{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"# 0, Introduction\nThe MNIST data are predicted using CNN and then the predicted results are explained using SHAP.\n\n**I won't go into how it works. Usage only.**\n\nIf you only want to see how to use SHAP, please refer to 5, as it is basically the same as usual until the model is created.","metadata":{}},{"cell_type":"code","source":"import os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:26.686371Z","iopub.execute_input":"2022-11-14T16:30:26.686698Z","iopub.status.idle":"2022-11-14T16:30:26.708688Z","shell.execute_reply.started":"2022-11-14T16:30:26.686662Z","shell.execute_reply":"2022-11-14T16:30:26.707789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1, import library\nThe only library required by shap is \"shap\". It's simple.","metadata":{}},{"cell_type":"code","source":"import shap\nimport tensorflow as tf\ntf.compat.v1.disable_v2_behavior()\n\n# from keras.layers import concatenate\nfrom tensorflow.keras.layers import Reshape, Lambda, RepeatVector, concatenate  # 注意全用tf.keras或keras\nfrom sklearn.metrics import r2_score as R2\nimport numpy as np\nfrom sklearn.metrics import mean_squared_error as mse\nfrom sklearn.metrics import mean_absolute_error as mae\nfrom tensorflow.keras import backend  # 用于在训练中添加metrics，如RMSE函数和R2\n\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Dense\nfrom tensorflow.keras.layers import LSTM  # LSTM\nfrom tensorflow.keras.layers import Conv2D, MaxPooling2D, Dense, Flatten, Activation\nfrom tensorflow.keras.layers import Conv3D, MaxPooling3D, Dropout\nfrom tensorflow.keras.layers import BatchNormalization\nfrom tensorflow.keras import optimizers\nfrom matplotlib import pyplot as plt\n\nfrom tensorflow import keras as keras\n# import keras\nfrom sklearn.metrics import r2_score as R2\nimport os\n","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:26.710506Z","iopub.execute_input":"2022-11-14T16:30:26.711108Z","iopub.status.idle":"2022-11-14T16:30:34.326284Z","shell.execute_reply.started":"2022-11-14T16:30:26.711067Z","shell.execute_reply":"2022-11-14T16:30:34.325463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cnn_module = Sequential([\n    BatchNormalization(),\n    # 注意第一层中变量的特征\n    Conv3D(filters=64, kernel_size=(3, 3, 1), padding='same', input_shape=(7, 7, 4, 18),\n           kernel_regularizer=None),  # 输出空间的维度   # regularizers.l2(0.01)\n    BatchNormalization(),\n    Activation('relu'),\n\n    Conv3D(filters=128, kernel_size=(3, 3, 1), padding='same', kernel_regularizer=None),  # 输出空间的维度\n    BatchNormalization(),\n    Activation('relu'),\n\n    Conv3D(filters=256, kernel_size=(3, 3, 1), padding='same', kernel_regularizer=None),  # 输出空间的维度\n    BatchNormalization(),\n    Activation('relu'),\n\n    Conv3D(filters=128, kernel_size=(3, 3, 1), kernel_regularizer=None),  # 输出空间的维度\n    BatchNormalization(),\n    Activation('relu'),\n\n    Conv3D(filters=64, kernel_size=(3, 3, 1), kernel_regularizer=None),  # 输出空间的维度\n    BatchNormalization(),\n    Activation('relu'),  # 这里的输出是 batch * 64 * 3 * 3 或者 batch * 64 * 3 * 3 * slen(4)\n\n    MaxPooling3D(pool_size=(3, 3, 1)),  # 这里的输出是 N * 1 * 1 *4 * 64\n])\n\nlstm_module = Sequential([  # 这个输入应该是 N * 4 * 64\n    LSTM(64, return_sequences=True,  # 这里的input—size默认为1  return_sequences=True,\n         input_shape=(4, 64), kernel_regularizer=None),  # 这里对应的是slen 以及cnn_output对应的特征维度\n    BatchNormalization(),\n    Activation('relu'),\n\n    LSTM(64, return_sequences=False, kernel_regularizer=None),  # 这此时的输出是N * num_directional * num_unit\n    BatchNormalization(),\n    Activation('relu'),\n\n])\n\ndnn_module = Sequential([\n    Dense(64, kernel_regularizer=None),\n    BatchNormalization(),\n    Activation('relu'),\n\n    Dense(1),\n])\n\ninput_layer = tf.keras.Input(shape=(7, 7, 4, 18))\n# input_at = attention_block_channel2(input_layer)  # 为了构建模型，input_at不能使用input_layer代替\n# cnn_out = cnn_module(input_at)\n\ncnn_out = cnn_module(input_layer)\nprint(cnn_out)\nlstm_input = tf.squeeze(cnn_out, axis=1)\nlstm_input = tf.squeeze(lstm_input, axis=1)\nprint(lstm_input, lstm_input.shape)\nlstm_output = lstm_module(lstm_input)\nprint(lstm_output)\ndnn_output = dnn_module(lstm_output)\nprint(dnn_output)\nmodel = tf.keras.Model(inputs=input_layer, outputs=dnn_output)","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:34.327662Z","iopub.execute_input":"2022-11-14T16:30:34.327944Z","iopub.status.idle":"2022-11-14T16:30:35.328279Z","shell.execute_reply.started":"2022-11-14T16:30:34.327907Z","shell.execute_reply":"2022-11-14T16:30:35.327421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.load_weights(\"../input/model-cnnlstm/best_weights.h5\")\n\n\n# if load model failed, retrain\n\ndef me_rmse(y_true, y_pred):  # 添加自定义的metric\n    return backend.sqrt(backend.mean(backend.square(y_pred - y_true), axis=0))  # 这里注意维度，由于预测值是N*1，这里开方的维度是0\n\n\ndef me_r2(y_test, pred_test):\n    return 1 - backend.sum(backend.square(pred_test - y_test)) / backend.sum(\n        backend.square(backend.mean(y_test) - y_test))  \n\nfrom sklearn.model_selection import train_test_split\nX = np.load('../input/fulldata-cnnlstm/feature_st.npy').transpose(1, 3, 4, 0, 2)\ny = np.load('../input/fulldata-cnnlstm/label_st.npy')\nX1000, X_test, y1000, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n# n_epochs = 30\n# opt = optimizers.Adam(lr=0.0001)\n\n# model.compile(optimizer=opt, loss='mean_squared_error', metrics=[me_rmse, me_r2]) \n# history = model.fit(x=X1000, y=y1000, epochs=n_epochs,\n#                     batch_size=64, validation_split=0.1) ","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:35.32985Z","iopub.execute_input":"2022-11-14T16:30:35.330145Z","iopub.status.idle":"2022-11-14T16:30:46.534691Z","shell.execute_reply.started":"2022-11-14T16:30:35.330108Z","shell.execute_reply":"2022-11-14T16:30:46.533908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# X_shap = shap.utils.sample(X1000, 1000,random_state=42)\nX_shap = X1000[:1000]","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:33:05.599055Z","iopub.execute_input":"2022-11-14T16:33:05.599359Z","iopub.status.idle":"2022-11-14T16:33:05.603532Z","shell.execute_reply.started":"2022-11-14T16:33:05.599326Z","shell.execute_reply":"2022-11-14T16:33:05.602539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"explainer = shap.DeepExplainer(model, X_shap)\n# 结果是1000*7*7*4*18\nshap_values = explainer.shap_values(X_shap)  # 计算每一个样本的shap值`","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:33:05.795292Z","iopub.execute_input":"2022-11-14T16:33:05.795564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.predict(X1000).mean()","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.940348Z","iopub.status.idle":"2022-11-14T16:30:46.940939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 原始的代码（参考）","metadata":{}},{"cell_type":"code","source":"# #shap!!\n# import shap\n\n# import pandas as pd\n# import numpy as np\n# import tensorflow as tf\n# import tensorflow.keras.backend as K\n# import matplotlib.pyplot as plt\n# plt.style.use('ggplot')\n\n# from PIL import Image\n# from sklearn.model_selection import train_test_split\n# from tensorflow.keras.preprocessing.image import ImageDataGenerator\n# from tensorflow.keras.utils import to_categorical\n# from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping\n# from tensorflow.keras.models import Model, Sequential\n# from tensorflow.keras.layers import Dense, Dropout, Flatten, Conv2D, MaxPooling2D","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.942256Z","iopub.status.idle":"2022-11-14T16:30:46.942652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2, Load data\nThe process is the same as usual.","metadata":{}},{"cell_type":"code","source":"# train_df = pd.read_csv(\"../input/digit-recognizer/train.csv\")\n# test_df = pd.read_csv(\"../input/digit-recognizer/test.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.945189Z","iopub.status.idle":"2022-11-14T16:30:46.945795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X = (train_df.iloc[:,1:].values).astype('float32')\n# Y = train_df.iloc[:,0].values.astype('int32')\n# test = test_df.values.astype('float32')","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.946955Z","iopub.status.idle":"2022-11-14T16:30:46.947553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X = X.reshape(X.shape[0], 28, 28, 1)\n# Y = to_categorical(Y)\n# test = test.reshape(test.shape[0], 28, 28, 1)","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.948622Z","iopub.status.idle":"2022-11-14T16:30:46.949213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# #Confirm the image.\n# plt.imshow(X[0].reshape(28,28))\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.950301Z","iopub.status.idle":"2022-11-14T16:30:46.950922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X = X.astype(\"float32\") / 255\n# test = test.astype(\"float32\") / 255","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.95199Z","iopub.status.idle":"2022-11-14T16:30:46.952598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X_train, X_val, Y_train, Y_val = train_test_split(X, Y, test_size = 0.1, random_state=0)","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.953694Z","iopub.status.idle":"2022-11-14T16:30:46.954277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X_train.shape, X_val.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.955346Z","iopub.status.idle":"2022-11-14T16:30:46.955944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# datagen = ImageDataGenerator(\n#         featurewise_center=False,  # set input mean to 0 over the dataset\n#         samplewise_center=False,  # set each sample mean to 0\n#         featurewise_std_normalization=False,  # divide inputs by std of the dataset\n#         samplewise_std_normalization=False,  # divide each input by its std\n#         zca_whitening=False,  # apply ZCA whitening\n#         rotation_range=10,  # randomly rotate images in the range (degrees, 0 to 180)\n#         zoom_range = 0.1, # Randomly zoom image \n#         width_shift_range=0.1,  # randomly shift images horizontally (fraction of total width)\n#         height_shift_range=0.1,  # randomly shift images vertically (fraction of total height)\n#         horizontal_flip=False,  # randomly flip images\n#         vertical_flip=False)  # randomly flip images\n\n# datagen.fit(X_train)","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.957056Z","iopub.status.idle":"2022-11-14T16:30:46.95766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X_train.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.958765Z","iopub.status.idle":"2022-11-14T16:30:46.959344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3, Define the CNN model\nBy the way, When I use the BATCH NORMARIZER, which is very popular nowadays, I get an error when I get a shap value...\n\n(⇒LookupError: gradient registry has no entry for: shap_FusedBatchNormV3)\n\nThat's why I don't use BATCH NORMARIZER.","metadata":{}},{"cell_type":"code","source":"# # 注意，BN层可能SHAP无法解释\n# model = Sequential()\n# model.add(Conv2D(32, kernel_size = 4, activation=\"relu\", input_shape=(28,28,1)))\n\n# model.add(MaxPooling2D(pool_size=(2, 2)))\n# model.add(Conv2D(64, kernel_size = 3, activation=\"relu\"))\n# model.add(MaxPooling2D(pool_size=(2, 2)))\n# model.add(Conv2D(64, kernel_size = 2, activation=\"relu\"))\n# model.add(MaxPooling2D(pool_size=(2, 2)))\n\n# model.add(Flatten())\n# model.add(Dense(128, activation=\"relu\"))\n# model.add(Dropout(0.25))\n# model.add(Dense(10, activation='softmax'))\n\n# model.compile(loss=\"mean_squared_error\", optimizer=\"rmsprop\", metrics=['accuracy'])\n\n# model.summary()","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.960457Z","iopub.status.idle":"2022-11-14T16:30:46.961073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# learning_rate_reduction = ReduceLROnPlateau(monitor='val_acc', \n#                                             patience=3, \n#                                             verbose=1, \n#                                             factor=0.5, \n#                                             min_lr=0.00001)\n# # 学习率的调度，factor是学习率的缩放因子，min_lr是最小的学习率。需要将Patience设置为小于早停的Patience才有效","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.96217Z","iopub.status.idle":"2022-11-14T16:30:46.96277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# batch_size = 100\n# epochs = 30\n\n# hist = model.fit_generator(datagen.flow(X_train,Y_train, batch_size=batch_size),\n#                               epochs = epochs, validation_data = (X_val,Y_val),\n#                               verbose = 2, steps_per_epoch=X_train.shape[0] // batch_size\n#                               , callbacks=[learning_rate_reduction])","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.963842Z","iopub.status.idle":"2022-11-14T16:30:46.96441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plt.figure(figsize=(8, 4))\n# plt.plot(hist.history[\"loss\"])\n# plt.title(\"Loss\")\n# plt.ylabel(\"Loss\")\n# plt.xlabel(\"Epoch\")\n# plt.legend([\"Train\", \"Test\"])\n# plt.show()\n\n# pred = model.predict(X_val)\n\n# test_loss, test_acc = model.evaluate(X_val, Y_val, verbose=0)\n# print(f'\\nTest accuracy: {test_acc}')","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.965526Z","iopub.status.idle":"2022-11-14T16:30:46.966128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4, Check error\nCheck the prediction results that were wrong.","metadata":{}},{"cell_type":"code","source":"# Y_pred_classes = np.argmax(pred,axis = 1) \n# Y_true = np.argmax(Y_val,axis = 1)","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.967191Z","iopub.status.idle":"2022-11-14T16:30:46.967791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Y_pred_classes.shape, Y_true.shape, pred.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.968862Z","iopub.status.idle":"2022-11-14T16:30:46.969445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# errors = (Y_pred_classes - Y_true != 0)\n\n# Y_pred_classes_errors = Y_pred_classes[errors]\n# Y_pred_errors = pred[errors]\n# Y_true_errors = Y_true[errors]\n# X_val_errors = X_val[errors]\n\n# def display_errors(errors_index,img_errors,pred_errors, obs_errors):\n#     \"\"\" This function shows 10 images with their predicted and real labels\"\"\"\n#     n = 0\n#     nrows = 2\n#     ncols = 5\n#     fig, ax = plt.subplots(nrows,ncols,sharex=True,sharey=True,figsize=(20, 8))\n#     for row in range(nrows):\n#         for col in range(ncols):\n#             error = errors_index[n]\n#             ax[row,col].imshow((img_errors[error]).reshape((28,28)))\n#             ax[row,col].set_title(\"Predicted label :{}\\nTrue label :{}\".format(pred_errors[error],obs_errors[error]))\n#             n += 1\n#     plt.show()\n    \n# Y_pred_errors_prob = np.max(Y_pred_errors,axis = 1)\n# true_prob_errors = np.diagonal(np.take(Y_pred_errors, Y_true_errors, axis=1))\n# delta_pred_true_errors = Y_pred_errors_prob - true_prob_errors\n# sorted_dela_errors = np.argsort(delta_pred_true_errors)\n# most_important_errors = sorted_dela_errors[-10:]\n\n# display_errors(most_important_errors, X_val_errors, Y_pred_classes_errors, Y_true_errors)","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.970576Z","iopub.status.idle":"2022-11-14T16:30:46.971168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5, SHAP value\nNow it's time for the SHAP, to calculate the SHAP value and see where the model sees it and what it expects.","metadata":{}},{"cell_type":"code","source":"# X_val_errors.shape, most_important_errors.shape, X_train.shape, X_val_errors[[i]].shape","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.97224Z","iopub.status.idle":"2022-11-14T16:30:46.972846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X_val_errors.shape, X_val_errors[[0]].shape, ","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.973906Z","iopub.status.idle":"2022-11-14T16:30:46.974486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# #When describing Deep Learning, we use DeepExplainer.\n# explainer = shap.DeepExplainer(model, (X[0:1000]))\n\n# #Check out the 10 data that the model has mistakenly predicted.\n# for i in most_important_errors:  \n    \n#     #Calculates the SHAP value.\n#     shap_values = explainer.shap_values(X_val_errors[[i]])\n    \n#     #The following two lines are extras. It works even if you turn it off.\n#     index_names = np.array([str(x) + \"\\n\" + '{:>7.3%}'.format(Y_pred_errors[i][x]) for x in range(10)]).reshape(1,10)\n#     print(\"Predicted label :{}\\nTrue label :{}\".format(Y_pred_classes_errors[i],Y_true_errors[i]))\n    \n#     #Displays the results.\n#     shap.image_plot(shap_values, X_val_errors[[i]] ,index_names ,show=False)\n#     plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.975593Z","iopub.status.idle":"2022-11-14T16:30:46.976173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The blue areas are pixels that were not what model would have expected as its label. The red areas are the opposite of that.\n\nFocusing on these areas may help to further improve the accuracy of the model.","metadata":{}},{"cell_type":"markdown","source":"# 6, submission\nI'm going to submit it. Thank you for taking the time to watch it to the end.","metadata":{}},{"cell_type":"code","source":"# results = model.predict(test)\n# results = np.argmax(results,axis = 1)\n# results = pd.Series(results,name=\"Label\")","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.977294Z","iopub.status.idle":"2022-11-14T16:30:46.978001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submission = pd.concat([pd.Series(range(1,28001),name = \"ImageId\"),results],axis = 1)\n# submission.to_csv(\"cnn_mnist_datagen.csv\",index=False)","metadata":{"execution":{"iopub.status.busy":"2022-11-14T16:30:46.979102Z","iopub.status.idle":"2022-11-14T16:30:46.979691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}