Tutorial - How to use BJ-RNN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### International Conference on Machine Learning (ICML), 2020.\n",
    "\n",
    "## Frequentist Uncertainty in Recurrent Neural Networks via Blockwise Influence Functions\n",
    "\n",
    "*Ahmed M. Alaa and Mihaela van der Schaar*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recurrent neural networks (RNNs) are instrumental in modelling sequential and time-series data. Yet, when using RNNs to inform decision-making,  predictions by themselves are not sufficient — we also need estimates of predictive uncertainty. The *blockwise Jackknife* method capitalizes on ideas from classical jackknife resampling to develop a frequentist approach for quantifying uncertainty in RNN predictions. This approach does not interfere with model training or compromise its accuracy, and applies to any RNN architecture. The blockwise Jackknife operates by estimating the variability of the (jackknife) sampling distribution of the RNN outputs, which is estimated by repeatedly deleting “blocks” of (temporally-correlated) training data, and collecting the predictions of the RNN re-trained on the remaining data. To avoid exhaustive re-training, we utilize *influence functions* to estimate the effect of removing training data blocks on the learned RNN parameters. This notebook provides a tutorial on how to use our method through a synthetic data example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 1: Import all necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline \n",
    "\n",
    "from models.sequential import *     # library containing the base RNN model implementation in torch\n",
    "from utils.make_data import *       # library for generating synthetic data  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Set all data and model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seq_len      = 10                                  # length of temporal sequences\n",
    "n_samples    = 1000                                # number of samples\n",
    "\n",
    "INPUT_SIZE   = 1                                   # feature dimensions\n",
    "EPOCH        = 10                                  # number of epochs\n",
    "N_STEPS      = 500                                 # number of steps \n",
    "BATCH_SIZE   = 100                                 # batch size \n",
    "MAX_STEPS    = seq_len \n",
    "mode         = \"RNN\"                               # RNN architecture: RNN, LSTM, GRU \n",
    "coverage     = 0.9                                 # Target coverage probability \n",
    "\n",
    "params       = dict({\"INPUT_SIZE\":INPUT_SIZE,      # RNN parameters\n",
    "                     \"EPOCH\":EPOCH, \n",
    "                     \"N_STEPS\":N_STEPS, \n",
    "                     \"BATCH_SIZE\":BATCH_SIZE, \n",
    "                     \"MAX_STEPS\":MAX_STEPS, \n",
    "                     \"mode\":mode})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3: Generate synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "noise_profile = [0.1 * k for k in range(seq_len)]       # variance of label noise within a sequence = [0.1, 0.2, ... , 0.9]\n",
    "\n",
    "X, Y          = create_autoregressive_data(n_samples=n_samples, noise_profile=noise_profile,\n",
    "                                           seq_len=seq_len, mode=\"time-dependent\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4: Instantiate and train a vanilla RNN model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RNN_model      = RNN(**params)\n",
    "    \n",
    "RNN_model.fit(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 5: Post-hoc uncertainty computation\n",
    "\n",
    "Our method computes uncertainty estimates in a post-hoc fashion. The function *RNN_uncertainty_wrapper* takes the trained RNN model as an input, and computes the influence functions per training point in order to evaluate confidence intervals on a new testing point. The object returned by *RNN_uncertainty_wrapper* is a super class containing the RNN object with extra attributes needed to computer confidence intervals on test data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RNN_model_     = RNN_uncertainty_wrapper(RNN_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 6: Generate a test sequence and compute confidence intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_test, Y_test = create_autoregressive_data(n_samples=1, noise_profile=noise_profile, seq_len=seq_len, mode=\"time-dependent\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Y_predicted    = RNN_model_.predict(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The predict function retrieves three outputs: the mean prediction (RNN output), and a lower and upper limits on the confidence intervals associated by this prediction over time. In what follows, we plot these predictions and show how they compare with the original test labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Prediction')"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEJCAYAAAB2T0usAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2de5yM5fvHP9fu2mVZ52NOi1CUwuqElEOIfElEmyjCr1QqldK5dNI3QiopOigJfYtSKJEOQgk5lKxT2HU+27W79++Pz0wzO4fdmd2ZeWZnrvfr9bx255lnnrln2Pu67+vwucQYA0VRFCW6ibF6AIqiKIr1qDFQFEVR1BgoiqIoagwURVEUqDFQFEVRAMRZPYDCUrlyZZOcnGz1MBRFUYoVa9asOWCMqeJ6vtgag+TkZKxevdrqYSiKohQrRGSHp/PqJlIURVHUGCiKoihqDBRFURSoMVAURVGgxkBRFEWBGgNFURQFagwURVEUFOM6A0VRlIjl11+BjAzg4EEe558PdOoU1LdUY6AoihJosrOBQ4cck/nBg0B8PNC1K59/7DFgw4a8z7doAcycCZw9C1x3HbBnj+N+3bsDF1wA1KgRtCGrMVAURcmP06eBAweAY8eApk15bsEC4LffHBP5gQNAmTLAhx9yMu/cGVixIu99GjUCKlcGzpwBliwB9u4FSpfmUa0aJ/pFiwBjgMGDgbg4oFw5HiJ8HzUGiqIoFjBvHtC/P5CVxUn7jz+AzEzgtdeAr74CEhOBpCQeVaoAc+Zw4m7ZEjj3XBqIpCT+rFgR2L2bk/zjjwOxsfw9xkPotnr1vI8PHAj6R5Xi2vYyJSXFqDaRoihBY+9eumYqVgSuuoqTeqtWnMSzsoCSJXnExfGIjaUhCAYHDgDnncfxFBERWWOMSXE9rzsDRVEUT9x3H3DyJPDww0CK29wZcagxUBRF8cTDD9NHf/HFVo8kJGidgaIoijPHjjEbaM8eoF07uoCiADUGiqIodnJzgf/8h6mchw8zkydKUGOgKIpiZ9Ik4LvvgNq1gapVrR5NSAkbYyAi5UVkjohsFpFNInK51WNSFCWK2LwZGD2aGUM9ejA7KIoIJ2fYqwC+MsbcICLxABKtHpCiKFFCdjYwcCCQkADceivTSKOMsNgZiEhZAFcCeBsAjDFZxpgj1o5KiTqys4FLLwXGj7d6JEqo2b4d2LULGDAAaNzY6tFYQrjsDOoD2A9guohcBGANgHuMMSedLxKRoQCGAkCdOnVCPkglwlm2DPjlF4qEtWlDd4ESHdSpw0VAmTKeK4KjgHD51HEAWgB43RjTHMBJAKNdLzLGTDXGpBhjUqpUqRLqMSqRzuzZQKlSwGWXMZNEiXwyM4GJE4E1a/i4dGlrx2Mh4WIMdgPYbYxZaXs8BzQOihIasrOpQ9OyJXVj9u8H/vnH6lEpweaJJ4B77qHwXJQvMMPCGBhj9gHYJSJ2Z10HABstHJISbcTFAf/7H3PMS5QAcnKAjh2BpUutHpkSLH78ERg3DujQgW7BKHUP2QmnT38XgJkisg7AxQCes3g8SrRRsiRQvz5/r1SJksG33AKcOGHtuJTAc/Ik/22rVePPUqWsHpHlhI0xMMastcUDmhljehpj1GmrhIazZ5lOuGiRo+K0dGlq0+zeDYwYYe34lMAzejSwbRtw220sMFPCxxgoimV8+y0wYwZ3As6FRs2bA716Ae++S5+yEjn078800rZtgyc7XcxQY6Aos2dzJ+BJpnjYMKBmTeDJJ9mBSine5OTw3zE7m60lExKsHlHYEC51BopiDVlZwKefMovIUzZJQgLw4otsVbh9O1CvXsiHqASQ226jW7BrV3UPuaA7AyW6+eYb1hRccol3LZqaNTlxLF/u3tdWKT7873/Ae+9xAVCtmtWjCTvUGCjRzcmTlB8oqJNVfDzw5pt0LezbF5qxKYEjIwMYOpR9ifv25b+nkgc1Bkp007Mn8PTTQOXKBV97++00HgMGaPygOGEMMHw4cOQIMHhw1ElT+4oaAyV6OXCAq/ycHN/kiuvVYwrqkiXcJSjFg7Q0/pvdcANjQ4pH1Bgo0cv99wOXX+6fy6BvX6BJE752+/agDU0JIHXrAq+/TmNQooTVowlb1Bgo0UlmJgOKDRoA5cv7/rrYWGDMGL5uy5bgjU8pOsYAn38ObNrEWoKKFa0eUVijxkCJThYtYuPzyy7zv6PVOedQ7vjgQeDQoeCMTyk6r79Oranp0zV7yAfUGCjRyezZQNmyQItCiuPGxrJ5+vXXc+WphBdbtwIPPMAq8nbtoq6FZWFQY6BEH2fOAJ99xmBipUqFv09CArBqFdCvH4PQSniQkwMMGkQV0sGDafSVAlFjoEQf8fHUG+ratWgrxipVqIW/bh3lKpTw4L//BX74Abj5ZuC886weTbFBjYESfcTEADVqsACpqHTqBLRuDbzwgqNblmItTZrQ0F97rbqH/ECNgRJdnD4N3Hcf8PPPgXEfiNA3Xbo0cOedRb9fuLBgATBhAnDqlNUj8Z8aNYCBA4GkJKtHUqxQY6BEF199xUygjIzArRrLlQOefZZNUvbuDcw9rWT8eMpu3HsvU2jfead4xESeeYb1H5s2RX0Ly8KgxkCJLmbPZl2BpyyiDRs4oT/1FOsQ/OGCC1ihvGIF8NdfgRmrVVx4Id1fjz7KHc/gwcCDD1o9qvxZtYr/br/+yqSAKG9hWRhUwlqJHk6dAubPZ9VxhQruz48fD+zaxaN2bcod+0OpUtwh7N4NbNzIibS4sG0bv5uePYH9+4G77wbKlAHat+duqlYt4PvvuZuKjy9Y2C+UnD5NI16pkrawLAJhYz5FZLuIrBeRtSKy2urxKBHIwoUUmrv0UncX0fbtnBDtvP8+J3V/6d0b2LmzeLXKXLoUaNUKePxxxgoqVaIhABgT6dqVO5/9+5k91aoV0KdP3u/LSh55BNi8mbpRyclWj6bYEjbGwMbVxpiLjTFhtOxQIoaMDK74mzd3f+7bb93PLVvm/3u0aMGq1xkzgC++8P/1ocQYYMoUuoSSkoAnnmAqpietJrucw+OP0+B9/jmlv+++m5XYVrF7Nz9Dly7AVVdpC8sioG4iJXoYMoQTmqsWkTFscuPMiy+y4U1hGD4c+OUXupm2bPFP+yiUjBwJTJzIlf4dd/i2qk5K4q7nhhuAt94CXnsNKFkSeP55a9I4a9akgmxsLMehFJpwMgYGwCIRMQDeNMZMdb1ARIYCGAoAderUCfHwlGLNyZPUs8/NdZ+0Nm8G9uzh76VLA3PmFG1iKVmSwdfnngNWrwY6diz8vYJJo0bcxdxyi/8ibtWrA489Btx4I3cSX3zh+H4HDAiNYVi/nsapRAnqRSlFIpzcRK2NMS0AdAVwp4hc6XqBMWaqMSbFGJNSRVPHFH8YNIjpkp6Ci867grZtA7PCbNKEK+eDB4ETJ4p+v0Cxdi3wySeOHcuwYUVT82zUiDuKuDgKw916K3DRRQw6B7MB0BdfAM2aAZMmMY1U3UNFJmyMgTFmj+1nBoBPARRyj64oLpw4wcmjShX3QrOcHAZQbdz888+Iad8eyf36YeaSJTx5+jTjDf5SqhRXyrffDqSnF+EDBIhPPmG19H33AT/9RCXPQGXeJCZyJ/TAAww0d+0KXH01jU+gOXiQLr/69enK0xaWASEsjIGIlBaRJPvvAK4BsMHaUSkRwxdfcEK/9FL3/PPjx4HzzkNOTAwyAMw6cgTGGOxIT8fQceOwfPJkVrM++2zhVrpnztDtdMst1rXKzM1l4LdvX67ix4xhw5e4AHuJRSgB8f77jJesWcPvPtBVzHfeSYMweLBKUweQsDAGAKoBWCEivwP4BcAXxpivLB6TEinMns10yYsvdn+ufHlg7Fi0qFQJPQE419lWycrC5XPncqW7fj3w9df+v3e9enRRLVpEt1Goyclh9s8zz7Bm4PHHmTEUTLdKfDzjBjNnUv9p/nxg9GgGrIva/+Hjj3n06sXAtxIwwsIYGGO2GWMush1NjTFjrR6TEiEcPw58+SUnjnyyetYfOICfXM7tAPCK84k33+T9/KVfP+D88+me2bHD/9cXhdhYBlcHDOBkHMqVdNmyfL9KlWhMJ06kcXzxRe6YCsPZs0zf7dNHW1gGmLAwBooSNOLjgVde4ao4H4mCOlWrejw/o2pVwP7ckSPAtGn+j8HeKjM72/+q5sLy7beMC/z+O3DFFZRztkq4LS6OsYQpU+jnHz2amkcLF/p/r5QUGtXKlQM/zihHjYES2SQkUGunUSP355x8+GOHDEFiQkKepxMTEvDo7bfnrSaeP5+pqP5SsyYrZbt1C26rTGOY+3/NNcBddwF//ME00HAIsp53HiU/nnuO/y5btgAHDgBZWQW/9v33qaL6668aJwgSagyUyOXYMWDsWEc+ujOnTwOpqXRdbNiA1A4dMHXUKNStVg0igrrVqmHqqFFI7dgRaNOGwWeAk+2ECYVT8bzySqBhQ+DHH4HDh4v++VzJymKq6IgRdKXcdRddROGm6X/55cDbb9Nl9PXXrFVo39575lFaGovipk9n9lOgA98KADUGSiQzfz6LvzIy3F1EK1ZQbvrTT+lGEkFqx47YPmsWcr/9FttnzaIhABhsvesuh496yxZq+BSGpCR24mrfPrCy0EePAh06MEjdsyddMXXrBu7+gUaE0t/nnMOYwqpVNGCpqXnjKrm5rF0whtlDngQGlYCgxkCJXGbPZm1Bs2buzzlrEbVvX/C9atbkRGVn2rTCr+4vu4yr4KeeKtzrPVG6NI877gCGDi1aIVkoEWE844MPgB49WAvRqBEfA8Crr1Ij6qabgKZNrR1rhKPGQIlMjh5lFWyrVlyBuj63apXjcYcOvt2zf3+H7MGJE8wuKgzXXMOg7vPP0wdeFD79lCvpX3/lpNqrV/GUcC5XjtlOM2bwuzl+nPpOjzzCf8Prrgs/d1eEoc43JTL5/HP60D0Vmi1b5nDRNGnCNom+EB9PCeeHHuIqvLDN1u2tMgcNooFZv97/AG9uLlVGn30W6N6du5aaNYu/LMM55/BzZWYCf//NzKFatSzJhJq5ZAnGTJuGnRkZqFO1KsYOGeJwHfpKTg7TYbOy6GZ0NdRbtzKhICvL/bC/LiuL30Fh/7/5iBoDJTL5+29OjoFwETlzySWUbW7XrmiumPLlgVGjgHHjgB9+oHSDrxw/zrqBzz7jrubmmyMvwyYhgUa6QgVLMqFmLlmCoS+/jPMzM3EHgNLp6fj2xRcBIK9BmDmTTX+8TebZ2Y5rb72VlejOTJsGrFxZ8IC6dCn6hyoANQZKZPLww6x+ddUi2r8fWLeOv8fE+DcJ2+nVq+jjA5ildP75jD2cOeObQN7OnZR82LyZMhm9e0d243eLZKnHTJuGppmZ+A5Aou3c29nZGDNtWl5jkJHBhAJf8JRC66uhO3vWt+uKgBoDJfLIznY0vHd1ES1d6qgvaN7c+kBrpUrAvn2Ug37yyYJbZSYlOYq4rr46POoHIpCY9HTMh8MQAEBJADtdBQt9+f7j43l4ink0bMiFgP2a+Hi6k5wfx8c7Ch+DiBoDJfLo1Yt/YHfd5f6cXYkU8D1wXBB//slq34EDC/f6gweBl1+m7/jtt92fNwaYNYtKoOvX0xBUr64B1WBx/DgWxcaimi2udBDA0wD+gIdK9d69+f/I22QeF5d/HGfAAN/GdOBAYT6JX6gxUCKLw4eZRdSli7v7ZNcu4K+/+HuJEuxdUBRyc4HJk+m7z81lMLow4mlNmzKt8p13OLlce63juawsqnROm0Yp7E6dGAtRgkNWFvD44zjXZggyAfwHwA9gRfrUIUPyXl+9Oo8IQFNLlcjif/+jm+iyy9xXZDa5agDMMrI3fS8sMTFMMc3N5eOJE32TVvDE//0fA6a33UYNJICurg4daAh69eLOQJs6BQ9juENzqoQeVa4cfnStSI9QdGegRBazZ3NSveAC9+eaNGE3rn/+8T5pr1pFMbVKlXx7v+HDKS9x8iSbs3/8se9bf2dKlqSY3V13sWhszBjm1u/fz0Ky7t2LZ/1AceLdd4HFix2Pb78dk266CZOsG1FI0Z2BEjkcPMiYQEqKexaRMzVrUhfHlf37ubr/80/f37NiRcok2PngA0c/ZX9p2pQ7hGbNmD2SkECjUFwLyYoTv/xCY2CnWzfWgEQRagyUyKFECWoRtW1buOKrJUu4uk9K4krfV3r0YFYIwB3HpEmF72rWpw93MFu3Uq6idWsVZgsFLVtyJwYw7jNyZPEv4PMTNQZK5FC2LJVBXSs1fZmYjaGCZqNGNAZPP+27kFxsLHDvvY7J4+efWUhWWCpWdIi4RdmEZBn2f8NRo1gBHYUGOKyMgYjEishvIlJISUglajl4kIHWbdvcA8Pr1rH6Mz8Xzl9/UeOna1emBP7yC1VPfeX88+nXtzN5MmWyC4u6hUKPCN1DBdV6RChhZQwA3ANgk9WDUIoh8+Yx9TIjw301/c03wPbtzOH/+GPPr1+8mKvBwYOpr3/xxTQe/rRnHDLEIYqXnu5Q3lTCj9OnKfJnzwRTwscYiEgtAN0AFKKvoBL1zJ5NMS9XmePsbArT2fFWaNasGXP8GzViyuj48dxt+DOhly3L5jJ25s6lQqoSXuTkUOBv4kT+LGw6cIQRNsYAwAQADwLwaqpFZKiIrBaR1fv37w/dyJTwZv9+is+lpLgXmq1ezY5nAEv6PaWcAjx/993M4AGAq65i4dq8ef61qezcmfdq2pTtJ13lsxVrsbcF/fFHPl66lDEeJTzqDESkO4AMY8waEbnK23XGmKkApgJASkpKIdM1lIhj3jxu9y+91LOLyM7VV7trFQEM9pYrRx19Z/77XzZb8WflGBPD1WZSkuf3Uqxl7ly6h+z068ekAyVsdgatAfQQke0AZgFoLyLqcFV849dfPbuIzpxhe0s7nlxEJ08yc2jBAqBy5bzPNWnC3gWxsf6pRpYrp4YgHPn+e2DKFMfjdu0YZ1IAhIkxMMY8bIypZYxJBtAPwLfGmJstHpZSXJgwgYqfrllEP/7oCADXqUNJa1eWL+fKPzXVs/BbyZLAxo3sO1AUipJZpBSdTZuAsWMdacZNm1LmXI32v+g3oRRvjGHmTunS7i4i1yY2nnL2Fy2i0Fi3bt7fIyuL2Ua//+7/+HJy6Grq25eFZEro2buXldyZmXx8zjl05dnjQwqAMDQGxpjvjDHdC75SUcCq0dGj3XcF9h66djy5iNLTKUrWoQM7anlj9GhqFb35pv+VxZMm0TVx4gQzlIpjKqM/1djhxvHj/Pc7fJiPy5YFXniBneaUPISdMVAUn9m3D/jyS/rzXQuFli93+PkbN2ZMwZUNG+gmuPXW/Ct9k5IoDbFpE7NP/OH66x3VrBs3AgsX+vf6UJGdDaSluZ9fsYIutMWLCy+xYSUvvcTucADlSp59Fqhd29oxhSlqDJTiy9y5nKA8yVUDju5Q3voct29PFdPLLiv4vYYOBRo0YOGarzIVAGMV/fo5Hk+dGn61BytXstju3nu5g7Fz8iTw6qsc73PPAQ8+SMXX4sTgwY7+0KNHAxdeaO14wphCGQMRiXE+Aj0oRfGJ2bOBunUpBeFKt27ARx+xsKhTJ/fnjeFk16CBb/IDJUrQTTRqlKNuwVdSUx0T0rFjwFtv+ff6YJGWxgl+9Giuno8ezVtkt3dv3gDr6tXst/DRR3kbvYczycl00z30kPdFgQLAD2MgIi1E5CcROQngrO3Itv1UlNCydy9TBVNSvE/mMTFcCXqKB0yZAjz+uOcMI2906AAMGsTMIH98/yVLsqDNzhdf0GVkFUeOMH4xZAj7N9hJTMzbx+Hcc4EZM4AbbnAYhaws7m6GDaPbrDhQsSILCJV88WdV/y6ApQBSANS3HfVsPxUltIiwsUzr1v4re2ZnU646Pp6NcPyhXDkK2HnqVZwfV1yRt6ht/Hj/3E2B4OxZ7qZuvhn4/HOHQYuJocje++9TQtuZUqXYdnPKlLyGc9s2np80CTh1KnSfoSDWrgW2bLF6FMUSf4xBXQBjjDGbjDE7nI9gDU5RvFK9Opu+2PsI2PElyLlqFVfHvXvTIPhDTAx3BvPm+d+k/K67HOmMW7eyd3IoMIaB4EGDGCNxzg5q0YIr/fvv5wraG40bA2+8QQNs/wzG8HsYNMjanY6dtDT2sxg50iE3ofiMP8bgUwDXBGsgiuIz6ensdbxvn7uL6KOP6Nf/8su8wVBnFi1ihlBhO1mNG8fdxdSp/r2uenWuyu288w7F8ILNmTPciTjLd9euzSKsl19m3MQXYmOp6Dp9OhvA2Dl1yvqm8IcOsYjs5El+3qL0o45S/DEGJQF8KiKLROQ95yNYg1MUj8yaxV3B/v15XUTGcKJfs4YT9sqV7q89cYJaRG3a+O8istOwIWUMvvkG+Ptv/17bt68jtTEpiZ8h2JQqxcAvwHqMO++km+uKKwrXPKdGDeDFF7kKL1+e8YP8dhXB5vRp4JFHuEgA+Hmffdb/XV+U449Q3UbboSjWMns2m9Y3bpz3/LZtbFADMGjrKjwHcPLr14+tKj3JT/jKk08C773H3cGLL/r+uvh4ujHWrgVuuonjDCSZmVThbNcu7/kuXega6949MEqqIgyoX3KJ5wD+Z59RFtxTr+lAYpejtscJYmLYqcyfxAAFgB/GwBjzVDAHoig+sWsX/cE33ug+CTkrlF5xheduYaVLc2L01tfAV6pWpbvkwAGuTP3pTNaiBY9AYgwL4qZO5Qp54sS8OfWxsUxxDTSukuEAJ+aJEzkx9+9P11iwVulTpuSND9xzD9VrFb/xq0ZARK4WkXdE5GvbT03cVULL3Ln86foHb0xeLSJPk31GBtM6ExICI0fQpw/lMOxSB1axcSOD088843CVTJ5sjfSFMSxUy81lXOX991n4tXZt4N9rzhwGsO3Yd3xKofCnzmAIgI8B7AMwD8BeAB+KiGrAKqHj228Z8GzUKO/5P/5wTIRJSXkDnHa+/poB04oVA9dovmRJNkv5+uui3efHHxkQ94eMDAaB77yTn99OuXLAtddaIx8hAjzwQF458d27Wd380kv+F+x5Y8UKlaMOMP7EDB4E0MkY8690o4h8DGAugDApqVQing8+4GozPxfRlVeyYtgZe3D5vPPo5w4UFSowk2X6dDbP8dcdsm8fc/V//JFB7WeeKfg1p08za2r2bIcSJ0ANpN696ZZxFe4LJfXq0U20YAHdVvZU1oULgZ9+ovHq0KHwBnnPHsYJVI46oPjz7VWCewB5CwAL0wiUqGP//rxVsgCDiN9953jsyUW0eTNXqN26+SY/4SuxseyIlp7Oydlf9u93+LxXrOBk6Q1jgK++AgYMoEF0NgRt27JaePhwaw2BnZgYumxmzMjbSezIEe5mHnqIVeSFoUYNfgeAylEHEH+MwQoAr4hIIgCISGkA4wBodYcSGnr1oquhbNm859es4SQD0FA0a+b+2kWLuFu49dbAj6trV054s2dTMtkfLryQr7czaVLeSd4ZEUpwONcmNGzIGoKnnwZq1vR/7MGmcmUqvj77LFCliuP8qlWFbxgkwmD4Y4+pHHUA8ccYDAfQDMBREUkHcATARQCGBWNgipKH7dtZaHbsGDV0nHEOHF99teeU0R07gJYt3SuWA4EI8MorNATvFaLsZtgwh4HbuxeYOdP7tcOH8/NVrMjV9RtvABdfXLhxh5LWrblLuP56fl+xscCIEUW7Z/v2KkcdQPxJLd0LoJ2I1AZQA8AeY8zuoI1MUZz55BP+vPxy9+dKleJx+rT3lNFHH+WuIVgpji1bspDLGGbRxPkRjitXjqJxr7zCx7NmcfJcuZJFas61CLVrc5V90UX+pbOGA4mJzHrq2BH46y/WijiTk8OqYdfPZQyzwK65RgvJgoiYfDIORESM7YL8pKqNMUXKYRORkgCWA0gADdQcY8wT+b0mJSXFrF69uihvqxQnWrWixPKECe47A4CuldWrPVfVZmfTN9+hg6PHQbD44w9g3Tr/q5tzcxlY3bw57/lBg4CBAwM2vLBmzhwa/XvuyVswOGMG8O67NIDPPOO5tiHSOXCAyQ8XXFDkW4nIGmNMiuv5gtxEzl047HLVzkegJKwzAbQ3xlwE4GIAXUTEh44jSlSQlsaJPiXFsyEAGED0pGB6/DgzbH7+2T3wHAyys7kL8Vc5MyaG6ZeuGTGzZjFbKdJJT+fOKiOD/YqffJKxka++oiEA2IN6zhxLhxnJFGQMnJKF/5Wrdj4CImFtiF1VrITtKIY99pSgkJvLCd2Ti6ggvvuOcYbLLy+a/ISv1K7NVdzrr/v/2kaNgJ49HY+rVqWaaDQESPfsyZsRtGwZd0Qvv+w4l5LiyCKKImYuWYLk//s/xDRrhuTkZMzML6ZUBPJ1bBpjdjk97GOMedn1GhG5D8ArRR2IiMQCWAPgXACvGWM8qIwpUUmDBhRaK0xe+uLFzLLp3Dnw4/JExYpc2T70EFNGPekj5ccdd7BVZlwc3VqB1i4KV5o35w7gjTe4GwDySm3Xr8/dgj+xmAhg5pIlGPryyzhlyzDbsWMHhg4dCgBIDbC8SL4xgzwXihwzxpT1cP6QMSZgtQYiUh6Uy77LGLPB5bmhAIYCQJ06dVru2KGtFCKeffuoDLpjB3PKnXniCa4mO3RgANd1otizhymIt9xCv3Ogqo4L4swZCqXFxVGmWouh/OPXXxlMt/dbrlyZVd7BjveEIcn9+mGHvbLeibp162L79u2FumdhYwYQkfY2DaJYmzZRe6djCAA/E6vzxxhzBMB3ANz61BljphpjUowxKVWcc5aVyGXGDFbmuur/HDrEIq3Fi9nD11NfgMWL+bOwu4rCUrIk8PzzNGD2VW4xYuaSJUju1w8x7dsjuV8/zFyyJLQDaNGC8YNhwyirMX58VBoCANiZkeH5/M6dAX8vX/Zc9v5+JQG843TeAEgHcFdRByEiVQCcNcYcEZFSADoC8EMXWIlYZs9mFoVrGuKyZQ4htgsvdDScd6ZtW/70pFMUbFJT2TshMZGpkT4ao5lLlmDMtGnYmZGBOlWrYuyQIUjt2DHIg837/nncEunpGGrz24dyHEhIoAHgQ4MAACAASURBVPBclFOnalWPO4M6deoE/L0K3BkYY+oZY+oBmGn/3XbUN8Zcboz5PADjqAFgqYisA7AKwGJjzIIA3Fcpzvz1F/DbbwwcuuaeO2sReastqFyZLiJvGUjBJCYG+L//Y5Gbj9lA9ol4R3o6jDH/TsShXJmPmTbtX0Ng51RmJsZMmxayMSgOxg4ZgkQXqY3ExESMHTs24O/ljzPzFVvB2b+ISG0RuaiogzDGrDPGNDfGNDPGXGCMebqo91QiAHuh2WUuWcb79jlUOmNi3Bu5AJSfWL3a+iYn+/axmY0PMtfhMBF7dUt4Oa8El9SOHTF11CjUrVwZIoK6deti6tSpAQ8eA/4Zgw/AlE9n4gG8H7jhKIoT8+cD55/v7iJylp9ISXFPvTx7lnr+33zj2X0USho2ZCD0nXcKvDQcJuI6Xnzz3s4rwSe1Y0dsf/115K5bh+3btwfFEAD+GYM6xphtzieMMX8DSA7oiBTFzrx51Kj310W0ciWLzW64wXr5ggsvpKT0V1+xS1s+hMNE7NEtkZCAsUOGhGwMijX4Ywx2i0ieXn22x3sCOyRFsXH0qHs6aVoaex0DnOjbtHF/3aJF1Pu58cbgj9EXnnuOBW9vvpnvZeEwEf/rlqhWjW6JatUwddSo0AaPFUvwp4JjPIDPROQlAH8DaABgFIDARzIUpV8/VvO2d+ms6uwiuvxy9+DwsWPsCdClC1C9evDH6Qu1alFv56WXgD//dO/SZsM+4VqZTWQfh07+0Yc/qqVvicgRAIMB1AawC8D9xhgVC1ECy6ZNwMcfU47A2UXkS5/jnTvZvGbAgPAq9nrkEaBx4wIrinUiVqzCr9puY8wnAD4J0lgUhXzyCfPyXZvenzhBt9G+fTQSrs8DVHWcNAno1Ck0Y/WVcuWoQLpwIXcvrg16FMVi8jUGIjLAGPO+7ffbvF1njCk4VUJRfGX2bKBJE/bSdSYpid2xDh1i7MA1OHzmDDXxK1UKT3G3mBgqb77zDuMHoRDOUxQfKWhn0B+O1FFvcoEGeSuTFaXw/PEHj0GDvLtUKlbk4cqsWWyCsnRpaOUn/KFhQ2otffops52U8OannyiJYm9YdOut3JFu20aDHhfHIzaWP/v0YUwoLQ1YsMD9+WuuYSxr9262a3V+Li6OGltJSey/sWuX47wIK/GDSEGqpdc6/X51UEeiKAA7XbVt615oVhDGMIuoWjX3uoRw4oYb+Af/4YdAt27Fr1tZtGAMeye8/jobFdWqxfqVM2fYcyE9nTvUnBzHkZ1NQ1G2LPtZfPWV43xODu9buzYn9p9/pvieK889ByQnA0uW5K1Nueoqh7xKkCio05lPEbiidjorDNrpLEIxhiuqEiX8k29evx64+25KR7/wQvDGFwi+/x648kqgf3/AJkeshBnGAI89xsXJ++8z+O/6vOuRm+v9cW4ujQJAY3D6NFOns7JoZLKyeNSuzf/3+/Zx92F/rkoVplFXqFDkj+ZNtbQgN1E2fGsyo85PpejYV1snTuRtG5mVRZ2fVq2YQXTuue5uoMWLKW52662hHXNhaNuWu4LPPgNuugkoU8bqESl2TpwATp3iYmTECP5/86SQLFJ0V2StWt6fq1075AKLBa38nbub3QVgGSgtfb7t51IAI4I5QCWKeO01oFkz/jE688svXCV9/DF7GLiSlcU4QatW4e0icubVV9nMxbmBS7hgDCdFC7FERvuff9iHevRoJid07+7ZEEQoBcUM/u0eY+tolmLrNwAAf4rIagCrARSix5+iOGGMI4vIdcXkLD/Rvr37iiw2li6i887jiq440KABO7AtWECXgdWxg5wcYMMG9ohYsYI7swcfZKC+RImQBuQtkdFev569q3Nz6bfv1Ml6KZMQ409VTjkArlrAibbzilI0li9n0K1Vq7y9cE+dYvtIO54KzWJjWV/Qxa0fUnhTsiQwdy7w1FPWjuPDD9ljeuRIuq5q1OB3edFF1Hm64w5H17EQEHL11kWL2Gs6MRH44APuDqLMEAD+FZ29C2CJiEwAq49rA7jbdl5RCo8xDPxWqeIuP/HDD3QDAXQBudYeHD0KfPQRjUSlSqEZbyBp3Jgpsb/9xj7Aweb4cWay/PgjcN99jsDlBRcAXbsyNbJuXYdBbtyYaZBDhgB33cVrgrxLCKl6a3Y23Y8NG9Jt17Jl+KYlBxl/jMGDALYCuBHAOQD2ApgM4K0gjEuJJjZs4GQ4cGDewDHg7iJyZelS/jH36RNe8hO+MmoUMGUKG8G/8UZwJqKjR/k9rVgBrF1Ll1D58sDmzfxOX3yRhtRTs/m+fYFLLuH3O24cjchDDzEXPkh47e4VSPXWzEy6hM6eZYygSxcGbaMYn/96jDG5xpg3jDEdjDHnG2Pa2x7nBHOAShRwwQVclXXunHcyPHqUDWrseDIGixfzj7i46vmUKQM8/TQF7AIZJN2xgzpNAF08r74K7NkD9OxJV8iff9I10rIlazM8GQI7ycncTTz2GH/OmcPdXJAIunrroUN0iT35JD9H375RbwgAP3YGIiIAhgDoB6CKMaaZiFwJoLoxZnawBqhEOHv3ckWfk+OuMrpsmaNYp2lT913Drl3Axo1seF+ctX6GDGHTd7u7qzA7nNxcrvTtAeBdu1jLMHw43T4LF1LltWzZwu0+YmNptK6/nnGctDRmHNWvH/CgfVDVW//+m6KBR48CDz/MimCrg/dhgj9uoqcBdAIwAcAbtnO7QWnrIhkDWzvN9wBUB5ALYKox5tWi3FMpBhw7xlTS9u2pMupKQS6ixYs5sQ0eXLz9vHFxjBukpXGS8rWwyBjH5x4xgmqvsbHMyOrZk/0cLroosHUMF1/M961QgVXilSsz3TfADdqDot7600/AM88wcP/mm0Bqav47oijDn29iEIDmxpgDImJPJU0DaxCKSjYoh/2riCQBWCMii40xGwNwbyVcGTcOOHCAbiLXCSsjA1i3jr/HxLAc35UzZ4AWLUITeA02zZszOP7FF5Tg9pbNcvo0M3xWrGD21SuvMAjarh0zgvr2ZTDUn+ptfxFhO9I33uDO4/bbmXHUo0f4GuUzZ4CXX6ZL7O23gauvDt+xWoQ/xiAWgL0Sxe4wLON0rtAYY/aCAWkYY46LyCYANQGoMYhU9uwB/vtflth7qrQ8cICTY1oaJ3xPwnSpqZSGjpRtfsmSNJANGzK90Zn165kCumYNg55JSez/XKECjekNN4S+xqJ/f1ZT33gjMGECV95PPZU3NdhqcnI46efkMPDdo0fQBd+KK/4Yg4UAXhGRe4F/YwjPAJgfyAGJSDKA5gBWenhuKIChAFAnwNtSJcQ88QRXtL17e17FNmlCoa60NEdqqTNHjnCVnOImsVJ8KVmSUhvz59Pfv2kTDWWNGtSq+ftvZr307Alcey1Tca2Wwa5Vi1pLzz9Pt93hw0DVquGR2XXiBI1TcjI75w0eXDzTj0NEvkJ1eS4UKQv69bsAKAHgDIBFAG4xxhwPyGBEyoCSF2ONMfPyu1aF6ooxJ05w9duiBbM6/F3RZmUxkNmlC4OuxaXq2Bf+/psrV7uoWWoqJ7LkZDb2qVAhfN0bmZl07S1fzqyj22+3rnhrzx4GinfvZn3EM8+oBpSNwgrV2V8sACoDuAFARQB1AewyxuwL4ABLAJgLYGZBhkAp5pQpw/TELVv+nchnLlnie/bIjz9S0+eqqyLLEACUqXj3XcYFevRgjn+ZMuFrAJxJSOBO5ssv+e/7yy9M33QtFAw269czDTY7mzGV4cOjsqLYX/zZGZwEkBQMuWqbsXkXwCFjzEhfXqM7g2LKnj0s+1+4kO6E2Fg3LRqAeeVTR43ybBDGjKELZePG8Gl6r+Rl3jymzJ44QZnu3r1DY9COHmUso1w5Brivuy48XFZhhLedgT/f0m8AGgVuSHloDXZSay8ia23HtQW9SClmGONw75Qo8a+/21mLpgWAXwAMy8zExDffdL/HkSNcNV95JY2JEp5cfz2N9RVXUI3244+D+372RW18PDOb5s/nzkoNgc/4E0D+DsBXIjID1Cb6d0tR1B7IxpgVAIrBPlgpEnPnciK//fY8gTxnzZmbALSyHe8eOOB+j6VLmRkycKD+oYc71avz3+v115nSmZ7OtNlA++4zMymp0aoVK6pHjdIdYyHwxxi0BusK2rmc1x7ISsGcPcuKz3r13GQn7Fo0Agpf2fm2fHkMdL1Pp0786akITQk/RLhSP3uWC4E+fThhjxgRmFqIQ4foNtyyhdXQnTvTRaT4TYFLKxFJFJHnwHqC5QC6GGOudjr0r1IpmKlTga1b6Tt2Se+za9G0BWDvZLAfQOfhw93vExND5Uz9gy9elCjBlft//sPCuiFD+P+hKPz9NzvgpaUxVXnSJP1/UQR82WdPBnAdgE0AegN4OagjUiKTOXOACy+kr9+F1I4dMXXUKAxzKh473LIlburcOe+FCxbwPo2CFbpSgkpCAoO68+ezInj4cKYG5xYiJyUjgymjZ89SWmLMmMgpPrQIX4xBVwDXGGMetP3ePbhDUiKSt9/mKs6L9HFqu3a4ySlNtNFAFwdRbi4bk69Z4y5YpxQvundnNlj79owpHDlS8GtcKV+eu8w5c4BbblGNoQDgizEobZOLgDFmF7SzmeIPBw/y+P33/PPNV6+mcB3AYGPTpnmfX7eOq8EePXQFGAlUqQJ8/TVlu0uVoqT299/n/5rsbGYmrVvH/1Njx3pug6oUCl/MaZyIXA1Hto/rYxhjvg3G4JQIYPRougVefDF/NU5nLf+rr3bPFFq0iAFH1x2DUnwRodrpOedQKuK99ygpPXKku8G3S0usXk1306uvqrREgPHFGGQgb7bQQZfHBoFRLlUijY0bqS/UpQv/4L2xdSvw3XeOx659jjMz2dvg0kupza9EFnFxwFtvMfg7ebKjgvj88/m8s7TEyJEqLREkCjQGxpjkEIxDiURGj2a1cc+e3mUjcnOpeGkPIjZvTkkGZw4dojbPjTdGnvyEQuLjgYkTWTF8880MDj/8MDWs7r6bLqL//pdxJ5WWCApataMEh+XL6R7q3p35394Q4SRfpQpXiCNHuvuAa9QAHn2UOepKZNOpE4PLN97IVpQlS1LQ8IMPaCDUEAQNDcErwWH+fE7w3bvnL7MsQk38lBS6B1ylyU+eZJvFxETPPQ2UyKNiRWDmTCYU/PILa1QaNNBAcZBRY6AEB3u1cX6xAmdKlaJCpytffUU5g6+/VvmJaKNsWSDQrS8Vr0TtX9fMmTORnJyM70Twv9Kl8cOdd9I3rRSNrCxWhP7+Ow1BUVdzixbRXeCpG5qiKAEjKo3BzJkzMXToUOzZsQP7AbQ+dQqtp0xBbuXKVFn83/+sHmLxZepUNmdZvZqFQZ44e5bZIb/8kv+9duxg/vk113CVqChK0IhKYzBmzBicOnUKZ22PSwH4BMAH8fHA3r0sftm8GVi7llrsn39O37WSP8eOAU8/DTRunH+f2Tlz2C/3oYdYf+CNRYvoGho8OPBjVRQlD1EZM9i5cycAoDaAXuCX0AdgPnuVKsx3/u03GoP332cOdEIC0K6dI/XN26o3mhk3Dti/n1kfXmQnsG8fi4vseKtKNoaFaBdeyENRlKASlcagTp062LFjB6oD2Aogzxp21SoetWoBvXrRGKxfz3O//cam35dcwpXv779TW791a81/37OHeeBt2uTv3580iSJlAFNOe/f2fJ0Idxnx8So/oSghICrdRGPHjkViYiK2AqgK4HcAy0WQpwHo7t2cuAYOBDZsYGPyDz+k4FpaGmV4R46kdELlysyBf/dd6udEI0uXcjXfu7d3nfoVK9i/2M699+afdlq2LPXpFUUJOmFjDETkHRHJEJENwX6v1NRUTJ06FRXq1MEjAGqJoK0xkDZt6AYqXdpx8alT7Of62mtcrdarRyG1GjVYMj9qFJt1fPMNMGgQcO21jvjCX38VTp63ONK7N78jV4E5O6dPU2rATrduwAUXeL72zBnq0Gzfrq0tFSVEhJObaAbYO+G9Aq4LCKmpqUhNTeWDw4cZzJw+nVWwr79OF9C8eYAtvuDRnVG2LCe1bt3oLlq/nnK88+fzPv37s4CmWzcWX3Xq5N2XXpzZuJEZQqVLe3eXvfce2x4CjMncfrv3+/3wA7WK+vVTaWJFCRFhszMwxiwHYE2if4UKTIn8/Xc23EhIoO97zBjgpZcom9yiRd7XZGVxJzBlCjOQYmOBiy8GrrqK/VdjY3mvBg3YDNze4Wv6dL4+UnYMy5dzNzB9Ot1lnkhLAz75xPF42LD8O1ItXszvqlevwI5VURSvFKtll4gMBTAUYBA44DRpAowfT7ncL79k79bGjSmU5Vo89e23wK5dPObMYX3C9ddTaE2EBqZ3bx6ZmQxAr1xJHfbFi9mzdfx47hr69qXxKW4YAzzwADOwWrXy7P83hkJ0OTl8fOGF+ccBDh3id9Wrl3fjoihKwAmbnYEvGGOmGmNSjDEpVapUCd4blSnDSXzcOGbJDB/ORhoHDjiucQ6EGkPXxv33Myd+wQJHxgzg2Gncfz/jC1lZdJmULcuWfW3bAvfcQ2XG4sScOSwc69nTu+yEvXAMoLG49978ZSW++Ya7pkGDVH5CUUKIGGMKvipEiEgygAXGGC+RRQcpKSlm9erVQR8Tjhyhzs60aZzUP/yQk3huLifCuXNZbeuKPZ7wn/8w4OyNY8cYo/jqK15bXKqfs7LoHsrJodHMr9FIRgaDx+ecQ8OaH19+SdfTokVadawoQUBE1hhjUtzOqzHwkc2bGVCuX59B0l27gIsuoktoxw7g008ppua8IwC4un32WeDyy/O//5df0tikptI1Fe4Kjb/+yvjIsGE0er6Qk5N/KinATKwSJShBoShKwAl7YyAiHwG4CkBlAOkAnjDGvO3t+pAbAzsnTnB3MGwYJ+177nF0ZDpxAli4kIZh716eS0wEZs/Om67qjZwcVuguWQI0a+ZZ2z9cyM2lm6hUqcBlSKWl0VhceSWb2SiKEnC8GYOwCSAbY/pbPQafKFOGcYHTp5kLf8cdbNM4bBgDqX36MJD888/cSdSr524ItmxhAVaPHnyNndhYupQ2bABmzKB+z/TpvhmSULJ+PQ1AdrZnQ5CTA/zzj3tvgvzIyQEefBA491wG1BVFCSkaoSsMsbHcEaSlAXfeSR/38OEOF1FsLCUq7G36XPnkE3Zu6t+f/Vz/+MPxXFwcxdtuvpkr75YtWbwWLuzZw17E99/vvdnMvHnAbbcxOH76tG/3/fRTBui7d/dewawoStBQY1AU7A28N2wAnnuOweb0dK767XUErj7ygwcdzd9zcpiiOmIEZlx7LWYuXszzdqXOsWO5wr7sMgaaw4Enn+SOoG1bz5P2/v3czeTkALNmUbajINLSWOfRvDl3WoqihBw1BoGgUSO6ia67Dvj7b8pU/N//sTLXlfLlgccfR7qLC2XQ6dM4+8IL+HDRIsfJyy/nJDlwIDOWTp1iGqtVbNxIbaZOnbzLTrz2mmM3ULcuM6Ty4+xZGtJSpagFVaFCYMesKIpPqDEIJKVLA6NHc7dw4ABdSE89lVe8zhYgvTQzExcD+Nzp5YNyc1Fy/HhHgRZADaQePbjrGDOGMtr794fqE+Vl9GgGxHv29Cw7sXIlsGyZ4/HIkQWruR4/zu/kgQcKzrhSFCVoqDEINDExNAJpaaxc/vFHBkbPns1z2c6MDPwO9lOY7nT++jNnmIrqWoBWsSInzZ9/ZqbRypXB/iR5OXyY7rDu3Zle60pmJjBxouNx586U5yiIMmXoeiqoGE1RlKCif33BomxZ4NVX6Vp58UXGCnbtYtpobi7q2NQ4cwEMBvCm82u/+447ChcDgu7dec+cHPrsJ08OnduoQgVKdPfo4blWYOZMBpcBZhgVVFx28iRlKv7+mzuCcMuYUpQoQ41BsGnQALjxRsYT1q1jUHj4cLzZpQsSExIAAAbAcABTnCfZuDjPK+Xzz2fntSZN2FHMHnQOJps2sb/Dzp0U4XNl507go48cj2+/veBOcJMns51obKx/KaiKogQFNQahonRp7hCmTAEOH0bnd9/FH/Xq4YrKlSEiqFutGso99BANh10x1Vu1brlywCuvAI8+yljC7t3ulc+BIiuLu4FevfgZXIvg7EJ0drdWkyYFVyQvX075jeuvB265JXwL6xQlilBjEEpiYphllJYGjByJ5G3b8EOtWshdtAjbZ81CaqdOzEp68smCdfxjYljsVq4cC9Tq1GEfhUAzdSqwdSulJzyt9r/5hu1A7WMqyPd/6BDrL+rVoxHRmgJFCQvUGFhBmTKUr960CXj+eWYeHT3K50TcdwTGsCjrxAn3eyUkMIUzMZEr+EceyZuNVBSOHWMf4mbNKBHhidq1mVoLcKV/7rn53/ONN5h6Om4c+0wrihIWhI02kb9Ypk0UDHbupGulc2cGiZ0xhoqpH34InHcem+14koA4cwZ44QWmdnboQD0kbxXCvvLYY8xseuaZ/Pst5ORQaK9DBxql/Ni1i0HjMWMKTjtVFCXgeNMm0p1BOFC9Ot09EyYAa9bkfW7TJhoCgMqp993n2EU4U7Ik8MQTTGtdtoyTeFEMvTFsMtOmDZDi9v8mL7GxDJDnZwgOH2b8ITaWKbdqCBQlrNCdQbhw+DC7haWnM9Ombl3Hc59/TreSneRk+t29rfz//JPGoVEjxhK8NZ4piM2bKZZXr17hXm8nO5uZT4mJ1GSyq7wqihJydGcQ7lSowIYuJUoADz3ELCE7PXqwcM2edbN9O6t7vVUiN2pEX/6GDfT333qrf9lGO3YA27ZRndSTX3/hQhooX/s4z5xJw3LttY74gqIoYYUag3Cifn1OsiVKcDJ2pmtXdlyzZ+rs2kWDkJ7u+V4iXNG3a8dso0sv5STvC3ffTXG8rCx3d86hQ9QfGj8eGDGiYGmMTZuA996ju+nBBwtubqMoiiWoMQg32rShm+eCCzjxOrvxOnViUNc+oe7ZQ4Ngb6TjSmws3TOPP04Z7IsuYn5/fixfToPUoYPnArMpU1g9DFBXqFw57/c6fZoidBUq0PUVqCY4iqIEHDUG4UipUpRoeOcd5vk7c9VVlKqwr9j37WNvhd27vd/v6qvZZzkpia/11mPAGK7eq1RhVpPrKn7NGtYV2Bk5EoiP9/6+Bw8y0+ixx+iuUhQlbFFjEK6UKcMir1mzgAUL8j7XujWzhewG4eRJrtLzo25dGpZhwxib2Lgxb1wCAObOpQBez57uQeesLOoi2Wnfno138qNSJWZIDR2qVcaKEuaEjTEQkS4iskVEtorIaKvHYzki3Bm0bu055fSSS1iwVqEC6wt8ydApVYqZSDExQO/eVBX9/XfH87/9xrhF587uk/esWYxTAJSlyK8JzdGjNDx799LtVaqUTx9ZURTrCAtjICKxAF4D0BVAEwD9RaSJtaMKA+LjKTGRnMwaAtcAcMuWzNS58EL/7puUxB3C0aMMLL/7Ls/ffTf1jipVynv9P/8wJdTO4MHu19gxhmmvn3xCI1DYtFZFUUJKWBgDAJcA2GqM2WaMyQIwC0ABLbKiBHvKadmyTBV1rQvxtOrevt1zlzVnLr6Ylc316gGDBnGCX7OG8QJnjGGfArucdqNGTHX1xtdfA99/D6SmcvehKEqxIFyMQU0Au5we77ady4OIDBWR1SKyer9V3b6soH59Sjj07Jm3a5ondu1is/oHHmCdQH5UqsRWkz16UIL6r78Yq3Bm+XLgl1/4uwgroL2lh+7dy/uddx4VWvMLLiuKElaEizHwFF10K402xkw1xqQYY1KquK5gI51SpRgnWLmS4nGeKseNoY7QoUPsl/zggw5FUW/ExVFpdM4c977GxjBWYKdHD6BxY+/3mjCBr3n1VaBaNd8/m6IolhMuxmA3gNpOj2sB2GPRWMKXEiWY1790qXvKKcCV+8MPO5rKnznDvsX2lX1+JCa6y2aL0P9/ww1A5crAkCH532PgQKqmdujg2+dRFCVsCAttIhGJA/AngA4A/gGwCsBNxpg/vL0m4rSJfMUYYMAABo7vv99d5RSgCur991MaG6AReeqpojWcP33ae1bQsWPUQjp4kJIT+RWiKYpiKWGtTWSMyQYwAsDXADYBmJ2fIYhqCko5BShON2GCw1Vz9iyrkJcvL/z7ejMEWVksehs3jkJ7aggUpVgSFsYAAIwxXxpjGhljGhhjxlo9nrDGnnJaty6NgSfBuJo1aRDsqZ3Z2dwdfPttwfc/dMj3sbz1FrOXOnViv2dFUYolYWMMFD+pUIHZQg89xCweT+6+6tVpEGrbwjG5ucDYsUz/tDFzyRIk9+uHmPbtkdyvH+Z89hlVTp96yuFm8saaNQw8d+5M0TqtMlaUYosag+JMYiILzo4epVHIynK/pkoVKoza+yPk5gJr1wKgIRj68svYkZ4OYwx2pKfj9MSJjAF89x2D0d5iSidOMH20enWmkxbU4UxRlLBGjUFxJyaGfvpVq7jq9zR52zWC6tenptCoUQCAMdOm4VRm5r+XtQEwwNnlNGSI99X+rl10PT33XMF9jxVFCXviCr5ECXtSU1kw9tRTDjE6V8qXZ/5/qVL/Fo3tdCpgKwHgdefrr7ySUhXesAvfXXeduocUJQLQnUGk8MQTQP/+LBL74gvP15Qpk6d6uE7VqgCAawDcC+AC2/kTIowBeGL/fr7HwYNsnJOQELCPoCiKdagxiBRE2NHsiiuAFSs8xw9cGDt4MCbHxuJrAM7pW1s6dHDXKAIYb3jhBb5PpUqer1EUpViibqJIIj6emUJ79jCGUKOGo02mB1Kzsth8Bo7/CIerVkXL0V4UxOfNA379ld3TrrkmwINXFMVKdGcQaZQpQ2XRSpWYYeTawMaZjh2pd+REBee2ms5s384YQfPmjE24SlcoilKsUWMQqSQlAevWUZvIm8soIYHC+5XwYQAACClJREFUdtdcw99vvZW9l10xhu6hkiWZRmrXPlIUJWIIC22iwhC12kT+MHMmcPPNzAx68sn8s36Myf/5n36iy+mBB/J1PSmKEt6EtTaREiRSU5lltHy5Z5VTZ7wZgtOnubNo0AC48041BIoSoehfdqRjTzn96SdWKvvDqVPsgPb226w5cG18oyhKxKDGINKxp5wuXkwJCXv7Sl+YPBnYt49Vy3XqBG2IiqJYjxqDaCA+HmjYkDpGzz7LfgcFsWIFsHAh0KsXm9ZolbGiRDRqDKKJsmWZYfTgg/mnnB46BLz8MiUnXnmFWUSKokQ0agyiieRk9kE4fDj/lNM//2Qxmt0gKIoS8agxiDbatGGntC1bvKucNm0KTJ8O/Oc/oR+foiiWYLkxEJE+IvKHiOSKiFvuqxIE7Cmna9cC27Y5zu/eTTmLY8cYNC5RwroxKooSUiw3BgA2ALgeQBEa9Cp+88QT1C8qU4a1BDk57E0wcSJQrx5QsaLVI1QUJYRYLjBjjNkEAKLZKqFFhBpGZcpQrvrQIWDTJuCRR/LvY6AoSkRiuTHwBxEZCmAoANTRvPfAkJQE/P473UVt2jCw7EmoTlGUiCYkbiIRWSIiGzwcfkUojTFTjTEpxpiUKqqlHxiSkoAlS4BBg4DXXuNjRVGijpDsDIwxHUPxPkohqVeP2UOKokQt4RBAVhRFUSzGcmMgIr1EZDeAywF8ISJfWz0mRVGUaMPyALIx5lMAn1o9DkVRlGjG8p2BoiiKYj1qDBRFURQ1BoqiKIoaA0VRFAVqDBRFURQAYjxJGBcDRGQ/gB0BuFVlAAcCcJ/ijn4PDvS7IPo9kEj7HuoaY9wkHIqtMQgUIrLaGBP10tn6PTjQ74Lo90Ci5XtQN5GiKIqixkBRFEVRYwAAU60eQJig34MD/S6Ifg8kKr6HqI8ZKIqiKLozUBRFUaDGQFEURUGUGwMR6SIiW0Rkq4iMtno8ViAitUVkqYhsEpE/ROQeq8dkJSISKyK/icgCq8diJSJSXkTmiMhm2/+Ny60ekxWIyL22v4sNIvKRiJS0ekzBImqNgYjEAngNQFcATQD0F5Em1o7KErIB3G+MOR/AZQDujNLvwc49ADZZPYgw4FUAXxljzgNwEaLwOxGRmgDuBpBijLkAQCyAftaOKnhErTEAcAmArcaYbcaYLACzAPjVkzkSMMbsNcb8avv9OPhHX9PaUVmDiNQC0A3ANKvHYiUiUhbAlQDeBgBjTJYx5oi1o7KMOAClRCQOQCKAPRaPJ2hEszGoCWCX0+PdiNJJ0I6IJANoDmCltSOxjAkAHgSQa/VALKY+gP0ApttcZtNEpLTVgwo1xph/ALwMYCeAvQCOGmMWWTuq4BHNxkA8nIvaPFsRKQNgLoCRxphjVo8n1IhIdwAZxpg1Vo8lDIgD0ALA68aY5gBOAoi6mJqIVAC9BfUAnAOgtIjcbO2ogkc0G4PdAGo7Pa6FCN4C5oeIlAANwUxjzDyrx2MRrQH0EJHtoMuwvYh8YO2QLGM3gN3GGPsOcQ5oHKKNjgDSjDH7jTFnAcwDcIXFYwoa0WwMVgFoKCL1RCQeDAx9bvGYQo6ICOgb3mSMecXq8ViFMeZhY0wtY0wy+H/hW2NMxK4C88MYsw/ALhFpbDvVAcBGC4dkFTsBXCYiiba/kw6I4EB6nNUDsApjTLaIjADwNZgl8I4x5g+Lh2UFrQEMALBeRNbazj1ijPnSwjEp1nMXgJm2hdI2ALdaPJ6QY4xZKSJzAPwKZt39hgiWplA5CkVRFCWq3USKoiiKDTUGiqIoihoDRVEURY2BoiiKAjUGiqIoCtQYKFGMiDwiIlGtQ6QodjS1VIlYROSE08NEAJkAcmyPhxljZoZ+VHkRke8AfGCMUaOkWErUFp0pkY8xpoz9d5vMxBBjzBLrRqQo4Yu6iZSoRUSetOsPiUiyiBgRuVVEdonIYREZLiKtRGSdiBwRkckur7/N1vjlsIh8LSJ1vbxPSRH5QEQO2u6zSkSqichYAG0BTBaRE/b7i8h5IrJYRA7Zmi/1dbrXDBF5w/b8cRFZ5u19FcUf1BgoSl4uBdAQwI2gpPUYULCsKYC+ItIOAESkJ4BHAFwPoAqA7wF85OWeAwGUA4URKwEYDuC0MWaM7XUjjDFljDEjbFLRiwF8CKAqgP4ApohIU6f7pQJ4BkBlAGsBWO7uUoo/agwUJS/PGGPO2HTrTwL4yBiTYdO2/x7s9wAAwwA8b4zZZIzJBvAcgIu9rNLPgkbgXGNMjjFmTT4y4d0BbDfGTDfGZNsaD80FcIPTNV8YY5YbYzJBY3W5iNT2dDNF8RU1BoqSl3Sn3097eGyPQ9QF8KrN7XMEwCGwR4anBknvg4KIs0Rkj4i8ZJMN90RdAJfa72u7dyqA6k7X/NuUyRhzwvbe5/j8CRXFAxpAVpTCsQvAWF8ykmxa+E8BeMrWTe5LAFtA6XDXdL5dAJYZYzrlc8t/dwG2pkQVEaW9OJTAoTsDRSkcbwB42O7LF5FyItLH04UicrWIXCgisQCOgW4je4prOthm0s4CAI1EZICIlLAdrUTkfKdrrhWRNjZ56WcArDTGOLdwVRS/UWOgKIXAGPMpgBdB188xABsAdPVyeXWwW9gxsDnKMgD2LmqvArjBlpE00RhzHMA1YIOdPQD22d4nwel+HwJ4AnQPtQTdSIpSJLToTFGKESIyA2xJ+ajVY1EiC90ZKIqiKGoMFEVRFHUTKYqiKNCdgaIoigI1BoqiKArUGCiKoihQY6AoiqJAjYGiKIoC4P8BG08ktVqwxXkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.fill_between(range(len(Y_test[0])), Y_predicted[1][0], Y_predicted[2][0], color=\"r\", alpha=0.25)\n",
    "plt.plot(Y_predicted[1][0], linestyle=\"--\", color=\"r\")\n",
    "plt.plot(Y_predicted[2][0], linestyle=\"--\", color=\"r\")\n",
    "plt.plot(Y_predicted[0][0], linestyle=\"--\", linewidth=3, color=\"r\")\n",
    "\n",
    "plt.scatter(range(len(Y_test[0])), Y_test[0], color=\"black\")\n",
    "\n",
    "plt.xlabel(\"Time step\", fontsize=12)\n",
    "plt.ylabel(\"Prediction\", fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the black dots correspond to the original test labels, the red line corresponds to the RNN predictions and the shaded area corresponds to the sequential confidence intervals. As we can see, because the noise variance increases progressively over the time steps, the RNN predictions deviates more and more over time, and the corresponding confidence intervals widen. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tensorflow",
   "language": "python",
   "name": "tensorflow"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}