StochasticGradientDescent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Gradient Descent\n",
    "\n",
    "\n",
    "Stochastic Gradient Descent (SGD) is an optimization method that is applicable to problems where\n",
    "the objective function is given as the sum of several terms \n",
    "$$\n",
    "E(w) = E_1(w) + E_2(w) + \\dots + E_N(w)\n",
    "$$\n",
    "\n",
    "Many models in machine learning lead to optimization objectives that have this form. For example, when $E_i(w)$ is the prediction error on example $i$ when using parameter $w$, the total error $E(w)$ is the sum of individual terms. \n",
    "\n",
    "SGD is an algorithm that can be used in this case. We make the following observations:\n",
    "\n",
    "* We can work with the average error\n",
    "$$\n",
    "\\frac{1}{N}E(w) = \\frac{1}{N}(E_1(w) + E_2(w) + \\dots + E_N(w))\n",
    "$$\n",
    "and the true average gradient is \n",
    "$$\n",
    "\\frac{1}{N} \\nabla E(w) = \\frac{1}{N}\\nabla  E_1(w) + \\frac{1}{N}\\nabla E_2(w) + \\dots + \\frac{1}{N}\\nabla E_N(w) \\equiv g_{[N]}(w) = g(w)\n",
    "$$\n",
    "\n",
    "* In a gradient descent algorithm, the computational burden will come from calculating the gradient for each data point in a big dataset. Remember that after taking a gradient descent, step the parameter will change and we will need to calculate the gradient again from scratch. If the data set is large, we could approximately calculate the true average gradient as\n",
    "\n",
    "$$\n",
    "\\frac{1}{N} \\nabla E(w) \\approx \\frac{1}{B}\\sum_{i\\in \\Omega}\\nabla  E_i(w) \\equiv g_\\Omega\n",
    "$$\n",
    "Here, $\\Omega$ is a random subset of $\\{1,\\dots,N\\}$ with $B$ elements.\n",
    "\n",
    "* We actually don't need to calculate the true average gradient $\\frac{1}{N} \\nabla E(w)$ exactly. In the extreme case, by just choosing $B=1$, that is just computing the gradient for a single data point, we get an unbiased but noisy estimate of the true gradient. Taking $B$ larger, reduces the variance of the gradient estimate.\n",
    "\n",
    "\n",
    "$$\n",
    "g = \\left\\langle{g_\\Omega}\\right\\rangle\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum Likelihood Estimation\n",
    "\n",
    "Suppose we wish to estimate the parameter $\\theta$, given conditionally independent observations\n",
    "$$\n",
    "p(x| \\theta) = \\prod_{i=1}^N p(x_i| \\theta)\n",
    "$$\n",
    "The maximum-likelihood estimator $\\theta^*$ maximizes the log likelihood function\n",
    "\\begin{align}\n",
    "\\theta^* = \\arg\\max_\\theta \\mathcal{L}(\\theta) \n",
    "\\end{align}\n",
    "\n",
    "The likelihood function can be written as the sum \n",
    "$$\n",
    "\\mathcal{L}(\\theta) = \\sum_{i=1}^N \\mathcal{L}_i(\\theta) \n",
    "$$\n",
    "\n",
    "$\\mathcal{L}_i(\\theta) \\equiv \\log p(x_i| \\theta)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Estimating the variance parameter of a Gaussian\n",
    "\n",
    "Assume for $i=1\\dots N$, we have conditionally independent observations \n",
    "$$\n",
    "x_i \\sim \\mathcal{N}(x; 0, \\theta)\n",
    "$$\n",
    "\n",
    "The likelihood contribution of a single data point is \n",
    "$$\\log p(x_i|\\theta) = \\mathcal{L}_i =  - \\frac{1}{2} \\frac{x_i^2}{\\theta} -\\frac{1}{2}\\log 2\\pi \\theta$$\n",
    "\n",
    "The total log-likelihood is the sum of individual data terms \n",
    "$$\\log p(x|\\theta) = \\mathcal{L}(\\theta) = \\sum_{i=1}^N \\mathcal{L}_i =  - \\frac{1}{2} \\frac{\\sum_{i=1}^N x_i^2}{\\theta} -\\frac{N}{2}\\log 2\\pi \\theta$$\n",
    "\n",
    "$$\\frac{\\partial \\mathcal{L}(\\theta)}{\\partial \\theta} = g(\\theta) = -\\frac{N}{2 \\theta} + \\frac{1}{2} \\frac{\\sum_i x_i^2}{\\theta^2}$$\n",
    "\n",
    "This problem has a solution in closed form:\n",
    "\n",
    "\\begin{align}                                                                                                                                                                                                                                                                                                                                                    g(\\theta) & = 0 = -\\frac{N}{2 \\theta} + \\frac{1}{2} \\frac{\\sum_i x_i^2}{\\theta^2} \\\\\n",
    "                                                                                                                                                                                                                                                                                                                                                    \\theta^* & = \\frac{\\sum_i x_i^2}{N}\n",
    "                                                                                                                                                                                                                                                                                                                                                    \\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recursive estimation\n",
    "\n",
    "Let's denote the estimate at $N$'th step as\n",
    "\\begin{align}\n",
    "\\theta^*_{N} & = \\frac{\\sum_{i=1}^N x_i^2}{N} \\\\\n",
    "N \\theta^*_{N} & = \\sum_{i=1}^{N-1} x_i^2 + x_{N}^2 \\\\\n",
    "\\theta^*_{N} & = \\frac{1}{N}\\sum_{i=1}^{N-1} x_i^2 + \\frac{1}{N} x_{N}^2 \\\\\n",
    " & = \\frac{N-1}{N} \\theta^*_{N-1} + \\frac{1}{N} x_{N}^2 \\\\\n",
    "& = \\theta^*_{N-1} + \\frac{1}{N} (x_{N}^2 - \\theta^*_{N-1})\n",
    "\\end{align}\n",
    "\n",
    "This recursion leads to the following algorithm\n",
    "\\begin{align}\n",
    "\\theta^*_{1}  & = x_1^2  \\\\\n",
    "\\text{for } & \\tau=2,3,\\dots & \\\\\n",
    "&\\theta^*_{\\tau}  = \\theta^*_{\\tau-1} + \\frac{1}{\\tau} (x_{\\tau}^2 - \\theta^*_{\\tau-1})\n",
    "\\end{align}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate sample data from the true model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy as sc\n",
    "import pandas as pd\n",
    "\n",
    "mpl.rc('font',**{'size': 20, 'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "mpl.rc('text', usetex=True)\n",
    "\n",
    "N = 1000\n",
    "x = np.random.randn(N)\n",
    "\n",
    "def g(th, N, s):\n",
    "    return -0.5/th + 0.5*s/N/th**2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exact recursive maximum likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAADOCAYAAAAuc4QFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd8W/W9//HX0fSQ94pHnMSOsxOy\nSCAQRgiUUWaBFrp7KWkp3YOOe1toy6W5ve3v3pbSBigtlxbKKKNAGFmQAFnOXs6wY8dxHO+9JZ3f\nH5Jl2bEdBRyP+P18PPSIdHQkH+crSx99zuf7+RqmaSIiIiIiIv2zDPUBiIiIiIiMBAqcRURERERC\noMBZRERERCQECpxFREREREKgwFlEREREJAQKnEVEREREQqDAWUREREQkBAqcRURERERCoMBZRERE\nRCQEtqE+gL4kJiaa48ePH+rDEBEREZFz2LZt2ypN00wKZd9hGziPHz+e3NzcoT4MERERETmHGYZR\nFOq+KtUQEREREQmBAmcRERERkRAocBYRERERCYECZxERERGREChwFhEREREJgQLnHs5/cDW3r9g4\n1IchIiIiIsPMsG1HN1QqGtqoaGgb6sMQERERkWFGGWcRERERkRAocBYRERERCYECZxERERGREChw\nFhEREREJgQJnEREREZEQKHAWEREREQmBAmcRERERkRAocBYRERERCYEC5z5c/t/vDPUhiIiIiMgw\nosC5D0crm4b6EERERERkGFHgLCIiIiISAgXOIiIiIiIhUOAsIiIiIhICBc4iIiIiIiFQ4NwP0zSH\n+hBEREREZJhQ4NwPr+JmEREREfFT4NwPt9c71IcgIiIiIsOEAud+KG4WERERkU4KnPuhjLOIiIiI\ndFLg3A+PipxFRERExE+Bcz/cCpxFRERExE+Bcz+UcRYRERGRTgqc+6GMs4iIiIh0UuDcD49HgbOI\niIiI+Chw7oe6aoiIiIhIJwXO/fBqyW0RERER8VPg3A/VOIuIiIhIJwXO/XCrxllERERE/BQ490Pt\n6ERERESkkwLnfqhUQ0REREQ6KXDuhzLOIiIiItJJgXM/1I5ORERERDopcO6HMs4iIiIi0kmBcz8U\nOIuIiIhIJwXO/VDgLCIiIiKdFDj3Q101RERERKSTAud+KON8bnr/SCW/fG2/xldERETOiG2oD2A4\nU8b53GKaJj99ZR9PbSoC4GhlEwWVTUSF2fjz588nKco5xEcoIiIiw5kyzv3wqB3dOeWdgxWBoHls\nfDhr8soprm5m34l67npyKy3tniE+QhERERnOlHHuh9ujjPO5ZPPRauxWg30PXI3darDreB1j48LZ\nVlTDsr9t49vP7uSPn5mLYRhDfagiIiIyDCnj3A/VwJ5bDpTWMzE5CofNgmEYzB4bS4LLyVXTx/Dj\na6by5r6TPLahoM/He7wmX/v7dj7++w38c9vxQTxyERERGQ6Uce6Hx1TgfK7o8HjZX1rP4omJvd5/\n1+IJbD9Ww/I3DzJ7bBzvHa7g0slJzBsXD0BZfSsL/3NNYP/vPr8Lr2ly2/yxg3L8IiIiMvSUcQ5i\n9giUlXE+d/z81f1UNLSxdFpKr/cbhsHyW2cxNi6c21ds5Hdrj/CFv2zlcFkDAM9uLQbgqmkpPPPl\nC1iUncD3X9jN1P94kwUPruY5//0fVn1rB09tLGRXcS0AmwqqeHvfyY/0nCIiIjKwlHEO0jPBrBrn\nc8OzW4/x1KYi7rp4AtfOTO1zv+gwO498eh7X/m4DAE6blS89uZWH75jLY+sLuGRSEo9+bj4A542N\n4Yn3jrLhcCWbj1bzg3/u5on3j5Kd7OLaGalcN6vvn9OTaZp86S9byS2qAWBxTiLvHanENOHfr5vK\nXYuzPsJvLyIiIgNFgXOQnmGyMs4jn9drsvzNgyyYEM/3Pjb5tPtPS4tmzXcvpba5A6vF4JMrNnLr\nnz7AwODBm2YE9otw2Lh3SQ73Lsmhud3NXU/m8kF+FXknG3h9dynHa6awKDuRH764m+WfmMWM9Jg+\nf+b7R6rILarhG1fkUNHQyjNbfNnrnGQXv3z9AAdPNvCja6cSH+n46P8hIiIi8qEpcA7Ss1RDfZyH\nvzsf24Tba/LIp+eS6Orqw/y5J7ZgtxhcPWMM1U3t3LFgLGF2a0jPmZ3kClz/f5+czT1/384tc9IY\nGx/R6/4RDhtPf/kCOjxemts8/Psre3nojTwy4yM4Vt3MLY98wOVTkpiUEsVlk5NJiHQwPjES8PWS\n/syfN2OzGNx9SRYup427L8nm/SOV3DY/g+VvHORvm4p4e38Z998wjZtmp6vrh4iIyBBR4Bzk1Iyz\n+jgPZ83tbj7IrwLgK09t45HPzCU5KoyWdg/rD1UAsCavHIfVwkV9TAo8nWtnpvLK1y4iO9l12n3t\nVgsxERb+3+3n0dbh4e39ZUxOiWJ8YgQbDlfy1r4yfr/2CAC/uGkG18wYQ25hNQD33zAdl9P35zgh\nMZIJ/sD6p9dP47b5Gfz0lb18+9ldPLu1mP++7Twy4noP4kVEROTsUeAc5JQaZ2Wch4Xtx2o4eLKB\nW+dlsHJPKUumJBMVZudkXSsAV0xJZk1eOQseXMPV08fw+UXjAfjJtVOxWw1mZsSSHBX2oX/+eWNj\nz2h/m9XC7+6Yw0MrD3DtzFQWZiXg8ZrsLK7hhW0lvLbrBP/x8l5+/WYe0eF2wuwW7liQ2efzTU2N\n5ukvX8AzW47x67cO8ok/fsBfv7iAqanRH/p3EhERkTOnwDmIibpqDEc/f3U/O4treXrzMfaU1AG+\n7ha3zE0H4N8WT+A7V03ijT0neXjdEd70d6O4esaYPssrzrYwu5UHbuyqibZaDOaNi2feuHh+ct1U\n/rzhKM9vK+Z4TQvxkQ6slv7LL+xWC5+7cDwLJyTw+Se2cOPD7/P5ReP44TVTT/tYERERGRgKnIP0\nzDh3qKvGkPjXrhOcrGvhMxeMI8Jho7XDtxR2Z9AM8Pb+Mt7eXwbAmOgwspJcTE+LIS02nA2HK0iJ\nDiMjLnxIjv90XE4b31yawz2XZ/O7NYeZlBIV8mMnj4nipa8tYvkbeTy24ShHyhv52PQx3DI3A4dN\n3SVFRETOJgXOQU4NnFXjPNhe3XWCbzyzA4C39pXx97sWUlrXym3zMthTUseczFg+dX4mL+8s4S/v\nFwIwJqarDOPOhZncubDvsofhxG618N2rTt/po6fUmHD+51NzmJEew3+9eZB1Byt4dfcJ/viZeUSH\n2c/CkYqIiAgocO6mZ6mGAuez4ytPbSPcYeW3t5+HYRiYpsn/rD5McXUzL+4oITspkm8uncQ3/7GD\nZU9to66lg6wkFw/dMhOrxcAwDKalRfPm3pPUt3QQ4RidL+O7FmfxqQWZvOqvmb7lkQ/4zW3nnVKT\n7fWaWFTOISIi8pGNzoijD8o4n32maQZqkA+ebCAjLjxQctFp6dQUbjgvjfL6Vn75+gEAUmPCsFm7\nShHsVgvrf3A59S0dg3fww5DLaeOOBZmMi4/g28/t5OZH3ufLi7P49pWTsFoMGlrdfPbPm3E5bTx8\n51ySopwcLmvAZrUEOneIiIhIaEIOnA3DuBWoNU1zdYj7zwXmAwVAFlAQ6mOHSs+K5na3apwHWm2z\nL9C1Ww32l9azv7Q+cN+nF2bictr4t4snAL6M6sm6Vh5/7yiZCadO8rNbLSQE9W4ezRZNTGTVdy7l\noZV5rFhfwIbDlSS4HGw4XBnY52P/s54Z6TGsP1SBYfi6kVyYncgXFo3XBEMREZEQhBQ4G4axFHgM\nuC3E/bOA5aZpXhm07XnDMApM0yz4UEc6CHougKKM84djmiamSa/lAcU1zQA8fOdcEiId5CRHUd7Q\nSn5FE1fPGHPK/j+5biq3zM1gamroE+hGq+gwOw/dMpMrpyXz7Wd3Bb6UPP65+WTEh/OL1/az/lAF\n6bHhLM5J5O39Zaw+UM6Kd/P51IJMvrFkYresvoiIiHTXb+DsD4DvA7YB1WfwvMuAFT22rQCWE2Lw\nPRR65pcVOH849/x9O8drWvjfT83GbrV0awl3vKYFgLFxEUxL8/Uhjomwk9NHZ4nOemYJ3ZIpKbz2\n9Yt5cXsJdy2eQKR/YZW/33UBRyubiIuwExvh4KFbTN7ce5J/bC3md2sOs3JPKePiI6hv7SAnJYpL\nJyVRWNnEuoPlOG1W7liQSU6Kq9vKiiIiIqNJv4GzPzu8DMAwjPvO4Hlv5dTAORdYdUZHN8h61ji3\nuxU4fxhv7PXVMC/5zbsA/POrFzJvXDwAxdW+jHNG/PBsFXeuGBsfwTeX5pyyPbiu2TAMrpmZyjUz\nU1m5p5Tlb+axJq+c6WnRPL35GE9vPgZAWkwYHtPkK3/bFrj9lcuyWZyTxI5jNcRHOpiVEUuY3UJV\nYzuJLif1rR0kRznPeHlwt8fLhsOVTEx2DVkPbhEZGYqrm4lwWFWyJ4NqwCcHGoYRi6+muVuG2jTN\nWsMwMAwja9iWa/QMnJVxPmO9fdn4j5f38eI9iwizWzle00JMuF1t04aZa2emcok/wzwjPYa9JXXs\nKK4lwm7l6hljsBgGb+8/SXl9G2/uO8lPX9l3ynM4bRba3F4sBnhNiHLamJ0Zy9eX5LBgQnxIx/G/\naw7z+7VHcNp8qynec1k2ydGhrfr4wZFKosPtTEuNVhcRkXOIaZocKG3g3UMVtHZ4mD02lqOVTfz8\ntf0AZCVFMn9cHPPHxTNvfByVDW3YrAYOq5UDJ+s5f3w84xMizviLvEhvzkZXjXjwBcp93J+Fb8Lg\nsDNc29GZphlYxXC41qA2tbn5yUt7qPN3ufjPm2eydGoyu4/Xcdf/5XLv0zv4zAWZFNc0D9uFSUY7\nl9PGjPQYAGakxwSud7pxtm+lxrsWT2BjfhVv7D3JtLRoaps7eGNvKTHhdpw2K6kxYbS5PbS7vWws\nqOL2FRtZOjWZOZlxeL0mM9JjuGxyErXNHaw/XEFshIOfvbKXBRPiyS2qIT02nAuzE/jbpiKeyy3m\nq5dmc9fiLMIdVsD3WsstquHCrITAoi81Te3c+fhmwJdV//LiLC6fksSY6DB9WMqgCJ7b8da+k7xz\nsIL54+K4anoKUWchUVBc3cyx6mamp0UTZrcSZrfi8ZojeqJvu9vL/a/uI9xuZaZ/IrPdaqGgspGt\nhTUAgS/mAEumJHP++Hi2Flbz1r4ynss93udzp8aEMTczjguy4lmck0R0uJ3oMNuw/UyV4etsBM6x\np99leOq5wvZwWTnw31/ey9/9p813/vRKYiMcQ3xEPqZpsq2ohnnj4vggv4qXd54I3JeVFElydBhL\np4Xxw2um8F9v5rH6QBlWi8GVU1OG8KjlozIMg0UTE1k0MTGw7auXZfe6b0u7hyfeP8qf3sln9YHy\nwPaZ6THYrQbbj/m+X9ssBkXVxzFNuO/qKXz1smzuvXwiv3ojj9+sOsSK9QXERthx2ix4TTha2USi\ny8Ht88cyMz2GXcd9q0reOi+DvSV1/PilPQDkJLu4cXYaU1OjSXQ52Xy0iulpMUxLjSY63D6igwwZ\nHg6U1vOTl/ZQUNmEx2syPiGSPSV1GAY8s+UY9hcNZqbHMCYmjN3H64gJt5OT7GL++Hg8XpNLJiWR\nGhOG02YJ6Uve7uO1PL7hKP/a1fV+azEg0mmjodVNZnwEV01L4dLJSczJjKOpzU1UmI0Ot0lMxPA9\n01fT1M6yp7axpbAah81Cu9uLzWIQ7rDS0Opm2aVZfOmiCbicNvaU1FHR0MbHpo/BYbPwVbLxek3y\nKxrJLarBa5okupyU1rYwIz2GA6X1bD5azfaiGl7fUxr4mVFOGxdmJzA1NZoEl4NF2QmE2a28c7CC\n8QmRzBsXh9vrJcxuxd4jwK5r6SDMbsFp832hN03zlPHrbds/thzDYjGYmOxiUkoULqe6Ao801vvv\nvz+kHR944IFvAq/df//9/WaLH3jggVRg2f333/9AL/fdDzzV13MYhnH3Aw88sOKBBx64u4GItJxL\nb2Z6WgwdHi93PrYJm9Vgamo0Le0ePvP4ZsLsViaPiaK+tYPP/XkLUWE2JiZHUd3Uzhf+soW4CAdZ\nSS7KG1r50l+3khTlZHxiJCdqW7jryVxSY8LJTIjgWFUzdz+Vy5joMF7d3fVHFeGw8vlF4zl4soGv\n/X07k1KiGBMTxr4TdXz96R1MS4smOSqMXcW1fPMfO5iVEUuiy8m2omq+/exO5o6LIz7SwaaCKr77\n3C4WTkggJsLOe4cr+f4Lu7hoYiJRYXbeOVjOff/czaWTk4h02li9v4wfvbiHK6amEO6w8m9P5gaO\naVpaDIfKGvjpK3v5+Kw07FYLL+04zv3/2sdNc9KxWgyezy3m56/u57b5YwHfm/ev3jjAJ+ZmAPDU\nxkJ+u+oQN8/xZRCfeO8of1h3hBv8GcVH1+fz2PqjfPy8NAAeeecIT35QyLUzU6ltbucnL+1l5Z5S\nnHYLn358C/+75jCv+t/E71iQyd6SOjxeM9Al4619JwNvEMXVLTjtFo6UN3LZ5GQA/uPlvWzMr+SS\nSUkA/OjF3WwrquFif2D2ved3sbekjguzfbe/9Y8dHCprZGFWAgD3Pr2dwqomzh/vKwf4ylPbOFHX\nEqirvuvJrVQ2tjMnMw6Azz+xhcbWjsBCIZ9+fBPtbi8z/RnWT67YCDCor72x8RFkxEWQX9HIV/+2\njaykSNJiw4f8tffm3lJ+8tJePjZjDE6blVd3nTij194L246zJq+M576yiJkZMZyXEUt5Qyser8ne\nE/WMiXYyJiac1d+5lI/PSmNPSR3Vze3cNCed2AgHx2qaMU24ICuBlOgwjlY2cbKula9dPhGbxcKz\nucW8vqeU3KIaLsiKJzspktgIBz+4egoz0mN452A5a/Mq+NeuE/xjazEbDlfy4vYSVqwv4In3jvLU\npiK2HK3m47NSMQxDr71z6LU3kO97AL9bc5i/bz5GToqLo5XNPPlBIY+sO8I/t5dQVNXMmJgwX9bX\nNLl25hjmj4sjJTqMJVNSKKhsYlNBFS6nnYRIJ1sLq3ljry8r/eQHhTy87giPrDvC/20sZGN+FU3t\nbh5dX0BeaT05KVF84S9beeiNPF7eWcLavHLeO+JrM3lJTiK3zM1gTmYsBRVNTE+LJjnayau7T/DC\nthIeeSefxzcc5Q/r8vnTu/n89YNC9vqDzu89t5OjVU2kRIex7mA5X396B4fKGmhsc7O72Pfls6nN\nTXlDG8XVTXz/+d2EO0J/7X3hL1tIcjmZEOJr787HNrGvtJ67Lp7Az66fzu7jtfzmttk8ePNMFk1M\n4KXtJZw/IZ70uAiqGtv5w7ojzMmMC7z2vvPcLpZOS+GSSUk0t3v44zv5LLskm6lp0TS0ull3sJw/\nf34+n7lgHBi+1qiLshPYWVzHqgNlrDtYwf9tLOKJ9wtZm1fOiztK+OM7R3jknXz+sO4If3wnn9zC\nGnYW1/KrN/J48PUD/OndfDYcruCP6/K5/9X9vLDtOI+uL+AfW4t5aGUev1t7hPWHKvjXzhP86o08\n/vpBIS/tKGHV/jKe3VrMivUFrD1QxnO5xfxtYxFzMuMob2jlgVf38cSGQiaP8U2af3Dlfp7efIzF\nOUm0dHj4/drDvLSjhKum+z5jl7+Zx+u7T7DUn5R68PX9rNpfxpIpvs/YB17dx7uHKs74M9eXbKjm\nZ//ax8aCKnJSogi3W/nWszs/9PtedVM7X/rrVjo83d/3PKZJZWMbJbXNfP/53YP+vrfvtT+X3n//\n/Y8SgrP2VccwjNh+yjV6ZZrmo8CjAPHjpg56uvfUrhrDI+M8ITGCo5W+SXXvHqzgsslJg/rzO9xe\nTtS2kHeynk888gFN7R7mZsayNq/8lH0fumUmEQ4rze3ubtttFgv/efNMrv/9eyREDo+MuQye+EgH\nH5+VxlMbC0mNCeevXzyfI+WNvHuogk0FVcRGOIiNcHD9ealsL6o95bEP3jwT8AUvBRWNgaXKf/na\nfvJONvDVy7K5MCuBX799EMOASyYlccmkJE7WtVDf6ua2eRlsPlrNgdJ6Wjo8LJwQT3F1Cy/tKOHt\n/WXc+qeNXDkthXa3l+Z2Dyvezfe/qQ+P9wAZei3tHlYfKOuW6e308xunU1rXSm1zOw/dMgvwBS9J\nUU5+eM0UfnjNFB741z4w4GfXT6fN7eEbz+zA5bQzd1wsdS0dvLDtOB6PSUFlEz95aS/gOxPz5MYi\n6lo6iHRYOVrZRGuHlx9dM4WdxbXMzIjhnssmAlBQ0cTccbHcfUk2jW1uPvP4ZiKdVq6YkkJ5Qxtv\n7i0l0mkjt7CG1/wJoiMVTfzl/cLA71FQ2RQ4uwkEzuR0Km9sY82BMuaNj8cdVMro9nhpaHXT6vb4\n9qtvZVtRDV9/ZgfT06JpbHNTWNXEYxvy2VZUAwYUVjXx5p5Syupb2VtSR2FVM4kuBz+5biqHyhqJ\ndNpwhdmwWIwBy8oahkFWkovLJyeTV9rAj6+bSpLLyRt7TvLohgJunJ2Gx2tiMQyezS3mwqwEnDYL\nBZVNbCus5kRdC+/nVxLhsJHocnDDeWlsL66lpcNDcpSTeePi8HhNdhbX4rBZuHlOOruO17ItvxKv\n6TsDdvHERKalRbO5oIrFOUlsKqhib0kdrW4v1/zvhm7He9Mf3u92e+4vfP0VDCDcYaW2uYM5mbEU\nVTXR4THZXFBFpNPGmVaZuj3mKe14dxbX8qf1Bd3mLb2wzVcKY7UYvHekktzCGiIcVjYdrQq0mq1u\n6uBIeSNeEw6VNVBU1cyB0npO1rWwp6SOV3eV0trh4XhtM/nljRyrbmZrYTW5RTV4vCbTUqPOSmnT\nQDJ6/mf1uaNh5APLTreIiX9yYA2QHTwJsK/tfZk/f76Zm5t7ut0GVEVDG+c/2PXrZSdFsua7lw3q\nMQDcvmIjF2Yl8O0rJwFwzf9uID02jHCHjY35VWz58RWDOvnpb5uK+PeXfW/kLqeNxrbuQfGqb1/C\n01uOkRkfwRcvmtDvc7V2eAizW8/asYqciQ6Pl79vKuKpTUXkVzQR5bTR5vYGJgbPyYzlnssm8pf3\nj1JW38oXLprADbPSCHdYqWpqIzVG9fqjwd6SOj7++/cA+NbSHHKSo3h7/0nGJ0SyKDuBBRPiB6yW\n3jRNDpc38va+k5yoa/V9Lo2P4+5LsnF7vBwqa2TymKgPXWZkmmagLWhJbQslNS1MSY0K1Ei3dniw\nWy3sPl6LgcGEpEjaOrxsKaxmU34VBZVNVDa2YRgwZ2ws09NiWLW/jJP1rThsFqanRVPZ2MbJulau\nnZnKgdJ6HDYL4XZf4F/Z2N7nsd1zWTY/uHrKh/q9Bkub24PdYjmjz+CWdg8NbR0kR/U+0blzzPed\nqMNutTA+IRLT9K17cKK2haykSGqaOgKfvRUNbRTXNLPneB0FlU2nPJ/dapDkchIdbic2wk5MuJ1x\nCZEsGB9PUpQTr2nS3O77krP+cAWPri/AahgkuBwkupy4PSYHyxq4eGIiyy7Nwmv6nrO0tpUTtS3U\nNHdQ1tBKfnkjlY1tjEvwZXZL61qBrjr0RJeDzPgI3F4Tt8ekuLqZmAg7181KZVthDbuP1+G0Wbh+\ndhpuj5f54+KZnRnLpD7a055NhmFsM01zfkj7DnTgHLTvbaZpbg/algVsM00zLpSfNxSBc3lDKwse\nXBO4PTY+nA0/WDKox+D1mmT9eCUAhb+6DoCLl6/l/PHxXDQxke89v4vXv3Ex09Ni+nuaAfXrt/L4\nw7p8ABbnJPKz66fz9Wd2cKC0nutmpfKHO+cO2rGInC17S+pYsb6AwsomfnTtFAoqmvjTu/mBICMt\nJowTda1EhdmYkOirY71uZirfWjqJicn997Yuq29lbV45dquFFe/mExth59tXTmJRdmK3/eqaO3hy\no+90/ZIpyUxMdpEeF05OskuTHAeQ12uyv7Se9YcrmDommotzEqlt7iAuwo7NamHNgTL+sO4IF09M\nZNKYKP76fiG5RTX88qYZvlP9o5hpmuwtqWdNXhlrDpSzp6SOMLsvuxrh8NUfl9W38pVLs7ljQeYp\nj29ze6hsbCc23E5lYxv5FY1MT4shyeVUN5wPoba5ncPljVQ1thHhsFHX0sG+E/VUNLRR19JBfUsH\nNc3tFFU199kpbMmUZKaMiaKysY2KhjZK61px2iw8u+zCM0p0VTa24faYxEbYcXvNU84UeLy+Rged\nk7rdHi8WwxgW434mgfPZKtVYjW+57e1B2+b6tw9fPScHDsGS29XNXd/Gm9rcgQkf0WE2Lpnk+5B9\n91DFoAbORVXNges3z0lnYrKLN765mJqmdmWP5ZwxIz2G398xJ3B7UXYinzp/LCv3nqSosomvXpbN\n/tJ6frfmCKsPlHHJpCTW5pWzck8pszJiSY8NZ05mLOMSImls66C2uYPCyia2HavhQGlDoDNOdJiN\n6qZ27nxsMznJLq6eMYaESAdJUWH86MXd1Le6mZAYyYMrDwSOJTbCzqLsBKanxXDVtJReFwxqd3ux\nWw0F2Pi+BJkmTE2N6tY1obndjcUw+NumIn75etf/r2H4+vgnupwsmZLE9mO1lNS0sLO4Fq8J4XYr\nv7hx+qgPmsFX7jAzI4aZGTF8a+mkMz6L6LRZSY/1namJdNoYlxB5mkdIf2IjHIFa407X++v0g7V2\neNh9vI7GNl/nqw6PSYTDimnCRRMTB2SidOJp+mlbLUa3nzNSO5qcSeAcTy8dM/wlGNvwLbHdWVh9\nH/A8/nplv2X+y7A1HFYOLKtvDVx/52AF184cQ2Obm6gwO8lRYUxPi+adgxWBurbBcKy6mcU5iTx8\nx1yiw7teMnGqVZZznM1q4YagD6FZGbE8/vn5lDe0kuRyUt3Uzp/fO0puYQ27S2q7zdjv5LBauHNB\nJtefl8b+E3VcOyuV6DA7/9x+nH/tPMHD644EFl+KCrOx8huLmZYWzaGyBv+p7TZ2Fdey7mAFK/ec\n5NdvHWRyShTR4Taiw+zcNj+D2WPjuOHh94gKs7FgQjz1rW5unZfBpTlJZz2b0+HxntJxYLB5vCaf\nfnwTEQ4b4Q4rr/treCMcVuZkxjJvXDzRYTb+682DWCzQ2uElJdrJX76wgBO1Laz0j1ur28Nb+8qo\na+ng5zdO58bz0jlW3Ux2ciSxwROrAAAaJklEQVQRDnU/6I2SJyNDmN0acj996d/pltyOBX6Er/dy\nLLDcMIwrgVWmab4QtGu30fAvdnKfYRjLga3+xy8ftguf+J2ycuAQBM7l9W2B689vK+ayyUl4vCZR\nYb6humJKMg+vO0J5fWvIC0OciYfXHsZqsbDskqzAB25RVTPXn5c6rFsZiQymzlrFBJezW03mwZMN\nFFc3ExVmIyU6jJYOD1NTu5aMD/7g+vTCcXx64TgKK5soq29lZ3Etk1KiAkvMT0qJCtT6fXqhL9NZ\n2djGc7nFvJNXQX5FI+0eL2v8k3QtBiRFOfnn9hKshsHru0vJSorkjvMzmZ4eTXaSi5QzfM/omU2s\naWqnorGNI+WNzMmM5TdvH+KFbcdJiXYyPS2G7KRIvCbYrAam6VvBcunUZP6xpZh1B8s5Ut7Iggnx\n3D5/LEumJIccdNU0tbOzuJaFWfGBAHbHsRr2najnkpwkCqua2FRQjctpo8PjJTspknuXTGTHsVpy\nC2t4eO3hQLvRT5yXgdvr5btXTiYzIYJpadEsndbVItPr9dUAZ8SFY7EYzIwYvLN7IjL8nW7J7Vp8\n2ePT7XNK3bK/vnn7qY8YvoIXQEl0OU6ZBDcYOjPO181MDUy4AAKzTG+ak87v1h7hlZ0n+PIlWQP6\ns1s7PPz324cAeHbrMZZ/YhZTxkRT19JBppY/FjmtyWOiAi2kQjU+MZLxiZGBFnf9SXQ5ueeyiYEz\nTh6vybq8cgoqGxmXEMnH/O2p2t1eVu4p5S/vHw2UfFgMuDgnCZvFwOM1SYsN48bZ6YxPiMQVZjul\nHnH38Vpu/MP7ZMSFc15GLNlJLh7bUBCYVNTp0klJJEQ6fDXDhypwe03fz/AvCPIfL3fte8N5aWwt\nrOaev2/HMGBuZhxXTx/DVdNTSHQ5ieyle4Jpmtz7zHbeP1KF02bhkklJTBkTxZ/fO9rtWMLsFjb/\n+AoiHNZAucrNc3yt6BpaOyhvaKO1w3PaMjeLxSAzQe93ItI7nXsK0pmR+OYVOXi8Jn98N3/Qj6HM\nn3G+bX4Gr+8p5a19JwECGeesJBfnjY3ln9uPD3jgvLPY1wbsgqx4NhVU88lHN/Hja33ZtMx41aGJ\nDDdWi+HPlnZfVMhhs3DTnHRumpNOaV0L+eVNbCyo5Nmtx3HaLFgssPloFc9sKQZ8M+YXTkjgsslJ\nZMSFk53k4unNx7AaBtNSo9l1vDbQwmzZJVlcPiWZvNJ62txevnjRhMBkH6/XDJypc3tNSmpaWJNX\nRlZiJBfnJOFy2vB4Td47UsnWo9WszSvnwZUHeHDlASwGzB8fz6WTkvggv5L6FjcLJ8QzISmS949U\n8eXFE+jwmLy17ySr9pdx/vg4vnPlZI6U+1peLZma3GvgDb7Ew3BvcSUiI4MC5yCdHUbSY8MprWsN\nzAAdzNXFimuaSYl2sig7kUiHlZd3lABdgTPAJ+am89NX9rGruDawkMJAyC2sBmDFZ+ez5kAZ33lu\nF/+5Mg+AccrAiIxIqTHhpMaEc3FOIt//WFdZSUu7hzf2llJW30ZNcztr88q7TZgDWDghnhWf9U00\nb2jtwGtCTLgvAL2glwy5xWIQZukqv+gtA2+1GFw6KYlLJyXxvY9NpqiqiXV55ZQ1tLH2QDm/fusg\naTFhjEuI5MmNhXT4Z+l/96rJhNmt/Oz6aTS0uYn2B8IXZp8+Uy8iMlAUOAcJ1DgbYLf5guUOjxer\nZXAmP7S0ezhQWs+4+EgcNgsXTUzk7f1lOGwWpgXVSd40J53lb+Tx5MZCfjt29oD9/D0ldWQlRhIT\nbueWuRm0dngDSxerVEPk3BLusHKLf1U9gB9fO5WS2hZO1rVyvKaZioa2bkHp2crYjkuI5Av+/u/3\nXT2FsvpW4iIcOGwWqpva+dfOEjITIgL10IZhBIJmEZHBpsC5Fwa+mfDgmyA4WLOGb37kffJONnDb\nPN+H2SfPH8uekjruXJDZbSJgdJid689L49VdJ2i9uf9WQKV1Lby2q5R/u3gC339hN1fPGMOV01L4\n3ZrDzMyI4XL/EpwAe0vqmTuuq1z99vkZ5BZVU9PU3ucpUBE5d6THhpMeG868cSG12z8rgicwxkc6\nAkG1iMhwoGgoSGfG2TCMQM1eh3twOmu4PV7yTjYA4PEfyBVTU7hiakqv+398Vhr/2FrMurxyrpmZ\n2ufz3vVkLvtO1LNqfxlbCqvZWVzDhdkJ/HaVbxLgzp9eiWEY3PfCbkpqW/jiReMDj7VZLfz29oHL\naIuIiIiMZCOz+/RZ0tlVw4BAX9IOz+AsgnKwrCFwvbP1VH8uyIonIdIRmLATbP+Jen752n4aWn0r\nCAFs8dcv51c0MeNnbwX2/fN7R1m1v4w3/ZMQl/YRqIuIiIiMdgqcg3RlnLsC5/ZByjgfKPUFzmu/\ne2lIp0ltVgvXzBzDmrwyGlo7ut33m7cP8vh7R3l8w1EAvtTHqc6FE+J5bEMBW45WATBvXBzjE9U9\nQ0RERKQ3KtUIEpgbaIDTHzi3uT19P2AAlTf4+jWPiQl9gYJb5mbwt03HeHVXKXcuzOTrz+wgKszG\n2oO+BRE62+ndMDuN731skm8J4KomNuZX0e7x8oVF47nqt+t5Lve4b+GEry4a+F9MRERE5ByhwDlI\nZzs6AwOnrTNwHpyMc0VDGy6n7YyWdZ0zNpYZ6dH8dtUhrpuVyqu7TgTu+/qSifx+7RHA1xEjwuF7\n7rTYcBZlJwb2+/7Vk/npK/uoaGg75flFREREpItKNYJ0yzj7O1UMXsa5jaQo5xk9xjAMfnnTTCob\n2/j9msOB7Z+9YBxfu3wiiS4HkQ4rcf0slf3J88cCBDp5iIiIiEjvlHEOEsg4G0EZ547ByzifaeAM\nMHtsLBdkxfP4e7565vuvn8bnF43HMAzuv2E6BRVNgeVne+O0WTnw86uxWwdvkRcRERGRkUiBc5DA\n5EAYklKNaWnRp9+xF19enMWmAl/XjEUTEwOB8sdnpYX0+HDH4PSpFhERERnJVKoRpFuphm1wSzUq\nGtpIcp15xhnotohJemz4QB2SiIiIiARRxjlIV8bZwGkfvIxzu9tLY5ub+EjHh3q8xWKw9ruXsqmg\nWiv8iYiIiJwlirKCBBZAMQgsY/1Ra5yLqppo7fAyeUxUn/vUtrQD9DuJ73SyklxkJbk+9ONFRERE\npH8q1QjSe43zRyvV+OXrB7j36e193l/V2EZRVTMAMREfLuMsIiIiImefMs5BglcOHKjJgTVN7eRX\nNNLa4QlksYN99W/bOVTuWzXwo2ScRUREROTsUsY5iNk1PTBocuBHC5wb29x4TThS3njKfV6vyd4T\nddQ2+5bMjg1XxllERERkuFLgHKQz42wxwG41MAxo6/hopRpN7W4A8k42dNvu9ZoUVTfT3N71/LHK\nOIuIiIgMWyrVCNJVqmEEFkE504xzh8fL3pI65mTGAdDU5guMD56s77bfd57bycs7T3TbFvchu2qI\niIiIyNmnjHOQQFcN/22nzUrrGWacn3jvKDc/8gEb86sAX6kGnJpxDg6aZ2XEYLMYRGohEhEREZFh\nS4FzkODJgcCHyjjX+OuV3z9SSYfHS7v/8Qd7BM7BvrU0J7BMtoiIiIgMTyrVCBK8ciCA037mgXNn\nZ4x9J+po8meb02PDKaltobyhleSoMLxeE4sBXhN+cPVklkxJYcmUlIH6NURERETkLFDGOYhpdpZq\n+CJnp83abx/nlnYP/7nyAHX+LDMQmOxXUNlEk//6RRMTANhVXAdARWMbXhN+ceN07rls4sD/IiIi\nIiIy4BQ4B+nMOBNcqtHPyoEbDlfw6PoCfvji7sC2Fn9N9LHqZqoa2wBYOCEBq8VgZ3ENAMdrfAue\nZMRHDOwvICIiIiJnjQLnIMErB4IvcG7tJ+Pc7vEF1W/sPYnX63tws7/9nGnCzuJaABJcDqaMiQrc\nPl7TAsDYuPCB/hVERERE5CxR4NyNv1TDX+QcG+Ggpqmj1z1f2VnCvU/vCNzefLQagOa2rkB7U4Gv\ns4bLaWP22Fh2F9fh9ZqBwDk9VhlnERERkZFCgXOQ4AVQABJdDqqa2nrd939WH+52+91DFYCvxjkn\n2UVKtJOVe05iGJASHcbssbE0tLkpqGzkeE0ziS4H4Wo/JyIiIjJiKHAO4g2Uavgi5wSXk6rG9kAZ\nRrCEoMVKspMiOVrpW1K7qd1NhD/DDHDXxRMYGx8RuL3jWC3Ha1pIj1O2WURERGQkUTu6IIGuGoGM\nsxO316S+tYPYiO6r+sUHBc5ZSS4KKpoAX6eNSIeVf79uGldMTeHWuRkAZCe5iHLa2HXcFzhPS4se\nhN9IRERERAaKMs5BAn2c/f8munzBcWVj+yn7JriCAufESIqqmvF4TZrbPUQ4rIyNj+D2+WOx+Os+\nLBaDWWNj2F5US0lNCxmaGCgiIiIyoihwDmL2iJwTXU4AKhtPrXMOt3cl67OTXbR7vBRWNdHc7ibc\n0Xsif/bYWPaX1tPu8ZKhUg0RERGREUWBcxCT7gugdGaVq3rJOLu9vlZ04XZroH559f4yCquaiexj\n0t95GbGB68o4i4iIiIwsqnEO1jk50J9xjgn3LZ9d33pqS7oOfw/n9T+4nPhIBy6njYfeyOtzf4DZ\nmV2Bs3o4i4iIiIwsyjgH6Vnj7HL6vlc0trpP2bfdbZIeG05SlBOrxeCCrPjAfbfNH9vr8ydHhfG9\nqyYxMdmlUg0RERGREUYZ5yCBlQP9KedIhw3DgIZeMshurxe71Qjc/tn106lqaudXt8xi8pioPn/G\nvUtyuHdJzsAeuIiIiIicdQqcg3TWOHcugGKxGLgcNhraTs04d3i82KxdCfux8RG8dM9Fg3KcIiIi\nIjL4VKoRxNujxhnAFWajoY9SDbtV/30iIiIio4UivyBmz350QFSYrdcaZ7fXiyOoVENEREREzm0K\nnIMEwubgjLPTRmMfpRrKOIuIiIiMHor8gnWWagRtigqz9zo5sMNtYlPGWURERGTUUOAcJLAASlDK\n2RXWx+RArzLOIiIiIqOJIr8gZi8Z5+g+Jgd2eLw4FDiLiIiIjBqK/IKYvXXVcPY+OVClGiIiIiKj\niwLnIF09NYK7athp6fAEltjupFINERERkdFFkV+QznZ0PTPOAE096pxVqiEiIiIyuijyC9LbAihR\nYb7AuWeds0o1REREREYXBc7d+DPOPRZAgV4CZ/VxFhERERlVFPkF6X1yoB3glF7OCpxFRERERhdF\nfkF6WzmwM+Pcc/XADo+JXaUaIiIiIqOGAucgXX2cuy+AAirVEBERERntFPkF6Vo5sGtboMY5KONs\nmiZur6nAWURERGQUUeQXpLeVA6N6qXHu8Ph2VKmGiIiIyOihwDlIbzXOYXYLNovRbfXAzsVQlHEW\nERERGT0U+QXpWgClK3I2DIO4SAcVDW2Bbe1uBc4iIiIio40ivyC9lWoAZCdFcqSiMXC7ucMDQITD\nOkhHJiIiIiJDTYFzkK7Jgd1D55zkKI6UNwYy0i3tvsA5XIGziIiIyKihwDlIXxnnickuGlrdlPvL\nNToD5wiHbRCPTkRERESGkgLnIL2tHAiQGR8BwPGaFgCa230TBVWqISIiIjJ6KHAOEuiq0SPnnBTl\nBKCy0Zdx7qxxDrMrcBYREREZLRQ4B+nqqtF9e6Kre+Dc2q7JgSIiIiKjTUhFuoZhzAXmAwVAFlBg\nmubq0zwmFrgbeAGoBuKBZcCq0z12qJh9bE9wOQCobGgHoFmBs4iIiMioc9rA2TCMLGC5aZpXBm17\n3jCMAtM0C/p5aDyw3H8BqAW+PFyDZiAQOffMONutFmIj7FQ0tgJdpRrhKtUQERERGTVCKdVYBqzo\nsW0FXQFxf64E4oBs0zTjTNN84QyPb1B1tqOz9Iyc8ZVrdGacW/yTA9WOTkRERGT0CCVwvhXY3mNb\nrn/7aZmmWXuazPSw4e0j4wyQ5HJS0djZjs63cqAyziIiIiKjR7+Bs79OOQtfjXKAaZq1/vuzzt6h\nDb6uPs6nRs7pceEUVzcD0NzhxmG1YNOS2yIiIiKjxulqnOOhK1DuRRa+CYN9yTIMozMzHQ9UD+dy\nja6VA0+9b3xCBC80tNHS7qGl3aMyDREREZFR5nSBc+xHeO5qgOBA2T+pkOEaPPe1ciBAZkIkAMeq\nm2lu96ijhoiIiMgoc9ZqDfy1zY/22NzvpELDMO42DCPXMIzcioqKs3VofQq0o+sj4wxQWNWkjLOI\niIjIKBRS4OyvdR4IBfjKN3p9PtM0HzVNc75pmvOTkpIG6Eeegc4FUHqJnMfF+zLOxdXN1La0ExNu\nH9RDExEREZGhdbrAubN+OT54Y1Dg22d9s2EYP+hlc+ckw2E5qTCw5HYvGefocBsup43jNS1UNbaT\nEOkY1GMTERERkaHVb+DsnxRYwKm1zvFAn23mOhdN6aXrRmcAPizb0/VX42wYBhlx4RyvaaGmuZ14\nBc4iIiIio0ooS26vxrfcdnAv57n+7b0yTbPAMIxlvQTWS4Ht/XTpGFKfvWAcn16YidXSW+iMP3Bu\nprqpnfhI5yAfnYiIiIgMpVBqnO8DbuuxbZl/O+Ar3TAMI98wjLuD9qkOzjj7yzuWAV/+CMd7Vlks\nBjarBaO3Wg0gIy6CvJMNdHhMlWqIiIiIjDKnzTibpllrGMZ9hmEsB7biq09e3ks2Ob7H414wDONW\nfx/nBHzlHreNlFUEe5MeGx64rlINERERkdEllFINTNPczqnLbgffXwvE9bJ9WPZr/rAy4oICZ5cC\nZxEREZHRRGtGn4GMuIjA9dSYsCE8EhEREREZbAqcz0BwxjknOWoIj0REREREBpsC5zMQG9G16Elf\nnTdERERE5NwUUo2z+BiGwc+un6Zss4iIiMgopMD5DH3xoglDfQgiIiIiMgRUqiEiIiIiEgIFziIi\nIiIiIVDgLCIiIiISAgXOIiIiIiIhUOAsIiIiIhICwzTNoT6GXhmGUQEUDcGPTgQqh+DnyuDSOI8O\nGufRQeM8OmicR4ehGOdxpmkmhbLjsA2ch4phGLmmac4f6uOQs0vjPDponEcHjfPooHEeHYb7OKtU\nQ0REREQkBAqcRURERERCoMD5VI8O9QHIoNA4jw4a59FB4zw6aJxHh2E9zqpxFhEREREJgTLOIiIi\nIiIhUOAsIiIiIhIC21AfwHBhGMZcYD5QAGQBBaZprh7ao5JQGYZxK75xy/b/u8I0zRd67HPaMdbr\nYOQwDCMWWG6a5rIe2zXO5wD/GH0SqAIS8P1NF/S4X+M8ghmGsRSY67+ZAOSbpvloj300ziOM//O4\ntrcxGKjxHNIxN01z1F/8/+mremx7Hsga6mPTJaTxuxWYG3Q7FsgH7j6TMdbrYGRdgBXA8z22aZzP\ngYv/b3pFz/HWOJ87F3wB89Jexl3v2yP4AiwFanqO7UCO51CPuUo1fJbh+xAOtgJYPgTHImcuyzTN\n7Z03TNOsxTd2wWMayhjrdTBCGIaRBcT3cpfGeYTzn0l4zAw6k+DPYC0N2k3jPPItM3tkCE3fWcIr\ng/dB4zwiGIaRZRjGCnxBbXUfuw3UeA7pmKurBmAYRj5wpdn9NGAsUGOapjF0Ryan4x+nNcAV/oC5\nc3sWvqxztmmaBaGMsV4HI4dhGHf7r15pmuZtQds1ziOcYRjLAUzTvK/H9qzOMdM4j3yGYWwDbgse\nH//2VaZpXum/rnEegfxjcsoXo4Eaz6Ee81Gfcfb/Z5/yDakzCPMHYDJM+ccpy3/pVShjrNfByOGv\ni3yul+0a53PD3cDWnhuDgmaN87lhNbDKX6sKBM4sPO+/rnE+hwzUeA6HMdfkQP/p3uBsZQ9Z+IrP\nZZgyTTOul81L8U1OKOj8QzrNGBPCPnodDA+xpmnWGsYpiYVQ/pYJYR+N89CKBQr8ZxWq8Y1rtdk1\n2VfjfA4wTfM+/5fgbYZh3Id/PMyuyYEa53PLQI1nKPuc1TFX4Ox7k5ZzzzLgIf/1UMZYr4MRwDCM\nW80e3VKCaJxHuKBsUVZQAIVhGMsNw4j3b9M4nyNM05xnGMYqfLWp24Ergu7WOJ9bBmo8h3zMR32p\nhpx7OjNVpmn+11Afiwwc/ym6vrIMcm7o60PxWTTZ65zjf6++D9+EwCx82WeVV8iwpsDZz/+hLCOc\n/013Wefkkh73nXaM9ToY1m7vOdmkNxrnEa2gx78A+LvmxAYHVRrnkc3fgWG1aZrb/X/XE/CN+7Ye\n+2mczyEDNZ5DOeYKnLveoLu1tgoaFNVHjSzL6X66D0IbY70OhjH/BKLc0+ymcR7hguoWQ6lf1DiP\nUJ3jENwVwTTNWn/Co8Bf+6xxPrcM1HgO+ZiP+hpn/ySjAk49RRiPf3LZEByWfAj+DMZ9PScNhDrG\neh0Ma/OBbMMwPhm0bS6Q5W9fttU0zRc0zueEwEpgvd2nv+dzwnx87UJ7s4KuCcAa53PEQI7nUI+5\nMs4+q/H9IQeb698uI4C/Vm55j76OS4NO7YYyxnodDFOmaT5qmuZ9wRdgFbDdf7tzwqDGeeRbQdcy\nzEDgjEPwh6LGeWQrALL7uC8W30RB0DifawZqPId2zD/q0oPnwgXfH2rP5RtXoSU7R8QF/zKt/j+c\nzstSui/Re9ox1utgZF3wleX0XHJb4zzCL/7xye9lfG7VOJ87F3xfkHouuZ2FLwGicR7BF3xLbt/a\ny/YBGc+hHnOtHOjnz2h8El/j/Sx8mSx9Yx3mOlcL6uPuAtM0s4P2Pe0Y63Uw/PnPItwH3I7vDfRR\nfF+Stvvv1ziPcEFjnI8vM/n8hxlDjfPw5j9TmA1U+TfVmkFtCP37aJxHAP9n8Y/w/f/fiu+swmp8\nAe4LQfsNyHgO5ZgrcBYRERERCYFqnEVEREREQqDAWUREREQkBAqcRURERERCoMBZRERERCQECpxF\nREREREKgwFlEREREJAQKnEVEREREQqDAWUREREQkBAqcRURERERCoMBZRERERCQE/x8mnjZEzhYQ\nWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119a634e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "th_star = np.mean(x**2)\n",
    "th = 0\n",
    "TH = []\n",
    "for n in range(N):\n",
    "    th = th + (x[n]**2 - th)/(n+1)\n",
    "    TH.append(th)\n",
    "\n",
    "plt.figure(figsize=(12,3))\n",
    "plt.plot(TH)\n",
    "plt.axhline(th_star, ls=':')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following panel illustrates the true gradient of the likelihood $\\mathcal{L}$ with \n",
    "estimates obtained on random subsets denoted by $\\Omega$. The expectation of the random gradients\n",
    "is the true gradient\n",
    "\n",
    "$$\n",
    "g(\\theta) = \\left\\langle{g_\\Omega(\\theta)}\\right\\rangle\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7b3d9469d16a45829d4812c2da0ade1f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in Jupyter Notebook or JupyterLab, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another notebook frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='B', max=500, min=1, step=10), IntSlider(value=10, description='NumOfBatches', max=400, min=10, step=10), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display, HTML\n",
    "from matplotlib import rc\n",
    "\n",
    "\n",
    "\n",
    "def plot_batch(B=1, NumOfBatches=10):\n",
    "    th = np.linspace(0.4,2)\n",
    "\n",
    "    for J in range(NumOfBatches):\n",
    "        s = np.sum(np.random.choice(x, size=B,replace=False)**2)\n",
    "        plt.plot(th, g(th, B, s), ls='-', alpha=0.3, color='blue')\n",
    "\n",
    "    s = np.sum(x**2)\n",
    "    plt.plot(th, g(th, N, s), ls=':', linewidth=4, color='red',label='true')\n",
    "\n",
    "    plt.gca().set_ylim([-0.2,1])\n",
    "    plt.grid()\n",
    "\n",
    "    plt.xlabel('$\\\\theta$/Parameter')\n",
    "    plt.ylabel('$g$/Gradient Estimate')\n",
    "    #plt.legend(['$g_\\Omega$','$g(\\\\theta)$'])\n",
    "    plt.title('True gradient versus mini-batch')\n",
    "    plt.show()\n",
    "    \n",
    "res = interact(plot_batch, NumOfBatches=(10,400,10), B=(1,500,10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SGD\n",
    "\n",
    "* Local search algorithms have the form :\n",
    "\\begin{align}\n",
    "w_0 & = \\text{some initial value} \\\\\n",
    "\\text{for}\\;\\; & \\tau = 1, 2,\\dots \\\\\n",
    "& w_\\tau = w_{\\tau-1} + \\eta_\\tau s_\\tau\n",
    "\\end{align}\n",
    "\n",
    "$\\eta_\\tau$ is a learning rate and $s_\\tau$ is a search direction. \n",
    "\n",
    "At each iteration $\\tau$, minibatch SGD selects random $i$ from $1\\dots N$. Set the search direction \n",
    "$$\n",
    "s_\\tau = \\nabla\\mathcal{L}_i(\\theta_{\\tau-1})\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "# Minibatch SGD\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "At each iteration $\\tau$, minibatch SGD selects a random subset $\\Omega_\\tau$ of size $B_\\tau$ of $1\\dots N$. Set the search direction \n",
    "$$\n",
    "s_\\tau = \\frac{1}{B_\\tau}\\sum_{i\\in \\Omega_\\tau} \\nabla\\mathcal{L}_i(\\theta_{\\tau-1})\n",
    "$$\n",
    "\n",
    "When the batch size is one ($B=1$), the minibatch SGD algorithm is SGD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "g_\\Omega(\\theta)/B & = -\\frac{1}{2 \\theta} + \\frac{1}{2} \\frac{\\sum_{i\\in \\Omega} x_i^2}{B\\theta^2} \\\\\n",
    "& = \\frac{1}{2\\theta^2}\\left(-\\theta + \\frac{1}{B}\\sum_{i\\in \\Omega} x_i^2\\right)\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "09c4d2dca8de4de3ae8043e9c774d859",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in Jupyter Notebook or JupyterLab, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another notebook frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='B', min=1), FloatSlider(value=1.0, description='eta', max=5.0, min=0.01, step=0.01), IntSlider(value=0, description='seed'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "th_true = np.sum(x**2)/N\n",
    "\n",
    "def plot_path(B=1, eta=1., seed=0):\n",
    "    np.random.seed(seed)\n",
    "    theta = 1.5\n",
    "\n",
    "    MAX_ITER = 1000\n",
    "    \n",
    "    TH = []\n",
    "    TH.append(theta)\n",
    "    for tau in range(MAX_ITER):\n",
    "        eta_tau = eta/(tau+1)\n",
    "        x_omega = np.random.choice(x, size=B)\n",
    "        theta = theta + eta_tau*(-0.5/theta + 0.5*np.sum(x_omega**2)/B/theta**2 )\n",
    "        TH.append(theta)\n",
    "\n",
    "    TH = np.array(TH)\n",
    "    plt.figure(figsize=(12,3))\n",
    "    plt.plot(TH)\n",
    "    plt.axhline(th_true, ls=':',color='r')\n",
    "    plt.ylabel('$\\\\theta_\\\\tau$/Estimate at $\\\\tau$')\n",
    "    plt.xlabel('$\\\\tau$/Iteration')\n",
    "    \n",
    "    plt.ylim([0.75, 1.5])\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "res = interact(plot_path, B=(1,100), eta=(0.01,5,0.01), seed=(0,100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A model for outliers\n",
    "\n",
    "$$\n",
    "x_i \\sim \\frac{1}{2}\\mathcal{N}(x; 0, \\theta) + \\frac{1}{2}\\mathcal{N}(x; 0, 100)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEECAYAAAA/L9PCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd4leXBBvD7yd57kD3ZkEAIG1SW\n4oKKoILFgQIKtrV+LWhbv9qv+lWotl9rtQY3CgqyilKRIcqQkcEKhJW9CAnZOzl5vj9yYmPIOiTn\nPGfcv+vi0uvNOcnN4eTc7/M+7/u8QkoJIiKi3rJSHYCIiEwLi4OIiHTC4iAiIp2wOIiISCcsDiIi\n0gmLg4iIdMLiICIinbA4iIhIJywOIiLSiY3qAPrg4+Mjw8PDVccgIjIpycnJJVJK354eZ5bFER4e\njqSkJNUxiIhMihAiuzeP46EqIiLSCYuDiIh0wuIgIiKdsDiIiEgnLA4iItIJi4OIiHTC4iAiIp2Y\n5XUcRNR3lfVNSPguHXvPF2FEkDte+clIONpZq45FRoAjDiK6gaZFYuWGFPzz23S4Othi+8l8LHzn\nGGobm1VHIyPAEQdRP6tv0uBqRT1KqhtQUt2Ayvpm1DdpUNeoQXOLhI2VgLWVgL2tNdwcbODmaAtv\nZzsEejjC29kOQgjVfwW8vuciDl0uwavzRuKhcaHYdaYQKzemYP3RbDx1a5TqeKQYi4PoJtU3aXDh\nahXOFVTgXEEl0q9VI/t6La5W1t/097SzsUKwpyMG+rlgkL8rhgxww6hQDwS6OxisUFLzK/DP79Lx\nYHwIHhoXCgC4OyYAnyX64J2DGXh0YjgPWVk4FgdRLzU0a5CYWYajGSU4nlGK03nlaNJIAICbgw0G\n+rticrQPwr2dEOjhCB9Xe/i42MHNwRZOdtZwtLOGlRBokRKaFom6Jg0q65pRWd+EkqoGFFbUo6C8\nDtnXa3HpWhX2pV2DpqX1+/u52iMu1BOTB/rgloE+CPN21svfUUqJ3+88By8nO/zm7qE/+tovZgzE\n/LePYsPxbDw5NVIvP59MA4uDqBsVtU34+vxV7E8rwqHLJaht1MDaSmBkkDuWTI7A6FAPDA90R7Cn\no84jAlcHW/i5dv31hmYNLl6twqnccpzMKceJzFLsPncVABDq5YTbh/lj9ogBiAv1hJVV/4xGtp/M\nR3J2GdbeHwN3R9sffS0+3AvjIryw4XgOnpgSYRSH1EgNFgdRBw3NGuw9X4QdJwvw3aVraNJIBLg7\n4L7RQZg+xA8TIr3hbK//Xx17G2vEBHsgJtgDj0xsHQ1kXa/F4cvF+ObCNaw/mo13D2fCz9Ue98UF\n4YH4EET5utz0z6uqb8KfvrqA2BAPzB8T3Olj5o0OwvPbziI1vxIjg91v+meRaWNxEGnlltZi44kc\nbE7MxfWaRvi72eORieGYExuImGB35XvYQghE+DgjwscZiyeGo6q+Cd9cuIYvThfi3UOZSPguA2PD\nPfFAfAjujgmAk51uv95vfHMFxVUNeOeR+C5HMLNHDMCL/0rFF2cKWBwWTEgpVWfod/Hx8ZL346De\nkFLiyJXr+OBIJr65eA0CwIyh/vjphDBMifaBdT8dAtK3a5X12HYyH5sTc5FRUgNXBxssHBeKRyeF\nI8jDscfnH8+4jkXvHsf9cUFYOz+228c+8WEi0gorcXj19H47REbGQQiRLKWM7/FxLA6yRFJKfHux\nGH/bfxmncsvh42KPheNCsHBcKAJ78UFrrKSUSMwqw0dHs7A7tXU+ZPbwAVgyJRxxoZ6djpquVtTj\n3n8chqu9Df71zGS4Otje8Jj2dpzMx7ObTmHLUxMRH+6lj78GKdLb4uChKrIoUkrsPV+EN765grP5\nFQjycMT/3jcS948Jgr2N6Z9iKoTAuIjWSez88jqsP5qFT4/nYNfZQsSGeGDxhDDMGuoPdydbSCnx\nffp1PLf5FGoamvHJE+N7LA0AmDHUD7bWAnvOF7E4LBRHHGQxTuaU4Y9fnkdKTjnCvJ2wclo07hsd\nBFtr815AobaxGVuT8/DBkSxklNTASgAD3BzQqGlBSXUjInyc8eaiOAwLdOv191z83nHkldXhm/+6\nVfncD/UfjjiItArK67Bm9wX861QBfF3tseb+kbg/Lhg2Zl4YbZzsbLB4YjgeHh+GM/kV+ObCNeSV\n1QIAJkR4466YALjoeJbY7cMH4MUdqUgvrkZ0d+cUk1licZDZqm/S4K1v05HwXToA4Jlp0Xjqtiid\nPyTNhZWVwKgQD4wK8ejz95o11B8v7kjFnvNFLA4LZJm/QWT2ErNK8fzWM0gvrsE9MQF4/s4hCPZ0\nUh3LbAxwd0BMsDv2ni/CituiVcchA2NxkFmpqm/C2t0X8fGxbAR5OOKjJeNw6yBf1bHM0qyh/nh9\n7yVcq6yHn5uD6jhkQJZxkJcswuHLJbjjrwfxyfFsPD45HHt+eQtLQ49mDfcHAOxLu6Y4CRkaRxxk\n8hqbW/DanotYdzAD0X4u2Pr0JMSFeqqOZfYG+7si1MsJe89fxaLxoarjkAGxOMikZZbU4OefnsTZ\n/Ar8dEIofnf3MDjYmv71GKZACIFZw/zx8bFs1DQ0G2T9LjIOPFRFJmtLch7u/vsh5JbVImHxGLz8\nk5EsDQObNcwfjc0t+OYCD1dZEhYHmZyGZg1e2HYWv/r8NEYGueOrX0zFHcMHqI5lkcaGe2GAmwN2\nnMxXHYUMiGNLMilFlfV4+pNkpOSU4+nbovCr2webzEKE5sjaSmDu6EC8dygT16sb4O1irzoSGQBH\nHGQykrNLcc8bh3HhahXeXBSH1bOHsDSMwLzRwWhukfjidIHqKGQgLA4yCZsSc/DQumNwsrPG9hWT\ncXdMgOpIpDV4gCuGBrhhS0oezHHtO7oRi4OMmpQSr319Eau3nsWESG/sXDkFgwdwiQtjs2h8KFLz\nK3Eis1R1FDIAFgcZrcbmFjy3+TT+ceAKHhobgvcfGwt3p56X/SbDmx8XDE8nW7xzKFN1FDIAFgcZ\npcr6Jjz2wQlsP5mPX90+CH+aN9Lslz83ZY521lg8IQz7LxQhvbhadRzSM/4mktG5WlGPBf88ihOZ\npfjLA7F4ZvpA3vPBBCyeGA4HG2v8efdF1VFIz1gcZFRyS2uxIOF75JfX4aMl4zAvLlh1JOolX1d7\nrLgtCrvPXcXR9Ouq45AesTjIaKQXV2PB20dRWdeMDU+Ox+RoH9WRSEdLb4lEkIcjXtp5Dg3NGtVx\nSE9YHGQU0gor8WDCUTS3tOCzZRMQ2w83GyLDc7C1xv/MHY6LRVV4fc8l1XFIT1gcpNzp3HI8tO4Y\nbK2tsGn5RAwN6P29r8n4zBjqj59OCMW6gxk4cJFrWJkjFgcpdSavHD999zjcHW2xeflERPm6qI5E\n/eC3dw3D0AA3/GzjSVy8WqU6DvUzFgcpc66gAovfOwEPZ1t8tmwCQrx4a1dz4Whnjfcfi4ezvTUe\n/+AErlXWq45E/YjFQUpcKqrC4vdOwNnOGhufnIBAD0fVkaifBbg74r1Hx6K8rglPfJSE2sZm1ZGo\nn+i1OIQQM4UQa4QQq4QQCUKIZb18XpwQYpn2+cuEEDP1mZMMK6O4GoveOQ4bK4ENSznSMGcjgtzx\nxsLROFdQgWc2nkSzpkV1JOoHeltWve3DXkq5ut22ZCGEh5RybTfPiwSwRko5q922z4UQGVLKDH3l\nJcPILa3FoneOQ0qJjcsmIMLHWXUk0rMZQ/3xh7kj8OKOVPxm+1msuT+GF3SaOH2OOJZ3sm1fF9s7\nPi+hw7YEAGv6IxSpU1LdgMXvHUddkwafPDke0X5crNBSLJ4Qhp9Pj8bmpDy8todXlps6fc9xzOpk\nW3kPz5kPIKXDtiTtdjJRNQ3NWPJhIq5W1uP9x8bylFsL9MtZg7BwXAjePJCOD49wMURTprdDVVLK\nBZ1sno8bRxM/EEJ4AIgE8KO1maWU5UIICCEiuzpcpZ0/WQYAoaGhN52b+l9jcwue+iQZ5woqsW7x\nGIwJ81QdiRQQQuCPc0egpLoRf/jyPHxc7XFPTKDqWHQTDHZWlfaDPaW7+Q0AXkBrUXTx9ciuniil\nXCeljJdSxvv6+vYhKfWnlhaJVVtO49DlEvxp3kjMGOqvOhIpZGNthTcWjkZ8mCee23Qa36eXqI5E\nN0HvxSGEmC+ESAAwpotRSHtcZ8LMvLr7AnacKsCv7xiMB+JDVMchI+Bga413HxmLcB8nLFufjHMF\nFaojkY70XhxSyi1SyuUAVmvPqorT988k4/DxsWysO5iBRyeGYcVtUarjkBFxd7LFR0vGwdXBBk98\nmIQiXiBoUnosDiFEpA5/uhwxaA8/JQDY34ufyZGHiTt0uRgv7TyH6UP88N/3Dufpl3SDAHdHvP/Y\nWFTVN+HJj5JQ18jVdE1Ft5PjbddU6PD9EgF0N4exD4CHEGKmlHJfJ19vm/j2Qruzr9oVCa/jMAFX\nrlVhxYYUDPRzwd8Xjoa1FUuDOjc0wA1/XzgaT65PwnObT+HNRXGw4vvF6HVbHNozmHqal7iBtnCS\nAcyQUnY8tRboYi5De/ZURidf9wJQzgsAjV9pTSOWfJgEextrvPtoPFzs9XbiHpmJGUP98du7huLl\nXWl4fe9F/PqOIaojUQ/0NcfhgdbRQccP+razojorkzb7AMR32Ban3U5GrKFZg6c+TsbVynq888gY\nBHtyKRHqnSemRGDhuFC8eSAdW5PzVMehHuilOLSjjE2dfGk1gLVtIwchhIcQIr3DGlarceMoZ7l2\nOxkpKSVe3JGKE1mleH1BLEaH8loN6j0hBP5n7nBMivLG89vOIDm7THUk6obezqrSXq/xQPtFDgF8\n3n7tKi2vDs8rR+sZWGu0p/KuQuvaVTxMZcQ2HM/B5qQ8/Gx6NO6N5UVdpDtbayv88+ExCHB3xIoN\nybhWxTOtjJWQUqrO0O/i4+NlUlKS6hgWIzm7DA+tO4rJ0T5479GxnAynPkkrrMS8t77HiCA3bFw6\nAbbWvPuDoQghkqWUHacKbsB/EeqTa1X1WLEhGQHujvjbgzyDivpuaIAbXr1/JBKzyvDKrjTVcagT\nPOWFblqTpgUrN6Sgsq4Z21aMg7uTrepIZCbmjgrCmbwKvHc4EzHB7pgXF6w6ErXDEQfdtFd2pSEx\nqwxr5sdwtVvqdy/cOQQTIr3wwrazSM3nsiTGhMVBN2Xn6QJ8+H0WnpwSgTmcDCc9sLG2wj8WxcHL\n2Q5Pb0hGRV2T6kikxeIgnWWW1OCFrWcwJswTq+/kxVqkPz4u9njz4TgUltfj+a1nYI4n85giFgfp\npL5Jg5UbUmBr07o8Ns94IX2LC/XEr+8YjK9Sr+LjY9mq4xBYHKSjV3al4XxhJV5fEItAD0fVcchC\nLJ0aiWmDffHyl2mc7zACLA7qtX+fLcTHx7KxdGoEb8hEBmVlJfD6A6Pg5WyHlRtTUFXP+Q6VWBzU\nK9nXa7B6yxmMCvHAqtmc1yDD83K2wxuLRiOvrA4vbDvL+Q6FWBzUoyZNC37+6UkIAc5rkFJjw73w\n3KxB+PJMITYl5qqOY7H4CUA9+vv+yzidV4E198cgxIsr3pJaT98ahcnR3vjDF+eRWVKjOo5FYnFQ\ntxKzSvHmgStYMCYYd44MUB2HCFZWAq8tiIWdjRWe3XQKTZoW1ZEsDouDulRV34RfbjqFYE8n/H7O\ncNVxiH4Q4O6I/71vJE7nluON/ZdVx7E4LA7q0u93nkNhRT3++uAo3smPjM7dMQG4Py4Y/zhwBcnZ\nparjWBQWB3XqyzMF2JaSj2emRWNMGG/KRMbppTnDEOTpiGc3neIpugbE4qAbFFbU4bfbUzEqxAM/\nmx6tOg5Rl1wdbPHXB0Yhv6wOL+08rzqOxWBx0I9IKbF661k0Nrfg/x4cBRueektGLj7cCyunRWNr\nSh72ni9SHcci8FOBfmRzUi4OXirGC3cNQbiPs+o4RL3ys+kDMTTADb/ZfhbltY2q45g9Fgf9IL+8\nDi9/mYYJkV746fgw1XGIes3OxgqvLYhBWU0jXtp5TnUcs8fiIACth6ie33oGGinx5/mxsOItYMnE\nDA90x8pp0dhxqgB7zl1VHcessTgIAPBZYi4OXS7BC3cO4dXhZLJWTovWHrJKRVkND1npC4uDkFdW\ni1d2pWFipDce5iEqMmFth6zKaxvxex6y0hsWh4VrPUTVutLo2vkxPERFJm94oDt+Nn0gdp4uwO5U\nHrLSBxaHhduSnIfDV0rwPA9RkRlZMS0KwwLc8N//SuW9yvWAxWHBSqob8Mq/0xAf5slDVGRWbK2t\n8Or9I1FS3YA1uy+ojmN2WBwW7H++OI+ahmb8ad5IHqIisxMT7IHHJ0dg4/EcnMjkWlb9icVhoQ5c\nvIadpwuw4rZoDPR3VR2HSC+emzUIQR6OeGHbGTQ0a1THMRssDgtU09CM321PRZSvM1ZMi1Idh0hv\nnO1t8PJ9I5BeXIO3DqSrjmM2WBwW6C97LyG/vA6v3h8Dextr1XGI9GraYD/MiQ3EW99eweWiKtVx\nzAKLw8KcySvHB0cysWh8KMaGe6mOQ2QQ/33vMDjb2+D5bWfR0iJVxzF5LA4LommR+M32s/B2scfq\n2UNUxyEyGB8Xe/z2rqFIzi7D5qRc1XFMHovDgmw8no3U/Er87u6hcHe0VR2HyKDmjwnG2HBPrNl9\ngcuR9BGLw0KUVDfgz19fxMRIb8yJDVQdh8jghBD4409GoLK+GWu/5rUdfcHisBB/+vcF1DVp8Mef\nDIcQvGaDLNOQAW54fFI4PkvMxcmcMtVxTBaLwwIkZpVia0oenpwaiWg/XrNBlu3ZWYPg52qPF/+V\nCg0nym8Ki8PMNWta8OKOVAR5OPL+4UQAXOxt8Lu7hyE1vxIbjmerjmOSWBxm7sPvs3DhahVevGcY\nnOxsVMchMgr3xARgSrQP/vz1RRRXNaiOY3JYHGbsWmU9/m/fZdw22Bd3DPdXHYfIaAgh8Ie5w1Hf\npMGrX3GiXFcsDjO2ZvdFNDa34KV7OSFO1FGUrwuenBqJrSl5nCjXEYvDTJ3KLcfWlDwsmRKBcB9n\n1XGIjNLKadHwc7XHH744zyvKdcDiMEMtLRIv7TwHX1d7PMMJcaIuudjbYPXsITiVW44dp/JVxzEZ\nLA4ztONUPk7llmP17CFwseeEOFF37hsdhNgQD7z61QVUNzSrjmMSWBxmpqahGa9+dQGxwe6YNzpI\ndRwio2dlJfD7e4fhWlUD3jpwRXUck8DiMDNvfXsF16oa8Ps5w3lXP6Jeigv1xLzRQXj3UCZyrteq\njmP0WBxmJOd6Ld45lIl5o4MQF+qpOg6RSVl95xDYWAu8vOu86ihGj8VhRv7332mwsRJYxSXTiXTm\n7+aAldOised8EY5cKVEdx6ixOMxEYlYpdp+7iqdujcIAdwfVcYhM0hNTIhDk4YhXdqXx9NxusDjM\ngJQSL+9Kg7+bPZZOjVQdh8hkOdhaY9XswThfWIntJ3l6bldYHGbgyzOFOJ1bjl/dPhiOdryHOFFf\n3BsTiNhgd7y25yLqmzSq4xglFoeJa2jWYO3XFzBkgCvmxQWrjkNk8qysBH5z11AUVtTjvcOZquMY\nJRaHifv4aDZyS+vw27uHwpqn3xL1i/GR3rh9mD/++W06Sqq5em5HLA4TVl7biL/vv4xbB/li6kBf\n1XGIzMrqO4egrkmDv+27rDqK0WFxmLA3vrmC6oZmvHAXT78l6m9Rvi54eHwoNp7IwZVr1arjGBUW\nh4nKuV6L9UezsGBMCIYMcFMdh8gs/WLGQDjaWvOeHR2wOEzU63svwtpK4LnbB6mOQmS2vF3s8fRt\nUdiXVoSkrFLVcYwGi8MEpRVWYufpAjw2KQL+brzYj0iflkyOgJ+rPdbsvgApeVEgwOIwSa99fREu\n9jZ4+tYo1VGIzJ6jnTV+PmMgErPK8O3FYtVxjAKLw8QkZZVi/4VreOrWKLg72aqOQ2QRHhwbgjBv\nJ6zZfYFLkYDFYVKklFj79UX4uNjj8cnhquMQWQxbayv81+2DceFqFb44U6A6jnIsDhPy3aVinMgs\nxc9nRMPJjnf2IzKke0YGYFiAG17fcwmNzS2q4yjF4jARLS0Sf/76IkK8HPHQ2FDVcYgsjpWVwKrZ\ng5FTWotNiTmq4yjF4jAR/04txLmCSvxy5iDY2fCfjUiFWwf5YnyEF/62/wpqGy33/uT8BDIBzZoW\n/GXPJQzyd8HcUbyPOJEqQrTeKK2kugEfHMlSHUcZvRaHEGKmEGKNEGKVECJBCLGsF8/x0D4+Uvv/\nkdrvMVOfWY3ZluQ8ZJTU4Fe3D+ZChkSKjQnzxIwhflh3MANV9U2q4yiht+Jo+6CXUq6WUq6VUi4H\nsFwIsaqHp3oBWAMgHUAZgGQAiVLKffrKaszqmzT42/7LGB3qgVnD/FXHISIAz84chIq6JnxooaMO\nfY44lneybV8X2zuaBcATQJSU0lNKuaVfk5mQT45lo7CiHr++YzCE4GiDyBiMDHbHrGH+eOdQBirq\nLG/Uoe85jlmdbCvvzROllOVSyox+zmNS6ho1ePu7DEyK8sakKB/VcYionWdnDkRlfTPet8CbPemt\nOKSUC6SUqztsng9gkz5+nhBimRAiSQiRVFxsHssCbDiejZLqBjw7kwsZEhmb4YHuuHPEALx/OBPl\ntY2q4xiUwc6q0k6Mp0gp1/bi4ZFCiPnaP8uEEPN7eoKUcp2UMl5KGe/ra/o3NWobbUyO9sa4CC/V\ncYioE7+YORBVDc1495BljTr0fvmx9kN/FtA6CunFU0q1j/1hXkMI8bkQApY019E22nhrRpzqKETU\nhSED3HB3TAA+OJKJJ6ZEwNPZTnUkg9D7iENKuUV7RtVqIUSyEKLbT0Lt3Ma6DpsT0HqmlUXgaIPI\ndDw7YyBqmzRYd8hypmR7HHEIISJ1+H6lUspOJ7+llOVCiAQA+9F6xpQuMtB6+Mqjq+9vTjjaIDId\nA/1dMSc2EB99n4UnpkTAx8VedSS967Y4tKWhy55+IoDu5jD2AfAQQszs6roMIcSqTuZB2m69FQkg\nRYc8JoejDSLT8/MZA/HF6QKsO5iB39w1VHUcveu2OLSnw/ZmXuJHtIWTDGCGlLKzD3qPbp63Rgix\npcOpuG2foGY/FuRog8j0RPm64CejgrD+aBaWTo2Er6t5jzr0NcfhgdYP+Y4f9G2HvTodNWjLYnkn\n12/MROsZWWZ9mKp1tJHO0QaRCfrZjIFo0ki8/V266ih6p5fi0I4yOrteYzWAtW3FoF2LKr3DGlal\n7edVhBAeaL3afKk+shqT1tFGI34xg9dtEJmaCB9nzB0ViI3Hc3C9ukF1HL3S2+m4Usq12mswogBc\nBxAF4PNOzpjy6vC8LW3XcADwRuvoZYG5X0XO0QaR6VtxWzS2n8zH+0cy8es7hqiOozd6vY6jk5Lo\n+PVydHKGlSVdr9GmbbTxFkcbRCYr2s8Fd40IwPrvs7Hslii4O9qqjqQXvB+HEWho1mDdwQxMjORo\ng8jUrZgWhaqGZnx8NEt1FL1hcRiBrcn5uFbVgGemR6uOQkR9NDzQHdOH+OG9w5moaTDPuwSyOBRr\n1rTg7e/SERvsjklR3qrjEFE/WDktGmW1Tfj0hHnem5zFodius4XIKa3FimnRvN8GkZkYE+aJSVHe\nSDiYgfomjeo4/Y7FoVBLi8RbB9Ix0M8Fs4by7n5E5uSZ6dEormrA58l5qqP0OxaHQvsvXMPFoiqs\nmBYFK95LnMisTIz0RlyoB97+Nh1NmhbVcfoVi0MRKSXePHAFwZ6OuDcmUHUcIupnQgg8Mz0a+eV1\n2HEyX3WcfsXiUORoxnWcyi3HU7dGwcaa/wxE5mjaYD8MC3DDP79Nh6ZFqo7Tb/iJpchbB9Lh62qP\n+WOCVUchIj1pG3VklNTgq9RC1XH6DYtDgVO55Th8pQRLp0bAwdZadRwi0qPZwwcgytcZbx5Ih5Tm\nMepgcSjw1oErcHe0xaLxYaqjEJGeWVkJLL81CmmFlTh8pUR1nH7B4jCwS0VV2HO+CI9OCoeLvd5v\n+U5ERmDuqED4udoj4TvzWKuVxWFgb3+bDic7azw+KVx1FCIyEHsbayyZEoHDV0qQml+hOk6fsTgM\nqLCiDjtPF+DBsSHwdLZTHYeIDGjR+FC42Ntg3UHTH3WwOAzowyNZkACWTI5QHYWIDMzNwRaLxodi\n19lC5JbWqo7TJywOA6mqb8LG4zm4a2QAQrycVMchIgUenxwOAeC9w5mqo/QJi8NAPjuRi6qGZiyd\nytEGkaUKcHfEnFGB2JSYi7KaRtVxbhqLwwCaNC14/0gmJkR6ISbYQ3UcIlJo2S2RqGvS4JNj2aqj\n3DQWhwHsOlOIwop6LLslUnUUIlJsyAA33DbYFx8dzTLZJddZHHompcS6gxmI9nPBbYP8VMchIiOw\n/JYolFQ3YmuKaS65zuLQs+/Tr+N8YSWWTo3g0ulEBADaw9bueOdghkkufsji0LOEgxnwcbHH3FFB\nqqMQkZEQQmD5LVHIul6Lveevqo6jMxaHHqUVVuLgpWI8NimMixkS0Y/MHjEAoV5OSDDBCwJZHHr0\n7qFMONpa42EuZkhEHVhbCTwxJQInc8qRklOmOo5OWBx6crWiHjtP5+OB+GAuL0JEnZo/JhiuDjZ4\n38QuCGRx6MmH32dB0yKxZAov+COizjnb2+ChsSH4KvUqCsrrVMfpNRaHHtQ2NmPj8WzcMXwAwryd\nVcchIiP2yMRwSCmx/qjpXBDI4tCDrSn5qKxv5miDiHoU4uWEO4YPwKcnclDb2Kw6Tq+wOPqZlBIf\nHsnEiCA3xId5qo5DRCbgiSkRqKhrwraUfNVReoXF0c8OXS5BenENHp8UASF4wR8R9WxMmCdigt3x\nwZFMtJjABYEsjn72wZFM+LjY457YANVRiMhECCGwZHIE0otr8N3lYtVxesTi6EcZxdU4cLEYD48P\nhb0NL/gjot67a2QA/FztTeLUXBZHP/ro+yzYWgs8PCFUdRQiMjF2NlZ4ZGIYDl0uwaWiKtVxusXi\n6CeV9U3YkpyHe2MC4efqoDoOEZmgRePDYG9jhQ+OZKmO0i0WRz/ZnJiLmkYNHuf9xInoJnk522Fe\nXBC2peQZ9R0CWRz9QNMi8dF6OHGWAAAKQklEQVTRLMSHeWJksLvqOERkwh6fHIGG5hZsPJGjOkqX\nWBz9YH9aEXJL6zjaIKI+G+TviqkDfbD+aBaaNC2q43SKxdEPPjiShUB3B9wx3F91FCIyA49NCkdR\nZQP2nCtSHaVTLI4+SiusxNGM61g8MRw21nw5iajvbhvsh2BPR6w/mqU6Sqf4SddHH32fBQdbKzw0\nNkR1FCIyE9ZWAosnhOF4ZikuXK1UHecGLI4+KK9txPaT+bhvdBDvuUFE/eqB+BDY21jhYyNcNZfF\n0QdbkvPQ0NyCxRPCVUchIjPj6WyHe2MDsf1kPirrm1TH+REWx01qaZH45Fg24sM8MSzQTXUcIjJD\nj04MR22jBluT81RH+REWx006fKUEWddrsXgi7ydORPoxMtgdo0I88PGxbEhpPKvmsjhu0vqj2fB2\ntsPsEQNURyEiM/bIxDBkFNfgyJXrqqP8gMVxE/LKavHNhSI8ODaEq+ASkV7dNTIA3s52RnVqLovj\nJnyqXQpg0XiugktE+uVga40Hx4ZgX1oR8svrVMcBwOLQWUOzBpsSczF9iD+CPZ1UxyEiC/DwhNa5\n1A3HjOPUXBaHjnanXkVJdSMnxYnIYII8HDFzqD8+S8xFfZNGdRwWh64+PpqNMG8nTI32UR2FiCzI\nIxPDUVrTiH+fLVQdhcWhi/MFlUjKLsNPx4fBykqojkNEFmRytDcifZ2x3giuJGdx6OCT49mwt7HC\ngvhg1VGIyMII0bp+1anccqTmVyjNwuLopcr6Juw4mY85sYHwcOK6VERkePPiguFga4UNx9Xe5InF\n0UvbU/JR26jhpDgRKePuaIt7YgKx81Q+qhualeVgcfSClBIfH8tGbLA7YoI9VMchIgu2aHwoaho1\n2HmqQFkGFkcvnMgsxZVr1T+cS01EpMroEA8MGeCKjSfUTZKzOHrhs8RcuDrY4N6YQNVRiMjCCSHw\n8PhQpOZX4kxeuZIMLI4elNc2YtfZQvxkVBAc7bguFRGpN3d0EBxtrbFR0SQ5i6MH21Ly0djcgoXj\nuC4VERkHNwdbzIkNxM7TBahScJMnFkc3pJT49EQOYkM8eLMmIjIqi8aHorZRgx0KJslZHN1Izi7D\n5WvVWDQuRHUUIqIfiQl2x/BAN2w8nmPwmzyxOLqx8UQOXOxtcA8nxYnIyAghsHBcKNIKK3Eq17CT\n5CyOLlTUNmHXmULMHRUIZ3sb1XGIiG4wd1QgnOwMP0lusOIQQngIIRJ6+dg4IcQyIcTMtv/qO19H\nO07lo4GT4kRkxFwdbDF3VCC+OFOAijrDTZIbcsSxBoBXTw8SQkQCWCOlXCel3CelXAdguXa7QbRN\niscEu2NEkLuhfiwRkc4WjQtDfVMLdpzMN9jPNEhxaD/0eywNreUAOo5MEtBaPAZxMrccF65W4aGx\nHG0QkXEbGeyOkUHuBp0kN9SIYyaAvb187HwAKR22JWm3G8Snx3PgZGeNOaM4KU5Exm/R+FBcLKpC\nSk6ZQX6e3otDOz+xuZeP9QAQCaC0/XYpZbn2610ertLOhSQJIZKKi4tvOm9lfRO+OFOAuaMC4cJJ\ncSIyAXNiA+FsZ41NibkG+XmGGHF4tH3w94IX8J+i6ESXxaGdE4mXUsb7+vrqmvEH/zqZj/omTooT\nkelwtrfBvbGB+OJ0oUGuJNdrcQgh5kspt+jwFKVrlkspsfFELoYHumEkJ8WJyIQ8ODYEdU0afHlG\n//ck11txaA87qVm68SadK6hEWmElFo4LhRC8pzgRmY5RIR4Y7O+Kz5P0f7iqx4P4Op4GW9ruMNMD\n2lNpdSaE0OXwVr8ZHuiGzcsnYmiAq6F/NBFRnwgh8OcFMRjg7qD3n9VtcbRdU6HD90sEsFYIEYfW\nM6F0laH9rxfajVa0o5f2X9cLIQTGRfT2rGEiIuNiqDuUdlscUsoMAAtu4vvGA4gSQjzYblscgEgh\nxBoAiZ3NfUgpy4UQGbhxrsMLQLk2DxERKaSX8007O0QlhFgFYKyUcnUPT9+H1uJpfy1HnHY7EREp\nZsglR7w7btCuX5UuhFjWbvNq3DjKWa7dTkREiun9CjftPMlqAA8AaFvoMEFK2Tai+NGkgvZw1eq2\nQ1povXZjDQ9TEREZB70Xh/YDf7n2T8evlQPw7GR7Cm5cdoSIiIwA78dBREQ6YXEQEZFOhKHvVWsI\nQohiANl9+BY+AEr6KY4l4OulG75euuHrpZu+vF5hUsoeF/szy+LoKyFEkpQyXnUOU8HXSzd8vXTD\n10s3hni9eKiKiIh0wuIgIiKdsDg6d1OLM1owvl664eulG75eutH768U5DiIi0glHHEREpBMWBxER\n6UTvS44QEaC9R008Wu8pEwkgQ0rJFZ87ob3/zjIAWwCUonU9u+UA9vI1+w8hxHy03m7ihtdE3+83\nFgeRnrXdEE1KOavdts+FEBlcvLNTXmi9gVzbTeTKASxlafyHEGImgHfQyf2SDPF+46GqLggh5mv/\ncTr7WpwQYpkQYmbbfw2dzxhpl8lfJYSI1P5/pBBiDV8fLAeQ0GFbAnS7u6almYXWBVCjpJSend34\nzRJpf6cS0DqKKO3iYXp/v/Gsqk5oP+g+B7Cg416Ots0TOrY5gNWWvveofW3S221q21O06F96IUQ6\ngFnt3x/awzFlUkqhLplx0r6PIjnC6J72fbW8k88ovb/fOOJox1ja3MRxT7Ed7S/sDe8n7S0F2j4k\nifqFod5vnONop929QyCE6OqOg/NxY3EkAdirx2gmRfsmLVedw0h4Af/5xe1EJFonMOnHIrWTv0Dr\na1hq6TshvWSQ9xuLQwfdtbkQAkKISEs/XEU38FAdwASVAkD7otBO7oLl0SODvN9YHLrh3mPPuKdI\nfaL9/eq4bEaC9g/fS0aAxaEb7j12j3uKXRBCeHSzw0E9y0DrTglfx17Q9+vEyXHqN1LKcillZ3uK\nlnziQNsI1Kv9Ru1hz/ZfJy0hxKpONrcdHubJBN0zyPvNLEccOp45UKprM5v7Xk8/v34Wvaeonf/K\nwI2jVS+0XvXL4min7eI1IcSWDq9N2wchX69uGOr9ZnbF0fbG0+EpiQDW9vKx7dv8hw9Bc9p77Mvr\nJ4RYJaXs+Fq231NM6XtCk7QPrcs/tP/7x2m3UztSygwhxPJOPuBmAkixxJ2Pm6D395vZFYf2DXfD\nZfj99L3Nfu/xZl8/7il2azVaLyhtfxhvufYP3ai0/RmK2h2z5QCWqo1ldLzQ+byr3t9vZlccBsC9\nx05wT7Fr2h2O1UKINWgdobWtJWTJZdolKeUW7ZI/8wF4o/XDcQFfrx9K9AW0voc80LqzNgutC0Bu\nAQzzfuOSI10QQpShk+UytP9wn3dYcmQvWi/9t+g3tvYXPaXDnuJ+tL6OlnqYisjssDja6dDm89F6\neGUf2rW59nFxAB7Ef9o8hevqtNKWRyT+s6fIPWsiM8PiICIinfA6DiIi0gmLg4iIdMLiICIinbA4\niIhIJywOIiLSCYuDiIh0wuIgIiKdsDiIiEgn/w8iLGww6cR1YwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11a460208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def log_gauss_pdf(x,mu,v):\n",
    "    return -0.5*np.log(2*np.pi*v) - 0.5*(x-mu)**2/v\n",
    "\n",
    "x = np.linspace(-10,10,200)\n",
    "v = 0.1\n",
    "v0 = 40.\n",
    "m = 3.\n",
    "\n",
    "y = np.log(0.006*np.exp(log_gauss_pdf(x,m,v)) + 0.994*np.exp(log_gauss_pdf(x,0,v0)))\n",
    "\n",
    "plt.plot(x,y)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.autograd\n",
    "from torch.autograd import Variable\n",
    "\n",
    "X_np = np.array([1, 2.5, -3.2, -1.1, 2.2, 2])\n",
    "X = Variable(torch.from_numpy(X_np).double(), requires_grad=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.62786465]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEECAYAAAAxqm/oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHuVJREFUeJzt3XlcVOX+B/DPMxuLgIDgSqjgvosE\noaaWWmnXXELTNC010LTdFl/d373dFstSyzUxbdNuqWRplt60zH0DUtw3zBRNXEBF9pnn98cMhcoy\nIDPnnJnP+/WaF87S+J1vx88cnvOc5wgpJYiISHt0ShdARERVwwAnItIoBjgRkUYxwImINIoBTkSk\nUQxwIiKNYoATEWkUA5yISKMY4EREGmVw5JsHBQXJRo0aOfKvICJyOcnJyRellMEVvc6hAd6oUSMk\nJSU58q8gInI5QohT9ryOQyhERBrFACci0igGOBGRRjHAiYg0igFORKRRDHAiIo1SbYBfzStUugQi\nIlVTZYAfz7iGmCk/Y9r/juB6fpHS5RARqZIqA9zP04j7WtfFnA3H0XP6Rnz3Wzp47U4iohupMsBr\n+3nig0c64JvxMQj29cBzS/dg8Pzt2J9+RenSiIhUQ5UBXqxTw0CsnNAFUx9ui98vXUe/OVvw6jep\nuJidr3RpRESKU3WAA4BOJ/DInaH4ZVIPjOnSGInJZ3DPtF+xcHMaCs0WpcsjIlKM6gO8mJ+nEf/8\nRyusfa4bIkID8NYPh/DAh5uw8egFpUsjIlKEZgK8WJPaPvjsiTuxaFQkzBaJUZ/swtjPd+P3i9eV\nLo2IyKk0F+AAIIRAz5Z18L/nu+HVPi2w/cQl3PfBJry75jCyOe2QiNyEJgO8mIdBj3Hdw7FhUg/0\na18f8zeewL3TfsXKPZx2SESuT9MBXqy2nyemD2mPb5/qjLo1PfHs13sw7OMdOHb+mtKlERE5jEsE\neLGOoQH49qkueHtgGxw6dw19Zm7GlB8PcViFiFySSwU4AOh1AsOjG2LDpB54OCIECzalodf0jVid\nepbDKkTkUlwuwIsF1jBhamw7rHiqM2r5mDDxv79hxKKdOJ6RrXRpRETVwmUDvFhEaABWTeyKN/q3\nRuqZK+gz0zpbJaeAwypEpG0uH+CAdVhlZEwjbJjUA/07NMD8jSfQa/pGrNl3jsMqRKRZbhHgxYJ8\nPDBtcHssHxcDPy8jxn+ZglGf8iQgItImYc8eqBCiF4DeAC4BCAeQLKVcUNF/FxkZKZOSkm67SEco\nMluweMcpTP/pKArMFjx9TxPEdw+HyeBW32lEpEJCiGQpZWSFr6sowG3hDSnl+pJvDmCplPK98v5b\nNQd4sfNX8/DG9wfxw75zaFLbB1MGtkVU40ClyyIiN2ZvgNuzuxlfymPry3hcc+r4eWLu8Ah88ngk\ncgvMGJKwHa8kpiLzeoHSpRERlcve8YLepTyWVZ2FKO3eFnWw7oVuiO8WhsSUM+g5YyNWpJzhQU4i\nUq0KA1xKOVhK+cpND8cCWOqYkpTjbTJgct+W+H5iV4QGeuOFZXsxfOFOpF3g3HEiUp9KH7ETQsQB\nSKlo/FvLWtX3w4rxnfHWgDbYl34FD8zcjJnrjyG/yKx0aUREf7FrFgoACCFiYRtKkVKWOf5tC/g4\nAAgNDe106tSpaihTORlX8/DG6oNYnXoOYcE1MGVgW9wVVkvpsojIhVXbLJRS3tgfwM8AnpRSppT3\nWi3MQrHXhiMZ+L/v9uNMZi6GRIbgtb6tUNPbqHRZROSCqnMWyg2klFkAEmANcbdxT/PaWPd8d8R3\nD8M3Keno9cFGrN1/TumyiMiNVfWslfUA/IvniLsLL5Mek/u0xMoJXRDs44FxS1IwfkkyMq7lKV0a\nEbmhcgNcCBEmhMgUQkSU8RJ/B9Skem0a1MTKiV3w0v3N8fPhDPSesQnLk05zyiEROVVFe+D+ANJs\nt5LCbD/LHQN3ZUa9DhPuaYI1z96NZnV88FJiKkZ+sgunL+coXRoRuYlyA9x2kLK0+d6vAHhPSnlz\nsLud8GAfLI2LwZsD2iDlVCbu+2ATPtlyEmYL98aJyLHsXcwqDtZFrFxmMStHSM/KxT+/3YcNRy6g\nwx3+eC+2HZrV8VW6LCLSGIdNI6wMdwtwAJBSYtXes3h91QFk5xdh4j1NMb4HVzkkIvs5bBohlU8I\ngf4dGmD9C93Rt209fLD+KB6aswX7068oXRoRuRgGuIPU8vHAzKEdsWhUJDJzCtB/7lbMWHcUBUUW\npUsjIhfBAHewni3r4KfnuqN/h/qY9fMx9J+7FQfPXlW6LCJyAQxwJ6jpbcSMIR2wcGQkLmbn46E5\nWzBz/TEUmrk3TkRVxwB3ol6t6mDd893wj3bWsfEBc7fi8J/cGyeiqmGAO5m/twkfDu2I+SM64fzV\nPPSbvQVzfjmGIu6NE1ElMcAV8kCbuvjp+e54oE09TPvpKAbO24Yjf15Tuiwi0hAGuIICa5gwe1hH\nfDQ8AmezctFv9hbM33iCZ3ESkV0Y4CrQp209/PR8N/RsWRvvrjmMYQt2cE0VIqoQA1wlavl4YN7w\nCMwY0h6Hzl3FAx9uwjKucEhE5WCAq4gQAoMiQrDmubvRpkFNvJyYirjFybiYna90aUSkQgxwFQoJ\n8MZXT96F1/q2xMYjF/DAh5uw/uB5pcsiIpVhgKuUTifwZLcwrHq6C4J8PDD2iyRMXpGK6/lFSpdG\nRCrBAFe5FnX9sHJiF4zrHo6vd59Gn5mbkXzqstJlEZEKMMA1wMOgx6t9WmBpXAwsUmLw/O14/3+H\neSo+kZtjgGtIVONArHn2bjwcEYK5G04g9qNt+P3idaXLIiKFMMA1xtfTiPcHt8e84RE4efE6Hpy1\nGd8kn+F0QyI3xADXqL5t62Htc93QukFNvLh8L55bugdX8wqVLouInIgBrmH1/b3w1ZN34cXezbA6\n9Rz6ztyM5FOZSpdFRE7CANc4vU7g6Z5NsSw+BgAwJGE7Zv98jOupELkBBriL6NQwAD8+ezcebFsP\n09cdxbCPd+BsVq7SZRGRAzHAXYifpxEzh3bA9MHtcSD9CvrM3Iy1+88pXRYROQgD3MUIIfBwpxD8\n8MzdaFTLG+OWpGDyilTkFZqVLo2IqhkD3EU1CqqBxPGdMb5HOL7adRoD5m7F8YxspcsiomrEAHdh\nRr0OrzzQAp+PjkLGNevFlL/97YzSZRFRNWGAu4HuzYLx4zPWJWqfX7oXrySmIreAQypEWscAdxN1\na3riv2OjMfGeJliWXDykwmtwEmkZA9yNGPQ6TLq/OT5/IgoXs/PRb/ZWfJPMIRUirWKAu6FuzYLx\n47N3o12I9TT8l5bv5ZAKkQYxwN1UHT9PfDk2Gk/f2wSJKWfQf+4WHDvPIRUiLWGAuzGDXocX72uO\nL0ZH4VJ2AR6asxUrUjikQqQVDHDC3U2DscY2pPLCsr147dt9yC/ikAqR2jHACQBQ2zakEt8tDF/u\n/ANDEnYgnWupEKkaA5z+YtDrMLlvS8wfEYETGdn4x6zN2HzsgtJlEVEZGOB0iwfa1MOqiV0Q7OuB\nkZ/swpxfjsHC5WmJVIcBTqUKC/bBdxO64KH29THtp6OIW5yEKzm84g+RmjDAqUzeJgM+fKQD/vNQ\na/x65AL6zdmCA2evKF0WEdkwwKlcQgiM6twIS+NjkF9kxqB525DIszeJVIEBTnbp1DAAPzxzNyJC\nAzBp+V68vuoACs0WpcsicmsMcLJbkI8HFo+JwugujfHZtt8xctEuXMrOV7osIrfFAKdKMeh1+Fe/\nVpg+uD2S/8jEQ3O2Yn86x8WJlMAApyp5uFMIEsfFwCIlYudvw8o96UqXROR2GOBUZe1C/LFqYle0\nbVATz369B+/8eAhmzhcnchoGON2WYF8PfDn2Loy4KxQJm9Lw+Ke7kJVToHRZRG6BAU63zWTQ4a0B\nbfHOoLbYkXYJ/eduxZE/uTQtkaMxwKnaDIsKxddxdyGnwIyB87Zi/cHzSpdE5NIY4FStOjUMxPcT\nuyI82AdPLk5CwsYTkJLj4kSOwACnale3pieWxcegb5t6eGfNYbyUmMr1xYkcwKB0AeSavEx6zB7W\nEU1q+2Dmz8fwx6UcfDQiArV8PJQujchlcA+cHEanE3i+dzPMGtYRe85kYcC8rTjK624SVRsGODnc\nQ+3rY2ncXcgrtGDQvG3YcDhD6ZKIXAIDnJyiY2gAVk7ogtBAb4z5fDcWbk7jwU2i28QAJ6ep7++F\nxPEx6N2qDt764RAmr9jHFQ2JbgMDnJzK22TAR8M7YcI94fh692k88eluXM3jlX6IqoIBTk6n0wm8\ndH8LvBfbDjvSLmHI/O04m5WrdFlEmsMAJ8UMibwDnz0RhfTMXAyYy2VpiSqLAU6K6to0CMvHx8Cg\nExiSsJ0zVIgqgQFOimtR1w/fTuiCxkE1MPaLJHy585TSJRFpAgOcVKGOn/X0+25Ng/Dat/vxzppD\nsHBtcaJy2RXgQohYIcTLQogEIcQ6IUSsowsj91PDw4CPR0ZieHQoEjam4emvf0NeIddQISpLhWuh\n2MI6TUqZaLvvDyBZCBEopVzg6ALJvRj0Orw1oA1CA73xzprDOH8lDx+PjERADZPSpRGpjj174GFS\nypTiO1LKLABTASQ4rCpya0IIxHcPx9xHI5CafgWx87chndMMiW5RboDb9rYfsf0sab3t+TBHFUb0\nYLt6WDw6ChnX8vHwvG28yg/RTcoNcNvedpjtRuR00WG1sHxcDCxSYvD8bdj9+2WlSyJSjQqHUKSU\nASWHUGx6AciSUqY5piyiv7Wo64dvxndGkI8HRizciZ8O/Kl0SUSqUNVphPEA3qnOQojKc0egNxLH\nd0aLen4YtyQZX+36Q+mSiBRX6QAXQsQBuCylfK+s54UQSUKIpAsXLtx2gUTFAmuY8NWT0ejWLBiT\nV+zDrJ+PcUlacmuVCnDbQct4KWXvsl4jpVwgpYyUUkYGBwffdoFEJXmbrHPFB0U0wIx1R/GvlQdg\n5gk/5KYqe03MqQB6OqIQInsZ9TpMH9wewb4eSNiYhkvX8zFjSAd4GvVKl0bkVHYHuBAiAcArtpkp\nRIoSQmByn5ao7euJN1cfRFbObnw8MhI1PHidbnIf9p5KHwdgaslZJ0KIXpwHTkob07UxPnikPXae\nvIzhC3ciK6dA6ZKInKbCAC+x7om/ECLCdusFYDCnEZIaDOwYgnnDI3Dw7FUMXbADGdfylC6JyClE\neUfxbWdgZpbxdJqUMry8N4+MjJRJSUm3UR6R/bYcu4gnv0hCHT8PLBkbjZAAb6VLIqoSIUSylDKy\notdVeCamlFKUcSs3vImcrWvTICwZG43L1wswZP52nLiQrXRJRA7F9cDJpXRqGICv42JQYLZgyPzt\nOHCWl2kj18UAJ5fTqr4flsXHwMOgw9AFO5B8iuunkGtigJNLCgv2wfK/1k/Zhc3HeFYwuR4GOLms\nBv5eWBYfg4a1vDHmsySs3c9FsMi1MMDJpQX7emBpXAxaN/DDhP+mYHXqWaVLIqo2DHByeTW9jVg8\nJhoRof545qvf8N1v6UqXRFQtGODkFnw8DPh8dBSiGgfi+WV7kJh8RumSiG4bA5zchrfJgE8fj0KX\n8CC8lLgXX3NNcdI4Bji5FS+THgtHRaJb02C8umIfFu84pXRJRFXGACe342nUY8HITujVsjb+77v9\n+HTrSaVLIqoSBji5JQ+DHvOGd8L9revgP98fxMebuC4baQ8DnNyWyaDDnEcj8GDbenj7x0OY9+tx\npUsiqhSufk9uzajXYebQDjDoBd5bewRFZolnejZVuiwiuzDAye0Z9DrMGNIBep3AjHVHoRPAxHsZ\n4qR+DHAiAHqdwPux7SElMO2no9DrdBjfgysmk7oxwIls9DqBaYPbw2yRmLr2MAw6gSe78aqBpF4M\ncKIS9DqBGUOsIf72j4eg1wmM7tpY6bKISsUAJ7qJQa/Dh0M7wGyReGP1QRj0AiNjGildFtEtOI2Q\nqBRGvQ6zhnVE71Z18K+VB7CEZ2ySCjHAicpgMugw99EI9GxRG//8bj/XTiHVYYATlcNk0GHeiAj0\naB6Myd/uw/Kk00qXRPQXBjhRBTwMeswf0QldmwTh5W9SsSKFS9GSOjDAiezgadTj45GRiAmrhUnL\n9+KH1HNKl0TEACeyl6fRuhRtp4YBePbr3/DL4fNKl0RujgFOVAneJgMWPX4nWtbzw7glKdh24qLS\nJZEbY4ATVZKfpxGfj45Co1reGPt5EpJPZSpdErkpBjhRFQTWMGHJmGjU9vXA45/uwoGzV5QuidwQ\nA5yoimr7eWLJ2Gj4ehgwctEuHM/IVrokcjMMcKLbEBLgjSVjoyGEwPCFO/DHpRylSyI3wgAnuk1h\nwT5YMjYK+UUWDF+0A39eyVO6JHITDHCiatCirh8+fyIKmdcLMXzhDlzMzle6JHIDDHCiatL+Dn8s\nGhWJ9KxcjFy0C1dyC5UuiVwcA5yoGkWH1ULCY5E4lnENYz7bjdwCs9IlkQtjgBNVs+7NgvHhIx2R\n/EcmnvoyGYVmi9IlkYtigBM5wIPt6uHtAW2x4cgFTFq+FxaLVLokckG8Ig+RgzwaHYrMnAK8/78j\n8Pcy4vWHWkMIoXRZ5EIY4EQO9FSPcGReL8DCLScRUMOE53o1U7okciEMcCIHEkLgtQdbIiu3EB+u\nP4YAbxNGdW6kdFnkIhjgRA4mhMC7g9riSm4h/r3qAPy9jejfoYHSZZEL4EFMIicw6HWYPawjohsH\n4sVle7HhcIbSJZELYIATOUnxBSFa1PPF+C+Tsfv3y0qXRBrHACdyIl9PIz57Igr1a3ph9Ge7cejc\nVaVLIg1jgBM5WZCPB74YEwUfDwMeW7QLpy5dV7ok0igGOJECQgK8sXhMFMwWC0Ys2omMq1zBkCqP\nAU6kkCa1ffHpE1G4lF2AxxbtwpUcLn5FlcMAJ1JQhzv88fHISJy8eB1PfLYLOQVFSpdEGsIAJ1JY\nlyZBmDm0A/aczsL4JSkoKOLiV2QfBjiRCvRpWw9TBrbFxqMX8OLyvTBz8SuyA8/EJFKJoVGhyMwp\nxNS1h+HvZcQb/bn4FZWPAU6kIuN7hCMrpwAJm9IQ4G3EC/c1V7okUjEGOJHKvNqnBTJzCjDrl+Pw\n9zZhdNfGSpdEKsUAJ1IZIQSmDLQufvXG6oPw9zZiUESI0mWRCvEgJpEKGfQ6zBzaEZ3Da+GlxFSs\nP3he6ZJIhRjgRCrladRjwchItK7vhwn/TcHOtEtKl0QqwwAnUjEfDwM+ffxONAjwwtjPk3Dg7BWl\nSyIVYYATqVwtHw8sGRMNX08DRn2yCycvcvErsmKAE2lAfX8vfDEmGhYJjFi4E39e4eJXxAAn0owm\ntX3w2RN3IiunACM/2YmsnAKlSyKFMcCJNKRdiHXxq98v5uDxT3fjej4Xv3JnDHAijencJAizhnVE\n6pksjFuSjPwis9IlkULsDnAhRKwQopcjiyEi+zzQpi7eHdQOm49dxAtLufiVu7LrTExbcH8MYLBj\nyyEiew258w5k5RZgyo+H4edlwJSBbbn4lZspN8CFEGEAXgGQDICX0CZSmbhu4cjKKcS8X0/Aw6DH\nv/u1Yoi7kXIDXEqZBiAeAIQQrzilIiKqlJfub478IgsWbTkJk0GHyX1aMMTdBBezItI4IQT++WBL\nFJotWLApDSa9DpPu5zK07oABTuQChBB4vV9rFBRZMGfDcZgMOjzTs6nSZZGDVXuACyHiAMQBQGho\naHW/PRGVQaezLkNbYLZgxrqjMBl0GNc9XOmyyIGqPcCllAsALACAyMhIzm0iciKdTuD92PYoNEu8\nu+YwjHodxvCCEC6LQyhELkavE5gxpD2KzBa8ufogzBYL4rpxT9wV8UxMIhdk1Oswa1hH/KNdPUz5\n8TBm/3xM6ZLIAbgHTuSijLar+pgMOkxfdxT5RRa8eF8zTjF0IQxwIhem1wlMi20Pk16HORuOo8Bs\n4TxxF1KZAA8E4O+oQojIMYpnp3gYdFiwKQ3X8orwZv/WMOg5gqp1FZ1K7w9gMoAwWMN7qhCiN4B1\nUspEJ9RHRNVApxN4/aHW8PU0Ys6G47hwLR+zh3WEl0mvdGl0Gyo6lT4L1rVQiEjjhBCYdH9z1Pbz\nwL9XHcCIRTuxaFQk/L1NSpdGVcTfoYjczMiYRpj3aAT2pV/Bwx9tw5nMHKVLoipigBO5oT5t62Hx\n6ChkXMvHoHnbsPd0ltIlURUwwIncVHRYLSSO6wyTQYfBCdvx3W/pSpdElcQAJ3Jjzev6YuWELuhw\nhz+eW7oH76w5xKv7aAgDnMjN1fLxwJIx0Xg0OhQJG9MwfOEOZFzNU7ossgMDnIhgMugwZWBbTBvc\nHntPX0HfWZux5dhFpcuiCjDAiegvsZ1CsGpiFwR4m/DYJzvx1uqDyCvkVe/VigFORDdoWscXKyd2\nwYjohli45ST6ztqM3/7IVLosKgUDnIhu4W0y4M0BbbBkTDTyCsx4+KNt+PfK/biSU6h0aVQCA5yI\nytS1aRDWPt8NI+5qiMU7TuHe6b9i2e7TnKmiEgxwIiqXn6cRb/Rvg++f7opGQTXw8jep6DNzE9bu\nPwcpGeRKYoATkV1a16+JxHExmD2sI4osEuOWpKDfnC1YnXoWRWaL0uW5JeHIb9DIyEiZlJTksPcn\nImUUmS34bs9ZzP7lGE5dykG9mp54LKYhBne6A8G+HkqXp3lCiGQpZWSFr2OAE1FVmS0SGw5n4JOt\nJ7HtxCXodQJdmwRhUEQD3NuiNnw9jUqXqEkMcCJyquMZ17AiJR0r95xFelYujHqBOxsF4t4WtdE5\nPAjN6/pCr+OVgOzBACciRVgsEsl/ZGL9ofPYcDgDR89nAwB8PQ2ICA1Au5CaaFbHFy3q+qJRUA0Y\neWWgWzDAiUgV0rNysevkJez+PRPJv2fi+IXsv6YhGvUCYUE+CAnwQn3/4psnavt6wt/biABvE/y9\njfA0uteVg+wNcF7UmIgcqoG/FwZ2DMHAjiEAgLxCM9IuXMeR81dx5M9sHM/IxtmsXCT/kYmsMk4U\n8jLq4e9thJ+nEZ4mPbyMOngZ9fAy6eFp1Fv/bNTDaNDBoBPQ6wQMOgGDvsR9/Y3PCQEIWH8C1isW\nCeCGx4vv44b7JV5ney1KPFcsyMeE1vVrOq6xYIATkZN5GvVoVd8Prer73fLc9fwinM3KxYXsfGTl\nFCIzp8D683oBsnILcTW3EHlFFuQVmHExuwC5hWbkFpiRX2RGToEZRWaJIosFajjPqG/bupg3vJND\n/w4GOBGpRg0PA5rW8UXTOr639T4Wi0SRRcJssQa6Ndit9wvNFpgtEhKAlMU/AUBCSvx1XxbfL/Fn\n3PzcLe/xN39vx8/AYYATkcvR6QRMf814cd3xcx7+JSLSKAY4EZFGMcCJiDSKAU5EpFEMcCIijWKA\nExFpFAOciEijHLoWihDiAoBTt/EWQQAuVlM57oD9qjz2rHLYr8qpar8aSimDK3qRQwP8dgkhkuxZ\n0IWs2K/KY88qh/2qHEf3i0MoREQaxQAnItIotQf4AqUL0Bj2q/LYs8phvyrHof1S9Rg4ERGVTe17\n4EREVAYGOBGRRqluPXAhRASASABpAMIApEkp1ytblToIIfwBxAFIBHAZQCCAeADrSvbI3XsohIgF\nkFXaZ7anN+7Wv7L6xe2tdLZ+hQEIt/1MkFIm3vQa52xnUkrV3GwfYt1Njy0HEKZ0bWq42fojS9wy\nAcSyhzd81l62vvQqo3/l9sbd+mdHv7i93fjZYgFElLjvD+AEgDgltjO1DaHEA0i46bEEAFMVqEWt\negMIABAupQyQN33zw017KIQIE0IkwPoP43IZL7OnN27RPzv7BXB7u1mYlDKl+I6UMgvWz1qyB87b\nzpT+RrvpG+gEbvoGgvUbTipdmxpusP5ju2VPiT0stQel7VFW2Bt37F85/eL2duvnSgbgX0qfZHEf\nnLmdqWYP3DbedsvegLR+w0EIEaZEXVrCHpbNnt6wf5Xjbv2yfa4w261Uzt7O1HQQMxD4+0OUIgzW\nwX53F2Y7iAJYe3ZZ/v1rLXtYNnt6Azte42794/ZWgpQyoJSHe8F6EDitOHydtZ2pKcD9lS5AAy4D\nQIl/QBBCLBdCFD/GHpbNnt6wfzfi9mafeADv2P7s1O1MNUMoVDEpZZaU8uZTc135gBEpiNtbxYQQ\ncbD+VvKeEn+/6gLcNj5E9kuD9dfcv/rGHpbNnt6wf+Xi9mZjGy6Jl1L2LuU5p2xnagrw4jGfwJIP\nlviQLjWWVhVCiJdLebj4QEjJcTP28Fb29Ib9K4HbW4WmAuh502NO3c5UE+C2Af003Do+FAjbAQLn\nV6Uetm/7qaUcoS7eCNLYw7LZ0xv272/c3spnm0P/ys0HIp29nakmwG3Ww3pqaUkRtsfdmu1/anwp\n/3N7AUgpsSGxh2WzpzfsH7i9lcc27j21ZG+EEL1KfNk5bztTenJ8KRPZbz69dB1c9LTcKvQntmQv\n8PeJBTef2uvWPUQpp3zb2xt37F85/eL2VnpP4mxhW3zrBet6KE7fzlS3HrhtgZdHAOyGdZwtRbrw\nwjiVVWIhnVqwbgQ37AnYXuN2PbSNH06G9fPGwvor6npY/5GUnAZXYW/coX+V6Be3NxtbzzLLeDpN\nShle4rVO2c5UF+BERGQftY2BExGRnRjgREQaxQAnItIoBjgRkUYxwImINIoBTkSkUQxwIiKNYoAT\nEWkUA5yISKMY4EREGvX/x6EV1prZqVUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x120510940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEECAYAAACBeYufAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAH5pJREFUeJzt3Xl8VNX9//HXSQj7EgLITiBsInsI\nKFZrVbAqahWxWlfUClbrVrfauvzqt1bRWre6ACpWRaGgqK1aC66AVUgI+xaSAGEPhpBAyDZzfn/M\ngCFkSAKZnFnez8djHiQzN8NnTm7mPfecc8811lpERERciHFdgIiIRC+FkIiIOKMQEhERZxRCIiLi\njEJIREScUQiJiIgzCiEREXFGISQiIs4ohERExJkGrgsIdW3btrXdu3d3XYaISNhIS0vbba1tV5Nt\nFULV6N69O6mpqa7LEBEJG8aYTTXdVt1xIiLijEJIREScUQiJiIgzCiEREXFGISQiIs4ohERExBmF\nkIiIOKPzhEREqlBc5mFzXhG7Ckr4YX8JuYUl7N5Xyt4DZRSXeThQ6uFAme9WXOahpMyLx1q8XovX\nWv/X+L72Wrz2x6+ttUf9v4/+aE02qH6T6moAWPXoudX/R8dJISQiUS+3sISlOfmkb97Diq17ycrd\nz9b8A0ds1yDGEN80jiYNY2kS57s1joulTbOGNGwQQ2yMIcaYQ//6vsb3dYwhxkCsMRhjqq2puk0M\ndfEc7imERCTqWGtZta2Auat3Mm/NTlZtKwB8IdO3QwuGd2/NL9t2pXvbprRv2Zi2zRvRtnlDWjWJ\nq1GASM0phEQkauwrKWdO+lamf7eJtTsKMQaGdWvNfef2ZUT3BAZ0bkXjuFjXZUYVhZCIRLy9RWW8\n8k0mb/1vE/tKyunfqSWPXTKAn/fvQNvmjVyXF9UUQiISsYrLPExbuJGXv9pAQXE5FwzqyI2n9WBI\n13h1q4UIhZCIRKTUjXncN3s5Wbv3c2bfdtzz877079TKdVlSiUJIRCLKgVIPT322jmnfZtOpVRPe\nunEEp/eu0aVtxAGFkIhEjI279zPxrTTW7Szk2pGJ3HfuiTRvpLe5UKbfjohEhC/W7uSOGUuJjTG8\necMIftpHRz/hQCEkImHNWssrX2fx5GdrOaljS165ehhdE5q6LktqSCEkImHL67U89skaXluQzYWD\nO/HUuEE6zyfMKIREJCyVebzcN3s5c9K3cv1PuvPQmJOIidG063CjEBKRsFPm8XLL9CXMXb2Te87p\nw61n9tJ5P2FKISQiYaXc4+XOGUuZu3onf7qoP9ed2t11SXIcdD0hEQkbHq/lnlnL+HjFdh4c008B\nFAEUQiISFqy1PPjBSj5Yuo17f96XX5+e5LokqQMKIREJCy98sYF3F23mlp/15NYze7kuR+qIQkhE\nQt6s1Bz+Nnc9Y5M7c+/P+7ouR+qQQkhEQtr8jFweeH8Fp/VqyxNjB2kWXIRRCIlIyNqwax+3vL2E\nXic056Wrk2nYQG9ZkUa/UREJSXuLypjwZioNG8Tw2vjhtGwc57okCQKdJyQiIcfjtdw2I52cPUW8\nc9MpdI5v4rokCRKFkIiEnEn/Wcs363N5fOxAhndPcF2OBJG640QkpHy6YjtTvsnimlMS+dWIbq7L\nkSBTCIlIyMjK3ce9s5czpGs8D11wkutypB4ohEQkJBwo9XDL9CXExRpevEoz4aKFxoRExLmDS/Ks\n21nItPHDNREhiuijhog498/UHN5bsoXbzurNz/qe4LocqUcKIRFxatW2vTz04SpO69WWO87u7boc\nqWcKIRFx5kCph9vfTSe+SRzPXTGEWF0ZNepoTEhEnPnLJ2vIzN3P2zeeTJvmjVyXIw7oSEhEnPh8\nzU7e+m4TN53eg9N6t3VdjjiiEBKRepdbWMJ9s5fTr2NL7tGlGaKauuNEpF5Za7lv9jL2lZTz7hVD\naNQg1nVJ4pCOhESkXr313Sa+XJfLH87vR5/2LVyXI44phESk3qzfWchjH6/hzL7tuHZkoutyJAQo\nhESkXpSU+6ZjN2/UgCfHDdYVUgXQmJCI1JO/fraOtTsKee26FNq10HRs8dGRkIgE3YKM3Uydn801\npyRydr/2rsuREKIQEpGg2rO/lLtnLaVnu2b84fx+rsuREKPuOBEJGmstv39/OXn7S3ntuuE0aajp\n2HI4HQmJSND8MzWHz1bt5N6f92VA51auy5EQpBASkaDI3r2fP/1rNaf2bMOvT0tyXY6EKIWQiNS5\nMo+XO2ekExcbw9O/HEyMVseWADQmJCJ17vnPM1i2ZS8vX5VMx1a6SqoEpiMhEalTi7LzePHLDVw2\nrAvnDezouhwJcQohEakzBcVl3DVzKV0TmvLIRf1dlyNhQN1xIlJnHv5gJTsKipl980iaN9Lbi1Sv\nXvYSY8w4IN9aOy/A48nA5cAPQBtgsrU2q5rnHAUk+79tA2Raa6dUeHwSMBdItdbmH/+rEJGj+XDp\nVj5Yuo3fje7D0G6tXZcjYSLoIeQPi6nAZQEeHweMttZOrHDfZGBiVdv7H08GsNY+WfF5jDETKgRR\nMnCf/7HKT5Flre1Z+1cjIlXJySviwTkrSUlszS0/05+W1FzQQsgYkwTcD6QBeQG2iQemWmtbV7hv\nHDCqmqefWDG0AKy1s40xs4CDIZQFjPb/3xWPhEYBqbV4KSJyFB6v5e5/LsMCz1w+hAaxGmqWmgta\nCPm70yYCGGPuD7DZA/wYGgd/brYxZkk1T59ijEmqossuvsLXaVV1/xljsNZW9/wiUkOvfJ3Joo15\nPHP5YLomNHVdjoQZ1x9ZJgCLK99Z3XgQMA+Ye7BbDg4dQc2q8BxTKv9Qpe66gIwxE4wxqcaY1Nzc\n3Oo2F4laS3PyeWbuei4c3ImLh3R2XY6EIdfTV+KBLGPMBHzdZglAnrV29tF+yFp7v3+sKc1/lJXl\nvz9gwPi3r1E3nP95pgCkpKTYmvyMSLTZX1LOnTPSad+yMX++eIAuUifHxFkI+ceMAJIqz2ozxiRU\nd8RirR1mjJkLTAKWAGdX818mV5zIICLH5//+vZpNeUXMuOkUWjWJc12OhCmX3XHxAe6fiS9Yjsp/\n9HQ/vskHSfiOiqpcJdHfVVddF5+I1NAnK7YzY3EOvzmjJycntXFdjoSxao+EAr2xB5BXi3Nysir9\nC4C1dokxJj7AxIODNU0GJh183BjTA994UBpQ1QkKD1D9kZKI1MC2/AP8/r3lDO4az12j+7guR8Lc\nUUPIH0DVHpVUsBioUZeXtTbf34ccKLSSqOLoxT+t+7DJC/7gG22MSTPGjKo4K86/fbJOWBU5fh6v\n5c6ZS/F4Lc9fMYQ4TceW43TUEPK/0Vd5kmkdySJA2AS4DyAFyAzw2GSO7OYbReCgE5FaePmrDSzK\nzuPpywaT2KaZ63IkArj+GDOZH5feAQ6thpB/lGnaWUCgU7Lj8U1SqGg4Gg8SOW5LNu/hmXkZXDS4\nE2OTNR1b6kZ9hVACVU9EmMKRy/NMAm46+I1/fCjTPxGBCuNAh62q4O86bFNFeCURYMUGEamZwuIy\n7piRTsdWjfnzJZqOLXUnmMv2xOObEJCEL4AmGWNGA3MPngfkHxca7Z9okInvCGdSFSsdJFT8xlo7\n0X9C6SR8i56C7+ipqpUZ8lAIiRyXhz9cxdY9B5h180haNtZ0bKk7wVy2Jx/fFOrqtju0vM9RnueI\nGW81WfnAv13A5xaR6s1J38Kc9K3cNaoPwxITqv8BkVpwPSYkIiFs8w9FPPTBKoZ3b82tZ2p1bKl7\nCiERqVKZx8vtM9IxRqtjS/C4XjtORELUM3PXszQnn79fOZQurbU6tgSHPtqIyBG+Xp/LS19lcsXw\nrlwwqJPrciSCKYRE5DA7C4r53cyl9G3fgkcu7O+6HIlwCiEROaTc4+X2d9MpKvXw4lVDadIw1nVJ\nEuE0JiQihzz/eQbf+5fl6XVCC9flSBTQkZCIALAgYzcvfLmBccO6cOmwLq7LkSihEBIRdhUUc+fM\ndHq1a86jv9A4kNQfdceJRDmP13LHjKXsKynnnZtOoWlDvS1I/dHeJhLlnv88g/9l/cCTlw6iT3uN\nA0n9UnecSBT7ct0unv8ig7HJnbksReNAUv8UQiJRKieviDtn+M4Heuzigbo8gzihEBKJQsVlHn4z\nPQ2vtUy+ZpjOBxJnNCYkEoUe+XAVK7cW8Np1KbpMtzilIyGRKDNj0WZmpubw2zN7cXa/9q7LkSin\nEBKJIiu27OXhj1Zxeu+23DW6j+tyRBRCItFiz/5Sbn47jXbNG/HcFUOJjdFEBHFPY0IiUaDcf4G6\n3MISZt08koRmDV2XJAIohESiwhOfrmV+xm4mXTqQwV3jXZcjcoi640Qi3HtpW3h1QTbjT+3O5cO7\nuS5H5DAKIZEIlr55Dw/MWcGpPdvwxzH9XJcjcgSFkEiE2rG3mIlvpdGhZWNevDKZuFj9uUvo0V4p\nEoGKyzxMfCuVfSXlTL02hdaaiCAhShMTRCKMtZYH3l/Bsi17eeXqYfTtoJWxJXTpSEgkwkydn8Wc\n9K3cNaoP5w7o4LockaNSCIlEkM9W7eDxT9dy/sAO3HZWL9fliFRLISQSIZZvyeeOGekM6hLP05cN\nIUYrIkgYUAiJRIAte4q48R+ptGnWiFevTdGlGSRsKIREwlxBcRk3vpFKcamHN64fTrsWjVyXJFJj\nmh0nEsbKPF5unb6EzNx9vHH9CHq310w4CS8KIZEwZa3lkY9WMT9jN0+MHchpvdu6Lkmk1tQdJxKm\npnyTxTvfb+Y3P+vJFSO0JpyEJ4WQSBj6IH0rj3+6ljEDO3LvOX1dlyNyzBRCImHmm/W53DNrGSf3\nSODpXw7WVGwJawohkTCyLCefm99Oo3f7Fky9LoXGcZqKLeFNISQSJrJy93H9G4tJaNaQf1w/nJaN\n41yXJHLcFEIiYWBXQTHXvr4IgDdvGMEJLRs7rkikbiiEREJcQXEZ101bTN7+UqaNH05Su+auSxKp\nMwohkRBWXOZhwpupZOws5JWrhzG4a7zrkkTqlE5WFQlRpeVefvN2Gt9n5/Hs5UP4aZ92rksSqXM6\nEhIJQeUeL3fOTOfLdbk8dvFAfjGks+uSRIJCISQSYrxey32zl/PJih08OKYfV56s1RAkcimEREKI\ntZaHPlzJ++lbuXt0H359epLrkkSCSiEkEiKstfzlkzVM/34zN5/Rk9/qyqgSBRRCIiHi2XkZTJ2f\nzXUjE7n/3L4Yo+V4JPIphERCwPOfZ/Dc5xlcNqwLj1zYXwEkUUNTtEUce25eBs/MW8/YoZ154tJB\nWpBUoopCSMShZ+au57nPM7g0uQtPjhtErAJIooxCSMQBay3PzMvgeX8X3BOXKoAkOimEROqZtZa/\nzV3PC19s4PKUrjw+dqC64CRqKYRE6pG1lr/+dx0vfpnJFcO78pdLFEAS3RRCIvXEWstjH6/h1QXZ\n/GpENx67eIACSKKeQkikHni8lj/OWcGMxTlcNzKRRy7srwASQSEkEnRlHi93zVzKv5dv57dn9uLu\nc/roPCARP4WQSBAVl3m4ZfoSvli7iwfOO5GJZ/R0XZJISFEIiQTJvpJybnxjMYs25vHYJQO46uRE\n1yWJhByFkEgQ7Nlfyvg3FrNy616evXyIrgckEoBCSKSObdlTxHWvLyJnzwFeuXoYo09q77okkZCl\nEBKpQ6u3FTB+2iIOlHl464YRnJzUxnVJIiFNISRSR77dsJsJb6XRonEDZt98Kn07tHBdkkjIUwiJ\n1IGPlm3j7n8upUfbZrxx/Qg6xTdxXZJIWFAIiRynV+dn8eeP1zCiewJTr02hVdM41yWJhA2FkMgx\n8ngtf/54NdMWbuS8AR145vIhNI6LdV2WSFhRCIkcg30l5dzxbjqfr93F+FO789AFJ+lSDCLHQCEk\nUkvb8g9wwxuLydi1j//7RX+uGdnddUkiYUshJFILy7fkc+M/Uiku9fD6+OGc0aed65JEwppCSKSG\n/rNyO3fOXEqbZo2YfsvJ9GmvKdgix0shJFINay0vfZXJU5+tY2i3eKZck0K7Fo1clyUSERRCIkdR\nVFrOvbOX8/Hy7Vw4uBNPjRukGXAidUghJBJATl4RE95KY+2OAn5/3olM/GmSrgMkUsfqJYSMMeOA\nfGvtvACPJwOXAz8AbYDJ1tqsap5zFJDs/7YNkGmtnVLbbUSq8u2G3dz6zhI8Xsu08cP5Wd8TXJck\nEpGCHkL+IJgKXBbg8XHAaGvtxAr3TQYmVrW9//FkAGvtkxWfxxgz4WDI1GQbkcqstUxbuJHHPllD\nUttmTLk2hR5tm7kuSyRixQTriY0xSf4wSQLyAmwTD0ytFEDjgFHVPP3EykdV1trZwOhabiNyyIFS\nD3fPWsaj/17NWSeewJxbf6IAEgmyoB0J+bvTJgIYY+4PsNkDwGFHJdba2caYJdU8fYoxJqmKLrv4\nWm4jAkBW7j5umb6EdTsLuXNUb24/qzcxWgFBJOiCdiRUQxOAxZXvrG48CJgHzD3Y5QaHjqBm1XKb\nKhljJhhjUo0xqbm5udVtLmHu0xXbuejvC9lRUMy08cO5c1QfBZBIPTHW2uD/J8ZkUkX3mDHGAsOA\nFHxddglAnr/brLrnTMM36eB+IAsOdbfVapvqpKSk2NTU1Nr8iISJMo+XJz5dy2sLshncNZ6Xrkqm\nsy7BIHLcjDFp1tqUmmzrbIq2MSbJ/2VSxYkCxphJxpiE6iYPWGuHGWPmApOAJcDZx7KNRKcde4v5\n7TtLSN20h2tHJvLHMf1o1EDn/4jUN5fdcYHGZmbiC42jMsZMwHeEMxrf5Ie0CsFW420k+ny5bhdj\nnp/P6u0FPHfFEB79xQAFkIgj1R4J1fJNO89am1/DbbMq/QuAtXaJMSY+wKSCgzVNBiYdfNwY0wPf\nWE8a0Lqm20h0KSn38OR/1vHagmxO7NCCv185lF4naP03EZeOGkL+AKr2qKSCxcCT1W4FWGvz/Wef\nBwqtJCoFlL+meP/PH3rMH3yjjTFp/vOSUqvbJtCJsxKZsnL3cfuMdFZuLeC6kYk8cH4/Lb8jEgKO\nGkL+N/EqTzKtI1kECJsA94FvEkNmgMcm4+vmq8k2EiXeS9vCQx+upGGDGKZcM4xz+ndwXZKI+Lme\noj2ZH5fVAQ6tdJB/lGnaWUDPAI/F45uAUJNtJMIVFpdx18yl3D1rGQM7t+LTO05XAImEmPoKoQSq\nPvqYwpHL80wCbjr4jX98KNM/yYAKYzyHrarg7zpsY63Nqsk2x/l6JMQtzcnnghcW8OHSrfxudB/e\nuekUOrbS9GuRUBO0Kdr+sZsH8HW3xQOTjDGjgbkHz9XxjwuN9k8iyMR39DKpivGahIrfWGsn+k8o\nnYRv0VPwHT3dX5ttJPKUeby88HkGL36VSfsWjZg5cSTDuydU/4Mi4kS9nKwaznSyavjI2FnIXf9c\nysqtBVya3IVHLjqJlo3jXJclEnXC4mRVkbri9VpeX5jNk5+to3mjBrxy9TDOHaCxH5FwoBCSsJaT\nV8Q9s5bxfXYeo/q15/GxA3XpbZEwohCSsGStZVbaFh7912oAnho3iHHDuujKpyJhRiEkYWfLniL+\nMGcl36zP5ZSkBP562WC6tG7quiwROQYKIQkbXq/l7e83MenTtVjgTxf155pTEnXZBZEwphCSsJCZ\nu4/fv7ecxRv3cHrvtjw+dqCOfkQigEJIQlqZx8vU+Vk8Oy+Dxg1iNPYjEmEUQhKyVm3by32zl7Nq\nWwHn9u/Aoxf354QWjV2XJSJ1SCEkIWdfSTnPzF3PtIXZJDRrxMtXJXPewI6uyxKRIFAISciw1vKf\nlTv4079Ws6OgmF+N6Mb95/YlvmlD16WJSJAohCQk5OQV8fCHK/lyXS79OrbkpauTSe6maw+KRDqF\nkDhVWu6bePDCFxnEGsODY/ox/tTuNIh1fZUREakPCiFx5tsNu3n4o1Vs2LWP8wZ04OELT9LlFkSi\njEJI6l1OXhF/+WQNn67cQdeEJkwbP5wzTzzBdVki4oBCSOrNgVIPL3+dyeSvM4kxhnvO6cOvT0+i\ncVys69JExBGFkASdtZaPV2znLx+vYdveYi4a3IkHzj9RXW8iohCS4FqzvYD/99Eqvs/Oo1/Hljx7\nxVBG9NCVTkXERyEkQbGroJin/7ueWWk5tGoSx2OXDOCK4d2I1WKjIlKBQkjqVFFpOVO+yWLy11mU\ne71c/5Me3HZWL51wKiJVUghJnfB4LbPTcnj6v+vZVVjC+QM7cP+5J5LYppnr0kQkhCmE5Lh9vT6X\nxz9Zw9odhQztFs/LVyczLFHjPiJSPYWQHLPlW/J56rN1zM/YTbeEprx4ZTLnD+ygyyyISI0phKTW\nNuwq5On/rufTlTto3TSOB8f045qRiTRqoPN9RKR2FEJSY1v2FPHsvAzeX7KFJnGx3HF2b359eg9a\nNI5zXZqIhCmFkFQrt7CEF7/cwPTvN2GM4Yaf9OCWM3uR0Ewz3kTk+CiEJKC9B8qY+k0Wry/MpqTc\nyy9TunDbWb3pFK+VDkSkbiiE5Ah7D5TxxsKNvLYgi4Lici4c3Im7RvUmqV1z16WJSIRRCMkhew+U\n8fqCbF5fmE1hcTmjT2rPHWf3ZkDnVq5LE5EIpRAS9haV8drCbKb5w+eck9pzu8JHROqBQiiK5ReV\n8tqCbN5YuJHCknLO7d+B287uRf9OCh8RqR8KoSj0w74Spi3cyBvfbmRfSTnnDejA7Wf3pl/Hlq5L\nE5EooxCKIlvzDzD1myxmLN5McZmX8wf6wufEDgofEXFDIRQFMnYW8vLXmXy0dBsAFw/tzM1nJNHr\nhBaOKxORaKcQimDpm/fw0leZzF29kyZxsVwzMpGbTk/SeT4iEjIUQhHGWsv8jN289NUGvsvKo1WT\nOG4/uzfjT+2uFQ5EJOQohCJESbmHfy3bzqvzs1i7o5D2LRvx4Jh+/GpEN5o10q9ZREKT3p3C3J79\npUz/fhP/+N8mcgtL6Nu+BU+OG8QvhnTSqtYiEvIUQmEqK3cfry/MZnbaForLvPy0Tzv+9ssenNar\nra7nIyJhQyEURqy1fJ+dx6vzs/h87S7iYmK4ZGhnbjitB307aKabiIQfhVAYOFDq4aNlW/nHt5tY\nvb2AhGYNue2s3lxzSiLtWjRyXZ6IyDFTCIWwnLwi3vpuEzMX57D3QBkndmjB42MHcsnQzjSO03iP\niIQ/hVCI8XotCzbs5s3/beTztbuIMYZz+3fg2pGJjOiRoPEeEYkoCqEQUVhcxntpW3jzf5vI2r2f\nts0bctuZvbjy5EQ6tGrsujwRkaBQCDm2dkcB07/bzPtLtrC/1MPQbvE8e/kQzhvYQVOsRSTiKYQc\nOFDq4d/Lt/HOos2kb86nYYMYLhzUietOTWRQl3jX5YmI1BuFUD1au6OAd77fzJz0rRQWl9OzXTMe\nHNOPS5O70FpL6ohIFFIIBVlRaTn/Xr6ddysc9YwZ2JFfjejG8O6tNdFARKKaQihI1mwv4N1Fm5mz\nZCuFJeX0OqE5D11wEmOHdtZRj4iIn0IoCDxey/hpi9hTVMaYgR258uRupCTqqEdEpDKFUBDExhhe\numoYPds1I76pjnpERAJRCAXJsMTWrksQEQl5Ma4LEBGR6KUQEhERZxRCIiLijEJIREScUQiJiIgz\nCiEREXFGISQiIs4Ya63rGkKaMSYX2HSMP94W2F2H5UQ6tVftqc1qR+1Ve8fSZonW2nY12VAhFETG\nmFRrbYrrOsKF2qv21Ga1o/aqvWC3mbrjRETEGYWQiIg4oxAKrimuCwgzaq/aU5vVjtqr9oLaZhoT\nEhERZ3QkJCIiziiERETEGV1PqI4ZY5KBFCALSAKyrLXz3FYVGowx8cAEYDaQByQAE4G5Fdso2tvQ\nGDMOyK/qNdekbaKx/QK1mfa5I/nbKgno6f93srV2dqVt6m8/s9bqVkc3/y9ibqX7ZgFJrmsLhZu/\nfWyF2x5gnNrwsNc6yt8uowK031HbJhrbrwZtpn3ux9c1Dkiu8H08kAlMcLWfqTuubk0EJle6bzIw\nyUEtoWo00Broaa1tbSt9AiNK29AYk2SMmYzvjzsvwGY1aZuoab8athlon6soyVq75OA31tp8fK+z\n4uuv3/3MdTJH0g3fJ4qkSvfF+5rZfX2ub/jeLI74tKo2rLINqvpUX23bRGv7HaXNtM8d/prSgPgq\n2sgebIP63s90JFRH/H3PR3wis75PGhhjklzUFU7UhoHVpG3UfrUXTW3mf01J/luVXOxnmphQdxLg\nx19EFZLwDeBFuyT/wCj42izP/tg9ojYMrCZtQw22icb20z7nZ61tXcXdo/BN6sg6GCD1uZ8phOpO\nvOsCwkAeQIU3AIwxs4wxB+9TGwZWk7ZR+x1J+1z1JgKP+7+u9/1M3XFSb6y1+dbaykuAROoAsIQA\n7XNHZ4yZgO/I8ElXNSiE6pi/v1RqLgtfd8mhdlMbBlaTtlH7VUv7HIfGbiZaa0dX8Vi97WcKobpz\nsA80oeKdFX5REdOvfKyMMfdVcffBwc2K/chqwyPVpG3UfpVonzuqScDZle6r9/1MIVRH/IN0WRzZ\nX5qAf9Cv/qsKHf5PXZOqmDlzcEfOUhsGVpO2UfsdTvtcYP7zq+6vPLnAxX6mEKpb8/AtY1FRsv/+\nqObfMSdWsYOOApZU+GNQGwZWk7ZR+/lpn6uafxxoUsV2McaMqhDW9bufuT6BKpJu+D4ZVF7KYi4R\nuPzHMbbPuIptwY8nz1VeRiSq25AqlpapadtEa/sdpc20zx3ZHhP8gXHwNgrf+nFO9jNdT6iO+Rf1\nuxxYjK/PeYmN0IUQj0WFxRPb4NuRD/tE5t8m6trQ35/+AL7XOw5fd8c8fH/oFacXV9s20dJ+tWgz\n7XMcaq89AR7Ostb2rLBtve1nCiEREXFGY0IiIuKMQkhERJxRCImIiDMKIRERcUYhJCIiziiERETE\nGYWQiIg4oxASERFnFEIiIuKMQkhERJz5/+HIRIpe6XmpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119f745f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Implementation \n",
    "v  = Variable(torch.DoubleTensor([3.3]), requires_grad=True)\n",
    "m0 = Variable(torch.DoubleTensor([0.]), requires_grad=False)\n",
    "m  = Variable(torch.DoubleTensor([3.]), requires_grad=False)\n",
    "v0 = Variable(torch.DoubleTensor([40.]), requires_grad=False)\n",
    "c  = Variable(torch.DoubleTensor([0.006,0.994]), requires_grad=False)\n",
    "\n",
    "def loglik(X, m, v, m0, v0, c):\n",
    "    return torch.log(c[0]*torch.exp(log_gauss_pdf(X,m,v)) \n",
    "                     + c[1]*torch.exp(log_gauss_pdf(X,m0,v0)) )\n",
    "\n",
    "def log_gauss_pdf(x,mu,v):\n",
    "#    return -torch.log(v) - (x-mu)*(x-mu)/v\n",
    "    return -0.5*torch.log(2*np.pi*v) - 0.5*torch.pow(x-mu,2)/v\n",
    "\n",
    "\n",
    "eta = 2\n",
    "MAX_ITER = 200\n",
    "V = []\n",
    "L = []\n",
    "for epoch in range(MAX_ITER): \n",
    "\n",
    "    # Compute the loglikelihood\n",
    "    LL = torch.sum(loglik(X, m, v, m0, v0, c)) \n",
    "    L.append(LL.data)\n",
    "    # Compute the gradients by automated differentiation\n",
    "    LL.backward()\n",
    "    \n",
    "    # The gradient ascent step\n",
    "    v.data.add_(eta*v.grad.data)\n",
    "    V.append(1*v.data.numpy())\n",
    "    \n",
    "    #print(v.grad.data)\n",
    "    # Reset the gradients, as otherwise they are accumulated in v.grad\n",
    "    v.grad.zero_() \n",
    "    \n",
    "print(v.data.numpy())\n",
    "\n",
    "plt.plot(V)\n",
    "plt.show()\n",
    "plt.plot(L)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "1.00000e-03 *\n",
      "  2.2835\n",
      "[torch.DoubleTensor of size 1]\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0\n",
       "[torch.DoubleTensor of size 1]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "v  = Variable(torch.DoubleTensor([0.05]), requires_grad=True)\n",
    "m0 = Variable(torch.DoubleTensor([0.]),  requires_grad=False)\n",
    "m  = Variable(torch.DoubleTensor([3.]),  requires_grad=False)\n",
    "v0 = Variable(torch.DoubleTensor([40.]), requires_grad=False)\n",
    "c  = Variable(torch.DoubleTensor([0.006,0.994]), requires_grad=False)\n",
    "\n",
    "def loglik(X, m, v, m0, v0, c):\n",
    "    return c[0]*torch.exp(log_gauss_pdf(X,m,v)) + c[1]*torch.exp(log_gauss_pdf(X,m0,v0)) \n",
    "\n",
    "LL = torch.sum(loglik(X, m, v, m0, v0, c)) \n",
    "\n",
    "LL.backward()\n",
    "\n",
    "print(v.grad.data)\n",
    "\n",
    "v.grad.zero_()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEECAYAAAAh5uNxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4VOXZ+PHvyUYWAlkhJASysATC\nFgj7LogIslgWBUSpC/RV31frRm3706ptfYFWX5VaoVbUiigiClhAQWSVLQmENUBICCQQQkjIStZ5\nfn+cSYnIksDMnFnuz3XlypwhmXOTSc59nu1+NKUUQgghXI+b0QEIIYQwhiQAIYRwUZIAhBDCRUkC\nEEIIFyUJQAghXJQkACGEcFGSAIQQwkVJAhBCCBclCUAIIVyUh9EB3EhISIiKiooyOgwhhHAoycnJ\n+Uqp0Jt9nV0ngKioKJKSkowOQwghHIqmaVkN+TrpAhJCCBclCUAIIVyUJAAhhHBRkgCEEMJFSQIQ\nQggXJQlACCFclCQAIYRwUZIAhBDiet5+G95//8rxgAHw6KNXjl94AT76yPZxWYgkACGEqPPxx/Dq\nq1eOv/oK1q+/cjx6NPTvf+X4++/h8GH9sVLw4ouQmmqbWC3ArlcCCyGE1Z04Ae3b649TUmDLFvj9\n78HNDTZsAI96l8mXXvrp9yYng8mkP87MhIUL9dfq3l1PCACaZv3/wy2SFoAQwnUtXQodOsDBg/rx\nggWwb59+8YefXvyvp+5rY2IgJwemT9ePP/8cBg2CvDzLx20hkgCEEK7l3Dk4elR/fM898PrrEBmp\nH3t63t5rN2sG3t76Y3d3/Tg4+PZe04okAQghXEdtLQwZAk89pR83bw6/+Q0EBFj+XFOmwLp1eiIo\nLYXJkyEtzfLnuQ0yBiCEcH65uRAWpl+M330XbF1m/sQJ2LFDb33Exdn23DcgLQAhhHNLTdUHZj/5\nRD++884rg762kpAAGRkwfLh+nJZ2ZZDYQJIAhBDOrUsXmDMHhg41Ng4fH/3zgQP6LKF33zU2HiQB\nCCGcUXo6TJoExcV6t89f/nJloNdoXbroaw3qZgsZSBKAEML5ZGfrfe4nThgdyc+5ucHcuRAYqA9K\nL1tmWHeQJAAhhPM4flz/PGyY3ufeq5eh4dzUsmV6S2DTJkNOLwlACOEcFi/Wu1f27dOPfX2Njach\nZsyAb7+FESMMOb1MAxVCOIepU+HiRejWzehIGk7TYNQo/XFmJpw+bdPBamkBCCEcV2mpvpK3tlZf\nzPXii/qgryN67DH45S+hutpmp5QWgBDCca1aBf/v/+mrewcONDqa2/P++3oiu91yFI0gLQAhhOOa\nMUMv5OboF3/QVyfHxuqPt261ycwgSQBCCMeSl6ev5q2b4tmpk7HxWNqWLXpVURuUkZYuICGEY8nP\n1y/+Z8/avqSDLdhwEFgSgBDCMVRVgZcXdO6sz/f38jI6IocnXUBCCPt36RL07Xulfo5c/C1CEoAQ\nwv75+OhllNu1MzoSpyJdQEII+3X5sj4bxtdXL5sgLEpaAEII+6SUvrp37Fh9frywOGkBCCHsk6bB\nAw9cKeksLE4SgBDCviil18Rp2xbuu8/oaJyadAEJIezLO+/oVT2PHTM6EqcnLQAhhH2ZPBkKC51z\nkZedkRaAEMI+1G2UHh4OL7+s75wlrEp+wkII4x0/Dj17woIFRkfiUqQLSAhhvPbt4bXXYOZMoyNx\nKZIAhBDGOX9e/9yyJTz7rLGxuCBJAEIIY9Qt9Cos1Pfxlbn+NicJQAhhDE2D+fP1fXzl4m8ISQBC\nCNtSCg4cgO7d9QqfwjAyC0gIYVsffqjP+PnxR6MjcXnSAhBC2NbUqVBUBP36GR2Jy7N6C0DTtABN\n017QNC3G/DhG07R5mqaNtPa5hRB25PBhfVcvPz94+mlZ6GUHbPEOBAHzgJNAIZAM7FVKbbTBuYUQ\n9qCwEIYMgSefNDoSUY+tuoDuBJKAIKVUho3OKYSwF4GBsHChDPraGZuNASilLgGXbHU+cW0mk+JC\naSVnL10mt6iCyhoTmgZumkaAryetA30JD/CmiYdMyxO3TinFxbIqcs5epPjYSQrbxKI6DUEVQfCJ\nC7Rq7kN4gDe+XjIMaST56Ts5pRSp2UXsSM9nd2YByacKKKu6+e5K0SF+9GwTSGJUIANig2kb7GeD\naIWjqqiuZVfGRZJOFZKUVcCB7CLKq2r58/p3uPvYjwyd8w+KvZv+7PvatWhK3+gg+sYEM7RDKM19\nPA2I3nVpSinrnkDTYoCRQIH5qSCgQCm14mbfm5iYqJKSkqwZntPKK6lgZUoOy5POkHGhDICOLf3p\nEx1EhzB/wpt706q5Dz5e7piUwmTS79iyCy+TXVjO4bPFJGcVUlBWBUCnVs0Y0yWMe7qHEx0iyUDo\nF/1NaXmsPXiOTWl5lFfV4u6mER/ejITIAKJC/Gh3uYCooylU3nc/bpqGAvJLKjlbdJnsgsskny4k\n6VQhpZU1NPFwY3SXMO5LjKRfTDBubprR/0WHpWlaslIq8aZfZ4MEEABMVUotrvfcF8Dn10oCmqbN\nBmYDtGnTpldWVpZV43M2F0oqWbjpBJ/uOU11raJ3VCBTEiMZ2aklQX5ejXotpRSZ+WX8cOwC6w6e\nIymrEIBB7UJ4oF9bRnZqgYe7zORwNacvlrN0dxbLk85QWF5NSFMvRsWHMTo+jMSoQL1bJyUFEhL0\n1b43UVNr4kBOEV+l5PD1/hxKKmqIC/PnhdEdGd6xBVoDXkP8lN0kgGueVJ8CukgpFXujr5MWQMNd\nrqrl3c3p/HN7JpU1Ju7rHcmjg6KJCf15s/tW5RZV8GVKNkt3ZXG2qILw5t78algsUxMj8faUMQNn\ndyiniIWb0vn2SC5umsZd8S2Z3qct/WODca9/t75jBwweDB98ALNmNeocFdW1/PvAOd7edIKsi+X0\njgrkd2M70yMywLL/GSdn7wkgBn1aaKB5cPiaJAE0THJWIc99kUpmfhnjuofzzJ0drNpNU1NrYlNa\nHou3ZpCUVUgL/yb8amgsM/q1kcFjJ3Qop4g3Nxzn+7Q8/L09eKh/FA/0a0tYc+9rf4PJBH/7Gzz6\nKPj43NI5q2tNLE86w1sbT5BfWsmcobE8PbK9/H41kN0kAE3TXlBKzb/quQD0NQG9lFIp1/teSQA3\nVlVj4s2Nx1m05SStmvuwYEo3BsSG2Oz8Sil2Zlzk7e9PsCujgMggH+aOjmNs11bSbHcCZy9d5i/f\nHmPlvhya+3jy6KBoHhwQdf2B2vR0fbpncLDFYiipqOaP3xzl86QzdGzpz5v39aBzeDOLvb6zsosE\nUO9OP7b+/H9pAdy+wrIqfvVJMrszC7i/dyS/G9sJf2/jZlBsPX6BP689SlpuCQltAnh1fBe6tm5u\nWDzi1lVU1/LuD+ks2pqBAh4eGM3jw2NpdqPfr5oa6NpVr+v/ww8N6vtvjE1p55n75UFKKqr565Qe\njO3WyqKv72zsIgGYA5ldfwC47jlgjlKq142+VxLAtZ04X8IjHyWRW1zB/EndmJgQYXRIANSaFF8m\nZzP/22MUlFUys19bnhnVUab2OZAf0vJ4afUhzhRcZnz3cF4Y3ZHWgb4N++Zt2/Qun8SbXnduSX5p\nJXP+lUxyViG/HtmB/xnRTlqa12FPCWAykFLXAjB3/3wPPHaj7h+QBHAtO09eZPbHSTTxdOcfD/Yi\noU2g0SH9TNHlat7ccJyPd54iyK8Jr06IZ0xXuWOzZ3nFFby8+jDrDuUSG+rHaxO7NKw7saoKkpOh\nf3/rBwlU1tTy4sqDrEzJYUKPcP46pbvMRLsGu0kA5mAmAzFAMBAAzGtISQhJAD+1Iz2fRz7aS2Sg\nLx893IfwgFsbYLOVQzlFvLjyIAdzihgdH8arE+Jp0ew6A4fCEEopViRn89o3R6ioMfHUiPY8NjgG\nL48GXlRfeglefx2OHYOYGOsGa6aU4t3NJ1nw7THGdA3jrfsT8JQk8BN2lQBulSSAK7aduMCjHyUR\nFezH0sf6EtK0idEhNUhNrYn3t2fyxobjeHu48cqEeCb2iJCmux3ILapg7pcH2HL8Ar2jAvnfSd2I\nbey04eJiWLsW7r/fOkHewPvbMvjjv48yOj6Mt6clNDxpuQBJAE7kx/R8Zn24l5gQP5Y+2pdgB7n4\n15dxoZTnVxwgOauQ0fFh/OneLg75/3AGSilW7T/LS6sOUV2rmDu6Iw/2j2rcytu1a+HOO8HT2PGd\nD7Zn8uo3R7grviXvzuj10/UILqyhCUBSpp07llvCnH8lExXsy6eP9XPYi2ZMaFOWz+nPb+6OY1Na\nHqPe3MrGI+eNDsvlFJZV8cSnKTz9+X7atWjK2qcGM2tgdOMu/qmpMHYsvPOO9QJtoIcHRfPSPZ35\n9vB5XllzGHu+obVHUgzOjp0vruCXS/bg4+XOh7/s0+hSDvbG3U3jV0NjGdYxlF9/nsqjHycxvW8b\nfj+2k1SFtIFtJy7w3BepFJRV8cLojswZEntrd8zdu8OqVTB6tOWDvAUPD4omt7iCxVsziAz05bEh\nthmLcAbSArBTZZU1PPzhXoouV/PBrN52P+DbGHFhzfj6iQHMGRLDsj2nueft7RzMLjI6LKdVUV3L\nq2uOMPOfe/D39uSrxwfy+LB2jb/4b94MaWn64/Hjwct+bkh+MzqOMV3D+NPao6w9eM7ocByGJAA7\npJTiuS9SScstYeGMnnSJcL4FVU083HlxTCeWPtqXy9W13PvuDv6++SS1JmnCW9Lx8yVM/NsOPtiR\nyYP927LmyUG39vtUXQ2PPAJPPGH5IC3AzU3jjak96NkmgF9/vp8jZ4uNDskhSAKwQx/sOMW6Q7nM\nNVdDdGYDYkNY99RgRsW3ZN76NGa8v4tzRZeNDsvhKaX4eOcpxr2znfzSSpbM6s2rE7rg43WLtXQ8\nPWH9evjkE4vGaUnenu4smplIcx9PHl+aTHFFtdEh2T1JAHYmOauA19ceZVTnljw22DX6MgN8vfjb\n9J4smNyNA9lFjP6/bayTZvwtyy+t5JGPknhp1WH6xwaz7qkhDI+7xRuJtDR47z39cfv20Mq+F/SF\n+jdh4fSenCm8zNwVB2RQ+CYkAdiRi6WVPLF0HxGBPiyY0t2l5sprmsaUxEjW/s9gooJ9+a+lKcxd\ncYCyyhqjQ3Mom4/lMfr/trE9PZ8/jOvMklm9CfW/jZlj77wDf/gDFBTc9EvtRZ/oIH4zOo51h3L5\nYMcpo8Oxa5IA7IRSiudXHKCgvIp3Z/R02fo5USF+rPivATwxPJblyWcY+/Y29p0uNDosu1dRXcvL\nqw4xa8legv28WPPkIGYNjL79m4i33oKdOyEoyDKB2sijg6O5K74lr689yqEcmWBwPZIA7MTne8+w\nKS2PF++OIz7c+QZ9G8PT3Y3n74rjs8f6UV2rmPzeTt7ccJyaWpPRodmlQzlF3PPOdj7amcXDA6NZ\n9eRAOob53/oLnj4N06frq3w9PCA62nLB2oimacyf1J3gpl48s3w/FdU33wfbFUkCsANnCsp57Zsj\n9I8J5qH+UUaHYzf6xgSz7unBjO8ezlvfn2DS338kPa/U6LDsRk2tiYWbTnDvuzsoqajmX4/04aVx\nnW9/d7aDB2HDBsjMtEygBmnu68m8Sd04fr6UNzccNzocuyQJwGAmkz7lU9M0FkzpJhthX6WZtydv\n3teDhdMTyCooZ+zb2/jn9kxMLj5dND2vlEl//5G/fHecu+LDWP/UEAa3D729F6013yWPHQsZGfqC\nLwc3rGMLpvVpw+JtGew95TjjGLYiCcBgS348xe7MAl4a17nhdddd0D3dwvnu10MY1C6E1745wv2L\nd5GZX2Z0WDZXU2ti0ZaTjH17G1kF5bwzLYGF03sSeLurxDMz9Qv+5s36sf9tdCHZmd+N7UTrQB+e\nXZ5KeZVMKqhPEoCBci5d5q/fHeOOuBZM6dXa6HDsXgt/b95/KJH5k7uRllvM6P/byqItJ11mbCAt\nt5hf/P1HXl+XxuD2oXz39BDGdQ+3zIv7+enbOTZtZDVQB9C0iQd/mdyd0wXlvP19utHh2BVJAAZ6\nZfVhTErxyvh4l5ryeTs0TWNqYiQbnhnK0A6hvL4ujYnv7iD1zHV3FnV4l6tqmb8+jXHvbCen8DLv\nTEvgHw/2sszeCqdP65u4t2gBW7dabTcvo/WNCWZqYmve35bBsdwSo8OxG5IADLLhyHm+O3Kep0Z0\nIDJIun4aq2UzbxbN7MXfpvckr7iSie/u4PdfH6TosnOt/vz+6HnufHML724+ybju4Wx4Zijjuodb\n5oYhOxt69IBXXtGPnfwm5Dd3d8Lf24Pff33Q5ceQ6kgCMEB5VQ1/WH2YDi2b8uhgx5tiZy80TWNs\nt1Z8/+xQZg2I4tPdpxn+l838a1eWw3cLpeeV8Msle3jkoyR8PN35fHY/3pjaw7IVYSMi4Pnn4aGH\nLPeadizIz4sXx3Ri76lCViRnGx2OXZANYQzwv+vSeG/LSZbP6U+faMdaYGPPDuUU8eo3R9iTWUCH\nlk357ZhODO0Q6lDdaxdLK3nr+xMs3X0aX093/ntEO2YNiLbsblcrV+pdPW3aWO41HYTJpLhv8U5O\n5JWy6dlhDl9i/XpkQxg7lXWxjA+2ZzKpZ2u5+FtYl4jmfD67H+890JOKahOzluxl6qKd7Dx50ejQ\nbupSeRXz16cxeP4PfLIri2l9Itn8/DBmD4m17MW/sBAeffRKt4+LcXPT+OPErpRU1PDWRlkbILtw\n2Nj/rkvD3U3jhdEdjQ7FKWmaxugurRge14Lle8+w8Id0pv1jF/1jgpkzNMbuWgR5xRV8+OMp/rUz\ni9KqGu7pFs5TI9rTroWFZ+MopffxBwbCpk0QF2fZ13cgHcP8mdYnkk92n2Zm/yjL/6wdiHQB2dCe\nzAKmLtrJr0d24KmR7Y0OxyVUVNeydPdpFm89yfniSjq0bMojg6K5p1s4fk2Mu/85lFPEhz+eYtX+\nHGpMitHxYTw1sj1xYc0sf7LycnjgAbj3Xpg50/Kv74AullYybMFm+kQH8c9ZvY0Ox+Ia2gUkLQAb\nMZkUf/z3EcKaeTNbtqyzGW9Pdx4ZFM3Mfm1Zk3qWf2zLYO6XB3l1zRHGdQ9nSmIkCZEBNlmBXVBW\nxer9OSxPyubIuWJ8PN2Z3qcNDw+Kpm2wn/VO7OkJJSVQJEXR6gQ3bcKTd7Tj9XVpbD+Rz6D2IUaH\nZAhpAdjIypRsnlmeyhtTu/OLnrLoyyhKKfaeKmR50hn+feAcl6trCWvmzZ2dW3Jn55b0jgq69U1T\nrnGuUxfL+SEtj++O5LInswCTgi4RzZiaGMmE7hE097Vi1dd9+6BDB32Rl8kEbjLkV19lTS0j39iC\nn5cH//6fwbe2P7KdamgLQBKADVTW1DJ8wWZC/Jvw9eMDpd6PnSipqOa7w+fZcOQ8W45f4HJ1LZ7u\nGl0imtM7Koi4MH/atWhKbGjTm3YXmUyKnEuXSc8rJT2vlP1nLrHnVAEXSioB6NjSn1HxLbm7Sys6\nh1uhm+dqFy7oVTxnzYKFC61/Pge19uA5Hl+awvzJ3ZiaGGl0OBYjXUB25NPdpzlbVMGCKd3l4m9H\n/L09mdSrNZN6taaiupadGRfZk1nA3swCPtxxiqp6awmaNvEgyM+LQD8vvNz191ApKKmo4WJZFYXl\nVT/ZzzgiwIdB7ULoHRXEgNhgokKs2MVTX91gb2goLFkCw4bZ5rwO6u4uYXRr3Zy3Np5gQo9wmnhY\npvXnKKQFYGXlVTUMmf8D7Vv4s2x2P6PDEQ1UXWsi62I56XmlZOSXcqGkkoKyKgrKfnqh9/fWE0OQ\nnxetA31p36Ip7Vo0JcDXgPnl587pg71//jP07Wv78zuorccv8OAHe3h1QjwPOkk5dmkB2IklO06R\nX1rFopky7dOReLq70c58MXcYTZrA+fN6IhANNrh9CH2ig3hnUzpTekVabAzIEciokBUVXa5m0ZaT\n3BHXgl5tA40ORzij6mq9q8dk0rdtTE2FiRONjsqhaJrG83d15EJJJR/tPGV0ODYlCcCK3t+WQXFF\nDc+O6mB0KMJZrVoFDz8M69frx+6uc/dqSb2jghjWMZT3tpykuMK5CgreiCQAK7lUXsUH2zMZ27WV\ny+/xKyxMKb2SJ8CkSXoZ5zFjjI3JCTw3qiOXyqv5cMcpo0OxGUkAVrJkxynKqmr57xHtjA5FOJsX\nXoBevSA/X5/xM3iw0RE5hS4RzRkR14IPdmRSVukaO4dJArCCkopqluzI5K74ltZZ2i9cj1JQY74o\nzZwJL76o9/kLi3ryjnZcKq9m6e4so0OxCUkAVvDxziyKK2p4crjU+xEWUF2td/H89rf6cbdu8PTT\nsrLXChLaBDK4fQiLt2ZSUV1rdDhWJ79BFlZeVcM/t2cyrGMoXVtL37+4DXVrdDw9oVMniIoyNBxX\n8eTwduSXVvL53jNGh2J1kgAs7NPdpykoq+K/75C+f3EbUlP1fv508ybmb7wBjz9ubEwuom9MMH2i\ngnhvy0kqa5y7FSAJwIIqa2pZvDWD/jHB9Gor/bPiFtTd9YeG6nP7L1wwNh4X9eQd7ThXVMFXKTlG\nh2JVkgAsaNW+s+SVVPJfw2KNDkU4oj/96Uq9/vBwvZpn//7GxuSiBrcPoUtEMxZvy3DqDeQlAViI\nyaRYvC2DTq2aMdhFa4uLW1Bz1XRDNzd90Bf0KZ7CEJqm8djgGDIulLEpLc/ocKxGEoCFbDl+gfS8\nUmYPibarLQeFHTtyBDp2hC1b9OPf/hY+/lgf9BWGG9u1FREBPizemmF0KFYjCcBCFm09Savm3tzT\nLdzoUIQ9M5n0gm2g1+uPi7tywZcbB7vi4e7Gw4Oi2XOqgH2nC40OxyokAVjAgexL7Moo4OGB0Xi6\ny49U3MDkyTB2rD7Y6+MD//43DBhgdFTiOu7rHUkzbw/+sc05WwFytbKAxVsz8G/iwf19nGdHIWEh\ntbWwcuWVvv4HH9QXcdnxPhziiqZNPJjRry3rD+WSdbHM6HAsziYJQNO0npqmzdY0bWTdZ1uc1xZy\nLl1m3aFcpvVtg7+39N2Kq6xfrxdsW7NGP544Ud+0RVbxOoxfDojC3U3jg+2ZRodicVb/LdQ0LQaY\np5RarJTaqJRaDMwxP+/w/rUzC6UUDw2IMjoUYQ9qa2HePH0wF+Duu+Gbb2DCBGPjEresRTNvxnUL\nZ0VyNiVOViraFrchc4BFVz23CJhng3Nb1eWqWj7be5q74sOICPAxOhxhpPJy/bO7u363v3mzfuzm\npvf5yx2/Q3toQBRlVbWsSM42OhSLssVv5WQg5arnkszPO7Sv9+dwqbyaWXL379refluf0VNaqh9/\n9x188IGxMQmL6h4ZQM82AXz04ymnWhhm1QSgaVoAEAMU1H9eKXXJ/O8O2w2klOLDHafo3KoZfaKl\n7INLMZn0u/y66Zx9++r9+lVV+rGvr3GxCav55cBoTl0sZ/Nx51kYZu0WQBBcueBfg8MmgJ0ZFzl2\nvoRZA6Nk4ZeryczU+/SXLNGP+/aFv/5V6vM7udFdwmjZrAlLnGjHMGsngIDGfoN5llCSpmlJF+y4\nENaHO04R5OfF+O6y8Msl/OlP+k5cALGxeh//c88ZGpKwLU93N2b2a8u2E/mk55UYHY5F2N3IlHm2\nUKJSKjE0NNTocK7pTEE5G4+eZ1qfSLw9ZRNup5Vdb8AvJwfOnLkyf3/IEPDwMCYuYZhpfdrg5eHG\nRz86x45htloH0OiWgD1btuc0ADP6tjU4EmE1n3wCbdvCsWP68cKFsGyZlGtwccFNm3BP11Z8tS/H\nKfYNtnYCqFs//ZPO0XoJweHWV1fVmFiedIY74loSLlM/nYdS+nz9/fv141Gj4OWX9br8INM4xX/M\n6NeW0soavt7v+HsFWPW32jz4m8HPxwKCgEtKKYdLAOsP55JfWsUD/doYHYqwpPJymDUL3nlHP27R\nAl56SQZ2xc/0bBNAp1bN+GTXaZSDl/SwxW3NRiDxqud6mp93OJ/syqJNkC9D2tvn+IRohNWr4aGH\n9Lt/Pz99YPe994yOStg5TdN4oF8bjp4rJuX09SY4OgZbJIC5wJSrnptjft6hHD9fwp7MAqb3bYOb\nm/QFOySTSf8AfZB33z4oMC9T6dJFavGLBpnYI4KmTTxYusuxB4OtngDM3UBzNU2bp2naZE3TXkCv\nDeRw3T9Ld2Xh5e7GlF6tjQ5F3IozZyAhAb7+Wj+ePVvv8w8ONjYu4XD8mnhwb0IE3xw8R2FZldHh\n3DKbjGwppVKUUnOVUiuUUvOVUg7X/VNeVcPKlBzGdA0juGkTo8MRjVFinrPdqhVERup1+EGfximD\nu+IWzejXhqoaE18knzE6lFsmv/0N9E3qOUoqa5jRT6Z+OpTf/Q569tTLNHh46DN97r7b6KiEE4gL\na0Zi20A+23PGYQeDJQE00LK9p2nXoimJbQONDkXcTGnplbo8Q4fC9Ol6mWYhLOz+Pm3IyC9jd2bB\nzb/YDkkCaIBjuSXsO32J+3tHSt0fe3f+vL7Pbt10zlGj4JVXrnT7CGFBY7qG4d/Eg8/3OmY3kCSA\nBvhs72k83TV+0VMGf+1WRYX+uWVLmDZN9tkVNuHr5cGEhHDWHjxHUbnjbRYjCeAmKqpr+WpfDqPi\nwwjy8zI6HHEty5bpBdryzGV6FyyA/v2NjUm4jPt7t6GyxsRX+xxvsxhJADfx7eFcLpVXM623rPy1\nO3UDbwkJel+/EAboEtGcrhHN+Wyv4w0GSwK4ic/2nCEyyIcBsTJX3G4opZdpeOYZ/TguDj79VC/f\nIIQB7usdSVpuCanZRUaH0iiSAG7gVH4ZOzMucl9ipKz8tSeaBsXFUFR0ZVWvEAaa0CMcH093PjNX\nCnYUkgBuYEVyNm4aTO4VaXQowmTS9949elQ/fuMNfd9dWcgl7IC/tydju7XimwPnKK9ynDLR8tdz\nHbUmxZcp2QzpEEpYc2+jwxEXL8Krr17ZbF0u/MLOTOnVmtLKGtYfyjU6lAaTv6Lr2J6ez7miCqbI\n3b+x0tL0Pv/QUNi7F+bPNzoiIa6pT3QQbYN9+SLJcWYDSQK4ji+SzhDg68nIzjKwaJjt2yE+Hj77\nTD+OjpYduYTd0jSNyT1bszPJJpu/AAAQGElEQVTjImcKyo0Op0EkAVxDUXk13x05z4Tu4TTxkD1/\nDdO/P/zxjzB2rNGRCNEgk3q1RtP08UNHIAngGlan5lBVY2JKonT/2NzRozBuHFy6BO7u8OKL0KyZ\n0VEJ0SDhAT4MahfCiuRsTCb7XxMgCeAavkjOJi7Mn/hwufDYXH6+vklLZqbRkQhxSyb3ak3Opcvs\nyrhodCg3JQngKsdySziQXcTURCn8ZjNKwcGD+uPBgyE9XV/dK4QDuis+jGbeHnzhAN1AkgCu8mVK\nNh5uGhN6hBsdiuv429/0mv2pqfqxt0y7FY7L29Odcd3DWXfoHKWV9r0mQBJAPTW1Jr7al8PwuBay\n65ctzZypF3Dr2tXoSISwiEm9WlNRbWLtwXNGh3JDkgDq2Z6ez4WSSib1jDA6FOd39iw89xzU1EDz\n5vD007K4SziNhMgAokP8WJli391A8hdXz8qUHAJ8PRkeJ3P/re6772DRoiulHYRwIpqm8YuECHZl\nFJBdaL9rAiQBmJVUVPPt4VzGdZO5/1ZVbd40Y9YsOH5cun2E07rX3JPwVUqOwZFcnyQAs3UHc6ms\nMfEL6f6xnh9+gI4d9Qs/QKtWxsYjhBW1DvSlX0wQK/fl2O0+AZIAzFakZBMT6kePyACjQ3FeERH6\nzl3+/kZHIoRN/KJnazLzy0g5fcnoUK5JEgBwpqCcPZkFTOrZWub+W0NdP3+HDrBhg9z5C5dxd5cw\nvD3d7HYwWBIA8PU+vY9uYoJ0/1jc5s3QpQssX250JELYnL+3J6Pjw1iTepbKmlqjw/kZl08ASim+\n2pdDv5ggIgJ8jA7H+QwapBd0u+ceoyMRwhATEyIorqjhh7QLRofyMy6fAA5kF5GRX8a9cvdvWStX\n6ts2enjoBd18fY2OSAhDDGoXQkjTJv/pabAnLp8AvtqXg5eHG6O7SL+0xWRlwf33w+uvGx2JEIbz\ncHdjfPdwNqXlUVRebXQ4P+HSCaC61sSa1LPc2aklzX08jQ7HebRtC99/D3/4g9GRCGEX7k2IoKrW\nxNpD9lUawqUTwPYT+Vwsq5LBX0tZuhS2bNEfDx4MTaSekhAAXSKaERvqZ3eLwlw6AXy1Ty/9MLRD\nqNGhOL7qapg3D/76V6MjEcLuaJrGvQkR7DlVYFfbRbpsAiitrOG7I7nc060VXh4u+2OwHE9Pvdtn\n2TKjIxHCLk3oofc0rE49a3AkV7jsle/bQ7lUVJtk9s/t2rQJnn0WTCYIDQU/P6MjEsIuRQb50icq\niJUp2XZTGsJlE8DX+3OIDPKhZ5tAo0NxbJs26ZU9S0uNjkQIuzchIZyTF8o4fLbY6FAAF00AeSUV\n7EjPZ0L3CCn9cLteew1+/FE2bheiAcZ0aYWnu8aq/fYxGOySCeCb1HOYFLLt4606dQruuEOf769p\nUtxNiAYK9PNiaIdQVqeepdZkfDeQSyaAValn6dyqGe1byoXrluTm6hd/6fYRotEm9IjgfHEluzMv\nGh2K6yWAzPwyUs9cYmKC3P03Wt3AVb9+kJYG8fHGxiOEAxrZqSV+Xu6s2mf8bCCXSwCr959F02Bc\nd0kAjVJVBWPGwEcf6ceesnJaiFvh4+XOXfFhrD10zvAKoS6VAJRSrNqfQ9/oIFo1l8qfjVJVpU/1\nFELctgkJEZTYQYVQD0PPbmMHc/TKn7OHxBgdimOprYWmTWHdOnBzqXsGIaxiYGwwIU29WLU/h9Fd\nwgyLw6X+mlftP4unu8bdUvmz4d58E8aPh8uX5eIvhIV4uLtxT7dwvk/Lo7jCuAqhLvMXXWtSrEk9\ny7COLWjuK/3XDebnp0/z9PIyOhIhnMr4HuFU1Zj49lCuYTFYLQFomhagadoLmqbFmB/HaJo2T9O0\nkdY6543szrxIXkmlzP1vqJoa/fPs2Xp9H3d3Y+MRwskkRAYQGeRjaG0ga7YAgoB5wEmgEEgG9iql\nNlrxnNe1ev9Z/LzcGRHX0ojTO5ajR6FzZ9izRz+W1dJCWJymaUzoHsGO9HwulFQaEoO1u4DuBAKB\nWKVUoFJqhZXPd02VNbWsO5TLqPgwfLzkTvamPD2hZUu9uJsQwmrG9wjHpGDtQWM2irH6GIBS6pJS\nKsPa57mRrcfzKbpczXjp/rmxum6fdu1g61aIjjY2HiGcXIeW/sSF+RtWG8juBoE1TZutaVqSpmlJ\nFy5YZo7sqv05BPl5MahdiEVezynV1MDEifC73+nH0u0jhE1M6BFByulLhmwUY+0EEKNp2mTzx2xN\n0ybf7BuUUouVUolKqcRQC3RBlFXWsPHoecZ0DcPT3e7ynX1p0wZatzY6CiFcyrju+rR0IwaDrbkQ\nrACgfr+/pmlfaJqGLccCNhw5T0W16T+78YhrqKkBDw94912jIxHC5bQO9CWxbSCr95/lieHtbHpu\nq90Sm/v+F1/19CL0mUE2s2p/DuHNveklG79c2/LlkJioV/gUQhhiQo9wjp0vIS3XthvF3DQBmOfv\nN/Qj4CYvl4HeLXSzr7OIwrIqtp3IZ1yPcNzcpE/7mgID9W6fAJu8JUKIaxjTtRXubhqr99u2G+iG\nXUCapsXQuDv2vcB88/e+oJSaf9W/F5g/xwApjXjdW7L20DlqTIrxUvnz56qr9emed96pfwghDBPc\ntAkD24Ww5sBZnr+ro812KrxhAjBP35zS2BetSxyapq24agpokPmzTaaFrt5/lthQPzq3ku0KfyIv\nD4YNg1degSmNfnuFEFYwvns4z32Ryr4zl2y2V7lVxgDMF/0515j/PxJIUUpdssZ56ztXdJk9pwoY\nL/v+/pynJ0RF6bN+hBB24a74lnh5uNm0G8ia8yILzC0BQK8NBMwBHrPiOf/jm9RzKIUs/qqvtlb/\nCAyEtWuhb1+jIxJCmPl7ezIirgXfHDhHTa1t9t6w5iygFUBPc0G4eehjCVOUUlbv+wd9Tm231s2J\nDvGzxekcwzPP6F0+dSt+hRB2ZXz3cPJLK9mVUXDzL7YAq24IY1Ttn4wLpRzMKeL3YzsZcXr7FROj\nl3X2cKl9gIRwGMPjWtC0iQerU3MY1N76lQuccmns6lR93997ukn3DwCV5kqDTz0FCxYYG4sQ4rq8\nPd0ZFd+SdYdybbJfsFMmgGO5JfSNDiKsubfRoRhv82bo0AEOHjQ6EiFEA0zoEYFScOJ8qdXP5ZR9\nAX9/oBdlldLPDUBYGMTHQ4SUwhDCEQyMDSbp9yPx9rR+6XqnTAAAfk2c9r/WMOXl4OsLcXH6jB8h\nhEPwcHfDw0bbljhlF5DLKyyEXr2kv18IcUOSAJxR06YweLDM8xdC3JCL95M4mcpKqKoCf39YfHUh\nViGE+ClpATiTRx6BO+7Qk4AQQtyEtACcybRpcOKEvthLCCFuQhKAo1MK0tOhfXsYO9boaIQQDkS6\ngBzdwoXQrRscPmx0JEIIByMtAEc3bRqUlkLnzkZHIoRwMNICcFRr14LJBCEh8OKLIHseCCEaSRKA\nI9q6Ve/v//BDoyMRQjgw50wAJhN88om++YkzGjwYvvgCHnrI6EiEEA7MORPA+vUwcyZ89ZXRkViO\nUvD663DqlN7dM3kyuNuoYIgQwik5ZwK4+27YuBEmTTI6EsvJzob58+Hjj42ORAjhJJxzFpCmwYgR\n+uOcHLh8Gdq1MzamW2UygZsbREbCvn3Qtq3REQkhnIRztgDqmEx6a+CBB/QuFEdTXAwjR8KSJfpx\nVJTM9hFCWIxztgDquLnB3/+uF0dzxAunt7de019KOwghrMC5EwDAwIFXHq9frxdLs/cL6u7d+sIu\nf39Ys8Yxk5cQwu45dxdQfUeOwJgx8MYbRkdyY3l5MHw4/Pa3+rFc/IUQVuL8LYA6nTvD11/DXXcZ\nHcm11Q32tmgBy5bBkCFGRySEcHKu0wIAGD8emjTR98t9+ml960R7kJ4OCQmwa5d+PGECBAYaG5MQ\nwum5VgKok5QEixbBzp1GR6ILCdG3caysNDoSIYQLcc0EMGQIZGToYwIAWVm2j2H/fnj8cb3rJyAA\ntm+HoUNtH4cQwmW5ZgIAaNVK/3zsGHTqpNfVt6X9++HLLyEzUz+WwV4hhI25bgKoExMDc+fqtXVA\nX3xlMln+PFVV8Oqr8Pnn+vFDD8Hx4xAba/lzCSFEA0gC8PSEl1+GsDD9+OGH9dXDllJcfOU8a9bA\njh36saZB8+aWO48QQjSSJICr3XvvT6eKLlqkb7R+K+bN0/fqrazUL/hbtsDbb1smTiGEuE2SAK42\nYwY884z+OD0dfvUr+O47/fj8eX2BVlqafpyXB8uX658BNmzQi87l5OjHiYnwxBNQXa0f+/ra7v8h\nhBA3IQngRtq1g7NnYfp0/TgzExYsgDNn9OO0NLjvPjh0SD8ODYUePa6sLxgxAl56SZ/iKYQQdkZT\ndlwlMzExUSUlJRkdxk/V7TLm7g5lZfoGLVFR4OdnZFRCCPEfmqYlK6USb/Z1rlMKwlLq78Ll5wfx\n8cbFIoQQt0G6gIQQwkVJAhBCCBclCUAIIVyUJAAhhHBRkgCEEMJFSQIQQggXJQlACCFclF0vBNM0\n7QJgQLF+uxMC5BsdhPgJeU/sk7wvurZKqdCbfZFdJwCh0zQtqSGr+oTtyHtin+R9aRzpAhJCCBcl\nCUAIIVyUJADHsNjoAMTPyHtin+R9aQQZAxBCCBclLQAhhHBRkgCEEMJFyX4AQjSCpmk9gUQgA4gB\nMpRSG42NyrVpmjYZ/b2INX9epJRaYWxUjkESgBANpGlaDDBPKXVnvee+0DQtQymVYWBoLst88c+o\nu+BrmhYAJGuaFqSUkgHhm5AuIAejadpkTdNe0DRtkaZpG8x/AMI25gCLrnpuETDPgFiELkYplVJ3\noJS6hP5+XP0+iWuQWUAOpN7dTor5OABIRr8rlbsdK9M07SRwZ/27ffN7UKiU0oyLzDWZf/bfAyPM\nF/6652OAk0CstMxuTFoAjkXudgxivtjEAAX1n6+78JgvOsKGzD/7GPOHuAWSAByE+QJ0n/lzfRvN\n/y5/BNYVBFcu+NcgP38DKKUC698UmY0ELsnd/81JAnAQcrdjuKsTr7Bfc4DXjQ7CEcgsIAeilAq8\nxtNytyOEmaZps4ECpdR8o2NxBNICcHxyt2ND1+iCE3bC3A06p/40XXFjkgAcmNzt2FRdCyuo/pP1\nEoK0wIw3DxhhdBCORLqADNDIAduCaw081rvb6WW5yMT1KKUuaZqWwc/HAoKQLjjDaZq2CJh7g0F6\ncQ2SAGysbjVpI75lL3CtO3y527G9jehlIOrPOulpfl4YxNwSnnfV+oyR6GtmJDHfgCwEc0Dmu515\n8sttW+buni+uKgWxAb0lJu+FAcyLI4OApHpPBwFTlFJzjInKcUgCcDDmu52NcrdjDHMxuPvQW2Yx\nQIoUgzNG3Srs6/xzhlIq1pbxOCJJAA5E7naEEJYkCcBByN2OEMLSJAEIIYSLknUAQgjhoiQBCCGE\ni5IEIIQQLkoSgBBCuChJAEII4aIkAQghhIuSBCCEEC5KEoAQQrio/w9Lwtf+FxIKNgAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119dc0cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0\n",
       "[torch.FloatTensor of size 100]"
      ]
     },
     "execution_count": 230,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def fun(v):\n",
    "    return torch.sin(v)*v**2 + 2*v\n",
    "\n",
    "v  = Variable(torch.linspace(-3,3,100), requires_grad=True)\n",
    "LL = torch.sum(fun(v))\n",
    "LL.backward()\n",
    "\n",
    "plt.plot(v.data.numpy(), v.grad.data.numpy())\n",
    "plt.plot(v.data.numpy(), fun(v).data.numpy(), 'r:')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "v.grad.zero_()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 0.9093\n",
       "[torch.DoubleTensor of size 1]"
      ]
     },
     "execution_count": 224,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.linspace(-3,3,100)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.1"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}