lda-advi-ae.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Automatic autoencoding variational Bayes for latent dirichlet allocation with PyMC3\n",
    "\n",
    "(c) 2016 by Taku Yoshioka\n",
    "\n",
    "For probabilistic models with latent variables, autoencoding variational Bayes (AEVB; Kingma and Welling, 2014) is an algorithm which allows us to perform inference efficiently for large datasets with an encoder. 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. Combining AEVB with ADVI (Kucukelbir et al., 2015), we can perform posterior inference on almost arbitrary probabilistic models involving continuous latent variables. \n",
    "\n",
    "I have implemented AEVB for ADVI with mini-batch on PyMC3. To demonstrate flexibility of this approach, we will apply this to latent dirichlet allocation (LDA; Blei et al., 2003) for modeling documents. 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. 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. I hope this notebook will attract readers, especially practitioners working on a variety of machine learning tasks, to probabilistic programming and PyMC3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/taku-y/anaconda/anaconda/envs/py35con/lib/python3.5/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.\n",
      "  \"`IPython.html.widgets` has moved to `ipywidgets`.\", ShimWarning)\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sys, os\n",
    "sys.path.insert(0, os.path.expanduser('~/git/github/taku-y/pymc3'))\n",
    "\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.advi import advi, sample_vp\n",
    "\n",
    "theano.config.compute_test_value = 'ignore'"
   ]
  },
  {
   "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 1000.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading dataset...\n",
      "done in 1.896s.\n",
      "Extracting tf features for LDA...\n",
      "done in 1.885s.\n"
     ]
    }
   ],
   "source": [
    "# The number of words in the vocaburary\n",
    "n_words = 1000\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 1000-dimensional term-frequency vector. Let's check the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAFVCAYAAADCLbfjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXucHFWZ93/dc5/pyWXCZEVkAwZZDbuKoq4sSmZBRCBG\nXDMCQsAXdnlx0Q8rrisXP6i4LpHPrq8vCi7ou2bFXVEuu0FEFzESTSAYgwRJSAIhCcnk1nPv7rl1\nd533j+6urnPqnFOnqqt7qqef7wcyXVWnznnO9TnPuVWMMcZAEARBEETkic+2AARBEARBmEFKmyAI\ngiDqBFLaBEEQBFEnkNImCIIgiDqBlDZBEARB1AmktAmCIAiiTjBS2kNDQ+jr68PevXu5++vXr8eq\nVatw6aWX4sEHH6yKgARBEARBFGj2cpDL5fDFL34R7e3trvtr1qzBI488gra2Nlx22WU499xz0dPT\nUzVhCYIgCKKR8bS0v/a1r+Gyyy7D4sWLuft79uzBkiVLkEgk0NLSgjPOOANbtmypmqAEQRAE0eho\nlfYjjzyCRYsW4ayzzoJ4cFo6nUZ3d7d93dXVhVQqVR0pCYIgCILwVtqbNm3C6tWrsXPnTnz+85/H\n0NAQACCRSCCdTttuM5kM5s2b5xkgnZpKEISTJ57djw99dh0+9Nl1sy0KQUQe7Zz2D37wA/v36tWr\ncfvtt2PRokUAgKVLl2L//v0YHx9He3s7tmzZgmuuucYzwFgshmSSLPJq0tvbTWlcAyidwyGVmrJ/\ni+lJaVx9KI1rQ29vt7cjAzwXopWIxWIAgMceewyTk5Po7+/HzTffjKuvvhqMMfT397vmvQmCIAiC\nCA9jpf39738fAHDyySfb9/r6+tDX1xe6UARBEARBuKHDVQiCIAiiTiClTRAEQRB1AiltgiAIgqgT\nSGkTBEEQRJ1ASpsgCIIg6gRS2gRBEARRJ5DSJgiCIIg6gZQ2QRAEQdQJpLQJgiAIok4gpU0QBEEQ\ndQIpbYIgCIKoE0hpEwRBEESdQEqbIAiCIOoEUtoEQRAEUSeQ0iYIgiCIOoGUNkEQBEHUCaS0CYIg\nCKJOIKVNEARBEHUCKW2CIAiCqBNIaRMEQRBEnUBKmyAIgiDqBFLaBEEQBFEnkNImCIIgiDqBlDZB\nEARB1AmktAmCIAiiTiClTRAEQRB1AiltgiAIgqgTSGkTBEEQRJ1ASpsgCIIg6gRS2gRBEARRJ5DS\nJgiCIIg6odnLgWVZ+MIXvoC9e/ciHo/jy1/+Mk455RT7+dq1a/HQQw+hp6cHAHD77bfjpJNOqprA\nBEEQBNGoeCrt9evXIxaL4Yc//CF++9vf4utf/zruuece+/n27dtx5513YtmyZVUVlCAIgiAaHU+l\n/f73vx/nnHMOAGBgYADz58/nnm/fvh333nsvkskk+vr6cO2111ZHUoIgCIJocDyVNgDE43HcdNNN\nePLJJ3HXXXdxzy666CJcfvnlSCQSuP7667FhwwYsX75c619vb3dwiQkjKI1rA6Vz5XR3t9u/ZelJ\naVx9KI3rByOlDQBr1qzB0NAQ+vv78fjjj6O9vVDRrrrqKiQSCQDA8uXLsWPHDk+lnUymKhCZ8KK3\nt5vSuAZQOodDKjVl/xbTk9K4+lAa14awOkaeq8fXrVuH++67DwDQ1taGeDyOeLzwWjqdxooVKzA5\nOQnGGDZv3ozTTjstFMEIgiAIguDxtLQ/8IEP4Oabb8YVV1yBXC6HW265BU888QQmJyfR39+PG2+8\nEatXr0ZbWxvOPPNMnH322bWQmyAIgiAaDk+l3dHRgW984xvK5ytXrsTKlStDFYogCIIgCDd0uApB\nEARB1AmktAmCIAiiTiClTRAEQRB1AiltgiAIgqgTSGkTBEEQRJ1ASpsgCIIg6gRS2gRBEARRJ5DS\nJgiCIIg6gZQ2QRAEQdQJpLQJgiAIok4gpU0QBEEQdQIpbYIgCIKoE0hpEwRBEESdQEqbIAiCIOoE\nUtoEQRAEUSeQ0iYIgiCIOoGUNkEQBEHUCaS0CYIgCKJOIKVNEARBEHUCKW2CIAiCqBNIaRMEQRBE\nnUBKmyAIgiDqBFLaBEEQBFEnkNImCIIgiDqBlDZBEARB1AmktAmCIAiiTiClTRAEQRB1AiltgiAI\ngqgTSGkTBEEQRJ1ASpsgCIIg6gRPpW1ZFm655RZcdtlluPzyy/HKK69wz9evX49Vq1bh0ksvxYMP\nPlg1QQmCIAii0fFU2uvXr0csFsMPf/hD3HDDDfj6179uP8vlclizZg3Wrl2L+++/Hz/60Y8wPDxc\nVYEJgiAIolHxVNrvf//78ZWvfAUAMDAwgPnz59vP9uzZgyVLliCRSKClpQVnnHEGtmzZUj1pCYLw\nhTU9jezwEAAgZ+UwODkExhhmjhwBY2yWpSNqTXp8CtmZ3GyLMaeYmc4hnZquWXjNJo7i8Thuuukm\nPPnkk7jrrrvs++l0Gt3d3fZ1V1cXUqlU+FISBBGIvbd+HvnRUZxyz3341x3fx0vDu/H57FmYevC/\ncNxfrULPhStmW0SiRuTzFu6/ZzPa2ptx9d+9d7bFmTOs/ebTyOcsfPKmvpqEZ6S0AWDNmjUYGhpC\nf38/Hn/8cbS3tyORSCCdTttuMpkM5s2b5+lXb2+3pxuiMiiNa0PU03n36CgAoCfRjJeGdwMApna8\nUPj7wu/Re9VlsyZbiUSi3f4tS8+op3G9MDNdsLCnp3KuNKU0Dk4+ZwGoXRp6Ku1169bh6NGjuPba\na9HW1oZ4PI54vDCqvnTpUuzfvx/j4+Nob2/Hli1bcM0113gGmkySNV5Nenu7KY1rQD2l8+BQuXOd\ny1loBpDNWZGQP52esn+L8tRTGked7Eze/u1MU0rjcPBKw7CUuqfS/sAHPoCbb74ZV1xxBXK5HG65\n5RY88cQTmJycRH9/P26++WZcffXVYIyhv78fixcvDkUwgiAIgiB4PJV2R0cHvvGNbyif9/X1oa+v\nL0yZCIIgCIKQQIerEARBEESdQEqbIAiCICqkVlsoSWkTBEEQRJ1ASpsgCIIg6gRS2gRBEARRJ5DS\nJgiCIIgKqdWpwKS0CYIgCKJOIKVNEATRENAHYqoLrR4nCIIgCMIBKW2CIIgGgL7EWl1oTpsgCIIg\nCA5S2gRBEARRJ5DSJgiCIIg6gZQ2QRBEA0Bz2lWG5rQJgiAIgnBCSpsgCKIhIFO7mjDap00QBEEQ\nhBNS2gTRCJCRRRDVhea0CYIgCIJwQkqbIBoBWjrc8FARqC61Sl5S2gRBEARRJ5DSJoiGgMwsgqgq\nNKdNEARBEIQTUtoE0QiQod3wMJrUrjK0T5sgCIIgCAektAmiISAriyCqCX1PmyAIgiAIDlLaBNEA\n0HQmQWVgbkBKmyAIgiDqBFLaBNEIkJlFEFWlVlWsWfcwl8vhlltuwcDAALLZLK677jqcc8459vO1\na9fioYceQk9PDwDg9ttvx0knnVRVgQmCIAiiUdEq7UcffRQLFy7EnXfeibGxMVx88cWc0t6+fTvu\nvPNOLFu2rOqCEgRRCWRpNzw02lJlapO+WqV9wQUX4IMf/CAAwLIsNDfzzrdv3457770XyWQSfX19\nuPbaa6snKUEQBEE0OFql3dHRAQBIp9O44YYb8JnPfIZ7ftFFF+Hyyy9HIpHA9ddfjw0bNmD58uXV\nk5YgiGCQkdXwUBGYG2iVNgAcPnwYn/rUp3DFFVfgwgsv5J5dddVVSCQSAIDly5djx44dRkq7t7c7\noLiEKZTGtSHq6by7+HfRoi77XnNzYf1pS3M8EvJ3d7fbv2XyREHGuUB7W4v9W0xTSuPKWbQogY7O\n1qqHo1Xag4ODuOaaa3DbbbfhPe95D/csnU5jxYoV+NnPfob29nZs3rwZq1atMgo0mUwFl5jwpLe3\nm9K4BtRTOg8Npe3fuayFZgDZbD4S8qdSU/ZvUZ56SuOok0lN27+daUppHA6Dg2m0d7Qon4fVMdIq\n7XvvvRfj4+O45557cPfddyMWi+FjH/sYJicn0d/fjxtvvBGrV69GW1sbzjzzTJx99tmhCEUQBEEQ\nhBut0r711ltx6623Kp+vXLkSK1euDF0ogiBChlYONzxUAuYGdLgKQRAEQdQJpLQJohEgM4ug0Zaq\nUqvvlZPSJgiCIIg6gZQ2QTQAjExtgpgTkNImCIIgiDqBlDZBNAJkaDc8NKVdXWqVvqS0CYIgCKJO\nIKVNEI0AmVkEMScgpU0QBEEQdQIpbYJoCMjSbnRqtY+4YaF92gRBEARBOCGlTRCNABlZBFFValXF\nSGkTBEEQRJ1ASpsgGgIytQmiqtA+bYIgCIIgnJDSJohGgFYONzxUBKoLzWkTBEEQBMFBSpsgGoAo\nW1m0f7hWUDpXFdqnTRAEQRCEE1LaBNEQkJXV6NCAxtyAlDZBELMK6RKCMIeUNkE0AqQZCaKq0Pe0\nCYIgCILgIKVNEI1GrPQ3pnVWM2gUgCCMIaVNEA0BacZGhxaizQ1IaRMEMauQLiEIc0hpE0QjQGYW\nQd2jqlKrQ4JIaRMEQRBEnUBKmyAagSgbWTQKUBsomecEpLQJgiAIok4gpU0QDQGZWY0OlYC5QbPu\nYS6Xwy233IKBgQFks1lcd911OOecc+zn69evxz333IPm5mZ89KMfRX9/f9UFJghibkHKhCDM0Srt\nRx99FAsXLsSdd96JsbExXHzxxbbSzuVyWLNmDR555BG0tbXhsssuw7nnnouenp6aCE4QhDn0+UuC\nekfVJRLHmF5wwQW44YYbAACWZaG5uazj9+zZgyVLliCRSKClpQVnnHEGtmzZUl1piaqxZ2cSySOp\nWZVhJDWN9c8dhCUp/RO7dyHzhxdCC+u11EE8dyw8/4IwkD6M3x19viZhpSdzymd7x17DtuSLyucz\nE4cwMbIDmfQ0/vtb/4kX1q8PRaZMdgJPHdyEnKWWTcdUbhpPHdyEydxkKPIQwdj60gbs+8W6yHcM\nx595GtOHBmZbjIrRWtodHR0AgHQ6jRtuuAGf+cxn7GfpdBrd3d32dVdXF1Ips0a/t7fb2xFREX7S\nmFkM3/7vpwAAt/3Lh6okkTe3/dtvcfBYGq/r7cbyd7yBe7bpr+8AAJy17uFQwrp+/V0AgPf3fwtN\n8abA/lRSlq9f/w8AgOWnvhOdrR2B/dGxu/j3oadeAd5U+N3cXOirtzTH0dvbjevXfwsA8ONLvi31\nY+sT3wUAHDiyAofTr0fm10mce0nldfjfN/0nfnvwefxZ93sALAAgT09VGq/9/c/x+O71ODp9BJ96\nzycqlmeuk53J2r/FNA1ajieyk2j7v2sxk2NoP/0MzFv2lopkrBbTyUHs/n/3AQivDRHp6elCz3Fd\nVfHbiVZpA8Dhw4fxqU99CldccQUuvPBC+34ikUA6nbavM5kM5s2bZxRoMjm7Ft1cp7e321caW1a5\nhzybeXPwWKE8vXZoDMkT50vdhC1fMpkKrLT9prOKI8lRdLcGszZNOXgsZSvtXM5CM4BszuLk94rL\n2EghfyZbwon3vuGDAICRqWGUlLbory6N9w8dKvozQG2KAUPD5fbamV6VlOPUTBqtuUL7MXTwGKZ7\n3+DxxuwwfWjI/l2tsjI8nEaeWcrnYRmr2uHxwcFBXHPNNfjc5z6Hj3zkI9yzpUuXYv/+/RgfH8fM\nzAy2bNmC008/PRShCIKIIKW+XcSHQQliNqhVtdBa2vfeey/Gx8dxzz334O6770YsFsPHPvYxTE5O\nor+/HzfffDOuvvpqMMbQ39+PxYsX10ZqImSoESa8oVJS70Tkq25ERWiV9q233opbb71V+byvrw99\nfX1hy0TUmEY2nFiDqKIwm2tq+uuUxijqcx46XIVoaKgdiwCN3GusKZTOcwFS2gTREITZYFPjX49Q\nrlWXSOzTJhqERq7NZOUZQ0lV31D2zQ1IaRMNM68ro1FiHgslogVPYiGnWqPkwWwT9cNP6h/6njZR\nKxq6Ljd05ImGgsr6XICUNkE0BEz6k2ggKN+rC81pE7WiketyI8ddxHP4lAl/ibqCsm1uQEqbaOja\n3JjzfPI4e61tYJJfRD1B+VZNapW6pLQJNHZlbuS4E41EQ/ZP5yCktAmiAXCeYha47WZuvwiCKEJz\n2kStaOQeeANH3YXXVAGlVb1DOTgXIKVNNDiz35DVZp88nYjW6FhV6J038hkPIrVKC1LaRGNb2hGI\nexRkAKgBJvwTlbLbSJDSJogGgOahCeqUVRma0yZqRyNX5ijEPQoyeEvB7IVo0ZCX8ElVso3KQq0h\npU00NFFocmpiATmDUJndpmOdUUg0wjfVOJOArPfaQ0qbaOh5KWp0ysxWSjRy+SMIv5DSJogGgIa0\nCaK60Pe0CaIWRECXRecoVcOzx6OQaIRvqpFr0Sm7jQMpbaKhK14jDo+rsttzIZr9PW2iHmnkel4b\naJ82QTQEUek4GH/liyCIWYOUNkE0AGHMaUelc0EEg1aPVxea0yZqRiOPmkWh0YlO+nucPU5z2oRA\ndMpu40BKm6grwrYWotHo1FoI1fe0vV6LRGIRAalOB7VeykS9yOkNKW2ivghdccx+ZZ59CUrQV74I\nIuqQ0ibIgmoAQlnxXTrGlMpLXVKdOe06oW4E9YaUNjGXyrNvojCnHZUc8Fw8TsqacFEvZaJe5PSG\nlDZRX+V5DiqO2kTJJBDa8jWXqYqlXS9logaC1qpTS0qbqC9CX4g2+61OLax9k2h6H64SMhV6GI1R\nEqIM5UctIKVNUFVrAMKY067WpzlJ+daGavRP6ybv6kRME4yU9rZt27B69WrX/bVr12LFihW48sor\nceWVV2Lfvn1hy0fUgghYm41NNNLfuwEOWc4KexIxOlA1Aph88zUKRKOOhUGzl4Pvfve7WLduHbq6\nulzPtm/fjjvvvBPLli2rinAEIcIYC7VpiIKlUAsJjNLMcyFaGJIQs0U1ynq9FIlaTINF5kS0JUuW\n4O6775Y+2759O+699158/OMfx3333Re6cAThpl6aCR/UokGpeggBoDntOQblRy3wVNrnnXcempqa\npM8uuugifPnLX8b3v/99bN26FRs2bAhdQKL6NLIF1Shxd+6tVkXZSwnSMaZ1TlUmtaks1BrP4XEd\nV111FRKJBABg+fLl2LFjB5YvX+75Xm9vdyXBEgb4SeMWR6csCnmTSLS55Nhd/HvcogSa2tpCC6tn\nURd6u4LHOYz0WtjThd551Un3UrohVh4gb24u9NVbmuOc/IsWdWF+u1uO14p/m+IFP2IsnHg3FeVo\n9ih/qrDaWlsAAC3NTZEot1HnwNEO+7eYXkHTb6Ytg7Hi73nzO7EoovmQHu+0y3G1ysqCBZ01KYfG\nSlucE0in01ixYgV+9rOfob29HZs3b8aqVauM/EomU/6kJHzR29vtK43HRyft31HIm3R6WinHYDKF\neNtMaGENDaURm2gN9K7fdFYxPJxG63R109254juXs9AEIJuzOPkHh9KYaVXPfufyFoAmACyUeOfz\nVlGevH1P9FeXxjMzOQBANpePRLmNOuNjE/ZvZ3pVUo6HM2mH/5OwIpoPU8MZ+3e1ysrISAZtnWqV\nGpZCN1basWJP/bHHHsPk5CT6+/tx4403YvXq1Whra8OZZ56Js88+OxShCKKRqMUAYzijmDHhb4XQ\nnPYcg/KjFhgp7RNOOAEPPPAAAGDFihX2/ZUrV2LlypXVkYwgakJjNDT83mrFV748d3xVJ60aIwdm\nn+rs064X6kdSL+hwFaK+mIOf5qzJqWwh7PkqPWVh7bmjfdp1D192I5wfUajoIUFKm+DKcxSO9dRD\np3EFwaQ5NV89HuHGmVAS/bpdPWoS9ajs0yaIKDEX252onD1ec2hOe45B+VELSGkTcFa2SDbuVaXh\nIhwcSqq6pjonotVLoZg7BxiR0ibqa3g87DntUH2LLjGDEW2vvA89rWhOe9aoTj2PcH5EvV3zASlt\nos4IufJFoDLXoqPEnYgWNDhW+hPhxplQwnXOQ6pHke/kVxnGWzw1CZOUNsER+UoYcfGCUC9rZOh7\n2nOH6tRzyo9aQEqb4DuIDVbvGia6BhH1VIINk1hzk2oo6ropEnPoozyktAnU1UK00AWMQoRrMDxu\nNKft8Zy+pz1nsELLy3rZpz3bAoQHKW1CYA6VbgOiENtaDPOaNadechR8YSY9ACLaRL53Xh/Mxigl\nKW2CK2yRr8uRFzAANYhSGMkWtUGOUmeHLG4znJ3D0BaiGRyPGw2iLJs/SGkTfLWbO2XbiMgvvAsN\ng7PHayNI6NCCNEOY9GdofkaZ6tXz2icAKW1CKHfRroVzsYGuyfB4GGePh32MKc1p1xTO0mZWFfwM\nxUvCA1LaBEfkK17Y27TnYCcgKLOV95Evc3OFaljanP+Nl5GzcTAVKW1C6C1HvOJFXb4A1KTjEGIQ\nURlapTltf/DTYA22T3sOtRuktIlI1zUiLLwz2fwrX9GCRktMCX/RWN2kfVQLbwBIaRMcUSzbrIr7\nKqIwslCTz2mHMKft8K0CScLzhixsfzgVbFj7tBv5YKbZgpQ2wRMBJaYl9DntKBCN4fFopAVRNaqu\nYBuwBM3CdllS2oSwmGL25FAyx7vz9RIjOxvCMnBpTnvWCG9Yu9qr20Iikg1bMEhpE+COMZ1FKWYH\ninH5QX2mRd3Mq84yzqmg8A5XaWxmo+yR0iYEQzaC1TCMz0qqvI5EsxONY0y9pGD2X5rTrnfCW8sR\n9WG6AlFYuxIWpLQJjsiX7fC19qxTmzQ3CcTDTQTSigiOM/uskA5X4f1vwAIyCwdTkdImOCLZI53z\nc9rRWIjmTUz4WyE0p11jwj+PoW6qZhTbtYCQ0iYamoa0DhSYDo9HDcpDM6qhtxo97WdjZpGUNsEv\nUIlgHeSsggjKVw+YfU+7xsPjtkzBPCYLOzgWnYhWt5DSJgSiXrhDPlwlAvGt/ZREZeHNfooRgWDh\nD4+r/G8YZiHKpLQJjmjWu/pYoRqc+pjTLnsR7px2UNFoTtsfTHMV2M+6GQWLtHC+IKVN1MHhKlX0\nOorxrQqSiDKxGa/PxKhXuWcTK7Qkq5O0r5qYtZ9aJKVN8ERRi1W1Nz/78a2FBGZnj+thrOhJWIYt\n7dOuKXw1Cut72uqrKBHJXTEBIaVNcES9aM9Nq6r6cYrkKGZIW74IU6psFTZgdvDpGKF92tu2bcPq\n1atd99evX49Vq1bh0ksvxYMPPhi6cERtYNVeoBImc+9slcgMbpjmfWgnohWp1GImi9sMfntS+Cei\nRbvtiLJs/mj2cvDd734X69atQ1dXF3c/l8thzZo1eOSRR9DW1obLLrsM5557Lnp6eqomLFEDoli2\nq9oYRDHC4eBsRMM8xjRsKrWYyeIOQniHq9RFl6kWaz2jMqe9ZMkS3H333a77e/bswZIlS5BIJNDS\n0oIzzjgDW7ZsqYqQRJWJ+kK0KgoYhfjWQumYNayztE87oL9kYfvD2YkLa582X3YjUJlURKGih4Sn\npX3eeedhYGDAdT+dTqO7u9u+7urqQiqVMg6YMYZfrNuBWDyG93/oLYjFYpjJZ7F2+3/iL098L960\ncCmyOQv33b8Vb+huw1vPOgk/2bQPf73iLehsbwEA5C0L9z66A+976/H4szcuAgBs2XkM214ZxDUX\nFfzcuusYtu5O4q9XLENcshrnxcGX8Nsjz+Gik8/DI6/8FB879WIs6lgIANi84VV0tMfR/Mx/YM+p\nC/ChC/4WsVgM03kLP371CJYfvxB/nOgAAGRS09jw8914z1++ET3HFUYlhn/2Uww+/CBGTu5FItGD\nRZd/Gt/7+S60xhg+NvIMev7yL9F12p9i+95hPPX8AN789jGkcuO4+JQLbfnWbdyLluY4LnzPEk7u\nmfwM7nz6u5gYacZbjp2Oo9k8+v5yKd755sVI7X4Z6x96FM81z8efv8HC1tFluPy8U9G7oEOaF795\n4ZAzZ1zPf/rMPuTzDCvfe7J9L/ngj9C8aBGSi/8Uhw6Moe+CP1Hm9eh0Fo/uT+ItC7uwazSDS5a+\nDi1xvr/4f3/9MJpfP4TcoVO4+9lsFt/89ydw4vHn4OTh57F5wx687jTg1NP+CACw4ee70BZ/GYPD\nx2EIGXTMb8LHLujDscPjePzBP6Dvgj/BSW86DgDwP/vWw2IWLjj5/bb/j+/9BT799r+xrx96ag/m\nJ1px3jtPBABMvbYfx/7zB1iwvA/zzjxLGUcnlmXh2z/5HTqOjuDIcAusP0pgVW47Tjj9NLzKjsfU\nVBZ/fvYbbfeP/HoP/nr5CVi8sNPIf6BQf478v/uQePs70H3Gu1SOsHNJGw7+USvSXSfitMM55Kde\nQDrfgZePPxcnT+xF9qVjOGHPWzHwxhfwxP5f4Zo/vcJ+ffjAzzG+bht+OfNmHHr9IHrtJzH8/LFt\n2Hvic/ir15+OePoV7PmVhdTIFGb++O3YHwfee/oJeOvSRWDMwtC+/8LOzfswuHcY1kmr8MGVy9DW\nNoPszDgA4Bj2APGT8TqrGVt+sxfvet/Jrqj8bNOzGHhtBPNftxTzu9rwgXf/cVlh5PP44X3/jdH5\nR/EX7zoHz20cw5uXn4SZkWkMJDNADMg2j+H3e8fwZ6e9Ede8u5z2YyMT2PTkHrzvA29C9/x2AMBP\nn3oWm2Z+hfaRd+GDy5ZhcOcg/uLcpVjQw+dPauvvcPjb38KCc96PxR+/AhNju3D0wIt4/MnFOG5+\nBz548TJsevIVjKemMTg+haHeIZy56ynkln4Ax5b0YKrnJax+yyXoain4+4N/fgIL2Tgu+LsP497d\nv8fOl3+HjsluWAtHcHz+zxCf6MFVF7wZrbEYfr5uO545mkJiQQc++6FTMPbDf8eiD38EPz+2A88c\nfRG39V2F9OFf4p7fduNNS47DFWe9h6vau15N4g1vP964vG184TBePTyOK88v1PNnXjyCXQdG8efv\nLtTjXUvasDH7PP4Xey9ixXZ2yyPfQebwEWxbvBJXnXsyXv72d/Fi4u24+PJ3Y2GxjTw8lMGP17+C\n1ef/CXrmtdvh5fMW/uv+55BuHsXoW17Bmw/9BeYv6MAD2w7g3d05nPKGDvz5ue/Di88NYHx0En9x\nzinYuutVW+GxAAAgAElEQVQYniu29bFYDJaVxdC+R9C9+D1gAAY734DD3UtxisUQj/O6gDGGp362\nCycsWYiXnvgJpqeBd/ztx7F7w69x+rEDOP7qv7bjVcJiFv7jpYcAvA4AsOH3A3jrI/dj6TlnofuM\nd+HIwBj+6/7fY9HiLqz6xBnGae2Fp9JWkUgkkE6n7etMJoN58+YZvdvb243x0Uns2ZkEAHyo/23o\nnteOX736NLYNbse2we348SXfxm9+P4D40QwOHc3gif3DmMxa2LwziUvOKxSc53Ydw+92Fv7/yb98\nGADw7TXrAQBXfeg0vP64BO4uXl9xwTIsOb7bJcu313+v4FfyBTDG0P5aK/7hvdcBAH7/zGvomRjA\n2w+9jDe/CDT3W+jpXIAnXj2Kl0YzeGk0g+9c+A4AwDO/3IP9e4YwmZnBdZ/rAwDsfrgwz79wbxJA\nEv/+kxfw8tEpdOUmMLlvKwae34qz1j2MfynK+AdrK5oWJvE3Z15iy7du495ifP6Uk3vP8H4czu4D\nEkB6x6noAvBfv9mLC963FE9/8m/wgyWXAQB2HAWAIbS07MGX/uZMaX785oXDeEtx0GXBgk709vLp\n9PCGVwEA13zkrfa93f/zMwDAL0/5hJ2H7R0tUv8ffu5V7BzLYOdYBgCwL5vDX7xhEedmd+5ZtLwB\nyB06BV1dbbYMT/z6GYwPd2Gwaz7G2nrRvDOJ7a/O4Ky+UzAzncOO5w8DSACYAtCECQC9V3bjh/f9\nFpMTWfxi3Q7c8rWLAACPrv85AODKd3/EDnfnyMtcfB/fvB8A8PELlgEAXvvlDky98jKShwewdOUH\nXXET0woAXnptGPGXJjCNNqTiORw9msajKYaLn70HTxfTa8VH32a7f/XQOH64/hX843VmnQIAmDhw\nEKnNzyC1+Rm8cd3DUjfMsvA/Z80HAHR3LcbB+BJ84t5NeOakd2Gi60RMTXcivW4HFuINSL5+D547\n9gJucsRn3y83Yv+eNsz80WIct28xGPIFfwHsfXEE29t24yLrAADguIOH8WLTh5HdM4zDYPjGy4P4\nyb98GOmRvTgwuh2Lf/kqFgP4ZXwML249hNPesguwsnZYTccN4MRjJ+N3m/bjwr8ql7NSGu/7zSSA\ndjxxcA/y+VZcftFpaGstNF9tB5swOrwAGF6AXx3ehbapBDb+ei9eGhCNiDg2Hd2Hf7jwrXYH/qc/\nfgGvvTqM9l8345Kr3w0A+J+ZR5FvzmJw5mX8609jeBfiYIzhE9fz+bP7298CAIyufxKn3fBJbP39\nj/DSjqWIT+UxPJXGr3/+Mg7uHwEA5MGQyP8KJx0aAw79GI8dfwJig1lsHtqMS//sw5gYn0Aq14oU\njkNiehz70/ORm7cPqWKT+goOYvIPH8TPfnsApzQ34ehrY3g9GP5wNIVXHnwEndueR/bga/jF+c1A\nB/CTp9fi9Hk5HHjtPTjw2gQ+c3E3urpabdmfe3QjLv1AOT6ycuzk3x4vtFGfvuTtaG1pwnceK1yf\nffZS5AD8/Kz5gHUA13YzLOwoCD3/8U2YD+D7pwzh+ZGdGM3Mw1RTDL/66S78779fDgD4px9sxSsH\nx7Dgmdfw91eUFdvOPxxG8kgaQDP29xxB2/ZhAEAKwHi6Fc8dzmPFpd34zRMvAwA+fMnby239hcvw\nx6/rxuDBZzE5tguTY7uwCB/CttcXOuvTEzn88cn8NG4mNY2dLxzBzheOACgYDv+2awCfWPcA0gDm\n/a8r0L54MffO3pEDeO7YC1hWVNqvvLQf79n/PA7veh5vXPcw/vXOpwAAQ8cyyIzNAH+kTWJjjJW2\nuMhg6dKl2L9/P8bHx9He3o4tW7bgmmuuMfIrmUwhNTZlXw8NpjE1ncXY+ATnZmSsfD2VLWxRGE9N\nIZksVMaREd69k+GhDFocMg8PZ9DZrB5OK8VvYnKK8yvm8GNoOI18pgnj6bLsJbeZzDQAYHo655Kl\nxOR03nWPc1vcUiN73xW/8bTLzcRUFslkCiybdT+bzCrlcjIyMoHmtibpM937g4NptLXLi9PEFC/P\n+Pik1q9MZtp+nkpPI84AxAAr3gxYM7YsM9M5pZzZbCGtLcZcYXldO+9linmdz0y43PX2dkvfHR3J\n2L9LpScX49NUfM80f0pMD5XzX/Ues+TbelhRFivmzme+PAKMm0HTDEfHgHzRv5KrZDKFqZS7nGYy\n05icmBTeL9czpwyqNE4mU5ieKeY/K8sYt4px0mxETiZTttKemiyUTWf6W/G8S6Ypj/yx37XKskxP\n8+U+7pSpGEYqXWhvJhxt3/BwGqqh5nRmGummQhxLIWUnC21PbmoapSZdzPlkMoV0sY0CgBjK9UKV\nxjIGB9NoaS7HcWx8Es7VTkNDaeTa3LOuM1PTsGKF+1PT5bScmCrkYWZimpNhdHTC5YeI073z9/Bw\nBh1NMaQcaZoaL5e3kZEMOhK8gTGRmdGGNTyURkucH6mUtcFOeZwfUjOJjynGSrs0NPDYY49hcnIS\n/f39uPnmm3H11VeDMYb+/n4sFnoipsyd2YbaUa0pmiisANXN8dIs5tyh0rwMOqc9+yVcBtNeRofI\nCtYwGCntE044AQ888AAAYMWKFfb9vr4+9PX1BQqYP/7OuyDEQMXFSbUWL0VAZ3tkdBQEDJlY1Xpg\n1fE3jJNaUHlOBq0DrDhyEyaVdnZZ3mGWWQzVKOdhLBpzRXMOVsegxKK0T7vaREJR1B1zN9HcMZu7\ncY0qciUU01zVE1UoT5Uqbc1VaITgrVtn10ndnENKZtaUdl0d6BFBwvy0nrPxDZ4Xtc/DmhSbei6b\nFclu8m59qu2q5GiIljazGLgJ0ZDgJAyYdWL7ECjWdVylRGZj2+EsKu3ZCnmOEXY6RiBfdB0HYaNG\ntUWpURj1QiH17RQJkDTMoYyCNHcx7nfAOe2qHOFZmacWt3AurD3UwZ8q3xKn3mtSBeujc1grKSOh\ntE0yvj6yrXawKvTEC/6GNNdVLbgRmhqEF97nkOoLbeJWUhuZ5FdQn7yH8PUShEfFc9qWY2cJY5H9\n3rVLaZt8eKRKbRUXhFfcqrZyt/wzVqNGcBbntOvkJJ2IUk6xcLszJjkhqyDVnOKgDlsY+M0fk8Wh\n/nOm8sU6UW0rKlXazotwFqIxph6XCr6WUBgeFzq11ZjjjoViac+SUq8CdWNpR4Eozb2XZAl94bGZ\n1ja5FVwEQ79qkR+qvc4hh1IdXytIH6PGN8jweKzC7mbM+TPo8HgV0rvCERkrH36D6PIlhEltMZpm\no6SV1iGD0RNPS9vrfR/iKJj7q8dpeFyJSaNSblRDtrQj0DExlSBsUWct7jUp3GHGLSb89SNGZXKE\n0TDORi57hikuzI1APQyCrA6JVr3fmJlY2pUml+x1lncfhhUFIrF6vF4KaJgrtnWYWDmzuk9bYn2G\nqvC0C9GqV26kvtVgPq56VDAFpUnbSlLdmX/B/Ck34IH3aQd6y8vTysqJ5dynXbXh8cqnJMU2ULyW\n5UnFc70mw+OehnZlprb6/bJs1TpuQSQa+7Qj6KdcCdVIaZsUcib8rYBQDL1Z6HeFHuTsZXkEMZjT\nDpI2ER8ym62FmJbLiAl/RCGUohxI91U6uuKNt1L2eCx5HnN2xCLUDkTC0japKH7resWrOSWv1+og\nAZNwrKKbsPcJGg3NS6yKMC1t0SdVDIMGqZJVNpIShemCwDDlhd+XeWKuHz6obMtXnJsTDzovW4U8\n9bC0PUMULO0wyp1oafMVJmjnRLC0DRaixSo1QUMZHvewpKUNvrMjFZ0Rt0gsRItSL6aMrAGfrZCD\nu6oK1Q5a63/lw+O+Ol8BwjCp3qziChBuJriaRa33ETeXK2C2zhayICq/aOJtvUuUdoU2v1HnjEl/\nVgQ3rG/gaeUL7syIxvB4FUqoa5Wjz/fllnZtMsXM2i26CfngAbOh+ep2aPQNgXOEJrwwC/559LYr\n8SdsfOeT1/igbv5TTpCtOM6GrdIZ6WrMaQfOu1BPRLNCKdzug1AcNwKfiKa/lo1iuJS2z6iZLUTT\nK1hxa5r7fVm4ztEP77Y/XqPObCSGx43mz/wHoAnP5H2jW1XCQGnP5kK0WRxGrma1kC9jCKK0w3Gj\n96CqzgtvRHFqIIQCEMkT0aogVHWi6TUcXoWG02fnUBqcZ/rqO+xm2VOb+hKJ4fFaHP7jX2fPnmIy\nqcD2Pu3QP1dk4KTK6aD3nxm6M/e/dC03tKultCuzevyfiOcllGT+U/FK+Xbt57RjXOjhf5ozrDLl\nJ0xAsASZZWTZecukOfQ3JEvbtXpcuuWLj4vfFHbms3p9i75d8NzGLXkeZ+KKfuEdMG4xZnyuK+1q\nU3HnTt5dq9BX49CLMqjDY64fAUMy2LIhecnkVlWohqVdirM87kFsVN+TMb7D8I8+DH/pWkEumExP\n6gpTKElVBf8r3tZUfp87aCVieJXtqkju19IO0HZXupDN1EkYRGJ4vBarx30PP0mc13qftr7tKj2s\nTI0xJn7lS+eWKR2Funpc8Ir7SETMX7mR+i8u+tFY2kFOujJ5peL0MgvE8dvLrViO1JZ2iUDHSzLv\nOW1t0jiCrOh72kHCDuqpAfm8s1xbCEMDMNfoScVeukcwjVaPVzZq4PdwFfkoqQ8PSuE6z6OQ1LeC\nPnCOAsx1S7vGw+P+35+93q6tG7UWQTjyucIwMUKqmjayMWqFy4Bi+FMWVYprxd5W25qvTrzN1Lw6\n7OqfQh1Ya1fmq8OBpZma8INLwSoWdPry0+tapvxq0JTyRqDUgeO35LF3CJ735rzSrvb3tEVlFIaV\nXOt92jqR7X3a4Y3KFa81Q/K1srQ1d7gTtYIqbXEURjc8HiAQszUJvr3l3zcaAfDRSIvzn4xphIxJ\n3zHB5Oxx07QpzXXGhbueMuieBS1TlZ49bjlXj7NQ2hq3Qq3YS88RTKmlXeGuG6elbVJmgiwolXc2\nLO1zMa5z/itf1Y6eK/18BljrL1lx4Rgo7TKVD4/7dStXbhWJIXikbVYrDtQSF/nMwkK0yhMs3LIo\nLUWqEQ5fRS5Aw6ZzEobZVpUBrArlcij9fGhbS6uwAqTCdjUYfleP67r9ijotu2V5LESbpcHY5tkJ\nFr6Hx/3Paeuvg1B579fs/fIcq87qDadiq1ZSy90Wf0isikrS12llIKb3y/nxAX8dDoeFLs5pQ5Pe\nFVraJpZBoBOjfFvz+uFBydueLoINdTPlE/u+Jm7c+gt7tKnyRt35NBAVrh7nOpJV2qcdhoZ1W9b6\n54DE0vYpht992l6db9Nqzh+u4m1px2t05lAkhscrWP2heeI9bKP1WdoZM/cjlGMIjZ5VaGl7Pefi\nwYS/zifB4yueBmWMH6XNDavLGx6pd1XqTVc8XRN6Nz9W9NaRGp5BBPmQg0Fnw9tX/3IY+j9bC9E4\nnR3S0LjWl5AUjHtaQKa0vVx4BeI9PO45AOe9Es11x1Np+xjODbO61vU+bT+rQP2GwTVerntGHkju\nmb1aatB18/5hfU+bMX6NqS7d7PpZQdxkcJa2h0dBV49rLe1iiymfC/M/omE2p13hSInRMLNXS+Z+\nbKeNdjFUaU5bcs/lreCJak6bE1U31OL8GdOErEbfbgQryJV/66BcHiztegJD/wC49mn7KQ8K3CdN\nytslJ5Uu0DLZj+8MV1r/PNZQyS1tkzlt51e+dEo7PK0d7X3ajnjKsk2bDC7lU7ml7SM4D1SFsHS/\npLTV/nvN27ieKSLktRCNC8fRnnv5o5PLde28IRkejykdy+E6ARLn4py27SSszojsHd23GwIEE/b6\nipjL0ga8pYoJP2VKVBJRr/6D78RQB+flXBuon96AWZ9IDXe4ijqRws31ALiMB+GxzGKNsWK2uxPU\nngrUBRPGUc0BEo4ziOQ9erX7KhKJ4fFq9G7dJ/UECoLrJfoa0pR8c9oU2UI09dxzuAvR9JZ2SWuH\n+5UvS1AUuqF2fvW43J0sn5xDjirrQLYAOIhFzM9pm3SWAqSdUXqr46x+g5X/Kl6RljzDpeAmZ4/r\nDW13nGLOZsxoxF5nEXm/L3+xwu9pOwpfGFu+LImxLvSNA+HZXngNM7sk8capDM1Wj8ssaS9LW1b5\n9QvRLMYERd1AlnYVprQrcVpwbyuoSnwJJoV8aL4iL3WhechSjTBF73w0eAYL0bwWlPmZiwozrnzH\nvQbD4wp0K8Wdow6eIQQRwUhb+GzUA4ihDjpgulZYTpyvW5Zl3MnS+SezbCvFbRGbGEeVViJ5PLht\ncsygrdQ89xxoMYiCNrVDbEciYWmbxEiWILqtkX5WRcs9KIXr31ophKdrlD0UJdyWnyh+WN/Tds1R\nadKt3JEJ19LOO4WIiVaGaKl549rSBT7vXHv47ZENw9rsGb7T0pYjTgn4xmhOW3nhwj08rrP2YsJf\naOKgNs24k+7kTiQhO9+PcX9N0c9p+/LK8aL397R1w6fMcHjcWByJpR2G4lAt4ixfSyxtIeJ+09j3\n7gAvS1pazSVyc50C7zZPJ2WYdk5dL0TTD0eYulS9L8t4Px4Y3uP8LzWcBh6GPKepDEb2O2xLWxcm\n4++YDI97TU2792mr3wsSWXXWcGdw1hjm+NeweVHJGJO9b7gQjes8eb+vgjHHNIrf9SoBn+nlCfhi\nEWdHL4/KlbZnelbp05wyuYP1SZ3lUFG2FNa1VzU2N77MjCsj93NieNxAaXseTl9NS7tIUEu7kkwq\nW9plP5SrNitcpGExcfW4ztKWu3E+C0JeXD3OlQ1xkE9eUZ1YkvUETLMSVLdPO8hJV6pyou04+A3D\nZM2ELD5ShQu7HJXXLXgrDs4CUg4pCO8YrNbR7tOWLVwPcctX8OFxSVoL19qoO8tZhaerAYXRI1fN\nMVhr4YXXanHpPu0AZd3pjeorX6o6LZ/Ss7TPpcPjXIXVj94BemUapo01iyeiSW043g2XcTI/dP77\ncOwRdvleOEpb3cTwq8d1SRTakaqe1r/zt9qyqWyhhS6jgzTP+gZUebiK855kXYEpJpZ27Q1tf9aF\nbiGaFNPhcYWnpsPjMm98TxFpAqjGqF8J7fiGuBCt0uFxeFjBASMaaASz0ql1lWHiw9jznI1lpY67\nc/GZx5SHjzScG0rbxNL2VCjmlc9v51WmoCo9XMXz/eJj2T5t1XdrK53T9lpIIs0nr/Fnn7i3YDni\nLdYbg33aeemJbd5pye/lth+qxFai6mw6rVze2wCJF/B72ipL257TtmUpKG1le8ksmLXG6gIVaJ92\nkcIIUXFOm3Ne6dnjAQuydOuB8FPjtbNMMmZV3KtjiLksbScx1YiLl7+a9qFw7S6XQax6flRAZWnL\n2wKv5inQSKGyPXfu0w6+jskPkRgeLxuWcgtI6YXmsa/5Bg0xmaAmVLTaxW35yV34E0nrj32tTjfZ\nryBiiG5dKptr7DSWtjJQaVe6/EvRSZDq0SBK2yQ1Ku16B35db7UEXrdhqgMMpvV1Isinq0IcHg+M\niaWtMzL42lXxSJpXvzqope0xPC7zNdCHNLhmVzWnLT+PQT48biqDw51idboK7UhKiIXO8+xxxhi+\n9KUvYdeuXWhtbcVXv/pVnHjiifbztWvX4qGHHkJPTw8A4Pbbb8dJJ53kGbCsl+RSFwqLReaH+xl/\n7d/Sdt/ztU87UC7xK3ilQ9Ou63D3aYuZYElkkM3zVjJUZFnqjNYN8an3abt7vCb7tJ3eWYwhjlgg\ny8tkeNzyKtyeYfid09Zb2hAtbVa0tDWT1dwoj3KuXOx0yeV2+mWS5JVMx1TUn1a+Z9Co65xww+PF\ns8crqNpM6PwWb5ZlCdgp8La03f4GOVufs7QVrzvbIXlbKdfkutXlTlm9jjEN48uRQfBU2k8++SRm\nZmbwwAMPYNu2bbjjjjtwzz332M+3b9+OO++8E8uWLQssRNmo8ZcIWku0QkO7kjlNuQAGXjHuj3y4\n1tgzM8z2WYrPAsRN45RXQGLPXdMxU92X6QVBIUvfk+ZZEEvb4P5sWdoeioAbHheG/4Lg53vtZTfe\njirs8+h8Dvia/j3Pqu9SOpWVD9nwuNPHwGnmyk55fZUtBlV4YYBqZ4JcDqn/HoHas6GOVOM6NnLz\nXZBSZ0SGp+A9lfbWrVvxvve9DwDwtre9DS+++CL3fPv27bj33nuRTCbR19eHa6+91ihg+QdD+Ih5\n9WT0lrbKMvVHzFWZDNG49VqIJj8RjXcp/Z52CFahbjVoWAvRvL51zsVbsOq5XrsiSFm50e7T1o1s\nBEpTR29dsRu58jltA8UmmdphiEmDizF+lKdoZiuJic81C9H4r7g500bxhjbcshs7bUP8ylfgttVg\nnlRncHL1rPIp7UKXyyv/AuCqWwol7uyIBxke5zplig+GMMVxxdLDqRy/Lemwa6k9ZeVAdAdllMIx\n7DzWdHg8nU6ju7u7/EJzMyzLQjxemA6/6KKLcPnllyORSOD666/Hhg0bsHz5cq2fvb3dSB5K2dfz\n5nWgt7cbXaNtnJtEot31bmdnG3p7C/J0d49w7p3Mm9/B3eue1+FyI6O1tQm9vd2I2eWhnNolP7vG\nM65w29taAABNTXH09nYj22phj+B3U5M7W2UyLVzYid753WBNTfa9np4uzE+U0ycxUvpd9rMpHkNv\nbzd2S+PVLA2LNTVxhS2RaOfctaSmHTIk0LuwA5lMJ/YL/sxf0KlM37b9fDETw5iITXHPOzpa7Odt\n7S2QNWG9vd3IjE277gPAwp5ONBW/kxeLFdy2TJX9WLCgk3Pfs6gLvV3dmMyX3Sw6rhttLU1Itjah\nlNuy+MnudXS0ue6JRuuChbwMrS1NRuWzRHxeOw5pZACAbKvE2vFYgLRoURfmt3cjM9aJAQb1cHrM\n/scBc/xbkKsp146k4Ka9owXt7S1KtdHT04UFPeX04eJXDHLRcQm0tjY5bxkzf0EnehcmAADNLQU/\nWtsK9YPbxuYQsFmSP856tmBBJ46NcI/R3MwvGZLJ2dHZit7ebnR0lOtIe2czkJanTnt7C9od7QIA\ntBTjEHd8G1LsmPT2dqO9rdl1T/ZbxwKhnnd0tnLP5xfbyFx2xtH+MbQ0NyOGHIByGwkAzc0F2dvb\nWjh/DySGHb7Kc7jHUYec9amkT9hEG0pZUmqfAaC7u90V3/HhUhsk71TOm9eOHuGd7ilRP5XfFf3v\n7nbrsqB4Ku1EIoFMpqyknAobAK666iokEoUKsHz5cuzYscNTaSeTKYyNTdrXY2OTSCZTSKWmODfO\n61ICTkxMI5ksKPzxcd69k7HRSe5eKQwvZmbySCZTGBnOuJ6Njk0g2ZxCJlNWFiU/p6ezAIB83kIy\nmUJuPO16P593N6IymYaHM2ifSWFotJxGg0NpzEzO2Nep9JTrvbzFlHHMZvPSZ84wACCVEtItUw5z\naCgN5HKYlqXNyIQy7JmZnBDGFOd2eLicVrEYw8TEjP18airLnz1c7LImkymMjk5IwxsaStsryFnR\n7dh0OYyhkZTLfWyiFcOOeCWTKbS1NGFmOsfdc9Lb2y2Nc2ZC3plwDlcODfPlQ5U/KtJj5bir3nOW\nQediO52SGxxKY6Y1humMO4+diLuAZbo9mUwhMz7pmlqfmspiairLuY07MnloKI1sPg9Alsbl/J+2\ny5U/tT06MoFkruBPLlsIR1fvS+50+SN7L5fL8zckeniyWNYnMuX0mJiYAdDqdgxgeiqLqSa+HckW\n42Bxc7y8LZ1MpjDpSPMYyuVGVY5ljIxkkEyWZXO2hQAwOjqBZFMKuZny/ViMIZct16NSG1n6DRTa\nT6cMXNuvGEkZHCyXb2f62/rE0UZOTs4A6AAAjEt0wVipPilWuo6NTSAvvjM+qVzZLvo/Ps63s5Xg\nuXr8He94BzZs2AAAeP7553Hqqafaz9LpNFasWIHJyUkwxrB582acdtppRgHLF6L5G9LWDZ9XOjxe\ncs2fwOVjm02Q5f+lIcrSvJDzkWIYihsaDDAG45JSM1xeXrDlMWnsgegyr5HblccG0ZVvt9McrqIZ\nHg90uAo3PK505Ntf7nXfchXj47UQzTU8rlOIBgvRwC+IisUU9cJ4iLskp1QKI7Q1M/DwuMExptrX\n+eFxP4KoBnvFT3PyVmSwiLqm6RTtkiUOj/ttf90zOwW/uLDlB6ZUshDNWTpiHgfeuPzRRLGmw+Pn\nnXceNm3ahEsvvRQAcMcdd+Cxxx7D5OQk+vv7ceONN2L16tVoa2vDmWeeibPPPjuwMG6l7flCkEdm\nsmg7BCYeqB95NTLS+VTFsk1ujsdALElgukteBHNvKkAoBWqdrZTIcy2EYvEM106w8lO/iL6X8bdC\n2k8ocicyNx5bvpTLeyR+B4iDPDXcXqsoKxvH6vUQG8RQ/fLhNzfnWtDa4QcvL+D+/PAMg3F/g+Ps\nYCgWoik6rp51V3Jbvs7MQAGpehe+/TLHU2nHYjF8+ctf5u6dfPLJ9u+VK1di5cqVvgOWr0NTW8e+\nT0TTLKgyE7AYrrOHJmncTcM3e1G9EE0sn5a05BkIJr6i6CmXn0t6qNJep58wxbxxLlSC0EEWPXbK\no/Dfo1dscoypqkyaoJRLIU+gD/H6lqsYR6/DVZwdRlZ4IncvPFPu+GJ8Pebqk9NDvSUkYrFyvvn/\nYIiuQx6wcTWoEzrr1nWSb9xcDlVT47a0nbIEw2sEs1yXylMDsRgThg1NwnFeqJS24hhTiSd83ZaM\nitiPHZY2H5hERj4itTrG1FNpVw1JgoqKyLRX5+F9IKQVwR7GD+iB5zulP97hyPZpB9pT7HlD9kja\nLQ0cpmv1OBcS/4yfrpD7L+vQ8JaM2FNxd5LKmRFibXOeiFapJWVWCCUy+AiA6dybfvjRPa4hI9xt\nW3pCzFFjXxk8+maOcmbZUxOVSeMtUTB/+Wu5P/z0RYDSzvVp5aVDtbhb3g4y6c/yLfdNroOp6BgF\n6W9XSt0dY+qxCr/8zPVeMMtEdgKTkXL0KjhSxMNVNNahbRF5BOmB+EEO95yVu+DKLVkfYQrXlmhm\nOB1U4qEAACAASURBVK8s9ZC4Kh9kH9Pge+IK64BT7PaL0jB0qD7N6Uzpyj+nHfAYU4/hcXv7m5d7\nxzsAtMqdn9M2kFST5qX3ufrse8tXsLC1fhq1CZpHFX2aU14fdVOuoVnaQjG0y4/lsLRlwvgKx3vL\nl/wYU7milhtkpV/OCOnPIXctxtTkWZj7tGfxgyGSe7oJVcmbviqfzzSTZbzPYucvQK0MEmUacGhQ\nF0bxjvJSe+BMBW0MdylsIxa9NWr0Pe65jzEt3Ze8F8boRUA3VQnEYyGal0XCOw9ylJti+JFzY+KP\nIIofP3TKM2jGGLyoTy1nh9FvY6W6FfN0Z+Q9MywTDgfclEgsgNJ2/FZZ2spOu0QGnQHkCtC+J9/+\nx/lT4VcWgzCLZ4+7E1G0gJwffnAeqiC+5+G99DoIvobHNR8QUDdWseKrPixtj3teeJ2IJl3BXmHY\n7sNN9J/mVPqjeCY7PEFraUuGx/lPVPpDYSTw/XLNKXBmgfi0tEtxVM5plyRxxJu5P+/o8FDwQDmp\nbTTHbzLt4XyLP9jC34oX48FaP+2x3nxTDqWWTxAr37Msf0pO1UmV1AL7l5/V406XrtXiinbJaQUX\n4u3P+FFZ2k4sS7V6XBKSh9KWGSTcB0Bk74iWdoC2KgiRGB4vtxPeETMZVi94ybTX3uEU+6qSxt7M\nEAjY1Xc8lySRQz73gyDFwrvjLKls0gIcPEzdrCcTmh/+PGDTEPgwVHPoTJaYVbK0K/5kstH7uhIk\nOhW3UpV+qJW26nvarukM7/ZXeENbse0/QdcFaGumrGEy8jSoLO6yZwlKLig6w8XXXCxTXri2z8nb\nCH+dEDEY9YlozgsfDaHmuWKbtrrNM4xWiDo7Gt/TVlmwflcB8+70157ySTvOjPurpYJWWX6MqVBZ\niiU2Bv0QjmdYnj1nyTNZfvmxDASn/AdDhPBdn+b0DtLzGFOVdaCJq+w9FcpdD6EuRAuoJDzmF8r7\n2YvK1t/2aYl/QoOtGiEw3LtYEp9XQH5Xjwd7pvfTuyMrtbRty9SptP1t+ZLJXFhdr9unbY6zPrma\nNbH9KI0ciKvHfXaGuDZB4Zwp9mmX5ZXXXfn3tt2dp5hH3RfXA2nTdC4obVkk3NsHykhm3PTpULHS\ndme8v5QPYKEJ1o50G4N9LUtA30F6Do/LOtnyAuw/bJUM0FQwk0+lSvdvOhselWXv1U4EiaTKT1nh\n9uOtz1Epr8Lhql+lvK50XbfQXps0cu5Onfd2m2rg6+tUHvJUueXwvK37gIffoLTtkqQDHGMIUHec\nHQzV6nEPP1XWt8agUk8byO4zTrFrF6KFqLUjMTwunYOAuncj80PE/RGKYInGz7W5C6QSqSNDS01q\nafNuynOuMcc9I+/5sNy1Trh09rIZ54ZbDe0jbLH5cA9XO5+J1qp3/luSM694S1tciGa55LCjaunn\ntaThK905V49XWImNlLa786O0tIUPhkg/7SiGzw2Py4cw4fJH3uDrPswjt4zK7qSNukZ2/QdDAuaL\nwXuyOU+ppZ1n/hp5WR9VWNBpCT76mtMWLE7dwIn0gyHw/wUUfuhbNTzubpsB9w4I8bm0ZZbooBi3\n0ECvi7yYE8PjfOUt/HZZQF4WnZ9E8yecfHjcj18a2YxPRNPck7kJejShbmRS3vcI3iGRuRVVqC6P\nS4272Hh4icd56R7Tk/ghvSkPUIPmcE+9AJ74U9rltzxLX/kv0yxEi5nv0+ZeCzQC5Xy/5Ksq/byl\n0koQuHH1KHRQLZdjLqd+92krOzXi1FkI/RHmCk88wdCttAvX/gLnO3TeJ6J5j77LFbi5fLK6JKzr\n0PcUQyMSlna5nRAac4/V4zpjpWJLu9SLl1krRq/rKrFK2/AnosksP9H/WEBrt4Tus5jOcLhn9o9g\nYbvDEE5E07gtae3CfXmgXmshlAvRpHFVn1muQlUunVZuxT1vA0tdqthiMcUiJH7ngi2gdg7c+0Q0\nMP5ENKUFya1VEPLHuRK57K3DQZhz2s4egg+jQFbmhN+61ePubVU+CogibKeXeYu3tf18LpMf8WPa\nNrhkePF5xjhPTII22afN12Pm/qUaHtf2x+VpJM9fsUOqJszpnFk8Ec3xs9ROiJa29DXDyIuKwVwy\nTib+nt2S+w7fnzD6BqBwzVwPQikXmqFqV/6YbJqWBeEKUqxRfK9YPg+q20qkr2DKw1WYxD3jXZoQ\n5vxVZWFIpglUhrM8EHM0Z7ZwysLhr+nBFBbE8sH7W+lZBU6C51xJrmAdCLci1PUsPG9IXjFXMLp3\nGa9/C/dksxNinfaZsNxiToPvaYsdCztc+xeTPNejOtTJGSb3BUIjXytnFi1tdyK6GnPZEbHSzJG4\n07znD3fGG/kV6EtYvLVjSdKoRPmZ4fCMAm9L2/3b68QxL0SX4tnj7gbMXR2cc5oinvu0FaMwMkub\nuYXxBTcS4vjNpXuQ2m6Q3vLtSx5z2uVhr+J/mn3aih6Afk7bUtR39SgEv7ug5Maj26JJH202Buik\ncZ5qXtFZ2k55LQat/CaHUBW8cJQ3r2XfGlzD48KoFVfKSnkrWNp+T/DjT4jzPlxF1k7xFo38PfFe\nDIo1LNI05jsjutGL0BYCIuInokm3zsjzweXepGBr5bMz0R2eVHbPGwFUqqzNNXQeKAzZY+2wVsDh\nXo+sYa5nznLg3Th6BWmpTkSTuQ/Q21O/olqIFiDnjJS25F5M38lTHXTkxhIWCCql4Bs2aahyGWTX\n3Je9AvfEfSjESjDJI9koj+Y1BllbI1cofJlngjvznqJreFwUiJs6d9emGH9pFqbjt2rHAa8HmeSX\n87GuHZO/FPPIFPeKD3WahliqZm94XG5pi5XV/Z5zqEy3RczVW/dbaJj4wymn7AXxffWXZNTtm7CC\n1+GnaBHb+7SFoSPfCz5ccqvzoLx3uRA2Uzn0ChNiXNQKTBUf3UI0USkXfHUMpSnCkMfVedOst6yU\nWfE7ED4UgjNEXmm7f3PTAkx3IhqgWj0uyumqpUJHrBC6ujcuKx/8Pb+rxzXPOG995JJdJxRWocI7\neV1XSylXOLLOV0w9sgOfq8e5Ntfdkefb3WL+uOa0LX1ZEsN0Dr4ZfeXLLYPS4NBY2sJN7XPGmJCn\nmjQNz9COyD5tu31UK2HZTVc6qut9ANTdMXl7GV5fSnpCl+hG8civFGbfYBKe2T9i7mdGYeplEPPR\nOcRc2jtbaNf8NGzO52KcJekt7bRJg9MFxeNUbJVadB4NSuG+eUthr8oXLSVlO6tugF1PhIatdOCK\nZuCdv1Z2gpiHJGb+hwETkk2GWp1LFKMyHMPtYEJg4jB24IVoYK5r6flOQkfLKzjXY0WnTLUQTa6g\n5d1kvWFTfhYzOODFCTeG5TKywit1kVg9rlRA0lWj6sQXM1G3YthUPuk+bUkOuryXmvYeNVuYV9Sl\nkapBDmJpx8QbDmQVoxR20NXQ7lEQR0WJ8Y2Sa9WmwwJWlhtmuQLRLkTTjGzwq3oNLW3nhcJaq7QS\nM4MKFMjSdmofJrrhfTMdHnfmYflwLEGJaFaPy0ZOPGcXNMmrHX42zBf39Jvb0hZ90u3TdkbIdYoc\n517ySOGv864sDU1xWdZC2vOXxXQQPwzDmH4ngqtNKP9WLkRTtO/SBcNOt5qvFDrb+5hH3WfMUlva\nmhHLSonYPm3vguVVV+1nYm/VZ6JJR0t8Wb7Bc0k6vKO0iIVlP77jaV64tNaE5kUvkdwNIH8pG8rT\nWiOSnOLPdpafmOw59Giats4Gx/Gj4tPFFGGo3cgaER8yaNK4sGDHYE+0bDjMs0Cob5TXtnhZnCG0\nkrrhcVWZVTRQ6rLqfq4dHpc9kiltxMCEPcRMVi4NcLejjLvkOoKlPrW4XZf5zS1n4nnv0+bWj8k8\nFXsartDc93hjRhK+1r1YPsLT2pGwtFUGqPQMac1wpfvkHlnvy1hCAIKlrWnc3XtLDXJZGbI7HJV1\nKm458K+09dfSHqxEARiOOBUv+Ru6Oe3CM8nwOFMHWvhKkhCmpizI9spKv6dtmLjqfQPO9KqwEjs7\nHoqVqeJhRUChgbVHrdwCSr5uZmYdxXyciCYd4uWSWSgDHt9vl2350upbbefbtJLKKw6vJ3g3+rPH\nRa8USlto16SywN1xFdPQ34loDn8Y46xgsfNs79N2nD0OZhUt79JIoiwMdfvJTY8p3MjPtJDrCv3G\nHke4XKYovqet6ghpRiwrZRYt7TJ2O6GpyKLbwm+mfia5DiITJIXBzK8gmaTeGajuoIgNlr9wtT1c\nlXeSBkofhjqfXM9j+nhzb6meMUmYBvIw2U1em6iEESWT3BVHRByNg69PLrnFMnMja8CcjSGvynVD\ntGVXzoVoBvI45WDC/ntd30DWDkCeRVDf0vpnl2lFmnl54ndLUzkExv0t+KXtcZje5G5bkjphjK4c\nMCY1SZ1vxMQtAx7tusuJwfA43++XZrD2sewef7CW7B0myCOXTelBQKKxT1tScAG+c1Oe03b6IfrJ\n/w5gJJVlYHy4OjkL/pv1fj0CByCfT1VZhyL+V8l7KFSh4nPvxMwsR/GJ61r4trRY+GWzsEzTCFmS\nsXO/lrb9U3G+sRbOWdm6kA0jBoX/wpH3+gY7jrGYVkO69mlr3JrpafewbEEudW1wj/YIK5Edfqgk\n0VvaWnF5YZUyivWmJJc8j5lCJvmomnpjnGVsaQP86nH+/G8/w+PcKDRjvFUrSCrdp11sjHXTQ9pt\nfkb7tMu/y1vmhQyw3arPmVDuYlCdiMalqVrhVPytAQeRHh73PtlK457xxcm0l+lebGboh1sTGYXH\nI1o7Ou+Z64FuP6TJNiTZDXealv+GVQzd5TmIcnTed1c8bk7bfdq5yz8mPhP91cAfpCBXLEY7JUxR\nvCxdI+LsaHGdiNJCNNFPzRpvl4KLud9gkvUFnpFVN+DlBXPhlT+nTObWqFdX1H1HlpKqUQRtyN5B\nu26Htv9c0NKipW0/Ejri4Ny5la1pGVHt0/basmVxnW+5366wvFY7MvXQvUtpmwVpxCxu+XJnnPuD\nIdrXtFaixeCae3H7pe4U6NwbDa9oS4a+nyvfZ8j7V0orcT7P/yp5XQMppKnrpqGlLTxyLwhXnz3u\nik+snDaqIMUvGonyqeJsycpWgOEaZ+h2NyzG75utuBJ7FW5BjnJHMAZd+XOfiKbDe3i8sEXI6UvR\n4hNkVg0zAoq1LZzmkO3T1pVHiX/Suq2JvevQ7dKqaWceu3rA7nDt1dYOr9yCOOSU+COdbxXbQ3V6\ne8GP2Ajpx/jkL+/TLruJF5W23tIWrk32aXNKWXZfkYbaExOddVdQKAIWmLAjQF0nw+ozARE7EU3X\nw5YNj7ve1vRAdQpa+gorhevdg5bdlSsxj5xjgitNuxxka5cqTF0FFiuoSzAv/+WutVdcNZDr7OJ9\nVaWUPXJUcGX4TsUuE85QaUsbZ4mMFWHQmZDdVuzTcm+pdHfMvD03c8ogX99v/3KVdbc7Vf6rGniN\nOOqHvhb76zOdeXjnbvtVSttsn3bBhXN4XFLPTDuhQjsqispfyxWml8y68qDap8195Utai52OneVL\nl1cOneNZ98WOUEz5LOiaBxmRGB4vGzVCRD2+rKL6WlPJL91XsmThFe5Z/DPOj9JeTAOlpSkY6uLL\n79OWWn62VyXP+AKtHAZX3Hdbveoeorin1PXpPwVeB7i4LCmtPirni1pXiU2Ju2xwz5jb2imnr2Pu\n2HBeiu/olPKUn0tmEIcPjbx2BiIPz4FsVX5hSFyyitd1RgCDbGjb6VtMcSKae5jQ0RCilG+iz+pR\nG25OO+YYHmeS8AyQNr/SNkhXqIUplpKpp3ld+z1tbo2CJljZc51VXsRieU6gGNOHI5Ox5CdXvC3G\nW9qSelP+ypdu9Th/zZ3doFiIxrXvMkvbaYl76gL3vZjHehbG1Ceiuc+XcPsflEgMj5d+6/pa0pti\nD5DJnSn98olsHlnpfwXjIV7DOzr/3erKyye9QpUlt10xZZVVFoJX+eXHwlxh8quMHWVFpbQ9tggp\nF27J3CsqpTGKfg2XJjF9+slgkvojceS+xckhnWXlfyoWopkoSnuxWMWVUdL5UJb0Ar4X5OuG0006\n6bZciqFcRcdMtp7C8xhThdIWV1K4Us1xIwYzi12UpNBpEBtabny88Ee6T9vHFgEuUEWaKpRq+aei\npGtHXZUKRvOOJAzhWYjr0KJmaas/nF5KEFWPSSxMJvu0ZftYxZXb0n3akvh4jOmp3+NeUVvarik0\nWz5xgROkK35VbZLqwLHytSQN7U6W0zqS+19wpw8jr+2bqBsv9eI6txXOL0QTOyp8njtlZHwXXi2o\nIFuJ8tJCYS7Zu7iYB6JJB1eADkubn5IuXJTrF4PLkds3pwcKEYW8YMUVB4LMspMHS3B7jPki6A7c\nYHhcthxAZmk7t+K5rVuh4+fRqDMoukh+DJaie/dcebG9EHSn7itfMbnIyjC58MV21xFueZ+2xNLm\n19lyuNohg33a6mNMS2nqvOfRVknueX4wRLS0nb8r+KqaF9HYpy38Fe+r3+NLD5d+wrXc/DTMvdKT\n0mITU6/MvZZ7JkaI80riBnpFZiSTKxzJI/tH0C1M+pwWF75I3eqGbr0aUDGOdkdN4l6TByqkW0pi\ngleuOBt5LX/Bh9LmF8NpZ1lddUr23EamtOFO1/KeXdFnXYdG0oArZJMsUJYgSxcvRKHEx5IyZNCg\nuaqU/VuVp8pH+q18TNKRNVXawjsme8pdI0FM1W0xEUAxp60x4FxeeMjsqQoURliM7/kqCXFKO2Jf\n+RIbU0lEVYtmLbEHKFx7rkAV7pWcx6RyuuVyyy45L7n0Q7m/uWjtSCw/18GbEsVZ6j37sbTFfdB6\nS1vtmX71uLojAAjHigrjmqJ8zgWJurlcbTx8WNqQ5H8QREubL9vmi4Icwsh/c04k97k5bbeFyo2m\nMEDXEsWlOSPJI0FWe9icl9bpRBmPkvVrOS1OafS13Q13WLq6Lbvv6viVMlTR+YjJu0jS4XHEoJLf\n2NIGn7+WuDpTU3/cYQrhi9eS1eOcpV1smHSh6XaxqKY6PE9EU9QRXadedfa48h3Hbb4dFUeNpVEI\nRDRWjzPpXS6mMck9XUUqVDR9Skk7V3abJdXMxfcC5oD9fhmuMyFGX3mhbpJkzaEjaKl7LZx4fANj\nepa2R3unvxaHO5nity5AaNLc8cyzTBpmO1PJyC1E03dkvMNQd0J0ghSSsySHpHMHA38d/tlohsdl\nzuw+gYHfymps35coSo3osvwpt/OKzPMsxN7h8kOpjPvrWv1skqUlf+285cuXy9gVBDYtc64tX4I8\nfIrx8Sr9LqlEdRjCNXelWj0uf1+2JsizjytLV1l+QXfPR3mpgEgsRLPncIWYSo/vlmZO4bfrTFzH\ntdwv+XmyznDkX/nSy1W4Id876bonschlIw+qnqj77HEGFndnq+mqcvfKXeez0l+rHKCH/4Ck/Lqu\nZYeRlMJnwkq0ctroRg/cDb8lPBfd8/fttQ1cAph+5cstGBM2D+qsSbNAPAo3hPpkl92yHM5GPiau\nHi9aEfoT0fRDg6V84KMmXz2uG+3hDquxy2DZER+0d0eSzwcm/JV75VVP5F/54ts4STEu55Erq9T1\nVVV2LE3/Vj6nbVbmXMpfNDQk6SR+TzvGmLYsudcxOGT1eSJa+adbxwAw3qftpend+7SdHlout2Hh\nqbQZY/jiF7+ISy+9FFdeeSUOHDjAPV+/fj1WrVqFSy+9FA8++KBxwNIOrStessZP8VtSkMWepbfv\nHv3PUkaYaG1ZJ0HqpbuVkM0FuhUd497hxJDOL8rRdRbFG+5yHZO4kgWie2iUdDa6FZqcf7psVJyI\nJq9/mkxQ+a90x6kmZTCGoRi8zNy/lJ/mLLqTNV68A83b/P2S7hHroXx4XBpE6Q2X5wULz63w/M6a\nStp32VNJtfCsOO7HmjB0dd0dtKCAHR0yzg9nvw7MlXem/US+z8VnqBgtWYqJHQRZu61NX4M5bXGU\nVfTCZIGiiNkOBGfH1xke76qmx5g++eSTmJmZwQMPPIDPfvazuOOOO+xnuVwOa9aswdq1a3H//ffj\nRz/6EYaHh40ClqW3qDC5b6rabtW9K69rtwzqHpd0iMXenyzzS3xfZkGXfHTObbprq8yi9+zpO+7L\n57QV7j3cyU8Jc1sVeiXpFYYzo/nGxS12OX+UnVxJ18v3QrTSRRBLW1auYqKlzRdu33Xa5EQ0Sc+Y\n26etWTxWapzFKRB7KFZ8orG0RaVasiQ1apu74vbsin4rXtV/VdNdFsT66wyMifedL9iXkr3+Ticx\nfqwlJoTnLhvyCFjMFXT5WVzs0jrntAWlrfFHhB+BEqf0+PamXG/cq8eZM0EF97rkVa4ed85pe5yI\nxu0eMf2eNtf2y0dElfu0/VgiPvFciLZ161a8733vAwC87W1vw4svvmg/27NnD5YsWYJEIgEAOOOM\nM7Blyxacf/75Sv/Wr/81xscmcWggZ987eHAAmzYewWhmEJ0zTQCATRufRnoojTa0AgBaMYMuzCAz\nNINNG9MAgKFDaXRhCgDw9DObkctZ6ELh2ZEDe/C73CH7+uiBV7Fp41FOlmkri95kN58gowybpp/G\n0OAUrNYU8s15ZONthfB2H8Om5NMYmQJail+e27QxCQAYnRyB1ZpDJpfCpo1PI3fsKLK9x3N+N+Uz\n6MIM2psn8Vrx2eFnn7VlbJ8EupLdOJDdB+vVLPYnZ+xnO3e8iNThVtuvzNgoOq0mtE032fc6MImt\nW7bA6nmd/V4JKz2FTRufduXHoeEM5reOwiqm85GhKWzaOGI/Hxguy/Dyzh2YHnwF2QOHkOs9HtNN\n7bBaxwAArx7YjYncQZf/ADAxAbQ4Cm1yPIlNh161r4+N7kcnK8SjJTeN9NAANm1MAQCm0im0oAtA\noRGaYhas1jH8buvvMD7o1hRW6xgGXh1H3uoC0ASL5bFp49M4NjWI3qFCXg9OHUVvqpzve7e9iqnO\nDPYenbbjuuMPf8Dga82Ylwfixfx/7rnnkZ+/z36vvbMNUxPT7vgOTduLtJqRQxdmEG/J49iC+XZ6\n7ds7id5sQQZregad2f3YtDHn8ktF+8ARtBfl2vLbrbA6O11uUukkOrPFdJ2aQEtTEsmOeWDZeDE9\ny3321pkmdGaasPsPuzDeMgrkD6N9sg1WrEnwtZCRTc3T6JxqwsFkOwAgm4fdeDcjjy5M4JlnnkWr\nlcTEeBs6i7ICeYxODOCV16bRMtmMzqKy68qXP+W486XdODL4SiGenW04nBwG0FKQM1+o4zte/AOa\nZ4DebDeyM2UZS41oLJd11YES+196CSh4j4mJLAAgnU5j08ZCve+MxcHiMeSnLLQgDbQwZGYyXP1h\n01OYctTvQ4cH0Zptx3QxvRksZPMTttwtTSnEp1qRjRcE7JqIw2oGMtNj2JR6GuPZCQCF8pCayaFj\nYgqM8WkfRxozo1MYi3cVry0c33oMxzo7kOk9HtNdbejMFMpj0yQwMVluD3f/YTvGBycALCymUxM2\nb34Wbc0xZTkuMZLJ2f4c3PsyrLF99vVMcgqHFyxGZ6aYDq8cwKaBpzF67BDaiunT1NKMdHcc0x0z\nsFrHkI3F7LTszg4g3ppD8+QoNm0sy5AcTQFoK8paruetyNr3Dx0estN3//69OL610L6PHp4stPUz\no4jlCuVzMpMGismZHBx0tYVHDk3Aas0Auby96jeWY3bbv+/AYewS3hkZH0TnhEN/tDTZ7n/3261g\nTTnE8oXwM5PyshiEGPOY2PjCF76A888/31bc55xzDp588knE43Fs3boV//Ef/4Gvf/3rAIC77roL\nr3/967Fq1Sqlf7d/9iehCU8QRG1psrLIx1tmWwyCqCviVhZf+D9/FYpfnpZ2IpFAJpOxry3LQry4\n0CmRSCCdLvcgMpkM5s2bp/Xvtn/5UFBZCYIgCKKh8ZzTfsc73oENGzYAAJ5//nmceuqp9rOlS5di\n//79GB8fx8zMDLZs2YLTTz+9etISBEEQRAPjOTzOGMOXvvQl7Nq1CwBwxx13YPv27ZicnER/fz+e\neuopfOtb3wJjDKtWrcJll11WE8EJgiAIotHwVNoEQRAEQUSDSJw9ThAEQRCEN6S0CYIgCKJOIKVN\nEARBEHUCKW2CIAiCqBNq9mlO5yr01tZWfPWrX8WJJ55Yq+DnFLlcDrfccgsGBgaQzWZx3XXX4ZRT\nTsFNN92EeDyON73pTfjiF78IAPjxj3+MH/3oR2hpacF1112Hvr6+2RW+zhgaGsJHP/pRfO9730NT\nUxOlcRW47777sH79emSzWXz84x/Hu971LkrnEMnlcvj85z+PgYEBNDc34ytf+QqV5RDZtm0b/vmf\n/xn3338/XnvtNeN0nZ6exuc+9zkMDQ0hkUhgzZo1WLhwoXeArEY88cQT7KabbmKMMfb888+zT37y\nk7UKes7x8MMPs3/6p39ijDE2NjbG+vr62HXXXce2bNnCGGPstttuY7/4xS9YMplkK1asYNlslqVS\nKbZixQo2MzMzm6LXFdlsll1//fXs/PPPZ6+++iqlcRV49tln2XXXXccYYyyTybBvfvOblM4h8+ST\nT7K/+7u/Y4wxtmnTJvbpT3+a0jgkvvOd77AVK1awSy65hDHGfKXr9773PfbNb36TMcbYT3/6U/aP\n//iPRmHWbHhcd4Y54Y8LLrgAN9xwAwAgn8+jqakJO3bswDvf+U4AwNlnn42nn34aL7zwAs444ww0\nNzcjkUjgpJNOsvfbE9587Wtfw2WXXYbFixeDMUZpXAU2btyIU089FX/7t3+LT37yk+jr66N0DpmT\nTjoJ+XwejDGkUik0NzdTGofEkiVLcPfdd9vX27dvN0rXnTt3YuvWrTj77LNtt88884xRmDVT2ul0\nGt3d5cPVm5ub+e/kEsZ0dHSgs7MT6XQaN9xwAz7zmc9wX5Xp6upCOp1GJpPh0ryzsxOpVGo2bXrw\nFgAAAmhJREFURK47HnnkESxatAhnnXWW9BvBlMbhMDIyghdffBF33XUXvvSlL+Hv//7vKZ1Dpqur\nCwcPHsQHP/hB3HbbbVi9ejW1FyFx3nnnoamp/GEX03Qt3S99bKvk1oSazWn///bunkVxKIwC8AmI\nVqK9lYigEQuxsRFsBMFKsNCoWNjaKEpAx85KsVSwE7T2P6ik8eMH2NhKQBjQQjCCWwybGZZdHWcz\nuwTO06a4l1PkcBPy5t4Mc3refr9HuVxGPp9HMplEt9vVr/2cAf+V2fD0ZjqdQhAEKIqC7XYLWZbx\n+vr+9zNmbAyn0wmPxwOLxQK32w2bzQZVff8bH3P+e6PRCNFoFJVKBaqqolAoQNM0/TozNs7HTnuU\n68dO/LXY765h7Jb/7N4Mc3rO4XBAqVRCvV5HKpUCAPj9fqxWKwDAfD5HOBxGMBjEZrPB5XLB6XTC\nbreD1+v9n1s3jclkgvF4jPF4DJ/Ph06ng2g0yowNFg6HsVgsAACqquJ8PiMSiWC5XAJgzkZwOBz6\nic5ut+N6vUIURWb8DURR/PQ9IhQK6Z04m830x+qP/LOTdjweh6IoyGQyAN5mmNPXDIdDHI9HDAYD\n9Pt9CIKAZrOJdrsNTdPg8XiQSCQgCAIKhQIkScLtdkO1WoXVan28AP2WLMtotVrM2ECxWAzr9Rrp\ndFr/wsTlcuHl5YU5G6RYLKLRaCCXy+F6vaJWqyEQCDDjb/DMPSKbzUKWZUiSBKvVil6v96k1OHuc\niIjIJPhSmYiIyCRY2kRERCbB0iYiIjIJljYREZFJsLSJiIhMgqVNRERkEixtIiIik/gB+rhwGZD6\nREgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x103fbdeb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(tf[:10, :].toarray().T);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split the whole documents into training and test sets. The number of tokens in the training set is 480K. Sparsity of the term-frequency document matrix is 0.025%, which implies almost all components in the term-frequency matrix is zero. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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 = 480287\n",
      "Sparsity = 0.0253837\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 $d$ 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}\n",
    "where $\\theta_{d}$ is the topic distribution for document $d$ and $\\beta$ is the word distribution for the $K$ topics. We define a function that returns a tensor of the log-likelihood of documents given $\\theta_{d}$ and $\\beta$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "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": [
    "## LDA model\n",
    "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 $\\theta$ and $\\beta$, 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 $\\theta_{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 $\\theta_{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": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "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": "markdown",
   "metadata": {},
   "source": [
    "## Mini-batch\n",
    "To perform ADVI with stochastic variational inference for large datasets, whole training samples are splitted into mini-batches. PyMC3's ADVI function accepts a Python generator which send a list of mini-batches to the algorithm. Here is an example to make a generator.  \n",
    "\n",
    "TODO: replace the code using the new interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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": "markdown",
   "metadata": {},
   "source": [
    "The ADVI function replaces the values of Theano tensors with samples given by generators. We need to specify those tensors by a list. The order of the list should be the same with the mini-batches sent from the generator. Note that `doc_t` has been used in the model creation as the observation of the random variable named `doc`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The value of doc_t will be replaced with mini-batches\n",
    "minibatch_tensors = [doc_t]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To tell the algorithm that random variable `doc` is observed, we need to pass them as an `OrderedDict`. The key of `OrderedDict` is an observed random variable and the value is a scalar representing the scaling factor. Since the likelihood of the documents in mini-batches have been already scaled in the likelihood function, we set the scaling factor to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# observed_RVs = OrderedDict([(doc, n_samples_tr / minibatch_size)])\n",
    "observed_RVs = OrderedDict([(doc, 1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Encoder\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. The `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 automativally exponentiated (thus 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 have never see any significant improvement with this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class LDAEncoder:\n",
    "    \"\"\"Encode (term-frequency) document vectors to variational 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": 11,
   "metadata": {
    "collapsed": false
   },
   "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": 12,
   "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": 13,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0 [0%]: ELBO = -3722434.15\n",
      "Iteration 300 [10%]: Average ELBO = -3541472.9\n",
      "Iteration 600 [20%]: Average ELBO = -3448505.31\n",
      "Iteration 900 [30%]: Average ELBO = -3384165.84\n",
      "Iteration 1200 [40%]: Average ELBO = -3327991.5\n",
      "Iteration 1500 [50%]: Average ELBO = -3287262.03\n",
      "Iteration 1800 [60%]: Average ELBO = -3270108.58\n",
      "Iteration 2100 [70%]: Average ELBO = -3258014.51\n",
      "Iteration 2400 [80%]: Average ELBO = -3246227.58\n",
      "Iteration 2700 [90%]: Average ELBO = -3238394.33\n",
      "Finished [100%]: ELBO = -3242038.97\n",
      "CPU times: user 43.6 s, sys: 866 ms, total: 44.5 s\n",
      "Wall time: 39.7 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x115135d68>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAFVCAYAAACD77etAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8FPX9/197ZHNt7stAuG+IHEk4lEPwKrYiKqGgVWmh\n1aJ8PQAN/CyotDZpEbwAW0QKJmo4vKrVFhEkcggkXBKuQCAc4UhISLKbkGSz8/tjs5vZ3ZmdY2eP\nbN5PHz7YzHzmM5/5zMzn8573532oGIZhQBAEQRAE4QK1rxtAEARBEIT/QwIDQRAEQRCCkMBAEARB\nEIQgJDAQBEEQBCEICQwEQRAEQQhCAgNBEARBEIJo5RzU0NCAefPmoba2FjqdDjk5OUhMTMShQ4fw\n17/+FVqtFrfffjvmzJkDAFixYgV27NgBrVaLhQsXYvDgwaiursb8+fPR2NiIxMREZGdnIzg4GNu2\nbcOqVaug1WoxZcoUTJ06FQzD4NVXX8XJkyeh0+nw+uuvo0uXLjh//jwWLFgAtVqNPn364JVXXlG0\ncwiCIAiCsCBLw7Bx40akpqYiLy8PkyZNwpo1awAAr776KpYvX46PP/4YR44cwYkTJ3Ds2DEUFhZi\n06ZNWL58OZYsWQIAWLlyJSZNmoS8vDz0798f+fn5MJlMyMnJwbp165Cbm4sNGzagqqoKW7duRVNT\nE/Lz8zFv3jxkZ2cDALKzszF37lzk5eXBbDZj69atCnULQRAEQRBsZAkMM2bMwOzZswEA5eXliIyM\nhMFgQHNzM1JSUgAAY8aMwa5du1BUVITRo0cDAJKTk2E2m1FVVYUDBw5g7NixAIBx48Zhz549OHPm\nDLp16wa9Xo+goCBkZGRg3759KCoqspUdMmQIiouLAQDFxcXIyMiwq4MgCIIgCOURXJLYvHkz1q9f\nb7ctOzsbqampmDFjBkpKSrB27VoYjUbo9XpbmfDwcFy4cAEhISGIjo62224wGGA0GhEREWHbVldX\nZ7cNAMLCwji3azQatLS0gB2k0loHQRAEQRDKIygwZGZmIjMzk3Pf+vXrUVpaiqeeegpffPEFDAaD\nbZ/RaERUVBSCgoJgNBpt2w0GAyIjI22CQ2xsrE0g0Ov1nHXo9Xq7OsxmMzQaDdRqtV3ZyMhIl9fC\nMAxUKpXQJRMEQRAE4YAso8fVq1cjKSkJkydPRlhYGDQaDcLDw6HT6XDhwgWkpKRg586dmDNnDjQa\nDd544w3MnDkTly9fBsMwiI6ORlpaGgoKCvDggw+ioKAAGRkZ6NmzJ8rKylBbW4uQkBAUFhZi1qxZ\nAIDt27dj4sSJOHToEPr27QsAGDhwIPbv34/hw4ejoKAAo0aNctlulUqFigrSQoghISGC+koE1E/i\nob4SB/WTOKifxJOQECFcSASyBIYpU6YgKysLmzdvBsMwyMnJAWAxepw/fz7MZjNGjx6NwYMHAwDS\n09Mxbdo0MAyDxYsXAwBmz56NrKwsbNy4ETExMVi2bJnNi2LmzJlgGAaZmZlITEzEPffcg127dmH6\n9OkAYDN6zMrKwqJFi9Dc3IxevXph4sSJbncIQRAEQRDOqDpatkqSSMVB0rs4qJ/EQ30lDuoncVA/\niUcpDQMFbiIIgiAIQhASGAiCIAiCEIQEBoIgCIIgBCGBgSAIgiAIQUhgIAiCIAhCEBIYCIIgCIIQ\nhAQGgiAIgiAEIYGBIAiCIAhBSGAgCIIgCEIQEhgIgiAIghCEBAaCIAiCIAQhgYEgCIIgCEFIYCAI\ngiAIQhASGAiCIAiCEIQEBoIgCMIr1N80gWEYXzeDkAkJDARBEITHuXajAXPeKsC/vjnh66YQMiGB\ngSAIgvA4peU1AICdP1/2cUsIuZDAQBAEQRCEICQwEARBEB5HBZWvm0C4CQkMBEEQhMdRkbzQ7iGB\nwcOUXanDyfPVvm4GQRAEQbgFCQwe5rV1+/G3jw/6uhkEQRAE4RYkMBAEQRCEj2kP8SlIYCAIgiAI\nH7NkXSH+/vEBXzfDJVpfN4AgCIIgOjplV+t83QRBSMNAEARBEIQgJDAQBEEQHkdFfpXtHhIYCIIg\nCKIds/PIZaz9z3GPG06SwEAQBEF4HNIveI613xzHzp8v42ZTi0fPQwIDQRAEQSjMqQs3UGtsUrze\nkos3sP/ENc59ag8v+5DAQBAEQQQsDMPg4+9O4USZ64i7NYZGmBVS6V+vuYmcjw7gT2v2KlIfm+y8\nA3jvi6PcOz2sxiGBgSC8yKGSSuRtOdkugrQQRCBQerkWW4su4u+f8EfcvVRpxAsrduH9r44pcs7a\neotmwdDQrEh9ovHwsEICA0F4kXc+PYJtBy6hsuamr5tCEH7H9ZqbooTpGmMTzpbXiKrTZDILljl3\nuRYAsPfYVVF1CuGr7wHGwxIDCQwE4QNIw0B0NISW1w+VVOLF93bj0x2lgnW98O5OPLvsBzSLEAbE\nvGqeXvuXwnk3Ajh5elghgYEgCILwOT+fvQ4AKDhcLvoYU4sIgUFEPSoPzoSLPtiLM5f4tSGFJ67h\nOkvjWFXXKPtcJDAQRCDiR180hPJINZ4rOFyumDrcH7lUYRCcCK1vhBTtm6iiIgpJ1TAYGppFCSsA\ncKnCyGsbceGaAau+OIo/faCUcSQDU4sZe4qvoP6mSaE62yCBgSCIgMTQ0Iwdhy6JHtiV4uI1A37/\nt+3YduCi6GPWfXsC//x3saTzVNY0eMyozmxm0NisnE//og/24ZOtJba/rUsJl68bUddqIKiSZeIv\nLAyIkSmkCAymFjOefftHLP5gn4tzOpy1tfpaYxO++akMTa19W2O0CFGNrPgJLS1txxoamtFsEn8f\nGADbDlzC+18dw5y3CnD8XJXoY8VAAgNBEAHJP748ivX/PYnn3tmpaL2VNxpw2oWKee9xi6Ygb8sp\n/GlNmzr6WnU9GhpNqKoVZ9gnxEvv7cGzb//odj1cvPqv/Zi9bAfv/gOnKvDX3CLbRGdoaMZ/9pxD\nQ6O4r9r3vjgKs5nBy+/vxfPvyr8/rnrx6Nnr+Oe/i/HRd6c491++brQJk2q1s8DwXeEFvLP5iNO9\nsgpSV6rqAQCHT1fi3U+PuBRMrbX/89/F2PzDGfx373kAgJnjkJWf/2z7/ezbP+LFVbt563WEYSza\nHCtL8w+JPlYMJDAQBBGQnG21fG9oNOFmE/9EJsZwjs1L/9iDv+YWwWwWnvTLK41Y8/Ux3GwyYcE/\nf8IzbxZg/qrd2LDttKRzCtFsasGBUxWSrwWwGNm99N5uW38BwMXWSefUhRu262QYBuev1sHUYsaK\nz37G6Us1OHfFckzelpP4dEcpPt1xBgBww9CIHw+X8y7NHDpdaZtgbUVaZ1U+Wer4uSocOXPdbhu7\nbHVdI9799Aiqai32AMs3HMbeY1dx+Xq9U11nymvw8vt78c8vLVodLg3DJ1tLcOh0pVP0RMf2vb35\nCA6WVOK4izgP1jwa1jLVBotmQczSVW19mxZp18+Xbb8//N9Jp7IMw3jUlZMEBoIgAoajpdclDZif\nFZTiqTd+wOXrRsnnEmunYGYYNDTaTzo/HLok+Xyu2PTDGaz47Gd8vfuc9GO3n0ZlzU18vNX5Szzn\nowP4vsiytHL0bBVe/dd+rGYtnVi/zK9VNwCAzV14Wf4h/OvbEzhwsoK3n+pZ2oja+ibbefhYmn8I\nb206bLeNLbS99+VRHCypxMbt/MKYVXA8f8XiiVB0qqL1OvjP66hh4BMUXT0OjvKIGHsNx30Xrxnw\nwX+O2/7+4aDzM1R44hoOllTyN8RNSGAgCCIgOHXhBpZvPIy/f+wcoIdL9QvANsEWn5W+1it2VcFs\ndp4wxKyZMwxjW+sWwrrsce6KdJc869cv3/WcvHADR0uv482Nlsm68GSFbd/yDYft6rAKB5cqLQLY\nqi+O4vd/285Z79wVu2y/v9x51vZbymINe1KtrrV8tWtczP47j1i+0FVOSxD896PsSp3NFuBo6XX8\nVHyFs1yzyYzcLSdx8ZqBc781mBMA26TO91wCzvdDjCCcu4V7+UUpSGAgCCIgqLhh+cq1qtPZA64Y\nbUB1XaOd8Zkw4qY2481mpy9jMame39x4GH9ctkPQLoA9aToa27WYzTh/tc7ll2ybwMBfZvnGw5zb\nrev51vmXEbFMw4V9cCXxdbz3ZTGr3YxdW7hgGOBqVb2gwMZewlqaf8hmC7B842Hks5aTXlnbZvj4\nw8GL2H7gEt7efNjpEtQqFTaxNB81xibsKb6Cq9XOyyUAsLXwAj783wnb36cu3EATj/HjVzK0SnLR\neu1MBEEQXoQ9ZgsJDGYGmLdyF/ShQXjnubHi6hc5r91sasH5q/Zfna4mNStHW7UeC1f/hLf+bwxu\nNpkQHKRxEjZmsb7gG26aUHmjAfHRoQCAT3eU4r97z+P39w/A0N4J2Hf8KsYN6QS1WoWKGw2IjQy2\naT/45nox/gPWL3aZ8gJuGPiTNF2qNMLA+jpnCzanLtxAZc1NJESHtt1vFw3+avc5fPJ9CQZ1j3HY\nY9/wv+YecDq2tLzWadsFljbhWqvA6rj8ZG1TXb29hsBVGOqPWR4lgGVpSKvh/r7/vEA40JVSyBIY\nGhoaMG/ePNTW1kKn0yEnJweJiYnYs2cP3n77bQQFBSE2NhZ///vfERwcjBUrVmDHjh3QarVYuHAh\nBg8ejOrqasyfPx+NjY1ITExEdnY2goODsW3bNqxatQparRZTpkzB1KlTwTAMXn31VZw8eRI6nQ6v\nv/46unTpgvPnz2PBggVQq9Xo06cPXnnlFaX7hyAINzC1mPH2psMYM7gT7r8jQtKxZy/XorS8Fnel\np8g6N1tbYDYzuHzdiKraRkAF/HzmOqbd2du2/8gZi4pYiv2DO34OarUKpeW16BQfhhCd62G41tiE\nGmMTXnh3J4b3T8TsB1N5y54pr8VL/9iDN/9vDKLCdShszWp4/Fw1fjhYjtOXaqBSAT2SI/Hqv/Zj\nxIBE29c2n4ZByL1ya+EFW+Ahud4fP5e2GTM6VrHIIYGT4/JRrbEJxpvi7pv1/hafazNQLLl4w+mc\nFyuclxX+8mGhy7orbrT2AceTcanCiEsV0u1k2HjbPZgLWUsSGzduRGpqKvLy8jBp0iSsWbMGALBk\nyRKsWrUKubm56NatGzZt2oRjx46hsLAQmzZtwvLly7FkyRIAwMqVKzFp0iTk5eWhf//+yM/Ph8lk\nQk5ODtatW4fc3Fxs2LABVVVV2Lp1K5qampCfn4958+YhOzsbAJCdnY25c+ciLy8PZrMZW7duVahb\nCIJQgtLyWhSfq5YcYwAA/ry+EB99dwo3Wi3KyyuNdpb8rqirt/9itbrwLdtwCMvyD2HL/gsoZ1nP\nHzvnOpPhum9P4EeHCITuuEbW1TfjLx8W4s2Nh0UJKeWtE5g1rbHQuZes249aY5NNGGhhGJsr6Pr/\nnsTuo5Z1+H3Hr9m0L/U3TWgxm/HtT2V2dR0VsO/4eGsJqluDMikRaZCBZXLc9fNl/GfPOaf9jl/q\n6/97AkvWFaLWaNm+99g1W3vEkJ13AG9vPmL7+8cj4iNNcsGpYfAQs/62TbCM2HdGDLIEhhkzZmD2\n7NkAgPLyckRGRgIAcnNzERsbCwAwmUwIDg5GUVERRo8eDQBITk6G2WxGVVUVDhw4gLFjLaq/cePG\nYc+ePThz5gy6desGvV6PoKAgZGRkYN++fSgqKrKVHTJkCIqLLYNPcXExMjIy7OogCG/x373nsWVv\nmXDBDkLZlTp8vfuc7In0h4OXOF3TrNX9ac1e/Hk9/1ce2zq8qtZ+wuCybHeVlOhShQHvbD6CGkMj\nmk0tKDhcjn99e8KuDPsyL183unTd5KPkYg02crhYlly8Yfe3Y5wAoR6urmvE8+/utKnJHa9/y/4L\ntt9WV8XKmpt4ZnkBNv1wRmzznVAivkRjUwv+u/c8PvjPcc68Eo7LSxdbv9ytX+CmFjPmrdzldJxY\n/vXNCeFCfoKY7nb1zkhFUGDYvHkzJk2aZPf/0aNHoVKpMGPGDHz00Ue4++67AQDx8fEAgC1btmDf\nvn2YPHkyDAYDIiLaVJHh4eEwGAwwGo227eHh4airq7PbBgBhYWGc2zUaDVpaWuweTmsdnuC7wgt2\nKjOCAICN20/j3Y3KBkZpz7y2bj8+KyhFGSt5jpBtn9nMoPJGAxiGwYf/O4mlHCmINWIW/GEJJmTl\nr3lF9ufhGFld2TX849/FOHS6Ep//WIrtB7m/OK1fbtV1jXj5/b34a24RZzkhdrJ8661k59mvoTuF\n+ZU4L7eINC5okhHHgc2pizWoMfLbI4jlMxfr8mKvxZew3R8DCUEbhszMTGRmZnLuW79+PUpLS/HU\nU0/hu+++AwCsW7cOW7ZswQcffACdTge9Xg+jsW3txmAwIDIy0iY4xMbG2gQCvV4Pg6Ft7choNCIq\nKsqpDrPZDI1GAzXLfcZoNNo0Ha5ISJC2jgrAFtL0q2WTJR/rznl9TXtssy+Q009xseFIiAv3QGt8\nT7g+xNYnlQZ79bFjX735yQFsK7yAnGfG8JZZ8cVR9Osaw7ufC8cARv/61jnITWRUKOexCQkRti+3\ngsOXnfZZeSP/EL5aNhmNrWUvVhiRkBABjVYj2D4++K7t3c9+tisjVZuhchVoQGE2uqGhEMO6b/1f\nA2CNAhloyDJ6XL16NZKSkjB58mSEhYVBo7G8IO+99x6OHz+OdevWQafTAQDS0tLwxhtvYObMmbh8\n+TIYhkF0dDTS0tJQUFCABx98EAUFBcjIyEDPnj1RVlaG2tpahISEoLCwELNmzQIAbN++HRMnTsSh\nQ4fQt29fAMDAgQOxf/9+DB8+HAUFBRg1apRg2ysq5GshfHWsL0hIiGh3bfYVcvrpepURGldO2ApS\ndqUOoSFaJEZzT5BK1MNeM66pabD1yY0a+4HTsa+2FVpU4wtW7rQrwx5wT1+4gdMXbtjtByyq5yvX\n62FoaEavzq4/Fk6ed17qqKriNkI7eaaCN/6BY/t/OnwRr3/Ypln4cvsp/PtH+VbrFRV1ggLR4n/s\nkhycp/C49xJbXZMRBItoH8gSGKZMmYKsrCxs3rwZDMMgJycH169fx8qVK5GamopZs2ZBpVLhl7/8\nJaZPn4709HRMmzYNDMNg8eLFAIDZs2cjKysLGzduRExMDJYtW2bzopg5cyYYhkFmZiYSExNxzz33\nYNeuXZg+fToA2Iwes7KysGjRIjQ3N6NXr16YOHGiQt1CEB7G03loWby2bj8AYO2CO23bvtt/Ad2T\nI9AnJdqteipuNOB/+87bhd9l+7izFxMuVxpFDziujNbMDAO1SoWnl++AqTVRz5hbk0XW3AZ7HZ/N\nfBex+zc7fD3/rzUngJU1X7uvir7EYaHPxpOR/JTgDIf7IREYyBIY4uLibJ4RbI4ePcpZfs6cOZgz\nZ46oOsaPH4/x48fbbVOpVHjttdecynbv3h25ubkSWk4Q/oEvV2Hr6pvwyfeWZTbr5M8wjJ1//+cF\npeiaFIH0fgm2/Vy898VRp+iCdtpvlsTwZPZWPD91CAb3ihNsoyuzBbOZwaEzlTZhAQD2n7wmWKcj\nRayIhWL5xsGDoFBGHa7Y9MNpfPvTeeGCBOEDKNKjAEpY/RKEEw6P1aGSSsEvS8VOzTr3xQoDZuZs\nswv+02xqwVe7z9llzVu+gdu409HFDXDUMNjP/Cc4lgYcOXXhBvK/588HsPfYVaxgrekDkBih0X8h\nYYHwZ0hgIAgfwJYXmppb8M6nR7Dog3285YUwtZhbsxXaT5xcngBsz4WPOVL/sk0rrJ4A7EA39TdN\nONVqU6DVOKsC2IGALlyz1z7sOFSO7/ZfwLUbDTjK43mU89EBO08LR4rPSc/7QBCE+1BoaAFIv0B4\nArbmSqyb2NGz1xEeEoQeyc4Gfl/tOoevdp/DvcO7YPpdfWzbHf3vm01mOw0D27+/xWyGRq22a8/y\nDYfw7vPj7OpYmn8QZVfqsGhGBjQc4WpPXbyBxWv3oUdypFPQmIZGEz75vsS2JCKHn4q9Z8BHEEQb\npGEgCD+j/mazLfsgm+UbDvMGYSltnZgdDfnYk/+hkko89cYP+K6wrQx7+eDLnWdx4ZoBLSwVg6N7\nImDxlgCAq9X10HIYG+S1ZsxTMsIcQRC+hzQMBOED2F/5jqsGOR8dxMUKA5bMGoGUBL1b52FrGN75\n1BL+9j972gz32BqGr3eX4evdZQgNbosj0GQy2+KQOFJV24jzPKl8CYIIPEjDIAStSRAegLH77Rjq\n1jIJX5UQ/MXxO7+0vBZ19U2Cyx1caX4dY+GzNRJsHF0MCYIIbEhgIAhfwFIrsDUMZSwXxWae7HQ1\nhka0mM34z55zqGzNFcCWGGoMjfjLh4V4+f29nDkU2AiFbiYIgrBCAoMAXKlKCcJd7DQMLInBGhwJ\nAEwm7mdv+8FL2HP0Kj7dUYql+c65F7YfvATAkspXKCUul40CQRAEFyQwEISbOKZSFgVLFuBTApy+\nZHFdvNlkshMqGMZicAgAFTdugmEY3GBFRvz3rnO23//6xnXkQaHUxQRBEFZIYBCA4jYRXLAn8Ofe\n2emipIXv9l/ALlZWwg+3nERV7U0s33gI+TwuhiE6La7daMDTywuQx4qXwICxM1z8bv8FW4pfR9jx\nEwiCINyBvCQIQgZSBMkPvj6GXUev2G07fbEGa785jmMuJnSVCihpDZC0/cAl23bHnFW7HeomiEBk\n7OBk/HjEORU44T1Iw0AQMpBi2+IoLFipMbheylBBxRmpkb2NjBaJjoI+NMjXTejwkIaBIGSw/4T0\nZEeOXKp0nQZ4+8FLaORIs9zCSrqkUaucfSoJIhCh59znkIbBz2lqbsHR0uuC7nGEd1n972MePweX\nsADYx0XQatROCZ4IIhDpKM/5rT3joAvyz6mZNAwCKGX06Jg+WCy5W05i189X8Pgv+mHCsM7KNIYQ\nZO1/jqN7cgTuTEsRVb6h0YTQYPvXqar2pqx7LoWbTS2oa5DhpUEQ7YhZvxqAKxICmbVnXvj1ENTV\nN6HoVAV+OHDJr6Kp+qcYQ9gobnV7K7tCcfm9BcMw2PnzZVtOBDH8e9dZ229Tixm19U2Yv2o35q3c\n5Ykm2lFV2yhciCDaMb06R/m6CV4hMToUABARpsP4oZ3xu18O8HGL7CGBQRD/WAog907vwe7qopMV\ndi6UfNxgGTD+NbcIz4twtSQIf0Kn9f108OsJvTm3a9UqRIXrPHZefzGoDHK4B91uiRA8ZtSgJE81\nxwnfPyEdBLnzvadV2oRrVn7+M4pOVgiWY9uYnGOFdyYIV3hzsBeCyyPH20wc2ZVzu0ajxh1DO+PB\nMT283CJpOE74kpEx3D88rqd755QACQwC+PodEvN1SyiMQ5dfrqrH2cu1LgfU05dq8Of1hbjoR+uN\nhP8jND9YVdTewFWiskUzMrzWDi40ahWCtGo8IFFgePKBgW6dt3/XaEnlX//9SLfOx0VwkEa4kJcg\ngaGdQIoG92hoNOGzgjOornNe7z9/tQ6zcrbh8OlKAM4xFj4vKMWf1xfi93/bjkIed8rqukacvVyL\ndf89oXzjfcBD43ri1d8Nd6uOHsmRCrVGmDuGdvLauZRESIPozVw2rr5NeiRHIiUhXHbdq18cL/tY\nAAjRyZs0Rw28RVQ5vg+zBKkCm8PtXPnCOIwZnCytDscqBcZ+b3qPkMDgLRT0tuiIMAyDTdtP49g5\nebkPvvmpDF/vLsM/vzzqtO9/+y6AAfD25iOcAgWbVV8cxftf8btUBpJgJ3mwdOA+lno5LjLY3ea4\nJCy4fTp8CT0v/vW6y3+43X0vdAp/ZWc9Ogx9UpQ3pGRP3j07RSI0WItgrfi2c3WToMDgxTGHBAYB\n/OV9ZRigtr4Js/62Hd/uLRM+IMC4Wt2Ab/eexxv5h2Qdb00Qdc2aDpoF+4VbtGav4CC9p5g/FLM6\ngCQGrca9a2F3xdjBHtYAiGxqz078Wo9nMwcLHj+sT7zYFolC6HnxlsDw8uPpHq3f8Su4bxfxqv7w\nEOWFQY1GjWh9mxCrVD+zb6erJR6JtSpUj/uQwOAl5KoW2SpL69f1pu1nFGlTe8JVmubdRy/jmTcL\nUGPg1g4wDGMzXDSbGbSYzagxtnk1sF/y+kYTytwwWgwYI1WGgUbj7vDQ1hee7hYlBLV+Iiax+2/v\n7rTtltgw2efka3ab+l/euLHsmdGSyrPdFoO0ajzxi35OZdzqYjeOlWq3IAbH5igpl2X0TwRgCbqn\nBHK7btLt3RUXNUhgEMJfVAwdHYf70NjcgpbWLExrvj6OhkYTCjm8GQwNzfi5tArGmyYAllTSyzcc\nxgvv7sQf/r4d+09cc5rkX88tkt3MU63JovwRjVqFWb8S79etpLYkvV+iy/2De8UBALom6Z32JcUI\nL414yy1OyyFExUeFyK6PT8C0bpc7/GjUbfVKXQ7KHN8L4zmCxPlIXvCp1k7qqVUqFYJbPSVsAoNC\nzb8rXVwQOSsPjesJnUzbDz5IYGhHdJTQqFw4Dpyzl+3AH/7+g51XgvXlLj5bhfzvS5C35SSefftH\nvLXpsK1Mi5nB8bJq2+/3vjjaYXo1MSYUo291zwBLCuzBNjRYi07x/EZzj97dBy8/no4Xpg5x2tez\nk/Bas9jrclf1rOZ6WCQ+QGzhRqVSYdXccbxVym0vu+/T+roW1hwRG49h6vheEtqjnP3DMw/dKnny\ndK7UcYNyX4aDesYCkN7vfAjbMHhvBGuflkJeRCkrZf8yXgocFq/dZ/utAvDt3jKXSzYNjSanbXw5\nGwINKQOLEo8r+2wqFdA7JRrlPAm3tBq1TS1+d0YKthZetO3jnKQd0KhV6N81GifOy9fwiOkerj6U\nKsi/+rvhmL9qd+uxQIjOeRh2dxJgHy+1Kh2fkZ5jPW400dWhsZHBdtFLHcum90tAz06R+L7oIuSi\ngsquX5Qcn0cOSEKnuHB0luVV4qJn/GAOIQ0D0S4Q4x2Su+WULPuOfcfdzzzZnnhzjvj17ZUvOH/9\nikZl//PpTGftASeOt9qHKqAhrUslVrgmX8kTMtvin+dYa518z31MhHteJ4t/6xxXYfaDqeiSqMdQ\nHsNOb2nQv4MAAAAgAElEQVQ4X5w+DL8Y0cWj5/DUR7lKZRHWuiZFQKNWZnq1PgKy2qywkEECgwCk\nGfAsn2wtwf/2nRcsx46e+OPhck82yecoEVOAa93fOt5E6YUnG+taPdeavWhY745KpUJ4aBCeeSiV\nsyj7i9jxlRMz8Ipf5+Z/ocVMiEqvp6t5zqluVas4jj+R4To8fm9f/L/HXHs1CDWz+y3O3iLD+yfi\ntZkjnJKoeZuk2DBMu7OPyzJKT/iungpX3H97Nzz9YNsz7Qk5xKrl5rd38cBJeSCBgfAp3xVewIZt\np12WuXajAeu+bQuI9J+fAtut9NaeccKFBOAcXHgt8u0NDdP6JsjOjNqL5bZoNwi3nluMsOKIRsya\nhEq+Gj8xJhQLfpPG2T9Oyg7x3cobspddB1+Tx7YG+3n4Doc6GAYT0lIQx2NoOWJAIv7fY+ke0QbE\nejiWhhRcXd3oW8UFa2IjN77NXWkpNq8IAJw3VMqd4HoezDwahp6dIoUjcCr8GJDAQPglJ89XY+Hq\nn3DkTCWqa2/a7esoNgfuIGXunDPlVtvvmIhgzHn4VttXppR6grRqvPwE9wAmVI2rAVuMwODOuPjS\nI8PQt0u0qDo4bRh4OonLBVPssb06RWHNSxOc4ldw9dKMiW0ukH+cnIreKVEe+eqcMbE/7hvFnevB\no0hcB5r1K+Fw0I6Hyw785VARZ6vcvBfWd8OxzVPH90KP5EivrtiRwOAl3F3a6GgrI3/7+CCuVtXj\nrU1HYGhotttXw8oMSXAj9mt7xIBExQYcV8+4tT2izuVQj1qMwKASJ9xY29g10dl9UwxKLEk4GoPy\nFeK6bq4+lqO1kUNkuA5Tx3Nnk5SKlPHMXccUrpDWjhqYub8eKqFGF3hg9m6zYRAwePECJDAQXudg\nSYWTEGBl18+XkbflpN22lZ87h3MOZLxpN3N7arLdGOfx0OM8Y5ur04pakpCohI9mGQ3ahBlRXhIS\nTsJfiwIlXEfiZAs2gRJLzA4J1/TSo2ku9yfFhiGFR4AUzuPgeazvpD9EkSW3SgHI6FE6/9lzDlV1\njXj8XudocYUnrmHVF0cxoFsMXnxkmG17/U0T3tx4CGfKa73YUs+R2jMWR0vl5b3o3dn9pE1cvvS8\nH7NeGIiE5nxX7ssaESGqVSoIju6/uq0bis+6uifCX/TcSxKCzeM/I/9N4dzMFuhWPD8OphYGpy5y\nuJL6fm5xiaTmcRSWMnnqQ4PcdrnlReXyT0UQ8pKgJYmApONIHp/uKMX2A5c4952+VAMAOFNeY7c9\nd8vJgBEWAKC3iGBDfETpg/HMQ7cKF3SBk7EcD2JV+e5jXZLgUzGwf9q/K2Ld04Qu4+70FLffQiX6\nys7okafVYk6jC9IgjCfPgrflhRXPu+F+6wUcBT27P118FQr1o+N+d+N0cAn6Zsa1l4Q3IYFBkI4z\n0bP5309luFpVL/k4vqUGK82tOSF0Wo1dEJ/zV+Xnb/BLeN7tob3FJS9K75cg+9QTR3RFZJhOdKOU\nG4hcDLytp9AFcQ85/JlCxAVusmgYxBSSvdtShvM47w3kXHMbd5s83hQ7+AQXTyH1+riW2vxhAmYz\npFccZt3vbLBpbTrve0BulYGHEksb3nq+z12pxYpNh7Bw9U8AgGaTmTdCHwCLJ0NrWugdh7g1C4BF\nq1DCyrXwpzV7bb9rjf5ryDisTzweudu1X7gjfLfKGhZYifXI2Q9yxzQAxA+GQs+lO83kqrszX3ho\nVmGnw0Q0QtX6n+syPNslXCP3F6QF3mtTELHDiF2kRz9cn7Dm35AXDdH9OVKlUilir+NScyGR56YO\n4UxkNqY17HlvnnTc3ry7ZMMgQEfUL1gTNVl597MjOFpahVd+Oxzdbomw23f5uhFvbToCAMh+cpSd\n+viHg5eQFBuGpZ8cxILfpCHnowO2fY6aCMdz+hP/N8WS9viTrSXiDxIIsqJSwe2Hq19X/uyKUr46\n7ULkutGeIBc5CKwCEp8gow/l0ohIQMzSisAaMNfxTrYV7prsO5wnMpz7uvmuRfEskj5iUI9YDOkd\nj75do/H8OztdlpW/aON9PNGqGff1wy9GdvWKQCoEaRh8BMMwKK802jIu+hPsh766rtFmvHeBlejJ\nyg2Wi+OyDYfsBtgP/3cSSz85CAB4I/+gZxrrp/CpD119+S+ZNUJSWuKQIGUy0SmhmlWpgPnThzls\nZYd6dH28uyptcX4UkC4RSZAX5Ahbd6WLD5B1d0YKRg5MctrOuUzhB/PpHycP4t+psqSB5l46E8bd\n67PYyIpb6xrF0eeSzuVmWzVqNaewoBIQwj0BCQxewvGdPnzmOv60Zi8+2nJKXgUeovhsFdZ+c9z2\n95sb2zI9MmBQWl6Lj747BbM1/BhrtKqsucmby8HU0hF1Nc64Wo9PSdBLyhGgC9LgngyemPsc5+nC\n5zrm4hxi1dnPZQ5Bj2R+7w5XtTi5B4qYpJ3qF7NsYaemt9vRuk1aHY7Hs1XcQn367vNjsXT27Qji\nS/TEgZRlLH9Ynx8xwL2J1qNI0Kr9YdJAhAZz3yc/6GY7kuOclzSUhAQGATzlVmldy9919IrLct5+\nIJdtOGSXKe5ihb1W4S8fFuL7oos4Unrduw3zEeOGyMvr0Lszz3qjBJ9/MUyd0Au/v3+A83k4yzoH\n3RFlLKgAOgnaELmvnNzLUDn9kFapHEPI8JAg3vDOgLRxR0q46vaI0q6sQKuGQcL5g3hyqvhbP/9i\nhGcjcZLA4C14BoD2FOfhzKU2t8em1vDM7aj5shjUI1bWcclx4Xjr2TFO21UuZyfpaDVq3J6aLKps\nMIeHAsP4R0AYudw3sqtNqBOapNUqGenqHYrzTRyA/bvssikyu9td48xAwn8MOYXbEa1gJM7f3tff\nuQXkJdEBEHmT/UmgKL/e5inBMEDZlTo0m/zPBkNJZFtSq+ByfdbjL7mb/uCsg7yPhD6fkNaZcxDl\nhrUkYR8MwWEvqykOfwfrnDUl1qpEezCIKcOpNfDNJPnIXdI8hLyC+24SXjkNANyVnqJALRbkajyV\ngrwkfI7IYcbDY0VDown7T1xzWcZmtwDg693ncMmFq2WgwJ67NGoVWszi7pf1q71roh7nWcai3hr0\npZzFzq5CQQHVG8Kuyk4IECgruN91AcH1YbYNgyduM2+dQhfGvVls9MMQnvV7X+Ju94qKNi4CMffZ\nlfeQIm1wsU8fEoTGJuWS9cm6koaGBjz99NN47LHHMHPmTFy7Zj/R/OMf/8DcuXNtf69YsQJTp07F\nI488giNHLC541dXVmDVrFh577DHMnTsXjY2WdfNt27YhMzMT06dPx6ZNmwBYvvJeeeUVTJ8+HU88\n8QQuXLgAADh//jweffRRPPbYY3jttdfkXIogSsXWd1SFWgc60dV7ePD96LtTdimkuShlRWLsCMIC\nIEOF7QJ26mdPig28bZY74QQSIgMfScEqaDAc27jLyz0P3x55z6hQjgV/gdtGxL067ZbhWrvvj5MH\n4ZmH+GObiMGPFMIAgBd+PURWum8+ZAkMGzduRGpqKvLy8jBp0iS8//77tn07duzAjh07bC/MsWPH\nUFhYiE2bNmH58uVYsmQJAGDlypWYNGkS8vLy0L9/f+Tn58NkMiEnJwfr1q1Dbm4uNmzYgKqqKmzd\nuhVNTU3Iz8/HvHnzkJ2dDQDIzs7G3LlzkZeXB7PZjK1bt7rbH15D7APvrfWpsivtK9LigG4x3jmR\n3BUJvvvmpfsp7fmSHxrXFdYYA6ESUwc7tkapd4A/DoWLEwhI9FoNS8ch3ou0XfDwOEt48UHd5dnx\neBb3epjrlo8YkIT0folu1esb+PuiU3y4qHTfYpElMMyYMQOzZ88GAJSXlyMqymIRXlZWhk2bNuHZ\nZ5+1lS0qKsLo0Rbf8uTkZJjNZlRVVeHAgQMYO3YsAGDcuHHYs2cPzpw5g27dukGv1yMoKAgZGRnY\nt28fioqKbGWHDBmC4uJiAEBxcTEyMjLs6lAaf5MYPUVH0RhIxWNpyf1kRnG8PiWf9z4pUfjdff3x\n2u+GC5RUpjOElhQ8sRz02swRtuYzoiUGpdvhmYfp/tu7Y03WBMRG2ntzTLmjJ35zT1+PnFMsbntJ\niJaoxbWjT2sUxmABbyC+0Oju4E2jR0HRf/PmzVi/fr3dtuzsbKSmpmLGjBkoKSnB2rVrUV9fjyVL\nlmDp0qUoKWmLiGcwGBAT0/Y1GB4eDoPBAKPRiIiICNu2uro6u20AEBYWxrldo9GgpaXFbrnAWoe/\n4k/Gix0drUYlOi4EexJQJLy37V/vSwx8Z3R1We5Y26tUKoyVYaTlhp2pYAElX8PbU2/BLbFhbYGb\n7OQF5ZckZGu7ZJ4O4Pag+dVt3QEA3+4tk12vpOefywDUbYGBu17XB3Bsbq0k6zdpMJnM0PJ40bz+\nh5FoNpnRNSkCM3O2SWytpCZ5FEGBITMzE5mZmZz71q9fj7Nnz+LJJ59EVlYWrl+/jueffx61tbWo\nqKjA+++/j4iICBiNbV+vBoMBkZGRNsEhNjbWJhDo9XoYDG0GYkajEVFRUdDr9XZ1mM1maDQaqFlh\niI1GIyIjhdMCJyRECJZhozO0xSSQeiyb+Hg9wkKCbH+HtVrQMwL1Wq8xOESLyMhQRdpy+FQFdEEa\nDJDpMqgEvbtE4/QF+elmdRzW6lYG947HkdOVvPsjw3V2sSZcERER0tbXKogetOPj9NCH6ewC82iD\nNAhtve9cCRil3lO+8qGhOsTFOQdpik+IcPoCiooKRTyrrFqtkvVsRUWFOh3HVQ/XNrXKfntoaJDd\n/vDwYPTvFoMTZdVOx8bF6ZEQY3kvdAJLHwkJEdC2GqAFs8rGx+sRwePREqRrK6fVqu3aufB3IwEA\nIa3vtYY1WahU/PcnPj5C0BAuNjbcdrxKZRFGwsJ0nHVGX29LEse1P4zj2sTcKyH04W3ugkL3OiEh\nAnMfTcPyjy3h4SMjQ0Q/H5ERzmVdGR87nhcAgoLsn434uAiEBLfdN77rDwu19Lmax0oyPkGPEJ3w\nkpuY/pV6D6Kjw5CQEIH6m21h9iMi2rRB7swRfMjykli9ejWSkpIwefJkhIaGQqPR4O6778bdd98N\nANi3bx82bNiAP/zhDyguLsYbb7yBmTNn4vLly2AYBtHR0UhLS0NBQQEefPBBFBQUICMjAz179kRZ\nWRlqa2sREhKCwsJCzJo1CwCwfft2TJw4EYcOHULfvhZ12MCBA7F//34MHz4cBQUFGDVqlGDbKyqk\naSHYSZGkHsumstJgt5bb0NBaL+O6XnNr6OibN5tRW9sgui1mhsFbmw6jpYXBvGlD7R74P/1zNwBg\n7YI7JV+HUqR2j3FLYGhyYfk7pFccgrVqXq8PsZ4OAFBbe9PW11IMYCuvG9BgDILJ1NZOU3OL7b5z\nNUHq88VXvqGhCVVVzktM1yvrnCIL1tQ0oN7QJoiazYys57ympt7pOMe/ExIiOOs2O7wDDY55RoyN\n+L+Hb8UzbxY4HXv9ugEwWfKQNDW6zkdyvbIOplY34KamtrLXrxtw0xjEeYzd/TOZ7dpp/d3Yet6W\nljYXYwb896eyso73S9RKVZURwa2vrFajRrPJDGN9E2edNTWuxwUtx1wndK/EYDC2Cd1cxzv2VSor\n9wn7vRJqR22dc1mzi3eR6x41N9s/G9XVRtxstDxnppYWl+9SRUWdnYcYm8pKg+AyhFjE3oNnHkrF\n7qNXEBumRUVFHW6ynuW6upuc9SklPMgSGKZMmYKsrCxs3rwZDMPYjBC5GDRoENLT0zFt2jQwDIPF\nixcDAGbPno2srCxs3LgRMTExWLZsGbRaLRYuXIiZM2eCYRhkZmYiMTER99xzD3bt2oXp06cDgO18\nWVlZWLRoEZqbm9GrVy9MnDhRzuX4FLFW+K5UeKYWZ1VY8dkqWw6IolMVKLtSBzPDYOr4XrYyZVfq\nkBwXhve/Oiaj5e7hJTdqtxGSERY8MRw5H+532i43c6ISMIxUbasK44d2wg+Hyt04q5JX5NzpUg0n\nuVGJTg9tZdTAJNt7xF+rBSWXHXnDWMvgzrQUbNh22s1anLFmVuzLk0XRk0jtE0c7DKnjh49tmO1I\n75fIa5xptaW4Y6hn4jXIegvj4uKwZs0a3v0jRozAiBEjbH/PmTMHc+bMEVXH+PHjMX78eLttKpWK\n022ye/fuyM3Nldh6aXjO9ECZR63sSh1eW7cf0+7sbRcWlO17u+/4VRSdrABgnyXyjfyDuCU2DGdY\n7pKBgEqlnEBhnx9AZSfgaTUqhDokTXo2czBKy2t4kyl5zXaB8zQizi171pP/prhzr6RFP5Rev06r\nQUb/RBS6ilHC4Sbh6lTeXnv2VByAob3j8WzmYI8LDFzvjKNtzZBecTh8hj9c/SN390FCdCi+3HmW\n83ghuiTpUSMgOPoDyXHhWPH8ON7cF+5CkR69BO847KZEcuCURRBw9QVRVdumptp55LLtt/GmySPC\nQr8u0Xj5iXSXZTwZulbJmjWsxEhitEFDe8fj4XG9BMsBwEuPOGZ39CEK3o/kuDCEu5F9UorMIjXn\nQu/OFjunrklsFa172iBbTBWx7fAXFxk3UalUGNo73s42y0pXnkRnnuD1P4zEMw/f6rJMeEgQJo/p\nYftbpRIZcbO11JOTBuGxe/uitwLC0bOZg/H0g+7Fe3BFWIjWY+MrRXoUwkPuDVJvJ99kpWV9PfxU\nfAU3DE2YONI+AcnZy971Hnnp0WHCLm4+HDOlnHp4f/6Me3IeDfZ18xlSKYHkmhV8zv/y+5FuVecp\nrZ4KKky7qw8GdI/F0N7xtq9NBSoGICGXhFw80DFdE/VoaHJt+yGVBb9JQ0qCczpmT5Ec53wuIU8o\ntveHmGdVHxqEO9NSUHalDqcv1shqp5WhvePdOt4RbwqgJDD4GKFnVehRYCfDWd1qi3BXemefxo+Q\nmmrYE+cfMzgZ+45zq5Gl9A1bnTt5dA98odAk4/H744EMf+JPrfLiuaQUtvjJD+9vv/7LVYd1whGj\nKZF8qSIOUPH+ocD5WbwiGCNDOn27RAsXkoCcZ2nlC+NgdpHmRu7zOe3O3kiMCcWnO0rlVdDOoSUJ\nAZQb2B1CQ0t8YHf9zJ0GW6txrsjbGgU5uG30KLA/tUcc/jl/vJtnseeBMT2QynJFDaiEgD68GKcz\ne2pJQkQZjVqF3ilRWDJrJKbf1Qf9WRFFpWlNvPjVJ+Bx4QqLcOeDmCAePmWQVsOZKMx2fgk5SNiE\nhQTZ4lB0REjD4Cdcvm7E5ev1SOubYLedzwOwqvYm8racQmJMqNO+nI8OYKKH86K7jQcHDGvVnjD2\ncjUIScLDKobAkmWUuRox9fxz/vhWo1mVzQtAuF7Lv2Ldbrla8db/jcHz7+60/c1VE9+yZN+u0ZiQ\n1hkjHDQnS2aNgE4hlz9/ZNVLd6KxXlw8FUfsUkm4jFzGvTkxJhTXqhvsbJx8hZL5boQggcFLCN3S\nl9/fCwB49/mxCGcZEfH5/+ZvO41DLoIT/Xffeclt9CZur7u5GvxlVC02c5+7k5e/ayXaQ0DS5Lgw\nXG4NVqREf7pvV+Js9Ci1lshwHW6JDcOVqnrhwg6oVSo8fm8/p+0pCd4zPPQFXZIiUFEh71jLPZf/\n8Pzl9yPR2NwCDVcEtgCmY12tDJSyBTslMkiR1VCn6OQ1ZP1jN2rYgaNutAVocenm5UF6d46yy7oo\nF3dT5rrlisexjctwSuhY9wz7PDs1c10ju8+sRmkJ0SG85X2FUN/cnsqTfU/mRbgrvNr6VWymep6H\nl314fFQIZ5n2xrvPj8WK58e5VYcnhGx377lWo7b7sHOXKD13pFExeNPokQQGL/Hupz/bfjMMg3/v\nOsdZTtP6hbPy86OouHHTbp8/GNpoNSo8MbG/pGO4LKbjI90bEH9zT1/07RKNR+/u47RPzgs0ZnAy\nnnxAOKsbe/CS85L7Ng5DGwsfS8fi32aIFpR8iT8JM1xwyQvuTHJjh3a2s0toz26Y4SFBdjFJnnko\nFT07RWJIL2FPgajWbKfWrKdK4k+avhXPj8PS2bf7uhmiIIHBB7hy99lx6BL+u9e3ywljBifz7lOr\nVeiSqMc9GV1E1/fMQ659pIXow+H7nBQThgW/SUNyvPOEJzQYcH3hqVUqjBooLW/8vGlDJZW3w8e6\n/9BgLbrf4r6mSBH8aPDmQnByaTNicFnM3WiVgZDALr1fIv70RIYoW6A/PZGBGRP7eSS9tkLJKhUh\nLETrluGqN+0oyIbBz/AHLYKrx89qU/HI3X3wixFdMH/VbsH6PPGF4A5SckI4YhU2kmJCkRwXjgtV\nDQJHOFYg+9QST8PhVunvM7MVP5sYhR4XsSsSy58ZbRf3X6g+/g3y6Zqox/23d1euQgVZPme03d9x\nUSG4Y2hnj5xLrfae66+n0WrUWDJzhFfGWRIYBHBncuHiYEkFdh/ldpH0F1wZ9rG1IyFKeQwojCcH\nAonL1bww8HA7ZdYdCF+xvkKo74J1rl39bPUo1B4u0volIKM/dx4CXxOtDxYu5CZvPTsGhvpmaNRq\nW9IoMdkm/Z0UL0XWbP891Y4ou1JnZ8tgpba+iaO0bxgxINHlZPMIy2ZAyGNg8pge+Ln0unKuiBx4\n/SNBogsdz+G+wdXJ3WxYIAsawksSXmmGW8REBKO6rhFhiiTxUpbsJ0ehsZk/+6ySRIbpENma7nvy\n2B5oaDLhgdE9BI4irJANg0KcKKvGrp8vuyyzp5hbs/D8Ozs5t3uLhY+l2X7fN7KbS0GgC0uSFRpI\nJ4/pgT89kWEXhlUOUuciOap3qWua7k6Q/jbBWgdRf7DO97OuEaTtmWAln3LjmfeE/PHSI8MwcWRX\nj2UxdIek2DCH3B7eITJMhycnDXIdb6MdCIPexP/EzXbK3z85CADI6JeIixUG9OocZZcxEgC27L/g\ni6a55MkHBqJPSlsoVymZHn0RIU4UMpolfgJ3N36Ee4eLPo3EezNxRFeYWswYL3PN2D03V/uDne6F\nnz5mVjxtG6JE7UmxYfj1hN4K1BT4zH4wFe99cdTXzfBLSMPgJrXGJlxjxUdY+fnPeD23CCfKqvHq\nv/b5sGXA1AnCWRO5PAPULoYo9sTgD+M41zwv1K5gN9YsRRrEi8D739Gu+iVYp8GUO3ohzg80DP6q\nY+BtlWLPBOEPDO+faEsP3W4Mhb0EaRgEcDUIHDlTibc2HbHbdvSsJWf6pUojrlZLtKBXmBC5YWFd\nBlFsBy+QQBPT+sSjxczI0vi0CUxybRjaQf+1A6aO74V9J64pYyinkNLI0/KCvwsko2+9xedjHuFZ\nSGBwg6KT/HFJP/rulBdbwo1cY8N7h3fB90UXOfexNQz+EEddDmq1ClMn9MbxsmpcuGaQdKx1wufL\n8cFxACdyB//Rqbdgl4CXDV/dnlxC8uRk5ihkqVQq3DeqG+4b1c2tep/LHIyyq3U2a3nZcAiRSnZ1\nZLgON5sa/NYrycqsXwkHPvMkf599m2IeD4/d0w/vf33MZUyajggtSQjgKkwtX54HfyFJZPIcNiqV\nCgnRzgmt2PutaNRq/HnWCFlt8ySe/Iof0N2SvdAxRTIvDo+Iv5p9+BNCb5VSrs5DescrYiE/uFcc\nAGDCsBS36rGGvB4+MMlu+/NTh2DckE6YONLPE8r5mPioUOhDlQnXfFvqLfgga4KdkTdBGga3MPu5\njtDVxA+4lxbXSme5CW4UmjiFciYo3YRRA5PQJVGP5LhWYcy/H4F2AVfG1fZEao84LJ8zGlHhOnzz\nUxkAec/W/bd1w22DkjCgdyIqKtpS1N8SG4bf3ictHDvhPn5r1O1DSMMghIsJQWxCKSXQyMiiJ+TO\neHeG8xeR114RVr8u+E0afzk/Q6VSISVB75MsdZ3iwwNOPnlgdHc8O2WwpGN8PZCn903AKActQLQ+\nWIFMpirER7Vv4YkIbDq0wHD2cq1kwzdTixkNjZbwrtdr5eVil8P7L02Q7CPPJTD4Y8rbvl2icd8o\n5dStvp5QPMGLjwzDy4+n+7oZivPg2J7OnhkCUpHS0Vel8szDt+LJBwa5LhSAzyBBdGiB4c/rC5H/\nfYmdW6Qj7KHp4+9O4cmlP+CZNwtw+HSl5xvowOLfDpdUXq22H7e6JOqxRMjmwFvjnBfCNwvhatrp\nFB+Op4QmBSUb44I3nr4dA7rFuExe1O0W7we+8RT+GrqYIDo6HVpgsNIsMizpVpbnwN7jVz3VHCd+\ndZvFGlyqQY9KpbJTnVs/zHTadnbbvfBBGR1h7573pyfSMdJB7cyJGzkbxBpnitGYLJ6RYbP2d2Wo\n2x5I65uAd54by7vfZxokCZoN0i8QgUg7mzn8Bzk2BXIZ1idB1nFqlUpyOzvSQGe91uS4METJzfQm\ncW52d67jO1ylUkHNutftXSOulLW7lft87GHgjyGZCUIqJDAI8MHXxzi37/rZ9xknbZb6PKjVsJtE\nnGjfH6IAuL82pU6WSTHS3U9Fw9PHcrUAAXDLfMJUb4dFdngGB3SL8e75CcIDkMAA14PwqYs1XmsH\nH3yTi5BqVq1S4d7hXaSdrJ19mvbtEoW0vgl45qFU1lb3r8FzsRzaV/9y0Sk+HAAQrZeplWnlxUeG\nKdEcjxOIRrQEIQcSGHzAr27rhszxwnkehGCPY1xfMCqVCun95C1ntBc0ajXmPHwr0vu1Gcr51fjO\n1xYJqgK/uh4AocEafJA1AQ+N7elWPYH81f2wm31DEP4ICQw+YModvfBLKWFt+UL9sn4PH8BtWc45\n17jMFeEjFNS1+9n8aoc/t00KKpXKa8sj/iYwiWHsELJZIAIPEhjgnUG8Z6dISeXnTx8qolRby7V8\ngYTa42jrLkpcsoe7Tcpk2x7uoD/YVni6n/zhGonAYUJaZ0y6vbuvmyEJEhjgnYHgwTHux6x3xC7V\ntKjRUvhK/VG+kGogGKRAyGtPIad/JT+fHk+b6OH6CaID8Pi9/fDQuPa1dEW5JNoDPJMM39zz8Lie\ntjBr4ukAACAASURBVGBUruYnf/LXl+3WyGLhY2nYU3wVA7vHKtAi/8Ffk5yxIy76k5yZHBeGy9fr\nFa/Xn66RIHwBCQwuqK4TH/r5yUkDsforbhdMuTyXORh7j19Fj2Tu5Qw+6+37WWouf9QYcDEhLQVN\nJjM+KyiVXUeflGj0SYlWpD3+1G1ykpx5o/1sl11PijRiryWiVeiM1gd7RGAgiI6O/+pufcy+41cx\nb+Uu0eUFU/LKaMOQ3vF4ctIg/iRSXprVvHGaIK3aTtAh2hCrYPC2kDNygIhImF7k4XE9ce/wLvj9\n/QN93RSCCEhIw8CBmWHwjy+LfdoGMXOE3AnCVYwBb/mcey7OgX8jp38ZDolBaiIyT6BrDUXtL4SH\nBGH6XX08Vr9/LgwRhPcgDQMH+44pnyfC1TTxsEzDl/ay3MCFt+wn4iLlTaz+1LfsJYnEGEv64y6J\nzllHR9+aDCCw4xsQBOE7SMMAOH06XLhm8OrptbKt+oVnNalftH40TyrCwsfScLCkEp8XlKK+NS25\nDR9crJRTPv6Lfth77Kpd6Or7RnZFRGgQhnMsB0y7qzfuTO+MpJgw1Bi8l3qdIIiOAWkYOPh273nJ\nx/hiohXjVsm4MJjj3NVaz/i0FPkNE4G3liRiI0NwV3qKX2kMxDJhWGcs+E2anXFhkFaDCWkpnMmZ\n1CqVZ/NitJIcF+7xc9jh45vXDh8dgvAIJDAAdiNCweFyziKpPZ1d9abf2ZbQxhNGj0Kwx1GXSaZk\nMO836W3X7IER09sundJtBzw0Tah4fit+Hs9V/sjdnrMTIAjCf6ElCQfWfXuCcztnul1PfvmImE/Z\nX+npfROR3q8Cd0nQDHA1n72J1zvDTeZOG+KRel3hStMSiESEBaFXp0gM7ROPT3fId1XlIjSYhg2C\n6IjQmw8ADHClqh7/3nXWaVdEWBDq6ps5D5Myn3pk6mVVGqRV45mHbpV0uK/m0NQecSg+W+WbkyvM\nrb3jJJX3lnpbrVLh5ScyAEBxgcHb+M2SgIj3ZfKYHgjSkuKWCExIYIDFCj07r4hTMLANVhyDhSe1\ny2KEEc8IIc61qqDyq6iQSsOlfRArDAZpNZg8pge+3OksbBKewW8ECA4meyAEPEH4CyQKt8IlLEy5\no6diyw7W5FMPjO4uqny/rsIRC/154CR4YD1P4SEkr8tBiujav2t0u4vXTxD+Co1YAG4Ympy2Bes0\n+NVt3bG18CIA7kFKiiFdWEgQPsiaIPoYDV/2SZnntz9Q2i5PaBeSWuMJ9E6JEi6swOnZfeUvHhPJ\nceF4cGwPfPEjaSc8xUuPprlfiQLPi7cCohGEJyGBAcBbmw47bdMq7HUAKD9oSK3OnxYV4qND8bc/\n3oZofbCvm+JVHG/ZqIFJJDBIhKZegvANsgSGhoYGzJs3D7W1tdDpdMjJyUFiYiLOnz+PV155BSaT\nCTqdDsuXL0dUVBRWrFiBHTt2QKvVYuHChRg8eDCqq6sxf/58NDY2IjExEdnZ2QgODsa2bduwatUq\naLVaTJkyBVOnTgXDMHj11Vdx8uRJ6HQ6vP766+jSpQvOnz+PBQsWQK1Wo0+fPnjllVcU65j0fomW\nHy5Gp0C0vPfmYJwQHcp5/sDrVUJJ6PkgCN8gy4Zh48aNSE1NRV5eHiZNmoQ1a9YAABYtWoQXXngB\nubm5mD59Os6dO4djx46hsLAQmzZtwvLly7FkyRIAwMqVKzFp0iTk5eWhf//+yM/Ph8lkQk5ODtat\nW4fc3Fxs2LABVVVV2Lp1K5qampCfn4958+YhOzsbAJCdnY25c+ciLy8PZrMZW7duVahbgMfu7QtA\nKD20b/GamtPXF6owVjnPF/ksBveyeFWItWUhAoNA/LggOh6yBIYZM2Zg9uzZAIDy8nJERkaisbER\nVVVV+P777/H444/j4MGDGDx4MIqKijB69GgAQHJyMsxmM6qqqnDgwAGMHTsWADBu3Djs2bMHZ86c\nQbdu3aDX6xEUFISMjAzs27cPRUVFtrJDhgxBcbElMVRxcTEyMjLs6lAKx3DNUqYWq7DRLlE5/+yI\nQ52nZLGuSRH4x7w78OBYiyFebGQI9KFBuG9UV8+cMAChJQmC8A2CSxKbN2/G+vXr7bZlZ2cjNTUV\nM2bMQElJCdauXYsbN26gpKQEixcvxgsvvICXX34Zn332GQwGA2Ji2pLhhIeHw2AwwGg0IiIiwrat\nrq7ObhsAhIWFcW7XaDRoaWmxk9qtdQiRkBAhWIZdTtMqOAQHOwdu0rPW3yMi29Tr0VHOqnZX59Vz\nrOOLaadO15YtkK/8TXPbb61GjYSECFjNM0JDdU7HxcXqkRBrCS88emhnHD5zHbcPTsbuI5dFt0uo\njKv9KpVzfAhtkAYDe8TiGCt2g9j72FYvyzshPBgJCRHQBlnurS5Y61RfQnyE7d4LER6m422XVtt2\nj7RBGs52f/KXX4o6j7tI7TMp9URFhUp+NsS2Rx8RbFfWev/EoMQ1h7S++2qNSnZ9YvvHilL3KtCh\nfvIuggJDZmYmMjMzOfetX78epaWleOqpp/D1118jPDwcw4cPBwBMmDABu3btQq9evWAwtCVzMhgM\niIyMtAkOsbGxNoFAr9fblTUajYiKioJer4fRaLRtN5vN0Gg0ULM8CYxGIyIjIwUvuKJCWKhglzO3\nphZubHR2uzSwEvzU1Ta0/a67Kem8Bo5EQWLa2dzcIli+urqt30wtZlRU1Nkm5IaGJqfjqqqMULW0\nICEhAum94/D6H0aisbnFJjCIaZdQGVf7ubQZpuYWmB2MUMXeR1u9bCmktR9MzRZpqqnR5FRfRWWd\nKE+VhIQIGOvbvGwc6zGZ2u6RqblFcruVRKlzc9VTc6NBsP6EhAi7MmLbYzA02pU1Ghslv8fucLP1\n3W9pfW6k0DclCqcu1iBcqxZ9rGM/EdxQP4lHKcFK1pLE6tWr8eWXXwKwaAE0Gg2Cg4PRo0cPFBUV\nAQD279+PPn36YNiwYdi1axcYhkF5eTkYhkF0dDTS0tJQUFAAACgoKEBGRgZ69uyJsrIy1NbWoqmp\nCYWFhRg6dCiGDRuGHTt2AAAOHTqEvn0tKv+BAwdi//79tjrS09Pd6w0X+KNaXq1W4d7hXfCH+wd6\npH6VSuX9RENeoEuSJTV05oReAIBf3d7dqYwn7BvCucKLEwHN3GlDsWTWCKRwpCMniPaGLC+JKVOm\nICsrC5s3bwbDMMjJyQEA/OUvf8GSJUtgNpvRuXNnvPjii9BqtUhPT8e0adPAMAwWL14MAJg9ezay\nsrKwceNGxMTEYNmyZTYvipkzZ4JhGGRmZiIxMRH33HMPdu3ahenTpwOAzegxKysLixYtQnNzM3r1\n6oWJEycq0Sd2uFzL5pEiBnS3JG26J6MLviu8oHib2Ey/S9lEQIHuLh4TEYxBrfdnUPdYrF1wp9fO\nPWNif6+dy+t48Lnx9SPpzvl1QRqkJJCwQAQGsgSGuLg4m2cEm/79++Pjjz922j5nzhzMmTNHVB3j\nx4/H+PHj7bapVCq89tprTmW7d++O3Nxcia3n5oVfD4HZzODtzUfstocFB6EKjQgOEq+MSYwOxeoX\nx4Nh4FmBwU21hz8abns6DPWwPvEeq9sVQ3rFISYigGNOeOCW9UiOxNnLtUiMsbcJ8vZj64evCUH4\nhA4bGrqh0WT3t06rdhqYAODph1IxalASJo9xDi9rN5A4fIZoNeqA/1r3FkoO2KLronvnc1749RD8\n35RbcWtPaQm+CILwDB020uPyjYfs/lapVJxxDW6JDcOTkwbB1GJ22hfGTvPri88Qb4Vh8OK1qQIs\ncpMvYj0ECvrQIAzrk+C03ds9SneQICx0SIHh/NU6nLlUa7dNrZI2tC98LM1JS+GIv2oYXLWLYt4r\nSyBn+SQIomPRIZckrtc4uz2q1NIm+D4pYrJJenjyFTEX+aOdQnsgkMSmzgnh0HggNwpBEB2LDqdh\nKDxxDau+OOq0Xa1SuZQYZH140xhN+AGvzRzR4UITq9DmPqsUHawLCcKJDicwcAkLgPQlCcBhAOE4\n2OPyguQGe6QVikIrIsojJAwHIqtfGm+5boIgFKNDLklwoVK5nn8FxQmOydjj9gAeEABojLWg1L0L\ndKPHHsmWCHJcHka+RKNWkz0OQShMh9MwuITGFz8gwNwkApyXHknDtRsN6BQfeBFBCYKwhzQMrZgZ\nxvXXYAAKE2TBT7hLsE6DLhT2mCA6BB1Kw3Czid8NssXMuKeO94VAEYBCTLckPc6U27u8krFZYPH2\ns2PQ1Owc18R/CcAXjSBk0KE0DJ//cMbu70UzMmy/rVkp+fDLIUPERKpWwJ3OmxP2nCmDvXcywidE\nhOkQFxXi62YQBCGRDiUwVFTX2/3dI7ktHbbZzLg0khI0oPLTr+AkTmM0aUKEN5cuosJ1iIsM4JwL\nBEEQ7ZQOJTCUltc4beudEgUASIj2LytvpVCpVEjv5xxe15/xU9lLFr8Y2QUAcMewzj5uCUEQhHt0\nKBuGMxedBYYXpw9DjbERsZEhqK5rlF852TB4lAnDOqOE4/75O6MG3oKMfonQajqUbE4QRADS4Uex\nIK0a8VHtVLvg7qd4O/qUHzXoFrw39w7Zx08a3R0AMKJ/okItEg8JCwRBBAIdSsNAtG+CdRrZx96T\n0QXjh3ZCkFZ+HUTHhOI/EYSFDvvp89QDg5y2dbR4+0qRFBumaH2O43OQVpnHlIQFgiAI+XRYDcPI\ngUlu1yFGvPjdL/sjJsJDVv8yv3wkHyZwoS89MgwHSyqQt+WUvAbxMKBbDIKDNJh2V29F6yUIgiCk\n02E1DFy4E3u+c7wl2l3/rvZpr8cO7oTUHnGCxz/QusYuCU/kkpBxTExEMO5MS1G8LbGRwXg2czCS\nYpTVYBCEFKxJrJSIaUIQ7ZkOq2FQmi6Jevz59yORGC0vIM2DY3tCrVbhix/PKtyyNqzyxS9GdMHn\nP57F0D7xHjsXQXgML8/bD43riaq6m/j1BNJ0ER0bEhgUpLObCXg0HvqCcax10ugeuCs9BWEhQaKO\nJ8sOwq/w8gMZExGM+dOHefekBOGH0JKEHyF5ScQN+UKssEAIM6hHLADg3uFdfNwSgiAIz0ECA4vQ\nYIsVfbRe5+OWiCRAP/1/eVt3AMBtg27xbUNE0rtzFN55biym3Ukqa69ApgQE4RNoSYJFiE6Lv//x\nNujD6OtbDktmjbAZiLnDhGGdMebWW9qVG6Q+lJ4ZrxGggjJB+DskMDgQ78OcEn774SRygE5J0Ct2\nyvYkLAQ60+/sDY1GjWaTGQO6xfi6OQRB+AgSGDoAv7ytGw6WVGI6xTMgZHDviK6+boI9fitZE0Rg\nQwJDB6D7LZF4/6UJvm4GQRAE0Y7pUEaP0+7u6+sm+D8cNggMLRoTBEF0eDqUwPDAuF6+boJrSNVK\nEARB+CkdSmBQPOscfXgTBEEQHYQOJTAEAlmPej/inNVl0FORKAlCDHGRliRu5MJKEL6hQxk9+vt0\npxLRwn5dY9CzUyRKy2u90CILyXHhePrBVHRPjnC7rmGUv4KQyUuPpuGn4ivtJqAXQQQaHUpgUH5N\nwjcwHlwKCdJwK50y+icqUv/YwZ0UqYfoeCREh2LS6B6+bgZBdFg6lsBA8JL95ChcqapHWAg9EgRB\nEIQzHWp2CAz9gmcUJUmxYUiKDVO+YoIgCCIg6FBGjwGyIkEQBEEQXqdDaRj8HbECjSdtGDo6i2Zk\noMbY5OtmEARB+B0kMBAEix7Jkb5uAkEQhF/SwZYklF2T6JQQDgAY1CNW0XqFoKUVgiAIwtuQhsEN\nEqNDseyZ0YgMp0AyBEEQRGDToQQGT3yYx0QEK1aX2PaRDQNBEAThbTrUkkTA+FUSBEEQhJfpWAKD\nQvxz/h2K1DOgWwwAoFO8tPgHZMNAEARBeBtZSxINDQ2YN28eamtrodPpkJOTg8TEROzevRvLli1D\nUFAQbrvtNjz33HMAgBUrVmDHjh3QarVYuHAhBg8ejOrqasyfPx+NjY1ITExEdnY2goODsW3bNqxa\ntQparRZTpkzB1KlTwTAMXn31VZw8eRI6nQ6vv/46unTpgvPnz2PBggVQq9Xo06cPXnnlFZftVsro\nMUirUaSeOQ/firOXa22CA0kCBEEQhL8iS8OwceNGpKamIi8vD5MmTcKaNWsAAEuXLsXSpUuRn5+P\nvXv3oqSkBMeOHUNhYSE2bdqE5cuXY8mSJQCAlStXYtKkScjLy0P//v2Rn58Pk8mEnJwcrFu3Drm5\nudiwYQOqqqqwdetWNDU1IT8/H/PmzUN2djYAIDs7G3PnzkVeXh7MZjO2bt3qst3+Nh2HBmsxsHus\n4t4bBEEQBKE0sgSGGTNmYPbs2QCA8vJyREZafNcHDhyI6upqNDU1oampCWq1GkVFRRg9ejQAIDk5\nGWazGVVVVThw4ADGjh0LABg3bhz27NmDM2fOoFu3btDr9QgKCkJGRgb27duHoqIiW9khQ4aguLgY\nAFBcXIyMjAy7OtozJDYQBEEQ/orgksTmzZuxfv16u23Z2dlITU3FjBkzUFJSgrVr1wIA+vTpgz/+\n8Y+IiYlBv3790LNnT2zZsgUxMTG2Y8PDw2EwGGA0GhEREWHbVldXZ7cNAMLCwji3azQatLS0gGG5\nC1jrcAX7Qz4hQX6qZneOdYVe3+Zx4eoc1iURnU7rsbYItUEuUVGhitTryeuWij+1xd+hvhIH9ZM4\nqJ+8i6DAkJmZiczMTM5969evR2lpKZ566il89tlnWL16Nb755hskJCRg6dKlWLt2LSIiImA0Gm3H\nGAwGREZG2gSH2NhYm0Cg1+thMBhsZY1GI6KioqDX6+3qMJvN0Gg0UKvVdmWtmg5+2iSGigrXwoUr\n3DnWFQZDo6hzNJtaAABNTSaPtSUhIcIjddfUNChSr6euWyqe6qdAhPpKHNRP4qB+Eo9SgpWsJYnV\nq1fjyy+/BGDRAmg0GgQHByM8PByhoaGtDUxAbW0thg0bhh9//BEMw6C8vBwMwyA6OhppaWkoKCgA\nABQUFCAjIwM9e/ZEWVkZamtr8f/bu/+gqKt/j+MvWBZIlx9q+v3q1Ythmqmlwja375jGLRw1tcla\nE71kXbl3lBlnLNRWTVEMk3HUf8KcqVExmAZd/KM/mkZzKhnNScEhC7M/6IeOOo2KJbsWC+65fzhu\ngNiHC7rrus/HX3A8fvZ83h7Yl+fz2c/x+/2qqanRuHHjNH78eB06dEiSVFdXpxEjRki6cQnk+PHj\nwWNkZmb2rBrhxjUJAMA9qlufknjppZfkdrtVVVUlY4xKSkoUHx8vt9utBQsWKCEhQcnJySopKVFS\nUpKcTqfmzJkjY4wKCwslSfn5+XK73dq7d6/69OmjLVu2BD9FsWDBAhlj5HK5NGDAAE2ePFlHjhxR\nTk6OJAVvenS73VqzZo1aWlo0bNgwTZ069W/Hzb2FAAB0T7cCQ79+/YKfjGgrOztb2dnZt7QvXrxY\nixcv7tIxsrKylJWV1a4tJiZGRUVFt/QdOnSoysvLuzxu8gIAAN3Dg5sAAICl6AoM9/g1iXt7dACA\naBZVgeFef0PmAU4AgHtVVO1W2VP/NXmEEux35rHQAABEkqgKDD39D/yzmYPvzEAAAIgw0XVJgiV/\nAAC6JaoCw73OHsc/BwDg3hRVlyTudf8a/Q/9cOaKns0cEu6hAADQDoHhHmKPs+l/Z44O9zAAALgF\na+ARzFh3AQDgjiAwAAAASwSGCMZnPgAAocI9DBEsGi9J/Ps/HIqP4+FZABBqBIYIFM0rC2tfeyLc\nQwCAqERgiEDRuLJwEw/fAoDw4B6GCMZbJwAgVAgMESyaVxoAAKFFYIhArCwAAEKNwBCBWFkAAIQa\ngSGCsdIAAAgVAkMEY6UBABAqBIYIxMoCACDUCAwRiJUFAECoERgiGCsNAIBQITBEMFYaAAChQmCI\nQKwsAABCjcAQgVhZAACEGoEhgrHSAAAIFQJDBGOlAQAQKlG3vfXAfr00IPWBcA+jRyJxZSGld7x+\n9/mV4ogP91AAAN0QdYGh+H/+QzExkfiW+5dIXFlY+99P6KfzV/XQwORwDwUA0A1Rd0ki0sNCpEp1\nJGj8iP7hHgYAoJuiLjDcD4g8AIBQIzAAAABLBIYIFIn3MAAAIhuBAQAAWCIwRCDuYQAAhBqBAQAA\nWCIwRCDuYQAAhBqBAQAAWCIwRCDuYQAAhBqBAQAAWCIwRCDuYQAAhFqPAkNDQ4OcTqf8fr8kqa6u\nTi+//LLmzZun0tLSYL/S0lLNnj1bc+fO1cmTJyVJV65cUV5ennJzc1VQUKDm5mZJ0ueffy6Xy6Wc\nnBx5PB5JkjFGa9euVU5OjubPn6+zZ89Kks6cOaN58+YpNzdXRUVFPTkVAADwN7odGLxerzZt2qSE\nhIRg27p167R161Z99NFHOnnypE6fPq1Tp06ppqZGHo9HW7du1fr16yVJ27Zt08yZM1VRUaGRI0eq\nsrJSra2tKikpUVlZmcrLy7Vnzx41Njbq4MGD8vv9qqys1NKlS7Vx40ZJ0saNG1VQUKCKigoFAgEd\nPHiwh+WIDNzDAAAItW4HhsLCQhUUFCgxMVHSjQDR0tKiwYMHS5KeeuopHTlyRLW1tZowYYIkaeDA\ngQoEAmpsbNSJEyc0ceJESdKkSZN09OhRNTQ0KC0tTQ6HQ3a7XU6nU8eOHVNtbW2w79ixY1VfXy9J\nqq+vl9PpbHcMAABw58VZdaiqqtLu3bvbtQ0aNEjTp0/XI488ImNuXFH3+XxyOBzBPr1799bZs2eV\nmJio1NTUdu1er1c+n09JSUnBtqampnZtktSrV69O2202m65fvx587bbHiAbcwwAACDXLwOByueRy\nudq1TZkyRVVVVfJ4PLp06ZLy8vK0fft2eb3eYB+fz6eUlBTZ7Xb5fL5gu9frVXJycjA49O3bNxgI\nHA5Hp8dwOBztjhEIBGSz2RQbG9uub3JysuUJ9++fZNnnXmePs0mS4u22u3o+90OtQoE6dR216hrq\n1DXUKbQsA0Nn9u/fH/z6mWee0c6dO2W32xUfH6+zZ89q8ODBOnz4sBYvXiybzabNmzdrwYIFunDh\ngowxSk1NVUZGhqqrq/XCCy+ourpaTqdT6enp+uWXX3T16lUlJiaqpqZGeXl5kqQvvvhCU6dOVV1d\nnUaMGCFJGjVqlI4fP64nnnhC1dXVevLJJy3HfvFi5K9CtLZelyT5W67ftfPp3z/pvqjV3Uaduo5a\ndQ116hrq1HV3Klh1KzC0FRMTE7w0UFRUpGXLlikQCGjChAl6/PHHJUmZmZmaM2eOjDEqLCyUJOXn\n58vtdmvv3r3q06ePtmzZori4OK1cuVILFiyQMUYul0sDBgzQ5MmTdeTIEeXk5EhS8KZHt9utNWvW\nqKWlRcOGDdPUqVN7ejoAAKATMabtjQBR4H5IpBs+rFHD+asa/VBfLZ0z7q68Bum9a6hT11GrrqFO\nXUOduu5OrTDw4CYAAGCJwAAAACwRGAAAgCUCAwAAsERgAAAAlggMAADAEoEBAABYIjBEoLR/3vhM\n7ZABDoueAADcGT1+0iNCb/Z/Pqxh/5aizBH9wz0UAECUIDBEoAS7Tf8a/c9wDwMAEEW4JAEAACwR\nGAAAgCUCAwAAsERgAAAAlggMAADAEoEBAABYIjAAAABLBAYAAGCJwAAAACwRGAAAgCUCAwAAsERg\nAAAAlggMAADAEoEBAABYIjAAAABLBAYAAGCJwAAAACwRGAAAgCUCAwAAsERgAAAAlggMAADAEoEB\nAABYIjAAAABLBAYAAGCJwAAAACwRGAAAgCUCAwAAsERgAAAAlggMAADAEoEBAABYIjAAAABLBAYA\nAGCJwAAAACz1KDA0NDTI6XTK7/dLko4ePaqcnBy98sorWrJkiZqbmyVJpaWlmj17tubOnauTJ09K\nkq5cuaK8vDzl5uaqoKAg2Pfzzz+Xy+VSTk6OPB6PJMkYo7Vr1yonJ0fz58/X2bNnJUlnzpzRvHnz\nlJubq6Kiop6cCgAA+BvdDgxer1ebNm1SQkJCsG39+vV67733VF5errS0NHk8Hp06dUo1NTXyeDza\nunWr1q9fL0natm2bZs6cqYqKCo0cOVKVlZVqbW1VSUmJysrKVF5erj179qixsVEHDx6U3+9XZWWl\nli5dqo0bN0qSNm7cqIKCAlVUVCgQCOjgwYM9LAcAAOhMtwNDYWGhCgoKlJiYGGwrLy9X3759JUmt\nra1KSEhQbW2tJkyYIEkaOHCgAoGAGhsbdeLECU2cOFGSNGnSJB09elQNDQ1KS0uTw+GQ3W6X0+nU\nsWPHVFtbG+w7duxY1dfXS5Lq6+vldDrbHQMAANx5cVYdqqqqtHv37nZtgwYN0vTp0/XII4/IGBNs\nf/DBByVJBw4c0LFjx/T6669rx44d6tOnT7BP79695fV65fP5lJSUFGxrampq1yZJvXr16rTdZrPp\n+vXr7V775jEAAMCdZxkYXC6XXC5Xu7YpU6aoqqpKHo9Hly5dUl5ensrLyyVJZWVlOnDggHbs2KH4\n+Hg5HA75fL7g3/V6vUpOTg4Gh759+wYDgcPhkNfrDfb1+XxKSUm55RiBQEA2m02xsbHt+iYnJ1ue\ncP/+SZZ9cAO16hrq1HXUqmuoU9dQp9Dq1iWJ/fv368MPP1R5ebkefPBB7dy5U5K0fft2nThxQmVl\nZUpJSZEkZWRk6PDhwzLG6Pz58zLGKDU1VRkZGaqurpYkVVdXy+l0Kj09Xb/88ouuXr0qv9+vmpoa\njRs3TuPHj9ehQ4ckSXV1dRoxYoQkadSoUTp+/HjwGJmZmT2rBgAA6FSMabuu3w3PPvusPv30UzU1\nNenpp5/WmDFjZLfbFRMTo+eee045OTkqLS1VdXW1jDFauXKlMjIydPnyZbndbl27dk19+vTRli1b\nlJiYqC+//FKlpaUyxsjlcmnu3LkyxmjdunX64YcfJN242fGhhx7Szz//rDVr1qilpUXDhg1TDIfY\nwwAABQBJREFUcXGxYmJi7khhAADAX3ocGAAAwP2PBzcBAABLBAYAAGCJwAAAACwRGAAAgCXL5zDc\nD9p+yiI+Pl4bNmzQkCFDwj2ssHvxxRflcDgkSYMHD9aiRYu0YsUKxcbGavjw4Vq7dq0kae/evdqz\nZ4/sdrsWLVqkrKysMI46dL755htt3rxZ5eXlOnPmTJdr09zcrOXLl+vy5ctyOBwqKSlp9/Cy+1Hb\nWn3//fdauHChhg4dKkmaO3eupk2bFtW1am1t1apVq3Tu3Dm1tLRo0aJFevjhh5lTHXRWp4EDBzKf\nOhEIBLR69Wr99NNPio2NVVFRkeLj4+/unDJR4MCBA2bFihXGGGPq6upMfn5+mEcUfs3NzWbWrFnt\n2hYtWmSOHz9ujDGmsLDQfPbZZ+bixYtmxowZpqWlxTQ1NZkZM2YYv98fjiGH1AcffGBmzJhh5syZ\nY4z5/9Vm165d5t133zXGGPPJJ5+Y4uLisJ1HKHSs1d69e82uXbva9Yn2Wu3bt8+88847xhhjfv/9\nd5OVlcWc6kTbOv32228mKyvLeDwe5lMnPvvsM7Nq1SpjjDFff/21yc/Pv+tzKiouSXTci+K7774L\n84jC7/Tp07p27Zry8vL02muv6ZtvvtGpU6fa7c3x1Vdf6eTJk8rMzFRcXJwcDoeGDh0afB7G/Swt\nLU3btm0Lft9x35Lb1eb06dOqra3VpEmTgn3v9z1OOqvVl19+qdzcXK1evVo+ny/qazVt2jQtWbJE\nknT9+nXZbLYu/7xFa50CgYDi4uJUX1+vL774gvnUQXZ2tt5++21J0vnz55WSknLX51RUBAav19tu\nL4q4uDgFAoEwjij8EhMTlZeXpx07dmjdunVatmzZLXtzdNzzQ/prf4/73eTJk2Wz2YLfd7U2N9tv\nXuq52fd+1rFWY8eO1ZtvvqmKigoNGTJEpaWlt/wMRlutHnjggeA5L1myRG+88QZzqhMd6/T666/r\n8ccfl9vtZj51IjY2VitWrFBxcbFmzJhx1+dUVASGzvaiaLsPRTQaOnSonn/++eDXqampunz5cvDP\nb+7N0dn+Hl3Zs+N+09m+JberTdv51vGHNRpkZ2dr1KhRwa9Pnz6tpKSkqK/VhQsX9Oqrr2rWrFma\nPn06c+o2OtaJ+fT3SkpKtH//fq1evVrNzc3B9rsxp6LiXTMjI6PTvSii2b59+1RSUiJJ+vXXX+X1\nejVhwgQdO3ZM0l97czz22GOqra2V3+9XU1OTfvzxRw0fPjycQw+LzvYtuV1t2u59cujQoeASYbTI\ny8vTt99+K0k6evSoRo8eHfW1urlJ3/LlyzVr1ixJ0qOPPsqc6qCzOjGfOvfxxx/r/ffflyQlJCQo\nNjZWY8aM6fLv8O7UKioeDW1usxdFNGtpadHKlSt1/vx5xcbGavny5UpNTdXq1atv2ZvD4/Foz549\nMsYoPz9f2dnZ4R5+SJw7d05Lly5VZWXlbfct6aw2f/75p9xuty5evKj4+Hht2bJF/fr1C/fp3FVt\na3Xq1Cm9/fbbstvt6t+/v9avX6/evXtHda02bNigTz/9VOnp6TLGKCYmRm+99ZaKi4uZU210Vqc3\n3nhDmzZtYj518Mcff2jlypW6dOmSWltbtXDhQqWnp3f5d3h3ahUVgQEAAPRMVFySAAAAPUNgAAAA\nlggMAADAEoEBAABYIjAAAABLBAYAAGCJwAAAACz9H4w46gW6AC/8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11303e320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def run_advi():\n",
    "    with model:\n",
    "        v_params = pm.variational.advi_minibatch(\n",
    "            n=3000, minibatch_tensors=minibatch_tensors, minibatches=minibatches, \n",
    "            local_RVs=local_RVs, observed_RVs=observed_RVs, encoder_params=encoder_params, \n",
    "            learning_rate=2e-2, epsilon=0.1, n_mcsamples=1 \n",
    "        )\n",
    "        \n",
    "    return v_params\n",
    "\n",
    "%time v_params = run_advi()\n",
    "plt.plot(v_params.elbo_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see ELBO increases as optimization proceeds. The trace of ELBO looks jaggy because at each iteration documents in the mini-batch are replaced. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extraction of characteristic words of topics based on posterior samples\n",
    "By using estimated variational parameters, we can draw samples from the variational posterior. To do this, we use function `sample_vp()`. Here we use this function to obtain posterior mean of the word-topic distribution $\\beta$ and show top-10 words frequently appeared in the 10 topics. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Topic #0: don think just know make going like people want sure\n",
      "Topic #1: year team game play win games players season period new\n",
      "Topic #2: edu information com mail list send available university 1993 email\n",
      "Topic #3: people state government gun world said years war states armenian\n",
      "Topic #4: god people believe does true jesus say question life way\n",
      "Topic #5: windows use thanks drive using window card file does work\n",
      "Topic #6: key use chip encryption government public keys used law clipper\n",
      "Topic #7: did time didn got said just day right thought let\n",
      "Topic #8: good like ve better really car probably lot know problem\n",
      "Topic #9: new power space 10 00 years used price 50 high\n"
     ]
    }
   ],
   "source": [
    "def print_top_words(beta, feature_names, n_top_words=10):\n",
    "    for i in range(len(beta)):\n",
    "        print((\"Topic #%d: \" % i) + \" \".join([feature_names[j]\n",
    "            for j in beta[i].argsort()[:-n_top_words - 1:-1]]))\n",
    "\n",
    "doc_t.set_value(docs_te.toarray().astype('float32')[:minibatch_size, :])\n",
    "\n",
    "with model:\n",
    "    samples = sample_vp(v_params, draws=100, local_RVs=local_RVs)\n",
    "    beta_pymc3 = samples['beta'].mean(axis=0)\n",
    "\n",
    "print_top_words(beta_pymc3, feature_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare these topics to those obtained by a standard LDA implementation on scikit-learn, which is based on an online stochastic variational inference (Hoffman et al., 2013). We can see that estimated words in the topics are qualitatively similar. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 11.8 s, sys: 11.7 ms, total: 11.8 s\n",
      "Wall time: 11.8 s\n",
      "Topic #0: government people law mr president gun state states public rights\n",
      "Topic #1: drive card scsi bit disk use mac hard memory does\n",
      "Topic #2: people armenian said armenians turkish did war killed saw russian\n",
      "Topic #3: year just good time team game car years like think\n",
      "Topic #4: 10 00 25 15 20 12 11 14 16 17\n",
      "Topic #5: windows window program dos file use version display ms application\n",
      "Topic #6: edu space com file information mail data available send ftp\n",
      "Topic #7: ax max g9v pl b8f a86 34u 145 1t 75u\n",
      "Topic #8: god people jesus believe does think say life don true\n",
      "Topic #9: don like know just think ve does use want good\n"
     ]
    }
   ],
   "source": [
    "from sklearn.decomposition import LatentDirichletAllocation\n",
    "\n",
    "lda = LatentDirichletAllocation(n_topics=n_topics, max_iter=5,\n",
    "                                learning_method='online', learning_offset=50.,\n",
    "                                random_state=0)\n",
    "%time lda.fit(docs_tr)\n",
    "beta_sklearn = lda.components_ / lda.components_.sum(axis=1)[:, np.newaxis]\n",
    "\n",
    "print_top_words(beta_sklearn, feature_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predictive distribution\n",
    "In some papers (e.g., Hoffman et al. 2013), the predictive distribution of held-out words was proposed as a quantitative measure for goodness of the model fitness. The log-likelihood function for tokens of the held-out word can be calculated with posterior means of $\\theta$ and $\\beta$. The validity of this is explained in (Hoffman et al. 2013). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def calc_pp(ws, thetas, beta, wix):\n",
    "    \"\"\"\n",
    "    Parameters\n",
    "    ----------\n",
    "    ws: ndarray (N,)\n",
    "        Number of times the held-out word appeared in N documents. \n",
    "    thetas: ndarray, shape=(N, K)\n",
    "        Topic distributions for N documents.         \n",
    "    beta: ndarray, shape=(K, V)\n",
    "        Word distributions for K topics. \n",
    "    wix: int\n",
    "        Index of the held-out word\n",
    "        \n",
    "    Return\n",
    "    ------\n",
    "    Log probability of held-out words.\n",
    "    \"\"\"\n",
    "    return ws * np.log(thetas.dot(beta[:, wix]))\n",
    "\n",
    "def eval_lda(transform, beta, docs_te, wixs):\n",
    "    \"\"\"Evaluate LDA model by log predictive probability. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    transform: Python function\n",
    "        Transform document vectors to posterior mean of topic proportions. \n",
    "    wixs: iterable of int\n",
    "        Word indices to be held-out. \n",
    "    \"\"\"\n",
    "    lpss = []\n",
    "    docs_ = deepcopy(docs_te)\n",
    "    thetass = []\n",
    "    wss = []\n",
    "    total_words = 0\n",
    "    for wix in wixs:\n",
    "        ws = docs_te[:, wix].ravel()\n",
    "        if 0 < ws.sum():\n",
    "            # Hold-out\n",
    "            docs_[:, wix] = 0\n",
    "            \n",
    "            # Topic distributions\n",
    "            thetas = transform(docs_)\n",
    "        \n",
    "            # Predictive log probability\n",
    "            lpss.append(calc_pp(ws, thetas, beta, wix))\n",
    "            \n",
    "            docs_[:, wix] = ws\n",
    "            thetass.append(thetas)\n",
    "            wss.append(ws)\n",
    "            total_words += ws.sum()\n",
    "        else:\n",
    "            thetass.append(None)\n",
    "            wss.append(None)\n",
    "    \n",
    "    # Log-probability\n",
    "    lp = np.sum(np.hstack(lpss)) / total_words\n",
    "    \n",
    "    return {\n",
    "        'lp': lp, \n",
    "        'thetass': thetass, \n",
    "        'beta': beta, \n",
    "        'wss': wss\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To apply the above function for the LDA model, we redefine the probabilistic model because the number of documents to be tested changes. Since variational parameters have already been obtained, we can reuse them for sampling from the approximate posterior distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "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_docs_te = docs_te.shape[0]\n",
    "doc_t = shared(docs_te.toarray().astype('float32'), name='doc_t')\n",
    "\n",
    "with pm.Model() as model:\n",
    "    theta = Dirichlet('theta', a=(1.0 / n_topics) * np.ones((n_docs_te, n_topics)).astype('float32'), \n",
    "                      shape=(n_docs_te, 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)\n",
    "\n",
    "# Encoder has already been trained\n",
    "encoder.p_corruption = 0\n",
    "local_RVs = OrderedDict([(theta, (encoder.encode(doc_t), 1))])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`transform()` function is defined with `sample_vp()` function. This function is an argument to the function for calculating log predictive probabilities. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def transform_pymc3(docs):\n",
    "    with model:\n",
    "        doc_t.set_value(docs)\n",
    "        samples = sample_vp(v_params, draws=100, local_RVs=local_RVs)\n",
    "    \n",
    "    return samples['theta'].mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of the log predictive probability is -7.00. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 1s, sys: 3.76 s, total: 4min 4s\n",
      "Wall time: 1min 45s\n",
      "Predictive log prob (pm3) = -7.104652449902257\n"
     ]
    }
   ],
   "source": [
    "%time result_pymc3 = eval_lda(transform_pymc3, beta_pymc3, docs_te.toarray().astype('float32'), np.arange(100))\n",
    "print('Predictive log prob (pm3) = {}'.format(result_pymc3['lp']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare the result with the scikit-learn LDA implemented The log predictive probability is significantly higher (-6.04) than AEVB-ADVI, though it shows similar words in the estimated topics. It may because that the mean-field approximation to distribution on the simplex (topic and/or word distributions) is less accurate. See https://gist.github.com/taku-y/f724392bc0ad633deac45ffa135414d3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 25.7 s, sys: 41.8 ms, total: 25.7 s\n",
      "Wall time: 25.7 s\n",
      "Predictive log prob (sklearn) = -6.041000209921997\n"
     ]
    }
   ],
   "source": [
    "def transform_sklearn(docs):\n",
    "    thetas = lda.transform(docs)\n",
    "    return thetas / thetas.sum(axis=1)[:, np.newaxis]\n",
    "\n",
    "%time result_sklearn = eval_lda(transform_sklearn, beta_sklearn, docs_te.toarray(), np.arange(100))\n",
    "print('Predictive log prob (sklearn) = {}'.format(result_sklearn['lp']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "We have seen that PyMC3 allows us to estimate random variables of LDA, a probabilistic model with latent variables, based on automatic variational inference. Variational parameters of the local latent variables in the probabilistic model are encoded from observations. The parameters of the encoding model, MLP in this example, are optimized with variational parameters of the global latent variables. Once the probabilistic and the encoding models are defined, parameter optimization is done just by invoking a function (`advi_minibatch()`) without need to derive complex update equations. \n",
    "\n",
    "Unfortunately, the estimation result was not accurate compared to LDA in sklearn, which is based on the conjugate priors and thus not relying on the mean field approximation. To improve the estimation accuracy, some researchers proposed post processings that moves Monte Carlo samples to improve variational lower bound (e.g., Rezende and Mohamed, 2015; Salinams et al., 2015). By implementing such methods on PyMC3, we may achieve more accurate estimation while automated as shown in this notebook. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# References\n",
    "* Kingma, D. P., & Welling, M. (2014). Auto-Encoding Variational Bayes. stat, 1050, 1.\n",
    "* Kucukelbir, A., Ranganath, R., Gelman, A., & Blei, D. (2015). Automatic variational inference in Stan. In Advances in neural information processing systems (pp. 568-576).\n",
    "* Blei, D. M., Ng, A. Y., & Jordan, M. I. (2003). Latent dirichlet allocation. Journal of machine Learning research, 3(Jan), 993-1022.\n",
    "* Hoffman, M. D., Blei, D. M., Wang, C., & Paisley, J. W. (2013). Stochastic variational inference. Journal of Machine Learning Research, 14(1), 1303-1347.\n",
    "* Rezende, D. J., & Mohamed, S. (2015). Variational inference with normalizing flows. arXiv preprint arXiv:1505.05770.\n",
    "* Salimans, T., Kingma, D. P., & Welling, M. (2015). Markov chain Monte Carlo and variational inference: Bridging the gap. In International Conference on Machine Learning (pp. 1218-1226)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  },
  "widgets": {
   "state": {},
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

@taku-y Thanks for sharing. I have not tried your code but reading it I think there is bug here:

    def encode(self, xs):
        if 0 < self.p_corruption:
            dixs, vixs = xs.nonzero()
            mask = tt.set_subtensor(
                tt.zeros_like(xs)[dixs, vixs], 
                self.rng.binomial(size=dixs.shape, n=1, p=1-self.p_corruption)
            )
            xs_ = xs * mask
        else:
            xs_ = xs

        w0 = self.w0.reshape((self.n_words, self.n_hidden))
        w1 = self.w1.reshape((self.n_hidden, 2 * (self.n_topics - 1)))
        hs = tt.tanh(xs_.dot(w0) + self.b0)
        zs = hs.dot(w1) + self.b1
        zs_mean = zs[:, :(self.n_topics - 1)]
        zs_std = zs[:, (self.n_topics - 1):]
        return zs_mean, zs_std

don't you need to enforce non-negativity for zs_std ? I might be wrong or have not carefully read the rest of the code.

Thanks for your comment. zs_std is actually log transformed std and tt.exp(zs_std) is used as std in advi_minibatch(), the inference program. I will add a comment to the notebook about this.

@taku-y Thanks for sharing this. I've been working on this problem for a while myself, so it is a relief to see it work.
but when i tried the code I couldn't replicate the result. I had to make the following edits in order to make it work, which might be causing the problem:

in block 8 above, minibatches = [create_minibatch(docs_tr.toarray().astype('float32'))] had to remove the braces otherwise _check_minibatches(minibatch_tensors, minibatches) was raising
_value_error(isinstance(minibatches, types.GeneratorType), 'minibatches must be a generator.')

had to change e = f(*next(minibatches)) to e = f(next(minibatches)) in line 450 in advi_minibatch.py, otherwise (after 1.) too many arguments was being raised

most importantly, had to change the learning rate. i tried values in the range [1e-2,1e-4] since keeping it to 5e-2 was resulting in elbo throwing nan after ~40 iterations. for example , the plot is for 1e-3 for 4000 iterations. when i run with 1e-2 and sometimes when it finishes without nans the result is slightly better but not even close to yours. for me 4e-3 with roughly 5000 iterations works the best but again not as good as in your case. seems like its too sensitive to the optimiser and the optimisation with regard to the learning rate (what algo is being used behind the scene, is it adam? )

As such my runs never reach the elbo values that are reported here. and consequently my topics are way worse:
`Topic #0: like just don make does know people new think did

Topic #1: people like just does new use good time know think

Topic #2: like just people know does don good new use think

Topic #3: people know just like don does good time use god

Topic #4: people know time like use don does just good god

Topic #5: just like don time think new say know people good

Topic #6: don just time know good think use people like make

Topic #7: just people know edu use like good don does think

Topic #8: people new does time just don good use like know

Topic #9: like don just new good people know does time use`

is there anything that pops out to you about my changes that might be causing this behaviour? in particular, it seems like most of the topics are quite similar in my case.

@akashgit Sorry for late to reply. I had been working on incorporation of autoencoding VB into the PyMC3 repo.

For 1. and 2. of your comments, I fixed the bug on the notebook. Your comments were helpful for the fix. Thanks.

1. The sensitivity might be due to the LDA model and/or optimization algo (adagrad with the gradients of the latest 10 updates), although it's not clear to me. I also encountered the same problem that ELBO was nan with learning rate of 5e-2. So in the notebook I set the learning rate to 2e-2 and n=3000. Try to run the updated notebook with the latest version of PyMC3. I will continue to examine this sensitivity issue.

hi @taku-y, I tried using this notebook after updating to the latest pymc3 version (instead of running it from your fork) and I keep getting test_compute_value related error on line doc = pm.DensityDist('doc', logp_lda_doc(beta, theta), observed=doc_t) about failure to downcast from float64 to 32 (my global floatX is set to 32)
I've double checked my .theanorc and even tried manually setting the flag to ignore/off but cannot seem to resolve it. any suggestions? just to clarify, this problem does not occur when using your old fork from my last comment with the edits but I am still struggling with setting the learning rate and had not been able to reproduce the topics of similar quality to yours.

Hi, thanks for your report. What version of theano and numpy are you using? I use theano 0.8.2 and numpy 1.10.4. Updating these libraries to latest versions might resolve the problem. Actually, I'm using a Docker container. If you want, I prepare a Docker file to build a container image on which the notebook works.

Hi @taku-y , I was using theano 0.8.2 and 1.10.4 before updating to latest pymc3 (when things were working). After upgrading, while debugging through the error i ended up on 0.9.0.dev2 and numpy 1.11.1 (though it wasn't helpful).

i found these two related issues mentioned on pymc3 repo
pymc-devs/pymc#1246
pymc-devs/pymc#1253
though not sure if they are resolved.

The interesting thing is that these problem only occur when i update to the latest pymc3 version. Running them on your old fork (from the time of my first comment with the suggested edits) doesn't cause any of those errors.

i haven't used docker before, but m willing to try it so it would be nice if you could send me the image. Thanks again.

I could not reproduce the test_compute_value error with the latest master branch of PyMC3. Could you show the error message? We might find something wrong.

Hi @taku-y, Here is the actual error. Somehow magically float64 are appearing even though my config file and flags are appropriately set.

Applied stickbreaking-transform to theta and added transformed theta_stickbreaking_ to model.

Applied stickbreaking-transform to beta and added transformed beta_stickbreaking_ to model.

---

TypeError Traceback (most recent call last)
in ()
10 beta = Dirichlet('beta', a=(1.0 / n_topics) * np.ones((n_topics, n_words)).astype('float32'),
11 shape=(n_topics, n_words), transform=t_stick_breaking(1e-9))
---> 12 doc = pm.DensityDist('doc', logp_lda_doc(beta, theta), observed=doc_t)

/IPC_MAP/lib/python2.7/site-packages/pymc3-3.0-py2.7.egg/pymc3/distributions/distribution.pyc in new(cls, name, _args, *_kwargs)
24 data = kwargs.pop('observed', None)
25 dist = cls.dist(_args, *_kwargs)
---> 26 return model.Var(name, dist, data)
27 elif name is None:
28 return object.new(cls) # for pickle

/IPC_MAP/lib/python2.7/site-packages/pymc3-3.0-py2.7.egg/pymc3/model.pyc in Var(self, name, dist, data)
304 self.named_vars[v.name] = v
305 else:
--> 306 var = ObservedRV(name=name, data=data, distribution=dist, model=self)
307 self.observed_RVs.append(var)
308 if var.missing_values:

/IPC_MAP/lib/python2.7/site-packages/pymc3-3.0-py2.7.egg/pymc3/model.pyc in init(self, type, owner, index, name, data, distribution, model)
579 self.missing_values = data.missing_values
580
--> 581 self.logp_elemwiset = distribution.logp(data)
582 self.model = model
583 self.distribution = distribution

in ll_docs_f(docs)
17 vfreqs = docs[dixs, vixs]
18 ll_docs = vfreqs * pm.math.logsumexp(
---> 19 tt.log(theta[dixs].astype('float32')) + 20 tt.log(beta.T[vixs].astype('float32')), axis=1).ravel().astype('float32')
21

/IPC_MAP/lib/python2.7/site-packages/theano/tensor/var.pyc in getitem(self, args)
502 TensorVariable, TensorConstant,
503 theano.tensor.sharedvar.TensorSharedVariable))):
--> 504 return self.take(args[axis], axis)
505 else:
506 return theano.tensor.subtensor.advanced_subtensor(self, *args)

/IPC_MAP/lib/python2.7/site-packages/theano/tensor/var.pyc in take(self, indices, axis, mode)
546
547 def take(self, indices, axis=None, mode='raise'):
--> 548 return theano.tensor.subtensor.take(self, indices, axis, mode)
549
550 # COPYING

/IPC_MAP/lib/python2.7/site-packages/theano/tensor/subtensor.pyc in take(a, indices, axis, mode)
2369 return advanced_subtensor1(a.flatten(), indices)
2370 elif axis == 0:
-> 2371 return advanced_subtensor1(a, indices)
2372 else:
2373 if axis < 0:

/IPC_MAP/lib/python2.7/site-packages/theano/gof/op.pyc in call(self, _inputs, *_kwargs)
610 for i, ins in enumerate(node.inputs):
611 try:
--> 612 storage_map[ins] = [self._get_test_value(ins)]
613 compute_map[ins] = [True]
614 except AttributeError:

/IPC_MAP/lib/python2.7/site-packages/theano/gof/op.pyc in _get_test_value(cls, v)
547 # ensure that the test value is correct
548 try:
--> 549 ret = v.type.filter(v.tag.test_value)
550 except Exception as e:
551 # Better error message.

/IPC_MAP/lib/python2.7/site-packages/theano/tensor/type.pyc in filter(self, data, strict, allow_downcast)
138 '"function".'
139 % (self, data.dtype, self.dtype))
--> 140 raise TypeError(err_msg, data)
141 elif (allow_downcast is None and
142 type(data) is float and

TypeError: For compute_test_value, one input test value does not have the requested type.

The error when converting the test value to that variable type:

TensorType(float32, matrix) cannot store a value of dtype float64 without risking loss of precision. If you do not mind this loss, you can: 1) explicitly cast your data to float32, or 2) set "allow_input_downcast=True" when calling "function".
[[ 0.09999998 0.1 0.1 ..., 0.1 0.10000001 0.10000001] [ 0.09999998 0.1 0.1 ..., 0.1 0.10000001 0.10000001] [ 0.09999998 0.1 0.1 ..., 0.1 0.10000001 0.10000001] ..., [ 0.09999998 0.1 0.1 ..., 0.1 0.10000001 0.10000001] [ 0.09999998 0.1 0.1 ..., 0.1 0.10000001 0.10000001] [ 0.09999998 0.1 0.1 ..., 0.1 0.10000001 0.10000001]]

Hi @akashgit, thanks for sending me detailed information. But I could not reproduce your result yet. I tested the notebook on theano-0.8.2 and 0.9.0dev2, python2 and 3. Under all of these environments, ADVI worked.

Here is my Dockerfile on which the notebook ran.

From jupyter/datascience-notebook
MAINTAINER Taku Yoshioka <contact@google.com>

ENV TERM xterm-color

USER jovyan
RUN pip install theano joblib

Save the above as Dockerfile, then you can build an image by typing

docker build -t pymc3-test .

pymc3-test is the name of the created image. Finally, start running a container with the above image

docker run -d -v /your/home:/home/jovyan/work -p 8889:8888 --name pymc3 pymc3-test

You can access the notebook server via localhost:8889. When using VirtualBox (DockerMachine), you need to specify the IP address on which DockerMachine is running. See https://docs.docker.com/machine/reference/env/. I recommend Linux or Mac, both of which have native support of Docker so you don't need to specify the IP address of the DockerMachine.

Here is the option of the command docker run:

• -d: Background execution
• -v: Mount volume /your/home on your host to /home/jovyan/work on the container.
• -p: Port forwarding
• --name: Image of the container

thanks for the file @taku-y. I will try to see what is wrong with my current configuration.