<jupyter_start><jupyter_text>SHAP on R2B5 NARVAL data<jupyter_code>import os
import sys
import gc
import time
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
import shap

from tensorflow.keras.models import load_model
from tensorflow import nn 

# Add path with my_classes to sys.path
sys.path.insert(0, '/pf/b/b309170/workspace_icon-ml/cloud_cover_parameterization/')

import importlib
import my_classes
importlib.reload(my_classes)

from my_classes import load_data
from my_classes import read_mean_and_std

seed = int(sys.argv[1]) #Usually 10 or 100
np.random.seed(seed)
root_path = '/pf/b/b309170'
# data_path = os.path.join(root_path,
#                          'my_work/NARVAL/data_var_vertinterp_R02B05')
model_path = os.path.join(root_path,
                          'workspace_icon-ml/cloud_cover_parameterization/grid_column_based/saved_models')<jupyter_output><empty_output><jupyter_text>We use SHAP on the R2B5 NARVAL data (Key Point 3)<jupyter_code># Load data
ORDER_OF_VARS_NARVAL = ['qv', 'qc', 'qi', 'temp', 'pres', 'rho', 'zg', 'fr_lake', 'clc']
# data_dict = load_data(source='narval', days='all', vert_interp=True, resolution='R02B05', 
#                     order_of_vars=ORDER_OF_VARS_NARVAL)
data_dict = load_data(source='narval', days='all', vert_interp=True, resolution='R02B05', 
                    order_of_vars=ORDER_OF_VARS_NARVAL)
# Are there any bad data points
ta_is_0 = np.where(data_dict['temp'] == 0)
for i in range(3):
    assert ta_is_0[i].size == 0

del ta_is_0
gc.collect()
(TIME_STEPS, VERT_LAYERS, HORIZ_FIELDS) = data_dict['clc'].shape
try:
    #Reshaping into nd-arrays of equaling shapes (don't reshape in the vertical)
    data_dict['zg'] = np.repeat(np.expand_dims(data_dict['zg'], 0), TIME_STEPS, axis=0)
    data_dict['fr_lake'] = np.repeat(np.expand_dims(data_dict['fr_lake'], 0), TIME_STEPS, axis=0)
except:
    pass
# One sample should contain a column of information
data_dict_reshaped = {}

for key in data_dict.keys():
    if data_dict[key].shape[1] == VERT_LAYERS:  
        # Removing data above 21kms
        for i in range(4, VERT_LAYERS):
            new_key = '{}{}{:d}'.format(key,'_',(i+17)) # Should start at 21
            data_dict_reshaped[new_key] = np.reshape(data_dict[key][:,i,:], -1)
    else:
        data_dict_reshaped[key] = np.reshape(data_dict[key], -1)
#Converting dict into a DataFrame-object 
df = pd.DataFrame.from_dict(data_dict_reshaped)
df.head()
#Modifies df as well
def split_input_output(dataset):
    output_df = pd.DataFrame()
    for i in range(21, 48):
        output_df['clc_%d'%i] = dataset['clc_%d'%i] # Should start at 21
        del dataset['clc_%d'%i]
    return output_df
output_df = split_input_output(df)

input_narval = np.float32(df)
output_narval = np.float32(output_df)<jupyter_output><empty_output><jupyter_text>Remove columns that were constant in at least one of the training folds<jupyter_code>remove_fields = [27, 162, 163, 164]
(samples_narval, no_of_features) = input_narval.shape

input_narval = np.delete(input_narval, remove_fields, axis=1)
no_of_features = no_of_features - len(remove_fields)<jupyter_output><empty_output><jupyter_text>Scale the data according to the model's mean and std<jupyter_code>mean_1, std_1 = read_mean_and_std(model_path+'/model_grid_column_based_final_1.txt')
input_narval = (input_narval-mean_1)/std_1<jupyter_output><empty_output><jupyter_text>Load the model<jupyter_code>os.listdir(model_path)
# custom_objects = {}
# custom_objects['leaky_relu'] = nn.leaky_relu

fold_1 = 'model_grid_column_based_final_1.h5'

model_fold_1 = load_model(os.path.join(model_path, fold_1))<jupyter_output><empty_output><jupyter_text>Setting up SHAP<jupyter_code># Load column-based R2B5 training data
path_train_data = '/pf/b/b309170/my_work/icon-ml_data/cloud_cover_parameterization/grid_column_based/based_on_var_interpolated_data'
input_train = np.load(path_train_data + '/cloud_cover_input_train_1.npy')
no_samples_train = 10000

rand_indices_train = np.random.randint(0, input_train.shape[0], no_samples_train)
# Load column-based R2B5 training data
input_train = np.float32(input_train)[rand_indices_train]
#initialize js methods for visualization
shap.initjs()

# It takes any combination of a model and masker (usually the training data!) and 
# returns a callable subclass object that implements the particular estimation algorithm 
# (e.g. 'kernel' or 'deep') that was chosen.
# Data is taken to compute the base value.
explainer_shap = shap.DeepExplainer(model=model_fold_1, data=input_train)
explainer_shap.expected_value
#I doubt we can look at much more than 10000 at a time
# for no_samples_narval in [10, 10**2, 10**3, 10**4, 10**5]:
for no_samples_narval in [10000]:
    rand_indices_narval = np.random.randint(0, input_narval.shape[0], no_samples_narval)
    t0 = time.time()
    # Estimate the SHAP values on a subset of the data (you can try all but then gets slower)
    shap_values = explainer_shap.shap_values(X=input_narval[rand_indices_narval])     
    elapsed_time = time.time() - t0
    
    # Print to a file: no_samples_narval, elapsed_time, mean shap values, 
    # variance shap values, mean absolute shap values
    with open('/pf/b/b309170/workspace_icon-ml/iconml_clc/additional_content/shap_values/averaged_shap_values/r2b4_column-based_on_narval_r2b5.txt', 'a') as file:
        file.write('Number of NARVAL samples: %d\n'%no_samples_narval)
        file.write('Elapsed time: %.3f\n'%elapsed_time)
        file.write('Mean SHAP values:\n')
        for i in range(27):
            file.write(str(list(np.around(np.mean(shap_values[i], axis=0), 3))))
        file.write('\nVariance SHAP values:\n')
        for i in range(27):
            file.write(str(list(np.around(np.var(shap_values[i], axis=0), 3))))
        file.write('\nMean absolute SHAP values:\n')
        for i in range(27):
            file.write(str(list(np.around(np.mean(np.abs(shap_values[i]), axis=0), 3))))
        file.write('\n\n')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11713055003313452, 'fraction_numerical': 0.026507620941020542, 'mean_word_length': 4.799231508165226, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '531210', 'license': 'MIT', 'n_tokens_mistral': 2313, 'n_tokens_neox': 2199, 'n_words': 469, 'path': 'n2_column_based_narval_r2b4/source_code/average_absolute_shap_values_on_narval_r2b5.ipynb', 'repo_name': 'agrundner24/iconml_clc', 'size': '6036'}