chapter12_time-series/issm-scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exponential Smoothing and Innovation State Space Model (ISSM)\n",
    "\n",
    "In this notebook we will illustrate the implementation of filtering in innovation state space model (ISSM, for short) using MXNet. Let us first briefy reivew the basic concepts.\n",
    "\n",
    "\n",
    "Time series forecasting is a central problem occuring in many applications from optimal inventory management, staff scheduling to topology planning. \n",
    "Given a sequence of measurements $z_1, \\ldots, z_T$ observed over time, the problem here is to predict future values of the time series $z_{T+1}, \\ldots, z_{T+\\tau}$, where $\\tau$ is referred as the time horizon.\n",
    "\n",
    "Exponential smoothing (ETS, which stands for *Error, Trend, and Seasonality*) is a family of very successful forecasting methods which are based on the key property that forecasts are weighted combinations of past observations ([Hyndman et. al, 2008](http://www.exponentialsmoothing.net/home)).\n",
    "\n",
    "For example, in simple exponential smoothing, the foreacast $\\hat{z}_{T+1}$ for time step $T+1$ is written as ([Hyndman, Athanasopoulos, 2012](https://www.otexts.org/fpp/7/1))\n",
    "\n",
    "$$ \\hat{z}_{T+1} = \\hat{z}_T + \\alpha (z_T - \\hat{z}_T) = \\alpha\\cdot z_T + (1 - \\alpha)\\cdot \\hat{z}_T, $$\n",
    "\n",
    "In words, the next step forecast is a convex combination of the most recent obseravtion and forecast. Expanding the above equation, it is clear that the forecast is given by the exponentially weighted average of past observations, \n",
    "\n",
    "$$ \\hat{z}_{T+1} = \\alpha z_T + \\alpha(1-\\alpha) z_{T-1} + \\alpha(1-\\alpha)^2 z_{T-2}+ \\cdots. $$\n",
    "\n",
    "Here $\\alpha > 0$ is a smoothing parameter that controls the weight given to each observation.\n",
    "Note that the recent observations are given more weight than the older observations.\n",
    "In fact the weight given to the past observation descreases exponentially as it gets older and hence the name **exponential smoothing**.\n",
    "\n",
    "General exponential smoothing methods consider the extensions of simple ETS to include time series patterns such as (linear) trend, various periodic seasonal effects. All ETS methods falls under the category of forecasting methods as the predictions are point forecasts (a single value is predicted for each future time step). \n",
    "On the other hand a statistical model describes the underlying data generation process and has an advantage that it can produce an entire probability distribuiton for each of the future time steps.\n",
    "Innovation state space model (ISSM) is an example of such models with considerable flexibility in respresnting commonly occurring time series patterns and underlie the exponential smoothing methods.\n",
    "\n",
    "The idea behind ISSMs is to maintain a latent state vector $l_{t}$ with recent information about level, trend, and seasonality factors.\n",
    "The state vector $l_t$ evolves over time adding small *innvoation* (i.e., the Gaussian noise) at each time step. \n",
    "The observations are then a linear combination of the components of the current state.\n",
    "\n",
    "Mathematically, ISSM is specified by two equations\n",
    "\n",
    "* The state transition equation is given by \n",
    "\n",
    "$$l_{t} = F_t l_{t-1} + g_{t}\\epsilon_t,\\quad \\epsilon_t\\sim \\mathcal{N}(0,1).$$\n",
    "\n",
    "Note that the innovation strength is controlled by $g_t$, i.e., $g_t\\epsilon_t \\sim \\mathcal{N}(0, g_t^2)$.\n",
    "\n",
    "* The observation equation is given by\n",
    "\n",
    "$$z_t = a_{t}^{\\top}l_{t-1} + b_t + \\nu_t, \\quad \\nu_t \\sim \\mathcal{N}(0, \\sigma_t^2)$$\n",
    "\n",
    "Note that here we allow for an additional term $b_t$ which can model any determinstic component (exogenous variables).\n",
    "\n",
    "This describes a fairy generic model allowing the user to encode specific time series patterns using the coefficients $F$, $a_t$ and thus are problem dependent. The innovation vector $g_t$ comes in terms of parameters to be learned (the innovation strengths). Moreover, the initial state $l_0$ has to be specified. \n",
    "We do so by specifying a Gaussian prior distribution $P(l_0)$, whose parameters (means, standard deviation) are learned from data as well.\n",
    "\n",
    "The parameters of the ISSM are typically learned using the maximum likelihood principle. \n",
    "This requires the computation of the log-likelihood of the given observations i.e., computing the probability of the data under the model, $P(z_1, \\ldots, z_T)$. Fortunately, in the previous notebook, we have learned how to compute the log-likelihood as a byproduct of LDS filtering problem. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtering\n",
    "\n",
    "We remark that ISSM is a special case of linear dynamical system except that the coefficients are allowed to change over time. The filtering equations for ISSM can readily be obtained from the general derivation described in LDS.\n",
    "\n",
    "\n",
    "Note the change in the notation in the following equations for filtered mean ($\\mu_t$) and filtered variance ($S_t$) because of the conflict of notation for the ISSM coefficient $F$. Also note that the deterministic part $b_t$ needs to be subtracted from the observations $[z_t]$. \n",
    "\n",
    "$$\\mu_h = F_t \\mu_{t-1} \\quad \\quad \\quad \\mu_v = a_t^{\\top}\\mu_h$$\n",
    "\n",
    "$$\\Sigma_{hh} = F_t S_{t-1}F_t^T + g_t g_t^T \\quad \\quad \\quad \\sigma^2_{v} = a_t^T\\Sigma_{hh}a_t + \\sigma_t^2$$\n",
    "\n",
    "$$K_t = \\frac{1} {\\sigma^2_{v}} \\Sigma_{hh}a_t $$\n",
    "\n",
    "$$\\mu_t = \\mu_h + K(z_t - b_t -\\mu_v) \\quad \\quad \\quad S_t = (I - K_t a_t^T)\\Sigma_{hh}(I-K_t a_t^T)^T +  \\sigma^2_t K_tK_t^T$$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import mxnet as mx\n",
    "from mxnet.ndarray import linalg_gemm2 as gemm2\n",
    "import mxnet.ndarray as nd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ISSM Filtering Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ISSM_filter(z, b, F, a, g, sigma, m_prior, S_prior): \n",
    "    \n",
    "    H = F.shape[0] # dim of latent state\n",
    "    T = z.shape[0] # num of observations\n",
    "    \n",
    "    eye_h = nd.array(np.eye(H))    \n",
    "\n",
    "    mu_seq = []\n",
    "    S_seq = []\n",
    "    log_p_seq = []\n",
    "        \n",
    "    for t in range(T):\n",
    "    \n",
    "        if t == 0:\n",
    "            # At the first time step, use the prior\n",
    "            mu_h = m_prior\n",
    "            S_hh = S_prior\n",
    "        else:\n",
    "            # Otherwise compute using update eqns.\n",
    "            F_t = F[:, :, t]            \n",
    "            g_t = g[:, t].reshape((H,1))    \n",
    "            \n",
    "            mu_h = gemm2(F_t, mu_t)\n",
    "            S_hh = gemm2(F_t, gemm2(S_t, F_t, transpose_b=1)) + \\\n",
    "                   gemm2(g_t, g_t, transpose_b=1)\n",
    "\n",
    "        a_t = a[:, t].reshape((H,1))\n",
    "        mu_v = gemm2(mu_h, a_t, transpose_a=1)\n",
    "\n",
    "        # Compute the Kalman gain (vector)\n",
    "        S_hh_x_a_t = gemm2(S_hh, a_t)\n",
    "        \n",
    "        sigma_t = sigma[t]\n",
    "        S_vv = gemm2(a_t, S_hh_x_a_t, transpose_a=1) + nd.square(sigma_t)\n",
    "        kalman_gain = nd.broadcast_div(S_hh_x_a_t, S_vv)\n",
    "\n",
    "        # Compute the error (delta)\n",
    "        delta = z[t] - b[t] - mu_v\n",
    "\n",
    "        # Filtered estimates\n",
    "        mu_t = mu_h + gemm2(kalman_gain, delta)\n",
    "\n",
    "        # Joseph's symmetrized update for covariance:\n",
    "        ImKa = nd.broadcast_sub(eye_h, gemm2(kalman_gain, a_t, transpose_b=1))\n",
    "        S_t = gemm2(gemm2(ImKa, S_hh), ImKa, transpose_b=1) + \\\n",
    "                nd.broadcast_mul(gemm2(kalman_gain, kalman_gain, transpose_b=1), nd.square(sigma_t))\n",
    "                \n",
    "        # likelihood term\n",
    "        log_p = (-0.5 * (delta * delta / S_vv\n",
    "                         + np.log(2.0 * np.pi)\n",
    "                         + nd.log(S_vv))\n",
    "                 )\n",
    "\n",
    "        mu_seq.append(mu_t)\n",
    "        S_seq.append(S_t)\n",
    "        log_p_seq.append(log_p)\n",
    "\n",
    "\n",
    "    return mu_seq, S_seq, log_p_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "We will use the [10 year US Government Bond Yields dataset](https://datahub.io/core/bond-yields-us-10y) to illustrate two specific instances of ISSM models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (12, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"https://datahub.io/core/bond-yields-us-10y/r/monthly.csv\", header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df.set_index(\"Date\")\n",
    "\n",
    "# get the time series \n",
    "ts = df.values[:,1]\n",
    "\n",
    "# Let us normalize the time series\n",
    "ts = np.array((ts - np.mean(ts)) / np.std(ts), dtype=np.double)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd829XVP/DP1ZZledvxSGLHcfbeCWGGsMsuPKSMDkYHLQXa5/lBSydtoaXQ3VJKKWUXWlbDDHsmIXsPJ3Fsx3vI1p7398d36Kst27I1fN6vV16RNW9i+fjq3HPPZZxzEEIIyR2qdA+AEEJIalFgJ4SQHEOBnRBCcgwFdkIIyTEU2AkhJMdQYCeEkBxDgZ0QQnIMBXZCCMkxFNgJISTHaNLxomVlZbyuri4dL00IIVlr69atPZzz8kT3S0tgr6urw5YtW9Lx0oQQkrUYY8eTuR+lYgghJMdQYCeEkBxDgZ0QQnIMBXZCCMkxFNgJISTHUGAnhJAcQ4GdEEJyDAV2ktXaLE5s2NeZ7mEQklEosJOs9sVHNuPGx7bA6fGneyiEZAwK7CSrtQ+4AABHum1pHgkhmYMCO8lqJSYdAOBYjz3NIyEkc1BgJ1nNH+AAQKkYQhQosJOsZnV5AQA+McATQiiwkyxndfsAAP5AIM0jISRzUGAnWY2LE3WasRMSRIGd5AQ/BXZCZBTYSdby+oPpF5qxExJEgZ1kLZc3WAlDM3ZCgkYc2BljBsbYZsbYTsbYXsbYT1IxMEIScSoC+31vHMSPXtqTxtEQkjlSMWN3A1jDOV8AYCGAcxljK1PwvITE5faGVsL889OkjoMkJOeN+DBrzjkHIO3n1op/6HMxGXXKGTshJCglOXbGmJoxtgNAF4ANnPNNqXheQuJxUWAnJKqUBHbOuZ9zvhDARADLGWNzw+/DGLuJMbaFMbalu7s7FS9LxrlobQTOvP89vLj9RBpGQ0jmSGlVDOfcAuBdAOdGue0hzvlSzvnS8vLyVL4sGaeipWKOdNtx6792pGE0hGSOVFTFlDPGisTLRgBnATgw0uclJBGXl9oIEBLNiBdPAVQB+CdjTA3hF8WznPP1KXheQmIadHnxtSe2pnsYhGSkVFTF7AKwKAVjISRpbRZnuodASMainackK/n8VFFLSCwU2ElWctDBGoTERIGdZCXpgA1CSCQK7CQr2cQDNgghkSiwk6xkdcUP7AHq9kjGMQrsJCtFm7FfuXSifNnti13j3mtz43ivfVTGRUgmoMBOslK0HPsvLp2HdcsnAQDcvtiLq6f/+j2cdt97ozU0QtKOAjvJSrYoqRiNWoX5E4sAxN+VKqVxPHFm9YRkMwrsJCtZxVTM+m+dHHK9XiO8pWN1flQG8xO0yYnkqFS0FCBkzNlcPsysNGNuTSGeumEF/FxYLDVo1QBi59gdnuBMv33AiSllptEfLCFjjAI7yUpWlw9mg/D2PamhTL7eoI0/Y1emaBxu2uREchMFdpJ1jnTb8OnRXuTp1BG36TXCdbEDe/B6Bx3UQXIU5dhJ1tnfPgggelsBacYeKxXjUlTLOD20yYnkJgrsJOswMADAQ9cuibgt8YxdkYqhfjMkR1FgJ1nHLs60Z1UVRNwmLZ66Ys3YlakYCuwkR1FgJ1lHOuvUpI9cIpLKHd1J5NijnZlKSC6gwE6yisvrx49e3gsAURdPE8/YKRVDch8FdpLxNh/rwwMbDgEAjvUEe7xIs3MlefE0xoxd2Wog2mHYhOQCCuwk413510/x+7cPAwAOdVrl6xljEfeVFk9jVcW4xRm7TqPCoJN6upPcRIGdZJVuqzvu7Vo1g4rFroqRZukN5flo7XekfHyEZAIK7CRr+AMcHn/8xl2MMRi06piBvd/hAQDMn1iI430U2EluosBOsobXH5BTKR/87xkx76fXqGKmYnpsbhTnaVFbaoLF4YWdTmIiOYgCO8kaPnHGrlOrMLk0L+b94s3Ye6welOXrUWLSAgjO4AnJJRTYSdbw+oQZuy5KNYySENijz9h77W6U5utQlKcDAFgctIBKcg8FdpI1vIEAPH5/1DJHJa2awRsjF99n96DEpEOJSSd/TUiuocBOsobPz5OasWtUKnj90Q+ztrl9MOu1KM6jVAzJXRTYSUbjPBigvf4APP5AUjN2XyD6jN0m9nGnVAzJZRTYSUbzB5SBPckZu1oFX5QZuz/AYff4kW/QwKQT+szYqXUvyUEU2ElG8wUiZ+yJUzHRc+w2sbTRbNDCoFWBMWoERnITBXaS0ZSB3efncPv8ctuAWDRqFvI4idUlpF3Meg0YYzDpNLDT8XgkB1FgJxnNr0ipeMQNSjp14sVTZWB/dXc7rn/0Mzz84TEAkM9KzdOp0dxnR90dr+CtfZ2jMHpC0oPOPCUZTbkI6hNTMfmG+G9brZrBp0jFfOPJbSG3VxToAQj93Dcf6wMA/PPTJqydPSFFoyYkvUY8Y2eMTWKMvcsY28cY28sY+3YqBkYIEJaKCXDYXL6ofdiVNKrQxdPwGX5DuRmAMGMfdAl5d7UqslMkIdkqFakYH4DvcM5nA1gJ4GbG2OwUPC8hIYHd4w+gfcCFqkJj3Mdo1AxexUxf2owkKRRr2KXKGEBYcCUkV4w4sHPO2znn28TLVgD7AdSM9HkJAUJz7D1WN5xeP6oKDXEfow0rd6wqCt5/5w/Pli/n6YMzf1WU3u6EZKuULp4yxuoALAKwKcptNzHGtjDGtnR3d6fyZUkOU868W8Q2u9VFCWbsqtAcu0fs9HjKtDJ5tg4A82sKUzlUQjJGygI7YywfwH8A3Mo5Hwy/nXP+EOd8Ked8aXl5eapeluQ45QalbptwyEZ4aiWcRq2CV/E4m9uHSxZW4/HrV4Tc79y5VSH3ISRXpCSwM8a0EIL6k5zz51PxnIQACEmpdFuFvi6FRm2suwOIrIqxunxRK2lmVJrly3baqERySCqqYhiAvwPYzzl/YORDIiRIWe7YI87YEwV2tYqF/EIQ+sNEPkatYthw26lYPqUETmotQHJIKmbsqwFcC2ANY2yH+Of8FDwvISFVMckGdq1aJefm3T6/UPuuj177Pm2CGTVFRjhoxk5yyIg3KHHOPwJAJQVkVChz7L02DzQqlkQdO5MfZ3VJ/WFiv9WNOjX1jCE5hVoKkIymbObl9PpRYNSCJShN1KiFfuycCxuaAMScsQNAnlZNM3aSUyiwk4zmD2vmZdTGn60DgFbcbOQP8JCOjrHk6dRwev0IRGkcRkg2osBOMlp4l0a9NvFbViO2EPAFuJyKiTdjN4o7UF0+mrWT3ECBnWS08AMzErXsFe4jvK1dXn+wVW+cHLuUs6d0DMkVFNhJRvOHHXGX6Fg8ACgSd5daHF45FRN/xi4EdlpAJbmC2vaSjGYLOwgjmcBeLO5M/fWbB+EW2wnQjJ2MJxTYScbyBzi++9xOAICKAQEOGJJYPC0WD6pev6tdvi5eD3epy6ODNimRHEGpGJKx9rcHWw6VmITDMZKaseeFVsDo1Kq4uXlKxZBcQ4GdZCy3okpFCuj6JGbsZfn6kK8TnbhEqRiSayiwk4zVb/fKl7VqoTY9mRm7Sa/BpYuCRwL02T1x7y8Hdi8FdpIbKLCTjNXnCAZkqTY9mcAOAPdcNg83njJFeGyC05GkOva/vHcEL+9sG85QCckotHhKMla/ONN+/Prl+MWrBwAkt3gq3e/7F8zGtAlmFMTZdQoILQUAIad/y9PbcdGC6hGMmpD0o8BOMpbN7QNjwMkNZUNKxShduXRSwvsYw5qKOTw+5OnoR4NkL0rFkIzl9gWgU6vAGEOACztQTXE2Gg1X+C+LHmv8nDwhmY4CO8lYHl8AOjHoinE97kaj4WKM4Usn1eGalZMBAFa3N8EjCMlsFNhJxnL7AnL9uUps1TsagR0AfnzRHJwvnoEqtfolJFtRYCcZy+MLyGkSlVjZYtbHXwgdCanenQ62JtmOAjvJWB5/MBUjrp0m3Gw0ElKjMCvN2EmWo8BOMpbH54dOrF9XizP20TyDUTqMw0ozdpLlKLCTjOXxBeSDNU6fUQEAqCo0jtrrSfl7yrGTbEfFuiRjSeWOAPD106bi0kU1qC4avcCu16igUTHYqCqGZDmasZOMpSx3VKnYqAZ1QCh7zDdoKMdOsh4FdpKxlIunY8Vs0MRNxXDO8dyWFrQPOMdwVIQMDQV2krE8ilTMWMnXa+Munh7ptuF//70LP3hxzxiOipChocBOMpaweJpc069UMevjz9i3NPUDiH6otj/A0dzrGLWxEZIsCuwkY7nTMWM3aOK2FOh3CLcVGCM3Sv1mwyGcet+7OGGhNA1JLwrsJOMc6rTil68fwIDTO2otBGIxGzQYcMYO7E7xXFQuNq9xevy48bEtaOyy4q39nQCAY9320R8oIXFQuSPJOGf/5gP5crlZH+eeqTe3uhAv7WjDkW4bppbnR9zuFE9Zkv7+qLEHG/Z1wuMLyP1sjvXYcPK0srEbNCFhaMZOMsrfPzoW8nV5/tgG9hX1JQBiz7qlc1GdHj/8AY5Xd7cDAHa1WrBPPHy7bcA1BiMlJDYK7CRjNPc6cPf6fSHXlZl1YzqGYFuB6OkY5Yz9sU+b8ML2EwCCuXcA6LW5R3eQhCRAgZ1kjI8aewAAD16zWL6uttQ0pmOQGoHFqoxxijP2Dw/3YMvx/pDb1syswOyqAvTa6KAOkl4pCeyMsUcYY12MMSruJcN2rMcGnUaFRZOL5eui5blHk7RYG6uWXZqxA8Bnx/pCbptTXYDSfB167JGBfcDhxVObmuVFV0JGU6pm7I8CODdFz0XGqeY+ByaX5KHCrMfXTpuKV285ZczHIPeLiTFjl3LsQDDI15bmARB+CZWYdOizh6ZiXF4/vv/ibnzvhd3Y1TowSiMnJCglgZ1z/gGAvoR3JCSO471CYGeM4Y7zZmJ2dcGYj4ExBnOcfjEurx/1ZUJ6yOryobrQgOtW1QEAZlSaYTZoYHcHg/8nR3ow8wevY/0uYZG1uY82MJHRR+WOJCNwznG814GTpqa/TFBoBBZ98XTQ6cW8iUXotrphdftQlKfDV1bXYe2sCtSWmmDSaWBXpHF2tFhCHn+UatzJGBizxVPG2E2MsS2MsS3d3d1j9bIkC2w+1ocpd74Kp9ePurK8dA8HRq0aLm8g4vpXd7ejqdeBsnwdik1CtU6xSQvGmLzIm6fTwO0LwOcXHt/aH7oLtdtGpZBk9I1ZYOecP8Q5X8o5X1peXj5WL0uywG/fOiRfnlyS/sCu16jh9vlDrgsEOL7x5DYAQFm+HgVG4cNucV5oOaZJL/SQcYj596PdtpDbqWKGjAUqdyRp51JUmtSNcXljNHqNCh5/6Izd7Qt+XZavkw/VDg/seToh4Dvcfvj8gYjFUgrsZCykqtzxaQCfApjBGGtljF2fiucl2ecHL+7BjY9tGdJjlL1ZaopH9zCNZOg0KrjDUjHKMsfKQiNMYr27lJKRSDN2u8eHln4nHB4/zphRjqpCA2pL89Bjp81LZPSlZPGUc74uFc9DspvL68fjG48P6TFWlxfHeoQFxe+cNR3aMe7mGI1eo4qoinF4gl+vnlqKf3wstD4oyQvt8qicsdvERdQbTqnH6oYy/PClPXh5Z9toDp0QAFQVQ1KoR7GV3u72ybPaeBq7bAhw4OHrlmLt7AmjObyk6TQqeHyhM3YpXfSHdYugUavQbRX+rXVloakjacZuc/vQOSgslE4oMAAASk16WBxeeP0B/GbDIdjdPvz4ojlgYvMwQlKFAjtJie+9sBtPbWqWv+62upMK7FJdt7TJJxNEWzx1eoRAbxQP/pAC/7QJ5pD7STn3hz88ivpyIehXFoqBPV+4rc/uwZ/fOwIA+M45M1BgiOztTshIUGAnKaEM6gDQOeiKmM1GI6VhJmVANYxEH2XGLqVijDohsP/p6sVYv6sd1WLQllSKs/O3D3Th7QOAQauS+8+UiYH9mc0t8v0tdi8FdpJy6U9okqzV0ufAa7vbwTmHVh2aTvisKfFGZM453j3QhWkV+TCM8RF48eg0qpAqGCC4eCoF9ukTzLj9rOkRaZSisJy7smqmTGxB/BtFeafFSVUyJPUosJNh+/Yz2/H1J7dh6/H+kL7ps6oK8NSmZnmTTiyv7u7AztYBnDkrM3LrEiEVEz3HbkzwCygy0AcD+8JJRRH3V7b7HYo+uydkQZcQJQrs48iRbhve2teZsueT8uO7WgdQkh8MYNeurEXbgAuHu2zYGtbaVum/O9tQlq/Hd86enrIxpUK0xVNnkoEdAF695RRcsrAaAGDSBe+vUatQUySUc9598RwAgMUxvBn74rs34KI/fjysx5LcR4F9HLnpsS244bEtckXHSPkCQgvagx3WkLrvKjHv/H//3oXL//IJtsRIy5ywODGnuiAjShyVpA1KgUCwxa7U1dGoSxzYZ1cXYI34KaQvrIWvNKGfWSU0OBvOhiW/OK7GLluCe5LxKrN+osiosLt9+MaTW3FEbEAVK9AOhdPjh0VMI7QNOGFz+6BRMfzl6sVyed/uE8Kuyzuf3w1vlLRM56ALEwrG9ui7ZOi1wo/Fkp9tQJdYsigdsJHsWsAZM8rRUJGPdcsnh1z/xy8sxjlzJmDhpCLoNCq5JHIolGWlFNxJNBTYx4GXdrTh1d0d8tcW5/DyukodioDUZnHC5vLh2lW1OG9eFaqLQitFDnfZsKs1tMuhP8DRY3PLvwQyiZRu6Xd4sa1ZSCVJOfa8JGbsgHDE3lu3n4YbT60PuX7hpCL89dql0KpVqCo0oH0Y56O2WYKNxd7en7rUGskdFNjHgabe0FaxlmEu2Cm1DwjBZU51AZr7HEILW6OQZy/K0+GuC2aF3P9IWLvalj4HAhyoLkp/C4Fw+Yr6e2kx1OHxQ6NiKU0bVRYY5P/HoVD+Muih81VJFBTYx4Hw7fHhJ/wMxy1P7wAALKkthtcv5HxDFlBX1cqXtWom16tLNh7tBQAsrS1GpjEr6sqlk5ScXn9SC6dD0VCRj/3tVnj9AZywJB/gpRl7gUGTsvUSklsosI8D4RUef/vwGAZjHCSRjJY+hzxTXKw4n7REUdqn16gxrSIfN58xFZNL8nC024ZAgMPnD4Bzjue2tqK+zISGirE90zQZ0rmnAOR+Ly6vH4Yk0zDJWj6lBDa3D/e8egCr731H/mWXSPuAC0atGg0V+eimGTuJgnaejgM2d2QQ39M6gJMakj+tiHOOzkE3KgsN2Nc+CAD43VUL5QMmAKAkrNPhhttPAwAc7LDhaLcdd720B09tasbqhlJsPd6PG0+ZkpF9UpSB/Ucv7wVjwuJpsvn1ZC2rKwEAPL6xCQDw7GctWFlfmvBx7QNOVBUZUGE24GgPLZ6SSDRjHwdsiqPaFkwsBABsHmJlzF8/OIqV97yN4712dIkf/1fWl2KGoleKdPhEuIaKfBztscttBz5uFGammdRGQCk/rMfND1/aC4cn9amY6iIjyvL1cirr4yM94JwneBTQZnGhutCIcrNe/l4QokSBfRywKXLsF8yvwqr6Uvlw5WS9uVeoqmmzuNA96AJjQKlJB6NOjUM/Ow+PfGkp5lQXRn3soslFcu11yPWTMi+/DgD6KAG8bcA5Km0PyhTrEp2D7oRVMrtbB7CjxYISkw7lZqFbZHjDMkIosI8Dyhl7qUmPqRWmiI0ziUibkawuLzoH3SjL10MjVojoNCqsmRm7LcDMyuCs/qun1WNCgR7vfvd0zJsY/RdBulUXGvC106aGXLfnxCAM2tT/uEgdH6UdqXvbBuPe/0/vNgIAvP4Ays3CHgA6lYmEo8Ce457Z3BxSaliar4PZoIXV5U3qY7/EJ6YLemwedFldqDAnv7FIKoMEgOtW1WHT99ZiShKdH9OFMYY7zpsZsWZw2aKJKX+tEpPw/zirSvjl1x/nF+7RbhteFz853fW52XInyaFU1JDxgQJ7jrvj+d0hX5sNGpgNGnj9PKLRVTzSztEemxudg+4hBfZ8xWJkaViwzGQv3bwav7p8vvz1OXMrU/4aUj5fOus1XrWS8vzUmiIjZoltCfYlmOWT8YcCe46bLC5Q/uXqxQCA+rJ8mMVgMpSSx15xJtnS58C+9sEh7RhVq4KVL5nUnjeRSSV5uHLZJNx9yVx8bn4VCo2p75supalMeg0YCz3/NZzUfuDRLy8DAEwo0KMoT4vDXdaUj4tkNyp3HAcuXliN8+ZVoeneCwAEN+DYXD5UmOM9UuDy+uWc/HNbWwEg4xp3jaZrV9bi2pW1ie84DFevmAxfgOOyRTV49JMmDMYJ7K39ThQYNDh9RgUAIWVUnKfDgDN0A5rPH0Cv3ZOR7RrI2Bg/P53jlNXlDanLBoJ12uE7UmM51Bk5I7wprAcKGR6NWoXrT56CYpMOhUYtOgfd+OXrB6J+mhpweiPy/gUGTcQvg0c/acKKX7yNPWITtqc2NeOyP1OL3/GEZuw5jHMOm9sXskUeCM7Ykw3sL25vg0GrwpLaYnzc2IspZaYh16DfcuY0aFSZtxkpkxQatfLiqNXlxc8umRdyu83tC1mvAIACo7AQriTtYP3cHz7CG7eeiu+9sDv4+CTOoSXZj2bsOayp1wGvn8eZsXvhD3A88OZBuT1tOJ8/gCc2HkdDRb5cdx7+fMm4/azpuOXMaUN+3HiiTJ3sViyUSmyuyMBsNmgwGPYLukdR/qg86KS515GqoZIMR4E9h93/5kEAQF7YgqUyFbOjxYLfv9OI//efXVGf49FPmuDxB+DxBXD1yslYUluML51UN6rjHq+qFAdjR2utbHX7kK8P/fRVYNBGpGJa+oIBXPkpqbkvtBFbPJ809uClHSeSvj/JLPS5LIdJC56XLwmtv5ZSMYMur9wdsKU/shaac453D3YBAH584RxUFRrxn6+fNJpDHtcqFYG93eKCP8BDKopsbi/MhtDVbiEVE5yxu7x+uYIJAN4/1C1fPj6EGfsXHt4EALhoQXVG9vMh8dGMPYc1dtlw+eKJETl26eO81eWTZ3eBsC3/dz6/G4vv3oAjXXZctrhmSA3DyPBMLA72pvf4A+gN69xod/sjUzF6DZxev9zKV2pJ8KvL50PFgFd2B1tHHO8beipmOAeBkPSjwJ6jvP4Aum1u1BRHHmShVjHk6zWwuX1yxUv4pOzpzc3od3jRMejC1PLMa62bi8L/n3vCWgXYXNEXTwHgpHvfwbEeuxzgJ5XkRRxi0jKMwP6Hdw6HtKQg2YECe47qsrrBeWjeVilfL5TJSbsZwxfglCmA+Rna0yXXSIFdahWgPB3J7fPD4w9EzNiVHTVP9DvlwF5TZERpfnB3cE2RcVjnqz69uQV3hu1eJpmPAnuO+vuHxwCE5m2VaoqNON7rQJt4NFv4ApxJ0Xt8xZTEPcLJyBl1arz73dPxD3Fnaa/ipCupQ2d4RVKBIs026PKizSIE7wmFennR/NJFNThzVgUOddqSOoovPC13jHq+Zx0K7Dmotd+BRz4WAnttjHrzudUF2Ns2IH/MdvsC8oHNnHM4xcuVBQboNPQ2GStTykxy+qzHGkzFSN+nyHLHYGDvtrrRPuBEWb4eeo0a/Q7h8WfOqpB7vl8tLorGE745SmoAR7IH/cTmoJ0twRpo5QlHSnNqCmH3+MG50KYWEBZTH/u0CVf+9VN4/RxfOqkOb9x26lgMmSiY9RroNCr0KGbsUuVLeGBXthLuHHTB4vCiOE8I9lJb35X1pVhWJ+xBONqduOSxIyxlIzWAI9mDyh1z0K5WCwDgulW1IblypbmKQzGqioxoG3Chc9CFH760V77+4oXVo9L4isTHGEOZSRd9xh6WiplWYcbK+hJsPNqHI9022Nw+eUH1N/+zEIc6rCjL1+PSRTU41mPHH95phMXhQVFe7C6bP3hxT8jXXpqxZ52UzNgZY+cyxg4yxhoZY3ek4jnJ8NjdPjy3tRULJhXhpxfPjXm/2tJgimaaeKD0bzYckq+7YslEzJ9YNHoDJXGVmfX4z7ZWvCqWK8o59rANSkadGs/ctArnz6vEwQ4rBp0+OQ9flq+Xy1QZY/LB441d8XPmnzX1h3zd3OeIKL0kmW3EgZ0xpgbwJwDnAZgNYB1jbPZIn5cMz84WC/rsHpwzJ/aJRgBCDmaWWse+fUDYjLT/p+fivisWxJztk9GnF9c1vvHkNgCxZ+ySuTWFaOp1oKnXHrFvQdIg/gI/HCewH+kO3ja1PJjGOxilERzJXKmYsS8H0Mg5P8o59wB4BsDFKXheMgSHOq0YdHnlXYdnzYof2JW7CesUpxl9fslEGHXZ0zM9V926drp8ucfmlncIh3d3lKyYUgJAyMXH6uVTU2REnk4dtVsnABzutOLM+98HAPzjS8vw7FdX4d7LhEZkJ6LsTCaZKxWBvQZAi+LrVvE6MkYCAY6zf/MB1j20Uf7IHCsARKMMBN9a05Dy8ZGhW91QhqdvXAlAOAe1td8Bs0ETc82jviy4uakgxoxdpWKYU12Ajxt7oh6L+MvXD8qXT5tejtJ8PS5bPBGMJT6LlWSWMauKYYzdxBjbwhjb0t3dnfgBJKEBpxf3v3lQbvW6t20QfXYPGEPcxbFwVYVG3LZ2OmqKjPKJSyT9ZlcLR99t2NeBf356PO7hJkV5wWBeYoq94H3OnEoc6rShO0rOvLJQqKK5ZGE1VGIaTqdR4cyZE/DU5uZh7Vwl6ZGKwH4CwCTF1xPF60Jwzh/inC/lnC8tLy9PwcuSN/Z04A/vNMp5WADosXtQnKcbUn68qtCAb6+dho/vWEMNnzJIoVGL2tI8PLGxGQBwzYrJMe+r/L5NmxD7WCzpF/fmY30hO1sBwC9uTLrnsvkh1991wSx4fAHc+9qBof0DSNqkIrB/BmAaY2wKY0wH4CoAL6fgeUkC0U7Z2Xi0Fw1J9nZ57mur8LfrllIwz2DS97Kq0IDbz54R976XLKwGAMyqLIh5nwqxXcE3n9qO1fe+I29KA4B+uxfTKvIj1ljqykxY3VCKExbKs2eLEdexc859jLFvAngDgBrAI5zzvQkeRlKgyxr5cfpotx3nzKlM6vHL6kpSPSSSYmViv5dk9hP8+ooFuOnUqTHbSADCAdgSty+AQ51Wuay1z+FBcYwKKs6bAAAbGklEQVS1mSKjDh0DlGfPFinJsXPOX+WcT+ecT+Wc/zwVz0kSk049KsvX4a/XLpGvXyLWK5PsV5ovBNqCJAK7Rq2S8/KxVJgNOLmhDN8WT7P68HCPfJvF4ZF3rYaTuoGS7EAtBbJYc58DK+tL8Nn31+IM8eR6AFhSS4E9V0gdGlO1pUCtYnjihhW4de00rKwvwdObm+Xb+h2Rh2VLzAYNrC4fLA4PXt/TkZrBkFFDgT2LHe2xo748H4wx6DQqPHjNEvz9i0tjfpwm2Ufq9+IPpHZbP2MMa2dNQGu/E12DLnDO0W+P3WrAbNDC4fHjrhf34GtPbI1ZC08yAwX2LBFed/zm3g5YHF551ygAnDu3Emcm2JhEssuamRU4ZVoZrl1Vl/LnniruRG21ONHv8MIX4CiJEdilHa9SyeMH4pF7r+9px9UPb4SPGoVlFArsWaDH5sacH72BB98/Iv8Avb6nAwUGDa5cOinBo0k2y9dr8Pj1K3DRgupReW5A6C/0nWd3AAC06ug5nzIx1y91fnxqk5DCeWDDIXzc2IvXKD2TUSiwZ4ED7VY4PH7c+9oB3CV23tveYsGqqaUwaGn7Pxkek04I7DaXD/vahYqXJbXRK6VOn1EBg1aFzkGhEqvV4oTL65eP33t2S0vUxw3V+l1t2LCvM+S6xi4brFFKe0lsFNizQEt/cMffv7e2wuHxoanXjtlVdGQdGT5pxm5z+1BZYMAp08owL8YxiIVGLc6eHSyj9fgCuOLBT/HeQSElk4pdqV1WF7751Hbc+NgWbG8Odphc+8D7uPiPH4/4+ccTCuxZoFnxQ5Nv0OBQpw2cAzOrYu8wJCQRKW9ud/vQMeiKeT6uROoOWSPO0nefCB7o0iEuwI7EtuMW+bLUm0baQHW0J/EBISSIAnuG45yjuc8hH3A8sdiIgx3Cm165cErIUJn0QhpvwOlDt9Utv8dikVoJl5n1Ebe5vAEMOEeWLjmqOFv1mBjIpeP9gNRXBuUyCuwZrLHLhvk/fhOv7GrH9Eozrlw6Ed1WNw50WJGnU2NSMTXsIsOn16ihVTP8+b1GBDhQWWiMe/+l4k7lyxeHNm9dK1ZihR+pN1Rdg26Y9RrMrDTLgb3XFgzsFkWQJ/FRYM9g/9nWCqu426++zIRysx49Ng+O9zowuSRP7sBHyHB5/Rxun1BpJXV3jGVJbTF2/PAsrFs+GVJ7oa+snoKvn14PAGgfGGFgt7pQUaBHfbkJTT12fNzYE5Jr77NTYE8WBfYkWF1ebG/ux7qHNuI18aiysaA83OCrp9WjPF8Pf4DjQPtg3H4ghAzHhASpGEBoB61Vq1BqklodaOTHdcYJ7Nua+zHrB6/HrZ7pGHChwmxAXakJR3vsuPrhTfiB4gzeXgrsSaPAnoSvP7ENl/75E3x6tDdlZV2J+PwBvLyzDQsmFWHPT85BVaER5WbhB6htwJUwH0rIUEmLo8koFjcymQ1aVJgNYCz+jH17swVOrx/viscvhuOc40i3HXVlJvzPsuh7M6I1vSPRUWBPwkeNwUZJvjFawHlH/AGoK82Ty9KUfTxoxk5S6b7Pz4dek/yeCJWYi5k/sRA6jQqFRm3IQmc4KT8erdW0zx/AG3s7MeD0YnaVGbWlJnlDFCCkIQFg87HepMc33lFgT4KyZepI84jJkg4c/vGFc+TrltQWy2M5b27VmIyD5Lbl8oLoxCE97t7L5+GalZPlTqJFRi0sjthVMVLQb7cIPz+7Wi24e/0+cM7x4PtH8LUntgIAZlYJ3Sn/dt1SNFTko640D/dfuQBmvQZPbGxGP6VjkjLifuy5zusPQK1iqCo0YP7EQnx0WDgvcrQPpzjYYUVVoSGkoZdOo8L2H5yFfodH7vpHyEg88uVl6BhwDnkhftHkYixStIcuNGphCSt3XPGLt7BmZgXuuWw++sWgf8LiBOccX318K9oHXLjhlCl4/1DwqMwZYgnvosnFeOv20+Trr1lVi7+8dwTtAy5qcpcEmrEn8OTG4+ize/D9C2ZhWV0J7B4/7nvj4Ig3Y8QTCHBsPNorl5cpqVSMgjpJmXy9Bg0VI98PUZinw0BYKqZz0I2nNwtrUlIqxu0LoMfmkc9vbeyyoVvMnV8wvyrmQdynThOO07Q4acaeDArscVhdXtz3xkEsryvB5+ZXyzvu/vzeEbnOdjT8d1cbuqxunD2bOjWS7FAUNmNXTnxuf3ZHSJrmhMUprxdtb7agqdeB60+egj99YXHs5xcPABmIk+4hQZSKieNwlw12jx83nDIFADC5NLghqKlX6IU+XDtbLLjz+d04Y2Y5JpfkobbUhHcPdGHA6cUzn7VgVlUBLphHeXSSHYryQnPsUm08ADy/TTjbfmalGQc6rLjtXzvkidEDGw4BSHyQiBTY+ymwJ4UCexxN4ptPCuCzqwrwt+uW4sbHtuBotx1rZg79OTnneGH7Cdz+7E4AkLvqSTQqBp1ahR9cMIs2IJGsUWTUYtDlhT/AoVYx2MWNdbesacDv32kEAMytKcSBDmvUT7trZsb/dCqVV/bZqeQxGRTY4zjSbYNaxTCpREjBMMZw1uwJKMrTRrw5ff4AnF4//ruzHa39DvzvOTOiLrDevX4/Hvn4GADg55fOhUbFUGjU4vU9HbjlzGkj+hRASLoU5unAuXBOwIJJhQiIE/ZJJXmYXVWAfe2DKDJqUVlgQMegCz+6cDae3dKK/e2D+P75s7Bqamnc5zdo1ZhZacb6Xe04a3YlzvntB1j/rZMxt4Y6nEYz7gN7IMBhdfsiToH/5ydN+NO7RzCnuiCivndKmQlHu0MD+0/+uw+Pbzwufz2j0oyNR3txxdJJWDy5GG6fH79/+zAe+fgYLlxQjV9fEVo3fC6VL5IsViT+/Nz81DYU5Wnx9I0rAQiLs6fPKMe+9kEYdWr891snY2eLBWtnT8C2Zgv2tw+iwJhcGDq5oQxPbmrGh4eFKpqHPjiK369bNDr/oCw3rgP7J409+MLDmwAAR39xfkjq40cvC1uZl0Y5GHpKmQkfHu7BX947ggsXVGFicV5IUAeA2/61AwEOPL25Bb++YgF++t+9GHQJH0+/taZhSJtBCMl0yomRxeHFeb/7EACQp9fItzk8fpSb9VgrFgV8a00DDnYMJn2cY7FJB6fXL/eMaU5BD/hcNa6rYu4XF24AoEeRu+OcQ6dWoShPizvPnxXxuPoyE7qtbvzy9QP46uNbsas12Ed6waQizK4qQIALJ8IX52nx3ed2ykF9y11rMX0CtdsluSVPF32ikq9Xy33fw09Bmj7BjDdvOw1lSZbvSnl2qQ/8/vbBmGetHu22YXfrQNTbxoNxG9i9/gD2nBjA5BKh0qVDsaP0ey/shscfwB3nzox69NyUsmAefG/bIC4ST3fZetdavHTzalyxVNjF99QNK/DIl5YBAC6YV4Wmey9I+k1MSDZZMKkIFy6oxlu3nxpyfZ5OI09kRrp+VCxWxuxoESZSbl8A33luZ9T7rrn/fVz4x49G9HrZbNwG9pY+B9y+ANbMrAAQbBWwv31Q3lSxJEoaBgBW1pfglGlluPviOSHXSxuHrltVh8++vxYr6kvFHXSn4t7L543WP4WQtDPpNfjDukVoqDCHtCfI12uwrK4EL968GjeeUj+i15B2nFpdPrmg4aUdbQCEn9uXdghllb9961D0JxhHxm1gPy7m56Tg3WMTUjFHuoUeLRcuqI7Z7a40X4/Hr1+Ba1fV4ZVbTgYAnDGjXL5drWIoV5wy01BhhjnGjjpCcs2XV9fJl6UUzcJJRVCPsHx34aQi+fIZMyrky5xzXPTHj/DtZ3bA6w/gt28dHtHr5IJxt3ja2GXFo580yemUOdVC0yGp7rZV7IF+z2XzkuoHM6e6EMfuOX/Ue8cQki2UC6kmfepCjEGrxtRyE45023HS1FL02j14ZVc7rG4fvH5hp+udz+8OeYzb54deo8b7h7qhU6sSllXminE3Y7/5ye14YmMzntp0HKUmHepKhZagNrdwaG5rvwNFeVq5VW4yKKgTEqRsLy2dk5oqHnGxdHKJCWtnCbP2rkG33Ob331tbAUDetT3oFCZsX3xkM9b9bWNKx5LJxl1gbx8QZuRHuu1YXFsMlYohX6+BTaxa6bF6UE4LnIQMm3KWnupJz88vmYeZlWbUl5tg1ktllL6IogQpxSrthpW8c6AzpePJVOMqFfPWvk657BAAPjdf+K2er9fA5hZKsfrsnpAZByEkc5w6vRynThfWs7TipwGvP4BBRQOyAoMGU8qFT+IWhwdt6uD89U/vHknYviAXjJvAzjnHDY9tCbnunDmVAACTXg27mIrpc3gwbQhHhBFC0kMnBmyPj8Pi9GLd8kkwG7RY3VCGCQXCDL59wIXjvcGNTOPl5LFxE9ilrnDnzJmAn148F1aXT65RzzdoYRUXT2nGTsjIvXnbqREbklJNpxHSPHa3Dw6PHzVFRnxzzTQAwIA4g2+3uNA24IRJp8bs6gK593uuG1GOnTF2BWNsL2MswBhbmqpBpdqGfZ247M/CJqLLFk/EhAJDSCljkVGLbqsbPn8AFodHPoGdEDI80yeYsaQ28qCYVJIO6+gWS5WV1TgFBg1MOjXaBpywuoReUBUFBgrsSdoD4DIAH6RgLEm5e/0+PPzh0SE95sbHtqBJ/DgmHZahNK+mEIc6hXaiAQ5UFkbehxCSWXRijr1HDNaFecEJGWMMVUVGtFtcsLt9MOk1qCky4oTFGbKYmqtGFNg55/s55wdTNZhEWvoc+PtHx/CzV/YP6XFSvg1A1E1HZ8wshz/Acbf4vDXFFNgJyXTSjL1LDOxFYR1aqwoNaB9wwiYG9obyfHh8AexU9HbKVWNW7sgYu4kxtoUxtqW7uzvxA6I43GWVL7ck6Oz28IdHsfaB9/GVRz9D56DwjZ9ckhe198uS2hKcNr0cH4iH6laPkwUWQrKZtHgqpVfCW29XFxrRNuCCze2D2aDBnBphM+Kf320c24GmQcLAzhh7izG2J8qfi4fyQpzzhzjnSznnS8vLyxM/IAppVygA3PvagZj38/oD+Nkr+9HYZcM7B7oAADedWi9v/4/mvLmV8uUpZaZhjY8QMnakGfvrezsAIKLooarIgB6bGxaHFyadBnOqC1GWr5d3qSbyzoFOXP/oZ3B6/CkZ79FuG+55bT8cHl/iO49QwqoYzvnaUR9Fklr7ndCpVVg2pRifNfXFvN++NuG4OcYA6Uzd1Q1lcfu1nDlrAmZWNuEnF82BRj3u9m0RknV0YbtaJ4alUKsLjeAcONZjx+LJwoalhgoTnN7Egfqzpj585VGhPPrlnSfwP8smj2isXn8Aa+5/HwCweHKxXGo9WrKq3HFaRT6uWCpUtXzc2AuX1x+RWhlwePH9F3eDMeDTO85EgVGDl3e04ZSGsrjPXW7W4/VbT417H0JI5tCqg7ta186qiNjlWlUUTKnm64U4YdCq5YM6YuGc44oHP5W/fnJTMz6/ZNKImpj12oTXPH1G+agHdWDk5Y6XMsZaAawC8Apj7I3UDCu6K5ZOws8vnSf3UP/nJ00R93li03HsOTGIn140B5WFBuTpNLhq+WQ6GJqQHKNVfLK+7azpEbdXKyrgpG6rRq06YWol/GSmXa0D+OsHR0YyVHRZhbbgV6+oHdHzJGukVTEvcM4ncs71nPMJnPNzUjWweOZPFA6w/VtY2eOBjkH87q3DWDOzAteuqhuLoRBC0kSnCOzRNhVOKTXJ95Eq3YxadcJUzMGOYJHGdatqoVUz/Or1gxhwDH/DVZdYwFFhHps+VFmZTK4vz8d3zpqOHpsHvbbghoPntwmN9u/7/Px0DY0QMkaUn8KL8yIDu0rFsHCy0MO9xCQEVINODVeCwH60Rzio/pqVk3HnebOwSMzPr7znbTT3Du+c1fZBYcY+oWBsKu6yMrADkBsB/VEsXeocdOFfn7VgRX2JfJIRIWR8iFbGDAiTvPPnVWJZnRCck0nF9FjdMGrV+Nkl82DUqfGry4WJotPrx9ef3JpwLG/s7UDdHa/g/N99KOfzm3rsMGrVNGNPZF5NIZbXleCfnzThuS0t+On6fXB4fPjRhXMSP5gQMi7Ulprw56uXIE8n1IkYtWq4fAFwHrvksc/uQWl+8BNAXZkJly6qAQBYxHQM51w+nCfcM5ubAQD72gfx6ZFe+AMcu1sHUFdmGrO1vqwN7CoVwzfXNCDAgf/99y68sqsdX1xVF/M4O0JI7vnHl5fhve+envT9jTo1/AEet5a9xx7ZL+rHFwkTxhMWJ7qtbjz6SRPm/OiNqKkZZUuS+944gHUPbcTmpj6sqh+705uyNrADwIr6YJMhnUaFO8+flcbREELG2hkzKlA3hA2FRjFlE2u2DQB9dndEOrfQqJU3Mb6+px2v7GoHAPzslX34xav7Q9I7FocHDRX5KDBo0NTrwGZxz82tZ01LepwjlVV17OH0GjVe+/Yp0GtUKDXpR3xYLiEkt0kplo5BFyxOLzy+AGZUmkPu02vzYFZlQcRjH7hyIV7b8zp2tAzAJv5ieHOfcCJTfZkJVy0XNjH1OzwoztPiSyfV4a4X9wAALphfhYIxPNA+qwM7AMyqivwGEEJINFJZ5Hm/+1C+runeC+TLnHP02j0oyY+ssjHq1Fg8uQjtA86Q9iaAkE+XWBxeTCrJw+eXTMTb+zuhYgw3nVKf6n9KXFkf2AkhJFmlpsiqlB6bWz4z1eb2weMLoCzK/QCgNF+PDeIs/ZtnNMhVeYc6g7Xv/Q4P5k8shEGrxj++vDzV/4SkZHWOnRBChiLaRqbGLhsA4Ei3DTc/tT3m/QCgOC+YTqkrM+HvX1yK2VUF6BFbBnDO0W/3Rq2rH0sU2Akh40aFWY9LFlajvjy44Cq1AP/dW4fl1t3RUjGAcIaqpMSkxZmzJmBpXTEau2zotbnh8Pjh8QdQnOZT2CiwE0LGDZWK4bdXLcI73zkdh352HhgDWsR8+dbj/fL9pH5U4aQukfdfsQBnzKgAAOjFLpPLfv6WvCFJObNPB8qxE0LGJZ1GhaoCA1r7HXj3YBdOWIILorUxAvu31jRg3fLJqFQcxuMQSx0DPNjWpCjNqRgK7ISQcWticR6e33ZCDshTy02oKc6LeSaDRq0KCeoA8O0zp8Ef4Hjmsxa8tkeob68co54wsVBgJ4SMWxNLjNjcFPz66RtXomKIQbmiwIB7L5+PjUd7caDDCq2aYWaVOfEDRxHl2Akh41a9uGt1bk0BHrxmyZCDulJDhRDMl9WVQK+J3pRsrFBgJ4SMW6eLC6BXLZuMc+eO7GQjo04I5uuWj+wYvVSgVAwhZNyaW1OID//vjIjzUofju2dPx8Rio9xTJp0osBNCxrVJMSpghqq21IT/d+7MlDzXSFEqhhBCcgwFdkIIyTEU2AkhJMdQYCeEkBxDgZ0QQnIMBXZCCMkxFNgJISTHUGAnhJAcwzjnY/+ijHUDOD7Mh5cB6EnhcFKJxjY8NLbhobENTzaPrZZzXp7oSdIS2EeCMbaFc7403eOIhsY2PDS24aGxDc94GBulYgghJMdQYCeEkByTjYH9oXQPIA4a2/DQ2IaHxjY8OT+2rMuxE0IIiS8bZ+yEEELiyKrAzhg7lzF2kDHWyBi7Iw2v/whjrIsxtkdxXQljbANj7LD4d7F4PWOM/V4c6y7G2OJRHtskxti7jLF9jLG9jLFvZ8r4GGMGxthmxthOcWw/Ea+fwhjbJI7hX4wxnXi9Xvy6Uby9brTGJr6emjG2nTG2PpPGJb5mE2NsN2NsB2Nsi3hd2r+n4usVMcb+zRg7wBjbzxhblQljY4zNEP+/pD+DjLFbM2Fs4uvdJv4c7GGMPS3+fKT2Pcc5z4o/ANQAjgCoB6ADsBPA7DEew6kAFgPYo7juVwDuEC/fAeCX4uXzAbwGgAFYCWDTKI+tCsBi8bIZwCEAszNhfOJr5IuXtQA2ia/5LICrxOsfBPB18fI3ADwoXr4KwL9G+f/udgBPAVgvfp0R4xJfpwlAWdh1af+eiq/3TwA3iJd1AIoyZWyKMaoBdACozYSxAagBcAyAUfFe+1Kq33Oj/h+bwv+QVQDeUHx9J4A70zCOOoQG9oMAqsTLVQAOipf/CmBdtPuN0ThfAnBWpo0PQB6AbQBWQNiIoQn//gJ4A8Aq8bJGvB8bpfFMBPA2gDUA1os/3Gkfl2J8TYgM7Gn/ngIoFAMUy7SxhY3nbAAfZ8rYIAT2FgAl4ntoPYBzUv2ey6ZUjPQfImkVr0u3CZzzdvFyB4AJ4uW0jVf8uLYIwsw4I8Ynpjt2AOgCsAHCpy8L59wX5fXlsYm3DwAoHaWh/RbA/wEIiF+XZsi4JBzAm4yxrYyxm8TrMuF7OgVAN4B/iGmshxljpgwZm9JVAJ4WL6d9bJzzEwB+DaAZQDuE99BWpPg9l02BPeNx4ddqWsuMGGP5AP4D4FbO+aDytnSOj3Pu55wvhDBDXg4g7YdDMsY+B6CLc7413WOJ42TO+WIA5wG4mTF2qvLGNH5PNRDSkn/hnC8CYIeQ3siEsQEAxDz1RQCeC78tXWMT8/oXQ/jFWA3ABODcVL9ONgX2EwAmKb6eKF6Xbp2MsSoAEP/uEq8f8/EyxrQQgvqTnPPnM218AMA5twB4F8LHzSLGmHSguvL15bGJtxcC6B2F4awGcBFjrAnAMxDSMb/LgHHJxBkeOOddAF6A8EsxE76nrQBaOeebxK//DSHQZ8LYJOcB2MY57xS/zoSxrQVwjHPezTn3Angewvswpe+5bArsnwGYJq4e6yB8xHo5zWMChDF8Ubz8RQi5ben668QV95UABhQfA1OOMcYA/B3Afs75A5k0PsZYOWOsSLxshJD73w8hwH8+xtikMX8ewDviDCulOOd3cs4ncs7rILyf3uGcX53ucUkYYybGmFm6DCFfvAcZ8D3lnHcAaGGMzRCvOhPAvkwYm8I6BNMw0hjSPbZmACsZY3niz6z0/5ba99xoL16keOHhfAjVHkcAfD8Nr/80hLyYF8KM5XoI+a63ARwG8BaAEvG+DMCfxLHuBrB0lMd2MoSPlrsA7BD/nJ8J4wMwH8B2cWx7APxQvL4ewGYAjRA+LuvF6w3i143i7fVj8L09HcGqmIwYlziOneKfvdJ7PhO+p+LrLQSwRfy+vgigOIPGZoIwsy1UXJcpY/sJgAPiz8LjAPSpfs/RzlNCCMkx2ZSKIYQQkgQK7IQQkmMosBNCSI6hwE4IITmGAjshhOQYCuyEEJJjKLATQkiOocBOCCE55v8DSt4H0QEqU1MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ts);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Level ISSM\n",
    "\n",
    "The simplest possible ISSM maintains a level component only. Abusing the notation and let $l_t$ denote *level*, the level ISSM can be written as\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "l_t = \\delta l_{t-1} + \\alpha \\epsilon_t.\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Or in ISSM terminology, \n",
    "$$\n",
    "  a_{t} = [ \\delta ],\\quad F_{t} = [ \\delta ],\\quad g_{t} = [ \\alpha ],\\quad \\alpha>0.\n",
    "$$\n",
    "\n",
    "The level $l_t \\in \\mathbb{R}$ evolves over time by adding a random innovation $\\alpha \\epsilon_t \\sim \\mathcal{N}(0,\\alpha^2)$ to the previous level, so that $\\alpha$ specifies the amount of level drift over time. At time $t$, the previous level $l_{t-1}$ is used in the prediction $z_t$ and then the level is updated. \n",
    "The damping factor $\\delta \\in (0, 1]$ allows the ``damping'' of the level.\n",
    "The initial state prior $P(l_0)$ is given by $l_0 \\sim N(\\mu_0, \\sigma_0^2)$. For Level-ISSM, we learn the parameters $\\alpha>0$, $\\mu_0$, $\\sigma_0>0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will fix the parameters for the illustration of filtering.\n",
    "Learning of the parameters will be discussed in another notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "latent_dim = 1\n",
    "T          = len(ts)\n",
    "\n",
    "# Set the coefficients of the ISSM\n",
    "delta      = 1.0\n",
    "F          = delta * nd.ones((1, 1, T))\n",
    "a          = delta * nd.ones((1, T))\n",
    "\n",
    "# Set the parameters of the ISSM\n",
    "alpha      = 0.5\n",
    "g          = alpha * nd.ones((1, T))\n",
    "\n",
    "m_prior    = nd.zeros((latent_dim, 1))\n",
    "S_prior    = nd.zeros((latent_dim, latent_dim))\n",
    "sigma      = 0.5 * nd.ones((T, 1))\n",
    "b          = nd.zeros((T, 1))\n",
    "z          = nd.array(ts).reshape((T, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu_seq, S_seq, _ = ISSM_filter(z, b, F, a, g, sigma, m_prior, S_prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the filtered mean and variance of observations\n",
    "\n",
    "\n",
    "Given $p(l_{t-1}|z_{1:t})=\\mathcal{N}(\\mu_t, S_t)$, we can compute the distribution of the reconstructed observations \n",
    "\n",
    "$$\n",
    "p(\\widehat{z_t}) = \\mathcal{N}(a_t^T\\mu_t, a_t^TS_ta_t + \\sigma_t^2).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "\n",
    "def reconstruct(mu_seq, S_seq):\n",
    "    a_np = a.asnumpy()\n",
    "    T = len(mu_seq)\n",
    "    sigma_np = sigma.asnumpy()\n",
    "    \n",
    "    v_filtered_mean = np.array([a_np[:, t].dot(mu_t.asnumpy()) \n",
    "                                for t, mu_t in enumerate(mu_seq)]\n",
    "                              ).reshape(T, )\n",
    "    \n",
    "    v_filtered_std = np.sqrt(np.array([a_np[:, t].dot(S_t.asnumpy()).dot(a_np[:, t]) + \n",
    "                                       np.square(sigma_np[t]) \n",
    "                                       for t, S_t in enumerate(S_seq)]).reshape((T,)))\n",
    "    \n",
    "    return v_filtered_mean, v_filtered_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "reconst_mean, reconst_std = reconstruct(mu_seq, S_seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forecast\n",
    "\n",
    "One advantage of the ISSM model is that one can obtain the complete probability distribution for each of the future time steps:\n",
    "\n",
    "$$\n",
    "p(\\widehat{z_{T+t}}) = \\mathcal{N}(a_{T+t}^T\\mu_{T+t}, a_{T+t}^TS_{T+t}a_{T+t} + \\sigma_{T+t}^2),\\quad t > 0 \\\\\n",
    "p(l_{T+t}) = \\mathcal{N}(F\\mu_{T+t-1}, FS_{T+t-1}F^T + g_{T+t} g_{T+t}^T)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def forecast(mu_last_state, S_last_state, F, a, g, sigma, horizon):\n",
    "    \n",
    "    forecasts_mean = []\n",
    "    forecasts_std = []\n",
    "    \n",
    "    mu_last_state = mu_last_state.asnumpy()\n",
    "    S_last_state = S_last_state.asnumpy()\n",
    "    F = F.asnumpy()\n",
    "    a = a.asnumpy()\n",
    "    g = g.asnumpy()\n",
    "    sigma = sigma.asnumpy()\n",
    "    \n",
    "    for t in range(horizon):\n",
    "        a_t = a[:, t]\n",
    "        forecast_mean = a_t.dot(mu_last_state)[0]\n",
    "        forecast_std = a_t.dot(S_last_state).dot(a_t) + np.square(sigma[t])[0]\n",
    "        \n",
    "        forecasts_mean.append(forecast_mean)\n",
    "        forecasts_std.append(forecast_std)\n",
    "                    \n",
    "        mu_last_state = F[:, :, t].dot(mu_last_state)\n",
    "        S_last_state = F[:, :, t].dot(S_last_state).dot(F[:, :, t].T)\n",
    "        \n",
    "    return np.array(forecasts_mean), np.array(forecasts_std)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let us use the same cofficients (constant over time) for the future as well\n",
    "forecasts_mean, forecasts_std = forecast(mu_seq[-1], \n",
    "                                          S_seq[-1], \n",
    "                                          F, a, g, sigma, horizon=13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the reconstruction as well as the forecasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_reconstruction_forecasts(v_filtered_mean, v_filtered_std, forecasts_mean, forecasts_std):\n",
    "\n",
    "    plt.plot(ts, color=\"r\")\n",
    "    plt.plot(v_filtered_mean, color=\"b\")\n",
    "    T = len(v_filtered_mean)\n",
    "    x = np.arange(T)\n",
    "    plt.fill_between(x, v_filtered_mean-v_filtered_std, \n",
    "                     v_filtered_mean+v_filtered_std, \n",
    "                     facecolor=\"blue\", alpha=0.2)\n",
    "    \n",
    "    plt.plot(np.arange(T, T+len(forecasts_mean)), forecasts_mean, color=\"g\")\n",
    "    plt.fill_between(np.arange(T, T+len(forecasts_mean)), forecasts_mean-forecasts_std, \n",
    "                     forecasts_mean+forecasts_std, \n",
    "                     facecolor=\"green\", alpha=0.2)\n",
    "    \n",
    "    plt.legend([\"data\", \"reconstruction\", \"forecasts\"]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VFX6xz93ZpJMek8oAUIHgSR0FLtiWyv2urhrW11/26y77uqq29RV19117eK6dl2xi2VBARUBBZQSIBBCei8zmT7398fJzZ3J1JCBFM7neXjIzK0zc+/3vud73vMeRVVVJBKJRDJ0MPT3CUgkEokktkhhl0gkkiGGFHaJRCIZYkhhl0gkkiGGFHaJRCIZYkhhl0gkkiGGFHaJRCIZYkhhl0gkkiGGFHaJRCIZYpj646A5OTlqYWFhfxxaIpFIBi0bNmxoVFU1N9J6/SLshYWFrF+/vj8OLZFIJIMWRVH2RrOetGIkEolkiBEzYVcUxagoyreKorwbq31KJBKJpPfEMmL/GbAthvuTSCQSyX4QE49dUZQC4AfAH4BfxmKfEokkdrhcLiorK7Hb7f19KpIoMJvNFBQUEBcXt1/bx6rz9GHgFiA1RvuTSCQxpLKyktTUVAoLC1EUpb9PRxIGVVVpamqisrKSsWPH7tc++mzFKIpyOlCvquqGCOtdoyjKekVR1jc0NPT1sBKJpBfY7Xays7OlqA8CFEUhOzu7T62rWHjsC4EzFUUpB14GjlcU5T89V1JV9QlVVeeoqjonNzdiGqZEIokxUtQHD339rfos7Kqq3q6qaoGqqoXARcD/VFW9rK/7lUgkEsn+IfPYJYMW2Q84uLnrrrt44IEHQi5ftmwZW7duPYhnNHSIqbCrqrpSVdXTY7lPiSQUbW39fQaSA4kU9v1HRuySQYnHI/5JBhd/+MMfmDRpEkceeSSlpaUAPPnkk8ydO5fi4mLOPfdcOjs7+eKLL3j77be5+eabKSkpoaysLOh6kuD0S60YiaSvuN2gqv19FoOUn/8cNm6M7T5LSuDhh8OusmHDBl5++WU2btyI2+1m1qxZzJ49m8WLF3P11VcDcMcdd/D0009z4403cuaZZ3L66adz3nnnAZCRkRF0PUkgUtglgxK3G7ze/j4LSW9YtWoV55xzDklJSQCceeaZAHz//ffccccdtLa2YrFYOPnkk4NuH+16EinskkGKFPY+ECGyPtgsWbKEZcuWUVxczNKlS1m5cmWf1pNIj10ySPF4pBUz2Dj66KNZtmwZNpuNjo4O3nnnHQA6OjoYPnw4LpeLF154oXv91NRUOjo6ul+HWk8SiBR2yaDE45ER+2Bj1qxZXHjhhRQXF3Pqqacyd+5cAO655x7mz5/PwoULmTJlSvf6F110Effffz8zZ86krKws5HqSQBS1H8KeOXPmqHKiDUlfqK0FqxXGjgWDDE8ism3bNqZOndrfpyHpBcF+M0VRNqiqOifStvKWkAxKtIi9pUWmPUokPZHCLhmUaMLudILL1d9nI5EMLKSwSwYlWuep2y3EXSKR6EhhlwxKtJGnbjc4HGCx9PcZSSQDBynskkGJJuhuN7S3Q2trf5+RRDJwkAOUJIMSraRAZycYjTKnXSLxRUbskkGHzSaEHfTI3WaT4i6JzNKlS6muro7Z/h5++GG/YmSnnXYarQOg+SiFXTLo6OgITHF0OmXa42BCVVW8/TDCLJywe/bjAuop7O+//z4ZGRn7fX6xQgq7ZNDR0YEIz3sIgxT2gU15eTmTJ0/miiuuYPr06Tz//PMcfvjhzJo1i/PPPx9LVw/4unXrOOKIIyguLmbevHl0dHRgt9u58sormTFjBjNnzmTFihWAEOrFixdzyimnMHHiRG655RZAiPSSJUuYPn06M2bM4KGHHuL1119n/fr1XHrppZSUlGCz2SgsLOTWW29l1qxZvPbaaxx77LFogycbGxspLCzs3t9NN93E9OnTKSoq4u9//zuPPPII1dXVHHfccRx33HEAFBYW0tjYCMCDDz7I9OnTmT59Og931ecpLy9n6tSpXH311UybNo2TTjoJm80W8+9aeuySQYfTCRN/cgKJe7axeXlN9/tS2KOjn6r2ArBz506ee+45JkyYwOLFi/nkk09ITk7mL3/5Cw8++CC33XYbF154Ia+88gpz586lvb2dxMRE/va3v6EoCt999x3bt2/npJNOYseOHQBs3LiRb7/9loSEBCZPnsyNN95IfX09VVVVfP/99wC0traSkZHBP/7xDx544AHmzNEHb2ZnZ/PNN98A8NhjjwU97yeeeILy8nI2btyIyWSiubmZrKwsHnzwQVasWEFOTo7f+hs2bODZZ59l7dq1qKrK/PnzOeaYY8jMzGTnzp289NJLPPnkk1xwwQW88cYbXHZZbGcTlRG7ZNDhdkPa+hXENdXqZjtS2AcDY8aMYcGCBXz11Vds3bqVhQsXUlJSwnPPPcfevXspLS1l+PDh3XVk0tLSMJlMrF69ulv8pkyZwpgxY7qF/YQTTiA9PR2z2cxhhx3G3r17GTduHLt37+bGG2/kww8/JC0tLeQ5XXjhhRHP+5NPPuHaa6/FZBKxcFZWVtj1V69ezTnnnENycjIpKSksXryYVatWATB27FhKSkoAmD17NuXl5RGP31tkxC4ZdLjdUEs+D/JLLt6zC+9EURAqkrB3dkJ8PJgO8au+P6v2JicnA8JjX7RoES+99JLf8u+++67X+0xISOj+22g04na7yczMZNOmTSxfvpzHHnuMV199lWeeeSbsOQGYTKZu799+gCbV7Xm+B8KKkRG7ZFChlRI4nv9xP7ewc4eeChOpL85qlaNUBwoLFixgzZo17Nq1CwCr1cqOHTuYPHkyNTU1rFu3DhClet1uN0cddVR3qd4dO3ZQUVHB5MmTQ+6/sbERr9fLueeey7333ttttfQsBdyTwsJCNmzYAMDrr7/e/f6iRYt4/PHHcXe1EJubm8Pu76ijjmLZsmV0dnZitVp58803Oeqoo6L+fvqKFHbJoMLtBoPNyjYOA8Dj0MP0SBG7y+Xn3Ej6kdzcXJYuXcrFF19MUVERhx9+ONu3byc+Pp5XXnmFG2+8keLiYhYtWoTdbuf666/H6/UyY8YMLrzwQpYuXeoX+fakqqqKY489lpKSEi677DL+9Kc/AWKyjuuuu66787QnN910E//617+YOXNmdycowFVXXcXo0aMpKiqiuLiYF198EYBrrrmGU045pbvzVGPWrFksWbKEefPmMX/+fK666ipmzpwZi68uKmTZXsmgwmKB3Z9XUvyDAgCW/nIz0y8pAqCgAPLzQ2+7dy+kpkIEe3RIIsv2Dj5k2V7JIYPFIiJ2DadDD0wiRexut6wEKTk0kMIuGVRYrWC06RW/HHbdWI/ksXs8UtglhwZ9FnZFUcyKonytKMomRVG2KIry+1icmEQSDJcLvB26sKcs+w+mpjogcsTu8cjOU8mhQSwidgdwvKqqxUAJcIqiKAtisF+JJAC3GyzNujp7m1oovOuH4u8oInYp7JJDgT4LuyrQQqi4rn+yHJPkgOBygbVF91P2MBbr93sAKewSiUZMPHZFUYyKomwE6oGPVVVdG2SdaxRFWa8oyvqGhoZYHFZyiNHYKMS7o1X3XO7ltyzs+JDkzV+iesPHE5rHLqtASoY6MRF2VVU9qqqWAAXAPEVRpgdZ5wlVVeeoqjonNzc3FoeVHGJoAwE72vzN9HLGMuVHR5D85n9CbququqDLDtT+45FHHmHq1Klceuml/X0qbNy4kffff7+/T+OAENOsGFVVW4EVwCmx3K9EArogd7QHX56w9ZuQ2/oOTJLC3n88+uijfPzxx92jSMPhPsCjyaSwh0FRlFxFUTK6/k4EFgHb+7pfiQT8BdnthsQdm1A++jBgvSI2UVabHPC+hq+3Lkef9g/XXXcdu3fv5tRTT+Wvf/0rZ599NkVFRSxYsIDNmzcDcNddd3H55ZezcOFCLr/8cjweDzfffDNz586lqKiIxx9/vHt/f/nLX5gxYwbFxcXcdtttADz55JPMnTuX4uJizj333O5a6a+99hrTp0+nuLiYo48+GqfTye9+9zteeeUVSkpKeOWVV/jss88oKSmhpKSEmTNnhi09MNCJRTmk4cBziqIYEQ+KV1VVfTcG+5VIcLn0ol1uN2R89hYtZAas9x1F3LxeYWWY/Wgc6lUgf/7hz9lYG9u6vSXDSnj4lPDVxR577DE+/PBDVqxYwe9//3tmzpzJsmXL+N///scVV1zBxq5awlu3bmX16tUkJibyxBNPkJ6ezrp163A4HCxcuJCTTjqJ7du389Zbb7F27VqSkpK6a7csXryYq6++GoA77riDp59+mhtvvJG7776b5cuXM3LkSFpbW4mPj+fuu+9m/fr1/OMf/wDgjDPO4J///CcLFy7EYrFgNptj+h0dTPos7KqqbgYOXhEEySGF0wmJieJvjwfwuNnLGJKw0ol/hG53hb6cfYVdRuz9z+rVq3njjTcAOP7442lqaqK9XXhsZ555JoldP/pHH33E5s2buwtytbW1sXPnTj755BOuvPJKkpKSAL2M7vfff88dd9xBa2srFouFk08+GYCFCxeyZMkSLrjgAhYvXhz0nBYuXMgvf/lLLr30UhYvXkxBQcGB+wIOMId4AVPJQMbj8Y+uvV7YUDOcx7kOgFFUsI/R3csTDKHNc18r5lCP2CNF1v2NbxldVVX5+9//3i3QGsuXLw+67ZIlS1i2bBnFxcUsXbqUlStXAqK1sHbtWt577z1mz57dXcHRl9tuu40f/OAHvP/++yxcuJDly5czZcqU2H2wg4gsKSAZsHi9/iLs8cC3dSO7X29gNmuZ1/3a5pYR+2DBtwzvypUrycnJCToZxsknn8y//vUvXF0/4I4dO7BarSxatIhnn32220PXrJiOjg6GDx+Oy+Xy66AtKytj/vz53H333eTm5rJv376AkrtlZWXMmDGDW2+9lblz57J9++DtKpQRu2TAEixiN6v65Ae5NJJLI+uZzQK+osOdGHJfUtgHFnfddRc/+tGPKCoqIikpieeeey7oeldddRXl5eXMmjULVVXJzc1l2bJlnHLKKWzcuJE5c+YQHx/Paaedxh//+Efuuece5s+fT25uLvPnz+8W7ptvvpmdO3eiqionnHACxcXFjB49mj//+c+UlJRw++23s3r1alasWIHBYGDatGmceuqpB/MriSmybK9kwGK1QmsrjBwpRP3bb+H1i1/nzzvP4wJe4RUuAqD+/Ou57bXZfGQ6jUrXsID9qCps2exhxK2XU3/+DShHLiTMHA1DElm2d/Ahy/ZKhiRer14mQIvcmy3xpNPaLeoA+279J/Gjh2HxJAXdz+7dYNjyHVnLX2Lcby7uHugk89klQxUp7JIBi9erjxbVBL6pM4k86tn58Ht+6yYneGhXUwI6RpubRdSftG0DT/MjPnAej9stOlPb2g7Ch5BI+gHpsUsGLKoaGLE3OlLISvfSfuRp7P3NE1inzgYgJdGDioGODsjIEOva7bBH1AfDUm/lKp6GFlgPNDUdetkxqqqiKEp/n4YkCvpqkcuIXTJg8bVinE4wtTTQZE8mK0nMVen84dXYpswCIDlJrNjaqm/v20la3+Qfw1RXH1pWjNlspqmpqc+CITnwqKpKU1NTnwZIyYhdMmDxjdhdLsh7+RHq+SlzbWtQlJkkJUHXmBaSu+z11mYPFBoBf2Gva43r/ttuB7P50MqOKSgooLKyEllZdXBgNpv7NEBKCrtkwOLrsbtc0GlKo5EczIfPxGgEg097MylFWAzt9XboGpHqa7XUt+kz2rc0uBg+Ku6QEva4uDjGjh3b36chOUhIK0YyYPG1YlwuWFM2DC9Gpv2gEJPJX9iTu4S9rc7R/Z6vcLda4vW/qzsDlkskQwkp7JIBi1/nqctL/YqtAEyfrmAygW8/YFKqaHy2NwYXdovd2P13U5VY51DrPJUcOkhhlwxYlLJdJK5dCUDCx+9S78kiATvJyQRE7EnpQtg7mnU19xN2h+6x11eLBVLYJUMV6bFLBiz5R04UfyxRMZWVUkc++dShKGMwGv0j9sRUIdwdLaGFfSSVNJBLfY2u6B4PGPVgXiIZEsiIXTIosO5r5nmuwEIKQEDnqSlVpIbZOkIIuzOB9DgbBVTSUKnbNTJqlwxFpLBLBgWl9WLU0fjEakAX9viuPtG4FCHsdktwYe9wxJGcolCgVFNXr78vhV0yFJHCLhnwqCo0dgjh/tmTMwC6rZiMDOG3G1NFInunxdMt6N3C7nbR4TKTlAQjklqp6Ujt3rcUdslQRAq7ZGCiqtzAPxhONS6nSoNFlOTNzhaLtYjdZBJRu5qcTDwOnJ0eysqEYGvCHtdYQz15ZGaqDEu1UO3I6s6Pl8IuGYpIYZcMTFwuHuUGahmO0+KkwSYGHWV2TXeqRewmE8TFAQlmEnDgdILFImrBaJgaaqlmBLl5CplpHlxqHF3zM+BwIJEMOaSwSwYkHrteyMXe7qTVnki6oaN7YmujUYi6yQRJSZCRF48ZO86uzaqq9H111HTgJIGcfCOpGSIFRitF0Nqqj26VSIYKUtglA5LmWn2S0tZ6J1ZXHCkmW/d7JhMkJIhoPTUVkjLiScCBwylyILWBTQBN1SIszxkRT2pmV1pkk9h/Rwfs2nWgP41EcnCReeySAYlvxN5Q7aLTFU+yUQi05q0DJCcLn1zxaBF7YKzSWCuM9OxRSZj3ieWW+k5ApNRYrQfwg0gk/YCM2CUDEq9DF/bqfW6s7gSS4kSUnZTU5auj++yJKUbhsbt9Ri253eDxUN9V0DB7TAqpWeKJYK1px9Qs8h49Hv8IXyIZ7EhhlwxIPDbdirG0ebB4EknuEva0NP9RpwBGk4JZceB065f07AVxFJ1WgOWzDQBk58eRniui9Li//5Xik/K71z2UarNLhj59tmIURRkF/BvIB1TgCVVV/9bX/UoObXytGJvFi8WTSG68yF9MSQm+TbziwuHyrw+wsmk6v+cusTwekoYnkEgnL7nPx4iDuV01BVwu4dlLJEOBWETsbuBXqqoeBiwAblAU5bAY7FdyCOPs1IXd1daJ1WsmKVlcrqYQ4YjZ4MTh7hJ2j4c68jibZQDckdYVa6SmMYFdrOFIfsJjeJrElEuyhK9kKNFnYVdVtUZV1W+6/u4AtgEj+7pfyaGNw6c0gFpdg4UUzBli9Gk4Ybd7xEJjRyun8CGdJHMO/+X8l88DwJOSzjS2dG9TsaUDkMIuGVrE1GNXFKUQmAmsDbLsGkVR1iuKsl5OzyWJhN2qDwl11zXSTBbJuaJsQChhTzHasLqEh27qaKGM8QD84r4RuPJGkpEB6phCSpJ2dG9TWSYybaSwS4YSMRN2RVFSgDeAn6uq2t5zuaqqT6iqOkdV1Tm5ubmxOqxkiOLo1IW9oR4cmMkpMAdMsOFLqtGGxd01AXBLMxZSuOGkHSQdvwAQHntcoolzD9vGZTwPQHuDHZDCLhlaxETYFUWJQ4j6C6qq/jcW+5Qc2jisutLubRVFu7JGJYWM1gFSTDasXcJurbOgYiA9W98gLk5E+/bb7+bv4x8CoE1Y7FLYJUOKPgu7oigK8DSwTVXVB/t+SpJDHa8XHHY9sbzCmgNAzoiE7jK9wUiNc2DxJOL1QkedGKWalqunumjC7hgzib3/fJ9U2mlpE+G/3X4APohE0k/EImJfCFwOHK8oysauf6fFYL+SQxSvF1w23YqpcuUBokRvWlro7VLihTrbbNDeILzz1LzE7uXdBcMAb3IqOTTS2iGyaKIRdjmISTJY6HMeu6qqq4EQrqdE0ns8HnA69Mpc1YwAhKibzaG3S40XA5hsNmjvmvs0bbie9G406hNzeM1JZNFMmzW1+5jhpslzu8UgpsTEwGVOJ2FbEhLJwUaOPJUMGFRVlNx1u8Hp0MNjJ8JOSUvTI+5gpCYIYbdYoLVFbJ+Woyuur7CjKKQZrVjs+g7D+ewOhxBw33PVqKqS0bxkYCGFXTJgqKyEtjYRcTvt/rV0kw227jK9ochLFtW8GhuhrV00ItPT9eV+wg6kmmxYHLqwhysr4HT6L7dYhNi7XOKcZV13yUBCCrtkQNDRAfX1QiTb2/2tGICMOAtGY/iIfcxIobzVlR5aLSZMuEhO1pf7zpEKkBZnp92pezvhIna3WzwwNCwWcc4VFcLCkZ2vkoGEFHbJgKCtTfxvs0FLCyhlOwFIQkThGfGd5OWFzmEHGFUYhxE3n977Je0VrWSb2v3WNxr1KfUAUuKddLh10zycsHs80NkpLBi3G5qbxXlqE3bIdEnJQEIKu2RAoE1VB2De9R3Gb9YBkIcorZsebw9rwwCYRuThwcQajuQZfkxugj5OzmDQHwpa1J9idtPuSYpq/lO3W4i6wyEidbtdiLrmrUthlwwkpLBLBgQ2fXIk4uurqGA0Cl4mKmJ6o3SzI6wNA0BBAUfxeffLUSkt3X/7btst7IluVAzdE22E6wDVRL+n167RU9i1h0VtbYRzlkgOAFLYJf2Ow+EvjMaOVvYwlgIqyTaKoaFpic6Iwu79wRm8xVnczH0AjEy3dC/z3VaL/MdnNgNQulUcPBphd7n8s2M0egp7Y6NYr7FRZsxIDj5S2CX9Ts+p6exVTbzOeUxiB2kGIc7pSa6QOeYaBnM8FR+Vcv35jSzgSw6fokfsvp2m2n7mj6ohHgfXXm+i8rOysFaMlvXidkcn7FYrNDWJ7WTGjORgI4Vd0q80NgrP2pd9FV5sJHHErUdhVrtKAyR7ujs9Q6Eo4M7Ko+nW+3jq5RRKbjmpe1mwiD3T08i/uQKAt//dEjay1sTZ5fK3jTR62jM2G9TU+G8rkRwspLBL+pWaGv80QoB2i7gsx05OwOYRoXZOpidyxO5zNdsnzEA16xkvvsKu7SehYgcX8iolfMu+GlPIiN1m062Y1tbgQu0r7KoqttF8dinskoNNn0sKSCQabW3+A4Ii4XQGtzXarOKyzMyE6w2PEe+1c+bCEVFF7KEIFrFX/vyvDH/6XoZtaKPKOjpkxO7bARpKpLWsGUXxF3UI/hklkgOJjNglMUFVdeshWnpaMBqtnSJKz8yEYtMW/sGNGLMzIwp7uOXBhL1z+jzKHnqbYckWahyZIYXdYgn+fk800e8p/r7WTTAbRyKJNVLYJTHB6RQdhuGG5fekZ6epRptN1IbJzATVIHwTb0ZWxP2Fitjj4/2Ld/UsJJafYqXRnRkyGg/XqerLrl3BR6FaLPoAp+bm6PYlkfQFKeySPuP16mLW0y8PR3vAPFuCVkcCiYpdCHCXIa5mRifswcR90iT/GjMJXSXaNbHPSBXKHUp0oxV2h0N8pp4PCFUV3nxNjSw9IDk4SGGX7BdaFApCsJqa9Pejobk5tF9tdcaTYuzEYADHeZcC0Qk7+NsxBgOMHasLue/7SUl6bfeMdP2cehKtqGu0twe3WyorQ3e8SiSxRgq7ZL/wTVN0uUTdFIjeQw5nSXS64kkx2jAYwPbnR9i0vBZjWnLoDXzwFfaEBMgK8TwYPlwX/LRMsVFLbaDq9lbYm5uDR+VanrsUdsnBQAr7IYSqRl4nWux23SP3HZzjcoksktpaEaWGIpS/DmBxJZBscghrJc6EOzs/cjmBLnxTIsPVlsnI0Jdr86J21FgChFx7bWxrJq6hOuLxvd7wI0293v3PkmlrE1G/RBIJKeyHCB6PHlXHArtdj857dpg2NAhLIlRU7naHL5pl8ZhJiXNgMOgReLTC7huxRywapg1UGi5C96ZqBw0N/utowl5yQjZFp46M7iQi4FvwrDdYrVLYJdEhhf0QwOkUQry/gtITh8M/+6OnsDudwmt3uYJbD5HOw+oxkxTnQlF6L+y+EXukAU2asGePSCQZC3vLVaqq/CNqtxsMnRY2UcQqjgzciaqSseJNFGf0Hku41ko4HA7Z+SqJDinshwBtbUJkGxpiU15W89YdDhFBBktx1GyfYCIesuWgqhhbm+jwJJES7wpaajcS0Voxvss9GVlMYTt7KsVBegq7sb2ZEjZxNKu634+rr8LY3kLy+s/45uYXyXjk7uhOkND5+5Gw26WwS6JDCvsQx+MREaLbLfzdWNgxviMxKyrCr+srYl6vOJ9QFk3TL+7liRNfocObTFKKwS9ijxR9a/TGitH26c4exij20dgonkYOB3412o1W/UNo748/bRLjL5rLlq86uIDX+PlHp0Z3gugWVmtr7zpn7Xaxvqz9LomELCkwxHG5hCBoGSDNzZCbu//7c7v97ZVIA5I022HHDkhNFaIWrHPR/OGbnLz6t92vs7Ib/Dz2SCKt0ZuIXWsRuHJHkJml0NQuInZt2rvCwq5SAT5PJ2uLk7zqjeTQyMz6bzntufeBM1jXMSW6E0R8frdbTAXodkNOTuRtnE79e3M4ov8+JIcm8vIY4mj+ula2dn8yMrQaKNB7K8Dp1GcdCmVBxNVWcO8d/jseO9HgZ8VEK2TBSgeEw2gU4pqV7qGpOR2PR+TkGwzic3s80FqnP8l2rG+j+NcnYaeVLzmCTRQD0OZOxusNX9bAF4dD2FTx8dEJu+/4AJsNv7lcJZKexMSKURTlGUVR6hVF+T4W+5PEDpdLiJM2p2hvhvxrtLfr2/U2D9vtjuwpO5Z/xstczJVH7sCI8BmmHJ2/X1aM72CkaIUdIDvVhRcjbW26oDc3i3+WZv1puOqeFbzDGd2vO0nmukmfYlMTezXqtqlJHCPa36PaJ9NS1puRRCJWHvtS4JQY7UsSQzQR0Jrx2mTMvcE3o2Z/Btj0TCHsyfqvhdF8wuXD2UQxH064geHF+RgMuvBGG7Hvr7CnZ4j/21r1ZP+9e4X4tjfrRvh2WyHX82j36wmG3cyeKGojtFT5p7skbv+W2XMUzLu+Cziu9hCIpgWltXp8X0sk4YiJsKuq+jkgyxsNMDye4KLa26jdN81uf1ImI22zvdxMvOJkQnEq3rfeJeeFRwC6PXZfSyYSiYm6TRFNJo0m7ElpYmVbsx4Oax2lljYh7PNYyzrmYSGVm49dx/ucyidzbiFzmHiatFYIgTd2tDLyujOw/Op35FHH1ufXBxxX23dgY68wAAAgAElEQVQ0Iq2Va9DYn1aX5NBCZsUMYazW4B2VvY3YOzuFsGv56bHEYO1ga9MwJqfVYDKBc+TYbrXV+gV8p7WLuD8DTJggIvdo/G5N2BMzxEHsjYG+UXuL+MImFehPqNFHjmLK3ZfS8sfHycgVTYPWOqHScd98xYT1LzG/7h0ayOOhr44IeXytIzUcPfPepbBLInHQhF1RlGsURVmvKMr6hkhtc0lMCBUN9kYYtHKzHR0iVbK3tVPCkf75OxQdk8Fa9yyKxwQOqdRsld4IOwgLZtiw6NZNSur6P70rYm/1/9LSVy5DXbESgPHD9Wg+e3w6zaddhicjm4xhog5wa6NQ6NIvW7CS0r2uzR6+uRHJ3ur5O4b6/aRFI9E4aMKuquoTqqrOUVV1Tm5f8u0kURNKAHoTsTc3C9vA4ej9RBqRqPrvV5jw0EEaJTP9L8XERN1K6a2wQ3SZJiCKhCUlQVKqiLo7O/y/nAk3nUMZ40kzdDB2hK7Aw8boBd7N2SmYsdHcJJpHe7/zr0e815Yf9hyczvClAnr+jlrnri9ut+gT0PpUbDYZ2R/KSCtmCBPK245W2D0e//rq4aL1pG0bGHPXEhR79Ckbr++Z0/339BPz/JZpJXUh+lGn+0NCgshXT+zy2K0dgd5VKZOZZNpNVr7eG+t7fmpKGnnU09IqbqeyPQpxODmDt1nER7R400PWngdhczU0BC81EKqzO9jk2e3tegaSxaJnQkkOPWKV7vgS8CUwWVGUSkVRfhyL/Ur2H7c7dBTo22QPNxK1pSX6ztL2q37Bce/+Ctv6KDNe3W5214qo9/a8p8mZ4i/svhktviJ6IEhIgKSMLivG4i/sVTkzWMNCphTFM268uF1ONXzgt44nJZ086mlqN2Fsb2GbczyH5Tdx3zO5XHi8eDJW7wk9AMBi0fsxehLqIdzzfW1b307uWNUGkgw+YpUVc7GqqsNVVY1TVbVAVdWnY7Ffyf4TrtCU5umqamDGhS/hosxuvF7SP3ubWxz38D0zWLEiuvSVhIqdfOWZy/kzd3Hu+z8OSHvxzVtPSeGAYjCAOS0BI246rf61jR833IAdM2ffNhVjdjoWknlx2K/81vEki4i9uSOe+Oo9bGEaE0Y5sRYdzrCJ4uRrt7UQd9PPWHnSHwMUvL09cESvRm+F3eEQy6zW3g8mi3XHuKT/kFbMECVctKYJiDYqtSdapB5NsaqsD19kwq/OohPRC7l+W3RDImu+rqCZbGbMDu6zHOwh86bURFLpwNLjgVhqLaAwsY7CQvCkZ5NMJ7ZTFvuto8YnMCq+ju8ah/P6L7+gjmGMnSieTAVTUkikk2/WeTly5T3c1Pxrar4N3lnRG2HvGeH7/qZ1dcKa6a2wNzXJDtihghT2IUq40YkulxAMLYWxJ1arqO0S0ov3eknc/g2KvRODzYqVJEqZDMDnewqiEpStK+oAmHpM8I7Fgy3sxtQk0min0+bfcmhzJZEeL75M+9ipbH1xI9XXBVZy/KnrIQDurP8pAIWHiQedaeQwFvExL342kg6Ep7R3W/AfJ9hvEapfo67O30PXtrXb9VaY78xW0aBNSB5uohDJ4EAK+xDEbo8cbXd2CvFX1UBB0aoIhiJ99Xu8ddlr7Lv9URosZlKw0kEaJxs/psWVytpl4dNnEr9by383jGF0Qh1jJpuDrhNtCYFYYUgWEbvV5n9LtHpSSEvQvyDbpOKgCfKTsxp5liXdrwtnCBF3jJrASKr81t1TFvzLjRSxm3d911333ev1t9tCFWbzrcQZCadTPCykJTP4kcI+BNmzJzDaHvbUvaR//k73a6dTF3RNFGw28UCINBHEnq/r+TV/4uxVN7F5l269/PjHKkbc7F1VwYh//gZjW/DByLuXfs7nHMNFVyaGHEQUbTGtWNEt7Hb/J0qbJ5X0pMj+ROnjK1h4dh7xOEilnfyRosmhxsVz9GGiA/WFYx5nFBXs3hc8f9PlCoyWNZE2tjYx7aIixtx7dfcyTYC1kszBsNmiG3ugpbTK6feGBlLYhxhut7+/rjjsJH2/lpzH7sH8y590v+9y+TffVRX27YtgwXSxfqsu5t9ViMj0N9zLxIvnMZlSPl2bzDXPLuCm07ZiaQ9UnFW7hmPEzakXhk53OdgRuzElkTTasTXbMe/eCoDictJKOqmJkZXRUTiFzp/8ivN5jaNY5dcXfOzcTqwkMef8sUxVtrGrNnQ/RM+xe5qwx1Xs4nr+yZPvj+guR6DNUBXOF1fV6GZsam/XUytlmuTgRwr7EKOnBTPutguYumQB5/IGo6jE0yHucqdTFw2rVXix0c7ss6NGT1NZsbsQgHOfOwtPaga/5o9sYTrvcgafOo7kkwc2Bmz/fdMIJqVUk5oa+hgHO2I3JcaRSgfrLVN45eJl4hysHbSQSVpKdMNt3dn53PFgFg894P8wq77299T8+d90zF/E5OQqdrXmhoywKyv9hVh7yNZvruVfXM/t/Jnt75d1L7dYInd4RmOt+E487nttSAYnUtiHGD2b0TtW1XIBr/AepwNQ872wR6xWXRBaWqC8PPpjlLYNY6ZJVCzcZJtEmqED77QZABz9C33Q0Qiq+Hazf+it2G1st49hUk74Xr2DHbHHxUErosTjPZ5fA+Bus+DATGov0i2tR59Gx7Fn+L2nxifQeuJ5oChMzGvH5k0IO4rX9zfUBLayVFf78o16SB3NCNNI1kpHR2AGjbRjBjdS2IcAvtZJzxv0PF7nNS7ofv3FO8Lvtdn07bxevdpgJDI+eIktzonMHFHLKMS8eMPjdS+9/tJf8BA/51oe45j0TayrHeW3b+u7/2M345lQnBTyGAc7Wgcwm+Fm7gdgOKL4uaVeZK+kpEdZWjIKJkwR3vv2tS1kvfd8UANcGz+gqrqttq9cD/H3lOvrRhOxh5q1SqO+Pvh70V4TkoGHFPYhQHV1YEeohhP/jrpnPxrZp2b26r+spp10JhyRzzm8CYDS4yo68955/Pa6BmaPbabenUXNt7Xk3n4V3g8+ZMXzos1/1OWFIY9xsKN1EMJ+Mh9xB/dQRz7mtZ9hbRZPydR0cUKxSMEcfUYRqbSz8/GVfHDnlzT/44WAdbRspaYmXZDL6xPJMHYwS/mW3TX6Q9FqJaoJPkKloGodpsHWj8X8uJL+QQr7IKa2Vtx8jY26teIXAKoqCnrYdTn/poE81nwURW+aD4W/vYzMj14Br5e3rSeQn9jO8T8v4rfcA8D5ye/7rd9yyiXUXvVbSmYLQXz9hv8x5uOnmPfbU7iz6lrmZu2ioDB0AZj+iNi1zs4J7MKLkek3HMOOlcIvSck0ER/vP4mHL70peWCfeThTDaU81XQON/AoF7xwdkDQrqWg+kbSu9tyGJfWwKSkfexsyfZbP5raPzU1Ilsq4HzsoSNzOVhp8CKFfRDT0iJuVlUVzffmHtmFBmsHdvQ88UVLRpKElW3Ltkd9DMVhx/vBckb++gpMzXWsV2czv7AWkwmq39/EnmOv5ML75gTdNv/E6aTRxr9dl/i9f/G8sqDrx8WJaov7U80xVsxAn+3ouVXjAEjNSiAhIfR8qsOG9SKaN8UxOkX/oVq9aWz9IjA0ttn0QWbpn/6Xze6pjM+3MiGziSpHblSZLr60torro2cUHk68pbAPXqSwD2J8o62mpsDBKO66RlrJ7H496uw5jDTU9GpuzrjqcvJo4FQ+QKmoYB+jGDNGLHPljaTpgWexz5gbdFvHhOndg3MuHr2Gl8xL+OiSpRxz53EB6yoKjB8v6sIkJgYsPihUPfYO40famY2Y8WiTXYymTcpLIS5O2DUaKSlipqbUVPHPHHycVVBGpwvvI18RIXnVFxUB63SnPbpdbL11KQ3kMXOSlfHDhOnem85uX3r66eFsOZkZM3gZdMIuhzsLHA7/7yJYje7sC08EII86jmElqSPTGGZqos4S/RT3naXCE1/B8dRva0LFwPBxUSqvotDS9WApOjqdiaueJeuXS1DjAkPyzEwhlJmZ0ddSjzWORaez98mPWc9cnkAMBMo1NjF81nBMJvCdRiA1VYj7iBHidSibJhjTRoiwOT+5AyNuamsCL2qtAzXvlX/wONeSSTNH3DCT0ZPEd1e5Pcrc1B70jPTDReVy0uzBy6ATdtlTL4hUeVGx2/iIkwBYxtms5DhQFIYltlLXGb0pXL1Vz3vbvE70wI2YlhX19imIJOojjjaFnbhUy2k3mfovYjeZRCvk+zdKWXjtdE7kY57+8Zfdk2rHx3dVgjRDdrYQeq3yZG86VhdOrOcwtnDPrGWMpIqahtAbq+V7+JLDOXlKBSnZCaTOnwZAx5bAKD8aVNU/Eg83e5NWKTJUSWHJwOUgl1rqO15v/2RNDDQizS5oam/mfm5m5vBaUv75PFs7xZMgP9lCbXsWqhrdBNHlPnVN3l6dhQEPE4qjj/jf5XS2M4W08c8SbphPUujsx4OG1mnrGDMJrprI347aKGrDoF9zZrP4ZzT6X4ehrkmTKUiJ3cuvYe3mc9j3o4cYtaqS6pbQMyx9unUELWRx5E/FwzRu0jjicdBS6x9qG9tbUA0GvCnpET+n2633F0Salq+tTUTuigIjR4r3bDbxuQ7kBCiSviEj9kFAY6N/SpvLFbmZbK9rZTfjObG4AefoCdimzAIgL9WGTTVHN8pUVdm2Q79EPuNYjjRv6JWfPJkdnMXbeNIyQ66jKP0XpfviJ86Kgm3KzG611yLyuLjgnbs9I/b4ePG5CgoC13Vn51P6zBd0Tp9PQUID1W3JjP3NJRhbAzs/alrElz19unjtzcginzqaWv1v3RE3nE3dsReSsHcHANMWT2bYU/cG/ZxaxN7eHtlHt1jEteZbpqKuLrYDmKSXH3sGnbAfih67xSLquPi+jkRThTBT84f7/8R5meIuamgAPB4MttDpFYbdO3mv5QhOzP62+71Lru5daL3j0U+o+sm9YZsHiYnRtR4ONOHSLHsr7NnZ4nNlZIS3aUYlt7DXMZynlhfQ9McnApY3diZjNjhJ7mokqXHx5BsaaWr3OQm3myu23c6pfMjN55Zh6LRwf8UFbH5sTdBjakIazXXU0SH+aYGE1ytsmd5m5YSirQ2qqiKvJ+kdUtgHAZ2d4nNrzeZobir1wQcByB3pr0K5OaLJ01DnYdRff87Mo1KCJ0KrKrsv+Q3VjOSMCxL5T9wSpvE904/LC1w3DB3zTqD2x78JudxgOPAzJEVLOGHXOkdNpuAdpT3FOzFRpEEajeEzZkbkCZW9lfv42RcXBiyvd6SSa+7we/DlxrXSaPWZTHvfTj7jGAA+5FRs3+/id9zDaXwQdr5UP99cVUne9EVAk9jh0H15rxcqKsT1GMnCiRbtweFLU1Ps9n+oMuiE/VCyYrxeoblatKTdiBGF3evlrbbjMOBh2OyRfoty8oXf0FRpo/LVNdzE/bDl+4DtzVu/4Vee+0mjjZLFYyl++de8+h8n5tG9E/ZIFBSIqHYgEMonVxRdzOPjgwt1T785KUlk+GjbhKJwjH5Bb7GP989S8XppcGWSk+w/HVZuYgcNdp8KarvL/EYYf3LfN91/l28NnEpLO4avnZfx6RtM+fFCst99LuS5arn1qhq7ztTOTv8S0iD2LTNy+sagE/ZDKWLftw+2+4wl0qKYSJMUx1eXs5yTOWlyBdmj/K2T7OFCAJqr7NzHrfyVm/h4mb5DU1Mds+cZqbvzUfZSyC9vdJOcmYBjzKRunz4U+2OnZGURtsrjwSRUxJ6ern+27OzgQu0r7Eajf1QfTtjHTPNvrtTt0X8Lo7WdOvLJSvHvKM1OstPgSiehdDMAtd834sXIr6a+B8A95Zd3r7vri7qAY1qtwl/3FWdjSwMvcAlVH4WejNxm06/BcDXge4Mm4HU+pymFve8MOmHvr4jdtyDTwaKtzb9J2tAQuaATgK10L3spZNq0wC8rPieddFppqrCwl9EA7KvRQ9WE3du4mBc5oVzMR37kmdGnNub1MphPShpYGU6hhN33c4Vax2jUl/XsCA6X426eNY2POZGXuAiAum36qFRjewt15JOd4Z9PlJ3mxIOJ/1wqSjns2yUukuLzJgDgRjxlTLj8qkJqdHbCrl3+7327O53LeIH/+/Ji8l54KOi5NjT4j5WI1OkZ6V71eHQXsLlZJAg0NIiHzsG+14Yag07Y+yNid7n0C+9gEawcq90Oe/dG3rb8c5HjPGFmYCjszsylgEoadraxg0kAVDbpSrR7q52XuRiAkvTdZGRGF4YbDL0X9oHirWuEeshEmwWktTx6Crn2OQ2GQNvJNn4ah0/rYPI5hwFQW6b3aCptLTSQS3aWv0JOyhepq3/hNuytNvbuEyc+bqaeeXRN5quMp4xvtiYGFdie75VXik6Cb5hNy0PPBv18PcU2UsmB5ubw6/h66263aKFWVIh7PJqMHUloBp2w90fE3twshnAfzCgiWClViOyvp618m5ffEwozYV52wHLbhBlMV7bwee0kmhHLK1v1B8DePSIkG041d/6oMmD7UCQkCMuhN/noydGnwx8UQgl7tK2KCROCe/Bms7BqUlKC2E6mOLY/txbDJRcTh5PqCr2301JjwYOJzBz/E5g/Xh/EsOOLJvY1J5Ed1465IIe7uJPR7OXKV0/jSp7l66bxrHurOuK51zbqXlIxm6P6vJGE12oNnXnjdAamTPoGbV6vjNr7wqAT9v6I2LXI4mCOvot2NqOerHl5H2+ymFTaycwO/HnVBDNT81uwIaL0MUoFey1ZGGxWLOf9kK8/EQf++BcfMOriI6M6ZmKiLmbRlgMwGIR3PZAIJuCK0rtqk+npwSP8/HwYPTp09O8ZPpIx7KW6Vj+JlhpxwWXl+/fMJjVV0txVqmHT1w6q7FmMTG0Hg4E7uZtyConLTOEXPEQGLaz+b6DP3pOaZv+OgMaqyBd7JGF3ufzF2TcoKy0NLFrXE1mEbP+JibArinKKoiiliqLsUhTltljsMxQHW9jb2/V61R6PaCoejGPub7rXqlJR0OSFG78KuU7uJN0PODpjE02udOpveYBjy59jqf0iZqeUYr/0x2EVLSdHtxxycvSaKdFG7Hl5A8tfh+Dn09tzzMgILt55eeL7CuW3q+YkJpvKWLV3VLeitdQL5cwY5m/a115+E8Y5M5nOd2zaqLLPW8CIDNGU2/7MF2x5oxSArSvqWRD3Dd+Uh+8nSVv9PuuaxlGIXtf309ciqC6hywVr96jL5d/C1MZiaPO0Rmp9y5TH/afPwq4oihH4J3AqcBhwsaIoh/V1v6E42MJe1yPYaWiITSThcokHhXZxezzCV7dYovPRg5H50Sus6JjDyYWljPjhSSHXyxqjm9vzx4lm/alf/haAW+f+j7v/FHmseHq6f50UTcwSE8XrSEPOB1q0DuI51jOzp7fCnpYWXNi1/WojUoNxSuJntLlT+NP5Il2xqVFc7BkFyX77cBROYedj/+Pw5M18Uz2MPYxleLa4KK1Fh4uSCIAnNYOZhc1s7xwTtgX4+T82sZ2pXJD4DtuZzHh2sXllZGEPFrHv2aN3zDqd+hiMzk594pBoR61GM4BKEpxYROzzgF2qqu5WVdUJvAycFYP9BuVgdqioavCLq6/eX3OzSGPctAm++Qa2bIHNm2HrVtFE7e2DI766HLxeKtbVsYdxzDxvXND1NKHNHakrz+yZ4smiYuCOszZz/r+OJ//w4Nv7YjbrHrmvgBsMMGqUiOC18r7BGAi1YYLR82EU61aFooS2Y845vo0zeJs3q+bS3uqluUkIe+YI8WVlZ4tBTxolh7lo96bQSTLD8oOHv0VF4mmwdU1woR721D28tUsUFjvt6fMYOT2TeXzNxroRwWbt86PnvejxiGu7o0PcI9qgptJS2LZN/N3Y6D9xdjg6Ow+tcSuxJBbCPhLwGfBOZdd7MeeGG2DmzAOx5+BYrcFbCOFybCN547t2iajGtylqt+9/SyS+Zi/qmWeS8/g9rCoVaSnHnBAYKsfH6wNmcgv1Xsv0I6ZxHf/idv7ID26ZHtUxtQE7mkD1HHWpRfOpqcGjU61K4kBk1Cj/czsQdlGoh1rtHf9kyWUevBj59u191FSpJCo20rrmXDWb/bedsVCv0pk/PPiJTjxxNAY8lH5YHrjQ66Xmsbd4hzO5qngdKZNGULr0K04qqqXOlcXqT8L77FrmlmbJ+FozTU36376BULSiDoGVKCXRc9BuL0VRrlEUZb2iKOsbIpUmDIHJFPsOlVARgccTOHGFRrDMlNZW4cXv2CGifLtd2Dja/js7RZQebH7J/SF585dMvH4Ryk9/QjGb+b//LKC8KZVMQ5tf3XAQYhAfr0fYhrwctjOZ9cymc9pc7i9+gavuHUtcQnSXg1bbRRP2ngNwjEaxjsEQvLjXQBmQFIyMjNAzJcWKkNlAisK4xSVk0syGF0v5vHYSM1N3dT8czWb/7zO7WK8wlj8hRO7onLkUxW1n4/eBH8TY3syLXIIBDxc/MLv7/QWXT2IM5Tz6p9aw0+65XPDdd3qtF991Q9ktvY3Atfs91L0oCU4shL0KGOXzuqDrPT9UVX1CVdU5qqrOye2pPFESHx/bJ7jHI8Q2WJOzri60CFssQtx9L96qKt1brKoS6ZGVlVBWJjqNtm2LXeEkAOPbb3LK179nxl4xSOVtx8nsac1gbHJgnmR6uj7tnNkMruxhTGYHo84/HIxGSp9eTcspF0d97Kyuvri4OFHKNVxUG0zEsgOzMAcUvg+qAyHs4R5snlGFlMRt4cXGk9jGYRx/hv4F9ux8dYydwv/xN+axlvHHjAqyN0BRmJu9m89bi7jvB5+x42W93ICpoZYnuIZTp5STnqlLgfW407l+8qeUWYbRsC384A1VFdPttbb6C3usAjCtdVxbKztTe0MshH0dMFFRlLGKosQDFwFvx2C/ASQkxFbYW1tFJB1McMNF1h6P8A3LykSTs6XFPxVSE35tP6Fy0vvCG7uK+JIj/N5b4VjIxDH+d1RcnIjytCg0Jwe8yals/qCKfTc/sl/H9hXm/NClxAG9g1XLc4+PH9gRO4iHnyagB8KKSUgQLZ4RI4JkySgKY9JFxJBrauGM/xunvd29nfbg8aak8atbTLxx3ScYTKFv5XPGi7z0V+uO4W9P6iG/taIJKykUFQVuM3zheACa9kbuwdS89QNhm2ijrz2e8AMEI/UHHGr0WdhVVXUDPwWWA9uAV1VV3dLX/QYjPj52NSpA7xjteUHW1UXuINWalOXl+z//ZF/YVSvU8TKe59Pz/9X9/rFn+c+OpM3VqQl7WtdiV+6I/TK6ExL8o9hI9WGSksQ6U6aIPO7Jk3t9yINORgYUFoq/D0TEriji90hP11s/vozJFhfm6OQmvwk+tO/a92HQcMEN1F0VunomwJRFBWwxTOf32Q+ztm0qm1aJkatNFeI42QWBvbk5I8TTo6EqujDZbj8wwt7Zqd+n9fWhrZyamtgfezATE49dVdX3VVWdpKrqeFVV/xCLfQZDi1RidQFpUXXP/fl2/ERDfwyaKmvNZk76Tn7xxQXk/PAH3MJfuMDwGrPP9J/ZQRN2LXLen1lvfMWtN3N7asfTUiDT08MXxBoopKUJC0mbDu9AkJjoP7DLl7njRQaLMU6/PX3TQ3ubUdR8+g+xfbmJM/94OKm088mzYjBGY6VoZmaPC8w9zR0lTqymwsWov9xAfGVZ2GM4HAdG2F0u3V/3LV3tS01N5BnFDjUGaG5CcLSbrK9em9stmm6afeLrDcayJOmBIm7bJja4i5hSaEONT8A1bDQ/euJwfvPeQgxG/xA6Lk5Eer45572twpiaqgtQb4XdaByYOeuR0KaCyww98VOfyM/3t1V8mTorkV/zB+479r3ucxk+XF/emxmsujEa8cyez5yE7ygrF0/qmhpxIeQVBnaEJORnUMgedm1xsPa1Clp/eodY4HaR+fGrAaOTvN7Y9iH54ns/apaLquoWZ1ubOL4cqaozKIW9L3nkFgtUVwu/TmvW+UYa0YyIO5gYOi2M/c0lJO7U63eUPvQBHaRRtEivumWZdbSwV3qQFmTe6t5G7b7ZGPsjKvvZV97vHMjz1voZggl7y5lLuOruMWRcfQEgInRf16wvrYi8lE7qbGkoDjsNX+/GjI3cvMAnvSc9izms5+3K2ZzJOxxd+RKNO5vJWLmMhtv/ivNPDwRsczBqu2itY5tNJCf49mf5PgC0Ou+HKoNS2PsSGTQ1CVH3nY7LtwUw0HreEzau5YfLL2baxUW0PPJvUFXe/G48ufFtHLl4WPhtE4JHnL31jbUJnGH/BhbJSY9Do7Wo/DAYaD7tMtzZome6ZxZRX77PvGQrte5sEnZvYz1zGMfuoF0tntQMFo3a5vde2f8q2PtlFQtYy6K3forL2ccISFXJeOMpHO9+7Pd2yjefoziD34haxO52iwCstFRf5nvv1tX13lIdSgxKYd/fIvxtbXqk7huV+z7pB5qwf7+mlXc5A4AHXx5BQsVOvnbNZPa45oiRW6gJonsrDElJ4l9y8sCryDjY0TpSw9GzldQnYU9z4FAT2LWuhc84luNPD/2knnzJXL/Xtf9dw7q3RC+llRR2rYtcdiAciaUb+cmfRrPwrkVUrxK5wnF1lUy+5hjG/vayoNv41qHpie993HMC7kONQSns+/uDhUph9PXbB1rzbcsWPZxrNeZg/3oTZUzgsNmRPZFQAtCbiF2L1lNTRQbHQJh0eqgR6QHd83cM1U+Slhb598nJEIr45XpxEcw7NXQ5zoyjRR5kNo2Mo4zlTXO4jb90L9+1PsqiLyFo+XQ9HyNqGn3zkQivDTWVHM+n/OrTU4NuE07YtYBP6yeTwj5I0G6A9vbeb+t2h7dwWlrE/wMlYjdY2sm/6zq2bDcyNr6KH2Uto9RRyPY14gaYelTkWS1CCUZvhF1LxzMaez+RhiQ6EhJETnuo7NNgv2PP35eihvQAACAASURBVDAhQaSTat59Wlpw2ywnSyjjmi2iwueY6aEHFbjzRrCVqWxjKtMSyvia+QCcwCck0kllWd+ioF3b9A7Y9VtEoLJ3WycrOJ5n+VHQvq5wwq7du9oE3D3nUj2UGJTC3tth+TabKLAV7gmu+XEHqme/t7ieeZ7D3r2PN9xnUZJXxZiRbpq9GazcPgwjbqZOjzxyJljHKUTflE9IGLwdn4OJxEQhyMFmlDIYggt+T2FPTBS/l2aVJSUFHwimdZSubitiTFxVeGtNUYh78D7qX17BlEyRc5iMhYdfH8UEdlFRFf4aVBz2sJkIpVUpKHg5K/VT1lSNRfngPWo36eVUa6sD84g1jz1YAKYlPvgu0wK2Q41BKewdHZFzZt1ukdva2ip6ziOt73CIgUbBnvBZ7z2P/Re389HpfztoQ9ze+yqbDtLIo54Lb8xjzGhxgzzdeDbT0iojZqckJITu6PSNAEcEJtJ0k5V1YAboSPxJTxe/VbBU0lDff8/3e6ajav0iPckcoz/tR5sjD4luO/oM7BOmM3GkiIqspOAonMyouDpqWkLnviZu/4ZxC4eR/ejdIdfZ1pjLuMRabpjxOS2eNGb/9gdc9Yle2qJsY2BFPS1iD9Wydjr9MzEP1SJig1LYe87MEoyGBjEpRllZ9BXlgvWiG1sbcd55L0eu+hO/rv0Zdd8d+GpEpsYavt2RxMSEvbz/VSYTTihkTrGTYYiOq5kFkW/IcJGYdvMbDMGjOq0u+UAf+j9UMJvFd94XYde21Yq9paQE319cYQFP8WNO4z1+dpr/jNa5uaE9+plzRHR+RaaoFpKX2kl9Z+gLJLFsC5m08qMXTwy63GBpY7NtIpOHtTLiFxehoEfnc40bUPCy9etAYXe5hHiHEvaeI2B7E4dFM1F8bwhXQO1AM2iFPZzPbrf7j0Try4/lKN3LFPScqm9XxKg8Yxja31/Dck5m/sQWMAnfRJkxgzu4FxMuFi6I/IFCZcSAf3QX7ObPyRE3+UCbbHqo0xthz8jwf/Bq25rN4vfTRhz3FGrnsNH8mGd4j9MpvOk8v2UpKaHHKSQfXoyVJH7zYxFc5KQ6aXSnhxROS60Q5RWOhXi9IkBKWy0K1qWs/ZTOY09jN+MpmWXAPnYqy0b/jEmUcjn/5vElX3Iyy3nzg8CTcThE2etQDk/P1ny04qqqIkUyVlUk3e4DUyMqWgatsDc0BG9meb2wc2fsmmC714iJgKcbRPmbXaUH3opZs8qLkwTO/b1enck2sYjjXrqWnb/4F7N+Mj/iPsLlm2vRoVYcrOfNbzaLuuQyA+bgEk0nqUZWlr+wa9vGxenzzhqNgVabOyOH6qt/x5aXNwf8wFq0H6wPpnP6PMre2kLHRdcCkJvpxovRz8Me/tjvGPnIrQBUVur73rfbSeFvr2D8z0/H1FjD5/d+zpGsAWD+YjF1Q8GLf+Gtpxr5zZNj8V53AyfwKQ3eHNpb/YMYmy18P1h7u39rPlphr6gQtm1NTfQzPIXDYumfUiMag0rYtSjU3tUnE2wy3H37YtcTbmxtZNnrHky4ePid8cxV1rF0/Qzq3l0XmwMEw+1ixZZcRifUMmqM/89jm1hE06X/F1Fxs7Mj2yi+dUp6Ror7NWRd0md6I+zgb7eFyoAKKDKmKNRc+3vsE2YEPX7PloAvzpFju//O7np4NNbpyjniqXsY9u/7MFjaqKjTT2jHF428sXMGJtzsXF7GsiYxSfpvx/2HkZPFwVRzEtaShVhmHgWKQurpxwJQu9PfjvGdTjIY2vSSGtEIu6oKG1ZrfVRU9D07rr/LkgwqYdeiUO1Lq6vz99Ds9tiONut48EledZ7Nlac3kJFv5kfq0wD86899fKSrKq7VX+Go0uuQKg47zpomvvrVa3zkPI6zj93/Y4wYETnajovTxUAT8sxMKCiQFkx/YTIFCnQkYdfqzYT6vbVSyZHQCp6lpQWvONmT7GEirG+p0IW3jTSOYSWnHmvjk+16MbpdX9TzfONpqBh4/9/NvO9cxKVT1nPWq8EHIQEUjBGefsX2vqWpRSPs9h7JOy6XCBD74pH3t7APqpwHLULRvjSXSwh5Xp4Q+LKy2NZ5eXndeBIUBxf8TLRnz7iugA8ee4tVzmNQ1f2zKjyt7Xx47lP8se16APKUfSSrFryKkR3qJOASilL3cP4dU/brnLW655GIj9eb3GazEJCCgsFRfXEoM2oU7N6tX8fhUlONRr0scjiSkyO3YrW0Vt+CceHIHCGazy2VVkDUrbiRv/M5xwDwevswRppqGefewdPrj+7e7tmmMwEoHhM+cBk7JYF0WtmwAU68PPL5hEIrPRDuOwomwm1twq+Ppgic1yuCTN9CbdFue6AYVBG7Juy+JQVqasRk0Nu2xfgp6fWytmkCc7LLu3+g2qvuYP4ChRZvBo3V/m01Y1szxt072HnbU7jqgw+1ti1fybWL9nBn2y8Zz27OyfqMubl7SVBcjM+3ctP8z7lv0cf848XMsJ2f4Yi2kmJcnC4aqakiypei3v9kZPj/hpHGHITyxH0JNZ7BF98JU4zGyHZc5mhhodSt20fHxp2gqnzLLE4atZVrMl4BIM/cwbXJz5OAneONn7GALwEYbazkqJsWhN2/p3A8M/mW8u0O0r5czribz93vVONIkXco26WyUn/AhvPLGxpEYUFtvuOamv4fGDWoIvaeVgyIH60vTSZj6VbUzEy8efrj1tBpIfPo6WxkF9dP2+q3/phpqfAV7F1TSe4FYpaZuG2bWHX5E7zNmXzEVSzevpafjv+AjWoxc+84mZQ4J/vuWcoN/1tMk5LDw9dt44glUzCYDhM79XhiNlVPNDcxCPHQBr5Eu43k4JCernfgRVMPKFLEHulhn5cXaPlkZ/sXyuuJMT+XfGr55zeH8/RVNpbft5YGCpmRWse8WUk88RY4PCZKnv0/dn5yH/VX3Y73ltu48/MW7vpFB8bMC8Oek3P4GEZmfs+KtmSUP9zDTbVXcsG7X5Jz1pHhP0ywfTnDP/xCibDTKTpqU1LEd5GdHdwu0zpzy8vFvTQQio8NKmHXfMD9LQLmh8dD5+XXcNaO+zk/ewXXLD+3e5H76w0s4CsMeDnxhyP9Nht/zEgMT3t4+qEOrJ/8h2OeuIyvn9rMT/ln9zr/rZzPfytF5sphp5aSpFpY7/05mTTzwNWlzLtqZuAH6wPJyfrFFW29dDnwaOCiPWi1ScHDkZoaXX+KyRQ6AArmqefmikg0lOi5MvN4lKv5I79mA3M45hYRgWemV1NYlAZvweT43dgnnNDdUWu+6UZeMP6S8tOXEk3CyMhMGzUtOdzv/BnPcD4Nz37LnWdFsWEP6upg3LjQy8NF1w6HEHaLRWTcKApMmOAv7loGntMZOH3fhuoNAfucPWJ2wHuxZlBZMSCi9lhYLubtG7lixx20kMUTTbqo43ax/dH/Uccw7v5ZEyOL/Gumxk0aRxGbWecq4VffXMZDP9zIM5+NR8HL3Ye9zJtLlpFKOxm0cMvsT9jqmcx672wuy/6AV19wMe+aHqLeBwwGcfMPH64LQG8nwpAMPOLjxW8bblSwRkJCdJODh3pAhKrYaTSGb8l50jJZzJusGekfeWdmwPCpGbzNGfw9/x6/Zc4Rhey+/794U6JrIk4e3oaKgWeazwb+v70zD4+iyhb47/aSXhJI0tkXkrDJIuKIKKCA4KAo7qKOPhecGWVUnNGnPvfP0TdvxnFw33Bc3htHxdFhRBF13HcBh0UW2TdDIBCIEMgC2e7741bRnaS700l30tWZ+/u+/lJd1V33pOvWqXPPPfcceLvsGKrXlLb5nK22mvqTJuH84J2g59m3zx8iHYxwETCHDik3TF2d0jvBskbG2+0SjISz29zu2FjsOz5Zy1YuZaCnjA11hdTUqA7ue+p3PLT5fHxUMvGSvLZfdDh4fOgz/HX1SOrsKbzyvVoCffcZy5hy38XQ2MAC8RgHTj4X55BJnPDQaxyoszPk7gvanitKiotVh01NVZM95opRTeITKpd+ZykuVsp61aqWrupwrp6wbiCbjVVz1lKfW8SXv76eccueBCAtw0Zd/2GMnlpIxYUzopJ57NC92L5uohEnHmqpw8uL1y1ixqdFpCz+jPrcPjRkFyKvmMYJNR9xzUOvc9Wpbc8jpVrb0tjYtjxj69wyrdm1S+WbCZXmG6yZtiDhFHusLPaVK9SVOv3IUjYsLmTvjlqSB3r55xdelvMT7jxtCQ5HcFMo+YUnuK7+ELKmlsIpf2Q5RzP58mHqoMPJwRm3YLr0Cm8O70uMBq/XP8mVlta5Ihgaa5KSEtsCJeZILiur5erKcG201/6hElWZPPnxB2Cc2peZ6wC7ndI7ZoX5ZmQ0jZ9I83PKTTmz79PcuuVaVtYon8qgayYCsOyP73HF1nsB2OQIXSndNAZra/2Kfc8ef9RMKJqb2+qbPXvUIjCHQyl1K1VcM0k4V0xKSnQZGHvdfwe9bruGZVvT8Nn3MaK/mqXaV6ZWNXxd3p8sVxXn/S60H0w6k2hO7oXMzuHCxbfzh8/H4hnQp/NCdQIhWrpdevf2rzjUJD6BUSqxpHUfCWeVR+rWa/b4fTlZBbHzBdYOOZbp/BmAofdfzjTHbDYxAJqa2E0mV/Msc56qYBXKh19zSD2JHD9WcMT0CSTt2Nr2nAFulL17VTRLRzl0SIWkgjXdMJCAFnt6eucrkov6Q1z+j3P5llG4qWNC0RZ8Oeon2FtehzhYy+b6Ao4o2IsQkVdgbk7u/mxZkcQvaxKXrporcblUvzGtzHCKvTMhtxlFsR02PsUMHuYm1g6owdmvD9XrPTiqKrmI1/mMibANMmx7OTV1IZ9Wq4pPae++wstLh3DEz64n7cv5Lc5XU6Os9P371YRoZ63tAwfUeaxSv6E1CWexp6d3PB+7ie3BBw4XCziIh1PPdZORp3wZe3Y04txTzmb6UZBt0cdwAGlp8ZZAk6gEuljCPUAcDuXq8/naD5mcy7n8jL/hyoxt7OyaN9dR+uJnAHjckgaSkLv3sAh/vqSxOevpn7qHnQ2Z1NbCwrVpXMcsJtXNp2pry9jD6mplbW/ZEn0ulwMHYpNXpitISIu9Mz+mq3QDN72hQrLOYD5HJa1j3OU3412zlxQOUFomaNy+k12cQEFh23Sh8cThaJsjQ7tdNJ3F6VQuhNbuvGCkpvo/YxpU6enKUg10a5yZ+hXnVr3FEm9sM1/VF/anvlCtF/G61bn3/bCPOrzc0Hcejpp9nDBjJD++2wBboWybpHSZX0Fs/PZHji3xz5U1N/sXEkVLaWl8U/OGIyrFLoS4ELgXGAIcL6VcHAuhwuHz+YdTkcZi2/ftYe+73/Ah07js6BXcdWMmh/KOpVFAU+80BrOWzTsK2LlRXfG8vtaKGczKUpbTli3qvcej49A1ncfsOz5f+0sokpOVYg80KtLTlVIPVOyrX1mKd+2yDvsHbTbVRiTK1mMo9h0bVcM5U0Zw6s9VThrP2u/gG9i+por1u3rjFA00SCfrV9Zz7EUdEilirKrUIXpXzCrgfOCLGMgSEX3qVXGASItnAAw6YwDLnlcLBc67fTA1R42mMVOFMjZk5jFKfMuSH7LYuEbFgeUeEbl/vTvwetXNZK4UDVvOTKNpB9MVE4kP3aydakZfeTzq1ToCqyG3iKoJHV895PO1nSgOVazblGHHVhVfmFboT2qTO0ytsir/bherGcrInDL6sYk1G2MYWpRARKXYpZRrpJTr2v9kjKis5OQnzgNUfpiIkJKph2ZzA48D0Kd/y9ki6fZyzuB1HGxy8sxXRwKQP9hait3M4mcOibVi10RD66ye4TAtertdPRAKCtT3YpXa2az2FIjPF/z85sNkc5n6BzJK/EELzgEl5FLO9g11rGYo/QsPMippGV9vzrNMHePuJLEmTxctYiir8VHJvz4NP3bbfdufmH3y8yy46BHeYwoA9w95MWhh4Pw7r6SYrayuLqZf0jbLTEy63S2TPJmdXSt2TTSYxVU6uu7B7fYbF7GK2vF61UMjsE+bRWBa4/EoM35thbLOC/r4/Uj1+X0ZwEZWb+/NflLJz2lmav5Cqpp6seWLbZEJY8WA9E7SrqdWCPERkBvk0F1SyrcibUgIMR2YDlBUVBSxgC1Ytw47zUzmfT77ZnLIj9nXfs/Uj2dQSzIYJfSeOvt9xtwyNWiOioNDRnBN7z9wx/47OeuY7UD3xqSHok+flsu6PR41gaULYWiiITVVWd4dXQAVWHDbTGcRbW13sy/n5qokWlKqfh6skI7bqxT7yn2F5NorcLuzDx+TSS76ecr5a7VaKZWRbcNVcjI8BQe2V9HePe3cvYO6qf9B7ZUz6P2LC6P7pwDf35/h0nUvUn1SNXhbegnGFY3j0dMejbqNcLRrsUspJ0kphwV5RazUjfM8K6UcKaUcmWUmf+4o69bRnJbO0MHNlNdnBM+i1tjIVw8vopZkphV8yGW8xKW8zLH/NZFmb+hE0+c+cAJrLrmPy584Puhxh6P78yu3tsyzs9WwVceva6LB6ezcAqi8vJZ9r73Mk+09OMycOKAeNoWFUFKiEpsFG5VmZyizbA9ZnJi7qc3xIZn+BS6+3CSSB6hkO1U/ho/U6bXwAxpOP4vRtZ8x7elRUacIEPWHeOuBtXxxcBRVP8Yma2tHSazYiltuoeqnF3DM3HJYC/NnlTLt7pbW/6bf/42bll7G0JRSZsw9BZtswlZXTbMnfC+sOW4CHDch6LHkZJVoyW6nRY3HriQtrW3Egt0ePkudRtOVtO6PTmf4vE0+n4p2aZ00yyRwNCpEyxBes7B24HL+7Ez/8ODnDx3Z5nyjB+4Bw+uSVpBCarrS0O3ds6vvepnzUMEV2yji1nPW8si7nSt0A+DcuY0HuI0J6xuY+fzlbQwxy2d3FEKcJ4QoA8YA7wgh3o+NWCEYMICGkyYxYGIxZzGP2fN6tczj8MXn/Nfb4+jj2cMTbxUpa8BupzklusnQ1FQVctjZ4hedIVTGvhilbddooqY9i9ycIwqGz9d+VsrW95sw/DYZ7CF3QNuFUAWT/Mo+b3gWjtRkUjjA3qowQ1wpeb36DABGs4DfJs/ky4rBfPTkmvDChaFqXTk7yWPyhINxG11HGxUzV0pZKKV0SSlzpJShHd8xwmaD6hHjmX7cd1Q2p/PP1/eZwvD3mVsppZh7HkyNuJJQewjhX+XpdhN08jXW2GztF6PWaOJNJK6WwAla058uBBQVtV+Cr7WrpznJzW4y2UT/oJ/ff8pUnmQGT3A9bq+NZm8KWexm74HQjglH5U42NxUzwrGcF2bVc+ZdwzmGpdz+lyEsuzcCb3NjI561y3Ds2Xl48rXse7WSq+So+N3EiRUVg1+x9r3xHI5lMa89e4DmJsmOp97gsfKLOKPfGoaNil3YSGam33IQAoYMidmpQ5KVpS1zjfUJptidTr+VHrgNKhgA1P0USf9uHSQgZDOZVOLOD1FtWwgmP3I6kx9VUXBN3hQy2cOP1aHdsEkV29nIAAaOTKXuuJMQk07hkXFvAHD1/HPaDZTx/fMVnJddxPWnbWD3m18BsGGlciP0P6GLMrlFQMIq9oODjuaXo1ax4WAfjh8lOPsvU8lOquKGWZ33jQWj9Tyv2x2b9LhpaWoyNi1NRRqkpanIgIIC9dJorE6w1c+DBkHfvmrb6fTHvOflKZ+606n6eSS0HrXW9R2KtNv54a7nQn6natyZ7B87xRDQSaaoZG9t6NjMg6UVVJJJYbHxz9hspN53M1eh2ljw5k6cFdsp+fUZJG1c3eb7ctcuTuJzvmIcLzwPvWfP4p3lhfRzlZGXHz/1mliTp7R0hYy6bQInnf8ZnzOBYazk3sfz8GXEzqkVqjSZxxN6Qqg9kpKU4k5P19EtmsSmtcXucPjDIdPS/Ip/6FB/Xy8sjDy6zOVS97uZrKsxM5elizq2jj/DsZ+VdaEnx3ZsVKuX8gf6rbWm3un85rc+Pr1vIy89eogxB5Zy1ILnOGvFp/z85DcRV07DUaysrwXrMignnwJ7OXN2jWPOwyrc8veXtY3a6U4S1mIHkEUlvPCnH1nz8Du88lYvSkbGNjNWqNJg4eLICwpCK+ykJGXR+HxaqWsSHzNazPSVB1rwgXVUA/t6sPqq4TD97J29X7K91ZTVZlC+4QD5T9+NrXp/i+NlP6inRt6glpNy9Wecx0Xe+SypGcwnn9soJ59nay5lzNt3Mus3/onVzTvUQ+PlGQsO7zs17VsmXxt8HqC7SGiLHWDfyedHdT6XK3hOZa83tEsk2KRPQYGy8LOyVAhYU5OKFS4tVZn0cnLaxgFrNImMzabizquqYNOmloo9VsELSUkq5LF/f9VGRxeH/qrhcR5jGm9Om8ur9UO4csdz+P7n5sPHt+1UQhcWt3L622zkjO5L8yd23lo3BDuNNBnq8rUdY7leqnt55z43XlFHymXnsa3uPjZs8+A8/6yo/udYkPCKPRpcLjVM/P77tpVQMjNDL5tOTlaFhp1OpbhtNn+pLFDDTSHU+4EDVRY4XbZO01Pp3Vv5zQNDj2N1n6akqFj43r075wI9onYFU3iXF+qvAGDFF5v4i3Gs18IPKCuVZNsr8Xrbxl72KVb/xDcHR3C0dz0zBz7Hs8uPZ468kMrdzWRm29hR3YsCdyXCVsiuX/2W2Gaj7zwJ7YqJlrQ0db7WSjcwxDEYQijrOzNTWSx5eS2tFafT/751yJdG09MQQingrkglbSp0IZQxFYpQRtiOq3/Lg9xCHqoG3vaDmUgJSetX8uvrJS/WXEBJckXQ7w49ys5ZzAPgouzP8b0wk7MvUA0tfFXl0N5el0Feyv6g348n/9aK3QzFat0hi4sjz6MRLO2oRvPvhssV2+LbJh6PPzomNTV4mKTTqeLig1H+q3s5sHg9b7/nYObwl9jbnErp1iZK5y7hA9Sym8yUg0G/29SnL/M4h+2DJjLlSRVpM3w4DGQ9j77ko/RfO9nelENeevDvx5OEc8XEMr47MN42kFCTphqNJjhJSV1Tp9Vma2mpJyX50xiYtVvz89U9G1jLtQ1Z2Rw12gsr4LP7FzJ/6XgAxvM590xcDhzT5isH+w5hxTvbaMgpPLzPVZjF25zFMFYx5w8rKWcsuZmdqIjdxSScxS5EdMrddL243f4Z90DFbrN1jeWh0fRkTHdMVxA4Sg/MD9/fCDwxLfr2HiyFI3NI4QBPLD2RLfRjtGsps6/5Eue1V4X8TqBSB6gZPgae+TOnJX/F7G3jkdjIHmg9SzDhLHbw1wDtDKmpyjceWNYq0uK+Go0mNO1le4wFppJ3u/3ph817NjBpWEaGCogILLnX0G8Io1nIF4znytS5XPzfQ9l54t0dlqF65ATGHvEa85ep98dOLQn52e5I+BWMhFTsnbXYbTZ/6ttAApW5VuwajXUxFbtpjAWuYg1cX5KerpR8oGJvSsvg6V63UVfdRMPHy6OS4/SJdSxY9ibHDz9EZv7PojpXV5CQir2zs+8+X/Dhotvt989pxa7RWBdTsQfTAYGK3eMJHmhRM+8TRGOUCdeBpnMv4Omah9l16U1Bi/fEm4RU7J1VvqEyJgqhJlKrq7tnOKnRaDqHOVoPNg9mKnabrWURj0CaesWm7mWzN4Xyq++Jybm6goSbPIWO1/y02dSsebjJHXOpsy47p9FYF1NZBzPAAn3toKz6rprQtToJqdg7mqvc61UrQMNZ4z6f6jS6ULRGY11MxR5s1O5w+K11k3/XugYJ6YpJSmpbNqs1Ho+aMW9s9GeHC4fdDkceqfOgazRWxrw/Q7ljWy+U6orVsIlAwv7bXm94xZ6d3fFERNq/rtFYm6Qkle4j1L3q8bRU5h0x1NLT1TxbtMWsTeKZSiQhXTEQPv1nJPUUNRpN4uFyhS9EYxbzMInUYnc6lc7o2zc2o3a3u+MpimNJwlrsqanq6dy6SrrbrXK96PS4Gk3Pw0wIForevVu6XiNV7AMG+C3szEzYtavzMkL8E/8lrMUO6gK0pqSkewpOazSa7qc9g83pbOl/j0Sx2+0tFXFhYfSTrvGOrktoFZiV1fIipqfrqBaNRuMnEsUeTAkXF3dsvUxrY1Jb7FEghFLuHo9S6CUl8ZZIo9FYCbu9/RF8MMXuckXuI3e5lCsnsM14x88nrI/dJCdHRcCA9qtrNJq2OJ3By1+ahLLMc3OhokIlHAwXXp2VpVw3OTlQWRm7CdhoiMpiF0LMFEKsFUKsEELMFULEZr1uh+XQSl2j0QTH4wl/PFToZOCCRbNiWnp628+ZFn9BgSpWb4V6DtG6Yj4EhkkphwPrgTuiF0mj0Whih5kBMpTxF279imnNJyerSdWiorb+czO8Uoj4T5qaRKXYpZQfSCnNzOYLgcJwn9doNJruJjlZ+byLioIr93CTpMnJSpG7XMq94nC0Vd5WLMwTSx/7L4DXQh0UQkwHpgMUhSpQqNFoNF3AwIHqb0WFf+2Lw6EUfTiLPSOj7SRqoGIXwpqKvV2LXQjxkRBiVZDXOQGfuQtoBF4JdR4p5bNSypFSypFZWVmxkV6j0WgiwGbzl8U0ycoKXQQ7kNZWfm6uP+rFqmlI2rXYpZSTwh0XQlwJnAn8VMqQpWQ1Go0m7ni9KnIFVM6ZzsSbC6FCq1etSmDFHg4hxGnArcBJUsra2Iik0Wg0XUNysnplZ0e3iMjlUi6Z9iJu4kW0PvYnARfwoVDjlYVSymuilkqj0Wi6gORk6N8/Nn7x9HRl9VuRqBS7lHJA+5/SaDQa6xCryc78/NicpytI6JQCGo1Go2mLVuwajUbTw9CKXaPRaHoYWrFrNBpND0Mrdo1Go+lhaMWu0Wg0PQyt2DUajaaHoRW7RqPR9DC0YtdoNJoehohH3i4hxG7gh05+PRPYE0NxYolVZdNydRyryqbl6jhWla0zchVLKdtNjxsXxR4NQojFUsqRqDHy6wAABORJREFU8ZYjGFaVTcvVcawqm5ar41hVtq6US7tiNBqNpoehFbtGo9H0MBJRsT8bbwHCYFXZtFwdx6qyabk6jlVl6zK5Es7HrtFoNJrwJKLFrtFoNJowJJRiF0KcJoRYJ4TYKIS4vZvb/l8hRIUQYlXAPp8Q4kMhxAbjb7qxXwghHjfkXCGEGNGFcvURQnwqhFgthPheCHGDhWRzCyG+FUIsN2S7z9jfVwixyJDhNSFEkrHfZbzfaBwv6SrZjPbsQohlQoj5VpFLCLFVCLFSCPGdEGKxsS/u19JoL00IMUcIsVYIsUYIMSbesgkhBhm/lfnaL4S4Md5yGW39p9HvVwkhXjXuh+7pY1LKhHgBdmAT0A9IApYDQ7ux/fHACGBVwL4/Abcb27cDDxjbU4D3AAGMBhZ1oVx5wAhjuxewHhhqEdkEkGJsO4FFRpuvAxcb+58BrjW2rwOeMbYvBl7r4mt6EzAbmG+8j7tcwFYgs9W+uF9Lo70XgauM7SQgzSqyGW3agZ1AcbzlAgqALYAnoG9d2V19rEt/6Bj/UGOA9wPe3wHc0c0ylNBSsa8D8oztPGCdsf1n4JJgn+sGGd8CTrGabIAXWAqMQi3KcLS+rsD7wBhj22F8TnSRPIXAx8DJwHzjRreCXFtpq9jjfi2BVENRCavJFtDGqcDXVpALpdi3AT6jz8wHJndXH0skV4z5Q5mUGfviSY6UstzY3gnkGNtxkdUYvh2DsowtIZvh7vgOqAA+RI269kkpG4O0f1g243gVkNFFoj0K3Ao0G+8zLCKXBD4QQiwRQkw39lnhWvYFdgP/Z7ivnhdCJFtENpOLgVeN7bjKJaXcDjwIlALlqD6zhG7qY4mk2C2NVI/auIUYCSFSgH8AN0op9wcei6dsUsomKeVPUBby8cDgeMgRiBDiTKBCSrkk3rIEYayUcgRwOjBDCDE+8GAcr6UD5YqcJaU8BqhBuTisIBuGr/ps4O+tj8VDLsOnfw7qgZgPJAOndVf7iaTYtwN9At4XGvviyS4hRB6A8bfC2N+tsgohnCil/oqU8g0ryWYipdwHfIoafqYJIRxB2j8sm3E8FajsAnFOBM4WQmwF/oZyxzxmAblMSw8pZQUwF/UwtMK1LAPKpJSLjPdzUIreCrKBehAulVLuMt7HW65JwBYp5W4pZQPwBqrfdUsfSyTF/i9goDGrnIQads2Ls0zzgGnG9jSUf9vcf4UxAz8aqAoYFsYUIYQAXgDWSCkftphsWUKINGPbg/L9r0Ep+AtCyGbKfAHwiWFtxRQp5R1SykIpZQmqH30ipbw03nIJIZKFEL3MbZTPeBUWuJZSyp3ANiHEIGPXT4HVVpDN4BL8bhiz/XjKVQqMFkJ4jXvU/L26p4915WRGF0xITEFFfWwC7urmtl9F+coaUNbLL1E+sI+BDcBHgM/4rACeMuRcCYzsQrnGooaZK4DvjNcUi8g2HFhmyLYKuMfY3w/4FtiIGjq7jP1u4/1G43i/briuE/BHxcRVLqP95cbre7OPW+FaGu39BFhsXM83gXQryIZyc1QCqQH7rCDXfcBao++/BLi6q4/placajUbTw0gkV4xGo9FoIkArdo1Go+lhaMWu0Wg0PQyt2DUajaaHoRW7RqPR9DC0YtdoNJoehlbsGo1G08PQil2j0Wh6GP8Pp+nWPdeRLzsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_reconstruction_forecasts(reconst_mean, reconst_std, forecasts_mean, forecasts_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Level Trend ISSM\n",
    "\n",
    "We can model a piecewise linear random process by using a two-dimensional latent state $l_{t}\\in \\mathbb{R}^2$, where one dimension represents the level (again with a slight abusing of notation, $l$) and the other represents the trend (slope) $b$. \n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "  l_t &= \\delta l_{t - 1} + \\gamma b_{t - 1} + \\alpha\\cdot\\epsilon_t\\\\\n",
    "  b_t &= \\gamma b_{t - 1} + \\beta\\cdot\\epsilon_t\n",
    "  \\end{split}\n",
    "$$\n",
    "\n",
    "\n",
    "In ISSM framework, such a (Damped) LevelTrend-ISSM is given by\n",
    "\n",
    "$$\n",
    "  a_{t} = \\left[\\begin{array}{c}\n",
    "    \\delta \\\\\n",
    "    \\gamma\n",
    "  \\end{array}\\right], \\quad F_{t} = \\left[\\begin{array}{cc}\n",
    "    \\delta & \\gamma \\\\\n",
    "    0 & \\gamma\n",
    "  \\end{array}\\right], \\quad g_{t} = \\left[\\begin{array}{c}\n",
    "    \\alpha \\\\\n",
    "    \\beta\n",
    "  \\end{array}\\right],\n",
    "$$\n",
    "\n",
    "where $\\alpha>0$, $\\beta>0$ and the damping factors $\\delta, \\gamma \\in (0, 1]$.\n",
    "Both the level and slope components evolve over time by adding innovations $\\alpha \\epsilon_t$ and $\\beta \\epsilon_t$ respectively, so that $\\beta>0$ is the innovation strength for the slope. The level at time $t$ is the sum of level at $t-1$ and slope at $t-1$ (linear prediction) modulo the damping factors for level $\\delta$ and growth $\\gamma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "latent_dim = 2\n",
    "T          = len(ts)\n",
    "\n",
    "# Set the coefficients of the ISSM\n",
    "damp_fact = 1.0\n",
    "damp_growth = 1.0\n",
    "\n",
    "# Set the parameters of the ISSM\n",
    "alpha      = 0.5 \n",
    "beta       = 0.1 \n",
    "g_t        = nd.array([alpha, beta])\n",
    "g          = nd.repeat(g_t, T).reshape((latent_dim, T))\n",
    "\n",
    "# F and a are constant over time\n",
    "F_t = nd.reshape(nd.array([damp_fact, damp_growth, 0, damp_growth]), (latent_dim, latent_dim))\n",
    "a_t = nd.array([damp_fact, damp_growth])\n",
    "F   = nd.repeat(F_t, T).reshape((latent_dim, latent_dim, T))\n",
    "a   = nd.repeat(a_t, T).reshape((latent_dim, T))\n",
    "\n",
    "m_prior    = nd.zeros((latent_dim, 1))\n",
    "S_prior    = nd.zeros((latent_dim, latent_dim))\n",
    "sigma      = 0.5 * nd.ones((T, 1))\n",
    "b          = nd.zeros((T, 1))\n",
    "z          = nd.array(ts).reshape((T, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu_seq, S_seq, _ = ISSM_filter(z, b, F, a, g, sigma, m_prior, S_prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let us use the same cofficients (constant over time) for the future as well\n",
    "forecasts_mean, forecasts_std = forecast(mu_seq[-1], \n",
    "                                          S_seq[-1], \n",
    "                                          F, a, g, sigma, horizon=13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the reconstruction as well as the forecasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VFX6xz93JpMy6aRASCAU6S00ERF7QVddRVd0bejatrC6u+q6/ljX1bXt2lbdVbHhFhHXwtrbAgqiCCgiPQQIJCGkJ5NJpt/fHyd37sxkWpIhjfN5Hh4mc+vM3Ps97/2e97xHUVUViUQikfQfDD19AhKJRCKJLVLYJRKJpJ8hhV0ikUj6GVLYJRKJpJ8hhV0ikUj6GVLYJRKJpJ8hhV0ikUj6GVLYJRKJpJ8hhV0ikUj6GXE9cdDs7Gx12LBhPXFoiUQi6bNs2rSpRlXVnEjr9YiwDxs2jI0bN/bEoSUSiaTPoihKaTTrSStGIpFI+hlS2CUSiaSfIYVdIpFI+hkx89gVRTECG4FyVVXPjdV+JRJJ13E6nZSVlWGz2Xr6VCRRkJiYSEFBASaTqVPbx7Lz9GZgB5AWw31KJJIYUFZWRmpqKsOGDUNRlJ4+HUkYVFWltraWsrIyhg8f3ql9xMSKURSlAPgB8Hws9ieRSGKLzWYjKytLinofQFEUsrKyuvR0FSuP/XHgdsATo/1JJJIYI0W979DV36rLwq4oyrlAlaqqmyKsd4OiKBsVRdlYXV3d1cNKJBKJJASxiNjnAOcrirIfeBU4VVGUfwWupKrqElVVZ6iqOiMnJ+LAKYkkIhZLT5+BpCvcfffdPPzwwyGXr1ixgu3bt3fjGfUfuizsqqr+TlXVAlVVhwGXAitVVb2iy2cmkYRBVcFq7emzkBxJpLB3HpnHLumTuFxC3CV9i/vuu4/Ro0dzwgknsGvXLgCee+45Zs6cyZQpU7joootoaWlh3bp1vP3229x2220UFRVRUlISdD1JcGJaK0ZV1dXA6ljuUyIJhssFHtlV3zluuQU2b47tPouK4PHHw66yadMmXn31VTZv3ozL5WLatGlMnz6d+fPnc/311wOwePFiXnjhBRYtWsT555/Pueeey8UXXwxARkZG0PUk7emRImASSVeREXvfY82aNVx44YWYzWYAzj//fAC2bt3K4sWLaWhooLm5mbPOOivo9tGuJ5HCLumjyIi9C0SIrLubhQsXsmLFCqZMmcLSpUtZvXp1l9aTSI9d0kdxu2XE3tc48cQTWbFiBa2trVgsFt555x0ALBYLeXl5OJ1O/v3vf3vXT01NxeKT+hRqPUl7pLBL+iRut4zY+xrTpk1jwYIFTJkyhbPPPpuZM2cCcO+99zJr1izmzJnD2LFjvetfeuml/OUvf2Hq1KmUlJSEXE/SHkXtgbBnxowZqpxoQ9IVKirAZoMRI3r6TPoGO3bsYNy4cT19GpIOEOw3UxRlk6qqMyJtKyN2SZ9Ei9hbW3v6TCSS3ocUdkmfRPPYm5rA6ezps5FIehdS2CV9Ei1id7nEP4lEoiOFXdIn0fLY3W4RscsMGYlERwq7pM+hqnoeu9stOlHr63v6rCSS3oMUdkmfw2YTwm63i/8tFmhu7umzkkh6D1LYJX0Oq1X32K1WaGkRIi+RRGLp0qVUVFTEbH+PP/64XzGyc845h4aGhpjtv7NIYZf0OSwW3VN3u8HhEFG8HLDUd1BVFU8P/GDhhN3tdnd4f4HC/v7775ORkdHp84sVUtglfY5guesOhxT23s7+/fsZM2YMV111FRMnTuSf//wns2fPZtq0afzoRz+iuc1P27BhA8cffzxTpkzh2GOPxWKxYLPZuOaaa5g0aRJTp05l1apVgBDq+fPnM2/ePEaNGsXtt98OCJFeuHAhEydOZNKkSTz22GO8/vrrbNy4kcsvv5yioiJaW1sZNmwYv/3tb5k2bRr/+c9/OPnkk9EGT9bU1DBs2DDv/m699VYmTpzI5MmTefLJJ3niiSeoqKjglFNO4ZRTTgFg2LBh1NTUAPDoo48yceJEJk6cyONt9Xn279/PuHHjuP7665kwYQJnnnkmrUdgMIYsAibpc7hcoDiFkqsJid733W6Ik1d0RHqoai8AxcXFvPzyyxxzzDHMnz+fTz/9lOTkZB566CEeffRR7rjjDhYsWMDy5cuZOXMmTU1NJCUl8de//hVFUfj+++/ZuXMnZ555Jrt37wZg8+bNfPvttyQkJDBmzBgWLVpEVVUV5eXlbN26FYCGhgYyMjJ46qmnePjhh5kxQx+8mZWVxTfffAPAM888E/S8lyxZwv79+9m8eTNxcXHU1dUxYMAAHn30UVatWkV2drbf+ps2beKll15i/fr1qKrKrFmzOOmkk8jMzKS4uJhly5bx3HPPcckll/DGG29wxRWxnZtIRuySPofLBXFXXMLrcx/3G5wU6Um6qSnyOpIjS2FhIccddxxfffUV27dvZ86cORQVFfHyyy9TWlrKrl27yMvL89aRSUtLIy4ujrVr13rFb+zYsRQWFnqF/bTTTiM9PZ3ExETGjx9PaWkpI0aMYO/evSxatIgPP/yQtLS0kOe0YMGCiOf96aefcuONNxLXFjkMGDAg7Ppr167lwgsvJDk5mZSUFObPn8+aNWsAGD58OEVFRQBMnz6d/fv3Rzx+R5HxjaRPoeWvn1HyLFUMZOJOO2MmJQCRRbu1FUwmSErqhhPtxfRk1d7k5GRAeOxnnHEGy5Yt81v+/fffd3ifCQkJ3tdGoxGXy0VmZibfffcdH330Ec888wyvvfYaL774YthzAoiLi/N6/zabrcPn0pnzPRJWjIzYJX0KbZRpFQMBUA+WeZdF8tgdDjlKtbdw3HHH8cUXX7Bnzx4ArFYru3fvZsyYMRw6dIgNGzYAolSvy+Vi7ty53lK9u3fv5sCBA4wZMybk/mtqavB4PFx00UX86U9/8lotgaWAAxk2bBibNm0C4PXXX/e+f8YZZ/Dss8/iaruA6urqwu5v7ty5rFixgpaWFqxWK2+99RZz586N+vvpKlLYJX0Kux0/BXe16F5MpIjd6ZTC3lvIyclh6dKlXHbZZUyePJnZs2ezc+dO4uPjWb58OYsWLWLKlCmcccYZ2Gw2fvazn+HxeJg0aRILFixg6dKlfpFvIOXl5Zx88skUFRVxxRVX8MADDwBiso6bbrrJ23kayK233srTTz/N1KlTvZ2gANdddx1Dhw5l8uTJTJkyhVdeeQWAG264gXnz5nk7TzWmTZvGwoULOfbYY5k1axbXXXcdU6dOjcVXFxWybK+kT7F/PzTsrWPqacLjXHLLdqZdMR6AwkII6MPyo7gYMjPDr9NfkWV7+x6ybK/kqMHhAIOl0fu3065H75EidrdbRuySowMp7JI+hdMJ9gZ9QIjDR9gjeezaYCaJpL/TZWFXFCVRUZSvFUX5TlGUbYqi/DEWJyaRBMPlguZq3RtV9pZ4FT1SxO7xSGGXHB3EImK3A6eqqjoFKALmKYpyXAz2K5H4UVvbVvSrTldn8+r3Gfz0YkBG7BKJRpeFXRVotfVMbf9kdWxJzNFKcjTX6Zkwz3E96948DETnscvZliRHAzHx2BVFMSqKshmoAj5RVXV9LPYrkfgKsfba0qAr+AaO5eLGFzBVHgw72YYm+toAJ4mkPxMTYVdV1a2qahFQAByrKMrEwHUURblBUZSNiqJsrK6ujsVhJUcBvtaJV9gb2yvz5HOHkvjJ2yH3E6yBkHQ/TzzxBOPGjePyyy/v6VNh8+bNvP/++z19GkeEmGbFqKraAKwC5gVZtkRV1Rmqqs7IycmJ5WEl/RRtpiQNLepuCjJwsII8kr5aHXJfUth7B3//+9/55JNPvKNIw+E6wrmpUtjDoChKjqIoGW2vk4AzgJ1d3a9Eok2m4fu3edsGtn3RgBH/mz6fCt7cMTbkvnzFXOay9ww33XQTe/fu5eyzz+aRRx7hggsuYPLkyRx33HFs2bIFgLvvvpsrr7ySOXPmcOWVV+J2u7ntttuYOXMmkydP5tlnn/Xu76GHHmLSpElMmTKFO+64A4DnnnuOmTNnMmXKFC666CJvrfT//Oc/TJw4kSlTpnDiiSficDi46667WL58OUVFRSxfvpzPPvuMoqIiioqKmDp1atjSA72dWBQBywNeVhTFiGgoXlNV9d0Y7FdylON2+3eIejywcdlu3iT4Y/yLu+fykxD76kgVyP7OLR/ewubK2NbtLRpUxOPzwlcXe+aZZ/jwww9ZtWoVf/zjH5k6dSorVqxg5cqVXHXVVWxuqyW8fft21q5dS1JSEkuWLCE9PZ0NGzZgt9uZM2cOZ555Jjt37uS///0v69evx2w2e2u3zJ8/n+uvvx6AxYsX88ILL7Bo0SLuuecePvroI/Lz82loaCA+Pp577rmHjRs38tRTTwFw3nnn8be//Y05c+bQ3NxMYmJi8A/SB+iysKuqugXoviIIkqOGYMJeXJcVeoMw+Y6+Xr2M2HuetWvX8sYbbwBw6qmnUltbS1NTEwDnn38+SW0lOD/++GO2bNniLcjV2NhIcXExn376Kddccw1msxnQy+hu3bqVxYsX09DQQHNzM2eddRYAc+bMYeHChVxyySXMnz8/6DnNmTOHX//611x++eXMnz+fgoKCI/cFHGFk2V5Jr8Xj0bVaVcVrE3ronU8Z5eg3X5wSWrFlxK4TKbLuaXzL6KqqypNPPukVaI2PPvoo6LYLFy5kxYoVTJkyhaVLl7J69WpAPC2sX7+e9957j+nTp3srOPpyxx138IMf/ID333+fOXPm8NFHHzF2bGh7rzcjSwpIei2+wq7939iqV/T7lqms51jv3y1OU8h9+UbpMmLveXzL8K5evZrs7Oygk2GcddZZPP300zjbWubdu3djtVo544wzeOmll7weumbFWCwW8vLycDqdfh20JSUlzJo1i3vuuYecnBwOHjzYruRuSUkJkyZN4re//S0zZ85k586+21UoI3ZJryWYsNc2xwOwmSnkUEMONaxjNsfzJc2u4J6ozSYm2Ujcux374OG43Uf5TBu9gLvvvptrr72WyZMnYzabefnll4Oud91117F//36mTZuGqqrk5OSwYsUK5s2bx+bNm5kxYwbx8fGcc8453H///dx7773MmjWLnJwcZs2a5RXu2267jeLiYlRV5bTTTmPKlCkMHTqUBx98kKKiIn73u9+xdu1aVq1ahcFgYMKECZx99tnd+ZXEFFm2V9Jrqa8X09kVFgpx3rYN/nj2V5RUp7INfajE9lc28+crvuNTw5mUOQf57UNVYetWcFbW8OCZKzl3xmGO/ecixo4VjYXhKHlmlWV7+x6ybK+kX+Lx6KNENfukttVMLlVYJ+gWTOvoKSQMGUiz2+y3vd0O27eLjtOWNRv5D5dw9cZFYrIORKMhkfRHpLBLei1ahynonZ819lQyBxjY+cJa3OmZVF55KwDJSR4sarJfuQCXS0T6AHXF9X7vW6167RmJpL8hPXZJr8U3Ync4QLHbqHJmMDetBEO8icptdVRWiuXmJA8ejLS0gJZU4dtJWlmnX+ouF1RUQHx8N32QXoKqqiiK0tOnIYmCrlrkMmKX9FoCI/aWF16hgUwKbbuIi/P3x9vSmWms95kP1UfYq+t1FW9sFDbM0ZQdk5iYSG1tbZcFQ3LkUVWV2traLg2QkhG7pNcS6LF/vlt0jJ54+VDi4sA3+EzWhP1wK4MLRMjum6/e3Kqv3FjnIisr7qgS9oKCAsrKypAF+PoGiYmJXRogJYVd0mvxjdhdLqhuTkLBQ+LF57aL2JNSxR+Nla1AsncbDWur0fu6sbIVRqUeVcJuMpkYPnx4T5+GpJuQVoyk1+Kbx+62u2jZvJtsaogzKe0idnOKEO6mKrv3Pb+I3abHMLXlokf1aBJ2ydGFFHZJr8W06SvSVvwDgJbX32cJN1JNLgBGY4AVky5GnTbV6kVhfIW72W4iG2FDVB0SC4720gKS/ou0YiS9loQL5mEjHe68in2bG/2WGY0BVkyaEHZLCGG3OBIooAwbidRUiNxJzeo5WgYpSY4e5CUt6bVM4TuGUQpATZUIr1dzEtA+Yk/KEDVkLPV6GO4XsTsSSElWyaecmnLdrok0AbZE0heRwi7ptRygEBCRdVWTSP1K/5WouK5F7GlpQuBNqWK5rVlXc1+rpcmRgDndxOC4KqpqlKDrSCT9BSnskl6PwwE1VjMmHLh+fBWgR+xmM8TFQVyqKOzV0qQLuzZa1dDcRKUnl5wBHgaZm6hs1qsISmGX9EeksEt6JU0+lkpDrZtmu4l0Q7PXftEidqMRTCaISxeJ7K3NLqxWYcNoNotSfpAqBpKda2BgWiuVjkxvfrwUdkl/RAq7pFdirdc7QesqHVgdcSQbbd73tIg9Lk4Iuyk5AQNubDY4eBCam/V9NZSKal8D8kxkZIALE1arWCZTHiX9ESnskl6Ju1UX9sYqOy1OE8lxQtgNBiHoRqP4l5QEmYMSSMCO3a5gtUJpqb6v2vJWALIHJ5CaKRLBLA0iVG9s9J82TyLpD0hhl/RKPHZ9LruGGhdWVwLJcUKBk5KEsCckiGg9LQ0S0+JJxIajbTPfSLzmkHgza0gyqRliIJO1WpR2rK2F4uJu+EASSTci89glvRJXix5GV1UIYTcni/dSUoSwg+g8dbtBwUQCdhzO9rFKTVt5lMxhaXgyxIbNlc3gNoPRKCN2Sb+jyxG7oihDFEVZpSjKdkVRtimKcnMsTkxydGOz6Grb3OjC6k4i2SQi79RUfT2DQUTtiUlKm7DrqYzjLp9G/s8v4OU1wzHgJmNgIilZIt/d/NRDTJ8lRN7jkV67pH8Ri4jdBfxGVdVvFEVJBTYpivKJqqrbY7BvyVGK3Scf3dnqwupJxJwgOkETEtqvHxcHiYoDh0sXdvOub3mYcylmNCD8+PRcke++83Am6zmWuLahp06n/hQgkfR1uhyxq6p6SFXVb9peW4AdQH5X9ys5umn1FfamViyeZBJTRNmAUAIcrziwa1aMy8Vi7uUP3APAK7m3oCiQlZ9IIq3cwx84jvUoDWJmJacz+D4lkr5ITDtPFUUZBkwF1sdyv5KjD0eLLuyug5XUkkX6oERvimMwEgxOHG7ROWpsbuA+FgNwj3IX4/9xB5mZYMzOZAy7vNscLm4Qx5BWjKQfETNhVxQlBXgDuEVV1XbTBCuKcoOiKBsVRdkoi/1LImGz6iOHKg/Y8WAkqzAlrF2SaHDicIlL2lhfSxxOLuMVLvzHfFzZg4iPB7VwGBOSSrzbHNghUiHlQCVJfyImwq4oigkh6v9WVfXNYOuoqrpEVdUZqqrOyMnJicVhJf0UVfWP2A/UiFGlkYQ93uDE7hIrOA7X48LE0B9Owz6uSCyPB6PJwKRj9IFOjW3122XELulPxCIrRgFeAHaoqvpo109JcrTj8YDDppddPGjJBCCzMI2kpNDbmY0OWlxiblNLpRh6mpKt97SaTMLGOenhc3lqziuAPkdqJGFXVSn+kr5DLCL2OcCVwKmKomxu+3dODPYrOUrxeMBp072Rww4h7OkD4sjICL1diqmVZpfIerFUCYslLVdvCbQRq/FZ6Rx/91nE4aShUWTRRMpldzjkCFVJ36HLCV6qqq4FlIgrSiRR4HaLf74R+2FVWHepqSLqDkWayUZzc5uw1wiLJWVQsne5VjAMwJOaRha1NFhEbNPaGv687HZ9Yu1AGhoI2+BIJN2NLCkg6TU0Nooh/nY7OO26sDsQdkpycgRhj7dhcQs/vqlW+Capg1K8y32FnTgTWUodDVbxRjjhBpEO6WvF2Gx69cjDh2Xnq6R3IYVd0iuwWqGkRETONhs47P4qm2Zs9nrkoUiJd+BQ43E4oKlRqG5auv4waTSKDlSNdGMzza1G79/hPHS3G5p8cr2am8U519eL13Z76G0lku5GCrukV1BfLyLmhgYRtbcT9rgWcnKEOIcira3UQEsLNFqEoKen68sDhT01zobFrr8RbpCSyyWeKEA/T6sVKivFe1LYJb0JKeySXoFvffTWVrDXiKyWVESYnG5qDVpKwJfkbOGv511/DtayRlIMLX7WjVZXRiMl3k6zU99ppIjd7dYFvrERqqtFIwJy5KqkdyGFXdIr8O28jC/bS8Pn32HCwWh2A5CWYAvrrwNk5Ynou3FfPfX2JDJM+mwbmoWjzbwEkJLgxOLUs2bCibPmodvt+nq+WTKhPHYtypdIuhMp7JIex273F8b40mIe5HdkUk+6wQJAWqIjorCPniCE/Tze4VumMjDJ4l3mu632OjXR6e1shfAdoFo073QGbwACo/2aGvHeoUPhz1kiORJIYZf0KG63/zR2AHuLhcKmYiHVIDya9CRHWH8dYMjVpwJQSzbfM5lBqfqOgwl7SpIbq2rG3Sz8FI+eiNMOW9tgVZcreD57oLA3N4PFIqwamf8u6W6ksEt6lMpKIYC+NFeKN/6WeodX2FOT3F4LJRRxaWa/v7Vp8MBf2DVbJi1FKLn95LNI2rU5pLCrqi7OTmfwjtLAaL+1FcrLxbayY1XS3Uhhl8SMcBFvMFRVdEDW1vq/rw0usixZRrLa1omaHFnYFQWWLYOFC8Gc5GHC/DHeZb7ZMFrkf0ry1wCcyBo2/GN7SCumulrPcXc69ejdF197RlXbOoDbBF1G7JLuRgq7JGZUVXVsfas1uK/d0CxC6gG5Juxu8To1zRDRijEYYNQo+MUv4LPPDZx9bviIfWzcHp7jOgA++L4gZMPk+7k07zwQX2FvbfUf7CSFXdLdSGGXxAS3G+rqOrZNoLeu0dg2GnTAAHAj1HzaBDtKhMIVvssD1/UVdq2BOHjrE5z1m0nMStpCpTU1ZMTuK+ShRqc6nfqyQOvFV9ilLSPpDqSwS2KC3S4i1Y7kc2u564E02uIx4cBshofif8+rLGD85MhljcJZNcGE3TmwgKrLbmaQuZFDrRkhI/ZoLSbNogkUb60PobW1442fRNIZpLBLOoWveNXX66JWUxPd9g5H6Bxvqz2ONIMVoxEGGmtYwGuoWdkR9xkqos/IwK/cb2Da5MAUK4ec2UEF3OMJX0PGl507xfcSaL1ojZ7v9ySRHEmksEs6RWWlHnFbLFBRIV6HisKDbR9KMK0OE8lGG4oCarao7KhmZEbcZ7CIPSEBhg3zX6aNYNUEPivNiUVN9Y4i9aUjxb08HtFYBRPvsjKR0y6tGEl3IIX9KCLayDMaWlv9h9NrgqWNIPV4wlsYTe0mT9RpccWTbLRhMEDT8g8p/9l9GHI7HrEnJsK4ce3ryyQminXT0sTfWsnd2qr2vaIdrdpYVxdcvLXP2xVhr6kJ/71JJBpS2I8CtNl/Ghpit0+bTRd2385Fh0OUsa2ogAMHgm/r8YQXOKsznuS4ts7S0aOpvPZOTPGRS/4HRuxxccGLhhkMQtQ1eyatLd+94UB71fQKe5StotUaPgvG5ep8iV+bLba/oaT/IoX9KMBqFUIaK3/X4RDipIlzYPpfVVVbhcUQHnqkqNXqSiTZ5MBg0KNw3zz0cPhG7eFK/BYU+AxUyhaeTGNZ+zQd7bN9u+ABPj7/yehOIgKRJvUIhc0mrRxJdEhh7+eoqvDAnU4R7cXCjtHsAJvNf1SmhsMhUhldruBCFKkwltWdiDneiaLoUXikOjEavlF7uLz3xEQfYR8oQvemwy3tUjDdblDsNq7feyd3ViwKuq/k77/qUBgebT9EIDab7HyVRIcU9n6O5n+7XCKKDpU73hG0jlKnUwybD+alaw1IMBELlTnjdqkc3t1Is8eMOd6FwaALdbjo2xdfMY+0jbY8e6goRVBR5qG62j9l0+2GuAb9hLXPFVdzCKOlAfXrDfz2msM0PvJcdCdI53xyrQH1zZeXSEIhhb2fowm7JlaBw/c7iiYuGpFGm2oNSUODsCAsltB2wh3nbeUHP07noGcwSclGFEW3ViKNOtXwXS/SNpqwmwtzyaeMkl0u6uv9n0BcLrDV6K2h1q/w8bxH2X7RYjaubORtfsjP3zk7uhP02cfBg9FH4Npo1mBPSBJJIF2ezFrSu9Fqm2gedWdsAKdTt0ICRTlS9Gi1ioh+3z6RZhjKX17+qy9ZVT1bHI948vPVTkXsvlZMtBG7MzefUSkbOFCW5C3aZbe3jXx1Q9Nh/aSryx0MPPA+t/MXqIMrXv8nADZ3lF4Regdqc7M4h7y8yNtUV+uvbTYiTjoiObqREXs/R5uEWXv8DzdLUCh8o+yOdt45HELMPZ7Qot767U4eXjOLSWzxvjdsVHyXI/ZIwu7bCGSn2Kh2iHn06utF3jmI789Sq4fVH79aR+sdd3v//hdXAlDhyO5QtovWmR1NxO52+9tXne18lRw9xETYFUV5UVGUKkVRtsZif5LYoFVPBF3QXa6Oe7Q2my4mHRV2lyuyr//dWyWoGFh8h5vj+JKBVDL2rMIj7rH7ds5mpzmodg8AhG3kdIqc9Lo6aK7VvacDb3/Lqaz028+fZq7AQTw11QFfrKpi3rEp6BdeWysau2hslcA+CdmBKolErCL2pcC8GO1LEiOs1uDRXUfn5/RNsws2OjMSkXz9nbsU4rFTeM543h96E6t/upzUIRl+09gdCY/dd50BaW5aSKbFoofdpaXif0u9aBXHsoP1zKKeAYxlBwDXJfyT7LYp+ZrL6r3bmqrK4Y03uO/KHdhXvN/uuJpYR/NbBNarl/OrSiIRE2FXVfVzQJY36mWEigY7ase0tOi+c2c8+kjWwfZDmYw1HyDOnMDeN7/D+pObAdEvoEXtkSo7avhG6dFE+V5hF8E6jWW6imrZPpZG8WIGG6lHrHjTdS6qyeau6ypIzRaGt+WQ+HLi9u7m9nO+58QHz+bfXMH9/ypsd1xt39FE7IHrSGGXREJ67P2YUKLR0QqMdruIGhsaOj9qMhiKpZHqPzzJhpYJFBW0z4FMSBCiHq0NA+GLfQVDE/bkTLFya5V/LqKxqZ6mfcLPGjop3ft+xoR8Dq4uoWrh7aTliIjdUiM8kj0f7+Uj5tFCMgCNYdIbVTXy7xGNsHs8IkNJc31kSuTRTbcJu6IoNyiKslFRlI3Vvl38kiNGLCL2+jZ3wWYTpQJiybqSL2axAAAgAElEQVQ73+Ps9xZhIY2pk/2T4bOydEGPdnASCGEfMEBsE02Ur1k95lRxsNZGf9UcfePJlJU4yDceYvBAvVXLGpGBJyUdFIXUHNGaWGrEF773O38lL7dmhD0Huz20xeXxtG9Mg/WTtLaK30ezbRobO2ebSY4sbk8MI6MwdJuwq6q6RFXVGaqqzsjJyemuwx61BGZS+BLNxBEgonXfNjhcZKm0WKn9eGOHzvHbg1ne1xNP9i/y5ZvOF205ARDCPngwZEYuBum3b3OaCN1bmvw/ZFLxFr5hGqPj95OZo7cwuYP0Wyd1UFtkXi9u2r27/Pdx0DEw7JOO3S6qXQbrFI1m8g8Q2zocupi3tsZmMJoktrg8nUhL6wTSiumnWK2hRdtXFLRRpMGor49+komPLn2Rs+6cwb7/7Y1uA1Xl+6qBANzJfWTOGOG32LfjsyMRO4hGoaAgunWzskRBMHOaOEhLk7+SLjVey1Ymcd7IbQzI01sYv+ybrHTSaKS2zoBia2VPUy5TM/ex/oz/48ERS3BhovpQ6BtaS3sM1hcR6ukqmLCDfyd3R9MiZeXII49b7UMRu6Ioy4AvgTGKopQpivKTWOxX0nnC3dTaze/xBI/qWlrEssBsjFDEl+/jjQoxuOi7ddH1rroPlLHZMZ5fjPmYy968BOL81dvXV+/MYJxoO1tTU6GwEJLS24Td4n/jfZJ6MfmUMffpKxg2Kp7beYi/Jd/ut44rNYPBVFBdH0dCxT52MpZhQ1WMD9xH/qx8AA7vqCPx+SfZ8fMn24XhdXX+k1/77TtKYfedONtqFf86mhZZV9e5cQ6S6OlTEbuqqpepqpqnqqpJVdUCVVVfiMV+JZ0nnL+qiYD2+B5IQ4MYoBONR5u27kMm/XAELYh6K9/sSIzq/PasOoiTeMafOQT70FHtlvsKe7S2SmcxmcCcIVqPFqv/I0qjI4kBKQ6UpEQ8qWk8xB0sONm/s0FNSCLPcJiakgZS/vYgZQyh8BjxAQaOSAGgbEcT8585kyvXL6J2S5nf9tpv0FFhDyyXrP1fV+dfIz9atMFkkiNHv/PYJd2H3R7eX9UiOZsteFEpu93fW2+HqrL9vb3UVNgxle3lfn7HDsYD8NW+QVGl8H29shkDbsaeHtwz6YoV01EUBVKyRYNkDfjeGl1m0uLFF9Y6Zip771tG6e+eabeDIRxkPcfx489uBGDIGOG7D5owgKGU8sU/StjNGAD2bw3+4wT73kIJe1WVf21238ZaGzfgdEb4HYMc32KRna5Hkj5lxUh6FyUl4fOjPR7d1w1WVCpS1Nb09mqu+sMIbr6ili/LhvB/3A/ABYb/UudMZcPbYYx7oGXtJl7aPoszMzeQlp8adJ1oByTFClO6mRQsWANErcllJi2xrTNUUag/61LUxKR22/8q9XkA1jEHgIKJYnomZ+EoTuN/fMxZ3nVL9wRX60gRu3nbBhSb+HGam/XG2+PxX8/X6YlUpE1DGwXb2Cg7XY8kfcqKkfQenM7gwjz4kd+Q8Nar3r9bW/UsF01QKipgx47Iwr7+c7HBrqbBlFXol9CiW4woeDiwsoQh915PXGVZ8O2X7qCJdK69b2TIYwSbv/RIYkgxk4qFllafA6sqjZ5UUpMi34ym555mxcmPiX3hZsgw8ZihJiRiGiyywObk7CKLGvYdDN5qBWuMtd/Ic7iaS6+O591rXvcu0waLhWvEbbbofPO6OtHIt7ToKa6S2KKqqrRiJJ0jcBILpboK96vLuWPZZCbddylul/BdfMvvaiUDDh+O7jG8uFT3Rr4sEZkt65hN5oWnMIK9vPt1Dqf+92Z+dF4rdRXtd7jhwEDy4qoZPjN02mt3C7sxJYk0mnDsPUjq1/8DQHHYaSCD1OTIqUG2EeMZcudVDKWU4YZSvxTN6848wA08yxN/rGe8spOSQ+ag+1DV9uUXNFE+tOEgW5jCH4qvZE+xOB+tjyTSAKdoInDf7KiWFjnA6UjQXTYMSGHvdwQK+78ve5dZDy/gH1wNQM0+keridOqCYLWKPOpoUxt3Vw8gG2HeflExnDicJH3wFp6kZG7jL+xiLNuYSLE6iv89uqXd9sVNgxibcShs5kp3WzFxiXFtEbvC6J+dDoDSbMFCKqkp0amcOiCL39/h5FeL/G9gw89+yi0vT0M99jhGpVSwpz70xNz79/tnI2m/Udn3uqF+YOUe7+vm5sjCHim7qbLSfx/hKnFKOk932TAghb1f4HsTBqa4LW34od/f+74Tyq/NqgTiMTzUYKZAVKeTrdZhnJ22jjicVLmzyI5rwJ0zCICTnr7Uu+5wZT+bA7JkPBWV7HKNYOTg0Ll4vuV6u4u4OEijiQ84h5dYCEBLdTMqBpJToz+ZqRcfw3FXBmT5GAy0TJgJwMhcC3WutLB2h2+nqPbbHijR/ZaDu/XvrqUlsrCHmwDb4wn+20frzUuip7tsGJDC3qdpbRXiXFHhn8fsix3/JPBlz4rwraGhczPxrP/pyzSQycSxTkZSAsAQs+4fWGeewm95kPP5Lydk72BD9TC/J4FtT62ihWSmnzMw5DG624YBIexOhMV0LS8B0FIjVDU13YDB0LERsKEonCJyNyve+JJRN51KXM2hdutoT11NTXpHaOmhBPINhxjJHvaX6bmgzc2R89W1yceDodUCCqSurvNzs0qCIyN2SVQcPAiHDgmRbmkRN3ignWJEv6MnsYUv6sdTXtYxAzVxz1aMjaJ454ffDyabak6/fRo/4+8AtDj8q3RdvuISHvxHPtPGWKnzZFD1/Nssm/Ewq8/9C+99HEdWXANTL2hf8dB7zt1sw4AYBLWGE71/j7zlXForxGdOTjeRkBCbtMsRP5qGERdrlmznnI338MGvPm63jt0ufsfKSv29PQ1ZjEyrZrRpP3ur9Uwiq1WIcCRaWoJbMqEad1WVwh5rpLBLosJm0x+Zm5qCZDM4HLSip+bdYvgrAGuW7iFqVJVNl/4Zx9XXY2i1stE9lVPHVqAOG8EFl5g4nU/4v9n+E084CkbQMn4Gk04Qxa9+uqSIR7iVWytv43V+xHnDvw9ZsTEhoWeEPSkJ/sxt3r/vXXsKllWbAFH5MT4+tLAPDP3w0Y74UcMYbyrmec9P+IIT+N2Oq2lpan/D+5ZIbiytZ6N9MuMHNzAitZq9zTl+DXg0fSN798Lu3e2983A2jiwPHFuksEsi4nL533h1de2rL1rKGnCi+weFSxaTy2EObY1+Rmvl8CGu5h/MKnsD975SKsnz1mGpvv0RHnl3DGPuXxh028x5s0ingYMM9b6XTgMLfjko6Pp5eWKUaU/M52kywU+nrud5RDWMR/kNL22YAEBydiJxcf6jYZOT9QYoO7tj0Xxhsv/3v3VFcbt1LBZdsJ+/5gsAZo9vYHhuM62exA574Fp/iu9TAEhh706ksIdBpmEJAlPYgpV3bbnjXr+/MyYNJS+umqqm6Ib9A9RtPuh9Xb1NqMmgkcne95yDhqKagpvPakoqWQgR++kpO9jz4Ov876tkcma3LyGQkSGqMiYlQWL0pxdT9i/5mMzf3uD9+3+I7JiU0fnExfmXNkhOhpQUGDRInG9HGqMhmcITGZe4D4DDxe09Ek24kzat4YumScxgA0VXTGLYUHEDlBZ3ooOE9h2p4coOdLQkgSQ8UtjDIIVdEHGouMvJg3t/xEAq2ccwasgCo5FBCfVUtaREfZzybbq/s+Vb0XIMnjgg6u0PkQfA2MnxNJx+cbtiXxppYqAmmZkdszZiSVxKIvkXHMupw/dydvbXAEzIOczAgniMRiHkICL3rCwh7tltmYsd6VgdO0j0jmaY7SRgo7KqfdaNJqrrn/+eAxTyk9FrcBYMJ3vmcABqN+3r1GcMDADCibeWz97UJP32WNCdwt6BuWl6Bx5Pz2RN9CT19f6RrNsdOTfZU1XLGk7k12Pfw3nT37HUCZ9moNnC5vpxUR+7ZLdu4L7/cRxGXAwrir4ql5kWWjEz9cTgpQM0ktseAhSlYzMmxZK4OIgzKfz5PyNorMil4vpKfnNXsneZlhmTnAxms/g9OjrZNsCUH0+AL+GyHzRR88+DVNaGbhUO1ovWZMIj1wJgniF+u6YDAQMWXC5ADdlw+uJ0ChvJ5Qov7B6PEHStnK/2Gx06JBq51PA/adS43T3Tr9LdyIg9DEdjxN7YCHt8+jubmyN/D3Wl4m7MnZBL0wnnUHv+NeLvlFaqXAOiK8+qqny3Rx8l+TkncYLxSxKTos/rXsmp/IvLMQ8JPShHUXrOfvHF1ydPH5zCC+8NYvwsoV6a8Ph2ovoGGIHCrq2fm9v+OCmzJ7FpnZ3jbz6W/PjDVNQlMOQvv8TQ3Nhu3bp6AwmKg6RBoiPamJtNGo3UBXSU77rmAX46Zwute0QZhyG3X0rmOy8H/Zy+4xciXUcNDaLD1XdEstUau3oy9fXh5wToT0hhD0O0oyP7E9pk0tojdDSPxbWlYqXsPH/Fyc10omIQGTQuJ8am0CNlrNv28X7jHG7Me9v73oIbw0/zFoj54Xs5/bLcsI9ZSUm94yks3Dn4TtMXrKM0UNhzckQjkJMTfLCVGi9M+YKkeiqaUtizfCPxjz/Ubr2aVjPZCU3efagJieQq1dQ2+hzQ5eLyHb9ng3s6cy8tQLE2c9zK+1n8x+DRu28piUhoxca0dTWhj9XIVIulfSPR0NA/O26lsIfhaIrY3W7xz7fMLkQn7Fv/tgqA3KH+oXBOtvgCqytdDH3o5xSdOgBcQe4iVeXvC7/GhYnTryvkHtO9pNPAhFM6ZoA3nvxDyn7zWMjlmrXRGwgn7NoTRVxccD890EowmyE9XWwXzn/Py7JTTgFzWMevV57Xbnm1PZ2cJP8fPMfUQG2z3lsbv2+X3/KytfvYxwiW8eOgVktg8TcA3G4yVr7ZLnKyWsVbTqc+YbbDEbuOVW2mJ9/DamMy+htS2MNwNAl7aSkUF+ufWbuZIkZLLhdvtp7DbNYx+PjhfouycsRPXnfAyptvKaTTQNXq7X7rJJTuxrbma17nYiaxhWHzxnHB61fw+bJykocH8Ra6QH6+yIjpDYTyeY1GXZxNpuC2UWDEnpQkPhuEz5gZnK+H8x83zfYWaQPA7abKPYCsFH+Vy06wUNOqd4A3bCkFYDpiztnV933hXVaypb1nog1K8r2O0t5/lS9uX0Hzs/8Mea6treEnBekMWunowNry/U3Yu1PUoQ8Ke09ZMdF0WMYai8U/OtfqgkR6TPXsLmY3o5l5cnI7FcrKFypTV97CMi6jiXTWr9IN1PiyvdRcdAP/+fWXuInjzr8VoCTE48gfjnPUhLDH7ejITINBZMFoGTE9TaiIXRNoCJ3a6PvZDQaxnq8vH4rCCf4ZSlrfCICxuZEqcslK8xeFnOQWalrNDHj/XwCUbxWquOBUUfTlsZabvOvu/7J9yQKbTXjbvtfR+p3pXM0/uPP1GSHP1bfUs/Y0GY5IJSucTn0fhw+La10bIdvfJvtwurvXW+pzwt4TEbvFIgZ2dGedaoulfR3tqqrozmHf19WoGBg9pf2EEBn5ySh4qN9VRR0ibbHisB6qNmwr5xRW81duIVWxMGJG9KmNg4KPOwpJamr3F/sKR6iIPT1df52ZGXw9X2FPCvjawwl7/lxRkz4ZEVlXb9fzWA31tVSRy4AB/hf9gFQnNeRQetfzoKoc2CsulMkX+k8IbsJBRXF7366pSYxE9aXisJCCLxvHMeaCsUHPNbACaKQAo64ufOTt6623tMCBAyJJwOnsf5N9yIg9Aj0RsTc2iou6O0uZBstTV1VRHyYcarOVfz4jwp2Rx7W3TTz5Q8ilivqdhylGDBQqr9Wj+pLt+p140ejvo05DM5k6bqn0Fm9dI1TEHs13EBenN1KBEb32t6LoufDefY8eyTfXPsGKX4ga8IeKdUWzVjbiwsSALP8TG5AhboJTWI1aX8+BShNmQys5E/WW9dFT32UEe9nwpRuP0z+0DhYcVVaLY6gYGFH2edDPGGi/RBJ2my10f1Bzc/vBUr6VKm22/iXuUtgj0N0Ru8ulR8nd5fupaudtn/X3fMi7rnkAZB/TXmnthWMoNJaz5vBobG11ZMqbREpf6lcfc6hENAovcg2//GNW1MdNSBCRaWC0Go5Y5UHHiq4IO+gNW6Cwp6UJUdc6VP1QFDw/+yXpZ4iyvpUHdP+ioVz8Fpm5/h7X7IF6uF3yZTWHmpIZbG5ETU3jcW7mFFZy/D3zuJFn2cR0vnpmc8Rzr6zVT7qaXFzOyDdaNFaLbzCkibaqiukbIxUv60z10d6K09MHrRhFUeYpirJLUZQ9iqLcEYt9hqK7I/bKSv0Cc7nE4IwjTV1ddNOZBeOL7UJd3sm7IbjNoShMGlxDKcMAGBdXzEHrAAzrv+SGX8Rz51c/JFWxMOXzp/AcMybkcXxtlIEDhWiBGJEZDSkp7aPXniaYgHdk4MyIEeIpJFDY4+LE95WbG7oj1ZA3iMGUU17pY4tVihA5Y7B/aznyGIWlbROnbPnCQoUzh0FpIjS+eu4+ls9/jbjEOOYvKiCXw6xdHV5U1MpKllWe4rWDAL5b2xRmC0Gwa9Ru1/1xp9PfKy9rmykx2un6+lNJgz4XsSuKYgT+BpwNjAcuUxRlfFf3G4rujtgDpyqrqIhNjq02aXBLiz5jTUODmPQgkt0SioSS7XxWOZb5maso+NeDIdcbMUaPAI/N3UedK523bvmMzziZfMpYdMJmVHN4nyQ7W4+4ExNFnRfQBT4SOaFnxesxuirsIKLzYOJdUAADBoQZiGUwUGiq4Js9aSSvfg+A+mohBhn5KX7nUnXpL5n5r1sYpezhm29VyihgUKaIPkoee5sDdz4DQM3Vv2FW5m6+rigIe84f3SFSYxfxJOchxiysW3Yg4mcNvA9UVcyZW1qqL/ctS1BfL671wHsqFP1pFqe+2Hl6LLBHVdW9qqo6gFeBH0bYptN0Z8TuO8uQL13psbfZRMdVSQns2iVuhO++g+3bxXulpZGzDUJR8urXVJDPjAUjcae37/TUUvJy8nRhnz5GRGn3Ou9gTuY23vq6gPmPzY14rKQk3XYxmXTR0QYbxceHT/PrbTYMxEbYs7ODi7f2XYX7Tk7O2MweRrH4VqFodbXiYk8vEF9WVlbbE1GcidaxUzkudy/rqkdTTj652cEvmumjmtjvLKD6YHCVNG7bwr+2FgFw+e+G8jY/5ATWsKMkcvGbwHvD5RLXrpbZov29a5eetltd3b4KaSi6OwvtSNLnInYgH/CNMcva3os5jz8OV1xxJPYcnFAeYDivff/+0E8Vra3iIq+v91+nK42VYmlk04m34ProUzbsFjbM7AVD/dbRhFbzd7OG6I/2RXP1XMOFiwuiGgGqKEKgNLHyzeGOixPRaV4eDB0afHujMTaTVsSaWAh7fHz4ujEGQ+jPfu6jp1JoKucNLqZpRxnVh1zE4SR9gDiJxER/+2rqDCNNpKNiYGBe8B9u6lyxwYZlQWrwqyq7r/4TOxnHw6OXYL/ocra8W0rR4Co2Nw6ntib847GvZw7+Qu977/h2oJaXh92lH1rD4HuMvkqf9NijQVGUGxRF2agoysbqiKUJg1NSAuvXx+6ctCyToFkClaFrWASL2PfvF5F4ba0Q7vp62LlT9xMPHhTReWe980ASSnez//an+ObWV7ix5XFuvDuPg7VJ5Bpr2kXDqalCcLRIMnu4voL5jOP5P/7Ec8e/yNSTAnv2Qhw7Qc/VhvZClpkpLJmUlODpjB3pYO1OYiHs0RDKrkoYN5KH7xEX1/d/XcmW/WmMN+/3nkNiov/TwOQ5usrnDAn+pebPP46Rhn18tbp9NGKwNnEPd5FHBXNfWAiIMswXn9WMkzheu31j2M9htYrrXhNr32s7cFJ1jY4KtNa/VVLSse16G33RiikHhvj8XdD2nh+qqi5RVXWGqqozcjppsMbHx7aGRH29yA0P5uWFujBBPCJaLLrwO526mIPoJCotFRd+cbGYuaaqKrZRx7dPreXilb/g+k0/BeA75wRKG9IpNPs3mgaDEFiTSXi8BgOkjczleX7CkrGPQpKZi1f9gqmPXR31sQe0uTzaIJzAPO24OCFeBkNwEe8tI02DkZnp3xgdiUqT4dI8h506ksHKIR7feAKrOYWZ0/WLJiHB38oZMPMY7+uBM0P46AkJjM8oZ09VGo133I/nm2+9i1wHD7GNiVw++XuMSfqPmPKzq7gk4xNe3DIdS2no1BW3WwQytbXifvEV9ljdp3a77tGH89x7e0TfFyP2DcAoRVGGK4oSD1wKvB1hm05hMsVW2DUPL7D3/fDh8Dm0TqcQ7EOHhIAHG7ihPUI6HEemE2jDwbx2731mO44xhf4fxmQSYqDVOElLA3daJse/fBPTX/gZAO7UjA6FplrFQqMRhg0LP8goMDpVlOgzZ3qCwNmQjkTEro20DdZoKEYDE1JL2ccI4hUHV9wzGtDtNJNJT8v0ZGbxi6I1nFOwhcIRob2t0Vl17GYMp316J3/9gy7UNXtEIvnA6UP8N1AUTjrbjIqB2pKAZPMguFxC3GP1NOpLY6Mu7oGzP/kSbYdsT+DyuFC7ueXpsrCrquoCfgF8BOwAXlNVdVtX9xsMLWKP1XekeX+BF2Q0TpF2DjU1dHiasliwtzqVgVTyEgtZefHfve/PvdBfNTVh18RKmxiiZcJM1ISO18r1HSoPkaNvzRbKyxOClpPTc/XWoyEtzf8J5EgIu/YUM3Ro8JTPwrYZlianlXq/P98o3/f8Fj4/l3tWTA57vAsugInKVgDePDSb2p3iAq/eJ26AnBHte7KzBouDVJdFF5XYbEemIqNvWQ0tqyaQ6uqOeffdTXfbMBAjj11V1fdVVR2tqupIVVXvi8U+g2EyCUHtbNYIiO21iFrrBPUVdo+nD+TPqirbmoYwa9ABJq97lgHXXsBPeJ5z4j9h6rn+j+SasGuRc2fqnvvaKR2dj9RsFuKYlwcjR8KQIZG36Wl8C3gdiUbIYBD7TUoK3jAOHSguzPh4PYLxLdfQ0d8w6ZLzWPr1BD58dBs2Evnwid0AHC4TgpM9uv3EKdltnn1taQvH3PwDknaFH+Rktx8ZYbfbRfAE4t7V7k2t0iQIwfftaO1tdLcNA31s5KkWqXTlAmpqElaLb2aK7/565Wi3gEeUmrfWsF8tZOpEB2p8Aq7cwfzqteP508fHYjD6+yLx8UJYwz3+RyI5WRf3jkwBB6JhSU4WFkxvqLkeDVlZeh31IzUBSFKSPlo3kNmniIMWDhQqplloGtGOFfCiKKAoZJ84gWkJ29i0S4T/leUuDLjJLmz/2JA5LJ0EbJR/X8crXwxl7/X3A2Bsqqfgsd9gtPhbNB7PkSvc5WuLahG7zSaidIdDX+4bkGkWTm/A4e5+UenFD8Xt0ewEh6NzmRUej2j9A/NjfYW9t1wMGnFVFSyfv4zWwvEseGIOpqw0/vNCI0m0cPx1erVF24j2Y8KMRj3S03xwo1EIbEdSLBMSRDTU2trx791o7J056+HQbA/NPjoSaA1HMGHP/NEZ/CthN0NPEJ2jgbMwdfSpyZf85AZ2WQswNtZRu6uGwVQQF9f+MUrNzqGIzSzfN5OnOQta4P1vyhh6YA3X//skxqx+n0v/+2O/bTrbl6Sq0ReD06JyrT78jh3+Za21Rq+6Wp9svKfps1ZMd9HViL2+XozuDCw56puX3tsi9uL3dvGA7Tc8vutslt2wEsVu4/PDYzlh8D6yjgk/92hycvAIvaM55Fp0qb3uKL25szQcee37p2NGZttPF+oJaOz5ozEPEF92YITelTEAuclWDjsGoO7cybucy5ik4CNM1UQzs0dWUafqA91K3/2e7x/9hHc4n4fLf0xDtf/N0uG+r9YWHjrzU06bbaW6XK/bMfT+m0jYvzPoJlpAommAr43qex/bbL1n5GpPROx9Sti1C7qzj3yhRrJpve7Q+4R9+zf6I8TWymzsGzZTzCiKpkbeNpRodNSO0ayY5OTOVWTsjYORegu+VSFDEfj9hfo+o7G6ctLsNKjp7N7QSC3Z/GDRiJDrHn+F/7Lyd7/hspYXvX8Xf9G1rIGad77iP3Wn0+RKZstbYgBVQulunngzn5VXvBh0G03Yg2XgaMKu9ZP1lpru0mOPgCZUHf3BPB597sZg+FZT7E1WTOP3B/h+r5nBlDOPD6g15LB3pSjEMeqkyIN7QwlAR4Q9PV0vYpWe3rvqp/cHQtkxvgQuDzXnal6ensUT6jfOyRSKuGGz2GnhjNCTjA89UZ996xilhD97bvVbfmBr5EJh4di+Tk+93LpJRFTl2xv4E7/nZtufg24TGLH7EjiFpN3eO+ZIlhF7BLQLOtpHLC0Kb2kRg4TCibbW895bWvn6N1dx2jVDef3wiUzNLqNwGBTbhrBzi7iiR82MPMonlG3SEWHXLAOTqXf4lf2R+Hj9ew5E6xOJ9F5iovDitf6MtLTgfRvZWeKm2LA5jngc5A8N/TgVl57MI+n3sPy0ZxmfXIoD4ce9fMISMqmjbH/XEtd3HzBjVlo4JX4tm/YNQHE6qCjRb+5ggZgm1MHuZe0931IHvSG/XQp7BLTIJZriQBaLEPPSUuGrR/L/GhrEv2AXjOlwGcXLNvDCjV+jerpnoMFHK/QW5qKzLIzIs9OimvmgYjIjkg5F7JA0mUKLhW/EFy5azMjw98dltH5kSEwU33MwmytUHn2wOVYNBt2P16yzQAbkiYybzzmJkYnlERv5k/53FyMfupHROSK6zjbWMeGx6xlp2M/ByjCNQu1hSn5yH/VvrQq5TkltBiOTKjijYCfbm4fy3OwXeOFfes/wvm3tZ+nQhD2YZar1nXU2y+1IRPdujxu3p/vzMPuUsGuCFE0rXF0tWvyamuiqyWnF/9vR2sILP3iDyx6ZydObjuXQziM/P55ia2Hz9ngms4XdT37E5F+dzigxAN8mTYMAACAASURBVJFNjslMyo38gcJ54b43c7govN2kEJIjQmKiEOJoJskO9b7Wua011GZziHIOo3TrZWiGv5WSHdqVYaSYwQ+nagJFYbC5gaqm0D3pxnWfs+C7/+OiB2YGX8HhYId1KCOzm7jwSnGxPsNP+YITvKvs+7b9qFe3O/xYk8CKrNHmtquqSJ+M1RO7lmdvd/eMt9unhF27aJuawldYLC+P3fyk1ZsO8AQ3e//etqomNjsOQ8Vra/mUMzjhmEM0zT4LgPGz0jidT8TrUZE7Y8Jlr2gioNV0CSQ9XTSifS1Nsa+SnS2u7WBPT50R9oyM4BN+ACSPymc6GxlIJTf8+Ri/ZVqxuGBMmJtFJnVcOUZU4ctOtVFlD50Lunmb2FGTJwW3GxRbK4n7dgBg3LmN9efdSxkFFE0zYD3vMi40rPBu++rxj5NNNauXtw9g3O72JTx8sVr9I/ZoyxwcPCie2Pfujc1AK4tFPC30hA0DfVTYXa7Qdow2P2ms+L5NyG/iaYy42LvzyLfA69eIi+GHj5/ifa955iksvvogt419hzNuDT+EHMIPYtEiw8TE4DdySoqYDagr+dKS6NH88mDfdyhhN5n8I3xt28REIexGo3gdaJ+5B+Xzv+PvYs2T31E43v+xLj5eL94WiPnsk9j4wKdc84Ko1Z+d4aTGPcBPBLPfeJbsN5cAUF6p7+TAXhcDH7iF+B+dj7Gxjvfv+Iyf194LwLEXiZHSd74xjbd+9Ar7bn6cUX+5iWt4iS8sk9tN0dfaGj5oCxyYFI2wu93iCd/hENvGokSItwPXJSP2iGhWjMslKigG+meqKmY6jxWu6jr++s4IRhv3cM2nVzCeHTz35SSaV22I3UHaHdTFp9sGMy5pH1mDfFRXUUhZdC0L/nUeadnhFTc9PXy0nZCg11Q3mdrf/AkJvW/auqOBjkTseXn+Vpm2rcmk94sYjXppBC8GA3ueeN/7JBh4/NTUEEGBwUD9GZegxotrTztGXbXudRQ+cBOm+++mbHsjh2v0zoGS9dUs/eIYRlPMlrdK+KhK5Or+Ke9JcscJ/0cZMpQhv/0xtVfegpqQSM7pRbiJo3a/fwTX2hr+aV2rvKoRjRUTuL/KyvYTbXcUrXGRVkwU+EbsHk9777yqKrZ56P9bvJJyz2Du/JUNY0YqP1DElGWP3tU1n0d1e9j+/Dr2fVHu/SzNpTWUfLqPtxa+xVeOaVx0Zuenjxk6NHzxKm1iZS3i87VmDIbO5apLuk5HIvaEBP13MhpD/96hOtCDHUebjCUaCy47V5xYQ6mugP/lfI5hDxdelcqKkskMpJJ47BSvLmdZw9kAvP/PGj53zOba8V8x7+1fhNz/oCHiAx0O6NOKJNQej79NE03EHqyhOHAg+uw7p7P9XMjeMgc9FLH3qZICmrBr4l1VJR47U1KEtxbrCm+f7R7MmMT9TLx0IgA3LIrn308c4Gt7Uaf3Wf75Hh65rZLP3XonkYIHhUw8ZAPDOTlnK+f+LrLdEoxQ9kogvhUfExLExZ2fLyIxmf3SM2gleX2FKVzWivZUFe731rz7SAHPwIHi/4QE0RhEmrQ9s22C7YbSRpgtwvf7uZNcqqhkEPucBRybtIWBrYd5YfMM73bLG4XATx7VEvZCGz5WRB3bNjuYfF74cwlHNH55sO/G6RSzQLV74gnC4cPCysnMFPdfQ4NoFNLSwOYK83hxBOlTEbvWIejroR08CFu36nMqxgq18hAbLGOZPli/wuuu+jULj9tJpSeXxgqfJFtVJW3dhzS+8xmPzfuIfav2tT8ZVWXXH17h4l8PZYN7Glenvcni0cu5ZcwH/CBpJb+cuoZnF/yPZy//jPtfHdHpqoLRTmKh1fYGPTUuO1uKek+iKEJIfH+DcKN2TSbx20Ua2RvJVjMY/CcXT0yMPIo1Y6jYqeutd8j88BVQVQ4zkJkj67kp9y1x3EQ3P8kWnaKXZnzAsYiO15EJB5n9mzlh9586ZSTj2c6Wr1rI+u+LFJ2UhmLreMpKNNVgQzV6hw/rEf+hQ6FLhttsojE+cEBk7O3fL973qJ4e6zztUxG79ujp+4jUlXoQ1j2HWPbgfibPTuHYn0zyvu8sq+SmCw5TTxEnneyvdIWT0uErqPhsN5nnjcSTkk7LB59x5V1D2YEoxLXrjzuYd+8TbFKncdXv8hngrubJ31fxLj9mTMI+nn4xgbQx8zt/4mGINpPFZNKjwdTUI1fsStIxcnPFY7zWQRhJtJOTIwc06emh5+8F0Zj4WjmKIhqDpjADSzOGiRoyN5fczNOLd/De+l9Rxf0MSK9m4vQseA0aHEmc9MaNbFz3XzjvfCx3P8yt7xl57DfVGM1nhz1nV/YgRmTuYU/DYPY+9R6nWjdz93NfMmHRaeE/bBCczvDWZKjUSVUVTkB6uj6Dk6qKyWV896c9FQT6+x6154a99qmIPZiwdxbF2szNl1XxzObZPPSsf5i7/b97+I4iLje+yvRr/W2X4SeJGZrvfySRf5+xFJelhU9eKPWKej5lbGwex5+abuYDy1wuu3MEZ/1+Fu9yHvOSP+eppyBtzOCufwAffIf6R5vJYjbrj/ChRilKegbtqUtRIhddM5sji3+48QjalImBaNZMSAYOZHDbDJg7GcfIdx6nFTOZGSpDJ4oDFhoP4soeBOf/EBSF9N/+lDdvXYf5h2dG2LkgP93KQcdAXnJfzT5G8MI7uSHXjd+6KWREr02cE6rTNZxNpdWRamkRjW1jY/tc9yNRh76rHLXCXvvhBjarUwA47MnxPq55tm7j/14aiRkrP/10frtRI+ljBpFPGTsYzyPOm7nklGruLb2aubk7+fz3H/OfJfXM4itOT/+alX/dwnD2kkUNL9/yDfeuPIHMqcMDT6XTZGQIIc/L01Pboq2XnpLSu2cyOprRGtmsrMiWSGpq5M5uozG0+JvNwa+DtLTwDYLHnMInk3/NsgUr/N7PHGBg0LQ83uIC/nLC2+22qb70l1FPS5Wf48CqJvNh02wAPq+bRNU/Pmi3nvX7vUxeOJ3Vt6xotwyEsNtsIpPOt/9CVcXf4YS9uVlYMtp2qurfQGgT9/Q2+pSwJyUJ8QqX7hQt29eJZ92rj1mLjUQa64Syf/HA51SSx9UD3iExNYhKKgoPX/YNFxRs5PaZqzhAIQCXXezA/MMziZ82iefeyePBj6eTNmcyy16L462XGphwxTQUY2y/7kGDRAZMcrK4Qc1m6ZH3B0wmoX3RzPmu5a1HIjU1+CjjcNF+2ElGFIXWF5cz6rYLePknn3nfzsyJwzloKGOW/o7GOx+KfGJhmDtV+Bo15DAWMbhpyRPi5h/1szPIe/ZujPuKefKaTQC8VTwp6H5UFbZvbx9tV1cL0Q5nZdXXt0/KaGrqXDVYm8vGXavu4kBjDHOyQ9CnhF17NO1KxG7533rq3lrFd8Vm0hQLM8eKTtDGMmEorisrJN94iBuWnRpyH6N+cz6LV8zgkqdP4ZlrvuL3RW8z/bKx3uWOvEJvVBI3YijmSceE2lWXSEzUvfH09L5b91zSHq3MQKwYPlx46YF9KZE6Z6Nhwg16hlf6INEatEychZrY0ame/Em76kLv6xsvbWIamyiJHweqSuvXW3jyuQQ++NteXudHACQmtIXVHg+ZH70KLj2U1oTYVzsaGiJn/wSjoUEkbUDHhP3Vra/y4uYX2Ve/r+MH7SB97mG8S8Kuqlzx23zKKSBXqWJG9j6ycoUAN5Y3w6RUSq3ZjMxrxp01Kqpdzvj5ccyIvFrM0aa804g2X1nSN8jIODJPX4MG+XeKhrPuoq6j73MhZuR1Tcx9UROTOIZi9jCKMZfPYPqX/+O1A8dgsFo4mw/4lmmwGpJo4ZSUDexuEgVt4t99kz/eY2feyn8w+qGf+O2zpUXYKtok2Z3NpKupgYKC6IW90dbIy9+9zGnDT+OkYSd17qAdoE9F7CDshs4Ke82bn1GOGMJcpeZy+klO0geK3saGShum2kr2MYz83F5omgUgOzv7N9HYMJ0hMPUxnLD/f3tnHidVdeXx76mlq7qq96Y3uulmU1kECYJC4qhBNGoSNEbjFhWXEDQzTmISlRhJnM/EiZOJo475RBwTkxh1NEYibnHBJYkKiorI3qDI2jQgm6DQwJ0/7ntWdXdVdVV3dS3t+X4+/elX772qe7rrvt8779xzz+3Oeq8lDentmLO+t4y/TPkt1XVeCsPCHhPCt20zS4gsCzmxbDmHl7TwwSc1HDwIb7xh+AOXcMHcy9nV0r5C5O7dNrSyalXPqjkaY+Pvyc5QvXfhvXy0/yNmHDej+42mwGdG2L1bN/Ot/7ArwhzOCs5kNl/87hgqB1gPo3XjQdrWt7CFavrX51agOpbn1mXWgpLXJDm+mDIi7QdLE2VRuaWAKyvbx/Fj2XYH/8JQmimoTG/ebOmFU2iYeZm1J3CIA/jZvXYb+wlwfuFsxvM6Z03aTmP1J7ThZ/Nm+CCq3O/789oXftm3L5Il01PcwmFd0fJRCw8vfZgvH/ZlhvUb1vUb0kCPhF1EzhGRJSJySEQyEpEoL0+9joN8spe1j77BOhqZ1v8Jnp3xEjf/rj/+oJdQbSk1tLBmg4/WlfaDawb20tL03aS6uv2F5U5MUZTuED0xLZHHLmJj8uFw+9h8fX3n0N/J917Iy9MeSPmO5PEk35cLgzZusrnZjosddelY7nplNBN+dBL1g+wdauPK3axaV0BAbHL6qsWdvcB0TWRMNonjrgV3ATB93PT0NJwEPY2xLwbOAmalwZakqPZ/yKK1RUCSeX3AobPO5v7WSwA4886T2NoY6UltlbUM4y2WrxnK5tW2w9QekeT0zQwRDtvYqHtD0wJdSk9wPfZkljosKeksvMGgnc0ZXWVxz6gJ7Bk1IWVb3AVG3NmaiQg6Txcb37NiXVJfjAlYJ6xuZDnMhs0L1rHUDOfz/dewdGMpzc0pm5RWmrc182Tzk3xz9DepLcrcEmQ98tiNMcuMMSvSZUyXbN1K09//yPZNKcRijOGC1tt4mHMBqBrQvpceChczqX45i7bV87c37LHaYbkl7G6usetpaZEupSdE1wjqCjdP3o2319RYoU/XE2M43NlRcevTdyQYsnehtevs74pBkUT7slGNBPiEjct2spxhDG5sY4x/CQuay5Kuyd4b3PnGnRQVFDH1qKkZbTdjMXYRmSYiC0RkwZbuBrnmzaOGzew0pXgWvR3/PGO4/6w/M3Hcfq47YR7N2OWHfjPh7pgeygk3HIeXA9y74Uv097VSUZkbMfZg0F6E7gXo/lZhV3qCK+zJDI66E998PvtTVmZ/p6tWfyhkPys67h8MxrbNvZk0b7Zxof4NkbBPW+MQBvMeS9eE+IRCaqoNZ9e9wpr9/Xn/hfeSssX3YRoKsUfx13ee45V1rzB1zFRKg5ldjqxLYReR50VkcYyfM1JpyBhztzFmnDFmXFV3h/yXLmUYywHY9Is/xj1t50tv8d9rv04bBczda2et3XfFS4z55cUxzy86ZiSnFv0DgPHDY6ygmyUaG2FkZPD/04tJ4+tKT3C94VT7UXRMPl3C7n5O9LJ8gUBs24IhK1fNW8up8GxvVzfe+AsYFGzh9Z1HAFBe5afyTLsoyK71XZfA3vbUPL5+yi5+cW561lpY+8PbuerxGykzVZw78ty0fGYqdBljN8ZMzoQhSbFyJSdVLSKwrY1nmgdxsTGdgoRt61u47d8/wsNBZp/9ADtXtfL+piDDvnUlJsHMz+/8cghjH3+LSd8fG/N4QYF9ZExUTCmdiFjPPHpKeWmpTbHqrYwJ5bNBKNS58Fey73O9fdeLTxTmCIUSryEaCET6d79+NmPFHbCNlYpYGLYnL2YUn6/7AGg/gttUvpN9m6yrX1oToKjS3i12bEtc3vGNHz/GlX+1furq1XDiyx8x/oTuD2R59uzmhpcmUbVjL1d/fxdBX+aTMfJrgtKsWXzynVa+cM5G5qz+Mj/43S20Xnp95PiBA9w1fSFP7DyFa4+fT/31F1EPTnmuxFQfPYAzjh4Q81hFhf05cCBzwl5R0blOSHl58ovzKko8ioq6NwDfv3/nksKJhL283J6/Z0/s49FzMQIBGDDA3ixE7LGOtenLyiKNn/HPDZ0+7+ghO8CZSVrSVE5J2Bq3Y3viNJinn4/MxDqC5Vz/w3qeebX7tZQ+XrGWZWYUP2lq49QvprGWeAr0NN3xayKyHpgIPCkiz6THrDh4vZjaOo47NcwaBvHsH9sXelh188P8vuVUph39JufcOjFtzRYVWW85kyGQeDNJE60kryi9SUdHo6uZqQUF8a+ZqqrOJTCil2n0+zuPJdXXRVT+pJM7P26MvDySlVN7RBmByiIK2Mf2XfHHzKRtP+sP1FJFK+/VHMvMUbPZeaiYm6e8hnd399bH2/DuNgAGjsxezLSnWTGzjTENxpiAMabGGNN5IcU0IwKnnN+PMfWt/HDnjRQ+ch8A3o1ruf3JoTQWtHDxHePT1l4wGOmA3ZmJ1x28Xp1ZquQ+XVUS9fvbr5/qbns8djp+V08NHT/fG7IB+UZiF9HyjxrOkZ4lHOlfQVERmHAx/djKjl3xXW//hjWsZjCnD3uPHQ89y6ifns3V3M6c1olc88W3EhsIBNY2c/DSy/m301+jdZld+H7tu7Zmw8BjszeLMO9mnopYgb16RhHbqeB/btnDwJkXM+eiR3jt4DFcelFb2gZ2wJbEdT0VjweGDEnfZ8ejpqbrcq2Kkm1ihSrKyiJPmx2F3V1mzp3R2hUdr2Pj9fEh5bx82s1x33PP3EHc86wNqR4qDNOPrWz/KP6jxd73W9hJGbVH9uNQUSn7mw7jnDtPBOBFJrF/n12CqfRvj7crKvYpTz/F8e/+ijmtE3ns9jV4d2xj3muHqPJso2l09lavyTv5cB/VRk8Icd4X1nEXV9LvqT/wk53X8NWjPuAr02PHybvbVseSqGVlPffcRewNo7bWPpK6y9LV10NTkz2mKLlOx1CM6/g0NNg+7gq732/7dEmJvRkk279DHeqJ7Rk1ge0/+zUf3nh73Pf4ikP4ip03ilDp2cH2vfE9vU0rbcZM/6GRxgonHMWfhtqxu4enz6Xghb9y7TVtvDD9YcKLXms30PVicwP7CFLCTmYtGMekycIj+6ZwxoiVWS2hnXfCHn2nv3xmPcVeOzJz1uhmbpjVlNZ/ZjzPoiex9rIym8LYv78V8sZGGDbMCnptrcbQlfyho8fuvi4osF6712sFfsQI29/B9vlEC3hE0ykcKcL2L52HKUj+kbwy8BHb9xSAMRQteKlTPYGNa6wXXjui/TJSh91xNafxFH9aMpJXn9rBbM7i2oUXMvayo3jq6shiHy2brUDMPuE2ALZTwQBZzxX/m/os3HSS18JeXulh9jNh5s+HH/32sLSuCOTzxS+0FctjLyiwto0cGUkj65hOVlkJgwenLwdYUbJJQUHEM4f2Qh/9pBu9PxXHxeNpf7PoDpVspfnjBvpdcxG7p/+AivvvaHd8/UbrCfYf1P6iPlDdn7HjfWw4VMcryyOi/zEhZr0ZGcPbvMNPhWcHZTdfxxNVU7mGX/Lzy1bi82d3kmNeCzvYDtTdvG6fz3rLsbzymprYa0FCZLQ++umgsRGGD48MttbVwahR1jsJhWw7AwfqCkdK3yEctn3ejZ1HC3i6Fkd3naAhQ7p37Uz++HEAvvn3bzOeBdz4wPBPjwVXvMOaDX5qvFtiPoU3DLAN3rflNI4pWsJMbgJgw4EaPt5rPf+W3UXUBT/EBILUPv07LljwfY64Mv4iPZki74W9J5SXR5aV60ii5caKiqzwDxli7QkGrYC7nnxDg3309Hqth37YYVoGQOl7uCualZfbCqTRwp6uSXTuJCZ36cd4xGvvhK9V8lXm8A/sLNRXt1th965ewQUXwoM7TmdgOHYpgcHDIoMI9QVbOPPla/j15EcAWD/Xlsja+HEFdcW5M1vdJe+EPZ2zLt10q46DQE1NiQdIvV4r3qWlcPjhdmJFNNGeRfTjpKL0RTyeyMBouikvj4xpdbzOXNxkhFisve5XXP3oiZx5vM0tX9dWy769B1ly/0IWYRezLwvFnmVVP6aauUzi4sBDXPGzQRwKF9M0yj6KzLzJw7731rHhYA21Zfu6/wf2Enkn7On02F0vuqOwJzu4435Guh47FSVfCQRSWEovBcLhyHUaDse+edTXJ4jd+3yUN5bw41srufPrL3AAPyuv/w23zLGe+wRe44aT58d86yeDR1D1l3u4+h/foGa8XbR+4JhyZjGNZg7ngWvfoZUaqqt6sBRTL/GZE3av13oA0VUTOz5C9kYHVZS+TCDQ/QHORPj9Nszj4rbh89kxK7DZM15v108MQ4+1g2aXvjqNdxnNMaF3eei6twlfNTXue/Y3DG73CL535Hj+6b5v8/nwQmatsfMxq4bHvqsI2RtQ+8wJe2mpjY9HF5iMFnLNWFGU1HGL1vUG0deke/0HAja5obg4EnuPPq+6urM9VccO/XR7QsUKZt5axpZzrsL4U7sj7R1+NMcNbeEAVjg+d+bAmOd5PV5C/vQt7p0KeRf97YmwuymMHQdhooW9N7wORfkskAmnyL3+3boyhx8eORYMRgqOlZba6z26AJkJF/HzMQ+yYUeYSx6Z0iM7TjrFy/3vrOecievpVxM/Z700WMretgQlLnuJvBN2d9JDd9YtrKyMPbIePVCqwq4ouUu0sHck+jouLIxd+nfyPeenxY7Kb0zmuSMXsHfEsQnPKwmUsGn3prS0mQp5F4qB7o++xys6FF0XWkMxipK7uNdpLA1whd0dJ0uUHtljRNg7cnyXyfUhfwjJwuSVvBT2VAc3fT4bc0tUTc7NW1dhV5TcxRX2mGuiBtv/LijIfiKERzwU+jJfvjcvhT1Zj9398kOhSBH/eLiDqb16l1cUpUe481hiCbu7L1rM480ezyTZGEDNuxg7WM97167E51RU2PzWAwdg//7kPrO+Pvt3eEVR4pPIY/d47PUbfQ2nErYNBKxWdGf8LhbuTajQrx57UlRWJg5tidjQS0GB9cATlQeIprY2PfYpitI7uOIdL2QaCLQX82SF3e+3138qkxO7+jy3iKCGYpLEFex4DB2qtVkUpS/i99sn63iOXVFR94R90CBbvmDQoPQsgVlUFHm6yMZi1nkp7BD/zlpdrVP8FaWvEgwmvr7daqouydSW8vkis1c9Hlt3pqeJLNFPFH6vH68njUWukiBvhb2urrPX7vFECvoritL3cFdkikc43D77LRmPvaOHXl5ukylSEfeO53b8zEx77Xkr7GBFvGNN9HRWf1QUJbfoSmw7Hk9G2GPF6xsakg/nhkI2hBNNx9TqgDezedQ9EnYR+YWILBeRRSIyW0SSHKZMD6Wl9u4aCNg7bGVlJltXFCXX6a6wuwkYLolKmVRVWR1yl/Krq+uctRPw5ZGwA88BRxpjRgMrgRk9Nyk1Bgywy9E1Nma6ZUVRch2RrsU9XhkR12kE68EXFibOn29qsnoUqzZ8XnnsxphnjTFulfp5QEPPTUoNn0+Xm1MUJT5dzU1JVB/KjZWHw7ZM8JAhnT/PfR0IWC8/lh4VeDNbhCqdMfbLgKe7PEtRFCWDNDjuZrxsmkTC7taRKiyMLM+XqDpsPDIt7F1GoETkeSDW1J0bjDGPOefcABwA7k/wOdOAaQCNGjdRFCVDlJTYmejV1bB3r52NDpGn/UTC7q6LHO2Fd6wPn0wcP+eE3RgzOdFxEZkKfAU4yZj4k3GNMXcDdwOMGzcuTZN2FUVRusbNWgmFIuVI3NnpiYhems+lshK2bYODB5Mv8y0i+L2Zq1fS06yYU4FrgSnGmMxXk1cURUmB6Pzy0tLuZdKFQjB4sN1OZf2GTHrtPS0CdicQAJ5zag7PM8ZM77FViqIovUA4bMW4X7+eVXItKbEhmGAK8478nsx57D0SdmPM0K7PUhRFyQ2Ki21qdDqKfRUX2xh8suSTx64oipI3+Hzpq+A4cGBqazD7vX4OHjqYnsa7IK9LCiiKomSLVEQdMuuxq7AriqJkgEzG2FXYFUVRMkDepDsqiqIoyaEeu6IoSh/D6/HikcxIrgq7oihKhshUOEaFXVEUJUP4PJnJMFdhVxRFyRCZirOrsCuKomQI9dgVRVH6GBpjVxRF6WNoVoyiKIrSLVTYFUVR+hgq7IqiKH0MFXZFUZQ+hgq7oihKH0OFXVEUpY+hwq4oitLHUGFXFEXpY6iwK4qi9DHEGJP5RkW2AB908+39gK1pNCed5Kptalfq5Kptalfq5Kpt3bGryRhT1dVJWRH2niAiC4wx47JtRyxy1Ta1K3Vy1Ta1K3Vy1bbetEtDMYqiKH0MFXZFUZQ+Rj4K+93ZNiABuWqb2pU6uWqb2pU6uWpbr9mVdzF2RVEUJTH56LEriqIoCcgrYReRU0VkhYisEpHrM9z2b0WkVUQWR+2rEJHnRKTZ+V3u7BcRucOxc5GIjO1FuwaIyIsislRElojIv+aQbUEReV1E3nFsu8nZP0hE5js2PCQiBc7+gPN6lXN8YG/Z5rTnFZG3ReSJXLFLRNaIyLsislBEFjj7sv5dOu2VicgjIrJcRJaJyMRs2yYiRzj/K/dnl4h8N9t2OW19z+n3i0XkQed6yEwfM8bkxQ/gBVYDg4EC4B1gRAbbPx4YCyyO2vefwPXO9vXALc726cDTgAATgPm9aFcdMNbZLgZWAiNyxDYBipxtPzDfafNh4Dxn/13Alc72VcBdzvZ5wEO9/J1eAzwAPOG8zrpdwBqgX4d9Wf8unfZ+D1zhbBcAZblim9OmF2gBmrJtF1APvA8URvWtqZnqY736j07zP2oi8EzU6xnAjAzbMJD2wr4CqHO264AVzvYs4PxY52XAxseAk3PNNiAEvAUci52U4ev4vQLPABOdbZ9znvSSPQ3AXGAS8IRzoeeCK9RSpAAAAv5JREFUXWvoLOxZ/y6BUkeoJNdsi2rjFOCVXLALK+zrgAqnzzwBfClTfSyfQjHuP8plvbMvm9QYYzY52y1AjbOdFVudx7fPYT3jnLDNCXcsBFqB57BPXTuMMQditP+pbc7xnUBlL5l2G3AtcMh5XZkjdhngWRF5U0SmOfty4bscBGwB7nXCV/eISDhHbHM5D3jQ2c6qXcaYDcB/AWuBTdg+8yYZ6mP5JOw5jbG32qylGIlIEfBn4LvGmF3Rx7JpmzHmoDFmDNZDPgYYlg07ohGRrwCtxpg3s21LDI4zxowFTgO+IyLHRx/M4nfpw4Yif22M+RywBxviyAXbcGLVU4A/dTyWDbucmP4Z2BtifyAMnJqp9vNJ2DcAA6JeNzj7sslmEakDcH63OvszaquI+LGifr8x5tFcss3FGLMDeBH7+FkmIr4Y7X9qm3O8FNjWC+Z8AZgiImuA/8OGY27PAbtcTw9jTCswG3szzIXvcj2w3hgz33n9CFboc8E2sDfCt4wxm53X2bZrMvC+MWaLMaYNeBTb7zLSx/JJ2N8ADnNGlQuwj11zsmzTHOASZ/sSbHzb3X+xMwI/AdgZ9ViYVkREgN8Ay4wxt+aYbVUiUuZsF2Jj/8uwAn92HNtcm88GXnC8rbRijJlhjGkwxgzE9qMXjDEXZtsuEQmLSLG7jY0ZLyYHvktjTAuwTkSOcHadBCzNBdsczicShnHbz6Zda4EJIhJyrlH3/5WZPtabgxm9MCBxOjbrYzVwQ4bbfhAbK2vDei+XY2Ngc4Fm4HmgwjlXgF85dr4LjOtFu47DPmYuAhY6P6fniG2jgbcd2xYDM539g4HXgVXYR+eAsz/ovF7lHB+cge/1RCJZMVm1y2n/HednidvHc+G7dNobAyxwvs+/AOW5YBs2zLENKI3alwt23QQsd/r+fUAgU31MZ54qiqL0MfIpFKMoiqIkgQq7oihKH0OFXVEUpY+hwq4oitLHUGFXFEXpY6iwK4qi9DFU2BVFUfoYKuyKoih9jP8H14VK1cucbIcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reconst_mean, reconst_std = reconstruct(mu_seq, S_seq)\n",
    "plot_reconstruction_forecasts(reconst_mean, reconst_std, forecasts_mean, forecasts_std)"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}