aevb.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced PyMC3: What the heck is ADVI?\n",
    "![PyMC3][logo]\n",
    "[logo]: https://raw.githubusercontent.com/pymc-devs/pymc3/master/docs/pymc3_logo.jpg \"Logo\"\n",
    "\n",
    "Peadar Coyle - Contributor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why Variational inference? \n",
    "* Scaling probabilistic modelling to increasingly compex problems on increasingly large datasets. \n",
    "\n",
    "**PyMC3**: Variational inference: ADVI for fast approximate posterior estimation as well as mini-batch ADVI for large data sets.\n",
    "\n",
    "**Uses** \n",
    "* Large scale topic models (LDA, Hierarchical Dirichlet Process)\n",
    "* Semi-supervised classification\n",
    "* Generative models of images (used in Deep Learning work)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Limitations**\n",
    "* One limitation: Variational Inference requires that intractable distributions be approximated by a class of known probability distributions. \n",
    "\n",
    "**Objection**\n",
    "Unlike MCMC, even in the asymptotic regime we are unable to recover the true posterior distribution.\n",
    "*There are solutions* \n",
    "We can use a method called **normalized flows** - this is out of the scope of this talk. There is work being done on this in PyMC3.\n",
    "\n",
    "**We'll assume that these limitations don't apply to our case** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "%matplotlib inline\n",
    "import sys, os\n",
    "sys.path.insert(0, os.path.expanduser('~/Code/Advanced_PyMC3_Talk'))\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.animation import ArtistAnimation\n",
    "from matplotlib.patches import Ellipse\n",
    "from collections import OrderedDict\n",
    "from copy import deepcopy\n",
    "import numpy as np\n",
    "from time import time\n",
    "from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import theano\n",
    "from theano import shared\n",
    "import theano.tensor as tt\n",
    "from theano.sandbox.rng_mrg import MRG_RandomStreams\n",
    "\n",
    "import pymc3 as pm\n",
    "from pymc3 import Dirichlet\n",
    "from pymc3.distributions.transforms import t_stick_breaking\n",
    "from pymc3.variational import advi, advi_minibatch, sample_vp\n",
    "\n",
    "theano.config.compute_test_value = 'ignore'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# configure pyplot for readability when rendered as a slideshow and projected\n",
    "plt.rc('figure', figsize=(8, 6))\n",
    "\n",
    "LABELSIZE = 14\n",
    "plt.rc('axes', labelsize=LABELSIZE)\n",
    "plt.rc('axes', titlesize=LABELSIZE)\n",
    "plt.rc('figure', titlesize=LABELSIZE)\n",
    "plt.rc('legend', fontsize=LABELSIZE)\n",
    "plt.rc('xtick', labelsize=LABELSIZE)\n",
    "plt.rc('ytick', labelsize=LABELSIZE)\n",
    "\n",
    "plt.rc('animation', writer='avconv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "blue, green, red, purple, gold, teal = sns.color_palette()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior distribution is equal to the joint distribution divided by the marginal distribution of the evidence.\n",
    "$$\\color{red}{P(\\theta\\ |\\ \\mathcal{D})}\n",
    "    = \\frac{\\color{blue}{P(\\mathcal{D}\\ |\\ \\theta)\\ P(\\theta)}}{\\color{green}{P(\\mathcal{D})}}\n",
    "    = \\frac{\\color{blue}{P(\\mathcal{D}, \\theta)}}{\\color{green}{\\int P(\\mathcal{D}\\ |\\ \\theta)\\ P(\\theta)\\ d\\theta}}$$\n",
    "For many useful models the marginal distribution of the evidence is hard or impossible to calculate analytically.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variational Inference\n",
    "* Choose a class of approximating distributions\n",
    "* Find the best approximation to the true posterior\n",
    "\n",
    "Variational inference minimizes the [Kullback-Leibler divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence)\n",
    "\n",
    "$$\\mathbb{KL}(\\color{purple}{q(\\theta)} \\parallel \\color{red}{p(\\theta\\ |\\ \\mathcal{D})}) = \\mathbb{E}_q\\left(\\log\\left(\\frac{\\color{purple}{q(\\theta)}}{\\color{red}{p(\\theta\\ |\\ \\mathcal{D})}}\\right)\\right)$$\n",
    "\n",
    "from approximate distributons, but we can't calculate the true posterior distribution.\n",
    "\n",
    "Minimizing the Kullback-Leibler divergence\n",
    "$$\\mathbb{KL}(\\color{purple}{q(\\theta)} \\parallel \\color{red}{p(\\theta\\ |\\ \\mathcal{D})}) = -(\\underbrace{\\mathbb{E}_q(\\log \\color{blue}{p(\\mathcal{D}, \\theta))} - \\mathbb{E}_q(\\color{purple}{\\log q(\\theta)})}_{\\color{orange}{\\textrm{ELBO}}}) + \\log \\color{green}{p(\\mathcal{D})}$$\n",
    "is equivalent to maximizing the Evidence Lower BOund (ELBO), which only requires calculating the joint distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example of Variational Inference\n",
    "In this example, we minimize the Kullback-Leibler divergence between a full-rank covariance Gaussian distribution and a diagonal covariance Gaussian distribution.\n",
    "\n",
    "* Thanks to Austin Rochford and Taku Yoshioka for these examples\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SIGMA_X = 1.\n",
    "SIGMA_Y = np.sqrt(0.5)\n",
    "CORR_COEF = 0.75\n",
    "\n",
    "true_cov = np.array([[SIGMA_X**2, CORR_COEF * SIGMA_X * SIGMA_Y],\n",
    "                     [CORR_COEF * SIGMA_X * SIGMA_Y, SIGMA_Y**2]])\n",
    "true_precision = np.linalg.inv(true_cov)\n",
    "\n",
    "approx_sigma_x, approx_sigma_y = 1. / np.sqrt(np.diag(true_precision))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApMAAAHBCAYAAAA1uDCcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3WlsXOd97/HfObNyhsPhNuIqUaIWarFWL9exEydIGwRF\njKQ12jpo6l6jDYI0RQMHdYE4MYLAceC+cJq2TtOmLepcFEHdBnbRi7YGnAJpHMexHcWyFkuULEqU\nSIqkhsuQwxnOfu4LybyWTFHU4fLM8v0ACSVzLP55og6/Pec8z7Ecx3EEAAAAuGCbHgAAAACVi5gE\nAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALi2ZEyyaxAAAACW4l3qk5ZlKR5PrtcsFSkWi3CM\nlsDxuTmO0dI4PjfHMVpaLBYxPQJQ1bjMDQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgE\nAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMC1JZ/NDQAAKl+pVNL09PS6fs2mpibZ\n9vLOWX35y1/Wv/3bv93w80899ZQeeOCB1Rpt2e677z498sgjeuCBB/TQQw/p0KFD+tKXvrTkv5PL\n5fTCCy/o05/+9KKfHx4e1q/8yq/opZdeUk9Pj/r6+vTss8/qnnvuueX5HMfRc889pwcffFC2bevL\nX/6yCoWCnn766Vv+s1aCmAQAoMpNT0/r//zfwwrXN6zL10vNzep/f/IOtbS0LOv1X/3qV/Unf/In\nkqTDhw/rkUce0SuvvLLw+UgksiZz3opnnnlGPp/vpq/7z//8T333u9+9YUx2dHTolVdeUXNz84pn\n+sUvfqGvf/3r+q3f+i3Ztq2vfvWrK/4z3SAmAQCoAeH6BkUamkyPsahIJLIQjNFoVJIUi8VMjvQ+\njY2Ny3qd4zhLft7j8aza93b91zIV3dwzCQAAyt4zzzyjz3/+83rooYd055136uWXX9ZHP/pR/fCH\nP1x4zeuvv66+vj4VCgVJ0tjYmL7whS/owIED+shHPqKnn35auVzuhl/jueee04c//GHdfvvt+t73\nvnfN5x566CF9+9vfliSNjo7qs5/9rA4dOqS77rpLjz32mFKplF5//XU99thjGh8fV19fn4aHh/XQ\nQw/piSee0Mc+9jF96EMf0vHjx9XX16cLFy4s/NmHDx/Wxz/+ce3fv19f/OIXlUgkFv1+pCu3BDz6\n6KMaHh7W7/3e70mS9uzZo9dff33hc+/68Y9/rN/4jd/Qvn379Gu/9mt68cUXr/l+/vqv/1p/8Ad/\noH379uljH/uYfvKTn9zy/y4SMQkAACrEj3/8Y3384x/XP/3TP+nQoUNLvtZxHP3RH/2RotGonn/+\neT399NP6n//5H/35n//5oq//6U9/qm9+85v60pe+pOeee05vvfWWxsfHF33tE088Ia/Xq+eff17/\n+I//qCNHjuhv//ZvdfDgQX3lK19RLBbTK6+8oo6ODknSCy+8oKeeekrf/e531dT0/rPD//zP/6yv\nfOUr+sEPfqDBwUE9+eSTNz0WHR0deuaZZyRJL7/8sg4ePHjN53/+85/rj//4j/WpT31K//7v/64H\nH3xQjz76qI4dO7bwmr/7u7/TJz7xCf3Hf/yHdu/erccff1zFYvGmX/t6XOYGAAAVobGxUb/7u7+7\nrNe+9tprGh4e1r/+67/K4/FIkr72ta/p93//9/Xoo4/K6702gX74wx/qE5/4hH79139dkvTNb35T\nH/7whxf9s0dGRtTX16euri75/X595zvfkWVZ8vv9ikQism37mkvZ9913n+644w5JVxbgXO8LX/jC\nwtd6/PHH9fDDD+trX/vakt+fx+NZuCWgpaXlfd/PD37wA/3qr/6qHn74YUnSli1bdPToUf3DP/yD\n/uqv/mphrncXNv3hH/6hPvWpT2l8fFydnZ1Lfu3rcWYSAABUhK6urmW/dmBgQLOzs7rjjjt08OBB\nHTx4UJ/73OeUz+d16dKlRV+/c+fOhd83Nzff8Ot97nOf04svvqi7775bX/ziF9Xf368tW7a4nnvv\n3r0Lv969e7eKxaIGBwdv8h0ubWBgQPv377/mnx08eFDnzp1b+P3GjRsXfl1fXy9J11xSXy7OTAIA\ngIoQCASW/Px7L9EWCgX19PS8795HSWpvb1/0379+QcuNVm/ff//9uueee/Tf//3fevnll/XYY4/p\nlVde0Z/92Z8t+nq/37/k3O/dQundGfx+vzKZzPteu9zYCwaD7/tnpVLpmmO02Pd3swVEi+HMJAAA\nqEg+n0+pVGrh90NDQwu/3rJli8bGxtTY2Kienh719PQoHo/rW9/61qLBtH37dh0/fnzh93Nzc9f8\nee/17W9/W2NjY/rt3/5tfec739GTTz6p//qv/5IkWZZ1y9/H6dOnF3597Ngx+Xw+bdy4cSH23vs9\nvvcy+VJfq7e3V0ePHr3mnx05cmTJM6huEZMAAKAi7d27Vy+88ILOnDmjN954Q88+++zC5z74wQ+q\nu7tbjz76qPr7+3XkyBE9/vjjsm170TOcn/nMZ/TSSy/pueee08DAgB5//HFls9lFv+65c+f0xBNP\n6OTJkzp37pxeeukl7dmzR5IUCoWUTCZ1/vz5ZZ9F/Mu//Eu9+uqrOnr0qJ588kk9+OCDCofD2r59\nu4LBoL73ve9paGhIzz77rE6ePLnw74VCIUnSyZMn3zfrww8/rB/96Ef6/ve/r8HBQX3/+9/Xj370\nI33mM59Z1ky3gpgEAKAGpOZmlZydXpf/pOZm1+V7euSRRxSNRvXAAw/oG9/4hh555JGFz3k8Hv3N\n3/yNPB6PPv3pT+vzn/+87rjjjhuulL7zzjv11FNP6e///u/1m7/5m2pra9OOHTsWfe3Xv/51tbW1\n6eGHH9YDDzygYrGob33rW5Kku+++W729vfrkJz+pU6dOLev7+OxnP7uw8ObAgQP60z/9U0lX7mP8\nxje+oRdffFH333+/Tpw4sbAdkCTt2LFDH/zgB/U7v/M779vWZ+/evXr66af1L//yL7r//vv1/PPP\n6y/+4i907733LmumW2E5N7k4Ho8nV/2LVpNYLMIxWgLH5+Y4Rkvj+Nwcx2hpsZj5p6eYVu6PU0Rl\nYwEOAABVzrbtZT/aELhV/L8MAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwC\nAADANWISAAAArhGTAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWIS\nAAAArhGTAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGT\nAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGTAAAAcI2Y\nBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGTAAAAcI2YBAAAgGvE\nJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4ZjmO45geAgAAAJXJe7MXxOPJ9ZijYsViEY7REjg+N8cx\nWhrH5+Y4RkuLxSKmRwCqGpe5AQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgE\nAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8Qk\nAAAAXCMmAQAA4JrX9AAAUMsKxZIkybIky7JkW5bhiQDg1hCTALCKcvmisvmicvnSlY+ForK5onKF\n0tWPRWXzJeXyReXyRRVLzvv+DNuyFuLyvR9ty1LA71E46FM46F34GKoPynEcWYQoAAOISQC4RY7j\nKJ0tKJnOK5nOLXycm88rXyit+M8vOY7kSFf/6xrz2YISyew1/yx8bkqZ+ZxC70ZmnU8tDUHFGuvk\n83I3E4C1RUwCwA2UHEep+fw10Tg3n1dyPq9iceXRuJqKJefqjDlJ0sDIjGzLUlNDQG1NIbU11Sla\nHzA8JYBqREwCgK6cbZxN5zWdzGh6Nqvpuazm5vMqLXIZulKUHEeTMxlNzmR0clAK+r3qaAlpc3uE\nsASwaohJADUply9qKpnV9GxGU8msEnPZVblEXc4yuYLOj87q/OismiIBbW5vUFcsLK+HS+EA3CMm\nAdSETK5w5Szd7JUzdbPpvByncs86rtR0MqvpZFynLkyrb1OjetojrCQH4AoxCaAqZfNFxafnNXE1\nIN+9lxDXyuQKOnp2QmdHZrSrp0ndsXrTIwGoMMQkgKoxM5fV2FRa49Pzmk5ma/rM461Kzed1uP+y\nRuIpHdjWqoDfY3okABWCmARQsQrFkuKJeY1PzWt8Oq35bMH0SBVvdDKlqWRGB7a1qqMlbHocABWA\nmARQUdKZvMam5jU+ldbEzPyim35jZbK5ol4/Oa7dm5u1Y2Oj6XEAlDliEkBZK5UcTSTmNT49r7Gp\nNPc+rqOTg1PK5Ara29vC03UA3BAxCaDslBxH8el5DcfnlMyOKjEzb3qkmnXu0qxy+ZJu74sRlAAW\nRUwCKAvO1Q22hydSGp1IKZsvSpLCYTbXNm04Pqdw0Ktdm5tNjwKgDBGTAIyaTmY1MjGnkXiKBTRl\n7PRQQuE6nza1RUyPAqDMEJMA1l0yndNIPKXh+Jzm5vOmx8EyHT07ocb6gBrCftOjACgjxCSAdTGf\nLSwEZGIua3ocuFAsOTryTlz37e/k/kkAC4hJAGsmmy/q0sSVgJyaZRPxajCdzGpwLKktHQ2mRwFQ\nJohJAKsunpjXhbGkLk2mVGIfyKrTf3Fam9rq5bFt06MAKAPEJIBVkc0XNTQ+p8GxWe6DrHLZXFFD\nl+e0uZ2zkwCISQAr9O5ZyNHJFE+jqSEDI7PqaYtw7yQAYhLAreMsJJLpnOIzGW1orDM9CgDDiEkA\nyzaRmNcgZyFx1dhkipgEQEwCWBpnIXEj41Pz0lbTUwAwjZgEsKjJmYzOj85yFhI3lMrkNZvOqSHE\nJuZALSMmASwoOY5GJ1IauDSrqdmM6XFQARLJLDEJ1DhiEoDyhZIujCd17tKs0hkuZWP5kmn+vgC1\njpgEalg6U9C50RldGEsqXyiZHgcVKJnOmR4BgGHEJFCDppNZDYzM6NJESiUecYgVYFEWAGISqBGO\n42hsKq2zIzOanOF+SKwOzmgDICaBKlcoljR0eU4DIzOcRcKqKxSJSaDWEZNAlcrkCjp/aVaDY0ll\n80XT46BKFUuOSo4jm8cqAjWLmASqTDqT15mhGQ1dTrI/JNZFqeTI9hCTQK0iJoEqkcrkdeZiQkOX\n51hUg3VjWZY8NiEJ1DJiEqhwc/N5nRlKaJiIhAFejyWLS9xATSMmgQqVTOd0ZiihkTjb+8Acr8c2\nPQIAw4hJoMLMpnM6czGhkYmUHCIShvm8xCRQ64hJoELMpnI6PZTQJSISZSQU5McIUOt4FwDK3Ewq\np9MXpzU6mSYiUXYaQn7TIwAwjJgEytTMXFanhxJEJMpahJgEah4xCZSZxFxW/RenNTaZNj0KcFOR\nkM/0CAAMIyaBMpHK5HVqcJqFNagYlmWpvo6YBGodMQkYlskWdGxgQoNjSZV4Yg0qSCjoZWsgAMQk\nYEqhWNLAyIxGE6NKzMybHge4ZY31AdMjACgDxCSwzkqOowtjSZ2+mFAmV1A4zA9kVKa2pjrTIwAo\nA8QksI5G4nM6dWFac/N506MAK2JbltqbQ6bHAFAGiElgHUwk5vX24JSmk1nTowCroiUalN/nMT0G\ngDJATAJraCaV08nBKY1Psc0PqktHS9j0CADKBDEJrIF0Jq9TFxIajs+xzQ+qjmVZ6mjhEjeAK4hJ\nYBXl8kWdGUro/OisimzzgyoVrferLsCPDwBX8G4ArIKS4+j8pVn1X5xWvlAyPQ6wpjpYeAPgPYhJ\nYIUuJ+Z1fGBSyXTO9CjAuuho5X5JAP8fMQm4lM7kdeL8lC5NpEyPAqyb+jqfGkJ+02MAKCPEJHCL\nCsWS3hme0dnhBPdFoub0tEdMjwCgzFgOS02BZbswNqsjp+NKZ9h0HLXH57X1qfu2sr8kgGvc9Mxk\nPJ5cjzkqViwW4RgtoVqOz2wqp+PnJhVPrP4ztMPhgFIpNjO/EY7Pza3XMdrWFdVMovL2TI3FOJsK\nrCUucwNLyOWL6r84rcHRpEqcxEcNsy1LvZ1R02MAKEPEJLAIx3E0OJZU/4VpZfNF0+MAxnW2hhUK\n8iMDwPvxzgBcZ3Imo+PnJpWY49Iq8K5t3ZyVBLA4YhK4aj5b0NvnpzQcnzM9ClBWYo11aqwPmB4D\nQJkiJlHzSo6jcyNXnl5TKPL0GuB6W7s4KwngxohJ1LTpZFZvnZ3QDJe0gUVFQn61NdWZHgNAGSMm\nUZPyhZJOXZjS+dGk2GoVuLFtXVFZlmV6DABljJhEzRmZSOn4wKQyuYLpUYCyFg761L2B53ADWBox\niZqRzhR0bGBCY1OVt+kyYMKeLc3y2LbpMQCUOWISVa/kODp3aVb9F1hgAyxXrLFOna2clQRwc8Qk\nqtrMXFZHzk4okWSBDbBctmXptt4W02MAqBDEJKpSsVRS/8WEBoZneAwicIs2tUcUDftNjwGgQhCT\nqDoTiXm9dXZCc/N506MAFcfntbWrp8n0GAAqCDGJqpEvFPX2+WldGGe7H8Ctvk1NCvg8pscAUEGI\nSVSF0cmUjp5lux9gJSIhv3o7GkyPAaDCEJOoaNl8UccGJjXC87SBFbttS7Nsmw3KAdwaYhIVa2wq\nrbfemeBsJLAK2ppDamsOmR4DQAUiJlFx8oWSTpyb1IXxpOlRgKpgW5Zu29JsegwAFYqYREWZSMzr\nzXcmlM6wUhtYLdu6o4qE2AoIgDvEJCpCsVTSycFpnbs0y0ptYBVF6wPauYmtgAC4R0yi7E0ns3rz\nTFzJdM70KEBV8diWDu2IsegGwIoQkyhbpZKj00MJvTOU4Ck2wBrY2dPEk24ArBgxibI0m8rpzTNx\nJeZ4pjawFloagtrWFTU9BoAqQEyirDiOo7MjM+q/MK1iibORwFrwemwd6ovJsri8DWDliEmUjVQm\nrzdPxzU5mzE9ClDVbuttVjjoMz0GgCpBTKIsnB+d1dvnp1QolkyPAlS1tuaQNrfzyEQAq4eYhFGZ\nXEFHzkxofDptehSg6gV8Hh3c3mp6DABVhpiEMePTab15Jq5srmh6FKAm7NvaoqCft30Aq4t3Fay7\nkuPo1OC0zo7MsAE5sE66Y/XqitWbHgNAFSImsa7SmbwOn45rikU2wLqpC3i1b2uL6TEAVCliEutm\nZCKlt96JK19gkQ2wXmzb0p07N8jv85geBUCVIiax5oqlko6fm9Lg6KzpUYCas2dLs5obgqbHAFDF\niEmsqZm5rF5+65JmUjxXG1hv3bF6be3kKTcA1hYxiTVzYSypgfGkZglJYN1FQn4dYBsgAOuAmMSq\nyxdKOjYwoaHLcwqHA6bHAWqOz2vrrl0b5PXYpkcBUAOISayqxFxWv+i/rNR83vQoQM06sD2mSMhv\negwANYKYxKoZGJnR24NTKpXYOxIwZfvGRnW1hk2PAaCGEJNYsVy+qDffiWtskkciAia1NYe0u6fJ\n9BgAagwxiRWZms3ocP9lpbMF06MANa2+zqc7+mKyLMv0KABqDDEJ1wbHZnVsYJLL2oBhPq+t27fF\n5POyMTmA9UdM4paVSo6ODUxqcIxNyAHTLMvSB/Z2KMAJSQCGsG8Ebsl8tqBXjo8SkkCZ2LmpUd0b\nIqbHAFDDODOJZZuYmdfh/rgyOe6PBMpBT3tEfZtYcAPALGISy3Lu0qxOnOf+SKBctDeHtH8bT7gB\nYB4xiSUVSyUdPTupi+NJ06MAuKopEtAdOzfIZuU2gDJATOKG0pmC3ugfVyKZNT0KgKvCdT7dvbud\nRyUCKBvEJBYVT8zrcP9lZfNF06MAuCrg9+gDe9oV8LMFEIDyQUzifc4Oz+jk4JRKDvdHAuXC67F1\n9+521df5TI8CANcgJrGgUCzprXcmNByfMz0KgPewLUt37tygpkjA9CgA8D7EJCRJqUxeb5wc10wq\nZ3oUANfZv61Vbc0h02MAwKKISSiemNcv+i8rx/2RQNnZualJPe1sSg6gfBGTNe7CWFJHBybYPxIo\nQz3tEe3sYVNyAOWNmKxRjuPo1IVpnRlKmB4FwCLYlBxApSAma1CxVNKbZyY0wkIboCyxKTmASkJM\n1phsrqjXT41rajZjehQAi2isD+gDe9iUHEDlICZrSDKd02tvjyuVyZseBcAiomG/7rmtXX4fm5ID\nqBzEZI2IJ+b1xqlx5Qsl06MAWERD2K979nYQkgAqDjFZAy6OJ/XWWVZsA+UqEvLr3ts6FCAkAVQg\nYrKKsWIbKH/1dT7du5fnbQOoXMRklSqWSjpyhkcjAuUsXOfTvXs7FPTzVgygcvEOVoWy+aJeP8mK\nbaCchYM+3Xtbh+oCvA0DqGy8i1WZZDqn106OKzXPim2gXIWCVy5th4K8BQOofLyTVZHJmYxePzXO\nM7aBMhYKeK+GpM/0KACwKojJKjE+ldYb/ZdVLLL1D1Cu6gJe3bO3Q2FCEkAVISarwHB8Tm+eibP1\nD1DGgn6v7t3bofo6QhJAdSEmK9z50VkdG5iU4xCSQLkK+D26Z287IQmgKhGTFezMUEInB6dMjwFg\nCaGAVx+4rV2RkN/0KACwJojJCnXi/KTODs+YHgPAEurrfLrntg5WbQOoarzDVZiS4+jo2QldGEua\nHgXAEhrrA/rAHp5sA6D6EZMVpFRydPj0ZV2aSJkeBcASWqJB3b27TT4vIQmg+lkOKzcqQr5Q0k/f\nGtHYJCEJlLPO1rA+eKBLXo9tehQAWBc3jcl4nMupS4nFImt+jHL5ol6r0McjhsMBpVJZ02OUNY7R\n0irp+HTH6nVoR0y2ba3r112P96FKFotFTI8AVDUuc5e5TK6gn58Y00wqZ3oUAEvY3NGg/VtbZFnr\nG5IAYBoxWcZSmbxePTHGc7aBMrdjY6N2b242PQYAGEFMlqnZdE6vHh9TJlcwPQqAJdy2pUXbuqOm\nxwAAY4jJMjSbzulnx0eVzRVNjwLgBizL0oFtrepp5348ALWNmCwzc/N5vXp8jJAEyphtW7q9b4O6\nWsOmRwEA44jJMjI3n9crx0a5tA2UMa/H1l27NmhDU8j0KABQFojJMpHK5PXqcUISKGcBv0f/a1eb\nmhuCpkcBgLJBTJaBdCavnx0fUzpLSALlKhLy6+49bQoHfaZHAYCyQkwals4UroRkhu1/gHIVa6zT\nXbs28HhEAFgEMWnQfLagn50YVYqQBMpWT1tE+7e1rvtTbQCgUhCThsxnC/rZ8VE2JAfKlGVZ2tXT\npB0bG02PAgBljZg0IJMr6NUTY5ojJIGy5LEtHdoRU1es3vQoAFD2iMl1ls0V9eqJMSXTPGsbKEes\n2AaAW0NMrqNsvqhXT4xqNkVIAuWIFdsAcOuIyXWSyxf18xNjmiEkgbLEim0AcIeYXAeFYkmvnRxX\nYi5rehQAi9jUFtEBVmwDgCvE5BorOY5+eTquqdmM6VEAXMeyLO3c1Ki+TU2mRwGAikVMrrHjA5Ma\nnUyZHgPAdTy2pYM7YupmxTYArAgxuYbODCV0fnTW9BgArhP0e3XXrg2s2AaAVUBMrpGL40mdHJwy\nPQaA6zQ3BHXnzg2qC/D2BwCrgXfTNXB5Oq23zk6YHgPAdTZ3NGhfbwsLbQBgFRGTqywxl9Uv+i+r\nVHJMjwLgKo9tad/WVvW0R0yPAgBVh5hcRelMXq+9Pa58oWR6FABX1QW8umtXm5oiAdOjAEBVIiZX\nSS5f1M/fHlcmVzA9CoCrWqJB3bWzTQE/G5EDwFohJlfBu5uS87xtoHz0dkZ1W2+zbIv7IwFgLRGT\nK+SwKTlQVjweWwe2tWrjBvaPBID1QEyu0C/7L7MpOVAmQkGf7tq1QY313B8JAOuFmFyBgZEZnRuf\nMz0GAEkbmup0e98GBXzcHwkA64mYdOlyYl5vn59SXchvehSg5m3vbtSuzU3cHwkABhCTLqQyef2y\n/7JKDntJAiZ5PbYObm9VF8/XBgBjiMlbVCiW9MbJcWXzRdOjADUtEvLrzp0b1BDm6gAAmERM3qIj\n70xoJsUWQIBJPW0R7d3aIq/HNj0KANQ8YvIWnBlKaCTOghvAFJ/X1v6trepm2x8AKBvE5DKNT6d1\n6sK06TGAmtUYCeiOvg2qr/OZHgUA8B7E5DKkM3m9eTouhwU3gBFbu6Las7lZts1qbQAoN8TkTRRL\nJb3Rf5kFN4ABAZ9H9x3sUoCGBICyxd3rN3H83JQSyazpMYCa0xIN6iMHu9S9IWJ6FADAEjgzuYSL\n40kNjs6aHgOoKZZlacfGRvVtamQTcgCoAMTkDcykcjo6MGl6DKCmBP1e3d4XU6yxzvQoAIBlIiYX\nUSiWdLj/sorFkulRgJrR1hTSoR0xBfw8WxsAKgkxuYiTg9NKptmYHFgPtm1pd0+ztnY1yOKyNgBU\nHGLyOpcT8zrPfZLAuggHfbq9L6bmhqDpUQAALhGT75EvFPXWGfaTBNbDpraI9vY2y+flsjYAVDJi\n8j2ODUwpnS2YHgOoagG/Rwe2taqjJWx6FADAKiAmr7o0kdLQ5aTpMYCq1tka1v6trSyyAYAqQkxK\nyuaKOjowYXoMoGr5vLb29rZoUxsbkANAtSEmJR05G1c2x+MSgbWwoalOB7fHVBfg7QYAqlHNv7tf\nGEtqbDJtegyg6ng8tm7b0qzN7RG2/AGAKlbTMZnO5HXiPE+5AVZbc0NQh3bEVF/nMz0KAGCN1WxM\nOo6jN89MKF/gKTfAarFtSzs3NWlbd5TnagNAjajZmDw3OquJmXnTYwBVI1of0KEdMUXDftOjAADW\nUU3GZDZXVP+FadNjAFXBtixt645q56Ym2TZnIwGg1tRkTJ66MM3lbWAV1Nf5dGgHj0MEgFpWczE5\nM5fVhXE2JwdWwrIs9XY2aFdPk7we2/Q4AACDai4mj5+b4tnbwAo0hP06sK2Vs5EAAEk1FpMjEykW\n3QAu2bbGsG00AAALJElEQVSlvo2N2t7dyL2RAIAFNROTxVJJJ89PmR4DqEjNDUEd2N6qhhArtQEA\n16qZmBwYmVUqkzc9BlBRfF5bu3qataWDp9gAABZXEzE5ny3ozFDC9BhARWlvDmnf1laFgjXxNgEA\ncKkmfkqcHJxWochWQMByBP1e7e1tVles3vQoAIAKUPUxOZ3Majg+Z3oMoOxZlqXN7RHt3twkn9dj\nehwAQIWo+pg8cW6SrYCAm4iG/drPdj8AABeqOiYnZzKanM2YHgMoWx6Prb6NjdrWHZXNAhsAgAtV\nHZNnR2ZMjwCUrbamkPZta1E46DM9CgCgglVtTM7N5zU2lTY9BlB26gJe7dnSrG4W2AAAVkHVxuTZ\nkRnulQTew2Nb2toV1Y6NjTxPGwCwaqoyJrO5ooYus4IbeFd7S0h7e7mkDQBYfVUZk+dHZ1VkX0lA\nkZBft/U2q60pZHoUAECVqrqYLBRLOj82a3oMwCif11bfxib1djbItlmlDQBYO1UXk0OX55TNFU2P\nARhhWZY2bqjX7s1NCvqr7v+8AQBlqKp+2jiOowG2A0KNaooEtLe3hY3HAQDrynKqaMnz8OWkXj4y\nYnoMYF0F/R7t3x5Tb1dUFhuPAwDW2U1jMh5PrtcsK/aL/ssaWefncIfDAaVS2XX9mpWE43Nzbo+R\nbVna0tmgnZsaq/pZ2rFYpKLeh0zgGC0tFouYHgGoalVzmbtYKunyNJuUozbEGuu0d2uLGkJ+06MA\nAGpc1cRkPJFRvsB2QKhu4Tqf9mxuVmdr2PQoAABIqqKYHJ1MmR4BWDNBv1d9mxrV0xZhqx8AQFmp\niph0HIfncKMqeT22tnVFta07yiMQAQBlqSpicmo2y96SqCq2bWlze0R9G5sU8Ffv4hoAQOWripgc\nneISN6qDZVnqag1rZ0+T6ut4jjYAoPxVR0xOcokblS/WWKfdm5vVFAmYHgUAgGWr+JicSeWUms+b\nHgNwrSkS0L4tTWprCpkeBQCAW1bxMcnekqhUoaBPu3qadHB3uyYm1nezfQAAVkvFx+RsirOSqCwB\nn0c7NjZqS0eDbNviEYgAgIpW8TGZTOdMjwAsi8dja1tng7Z1N8rnZZsfAEB1qOiYdBxHSe6XRJnz\nemxt7ohoe1cj2/wAAKpORcdkOltQscgjFFGevB5bWzobtK0zSkQCAKpWRcdkMs1ZSZQfn9fWlo4G\nbe2KKuAjIgEA1a3CY5L7JVE+fF5bvZ1Rbe1skJ+IBADUiIqOyflswfQIgHxeW1s7o9ra1SCfl4gE\nANSWio7JfMExPQJqmN/n0dbOBvV2EpEAgNpV0TFZYPENDAj4PNraFdWWjga2+AEA1DxiElimdyOy\nt7NBXg8RCQCAVOExWSpxmRtrry7g1dbOqDZ3RIhIAACuU9Ex6fHwGDqsncb6gLZ2RdXVGpZt83cN\nAIDFVHRMcpYIq82yLLU11WlbV1StjXWmxwEAoOxVdEx6bGISq8PjsbVxQ722djYoEvKbHgcAgIpR\n0TEZ5BF1WKGA36PejgZt7mjgaTUAALhQ0TFZX+czPQIqVCTk17auqLo3hDnDDQDAClR0TEZCxCRu\nTayxTtu6o9rQWCfLYlENAAArVeExyb1tuDnbttQdq9fWrqiiYf7OAACwmio6Jn1eW9H6gGbmsqZH\nQRkKBX3qaavXpraI6gIV/VcdAICyVfE/Ydub6ohJLLAtS+0tIfW0R7iUDQDAOqj8mGwJ6/RQwvQY\nMIyzkAAAmFHxP3Ub6/0KBbxKZwumR8E64ywkAADmVXxMWpal3q6oTpybND0K1glnIQEAKB9V8ZO4\npy2i0xenlS+UTI+CNcJZSAAAylNVxKTPa2tzR4Pe4d7JqhMO+rSJs5AAAJStqvnpvKM7qqHxOWVy\n3DtZ6YJ+rzpbw+qOhdXcEDQ9DgAAWELVxKTP69Ftvc063H/Z9Chwwee11dkSVteGerVGg7K5jA0A\nQEWompiUpO5YvUYn0xqJz5keBcvg8dhqbw6pOxbWhqY6npENAEAFqqqYlKSD21s1N59nI/MyZVuW\nYk116o7Vq705JJ+XgAQAoJJVXUx6Pbbu3t2mn7x1ifsny4RlWWpuCKg7Vq/O1rACPo/pkQAAwCqp\nupiUpLqAVx/Y06ZX3x5TNlc0PU5Nsi1LzQ1B7extVchrKRSsyr9qAADUvKr9CR+tD+hD+zr187fH\nlJrPmx6nJgR8Hm1oCqmtuU5tTXXyeT2KxSKKx5OmRwMAAGukamNSkurrfLpvX6deOzmm6ST3UK6F\naH1A7U11amsOqSkSYDNxAABqTFXHpCQF/B59aF+nTl2c1sDwjEqOY3qkiub12Io11l09+xhiI3EA\nAGpcTZSAbVvas7lZHc0h/fJMnMvetygc9Kmt+crl69ZokC18AADAgpqIyXc1NwT10UNdOn8pqXeG\nE8rmWZyzmFDQp5aGoFqjQbVEg6qv85keCQAAlKmaiklJ8ti2tnVH1dMe0dmRGZ27NKN8oWR6LKMi\nIb9aGq6EY0tDkJXXAABg2Wq2GnxeW7t6mrS9O6qReErnx2aVqIFFOj6vrcb6gJojATU3BNUUCcjP\nvo8AAMClmo3Jd3k9tnraI+ppjygxl9VIPKXx6bRmUznTo61YwOdRfcinSMivxnq/miJBNYR8rLgG\nAACrpuZj8r0a6wNqrA9oz5ZmzWcLujw9r/HptOKJ+bK+FB70exW5Go1XPl75NU+aAQAAa42YvIG6\ngHfhjGXJcZRIZjU3n1dqPq9UpqBU5srH9WDblgI+j/w+j+r8V8821v3/cPR5iUYAAGAGMbkM7z4a\nsLkh+L7PRRtDGhyaViqTVzpTUDpbUKnkyHEcOY7kSHIcR6V3f3/1Y+nqJ71e+2ooXv3ofc+vfR4F\nfB75vGzFAwAAyhMxuUJ+n0dNkYCaIgHTowAAAKw7TnkBAADANWISAAAArhGTAAAAcI2YBAAAgGvE\nJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGTAAAAcI2YBAAAgGvEJAAAAFwj\nJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArlmO4zimhwAAAEBl8t7sBfF4cj3mqFixWIRj\ntASOz81xjJbG8bk5jtHSYrGI6RGAqsZlbgAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8Qk\nAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMm\nAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4Box\nCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeI\nSQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhG\nTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1\nYhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACu\nEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABw\njZgEAACAa8QkAAAAXCMmAQAA4JrlOI5jeggAAABUJs5MAgAAwDViEgAAAK4RkwAAAHCNmAQAAIBr\nxCQAAABcIyYBAADg2v8DSm3DUEfeZiIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1115399e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax  = plt.subplots(figsize=(8, 8))\n",
    "ax.set_aspect('equal');\n",
    "\n",
    "\n",
    "var, U = np.linalg.eig(true_cov)\n",
    "angle = 180. / np.pi * np.arccos(np.abs(U[0, 0]))\n",
    "\n",
    "e = Ellipse(np.zeros(2), 2 * np.sqrt(5.991 * var[0]), 2 * np.sqrt(5.991 * var[1]), angle=angle)\n",
    "e.set_alpha(0.5)\n",
    "e.set_facecolor(blue)\n",
    "e.set_zorder(10);\n",
    "ax.add_artist(e);\n",
    "\n",
    "ax.set_xlim(-3, 3);\n",
    "ax.set_xticklabels([]);\n",
    "\n",
    "ax.set_ylim(-3, 3);\n",
    "ax.set_yticklabels([]);\n",
    "\n",
    "rect = plt.Rectangle((0, 0), 1, 1, fc=blue, alpha=0.5)\n",
    "ax.legend([rect],\n",
    "          ['True distribution'],\n",
    "          bbox_to_anchor=(1.5, 1.));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApMAAAHBCAYAAAA1uDCcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3WlsXOd97/HfObNyhsPhNuIqUaIWarFWL9exEydIGwRF\njKQ12jpo6l6jDYI0RQMHdYE4MYLAceC+cJq2TtOmLepcFEHdBnbRi7YGnAJpHMexHcWyFkuULEqU\nSIqkhsuQwxnOfu4LybyWTFHU4fLM8v0ACSVzLP55og6/Pec8z7Ecx3EEAAAAuGCbHgAAAACVi5gE\nAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALi2ZEyyaxAAAACW4l3qk5ZlKR5PrtcsFSkWi3CM\nlsDxuTmO0dI4PjfHMVpaLBYxPQJQ1bjMDQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgE\nAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMC1JZ/NDQAAKl+pVNL09PS6fs2mpibZ\n9vLOWX35y1/Wv/3bv93w80899ZQeeOCB1Rpt2e677z498sgjeuCBB/TQQw/p0KFD+tKXvrTkv5PL\n5fTCCy/o05/+9KKfHx4e1q/8yq/opZdeUk9Pj/r6+vTss8/qnnvuueX5HMfRc889pwcffFC2bevL\nX/6yCoWCnn766Vv+s1aCmAQAoMpNT0/r//zfwwrXN6zL10vNzep/f/IOtbS0LOv1X/3qV/Unf/In\nkqTDhw/rkUce0SuvvLLw+UgksiZz3opnnnlGPp/vpq/7z//8T333u9+9YUx2dHTolVdeUXNz84pn\n+sUvfqGvf/3r+q3f+i3Ztq2vfvWrK/4z3SAmAQCoAeH6BkUamkyPsahIJLIQjNFoVJIUi8VMjvQ+\njY2Ny3qd4zhLft7j8aza93b91zIV3dwzCQAAyt4zzzyjz3/+83rooYd055136uWXX9ZHP/pR/fCH\nP1x4zeuvv66+vj4VCgVJ0tjYmL7whS/owIED+shHPqKnn35auVzuhl/jueee04c//GHdfvvt+t73\nvnfN5x566CF9+9vfliSNjo7qs5/9rA4dOqS77rpLjz32mFKplF5//XU99thjGh8fV19fn4aHh/XQ\nQw/piSee0Mc+9jF96EMf0vHjx9XX16cLFy4s/NmHDx/Wxz/+ce3fv19f/OIXlUgkFv1+pCu3BDz6\n6KMaHh7W7/3e70mS9uzZo9dff33hc+/68Y9/rN/4jd/Qvn379Gu/9mt68cUXr/l+/vqv/1p/8Ad/\noH379uljH/uYfvKTn9zy/y4SMQkAACrEj3/8Y3384x/XP/3TP+nQoUNLvtZxHP3RH/2RotGonn/+\neT399NP6n//5H/35n//5oq//6U9/qm9+85v60pe+pOeee05vvfWWxsfHF33tE088Ia/Xq+eff17/\n+I//qCNHjuhv//ZvdfDgQX3lK19RLBbTK6+8oo6ODknSCy+8oKeeekrf/e531dT0/rPD//zP/6yv\nfOUr+sEPfqDBwUE9+eSTNz0WHR0deuaZZyRJL7/8sg4ePHjN53/+85/rj//4j/WpT31K//7v/64H\nH3xQjz76qI4dO7bwmr/7u7/TJz7xCf3Hf/yHdu/erccff1zFYvGmX/t6XOYGAAAVobGxUb/7u7+7\nrNe+9tprGh4e1r/+67/K4/FIkr72ta/p93//9/Xoo4/K6702gX74wx/qE5/4hH79139dkvTNb35T\nH/7whxf9s0dGRtTX16euri75/X595zvfkWVZ8vv9ikQism37mkvZ9913n+644w5JVxbgXO8LX/jC\nwtd6/PHH9fDDD+trX/vakt+fx+NZuCWgpaXlfd/PD37wA/3qr/6qHn74YUnSli1bdPToUf3DP/yD\n/uqv/mphrncXNv3hH/6hPvWpT2l8fFydnZ1Lfu3rcWYSAABUhK6urmW/dmBgQLOzs7rjjjt08OBB\nHTx4UJ/73OeUz+d16dKlRV+/c+fOhd83Nzff8Ot97nOf04svvqi7775bX/ziF9Xf368tW7a4nnvv\n3r0Lv969e7eKxaIGBwdv8h0ubWBgQPv377/mnx08eFDnzp1b+P3GjRsXfl1fXy9J11xSXy7OTAIA\ngIoQCASW/Px7L9EWCgX19PS8795HSWpvb1/0379+QcuNVm/ff//9uueee/Tf//3fevnll/XYY4/p\nlVde0Z/92Z8t+nq/37/k3O/dQundGfx+vzKZzPteu9zYCwaD7/tnpVLpmmO02Pd3swVEi+HMJAAA\nqEg+n0+pVGrh90NDQwu/3rJli8bGxtTY2Kienh719PQoHo/rW9/61qLBtH37dh0/fnzh93Nzc9f8\nee/17W9/W2NjY/rt3/5tfec739GTTz6p//qv/5IkWZZ1y9/H6dOnF3597Ngx+Xw+bdy4cSH23vs9\nvvcy+VJfq7e3V0ePHr3mnx05cmTJM6huEZMAAKAi7d27Vy+88ILOnDmjN954Q88+++zC5z74wQ+q\nu7tbjz76qPr7+3XkyBE9/vjjsm170TOcn/nMZ/TSSy/pueee08DAgB5//HFls9lFv+65c+f0xBNP\n6OTJkzp37pxeeukl7dmzR5IUCoWUTCZ1/vz5ZZ9F/Mu//Eu9+uqrOnr0qJ588kk9+OCDCofD2r59\nu4LBoL73ve9paGhIzz77rE6ePLnw74VCIUnSyZMn3zfrww8/rB/96Ef6/ve/r8HBQX3/+9/Xj370\nI33mM59Z1ky3gpgEAKAGpOZmlZydXpf/pOZm1+V7euSRRxSNRvXAAw/oG9/4hh555JGFz3k8Hv3N\n3/yNPB6PPv3pT+vzn/+87rjjjhuulL7zzjv11FNP6e///u/1m7/5m2pra9OOHTsWfe3Xv/51tbW1\n6eGHH9YDDzygYrGob33rW5Kku+++W729vfrkJz+pU6dOLev7+OxnP7uw8ObAgQP60z/9U0lX7mP8\nxje+oRdffFH333+/Tpw4sbAdkCTt2LFDH/zgB/U7v/M779vWZ+/evXr66af1L//yL7r//vv1/PPP\n6y/+4i907733LmumW2E5N7k4Ho8nV/2LVpNYLMIxWgLH5+Y4Rkvj+Nwcx2hpsZj5p6eYVu6PU0Rl\nYwEOAABVzrbtZT/aELhV/L8MAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwC\nAADANWISAAAArhGTAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWIS\nAAAArhGTAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGT\nAAAAcI2YBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGTAAAAcI2Y\nBAAAgGvEJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGTAAAAcI2YBAAAgGvE\nJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4ZjmO45geAgAAAJXJe7MXxOPJ9ZijYsViEY7REjg+N8cx\nWhrH5+Y4RkuLxSKmRwCqGpe5AQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgE\nAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8Qk\nAAAAXCMmAQAA4JrX9AAAUMsKxZIkybIky7JkW5bhiQDg1hCTALCKcvmisvmicvnSlY+ForK5onKF\n0tWPRWXzJeXyReXyRRVLzvv+DNuyFuLyvR9ty1LA71E46FM46F34GKoPynEcWYQoAAOISQC4RY7j\nKJ0tKJnOK5nOLXycm88rXyit+M8vOY7kSFf/6xrz2YISyew1/yx8bkqZ+ZxC70ZmnU8tDUHFGuvk\n83I3E4C1RUwCwA2UHEep+fw10Tg3n1dyPq9iceXRuJqKJefqjDlJ0sDIjGzLUlNDQG1NIbU11Sla\nHzA8JYBqREwCgK6cbZxN5zWdzGh6Nqvpuazm5vMqLXIZulKUHEeTMxlNzmR0clAK+r3qaAlpc3uE\nsASwaohJADUply9qKpnV9GxGU8msEnPZVblEXc4yuYLOj87q/OismiIBbW5vUFcsLK+HS+EA3CMm\nAdSETK5w5Szd7JUzdbPpvByncs86rtR0MqvpZFynLkyrb1OjetojrCQH4AoxCaAqZfNFxafnNXE1\nIN+9lxDXyuQKOnp2QmdHZrSrp0ndsXrTIwGoMMQkgKoxM5fV2FRa49Pzmk5ma/rM461Kzed1uP+y\nRuIpHdjWqoDfY3okABWCmARQsQrFkuKJeY1PzWt8Oq35bMH0SBVvdDKlqWRGB7a1qqMlbHocABWA\nmARQUdKZvMam5jU+ldbEzPyim35jZbK5ol4/Oa7dm5u1Y2Oj6XEAlDliEkBZK5UcTSTmNT49r7Gp\nNPc+rqOTg1PK5Ara29vC03UA3BAxCaDslBxH8el5DcfnlMyOKjEzb3qkmnXu0qxy+ZJu74sRlAAW\nRUwCKAvO1Q22hydSGp1IKZsvSpLCYTbXNm04Pqdw0Ktdm5tNjwKgDBGTAIyaTmY1MjGnkXiKBTRl\n7PRQQuE6nza1RUyPAqDMEJMA1l0yndNIPKXh+Jzm5vOmx8EyHT07ocb6gBrCftOjACgjxCSAdTGf\nLSwEZGIua3ocuFAsOTryTlz37e/k/kkAC4hJAGsmmy/q0sSVgJyaZRPxajCdzGpwLKktHQ2mRwFQ\nJohJAKsunpjXhbGkLk2mVGIfyKrTf3Fam9rq5bFt06MAKAPEJIBVkc0XNTQ+p8GxWe6DrHLZXFFD\nl+e0uZ2zkwCISQAr9O5ZyNHJFE+jqSEDI7PqaYtw7yQAYhLAreMsJJLpnOIzGW1orDM9CgDDiEkA\nyzaRmNcgZyFx1dhkipgEQEwCWBpnIXEj41Pz0lbTUwAwjZgEsKjJmYzOj85yFhI3lMrkNZvOqSHE\nJuZALSMmASwoOY5GJ1IauDSrqdmM6XFQARLJLDEJ1DhiEoDyhZIujCd17tKs0hkuZWP5kmn+vgC1\njpgEalg6U9C50RldGEsqXyiZHgcVKJnOmR4BgGHEJFCDppNZDYzM6NJESiUecYgVYFEWAGISqBGO\n42hsKq2zIzOanOF+SKwOzmgDICaBKlcoljR0eU4DIzOcRcKqKxSJSaDWEZNAlcrkCjp/aVaDY0ll\n80XT46BKFUuOSo4jm8cqAjWLmASqTDqT15mhGQ1dTrI/JNZFqeTI9hCTQK0iJoEqkcrkdeZiQkOX\n51hUg3VjWZY8NiEJ1DJiEqhwc/N5nRlKaJiIhAFejyWLS9xATSMmgQqVTOd0ZiihkTjb+8Acr8c2\nPQIAw4hJoMLMpnM6czGhkYmUHCIShvm8xCRQ64hJoELMpnI6PZTQJSISZSQU5McIUOt4FwDK3Ewq\np9MXpzU6mSYiUXYaQn7TIwAwjJgEytTMXFanhxJEJMpahJgEah4xCZSZxFxW/RenNTaZNj0KcFOR\nkM/0CAAMIyaBMpHK5HVqcJqFNagYlmWpvo6YBGodMQkYlskWdGxgQoNjSZV4Yg0qSCjoZWsgAMQk\nYEqhWNLAyIxGE6NKzMybHge4ZY31AdMjACgDxCSwzkqOowtjSZ2+mFAmV1A4zA9kVKa2pjrTIwAo\nA8QksI5G4nM6dWFac/N506MAK2JbltqbQ6bHAFAGiElgHUwk5vX24JSmk1nTowCroiUalN/nMT0G\ngDJATAJraCaV08nBKY1Psc0PqktHS9j0CADKBDEJrIF0Jq9TFxIajs+xzQ+qjmVZ6mjhEjeAK4hJ\nYBXl8kWdGUro/OisimzzgyoVrferLsCPDwBX8G4ArIKS4+j8pVn1X5xWvlAyPQ6wpjpYeAPgPYhJ\nYIUuJ+Z1fGBSyXTO9CjAuuho5X5JAP8fMQm4lM7kdeL8lC5NpEyPAqyb+jqfGkJ+02MAKCPEJHCL\nCsWS3hme0dnhBPdFoub0tEdMjwCgzFgOS02BZbswNqsjp+NKZ9h0HLXH57X1qfu2sr8kgGvc9Mxk\nPJ5cjzkqViwW4RgtoVqOz2wqp+PnJhVPrP4ztMPhgFIpNjO/EY7Pza3XMdrWFdVMovL2TI3FOJsK\nrCUucwNLyOWL6r84rcHRpEqcxEcNsy1LvZ1R02MAKEPEJLAIx3E0OJZU/4VpZfNF0+MAxnW2hhUK\n8iMDwPvxzgBcZ3Imo+PnJpWY49Iq8K5t3ZyVBLA4YhK4aj5b0NvnpzQcnzM9ClBWYo11aqwPmB4D\nQJkiJlHzSo6jcyNXnl5TKPL0GuB6W7s4KwngxohJ1LTpZFZvnZ3QDJe0gUVFQn61NdWZHgNAGSMm\nUZPyhZJOXZjS+dGk2GoVuLFtXVFZlmV6DABljJhEzRmZSOn4wKQyuYLpUYCyFg761L2B53ADWBox\niZqRzhR0bGBCY1OVt+kyYMKeLc3y2LbpMQCUOWISVa/kODp3aVb9F1hgAyxXrLFOna2clQRwc8Qk\nqtrMXFZHzk4okWSBDbBctmXptt4W02MAqBDEJKpSsVRS/8WEBoZneAwicIs2tUcUDftNjwGgQhCT\nqDoTiXm9dXZCc/N506MAFcfntbWrp8n0GAAqCDGJqpEvFPX2+WldGGe7H8Ctvk1NCvg8pscAUEGI\nSVSF0cmUjp5lux9gJSIhv3o7GkyPAaDCEJOoaNl8UccGJjXC87SBFbttS7Nsmw3KAdwaYhIVa2wq\nrbfemeBsJLAK2ppDamsOmR4DQAUiJlFx8oWSTpyb1IXxpOlRgKpgW5Zu29JsegwAFYqYREWZSMzr\nzXcmlM6wUhtYLdu6o4qE2AoIgDvEJCpCsVTSycFpnbs0y0ptYBVF6wPauYmtgAC4R0yi7E0ns3rz\nTFzJdM70KEBV8diWDu2IsegGwIoQkyhbpZKj00MJvTOU4Ck2wBrY2dPEk24ArBgxibI0m8rpzTNx\nJeZ4pjawFloagtrWFTU9BoAqQEyirDiOo7MjM+q/MK1iibORwFrwemwd6ovJsri8DWDliEmUjVQm\nrzdPxzU5mzE9ClDVbuttVjjoMz0GgCpBTKIsnB+d1dvnp1QolkyPAlS1tuaQNrfzyEQAq4eYhFGZ\nXEFHzkxofDptehSg6gV8Hh3c3mp6DABVhpiEMePTab15Jq5srmh6FKAm7NvaoqCft30Aq4t3Fay7\nkuPo1OC0zo7MsAE5sE66Y/XqitWbHgNAFSImsa7SmbwOn45rikU2wLqpC3i1b2uL6TEAVCliEutm\nZCKlt96JK19gkQ2wXmzb0p07N8jv85geBUCVIiax5oqlko6fm9Lg6KzpUYCas2dLs5obgqbHAFDF\niEmsqZm5rF5+65JmUjxXG1hv3bF6be3kKTcA1hYxiTVzYSypgfGkZglJYN1FQn4dYBsgAOuAmMSq\nyxdKOjYwoaHLcwqHA6bHAWqOz2vrrl0b5PXYpkcBUAOISayqxFxWv+i/rNR83vQoQM06sD2mSMhv\negwANYKYxKoZGJnR24NTKpXYOxIwZfvGRnW1hk2PAaCGEJNYsVy+qDffiWtskkciAia1NYe0u6fJ\n9BgAagwxiRWZms3ocP9lpbMF06MANa2+zqc7+mKyLMv0KABqDDEJ1wbHZnVsYJLL2oBhPq+t27fF\n5POyMTmA9UdM4paVSo6ODUxqcIxNyAHTLMvSB/Z2KMAJSQCGsG8Ebsl8tqBXjo8SkkCZ2LmpUd0b\nIqbHAFDDODOJZZuYmdfh/rgyOe6PBMpBT3tEfZtYcAPALGISy3Lu0qxOnOf+SKBctDeHtH8bT7gB\nYB4xiSUVSyUdPTupi+NJ06MAuKopEtAdOzfIZuU2gDJATOKG0pmC3ugfVyKZNT0KgKvCdT7dvbud\nRyUCKBvEJBYVT8zrcP9lZfNF06MAuCrg9+gDe9oV8LMFEIDyQUzifc4Oz+jk4JRKDvdHAuXC67F1\n9+521df5TI8CANcgJrGgUCzprXcmNByfMz0KgPewLUt37tygpkjA9CgA8D7EJCRJqUxeb5wc10wq\nZ3oUANfZv61Vbc0h02MAwKKISSiemNcv+i8rx/2RQNnZualJPe1sSg6gfBGTNe7CWFJHBybYPxIo\nQz3tEe3sYVNyAOWNmKxRjuPo1IVpnRlKmB4FwCLYlBxApSAma1CxVNKbZyY0wkIboCyxKTmASkJM\n1phsrqjXT41rajZjehQAi2isD+gDe9iUHEDlICZrSDKd02tvjyuVyZseBcAiomG/7rmtXX4fm5ID\nqBzEZI2IJ+b1xqlx5Qsl06MAWERD2K979nYQkgAqDjFZAy6OJ/XWWVZsA+UqEvLr3ts6FCAkAVQg\nYrKKsWIbKH/1dT7du5fnbQOoXMRklSqWSjpyhkcjAuUsXOfTvXs7FPTzVgygcvEOVoWy+aJeP8mK\nbaCchYM+3Xtbh+oCvA0DqGy8i1WZZDqn106OKzXPim2gXIWCVy5th4K8BQOofLyTVZHJmYxePzXO\nM7aBMhYKeK+GpM/0KACwKojJKjE+ldYb/ZdVLLL1D1Cu6gJe3bO3Q2FCEkAVISarwHB8Tm+eibP1\nD1DGgn6v7t3bofo6QhJAdSEmK9z50VkdG5iU4xCSQLkK+D26Z287IQmgKhGTFezMUEInB6dMjwFg\nCaGAVx+4rV2RkN/0KACwJojJCnXi/KTODs+YHgPAEurrfLrntg5WbQOoarzDVZiS4+jo2QldGEua\nHgXAEhrrA/rAHp5sA6D6EZMVpFRydPj0ZV2aSJkeBcASWqJB3b27TT4vIQmg+lkOKzcqQr5Q0k/f\nGtHYJCEJlLPO1rA+eKBLXo9tehQAWBc3jcl4nMupS4nFImt+jHL5ol6r0McjhsMBpVJZ02OUNY7R\n0irp+HTH6nVoR0y2ba3r112P96FKFotFTI8AVDUuc5e5TK6gn58Y00wqZ3oUAEvY3NGg/VtbZFnr\nG5IAYBoxWcZSmbxePTHGc7aBMrdjY6N2b242PQYAGEFMlqnZdE6vHh9TJlcwPQqAJdy2pUXbuqOm\nxwAAY4jJMjSbzulnx0eVzRVNjwLgBizL0oFtrepp5348ALWNmCwzc/N5vXp8jJAEyphtW7q9b4O6\nWsOmRwEA44jJMjI3n9crx0a5tA2UMa/H1l27NmhDU8j0KABQFojJMpHK5PXqcUISKGcBv0f/a1eb\nmhuCpkcBgLJBTJaBdCavnx0fUzpLSALlKhLy6+49bQoHfaZHAYCyQkwals4UroRkhu1/gHIVa6zT\nXbs28HhEAFgEMWnQfLagn50YVYqQBMpWT1tE+7e1rvtTbQCgUhCThsxnC/rZ8VE2JAfKlGVZ2tXT\npB0bG02PAgBljZg0IJMr6NUTY5ojJIGy5LEtHdoRU1es3vQoAFD2iMl1ls0V9eqJMSXTPGsbKEes\n2AaAW0NMrqNsvqhXT4xqNkVIAuWIFdsAcOuIyXWSyxf18xNjmiEkgbLEim0AcIeYXAeFYkmvnRxX\nYi5rehQAi9jUFtEBVmwDgCvE5BorOY5+eTquqdmM6VEAXMeyLO3c1Ki+TU2mRwGAikVMrrHjA5Ma\nnUyZHgPAdTy2pYM7YupmxTYArAgxuYbODCV0fnTW9BgArhP0e3XXrg2s2AaAVUBMrpGL40mdHJwy\nPQaA6zQ3BHXnzg2qC/D2BwCrgXfTNXB5Oq23zk6YHgPAdTZ3NGhfbwsLbQBgFRGTqywxl9Uv+i+r\nVHJMjwLgKo9tad/WVvW0R0yPAgBVh5hcRelMXq+9Pa58oWR6FABX1QW8umtXm5oiAdOjAEBVIiZX\nSS5f1M/fHlcmVzA9CoCrWqJB3bWzTQE/G5EDwFohJlfBu5uS87xtoHz0dkZ1W2+zbIv7IwFgLRGT\nK+SwKTlQVjweWwe2tWrjBvaPBID1QEyu0C/7L7MpOVAmQkGf7tq1QY313B8JAOuFmFyBgZEZnRuf\nMz0GAEkbmup0e98GBXzcHwkA64mYdOlyYl5vn59SXchvehSg5m3vbtSuzU3cHwkABhCTLqQyef2y\n/7JKDntJAiZ5PbYObm9VF8/XBgBjiMlbVCiW9MbJcWXzRdOjADUtEvLrzp0b1BDm6gAAmERM3qIj\n70xoJsUWQIBJPW0R7d3aIq/HNj0KANQ8YvIWnBlKaCTOghvAFJ/X1v6trepm2x8AKBvE5DKNT6d1\n6sK06TGAmtUYCeiOvg2qr/OZHgUA8B7E5DKkM3m9eTouhwU3gBFbu6Las7lZts1qbQAoN8TkTRRL\nJb3Rf5kFN4ABAZ9H9x3sUoCGBICyxd3rN3H83JQSyazpMYCa0xIN6iMHu9S9IWJ6FADAEjgzuYSL\n40kNjs6aHgOoKZZlacfGRvVtamQTcgCoAMTkDcykcjo6MGl6DKCmBP1e3d4XU6yxzvQoAIBlIiYX\nUSiWdLj/sorFkulRgJrR1hTSoR0xBfw8WxsAKgkxuYiTg9NKptmYHFgPtm1pd0+ztnY1yOKyNgBU\nHGLyOpcT8zrPfZLAuggHfbq9L6bmhqDpUQAALhGT75EvFPXWGfaTBNbDpraI9vY2y+flsjYAVDJi\n8j2ODUwpnS2YHgOoagG/Rwe2taqjJWx6FADAKiAmr7o0kdLQ5aTpMYCq1tka1v6trSyyAYAqQkxK\nyuaKOjowYXoMoGr5vLb29rZoUxsbkANAtSEmJR05G1c2x+MSgbWwoalOB7fHVBfg7QYAqlHNv7tf\nGEtqbDJtegyg6ng8tm7b0qzN7RG2/AGAKlbTMZnO5HXiPE+5AVZbc0NQh3bEVF/nMz0KAGCN1WxM\nOo6jN89MKF/gKTfAarFtSzs3NWlbd5TnagNAjajZmDw3OquJmXnTYwBVI1of0KEdMUXDftOjAADW\nUU3GZDZXVP+FadNjAFXBtixt645q56Ym2TZnIwGg1tRkTJ66MM3lbWAV1Nf5dGgHj0MEgFpWczE5\nM5fVhXE2JwdWwrIs9XY2aFdPk7we2/Q4AACDai4mj5+b4tnbwAo0hP06sK2Vs5EAAEk1FpMjEykW\n3QAu2bbGsG00AAALJElEQVSlvo2N2t7dyL2RAIAFNROTxVJJJ89PmR4DqEjNDUEd2N6qhhArtQEA\n16qZmBwYmVUqkzc9BlBRfF5bu3qataWDp9gAABZXEzE5ny3ozFDC9BhARWlvDmnf1laFgjXxNgEA\ncKkmfkqcHJxWochWQMByBP1e7e1tVles3vQoAIAKUPUxOZ3Majg+Z3oMoOxZlqXN7RHt3twkn9dj\nehwAQIWo+pg8cW6SrYCAm4iG/drPdj8AABeqOiYnZzKanM2YHgMoWx6Prb6NjdrWHZXNAhsAgAtV\nHZNnR2ZMjwCUrbamkPZta1E46DM9CgCgglVtTM7N5zU2lTY9BlB26gJe7dnSrG4W2AAAVkHVxuTZ\nkRnulQTew2Nb2toV1Y6NjTxPGwCwaqoyJrO5ooYus4IbeFd7S0h7e7mkDQBYfVUZk+dHZ1VkX0lA\nkZBft/U2q60pZHoUAECVqrqYLBRLOj82a3oMwCif11bfxib1djbItlmlDQBYO1UXk0OX55TNFU2P\nARhhWZY2bqjX7s1NCvqr7v+8AQBlqKp+2jiOowG2A0KNaooEtLe3hY3HAQDrynKqaMnz8OWkXj4y\nYnoMYF0F/R7t3x5Tb1dUFhuPAwDW2U1jMh5PrtcsK/aL/ssaWefncIfDAaVS2XX9mpWE43Nzbo+R\nbVna0tmgnZsaq/pZ2rFYpKLeh0zgGC0tFouYHgGoalVzmbtYKunyNJuUozbEGuu0d2uLGkJ+06MA\nAGpc1cRkPJFRvsB2QKhu4Tqf9mxuVmdr2PQoAABIqqKYHJ1MmR4BWDNBv1d9mxrV0xZhqx8AQFmp\niph0HIfncKMqeT22tnVFta07yiMQAQBlqSpicmo2y96SqCq2bWlze0R9G5sU8Ffv4hoAQOWripgc\nneISN6qDZVnqag1rZ0+T6ut4jjYAoPxVR0xOcokblS/WWKfdm5vVFAmYHgUAgGWr+JicSeWUms+b\nHgNwrSkS0L4tTWprCpkeBQCAW1bxMcnekqhUoaBPu3qadHB3uyYm1nezfQAAVkvFx+RsirOSqCwB\nn0c7NjZqS0eDbNviEYgAgIpW8TGZTOdMjwAsi8dja1tng7Z1N8rnZZsfAEB1qOiYdBxHSe6XRJnz\nemxt7ohoe1cj2/wAAKpORcdkOltQscgjFFGevB5bWzobtK0zSkQCAKpWRcdkMs1ZSZQfn9fWlo4G\nbe2KKuAjIgEA1a3CY5L7JVE+fF5bvZ1Rbe1skJ+IBADUiIqOyflswfQIgHxeW1s7o9ra1SCfl4gE\nANSWio7JfMExPQJqmN/n0dbOBvV2EpEAgNpV0TFZYPENDAj4PNraFdWWjga2+AEA1DxiElimdyOy\nt7NBXg8RCQCAVOExWSpxmRtrry7g1dbOqDZ3RIhIAACuU9Ex6fHwGDqsncb6gLZ2RdXVGpZt83cN\nAIDFVHRMcpYIq82yLLU11WlbV1StjXWmxwEAoOxVdEx6bGISq8PjsbVxQ722djYoEvKbHgcAgIpR\n0TEZ5BF1WKGA36PejgZt7mjgaTUAALhQ0TFZX+czPQIqVCTk17auqLo3hDnDDQDAClR0TEZCxCRu\nTayxTtu6o9rQWCfLYlENAAArVeExyb1tuDnbttQdq9fWrqiiYf7OAACwmio6Jn1eW9H6gGbmsqZH\nQRkKBX3qaavXpraI6gIV/VcdAICyVfE/Ydub6ohJLLAtS+0tIfW0R7iUDQDAOqj8mGwJ6/RQwvQY\nMIyzkAAAmFHxP3Ub6/0KBbxKZwumR8E64ywkAADmVXxMWpal3q6oTpybND0K1glnIQEAKB9V8ZO4\npy2i0xenlS+UTI+CNcJZSAAAylNVxKTPa2tzR4Pe4d7JqhMO+rSJs5AAAJStqvnpvKM7qqHxOWVy\n3DtZ6YJ+rzpbw+qOhdXcEDQ9DgAAWELVxKTP69Ftvc063H/Z9Chwwee11dkSVteGerVGg7K5jA0A\nQEWompiUpO5YvUYn0xqJz5keBcvg8dhqbw6pOxbWhqY6npENAEAFqqqYlKSD21s1N59nI/MyZVuW\nYk116o7Vq705JJ+XgAQAoJJVXUx6Pbbu3t2mn7x1ifsny4RlWWpuCKg7Vq/O1rACPo/pkQAAwCqp\nupiUpLqAVx/Y06ZX3x5TNlc0PU5Nsi1LzQ1B7extVchrKRSsyr9qAADUvKr9CR+tD+hD+zr187fH\nlJrPmx6nJgR8Hm1oCqmtuU5tTXXyeT2KxSKKx5OmRwMAAGukamNSkurrfLpvX6deOzmm6ST3UK6F\naH1A7U11amsOqSkSYDNxAABqTFXHpCQF/B59aF+nTl2c1sDwjEqOY3qkiub12Io11l09+xhiI3EA\nAGpcTZSAbVvas7lZHc0h/fJMnMvetygc9Kmt+crl69ZokC18AADAgpqIyXc1NwT10UNdOn8pqXeG\nE8rmWZyzmFDQp5aGoFqjQbVEg6qv85keCQAAlKmaiklJ8ti2tnVH1dMe0dmRGZ27NKN8oWR6LKMi\nIb9aGq6EY0tDkJXXAABg2Wq2GnxeW7t6mrS9O6qReErnx2aVqIFFOj6vrcb6gJojATU3BNUUCcjP\nvo8AAMClmo3Jd3k9tnraI+ppjygxl9VIPKXx6bRmUznTo61YwOdRfcinSMivxnq/miJBNYR8rLgG\nAACrpuZj8r0a6wNqrA9oz5ZmzWcLujw9r/HptOKJ+bK+FB70exW5Go1XPl75NU+aAQAAa42YvIG6\ngHfhjGXJcZRIZjU3n1dqPq9UpqBU5srH9WDblgI+j/w+j+r8V8821v3/cPR5iUYAAGAGMbkM7z4a\nsLkh+L7PRRtDGhyaViqTVzpTUDpbUKnkyHEcOY7kSHIcR6V3f3/1Y+nqJ71e+2ooXv3ofc+vfR4F\nfB75vGzFAwAAyhMxuUJ+n0dNkYCaIgHTowAAAKw7TnkBAADANWISAAAArhGTAAAAcI2YBAAAgGvE\nJAAAAFwjJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArhGTAAAAcI2YBAAAgGvEJAAAAFwj\nJgEAAOAaMQkAAADXiEkAAAC4RkwCAADANWISAAAArlmO4zimhwAAAEBl8t7sBfF4cj3mqFixWIRj\ntASOz81xjJbG8bk5jtHSYrGI6RGAqsZlbgAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8Qk\nAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMm\nAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4Box\nCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeI\nSQAAALhGTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhG\nTAIAAMA1YhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1\nYhIAAACuEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACu\nEZMAAABwjZgEAACAa8QkAAAAXCMmAQAA4BoxCQAAANeISQAAALhGTAIAAMA1YhIAAACuEZMAAABw\njZgEAACAa8QkAAAAXCMmAQAA4JrlOI5jeggAAABUJs5MAgAAwDViEgAAAK4RkwAAAHCNmAQAAIBr\nxCQAAABcIyYBAADg2v8DSm3DUEfeZiIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1115399e8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Approximate the true distribution using a diagonal covariance Gaussian from the class\n",
    "$$\\mathcal{Q} = \\left\\{\\left.N\\left(\\begin{pmatrix} \\mu_x \\\\ \\mu_y \\end{pmatrix},\n",
    "                              \\begin{pmatrix} \\sigma_x^2 & 0 \\\\ 0 & \\sigma_y^2\\end{pmatrix}\\ \\right|\\ \n",
    "                       \\mu_x, \\mu_y \\in \\mathbb{R}^2, \\sigma_x, \\sigma_y > 0\\right)\\right\\}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vi_e = Ellipse(np.zeros(2), 2 * np.sqrt(5.991) * approx_sigma_x, 2 * np.sqrt(5.991) * approx_sigma_y)\n",
    "vi_e.set_alpha(0.4)\n",
    "vi_e.set_facecolor(red)\n",
    "vi_e.set_zorder(11);\n",
    "ax.add_artist(vi_e);\n",
    "\n",
    "vi_rect = plt.Rectangle((0, 0), 1, 1, fc=red, alpha=0.75)\n",
    "\n",
    "ax.legend([rect, vi_rect],\n",
    "          ['Posterior distribution',\n",
    "           'Variational approximation'],\n",
    "          bbox_to_anchor=(1.55, 1.));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqkAAAHBCAYAAABKaSJtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcHGWdP/BPHX13z91zZyaTgySQkATBcLoYRJAF5RAh\nASGAIgoYEVzDAnJG3HCpwUU0i/6ICpIlgCggRxZRTpUzhFyTa+57pu+uruP3x5A2Q+ae7q7q7s/7\n9drNzHR31XdqpOczTz3P9xEMwzBARERERGQhotkFEBERERF9EkMqEREREVkOQyoRERERWQ5DKhER\nERFZDkMqEREREVkOQyoRERERWc6oIZXdqYiIiIjIDPJoDwqCgK6uYKZqyUp+v4/XaBS8PmPjNRod\nr8/YeI1G5/f7zC6BiCaBt/uJiIiIyHIYUomIiIjIchhSiYiIiMhyGFKJiIiIyHIYUomIiIjIchhS\niYiIiMhyGFKJiIiIyHIYUomIiIjIchhSiYiIiMhyGFKJiIiIyHIYUomIiIjIcmSzCyAiIso1uq6j\nr68vo+csLi6GKI5/7Gnp0qVoaWlJfi7LMiorK3Heeefh8ssvn3I9TU1NaGxsxIknnjip13/1q1/F\nEUccgWuuuWbKtQxn2bJlOPbYY3H11Vdj1apVUFUVd99996ivMQwDjz76KM4777wRr/WcOXPwq1/9\nCsceeyyWLl2Kb37zmzj33HMnVeNzzz2HT33qU/D7/Vi7di1ee+01PPLII5M6VjZiSCUiIkqxvr4+\n/L8//AMeb0FGzhcOBXDxF49EaWnphF63atUqnH766QAAVVXxxhtv4IYbbkB5eTnOPPPMKdX0n//5\nnzjiiCMmHVLXrl0Lm802pRrG64YbbhjX8/7+97/jlltuwbnnnjtiSP3b3/6GwsLCKdfU0tKClStX\n4vnnnwcAXHrppfjqV7865eNmE4ZUIiKiNPB4C+ArKDa7jFF5vV74/f7k52eddRb++Mc/4vnnn59y\nSJ2qoqKijJ3L5/ON63mGYYz5nAOv51R88lwejyclx80mnJNKRERESbIsJ0cwdV3HunXr8LnPfQ6H\nH344LrzwQmzdujX53Oeeew6nnXYaFixYgM9//vN4/PHHAQyO0L711lv4+c9/nhz9a29vx7e+9S0s\nWrQIJ554Iu6++24oigIA2LhxI77yla/g29/+Nj71qU9hw4YN+OpXv4r77rsvea6NGzfitNNOw+GH\nH46zzz4bb775ZvKxpUuXYs2aNTj++ONx2mmnQVXVg76vF154AaeccgoWLVqEO+64Y0gIXLVqFa67\n7joAQDAYxHe+8x18+tOfxhFHHIGrrroKXV1daG5uxkUXXQQAOOyww/Dmm29i1apV+P73v48zzzwT\nS5YswbZt2zBnzhy89tpryWPv3LkTZ511FhYsWIBLLrkEzc3NAIDm5mbMmTMHe/fuTT537dq1WLZs\nGQDgpJNOAgB8/vOfx8aNG4c8BgDvvPMOli1bhkWLFmHp0qX47W9/O+T7ueOOO/Dd734XixYtwmc+\n8xls3LhxrB+95TCkEhERERKJBJ5//nm8+uqryYD0s5/9DA899BCuv/56PPHEE6itrcXXvvY1hEIh\n9PT04LrrrsOKFSvw3HPP4Rvf+AZuvPFGNDY24oYbbsDixYtx8cUXY+3atTAMA1deeSUKCwvx+OOP\n4+6778bLL7+Me++9N3n+9957D/X19diwYQM++9nPDqlt48aNuO2223D55ZfjqaeewnHHHYfLL78c\nra2tyef84Q9/wLp163DPPfdAlofeKN65cye+853vYNmyZXj88cehKAreeeedYa/DT37yE7S0tGD9\n+vV47LHH0NPTgzvvvBNVVVVYu3YtAOCVV17B4sWLk+e98sor8ctf/hKzZ88+6Hi///3vcdlll+Hx\nxx+Hruv43ve+N66fx4YNG5KvP+2004Y81tjYiIsvvhhHHXUUnnjiCVx99dW466678Oyzzyaf8+ij\nj2LevHl4+umnccopp+CWW25Bf3//uM5tFbzdT0RElKduu+02/PCHPwQAxGIxOJ1OXHzxxfjiF78I\nwzDwm9/8BitXrkyG1ttvvx0nn3wynnrqKSxevBiJRAIVFRWoqanBOeecg+rqapSVlcHn88Fms8Hl\ncqGoqAivv/46mpub8dhjj0GSJADAD37wA1x66aXJEUwAuOKKK4a9rb1+/XpccMEFySkI1157Ld56\n6y2sX78e3//+9wEAZ5xxBubOnTvs9/n444/jiCOOwIoVKwAAN910EzZt2jTsc1taWuB2u1FbWwuP\nx4M1a9YgEAhAkqTkXNPS0tJkEJ43bx5OPvnkEa/x+eefn5z3u3r1apx00knYvn073G73iK8BgJKS\nEgCDC+KcTueQxx577DHMmTMH3/3udwEADQ0NaGxsxLp16/CFL3wBAHDIIYfg61//OgBg5cqVePjh\nh7Fjxw4cddRRo57XShhSiYiI8tRVV12FU089FQDgcDjg9/uTIbKnpwf9/f1YuHBh8vk2mw3z589H\nY2Mjli9fjqVLl+Lyyy9HXV0dPvvZz+Lss88edtFQY2MjAoEAjjzyyOTXDMNAIpFIjoYWFRWNOO+y\nsbER3/zmN4d8bdGiRdi1a1fy85qamhG/z8bGRsyZM2fI93Hg5wdasWIFvvnNb+KYY47BkiVLcPLJ\nJ486P7e2tnbExwBgwYIFQ55bVFSExsbGIV+fqMbGxiE/FwBYvHjxkFv+06ZNS37s9XoBYNhpEFbG\nkEpERJSnSkpKUF9fP+xjnxy920/TNGiaBkEQ8MADD+DDDz/Epk2bsGnTJvzud7/Dz3/+cxx//PFD\nXqOqKurr6/Hggw8edLzKykoAgyF5JMPVsr+O/UZ7/XBG6hywZMkSvPLKK9i0aRP+8pe/4Ec/+hGe\nfvpprF+/ftjn2+32Uc8jCMKQz3Vdh81mO+jrwPhD5HDXQ9f1IddjuO9vPAu/rIRzUomIiOgg+1f+\nv/fee8mvJRIJfPjhh8nbyz/60Y9w2GGH4eqrr8YTTzyBI488Ei+88MJBx2poaEB7ezuKiopQX1+P\n+vp6dHV14Z577hlXcJoxY8aQOoDBOawNDQ3j+l5mz56NDz74IPm5pmnYtm3bsM/99a9/jffeew9f\n/OIXcc899+AXv/gF3nrrLXR3dw8bLMeyffv25Md79uxBIBDAzJkzkyEyHA4nH9+/qAo4ONweaLjr\n8c4774z7emQLhlQiIiIa1qWXXor7778fL730EhobG/GDH/wA8Xgcp59+OgoKCvDoo49i7dq1aGpq\nwhtvvIFt27Zh/vz5AAZbJu3btw89PT04/vjjUVtbi+uuuw5bt27FO++8gxtvvBGiKI5rBPSSSy7B\n7373Ozz55JPYvXs37rnnHmzduhVf+cpXxvV9nHvuudiyZQvuv/9+7Nq1C3feeSfa29uHfW57eztu\nv/12vP3222hqasLTTz+N6upqFBcXJ+eRbtmyBfF4fFznfvjhh/HnP/8ZW7duxfXXX4/PfvazaGho\nQFlZGaqqqvDQQw+hqakJTz75JF5++eXk6/afa+vWrUOCLAAsX74c27dvx7333ovdu3fjySefxO9+\n9ztceOGF46opW/B2PxERURqEQ4GsP9eKFSsQCoVw8803IxgMYtGiRXj44YdRVlYGYLBl0j333INf\n/vKXKCwsxLJly/DlL38ZAHDeeefh+9//Pr72ta/hiSeewAMPPIDVq1fj/PPPh8PhwMknn4xVq1aN\nq45TTjkFXV1d+OlPf4quri7MmzcP//M//zPsavrhTJ8+HT//+c9x55134he/+AVOPvlknHDCCcM+\nd+XKlQiFQrjyyisRDoexcOFCPPDAA5AkCYcccgiOP/54LF++fEhngtFcdtllWLt2Lfbt24cTTjgB\nt99+OwBAFEWsXr0at99+O0477TQsWbIE3/rWt/Diiy8CGFwwdfbZZ+Paa68dsrgMGJwi8eCDD2LN\nmjV46KGHUF1djVWrVk16ZyurEowxxtm7uoKZqiUr+f0+XqNR8PqMjddodLw+Y+M1Gp3fP75G7amU\nDduiElkdR1KJiIhSTBTFCW9RSkRD8U8uIiIiIrIchlQiIiIishyGVCIiIiKyHIZUIiIiIrIchlQi\nIiIishyGVCIiIiKyHIZUIiIiIrIchlQiIiIishw28yciIkoxXdfR2tqa0XNWV1ePe8ep5cuXo7y8\nHD/+8Y8Peuzll1/GVVddhb/+9a8oLi6eUA0bNmzAAw88gE2bNo35XF3X8cgjj+CCCy4AAFx33XWQ\nZRk/+tGPJnTOiXrttddwySWXYNu2bWk9T7ql83q9/vrrKC8vx8yZMyf0M001hlQiIqIUa21txV9v\nugVlXm9GztcdCuGE229BbW3tuJ5/xhlnYM2aNYjH43A4HEMee+aZZ3D88cdPOKDuP+5JJ500rue+\n+eabuO2225Ih9eabb57w+fJZuq6XpmlYsWIFHn74YcycOXNCP9NUY0glIiJKgzKvF5UFBWaXMaxT\nTz0Vq1evxl//+ld87nOfS35dURRs2rQJt91226SO63Q64XQ6x/VcwzCGfO7z+SZ1znyVruv1yZ/L\nRH6mqcY5qURERHmmuLgYxx9/PJ577rkhX3/llVeg6zqWLl0KAAiFQrj++utxzDHHYP78+Tj11FPx\n4osvAgBUVcWcOXPwk5/8BEuWLMFll12GDRs2JF8LAC+99BLOPPNMLFiwAEceeSS++93vIhwOY+/e\nvbjkkksAAHPmzME//vEPXHfddVi1atVBrz388MNx2mmn4fnnn08+tmzZMjzwwAO47LLLcPjhh+OU\nU07BX//61+TjO3fuxGWXXYbFixdjwYIFWL58ORobG8d1bUaqGQDuu+8+XHPNNbj++uuxcOFCnHLK\nKXjppZeG1HX//fdj2bJlWLhw4ZDzDne9AOCf//wnzj//fCxatAhLly7FI488AgCIxWI4+eSTcf31\n1yePf+ONN+L000+HoihDrtd9992H733ve7j11luxePFinHTSSXjzzTfx61//GkcffTSOOeYYPPro\no+O6Pvt/fhdddBH++7//+6Cf6Y4dO3DZZZfhiCOOwAknnID7778fuq4n67j22mtx880344gjjsAx\nxxyDX/7yl+O67sNhSCUiIspDp59+Ol5++WUoipL82rPPPovPf/7zyZGzO+64A3v37sWvfvUr/PGP\nf8TixYtxww03DHnN//3f/+GRRx4ZEjABYO/evVi5ciUuuOACPPvss7jvvvvw6quv4rHHHkNtbW1y\nPuzf/vY3LFy4cMhr//a3v2HlypU4++yz8dRTT+Gcc87BNddcg82bNyef8+CDD+KMM87AH//4Rxxy\nyCG46aaboOs6dF3HFVdcgbq6Ojz11FN45JFHoCgK7rrrrjGvyWg17/fCCy8AADZu3Iizzz4b3/72\nt7Fz587k47/4xS/whS98AY8//jjKyspw+eWXj3i9tm/fjhUrVuDoo4/Gxo0bceWVV+LOO+/ECy+8\nAKfTiVtvvRVPPvkk3n33Xbz55pt44okncOedd8Jutx9U+7PPPovCwkI89dRTmDNnDq666ir84x//\nwG9/+1uce+65WL16NQKBwJjX53//938BAGvXrsWKFSuGnKOnpwcXXnghqqqqsGHDBvzgBz/A+vXr\n8fDDDyef8+c//xlutxtPPPEELrnkEtx9993YvXv3mNd+OAypREREeeikk06Cpml47bXXAADxeByb\nNm3CGWeckXzOkUceiVtvvRVz587F9OnTcemll6K/vx89PT3J55x33nmYMWMGZs+ePeT4mqbhpptu\nwrnnnova2lqccMIJWLJkCXbs2AFJklBYWAgA8Pv9sNlsQ177m9/8BqeccgouuugiNDQ04LLLLsNJ\nJ52Ehx56KPmcE088EWeeeSbq6upwxRVXoK2tDV1dXYhEIjjvvPPwH//xH6irq8P8+fNx1llnDQmS\nIxmt5v2Kiopw6623YubMmfjGN76BBQsW4PHHH08+/m//9m+46KKLMGvWLKxevRo9PT149dVXh71e\njz32GA477DB85zvfwYwZM3DOOefgggsuwLp16wAAxx57LL70pS/hjjvuwE033YRLL70UCxYsGLb2\n4uJirFy5EnV1dTjrrLMQCARwww03YObMmVixYgUURUFzc/OY16ekpAQAUFhYCLfbPeQcf/jDH+B2\nu3HLLbdg5syZOPnkk3H11Vcn691/fb73ve+hvr4el19+Obxe75A/LiaCc1KJiIjykMvlwkknnYQ/\n//nPOPHEE/GXv/wFHo8HRx99dPI5Z599Np5//nk8+uij2LVrFz788EMAg2Fuv5qammGPP2PGDDgc\nDjzwwAPYsWMHduzYgZ07d+JLX/rSmLXt2rULF1544ZCvLV68GE8++WTy87q6uuTH3o8XqKmqCq/X\ni2XLlmHjxo3YvHkzdu/ejQ8//BB+v3/M846n5vnz5w8ZyZw/f/6QqQSLFy9Ofuzz+VBXV4fGxkac\ncMIJAIZer8bGxoNGkRcvXowNGzYkP1+1ahVOOeUUeL1eXH311SPWXltbC0EQAAAOhwOiKKKysjL5\nOTA453gq16exsRGHHXYYZPlf8XHx4sXo6upCIBBI1nFglwmPxwNVVcc89nA4kkpERJSnzjjjDGza\ntAmqquKZZ57BaaedBkmSko9fe+21uOuuu1BYWIjly5fjgQceOOgYn+wOsN+HH36If//3f8fOnTtx\n5JFH4oc//CFOPfXUcdU13DE1TUvOfQRw0OgrMLjoJxQK4ZxzzsGf/vQnzJw5EytXrsR11103rvOO\np+YDrw8w2ErrwFA21uMHfm8jfZ8H/hHQ0tKCcDiMjo4O7NmzZ8TaDwyOACAIQjK0Hmgq12e4BVT7\na93/70g/l8ngSCoREVGeOu644yCKIl5//XX85S9/wfr165OPDQwM4JlnnsGGDRtw+OGHA0BykdB4\nQseTTz6JT3/607jnnnuSX9uzZw8OOeQQABg2QO03Y8YMvPfee0O+9u6776KhoWHM877xxhvo7u7G\nn/70p2Rwe/nll1NSMwBs27ZtSPDcvHkzjjnmmOTjW7duTX48MDCApqYmzJkzZ8Tv8+233x7x+9Q0\nDTfeeCO+/OUvIxqN4sYbb8Sjjz467n64wxnr+oz2c2loaMBLL70EVVWTr3333XdRUlKCoqKiSdc0\nEo6kEhER5SlZlvGFL3wB99xzDyoqKjB//vzkY06nEy6XC88//zyam5vxyiuvYPXq1QAwZCHQSIqK\nirBt2za8//772L17N1avXo0tW7YkX7t/vuPmzZsRj8eHvPaSSy7Bc889h/Xr12PPnj146KGHsGnT\nJixfvnxc541EInjxxRfR3NyM3//+98nFQVOtGQCamppw1113YdeuXfjZz36GrVu34stf/nLy8aef\nfhpPPvkkdu7ciRtuuAG1tbVYsmTJsOe74IILsGXLFvz4xz/G7t27sXHjRjzyyCPJqQ6//vWv0d7e\njmuuuQbXXXcddu3ahd/85jdjfh9TuT6SJMHhcGDHjh0IBoNDXvulL30J8Xgct9xyCxobG/Hiiy/i\nZz/7GZYvXz5quJ0sjqQSERGlQXcolBXnOuOMM/Db3/4W3/72t4d83eFwYM2aNVizZg3Wr1+P2tpa\nXHnllbj33nvx0Ucfob6+ftTjXnzxxdi2bRsuvvhiOJ1OHHXUUfjWt76FZ555BgAwb948HHfccTj/\n/PPxk5/8ZMhrFy1ahDVr1mDt2rX4r//6L8yYMQM//elPh8yXHcmRRx6JK664Arfeeivi8Tjmzp2L\nm2++GTfccAM6OjqmVDMALFy4EL29vTjzzDPR0NCAdevWYdq0aUOu56OPPoqPPvoIRx11FNatWwdZ\nloedl1lTU4MHH3wQa9aswbp161BTU4Mbb7wRZ599NpqamrB27VrcdNNNyUVmK1euxL333jukt+1E\njXV9KioqsGLFCqxZswYtLS2YMWNG8rVerxfr1q3D6tWrceaZZ6K0tBSXXnopvv71r0+6ntEIxhjj\n311dwdEeznt+v4/XaBS8PmPjNRodr8/YeI1G5/dnvkm81bdFpcm577778Pbbbw+ZFnGgZcuW4dhj\njx11gRONH0dSiYiIUkwUxXFvUUpEw+OfXERERERkORxJJSIiIhqHa665ZtTH929pSqnBkVQiIiIi\nshyGVCIiIiKyHIZUIiIiIrIchlQiIiIishyGVCIiIiKyHIZUIiIiIrIchlQiIiIishyGVCIiIiKy\nHIZUIiIiIrIchlQiIiIishyGVCIiIiKyHIZUIiIiIrIchlQiIiIishyGVCIiIiKyHIZUIiIiIrIc\nhlQiIiIishyGVCIiIiKyHIZUIiIiIrIchlQiIiIishyGVCIiIiKyHIZUIiIiIrIchlQiIiIishyG\nVCIiIiKyHMEwDMPsIoiIiIiIDiSP9YSurmAm6shafr+P12gUvD5j4zUaHa/P2HiNRuf3+8wugYgm\ngbf7iYiIiMhyGFKJiIiIyHIYUomIiIjIchhSiYiIiMhyGFKJiIiIyHIYUomIiIjIchhSiYiIiMhy\nGFKJiIiIyHIYUomIiIjIchhSiYiIiMhyGFKJiIiIyHIYUomIiIjIchhSiYiIiMhyGFKJiIiIyHIY\nUomIiIjIchhSiYiIiMhyGFKJiIiIyHJkswsgIspnqqYDAAQBEAQBoiCYXBERkTUwpBIRpZCS0BBP\naFAS+uC/qoa4okFR9Y//1RBP6FASGpSEBk03DjqGKAjJ0Hrgv6IgwGGX4HHa4HHKyX/dXicMw4DA\ngEtEOYQhlYhoggzDQCSuIhhJIBhRkv+GogkkVH3Kx9cNAzCAj//fENG4iv5gfMjXPLt6EYsqcO8P\nry4bSguc8Be5YJM5q4uIshNDKhHRCHTDQDiaGBJGQ9EEgtEENG3qYTSVNN34uEYFANDYMgBREFBc\n4EBFsRsVxS4Ueh0mV0lENH4MqUREGBwdDUQS6AvG0BeIoy8URyiagD7M7fhsoRsGegZi6BmIYcse\nwGmXUVXqxvRKHwMrEVkeQyoR5SUloaE3GEdfIIbeYBz9oXhKbtVbWUxRsbstgN1tART7HJheWYAa\nvweyxCkBRGQ9DKlElBdiijo4qhgYHFkMRBIwjOwdJZ2qvmAcfcEufLS3D3PqilBf6WNnASKyFIZU\nIspJ8YSGrr4ouj8OpvvnatJQMUXFezu7sbNlAPPqi1Hr95pdEhERAIZUIsohA6E42nsj6OiLoi8Y\nz+uR0okKRxP4x9ZOtHSFsWhWGRx2yeySiCjPMaQSUdZSNR1d/VF09EbR0RdBNK6aXVLWa+sJozcY\nw6JZZagq9ZhdDhHlMYZUIsoqkVgC7b1RdPRG0D0QHbYZPk1NXNHw5pYOHDq9BIdMKzK7HCLKUwyp\nRGRpum6guz+Kjr4o2nsjnFuaQVv29CKmqFgwo5S7WRFRxjGkEpHl6IaBrr4omrtCCMbb0D8QNbuk\nvLWrNQAloeNTc/wMqkSUUQypRGQJxseN55u7w2jrDiOe0AAAHg+bzputuSsEj1PGvOklZpdCRHmE\nIZWITNUXjKOlO4SWrjAXPlnYtqZ+eFw21FX4zC6FiPIEQyoRZVwwoqClK4zmrhBC0YTZ5dA4vbez\nG0VeBwo8drNLIaI8wJBKRBkRjavJYNofiptdDk2Cpht4Z0cXPrOwmvNTiSjtGFKJKG3iCQ2t3YPB\ntDfA5vq5oC8Yx572IBqqCswuhYhyHEMqEaVcV38Ue9uDaO0JQ2cf05yzdV8f6iq8kETR7FKIKIcx\npBJRSsQTGpo6QtjTHuA80xwXVzQ0dYYwvZKjqUSUPgypRDQl+0dN23rC3P0pjzS2BFBf4ePcVCJK\nG4ZUIpowjppSMKKgayCG8iKX2aUQUY5iSCWicevuj2IPR03pY+09YYZUIkobhlQiGhVHTWkkHb1R\nYKbZVRBRrmJIJaJh9QzEsLstwFFTGlE4lkAgoqDAzeb+RJR6DKlElKQbBtq6w2hsDaA3EDO7HMoC\n/cE4QyoRpQVDKhEhoerY2xHErtYAIjHe0qfxC0b4vxciSg+GVKI8Fomp2NU2gL3tQSRU3exyKAsF\nI4rZJRBRjmJIJcpDfcE4GlsG0Nodhs6tSmkKuJiOiNKFIZUoTxiGgfbeCHa2DKBngPNNKTU4Ak9E\n6cKQSpTjVE1HU2cIjS0DHPWilFM1hlQiSg+GVKIcFVNU7G4NYE97EPGEZnY5lKM03YBuGBC5PSoR\npRhDKlGOicQS2N40gKbOIPubUkbougFRYkglotRiSCXKEeFYAtv39aOpM8TFUJQxgiBAEhlQiSj1\nGFKJslwomsD2pn40M5ySCWRJgMBb/USUBgypRFkqGFGwvakfLV1sI0XmkSXR7BKIKEcxpBJlmUBE\nwfZ9/WjpDsNgOCWT2WSGVCJKD4ZUoiwRCCvY1tSPVoZTshC3k79GiCg9+O5CZHEDYQXb9vWhrSfC\ncEqWU+C2m10CEeUohlQiixoIxbGtqZ/hlCzNx5BKRGnCkEpkMf2hOLbu60N7T8TsUojG5HPbzC6B\niHIUQyqRRYRjCXy0p48LoihrCIIAr4shlYjSgyGVyGSxuIr3G7uxpz0InTtEURZxO2W2oCKitGFI\nJUoRQ9dhJBKD/6eqMNR/fawnVMDQAcMADAAwoKoaWrvD6AkpCIXjcEAA9jdFFwQYogRdkmFIMgxR\nGvxXkqCLMiBJZn6rRACAIq/D7BKIKIcxpBKNg6Gq0KIR6NEY9GgEejQKPRYbDKOJxGAgHecoqAED\nvQMxtPdGkdA0OOwy7Io6wYoEGNJgcNVFGYZsg2ZzQLc7oNsc0OxOGDJvw1J6VRS7zC6BiHIYQyrR\nxwxNgxYODwbQaBR6LPpxII1BTyRSco7+UAxtPRHEE9pUq4WgqRA0Fftvth4USQXxoODKAEupIgoC\nKkvcZpdBRDmMIZXykmEY0MNhaMEA1GAQWiAALRL++FZ86oWiClq7w4jEJzpiOgWGDkmJQlKiBz2k\ny3ZoTg9UpweqywPN4QZEzi2k8SstdMJu47QTIkofhlTKC3o8PhhIA0FowQC0YBCGNtXRzLFFFRVt\n3WEEIkqHXhiSAAAgAElEQVTazzURoqpADCmwhfo+/ooAzeGC6vQMhleXB7rdaWqNZG1VpR6zSyCi\nHMeQSjnJUFUk+nqh9vRAHeiHHotn9PxKQkN7bwR9wTiMdA3PppQBKR6BFI8AA12DXxElqC4vEp5C\nJDyFMGxs2k6DBEFAVSlv9RNRejGkUs7QImGoPb1I9PZAGxgwpdeopuno6IugeyAGPct7nQq6Blt4\nALbwAABAs7sGA6u3EJrT869OBJR3Cr12uBz89UFE6cV3Gcpahq5D7e+H2tuDRG8P9GjMvFpgoLs/\nhvbeMLQc7XW6f36rs68dhihD9fj+Ncoq8a0kn1RxwRQRZQB/s1BWMXQdak8PlM4OqH29MDTd7JIQ\njCpo6QojNuE2UtlL0FXYgn2wBfsACFCdHiR8xVAKShhY80BVGeejElH68bcJZQUtFILS3oZEZ2fK\n2kFNlZLQ0NoTRn8os/NdrceAHAtBjoXg6m5BwlOIeEEpVE8BpwTkIK/LhgI35ycTUfoxpJJl6QkF\nic5OKO1t0EJhs8tJ0nUdnf1RdPZFs37eacoZOmyhPthCfTAkG5SCEsQLy9gpIIfUV/rMLoGI8oRg\nmLG6hGgEhq5D6elBtKUVsc6uwa1ELaRnIIamjiCUKTfjzy+ay4tEURnUwhKA0wGylk0W8aXPzGR/\nVCLKiDF/W3R1BTNRR9by+328RqMY7/UxVBVKWyvizc3QFWv1FAWAmKKipSuMYDT1tTnsMuK5Pp9V\n6Ycw0A+bIEIpKEGsuBK6fXz7vns8DoTD+T6lYnSZukazagox0B9J+3lSze/n6C9RNuKQBplKTyhQ\nmpsRb22BoVpvdFLTdLT3DraUyo5+pxZn6LAPdMM+0APFV4xYSSV0B/d/zwaiIGBGdaHZZRBRHmFI\nJVPosRjiTU1QOtossUL/YAZ6BmJo641AtWR92c6APdgLe7AXCU8hYiWV0Fxes4uiUVSXeeB28lcG\nEWUO33Eoo7RwGPGmfUh0dprSbH88wtEEWrpDiMRz/Ba8RezfMEB1ehErqYTq5WidFc2q5c+FiDKL\nIZUyQgsGEdu7B4meHrNLGVFC1dDaHUZf3reUMoccC8HbuhOa3YVYaRUSvmKzS6KP+YtcKPKObw4x\nEVGqMKRSWmnxOCLbt0Fpb4NVp3QaMNDd9/FuURYd3c0nkhKFp20X1D4vouXTAA/Dkdlm1nAUlYgy\njyGV0sLQdSgtLejua4fSZ50ep58UiSXQ1BVClLf2LUeOheDbtxVivBIRjx+GbDO7pLzkc9tRUczF\nbUSUeQyplHKJnh7EGndCi0Yh+6w5CqbpOtp7uGrf+gzY+rpQ0NWJWEkV4sXl3MUqw2bVFELgNSci\nEzCkUspokTBijY1I9PaaXcqo+kNxtHSFkdCs1/KKhifoGlzdzXAMdCPir+XiqgzxOG2oLfeYXQYR\n5SmGVJoyQ1UR27sHSkuLZVfsA4CS0NDSFcJAxHqbBdD4iIkYvK07oboLECmfxu1W0+ywhhJIomh2\nGUSUpxhSaUrU/j5Etm2FHrPuingDBrr7Y2jv4cKoXCFHAijY+xGiZdWIF1eYXU5O8he5UF3GUVQi\nMg9DKk2KoWmI7d6FeEuL2aWMKhpPoKmTPU9zkqHD1dUMW6gfkcrp0G3WnP+cjURBwPwZpWaXQUR5\njiGVJkwNBhDduhVaxLp7eBuGjrbeKLr6olwYlePkaAi+PR8hWl4LpbDM7HJyQl2lD4Ueu9llEFGe\nY0ilCYk37UNs925Lzz0NRRU0dYYQT3BhVL4QDA3ujr2whQMIV9QDkmR2SVnLJouYV8+NFIjIfAyp\nNC66oiCy9SOofX1mlzIiTdfR2h1Bb4BtpfKVLdSHglgE4aoGaC7Op5yMOXXFcNgY8onIfAypNCa1\nvw+Rj7ZAVxJmlzKigXAczZ1sK0WAqMbha9o2uKiqpNLscrKKz23HjKoCs8sgIgLAkEpjUNpaEd2x\nw7K391VNR3NXCP0h63YXIDMYcHW3QFJiiFTUcwOAcZrfUAJR5LUiImtgSKURRXftQrxpn9lljCgQ\nHpx7ytFTGok90AMxoSBUPZPzVMdQUeJGRYnb7DKIiJLYpZkOYug6Ih9tsWxA1XQdTR1B7GobYECl\nMcnRIAqatkJMcLR9JKIgYH5DidllEBENwZBKQ+gJBeH33oXS2Wl2KcMKRRVs29ePnmDM7FIoi4hK\nDL592yBFw2aXYkmzagvhc7PlFBFZC2/3U5IWiSC8+X3oUesFQMPQ0doTQXc/V+7T5AhaAr7m7QhX\nTkfCxxZL+xV6HZhbx+tBRNbDkEoAAHWgH+HNm2Go1tuZKRJLYF9nCDHFerVRljF0eNp2I5pQEC/h\ndqqSKOCIQ/xcLEVElsSQSlD7+xHe/D4MTTe7lCEMw0BHXwQdvdw1ilLJgKu7GYKhI1ZaZXYxpppb\nX8ydpYjIshhS85waDCC8+QPLBdSYomJfRxCROEdPKT2cPa0wRAnx4nKzSzFFaYETs2oKzS6DiGhE\nDKl5TAuFEH7/fRiWWiFvoLMvivbeCHSL9mal3OHqaoIhilAKy8wuJaNkScQRc/wQ2D+WiCyMITVP\nadEowh+8b6k5qEpCw76OIEIx6+5sRbnH3bEPhijl1WKq+TNK4HHazC6DiGhUDKl5SI/HEX7/PeiK\nYnYpST0DUbR2h6Fx9JQyzoCnbTdCogjVk/u3vytK3Jheya1Picj62Cc1z+gJZTCgxqzRZkrVNOxq\nHUBTV4gBlUxkwNu6C3IkaHYhaeWwSVg8O7+mNhBR9mJIzSOGqiL8/vvQIhGzSwEABCMKtu7rRyBi\nnRFdymOGDk9rI6SYNf77SIfDZ5bCaecNNCLKDgypeSSyfRu0UMjsMmDAQGt3GLtaA1At1lWA8pug\na/C0NkLQrDNXO1Vq/V7U+L1ml0FENG4MqXki3tKCRFeX2WVASWjY2TyAzv4Ie5+SJYmqAnf7HrPL\nSCmXQ8bhM0vNLoOIaEIYUvOAGgwg1rjT7DLQH4pjW1Mfwly9TxZnCw/A0dtudhkpIYoCjppbDrtN\nMrsUIqIJ4eSkHGeoKiJbtsAwcVGSYeho7o6gZyBqWg1EE+XqboXq8kJzZfct8sMaSlBS4DS7DCKi\nCeNIao6LbNtq6kr+aFzF9qYBBlTKQoOtqbJ5fmqt34uZ1bnfVouIchNDag6LNzch0d1t2vl7AzF8\nuKsHUSV7f8lTfhNVBZ623WaXMSk+tx2L2G6KiLIYQ2qOUgMBxHbtMuXcmq5jX0cA+zqD0HUujqLs\nJkcCWTc/1SaL+PS8csgS3+KJKHtxTmoOMnQd0e1bTZmHGo0nsKc9iHhCy/i5idLF1d2KhLcIuj07\n5nYumu2Hz203uwwioinhn9k5SGlrhRbOfEPyrv4otjcPMKBSDjLg7mwyu4hxmT2tCDVlHrPLICKa\nMo6k5hg9oSC2J7Nz6DRNx77OIAbC3DmKcpccCUAODUD1WnchUkWJG4fWF5tdBhFRSnAkNcfEdu+G\noWZuJDMSS2BbUz8DKuUFd1cToFtzlzSvy4Yj5/ghCILZpRARpQRDag5RgwEo7W0ZO19PIIodLQNQ\nMhiKicwkJuJw9HeaXcZBbLKIJYdWwCazYT8R5Q7e7s8hsZ07kYmdRg3DQHNXGD0B9j6l/OPqaYPi\nK4Fhs8bCJEEQcMyCKjg4gEpEOYYjqTlCaW+HGgik/TwJVcPOlgEGVMpfhg5Xd4vZVSTNrStCbbnP\n7DKIiFKOITUHGLqO2O7090QNRRVsbxpAOJZI+7mIrMwe7IUUC5tdBuorfZhTx4VSRJSbGFJzQKKr\nE7qS3oVL3QNRNLYGkNA4/5QIABx95s5NrSxxY+Es7ihFRLmLITUHxFvSd+vRMHTs6wiiuStkyuYA\nRFZlD/ZBUM25q1Dsc+DIueUQuZKfiHIYQ2qWUwcGoAWDaTm2ktCwo3kAvcFYWo5PlN0MOAa6Mn5W\nj8uGow+t5JanRJTz+C6X5ZTW9IyihqIKtjf3IxJX03J8olzg6O8GMniHwWGXcMxhlXDY2WqKiHIf\nQ2oW0+NxJLpSP5LT1RdFY0sAqmbNpuVEViFoCdiDvRk5lyyJOPrQSnhdtoycj4jIbAypWUxpa03p\nPFFd17G3PYCWnhCMTDRcJcoBjr703/IXBQFHzS1Hsc+R9nMREVkFm/lnKUPXobS1pux4SkLD7rYA\nogpv7xNNhBQPQ4qGobk8aTvHwlllqChxp+34RERWxJHULKX29kBXUrOyeP/8UwZUosmxB3rSduy5\ndcWor2SzfiLKPxxJzVKJ3tTMg+sNxNDE9lJEU2ILB5COPdjqK32YW89m/USUnxhSs5TaN9WQaqCt\nJ4KOvkhK6iHKZ6Iah6jEoNudKTsmm/UTUb7j7f4spIXD0GPxSb/eMHTsaQ8yoBKlkC08kLJjsVk/\nERFDalZSp3CrX9V07GwJoD80+ZBLRAezhQMpOU6R14FjDmOzfiIi3u7PQom+yS3SiCsqdrUGEFe1\nFFdERHI0BOg6IE4+XBZ67Dh2fiXsNjbrJyLin+pZxtA0aAMTv624fwU/AypRmhg65Ojktygu8Nhx\n7IIqBlQioo8xpGYZtb8fhj6xlfi9gRgaWwPQJvg6IpqYyd7y97ntOG5+FRwMqERESbzdn2W0wERG\nUbmCnyiT5Gh4wq/xumw4bkElHHYGVCKiAzGkZhktOr5ujIahY19HCH1cIEWUMWIiNqHne1w2HLeg\nCk4734qJiD6J74xZRh9HSFU1HbvbAgjHUrMjFRGNj6BrEDQVhjT2W6vHacNx86vgcvBtmIhoOHx3\nzDJ6bPSQGldU7GoLIJ7gAikiM4hKHJpr9LdWt3PwFr/bybdgIqKRcOFUFtETCoxRVueHownsaBlg\nQCUykTTGLX+3Q/44oNoyVBERUXbin/FZRI+MPIoaDCvY3R6AbnAFP5GZxMTI88BdDhnHLqiChwGV\niGhMDKlZZKRb/f2hGPZ2hGAwoBKZTlSGD6lOu4zjFlTB62JAJSIaD4bULDLcoqmegSiau8IwwIBK\nZAXSMCOpDruEYxdUMqASEU0AQ2oW0RVlyOedfRG09ky8LyMRpY+gDu2q4XbIOGZ+JXxuu0kVERFl\nJ4bUbKLryQ9bu0Po7B9fz1QiyhzB+Nd/p16XDcfOr+IqfiKiSeA7ZxYxdB0GDDR3htATmFjTcCLK\nkI/nhhd5HTjmMO4kRUQ0WQypWcTQdOxtD6Kfu0gRWZZgGCgtdOLoQytgkxlQiYgma8yQ6vf7MlFH\nVsvENUqoOjZ3hRBVNDiybAvFbKvXDLxGo8um61Poc+K4E2dDljLbhprv1USUa8Z85+/qCmaijqzl\n9/vSfo2UhIY3tnQg3hmCTVHTeq5Uc9hlxLOs5kzjNRpdNl2fYq8D1SVu9PVmdkFjJt6HshkDPFF2\nyp7hiTwVU1S8vrkdA2EFHlEwuxwiGkFpoQvT/B4IEm/xExGlAkOqhYVjCby2uR3h6GBLG0PgLrZE\nVlRR7EZVqQcAIPCPSSKilGBItahARMFrH7QjdsBtToZUIuupKfXCX+z61xdEjqQSEaUCQ6oFBSIK\nXv2gDXFFG/J13cZm4ERWIUDAtHIvSgqcQ74uOp0jvIKIiCaCIdViQtEEXvug/aCACgC6zWFCRUT0\nSYIgoL7ChyLvwf9Nii7XMK8gIqKJYki1kFA0gb+93zbkFv+BGFKJzCcJAqZXFYy4zanEkEpElBIM\nqRYRjiXw2gcjB1QA0Oy8jUhkJlkSMaOqAG6nbcTniE6GVCKiVGBItYBILIFXP2hHJD5GL0hRhCHZ\nIGiJzBRGRElOu4wZVQWw20ZfGMXb/UREqcGQarJITB0MqLHxBU/N5oDMkEqUUT6XHdOrfJDEsTts\nMKQSEaUGexqZKBpX8ermNoTHGVABQLdzXipRJpX6nJhRXTC+gGq3s5k/EVGKcCTVJNG4ilc/aEs2\n6h8vzkslygwBAqpK3Sgvdo/7NaJ7/M8lIqLRMaSaIKaoeG1zO0ITDKgAoLq4BzVRuomCgLoKL4q8\nE/ujUC4qSlNFRET5hyE1w+KKhtc2tyMYUSb1es3phiHJELQxFlkR0aSMZwX/iK8tKUlDRURE+Ylz\nUjMontDw2uY2BMKTC6gAAEGA6uZoKlE6OO0yDqktmlRAFW02yL6CNFRFRJSfGFIzREloeH1zOwam\nElA/lnAXpqAiIjqQz2XH7NrCMVtMjUQuLk5xRURE+Y23+zNA1XS8saUD/aF4So6X8HC0hiiVSnxO\nTCv3QhCESR+Dt/qJiFKLI6lpphsG/rmtC72BWMqOacg2aA6uIiaaKgECqkrcqKvwTSmgAgypRESp\nxpCaZh809qCtJ5zy43I0lWhqREFAfaUXFSWeKR9L8vkg2uwpqIqIiPZjSE2j7U392N0WSMuxEx7O\nSyWaLJskYVZN4YRbTI14vJLSlByHiIj+hXNS02RfRxBb9vSm7fiaywvN7oKkRNN2DqJc5HHaML3S\nB5ucmp2hBEGAvbIyJcciIqJ/4UhqGnT2RfDuzu60nydeXJ72cxDlktJCF2bVFKYsoAKAXFYG0cmd\n4IiIUo0hNcX6Q3H8fWsndN1I+7kUXwkMiYPhRGMRBQF15T5M809tBf9wHDU1KT0eERENYkhNoUgs\ngTc+7EBC1TNzQlFEvKAsM+ciylJ2eXD+aUlB6kc7Ja8XciG3QiUiSgeG1BRREhpe/7ADMSWz25XG\ni/wAUjsyRJQrvC4bDpk2uR2kxoOjqERE6cOQmgL7m/UHI1PfTWqiDJsdCS9X+hN9kr/QhZk1hZCl\n9LzNiTYbbOUVaTk2ERExpE6ZkYZm/RMVL+ICKqL9REFAfYUPNX4vhDTeZbBXVUEQ+RZKRJQuXHUz\nRf/c2pmWZv0Tobp9UJ1eyLGQqXUQmc0uS2io8sHlSM/t/f0ESYS9mrf6iYjSicMAU9DYMoDt+/rM\nLgMAEC2vNbsEIlP53HYcMq0o7QEVABx19RAdjrSfh4gonzGkTlJnfxQf7k5fs/6J0pweKAXc9Yby\nU3mRGzOqC9I2//RAossJR+20tJ+HiCjfMaROQjiWwD+3dkI30t8LdSKiZTUwxNQ1KSeyOkkQML3S\nh+oyT1rnnx7INWMm56ISEWUA32knSNV0vLWlA/GEZnYpBzFkG2IlVWaXQZQRTruM2dOKUOTN3G5P\ncnExbGX+jJ2PiCifMaRO0Ds7ujEQznyrqfGKF5dDt3OLRsptpT4nDqkthNOeubWfgiDANWtWxs5H\nRJTvGFInYHtTP1q6LL6CXhAQ8XO+HOUmSRQwvcKHaRU+iBm+5W6vqYHk9mT0nERE+YwtqMapoy+C\nj/ZaYyX/WFRPARKeItjC/WaXQpQyboeM+soCOGyZn3ct2u1w1k/P+HmJiPIZR1LHIRJL4O1tXTAs\ntlBqNJGKOhhS+lvxEGWCv8iF2bVFpgRUCIBrzlwIMv+mJyLKJIbUMWi6jre2dlpyodRoDNmGcFUD\nkKEVz0TpIEsiZtcVo6bMC0Ew53/Lzrp62EpKTDk3EVE+Y0gdwwe7etEfjJtdxqSobh9ipVztT9nJ\n67JhzrQiFPvMa5ovFxXBwdv8RESm4P2rUezrCGJPW8DsMqYkVloFORqCHMnu74PyhwABFSUuVJS4\nM9b7dDii3Qb33HmmjeASEeU7jqSOYCCs4L3GHrPLSIlwVQN02W52GURjskkSZtYUoLIkc835hyUA\n7nmHcutTIiITMaQOQ9V0/GNrJzRNN7uUlDAkmfNTyfIK3HbMqSuC12X+H1TO+umQi4rNLoOIKK8x\npA5jy54+BCPWbdg/GZrLi2hZtdllEB1EEATUlHoxo7oAsmT+W5JcXAxHXb3ZZRAR5T3OSf2Ezv4o\ndmf5PNSRxEsqISlx2APdZpdCBABwyBLqK31wO63RLk3yeuCedyjnoRIRWQBD6gESqoZ3t2dXP9SJ\nilTUQTA02ILZsTEB5a4SnxM1fg+kDO8cNRLJ5YJnwUKINmsEZiKifMeQeoD3G3sRiatml5FegoBw\nZQM8ug5beMDsaigPyZKIaeVeFHqssyhJdDrgOXwhRLv582GJiGiQNYYwLKC1O4ymzqDZZWSGICBc\nNQOqy2d2JZRnirwOzK0rtlZAtdsGR1CdTrNLISKiAzCkAogrGt5rzLN5mqKIUM1MaE6P2ZVQHpBE\nAXXlPkyvtMbiqP0EWYZnwUJIbrfZpRAR0SdY57eFid7Z2YW4kl3bnqaEKCFUMwua3WV2JZTDfG47\n5tYVo6TAWiOVgiTBs+BwSF6v2aUQEdEw8j6k7m0Por0nYnYZpjEkGaHa2dDt1goQlP1EQcA0vxcz\nqwtgkyWzyxlCkER45s+HXFBgdilERDSCvA6pkVgCm3fnxq5SU2HINgSnzYHq5IgSpYbHacOcumKU\nFrpgtU0kRLsNnoWL2KyfiMji8nZ1v2EYeHt7NxJqbuwqNVX7R1TdHXthD/aaXQ5lKUEQUFXihr/Y\nZe62piOQPG545h/ORVJERFkgb0PqrrYAugeiZpdhLaKISFUDdJsdzt52s6uhLONyyKir8MFlt+bb\nilxcBM+h8yHI1qyPiIiGyst367iiYeteNrMfSaysBrrNAXfHPgC5u7EBpYYAAeXFLlSWuC27U5O9\nshKu2YdAsMjGAURENLa8DKkf7e3jbf4xKIVl0GU7PG27IOh52PmAxsVhk1BfYZ1tTYfjbGiAs67e\n7DKIiGiC8m5YYSAUx96OPGnaP0WqpwDBaXOgy9yFh4YSIMBf5MKcaUWWDaiCKMA9dx4DKhFRlsq7\nkdQPdvXCMHgLe7x0hwvB+nlwt+/hNqoEAHDZZUwr91o2nAKA6HLCPe9QyD62mCIiylZ5FVJbusNc\nLDUJhiQjXDMLjt4OuLpbwHmq+UkQBFQWu1BebN25pwBgLy8fnH/KBVJERFktb97FNV3Hlt1srTQV\n8ZIKqG4fPG27ICbiZpdDGeRx2jCt3AunRVfuA4MN+l2zZsNeWWV2KURElALW/Y2TYo0tAYRjCbPL\nyHqa041A/Ty4u5phH+g2uxxKM0kUUFXqQVmhE1Zryn8gucAH15x5kNxus0shIqIUyYuQGo2r2N7U\nb3YZuUOUEKmoh+Itgqd9LwSN4T8XFbrtqPF7YbdZa0vTAwmCAMf06XBMq7P0FAQiIpq4vAipW/b0\nQdXYcirVVE8hAtMPhauzibtU5RCbJKHG70aR19q7MkkeD9xz50HycjtfIqJclPMhtS8YR3NXyOwy\ncpYhyYhUNUApKIWrqxmSwoVp2UqAgNJCJ6pK3ZAs3PRekGU4p0+HvaqazfmJiHJYzofUzbt62HIq\nA1RPAYLueXD0d8HZ0wZBV80uiSYgG9pKQQAcVdVwTJ8O0cbevUREuS6nQ2rPQAw9gZjZZeQPQUC8\nuBxKQQmc3a1wcGGV5YmCgMoSN/zFLghWXhhVVATXzFm8tU9ElEdyOqTubGHzeTMYkoxoRR3iRX6U\nBNoBhfNVrajAbUetxRdGiU4nihYdDknkqn0ionyTsyE1FE2gvTdidhl5TXe4EJ0+F4qzHa6uFogq\ne6tagV2WUF1m7YVRgiTCMa0Ojml1cFYUItjFrYyJiPJNzobUnS0DnItqEQlfMRKeQtgDPXD2djCs\nmkQUBPiLXKgodkG06IIjQZJgr66Go6YWosNhdjlERGSinAypcUVDUydX9FuKKEIp8kMpLIMt2Adn\nbzs7AWRQoceOmjLr3toXbTbYa2vhqK7hdqZERAQgR0Pq7rYANPZFtSZBQKKgBImCEsihATh72yHH\n+AdFujjtMmrKPPC5rbkaXnQ64KidBntlFQTJmgGaiIjMkXMhVdV07G4PmF0GjYPqLUTIWwgpGoKr\npw1yhD+3VJFEAZXFHpQVOS25E5PkccNRWwdbeTl7nRIR0bByLqQ2dYYQVzSzy6AJ0FxehGpnQ4pF\n4OjrgD3UDxgcCZ8MAQKKfQ5Ul7khW3BkUi4qgqOmBrYyv9mlEBGRxeVUSDUMA41sO5W1NKcbkaoG\nRDQN9lAf7AM9nAowAW6HjFq/9Rryi04n7BWVsFdWQHS6zC6HiIiyhGDk0BL45s4gXnmnxewyKIWE\neBS2/m7YBnogJBSzy7EkmyyittyHsiIXrHJnX5AkOCsq4KyuhqO0xOxyiIgoC40ZUruyqD/h37d2\noqUrsyNvHo8D4TBbKo0kZdfHMCBHAnAEemALDeTUdACHXUZcmfg2sgIElBU5UVnihmSReZ1yYSHs\nlZWwlflTtkrf7/dl1fuQGXiNRuf3+8wugYgmIWdu92u6js4+Nu/PWYIA1VMI1VMIQVNhC/bBFh6A\nLRLMqcA6Xj6XHTV+D5x2k/8TFgDZ54NcUgpbeQUkF2/nExFRauRMSO3qjyGh5l9YyUeGJA/2XC3y\nA7oOORIcDKzhAYhqbk8JcNgkVJd6UOg1r9G9IEuQi0tgKymFXFoC0WbN9lZERJTdciaktvWEzS6B\nzCCKUL2FUL2FiAIQ49FkYJWjYQC5MeXaJkmoLHGhpMCcllKS2w25pAS20lJIBYVsG0VERGmXEyHV\nMAy09/JWPwG6w4W4w4V4SSWgaYOBNRqEFItAikeRbaFVEgT4i10oL8rsVqai0wHJVwC5sAi20hKu\nyicioozLiZDaG4izNyodTJKSu1sBGJwaEAtDioU//jdi2ekBgiCgrMCJihI3ZCm94VSQJUg+H2Rf\nAaSCAki+Aoh23sInIiJz5URIbevlrX4aB1GE6vZBdfuwv9+AoCYGA2v04+Aaj0LQJ77SPlUECCjy\n2lFZ6oHDlvpm/IIkQnS5IRcUQPL5BgOp223JXamIiCi/5UZI7eGtfpocQ7Yh4S1CwluU/JqgqRCV\nOMREHFIiDlGJDf6biEPQ0hdgfS47qkrdU27GL0giRKcLossF0eWG5HImPxYd5i24IiIimoisD6kD\nYQXhaMLsMiiHGJIMzSVDc3lw0P+yNC0ZXEVVgaipEDQNgqZC0D/+WNcgauq4W2O5nTbUlnngc498\niyW6DBMAAAiBSURBVF0QBAg2GYJsgyDL//rYNvi56HAMBlGni0GUiIhyQtaHVPZGpYySJGiSG5rT\nPfZzdR2Cvj/AahB0HYABGIAAAy6HjJnVBThsth99fRFAEAEBAAQIogBBkpMhNFWN8YmIiLJF1v/m\nC4Q5ikoWJYowRBGGPPT2vcMm4ZBpRWioKoAoCnD5fQhJ3C2IiIjoQFkfUoMRa67OJvokSRIxq7oA\ns2qLYJPZZ5SIiGg0WR1SDcNAkPNRyeJkScT0Kh9m1xTBYU/9in0iIqJclNUhNRJXoWncCpWsSZZE\nNFQXYFZ1IcMpERHRBGV1SA1GOIpK1mOTRTRUFWBmTWFaep0SERHlgywPqZyPStZhk0XMqC7EzOoC\n2BlOiYiIpiSrQ2o0bt7OQET72WQRM6sLMbOmADaZ4ZSIiCgVsjqkJlTD7BIoj9ltEmZWF2BGNcMp\nERFRqmV1SFW5aIpM4LBJmFlTiIaqAraSIiIiShOGVKJx2h9OZ1QXQJYYTomIiNIpq0OqrvN2P6Xf\n4PalhZhe5WM4JSIiypCsDqn/v727+W0jLwM4/njs8diOX2InTtIkbErbRUBRJTjBYTnwn3MBDiBx\nYIFlYRdtt6XdaJV0m5emCTGHpkulhXZf2s5j+/O5jCPn8Gg0kr/6zVuz2ah7BBbYar+Kmzuj2Flf\niaJwrAHA2zTXkWpVi9et0WjE5rgbt3ZGsb7arXscAFhacx2pzUKk8no0m0V8b6MfN7eHMei16x4H\nAJbeXEdqx6sm+Y6qdjNuXBvG9WtDb4cCgETmOlL73bLuEZhTg147bu2MYndjxYo8ACQ015E66IlU\nvpnpajdu7Y5iY7UbjYaboQAgqzmPVNcO8mpF0YjdaT9u7oxitOKYAYB5MNeRWraKGPWreHR0Vvco\nJNTrlLG32Y93NgfRreb6UAeApTP3v9xb465I5UtFoxFba73Y2xo4pQ8Ac2z+I3VtJT64e1j3GNTM\nqikALJa5/zVf7bejV7Xi5Oyi7lF4y6yaAsDimvtIbTQacWNnFH/66PO6R+EtsWoKAItvIX7h9zYH\n8cEnB3F+cVn3KLwhVk0BYLksRKSWrSKuXxvGh65NXTgrnTLesWoKAEtnYX71f7A7irsPj+LJU9em\nzrtOuxXb6yuxO12JybBT9zgAQA0WJlLLVjN+cmMSv//rZ3WPwrdQtorYXluJnY1+rI86UTidDwBL\nbWEiNSJid9qPf31+Evf2j+oeha+h2Sxia9KL3elKbIy70SyKukcCAJJYqEiNiPjpu+txdHruAf9J\nFY1GTMfd2J32Y2vSi7IlTAGAr1q4SG01i/j5jzfj13+87/rUJBqNRkyGVexO+7G9vhJV2ax7JAAg\nuYWL1IiIbtWKX9zejN+8/yDOnv677nGWUtFoxGTYiR/eWI9eqxG9zkIeagDAG7Kw5TDqV/Hene34\n7fsP4vj0vO5xlkJVNmNj3IvNSTc2x90oW82YTgexv/+47tEAgDmzsJEaEdHvlvHLO9vxuz8/iIPH\nrlF9E0b9KrbG3dic9GI8qDxkHwB4LRY6UiMiqnYz3ruzHX/55CD+8emjuJzN6h5prrWaRUxXu1er\npT0P2AcA3oilKIyiaMTt65O4NunFH/627/T/N7TSKWNz8uw0/vqo41FRAMAbtxSR+txk2Ilf/Wwn\nPr7/OD789DDOzt1U9b/0OmWsDTuxPurE2qgT/W5Z90gAwJJZqkiNiGgWRdzaHcXe1iD+fu9RfHT/\nUZxfXNY9Vq0GvXasDZ8F6dqw4058AKB2S1sjZauIH+2N493dUdzbP46PH3wRh0twc1XZKmK1X8Vk\nUMVk2InxoIq255YCAMksbaQ+12oWsbc1iL2tQRwencW9/eN4eHASXxw/rXu076wqm9HvlTHotWO1\n347xoBPDXukOfAAgvaWP1Bet9qtY7Vdx+/uTOD27iM8OTuPhwUnsH56mviSg027F4CpGn22fffZm\nJwBgXonU/6Nbtb5cYb2czeLw8VkcnZ7H8el5HD+5iOMnz7ZvQ1E0oiqb0S6b0W1frY52/xukZUuM\nAgCLRaR+Dc9f8TkZdr7y3Wi1F/+8exDHT87j5MlFnJxdxOXlLGazWcxmEbOImM1mcfn876vt5dWX\nrVZxFaBX29YLn8tmVGUzypZHPgEAy0WkfkftshnjQRXjQVX3KAAAC8MSHQAA6YhUAADSEakAAKQj\nUgEASEekAgCQjkgFACAdkQoAQDoiFQCAdEQqAADpiFQAANIRqQAApCNSAQBIR6QCAJCOSAUAIB2R\nCgBAOiIVAIB0RCoAAOmIVAAA0mnMZrNZ3UMAAMCLWq/6h/39x29jjrk1nQ7so5ewf17NPno5++fV\n7KOXm04HdY8AfAtO9wMAkI5IBQAgHZEKAEA6IhUAgHREKgAA6YhUAADSEakAAKQjUgEASEekAgCQ\njkgFACAdkQoAQDoiFQCAdEQqAADpiFQAANIRqQAApCNSAQBIR6QCAJCOSAUAIB2RCgBAOiIVAIB0\nRCoAAOmIVAAA0hGpAACkI1IBAEhHpAIAkI5IBQAgHZEKAEA6IhUAgHREKgAA6YhUAADSEakAAKQj\nUgEASEekAgCQjkgFACAdkQoAQDoiFQCAdEQqAADpiFQAANIRqQAApCNSAQBIR6QCAJCOSAUAIB2R\nCgBAOiIVAIB0RCoAAOmIVAAA0hGpAACkI1IBAEhHpAIAkI5IBQAgHZEKAEA6IhUAgHREKgAA6YhU\nAADSEakAAKQjUgEASEekAgCQjkgFACAdkQoAQDoiFQCAdEQqAADpiFQAANIRqQAApCNSAQBIR6QC\nAJCOSAUAIB2RCgBAOiIVAIB0RCoAAOmIVAAA0hGpAACk05jNZrO6hwAAgBdZSQUAIB2RCgBAOiIV\nAIB0RCoAAOmIVAAA0hGpAACk8x+tLWgK6zcqWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1115399e8>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pros ###\n",
    "- A principled method to trade complexity for bias\n",
    "- Optimization theory is applicable\n",
    "- Assesment of convergence\n",
    "- Scalability\n",
    "\n",
    "### Cons###\n",
    "- Biased estimate of the true posterior\n",
    "- Better for prediction than interpretation\n",
    "- Model-specific algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean field Variational Inference\n",
    "Assume the variational distribution factors independently as $q(\\theta_1, \\ldots, \\theta_n) = q(\\theta_1) \\cdots q(\\theta_n)$\n",
    "\n",
    "The variational approximation can be found by *coordinate ascent*\n",
    "$$\\begin{align*}\n",
    "q(\\theta_i)\n",
    "    & \\propto \\exp\\left(\\mathbb{E}_{q_{-i}}(\\log(\\mathcal{D}, \\boldsymbol{\\theta}))\\right) \\\\\n",
    "q_{-i}(\\boldsymbol{\\theta})\n",
    "    & = q(\\theta_1) \\cdots q(\\theta_{i - 1})\\ q(\\theta_{i + 1}) \\cdots q(\\theta_n)\n",
    "\\end{align*}$$\n",
    "\n",
    "#### Coordinate Ascent Cons\n",
    "- Calculations are tedious, even when possible\n",
    "- Convergence is slow when the number of parameters is large\n",
    "\n",
    "### Variational Inference in Python\n",
    "- Maximize ELBO using gradient ascent instead of coordinate ascent\n",
    "- Tensor libraries calculate ELBO gradients automatically\n",
    "- In Edward you can get BBVI\n",
    "- In PyMC3 we have ADVI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mathematical details\n",
    "\n",
    "* Monte Carlo estimate of the ELBO gradient\n",
    "    * For samples $\\tilde{\\theta}_1, \\ldots, \\tilde{\\theta}_K \\sim q(\\theta)$\n",
    "\n",
    "    $$\n",
    "    \\begin{align*}\n",
    "    \\nabla \\textrm{ELBO}\n",
    "        & = \\mathbb{E}_q \\left(\\nabla\\left(\\log p(\\mathcal{D}, \\theta) - \\log q(\\theta)\\right)\\right) \\\\\n",
    "        & \\approx \\frac{1}{K} \\sum_{i = 1}^K \\nabla\\left(\\log p(\\mathcal{D}, \\tilde{\\theta}_i) - \\log q(\\tilde{\\theta}_i)\\right)\n",
    "    \\end{align*}\n",
    "    $$\n",
    "\n",
    "* Minibatch estimate of joint distribution\n",
    "    * Sample data points $\\mathbf{x}_1, \\ldots, \\mathbf{x}_B$ from $\\mathcal{D}$\n",
    "\n",
    "    $$\\log p(\\mathcal{D}, \\theta) \\approx \\frac{N}{B} \\sum_{i = 1}^B \\log(\\mathbf{x}_i, \\theta)$$\n",
    "\n",
    "BBVI and ADVI arise from different ways of calculating $\\nabla \\left(\\log p(\\mathcal{D}, \\cdot) - \\log q(\\cdot)\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic Differentiation Variational Inference (ADVI)\n",
    "\n",
    "* Only applicable to differentiable probability models\n",
    "* Transform constrained parameters to be unconstrained\n",
    "* Approximate the posterior for unconstrained parameters with mean field Gaussian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LDA and AEVB\n",
    "Taku Yoshioka example\n",
    "\n",
    "#### Problem\n",
    "Let us consider the problem of probabilistic models with latent variables (topic models are a great example), for **Big DATA**\n",
    "\n",
    "* One solution is autoencoding variational Bayes (AEVB; Kingma and Welling 2014).\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- In AEVB, the encoder is used to infer variational parameters of approximate posterior on latent variables from given samples. By using tunable and flexible encoders such as multilayer perceptrons (MLPs), AEVB approximates complex variational posterior based on mean-field approximation, which does not utilize analytic representations of the true posterior.\n",
    "- Combining AEVB with ADVI (Kucukelbir et al., 2016), we can perform posterior inference on almost arbitrary probabilistic models involving continuous latent variables.\n",
    "\n",
    "It is worth highlighting that there's been a lot of development in the Deep Learning community on encoders. \n",
    "\n",
    "- In the LDA model, each document is assumed to be generated from a multinomial distribution, whose parameters are treated as latent variables. By using AEVB with an MLP as an encoder, we will fit the LDA model to the 20-newsgroups dataset.\n",
    "\n",
    "- In this example, extracted topics by AEVB seem to be qualitatively comparable to those with a standard LDA implementation, i.e., online VB implemented on scikit-learn. Unfortunately, the predictive accuracy of unseen words is less than the standard implementation of LDA, it might be due to the mean-field approximation. \n",
    "\n",
    "- One solution is the work on **normalized flows** which is currently WIP for PyMC3\n",
    "\n",
    "However, the combination of AEVB and ADVI allows us to quickly apply more complex probabilistic models than LDA to big data with the help of mini-batches. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset\n",
    "Here, we will use the 20-newsgroups dataset. This dataset can be obtained by using functions of scikit-learn. The below code is partially adopted from an example of scikit-learn (http://scikit-learn.org/stable/auto_examples/applications/topics_extraction_with_nmf_lda.html). We set the number of words in the vocabulary to 2000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading dataset...\n",
      "done in 1.989s.\n",
      "Extracting tf features for LDA...\n",
      "done in 2.143s.\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# The number of words in the vocabulary\n",
    "n_words = 2000\n",
    "\n",
    "print(\"Loading dataset...\")\n",
    "t0 = time()\n",
    "dataset = fetch_20newsgroups(shuffle=True, random_state=1,\n",
    "                             remove=('headers', 'footers', 'quotes'))\n",
    "data_samples = dataset.data\n",
    "print(\"done in %0.3fs.\" % (time() - t0))\n",
    "\n",
    "# Use tf (raw term count) features for LDA.\n",
    "print(\"Extracting tf features for LDA...\")\n",
    "tf_vectorizer = CountVectorizer(max_df=0.95, min_df=2, max_features=n_words,\n",
    "                                stop_words='english')\n",
    "\n",
    "t0 = time()\n",
    "tf = tf_vectorizer.fit_transform(data_samples)\n",
    "feature_names = tf_vectorizer.get_feature_names()\n",
    "print(\"done in %0.3fs.\" % (time() - t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each document is represented by 2000-dimensional term-frequency vector. Let’s check the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAFpCAYAAADUeIJKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFPWdN/BPT8/BAMMw6EA8EhcQPAABWTCCMcujqyKa\na2PWx/V6ErO6UUyMWfEC3ScvF+IGjxAP8A4mHjGn6IasBoKPgJCRSw6BQWa4hum5e7p7+vw9f/RM\nT1d3VXddv+7q7s/79UKnu6t+9buqvlW/ulxCCAEiIiKSpizfGSAiIip2DLZERESSMdgSERFJxmBL\nREQkGYMtERGRZAy2REREkpXLStjj8dqaXl3dUHR2+m1NM19YFucplnIALItTsSzOY3c56utrNH8r\nmCPb8nJ3vrNgG5bFeYqlHADL4lQsi/PkshwFE2yJiIgKFYMtERGRZAy2REREkjHYEhERScZgS0RE\nJBmDLRERkWQMtkRERJIx2BIREUmmK9iuXr0aZ511luLf9773Pdl5IyIiKgq6Hte4f/9+/OM//iMe\neuihxHdVVVXSMkVERFRMdAXbxsZGnHXWWaivr5edHyIioqKjaxj5wIEDGDt2rOy8kEWfHGxHjy+U\n72wQEVEKlxBCZJogFAph+vTpuOKKK7Bjxw4IIXDFFVfgzjvvRGVlpeZ8kUi0aB5WXQg+O9aNO5et\nw8m1Q/DS4svznR0iIkqSdRi5qakJkUgEQ4cOxfLly9Hc3IxHHnkEPp9PcQ43ld2vX6qvr7H9tX35\nIqMsjU0dAIC27r6c1lOxtEuxlANgWZyKZXEeu8uR6RV7WYPthAkTsGnTJtTV1QEAzj77bAghcPfd\nd+OBBx5Aebm0V+ISEREVBV3nbAcC7YDx48cjHA6jo6NDSqaIiIiKSdZg++c//xmzZ89GKDR44c3u\n3bsxYsQIXp1MRESkQ9ZgO3PmTAghsHjxYnz22WdYt24dHn30UXznO9+By+XKRR6JiIgKWtYTrnV1\ndXjhhRewZMkSfOMb38Dw4cNx7bXX4tZbb81F/oiIiAqerqubzj33XKxatUp2XoiIiIoSX0RAREQk\nGYMtERGRZAy2REREkjHYEhERScZgS0REJBmDLRERkWQMtkRERJIx2BIREUnGYEtERCQZgy0REZFk\nDLZERESSMdgSERFJxmBLREQkGYMtERGRZAy2REREkjHYEhERScZgS0REJBmDLRERkWQMtkRERJIx\n2BIREUnGYEtERCQZgy0REZFkDLZERESSMdgSERFJxmBLREQkGYMtERGRZAy2REREkjHYEhERScZg\nS0REJBmDLRERkWQMtkRERJIx2BIREUnGYEtERCQZgy0REZFkDLZERESSMdgSERFJxmBLREQkGYNt\nkRD5zgAREWlisCUiIpKMwbZIuPKdASIi0sRgS0REJBmDLRERkWQMtkRERJIx2BIREUnGYEtERCQZ\ngy0REZFkDLZERESSGQq2Dz74IG644QZZeSEiIipKuoPtxo0b8etf/1pmXoiIiIqSrmDr9/uxaNEi\nnH/++bLzQ0REVHR0BdvHH38cs2bNwqxZs2Tnh4gKTDgSQ0yYexVGLByyOTdEzpQ12G7duhV/+tOf\nsHDhwlzkh4gKiBACt/50Hf5zVYPhecNtHhz4t39F6+u/kpAzImcpz/RjKBTCAw88gPvvvx+1tbWG\nEq6rG4rycrelzKWqr6+xNb18srsstR6ftLSzKZZ2KZZyALkri+g/oj14rMfwMls/iQforvf+jEkL\nbtWcju3iTMVSllyVI2Owfeqpp3DGGWdg3rx5hhPu7PSbzpSa+voaeDxeW9PMFxll6e4OJP7OZT0V\nS7sUSzmA3JZFJA0fG12m19uXdV62izMVS1nsLkemwJ0x2L799tvweDyYPn06ACAcDiMajWL69OnY\nunWrbRkkIiIqZhmD7apVqxCJRBKfX375ZXzyySf46U9/Kj1jRERExSJjsD3ttNMUn0eMGIEhQ4bg\njDPOkJopIioM5q5BJio9fFwjERGRZBmPbFPdddddsvJBRERUtHhkS0REJBmDLRGZx5O2RLow2BIR\nEUnGYEtERCQZgy0REZFkDLZERESSMdgSkWmCV0gR6cJgS0REJBmDLRERkWQMtkRERJIx2BKRacLK\nKVue7qUSwmBLREQkGYMtEeWHK98ZIModBlsiIiLJGGyJiIgkY7AlIiKSjMGWiIhIMgZbIiIiyRhs\niYiIJGOwJSIikozBlohMs/QEKaISwmBLREQkGYMtERGRZAy2REREkjHYEpEFPGlLpAeDLRERkWQM\ntkRERJIx2BIREUnGYEtERCQZgy0RmcaHWhDpw2BLREQkGYMtERGRZAy2REREkjHYEpFpPGVLpA+D\nLRERkWQMtkWCRxhERM7FYEtERCQZg22RcOU7A0REpInBlojM4/kLIl0YbImIiCRjsCUiIpKMwZaI\niEgyBlsiMk3wpC2RLgy2REREkjHYEhERScZgS0REJBmDLRERkWS6gm1jYyNuvvlmTJ8+HXPnzsXz\nzz8vO19EVAAEr48i0iVrsA2Hw/jud7+LU045Bb///e+xePFiPP300/jjH/+Yi/wREREVvKzB9sSJ\nEzjvvPPw0EMP4YwzzsDcuXMxe/ZsbNmyJRf5IyIiKnhZg+3pp5+OJ554AkOGDIEQAg0NDdiyZQsu\nvPDCXOSPiIio4JUbmfjiiy9Ga2sr5s6di8svv1xWnopCuM2D1l+9ivp/vg6VY8bkOztSiVgMLS88\nh+Hnn4+aGTPznZ2S0Pbbt1A2bBg2TSyH2+XG7I5a9P5tCz53y7/mO2uIxWJ47497cNbkz+GMM0/S\nnpDne03bvuUwAr4wvvgP4/KdFcfbvP4zVFS5Mf2CL+Q1H4aC7dNPP43W1lY8/PDDWLJkCR588EHN\naevqhqK83G05g8nq62tsTU+m3c8uh2/HdrhjEUx55P+m/W53WUZ4fNLSzmZ4yAvvRxvh/Wgjxv3h\nNzldtp0KqX/te3c1AOCd60YDAMb+qhUAMO5/fxNAbc7K4u8LJ/4eWOahxjY07vWgca8Hi5ddrTmv\nqBmClpR51RRSu2RjV1k2vN8IALj6mqm2pGdGobRLw4YmAMBlV01S/T1X5TAUbKdMmQIA6Ovrw8KF\nC3HPPfegsrJSddrOTr/13CWpr6+Bx+O1NU2Zgr4+AEAoEEzLt4yydHcHEn/nsp7q62vQ0d6bl2Xb\nqdD6l5bOTh+GI3ftEAhGEn8PLLMrad3PlA9vb1/W6YqlXQA5ZclX3RRiu6jl1+5yZArcui6Qev/9\n9xXfjR8/HuFwGL29vRpzUa7x5fFERM6VNdg2NjZiwYIFaG9vT3y3a9cujBo1CqNGjZKaOSIiomKQ\nNdjOnDkT48ePx7333ovGxkasXbsWy5Ytw2233ZaL/BGRg/GhFkT6ZA22FRUVWLlyJdxuN6655hos\nXrwYN910E2688cZc5I+IiKjg6bpA6pRTTsGzzz4rOy9ERERFiS8iICIikozBlogssHDSlud7qYQw\n2BIREUnGYEtE+cGbw6mEMNgSERFJxmBLREQkGYMtEZnGa5yI9GGwJSIikozBloiISDIGWyIiIskY\nbIsEz50RETkXgy0Rmca3/hDpw2BbJPh8ACIi52KwJSIikozBloiISDIGWyIiIskYbKXJ7ZUj+b1O\nhVfJ0CBeNEWUjsGWiIhIMgZbaXJ7fXB+r0bmtdA0yMXuQJSGwZaIiEgyBlsiMk3wBC2RLgy2RERE\nkjHYEhERScZgS0REJBmDLRGZZumMLU/3UglhsCUiIpKMwZaI8oP341IJYbAtEhyRIyJyLgZbIiIi\nyRhsiwRH5CgvOKRCpAuDLRERkWQMtkRERJIx2BIREUnGYCsNXx5PxU+t5XW/m4DdhkoIgy0REZFk\nDLbS8OXxVJp0vzye3YZKCIMtERGRZAy2REREkjHYEpF5uq+GIiptDLZERESSMdgSERFJxmBLREQk\nGYMtEZlm6YwtT/dSCWGwJSIikozBlojygw+1oBLCYEtERCQZgy0REZFkuoJtc3MzbrvtNsycORMX\nX3wxli5dimAwKDtvRORwfKYFkT7l2SYIhUK47bbbcOaZZ+L1119He3s77r//fgDAvffeKz2DRERE\nhS7rke2OHTvQ3NyMJUuWYPz48Zg1axa+//3v4+23385F/oiIiApe1mA7btw4rFy5EsOGDUt853K5\n0NPTIzVjhY8vj6fSxKFlonRZg+2oUaMwe/bsxOdYLIZXX31V8V2uCCGwZ/tx9Hqdeb7Y3xfBuq1H\nEQxH850Vw6KRGHZtPYZgXzjny27p8GPT7hbV33p3bEfg4EFbl+cP+/HB0U0IR3NX1m2eT3C093jO\nlrfdswuHvcdUf+vzfoa+3iZsfe8D/PeKV0wvw+Nvx9a2babn19pHO+47gY9bd5hP12GChw/D+3GD\nlLT93Z8i5Lfer7p6O9Dw+xcR7vVaSicajeBvf3gRnW1HLecJAPqam9C7bastaeVb1nO2qZYsWYI9\ne/bgrbfeyjhdXd1QlJe7TWdMTafHj3X//SlGjqrGnQ9camvadlj2qwasazgCXyiKv6+MV21FhRv1\n9TVp06p9Z8VIj89S2h+8tx/r1+xDy5FuXPvtWYbmHTVqGJosLPvbS/8CAJg5+VR87qRhit/2/exx\nAMCcP/zGcLpaHv1/v8Tfjm4HKiP4xrnzEt/b3SYDItEInvvLLwAAb/7zM7akuU/j+7q6oRBCYOXO\nVzSX17B1FQBg098uBnAG2g/sxNkXGt95vuONhRAQcFXPgQjUJOqvt3twZzhTnYqaIWhRme72v9wD\nAJjz9Z9mTaMQfHjLIgDA382dbXtZ2g6+AQCYcdl/WUrn9yvvw7mbj+Ng83Fc9OOluuZRK8tfX38O\nI95ej/0f/Q3zn1tlKU8A8OEtDwGwZ/3Xqvtc9S/dwVYIgUceeQSvvfYannzySUyYMCHj9J2dfsuZ\nS1ZfX4PjR7sBAF0dAXg81vbAZDjQ3AUAOHikC1ND8aPbcDialtf6+hrb89/dHUj8bSbtY0fieT92\nuMvQ/PX1NejoGGxrK+U6erwb7lhM9Tc766uxLb5rcMhzLJGujDYZEEo6gpbdbzs7/Rg2dvCznuW1\nHm3BSSbyJfoPTV3lYYikZXV36esP3t6+jNO1tHbhzNOHOXJdN0PEYtLKYjldTycAINJ8XFdaWutL\n95GjOAnAyFa/rWW1Iy21NOxe7zMFbl3BNhaL4YEHHsDbb7+Nxx9/HJde6ryjSiKKEzyHXjKEAFx8\nEldB0BVsly5dirfffhvLly/H3LlzZeeJiIioqGQNttu2bcMrr7yCu+++G5MnT4bH40n8Vl9fLzVz\nRGSC0QNbywfCPJImyiZrsF2zZg0AYNmyZVi2bJnit127dqG83PA1ViWilDZApVRWyoa3/hClyxop\nFy5ciIULF+YiL0RUkHjSkCgbvohAmlLaAJVSWZ3P6AVSdh+I8oIdspsoguESBlsiIiLJGGyJikzh\nHwMQpeCRLREREWXDYEskXY73yovgKICo2DDYEhEVKO5XFQ4GW1sVcM/nWlswiuHKTAAFvboUG5fT\n26II+jyDLZFkud5MGL71p/C3Y4WBFV3SGGxtlXyDYYG9PN7SzZHciNAg3TGF9+PawJ5KFGwL6Rhs\niSTL9bAvd32ciQe2FhRB5THYSpPbXUXumJYQh294+AQponQMtkTS8dYfAjjmUNoYbMk6btyJiDJi\nsCUii7izlS9FcxtYFsVQTgbbIlH4XbF45f7WH3KkIggYZB6DrTQltGKVUFFJjfKKKN0xhf3GMlZh\n4WCwLRK8ANTJcn3rDx9q4UgS6rl0ms5cSZ00/MxgKw3DH0nioA2IGt23/nAVoRLCYEuO5aS9Uity\nXowiqbfiw3YxrQiqjsGWbFAEawJRAeJ+VeFgsCUqNJm2sHnZ+nKLny8lU/NFsFfBYEvOVQQrWByf\njUzFc1qEzGGwtZXQ+LsAWNgQcCNS6kze+kM2sOcqM8e/z7YIMNgSSeb499lKyodVRsvheLz1J+ec\ntOPHYGsrl8bfBYCvaikODti46O5KDsgrxfF9tvIx2ErD83SWOWm31IKcH6EZXZzN2SuSZpPA/oop\nmaougk7FYEtE+cGjKZLOOUGawVYavjzesiLYmwUg4cgxU4IiD+c6lcuz64xE0V14J6E8RVZDGRR+\nSRlsiYiIJGOwJeuK7QiEqFCUyKpndhPjpE0Tgy2RZLke1i26W2aKhJM2/JR7DLa2yt/aZHnJ3BIU\nEKc9rrEorxgoCHa1tuMfalEE2ycGW2kKv3PoJqmoRXeBTK4YrDa7a9muZiu+I3QJF0gVWxUVMQZb\nWxXwHj4falFSHLGRdkIeCpxdVciHWsjHYCsNe691xbE1tv0ILcsosuHlWT8Hofhk360/9qTjGEVX\noFwq/LpjsCWi/OD+qGWFH4JKB4Mt2YCrfKEpjIOsgshknpXIHovZruCgLsRgWywc1KlsUyRlynVg\n07e4EtlIO4mUt/4UyUpSAhhsiQpNlsc1FoviKQlZVhhDMRkx2EpT+J1DN1krQhGsYHG5PrTNb70V\nTbPZjLeylTYGW1sV8MrEDUFJsbe1OSSdLyXzUAuTnDTMzmArDTdAFJfTFV7ofOuPUP3T7EIVn+y7\nZds5G0p7FFt5yAgGW1vlL8Ba3qDzoRYlxgHtzdhjGR9qUTgYbMkyjkCTDOxXOpRKHZl+7Y+92bCC\nwZaci1tbU/TUmqJq81XNpXY0JeXWHyoUDLZEhaZEdkKcdHEL5Vcx9AUGW7KBrBWh8FcwIA+3fBhd\nnt2Pbi6OZrOflIopteEBY5zUFRlsbeWkpjWIW8gSY+dG2mRa7HKO4fhbf5yePx0MBdtQKISrrroK\nGzZskJUfMqsIOmMqxn8tmV/7o+ucrWIWqxXNkQ19+D7bUqY72AaDQfzwhz/E/v37ZeaniJTQWmDT\nGs8NRw45oa45Akp6FcHGQVewPXDgAL71rW+hublZdn6IiAA4Y3/A6VhHWTiognQF282bN+OCCy7A\nG2+8ITs/GSUPd7V298IX9iMmBI61eNHR0wchBIQQ6OoNAgCCoSj8fZG0dMKRGHoD4YzL6gtF8Glb\nE2Iihi5vAMfa2xS/R6Mx+L1+dLYdTcpXfFc9JgQAF3qqy9EHdyLvXY3NCPt9OHb8EGLhMHzeIDp6\n+tDhDaKzsxdRnw+RSBQBXxCRri4IIdDW44M31Du43FgMPb6QVg3BVeVHT7cfzce6IYRAJBpD+7FW\nHD58CJFwGF5/CL6+9LKHo7H4/0PRtN+C4Sj8KvNEursholEE+qKIJR2mhENRBJPqvTsUhjccSdRN\nV28Qnd4ghBCIxWJARR8AwNcbRHdS2VpaTiDorgYA+H0hxGIiMV000odYNIiAP4TWjg4E+nyIRQIA\ngPYTXrR3BeDt7ku0Tbu3G4FAcKCa4A5XojvUg86+LsREf9kj0bR+ISIR9LYcQSDSp1HnynravbcZ\nu/acQJc3gEhXJ3y9wUT6A8vQQwiRmF9tmDdY4ULIDVSHhsAVE4jBhZB7CPr6Iuj1BlU3MkII9LUd\nwdH24egJuRPfhyJRHOlsRzjUg47eThw7eBSBQEhRF9GwD8ePHsKJE13pCbvjbevzBtN+ikZjCPgH\n27StuwvBcPxzIJi+bibEgAOfHUA4GkFLhy+xbnd6g4jFBDq9ATS3t6Orv1+GPR6IaFS1vqLRGBpb\nD8PT3Qtvd58iPwAQDQTQcfQgorFovN7DXvj8IbS19SLYF4a/Nwi/P4jdB0+gtbMd3cfb0OMLoqfH\nn1g3A/4QotEYWvYdQSQYrwcRi8DX14MDo0fDV+lGa7cXnd4+BMPhxHy+3iA6Ovw42uZDtLcXsfBg\n3g57jiAY6U9LRHG07RiC4czbrWBfWHUdFtEoulraEInG0NnfTjEhcLTLj95ACMIV71Nht3K+WCyG\nE0cPobXdh95wBH09XnhOdCPgD8XX/9hg3+72BRGJDH6OhcNoPd6DLn8PegK92Nd4HO2dvejubEcs\nFoMQItFeka5ORGOxxPovRAzRcC8EkNgGxL8fnMfnDSIY8OPI7l3o7vUh6velldsb9A/WTUd7oh81\nNh2D35feX2Uq1zPRddddJzsfuhzc6wEAdEHg4U1L4aoIYfqJ/41wUzcOI4YLvzQWwXAM725qwl3f\nmorH39wOAHjx3v+lSOe+lRvR0RPE8wvnokzjyUkLVryDykkfYtyIM1Dxl4moiFXgq3ecjVOHfw4A\n8JuXGzDu49cxItiOv/7gq/ja5K8n5t3b7UdnZRQvfX0UajorsRTA2tc24NPmMGYc+TlG9rXiszMu\nwo6K8diTtFW87dBvsOvc69DnD+OSAy9jw7w7sH5/DyrP2YSfXfUjVLor8V+vbcO+w114YsFFGDGs\nUpHn8lMbUXH6AbywsgpDohU469Lx+KDhEJo64yvoP09bjd9+MhlDq8rx+IKLEvN19Qbx0Z5W1MOl\nCJID7nh8PaIxoajHcGcnPvv3u3B86vlY7TsPtaddgb8/+t8AgOcf+wAA8G/3/gM2tXbhj03xdvvB\n5DPw7rpGrN9+HADwtS+NxZ7gRlRP34bgpzPw7ls7cRwCD94+B71eD159rREjx/4zLmj6HV5ZvgGn\n/10dzpt5Ot799U6MH9eMCeOa8Kf3vgQAuOIfP4C7TMA18g6sfmNHIp8Xzh2PCeefhDef2prI08k7\nJuG0nlrsmf4eHuz4T1zwuRm4e/QtuOeZjej2hRTlbHnlRXg3bsDrl9dh8TWPq/aVAU/89zaM2O1F\nuzuMg1E3/s/R/8Enp12OCdNOBvqb6p5nN+LxOy7KmA4AdLy7Gp+9uxZbPv8VnD3lc5g7/+zEb0IA\nz15TDwCorbkBc9b9Ch+f9vforh6Dob/ZC//qEzi1fgqOjd2Jjce24MJTZwIA2vf8Fut/2Y0TNecr\nlvXxzk78afh/4eYRQ9HZ5MOGvfMQxX58jBhevPd/IRr24cDmn8D1UjNO1A1B+Pb/RH394EasauJW\nDN98BX7x1EZcfPlE1NYNbhzferkBHR4fbvnhlxBFFL9+Zhsi1X1Y8P0r8Oqf9+FqjfLPbp6BHX8L\nYuva38DdOxQnfenzqK2owBt/OaCYrurkIbh/5nAElj+GwKx52NAxBudf+AVc8OVxiWmeffct7B72\nN4QOTsHn205HPVz4t3v/IfH77h9+D89eU4/TGkfh9olz0H18LVav+RJcKePcO6b+Bf9nzVGM7I3i\n0fH/gtrJH6FvWA9+MvthrPrZ5sR0Q12f4KaFV+DorifxdMfZ6L4UAE4CGv4Lgc1XYMS0zQhXduCe\n8T/Cu6/tBgA0IoZ/PfALuIfXYPwTy3Gw9RCWffI0qgJ1eGz+fTi08wX8+N0zMWT4Vjx9x/zEsg42\nD8PJE72Jzy8+8SEAKMoHAIeXPYq+fZ/iwb+7Bv7yanxn/jlo2N+Gbfs8qJk4EnMx2KeeSppvy2s/\nx8fNkwHXIRy+5DQMCfhw6vrjCJdX45IDL2PopMk4/a4fIRCM4KfLN2BauAdj+ufd+OOfY8fQ6YhU\nBOESLrgjlWh0xTBelGHajDII96nYvvkIvvz5bpSv/R02Tb0a63x1eOyOOQi3/gGB7r3YPeZLaKkZ\njzHtfoyurcWG9xux429HMOHc0di/uxXDwifgqxiDiz5bjqpoHyY+/7Ki3C+8tR3D+/9uuudujLry\nKvz1tJFoXRsPfYuXafVA++kKtmbU1Q1Febk7+4QGHG2O71WHIeCqiO8BtTR14SS4cDJc+HDXCfT2\n77U2tgx2wPr6GkU6HT3xPZqTThqOcrfGwf3QTgDAwZ4mTI5Nis8nPJhaPwEA0O7xYVqwHQCw7eBG\nfHfujSgvH0yrvSa+h+etO4H6+hp82hwPeJ3Vn8PIvlYcdZ2CwRzGtVWORJ9/cM91/f4eAEDMOwpD\na8tRV12DfYfjdRBxuRTlGnGiF+76owCAIdEKAMC+w92JQAsAn7aehHAkhu5ISDFvS49yDy+1vqL9\nR5TJ33cd+wwAcGJvE/D589BdPSZtmvr6Guw6cCzx2et2JQItAHz4SQt6x30CuAF3bRvQPRq1cCEo\ngOPt7RhZFu8/3iEnAQCOHOrEKafVAgCaD5+CsWccTaTlLovnse344CgAEN9BmzBnhCJPVT3xNCqD\n1QhUhPBRSwOAWxJ71cll2LcxfjFgrTeaVi+p3J3xI4rh0XIAAieGng4A2L+tDZgVn6a7N5Q1HQBo\n3vj/0D1kNABg784WfOvmmYnfIn7lehWq+Ty6o/H6D7uHAABGeT6PY2N34pOuXfjK1PjOw6FDDThR\nc1nasmIATutfV+vOGAbsBQaWUF9fA19XJ3zeMIYDGNPZh/YTPpxy8glFGqP6A1PTgXZ86dIJie87\nPPGjjWFDqxBEfHSgPDAkrQ5SP/e2nQwAcPcOjafb2I6ASF9Xg2196D3eDTeA5oMdwMgx2LP9OK76\n5tTENAexL55WXQs62k5DPZTrjn9IPN2jkQ4EOrcBQFqgBYCyqj6M7I23cUUsir5h8fVTVCqPNv0i\nXr7miA++YHfaVjZc2QEAaG3tTnxX17+8aK8X9fU1WPvp4Xj5qjtRX1+Dbd74et/XW63Ie+fxcmBi\nev2lft6371MAwMhwL/zl1fjkUCe27YvvBJeVl2nOW97wCTD6vMTnvuphCJcP7kz5d32C+voaHPP0\nYnRKfXnCw+JphKsGy9n/XMjGfQH4/fFtQ/PeFowDUHFoH1B/AULChb7uvQCAlprx8eX44uvfzoYj\nAIDP9sVHGn0V8X7fV1GDqmhfWrmHn1AePHT/dS0aZ12IGpyuWVeySAu2nZ3+7BMZoKdCYtFYYggp\nkBS0PJ7UsDb4vWawVeH19qmnJeJpJQ+hJP/d2tqjexla2tt7EakazGtnlx+eoRWJz93dgbR5QiFl\nR0seXEsuR3eXsq0y1dcAv8ryUqfxeLwIhweHtLxe5TyxaAxqurv8CCiGcwc3fIEsw//+gHKIMBKJ\norNzcHhJq2zJWlt74FIZ8dAzr4LGBUB60omm1E3yPNGAet2rCYWixvOdstygTzk85/eH4NMYgguH\no+juTl+iBVXUAAAgAElEQVT3Ozp6EcTgPKnrhJ48Jq9TyYJ9YQxN+hyLiazpaf0+sGNpRFdXenn1\nlCcQUD8d5PF40Zf0W2paamnrmSZZMMMQvnLe7PXh8XjR3mFsey9E0oMqUob9u7r8GJIyfW9vUJEv\nrSvos5V7YBjZyDxGZIpTvM+WqNAUwZWZRKWGwdZWTnjgrDnOvgsjt3Vp29IKqwuoMnMPrp5Z4pcQ\n2sV4SnyRu5KZh1okzyJ/+1H47VUEwdZ8Ixh+qp3JFdTMbDK6lnB4SHWMwl+vbSOtKljHSgz+Ra8I\ngq0Vdndwl8bfDlmRMj14KHe5cLxieOi5bSRVhc5X3EtTLLHNrlo09z7b3O28F0N7Gb5A6tNPP5WR\nj7ww3IBmX6loaj5jHVnv0J2TOSV/zl+x5WUwteyyQmL+69jegexSlNMazH+Hsaykj2yd3HxOzhsV\nF1eG3lYE2zhVxVouZ+HOTLKSDrY5i2gOWbELawOTv8wWVj0VJtZxMchlMC38DlPSwdboEJnZDYS5\noTgZHblQ9zRznW+bVmxZ2bb7/bPJf6cOIxfEFVL5uAKx8Df+dtJbG6ZXiSKo7qIKtvnfW+atP2lM\nn+c2e+W3kfnUS53/fuR8mepZ66e8325TqPuaGVkvlBAmb/3ReNSt/gSUH10m+lQhKapga1TOGtAh\nHcUh2SBHS92AFucFUlYv/FLca2BDYcyHLeeu1c7NWX4UVbC1uqNlr8HMmOl0xvca1V7zYmLBdjPZ\nJi49t06pfK32mEVt6uk6/qEWOYxUWosyVs8DaQ0mlpeHWtg68OSElcsedmw37bhzKPM2r/Dru6iC\nrVFyb/0Rqn/mU76PJozIa1YLqJ7sllp0k4P5tkwik9XF238XfQl3Oh2KoXZKOtjKfaiFvcux535H\n7T3HQuzM0g4cC7I2zHNl2jG0sSpiiiPbPNdysTSxTcMDTn+oRTG0VxEEWwuPazQ6vYEZkh+NmIuH\nWuhRBP01hZyVvZBGAOynrFM7Q2Jyvea/ji2esxUZdlBKRPZip6+fjjrTl2NFEGzNy/8Kry3tAMOO\nzCo2dg4uPFLLr7GKqhXBSeXKy5bF1BUC2r/YVJ1C2JhYeurG57AxK/ntcc7o79JzYeudYvkJ+SUd\nbHPFKdt/p9+YJDQ/5JZT2isfbNmn00gjltSohV7His21DWUpzCO+HD4bWWclO/kgoqiCrYPrWTrH\nlt10vszeZ2tkannn2AtLrl5Rx/tsncj0tsNqfRq4zzbfXccORRVsjTK8wcg6edKeu+VdXxnnzpLS\nTO3oOubO6V6j6Qta7bgYzS6Stu4Z2sHuJdo3jCwgYs4ZW7FeriLY+icx9z7bXL71R2cGHXvUUQzB\n1kJ7O7dZkH4PmqkDCeVMyk2difNcJrIgVZ7OvRS3HD3UQkqqucuBcl/a3uspVH/Ws/MpO9AYTd9h\nq2e++1zhB9skdjwHwsbJB+ezZa/R7gukrM2vm9kVzuTIgLE+oHFsrPK1mWI4+fyRFr13/ph6qEXy\n3yLPD7WwUaHfKpbLhwG5ND8g80Mt9ObRwetcUQVboww3i+YMaj/Y24Pt6EJW08jpRsX0MLIkzl2H\n7ZVSTjuvRnbSzofizh0Tq2ryvcn2FCtzInoWkc/qVWtbW7OjN7Hke7kd1N+AEgm2mpVuU2O4VP5S\n7sU7IVQCivylblT15CBvfTe341G2raSSDhliImZzirloWJFyn62td/DalpKZRZotibG5NF6aYSFF\nxZy6LtpI/SLbTPkdR3baiENpBFuD39uSkOV+lnKBlO2nhZzVEVPpyp3qFkLOBVKZXrCunY6kh25I\nSXVAar8zdQ4k/Suh/CHfvS/1Oc1GKY5sYyZLY2A2fedszWXDrNTTAnIXpm8BRvtrLo9+SyLY5ofD\nrg4oOPkcE8vfovMjNwV22Kge2ahgmzaHGS+qYKu5MssdRc66IDsODuweEknNk75bf4wuQ5juzGbb\nxo77bFWTcNLWRGLU0nt6If2IIHueFH3Y1gukdNJaoOXXslpvj3zumpt9n63iNIkNG7nM77PVfdI2\n+zx8gpQ8ssfuDWyyLbBheFS5rctLHuyksWujc7pMCafP4aRxinxcGGbsueAqbSCUw6357klW1wU7\n+oOR7ZK+0yqms2JqVqHxdybS16OMGUlfei77YREEWx171JpHtkbH9/VnQXG+zsxeX8rFNaY6RaaZ\nzOyIGj6yNTi9sckNzGCwnY3mQ3I66QnbnXIu3gaV8p4fWy+QyofB3MdyMj6ePUxJr9FMF0jJXrSJ\nh1po1Yft90jrVATBNv+yrQZW9xrjX2RORc8yRNLwSW4eauGMzaktw3z5OvxRIaVWB45gtZ6mYnWj\npHN4OmeSN8qW20l+aXRVv+RsOGl0Rw8HVJlCUQVbrTst7LvzR/5Vh4NsOLLNIBcPtRACUtdQtat9\nVfuAZr7tO8eeU7bf+pOcdsaPCRkfaqE60qPc6bH3oRY6aQ4dG8+J4uXxJjuMU040uVx6581wTlU9\n5SwLTknDjlvlMrzHMd/XYhRVsNWmdXbPpnFOjY1LttksLdz69Qgm5lc7H+r8/V3tdtafd1NPkMrD\nrT+ydsospyscthNjeWWw4/yz/q2E7HO2li8Sc0zjGmuXXJ7MKIlga/QqZaPUHmoBq+dsU6SmYOO2\nAoCcq5EtZdL0GJDKDoHm0IbG12rTO2VbAhkbiGwbKGG5/CIl2tr7UAsz+UnKh6knSCWlZcP6nTUL\nGpkUivMb2fOhdb+4mYdamNmZlLVrPtieat/qmDEHSiLY5oXzD/ioYGlvIZzc7eS9O95Be0IlyikP\ntTCcrJRU1RVVsDV4AGO4orVP/an9Yq0ZU/ca0/ae7RjXs0WmK1nNHxHpemKOzvm0r0rUOFowcsRb\nSHRlX6XfqTRjprrQGnDI+/XHysNR5f8tMN0vDMyma6BH4rqmSmX1MXScnTq6Zqkt+jNjy/C+HEUV\nbDUZDMJG01F/IrL2c4jtkLrh0rPCK69GNrFMtZkyHUpJ7+36juOMXk+kFhRcGt9nTkeOzG1trkFc\nGSexPoycmrStF0hZXMHMzK0YjjUda3MfDjQfO2oqK+l3N9hy/7HaSw10r8PKUxU6FqY3YcsKPtjq\nuYBNs1PLfBKP1W6XeqVeWnDVk4mUeTR+0l8NKiuBoal1LsVSs6jl0WCC9g9UFADto3w7ip7+UIv8\nVmiij9kQ9WM5KIvMI1sApupAa3ti82J0EWl/pH1I+lrugZCWgg+2ehgdXtZOSP1r1Xswk740N8wk\noVtarAetwVijyzO3HPMTGn6Mp8p3LphpR0lXIxsdwtX1o/ok8aBoMOJqHR0bzIv+5VlIzOSFWvYc\n2Rr4TXNil56JVKc2LC15Wf1bbdHZyqZnGDm/VzSURLDVJHGvxmrS2c6RmQngQjG0rTy0NXo1sr7l\nWzmJpGNYXPe3WmkZedCmmfp2Fl35Sb0sVQz8z+JQrVD2mXzXjaL/Wn42sv33Paet77YFCrt6d8r8\nqgmYy7N9fUPPqTXbFpZVEQRb8xVqtJ6NTW915ZD8UAurafdXasa9c5Ojt/pn03nO1obKc4n8D30O\nMJ0PE2P+ApCyRXJGTZoddVKOZpmvnsEZ03qyzkSTLyrSM4upJ6ENzCt9hn5mbr1zqRzZ6tu71Jsr\nywo+2CZXlfb5W429OZse3putA9uzrUo9srWYgolRMPX6MvaUGF3LsTCh2pONYooVMMMGbnAia/mS\nTOZV0enDyOo/ZHyClGbatkQolYT1nk9Q/9JcTpKP0k3sxaT9ovyUeh5Y5hGY/qbMVBYTGdR4gpSJ\nrUzSjrqxwwge2RYsJ9/lSHo5Jajmjsxxk6RvS69ipZCxs2XmwR7Kl63Ylxf1ZRW+gg+2arfOpcpd\nQ2ktyXoO0q5Gtpye8cSU8+jcazSTUQsHKuobImPnbDVP+xrcyOXu+FPibAONmDISYvR9tkIIeeds\nTVwKO3g1MkztI+vZ7lihO01FmbLPpHbrj+4LslNntRprU0dLMp6aMr5RcFqALshgq/UwAi3aDymw\nITPQ6qjWLi9Pf/tK6gTGE1XMotjwmYq28f9lHIMyO0AnbN17V9xyomOUSWsYyykrr/l86Jgz5RBH\nY+TV+vLtjbbm57HlsZHmhpEz/5S6c21gx9AEy8kY398xlFb2zYHa1ch6Dghyt1YXZLC1i+Fq1jxw\nVdlbNNF9FStUtvts9SWYQuOhFjr3bNVirZSuqjtRva/4MXoex64Ritzf+mNepo25yP7M2SyPAo/f\nr6s8stV8wIJBZi76UebRzDN+BxMw+z5b5TnsDL9lSSXxl7mNggXWXtmpxVpaGerUjrcKWVASwVb7\nQfS52auR8wQpa/MI7T6pq770PbFKf94U89k8YfKNGWY3si5hvM6tvydVI10ZM2rsPAnlf6wtP6kh\n7N1hsJCWMPciAt3Jm55RZPqYYHSHxc77bLW2JxqTm16OnsRE2h8ZgjYfaiGPXccpWtOrHl+5kqtW\n35KUe16pw3kZejrUO0367QQayzUzipz4ZP+WKn5PptH8ZPhWIzGtnGvvazhkIFlmNtL6lYgP61tN\nVgh5T1qydGRrLk/JQc7s1ciZZExRxsGDwXO/aTOZ+l2d+tKz3cs8sCyDdcBhZJvJjraq32c4dDSx\nqNQ+YW6YOjk97aFVXeuyjoBkPjjZuwJo5lurDCrfad+UkHvC7MvjbRqdM51O8iph64Gt+XEeO7qo\nHf08NSxlOmer7M/Gnneu9bB/YXFIWG8dSHmVZ1LCyhG77DvZuVyjSyLY2hVrtai/z9b4cpTTpa9+\nphLNmGZ/Umlp6TmsFOazoS9pHRPqfKiF7qOQjInYtm/mPCLpuCD1rT/xb21YRKaPuTeQAVuGkY2U\nRt9wvKngZWU03VF3LdoTbWU+5M6Mgg+2yXWV5/PfKRyVmayk5dZkwnltS6cMFzuY5Yda5EORNKvd\n9Wh2XbOcC42HWqguyzlnIEwr+GCbbmDAL/uYlfFXlWldsany1CJTYxWaZxHVL1xJn0xfckgdbtH+\nTfm99gCrarpZy609EqB5QUiWe/vS8ujKNPxtdCQic4Fs3gRmWE7KUH2mirY4Zqb31KaeJ/xo3Xqm\niytDUFDcNKv3TWDpfxnLzuB8sYxJaP+Y6SYr/W/5MlaSgWVmqiPl/dBJ43Zp25f0oe3UUuiSGnQt\nnXbLdr7DpVwgz9lmpnUZkWLIK0Mw0Udlrqy3hKg3tN49Ue1ukn6Xp64nJ6YO3WmdzzB1ykuoL9PQ\neKHGqpkhfhheN0TqMLK+eVK5dCw7fftl5fA8U8hP3aHItOW0uNj+YJ5ahxnfPKS+F6Qc4jOYrUx9\nYjA1l/byB39WzmLi9ED64jMOCGeYzWAdWsqHMj96d4r1Xu2sbAGoftKcaWBTkuHUlJlnAajPY0eD\nm1OQwVZB4yY7PZ3JaIfWsw4nprV6Fj5lA5rtvKpqx0q9fUDjQ7YrndWTU1+9Mh296L/qWWhnIvvM\n6WJZfteRRDzzxvIkb53OtPdu9iEj6peAqddFhmQyLFmr/1lmKj8DO+XC+hOkMmUgw3pr5IBY+60/\nxgqv7376bAcVaj/Z2KAaMTIzc1cj5zLuFn6w1SDl/JChPmjmqEbf0CFgtm9rP9Qi07LUvlcfNkqd\n3qSMRzFJk6kM36sdLWg9nNylcVWIWhqOeoKU2Yxozqd9y1l8HuXRgKl1S2WEwa6HWpiO/ubnVsxl\nqD10bhYy7QBrHRELHdPoW7b63+nTGbsSGrAQ7PUuQMfdEHZdVGaUrmAbCoWwaNEizJw5E3PmzMFz\nzz0nO1+WKQ+wsgcPK7I9x0j3EZ2B+Gx1lCnjCqW10uopVIaEda8vSN0wG6kYtWCb8Wf9STsl3Oo9\nR2vHomycS8h6qIWVpNT6i470dDzU0tBv2W/9Sf6QvMOYIXkV6js4qSuqjnVb10+5u/JKpP2hRmXn\nPIfrdLmeiR599FFs27YNL730ElpaWnDPPffg1FNPxfz582XnzzxdR0Za32sFZ/0LM/qqp7jkPcXU\n4VmTkUtLyiXyWqPe2qUY+KT98A3TB2A6j2xV51VNz+DQksrkLqGRuMG82CHDZtjIjAkuiMGNdtp2\nt/+hFhYLM5COHMbTHRyZMXnrT9IiM18gpS+NtJ/03/+mOo8Qxq4y1l7PMyWS5WpFk4xufZPzous2\nPyEn39lkPbL1+/148803cd9992Hy5Mm49NJLccstt+DVV1/NRf5Ms3I0ozm51rCE2tcqV+qZXq7K\nr5YfapFpyRoTCq1pNOc1GxYsvIhA9UhFZPpZf9KG67wwbv+Kl8rAYzts2EDZ+1ALCzMJwOqzkY0V\nRt/OaOpv1tb2QVoPtbA6tqq2bljZ2dbznYLKKVt9p6L05sq6rMF27969CIVCmDFjRuK7GTNmYOfO\nnYhGo1IzV3gKY+NKRJTMrodaOOREiwG5y7FLZDmEWLNmDRYvXoyPPvoo8V1jYyOuvPJKfPDBBxg9\nerTqfB6P17ZMPv/Yc0DgVITdwwAAZbEw3CIMAAi7hyamq4j6bVtmsoFllEf74Eq6tNVlbefY0LLM\npOGOBVEmCneHKFGOaAhliCi+A4DyaBARdxWAwbaPucoRLatUfCfgQsRd3f9dAOH+v/XUj5E2trMv\nugQQLatAtKxCV3qJZSeNH6bNI4Bw+VCkckeDKEO0f5IyRNxD0udPOtWgtmFOLrs7Fkprg3gSbmV7\nZUgzOb0BmnXQn06o3K76dyX6SFksjFhSGyT6g0odRF0VimmziZRVQbjcJvI8mL/k7YSy7QJQnBM2\nva1SX5bm1MnD667B/GSjVfbBbVkIZSKSto4Nbi/j63KmflQRCcAFgZhrsB9WR5tx8wM36sqjHvX1\nNZq/ZT1nGwgEUFlZqfhu4HMoFNKcr65uKMrL3Zq/G+FyAy4RSVRuRWywI1VGvAiV16A8GkCZiNiy\nvEHK/RC3UJbX3X+yJlrm0thBUj+vqUfqsvQbzEh5LGgyjVwTUK0jEQNcZaiKBuJ1DKAq0oNg+Qi4\nY0G4RRARxFeagbYvE5HEhj65P0TQv/EUYaD/73j9ZG4jlxAoE5naeJCdfdElBMojYfgra+ESUZX0\n1PNdFfUiWD4CLhHpL6vyd3csnAjgA8qFsp9EEN9Api5T0d8VeQAqoj6E3cNQEfWhTERV2wCIpLWX\nOxYfgIyVpeYziGhZFdyxPkTLhqAy4s149bI7JlAW7kFfxYj+nSi1+tJzW9TANP07Z7EAgv31VSbC\ncPfHmngdKOcpExEEyypU2is534PzVUYjCJaPiP+tWr5MfbN/hzFlOzHYdmHF99m3VZmWp74sLe5Y\n/Dy4C65EfgYkryPC5Ua0rBJVkR7NtFwiBuEqQ3msDwBQFgshpjJPZTT+ezSlH5XFIoiVxcNcRSyE\nmAtwJfVDuGIZA6SdsgbbqqqqtKA68Lm6ulpzvs5O+44yv/P976K+vsbWo+V8Ylmcp1jKAbAsTsWy\nOI/d5cgUuLOesx0zZgx6enoUAdfj8aCyshK1tbX25JCIiKiIZQ2255xzDioqKrB169bEdw0NDZg0\naRLKy3XdOURERFTSsgbb6upqfO1rX8N//Md/YMeOHXj//ffx4osv4sYb7TupTEREVMx0HZred999\nePjhh3HTTTdh2LBhuP3223HllVfKzhsREVFR0BVsq6ur8ZOf/AQ/+clPZOeHiIio6BTtiwiIiIic\ngsGWiIhIMgZbIiIiyRhsiYiIJGOwJSIikozBloiISDIGWyIiIskYbImIiCTL+j5bIiIisoZHtkRE\nRJIx2BIREUnGYEtERCQZgy0REZFkDLZERESSMdgSERFJ5vhgGwqFsGjRIsycORNz5szBc889l+8s\naWpubsZtt92GmTNn4uKLL8bSpUsRDAYBAIsXL8ZZZ52l+Pfyyy8n5t20aROuvvpqTJ06FTfccAOa\nmpryVIq41atXp+X3e9/7HgDg6NGj+Pa3v41p06Zh3rx5+Otf/6qY1yll+e1vf5tWhoF/x44dK5g2\nCYVCuOqqq7Bhw4bEd1bbYNWqVbj44osxffp03HffffD7/Xkry8aNG/FP//RPmD59Oi6//HL8+te/\nVsxz+eWXp7XTnj17HFkWq33KKWW59957VdebSy65JDGPk9ol07bXMeuKcLgf//jH4qqrrhI7d+4U\n//M//yOmT58uVq9ene9spQkGg2LevHliwYIF4sCBA+Kjjz4Sl1xyiViyZIkQQohrr71WPP/886K1\ntTXxz+/3CyGEOHbsmJg2bZpYuXKl2L9/v/jBD34grrzyShGNRvNWnscee0zcfvvtivx2d3eLWCwm\nvvKVr4i77rpL7N+/X6xYsUKcd955orm52XFlCQQCivy3tLSIr33ta2LBggVCiMJok76+PnH77beL\niRMnig8//FAIISy3wZo1a8T5558v3nvvPbFjxw4xf/58sWjRoryU5bPPPhNTpkwRzzzzjDh06JD4\nwx/+ICZPnizef/99IUR8vTrnnHNEQ0ODop3C4bDjyiKEtT7lpLL09PQoyrBnzx4xffp0sWrVKiGE\ns9ol07bXSeuKo4Otz+cTU6ZMUXTmp556Slx77bV5zJW6LVu2iEmTJone3t7Ed3/84x/F7NmzhRBC\nzJo1S2zatEl13ieeeEJRJr/fL6ZPn64od67dfvvt4mc/+1na9xs2bBBTpkwRXq838d1NN90kHnvs\nMSGEM8syYNWqVeKCCy4QXV1dQgjnt8n+/fvFV77yFXH11VcrNoRW2+C6665LTCtEvO9OnjxZ0Xdz\nVZannnpKfOtb31JM++CDD4of/OAHQggh9uzZI84991wRCoVU03VSWYSw1qecVpZkCxYsEDfffHPi\ns5PaJdO210nriqOHkffu3YtQKIQZM2YkvpsxYwZ27tyJaDSax5ylGzduHFauXIlhw4YlvnO5XOjp\n6YHH40FXVxfGjh2rOu/27dsxc+bMxOfq6mpMmjQJW7dulZ5vLQcOHFDN7/bt23Huuedi+PDhie9m\nzJiBbdu2JX53WlkAoLe3Fz//+c9x5513ora2tiDaZPPmzbjgggvwxhtvpOXNbBtEo1Hs3LlT8fu0\nadMQjUYVQ4C5Ksu8efOwaNEixXcD6w0ANDY24vTTT0dFRUVamk4ri5U+5bSyJNu6dSvee+893Hff\nfYnvnNQumba9TlpXys0ULlc8Hg9qa2tRVVWV+O7kk09GOBxGe3s7Ro8encfcKY0aNQqzZ89OfI7F\nYnj11Vcxe/ZsHDhwAOXl5XjyySexfv161NXV4eabb8Y3vvENAPFyppblpJNOwokTJ3JahgGhUAiH\nDx/G2rVr8eSTT0IIgSuuuAJ33nmnZl5bWloAOK8sA9544w1UVlbimmuuAYCCaJPrrrtO9XsrbdDT\n04NgMKj4vby8HCNHjkzML4NWWVIDU1tbG955553E9QEHDhyA2+3GLbfcgj179mDs2LH493//d0yd\nOtVxZbHSp5xWlmTPPvssLrvsMkycODHxnZPaJdO210nriqODbSAQQGVlpeK7gc+hUCgfWdJtyZIl\n2LNnD9566y1s3rwZAHD22WfjhhtuwObNm7F48WJUV1dj3rx5muXMVxmbmpoQiUQwdOhQLF++HM3N\nzXjkkUfg8/kQDAbT9mYrKysRDocBaLdZPttLCIE33ngD119/fSLvBw8eBFA4bZIsEAiYboO+vr7E\nZ7Xf88nv9+OOO+7A6NGjE0GgsbERPT09+NGPfoQxY8bgzTffxE033YTVq1fD7XYDcE5ZrPQpp7bL\n0aNHsX79erz++uuK753cLsnb3pdeeskx64qjg21VVVVaoQY+V1dX5yNLWQkh8Mgjj+C1117Dk08+\niQkTJuDMM8/E/PnzMXLkSADxlbGpqQmvvfYa5s2bp1nOgelzbcKECdi0aRPq6uoAxPMrhMDdd9+N\na665Br29vYrpQ6EQhgwZAkC7zfJVFgDYtWsXmpub8dWvfjXx3XXXXVdQbZKsqqrKdBsMjBKp/T4w\nfz54vV7ceuutOHLkCH71q18l1u9ly5YhGAwmhgEffvhhfPzxx/j973+fCMhOKYuVPuXUdlmzZg2+\n8IUvYOrUqYrvndguatteJ60rjj5nO2bMGPT09CgK6/F4UFlZidra2jzmTF0sFsP999+P119/HY8/\n/jguvfRSAPHzB6kb6XHjxiWGJMeMGQOPx6P4va2tDfX19bnJuIqBQDtg/PjxCIfDGD16dMa8OrEs\n69evx9SpUzFmzJjEd4XYJgOy5S3T7wMbkba2tsRvkUgEXV1deTst09HRgRtvvBGHDx/GL37xC3zh\nC19I/FZRUaE43+ZyuTBu3Di0trY6rixW+pTTyjJg/fr1uOyyy9K+d1q7aG17nbSuODrYnnPOOaio\nqFBclNLQ0IBJkyahvNx5B+VLly7F22+/jeXLlys66NKlS3Hrrbcqpt2zZw/GjRsHAJg6dSo+/vjj\nxG+BQAC7d+/GtGnTcpPxFH/+858xe/ZsxU7O7t27MWLECEybNg179+5V3GvW0NCQyKvTygKkXwQB\nFF6bJJs6darpNigrK8OUKVPQ0NCQ+H3btm1wu90455xzcleIfqFQCLfddhs6Ozvxy1/+MlH/A775\nzW9i5cqVic+xWAyffvopxo0b57iyWOlTTisLED9S3LFjR9q6AzivXbS2vY5aVwxfv5xjixYtEvPm\nzRPbt28X7733njj//PPFO++8k+9spdm6dauYOHGiWLFiheK+s9bWVrFp0yZx9tlni1deeUU0NTWJ\nVatWiUmTJoktW7YIIYQ4fPiwmDJlinj66afF/v37xV133SXmz5+ft/tsOzo6xBe/+EWxcOFCcfDg\nQbF27VoxZ84c8cwzz4hIJCKuvPJKsWDBArFv3z6xYsUKMXXqVHH48GFHlkUIIebOnSt+97vfKb4r\ntDZJvi3DahusXr1aTJs2TaxZs0bs2LFDXHXVVeKhhx7KS1lWrFghzj33XLFhwwbFOtPZ2SmEEGL5\n8qh7/uIAAAF4SURBVOVi1qxZYt26daKxsVEsWrRIfPGLXxQ9PT2OK4vVPuWksgzkd+LEieLYsWNp\n0zqpXTJte520rjg+2Pr9fnHPPfeIadOmiTlz5ogXXngh31lStXTpUjFx4kTVf+FwWLzzzjti/vz5\nYvLkyWLevHlizZo1ivnXrVsnLr/8cnHeeeeJG264QTQ1NeWpJHG7du0S119/vZg2bZq46KKLxPLl\ny0UsFhNCCHHo0CHxL//yL2Ly5MniyiuvFB988IFiXqeVZcqUKWLt2rVp3xdSm6RuCK22wYoVK8SF\nF14oZsyYIe69914RCARyUg4hlGX5+te/rrrODNz7GIlExJNPPim+/OUviylTpojrr79e7N2715Fl\nEcJ6n3JSWbZt2yYmTpwofD5f2rROapds216nrCsuIYQwcdROREREOjn6nC0REVExYLAlIiKSjMGW\niIhIMgZbIiIiyRhsiYiIJGOwJSIikozBloiISDIGWyIiIskYbImIiCT7/z520t8xaV8vAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1115e17f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(tf[:10, :].toarray().T);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create training set\n",
    "We split the whole documents into training and test sets. The number of tokens in the training set is 610K. Sparsity of the term-frequency document matrix is 0.016%, which implies almost all components in the term-frequency matrix is zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of docs for training = 10000\n",
      "Number of docs for test = 1314\n",
      "Number of tokens in training set = 610099\n",
      "Sparsity = 0.01677095\n"
     ]
    }
   ],
   "source": [
    "n_samples_tr = 10000\n",
    "n_samples_te = tf.shape[0] - n_samples_tr\n",
    "docs_tr = tf[:n_samples_tr, :]\n",
    "docs_te = tf[n_samples_tr:, :]\n",
    "print('Number of docs for training = {}'.format(docs_tr.shape[0]))\n",
    "print('Number of docs for test = {}'.format(docs_te.shape[0]))\n",
    "\n",
    "n_tokens = np.sum(docs_tr[docs_tr.nonzero()])\n",
    "print('Number of tokens in training set = {}'.format(n_tokens))\n",
    "print('Sparsity = {}'.format(\n",
    "    len(docs_tr.nonzero()[0]) / float(docs_tr.shape[0] * docs_tr.shape[1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log-likelihood of documents for LDA\n",
    "For a document dd consisting of tokens *w*, the log-likelihood of the LDA model with *K* topics is given as\n",
    "$$\\begin{eqnarray}\n",
    "    \\log p\\left(d|\\theta_{d},\\beta\\right) & = & \\sum_{w\\in d}\\log\\left[\\sum_{k=1}^{K}\\exp\\left(\\log\\theta_{d,k} + \\log \\beta_{k,w}\\right)\\right]+const,\n",
    "\\end{eqnarray}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def logp_lda_doc(beta, theta):\n",
    "    \"\"\"Returns the log-likelihood function for given documents.\n",
    "\n",
    "    K : number of topics in the model\n",
    "    V : number of words (size of vocabulary)\n",
    "    D : number of documents (in a mini-batch)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    beta : tensor (K x V)\n",
    "        Word distributions.\n",
    "    theta : tensor (D x K)\n",
    "        Topic distributions for documents.\n",
    "    \"\"\"\n",
    "    def ll_docs_f(docs):\n",
    "        dixs, vixs = docs.nonzero()\n",
    "        vfreqs = docs[dixs, vixs]\n",
    "        ll_docs = vfreqs * pm.math.logsumexp(\n",
    "            tt.log(theta[dixs]) + tt.log(beta.T[vixs]), axis=1).ravel()\n",
    "\n",
    "        # Per-word log-likelihood times num of tokens in the whole dataset\n",
    "        return tt.sum(ll_docs) / tt.sum(vfreqs) * n_tokens\n",
    "\n",
    "    return ll_docs_f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the inner function, the log-likelihood is scaled for mini-batches by the number of tokens in the dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the log-likelihood function, we can construct the probabilistic model for LDA. doc_t works as a placeholder to which documents in a mini-batch are set.\n",
    "\n",
    "For ADVI, each of random variables $θ$ and $β$, drawn from Dirichlet distributions, is transformed into unconstrained real coordinate space. To do this, by default, PyMC3 uses a centered stick-breaking transformation. Since these random variables are on a simplex, the dimension of the unconstrained coordinate space is the original dimension minus 1. For example, the dimension of $θ_d$ is the number of topics *(n_topics)* in the LDA model, thus the transformed space has dimension *(n_topics - 1)*. It shuold be noted that, in this example, we use *t_stick_breaking*, which is a numerically stable version of stick_breaking used by default. This is required to work ADVI for the LDA model.\n",
    "\n",
    "The variational posterior on these transformed parameters is represented by a spherical Gaussian distributions *(meanfield approximation)*. Thus, the number of variational parameters of $θ_d$, the latent variable for each document, is *2 * (n_topics - 1)* for means and standard deviations.\n",
    "\n",
    "In the last line of the below cell, *DensityDist* class is used to define the log-likelihood function of the model. The second argument is a Python function which takes observations (a document matrix in this example) and returns the log-likelihood value. This function is given as a return value of *logp_lda_doc(beta, theta)*, which has been defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Applied stickbreaking-transform to theta and added transformed theta_stickbreaking_ to model.\n",
      "Applied stickbreaking-transform to beta and added transformed beta_stickbreaking_ to model.\n"
     ]
    }
   ],
   "source": [
    "n_topics = 10\n",
    "minibatch_size = 128\n",
    "\n",
    "# Tensor for documents\n",
    "doc_t = shared(np.zeros((minibatch_size, n_words)).astype('float32'), name='doc_t')\n",
    "\n",
    "with pm.Model() as model:\n",
    "    theta = Dirichlet('theta', a=(1.0 / n_topics) * np.ones((minibatch_size, n_topics)).astype('float32'),\n",
    "                      shape=(minibatch_size, n_topics), transform=t_stick_breaking(1e-9))\n",
    "    beta = Dirichlet('beta', a=(1.0 / n_topics) * np.ones((n_topics, n_words)).astype('float32'),\n",
    "                     shape=(n_topics, n_words), transform=t_stick_breaking(1e-9))\n",
    "    doc = pm.DensityDist('doc', logp_lda_doc(beta, theta), observed=doc_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def create_minibatch(data):\n",
    "    rng = np.random.RandomState(0)\n",
    "\n",
    "    while True:\n",
    "        # Return random data samples of a size 'minibatch_size' at each iteration\n",
    "        ixs = rng.randint(data.shape[0], size=minibatch_size)\n",
    "        yield [data[ixs]]\n",
    "\n",
    "minibatches = create_minibatch(docs_tr.toarray().astype('float32'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The value of doc_t will be replaced with mini-batches\n",
    "minibatch_tensors = [doc_t]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# observed_RVs = OrderedDict([(doc, n_samples_tr / minibatch_size)])\n",
    "observed_RVs = OrderedDict([(doc, 1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encoder\n",
    "\n",
    "Given a document, the encoder calculates variational parameters of the (transformed) latent variables, more specifically, parameters of Gaussian distributions in the unconstrained real coordinate space. \n",
    "`encode()` method is required to output variational means and stds as a tuple, as shown in the following code. As explained above, the number of variational parameters is `2 * (n_topics) - 1`. Specifically, the shape of `zs_mean`(or `zs_std`) in the method is `(minibatch_size, n_topics - 1)`. It should be noted that `zs_std` is defined as log-transformed standard deviation and this is automatically exponentiated (or bounded to be positive) in `advi_minibatch()`, the estimation function. \n",
    "\n",
    "To enhance generalization ability to unseen words, a Bernoulli corruption process is applied to the inputted documents. Unfortunately, I've never seen any significant improvement with this. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class LDAEncoder:\n",
    "    \"\"\"Encode (term-frequency) document vectors to variational\n",
    "    means and (log-transformed) stds.\n",
    "    \"\"\"\n",
    "    def __init__(self, n_words, n_hidden, n_topics, p_corruption=0,random_seed=1):\n",
    "        rng = np.random.RandomState(random_seed)\n",
    "        self.n_words = n_words\n",
    "        self.n_hidden = n_hidden\n",
    "        self.n_topics = n_topics\n",
    "        self.w0 = shared(0.01 * rng.randn(n_words, n_hidden).ravel(), name='w0')\n",
    "        self.b0 = shared(0.01 * rng.randn(n_hidden), name='b0')\n",
    "        self.w1 = shared(0.01 * rng.randn(n_hidden, 2 * (n_topics - 1)).ravel(), name='w1')\n",
    "        self.b1 = shared(0.01 * rng.randn(2 * (n_topics - 1)), name='b1')\n",
    "        self.rng = MRG_RandomStreams(seed=random_seed)\n",
    "        self.p_corruption = p_corruption\n",
    "\n",
    "    def encode(self, xs):\n",
    "        if 0 < self.p_corruption:\n",
    "            dixs, vixs = xs.nonzero()\n",
    "            mask = tt.set_subtensor(\n",
    "                tt.zeros_like(xs)[dixs, vixs],\n",
    "                self.rng.binomial(size=dixs.shape, n=1, p=1-self.p_corruption)\n",
    "            )\n",
    "            xs_ = xs * mask\n",
    "        else:\n",
    "            xs_ = xs\n",
    "\n",
    "        w0 = self.w0.reshape((self.n_words, self.n_hidden))\n",
    "        w1 = self.w1.reshape((self.n_hidden, 2 * (self.n_topics - 1)))\n",
    "        hs = tt.tanh(xs_.dot(w0) + self.b0)\n",
    "        zs = hs.dot(w1) + self.b1\n",
    "        zs_mean = zs[:, :(self.n_topics - 1)]\n",
    "        zs_std = zs[:, (self.n_topics - 1):]\n",
    "        return zs_mean, zs_std\n",
    "\n",
    "    def get_params(self):\n",
    "        return [self.w0, self.b0, self.w1, self.b1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To feed the output of the encoder to the variational parameters of $\\theta$, we set an OrderedDict of tuples as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "encoder = LDAEncoder(n_words=n_words, n_hidden=100, n_topics=n_topics, p_corruption=0.0)\n",
    "local_RVs = OrderedDict([(theta, (encoder.encode(doc_t), n_samples_tr / minibatch_size))])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "theta is the random variable defined in the model creation and is a key of an entry of the OrderedDict. The value (encoder.encode(doc_t), n_samples_tr / minibatch_size) is a tuple of a theano expression and a scalar. The theano expression encoder.encode(doc_t) is the output of the encoder given inputs (documents). The scalar n_samples_tr / minibatch_size specifies the scaling factor for mini-batches.\n",
    "\n",
    "ADVI optimizes the parameters of the encoder. They are passed to the function for ADVI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "encoder_params = encoder.get_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AEVB with ADVI\n",
    "`advi_minibatch()` can be used to run AEVB with ADVI on the LDA model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'v_params' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-705c6e539bd9>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;31m#%timeit v_params = run_advi()\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv_params\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0melbo_vals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'v_params' is not defined"
     ]
    }
   ],
   "source": [
    "def run_advi():\n",
    "    with model:\n",
    "        v_params = pm.variational.advi_minibatch(\n",
    "            n=1000, minibatch_tensors=minibatch_tensors, minibatches=minibatches,\n",
    "            local_RVs=local_RVs, observed_RVs=observed_RVs, encoder_params=encoder_params,\n",
    "            learning_rate=1e-2, epsilon=0.1, n_mcsamples=1\n",
    "        )\n",
    "\n",
    "    return v_params\n",
    "\n",
    "#%timeit v_params = run_advi()\n",
    "\n",
    "plt.plot(v_params.elbo_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'zip' object is not subscriptable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-25-aed8ca2df4b0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mv_params\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrun_advi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-23-705c6e539bd9>\u001b[0m in \u001b[0;36mrun_advi\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m             \u001b[0mn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mminibatch_tensors\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mminibatch_tensors\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mminibatches\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mminibatches\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m             \u001b[0mlocal_RVs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlocal_RVs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobserved_RVs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mobserved_RVs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mencoder_params\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mencoder_params\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m             \u001b[0mlearning_rate\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepsilon\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_mcsamples\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m         )\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/peadarcoyle/anaconda/envs/pymc3_dev/lib/python3.5/site-packages/pymc3-3.0rc2-py3.5.egg/pymc3/variational/advi_minibatch.py\u001b[0m in \u001b[0;36madvi_minibatch\u001b[0;34m(vars, start, model, n, n_mcsamples, minibatch_RVs, minibatch_tensors, minibatches, local_RVs, observed_RVs, encoder_params, total_size, optimizer, learning_rate, epsilon, random_seed, verbose, local_NFs, global_NFs)\u001b[0m\n\u001b[1;32m    329\u001b[0m     \u001b[0mlogp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtheano\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclone\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlogpt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreplace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstrict\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    330\u001b[0m     elbo = _elbo_t(logp, uw_g, uw_l, inarray_g, inarray_l,\n\u001b[0;32m--> 331\u001b[0;31m                    n_mcsamples, random_seed, local_NFs, global_NFs)\n\u001b[0m\u001b[1;32m    332\u001b[0m     \u001b[0;32mdel\u001b[0m \u001b[0mlogpt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    333\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/peadarcoyle/anaconda/envs/pymc3_dev/lib/python3.5/site-packages/pymc3-3.0rc2-py3.5.egg/pymc3/variational/advi_minibatch.py\u001b[0m in \u001b[0;36m_elbo_t\u001b[0;34m(logp, uw_g, uw_l, inarray_g, inarray_l, n_mcsamples, random_seed, local_NFs, global_NFs)\u001b[0m\n\u001b[1;32m    184\u001b[0m         logps, _ = theano.scan(fn=lambda z_g, z_l: logp_(z_g, z_l),\n\u001b[1;32m    185\u001b[0m                                \u001b[0moutputs_info\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 186\u001b[0;31m                                sequences=zip(zs_g, zs_l))\n\u001b[0m\u001b[1;32m    187\u001b[0m         \u001b[0melbo\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mtt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlogps\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    188\u001b[0m             \u001b[0mtt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw_g\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m0.5\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0ml_g\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2.0\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/peadarcoyle/anaconda/envs/pymc3_dev/lib/python3.5/site-packages/Theano-0.8.2-py3.5.egg/theano/scan_module/scan.py\u001b[0m in \u001b[0;36mscan\u001b[0;34m(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict)\u001b[0m\n\u001b[1;32m    473\u001b[0m                 \u001b[0;31m# If not we need to use copies, that will be replaced at\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    474\u001b[0m                 \u001b[0;31m# each frame by the corresponding slice\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 475\u001b[0;31m                 \u001b[0mactual_slice\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mseq\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'input'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mmintap\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    476\u001b[0m                 \u001b[0m_seq_val\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtensor\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_tensor_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseq\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'input'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    477\u001b[0m                 \u001b[0m_seq_val_slice\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_seq_val\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mmintap\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: 'zip' object is not subscriptable"
     ]
    }
   ],
   "source": [
    "v_params = run_advi()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Variational Inference\n",
    "\n",
    "Angelino, Elaine, Matthew James Johnson, and Ryan P. Adams. \"[Patterns of Scalable Bayesian Inference](https://arxiv.org/pdf/1602.05221v2.pdf).\" arXiv preprint arXiv:1602.05221 (2016).\n",
    "\n",
    "Blei, David M., Alp Kucukelbir, and Jon D. McAuliffe. \"[Variational inference: A review for statisticians](http://arxiv.org/pdf/1601.00670).\" arXiv preprint arXiv:1601.00670 (2016).\n",
    "\n",
    "Kucukelbir, Alp, et al. \"[Automatic Differentiation Variational Inference](https://arxiv.org/pdf/1603.00788.pdf).\" arXiv preprint arXiv:1603.00788 (2016).\n",
    "\n",
    "Ranganath, Rajesh, Sean Gerrish, and David M. Blei. \"[Black Box Variational Inference]((http://www.cs.columbia.edu/~blei/papers/RanganathGerrishBlei2014.pdf).\" AISTATS. 2014."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Models and Methods\n",
    "Diederik P Kingma and Max Welling \"[Auto-Encoding Variational Bayes](https://arxiv.org/abs/1312.6114).\"\tarXiv preprint arXiv:1312.6114 (2014).\n",
    "\n",
    "Blei,David M. et al. \"[Latent Dirichlet Allocation](https://www.cs.princeton.edu/~blei/papers/BleiNgJordan2003.pdf).\" Journal of Machine Learning Research 3 (2003) 993-1022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:pymc3_dev]",
   "language": "python",
   "name": "conda-env-pymc3_dev-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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}