araucana_AI_special_issue.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/detsutut/AraucanaXAI/blob/master/araucana_XAI_special_issue.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B2E--xGOJYum"
      },
      "source": [
        "#**Araucana XAI Special Issue Notebook**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Rkxkcaf7JVXo"
      },
      "source": [
        "## **Init**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "na_QGmtfwdda"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "!pip3 install shap\n",
        "!pip3 install araucanaxai\n",
        "!pip3 install lime"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XV7J-kvWcJfa"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "from google.colab import files\n",
        "import numpy as np\n",
        "from itertools import combinations\n",
        "import pandas as pd\n",
        "import io\n",
        "import time\n",
        "import warnings\n",
        "from tqdm import tqdm\n",
        "from matplotlib import pyplot as plt\n",
        "from IPython.display import display\n",
        "import collections\n",
        "import random\n",
        "\n",
        "import araucanaxai\n",
        "import shap\n",
        "import lime\n",
        "\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier\n",
        "from sklearn.linear_model import LogisticRegression\n",
        "from sklearn.neural_network import MLPClassifier\n",
        "from sklearn.preprocessing import MinMaxScaler\n",
        "from sklearn.metrics import precision_recall_curve, f1_score, accuracy_score, precision_score, recall_score, confusion_matrix\n",
        "\n",
        "# still have to fix this (the warning is just annoying but does not affect the pipeline in any way)\n",
        "warnings.filterwarnings(\"ignore\", message=\"X does not have valid feature names, but LogisticRegression was fitted with feature names\")\n",
        "warnings.filterwarnings(\"ignore\", message=\"X does not have valid feature names, but RandomForestClassifier was fitted with feature names\")\n",
        "warnings.filterwarnings(\"ignore\", message=\"X does not have valid feature names, but DecisionTreeClassifier was fitted with feature names\")\n",
        "warnings.filterwarnings(\"ignore\", message=\"X does not have valid feature names, but GradientBoostingClassifier was fitted with feature names\")\n",
        "warnings.filterwarnings(\"ignore\", message=\"X does not have valid feature names, but MLPClassifier was fitted with feature names\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t4eEbOQGJP7Y"
      },
      "source": [
        "## **Data**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ggAFbLC_cMsd"
      },
      "outputs": [],
      "source": [
        "#UPLOAD DATA\n",
        "uploaded = files.upload()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kPPBA2d4yP30"
      },
      "source": [
        "### MIMIC Dataset"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ti8bX18LLpF5"
      },
      "outputs": [],
      "source": [
        "#MIMIC PREPROCESSING\n",
        "df_raw = pd.read_csv(io.BytesIO(uploaded['MIMIC_III_dataset.csv']))\n",
        "#df_raw = df_raw.groupby('In-hospital_death', group_keys=False).apply(lambda x: x.sample(frac=0.1))  # ********* <<< DEBUG, remove this row >>> **********\n",
        "\n",
        "# Removal of features containing more than 90% of NaN\n",
        "original_shape = df_raw.shape\n",
        "print(\"DataFrame shape before NaN removal:\", original_shape)\n",
        "df = df_raw.dropna(thresh=round(df_raw.shape[0] * 0.90), axis=1)\n",
        "\n",
        "# NaN removal by removing rows with at least one NaN value\n",
        "# Check for NaN values\n",
        "if df.isnull().any().any():  # If there's at least one column with NaN, the output is \"True\"\n",
        "    old_row_size = original_shape[0]\n",
        "    df = df.dropna(axis=0, how='any')\n",
        "    print(\"DataFrame shape after NaN removal:\", df.shape)\n",
        "    new_row_size = df.shape[0]\n",
        "    print(\"Relative percentage of removed rows (wrt the old row size): %.2f\" % (\n",
        "            (old_row_size - new_row_size) / old_row_size * 100))\n",
        "\n",
        "# NaN values have now been removed, so the output should be \"False\"\n",
        "print(\"Is there any NaN after NaN removal?\", df.isnull().any().any())\n",
        "df = df.drop(columns=[\"recordid\"])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vqvlOovnMjZi"
      },
      "outputs": [],
      "source": [
        "#DATA PREPARATION\n",
        "pred_col = 'In-hospital_death' #edit this line when using another dataset\n",
        "\n",
        "X = df\n",
        "y = X[pred_col]\n",
        "X_feat = X.drop(columns=[pred_col])\n",
        "feature_names = X_feat.columns\n",
        "\n",
        "#check which features has less than 10 unique values to identify a subset of probable categorical features\n",
        "print(feature_names[np.where([len(np.unique(X_feat.iloc[:,i]))<10 for i in range(X_feat.shape[1])])])\n",
        "categorical_feat = ['CCU', 'CSU', 'SICU'] # edit this line when using another dataset\n",
        "iscat = [x in categorical_feat for x in feature_names]\n",
        "\n",
        "X_train, X_test, y_train, y_test = train_test_split(X_feat, y, test_size=0.30, random_state=42)\n",
        "\n",
        "# Add 100 duplicates\n",
        "X_dup = X_test.sample(n=100,replace=False, ignore_index=False)\n",
        "y_dup = y_test[X_dup.index]\n",
        "X_test = pd.concat([X_test,X_dup], axis=0)\n",
        "y_test = pd.concat([y_test,y_dup], axis=0)\n",
        "\n",
        "X_test_nd = X_test.to_numpy(dtype=float)\n",
        "X_train_nd = X_train.to_numpy(dtype=float)\n",
        "\n",
        "# Normalizing data ********* <<< DEBUG, remove next rows if you don't want normalization >>> **********\n",
        "scaler = MinMaxScaler()\n",
        "\n",
        "X_train_nd = scaler.fit_transform(X_train_nd)\n",
        "X_train = pd.DataFrame(X_train_nd, columns=feature_names)\n",
        "\n",
        "X_test_nd = scaler.fit_transform(X_test_nd)\n",
        "X_test = pd.DataFrame(X_test_nd, columns=feature_names)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "srARwMTh4Rf1"
      },
      "source": [
        "### HEPAR Dataset"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UdTHIz4-4Yxm"
      },
      "outputs": [],
      "source": [
        "df_raw_iid = pd.read_csv(io.BytesIO(uploaded['HEPAR_simulated_patients.csv']))\n",
        "df_raw_ood = pd.read_csv(io.BytesIO(uploaded['HEPAR_simulated_patients_ood.tsv']))\n",
        "df_raw_iid.drop(\"Unnamed: 0\", axis=1, inplace=True)\n",
        "df_raw_ood.drop(\"Unnamed: 0\", axis=1, inplace=True)\n",
        "display(df_raw_iid)\n",
        "display(df_raw_ood)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "P0AFCH5zgZSn"
      },
      "outputs": [],
      "source": [
        "def preproc_wrapper(df):\n",
        "  # Categorical variables handled\n",
        "  df_num = df.copy()\n",
        "  \"\"\"Handling categorical and nominal data\"\"\"\n",
        "\n",
        "  nominal_feat = ['age', 'triglycerides', 'bilirubin', 'phosphatase', 'proteins', 'platelet', \n",
        "                  'inr', 'urea', 'ESR', 'alt', 'ast','amylase', 'ggtp', 'cholesterol', 'albumin']\n",
        "\n",
        "  # preprocessing: data transformation\n",
        "  nominal_dict = {}\n",
        "  for n in nominal_feat:\n",
        "    unique_val = df[n].unique().tolist()\n",
        "    unique_num = [int(x.split('_')[1]) for x in unique_val]\n",
        "    val2num = dict(zip(unique_val, unique_num))\n",
        "    num2cat = dict(zip(sorted(unique_num), range(1, len(unique_num)+1)))\n",
        "    dict_n = {}\n",
        "    for k,v in val2num.items():\n",
        "      dict_n[k] = num2cat[v]\n",
        "    nominal_dict[n] = dict_n\n",
        "  # print(nominal_dict) ***DEBUG to check that it is consistent between the two datasets (it is, but it might not if using different datasets --> do not copy paste)\n",
        "\n",
        "  special_feat = ['ChHepatitis', 'sex', 'Cirrhosis']\n",
        "\n",
        "  dict_chhepa = {'absent':0, 'active':1, 'persistent':2}\n",
        "  dict_sex = {'female':1, 'male':2}\n",
        "  dict_cirr = {'absent':0, 'decompensate':1, 'compensate':2}\n",
        "\n",
        "  df_num['sex'] = [dict_sex[x] for x in df['sex']]\n",
        "  df_num['ChHepatitis'] = [dict_chhepa[x] for x in df['ChHepatitis']]\n",
        "  df_num['Cirrhosis'] = [dict_cirr[x] for x in df['Cirrhosis']]\n",
        "\n",
        "  categorical_feat = [x for x in df.columns if x not in nominal_feat+special_feat]\n",
        "  dict_cat = {'absent':0, 'present':1}\n",
        "\n",
        "  #print(categorical_feat)\n",
        "  for c in categorical_feat:\n",
        "    #print(c)\n",
        "    df_num[c] = [dict_cat[x] for x in df[c]]\n",
        "\n",
        "  for n in nominal_feat:\n",
        "    newcol = [nominal_dict[n][x] for x in df[n]]\n",
        "    df_num[n] = newcol\n",
        "\n",
        "  pred_col = 'hospital'\n",
        "\n",
        "  df_num[pred_col].value_counts()\n",
        "\n",
        "  np.random.seed(1)\n",
        "\n",
        "  dataset_class_1 = df_num[df_num[pred_col]==1]\n",
        "  dataset_class_0 = df_num[df_num[pred_col]==0]\n",
        "\n",
        "  # selecting a subpopulation of X% that will be used for training and testing ***DEBUG\n",
        "  # n_subpop = int(0.5*df_num.shape[0])\n",
        "  # i_class_yes = np.random.randint(0, high=dataset_class_1.shape[0], size=int(n_subpop/2))\n",
        "  # i_class_no = np.random.randint(0, high=dataset_class_0.shape[0], size=int(n_subpop/2))\n",
        "  # X = dataset_class_1.iloc[i_class_yes].append(dataset_class_0.iloc[i_class_no])\n",
        "\n",
        "  X = df_num\n",
        "  y = X[pred_col]\n",
        "  X_feat = X.drop(columns=[pred_col])\n",
        "  feature_names = X_feat.columns\n",
        "  iscat = [x in categorical_feat for x in feature_names]\n",
        "\n",
        "  return X_feat, y, feature_names, iscat\n",
        "\n",
        "X_iid, y_iid, feat_names_iid, iscat_iid = preproc_wrapper(df_raw_iid)\n",
        "X_ood, y_ood, feat_names_ood, iscat_ood = preproc_wrapper(df_raw_ood)\n",
        "\n",
        "feature_names = feat_names_iid # = feat_names_ood\n",
        "iscat = iscat_iid # = iscat_ood"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iCNUSywBwC83"
      },
      "outputs": [],
      "source": [
        "X_train, X_test, y_train, y_test = train_test_split(X_iid, y_iid, test_size=0.05, random_state=6)\n",
        "#to test on iid\n",
        "# X_test = X_ood\n",
        "# y_test = y_ood\n",
        "\n",
        "# Add 100 duplicates\n",
        "X_dup = X_test.sample(n=100,replace=False, ignore_index=False)\n",
        "y_dup = y_test[X_dup.index]\n",
        "X_test = pd.concat([X_test,X_dup], axis=0)\n",
        "y_test = pd.concat([y_test,y_dup], axis=0)\n",
        "\n",
        "X_test_nd = X_test.to_numpy(dtype=float)\n",
        "X_train_nd = X_train.to_numpy(dtype=float)\n",
        "\n",
        "# Normalizing data ********* <<< DEBUG, remove next rows if you don't want normalization >>> **********\n",
        "#scaler = MinMaxScaler()\n",
        "\n",
        "#X_train_nd = scaler.fit_transform(X_train_nd)\n",
        "X_train = pd.DataFrame(X_train_nd, columns=feature_names)\n",
        "\n",
        "#X_test_nd = scaler.fit_transform(X_test_nd)\n",
        "X_test = pd.DataFrame(X_test_nd, columns=feature_names)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NqUWeYBzJMoy"
      },
      "source": [
        "## **Classifiers**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tTlhyS0zMzD8"
      },
      "outputs": [],
      "source": [
        "#CLASSIFIERS\n",
        "def get_best_threshold_f1(X,y,y_pred):\n",
        "  prec, recall, threshold = precision_recall_curve(y, y_pred)\n",
        "  f1 = 2*prec*recall/(prec+recall)\n",
        "  #return(threshold[np.where(f1==np.nanmax(f1))])\n",
        "  return(0.5)\n",
        "\n",
        "def print_metrics(y_true,y_pred):\n",
        "  tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()\n",
        "  specificity = tn / (tn+fp)\n",
        "  print(f\"- Acc: {accuracy_score(y_true,y_pred)}\\n- Rec: {recall_score(y_true,y_pred)}\\n- Prec: {precision_score(y_true,y_pred)}\\n- Spec: {specificity}\\n- F1: {f1_score(y_true,y_pred)}\")\n",
        "\n",
        "print(\"True Labels\")\n",
        "print(collections.Counter(y_test))\n",
        "print()\n",
        "\n",
        "#Gradient Boosting\n",
        "clf_gb = GradientBoostingClassifier(n_estimators=100, learning_rate=1.0, max_depth=1, random_state=0).fit(X_train, y_train)\n",
        "ypred_proba = clf_gb.predict_proba(X_test)[:,1]\n",
        "sel_thr_gb = get_best_threshold_f1(X_test,y_test,ypred_proba)\n",
        "ypred_gb = (clf_gb.predict_proba(X_test)[:,1] >= sel_thr_gb).astype(int)\n",
        "\n",
        "print(\"\\nGradient Boosting\")\n",
        "print(collections.Counter(ypred_gb))\n",
        "print_metrics(y_test, ypred_gb)\n",
        "\n",
        "#Random Forest\n",
        "clf_rf = RandomForestClassifier(random_state=1).fit(X_train, y_train)\n",
        "ypred_proba = clf_rf.predict_proba(X_test)[:,1]\n",
        "sel_thr_rf = get_best_threshold_f1(X_test,y_test,ypred_proba)\n",
        "ypred_rf = (clf_rf.predict_proba(X_test)[:,1] >= sel_thr_rf).astype(int)\n",
        "\n",
        "print(\"\\nRandom Forest\")\n",
        "print(collections.Counter(ypred_rf))\n",
        "print_metrics(y_test, ypred_rf)\n",
        "\n",
        "#Logistic Regression\n",
        "clf_lr = LogisticRegression(random_state=1, penalty='l1', solver='liblinear').fit(X_train, y_train)\n",
        "ypred_proba = clf_lr.predict_proba(X_test)[:,1]\n",
        "sel_thr_lr = get_best_threshold_f1(X_test,y_test,ypred_proba)\n",
        "ypred_lr = (clf_lr.predict_proba(X_test)[:,1] >= sel_thr_lr).astype(int)\n",
        "\n",
        "print(\"\\nLogistic Regression\")\n",
        "print(collections.Counter(ypred_lr))\n",
        "print_metrics(y_test, ypred_lr)\n",
        "\n",
        "#Multi-layer perceptron (neural network)\n",
        "clf_nn = MLPClassifier(random_state=1,hidden_layer_sizes=(100, 50, 20), max_iter=500).fit(X_train, y_train)\n",
        "ypred_proba = clf_nn.predict_proba(X_test)[:,1]\n",
        "sel_thr_nn = get_best_threshold_f1(X_test,y_test,ypred_proba)\n",
        "ypred_nn = (clf_nn.predict_proba(X_test)[:,1] >= sel_thr_nn).astype(int)\n",
        "\n",
        "print(\"\\nMulti-layer Perceptron\")\n",
        "print(collections.Counter(ypred_nn))\n",
        "print_metrics(y_test, ypred_nn)\n",
        "\n",
        "classifier_dict = {'NN': clf_nn, 'GB': clf_gb, 'LR': clf_lr, 'RF': clf_rf}\n",
        "classifier_thr = {'NN': sel_thr_nn, 'GB': sel_thr_gb, 'LR': sel_thr_lr, 'RF': sel_thr_rf}\n",
        "\n",
        "print(\"\\nThresholds\")\n",
        "print(classifier_thr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OWkLz3ZP-qq3"
      },
      "source": [
        "## **Explainers: Overall Metrics**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xSFr8WHHHKgz"
      },
      "source": [
        "### **Araucana**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "n60eAV9Jh_XP"
      },
      "outputs": [],
      "source": [
        "#ARAUCANA EXPLAINER\n",
        "\n",
        "def run_auracana(X_test, feature_names, cat_list, X_train, nsize, oversampling, oversampling_size, classifiers):\n",
        "    xai_dict = dict(zip(list(classifiers.keys()), [[], [], [], []]))\n",
        "    fidelity_dict = {}\n",
        "    times = {}\n",
        "    # for each classifier\n",
        "    for k_clf, v_clf in classifiers.items():\n",
        "        xai_start = time.time()\n",
        "        fidelity_list = []\n",
        "        fidelity_fail_count = 0\n",
        "        # for each instance of the test set\n",
        "        for i in tqdm(range(X_test.shape[0])):\n",
        "            # declare the instance from the test set we want to explain\n",
        "            instance = X_test[i, :].reshape(1, X_test.shape[1])\n",
        "            # save the predicted label for the instance according to the current classifier\n",
        "            instance_pred_y = (v_clf.predict_proba(instance)[:,1] >= classifier_thr[k_clf]).astype(int)\n",
        "            # build xai tree to explain the instance classification\n",
        "            try:\n",
        "                xai_tree = araucanaxai.run(x_target=instance, y_pred_target=instance_pred_y,\n",
        "                                           x_train=X_train, feature_names=feature_names,\n",
        "                                           cat_list=cat_list, oversampling=oversampling, oversampling_size = oversampling_size,\n",
        "                                           predict_fun= lambda instance:(v_clf.predict_proba(instance)[:,1] >= classifier_thr[k_clf]).astype(int), \n",
        "                                           neighbourhood_size=nsize)\n",
        "                xai_dict[k_clf].append(xai_tree)\n",
        "                #fidelity online check \n",
        "                if xai_tree['tree'].predict(instance) != instance_pred_y:\n",
        "                  fidelity_fail_count += 1\n",
        "                  if fidelity_fail_count == 100:\n",
        "                    warnings.warn(\"Warning...........divergent predictions between explainer and original model\")\n",
        "                fidelity_list.append((xai_tree['tree'].predict(instance) == instance_pred_y).astype(int))\n",
        "            except Exception as e:\n",
        "                print(e)\n",
        "                xai_dict[k_clf].append(None)\n",
        "                fidelity_list.append(None)\n",
        "                continue\n",
        "        if len([x for x in fidelity_list if x != None]) != 0:\n",
        "            fidelity_dict[k_clf] = sum([x for x in fidelity_list if x != None]) / len([x for x in fidelity_list if x != None])\n",
        "        else:\n",
        "            fidelity_dict[k_clf] = None\n",
        "        times[k_clf]=time.time() - xai_start\n",
        "    return {\"fidelity\": fidelity_dict, \"test\": xai_dict, \"time\": times}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "n-izBuE-kVOo"
      },
      "outputs": [],
      "source": [
        "# Function to compare explaination\n",
        "def exp_similarity_trees(tree_a, tree_b):\n",
        "    feat_a = tree_a.tree_.feature.copy()\n",
        "    feat_b = tree_b.tree_.feature.copy()\n",
        "    n = min(len(feat_a), len(feat_b))       # if different len, trim the longer\n",
        "    feat_a = feat_a[0:n]\n",
        "    feat_b = feat_b[0:n]\n",
        "    return sum((feat_a == feat_b).astype(int)) / n\n",
        "\n",
        "# FIDELITY: concordance of the predictions between the XAI proxy model and the complex model\n",
        "\n",
        "def check_fidelity(araucana_xai):\n",
        "  for c in classifier_dict.keys():\n",
        "    print('- Classifier ', c, ': fidelity = ', araucana_xai[\"fidelity\"][c])\n",
        "\n",
        "# IDENTITY: if there are 2 identical instances, they must have identical explanations\n",
        "def get_equal_instances(X_test):\n",
        "  iDup = np.where(X_test.duplicated(keep=False))[0]  # instances in test which are duplicated\n",
        "  xai_dup = []\n",
        "  dup = []\n",
        "  for couple in list(combinations(iDup, 2)):\n",
        "    check = (X_test.iloc[couple[0],] == X_test.iloc[couple[1],]).all()\n",
        "    if check:\n",
        "      dup.append(couple)\n",
        "  return(dup)\n",
        "\n",
        "def check_identity(X_test, araucana_xai):\n",
        "  dup = get_equal_instances(X_test)\n",
        "  for c in classifier_dict.keys():\n",
        "    if len(dup)!=0:\n",
        "      right = 0\n",
        "      for d in dup:\n",
        "        xai1 = araucana_xai[\"test\"][c][d[0]]['tree']\n",
        "        xai2 = araucana_xai[\"test\"][c][d[1]]['tree']\n",
        "        s = exp_similarity_trees(xai1, xai2)\n",
        "        if s==1:\n",
        "          right += 1\n",
        "      score = right/len(dup)\n",
        "      print('- Classifier', c, ': identity = ', score)\n",
        "    else:\n",
        "      print('- Classifier', c, ': identity = na')\n",
        "\n",
        "# SEPARABILITY: if there are 2 dissimilar instances, they must have dissimilar explanations\n",
        "\n",
        "def check_separability(X_test, araucana_xai):\n",
        "  dup = get_equal_instances(X_test)\n",
        "  iNotDup = [i for i in range(len(X_test)) if i not in [dup[i][1] for i in range(len(dup))]] # instances in test which are NOT duplicated\n",
        "  combs = list(combinations(iNotDup, 2))\n",
        "  for c in classifier_dict.keys():\n",
        "    wrong = 0\n",
        "    for couple in combs:\n",
        "        xai1 = araucana_xai[\"test\"][c][couple[0]]['tree']\n",
        "        xai2 = araucana_xai[\"test\"][c][couple[1]]['tree']\n",
        "        s = exp_similarity_trees(xai1, xai2)\n",
        "        if s==1:\n",
        "            wrong += 1\n",
        "    total = len(combs)\n",
        "    score = wrong/total\n",
        "    print('- Classifier', c, ': separability = ', score, \"(\",wrong, 'equals out of',total,\")\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Fwtv6yVII4bD"
      },
      "source": [
        "#### Single Run"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zzrF2zieUmlK"
      },
      "outputs": [],
      "source": [
        "os_method = None # Possible values = [None,\"smote\", \"uniform\", \"non-uniform\"]\n",
        "\n",
        "#following parameters are not used if os_method = None\n",
        "os_size = 200 # number of new samples to be generated with oversampling\n",
        "nsize = 200\n",
        "\n",
        "araucana_xai = run_auracana(X_test=X_test_nd, cat_list=iscat, feature_names=feature_names,\n",
        "                            X_train=X_train_nd, nsize=nsize, oversampling=os_method, oversampling_size = os_size,\n",
        "                            classifiers = classifier_dict)\n",
        "\n",
        "print('\\ntime')\n",
        "for c in classifier_dict.keys():\n",
        "  print('- Classifier ', c, ': time = ', araucana_xai[\"time\"][c])\n",
        "\n",
        "print('\\nfidelity')\n",
        "check_fidelity(araucana_xai)\n",
        "\n",
        "print('\\nidentity')\n",
        "check_identity(X_test, araucana_xai)\n",
        "\n",
        "print('\\nseparability')\n",
        "check_separability(X_test, araucana_xai)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gdgjpAJ1y2fp"
      },
      "source": [
        "#### Hyperaparameter Search"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RHmZoAJfXh4L"
      },
      "source": [
        "Oversampling Methods"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3nw76KZny117"
      },
      "outputs": [],
      "source": [
        "os_methods = [\"smote\", \"uniform\", \"non-uniform\"]\n",
        "os_size = 100 # number of new samples to be generated with oversampling\n",
        "nsize_dict = {\"smote\": 800, \"uniform\":35, \"non-uniform\":35}\n",
        "\n",
        "for os_method in os_methods:\n",
        "  print(str(os_method))\n",
        "  print(nsize_dict[os_method])\n",
        "  araucana_xai = run_auracana(X_test=X_test_nd, cat_list=iscat, feature_names=feature_names,\n",
        "                              X_train=X_train_nd, nsize=nsize_dict[os_method], oversampling=os_method, oversampling_size = os_size,\n",
        "                              classifiers = classifier_dict)\n",
        "  # TIME: \n",
        "  print('\\ntime')\n",
        "  for c in classifier_dict.keys():\n",
        "    print('- Classifier ', c, ': time = ', araucana_xai[\"time\"][c])\n",
        "\n",
        "  # FIDELITY: concordance of the predictions between the XAI proxy model and the complex model\n",
        "  print('\\nfidelity')\n",
        "  for c in classifier_dict.keys():\n",
        "    print('- Classifier ', c, ': fidelity = ', araucana_xai[\"fidelity\"][c])\n",
        "\n",
        "  # IDENTITY: if there are 2 identical instances, they must have identical explanations\n",
        "  print('\\nidentity')\n",
        "  check_identity(X_test, araucana_xai)\n",
        "\n",
        "  # SEPARABILITY: if there are 2 dissimilar instances, they must have dissimilar explanations\n",
        "  print('\\nseparability')\n",
        "  check_separability(X_test_nd, araucana_xai)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4VcQkQrJX5VB"
      },
      "source": [
        "Oversampling Neighborhood Size"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Pl0CkZFDX4Wb"
      },
      "outputs": [],
      "source": [
        "nratios = [0.25, 0.10, 0.05, 0.01, 0.005] # neighborhood size expressed as a portion of training size\n",
        "os_method = \"non-uniform\" # Possible values = [None,\"smote\", \"uniform\", \"non-uniform\"]\n",
        "os_size = 50 # number of new samples to be generated with oversampling\n",
        "\n",
        "for nratio in nratios:\n",
        "  nsize = int(X_train.shape[0] * nratio)\n",
        "  print(nsize)\n",
        "  araucana_xai = run_auracana(X_test=X_test_nd, cat_list=iscat, feature_names=feature_names,\n",
        "                              X_train=X_train_nd, nsize=nsize, oversampling=os_method, oversampling_size = os_size,\n",
        "                              classifiers = classifier_dict)\n",
        "  # TIME: \n",
        "  print('\\ntime')\n",
        "  for c in classifier_dict.keys():\n",
        "    print('- Classifier ', c, ': time = ', araucana_xai[\"time\"][c])\n",
        "\n",
        "  # FIDELITY: concordance of the predictions between the XAI proxy model and the complex model\n",
        "  print('\\nfidelity')\n",
        "  for c in classifier_dict.keys():\n",
        "    print('- Classifier ', c, ': fidelity = ', araucana_xai[\"fidelity\"][c])\n",
        "\n",
        "  # IDENTITY: if there are 2 identical instances, they must have identical explanations\n",
        "  print('\\nidentity')\n",
        "  check_identity(X_test, araucana_xai)\n",
        "\n",
        "  # SEPARABILITY: if there are 2 dissimilar instances, they must have dissimilar explanations\n",
        "  print('\\nseparability')\n",
        "  check_separability(X_test_nd, araucana_xai)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QKuVT8I6jK6b"
      },
      "source": [
        "Oversampling Sample Size"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FUKG406jjKda"
      },
      "outputs": [],
      "source": [
        "nratio = 0.01 # neighborhood size expressed as a portion of training size\n",
        "os_method = \"non-uniform\" # Possible values = [None,\"smote\", \"uniform\", \"non-uniform\"]\n",
        "os_sizes = [1, 10, 50, 100, 500] # number of new samples to be generated with oversampling\n",
        "nsize = int(X_train.shape[0] * nratio)\n",
        "\n",
        "for os_size in os_sizes:\n",
        "  print(os_size)\n",
        "  araucana_xai = run_auracana(X_test=X_test_nd, cat_list=iscat, feature_names=feature_names,\n",
        "                              X_train=X_train_nd, nsize=nsize, oversampling=os_method, oversampling_size = os_size,\n",
        "                              classifiers = classifier_dict)\n",
        "  # TIME: \n",
        "  print('\\ntime')\n",
        "  for c in classifier_dict.keys():\n",
        "    print('- Classifier ', c, ': time = ', araucana_xai[\"time\"][c])\n",
        "\n",
        "  # FIDELITY: concordance of the predictions between the XAI proxy model and the complex model\n",
        "  print('\\nfidelity')\n",
        "  for c in classifier_dict.keys():\n",
        "    print('- Classifier ', c, ': fidelity = ', araucana_xai[\"fidelity\"][c])\n",
        "\n",
        "  # IDENTITY: if there are 2 identical instances, they must have identical explanations\n",
        "  print('\\nidentity')\n",
        "  check_identity(X_test, araucana_xai)\n",
        "\n",
        "  # SEPARABILITY: if there are 2 dissimilar instances, they must have dissimilar explanations\n",
        "  print('\\nseparability')\n",
        "  check_separability(X_test_nd, araucana_xai)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K7Z8hFtSHg5z"
      },
      "source": [
        "### **SHAP**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dCfxCoFXNDy1"
      },
      "outputs": [],
      "source": [
        "#SHAP EXPLAINER\n",
        "\n",
        "def run_shap(X_test, classifiers):\n",
        "    xai_dict = dict(zip(list(classifiers.keys()), [[], [], []])) \n",
        "    fidelity_dict = {}\n",
        "    times = {}\n",
        "    # for each classifier\n",
        "    for k, v in classifiers.items():\n",
        "      start = time.time()\n",
        "      #For Logistic Regression, we need to use shap.Explainer\n",
        "      if k=='LR' or k=='NN':\n",
        "        explainer = shap.Explainer(v.predict_proba, np.median(X_test, axis=0).reshape(1, X_test.shape[1]))\n",
        "        shap_val = explainer(X_test) \n",
        "        shap_values = shap_val.values\n",
        "      #For Random Forest and Gradient Boosting, we can use shap.TreeExplainer\n",
        "      else:\n",
        "        explainer = shap.TreeExplainer(v)\n",
        "        shap_values = explainer.shap_values(X_test)\n",
        "      if k=='RF':\n",
        "        shap_values = shap_values[1]\n",
        "      xai_dict[k] = shap_values\n",
        "      times[k] = time.time()-start\n",
        "      fidelity_dict[k] = 1 #by default, not a local explainer\n",
        "    return {\"fidelity\": fidelity_dict, \"test\": xai_dict, \"time\": times}\n",
        "\n",
        "shap_xai = run_shap(X_test=X_test_nd, classifiers = classifier_dict)\n",
        "\n",
        "print('\\ntime')\n",
        "for c in classifier_dict.keys():\n",
        "  print('- Classifier ', c, ': time = ', shap_xai[\"time\"][c])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uxA6S5W3jfOw"
      },
      "outputs": [],
      "source": [
        "for c in classifier_dict.keys():\n",
        "  print('- Classifier ', c, ': time = ', shap_xai[\"time\"][c])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0kRZBAaWLx5Z"
      },
      "outputs": [],
      "source": [
        "# Function to compare explaination (oversimplified, we only need to check if they are equal or not at the moment)\n",
        "def exp_similarity_shap(shap_vals_a, shap_vals_b):\n",
        "    #EXACT MATCH VALUE BASED\n",
        "    #is_equal = np.array_equal(shap_vals_a,shap_vals_b)\n",
        "    #EXACT MATCH FEATURE RANKING\n",
        "    if len(shap_vals_a.shape)>1:\n",
        "      shap_vals_a = shap_vals_a[:,1]\n",
        "    z1 = zip(feature_names, abs(shap_vals_a))\n",
        "    z1 = list(z1)\n",
        "    z1 = sorted(z1, key = lambda x: x[1])\n",
        "    z1 = [z1[i][0] for i in range(len(z1))]\n",
        "\n",
        "    if len(shap_vals_b.shape)>1:\n",
        "      shap_vals_b = shap_vals_b[:,1]\n",
        "    z2 = zip(feature_names, abs(shap_vals_b))\n",
        "    z2 = list(z2)\n",
        "    z2 = sorted(z2, key = lambda x: x[1])\n",
        "    z2 = [z2[i][0] for i in range(len(z2))]\n",
        "    return 1 if z1==z2 else 0\n",
        "\n",
        "# FIDELITY (1 by default, not a local explainer)\n",
        "\n",
        "print('\\nfidelity')\n",
        "for c in classifier_dict.keys():\n",
        "  print('- Classifier ', c, ': fidelity = ', shap_xai[\"fidelity\"][c])\n",
        "\n",
        "# IDENTITY\n",
        "  \n",
        "print('\\nidentity')\n",
        "iDup = np.where(X_test.duplicated(keep=False))[0]  # instances in test which are duplicated\n",
        "xai_dup = []\n",
        "dup = []\n",
        "for couple in list(combinations(iDup, 2)):\n",
        "  check = (X_test.iloc[couple[0],] == X_test.iloc[couple[1],]).all()\n",
        "  if check:\n",
        "    dup.append(couple)\n",
        "\n",
        "for c in classifier_dict.keys():\n",
        "  if len(dup)!=0:\n",
        "    right = 0\n",
        "    for d in dup:\n",
        "      s = exp_similarity_shap(shap_xai[\"test\"][c][d[0]], shap_xai[\"test\"][c][d[1]])\n",
        "      if s==1:\n",
        "        right += 1\n",
        "    score = right/len(dup)\n",
        "    print('- Classifier', c, ': identity = ', score)\n",
        "  else:\n",
        "   print('- Classifier', c, ': identity = na')\n",
        "\n",
        "# SEPARABILITY\n",
        "\n",
        "print('\\nseparability')\n",
        "dup = get_equal_instances(X_test)\n",
        "iNotDup = [i for i in range(len(X_test)) if i not in [dup[i][1] for i in range(len(dup))]] # instances in test which are NOT duplicated\n",
        "combs = list(combinations(iNotDup, 2))\n",
        "for c in classifier_dict.keys():\n",
        "  wrong = 0\n",
        "  for couple in combs:\n",
        "      s = exp_similarity_shap(shap_xai[\"test\"][c][couple[0]], shap_xai[\"test\"][c][couple[1]])\n",
        "      if s==1:\n",
        "          wrong += 1\n",
        "  total = len(combs)\n",
        "  score = wrong/total\n",
        "  print('- Classifier', c, ': separability = ', score, \"(\",wrong, 'equals out of',total,\")\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5Nuy7-YZQJn3"
      },
      "source": [
        "### **LIME**"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def return_weights(exp):\n",
        "   \n",
        "    exp_list = exp.as_map()[1]\n",
        "    #print(exp_list) = [(7, -0.539), (58, 0.215), etc]\n",
        "    exp_list = sorted(exp_list, key=lambda x: x[0]) \n",
        "    #print(exp_list) = [(0, 0.014), (1, -0.061), (2, 0.046), etc]\n",
        "    #sort so that we can compare with feature list (0 = feature at position 0 in the feature list)\n",
        "    exp_weight = [x[1] for x in exp_list]\n",
        "    \n",
        "    return exp_weight\n",
        "\n",
        "explainer = lime.lime_tabular.LimeTabularExplainer(X_train_nd, \n",
        "                                                    feature_names=feature_names, \n",
        "                                                    verbose=False, \n",
        "                                                    mode='classification')\n",
        "\n",
        "for k, v in classifier_dict.items():\n",
        "  print(f\"\\n\\nClassifier: {k}\")\n",
        "  trs = classifier_thr[k]\n",
        "\n",
        "  weights = []\n",
        "  fidelity = 0\n",
        "  start = time.time()\n",
        "  #Iterate over rows in feature matrix\n",
        "  for i in tqdm(range(len(X_test_nd))):\n",
        "      \n",
        "      #Get explanation for ith row\n",
        "      exp = explainer.explain_instance(X_test_nd[i, :], \n",
        "                                      v.predict_proba,\n",
        "                                      num_features = X_test_nd.shape[1],\n",
        "                                      num_samples = 200)\n",
        "      #Fidelity check\n",
        "      instance = X_test_nd[i, :].reshape(1, X_test_nd.shape[1])\n",
        "      instance_pred_y = int((v.predict_proba(instance)[:,1] >= trs))\n",
        "\n",
        "      ###DEBUG\n",
        "      if int((exp.local_pred >= 0.5)) == instance_pred_y:\n",
        "        fidelity +=1  \n",
        "\n",
        "      # if (exp.predict_proba[1] >= 0.5).astype(int) == instance_pred_y:\n",
        "      #   fidelity +=1\n",
        "\n",
        "      #Get weights\n",
        "      exp_weight = return_weights(exp)\n",
        "      weights.append(exp_weight)\n",
        "  print(f\"\\nTime = {time.time()-start}\")\n",
        "      \n",
        "  #Create DataFrame\n",
        "  lime_weights = pd.DataFrame(data=weights,columns=feature_names)\n",
        "  #display(lime_weights)\n",
        "\n",
        "  #Fidelity\n",
        "  print(f\"\\nFidelity = {fidelity/len(X_test_nd)}\")\n",
        "\n",
        "  #Identity\n",
        "  iDup = np.where(X_test.duplicated(keep=False))[0]  # instances in test which are duplicated\n",
        "  couples = list(combinations(iDup, 2))\n",
        "  equal_count = 0 \n",
        "  for couple in tqdm(couples):\n",
        "    check = (X_test.iloc[couple[0],] == X_test.iloc[couple[1],]).all()\n",
        "    if check:\n",
        "      exp_names1 = abs(lime_weights.iloc[couple[0]]).sort_values(ascending=False).index.tolist()\n",
        "      exp_names2 = abs(lime_weights.iloc[couple[1]]).sort_values(ascending=False).index.tolist()\n",
        "      if exp_names1==exp_names2:\n",
        "        equal_count += 1\n",
        "  if len(couples)==0:\n",
        "    print(\"\\nIdentity = n.a.\")\n",
        "  else: \n",
        "    print(f\"\\nIdentity = {equal_count/len(couples)}\")\n",
        "\n",
        "  #Separability\n",
        "  dup = get_equal_instances(X_test)\n",
        "  iNotDup = [i for i in range(len(X_test)) if i not in [dup[i][1] for i in range(len(dup))]] # instances in test which are NOT duplicated\n",
        "  couples = list(combinations(iNotDup, 2))\n",
        "  equal_count = 0\n",
        "  for couple in tqdm(couples):\n",
        "    exp_names1 = abs(lime_weights.iloc[couple[0]]).sort_values(ascending=False).index.tolist()\n",
        "    exp_names2 = abs(lime_weights.iloc[couple[1]]).sort_values(ascending=False).index.tolist()\n",
        "    if exp_names1==exp_names2:\n",
        "      equal_count += 1\n",
        "\n",
        "  print(f\"\\nSeparability = {equal_count/len(couples)}\")"
      ],
      "metadata": {
        "id": "vyM0ZZQUj2nv"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3jGF2tS02v6e"
      },
      "source": [
        "## Explainers: Single Instance (i.e. Local) Explanations"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Md75F3KVLBWl"
      },
      "outputs": [],
      "source": [
        "# declare the classifier\n",
        "k_clf = \"NN\"\n",
        "v_clf = classifier_dict[k_clf]\n",
        "# it could be usefult to separate right predictions and wrong predictions\n",
        "right = np.where((y_test-(v_clf.predict_proba(X_test)[:,1] >= classifier_thr[k_clf]).astype(int) == 0))[0]\n",
        "wrong = np.where((y_test-(v_clf.predict_proba(X_test)[:,1] >= classifier_thr[k_clf]).astype(int) != 0) & (y_test==0))[0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "u8730lIlBKAy"
      },
      "outputs": [],
      "source": [
        "# declare the instance from the test set we want to explain\n",
        "ind = wrong[0]\n",
        "instance = X_test_nd[ind, :].reshape(1, X_test.shape[1])\n",
        "\n",
        "# save the predicted label for the instance according to the current classifier\n",
        "instance_pred_y = (v_clf.predict_proba(instance)[:,1] >= classifier_thr[k_clf]).astype(int)\n",
        "\n",
        "# cumbersome line to print the example nicely\n",
        "print(pd.DataFrame(data=[np.concatenate([instance.flatten(),instance_pred_y],dtype=object)],\n",
        "                     columns=np.concatenate([feature_names,[\"pred_label\"]],dtype=object)))\n",
        "\n",
        "print(f\"Predicted class:{instance_pred_y}\")\n",
        "print(f\"True class:{y_test.iloc[ind]}\")\n",
        "\n",
        "## ARAUCANA\n",
        "print(\"\\n\\n ==================== ARAUCANA ====================\\n\\n\")\n",
        "\n",
        "# build xai tree to explain the instance classification\n",
        "xai_tree = araucanaxai.run(x_target=instance, y_pred_target=instance_pred_y,\n",
        "                           x_train=X_train_nd, feature_names=feature_names,\n",
        "                           cat_list=iscat, oversampling=None, max_depth=4, neighbourhood_size=300,\n",
        "                           predict_fun= v_clf.predict)\n",
        "from sklearn import tree\n",
        "fig, ax = plt.subplots(figsize=(10, 10))\n",
        "tree.plot_tree(xai_tree['tree'], feature_names=feature_names, filled=True, class_names=[\"None\",\"Death\"])\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "## SHAP\n",
        "print(\"\\n\\n ==================== SHAP ====================\\n\\n\")\n",
        "\n",
        "explainer = shap.Explainer(v_clf.predict_proba, np.median(X_test, axis=0).reshape(1, X_test.shape[1]))\n",
        "shap_val = explainer(instance) \n",
        "base_values = shap_val.base_values\n",
        "shap_values = shap_val.values\n",
        "\n",
        "labels = feature_names\n",
        "x = np.arange(len(labels))\n",
        "w=0.9\n",
        "fig, ax = plt.subplots(figsize=(8,10))\n",
        "vals = shap_val[0].values[:,1]\n",
        "bars = plt.barh(x, vals, label='SHAP', height=w, color= ['#449F89' if i else '#CA7952' for i in vals > 0])\n",
        "plt.axvline(x=0, ymin=0, ymax=1, color='grey')\n",
        "ax.set_yticks(x)\n",
        "ax.set_yticklabels(labels, fontsize=12)\n",
        "plt.title('SHAP values')\n",
        "plt.show()\n",
        "\n",
        "## LIME\n",
        "print(\"\\n\\n ==================== LIME ====================\\n\\n\")\n",
        "\n",
        "explainer_lime = lime.lime_tabular.LimeTabularExplainer(X_train_nd, feature_names = feature_names)\n",
        "exp = explainer_lime.explain_instance(data_row=X_test_nd[ind, :], \n",
        "                                      predict_fn=v_clf.predict_proba,\n",
        "                                      num_features=X_test_nd.shape[1],                \n",
        "                                      top_labels=1)\n",
        "\n",
        "exp_list = exp.as_list()\n",
        "vals = [x[1] for x in exp_list]\n",
        "names = [x[0] for x in exp_list]\n",
        "vals.reverse()\n",
        "names.reverse()\n",
        "colors = ['#4986E7' if x > 0 else '#D93025' for x in vals]\n",
        "x = np.arange(len(labels))\n",
        "w=0.9\n",
        "fig, ax = plt.subplots(figsize=(8,10))\n",
        "bars = plt.barh(x, vals, label='LIME', height=w, color= colors)\n",
        "plt.axvline(x=0, ymin=0, ymax=1, color='grey')\n",
        "ax.set_yticks(x)\n",
        "ax.set_yticklabels(names, fontsize=12)\n",
        "plt.title('LIME')\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "gdgjpAJ1y2fp"
      ],
      "name": "araucana_XAI_special_issue.ipynb",
      "toc_visible": true,
      "provenance": [],
      "authorship_tag": "ABX9TyPC15e6xcQZpKDYGOeZ1QmF",
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}