naive_bayes_ml.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Machine Learning - HW1\n",
    "mr3523"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive Bayes Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import scipy.io as sio\n",
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Load data.\n",
    "data = sio.loadmat('data/mnist_mat.mat')\n",
    "ytrain = data['ytrain']\n",
    "Xtrain = data['Xtrain']\n",
    "ytest = data['ytest']\n",
    "Xtest = data['Xtest']\n",
    "\n",
    "## Define parameters.\n",
    "a = 1; b = 1; c = 1; d = 1; e = 1; f = 1\n",
    "d = 1/a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Prior function.\n",
    "def y_prior(y_s, e, f, ytrain):\n",
    "    ind = np.sum(ytrain == y_s)\n",
    "    N = len(ytrain[0])\n",
    "    prior = float(e + ind)/(N + e + f)\n",
    "    return prior"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "## Likelihood function.\n",
    "def log_likelihood(x_s, y_s, b, c, d, Xtrain, ytrain):\n",
    "    \n",
    "    ## First term (numerator).\n",
    "    N = Xtrain[np.where(ytrain == y_s)].shape[0]\n",
    "    alpha1 = N/2. + b + 0.5\n",
    "    n1 = alpha - 1\n",
    "    numerator = n1 * np.log(n) - n1 + np.log(n*(1+4*n*(1+2*n)))/6 + np.log(np.pi)/2\n",
    "    \n",
    "    ## Second term (denominator).\n",
    "    dat = Xtrain[:, np.where(ytrain[0] == y_s)].T\n",
    "    xi_sum = np.sum(dat, axis = 0)\n",
    "    beta = ((x_s**2 + xi_sum**2) / 2) - (x_s + xi_sum)**2 / (2*(1 + N + d)) + c\n",
    "    denominator = alpha * np.sum(np.log(beta))\n",
    "    \n",
    "    ## Output.\n",
    "    result = dat.shape[2] * (numerator - 0.5*np.log(1 + N + d)) - denominator\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Likelihood function.\n",
    "def log_likelihood(x_s, y_s, b, c, d, Xtrain, ytrain):\n",
    "    \n",
    "    ## Data and sample statistics.\n",
    "    dat = Xtrain[:, np.where(ytrain[0] == y_s)].T\n",
    "    s = np.sqrt(np.var(dat, axis = 0))\n",
    "    x_mean = np.mean(dat, axis = 0)\n",
    "    \n",
    "    ## First term - numerator.\n",
    "    N = Xtrain[np.where(ytrain == y_s)].shape[0]\n",
    "    alpha1 = N/2. + b + 0.5\n",
    "    n1 = alpha1 - 1\n",
    "    num1 = n1 * np.log(n1) - n1 + np.log(n1*(1+4*n1*(1+2*n1)))/6 + np.log(np.pi)/2\n",
    "    \n",
    "    ## Second term - numerator.\n",
    "    num2 = 0.5 * np.log(d + N)\n",
    "    \n",
    "    ## Third term - numerator.\n",
    "    num3 = (N/2 + b) * np.log(c + 0.5 * (N * s + (d * N * x_mean**2) / (N + d)))\n",
    "    \n",
    "    ## Combining them.\n",
    "    numerator = num1 + num2 + num3\n",
    "    \n",
    "    ## First term - denominator.\n",
    "    alpha2 = N/2. + b\n",
    "    n2 = alpha2 - 1\n",
    "    dem1 = n2 * np.log(n2) - n2 + np.log(n2*(1+4*n2*(1+2*n2)))/6 + np.log(np.pi)/2\n",
    "    \n",
    "    ## Second term - denominator.\n",
    "    dem2 = 0.5 * np.log(2 * np.pi * (1 + N + d))\n",
    "    \n",
    "    ## Third term - denominator.\n",
    "    xi_sum = np.sum(dat, axis = 0)\n",
    "    beta = ((x_s**2 + xi_sum**2) / 2) - (x_s + xi_sum)**2 / (2*(1 + N + d)) + c\n",
    "    #dem3 = alpha1 * np.sum(np.log(beta))\n",
    "    dem3 = alpha1 * np.log(beta)\n",
    "    \n",
    "    ## Combining them.\n",
    "    denominator = dem1 + dem2 + dem3\n",
    "    \n",
    "    ## Output.\n",
    "    result = np.sum(numerator - denominator)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def naive_bayes(x_s, b, c, d, e, f, Xtrain, ytrain):\n",
    "    prob0 = log_likelihood(x_s, 0, b, c, d, Xtrain, ytrain) + np.log(y_prior(0, e, f, ytrain))\n",
    "    prob1 = log_likelihood(x_s, 1, b, c, d, Xtrain, ytrain) + np.log(y_prior(1, e, f, ytrain))\n",
    "    return prob1 < prob0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ypred = []\n",
    "for i in range(ytest.shape[1]):\n",
    "    ypred.append(int(naive_bayes(Xtest[:, i], b, c, d, e, f, Xtrain, ytrain)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confusion Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   0,  982],\n",
       "       [   0, 1009]])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "confusion_matrix(ytest[0], ypred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Misclassified Digits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "total size of new array must be unchanged",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-128-a6d90b1d5ec5>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mimg1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mXtest\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwhere\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mytest\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m15\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m982\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'Q'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m28\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m28\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mimg2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mXtest\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwhere\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mytest\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m15\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m982\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'Q'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m28\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m28\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mimg3\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mXtest\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwhere\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mytest\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m15\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m982\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'Q'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m28\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m28\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: total size of new array must be unchanged"
     ]
    }
   ],
   "source": [
    "img1 = np.dot(Xtest[:, np.where(ytest[0] == 0)].reshape(15, 982)[:,0], data['Q'].T).reshape(28,28)\n",
    "img2 = np.dot(Xtest[:, np.where(ytest[0] == 0)].reshape(15, 982)[:,0], data['Q'].T).reshape(28,28)\n",
    "img3 = np.dot(Xtest[:, np.where(ytest[0] == 0)].reshape(15, 982)[:,0], data['Q'].T).reshape(28,28)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD7CAYAAABKWyniAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEptJREFUeJzt3W2MXdV1xvH/ih0b/DaAiccYO+DwmjbB05aSJtBmIqeU\nUiRIlJLUiWQSQEillChqFAIfGNpKDamEShWIIiDpGEFDIKIQFQFJ0CgCwksh4xgwBhJM7BiPjcE2\nxmBje/XDXOyxmVnres59w/v5SZbv3Ofcc7aP75pz7t1n72PujoiU5X3tboCItJ4KX6RAKnyRAqnw\nRQqkwhcpkApfpEAtK3wzO8PMnjWz58zsG63abr3MbKWZLTWzX5nZYx3QnpvMbMjMfj3iuUPN7H4z\nW2Fm95lZV4e170ozW21mT9b+nNHG9s01swfM7GkzW2Zm/1h7viP24Sjtu6T2fEv2obWiH9/M3gc8\nBywE1gCPA19w92ebvvE6mdlvgT9x99fa3RYAMzsN2AIscfeTas9dDWxw92/Xfnke6u6XdVD7rgRe\nd/dr2tGmkcxsNjDb3QfNbBrwBHA28GU6YB8G7fs8LdiHrTrinwI87+4vufvbwA8Z/kd2EqODPvq4\n+4PAvr+Ezgb6a4/7gXNa2qgRxmgfDO/HtnP3te4+WHu8BVgOzKVD9uEY7TuyFjd9H7bqjX4ksGrE\nz6vZ84/sFA7cZ2aPm9mF7W7MGGa5+xAMv3GAWW1uz2guNrNBM7uxnR9FRjKzo4Ee4BGgu9P24Yj2\nPVp7qun7sGOOcB3gVHc/GTiT4R1/WrsbVIdOu976euAYd+8B1gKdcMo/DbgDuLR2ZN13n7V1H47S\nvpbsw1YV/u+BD474eW7tuY7h7i/X/l4P3Mnwx5NOM2Rm3bD7M+K6NrdnL+6+3vd8aXQD8KftbI+Z\nTWS4qG5297tqT3fMPhytfa3ah60q/MeBY83sKDObBHwBuLtF206Z2ZTab17MbCpwOvBUe1sFDH/W\nG/l5727gvNrjxcBd+76gxfZqX62Q3vFZ2r8Pvw884+7Xjniuk/bhu9rXqn3Ykm/1Ybg7D7iW4V82\nN7n7t1qy4TqY2XyGj/IOTARuaXf7zOxWoBeYCQwBVwL/A9wOzANeAs51940d1L5PMfxZdRewErjo\nnc/TbWjfqcAvgGUM/786cDnwGPAj2rwPg/YtogX7sGWFLyKdQ1/uiRRIhS9SIBW+SIFU+CIFqlT4\nnT7wRkRGN+5v9esdeGNm6jYQaRN3H/W6/4kV1rl74A2Amb0z8OZdI+5+7h/f/bi/bxWL++ZV2Gxz\nqX3VdHL7Orlt0Pj2LbRfjplVOdV/Lwy8EZFR6Ms9kQJVOdWve+BNf9+eE4Oph0yosMnmW9A7o91N\nCKl949fJbYPq7Rsc2MTSgc11LVvly70JwAqGv9x7meFroP/O3Zfvs5yP/IwvIq2x0H7Z+C/33H2n\nmf0DcD97Bt4sT14mIh2gyqk+7n4vcEKD2iIiLaIv90QKpMIXKZAKX6RAKnyRAqnwRQqkwhcpkApf\npEAqfJECqfBFCqTCFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAlcbji3S6Cexs6vp30tlTyY1FR3yR\nAqnwRQqkwhcpkApfpEAqfJECqfBFCqTCFylQ8f34Vfths37iPN8R5hMrrj+zjUlhvjN5i2T7L8ur\n7r/MdF4P88lsq7T+rUwJ82z/bmdype2Pl474IgVS4YsUSIUvUiAVvkiBVPgiBVLhixRIhS9SoEr9\n+Ga2EtgE7ALedvdTGtGoRqraz5zJ+pmnsDXM57AmzKfvjPuhd0yI27+RQ8J8Hd1hnvUzb+XgMH8z\n6efekez/7DqGmWwI826GwnzWplfDPLOhqyvM1zAnzKteJzFeVS/g2QX0uvtrjWiMiLRG1VN9a8A6\nRKTFqhatA/eZ2eNmdmEjGiQizVf1VP9Ud3/ZzD4A/NTMlrv7g/su1N+3avfjBb0z6OmNPxeJyP4b\nHNjE0oHNdS1bqfDd/eXa3+vN7E7gFOBdhb+4b16VzYhIHXp6u/Y6qC65avWYy477VN/MppjZtNrj\nqcDpwFPjXZ+ItE6VI343cKeZeW09t7j7/Y1plog007gL391fBHoa2JZxqdpPn/UjZ/2sB/NmmM/k\nlTA/8tmkH/m3ccyst8N4+8lx+17g2DB/mE+E+WDyFniRo8N8SrL/PsajYb6IW8L8yHXJ/n0xjpka\nx9MPSsb7T47H+7/O9KQBzaGuOJECqfBFCqTCFymQCl+kQCp8kQKp8EUKpMIXKdABP69+1k+fjTev\nOh56/qtr4wXuS1awKsnPi+Pb+HyYX8ffh/nzKxbEG3ghjtP8oDjeeEE8n8BnJtwZr2B5sv14OgQ4\nLo53TozfH9vaNG9+Rkd8kQKp8EUKpMIXKZAKX6RAKnyRAqnwRQqkwhcp0AHfj5+Np696f/js/vY2\nGMaQTV3ysTj+zkfOD/OvPvG9eAXvmihtH7OTfFqSv5XkY88OBcDWLfG8/bO61sUryMbbx/99kEwP\n+cqEmWGe3XegWfPmZ3TEFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAKnyRAhXQjx/3k2b3b89e303S\nj/zzOE5u787mK94f5pf85oZ4Bd9Jtn9aHHedE88nMGny9jBfv/GD8QaSfv4/6oovhPiLVx+LV5Ds\nXz4ax68eF08YsI7uMM/me2gXHfFFCqTCFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAaT++md0EnAUM\nuftJtecOBW4DjgJWAue6+6YmtnNMWT97Nt5+a9KPP534/ufHrUoGlD8Sx/xZHH99wr/HC/yHxfmJ\ncfyh858O8z/nF2G+nD8I8/Vbkn78ZLz/57gjzO32+PWpj8TxSuaHedX7209gZ6XXj1c9R/wfAH+1\nz3OXAT9z9xOAB4BvNrphItI8aeG7+4PAa/s8fTbQX3vcD5zT4HaJSBON9zP+LHcfAnD3tcCsxjVJ\nRJqtUdfqexT29+25AdyC3hn09CYTmYnIfhsc2MTSgc11LTvewh8ys253HzKz2RCPVFncN2+cmxGR\nevX0du11UF1y1dhfPNd7qm+1P++4mz33aV0M3LVfLRSRtkoL38xuBR4Gjjez35nZl4FvAX9pZiuA\nhbWfReQ9Ij3Vd/dFY0SfbnBbmiIbD51dBzCTV+INZPPSZ/e3vzqObxy6IF4gub/8Yf/0+zBfxC1h\nPpl4vP09/E3cgJ/FMV+N4wt23hgvkO3/v43jp+d8KMzXcESygVi2/9pFV+6JFEiFL1IgFb5IgVT4\nIgVS4YsUSIUvUiAVvkiB3vPz6u9I+uG3J+Pxs/HQ81kZN+CJOGZqHD9y8oIw3/XjZAXJvPhnTrgn\nzLP7AjzMJ8J8/dXJePukH//Dtz8Z5jOueztewbY4fmthnP86mVg/uw7kEDbGG1A/voh0ChW+SIFU\n+CIFUuGLFEiFL1IgFb5IgVT4IgXq+H78bLx8VVPYGuZHsCZeQXb/9Z44vu9dM5fvY2Wy/mRe+Oc4\nPsyXJf3YS69LJv6/7MU4nxjPS7+IW+PX/18c88k4vm/q6WG+jJPC/PBkPoZJyYUEk7MLDSoa77z8\nOuKLFEiFL1IgFb5IgVT4IgVS4YsUSIUvUiAVvkiBOr4fv6qDeTPMs378Q199K95AMlw+6SbnmeT+\n8jySrD+Zdv6xiUlH91PJ+vmvJF8Zxxf0hfHnuCN+fbJ/3zovzrPrJDYwM8wnsCPMs/H403k9WX/c\nDz/efvqMjvgiBVLhixRIhS9SIBW+SIFU+CIFUuGLFEiFL1KgtB/fzG4CzgKG3P2k2nNXAhfC7knZ\nL3f3e5vRwKwfc2KSZ/30WT+sJd34dCV5PC07x7MiXiAZz591g8PzSR5f5wB/HcfHdofxh78bz5t/\n4t0vxet/I45XTD0uzLN++qrzPTSrn73Z6jni/wBGvQriGnf/49qfphS9iDRHWvju/iDw2iiRNb45\nItIKVT7jX2xmg2Z2o5llJ7wi0kHGe63+9cA/u7ub2b8C1wDnj7Vwf9+q3Y8X9M6gp1e/J0QabXBg\nE0sHNte17LgK393Xj/jxBuAn0fKL++aNZzMish96erv2OqguuWr1mMvWe6pvjPhMb2azR2SfpY4x\nXiLSOerpzrsV6AVmmtnvgCuBT5lZD7CL4XGZFzWxjSLSYGnhu/uiUZ7+QRPaMi5ZP2o2HjrLU9l4\n/GRe+H85+d/CfP4VK8P8nivODPNXODzMs+schog/M65L8m/z9TDnG3GcrJ7u3ZeSjG4eq8I8Mye5\nr0L2/smuM2kXXbknUiAVvkiBVPgiBVLhixRIhS9SIBW+SIFU+CIFOuDn1Z/E9iRP7l+e7aFZST6Y\n5BfE8Vcm/necT41z5iTbXxjH934pnpf/h3w+zM969oEw/+6z8fZPiWNO2rQpzOd0xf3wmcPZEObZ\ndRDZ+69ddMQXKZAKX6RAKnyRAqnwRQqkwhcpkApfpEAqfJECvef78avOa74z2QVrZ8XzA85eGPcj\nc1DSgEeS/MUkz6aF/2Acv/WZOL8xudDg13w0XkHcDc5QHHNccl+CdV2HhfnrTA/zqvMxHMjz6ovI\nAUaFL1IgFb5IgVT4IgVS4YsUSIUvUiAVvkiB3vP9+Jns/udbOTjM1yQD2jfOPyTMZ86PO7IPP3NL\nmNtbYcz6OdPCfJCeML8tGU//44e+FDdg7ttxnrzDZsQxM5LrEF6gO8y3MiXMs/HyO5L3T/b+6tR+\nfh3xRQqkwhcpkApfpEAqfJECqfBFCqTCFymQCl+kQGk/vpnNBZYA3cAu4AZ3/08zOxS4DTgKWAmc\n6+7J4PTOs514wHeWZ+O9h5J+5p2Hxf3AG4mvE1jJ0WG+ghPC/H85M8zT+QLOSd5CyTsivooCiKdD\nSPvRM9nrs///rJ9/Qof289dzxN8BfM3d/xD4OHCxmZ0IXAb8zN1PAB4Avtm8ZopII6WF7+5r3X2w\n9ngLsByYC5wN9NcW6wfOaVYjRaSx9uszvpkdDfQwfALY7e5DMPzLgfxmUiLSIeq+Vt/MpgF3AJe6\n+xYz830W2ffn3fr7Vu1+vKB3Bj29yQc3EdlvgwObWDqwua5l6yp8M5vIcNHf7O531Z4eMrNudx8y\ns9nAurFev7hvXl2NEZHx6+nt2uuguuSq1WMuW++p/veBZ9z92hHP3Q2cV3u8GLhr3xeJSGeqpzvv\nVOCLwDIz+xXDp/SXA1cDPzKzrwAvAec2s6Ei0jhp4bv7Q4w9e/unG9ucxqvaz5vZzqQw31bxOoA1\nHBHmLyfzBWTzCax9KfkYtjGOu+YmM+P/Mo7jfz0kuy+9P32WZ7L3T3ZfBpLx/u2iK/dECqTCFymQ\nCl+kQCp8kQKp8EUKpMIXKZAKX6RAB/y8+s2WjcfO5u3Pxttv5NBk/fG88RuYGeasfn+cJ+ZNXhUv\n8EYcz620dTgkudAgmze/VDriixRIhS9SIBW+SIFU+CIFUuGLFEiFL1IgFb5IgdSPX1E2Hjuftz8e\nz5+NB9+WvD67ToBX4phpcTyHNfECU+P42GTz2etnvvFqmB8y9bUwz+ZD6NT721elI75IgVT4IgVS\n4YsUSIUvUiAVvkiBVPgiBVLhixRI/fiJqvPyZ/3A2Xjx6bwe5mk/fSbpp+fEOM7al71+7unJ9k+O\n461TDwrz6WwJ8+w6DPXji8gBQ4UvUiAVvkiBVPgiBVLhixRIhS9SoLTwzWyumT1gZk+b2TIzu6T2\n/JVmttrMnqz9OaP5zRWRRqinH38H8DV3HzSzacATZvbTWnaNu1/TvOZV1+x+2AnsCPPJbAvzrB88\ne312ncEx/CbMty6M5+WflGx/HvG8+i+d/IEwP+ry9WFOdxy/zJwwn8LWMM/2b3adRfb6Tr0OIC18\nd18LrK093mJmy4Eja7E1sW0i0iT79RnfzI4GeoBHa09dbGaDZnajmXU1uG0i0iR1F37tNP8O4FJ3\n3wJcDxzj7j0MnxF09Cm/iOxR17X6ZjaR4aK/2d3vAnD3kR/ObgB+Mtbr+/v2fA5c0DuDnl6dHIg0\n2uDAJpYObK5r2XoH6XwfeMbdr33nCTObXfv8D/BZ4KmxXry4b16dmxGR8erp7drroLrkqtVjLpsW\nvpmdCnwRWGZmvwIcuBxYZGY9wC5gJXBRpVaLSMvU863+QzBqn9G9jW+OiLSCxuMnsn7YyWk/b5xn\n93fPZP3ox7MizLck88pnZrEuzLN+9jc/GV9HkN03YAOHh3nWj5/J/v87tZ8+o0t2RQqkwhcpkApf\npEAqfJECqfBFCqTCFymQCl+kQOrHr6jd/bjZ9ruTfvYsr2orcT99lld1oPbDV6UjvkiBVPgiBVLh\nixSo5YU/OLCp1ZvcL2pfNZ3cvqUD1cZFNFsr913LC7/eiQLaRe2rppPbt7SDfylBa/edTvVFCqTC\nFymQuXtzN2DW3A2IyJjcfdQp8Jte+CLSeXSqL1IgFb5IgVT4IgVS4YsUSIUvUqD/B/OL1pQWhW2B\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa9092f72d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.matshow(img1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD7CAYAAABKWyniAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEptJREFUeJzt3W2MXdV1xvH/ih0b/DaAiccYO+DwmjbB05aSJtBmIqeU\nUiRIlJLUiWQSQEillChqFAIfGNpKDamEShWIIiDpGEFDIKIQFQFJ0CgCwksh4xgwBhJM7BiPjcE2\nxmBje/XDXOyxmVnres59w/v5SZbv3Ofcc7aP75pz7t1n72PujoiU5X3tboCItJ4KX6RAKnyRAqnw\nRQqkwhcpkApfpEAtK3wzO8PMnjWz58zsG63abr3MbKWZLTWzX5nZYx3QnpvMbMjMfj3iuUPN7H4z\nW2Fm95lZV4e170ozW21mT9b+nNHG9s01swfM7GkzW2Zm/1h7viP24Sjtu6T2fEv2obWiH9/M3gc8\nBywE1gCPA19w92ebvvE6mdlvgT9x99fa3RYAMzsN2AIscfeTas9dDWxw92/Xfnke6u6XdVD7rgRe\nd/dr2tGmkcxsNjDb3QfNbBrwBHA28GU6YB8G7fs8LdiHrTrinwI87+4vufvbwA8Z/kd2EqODPvq4\n+4PAvr+Ezgb6a4/7gXNa2qgRxmgfDO/HtnP3te4+WHu8BVgOzKVD9uEY7TuyFjd9H7bqjX4ksGrE\nz6vZ84/sFA7cZ2aPm9mF7W7MGGa5+xAMv3GAWW1uz2guNrNBM7uxnR9FRjKzo4Ee4BGgu9P24Yj2\nPVp7qun7sGOOcB3gVHc/GTiT4R1/WrsbVIdOu976euAYd+8B1gKdcMo/DbgDuLR2ZN13n7V1H47S\nvpbsw1YV/u+BD474eW7tuY7h7i/X/l4P3Mnwx5NOM2Rm3bD7M+K6NrdnL+6+3vd8aXQD8KftbI+Z\nTWS4qG5297tqT3fMPhytfa3ah60q/MeBY83sKDObBHwBuLtF206Z2ZTab17MbCpwOvBUe1sFDH/W\nG/l5727gvNrjxcBd+76gxfZqX62Q3vFZ2r8Pvw884+7Xjniuk/bhu9rXqn3Ykm/1Ybg7D7iW4V82\nN7n7t1qy4TqY2XyGj/IOTARuaXf7zOxWoBeYCQwBVwL/A9wOzANeAs51940d1L5PMfxZdRewErjo\nnc/TbWjfqcAvgGUM/786cDnwGPAj2rwPg/YtogX7sGWFLyKdQ1/uiRRIhS9SIBW+SIFU+CIFqlT4\nnT7wRkRGN+5v9esdeGNm6jYQaRN3H/W6/4kV1rl74A2Amb0z8OZdI+5+7h/f/bi/bxWL++ZV2Gxz\nqX3VdHL7Orlt0Pj2LbRfjplVOdV/Lwy8EZFR6Ms9kQJVOdWve+BNf9+eE4Oph0yosMnmW9A7o91N\nCKl949fJbYPq7Rsc2MTSgc11LVvly70JwAqGv9x7meFroP/O3Zfvs5yP/IwvIq2x0H7Z+C/33H2n\nmf0DcD97Bt4sT14mIh2gyqk+7n4vcEKD2iIiLaIv90QKpMIXKZAKX6RAKnyRAqnwRQqkwhcpkApf\npEAqfJECqfBFCqTCFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAlcbji3S6Cexs6vp30tlTyY1FR3yR\nAqnwRQqkwhcpkApfpEAqfJECqfBFCqTCFylQ8f34Vfths37iPN8R5hMrrj+zjUlhvjN5i2T7L8ur\n7r/MdF4P88lsq7T+rUwJ82z/bmdype2Pl474IgVS4YsUSIUvUiAVvkiBVPgiBVLhixRIhS9SoEr9\n+Ga2EtgE7ALedvdTGtGoRqraz5zJ+pmnsDXM57AmzKfvjPuhd0yI27+RQ8J8Hd1hnvUzb+XgMH8z\n6efekez/7DqGmWwI826GwnzWplfDPLOhqyvM1zAnzKteJzFeVS/g2QX0uvtrjWiMiLRG1VN9a8A6\nRKTFqhatA/eZ2eNmdmEjGiQizVf1VP9Ud3/ZzD4A/NTMlrv7g/su1N+3avfjBb0z6OmNPxeJyP4b\nHNjE0oHNdS1bqfDd/eXa3+vN7E7gFOBdhb+4b16VzYhIHXp6u/Y6qC65avWYy477VN/MppjZtNrj\nqcDpwFPjXZ+ItE6VI343cKeZeW09t7j7/Y1plog007gL391fBHoa2JZxqdpPn/UjZ/2sB/NmmM/k\nlTA/8tmkH/m3ccyst8N4+8lx+17g2DB/mE+E+WDyFniRo8N8SrL/PsajYb6IW8L8yHXJ/n0xjpka\nx9MPSsb7T47H+7/O9KQBzaGuOJECqfBFCqTCFymQCl+kQCp8kQKp8EUKpMIXKdABP69+1k+fjTev\nOh56/qtr4wXuS1awKsnPi+Pb+HyYX8ffh/nzKxbEG3ghjtP8oDjeeEE8n8BnJtwZr2B5sv14OgQ4\nLo53TozfH9vaNG9+Rkd8kQKp8EUKpMIXKZAKX6RAKnyRAqnwRQqkwhcp0AHfj5+Np696f/js/vY2\nGMaQTV3ysTj+zkfOD/OvPvG9eAXvmihtH7OTfFqSv5XkY88OBcDWLfG8/bO61sUryMbbx/99kEwP\n+cqEmWGe3XegWfPmZ3TEFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAKnyRAhXQjx/3k2b3b89e303S\nj/zzOE5u787mK94f5pf85oZ4Bd9Jtn9aHHedE88nMGny9jBfv/GD8QaSfv4/6oovhPiLVx+LV5Ds\nXz4ax68eF08YsI7uMM/me2gXHfFFCqTCFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAaT++md0EnAUM\nuftJtecOBW4DjgJWAue6+6YmtnNMWT97Nt5+a9KPP534/ufHrUoGlD8Sx/xZHH99wr/HC/yHxfmJ\ncfyh858O8z/nF2G+nD8I8/Vbkn78ZLz/57gjzO32+PWpj8TxSuaHedX7209gZ6XXj1c9R/wfAH+1\nz3OXAT9z9xOAB4BvNrphItI8aeG7+4PAa/s8fTbQX3vcD5zT4HaJSBON9zP+LHcfAnD3tcCsxjVJ\nRJqtUdfqexT29+25AdyC3hn09CYTmYnIfhsc2MTSgc11LTvewh8ys253HzKz2RCPVFncN2+cmxGR\nevX0du11UF1y1dhfPNd7qm+1P++4mz33aV0M3LVfLRSRtkoL38xuBR4Gjjez35nZl4FvAX9pZiuA\nhbWfReQ9Ij3Vd/dFY0SfbnBbmiIbD51dBzCTV+INZPPSZ/e3vzqObxy6IF4gub/8Yf/0+zBfxC1h\nPpl4vP09/E3cgJ/FMV+N4wt23hgvkO3/v43jp+d8KMzXcESygVi2/9pFV+6JFEiFL1IgFb5IgVT4\nIgVS4YsUSIUvUiAVvkiB3vPz6u9I+uG3J+Pxs/HQ81kZN+CJOGZqHD9y8oIw3/XjZAXJvPhnTrgn\nzLP7AjzMJ8J8/dXJePukH//Dtz8Z5jOueztewbY4fmthnP86mVg/uw7kEDbGG1A/voh0ChW+SIFU\n+CIFUuGLFEiFL1IgFb5IgVT4IgXq+H78bLx8VVPYGuZHsCZeQXb/9Z44vu9dM5fvY2Wy/mRe+Oc4\nPsyXJf3YS69LJv6/7MU4nxjPS7+IW+PX/18c88k4vm/q6WG+jJPC/PBkPoZJyYUEk7MLDSoa77z8\nOuKLFEiFL1IgFb5IgVT4IgVS4YsUSIUvUiAVvkiBOr4fv6qDeTPMs378Q199K95AMlw+6SbnmeT+\n8jySrD+Zdv6xiUlH91PJ+vmvJF8Zxxf0hfHnuCN+fbJ/3zovzrPrJDYwM8wnsCPMs/H403k9WX/c\nDz/efvqMjvgiBVLhixRIhS9SIBW+SIFU+CIFUuGLFEiFL1KgtB/fzG4CzgKG3P2k2nNXAhfC7knZ\nL3f3e5vRwKwfc2KSZ/30WT+sJd34dCV5PC07x7MiXiAZz591g8PzSR5f5wB/HcfHdofxh78bz5t/\n4t0vxet/I45XTD0uzLN++qrzPTSrn73Z6jni/wBGvQriGnf/49qfphS9iDRHWvju/iDw2iiRNb45\nItIKVT7jX2xmg2Z2o5llJ7wi0kHGe63+9cA/u7ub2b8C1wDnj7Vwf9+q3Y8X9M6gp1e/J0QabXBg\nE0sHNte17LgK393Xj/jxBuAn0fKL++aNZzMish96erv2OqguuWr1mMvWe6pvjPhMb2azR2SfpY4x\nXiLSOerpzrsV6AVmmtnvgCuBT5lZD7CL4XGZFzWxjSLSYGnhu/uiUZ7+QRPaMi5ZP2o2HjrLU9l4\n/GRe+H85+d/CfP4VK8P8nivODPNXODzMs+schog/M65L8m/z9TDnG3GcrJ7u3ZeSjG4eq8I8Mye5\nr0L2/smuM2kXXbknUiAVvkiBVPgiBVLhixRIhS9SIBW+SIFU+CIFOuDn1Z/E9iRP7l+e7aFZST6Y\n5BfE8Vcm/necT41z5iTbXxjH934pnpf/h3w+zM969oEw/+6z8fZPiWNO2rQpzOd0xf3wmcPZEObZ\ndRDZ+69ddMQXKZAKX6RAKnyRAqnwRQqkwhcpkApfpEAqfJECvef78avOa74z2QVrZ8XzA85eGPcj\nc1DSgEeS/MUkz6aF/2Acv/WZOL8xudDg13w0XkHcDc5QHHNccl+CdV2HhfnrTA/zqvMxHMjz6ovI\nAUaFL1IgFb5IgVT4IgVS4YsUSIUvUiAVvkiB3vP9+Jns/udbOTjM1yQD2jfOPyTMZ86PO7IPP3NL\nmNtbYcz6OdPCfJCeML8tGU//44e+FDdg7ttxnrzDZsQxM5LrEF6gO8y3MiXMs/HyO5L3T/b+6tR+\nfh3xRQqkwhcpkApfpEAqfJECqfBFCqTCFymQCl+kQGk/vpnNBZYA3cAu4AZ3/08zOxS4DTgKWAmc\n6+7J4PTOs514wHeWZ+O9h5J+5p2Hxf3AG4mvE1jJ0WG+ghPC/H85M8zT+QLOSd5CyTsivooCiKdD\nSPvRM9nrs///rJ9/Qof289dzxN8BfM3d/xD4OHCxmZ0IXAb8zN1PAB4Avtm8ZopII6WF7+5r3X2w\n9ngLsByYC5wN9NcW6wfOaVYjRaSx9uszvpkdDfQwfALY7e5DMPzLgfxmUiLSIeq+Vt/MpgF3AJe6\n+xYz830W2ffn3fr7Vu1+vKB3Bj29yQc3EdlvgwObWDqwua5l6yp8M5vIcNHf7O531Z4eMrNudx8y\ns9nAurFev7hvXl2NEZHx6+nt2uuguuSq1WMuW++p/veBZ9z92hHP3Q2cV3u8GLhr3xeJSGeqpzvv\nVOCLwDIz+xXDp/SXA1cDPzKzrwAvAec2s6Ei0jhp4bv7Q4w9e/unG9ucxqvaz5vZzqQw31bxOoA1\nHBHmLyfzBWTzCax9KfkYtjGOu+YmM+P/Mo7jfz0kuy+9P32WZ7L3T3ZfBpLx/u2iK/dECqTCFymQ\nCl+kQCp8kQKp8EUKpMIXKZAKX6RAB/y8+s2WjcfO5u3Pxttv5NBk/fG88RuYGeasfn+cJ+ZNXhUv\n8EYcz620dTgkudAgmze/VDriixRIhS9SIBW+SIFU+CIFUuGLFEiFL1IgFb5IgdSPX1E2Hjuftz8e\nz5+NB9+WvD67ToBX4phpcTyHNfECU+P42GTz2etnvvFqmB8y9bUwz+ZD6NT721elI75IgVT4IgVS\n4YsUSIUvUiAVvkiBVPgiBVLhixRI/fiJqvPyZ/3A2Xjx6bwe5mk/fSbpp+fEOM7al71+7unJ9k+O\n461TDwrz6WwJ8+w6DPXji8gBQ4UvUiAVvkiBVPgiBVLhixRIhS9SoLTwzWyumT1gZk+b2TIzu6T2\n/JVmttrMnqz9OaP5zRWRRqinH38H8DV3HzSzacATZvbTWnaNu1/TvOZV1+x+2AnsCPPJbAvzrB88\ne312ncEx/CbMty6M5+WflGx/HvG8+i+d/IEwP+ry9WFOdxy/zJwwn8LWMM/2b3adRfb6Tr0OIC18\nd18LrK093mJmy4Eja7E1sW0i0iT79RnfzI4GeoBHa09dbGaDZnajmXU1uG0i0iR1F37tNP8O4FJ3\n3wJcDxzj7j0MnxF09Cm/iOxR17X6ZjaR4aK/2d3vAnD3kR/ObgB+Mtbr+/v2fA5c0DuDnl6dHIg0\n2uDAJpYObK5r2XoH6XwfeMbdr33nCTObXfv8D/BZ4KmxXry4b16dmxGR8erp7drroLrkqtVjLpsW\nvpmdCnwRWGZmvwIcuBxYZGY9wC5gJXBRpVaLSMvU863+QzBqn9G9jW+OiLSCxuMnsn7YyWk/b5xn\n93fPZP3ox7MizLck88pnZrEuzLN+9jc/GV9HkN03YAOHh3nWj5/J/v87tZ8+o0t2RQqkwhcpkApf\npEAqfJECqfBFCqTCFymQCl+kQOrHr6jd/bjZ9ruTfvYsr2orcT99lld1oPbDV6UjvkiBVPgiBVLh\nixSo5YU/OLCp1ZvcL2pfNZ3cvqUD1cZFNFsr913LC7/eiQLaRe2rppPbt7SDfylBa/edTvVFCqTC\nFymQuXtzN2DW3A2IyJjcfdQp8Jte+CLSeXSqL1IgFb5IgVT4IgVS4YsUSIUvUqD/B/OL1pQWhW2B\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa9092c4710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.matshow(img2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD7CAYAAABKWyniAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEptJREFUeJzt3W2MXdV1xvH/ih0b/DaAiccYO+DwmjbB05aSJtBmIqeU\nUiRIlJLUiWQSQEillChqFAIfGNpKDamEShWIIiDpGEFDIKIQFQFJ0CgCwksh4xgwBhJM7BiPjcE2\nxmBje/XDXOyxmVnres59w/v5SZbv3Ofcc7aP75pz7t1n72PujoiU5X3tboCItJ4KX6RAKnyRAqnw\nRQqkwhcpkApfpEAtK3wzO8PMnjWz58zsG63abr3MbKWZLTWzX5nZYx3QnpvMbMjMfj3iuUPN7H4z\nW2Fm95lZV4e170ozW21mT9b+nNHG9s01swfM7GkzW2Zm/1h7viP24Sjtu6T2fEv2obWiH9/M3gc8\nBywE1gCPA19w92ebvvE6mdlvgT9x99fa3RYAMzsN2AIscfeTas9dDWxw92/Xfnke6u6XdVD7rgRe\nd/dr2tGmkcxsNjDb3QfNbBrwBHA28GU6YB8G7fs8LdiHrTrinwI87+4vufvbwA8Z/kd2EqODPvq4\n+4PAvr+Ezgb6a4/7gXNa2qgRxmgfDO/HtnP3te4+WHu8BVgOzKVD9uEY7TuyFjd9H7bqjX4ksGrE\nz6vZ84/sFA7cZ2aPm9mF7W7MGGa5+xAMv3GAWW1uz2guNrNBM7uxnR9FRjKzo4Ee4BGgu9P24Yj2\nPVp7qun7sGOOcB3gVHc/GTiT4R1/WrsbVIdOu976euAYd+8B1gKdcMo/DbgDuLR2ZN13n7V1H47S\nvpbsw1YV/u+BD474eW7tuY7h7i/X/l4P3Mnwx5NOM2Rm3bD7M+K6NrdnL+6+3vd8aXQD8KftbI+Z\nTWS4qG5297tqT3fMPhytfa3ah60q/MeBY83sKDObBHwBuLtF206Z2ZTab17MbCpwOvBUe1sFDH/W\nG/l5727gvNrjxcBd+76gxfZqX62Q3vFZ2r8Pvw884+7Xjniuk/bhu9rXqn3Ykm/1Ybg7D7iW4V82\nN7n7t1qy4TqY2XyGj/IOTARuaXf7zOxWoBeYCQwBVwL/A9wOzANeAs51940d1L5PMfxZdRewErjo\nnc/TbWjfqcAvgGUM/786cDnwGPAj2rwPg/YtogX7sGWFLyKdQ1/uiRRIhS9SIBW+SIFU+CIFqlT4\nnT7wRkRGN+5v9esdeGNm6jYQaRN3H/W6/4kV1rl74A2Amb0z8OZdI+5+7h/f/bi/bxWL++ZV2Gxz\nqX3VdHL7Orlt0Pj2LbRfjplVOdV/Lwy8EZFR6Ms9kQJVOdWve+BNf9+eE4Oph0yosMnmW9A7o91N\nCKl949fJbYPq7Rsc2MTSgc11LVvly70JwAqGv9x7meFroP/O3Zfvs5yP/IwvIq2x0H7Z+C/33H2n\nmf0DcD97Bt4sT14mIh2gyqk+7n4vcEKD2iIiLaIv90QKpMIXKZAKX6RAKnyRAqnwRQqkwhcpkApf\npEAqfJECqfBFCqTCFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAlcbji3S6Cexs6vp30tlTyY1FR3yR\nAqnwRQqkwhcpkApfpEAqfJECqfBFCqTCFylQ8f34Vfths37iPN8R5hMrrj+zjUlhvjN5i2T7L8ur\n7r/MdF4P88lsq7T+rUwJ82z/bmdype2Pl474IgVS4YsUSIUvUiAVvkiBVPgiBVLhixRIhS9SoEr9\n+Ga2EtgE7ALedvdTGtGoRqraz5zJ+pmnsDXM57AmzKfvjPuhd0yI27+RQ8J8Hd1hnvUzb+XgMH8z\n6efekez/7DqGmWwI826GwnzWplfDPLOhqyvM1zAnzKteJzFeVS/g2QX0uvtrjWiMiLRG1VN9a8A6\nRKTFqhatA/eZ2eNmdmEjGiQizVf1VP9Ud3/ZzD4A/NTMlrv7g/su1N+3avfjBb0z6OmNPxeJyP4b\nHNjE0oHNdS1bqfDd/eXa3+vN7E7gFOBdhb+4b16VzYhIHXp6u/Y6qC65avWYy477VN/MppjZtNrj\nqcDpwFPjXZ+ItE6VI343cKeZeW09t7j7/Y1plog007gL391fBHoa2JZxqdpPn/UjZ/2sB/NmmM/k\nlTA/8tmkH/m3ccyst8N4+8lx+17g2DB/mE+E+WDyFniRo8N8SrL/PsajYb6IW8L8yHXJ/n0xjpka\nx9MPSsb7T47H+7/O9KQBzaGuOJECqfBFCqTCFymQCl+kQCp8kQKp8EUKpMIXKdABP69+1k+fjTev\nOh56/qtr4wXuS1awKsnPi+Pb+HyYX8ffh/nzKxbEG3ghjtP8oDjeeEE8n8BnJtwZr2B5sv14OgQ4\nLo53TozfH9vaNG9+Rkd8kQKp8EUKpMIXKZAKX6RAKnyRAqnwRQqkwhcp0AHfj5+Np696f/js/vY2\nGMaQTV3ysTj+zkfOD/OvPvG9eAXvmihtH7OTfFqSv5XkY88OBcDWLfG8/bO61sUryMbbx/99kEwP\n+cqEmWGe3XegWfPmZ3TEFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAKnyRAhXQjx/3k2b3b89e303S\nj/zzOE5u787mK94f5pf85oZ4Bd9Jtn9aHHedE88nMGny9jBfv/GD8QaSfv4/6oovhPiLVx+LV5Ds\nXz4ax68eF08YsI7uMM/me2gXHfFFCqTCFymQCl+kQCp8kQKp8EUKpMIXKZAKX6RAaT++md0EnAUM\nuftJtecOBW4DjgJWAue6+6YmtnNMWT97Nt5+a9KPP534/ufHrUoGlD8Sx/xZHH99wr/HC/yHxfmJ\ncfyh858O8z/nF2G+nD8I8/Vbkn78ZLz/57gjzO32+PWpj8TxSuaHedX7209gZ6XXj1c9R/wfAH+1\nz3OXAT9z9xOAB4BvNrphItI8aeG7+4PAa/s8fTbQX3vcD5zT4HaJSBON9zP+LHcfAnD3tcCsxjVJ\nRJqtUdfqexT29+25AdyC3hn09CYTmYnIfhsc2MTSgc11LTvewh8ys253HzKz2RCPVFncN2+cmxGR\nevX0du11UF1y1dhfPNd7qm+1P++4mz33aV0M3LVfLRSRtkoL38xuBR4Gjjez35nZl4FvAX9pZiuA\nhbWfReQ9Ij3Vd/dFY0SfbnBbmiIbD51dBzCTV+INZPPSZ/e3vzqObxy6IF4gub/8Yf/0+zBfxC1h\nPpl4vP09/E3cgJ/FMV+N4wt23hgvkO3/v43jp+d8KMzXcESygVi2/9pFV+6JFEiFL1IgFb5IgVT4\nIgVS4YsUSIUvUiAVvkiB3vPz6u9I+uG3J+Pxs/HQ81kZN+CJOGZqHD9y8oIw3/XjZAXJvPhnTrgn\nzLP7AjzMJ8J8/dXJePukH//Dtz8Z5jOueztewbY4fmthnP86mVg/uw7kEDbGG1A/voh0ChW+SIFU\n+CIFUuGLFEiFL1IgFb5IgVT4IgXq+H78bLx8VVPYGuZHsCZeQXb/9Z44vu9dM5fvY2Wy/mRe+Oc4\nPsyXJf3YS69LJv6/7MU4nxjPS7+IW+PX/18c88k4vm/q6WG+jJPC/PBkPoZJyYUEk7MLDSoa77z8\nOuKLFEiFL1IgFb5IgVT4IgVS4YsUSIUvUiAVvkiBOr4fv6qDeTPMs378Q199K95AMlw+6SbnmeT+\n8jySrD+Zdv6xiUlH91PJ+vmvJF8Zxxf0hfHnuCN+fbJ/3zovzrPrJDYwM8wnsCPMs/H403k9WX/c\nDz/efvqMjvgiBVLhixRIhS9SIBW+SIFU+CIFUuGLFEiFL1KgtB/fzG4CzgKG3P2k2nNXAhfC7knZ\nL3f3e5vRwKwfc2KSZ/30WT+sJd34dCV5PC07x7MiXiAZz591g8PzSR5f5wB/HcfHdofxh78bz5t/\n4t0vxet/I45XTD0uzLN++qrzPTSrn73Z6jni/wBGvQriGnf/49qfphS9iDRHWvju/iDw2iiRNb45\nItIKVT7jX2xmg2Z2o5llJ7wi0kHGe63+9cA/u7ub2b8C1wDnj7Vwf9+q3Y8X9M6gp1e/J0QabXBg\nE0sHNte17LgK393Xj/jxBuAn0fKL++aNZzMish96erv2OqguuWr1mMvWe6pvjPhMb2azR2SfpY4x\nXiLSOerpzrsV6AVmmtnvgCuBT5lZD7CL4XGZFzWxjSLSYGnhu/uiUZ7+QRPaMi5ZP2o2HjrLU9l4\n/GRe+H85+d/CfP4VK8P8nivODPNXODzMs+schog/M65L8m/z9TDnG3GcrJ7u3ZeSjG4eq8I8Mye5\nr0L2/smuM2kXXbknUiAVvkiBVPgiBVLhixRIhS9SIBW+SIFU+CIFOuDn1Z/E9iRP7l+e7aFZST6Y\n5BfE8Vcm/necT41z5iTbXxjH934pnpf/h3w+zM969oEw/+6z8fZPiWNO2rQpzOd0xf3wmcPZEObZ\ndRDZ+69ddMQXKZAKX6RAKnyRAqnwRQqkwhcpkApfpEAqfJECvef78avOa74z2QVrZ8XzA85eGPcj\nc1DSgEeS/MUkz6aF/2Acv/WZOL8xudDg13w0XkHcDc5QHHNccl+CdV2HhfnrTA/zqvMxHMjz6ovI\nAUaFL1IgFb5IgVT4IgVS4YsUSIUvUiAVvkiB3vP9+Jns/udbOTjM1yQD2jfOPyTMZ86PO7IPP3NL\nmNtbYcz6OdPCfJCeML8tGU//44e+FDdg7ttxnrzDZsQxM5LrEF6gO8y3MiXMs/HyO5L3T/b+6tR+\nfh3xRQqkwhcpkApfpEAqfJECqfBFCqTCFymQCl+kQGk/vpnNBZYA3cAu4AZ3/08zOxS4DTgKWAmc\n6+7J4PTOs514wHeWZ+O9h5J+5p2Hxf3AG4mvE1jJ0WG+ghPC/H85M8zT+QLOSd5CyTsivooCiKdD\nSPvRM9nrs///rJ9/Qof289dzxN8BfM3d/xD4OHCxmZ0IXAb8zN1PAB4Avtm8ZopII6WF7+5r3X2w\n9ngLsByYC5wN9NcW6wfOaVYjRaSx9uszvpkdDfQwfALY7e5DMPzLgfxmUiLSIeq+Vt/MpgF3AJe6\n+xYz830W2ffn3fr7Vu1+vKB3Bj29yQc3EdlvgwObWDqwua5l6yp8M5vIcNHf7O531Z4eMrNudx8y\ns9nAurFev7hvXl2NEZHx6+nt2uuguuSq1WMuW++p/veBZ9z92hHP3Q2cV3u8GLhr3xeJSGeqpzvv\nVOCLwDIz+xXDp/SXA1cDPzKzrwAvAec2s6Ei0jhp4bv7Q4w9e/unG9ucxqvaz5vZzqQw31bxOoA1\nHBHmLyfzBWTzCax9KfkYtjGOu+YmM+P/Mo7jfz0kuy+9P32WZ7L3T3ZfBpLx/u2iK/dECqTCFymQ\nCl+kQCp8kQKp8EUKpMIXKZAKX6RAB/y8+s2WjcfO5u3Pxttv5NBk/fG88RuYGeasfn+cJ+ZNXhUv\n8EYcz620dTgkudAgmze/VDriixRIhS9SIBW+SIFU+CIFUuGLFEiFL1IgFb5IgdSPX1E2Hjuftz8e\nz5+NB9+WvD67ToBX4phpcTyHNfECU+P42GTz2etnvvFqmB8y9bUwz+ZD6NT721elI75IgVT4IgVS\n4YsUSIUvUiAVvkiBVPgiBVLhixRI/fiJqvPyZ/3A2Xjx6bwe5mk/fSbpp+fEOM7al71+7unJ9k+O\n461TDwrz6WwJ8+w6DPXji8gBQ4UvUiAVvkiBVPgiBVLhixRIhS9SoLTwzWyumT1gZk+b2TIzu6T2\n/JVmttrMnqz9OaP5zRWRRqinH38H8DV3HzSzacATZvbTWnaNu1/TvOZV1+x+2AnsCPPJbAvzrB88\ne312ncEx/CbMty6M5+WflGx/HvG8+i+d/IEwP+ry9WFOdxy/zJwwn8LWMM/2b3adRfb6Tr0OIC18\nd18LrK093mJmy4Eja7E1sW0i0iT79RnfzI4GeoBHa09dbGaDZnajmXU1uG0i0iR1F37tNP8O4FJ3\n3wJcDxzj7j0MnxF09Cm/iOxR17X6ZjaR4aK/2d3vAnD3kR/ObgB+Mtbr+/v2fA5c0DuDnl6dHIg0\n2uDAJpYObK5r2XoH6XwfeMbdr33nCTObXfv8D/BZ4KmxXry4b16dmxGR8erp7drroLrkqtVjLpsW\nvpmdCnwRWGZmvwIcuBxYZGY9wC5gJXBRpVaLSMvU863+QzBqn9G9jW+OiLSCxuMnsn7YyWk/b5xn\n93fPZP3ox7MizLck88pnZrEuzLN+9jc/GV9HkN03YAOHh3nWj5/J/v87tZ8+o0t2RQqkwhcpkApf\npEAqfJECqfBFCqTCFymQCl+kQOrHr6jd/bjZ9ruTfvYsr2orcT99lld1oPbDV6UjvkiBVPgiBVLh\nixSo5YU/OLCp1ZvcL2pfNZ3cvqUD1cZFNFsr913LC7/eiQLaRe2rppPbt7SDfylBa/edTvVFCqTC\nFymQuXtzN2DW3A2IyJjcfdQp8Jte+CLSeXSqL1IgFb5IgVT4IgVS4YsUSIUvUqD/B/OL1pQWhW2B\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa909176b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.matshow(img3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}