whale-recognition-model-with-score-0-78563.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "baff2752922fc86a6d62d3c8e02dd5d868a9d296",
    "collapsed": true
   },
   "source": [
    "# Whale Classification Model\n",
    "This notebook describes the strategy behind the 0.78563 submission to the Humpack Whale identification Challenge.\n",
    "\n",
    "It should be studied in conjunction with the [Bounding Box Model](http://www.kaggle.com/martinpiotte/bounding-box-model) notebook which describes separately the strategy for image cropping.\n",
    "\n",
    "To speed things up, the results of some slow computations are included as a dataset instead of being recomputed here. However, the code is still provided in the notebook as reference, even if it is not executed by default."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "3c7ec1841dc33ba31d374ab9f45cab51ba5b9886"
   },
   "source": [
    "# Abstract\n",
    "The approach used for this submission is essentially a [Siamese Neural Network](http://www.cs.utoronto.ca/~gkoch/files/msc-thesis.pdf), with a few modifications that will be covered in details later. The element that generated the largest accuracy improvement is the procedure used to generate image pairs during training. Each training epoch is composed of a sequence of image pairs (A, B) such that:\n",
    "\n",
    "* Exactly 50% of the pairs are for matching whales, and 50% for different whales;\n",
    "* Each image from the training set is used exactly 4 times per epoch: A and B images of matching whales, A and B images of different whale pairs;\n",
    "* Pairs of images of different whales are selected to be difficult for the network to distinguish at a given stage of the training. This is inspired from adversarial training: find pairs of images that are from different whales, but that are still very similar from the model perspective.\n",
    "\n",
    "Implementing this strategy while training a Siamese Neural Network is what makes the largest contribution to the model accuracy. Other details contribute somewhat to the accuracy, but have a much smaller impact."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "bf399701be9c6b8d72fc05cccdefdc4f5eeb2a89"
   },
   "source": [
    "# Overview\n",
    "This notebook describes all the different elements of the submission. Obviously, to cover everything, it has to be fairly long. I encourage everyone to skip ahead directly  to whatever you are most interested in, without  necessarily going through everything.\n",
    "## Content\n",
    "1. Duplicate image identification (not much to see here -- keep moving)\n",
    "1. Image preprocessing (just the regular stuff)\n",
    "1. Siamese Neural Network architecture (some interesting thoughts)\n",
    "1. Training data construction (most of the secret sauce is here)\n",
    "1. Training procedure (zzzzz.....)\n",
    "1. Generating the submission file (re-zzzzz.....)\n",
    "1. Bootstrapping and ensemble (classic but short)\n",
    "1. Visualization (everyone's favorite!)\n",
    "1. Off topic (why add this unless it is interesting?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "7afad70768ff847794be2ab8312d7fa04ca1c571",
    "collapsed": true
   },
   "source": [
    "# Duplicate image identification\n",
    "This section describes the heuristic used to identify duplicate images. The fact that the training and test set have duplicate images has already been well documented. Some images are perfect binary copies, while other have been altered somewhat: contrast and brightness, size, masking the legend, etc. \n",
    "\n",
    "Two images are considered duplicate if they meet the following criteria:\n",
    "\n",
    "1. Both images have the same [Perceptual Hash](http://www.phash.org/) (phash); or\n",
    "1. Both images have:\n",
    "    1. phash that differ by at most 6 bits, and;\n",
    "    1. have the same size, and;\n",
    "    1. the pixelwise mean square error between the normalized images is below a given threshold.\n",
    "\n",
    "The *p2h* dictionary associate a unique image id (phash) for each picture. The *h2p* dictionary associate each unique image id to the prefered image to be used for this hash.\n",
    "\n",
    "The prefered image is the one with the highest resolution, or any one if they have the same resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "_uuid": "34d1be7393ae2c49babbb540fc766fbc3274c579"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25460\n",
      "1624\n"
     ]
    }
   ],
   "source": [
    "# Read the dataset description\n",
    "from pandas import read_csv\n",
    "\n",
    "tagged = dict([(p,w) for _,p,w in read_csv('./train.csv').to_records()])\n",
    "submit = [p for _,p,_ in read_csv('./sample_submission.csv').to_records()]\n",
    "join   = list(tagged.keys()) + submit\n",
    "len(tagged),len(submit),len(join),list(tagged.items())[:5],submit[:5]\n",
    "join2 = []\n",
    "cc = 0\n",
    "tagged2 = {}\n",
    "print(len(join))\n",
    "for j in join:\n",
    "    if expand_path(j) != None:\n",
    "        join2.append(j)\n",
    "        if j in tagged:\n",
    "            tagged2[j] = tagged[j]\n",
    "print(len(join2))\n",
    "join = join2\n",
    "tagged = tagged2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "_uuid": "ab6d0c498190ba97c62b96cca930a1ba2201f987"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "90b98d12ef9848a69f3bc8f073b9f643",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1624), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1624,\n",
       " [('997cc007.jpg', (1050, 459)),\n",
       "  ('90d65ce2.jpg', (1050, 600)),\n",
       "  ('f976672a.jpg', (1050, 450)),\n",
       "  ('fa9c24f0.jpg', (850, 325)),\n",
       "  ('91fae2d5.jpg', (615, 292))])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determise the size of each image\n",
    "from os.path import isfile\n",
    "from PIL import Image as pil_image\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "def expand_path(p):\n",
    "    if isfile('./train/' + p): return './train/' + p\n",
    "    if isfile('./test/' + p): return './test/' + p\n",
    "    return None\n",
    "\n",
    "p2size = {}\n",
    "for p in tqdm_notebook(join):\n",
    "    size      = pil_image.open(expand_path(p)).size\n",
    "    p2size[p] = size\n",
    "len(p2size), list(p2size.items())[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: imagehash in /home/vahid/.virtualenvs/tensorflow/lib/python3.5/site-packages (4.0)\n",
      "Requirement already satisfied: six in /home/vahid/.virtualenvs/tensorflow/lib/python3.5/site-packages (from imagehash) (1.11.0)\n",
      "Requirement already satisfied: numpy in /home/vahid/.virtualenvs/tensorflow/lib/python3.5/site-packages (from imagehash) (1.14.1)\n",
      "Requirement already satisfied: pywavelets in /home/vahid/.virtualenvs/tensorflow/lib/python3.5/site-packages (from imagehash) (1.0.1)\n",
      "Requirement already satisfied: pillow in /home/vahid/.virtualenvs/tensorflow/lib/python3.5/site-packages (from imagehash) (5.1.0)\n",
      "Requirement already satisfied: scipy in /home/vahid/.virtualenvs/tensorflow/lib/python3.5/site-packages (from imagehash) (1.1.0)\n",
      "\u001b[33mYou are using pip version 18.0, however version 19.0.2 is available.\n",
      "You should consider upgrading via the 'pip install --upgrade pip' command.\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "!pip install imagehash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "_uuid": "885c2c1714e2ef38c2c30ec441684ec853d0e1a0"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "THERE IS'NT THE FILE\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8c167cf1cb164537956110b547743d2c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1624), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "daf688dd08fd4389a5607b7eb4bfd9d6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1612), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1624,\n",
       " [('997cc007.jpg', 'bfcc926c633c093a'),\n",
       "  ('90d65ce2.jpg', 'aa15c62d1ad4b765'),\n",
       "  ('f976672a.jpg', 'bfc0d025685e476b'),\n",
       "  ('fa9c24f0.jpg', 'ae4fc0b1248dcf72'),\n",
       "  ('91fae2d5.jpg', 'be4ad2e08fb334b0')])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read or generate p2h, a dictionary of image name to image id (picture to hash)\n",
    "import pickle\n",
    "import numpy as np\n",
    "from imagehash import phash\n",
    "from math import sqrt\n",
    "\n",
    "# Two phash values are considered duplicate if, for all associated image pairs:\n",
    "# 1) They have the same mode and size;\n",
    "# 2) After normalizing the pixel to zero mean and variance 1.0, the mean square error does not exceed 0.1\n",
    "def match(h1,h2):\n",
    "    for p1 in h2ps[h1]:\n",
    "        for p2 in h2ps[h2]:\n",
    "            i1 =  pil_image.open(expand_path(p1))\n",
    "            i2 =  pil_image.open(expand_path(p2))\n",
    "            if i1.mode != i2.mode or i1.size != i2.size: return False\n",
    "            a1 = np.array(i1)\n",
    "            a1 = a1 - a1.mean()\n",
    "            a1 = a1/sqrt((a1**2).mean())\n",
    "            a2 = np.array(i2)\n",
    "            a2 = a2 - a2.mean()\n",
    "            a2 = a2/sqrt((a2**2).mean())\n",
    "            a  = ((a1 - a2)**2).mean()\n",
    "            if a > 0.1: return False\n",
    "    return True\n",
    "\n",
    "if isfile('./1p2h.pickle'):\n",
    "    print(\"THERE IS THE FILE\")\n",
    "    with open('./p2h.pickle', 'rb') as f:\n",
    "        p2h = pickle.load(f)\n",
    "        # Find all images associated with a given phash value.\n",
    "    h2ps = {}\n",
    "    for p,h in p2h.items():\n",
    "        if h not in h2ps: h2ps[h] = []\n",
    "        if p not in h2ps[h]: h2ps[h].append(p)\n",
    "\n",
    "else:\n",
    "    print(\"THERE IS'NT THE FILE\")\n",
    "    # Compute phash for each image in the training and test set.\n",
    "    p2h = {}\n",
    "    for p in tqdm_notebook(join):\n",
    "        img    = pil_image.open(expand_path(p))\n",
    "        h      = phash(img)\n",
    "        p2h[p] = h\n",
    "\n",
    "    # Find all images associated with a given phash value.\n",
    "    h2ps = {}\n",
    "    for p,h in p2h.items():\n",
    "        if h not in h2ps: h2ps[h] = []\n",
    "        if p not in h2ps[h]: h2ps[h].append(p)\n",
    "\n",
    "    # Find all distinct phash values\n",
    "    hs = list(h2ps.keys())\n",
    "\n",
    "    # If the images are close enough, associate the two phash values (this is the slow part: n^2 algorithm)\n",
    "    h2h = {}\n",
    "    for i,h1 in enumerate(tqdm_notebook(hs)):\n",
    "        for h2 in hs[:i]:\n",
    "            if h1-h2 <= 6 and match(h1, h2):\n",
    "                s1 = str(h1)\n",
    "                s2 = str(h2)\n",
    "                if s1 < s2: s1,s2 = s2,s1\n",
    "                h2h[s1] = s2\n",
    "\n",
    "    # Group together images with equivalent phash, and replace by string format of phash (faster and more readable)\n",
    "    for p,h in p2h.items():\n",
    "        h = str(h)\n",
    "        if h in h2h: h = h2h[h]\n",
    "        p2h[p] = h\n",
    "\n",
    "len(p2h), list(p2h.items())[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "_uuid": "61e410583127029513732341a99145151e329475"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1612,\n",
       " [('ea9681fae7607849', ['9931ee1b.jpg']),\n",
       "  ('b4918ffb2694c0ad', ['93d859d2.jpg']),\n",
       "  ('b6c8d2b7ad289c0b', ['90c49192.jpg']),\n",
       "  ('bfc0d025685e476b', ['f976672a.jpg']),\n",
       "  ('ba4fd0f6a230e589', ['f8bcfd01.jpg'])])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For each image id, determine the list of pictures\n",
    "h2ps = {}\n",
    "for p,h in p2h.items():\n",
    "    if h not in h2ps: h2ps[h] = []\n",
    "    if p not in h2ps[h]: h2ps[h].append(p)\n",
    "# Notice how 25460 images use only 20913 distinct image ids.\n",
    "len(h2ps),list(h2ps.items())[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "_uuid": "7e768fbcafade78a656c42a1db326cdd7c24478d"
   },
   "outputs": [],
   "source": [
    "# Show an example of a duplicate image (from training of test set)\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def show_whale(imgs, per_row=2):\n",
    "    n         = len(imgs)\n",
    "    rows      = (n + per_row - 1)//per_row\n",
    "    cols      = min(per_row, n)\n",
    "    fig, axes = plt.subplots(rows,cols, figsize=(24//per_row*cols,24//per_row*rows))\n",
    "    for ax in axes.flatten(): ax.axis('off')\n",
    "    for i,(img,ax) in enumerate(zip(imgs, axes.flatten())): ax.imshow(img.convert('RGB'))\n",
    "\n",
    "for h, ps in h2ps.items():\n",
    "    if len(ps) > 2:\n",
    "        print('Images:', ps)\n",
    "        imgs = [pil_image.open(and_path(p)) for p in ps]\n",
    "        show_whale(imgs, per_row=len(ps))\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "_cell_guid": "79c7e3d0-c299-4dcb-8224-4455121ee9b0",
    "_uuid": "d629ff2d2480ee46fbb7e2d37f6b5fab8052498a"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1612,\n",
       " [('ea9681fae7607849', '9931ee1b.jpg'),\n",
       "  ('b4918ffb2694c0ad', '93d859d2.jpg'),\n",
       "  ('b6c8d2b7ad289c0b', '90c49192.jpg'),\n",
       "  ('bfc0d025685e476b', 'f976672a.jpg'),\n",
       "  ('ee9994a469781e4b', '96aaafd8.jpg')])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For each images id, select the prefered image\n",
    "def prefer(ps):\n",
    "    if len(ps) == 1: return ps[0]\n",
    "    best_p = ps[0]\n",
    "    if not best_p in p2size:\n",
    "        return ps\n",
    "    best_s = p2size[best_p]\n",
    "    for i in range(1, len(ps)):\n",
    "        p = ps[i]\n",
    "        if not p in p2size:\n",
    "            continue\n",
    "        s = p2size[p]\n",
    "        if s[0]*s[1] > best_s[0]*best_s[1]: # Select the image with highest resolution\n",
    "            best_p = p\n",
    "            best_s = s\n",
    "    return best_p\n",
    "\n",
    "h2p = {}\n",
    "for h,ps in h2ps.items(): h2p[h] = prefer(ps)\n",
    "len(h2p),list(h2p.items())[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "bd6640560036cf8e6d581a204e4a78a96f465d36",
    "collapsed": true
   },
   "source": [
    "# Image preprocessing\n",
    "Training is performed on images subjected to the following operations:\n",
    "\n",
    "1.  Rotate the image if it is in the rotate set;\n",
    "1. Transform to black and white;\n",
    "1. Apply an affine transformation ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "aec39ab5192624f08c55a330cf04804b0d63b7c1"
   },
   "source": [
    "## Image rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "5275aba3332f8600731d19e9b9a5c3cea696e415"
   },
   "source": [
    "I noticed that some pictures have the whale fluke pointing down instead of up as usual. Whenever I encountered such instance in the training set (not in the test set), I would add it to a list. During training, these images are rotated 180 degrees to normalize them with the fluke pointing up. The list is not exhausitve, there are probably more case that I have not noticed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "_uuid": "237fba51c8e460e2c362a8876ef963a37b60bbab"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'2b792814.jpg',\n",
       " '2bc459eb.jpg',\n",
       " '3401bafe.jpg',\n",
       " '56fafc52.jpg',\n",
       " 'a492ab72.jpg',\n",
       " 'd1502267.jpg',\n",
       " 'e53d2b96.jpg',\n",
       " 'ed4f0cd5.jpg',\n",
       " 'f2ec136c.jpg',\n",
       " 'f966c073.jpg'}"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with open('./rotate.txt', 'rt') as f: rotate = f.read().split('\\n')[:-1]\n",
    "rotate = set(rotate)\n",
    "rotate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "_uuid": "b4687ab632f33f5afc7d2f813a3be441927cd9ae"
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'NoneType' object has no attribute 'read'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m~/.virtualenvs/tensorflow/lib/python3.5/site-packages/PIL/Image.py\u001b[0m in \u001b[0;36mopen\u001b[0;34m(fp, mode)\u001b[0m\n\u001b[1;32m   2551\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2552\u001b[0;31m         \u001b[0mfp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mseek\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2553\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mAttributeError\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mUnsupportedOperation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'seek'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-47-f63e625a249e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mp\u001b[0m    \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrotate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mimgs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mpil_image\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpand_path\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mread_raw_image\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m \u001b[0mshow_whale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimgs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.virtualenvs/tensorflow/lib/python3.5/site-packages/PIL/Image.py\u001b[0m in \u001b[0;36mopen\u001b[0;34m(fp, mode)\u001b[0m\n\u001b[1;32m   2552\u001b[0m         \u001b[0mfp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mseek\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2553\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mAttributeError\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mUnsupportedOperation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2554\u001b[0;31m         \u001b[0mfp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mBytesIO\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2555\u001b[0m         \u001b[0mexclusive_fp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2556\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'read'"
     ]
    }
   ],
   "source": [
    "def read_raw_image(p):\n",
    "    img = pil_image.open(expand_path(p))\n",
    "    if p in rotate: img = img.rotate(180)\n",
    "    return img\n",
    "\n",
    "p    = list(rotate)[0]\n",
    "imgs = [pil_image.open(expand_path(p)), read_raw_image(p)]\n",
    "show_whale(imgs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "b9fb058e68b8d9c87426e46670e06e29fab8a85a"
   },
   "source": [
    "The left side image is the original image (fluke pointing down). The right side image is rotated 180 degrees."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "53bfdd79fc980891a85a6bd93dc54ef5b25a8d98"
   },
   "source": [
    "## Convert to black and white\n",
    "In my early experiments, I noticed that my models achieved approximately the same accuracy when comparing two colored images, or two black and white images. However, comparing a colored image with a black and white image resulted in much lower accuracy. The simplest solution was to convert all images to black and white, which did not reduce the accuracy even when comparing originally colored images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "c44ac944607ebb8bfed66324e60440fbda4bb8fc"
   },
   "source": [
    "## Affine tranformation\n",
    "The [affine transformation](https://en.wikipedia.org/wiki/Affine_transformation) maps a rectangular area of the original image to a square image with resolution 384x384x1 (only one channel for black and white). The rectangular area has a width over height aspect ratio of 2.15, close to the average picture aspect ratio. The rectangle is taken to be slightly larger than the computed [bounding box](https://www.kaggle.com/martinpiotte/bounding-box-model), as computed in the model defined in a separate kernel, The idea is that clipping the edges of the fluke is more harmful than the gain obtained by fitting it exactly, thus a margin is preferred.\n",
    "\n",
    "During training, data augmentation is performed by adding a random transformation that composes zoom, shift, rotation and shear. The random transform is skipped when testing.\n",
    "\n",
    "Finally, the image is normalized to zero mean and unit variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "_uuid": "b1b9bff9cbec6d6120abd1d11a94929481780acb"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('28781a20.jpg', (44, 135, 957, 396)),\n",
       " ('5775c9f2.jpg', (15, 183, 1045, 356)),\n",
       " ('0845b7cb.jpg', (22, 139, 1013, 418)),\n",
       " ('6c13f3de.jpg', (5, 135, 662, 294)),\n",
       " ('7162b783.jpg', (35, 129, 653, 236))]"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read the bounding box data from the bounding box kernel (see reference above)\n",
    "with open('./bounding-box.pickle', 'rb') as f:\n",
    "    p2bb = pickle.load(f)\n",
    "list(p2bb.items())[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "_uuid": "97c4192e165c301719289053cd8e66ce1cd2367d"
   },
   "outputs": [],
   "source": [
    "# Suppress annoying stderr output when importing keras.\n",
    "import sys\n",
    "import platform\n",
    "old_stderr = sys.stderr\n",
    "sys.stderr = open('/dev/null' if platform.system() != 'Windows' else 'nul', 'w')\n",
    "import keras\n",
    "sys.stderr = old_stderr\n",
    "\n",
    "import random\n",
    "from keras import backend as K\n",
    "from keras.preprocessing.image import img_to_array,array_to_img\n",
    "from scipy.ndimage import affine_transform\n",
    "\n",
    "img_shape    = (384,384,1) # The image shape used by the model\n",
    "anisotropy   = 2.15 # The horizontal compression ratio\n",
    "crop_margin  = 0.05 # The margin added around the bounding box to compensate for bounding box inaccuracy\n",
    "\n",
    "def build_transform(rotation, shear, height_zoom, width_zoom, height_shift, width_shift):\n",
    "    \"\"\"\n",
    "    Build a transformation matrix with the specified characteristics.\n",
    "    \"\"\"\n",
    "    rotation        = np.deg2rad(rotation)\n",
    "    shear           = np.deg2rad(shear)\n",
    "    rotation_matrix = np.array([[np.cos(rotation), np.sin(rotation), 0], [-np.sin(rotation), np.cos(rotation), 0], [0, 0, 1]])\n",
    "    shift_matrix    = np.array([[1, 0, height_shift], [0, 1, width_shift], [0, 0, 1]])\n",
    "    shear_matrix    = np.array([[1, np.sin(shear), 0], [0, np.cos(shear), 0], [0, 0, 1]])\n",
    "    zoom_matrix     = np.array([[1.0/height_zoom, 0, 0], [0, 1.0/width_zoom, 0], [0, 0, 1]])\n",
    "    shift_matrix    = np.array([[1, 0, -height_shift], [0, 1, -width_shift], [0, 0, 1]])\n",
    "    return np.dot(np.dot(rotation_matrix, shear_matrix), np.dot(zoom_matrix, shift_matrix))\n",
    "\n",
    "def read_cropped_image(p, augment):\n",
    "    \"\"\"\n",
    "    @param p : the name of the picture to read\n",
    "    @param augment: True/False if data augmentation should be performed\n",
    "    @return a numpy array with the transformed image\n",
    "    \"\"\"\n",
    "    # If an image id was given, convert to filename\n",
    "    if p in h2p: p = h2p[p]\n",
    "    size_x,size_y = p2size[p]\n",
    "    \n",
    "    # Determine the region of the original image we want to capture based on the bounding box.\n",
    "    x0,y0,x1,y1   = p2bb[p]\n",
    "    if p in rotate: x0, y0, x1, y1 = size_x - x1, size_y - y1, size_x - x0, size_y - y0\n",
    "    dx            = x1 - x0\n",
    "    dy            = y1 - y0\n",
    "    x0           -= dx*crop_margin\n",
    "    x1           += dx*crop_margin + 1\n",
    "    y0           -= dy*crop_margin\n",
    "    y1           += dy*crop_margin + 1\n",
    "    if (x0 < 0     ): x0 = 0\n",
    "    if (x1 > size_x): x1 = size_x\n",
    "    if (y0 < 0     ): y0 = 0\n",
    "    if (y1 > size_y): y1 = size_y\n",
    "    dx            = x1 - x0\n",
    "    dy            = y1 - y0\n",
    "    if dx > dy*anisotropy:\n",
    "        dy  = 0.5*(dx/anisotropy - dy)\n",
    "        y0 -= dy\n",
    "        y1 += dy\n",
    "    else:\n",
    "        dx  = 0.5*(dy*anisotropy - dx)\n",
    "        x0 -= dx\n",
    "        x1 += dx\n",
    "\n",
    "    # Generate the transformation matrix\n",
    "    trans = np.array([[1, 0, -0.5*img_shape[0]], [0, 1, -0.5*img_shape[1]], [0, 0, 1]])\n",
    "    trans = np.dot(np.array([[(y1 - y0)/img_shape[0], 0, 0], [0, (x1 - x0)/img_shape[1], 0], [0, 0, 1]]), trans)\n",
    "    if augment:\n",
    "        trans = np.dot(build_transform(\n",
    "            random.uniform(-5, 5),\n",
    "            random.uniform(-5, 5),\n",
    "            random.uniform(0.8, 1.0),\n",
    "            random.uniform(0.8, 1.0),\n",
    "            random.uniform(-0.05*(y1 - y0), 0.05*(y1 - y0)),\n",
    "            random.uniform(-0.05*(x1 - x0), 0.05*(x1 - x0))\n",
    "            ), trans)\n",
    "    trans = np.dot(np.array([[1, 0, 0.5*(y1 + y0)], [0, 1, 0.5*(x1 + x0)], [0, 0, 1]]), trans)\n",
    "\n",
    "    # Read the image, transform to black and white and comvert to numpy array\n",
    "    img   = read_raw_image(p).convert('L')\n",
    "    img   = img_to_array(img)\n",
    "    \n",
    "    # Apply affine transformation\n",
    "    matrix = trans[:2,:2]\n",
    "    offset = trans[:2,2]\n",
    "    img    = img.reshape(img.shape[:-1])\n",
    "    img    = affine_transform(img, matrix, offset, output_shape=img_shape[:-1], order=1, mode='constant', cval=np.average(img))\n",
    "    img    = img.reshape(img_shape)\n",
    "\n",
    "    # Normalize to zero mean and unit variance\n",
    "    img  -= np.mean(img, keepdims=True)\n",
    "    img  /= np.std(img, keepdims=True) + K.epsilon()\n",
    "    return img\n",
    "\n",
    "def read_for_training(p):\n",
    "    \"\"\"\n",
    "    Read and preprocess an image with data augmentation (random transform).\n",
    "    \"\"\"\n",
    "    return read_cropped_image(p, True)\n",
    "\n",
    "def read_for_validation(p):\n",
    "    \"\"\"\n",
    "    Read and preprocess an image without data augmentation (use for testing).\n",
    "    \"\"\"\n",
    "    return read_cropped_image(p, False)\n",
    "# p = list(tagged.keys())[312]\n",
    "# imgs = [\n",
    "#     read_raw_image(p),\n",
    "#     array_to_img(read_for_validation(p)),\n",
    "#     array_to_img(read_for_training(p))\n",
    "# ]\n",
    "# show_whale(imgs, per_row=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "e671c843a32506b635f92e75a7574346fe339254"
   },
   "source": [
    "The left image is the original picture. The center image does the test transformation. The right image adds a random data augmentation transformation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "4c35fc67c23a1677d37fe4766cd3fb3e14b840d9"
   },
   "source": [
    "# Siamese Neural Network architecture\n",
    "The Siamese Neural Network compares two images and decides if the two images are taken from the same whale, or different whales. By testing each image from the test set against every image from the training set, the most likely whales can be identified by sorting the pictures in likelihood for a match.\n",
    "\n",
    "A Siamese Neural Network is composed of two parts. A Convolutation Neural Network (CNN) transforms an input image into a vector of features describing the whale. The same CNN, with the same weights, is used for both images. I call the CNN the __branch model__. I used a custom model mostly inspired from [ResNet](https://www.cv-foundation.org/openaccess/content_cvpr_2016/papers/He_Deep_Residual_Learning_CVPR_2016_paper.pdf).\n",
    "\n",
    "What I call the __head model__ is used to compare the feature vectors from the CNN and decide if the whales match or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "f7fbe3e42597f837363c6cdebb5a8853d5c68399"
   },
   "source": [
    "## Head model\n",
    "The head model compares the feature vector from the branch model to decide if the pictures show the same or different whales.  The typical approach is to use a distance measure (e.g. \\\\( L_1 \\\\)) with a contrastive loss function, but there are a few reasons to try something different:\n",
    "\n",
    "* A distance measure will consider two features with value zero as a perfect match, while two features with large but slightly different values will be seen as good, but not quite as good since they are not exactly equal. Still, I feel there is more postive signal in the active features than in the negative ones, especially with ReLU (Rectified Linear Unit) activation, concept that is lost by the distance measure.\n",
    "* Also, a distance measure does not provide for features to be negatively correlated. Consider a case where, if both images have feature X, they must be the same whale, unless they also both have feature Y, in which case X is not as clear.\n",
    "* At the same time, there is this implicit assumption that swapping the two images must produce the same result: if A is the same whale as B, B must be the same whale as A. \n",
    "\n",
    "To address these concerns, i proceed as follow:\n",
    "\n",
    "1. For each feature I compute the sum, the product, the absolute difference and the difference squared (\\\\(x + y, x y, |x - y|, (x - y)^2\\\\)).\n",
    "1. The four values are passed through a small neural network, which can learn how to weigh between matching zeros and close non-zero values. The same neural net with the same weights is used for each feature.\n",
    "1. The output is a weighted sum of the converted features, with a sigmoid activation. The value of the weight is redundant, because the weight is just a scaling factor of the feature, which could be learned by another layer, however, it allows for negative weights, which cannot be produced otherwise when using ReLU activation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "0ba49143534d8fee4d9d7adeee7a912736c91e45"
   },
   "source": [
    "## Branch model\n",
    "The branch model is a regular CNN model. Here are the key elements of its design:\n",
    "* Because the training dataset is small, I tried to keep the number of learned parameter relatively small, while keeping the model expressive enough. A ResNet like architecture is more economical than a [VGG](https://arxiv.org/abs/1409.1556) like network for example.\n",
    "* The problem turned out to be memory bound, with most memory being taken to store activations from the feedforward pass, used to compute the gradient during backpropagation. With Windows 10 and a GTX 1080 card, there is about 6.8GB VRAM available, and this limit has constrained the model choice. \n",
    "\n",
    "The branch model is composed of 6 blocks, each block processing maps with smaller and smaller resolution,, with intermediate pooling layers.\n",
    "\n",
    "* Block 1 - 384x384\n",
    "* Block 2 - 96x96\n",
    "* Block 3 - 48x48\n",
    "* Block 4 - 24x24\n",
    "* Block 5 - 12x12\n",
    "* Block 6 - 6x6\n",
    "\n",
    "Block 1 has a single convolution layer with stride 2 followed by 2x2 max pooling. Because of the high resolution, it uses a lot of memory, so a minimum of work is done here to save memory for subsequent blocks.\n",
    "\n",
    "Block 2 has two 3x3 convolutions similar to VGG. These convolutions are less memory intensive then the subsequent ResNet blocks, and are used to save memory. Note that after this, the tensor has dimension 96x96x64, the same volume as the initial 384x384x1 image, thus we can assume no significant information has been lost.\n",
    "\n",
    "Blocks 3 to 6 perform ResNet like convolution. I suggest reading the original paper, but the idea is to form a subblock with a 1x1 convolution reducing the number of features, a 3x3 convolution and another 1x1 convolution to restore the number of features to the original. The output of these convolutions is then added to the original tensor (bypass connection). I use 4 such subblocks by block, plus a single 1x1 convolution to increase the feature count after each pooling layer.\n",
    "\n",
    "The final step of the branch model is a global max pooling, which makes the model robust to fluke not being always well centered."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "9f11d92483b771073da6c1a2f0eb62815e92aee0"
   },
   "source": [
    "## Code\n",
    "The following is the Keras code for the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "_uuid": "7f90b807f6236490117f67b485bb2268b32eb8b4",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_12 (InputLayer)           (None, 512)          0                                            \n",
      "__________________________________________________________________________________________________\n",
      "input_13 (InputLayer)           (None, 512)          0                                            \n",
      "__________________________________________________________________________________________________\n",
      "lambda_11 (Lambda)              (None, 512)          0           input_12[0][0]                   \n",
      "                                                                 input_13[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "lambda_9 (Lambda)               (None, 512)          0           input_12[0][0]                   \n",
      "                                                                 input_13[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "lambda_10 (Lambda)              (None, 512)          0           input_12[0][0]                   \n",
      "                                                                 input_13[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "lambda_12 (Lambda)              (None, 512)          0           lambda_11[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_3 (Concatenate)     (None, 2048)         0           lambda_9[0][0]                   \n",
      "                                                                 lambda_10[0][0]                  \n",
      "                                                                 lambda_11[0][0]                  \n",
      "                                                                 lambda_12[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "reshape1 (Reshape)              (None, 4, 512, 1)    0           concatenate_3[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_170 (Conv2D)             (None, 1, 512, 32)   160         reshape1[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "reshape_3 (Reshape)             (None, 512, 32, 1)   0           conv2d_170[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_171 (Conv2D)             (None, 512, 1, 1)    33          reshape_3[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "flatten (Flatten)               (None, 512)          0           conv2d_171[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "weighted-average (Dense)        (None, 1)            513         flatten[0][0]                    \n",
      "==================================================================================================\n",
      "Total params: 706\n",
      "Trainable params: 706\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "from keras import regularizers\n",
    "from keras.optimizers import Adam\n",
    "from keras.engine.topology import Input\n",
    "from keras.layers import Activation, Add, BatchNormalization, Concatenate, Conv2D, Dense, Flatten, GlobalMaxPooling2D, Lambda, MaxPooling2D, Reshape\n",
    "from keras.models import Model\n",
    "\n",
    "def subblock(x, filter, **kwargs):\n",
    "    x = BatchNormalization()(x)\n",
    "    y = x\n",
    "    y = Conv2D(filter, (1, 1), activation='relu', **kwargs)(y) # Reduce the number of features to 'filter'\n",
    "    y = BatchNormalization()(y)\n",
    "    y = Conv2D(filter, (3, 3), activation='relu', **kwargs)(y) # Extend the feature field\n",
    "    y = BatchNormalization()(y)\n",
    "    y = Conv2D(K.int_shape(x)[-1], (1, 1), **kwargs)(y) # no activation # Restore the number of original features\n",
    "    y = Add()([x,y]) # Add the bypass connection\n",
    "    y = Activation('relu')(y)\n",
    "    return y\n",
    "\n",
    "def build_model(lr, l2, activation='sigmoid'):\n",
    "\n",
    "    ##############\n",
    "    # BRANCH MODEL\n",
    "    ##############\n",
    "    regul  = regularizers.l2(l2)\n",
    "    optim  = Adam(lr=lr)\n",
    "    kwargs = {'padding':'same', 'kernel_regularizer':regul}\n",
    "\n",
    "    inp = Input(shape=img_shape) # 384x384x1\n",
    "    x   = Conv2D(64, (9,9), strides=2, activation='relu', **kwargs)(inp)\n",
    "\n",
    "    x   = MaxPooling2D((2, 2), strides=(2, 2))(x) # 96x96x64\n",
    "    for _ in range(2):\n",
    "        x = BatchNormalization()(x)\n",
    "        x = Conv2D(64, (3,3), activation='relu', **kwargs)(x)\n",
    "\n",
    "    x = MaxPooling2D((2, 2), strides=(2, 2))(x) # 48x48x64\n",
    "    x = BatchNormalization()(x)\n",
    "    x = Conv2D(128, (1,1), activation='relu', **kwargs)(x) # 48x48x128\n",
    "    for _ in range(4): x = subblock(x, 64, **kwargs)\n",
    "\n",
    "    x = MaxPooling2D((2, 2), strides=(2, 2))(x) # 24x24x128\n",
    "    x = BatchNormalization()(x)\n",
    "    x = Conv2D(256, (1,1), activation='relu', **kwargs)(x) # 24x24x256\n",
    "    for _ in range(4): x = subblock(x, 64, **kwargs)\n",
    "\n",
    "    x = MaxPooling2D((2, 2), strides=(2, 2))(x) # 12x12x256\n",
    "    x = BatchNormalization()(x)\n",
    "    x = Conv2D(384, (1,1), activation='relu', **kwargs)(x) # 12x12x384\n",
    "    for _ in range(4): x = subblock(x, 96, **kwargs)\n",
    "\n",
    "    x = MaxPooling2D((2, 2), strides=(2, 2))(x) # 6x6x384\n",
    "    x = BatchNormalization()(x)\n",
    "    x = Conv2D(512, (1,1), activation='relu', **kwargs)(x) # 6x6x512\n",
    "    for _ in range(4): x = subblock(x, 128, **kwargs)\n",
    "    \n",
    "    x             = GlobalMaxPooling2D()(x) # 512\n",
    "    branch_model  = Model(inp, x)\n",
    "    \n",
    "    ############\n",
    "    # HEAD MODEL\n",
    "    ############\n",
    "    mid        = 32\n",
    "    xa_inp     = Input(shape=branch_model.output_shape[1:])\n",
    "    xb_inp     = Input(shape=branch_model.output_shape[1:])\n",
    "    x1         = Lambda(lambda x : x[0]*x[1])([xa_inp, xb_inp])\n",
    "    x2         = Lambda(lambda x : x[0] + x[1])([xa_inp, xb_inp])\n",
    "    x3         = Lambda(lambda x : K.abs(x[0] - x[1]))([xa_inp, xb_inp])\n",
    "    x4         = Lambda(lambda x : K.square(x))(x3)\n",
    "    x          = Concatenate()([x1, x2, x3, x4])\n",
    "    x          = Reshape((4, branch_model.output_shape[1], 1), name='reshape1')(x)\n",
    "\n",
    "    # Per feature NN with shared weight is implemented using CONV2D with appropriate stride.\n",
    "    x          = Conv2D(mid, (4, 1), activation='relu', padding='valid')(x)\n",
    "    x          = Reshape((branch_model.output_shape[1], mid, 1))(x)\n",
    "    x          = Conv2D(1, (1, mid), activation='linear', padding='valid')(x)\n",
    "    x          = Flatten(name='flatten')(x)\n",
    "    \n",
    "    # Weighted sum implemented as a Dense layer.\n",
    "    x          = Dense(1, use_bias=True, activation=activation, name='weighted-average')(x)\n",
    "    head_model = Model([xa_inp, xb_inp], x, name='head')\n",
    "\n",
    "    ########################\n",
    "    # SIAMESE NEURAL NETWORK\n",
    "    ########################\n",
    "    # Complete model is constructed by calling the branch model on each input image,\n",
    "    # and then the head model on the resulting 512-vectors.\n",
    "    img_a      = Input(shape=img_shape)\n",
    "    img_b      = Input(shape=img_shape)\n",
    "    xa         = branch_model(img_a)\n",
    "    xb         = branch_model(img_b)\n",
    "    x          = head_model([xa, xb])\n",
    "    model      = Model([img_a, img_b], x)\n",
    "    model.compile(optim, loss='binary_crossentropy', metrics=['binary_crossentropy', 'acc'])\n",
    "    return model, branch_model, head_model\n",
    "\n",
    "model, branch_model, head_model = build_model(64e-5,0)\n",
    "head_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "_uuid": "7784a525bdebc38d2b89b106f33c2893e58ebf04"
   },
   "outputs": [],
   "source": [
    "# from keras.utils import plot_model\n",
    "# plot_model(head_model, to_file='head-model.png')\n",
    "# pil_image.open('head-model.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "_uuid": "0efc2c31c8d2a1aa8aef5d135439d8cf2e5489d2"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_11 (InputLayer)           (None, 384, 384, 1)  0                                            \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_115 (Conv2D)             (None, 192, 192, 64) 5248        input_11[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling2d_11 (MaxPooling2D) (None, 96, 96, 64)   0           conv2d_115[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_109 (BatchN (None, 96, 96, 64)   256         max_pooling2d_11[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_116 (Conv2D)             (None, 96, 96, 64)   36928       batch_normalization_109[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_110 (BatchN (None, 96, 96, 64)   256         conv2d_116[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_117 (Conv2D)             (None, 96, 96, 64)   36928       batch_normalization_110[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling2d_12 (MaxPooling2D) (None, 48, 48, 64)   0           conv2d_117[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_111 (BatchN (None, 48, 48, 64)   256         max_pooling2d_12[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_118 (Conv2D)             (None, 48, 48, 128)  8320        batch_normalization_111[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_112 (BatchN (None, 48, 48, 128)  512         conv2d_118[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_119 (Conv2D)             (None, 48, 48, 64)   8256        batch_normalization_112[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_113 (BatchN (None, 48, 48, 64)   256         conv2d_119[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_120 (Conv2D)             (None, 48, 48, 64)   36928       batch_normalization_113[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_114 (BatchN (None, 48, 48, 64)   256         conv2d_120[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_121 (Conv2D)             (None, 48, 48, 128)  8320        batch_normalization_114[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_33 (Add)                    (None, 48, 48, 128)  0           batch_normalization_112[0][0]    \n",
      "                                                                 conv2d_121[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_33 (Activation)      (None, 48, 48, 128)  0           add_33[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_115 (BatchN (None, 48, 48, 128)  512         activation_33[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_122 (Conv2D)             (None, 48, 48, 64)   8256        batch_normalization_115[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_116 (BatchN (None, 48, 48, 64)   256         conv2d_122[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_123 (Conv2D)             (None, 48, 48, 64)   36928       batch_normalization_116[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_117 (BatchN (None, 48, 48, 64)   256         conv2d_123[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_124 (Conv2D)             (None, 48, 48, 128)  8320        batch_normalization_117[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_34 (Add)                    (None, 48, 48, 128)  0           batch_normalization_115[0][0]    \n",
      "                                                                 conv2d_124[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_34 (Activation)      (None, 48, 48, 128)  0           add_34[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_118 (BatchN (None, 48, 48, 128)  512         activation_34[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_125 (Conv2D)             (None, 48, 48, 64)   8256        batch_normalization_118[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_119 (BatchN (None, 48, 48, 64)   256         conv2d_125[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_126 (Conv2D)             (None, 48, 48, 64)   36928       batch_normalization_119[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_120 (BatchN (None, 48, 48, 64)   256         conv2d_126[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_127 (Conv2D)             (None, 48, 48, 128)  8320        batch_normalization_120[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_35 (Add)                    (None, 48, 48, 128)  0           batch_normalization_118[0][0]    \n",
      "                                                                 conv2d_127[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_35 (Activation)      (None, 48, 48, 128)  0           add_35[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_121 (BatchN (None, 48, 48, 128)  512         activation_35[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_128 (Conv2D)             (None, 48, 48, 64)   8256        batch_normalization_121[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_122 (BatchN (None, 48, 48, 64)   256         conv2d_128[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_129 (Conv2D)             (None, 48, 48, 64)   36928       batch_normalization_122[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_123 (BatchN (None, 48, 48, 64)   256         conv2d_129[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_130 (Conv2D)             (None, 48, 48, 128)  8320        batch_normalization_123[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_36 (Add)                    (None, 48, 48, 128)  0           batch_normalization_121[0][0]    \n",
      "                                                                 conv2d_130[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_36 (Activation)      (None, 48, 48, 128)  0           add_36[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling2d_13 (MaxPooling2D) (None, 24, 24, 128)  0           activation_36[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_124 (BatchN (None, 24, 24, 128)  512         max_pooling2d_13[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_131 (Conv2D)             (None, 24, 24, 256)  33024       batch_normalization_124[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_125 (BatchN (None, 24, 24, 256)  1024        conv2d_131[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_132 (Conv2D)             (None, 24, 24, 64)   16448       batch_normalization_125[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_126 (BatchN (None, 24, 24, 64)   256         conv2d_132[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_133 (Conv2D)             (None, 24, 24, 64)   36928       batch_normalization_126[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_127 (BatchN (None, 24, 24, 64)   256         conv2d_133[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_134 (Conv2D)             (None, 24, 24, 256)  16640       batch_normalization_127[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_37 (Add)                    (None, 24, 24, 256)  0           batch_normalization_125[0][0]    \n",
      "                                                                 conv2d_134[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_37 (Activation)      (None, 24, 24, 256)  0           add_37[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_128 (BatchN (None, 24, 24, 256)  1024        activation_37[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_135 (Conv2D)             (None, 24, 24, 64)   16448       batch_normalization_128[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_129 (BatchN (None, 24, 24, 64)   256         conv2d_135[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_136 (Conv2D)             (None, 24, 24, 64)   36928       batch_normalization_129[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_130 (BatchN (None, 24, 24, 64)   256         conv2d_136[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_137 (Conv2D)             (None, 24, 24, 256)  16640       batch_normalization_130[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_38 (Add)                    (None, 24, 24, 256)  0           batch_normalization_128[0][0]    \n",
      "                                                                 conv2d_137[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_38 (Activation)      (None, 24, 24, 256)  0           add_38[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_131 (BatchN (None, 24, 24, 256)  1024        activation_38[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_138 (Conv2D)             (None, 24, 24, 64)   16448       batch_normalization_131[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_132 (BatchN (None, 24, 24, 64)   256         conv2d_138[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_139 (Conv2D)             (None, 24, 24, 64)   36928       batch_normalization_132[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_133 (BatchN (None, 24, 24, 64)   256         conv2d_139[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_140 (Conv2D)             (None, 24, 24, 256)  16640       batch_normalization_133[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_39 (Add)                    (None, 24, 24, 256)  0           batch_normalization_131[0][0]    \n",
      "                                                                 conv2d_140[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_39 (Activation)      (None, 24, 24, 256)  0           add_39[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_134 (BatchN (None, 24, 24, 256)  1024        activation_39[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_141 (Conv2D)             (None, 24, 24, 64)   16448       batch_normalization_134[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_135 (BatchN (None, 24, 24, 64)   256         conv2d_141[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_142 (Conv2D)             (None, 24, 24, 64)   36928       batch_normalization_135[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_136 (BatchN (None, 24, 24, 64)   256         conv2d_142[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_143 (Conv2D)             (None, 24, 24, 256)  16640       batch_normalization_136[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_40 (Add)                    (None, 24, 24, 256)  0           batch_normalization_134[0][0]    \n",
      "                                                                 conv2d_143[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_40 (Activation)      (None, 24, 24, 256)  0           add_40[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling2d_14 (MaxPooling2D) (None, 12, 12, 256)  0           activation_40[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_137 (BatchN (None, 12, 12, 256)  1024        max_pooling2d_14[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_144 (Conv2D)             (None, 12, 12, 384)  98688       batch_normalization_137[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_138 (BatchN (None, 12, 12, 384)  1536        conv2d_144[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_145 (Conv2D)             (None, 12, 12, 96)   36960       batch_normalization_138[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_139 (BatchN (None, 12, 12, 96)   384         conv2d_145[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_146 (Conv2D)             (None, 12, 12, 96)   83040       batch_normalization_139[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_140 (BatchN (None, 12, 12, 96)   384         conv2d_146[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_147 (Conv2D)             (None, 12, 12, 384)  37248       batch_normalization_140[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_41 (Add)                    (None, 12, 12, 384)  0           batch_normalization_138[0][0]    \n",
      "                                                                 conv2d_147[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_41 (Activation)      (None, 12, 12, 384)  0           add_41[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_141 (BatchN (None, 12, 12, 384)  1536        activation_41[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_148 (Conv2D)             (None, 12, 12, 96)   36960       batch_normalization_141[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_142 (BatchN (None, 12, 12, 96)   384         conv2d_148[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_149 (Conv2D)             (None, 12, 12, 96)   83040       batch_normalization_142[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_143 (BatchN (None, 12, 12, 96)   384         conv2d_149[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_150 (Conv2D)             (None, 12, 12, 384)  37248       batch_normalization_143[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_42 (Add)                    (None, 12, 12, 384)  0           batch_normalization_141[0][0]    \n",
      "                                                                 conv2d_150[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_42 (Activation)      (None, 12, 12, 384)  0           add_42[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_144 (BatchN (None, 12, 12, 384)  1536        activation_42[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_151 (Conv2D)             (None, 12, 12, 96)   36960       batch_normalization_144[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_145 (BatchN (None, 12, 12, 96)   384         conv2d_151[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_152 (Conv2D)             (None, 12, 12, 96)   83040       batch_normalization_145[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_146 (BatchN (None, 12, 12, 96)   384         conv2d_152[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_153 (Conv2D)             (None, 12, 12, 384)  37248       batch_normalization_146[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_43 (Add)                    (None, 12, 12, 384)  0           batch_normalization_144[0][0]    \n",
      "                                                                 conv2d_153[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_43 (Activation)      (None, 12, 12, 384)  0           add_43[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_147 (BatchN (None, 12, 12, 384)  1536        activation_43[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_154 (Conv2D)             (None, 12, 12, 96)   36960       batch_normalization_147[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_148 (BatchN (None, 12, 12, 96)   384         conv2d_154[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_155 (Conv2D)             (None, 12, 12, 96)   83040       batch_normalization_148[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_149 (BatchN (None, 12, 12, 96)   384         conv2d_155[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_156 (Conv2D)             (None, 12, 12, 384)  37248       batch_normalization_149[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_44 (Add)                    (None, 12, 12, 384)  0           batch_normalization_147[0][0]    \n",
      "                                                                 conv2d_156[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_44 (Activation)      (None, 12, 12, 384)  0           add_44[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling2d_15 (MaxPooling2D) (None, 6, 6, 384)    0           activation_44[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_150 (BatchN (None, 6, 6, 384)    1536        max_pooling2d_15[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_157 (Conv2D)             (None, 6, 6, 512)    197120      batch_normalization_150[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_151 (BatchN (None, 6, 6, 512)    2048        conv2d_157[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_158 (Conv2D)             (None, 6, 6, 128)    65664       batch_normalization_151[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_152 (BatchN (None, 6, 6, 128)    512         conv2d_158[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_159 (Conv2D)             (None, 6, 6, 128)    147584      batch_normalization_152[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_153 (BatchN (None, 6, 6, 128)    512         conv2d_159[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_160 (Conv2D)             (None, 6, 6, 512)    66048       batch_normalization_153[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_45 (Add)                    (None, 6, 6, 512)    0           batch_normalization_151[0][0]    \n",
      "                                                                 conv2d_160[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_45 (Activation)      (None, 6, 6, 512)    0           add_45[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_154 (BatchN (None, 6, 6, 512)    2048        activation_45[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_161 (Conv2D)             (None, 6, 6, 128)    65664       batch_normalization_154[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_155 (BatchN (None, 6, 6, 128)    512         conv2d_161[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_162 (Conv2D)             (None, 6, 6, 128)    147584      batch_normalization_155[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_156 (BatchN (None, 6, 6, 128)    512         conv2d_162[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_163 (Conv2D)             (None, 6, 6, 512)    66048       batch_normalization_156[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_46 (Add)                    (None, 6, 6, 512)    0           batch_normalization_154[0][0]    \n",
      "                                                                 conv2d_163[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_46 (Activation)      (None, 6, 6, 512)    0           add_46[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_157 (BatchN (None, 6, 6, 512)    2048        activation_46[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_164 (Conv2D)             (None, 6, 6, 128)    65664       batch_normalization_157[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_158 (BatchN (None, 6, 6, 128)    512         conv2d_164[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_165 (Conv2D)             (None, 6, 6, 128)    147584      batch_normalization_158[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_159 (BatchN (None, 6, 6, 128)    512         conv2d_165[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_166 (Conv2D)             (None, 6, 6, 512)    66048       batch_normalization_159[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_47 (Add)                    (None, 6, 6, 512)    0           batch_normalization_157[0][0]    \n",
      "                                                                 conv2d_166[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_47 (Activation)      (None, 6, 6, 512)    0           add_47[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_160 (BatchN (None, 6, 6, 512)    2048        activation_47[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_167 (Conv2D)             (None, 6, 6, 128)    65664       batch_normalization_160[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_161 (BatchN (None, 6, 6, 128)    512         conv2d_167[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_168 (Conv2D)             (None, 6, 6, 128)    147584      batch_normalization_161[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_162 (BatchN (None, 6, 6, 128)    512         conv2d_168[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_169 (Conv2D)             (None, 6, 6, 512)    66048       batch_normalization_162[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "add_48 (Add)                    (None, 6, 6, 512)    0           batch_normalization_160[0][0]    \n",
      "                                                                 conv2d_169[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "activation_48 (Activation)      (None, 6, 6, 512)    0           add_48[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "global_max_pooling2d_3 (GlobalM (None, 512)          0           activation_48[0][0]              \n",
      "==================================================================================================\n",
      "Total params: 2,692,096\n",
      "Trainable params: 2,674,304\n",
      "Non-trainable params: 17,792\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "branch_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "_uuid": "16cb29a5f2fe1e1923f51ec5207297dfeb69c216"
   },
   "outputs": [],
   "source": [
    "# # Oops, this is HUGE!\n",
    "# plot_model(branch_model, to_file='branch-model.png')\n",
    "# img = pil_image.open('branch-model.png')\n",
    "# img.resize([x//2 for x in img.size])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "1f6cb609b15e172955f634c8f68f15a98d6944b0",
    "collapsed": true
   },
   "source": [
    "# Training data construction\n",
    "As highlighted in the abstract section, this is the part that makes a big difference in the model accuracy.\n",
    "\n",
    "We want the Siamese Neural Network to pick the one correct whale from all the possible whales from the training set. While scoring the correct whale high, it must simultaneously score **all** other whales lower. It is not enough to have a random whale score low. To force all the other whales to a low probability, the training algorithm presents pairs of pictures with increasing difficulty, as evaluated by the model at any given time. Essentially, we forcus the training on pairs that he model is getting wrong, as a form of adversarial training.\n",
    "\n",
    "At the same time, we want the model to recognize **whales** and not **pictures**. Given the small number of pictures in the training dataset, it is not irrealistic to imagine the model recognizing a specific picture using the shape of a wave, or a bird flying by. To prevent this, the data presented to the model must be unbiased. If a picture is used more often in negative examples, the model risks simply learning to guess a mismatch whenever this picture is present, without learning how to compare the whales correctly. By presenting each image an equal number of times, with 50% positive and 50% negative examples, the model has no incentive in learning to recognize specific pictures, and thus focuses on recognizing whales as desired."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "7e04e1c924396cda1834e675b2de8a447267116d"
   },
   "source": [
    "## Image selection\n",
    "To begin, we reduce the number of images from the training set:\n",
    "\n",
    "* Images from the blacklist are removed;\n",
    "* Duplicate images are removed;\n",
    "* All 'new_whale' images are removed;\n",
    "* All whales with a single image are removed.\n",
    "\n",
    "The blacklist was constructed manually by spotting images unhelpful to training. Reasons could be the underside of the fluke is not visible, or we see only dead fluke fragments on the a beach, there are two whales in the picture, etc. The list is in no way exhaustive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "_uuid": "a3403d1d59eb66f496f042c0a08b3286814b2117"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "34"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with open('./exclude.txt', 'rt') as f: exclude = f.read().split('\\n')[:-1]   \n",
    "len(exclude)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "_uuid": "edca49e4b9e795a202c426804104d0ac917905c0"
   },
   "outputs": [],
   "source": [
    "# show_whale([read_raw_image(p) for p in exclude], per_row=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "4edecade09fba068c488c50b86f08a59417cba74",
    "collapsed": true
   },
   "source": [
    "## Matching whale examples\n",
    "Half the examples used during training are for pair of images. For each whale of the training set, compute a [derangement](https://en.wikipedia.org/wiki/Derangement) of its pictures. Use the original order as picture A, and the derangment as picture B. This creates a random number of matching image pairs, with each image taken exactly two times."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "0fd59cc7ab45460467f1df0fc32c75e089d351f8",
    "collapsed": true
   },
   "source": [
    "## Different whales examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "ecfa5542699ab68b70078ef65feaba5d4c0fcac2"
   },
   "source": [
    "Different whales examples are generated by computing a derangement of all pictures from the training set, subject to:\n",
    "\n",
    "* The image pairs must belong to different whales;\n",
    "* The pairs formed must be difficult for the model to distinguish.\n",
    "\n",
    "The following algorith is used to generate the pairs:\n",
    "\n",
    "1. The similarity between each pair of image is computed using the current model state. This has complexity \\\\( \\frac {n(n-1)} 2  \\\\), where n is the size of the training set. Fortunately, only the head model must be computed for all pairs, and it is very fast. The 512-feature vectors can be pre-computed once for each image, i.e. O(n) complexity.\n",
    "1. Entry that correpond to pair of images from the same whale are set with similarity \\\\( -\\infty \\\\).\n",
    "1. [Linear sum assignment](https://en.wikipedia.org/wiki/Assignment_problem) algorithm is used to find the most difficult matching. \n",
    "1. To randomize the selection, and control the matching difficulty, we add a random matrix to the cost matrix from step #1 . The random matrix has values uniformly distributed between 0 and K. The larger K, the more random the matching. The smaller K, the more difficult the pairing is for the model.\n",
    "1. To produce different matching for successive epochs, the selected entries in the matrix are overwritten with \\\\( -\\infty \\\\) to force an alternate selection for the next matching."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "70ba1841da5b8d6d51ce06734d4405ef4c129d6f"
   },
   "source": [
    "## Code\n",
    "The described logic is essentially implemented by the TrainingData class, that performs just in time data augmentation as well as computing the matching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "_uuid": "4b0b2120d9f3d38bcc0971e0c0baab3e07dcbbcc"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "605"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Find all the whales associated with an image id. It can be ambiguous as duplicate images may have different whale ids.\n",
    "h2ws = {}\n",
    "new_whale = 'new_whale'\n",
    "for p,w in tagged.items():\n",
    "    if w != new_whale: # Use only identified whales\n",
    "        h = p2h[p]\n",
    "        if h not in h2ws: h2ws[h] = []\n",
    "        if w not in h2ws[h]: h2ws[h].append(w)\n",
    "for h,ws in h2ws.items():\n",
    "    if len(ws) > 1:\n",
    "        h2ws[h] = sorted(ws)\n",
    "len(h2ws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "_uuid": "8844f37d5d2c40eb9347dfb29b2383ecdf417557"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f43183739a8f53c4\n",
      "9467679c9b638c98\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "561"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For each whale, find the unambiguous images ids.\n",
    "w2hs = {}\n",
    "for h,ws in h2ws.items():\n",
    "    if len(ws) == 1: # Use only unambiguous pictures\n",
    "        if h2p[h] in exclude:\n",
    "            print(h) # Skip excluded images\n",
    "        else:\n",
    "            w = ws[0]\n",
    "            if w not in w2hs: w2hs[w] = []\n",
    "            if h not in w2hs[w]: w2hs[w].append(h)\n",
    "for w,hs in w2hs.items():\n",
    "    if len(hs) > 1:\n",
    "        w2hs[w] = sorted(hs)\n",
    "len(w2hs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "_uuid": "e1420dd75317af6d7876f86187713e28e4cd385d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(79, 37)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Find the list of training images, keep only whales with at least two images.\n",
    "train = [] # A list of training image ids\n",
    "for hs in w2hs.values():\n",
    "    if len(hs) > 1:\n",
    "        train += hs\n",
    "random.shuffle(train)\n",
    "train_set = set(train)\n",
    "\n",
    "w2ts = {} # Associate the image ids from train to each whale id.\n",
    "for w,hs in w2hs.items():\n",
    "    for h in hs:\n",
    "        if h in train_set:\n",
    "            if w not in w2ts: w2ts[w] = []\n",
    "            if h not in w2ts[w]: w2ts[w].append(h)\n",
    "for w,ts in w2ts.items(): w2ts[w] = np.array(ts)\n",
    "    \n",
    "t2i = {} # The position in train of each training image id\n",
    "for i,t in enumerate(train): t2i[t] = i\n",
    "\n",
    "len(train),len(w2ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "_uuid": "bce8abd37c4ccdcc0059504c03e2d4f3371fd170"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Module lap not found, emulating with much slower scipy.optimize.linear_sum_assignment\n"
     ]
    }
   ],
   "source": [
    "from keras.utils import Sequence\n",
    "\n",
    "# First try to use lapjv Linear Assignment Problem solver as it is much faster.\n",
    "# At the time I am writing this, kaggle kernel with custom package fail to commit.\n",
    "# scipy can be used as a fallback, but it is too slow to run this kernel under the time limit\n",
    "# As a workaround, use scipy with data partitioning.\n",
    "# Because algorithm is O(n^3), small partitions are much faster, but not what produced the submitted solution\n",
    "try:\n",
    "    from lap import lapjv\n",
    "    segment = False\n",
    "except ImportError:\n",
    "    print('Module lap not found, emulating with much slower scipy.optimize.linear_sum_assignment')\n",
    "    segment = True\n",
    "    from scipy.optimize import linear_sum_assignment\n",
    "\n",
    "class TrainingData(Sequence):\n",
    "    def __init__(self, score, steps=1000, batch_size=1):\n",
    "        \"\"\"\n",
    "        @param score the cost matrix for the picture matching\n",
    "        @param steps the number of epoch we are planning with this score matrix\n",
    "        \"\"\"\n",
    "        super(TrainingData, self).__init__()\n",
    "        self.score      = -score # Maximizing the score is the same as minimuzing -score.\n",
    "        self.steps      = steps\n",
    "        self.batch_size = batch_size\n",
    "        for ts in w2ts.values():\n",
    "            idxs = [t2i[t] for t in ts]\n",
    "            for i in idxs:\n",
    "                for j in idxs:\n",
    "                    self.score[i,j] = 10000.0 # Set a large value for matching whales -- eliminates this potential pairing\n",
    "        self.on_epoch_end()\n",
    "    def __getitem__(self, index):\n",
    "        start = self.batch_size*index\n",
    "        end   = min(start + self.batch_size, len(self.match) + len(self.unmatch))\n",
    "        size  = end - start\n",
    "        assert size > 0\n",
    "        a     = np.zeros((size,) + img_shape, dtype=K.floatx())\n",
    "        b     = np.zeros((size,) + img_shape, dtype=K.floatx())\n",
    "        c     = np.zeros((size,1), dtype=K.floatx())\n",
    "        j     = start//2\n",
    "        for i in range(0, size, 2):\n",
    "            a[i,  :,:,:] = read_for_training(self.match[j][0])\n",
    "            b[i,  :,:,:] = read_for_training(self.match[j][1])\n",
    "            c[i,  0    ] = 1 # This is a match\n",
    "            a[i+1,:,:,:] = read_for_training(self.unmatch[j][0])\n",
    "            b[i+1,:,:,:] = read_for_training(self.unmatch[j][1])\n",
    "            c[i+1,0    ] = 0 # Different whales\n",
    "            j           += 1\n",
    "        return [a,b],c\n",
    "    def on_epoch_end(self):\n",
    "        if self.steps <= 0: return # Skip this on the last epoch.\n",
    "        self.steps     -= 1\n",
    "        self.match      = []\n",
    "        self.unmatch    = []\n",
    "        if segment:\n",
    "            # Using slow scipy. Make small batches.\n",
    "            # Because algorithm is O(n^3), small batches are much faster.\n",
    "            # However, this does not find the real optimum, just an approximation.\n",
    "            tmp   = []\n",
    "            batch = 512\n",
    "            for start in range(0, score.shape[0], batch):\n",
    "                end = min(score.shape[0], start + batch)\n",
    "                _, x = linear_sum_assignment(self.score[start:end, start:end])\n",
    "                tmp.append(x + start)\n",
    "            x = np.concatenate(tmp)\n",
    "        else:\n",
    "            print('1',time.time())\n",
    "            _,_,x = lapjv(self.score) # Solve the linear assignment problem\n",
    "            print(time.time())\n",
    "        y = np.arange(len(x),dtype=np.int32)\n",
    "\n",
    "        # Compute a derangement for matching whales\n",
    "        for ts in w2ts.values():\n",
    "            d = ts.copy()\n",
    "            while True:\n",
    "                random.shuffle(d)\n",
    "                if not np.any(ts == d): break\n",
    "            for ab in zip(ts,d): self.match.append(ab)\n",
    "\n",
    "        # Construct unmatched whale pairs from the LAP solution.\n",
    "        for i,j in zip(x,y):\n",
    "            if i == j:\n",
    "                print(self.score)\n",
    "                print(x)\n",
    "                print(y)\n",
    "                print(i,j)\n",
    "            assert i != j\n",
    "            self.unmatch.append((train[i],train[j]))\n",
    "\n",
    "        # Force a different choice for an eventual next epoch.\n",
    "        self.score[x,y] = 10000.0\n",
    "        self.score[y,x] = 10000.0\n",
    "        random.shuffle(self.match)\n",
    "        random.shuffle(self.unmatch)\n",
    "        # print(len(self.match), len(train), len(self.unmatch), len(train))\n",
    "        assert len(self.match) == len(train) and len(self.unmatch) == len(train)\n",
    "    def __len__(self):\n",
    "        return (len(self.match) + len(self.unmatch) + self.batch_size - 1)//self.batch_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "_uuid": "38a2392b38554486c1396ecf9614c3785c03c08c"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<__main__.TrainingData object at 0x7ff7e5b62160>\n"
     ]
    }
   ],
   "source": [
    "# Test on a batch of 32 with random costs.\n",
    "score = np.random.random_sample(size=(len(train),len(train)))\n",
    "data = TrainingData(score)\n",
    "print(data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "index 1 is out of bounds for axis 0 with size 1",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-61-092e5475a53d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-59-c9d0da2a885e>\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, index)\u001b[0m\n\u001b[1;32m     43\u001b[0m             \u001b[0mb\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m  \u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mread_for_training\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatch\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m             \u001b[0mc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m  \u001b[0;36m0\u001b[0m    \u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1\u001b[0m \u001b[0;31m# This is a match\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 45\u001b[0;31m             \u001b[0ma\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mread_for_training\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munmatch\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     46\u001b[0m             \u001b[0mb\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mread_for_training\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munmatch\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m             \u001b[0mc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m    \u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m \u001b[0;31m# Different whales\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: index 1 is out of bounds for axis 0 with size 1"
     ]
    }
   ],
   "source": [
    "(a, b), c = data[0]\n",
    "a.shape, b.shape, c.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "54c10e3e139a93c59a125a50329d59ef6059083b"
   },
   "outputs": [],
   "source": [
    "# First pair is for matching whale\n",
    "imgs = [array_to_img(a[0]), array_to_img(b[0])]\n",
    "show_whale(imgs, per_row=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "63da9a872cbea6641fb7a8740696b5b559fa38f9"
   },
   "outputs": [],
   "source": [
    "# Second pair is for different whales\n",
    "imgs = [array_to_img(a[1]), array_to_img(b[1])]\n",
    "show_whale(imgs, per_row=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "6571f015c52491684f8af4cf88bdaebbd09b018c",
    "collapsed": true
   },
   "source": [
    "# Training procedure\n",
    "This scection describes the procedure used to train the model. Training lasts 400 epochs, with the following quantities changing as the training progresses:\n",
    "\n",
    "* Learning rate\n",
    "* L2 regularization\n",
    "* Constant K measuring the scale of the random component of the score matrix used to match similar images to construct tough training cases.\n",
    "\n",
    "The procedure itself has evolved from many previous experiments, trial and errors on earlier versions of the model.\n",
    "\n",
    "Training the large model from random weights is difficult. In fact, if the model is initially fed examples that are too hard, it does not converge at all. In the context here, hard examples are similar images belonging to different whales. Pushed to the extreme, it it possible to construct a training dataset for which pairs of pictures of different whales appear (to the model) more similar than pairs of pictures from the same whale, making the model learn to classify similar images as different whales, and dissimilar images as same whales!\n",
    "\n",
    "To prevent this, early training is executed with a large value of K, making the negative examples essentially random pairs of pictures of different whales. As the model ability to distinguish between whales increases, K is gradually reduced, presenting harder training cases. Similarly, training starts with no L2 regularization. After 250 epochs, trainings accuracy is fantastic, but it also grossly overfits. At this point, L2 regularization is applied, learning rate is reset to a large value and training continues for an additional 150 epochs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "0a3e377c8cb687a83799c0f1a123b7216c08495c"
   },
   "source": [
    "The following table shows the exact schedule used for the learning rate (LR), L2 regularization (L2) and randomized score matrix (K). \n",
    "\n",
    "Also note that the score matrix for the Linear Assignment Problem is computed at every 5 epochs starting with epoch 10.\n",
    "\n",
    "Epochs | LR | K | L2\n",
    "---- | ---: | ---: | ---:\n",
    "1-10 | \\\\( 64.10^{-5} \\\\)  | \\\\( +\\infty \\\\) | 0\n",
    "11-15 | \\\\( 64.10^{-5} \\\\) | 100.00 | 0\n",
    "16-20 | \\\\( 64.10^{-5} \\\\) | 63.10 | 0\n",
    "21-25 | \\\\( 64.10^{-5} \\\\) | 39.81 | 0\n",
    "26-30 | \\\\( 64.10^{-5} \\\\) | 25.12 | 0\n",
    "31-35 | \\\\( 64.10^{-5} \\\\) | 15.85 | 0\n",
    "36-40 | \\\\( 64.10^{-5} \\\\) | 10.00 | 0\n",
    "41-45 | \\\\( 64.10^{-5} \\\\) | 6.31 | 0\n",
    "46-50 | \\\\( 64.10^{-5} \\\\) | 3.98 | 0\n",
    "51-55 | \\\\( 64.10^{-5} \\\\) | 2.51 | 0\n",
    "56-60 | \\\\( 64.10^{-5} \\\\) | 1.58 | 0\n",
    "61-150 | \\\\( 64.10^{-5} \\\\) | 1.00 | 0\n",
    "151-200 | \\\\( 16.10^{-5} \\\\) | 0.50 | 0\n",
    "200-240 | \\\\( 4.10^{-5} \\\\) | 0.25 | 0\n",
    "241-250 | \\\\( 10^{-5} \\\\) | 0.25 | 0\n",
    "251-300 | \\\\( 64.10^{-5} \\\\) | 1.00 | \\\\( 2.10^{-4} \\\\)\n",
    "301-350 | \\\\( 16.10^{-5} \\\\) | 0.50 | \\\\( 2.10^{-4} \\\\)\n",
    "351-390 | \\\\( 4.10^{-5} \\\\) | 0.25 | \\\\( 2.10^{-4} \\\\)\n",
    "391-400 | \\\\( 10^{-5} \\\\) | 0.25 | \\\\( 2.10^{-4} \\\\)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "d8878aaebeb16efc57a1628dc47c0a811b7c9a6e"
   },
   "outputs": [],
   "source": [
    "# A Keras generator to evaluate only the BRANCH MODEL\n",
    "class FeatureGen(Sequence):\n",
    "    def __init__(self, data, batch_size=64, verbose=1):\n",
    "        super(FeatureGen, self).__init__()\n",
    "        self.data       = data\n",
    "        self.batch_size = batch_size\n",
    "        self.verbose    = verbose\n",
    "        if self.verbose > 0: self.progress = tqdm_notebook(total=len(self), desc='Features')\n",
    "    def __getitem__(self, index):\n",
    "        start = self.batch_size*index\n",
    "        size  = min(len(self.data) - start, self.batch_size)\n",
    "        a     = np.zeros((size,) + img_shape, dtype=K.floatx())\n",
    "        for i in range(size): a[i,:,:,:] = read_for_validation(self.data[start + i])\n",
    "        if self.verbose > 0: \n",
    "            self.progress.update()\n",
    "            if self.progress.n >= len(self): self.progress.close()\n",
    "        return a\n",
    "    def __len__(self):\n",
    "        return (len(self.data) + self.batch_size - 1)//self.batch_size\n",
    "\n",
    "# A Keras generator to evaluate on the HEAD MODEL on features already pre-computed.\n",
    "# It computes only the upper triangular matrix of the cost matrix if y is None.\n",
    "class ScoreGen(Sequence):\n",
    "    def __init__(self, x, y=None, batch_size=2048, verbose=1):\n",
    "        super(ScoreGen, self).__init__()\n",
    "        self.x          = x\n",
    "        self.y          = y\n",
    "        self.batch_size = batch_size\n",
    "        self.verbose    = verbose\n",
    "        if y is None:\n",
    "            self.y           = self.x\n",
    "            self.ix, self.iy = np.triu_indices(x.shape[0],1)\n",
    "        else:\n",
    "            self.iy, self.ix = np.indices((y.shape[0],x.shape[0]))\n",
    "            self.ix          = self.ix.reshape((self.ix.size,))\n",
    "            self.iy          = self.iy.reshape((self.iy.size,))\n",
    "        self.subbatch = (len(self.x) + self.batch_size - 1)//self.batch_size\n",
    "        if self.verbose > 0: self.progress = tqdm_notebook(total=len(self), desc='Scores')\n",
    "    def __getitem__(self, index):\n",
    "        start = index*self.batch_size\n",
    "        end   = min(start + self.batch_size, len(self.ix))\n",
    "        a     = self.y[self.iy[start:end],:]\n",
    "        b     = self.x[self.ix[start:end],:]\n",
    "        if self.verbose > 0: \n",
    "            self.progress.update()\n",
    "            if self.progress.n >= len(self): self.progress.close()\n",
    "        return [a,b]\n",
    "    def __len__(self):\n",
    "        return (len(self.ix) + self.batch_size - 1)//self.batch_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install keras_tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "f0835cea78c2d2c63ef290d5ff826f3c2a866504"
   },
   "outputs": [],
   "source": [
    "from keras_tqdm import TQDMNotebookCallback\n",
    "\n",
    "def set_lr(model, lr):\n",
    "    K.set_value(model.optimizer.lr, float(lr))\n",
    "\n",
    "def get_lr(model):\n",
    "    return K.get_value(model.optimizer.lr)\n",
    "\n",
    "def score_reshape(score, x, y=None):\n",
    "    \"\"\"\n",
    "    Tranformed the packed matrix 'score' into a square matrix.\n",
    "    @param score the packed matrix\n",
    "    @param x the first image feature tensor\n",
    "    @param y the second image feature tensor if different from x\n",
    "    @result the square matrix\n",
    "    \"\"\"\n",
    "    if y is None:\n",
    "        # When y is None, score is a packed upper triangular matrix.\n",
    "        # Unpack, and transpose to form the symmetrical lower triangular matrix.\n",
    "        m = np.zeros((x.shape[0],x.shape[0]), dtype=K.floatx())\n",
    "        m[np.triu_indices(x.shape[0],1)] = score.squeeze()\n",
    "        m += m.transpose()\n",
    "    else:\n",
    "        m        = np.zeros((y.shape[0],x.shape[0]), dtype=K.floatx())\n",
    "        iy,ix    = np.indices((y.shape[0],x.shape[0]))\n",
    "        ix       = ix.reshape((ix.size,))\n",
    "        iy       = iy.reshape((iy.size,))\n",
    "        m[iy,ix] = score.squeeze()\n",
    "    return m\n",
    "\n",
    "def compute_score(verbose=1):\n",
    "    \"\"\"\n",
    "    Compute the score matrix by scoring every pictures from the training set against every other picture O(n^2).\n",
    "    \"\"\"\n",
    "    features = branch_model.predict_generator(FeatureGen(train, verbose=verbose), max_queue_size=12, workers=6, verbose=0)\n",
    "    score    = head_model.predict_generator(ScoreGen(features, verbose=verbose), max_queue_size=12, workers=6, verbose=0)\n",
    "    score    = score_reshape(score, features)\n",
    "    return features, score\n",
    "\n",
    "def make_steps(step, ampl):\n",
    "    \"\"\"\n",
    "    Perform training epochs\n",
    "    @param step Number of epochs to perform\n",
    "    @param ampl the K, the randomized component of the score matrix.\n",
    "    \"\"\"\n",
    "    global w2ts, t2i, steps, features, score, histories\n",
    "    \n",
    "    # shuffle the training pictures\n",
    "    random.shuffle(train)\n",
    "    \n",
    "    # Map whale id to the list of associated training picture hash value\n",
    "    w2ts = {}\n",
    "    for w,hs in w2hs.items():\n",
    "        for h in hs:\n",
    "            if h in train_set:\n",
    "                if w not in w2ts: w2ts[w] = []\n",
    "                if h not in w2ts[w]: w2ts[w].append(h)\n",
    "    for w,ts in w2ts.items(): w2ts[w] = np.array(ts)\n",
    "\n",
    "    # Map training picture hash value to index in 'train' array    \n",
    "    t2i  = {}\n",
    "    for i,t in enumerate(train): t2i[t] = i    \n",
    "\n",
    "    # Compute the match score for each picture pair\n",
    "    features, score = compute_score()\n",
    "    \n",
    "    # Train the model for 'step' epochs\n",
    "    history = model.fit_generator(\n",
    "        TrainingData(score + ampl*np.random.random_sample(size=score.shape), steps=step, batch_size=32),\n",
    "        initial_epoch=steps, epochs=steps + step, max_queue_size=12, workers=6, verbose=0,\n",
    "        callbacks=[\n",
    "            TQDMNotebookCallback(leave_inner=True, metric_format='{value:0.3f}')\n",
    "        ]).history\n",
    "    steps += step\n",
    "    \n",
    "    # Collect history data\n",
    "    history['epochs'] = steps\n",
    "    history['ms'    ] = np.mean(score)\n",
    "    history['lr'    ] = get_lr(model)\n",
    "    print(history['epochs'],history['lr'],history['ms'])\n",
    "    histories.append(history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "f84b80385ff2adbac2528fe5642f55c5bb50ea29"
   },
   "outputs": [],
   "source": [
    "model_name = 'mpiotte-standard'\n",
    "histories  = []\n",
    "steps      = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "2563e2e6b15e7ce23f94ad5a1093ad2e70861de4"
   },
   "outputs": [],
   "source": [
    "if isfile('../input/humpback-whale-identification-model-files/mpiotte-standard.model'):\n",
    "    tmp = keras.models.load_model('../input/humpback-whale-identification-model-files/mpiotte-standard.model')\n",
    "    model.set_weights(tmp.get_weights())\n",
    "else:\n",
    "    # epoch -> 10\n",
    "    make_steps(10, 1000)\n",
    "    ampl = 100.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for _ in range(10):\n",
    "    print('noise ampl.  = ', ampl)\n",
    "    make_steps(5, ampl)\n",
    "    ampl = max(1.0, 100**-0.1*ampl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# epoch -> 150\n",
    "for _ in range(18): make_steps(5, 1.0)\n",
    "# epoch -> 200\n",
    "set_lr(model, 16e-5)\n",
    "for _ in range(10): make_steps(5, 0.5)\n",
    "# epoch -> 240\n",
    "set_lr(model, 4e-5)\n",
    "for _ in range(8): make_steps(5, 0.25)\n",
    "# epoch -> 250\n",
    "set_lr(model, 1e-5)\n",
    "for _ in range(2): make_steps(5, 0.25)\n",
    "# epoch -> 300\n",
    "weights = model.get_weights()\n",
    "model, branch_model, head_model = build_model(64e-5,0.0002)\n",
    "model.set_weights(weights)\n",
    "for _ in range(10): make_steps(5, 1.0)\n",
    "# epoch -> 350\n",
    "set_lr(model, 16e-5)\n",
    "for _ in range(10): make_steps(5, 0.5)    \n",
    "# epoch -> 390\n",
    "set_lr(model, 4e-5)\n",
    "for _ in range(8): make_steps(5, 0.25)\n",
    "# epoch -> 400\n",
    "set_lr(model, 1e-5)\n",
    "for _ in range(2): make_steps(5, 0.25)\n",
    "model.save('mpiotte-standard.model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "a02a1fa6e97259186174afc82a4599e03d047d2b"
   },
   "source": [
    "#  Generating the submission file\n",
    "The basic strategy is this. For each picture from the test set:\n",
    "\n",
    "1. If the image duplicates one or more training set images, add the whales (possibly more than one) from the training image as the top candidates.\n",
    "1. For each image not new_whale from the training set, compute the image score, which is the model score for the image pair.\n",
    "1. For each whale from the training set, compute the whale score as the maximum image score for this whale.\n",
    "1. Add new_whale with a fixed new whale score of 'threshold'.\n",
    "1. Sort the whales in decreasing score.\n",
    "\n",
    "Duplicate image are free answers, assuming there is no tagging error. For new_whale, the algorithm will chose high confidence prediction first, then insert new_whale, then low confidence predictions. 'threshold' is selected through trial and error, although most model variants perform best with a value of 'threshold' that results in ~7100 test images with new_whale as first choice, something that can be measured without submitting a prediction to Kaggle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "74332a27926d3bf1b72f7b0d027baf6da053a52b"
   },
   "outputs": [],
   "source": [
    "# Not computing the submission in this notebook because it is a little slow. It takes about 15 minutes on setup with a GTX 1080.\n",
    "import gzip\n",
    "\n",
    "def prepare_submission(threshold, filename):\n",
    "    \"\"\"\n",
    "    Generate a Kaggle submission file.\n",
    "    @param threshold the score given to 'new_whale'\n",
    "    @param filename the submission file name\n",
    "    \"\"\"\n",
    "    vtop  = 0\n",
    "    vhigh = 0\n",
    "    pos   = [0,0,0,0,0,0]\n",
    "    with gzip.open(filename, 'wt', newline='\\n') as f:\n",
    "        f.write('Image,Id\\n')\n",
    "        for i,p in enumerate(tqdm_notebook(submit)):\n",
    "            t = []\n",
    "            s = set()\n",
    "            a = score[i,:]\n",
    "            for j in list(reversed(np.argsort(a))):\n",
    "                h = known[j]\n",
    "                if a[j] < threshold and new_whale not in s:\n",
    "                    pos[len(t)] += 1\n",
    "                    s.add(new_whale)\n",
    "                    t.append(new_whale)\n",
    "                    if len(t) == 5: break;\n",
    "                for w in h2ws[h]:\n",
    "                    assert w != new_whale\n",
    "                    if w not in s:\n",
    "                        if a[j] > 1.0:\n",
    "                            vtop += 1\n",
    "                        elif a[j] >= threshold:\n",
    "                            vhigh += 1\n",
    "                        s.add(w)\n",
    "                        t.append(w)\n",
    "                        if len(t) == 5: break;\n",
    "                if len(t) == 5: break;\n",
    "            if new_whale not in s: pos[5] += 1\n",
    "            assert len(t) == 5 and len(s) == 5\n",
    "            f.write(p + ',' + ' '.join(t[:5]) + '\\n')\n",
    "    return vtop,vhigh,pos\n",
    "\n",
    "if False:\n",
    "    # Find elements from training sets not 'new_whale'\n",
    "    h2ws = {}\n",
    "    for p,w in tagged.items():\n",
    "        if w != new_whale: # Use only identified whales\n",
    "            h = p2h[p]\n",
    "            if h not in h2ws: h2ws[h] = []\n",
    "            if w not in h2ws[h]: h2ws[h].append(w)\n",
    "    known = sorted(list(h2ws.keys()))\n",
    "\n",
    "    # Dictionary of picture indices\n",
    "    h2i   = {}\n",
    "    for i,h in enumerate(known): h2i[h] = i\n",
    "\n",
    "    # Evaluate the model.\n",
    "    fknown  = branch_model.predict_generator(FeatureGen(known), max_queue_size=20, workers=10, verbose=0)\n",
    "    fsubmit = branch_model.predict_generator(FeatureGen(submit), max_queue_size=20, workers=10, verbose=0)\n",
    "    score   = head_model.predict_generator(ScoreGen(fknown, fsubmit), max_queue_size=20, workers=10, verbose=0)\n",
    "    score   = score_reshape(score, fknown, fsubmit)\n",
    "\n",
    "    # Generate the subsmission file.\n",
    "    prepare_submission(0.99, 'mpiotte-standard.csv.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "eb8f53368668723f327040bb3b058f005d38ba77"
   },
   "source": [
    "# Bootstrapping and ensemble\n",
    "The mpiotte-standard.model is good for a 0.766 score.\n",
    "\n",
    "Because the training data set is small, and the test set is larger, bootstrapping is a good candidate to improve the score. In this context, bootstrapping means using the model to automatically generate additional training example, and retrained the model over this larger dataset. For this experiment, I selected all images from the test set for which the mpiotte-standard model provides a single whale prediction with a score > 0.999999 (score is just a number use to rank whales, it is **not** a probability).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "b832dd2d0d837eea499b129184f6633e193ade48"
   },
   "source": [
    "## Bootstrapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "fcef7ef0b197305f6944e99346451cd3bca63fce"
   },
   "outputs": [],
   "source": [
    "with open('../input/humpback-whale-identification-model-files/bootstrap.pickle', 'rb') as f:\n",
    "    bootstrap = pickle.load(f)\n",
    "len(bootstrap), list(bootstrap.items())[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "42062beebd3413eeced9706ddf629406dba773f3"
   },
   "source": [
    "Submitting these 1885 pictures as a submission show that this set if over 93% accurate.\n",
    "\n",
    "Adding these files to the training set and re-running the training from scratch generates the mpiotte.bootstrap.model, also included in the dataset. This model as a slightly better score of 0.744 with a threshold=0.989."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "664531df0a9fd9cf7f59d02b47cf0a612a3d9522"
   },
   "source": [
    "## Ensemble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "f294a99f12859710df82c9c6f663d5d0d7f3aaa9"
   },
   "source": [
    "The best score is obtained by an ensemble of the mpiotte-standard.model and mpiotte-bootstrap.model. Both of these models make different errors because of their nature, which make them good candidate for ensembling:\n",
    "\n",
    "* The standard model is trained on the smallest training set, and thus has more potential for overfitting.\n",
    "* The bootstrap model is trained on more data, however the tagging accuracy is lower since the bootstrap data is only 93% accurate.\n",
    "\n",
    "The assembly strategy consist in compute a score matrix (or dimension test by train) that is a linear combination of the standard and bootstrap model. Generation of the submission using the score matrix is unchanged. Trial and error suggest a weight of 0.45 for the standard model and 0.55 for the bootstrap model.\n",
    "\n",
    "The resulting ensemble as an accuracy of **0.78563** using a threshold of 0.92. It is interesting to note how the 'threshold' value for the ensemble is much lower, which is consistent with the fact that both models make different errors, and thus the ensemble scores are typically lower than the individual models, which are grossly optimistic about they guesses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "69020ecc39f36ad16c670da837e6c336580525e8"
   },
   "source": [
    "# Visualization\n",
    "This section explores the model through some visualizations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "ca8713813b6e8d8fff5ce8c031f4d2cd05f48561"
   },
   "source": [
    "## Feature weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "3ec585afe83f6fa115c2867c8ff60f0f3f066b0a"
   },
   "source": [
    "As was discussed in the model description, the head model makes a weighted sum of the features, allowing for negative weights. We can verify that we see a combination of positive and negative weights, that confirm that some features, when matched, reduce the probability of matching whales. This could be that we match uniform, unicolor flukes, which is less likely to be correct that mathing flukes with multiple caracteristic markings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "a67b17eb5b0472789036dff08b064cb709d13389"
   },
   "outputs": [],
   "source": [
    "w = head_model.layers[-1].get_weights()[0]\n",
    "w = w.flatten().tolist()\n",
    "w = sorted(w)\n",
    "fig, axes = plt.subplots(1,1)\n",
    "axes.bar(range(len(w)), w)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "7dee1a76306cba30215867133189f58f56a6c90e"
   },
   "source": [
    "We can also check how 'per feature' network behave for different feature values.\n",
    "\n",
    "What we expect to see is that equal zero feature should produce a smaller output than similar large values. At the same time, very dissimilar values must be penalized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "2af9ec4e2c5df9faa21aba4ed37833486ade6b55"
   },
   "outputs": [],
   "source": [
    "# Construct the head model with linear activation\n",
    "_, _, tmp_model = build_model(64e-5,0, activation='linear')\n",
    "tmp_model.set_weights(head_model.get_weights())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "e937b7a76cc7e2b6e7474592a993e9d837437c03"
   },
   "outputs": [],
   "source": [
    "# Evaluate the model for constant vectors.\n",
    "a = np.ones((21*21,512),dtype=K.floatx())\n",
    "b = np.ones((21*21,512),dtype=K.floatx())\n",
    "for i in range(21):\n",
    "    for j in range(21):\n",
    "        a[21*i + j] *= float(i)/10.0\n",
    "        b[21*i + j] *= float(j)/10.0\n",
    "x    = np.arange(0.0, 2.01, 0.1, dtype=K.floatx())\n",
    "x, y = np.meshgrid(x, x)\n",
    "z    = tmp_model.predict([a,b], verbose=0).reshape((21,21))\n",
    "x.shape, y.shape, z.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "0010e75e700bfa3501d7ff8bf73c7be764ad62e0"
   },
   "source": [
    "## Pseudo-distance function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "ab38b8c09129211f8d79cfd569ba4740cb56e540"
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot_surface(x, y, z, cmap=cm.coolwarm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "1ddf0cfc7f773b64bf7c2ba1fc490b10dc74e793"
   },
   "source": [
    "Not very easy to see, but still the largest output (best whale match) occurs for matching features with large values. Matching zeros get llower value.\n",
    "\n",
    "Just the colormap is probably easier to see. This confirms are initial assumptions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "1b76640bfb5f4c466a8077d38fc928246c8949df"
   },
   "outputs": [],
   "source": [
    "from matplotlib.colors import BoundaryNorm\n",
    "from matplotlib.ticker import MaxNLocator\n",
    "\n",
    "levels = MaxNLocator(nbins=15).tick_values(z.min(), z.max())\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "cf = ax.contourf(x, y, z, levels=levels, cmap=cm.coolwarm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "ef170b3c74a41e52c64750b807d3ca56ecdda144"
   },
   "source": [
    "## Feature activation\n",
    "This section attempts to reconstruct image that maximally activate a feature. This provides some insight on the filtering process.\n",
    "\n",
    "The code to generate images is modified from examples found in [Deep Learning with Python](https://www.manning.com/books/deep-learning-with-python) by Francois Chollet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "ced4199e3979bc63c12949047ca88af4fb9c1538"
   },
   "outputs": [],
   "source": [
    "from scipy.ndimage import gaussian_filter\n",
    "\n",
    "def show_filter(filter, blur):\n",
    "    np.random.seed(1)\n",
    "    noise   = 0.1 # Initial noise\n",
    "    step    = 1 # Gradient step\n",
    "    \n",
    "    # Construct the function\n",
    "    inp     = branch_model.layers[0].get_input_at(0)\n",
    "    loss    = K.mean(branch_model.layers[-3].output[0,2:4,2:4,filter]) # Stimulate the 4 central cells\n",
    "    grads   = K.gradients(loss, inp)[0]\n",
    "    grads  /= K.sqrt(K.mean(K.square(grads))) + K.epsilon()\n",
    "    iterate = K.function([inp],[grads])\n",
    "    img     = (np.random.random(img_shape) -0.5)*noise\n",
    "    img     = np.expand_dims(img, 0)\n",
    "\n",
    "    # Use gradient descent to form image\n",
    "    for i in range(200):\n",
    "        grads_value = iterate([img])[0]\n",
    "        # Blurring a little creates nicer images by reducing reconstruction noise\n",
    "        img = gaussian_filter(img + grads_value*step, sigma=blur)\n",
    "\n",
    "    # Clip the image to improve contrast\n",
    "    avg  = np.mean(img)\n",
    "    std  = sqrt(np.mean((img - avg)**2))\n",
    "    low  = avg - 5*std\n",
    "    high = avg + 5*std\n",
    "    return array_to_img(np.minimum(high, np.maximum(low, img))[0])\n",
    "\n",
    "# Show the first 25 features (of 512)\n",
    "show_whale([show_filter(i, 0.5) for i in tqdm_notebook(range(25))], per_row=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "38ac924d4edbedd64d312c99e353828e409e8ecb",
    "collapsed": true
   },
   "source": [
    "# Off topic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "596d0a17807e10828730542ea6fd1e90a42f60af"
   },
   "source": [
    "## Training Scalability\n",
    "As described, training the base model takes about 2 days, and the bootstrap version a little under 3 days, using a i7-8700 CPU and GTX 1080 GPU. More than 50% of the time is spent solving the Linear Assignment Problem because the algorithm used has complexity \\\\( O(n^3) \\\\) and provides an exact solution. However, the score matrix is randomized, so investing a lot of time to compute an exact solution to a randomized input is wasteful. In the context of this competition, with no constraint on runtime and a small dataset, it was a pragmatic choice. However, to scale this approach, a less costly randomized matching heuristic would be more effective.\n",
    "\n",
    "Another approach to training scalability would be to partition the training data into different subset, each subset being processed separately to match image pairs. The subset can be reconstructed randomly each time the cost matrix is computed. This would be effective not only for the Linear Assignment Problem part, but also when computing the cost matrix, which still has complexity \\\\( O(n^2) \\\\). By fixing the subset size to a reasonable value, the complexity then grows linearly with the number of subsets, allowing larger training datasets..  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "ac0cbe07b888d80c4ffb7121bfd9a4f658816b65"
   },
   "source": [
    "## Interesting results and scores\n",
    "Score | Description\n",
    "---:|-------\n",
    "0.786 | best score obtained by a linear combination of the standard model with the bootstrap model\n",
    "0.774 | bootstrapped model\n",
    "0.766 | standard model\n",
    "0.752 | VGG like CNN, trained like standard model\n",
    "0.728 | standard model without L2 regularization (result after 250 epochs)\n",
    "0.714 | standard model without exclusion list, rotation list and bounding box model (i.e. no manual judgement on the training set)\n",
    "0.423 | duplicate images and new_whale submission\n",
    "0.325 | new_whale submission\n",
    "0.107 | duplicate images only\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "1ad0edb6a1fca52effac15d40addd8560e854f28"
   },
   "source": [
    "## Validation set\n",
    "I have not discussed validation data until now. During development, I used a validation set composed of 570 images from the training set to test ideas and to tune the training procedure. However, greater accuracy can be achieved by retraining the model using all the data,  repeating he procedure that was successful on the validation set. This notebook essentially describes this final retraining, and thus there is no validation set involved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_uuid": "4b8ee534a637334e163acb2031cce28475a71d7c"
   },
   "outputs": [],
   "source": [
    "!pip install kaggle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}