notebooks/xnn_notebook_hmda.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train an XNN model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to implement using Keras (with TensorFlow backend) an Explainable Neural Network as described in [Explainable Neural Networks based on Additive Index Models](https://arxiv.org/pdf/1806.01933.pdf).\n",
    "\n",
    "The architecture of the network is as follows:\n",
    "\n",
    "![XNN Architecture](xnn_arch.png)\n",
    "[Explainable Neural Networks based on Additive Index Models](https://arxiv.org/pdf/1806.01933.pdf)\n",
    "\n",
    "And consists of three layers:\n",
    "\n",
    "(i) The projection layer (first hidden layer) using linear activation function\n",
    "\n",
    "(ii) Subnetworks, which learn a potentially nonlinear transformation of the input\n",
    "\n",
    "(iii) Combination layer calculates a weighted sum the output of the ridge functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import shap\n",
    "import subprocess\n",
    "import sys\n",
    "import pydot\n",
    "\n",
    "import keras\n",
    "from keras import backend\n",
    "from keras.layers import Activation, Add, Dense, Dropout, Input, Lambda, Concatenate\n",
    "from keras.models import Model\n",
    "from keras.utils import plot_model\n",
    "\n",
    "import plotly.plotly as py\n",
    "import chart_studio.plotly as py\n",
    "import plotly.tools as tls\n",
    "import matplotlib.pyplot as plt\n",
    "    \n",
    "from timeit import default_timer as timer\n",
    "import tensorflow as tf\n",
    "from keras import backend as K\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 12345\n",
    "\n",
    "np.random.seed(seed)\n",
    "\n",
    "my_init = keras.initializers.RandomUniform(seed=seed)\n",
    "\n",
    "output_label = \"BLDS_run100\"\n",
    "out_dir = \"output100/\"\n",
    "\n",
    "output_to_files = False\n",
    "\n",
    "def projection_initializer(shape, dtype=None):\n",
    "   \n",
    "    inps = shape[0]\n",
    "    subs = shape[1]\n",
    "    if subs > pow(inps, 2) - 1:\n",
    "        raise ValueError(\"Currently we support only up to 2^features - 1 number of subnetworks.\")\n",
    "    \n",
    "    weights = []\n",
    "\n",
    "    for i in range(subs):\n",
    "        w = [0] * inps\n",
    "        w[i] = 1\n",
    "        weights.append(w)\n",
    "    return weights\n",
    "\n",
    "\n",
    "\n",
    "def alpha_beta(alpha, beta, X , R):\n",
    "    \"\"\" Calculate the layerwise backpropagation function \"\"\"\n",
    "    \n",
    "    positive_values = [item for item in X if item > 0]\n",
    "        \n",
    "    negative_values = [item for item in X if item < 0] \n",
    "        \n",
    "    ans = np.array([0.0]*len(X))\n",
    "        \n",
    "    \n",
    "    if len(positive_values) > 0:\n",
    "           \n",
    "        ans += alpha*np.array([item / float(sum(positive_values)) if item > 0 else 0 for item in X])\n",
    "\n",
    "    if len(negative_values) > 0:\n",
    " \n",
    "        ans += -beta * np.array([item / float(sum(negative_values)) if item < 0 else 0 for item in X]) \n",
    "\n",
    "    return ans*R\n",
    "\n",
    "\n",
    "def deep_lift(X_bar, X , R):\n",
    "    \n",
    "    \"\"\" Deep lift backpropagation function\"\"\"   \n",
    "    \n",
    "    ans =  np.array(X) - np.array(X_bar)\n",
    "    ans = ans / (sum(X) - sum(X_bar))     \n",
    "    \n",
    "    return ans*R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XNN Class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class XNN:\n",
    "    # define base model\n",
    "    def __init__(self, features, ridge_functions=3, arch=[20,12], bg_samples=100, seed=None, is_categorical=False):\n",
    "        self.seed = seed\n",
    "        self.bg_samples = bg_samples\n",
    "        self.is_categorical = is_categorical\n",
    "        \n",
    "        #\n",
    "        # Prepare model architecture\n",
    "        #\n",
    "        # Input to the network, our observation containing all the features\n",
    "        input = Input(shape=(features,), name='main_input')\n",
    "\n",
    "        # Input to ridge function number i is the dot product of our original input vector times coefficients\n",
    "        ridge_input = Dense(ridge_functions,\n",
    "                            name=\"projection_layer\",\n",
    "                                activation='linear')(input)\n",
    "        \n",
    "        self.ridge_networks = []\n",
    "        # Each subnetwork uses only 1 neuron from the projection layer as input so we need to split it\n",
    "        ridge_inputs = Lambda( lambda x: tf.split(x, ridge_functions, 1), name='lambda_1' )(ridge_input)\n",
    "        for i, ridge_input in enumerate(ridge_inputs):\n",
    "            # Generate subnetwork i\n",
    "            mlp = self._mlp(ridge_input, i, arch)\n",
    "            self.ridge_networks.append(mlp)\n",
    "                    \n",
    "        added = Concatenate(name='concatenate_1')(self.ridge_networks)\n",
    "        \n",
    "        # Add the correct output layer for the problem\n",
    "        if self.is_categorical:\n",
    "            out = Dense(1, activation='sigmoid', input_shape= (ridge_functions, ), name='main_output')(added)\n",
    "        else:\n",
    "            out = Dense(1, activation='linear', input_shape= (ridge_functions, ), name='main_output')(added)\n",
    "            \n",
    "        self.model = Model(inputs=input, outputs=out)\n",
    "        \n",
    "        optimizer = keras.optimizers.Adam(lr=0.001, beta_1=0.9, beta_2=0.999, decay=0.0, amsgrad=True)\n",
    "        \n",
    "        # Use the correct loss for the problem\n",
    "        if self.is_categorical:\n",
    "            self.model.compile(loss={'main_output': 'binary_crossentropy'}, optimizer=optimizer)\n",
    "        else:\n",
    "            self.model.compile(loss={'main_output': 'mean_squared_error'}, optimizer=optimizer)\n",
    "\n",
    "        self.explainer = None\n",
    "                \n",
    "        \n",
    "    def _mlp(self, input, idx, arch=[20,12], activation='relu'):\n",
    "        if len(arch) < 1:\n",
    "            return #raise exception\n",
    "        \n",
    "        # Hidden layers\n",
    "        mlp = Dense(arch[0], activation=activation, name='mlp_{}_dense_0'.format(idx), kernel_initializer=my_init)(input)\n",
    "        for i, layer in enumerate(arch[1:]):\n",
    "            mlp = Dense(layer, activation=activation, name='mlp_{}_dense_{}'.format(idx, i+1), kernel_initializer=my_init)(mlp)\n",
    "         \n",
    "\n",
    "        # Output of the MLP\n",
    "        mlp = Dense(1, \n",
    "                    activation='linear', \n",
    "                    name='mlp_{}_dense_last'.format(idx), \n",
    "                    kernel_regularizer=keras.regularizers.l1(1e-3),\n",
    "                    kernel_initializer=my_init)(mlp)\n",
    "        \n",
    "        return mlp\n",
    "    \n",
    "    def print_architecture(self):\n",
    "        self.model.summary()\n",
    "    \n",
    "    def fit(self, X, y, epochs=5, batch_size=128, validation_set=0.0, verbose=0):\n",
    "        inputs = {'main_input': X}\n",
    "\n",
    "\n",
    "        if validation_set == 0:\n",
    "            self.model.fit(inputs, y, epochs=epochs, batch_size=batch_size, validation_split=validation_set, verbose=verbose)\n",
    "        else:\n",
    "            self.model.fit(inputs, y, epochs=epochs, batch_size=batch_size, validation_data=validation_set, verbose=verbose)\n",
    "            \n",
    "        #\n",
    "        # Prepare the explainer\n",
    "        # \n",
    "        np.random.seed(self.seed)\n",
    "        if isinstance(X, pd.DataFrame):\n",
    "            background = X.iloc[np.random.choice(X.shape[0], self.bg_samples, replace=False)]\n",
    "        else:\n",
    "            background = X[np.random.choice(X.shape[0], self.bg_samples, replace=False)]\n",
    "\n",
    "        # Explain predictions of the model on the subset\n",
    "        self.explainer = shap.DeepExplainer(self.model, background)\n",
    "                    \n",
    "        \n",
    "    def predict(self, X, pred_contribs=False):\n",
    "        pred_start = timer()\n",
    "        preds = self.model.predict(X)\n",
    "        pred_end = timer()\n",
    "        print(\"Predictions took {}\".format(pred_end - pred_start))\n",
    "\n",
    "        if pred_contribs:\n",
    "            explainer_start = timer()\n",
    "\n",
    "            self.shap_values = self.explainer.shap_values(X)\n",
    "\n",
    "            explainer_end = timer()\n",
    "            print(\"Explainer took {}\".format(explainer_end - explainer_start))\n",
    "\n",
    "            concat_start = timer()\n",
    "\n",
    "            preds = np.concatenate((preds, self.shap_values[0], preds), axis=1)\n",
    "            preds[:,-1] = self.explainer.expected_value\n",
    "\n",
    "            concat_end = timer()\n",
    "            print(\"Concat took {}\".format(concat_end - concat_start))\n",
    "        return preds\n",
    "    \n",
    "    def plot_shap(self, X):\n",
    "        shap.summary_plot(self.shap_values, X)\n",
    "        \n",
    "    def save(self, filename):\n",
    "        self.model.save(filename)\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the dataset\n",
    "xnn_data_dir = '~/article-information-2019/data/xnn_output/'\n",
    "\n",
    "DATA = pd.read_csv(xnn_data_dir + 'train_transformed.csv')\n",
    "TEST = pd.read_csv(xnn_data_dir + 'test_transformed.csv')\n",
    "\n",
    "\n",
    "\n",
    "# Select features and split into target and feature sets\n",
    "selected_vars = ['term_360', 'conforming']\n",
    "selected_vars += ['debt_to_income_ratio_missing','loan_amount_std', 'loan_to_value_ratio_std']\n",
    "selected_vars += ['no_intro_rate_period_std', 'intro_rate_period_std']\n",
    "selected_vars += ['property_value_std', 'income_std', 'debt_to_income_ratio_std']\n",
    "\n",
    "\n",
    "target_var = \"high_priced\"\n",
    "\n",
    "X=DATA[selected_vars].values\n",
    "Y=DATA[target_var].values\n",
    "TEST_X = TEST[selected_vars].values\n",
    "TEST_Y = TEST[target_var].values\n",
    "features = X.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit XNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the XNN\n",
    "is_cat = True\n",
    "xnn = XNN(features=features, ridge_functions=features,arch=[12, 8], is_categorical= is_cat)\n",
    "\n",
    "#plot_model(xnn.model, to_file='model_regression.png')\n",
    "xnn.print_architecture()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_test_preds = {}\n",
    "cv_train_preds = {}\n",
    "\n",
    "epoch = []\n",
    "fold_list = list(set(DATA['cv_fold']))\n",
    "epochs_max = 15000\n",
    "\n",
    "# Make predictions on the CV folds\n",
    "for cv_fold in fold_list:\n",
    "    TRAIN = DATA[DATA['cv_fold'] != cv_fold]\n",
    "    VALID = DATA[DATA['cv_fold'] == cv_fold]\n",
    "    \n",
    "    X = TRAIN[selected_vars].values\n",
    "    Y = TRAIN[target_var].values\n",
    "    \n",
    "    X_VALID = VALID[selected_vars].values\n",
    "    Y_VALID = VALID[target_var].values\n",
    "    \n",
    "    xnn.fit(X, Y, epochs=epochs_max, batch_size=1024, validation_set=(X_VALID, Y_VALID), verbose=1)\n",
    "    if output_to_files:\n",
    "        xnn.save(out_dir + 'cv_' + str(cv_fold) + '_hmda_model.h5')\n",
    "    \n",
    "    # CV predictions\n",
    "    cv_test_preds[cv_fold] = xnn.predict(TEST_X, pred_contribs=True)\n",
    "    cv_train_preds[cv_fold] = xnn.predict(X_VALID, pred_contribs=True)\n",
    "    \n",
    "    if output_to_files:\n",
    "        pd.DataFrame(pd.concat([TEST, pd.DataFrame(cv_test_preds[cv_fold])], axis=1)).to_csv(out_dir + \"test_output_\" + str(cv_fold) + \"_\" + str(epochs_max) + \"_\" + output_label +\".csv\" , index=False)\n",
    "        pd.DataFrame(pd.concat([VALID, pd.DataFrame(cv_train_preds[cv_fold])], axis=1)).to_csv(out_dir + \"valid_output_\" + str(cv_fold) + \"_\" + str(epochs_max) + \"_\" + output_label + \".csv\" , index=False)    \n",
    "      \n",
    "# Run the model on the full training set and make predictions on the test set  \n",
    "X=DATA[selected_vars].values\n",
    "Y=DATA[target_var].values\n",
    "xnn.fit(X, Y, epochs=epochs_max, batch_size=1024, validation_set=0, verbose=1)\n",
    "\n",
    "test_preds = xnn.predict(TEST_X, pred_contribs=True)\n",
    "\n",
    "if output_to_files:\n",
    "    pd.DataFrame(pd.concat([TEST, pd.DataFrame(test_preds)], axis=1)).to_csv(out_dir + \"main_\" + str(epochs_max) + \"_\" + output_label + \".csv\" , index=False)\n",
    "    xnn.save(out_dir + 'final_hmda_model_100.h5')\n",
    "    \n",
    "     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Record layer information\n",
    "# Plot projection layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Record the inputs, outputs, weights, and biases\n",
    "import scipy as sp\n",
    "\n",
    "int_output = {}\n",
    "int_output2 = {}\n",
    "int_weights = {}\n",
    "int_bias = {}\n",
    "int_input = {}\n",
    "\n",
    "original_activations = {}\n",
    "\n",
    "\n",
    "x_labels = list(map(lambda x: 'x' + str(x+1), range(features)))\n",
    "\n",
    "intermediate_output = []\n",
    "\n",
    "# Record and plot the projection weights\n",
    "# \n",
    "weight_list = []\n",
    "for layer in xnn.model.layers:\n",
    "    \n",
    "    layer_name = layer.get_config()['name']\n",
    "    if layer_name != \"main_input\":\n",
    "        print(layer_name)\n",
    "        weights = layer.get_weights()\n",
    "        \n",
    "        \n",
    "        # Record the biases\n",
    "        try:\n",
    "            bias = layer.get_weights()[1]\n",
    "            int_bias[layer_name] = bias\n",
    "        except:\n",
    "            print(\"No Bias\")\n",
    "            \n",
    "                       \n",
    "        # Record outputs for the test set\n",
    "        intermediate_layer_model = Model(inputs=xnn.model.input, outputs=xnn.model.get_layer(layer_name).output)\n",
    "        if (is_cat) and (layer_name == 'main_output'):\n",
    "            int_output[layer_name] = sp.special.logit(intermediate_layer_model.predict(TEST_X))\n",
    "            int_output[layer_name + \"_p\"] = intermediate_layer_model.predict(TEST_X)\n",
    "        else:\n",
    "            int_output[layer_name] = intermediate_layer_model.predict(TEST_X)\n",
    "        \n",
    "        # Record the outputs from the training set\n",
    "        if is_cat and (layer_name == 'main_output'):\n",
    "            original_activations[layer_name] = sp.special.logit(intermediate_layer_model.predict(X))   \n",
    "            original_activations[layer_name + \"_p\"] = intermediate_layer_model.predict(X)\n",
    "        else:\n",
    "            original_activations[layer_name] = intermediate_layer_model.predict(X)        \n",
    "\n",
    "\n",
    "        # Record other weights, inputs, and outputs\n",
    "        int_weights[layer_name] = weights\n",
    "        int_input[layer_name] = layer.input\n",
    "        int_output2[layer_name] = layer.output\n",
    "               \n",
    "    # Plot the projection layers    \n",
    "    if \"projection_layer\" in layer.get_config()['name']:\n",
    "        \n",
    "        print(layer.get_config()['name'])\n",
    "        \n",
    "        # Record the weights for each projection layer\n",
    "        weights = [np.transpose(layer.get_weights()[0])]\n",
    "\n",
    "        weight_list2=[]\n",
    "        for i, weight in enumerate(weights[0]):\n",
    "            weight_list.append(weight)\n",
    "            weight_list2.append(list(np.reshape(weight, (1,features))[0]))\n",
    "        \n",
    "            print(weight)\n",
    "\n",
    "            # Plot weights\n",
    "            plt.bar(x_labels, abs(np.reshape(weight, (1,features))[0]), 1, color=\"blue\")\n",
    "            plt.xlabel(\"Subnetowork {} coefficient\".format(i))\n",
    "            plt.ylabel(\"Weight value\")\n",
    "            plt.show()\n",
    "\n",
    "    if \"main_output\" in layer.get_config()['name']:\n",
    "        weights_main = layer.get_weights()\n",
    "        print(weights_main)\n",
    "        \n",
    "if output_to_files:\n",
    "    pd.DataFrame(weight_list2).to_csv(\"wp_\" + output_label + \".csv\", index=False)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate the feature importances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Calculate ridge and input function local feature importances   \n",
    "item = 0\n",
    "\n",
    "feature_output = []\n",
    "feature_output2 = []\n",
    "feature_output3 = []\n",
    "\n",
    "# Find the average outputs\n",
    "S_bar = sum(original_activations[\"main_output\"])/len(original_activations[\"main_output\"])\n",
    "# original_activations[layer_name]\n",
    "output_weights = np.array([int_weights[\"main_output\"][0][ii][0] for ii in range(features)])\n",
    "output_Z_bar = sum(original_activations[\"concatenate_1\"]*output_weights)/len(original_activations[\"concatenate_1\"])\n",
    "\n",
    "\n",
    "# For each ridge function calculate the average input activation\n",
    "input_Z_bar = {}\n",
    "for ridge_num in range(features):   \n",
    "    input_weights = np.array([int_weights[\"projection_layer\"][0][ii][ridge_num] for ii in range(features)])\n",
    "    input_Z_bar[ridge_num] = sum(X*input_weights)/len(X)\n",
    "    \n",
    "    \n",
    "# For each test instance, calculate the feature importance scores    \n",
    "for test_num in range(len(TEST_X)):\n",
    "    \n",
    "    # Calculate the output activations\n",
    "    activation_list=[int_weights[\"main_output\"][0][ii][0]*int_output[\"concatenate_1\"][test_num][ii] for ii in range(features)]\n",
    "    \n",
    "    \n",
    "    # Calculate layerwise backpropagaiton to the ridge functions\n",
    "    # For classification, change this to the inverse sigmoid of the output\n",
    "    features_ab = alpha_beta(2, 1, activation_list , int_output[\"main_output\"][test_num][0])\n",
    "    features_ab2 = alpha_beta(2, 1, activation_list , int_output[\"main_output\"][test_num][0]-S_bar)\n",
    "    \n",
    "    # Calculate deep lift backpropagation to the ridge functions\n",
    "    features_dl = deep_lift(output_Z_bar, activation_list , int_output[\"main_output\"][test_num][0]-S_bar)\n",
    "      \n",
    "        \n",
    "    # Calculate the deep lift and layerwise information scores for the input layer\n",
    "    input_scores = []\n",
    "    input_scores_dl = []\n",
    "    input_scores2 = []\n",
    "    input_scores_dl2 = []\n",
    "    for ridge_num in range(features):\n",
    "        weights = int_weights[\"projection_layer\"][0][ridge_num]\n",
    "        output = TEST_X[test_num,:]\n",
    "        \n",
    "        # Calculate the activations from the projection layer\n",
    "        act = TEST_X[test_num,:]*np.array([int_weights[\"projection_layer\"][0][ii][ridge_num] for ii in range(features)])\n",
    "    \n",
    "        # Input relevance scores for a single ridge function\n",
    "        input_scores += list(alpha_beta(2, 1, act, features_ab[ridge_num]))\n",
    "        input_scores_dl += list(deep_lift(input_Z_bar[ridge_num], act, features_dl[ridge_num]))\n",
    "        input_scores2 += list(alpha_beta(2, 1, act, features_ab2[ridge_num]))\n",
    "\n",
    "    # Sum the contribution of the variable importance from each of the projections\n",
    "    input_sum = [sum(input_scores[ii+features*jj] for jj in range(features)) for ii in range(features)] \n",
    "    input_sum2 = [sum(input_scores2[ii+features*jj] for jj in range(features)) for ii in range(features)] \n",
    "    input_sum_dl = [sum(input_scores_dl[ii+features*jj] for jj in range(features)) for ii in range(features)] \n",
    "    input_abs_sum = [sum(abs(input_scores[ii+features*jj]) for jj in range(features)) for ii in range(features)] \n",
    "    \n",
    "    # Recored the feature importance information for this instance\n",
    "    feature_output.append(input_sum+input_abs_sum+[int_output[\"main_output\"][test_num][0]]+list(features_ab)+input_scores)\n",
    "    feature_output2.append(input_sum+list(features_ab)+input_sum_dl + list(features_dl))\n",
    "    feature_output3.append(input_sum2+list(features_ab2)+input_sum_dl + list(features_dl))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find and plot the ridge functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "intermediate_output = []\n",
    "\n",
    "for feature_num in range(features):\n",
    "    intermediate_layer_model = Model(inputs=xnn.model.input,\n",
    "                                 outputs=xnn.model.get_layer('mlp_'+str(feature_num)+'_dense_last').output)\n",
    "    intermediate_output.append(intermediate_layer_model.predict(X))\n",
    "\n",
    "\n",
    "# Record and plot the ridge functions\n",
    "ridge_x = []\n",
    "ridge_y = []\n",
    "for weight_number in range(len(weight_list)):\n",
    "    \n",
    "    ridge_x.append(list(sum(X[:, ii]*weight_list[weight_number][ii] for ii in range(features))))\n",
    "    ridge_y.append(list(intermediate_output[weight_number]))\n",
    "\n",
    "    plt.plot(sum(X[:, ii]*weight_list[weight_number][ii] for ii in range(features)), intermediate_output[weight_number], 'o')\n",
    "    plt.xlabel(\"x\")\n",
    "    plt.ylabel(\"Subnetwork \" + str(weight_number))\n",
    "    plt.legend(loc='lower right')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "if output_to_files:\n",
    "    pd.DataFrame(ridge_x).to_csv(\"ridge_x_\"+ output_label +\".csv\", index=False)\n",
    "    pd.DataFrame(ridge_y).to_csv(\"ridge_y_\" + output_label + \".csv\", index=False)     \n",
    "    pd.DataFrame(feature_output2).to_csv(\"feature_output2_\" + output_label + \".csv\", index=False)\n",
    "    pd.DataFrame(feature_output3).to_csv(\"feature_output3_\" + output_label + \".csv\", index=False)   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make predictions on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds = xnn.predict(TEST_X, pred_contribs=True)\n",
    "\n",
    "if output_to_files:\n",
    "    pd.DataFrame(preds).to_csv(\"preds_\" + output_label + \".csv\", index=False)\n",
    "    pd.DataFrame(TEST).to_csv(\"TEST_\" + output_label + \".csv\", index=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the Shapley values.\n",
    "shap.initjs()\n",
    "shap.summary_plot(xnn.shap_values, X)\n",
    "\n",
    "y=xnn.shap_values\n",
    "ind=1\n",
    "\n",
    "\n",
    "# Calculate the average absolute feature importance across the dataset\n",
    "\n",
    "layerwise_average_input=np.array([0.0]*features)\n",
    "layerwise_average_input2=np.array([0.0]*features)\n",
    "layerwise_average_ridge=np.array([0.0]*features)\n",
    "layerwise_average_ridge2=np.array([0.0]*features)\n",
    "layerwise_average_shap=np.array([0.0]*features)\n",
    "lift_average_input=np.array([0.0]*features)\n",
    "lift_average_ridge=np.array([0.0]*features)\n",
    "\n",
    "for ii in range(len(feature_output2)):\n",
    "    layerwise_average_input += np.abs(np.array(feature_output2[ii][0:features]))\n",
    "    layerwise_average_ridge += np.abs(np.array(feature_output2[ii][features:(2*features)]))\n",
    "    layerwise_average_input2 += np.abs(np.array(feature_output3[ii][0:features]))\n",
    "    layerwise_average_ridge2 += np.abs(np.array(feature_output3[ii][features:(2*features)]))\n",
    "    lift_average_input += np.abs(np.array(feature_output2[ii][(2*features):(3*features)]))\n",
    "    lift_average_ridge += np.abs(np.array(feature_output2[ii][(3*features):(4*features)]))\n",
    "    layerwise_average_shap += np.abs(np.array(y[0][ii]))\n",
    "     \n",
    "layerwise_average_input = layerwise_average_input/len(feature_output2)\n",
    "layerwise_average_ridge = layerwise_average_ridge/len(feature_output2)\n",
    "layerwise_average_input2 = layerwise_average_input2/len(feature_output2)\n",
    "layerwise_average_ridge2 = layerwise_average_ridge2/len(feature_output2)\n",
    "layerwise_average_shap = layerwise_average_shap/len(feature_output2)\n",
    "lift_average_input = lift_average_input/len(feature_output2)\n",
    "lift_average_ridge = lift_average_ridge/len(feature_output2)\n",
    "\n",
    "\n",
    "SCORES = [list(layerwise_average_input), list(layerwise_average_ridge),\n",
    "          list(layerwise_average_input2), list(layerwise_average_ridge2),\n",
    "          list(layerwise_average_shap), list(lift_average_input),\n",
    "          list(lift_average_ridge)]\n",
    "\n",
    "if output_to_files:\n",
    "    pd.DataFrame(SCORES).to_csv(\"scores_\" + output_label + \".csv\", index=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot feature importance scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.bar(x, abs(np.reshape(y[0][ind], (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Shap Score Example \" + str(ind))\n",
    "plt.ylabel(\"\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(feature_output2[ind][0:features], (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Input Layerwise Propagation Score Example \" + str(ind))\n",
    "plt.ylabel(\"\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(feature_output2[ind][features:(2*features)], (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Ridge Layerwise Propagation Score Example \" + str(ind))\n",
    "plt.ylabel(\"Weight value\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(feature_output2[ind][2*features:(3*features)], (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Deep Lift Input Score Example \" + str(ind))\n",
    "plt.ylabel(\"Weight value\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(feature_output2[ind][3*features:(4*features)], (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Deep Lift Ridge Score Example \" + str(ind))\n",
    "plt.ylabel(\"Weight value\")\n",
    "plt.show()\n",
    "      \n",
    "plt.bar(x, abs(np.reshape(layerwise_average_input, (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Input Layerwise Propagation Score Average\")\n",
    "plt.ylabel(\"\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(layerwise_average_ridge, (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Ridge Layerwise Propagation Score Average\")\n",
    "plt.ylabel(\"Weight value\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(layerwise_average_input2, (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Input Layerwise Propagation Score Average 2\")\n",
    "plt.ylabel(\"\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(layerwise_average_ridge2, (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Ridge Layerwise Propagation Score Average 2\")\n",
    "plt.ylabel(\"Weight value\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(lift_average_input, (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Input Lift Score Average\")\n",
    "plt.ylabel(\"\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(lift_average_ridge, (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Ridge Lift Score Average\")\n",
    "plt.ylabel(\"Weight value\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(x, abs(np.reshape(layerwise_average_shap, (1,features))[0]), 1, color=\"blue\")\n",
    "plt.xlabel(\"Shapley Score Average\")\n",
    "plt.ylabel(\"Weight value\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}