{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "7fd843d4-25cd-fd9e-d5a4-7bf948d78f4a"
      },
      "source": [
        "Simple occlusion and saliency maps (WIP)\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "b8f23803-f7c9-c890-22dd-eab6672e749c"
      },
      "source": [
        "TODO: \n",
        "\n",
        " 1. set backend to theano and change array shapes, as the backprop stuff\n",
        "    is theano-only for now. Problem: theano runs too slow for a kaggle\n",
        "    kernel.\n",
        " 2. animated gif stuff maybe"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "19d8eb6b-9b71-029c-06e1-7e74975ae669"
      },
      "outputs": [],
      "source": [
        "#import os\n",
        "#os.environ['KERAS_BACKEND'] = 'theano'\n",
        "\n",
        "#from keras import backend as K\n",
        "#K.set_image_dim_ordering('th')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "76adcf6f-d109-e47d-d8f0-c7a26e986df5"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "\n",
        "from keras.utils import np_utils\n",
        "from keras.models import Sequential\n",
        "from keras.layers import Dense, Dropout, Flatten, Activation, Conv2D, MaxPooling2D\n",
        "from keras.optimizers import Adam\n",
        "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
        "from keras.preprocessing.image import ImageDataGenerator\n",
        "from keras.activations import relu\n",
        "\n",
        "%matplotlib inline\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "249daad0-7a0d-ad43-2bf9-4861ffac0d90"
      },
      "outputs": [],
      "source": [
        "train_file = \"../input/train.csv\"\n",
        "test_file = \"../input/test.csv\"\n",
        "\n",
        "mnist_dataset = np.loadtxt(train_file, skiprows=1, dtype='int', delimiter=',')\n",
        "\n",
        "val_split = 0.125\n",
        "n_raw = mnist_dataset.shape[0]\n",
        "n_val = int(n_raw * val_split + 0.5)\n",
        "n_train = n_raw - n_val\n",
        "\n",
        "np.random.shuffle(mnist_dataset)\n",
        "x_val, x_train = mnist_dataset[:n_val,1:], mnist_dataset[n_val:,1:]\n",
        "y_val, y_train = mnist_dataset[:n_val,0], mnist_dataset[n_val:,0]\n",
        "\n",
        "x_train = x_train.astype(\"float32\")/255.0\n",
        "x_val = x_val.astype(\"float32\")/255.0\n",
        "y_train = np_utils.to_categorical(y_train)\n",
        "y_val = np_utils.to_categorical(y_val)\n",
        "\n",
        "n_classes = y_train.shape[1]\n",
        "x_train = x_train.reshape(n_train, 28, 28, 1)\n",
        "x_val = x_val.reshape(n_val, 28, 28, 1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "363e8bc8-a731-01da-9a39-d3914ea7a5c8"
      },
      "outputs": [],
      "source": [
        "def create_model(activation=relu):\n",
        "    model = Sequential()\n",
        "\n",
        "    model.add(Conv2D(filters = 32, kernel_size = (3, 3), input_shape = (28, 28, 1), padding='same'))\n",
        "    model.add(Activation(activation))\n",
        "    model.add(Conv2D(filters = 32, kernel_size = (3, 3)))\n",
        "    model.add(Activation(activation))\n",
        "    model.add(MaxPooling2D(pool_size=(2,2)))\n",
        "    model.add(Dropout(0.25))\n",
        "\n",
        "    model.add(Conv2D(filters = 64, kernel_size = (3, 3)))\n",
        "    model.add(Activation(activation))\n",
        "    model.add(Conv2D(filters = 64, kernel_size = (3, 3)))\n",
        "    model.add(Activation(activation))\n",
        "    model.add(MaxPooling2D(pool_size=(2,2)))\n",
        "    model.add(Dropout(0.25))\n",
        "\n",
        "    model.add(Flatten())\n",
        "    model.add(Dense(128))\n",
        "    model.add(Activation(activation))\n",
        "    model.add(Dropout(0.5))\n",
        "    model.add(Dense(10))\n",
        "    model.add(Activation('softmax'))\n",
        "\n",
        "    return model\n",
        "\n",
        "model = create_model()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "17cdb0ec-ba92-6663-b329-0612f6c1c580"
      },
      "outputs": [],
      "source": [
        "datagen = ImageDataGenerator(zoom_range = 0.1,\n",
        "                            height_shift_range = 0.1,\n",
        "                            width_shift_range = 0.1,\n",
        "                            rotation_range = 20)\n",
        "\n",
        "model.compile(loss='categorical_crossentropy', optimizer = Adam(lr=1e-3), metrics = [\"accuracy\"])\n",
        "\n",
        "callbacks = [EarlyStopping(monitor='val_loss', min_delta=0, patience=2, verbose=2, mode='auto'),\n",
        "            ModelCheckpoint('mnist.h5', monitor='val_loss', save_best_only=True, verbose=0)]\n",
        "\n",
        "hist = model.fit_generator(datagen.flow(x_train, y_train, batch_size = 64),\n",
        "                           steps_per_epoch = n_train/100, #Take away 100 when not on Kaggle kernel\n",
        "                           epochs = 5, #Increase this when not on Kaggle kernel\n",
        "                           verbose = 2,  #verbose=1 outputs ETA, but doesn't work well in the cloud\n",
        "                           validation_data = (x_val, y_val),\n",
        "                           callbacks = callbacks)\n",
        "\n",
        "mnist_testset = np.loadtxt(test_file, skiprows=1, dtype='int', delimiter=',')\n",
        "x_test = mnist_testset.astype(\"float32\")/255.0\n",
        "n_test = x_test.shape[0]\n",
        "x_test = x_test.reshape(n_test, 28, 28, 1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "eadf5435-e7e1-480b-b763-7019bcc8e8a0"
      },
      "source": [
        "## Occlusion Map ##\n",
        "\n",
        "roughly after Zeiler et al. (2013): \"Visualizing and Understanding Convolutional Networks\", https://arxiv.org/abs/1311.2901"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "580a57a8-c8b2-cb5b-9962-a743f50eb5a4"
      },
      "outputs": [],
      "source": [
        "def iter_occlusion(image, size=8):\n",
        "\n",
        "    occlusion = np.full((size * 5, size * 5, 1), [0.5], np.float32)\n",
        "    occlusion_center = np.full((size, size, 1), [0.5], np.float32)\n",
        "    occlusion_padding = size * 2\n",
        "\n",
        "    # print('padding...')\n",
        "    image_padded = np.pad(image, ( \\\n",
        "                        (occlusion_padding, occlusion_padding), (occlusion_padding, occlusion_padding), (0, 0) \\\n",
        "                        ), 'constant', constant_values = 0.0)\n",
        "\n",
        "    for y in range(occlusion_padding, image.shape[0] + occlusion_padding, size):\n",
        "\n",
        "        for x in range(occlusion_padding, image.shape[1] + occlusion_padding, size):\n",
        "            tmp = image_padded.copy()\n",
        "\n",
        "            tmp[y - occlusion_padding:y + occlusion_center.shape[0] + occlusion_padding, \\\n",
        "                x - occlusion_padding:x + occlusion_center.shape[1] + occlusion_padding] \\\n",
        "                = occlusion\n",
        "\n",
        "            tmp[y:y + occlusion_center.shape[0], x:x + occlusion_center.shape[1]] = occlusion_center\n",
        "\n",
        "            yield x - occlusion_padding, y - occlusion_padding, \\\n",
        "                  tmp[occlusion_padding:tmp.shape[0] - occlusion_padding, occlusion_padding:tmp.shape[1] - occlusion_padding]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6436d590-9b96-3e41-ba3d-080b3f3af2a3"
      },
      "outputs": [],
      "source": [
        "i = 23  # for example\n",
        "data = x_train[i]\n",
        "correct_class = np.argmax(y_train[i])\n",
        "\n",
        "# input tensor for model.predict\n",
        "inp = data.reshape(1, 28, 28, 1)\n",
        "\n",
        "# image data for matplotlib's imshow\n",
        "img = data.reshape(28, 28)\n",
        "\n",
        "# occlusion\n",
        "img_size = img.shape[0]\n",
        "occlusion_size = 4\n",
        "\n",
        "# preview\n",
        "_ = plt.imshow(img)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "935a1b04-e5f7-1479-437d-702604d087fe"
      },
      "outputs": [],
      "source": [
        "print('occluding...')\n",
        "\n",
        "heatmap = np.zeros((img_size, img_size), np.float32)\n",
        "class_pixels = np.zeros((img_size, img_size), np.int16)\n",
        "\n",
        "from collections import defaultdict\n",
        "counters = defaultdict(int)\n",
        "\n",
        "for n, (x, y, img_float) in enumerate(iter_occlusion(data, size=occlusion_size)):   \n",
        "\n",
        "    X = img_float.reshape(1, 28, 28, 1)\n",
        "    out = model.predict(X)\n",
        "    print('#{}: {} @ {} (correct class: {})'.format(n, np.argmax(out), np.amax(out), out[0][correct_class]))\n",
        "    #print('x {} - {} | y {} - {}'.format(x, x + occlusion_size, y, y + occlusion_size))\n",
        "\n",
        "    heatmap[y:y + occlusion_size, x:x + occlusion_size] = out[0][correct_class]\n",
        "    class_pixels[y:y + occlusion_size, x:x + occlusion_size] = np.argmax(out)\n",
        "    counters[np.argmax(out)] += 1\n",
        "\n",
        "pred = model.predict(inp)\n",
        "print('Correct class: {}'.format(correct_class))\n",
        "print('Predicted class: {} (prob: {})'.format(np.argmax(pred), np.amax(out)))\n",
        "\n",
        "print('Predictions:')\n",
        "for class_id, count in counters.items():\n",
        "    print('{}: {}'.format(class_id, count))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "b0397aaa-115e-7aa6-8f80-0f88cd93eac2"
      },
      "source": [
        "heatmap (low value = important field) and predicted-class map"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "40f7c027-b751-74cc-94cd-9ba1c6b41e30"
      },
      "outputs": [],
      "source": [
        "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
        "from matplotlib.colors import BoundaryNorm\n",
        "from matplotlib.colorbar import ColorbarBase\n",
        "from matplotlib import cm\n",
        "\n",
        "fig = plt.figure(figsize=(8, 8))\n",
        "\n",
        "ax1 = plt.subplot(1, 2, 1, aspect='equal')\n",
        "hm = ax1.imshow(heatmap)\n",
        "\n",
        "ax2 = plt.subplot(1, 2, 2, aspect='equal')\n",
        "\n",
        "\n",
        "vals = np.unique(class_pixels).tolist()\n",
        "bounds = vals + [vals[-1] + 1]  # add an extra item for cosmetic reasons\n",
        "\n",
        "custom = cm.get_cmap('Vega20', len(bounds)) # discrete colors\n",
        "\n",
        "norm = BoundaryNorm(bounds, custom.N)\n",
        "\n",
        "cp = ax2.imshow(class_pixels, norm=norm, cmap=custom)\n",
        "\n",
        "divider = make_axes_locatable(ax1)\n",
        "cax1 = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
        "cbar1 = plt.colorbar(hm, cax=cax1)\n",
        "\n",
        "divider = make_axes_locatable(ax2)\n",
        "cax2 = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
        "cbar2 = ColorbarBase(cax2, cmap=custom, norm=norm,\n",
        "                         # place the ticks at the average value between two entries\n",
        "                         # e.g. [280, 300] -> 290\n",
        "                         # so that they're centered on the colorbar\n",
        "                         ticks=[(a + b) / 2.0 for a, b in zip(bounds[::], bounds[1::])],\n",
        "                         boundaries=bounds, spacing='uniform', orientation='vertical')\n",
        "\n",
        "cbar2.ax.set_yticklabels([n for n in np.unique(class_pixels)])\n",
        "\n",
        "fig.tight_layout()\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "6d9a5539-bce1-ed93-b16e-2045fe322730"
      },
      "source": [
        "heatmap overlay"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "d7e0b2d1-1d47-0eb7-4d0e-8ca529c27916"
      },
      "outputs": [],
      "source": [
        "import cv2\n",
        "\n",
        "plt.figure(figsize=(6, 6))\n",
        "\n",
        "plt.imshow(img, cmap=cm.gray)\n",
        "\n",
        "plt.pcolormesh(heatmap, cmap=plt.cm.jet, alpha=0.50)\n",
        "plt.colorbar().solids.set(alpha=1)\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "4fb77903-cd36-6611-f6e7-58b11d5e11c7"
      },
      "source": [
        "## Saliency Maps ##\n",
        "\n",
        "after https://github.com/Lasagne/Recipes/blob/master/examples/Saliency%20Maps%20and%20Guided%20Backpropagation.ipynb"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "fdce84d9-bcc6-748a-e05c-02fd26a56ce3"
      },
      "outputs": [],
      "source": [
        "import theano\n",
        "from keras import backend as K\n",
        "\n",
        "def compile_saliency_function(model):\n",
        "    \"\"\"\n",
        "    Compiles a function to compute the saliency maps and predicted classes\n",
        "    for a given minibatch of input images.\n",
        "    \"\"\"\n",
        "    inp = model.layers[0].input\n",
        "    outp = model.layers[-1].output\n",
        "    max_outp = K.T.max(outp, axis=1)\n",
        "    saliency = K.gradients(K.sum(max_outp), inp)\n",
        "    max_class = K.T.argmax(outp, axis=1)\n",
        "    return K.function([inp, K.learning_phase()], [saliency, max_class])\n",
        "\n",
        "def show_images(img_original, saliency, max_class, title):\n",
        "    # get out the first map and class from the mini-batch\n",
        "    saliency = saliency[0]\n",
        "    saliency = saliency.reshape(28, 28, 1)\n",
        "    #saliency = saliency[::-1].transpose(1, 2, 0)\n",
        "    \n",
        "    max_class = max_class[0]\n",
        "\n",
        "    # plot the original image and the three saliency map variants\n",
        "    plt.figure(figsize=(10, 10), facecolor='w')\n",
        "    plt.suptitle(str(max_class) + \" Saliency: \" + title)\n",
        "    plt.subplot(2, 2, 1)\n",
        "    plt.title('input')\n",
        "    plt.imshow(img_original)\n",
        "    plt.subplot(2, 2, 2)\n",
        "    plt.title('abs. saliency')\n",
        "    plt.imshow(np.abs(saliency).max(axis=2), cmap='gray')\n",
        "    plt.subplot(2, 2, 3)\n",
        "    plt.title('pos. saliency')\n",
        "    x = (np.maximum(0, saliency) / saliency.max()).reshape((28, 28))\n",
        "    plt.imshow(x)\n",
        "    plt.subplot(2, 2, 4)\n",
        "    plt.title('neg. saliency')\n",
        "    x = (np.maximum(0, -saliency) / -saliency.min()).reshape((28, 28))\n",
        "    plt.imshow(x)\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "92ae01d9-a73f-2c5c-3783-69494167a557"
      },
      "source": [
        "Simonyan et al. (2013): Plain Gradient\n",
        "---------------------------------------------\n",
        "\n",
        "**needs theano**\n",
        "\n",
        "Simonyan et al. (2013): \"Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps\", http://arxiv.org/abs/1312.6034"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "76b82e99-8035-da9f-c4df-597000f1b8d9"
      },
      "outputs": [],
      "source": [
        "saliency, max_class = compile_saliency_function(model)([inp, 0])\n",
        "show_images(img, saliency, max_class, \"plain gradient\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "c4c9dbd6-aea1-31b0-8db9-eb6597be1614"
      },
      "source": [
        "Springenberg et al. (2015): Guided Backpropagation\n",
        "--------------------------------------------------------------\n",
        "**needs theano**\n",
        "\n",
        "Springenberg et al. (2015): \"Striving for Simplicity - The All Convolutional Net\", http://arxiv.org/abs/1412.6806"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a1b67fc1-37fe-ca39-638e-a7677ed83db7"
      },
      "outputs": [],
      "source": [
        "class ModifiedBackprop(object):\n",
        "\n",
        "    def __init__(self, nonlinearity):\n",
        "        self.nonlinearity = nonlinearity\n",
        "        self.ops = {}  # memoizes an OpFromGraph instance per tensor type\n",
        "\n",
        "    def __call__(self, x):\n",
        "        # We note the tensor type of the input variable to the nonlinearity\n",
        "        # (mainly dimensionality and dtype); we need to create a fitting Op.\n",
        "        tensor_type = x.type\n",
        "        # If we did not create a suitable Op yet, this is the time to do so.\n",
        "        if tensor_type not in self.ops:\n",
        "            # For the graph, we create an input variable of the correct type:\n",
        "            inp = tensor_type()\n",
        "            # We pass it through the nonlinearity.\n",
        "            outp = self.nonlinearity(inp)\n",
        "            # Then we fix the forward expression...\n",
        "            op = theano.OpFromGraph([inp], [outp])\n",
        "            # ...and replace the gradient with our own (defined in a subclass).\n",
        "            op.grad = self.grad\n",
        "            # Finally, we memoize the new Op\n",
        "            self.ops[tensor_type] = op\n",
        "        # And apply the memoized Op to the input we got.\n",
        "        return self.ops[tensor_type](x)\n",
        "\n",
        "class GuidedBackprop(ModifiedBackprop):\n",
        "    __name__ = 'GuidedBackprop'\n",
        "\n",
        "    def grad(self, inputs, out_grads):\n",
        "        (inp,) = inputs\n",
        "        (grd,) = out_grads\n",
        "        dtype = inp.dtype\n",
        "        return (grd * (inp > 0).astype(dtype) * (grd > 0).astype(dtype),)\n",
        "\n",
        "\n",
        "modded_relu = GuidedBackprop(relu)\n",
        "model_guided = create_model(activation=modded_relu)\n",
        "model_guided.compile(optimizer=Adam(lr=1e-3), loss='categorical_crossentropy')\n",
        "\n",
        "saliency, max_class = compile_saliency_function(model_guided)([inp, 0])\n",
        "show_images(img, saliency, max_class, \"guided backprop\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "9c0a6f39-0693-b804-cd3e-460379731e0a"
      },
      "source": [
        "Zeiler et al. (2013): DeconvNet\n",
        "-------------------------------\n",
        "**needs theano**\n",
        "\n",
        "Zeiler et al. (2013): \"Visualizing and Understanding Convolutional Networks\", http://arxiv.org/abs/1311.2901"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "2aa3ce58-fbfb-b7e9-1bdf-653d4877cbcd"
      },
      "outputs": [],
      "source": [
        "class ZeilerBackprop(ModifiedBackprop):\n",
        "    __name__ = 'ZeilerBackprop'\n",
        "\n",
        "    def grad(self, inputs, out_grads):\n",
        "        (inp,) = inputs\n",
        "        (grd,) = out_grads\n",
        "        #return (grd * (grd > 0).astype(inp.dtype),)  # explicitly rectify\n",
        "        return (self.nonlinearity(grd),)  # use the given nonlinearity\n",
        "\n",
        "\n",
        "modded_relu = ZeilerBackprop(relu)\n",
        "model_zeiler = create_model(activation=modded_relu)\n",
        "model_zeiler.compile(optimizer=Adam(lr=1e-3), loss='categorical_crossentropy')\n",
        "\n",
        "saliency, max_class = compile_saliency_function(model_zeiler)([inp, 0])\n",
        "show_images(img, saliency, max_class, \"deconvnet\")"
      ]
    }
  ],
  "metadata": {
    "_change_revision": 0,
    "_is_fork": false,
    "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.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}