examples/cem_mnist.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Contrastive Explanations Method (CEM) applied to MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Contrastive Explanation Method (CEM) can generate black box model explanations in terms of pertinent positives (PP) and pertinent negatives (PN). For PP, it finds what should be minimally and sufficiently present (e.g. important pixels in an image) to justify its classification. PN on the other hand identify what should be minimally and necessarily absent from the explained instance in order to maintain the original prediction.\n",
    "\n",
    "The original paper where the algorithm is based on can be found on [arXiv](https://arxiv.org/pdf/1802.07623.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "import keras\n",
    "from keras import backend as K\n",
    "from keras.layers import Conv2D, Dense, Dropout, Flatten, MaxPooling2D, Input, UpSampling2D\n",
    "from keras.models import Model, load_model\n",
    "from keras.utils import to_categorical\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import tensorflow as tf\n",
    "from alibi.explainers import CEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and prepare MNIST data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train shape: (60000, 28, 28) y_train shape: (60000,)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f36b0f394e0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADXVJREFUeJzt3W+oXPWdx/HPZ00bMQ2Su8FwScPeGmUlBDfViygb1krXmI2VWPxDQliyKr19UGGL+2BFhRV1QWSbpU8MpBgal27aRSOGWvpnQ1xXWEpuJKvRu60xpCQh5o9paCKBau53H9wTuSZ3ztzMnJkzc7/vF1zuzPmeM/PlJJ/7O2fOzPwcEQKQz5/U3QCAehB+ICnCDyRF+IGkCD+QFOEHkiL8QFKEH0iK8ANJzermk9nm7YRAh0WEp7NeWyO/7ZW2f2N7n+1H23ksAN3lVt/bb/sySb+VdLukQ5J2SVobEe+VbMPID3RYN0b+myTti4j9EfFHST+WtLqNxwPQRe2Ef6Gkg5PuHyqWfY7tEdujtkfbeC4AFev4C34RsUnSJonDfqCXtDPyH5a0aNL9LxfLAPSBdsK/S9K1tr9i+4uS1kjaXk1bADqt5cP+iPjU9sOSfiHpMkmbI+LdyjoD0FEtX+pr6ck45wc6ritv8gHQvwg/kBThB5Ii/EBShB9IivADSRF+ICnCDyRF+IGkCD+QFOEHkiL8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kBThB5Ii/EBShB9IquUpuiXJ9gFJpyWdk/RpRAxX0RQ+74Ybbiitb9u2rWFtaGio4m56x4oVK0rrY2NjDWsHDx6sup2+01b4C7dFxIkKHgdAF3HYDyTVbvhD0i9t77Y9UkVDALqj3cP+5RFx2PZVkn5l+/8i4o3JKxR/FPjDAPSYtkb+iDhc/D4m6RVJN02xzqaIGObFQKC3tBx+23Nszz1/W9IKSXuragxAZ7Vz2L9A0iu2zz/Ov0fEzyvpCkDHtRz+iNgv6S8q7AUN3HHHHaX12bNnd6mT3nLXXXeV1h988MGGtTVr1lTdTt/hUh+QFOEHkiL8QFKEH0iK8ANJEX4gqSo+1Yc2zZpV/s+watWqLnXSX3bv3l1af+SRRxrW5syZU7rtxx9/3FJP/YSRH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeS4jp/D7jttttK67fccktp/bnnnquynb4xb9680vqSJUsa1q644orSbbnOD2DGIvxAUoQfSIrwA0kRfiApwg8kRfiBpBwR3Xsyu3tP1kOWLl1aWn/99ddL6x999FFp/cYbb2xYO3PmTOm2/azZflu+fHnD2uDgYOm2x48fb6WlnhARns56jPxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kFTTz/Pb3izpG5KORcTSYtmApJ9IGpJ0QNL9EfH7zrXZ35544onSerPvkF+5cmVpfaZeyx8YGCit33rrraX18fHxKtuZcaYz8v9Q0oX/+x6VtCMirpW0o7gPoI80DX9EvCHp5AWLV0vaUtzeIunuivsC0GGtnvMviIgjxe0PJS2oqB8AXdL2d/hFRJS9Z9/2iKSRdp8HQLVaHfmP2h6UpOL3sUYrRsSmiBiOiOEWnwtAB7Qa/u2S1he310t6tZp2AHRL0/Db3irpfyT9ue1Dth+S9Kyk222/L+mvi/sA+kjTc/6IWNug9PWKe+lb9957b2l91apVpfV9+/aV1kdHRy+5p5ng8ccfL603u45f9nn/U6dOtdLSjMI7/ICkCD+QFOEHkiL8QFKEH0iK8ANJMUV3Be67777SerPpoJ9//vkq2+kbQ0NDpfV169aV1s+dO1daf+aZZxrWPvnkk9JtM2DkB5Ii/EBShB9IivADSRF+ICnCDyRF+IGkuM4/TVdeeWXD2s0339zWY2/cuLGt7fvVyEj5t7vNnz+/tD42NlZa37lz5yX3lAkjP5AU4QeSIvxAUoQfSIrwA0kRfiApwg8kxXX+aZo9e3bD2sKFC0u33bp1a9XtzAiLFy9ua/u9e/dW1ElOjPxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kFTT6/y2N0v6hqRjEbG0WPakpG9JOl6s9lhE/KxTTfaC06dPN6zt2bOndNvrr7++tD4wMFBaP3nyZGm9l1111VUNa82mNm/mzTffbGv77KYz8v9Q0soplv9rRCwrfmZ08IGZqGn4I+INSf079ACYUjvn/A/bftv2ZtvzKusIQFe0Gv6NkhZLWibpiKTvNVrR9ojtUdujLT4XgA5oKfwRcTQizkXEuKQfSLqpZN1NETEcEcOtNgmgei2F3/bgpLvflMTHq4A+M51LfVslfU3SfNuHJP2TpK/ZXiYpJB2Q9O0O9gigA5qGPyLWTrH4hQ700tPOnj3bsPbBBx+UbnvPPfeU1l977bXS+oYNG0rrnbR06dLS+tVXX11aHxoaaliLiFZa+sz4+Hhb22fHO/yApAg/kBThB5Ii/EBShB9IivADSbndyy2X9GR2956si6677rrS+lNPPVVav/POO0vrZV8b3mknTpworTf7/1M2zbbtlno6b+7cuaX1ssuzM1lETGvHMvIDSRF+ICnCDyRF+IGkCD+QFOEHkiL8QFJc5+8By5YtK61fc801XerkYi+99FJb22/ZsqVhbd26dW099qxZzDA/Fa7zAyhF+IGkCD+QFOEHkiL8QFKEH0iK8ANJcaG0BzSb4rtZvZft37+/Y4/d7GvF9+5lLpkyjPxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kFTT6/y2F0l6UdICSSFpU0R83/aApJ9IGpJ0QNL9EfH7zrWKflT23fztfm8/1/HbM52R/1NJ/xARSyTdLOk7tpdIelTSjoi4VtKO4j6APtE0/BFxJCLeKm6fljQmaaGk1ZLOf03LFkl3d6pJANW7pHN+20OSvirp15IWRMSRovShJk4LAPSJab+33/aXJL0s6bsR8YfJ52sREY2+n8/2iKSRdhsFUK1pjfy2v6CJ4P8oIrYVi4/aHizqg5KOTbVtRGyKiOGIGK6iYQDVaBp+TwzxL0gai4gNk0rbJa0vbq+X9Gr17QHolOkc9v+lpL+V9I7t858tfUzSs5L+w/ZDkn4n6f7OtIh+VvbV8N382nhcrGn4I+JNSY0uyH692nYAdAvv8AOSIvxAUoQfSIrwA0kRfiApwg8kxVd3o6Muv/zylrc9e/ZshZ3gQoz8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU1/nRUQ888EDD2qlTp0q3ffrpp6tuB5Mw8gNJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUlznR0ft2rWrYW3Dhg0Na5K0c+fOqtvBJIz8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5CUm82RbnuRpBclLZAUkjZFxPdtPynpW5KOF6s+FhE/a/JYTMgOdFhEeDrrTSf8g5IGI+It23Ml7ZZ0t6T7JZ2JiH+ZblOEH+i86Ya/6Tv8IuKIpCPF7dO2xyQtbK89AHW7pHN+20OSvirp18Wih22/bXuz7XkNthmxPWp7tK1OAVSq6WH/ZyvaX5L0X5L+OSK22V4g6YQmXgd4WhOnBg82eQwO+4EOq+ycX5Jsf0HSTyX9IiIu+jRGcUTw04hY2uRxCD/QYdMNf9PDftuW9IKkscnBL14IPO+bkvZeapMA6jOdV/uXS/pvSe9IGi8WPyZpraRlmjjsPyDp28WLg2WPxcgPdFilh/1VIfxA51V22A9gZiL8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIrwA0kRfiApwg8k1e0puk9I+t2k+/OLZb2oV3vr1b4kemtVlb392XRX7Orn+S96cns0IoZra6BEr/bWq31J9NaqunrjsB9IivADSdUd/k01P3+ZXu2tV/uS6K1VtfRW6zk/gPrUPfIDqEkt4be90vZvbO+z/WgdPTRi+4Dtd2zvqXuKsWIatGO2905aNmD7V7bfL35POU1aTb09aftwse/22F5VU2+LbO+0/Z7td23/fbG81n1X0lct+63rh/22L5P0W0m3SzokaZektRHxXlcbacD2AUnDEVH7NWHbfyXpjKQXz8+GZPs5SScj4tniD+e8iPjHHuntSV3izM0d6q3RzNJ/pxr3XZUzXlehjpH/Jkn7ImJ/RPxR0o8lra6hj54XEW9IOnnB4tWSthS3t2jiP0/XNeitJ0TEkYh4q7h9WtL5maVr3XclfdWijvAvlHRw0v1D6q0pv0PSL23vtj1SdzNTWDBpZqQPJS2os5kpNJ25uZsumFm6Z/ZdKzNeV40X/C62PCJukPQ3kr5THN72pJg4Z+ulyzUbJS3WxDRuRyR9r85mipmlX5b03Yj4w+Ranftuir5q2W91hP+wpEWT7n+5WNYTIuJw8fuYpFc0cZrSS46enyS1+H2s5n4+ExFHI+JcRIxL+oFq3HfFzNIvS/pRRGwrFte+76bqq679Vkf4d0m61vZXbH9R0hpJ22vo4yK25xQvxMj2HEkr1HuzD2+XtL64vV7SqzX28jm9MnNzo5mlVfO+67kZryOi6z+SVmniFf8PJD1eRw8N+rpa0v8WP+/W3ZukrZo4DPxEE6+NPCTpTyXtkPS+pP+UNNBDvf2bJmZzflsTQRusqbflmjikf1vSnuJnVd37rqSvWvYb7/ADkuIFPyApwg8kRfiBpAg/kBThB5Ii/EBShB9IivADSf0/fhI1ni26LDgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "(x_train, y_train), (x_test, y_test) = keras.datasets.mnist.load_data()\n",
    "print('x_train shape:', x_train.shape, 'y_train shape:', y_train.shape)\n",
    "plt.gray()\n",
    "plt.imshow(x_test[4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare data: scale, reshape and categorize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train shape: (60000, 28, 28, 1) x_test shape: (10000, 28, 28, 1)\n",
      "y_train shape: (60000, 10) y_test shape: (10000, 10)\n"
     ]
    }
   ],
   "source": [
    "x_train = x_train.astype('float32') / 255\n",
    "x_test = x_test.astype('float32') / 255\n",
    "x_train = np.reshape(x_train, x_train.shape + (1,))\n",
    "x_test = np.reshape(x_test, x_test.shape + (1,))\n",
    "print('x_train shape:', x_train.shape, 'x_test shape:', x_test.shape)\n",
    "y_train = to_categorical(y_train)\n",
    "y_test = to_categorical(y_test)\n",
    "print('y_train shape:', y_train.shape, 'y_test shape:', y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, xmax = -.5, .5\n",
    "x_train = ((x_train - x_train.min()) / (x_train.max() - x_train.min())) * (xmax - xmin) + xmin\n",
    "x_test = ((x_test - x_test.min()) / (x_test.max() - x_test.min())) * (xmax - xmin) + xmin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define and train CNN model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cnn_model():\n",
    "    x_in = Input(shape=(28, 28, 1))\n",
    "    x = Conv2D(filters=64, kernel_size=2, padding='same', activation='relu')(x_in)\n",
    "    x = MaxPooling2D(pool_size=2)(x)\n",
    "    x = Dropout(0.3)(x)\n",
    "    \n",
    "    x = Conv2D(filters=32, kernel_size=2, padding='same', activation='relu')(x)\n",
    "    x = MaxPooling2D(pool_size=2)(x)\n",
    "    x = Dropout(0.3)(x)\n",
    "    \n",
    "    x = Flatten()(x)\n",
    "    x = Dense(256, activation='relu')(x)\n",
    "    x = Dropout(0.5)(x)\n",
    "    x_out = Dense(10, activation='softmax')(x)\n",
    "    \n",
    "    cnn = Model(inputs=x_in, outputs=x_out)\n",
    "    cnn.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])\n",
    "    \n",
    "    return cnn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_1 (InputLayer)         (None, 28, 28, 1)         0         \n",
      "_________________________________________________________________\n",
      "conv2d_1 (Conv2D)            (None, 28, 28, 64)        320       \n",
      "_________________________________________________________________\n",
      "max_pooling2d_1 (MaxPooling2 (None, 14, 14, 64)        0         \n",
      "_________________________________________________________________\n",
      "dropout_1 (Dropout)          (None, 14, 14, 64)        0         \n",
      "_________________________________________________________________\n",
      "conv2d_2 (Conv2D)            (None, 14, 14, 32)        8224      \n",
      "_________________________________________________________________\n",
      "max_pooling2d_2 (MaxPooling2 (None, 7, 7, 32)          0         \n",
      "_________________________________________________________________\n",
      "dropout_2 (Dropout)          (None, 7, 7, 32)          0         \n",
      "_________________________________________________________________\n",
      "flatten_1 (Flatten)          (None, 1568)              0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 256)               401664    \n",
      "_________________________________________________________________\n",
      "dropout_3 (Dropout)          (None, 256)               0         \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 10)                2570      \n",
      "=================================================================\n",
      "Total params: 412,778\n",
      "Trainable params: 412,778\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "cnn = cnn_model()\n",
    "cnn.summary()\n",
    "cnn.fit(x_train, y_train, batch_size=64, epochs=5, verbose=0)\n",
    "cnn.save('mnist_cnn.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate the model on test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy:  0.9892\n"
     ]
    }
   ],
   "source": [
    "score = cnn.evaluate(x_test, y_test, verbose=0)\n",
    "print('Test accuracy: ', score[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define and train auto-encoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ae_model():\n",
    "    x_in = Input(shape=(28, 28, 1))\n",
    "    x = Conv2D(16, (3, 3), activation='relu', padding='same')(x_in)\n",
    "    x = Conv2D(16, (3, 3), activation='relu', padding='same')(x)\n",
    "    x = MaxPooling2D((2, 2), padding='same')(x)\n",
    "    encoded = Conv2D(1, (3, 3), activation=None, padding='same')(x)\n",
    "    \n",
    "    x = Conv2D(16, (3, 3), activation='relu', padding='same')(encoded)\n",
    "    x = UpSampling2D((2, 2))(x)\n",
    "    x = Conv2D(16, (3, 3), activation='relu', padding='same')(x)\n",
    "    decoded = Conv2D(1, (3, 3), activation=None, padding='same')(x)\n",
    "\n",
    "    autoencoder = Model(x_in, decoded)\n",
    "    autoencoder.compile(optimizer='adam', loss='mse')\n",
    "    \n",
    "    return autoencoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_2 (InputLayer)         (None, 28, 28, 1)         0         \n",
      "_________________________________________________________________\n",
      "conv2d_3 (Conv2D)            (None, 28, 28, 16)        160       \n",
      "_________________________________________________________________\n",
      "conv2d_4 (Conv2D)            (None, 28, 28, 16)        2320      \n",
      "_________________________________________________________________\n",
      "max_pooling2d_3 (MaxPooling2 (None, 14, 14, 16)        0         \n",
      "_________________________________________________________________\n",
      "conv2d_5 (Conv2D)            (None, 14, 14, 1)         145       \n",
      "_________________________________________________________________\n",
      "conv2d_6 (Conv2D)            (None, 14, 14, 16)        160       \n",
      "_________________________________________________________________\n",
      "up_sampling2d_1 (UpSampling2 (None, 28, 28, 16)        0         \n",
      "_________________________________________________________________\n",
      "conv2d_7 (Conv2D)            (None, 28, 28, 16)        2320      \n",
      "_________________________________________________________________\n",
      "conv2d_8 (Conv2D)            (None, 28, 28, 1)         145       \n",
      "=================================================================\n",
      "Total params: 5,250\n",
      "Trainable params: 5,250\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "ae = ae_model()\n",
    "ae.summary()\n",
    "ae.fit(x_train, x_train, batch_size=128, epochs=10, validation_data=(x_test, x_test), verbose=0)\n",
    "ae.save('mnist_ae.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare original with decoded images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAwAAADnCAYAAACaL8QkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XuwXVV9OPB9Awkh4RkIL0HCY0AjLyGoYIBCGUCa8BAQBK1KxSoyUqHig7YCtjClCkPFEHCGUWAmZArKQHmoLSAwVV7yKPIooATBIISQQCDhkdzfH7/p6lpr56y7z+Hcc09uPp+/vnu+55y9bjjr7HM267u+A4ODgxUAAABAbMxIDwAAAADoP24YAAAAADVuGAAAAAA1bhgAAAAANW4YAAAAADVrlpIDAwNaKPSPBYODg5NHehD0B3Ozr5ibBOZmXzE3CczNvmJuEpibfWWlc9MKg1XHvJEeALBS5ib0J3MT+pO5Cf1ppXPTDQMAAACgxg0DAAAAoMYNAwAAAKDGDQMAAACgxg0DAAAAoMYNAwAAAKDGDQMAAACgxg0DAAAAoMYNAwAAAKBmzZEeQFN/+7d/mxyvvfbaId5ll12S3NFHH93ydS655JLk+Fe/+lWIr7zyynczRAAAABg1rDAAAAAAatwwAAAAAGr6uiRh7ty5IS6VGeRWrFjRMvfXf/3XyfGBBx4Y4l/+8pdJ7tlnn218TqC7dthhhxA//vjjSe7UU08N8fe///2ejQlGi4kTJybH//Iv/xLi/Dp5//33J8fHHHNMiOfNmzcMowMA+oUVBgAAAECNGwYAAABAjRsGAAAAQE1f7WEQ71lQVc33Lcjrm3/2s5+FeNttt01yM2fOTI632267EJ9wwglJ7rzzzmt0fqD7PvjBD4Y435fkueee6/VwYFTZfPPNk+OTTjopxPl822OPPZLjGTNmhPgHP/jBMIwORrfdd989Of7JT34S4ilTpgz7+Q866KDk+LHHHgvxH/7wh2E/P6xu4t+f119/fZI75ZRTQjx79uwkt3z58uEdWENWGAAAAAA1bhgAAAAANSNekjBt2rQQH3nkkS0f99vf/jY5Puyww0K8YMGCJLdkyZIQjxs3Lsn9+te/To533XXXEG+00UYNRgz0wm677Rbi119/Pcn99Kc/7fVwYJU3efLkEP/4xz8ewZHA6u3ggw9Ojtdaa62enj8vzz3xxBNDfNxxx/V0LDAa5b8pZ82a1fKxF198cYgvv/zyJLd06dLuDqxDVhgAAAAANW4YAAAAADVuGAAAAAA1I76HQdzaaWBgIMnF+xbk9V7z589v9Pqnn356cjx16tSWj73xxhsbvSbQfTvttFNyHLeZufLKK3s9HFjlfeUrX0mOjzjiiBB/6EMf6vh199133xCPGZP+f4eHHnooxHfccUfH54DRZs01/+8r96GHHjqCI6mq+++/Pzk+7bTTQjxx4sQkl+8hBAwtvk5WVVVtueWWLR87Z86cEC9btmzYxvRuWGEAAAAA1LhhAAAAANSMeEnCDTfcEOLtt98+yb322mshXrhwYUevn7eHGTt2bEevAwyv973vfclxvCxy7ty5vR4OrPIuvPDC5HjFihVded2Pf/zjK42rqqrmzZsX4mOPPTbJ5cugYXWy//77h3ivvfZKcueff35Px7Lhhhsmx3G57oQJE5KckgQYWt4a9cwzz2z83LjsdnBwsGtj6iYrDAAAAIAaNwwAAACAGjcMAAAAgJoR38MgFtc+vhtf+9rXQrzDDjsUH3v33XevNAZ664wzzkiO48+D++67r9fDgVXSTTfdFOK85WGnXn755eR4yZIlId56662T3DbbbBPie+65J8mtscYaXRkPrAryVsFx67Snn346yZ177rk9GdP/Ovzww3t6Phjtdt555+R4jz32aPnYd955Jzm++eabh2VM3WSFAQAAAFDjhgEAAABQ01clCZ2aMWNGcnzOOeeEeNy4cUnuxRdfTI6/+c1vhviNN94YhtEBKzNlypTkeNq0acnx//zP/4RYWydYuf322y853nHHHUOct1Fs2lZx9uzZyfHPf/7z5Hjx4sUhPuCAA5JcqZXUl770pRBfcskljcYCq6q/+7u/S47jVsGHHHJIkovLfIbLpEmTQpx/bnSr5Sqsro466qjGj82vqasCKwwAAACAGjcMAAAAgBo3DAAAAICaUbGHQV77nO9bEJs7d25y/Mtf/nJYxgSU5TWUuZdeeqlHI4FVS7z/x9VXX53kNt5440avkbcxvvbaa0N89tlnJ7nS/j7563zhC18I8eTJk5Pc+eefH+Lx48cnuYsvvjjEb7/9dsvzQT87+uijQ3zooYcmuaeeeirEI9EqON5fJN+z4Pbbbw/xokWLejUkGDX23XffYv6tt94KcWmvn35lhQEAAABQ44YBAAAAULPKliRcd911IT7ooINaPu6KK65IjvM2N8DI2HnnnYv5ePky8H/WXPP/Lt1NSxCqKi3BO+6445LcggULOhpLXpJw3nnnhfiCCy5IchMmTAhxPr+vv/76ED/99NMdjQVG2jHHHBPi+P1eVVU1a9asno4lb118wgknhHj58uVJ7h//8R9DrCQImtl7771XGq9M3B78wQcfHLYxDRcrDAAAAIAaNwwAAACAGjcMAAAAgJpVZg+DzTffPDmOa0XWWmutJBfXYsZ1WVVVVUuWLBmG0QFNfOQjHwnx5z73uST3wAMPJMe/+MUvejImGK3y1m0nnnhiiDvds2Ao8V4Ecc10VVXVnnvuOSznhJGy/vrrJ8fxNS53ySWXDPdwEnGL06pK9zt57LHHktxtt93WkzHBaNLONa3X87/brDAAAAAAatwwAAAAAGpWmZKEa6+9NjneaKONWj72qquuCrH2TNA/DjzwwBBPmjQpyd1yyy3J8bJly3oyJliVjRnT+r7/hz/84R6O5P8bGBgIcT620ljPOuusEH/605/u+rhgOOQlse95z3tCPGfOnF4PJ7Hddtu1zD3yyCM9HAmMTtOmTWuZW7RoUXKsJAEAAAAYddwwAAAAAGrcMAAAAABq+noPg8MOOyzEu+++e8vH3X777cnxt7/97eEaEvAu7LrrriEeHBxMctdcc02vhwOrpC9+8YshXrFixQiOpG7mzJkh/uAHP5jk4rHm4473MIBVxWuvvZYcP/jggyHeZZddkly8b8/ChQuHZTybbLJJiI8++uiWj7vrrruG5fwwmk2fPj05Pv7441s+dvHixcnxc889Nyxj6hUrDAAAAIAaNwwAAACAGjcMAAAAgJq+2sNgo402So6/9a1vhXjs2LEtnxfXjFVVVS1ZsqS7AwM6stlmmyXH++yzT4ifeOKJJPfTn/60J2OCVV28T8BImDx5coinTp2a5OLrdslLL72UHL/99tvvfmDQY0uXLk2On3766RAfddRRSe7GG28M8QUXXNDR+XbaaafkeNttt02Op0yZEuJ8n6BYv+19AquC/HfqmDGt/7/7L37xi+EeTk9ZYQAAAADUuGEAAAAA1PRVScLpp5+eHO+5554tH3vdddeFWBtF6E+f/exnk+O45dPNN9/c49EA3XDmmWeG+Mtf/nLj5z3zzDMh/sxnPpPknn322Xc9Lhhp8ffRgYGBJPcXf/EXIZ4zZ05Hr79gwYLkOC872HjjjRu9zo9+9KOOzg+rs1Kr0kWLFiXHl1566XAPp6esMAAAAABq3DAAAAAAatwwAAAAAGr6ag+D0047rfFjTznllBBrowj9aeutt26Ze+WVV3o4EqBTN910U3K84447dvQ6jz76aIjvuuuudzUm6EePP/54iD/xiU8kud122y3E22+/fUevf8011xTzP/7xj0N8wgkntHxc3g4SWLktt9wyxMcff3zLxz333HPJ8X333TdsYxoJVhgAAAAANW4YAAAAADV9VZLQjkmTJoX47bff7vh1Fi9e3PJ1xo4dG+L111+/5WtssMEGyXHT0orly5cnx1//+tdD/MYbbzR6DehnM2bMaJm74YYbejgSGD3idm1jxrS+7/+xj32sZe6yyy5LjrfYYouWj83PsWLFiqGGuFIzZ87s6HkwGjz44IMrjbvpd7/7XaPH7bTTTsnxI488MhzDgVXe3nvvHeLS9fa6667rxXBGjBUGAAAAQI0bBgAAAECNGwYAAABAzSq7h8HDDz/cldf5t3/7txDPnz8/yW266aYhPvbYY7tyvpIXXnghxP/0T/807OeD4TB9+vQQb7bZZiM4EhidLrnkkhCff/75LR/37//+78lxae+BdvYlaPrY2bNnN35N4N2L9zeJ45w9C6CZjTbaqGVuwYIFIb7ooot6MZwRY4UBAAAAUOOGAQAAAFDTVyUJN910U3J8+OGHD/s5jznmmI6e984774S4tDzz+uuvT47vu+++lo+98847OxoL9JMjjzwyxGussUaSe+CBB0J8xx139GxMMJr85Cc/CfHXvva1JDd58uRhP/9LL70U4sceeyzJfeELXwhxXuYHDK/BwcGVxkBnDj744Ja5Z599NsSLFy/uxXBGjBUGAAAAQI0bBgAAAECNGwYAAABATV/tYfDxj388OT7jjDNCPHbs2Mav84EPfCDE7bRDvPzyy5PjZ555puVjr7322hA//vjjjc8Bo82ECROS40MPPbTlY6+55poQL1++fNjGBKPZvHnzQnzccccluSOOOCLEp5566rCcP277+4Mf/GBYzgG0b/z48S1zS5cu7eFIYNWU/97cbrvtWj522bJlIX777beHbUz9wAoDAAAAoMYNAwAAAKCmr0oScueff/67fo3jjz++CyMBWsmXYb3yyishztuKXnTRRT0ZE6wu8vak8fHPf/7zJBe3PJw5c2aSi+fqZZddluQGBgaS40cffbSzwQLD6nOf+1yIFy1alOS+853v9Ho4sMpZsWJFcnzfffeFeKeddkpyTz31VE/G1A+sMAAAAABq3DAAAAAAatwwAAAAAGr6eg8DoP/lexjsvffeIzQSIHbLLbcUj4HR5d577w3xBRdckORuu+22Xg8HVjl5y+8zzzwzxIODg0nu/vvv78mY+oEVBgAAAECNGwYAAABAzUC+vCJJDgy0TtJr9w8ODk4b6UHQH8zNvmJuEpibfcXcJDA3+4q5SWBu9pWVzk0rDAAAAIAaNwwAAACAGjcMAAAAgBo3DAAAAIAaNwwAAACAGjcMAAAAgBo3DAAAAIAaNwwAAACAGjcMAAAAgBo3DAAAAICaNYfIL6iqal4vBsKQth7pAdBXzM3+YW4SMzf7h7lJzNzsH+YmMXOzf6x0bg4MDg72eiAAAABAn1OSAAAAANS4YQAAAADUuGEAAAAA1LhhAAAAANS4YQAAAADUuGEAAAAA1LhhAAAAANS4YQAAAADUuGEAAAAA1LhhAAAAANS4YQAAAADUuGEAAAAA1LhhAAAAANS4YQAAAADUuGEAAAAA1LhhAAAAANS4YQAAAADUuGEAAAAA1KxZSg4MDAyOGbP63FMYHBxMjgcGBkbs/Pm/+/LlyxcMDg5O7umA6Fv53Mzfu9ljGz1uOORzqOk4c+08r9O/sem/U36+FStWmJsE7cxNusvcpGR1+0470nynpSnXzZHT9LpZvGEwZsyYap111gnHy5cvb3mCWOk/dKc/yvPHNX1efr4VK1a0fGyeW3PN//vn6caPkPw4P198vNZaayW5RYsWzetoAIxKY8aMqdZee+1wXHpfD/cNg9Kcfjc3DOLj+LNnqHOU/i2anq80N8ePH5/kXn31VXOTYMyYMcl7JH/vxoZ7bpbe4+/mOl2aK8Mh/iLZznVzyZIl5ibBmDFjqokTJ4bjblw3S9e4Tq933fos6Nb/BItfp1vXTd9piZWum936n0lNn9eL12nnu3Gnnwfdvm661QoAAADUuGEAAAAA1BRLEgYHB6u33norHJdqv0pLNpouH25niUb8vHaWRMdlBlWVLnvJX6dUf1VadtZ0+Ur+uPgcS5cubfkaMDg4WL399tvhOH7vdFrT386yq6ZLq0rPy5dI5Y+N/6Z8/pWWWpX+LUpLUEvzPX4dc5OSFStWtLxutlM/3WnZXWlpc+n60+k5Sn9Tt5ZWl74nuG7S1HBcN/NrStO50Y5OSyJK5QOdXrdz8d/vukmnStfNpr8v29GN8px2X6f0vbU0x0vPK42l2783rTAAAAAAatwwAAAAAGqGaqvYlU4BTZdarOz8/yvfYTVeVlY63xprrFF8XmmnyPix7SyJKS0DafpvMW7cuCT35ptvtnweq5/S3Ox0yVQ787vpY7u1E3uu6bLrUilD6XmlseVzs/RZxOqnnblZMhxzrGm5Tn7dzMWPbackqBs7QbczNy2DJjYc1812xOcolTLk869UgtupXnQ3Kc1N32mJDQwMVGPHjg3Hnb7nO+1u0qnS7838/O+8807L8zct7Xs3pYStztf0ummFAQAAAFDjhgEAAABQ44YBAAAAUFPcwyDXjdqQUi6u76iqqpo4cWKIp0+fnuT+/M//PMTvf//7k9yWW24Z4rgmpqqq6uWXX06OzznnnBDffvvtSS6uRynVgnZLt2rTGP0GBwc7rmsajrHEOm2jmu8hEtc4xp8FVZW2Q81rIfPWqa3O3w5zk3a0mpvd2geo1Lqt6bU5f1w8b4ba+6fpPj1N27EN9dimzFOG0sl1s53517Tldz7H4vnXTuu0Ttsjluqt8+/ine6v0I3PPlYf3Zib3fhe3M6+e/H7P/5eurLzx8/N51H83NK+QJ1+3pT+nZqywgAAAACoccMAAAAAqGmrJKEbS4pKrzFhwoTk+MILLwzxEUcc0fJ18mUf66yzTojz5cpTpkxJji+++OIQz5gxI8k9+eSTIR6qzVRTw9HiDkby/dK0HUxVpfMonwtvvfVWcrz11luH+I477khy5513XohnzZqV5EplD6WWi91oYwe54W7tVGodmi8tjuUlQLNnzw7xYYcdluT++7//Ozn+4he/GOL4OllVaflCpy0Xh2MpKeSavl9K78dOl/3ncyN+bL60OX7sUOWwpTnWapz5cena2E77yU5bykJTw1GO2+l8L7VRrKp0XuVlfvGcz8tqu/Edop2/qRUrDAAAAIAaNwwAAACAGjcMAAAAgJoh9zBoVdvQrRrDDTbYIMTf+973ktzhhx8e4rzeY8mSJSF++OGHk9xvfvObEG+88cZJLq/NnDx5coiPPPLIJPf9738/xHktSlxvMlSbGxjN2qmvKtV75XsYbL/99iGeNGlSkovbo+b1lnEtaKmGTE0lI61x7WD2Ho+vP3m9c3ycv358Hd10002TXLyHT369nTp1anL8gQ98IMRPP/10khs3blyIly1bVrVSqpMutVx0faXX2mkdXHofx22+d9xxxyT3z//8zyE+4YQTktwrr7zSeGyx0n4Hm222WXIcf27Mnz+/5eu00w5VW0Xa0e3fm71W2rOrqqrq+OOPD/HZZ5+d5ObMmRPiM888s/E5eskKAwAAAKDGDQMAAACgZsiShG60ASzl4tZpeevE119/PcSPPPJIkjv11FNDPG/evCQXL6fcaaedktwhhxySHMfLsCZOnJjk4uWUeZuLdlrJxVaVpTXQDaW5kC9tjJcyV1X6efDiiy8muZ/97GchLrU8Nd/oJ6XlyyWlZcDttHnacMMNQ3zRRRclufj6l5cSxNfiqqqqLbbYIsR5KdH48eNDnF8387+j1biHah0H/app67T99tsvye28884hzq+FsaE+M0rtSeNWqgcccECS23///UP8+c9/PsnFc7ydlm/KEGhHv3xfa+c3XOlatf766yfH3/jGN0IcX4urKi2XP+uss5JcXFrbTqvUbs8/V2UAAACgxg0DAAAAoMYNAwAAAKCma20V4+O8TjGu8YhrqKoqrYVcvHhxknviiSdCfPLJJye5J598MsRxq5p8LOuuu26Si9s45o996qmnklzT+uu8hkUrGVZnTedNXguZt3mLW9DMnj07yS1atCjEpfnX6V4jMBzaqY1s2i6ttIdB3kb461//eojf+973JrmlS5eGON8XJN6XoKqqap999gnxwoULk9x//ud/hjhvzxaPLW95Gn8elP4tzFt6rZ0WoPFj82tcLJ+b8XfafE7F5xiqzrv0XfzVV19t+by99torxPl36ngPg7yNa+ka2y816awaSvtvDOe58vO10yq11WtUVX2fkK222irE+Ty65557Wp4j3tMkz/Xy2miFAQAAAFDjhgEAAABQM2RJQqy0RKTpUpK8zcx//Md/hDhvgRi3diotX8yXYcRLNr70pS8luXxpZVwGcddddyW50jKQpstXYHUTL1HM50a8DCvPHX300clxXL70wx/+sNH58uPS8jHotabX0KEeG8vf47vsskuIZ82aleTi1olvvvlmkitdb/M59tGPfjTEf/Znf5bk4haoJ510UpK77bbbQpwv146/G5RapUIvlJb2x7lS69B8Dh944IEhzr/vnnvuuSHO25qWynVKS5RLtt122+Q4bwEXi+e/Mj9Gg3ZKHkq/9+L5l1+38utf/JsyLgGsqqq69NJLW54/vh7nc6yXLYitMAAAAABq3DAAAAAAatwwAAAAAGra2sOgU3HNRd5KIq7/yOsmS+3R4lqRvKbkxBNPDPEhhxyS5PK6zccffzzECxYsWPkfsBJN60a0lWG0K+3hkefidk15vdenPvWp5Diej3FddFWV51/+GVMaKwy3pm1223lvxu/xuOVZVVXVZz7zmRDHexbk54/rKauqXIuZz6k33nij5TmmTJkS4r//+79Pcg888EDLccd7GLRTp6lumuFQuo6V5nSce8973pPkvve974X4scceS3JXXXVViEt7eJTGkh+XWg5/6EMfSnLx37H22msnuSVLloS4tGdDztykX3XaVrE033bdddckN23atOQ4vo7mvzcfeuihlmON51hpTmurCAAAAPScGwYAAABATVslCU2XZZTaTuRLJuLj0jLEvFwhXk658847J7lvfOMbIc5LEPKxxa1s4mVXVZUun86VWufAaNNpO7j8efE8njRpUpLbZpttkuPnn38+xK+//nqSa2fJZmmsMNyavufauY7Er5mX8hx//PEhjksHqiq9ppXm7Zw5c5LcddddlxzHSyvj81VVVR1xxBEhjtsvVlVVffKTnwzx7Nmzk1zcRrX0PWGoJdnQS/mS/HiOnX322Ulus802C/GXv/zlJLdw4cIQ58uOmy5Jzh+bf29+73vfG+LddtstycUlQfnnRtPv8OYio0Hp2pjn4rnxiU98IsnlvyHjx956660tz5GLv++W2qh22pq5KSsMAAAAgBo3DAAAAIAaNwwAAACAmrb2MGhaD1Hai6BU/5QrtU6M/c3f/E1yvOGGG4Y4bwcV71lQVVV15513rvR8VVVuF6NWi9VJp+3g8jm0dOnSEO+zzz5JLq5hrqqqevLJJ0Oc12KW2qyVPm+gX5Wujfn7fZ111gnx6aefnuRKe4jEc+yPf/xjkotbvl1xxRVJLr+OxuN5+OGHk1y8p1Bpf6FFixYluRtuuCHE+XwvXYtj5jvdUnovld6PcbvCmTNnJrm77747xPfcc0+Sa/qdcqi5EH+Pzffw+vSnPx3i/Hobjyd/Xmnvk06/30M/Ke3FEV//8mtxPI+mT5+e5OJ9QfJzxG1Uqyqdc6W2xqXvvsP9u9QKAwAAAKDGDQMAAACgpq2ShE41bfvQzjKQM888M8QHH3xwkovbI958881JLm/l1HQJR/64XraygH5Wev/nyycnTJgQ4rjFWlXVl0FeeeWVLc9ROn/T5cvQC6WWTJ0uQ47n0SabbJLk4mtlvrQxXvb/la98JcnFbd1ypTZvr732WpKbNWtWiM8777wkF7eV+/a3v53kbrnllhDn/07xdwHXV3qhaWlbPjfi61r+Xr3gggtCXGoVXFp2XGqrlufzVsWf/exnQ5y3Tvzud78b4rwkaM01/++nQn5+85HRoPQ+jnP53JgyZUqIt9tuu+JrPvfccyH+zW9+k+Sazv/Sb9HhZoUBAAAAUOOGAQAAAFDjhgEAAABQ0/EeBu3UYpbaVZReMz7eZZddktzJJ58c4rfeeivJvfrqqyH+zne+k+TyxzZtnWifAlY3pfd/aU7Hcyp/3lZbbRXij370o0kuru+qqqq69dZbQzx+/PhG48xp68RI6Mb1obRPTlxTnItrIefMmZPkzjjjjBDnew+UXjNvDzVu3LgQ5/WW119/fYhPO+20JDdx4sQQb7HFFkkuvsbH7Y7z89mjhF5oeo3Lr01xu+C8dei9994b4rytYfyaeRvTVuNamXgeH3nkkUluvfXWC/HTTz/dcmz5Z0Hp38L1l37Vzp55pefF17h8D4ODDjooxPlnQX6tuuaaa0Ic/06tqnRfotI1biTnlBUGAAAAQI0bBgAAAEBNxyUJ3VqGVHqdsWPHhjhf2hgv2cpbR11xxRUh/tOf/tTyeflY21lG2rTNm9IFRoPS0q58SXI8x/KlVQcccECI119//SQ3d+7c5Dhu+1Q6v3Ih+k2r60PpvVlql5bn4nZpG264YZJbtmxZiM8555wkFy+njJf552PL53R+jS1d4+LXjZdZ5vKlnVOnTg3xr3/965Zjy5U+CyDX6r1Ueu/kz4nnx9KlS5NcXGrw1FNPJbmXX345xHF5TlWl86G07D/P5XMztueeeybH8Xfq3//+90muVAZRapUcn9/8o590+l0wf17p9+YhhxzSMpeXwMe/TfNrY9P5pyQBAAAA6CtuGAAAAAA1bhgAAAAANV1rqxgr1V+V2rPlNRwzZswI8YEHHpjk4rYv8+fPT3I//OEPQ5zXkOS1mflx03GXdFpfHdeGlVpcQS+U6pTjeZO/x+M6rnx+zZw5M8T53Lz22muT4zift6CKPyvsWUA/63S/n/h5+fVg//33b/m8N998M8R566a4TjrfQyA+zuur83HG+xTkNdx77713iCdPntzydfLPhgULFoS4ndZt8etouchQWl3XSq0T8++m8V4AeW7hwoUhXnfddZPcBhtsEOK8rWl8vtKeIUN9h912221DvPXWW1et5G2M8+txq7GVvjNDP2mnrWLp+248H/M5tfvuu4c4/yzI59gLL7wQ4rz0qQFlAAAMQklEQVRVcdPz53q5Z56ZDwAAANS4YQAAAADUuGEAAAAA1HStUL5pbUguzq233npJ7qtf/WqI8xrOuIbyX//1X5Pcn/70p5U+bmVKtdjxcbd6XzatN9HPlpHW9P2ZPy6uhd5qq62S3B577BHi559/Psndc889yXFcJ5prWtMFI6FpLWHpcaXez3F/9bxuMp5/W265ZZKLayrzc8fX2Hx/g7xuOa6/POigg5Lcd7/73ZZji8/58MMPJ7lbb7210flKtebmPsOhNE/zXDzHPvaxjyW5yy67LMR33HFHkvvDH/4Q4vx9HM+HDTfcMMlNnTo1OY6vsXm9dTyP/uu//ivJlfZQiD8b2vl+D/2q6e/Sqkr395g+fXqSy3+3xu68886WudJeIO38Fu3lHl5WGAAAAAA1bhgAAAAANT3p3VdqCRMvpzj55JOT3Pbbbx/ifInU1VdfHeIrr7wyycWPzZdWttMCrvS84W7HCCOtVHbQtCQhb4c6YcKEEMdLkKuqqt54443kOF4GWVqilTOPWFWU5lGpzG+LLbYIcd5yNF4+mbdKi0sESm0Vc/G8raqq+tSnPhXis846K8mts846Ic5LCeMWjN/61reSXNwCsvQ9QVs3uqXp/MvbdcZL+/OymzPOOCPE+dzca6+9QpyX8sTyuRjP43ws+fH48eNXGueefPLJ5LhUvltqV9rp910Ybu18by21Do3LY/PvtPH8z1slxr9TqypteZxfG5u2Ch/J77euvAAAAECNGwYAAABAjRsGAAAAQE3X9jDotO1DXOP1l3/5l0kurr/6/e9/n+TOO++8EC9ZsiTJldpD5eK6lbxOq2kdSb6/QnzO/DXjx+a1b3GdzFDjhn6Rv4/j93he77X22muH+KGHHkpyeb11PB/amZvqKOlX3drTJr425tefeC+A+fPnJ7m4NjOvt4zPt9lmmyW5U045JTmOWx7nddrLli0LcT5v4/2G7r777pbnz//2eNyl2lMYStPW1qX3XGlPjVJbxW233TbEO+ywQ5LbdNNNW47l2WefDfEzzzyT5PL2w/GeIkcffXSSW7hwYYjnzZuX5ErfTV1vGQ1K15jSe3zSpEkh3m+//ZJcPP/zufnggw+2fGyu9HnTL3PMCgMAAACgxg0DAAAAoKbjkoR2lkyU2irGuY022ijJxUuk8vNNnjw5xIsWLUpycQuovFXMCy+8kBzHyyfzpZ2lpY5x66idd945ye22224hnjZtWpKLx/PYY48luUsvvTTE+dIW6IZOlzrljyu1oImXS+dzI5bPxbzNTOn88Tn7ZbkW/K9Wy55L8y+/3sRLjeN2TFWVlujtuOOOSS5uuXjMMcckuddffz3EJ510UpLbd999Q5yXK7TTgjH+O/Jr84UXXhjifCl1/LzSnG6nPAlyTd8vTVsJ5te/eO7k5TpxiVBculBVaRlA/l00ls+b/O+58cYbQ3zUUUcludJ1Mz5/0+/z0M9K5fGl93F+ffvIRz4S4vi3X27u3LnJcVweWFVpSe6qOMesMAAAAABq3DAAAAAAatwwAAAAAGq61laxpFSPEddULV68OMltsMEGIY7rMquqqq666qoQx3sWVFVaC52fOz/Hb3/72xC/9NJLSW7zzTcP8TbbbJPk4v0W8pqy+G/Ka9ji2pgDDjggycW1b9dee22Sy1vQQSc6betWqufMX3Pq1KkhzltFxS1QH3300SSX7zdSqqnUSo3RoJ1a/di9994b4g9/+MNJLq6xPPfcc5Nc6dr0xhtvhDivoc7rpOO9f/J5Gz/2ggsuSHILFiwIcd66rbTXEQy3duqdm14P23mPl85XarkW70tSVenczHPxuEtt5fq1hhq6pbSHRz5v43ao+fOWLl0a4vh3aVVV1cSJE1ueo/R5069clQEAAIAaNwwAAACAmmEpSSgtrciXYcXtou65554kN2PGjJavE5cd5OeLl13lS0vWX3/95Dhul5EvH4nLB/I2U3EZwltvvZXknn322RA/8MADSe7uu+8Ocd5W58EHHwxxvMwFRkJp+VS8ZDkvyfnkJz8Z4rxc6IknnlhpPJTS8smc5ZT0k1YtFqsqvR7m16r4OJ9jl19+eYgPPvjgJPf+978/xOuuu27LseTnKy2zzssX4uvhnXfemeTOOuusEOdzPL7e598F4s+UdpZrmu90QztLgkvL90stEWPtlPnF30Xz5+WPjVsZ558bcZvTvFVr6XPKHGNV1M77OP5Nmc/h7bbbLsT5tTC+xr344otJLr/GlUqCVoVW4VYYAAAAADVuGAAAAAA1bhgAAAAANV3bwyCuxyjVjeT1V3Ed5Ve/+tUkd9ddd4U4bhVVVWn9Vb4vQdzWbY899khy6623XnL8/PPPh/h3v/tdknvyySdDnLc1jI/nzZuX5OL6zlItSl5fFrfDyetkoBvaqQUutXKKH5u/j4888siWubjtTNzGrarq7dlKrdV6XeOl5RTDrdRmLa+pnD9/foj/6q/+KsmdfvrpIT700EOTXLynSH4tjudffr58rs6aNSvEl112WZKL2xPnrxOfM/97S/MqzrXT/g6aaue90/S92o6m36FL1+mqqqpNNtmkZS5us5i/TrwXV34tbjXOnLnJaBDvZ1BVVfW+970vxPneH4sXLw5xaV+CXD4fRnJ+NN2/xQoDAAAAoMYNAwAAAKCmayUJpeUU8XKHuD1MVaUlCQsXLkxy8bLHfMlEvNQxX/ZRWkqdL7WKl2Hl7RGbtsMqLTspLS3L/y3aaesDnWg6T4d6Xjz/8pajjzzySIjvu+++JDd79uyVvkZV1ZdIl9rM9Hqpo+WUdKobLQHzuRG/Zl4S9/nPfz7E06dPT3Jx2d9WW22V5P74xz+G+Kabbkpy+fELL7wQ4rw9Y1yGlF8b41w3lm6/m9eBflJ6H5euhXn56jbbbBPiJUuWJLkf/ehHIX7ttdeSXFyu1M71VjtGeq3p99jStbf0fTO/psWl6/vvv3+Si7//5qXzr776asvzl3435pr+ve20I276+jErDAAAAIAaNwwAAACAGjcMAAAAgJoh9zBoVQPRaW1SXrdcqs0q1ULGNSb5GOPalPw1830DOm3z1GosQ50/zuX/FvHr5OOE4VaqfyrVQi1dujQ5PvbYY1s+L35fD1XvmNeYjSS1mHSqVOMbK82H/BoTPzbfQyS+jv3qV79KcnGr4tK1uFTfmZ8/z8VzPD9Hq9fIj/tp7kMvlD4nSt9T87rphx9+OMRXX311kps7d26I4z0L8nO009ax6f4GMJRu7+PW6Z5d+TX1m9/8Zojza1q8v0/+vLw9Y2lsnc6VXj7PCgMAAACgxg0DAAAAoKbjtoqdLh0ptYcqvWa+JLO0XDPWTuuKfGydLmcptXUsvUa8fEWLRUZa6f1fKq0pLSculSCV5ngv5kM32tNAt3Tayqy0DDIu82s6T1d2jvi5Q5UvtNLOkshutKaEftZpO7jFixcnx//wD/8Q4jfffDPJxS0Y88+JUiltrJ3SXddN+knT77T5e/zll18O8Yknnpjk4hK8+PpaVeXvxu2U/TQ13PPNCgMAAACgxg0DAAAAoMYNAwAAAKCm4z0Mct2oN26nbjiuDclrn9vZtyCuOSm9TmmfhFyp9rST14B2DNWusBtK9Y5N58pQ87TXewqUPsPiXDufLzAc79WmNf2lVqXttCMunT+v0yy1RyzNndI5m36/MDfpVC+umyVNrz/ttByN9yyoqvIeQp22Tmz1uPyx5iadKr2v2rm+Nt0bpLRPR/7+L+0Z1M7vxqaG4/tE07lpBgMAAAA1bhgAAAAANUOWJMTLH+L2EaUlG+0smWi67Ku07LmdXD7u0mO7sSStnX+LdpaasXobGBhI3svx3Ox0WWA7bV7iZcjttHxr9bihxtPpMshu6UbLG1YfreZmPv9KyyCbKl3j2pmbpbZuedlB/Drx35efv9OlpO3MMXOTdnT7utktpetdaU7n3xvjJdL5kuj4se0spS61/G5aSmRuMpROfm+WXqOd77Sl8vTYuHHjkuO8dWmsndalI1lm23RuWmEAAAAA1LhhAAAAANS4YQAAAADUDAxRC/JSVVXzejccCrYeHBycPNKDoD+Ym33F3CQwN/uKuUlgbvYVc5PA3OwrK52bxRsGAAAAwOpJSQIAAABQ44YBAAAAUOOGAQAAAFDjhgEAAABQ44YBAAAAUPP/ADD8REr0AL6vAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x288 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "decoded_imgs = ae.predict(x_test)\n",
    "n = 5\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(1, n+1):\n",
    "    # display original\n",
    "    ax = plt.subplot(2, n, i)\n",
    "    plt.imshow(x_test[i].reshape(28, 28))\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "    # display reconstruction\n",
    "    ax = plt.subplot(2, n, i + n)\n",
    "    plt.imshow(decoded_imgs[i].reshape(28, 28))\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate contrastive explanation with pertinent negative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explained instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = 4\n",
    "X = x_test[idx].reshape((1,) + x_test[idx].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CEM parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 'PN'  # 'PN' (pertinent negative) or 'PP' (pertinent positive)\n",
    "shape = (1,) + x_train.shape[1:]  # instance shape\n",
    "kappa = 0.  # minimum difference needed between the prediction probability for the perturbed instance on the\n",
    "            # class predicted by the original instance and the max probability on the other classes \n",
    "            # in order for the first loss term to be minimized\n",
    "beta = .1  # weight of the L1 loss term\n",
    "gamma = 100  # weight of the optional auto-encoder loss term\n",
    "c_init = 1.  # initial weight c of the loss term encouraging to predict a different class (PN) or \n",
    "              # the same class (PP) for the perturbed instance compared to the original instance to be explained\n",
    "c_steps = 10  # nb of updates for c\n",
    "max_iterations = 1000  # nb of iterations per value of c\n",
    "feature_range = (x_train.min(),x_train.max())  # feature range for the perturbed instance\n",
    "clip = (-1000.,1000.)  # gradient clipping\n",
    "lr = 1e-2  # initial learning rate\n",
    "no_info_val = -1. # a value, float or feature-wise, which can be seen as containing no info to make a prediction\n",
    "                  # perturbations towards this value means removing features, and away means adding features\n",
    "                  # for our MNIST images, the background (-0.5) is the least informative, \n",
    "                  # so positive/negative perturbations imply adding/removing features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate pertinent negative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# initialize TensorFlow session before model definition\n",
    "sess = tf.Session()\n",
    "K.set_session(sess)\n",
    "sess.run(tf.global_variables_initializer())\n",
    "\n",
    "# define models\n",
    "cnn = load_model('mnist_cnn.h5')\n",
    "ae = load_model('mnist_ae.h5')\n",
    "\n",
    "# initialize CEM explainer and explain instance\n",
    "cem = CEM(sess, cnn, mode, shape, kappa=kappa, beta=beta, feature_range=feature_range, \n",
    "          gamma=gamma, ae_model=ae, max_iterations=max_iterations, \n",
    "          c_init=c_init, c_steps=c_steps, learning_rate_init=lr, clip=clip, no_info_val=no_info_val)\n",
    "explanation = cem.explain(X, verbose=False)\n",
    "\n",
    "sess.close()\n",
    "K.clear_session()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Original instance and prediction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original instance prediction: 4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f36453d0828>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADXVJREFUeJzt3W+oXPWdx/HPZ00bMQ2Su8FwScPeGmUlBDfViygb1krXmI2VWPxDQliyKr19UGGL+2BFhRV1QWSbpU8MpBgal27aRSOGWvpnQ1xXWEpuJKvRu60xpCQh5o9paCKBau53H9wTuSZ3ztzMnJkzc7/vF1zuzPmeM/PlJJ/7O2fOzPwcEQKQz5/U3QCAehB+ICnCDyRF+IGkCD+QFOEHkiL8QFKEH0iK8ANJzermk9nm7YRAh0WEp7NeWyO/7ZW2f2N7n+1H23ksAN3lVt/bb/sySb+VdLukQ5J2SVobEe+VbMPID3RYN0b+myTti4j9EfFHST+WtLqNxwPQRe2Ef6Gkg5PuHyqWfY7tEdujtkfbeC4AFev4C34RsUnSJonDfqCXtDPyH5a0aNL9LxfLAPSBdsK/S9K1tr9i+4uS1kjaXk1bADqt5cP+iPjU9sOSfiHpMkmbI+LdyjoD0FEtX+pr6ck45wc6ritv8gHQvwg/kBThB5Ii/EBShB9IivADSRF+ICnCDyRF+IGkCD+QFOEHkiL8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kBThB5Ii/EBShB9IquUpuiXJ9gFJpyWdk/RpRAxX0RQ+74Ybbiitb9u2rWFtaGio4m56x4oVK0rrY2NjDWsHDx6sup2+01b4C7dFxIkKHgdAF3HYDyTVbvhD0i9t77Y9UkVDALqj3cP+5RFx2PZVkn5l+/8i4o3JKxR/FPjDAPSYtkb+iDhc/D4m6RVJN02xzqaIGObFQKC3tBx+23Nszz1/W9IKSXuragxAZ7Vz2L9A0iu2zz/Ov0fEzyvpCkDHtRz+iNgv6S8q7AUN3HHHHaX12bNnd6mT3nLXXXeV1h988MGGtTVr1lTdTt/hUh+QFOEHkiL8QFKEH0iK8ANJEX4gqSo+1Yc2zZpV/s+watWqLnXSX3bv3l1af+SRRxrW5syZU7rtxx9/3FJP/YSRH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeS4jp/D7jttttK67fccktp/bnnnquynb4xb9680vqSJUsa1q644orSbbnOD2DGIvxAUoQfSIrwA0kRfiApwg8kRfiBpBwR3Xsyu3tP1kOWLl1aWn/99ddL6x999FFp/cYbb2xYO3PmTOm2/azZflu+fHnD2uDgYOm2x48fb6WlnhARns56jPxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kFTTz/Pb3izpG5KORcTSYtmApJ9IGpJ0QNL9EfH7zrXZ35544onSerPvkF+5cmVpfaZeyx8YGCit33rrraX18fHxKtuZcaYz8v9Q0oX/+x6VtCMirpW0o7gPoI80DX9EvCHp5AWLV0vaUtzeIunuivsC0GGtnvMviIgjxe0PJS2oqB8AXdL2d/hFRJS9Z9/2iKSRdp8HQLVaHfmP2h6UpOL3sUYrRsSmiBiOiOEWnwtAB7Qa/u2S1he310t6tZp2AHRL0/Db3irpfyT9ue1Dth+S9Kyk222/L+mvi/sA+kjTc/6IWNug9PWKe+lb9957b2l91apVpfV9+/aV1kdHRy+5p5ng8ccfL603u45f9nn/U6dOtdLSjMI7/ICkCD+QFOEHkiL8QFKEH0iK8ANJMUV3Be67777SerPpoJ9//vkq2+kbQ0NDpfV169aV1s+dO1daf+aZZxrWPvnkk9JtM2DkB5Ii/EBShB9IivADSRF+ICnCDyRF+IGkuM4/TVdeeWXD2s0339zWY2/cuLGt7fvVyEj5t7vNnz+/tD42NlZa37lz5yX3lAkjP5AU4QeSIvxAUoQfSIrwA0kRfiApwg8kxXX+aZo9e3bD2sKFC0u33bp1a9XtzAiLFy9ua/u9e/dW1ElOjPxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kFTT6/y2N0v6hqRjEbG0WPakpG9JOl6s9lhE/KxTTfaC06dPN6zt2bOndNvrr7++tD4wMFBaP3nyZGm9l1111VUNa82mNm/mzTffbGv77KYz8v9Q0soplv9rRCwrfmZ08IGZqGn4I+INSf079ACYUjvn/A/bftv2ZtvzKusIQFe0Gv6NkhZLWibpiKTvNVrR9ojtUdujLT4XgA5oKfwRcTQizkXEuKQfSLqpZN1NETEcEcOtNgmgei2F3/bgpLvflMTHq4A+M51LfVslfU3SfNuHJP2TpK/ZXiYpJB2Q9O0O9gigA5qGPyLWTrH4hQ700tPOnj3bsPbBBx+UbnvPPfeU1l977bXS+oYNG0rrnbR06dLS+tVXX11aHxoaaliLiFZa+sz4+Hhb22fHO/yApAg/kBThB5Ii/EBShB9IivADSbndyy2X9GR2956si6677rrS+lNPPVVav/POO0vrZV8b3mknTpworTf7/1M2zbbtlno6b+7cuaX1ssuzM1lETGvHMvIDSRF+ICnCDyRF+IGkCD+QFOEHkiL8QFJc5+8By5YtK61fc801XerkYi+99FJb22/ZsqVhbd26dW099qxZzDA/Fa7zAyhF+IGkCD+QFOEHkiL8QFKEH0iK8ANJcaG0BzSb4rtZvZft37+/Y4/d7GvF9+5lLpkyjPxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kFTT6/y2F0l6UdICSSFpU0R83/aApJ9IGpJ0QNL9EfH7zrWKflT23fztfm8/1/HbM52R/1NJ/xARSyTdLOk7tpdIelTSjoi4VtKO4j6APtE0/BFxJCLeKm6fljQmaaGk1ZLOf03LFkl3d6pJANW7pHN+20OSvirp15IWRMSRovShJk4LAPSJab+33/aXJL0s6bsR8YfJ52sREY2+n8/2iKSRdhsFUK1pjfy2v6CJ4P8oIrYVi4/aHizqg5KOTbVtRGyKiOGIGK6iYQDVaBp+TwzxL0gai4gNk0rbJa0vbq+X9Gr17QHolOkc9v+lpL+V9I7t858tfUzSs5L+w/ZDkn4n6f7OtIh+VvbV8N382nhcrGn4I+JNSY0uyH692nYAdAvv8AOSIvxAUoQfSIrwA0kRfiApwg8kxVd3o6Muv/zylrc9e/ZshZ3gQoz8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU1/nRUQ888EDD2qlTp0q3ffrpp6tuB5Mw8gNJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUlznR0ft2rWrYW3Dhg0Na5K0c+fOqtvBJIz8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5CUm82RbnuRpBclLZAUkjZFxPdtPynpW5KOF6s+FhE/a/JYTMgOdFhEeDrrTSf8g5IGI+It23Ml7ZZ0t6T7JZ2JiH+ZblOEH+i86Ya/6Tv8IuKIpCPF7dO2xyQtbK89AHW7pHN+20OSvirp18Wih22/bXuz7XkNthmxPWp7tK1OAVSq6WH/ZyvaX5L0X5L+OSK22V4g6YQmXgd4WhOnBg82eQwO+4EOq+ycX5Jsf0HSTyX9IiIu+jRGcUTw04hY2uRxCD/QYdMNf9PDftuW9IKkscnBL14IPO+bkvZeapMA6jOdV/uXS/pvSe9IGi8WPyZpraRlmjjsPyDp28WLg2WPxcgPdFilh/1VIfxA51V22A9gZiL8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIrwA0kRfiApwg8k1e0puk9I+t2k+/OLZb2oV3vr1b4kemtVlb392XRX7Orn+S96cns0IoZra6BEr/bWq31J9NaqunrjsB9IivADSdUd/k01P3+ZXu2tV/uS6K1VtfRW6zk/gPrUPfIDqEkt4be90vZvbO+z/WgdPTRi+4Dtd2zvqXuKsWIatGO2905aNmD7V7bfL35POU1aTb09aftwse/22F5VU2+LbO+0/Z7td23/fbG81n1X0lct+63rh/22L5P0W0m3SzokaZektRHxXlcbacD2AUnDEVH7NWHbfyXpjKQXz8+GZPs5SScj4tniD+e8iPjHHuntSV3izM0d6q3RzNJ/pxr3XZUzXlehjpH/Jkn7ImJ/RPxR0o8lra6hj54XEW9IOnnB4tWSthS3t2jiP0/XNeitJ0TEkYh4q7h9WtL5maVr3XclfdWijvAvlHRw0v1D6q0pv0PSL23vtj1SdzNTWDBpZqQPJS2os5kpNJ25uZsumFm6Z/ZdKzNeV40X/C62PCJukPQ3kr5THN72pJg4Z+ulyzUbJS3WxDRuRyR9r85mipmlX5b03Yj4w+Ranftuir5q2W91hP+wpEWT7n+5WNYTIuJw8fuYpFc0cZrSS46enyS1+H2s5n4+ExFHI+JcRIxL+oFq3HfFzNIvS/pRRGwrFte+76bqq679Vkf4d0m61vZXbH9R0hpJ22vo4yK25xQvxMj2HEkr1HuzD2+XtL64vV7SqzX28jm9MnNzo5mlVfO+67kZryOi6z+SVmniFf8PJD1eRw8N+rpa0v8WP+/W3ZukrZo4DPxEE6+NPCTpTyXtkPS+pP+UNNBDvf2bJmZzflsTQRusqbflmjikf1vSnuJnVd37rqSvWvYb7/ADkuIFPyApwg8kRfiBpAg/kBThB5Ii/EBShB9IivADSf0/fhI1ni26LDgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('Original instance prediction: {}'.format(explanation['X_pred']))\n",
    "plt.imshow(explanation['X'].reshape(28, 28))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pertinent negative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pertinent negative prediction: 9\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f36453b5400>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADrpJREFUeJzt3X+MVfWZx/HPI1CMgEYkTIDiTi24BpG1OiEkO1lrFiurjdAoBv7YsFFLTWqyTVaz6JpIYjYxRtg0JjahkRQ3lbIEjCTW3XaJrmuykgEjisy2umQqv6cEtVND6Mg8+8ecMSNyvudy77n33Jnn/Uomc+957rn38Tgfzrn3e+75mrsLQDyXVN0AgGoQfiAowg8ERfiBoAg/EBThB4Ii/EBQhB8IivADQU1s5YuZGacTAk3m7lbL4xra85vZMjP7jZl9aGbrGnkuAK1l9Z7bb2YTJP1W0m2SjkjqkbTa3Q8m1mHPDzRZK/b8iyV96O6H3P1Pkn4haXkDzweghRoJ/xxJh0fdP5It+xIzW2tme81sbwOvBaBkTf/Az903SdokcdgPtJNG9vxHJc0ddf/r2TIAY0Aj4e+RNN/MvmFmX5O0StKuctoC0Gx1H/a7++dm9pCk/5A0QdJmd3+/tM4ANFXdQ311vRjv+YGma8lJPgDGLsIPBEX4gaAIPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8ERfiBoAg/EBThB4Ii/EBQhB8IivADQRF+ICjCDwRF+IGgCD8QFOEHgmrpFN1ojtmzZ+fWjh071sJOWuuSS9L7rqGhoRZ1Mjax5weCIvxAUIQfCIrwA0ERfiAowg8ERfiBoBoa5zezPkkDks5J+tzdu8poCl920003Jes7d+7MrS1ZsiS57okTJ5L1adOmJesDAwPJejMtXbo0We/t7c2tHT58uOx2xpwyTvK51d1PlfA8AFqIw34gqEbD75J+ZWb7zGxtGQ0BaI1GD/u73f2omc2U9Gsz+193f2P0A7J/FPiHAWgzDe353f1o9rtf0kuSFl/gMZvcvYsPA4H2Unf4zWyKmU0buS3pO5IOlNUYgOZq5LC/Q9JLZjbyPC+6+7+X0hWApqs7/O5+SNJflNgLctx+++3J+uTJk3Nr8+bNS65bNM7fzHH8iRPTf35F5zesWLEiWb/vvvtya6tWrUquGwFDfUBQhB8IivADQRF+ICjCDwRF+IGguHR3Gyga8lq5cmXdz/3mm2/WvW4ZFi/+ykmfX+jp6UmuW3TZ8UOHDiXr3d3dubUpU6Yk1/3ss8+S9fGAPT8QFOEHgiL8QFCEHwiK8ANBEX4gKMIPBMU4fxu49dZbk/VFixYl608//XRubfXq1cl1L7vssmR9z549yfqBA+nrt/T39+fW3D257pEjR5L1osuKL1iwILdW9N/NOD+AcYvwA0ERfiAowg8ERfiBoAg/EBThB4JinL8FFi5cmKxv3bo1Wf/444+T9f3799e97osvvpisF8nmbcg1c+bM3FpfX19Dr110SXOksecHgiL8QFCEHwiK8ANBEX4gKMIPBEX4gaAKx/nNbLOk70rqd/eF2bLpkrZJ6pTUJ+led08PKI9jl156abL++OOPJ+tF15B/5JFHkvVt27bl1q666qrkuvv27UvWb7755mT94YcfTtZTvT/wwAPJdTds2JCsz58/P1kfGhpK1qOrZc//M0nLzlu2TtJud58vaXd2H8AYUhh+d39D0unzFi+XtCW7vUXSipL7AtBk9b7n73D349ntE5I6SuoHQIs0fG6/u7uZ5V6MzczWSlrb6OsAKFe9e/6TZjZLkrLfuVdpdPdN7t7l7l11vhaAJqg3/Lskrclur5H0cjntAGiVwvCb2VZJ/yPpz83siJndL+kpSbeZ2QeSlmb3AYwhVnTt9FJfLPHZwFh2zz33JOvbt29P1ouufX/DDTdcdE/tYs2aNbm166+/Prlu0fkPDz74YLL++uuv59aWLTt/9PrLBgcHk/V25u7piyxkOMMPCIrwA0ERfiAowg8ERfiBoAg/EBSX7i7BypUrk/Vz584l688991yZ7bTU0qVLk/UtW7bk1jo7O5PrvvXWW8l60XZ95plncmtjeSivLOz5gaAIPxAU4QeCIvxAUIQfCIrwA0ERfiAoxvlrdMUVV+TWlixZ0tBz79q1K1kv+mrrxIn5/xs//fTTunoacddddyXrRV+r7e/PvciTVq1alVx3xowZyXpvb2+y/uqrrybr0bHnB4Ii/EBQhB8IivADQRF+ICjCDwRF+IGgGOev0eTJk3Nrc+bMSa67Y8eOZP3o0aN19TQi1VuRyy+/PFkvOgehqH7LLbfk1h599NHkukXf1y+65DnS2PMDQRF+ICjCDwRF+IGgCD8QFOEHgiL8QFCF4/xmtlnSdyX1u/vCbNl6Sd+X9PvsYY+5+y+b1WQ7GBgYyK0Vfa/8uuuuS9anT5+erJ8+fTpZP3v2bLKecvXVVyfrjY6lF22blAkTJiTrjV6rILpa9vw/k3Shycz/xd1vzH7GdfCB8agw/O7+hqT0rgfAmNPIe/6HzOxdM9tsZleW1hGAlqg3/D+R9E1JN0o6LmlD3gPNbK2Z7TWzvXW+FoAmqCv87n7S3c+5+5Ckn0panHjsJnfvcveuepsEUL66wm9ms0bd/Z4kvl4FjDG1DPVtlfRtSTPM7IikJyR928xulOSS+iT9oIk9AmiCwvC7++oLLH6+Cb20tTNnzuTWDh48mFz37rvvTtZfeeWVZH3mzJnJ+rp165L1lKlTpybrK1euTNavueaaZL2zszO3VvR9/SL79+9vaP3oOMMPCIrwA0ERfiAowg8ERfiBoAg/EJS5e+tezKx1L9ZCRV/ZffbZZ5P17u7uZH3SpEkX3VNZTp06lawX/f2kptk2s7p6GjFt2rRkPTU8O565e00blj0/EBThB4Ii/EBQhB8IivADQRF+ICjCDwTFOH8LFE2DXfS12Hnz5tX92tu3b697Xal4LL7o7+fYsWO5tY6Ojrp6GnHttdcm64cOHWro+ccqxvkBJBF+ICjCDwRF+IGgCD8QFOEHgiL8QFCM86Opnnjiidza+vXrk+sW/W0uWrQoWW90evGxinF+AEmEHwiK8ANBEX4gKMIPBEX4gaAIPxBU4RTdZjZX0guSOiS5pE3u/mMzmy5pm6ROSX2S7nX3j5vX6vg1cWL6f8OECROS9bNnz5bZTqkauTZ/0Th/1HH8stSy5/9c0j+4+wJJSyT90MwWSFonabe7z5e0O7sPYIwoDL+7H3f3t7PbA5J6Jc2RtFzSluxhWyStaFaTAMp3Ue/5zaxT0rck7ZHU4e7Hs9IJDb8tADBGFL7nH2FmUyXtkPQjd//D6Pdy7u555+2b2VpJaxttFEC5atrzm9kkDQf/5+6+M1t80sxmZfVZkvovtK67b3L3LnfvKqNhAOUoDL8N7+Kfl9Tr7htHlXZJWpPdXiPp5fLbA9AstRz2/6Wkv5X0npm9ky17TNJTkv7NzO6X9DtJ9zanxfFvaGgoWZ89e3ayfuLEidxa0XDZ4OBgst6ozs7Outdt5dfNIyoMv7u/KSlvsPavy20HQKtwhh8QFOEHgiL8QFCEHwiK8ANBEX4gqJpP70XzFI3zf/TRRy3qpHznzp2rqyZJZ86cSdaLvupc9PzRsecHgiL8QFCEHwiK8ANBEX4gKMIPBEX4gaAY50dT3Xnnnbm1Tz75JLnuk08+mawzjt8Y9vxAUIQfCIrwA0ERfiAowg8ERfiBoAg/EBTj/GjIvHnzkvWenp7c2saNG3NrkvTaa6/V1RNqw54fCIrwA0ERfiAowg8ERfiBoAg/EBThB4KyojnQzWyupBckdUhySZvc/cdmtl7S9yX9PnvoY+7+y4LnYsL1cWby5MnJempOgsHBwbLbgSR3t1oeV0v4Z0ma5e5vm9k0SfskrZB0r6Q/uvsztTZF+Mcfwt9+ag1/4Rl+7n5c0vHs9oCZ9Uqa01h7AKp2Ue/5zaxT0rck7ckWPWRm75rZZjO7MmedtWa218z2NtQpgFIVHvZ/8UCzqZL+S9I/u/tOM+uQdErDnwM8qeG3BvcVPAeH/eMMh/3tp9bD/pr2/GY2SdIOST93953ZC5x093PuPiTpp5IW19ssgNYrDL+ZmaTnJfW6+8ZRy2eNetj3JB0ovz0AzVLLp/3dkv5b0nuSRo7hHpO0WtKNGj7s75P0g+zDwdRzcdgPNFlpQ31lIvxA85X6nh/A+EP4gaAIPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8ERfiBoAg/EBThB4Ii/EBQhB8IqtVTdJ+S9LtR92dky9pRu/bWrn1J9FavMnv7s1of2NLv83/lxc32untXZQ0ktGtv7dqXRG/1qqo3DvuBoAg/EFTV4d9U8euntGtv7dqXRG/1qqS3St/zA6hO1Xt+ABWpJPxmtszMfmNmH5rZuip6yGNmfWb2npm9U/UUY9k0aP1mdmDUsulm9msz+yD7fcFp0irqbb2ZHc223TtmdkdFvc01s9fM7KCZvW9mf58tr3TbJfqqZLu1/LDfzCZI+q2k2yQdkdQjabW7H2xpIznMrE9Sl7tXPiZsZn8l6Y+SXnD3hdmypyWddvensn84r3T3f2yT3tbrImdublJveTNL/50q3HZlznhdhir2/Islfejuh9z9T5J+IWl5BX20PXd/Q9Lp8xYvl7Qlu71Fw388LZfTW1tw9+Pu/nZ2e0DSyMzSlW67RF+VqCL8cyQdHnX/iNprym+X9Csz22dma6tu5gI6Rs2MdEJSR5XNXEDhzM2tdN7M0m2z7eqZ8bpsfOD3Vd3ufpOkv5H0w+zwti358Hu2dhqu+Ymkb2p4GrfjkjZU2Uw2s/QOST9y9z+MrlW57S7QVyXbrYrwH5U0d9T9r2fL2oK7H81+90t6Se03+/DJkUlSs9/9FffzhXaauflCM0urDbZdO814XUX4eyTNN7NvmNnXJK2StKuCPr7CzKZkH8TIzKZI+o7ab/bhXZLWZLfXSHq5wl6+pF1mbs6bWVoVb7u2m/Ha3Vv+I+kODX/i/3+S/qmKHnL6ukbS/uzn/ap7k7RVw4eBgxr+bOR+SVdJ2i3pA0n/KWl6G/X2rxqezfldDQdtVkW9dWv4kP5dSe9kP3dUve0SfVWy3TjDDwiKD/yAoAg/EBThB4Ii/EBQhB8IivADQRF+ICjCDwT1/0I6xONf3WLaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('Pertinent negative prediction: {}'.format(explanation[mode + '_pred']))\n",
    "plt.imshow(explanation[mode].reshape(28, 28))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate pertinent positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 'PP'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# initialize TensorFlow session before model definition\n",
    "sess = tf.Session()\n",
    "K.set_session(sess)\n",
    "sess.run(tf.global_variables_initializer())\n",
    "\n",
    "# define models\n",
    "cnn = load_model('mnist_cnn.h5')\n",
    "ae = load_model('mnist_ae.h5')\n",
    "\n",
    "# initialize CEM explainer and explain instance\n",
    "cem = CEM(sess, cnn, mode, shape, kappa=kappa, beta=beta, feature_range=feature_range, \n",
    "          gamma=gamma, ae_model=ae, max_iterations=max_iterations, \n",
    "          c_init=c_init, c_steps=c_steps, learning_rate_init=lr, clip=clip, no_info_val=no_info_val)\n",
    "explanation = cem.explain(X, verbose=False)\n",
    "\n",
    "sess.close()\n",
    "K.clear_session()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pertinent positive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pertinent positive prediction: 4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f363fa95588>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAC5pJREFUeJzt3V+sXWWZx/HvM6jhj70AT+bQIp06QrzpBU5OGgLNhLnQMMSkeEOEm5oRjxeSjAkXEuZCkomJMerEK5NKG6sZUBN0aMhkUMjMYEgQDkShgGjHlNhy2tJUYrkgKn282KvjsXSvfbr32nvt9vl+kpOz9npX13qy0t9Zf9611xuZiaR6/qrvAiT1w/BLRRl+qSjDLxVl+KWiDL9UlOGXijL8UlGGXyrqXbPcWET4OKE0ZZkZ61luoiN/RNwcEa9ExIGIuGeSdUmarRj32f6IuAj4JfAR4BDwDHB7Zr7U8m888ktTNosj/zbgQGb+OjN/D3wX2DHB+iTN0CThvwr4zZrPh5p5fyEiliNiJSJWJtiWpI5N/YZfZu4CdoGn/dI8meTIfxi4es3n9zfzJJ0HJgn/M8C1EfGBiHgP8AlgXzdlSZq2sU/7M/OPEXEX8ChwEbAnM1/srDJJUzV2V99YG/OaX5q6mTzkI+n8Zfilogy/VJThl4oy/FJRhl8qyvBLRRl+qSjDLxVl+KWiDL9UlOGXijL8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqkowy8VZfilogy/VJThl4oy/FJRYw/RDRARB4GTwNvAHzNzqYuiJE3fROFv/ENmHu9gPZJmyNN+qahJw5/AjyLi2YhY7qIgSbMx6Wn/9sw8HBF/Dfw4In6RmU+sXaD5o+AfBmnORGZ2s6KI+4A3M/MrLct0szFJQ2VmrGe5sU/7I+KyiNhwehr4KLB/3PVJmq1JTvsXgR9GxOn1PJCZ/9VJVZKmrrPT/nVtzNP+sezevbu1/Y033hjadvfdd3ddjubc1E/7JZ3fDL9UlOGXijL8UlGGXyrK8EtF2dUnXWDs6pPUyvBLRRl+qSjDLxVl+KWiDL9UlOGXijL8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1SU4ZeK6mKUXo1w6aWXtrbfcMMNre2PPfZYl+WUceWVVw5tO3LkyAwrmU8e+aWiDL9UlOGXijL8UlGGXyrK8EtFGX6pqJHv7Y+IPcDHgGOZubWZdwXwPWALcBC4LTN/O3Jjvrdfmrou39v/LeDmM+bdAzyemdcCjzefJZ1HRoY/M58ATpwxewewt5neC9zacV2Spmzca/7FzFxtpo8Aix3VI2lGJn62PzOz7Vo+IpaB5Um3I6lb4x75j0bERoDm97FhC2bmrsxcysylMbclaQrGDf8+YGczvRN4uJtyJM3Kerr6HgRuAhaAo8AXgP8Avg9sBl5l0NV35k3Bs63Lrj5pytbb1Tcy/F2qGv5t27a1tj/99NMzqkQVdNnPL+kCZPilogy/VJThl4oy/FJRhl8qyld3z8DFF1/cdwnnpWuuuaa1/cCBAzOq5MLkkV8qyvBLRRl+qSjDLxVl+KWiDL9UlOGXivIrvTpv3Xnnna3t999//4wqmS9+pVdSK8MvFWX4paIMv1SU4ZeKMvxSUYZfKsp+fvVmcbF9iMejR4+2tm/atKm1/bXXXjvnmi4E9vNLamX4paIMv1SU4ZeKMvxSUYZfKsrwS0WN7OePiD3Ax4Bjmbm1mXcf8Gng9WaxezPzP0durGg//9LSUmv7ysrKjCqZL6+//npr+8LCQmv7JZdc0tr+1ltvnXNNF4Iu+/m/Bdx8lvn/lpnXNT8jgy9pvowMf2Y+AZyYQS2SZmiSa/67IuL5iNgTEZd3VpGkmRg3/N8APghcB6wCXx22YEQsR8RKRNS8sJXm1Fjhz8yjmfl2Zp4Cvglsa1l2V2YuZWb7XS9JMzVW+CNi45qPHwf2d1OOpFkZOUR3RDwI3AQsRMQh4AvATRFxHZDAQeAzU6xR0hT4fX5N1Y033ji07cknn2z9t9dff31r+1NPPTVWTRc6v88vqZXhl4oy/FJRhl8qyvBLRRl+qSi7+i5wGzZsaG0/efLkjCrRrNjVJ6mV4ZeKMvxSUYZfKsrwS0UZfqkowy8VZT+/pmrr1q1D2/bv9x0w02A/v6RWhl8qyvBLRRl+qSjDLxVl+KWiDL9U1Mj39kuT2LRp09A2+/n75ZFfKsrwS0UZfqkowy8VZfilogy/VJThl4oa+X3+iLga+DawCCSwKzO/HhFXAN8DtgAHgdsy87cj1uX3+aUpW+/3+dcT/o3Axsx8LiI2AM8CtwKfBE5k5pci4h7g8sz8/Ih1GX5pyjp7mUdmrmbmc830SeBl4CpgB7C3WWwvgz8Iks4T53TNHxFbgA8DPwUWM3O1aTrC4LJA0nli3c/2R8R7gYeAz2Xm7yL+fGaRmTnslD4iloHlSQuV1K11vcAzIt4NPAI8mplfa+a9AtyUmavNfYH/ycwPjViP1/zSlHV2zR+DQ/xu4OXTwW/sA3Y20zuBh8+1SEn9Wc/d/u3AT4AXgFPN7HsZXPd/H9gMvMqgq+/EiHV55C9mYWFhaNvx48dnWEkdnXX1dcnw12P4Z8/39ktqZfilogy/VJThl4oy/FJRhl8qyld3a6pOnTo1eiH1wiO/VJThl4oy/FJRhl8qyvBLRRl+qSjDLxVlP7+mavPmzUPbTpxoff0DDzzwQGv7HXfcMVZNGvDILxVl+KWiDL9UlOGXijL8UlGGXyrK8EtF+epu6QLjq7sltTL8UlGGXyrK8EtFGX6pKMMvFWX4paJGhj8iro6I/46IlyLixYj452b+fRFxOCJ+1vzcMv1yJXVl5EM+EbER2JiZz0XEBuBZ4FbgNuDNzPzKujfmQz7S1K33IZ+Rb/LJzFVgtZk+GREvA1dNVp6kvp3TNX9EbAE+DPy0mXVXRDwfEXsi4vIh/2Y5IlYiYmWiSiV1at3P9kfEe4H/Bb6YmT+IiEXgOJDAvzK4NPinEevwtF+asvWe9q8r/BHxbuAR4NHM/NpZ2rcAj2Tm1hHrMfzSlHX2xZ6ICGA38PLa4Dc3Ak/7OLD/XIuU1J/13O3fDvwEeAE4Pd7yvcDtwHUMTvsPAp9pbg62rcsjvzRlnZ72d8XwS9Pn9/kltTL8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqkowy8VNfIFnh07Dry65vNCM28ezWtt81oXWNu4uqztb9a74Ey/z/+OjUesZOZSbwW0mNfa5rUusLZx9VWbp/1SUYZfKqrv8O/qeftt5rW2ea0LrG1cvdTW6zW/pP70feSX1JNewh8RN0fEKxFxICLu6aOGYSLiYES80Iw83OsQY80waMciYv+aeVdExI8j4lfN77MOk9ZTbXMxcnPLyNK97rt5G/F65qf9EXER8EvgI8Ah4Bng9sx8aaaFDBERB4GlzOy9Tzgi/h54E/j26dGQIuLLwInM/FLzh/PyzPz8nNR2H+c4cvOUahs2svQn6XHfdTnidRf6OPJvAw5k5q8z8/fAd4EdPdQx9zLzCeDEGbN3AHub6b0M/vPM3JDa5kJmrmbmc830SeD0yNK97ruWunrRR/ivAn6z5vMh5mvI7wR+FBHPRsRy38WcxeKakZGOAIt9FnMWI0dunqUzRpaem303zojXXfOG3zttz8y/A/4R+GxzejuXcnDNNk/dNd8APshgGLdV4Kt9FtOMLP0Q8LnM/N3atj733Vnq6mW/9RH+w8DVaz6/v5k3FzLzcPP7GPBDBpcp8+To6UFSm9/Heq7n/2Xm0cx8OzNPAd+kx33XjCz9EPDvmfmDZnbv++5sdfW13/oI/zPAtRHxgYh4D/AJYF8PdbxDRFzW3IghIi4DPsr8jT68D9jZTO8EHu6xlr8wLyM3DxtZmp733dyNeJ2ZM/8BbmFwx///gH/po4Yhdf0t8PPm58W+awMeZHAa+AcG90Y+BbwPeBz4FfAYcMUc1fYdBqM5P88gaBt7qm07g1P654GfNT+39L3vWurqZb/5hJ9UlDf8pKIMv1SU4ZeKMvxSUYZfKsrwS0UZfqkowy8V9Sfh5A9Nh9Nu1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('Pertinent positive prediction: {}'.format(explanation[mode + '_pred']))\n",
    "plt.imshow(explanation[mode].reshape(28, 28))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clean up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.remove('mnist_cnn.h5')\n",
    "os.remove('mnist_ae.h5')"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}