vrnn_tf_v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "from functools import reduce\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "GLOBAL_SEED=1\n",
    "VERBOSITY=1\n",
    "TESTING=0\n",
    "\n",
    "np.random.seed(GLOBAL_SEED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# useful function for recovering functions in deleted cells\n",
    "def rescue_code(function):\n",
    "    import inspect\n",
    "    get_ipython().set_next_input(\"\".join(inspect.getsourcelines(function)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# wrapper around tf.get_variable that checks if the variable has already been defined.\n",
    "# useful for playing around in jupyter. allows you to call cells multiple times\n",
    "def get_variable_wrap(*args, **kwargs):\n",
    "    try:\n",
    "        return tf.get_variable(*args, **kwargs)\n",
    "    except ValueError:\n",
    "        tf.get_variable_scope().reuse_variables()\n",
    "        return tf.get_variable(*args, **kwargs)\n",
    "\n",
    "## do matrix multiplication with weights and add bias (as in fully connected layer)\n",
    "def fc_act(x, next_layer_size, act=None, name=\"fc\"):\n",
    "    nbatches = x.get_shape()[0]\n",
    "    prev_layer_size = x.get_shape()[1]\n",
    "    with tf.name_scope(\"fc\"):\n",
    "        w = get_variable_wrap(\"weights\", [prev_layer_size, next_layer_size], dtype=tf.float32, initializer=tf.random_normal_initializer())\n",
    "        b = get_variable_wrap(\"bias\", [next_layer_size], dtype=tf.float32, initializer=tf.constant_initializer(0.0))\n",
    "        o = tf.add(tf.matmul(x, w), b)\n",
    "        if act: return act(o)\n",
    "        else: return o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sample Data\n",
    "I want to create 100 bi-variate time-series sampled from the same distribution and split them into train/validation/test 70/20/10, i.e. split into groups of time-series size 14/4/2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gen_time_series(n_points):\n",
    "    time_series = np.array([[np.random.normal(mu*i, sigma*i), \n",
    "                        np.random.normal(mu*i, sigma*i/4)] \n",
    "                       for i in range(n_points)])\n",
    "    return time_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 100, 2)\n",
      "20000 data values\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEFCAYAAAAc33cJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXmcHFd59/utpdfpnn0kjVbLlnwk71hgDMZgCCE3QIBL\nEpKXQEjIQt4kF7LcNySBLC+XvGSD3JAE8sYECNyQsBPsAGYJ2MbGm2xsy5aOFksaSaPZu6e7p/eq\nun/U0tUzPaMezabRnO/n4497uqurTtWMfvXU7zzPczTHcVAoFArF5YG+1gNQKBQKxfKhRF2hUCgu\nI5SoKxQKxWWEEnWFQqG4jFCirlAoFJcRStQVCoXiMkKJ+gZBCHGrEOK7QoinhBCHhBBfF0Jcu0rH\ndoQQ/YvY/g4hxKF5PnufEOLnl290Tft+nRDiw21sd6cQ4sAi9/0CIcQ/eq+fL4T4wsWOc6URQnxM\nCPHKtR6H4uIw13oAipVHCBED7gZeJaV83HvvLcDXhRC7pZTWmg5wEUgp/3gF9/1V4KttbPqjwP9e\n5O6vBbZ7x3kM+KlFfn/VkFL+8lqPQXHxKFHfGCSBbiAVeu9fgRxgCCEc4G+AW4E0oAG/LKV8QAjx\nSaAEvADYAnwOGAd+wvv5l6WU/+Vt5wD7gQHgm8A7pZS18ECEEL8E/DruU+Ik8JtSyiMtxpzyotk9\nQBb4VSnlUe84h6SUfz1rv78KvE5K+Vrv533Ad4CdwNuAdwBRoBf4cynlR4UQvwD8EtABTAP/AvyU\nlPK1Qohbgb8EYsAg8C0p5S8JIf4M2Ar8q/fEcAT4W+B6IOId839IKeuhse0A3gd0CSE+4R3n76WU\n1y3i+kaBvwBeBhjAE971zc26DvuAfwbi3u/xY1LKj3ifvQf4Se/anwJ+XUo5LIT4HjAF7AM+6m3z\n91LKLwghXuwdtwOwgT+VUt4thNgCfArwn8D+U0r5Ry1+j4pVRtkvGwApZQb4PeAbQojnhBCfBn4R\n+LaUsgq8EFeoXiSlvAZXdH4/tIvnAS8Cng/8NlCQUr4YV8zC290IvBK4xvvvHeFxCCFehiuwt0sp\nn4crml+aZ9g7gA9JKW8CPgN8+gKn+W/ASzyxwTu/TwAJ4FeAV3vH/BnvuD7XAndIKV8+a3/vAv5Y\nSvlC71xeJ4Q4IKV8DzAM/JyU8mHcm+FBKeUB3OvUD/xOeEdSyjPAHwP3Syl/scXY27m+vw/UgQNS\nyhu9Mfx5i339D+AubzyvBl4qhNC9G9D1wC3eNf0a8LHQ9zJSymuklH/nvyGE6PGu4VullDcDrwM+\nKoTY6V3T57z3bwf2CiG6WoxHscqoSH2DIKX8kBDiTtxI76XAu4F3CyFukVL+QAjxXuAdQoirgDuA\nfOjrd3kR94gQYgb4hvf+CdzI1+eTUsoCgBDiU8AbgL8Pff4a3Mj7QSGE/16vEKJXSjk1a8hPSSkf\n9PeLKybzioaUMu9F9m8RQvwN8BbgJVLKghDitcBrhBB7gZtofmJ5ana06/E24NVCiD/EjWCTs77n\n81rgFu8JBNybyGJp5/q+Fvdp60e9axcFxlrs68vAp4QQtwDfxo3mbe8a3AI85n3f8M7J5/4W+3oR\n7lPKV0K/Lwe4wRvj1zyB/zbw+1LK6cWeuGL5UaK+ARBC3Aa8WEr5V7je+t2eWD2NKxIl3Kjwg8B/\n4FoKbwntojJrlzVaUw+91oHZXr0BfFpK+W5vXDruE0Kmxb5mf9cJH1cI8TpcSwNgWEr5atzI85+A\nw8CzUsqTQojtwA+8978PfAFXIH0K85zL/cCTuOL1OdynGa3Fdgbw01LKw964ur2xLoZ2rq8BvEtK\n+XXvOClci6UJzxrZi+v7/wjwJ56FYgB/IaX8qPf9GNAT+mqr62AAh72nFbzvbQXGpZQ1IcRu3Cez\nVwCPCCHeELoRK9YIZb9sDMaB9wohXhJ6bxDXJ30aVwDu8v7BP4obYRsXcZyfEULEhBBx3Ej3rlmf\nfxP4b0KIQe/nX8P1oFtxoxDiJu/1O4DvSymL/odSyq9KKW/y/nu1995DuML7x8Cd3qbPxz3/90sp\n78ETdCHEvOfn2Q7PB94tpfwSsA33CcP/Th3XPwe4B/htIYTmCeVXgd9ssdvwdy6Ge4DfFEJEvZvh\nncAHWoz9M8DPSCn/HXfuIodrZd0D/LIQotPb9H1c2NJ6CNdWeam375uAY8BWIcSfA38kpfwKrlX1\nDHD1Es5PsUwoUd8ASCmP4gr1//I89Wdxo89flVJK4B+BlwkhnsKNak8Auz3xWAxF3Aj3ae//n5g1\njntwJ92+5R3rzcAbpZStItvDuFHmk7he7tvaHMOdwJXAV7yfvwmcBaQQ4gncidNxXJFuiTcH8QHg\ncSHEY8AfAA+EvvMV4LNCiFcB76Rxc3zK+/9fztmpe133CSG+3OZ5zOb/wZ3cfAJ4Fvfm9bvzbPdz\n3nV7GNeOuRf3KeZu4CEhxDO4FsovLHRAKeU47qTpX3n7+zSuv34a+H+Bm7zU08eAk7jzGoo1RlOt\ndxXLwXxZKQqFYnVRkbpCoVBcRqhIXaFQKC4jVKSuUCgUlxFK1BUKheIyYk3z1MfH80vyfnp6kmQy\nxQtveBmxEc8ZNuZ5q3PeOCz2vAcG0q1qJoB1Hqmb5sWkUq9vNuI5w8Y8b3XOG4flPO91LeoKhUKh\naEaJukKhUFxGKFFXKBSKywgl6gqFQnEZoURdoVAoLiOUqCsUCsVlhBJ1hUKhuIxQoq5QKC477n9q\nmCeOja/1MNaEC1aUeosJ3AkI3BVdfg0o4y4x5gCHgN/wlsz6FdwFDeq4ixLcvULjVigUinn57HeO\n09sZ53l7B9Z6KKtOO5H6TwBIKW8D3gv8GfAh4L1Syttxm/W/3lvw953AbcCPAR/wVoJRKBSKVaVm\n2dQse62HsSZcMFKXUn5FCOFH3LuALO66hPd6730deBXumpIPSCkrQEUIcRx3dZVH59t3T09yyeWx\nAwPpJX1/PbIRzxk25nmrc744LNvBcZx1df2Wa6xtNfSSUtaFEP8C/J/ATwE/GlqCLA90AZ1AeDVx\n//15WWrjnoGBNOPj+QtveBmxEc8ZNuZ5q3O+OBzHwbYdKlVr3Vy/xZ73QjeAtidKpZRvw11Y9k4g\nEfoojRu957zXs99XKBSKVcOy3XizvkHtlwuKuhDirUKIP/B+LAI28JgQ4g7vvR/HXWT4EeB2IURc\nCNEF7MedRFUoFIpVw7J8Ud+Yq7q1Y798CfiEEOI+IAL8Fu5K73cKIaLe6y9IKS0hxIdxBV4H3iOl\nLK/QuBUKhaIldduN0DdqpN7OROkM8KYWH72sxbZ34tozCoVCsSb4kbplO9iOg67Nu57EZYkqPlIo\nFJcVvqcOYG3AaF2JukKhuKwIC3mtvvF8dSXqCoXisiIcqfv++kZCibpCobisqIdFva5EXaFQKNY1\nYftlI2bAKFFXKBSXFU32yyrnqs+Ua7znzod45PDoqh43jBJ1hUJxWWFZYVFf3Uj97FiB85NFjgyt\nXTG9EnWFQnFZYYUmR1e7U2OhVAegUq2v6nHDKFFXKBSXFfWmPPXVtV8KpSoA5aq1qscNo0RdoVBc\nVjTlqa96pF4DoFJToq5QKBTLQpOnvsopjTOB/aJEXaFQKBZF3bJ5z50P8Y2Hh5reX8vsl7xvv6hI\nXaFQKBZHtlDh/GSRE+emm94PV5GudvaLitQVCoXiIilXXOGcLdxrmdLoe+pqolShUCgWiS+c4WwX\naLZf1ESpQqFQrBNKXi747Pa64Z9XP6XRFfVa3W7Kl19NlKgrFIp1SaniivrsydBw5F5bxewX23GY\nKdeCnytVJeoKhULRNoH9col46sVyHSd0f1krC0aJukKhWJfMF6lba5T94lsvPuU1ahVwwTVKhRAR\n4OPAFUAMeD9wBrgbOOZt9lEp5WeFEL8CvAOoA++XUt69EoNWKBQKP1Kf7V03R+qr56nPFvW1itQv\nKOrAW4BJKeVbhRC9wA+B9wEfklJ+0N9ICLEFeCfwfCAOfF8I8S0pZWUFxq1QKDY4jUi9WdSbFslY\ng0g9aupU6/aa5aq3I+qfB77gvdZwo/ADgBBCvB43Wv8t4BbgAU/EK0KI48ANwKPLPmqFQrHh8e2N\nS8V+mfFEvbczzshUcc1y1S8o6lLKAoAQIo0r7u/FtWE+JqU8KIR4D/AnuBF8uLQrD3QttO+eniSm\naVzk0F0GBtJL+v56ZCOeM2zM81bnPD82GgCO0/ydWCwSvDYj5qpdQ0d3F8bYOpBiZKpILB5d1LGX\na5ztROoIIXYAXwY+IqX8jBCiW0rpd4H/MvB3wH1AeFRpYMFO8ZlMcfEjDjEwkGZ8PL+kfaw3NuI5\nw8Y8b3XOCzOdd53das1q+k6+UGl6vVrXcHSiAEA64crq2GSh7WMv9ne90A3ggtkvQojNwDeBd0sp\nP+69fY8Q4hbv9Y8AB4FHgNuFEHEhRBewHzjU9igVCoViEfjFR/U5E6Wh4qNVLADyPfW+zjiwdv1f\n2onU/xDoAf5ICPFH3nu/A/yNEKIGjAC/KqXMCSE+DNyPe7N4j5SyvBKDVigUinLFryi9NIqPCkVP\n1LtcUV+rTo3teOrvAt7V4qPbWmx7J3DnMoxLoVAoFqRU8VMaHWzHQddcj30tUxo1oDcdA9YuUlfF\nRwqFYl0SLu4JC/maFR+VayTjJsm4O1GrRF2hUCjaxHGcppTBepOPvkZ56sUaHYkIsaib0bdWFaVK\n1BUKxbrD7YIYjs5bLza9WqLuOA6FUo10IkI84om66v2iUCgU7VGaZW2ExdvPhjF0bdU89XLVwrId\nOhIR4l6kruwXhUKhaBM/88Wn3qKHejxqrFqk7leTphIRIqaOpqlIXaFQKNpmdgl+80Sp+zoaMVYt\npTEfEnVN04hHDRWpKxQKRbuUForUbRtD14gYepPXvpKEI3WAWESJukKhULRNqTpb1Jtz0w1DwzC0\nVYvUC7NFPWoq+0WhUCjapewVHnXE3frJcKsAy3IwdJ2Ioa+ap56fJepxFakrFApF+/g54L6Izi4+\nMnQN09RXLftljv0SNajULGxndRe+BiXqCoViHeKnNPoiOrv4yDA0TF2jbtk4qyCss+0XP62xugYW\njBJ1hUKx7vAnShui3lx8ZHqROrAqk6W+qHeEJkphbXLVlagrFIp1h++pp5K+/dJcfGToOqbhyttq\n+OpzJ0rXrqpUibpCoVh3+Nkv6UQUaG63a3nZLw1RX51IPRY1iHhPB3EVqSsUCkX7+MVHHd4qQ3M8\ndV3HNNxWvKuR1jhTqpGKN5bRazT1UqKuUCgUF2S2pz4n+8Vwi4/cz1Ze1POlWjAWaEyUVpT9olAo\nFBemXK0TMfUgIp6dp27qGoYn6rUVFvVa3aJaswN/H9REqUKhUCyKUsUiETUwdT8adyN1x3E8+6UR\nqa+0p14oNT81gLJfFAqFYlGUq3XiUXNOhoufvmgYOqapNX3mv17u3PEg8yUetl9cr1/ZLwqFQtEG\npapFPGYEk6FzRV1rmdL4sbuf5X9+8tFlHUsjR72x5LNvv6zF6kcXXHhaCBEBPg5cAcSA9wPPAp8E\nHOAQ8BtSSlsI8SvAO4A68H4p5d0rM2yFQrFRsW2HStUiETUD39y3X/z/m+E89VD2y7nxGc5PFilV\n6iRiF5S/tpiZVXgEl/5E6VuASSnl7cD/Afw98CHgvd57GvB6IcQW4J3AbcCPAR8QQsRWZtgKhWKj\n4vvUiZiJoXuRuu1H6o1Vj4IoPpTD7n93eqa6bOMpepk4fnMxCEfql6aofx74I++1hhuFHwDu9d77\nOvBK4BbgASllRUo5DRwHblje4SoUio2Ob2nEo8acAiP//032SyhS9787Xags23iKZXefyViLSH0N\nRP2Czx9SygKAECINfAF4L/DXUkr/9pcHuoBOYDr0Vf/9eenpSWKaxkUMu8HAQHpJ31+PbMRzho15\n3uqc51LyhLunK8FAfwqAaNRkYCCNY7h60pGM0tOVACDREQv2Wal5Eb2mL9u1dbynha1bOoN9Gr7A\n6+0fZ7nG05apJITYAXwZ+IiU8jNCiL8MfZwGskDOez37/XnJZIqLG+0sBgbSjI/nl7SP9cZGPGfY\nmOetzrk1wyM594Vtk8uVAMgXKoyP5xnzNKVWsyiXXa87kykyPp6nbtnBpOmZ89OMjy8Yc7bNxJR3\nzHI1GLtfHDWdL7f1O1zs73qhG8AF7RchxGbgm8C7pZQf995+Qghxh/f6x4H7gUeA24UQcSFEF7Af\ndxJVoVAolg2/70s8Zs7NfgkmSrXAb/eLj8L+dnY57RdPwJPxucVHa+GptxOp/yHQA/yREML31t8F\nfFgIEQUOA1+QUlpCiA/jCrwOvEdKWV6JQSsUio2L36ExHio+8r30IKVR14PmWvVA1BvphdOF5Zso\nnfGeCJKhiVJd14hG9EvWU38XrojP5mUttr0TuHMZxqVQKBQt8SP1RLQRqftZL76At+rSWFmpSL1c\nxzQ0omaz8RFfo3VKVfGRQqFYV/iReiJmBHnqcyJ1Q5tjzYQFNruskXqdZDyCpmlN77vrlK5+8ZES\ndYVCsa5Y2FP389TnFh+F/e3pmeWL1EvlWlOOuo+/Tulqo0RdoVCsK5o89Vntdf1I3dRD9otnzYTt\nl1LFWha/23EcN1JvUZ0aixqUq9aqrJEaRom6QqFYV5RDnnpQUdrSfvEjdafpe75Lkl2GaL1as7Fs\npynzxSceMXCc1VmkI4wSdYVCsa4ohdoEaJqbulifPVEaXvnIao7U+7viAGTzSxd1P/NlPvsFVn+d\nUiXqCoViXeEX9vil+KahNyL1UJuAOSmNnrhu7k0Cy9P/xc9RT7QQ9YXWKf3Xbx3lOwfPLvn4rVCi\nrlAo1hVlT0hjgahrC3vqvqh7XvyWHlfUlyNS9/u+LBSpzxb1qVyZ7xw8y7OnppZ8/FYoUVcoFOuK\nUtUiHjXQPXPcCEfqfpdGQ5+bpz4rUs8uQ6QeFB7F5nrq89kv8ozbPeXqHd1LPn4rlKgrFIp1Rbna\n3AvdNLQ5bQKaWu/WmytKt/iivgwFSAtF6vPZL3IoA8C+nT1LPn4rlKgrFIp1RaliBX46uAti+LZL\no01AC/vFE9fNPW73xuVoFRC03W1pv5hNx/U5MpQlETPZsSm15OO3Qom6QqFYV/jrk/oYoUi90SZA\nnyPqfsTckYiQSkSWJVJv9H1pkdIYrH7UqCrN5CuMZUpcvb0LXdfmfGc5UKKuUCjWDbW6Td1ySMRC\nkXqTpx6eKPVTGv08dVfUYxGD7lR0WVoFLDhR2sJ+8a0XsULWCyhRVyjWLZZtk1vGZdnWA+HCI59W\n2S+GoaFpWtNn5ZpFNKKj6xpdqRilSn3JZfyNtrvt5akfGXInSfftWplJUlCirlCsWz7/3RP83kcf\nDPK2NwJ+4VE8FKk3Zb+Eio/8z8L91H3bpjsVBZa+rF2rpex8Wk2UyqEMiZjBzk0rt6KVEnWFYp0y\nlStTrdsUvNXsNwK+CKeT0eA9U9ewHQfbdpomSgEiIcGvVOuB0HanYsDSuzXOlGtoWvNNxieI1D1R\nz+QrjGZK7N3evWJ+OihRVyjWLb5YrXZvkbVkIuuuuzPglfoDjaZett208LT7mdaU/eJPXjZEfemR\nejJmBjnzYfxj5YrujUOeWdlURp+21ihVKBSXHr5YbSRRH8+6a5IOdCeC98JFRn7xkf+eO4lq4zgO\nlaoVRM9dHW6kv9RIvVip09Ei8wWgJx2jI27y0DOj7ji835PYuXJ+OqhIXaFYtwSibm0gUZ+eK+pG\nqKd6uPgICMS0WrdxoOGpp91Ifame+ky51rLvC96x/uAtB9i5OcX3nzrPQ8+Oun765pXJT/dRoq5Q\nrFMudfulUKpxfnJmWfc5ni2jAb2dc+0XN1JvIeqW00hn9O2XIFK/eFGvWzbVmt0yndFna38H7/35\n5/PqW3ehAddf2RdM4q4Uyn5RKNYptUvcfvn8d4/z8OFRPvQbL2mZ8ncxjGdL9HTGgg6M4E6Ugpv5\nEu79Ag1P3V9Wzp8o7VqGidKZoJq0tf0SjM/Q+ak7ruIVN29ram+wUrR9BCHEC4G/kFLeIYR4HnA3\ncMz7+KNSys8KIX4FeAdQB94vpbx72UesUCiARvrepSrqk7ky1ZrNeLbEri1LT+Gr1W2y+Qp7ZzXC\nCtYptZ3GRKkfqZtuSqMfqfuTlxFTpyNuLilSLwbNvNqT0fDTxUrS1miEEL8HvBXwn6UOAB+SUn4w\ntM0W4J3A84E48H0hxLeklMu3GKBCoQjwKyVr1uqvg9kOJa/V7cT08oj6ZK6MAwx0N4ujXznqRupe\nRakn9BFDx3EaPdhjoZ4x3ekYmdxSRH3+atK1pN3RnADeCHza+/kAIIQQr8eN1n8LuAV4wBPxihDi\nOHAD8Oh8O+3pSWKac/M7F8PAwMol8V+qbMRzho153guds7/yZSIRvSSvjW8PlerOosY337ZnJt1J\n0iu2dTdtk065Ip/uTAR6smkgTXc6RiLhWiN6xJW6vp5k8N2B7iTnxmfo6k4SjSxeh05PFN399HUs\ny/Vfrt9hW6IupfyiEOKK0FuPAB+TUh4UQrwH+BPgh8B0aJs80LXQfjOZ4uJGO4uBgTTj4/kl7WO9\nsRHPGTbmeV/onH2feCpTvCSvTcHLzz59brrt8S10zsdPTwKQNPWmbapV1wYZnyxQLLnHzGZnqJWr\nON6NZXg0B4BVs4Lvat5tcXhket60xIU47+0Ty17y9V/s3/dCN4CLnYb9spTyoP8aeB6QA8JHSgPZ\ni9y/QqG4ANYlnv3i+9gTXhriUhn3C49C6Yzgtt4F93q0yn4BgqrbcMtef7K13et3/Nw0v/sPD3B2\nvAAs3PdlLblYUb9HCHGL9/pHgIO40fvtQoi4EKIL2A8cWoYxKhSKFtQu4Tx122mkEU7kysuyTz9H\nvX8eT71u2U0LT4c/myl5nnrk4kVdDmXI5CsclOPuPhfopb6WXOxo/jvwd0KIGjAC/KqUMieE+DBw\nP+7N4j1SyuX5bSoUijlcytkv4SZWE9kyjuOgtSilXwzj2RIRUw+qQX2McJ76nDYB7md5z5ZpjtTd\n1+1eP3+h6hPnXJfZz365GOtmJWlb1KWUp4BbvdePA7e12OZO4M7lGpxCoXAbQVVqVrAMG4DjOJd0\n8VG4c2SlZlEo1ZqacF0ME9ky/V3xOTeH5jx1B13Tgl4spheN+5F6eHGNiLG4SD1fdEX8xHAO23EW\nXPVoLVEVpQrFJc6ddz3DX3zm8ab3fO8YLlFRn7WE28T00h7aZ8o1ipX6HD8dmvPULdsOonRo+O2+\npx5rw1Ov1W0msnPnAfze9aVKnfOTxVDbXSXqCoViEUzmyhSKze11w0J0KXrq/mIWUU84lyrqE/NM\nkkKzp25ZTjBJCmCa7utWE6XRQNSbb0DffHSIP/inh+YIu99tEVwLprGUnRJ1hUKxCGZKdSzb7Rfu\nc6lH6mWv8GjbQAew9AyYoDtj19yqzPBapHW7WdQjRhuR+qyb4sR0Gct2OD/VnHKdn6kG+z5xbppi\nuU48aqx4L5fFcmmNRqFQNGHbTpA6Fxafpkh9hUTdth3+16cPctcDJxf9Xd9T3z7gdiRcaqTequWu\njxFUlDpYlh2IvPuZ56l7UXU8lP1izmO/VGvuz9l8o9rUth3ypRpXDKaJRQxODOe8truXVpQOStQV\nikuaYmjCsR4SdWsegV9OpmeqHD83zX1PDi/6u6XqLFHPLlHUvZtCfyv7RW9E6pbtNHnqfqTuOKBr\nWlMjsPk9dfcpIxta/7VQquE47uIauwfTDE/MkC1USbRYxm6tUaKuUFzC+BEmzO+jr5Sn7h97Mldp\nOXG4EH6Oem+nu1DEctkv/QvaL27xUZOnHhL4WNRoypyZL/vF/znc7Mv30zuTUa7a1uUdb+G2u2uF\nEnWF4hLGT8UDgpVzoFFNCisXqc+E1j6VZxZXHF72njDiMZP+rgST026u+sUykS2RTkZatq4NGnrZ\nNpZlN3ncYSsmPEkKjUi9Ott+qc+1X/zMl3QywlVbG91PLrVJUlCirlBcMnzpvhN85+DZpveaIvV5\novPlEPXnhnNBxopPYQmiXgq1uu3vilOt20Ge92KxbYeJ6TL9XXOtF2guPqpbzfbLQqIenaf4qBGp\nN+wXP1Lv6ohy5bbO4P2LLTy67+wPeGr8mYv67oVQoq5QXAI4jsPXfjDEtxcS9RWK1E+ez/H+Tz3G\nNx4eano/LOpHhy4uUk9ETfo8y2T8Ii2YTL6CZTtzWu76NKU0zrZfzMbr+SL12fZV1ffUQ/ZLfsa9\nFulklM5klE097g2mVaT+zOQRCrX5V3zKlLN89uiX+cH5x+bdZikoUVcoLgGKlTq24wSdDX2a7Jew\nkIeEqL5ET/3pE273w/FZk5l+b5NoRGcsWyKTb7/3uB+pJ2JmkLEyeZEZMP6SeL6QzsaPxi1v4elw\ndG6GrJhYZB5Rn2O/WGBWyc1UsT3LKPDUvRYFvgUzW9SPZ0/ykSc/zrdP3zvv+Twx9hQA1/SJebdZ\nCkrUFYpLAL+4aKZcbxLp5ki9USRzMdkvZ8cKvOvD9yOHMk3vP3tqCmj0RwnG5EXq11/ZBzDnewsR\neOpRI4jULzat8bB33Ku3d7f83I/MWxcfhe2XZgH2Rb0+6/oVO48Qv+l72GYx+L34nrov6nu2u6Ke\nTjTbL09NuJZKpjL/k83jY0+hoXHTwHXzbrMUlKgrFJcAYb85bHu0E6m3m/1yciRHvljjvx4/F7xX\nrtY5MZybM4bwOA5cPQAszlcvVS20ZI7HJw4GGSuLzaDxOXI6g6Fr7J1H1P3IvFq3caBl8RG0sF/m\nyX6pd5xH0220ZD6wYPxr05l0Rfy267bwU3dcxQuv2dL03WcmjrjbVwstxzpZynAyN8TVPVeRjqbm\nP+kloERdobgECEfJYXEtzjNRejGeul9U8+SJCSo1N+o/eiYbVKfOtX7cY19zRS+xqMHRRYh6uVIn\ntv0En5FfJBJ393sxkXqxXOPUSJ4rt3Y2VYOG8T11/5yMsP0yK6UxTMNTbzwBVawqTsLtwqjFioGo\nT89UMQ2sOQtbAAAgAElEQVQtyL6JRgxefeuuJvtlojTJSHEMmF/Unxh3rZebN92w4HkvBSXqCsUl\nQLi3Sz4krr6vDc02wcVkv1Q90avW7MBHf/aUa20Yuka+WGtKO5wp1dA0SCUj7N3WxfnJYtB+9kKU\nqnX0qDv2vDVNKhG5KFGXZ7I4Duzf1TPvNn6k7rf7bc5TX1xK46npIdDca6DHSkEGTL5YJZ2MLtg+\n+JAXpQPka/OI+tjT6JrOjStkvYASdYXikiAfslzCjaMK80Tq9Vmv28kB9yNZgMekG1E+eypDxNTZ\nu72Lat1u2qZQrtMRj6BrGmKna320G62Xqxaa4RUvlabo74ozMV0OJh7b5fBp96azkKj7kbk/9qaJ\nUmNxE6XHss8Fr7VYiWy+guM45GaqdF6gdfChycMA9MV7mKkVsZ3mm+1kaYpTuSGu7l456wWUqCsU\nlwTNkXrYfmlE6mHxCfvr7s8XjtZ9+8XQNZ48PslEtsTZ8QJXb++iJx2fc+xCqRZUTIodrqi2O1la\nqlhguGOfKLuiXrdspgvtRfo+R05niJo6V26df7lj32Lxn0QWO1Eavq7HM16fG1sP7JdKzaJat4NJ\n0lZUrCrHss+xLTXI9tRWbMemWG+eQ3hi/GkAbt68ctYLKFFXrBBnxwocO6uWqG2XsOUSFtZwVefs\n6DxMOxZMxcueuf7KPio1i8999zjgeuZpbwLQP7bjOMyUaqS87I4rBtNETJ3j56Zb7LmZWt1dVs7R\n3XOaLE0FPVsmF7G0XW6mytnxGfZu72rq2TIbP22x4am3bhMwb566d+0s2+JUfgi7mCJipdFiJTKF\nCrlZk6StkFPHqNt1ru3bR8qLwguzfPXHR59acesFlKgrVoiP3f0sf/fFp9d6GOuGsP3iC7zjOPMW\nH/mi7lsK7Yh61fOcb7vezdh4zFtrc/8VPSFRd49drlpYtkOHJ+qmobO1v4PhiSKWvfCx3MpUG0f3\nInXPfoHFZcAc8Z4K9i1gvUBDxMstPPVw9svsiVJD19E1LbC1zhaGqdk17HwPCa0TzbDIFPNz0hlb\ncWjS9dOv69sfWCvhydJSvcTp/Bn2dO0mFelY8HyWihJ1xbJjOw7Dk0VmSrUl9fvYSORb2C+VmkXd\ncoKl2Zojdfe6xmPti3rF22bP9u6gOrMjbrJzUzpYas4/tv+EkArlYW8f6KBu2YxlFhbmUtUCs2Eb\nTZanghL/8UVMljb89N4Ft/NFvFKvgNZ+7xdwo3X/2h3PutaLXeghQRqAbDVDPuj70lrUHcfhmckj\ndJhJdnftJB3xRD1UVZopu084mzoGFjyX5aBtURdCvFAI8T3v9R4hxPeFEPcLIT4qhNC9939FCPGY\nEOIhIcRrV2jMikucqVzZffxmbrOkjYhtO3zqHsnxs/NbF4VSla6OKJrWiJZ9n72zwxXWpkjde+0v\npdZOrrrvOcciOs/ftwlwJyB1XWtE6l5qZaHFosp+G92z4/OXwIObzqgZDVGfruToTrvjnFxEq4DD\npzMkYga7tiw8qahpGqYBtb3fJbL76SbLpcl+iSws6id8Uc/3kNRcD3/GyQUteP3fw2yG8mfJVqa5\npk+gazrpqBuJhyP1TMX93ffE5p8bWC7aEnUhxO8BHwP85gsfAt4rpbwd0IDXCyG2AO/EXZD6x4AP\nCCFiyz9kxaXO6FTjH25l1lqVG5Gz4wW+98Q5vvHI0LzbFEo1OjuipBKRIFr2xb0r5f4zasp+8SwQ\nf/KvLfvFE/VoxOD2G7bS2xnjtusHAeZE6oUgUm9MLvqifmasdbqeT7lqgdF48nBw0GPu38TsVgTz\nMTldZixTQuzoaWtlISNWR4uWMXrGCG/eFKm36PDoirqF4zicmD5F2uzEqSboMNymXVq0yNlx93zn\ns1++M3QfAM/ffBNAS/sl61WYdl8qog6cAN4Y+vkA4Dc3+DrwSuAW4AEpZUVKOQ0cB1Z2mldxSTKa\naSwDFk6RWwumcuWmCs21wL8GJ8/nWn5eq9uUKhapRIR0MjonUu/2xKReb1hZ/uvkYuyXmkXEdH3k\nLb1J/vrXb+PGPf1AYxLQP7ZfyTrbfgE4N76wqJcqdTTPfono7ven69N0dkTb7v9yasS9Vnt3tCeC\nRtQdt2ZYlPVGhs5CKY3QiNRHi+MUajMMxrcDkDbd42qxEmdGPVFPRqnbzZ0sR4vjPD72FNtSg1zb\ntw+AlGe/hJt6ZT37ZTVEva1mwFLKLwohrgi9pUkp/b+wPNAFdALh50v//Xnp6Ulimq2rxNplYCC9\npO+vRy71c86FStuTqfiyjXex+7Fth9/5+wcQu3p479tfuCxjuBjOeE8umXwFPWrSN6uFrG9JDHj/\nHoYnZujp7eD4SB6ALQMpnjwxiRkxgmsQ8fzhrk7PG2/jOluOG9m32q7DS2ms1B0GBtJohjuJOri5\nM9i+vz9FZ0eU81PFBY8VOTMdROpXdG/j2NQpKsYMg/0dnDibpbcv1TSZOZuBgTTW0Qn3+9u62/q9\nm/Ea/q27HJto+o6ua9i2w7bBrjkrJyViJqVKnVHrPAC7unbzJBZbujbBpCvqZ4ddUZ+OneOv7/sX\nfvXAm3n5lS8G4POPfBkHhzfd8Bo2bXKj+2jalcaqVg7GUT7pBjpXDW5loLP1+SzXv5OL7fAeDgvS\nQBbIea9nvz8vmUxxoY8vyMBAmvHx/JL2sd5YD+d8arhxbx8ZzdFhzv8PuF0u5ryzhQrZQoWhkdya\nXrOx0LEfO3Sem69uniwbGnU/j+gasYgbWZ4amgqskJjnC+cLleA88l75uu6FVuOTBca7FnY7i6Ua\nEVNreS0cx8E0NCazJcbH84x40bhVqzdtv62/g8OnM5w5l5mT9x2c72Qh8NS3JAY5xilOTQzT3bGb\nuuVw/OQEvZ2t2+j6v+fhUe+ppm619bvTzMYTwET1XNN3TEOjajvM5Es4teZIW8N9gjk4dAiAVG0A\nGMGpasT1JMVYMciqefTs41i2xf9+7F/Rq1E2Jzdx36mH2ZzcxJWxPcExbcdBQ2Minw3eG5l2K3it\nGZPxytzzWezf90I3gIvNfnlCCHGH9/rHgfuBR4DbhRBxIUQXsB84dJH7V6xjRkOrsJfX0H7xW8Ve\n7OIMy0XYgmplwQT+dTISTFjmirWgF4vvqbfKfvF7kbTrqbeyIMCdbGyyfvwxzVoEYptvwUzMP1la\nrlhopvv9HamtgNvIyk9rHG8jrbGdNMKm8Ucaef5Z53xT1pWf1tiqd0zE1KnZVZ6aeJZNiX6S9ATv\nd0W60aJlwCEZNziaPU7CjKNrOh879Gn+/eiXsB2bH9v1cnStIaW6ptMRSVKohT31aRJmgri58tOM\nFyvqvwv8TyHED4Ao8AUp5QjwYVyB/y/gPVLKpa02q1h31C27aTKsuoYTpVNeoctMuYZtr0xqZbla\n5+sPnQ4mIVsRnixuJer+TSedjAStXPPFaiCs3SlX2FrlqScW5anbwWo/rUgnIkG+vJ8fn5rVWjbI\ngFlgsrRUqQfVpJuSA0T0CJOlyUauehu+elDw06ao44m6XYlTochkeSr4yDB0IqbecsI1YuroPSPU\n7Bq3bDkQ3Cwjpk5/og9Nd9CiZVI9ZfLVAtf1XcPbrvlZqlaNZyclffHeYII0TDqampP9shqZL7AI\n+0VKeQq41Xt9FHhZi23uBO5crsEp1h+TXn8PTXNXcF/LSH0q50bqjuMuQjFboJaDRw6P8fnvnSCd\njPKSGwZbblOpNQT35Pk8ttPIPYdwpkljfPliLRD7I+VHiVwlqVovDz5viHp7kbrjOF6kPn8cl05G\nGBorUK1ZFFpMlELrtMZ/uusZKlWL/+sn3byIUrUeROrJSIK+RC8T5Qz9A66f3Zaoz1TRNS0ofrog\nnv1iTW1BHzzFiewp+hNuH/iIoWHP84QSNQ2M/mEAbtnyPH7oPWVGIwZbYn08k3W7NRqd7vXY17uH\nmzfdwPTeHF86fjevvfJVGPrcfaciHZyfGcWyLWp2nVK9xBWdO9o7lyVy6a2aqljXjHj/KAb7Ohie\nmFnT7JfwSj35YnVFRN0v0plaoPzdX/uzOxUlW6gyOlVksK9RVehbHuEugLliNXj/0PQPMfuyZKdP\nAG5UOMd+uUCees3rNR6dR9z847vjqVEo1TANjeism8C2/g40Ghkwp0fyPPTMKJp3jIipU640UhqT\nZoL+eA8jM6N0pNwxt1NVmpupkk5Gmm5+C+GY7u/azmyGwVMcz57khYMHANi3s2fe62NHihjpKXan\nr6Av0Uu17s1vmDqbvZuCFitRT7oZNaJnDwAv3/ESXjT4gnntFD+tsVCboVx3/zZWK1JXFaWKZWXU\nqzbcudn9o17LPPWpfENoV8pXL3or/GQK8y/15jfS8jsNzrZgfMsjnQjbL56nrtlMV92J56nUU1Qs\nV+gXG6n7RWDzeergevrueKrMlGp0JCJzWs3GogYD3QnOjs/gOA73PHoavWscOrJMeFk85WojpTFh\nxunzxNE2i2i0F6lPF6strZdcNc8HD36Eg6NPNr1vG2WcegS70IVJhBPTp4LPfum11/Brr2/db6UQ\nc7e7ofdGAGre7ypq6gwkPVGPz1A0R9mc3ERPvLFQx0L+eDhX3S88Wo10RlCirlhm/Bz1nZvc2fm1\njNSnQpH6SuWqlyvu+S20fqdvQe3b6Yt6c5ZDIeype0JW8Dz1WLKGg4Nj69hGiW+f/h7giroGgZ0S\nXuquFf7NdXbkHSYcqc+U3WZeM7UitVm52dsGOihUZ/j8M9/kh+bniYmDxPY+HtzQS17xkYZGzIjR\nH3fPO1vL0J2OBeI/71hrFpWq1VLUv3TsP3lu+lSwzqePbZRxalFApz+yldHi2LwLVfg4jkPGPIFj\n61zduR9o3Pwipk5f3G1PYPSdx9bq7Ovdu+D+wjRaBRTI+qIeV6KuWIf4mS87LoFIPROyRFZK1P1I\nPbuAqPvXYM/2LgxdmxupezZLRyLS1C0xX6wRS7n71SZ3odXjfGvoXrKVaeqWg2nqRMz2Jkqrdb9F\nwEL2i1coVKhSLNdJJhz+9Ad/wT8f+nTTdoP9cWLXPsS9Y98Bs4puR9GiVc5PuTersld8lDATaJoW\nROp+X/WpfGXBVsF+r5XZXRHl1HEeHX3c3VdoItSyLWy9ilNzI+dN0W0APBeK1lsxlD9LWZ/GymzC\ncNwbiH+doqbhRtaOhh5z/472edZLO6RCkXog6rHWy/EtN8pTVywro1MlulNRuryob60iddt2yOTD\n7WwX18e7XUpt2C/+NUgno2wfSDE0mqduNVa9L5RqJGMmpqGTikfQcD31QrFKcmuZGmBUu4hM9lDe\n/DhfPfENLGs3pqGFlmRbWNQroRYB85FOuL+z0UwRB7A6z1Osl3h64jDPTEqu7RMAZJNH0K0i9fFt\nRMeuR9x6jsP5pzmXnQR2U6pa6EadpJkEoD/hRryT5Qz9XZs4dnaaqXyFTbMKgXymi3PTGWt2nc8e\n/TIaGnEzzmS5UTXqrzLkRuowGNvOUzPw4PAjHJk6xtHsc8zUZugwkyQjSUzdpGbVgsWhrYmtgZjX\nQ5G6oRvESVEmj4bG3p4rF7zGTdfS68RYCNkvylNXrDtqdYupXJnNPckgJ3itIvXpmSq24wRpdCvl\nqfuToPlibd5ouRJqpLV7ayd1y2nqn5Iv1gI/W9c1UskI0zNutGzEXavCrKfQMtvZlhrkkZHHKWs5\nTEOfd/Hk2fi+/kKi7jesGpl0n7Zm4qcB0ND40vG7sWyL6UqOQzMP49Si1E7v5+U37mJzp2uvjOUz\njWti1khGXNHu8+yXiVBa46Q3Wfrk8Ql+9x8eYCxUiJifmZvO+O3T9zJaHOel21/M7s6dzNSKwQRk\nzi/m8SL1zbGtGJrBockj3HfuB0yVMySMOPlagZPTpzmaOc7p/BnK9Qrd2iB2rj8Qc99+iXo3y37P\nV9+R2k7CbH0TakUQqddmVrVFAKhIXbGMjGVKOMDm3mTwmL9Wkbo/Sbprc5qJ6ZXr/+LbLwDThcqc\nMnRwb2y6pmEaOrsH03zvCTh1PsfuwU4cx6FQqgViB25E7wsrUVf8ok6KsgWv2vVyPvHMZyh3nsAs\nXhdE6vl6hm8P3csrdtzeVAgTjCF0Y5kP31M/P1UEs0JOP88VnTvZ2rGFB88/woPnH+F07qzbc3z4\nOgwivOLm7TyecVMCJ0uueJUqNdCtQATjZpxUpIPJ8hQ3ei1/x6fL7HMcvnL/STL5CnIoy6YeN7L3\nl/PrTEZxHIdHR5/gntPfoTOa5ieufBVfOfF193jlDNtSg+Sqrqj79kssEuXn97+JifIUV/dcxa70\njiDt0Hbcpf/8n//zB6f4ovNccFMMPHXv73dn9ybOlk5xbb+Y97q1vJahhTKylWmiRpSE2bqKdrlR\noq5YNka8HiebexONSL124aKYlSDj5ajv3Jzi4NHxFZ8oBdeCaSnqNYtYVEfTNHYPuv1BTgznePnN\nrn1j2U5Tr+50IsKwVxFpmQUMzSDidFCoV3jewPV8KdrJdOdpopP7PFG3OWr+F08dz7AzvZ2re66a\nM4ZqO/aL97QwOlXE6B8FHA5svpEDm27i4NgP+Y8T36BcL7O1YwsvuunlRHSD7lSMrhl3UjxfK1C3\nbMpWmTg0RbZ98V7OFYZxerPo3WMcyWgMnk1y2muRMB6aPPUXt7ajeT78xFc5mj1BRI/w5n0/ScJM\n0O9NYE6UpjxRb7ZfTF3jpi3Pa3mOuqa7vQE8Zj/p1Lzr5L+/u2sXD408xo0D18573Vpey1kTpT2x\nrgUXrV5OlKgrlg3/EXpzT5Koqbt9Nar1hb+0QviZL4N9HZiGvuIpjTB/BkwlVJ6/tb+DjrjJkaEM\njuME6YzhHPp0aIKwqhfojXdjmAY1y8bQDV66/UXc9dw92N1niJi3Ymw6G3QmPJZ9rknUn5mUfOzp\nT/GKrp8CFp4oTcZMDF3Dsh2ifW6Dq5s33UBXLM2rdr2Cu577BgBv3Pta9vc2Cmk6o14fErPM8MRM\nKEe9EZn2J3o5nT/DZ89+ktjV8KQNuSctfAc4XIWcn6miJXN8fvgTWI7FdX37edPVr6fP8+b7Ao/e\nnSydHakbRvuu8uw5icB+8Z5obh08wA0D1yx6taKEGcfQDKbKWQq1GbalWhemrQTKU1csG3464+be\nJJqmEY0aaxap+8VAPZ0x0snIkidKHceZ02rAcZxgohTmz4CpVBuirmsa+3b1MJWrMJYpNbUI8PHT\nGtHr1CjTF+8lYmjU6651cNvWF4KtU+l6jopTIrL9GLoTQUPjeOa5pmM/fP4xqnaN48VngYVTGjVN\nI5WIoEVLGOkMg7EdgQ/8ih23sy01yPM338T+3qubvtcZ83uPVzg9kg9WPQpH6nfseAkHNt3Iy7a9\nGGt0JwDHc0fZ2t+BoWtN/WByxSrmpiEsx+LN+36SX7vhFwIhh4ZHP1Vyb2RzRH2BDpCzMWetU1qz\nbDStsQ9d0y9q+TlN00hHU5yfGQVWz08HJeqKZWRkqoQGbPJ801jEWLM2AX7U3JuOk05Elmy/3Pfk\nML/+oXubovG6ZWPZTlAANF8GjGu/NCLka7wipMOnM0GOeios6l7UrnkLS/QlejFN3c1IsR3S0RT2\n1FbsSIF/fvaTaGaNnpnr2Jrawsnc6SCv3HZsjmSOATBcPQk4xObp/VKul7Edm3QygtE7AsD+robl\nEDUi/OEtv80vXvvmOd/t8iJ1LVJhaLSA5kXqYVG/smsXb7/u53iTeAPp6RtxLAOta4xXHthOf1e8\nSdSzM0WMvvN0x7p40eAL5tgWvsBPzInU3ZtheOHpCzF78ema1x9nOaySVKQD23H3u1qZL6BEXbFM\nOI7D8MQMfV3xIHc6FtEXbHS1kkzlyxi6RldHlFQyQrlqtdX0aj7kmSzVuu3aCx5Fz0/f2udO8LWy\nXxzHaYrUAfZf4YrSs6czwfJxfjohNCYsfVHvT/TO9X5H3Gj33MwwdqmDjsIe9nRfSc2uczp3BoAz\n+XPM1Nynp6KTQ0sUiLboVDhenOTd338fv3//+ygPPoqx6QyOo3FD//VtXZu4GSeiRdAiVdcj95p5\nJSLztNft7MCe7kNPFNlzpclAd4J8sRY89Uxqp9EMi1u23Nxy0rfDTBIzokyWXFHPV/3sF99Tb1/W\n/AZn1WCi1AqEfqn4k6WweoVHoERdsUwMTxYplGrs2d74441FzKAX9WozlavQnYq6KYJe5BuO1o+f\nnebwqan5vj6HiWyj46NP2ROhLb1JNFrbL37PlXCkvrknQU86xpHTmcB+CUfqfiqfFnMFuS/eG9gE\ndcvGth3sYiexqrvOqHVmP/W6xp7u3e65ZV0L5tnJowBc338NAEb3eMtI/dDkYep2HcuxmIkPoceL\n2NN9bEp1tn19OmNpN1Ifyzeaec2TAtjflcCadnvKH80dZaC7udFXMXkKgFu3HGj5fU3T6Iv3Mlme\nwnEcctU8EeL4crYY+6URqVve/+1lE3V/BSRQ9otiHXJ0yGt4tKNRNReLupF6uLf1amDZNtOFKj3e\nQgx+FBz21T/2n8/yD18+tODYTufO8J8nv0XNrgel7eEbgz9J2pGI0NkRbWm/lGtzKzk1TeOaXT0U\nSrXgxpJORjiaOcFQ/mxgv+jhSD1kE/iTelsKL+bXb3w7ZnEztbrN3m63OOa4t4Dys1MSDY037nmt\nWxnZPdbSU5eZ4wD84S2/w431N1I9dQ21U9e03yER6I51okUqVGtWsEDGfCl8N+7pp19znzQOTRwO\nRH08W2KsMIWTGida7WNzx6Z5j9eX6KViVZmpF8lVCsS0ZPDZouyXWU9A1bod5KgvFX8Bali9alJQ\noq5YJuQZtzrv6pCoxyMGlu0EHQVXi+mCW3jUm3YnzvxJSF+Qa3Wb8UyJYqVOttB6AnW6kucjT36c\nr538Fp85/EWynmDPhETdtwsSMZPudIxMvjrnJuH3k5+9kv3+K3xf3b1uHXGdjz71Cf7pqU8Fiz0H\nkXqiN6g+rVs2lt/MS0txbZ+b1lizbNLRFFuSmzgxfYpCdYZTuSGu6NzJpmQ/KWcAPZXF0pqbaVm2\nxbHMc/Qn+uhL9LAluRlrbCcx0k3re16IzlinmyoYqYA511MPc0AM8IG3v5wd6W0cz56kq9MV4fFs\niW8fexBNg35r4ZJ8P61xZGbMTaHUGgK6lOwXN1Jf2hKbPmH7RXnqinWF4zjIM1k6O6Js6W1ETNE1\nKkAKT5JCw9rwRX0s6xZJAYxMzl3Bx3ZsPn34sxRqM6QiHTwyehBj82lvH41sl5LnqSdiJj2pGHXL\nZqbcnMLpR+qzvez9u3q9Y7kjKZGjalXJVLJkbDdjQouViBkxOsxkU0RZ826Svuj6iycD7Om5kqpV\n5dtD92I7Nvv73EyVVG07mgYnZ040jWMof46yVQ5ayvo3wMW2KQ5PlvqR+nz2i891ffuwHIuC6aZP\njmWLPHjmURxbZ3v06gW/60+Wnpx2fy9xvfF3Z16U/eKLurVghtBi8O0XUzfpiCQvsPXyoURdsWTG\nsiWmC1XEju6mrIH4GrUK8HPUezrdSD0VamcLMBZabu/81Nx1cr939gEOTx3lml7B77/gXSSNDiI7\nj6B3jTd56o1I3aDHeyrwfXXHcTidO0Ox6j4JzI7Ue9Kx4AZo6BoT1bHgM5l7Fg0HLVaiP96LpjV6\nvNQtJ4jUTc9miBgNUd/btTs4B4BrvPTDeMldVu5Y7mjTOHzrRXi57f4kbUd8saLupTVGqkGe+oUq\nKK/rdzsjnqu4mTkni5Kx0hhWZhO9HQsvwuynNZ7MDbnHMkKR+kWIer3uzlXULWfZ7ZfuVSw8AiXq\nimVADs21XoA1axXg56gH9susiVK/RSzA+YlmUT+TH+Y/jn+NdCTFW695Ez3xbm7teC04OtGrnmSs\ndjbYtuQVViWirv0CblqjZVt8+vDn+MvH/o77R78HtC768S2YVDLC8MxI8P4T40+TTDtohhU0w/Kj\n8prV8NTni9QBanaNpJlgl7/aTjWFXU4gs8eoh1rp+qJ+9ZxIfXF1ielYKFL38tSTF4hOd6a3k46k\nOJw5QmL/44x2PoCGhjW6i64LLGPnR+qnvEg9aSzNfqnW7eAamsuc/bKa1gsoUVcsA0c9P13snCXq\n0dUXdduxOZj/HnrPSLBifSrZPFE6mglH6g37xbIt/r/Dn6PuWLxl/083KiVnuqk9dz3oFuc6v8N3\nhu5rKjxKxF37BWBsOsdHnvw4D48cBOBI7mnAaVme7+erpxNRzhbc/inX9+8nW5lm77Vu6bsvXmGb\noN7KfvGEvjvWxYDX6nZf794gJbBas3GmN1GxKoGQ16waz02fYltqMBCgIFK/WPslWgGjho5OVF94\nH7qmc23fPgq1GUiPY+f6+LGet2IXepqKsVrhR+rTXjpj8qIj9Ubr4nDb3eWgJ9aNhsZAon9Z9tcu\nqk2AYsnIoSwdcZOt/c2Vd0Gkvor2ywPDD3NOO0R0t0kk9ipg7kSp3/M9nYxwfrIh8Peee5CzhWFu\n3fL8wBoAdwLPmhrEqcWJ7/0hXzp+N6dzZ0hU3JS7RNSENGDU+Fb28+SdCa7r208ykuCRkcfROyeJ\nR/fNGeu+XT2YhkZ/V5xzebfY5iVbb+XpicOMRZ6FCsFCDb7VUrfsoKNgIOqe/eI4Dpqmsbf7SsZL\nk+zvbTShqtZsjJntsPk0Xzx2F3u7r+Tk9BB1ux746QBb+pK86gU7eN7exQlRl19V6nnqcTPeluXw\nip23U6jNkDu7GXkkxqRn+1xowem4Gacjkgzy8DvMFODeCM2LzH7xI/Xlsl+6Yp381s2/xubkwLLs\nr12WJOpCiMcBv+P/SeDPgE8CDnAI+A0p5drUiStWhYnpEpO5Ms/b2z9nPUk/Ul+tqtJCbYa7TtwD\nDmhmne+MfJu39f7sHE99NFOit9P1tJ89laFUqVNhhrufu4ekmeANe17dtN/JabeQqT+yjcKxTnbf\nem0UexAAACAASURBVIyDY0+yq9YBdJKIGcSiBuamIfLOBC8afAH/TbyRk7khHhl5HKN/uKX90hGP\n8AdvOYAZrfHnT+a4tm8f+3r3kjDjQa/vvoQbkYbFp27P8tRDfnvE1Nhl3MjB3AR7Uo0bSaVmEa33\ncduOl/DdM9/nKye+Ttxwny7CvWJ0TeNnf6T9FX58/KcaM1bBMWskzNQFvuGyLTXIf7/xF/li9gSS\n0xw57aZ4diUXFnVwb3i+qLv9y11RNxZRfBTOfvFFfbny1IGgdmA1uejRCyHigCalvMP77xeBDwHv\nlVLejpvg9PplGueG5eT5HGfHFl6Wa7mYypX5vY8+yEE5duGNPQLrZcfcPFxfyFaiqrRUqXP4dKbp\nKeCuE99gpl7EHLsGvdzFIyOPcyJ7CtPQScQMCqUalZpFJl9hc0+SwV73yWI0U+QLx+6iYlV5w55X\nN6WigVsU0+v1kCkWDH7x2jejazqjugTcRlhdHRGMgbNojsFP7v0JDN3gqq4rSGppjJ5RdLP1Ndg9\n2ElRc3P8t6UGMXWTG/ob5fn9npUSCRUfzYnUZ61+dPaMRvaIYHi0kTdfrVnETIPXXfnjbElu4t6z\nD/Dg+UfQNZ093e0v/jAfHZEkuqZjJupg1OiItN97HAhy1ce8+Y70BSJ1oKkfTEfEs39we9K3i39j\ndO0XP1JfHvtlrVjKLelGICmE+KYQ4r+EELcCB4B7vc+/DrxyqQPc6PztF57iH7/6zKoc65lTU0xM\nl7n7wdNtf6fhp/fM+cwX9fmqSnMzVZ44Nt72sRzH4f4nh/ngvz/BO//2fv7q357gaw+5Yx3Kn+WB\n4UfYnNxE4cx2+vLPB+BzR7+C7dikEm5TL180Nvck6Os10BJ57j3zEE+MPcXuzl28aPAFTces1iym\nZ6r0dyVIxSM4DsTo4Lq+/ZSNSbTkNPGYyVDxJHq8RKywI8j60DSNQf1qNMNiuNbcaCvMOc9P9zv5\n3bzphuAz3zs2myJ131N3BcmclWvtP5GEi60qNZtoxCBqRHjbtT+LrunkqwV2pXcsS59vXdPpjKaJ\nJStohn3BdMbZDIT6yceixoLdJH38a+M23fIyiRZhvQBBZlHYU1/OSH0tWIr9UgT+GvgYsBdXxDUp\npZ8CnAcWnPbt6UliLvGuODCwcOrTeqZcqZObcZc16+pOBpNtK3XOk97yb6dH80xXLPZsb10F99jh\nUQ6dmGBkssgTR8dJxk2ed+3gnAmqgX53EisSM+eMuVSp86effJShkTx/81svY48X6duOzQPPnOIf\nvvg4f/ZLr+Sq0BiOjRT4xNePAHDl1i6eG55mJFuiqzfOl568CweHnxRv4IPfG+LK3p1cc8Wt3Hvq\nIX44/UN6OxOcODdNyXIAh+H093k4d5T49fBI3hWGX3/RW9jc3fwne8br9719czrIKY8lYrzmmjt4\n6r5niGw+y7bBLj536gkA6hPbm851e2w/J4oHOVU9wsBAs60TXPfnJgG4fuceBjrT3N57M586/O/E\nzTjbtriReq/3VBFPROnocG2Trs4EAwPpIKrt7Eww0Juk4kWctq4HY6nVLTqSEQYG0gwM7OenS6/h\ns4fu4sCO65bt76kv2c2JjHuT7e5IL2q/ttHQgd50vK3vXjE9CEPQFU8z0OduHzH1RZ9PNGLgAB0d\n7o2lq6u94y83y3XMpYj6UeC4J+JHhRCTuJG6TxrILrSDTGZujvBiGBhIMz6ev/CG65QRb0LPduCZ\no2Ns35Ra0XM+errRC+Wr9x7nra+au9pLuVrn/R9/GMuLFiOmzisPbGdqssCDw4+gaTovGnSj5KrX\nrGpyqtg0Zsdx+OhXDjE04r734JNnOV98js8c+QKZyrTb2W4f/PN9dX77R18DuL/rbz/yHMamIX7u\n9gO8dM+1/NaHH+Dk1Gn+7699nZHiGM/bdAO1yW5giN6OKC/d/qM8fOaH/OuTX2FL7MepWzZPHx3D\n6DvPmcpRNsUHGB6Ks7Wzl5+/7SUka11zru2xU67gpuImRS9HfehclisGd6DXktA7zJEzp3n07JNE\nal3kxjsYPj/d8GpzcexyFyc4zolzw42MmhAnJk4T0U3MUoJxb2m2X7vh7XR3J4LxlGZcKyWTLYIX\nkVfKNcbH81ieiI+M5dAsK1gqbmQ8z/i4ux5q3XLQIdjfSwZuo/eGfvb2XLVsf09Jo5HCqFvmovbr\n2HbQy70j3t53o3X3RpcyOih4K13pmrbo8zENjVK5xtiEa3PWq/VV15XF/rte6AawlOeMtwMfBBBC\nbAU6gW8KIe7wPv9x4P4l7H/DE24QNdyi8nE5cRx33cyB7jjdqSgPPTPaMhVxLFPCsh1eeM1mPvSb\nt/HR330ZP/3yPeSqef5NfonPHf0Klt28yk54P0+NP8Nf3PcvPCbH2LXF/cOUQ1keGXmcyXKG7amt\n6LlBHMvguPMg2ZI7D1+rWxyq3kf0imf5/JlP88cP/jnRq56ksP1eRopjvGz7i/n5/W8KloEb7Oug\nM5rmjXtfQ9kqM9X5KOBwbHiCyA6JqZn8xk1vxxi+gfrwHnZ37Wp5XSY8gezvijc1BtM1HT2zEwyL\nf3r6U1iOxSb7akBjOtQDplKzqE9uxcHh++cemrN/y7Y4PzPKYMfmYIk1/v/2zjxMjru885+qvnt6\neu771EiakizJo1uyZFsGW7INGLAdGwPGj7lyLCHgZxMS2LALBDYhiwmQhNP4MQHbrIE4m8gYH0LG\ntiTLuqzLUumwpBlpLs1o7u7p6aP2j+qqrp7pubslTev3+Ufqq7p+1TVvvfV9L/RWtUpxIoBplVgi\no4uPRlVFGuPgjLa+xnxSq6QhSzJLixfjsk2uXU8V6wVrvA6N42GTZYqMXj2TpDMaFMfll1xXrim7\nTCed0cBhk0cFSq9dTf2nQL6iKK8D/xfdyH8O+KqiKLsAJ/Dr2e/itYu1QZS15WtGvmsgxNBwhLqy\nXG68voJgKJIyYGpo0vMq/OT7XGbGy+62fcS0GCPRETPnenRFaTQW5cm3n6Ul+jb5JcM8cl8TZQUe\nTp7vpXngPE6bk4caPsHQ8SbCLY1gD/P4oWfQNI1f7H4JqaiFHK2QdeWrCEaCBL3NEHHy4XkPcn/j\nB3HanObFryLeDndDxVqUggX0285jK2qjWdqP5Ayxpe5dFHuLqCjy0tETIBpLnaRldA4sznObVZZG\nVWm4swo0iQuDbdhlO/VuPQ3S+rsNh6NEuyvw2r389szLvNG2N2n7HYGLRLQoVb7KCX8fa+VjJEXx\nEegGX9O0hKYeT+EMmaPsMqsV+52Jro7T1dQBSuJ9+CcrPDLf7y3mhoo1bKxYax6L6RQeGRiaetjM\nU79GNXVVVUeAsR3zYdPMd0dgxeqpW/OpM0Fzh37rWVPqY92ScrbuPMerB9vYsFQP3nUMdVLgzqcz\n7rmWWmZxaprGzrY3zcene89Q568ZU1F6pPsYg1Hd8964zo0/x4lSW8Crh5tpH+qkIa8OtUUfXrws\nbyVv97dzmhM8d+YlXji/HS1i574FD7CmoZ6R6D38avdetu0dwNuQGBXW3h1AliRKC/T9kySJjy76\nI76661s46t4GWxRpJIctdbcAUF6Yw5m2Abr6hikrGFsBmTDqHvrizb8Gg2FiMY3hgIPCUCVB9wWW\nlyylNJgHdCb1VR8JRyHi5M+WfpIfHHmMXxz7FTbJxpr4DM0Lg3rfk8nGnTlsqTx1Oem1SCTG8EjU\nfN0IlJpFNVMIPs6GPJfFU5+RUfcAPZPmqBvIksyDi+8DoC1+MZ+Rp24GStOf0nglmNt7n+X0XEb5\npaVT1/NqSnMpzfewuK6AEy29tHT18cyJ/+Bru7/FT488ac4hNYwm6K1eOwNdiX7efWeBsfLLq+d3\nmZ8ZkPSsF6U2H9k7gIZGrb+a4+f09L57b56PvbUJYjLPn32ZGFHsF1ayql6XSZw2B9eVLICY3Yw9\naJpGW/cQpQWepA6DRZ5CmnJuQrJHkCSNsuBqHDbd6zY8+vEuml19Qew2iTyfkxy37gMNBcNmRk9J\naBll3lJuq900pv8LJDJ/6guq+ezyT+O2u/jZ279kR+tuNE2bulE3PfVE10tDkrHKL9aMF8NjN+6U\nppJRMhusnvpMMmqMtMapGnUrtrR46sKoCzKMcRtf5HfRcWl8iSAdtMRz4WvL9Bztm5sqkTwD/POh\nH/CH8zuRkDjSfYzmwWYkErfKgOmlv2/eFgpc+ZzuPYOmaUnyS0fgIsd7TiIHiiBq59yA3ohJqclH\nytG98xpfFcfO9ZDvc1JZnMPKunmMtOgNqcKtDawsX5qUg2w0xDK8tIFAmKHhiGmorazIX0Wks5pw\n2zzm5SQqKCuK9GBb+7hGfZgivxtZkszS+aFgxGwRUGgv43+u/0tqcqvIi2el9CWlEupTj2RJotZf\nzWeaPoXL5uKp47/hh4eeMAdaTGbU7ak8dXmspt5vGbBtyC+jhylnCqumPhP5ZfnCYpTaApbOK5z8\nzaMwjsV0OjQaOGwy0ZjGcLyXz7Wcpy7IML0DIWyyRGNNPpGoZk7fmYjdbfv4zv4fMjAyvYKlls5B\nvL4Ir3S8zPcO/Jj/uPQT3Mt2MMQlbq66gT9f/ikAOt1vUeB3msGkQDjAgc5DlHqKWZDfwPz8egbD\nQ3QGLibJL6/FvfSR9hqcI0V0BroYCgco9LvxFuhG2RkuZCAQZnFdAZIksXpRCdGOeiJH3kXkfCMr\nG5PLrUvyPdhkyfTUDeNensKo5+Y4CZ9dSqRFoczSHjjhqY+9EwqNRBkIhCmOe5BGoHRoOGw283K7\nEgqmP0d/vd/So10fZZf4M5uXV8uX1n6exoIFHOk+xpn+ZvJdeZO2ZrUWyZhGfbSnHk321EMjUcKR\nqHmnlGlP3Sq/eKdZfAT6BfZbn7uZ0hQy2GSYmvoMjLpxR2kMPXFk+OKXaeb23mc5PYMh02uFyYOl\nmqbxwrntnOx9hyeOPm0OvQVdX915pI1YbOzAiuFQhG7bSaRFf2Bby6uoPaewSTK2wVJs59bwIeVu\nFhUuZHFBI1pON/6yROrVno63CMcibKhciyRJzM8zJJgzyLJe2BGMhHijfS+5Dh/h7lJy0SfanI3P\n0rT7BtCiNvYd0o3zonijq8V1hXhcdsIBFx6XncV1yQVOdptMcb6H9u6ALr3EjXtl0djp79b+4Fbt\nvLTAgyxJKVvwGtOOiuOFMTmW7BfDU/cmGXVdNhjtqY/Wsos8hXx2+af4UDy4e13h2NTR0SRVlBry\ni+Gp26zyi+6dGwHsgUDYrOjNtAealP0yA099NhgXvekWH0Hi+AWGDU99bpvFub33WUxM0+gbHCE/\n12Uaqcl09dahdjoCnUhIHO85yQtntwN65eY3nzrAY1uPceRMd9JnBsNDfPfAj3E0HAEJ7mv8AI/e\n/Hd8feOXWBjewmBHkZmit77oZv0zeYfRNI0TPafZ1vwqsiSzrkIvUTB09dO9ZwHdOxx0nSMYGaap\nYCVoMsX2ckBvmzocCTEs9xIb8rPrqD4cYnG8OtVhl1m+QG8stWZxWUqts6LQy9BwhIFg2Gyjm9JT\nt/QSKStMGBy7TaakwENb15BZXGRgzXwB/Y/dbpN1Tz1u1N2W4Rdupx2XwzbGU3enGPYsSzI3V2/g\nWzd9lQeUu8e8Ppqk3i/jeeoWTd2IeQwEwglPPcV+pBPrMIh0VKlOB1NTn0bfFwPj+BkDToSmLsgI\nA0MjRGMaBT6X6albg3lPH/8N/7jnn5N6Y+/vPATA/Y0fpMCVz3NnXuSttuP8n18eML38CxZvP6bF\neOLo0zQHzhLtLeG9hQ9xS/VG3HZdG64p1fX15rjebhsuIHqpjCG5i2+8+W2+e+BHdA9f4pbqjaaX\nVp5Titfu4XR8TqbTqRHMO44sycx369PpK7xVgO6pG+mPWsCPpukGtNiSWbNpeSV2m8yWdanzyA0D\n3t4dMNvoGj1drHjddiRJvz0vzks2OI3VeQwNRzjyTvIFz5r5Anomjc9jj3vquqG0euqgp+P1DY3V\n1MfDJtuS8tPHI9lTHyelMRKjf0j31I1zZiA4YuapXw4P1BiWMRNNfTaYLRNm4KkbF0ejuOxazlMX\nZBAjSJqf66I4343dJpm6b1+on51tezg30MKb7Xp5uqZp7O88iF2yc/poLg3hTYDEY4eforW/k6b5\nerm5NSD4/NltHLt0gvxYNSMnVrK4IjlXuqZUN9RGELWjJ0j4wgIkJNqGOlhatJi/Wv3n3LvwLvMz\nsiTTkFdP1/AlekN9xMqPoTkDbK69heiwfrEo8ek9v8/2N9Mcl2ByNN0jHy2xNNbk8+O/uoWmxtTt\nS41gafulAO3dAfJ8TrzusZm6sqQb88rinDHe3K2rqgF4ae/5pOcv9ibLL6DLOIOWQKlnlFH3+5wM\nBMLxKToxojEtLR5ycqB0bD9147WBeBVvZbF+XJI89Qxr6gCLCxtpyKvDMUkv9XRjk2Xuu2U+d6yr\nnfZnR3vqc11+Ef3Ur1KMdMaCXBc2Waas0EtrXDs2Cn0AXjj3e9aVr6R1sIPOQBeR7nJeO60XDdnK\nGnHWHcd3/R7es2Iph97pNoOKb3erPH/mZQrdBbjeWY0shUxDYGBkwjTH+59c7A2iBXN5sOFhaorz\nx83YWJA/jyPdx3jx3HZCeafRgj7urL+Vl/fqXnmez0m9s5Y9HQfY23lQ/0xhHXtbAqaePlUMo97c\nMUB33/CYQR1WHrl/ecoOfrVluTTW5HP0zCVau4aoLM6hPzDCjsNtuJ02M0MG9Ha55y8OmQVI7tGe\nutdJTNMYDIZNfTcdxtSesvgohaYev0swJLuBQDhRVHMZjPo9C9+X8e8YjzvXp76bmwzj+AWFUReM\nJhKNTWsC+0QYuc7GRJ3KohwuXBziYk+QXW17cMh2mkqWsrfjLXa07OOFQ8cgB3KGa/jkA8uxyxL9\ngaWcGq5lx6WX+NdDj5Ff3URrqIeXmwd58dx2bJLMJ5c8yDd3vUNFkXfMbWdxnhuPy2Z66kaO+orK\nxgm9z/n59QD84fxO0GDkzBJkyWYW7+T7XNRHdaN+rr8Ft83FB1YvJc/eOibDZTIM+eXgqS40SDLA\nY95bOH5WxW2rqjnR0su2/ef52BaFX20/xdBwhA/fujDJ8zcCrt1xaWa0/OL36dp9/9CI+bl0eOqy\nJGGTpdTFR5bWu/2BMG6nzZz6NBgcQW9IS1IWjiCBkeo5lCXyizDqs6SjJ8CbxzrZc6yDtu4Af/nA\n8pRtaKdLz2AIydvH7/t+zfGjxZQX6k2ydp45QmewizVlK7mrYQv7Ow/x7IkXGbZFsGk2vnzPeyjI\nSRi2NdzG0q5KHj/yJMMVeon6s/o0Mx5Q7sETKyIUPmnq51YkSaKmxMfJC32MhKN09ATJ8zknNVI1\nudU4ZDvhWISC0CJaBwsIjcRMrTkvx8k8ai3vr6KqxMdHN088QT4VuR4HOW473f36RTBVjvpUWNFY\nTJHfxc7D7SybV8SOw+3Ulvp496qqpPflxGd3GtLM6CCoUeLeNzSSVk8ddG99vHF2kAiU5nodSYNB\njO+/HJ76XGR09osIlF7DvHmsgy/+6A2effUd2i8FiMY0nn3tzKy3OxwJcTj4Oq4lu7gw3Myejv10\nuvUg6K4LuwHYULmGIk8hq0pWMGLrR3YHWF62JMmgGywrvo5HVv4Z5doiwucX8N7Ke/jK+r/mpqr1\ntFjaA6SipjQXTdPb8Xb3D1OWP3kAzCHbWV5yPdW+SirDKwE9YGhk0fhznOZACNAHEM8USZKSsl0m\n8tQnwibLvHtlNaFwlH999jAS8LE7lDH6u5HWeHE8T9006qG0a9kOm0wkqo2RX+xJeeph/F6n2RTL\nmtJ4OTT1uYjhmWvo9zQzCbZeTQijPgteP6yXeD985yK+89kbuX5+ESdaelGbe2a0vbahDv791Fa+\nsuubdDqOooW8fHzxRyn2FHFo8A1sxec5FzxBsaeIhfFpNbUsx8jEsw5XGE2tv5qbCm8n0roA/0gd\nJV49cHq6Va/mbKj0p/xcTVxXP3CiC01jyoUhDy95gL9Z8zm8Tl0+CoWj9A6O4PM4sNtk7LKdGl+V\nuW+zoTxFMdFMuKmpEqddry7ctKKK+ZVjxwH43Mnyi8c9NvsFoH8onPbyfIdd1jX1MZOP5Ph36hlT\nuV4nOR4HkqTXJ4TCl6eidK5i9cwdDnlKs1WvZsSvPEMi0RgnW/ooL/Ryc1MlXreD922oB2DrzrNJ\n7+0Y6kxKPRxNV/AS/7T/B3x996Nsa36VmBbD3bMI28lNrK5o4o+XPYRTduBsOEKMCDdUrDFPvNNn\nIkQ7a8lz5LO0aOxwYyuGwWu3FNqcvtCPLEnUl49j1OMe/N54x0Zrz5fJkCQp0f9lRJ8glO9L5Isv\nLV6MU3YwP69+yttMhWHUXQ6b2X9lJvg8Dm5fW0t1SQ73bko94s3w1I1+8qPll1Seeqo89Zlgt8U1\n9dgo+SX+b3e/fqHJ9TqQJSmeqSM89cmwGvW53iIAhKaexCtvXaBvcIT3b6yf9Gp9rn2AUDialK2x\noCqPxXUFHD3bw+nWPuZX5vH8mW1sPfMCfmcuN1auY2PVOvJdCQ/waPdxnjj6NIFIkMWFjWyoXMuy\n4uv4/Hd3Uhg3UFW+Ch5cfD+PH30STYNl+bpHHonGOHiqmxzncr7xwMZJ99ma/md8/mz7ANWlOePq\n5FXFOciSZOZsT8eoQ8KgDQRHCIYi5PkSF48tdbewqXrDrAtVyuN56eVF3ll7WXff3MDdN48/s9Na\nmepy2MbIM1ZN3WjmlS4t22G3EQyNmJ66bVQ/9R7TqDvNf/uHRgjlXb7sl7mIw5LcMNf1dLhGjXrf\n0AhOu5yUY9zVG+TJF08QjWmUFXpYf135mM/FtBhHu49zsvcdWlpD2IqH8Jf5iMaiZgHJXRvqOXau\nh607znL9+gHToIdjYX579mV+d+73VOSUUe4txWVzsqttLzbZxkcW3cvGynWA7tUGQxHyLZLIqrIm\n9p1qZY/azk65l/veVcyJll4CoQg3LKmekjHz5zjxuGymUW/uGCQSjTG/avypg06HjfIir1m8lKo9\n7UQY3qHRhz3f0oFPluS0VB5Wl+pGvaZkahPsZ0OORW7xuMYayYT8MpJBT12fEmS0AjDll3iLAH9c\nT/d5HLR1DTEcyo4+4ZkiSX7JgmN0zRn1CxcH+frP95Hvc/GVh9eYHup/7TxLNKYhAU++eILFdYXm\nH2gkFmFP+wFebtan7Bg4G+ClS0fY98Y2NlVv0Acy1OazsDqPI32HOHHyMG7Zy4MNDzO/tJQ9HQd4\no20vrYNtZsvVQncBn176sSRd2Sg8MtIZDT6+bgsnD+1m2/7zbFlby4ETXYCeuTEVJEmivNBLS+cg\n0ViM0xd0PX3+OHq6QU2pzzTqJVMIlFpxjjLqeb6ZyyPjUVbg5fP3NZl59ZnE6qmPLjwC3Zv2uOz0\nDY1kUFPXklJnHaPSaBOeugMNuDQwjDMLtOJMkSy/CKN+VaNpWtKJHBiO8C//flhvBXspwG/+cJqP\nbG6ksyfAjsPtVBR52bS8il9uO8nPX1D5zN1LuTDQwfcP/Bt90S69x0n5KtaUruS7zx4gNxdWrrDz\nZts+nj31HP/v9PPYJBtalYZTi6BFHPQeW8Gjb5xk8+phHrh1HTdVrSemxegZ7qV7+BI1udVjvFUj\nRz1/lD7ssNu4/7ZGfvCbQzy36ywHTl3E67LTWDN+wc1orEMhjCDpRJ46QG2pj91vd5DrdaSs1pwI\nd9ygdcTvDqY61Wa6XB+vmM00OZMYddDXaPXU0ya/xFvEhqOxpAwNWZbM+Z4AufFukYZx7xkImVOb\nBGNJ9tTnvkSVtUb92LkeHtv6NhVFXu7dNJ+68lwe2/o2HT1BbltVzdGzl3h533lWNpbw+uE2YprG\nB26cx+pFpRw4cZH9Jzr59ktbeUfeCXKUSGc19y+5k1uvW8jJ872MdBfTVFfFhxWF9zfcwc7WNzl4\n8ahe6SmBS3ayJm8TwSIfv99/npf2tiBJ8KF3L0CWZIo8hRR5UveNNj31FEG/zWvreOalE2zbdx5N\ng/VLyqZV8GTtlXL6Qh8+jyNpilEqjGDpdPV0SBTeGBOT8nyZMeqXC6tx9Iwjq/hznHRcCpitXNMn\nv+i/8/BIZMxv7rDLRON3BrmeuKcevwBpmgiSToTw1K9yNE3j9/sv8PTLJ9E0jZ6BEG+f3UttmY/m\njkEW1xXwoVsXcLZ9gP/98738cPsrBLReSup9uIu7ebO9hdJlZziXp/KOvR8tYqcx+m6ONrt4sfci\nm66bb07nWRQvMspxeNlcdwub4yPSRrN6USn/+NR+XtzTgixLvPeGOoaCYQaDEfoGQ/QOhugbGqFp\nQTHzKvxjqkmtOOwyd22s54nnjwOwcuH0KjAr4sFStbmX7v4QTfOLJr0tr6/w43baaKiY2KNPhWFM\njGKd/AzIL5cTh13G5bARCkcn9NQ1MPvfp1N+AQiGomPy4x122QzMGhk4PssAZ5HOOD6jUxrnOlll\n1IOhCE9vO8nrh9rwex38t7uXMRIN88wb+2gdPoy/xs09WxQ9wOTtoWzNQfpoxwEMAj86vMfcls0r\nUyDV8PFV9zGvqJwnpRNs23+eVw+2cry5F2DCPiNW8nKcfOHDK/jmUwf43e5mfre7OeX7th+4wNc+\nsTap70sqNiwt57dvnKN3MMSSaU6JMTJgdh1tByaXXkDXkf/hT24Y14hNhGHQjCrIue6pA/g89kmN\nOiTuTtLV8tb01EMRMxhqYDVMhu6f603O1BGkxjFBfGIuklajriiKDHwfaAJCwKdUVT2Vzu9IRUzT\n2HWknV+/cpq+oRCV1TFWr4nw24u/5Ex/M5HSCA4gDHz74EFy7F6GIrrG6wpU4QvVsmV9BcHIMA6b\ng3p/DVW+Shxy4vC8b2M9rx9u4z93nCUYilBdkpPUo3sy8nwuvvCRFTyz/RShkSg5bl2fzve5TM3B\n5QAACuJJREFUyPc5udA1xHO7zvHYc8dMHXq0pm5gt8l84cMrCAxHpm1oywo9SGCW7E/FqMPM5kbC\nWIOWKU39cpLjdtDdHxr32BvH6mI8OJxuT12DsfJL/LHHZTffZz0/RTrj+CTJL1lwnNLtqX8QcKuq\neoOiKOuBR4EPpPk7CI6M8NT+bXSNdNEbusRQdICwFkZaqJHrjNHDCC9dAAmJSl85C/Ln0ZBXTyAc\n5HTfGc70naPKV8Fd829nnr9uSlkBeTlObl9bw3/uOAskpJfpkO9z8cd3LUn5mqZpNHcMcvidbmzx\nwFeud/zgVqHfTeHESSspcdhtFOW56eobRpJgXkXu5B+aBVaD5nLacDvn/s2hESydzFM3NPW09X6x\nGPLRRt1uGvLEOZPrEZ76VLAGR0VK41huBH4HoKrqG4qirJ7ozQUFXuwziDa/cOgA+wP6VB9soGkO\nXDYneV43LoeDhoJalpcvoal8MX53+ozWR99zHa+81Ur/0Ahrl1VSUpJeg/iFh9bw2Ue30zsQoqTA\nQ1lpaqs92++trfDT1TdMfYWfmqrZNx+bEHviFCvyu2e17+k+3jOlKN8D53ooLsxJuU81o9oLVFfm\nzWjKPSSv2W+JR7jd9qTXPPEAblGex3xetlxAc33Oq+b4Tcbl3k93TuK45uXO7hydDen63nQbdT/Q\nZ3kcVRTFrqpqyhr5np7UE9wno6msgT8a+BjFeX78ci6lfn9Kryk0ABcHBlJsYeZ85LaFvHqwlepC\nNxcvpnfbAJ+4cxHffuYgBT5Xyu2XlOTO+nsL47p2XdnstzUZg8HEdHufxzHj70vHutOFMRs0Fo6k\n3qdoNPFem8ylSxOPIRyP0WuOhBPbJaYlvSbFGwC5HbL5fDgSG/f9VytX4ncOWY5rNBy9Isdpuuue\n6AKQbqPeD1i/TR7PoM8GmyzzrsZlV+QEWLu4jLWLyzK2/aUNRTxyf9OsephMRl2Z/hMtmmKgdzZY\n0/nysyBICnqgFCaSXxK/XTp7mNvtCZlwdCfBVDq6wy7jcdkIhsYOvxYksAZH7UJ+GcMO4C7gmbim\nfjjN278mWNaQ2UKaG5aUU+h3XxajbrfJZmHMTIOtVxuNNfn4va3Ulaf2lqy6drpy1GGU8RknUOrP\nSY7D5HqcBENBMSBjAqzFWyJPfSzPApsVRdmJ3pr442neviANyLI0ZhZoJnE5bARCkTmfo26wdF4R\n3/mLm8Z93W6TzQ6J6fSQrV5kquIjSBQeGfi8Djp7g1nRfTCTGMVbwqiPQlXVGPCn6dymYO7jcupG\nPRvSGadKXo6TwWA4rZ66fQKZwDTqYzx1/XG6cuWzFWe8eMuRBTLV3L8sCa56jHS6bPHUp4IhNaUz\nldCabmeXJ9fUrY+zwQPNJMbxy4bjNPdXILjqMQxbNlSTThXjriSd8stEAT1/jgsJfVi4FUPfF3nq\nE2OkVos8dYFgChi3/teS/GJ46mkNlCZ56snG5461tTQtKBrT797w1IX8MjGOUaMB5zLCqAsyziql\nBL9lwv21QCY89WRNPVl+8brtKWeqrruujPZLgYxnVM11EvLL3L/4CaMuyDibV9eweXXNld6Ny4rp\nqWfKqE+xQrUg18XDd048u1aQMOrZ4KnP/RUIBFchRvzAnWLk3UxxTJDSKJgdRoA0G1oUz/0VCARX\nIYtqC7hzXS03LqtI2zaTi4/EaLp0YnrqWXCxFPKLQJAB7DaZ+961IK3bFJ565jCNehZkCQmjLhDM\nEazeuTDq6WVZQxG9gyMUZrDn0uVCGHWBYI6Q3CZAyC/pZOOyCjamUSq7kojLvUAwR5iooZdAYCDO\nDIFgjjBRQy+BwECcGQLBHEFkvwimgjDqAsEcQWS/CKaCODMEgjmCTRbZL4LJEWeGQDBHkCTJ9NaF\n/CIYD2HUBYI5hOGhC09dMB7izBAI5hAJT1386QpSI84MgWAO4YjLLkJ+EYyHMOoCwRxCyC+CyZhR\nmwBFUSTgPHAy/tQuVVW/qCjKeuC7QAR4UVXVr6ZnNwUCAQj5RTA5M+39Mh/Yr6rqXaOe/yFwL/AO\n8JyiKCtUVT0wmx0UCAQJEp66kF8EqZmpUV8FVCmKsh0IAo8AbYBLVdXTAIqivADcBoxr1AsKvObA\n15lSUpI7q8/PRa7FNcO1ue7Ra/bGRwKWl/nxebNz5uu1+DtD+tY9qVFXFOWT6EbbymeAv1dV9VeK\notwI/AK4G+i3vGcAaJho2z09gent7ShKSnK5eHFgVtuYa1yLa4Zrc92p1qzFNAB6ewIEh0JXYrcy\nyrX4O8P01z3RBWBSo66q6k+Bn1qfUxTFi66bo6rq64qiVKIbces35QK9U95LgUAwKTdeX0FZoReX\nc+4PcxBkhplGW/4X8HkARVGagBZVVfuAEUVR5scDqbcDr6VnNwUCAcANS8p56HblSu+G4Cpmppr6\nPwC/UBTlvege+8Px5/8UeBKwoWe/7J71HgoEAoFgyszIqKuq2gO8N8XzbwDrZ7tTAoFAIJgZItlV\nIBAIsghh1AUCgSCLEEZdIBAIsghh1AUCgSCLEEZdIBAIsghh1AUCgSCLkDRNu9L7IBAIBII0ITx1\ngUAgyCKEURcIBIIsQhh1gUAgyCKEURcIBIIsQhh1gUAgyCKEURcIBIIsQhh1gUAgyCJm2k/9iqEo\nigx8H2gCQsCnVFU9dWX3KjMoiuIAHgfqARfwdeBt4AlAA44An1FVNXaFdjFjKIpSCuwDNqP37H+C\n7F/zF4H3A070c/wPZPG64+f3z9DP7yjwabL8t1YUZR3wTVVVb1EUZQEp1qooyqeBP0E/Fl9XVXXr\ndL5jLnrqHwTcqqreAPwN8OgV3p9M8iDQrarqTcAdwL8A3wb+Nv6cBHzgCu5fRoj/sf8Ifag5XBtr\nvgXYAGwENgE1ZP+63wPYVVXdAHwN+AZZvGZFUb4APAa440+NWauiKOXAX6CfB7cDf68oims63zMX\njfqNwO/AHMqx+sruTkb5FfDl+P8l9Cv3KnQPDuB54LYrsF+Z5lvAD4HW+ONrYc23A4eBZ4H/AraS\n/es+Adjjd99+IEx2r/k0cI/lcaq1rgV2qKoaio8IPQVcP50vmYtG3Q/0WR5HFUWZczLSVFBVdVBV\n1QFFUXKBXwN/C0iqqhq9HQaAvCu2gxlAUZSHgYuqqr5geTqr1xynGN1BuY/EWEg5y9c9iC69HAd+\nAnyPLP6tVVX9DfqFyyDVWkfbt2kfg7lo1PuBXMtjWVXVyJXamUyjKEoNsB34uaqqTwFWfTEX6L0i\nO5Y5PgFsVhTlFWA58G9AqeX1bFwzQDfwgqqqI6qqqsAwyX/M2bjuR9DX3IgeI/sZejzBIBvXbCXV\n3/Jo+zbtYzAXjfoOdC0ORVHWo9+yZiWKopQBLwJ/rarq4/GnD8T1V4A7gdeuxL5lClVVb1ZVdZOq\nqrcAbwEPAc9n85rjvA7coSiKpChKJZADbMvydfeQ8EovAQ6y/PweRaq1vgncpCiKW1GUPGAxehB1\nysxF2eJZdE9uJ7rO/PErvD+Z5EtAAfBlRVEMbf1zwPcURXECx9BlmWznvwM/yeY1q6q6VVGUm9H/\nqGXgM8AZsnvd/wQ8rijKa+ge+peAvWT3mq2MOa9VVY0qivI9dAMvA/9DVdXh6WxUtN4VCASCLGIu\nyi8CgUAgGAdh1AUCgSCLEEZdIBAIsghh1AUCgSCLEEZdIBAIsghh1AUCgSCLEEZdIBAIsoj/D7tf\nvOsaO+whAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1146cb160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu=1\n",
    "sigma=1\n",
    "n_series=100\n",
    "n_points = 100\n",
    "n_features = 2\n",
    "# multi-variate time-series: \n",
    "# x_t ~ [Gaussian(mu*t,sigma*t), Gaussian(mu*t,sigma*t/4)]\n",
    "# Points X Time X Features\n",
    "x_data = np.zeros((n_series,n_points,n_features), dtype=np.float32)\n",
    "for i in range(n_series):\n",
    "    x_data[i] = gen_time_series(n_points)\n",
    "\n",
    "if TESTING: x_data_test=x_data[:10]\n",
    "plt.title(\"Sample bi-variate time series\")\n",
    "for i in range(n_features):\n",
    "    plt.plot(x_data[0,:,i]);\n",
    "print(x_data.shape)\n",
    "print(reduce(lambda x,y: x*y, list(x_data.shape)), \"data values\")\n",
    "\n",
    "x_train = x_data[:70]\n",
    "x_validate = x_data[70:90]\n",
    "x_test = x_data[-10:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building Variational RNN\n",
    "\n",
    "note that all gaussian distribtutions: p(z), q(z|x), p(x|z) should have the same dimensionality, i.e. share z_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First I'm going to define the loss functions necessary because they don't depend on anything\n",
    "def negative_log_likelihood_gaussian(y, mu, sigma):\n",
    "    \"\"\"\n",
    "                a                  b              c\n",
    "    (1/2)*[ ((y-mu)/sigma)^2 + 2*log(sigma) + log(2*pi) ]\n",
    "    learn more: https://www.statlect.com/fundamentals-of-statistics/normal-distribution-maximum-likelihood\n",
    "    \"\"\"\n",
    "#     sigma might be really small. in that case, a & b blow up\n",
    "    sigma=tf.clip_by_value(sigma, 1e-3,1e10)\n",
    "    a = tf.square(tf.div(tf.subtract(y,mu), sigma))\n",
    "    b = 2*tf.log(sigma)\n",
    "    c = tf.log(2*np.pi)\n",
    "    return tf.add(tf.add(a,b),c)\n",
    "\n",
    "def kl_gaussian_gaussian(mu1, sigma1, mu2, sigma2):\n",
    "    \"\"\"\n",
    "            a                  b         c           d\n",
    "    log(sigma2/sigma1) + (sigma1^2 + (mu1-mu2)^2)/2sigma2^2 - 1/2\n",
    "    \"\"\"\n",
    "    # Note: if sigma1 or sigma2 go to 0, this explodes\n",
    "    sigma1=tf.clip_by_value(sigma1, 1e-3,1e10)\n",
    "    sigma2=tf.clip_by_value(sigma2,1e-3,1e10)\n",
    "    a = tf.log(tf.divide(sigma2,sigma1))\n",
    "    b = tf.square(sigma1)\n",
    "    c = tf.square(tf.subtract(mu1,mu2))\n",
    "    d = 2*tf.square(sigma2)\n",
    "    e = tf.divide(tf.add(b,c), d)\n",
    "    return tf.add(a,e) - 0.5\n",
    "\n",
    "# I will also define a simple wrapper around reshape to flatten the 1st two dimensions of a tensor\n",
    "def flatten(x):\n",
    "    return tf.reshape(x, [int(x.shape[0])*int(x.shape[1]), -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class VRNN(object):\n",
    "    \"\"\"class for VRNN.\"\"\"\n",
    "    def __init__(self, input_shape, rnn_size, x_size, x_1_size, \n",
    "               z_size, z_1_size, phi_size, prior_size, theta_size, lr=.01):\n",
    "        self.batch_size = input_shape[0]\n",
    "        self.ntime_steps = input_shape[1]\n",
    "        self.nfeatures = input_shape[2]\n",
    "        \n",
    "        self.rnn_size = rnn_size\n",
    "        self.x_size = x_size\n",
    "        self.x_1_size = x_1_size\n",
    "        self.z_size = z_size\n",
    "        self.z_1_size = z_1_size\n",
    "        self.phi_size = phi_size\n",
    "        self.prior_size = prior_size\n",
    "        self.theta_size = theta_size\n",
    "        self.learning_rate=lr\n",
    "        \n",
    "        self.inputs = self._define_input_placeholder(input_shape)\n",
    "        self._build()\n",
    "\n",
    "    def _define_input_placeholder(self, input_shape):\n",
    "        return tf.placeholder(tf.float32, shape=input_shape, name=\"x_t\")\n",
    "\n",
    "    def _build(self):\n",
    "        self._init_cell()\n",
    "        # transpose input so T X N X F instead of N X T X F\n",
    "        # this allows you to do one time-step at a time for all batches in tf.scan\n",
    "        self.inputs_t = tf.transpose(self.inputs, [1, 0, 2])\n",
    "        self.outputs = self._scan_sequence(self.inputs_t)\n",
    "        \n",
    "        with tf.name_scope(\"negative-log-likelihood\"):\n",
    "            self.nll = negative_log_likelihood_gaussian(self.inputs_t, \n",
    "                                                        self.outputs[-2], self.outputs[-1])\n",
    "            self.nll_r = flatten(self.nll)\n",
    "            tf.summary.scalar(\"\", xent)\n",
    "        with tf.name_scope(\"kl-Divergence\"):\n",
    "            self.klgg = kl_gaussian_gaussian(self.outputs[5], self.outputs[6], \n",
    "                                             self.outputs[8], self.outputs[9])\n",
    "            self.klgg_r = flatten(self.klgg)\n",
    "\n",
    "        with tf.name_scope(\"objective-function\"):\n",
    "            self.objective = tf.add(tf.reduce_mean(self.nll_r),tf.reduce_mean(self.klgg_r))\n",
    "        \n",
    "        with tf.name_scope(\"train\"):\n",
    "            self.train_step = tf.train.AdamOptimizer(self.learning_rate).minimize(self.objective)\n",
    "\n",
    "    def _init_cell(self):\n",
    "        with tf.variable_scope(\"lstm\"):\n",
    "            try:\n",
    "                self.cell = tf.contrib.rnn.LSTMCell(rnn_size, state_is_tuple=True)\n",
    "            except:\n",
    "                self.cell = tf.contrib.rnn.LSTMCell(rnn_size, state_is_tuple=True, reuse=True)\n",
    "        self.state = self.cell.zero_state(batch_size, tf.float32)\n",
    "    \n",
    "    def _calculate_prior(self, hidden_state):\n",
    "        \"\"\"\n",
    "        calculate your prior based on your previous state information, i.e. z ~ p(z)\n",
    "        \"\"\"\n",
    "        with tf.variable_scope(\"prior\"):\n",
    "            prior = fc_act(hidden_state, prior_size, act=tf.nn.relu, name=\"fc_prior\")\n",
    "        with tf.variable_scope(\"prior_mu\"):\n",
    "            prior_mu = fc_act(prior, z_size, name=\"fc_prior_mu\")\n",
    "        with tf.variable_scope(\"prior_sigma\"):\n",
    "            prior_sigma = fc_act(prior, z_size, act=tf.nn.softplus, name=\"fc_prior_sigma\")\n",
    "        return [prior, \n",
    "                prior_mu, \n",
    "                prior_sigma]\n",
    "        \n",
    "    def _encode_x(self, x_t, hidden_state):\n",
    "        \"\"\"\n",
    "        infer the latent values given z, i.e. z ~ p(z|x)\n",
    "        \"\"\"\n",
    "        with tf.variable_scope(\"x_1\"):\n",
    "            x_1 = fc_act(x_t, self.x_1_size, act=tf.nn.relu, name=\"fc_x_1\")\n",
    "        with tf.variable_scope(\"phi\"):\n",
    "            phi = fc_act(tf.concat([x_1, hidden_state], axis=1), self.phi_size, act=tf.nn.relu, name=\"fc_phi\")\n",
    "        with tf.variable_scope(\"phi_mu\"):\n",
    "            phi_mu = fc_act(phi, self.z_size, name=\"fc_phi_mu\")\n",
    "        with tf.variable_scope(\"phi_sigma\"):\n",
    "            phi_sigma = fc_act(phi, self.z_size, act=tf.nn.softplus, name=\"fc_phi_sigma\")\n",
    "        epsilon = tf.random_normal(shape=[x_t.get_shape().as_list()[0], self.z_size], seed=GLOBAL_SEED)\n",
    "        # z = mu + epsilon*sigma\n",
    "        z_t = tf.add(phi_mu, tf.multiply(epsilon, phi_sigma))\n",
    "        return [x_1, \n",
    "                phi,\n",
    "                phi_mu,\n",
    "                phi_sigma,\n",
    "                z_t]\n",
    "    \n",
    "    def _decode_z(self, z_t, hidden_state):\n",
    "        \"\"\"\n",
    "        generate the mean and variance for your data from latent values z, i.e. x ~ q(x|z)\n",
    "        \"\"\"\n",
    "        with tf.variable_scope(\"z_1\"):\n",
    "            z_1 = fc_act(z_t, self.z_1_size, act=tf.nn.relu, name=\"fc_z_1\")\n",
    "        with tf.variable_scope(\"theta\"):\n",
    "            theta = fc_act(tf.concat([z_1, hidden_state],axis=1), self.theta_size, act=tf.nn.relu, name=\"fc_theta\")\n",
    "        with tf.variable_scope(\"theta_mu\"):\n",
    "            theta_mu = fc_act(theta, self.x_size, name=\"fc_theta_mu\")\n",
    "        with tf.variable_scope(\"theta_sigma\"):\n",
    "            theta_sigma = fc_act(theta, self.x_size, act=tf.nn.softplus, name=\"fc_theta_sigma\")\n",
    "        return [z_1,\n",
    "                theta,\n",
    "                theta_mu,\n",
    "                theta_sigma]\n",
    "    def _rnn_recurrence(self, x_1, z_1, cell_state, hidden_state):\n",
    "        \"\"\"\n",
    "        use x_1 and z_1 and the previous hidden state to update lstm cell\n",
    "        \"\"\"\n",
    "        try:\n",
    "            return self.cell(tf.concat([x_1, z_1], axis=1), [cell_state, hidden_state], \"rnn\")\n",
    "        except:\n",
    "            tf.get_variable_scope().reuse_variables()\n",
    "            return self.cell(tf.concat([x_1, z_1], axis=1), [cell_state, hidden_state], \"rnn\")\n",
    "\n",
    "        \n",
    "    def _vrnn_step(self, vals, x_t):\n",
    "        \"\"\"\n",
    "        A single step on the time-series (or batch) using the VRNN. \n",
    "        \"\"\"\n",
    "        # I chose to store them all but you don't HAVE to do this. \n",
    "        # This is a convenient thing to do however so you can check values later on\n",
    "\n",
    "        prev_hidden = vals[1]\n",
    "        # mu_z, sigma_z = f(x, ht-1)\n",
    "        encoding_vals = [x_1, \n",
    "                phi,\n",
    "                phi_mu,\n",
    "                phi_sigma,\n",
    "                z_t] = \\\n",
    "                    self._encode_x(x_t, prev_hidden)\n",
    "        \n",
    "        # mu_phi, sigma_phi = f(ht-1)\n",
    "        prior_vals = [prior, \n",
    "                prior_mu, \n",
    "                prior_sigma] = \\\n",
    "                    self._calculate_prior(prev_hidden)\n",
    "\n",
    "        # mu_x , sigma_x = f(z, ht-1)\n",
    "        decoding_vals = [z_1,\n",
    "                theta,\n",
    "                theta_mu,\n",
    "                theta_sigma] = \\\n",
    "                    self._decode_z(z_t, prev_hidden)\n",
    "            \n",
    "        output, state = self._rnn_recurrence(x_1, z_1, vals[0], vals[1])\n",
    "\n",
    "        return [\n",
    "            state[0], # vals[0],\n",
    "            state[1], # vals[1],\n",
    "            output, # vals[2],\n",
    "            x_1, # vals[3],\n",
    "            phi, # vals[4],\n",
    "            phi_mu, # vals[5],\n",
    "            phi_sigma, # vals[6],\n",
    "            prior, # vals[7],\n",
    "            prior_mu, # vals[8],\n",
    "            prior_sigma, # vals[9],\n",
    "            z_t, # vals[10],\n",
    "            z_1, # vals[11],\n",
    "            theta, # vals[12],\n",
    "            theta_mu, # vals[13],\n",
    "            theta_sigma # vals[14]\n",
    "        ]\n",
    "   \n",
    "    def _scan_sequence(self, sequence):\n",
    "        with tf.name_scope(\"scan_input\"):\n",
    "            return tf.scan(self._vrnn_step, sequence, initializer=[\n",
    "                    self.state[0], # cell_state - 0,\n",
    "                    self.state[1], # hidden_state - 1,\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.rnn_size), dtype=np.float32)), # cell_output - 2\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.x_size), dtype=np.float32)), # x_1 - 3\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.phi_size), dtype=np.float32)), # phi - 4\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.z_size), dtype=np.float32)), # phi_mu - 5\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.z_size), dtype=np.float32)), # phi_sigma - 6\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.prior_size), dtype=np.float32)), # prior - 7\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.z_size), dtype=np.float32)), # prior_mu - 8\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.z_size), dtype=np.float32)), # prior_sigma - 9\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.z_size), dtype=np.float32)), # z_t - 10\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.z_size), dtype=np.float32)), # z_1 - 11\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.theta_size), dtype=np.float32)), # theta - 12\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.x_size), dtype=np.float32)), # theta_mu - 13\n",
    "                    tf.Variable(np.zeros((self.batch_size, self.x_size), dtype=np.float32)) # theta_sigma - 14\n",
    "            ])\n",
    "    \n",
    "    @property\n",
    "    def loss(self): return self.objective\n",
    "    \n",
    "    @property\n",
    "    def x(self): return self.inputs\n",
    "\n",
    "    @property\n",
    "    def reconstructions(self): \n",
    "        \"\"\"\n",
    "        mean and variance to reconstruct x\n",
    "        transposed from T X N X F to be N X T X F\n",
    "        \"\"\"\n",
    "        return tf.transpose(self.outputs[-2], [1, 0, 2]), tf.transpose(self.outputs[-1], [1, 0, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# test computing latent variable\n",
    "batch_size=10\n",
    "\n",
    "rnn_size=20\n",
    "x_size=x_data.shape[2]\n",
    "x_1_size=x_data.shape[2]\n",
    "\n",
    "z_size=15\n",
    "z_1_size=15\n",
    "\n",
    "phi_size = 25   # N ~ q(z|x)\n",
    "prior_size = 25 # N ~ p(z)\n",
    "theta_size = 25 # N ~ p(x|z)\n",
    "\n",
    "tf.reset_default_graph()\n",
    "vrnn = VRNN([batch_size, n_points, n_features], rnn_size, x_size, x_1_size, \n",
    "       z_size, z_1_size, phi_size, prior_size, theta_size)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define session and model saver\n",
    "def get_new_session():\n",
    "    sess = tf.Session()\n",
    "    sess.run(tf.global_variables_initializer())\n",
    "    return sess\n",
    "sess = get_new_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import os.path\n",
    "import shutil\n",
    "LOGDIR = \"/tmp/test1/\"\n",
    "if os.path.exists(LOGDIR):\n",
    "    shutil.rmtree(LOGDIR)\n",
    "writer = tf.summary.FileWriter(LOGDIR)\n",
    "writer.add_graph(sess.graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO: \n",
    "# scalar summaries: train loss, validation loss, test loss\n",
    "# histogram summaries: p(z), p(z|x), p(x|z)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Optimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Optimizer(object):\n",
    "    def __init__(self, loss, initial_learning_rate, num_steps_per_decay, decay_rate, max_global_norm=1.0):\n",
    "        \"\"\" Create a simple optimizer.\n",
    "        \n",
    "        This optimizer clips gradients and uses vanilla stochastic gradient\n",
    "        descent with a learning rate that decays exponentially.\n",
    "        \n",
    "        Args:\n",
    "            loss: A 0-D float32 Tensor.\n",
    "            initial_learning_rate: A float.\n",
    "            num_steps_per_decay: An integer.\n",
    "            decay_rate: A float. The factor applied to the learning rate\n",
    "                every `num_steps_per_decay` steps.\n",
    "            max_global_norm: A float. If the global gradient norm is less than\n",
    "                this, do nothing. Otherwise, rescale all gradients so that\n",
    "                the global norm because `max_global_norm`.\n",
    "        \"\"\"\n",
    "        trainables = tf.trainable_variables()\n",
    "        grads = tf.gradients(loss, trainables)\n",
    "        grads, _ = tf.clip_by_global_norm(grads, clip_norm=max_global_norm)\n",
    "        grad_var_pairs = zip(grads, trainables)\n",
    "        \n",
    "        global_step = tf.Variable(0, trainable=False, dtype=tf.int32)\n",
    "        learning_rate = tf.train.exponential_decay(\n",
    "            initial_learning_rate, global_step, num_steps_per_decay,\n",
    "            decay_rate, staircase=True)\n",
    "        optimizer = tf.train.GradientDescentOptimizer(learning_rate)\n",
    "        self._optimize_op = optimizer.apply_gradients(grad_var_pairs,\n",
    "                                                      global_step=global_step)\n",
    "    \n",
    "    @property\n",
    "    def optimization_op(self):\n",
    "        \"\"\" An Operation that takes one optimization step. \"\"\"\n",
    "        return self._optimize_op"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# writer = tf.summary.FileWriter(\"/tmp/tensorflow/\", graph=tf.get_default_graph())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tf.summary.scalar('dropout_keep_probability', keep_prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define optimizer\n",
    "def get_optimizer():\n",
    "    train_op = tf.train.AdamOptimizer()\n",
    "    try:\n",
    "        optimizer = train_op.minimize(tf.reduce_mean(vrnn.loss))\n",
    "    except ValueError:\n",
    "        with tf.variable_scope(tf.get_variable_scope(),reuse=False): \n",
    "            optimizer = train_op.minimize(tf.reduce_mean(vrnn.loss))\n",
    "    return optimizer\n",
    "optimizer = get_optimizer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_loss_on_data(sess, data, model, batch_size):\n",
    "    losses = []\n",
    "    nbatches=max(data.shape[0]//batch_size,1)\n",
    "    for j in range(nbatches):\n",
    "        batch = x_data[j*batch_size:(j+1)*batch_size]\n",
    "        loss = sess.run([model.loss], feed_dict={model.x: batch})\n",
    "        losses.append(loss)\n",
    "    return np.array(losses).mean()\n",
    "\n",
    "def train(sess, data, model, optimizer, batch_size):\n",
    "    losses = []\n",
    "    nbatches=max(data.shape[0]//batch_size,1)\n",
    "    for j in range(nbatches):\n",
    "        batch = x_data[j*batch_size:(j+1)*batch_size]\n",
    "        _, loss = sess.run([optimizer, model.loss], feed_dict={model.x: batch})\n",
    "        losses.append(loss)\n",
    "    return np.array(losses).mean()\n",
    "\n",
    "def reconstruct(sess, data, model, batch_size):\n",
    "    mus = []\n",
    "    sigmas = []\n",
    "    nbatches=max(data.shape[0]//batch_size,1)\n",
    "    for j in range(nbatches):\n",
    "        batch = x_data[j*batch_size:(j+1)*batch_size]\n",
    "        reconstructions = sess.run([model.reconstructions], feed_dict={model.x: batch})\n",
    "        mu, sigma = reconstructions[0]\n",
    "        mus.append(mu)\n",
    "        sigmas.append(sigma)\n",
    "    mus=mus[0]\n",
    "    sigmas=sigmas[0]\n",
    "    return mus, sigmas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "def train_model(N, sess, x_train, x_validate, vrnn, optimizer, batch_size):\n",
    "    train_losses = []\n",
    "    validate_losses = []\n",
    "    # train the vrnn for N epochs and see the training and validation losses\n",
    "    for _ in range(N):\n",
    "        train_loss = train(sess, x_train, vrnn, optimizer, batch_size)\n",
    "        validate_loss = compute_loss_on_data(sess, x_validate, vrnn, batch_size)\n",
    "        train_losses.append(train_loss)\n",
    "        validate_losses.append(validate_loss)\n",
    "    return train_losses, validate_losses\n",
    "\n",
    "train_losses, validate_losses = train_model(100, sess, x_train, x_validate, vrnn, optimizer, batch_size)\n",
    "plt.plot(train_losses);\n",
    "plt.plot(validate_losses);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# see how well the vrnn is making predictions\n",
    "mus, sigmas = reconstruct(sess, x_test, vrnn, batch_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_reconstruction(mus, sigmas, data, index):\n",
    "    mu_series=mus[index]\n",
    "    sigma_series=sigmas[index]\n",
    "    actual = x_test[index]\n",
    "\n",
    "    reconstruction = np.random.normal(mu_series, sigma_series)\n",
    "    print(\"reconstruction\", reconstruction.shape)\n",
    "    print(\"actual\", actual.shape)\n",
    "    plt.plot(actual[:, 0], label=\"real\");\n",
    "    plt.plot(reconstruction[:,0], label=\"reconstruction\");\n",
    "    plt.legend()\n",
    "    plt.show();\n",
    "    plt.plot(actual[:, 1], label=\"real\");\n",
    "    plt.plot(reconstruction[:,1], label=\"reconstruction\");\n",
    "    plt.legend()\n",
    "    plt.show();\n",
    "\n",
    "plot_reconstruction(mus, sigmas, x_test, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#train for 50 more and see improvement\n",
    "train_losses, validate_losses = train_model(50, sess, x_train, x_validate, vrnn, optimizer, batch_size)\n",
    "plt.plot(train_losses);\n",
    "plt.plot(validate_losses);\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mus, sigmas = reconstruct(sess, x_test, vrnn, batch_size)\n",
    "plot_reconstruction(mus, sigmas, x_test, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#train for 1000 more and see improvement\n",
    "train_losses, validate_losses = train_model(50, sess, x_train, x_validate, vrnn, optimizer, batch_size)\n",
    "plt.plot(train_losses);\n",
    "plt.plot(validate_losses);\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}