lab3_unsolved.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab  3 - Implementing Logistic Regression with Gradient Descent using Theano\n",
    "\n",
    "## Learning Aims\n",
    "\n",
    "- Using Theano\n",
    "- Logistic regression\n",
    "- Supervised learning\n",
    "- Automatic differentiation\n",
    "- Gradient descent\n",
    "$\\newcommand{\\Ex}{\\mathbb{E}}$\n",
    "$\\newcommand{\\Var}{\\mathrm{Var}}$\n",
    "$\\newcommand{\\Cov}{\\mathrm{Cov}}$\n",
    "$\\newcommand{\\SampleAvg}{\\frac{1}{N({S})} \\sum_{s \\in {S}}}$\n",
    "$\\newcommand{\\indic}{\\mathbb{1}}$\n",
    "$\\newcommand{\\avg}{\\overline}$\n",
    "$\\newcommand{\\est}{\\hat}$\n",
    "$\\newcommand{\\trueval}[1]{#1^{*}}$\n",
    "$\\newcommand{\\Gam}[1]{\\mathrm{Gamma}#1}$\n",
    "\n",
    "## Installing Theano\n",
    "\n",
    "We shall be using Theano 0.6 in this course:\n",
    "http://deeplearning.net/software/theano/index.html\n",
    "\n",
    "### Installation\n",
    "\n",
    "In Anaconda:\n",
    "\n",
    "$ pip install Theano\n",
    "\n",
    "If you are not using Anaconda, additional installation instructions can be found here:\n",
    "\n",
    "http://deeplearning.net/software/theano/install.html\n",
    "\n",
    "### Theano at a glance\n",
    "\n",
    "Theano is a Python library that lets you to define, optimize, and evaluate mathematical expressions, especially ones with multi-dimensional arrays. Some of its advantages over other libraries, such as numpy are:\n",
    "\n",
    "- Execution speed optimizations: Theano can use g++ or nvcc to compile parts your expression graph into CPU or GPU instructions, which run much faster than pure Python.\n",
    "- Symbolic differentiation: Theano can automatically build symbolic graphs for computing gradients."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Theano\n",
    "One powerful feature of Theano is that once you define symbolic variables, you can write expressions for those variables and convert these expressions into callable objects that can evaluate the expression for any value of the variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(5.8)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import theano\n",
    "from theano import *\n",
    "import theano.tensor as T\n",
    "# Many of the types and symbols that we use are within the Theano subpackage theno.tensor.\n",
    "\n",
    "# declare two symbolic floating-point scalars\n",
    "a = T.dscalar('a')\n",
    "b = T.dscalar('b')\n",
    "\n",
    "# create a simple expression\n",
    "c = a + b\n",
    "\n",
    "# convert the expression into a callable object that takes (a,b)\n",
    "# values as input and computes a value for c\n",
    "f = function([a,b], c)\n",
    "# Behind the scenes, when you call theano.function, f is compiled into C code.\n",
    "\n",
    "# bind 1.5 to 'a', 2.5 to 'b', and evaluate 'c'\n",
    "f(1.3,4.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same can be do with vectors or matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 11.,  22.],\n",
       "       [ 33.,  44.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# declare two symbolic floating-point matrices\n",
    "x = T.dmatrix('y')\n",
    "y = T.dmatrix('y')\n",
    "\n",
    "# expression\n",
    "z = x + y\n",
    "\n",
    "# convert into callable\n",
    "f = function([x,y], z)\n",
    "\n",
    "# Evaluate\n",
    "f([[1, 2], [3, 4]], [[10, 20], [30, 40]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are other available Theano types, and all of them we can preform similar operations.\n",
    "    \n",
    "- byte: bscalar, bvector, bmatrix, brow, bcol, btensor3, btensor4\n",
    "- 16-bit integers: wscalar, wvector, wmatrix, wrow, wcol, wtensor3, wtensor4\n",
    "- 32-bit integers: iscalar, ivector, imatrix, irow, icol, itensor3, itensor4\n",
    "- 64-bit integers: lscalar, lvector, lmatrix, lrow, lcol, ltensor3, ltensor4\n",
    "- nfloat: fscalar, fvector, fmatrix, frow, fcol, ftensor3, ftensor4\n",
    "- double: dscalar, dvector, dmatrix, drow, dcol, dtensor3, dtensor4\n",
    "- complex: cscalar, cvector, cmatrix, crow, ccol, ctensor3, ctensor4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic function\n",
    "Now suppose we want to evaluate the logistic function on a matrix, on an element-wise basis. The one-dimensional logistic curve is given by:\n",
    "\n",
    "$$s(x) = \\frac{1}{1 + e^{-x}}$$\n",
    "\n",
    "We can build a Theano graph that represents the logistic function, and then evaluate it in a 2x2 matrix as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.5       ,  0.73105858],\n",
       "       [ 0.26894142,  0.11920292]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = T.dmatrix('x')\n",
    "s = 1 / (1 + T.exp(-x))\n",
    "logistic = function([x], s)\n",
    "logistic([[0, 1], [-1, -2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also avaluate it on a vector and plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/python3/lib/python3.6/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['generic', 'test', 'logistic', 'f', 'dot', 'gradient']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n",
      "  \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10be77cc0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9xvHPNzuBJGxhDSHsiywCEbfrVlFB69LFKq1t\nba203tpaq/ZqF9vrrb12sWpd6tLaarVat3uLBRUVrdYqsij7EggYwhoSSCAh28z3/pGBmyKQkExy\nMpPn/XrlNTNnDpnnQObhl3PO/I65OyIiEl8Sgg4gIiLRp3IXEYlDKncRkTikchcRiUMqdxGROKRy\nFxGJQyp3EZE4pHIXEYlDKncRkTiUFNQL9+7d2/Py8oJ6eRGRmLR48eJd7p7d1HqBlXteXh6LFi0K\n6uVFRGKSmX3UnPW0W0ZEJA6p3EVE4pDKXUQkDqncRUTiUJPlbmaPmtlOM1txhOfNzH5jZuvNbJmZ\nTY5+TBERORbNGbn/EZh+lOdnACMiX7OA37Y+loiItEaT5e7ubwFlR1nlYuBxb/Ae0N3M+kcroIiI\nHLtonOc+ENjc6HFxZNm2KHxvEZEOwd2pCzk19SFq6sNU1zXc1tSFqa4PUVMXpqY+RHXktrY+TF3I\nqQ+Hqa0PUx926urD1IXCnD2mLxMHdW/TvO36ISYzm0XDrhtyc3Pb86VFpBOrqQ+xp6qO3VW17Kmq\nY09VLbur6thbXce+mhCVNfVU1tSzL3JbWRNi78H79VTVhqipDxGO0iWn+2alxUS5bwEGNXqcE1n2\nMe7+MPAwQH5+vq7MLSItVhcKU7K3hp17a9hRUc3OiuqD93dU1LBrX83BQq+qDR31e3VNSaRrahLd\nUpPomppE19REBnZPi9xPIj05kbTkRNKSE0hNSiQ1OYG0yG1qUgKpyYmkJiWQFrlNTWq4TU5MIDnR\nSEpMICUxgaREIynBMLM2//uJRrnPBq41s6eBE4Fyd9cuGRFptbLKWjbuqqR4dxVFpVUUlVWxeXcV\nm8v2s618/8dG0okJRu9uKfTNTKNfZhqj+2XSIz2ZHl1TyOqSTI/0FHqkJ9M9PYXu6clkdkkmPTmR\nhIS2L9v21mS5m9lTwJlAbzMrBn4MJAO4+4PAXOB8YD1QBXylrcKKSHwq31/H6m0VFOzYS8HOfazb\nsZeCHfsoraz9l/X6ZKSS2zOdE4f0JKdnOv2z0uibmUqfjDT6ZKbSq2sqiXFY1C3RZLm7+8wmnnfg\nm1FLJCJxraq2npVbK1i6eQ/ListZVryHTaVVB5/vlprEiL7dmDamLyP6dmNodldye6aT0yOdtOTE\nAJPHlsBmhRSRzmFfTT2LNpWxYGMZ7xWWsqy4nFBkf0r/rDTGD8zi0vxBHDcgk1H9MuiXmdYu+6Tj\nncpdRKLK3Vm9bS/z1+xg/pqdLI2UeVKCMXFQd75++lAm5/ZgQk4WfTLTgo4bt1TuItJq9aEw7xaW\n8vKK7cxfs5Nt5dUATMjJ4htnDOWkob2YMrgH6SmqnPaiv2kRaRF3Z0nRbmZ/uJU5y7exa18t6SmJ\nnDaiN9dPG8mZo7Ppk6GReVBU7iJyTEr21vDc4mKeXljER6VVpCYlcPaYPlw0cQBnjuqjg54dhMpd\nRJrk7izYWMaf3v2Ieau2UxdyThzSk29/YgTnHteXjLTkoCPKIVTuInJEobAzb+V2HnyrkKWb95DV\nJZkvnZzHzKm5DO/TLeh4chQqdxH5mPpQmBeWbOG3f9/Axl2VDO6Vzk8vGcdnp+Rot0uMULmLyEHh\nsDN3xTZ+PW8dhbsqmZCTxQNfmMx5x/XTJz9jjMpdRAD454Zd3D5nNSu3VjCqbwaPfCmfaWP66ANF\nMUrlLtLJbd2zn9vnrmbOsm3k9OjC3Zcdz4UTB2ikHuNU7iKdVF0ozCNvF3Lv6+sJu3P9tJF8/Yyh\n2qceJ1TuIp3Qmu0V3PjsUlZsqeDcsX350SfHMqhnetCxJIpU7iKdSF0ozINvbuA38wvITEvmwSsm\nM32cLnkcj1TuIp3E5rIqvvXUB3y4eQ+fnNCf2y4eR8+uKUHHkjaichfpBF5ZuZ2bnl2KO9w7cxIX\nThwQdCRpYyp3kThWFwrzs7mr+cM7m5iQk8V9MyeT20v71jsDlbtInNpdWcu/P7mEdwtLufKUPG45\nfzSpSToTprNQuYvEofU793LVY4vYtqeaOy+dyGem5AQdSdqZyl0kzrxdUMK/P7GE1OREnpp1ElMG\n9wg6kgRA5S4SR2Yv3coNz3zIsOxuPHrlCQzo3iXoSBIQlbtInHjsn5v4yYsrOSGvJ498KZ+sLppj\nvTNTuYvEOHfnntcLuPu1As4Z25d7Z07SFAKicheJZe7OXa8V8JvXC/jslBzu+PR4khITgo4lHYDK\nXSSGHSj2z+XncMenJ5CgmRwlQv/Fi8Sou15dp2KXI1K5i8Sg371dyD2vF3DpFBW7HJ7KXSTGvLCk\nmJ/OWc354/txx2dU7HJ4KneRGPLGmp3c9NwyThnWi7suO15XS5IjUrmLxIglRbu55snFjOmfwUNf\nnKJ5YuSoVO4iMaB4dxWzHl9En4w0/nDlVDLS9AElObpmlbuZTTeztWa23sxuPszzuWb2hpl9YGbL\nzOz86EcV6Zz21dTztccWUVMX5vdfzic7IzXoSBIDmix3M0sE7gdmAGOBmWY29pDVfgg84+6TgMuB\nB6IdVKQzCoWd7zz9Iet27OW+L0xmRN+MoCNJjGjOyH0qsN7dC929FngauPiQdRzIjNzPArZGL6JI\n5/XLV9by2uod3PrJsZwxMjvoOBJDmlPuA4HNjR4XR5Y19hPgCjMrBuYC3zrcNzKzWWa2yMwWlZSU\ntCCuSOfx8optPPj3DXz+xFy+fEpe0HEkxkTrgOpM4I/ungOcD/zJzD72vd39YXfPd/f87GyNQkSO\nZOOuSm56dhkTB3XnxxeOxUynPMqxaU65bwEGNXqcE1nW2FXAMwDu/i6QBvSORkCRzmZ/bYhrnlhM\nUqLxwBcm65RHaZHmlPtCYISZDTGzFBoOmM4+ZJ0i4GwAMxtDQ7lrv4vIMXJ3fvi/K1i7Yy93Xz6J\ngbrYhrRQk+Xu7vXAtcArwGoazopZaWa3mdlFkdVuAK42s6XAU8CV7u5tFVokXj27qJjnlxRz3dkj\ndABVWqVZU/66+1waDpQ2XnZro/urgFOjG02kcyks2cePZ6/k1OG9+PYnRgQdR2KcPqEq0gHUhcJ8\n5y8fkpqcwJ2XHq/JwKTVdLEOkQ7g7tfWsay4nN9+YTL9stKCjiNxQCN3kYAtKCzlgTc3cFn+IGaM\n7x90HIkTKneRAFVU13H9Xz5kcM90br3w0Fk9RFpOu2VEAvSzOavZXlHN89ecQtdUvR0lejRyFwnI\n2wUlPL1wM1efPpRJuT2CjiNxRuUuEoB9NfXc/PxyhmZ35fppI4OOI3FIvweKBOCOl1aztXw/z33j\nZNKSNb2ARJ9G7iLt7J8bdvHEe0V89dQhTBncM+g4EqdU7iLtaH9tiJufX05er3RuPHdU0HEkjmm3\njEg7und+AUVlVTx19Ul0SdHuGGk7GrmLtJOCHXt55O1CPjM5h5OH9Qo6jsQ5lbtIOzgwlW96ShLf\nP3900HGkE1C5i7SDF5ZsYcHGMm6eMZpe3VKDjiOdgMpdpI3tqarl9rmrmZzbncvyBzX9B0SiQOUu\n0sZ+/vIayvfXcfunxmsqX2k3KneRNrR08x6een8zXzkljzH9M4OOI52Iyl2kjbg7//niSnp3S+G6\nabqykrQvlbtIG5m9dCtLivbwvfNGk5GWHHQc6WRU7iJtoKq2nv+eu4bxA7P47JScoONIJ6RyF2kD\nD765ge0V1fz4wrE6iCqBULmLRNnmsioeequQiyYOID9PE4NJMFTuIlF2x0trMIObZ+iTqBIclbtI\nFC0oLGXO8m1cc8ZwBnTvEnQc6cRU7iJREg47P52zmgFZacw6fWjQcaSTU7mLRMmLy7ayfEs5N543\nStP5SuBU7iJRUFMf4lfz1jKmfyaXHD8w6DgiKneRaHjivSI2l+3nlhmjdeqjdAgqd5FWKt9fx73z\nCzhtRG9OH5kddBwRQOUu0mq/fXMD5fvr+I/pOvVROg6Vu0grbN2zn0ff2cinjh/IuIFZQccROahZ\n5W5m081srZmtN7Obj7DO58xslZmtNLM/RzemSMf061fXgcN3zx0ZdBSRf5HU1ApmlgjcD5wDFAML\nzWy2u69qtM4I4BbgVHffbWZ92iqwSEexelsFzy8p5urThpLTIz3oOCL/ojkj96nAencvdPda4Gng\n4kPWuRq43913A7j7zujGFOl4fv7yGjLTkvnmmcODjiLyMc0p94HA5kaPiyPLGhsJjDSzd8zsPTOb\nfrhvZGazzGyRmS0qKSlpWWKRDmDhpjLeXFvCNWcOIytdc7VLxxOtA6pJwAjgTGAm8IiZdT90JXd/\n2N3z3T0/O1unjElscnd++cpasjNS+fLJeUHHETms5pT7FqDxJdtzIssaKwZmu3udu28E1tFQ9iJx\n5+2CXby/sYxrzxquaQakw2pOuS8ERpjZEDNLAS4HZh+yzv/SMGrHzHrTsJumMIo5RToEd+dX89Yy\nsHsXLp86qOk/IBKQJsvd3euBa4FXgNXAM+6+0sxuM7OLIqu9ApSa2SrgDeAmdy9tq9AiQXll5Q6W\nFZdz3bQRpCZp1C4dl7l7IC+cn5/vixYtCuS1RVoiFHZm3PMW9WFn3ndOJylRnwGU9mdmi909v6n1\n9NMp0kwvLt3Kuh37uH7aSBW7dHj6CRVphrpQmLteW8eY/plcML5/0HFEmqRyF2mG5xYX81FpFTec\nM1JT+kpMULmLNKG6LsRvXi9gUm53zh6jmTUkNqjcRZrw5IIitpVXc9O5ozDTqF1ig8pd5Cgqa+p5\n4I31nDKsF6cM7x10HJFmU7mLHMUf/7mJ0spabjxvVNBRRI6Jyl3kCMqr6njo7xs4e3QfJuf2CDqO\nyDFRuYscwSNvF1JRXc8N52rULrFH5S5yGKX7anj0nY1cMKE/YwdkBh1H5Jip3EUO46G3CqmuC3H9\nNF0+T2KTyl3kEDsrqnnsn5u45PiBDO/TLeg4Ii2ichc5xANvbqA+7Fw3TZckkNilchdpZMue/fx5\nQRGXTslhcK+uQccRaTGVu0gj981fj+Nc+wld9Fpim8pdJKKotIpnF21m5tRccnqkBx1HpFVU7iIR\n97xeQGKC8c2zNGqX2KdyFwE2lOzjfz4o5oqTBtM3My3oOCKtpnIXAe55rYDUpESuOXNY0FFEokLl\nLp3e2u17eXHZVq48NY/e3VKDjiMSFSp36fTuenUdXVOSmHXa0KCjiESNyl06tRVbynl55Xau+rch\n9OiaEnQckahRuUun9utX15HVJZmrThsSdBSRqFK5S6e1pGg389fsZNbpQ8lMSw46jkhUqdyl0/r1\nvHX07JrClafkBR1FJOpU7tIpLSgs5R/rd3HNGcPompoUdByRqFO5S6fj7tw5bx3ZGalccdLgoOOI\ntAmVu3Q6bxfs4v1NZXzzzGF0SUkMOo5Im1C5S6cSDju/eGUNA7t3YeaJuUHHEWkzKnfpVOau2MaK\nLRV895yRpCZp1C7xq1nlbmbTzWytma03s5uPst5nzMzNLD96EUWioy4U5s556xjZtxuXTBoYdByR\nNtVkuZtZInA/MAMYC8w0s7GHWS8DuA5YEO2QItHw3OJiNu6q5MZzR5GYYEHHEWlTzRm5TwXWu3uh\nu9cCTwMXH2a9/wJ+DlRHMZ9IVFTXhbjntQIm5XbnnLF9g44j0uaaU+4Dgc2NHhdHlh1kZpOBQe4+\nJ4rZRKLm8Xc3sb2imu+dNxozjdol/rX6gKqZJQC/Bm5oxrqzzGyRmS0qKSlp7UuLNEtFdR0PvLmB\n00dmc/KwXkHHEWkXzSn3LcCgRo9zIssOyADGAW+a2SbgJGD24Q6quvvD7p7v7vnZ2dktTy1yDB55\nq5A9VXV877xRQUcRaTfNKfeFwAgzG2JmKcDlwOwDT7p7ubv3dvc8d88D3gMucvdFbZJY5BiU7K3h\n9//YyAUT+jNuYFbQcUTaTZPl7u71wLXAK8Bq4Bl3X2lmt5nZRW0dUKQ17ptfQE19mBvOGRl0FJF2\n1awZk9x9LjD3kGW3HmHdM1sfS6T1Ckv28eSCIi47YRBDs7sFHUekXekTqhK3fv7yGlKTErh+mkbt\n0vmo3CUuvb+xjFdW7uCaM4eRnaGLXkvno3KXuBMOO7fPWUW/zDSu+jdd9Fo6J5W7xJ2/Ld/G0uJy\nbjxvlKb0lU5L5S5xpbouxM9fWsPY/pl8SpODSSemcpe48vi7m9iyZz8/uGCMJgeTTk3lLnGjrLKW\ne+ev56xR2Zw6vHfQcUQCpXKXuHHnvLVU1Ya45fwxQUcRCZzKXeLCii3l/Pn9Ir508mBG9s0IOo5I\n4FTuEvPcnf98cSU90lP4jj6wJAKo3CUOzF66lYWbdvO980aR1SU56DgiHYLKXWJaZU09/z13DeMH\nZnFp/qCm/4BIJ9GsicNEOqoH3lzP9opq7v/CJJ36KNKIRu4Ssz4qreSRtzby6UkDmTK4Z9BxRDoU\nlbvEJHfnR39dSXKi8R8zRgcdR6TDUblLTPrbsm28ta6Em84bRd/MtKDjiHQ4KneJOeX76/jPF1cx\nISeLL56cF3QckQ5JB1Ql5vzi5TWUVdbwx6+coIOoIkegkbvElMUf7ebJBUV89dQhuuC1yFGo3CVm\n1IXCfP+F5QzISuN6XfBa5Ki0W0ZixkN/38DaHXv53Zfy6ZqqH12Ro9HIXWLCmu0V3PN6ARdM6M+0\nsX2DjiPS4ancpcOrC4W58dmlZHVJ5r8uHhd0HJGYoN9tpcP77ZsbWLGlggevmEzPrilBxxGJCRq5\nS4e2amsF984v4MKJA5g+rn/QcURihspdOqza+gO7Y1K47aLjgo4jElO0W0Y6rDvnrWXVtgoe/uIU\nemh3jMgx0chdOqS3C0p46K1CPn9iLuce1y/oOCIxR+UuHU7pvhq++8xShvfpxo8uGBt0HJGYpN0y\n0qG4Ozc9t4zyqjoe+8pUuqQkBh1JJCZp5C4dymP/3MT8NTu5ecZoxg7IDDqOSMxqVrmb2XQzW2tm\n683s5sM8/10zW2Vmy8zsdTMbHP2oEu8+KNrN7XNXc9aobL5yal7QcURiWpPlbmaJwP3ADGAsMNPM\nDt0R+gGQ7+4TgOeAX0Q7qMS30n01/PuTS+ibmcZdlx2PmabyFWmN5ozcpwLr3b3Q3WuBp4GLG6/g\n7m+4e1Xk4XtATnRjSjwLhZ1vP/0BpZW1PHjFFLqn67RHkdZqTrkPBDY3elwcWXYkVwEvtSaUdC53\nzlvLO+tL+ekl4zRHu0iURPVsGTO7AsgHzjjC87OAWQC5ubnRfGmJUS8u3coDb25g5tRBfC5/UNBx\nROJGc0buW4DG77qcyLJ/YWbTgB8AF7l7zeG+kbs/7O757p6fnZ3dkrwSRz4o2s2Nzy7lhLwe/ETT\nC4hEVXPKfSEwwsyGmFkKcDkwu/EKZjYJeIiGYt8Z/ZgSb7bs2c/Vjy+mb2YaD30xn9Qknc8uEk1N\nlru71wPXAq8Aq4Fn3H2lmd1mZhdFVvsl0A141sw+NLPZR/h2IuyrqeeqPy6kpj7Eo1fmaxpfkTbQ\nrH3u7j4XmHvIslsb3Z8W5VwSp2rqQ1zzxGIKdu7jD1eewPA+GUFHEolL+oSqtJtQ2PnuX5bydsEu\n7vj0eE4fqeMuIm1F5S7twt259a8rmLN8Gz84fwyX6swYkTalcpc25+7cOW8dTy4o4htnDOPq04cG\nHUkk7qncpU25O3e9uo773ljP5ScM4j+mjwo6kkinoCl/pc0cGLEfKPaffWq85owRaScqd2kT7s4v\nX1l78NOnt18ynoQEFbtIe1G5S9SFws6PZ6/gifeK+PyJufz04nEqdpF2pnKXqKquC3H9Xz7kpRXb\n+foZQ7l5+mjtihEJgMpdoqaiuo5Zjy/ivcIyfnjBGL52ms6KEQmKyl2iYuOuSr722EI+Kq3i7suO\n55JJR5sVWkTamspdWu2tdSVc++clJCYYf7rqRE4e1ivoSCKdnspdWszd+f0/NvKzuasZ2TeDR76U\nz6Ce6UHHEhFU7tJCuytruem5pby2eifTj+vHnZ+bSNdU/TiJdBR6N8oxe39jGdc9/QGl+2r58YVj\nufKUPJ0RI9LBqNyl2arrQtz9WgEPv7WBQT3Tef6aUxifo2ueinREKndpliVFu7np2aVsKKnksvxB\n/PCTY8hISw46logcgcpdjmpvdR13v1bAH97ZSL/MNB776lTO0DzsIh2eyl0OKxx2XvhgC3e8tIbS\nyhpmTs3llhmjNVoXiREqd/mYxR/t5qdzVvFB0R6OH9Sd3385n4mDugcdS0SOgcpdDlq5tZw7561j\n/pqd9O6Wyq8uncinJw3UpF8iMUjlLqzYUs5v39zAnOXbyExL4qbzRnHlKXk6b10khund20m5O39f\nV8IjbxfyzvpSuqUm8a1PDOdrpw0lq4v2q4vEOpV7J1O+v46/friFJ98rYu2OvfTNTOXmGaP5/Im5\nZOpgqUjcULl3Au7O4o9289T7m5mzfCvVdWGOG5DJLz87gYuPH0hKki6lKxJvVO5xyt1ZubWCvy3b\nxpzlW9lctp9uqUl8enIOM0/I1SdLReKcyj2O1IfCLCnawxtrd/Lyiu1s3FVJYoJx6vDefPsTIzh/\nfH8dJBXpJPROj3Fb9uzn3Q2lvLl2J2+tK6Giup6kBGPqkJ5cfdpQpo/rR8+uKUHHFJF2pnKPIe7O\nxl2VLNxUxoLCMhZsLGPLnv0AZGekMn1cP84a1YdTR/TWwVGRTk7l3kGFw05RWRXLt5SzYkv5wduK\n6noAendLiYzOhzB1SC9G98vQh41E5CCVe8DqQ2E+Kqti/c59rN+5jw0797G+pOG2sjYEQEpiAqP7\nZ3DBhAFMyMnihLyeDMvuqjnUReSIVO5trLouRMneGrbu2c/m3fsp3l1FcaPbbeXVhMJ+cP1+mWkM\n79ONS/MHMbpfBuMGZjGyb4ZOVxSRY9Kscjez6cA9QCLwO3e/45DnU4HHgSlAKXCZu2+KbtTguTtV\ntSHK99dRvr+OPVUNtxX769hVWUPJ3kZf+xpu90Z2oxxgBn0z0sjp0YX8wT0Y2KMLQ3p3Y3ifbgzL\n7qpZF0UkKposdzNLBO4HzgGKgYVmNtvdVzVa7Spgt7sPN7PLgZ8Dl7VF4KNxd2pDYWrrw9SFnNr6\nhvu1oRC19f//XFVtPVW1ochX5H5N5Lau0f3aEJW19Q1lHiny+kaj7ENlpCaRnZFK74xUxvTP5PQR\nqWRnpJLdLZUB3buQ06ML/bunkZqU2I5/KyLSGTVn5D4VWO/uhQBm9jRwMdC43C8GfhK5/xxwn5mZ\nux+5CVvoLwuLePitwoNFfeCrLtRQ3q2RnpIY+UoiPSWRLimJdE1JYkBWF7LSk8nq0vDVvcv/3z+w\nvFfXVLqkqLRFpGNoTrkPBDY3elwMnHikddy93szKgV7ArsYrmdksYBZAbm5uiwL3SE9hdL9MUpIS\nSElMaLiNfCUnJpB6yPLkA/cPPJeUQJdIiXdNSTp4Py0pUWebiEjcaNcDqu7+MPAwQH5+fotG9ece\n149zj+sX1VwiIvGmOadgbAEGNXqcE1l22HXMLAnIouHAqoiIBKA55b4QGGFmQ8wsBbgcmH3IOrOB\nL0fufxaY3xb720VEpHma3C0T2Yd+LfAKDadCPuruK83sNmCRu88Gfg/8yczWA2U0/AcgIiIBadY+\nd3efC8w9ZNmtje5XA5dGN5qIiLSUPvYoIhKHVO4iInFI5S4iEodU7iIicciCOmPRzEqAj1r4x3tz\nyKdfY5i2peOJl+0AbUtH1ZptGezu2U2tFFi5t4aZLXL3/KBzRIO2peOJl+0AbUtH1R7bot0yIiJx\nSOUuIhKHYrXcHw46QBRpWzqeeNkO0LZ0VG2+LTG5z11ERI4uVkfuIiJyFDFd7mb2LTNbY2YrzewX\nQedpLTO7wczczHoHnaUlzOyXkX+PZWb2P2bWPehMx8rMppvZWjNbb2Y3B52npcxskJm9YWarIu+P\n64LO1BpmlmhmH5jZ34LO0hpm1t3Mnou8T1ab2clt9VoxW+5mdhYNl/eb6O7HAb8KOFKrmNkg4Fyg\nKOgsrfAqMM7dJwDrgFsCznNMGl0veAYwFphpZmODTdVi9cAN7j4WOAn4ZgxvC8B1wOqgQ0TBPcDL\n7j4amEgbblPMljtwDXCHu9cAuPvOgPO01l3A94CYPQji7vPcvT7y8D0aLuwSSw5eL9jda4ED1wuO\nOe6+zd2XRO7vpaFEBgabqmXMLAe4APhd0Flaw8yygNNpmCIdd6919z1t9XqxXO4jgdPMbIGZ/d3M\nTgg6UEuZ2cXAFndfGnSWKPoq8FLQIY7R4a4XHJOF2JiZ5QGTgAXBJmmxu2kY+ISDDtJKQ4AS4A+R\nXUy/M7OubfVi7XoN1WNlZq8Bh7tg6g9oyN6Thl85TwCeMbOhHfUKUE1sy/dp2CXT4R1tO9z9r5F1\nfkDDboEn2zObfJyZdQOeB77j7hVB5zlWZvZJYKe7LzazM4PO00pJwGTgW+6+wMzuAW4GftRWL9Zh\nufu0Iz1nZtcAL0TK/H0zC9MwX0NJe+U7FkfaFjMbT8P/6EvNDBp2ZSwxs6nuvr0dIzbL0f5NAMzs\nSuCTwNkd9T/ao2jO9YJjhpkl01DsT7r7C0HnaaFTgYvM7HwgDcg0syfc/YqAc7VEMVDs7gd+g3qO\nhnJvE7G8W+Z/gbMAzGwkkEIMTirk7svdvY+757l7Hg0/AJM7YrE3xcym0/Dr80XuXhV0nhZozvWC\nY4I1jBR+D6x2918Hnael3P0Wd8+JvDcup+H6zLFY7ETe05vNbFRk0dnAqrZ6vQ49cm/Co8CjZrYC\nqAW+HIMjxXhzH5AKvBr5LeQ9d/9GsJGa70jXCw44VkudCnwRWG5mH0aWfT9yyUwJzreAJyODh0Lg\nK231QvqEqohIHIrl3TIiInIEKncRkTikchcRiUMqdxGROKRyFxGJQyp3EZE4pHIXEYlDKncRkTj0\nfwN9aSXmuHWfAAAAAElEQVQ1rwYeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bdc54a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = T.dvector('x')\n",
    "s = 1 / (1 + T.exp(-x))\n",
    "logistic = function([x], s)\n",
    "\n",
    "# Build a numpy vector for the support of the function\n",
    "x_eval = np.linspace(-6,6,1000)\n",
    "\n",
    "# Evaluate the function. Note that Thenao understands numpy arrays.\n",
    "y = logistic(x_eval)\n",
    "\n",
    "# Plot\n",
    "plt.plot(x_eval,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shared variables\n",
    "Shared variables are hybrid symbolic and non-symbolic variables whose value may be shared between multiple functions. They are symbolic because they can be used just as the elements returned by `dvectors` or `dmatrices`, but they also adopt an internal value that defines the value the symbolci variable will have in all functions that share it. As we will see, shared values are useful when performing training in supervised learning algorithms, because we want these values to remain around across many executions of a Theano training or testing function. Often, the purpose of a Theano training function is to update the weights stored in a shared variable. And a testing function needs the current weights to perform the network's forward pass. A good example of how to use shared variables is by building an accumulator function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from theano import shared\n",
    "\n",
    "# We create a shared valuable and initiazlize it with value 0.\n",
    "state = shared(0)\n",
    "\n",
    "# Define a scalar variable with the increment\n",
    "inc = T.iscalar('inc')\n",
    "\n",
    "# Everytime the function is executed, our shared valuable is updated and returns old value\n",
    "accumulator = function([inc], state, updates=[(state, state+inc)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array(0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let us run it:\n",
    "print(state.get_value())\n",
    "accumulator(1)\n",
    "# get_value extracts the current value of the shared variable.\n",
    "# you can also reset this value with set_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array(1)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(state.get_value())\n",
    "accumulator(300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "301\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array(301)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(state.get_value())\n",
    "accumulator(300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiation in Theano\n",
    "One reason why Theano is so powerful is because it uses Automatic Differentiation (autdiff) in order to compute gradients. Autdiff is a set of techniques to numerically evaluate the derivative of a function defined by a computer program. This implies that Theano can compute the derivative of any function for which you can write code (not only for parametric functions, but also data strcutures). Autdiff is different from numerical differentiation (which uses finite differences) and also from symbolic differentiation (which is a way to write new computer expressions for the derivative of a given function). Both numerical and symbolic differentiation can be very inefficient when calculating partial derivatives with respect to many inputs. Autdiffs solves for these problems.\n",
    "\n",
    "At the core of autdiff is the concept of dual numbers. We can extend any natural number to a dual number:\n",
    "\n",
    "$$x \\rightarrow x + \\epsilon x'$$\n",
    "\n",
    "$\\epsilon$ is just a symbol that distinguishes the second component, and we define its arithmetic as: $\\epsilon^2 = 0$. With this, it is easy to see that:\n",
    "\n",
    "$$(x + \\epsilon x') + (y + \\epsilon y') = x + y + \\epsilon(x' + y')$$\n",
    "$$(x + \\epsilon x')(y + \\epsilon y') = xy + \\epsilon(xy' + x'y)$$\n",
    "$$\\frac{1}{x+\\epsilon x'} = \\frac{1}{x}-\\epsilon\\frac{x'}{x^2}$$\n",
    "\n",
    "We can then apply these rules to a polinomial:\n",
    "\n",
    "$$p(x) = p_0 + p_1x + p_2x^2 + ... + p_nx^n$$\n",
    "\n",
    "to obtain:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "p(x + \\epsilon x') & = p_0 + p_1(x + \\epsilon x' ) + ... + p_n(x + \\epsilon x')^n \\\\\n",
    "                   & = p_0 + p_1x + p_2x^2 + ... + p_nx^n \\\\\n",
    "                   &   \\, \\, \\, \\, + \\epsilon p_1x + 2\\epsilon p_2 xx' + ... + n\\epsilon p_nx^{n−1}x' \\\\\n",
    "                   & = P(x) + \\epsilon P'(x)x' \\\\\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "We are free to choose $x' = 1$ and then **a function applied to a dual number returns its derivative in the second dual component**. Dual numbers are data structures that allow us to carry the derivatives around without too much additional computation or errors. Theano profits from autdiff to compute complicated gradients in many dimensions using `T.grad`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derivative of the logistic function\n",
    "We can also compute the gradient of complex expressions such as the logistic function defined above. It turns out that the derivative of the logistic is: $ds(x)/dx = s(x) \\cdot (1 - s(x))$. But we do not even need to know its symbolic expression to evaluate it at a given point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.25      ,  0.19661193],\n",
       "       [ 0.19661193,  0.10499359]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a Theno graph that computes the logistic function, \n",
    "# and plot the derivtaive for the range x in [-6,6]\n",
    "x = T.dmatrix('x')\n",
    "s = T.sum(1/(1+ T.exp(-x)))\n",
    "gs = T.grad(s, x)\n",
    "\n",
    "# create the callable function \n",
    "dlogistic = theano.function([x], gs)\n",
    "dlogistic([[0,1],[-1,-2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10c075b38>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJ5PJvm9khbAEJOwSoIDgggsIBVuh4tKq\ntbXaatvb/tra621rvUtr7aat916pS29brYprWFxQqFoVTFgSCGtYsy8kJCEh+/f3RyY2poFMkknO\nLJ/n45GHyZkzM+8xzHtOvud7zhFjDEoppXyDn9UBlFJKjRwtfaWU8iFa+kop5UO09JVSyodo6Sul\nlA/R0ldKKR+ipa+UUj5ES18ppXyIlr5SSvkQf6sD9BYXF2fS09OtjqGUUh5l586d1caY+P7Wc7vS\nT09PJzc31+oYSinlUUTkpDPr6fCOUkr5EC19pZTyIVr6SinlQ7T0lVLKh2jpK6WUD3Gq9EVkqYgc\nEpFCEbmvj9u/IyL7RSRfRN4RkTE9busQkT2Or2xXhldKKTUw/U7ZFBEb8BhwFVAM5IhItjFmf4/V\ndgNZxpgmEbkb+AVwg+O2c8aYmS7OrZRSahCcmac/Fyg0xhwDEJHngFXAJ6VvjNnWY/3twC2uDKmU\nlcrqzvHR0dOU1TVjtwnj4sJYMCGWkAC3O8xFqX458682BSjq8XMxMO8C698BvN7j5yARyQXagZ8b\nY17tfQcRuRO4E2D06NFORFJq+B0oq+fhNw+x9WDlP90WbLdx07zR3HP5BKJDAyxIp9TguHRTRURu\nAbKAS3ssHmOMKRGRccBWEdlrjDna837GmHXAOoCsrCy9UruylDGG328t5LfvHCEiyJ9vLclg6dRE\nxsaF0trRyb7iOl7cWcwfPzxBdl4pv/7CDBZl9Hv0u1JuwZnSLwHSevyc6lj2KSJyJXA/cKkxpqV7\nuTGmxPHfYyLyN2AWcLT3/ZVyB63tnXx3fR4b8kpZOSOZB1dNISrkH1vyQXYbCybEsWBCHHcsGsu/\nPL+H257O4Wefm8YX5qRd4JGVcg/OzN7JATJEZKyIBABrgU/NwhGRWcDjwEpjTGWP5dEiEuj4Pg5Y\nSI99AUq5k/aOTr79/G425JXy/aWTeGTtzE8Vfm9TkiN56e4FLBgfy/dfyuf5nFMjmFapwem39I0x\n7cA9wJvAAeAFY0yBiDwoIisdqz0MhAHre03NnAzkikgesI2uMX0tfeWWHthQwOa95fzb8sl8/bIJ\niEi/9wkPsvPkrXNYPDGeH768l7f3V4xAUqUGT4xxryH0rKwso2fZVCPtxZ3F/L/1eXxt8Th+eO3k\nAd+/saWdG/+wneNVjWy49xLS40KHIaVS5yciO40xWf2tp0fkKp93oKye+1/Zy/xxsXzvmkmDeozQ\nQH/+55bZ2GzC3c/sormtw8UplXINLX3l09o6OvnuC3mEB9l59MZZ+NsG/5ZIiQrmV2tmcKCsnkff\nOeLClEq5jpa+8mnr3jvG/rJ6/uO6qcSHBw758ZZMHsXq2ak8/t4x9pXUuSChUq6lpa981tGqszzy\n9hGunZbI0qmJLnvcHy3PJCY0gB+8lE97R6fLHlcpV9DSVz7rvzYdINDfj5+unOrSx40MsfOTz2ZS\nUFrPC7nFLn1spYZKS1/5pA8Lq3nnYCVfv3yCS4Z1els+LYk56dH8esshGprbXP74Sg2Wlr7yOR2d\nhv/YdICUqGBuX5g+LM8hIvzb8kyqz7by2DY9AF25Dy195XNe21PC/rJ6vr90EkF227A9z4y0KD43\nK4WnPzhOZUPzsD2PUgOhpa98SntHJ7/bWkhmUgQrZyQP+/N9a0kG7Z2Gx989NuzPpZQztPSVT9mY\nX8bx6ka+uSTDqdMsDFV6XCifm5XCX7af1K195Ra09JXP6Og0/G7rES5KDOfqzFEj9rz3XD5Bt/aV\n29DSVz5j894yjlY1cu8VGfj5Df9Wfrf0uFCum9m1tV/T2Dpiz6tUX7T0lU8wxvD4e0cZHx/KMhce\niOWsuy4dR0t7J3/ZfnLEn1upnrT0lU/4+HgN+0rqueOScSO6ld8tY1Q4l02K508fndCTsSlLaekr\nn/Dk348TFWLnc7NSLMvw1UXjqD7bSnZeqWUZlNLSV17v5OlGthyo4OZ5owkOGL55+f1ZMD6WixLD\nefL947jbdSyU79DSV17v6Q9O4O8nfGl+uqU5RISvLBrHoYoG3j9SbWkW5bu09JVXa2huY31uESum\nJzMqIsjqOKyckUxcWAB/1h26yiJa+sqrZeeV0tjawZfmj7E6CgAB/n6syUpj68FKyurOWR1H+SAt\nfeW1jDE8u+MUk5MimJkWZXWcT9w4ZzQdnYbnc4qsjqJ8kJa+8lr5xXUUlNZz09y0ETnlgrNGx4aw\neGI8z+cU6UVW1IjT0lde668fnyLYbmOVhdM0z+emuaMpq2tm26Eqq6MoH6Olr7xSQ3Mb2XmlrJyR\nTESQ3eo4/2TJ5ARGRQTyzA7doatGlpa+8kqv7SmlqbWDG+eNtjpKn+w2P27ISuPdw1UU1zZZHUf5\nEC195ZWey+nagTsjNdLqKOe1JisNY+DlXSVWR1E+REtfeZ1D5Q3sK6nnC1mpbrUDt7e0mBDmj4vl\npV3FeoSuGjFa+srrvLyrGH8/GZErYw3V9bNTOXm6idyTtVZHUT5CS195lY5Owyu7S7hsUgKxYYFW\nx+nXsqmJhATYeGlnsdVRlI/Q0lde5e+F1VQ2tHD9xe43TbMvoYH+XDstiY35ZZxr1VMuq+Gnpa+8\nysu7iokMtnPF5ASrozjt+otTOdvSzlv7y62OonyAlr7yGg3NbbxZUM5nZyQR6G/dKZQHat7YGFKj\ng3lRh3jUCHCq9EVkqYgcEpFCEbmvj9u/IyL7RSRfRN4RkTE9brtVRI44vm51ZXilenp9bznNbZ1c\nf3Gq1VEGxM9P+PzFqfy9sFpPwqaGXb+lLyI24DFgGZAJ3Cgimb1W2w1kGWOmAy8Cv3DcNwb4CTAP\nmAv8RESiXRdfqX94cVcx4+JC3erkas76/KwUjIENelUtNcyc2dKfCxQaY44ZY1qB54BVPVcwxmwz\nxnQfVrgd6N7UugbYYoypMcbUAluApa6JrtQ/FNc28fHxGj5/cYpbz80/n/S4UGakRuqlFNWwc6b0\nU4Ce54Atdiw7nzuA1wdyXxG5U0RyRSS3qkpPQKUGbvPeMgBWzvCMWTt9WTkzhX0l9RytOmt1FOXF\nXLojV0RuAbKAhwdyP2PMOmNMljEmKz4+3pWRlI/YlF/G9NRIRseGWB1l0FZMT0IEsvfo1r4aPs6U\nfgmQ1uPnVMeyTxGRK4H7gZXGmJaB3FepoSiqaSKvuI7l05KsjjIkoyKC+MzYWLLzSvW0DGrYOFP6\nOUCGiIwVkQBgLZDdcwURmQU8TlfhV/a46U3gahGJduzAvdqxTCmX2eQY2rnWw0sfYNXMZI5XN7Kv\npN7qKMpL9Vv6xph24B66yvoA8IIxpkBEHhSRlY7VHgbCgPUiskdEsh33rQH+na4PjhzgQccypVxm\nU34ZM9KiSIvx3KGdbsumJmG3Ca/t0T+I1fDwd2YlY8xmYHOvZT/u8f2VF7jvU8BTgw2o1IWcPN3I\n3pI67r92stVRXCIyxM6lExPYkF/KD6+djM3P82YiKfemR+Qqj9Y9tLNsWqLFSVxn5cxkKupb+Pi4\n/lGsXE9LX3m0TfllzBodRWq05w/tdLtycgIhATads6+GhZa+8ljHqxspKK1nxXT3P2/+QIQE+HNV\n5ihe31dGW0en1XGUl9HSVx5r8yezdrxnaKfb8mlJnGlq48Ojp62OoryMlr7yWBvySskaE01SZLDV\nUVxu8cR4wgL92ZxfZnUU5WW09JVHKqw8y8HyBpZP9/y5+X0Jstu4cnICb+4v1yEe5VJa+sojbd5b\nhkjXvHZvda0O8ahhoKWvPNKm/DLmjIkhMTLI6ijDpnuIZ1O+zuJRrqOlrzzOkYoGDlV479BOt+4h\nnrf2V+gQj3IZLX3lcTZ9MrTjfbN2els+PVmHeJRLaekrj7Mpv4y56TEkRHjv0E63RRlxOsSjXEpL\nX3mUwxUNHKk8ywovH9rpFmS3cVXmKN4s0CEe5Rpa+sqjbMwvw0/gGh8Y2ul27bQk6s618UFhtdVR\nlBfQ0lcewxjDpvxS5o2NJSHc+4d2ui3KiCM80P+TI5CVGgotfeUxDlU0cLSq0etn7fQWZLdxpQ7x\nKBfR0lceY5NjaGepDw3tdFuuQzzKRbT0lUfoGtopY/74WOLCAq2OM+IWTdQhHuUaWvrKIxwoa+BY\ndSPLp3nXaZSdFeivs3iUa2jpK4+waW8pNj/hmimjrI5iGZ3Fo1xBS1+5ve6hnQXjY4n1waGdbt1D\nPJv0dMtqCLT0ldsrKK3nxOkmlk/zrVk7vXUP8by1v4LWdh3iUYOjpa/c3qa9ZY6hHd+btdPb8umO\nIZ6jOsSjBkdLX7m17qGdhRPiiA4NsDqO5S7JiCM8SK+opQZPS1+5tX0l9ZyqaWKFjw/tdPvHLJ5y\nHeJRg6Klr9zaxvxS/P2Eq3141k5vK6YnUd/crrN41KBo6Su3ZYxhY34Zl2TEERWiQzvdLpkQT3iQ\nP5v0QC01CFr6ym3lFddRcuacz8/a6S3A34+rMxN1iEcNipa+club8kux24SrM3XWTm8rpifR0NzO\n3wurrI6iPIyWvnJL3bN2FmXEExlitzqO21k4IY6IIH825ZdbHUV5GC195ZZ2F52htK5Zh3bOI8Df\nj6unJPLW/nJa2jusjqM8iFOlLyJLReSQiBSKyH193L5YRHaJSLuIrO51W4eI7HF8ZbsquPJum/LL\nCLD5cZXO2jmv5d1DPEd0Fo9yXr+lLyI24DFgGZAJ3Cgimb1WOwXcBjzbx0OcM8bMdHytHGJe5QM6\nOw2b95axeGI8EUE6tHM+C8fHERls13PxqAFxZkt/LlBojDlmjGkFngNW9VzBGHPCGJMP6FQCNWS7\nTtVSVtfM8um6A/dCumbxjGLL/god4lFOc6b0U4CiHj8XO5Y5K0hEckVku4hcN6B0yidtzC8jwN+P\nKyfr0E5/lk9PoqGlnfcP6xCPcs5I7MgdY4zJAm4Cfisi43uvICJ3Oj4YcquqdAqaL+voNGzaW8bl\nk+IJ16Gdfi2c4Bji0QO1lJOcKf0SIK3Hz6mOZU4xxpQ4/nsM+Bswq4911hljsowxWfHx8c4+tPJC\nHx+voaqhhRXTffMKWQNlt/lxzZRRvL2/guY2HeJR/XOm9HOADBEZKyIBwFrAqVk4IhItIoGO7+OA\nhcD+wYZV3m9jfinBdhtLJidYHcVjLJ+e3DXEo7N4lBP6LX1jTDtwD/AmcAB4wRhTICIPishKABGZ\nIyLFwBrgcREpcNx9MpArInnANuDnxhgtfdWn9o5O3thXzhWTEwgJ8Lc6jsdYMD6WqBA7m/JLrY6i\nPIBT7yxjzGZgc69lP+7xfQ5dwz697/chMG2IGZWP+OjYaU43tvLZ6XpA1kDYbX5ck5nIpr1lNLd1\nEGS3WR1JuTE9Ile5jY15ZYQG2Lhskg7tDNTy6UmcbWnnvcM6EUJdmJa+cgut7Z28UVDOVZmjdEt1\nEOaPjyU6RGfxqP5p6Su38EFhNXXn2nTWziB1zeJJ1Fk8ql9a+sotbMgvJTzIn0UT46yO4rGWT0+i\nsbWDd3WIR12Alr6yXHNbB1sKKrhmSiKB/jq0M1jzxzmGePRcPOoCtPSV5d47XEVDSzsrdNbOkPjb\n/Fg6NYm3D1TQ1NpudRzlprT0leU25pcRHWJn4QQd2hmqVTOTaWrt4O0DlVZHUW5KS19Z6lxrB28f\nqGDp1ETsNv3nOFRz02NIigzitd1OnylF+Rh9lylLdQ1FdPBZnbXjEn5+wsoZybx7uIqaxlar4yg3\npKWvLPXq7hISI4L4zLhYq6N4jVUzU2h3nK1Uqd609JVlahpbefdwFStnJuPnJ1bH8RqTk8LJSAgj\ne48O8ah/pqWvLLNpbxntnYbrZg7kmjyqPyLCdbNSyDlRS3Ftk9VxlJvR0leWeXV3CRNHhTE5Kdzq\nKF5n5YyufSSv7dEzb6pP09JXljh1uomdJ2tZNTMFER3acbW0mBCyxkSTraWvetHSV5Z4zTHevGqm\nztoZLqtmJnOoooEDZfVWR1FuREtfjThjDK/uKWFuegyp0SFWx/Fay6cn4+8nvKo7dFUPWvpqxBWU\n1nO0qpHrZukO3OEUExrA4onxbNhTSmensTqOchNa+mrEvbK7BLtNuHZaotVRvN6qmcmU1jWTc6LG\n6ijKTWjpqxHV0WnYkFfKZZMSiAoJsDqO17sqcxShATZe3qVDPKqLlr4aUe8fqaKyoYXP6dDOiAgJ\n8OfaaUlszC/VM28qQEtfjbD1O4uJCrGzZLJeB3ekrMlKo7G1gzf2lVsdRbkBLX01Ys40tbKloILr\nZqboxVJG0Jz0aMbEhrA+t9jqKMoNaOmrEZOdV0prRyerZ6daHcWniAirL07lo2OnKarR0zL4Oi19\nNWJe3FnM5KQIpqZEWh3F53x+dioi8NIu3dr3dVr6akQcKm8gv7iONbqVb4mUqGAWjo/jxZ3FOmff\nx2npqxGxPrcIu030gCwLrZ6dSnHtOXYc1zn7vkxLXw27to5OXt1TwpKLRhETqnPzrXLNlETCA/1Z\nv7PI6ijKQlr6atj97VAV1WdbWZOlQztWCg6wsWJGEq/vLedsi87Z91Va+mrYvZBbRFxYIIsnxlsd\nxeetnp3GubYONubpKZd9lZa+Glbldc1sPVjJ6tmp2G36z81qF4+OIiMhjL9+fMrqKMoi+i5Uw+r5\nnCI6Og03zR1tdRRF15z9m+aNJq+4jn0ldVbHURZwqvRFZKmIHBKRQhG5r4/bF4vILhFpF5HVvW67\nVUSOOL5udVVw5f7aOzp5LucUizLiGB2r5813F5+flUqQ3Y9ndujWvi/qt/RFxAY8BiwDMoEbRSSz\n12qngNuAZ3vdNwb4CTAPmAv8RESihx5beYK/HaqirK6Zm+fpVr47iQyx89npyWTvKdEduj7ImS39\nuUChMeaYMaYVeA5Y1XMFY8wJY0w+0NnrvtcAW4wxNcaYWmALsNQFuZUHePbjU8SHB7Jk8iiro6he\nbpo3msbWDl7dradc9jXOlH4K0HNib7FjmTOGcl/lwYprm9h2qJIbstJ0B64bmpkWRWZSBM/sOIUx\neoSuL3GLd6OI3CkiuSKSW1VVZXUc5QLP53R91q+dm2ZxEtWX7h26B8rq2VN0xuo4agQ5U/olQM93\nbqpjmTOcuq8xZp0xJssYkxUfr3O5PV1bRyfP5xRx6cR4vfC5G7tuVgqhATae1R26PsWZ0s8BMkRk\nrIgEAGuBbCcf/03gahGJduzAvdqxTHmxN/aVU9nQwi3zxlgdRV1AWKA/K2emsCG/lDNNrVbHUSOk\n39I3xrQD99BV1geAF4wxBSLyoIisBBCROSJSDKwBHheRAsd9a4B/p+uDIwd40LFMebGnPzjOmNgQ\nrrhIr47l7r40fwzNbZ08l6Pn4/EV/s6sZIzZDGzutezHPb7PoWvopq/7PgU8NYSMyoPkFZ1h16kz\n/GhFJn5+YnUc1Y/JSRHMHxfLnz48wVcuGYu/7nT3evobVi71xw9PEBpg05OreZDbF6ZTWtfMmwUV\nVkdRI0BLX7lMZX0zG/NLWZOVRkSQ3eo4yklLJo9idEwIT39w3OooagRo6SuXeWbHKdo6DF+arztw\nPYnNT7h1QTq5J2vJL9bpm95OS1+5REt7B8/sOMnlk+IZFx9mdRw1QF/ISiUs0J+nPzhhdRQ1zLT0\nlUtsyCuj+mwrty8ca3UUNQjhQXZWz05lY34plfXNVsdRw0hLXw1ZZ6dh3XtHmTQqnEUZcVbHUYN0\n24J02jsNf/zwhNVR1DDS0ldDtu1QJYcrzvK1S8chotM0PVV6XCjLpiby5+0naWhuszqOGiZa+mrI\n/udvR0mJCuazM5KtjqKG6O5LJ9DQ3K7n2vdiWvpqSHJP1JB7spavLBqrZ9P0AtNSI1mUEceTfz9O\nc1uH1XHUMNB3qRqS/333KNEhdm6Yo2fT9BZ3XzqeqoYWXtpVbHUUNQy09NWgHalo4O0Dldy6IJ2Q\nAKfO6KE8wPzxscxIjeTxd4/R3tH7ukjK02npq0H7/bZCgu02bp2fbnUU5UIiwt2XTeBUTROb95Vb\nHUe5mJa+GpTCyrNk55XypQVjiA4NsDqOcrGrM0cxPj6U/95WSGenXlnLm2jpq0F59J0jBNtt3Llo\nnNVR1DDw8xO+uSSDg+UNbN5XZnUc5UJa+mrAjlQ0sCG/lC/NTyc2LNDqOGqYrJieTEZCGL99+wgd\nurXvNbT01YA9urVrLP/OxbqV781sfsK3r5xIYeVZNuSVWh1HuYiWvhqQwxUNbMwv5dYF6cToWL7X\nWzY1kYsSw3nknSM6k8dLaOmrAfnVW4cIsdv4qo7l+wQ/P+FfrprI8epGXtldYnUc5QJa+sppO0/W\n8GZBBXcuHq9b+T7k6sxRTEuJ5LdvH9GjdL2Alr5yijGGn20+SHx4IF9ZpKdP9iUiwn3LLqLkzDn+\n9NEJq+OoIdLSV07Zsr+C3JO1fPvKDEID9ehbX7NwQhyXTYrnd1sLqW1stTqOGgItfdWv9o5OHnrj\nIOPiQ7khS8+x46t+uGwyjS3tPLr1iNVR1BBo6at+PZ9bxNGqRn6w9CL89UyaPmtSYjg3zEnjL9tP\ncqK60eo4apD0Hawu6ExTK7988xBz02O4OnOU1XGUxf7lyonYbX489MZBq6OoQdLSVxf06y2HqTvX\nxgMrp+hVsRQJEUHcfel4Xt9Xzt+PVFsdRw2Clr46r/2l9fxl+0m++JkxZCZHWB1HuYmvLh7HmNgQ\nfvzaPlradQqnp9HSV30yxvCT7H1EhQTwnasmWR1HuZEgu40HVk7hWHUjT7x/3Oo4aoC09FWfXt1T\nQs6JWn6wdBKRIXar4yg3c/mkBK6ZMorfbT1CyZlzVsdRA6Clr/5JTWMr/77xADPTolgzW6doqr79\naEUmAD/NLrA4iRoILX31Tx7cUEBDcxsPXT8dPz/deav6lhodwreWTOSt/RVs3qvn3PcUWvrqU7Ye\nrODVPaV84/IJTEoMtzqOcnNfXTSWaSmR/OjVfdTokboewanSF5GlInJIRApF5L4+bg8Ukecdt+8Q\nkXTH8nQROSciexxf/+va+MqVGprbuP+VfUwcFcbXL5tgdRzlAfxtfjy8Zjr1zW08oMM8HqHf0hcR\nG/AYsAzIBG4Ukcxeq90B1BpjJgC/AR7qcdtRY8xMx9ddLsqthsF/bT5IeX0zD10/nQB//SNQOeei\nxAjuuTyD7LxS3irQC6m7O2fe2XOBQmPMMWNMK/AcsKrXOquA/3N8/yKwRPRIHo+yZX8Ff/34FHcu\nGses0dFWx1Ee5u7LxnNRYjj3v7qP02dbrI6jLsCZ0k8Binr8XOxY1uc6xph2oA6Iddw2VkR2i8i7\nIrJoiHnVMKisb+YHL+UzJTmC71w90eo4ygMF+PvxmxtmUneuje+9mI8xek1ddzXcf8OXAaONMbOA\n7wDPisg/HdopIneKSK6I5FZVVQ1zJNVTZ6fhu+vzaGpt55G1swj0t1kdSXmoyUkR/Ouyi9h6sJI/\nfnjC6jjqPJwp/RKg52TtVMeyPtcREX8gEjhtjGkxxpwGMMbsBI4C/7QpaYxZZ4zJMsZkxcfHD/xV\nqEF76oPjvH+kmh+tyGRCQpjVcZSHu3VBOksuSuBnmw+yv7Te6jiqD86Ufg6QISJjRSQAWAtk91on\nG7jV8f1qYKsxxohIvGNHMCIyDsgAjrkmuhqq3BM1/Pz1g1yVOYqb5o62Oo7yAiLCw2tmEBVi556/\n7uJsS7vVkVQv/Za+Y4z+HuBN4ADwgjGmQEQeFJGVjtWeBGJFpJCuYZzuaZ2LgXwR2UPXDt67jDE1\nrn4RauAqG5r5+jO7SI0O5pdrZugZNJXLxIQG8MjaWZw83cR3X9hDZ6eO77sTcbcdLllZWSY3N9fq\nGF6traOTm/+wg/ySM7z6jYVclKhn0FSu98T7x/iPTQf43jWT+MbletzHcBORncaYrP7W08nYPug/\nNx3g4xM1/Pzz07Xw1bC545KxrJqZzC/fOsS2g5VWx1EOWvo+5k8fneCPH57g9oXpXDer98xbpVxH\nRD7ZsPjmc7s5UtFgdSSFlr5PeedABQ9kF3Dl5AT+bXnvg6qVcr3gABvrvjibILuN257OoaK+2epI\nPk9L30fsK6nj3r/uJjM5gkfWzsKmZ89UIyQtJoSnb5tDbVMrtz+dozN6LKal7wOOVzdy29M5RAXb\neerWOYQG+lsdSfmYqSmRPHbzxRyqaODrz+yitb3T6kg+S0vfyxXXNnHzH7bTaQx/umMuCRFBVkdS\nPurySQn81+em8t7hKr751920dWjxW0FL34tV1Ddz8xM7ONvSzp/vmMuEBD0/vrLWDXNG8+MVmbxR\nUM53XsijQ+fwjzj9O99Lldc1c/MT26luaOHPX5nHlORIqyMpBcCXLxlLa0cnP3/9IHab8PDqGbqP\naQRp6XuhU6ebuPnJ7dQ2tvH07XO5WE+VrNzMXZeOp629k19tOUxLWye/uWGmXsNhhGjpe5kjFQ3c\n/MQOWjs6efar85ieGmV1JKX6dO+SDILsNv5z8wHqm9v431tm6ySDEaAfrV5k+7HTrHn8Iwzw/J3z\ntfCV2/vq4nH8YvV0Piis5pYnd+gFWEaAlr6XeGlnMV98cgexoQG8dNcCvai58hhfyErjv2+ezf7S\nelY99gGHyvXI3eGkpe/hOjoND795kO+uz2NOegwv372Q0bEhVsdSakCWTk3k+a/Np6W9k+v/50O2\nHqywOpLX0tL3YNVnW7j1qY95bNtR1s5J4/++PJfIELvVsZQalJlpUWTfs5AxsSHc8X+5/GbLYZ3S\nOQy09D1Uzokalj/6Pjknanjo+mn87PPTsNv016k8W1JkMOvvms/nZqbwyDtHuOWJHVTq+XpcSlvC\nw7R1dPKbLYdZu247wXYbr3x9ITfMGa0XQVFeIyTAn199YQa/WD2d3UW1XPvo+2w7pKdmdhUtfQ9y\nsLye6x6XEAw6AAAKkklEQVT7gEfeOcLKGclk33sJmcl6PnzlfUSEL2SlkX3PJcSEBnD70zl8b30e\ndefarI7m8XRSrAdobuvg8XeP8fttR4gMtvP4F2dzzZREq2MpNewmjgon+55LePSdIzz+3jHeO1LF\nf143jSszR1kdzWPp5RLd3LaDlTywoYCTp5tYMT2Jn66cQmxYoNWxlBpx+cVn+N76fA5VNHDFRQn8\neEUm6XGhVsdyG85eLlFL300drmjgF28c5O0DlYyPD+WnK6dySUac1bGUslRreyd//PA4j7x9hLYO\nw5cvGcvXLx9PRJDOWtPS91BFNU385u3DvLK7hLAAf75xxQS+vHCsnpdEqR4q65t56I1DvLSrmMhg\nO1+7dBy3LUgnJMB3R6y19D3MiepG1r1/jPW5RfiJcNuCdO66dDzRoQFWR1PKbe0rqePXWw6z9WAl\ncWEB3HXpeNbOHU2YD57DR0vfQ+wpOsO6947y+r5y7H5+rM5K5ZtXZJAYqRc7UcpZO0/W8qu3DvHh\n0dOEB/lz09zR3LYwnaTIYKujjRgtfTfW2NLOxvxSnv24iLyiM4QH+fPFz4zhtoXpJIRr2Ss1WHuK\nzvCH94/x+t4y/ERYNi2JtXPSmD8uFj8vP2e/lr6b6ew07C6q5ZXdJby6u5SzLe1kJIRx07zRrMlK\n88k/R5UaLkU1TTz9wQle3FlEfXM7aTHBrJmdxvWzU0mJ8s6tfy19N2CMoaC0ng15pWzML6PkzDkC\n/f1YPj2Jm+aOZvaYaD2SVqlh1NzWwRv7ynk+p4iPjp0G4OLRUVw7LYlrpyWR7EUfAFr6FmlqbefD\nwtNsO1TJ3w5VUXLmHHabsDgjnhUzkrhy8ijCdXqZUiPu5OlGNuSVsmlvOQfK6gGYkRbFFZMSuHRS\nPNNSIj36so1a+iOkpb2D/OI6Pj5ew/Zjp9lxvIbW9k5CAmxcMiGOJZMTuGZKIlEhOgtHKXdxvLqR\nzXvLeKugnPySOoyBqBA7izLiuWRCLFnpMYyLC/Wov8S19IeBMYby+mb2ldSTV3SGj0/UsKfoDK3t\nnQBkJISxeGI8l09KYM7YaAL9bRYnVkr1p6axlfePVPHu4SreO1xNtePqXTGhAWSNiWZOegwz0qLI\nTI5w631vWvpD1NzWwfHqRo5WneVAWT37SurZV1LH6cZWAPwEpqZEMic9hrljY5iTHkOMzqlXyqMZ\nYzha1UjuiRpyTtSSe7KGk6ebPrl9bFwomckRTEmOIDMpgvHxYaREBbvFzCAtfSc0t3VQXHuOkjPn\nKK5t4kR1I0erGimsPEtRbRPd/2tsfkJGQhhTUyKZmhzBtNRIJidF+PTRf0r5isr6ZvaV1rGvpJ6C\n0joKSusprj33ye1Bdj/SY0MZnxDG+LhQxsWHkRodTHJUMKMigkZsP4FLS19ElgKPADbgCWPMz3vd\nHgj8CZgNnAZuMMaccNz2Q+AOoAP4pjHmzQs9lytKv72jk5rGVqrOtlDV0EL12Vaqz7ZQ3dBCWV0z\nxbVNlJw5R/XZ1k/dL9Dfj3HxYYyPD2V8fFjXL9HxfZBdh2qUUl3qmto4UF7PsapGjlWd5ZhjVKCo\npomeF/vy9xMSI4NIjgomNSqYxMgg4sMDiQvr+ooPDyAuLJDIYPuQ9x84W/r9bqqKiA14DLgKKAZy\nRCTbGLO/x2p3ALXGmAkishZ4CLhBRDKBtcAUIBl4W0QmGmM6Bv6SLqyqoYVbnthB9dkWappa6euz\nLNhuIykyiJToYDKTI0iJCiY1OoSU6GBSRvhTWSnluSJD7HxmXCyfGRf7qeUt7R0U1TRRXHuO0jPN\nlJxposQxmrD92GkqGlr6vARkgM2P2LAAstJj+N2Ns4Y1uzPjE3OBQmPMMQAReQ5YBfQs/VXAA47v\nXwR+L10fW6uA54wxLcBxESl0PN5Hron/D+FB/oyJDWF2ejTxYYHEhQcSHxbg+DTt+lQNdeOdMEop\nzxfob2NCQjgTEsL7vL2z03DmXNsnIw+9RyMSwof/tOnOtGAKUNTj52Jg3vnWMca0i0gdEOtYvr3X\nfVMGnfYCguw21n2p379slFLKMn5+QkxoADGhAUwc1fcHw7BnsORZexGRO0UkV0Ryq6qqrI6jlFJe\ny5nSLwHSevyc6ljW5zoi4g9E0rVD15n7YoxZZ4zJMsZkxcfHO59eKaXUgDhT+jlAhoiMFZEAunbM\nZvdaJxu41fH9amCr6ZoWlA2sFZFAERkLZAAfuya6Ukqpgep3TN8xRn8P8CZdUzafMsYUiMiDQK4x\nJht4EvizY0dtDV0fDDjWe4Gunb7twDeGY+aOUkop5/j0wVlKKeUtnJ2n7xY7cpVSSo0MLX2llPIh\nWvpKKeVD3G5MX0SqgJNDeIg4oNpFcazkLa8D9LW4K295Ld7yOmBor2WMMabfOe9uV/pDJSK5zuzM\ncHfe8jpAX4u78pbX4i2vA0bmtejwjlJK+RAtfaWU8iHeWPrrrA7gIt7yOkBfi7vyltfiLa8DRuC1\neN2YvlJKqfPzxi19pZRS5+GVpS8i94rIQREpEJFfWJ1nqETkuyJiRCTO6iyDJSIPO34n+SLyiohE\nWZ1pIERkqYgcEpFCEbnP6jyDJSJpIrJNRPY73h/fsjrTUImITUR2i8hGq7MMhYhEiciLjvfJARGZ\nPxzP43WlLyKX03XFrhnGmCnALy2ONCQikgZcDZyyOssQbQGmGmOmA4eBH1qcx2k9Lhm6DMgEbnRc\nCtQTtQPfNcZkAp8BvuHBr6Xbt4ADVodwgUeAN4wxFwEzGKbX5HWlD9wN/NxxiUaMMZUW5xmq3wDf\nBzx654sx5i1jTLvjx+10XVvBU3xyyVBjTCvQfclQj2OMKTPG7HJ830BXsQzL1exGgoikAsuBJ6zO\nMhQiEgkspuuMxRhjWo0xZ4bjubyx9CcCi0Rkh4i8KyJzrA40WCKyCigxxuRZncXFvgy8bnWIAejr\nkqEeW5TdRCQdmAXssDbJkPyWro2iTquDDNFYoAp42jFU9YSIhA7HE3nklcJF5G0gsY+b7qfrNcXQ\n9afrHOAFERln3HSaUj+v5V/pGtrxCBd6LcaY1xzr3E/XEMMzI5lNfZqIhAEvAd82xtRbnWcwRGQF\nUGmM2Skil1mdZ4j8gYuBe40xO0TkEeA+4EfD8UQexxhz5fluE5G7gZcdJf+xiHTSdT4Lt7z47vle\ni4hMo+vTP09EoGs4ZJeIzDXGlI9gRKdd6PcCICK3ASuAJe76IXweTl3201OIiJ2uwn/GGPOy1XmG\nYCGwUkSuBYKACBH5izHmFotzDUYxUGyM6f6r60W6St/lvHF451XgcgARmQgE4IEnYzLG7DXGJBhj\n0o0x6XT9o7jYXQu/PyKylK4/w1caY5qszjNAzlwy1CNI1xbEk8ABY8yvrc4zFMaYHxpjUh3vj7V0\nXabVEwsfx/u6SEQmORYtoeuKgy7nkVv6/XgKeEpE9gGtwK0etlXprX4PBAJbHH+5bDfG3GVtJOec\n75KhFscarIXAF4G9IrLHsexfjTGbLcykutwLPOPYsDgG3D4cT6JH5CqllA/xxuEdpZRS56Glr5RS\nPkRLXymlfIiWvlJK+RAtfaWU8iFa+kop5UO09JVSyodo6SullA/5/94eoZzbCdjrAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c016278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = T.dvector('x')\n",
    "s = T.sum(1/(1+ T.exp(-x)))\n",
    "gs = T.grad(s, x)\n",
    "dlogistic = theano.function([x], gs)\n",
    "x_eval = np.linspace(-6,6, 1000)\n",
    "y = dlogistic(x_eval)\n",
    "plt.plot(x_eval, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent\n",
    "Some parts of this section were taken from Sebastian Ruder's excellent blog on GD.\n",
    "\n",
    "As we have seen, the key aspect of training a neural network has to do with minimizing a cost function $J(\\theta)$. Gradient descent is a method that allows to do exactly that in a very effective way by updating the parameters $\\theta$ in the opposite direction of the gradient of the objective function $\\nabla_{\\theta}J(\\theta)$ w.r.t. to the parameters. In other words, you try to go down the slope in the direction where the slope is steeper. The learning rate $\\eta$ determines the size of the steps we take to reach a (local) minimum. In other words, we follow the direction of the slope of the surface created by the objective function downhill until we reach a valley.\n",
    "\n",
    "![](images/GD1.png){:height=300 width=300}![](images/GD2.png){:height=300 width=300}\n",
    "\n",
    "\n",
    "The idea is that the direction we should pick to go down the slope to maximize the change in $J(\\theta)$ is that of the gradient $\\nabla_{\\theta}J(\\theta)$. This can be easily proved if we calculate the difference in $J$ for two different epochs: $\\Delta J = J(\\theta(t+1) - J(\\theta(t)))$ (try it!) \n",
    "\n",
    "### Batch gradient descent\n",
    "\n",
    "Vanilla gradient descent, aka batch gradient descent, computes the gradient of the cost function w.r.t. to the parameters $\\theta$ for the entire training dataset:\n",
    "\n",
    "$$\\theta = \\theta−\\eta\\nabla_{\\theta}J(\\theta;x;y)$$\n",
    "\n",
    "Keep in mind that the cost function $J$ is the product of the individual likelihoods for each example in the training set, and so in order to produce one single update to the parameters, you need to use the entire training set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classifying the Iris data set\n",
    "\n",
    "The Iris flower data set is a multivariate data set introduced by Ronald Fisher in his 1936 to describ morphological properties of three related species of the Iris flower. The data set consists of 50 samples from each of three species of Iris (Iris setosa, Iris virginica and Iris versicolor). Four features were measured from each sample: the length and the width of the sepals and petals, in centimetres.\n",
    "\n",
    "The goal of this part of the lab is to perform classification of this dataset using the logistic regression, and batch gradient descent to perform the minimization of the cost function (training)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/python3/lib/python3.6/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['logistic']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n",
      "  \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import six.moves.cPickle as pickle\n",
    "#import gzip\n",
    "#import os\n",
    "#import sys\n",
    "import timeit\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import theano\n",
    "import theano.tensor as T\n",
    "#import matplotlib.pyplot as plt\n",
    "#import csv\n",
    "#import pandas as pd\n",
    "import six.moves.cPickle as pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#iris_dataset = pd.read_csv('iris.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#feats = []\n",
    "#species = []\n",
    "#for i in arange(len(iris_dataset['Sepal.Length'])):\n",
    "#    feats.append(np.array([iris_dataset['Sepal.Length'][i],iris_dataset['Sepal.Width'][i],\\\n",
    "#                 iris_dataset['Petal.Length'][i],iris_dataset['Petal.Width'][i]]))\n",
    "#    species.append(iris_dataset['Species'][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#iris_data = {'features': feats, 'species': species}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# pickle_file = open(\"iris_dataset.pickle\", \"wb\")\n",
    "# pickle.dump(iris_data, pickle_file)\n",
    "# pickle_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the dataset\n",
    "\n",
    "The dataset is saved as a pickle file, and you can easily separate it into features and classes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pkl_file = open('./iris_dataset.pickle', 'rb')\n",
    "iris_data = pickle.load(pkl_file)\n",
    "features = iris_data['features']\n",
    "features = np.array(features)\n",
    "species = iris_data['species']\n",
    "species = np.array(species)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us explore the dataset a bit by plotting the features color-coded according to their class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2]\n"
     ]
    }
   ],
   "source": [
    "uni = np.unique(species)\n",
    "\n",
    "species[species == uni[0]] = 0\n",
    "species[species == uni[1]] = 1\n",
    "species[species == uni[2]] = 2\n",
    "species = species.astype(np.int)\n",
    "print(species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Shuffle entries and separate into training, validation and test sets:\n",
    "\n",
    "indices = np.arange(150)\n",
    "np.random.shuffle(indices)\n",
    "\n",
    "species = species[indices]\n",
    "features = features[indices]\n",
    "\n",
    "train_set = [features[0:110],species[0:110]] #110\n",
    "valid_set = [features[110:130],species[110:130]] #20\n",
    "test_set = [features[130:150],species[130:150]] #20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x10c344ef0>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+QHGd557/PzK6MVkIrtFLlbImZMTHlKwVxGG+RgFOc\nj10SkI2hqMsVZOwokKtJdiFlY2IO2KrY8t1Sdxxl5D9uNzVlLGx2AkWMncRYJOAFJ7HqYrIyBmGZ\nJGC0i4TB+hGvLK3O2p197o+e3p3ped/u953pnu7peT5VXbvT887bz8y8/cz7Pu/zg5gZgiAIQveQ\niVsAQRAEwQ5R3IIgCF2GKG5BEIQuQxS3IAhClyGKWxAEocsQxS0IgtBliOIWBEHoMkRxC4IgdBmi\nuAVBELqMvig63b59OxcKhSi6FgQcOXLkNDPv6PR1ZVwLUWIzriNR3IVCAXNzc1F0LQggovk4rivj\nWogSm3EtphJBEIQuQxR3G1QqQKEAZDLO30olbokEQegFRHG3SKUClErA/DzA7PwtlUR5C0LSqRyt\noHCggMz+DAoHCqgcrTQ9R/sJfXf3gfZTU5skEImNuxeYmACWlhrPLS0554vFeGQSBMGfytEKSo+W\nsLTs3Lzzi/MoPVpae77+uSpXm9oU9yTj5hbF3SILC3bnBUGIn4nZiTXF7LK0vISJ2Ym1/1W4bURx\ndzm5nGMeUZ0XBCGZLCyqZ1a687ZtOoXYuFtkchIYGGg8NzDgnBcEIXr8bNVexh8bR9/dfWCoK37l\nBnPIDfrPuoKe7yQy424R1449MeGYR3I5R2mLfVsQosfPVu01Z4w/No7puWltXwP9A5gccWZc9X3q\n2iQBmXH7EOTuVywCx48Dq6vOX1HagtAZgmzV9ZSPlLX95AfzKL+njOKeIop7iii/p4z8YB4AkKVs\nU5ukIDNuDa67n+s54rr7AaKgBSFudPbm+cV5bP/sdpy5eAZZyq55hvhxeOEwJmYnsLC4gNxgbm1m\n7Z5rh8rRSlPfYfwAUBRV3oeHh7nbQ4MLBfXmYz7vzK6F+CCiI8w83OnrpmFcpwVXOUdBf6YfRIRL\n1Utr5wb6B6xn3V5zTlA/NuNaTCUaxN1PEHqT5dXlBqUN6M0wftiYc2zpOcVtGqauc+sTdz9BiJ+z\nF892/Jq2ZpN2XA+D6CnFbROmLu5+gpBcWnHNczcbo7xmvYtihtTqNQy3wp5S3H5h6l6KRaBcdmza\nRM7fclk2JgUhCUyOTGKgfyC4YR0mG5WAY+PekN3QcM7EHdC1ac8vzoPByuuF5VZorLiJKEtE3yOi\nr7d91ZiwtVvbuvtJtkBB6Aw6172hjUMY2jjUcM5mpp2lLA6+7yDuf+/9yA/mQSBjd0CVTdvt06Yf\nE2zcAW8F8ByALW1fNSaiDFMX90FB6Cyu77UK1w1vfnHeeKYNrM/K/fr2u5aKVV7F6p2rxjKYYDTj\nJqJdAG4AcF+oV+8wUdqtbcwwgiBER73JohVKj5aM07iaXCuKUHlTU8kBAJ8AoP3ZIKISEc0R0dyp\nU6dCES5sorRbi/ugICQDncnCFBuXvaBrRRUqH6i4iehGAC8y8xG/dsxcZuZhZh7esaPjdVyNCbJb\nj48DfX2OYu/rcx6bIO6D3QcRvZaIvkNEx4joWSK6NW6ZhPYJw93OtA+/dlGGypvYuK8DcBMR7QXw\nKgBbiGiGmW8OXZqYGR8Hputy0VSr64+npvxfOznZaOMGxH2wC1gB8HFmfpqIXg3gCBF9i5mPxS2Y\n0Eh96Pi2jdsAoCGs3f2bH8xj28ZtbUdV+pk36mXJUMbKhh4WgTNuZv4UM+9i5gKADwD4dhqVNuCY\nTWzO1yPug90HM7/AzE/X/n8Zzub7znilErx43ezOXDyzpphdpVlfreblSy+jP9Pf8vX8zBsmLn8u\nbsbCKMqepd6Pe3TUUaTuMTqqb1vVfAe68173P0CyBeoIyp1sk1s5CoioAOAaAE919MJCILY260vV\nS9hy2ZYmV0Gvy6Drojc2PGbs+ufn8qcirBB3L1bZAZn5CQBPhC5FRIyOArOzjedmZ53zjz/e3D6b\nVSvprOI7Efc/c4JyJ9vkVo4CItoM4GsAbmPmc57nSgBKAJCTDYtYaMVmffbiWZz+xOmOybLKqyCQ\nslBDFJVzUj3j9irtoPOlkvl5cf8zJyjZTpTJeIIgon44SrvCzA97n++WTfc004o7HYPbWrnpVoA6\nWfwq6MTpDtgTTE0BY2PrM+xs1nms2pgU9z9zgpLtRJmMxw8iIgBfAPAcM98T6cWElmklvB1o3cbs\ntWPX96OSxbWJ+z0XNj2vuFV26l27HHv4rl3AddepXyfuf+YEzUQ6OVPxcB2AWwC8g4ieqR17o76o\nYIc3vN2GsNOx1svitYn7PRc2qa6As3s3cEzh2LV7t/NXZaeudwf0s1uL+585kyOTyoTy7kwk6Pmo\nYOYnAVCkFxHaYvyeUZQXZ1F1v6UWvq35xXnQ/vUX5gfzTVVu6ivf6KIg5xfnUThQWGv/pfd/qUkp\n24TKt0OqK+AEVbHRPa9r76VSkWLBpgSVcLIp8SQVcHqD8XtGMX1uNpKfVlWVG9W5erybj61UxfHD\nZlynWnFnMk7ebS9Ejsue7nldeyEZiOLuDfruJFQTYszVeYzkB/M4ftvxcK6RhtJlNilSvW3Hx52/\nOqXs2qFN7dG9areO27da6G2qLc60XRtzmKiUNhD9BrqORCpum0o1qrbT03oTSL0dWpUt0K99L+G3\nsy4IYaKbIGRbMAa4M+DVO1fbrnhT36duY7QDG+hKEqm4bXykVW11eMPQVWHqY2MStg7E61st9A5+\nE4TS4Ag0E12AgawnWM67oV26VhOYAXWVG7/KN5109TMhkYrbxkfaxm/6xAng4EH/MPWpqfXHk5PO\nD0MvVrSJy7da6C38JghTtz+OsS0jyK7CUeC1I1sFxr4LPPCXQP4lgBjIn882bRRO3TCFseGxppl3\nfjCvrHLjV/mmk65+JiTSHdCmUo2urYpqtTFq0s/dr9dD2nODOaVbVFxLQyGd6CcI80ChgKmFBUzl\n8usuWx6PguJR978q8CYAexr7mbphClM3qFN7qsx+fu58nXL1MyGRM26bSjUmdmo/bEwwvRTSnrSl\noZBOcn3b1OcXod7k8vMU0G2EKej2PZxEKm6bFKnetq1gY4LplZD2pC0NhXQy+Tgw4HGbHrjknG/A\nnTX5zdQsZlbdvoeTSMUdhLdKzeHD63ZpVSa/IOp/xN2+g1wJrQnwb6zfWd/+2e3Y/tntsbvhFfcU\n13boj992XJS2EDrFvzuL8qN1tuqXgPKj9SaQOhYW1mdqOgxnVt2+h5NIG7efffnwYf8qNVdfrQ5z\n92NvLTuFtwKOl5ZdAwMM5t60pvXVOzqd4lQQIkEXZpzLoXh0Xq2ovTADRKi8kTDxMcL8FkZ2Fahm\ngPwiMDkLFM+Zzay6fQ8nkZGTfqHqJ07oc2avrDizZV3hAx1uSLvfa/P5NkLaA2LvCwcKgRWpw4zQ\n6nYkcrLL8E5cAGcWVC4DH/0o8NJL5l3tAUrvAZY2ND83sAyUd46hOBZQZxDNOeKB8EPYben6yEk/\n+3JQlRpbpV1/Pb/XtlXRJsBgbrI8M13CjT82jr67+0D7CX1392H8McNqx4IQFX47/RZKGwAmRtRK\nGwCW+oGJVw4Z9dPteziJNJX4uQP6zbgBx+5tu4hw7dY2FXCsL+Dj36hbttWzLbMp8DLjj41jem7d\n1lPl6tpjnUuUIEROiDv9C4MBz1vYqJPk3mdLImfcfu6AQVVqNgXrN2W/9X3o+m6ZAP9Go0TxSxcC\nL1M+ot600Z0XhFDRbcDrdvRb8OPNLQY8r7FRpy3vTiIVt587YFCVmgvB+g2A2s3QpgJOaG8InkTx\nmtXC2cuClxG6itN+lagFIRT8EgypJi6ZjPnNWsfkbLP7oIsuzqDbfbZVxKa4bbL/nT4N/N7vrbv/\nAc5GJLPzt16xmrjr5fP6SuxTU/q+rbAsAe+63uUvqO0ymVUgcxehcEcfxu8ZVSfl8Umq47Ydf2y8\n8bXT49ovIm2zFCFC/OzYqolLixSPOu6C+epmAOtj3s9G3e0+2ypi8Srx22QuFtXPe9HNhE1c+iJP\nHBX0Bv1eOj2O0slpLPXXnWQ0JpP3PHZ3ww8vHG6wcZswsAyU/7rOb7YmZ+WNSNyuu4t4lcSIzq3P\nNLk90NpGlPf1FgnyM/szyrSsBMLqnclJtJ/4QgphVKZx3f9M+3b770iVmqA3GEBlehwTz5exsKmK\nzCpQNdgcdd0Fxx8bR/lI2co8kn8JOH6gUc7CbeoSTklwSxTFHRN+E5KJCfOkQe1ieB+56NxtkzCW\n60m84g6rMo2qTVDfHSFEITJ3EdgglF81e9DNNJpey8Dq/kY5M3eqk8cnYZYiijsm/CYkqiKsUdDC\nkjmJPtsqEu/HHVQh3TSs3DXLjo6uh8CHHqreCiGWgM9pbN5N7Wq76ZXbR1H4GCFzFyFTNftRbtqp\nz+XirLwuJBU/t756O3ZUtJggv9t9tlXEoriDsv+5IehBuJvXs7P+wTMdr2Jjk94wgL2XX9/saeJ5\n7O6mV24fRWnjLOa3Akw1E0vAa8HAVaeb5ZTsgEITIU5IrMjngZmZtqLg0pZ3JxbFHZT975BZ8FMg\nsVWxsUlvGMChzI+bq1yTs5vunT1M8GxzVBk5iecJ1LzJWXv+iSvRJGcaZylCm/hNSOrdAcPGr3Zh\njxKouInoVUT0XSL6PhE9S0T7g16jwusdd/iwvq1uRUbUeupWY/z8FL1pCUdH9a505w6j8DsnkPlT\nRuF3TmD8lwcbXOvG7xlF4Y6+NRe/yrQ6NF0XCVbl6ppP6sHHPgMUCtqosqBK2dUM1twUK2/Empy3\nfuNWnL903v/FGryuhE1uiOJa2H34TUhsagi2Qi8lwzcgcHOSiAjAJmY+T0T9AJ4EcCsz/6PuNd5N\nHBP3vvo9B789kNOn7f32jfcz/HbNvWkJfS5UOXfY3qVPkyDHJAEVGBj5CfDM5cAZy8hRAMiuAiv7\nWbmJU4/phk5QPzZ9qZDNyQRi4w7ohzszi93DoPOEujnJDu60q792WH1DJj/G9T+ofiuyixdtrtzc\nt7Wg7ov9cgB72k48X25U2oDSRNHw0n5g4vnmaxiFwxMw+6vB4ilhOEVZoQ5UaJDRMGghqB+bvoQI\nsYmCC+onE5LVNZeLz5beRRh92kSUJaJnALwI4FvM/JSiTYmI5oho7tSpUw3PmeaScdv5rcha/cE1\nkqGVtISKtgubWgsxV73Oa2v246xp6ge36OoqMLZlBFO3O+VGwspSaJrop1uS1qcSvxD1VvqxTcu5\nYQPQ75nduLOzEDf304qR4mbmKjO/CcAuAG8hojco2pSZeZiZh3fs2NHwnOkPZX27YlEdId5qpj4j\nGfx+6U0vnMsZu/A1vXQRQCaDyn/ajsLkegUcPHkYxw/UfK191jpBCXhc8lvz4LsYK/t5TWkDZq5+\nGcoE2qlNXQbjdi0kovuJ6EUi+mGsgsRBWEVVTZbT+bwT6lw/E7v/fuDgQfXsLMTN/bRitb5h5pcA\nfAfAu2xeZ1rQ18QNsJVMfcY/1n6/9Ndfb3axvXsx+boSBpY954Nc+i45CXQqb2CU3nYG8ytn1hPi\nnJxGZYszMxr5ibqvkZ9oEvBoXAdVmJhl6jdFdYl6TPpJiGvhF2E5llNDUKpV70b8r/1a4+PxcWe2\nbeJFcuLE+v7Ql760PhPTzc4A/+cEI6+SHUS0tfb/RgDvBPAjm4t4f0B1k1cTN0BVBr+RkcYfZ++P\nu/GPtd8v/Y9/bPZmDx1CcWwK5Z1jyJ/PrtXRG/tuY10972O3zp4qUfxSv3MeAB6fQaPyrintx2fq\nEvC4/Z7PYmzLiLFLn9csM7RxCEMbh0AgZQIrnZ1a5Uo4NjyWONdCZv57AGdjFSIu/FaXbsKf+uok\nx441Pp6edjK/meC+Ttz6QsPEq+SNAB4AkIWj6L/KzHf7vabdkPdEksmg8gbGxIiTzD3n1rirq5VX\n2eMo2IWthNxgDpMjk45yMthxd187P4jmjUxgTVFnV4HSHDD1DfiHirofpi4pkO3b9wmfJ3jer+r9\nHa1gYnYCC4sLgW2DCNOrhIgKAL7OzE3mPy+p8irx86Dat6+1UlKmWOYa6RVsxnVgBRxm/gGAa9qW\nqo6AgjCJpPIft6H0tjNrs+H5rU7tO8BR3o218LixyK/uDbt9+9TRW6OmzKtZYPotzv9T3/AzeOcC\nixTb4Felh73v16OQve6B3VQAmYhKAEoAkEvyALXF/f5VP+o33xzttVuofCM0ksiQ9yQyMaowYWxY\nN2EoTRyuKSHAyO9XR08JAWW/32X3wwxrAwpmdmud6aSb8yH7bbp3PWF7AADObDooX0mafgBjIpEh\n70lkYUVtCnWjFXVRiwuLC81vWNOHDcpoSO+HGWKtP1O3RJWLn87tT9wBE4qpB4DXd9vPnc/bRmiL\n2CrgxLZp3GLQgTZb3mLjX9PXqfqwIevdC1CV9Qk5kKE+UU9+UD2rUr3fpGYaJKIvA/i/AK4mohNE\n9AexCpQUVB4Au3c31/R78MFgdz63PdAdM7QuIZE1JyOjjaCDycv2Nrn4uS58QC3Dnirz3ovLzdf1\n9u1TR08JA9f/tF4QzSwmQpuUTfbApGYaZOYPMvPlzNzPzLuY+QuxCpQkvDX8nn22uaafiTuf255Z\n3PpCpLcUdxs23+L/OoTyX6td+AA3w57nRQQ80ffzwCCF4rEsyhdHtLPYJgj48b/LBtuZIrRJ2WQP\nlEyDghAusVTAiY12/BADXkt3kdaNj+/WuO0prmtctSYBlWjiQpJMCWkk8RVwYsPW5ltvD9cl0ak9\n12RzrpFdVfdf2QMUPt4cPt4t4eK2SMV4QQiP3lLcNjZfr11aF5BQrQLMKM1BaeMu/eKKputW9gCl\nm4D5zc3h410ULm6M68c9vzgfGC4vCEIwvaW4bWy+Ort0NquM25/6hhPGnq3CybxXdR5PfbO/6boT\nv51tSvvq+jV3S7i4Dd3sxy0ISSQVitt3Ge51/wPM/BB1vs51+Roqe4DCbUDmTufvdSeAlf8O8H7n\n79Q30Jirtnbdhc1qu8qaX/OTh4GTJ5yZ/skTuG4BxvXykmiSED9uQQiXwJD3pOMbTv0DtB7yvW0b\ncOaM+jnmpjB1bwj8Ggr7ti58PDeYQ2V63Kmes7nW7+YqSiengWk0VcfxktTQcr/3KwiCPV0/4/Zd\nhocY8t10XVWIe10IPACt/dzPr1lVPUdXHadJpoSaJJLqxx0nYRWfEXqTrlfcvsvwdkK+z/pn+9SG\nuA8i0H7u59esq55jUlUnqSYJ8eNuJKziM0Lv0vWmEt9leA6tpyEMyOiXW3TMI03nt+aB1eOB3Rf3\nFJWKK3chi/nNzUrapKpOkk0Suvfbi/gtBCWwUDCh62fcylD0Zed8WyHfV13lf91ZNF83hOW/qnrO\nwLJzPvC1YpLoCnQLvvl5MZsIZnS94laGov+1c76tkO8nntA/R4TiubxT5Sbk5X9T9ZzzWZR3jgVu\nTAJikugW/BZ8YjYRjGDm0I9rr72WbZmZYc7nmYmcvzMzhi8kYnZMhc2HdWd16PoE7PuypeUPozcA\nMMcRjNugo5VxrWJmhnlgwH+I5fOhXEroImzGdSJm3G1t1vhNX9rZ+VHkzfY9Hxayc5V6vFlPVUiR\nGMGPRCjutrz2TErIt+ICuGmT3fmwiNCFUegsKpe/0VHnt//mm/0LpLvpccTmLahIhFdJW4VavLXz\ndNkObacwFy7YnQ+LEKvWCPGhKvdpU8rRWxgdEI8TYZ1EzLjbLtRSn9Bdt/60rfoScvWYyvQ4Cnf0\nIXMXoXBHHyrT454GtemZ7odH6vR1FQEp2K2QBZfgJRGKO9RCLWF1FqJQbhj7/OYqmNbD2NeUd71d\nW4XU6es6wl4gyYJLqCcRijvUQi1hdRaiUIFh7H7TM6nT15UEbbvYIgsuoZ7eqoATE5m7CKxwRiEG\nVu/i9irz9CDdUAEnmw3vqxsYkN/uXkAq4CQMXbj62vmQ7elC/IShtEMuEyqkCFHcHSAwjD3CauxC\n9Hjd/kZH2+8zG5yaRuhhAhU3Eb2WiL5DRMeI6FkiurUTgqWJwDD2CKuxC9GiipeanW2/31pFPIm/\nEpQE2riJ6HIAlzPz00T0agBHALyPmY/pXhO7jbtSWffrzuWcmasowdSQJBt3oeAfSBMW+bzj8Sqk\nF5txHRiAw8wvAHih9v/LRPQcgJ0AtIo7VlSRDxLBIEREp9z0xB1QqMfKxk1EBQDXAHgqCmFCQULG\nhQ7Sqf1j2acW6jFW3ES0GcDXANzGzOcUz5eIaI6I5k6dOhWmjHZIyLjQQUxS5bSL7FMLXowUNxH1\nw1HaFWZ+WNWGmcvMPMzMwzt27AhTRjvEtU7oIO6+crteIETr+9JjY7JPLfgTaOMmIgLwBQDPMfM9\n0YvUJpOTjTZuQKYsghYieheAewFkAdzHzP/Tto9iEbjlFl3/+vQz9UiclWCDyYz7OgC3AHgHET1T\nO/ZGLFfriGudYAgRZQH8HwDvBrAbwAeJaHcrffkt9IJm4+KzLdgSqLiZ+UlmJmZ+IzO/qXYc6oRw\nLVOfLfD4cVHago63APgxMz/PzJcAfAXAe1vpyC+GqhRQLjToeUHwIpGTQi+zE8DP6h6fqJ1bw3TT\n3W+hNzXl2K29M+ts1jk/FVxOVBAaSEQhBUFIKsxcBlAGnAAcv7bFon5xNzUlCloIj0gU95EjR04T\nUQfiyQLZDuB03EJ4SKJMQDLl0snkU63RipMAXlv3eFftnJKYx3USvx8TRG5zjMd1JGldkwIRzcUR\nGu1HEmUCkilX1DIRUR+AfwEwAkdh/xOA32XmZ6O6Zqsk8fsxQeSOBjGVCD0LM68Q0UcB/C0cd8D7\nk6i0BcGLKG6hp6l5SCXbS0oQPKTdq6QctwAKkigTkEy5kihTXHTrZyFyR0CqbdyCIAhpJO0zbkEQ\nhNQhilsQBKHLSKXiJqLjRHS0llclMeXmiWgrET1ERD8ioueI6K0xy3N1Xf6ZZ4joHBHdFqdMLkT0\nsVqpvB8S0ZeJ6FVxy9Rpur1sIBFlieh7RPT1uGUxJWn3qI5U2riJ6DiAYWZOlOM/ET0A4B+Y+T4i\n2gBggJlfilsuYC3h0kkAv87MsQZPEdFOAE8C2M3MF4noqwAOMfMX45Sr07RSNjBJENHtAIYBbGHm\nG+OWx4Qk36P1pHLGnUSIaBDA2+GkyAUzX0rYgBgB8JO4lXYdfQA21oJkBgD8PGZ5Og4zv8DMT9f+\nfxmAWzYw8RDRLgA3ALgvbllM6YJ7dI20Km4G8E0iOkJEScm9diWAUwAO1paP9xHRpriFquMDAL4c\ntxAAwMwnAXwOwAKceqeLzPzNeKWKl64oG9jIAQCfANBNmcaTfo+ukVbF/ZvM/GY4eZY/QkRvj1sg\nODPINwOYZuZrAFwA8Ml4RXKoLQlvAvAXccsCAET0GjjpVa8EcAWATUR0c7xSxUdQ2cCkQUQ3AniR\nmY/ELYslib1HvaRScddmbGDmFwE8AifvctycAHCCmd0Z00NwBkkSeDeAp5n5l3ELUmMUwE+Z+RQz\nLwN4GMDbYpYpFkzKBiaQ6wDcVNtr+gqcIiwz8YpkRJLv0QZSp7iJaFNtIwe1Zc5vAfhhvFIBzPwL\nAD8joqtrp0YAJGWT6YNIiJmkxgKA3yCigVrpvBE49t2eouvKBtZg5k8x8y5mLsAxwX2bmRO/Ykr4\nPdpAGnOV/AqAR5wxjz4Af87MfxOvSGv8MYBKzTTxPIAPxSyP++P2TgB/GLcsLsz8FBE9BOBpACsA\nvoeEhyBHhFs28CgRPVM79+nEV6DqbhJ3j6pIpTugIAhCmkmdqUQQBCHtiOIWBEHoMkRxC4IgdBmR\nbE5u376dC4VCFF0LAo4cOXKamXd0+royroUosRnXkSjuQqGAubnE5HYSUkZcBXtlXAtRYjOuxVQi\nCILQZYjiFgRB6DJEcQsdoXK0gsKBAjL7MygcKKBytBK3SEIMjD82jr67+0D7CX1392H8sXGr18s4\nckhj5KSQMCpHKyg9WsLS8hIAYH5xHqVHnaSNxT3FOEUTOsj4Y+OYnptee1zl6trjqRumAl8v42gd\nmXELkTMxO7F2s7ksLS9hYnYiJomEOCgfUWct0J33IuNoHVHcQuTML6o3y3Xnhc7RrunB5vVVrmrP\nm/SxsLhgdT7NiOIWIqH+hhaSiWt6mF+cB4PXTA+mytv29VnKavsy6SM3mFO+Vnc+zchdJYSO94YW\nkkm7pgfb15euNStGpetjcmQSA/0DDecG+gcwOTJp1G+aEMUthI7qhlbhNwMToqdd04Pt66dumMLY\n8Nja9+73/av6KO4povyeMvKDeRAI+cE8yu8p99zGJCCKW4gA0xu/fgYmbl7hYPM52poevK583tlv\n/etHHxwF7ae1Y/TBUQCO8l750xXwnYyVP11BfjBvJUNxTxHHbzuO1TtXcfy24z2ptAFR3IIFpkoh\nyOaYpSzGhsfWXMDatbUKDrafo43pwXXlczcYq1zFheULyn6Xq8uY/elsw7nZn86uKe969r5+r7IP\n3XnBQRS3YISNUtAphJn3z6zNtOr9dsXNKxxsP0cb04Opyx4A/Pz8z5XnvcocAA79q7qYj+684CCK\nWzDCRinY2iLFzSscovwcda587RL1d59WE5wobkFL/aD388VW3RReWyQA7Q0kbl7hYPs52qyiMhSN\nqojyu0+zCc7o2yCirUT0EBH9iIieI6K3Ri2YEC82Ln1BN0XQDSRuXuFg+znarKI29m00kmGgfwC7\nt+9WPjdy5UjbMtuQZhOc6c/ovQD+hpn/PYD/AOC56EQSkrC8M3Xpq0d3UwTdQOLmFQ5Rmqj8xoL3\nes9+5NkmJT1y5Qge/73H25bZhjSb4AKTTBHRIIC3A/h9AGDmSwAuRStW7xJ2Ip3K0QomZiewsLiA\n3GAOkyOTvv247f3C0QmknYWrbgqTG6i4pyiKOgRsPsfcYE75PavMFLq2QxuHlH2rlDSgH49RfPc2\n76/bMJkhi5QGAAAYmElEQVRxXwngFICDRPQ9IrqPiDZFLFfPEubyztbGV99eR5ayWL1z1cr/VmzY\nycTGFW9Tv/qWP3PxTEvjqxM25zSb4EwUdx+ANwOYZuZrAFwA8ElvIyIqEdEcEc2dOnUqZDF7B5PZ\nqZ8ppf65fY/ss/oRMDGPuN4FNjdFmm+gbsbGFe/Y6WNGfdqOryhtzmk2wZnk4z4B4AQzP1V7/BAU\nipuZywDKADA8PCwJKlokaHnnZ0oBgH2P7GsIklDheoJ4TScmtj93aey+xsQMY9O2kxDRawE8COBX\nADCAMjPfG6tQHSQqG7Btv1HanNNqggtU3Mz8CyL6GRFdzcz/DGAEgNnPr2DN5Mhkg2IGGmenfrOW\n00unjf1t65eqgDPAdT8aOmxuioTeQCsAPs7MTxPRqwEcIaJvMXNXjG/b/QsvUdmAc4M5pWxx2Jzb\n/YySiqlXyR8DqBDRDwC8CcBnohOpN9CZO4p7injrrkZvy62XbcW+R/aB9pNWsS4sLmhDkP2oX6qq\nTBpezl48a32NpMLMLzDz07X/X4bjLbUzXqnMCMNefNW2q4zP61z8vAz0D2Dv6/cqZdv7+r0dNZn1\nvB83Mz/DzMPM/EZmfh8z/1vUgqUZvwE1/th4U2jwz8//PHAm3c6sxV2q1tsEo7hOkiGiAoBrADzl\n3zIZhGEvfuL4E8bndZOCoY1DTTbkQ/96SCnboX891FGbc5r9uKXmZAz4DagT50601Ofe1+9tqOdn\nQ70ydk0aXls6kN4NRSLaDOBrAG5j5nOe50oASgCQyyXnR8vWXjz64GjDhGDkyhHfijSm/Z69eBan\nP3G64dwtD9+ila2TJrM0+3FLyHsM6AbO/OJ8yzkhWk3KsyG7QamM07wjXw8R9cNR2hVmftj7PDOX\na6vN4R07dnReQA3bNm4zPu9V2oA64ZOLKrzdxqUzKe6fSZEjCkRxx0AUA2dhccHXxKGDWe8AlPbc\nx0REAL4A4DlmvidueaLCT0mrUIW3d6P7Z1LkiAJR3DFgsgloi7tjbtvv8upyKmx+LXIdgFsAvIOI\nnqkdXZEIWrdJHMbmscqX32YFlpTVWlLkiAKxcceA16+53bqM/Zn+NTenwwuHUT5SRpWrIBCymSxW\nVld8X58Gm18rMPOTAChuOVohSte6gf4B9N3dhypXkaUsSteWMHXDFD7z959Zu+b84jw+8/eOc5nS\n3e7Jw8DJE8Amdv4+eRiIQWEm1A21bWTGHSLeck2v+h+vWnf3mx4HCgUgkwEKBRR/gDUzRCsmjnpc\nxV85WsED339gzU7OYBAIQxuHHCWuqfGXBptfr2ETrq5z5dt62Vbl+QvLFxqCuKbnpvGa//mapujJ\nY6eP4eaHb27yjhq/ZxSlk9OY31wFEzC/uYrSyWnnHhBCQRR3SKg2gF6pvrI+oE9Oo7JlHmAG5ueB\nUgmo6NOa2rCyuoJbv3Gr0ltleXUZmzdsxuqdq9oq21ImqvuwCVfXufK9fOll4+u99MpLRu2WlpdQ\nXpzFUr/nfD8w8bx5FR3BH1HcbVAfRBO0AbTUD+x7H5C5EyjcBlR+dQmYaE5r2ipnLp4JdH+SMlHp\nwcbVTdc2qqo2VY3xaWFTNNfrRURxW1CvqLd/djs+/FcfNio04FLNwlk6bgU+/F5g+39Zrx4DYK1S\nTKvk+tQuYu75NPu19ho27oA6U5jOdNYuWc3tkLsQzfV6EVHchnijHc9cPINL1dbTkl/qA85saq4e\nk9HslWVA4DtZm/94aOMQJh8HBjwiDVwCJmupkYMUu5BOdG5xpWtLTeczGpWw9SLQND9hAKuNpwb6\nB1AaHMHAsuf8MjD5OrWpTrBHFLchrVSEsWFpeQn7HtmHVY1ftXv+3nffi/5MowGxP9OPe999L4p/\ndxblR4H8SwCx87f8KFB84gyQyWDyq2d9FbvQPdi4A+rc4qZumGo6/+D7H8TY8NjabDxLWYwNj+Hf\nPkvY/Us4yrp27P4lMPNIcwWcqdsfR3nnGPLns844PJ9FeecYimNT0X0gPQb5BWC0yvDwMM/NzYXe\nb5xk9mfs3PYYuGzFmVkTA6vt/kQywHete48oXbAKBWfj04fKHmBiBFgYBHKLwOQsUDwKgAjI5YDJ\nSaCYbPcpIjrCzMOdvm6SxnXhQEHpDpgfzKtNbpWKs6eysNDa91wooLJlvnnsnMsDxxXX06GSA2hP\ntpRgM67Fj9sQk5SnA/0DjoP/eyYaFOj2OxyzSDsM/b91E4rWN3XvXmDaP19J8WhNUXup93YBevLG\n6SaC0v82UKk43+tSrW0L33Plv+1F6eT0mrfI/FagdBOAnXthPFJUcnz4w87YW15uWbZeREwlhqjs\nhP2Z/jUf6YaorIXGzb6zbQZJblgB7v27uk4qlQafcIzXfMQDlLYRS+veLkJysYoKnJhYV5Yult/z\nxCuH1C5+r1h4JKnkuHRpXWm3KFsvkgrFHWpVdK9SrKznyS6/Zl+D3e7gv70dp6c2Y3U/cPwAUPxB\nrQ9PFrncovpS2apjRsnqvKRqdur7/woo/uMFR6bt24EPfciZmbiz5OnpQBNJ44WzjmlEx4J4mXQD\nxrlkdN+nxfccikeSzbiSMehL1yvuUJOlu0u5eUWgTKWC4p88gOOfqzqK+nNVFD8/q247OQkMrM+Q\nJ2fV3h4P/CWwut/5q3p+5uHaD4Jr2mAGzpxpnqHYMDAAPPAAsLoK5DV+4wlKX5pWKtPjKNzRh8xd\nhMIdfa1FFWomGU3ovs9cTi2Hol9tpr3qJqCvz5kI9PU5qz+dbDbjSsagP8wc+nHttddyp8h/Ps+4\nC01H/vP5FjrLMzvqsfnIZvXP1R/52nVnZpz/iZgBntkDzt8GpjudvzN7Gl8X9HzLR18f89CQI0c+\n78jlMjPDPDDQ2H5goLFNAgEwxxGM26AjrHE9MzXGAxON43VgAjwzNWbRicV3p2k787GRZjk+DZ55\nU7a57W9fwQOfVrRVjdOREbVsY2PN5zdsYO7v77oxGAU247rrvUp03h4Ewuqdq4pX+HWWcYZOOxA5\ns1nvubjo7wcOHtRv9LTrbRAD3e5VUrijD/Obm+1j+fNZHP/f/gnB1jspqM1jeY2Xh+J7LjyzTy3H\nS85Kr6kLnUeSKfm8M77Eq0RJT3mVhJolLZezsxXr+gAab5Q4WV525NDdCMViT94kcaIL/bYKCdeM\nq8qWeUwcKBgVx9XKMai+pNYjyZSFBf14kzFoRdfbuENNlu6xTbfEVVc128p1bNrkzMYzLXwNROYz\n+bh/PIQGdKHfViHhChtwZQ9Qei817/dMjyv3bnLn1ONHt5neNtskQjcsul5xh5osvVgEymX9pp0J\ns7PAzTc3uz3Vk80CY2PA+fOOWaVV84zfBmM9stGTKCZfV2o/JFwxyZgYJSz1NY6lpeUlJyufwh1w\n8nFWR9LaFcwRYsBYcRNRloi+R0Rfj1KgVgi1xFax6NgIZ2ZamwkHQQTs2uX87+66t3IdVxkHrRIG\nBtbtiEIiKI5NtR8SXj/JIALyeSwMqicAOpNI8SjUKRLaMYf4cbb96jyCg42N+1YAzwHYEpEsyeLw\n4eZNxjBwl6r1wTJVxY1F5LQdGgJeftkJVHCpV8aubdC1p7vL0bNne3qjJ+kUx6ZQhEJRj487Crla\ndVZmpRIwZabQc31DmF8503z+QhaAYoxlsygerYavqLNZ9ZiWlV9oGE31iGgXgBsA3BetOB3A1Pe1\nHEPSdzcwJp8HvvQlR3GfPg3cf3/DzArlcqMydlcJq6tO+9Onnf+PHxel3U2Mjzs/6K7Sq1adx+MK\nH29FzMHkI+cwQBsamg30D2Cyer36eq9+tblsulWhd59lYMCRy7sKlJVfuJj4DAJ4CMC1AK4H8HVN\nmxKAOQBzuVyuQ56Pltj4vobhQ217EHX+M+lC0OV+3Fp0sQLZbHNbTczBzPVDnP98nuku4vzn8zzz\ngxn/+IR2j6Gh9XiF+jiB+jgGb/yAoMRmXAf6cRPRjQD2MvM4EV0P4E+Y+Ua/1yQpi1oDNr6v2Ww4\nphI3656Jm6HOB1dooNv9uLX4eQl571O/mAPXVOGaWv7sz9qPT9ChilsQWsJmXJuYSq4DcBMRHQfw\nFQDvIKKZNuSLD5ucDRs3mvcbdMOZuBnKUlLIatwBVef97MVeU8uGDfq27SIufrEQqLiZ+VPMvIuZ\nCwA+AODbzHxz5JJFgU/OhiYuqAusNpHP+7vl5fNKDwCMjfnbrYXeo6RxB1Sdt4k5eOWV1mUSEknX\n+3FboRrsupmubvbjZe9es77rNxCPH3c8Beofi9IWpqacH3R37Ln+/iqvEtVkICquuEK/qhQXv1iw\nUtzM/ESQfTsR6DxHVINdN9NVuTOpKJed60xMAPv2ySxaaI+pKWBlxTGxraw4j3Xj+fBh4MQJp+2J\nE/bX0k1O8vnGLciTJ+1Wq0L0mO5i2hydzA7YRFgZ71rZie/RrGadBmn1KlGhG88jI+bj0iZbn2WG\nQRnv4WEzrtMzwF10CtdNt2qKaqCaHLbXEazpKcVtO4EgWncrzGYd5cysd8+zcdsTF79IsRnXXZ/W\nFUBjJj7d+2nFbam+323bgHPngosYiHtU5ITlDkhE9wO4EcCLzPyGoPaxuLm2kmo4gntaiJ6w3QGT\njWkmvlbclrwRiQcPrtuwdfZBsfl1E18E8K7QejONyvVr6z1vO27d1wX16yebkHxMp+Y2R0eXlKZL\nyaGhcK8rNr/YQIimEgAFAD80aes7rkOoSKO0OdcqKDUdV1yhPu+NvrS1ZQuxYTOuu19x6wa2yvYX\nNmLzi4VEKm6bvRVdW9PyePX263p79qZNdv3KfkyisBnX3W8qMTVNMIe/RPT6ZovrX+ogohIRzRHR\n3KlTp/QNbaJydW1NXVDdtl7XQV0OeF2/UmCja+l+xW0TQVZfiV0QDGDmMjMPM/Pwjh079A1t/Jx1\nbU2DvnRtbfuV/ZiupfsVt184uYqlJcdTRBDCxCYqV9dWlQ61T5My3zQMXtKsphNTm4rNEasfdz06\n+7ekT+1qEJKNG8CXAbwAYBnACQB/4Nc+cFyH4ROtOu+1Zbu+2e32KyQKm3GdDj9uHTZpXIWuIbVp\nXYWeprf8uP2wWb4KgiB0CelW3DZJpQRBELoEm2LB3UmxKIpaEIRUke4ZtyAIQgoRxS0IgtBliOIW\nBEHoMkRxC4KQSCShoZ70b04KgtB1uNma3fQrbrYKQHwNAJlxC4KQQCYmmnNmSbaKdQIVNxG9loi+\nQ0THiOhZIrq1E4IJgtC72CRb7EVMZtwrAD7OzLsB/AaAjxDR7mjFEgQhjZjaraWovD+BipuZX2Dm\np2v/vwzgOQA7oxZMEIR04a0y6JdlWbJV+GNl4yaiAoBrADwVhTCmhFHaTxCEzmJjt5ZsFQGYphEE\nsBnAEQDv1zxfAjAHYC6Xy0WW+jCM0n6S0bK7QYily2yOxKQr7lL8qgyaZq3VEVWW3DAwzahrM65N\nlXY/gL8FcLtJ+ygHeBil/fxK7YWR+liIFlHc3YlpXW/ATpna1F7u62v/ejbYTB5DVdwACMCDAA6Y\ndhrlALepjeA3MFSMjZl/qTKbjw9R3N2J6p7xm4GbElbt5SiwmTzajGsTG/d1AG4B8A4ieqZ27A3N\nVmNAvZ06o5G4frd5dNSxi+nQleArl83Pi5+pINihslvrsKmbHFbt5SiIyq3RxKvkSWYmZn4jM7+p\ndhxq77LBuMqaCLjllvWdaN0HPD/vtCUCZmf9+65WnXZ9fY6Sd38UdH2rzoufqRAFSd5MD0O2YtEp\nPrW66vzVTaKyWWB83LlH3Xt1fFzdNqray2EQmVuj6dTc5mh3SWmzpOrU4WVoSN1uaKitty4YgJSa\nSpJsfotKtpER9X10xRXq8zZmy927ze/v1Nm4WzlMBrjfxp7NJkanDq+MorjjI62Ku5XN9E4RlWy2\n97rOFq3SJ3427m73KomlWHClAvz+7wMrK+vnMhngNa8Bzp51PtoksmEDcP/9znIvk1HLSeQsA4Xo\nSGux4CSPqahk0/Xrh2l7v32uJOqYxBcL/qM/alTagPPlnzmTzA/U5dIl4NZaphYJyRXCJiljSmXL\n1slw2WXre0tEzp4RAOzc2Xh+505139u22cnm2ru9dm+VPdzPfm76vnXEvhdhOjW3OYKWlHGbPdo9\n3OWP13bV3++YSsSvO1qQUlNJEmzcOhl0NmfV0d+vPr91azR7Vzpbtk7mkZH2Pvuoviebcd2xAV5v\n54lb8XqPbNZOLtV7Ghpi3rAh/C9TaCatips5/qCuJO4vhX20G6wXlb0/cYo7iV4iKuVqMmh1m49J\n3lhKG2lW3HGTxIlV2IcqWM8vsM/7Y2rTrw0247ojNm5VsEoQfhsL9YE42Sywe7eZH2Y2C4yM6BPX\nqDKS1dPfD9x7r/o58esW0kAv7NGo3qPufW/b1pzRUKebOvnZdURxBykvr0IdGnI2GOpxP6x8Hnjw\nQScohtnZ5Hz2Wecvsz4aK5932jz+eGMAQH22MW9k19CQc7hK/uBBfXaypGwsCUI76NKphsHWreZ9\n6SKkVXh1hcvIiHlqWL/37Z10Mjcr746nnDWdmtsc3iWlrRmhHbNDXBs8SdhY6hUQoqkEwLsA/DOA\nHwP4pF/bXjCVMKvt7LamCO/G4BVX6PveurWx7datTltdcI7qerrkcDZ7Bqq2fqajsPcibMZ1RxS3\nrVKzSSSlIq4Nnrg3lnqFsBQ3gCyAnwB4HYANAL4PYLeufa8obhU2CZui2teJYx+pk9e0GdcdMZXY\nJkVv1+zgzYfQqeTrcV1XaJm3APgxMz/PzJcAfAXAe2OWKZG4Fda9ePeWojQZxFEVJ6mVeDoWgGOj\n1JL6YQmpYyeAn9U9PgEpy6dkagoYG1tX1Nms8/iBBzpXpSaOqjhJrcQTS8i7CZWK442ysODMtCcn\n4/+whGQQVsg7Ef1nAO9i5v9ae3wLgF9n5o/WtSnBqe6EXC537fz8fLuXFQQlNuNasx8bP8WiKGoh\nck4CeG3d4121c2swcxlAGXAmJJ0TTRD0RDLjJqJTAJIwNdkO4HTcQnhIokxAMuXSyZRn5h3tdk5E\nfQD+BcAIHIX9TwB+l5mf1bQ3HddJ/CzDRN5fNBiP60hm3GHcVGFARHNhLKnDJIkyAcmUK2qZmHmF\niD4Kp55qFsD9OqVda280rpP4WYaJvL/4SaypRBA6ATvVnCKv6CQIYRJLWldBEAShddKuuDXlf2Ml\niTIByZQriTKZ0K1ymyLvL2Yi2ZwUBEEQoiPtM25BEITUkUrFTUTHiegoET1DRNEVCbSEiLYS0UNE\n9CMieo6I3hqzPFfXPiP3OEdEt8UpkwsRfYyIniWiHxLRl4noVXHL5AcRvZaIvkNEx2py3xq3TGFC\nRK8iou8S0fdr729/3DJFARFlieh7RPT1uGXxI5WmEiI6DmCYmRPla0pEDwD4B2a+j4g2ABhg5pfi\nlgtwBiwcX+ZfZ+ZYffCJaCeAJ+EkfLpIRF8FcIiZvxinXH4Q0eUALmfmp4no1QCOAHgfMx+LWbRQ\nICICsImZzxNRP5zv51Zm/seYRQsVIrodwDCALcx8Y9zy6EjljDuJENEggLcD+AIAMPOlpCjtGiMA\nfhK30q6jD8DGWpDMAICfxyyPL8z8AjM/Xfv/ZQDPIUV5T2oJ7M7XHvbXjlTN+ohoF4AbANwXtyxB\npFVxM4BvEtGRWq6JJHAlgFMADtaWYvcR0aa4harjAwC+HLcQAMDMJwF8DsACgBcALDLzN+OVyhwi\nKgC4BsBT8UoSLjUzwjMAXgTwLWZO1fsDcADAJwCsxi1IEGlV3L/JzG8G8G4AHyGit8ctEJwZ5JsB\nTDPzNQAuAPhkvCI51Mw2NwH4i7hlAQAieg2c9KpXArgCwCYiujleqcwgos0AvgbgNmY+F7c8YcLM\nVWZ+E5ycLm8hojfELVNYENGNAF5k5iNxy2JCKhV3bcYGZn4RwCNw8i7HzQkAJ+pmKQ/BUeRJ4N0A\nnmbmX8YtSI1RAD9l5lPMvAzgYQBvi1mmQGq2368BqDDzw3HLExU1E9934FQPSgvXAbiptj/2FQDv\nIKKZeEXSkzrFTUSbaptDqJkifgvAD+OVCmDmXwD4GRFdXTs1AiApG1cfRELMJDUWAPwGEQ3UNsVG\n4NiME0tNzi8AeI6Z74lbnrAhoh1EtLX2/0YA7wTwo3ilCg9m/hQz72LmAhyz4beZObGrvDTmKvkV\nAI849xH6APw5M/9NvCKt8ccAKjXTxPMAPhSzPO6P2zsB/GHcsrgw81NE9BCApwGsAPgekh/Ndh2A\nWwAcrdmBAeDTtVwoaeByAA/UvI8yAL7KzIl2mUszqXQHFARBSDOpM5UIgiCkHVHcgiAIXYYobkEQ\nhC5DFLcgCEKXIYpbEAShyxDFLQiC0GWI4hYEQegyRHELgiB0Gf8fVA9c1mFoHQ8AAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c1617b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let us plot some of the features agains each other for different species\n",
    "\n",
    "f, axarr = plt.subplots(2, 2)\n",
    "axarr[0, 0].scatter(features[:,0][species == 0],features[:,1][species == 0],color = 'blue')\n",
    "axarr[0, 0].scatter(features[:,0][species == 1],features[:,1][species == 1],color = 'red')\n",
    "axarr[0, 0].scatter(features[:,0][species == 2],features[:,1][species == 2],color = 'green')\n",
    "#axarr[0, 0].set_title('Axis [0,0]')\n",
    "axarr[0, 1].scatter(features[:,2][species == 0],features[:,3][species == 0],color = 'blue')\n",
    "axarr[0, 1].scatter(features[:,2][species == 1],features[:,3][species == 1],color = 'red')\n",
    "axarr[0, 1].scatter(features[:,2][species == 2],features[:,3][species == 2],color = 'green')\n",
    "#axarr[0, 1].set_title('Axis [0,1]')\n",
    "axarr[1, 0].scatter(features[:,0][species == 0],features[:,2][species == 0],color = 'blue')\n",
    "axarr[1, 0].scatter(features[:,0][species == 1],features[:,2][species == 1],color = 'red')\n",
    "axarr[1, 0].scatter(features[:,0][species == 2],features[:,2][species == 2],color = 'green')\n",
    "#axarr[0, 1].set_title('Axis [1,1]')\n",
    "axarr[1, 1].scatter(features[:,1][species == 0],features[:,3][species == 0],color = 'blue')\n",
    "axarr[1, 1].scatter(features[:,1][species == 1],features[:,3][species == 1],color = 'red')\n",
    "axarr[1, 1].scatter(features[:,1][species == 2],features[:,3][species == 2],color = 'green')\n",
    "#axarr[1, 1].set_title('Axis [1,1]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Theano to do classification with the logistic regression\n",
    "\n",
    "We will use Theano to train a logistic regression classifier using this dataset. We leave validation and testing for the homework, but this part of the lab should form the basis of the code for your homework.\n",
    "First, let us define symbolic variables for the features and the classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Declare Theano symbolic variables for the features and the classes\n",
    "x = T.dmatrix('x') # vector of features \n",
    "y = T.lvector('y') # vector of outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create shared variables that contain the matrix of weights and the bias terms and initialize it with normally distributed values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# initialize the weight vector w randomly\n",
    "\n",
    "# this and the following bias variable b\n",
    "# are shared so they keep their values\n",
    "# between training iterations (updates)\n",
    "w = theano.shared(np.random.randn(len(features[0]),len(np.unique(species))), name='w', borrow=True)\n",
    "\n",
    "# w dimensions nfeatures x numclasses  for the softmax formulation each theta is the probability of each of the classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "No suitable SharedVariable constructor could be found. Are you sure all kwargs are supported? We do not support the parameter dtype or type. value=\"[ 0.  0.  0.]\". parameters=\"{'dtype': 'float64', 'borrow': True}\"",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m-----------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-35-bf87f2b69765>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# initialize the bias term\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtheano\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshared\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspecies\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtheano\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloatX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'b'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mborrow\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m//anaconda/envs/python3/lib/python3.6/site-packages/theano/compile/sharedvalue.py\u001b[0m in \u001b[0;36mshared\u001b[0;34m(value, name, strict, allow_downcast, **kwargs)\u001b[0m\n\u001b[1;32m    265\u001b[0m                     \u001b[0;34m' We do not support the parameter dtype or type.'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    266\u001b[0m                     \u001b[0;34m' value=\"%s\". parameters=\"%s\"'\u001b[0m \u001b[0;34m%\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 267\u001b[0;31m                     (value, kwargs))\n\u001b[0m\u001b[1;32m    268\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    269\u001b[0m \u001b[0mshared\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconstructors\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: No suitable SharedVariable constructor could be found. Are you sure all kwargs are supported? We do not support the parameter dtype or type. value=\"[ 0.  0.  0.]\". parameters=\"{'dtype': 'float64', 'borrow': True}\""
     ]
    }
   ],
   "source": [
    "# initialize the bias term\n",
    "b = theano.shared(value=np.zeros((len(np.unique(species)),)), dtype=theano.config.floatX, name='b', borrow=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"Initial model:\")\n",
    "print(w.get_value())\n",
    "print(b.get_value())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next step is the creation of the computational graph that computes the cost function for using all examples in the dataset. The likelihood is is built using the softmax formulation. In this formulation, we project the input vector into a set of hyperplanes that represent the different classes, and the distance from the vector to each of these hyperplanes relates to the probability of the input belonging to a given class. Write a graph that calculates the likelihood, the loss and the cost function, as well as the gradients of the cost function with respects to the weights and the bias terms. Use an L2 regularization approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 613,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Construct Theano expression graph\n",
    "p_y_given_x =                                   # Probability that target belongs to class i (softmax)\n",
    "prediction =                                    # The prediction of the model (class whose probability is maximal)\n",
    "loss =                                          # Loss function\n",
    "cost =                                          # The cost to minimize with regularization\n",
    "gw = \n",
    "gb =                                            # Compute the gradient of the cost\n",
    "                                                # w.r.t weight vector w and\n",
    "                                                # bias term b\n",
    "                                                # (we shall return to this in a\n",
    "                                                # following section of this tutorial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now write a Theano function that uses the graph to evaluate the predictions of the classiffier and the loss function for a single epoch, and evaluates the gradient using all examples of the training set in order to update the parameters (batch gradient descent). Leave the learning rate as an adjustable parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 614,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compile. What is a good learning rate?\n",
    "eta = 0.1 # Learning rate\n",
    "train = theano.function(\n",
    "          inputs=[x,y],\n",
    "          outputs=[prediction, loss],\n",
    "          updates=((w, w - eta * gw), (b, b - eta * gb)),name='train')\n",
    "validate = \n",
    "\n",
    "test = \n",
    "\n",
    "predict = theano.function(inputs=[x], outputs=prediction, name='predict')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now train your classifier using a certain number of epochs. Plot the evolution of the loss function for different values of the learning rate and the regularization parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 615,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_epochs = 50000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 616,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Train and every certain number of empochs (let us say 500), evaluate\n",
    "# the performance of the classifier on the validation set. Save the\n",
    "# validation loss each time you do so.\n",
    "errors = []\n",
    "errors_val = []\n",
    "for i in range(n_epochs):\n",
    "    pred, err =                      # train function applied to the training set\n",
    "    errors.append(err)\n",
    "    \n",
    "    if (i%500 == 0):\n",
    "        pred_val, err_val =          # validate function applied to the validation set\n",
    "        errors_val.append(err_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the validation losss\n",
    "plt.plot(500*arange(len(errors_val)),np.array(errors_val),marker='o')\n",
    "plt.yscale('log')\n",
    "plt.axis('tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now estimate the number of errors made by the classifier on the training set itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"Final model:\")\n",
    "print(w.get_value())\n",
    "print(b.get_value())\n",
    "print(\"target values for D:\")\n",
    "print(test_set[1])\n",
    "print(\"prediction on D:\")\n",
    "print(predict(test_set[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 619,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Estimate error\n",
    "# Function to count errors\n",
    "aa = T.dvector('aa')\n",
    "bb = T.dvector('bb')\n",
    "cc = T.mean(T.neq(aa, bb))\n",
    "f_err = function([aa, bb], cc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('Test error = ',100.*f_err(test_set[1],predict(test_set[0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [python3]",
   "language": "python",
   "name": "Python [python3]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}