5_intro_gans-full.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GANs for the MNIST dataset\n",
    "\n",
    "In this notebook we will implement a generative adversarial network (GAN) for generating digits resembling those from the MNIST data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torch.optim as optim\n",
    "import numpy as np\n",
    "from torchvision import datasets, transforms\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GANs on a nutshell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GANs are relevant tools for generating data from complex probability distributions (especially, probability distributions over images). They view the sampling process as a game between two agents: a **generator** (*G*) and a **discriminator** (*D*).\n",
    "\n",
    "*G* generates samples of data intended to come from the same distribution as the real data $p_{data}$. *D* is trained to distinguish real data from fake data (generated by *G*). This is a simple classication task. In turn, *G* is trained to fool *D* in its task. Each player is seen as a diffentiable function with some parameters $\\theta^{(G)}$ (for the generator) and $\\theta^{(D)}$ (for the discriminator). The discriminator is simply a function that maps a given image into a scalar value, representing the probability of the image being real. The generator, maps a prior noise (usually belonging to a low dimension latent space) into an image of the same size as the real images.\n",
    "\n",
    " <img src=\"img/gan.jpeg\">\n",
    "         \n",
    "#### Cost functions\n",
    "\n",
    "One popular choice for the discriminator's cost function is the Binary Cross Entropy. In that case, assuming that the data is balanced (same number of fake and real images), the cost function is\n",
    "\n",
    "$$\n",
    "C(\\theta^{(D)}, \\theta^{(G)}) = -\\frac{1}{2} \\mathbb{E}_{x \\sim p_{data}}[\\log(D(x))] - \\frac{1}{2} \\mathbb{E}_{z}[\\log(1- D(G(z)))] \n",
    "$$\n",
    "\n",
    "The discriminator will choose $\\theta^{(D)}$ minimizing these function.\n",
    "\n",
    "Under the minimax setting, the loss for the generator is just $-C(\\theta^{(D)}, \\theta^{(G)})$. Thus, the optimal value (equilibrium) of the generator's parameters is\n",
    "\n",
    "$$\n",
    "\\theta^{*(G)} = \\arg\\min_{\\theta^{(G)}} \\left \\lbrace \\max_{\\theta^{(D)}} [-C(\\theta^{(D)}, \\theta^{(G)})] \\right \\rbrace\n",
    "$$\n",
    "\n",
    "As this objetive is intractable to optimize exactly (as requieres solving a bilevel optimization problem), the usual approach is to perform silmultaneous stochastic gradient descent: at each iteration, we first uptade the parameters of one of the players, and then the other's. \n",
    "\n",
    "Now we will try to implement this approach to generate digits from MNIST."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GANs on MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# De-normalization\n",
    "def denorm(x):\n",
    "    out = (x + 1) / 2\n",
    "    return out.clamp(0, 1)\n",
    "\n",
    "def plot_result(generator, noise, num_epoch=1, save=False, save_dir='MNIST_GAN_results/', show=True, fig_size=(5, 5)):\n",
    "    generator.eval()\n",
    "\n",
    "    noise = noise.to(device)\n",
    "    gen_image = generator(noise)\n",
    "    gen_image = denorm(gen_image)\n",
    "\n",
    "\n",
    "    n_rows = np.sqrt(noise.size()[0]).astype(np.int32)\n",
    "    n_cols = np.sqrt(noise.size()[0]).astype(np.int32)\n",
    "    fig, axes = plt.subplots(n_rows, n_cols, figsize=fig_size)\n",
    "    for ax, img in zip(axes.flatten(), gen_image):\n",
    "        ax.axis('off')\n",
    "        ax.set_adjustable('box')\n",
    "        ax.imshow(img.cpu().data.view(image_size, image_size).numpy(), cmap='gray', aspect='equal')\n",
    "    plt.subplots_adjust(wspace=0, hspace=0)\n",
    "    title = 'Epoch {0}'.format(num_epoch+1)\n",
    "    fig.text(0.5, 0.04, title, ha='center')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "image_size = 28\n",
    "batch_size = 128\n",
    "Gn_hidden = [256, 512, 1024]   # Sizes of the hidden layers\n",
    "Dn_hidden = Gn_hidden\n",
    "dim_z = 100\n",
    "lr = 0.0001\n",
    "n_out = 28*28\n",
    "num_epochs = 100\n",
    "device=\"cuda\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We are going to load the dataset. Since the MNIST dataset is very popular, PyTorch already has it included into datasets.MNIST, and will automatically download it the first time. We can also specify some transformations, like converting the matrix to a tensor and normalize the pixels\n",
    "\n",
    "Note that we create two data loaders: one for the training set and another one for the testing set. Since this is just a quick experiment, we are not going to validate the initial hyperparameters, so we don't create another data loader for a validation set.\n",
    "\n",
    "The main benefit of the DataLoader class is that it creates an object which can be used in a for loop, as we did in the previous train function, to iterate over the minibatches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_loader = torch.utils.data.DataLoader(\n",
    "    datasets.MNIST('../data', train=True, download=True,\n",
    "                   transform=transforms.Compose([\n",
    "                       transforms.ToTensor(),\n",
    "                       transforms.Normalize((0.5,), (0.5,))\n",
    "                   ])),\n",
    "    batch_size=batch_size, shuffle=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize some examples together with the true labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([512, 1, 28, 28])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "examples = enumerate(train_loader)\n",
    "batch_idx, (example, label) = next(examples)\n",
    "example.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADhCAYAAAD/Ec//AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAHUhJREFUeJzt3XmUVMXZx/HvI6vBDYSgrAqvCyhiiAbcCFFcMKJGgzHH/TXGJUaNuKLEqFERd4JGjcYo7uKOIsQFEY1yVNQXEzGI4K6sKoqKcN8/eqqmeqZnpqeZ7lvd8/ucM4eH6r7d1c/0VFfXrVtlSZIgIiLpWyvtCoiISIYaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUhE2yCb2TQz+02pj20OlNviUn6Lp9JzW/QG2czmm9nQYj9PoczsSDNbZWbLg58hadcrH2WQ263NbIqZLTKzspvwXgb5bWNmV5nZR2a21MyuM7NWadcrH8ptbtH2kEvsX0mSrBP8TEu7QhViJXAvcHTaFalQZwHbAVsDmwMDgHNTrVHlSCW3qTXIZtbezCaZ2cKqT6BJZtatxt16m9lMM/vCzB42sw7B8YPM7AUzW2Zmr5dLr7YUYsltkiRzkiS5GXhzDV5OdGLJLzAcGJckyZIkSRYC44D/LfCxotDcc5tmD3kt4BagJ9ADWAGMr3Gfw8kkYWPgezJJwcy6Ao8BfwY6AKcB95tZp5pPYmY9qn45Peqpy4+qvla/bWajzazlmr201MWU20oUU36tRtzNzNYv5EVFonnnNkmSov4A84GhedxvW2Bp8P9pwJjg/32B74AWwJnAhBrHTwGOCI79TZ716wVsSuaN0A/4N3B2sfPSHHIbHP8/mbda+jmrpPySaXieBzoBGwEvAQmwcdq5U24L+0lzyOIHZnaDmS0wsy+A6cAGZtYiuNv7QbwAaAV0JPPpOaLqE26ZmS0DdibzidkoSZLMS5Lk3SRJVidJ8n/ABcAvC31dMYglt5UqovxeBMwCXgNeAB4iM27/aQGPFYXmnts0hyxGAlsAA5MkWQ8YXFUefk3oHsQ9yCRkEZlfyIQkSTYIftolSTKmCeqV1KhDOYo1t5UiivwmSbIiSZITkyTpmiRJL2Ax8EqSJKsLeVGRaNa5LVWD3MrM2gY/LYF1yYwPLasalD8vx3GHmllfM/sBmZ7rxCRJVgG3A8PNbE8za1H1mENyDP43yMyGmVnnqnhLYDTwcIGvMw0x59bMrC3Quur/bc2sTaEvNCUx57ermXWpyvMgMu/dXHWJlXJbU4nGipIaP38GupAZ01kOvA0cW3Vby2C85xJgJvAF8CjQMXjcgcCzwBJgIZnB/B41x4rIfIIud7flqN/lZL6GfAXMI/MLblXsvDST3G6So37z085bBeV3cFUdvwbmAIeknTPlds1+rOrJRUQkZbowREQkEmqQRUQioQZZRCQSapBFRCKhBllEJBKNWrPBynAJxVJLkqSgi0qU27wsSpKk1roEDVFu81JQbkH5zVNe+VUPWcrJgrQrUMGU2+LKK79qkEVEIqEGWUQkEmqQRUQioQZZRCQSapBFRCKhBllEJBJqkEVEIqEGWUQkEtHvrtyyZaaKa6+9dt7HjBo1CoBwred99tnHx/369cvrcV588UUf77777gAsX74873qUs5EjR/r48ssv9/HGG2e2J/vkk09KXieRSqcesohIJNQgi4hEIvohi8suuwyAk08+Oe9jzDLr+4RDFp9//rmPP/74YwDef796N/F58+YB8MILL/iy2267zcfNZajC6dOnj49Xr67eaPeOO+4AYLfddit5naR5W2ut6v7jqaeeClS3D6F3333Xx7vuuquP58+fX7zKNRH1kEVEIqEGWUQkElEOWXTt2tXH66+/PgDPP/+8L+vcuXOtY8KhhmeffRbIHrKYMWOGj+fOndt0lW1mttlmGwB69erly9xwj1Tr2bMnAIMHD/ZlLu7bt68v23HHHWsd+8EHH/j44osv9vFTTz0FwNtvv920lS0TYd7Gjh0LZP+NO5tssomPn3zySR9fcMEFQPZQZGzUQxYRiUSUPeSbbrrJx/379weyexLlMDhfrtyJk7CXEerQoQMA3bp182XNuYccnmj6xS9+4eNx48YB1fO2obo39/XXX/uyyZMn+9i91zfYYANfdvXVV/v4rbfeAuDYY4/1ZeFc+Uq0zjrr+Di8liCX77//vlZZ+E1u9OjRAEydOtWXxTafXj1kEZFIqEEWEYlElEMWoTZt2gDVl1BLcbmv1eG8balt6623BuCcc87xZb/61a98/M477wBw4403+rKHHnoIgCeeeKLRzwNw6aWXAjBhwgRfttlmmzWm2mXnuuuu8/Ghhx5a6/aVK1f6+MwzzwSyhzSvvPJKH/fu3RuADTfc0JdpyEJERHKKptu5+eab+/gnP/mJj91JC01VK41WrVoBsP/++6dck/gMGzbMx/feey+QvejV7bff7mO3ONPChQvX6Dlnz57t44MOOgjI7mG7aaGtW7f2ZcuWLfNx2IMsJx07dgRghx12yHn7d999B8DZZ5/ty8IToE54lemFF14IwPnnn+/Lpk+f7mPXG891crBU1EMWEYmEGmQRkUikPmSx7rrrAtVX3gC0b9/ex3//+99LXiep2+LFiwF47733Uq5J6fz85z8H4O677/Zl7dq1A+BPf/qTL3NXghWL+1txJwyheigvPFE1aNAgH8+cObOodSoWN//YnYir6c033wTgqquuqvdxHnnkER8feOCBABx22GG+7IADDqh1jJtDngb1kEVEIqEGWUQkEqkPWbiFgvbdd9+ct3/66aelrI40wP0+Kv3y9XBGxX333QfAqlWrfNlRRx0FZM+sKIYhQ4b4+JprrgFyb0HmLqsGeO2114pap1IIL0PPZeLEiXk9zp577unj4cOH13vffLd2Kyb1kEVEIpF6D3n77bevVRbu5BEuRSjFF15t1tyEJ8auv/56H7dt2xaonscKcOuttxatHiNGjPBxuPymq8ecOXN82RZbbAFkL8np5uiWGzcHHrKXLc3lP//5T523rbfeej4OdxQJF22KlXrIIiKRUIMsIhKJ1IcsXn/99Vpl3bt39/G///1voHreIcADDzxQ65grrrjCx+X6lS0G4RzwXKZNm1aaiqQgHJII34OXX345kH3JbVPp0qWLj48//nigepEcgMcff9zHbvgiXM/XOeKII5q8bqXmLpcG2G+//Rp9vBuquPnmm31ZuDhTOVAPWUQkEqn3kN3yd+EyhgcffLCPt9xySyD7k2677bar9TjhyY+99toLgGeeecaXqdfcNMqtx9EY4QI9oVmzZgHZC9UUwp00POOMM3zZkUce6eNOnToB8OGHH/qyUaNG+dgtwBWetHJb3q9YsWKN6haDJUuW+NgtoOT+lms67rjjAFi6dKkvO/HEE4HqK/LKkXrIIiKRUIMsIhKJ1Ics3NeUSy65xJeFsRueCNd4dVdRnXTSSb4s3MzQfd0JT4iEczvDTSYl2zHHHFPv7eEc2EqTa0t5qP7a3JjXPnDgQCD76q+dd94ZADPzZeedd56P3dVn4dfwcLNZt5BOOPzm5o1XwpDct99+6+PnnnsOqHvIYo899sj6t1KohywiEgk1yCIikUh9yMJxm5lC9oamL7/8cq37/uUvf8n6F7I3Q3TzOcMFYu666y4fu80Sv/zyyzWtdsVZa636P6PD+eCVJpyHHM6Jde+XcB3dfIWX+N5///0AjBkzxpc1NHwWPmfPnj0BmDRpki/L9fdRCdzmpLvssosvq2v4oia37RtkXzrthkLD7eJiox6yiEgkUu8hb7vttkD2J1l41V2+W6afcMIJPnYnBNxyhZC99N5tt90GZF/d9MUXXzSm2hXFLVoD2Qu85FLJO7iEu6CEyz+6hW7cnHio7jXfeeedvixXbzd8/3722Wd51WPo0KE+dpulQvUJ8GLvTBIDd4Lv3HPP9WXhe9PtKPLRRx/5Mvf7Gz16tC9bvny5j8PyXMIlTNOiHrKISCTUIIuIRMLqmnuZ885m+d+5Hv379/exO0HxzTff+LIdd9zRxwsXLiz4eU4//XQfX3rppbVu33vvvX2c79BIQ5IksYbvVVtT5bYQP/3pT3389NNP13tfN7yxcuXKotapDq8kSVL7uvkGpJnbxnCbmM6YMcOXhfOY3c4lRVqzuqDcQunyGw6tucvcGzPU+OqrrwLVw6Q1uTnfRdrAN6/8qocsIhIJNcgiIpFIZZbF0Ucf7eOuXbsC2RszrskwRSicZRGub+tWherdu3eTPE+5GzBgQNpVEOCee+4BsocpZs+e7eNTTjml5HWKSTisGcb1CZdX6NOnT63b582b5+OvvvpqDWrXNNRDFhGJRCo95FwbFLrNGiG7Bx2u/t9Y4YIr48aN87HrIYcnR6699tqCn6fc5dpoVkqjR48ePs71e/jjH//o448//rgkdSp37ls3wMknn+zj8GpgZ/z48T5evHhxcSuWB/WQRUQioQZZRCQSqQxZ/O1vf/OxW8Ql3EAyXDTorLPOArKHOR588MF6H/+WW24BYNCgQb5syJAhhVdYpAmFCziFCxq5LZ7C9/eUKVNKV7EKcfjhh/t40003rfe+4UJNMVAPWUQkEqn0kL///nsfX3/99UD2lJRwk1M3NS2corbPPvvU+/humb127dr5sjB2wmlxUrfp06f7eNWqVSnWpDK0b9/ex2Fvzk3lCjf8rYTNS0st32U6Y6QesohIJNQgi4hEIvX1kN1VeeHOCOFGm648vNLuuOOOA7K/+oU6d+4MZG9aGW4c6U4gPvzww2tU9+Zi7ty5Pl69enWKNSlvHTp0AODRRx/1ZeF79IEHHgDiWJe3nIU7soQ7jpQD9ZBFRCKhBllEJBKprIdcycpxPeQyUtbrIbutmaZOnerLwgVtfvSjHwHZQ0QlFP16yPnafffdfZxrHne44FC4ZEJj2sICaD1kEZFykvpJPZHmomXLzJ9b2BNzm6VCaj3jivPPf/7Tx+FVkeWgvGorIlLB1CCLiERCQxYiJeI20m3RokXKNZFYqYcsIhIJNcgiIpFQgywiEgk1yCIikWjsSb1FwIJiVKRC9FyDY5XbhhWaX+W2YXrvFlde+W3UpdMiIlI8GrIQEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSETbIJvZNDP7TamPbQ6U2+JSfoun0nNb9AbZzOab2dBiP8+aMLNeZjbJzL40s0VmNjbtOuVDuS2u2PNrZlub2ZSqvCZp16cxyiC3bczsKjP7yMyWmtl1Ztaq2M8bbQ+5VMysNfBP4GlgI6AbcHuqlaoQym3RrQTuBY5OuyIV6CxgO2BrYHNgAHBusZ80tQbZzNpX9ZwWVn0CTTKzbjXu1tvMZprZF2b2sJl1CI4fZGYvmNkyM3vdzIYUWJUjgY+SJLkySZKvkiT5JkmSNwp8rCgot8UVS36TJJmTJMnNwJtr8HKiEktugeHAuCRJliRJshAYB/xvgY+VtzR7yGsBtwA9gR7ACmB8jfscTiYJGwPfk0kKZtYVeAz4M9ABOA2438w61XwSM+tR9cvpUUc9BgHzzWxy1Ve/aWbWb41fXbqU2+KKJb+VKKbcWo24m5mtX8iLyluSJEX9AeYDQ/O437bA0uD/04Axwf/7At8BLYAzgQk1jp8CHBEc+5s86zeVzFe/YUBr4HRgHtC62LlRbuP+iT2/wfH/k/lTTj9nlZJbMo3680AnMsNtLwEJsHEx85LmkMUPzOwGM1tgZl8A04ENzKxFcLf3g3gB0AroSObTc0TVJ9wyM1sG7EzmE7OxVgAzkiSZnCTJd8DlwIZAnwIeKwrKbXFFlN+KE1FuLwJmAa8BLwAPkelcfFrAY+UtzSGLkcAWwMAkSdYDBleVh18TugdxDzIJWUTmFzIhSZINgp92SZKMKaAeb5D55Kskym1xxZLfShRFbpMkWZEkyYlJknRNkqQXsBh4JUmS1YW8qHyVqkFuZWZtg5+WwLpkelDLqgblz8tx3KFm1tfMfgBcAExMkmQVmTP1w81sTzNrUfWYQ3IM/ufjdmCQmQ2t+hQ+hcwv9z+FvNAUKLfFFW1+LaMtmeEgqh6rTaEvNAUx57armXWpyvEgYHQddWlaJRorSmr8/BnoQmZMZznwNnBs1W0tg/GeS4CZwBfAo0DH4HEHAs8CS4CFZAbze9QcKyLzCbrc3VZHHQ8A5lY9zzRgq2LnRbmN/yf2/AKb5Kjf/LTzViG5HVxVx6+BOcAhpciLVT25iIikrNlfGCIiEgs1yCIikVCDLCISCTXIIiKRUIMsIhKJlo25s5XZEn9pSJLEGr5XbcptXhYlSVJrXYKGKLd5KSi3oPzmKa/8qocs5WRB2hWoYMptceWVXzXIIiKRUIMsIhKJRo0hS/M0YsQIH99999313rdFixb13i4idVMPWUQkEmqQRUQiUTFDFm3btvXxZZddBsDvfvc7X2ZWPRtt1KhRAFxyySUlql1522abbXycazGqCy+8sJTVEalY6iGLiESirHvIa6+9to+ffPJJHw8YMADI7gG//fbbPh4/PrNn4qOPPurLZs+eXbR6lqNf//rXPh45cmSt22fOnOnjG2+8sSR1qnSPPPIIAMOHD/dlixYt8nH//v0B+Oijj0pbMSkZ9ZBFRCKhBllEJBJlPWRx5513+rhv374+3n///QGYMmVKzuO++eYbAHr16uXLmvOQRe/evX188cUXA7D55pv7statW/v422+/BbJP5H388cfFrmKz4E6Yrlq1ype1b9/ex/fddx8AO+20U2krViE6duwIwNSpU33Ziy++6OMTTjih5HWqST1kEZFIqEEWEYlEWQ5ZXHnllQDsu+++vuzggw/2cV1DFc7EiROB3HNqmyM3TAFw4IEHAtnztsM8nX/++QBMnjy5RLWrbIMGDfJxjx49UqxJZdp999197NqNrbbaype1bBlXE6gesohIJOL6eKhH+El38sknA/Dggw/6skmTJuX9WOFJk+bsrLPOAuCXv/xlrdvWWqv6s/qee+7x8dixY4tfsWYk/GbXr1+/FGtSnsIerjspt9tuu/myPfbYw8fuZP6sWbN8Wc+ePX3cpUsXIN153uohi4hEQg2yiEgkoh+ycJdHX3PNNb7MXbZ75JFH+rIVK1aUtF7l6pxzzqkV5zq5+dZbb/nYDW1I0xg2bJiP99xzzxRrUr5++MMfAtnDE+eddx4AHTp08GWjR4/28fXXX591LMArr7zi41133RWA22+/vQg1zo96yCIikYi+h+yW0Nx000192UEHHQTA8uXLU6lTuQkXCgp7yOEVeM78+fMB2GuvvXzZggX178+4ySab+HijjTYCYMmSJb4sXNhJsq+C3GyzzVKsSXnp3r27j2+44QYg+326bNkyAI466ihfdtttt/l49erVQPaCTR988IGP3Xs3Teohi4hEQg2yiEgkohyyCHf/OPbYYwGYO3euLwvXMZaGnXHGGT7ONUyxePFiH7s1pBsapggXJArngLuv4++//74vc/OcX3755cZUW8g+4docT1zvvPPOPv7rX//qY3e13fTp033ZYYcdBmS/93IJ25dWrVr5OFygLC3qIYuIREINsohIJKIcsnAL3ED1V+MJEyakVZ2ytd9++wENb1IabgbrFl6qi5tRccwxx/iyXDMFunXr5uOrrroKgF122SWPWleu9dZbD2jc2fyFCxf6eOjQoU1ep1i5S5rDedpbbrmlj2+66SYATj/9dF/2+eef5/XYI0aM8HG4oJMrDxfbCodKS0E9ZBGRSETZQw57dI4WtcnPPvvs4+NTTjkFyF4oyM3FhOpFVP773//W+5jhvNknnngCyF6UJeSeK3yeHXfcEchexKihnnilcL1iqO7Nhb26hoSb91a68D3l3scnnXSSL3O9YoDjjz++0Y/vrtCrK/9u55tS94pD6iGLiERCDbKISCSiHLLo37+/j5cuXQo0PC82XBfVzSfcYostfNkbb7zhYzfvNryEslK4uZhQfRItHD547733fOxOYrz++uu1Hiccpnjsscd87E6C1LXbiltrNtewU3MUnjQ6++yz8zomXNwmPHla6cL5w6eeeiqQfVIzPPns/t7DecQbb7wxkL0Y1r/+9S8ff/LJJwB89913OZ//2WefLbjuTUU9ZBGRSETZQ3YngQDeffddIPdCQm6Ff4DrrrvOx26vPbfYCGT36Nyn6qhRo3zZ+PHj17TaJbfhhhv62O2o8rOf/azeY8LdP+q7cs6dvIPce725JVAB/vCHP/j4rrvuqvf5m4s2bdoAcOaZZzb62Msuu8zHK1eubLI6xc7tBATZi4k5YS7dUplDhgzxZe7b2dNPP+3LwqU43Y4hnTt39mVhG3HzzTcXWvUmox6yiEgk1CCLiEQiyiGLOXPm+DjXYjg//vGPgeyvx25AH6o3OwxvD4cs3MaSV1xxhS97+OGHgYYXJomJ+9oG9e9y8OWXX/r46quvzuuxwzmhuU7guS3VIfv35Y4Lj3HrIb/44ot5PXclcMNi4VrUUr/w6s5cLrroIh/PmDEDyB46c+/t++67L+fxbh5yu3btfFm4HnIM70/1kEVEIqEGWUQkElEOWUybNs3HJ554IpA9r3Xy5MlA9jZB+++/v4+feuqpeh//xhtvBOC3v/2tL3MzM6699toCax0vt8gQVM/FXFPhfOeGZhL84x//ALK/HlY69x5tjN///vcAvPPOO01dnbIQDn1NmTIFyJ5dFc4fdkMV4eyoVatW1fv4ffr0AbLnLsdGPWQRkUhE2UMOr6pz8zndJyZAixYtgOzFambPnt3o53n++ed97ObvllMP2cxyxo678qiQK5DqWpDICRcxynWcO7EK2Ts9VJr111/fx48//riPBw4cWOcx4c4f4W4u7ptbQz29SuVef824qbjNkcOTeg888ECTP8+aUA9ZRCQSapBFRCIR5ZBFODzhTgSFcxTdfMRChimg+uu9m5cIcOuttxb0WGnq16+fj3PNFZ4/f369x7tLr8P5zE44TFHXQkK5uN9TU508jF04HFPfMEUo/DpeycM5MQgXGHPXH4ReeumlUlanQeohi4hEQg2yiEgkohyy+Oyzz3z86quvAtlDFu4saXhZda41TsPb3RrJUL39SzhT4JxzzlnTapfc4MGD67398MMPB7I3hwyHH9ZZZx0AttpqqzWqh7vsHKrneDbXmQISl3Bt9fbt2wPZs7hi2yJLPWQRkUhE2UMO5VpwxG3eecghh/iyO+64w8fuk3Dvvff2ZR07dvTxa6+9BsC2227ry+bNm9dENS6d8CqlnXbaqc77hSebGnOCLl/hlYBu7eRwy/p8t2evdG4z2dh6ZZUsvFbBcTsGQXzrTauHLCISCTXIIiKRiH7Iwm3KOWDAAF/mNj7t1KmTL3PDGFD9tTy8lDW8JDrcnqichZuTjhs3DoCTTjqpSR575MiROcvdAjjhesmh0047DdAwRS5vvfUWUNjCQ5K/YcOG+fiAAw7w8bfffgtkb5EVG/WQRUQiEX0PecyYMQBsv/32vqxr164AfPjhh75s+vTpPh47diyQe3v7SuJ24oDqbdPdv8WS744jzZ1bQMgt6wrZW9pL8fTq1cvH4SJZzz33HBD3N2T1kEVEIqEGWUQkEtEPWbidAbp3755yTUTy5xZneuaZZ1KuSfMT7mYTmjhxYolr0njqIYuIREINsohIJKwxl9KaWdNfd1thkiSpvZdSHpTbvLySJMl2jT1Iuc1LQbmF+PIbruU9a9YsH++www5A7oXISiCv/KqHLCISiehP6omINEY497jclG/NRUQqjBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJRGPnIS8CFhSjIhUi9zYa+VFuG1ZofpXbhum9W1x55bdRl06LiEjxaMhCRCQSapBFRCKhBllEJBJqkEVEIqEGWUQkEmqQRUQioQZZRCQSapBFRCKhBllEJBL/D40KFTN0qR10AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADhCAYAAAD/Ec//AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAHUhJREFUeJzt3XmUVMXZx/HvI6vBDYSgrAqvCyhiiAbcCFFcMKJGgzHH/TXGJUaNuKLEqFERd4JGjcYo7uKOIsQFEY1yVNQXEzGI4K6sKoqKcN8/eqqmeqZnpqeZ7lvd8/ucM4eH6r7d1c/0VFfXrVtlSZIgIiLpWyvtCoiISIYaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUhE2yCb2TQz+02pj20OlNviUn6Lp9JzW/QG2czmm9nQYj9PoczsSDNbZWbLg58hadcrH2WQ263NbIqZLTKzspvwXgb5bWNmV5nZR2a21MyuM7NWadcrH8ptbtH2kEvsX0mSrBP8TEu7QhViJXAvcHTaFalQZwHbAVsDmwMDgHNTrVHlSCW3qTXIZtbezCaZ2cKqT6BJZtatxt16m9lMM/vCzB42sw7B8YPM7AUzW2Zmr5dLr7YUYsltkiRzkiS5GXhzDV5OdGLJLzAcGJckyZIkSRYC44D/LfCxotDcc5tmD3kt4BagJ9ADWAGMr3Gfw8kkYWPgezJJwcy6Ao8BfwY6AKcB95tZp5pPYmY9qn45Peqpy4+qvla/bWajzazlmr201MWU20oUU36tRtzNzNYv5EVFonnnNkmSov4A84GhedxvW2Bp8P9pwJjg/32B74AWwJnAhBrHTwGOCI79TZ716wVsSuaN0A/4N3B2sfPSHHIbHP8/mbda+jmrpPySaXieBzoBGwEvAQmwcdq5U24L+0lzyOIHZnaDmS0wsy+A6cAGZtYiuNv7QbwAaAV0JPPpOaLqE26ZmS0DdibzidkoSZLMS5Lk3SRJVidJ8n/ABcAvC31dMYglt5UqovxeBMwCXgNeAB4iM27/aQGPFYXmnts0hyxGAlsAA5MkWQ8YXFUefk3oHsQ9yCRkEZlfyIQkSTYIftolSTKmCeqV1KhDOYo1t5UiivwmSbIiSZITkyTpmiRJL2Ax8EqSJKsLeVGRaNa5LVWD3MrM2gY/LYF1yYwPLasalD8vx3GHmllfM/sBmZ7rxCRJVgG3A8PNbE8za1H1mENyDP43yMyGmVnnqnhLYDTwcIGvMw0x59bMrC3Quur/bc2sTaEvNCUx57ermXWpyvMgMu/dXHWJlXJbU4nGipIaP38GupAZ01kOvA0cW3Vby2C85xJgJvAF8CjQMXjcgcCzwBJgIZnB/B41x4rIfIIud7flqN/lZL6GfAXMI/MLblXsvDST3G6So37z085bBeV3cFUdvwbmAIeknTPlds1+rOrJRUQkZbowREQkEmqQRUQioQZZRCQSapBFRCKhBllEJBKNWrPBynAJxVJLkqSgi0qU27wsSpKk1roEDVFu81JQbkH5zVNe+VUPWcrJgrQrUMGU2+LKK79qkEVEIqEGWUQkEmqQRUQioQZZRCQSapBFRCKhBllEJBJqkEVEIqEGWUQkEtHvrtyyZaaKa6+9dt7HjBo1CoBwred99tnHx/369cvrcV588UUf77777gAsX74873qUs5EjR/r48ssv9/HGG2e2J/vkk09KXieRSqcesohIJNQgi4hEIvohi8suuwyAk08+Oe9jzDLr+4RDFp9//rmPP/74YwDef796N/F58+YB8MILL/iy2267zcfNZajC6dOnj49Xr67eaPeOO+4AYLfddit5naR5W2ut6v7jqaeeClS3D6F3333Xx7vuuquP58+fX7zKNRH1kEVEIqEGWUQkElEOWXTt2tXH66+/PgDPP/+8L+vcuXOtY8KhhmeffRbIHrKYMWOGj+fOndt0lW1mttlmGwB69erly9xwj1Tr2bMnAIMHD/ZlLu7bt68v23HHHWsd+8EHH/j44osv9vFTTz0FwNtvv920lS0TYd7Gjh0LZP+NO5tssomPn3zySR9fcMEFQPZQZGzUQxYRiUSUPeSbbrrJx/379weyexLlMDhfrtyJk7CXEerQoQMA3bp182XNuYccnmj6xS9+4eNx48YB1fO2obo39/XXX/uyyZMn+9i91zfYYANfdvXVV/v4rbfeAuDYY4/1ZeFc+Uq0zjrr+Di8liCX77//vlZZ+E1u9OjRAEydOtWXxTafXj1kEZFIqEEWEYlElEMWoTZt2gDVl1BLcbmv1eG8balt6623BuCcc87xZb/61a98/M477wBw4403+rKHHnoIgCeeeKLRzwNw6aWXAjBhwgRfttlmmzWm2mXnuuuu8/Ghhx5a6/aVK1f6+MwzzwSyhzSvvPJKH/fu3RuADTfc0JdpyEJERHKKptu5+eab+/gnP/mJj91JC01VK41WrVoBsP/++6dck/gMGzbMx/feey+QvejV7bff7mO3ONPChQvX6Dlnz57t44MOOgjI7mG7aaGtW7f2ZcuWLfNx2IMsJx07dgRghx12yHn7d999B8DZZ5/ty8IToE54lemFF14IwPnnn+/Lpk+f7mPXG891crBU1EMWEYmEGmQRkUikPmSx7rrrAtVX3gC0b9/ex3//+99LXiep2+LFiwF47733Uq5J6fz85z8H4O677/Zl7dq1A+BPf/qTL3NXghWL+1txJwyheigvPFE1aNAgH8+cObOodSoWN//YnYir6c033wTgqquuqvdxHnnkER8feOCBABx22GG+7IADDqh1jJtDngb1kEVEIqEGWUQkEqkPWbiFgvbdd9+ct3/66aelrI40wP0+Kv3y9XBGxX333QfAqlWrfNlRRx0FZM+sKIYhQ4b4+JprrgFyb0HmLqsGeO2114pap1IIL0PPZeLEiXk9zp577unj4cOH13vffLd2Kyb1kEVEIpF6D3n77bevVRbu5BEuRSjFF15t1tyEJ8auv/56H7dt2xaonscKcOuttxatHiNGjPBxuPymq8ecOXN82RZbbAFkL8np5uiWGzcHHrKXLc3lP//5T523rbfeej4OdxQJF22KlXrIIiKRUIMsIhKJ1IcsXn/99Vpl3bt39/G///1voHreIcADDzxQ65grrrjCx+X6lS0G4RzwXKZNm1aaiqQgHJII34OXX345kH3JbVPp0qWLj48//nigepEcgMcff9zHbvgiXM/XOeKII5q8bqXmLpcG2G+//Rp9vBuquPnmm31ZuDhTOVAPWUQkEqn3kN3yd+EyhgcffLCPt9xySyD7k2677bar9TjhyY+99toLgGeeecaXqdfcNMqtx9EY4QI9oVmzZgHZC9UUwp00POOMM3zZkUce6eNOnToB8OGHH/qyUaNG+dgtwBWetHJb3q9YsWKN6haDJUuW+NgtoOT+lms67rjjAFi6dKkvO/HEE4HqK/LKkXrIIiKRUIMsIhKJ1Ics3NeUSy65xJeFsRueCNd4dVdRnXTSSb4s3MzQfd0JT4iEczvDTSYl2zHHHFPv7eEc2EqTa0t5qP7a3JjXPnDgQCD76q+dd94ZADPzZeedd56P3dVn4dfwcLNZt5BOOPzm5o1XwpDct99+6+PnnnsOqHvIYo899sj6t1KohywiEgk1yCIikUh9yMJxm5lC9oamL7/8cq37/uUvf8n6F7I3Q3TzOcMFYu666y4fu80Sv/zyyzWtdsVZa636P6PD+eCVJpyHHM6Jde+XcB3dfIWX+N5///0AjBkzxpc1NHwWPmfPnj0BmDRpki/L9fdRCdzmpLvssosvq2v4oia37RtkXzrthkLD7eJiox6yiEgkUu8hb7vttkD2J1l41V2+W6afcMIJPnYnBNxyhZC99N5tt90GZF/d9MUXXzSm2hXFLVoD2Qu85FLJO7iEu6CEyz+6hW7cnHio7jXfeeedvixXbzd8/3722Wd51WPo0KE+dpulQvUJ8GLvTBIDd4Lv3HPP9WXhe9PtKPLRRx/5Mvf7Gz16tC9bvny5j8PyXMIlTNOiHrKISCTUIIuIRMLqmnuZ885m+d+5Hv379/exO0HxzTff+LIdd9zRxwsXLiz4eU4//XQfX3rppbVu33vvvX2c79BIQ5IksYbvVVtT5bYQP/3pT3389NNP13tfN7yxcuXKotapDq8kSVL7uvkGpJnbxnCbmM6YMcOXhfOY3c4lRVqzuqDcQunyGw6tucvcGzPU+OqrrwLVw6Q1uTnfRdrAN6/8qocsIhIJNcgiIpFIZZbF0Ucf7eOuXbsC2RszrskwRSicZRGub+tWherdu3eTPE+5GzBgQNpVEOCee+4BsocpZs+e7eNTTjml5HWKSTisGcb1CZdX6NOnT63b582b5+OvvvpqDWrXNNRDFhGJRCo95FwbFLrNGiG7Bx2u/t9Y4YIr48aN87HrIYcnR6699tqCn6fc5dpoVkqjR48ePs71e/jjH//o448//rgkdSp37ls3wMknn+zj8GpgZ/z48T5evHhxcSuWB/WQRUQioQZZRCQSqQxZ/O1vf/OxW8Ql3EAyXDTorLPOArKHOR588MF6H/+WW24BYNCgQb5syJAhhVdYpAmFCziFCxq5LZ7C9/eUKVNKV7EKcfjhh/t40003rfe+4UJNMVAPWUQkEqn0kL///nsfX3/99UD2lJRwk1M3NS2corbPPvvU+/humb127dr5sjB2wmlxUrfp06f7eNWqVSnWpDK0b9/ex2Fvzk3lCjf8rYTNS0st32U6Y6QesohIJNQgi4hEIvX1kN1VeeHOCOFGm648vNLuuOOOA7K/+oU6d+4MZG9aGW4c6U4gPvzww2tU9+Zi7ty5Pl69enWKNSlvHTp0AODRRx/1ZeF79IEHHgDiWJe3nIU7soQ7jpQD9ZBFRCKhBllEJBKprIdcycpxPeQyUtbrIbutmaZOnerLwgVtfvSjHwHZQ0QlFP16yPnafffdfZxrHne44FC4ZEJj2sICaD1kEZFykvpJPZHmomXLzJ9b2BNzm6VCaj3jivPPf/7Tx+FVkeWgvGorIlLB1CCLiERCQxYiJeI20m3RokXKNZFYqYcsIhIJNcgiIpFQgywiEgk1yCIikWjsSb1FwIJiVKRC9FyDY5XbhhWaX+W2YXrvFlde+W3UpdMiIlI8GrIQEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJhBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSETbIJvZNDP7TamPbQ6U2+JSfoun0nNb9AbZzOab2dBiP8+aMLNeZjbJzL40s0VmNjbtOuVDuS2u2PNrZlub2ZSqvCZp16cxyiC3bczsKjP7yMyWmtl1Ztaq2M8bbQ+5VMysNfBP4GlgI6AbcHuqlaoQym3RrQTuBY5OuyIV6CxgO2BrYHNgAHBusZ80tQbZzNpX9ZwWVn0CTTKzbjXu1tvMZprZF2b2sJl1CI4fZGYvmNkyM3vdzIYUWJUjgY+SJLkySZKvkiT5JkmSNwp8rCgot8UVS36TJJmTJMnNwJtr8HKiEktugeHAuCRJliRJshAYB/xvgY+VtzR7yGsBtwA9gR7ACmB8jfscTiYJGwPfk0kKZtYVeAz4M9ABOA2438w61XwSM+tR9cvpUUc9BgHzzWxy1Ve/aWbWb41fXbqU2+KKJb+VKKbcWo24m5mtX8iLyluSJEX9AeYDQ/O437bA0uD/04Axwf/7At8BLYAzgQk1jp8CHBEc+5s86zeVzFe/YUBr4HRgHtC62LlRbuP+iT2/wfH/k/lTTj9nlZJbMo3680AnMsNtLwEJsHEx85LmkMUPzOwGM1tgZl8A04ENzKxFcLf3g3gB0AroSObTc0TVJ9wyM1sG7EzmE7OxVgAzkiSZnCTJd8DlwIZAnwIeKwrKbXFFlN+KE1FuLwJmAa8BLwAPkelcfFrAY+UtzSGLkcAWwMAkSdYDBleVh18TugdxDzIJWUTmFzIhSZINgp92SZKMKaAeb5D55Kskym1xxZLfShRFbpMkWZEkyYlJknRNkqQXsBh4JUmS1YW8qHyVqkFuZWZtg5+WwLpkelDLqgblz8tx3KFm1tfMfgBcAExMkmQVmTP1w81sTzNrUfWYQ3IM/ufjdmCQmQ2t+hQ+hcwv9z+FvNAUKLfFFW1+LaMtmeEgqh6rTaEvNAUx57armXWpyvEgYHQddWlaJRorSmr8/BnoQmZMZznwNnBs1W0tg/GeS4CZwBfAo0DH4HEHAs8CS4CFZAbze9QcKyLzCbrc3VZHHQ8A5lY9zzRgq2LnRbmN/yf2/AKb5Kjf/LTzViG5HVxVx6+BOcAhpciLVT25iIikrNlfGCIiEgs1yCIikVCDLCISCTXIIiKRUIMsIhKJlo25s5XZEn9pSJLEGr5XbcptXhYlSVJrXYKGKLd5KSi3oPzmKa/8qocs5WRB2hWoYMptceWVXzXIIiKRUIMsIhKJRo0hS/M0YsQIH99999313rdFixb13i4idVMPWUQkEmqQRUQiUTFDFm3btvXxZZddBsDvfvc7X2ZWPRtt1KhRAFxyySUlql1522abbXycazGqCy+8sJTVEalY6iGLiESirHvIa6+9to+ffPJJHw8YMADI7gG//fbbPh4/PrNn4qOPPurLZs+eXbR6lqNf//rXPh45cmSt22fOnOnjG2+8sSR1qnSPPPIIAMOHD/dlixYt8nH//v0B+Oijj0pbMSkZ9ZBFRCKhBllEJBJlPWRx5513+rhv374+3n///QGYMmVKzuO++eYbAHr16uXLmvOQRe/evX188cUXA7D55pv7statW/v422+/BbJP5H388cfFrmKz4E6Yrlq1ype1b9/ex/fddx8AO+20U2krViE6duwIwNSpU33Ziy++6OMTTjih5HWqST1kEZFIqEEWEYlEWQ5ZXHnllQDsu+++vuzggw/2cV1DFc7EiROB3HNqmyM3TAFw4IEHAtnztsM8nX/++QBMnjy5RLWrbIMGDfJxjx49UqxJZdp999197NqNrbbaype1bBlXE6gesohIJOL6eKhH+El38sknA/Dggw/6skmTJuX9WOFJk+bsrLPOAuCXv/xlrdvWWqv6s/qee+7x8dixY4tfsWYk/GbXr1+/FGtSnsIerjspt9tuu/myPfbYw8fuZP6sWbN8Wc+ePX3cpUsXIN153uohi4hEQg2yiEgkoh+ycJdHX3PNNb7MXbZ75JFH+rIVK1aUtF7l6pxzzqkV5zq5+dZbb/nYDW1I0xg2bJiP99xzzxRrUr5++MMfAtnDE+eddx4AHTp08GWjR4/28fXXX591LMArr7zi41133RWA22+/vQg1zo96yCIikYi+h+yW0Nx000192UEHHQTA8uXLU6lTuQkXCgp7yOEVeM78+fMB2GuvvXzZggX178+4ySab+HijjTYCYMmSJb4sXNhJsq+C3GyzzVKsSXnp3r27j2+44QYg+326bNkyAI466ihfdtttt/l49erVQPaCTR988IGP3Xs3Teohi4hEQg2yiEgkohyyCHf/OPbYYwGYO3euLwvXMZaGnXHGGT7ONUyxePFiH7s1pBsapggXJArngLuv4++//74vc/OcX3755cZUW8g+4docT1zvvPPOPv7rX//qY3e13fTp033ZYYcdBmS/93IJ25dWrVr5OFygLC3qIYuIREINsohIJKIcsnAL3ED1V+MJEyakVZ2ytd9++wENb1IabgbrFl6qi5tRccwxx/iyXDMFunXr5uOrrroKgF122SWPWleu9dZbD2jc2fyFCxf6eOjQoU1ep1i5S5rDedpbbrmlj2+66SYATj/9dF/2+eef5/XYI0aM8HG4oJMrDxfbCodKS0E9ZBGRSETZQw57dI4WtcnPPvvs4+NTTjkFyF4oyM3FhOpFVP773//W+5jhvNknnngCyF6UJeSeK3yeHXfcEchexKihnnilcL1iqO7Nhb26hoSb91a68D3l3scnnXSSL3O9YoDjjz++0Y/vrtCrK/9u55tS94pD6iGLiERCDbKISCSiHLLo37+/j5cuXQo0PC82XBfVzSfcYostfNkbb7zhYzfvNryEslK4uZhQfRItHD547733fOxOYrz++uu1Hiccpnjsscd87E6C1LXbiltrNtewU3MUnjQ6++yz8zomXNwmPHla6cL5w6eeeiqQfVIzPPns/t7DecQbb7wxkL0Y1r/+9S8ff/LJJwB89913OZ//2WefLbjuTUU9ZBGRSETZQ3YngQDeffddIPdCQm6Ff4DrrrvOx26vPbfYCGT36Nyn6qhRo3zZ+PHj17TaJbfhhhv62O2o8rOf/azeY8LdP+q7cs6dvIPce725JVAB/vCHP/j4rrvuqvf5m4s2bdoAcOaZZzb62Msuu8zHK1eubLI6xc7tBATZi4k5YS7dUplDhgzxZe7b2dNPP+3LwqU43Y4hnTt39mVhG3HzzTcXWvUmox6yiEgk1CCLiEQiyiGLOXPm+DjXYjg//vGPgeyvx25AH6o3OwxvD4cs3MaSV1xxhS97+OGHgYYXJomJ+9oG9e9y8OWXX/r46quvzuuxwzmhuU7guS3VIfv35Y4Lj3HrIb/44ot5PXclcMNi4VrUUr/w6s5cLrroIh/PmDEDyB46c+/t++67L+fxbh5yu3btfFm4HnIM70/1kEVEIqEGWUQkElEOWUybNs3HJ554IpA9r3Xy5MlA9jZB+++/v4+feuqpeh//xhtvBOC3v/2tL3MzM6699toCax0vt8gQVM/FXFPhfOeGZhL84x//ALK/HlY69x5tjN///vcAvPPOO01dnbIQDn1NmTIFyJ5dFc4fdkMV4eyoVatW1fv4ffr0AbLnLsdGPWQRkUhE2UMOr6pz8zndJyZAixYtgOzFambPnt3o53n++ed97ObvllMP2cxyxo678qiQK5DqWpDICRcxynWcO7EK2Ts9VJr111/fx48//riPBw4cWOcx4c4f4W4u7ptbQz29SuVef824qbjNkcOTeg888ECTP8+aUA9ZRCQSapBFRCIR5ZBFODzhTgSFcxTdfMRChimg+uu9m5cIcOuttxb0WGnq16+fj3PNFZ4/f369x7tLr8P5zE44TFHXQkK5uN9TU508jF04HFPfMEUo/DpeycM5MQgXGHPXH4ReeumlUlanQeohi4hEQg2yiEgkohyy+Oyzz3z86quvAtlDFu4saXhZda41TsPb3RrJUL39SzhT4JxzzlnTapfc4MGD67398MMPB7I3hwyHH9ZZZx0AttpqqzWqh7vsHKrneDbXmQISl3Bt9fbt2wPZs7hi2yJLPWQRkUhE2UMO5VpwxG3eecghh/iyO+64w8fuk3Dvvff2ZR07dvTxa6+9BsC2227ry+bNm9dENS6d8CqlnXbaqc77hSebGnOCLl/hlYBu7eRwy/p8t2evdG4z2dh6ZZUsvFbBcTsGQXzrTauHLCISCTXIIiKRiH7Iwm3KOWDAAF/mNj7t1KmTL3PDGFD9tTy8lDW8JDrcnqichZuTjhs3DoCTTjqpSR575MiROcvdAjjhesmh0047DdAwRS5vvfUWUNjCQ5K/YcOG+fiAAw7w8bfffgtkb5EVG/WQRUQiEX0PecyYMQBsv/32vqxr164AfPjhh75s+vTpPh47diyQe3v7SuJ24oDqbdPdv8WS744jzZ1bQMgt6wrZW9pL8fTq1cvH4SJZzz33HBD3N2T1kEVEIqEGWUQkEtEPWbidAbp3755yTUTy5xZneuaZZ1KuSfMT7mYTmjhxYolr0njqIYuIREINsohIJKwxl9KaWdNfd1thkiSpvZdSHpTbvLySJMl2jT1Iuc1LQbmF+PIbruU9a9YsH++www5A7oXISiCv/KqHLCISiehP6omINEY497jclG/NRUQqjBpkEZFIqEEWEYmEGmQRkUioQRYRiYQaZBGRSKhBFhGJRGPnIS8CFhSjIhUi9zYa+VFuG1ZofpXbhum9W1x55bdRl06LiEjxaMhCRCQSapBFRCKhBllEJBJqkEVEIqEGWUQkEmqQRUQioQZZRCQSapBFRCKhBllEJBL/D40KFTN0qR10AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "for i in range(8):\n",
    "  plt.subplot(2,4,i+1)\n",
    "  plt.imshow(example[i][0], cmap='gray', interpolation='none')\n",
    "  plt.title(\"Label: {}\".format(label[i]))\n",
    "  plt.xticks([])\n",
    "  plt.yticks([])\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generator model\n",
    "class Generator(torch.nn.Module):\n",
    "    def __init__(self, hidden_dims, input_dim, output_dim):\n",
    "        super(Generator, self).__init__()\n",
    "        self.lin1 = torch.nn.Linear(input_dim, hidden_dims[0])\n",
    "        self.lin2 = torch.nn.Linear(hidden_dims[0], hidden_dims[1])\n",
    "        self.lin3 = torch.nn.Linear(hidden_dims[1], hidden_dims[2])\n",
    "        self.output_layer = torch.nn.Linear(hidden_dims[2], output_dim)\n",
    "    \n",
    "    def forward(self, x):\n",
    "        h = F.leaky_relu(self.lin1(x), negative_slope=0.2)\n",
    "        h = F.leaky_relu(self.lin2(h), negative_slope=0.2)\n",
    "        h = F.leaky_relu(self.lin3(h), negative_slope=0.2)\n",
    "        out = F.tanh(self.output_layer(h))\n",
    "        return out\n",
    "\n",
    "\n",
    "\n",
    "# Discriminator model\n",
    "class Discriminator(torch.nn.Module):\n",
    "    def __init__(self, hidden_dims, n_out):\n",
    "        super(Discriminator, self).__init__()\n",
    "        input_dim = n_out\n",
    "        self.Dropout = torch.nn.Dropout(0.3)\n",
    "        self.lin1 = torch.nn.Linear(input_dim, hidden_dims[0])\n",
    "        self.lin2 = torch.nn.Linear(hidden_dims[0], hidden_dims[1])\n",
    "        self.lin3 = torch.nn.Linear(hidden_dims[1], hidden_dims[2])\n",
    "        self.output_layer = torch.nn.Linear(hidden_dims[2], 1)\n",
    "        \n",
    "        \n",
    "    def forward(self,x):\n",
    "        h = F.leaky_relu(self.lin1(x), negative_slope=0.2)\n",
    "        h = self.Dropout(h)\n",
    "        h = F.leaky_relu(self.lin2(h), negative_slope=0.2)\n",
    "        h = self.Dropout(h)\n",
    "        h = F.leaky_relu(self.lin3(h), negative_slope=0.2)\n",
    "        h = self.Dropout(h)\n",
    "        out = F.sigmoid(self.output_layer(h))\n",
    "        return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GAN Loss Definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "gen = Generator(Gn_hidden, dim_z, n_out).to(device)\n",
    "disc = Discriminator(Dn_hidden, n_out).to(device)\n",
    "##\n",
    "gen_optimizer = optim.Adam(gen.parameters(), lr=lr)\n",
    "disc_optimizer = optim.Adam(disc.parameters(), lr=lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Minimax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "disc_avg_losses = []\n",
    "gen_avg_losses = []\n",
    "num_test_samples = 5*5\n",
    "fixed_noise = torch.randn(num_test_samples, dim_z)\n",
    "bce = torch.nn.BCELoss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disc_avg_losses = []\n",
    "gen_avg_losses = []\n",
    "num_test_samples = 5*5\n",
    "fixed_noise = torch.randn(num_test_samples, dim_z)\n",
    "bce = torch.nn.BCELoss()\n",
    "##\n",
    "for epoch in range(num_epochs):\n",
    "    \n",
    "    disc_losses = []\n",
    "    gen_losses = []\n",
    "\n",
    "    # minibatch training\n",
    "    for i, (images, _) in enumerate(train_loader):\n",
    "\n",
    "        # image data\n",
    "        batch_size = images.size()[0]\n",
    "        x_ = images.view(-1, n_out).to(device)\n",
    "\n",
    "        # labels\n",
    "        y_real_ = torch.ones(batch_size, 1).to(device)\n",
    "        y_fake_ = torch.zeros(batch_size, 1).to(device)\n",
    "\n",
    "        # Train discriminator with real data\n",
    "        disc_real_decision = disc(x_)\n",
    "        disc_real_loss = bce(disc_real_decision, y_real_)\n",
    "        \n",
    "        # Train discriminator with fake data\n",
    "        zmodel = torch.randn(batch_size, dim_z).to(device)\n",
    "        gen_image = gen(zmodel)\n",
    "        disc_fake_decision = disc(gen_image)\n",
    "        disc_fake_loss = bce(disc_fake_decision, y_fake_)\n",
    "        \n",
    "        disc_loss = disc_real_loss + disc_fake_loss\n",
    "        ## Take grads\n",
    "        disc.zero_grad()\n",
    "        disc_loss.backward(retain_graph=True)\n",
    "        disc_optimizer.step()\n",
    "\n",
    "        # Train generator\n",
    "        zmodel = torch.randn(batch_size, dim_z).to(device)\n",
    "        gen_image = gen(zmodel)\n",
    "        disc_fake_decision = disc(gen_image)\n",
    "        disc_fake_loss = bce(disc_fake_decision, y_fake_)\n",
    "\n",
    "        gen_loss = -( disc_fake_loss + disc_real_loss )\n",
    "\n",
    "        # Back propagation\n",
    "        disc.zero_grad()\n",
    "        gen.zero_grad()\n",
    "        gen_loss.backward()\n",
    "        gen_optimizer.step()\n",
    "\n",
    "        # loss values\n",
    "        disc_losses.append(disc_loss.data)\n",
    "        gen_losses.append(gen_loss.data)\n",
    "        \n",
    "        if i%1000==0:\n",
    "            print('Epoch [%d/%d], Step [%d/%d], D_loss: %.4f, G_loss: %.4f'\n",
    "                  % (epoch+1, num_epochs, i+1, len(train_loader), disc_loss.data, gen_loss.data))\n",
    "\n",
    "    disc_avg_loss = torch.mean(torch.FloatTensor(disc_losses))\n",
    "    gen_avg_loss = torch.mean(torch.FloatTensor(gen_losses))\n",
    "\n",
    "    # avg loss values for plot\n",
    "    disc_avg_losses.append(disc_avg_loss)\n",
    "    gen_avg_losses.append(gen_avg_loss)\n",
    "\n",
    "\n",
    "    # Show result for fixed noise\n",
    "    plot_result(gen, fixed_noise, epoch, save=False, fig_size=(5, 5))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the discriminator is learning fast, causing the generator's loss to be zero, (and its gradientes consequently).\n",
    "\n",
    "An heuristic solution to this problems, changes the generators loss function to be\n",
    "\n",
    "$$\n",
    "- \\frac{1}{2} \\mathbb{E}_{z} \\left[ \\log(D(G(z)))\\right]\n",
    "$$\n",
    "\n",
    "Try this new loss and see if things improve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Heuristic Loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gen = Generator(Gn_hidden, dim_z, n_out).to(device)\n",
    "disc = Discriminator(Dn_hidden, n_out).to(device)\n",
    "##\n",
    "gen_optimizer = optim.Adam(gen.parameters(), lr=lr)\n",
    "disc_optimizer = optim.Adam(disc.parameters(), lr=lr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "disc_avg_losses = []\n",
    "gen_avg_losses = []\n",
    "num_test_samples = 5*5\n",
    "fixed_noise = torch.randn(num_test_samples, dim_z)\n",
    "bce = torch.nn.BCELoss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for epoch in range(num_epochs):\n",
    "    \n",
    "    disc_losses = []\n",
    "    gen_losses = []\n",
    "\n",
    "    # minibatch training\n",
    "    for i, (images, _) in enumerate(train_loader):\n",
    "\n",
    "        # image data\n",
    "        batch_size = images.size()[0]\n",
    "        x_ = images.view(-1, n_out).to(device)\n",
    "\n",
    "        # labels\n",
    "        y_real_ = torch.ones(batch_size, 1).to(device)\n",
    "        y_fake_ = torch.zeros(batch_size, 1).to(device)\n",
    "\n",
    "        # Train discriminator with real data\n",
    "        disc_real_decision = disc(x_)\n",
    "        disc_real_loss = bce(disc_real_decision, y_real_)\n",
    "        \n",
    "        # Train discriminator with fake data\n",
    "        zmodel = torch.randn(batch_size, dim_z).to(device)\n",
    "        gen_image = gen(zmodel)\n",
    "        disc_fake_decision = disc(gen_image)\n",
    "        disc_fake_loss = bce(disc_fake_decision, y_fake_)\n",
    "        \n",
    "        disc_loss = disc_real_loss + disc_fake_loss\n",
    "        ## Take grads\n",
    "        disc.zero_grad()\n",
    "        disc_loss.backward()\n",
    "        disc_optimizer.step()\n",
    "\n",
    "        # Train generator\n",
    "        zmodel = torch.randn(batch_size, dim_z).to(device)\n",
    "        gen_image = gen(zmodel)\n",
    "\n",
    "        disc_fake_decision = disc(gen_image)\n",
    "        gen_loss = bce(disc_fake_decision, y_real_)\n",
    "\n",
    "        # Back propagation\n",
    "        disc.zero_grad()\n",
    "        gen.zero_grad()\n",
    "        gen_loss.backward()\n",
    "        gen_optimizer.step()\n",
    "\n",
    "        # loss values\n",
    "        disc_losses.append(disc_loss.data)\n",
    "        gen_losses.append(gen_loss.data)\n",
    "        \n",
    "        if i%1000==0:\n",
    "            print('Epoch [%d/%d], Step [%d/%d], D_loss: %.4f, G_loss: %.4f'\n",
    "                  % (epoch+1, num_epochs, i+1, len(train_loader), disc_loss.data, gen_loss.data))\n",
    "\n",
    "    disc_avg_loss = torch.mean(torch.FloatTensor(disc_losses))\n",
    "    gen_avg_loss = torch.mean(torch.FloatTensor(gen_losses))\n",
    "\n",
    "    # avg loss values for plot\n",
    "    disc_avg_losses.append(disc_avg_loss)\n",
    "    gen_avg_losses.append(gen_avg_loss)\n",
    "\n",
    "\n",
    "    # Show result for fixed noise\n",
    "    plot_result(gen, fixed_noise, epoch, save=False, fig_size=(5, 5))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "torch.save(gen.state_dict(), 'gen.pt')\n",
    "torch.save(disc.state_dict(), 'disc.pt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We load some pretrained weights\"\n",
    "gen.load_state_dict(torch.load('gen.pt'))\n",
    "disc.load_state_dict(torch.load('disc.pt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolating in the latent space\n",
    "\n",
    "Given two random samples $z_1, z_2 \\in \\mathbb{R}^d$ with $d=100$ as in the previous code, we can interpolate in this latent space and generate the corresponding images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "torch.manual_seed(16)\n",
    "z1 = torch.randn(dim_z).to(device)\n",
    "z2 = torch.randn(dim_z).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "alpha = np.linspace(0, 1, 9)\n",
    "alpha = np.atleast_2d(alpha).T\n",
    "alpha = torch.tensor(alpha, dtype=torch.float).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "zs = alpha*z1 + (1-alpha)*z2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/dist-packages/torch/nn/functional.py:1339: UserWarning: nn.functional.tanh is deprecated. Use torch.tanh instead.\n",
      "  warnings.warn(\"nn.functional.tanh is deprecated. Use torch.tanh instead.\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAFBCAYAAAD9m4PQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XmQVNXZP/DnRgXZl2ER2RFQkCUIKIsaZFFCAhEThF+9BhdilNUQIQlQUYkRiEtpSRT8aUJpRK0kxCBBE5VFooIKCAQQkFX2dYZFlmG57x/v8OR7jn2bnu7b3bdPfz9Vlk9P9+m5faY5dZ5zz+L5vi9ERC75VrYvgIgobGzYiMg5bNiIyDls2IjIOWzYiMg5bNiIyDls2IjIOWzYiMg5F2fyl3mex9nAGeL7vifCOs8k1nnmna9zG3tsROQcNmxE5Bw2bETkHDZsROQcNmxE5Bw2bETkHDZsROQcNmxE5Bw2bETkHDZsROQcNmxE5Bw2bETkHDZsROQcNmxE5Bw2bETknIzux5YvNm3aZDyeO3euxg899JDGRUVFGbsmF3nef7fi+uMf/2g899e//lXjTz75ROMDBw6k/8Ic9q1v/bcv1K9fP+O55cuXa3zy5EmNCwsLNT5z5oxRJl0HtrPHRkTOYcNGRM7x0tUVjPnLHNsyedmyZRpfc801CZXBrnw6697Vbaq3bt2qccOGDRMqU7ZsWY2Li4vDviTlUp1jne3Zs0fjqlWrBpY5ffq0xvi3wfJhf+e5NTgR5Q02bETkHDZsROQcTvcoJRwja9u2banLFxQUaMypB7HhNA4Rka5du2rcoEGDhN4Dx3LatWunMU79yDb8nOkcb8XvbKNGjTTeuXOnxjfffLNRBqfPxBtXQ/h5Ro8erfEjjzyi8fHjxxN6r1Sxx0ZEzmHDRkTO4XSPUipTpozGK1eu1Piqq64KLINTDPA2ejrlwtQDTIu++uorje1UtG7duhovXLhQ48aNGwe+94kTJzSuWLGixufOnUvmUhMS1TqvUaOGxp9//rnGo0aN0njbtm1Gmc6dO2s8ZswYjfFvZjt27JjGmP5v3LixdBdcCpzuQUR5gw0bETmHd0VLdOjQQeMVK1ZobKcubdq00TheKoQylX7mAkyLpk+frvFPfvITjffv32+UOXz4sMaJppK1a9cudRlX2Kn8jBkzNMY7nDh8Mn/+fKMMDrNgKv+b3/xGYxyWERF58sknNcahgGxgj42InMOGjYick9ep6CWXXKJx+fLlNT579qzGI0aMMMq0bNlS46AUM5N3mqPuoosuMh5Xr15d4x07dmi8d+9ejdu3b2+UKVeunMaJLnw/evRoqa7TJb179zYeb9++XWOcFP7CCy9ojN9/EXNSLw4T4L8Z28yZMzXetWtXKa44fOyxEZFz2LARkXPYsBGRc7jyoASOBV155ZUa2wuAP/roo5jlcS93vD0uInLq1KkwLrFUojIL3p56gI/r16+v8a233qqxXceLFy/W+OKL/zssjN/d22+/3SiDZx5kSlTqvGnTpsZjrLP7779f4y5dumg8bNgwo8yf/vQnjYNW1dhne9i/t7SS2RSAKw+IKG+wYSMi5+T1dA+EUzxw0foHH3yQUHlcUJyN1DOq7JQCH+PCa5xu8OGHHxplMJVCX375pcbZSD2jyl50jsMszzzzjMY/+MEPNF60aJFR5tJLL4353vjvpHnz5ildpy3MYTH22IjIOWzYiMg5TEVjwCPegtIgEZF169ZpPG3atHRekpNwdjvuExZv0wA8yg3vXlMwTB9//OMfa4wrOjC24XbeeMc/yits2GMjIuewYSMi5zAVjWHixIkJva5///5pvhK34T5pDz74YEJl7rrrrjRdTX6YNGmSxvYGD0HwlKkop5+IPTYicg4bNiJyDhs2InIOF8GXwNnZuHLA3igR2Qu8oyQqC7Ljwak0u3fv1hjPRbCxzlODU2neeecdjbt162a8DusZV4Vk+ywDGxfBE1HeYMNGRM7hdI8SeJRevPQT94+n1NSpU0djTD/t4RFuKpAaTCvxWMJOnTppHG+zgpMnT6bx6tKDPTYicg4bNiJyDlPREhs2bIj5c7sb3qBBg0xcjrNw4bu9tXQQ+2g4+iZMN+0T2itVqqTxqlWrNI638P22227TOFdWGyD22IjIOWzYiMg5eZWK2pM7jxw5Evjcea1bt07rNbnO3ltt3759GgedKv7cc88Zj3MxFco0rKMOHToYz7399tsaY1qKcIK0iMi8efNCvLrMY4+NiJzDho2InMOGjYick1djbPb5BfaJ7bFEbdFvrrHH2ILGeFBhYWG6LicvHDt2zHgc7wyJ806fPm08Pnr0aKjXlGnssRGRc9iwEZFz8ioV7dOnj/EY99zHGfF4XNnOnTvTf2EOa9GihfG4qKhI42rVqmmMdY577FPpYV2KmNNnRo8erXFxcbHGI0eONMrk+hQb9tiIyDls2IjIOXmVitasWdN4jKsNsOsdb2tqKh07pcEF7Zgy2cMEVDr4Xbbvivbo0UNjXG2zYMECjTdv3pzGq8s89tiIyDls2IjIOc6noi1bttR448aNxnNBC9+x6z5r1qz0XJjDmjZtqvGZM2eM54Imiw4ZMkTjd999Nz0X5hi8k4+bNVSvXt14Hf4bwI0HunbtqvGaNWvScYlZwx4bETmHDRsROYcNGxE5x8kxNhx7mDx5ssb9+vVLqHzQBogUDPfPnzJlisa4d348QWdOULBWrVppPH36dI2vu+66wDI4/WbJkiUxf+4C9tiIyDls2IjIOV4mu6Ce52W8v4tp6RdffGE8d8UVV2iMs+AT2b8q6nzf90SyU+dVqlTR2J7RXqFCBY1xD7BE9mmLumzW+S233KLx7NmzjedwsTseJ1mnTh2N7YXzueJ8ndvYYyMi57BhIyLnOHlXFF100UUa23urPfzwwxrv2LEjY9fkOlx58NVXXxnPPfvssxrv378/Y9cUdUEbMsR7HQ6zDBw4UONPP/3UKIMzA3Cr+1xNPxPBHhsROYcNGxE5x/m7oqigoMB4fPDgwSxdSfpl8w4dqlWrlvEYT4J3TVTqvHLlysZj3IPNNbwrSkR5gw0bETmHDRsROSevxtjySVTGe/IJ6zzzOMZGRHmDDRsROcf5lQeUmERnvlN4cOXAuXPnsngl0ZPq95E9NiJyDhs2InIOU9ELsI/oczVNi9LnwhRNxN00LUqfKwp1HuZwCHtsROQcNmxE5BymohcQpRQtV6SaUkQpRcsVuO9gMvusRaHOw/y3xh4bETmHDRsROYcNGxE5h2NsOSDXVgXkwjVeSK6tCnDh/AJO9yAiioMNGxE5h6loDnAhtUuHdKaLqbxfrg0dlEbNmjU1Dvv4RE73ICKKgw0bETmHW4OHpEKFCho3a9bMeC7olHmcLX748GHjueLi4phlEk2RXNqmuly5chqfPHlS465du2p85swZo8z69es1Pn36tMaNGjXSeOPGjUaZSy65ROPjx49rjHcc8W9mPxfVOg9K2eOl8nfffbfGr776qsazZs3S+MMPPzTKvPzyyxoXFRVp3KJFC43Xrl1rlGncuLHG+PfAeq1Ro4ZR5sCBAxpza3Aiyhts2IjIOUxFY8C7WhdfbN44Hj58uMaPPvqoxuXLl49ZXsRMM++44w6NFyxYoPGpU6eMMpgaJPM3impaFKR69eoa33nnncZz48eP17hixYoaY+po7yd26NAhjefMmaPxhAkTND569KhRBtNcTF8TlQt1jvV07733aoz1IiJSq1YtjTEtLFu2bMz3EjHrb8+ePRp3795dY0xR7TIYJ4qpKBHlDTZsROQcNmxE5JxIjrGla0a5PfbVoEEDjZ9//nmNd+3apfGtt95qlMGxIHuMIQhORdi0aZPGv/vd7zSeMWNGQu+VqFjjPfFmxIc5Wz7eGOXUqVM1HjRokMY4vaVq1apGGXwP+28YBD8DjpetXr1a42uvvdYok+pC8tLWeZhwGgqOQ4qIbNmyRePKlStrjN9LexoLfrfxMyRT/3v37tUY/82JJDeWaf0ejrERUX5gw0ZEzolkKpoudlr01ltvady7d2+NseuOrxERqV27tsY4o/rb3/62xnadHjt2TOPZs2drPHjw4MAyqbrQ1IOwj1sLSssvvfRS4zFOA6hUqZLGOPXi/fffN8pgylSlShWNb7rpJo3tNBJTJkx3cBa7PcUGp49gmUTTyQvVebxVC8nAz4LXZdf5tm3bYl4Dfv6//OUvRhmclnH11Vdr3LFjR43tf084FQRVq1ZNY3uFDa7Y+frrrzUubZ3b2GMjIuewYSMi5+TVfmwFBQXG427dummM3V1MRadMmWKU2bdvn8a4gPell17S+Dvf+Y5RBlPRyy+/PObvzLR0bneN743pooiZeuDnX7ZsmcZDhw41yuBdNbxj2rlzZ43//Oc/G2XKlCmjMabJdvqJgu7QhfV3Cnv7brxeTPFGjhxpvC5omOCdd97R+MEHHzSew+85ppI4FLNkyRKjDKaiWGd2+okw/UTcGpyIyMKGjYic43wqinduPvvsM+M5vEOHXez77rtPY0wj7ffDLn+nTp00trv+u3fv1hgXwceTa9tL4zXihgBvvvmm8brCwkKN8a7o6NGjY75GxLyTh/uE4SJu/J0i5oLqn/70pxf+ABlgT25N9e965MgRjTENfPrpp43X4f5lmFY+9dRTGtt1jneIn3jiCY1xJoBd55hq20M4mcYeGxE5hw0bETmHDRsROcf5lQc4JoYbO4qIHDx4UOM2bdpojLPb+/bta5TBqRyTJk2K+Ttxv3y7DE5rSLTukxlvy+amhzg+g+OYIub149QPHJe0Nx7AI99wX3104sQJ4zHux//9739f46CzJJKFf5tz585lrc6x/uzxMlyJYK8WOG/IkCHGY9xoEjdrQPbxe5s3b9b4+uuv19g+jyJV+FnPnj3LlQdElB/YsBGRc5xPRZHdDe/fv7/GeEsbVwfYt+jxsb2o+TycES/yzRnamZDNVBTrqG7dusZzc+fO1bh169Yxy9iCjoxD06ZNMx4PGzYssYsNUVTOPLCHT1577TWN7b3aUmHXsf03yAQugieivMGGjYick1epqA1TU1xhgAuo46VIONMat1y274pmQ1TSIrv+brjhBo0/+OCDUr8ffl+bNGmi8datW0t/cSGIdYcu23WO318RkY8++kjjDh06pPTeuPJg3bp1Kb1XsvCue3FxMVNRIsoPbNiIyDl5nYoiPD0HT5IKmtAoYk4Wveuuu9JyXcmKSioaz8CBAzV+4403EiqDE28xrY2CXKjzxx57TOPx48cnVAYn4uJmEVHYnIF3RYkob7BhIyLnsGEjIuc4v9FkonBB/Jw5czS2F2Tj9IX58+en/8Ik9zadTNTbb7+t8bx58zTu0aNHYJmg1R5hCzqKL9fhxg233367xk2bNg0sg9NH0vn9CzqOMRnssRGRc9iwEZFz8joVxbSmfv36Grdr105j+5g6LHPo0KE0Xt1/uZR+4vQZrHOsZ3v/LiyD0w3SyaX0E+sPNyXAqTP16tUzyuAebvbJ8umSavqJ2GMjIuewYSMi5ziZiuJdRNzm214A3KxZM43HjRuncZ06dTS278LhqdYbN25M/WIdgXXetm1bjX/xi18Yr2vevLnGeHxb9erVNbZXe2CaiieH4wL0dJ5sH1VY5z179tT4oYceMl6Hdb53716Nd+3apbG9cB7rE9Py9u3ba4zb3KdTMscWssdGRM5hw0ZEznEyFUV4IjguABYxJyXinR9MceyT4Ldv364xduUzJaqTdbH+qlatqrE9wbmoqEjjnTt3aozpp53+4+nl+Lp8T0XxrvKNN96ocZcuXYzX7du3T+PZs2drvHTpUo13795tlBk8eLDG+LetUaOGxmGfbB8kmfdlj42InMOGjYicw4aNiJzj5Bgb5uQ4JrBw4ULjdXjMHu6Zj2MX5cqVM8o888wzGuPUg0wJaxzDHsfC8xsShWMsON41ceJEje0pMTiuNnz4cI1xSsHy5cuNMjjOuWLFipSuOZvatGljPF61alWp3wPrGY8vHDNmjMYHDhwwyrz33nsaP/LIIxpj/d1zzz1GGXwOxy9PnTqlcZTGeG3ssRGRc9iwEZFznD/zAFOujh07Gs/hdI0ZM2Zo3Lhx45jlRczu/5EjR0K7zrBlc/99XFDdq1cv47lXX31V47Jly2rcvXt3jV9//XWjzMcff6zxrFmzNH7hhRc0jkJalOk6x7R0wIABGv/oRz8yXodnS2BaifWPK2pEzPrE34OrGLZt25bMZYeKZx4QUd5gw0ZEznHyrijCO3eFhYXGc3369NEYT2/H2e149JjIN1ciRF0yKxUSLYOvw620cV+ttWvXGmXwJHFcoYBDAfYieLz7jLPloyqddY51gxsHrF+/XmP7TvRNN92k8YYNGzTGk9ztRfB4DcXFxTHjKGOPjYicw4aNiJzjzF1R7MoXFBRoPGLECI3tbaUHDRqkMaZFuDW1vYfbypUrU7/YDAjrDh3Wq50iVqhQQWOsywkTJsR8jf0eeCoRsocMmjRpojGmpVHbvjsddW5vy33llVdq/PLLL2vcqlWrmOVFzDuhiZ7yhZOHcVL6K6+8onEUJkjzrigR5Q02bETkHDZsROQcZ8bYcLoBnuTeuXNnjStXrmyUwc+O8ZYtWzTGmdYiubOhYTpmweMMdBGRyy67TGOcboBnGdhlENYlbgiKqztERDZv3lz6i82CdNS5PV6G50ngmQPx6hnh9xzHyF588UXjdSNHjtQ46DsfpdUeNvbYiMg5bNiIyDk5m4rG208MpxRMmTJF4wceeMAog/vv44J2vHV+4sSJ1C82C9KRFtmpPK4wwFnwkydP1vjuu+82yuB+bHh+xHe/+12Nc211x3npqHM870HE/J7ivoG9e/fW+IknnjDK4MYB//nPfzSeOnWqxplaURD2OQlMRYkob7BhIyLnRD4Vbdiwocbx9n8Kmq39t7/9TePFixcbZV566SWNcbUBHleWq1JJiyZNmqTx+PHjA18XtAh+7NixGs+cOdMog+knfvdy5W5zPKnUOaaI9l3hIHgnFNNSrGMRN+o2CFNRIsobbNiIyDmRT0XDhOmSSPQWUYcpm1uD5yvWeeYxFSWivMGGjYicw4aNiJyTU2ceJLOXPHJhTC3VOigtXOERhY0FswGPqcOT0NOFdW5u/IorhBLFHhsROYcNGxE5J5LTPTKdbuWSROumtFMPWOfBcIZ/vFn8nO4RnkTTcU73IKK8EcmbB1HqMYS9zUqq15Cu3x+lOo8CrHOX11pGCdZ5qjdN2GMjIuewYSMi50QyFQ1TogO/QaKQokXhGkrDhRsRuXrduSzMOmePjYicw4aNiJzDho2InOP8GBtv1Wcex6co29hjIyLnsGEjIuc4n4omKpkpCljm6quv1nj16tXhXZjDcD1goidWYZ1fd911Gi9ZsiTkq3MTnuCGU6GOHz+eUPmWLVtqvHbt2vAuLGTssRGRc9iwEZFzIrltUabUrFlT40OHDmlcUFCg8cGDB40y//73vzXu1auXxsuWLdP4rbfeMspMnjxZ4yNHjmjcqFEjjbdu3WqUufHGGzVetGiRxrg4GH+/iMh7772ncVS30Bk8eLDGf//73zUeOHCgxq+99ppRBg/Kbt68ucZ4yPCsWbOMMuPGjdMY06waNWpofODAAaPMzTffrDHWJabGU6dONcqMHDlS46jW+Zw5czQeM2aMxlj/ffv2NcqsWbNG4yZNmmi8dOlSjfHAcRGRRx99VOPi4mKNy5cvr7Gd8rZq1UpjTG2xzvFaRMxhH25bRER5gw0bETknr1NRvCuE6ee7776rMXZ7bSdPntQY7zZdfLF5sxm75XigDL43pqgiImfOnNH42LFjgdcQJJtpUbw7zFWqVNG4f//+Gj/77LMaV6hQwSiD6TfWebly5TTGO6wiZv2hyy+/XON4h/scPnw48Lkg2axzPAzc/uyDBg3S+KGHHtIY03r8tyBi1g0eYIPfc7vO8e8U9G8r3j5rR48eDXwuCFNRIsobbNiIyDls2IjIOVkbY8vUZoQ4JnDVVVcZz+EUjaANKe1Z8PYZCOfZ4w2JKCws1LhOnTrGczgul0z9ZHO8Bz8LTu8QEZkyZcoFy8cbh8H6t8eFEoFTPOrXr288h2NJuVbnvXv31tiu47Zt216wfLzvOdZFMnWOU5nsf4Nhfc9t7LERkXPYsBGRc7K2CD5e9zbVPdSGDh2qcZkyZTRev3698bqgbjX+/jfeeMN47ssvv9R4wIABGuMqAkx/Rb45/eO8Fi1aaIxpkIg5LcKeCnJeKil82McK4oz0SpUqafz+++8nVB7r/M033zSeW7lypcZjx47VGGe020MBQUMG7dq109iu84oVK2qMU2yiur8cfi/wO2evYgmCn+uzzz4znlu4cKHGuFoBywTVse2GG27QGFNPEfO6T5w4kdD7JYI9NiJyDhs2InJOJPZjC3v7buyK//Of/9QYZ7fHg6kQpj4iInv37tX4lVde0RjTzc8//9wogylO0HvZkpn5Xhphp1crVqzQ+L777tP4H//4R0LlceH7r371K+O53bt3a4x3svft26exvR9bUPq/Y8eOwGtIZuZ7Nm3atEnja6+9VuPly5cnVP7111/X+Ne//rXxHNYTpvn/+te/NMZ/W/brgt7LFmb6idhjIyLnsGEjIudEIhUN+64oLmLHxcHDhg0zXod7sOFdSEyL7D278P2mT58es3y8z/Pcc89d+ANkgJ2qBS0ajwfvit1///0a4+fHPc5EzBQb755iioOpp4iZ4jz22GMxf46L40XMhfRRqXP87ojEX4QfBOsW00/8jtkbN2C6h3chcW/BLVu2GGXw+/G9731P4549e2pcVFRklKlWrZrGDz/8cJxPkX7ssRGRc9iwEZFz2LARkXOc32gSxyTsKRQ4cx1fh3WCUxdERGrVqqXxb3/725i/84svvjAe42L366+/PubvCQOOOZ05cyZrC7KxLr/++mvjOXtVRix2nderV09je1rCeTgNRMQcM+zUqdMFf2eyoljn9nijPbYXywMPPGA8xpU0o0ePjlnGXiGC/55wUX7Y8LOePXuWi+CJKD+wYSMi5zifiqLnn3/eeIwpTzL7TAXp3Lmz8Tgbp5RH5Sg4PCJPxDxuLUzNmjUzHm/cuDEtvyeeqNS5fcQdng0RJjzKUOSbR1VmAvdjI6K8wYaNiJwTiZUHmWJvmXzvvfdqnGoqetttt2mcjdRTJLE7jpk2c+ZM4/HkyZNDe288yTwbqaeIeScwKuzjGsNMRTHdzEbqKZLY52GPjYicw4aNiJyTV3dFbThxcd68eRrjVsbx4CRQTEmSWdwctqjcobPhKUWLFy/WuGrVqim9l73tezZEtc4//fRTjdu3b69xMsMvuL9e3759U7uwEPCuKBHlDTZsROQcNmxE5Jy8HmNDuDHhmjVrNG7YsGFgGTyxPGiP/TA0adJE482bNydUJqrjPQhv2+MKhSuuuCKwDH5fcQF62N9jXMRt7+0fJBfqHKcE4ZF7ia4Iwc1B7WklqRo1apTGiZ5PwjE2IsobbNiIyDl5nYri7e46depo3K1bN43tmfL169fXGPeZx1Q0CieHRzUtwnMSKleurDEuqP7lL39plMEVIljnmMraJ4zj3zbs4x2DRLXOUZkyZTTGOsL910TM4RisPxwmsPc3xI0Ili5dmvK1JoKpKBHlDTZsROQc5xfB41bcw4cPN57r3r27xqtWrdL4D3/4g8Z4pJiIudoAt2B+/PHHNR43bpxRBu+ehpmmhn1sYViuueYajW+55RbjuREjRmi8cOFCjSdMmKCxfXwcntCOdY4nmdvHve3Zs0dj+wjFVIRxbGFYMK3HIZLLLrvMeN2LL76oMa4cwKMMGzRoYJT5+OOPNcb93Z566imNa9asaZTZv3+/xgMGDLjwB0gQps8i3xx2iIU9NiJyDhs2InKOk3dF8VSie+65R+MhQ4YYr1u9erXGP/vZzzTGOzp2F33lypUaYyowd+5cjX/+858bZVzdpho//w9/+EON77jjDo379etnlJk/f77GWOd4spd9qhLu+4UTTDds2KDxggULjDLDhg278AcIWSbqHIcf7rzzTo1xP8A+ffoYZZ5++mmNp02bpnG8yd7r1q3TuHHjxhpjim+niDizIFPtCu+KElHeYMNGRM5hw0ZEzonEdA+cXiHyzbGwROCCaBy7wdvO+/btM8rgnvn2GM15ePK7iHl7H8eC8HWbNm1K9LKzZvv27cZjnC6QKKzzJ598UmPcOGDv3r1GGRy/xDFOFG+KDY7r1a5dW+MxY8YketlZYx+Ll8x5Cfidw4XiFStW1HjXrl1GGazPoO+mPV6G01rw74xTPHr16mWUicKKm/PYYyMi57BhIyLn5Ox0D0xJRMxuMHbXMUXq2bOnUaZ169Ya46x9fG97zynssmN3HU+Zt1c4ZEM6ph7EW+mAC9pxdnubNm2MMrjYPWjWvr24GvfKw7Ro7dq1GturFbIhHXUeb6VDvXr1NMahFHtBe0FBgcZHjhyJ+XvsIYPq1avHvIZDhw7FfN9s4XQPIsobbNiIyDmRuCuaDDuFxvQTZ0rj4nZc9C4iMnHiRI1///vfa4x3DO1Z8EFw5ryr7EX2mKLgCo0tW7Zo3KFDB6MMpo84Qx6Pz0v0WDhceeAqO13HusEjC4uKijTGdF3EXJzeo0cPjRctWqSxPbQTJMwNBdKJPTYicg4bNiJyTk7dFcXusp0idu3aVeOZM2dqjAtzw4DpGC4ixjusp06dMspkY+JiOu7Q2SlOx44dNcbJzjhxNmw4yRX3HbMnv+IeeJmSjjq3U0TcfhvvPuPPw4bpcN26dTW2ZwycOHFCYy6CJyIKGRs2InIOGzYick5OTffAvN3e9xxnq6c6robjaDjTWsScSmKP67jOHrfCWeypjqvh39b+PTiVBKcuIHuGfjbG2NLBHqvC732Y42r27+nSpYvGuLkq/tvAlQ8i8TeuzDT22IjIOWzYiMg5OTXdA+HCaBFzigWeVl2lShWN8bg2EZFRo0ZpXFhYqPEnn3yicaZSmrCP0kvH1ANNBnqIAAADfklEQVR7zy68RtyDDYcFpkyZYpQZO3asxjh9BM81yNR30p6+kurfOhMbD+BjTP8xRcepTyLmCgNM2e1pSZkQ9vGFnO5BRHmDDRsROSfyqejOnTs1xlnP8WB3HY9ry6e7mKmkRbigGhdax4N1julGIqd2uyKVOsc7zLi3XTx2Kn2eK3eEE8FUlIjyBhs2InJO5FNRSk4mTiUnE+s885iKElHeYMNGRM5hw0ZEzsmpRfA4pSDVmfm5Ck/1xtUS6YIbHUbppO9MwukXQcfXUbjKli2rcTIrJNhjIyLnsGEjIudkdLoHEVEmsMdGRM5hw0ZEzmHDRkTOYcNGRM5hw0ZEzmHDRkTOYcNGRM5hw0ZEzmHDRkTOYcNGRM5hw0ZEzmHDRkTOYcNGRM5hw0ZEzmHDRkTOYcNGRM5hw0ZEzmHDRkTOYcNGRM5hw0ZEzmHDRkTOYcNGRM5hw0ZJ8TzvrOd5K+C/X4X43o08z1udwOt+7nneWs/zVnmeN8/zvIZhXQPltouzfQGUs074vv/tLF/D5yLSwff9457nDRWRx0VkYJaviSKAPTYKled5Wz3Pe9zzvP94nvep53lNS37eyPO8+dC7alDy89qe573ped7Kkv+6lLzVRZ7nveh53hrP8971PK+c/bt831/g+/7xkodLRKReRj4kRR4bNkpWOSsVxZ7SYd/3W4vI70XkmZKfTRWRl33fbyMiM0Xk2ZKfPysiH/i+31ZErhGRNSU/byYiz/m+f7WIFInIDy9wPUNE5J2UPxU5wfN9P9vXQDnI87xjvu9XjPHzrSLS3ff9zZ7nXSIie3zfL/A874CI1PF9/3TJz3f7vl/D87z9IlLP9/1T8B6NROQ93/eblTz+pYhc4vv+bwOu5Q4RGSEi38H3ofzFMTZKBz8gLg1soM6KyDdSURERz/N6isgEYaNGgKkopcNA+P/ikvhjERlUEv+PiPy7JJ4nIkNFRDzPu8jzvCqJ/hLP89qJyAsi0s/3/X2pXjS5gz02SlY5z/NWwON/+r5/fspHNc/zVsn/9br+X8nPRorIDM/zxorIfhG5u+TnD4jI//c8b4j8X89sqIjsTvAanhCRiiLyF8/zRES+8n2/X7IfiNzBMTYKVckYWwff9w9k+1oofzEVJSLnsMdGRM5hj42InMOGjYicw4aNiJzDho2InMOGjYicw4aNiJzzvwLuc/xPkS/AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 9 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_result(gen, zs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}