ml1_labs/lab2/lab02_mul_gallegoposada_menyah.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 2: Classification\n",
    "\n",
    "### Machine Learning 1, September 2016\n",
    "\n",
    "* The lab exercises should be made in groups of two people.\n",
    "* The deadline is October 9th (Sunday) 23:59.\n",
    "* Assignment should be sent to your teaching assistant. The subject line of your email should be \"lab\\#\\_lastname1\\_lastname2\\_lastname3\".\n",
    "* Put your and your teammates' names in the body of the email.\n",
    "* Attach the .IPYNB (IPython Notebook) file containing your code and answers. Naming of the file follows the same rule as the subject line. For example, if the subject line is \"lab01\\_Kingma\\_Hu\", the attached file should be \"lab01\\_Kingma\\_Hu.ipynb\". Only use underscores (\"\\_\") to connect names, otherwise the files cannot be parsed.\n",
    "\n",
    "Notes on implementation:\n",
    "\n",
    "* For this notebook you need to answer a few theory questions, add them in the Markdown cell's below the question. Note: you can use Latex-style code in here.\n",
    "* Focus on Part 1 the first week, and Part 2 the second week!\n",
    "* You should write your code and answers below the questions in this IPython Notebook.\n",
    "* Among the first lines of your notebook should be \"%pylab inline\". This imports all required modules, and your plots will appear inline.\n",
    "* If you have questions outside of the labs, post them on blackboard or email me.\n",
    "* NOTE: Make sure we can run your notebook / scripts!\n",
    "\n",
    "$\\newcommand{\\bx}{\\mathbf{x}}$\n",
    "$\\newcommand{\\bw}{\\mathbf{w}}$\n",
    "$\\newcommand{\\bt}{\\mathbf{t}}$\n",
    "$\\newcommand{\\by}{\\mathbf{y}}$\n",
    "$\\newcommand{\\bm}{\\mathbf{m}}$\n",
    "$\\newcommand{\\bb}{\\mathbf{b}}$\n",
    "$\\newcommand{\\bS}{\\mathbf{S}}$\n",
    "$\\newcommand{\\ba}{\\mathbf{a}}$\n",
    "$\\newcommand{\\bz}{\\mathbf{z}}$\n",
    "$\\newcommand{\\bv}{\\mathbf{v}}$\n",
    "$\\newcommand{\\bq}{\\mathbf{q}}$\n",
    "$\\newcommand{\\bp}{\\mathbf{p}}$\n",
    "$\\newcommand{\\bh}{\\mathbf{h}}$\n",
    "$\\newcommand{\\bI}{\\mathbf{I}}$\n",
    "$\\newcommand{\\bX}{\\mathbf{X}}$\n",
    "$\\newcommand{\\bT}{\\mathbf{T}}$\n",
    "$\\newcommand{\\bPhi}{\\mathbf{\\Phi}}$\n",
    "$\\newcommand{\\bW}{\\mathbf{W}}$\n",
    "$\\newcommand{\\bV}{\\mathbf{V}}$\n",
    "$\\newcommand{\\bA}{\\mathbf{A}}$\n",
    "$\\newcommand{\\bj}{\\mathbf{j}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "%pylab inline\n",
    "import gzip, cPickle\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1. Multiclass logistic regression\n",
    "\n",
    "Scenario: you have a friend with one big problem: she's completely blind. You decided to help her: she has a special smartphone for blind people, and you are going to develop a mobile phone app that can do _machine vision_ using the mobile camera: converting a picture (from the camera) to the meaning of the image. You decide to start with an app that can read handwritten digits, i.e. convert an image of handwritten digits to text (e.g. it would enable her to read precious handwritten phone numbers).\n",
    "\n",
    "A key building block for such an app would be a function `predict_digit(x)` that returns the digit class of an image patch $\\bx$. Since hand-coding this function is highly non-trivial, you decide to solve this problem using machine learning, such that the internal parameters of this function are automatically learned using machine learning techniques.\n",
    "\n",
    "The dataset you're going to use for this is the MNIST handwritten digits dataset (`http://yann.lecun.com/exdb/mnist/`). You can load the data from `mnist.pkl.gz` we provided, using: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def load_mnist():\n",
    "    f = gzip.open('mnist.pkl.gz', 'rb')\n",
    "    data = cPickle.load(f)\n",
    "    f.close()\n",
    "    return data\n",
    "\n",
    "(x_train, t_train), (x_valid, t_valid), (x_test, t_test) = load_mnist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tuples represent train, validation and test sets. The first element (`x_train`, `x_valid`, `x_test`) of each tuple is a $N \\times M$ matrix, where $N$ is the number of datapoints and $M = 28^2 = 784$ is the dimensionality of the data. The second element (`t_train`, `t_valid`, `t_test`) of each tuple is the corresponding $N$-dimensional vector of integers, containing the true class labels.\n",
    "\n",
    "Here's a visualisation of the first 8 digits of the trainingset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAADiCAYAAACFkJabAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFupJREFUeJzt3Xu0lFUZx/EfgmBeMGV5w0yQNM1LhoaauhaBV1QUSSXD\n5OINy4RQ6Srq0aJcBhKKYiiYLTNNa2ku81Z4ByRBFAXTRMAI0KRQQBH7o7WfeYaZc86cOTNz5pn5\nfv7hWfucmdm85/Cw3733++x2n3zyiQAAMWzW1h0AABSOpA0AgZC0ASAQkjYABELSBoBASNoAEAhJ\nGwACIWkDQCAkbQAIpEMFPoNHLgvTrojXcG0LU8y1lbi+heDallfO9WWkDQCBkLQBIBCSNgAEQtIG\ngEBI2gAQCEkbAAIhaQNAICRtAAiEpA0AgVTiiUjUqSVLlkiSrr/+emsbP368xaNGjZIkXXzxxda2\n2267Vah3QEyMtAEgkHYVOI29TWsMbNy4UZK0fv36Jr9v+vTpFr///vsWL1iwQJI0YcIEa/vBD34g\nSZo0aZK1fepTn7L4uuuukySNGDGiJV2tidojy5Yts/iLX/yiJOm9995r8jXbbbedxStXrixHt6iP\nIemVV16x+KijjpIkzZ0719p22GGHYt62Lq/tLbfcIkm64IILrC3lmoULF1rbXnvt1dqPovYIAERG\n0gaAQEIuRK5evdrijz/+WJI0b948a3v44YctTrfmU6ZMKeqzunXrJkkaPXq0tU2dOlWStO2221rb\nkUceaXGfPn2K+qyoFi9ebHHv3r0t/ve//y1Jatcuc4fnr1mnTp0kSStWrLC2N954Q5K0++67W1v7\n9u1L2+ESe+211yxOf2dJ6tWrV1t0p1EzZ860uG/fvm3Yk5gee+wxi7/73e9KkjbbLHfc63/fy4GR\nNgAEQtIGgEDCTI8sXbrU4gMPPNBifztaKv6WJ02F+N0hw4cPlyTtuOOO1rb11ltbXOQqfAgfffSR\nxWla5LjjjrO2tDe7Mf5nd80110iSjjjiCGvbc889JWVPZ6XrXa38bfOrr75qcbVMj6QdYn4aZ9Gi\nRW3VnbD8NVu3bl2b9YORNgAEEmak3aVLF4t32mkni4sZaR9zzDE573vvvfdaW1ogk7IX1iBdeuml\nFvt96oWaMWOGxWk//IABA6wt/RxeeOGFYrtYcRMnTrTY/25VizVr1kiSfvrTn1pbegq1lu8KSyE9\npyFJV1xxRc7Xe/bsaXHaALHVVluVtU+MtAEgEJI2AAQSZnrELwROmzbN4nvuuUeSdNhhh1nbwIED\nc17vF7v++Mc/WtyxY0dJ0vLly63NFzhC9uLiHXfcYXG+Egh+qiP9HAYPHmxtviDUPvvsI0kaM2aM\ntaWfZwXKK5RMelagWvlHrZN07ZHf3//+d0lSv379rO3dd9/N+b5x48ZZ7J9BKCdG2gAQSJiRtvfl\nL3/Z4gMOOEBSZsQsSZdddpnFP//5zyVJDQ0N1ua/N9l5550t9gs29SwVf/rSl75kbb74U3ry6xvf\n+Ia1pUI6UmYRx7cNGjTI4i233FKS1LVrV2tL2y1//etfW9v3vvc9i6updOvbb78tKbtIVjXKN0I8\n+uij26AncfzqV7+S1PgW1lNPPVWS9NWvfrVifUoYaQNAICRtAAgk5PSI5/dUJ74+c+L30vriTuUu\n7hLNqlWrLP7Zz34mKXsvvN8j3717d0nZdcP91FN6+tE/BVmoDz74wOJrr73WYv9zbGtpX67va7Xw\nNeHnz5+f83X/3AP+L9/vnH862l8zP91aaYy0ASAQkjYABBJ+eiSfkSNHWjxr1ixJ0n333WdtL7/8\nssX77bdf5TpWpTZs2GDxJZdcYnHak+33n/75z3+2+HOf+5yk7CJS5fCPf/yjrO9frJdeeimnrZip\noHL44Q9/aHHa5ZJ2Wkn5d1DVI78b6uSTT27ye/1j7HvvvXe5utQsRtoAEEhNjrT9KCKV+PTlM/3/\nqKeccook6fDDD7c2/1RfPSxUvvXWWxb7Jx6T5557zuJ8B5X6p1Xr3SGHHFKRz/EHVc+ZM0dSdjnb\nu+66K+c1fhF3iy22KGPv4njyySctfuaZZ3K+ftppp1k8ZMiQSnSpWYy0ASAQkjYABFKT0yPe9ttv\nLyl7Ac2ftDJhwoSsPyXp1ltvtTgVPfIn09Sab33rWxb7Qk1pmijflEi5bNy4UVL2/thIxaP8wlZT\n0uKglPk7S5l6437x9cMPP5Qk/fKXv7Q2X6Qq1W/2tbz99EdaKKZIVMbs2bMlSWeffXber5900kmS\nskswVMuUEiNtAAiEpA0AgdT89EjiD1n1+7RHjRolSbr77rutbdiwYRa//vrrkrKP2dpmm23K1s9K\nSkd6PfHEE9bmd8v4lfNKSdMivh8HH3xwxftRiFSl0Pe1f//+Fn/+859v9LXPPvusxX76p0OH//+T\n9NNxaUeK30PvSzGkveH+mCtfDTE90l7vR4v5qatDDz20ye9NzyCU++iwYjDSBoBA2lVgkaeqV5HW\nrVsnKXsv8lFHHWVxuj5f+9rXrC3fHtgSKGZDeKuubRrt+VGbr22d6mGXaxE2PYnp9w+nOxo/yr/9\n9tstLvJJvmI32xd0fadPn27xX//61xZ/yJlnnmlxGuGlYlwt8eCDD1p84oknWpye3vOH1JZQWa9t\nKf34xz+22J84k09aKK6Cu5Oc68tIGwACIWkDQCB1sxDZmLT3snfv3tbWvn17i9Mt/B/+8AdrW7hw\noaSmF5qi8ntRyzEt4otTTZ48WVL28XDdunWTlF3wqNqLG/m9vo3t+62EBx54IG+7X1ivN/4ouHRo\ndGOGDh1qcRVMizSKkTYABELSBoBA6nJ6xD9CfO+990rK3jfrb+ETfwJ8JR/rrrSzzjqr5O/pb1HT\nEWaSdOONN0rKvi31jw2jNNLJ4fXI7/H3R+klxx57rMWTJk2qSJ9ai5E2AARS8yPtlStXSpJuuOEG\na7vtttssXrp0aZOvT4uSaYFMqp0a22kPut+rP23aNIv9vtZi3HnnnZKkiy66yNr8IcHf+c53JEnj\nx49v1ecAjVmxYoXFvghZMmbMGIurfcE7YaQNAIGQtAEgkJqZHlmzZo3F999/v8VXXXWVJGnRokUF\nv1efPn0sTo+7HnTQQa3tYtVJ0zx+usdPF6VrN3z4cGvzxbJS4a2bb77Z2vzxTW+++aYkqUePHtY2\naNAgi9P0CErPT3ktXrxYkrTHHnu0VXcqLhXX8rXK8/GHHUfBSBsAAiFpA0AgIadHUn1gSVqyZIkk\nafDgwdaW6kQXIh3RdOWVV1qb35NdKztFCuWPsUrTI1OnTrW2dHybJM2fP7/J9zr++OMlZR/v9u1v\nf7sk/UTT/O9tc1MEtSLfI+t+x0inTp0sHjt2rKTqrJfdHEbaABBI1Y+0165dK0kaOXKktT311FMW\nv/rqqwW9T79+/Sy+/PLLLU6nfmy++eat6mdE++67r6Ts+uGPPvpozvf5xUk/mkl23HFHi0eMGGFx\na/d5ozQef/xxSVLfvn3buCfl5Tcj5Ps99c9a+P3Z0TDSBoBASNoAEEhVTY+kfb0/+clPrC3drqe9\npoVIB65KUkNDgyTpwgsvtLYoj6uWW+fOnSVl1xn2R3s1t4/66quvliSde+651talS5dSdhFFqsAx\ngmgjjLQBIJCqGmn//ve/l5S9xSyfnj17Wvz1r39dktShQ+avct5551nsT2JBfv6EGn9H4mNUv4ED\nB1p80003tWFP2sauu+5q8QknnCAp++noWsFIGwACIWkDQCDtKrBgwYpIYYp59JJrW5hiH2vl+jaP\na1teOdeXkTYABELSBoBASNoAEAhJGwACIWkDQCAkbQAIhKQNAIFUYp82AKBEGGkDQCAkbQAIhKQN\nAIGQtAEgEJI2AARC0gaAQEjaABAISRsAAiFpA0AgJG0ACISkDQCBkLQBIBCSNgAEQtIGgEBI2gAQ\nCEkbAAIhaQNAICRtAAiEpA0AgZC0ASAQkjYABELSBoBASNoAEAhJGwACIWkDQCAkbQAIhKQNAIGQ\ntAEgEJI2AARC0gaAQEjaABAISRsAAiFpA0AgJG0ACISkDQCBkLQBIBCSNgAEQtIGgEBI2gAQCEkb\nAAIhaQNAICRtAAiEpA0AgZC0ASAQkjYABELSBoBASNoAEAhJGwAC6VCBz/ikAp9RC9oV8RqubWGK\nubYS17cQXNvyyrm+jLQBIBCSNgAEQtIGgEBI2gAQCEkbAAIhaQNAICRtAAiEpA0AgZC0ASCQSjwR\nCaCFGhoaLL788sslSb169bK2hx9+2OJtt922ch1Dm2OkDQCBtPvkk7KXAKDGQGFqrvbI+vXrJUkf\nffSRtT311FMWL1u2TJJ09tlnW1uHDmW5+QtTH+O9996TJO25557W9u6770qS2rXL/DVeeOEFi/ff\nf/8K9S6vMNd21apVFm/YsEGSNGvWLGs7+eSTLd5ss5aPZ4cOHSpJuvnmm62tffv2LX6fTVB7BAAi\nI2kDQCAsRKLV0i29JF133XUWP/7445KkmTNnNvn6NE0iZRbd6tWWW24pSerfv7+1TZs2rY16E9fy\n5cslSbfffru1TZkyxeKNGzdKkt566y1r81MifiqqUOnntN1221nb1VdfbXGnTp1a/J75MNIGgEBI\n2gAQSM3vHnnzzTclZd9iPvTQQxbPnj075zW/+c1vLN5tt90kSY888oi1DRkyRJLUrVu30nU0yO6R\nlStXWnz99ddn/SlJa9eutTj9bnXv3t3aunTpYvGcOXMkSTvttJO1zZ07V5K0ww47lLLbYXY4JP62\neuzYsZLYPdIS6d/oHXfcUfBrfC4sZnokn4ULF1rco0ePYt6C3SMAEFlNLkQ+/fTTFp9++umSpH/9\n61/W5v9HPfXUUyVJS5YssbbBgwfnvKd/TRpt3nDDDSXqcXVat26dpOxR3+TJky1evXp1k69PI8AZ\nM2ZYW9ofK2VG2P5nk96zxCPtMNI19yNptNxJJ50kqfGRdteuXSVJl1xyibWlxUkp/z7tJ5980uL7\n7ruvJP0sBiNtAAiEpA0AgYSfHkm3NGnBUZJOOOEEi9esWSNJOuWUU6zN3+6nx4U//vhjaxs2bJjF\nv/3tb3M+8ytf+Uorex1DmmYaN25cwa/5whe+YPETTzwhSercubO1vfPOOyXqXW1Kj/wvWLCgye97\n7rnnLP7sZz8ricJR3oABAyRlSgBsKk1/bL311gW/5/nnn2/xPvvsIyl7n3fi88fuu+9e8PsXipE2\nAAQSfqT9l7/8RZJ07LHH5v36GWecIUm69dZbrS3fk0m+kFG+0bXf3pf+F691zT2Jt9dee0mS+vTp\nY23XXHONxX6EnSxevLg0natR22yzjSRp1KhR1jZixIic7/NtaRtlWlRHZiSd73ewWH/7298s9sWn\nNpXufKTyFEBjpA0AgZC0ASCQkNMjEydOtDjdRvonmHzRoTFjxkhqvljLyJEjm/z6XXfdZXEq6lPr\nbrzxRknSYYcdZm3HHXecxWmf9VZbbVXwe65YsaJEvatt5513nsX5pkdQGX7a1D/5+8EHHzT6mksv\nvbSsfWKkDQCBkLQBIJAw0yM33XSTxX5lPU17DBo0yNq+//3vW7z55pvnvJd/lHrevHmSpNdee83a\n/CPraSrm4IMPLrrvUaWdDBdeeGHJ3jPV2Ebh0rMIxRyBhcKl5wokafTo0ZKkl19+2do+/PDDJl9/\n5JFHSir/z4nfAgAIpOpH2qmATkNDg7X5Rcc0wvb7sPPxT0alvdtSZp+35598Ovfcc1vY4/pzzz33\nWPyf//zH4nTH4n9eqRyr559g3WOPPcrRxdDSyK1U5ULrQTpN6Xe/+521Pfjgg02+5v7777e4uWv9\n6U9/WlL2yThHHHGEpPx396XESBsAAiFpA0AgVT89kgo5+ZrL3vjx4yVJ77//vrX52/W0v/rZZ5+1\nNn8Ln26D/O3QOeecY3HHjh2L7nutSEWMJOntt9+2OO2Hb6xmcXMLaOlUoNtuu83aWGxDsf75z39a\n3Lt3b0nS66+/XpbPSvW6+/XrV5b3bwr/QgAgEJI2AARS9dMj7du3lyTtvPPO1rZ8+XKLt99+e0nN\nr/b6yltp5VfKHDPmD5ft2bNnK3ocm68rvnTpUkmZW00p+1i29Dh/muaQpOOPP97iO++8U1Kmpvmm\n0n75P/3pT9Z25plnSsr83IFipJ1LLTm4vLnjxry0a+Tiiy+2tgMPPLAlXSwaI20ACKTqR9pbbLGF\npOzCLYceeqjF6ZBdf2LKWWedZfE3v/lNSdlFjfzX08ixnovy+NH13LlzLT7kkENyvjcVkZKkvn37\nSpJ69OhhbWvXrrX4xRdflCTNnDkz7+emO6ahQ4daW9qn7T+7HDWJI2luQfeRRx6RRD3tXXbZxeLZ\ns2dLku6++25rO+aYYywudIPB1KlTLR47dmxru1gSjLQBIBCSNgAE0q4lE/VFKvsHFMIXhErHZEmZ\nW07/uOvAgQMr17GMYp5RbtW1TdMivk7wZZddlvN9aXFQkqZMmWJxmrrytYVPPPFEi2fMmCEpu5b5\ntddea3GaivH7tJPTTz/dYl8fPd9BrJ/5zGdy2jZR7PPfVfG7mxZlm1tsX7ZsmcV+Yb3MQl/b5qQy\nGlL+373nn3/e4jItROZcX0baABAISRsAAqmbZXl/m+NX4dMtp99fXMv8XtQJEyZIyhzJJmVqaEuZ\n09j9SfdpSkTKnKzuKyH6msT777+/pOzT7ffee2+L169fL0m66KKLrC1Va5w+fbq1+amrxFcDXLRo\nUc7Xa8mPfvQjSdkn3edzyy235LwGreNPYK8WjLQBIJC6GWmnUV+9e+CBByxOI2y/wOJrCh900EGS\npIULF1qbP0EoFYrye7MnTZpkcVrA7Ny5c96+pAXKAw44wNrS6N8vBvsRZJIKhdUDf33wf2kRff78\n+da27777WtzamtZp7/tpp53WqvcpB0baABAISRsAAqmbfdr+Nsrvp0wLkb7GdiqEVGEV2aft9zSn\nx8j94mKaEpGk1atXS5JeeumlJt9z8uTJFg8fPtziKqqNXRN7if0U34IFC3K+7heZ33nnHYtTUbUy\nqdi19c9aXHHFFZIy9fKl7CMFG5uS25Sf2ps1a5bFqSRA+jewqZQj/Gv8InsJsU8bACIjaQNAIHWz\ne+SNN95o6y5UhW7dulmcpkf8Hvann3465zWDBw+2+Oijj7Y47W339cmraEqk5vTq1cviV155Jefr\ntX7thwwZYnG+ypF+R1Gh0yN+t1QquyDlLxngqyiOHj1aUtmmRJpU2z9lAKgxdbMQ6Q/97Nq1q8Vp\ndPLf//7X2mp5ITI9hShlDjv2o2tfk/iMM86QlL1QGfREmZpYiJw3b57FfsE48f+WU515qXYWIg8/\n/HCLG6vR3hr++u26666SsmvvX3nllRZXsMY7C5EAEBlJGwACqZvpEc/vd00LOn4PaPfu3SveJ7VB\nPe06UhPTI37PsD86a86cOZJqf3okHTQtSRMnTpQk/eIXvyjqw9PxhH7B0l/TVATNTxe2EaZHACCy\nuhxpP/bYYxansqMDBgywNl/0qMpPAKm6a1ulamKkXaXa5Npu2LBBkvTQQw9Z2znnnGPxqlWrJEnD\nhg2ztv79+1vcu3dvSflPo6kyjLQBIDKSNgAEUpfTI36v8tChQyVln47iT2JJh9527Nix3N1ieqR8\nmB4pH65teTE9AgCRkbQBIJC6nB7x0lTJuHHjrK2hocHiZcuWSarILhKmR8qHW/jy4dqWF9MjABBZ\n3Y+0qwgj7fJhNFg+XNvyYqQNAJGRtAEgkEpMjwAASoSRNgAEQtIGgEBI2gAQCEkbAAIhaQNAICRt\nAAiEpA0AgZC0ASAQkjYABELSBoBASNoAEAhJGwACIWkDQCAkbQAIhKQNAIGQtAEgEJI2AARC0gaA\nQEjaABAISRsAAvkfRU1fuKzaqHIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a35e03ad0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_digits(data, numcols, shape=(28,28)):\n",
    "    numdigits = data.shape[0]\n",
    "    numrows = int(numdigits/numcols)\n",
    "    for i in range(numdigits):\n",
    "        plt.subplot(numrows, numcols, i+1)\n",
    "        plt.axis('off')\n",
    "        plt.imshow(data[i].reshape(shape), interpolation='nearest', cmap='Greys')\n",
    "    plt.show()\n",
    "    \n",
    "plot_digits(x_train[0:8], numcols=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In _multiclass_ logistic regression, the conditional probability of class label $j$ given the image $\\bx$ for some datapoint is given by:\n",
    "\n",
    "$ \\log p(t = j \\;|\\; \\bx, \\bb, \\bW) = \\log q_j - \\log Z$\n",
    "\n",
    "where $\\log q_j = \\bw_j^T \\bx + b_j$ (the log of the unnormalized probability of the class $j$), and $Z = \\sum_k q_k$ is the normalizing factor. $\\bw_j$ is the $j$-th column of $\\bW$ (a matrix of size $784 \\times 10$) corresponding to the class label, $b_j$ is the $j$-th element of $\\bb$.\n",
    "\n",
    "Given an input image, the multiclass logistic regression model first computes the intermediate vector $\\log \\bq$ (of size $10 \\times 1$), using $\\log q_j = \\bw_j^T \\bx + b_j$, containing the unnormalized log-probabilities per class. \n",
    "\n",
    "The unnormalized probabilities are then normalized by $Z$ such that $\\sum_j p_j = \\sum_j \\exp(\\log p_j) = 1$. This is done by $\\log p_j = \\log q_j - \\log Z$ where $Z = \\sum_j \\exp(\\log q_j)$. This is known as the _softmax_ transformation, and is also used as a last layer of many classifcation neural network models, to ensure that the output of the network is a normalized distribution, regardless of the values of second-to-last layer ($\\log \\bq$)\n",
    "\n",
    "Warning: when computing $\\log Z$, you are likely to encounter numerical problems. Save yourself countless hours of debugging and learn the [log-sum-exp trick](https://hips.seas.harvard.edu/blog/2013/01/09/computing-log-sum-exp/ \"Title\").\n",
    "\n",
    "The network's output $\\log \\bp$ of size $10 \\times 1$ then contains the conditional log-probabilities $\\log p(t = j \\;|\\; \\bx, \\bb, \\bW)$ for each digit class $j$. In summary, the computations are done in this order:\n",
    "\n",
    "$\\bx \\rightarrow \\log \\bq \\rightarrow Z \\rightarrow \\log \\bp$\n",
    "\n",
    "Given some dataset with $N$ independent, identically distributed datapoints, the log-likelihood is given by:\n",
    "\n",
    "$ \\mathcal{L}(\\bb, \\bW) = \\sum_{n=1}^N \\mathcal{L}^{(n)}$\n",
    "\n",
    "where we use $\\mathcal{L}^{(n)}$ to denote the partial log-likelihood evaluated over a single datapoint. It is important to see that the log-probability of the class label $t^{(n)}$ given the image, is given by the $t^{(n)}$-th element of the network's output $\\log \\bp$, denoted by $\\log p_{t^{(n)}}$:\n",
    "\n",
    "$\\mathcal{L}^{(n)} = \\log p(t = t^{(n)} \\;|\\; \\bx = \\bx^{(n)}, \\bb, \\bW) = \\log p_{t^{(n)}} = \\log q_{t^{(n)}} - \\log Z^{(n)}$\n",
    "\n",
    "where $\\bx^{(n)}$ and $t^{(n)}$ are the input (image) and class label (integer) of the $n$-th datapoint, and $Z^{(n)}$ is the normalizing constant for the distribution over $t^{(n)}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Gradient-based stochastic optimization\n",
    "\n",
    "### 1.1.1 Derive gradient equations (20 points)\n",
    "\n",
    "Derive the equations for computing the (first) partial derivatives of the log-likelihood w.r.t. all the parameters, evaluated at a _single_ datapoint $n$.\n",
    "\n",
    "You should start deriving the equations for $\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log q_j}$ for each $j$. For clarity, we'll use the shorthand $\\delta^q_j = \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log q_j}$.\n",
    "\n",
    "For $j = t^{(n)}$:\n",
    "$\n",
    "\\delta^q_j\n",
    "= \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log p_j}\n",
    "\\frac{\\partial \\log p_j}{\\partial \\log q_j}\n",
    "+ \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log Z}\n",
    "\\frac{\\partial \\log Z}{\\partial Z} \n",
    "\\frac{\\partial Z}{\\partial \\log q_j} \n",
    "= 1 \\cdot 1 - \\frac{\\partial \\log Z}{\\partial Z} \n",
    "\\frac{\\partial Z}{\\partial \\log q_j}\n",
    "= 1 - \\frac{\\partial \\log Z}{\\partial Z} \n",
    "\\frac{\\partial Z}{\\partial \\log q_j}\n",
    "$\n",
    "\n",
    "For $j \\neq t^{(n)}$:\n",
    "$\n",
    "\\delta^q_j\n",
    "= \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log Z}\n",
    "\\frac{\\partial \\log Z}{\\partial Z} \n",
    "\\frac{\\partial Z}{\\partial \\log q_j} \n",
    "= - \\frac{\\partial \\log Z}{\\partial Z} \n",
    "\\frac{\\partial Z}{\\partial \\log q_j}\n",
    "$\n",
    "\n",
    "Complete the above derivations for $\\delta^q_j$ by furtherly developing $\\frac{\\partial \\log Z}{\\partial Z}$ and $\\frac{\\partial Z}{\\partial \\log q_j}$. Both are quite simple. For these it doesn't matter whether $j = t^{(n)}$ or not.\n",
    "\n",
    "Given your equations for computing the gradients $\\delta^q_j$ it should be quite straightforward to derive the equations for the gradients of the parameters of the model, $\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial W_{ij}}$ and $\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial b_j}$. The gradients for the biases $\\bb$ are given by:\n",
    "\n",
    "$\n",
    "\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial b_j}\n",
    "= \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log q_j}\n",
    "\\frac{\\partial \\log q_j}{\\partial b_j}\n",
    "= \\delta^q_j\n",
    "\\cdot 1\n",
    "= \\delta^q_j\n",
    "$\n",
    "\n",
    "The equation above gives the derivative of $\\mathcal{L}^{(n)}$ w.r.t. a single element of $\\bb$, so the vector $\\nabla_\\bb \\mathcal{L}^{(n)}$ with all derivatives of $\\mathcal{L}^{(n)}$ w.r.t. the bias parameters $\\bb$ is: \n",
    "\n",
    "$\n",
    "\\nabla_\\bb \\mathcal{L}^{(n)} = \\mathbf{\\delta}^q\n",
    "$\n",
    "\n",
    "where $\\mathbf{\\delta}^q$ denotes the vector of size $10 \\times 1$ with elements $\\mathbf{\\delta}_j^q$.\n",
    "\n",
    "The (not fully developed) equation for computing the derivative of $\\mathcal{L}^{(n)}$ w.r.t. a single element $W_{ij}$ of $\\bW$ is:\n",
    "\n",
    "$\n",
    "\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial W_{ij}} =\n",
    "\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log q_j}\n",
    "\\frac{\\partial \\log q_j}{\\partial W_{ij}}\n",
    "= \\mathbf{\\delta}_j^q\n",
    "\\frac{\\partial \\log q_j}{\\partial W_{ij}}\n",
    "$\n",
    "\n",
    "What is $\\frac{\\partial \\log q_j}{\\partial W_{ij}}$? Complete the equation above.\n",
    "\n",
    "If you want, you can give the resulting equation in vector format ($\\nabla_{\\bw_j} \\mathcal{L}^{(n)} = ...$), like we did for $\\nabla_\\bb \\mathcal{L}^{(n)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer:\n",
    "\n",
    "---\n",
    "\n",
    "First, $\\frac{\\partial \\log Z}{\\partial Z} = \\frac{1}{Z}$.\n",
    "\n",
    "Recall $Z = \\sum_j \\exp(\\log q_j)$. Then, $\\frac{\\partial Z}{\\partial \\log q_j} = \\exp(\\log q_j) = q_j$.\n",
    "\n",
    "Finally, $\\delta^q_j = \\mathbb{I}_{j,t^{(n)}} - \\frac{\\partial \\log Z}{\\partial Z} \n",
    "\\frac{\\partial Z}{\\partial \\log q_j} = \\mathbb{I}_{j,t^{(n)}} - \\frac{q_j}{Z}$. Where $\\mathbb{I}$ is the identity matrix.\n",
    "\n",
    "---\n",
    "\n",
    "$\\frac{\\partial \\log q_j}{\\partial \\bW_j} = \\frac{\\partial \\bw_j^T \\bx +\\bb_j}{\\partial \\bW_j} = \\bx \\rightarrow \\frac{\\partial \\log q_j}{\\partial \\bW_{ij}} = \\bx_i$\n",
    "\n",
    "$\\nabla_{\\bw_j} \\mathcal{L}^{(n)} = \\mathbf{\\delta}_j^q \\bx \\rightarrow \\nabla_{\\bW} \\mathcal{L}^{(n)} = \\bx {\\mathbf{\\delta}^q}^T$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1.2 Implement gradient computations (10 points)\n",
    "\n",
    "Implement the gradient calculations you derived in the previous question. Write a function `logreg_gradient(x, t, w, b)` that returns the gradients $\\nabla_{\\bw_j} \\mathcal{L}^{(n)}$ (for each $j$) and $\\nabla_{\\bb} \\mathcal{L}^{(n)}$, i.e. the first partial derivatives of the log-likelihood w.r.t. the parameters $\\bW$ and $\\bb$, evaluated at a single datapoint (`x`, `t`).\n",
    "The computation will contain roughly the following intermediate variables:\n",
    "\n",
    "$\n",
    "\\log \\bq \\rightarrow Z \\rightarrow \\log \\bp\\,,\\, \\mathbf{\\delta}^q\n",
    "$\n",
    "\n",
    "followed by computation of the gradient vectors $\\nabla_{\\bw_j} \\mathcal{L}^{(n)}$ (contained in a $784 \\times 10$ matrix) and $\\nabla_{\\bb} \\mathcal{L}^{(n)}$ (a $10 \\times 1$ vector).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def logreg_gradient(x,t,W,b):\n",
    "    lq = dot(W.T,x) + b\n",
    "    Z = sum(np.exp(lq))\n",
    "    \n",
    "    delta_q = - exp(lq)/Z\n",
    "    delta_q[t] += 1\n",
    "    \n",
    "    grad_b = delta_q\n",
    "    grad_W = dot(x,delta_q.T)\n",
    "    \n",
    "    return grad_W, grad_b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### 1.1.3 Stochastic gradient descent (10 points)\n",
    "\n",
    "Write a function `sgd_iter(x_train, t_train, w, b)` that performs one iteration of stochastic gradient descent (SGD), and returns the new weights. It should go through the trainingset once in randomized order, call `logreg_gradient(x, t, w, b)` for each datapoint to get the gradients, and update the parameters using a small learning rate (e.g. `1E-4`). Note that in this case we're maximizing the likelihood function, so we should actually performing gradient ___ascent___... For more information about SGD, see Bishop 5.2.4 or an online source (i.e. https://en.wikipedia.org/wiki/Stochastic_gradient_descent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def sgd_iter(x_train, t_train, W, b):\n",
    "    N, M = shape(x_train)\n",
    "    perm = np.random.permutation(N)\n",
    "    eta = 0.0001\n",
    "    for ix in perm:\n",
    "        x = reshape(x_train[ix],(M,-1))\n",
    "        t = t_train[ix]\n",
    "        grad_W, grad_b = logreg_gradient(x,t,W,b)\n",
    "        W += eta * grad_W\n",
    "        b += eta * grad_b\n",
    "    return W, b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2. Train\n",
    "\n",
    "### 1.2.1 Train (10 points)\n",
    "Perform a handful of training iterations through the trainingset. Plot (in one graph) the conditional log-probability of the trainingset and validation set after each iteration.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJ5ONJYSwhT2ACoKIMONaXOJSrdTt8au7\npXa5tr/eLtbb3l7p7+HP5fZe6229vXrb/lqrVYpdXHpVWqQgV2OrVoWERUXltiIgqyhbBLJMPr8/\nZhKSkMlMSGbOTPJ+Ph7nMWe+58z3fDJi3jnnexZzd0RERDqSF3QBIiKSvRQSIiKSkEJCREQSUkiI\niEhCCgkREUlIISEiIgmlNSTM7AEz225ma1q1lZnZUjN728yWmFlpq2XzzOx/zOxNMzs/nbWJiEhy\n6d6TeBC4oF3bzcAyd58CPAvMAzCzacCVwFTgQuAnZmZprk9ERDqR1pBw9xeAXe2aLwXmx+fnA5fF\n5y8Bfuvuje7+LvA/wMnprE9ERDoXxJjECHffDuDu24AR8fYxwKZW622Ot4mISECyYeBa9wUREclS\n+QFsc7uZlbv7djMbCeyIt28GxrVab2y87TBmpmARETkC7t6lsd5M7ElYfGq2EPhsfP564KlW7Veb\nWaGZTQSOBl5N1Km7a3Ln1ltvDbyGbJn0Xei70HfR+XQk0ronYWa/BiqBoWa2EbgV+B7wmJl9HthA\n7Iwm3H2tmT0KrAUagL/3I/2pRESkR6Q1JNz92gSLzkuw/p3AnemrSEREuiIbBq6lGyorK4MuIWvo\nuzhE38Uh+i66x3LxiI6Z6UiUiEgXmRnexYHrIM5uEpEcN2HCBDZs2BB0GZJARUUF7777bo/0pT0J\nEemy+F+kQZchCST673MkexIakxARkYQUEiIikpBCQkREElJIiIh0oqmpiZKSEt57770eXTdXaOBa\nRLosmweuS0pKaH4UzUcffURRURGhUAgz42c/+xnXXHNNwBWmX08OXCskRKTLsjkkWps0aRIPPPAA\nZ599dsJ1otEooVAog1Wln85uEhFJQUc3trvlllu4+uqrufbaayktLeVXv/oVL7/8MqeddhplZWWM\nGTOGG2+8kWg0CsRCJC8vj40bNwIwd+5cbrzxRubMmcOgQYOYPXt2yzUjXVkXYPHixUyZMoWysjK+\n/vWvc/rpp/PLX/4yE19NyhQSItLnPPnkk3z6059mz549XHXVVRQUFHDvvffy4Ycf8uKLL7JkyRJ+\n9rOftazf/knKv/nNb/iXf/kXdu3axbhx47jlllu6vO6OHTu46qqruPvuu9m5cycTJ05k+fLlafyp\nj4xCQkTSwqxnpnQ4/fTTmTNnDgBFRUVEIhFOOukkzIwJEyZwww038Pzzz7es335v5PLLL2fWrFmE\nQiGuu+46Vq1a1eV1Fy1axKxZs7jooosIhULcdNNNDB06ND0/cDfothwikhbZPGQxbty4Nu/ffvtt\nvvnNb1JdXc3+/fuJRqOccsopCT8/cuTIlvn+/ftTW1vb5XW3bNlyWB1jx47t0s+RCdqTEJE+p/0h\noS996Uscf/zxvPPOO+zZs4fbb7897QPzo0aNYtOmTW3aNm/u8GGcgVJIiEift2/fPkpLS+nXrx9v\nvvlmm/GIdLnoootYuXIlixYtIhqN8h//8R/s3Lkz7dvtKoWEiPRa7fcYErn77rt56KGHGDRoEF/+\n8pe5+uqrE/aTrM9U1x0xYgSPPPIIN910E8OGDWP9+vXMmjWLoqKilGrOFF0nISJdlivXSeSSpqYm\nRo8eze9+9ztmz57drb50nQSwt7Yh6BJERLplyZIl7Nmzh7q6Ou644w4KCws5+eSTgy6rjZwNiYV/\nWRt0CSIi3fLCCy8wadIkysvLeeaZZ3jyyScpKCgIuqw2cvZw03U/eICHv/n5oEsR6ZN0uCm76XAT\nUL21OugSRER6vZwNiQ31NUGXICLS6+VsSBwoWcPe2sagyxAR6dVyNiSK6sey8KW3gi5DRKRXy9mQ\nGBcK88fVGpcQEUmnnA2JmeURqrdqXEJEes6GDRvIy8ujqakJgDlz5rBgwYKU1u2qO++8ky9+8YtH\nXGum5GxInD89woZ67UmISFsXXnght91222HtTz31FKNGjUr6S731rTSefvpp5s6dm9K6nXn++ecP\nu+PrvHnzuO+++1L6fJByNiQuO3UWB0pXs2dvNOhSRCSLXH/99Tz88MOHtT/88MPMnTuXvLzM/9pz\n95QDJdvkbEgMLxlMUUM5C19aF3QpIpJFLrvsMj744ANeeOGFlrbdu3fzhz/8gblz5/L0008TDocp\nLS2loqKC22+/PWFfZ599Nr/4xS+A2L2VvvWtbzF8+HCOPvpoFi1a1Gbdhx56iGnTpjFo0CCOPvro\nlr2E/fv3M2fOHLZs2UJJSQmDBg1i27Zt3H777W32UhYuXMj06dMZMmQI55xzDm+9dejEnIkTJ3L3\n3XdzwgknUFZWxjXXXEN9fX2PfF/J5GxIAIzV4LWItFNcXMwVV1zR5lnRjzzyCFOnTuX4449nwIAB\nLFiwgD179rBo0SJ++tOfsnDhwqT93nfffTz99NOsXr2aFStW8Pjjj7dZXl5eztNPP83evXt58MEH\nuemmm1i1ahX9+/dn8eLFjB49mn379rF3796WBxE1712sW7eOa6+9lnvvvZf333+fCy+8kIsvvpjG\nxkOn+T/22GMsXbqU9evXs3r1ah566KEe+LaSy+kn080aGaF6fQ3w6aBLEZF27PaeObzit3b99h/X\nX389F110ET/60Y8oLCxkwYIFXH/99QCcddZZLetNnz6dq6++mueff55LLrmk0z4fe+wxvvGNbzB6\n9GggNqbQ+hGnF154Ycv8GWecwfnnn8+f//xnZs6cmbTeRx99lIsuuohzzjkHgG9961vcc889vPTS\nS5x55pkA3HjjjZSXlwNw8cUXt3lkajrldEh8/Lgwi9Z9N+gyRKQDR/LLvafMnj2b4cOH8+STT3Li\niSeyfPlynnjiCQBeeeUV5s2bx+uvv059fT319fVcccUVSfts/7jRioqKNssXL17MHXfcwbp162hq\nauLAgQPMmDEjpXq3bNnSpj8zY9y4cW2eVNccEBB7DOrWrVtT6ru7cvpw06UnhzlQuorde47sFDQR\n6b3mzp3L/Pnzefjhh7ngggsYPnw4ANdddx2XXXYZmzdvZvfu3XzpS19K6WaF7R83umHDhpb5+vp6\nLr/8cr797W/z/vvvs2vXLi688MKWfpMNWo8ePbpNfwCbNm3Kimde53RIlA8aSlHTEH7/4l+DLkVE\nssxnPvMZli1bxv33399yqAmgtraWsrIyCgoKePXVV/n1r3/d5nOJAuPKK6/k3nvvZfPmzezatYu7\n7rqrZVnzHsmwYcPIy8tj8eLFLF26tGV5eXk5H3zwAXv37k3Y96JFi3juuedobGzkBz/4AcXFxZx2\n2mnd+Qp6RE6HBMQGr5es0UV1ItJWRUUFH/vYx9i/f3+b8Yaf/OQn3HLLLZSWlvLd736Xq666qs3n\nEj1+9IYbbuCCCy7ghBNO4MQTT+RTn/pUy7KBAwdy7733csUVVzBkyBB++9vfcumll7YsnzJlCtdc\ncw2TJk1iyJAhbNu2rc02J0+ezMMPP8xXv/pVhg8fzqJFi/j9739Pfn7+YXVkWmDPkzCzm4AvAE3A\na8DngAHAI0AF8C5wpbvv6eCzLY8vveLef+W1v+7irXu/n6HKRUTPk8huOf88CTMbDXwNCLv7DGID\n6NcANwPL3H0K8CwwL1lf500Ps7FBp8GKiKRDkIebQsAAM8sH+gGbgUuB+fHl84HLknVyyYlhDgyu\n4cMP9VeNiEhPCyQk3H0LcDewkVg47HH3ZUC5u2+Pr7MNGJGsr1GDRlBICYteeiedJYuI9EmBXCdh\nZoOJ7TVUAHuAx8zsOqD97kDC3YPWN/Aa9l4Ff1xTw9yLjur5YkVEclRVVRVVVVXd6iOQgWszuxy4\nwN1viL+fC5wKnANUuvt2MxsJPOfuUzv4vLeu+/L/vIPX397PWz/6XmZ+AJE+TgPX2S3nB66JHWY6\n1cyKLXZu17nAWmAh8Nn4OtcDT6XS2cenR/TMaxGRNAjkcJO7v2pmjwMrgYb4631ACfComX0e2ABc\nmUp/F0cifHlJNTt3OsOG5ebteEVySUVFRc7e+rovaH/LkO4I7DqJ7mh/uAmg6Dujuf+0vzD34p77\nckREepNcOtzU48aGIrptuIhID+s1ITGrPMyKLQoJEZGe1GtC4uPTI2zU4LWISI/qNSExZ1aYuqHV\nbN+ee2MsIiLZqteExNhBY8jPN5a+vDn5yiIikpJeExJmxrhQhMUavBYR6TG9JiQAZpaHqd6icQkR\nkZ7Sq0Li/OkRNtRrT0JEpKf0qpC4YEaY+mHVZOj54CIivV6vComKwePJL2zgmZeVEiIiPaFXhYSZ\nMTYU1pXXIiI9pFeFBMCs8ogGr0VEekivC4nzpofZUF9NDt63UEQk6/S6kDh/eoSG4TVs1jV1IiLd\n1utCYlLZRPKKa/nvl3cEXYqISM7rdSHRPHi9ZI3GJUREuqvXhQTEBq9XbNYZTiIi3dUrQ+Lj08Ns\naNDgtYhId/XKkDhvWoToiBo2bgy6EhGR3NYrQ+KoIUdh/Xfx3MsfBF2KiEhO65UhkWd5jAnN0uC1\niEg39cqQgNhtw/XMaxGR7um1IfHx4yJsqK/R4LWISDf02pA4d2qEppHVrF8fdCUiIrmr14bEMUOO\nwQbsoOqVXUGXIiKSs3ptSITyQowJzWTpmpVBlyIikrN6bUiABq9FRLqrV4fEecdF2FhfQ1NT0JWI\niOSmXh0SlVPC+Khq/va3oCsREclNvTokjh12LD5wC396ZW/QpYiI5KReHRL5efmMyZ/BEg1ei4gc\nkV4dEhAbvK7RM69FRI5Irw+J846LsKGhWoPXIiJHoNeHxJlHh7HR1axbF3QlIiK5p9eHxLTh02ga\ntIEXXq0NuhQRkZwTWEiYWamZPWZmb5rZG2Z2ipmVmdlSM3vbzJaYWWl3t1MQKmB0/nSWrF7VE2WL\niPQpQe5J3AM87e5TgROAt4CbgWXuPgV4FpjXExuaWR6hZqsGr0VEuiqQkDCzQcAZ7v4ggLs3uvse\n4FJgfny1+cBlPbG986aF2dhYTTTaE72JiPQdQe1JTAR2mtmDZlZjZveZWX+g3N23A7j7NmBET2zs\n9KMi5I2p4a23eqI3EZG+I6iQyAfCwI/dPQx8ROxQU/tHBPXII4Omj5hOdNDfePHV/T3RnYhIn5Ef\n0HbfAza5+4r4+98RC4ntZlbu7tvNbCSwI1EHt912W8t8ZWUllZWVCTdWGCpkVP5Unlmzhi9yag+U\nLyKS/aqqqqiqqupWH+ZdeL6nmQ0ADrp7t4/um9nzwA3uvs7MbgX6xxd96O53mdk/AWXufnMHn/Wu\n1A1w8X1f5PX/PoH1j3ylu6WLiOQkM8PdrSuf6XRPwszygKuB64CTgDqgyMx2AouAn7n7X4+w3q8D\nvzKzAuAd4HNACHjUzD4PbACuPMK+D3PutDB/XPoqjY2QH9T+k4hIjkn26/I5YBmxU1Ffd/cmADMb\nApwN3GVmT7j7w13dsLuvJhY87Z3X1b5SMXtihNDYn7J2LcyYkY4tiIj0Pp0ebjKzAndv6LSDFNbp\naUdyuOlg40EG/vMQfjz+Q770heI0VSYikr2O5HBTp2c3tf/lb2bFZvZ3ZvY1Mxva0TrZqji/mPL8\nyTyz+rWgSxERyRldPQX2HqAe2AU82fPlpNeskRE981pEpAs6DQkz+42ZHdWqaQjwGLFTVsvSWVg6\nnDM1zHtNNdTXB12JiEhuSLYn8X+Afzazu81sMPAD4AlgMXBbmmvrcR+bEKFgfDVvvBF0JSIiuaHT\ns5vc/R3gWjM7HXiE2Gmvn+yJ6ySCMKN8Bg2lb/KX5XXMmlUUdDkiIlkv2eGmMjP7CjANuILYWMQS\nM7s4E8X1tP4F/RlRcBT/vUa7EiIiqUh2uOlJYDexeygtcPcFwMXALDP7fbqLS4eZ5WFWbNbgtYhI\nKpKFxFDgcWKD1WMA3P2Au98BfDHNtaXFuVMjbPYa6uqCrkREJPslC4lbgT8SC4o291By963pKiqd\nThkfprCimtd0uYSISFLJBq5/R+x0115j5siZNAx+g1eWN3DiiQVBlyMiktWSDVz/3MymJ1g2wMw+\nb2bXpae09BhYOJBhBRUsW7M26FJERLJeshv8/Rj4v2Z2PPA68D5QDBwDDAJ+AfwqrRWmwczyMNWv\n1hB7tLaIiCSS7HDTKuBKMxsInAiMAg4Ab7r72xmoLy3OOTbCM7+v5sCBz9GvX9DViIhkr5SerODu\ntUBVekvJnJPHhSme+Chr1sAppwRdjYhI9grqGdeBmjVqFnWD1/DK8sagSxERyWp9MiQGFQ1iSMEY\nnl3zVtCliIhktaQhYWYhM/tBJorJpFnlEVZsrgm6DBGRrJY0JOI38zs9A7VkVOWxEbblVfPRR0FX\nIiKSvVIauAZWmtlCYrfnaPm16u7/lZaqMuDksWH6TXqKVatg9uygqxERyU6phkQx8AFwTqs2B3I2\nJMKjwtQNXsXyFVFmzw4FXY6ISFZK9RTYz6W7kEwbXDyYwQUjeG7NOr7B1KDLERHJSimd3WRmY83s\nCTPbEZ9+Z2Zj011cus0sj7BiiwavRUQSSfUU2AeBhcDo+PT7eFtOO2tymB2havbtC7oSEZHslGpI\nDHf3B929MT49BAxPY10ZcfLYCP2OqmHlyqArERHJTqmGxAdm9un4NRMhM/s0sYHsnBYbvF7J8hVN\nQZciIpKVUg2JzwNXAtuArcDlQM4PZg/tP5SSgjKq1vwt6FJERLJS0rObzCwE/C93vyQD9WTcrPII\nK/5UTezu5yIi0lqqV1xfk4FaAnHmMWE+KKhhz56gKxERyT6pHm560cx+ZGZnmFm4eUprZRly0tgI\n/Y+qpkZnwoqIHCbVK65nxl/vaNXmtL0COyeFR4U5WFbD8uXO2Wdb0OWIiGSVVMYk8oD/5+6PZqCe\njBsxYAQDCgbyp9fW820mBV2OiEhWSWVMogn4dgZqCcwJ5WFWbK4OugwRkayT6pjEMjP7lpmNM7Mh\nzVNaK8ugM4+OsKu4hl27gq5ERCS7pBoSVwFfAf4EVMenFekqKtNOGhOh/zHVVGtnQkSkjVTvAjsx\n3YUEKTwqzMHBscHr887T4LWISLNO9yTM7Nut5q9ot+xfu7txM8szs5r4A40wszIzW2pmb5vZEjMr\n7e42UjGqZBTFBYX8+bWNmdiciEjOSHa46epW8/PaLftED2z/RmBtq/c3A8vcfQrwbAfbTJuZGrwW\nETlMspCwBPMdve+S+PMo5gD3t2q+FJgfn58PXNadbXTFGUdF2NO/hp07M7VFEZHslywkPMF8R++7\n6ofAP7brp9zdtwO4+zZgRDe3kbLI6DADJ2vwWkSktWQhcYKZ7TWzfcCM+Hzz++OPdKNm9klgu7uv\novM9ku4GUcoioyMcLKtm+fKMbVJEJOt1enaTu4fStN3ZwCVmNgfoB5SY2QJgm5mVu/t2MxsJ7EjU\nwW233dYyX1lZSWVlZbcKGlMyhvx8eGHNZiDnn8wqIkJVVRVVVVXd6sPcg/3L2czOAr7p7peY2b8B\nH7j7XWb2T0CZu9/cwWc8HXWf9fMLWTv/y7z/Qq+8K7qI9HFmhrt3aTw51YvpMuV7wMfN7G3g3Pj7\njDl9UoTaQdVs357JrYqIZK/AQ8Ldn29+oJG7f+ju57n7FHc/3913Z7KWyOgwA4+p0eC1iEhc4CGR\nTSKjYoPXK3rNDUdERLpHIdHK+NLxkF/HC6u2Bl2KiEhWUEi0YmbMKo9QvUWPqRMRAYXEYU6bEOZA\nWTVbtgRdiYhI8BQS7Zw4OkLJZA1ei4iAQuIwkdERDmjwWkQEUEgcZuLgiTTl1/LiqoQXe4uI9BkK\niXbMjJkjwqzYXEPAF6OLiAROIdGBUyeEiY6oYfPmoCsREQmWQqIDJ46OMOAYjUuIiCgkOhAeFdaV\n1yIiKCQ6dPSQo6kPfchLqz4IuhQRkUApJDqQZ3mcMGKWBq9FpM9TSCRw2oQIjKph48agKxERCY5C\nIoHIqAglkzUuISJ9m0IigdjgdY1CQkT6NIVEApOHTuZA3nb+sjKjzz0SEckqCokEQnkhjh9xggav\nRaRPU0h04tTxEfLH1fDOO0FXIiISDIVEJ8KjwpRMrtZtw0Wkz1JIdCIyOqLBaxHp0xQSnTh22LHU\n5m3m5ZV7gy5FRCQQColO5OflM3348VRvXkVTU9DViIhknkIiiVPGRyiaWM1f/xp0JSIimaeQSCI2\neK1xCRHpmxQSSURGRThYpjOcRKRvUkgkMW34NPbaBl5ZWRt0KSIiGaeQSKIgVMC04cdRs2U10WjQ\n1YiIZJZCIgUnjQ3T/+hq1q0LuhIRkcxSSKQgMipCyTEavBaRvkchkYLYldcavBaRvkchkYLjhh/H\nLvsbr9TsD7oUEZGMUkikoCi/iGOHTWXV1jU0NgZdjYhI5igkUnTSmDCDplTz1ltBVyIikjkKiRRF\nRkcYqCuvRaSPUUikKPbM62qFhIj0KYGEhJmNNbNnzewNM3vNzL4eby8zs6Vm9raZLTGz0iDq68iM\n8hl8wDqWrzwYdCkiIhkT1J5EI/AP7n4ccBrwFTM7FrgZWObuU4BngXkB1XeY4vxiJg+dzOptr9HQ\nEHQ1IiKZEUhIuPs2d18Vn68F3gTGApcC8+OrzQcuC6K+RE4cE6ZsWg1r1wZdiYhIZgQ+JmFmE4CZ\nwMtAubtvh1iQACOCq+xwkVERBk7WuISI9B35QW7czAYCjwM3unutmXm7Vdq/b3Hbbbe1zFdWVlJZ\nWZmOEtsIjwpTV/YgK1bAF76Q9s2JiHRLVVUVVVVV3erD3BP+Hk4rM8sH/gAsdvd74m1vApXuvt3M\nRgLPufvUDj7rQdS9v2E/Q743jOl/2M2KVwozvn0Rke4wM9zduvKZIA83/QJY2xwQcQuBz8bnrwee\nynRRnelf0J9JQybxxvuvU18fdDUiIukX1Cmws4HrgHPMbKWZ1ZjZJ4C7gI+b2dvAucD3gqivMyeN\niTDkuBpefz3oSkRE0i+QMQl3fxEIJVh8XiZr6arwyDCvTK5mxYq/IxwOuhoRkfQK/OymXBMZHeHg\nEN2eQ0T6BoVEF80cOZPt/jrLq3VFnYj0fgqJLhpYOJCK0vG89cGbHNQdOkSkl1NIHIETx0QYPqOa\n738fli2Dd95Bt+oQkV4p0IvpclV4VJh9n6xhw/LP8dxzsZDYuhXGjIFJk2DixNhr62nIELAunZ0s\nIhI8hcQRiIyK8Pjax3nq/kNt9fWwcWMsMJqnxx47NO9+eHA0TxUVUKhr80QkCwV2xXV3BHXFdbM9\nB/cw5t/HsPvm3eTnpZazu3a1DZDW03vvQXl54hAZPlx7ISLSfUdyxbVC4ghN/s/J/PCCH3JmxZmU\nFJV0q6/GRti0qW1wrF9/aP7gwbah0fpw1vjx0L+/QkREklNIZNA9L9/Dz2t+zvrd6ynOL2bC4AlM\nHDyx7WvZRCpKKxhQOKBb29qzp21otJ42bYJoFAYPPjSVlrZ9335qv3zAAIWMSF+gkAiAu7Nz/07W\n717Pu7vfZf2u2Ou7e2LzG/ZsoKSwhIllbQOkeb5icAXF+cXdqqGuLhYku3d3PCVbVlfXNjg6C5mO\nlg0cCHk6T04k6ykkslCTN7G9dnssQFoHSTxENu3dxNB+Q1v2PCaUTjg0P3gC40vHUxhK76h2Q8Ph\nQdKV0Nm/P7Y3UlISmwYOPDSfyvv2bTp8JpIeCokcFG2KsrV2a8seSEuQxF+37NvCiAEjOjycNWHw\nBMYOGktBqCDQn6GxEWprY9O+fW2n9m2prFNX1zZ0Ug2egQOhqCh2plhRUdupozbt/Uhfo5DohRqb\nGnlv73ttDmW1DpIt+7ZQGCpkYOHADqeSwpIjWhZk8ESjycOlo7D56KNYwNTXx15bTx21hUKphUlX\ngqd1W2HhofedzbdvKyjQnpSkh0KiD3J3DjQeoLa+9rBpX92+Dttr62upbUi8zr76feTn5XcpYPoX\n9Kd/QX/6FfSLveb3a/O+dVtxfjEW8G9B99geUPvg6ErIdNTW3N48tV6vfVui+YaGWFCkEiiJ5gsK\nYlN+/qHXRPPJlnfnc/n5CrxsopCQHuHu1EXrkgdNPFBq62vZ37Cf/Q37OdB4IPbacKDD9/sb9lMf\nrac4v7jDADksZPL7JQ2g1u8LQ4UU5BVQGCpsMxWECijIKwg8nFLhnnqgdBZE0WgscBobY1PzfEdt\niea7sm77zzU0QFNTx+GSaD7Z8q6s2z6wQqHD53vqtX1bKJSdhzOPJCR0xbUcxswozi+mOL+YYf2H\n9Xj/Td7EgYYDCQMlUcjsq9vHjo92dBpIDdEG6qP11EfraWg6NF8fraexqbFNgBSEOgiTBAHT8j4v\n+ecKQgXk5+VTkFfQZj4/L78lrJrnk67XP5+BJW2X5ULQNWtqahtWzeHRUaD0VNvBg4e3RaOHgqx5\nPtFrT61j1jY88vLaBkjr12RtR/KZjtqOhPYkpM9w98OCo3WoJAqXrqzXHEYN0QYamhpi882vrdui\nh5a1nk+2XtSjhCyUUsAU5RdRFCqiKL+IwlBhy3xRqJP2+GthqPCwtmTLmvvKtSBLB/dDAdkcKs3v\n278mmu/ptmgUbrhBh5tEejV3p7GpMWnANEQbqIvWUddYR120jvpofct869f6aP1hbS3rJ2rvoJ/W\ny5q8qSU8QnkhQhbq9DXP8pKuE7L4esn6S7Befl4+obz4a7y9eb55WbrXa/9ztq+zeXk6A1aHm0R6\nOTOL7SmECuhHv6DL6VC0KdoSHlGPEm2Kdvra5E09sk60Kb5egmWNTY1EPfZa31DfMt96Wfv1OmyL\n95n0s+3et/4Z2v88zXU3eROGJQyQjgKmK+seCe1JiIhkCXfH8U6DL5WwSbS8cmKlDjeJiEjHjuRw\nUxaepCUiItlCISEiIgkpJEREJCGFhIiIJKSQEBGRhBQSIiKSkEJCREQSUkiIiEhCCgkREUlIISEi\nIgkpJESB0/6MAAAFk0lEQVREJCGFhIiIJKSQEBGRhBQSIiKSUFaGhJl9wszeMrN1ZvZPQdcjItJX\nZV1ImFke8CPgAuA44BozOzbYqrJXVVVV0CVkDX0Xh+i7OETfRfdkXUgAJwP/4+4b3L0B+C1wacA1\nZS39D3CIvotD9F0cou+ie7IxJMYAm1q9fy/eJiIiGZaNISEiIlki655xbWanAre5+yfi728G3N3v\narVOdhUtIpIjuvqM62wMiRDwNnAusBV4FbjG3d8MtDARkT4oP+gC2nP3qJl9FVhK7HDYAwoIEZFg\nZN2ehIiIZI+cG7jWhXYxZjbWzJ41szfM7DUz+3rQNQXJzPLMrMbMFgZdS9DMrNTMHjOzN+P/Pk4J\nuqagmNlNZva6ma0xs1+ZWWHQNWWKmT1gZtvNbE2rtjIzW2pmb5vZEjMrTdZPToWELrRroxH4B3c/\nDjgN+Eof/i4AbgTWBl1ElrgHeNrdpwInAH3ycK2ZjQa+BoTdfQaxw+tXB1tVRj1I7HdlazcDy9x9\nCvAsMC9ZJzkVEuhCuxbuvs3dV8Xna4n9IuiT15OY2VhgDnB/0LUEzcwGAWe4+4MA7t7o7nsDLitI\nIWCAmeUD/YEtAdeTMe7+ArCrXfOlwPz4/HzgsmT95FpI6EK7DpjZBGAm8EqwlQTmh8A/Ahpgg4nA\nTjN7MH747T4z6xd0UUFw9y3A3cBGYDOw292XBVtV4Ea4+3aI/aEJjEj2gVwLCWnHzAYCjwM3xvco\n+hQz+ySwPb5XZfGpL8sHwsCP3T0M7Cd2iKHPMbPBxP5yrgBGAwPN7Npgq8o6Sf+wyrWQ2AyMb/V+\nbLytT4rvQj8OLHD3p4KuJyCzgUvM7B3gN8DZZvbLgGsK0nvAJndfEX//OLHQ6IvOA95x9w/dPQr8\nF/CxgGsK2nYzKwcws5HAjmQfyLWQWA4cbWYV8bMUrgb68tksvwDWuvs9QRcSFHf/jruPd/dJxP49\nPOvunwm6rqDEDyVsMrPJ8aZz6bsD+huBU82s2MyM2HfR1wbx2+9dLwQ+G5+/Hkj6x2XWXUzXGV1o\nd4iZzQauA14zs5XEdhu/4+5/DLYyyQJfB35lZgXAO8DnAq4nEO7+qpk9DqwEGuKv9wVbVeaY2a+B\nSmComW0EbgW+BzxmZp8HNgBXJu1HF9OJiEgiuXa4SUREMkghISIiCSkkREQkIYWEiIgkpJAQEZGE\nFBIiIpKQQkL6PDPbF3+tMLNrerjvee3ev9CT/Yukm0JC5ND9ayYCXbq3T/xxu535TpsNuZ/elf5F\ngqaQEDnkTuD0+N1Tb4w/yOjfzOwVM1tlZjcAmNlZZvYnM3sKeCPe9oSZLY8/AOrv4m13Av3i/S2I\nt+1r3piZfT++/mozu7JV38+1emjQggx/ByJt5NRtOUTS7Gbgm+5+CUA8FHa7+ynxe4W9aGZL4+vO\nAo5z943x959z991mVgwsN7Pfufs8M/tK/G6szTze96eAGe5+vJmNiH/m+fg6M4FpwLb4Nj/m7i+l\n8wcXSUR7EiKJnQ98Jn5vrFeAIcAx8WWvtgoIgG+Y2SrgZWJ3Jz6Gzs0mdtda3H0HUAWc1KrvrR67\nZ84qYEL3fxSRI6M9CZHEDPiauz/TptHsLOCjdu/PAU5x9zozew4obtVHqttqVtdqPor+P5UAaU9C\n5NAv6H1ASav2JcDfx5/bgZkdY2b9O/h8KbArHhDHAqe2Wlbf/Pl22/ozcFV83GM4cAbwag/8LCI9\nSn+hiBw6u2kN0BQ/vPSQu98TfzRsTfx5BDvo+JnAfwT+t5m9AbwN/KXVsvuANWZW7e5zm7fl7k+Y\n2anAaqAJ+Ed332FmUxPUJhII3SpcREQS0uEmERFJSCEhIiIJKSRERCQhhYSIiCSkkBARkYQUEiIi\nkpBCQkREElJIiIhIQv8fvyd0tVHfA+EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ab574ab10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEPCAYAAADiVdsmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FdW9///XJ4Fwk2CQm4CiKFgtakXBC9qmtIpyUGnl\nphXwlJ+1tSp6bE/L8VihrV9brY+j2G+tfqsiaMVbRSwg6JGLqFwU8YYCfSi3oCK3BAUhJJ/fHzMJ\nOyGXnWTvzJD9fj4e89gza2bW/uwt7k9mrTVrzN0RERGJg6yoAxARESmjpCQiIrGhpCQiIrGhpCQi\nIrGhpCQiIrGhpCQiIrGR8UnJzC40s4/MbI2Z/SrqeEREMpll8n1KZpYFrAG+B2wGlgOj3P2jSAMT\nEclQmX6l1B9Y6+7r3b0YmA5cGnFMIiIZK9OTUjdgY8L2prBMREQikOlJSUREYqRZ1AFErAA4OmG7\ne1hWgZllbsebiEgDuLvV5fhMv1JaDhxvZj3MLAcYBcys6kB31+LObbfdFnkMcVn0Xei70HdR81If\nGX2l5O4lZnYdMI8gQT/k7h9GHJaISMbK6KQE4O4vAidEHYeIiKj5TuooPz8/6hBiQ9/FAfouDtB3\n0TAZffNssszM9T2JyKHC3XGcUi8tX6/tNR3H9u7QG6/jQIeMb74TOdS5Bz8SJV5CSWlJteslHm5X\nsd7Q8yovZccls1RXR33rLfvRLPvhTLa8LscmltdUVt2PeNkx9Vmvqe5SLy3/d2EYZpbUa5ZlpeXY\n+tCVUhJ0pXToK/VS9pXsS+lSXFJMiZewv3T/QUtJaVjuNeyr6rxq6qvpPMcxjOysbLItmyzLKl/P\nzgq3q1hvyLFVrWdZVoWlbF8yS52OreK9siyr/AezfNusyn1Vldfl2MTyymXVbSf+WJedV5/1yvXV\n9D5xYGa6UpJ4c3f2luzlq31f8VXxV+Wvu4t3H1T21b6wPLGs+Cv2FO+pcwIp8RJaZLcgJzsnJUvz\nrOY0z25Os6xmtGzWkmZZzcqXbMuusF1hX1YN+6o5r7Zzyn6MRJoCXSklIdOvlHYX76agqICCXQXs\n/HpnrUmjyvKE5NMsqxltmrehTU4bWjdvXb5e4bV5uK9yeU4bWjVrRYtmdUsw2ZatH26RRlafKyUl\npSQ01aTk7uz4egebijZRUFQQvO6q9FpUwO7i3XTL7Ua3tt04vOXhFRJHYrKoLsFUTi7NsnSBLpIJ\nlJTS5FBMSiWlJXz25WcVkktVSadFdgu65Xaje253urWt9BqWH9HqCF1liEidKSmlSdyS0tf7vy5v\nTqsu4Xz+5ee0b9X+QHJp2/3g5JPbjcNyDov644hIE6WklCZRJ6UVn65g4oKJbCjcwKaiTezat4uu\nbbtWe2XTrW03jmx7JDnZOZHFLCKi0XdN0PMfPc/VL1zN7wf+njO6nkH33O50aN2BLNNkHCLS9Cgp\nxZS7c8+Se/jTG39i1hWz6NetX9QhiYiknZJSDO0v3c/1s6/ntY2v8ca4Nzi63dG1nyQi0gQoKcVM\n0d4iRjw9AjNj8Y8Xk9siN+qQREQajTomYmT9zvUMeHgAx+UdxwuXv6CEJCIZR0kpJpYVLOOch89h\n3Gnj+PPgP+sGUxHJSPrli4FnVz3LT2f9lIcueYhLTrgk6nBERCKjpBQhd+eu1+/ivmX3MffKufQ9\nsm/UIYmIREpJKSLFJcX8bNbPeOvTt3hj3Bt0z+0edUgiIpFTUorAzq93MuypYbRq3opX//1VTfUj\nIhLSQIdG9vGOjzn7obPp06kPM0bOUEISEUmgpNSI3tj4BgMeHsB1/a7jngvvITsrO+qQRERiRc13\njeTJ95/k+jnXM2XoFAb3Ghx1OCIisaSklGbuzv959f/w4IoHeXnMy5zS+ZSoQxIRiS0lpTTaV7KP\nn7zwE97f8j5Lxi3hyLZHRh2SiEisKSmlyfY92/nhkz8kr1UeC69aSJucNlGHJCISexrokAb/2v4v\nzvrbWfTr2o9nRzyrhCQikiQlpRR7df2rnPvwudx89s3cdcFdehifiEgdxPIX08zuNLMPzWylmT1r\nZrkJ+yaY2dpw/wUJ5X3N7F0zW2Nm9ySU55jZ9PCcN8zs6IR9Y8PjV5vZmIbG/di7j3HZU5cx9QdT\nueaMaxpanYhIxollUgLmAd90928Ba4EJAGZ2EjACOBG4CPiLmZU9//1+YJy79wZ6m9mgsHwcsN3d\newH3AHeGdeUBvwH6AWcCt5lZu/oE6+5MXDCRW+ffyitjX+GC4y6o/SQRETlILJOSu7/s7qXh5hKg\nbGK4S4Dp7r7f3dcRJKz+ZtYFaOvuy8PjpgJDw/VLgUfD9WeAgeH6IGCeuxe6+06CRHhhXWP9ev/X\nXPnclbz4rxdZMm4JfTr1qWsVIiISimVSquTHwOxwvRuwMWFfQVjWDdiUUL4pLKtwjruXAIVm1r6G\nupK2dfdWvj/1++wr2cf8sfPpfFjnupwuIiKVRDYk3MxeAhJ/xQ1w4BZ3fyE85hag2N2fSOVbp6KS\n1VtX829//zeGnzSc2793uwY0iEhslJbCvn3Bsndv1euNsV0fkSUldz+/pv1mdhUwmAPNbRBczRyV\nsN09LKuuPPGczWaWDeS6+3YzKwDyK50zv7p4Jk6cWL7e7oR2/KHgD9zxvTv48Wk/ruljiEgGKSmB\nr78Olr17o3vdvx9atICcnGBJXK/rdtl669Zw+OE1H7tq1QI++GAB2dmQnQ3vvFP379DcPfX/ZRrI\nzC4E7ga+7e7bEspPAh4nGJjQDXgJ6OXubmZLgBuA5cAsYLK7v2hm1wJ93P1aMxsFDHX3UeFAhzeB\nvgTNmG8Cp4f9S5Xj8bLvacrKKfzq5V/xxGVPMPDYgZUPFZEYKS6GPXtqXnbvTt0xJSXQsmWwtGgR\n3Wvz5mApaRNqGDPD3esUSVxndLgPyAFeCgfXLXH3a919lZk9BawCioFr/UBW/TkwBWgJzHb3F8Py\nh4BpZrYW2AaMAnD3HWb2O4Jk5MCkqhJSmVIv5dZXbmX6B9NZeNVCvtHhGyn+yCKZxz34y/7LL+Gr\nr+r2Wrnsq68OThSlpdCqVfBXfqtWNS+Vjzn88NqPqbzk5MQjGRzKYnmlFDdm5iOeHsGmok3MGDmD\njm06Rh2SSCT274fCwoOXqpJEsgmlWTM47DBo06Zur5XLWrc+OGnE5YohUzWlK6XYybIs/nfM/9Ky\nWcuoQxGpl337qk4ohYVQVFT9vsRl717IzYV27Q4subnQtm3FJJGXB0cdVXtCadMmSBwiZXSllAQz\n89LSUkx/ckmE3IPEsHUrfPFFsGzfDjt3JpdQiosrJpP6LG3a6MpDklefKyUlpSQkDnQQSZX9+4ME\nk5hkalrfujXoyO7YMVg6dIAjjgj6PipfuVSVUFq3VkKRxqWklCZKSpKMr76qPbEkru/aFTRzJSaZ\nmtY7dAhGVokcKpSU0kRJKbPt3QsbN8L69bBuXfC6fj18/nnFJONedVKpLsnk5UGW7rmWJkxJKU2U\nlJq2PXsOJJqypJP4unUrdO0KxxwDPXocWI48smKyUfOYSEVKSmmipHRo27Wr5qRTWBiMFCtLOpVf\nu3YNhi2LSN0oKaWJklK87dxZdbIpe92zp+pkU/bapYua0UTSQUkpTZSUoldaCh9+CEuXwsqVFZOO\ne81Jp2NHNauJREFJKU2UlBrfli1BAlq6FJYsgeXLoVMnOPNM6NsXevY8kHgOP1xJRySOlJTSREkp\nvfbuDa5+liw5kIR27ID+/YMkdNZZwXqHDlFHKiJ1oaSUJkpKqeMOn3xyIPksXQrvvQe9ewfJpywJ\n9e6tfh6RQ52SUpooKdVfYWHQ9JaYhHJyDiSfM8+E008Ppq8RkaZFSSlNlJSSU1ICH3xQsRlu/fqg\nDygxCXXvHnWkItIYlJTSREmpap9+WvEK6K23gnt6Epvh+vTRLNAimUpJKU2UlALvvQcvvXQgCX35\nZcUroP79g6lzRERASSltMj0prVkD//3f8Npr8IMfHEhCxx+vodgiUj095E9SavNmmDQJ/vEPuPlm\nmDIlmN9NRCRdNOhWDrJjB/z613DyycGNqatXB9tKSCKSbkpKUm73bvjjH4N7hLZvh3feCbbbt486\nMhHJFGq+E4qL4ZFH4Le/hbPPhsWL4YQToo5KRDKRklIGKy2FZ58NBjF07w7PPQf9+kUdlYhkMiWl\nDPXyy0E/kTv8+c9w/vlRRyQioqSUcd58M0hGGzbA7bfDZZdpjjkRiQ/9HGWI1ath+HC49FIYMSKY\nDmj4cCUkEYkX/SQ1cQUF8JOfwLnnwhlnwNq1wbam/hGROFJSaqK2b4df/QpOOSUY0r1mTbCte41E\nJM6UlJqY3bvhD38IhnTv3Anvvhtsa046ETkUKCk1EcXF8Ne/Qq9esGJFME/dAw9At25RRyYikrxY\nJyUzu9nMSs2sfULZBDNba2YfmtkFCeV9zexdM1tjZvcklOeY2fTwnDfM7OiEfWPD41eb2ZjG+2Sp\nU1oKTz4JJ50U3HP0/PPw1FPBrAwiIoea2A4JN7PuwPnA+oSyE4ERwIlAd+BlM+sVTuF9PzDO3Zeb\n2WwzG+Tuc4FxwHZ372VmI4E7gVFmlgf8BugLGPCWmT3v7oWN+Tnryz14jMSECcFM3fffD9//ftRR\niYg0TJyvlP4H+GWlskuB6e6+393XAWuB/mbWBWjr7svD46YCQxPOeTRcfwYYGK4PAua5e6G77wTm\nARem5ZOk2LJl8L3vwfXXB0lp+XIlJBFpGmKZlMzsEmCju79XaVc3YGPCdkFY1g3YlFC+KSyrcI67\nlwCFYXNgdXXF1kcfwbBh8MMfwuWXB/caDRumZxqJSNMRWfOdmb0EdE4sAhz4b+C/CJru0vLW9Tlp\n4sSJ5ev5+fnk5+enKJzaff11cFU0Ywb88pcwdaqGdotI/CxYsIAFCxY0qI7IkpK7V5l0zKwPcAzw\njpkZQd/RCjPrT3A1c3TC4d3DsgLgqCrKSdi32cyygVx3325mBUB+pXPmVxdvYlJqbHPmBEO716zR\n0G4Ria/Kf7BPmjSpznXErvnO3d939y7u3tPdjyVoijvN3bcAM4GR4Yi6Y4HjgWXu/hlBs1z/MJGN\nAZ4Pq5wJjA3XhwOvhOtzgfPNrF046OH8sCx2Xn0Vhg5VQhKRpi+2o+8SOGGTm7uvMrOngFVAMXBt\nOPIO4OfAFKAlMNvdXwzLHwKmmdlaYBswKqxrh5n9DngzfI9J4YCH2Fm0CO69N+ooRETSzw78plex\n0+yHNZ3s7v9IeUQxZGZe0/eUTkVF0LUrbNsGLVpEEoKISL2YGe5ep3782q6ULg5fOwHncKDp67vA\n60BGJKUovf568OA9JSQRyQQ1JiV3/3cAM5sHnOTun4bbRxI0lUmaLVoE550XdRQiIo0j2YEOR5Ul\npNDnVBwFJ2ny6qvw7W9HHYWISONIdqDD/5rZXOCJcHsk8HJ6QpIye/bA22/D2WdHHYmISONIKim5\n+3Vm9gOg7G/2B939ufSFJRBMJ9SnD7RpE3UkIiKNoy5Dwl8H9hMMn16WnnAk0aJFaroTkcySVJ+S\nmY0gSETDCGbpXmpmw9IZmGiQg4hknhrvUyo/yOwd4PxwVgXMrCPwsrufmub4YiGK+5SKi4PHmG/Y\noJkcROTQVJ/7lJIdfZdVlpBC2+pwrtTD229Dz55KSCKSWZLtU3qxitF3s9MTkoD6k0QkMyU7+u6X\n4ZRD54ZFGn2XZosWwejRUUchItK4kupTAjCzzkB/wtF3lZrzmrTG7lMqLYUOHWDVKujSpdHeVkQk\npdLWp6TRd43r/feDpKSEJCKZJtk+pVuAfpVH3wHPpCuwTKaphUQkU2n0XQxpkIOIZKpkE8uLZjbX\nzK4ys6uAWWj0XVq4KymJSOaqy0CHy4AB4earmTT6rjEHOqxdCwMHBjfNWp26B0VE4iUdD/kr5+7P\nAs/WOSqpk7KrJCUkEclEyY6++6GZrTWzQjMrMrNdZlaU7uAykQY5iEgmS3buu38BF7v7h+kPKX4a\ns/muZ0+YNQtOPLFR3k5EJG3SOffd55makBrTxo2waxd84xtRRyIiEo0a+5TCqYUA3jSzJ4EZwN6y\n/e7+jzTGlnHKmu7UnyQimaq2gQ4XJ6zvBi5I2HZASSmF9PwkEcl0SQ8Jz2SN1af0zW/CtGnQt2/a\n30pEJO3q06dUY1Iys/909zvN7D6CK6MK3P2Guod56GmMpPTFF9CrF2zbBtnZaX0rEZFGkY77lMoG\nN7xZv5AkWYsXwznnKCGJSGarMSm5+wvh66ONE07m0tRCIiK1j757gSqa7cq4+yUpjyhDLVoEkydH\nHYWISLRq61P6Tk0nu/vClEd04L2vB64F9gOz3P3XYfkE4Mdh+Xh3nxeW9wWmAC2B2e5+Y1ieA0wF\nTge2AiPdfUO4byzBYzkcuN3dp1YTS1r7lIqKoGvXoD+pRYu0vY2ISKNKeZ9SYtIxs1bA0e6+up7x\nJc3M8gmGo5/s7vvNrENYfiLBQwZPBLoDL5tZrzBj3A+Mc/flZjbbzAa5+1xgHLDd3XuZ2UjgTmCU\nmeUBvwH6Aga8ZWbPu3thuj9fZa+/Dv36KSGJiCQ7993FwErgxXD7W2Y2M41x/Qz4g7vvB3D3rWH5\npcB0d9/v7uuAtUB/M+sCtHX35eFxU4GhCeeU9Yk9AwwM1wcB89y90N13AvOAC9P4maql/iQRkUCy\n0wxNBPoDOwHcfSVwbJpiAugNfNvMlpjZfDM7PSzvBmxMOK4gLOsGbEoo3xSWVTjH3UuAQjNrX0Nd\njU5JSUQkkOyjK4rdvdAqzn/ToE4WM3sJ6JxYFNb532Fcee5+lpn1A54Gejbk/Sq9T51NnDixfD0/\nP5/8/PyUBLNnD7z9Npx1VkqqExGJzIIFC1iwYEGD6kg2KX1gZlcA2WbWC7gBeL0hb+zu51e3z8x+\nSjiFUdhHVGJmRxBczRydcGj3sKwAOKqKchL2bTazbCDX3bebWQGQX+mc+dXFlJiUUmnpUjj5ZGjT\nJi3Vi4g0msp/sE+aNKnOdSTbfHc98E2CyVj/DhQBN9b53ZI3g7Dvx8x6Aznuvg2YCYw0sxwzOxY4\nHljm7p8RNMv1t+BybgzwfFjXTGBsuD4ceCVcnwucb2btwkEP54dljUrPTxIROSDZK6XO7n4LwfBp\nAMJmteXVn9IgjwAPm9l7BIlwDIC7rzKzp4BVQDFwbcJY7Z9TcUj4i2H5Q8A0M1sLbANGhXXtMLPf\nEcxW4cCkcMBDo1q0CMaPb+x3FRGJp2Qf8reC4CF/BeH2t4H/6+4npzm+WEjXfUrFxdC+PWzYAHl5\nKa9eRCRS6XzI3zXADDPrYmaDgfuAwXUNUCpasSJ40qwSkohIIKnmu3CwwQ0E9/J8DXzf3b9Ia2QZ\nQEPBRUQqquvcd62BQuCh8LJMc981wKuvwujRUUchIhIfsZ37Lk7S0adUWgodOsCqVdClS0qrFhGJ\nhbTOfSep9f770LGjEpKISKLamu8Wu/u5ZraLis14Bri756Y1uiZs0SI477yooxARiZfarpTODV/b\nNk44mWPRIhgyJOooRETiJdkh4Qcxsw2pDCSTuGsmBxGRqtQ7KVHPiU0F/vUvaN4cevSIOhIRkXhp\nSFJK36NYm7iy+5NMaV1EpILaBjr8R3W7gMNSH05m0CAHEZGq1Xal1Laa5TDg3vSG1nRpJgcRkaol\nNSFrpkvlzbMbN0LfvrBli5rvRKRpS+eErIlvsqKu58gBZaPulJBERA5Wn4EO+jltADXdiYhUrz5J\naVbKo8ggGuQgIlI99SklIVV9Sl98AccfD9u3Q3Z2CgITEYmxlE/ImlBx5bnvIHiExZvAze7+cV3e\nNFMtXgwDBighiYhUJ6mkBNwDbAL+TtCnNAo4DlgBPAzkpyO4pkb9SSIiNUu2T+kSd3/A3Xe5e5G7\nPwgMcvcnAT3MO0lKSiIiNUs2Ke02sxFmlhUuIwgeiw6abigphYWwejWccUbUkYiIxFeySelHwGhg\nS7iMBq40s1bAdWmKrUl5/XXo1w9ycqKOREQkvpLqUwoHMlxcze7FqQun6VLTnYhI7ZK6UjKz7mb2\nnJltCZdnzax7uoNrSvT8JBGR2iXbfPcIMBPoGi4vhGWShD17YOVKOOusqCMREYm3ZJNSR3d/xN33\nh8sUoGMa42pSli6FPn2gTZuoIxERibdkk9I2M7vSzLLD5UpgWzoDa0rUnyQikpxkk9KPgRHAZ8Cn\nwDDgqjTF1OQoKYmIJCeppOTu6939Enfv6O6d3H0ocFm6gjKzU83sDTN728yWmdkZCfsmmNlaM/vQ\nzC5IKO9rZu+a2RozuyehPMfMpofnvGFmRyfsGxsev9rMxqTjsxQXB813Awako3YRkaalPrOEl6nu\nUempcCdwm7ufBtwG3AVgZicRXLGdCFwE/MWs/MlE9wPj3L030NvMBoXl44Dt7t6LYLqkO8O68oDf\nAP2AM4HbzKxdqj/IihVw3HGQp3kvRERq1ZCklM7nKpUCZQnicKAgXL8EmB4OtlgHrAX6m1kXoK27\nLw+PmwoMDdcvBR4N158BBobrg4B57l7o7juBecCFqf4gelSFiEjykp2QtSrpnF7oJmCumd1NkPzO\nCcu7AW8kHFcQlu0nmDC2zKawvOycjQDuXmJmhWbWPrG8Ul0ptWgRjElLw6CISNNTY1Kq5pEVECSK\nVg15YzN7CehcqU4HbgG+D4x39xlmNoxgJvLzG/J+ld6nziZOnFi+np+fT35+fq3nlJYGj6v4f/+v\nPu8oInJoWbBgAQsWLGhQHbF8yJ+Z7XT3wytvm9mvAXf3P4blLxL0Oa0H5rv7iWH5KOA77v6zsmPc\nfamZZQOfunun8Jh8d/9peM5fwzqerCKeej3k7913YfjwYCJWEZFMU5+H/DWkTymdCszsOwBm9j2C\nviMIZpUYFY6oOxY4Hljm7p8BhWbWPxz4MAZ4PuGcseH6cOCVcH0ucL6ZtQsHPZwflqWMhoKLiNRN\nQ/qU0ulqYHJ4ZfM18BMAd19lZk8Bq4Bi4NqES5ifA1OAlsBsd38xLH8ImGZmawlu+B0V1rXDzH5H\n8PRcByaFAx5SZtEiGDIklTWKiDRtsWy+i5v6NN+5w5FHwpIlcMwx6YlLRCTOmlLz3SFv7drg2Uk9\nekQdiYjIoUNJKU3KHlVh6bybS0SkiVFSShMNchARqTslpTTRTA4iInWnpJQGGzbAl1/CN74RdSQi\nIocWJaU0UH+SiEj9KCmlQVlSEhGRulFSSgMNchARqR/dPJuEutw8u2UL9O4N27ZBdnaaAxMRiTHd\nPBsDixfDOecoIYmI1IeSUoqp6U5EpP6UlFJMSUlEpP7Up5SEZPuUCguhe/egPyknpxECExGJMfUp\nRez116FfPyUkEZH6UlJKIU0tJCLSMEpKKaT+JBGRhlGfUhKS6VPaswc6dAjuU2rTppECExGJMfUp\nRWjpUjjlFCUkEZGGUFJKETXdiYg0nJJSimiQg4hIw6lPKQm19Snt2wdHHBE8RykvrxEDExGJMfUp\nRWTFCjjuOCUkEZGGahZ1AE2Bnp8kkl7HHHMM69evjzoMqUaPHj1Yt25dSupSUkqBRYtg7NiooxBp\nutavX4+6GuLLUviYbfUpJaGmPqWSkuD+pI8+gs6dGzkwkQwR9k1EHYZUo7r/PupTisD770OnTkpI\nIiKpoKTUQLo/SUQkdZSUGkiDHEREUieypGRmw8zsfTMrMbO+lfZNMLO1ZvahmV2QUN7XzN41szVm\ndk9CeY6ZTQ/PecPMjk7YNzY8frWZjUkoP8bMloT7njCzOg/6cNeVkoikTmlpKW3btmXTpk0pPfZQ\nEuWV0nvAD4CFiYVmdiIwAjgRuAj4ix0Y2nE/MM7dewO9zWxQWD4O2O7uvYB7gDvDuvKA3wD9gDOB\n28ysXXjOH4G7w7p2hnXUydq1wbOTevSo65ki0hS0bduW3NxccnNzyc7OpnXr1uVlTzzxRJ3ry8rK\nYteuXXTv3j2lxx5KIktK7r7a3dcClUdmXApMd/f97r4OWAv0N7MuQFt3Xx4eNxUYmnDOo+H6M8DA\ncH0QMM/dC919JzAPuDDcNxB4Nlx/lCBB1omukkQy265duygqKqKoqIgePXowa9as8rLLL7/8oONL\nSkoiiPLQEsc+pW7AxoTtgrCsG5B4nbopLKtwjruXAIVm1r66uszsCGCHu5cm1NW1roEqKYlIGXc/\naFj0rbfeyqhRo7jiiito164djz/+OEuWLOHss88mLy+Pbt26MX78+PJkVVJSQlZWFhs2bABg9OjR\njB8/nsGDB5Obm8uAAQPKbyKuy7EAc+bM4YQTTiAvL48bbriBc889l6lTpzbGV1MnaU1KZvZS2AdU\ntrwXvl6czvfl4Kuv+h5TIw1yEJHazJgxgyuvvJLCwkJGjhxJ8+bNmTx5Mtu3b+e1115j7ty5PPDA\nA+XHV74R9YknnuD2229nx44dHHXUUdx66611PnbLli2MHDmSu+++m61bt3LssceyfPly4iitMzq4\n+/n1OK0AOCphu3tYVl154jmbzSwbyHX37WZWAORXOme+u28zs3ZmlhVeLSXWVaWJEyeWr+fn59Oz\nZz5ffQUnnFD3DygiqZWqCQXScX/uueeey+DBgwFo0aIFp59+evm+Y445hquvvpqFCxdy7bXXhjFU\nDGLYsGGcdtppAPzoRz/illtuSYg3uWNnzZrFaaedxpAhQwC46aabuOuuu1L5MQFYsGABCxYsaFAd\ncZlmKPGf1EzgcTP7H4Lmt+OBZe7uZlZoZv2B5cAYYHLCOWOBpcBw4JWwfC5wezi4IQs4H/h1uG9+\neOyT4bnP1xRgYlICePzx4FEVKZxdQ0TqKc6TPRx11FEVtlevXs3NN9/MW2+9xe7duykpKeHMM8+s\n9vwuXbqUr7du3Zovv/yyzsdu3rz5oDjSMUAiPz+f/Pz88u1JkybVuY4oh4QPNbONwFnAP81sDoC7\nrwKeAlYBs4FrE+b4+TnwELAGWOvuL4blDwEdzGwtcCNh4nH3HcDvgDcJEtakcMAD4TH/YWZrgPZh\nHUlTf5LRs/WUAAAN40lEQVSIJKNyE9s111zDySefzMcff0xhYSGTJk1K+xRKRx55JBs3bqxQVlBQ\nY+NQZCK7UnL3GcCMavbdAdxRRflbwMlVlO8lGEZeVV1TgClVlH9CMEy8XhYtgp/+tL5ni0im2rVr\nF+3ataNVq1Z8+OGHPPDAA2kf1j1kyBBuvPFGZs2axYUXXsh9993H1q1b0/qe9RXH0Xext2ULfPop\nnHJK1JGISFwkO1P23XffzZQpU8jNzeVnP/sZo0aNqrae2upM9thOnTrx5JNPctNNN9GhQwc++eQT\nTjvtNFq0aJFUzI1Js4QnofIs4f/4Bzz0EMyaFWFQIhlEs4SnVmlpKV27duXZZ59lwIABDa5Ps4RH\nbNGiYJCDiMihYu7cuRQWFrJ3715++9vfkpOTQ//+/aMO6yBKSvWgQQ4icqhZvHgxPXv2pHPnzrz0\n0kvMmDGD5s2bRx3WQdR8l4TE5rvCQujWDbZvD+a9E5H0U/NdvKn5LkKvvQb9+yshiYikg5JSHWlq\nIRGR9FFSqiMNchARSR/1KSWhrE9p927o2DG4T6lNm6ijEskc6lOKN/UpRWTp0uCGWSUkEZH0UFKq\nAw0FF5FUWb9+PVlZWZSWBo91Gzx4MNOmTUvq2Lq64447+MlPflLvWBuTklIdaJCDiCS66KKLDnqC\nAMDzzz/PkUceWWsSSZwaaPbs2YwePTqpY2uycOHCg2YEnzBhAg8++GBS50dNSSlJ+/YFzXcpmJFD\nRJqIsWPH8thjjx1U/thjjzF69Giyshr/J9bdk05gcaSklKQVK+C44+Dww6OORETiYujQoWzbto3F\nixeXl+3cuZN//vOfjB49mtmzZ9O3b1/atWtHjx49any+0He/+10efvhhIJib7he/+AUdO3bk+OOP\nZ1aliTanTJnCSSedRG5uLscff3z5VdDu3bsZPHgwmzdvpm3btuTm5vLZZ58xadKkCldhM2fOpE+f\nPrRv356BAwfy0Ucfle879thjufvuuzn11FPJy8vj8ssvZ9++fSn5vpKhpJQk9SeJSGUtW7Zk+PDh\nTJ06tbzsySef5MQTT+Tkk0+mTZs2TJs2jcLCQmbNmsVf//pXZs6cWWu9Dz74ILNnz+add97hzTff\n5Jlnnqmwv3PnzsyePZuioiIeeeQRbrrpJlauXEnr1q2ZM2cOXbt2ZdeuXRQVFZU/+K/s6mnNmjVc\nccUVTJ48mS+++IKLLrqIiy++mP3795fX//TTTzNv3jw++eQT3nnnHaZMmZKCbys5cXnybOwtWgRX\nXRV1FCJSFZuUmuYqv63uw87Hjh3LkCFD+POf/0xOTg7Tpk1j7NixAHznO98pP65Pnz6MGjWKhQsX\ncskll9RY59NPP82NN95I165dgaBPaOHCheX7L7roovL18847jwsuuIBXX32Vb33rW7XG+9RTTzFk\nyBAGDhwIwC9+8QvuvfdeXn/9db4d/uU9fvx4OnfuDMDFF1/MypUrk/kqUkJJKUmvvRY8rkJE4qc+\nySRVBgwYQMeOHZkxYwZnnHEGy5cv57nnngNg6dKlTJgwgffff599+/axb98+hg8fXmudlR9f3qNH\njwr758yZw29/+1vWrFlDaWkpe/bs4ZQkH/C2efPmCvWZGUcddVSFJ9GWJSQIHqv+6aefJlV3Kqj5\nLkmdOkHCfycRkXKjR4/m0Ucf5bHHHmPQoEF07NgRgB/96EcMHTqUgoICdu7cyTXXXJPUTcCVH1++\nfv368vV9+/YxbNgw/vM//5MvvviCHTt2cNFFF5XXW9sgh65du1aoD2Djxo1pf/ptspSUkqT+JBGp\nzpgxY3j55Zf529/+Vt50B/Dll1+Sl5dH8+bNWbZsGX//+98rnFddghoxYgSTJ0+moKCAHTt28Mc/\n/rF8X9kVV4cOHcjKymLOnDnMmzevfH/nzp3Ztm0bRUVF1dY9a9Ys5s+fz/79+/nTn/5Ey5YtOfvs\nsxvyFaSMklKSlJREpDo9evTgnHPOYffu3RX6i/7yl79w66230q5dO37/+98zcuTICudV9zjzq6++\nmkGDBnHqqadyxhlncNlll5XvO+yww5g8eTLDhw+nffv2TJ8+nUsvvbR8/wknnMDll19Oz549ad++\nPZ999lmF9+zduzePPfYY1113HR07dmTWrFm88MILNGvW7KA4oqC575JgZr5unVOpWVdEGonmvou3\nVM59p6SUhMSH/IlI41NSijdNyCoiIk2SkpKIiMSGkpKIiMSGkpKIiMSGkpKIiMSGphkSkdjr0aNH\n5PfPSPUqT4PUEJFdKZnZMDN738xKzKxvQvn3zexNM3vHzJab2XcT9vU1s3fNbI2Z3ZNQnmNm081s\nrZm9YWZHJ+wbGx6/2szGJJQfY2ZLwn1PmJkStEhMrVu3DnfXEtNl3bp1KftvHWXz3XvAD4CFlcq/\nAIa4+6nAVUDi84HvB8a5e2+gt5kNCsvHAdvdvRdwD3AngJnlAb8B+gFnAreZWbvwnD8Cd4d17Qzr\nkFosWLAg6hBiQ9/FAfouDtB30TCRJSV3X+3uawGrVP6Ou38Wrn8AtDSz5mbWBWjr7svDQ6cCQ8P1\nS4FHw/VngIHh+iBgnrsXuvtOYB5wYbhvIPBsuP4oQYKUWuh/uAP0XRyg7+IAfRcNE+uBDmY2DFjh\n7sVAN2BTwu5NYRnh60YAdy8BCs2sfWJ5qADoZmZHADvcvTShrq5p+yAiIpKUtPajmNlLQOIDHwxw\n4BZ3f6GWc78J3AGcX5+3TtExIiLSmKLuIAPmA30rlXUHVgNnJZR1AT5M2B4F3B+uvwicGa5nA1sS\njvlrwjl/BUaG61uArHD9LGBODTG6Fi1atGip+1LXnBCXEWflVy3hQIR/Ar9y9yVl5e7+mZkVmll/\nYDkwBpgc7p4JjAWWAsOBV8LyucDtYZ1ZBFddvw73zQ+PfTI89/nqgvM6TigoIiL1E9ks4WY2FLgP\n6EAw+m2lu19kZrcQJI6yQRAOXODuW83sdGAK0BKY7e7jw7paEIzSOw3YBoxy93XhvquAW8J6fu/u\nU8PyY4HpQB7wNnBl2HclIiIR0aMrREQkNmI9+i5qZnahmX0U3mD7q6jjiYqZdTezV8zsAzN7z8xu\niDqmqJlZlpmtMLOZUccSJTNrZ2ZPm9mH4b+PM6OOKSpmdlM4IcC7Zva4meVEHVNjMbOHzOxzM3s3\noSzPzOaFExfMTbhHtEZKStUwsyzgzwT3On0TuNzMvhFtVJHZD/yHu38TOBv4eQZ/F2XGA6uiDiIG\n7iVoSj8ROBX4MOJ4ImFmXYHrCQZtnUIwsnlUtFE1qkcIfisT/Rp42d1PIOjnn5BMRUpK1esPrHX3\n9WFf03SCm3Qzjrt/5u4rw/UvCX54utV8VtNlZt2BwcDfoo4lSmaWC5zn7o8AuPt+dy+KOKwoZQNt\nwinLWgObI46n0bj7YmBHpeLESQ0e5cBkBzVSUqpe5RtvE2/WzVhmdgzwLYKRjpnqf4BfEgyeyWTH\nAlvN7JGwKfNBM2sVdVBRcPfNwN3ABoKb9He6+8vRRhW5Tu7+OQR/2AKdkjlJSUmSZmaHEUzjND68\nYso4ZvZvwOfhlaOR2TdhNwP6Av/X3fsCuzlwy0VGMbPDCa4MehDMDnOYmV0RbVSxk9QfcUpK1SsA\njk7Y7h6WZaSwSeIZYJq7V3tPVwYYAFxiZh8DTwDfNbOpEccUlU3ARnd/M9x+hiBJZaLvAx+7+/Zw\nqrN/AOdEHFPUPjezzgDh3KVbkjlJSal6y4HjzaxHOIpmFMFNupnqYWCVu98bdSBRcvf/cvej3b0n\nwb+JV9x9TG3nNUVh08xGM+sdFn2PzB38sQE4y8xaWvDgp++ReYM+KrcczCR40gPUMkFBorjM6BA7\n7l5iZtcRzCyeBTzk7pn2jwwAMxsA/Ah4z8zeJrgM/y93fzHayCQGbgAeN7PmwMfAv0ccTyTcfZmZ\nPUNwI35x+PpgtFE1HjP7O5APHGFmG4DbgD8AT5vZj4H1wIik6tLNsyIiEhdqvhMRkdhQUhIRkdhQ\nUhIRkdhQUhIRkdhQUhIRkdhQUhIRkdhQUhJpRGa2K3ztYWaXp7juCZW2F6eyfpHGoKQk0rjKbgw8\nFqjT3Ghmll3LIf9V4Y3cz61L/SJxoKQkEo07gHPD2bXHhw8NvNPMlprZSjO7GsDMvmNmi8zseeCD\nsOw5M1sePnDx/wvL7gBahfVNC8t2lb2Zmd0VHv+OmY1IqHt+wkP6pjXydyByEE0zJBKNXwM3u/sl\nAGES2unuZ4ZzLb5mZvPCY08DvunuG8Ltf3f3nWbWElhuZs+6+wQz+3k4W3cZD+u+DDjF3U82s07h\nOQvDY74FnAR8Fr7nOe7+ejo/uEhNdKUkEg8XAGPCuQWXAu2BXuG+ZQkJCeBGM1sJLCGYvb4XNRtA\nMKM57r4FWAD0S6j7Uw/mG1sJHNPwjyJSf7pSEokHA65395cqFJp9B/iq0vZA4Ex332tm84GWCXUk\n+15l9iasl6DfBImYrpREGldZQtgFtE0onwtcGz63CjPrZWatqzi/HbAjTEjfAM5K2Lev7PxK7/Uq\nMDLst+oInAcsS8FnEUk5/VUk0rjKRt+9C5SGzXVT3P3e8FHzK8Ln8WwBhlZx/ovAT83sA2A18EbC\nvgeBd83sLXcfXfZe7v6cmZ0FvAOUAr909y1mdmI1sYlERo+uEBGR2FDznYiIxIaSkoiIxIaSkoiI\nxIaSkoiIxIaSkoiIxIaSkoiIxIaSkoiIxIaSkoiIxMb/DwKOifWPGfMwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ab6369e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def logreg(x,t,W,b):\n",
    "    lq = dot(W.T,x) + b\n",
    "    Z = sum(np.exp(lq))\n",
    "    a = max(lq)\n",
    "    lZ = a  + log(sum(np.exp(lq - ones(shape(lq))*a)))\n",
    "    lp = lq - ones(shape(lq))*log(Z)\n",
    "    lpt = lp[t]\n",
    "    return lpt, argmax(lp)\n",
    "\n",
    "def pred_iter(x_set, t_set, W, b):\n",
    "    N_set, M = shape(x_set)\n",
    "    L = 0\n",
    "    misclass = 0\n",
    "    for ix in range(N_set):\n",
    "        x = reshape(x_set[ix],(M,-1))\n",
    "        t = t_set[ix]\n",
    "        lpt, lpix = logreg(x,t,W,b)\n",
    "        L += lpt\n",
    "        misclass += lpix != t\n",
    "    return L, misclass*1.0/N_set\n",
    "\n",
    "\n",
    "N_tr, M = shape(x_train)\n",
    "N_vl, M = shape(x_test)\n",
    "\n",
    "K = 10\n",
    "max_iter = 10\n",
    "\n",
    "W = np.zeros((M,K))\n",
    "b = np.zeros((K, 1))\n",
    "\n",
    "L_tr = []\n",
    "E_tr = []\n",
    "L_vl = []\n",
    "E_vl = []\n",
    "\n",
    "L, err = pred_iter(x_train, t_train, W, b)\n",
    "L_tr.append(L)\n",
    "E_tr.append(100*err)\n",
    "L, err = pred_iter(x_valid, t_valid, W, b)\n",
    "L_vl.append(L)\n",
    "E_vl.append(100*err)\n",
    "    \n",
    "for it in range(max_iter):\n",
    "    #To visualize weights after each iteration uncomment the following line\n",
    "    #plot_digits(W.T, numcols=10)\n",
    "    W, b = sgd_iter(x_train, t_train, W, b)\n",
    "    L, err = pred_iter(x_train, t_train, W, b)\n",
    "    L_tr.append(L)\n",
    "    E_tr.append(100*err)\n",
    "    L, err = pred_iter(x_valid, t_valid, W, b)\n",
    "    L_vl.append(L)\n",
    "    E_vl.append(100*err)\n",
    "    \n",
    "plot(range(max_iter + 1), E_tr, label='Training')\n",
    "plot(range(max_iter + 1), E_vl, label='Validation')\n",
    "xlabel('Iteration')\n",
    "ylabel('Error (%)')\n",
    "legend(loc=1,ncol=1)\n",
    "show()\n",
    "\n",
    "plot(range(max_iter + 1), L_tr, label='Training')\n",
    "plot(range(max_iter + 1), L_vl, label='Validation')\n",
    "xlabel('Iteration')\n",
    "ylabel('Log-Likelihood')\n",
    "legend(loc=4,ncol=1)\n",
    "show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.2 Visualize weights (10 points)\n",
    "Visualize the resulting parameters $\\bW$ after a few iterations through the training set, by treating each column of $\\bW$ as an image. If you want, you can use or edit the `plot_digits(...)` above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAADTCAYAAABOQ5KuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnUmPLEfVht+L53n29ezrGczMwhuwkADBr2DPjiVihVkh\nluzY8yPYWUgIJGQZsI2N8XA9z/M829/i09Nx6q3I6u7blVUdvu+zya7qrMyMyMiM95w4ceLIF198\noRBCCGPxlW1fQAghhP2Tl3cIIQxIXt4hhDAgeXmHEMKA5OUdQggDkpd3CCEMSF7eIYQwIHl5hxDC\ngOTlHUIIA3Lq3Cf45S9/edJM4fzDH/5wZC/7/fWvfz1p6uTOO+/cU53cddddG6+TOrv4yJE9XeZa\nuOuuu/Z8st/85jcnRVv57W9/u+c6+f3vf39S1Ikk/epXv5qsl9lf3puAh3DVA7iXfcLJRW0LtA+2\nX/lKjNIvI3t9/tlvP+lDNp1qJC00hBAGZCjlTW/ovad/X/8/1YP6588//3zh+1ESdu2mJHrlcCvE\nt1NlH6VOvHx+3avaCdtTTjllYQufffaZpOX2wudV5zusTF3fftXnKO1Dmr7ffD711FMXPtey+Tti\nqk3s5V1ykDqL8g4hhAE5lMp7SmHTO+KP5P987qmlKXXJlt6S3vOTTz6RJH366acLW+lwqHK/fso+\n5aNdpby9Hl1luKp0ZbHbeeZiSvl8/PHHkqT3339fkvTBBx9Ikj766CNJ7V7yuzPOOGPnmPzN9txz\nz5UkXXjhhZKk008/XVJTZK62+Fz/3mSdcB2r1LLXF22drdcrZfXn6rTTTlv4P7iVUq9r26p8SmlT\nBso0ta3X7++Gqfrk+9o2/FgHIco7hBAG5FAp790UNb0k39Or0wN++OGHkhZ7Or6b6lH9mH4tVV1M\n+bbmZDfrw+uMa0R9vvvuu0vHREWiMvkt9chvULLsf+aZZ0qSzjrrrJ1j+fk3USfUP4oahf3WW29J\nkp577jlJ0iuvvCJJevXVVyVJb7/99sJx6r09++yzJUkXXXSRJOmWW26RJF1//fWSpEsvvVSSdM45\n50hqdUd5uZb6d7XaDspu4xDcK+4/dSK1+8mz8M477yxcJ88CZbvgggskNeuD77n/3Gs+ezuq10Md\ncD2uQreF+7Zp026B0fYrtD/KxJbvXYFTB24d1fOfSHRTlHcIIQxIXt4hhDAgW3ebVBPC3SSYLO7i\nwAx57733FrYvv/yyJOm1117bOSbuA8A8uvjiiyVJ559/vqRmKvq5+FzP6ybqJvBBI68b8EE7TONa\nz5i7bP23Xq+YkNyfOtC3KbdJPS7mKOYq9/2pp56SJD3//POSpNdff11Sq4upAVmpuRaoL8qFywAX\nAu4VH/StdTk10HsQpq6d6+X633jjDUmLzwDf4TZyVyKDslxnvb/1XNSBtz32r3VAvdCGgHuxjjrZ\nC34PuC6/R9QvZWI/f76kZXed1y/t098TfF9dJLx/eC/tx30S5R1CCANyKJU3PfvUwBo9nyswlHE9\nJt+hwugFUSZXXHHFwjFRVmzPO++8nWP5IFVvgsY6WGWNcA301NQR+1EnriRqGBdWBkrJQ55QmxwT\nxUVd1EEcHzyeK5yyF6rlahNFSbmuuuoqSa0O/PrrffOyUkeoUm+D/LZnobnCWkf78Pp0K4k6YLD2\nzTff3NkXRchv3IKiPl988UVJbeDXlTmqmWcD1cjx2E+SLrvsMkntXnBveM7mVN712B5kQFl5Htz6\n5LmibFjotWy0AR8gf+aZZyQ1S9CDBxgors9iLwhgr0R5hxDCgGxNeXs4oDQ9GYB9UNgeBI8qcSUp\nLfs7+cwxUGGoTc5NT10VFcoNtcUx5sTrwq0RD2ly317PF+316t9TTsLm/Bz1d5xv3eMAPsW9qinu\nLyoO1XLNNdcsfPZJI5TDJ5/U3/A/1KWHSVJ2n/DD/vU8XifrtEaoD5Q3KpA2WS3GSy65RNKy396n\ngFOfKEeOSd14qCX/px7r88Bv3NphX6+/ddKb7u8htu6vx9JypY0FwVZqZaE+PHQVsH48dLCOKdCu\nvK3upV6ivEMIYUC2rrxrj+P+WVcE9NZTU9h7CscjK3xSA2qeXhQ14v5B6cRGhPfDVFqA3v+mpvn6\n1Gbqhokq0vKkDAfV4WqUcld17eMPm4AyooC4TleUWFnUEf+n/NxjaXm8hM+cg89u0fgU8fqbdbYT\nbxM+1Z0y4efv+Z+xpDyCBuXIONDDDz8sSXrsscckLU+CcyvWI6CkpmBRtF5v62TVMd3K8uggn7Tn\nKrnXrl1p8+7Ap+3jTdQz56z1RJv0ZyvKO4QQvqRsXHnvJX2rTz/1ng4FwNZHc+vI7ZVXXilJOnr0\nqKTWSz777LMLx+B7VzT8X2qKxWN9N4lfn4+mex3id7v//vt3jkG9MvX72LFjkppKYqQd5eBT76vy\nnitVwCo1xf/w5aIsmcLuUUgoHRQQKqdGZLz00kuSmg/ZfffUCbiS60U4zNk+OK9bHZSReypJX/3q\nVyW1ekKV8xva/n333SdJevTRRyU1a436JDILKGfPiqbs3IM5/f/Qs3RcBXu74rqwOmgTbpVWSwZL\n5L///a8k6d5771343iOTvN3VCBcsA95t+xkzivIOIYQB2bjyXpWo3tNr0hvxG5/hiNKmd6d3vPzy\ny3eOefXVV0tqvR2xrCgWFDjqHRWCAsN3JzWfYS9Zzdx43bifn//7qDXx7VV5U0aUF0mYUAg+o5L6\n5n7UUXVPujMXNeqH+4ylcNttty1871FL3Gu2XCttQZIeeeQRSdITTzwhqbUHj5/nt9SdR09IzVrz\nVKvrhOeGOkAp0t5vv/32nX1vuOEGSa39Ug8eL//4449Lam2FGavu36fMfh+wfOo5fDbinLOS9xJl\nMhWl4+NqlAlrg3eN1CySe+65R1JrO1hnWISexKwX5XUiihuivEMIYUAO1QxLn5FEb+2RFPwfRYii\nuPbaayVJ3/zmN3eOeeONN0pqyhulTW/IOZ988klJrSdmxlpVlPTWnG/OWYTguRlc4Xq8vM/4Q2Xj\nl5Oauvjxj38sqcXiuq/OI3JcUUrzp8nt5ZhA4dx6662Smu/eY4qncuNwzTXeljKiQl1pY8GgpjhW\nb7amR0PNgZcJ1cszwBiP1OoFFYraI2Xu3/72N0nS3//+d0nSgw8+KKk9Z/jKORd1wPeMK3FuaTki\nxdvJHPQWGec77gnXQ9mJaccKo4xf//rXJfVztBw/flyS9J///GfhWNQv1rrPK/A5AtLyONJ+np8o\n7xBCGJCtK+8K6ocRX4+1RTGggvCroTq+9a1vSZK+973v7RyT/7n/k2N4on56UVRnXczAIxE2ubDs\nbsrbFTeWBuXELyc1xUpMMErWc11QB55Yv17DphamqPHU3FP8kfgYXen4vIGp/BZSqwOiNLjXrqr8\n2KjQqqaw2qZiw9eBXwefPf+OtDye4+Mgf/7znyVJ//jHPyS1Nu+zDvGn025uuukmSU3lM2ZSz7kq\nk+O66Slv9+ujtLHAUdF8T9tiHIU6q1YaWSuJUGIf6sXno/is59pWDhKFE+UdQggDsnXlXUdZGdFF\nedMroSL5HoVOVAnREihvesD6W89J7DPAPJKktzyRq8xNKm/wc7oC8xwSKIxaz9ddd93CljoAHwFf\n5afcVB30VDJqyPOrsPWc16hpVGGN4fdcNyg1z6zH//Ej92Zreiw+zDUztx7byyq1slFe/LsPPPCA\nJOmhhx5a2I9xItQnx+aZ8YyWUJde223h7zkVeFW21ANlJpro6aefltTUs+dx4X5ipdRjYq3zG/bh\nM9YYx+Cd03unHKQeorxDCGFA8vIOIYQBOVRuE8wOzBLMMExm/u+TApiQwMBJnX6KWeuDcJ4aFtOb\n3/qyaNLqxFFz4yGBu5ngDKJgFtbpvbhLCPHyCRV+jjmTCu2VXije1PJTnjCIOmB/9qvuItxmnpKB\nLaavu6fcZVP/9kEoX15sP7irzl0znBPXSHWbTP2P3zLgi7uR58pDYqk33JtM4uG5rOWjbjeRQsIH\nQ6s7jGujDTARie95L/Dc49rg/7gd63R5Fl9g6/XBMalXPnuaYWl5Ut1+3EpR3iGEMCBbTwlbBwLo\n7ejJ2PqkByYHELpEyBi9Y+3ZfNCKQU96Te+BUQyoktozcvxNKlBfBm1KBXO91BlLMqE+a3L+m2++\nWVIr61TvvxcFPmWNrEtxoWzrQtJMqGIQivbBICNb7i1JhxhY4lpr8ibUpi/CMZXO09tNXagAhe+L\nRawDX3qP83tytnq/uXbKwm+wwFCb7IfiBo5JOB31+MILL0hanrQjtWACrmPKclgnvWP6wh0e5otF\nyiA4W/YjVJQ2J7WQSp4xDxf1JfJ80YveMoeQxFQhhPAlZ+vKu/oy6R3p6d1H50uW9RbElRanb/M3\nPSiTVfBleSrYupyVf+bvucOdesree3GfeOJq0/1v+LelFgLmaQemrmNVGt+5cYtCasrP/dOUg/aC\n4sbKQjlicVRFzG9oL7QxwuM8RSj7+0LA9brAF7k9CJ60jeulPfhCv/V/njgJa4Nr94WJqWemgRNe\nx5gUFgb1WFWjTx5y9TkHvZSrPLMobPz63COvC08eRRsiTa7U3iG0O0+IhycAC4Z677UDHxeJzzuE\nEL7kbF15V3yE31NqTiVSd19nVeIoUZZ2YkICqoJzoTZQKXyuURrr9F326ClbT4DkSeXdWnG16alC\npaYIfOkv96dOTc5Z5fNetzXikyXq3750HcoG/6/7vFGYqDGUo9TqzyMCULJTqURpc1VtUweouXWq\nTR8D8GnolKNaAvwGSxXfNz5qT6Dl6t4jNlzd9xYo5jqm0tDOgS+fKDVlDb70G599ohr3Fb9+HXOh\nTGy/+93vSmoLYKPufQyO66rRMNThiSTsivIOIYQB2Zry9igKqakeFBS9ki8GzPf43ugde5BEBr8d\nv0Fd+HRylAE9ck22s07f5SrqCPSU8ub6UVhE0eDXxWohoqBGAvjSS74o827JpnpKct2Jh3xMpCpJ\nfK+oO/ZBZfk4B6rJEwZVvzAqinbg6YfZoshoV34OqVlrjDOsc/EOrovrcP8/97aO1XDvuUba8dSc\nAV8ajKgKyso5fWGQqqq5zt4CBHNBm6nx5ih+xi7c2uT5pu1Q9hpdIi1eP+2HY5Keg3kmUwt787te\nKoUTSfAW5R1CCAOytQWIdy6gpPpEseBj85hHPtPzE1FB784ob/0d6gF/nSsUVJEvd8X3PbU995Jf\nq+JAOTd+M/e5enQJI9017pf6oZ5RaXzvs05dDdTyzz2DrudLRx37YhFT/kNPLuWJuaSm0Dzm3ZeA\nox0RrYTFU9U1x6qzc9eFz+jEksQSoJ5qu3Erl996lAj1yDGJqvB5F/zeLYraxqhrTww3J72xMdo2\nyprnAXVMHfA+4LlCkVPGuhyijyH4smf8xueYeGx+/e5E6ifKO4QQBmRjytt7ll5PQy+Njw7FQE/G\naDaKkdhLnw1V1bznMHGfIJ896sSjOup55qKnMqcWY0Zx+5JwjJqjvHu5N/wYnize1cCqpdjmWozB\nfei9ma60C+6tz4r12XMsxIuFxvdSq3vaGHWD2mQ2JxEq1DvtpB6LKBPU2zpzw3gbQSliAfgzITVF\nzbUzw5b64VgobqwKttQFKtXjqVGg1KvUlKqPJW1Sgddr5bnmufel/4D2TNtiabcajUKde34Zvvex\nAE9V3FtGMMo7hBBOEjamvF119PzGPoPSY0JRnygBRsEBhdHLBEgP674szw7n6nOVf2rdizK40q1l\ncp+2L+xK+YiaoJw+TlB/i8r0zGqed2bVIqlzLW/lvtsa9cOoPlaTl8P3Y+Ybse6ed6P+1hW4L27A\nOVHaqNEaT+zLsq0zzptjoSB5VnyJs3o/uHaUIfef6/RZuqh4ttwDj+ahDrA06izeer+kgynMveLZ\nD6XljIrcT49A8qgmt9RredjXjzVlyUJvcZODjJ9FeYcQwoBsLc6bHqf6Mn0JIlQFqhNlSC+K2njs\nscckNR8diktq/nNUgytxjwhwf1X1HbpqmCtPQy/fi8eI+og//kt8c9Qdv+/NJmTrESme82TKB16/\nWzfUrWf3k9q947o9fzZtit9gjeCTRSHWeNwpK4+2xnV4znPabPV5cw98Zt06cOXtEVbE+dfZjlw7\nCpG2QL1QD5690zPvUa8cj7JznGrJgPt554zU4tj1mfXc7pSJ6/IcQa7IfW6A1J4b3ksobh974X3l\n9bZqnsR+iPIOIYQB2bryrr0QPZOvYoOCQW2gqFAQ9I4oiKrSUA2Mfnt+ZnC/7ir/7tz0lC1KwVWR\nL0AMrtBrPXNMVxceTTMV372pepCWM+JJTfVSZu47W5S5j6H4tuLKm/OiLlFXfh88mqFe1xxRJtxv\nV/x+L2sWRlfUnhPIZ5eioHl2fB4A95/9+H9tF36uqbY1h/Xas1x9tSXqZyqvj/ur6wxfn3XL1lf+\n4pgeS74uorxDCGFA8vIOIYQB2ZrbxAPcpeXBRMxbT/XoZnsdTJAWp+266wV8KrgPTGzDRdCb2uwh\njdQFn6krn5JNXbFfXdqK33iYlN+TqTDATdZJL3kSZfTlvTDjpxbp8LCxOnGJusD0pYy4JTiHJ7Xq\nudd8UG6dLgIfyMXF6O6zOmA5NdhNKKBPTPFJOJ7EyRcs6Q3qeVjcJgb7e+3Sp/77Z1/yzxci7w0y\n8vz4wtI+eccXpFjXQCVEeYcQwoBsTXlDr+fxadiozakllValm3SF7aGArrR9u0mVCb3pvZ5cyFOm\nevpcYL9q4UwNpLnS3qbi9nNWNeuLdLhqJhnZ1LF6x/Tp21N14AvYumKr381ZX67qsLwYZKz32wcz\npwYRUd4obCxfP5eHr8Iq62OTrJpE5laQW7CuuH3QUWrWjnsCGLj0Z9Gt+3XVTZR3CCEMyNaVd8UV\nkatlerQpf1lvGbGpnne3JOjbUJk9ppJDuYryOlm1PNnUNP+pcx8Gesm6POmWLwkGe/Gv7nb/T6R9\nzOHXnSqbj4309tmtjUzt59PO99p+DgP+/LjVwXPkU92nLN36W3DLxrdzWfNR3iGEMCCHSnnDKt/k\nKvajdA6zWuhxWC2DbZI62RsjKOS5mLK492rNw6r/b8tyjfIOIYQBOXIy9sYhhDA6Ud4hhDAgeXmH\nEMKA5OUdQggDkpd3CCEMSF7eIYQwIHl5hxDCgOTlHUIIA5KXdwghDEhe3iGEMCB5eYcQwoDk5R1C\nCAOSl3cIIQxIXt4hhDAgeXmHEMKA5OUdQggDkpd3CCEMSF7eIYQwIHl5hxDCgOTlHUIIA5KXdwgh\nDEhe3iGEMCB5eYcQwoDk5R1CCAOSl3cIIQxIXt4hhDAgeXmHEMKA5OUdQggDkpd3CCEMSF7eIYQw\nIHl5hxDCgOTlHUIIA5KXdwghDEhe3iGEMCB5eYcQwoDk5R1CCAOSl3cIIQxIXt4hhDAgeXmHEMKA\n5OUdQggDkpd3CCEMSF7eIYQwIHl5hxDCgOTlHUIIA5KXdwghDEhe3iGEMCB5eYcQwoDk5R1CCAOS\nl3cIIQxIXt4hhDAgeXmHEMKA5OUdQggDkpd3CCEMSF7eIYQwIHl5hxDCgOTlHUIIA5KXdwghDEhe\n3iGEMCB5eYcQwoDk5R1CCAOSl3cIIQxIXt4hhDAgeXmHEMKA5OUdQggDcurcJ/jTn/70xdznOCz8\n/Oc/P7KX/f74xz+eNHXyi1/8Yk918rvf/e6kqZNf//rXe6oTSbr77rtPinr50Y9+lDrpsKpeZn95\nh/XwxRdf7OvzKo4c6beHqe9PdL8QRmEdbXo/z+A6iNskhBAGZGjl/ZWv/H/f471mrxfdTanyG77/\n/PPP13ad64Tr+uyzzxa2H3/8sSTp008/lSR99NFHC9+zn7RcP9TjKaecIkk69dT/bxZnnnmmJOm0\n006TJJ1++ukL+9XfHiY1flDLojKlpjatstbB1HMytZ3Cy06brN8f1vrZrWz7rYsKZfbt1H4HJco7\nhBAG5FArb1Qd0AvyPQqR71GEVTVPKWgU6m69ZA/23YTa9OujPK7AUdi+fe+99yRJ77zzzs4xUeUf\nfvjhwv8++eQTSdI555wjSTp69Kgk6ZJLLpEkXXDBBZKks88+e+dYZ5xxhqTle7FJdlOQJ6KmplQU\n9e6W2mFTmrWM/ty4hcWW77Gw3KqijLQ52otbedJyO4Vt1dN+rQ7K7tu6z5QFQv24RbKq7CdSL1He\nIYQwIIdKedPz07O5rxWF4Mrb/18VwPvvvy+pqUtUAr0ix+IzShVVUXtczoMveBP+3t2UQfU/S9J5\n55238Bl1/cEHH+x8R5288sorklqdUGasElfx/H/b4wFTPnuvE9rPlA+/1h3HpGyUla3XgddFrZNt\nqHG3PmvZKD/l9vL7+Mkbb7yx8Bnr6txzz134fNZZZ0larpP6N+3Px142Zb3uprDd6vB3jluW9W+O\nMfW87NZm/G9pf20myjuEEAYkL+8QQhiQrbtNqnmHqcIWEwZTws0TTDDMYX5XTTF3BeA+wXXAsd5+\n+21Jzb2AuccxJeniiy+WJF1++eWSmhnpA6vrhPpxF42bW5iwDCoCZl81ad966y1J0mWXXSZJuvrq\nqyW1snMOBi7Z+n3pXdcmmXKTsOW6uU8XXnjhwvf13k65TWgHPvDLlv/X+qWtuYtgDrwO3Nz3v6VW\nVtxmzzzzjCTp6aefliQ9++yzklqZgXq85ZZbJEnXX3+9pDawXesT/BmmTjbhapR2byPUjQ/S8k7x\nbd3Xy8a7BDck9efvFLb1b3ep7KXNRHmHEMKAbE1508NVFefK23u0119/XVJTjvRWrrQq559//sL5\n6JX9mCipd999V1IbtKEXlVrIHNfFwEUNnTsIPRXrisEH59iXa0BJsEUx1J78oosukiRdeeWVC+Xw\nUMvdwhDr36sG8NYNytbPwWfKw71EAVEHrr6kVnb2oX5pi5SPeqVd+P7SsqqcQ4H7wLUryVo2v0fU\nx/PPPy+pPQPUARamD1h6YMCbb765UK5aB654sQx933XUydQAttTqgefAnxO/7yjhl156SZL03HPP\nSWpWidSUNGW67rrrJDVLlvcE/+ecXANWf732/SjunXLuec8QQgiHhq0r755vjt6QHgplwBZVjJIg\nPA7ljbKUmj+OfT3c0P3oKHDUyfHjx3eORQ/K+fB9r5uqJDzMy0Oc2NLLo464Vn5XfbLsy/VjnaAM\nUFyuXPnMVlpWX1WVH4SpyUhSU0e0D67P/YfcY+qC8QCUJT7wev3uy/Z6dp+o35e6r5dhDuXt0J6r\nvxpLFcXI/7j2a6+9VlJrB9QT+3P9PJevvvqqpFYXjBfRFupvgOeQup9TcVerg3vPlrbvipv7jjVy\n3333SZLuueceSdLdd9+9c0xUOfX205/+VJJ0xx13SJJuu+02Se35oj7dOpKW/eT7sdKivEMIYUC2\nrrx7fkeg5/dIEA+0R22iiPE5Sa3HRQn0/J3SctQJx66TW1As7LPuKAL3Y0rLym+3adkemePRNlKr\nJ1egHpni56QO67FQoPzmoFEEU/72OkKP4uY+oASxCDz1AdeE6sJCq1YC1hzH4DdXXHGFJOmqq66S\n1OqMNsaxqpryyWabgLJQDiJJpFYmFDjgo3W/MGXlueIZqPddkp588klJ7bl47bXXdv7HuTgmbYm2\nts4Iram0GdKyZUSboEy0W+oLHzd1RoTND37wg51j0jYZM7r11lsltbbh7w63jnsTuk6EKO8QQhiQ\nrcd5V18m4P+h96ZHR3Hxf3pHHxWvKoie32N7fV+fJs+2Kj5UjU+PX7cCX5UExyMHavyp1MpH3bF/\njcFFNaIqUAScaypypKe854oq4bjc6970fq7D45g99YFP/8Y3W9Up6t0jLzzWGRXl4zO1fjfp66YO\naKeUCd9t/Y7nxcc2psZLvH1Qb+zHOakz6lBq98h9uJ7Yax0KfFXiMU+eRRm4duoExU0ZeF/88Ic/\nlLT4nFEfWLBTap6242Wsz8yJxHdDlHcIIQzI1pQ3PUyNowZ6LNQjI8EoX3ydHn3C72rPRi/JSDqq\nk96Q2YXHjh2TJH3ta1+T1JTLo48+unMsvqNX3oRvk7L4LC0UF9uqTKVWFyhvRsClpqhQkT7SjWqi\nXqlnVAv1X3+zLgW+nwT2bj15BIGrO66R66ce6jGoXy8P5dwtaVfvfHPC+f2esa1/uxqmrKjOSy+9\nVFLz5bpVx+94HmkPvXbhcyF8DoHfq3Xg4z1Sey54l3CNfKbt08bxkX/jG9+QJN18880L5ajHR7Vj\n2RCFwjl5TxB1QhRc71g9D8RuRHmHEMKAbF15V8XiPTs9mCd7p9d84oknJDWlQ0+IcpCWR3rxU3m8\nJ9fBSDt5G4jllKTHHntMUuu93de6Lqpi87J77K3HI/MZtYXFwbZet89UdT8gSuKpp55a+L5eH+rV\nZ8WeKFNJ86vfEEUzldfD93Ofvi8ZJzX1yb4vv/zywmeiJDydas9nu4lZpkCZeHZ6qYJ9+ToUNjz0\n0EOS2n3GavvOd74jqbUt1Cl14+2m1j/1RLvznDG9yJCD4r7neq2oY8pAG3AL4IYbbpDUFDfXX+c2\ncExywGCdc2zKyG+pH88ZVL/b6xJqlSjvEEIYkK0p754qmVoUwP1B9HhkQ0Md0Wv2ok08xpLe0BPS\n0+OhIupszWuuuUZSU2mbSCjP9XqkjdcJ18JoOf/Hp0/5paZCPH+L+7ypZ8rLflXVca98LOGgTC3I\nUc/PvfVlvVCW5JxANXv2yTre8u9//1tSU22eAwMVhaXBtjfDck7l7cqMMqPmsDprdklvK9yr//3v\nfwtb1CX1xfV7PbPlnL5MXj0Gli7155bLOp4dz1hYfe8+XkPb9hxBPOe0GeqVZ6XmNnn88cclSS+8\n8IKk9h7CYqHteK6l3vjQQSLWorxDCGFADpXyhqnMePia6D35nl4elVnVssdi4uflGK6e6XHpJas6\nY19X6+uil51taiklt0581iFqCBVa806gqFElbp2gwFDxbLmWatm4RbMu5e2+23pOzyKJ0kHdMU7B\n+AWj/Z4rb7GHAAALCUlEQVSPpC7KTBmZJUe9elY8j2yhvVTltElftyvvnqXiOTW4VnzcRFG5cr33\n3nslNTWKD5xZp9RzL8rLn5WpPPTrxCOzpHYffTY218X7gTJ5rnPGfRhfq99Rb+xL/dKGqB/P/1/b\nci93z16J8g4hhAHJyzuEEAZkdrfJboN61czGhGHggy2mIS4PwvhuuukmSS3BDhNsWKZJaqaKD1Ti\nMvDV5DHBPQSr7js3qwZxuE5cGLg4MHUxEyk3x3rkkUd2jsHgC2XDveCJqXDF+MSg3vJh6xq8dbeL\np2Kt3/l0bsxVtrjPfPCxd29xK3kYISY45cN9wn5e/l4Z5gTz3+ughoZyzbhPKKunWKB94wbwCW4+\niccn2tVwOp9U5mG/65icM5WUbVVYqS/k4QOrPukJavpnzsOAJXVPeOG3v/1tSe09xbF7E5MO8k6J\n8g4hhAGZXXnvpsR6ypueDGWFmkCJo7D9d3xfJ+mgHnxgj2RXDLag5uh5USk1bahPQV83vQkpvvSU\nJ6hCcVMeVBK9PN+jEqSmuFBSU8rAF4tFudYQKE/sta6wSQ+hquGJHvLnE7D4rasoH+Stg1CEyzGI\niZXHebFkqBOO5QvXVnqqfF14EiYf3K8pj31yki9m4guj0IYIjfXp3Tyf/jtP6SwtD/x66oA5knbV\naf08N7w7qANPGsV9p06439RBPebDDz8sqbUZBjsZ0L399tsltfa4KgjCLbj9EOUdQggDsvWUsD3o\npenp8cmS9NyTMrnvrk65RXX51Hp6WraubOmRe9P356YqGJ/2jtL2FJcobFedlK/6qT3Ui98QQlgn\nOUjLCqIuH+bLra2bnjXiizKzj/up+Y377FGehMJJTYVTb1h77vdFwVGfnh5YWlaVcypvV7S+hF39\nmzbCEl//+te/JLXxnzvvvFOS9P3vf19SU53cY/cXrwr/83GFTUxo8yXqpOWQTtqCq2DqiLK6L7ze\nX37LPlgmN954o6T23vJnwuvN/94vUd4hhDAgW1feVZWgAAiCp7f0KBS+dwWOyq6qA+WJ6uAcvliB\n+1jx69bFGHy68LrpLQ3nEysoB+qSsr744ouSmpJAXaOeajIiJl2gMpjey5ZoA194AD9gjWTg77nS\noPpSbFKrH1QxWx8PcH8mivvBBx+U1HyWUqtP1BTl8igOv/e9SRacfw7FDT5Jh7J6ilOpTf3/5z//\nKalFG7G4LsfAV8vYEXXgkSL1+ZKWF37ufccz6hbCOuHYNRUEbZcy+pgAvm/K5ukPuIe1PmlHPEds\neU9NJeHyBSn87/0S5R1CCAOy9QWIaw+Mr5XEL6jMO+64Q1JTkd7D+QKoNRE9Cpt9UGGeJMZ7S/ar\naUPpredSVD3/LsrPrRB8ru7LRxUxPR7/G8pRapYN9YUSY4o46sSXkHNVKjWVQ73NNS7QS3HKuV1N\nefoC6oYFc/Fv1yXOsOpoYygvvq9llpYXRq4RBL4Y8xwRFZ4G16NgaiKlBx54QFJLXUrkEb5afNw/\n+clPJLU64Jjuu2WcxRd6ruWkbj1iyX3g64Tnpo7v+BJvPp7mz7+nY+C+1lQK/kzxuZfwrH7vy/HV\n8yYxVQghnCRsTXn3VCY9FooQJchnZk4ePXpU0vKCDj6SXPcBzod/3ONlOZYna6/7bBKu11N++vgA\n+/E9C0cw+7T67n0mHRaPjy2g+lEWHucrTc9yWxe9RWU9nedUmlHKjGryBFvVN0q9Ul/UgUcO+IK/\nvdjdOWKYHY8ymVqaq+5D2YiKYGk8ZgTi68bq8GfEUwajrnvPhS8n5zHgm4iBl5bVuCe98zh9f6fw\n7qlzG3xh5vo8SMvzQfzdUtuKWy/7Ico7hBAGZOvKu/rT6Mno2VCVxKUSh0zkBL0hStxjcaXWs6FI\nUSY++u0zLFEZtVf1RVTXTW80mu+oJ6JIuH5X0cz+In3uX/7yF0mLURKoDeqZOHr85J5K1WObq2Xj\n0QTrqhtXc9VCc5UyZT2hgLAsuLe+PJq0PIvXZ/sC5fOFQlZZe3Pg9UOZfeFpqaVKxnLlvpKLw5ey\nc8vBUwVjyaxS3uDW0Jyx7z08ttrj890X7j55vq9zG3hXeCSUj/t4PXo6596++7Huo7xDCGFAth7n\nXUH9HDt2bOF7VCW+NlS052ygh6sjw/4b1JgvfOp+KtRG3W/ubHG9RZndP4tiYHakR4TwPREGRFrU\nEXiynRGfir8XCwZ1z7E9e1y9Pv6eewGCenz3HXqsrmdF5P+em6Umxad+3IKoZZVaO5nKuChtRlW6\nqkPB0V56EUEeqURUFs/G8ePHF/bjd7R7LBjmFFBmXxCi/gaLcRPjAFCtILeMXOn6/fa5ApSt+ux9\n7oVb65750zMp9qJNTuTdEuUdQggDsnXlXXshemmfAUhcal3KS2qqiAVy8ZHXXpLj08O6mqTnm/KF\nzZWzo4fnHJea2vGlvzzHMrnNf/azn0lazppYe3YUFcea8t3xW8/K11MOc48D1HN6ju1qadV93R9P\nuXt5lSmr57ZxH63PNuzVySbUpc8u5h7iv654RAWRSPfff7+kFv/NsbDEPFeHq1ePQqpWii+Ptwlf\nd2+2q+c/chXsFrhHD/VmPVOfHuPvFqAvD9eLnOoto7dXorxDCGFADtUCxK56UJf0hvR4qFH3aXrk\niNSUiCvWqYVzV+Um9utcd9x3Tzm4InQfM9dNlATluuGGGxaOVevbI1V88WKULCrTfXe9aJh148et\nPmXPseEzGvnMMWg3Hi1TIzIok+e0oW74rSu4g8TpnghuMVIG/NRYYjXiinqiLFw7EUlkV+SzK2+i\nVRgLIV7clXjvefAZt3PWU0/d+/iHj2G4T97nn9RxET8PzwtjBnz2lYp8Nmwv3/iJEOUdQggDkpd3\nCCEMyNbdJr3JFz7IwcClu0/cHORYddDGkyntNbnUqv/PPU2+N/iF+8QHEwkF9KXjpiZHSM3U9+RK\nU4NxPiGlHmvuEEHOVd0mvrwZbh4P2/MJTtQJ5av1jPlKO2FL+Kov9+Yug1onmwwVdDdObxk0H6Rn\ncJt9mCZPPfqUdn7P5C0+0+Z43qqLYZOTctyN2Qsr9bbsn2k7tANPN1xDbX1inKdjdvfvqgWIfRBz\nP/UU5R1CCAOyNeXdmwrui9h6uB49uwfYe09XBwT8mFOJlDwkbVMDULvhynZqIoonKuqFJfkxvcw+\nYDoVBji32t4NH1TySUTUiZePbW/xAF/qy8/hU8hXKcpNtB23SLBCsUrqM+ApfVHQJKJyFe8hb/6s\n+jPSq4ttPj+9xQ5cYfvCLQw2VgtPWk4RW4/Jb6knt1RoU/y2p+IPkqAryjuEEAZk64mpKq7oXD16\nWE3twep+vUk6U4p7Sl3OPRV+v/h1Ti0ndSL++N3Sum5baUuL935qQQJ8i6hP92v6dObKlLXnCtf9\nlgeZ3rwOvF24dSE1VempdMHbjCtpP8deJt5sQ3mvuh7umyclwxrxNuLhp/UZ8OnxU4sse/tku64x\ns8P1hgohhLAntj49vsdBlV5vKa7DoB7XyZetPPsBRcMWVe7pW6csiJ615WrIfd69SIHDgF/POpah\nm1Lio1Cvf0rluirey7F2Y9PpbqO8QwhhQI6M1quGEEKI8g4hhCHJyzuEEAYkL+8QQhiQvLxDCGFA\n8vIOIYQBycs7hBAGJC/vEEIYkLy8QwhhQPLyDiGEAcnLO4QQBiQv7xBCGJC8vEMIYUDy8g4hhAHJ\nyzuEEAYkL+8QQhiQvLxDCGFA8vIOIYQBycs7hBAGJC/vEEIYkLy8QwhhQP4Pinba/Yr4skkAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ab6036850>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Final weights\n",
    "plot_digits(W.T, numcols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.3. Visualize the 8 hardest and 8 easiest digits (10 points)\n",
    "Visualize the 8 digits in the validation set with the highest probability of the true class label under the model.\n",
    "Also plot the 8 digits that were assigned the lowest probability.\n",
    "Ask yourself if these results make sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hardest 8 Examples\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAADiCAYAAACFkJabAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF3VJREFUeJzt3Xm0VWUZx/EfoiGIU5DScgpFnNBcLlIZaqmlqGlLRcUK\nG0wlwEzJYYkaBrkUjNREJA1LK00srbBMbeGA4VQ4rURXkAOWESWkRE5of7Ce5z7Hs++Z97nnPff7\n+Ydnvefsffbd99yHd79jj/fee08AgDRs0NUXAACoHEkbABJC0gaAhJC0ASAhJG0ASAhJGwASQtIG\ngISQtAEgISRtAEjIhk34DKZcVqZHDcdwbytTy72VuL+V4N7mq+j+UtMGgISQtAEgISRtAEgISRsA\nEkLSBoCEkLQBICHNGPLXpRYtWiRJ+va3v+1lF1xwgcdDhgyRJG222WbNvbBu4O2335YkzZgxw8su\nvPDCksfYphzLli3zsh133DGHqwPSRE0bABJC0gaAhPRowh6RXTrzqWfPnpKkDTbo+P/p3Xff9fiw\nww6TJG2//fZedvbZZ3s8cODAvC/RtN2MyJkzZ0qSzjnnnKqP3WGHHTyeP3++x3vssUctl8Ksvfxw\nb/PFjEgASFnb17R79Fj/H1VnNW0rzyqTpL/85S+SmtIZ1hY17Tlz5ng8adIkSdKbb75Z1zmvvPJK\nj0877bRaTkFtMD/c23xR0waAlJG0ASAhbT9O25o6YpPHH/7wB4+t+WT48OFFx0gdnWlXX311rteZ\nstdee83jBQsWeFxPs0j//v093nPPPWs+D2DWrl3r8TXXXCNJmjdvnpc99thjRcfE5mPLFZL03e9+\nV5L01a9+1cs23njjxl1sCdS0ASAhJG0ASEi3HD2ybt26ovf99a9/9fjggw8uKp8+fbqX1TLuuALJ\njh5ZsWKFx2PHjvU4NpVksbHYl1xySdFr2267rccjRoyo9xJbZoSD/b1lfQfLueeeezxeunRpyfcu\nXrzY4x/96EdFr1966aWSpLPOOsvLbE5DlVrm3kbWNPenP/3Jyw488ECP33nnnYZ8zqabburxgw8+\nKKljaYwGYfQIAKSs7WvaWTMib7zxRo8/+9nPFh0zceJEj23ccTzeFkJqsGRr2rE2c/jhh3v8r3/9\nq+RxI0eOlCTdf//9+VxYh5apDb7xxhuSpFGjRnnZwoULG/0xFXv00Uc9Hjp0aC2naJl7G918882S\nCp/88mZPhA888EAjT0tNGwBSRtIGgIS0/TjtQw45RJJ09913e9nnPvc5j+3x6eKLL/ayOBXbprfb\nwlIoFjtjtthiC4/LNY8ceuihuV1Tq6ulI7Iacczwhhuu/zNfs2ZN0ftih1xnY5JTEcdhT5s2raJj\n4nf3iCOO8Hjw4MGSCheSs7X5Jen111+XVDjO++mnn5Ykvfjii14WFz5rFGraAJCQtq9p33LLLZIK\na9d33nmnx9bBeP755xeVSR017JtuuinX60yZ1UqkwtrKFVdcUfK47rhbkNWADzroIC+LnYEf/vCH\ni46xRbJ22WWXij9n55139thm782dO7fofbFWuO+++3qcYk07DoN87rnnKjom3qef/OQnJd/7pS99\nyeO33npLUmFN22YGf+c73/Gyq666qqLrqAY1bQBICEkbABLS9s0jJq6XHTtcrDyrTOroiOiOj/K1\niAs92azGl19+OfO9NrN0k0028bL4CNrOJk+e7PGnP/1pj/fbb7+GnN8e3yWpV69enb5v11139Tg2\nC6Zo99139/j000+XVNj0ZIMSouOOO67h1/H5z3++4eeM0v4tAUA3Q9IGgIS0/TR2m1YdF9upZrsx\nax6JY4pzWlu76dPY7WeOPe277bZbPacsYOOCjz/+eC+76667it4XlxKYPXu2pIY3R7XkVOtGi+uX\njxs3zuO4bIOxESmXX365l7XTglH23Y5/1zZevRGefPJJSdI+++xT9FrMD3G97RoxjR0AUtY2Ne04\nrjTWlO3ni6/HWoiVx/sQa4PPP/980fHWoWFjwKWG1AybXtO2GoGN45Wko446yuMpU6ZIqv9nmz9/\nfub5s9jiU3vvvXddn/k+LVkbbLRf/OIXHmd1sMWlSX/3u99JkjbaaKN6PzaZexs7Z23Rt969e3tZ\nVq08Lg5nMx6ljg7jrEELTzzxhJdtvvnm9V42NW0ASBlJGwASkvw4bVs3NzaJxNgeX2ynDkk6++yz\nS57TmkQkadCgQUXntMWn4tT4O+64o+pr72p2zS+88IKXxannNt40TvWNC+zkwZYNWL58uZc1sgOp\nXdi63FLHYmezZs3KfO9WW20lSfrmN7/pZamPyS7n8ccfl1T4PYrTy21nnzFjxniZLQIldTQNPvLI\nI152yimnlPxMazJpQJNISe39mwOANkPSBoCEJDN6JG68u9NOO3lsU4C32247L4uPfvWOqbbPnTlz\nppddc801krJHlEgdo0qqHHXR9NEj119/vaTyj30HHHCAx2eccYbHtgFyXLs5SzWjR0wcc9yA5pFk\nRjiUEjfrjU189957b9F7t956a49tOzNr6muwlrm3cds7+3tcvXp1oz+mgG2ZJ3U0Nza4CZHRIwCQ\nspavaY8fP15S4YzGuO6tLRLTzAWdsjYLjuM1rTOtys7Jpte0bQxquZpyZ2xRnjjrKz592Drb1LTr\nc9lll0mSpk+f7mWvvvpq0fsGDBjgcdwsOXYk56Bl7m0cm37bbbc1+vQFvvKVr0gqnOPQt2/fPD6K\nmjYApIykDQAJaakBsDaWOm79Zc0OsSli//33b+6Fvc+qVaskSSeccIKX2bRgqXA7s1YTN3ctt/Fu\nOTYu+Ic//GHm60OGDJFUOBYWlbEmEamjWSSrSSSKnZM5N4l0S0OHDvXYNv/uivHu1LQBICEkbQBI\nSEs1j1izSNYjR5yG3tWsWSSOaOlslcFWY2OzJenMM8+s61zWdBWn/0YPPfRQwb8oz6ZalxspEi1Y\nsECSNGLEiPwuDDr66KM97sq/8dbNLgCAIi01Ttv+98qaaRgXuxk2bFjFH15LLW/atGkeWwdjvKas\nNbrj2PG4E0sVmjJOe+XKlR6PGjVKUscuHK2k3cdpx7WdYweifY+sszuaNGmSx3Fmqs1+bMDa2LVo\nmXt73333eTxhwgRJ0iuvvJL53mOOOUZS4byBj3/84x7vu+++kqRly5Z5Waxd24JdNe72Uw3GaQNA\nykjaAJCQlmoeKTU9PJadfPLJHpfbmDdrPGU1G/tmfX4ta3RXoOnT2M8991xJhesMd4WpU6d6bNcU\nHztjM1SNWuYR3ppFJk6c6GVz584teYxNT4/bWNka2S2gZe5tIx155JGSpN/+9reZr9vvkeYRAEBJ\nJG0ASEhLNY/Yin62nZeUvRt6vOas3dTLjfTIeu8OO+zgZbZKX2fH17tGdyea3jyybt06SYUjEWbP\nnl3PKauyzTbbSCpcDXGvvfbK46O69BE+bll14oknSpKWLl1a8pj4O7noooskNXclyyq0ZfPIDTfc\nIEk66aSTMl+neQQAUJGWmhFpO8LEjXVXrFghKf+aduzYGThwYB0/RTqslhBn38UOstGjR0uSnn32\n2czjbSZe3CWlGr1795ZU+JTTjnr16uWxdSrGmna/fv08tg1lx40b52VdNP66W8t7A+t6UNMGgISQ\ntAEgIS3VEdnNNb0jshvp0s6yuIa5dSZuueWWXvbwww97nOA62G3TEbl27VqPbcPezpZ4oCMSAFCR\nluqIBFIXF4KyDY/jQkYzZsyQVLhQ0aBBg5pzcShiCz9JHcMxpewati0iJTVklm7NqGkDQEJI2gCQ\nEDoiWwcdkfnJtbNsyZIlHsc1mW0mo60JL0k77bSTpKZ0YDVLl3ZELl++3OO48FbWXIu4Ibh1/t56\n661eVm5D7oULF3o8fPjw6i+2NnREAkDKSNoAkBCaR1oHzSP5yfUR/plnnvF4yJAhHvfp06fo9e23\n377GS2lZbTNOu0XRPAIAKaOm3TqoaeeH2mB+uLf5oqYNACkjaQNAQprRPAIAaBBq2gCQEJI2ACSE\npA0ACSFpA0BCSNoAkBCSNgAkhKQNAAkhaQNAQkjaAJAQkjYAJISkDQAJIWkDQEJI2gCQEJI2ACSE\npA0ACSFpA0BCSNoAkBCSNgAkhKQNAAkhaQNAQkjaAJAQkjYAJISkDQAJIWkDQEJI2gCQEJI2ACSE\npA0ACSFpA0BCSNoAkBCSNgAkhKQNAAkhaQNAQkjaAJAQkjYAJISkDQAJIWkDQEJI2gCQEJI2ACSE\npA0ACSFpA0BCSNoAkBCSNgAkhKQNAAkhaQNAQkjaAJAQkjYAJISkDQAJIWkDQEI2bMJnvNeEz2gH\nPWo4hntbmVrurcT9rQT3Nl9F95eaNgAkhKQNAAkhaQNAQkjaAJAQkjYAJISkDQAJacaQv5a2YsUK\nSdJ+++3nZcuXLy9637vvvuvxBhus/79u5cqVXvbBD34wr0sEAEdNGwASQtIGgIR0y+YRaxKRpOOO\nO06S9Le//c3LrPmjM+VeB9A9PPTQQx7Pnj1bkvStb33Ly3bccceGfybZBwAS0i1r2v/4xz88jv9T\nlrLnnnt63KPH+uUAevbs2dgLA0qIHeTz58/3+Bvf+IYk6Y033vAy+45K0oQJEyRJU6dO9TI6zqu3\nePFijydPnixJ+v3vf+9lNljh9ttv97I1a9Y0/DqoaQNAQkjaAJCQbtk8Uosnnniiqy+hZb399tse\n/+c///H4qquuklT4WD5v3jyPrRMYpd15552SpNGjR3tZbAoxm2++ucf9+vXz+KabbpIk3XHHHV52\nxhlnSJJOP/10L0uxg/3CCy/0OH4PL7300rrOu3btWknSD37wAy8766yzPH7nnXc6PfbQQw+t67PL\nSe+3BADdGEkbABLSLZtHJk2aVPF758yZk+OVpM16060nXZLuuecejwcPHiyp8LH7nHPO8XibbbaR\nJA0fPjzX60zJww8/LEm6/PLLvcxGI8RH8q9//esen3rqqZKk/v37e9mHPvQhj1999VVJ0qhRo7zM\n/gYOP/xwL7PfVwpsBNjVV1/tZdddd13V5/njH//o8aJFizy2sdarVq2q+FwbbbRRwbF5oaYNAAnp\nljXtBQsWeJzV+RJrHIccckhTrqkVvffe+m38nnrqKS8777zzPLZatb1PKnwyGTt2rCTptNNO87Ib\nbrjBY+sEi7Wd7uh73/uex2eeeaakwnvat29fSYVPiFOmTPE4jsnOYmOyf/3rX3uZPeXcd999XpZS\nTfuRRx6RJK1evdrLfv7zn3ts96wzM2fOlCQtXLjQy958882i9x177LEeX3DBBR5fe+21kjpmQUpS\nnz59JEl77LFH+R+gDtS0ASAhJG0ASEi3bB6JTSJZzSNf+MIXPN5uu+2ack2tIj4izpgxQ5J00UUX\nlTzm4IMP9viUU04pen3u3Lke77zzzh7bI2Z3YlOdv//973tZbPawZpGRI0d62axZsyRJe+21V12f\nPWDAAI/HjBlT17m6mn2PrElCkm655ZbMuFKx89aWBrDx7O8XF5gzn/jEJ6r+zFpQ0waAhHTLmjY6\nF2d+lqthjxgxQlJ1Q61izag7shr2xIkTM18/6aSTJElXXnmll22yySYN+exf/epXHltN9MADD2zI\nuZtt9913lyS99NJLXmadi1LHcNQlS5Z4WRzeOGjQIEnSiSee6GUf+MAHPI6zS03s9LR7GRfeqqV2\nXwtq2gCQEJI2ACSE5hFUJc5A+/KXvyxJ6tWrV8lj4ky+Zj1CtpI4qy4ucGRiE4UtstW7d++GX8dz\nzz1XVPbggw96bDMrUxKbJy6++GKPbUGt2LGe1eRRjazZ0dOmTfN44403ruv8laKmDQAJIWkDQEJo\nHkGBXXfd1eO77rpLkvTYY495WextL9csYuLiRrYgklT/msepiGtW2+JN8VE6Pnbn0Sxi4nrTZujQ\nobl9Xley+1tvk8Wf//xnj+MiXnbeYcOG1XX+WlDTBoCE9IgL0+Skrg+wBV0OOOCAzNet/LbbbvOy\nch0OcYGdrBmRcaeV888/v9JLrVfpVX+y5f7Lq4d1Bg0cONDL/vnPfxa9bkta5qiWeyvVeX+fffZZ\nSdKQIUO8zGZE/uY3v/Gyww47rJ6PKcs64/bff38vs9/D0qVLvazGWn6X3NtmGT9+vMdxFuuRRx4p\nqXDse06K7i81bQBICEkbABLS8h2RNt23s01HbT3ge++918uOOuqoovfFtYTLLRiV4ganebCms3//\n+99eFpuhrDPm+OOPzzz+ox/9qKTCJhFbiEeSNtyw5b9+VYsb7h5zzDGSOppEJGmfffaRVLiLTN5s\nc9onn3zSy+x3lmfHZ8reeustSdLKlSu9LHZq5r07TSlkJwBICEkbABLSNs+nf//73z1et26dxz17\n9pRUuBURKmPjerfeeuuS77Pp7J3ZbLPNPI4b+5bbJitFa9eu9dhGj0S33nqrpPyb4GKTVlwx0MQR\nLShmY7Jjc2BcC37vvfdu+jUZatoAkJCWr2nbDLzJkyeXfF+cdXfCCSd4HBeUQdewjjBJ6t+/fxde\nSf4uu+yyorJYK7MNdfMQOzxtXW6pYyz2Lrvs4mVxlibWizNGH3300aLXv/a1rzXzcjpFTRsAEkLS\nBoCEtHzziD3mlWseieIYSuuEiWvt/vSnP23Q1aESW265ZVdfQq5sTK+UvV74FVdc4XHc0qpRbL3y\n66+/3svmz59f9L7Ro0d7HDuHsV4czHD77bdLKlwUbfjw4U2/pizUtAEgISRtAEhIyzeP9OvXT5J0\n8803e9mYMWM8jj3mZtasWR6vWbNGkvTCCy+UPCYq93p3YeOo4xrbr7zyisc2MueTn/ykl33xi1/0\neNy4cZKksWPHetkzzzzj8RZbbNHgK+4a8bsVYzNy5MiGf2b8jp577rmSCtd7jmxM+LHHHtvw62gn\nNoY+it99W4Kgq1HTBoCEtHxN28ROgI997GMe264qnc0wu/HGG4vKWDCqMrbOddy55n//+5/H1knT\nt2/fzOMfeOABSdKIESO8LHaWTZo0qXEX2w3EGZZWu5ayOx2j8847T1LXzuJrVcuXL/c4blptjjji\nCI9ff/31kueyJ9PO/h4ahewEAAkhaQNAQlp+u7Es9rgnSTNmzJBUXZNG7MSx4+L42bjAzsknn1zz\ndVap7bYbM7YmuiRdd911Hv/3v/+VlP52Yy+++KLHgwcP9timRcep7dYZ2KdPHy/72c9+5rGNuX7+\n+ee97Je//KWkjsXPJGnKlCke2/Rqu5/vZ5sJ59Txm/t2Y3YfszYmXrZsmcfz5s2r+iJ+/OMfe/zS\nSy9VfJzlzbjomW1zuGrVqqqvowS2GwOAlJG0ASAhSTaPRDY9/e677/ayRYsWlTwmq3kkri/8+OOP\nN/ISK9UtmkfmzJnj8YoVKyQ1ZeW/pu0YHtdfti29ahn3H7dis93a4/yDbbfd1uMBAwZIKtwaK47M\nmT59uqTC5pUGyv3eHn300ZKasvN5kdjsuummmxa9HuclHHTQQZKkz3zmM428BJpHACBlyde0Tdw8\n9qmnnvLYailRVk37/vvv97IuWhimbWva1157rcfjx4/32DpsmrB4UdNq2tGSJUskSVOnTu04YZm/\nN1sL3mptUvb9Wbx4scdDhw6VVLjxbJy5ah1kOcn93trPn9URWSvb3co6fiXpIx/5iMe26NxWW23l\nZc3cjDmgpg0AKSNpA0BCkpnGXk58jPnUpz7lcSMfqVCbp59+OrPcOpZsS7l2s9tuu0kqXOwsT7HD\nN+cmkaZ67bXXGn7O1atXS5KGDRvmZaeeeqrHrfydpKYNAAlpm5o2Wo918rz88suZr7MEbu2sFi9J\nEyZMkNTRIYnybHaodRanhJo2ACSEpA0ACWmbcdptoC3GacdZeTZGvrMZptYZlDXTrMG6ZJx2N8G9\nzRfjtAEgZSRtAEgIzSOtoy2aR1oUj/D54d7mi+YRAEgZSRsAEkLSBoCEkLQBICHN6IgEADQINW0A\nSAhJGwASQtIGgISQtAEgISRtAEgISRsAEkLSBoCEkLQBICEkbQBICEkbABJC0gaAhJC0ASAhJG0A\nSAhJGwASQtIGgISQtAEgISRtAEgISRsAEkLSBoCEkLQBICH/B+Vs9GTA/QERAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ab5771c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Easiest 8 Examples\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAADiCAYAAACFkJabAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGIpJREFUeJzt3Xm4VVUZx/EfDqiQMSpq+kSmgRimaZqJIoIkJUFlagw5\nFYaZqPhQaqilDYZhAilQqQ+KYJCkJZA9oKApJpCmOCUICEkCoTI52x8+673v8e57zz7zXvd+P3+9\nzzpnD2zveV3TXqvF+++/LwBAHHao9Q0AANIjaQNAREjaABARkjYARISkDQARIWkDQERI2gAQEZI2\nAESEpA0AEdmpCtfglct0WhRxDM82nWKercTzTYNnW1n1ni81bQCICEkbACJC0gaAiJC0ASAiJG0A\niAhJGwAiUo0pf2jiVq9ebfFHPvIRi3fY4YM6waZNm6ysU6dOFrdq1aoKdxen7du3Wzxs2DBJ0rZt\n26zMP/OlS5dKki666CIr8895xIgRkqSWLVtW5mYzbOvWrYnlu+66qyRp5cqViZ/vu+++kqRddtml\nIvdVCmraABARkjYARKRFFfaI5M2ndKJ4I/KJJ56w+Nxzz5UkLV682Mratm1r8Y477ihJWr9+vZUN\nGTLE4ltuuSXnexWU6bf2fBP+pZdekiT17NnTyjZu3PjBzbjfaosW6f9JJ598siRp+vTpVha6B8og\nM8/20Ucftfj555+XJC1atMjK/DPr1q2bpLquow8bPXq0JOmKK64o920WijciASBm1LSzI3M17dde\ne02SdPPNN1vZj370I4v9YFkxQm3y/vvvL+k8KWSmNhiMGzfO4tDikKR169ZJym2dBN27d7c41BS9\nOXPmWPz6669bHH7j8+fPtzJfky9RZp7tOeecY/GUKVPKcs6hQ4cmxr169SrL+VOgpg0AMSNpA0BE\nmmT3yNtvv23xe++9Jym3abhs2TKLu3TpIkl68803rcw3fdq1ayepbs5xBWWie+Tdd9+1+Pzzz5ck\nTZo0KfG7PXr0kCR99atfTfz8gQcekCTNnj078fzhmT788MNWduSRRxZx13llpgkf/g5PPPFEK0sa\nVDzooIMsvueeeyTln+PuBxr9gG/4jZ966qlWNm3atILvvQE1ebbhNz558mQr+8EPfmCx/z2XIuQP\nSdp9990tnjFjhiSpd+/eVlahHEH3CADEjKQNABGJvnskNF/mzZtnZSNHjrT4qaeeKun8YR5nhw4d\n6pX55lIZZKJ7xHcj9enTp97nu+22m8XPPfecpLpXfhty3XXXWTxq1Kh6n++1114Wv/DCC5LK/op7\nZrpH3nrrLUm5z9F3jxx44IGSpL///e9W1r59+1TnDt1RUu5/u/AbP+aYY6xs4cKFBdx1o2rybJcv\nXy5J6tq1aymnyct3jyR1f/hZOv6/aRnRPQIAMYtywSj/f78777xTkjR48OCKXOuGG26oV/anP/1J\nUu784jLXumsm/NsaEhYvkvLXsIMLLrjA4gULFlh87733SqqbmyxJ77zzTqpzxios2vT0009b2dq1\nay0+4YQTynId34IOvxf/u4ldaJ2Ft3KlhgfMmxpq2gAQEZI2AEQkmu6RMIAj5Q40/uY3v2n0uP33\n31+SdMkll1jZ2LFjLQ6Din5e58SJEy1etWqVpNy532H9Yj+fu6l0lVx66aUWh4EVPxDmn31afh3n\n66+/3uLQPeI99NBDkqQvfelLBV8nJuH9gA/HpXj55Zct9oObYQDNz9OOXevWrSVJY8aMsTKfI/y7\nAYH/20vqhvJzvkNe8c80ybe//W2Lp06dmu+2y4KaNgBEJPNT/jZv3iwp920nXxPO59hjj5WUOwBW\niM997nOSpCVLljT6vUMPPdTicK0Ca9yZmPJXaWFKnyR96lOfqvf5hg0bJKWf5pZSZqb8VUJYuOvw\nww+3srA0qdR8FoyqBL9scNKUP/+7f+yxxypxC0z5A4CYkbQBICKZHIj0XTZXX321pPxdIn4DzrBT\nhyRdeeWVJd1L2CzVr9WbtBjN448/bnFoevodXaqw4FSTUIVdbJqcsKCU7xLx9txzT0nSUUcdVbV7\nQuWQSQAgIiRtAIhIJrtH/GI2frGhxlx77bUW+9emSzVo0CBJ0iuvvGJlF198caPHhK6SKszMaXIK\n2bAWH7j77rsb/Xzu3LmSyrqZb5Pkn6PfYq8xd9xxR6Vup0HUtAEgIpmsaQ8fPjz1d3/yk59Iks46\n66xK3Y4k6TOf+UzBx/gWQxU3As0c3+J45pln6n3euXNni3feeedq3FKU/IJPftAxLJrmhTcGpeL+\ndpu6sIPS7bffbmVnn312ve81tMhWmJMdltKtJmraABARkjYARCST3SO+CZ1vYCq8Zl7pRZo6duxo\ncRjQeeONNxo95rzzzrM4qVsg68KuHLfddpuVrVixwuJPfvKTkqS+ffsmHh8+D0sRSNKAAQPqfW/g\nwIEWV2j3jyZh27ZtFn/605+2OPxG/N9oWHgLdZ544gmLw85AYcE4Kf+7FP369bM47VrylUBNGwAi\nQtIGgIhksnski7p3727xPvvsIym3qyBJ2Pg26/zsjvBKtCR997vflST997//Leq8YQ3zTZs2JX4e\nlh7wa52jYTfddFOjn19xxRUWH3DAAZW+nSj4bd38euL5frtJevfubXFYGqAWqGkDQEQyWdPu1q2b\nxVkZwAuDclLugFBjDj744ErdTln99Kc/tdjX1kqV723WIUOGSKpruTQFfr3wn/3sZ5KSF3LyrZuk\nwfbzzz+/3jn94mf+mG984xuS6p5nc+V/o0cffbSk3J1n/IB4WuE8Ut1zlqT169enOt7X7pNa3rNm\nzbI47YJe1LQBICIkbQCISCa7R8IAmJQ7j7KW/AIy69atS3VMmAuaVcuWLZNUt2a5JH32s5+1OGxu\n6uf/9u/f3+IXX3xRUvL64mm8+uqrRR2XNX5D2WHDhlnslzH4sHzdI4sWLUp9/TAffqedMvlzroit\nW7dKkh5++GErC91RUsNrixfqkUcesfgTn/hEWc7p9ejRw2K/eXhjqGkDQERI2gAQkUzuxu5fwT3u\nuOMa/W5oEv3whz8s9DIN8qPQYfXA2bNnW1na7gDf3EmxjVbVd2MPzXLfPP/d735ncdKqZ9706dMl\nSWeccYaVpW3iSXXN+dNPP93KwrZyrVq1Sn2eFCq6Y7ifv/vAAw/U+9y/mj9q1ChJ0oIFCxo9Jolf\ncS7fK9ejR4+2OMyD93+DYTVFvx60fxfhsMMOS3VPqtFu7MuXL5ckde3atZTT5FXIMy9VA78ddmMH\ngJhlsqbt7+myyy6TlLszjRdqaxMmTLAyPxiU1syZMy32O9OsWbMm1fHTpk2zuE+fPpKk9u3bW1mK\nHVkyUdP2g1n+/pOENx0LqV3nE67Z0KDaIYccIkk66aSTEj9vYFehitYGfQ026b9zly5dLF69erWk\n3Ln+Scf4wd9QUw/HNnRMQ8LvyZ8zrO/u/+5/8YtfWFzAW6o1qWmHNz5XrVpVymnyqnRN27e4vvCF\nLyR9hZo2AMSMpA0AEclk94j34IMPSpJ69uzZ6Pd8c9rPfQwDP36BmLFjx9Y7fu3atRanHWj0TbOP\nfexjFhfZjKp690hYLuDZZ58t5TQ5fFdA6NL6wx/+YGWV2AjVD2Q2cP6ado8knriBedrh7/zGG2+0\nsk6dOknKfUXez+MOc4nnzJljZX4wPVwr373F1D0SnnmlBwfL2T2y3377ScoduPZzy/fYY4+kw+ge\nAYCYkbQBICKZ7x4J9+ebfieffHJpd1QEv6P1fffdJ0nq0KGDlZWhmVb17pHQhPbNtSVLltT7nv93\nhjWwpbptxoYOHWplfl59aML6JqZvtofm+OTJk60s6dV2f83BgwdLyl2NMDQ7pQa7ADLdPXLOOedY\nPH78eEm5/+a0Nm7caPHSpUstDivJ+eecJKbukR//+MeSpGuuuSb1MX55DL96X2PyLTdQiLBze9gi\nMSW6RwAgZpmvadtJ3H362loYsJkxY4aVPf744yVd66CDDrI4vK3Wpk0bK2vZsmVJ529A1WvawZYt\nWyz2tbXA/9t9rbJcmyn7a/p7CcLbe1LRa29XtDZ41VVXWewXMJo/f76k3LnjYceT4cOHW5mvVVdi\n0afQ0tm+fbuV/e1vf5OUO7jp76l169ZpT1+TmnaYLOAnEOQTBnSlgv59tUZNGwBiRtIGgIhE0z2S\nj99K6M9//rPFp5xySr3vJm00265dO4v9oGIV1yiuWfdIM1CTJnwzwbOtLLpHACBmTaam3QRQ064c\naoOVw7OtLGraABAzkjYARISkDQARIWkDQERI2gAQEZI2AESEpA0AESFpA0BESNoAEBGSNgBEpBqv\nsQMAyoSaNgBEhKQNABEhaQNAREjaABARkjYARISkDQARIWkDQERI2gAQEZI2AESEpA0AESFpA0BE\nSNoAEBGSNgBEhKQNABEhaQNAREjaABARkjYARISkDQARIWkDQERI2gAQEZI2AESEpA0AESFpA0BE\nSNoAEBGSNgBEhKQNABEhaQNAREjaABARkjYARISkDQARIWkDQERI2gAQEZI2AESEpA0AESFpA0BE\nSNoAEBGSNgBEhKQNABEhaQNAREjaABARkjYARISkDQARIWkDQERI2gAQEZI2AESEpA0AESFpA0BE\nSNoAEJGdqnCN96twjaagRRHH8GzTKebZSjzfNHi2lVXv+VLTBoCIkLQBICIkbQCICEkbACJC0gaA\niJC0ASAi1ZjyV3aPPPKIxT169CjpXEceeaQkac6cOVbWtm3bks7Z1GzevNniGTNmWHzbbbdJkhYu\nXJh43HvvvSdJ6tevn5XNmjXL4l122aWs9xmTpUuXWjxmzBiL77zzznrfff/9D2bHnXnmmVZ2yy23\nVO7mIhWek1SXI/zv+pprrrG4c+fOkqSRI0cmnmvw4MGScnNBixbFzm4sL2raABARkjYARKSFb1JU\nSEkXWLlypSTpwQcftLILLrjA4i1btpRyemvCf/3rX7eyiRMnWty+ffuSzl+AzL4ROXPmTItPP/30\n1MeFvy3frOzfv7/F06ZNkyTtuuuupd5iPjV5a+/tt9+WJI0fP97K5s6dK0lavHixlb322mupzuf/\nFg877DCLv/a1r0mSTjvtNCtr165dEXdclJq+Eenz17333mvxgAEDynF6TZkyxWL/t7/DDh/Ud6vQ\nZcIbkQAQs0zWtF955RWLv/KVr0iSlixZUr47ckJNO/yfU5L+8Y9/WOxrNBWWuZp2qGH7AbA33ngj\n9fFJNW0v1BBvvfVWK2vVqlWBd5lK1WqDy5cvt/jXv/61JOnGG28s8vKF8TXBqVOnVuWaqnFN+3//\n+5/Fe+yxRzlOmaNbt24WH3rooRZfe+21kqR99tmn7Nf8EGraABAzkjYARCST87T9vOBKdYs05sQT\nT7R4zZo1kqoyWJY5kydPllRYl0gh7rrrLknSXnvtZWXjxo2ryLWqJTwzqXrdIsGKFSss9oObbdq0\nqep9VNMNN9xQ0fM//fTTifHuu+8uSZowYYKV+S7WSqKmDQARIWkDQEQy2T3i55j27t1bkjRv3ryq\nXd83Laswuyazwmvq+bpHwnxkSRo9erTFTz31lCTpmWeeafT49evXF3uLmVPFWRv1+FlP/vcSZuk0\nJeFvzi9pkY9fNmG//faTJG3fvt3K1q5dm/pckyZNkiT17dvXygYOHJj6+FJQ0waAiGRynrYXantn\nn312WW7mw5LmaXuvv/66JGm33XaryPWdzM3TLtWrr74qSerVq5eVPfnkk/W+17FjR4vDm4L77rtv\nOW+lanOJ//Of/1h83HHHSZJefPHFIi//gUGDBknKXXjLD8D5tysD//ZkhVsyNZ2n7QcCf/nLX1oc\nauJHHHGElV199dUWhznXW7dutbKQayTpkksukZRbE08yatQoi3/+858XdO8pMU8bAGJG0gaAiETf\nPRK6LX7/+99b2QknnFDve37+rF9XNwyy0T1SOZs2bbLYvxac1GxfsGCBJOmYY44p5y3UpAkf/n1X\nXXWVlX3nO9+RJP32t7+1sqFDh1oclk34y1/+YmVh8aOddqqbN+AHf8Og2IgRI6zMD7rNnz9fkvT5\nz3++yH9Jo2raPeKfg1887s0335SU+w5AIcKc7AsvvNDKkiZDdOnSxeKwqF2HDh2KumYD6B4BgJhF\nU9P2O3rsvffeFn/5y1+WVNh0G79sY6jRUNOujv3339/i1atX1/u8U6dOkqRFixZZWZieVYKa1gYr\nIQygS9J1110nSbr00kutrLkMRFbas88+a7Ef1EwaoAzTXX3LqgyoaQNAzEjaABCRTL4R6YW5rqEb\nRKrqbjK6/fbbJdUNIKGyNmzYIEnauHGjlZWhe6TJ2bZtm8W+WyRozm/yllPXrl0tDmv7S8kbMFcL\nNW0AiAhJGwAikvnukY9//OM1vf7xxx9f0+s3N2GOvd/aCfVdfvnljX5ehQ1nm53zzjvPYrpHAACp\nZL6mXS5hLquU+7aZn++apMwLF0Vp4cKFFr/wwgv1Pvc7/eQbNPTPO2mw7Prrry/mFpu08HafX5xo\n1qxZjR5T5rnCzZZfUOpXv/pVo98966yzKn07kqhpA0BUSNoAEJEm2T0yd+5ci8MrvI8++qiVJb2y\nXq1NOWNy8803S5KGDx9uZe+++26977Vt29bif/7znxaHrpLZs2dbmd8dhMGyhvndgoYMGSIpf5eI\nd/jhh5f9npqT0C0Snr0k3XPPPfW+d+aZZ1pcra5UMhUARISkDQARaTLdI3/84x8tHjZsmMVhlb5C\n+K2cWrZsWdqNRcBvjjp48GCL161bJym5S8QL24pJdetB+3P5/zZJM0b8NlAHHHBA2tuO2ltvvWXx\n888/b3FYW9s/06TVEJP4mT1NfdZT2Nbtueees7KkGTN+rfJDDjkk9fnD1mFJXSJe2HhcknbcccfU\n5y8FNW0AiEjm19NO4hcT+utf/ypJOvfcc63MD+Lk06dPH0nSXXfdZWW+dl2t/3uqButphx1lfG3h\nX//6VymnzBH+tvINOC5fvtzisJ72ihUrrMzvdlOkzOyuElotY8aMsbK0c6r93+Vll11mcah1jhs3\nzsp23nnnou61CFV7tr7FFwYA77jjjiIvX5ozzjhDUu4ORBXKFaynDQAxI2kDQEQyPxAZmn7Lli2z\nstNOO83izZs3l3T+cK5CthP797//LUlauXJl6mP8q95Z0bdvX0nl7RIpxqmnnlqvLGysKkkTJ060\nuH///pKkj370o5W/sRKFTYr9Egp+znqh/ECa3zw2/DYWL15sZUcffXTR18mqsPWgVLtukSB00fol\nHoJjjz3WYr8Zc7lQ0waAiJC0ASAimZw94meHfPOb35Qk3X///eW7I+fiiy+WVNgr1ffdd58k6ckn\nn0x9jJ9B0ICqzB7xzbnQZZNvHrafe33llVdaPHbs2HrnzLm5lLNHChHmfA8YMKCQw2oyeyR0P82b\nN6+U06TWunVri/2u92n52Sd77723JOnAAw/Md1jFn234+9xzzz2tzM9jz5rOnTtb7GeUhFlCgwYN\nKuR0zB4BgJhlciDyqKOOsnjVqlUVvVYYJGouC0b5t+by1bCDl156yeLvf//7ieXVMmnSJEkF17Sb\nBb/2cyGtwKBXr14WX3TRRZJyB1FrJQz4Z7l27TU0QSG8F1Gq5pGpAKCJIGkDQEQy2T3SvXt3iyvd\nPYL8NmzY0Ojnvmtp5MiRFoetxW699VYr84PMxfCLK6FywrZvWegeKWShp+aAmjYARISkDQARyWT3\nyJQpUywOsxWmTp1aq9tpUvzr+mH+dCFz9Vu1amXx9773PUlSz549reykk06qd8yFF15osZ8LHFZI\ne+edd6xsy5Yt9Y733S8TJkxIfa+1Nn78eEnSiBEjUh8TXj/3M0HCMgNHHHGElfn/ZtOnT5dU2LIK\n+XTp0qVs5ypVWA7AzypLy28HdvDBB1sclq/w2+P5rrewVMXkyZMLvmalUdMGgIhk8o1ILywIFTbo\nlaSXX37Z4rD57OWXX25l3/rWtyz2tfZihGsdf/zxBR/rd73o169fvq9XfT3tadOmSZIeeughKwvz\noD3/hteiRYss7tixYymXN37+7d13313v8zZt2lg8cODAYi5R0/W0K23NmjWSchdV8y2WpAW5glNO\nOcXixx57zOKbbrpJkvTFL34x3+Ur/mxDjgoD21Jd68PPR/e7LgW+lVbIm7lJ1/RmzpwpKXeQPuSg\nhhaxC63M0EJNiTciASBmJG0AiEjmu0eakap3jzQjTbp7pMZ4tpVF9wgAxIykDQARIWkDQERI2gAQ\nEZI2AESEpA0AESFpA0BESNoAEBGSNgBEhKQNABEhaQNAREjaABCRaiwYBQAoE2raABARkjYARISk\nDQARIWkDQERI2gAQEZI2AESEpA0AESFpA0BESNoAEBGSNgBEhKQNABEhaQNAREjaABARkjYARISk\nDQARIWkDQERI2gAQEZI2AESEpA0AESFpA0BE/g+32Zmd236FlQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ab5979ed0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "P = []\n",
    "for ix in range(N_vl):\n",
    "    x = reshape(x_valid[ix],(M,-1))\n",
    "    t = t_valid[ix]\n",
    "    lpt, lpix = logreg(x,t,W,b)    \n",
    "    P.append(lpt[0])\n",
    "\n",
    "P = np.array(P)\n",
    "ix_hard = P.argsort()[:8]\n",
    "ix_easy = P.argsort()[-8:]\n",
    "\n",
    "print ('Hardest 8 Examples')\n",
    "plot_digits(x_valid[ix_hard], numcols=4)\n",
    "\n",
    "\n",
    "print ('Easiest 8 Examples')\n",
    "plot_digits(x_valid[ix_easy], numcols=4)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2. Multilayer perceptron\n",
    "\n",
    "\n",
    "You discover that the predictions by the logistic regression classifier are not good enough for your application: the model is too simple. You want to increase the accuracy of your predictions by using a better model. For this purpose, you're going to use a multilayer perceptron (MLP), a simple kind of neural network. The perceptron wil have a single hidden layer $\\bh$ with $L$ elements. The parameters of the model are $\\bV$ (connections between input $\\bx$ and hidden layer $\\bh$), $\\ba$ (the biases/intercepts of $\\bh$), $\\bW$ (connections between $\\bh$ and $\\log q$) and $\\bb$ (the biases/intercepts of $\\log q$.\n",
    "\n",
    "The conditional probability of the class label $j$ is given by:\n",
    "\n",
    "$\\log p(t = j \\;|\\; \\bx, \\bb, \\bW) = \\log q_j - \\log Z$\n",
    "\n",
    "where $q_j$ are again the unnormalized probabilities per class, and $Z = \\sum_j q_j$ is again the probability normalizing factor. Each $q_j$ is computed using:\n",
    "\n",
    "$\\log q_j = \\bw_j^T \\bh + b_j$\n",
    "\n",
    "where $\\bh$ is a $L \\times 1$ vector with the hidden layer activations (of a hidden layer with size $L$), and $\\bw_j$ is the $j$-th column of $\\bW$ (a $L \\times 10$ matrix). Each element of the hidden layer is computed from the input vector $\\bx$ using:\n",
    "\n",
    "$h_j = \\sigma(\\bv_j^T \\bx + a_j)$\n",
    "\n",
    "where $\\bv_j$ is the $j$-th column of $\\bV$ (a $784 \\times L$ matrix), $a_j$ is the $j$-th element of $\\ba$, and $\\sigma(.)$ is the so-called sigmoid activation function, defined by:\n",
    "\n",
    "$\\sigma(x) = \\frac{1}{1 + \\exp(-x)}$\n",
    "\n",
    "Note that this model is almost equal to the multiclass logistic regression model, but with an extra 'hidden layer' $\\bh$. The activations of this hidden layer can be viewed as features computed from the input, where the feature transformation ($\\bV$ and $\\ba$) is learned.\n",
    "\n",
    "## 2.1 Derive gradient equations (20 points)\n",
    "\n",
    "State (shortly) why $\\nabla_{\\bb} \\mathcal{L}^{(n)}$ is equal to the earlier (multiclass logistic regression) case, and why $\\nabla_{\\bw_j} \\mathcal{L}^{(n)}$ is almost equal to the earlier case.\n",
    "\n",
    "Like in multiclass logistic regression, you should use intermediate variables $\\mathbf{\\delta}_j^q$. In addition, you should use intermediate variables $\\mathbf{\\delta}_j^h = \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial h_j}$.\n",
    "\n",
    "Given an input image, roughly the following intermediate variables should be computed:\n",
    "\n",
    "$\n",
    "\\log \\bq \\rightarrow Z \\rightarrow \\log \\bp \\rightarrow \\mathbf{\\delta}^q \\rightarrow \\mathbf{\\delta}^h\n",
    "$\n",
    "\n",
    "where $\\mathbf{\\delta}_j^h = \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\bh_j}$.\n",
    "\n",
    "Give the equations for computing $\\mathbf{\\delta}^h$, and for computing the derivatives of $\\mathcal{L}^{(n)}$ w.r.t. $\\bW$, $\\bb$, $\\bV$ and $\\ba$. \n",
    "\n",
    "You can use the convenient fact that $\\frac{\\partial}{\\partial x} \\sigma(x) = \\sigma(x) (1 - \\sigma(x))$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer:\n",
    "\n",
    "We are given the following matrices (vectors) with their corresponding dimensions: $\\bW_{L\\times K}$, $\\bV_{M \\times L}$, $\\bx_{M\\times 1}$, $\\bh_{L\\times 1}$,  $\\bb_{K\\times 1}$, $\\ba_{L\\times 1}$. The notation $\\bA_{i\\cdot}$ and $\\bA_{\\cdot j }$  represent the $i$-th row and $j$-th column of the matrix $\\bA$, respectively. Let $\\bj$ denote a $L$ dimensional column vector of ones and $\\bA \\odot \\mathbf{B}$ the Hadamard product between $\\bA$ and $\\mathbf{B}$ of dimensions $m\\times n$.\n",
    "\n",
    "In this case, the model is summarized by the following equations:\n",
    "  * $\\bh = \\sigma( \\bV^T \\bx + \\ba)$\n",
    "  * $\\log \\bq = \\bW^T \\bh + \\bb$\n",
    "  * $\\log q_k = \\bW_{\\cdot k}^T \\bh + \\bb_k$\n",
    "  * $ Z = \\sum_k q_k = \\sum_k \\log(\\exp(q_k))$\n",
    "\n",
    "---\n",
    "\n",
    "Note that $\\nabla_{\\bb} \\mathcal{L}^{(n)} = \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log q_j} \\frac{\\partial \\log q_j}{\\partial b_j} = \\mathbf{\\delta}_j^q \\cdot 1 = \\mathbf{\\delta}_j^q$. Therefore, it is the same as in the previous case.\n",
    "\n",
    "$$\\nabla_{\\bb} \\mathcal{L}^{(n)} = \\mathbf{\\delta}^q = \\mathbb{I}_{\\cdot t^{(n)}} - \\frac{1}{Z} \\bq$$\n",
    "\n",
    "Besides, $\\nabla_{\\bW_{\\cdot j}} \\mathcal{L}^{(n)} = \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\log q_j} \\frac{\\partial \\log q_j}{\\partial \\bW_{\\cdot j}} = \\mathbf{\\delta}_j^q \\frac{\\partial \\bW_{\\cdot j}^T \\bh + b_j}{\\partial \\bW_{\\cdot j}} = \\mathbf{\\delta}_j^q \\bh$. Thus, this is almost the same as the previous case, replacing $\\bx$ by the new input of the logistic regression $\\bh$.\n",
    "\n",
    "$$\\nabla_{\\bW} \\mathcal{L}^{(n)} = \\bh{\\mathbf{\\delta}^q}^T$$\n",
    "\n",
    "---\n",
    "\n",
    "$ \\mathbf{\\delta}^h = \\nabla_{\\bh} \\mathcal{L}^{(n)} = \\frac{\\partial \\log q_{t^{(n)}}}{\\partial \\bh} - \\frac{\\partial \\log Z}{\\partial \\bh} = \\frac{\\partial \\bW_{\\cdot{t^{(n)}}}^T \\bh + \\bb_{t^{(n)}} }{\\partial \\bh} - \\frac{\\partial \\log \\sum_k exp(\\bW_{\\cdot k }^T \\bh + \\bb_k)}{\\partial \\bh} = \\bW_{\\cdot{t^{(n)}}} - \\frac{1}{\\sum_k exp(\\bW_{\\cdot k }^T \\bh + \\bb_k)} \\frac{\\partial \\sum_k exp(\\bW_{\\cdot k}^T \\bh + \\bb_k)}{\\partial \\bh}$\n",
    "$ \\hspace{5mm} = \\bW_{\\cdot{t^{(n)}}} - \\frac{1}{\\sum_k exp(\\bW_{\\cdot k }^T \\bh + \\bb_k)} \\sum_k exp(\\bW_{\\cdot k }^T \\bh + \\bb_k) \\bW_{\\cdot k} = \\bW_{\\cdot{t^{(n)}}} - \\sum_k \\frac{q_k}{Z} \\bW_{\\cdot k} = \\sum_k (\\mathbb{I}_{k,t^{(n)}} - \\frac{q_k}{Z}) \\bW_{\\cdot k} = \\bW (\\mathbb{I}_{\\cdot t^{(n)}} - \\frac{1}{Z} \\bq)$\n",
    "\n",
    "$$ \\mathbf{\\delta}^h = \\bW \\mathbf{\\delta}^q$$\n",
    "\n",
    "---\n",
    "\n",
    "Note that $\\frac { \\partial \\bh_l}{ \\partial \\bV_{ ij } } = \\bh_l (1-\\bh_l) \\frac{\\partial \\bV_{\\cdot l}^T \\bx + \\ba_l}{\\partial \\bV_{ij}} = \\bh_l (1-\\bh_l) \\frac{\\partial \\sum_k \\bV_{kl} \\bx_k + \\ba_l}{\\partial \\bV_{ij}} = \\begin{cases} \\bh_l (1-\\bh_l) \\bx_k \\hspace{5mm} \\text{if} \\hspace{2mm} j=l \\hspace{2mm} \\text{and} \\hspace{2mm} k=i\\\\ 0 \\hspace{5mm} \\text{else}\\end{cases} $\n",
    "\n",
    "$\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\bV_{ij}} =  \\left( \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\bh} \\right)^T \\frac{\\partial \\bh}{\\partial \\bV_{ij}} = {\\mathbf{\\delta}^h}^T \\begin{bmatrix} \\frac { \\partial \\bh_{ 1 } }{ \\partial \\bV_{ ij } }  \\\\ \\vdots  \\\\ \\frac { \\partial \\bh_L}{ \\partial \\bV_{ ij } }  \\end{bmatrix}= \\sum_l \\mathbf{\\delta}_l^h \\frac { \\partial \\bh_l}{ \\partial \\bV_{ ij } } = \\mathbf{\\delta}_j^h \\bh_j (1-\\bh_j) \\bx_i$\n",
    "\n",
    "\n",
    "$$\\nabla_{\\bV} \\mathcal{L}^{(n)} = \\bx (\\mathbf{\\delta}^h \\odot \\bh \\odot (\\bj - \\bh))^T$$\n",
    "\n",
    "---\n",
    "\n",
    "Note that $\\frac { \\partial \\bh_l}{ \\partial \\ba_{ j } } = \\bh_l (1-\\bh_l) \\frac{\\partial \\bV_{\\cdot l}^T \\bx + \\ba_l}{\\partial \\ba_j} = \\begin{cases} \\bh_l (1-\\bh_l) \\hspace{5mm} \\text{if} \\hspace{2mm} j=l\\\\ 0 \\hspace{5mm} \\text{if} \\hspace{2mm} j \\ne l\\end{cases}$\n",
    "\n",
    "$\\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\ba_j} =  \\left( \\frac{\\partial \\mathcal{L}^{(n)}}{\\partial \\bh} \\right)^T \\frac{\\partial \\bh}{\\partial \\ba_j} = {\\mathbf{\\delta}^h}^T \\begin{bmatrix} \\frac { \\partial \\bh_{ 1 } }{ \\partial \\ba_{ j } }  \\\\ \\vdots  \\\\ \\frac { \\partial \\bh_L}{ \\partial \\ba_{ j } }  \\end{bmatrix}= \\sum_l \\mathbf{\\delta}_l^h \\frac { \\partial \\bh_l}{ \\partial \\ba_{ j } } = \\mathbf{\\delta}_j^h \\bh_j (1-\\bh_j)$\n",
    "\n",
    "$$\\nabla_{\\ba} \\mathcal{L}^{(n)} = \\mathbf{\\delta}^h \\odot \\bh \\odot (\\bj - \\bh) $$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 MAP optimization (10 points)\n",
    "\n",
    "You derived equations for finding the _maximum likelihood_ solution of the parameters. Explain, in a few sentences, how you could extend this approach so that it optimizes towards a _maximum a posteriori_ (MAP) solution of the parameters, with a Gaussian prior on the parameters. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer:\n",
    "\n",
    "Recall that $posterior = \\frac{likelihood \\times prior}{evidence} \\rightarrow \\log posterior = \\log likelihood + \\log prior - \\log evidence.$\n",
    "\n",
    "Therefore, $\\theta_{MAP} = \\text{argmax}_{\\theta} \\log posterior = \\text{argmax}_{\\theta} \\log likelihood + \\log prior$. In this case it would be equivalent to finding the parameters that maximize a function $\\mathcal{L} - \\lambda \\mathcal{R}$, where $\\mathcal{L}$ corresponds to the log-likelihood (as before) and $-\\lambda \\mathcal{R}$ corresponds to the regularization term for a positive hyperparameter $\\lambda$ associated with the parameters of the prior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3. Implement and train a MLP (15 points)\n",
    "\n",
    "Implement a MLP model with a single hidden layer, and code to train the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.special import expit\n",
    "\n",
    "def mlp_gradient(x,t,W,b,V,a):\n",
    "    \n",
    "    h = expit(np.dot(V.T,x) + a)\n",
    "    lq = np.dot(W.T,h) + b\n",
    "    Z = sum(np.exp(lq))\n",
    "    \n",
    "    delta_q = - np.exp(lq)/Z\n",
    "    delta_q[t] += 1\n",
    "    \n",
    "    delta_h = dot(W,delta_q)\n",
    "    \n",
    "    grad_b = delta_q\n",
    "    grad_W = dot(h,delta_q.T)\n",
    "    \n",
    "    grad_a = delta_h * h * (ones(shape(h)) - h)\n",
    "    grad_V = dot(x,grad_a.T)\n",
    "    \n",
    "    return grad_W, grad_b, grad_V, grad_a\n",
    "\n",
    "def mlp_sgd_iter(x_train, t_train, W, b, V, a, eta):\n",
    "    N, M = shape(x_train)\n",
    "    perm = np.random.permutation(N)\n",
    "    for ix in perm:\n",
    "        x = reshape(x_train[ix],(M,-1))\n",
    "        t = t_train[ix]\n",
    "        grad_W, grad_b, grad_V, grad_a = mlp_gradient(x, t, W, b, V, a)\n",
    "        W += eta * grad_W\n",
    "        b += eta * grad_b\n",
    "        V += eta * grad_V\n",
    "        a += eta * grad_a\n",
    "    return W, b, V, a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter: 0 Error: 92.13%\n",
      "Iter: 1 Error: 4.162%\n",
      "Iter: 2 Error: 2.534%\n",
      "Iter: 3 Error: 2.072%\n",
      "Iter: 4 Error: 1.656%\n",
      "Iter: 5 Error: 1.472%\n",
      "Iter: 6 Error: 0.646%\n",
      "Iter: 7 Error: 0.5%\n",
      "Iter: 8 Error: 0.226%\n",
      "Iter: 9 Error: 0.138%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEPCAYAAABGP2P1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHFJJREFUeJzt3Xt4VfW95/H3NwnhIiFyUW7halHUSklbb4VOY9tRZHzq\nmU6ngNVO7ZzK9NTK2M5R2+fhGPucjqdzjs+xjr3ojKO1jy1OpVrPSA/SqWlHZ1RquYmIdwhyMwIR\nRCEm3/ljrQ17b3ayL9lrr335vJ5nP9lr5cfa3wTIJ7/1Xeu3zd0RERFJqIu7ABERKS8KBhERSaFg\nEBGRFAoGERFJoWAQEZEUCgYREUkRaTCY2T1mtsfMNg4w5g4ze9nM1pvZ3CjrERGR7KKeMdwLXNLf\nJ83sUuA0d58FLAV+GnE9IiKSRaTB4O5PAvsHGHI5cH849hmg2czGR1mTiIgMLO4ew2SgM2n7zXCf\niIjEJO5gEBGRMtMQ8+u/CUxJ2m4J953AzLSok4hIAdzd8hlfihmDhY9MHgW+DGBmFwAH3H1Pfwdy\n97J63HzzzbHXUCl1qSbVVAt1lWNNhYh0xmBmvwDagLFmth24GWgE3N3vdvdVZrbQzF4B3gWujrIe\nERHJLtJgcPcrchhzbZQ1iIhIftR8HoS2tra4S8ioHOtSTblRTbkrx7rKsaZCWKHnoErNzLxSahUR\nKRdmhufZfI77qiQRqULTp09n27ZtcZdRU6ZNm8Ybb7xRlGNpxiAiRRf+lhp3GTWlv+95ITMG9RhE\nRCSFgkFERFIoGEREJIWCQURkEPr6+mhqamLHjh1FHRsnNZ9FpOjKufnc1NSEWdCLfffddxk6dCj1\n9fWYGXfddRdLliyJucLCFLP5rGAQkaIr52BINnPmTO655x4uuuiifsf09vZSX19fwqoKo6uSRESK\nINNCc8uXL2fx4sVcccUVNDc388ADD/D0009z4YUXMnr0aCZPnsyyZcvo7e0FguCoq6tj+/btAFx1\n1VUsW7aMhQsXMmrUKObNm3fsno58xgL89re/5YwzzmD06NFcd911zJ8/n/vvvz/y74uCQUQkzSOP\nPMKVV15Jd3c3ixYtYsiQIdxxxx3s27ePp556itWrV3PXXXcdG584NZXwy1/+ku9///vs37+fKVOm\nsHz58rzH7t27l0WLFnHbbbfR1dXFjBkzWLt2bYRf9XEKBhGJhVlxHlGYP38+CxcuBGDo0KF87GMf\n49xzz8XMmD59Ol/72tf4wx/+cGx8+qzjC1/4Aq2trdTX1/OlL32J9evX5z32scceo7W1lcsuu4z6\n+nquv/56xo4dG80XnEZLYohILMq5BTFlypSU7a1bt/Ltb3+b5557jsOHD9Pb28v555/f75+fMGHC\nsecjRozg0KFDeY/duXPnCXW0tLTk9XUUqqJmDOEpPRGRSKWf7lm6dCnnnHMOr732Gt3d3dxyyy2R\nN9cnTpxIZ2dnyr4338z4BpdFV1HB8OqrcVcgIrXo4MGDNDc3M3z4cLZs2ZLSX4jKZZddxrp163js\nscfo7e3l9ttvp6urK/LXhQoLhnXr4q5ARKpJ+sygP7fddhv33Xcfo0aN4utf/zqLFy/u9zjZjpnr\n2FNPPZUHH3yQ66+/nnHjxvH666/T2trK0KFDc6p5MCrqPoabbnJuvTXuSkQkm0q5j6GS9PX1MWnS\nJFauXMm8efNO+HzN3segGYOI1JLVq1fT3d3NkSNH+N73vkdjYyPnnXde5K9bUcGQdMWXiEjVe/LJ\nJ5k5cybjx49nzZo1PPLIIwwZMiTy162oU0ljxzqbNsHEiXFXIyID0amk0qvZU0lz52rWICIStYoK\nhtZW9RlERKJWUcGgGYOISPQqakmM1lZob4+7ChHJZtq0aTnfIyDFMW3atKIdq6Kazz09TnMz7N4N\nTU1xVyQiUv6qvvnc0AAf/jBs2BB3JSIi1auiggHUZxARiVrFBYOuTBIRiVbFBYNmDCIi0aqo5rO7\nc/gwjBsH3d1QgjvDRUQqWtU3nwFGjIBp02DLlrgrERGpThUXDKA+g4hIlCoyGNRnEBGJTkUGg2YM\nIiLRqbjmM8Bbb8GsWbB/P+iuexGR/pVl89nMFpjZi2b2kpndmOHzo8zsUTNbb2abzOwr2Y55yikw\nciS88UYUFYuI1LZIg8HM6oA7gUuAs4ElZjY7bdg3gM3uPhe4CLjNzLIu7qc+g4hINKKeMZwHvOzu\n29y9B1gBXJ42xoHEknhNwNvu/kG2A6vPICISjaiDYTLQmbS9I9yX7E7gLDPbCWwAluVyYM0YRESi\nUQ7vx3AJsM7dP21mpwFrzGyOux9KH9ie9GYMp5/exrp1bSUrUkSkEnR0dNDR0TGoY0R6VZKZXQC0\nu/uCcPsmwN39B0lj/hdwq7s/FW7/b+BGd/9T2rE8uda+Phg9Gl59NVgiQ0RETlSOVyWtBT5kZtPM\nrBFYDDyaNmYb8FkAMxsPnA68lu3AdXU6nSQiEoVIg8Hde4FrgceBzcAKd99iZkvN7Jpw2N8CnzCz\njcAa4AZ335fL8RUMIiLFV5E3uCXcdx+sWQMPPBBPTSIi5a4cTyVFqrVVMwYRkWKr6BnD0aNw8snQ\n1RUsxy0iIqlqbsbQ2AizZ8Pzz8ddiYhI9ajoYICgAa07oEVEiqfig0F9BhGR4qr4YNCMQUSkuCq6\n+QzwzjswaRJ0d0N9fQyFiYiUsZprPgOMGgUTJsBLL8VdiYhIdaj4YAD1GUREiqkqgkF9BhGR4qmK\nYNCMQUSkeKoiGBIzhgrpo4uIlLWqCIaJE4NluHfujLsSEZHKVxXBYKY+g4hIsVRFMID6DCIixVI1\nwaAZg4hIcVRNMGjGICJSHBW/JEZCb2/w3gw7dkBzcwkLExEpYzW5JEZCfT2ccw5s2BB3JSIila1q\nggHUZxARKYaqCgb1GUREBq+qgkEzBhGRwaua5jPAe+/B2LFw4EDwftAiIrWuppvPAMOHw8yZsHlz\n3JWIiFSuqgoGUJ9BRGSwqi4Y1GcQERmcqgsGzRhERAanqprPAPv2wYwZsH9/sBS3iEgtq/nmM8CY\nMcHSGK+/HnclIiKVqeqCAdRnEBEZjKoMBvUZREQKV5XBoBmDiEjhqjIYNGMQESlcVQbD1KnB8hh7\n98ZdiYhI5anKYDALTidp1iAikr/Ig8HMFpjZi2b2kpnd2M+YNjNbZ2bPm9kTxXhd9RlERArTEOXB\nzawOuBP4DLATWGtmv3H3F5PGNAM/Ai529zfNbFwxXru1FVatKsaRRERqS9QzhvOAl919m7v3ACuA\ny9PGXAGsdPc3Ady9qxgvrBmDiEhhog6GyUBn0vaOcF+y04ExZvaEma01s6uK8cKzZ8P27fDuu8U4\nmohI7SiH5nMD8FHgUmABsNzMPjTYgw4ZAmedBRs3DvZIIiK1JdIeA/AmMDVpuyXcl2wH0OXu7wPv\nm9kfgY8Ar6QfrL29/djztrY22traBnzxxP0MF15YSOkiIpWno6ODjo6OQR0j0tVVzawe2ErQfN4F\nPAsscfctSWNmA/+VYLYwFHgGWOTuL6QdK6fVVZP96EewYQPcffegvgwRkYpVdqurunsvcC3wOLAZ\nWOHuW8xsqZldE455EVgNbASeBu5OD4VC6Q5oEZH8Vd37MSQ7dAjGj4fubmiI+qSZiEgZKrsZQ9xG\njoSWFti6Ne5KREQqR1UHA+h+BhGRfFV9MKjPICKSn6oPBs0YRETyU9XNZ4A9e4Ib3bq6glVXRURq\niZrPGYwfD42N0NmZfayIiNRAMID6DCIi+aiJYFCfQUQkd1mDwczqzewfSlFMVDRjEBHJXdZgCJe1\nmF+CWiKjGYOISO5yuirJzH5C8D4KvwKOvcOBu/86utJOqKGgq5IA+vrg5JNh2zYYPbrIhYmIlLFC\nrkrKdQWhYcDbwKeT9jlQsmAYjLo6mDMnOJ100UVxVyMiUt5yCgZ3vzrqQqKW6DMoGEREBpbTVUlm\n1mJmD5vZ3vCx0sxaoi6umNRnEBHJTa6Xq94LPApMCh//FO6rGLoySUQkN7k2n9e7+9xs+6I0mOYz\nwJEjQeN53z4YNqyIhYmIlLEol8R428yuDO9pqDezKwma0RVj6FCYNQs2b467EhGR8pZrMHwV+CKw\nm+C9m78AVFxDWn0GEZHssl6VZGb1wOfd/XMlqCdS6jOIiGSX653PS0pQS+Q0YxARyS7X5vM/AkOA\nB0m98/nP0ZV2Qg2Daj4DHDgAU6ZAd3dw05uISLUrpPmcazA8kWG3u/unM+yPRDGCAWDGDFi9Gk4/\nvQhFiYiUuUiWxDCzOuAn7v4/C66sjCT6DAoGEZHMcukx9AE3lKCWklCfQURkYLmeaf+dmf0nM5ti\nZmMSj0gri4iuTBIRGViuPYbXM+x2d59Z/JL6raEoPYbOTjj3XNi9uwhFiYiUuciW3Xb3GYWVVH5a\nWuCDD4JgmDAh7mpERMrPgKeSzOyGpOf/Nu1z/zmqoqJkpj6DiMhAsvUYFic9/07a5xYUuZaSUZ9B\nRKR/2YLB+nmeabtiaMYgItK/bMHg/TzPtF0xNGMQEenfgFclmVkvwRIYBgwHDic+BQxz9yGRV3i8\nlqJclQRB87m5OWhANzUV5ZAiImWp6O/H4O717j7K3ZvcvSF8ntguWSgUW0MDfPjDsHFj3JWIiJSf\nml1KTn0GEZHMajYY1GcQEcmsZoNBMwYRkcwiDwYzW2BmL5rZS2Z24wDjzjWzHjP7fNQ1AcyZA1u2\nQE9PKV5NRKRyRBoM4ZLddwKXAGcDS8xsdj/j/g5YHWU9yUaMgGnTgnAQEZHjop4xnAe87O7b3L0H\nWAFcnmHcN4GHgL0R15NCfQYRkRNFHQyTgc6k7R3hvmPMbBLwF+7+E0p8N7X6DCIiJyqH5vPtQHLv\noWThoBmDiMiJclp2exDeBKYmbbeE+5J9HFhhZgaMAy41sx53fzT9YO3t7ceet7W10dbWNqji5s4N\ngsE9WHVVRKTSdXR00NHRMahj5PRGPQUf3Kwe2Ap8BtgFPAsscfeMLV8zuxf4J3f/dYbPFW1JjGQt\nLfDkkzB9etEPLSISu6IviTFY7t4LXAs8DmwGVrj7FjNbambXZPojUdaTifoMIiKpIp0xFFNUM4bl\ny6GuDm65peiHFhGJXdnNGCqBZgwiIqlqPhh0ZZKISKqaD4bp06G7G95+O+5KRETKQ80HQ13d8ctW\nRUREwQCozyAikkzBgPoMIiLJFAxoxiAikqzm72MAOHoUTj45aEAPHx7JS4iIxEL3MRSosRHOOAM2\nbYq7EhGR+CkYQuoziIgEFAwh9RlERAIKhpBmDCIiATWfQ++8A5MmBXdB19dH9jIiIiWl5vMgjBoF\nEybAyy/HXYmISLwUDElaW9VnEBFRMCTRmkkiIgqGFJoxiIgoGFIkZgwV0o8XEYmEgiHJxIlgBjt3\nxl2JiEh8FAxJzNRnEBFRMKRRn0FEap2CIY1mDCJS6xQMaTRjEJFapyUx0vT2Bu/NsGMHNDdH/nIi\nIpHSkhhFUF8P55wDGzbEXYmISDwUDBmozyAitUzBkIH6DCJSyxQMGWjGICK1TM3nDN57D8aOhQMH\ngveDFhGpVGo+F8nw4TBzJrzwQtyViIiUnoKhH+oziEitUjD0Q30GEalVCoZ+aMYgIrVKzed+7NsH\nM2bA/v1Qp/gUkQql5nMRjRkTLI3x+utxVyIiUloKhgGozyAitSjyYDCzBWb2opm9ZGY3Zvj8FWa2\nIXw8aWbnRF1TrtRnEJFaFGkwmFkdcCdwCXA2sMTMZqcNew34F+7+EeBvgf8WZU350IxBRGpR1DOG\n84CX3X2bu/cAK4DLkwe4+9Pu3h1uPg1MjrimnGnGICK1KOpgmAx0Jm3vYOAf/H8J/DbSivIwdWqw\nPMbevXFXIiJSOg1xF5BgZhcBVwPz+xvT3t5+7HlbWxttbW0R13T8dNLFF0f6UiIiRdHR0UFHR8eg\njhHpfQxmdgHQ7u4Lwu2bAHf3H6SNmwOsBBa4+6v9HKuk9zEkfOtbMGEC3HBDyV9aRGTQyvE+hrXA\nh8xsmpk1AouBR5MHmNlUglC4qr9QiJP6DCJSayINBnfvBa4FHgc2AyvcfYuZLTWza8Jhy4ExwI/N\nbJ2ZPRtlTfnSlUkiUmu0JEYWPT3BHdB798JJJ5X85UVEBqUcTyVVvCFD4MwzYdOmuCsRESkNBUMO\n1GcQkVqiYMjB3LkKBhGpHQqGHLS2qgEtIrVDzeccHDoE48dDdzc0lM0tgSIi2an5HJGRI6GlBbZu\njbsSEZHoKRhypD6DiNQKBUOO1GcQkVqhYMiRZgwiUivUfM7Rnj1w1lnQ1RWsuioiUgnUfI7Q+PHQ\n2AidndnHiohUMgVDHtRnEJFaoGDIg/oMIlILFAx50IxBRGqBgiEPmjGISC1QMOThtNNg3z7Yvz/u\nSkREoqNgyENdHcyZo9NJIlLdFAx5Up9BRKqdgiFP6jOISLVTMORJMwYRqXZaEiNPR47A6NFBE3rY\nsLirEREZmJbEKIGhQ2HWLNi8Oe5KRESioWAogPoMIlLNFAwFUJ9BRKqZgqEAmjGISDVT87kABw7A\nlCnQ3R3c9CYiUq7UfC6Rk0+GcePglVfirkREpPgUDAVSn0FEqpWCoUDqM4hItVIwFEgzBhGpVmo+\nF2j3bpg9G5qaYOLEgR/jx0NDQ9wVi0gtKqT5rGAYhCNHgoDYtWvgR1cXjBmTPUAmTtQyGyJSXAqG\nMtXbC3v3Zg+Q3bthxIjcAqSpCSyvv2oRqUUKhgrnHizOly1Adu0KxmYKjAkTgnAZOjR4DBt2/PlA\nD4WMSHVSMNSQgwf7n3UcPhyc5sr2eP/94GNPDzQ25hYguQZNpnGJ12hsTH3e3z7dPCgyeGUZDGa2\nALid4Aqoe9z9BxnG3AFcCrwLfMXdT7jeR8EQnb4+OHo0txDJJ3CSH4njHz2a+nygfXV1+YdJPvsa\nGqC+/sRHXV3u+4s5VkEoUSgkGCK9VsbM6oA7gc8AO4G1ZvYbd38xacylwGnuPsvMzgd+ClwQZV3F\n0tHRQVtbW9xlnCDfuurqgt/wo2x851uTe9CbKTRUMu17//1gGZPEvs7ODk49tY3eXo49+vpI2S5k\nf6HHAGho6GDUqDaamjj2GDmy8OfFOE1YLf/OS6EcaypE1BdRnge87O7bAMxsBXA58GLSmMuB+wHc\n/Rkzazaz8e6+J+LaBq1c/xGUY1351mQW/Ebf0AAnnRRNTe3tHbS3515T1Pr64G/+poNly9o4dCg4\nXXjwIP0+37Xr+PP+xvX1DS5Ymppg5coOpk5to64udXaT/Ejfl74dRQ+rGv6dl6uog2Ey0Jm0vYMg\nLAYa82a4r+yDQaSY6uqCIDzllOBRDEeP9h8s6c+7ujJ/rrMTHnvs+Mymry/1kb4v0xiz/MMkWwC9\n9RasWgVDhhz/JSLxSN9XqjE7d8Jzzx0Pw8TXnfyxlM8LPT2p265EqlhjY3APzZgxhR+jvT14FMo9\neOQTJtkCqLcXfvxj+OpX4YMPjj96elK3M+3rb8x77+V3nPR9PT2wYwf86U/Hv97kr72/fVE/L0Sk\nzWczuwBod/cF4fZNgCc3oM3sp8AT7v5guP0i8Kn0U0lmps6ziEgByqr5DKwFPmRm04BdwGJgSdqY\nR4FvAA+GQXIgU38h3y9MREQKE2kwuHuvmV0LPM7xy1W3mNnS4NN+t7uvMrOFZvYKweWqV0dZk4iI\nDKxibnATEZHSKPtbaszsHjPbY2Yb464lwcxazOz3ZrbZzDaZ2XVlUNNQM3vGzNaFNd0cd00JZlZn\nZn82s0fjriXBzN4wsw3h9+vZuOsBCC/V/pWZbQn/bZ0fcz2nh9+fP4cfu8vk3/r1Zva8mW00swfM\nrLEMaloW/r+L9edBpp+XZjbazB43s61mttrMmrMdp+yDAbgXuCTuItJ8AHzL3c8GLgS+YWaz4yzI\n3Y8AF7l7KzAXuNTM0i8Njssy4IW4i0jTB7S5e6u7l8v36YfAKnc/E/gIsCXOYtz9pfD781HgYwSn\neh+OsyYzmwR8E/iou88hOB2+OOaazgb+PfBxgv97l5nZzJjKyfTz8ibgd+5+BvB74DvZDlL2weDu\nTwL7464jmbvvTizb4e6HCP4DT463KnD3w+HToQT/YWI/T2hmLcBC4L/HXUsao4z+/ZvZKOCT7n4v\ngLt/4O7vxFxWss8Cr7p7Z9aR0asHTjKzBmAEwaoKcToTeMbdj7h7L/BH4PNxFNLPz8vLgZ+Fz38G\n/EW245TNf4xKZWbTCX5LeCbeSo6dslkH7AbWuPvauGsC/hH4a8ogpNI4sMbM1prZ1+IuBpgBdJnZ\nveGpm7vNbHjcRSVZBPwy7iLcfSdwG7Cd4GbYA+7+u3ir4nngk+EpmxEEvwhNibmmZKcmrvR0993A\nqdn+gIJhEMxsJPAQsCycOcTK3fvCU0ktwPlmdlac9ZjZvwL2hLMrCx/lYl54imQhwanA+THX0wB8\nFPhRWNdhglMAsTOzIcDngF+VQS0nE/wGPA2YBIw0syvirClc++0HwBpgFbAO6I2zpiyy/pKmYChQ\nOI19CPi5u/8m7nqShacgngAWxFzKPOBzZvYawW+bF5nZ/THXBIC77wo/vkVw3jzuPsMOoNPd/xRu\nP0QQFOXgUuC58HsVt88Cr7n7vvC0za+BT8RcE+5+r7t/3N3bgAPASzGXlGyPmY0HMLMJwN5sf6BS\ngqHcftsE+B/AC+7+w7gLATCzcYmrDcJTEP+S1MUKS87dv+vuU919JkGD8Pfu/uU4awIwsxHhbA8z\nOwm4mOB0QGzCqX6nmZ0e7voM5dOwX0IZnEYKbQcuMLNhZmYE36dYm/QAZnZK+HEq8K+BX8RZDqk/\nLx8FvhI+/3dA1l9ky36tJDP7BdAGjDWz7cDNiQZdjDXNA74EbArP6TvwXXf/5xjLmgj8LFzqvA54\n0N1XxVhPORsPPBwus9IAPODuj8dcE8B1wAPhqZvXKIObPcNz5p8From7FgB3f9bMHiI4XdMTfrw7\n3qoAWGlmYwhq+qu4LhzI9PMS+DvgV2b2VWAb8MWsx9ENbiIikqxSTiWJiEiJKBhERCSFgkFERFIo\nGEREJIWCQUREUigYREQkhYJBao6ZHQw/TjOz9HcUHOyxv5O2/WQxjy9SCgoGqUWJm3dmAHmts2Nm\n9VmGfDflhdzjXoNJJG8KBqlltwLzw9VMl4Wr0/6X8A2P1idWXTWzT5nZH83sN8DmcN/D4cqsm8zs\nL8N9twLDw+P9PNx3MPFiZvb34fgNZvbFpGM/kfQGPT8v8fdA5ARlvySGSIRuAr7t7p8DCIPggLuf\nH74r2FNmllgqoxU42923h9tXu/sBMxsGrDWzle7+HTP7Rrg6aoKHx/43wBx3P8fMTg3/zB/CMXOB\nswiWS3/KzD7h7v83yi9cZCCaMYgcdzHw5XD9q2eAMcCs8HPPJoUCwH80s/XA0wTLnM9iYPMIF6Jz\n971AB3Bu0rF3ebA+zXpg+uC/FJHCacYgcpwB33T3NSk7zT5F8LaWydufBs539yNm9gQwLOkYub5W\nwpGk573o/6XETDMGqUWJH8oHgaak/auBvwrfawMzmxWuLpquGdgfhsJs4IKkzx1N/Pm01/o/wKKw\nj3EK8Eng2SJ8LSJFp99MpBYlrkraCPSFp47uc/cfhm/V+udwrf+9ZH5/3H8G/oOZbQa2Av8v6XN3\nAxvN7Dl3vyrxWu7+sJldAGwA+oC/dve9ZnZmP7WJxEbLbouISAqdShIRkRQKBhERSaFgEBGRFAoG\nERFJoWAQEZEUCgYREUmhYBARkRQKBhERSfH/AZeqZALGQhN0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a3409bc10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def mlp_class(x,t,W,b,V,a):\n",
    "    h = expit(np.dot(V.T,x) + a)\n",
    "    lq = np.dot(W.T,h) + b\n",
    "    Z = sum(np.exp(lq))\n",
    "    offset = max(lq)\n",
    "    lZ = offset  + log(sum(np.exp(lq - ones(shape(lq))*offset)))\n",
    "    lp = lq - ones(shape(lq))*lZ\n",
    "    return argmax(lp)\n",
    "\n",
    "\n",
    "def mlp_pred_iter(x_set, t_set, W, b, V, a):\n",
    "    N_set, M = shape(x_set)\n",
    "    misclass = 0\n",
    "    for ix in range(N_set):\n",
    "        x = reshape(x_set[ix],(M,-1))\n",
    "        t = t_set[ix]\n",
    "        misclass += mlp_class(x,t,W,b,V,a) != t\n",
    "    return misclass*1.0/N_set\n",
    "\n",
    "\n",
    "N_tr, M = shape(x_train)\n",
    "N_vl, M = shape(x_test)\n",
    "\n",
    "K = 10\n",
    "L = 200\n",
    "max_iter = 10\n",
    "\n",
    "W_mlp = np.random.normal(0,0.2,(L, K))\n",
    "b_mlp = np.random.normal(0,0.2,(K, 1))\n",
    "V_mlp = 2*np.random.random_sample((M, L)) - 1\n",
    "a_mlp = 2*np.random.random_sample((L, 1)) - 1\n",
    "\n",
    "E_tr = []\n",
    "\n",
    "\n",
    "for it in range(max_iter):\n",
    "    eta = 0.1\n",
    "    err = mlp_pred_iter(x_train, t_train, W_mlp, b_mlp, V_mlp, a_mlp)\n",
    "    print ('Iter: ' + str(it) + ' Error: '+ str(100*err) + '%')\n",
    "    E_tr.append(err)\n",
    "    W_mlp, b_mlp, V_mlp, a_mlp = mlp_sgd_iter(x_train, t_train, W_mlp, b_mlp, V_mlp, a_mlp, eta)\n",
    "\n",
    "plot(range(1,max_iter+1), E_tr, label='Training')\n",
    "xlabel('Iteration')\n",
    "ylabel('Error')\n",
    "legend(loc=1)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3.1. Less than 250 misclassifications on the test set (10 bonus points)\n",
    "\n",
    "You receive an additional 10 bonus points if you manage to train a model with very high accuracy: at most 2.5% misclasified digits on the test set. Note that the test set contains 10000 digits, so you model should misclassify at most 250 digits. This should be achievable with a MLP model with one hidden layer. See results of various models at : `http://yann.lecun.com/exdb/mnist/index.html`. To reach such a low accuracy, you probably need to have a very high $L$ (many hidden units), probably $L > 200$, and apply a strong Gaussian prior on the weights. In this case you are allowed to use the validation set for training.\n",
    "You are allowed to add additional layers, and use convolutional networks, although that is probably not required to reach 2.5% misclassifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test error: 2.35%\n"
     ]
    }
   ],
   "source": [
    "err = mlp_pred_iter(x_test, t_test, W_mlp, b_mlp, V_mlp, a_mlp)\n",
    "print ('Test error: '+ str(100*err)+'%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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
}