notebooks/source/time_series_forecasting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time Series Forecasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will demonstrate how to build a model for time series forecasting in NumPyro. Specifically, we will replicate the **Seasonal, Global Trend (SGT)** model from the [Rlgt: Bayesian Exponential Smoothing Models with Trend Modifications](https://cran.r-project.org/web/packages/Rlgt/index.html) package. The time series data that we will use for this tutorial is the **lynx** dataset, which contains annual numbers of lynx trappings from 1821 to 1934 in Canada."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q numpyro@git+https://github.com/pyro-ppl/numpyro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from IPython.display import set_matplotlib_formats\n",
    "\n",
    "import jax.numpy as jnp\n",
    "from jax import random\n",
    "\n",
    "import numpyro\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.contrib.control_flow import scan\n",
    "from numpyro.diagnostics import autocorrelation, hpdi\n",
    "from numpyro.infer import MCMC, NUTS, Predictive\n",
    "\n",
    "if \"NUMPYRO_SPHINXBUILD\" in os.environ:\n",
    "    set_matplotlib_formats(\"svg\")\n",
    "\n",
    "numpyro.set_host_device_count(4)\n",
    "assert numpyro.__version__.startswith(\"0.12.1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, lets import and take a look at the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of time series: 114\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAD4CAYAAADB2L5nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABhuklEQVR4nO29eZwcZ33n/3767pnuuUfn6LYOS7JlW0L4AAN2jE3iYCcE1iRgB0i8AbIbsksIbLJssll+YXMSEnDwQsAkgOOYy3EwxgiMwTaWJUuWLVnH6JzR3EfP9PR9PL8/qqq7NdN3V3XPqJ/36zWv6amumq6p6a7v870+XyGlRKFQKBQKxeLH1ugTUCgUCoVCUR7KaCsUCoVCsURQRluhUCgUiiWCMtoKhUKhUCwRlNFWKBQKhWKJ4Gj0CZSip6dHrl+/vtGnoVAoFApFXTh48OCElLI333OL3mivX7+eAwcONPo0FAqFQqGoC0KI84WeU+FxhUKhUCiWCMpoKxQKhUKxRFBGW6FQKBSKJYIy2gqFQqFQLBGU0VYoFAqFYolQ0mgLIbYKIQ7nfM0KIT4ihOgSQjwlhDilf+/MOeYTQoh+IcQJIcTtOdt3CyFe0Z/7rBBCWPWHKRQKhUJxuVHSaEspT0gpr5FSXgPsBsLAt4GPA/uklJuBffrPCCG2A/cAO4A7gM8LIez6r3sAuB/YrH/dYepfo1AoFArFZUyl4fFbgdNSyvPAXcBD+vaHgLv1x3cBD0spY1LKs0A/sFcIsRJok1I+L7V5oF/NOUahUCgUwGvDs7xwZrLRp6FYpFRqtO8BvqE/Xi6lHAbQvy/Tt68GBnKOGdS3rdYfz9++ACHE/UKIA0KIA+Pj4xWeokKhUCxdPvPDk/y3R15u9GkoFillG20hhAt4O/BvpXbNs00W2b5wo5QPSin3SCn39PbmVXJTKBSKy5JQLMXFQIRAON7oU1EsQirxtN8GvCSlHNV/HtVD3ujfx/Ttg8CanOP6gCF9e1+e7QqFQqHQiSZSABwbnm3wmSgWI5UY7XeTDY0DPAbcpz++D/huzvZ7hBBuIcQGtIKz/XoIPSiEuF6vGr835xiFQqFQANGkZrRfGw42+EwUi5GyBoYIIVqA24D/nLP508AjQogPABeAdwJIKY8KIR4BjgFJ4MNSypR+zAeBrwBe4An9S6FQKBQ60UQagGNDytNWLKQsoy2lDAPd87ZNolWT59v/U8Cn8mw/AOys/DQVCoWiOVDhcUUxlCKaQqFQLCIMo90/FiSeTDf4bBSLDWW0FYo6Ek2k+K2HDnBuItToU1EsUqKJND0+N4mUpH9srtGno1hkKKOtUNSRc5MhfvjaKPvPTTX6VBSLlGgixXVrOwAVIlcsRBlthaKOGEVGoViywWeiWIwkU2mSacmVK9vwOG2qGE2xAGW0FYo6YuQrldFW5COq57B9bgdbV7TxmvK0FfNQRluhqCMR3WjPxVIl9lQ0I8aizuO0sX2ln2PDs2ijGhQKDWW0FYo6ElOetqIIkbj2/nA77Wxf2cZMJMHQTLTBZ6VYTCijrVDUkaynrYy2YiGxpOFp29m+qg2A11ReW5GDMtoKRR0xCtGU0Vbkw3h/eBw2tq5oQwhVQa64FGW0FYo6YoQ/VXhckY9sTtuOz+1gXVeLqiBXXIIy2gpFHTGGQSijrciH4Wl7XXYAtq9q47URZbQVWZTRVijqSDSuctqKwmQ8bYdutFe2cX4yTDCaaORpKRYRymgrFHXE6MMNqZYvRR4iOS1fAFeu1IrRjo+oMZ0KDWW0FYo6onLaimLk5rQBtq7wA3BqVGmQKzSU0VYo6ohxU56LJ5VohmIBRiTGrXvaXa0uAGZVeFyho4y2QlFHjPCnlBCOqxC54lJi8zxtr9OOECoyo8iijLZCUUeM6mBQN2LFQoxIjFc32kIIWl0OVbioyKCMtkJRR4ybMqgKcsVCook0dpvAac/emlvddsKqcFGho4y2QlFHco22qiBXzCeSSOFxXHpbbnU7mIurBZ5CoyyjLYToEEI8KoQ4LoR4TQhxgxCiSwjxlBDilP69M2f/Twgh+oUQJ4QQt+ds3y2EeEV/7rNCCGHFH6VQLFYiiRQtunCG8rQV84kmUpl8toHP7VCpFEWGcj3tvwO+L6XcBuwCXgM+DuyTUm4G9uk/I4TYDtwD7ADuAD4vhDDehQ8A9wOb9a87TPo7FIolQTSRosfnBlROW7GQaCK9wGi3upTRVmQpabSFEG3AzcCXAKSUcSllALgLeEjf7SHgbv3xXcDDUsqYlPIs0A/sFUKsBNqklM9LrdflqznHKBRNQTSRpsentfGEVMhTMY9oMpVp9zJoddvV/HVFhnI87Y3AOPBlIcQhIcQXhRCtwHIp5TCA/n2Zvv9qYCDn+EF922r98fztCxBC3C+EOCCEODA+Pl7RH6RQLGZyPW0VHlfMJ5ZIZSrHDVpVeFyRQzlG2wFcBzwgpbwWCKGHwguQL08ti2xfuFHKB6WUe6SUe3p7e8s4xcXDz05N8OVnzzb6NBSLlEgiRbcKjysKkDc87nYQVlEZhU45RnsQGJRSvqD//CiaER/VQ97o38dy9l+Tc3wfMKRv78uz/bLiXw8M8Nc/OKnUrhQLkFISTaToanUCMBdVN2LFpUQSqYzuuIHPrfq0FVlKGm0p5QgwIITYqm+6FTgGPAbcp2+7D/iu/vgx4B4hhFsIsQGt4Gy/HkIPCiGu16vG78055rJhNpJgLpZkJqJkBxWXEk+lSUtocTlodak8pWIh0UQqM+HLoNXlIJpIk0ylCxylaCYcZe73X4CvCSFcwBngfWgG/xEhxAeAC8A7AaSUR4UQj6AZ9iTwYSmlcXf6IPAVwAs8oX9dVhgawQNTETpaXA0+G8ViwlBD8zjtKk+pyEu+lq9Wt/ZzKJ6i3aukNZqdsoy2lPIwsCfPU7cW2P9TwKfybD8A7Kzg/JYcs7qHPTgd5qq+9gafjWIxEc0Zu+hTghmKPBTKaYNWA9HudTbitBSLCLVsM5lZPU85MB1u8JkoFhu5utLK01bkI5ZcmNPONdoKhTLaJpP1tCMNPhPFYiOSM8FJqVwp8pHP0/blhMcVCmW0TSSaSBHT5+EOTClPW3EpRk7b8LRVIZpiPvmqx1tdytNWZFFG20SCOS08A8rTVswjontKbqcNn9uubsKKS0ik0qTScmH1uB4eV21fClBG21SMyvFev5vB6bDq1VZcQjSpctqKwkRz0ie5+FROW5GDMtomYuSzd6xqI5pIMzEXb/AZKRYT0filOe2gugkrcsi0BLouNdotRk5bvV8UKKNtKkbl+PaVbYDW9qVQGMz3tOPJNAklmKHQyXjajoWKaICqgVAAymibStbT1vqzVV5bkUskfqm4CijvSZEllswfHvc67dgESn9cASijbSpGTnv7KuVpKxaS26dttPGo4iKFQe6iLhchBK0upT+u0FBG20RmI9qHamW7h65WFwNTytNWZDH6tN1OW46nrUKeCo1oMquYNx9VuKgwUEbbRGajCVx2G26Hjb5Or/K0FZcQS6QQAtwOW06eUt2IFRqFqsdB0x9XCzwFKKNtKrORBG1eB0II1nS2KFU0xSVE9AlOQgjVxqNYQK74znxa1XhOhY4y2iYyG03S5tEE/fu6vFycjpBOq15thUY0kcart/OoQjTFfHIHysyn1aXC4woNZbRNZDaSwK9P4enrbCGeSjMWjDX4rBSLBc3T1j5yKjyumI9htN2O/J620h5XgDLapjIbTdDm0W7Gazq9gJr2pcgSTaQywhlKmlIxn2I5bSV7qzBQRttEtJx21tMG1falyBLVc9qgFRaBCo8rsmQU0VT1uKIIymibyEwkJ6dteNqq7Uuho41d1D5ybocdp10olStFhuKetipEU2goo20is1Gtehy0D54xOEShAC2n7c3RlVbekyKXaDKFwyZw2vN72rFkmqSSvW16lNE2iWgiRTyZznjaoOW1laetMMgNj4OqCFZcihaJWehlA7S4jHSKisw0O8pom4QhYWrktAHWdLUwGFCetkIjklOIBuD3qJCnIks0kcqbz4ac8ZxKf7zpKctoCyHOCSFeEUIcFkIc0Ld1CSGeEkKc0r935uz/CSFEvxDihBDi9pztu/Xf0y+E+KwQQpj/JzUGQ8K0Pcdo93V6GQpEVUhLAUAskb7U03Y71E1YkSGSSOVt9wLV16/IUomn/RYp5TVSyj36zx8H9kkpNwP79J8RQmwH7gF2AHcAnxdCGO/EB4D7gc361x21/wmLg4ynrbd8AazpbCGVlgzPRBt1WopFhJbTzn7kNJUrFe5UaMRyChXno/r6FQa1hMfvAh7SHz8E3J2z/WEpZUxKeRboB/YKIVYCbVLK56WUEvhqzjFLHmMsZ9slnrbR9qXy2oqFOW2f286cvthTKLTweClPWy3ymp1yjbYEfiCEOCiEuF/ftlxKOQygf1+mb18NDOQcO6hvW60/nr99AUKI+4UQB4QQB8bHx8s8xcYyG9VWwJcUonUpgRWFhpRyYfW4y6FuwooM0WQqr+44ZPv6laetcJTeBYCbpJRDQohlwFNCiONF9s2Xp5ZFti/cKOWDwIMAe/bsWRLi3VlPO3tJV7R7ABhV4fGmJ55KI+WlPbiq5UuRSzSRLmy0XSqnrdAoy9OWUg7p38eAbwN7gVE95I3+fUzffRBYk3N4HzCkb+/Ls/2yIJvTznrabocdl8PGnCo2anqyale54XGtEE3LFimanWLV40Z4PKzuJU1PSaMthGgVQviNx8BbgVeBx4D79N3uA76rP34MuEcI4RZCbEArONuvh9CDQojr9arxe3OOWfLMRpK4HLYFOSm/28FcVH3Qmp18E5xa3Q7SUitQUygiiRTuAp52thBNvVeanXLC48uBb+vdWQ7g61LK7wshXgQeEUJ8ALgAvBNASnlUCPEIcAxIAh+WUhrvtA8CXwG8wBP612WBNizEuWC7moOrgKzR9jovLUQDLU/Z4io3U6W4XJnfEpiLx2nDJlR4XFGG0ZZSngF25dk+Cdxa4JhPAZ/Ks/0AsLPy01z8aMNCFl5On8pbKsh605eExz05FcH+hpyWYhFRLDwuhFAOgAJQimimMRtN5vW0fR4HQRUeb3qMnHaup62KixS5RBOFq8dBOQAKDWW0TSJ3LGcuPqV6pQAicc3Tdud4UkowQ5FLNFlYexw0/XF1L1Eoo20SWk47f3hcFaIposmFOW2jIli9PxSJVJpUWhYMj4PhaatCtGZHGW2TmI0k83raKg+lAIjGF+a0W9UQCIVOvpqH+ai+fgUoo20aharH1SQnBeT3tFV4XGFgdBcUavkC5QAoNJTRNoHMLO0C1ePRRJqEmvTV1ETiC8VVDGlK5T0pYob4jqNEeFxFZZoeZbRNIJ8amoEaqaeA/H3aRvW4EsxQZN4frmKetl3ltBXKaJuBMUs7X07br0KgCrI5y9zqcZtNaBXB6r3R9GRkbguIq4C2yFP3EYUy2iaQb5a2gSGgoT5szU0skUIIcM8Lf6reWwVkax5KFaLFkyrV1uwoo20CM3lmaRuo8LgCNE/b47CjywFn8KniIh585jQ/PzPZ6NNoKJH4Qm36+WSGhqgQeVOjjLYJZMZy5lNE0z9oShWtuYkm0nlvyKqNB/72qVN88+Bgo0+joUTLaPnKaNWrYrSmRhltE5iNGjnt/NXjoMLjzU6kgERlq9ve1O+NSDxFJJHKpJialWjS6C4o7Wk3+yKv2VFG2wSKetoe9UFTGMMgFhptLTzevOHO6XAcyKaYmpVyPO1W5QAoUEbbFGajCdx5ZmmDCo8rNAoZ7WYPj0+FDKPdvNcAtEJFKGG01YAZBcpom0IhCVOAVld2ZrKieVE57fwYnvZs03vaC8V35pMV42neyIxCGW1TKDQsBMBht+F1ql7cZieSSOUVzmh12QnHm/cmnPW0m9toZ7THSyiigfK0mx1ltE2g0FhOA5/SH296onrL13w8TjvRZAopZQPOqvFM60Z7LpYk2cT9x9FECodN4LCXUYimqsebGmW0TWA2msxbhGbQ7MVGCr1PO4+n7XHakRLiTWqwpsJZD3u2ies+tPRJ4dA4qE4UhYYy2iYQLOVpux3MNXlLS7MTS6TzetqGQpqR02w2DE8bmjuvHU3mL1TMxe2wYbcJFR5vcso22kIIuxDikBDicf3nLiHEU0KIU/r3zpx9PyGE6BdCnBBC3J6zfbcQ4hX9uc+K+fJQS5RiOW1QvbgKI6e98ONm3KiN6uFmwyhEg+bOa2vdBcVvx0IYWvXN+V5RaFTiaf8e8FrOzx8H9kkpNwP79J8RQmwH7gF2AHcAnxdCGEvIB4D7gc361x01nf0iQEpZtHocwOd2qvB4k1Msp60936SedjiO3aat3ReL0T42NFt3bzZWRngclFa9okyjLYToA34J+GLO5ruAh/THDwF352x/WEoZk1KeBfqBvUKIlUCblPJ5qVXdfDXnmCVLLJkmnkoXzWn7PQ7mYovjhqSoP1LKgtXjhncVaVJPeyqUYE2nF1gcRjuVlvzqA8/yuR/31/V1I2V42qC3CKpCtKamXE/7M8DHgFx3YLmUchhA/75M374aGMjZb1Dftlp/PH/7AoQQ9wshDgghDoyPj5d5io0ho4aWR8LUoNVtZ66Ji2yanXgqjZT5e3AN7zvapEZ7OhRnfU8rsDiMdjieJJpI132ASaFIzHxaVVFr01PSaAsh7gTGpJQHy/yd+fLUssj2hRulfFBKuUdKuae3t7fMl20M2bGcxcPjKg/VvETjhYUzsuHx5nt/SCmZCsdZ372YjLb2f3jl4kxd/yeFFPPm43MrzYdmpxxP+ybg7UKIc8DDwC1CiH8BRvWQN/r3MX3/QWBNzvF9wJC+vS/P9iXNVEi70XS0FA+Px1NpYsnmuzErcmcl5ytE06vHk82X0w7HU8STaVa2e3A5bIuietww2omU5PBAoG6vW07LF2hSpspoNzcljbaU8hNSyj4p5Xq0ArMfSSnfAzwG3Kfvdh/wXf3xY8A9Qgi3EGIDWsHZfj2EHhRCXK9Xjd+bc8ySZXgmAsDKdk/BfTJSpipE3pQYs5LzTflqZk/bUEPrbHXR5nEuiklfuQbxwLmpur2u1vJV2ofyqZx201M4EVuaTwOPCCE+AFwA3gkgpTwqhHgEOAYkgQ9LKY070geBrwBe4An9a0kzFIgCsLLdW3Afnx46D8VSdPvqclqKRUTW0y5ciNaMRtto9+pqcdHudSyq8DjA/nPTdXvdcqvHW9yq5avZqchoSymfBp7WH08CtxbY71PAp/JsPwDsrPQkFzPDMxHaPI6MxGA+jOH1QVVB3pQU87TdDqNPu/nC47medrvXuSiMtuHF7ljVxkvnp0mlZaYlzUoqqR63WvPhwmSYv3nqBJ9+x9UlFxIXJsNcmArzhs09lp6TIotSRKuRoUCUVR2FvWzQCtFATedpVowebHfenLYeHm/CeoeMp72IjHZY/4y+aUsvc7Ekx0dm6/K65VaP+90O4klr62N+dHyU7xweKiun/5l9J/nIvx6y7FwUC1FGu0aGApHSRttjaAY3/qakqD9G6Dufp230bjdjeNwo4tTC44vDaBue9pu2aF0rL561Pq8tpSy7erxdF3GatXD++NCMlvI7MRIsue/psTmCqlanriijXSPDM5GiRWiQEx5Xb+6mxDDI+fu0m1d7fDqkqaH5PQ7NaIcbb7SNVMbm5X5WtXt48bz1ee1ESpKW5BXfmY+hvGhl0d5QQCuuLRVlkFJyejxELJkmnW7OKXWNQBntGojEU0yHEyo8rihKpIjRdthtOGyiOT3tcJzOFic2m6DN6yQYSzb85m942i0uO3vWd3Hg3JTlY1ON1Ii7yCxtA8NoWxmVyBrt4p72WDCWya83q6JfI1BGuwbKafcCFR5vdgwvOl94HPSZ2k3qaXe2uAAt7CslBBvcgxyOpbDbBG6Hjddt6GJ0NsbAVMTS1ywWiZmPIeJkZU+70RFzYiRYdBHVPzaXeayMdv1QRrsGhvXcTylPu8Wp+rSbmaynnf/j5nHamrIQbSoUp7NVM9qZsG+D89qheJIWpx0hBK9brw0ufNHifu1iinnzabfY006k0owGo6xo8xCOpxiYDhfc9/R4jtGON+b9K6WsSzRkMaGMdg1c1MNIq4r0aAPYbEKbqa3C401JKU/K7bA3ZXh8OhynK8fThsZLmYZjKVr0GpQty/y0eRwcOG+x0S6imDcfY8bBrEUOwOhsFCnhLdu0URLFQuSnF4Gn/fyZSX7tH5/n52fqJ4TTaJTRroFhPYy0vN1dcl9tprYKjzcj0UQKIQrnLD1OW5P2aScynvaiMdqJFK0uzTDabII967vYb3EFeWZRV0bLV7vFEQkjNP7mrb0IAceHixjt8VDmcaM87QuTWiTg0ED9hHAajTLaNTA8E6HH584IZBRDm4PbfN6UItuDq6n3LkTLaTfXe0NKSSAcp6tVM0KLxmjHkhlPG2D3uk5Oj4csPa9MzUMZ1eNuhx2P0zqddqMI7YplPtZ1tXBitHAF+enxuUw9T7hBRttIUR69WJ9++sWAMto1MDQTZVVH8SI0A5/H2fAiG0VjKKV25XHamy6nHYwlSaZlphCtHlXR5RCKJ2lxZdUNe/1aFM3KXHupmof5tHms62nPTfltW9FW0NOeiyUZnomyc3U70DidAaMY+JWLMw15/UagjHYNDAUiJfPZBj63nbkGDET40fFRHj+y5IepLWmiiXTBynHQbtaNCi82iulQVg0NrA/7lks4nqIlx+M1QuVWepJB/b5gtIaWot1r3XCV4ZkInS1OvC4721b6OTsZyvvePKMXoV2lG+1G5bQNT/vCVLhhff4/PjHGuYlQ6R1NQhntKpFSMhyIsLJcT7tB4fHP/PAUn913qu6vq8gSKaF25XE0X8tXru44aJPw7DbReE87lswYaiBjwMMWTtYyRJf8nvJGQbRZqB6XK8u8bYUfKeHU2EJv+/Q8o92o8PjITJQ2/bodHaq/t51OSz70Ly/x0PPn6vaaymhXyWw0SSieqsDTdlou9D+fdFpycjSYkYtUNIZYKaPdhOHx3AlfAEKIRSFlGpnnaRt5ZisjIYanXa7Rbvc6LZMxHQpEMhMLt61oA/IXo50eC2G3Cbau8AON9bSNSvdGhMgHpyNEEim2LvfX7TWV0a6SjLBK2Z62PfPhrBcXpsJEE2mmw/Gm6mNcbJTKabubsHrcWEgaOW1gURjtUDx1ycS+rKdtpdFOIgSXePjFaPNYN8b0YiDCav2etrarBa/Tnrft6/T4HOu6WjK1CNEGeNrBaIK5WJIdq9pY3eHl1aH6F6OdGNWuzWZltBc/w2XM0c7F53EQiqfqajyND1sqLS3r61SUZjaSxO8pnK9sxurx6Ux4PHtdrAz7lks4nrzE084YbQv/P8FoEp/bga3MEaBW5bSD0QTBaDITHrfZBFtW+PNqkJ8en2Njry9Tq9GI8PiIns9e0e5l5+o2Xm2Ap31SN9pblvvq9prKaFdJpsqyTE+71e0glZZ1zV0abyjI5hAV9Wda19guhKcJxVWmwnGcdk10yKDN42hoIVo8mSaRkvPC49r5RSzOabcVWdTNp83rZDaSMF2nPZ/C45Ur/BwfCV7ibCRTac5OhNi0rBW7TeBy2BoSHjemka1s97BzVTtnJ0J1j2aeHA2yusNbdFFuNspoV8nwTAS7TbDMX57R9us3p2AdBVZyR+spo904AuEEHTlh4PloMqbNFR43dMdze9cbHR4PZ4aF5ITHdU/SyiLSYDRRdj4btOuUltnhJmaRzxHZusLPVCjO+Fwss21gOkIiJbmiV/MuvU67pYuaQozkzH7Y2acVxB2tc4j8xEiwrl42KKNdNcMBTZ/XXmZIyxgaUs8K8uMjsxnxA2W0G0M8mWYulsy0NuXD47STSksSqeYx3FOh+IJrooV9G5fGMUK8re48hWgWh8crMdqGV272AmcoY7Sznna+YjRDvnTTMs1YtbjsDfG0h2eiCAHL/JqnDdQ1RJ5MpTkzHmJLHfPZoIx21QyVMUc7F6PIpF5DQ6KJFOcmw9ywsRvI5hAV9SUQ0XO3RcLjRl6wmULkWspgodGeiSQaVjSZz9N2O2zYhMUtX7FEReHV7HAVc89pOBBdED3cpleH50btjHavTT05nnYDCimHA1F6fG5cDhu9fjcr2jx1NdrnJsPEU2lltJcKwzNRVpaY7pVLdjxnfYz26fE5UmnJ9Zs0oz2pjHZDCOiCD6XC40BT9WoX8rRTaUmoQT2/RhQs19MWQtDqclhePV6Rp60PDbHC054fPexsdbG8zc3hgUBm2+nxOXp8btr1hajX1Zjw+PBslFU5jtPO1e11bfsyaoaMtrd6UdJoCyE8Qoj9QoiXhRBHhRB/qm/vEkI8JYQ4pX/vzDnmE0KIfiHECSHE7TnbdwshXtGf+6woJMa8yEmnJcOBS98wpfDrakf1MtrGyvjaNR14nLZMX6yivmSqpIsYbXdTetqJSyrHofH640aO2Ou81IBqRmnxGO2MepzJRVcXA5G8hbVv3b6C/3hlmD957CjJVJrT4yE29bZmntc87UZUj0dYcYnRbuPMRIhQHe+xQmg67fWkHE87BtwipdwFXAPcIYS4Hvg4sE9KuRnYp/+MEGI7cA+wA7gD+LwQwli6PgDcD2zWv+4w70+pH5OhOPFUuuQc7VyM1Xu9Jn2dGA3isttY39NKV4tL5bQbhLFY6ihWPa4b7ViTCKyk0vqwkDzhcaBhcpThPJ42aDlbqzxtKaVeiFZBeNyqnPZMJO897U/evoPfesMGvvLcOd7/0AFOjQYz+WywflFTiOFA9JKW26tWtyMlHBuuTzHaqbEg67paypqDbiYljbbUMAanOvUvCdwFPKRvfwi4W398F/CwlDImpTwL9AN7hRArgTYp5fNSS1p9NeeYJcVwTtViuWTD4/V5c58Y0T5YTruNzlZltBvFtG6AOosVojmaKzw+G0mQlguvSaOHhhi92C2u+Z62deHxmN5mVpGn3WK+Tns6LRmZieY12nab4I/v3M6nf/UqnuufYDaazFSOg+Zp17tPOxhNEIwlL7kHG5Kq9cpra5Xj9Q2NQ5k5bSGEXQhxGBgDnpJSvgAsl1IOA+jfl+m7rwYGcg4f1Let1h/P357v9e4XQhwQQhwYHx+v4M+pD8bM2Uo8baMftV6FaCdHgmzVWxG6lNFuGPPlOvPhabLw+FT40mEhBlaFfcslrIdV83nakYQ1n9vZjIRp+Z62z+VACHON9sRcjERKFk353bN3Lf/8gddz3doO3rS1N7Pd66q/zsDorCGskj3fZW0eev3uuuS1jULfeuezoUyjLaVMSSmvAfrQvOadRXbPl6eWRbbne70HpZR7pJR7ent78+3SUKrxtL1OOzZRn/D4TCTB0EyUrXq7hjLajSMQTuB22IrOSs4a7ebwtAvl+Ruf087vaVsZHjeGhbRV4GnbbML08ZwX87R75eOGTd1860M3sanBnvZQAUXKnavaeK3AOFEzOTMeIpWWi9fTNpBSBoCn0XLRo3rIG/37mL7bILAm57A+YEjf3pdn+5JjKBDB7bAV7b2djxCibpO+slWNWU9btXxpzEYTHDw/VbfXM0REipGtHm8STzuU39Nua/B4TsPTbpm3wPI67Zl8t9lUOuHLoM3rMLWnvZrooYG3AX3aIzlqaLms7PAypnvhVmJMPluURlsI0SuE6NAfe4FfAI4DjwH36bvdB3xXf/wYcI8Qwi2E2IBWcLZfD6EHhRDX61Xj9+Ycs6QY0nM/lRa/+9yOzIfUSozK8Yyn3eIiGEsSbzLVrflEEynu+6f9vOsLP6+bNzcdThQtQoMcT7tJCtECBfL8frcW9m2kp+2y23DaL70ttrjshC0KjwerCI+D+epx+YRVysXbAO18Q3J1edulRrvH52YqHCdpsVDRiZEgDptgQ09r6Z1NphxPeyXwYyHEEeBFtJz248CngduEEKeA2/SfkVIeBR4BjgHfBz4spTT+ox8EvohWnHYaeMLEv6VuDAcqE1Yx8HkcdWlHODESxO92ZPJTxs2xmdu+0mnJ7//rYQ5dCJBKS/rzzAi2gkAeEZH5eBzWj39cTEwVyPNbEfathEg8SYt7YRrD63JY9r+p2tP2OE2NSAzNRGh12SsK0xu0uOwkUvVV9BueiWSEVXLp9bmQMvses4qTo0E29rYueP16UPI/JKU8AlybZ/skcGuBYz4FfCrP9gNAsXz4kmAsGON167sqPs7ndtSlT/vESJAtK/yZSEC3brSnQvEFK9Nm4f9+/zhPvDrCb964nq88d46To3PsXlf5/7BSpsLxjKpUITLh8SaJhFycjuD3OPLm+du81o2dLEUonso7HrPV0px29Z52/9hc6R3LZCgQqSp6CNlIUSSRWhClsIrhmWhex6nH5wZgIhgvey5ENZwYDbKrr8Oy318MpYhWIVJKxoMxev3uio9tdTsIWmy0pZScGA1eUtXYmWO0m5F/+fl5vvDMGe69YR2fvHM7Xqf9kgloVlJqWAhkxVViTZLTPjVWuFWmkUND5o/lNDC0ta2QV63J0zaxyn4okL/dqxwy+ux1jBSNFDDaxn15ImfAidmE40kGpiJsbUA+G5TRrphgLEksmabXV7nR9tchPD46G2MmkrjEu+tuYqMdjCb4038/ypu39vLJO7djswmuWObj1Kh5Xkoh0gVERObTbIVo/WNzbC6gItXuNTfsWwmhWCqv0fa6HEhpTXX/bDSJEFobVyW0t5if0y53zPB8WhpgtAvNfjA87fGgdUbbuHdsVkZ7aWC8GarytF0Oy/u0jw5pPYq5nkwze9oDU9oYwXfuXoNDD91tXu6ri6cdjCZJy+JqaAAuuw0hmqPlayoUZ2IuXlD6sfGe9kLjaRglK4aGBKMJfC4HtjKnBRq0e51EE2lTVPTiyTSTNaTOvDnh8XowF0sSjCZZ0b4wMtBTB0/7RIM0xw2U0a6QWoy2z2N9TvuHr43R6rJzzZqOzLYOvZWmGY220VOf60VsWe5nLBizXC7TKPwrVYgmhMDjqH8FbiM4pd/wCnkpmtFuzHjOcDy1QFgFsuFfK/LaleqOGxgFY2ZM+poMafe0anPARk67Xr3aRrtXvshAq8uOx2mz1GifHpvD5bCxtqvFstcohjLaFVKL0fa7HYTiSdJpa0YPptOSp46N8uatyy7Rw3XYbXS0OJuyetxoZVmdk68zhtaftLiCPGO0W0sXGXmctqZo+TqlF08VCo+36eHxRoznDMdTRT1tKzzJSnXHDdpMVI8bm63+ngZZMZp6LTqNhfiKPJEBIQQ9PjcTc9bd6wYDEVZ3eC+ZhlZPlNGukIzRriKn3erWcmMhi8bYHRqYZmIuxlt3LF/wXFeLqynHc14MRHHaRSbXBbB5meblWR0iL2csp4HW63r5h8f7x+ZoddkLtky2eZzEU+mGXItQLJnX0zaMthX1KFV72iaqxxn3tGVVGu1MeLxOnvbwTH41NAPNaFvnadeS/zcDZbQrZHwuhtMuMpKLlWB80KwKkT95dBSnXfCWbcsWPLdYVNH+9qmTfPX5c3V7vaGANr4vN2e4usNLi8tueTFaueFx0EKMTREeHwtyxXJ/wdaiRkqZFvK0jVGdVhilao12u4nqceNztXnamfRBnd6/Rnh8eXv+8+3xuS0tRLs4HbkkcldvlNGukPFgjB6fu+LCEci2dZiRh5qPlJInj45w46aezOi+XBbDpK9QLMkDPznN40eG6/aawzMRVs1bkdtsgs3LrC9Gy8h1lmG03U3iaZ8aLVw5Do0z2lJKQkVavsCqnHaV4XETx3Ma4fGeKqKHkDXa0bp52hF6fC7cjvx6/r1+68LjsWSKsWCs6vY4M1BGu0LGg7Gqw0jGB82KKUYnRoOcnwznDY0Di2Km9jMnx4kn03X1+IcC0byr4s3L/Zy02NMOhBPYRHk9uB6n7bKfpz0TTjAWjBU12kY0KljnSV/RRBopFw4LgezULys8yZo9bRO6UcbnonS2OKtW9/I6rauuz8fwTPSS6V7z6fW5mArFSFlQOzQ6oy1wlNFeQlQrrALW3pCefHUUIeC27QWMts/FdDjekAIfgx8cGwXqV8WeTKUZmc0vGrFluY+JuZilC4jpcJyOFldZUZlmqB7vHzcqxwsbbcOA1UOjPxfD4OSvHjfC41bltKtJtRlRO3M87WrvaZBbqFefSNHITJQVbYWNZo/fTVpac58ZDIQB6FNGe+kwPleD0bYwPP7k0RGuW9tZsG2jq8VFIiUtV2QrRCKV5kfHtUFw0+G4ZRX0uYwFtdX2yjxFI0bLkZUh8kAZw0IMPE7bZR8ez4hSLCvc35r5jNTZ0w4XGMsJ0GJRS1M0kSKeSlflabsdWmuTWTntWiQ/3bqHXq8+7dHZKCsK5LMhR8rUgmK0WqahmYUy2hWQSksm52JVVY5DVl/Y7BvSwFSYY8Oz3F4gNA7ZMYiNKkZ78dwUM5EEr9/QRVrWJ2eZ7dHO52nrRttE/eb5TJcxLMSgGQrRTo3N4XHaihbxZD8j9V1cGh0drXkV0awx2tXM0s7FrOEqtUQPQWuz8jrtlkQi5hNPppkOJ4ouMqw12nq7WRUDo8xCGe0KmArFScvqqyyzhWjmGqwnj44AcPuOFQX3MYx2o9q+fnB0FLfDxt3Xrgasn8IDWrsXkNdIrGr34HM76LfQ054OJ+gs29O2X/Z92qfG5rhima9ousCo+6h3TtuYc59viInbYcMmzK8er3ZYiEG7t3b9cSklYzUabcjqs1tNOZXuPT7tXmeV0e71uy/Rwag3ymhXQC3CKqDdmN0Om+n5uh8cG2XbCj/rugvPdu1soKctpSb68oYrejIGtB55bWNVnK8nWAhNg9zKYrRyxnIaNEN4vH80yJYioXHQroPDJhqY017o9QohaHE5LPO0qwmPg1YjU6unPRtNEk+mqy6uNfA47UTi1r9/x2a1hXix881ImQbNv8dc1KehNRJltCug1n5G0BWfTPYiTo0G2bO+s+g+3Q30tI8Nz3IxEOGtO5ZnPP56GW2/x1HQk9my3McpC1XRpkLxzGKpFO46FqJ9+dmzPPFK/druQPMqh2aiXFGkCA00A+n3OOruaWdz2vk9KK/Lbnp1dNZo1+Bp11gfMx7UjKA5nrb1C62xjBBM4fC03+3A5bBGyvRiIMLqBgqrgDLaFZFVQ6v+n+b3OEwtREunJTORREmPrpGe9lPHtMr2W7bV22jnb/cy2LLcz8Rc3JJzicRTxJLpCgrR6mO0U2nJXz15gi8/d87y18qlf6x0EZqB3+NsnKddYNpWiwUztbPh8Wpz2rXPHh+rMXpo4HXZ66KIljHabYXPVwhBr8+dcbLMQkrJUKCxwiqgjHZFGEa7x1+e95QPs+fgzsW1SVKlFNpaXXZcDltdcsnz+cHRUXav7aTX7667p10slGVlBXklamighYUTKWlJb2kup8aChOIpzoyHLH2dha9bXHM8lzavo+7jOY2cdkueli/QepEXW3jcjJx2rRKmBh4Lrk8+xoMxhMhGDgvRY4HAylQoTjSRVuHxpcRYMIrP7cjbFlIuWnjcPC/CmFTVVsJoCyE0gRULhfTzMTitVbYb/eMep50Wl70+RnumuEawMTjklKVGu3xPG6wfunDoQgDQinTqqTrWr09GWlPGZCS/u3GedqHPdqvbYXr4d7bGQjRjuEot7ZNmRA9Bi0TUI1I0HozS3erOjNktRK/PZbqU6WJo9wJltCui1tYI0EJaQRNvlsaNtxwt9M5WV90nfT3XPwnArVdm9dDroYMejicJhBMFhwqANiXI73FkvEAzqWRYCIBH73W13mhPZx6fGbdWES6XU6NBNvX6ypqMpOW069zyZVSPF6gKtiY8rv2NvjzFb+XQ7nWSrnEA0Xgwhsthy4i1VIvXWZ/q8bHZ8hQprRgacjHPxMBGUNJoCyHWCCF+LIR4TQhxVAjxe/r2LiHEU0KIU/r3zpxjPiGE6BdCnBBC3J6zfbcQ4hX9uc+KQlMDFinjwep7tA38JofHDaPdUYbR7m6t/6SvwekwNgHrcyrbu+pwHkNF2r0MhBCsavdmBhCYSeXhcd3TTlpbgXvoQoD13Zq3W88Q+amx4prjuWg57fqGxyOJFF6nveCiQutDNt9ot7oKv2YpzNAfH9PvabXeiq1IH+RjLBgrms826PG5tRZdE9NNS8ZoA0ngv0sprwSuBz4shNgOfBzYJ6XcDOzTf0Z/7h5gB3AH8HkhhLF8fQC4H9isf91h4t9iObWooRm0eR3mhscNT7uMMGxnAyZ9XQxEWd7muSSc1VUHj99o9yoVyurxuyypMp3WPe1yZmlDfcLjM5EEp8bmePuuVdhtgjMT9fG0w/Ekg9ORso222Z+Rcig0ltPAqkK0akPjkDNTu4bC1vEyjWApvHUKj48Fo2U5Tj0+F6m0NPU+MxSI4HXayy4utYqSRltKOSylfEl/HAReA1YDdwEP6bs9BNytP74LeFhKGZNSngX6gb1CiJVAm5TyeakJYH8155glgTnhcSfxZNq0N3gl4fHuBkz6Gp6JLOiT7mpxMWlxbj1rtIvn6np8bku8fmNx1OGt0NO28Mb38kAAgL0bulnX1VI3T9v4X6ztLp3PBs3TnoslLS/Ky6XQWE4Dr0V92tUWoUFWf7w2T7s8I1gKKyIR80mlJRNz8fI8baNX28T7jDFHu9EB4opy2kKI9cC1wAvAcinlMGiGHTCSlquBgZzDBvVtq/XH87fne537hRAHhBAHxsfHKzlFy4gmUgSjSRM8bXOlTCvKabe4mI0mSaTqJ+KRr4K7Xp62ELC8rbjR7m51M2HB7N3pcByf3i9aDh6nkdO27n9z6EIAIWDXmnY29rZyuk457UrbigxZT6vmzucjFMs/ltOgxYo+7ViiJqPdbsK9xAxHBPQ+9kTK0oFEU6E4qbQsSye91wIp06FAhNWd5S08raRsoy2E8AHfBD4ipZwttmuebbLI9oUbpXxQSrlHSrmnt7e33FO0lGyVZe2FaGDe0JBAOIHLbitYQJNLlx6qrVcxmpSSoZmFU7Y6W12E4ylLvcqhmSjL/R6cJapMe/wuQvGU6V5CJcNCIOtpxyy8JocGptmyzI/f42Rjr49zk+G6eLPjZQhi5JKd9FW/vLbmaRc32hGTjVK1E74Mas1pl6PjXS5elx0pIWZhTUYl7WlZT9s8o70YhFWgTKMthHCiGeyvSSm/pW8e1UPe6N/H9O2DwJqcw/uAIX17X57tSwIz1NDAfG3lmUiCNq+zrJBNV6t27vUKkU+G4sSTaVbNC49316FX2whllaKn1ZrhApUMC4HcQjRrjLaUkkMXAly7tgOAjT2txJNpLk5HLHm9XMZmK/W0a8/VVko4nswrYWpgGCUzIyG1hseNOpZqe9onTLqnQbbq3soQ+Ziu3lZuIRpgWttXNJFiYi7OqiLdKPWinOpxAXwJeE1K+Tc5Tz0G3Kc/vg/4bs72e4QQbiHEBrSCs/16CD0ohLhe/5335hyz6KlVd9wgMwfXpEKb2UiC9jLbNYyiqHoZ7WG9gntlHk/b6vMoJaxiYAjlmJ3Xnq7Y07Y2PH5mIsRMJJE12r1aUdjpOhSjjc/FcDtsZU+z8jdgaEhJTzszntO8hUSthWg+lwObqN5omyWsArkzta002uVHbNo8Dlx2m2mqaMN6h8nqziVgtIGbgPcCtwghDutfvwh8GrhNCHEKuE3/GSnlUeAR4BjwfeDDUkrjP/lB4ItoxWmngSfM/GOsxDSj7altdTyfmUiirHw2UFc1MijcImG1p10oLJ+PbsPTNjmvHQjHM9e7HDwOawvRDFGVa9dqnZmberUWvHoUo43NRlnWVn5bUTY8XsecdjxZUMIUoEX3ws0sRpuNJqseywlgswn8NYznNEvCFLKRIiuNdiX3YCEEPT6XaUNDjIhUo4VVAEq+Y6SUPyN/Phrg1gLHfAr4VJ7tB4CdlZzgYsGQz6vkRpwPs2dqByLxsvPs3XUOjxvzrOdXj1vtaRcKy+fDyH1NhkwOj4eqDI9b5GkfujCN3+3gCt3D7mp10e511qUYbaxCfYOM0Y7V0dOOpfKO5TQw25OMJVPEk+mawuOg/R+nwrV52ksmPD4bpc3jKHsspiZlas7nemiR9GiDUkQrm/G5GF0trpKFTaUwwuNmeREzkUTZqludLU6EMN+rLMRQIILbYVuw0LHa0y63Rzv3XMxsDUmm0sxGk1WGx63ztK9Z25GZZS2EYGNva11U0caCsYqKnYwOi7p72kVy2obRNsvTrnXCl4HWxlnd59nIEfeY0PJltMtZHR6vZIFhpiraRb0bZUUZjoDVKKNdJqa1RjjtOGzCvPB4uPzwuMNuo6vFxXid9MeNEPX8sGibx4ndJhaF0fY47fjdDlN1igP6/3axFKKF40mOj8xy7ZqOS7Zv7PHVJTxeqYCHP9NhUR9PO5WWRBPpojltr9MIj5uzkKh1WIhBV2v1mgfjwRidLc6y2xKL4XVpv8NKVbRKF3+9JhvtcrpR6kHjz2CJYJbRFkKYNlM7lZYEY8mSw0JysUKTtxCFKrhtNkFni9OyiWOVCvt3+8yVVQ3of1clnrbbYV0h2uGBAGmZzWcbbOxtZSwYs7TgK5pIMRNJVFTs5HZoE+nq5Wkb3mHRnLbL3PBvsMZhIQbdvuqnWY1XaASL4alDeLzSxV+PX1vQmCFlWm43Sj1QRrtMzNAdNzBrpnYwmkCWMZYzF6tkO/MxHIgWHNjRaeHEsaFABI/TVvaErR6fuQIrGQnTCjxtIQRuh83UPm0pJY8fGeK/fP0QLS47180z2pv0/PbZCeu87Wrzpm2e+kmZhnURl0JjOSFrtEOmh8dr87R7fJpQUTWGqdJwczGM8LhV6R0pJWPBaEWLvx6fm2RamjLNbrEIq4Ay2mUhpdR0x03Q6AUtPGyGd1PJsBCDHp/bcglRgEQqzWiwcAV3Z6vLMk/7ot7uVW61siZlaqLRDlU2LMTAY+KkpNHZKPf/80F+9+uHWN3p5ZsfvHGBPr1RQW5lMVolbTq5mPUZKQfDEBcNj2c8bbPC44anXXt4PFWlYTIregjZQjSrwuPBWJJoIl3R+yjTq12jk5JOS4YC0UXjadf2jmkSZqNJ4sm0aZ62WQMRKpEwNahXeHx0NoqUFKzg7m51WTISE+D8ZJh1ZcxtzpyLz8X+c+YtIEZntfD88goXeR6nzRRPJZlKc/fnnmUqFOePfvFK3nfT+rzzh9d2t2AT1rZ9Vetp++voaYcMT7toeNzcli/jb2szITwOWvdDZwWdLVJKPTxuktG2uE/bEOipKDxuSJkGY2xZ7q/6tSfmYsRT6UVROQ7K0y4LM1sjQPugmlFkU8mEL4NunyYharaO8nxK5ZWtmjgmpeT8ZIh1OaNAS9HjczMdjpM0SZP9YiCK0y4qrsr1OO2m5LQHpiMMz0T507fv4Ldv3pjXYIOWO15j8eCQ8QpUrHKp53jOSnLa5leP1xge1w11pdGz2UiSeCptuqdtVXjcqHSv5Hx7deGkWj3to0Oaanctht9MlNEuA7ONtuZF1H5DCoSr87QB00QHCmH0aBcKKXW3Vp+LK8bEXJxQPJWZGV0OPT4XUmJauH4oEGFluzfTXlUuHoc54w1P6xGMLStK32Q29lg7OGQsGMMmshoB5eL3OOpWiBYqI6ftdtiwCfML0XxF2szKoctXnaLf+FzlRrAYTrvAbhOWOQOV6tcD9Or7Gl56tRwaCGATcNXq9pp+j1koo10GxkrNrFCSlq8zLzxeSU6716Q8TykMNbRihWhpWdtYwXxcmNK8xko9bajcWynEUCBSVSjN47QRNWHggmGEN/WUnl+9qdfH2YmQ6Ysng7HZGN0+N/YKFzD1zGkb3nMxT1sIQYuJ4zmDUW2qWKEoSLkYi6HJCj/PlerBl0IIQYvTTiRujThQNefb7nXi9zgy96JqOTwQYMtyf9E+/nqijHYZnB6bwyYKG6BKafM6CcdTNY/INAxepS1fYP6AjPkMB6K0e50F3+jdVXoIpTg3EQZgXQWedrfJ16Rc3fP5aOFxEzzt8Tl6fO6y0iYbe33Ekumab2yFGJ+rLm9qVodFOWRz2sWVtrwuO5GEeYVotYbGIUcwqcIFp9mOCIDHxOszn0r16w36OlsYnA5X/bpSSl4eCHDNPI2DRqKMdhkcGQyweZl5K602k7SVZyMJ3A5b2bJ+kDMgw+IK8lKGy6isNntM6PnJEDahfVjLpcdn3jVJpNKMzEarGuHncdpNafk6PR7KVIaXwtiv36IQeaVtOgZ+j5NIovaFbTmEy6geN54309OutUcbNMGkDq+zYqGirOdqXkW012m3rE+7Uv16g75OL4M1TLI7qw/aUUZ7CSGl5OXBGa7uMy+f4TdpaEglw0IMui0aRTmfoZloUe3vrioLaEpxfirMqg5vRSpPZnrao7NR0rK6wQJa9XhtRkpKSf/YHJuWlQ6NA2xb2QbAa8OzNb1uIcZmq2srMrzQuTrktUN6HrbUotzrNNtom+MEdFfRsjg0o2kZ1DKwZD7GzHErqFQNzcAw2tXOQX95MADANfp0vMWAMtolGJyOMBWKs8vElZZZ2sqBCsc/ArgcNtq9TuuNdglP2zDaZnva5ybDrK8gnw3ZMX5m6I9XqsaWi8dpr1nGdCoUZyaSyAinlKLd66Sv08uxIfONdiotmZir7mZbT/3xSDyFTWRV6Qqhedpmhsdr97RBK+qs9L17ZHCGHavaK/Zci+ExcVEzn7Eq29P6OluYiyWrrp05fCFAi8vO5mWLo3IclNEuibHS2tXXYdrvNFa3tVaQV+Npg5ZPttJoh/QPycoiIWKrxoRq7V6VKRcJIUy7JplpQFXM3TWjevy03r5VbngcYPvKNkuM9lQoTlpW3u4FOfrjdShGC8VStLocJQ2YmYVos6Z62q6KPkfxZJpXL84s0KKvlRaXOTUZ+RibrS7N0qd/DqsNkR8eCHDV6vaKCymtRBntErw8EMDlsLG1jPaZcjG8iEaEx8GQ7bQup220exWroPY47bS47KYa7ZlwgkA4UbHRBkMprnajbRR0raqiaNGM8HimcrxMTxtg+6o2zk6GMgVZZpHpra1ClKieRjsQjpdVzOl1mZOznYslOT8ZYlNPZRGhQnS3VvbePT4ySyyZNj3k6zVR0S+XaCLFbDTJsrbqwuMAA1OVF6NFEymODc8uqtA4KKNdkpcHZ9i+ss2USTgGZt2QZiKJiirHDcycfpMPI0Rcqtq+q7UyD6EU56to9zLQPO3az+ViIEJXq6vobOZCmFE9fnpsDrfDVlHL2Y5V7UgJx0eCNb32fDISplV42oZSWD3C4+enwqzpKn29zCpEO3xBG+Cye31Xzb8LtPfudDhRtjjQoQsBYOEAmVrxmFiol0stOhlGQWo1nvax4VkSKWl6RKJWlNEuQiotefXijOmVg2bl62YiCTq8lelbg1YtbWWfdnY0ZvGVsdlG+9yktpquNKcN5nnatUwDcjvtxJLpqotmQPO0N/b6KhJ22b5KK0Y7ZnIx2rghPVlNTrueRrvMOgizwuMHzk8hBFxrkgeXmU9fZn3I4YEAy/zuooWi1dDitBO1wGiP1WC0jV7tatq+DuuLm2vWmLu4qRVltIvQPzZHOJ4ytXIcwOdyIERt4fFkKs1cLFl1eDwYTRKzYHYzaJXjQsDyEuEss432hUnN015bge64QY8+4rAWgwm60a6yn9/j1D6OsRoEVs5MlN/uZbCq3UO712l6XttYGNZSPW61wMpcLMnEXIy1ZaRUWlx2UwaGHDw/zdbl/pp1xw2M7odyP0uHLkxz7doOU4vQQEsfhC0Ij2ekcKvsKdd6tSv3tF8eDLCizcMKkxc3taKMdhFeHggAmFo5Dto8ab+7toEIxrHt3sqLWXr85iqAzWeozIHxXS3me9rL29xVhaZ7fC7iqXRN/xMpJRenI1UVoYFWiAbV6zdHEykGpsJsrCCfDVoh3o5VbRwbmqnqdQsxNhvF73FUpCNg4DNSSBYLrFzQozPrusrxtDWjVMvCLpWWHLoQYPc687y37graJ6dCcc5Nhi3xHq3q0652UpxBtb3ahwcC7FqzOKRLcylptIUQ/ySEGBNCvJqzrUsI8ZQQ4pT+vTPnuU8IIfqFECeEELfnbN8thHhFf+6zwuxlngW8PBjA73awoYpwayn8NQ4NqWZYiIHVqmjDM5GileMGpue0KxwUkktWyrT6azIbTRKKp6qeBuTJDF2oztM+PxkmLSurHDfYvrKN4yNB04amQPVtOgBOuw2v0265p52VvS3taXtddqSsLRJyYiTIXCzJnvUmGm1dHKicz/PhgWnAvNB8Ll6Xlt4xWxL3zHiIVpc9I4JUKZrRDle02JoKxTlv0eKmVsrxtL8C3DFv28eBfVLKzcA+/WeEENuBe4Ad+jGfF0IYy+wHgPuBzfrX/N+56Hh5MMDVa9orHvxQDm1eZ01eXUDPX1Wb0wbrjLY2e7a04epsdRFJpExbnVc6kjOX7I2v+kXExWkjl19beLxaT7uaynGD7avaiCXTnJ0wb+LXeJWCGAZtXuuHhhh1EOUY7RYTZkYfPD8FwJ515hShQVYwqZwF8OEL2vALs1N+kJ30ZXYF+YmRIFtW+KsO5/d1thCKpzIDlsrBiLIuJiU0g5JGW0r5DDA1b/NdwEP644eAu3O2PyyljEkpzwL9wF4hxEqgTUr5vNSWO1/NOWZREk2kOD4c5GoT+7Nzaatx0lc1uuMGVk76SqTSZQ/MqLSAphjheJKxYIz1VbbRmOFpZwvwavO0q73pGdO9NlbjaVtQjDYWjFVVOW7g9zgJxqz1tM9PhulqdZUldGLM1K6lNe7A+WmW+d2ZViQzaPc6sdtEWeHxQwMBtq1oKzo7vFpaLJqpfXI0yNYaxmJW06ttTPayYnFTK9XmtJdLKYcB9O/L9O2rgYGc/Qb1bav1x/O3L1qODc+STEtTRVVyafOaFB6vwWhbUUF+ZHBG6wEtY4XaaaiimRAiP1+Bx5SPSkKMhRgqMY60FGZ42qs7vFXdkDf1+nA5bJnZwbUipWQsGK2qR9ugHuM5L0yVL8bjNcEoHTg3zZ71naYWgdlsgs4WV0kp03RacvhCwLK+48yi08S89ngwxmQoXtMs66zRLr+CfP/ZyUU12SsXswvR8r0TZZHt+X+JEPcLIQ4IIQ6Mj4+bdnKVcCRThGbNSqvWG9JsDUbb67LT6rJbUoj2XP8EANdv7C65r+FpmxGmP19BQVE+ulpcVU1LyuViIILLYaOnwtnRBtlCtOpypqfHQ1V52aDlkLcu95tWQR6MJYkm0jV72rUKEJXi3ET5KRXDk6w2PD4yE+ViIMJuE0PjBj0+V8nP8+nxOYKxpGV9x2YsauZzclTTDthWg7iV0as9UKbRvhiI8MLZKe7YuaLq17SSao32qB7yRv8+pm8fBNbk7NcHDOnb+/Jsz4uU8kEp5R4p5Z7e3t4qT7E2Xh6cYZnfzYoqVHjKoc3jrCk8buRnqjHaoFWQW5HTfu70JNtXtmVkSouxQQ9lmyHqcd5o96rS03bYbXS21CZlOhTQhqRUWwPhNgrRqmjFk1Jyenyuqny2wfaVbRwbnq257Q2yghg15bQt9rTjyTTDMxHWllm86M0Y7erO6eB5rQjMzMpxg26fq+SY20O6I2K2qIpBJjxuoqd9Qr83bKnBaGd7tcsLj3/7pUGkhHdc11d65wZQrdF+DLhPf3wf8N2c7fcIIdxCiA1oBWf79RB6UAhxvV41fm/OMYuSlwcCXN1nfi+jQZvXyVwsWXWl5UwkQYvLXrVSW48FqmjRRIqDF6a5cVNpLxu0/tI1Xd5M0UctnNNzk9UuYqA8b6UY1c7RNsj0aVfhqYzMRgnHU2VP98rHjtVtTIXijM7W/r4Yywir1OhpW2i0B6e1avv1ZS70WvW0Q7VG6cD5KTxOGzv0+gEz6SpDyvTQhQBtHgcbTZJPnY/HhEK9+ZwcDdLd6sqk9Kql3F5tKSWPHhzk+o1drKmyqNVqymn5+gbwPLBVCDEohPgA8GngNiHEKeA2/WeklEeBR4BjwPeBD0spjf/gB4EvohWnnQaeMPlvMY2To0HOTIR4/Qbzw1gGbR4HUmphxGqoVnfcoMeCoSEHz08TT6a56Yqeso+5Zk2nKUb7wlSoKlGVXLpba1vIXJyuzWh7a2j5Oj1W+aCQ+WzXx3QeNaFfO6M7XoPR1jxt68LjldZB1BoeP3h+ml19HSX1C6qhu7X0gvPQhWl2remwpBsGct+/Jnrao8Ga8tkGa/S2r1IcPD/NucnwovWyobzq8XdLKVdKKZ1Syj4p5ZeklJNSylullJv171M5+39KSrlJSrlVSvlEzvYDUsqd+nO/K82IwVnEP/3sLB6njV/bbd0/LivTWN1NqXaj7TZFazuX505PYLcJXlfBYmdXXztDM1HGZqM1vfa5iXDZHlMhevzukiHGQiRSaUaD5bW6FcJTw03PaPe6oobwuDFb24y8thnhcb/HQSyZtky5L5NSKbMOwltD+DccT3J0aNbU/uxcenwugrHCKoehWJKTo0HLQuOQra43K6ctpeTkSNCUYU2Gp13K7HzzpUFaXHZ+8aqVNb+mVShFtHlMzsX41qGL/Op1fZnqZito89am+BSocliIQbfPzXQ4bqqYxrP9k+zqa8dXQcWlUWV+uAZvO5ZMMTQTqVpYxaC71cVEsDpPe2QmipTQ1yCjfXI0iN/tqMmz9bkdrO9uMaXtazwYw+WwZd7n1WD1TO3zU2FaKhDtMIxSNTntwwMBUmlpan92Ll0lerWfPz1JWsJek4aU5MNrcnj8YiBCKJ4yxdPu6/QSjqeYLtKrHU2kePzlYd62c+WirBo3UEZ7Hv/y8wvEk2nef9MGS1/H8LSrLUabjSToqMFo9/pcSGnePOvZaIIjg4GKQuMAO/VZtcbc8mrQVtDVt3sZ9PrdetVz5TedWnu0Iaflq0LFLSklz5waZ++GrpprMHasbuflgUDNxWiGGlot55PVH7fIaE+GWdfdWvY5ZsLjVbw/nj89iU3AdRZ5ukbLYqEQ+b7jY7S67Oy1MOVndvW4UYS2dUX10SODctq+njw6QjCWtDTCagbKaOcQS6b455+f581be7mihoKecqh1prYZ4XGorcUpl/1npkhLuKHMIjQDj9POthV+Xh6oPo968JxWlVvritzwuKoJkdfaow3Va4+fHp9jYCrCW7YtK71zCW7Y2M3QTJQzNSijheNJnj4xVnPBld9dWwqpFOcnQxUp6LkdNoSAcKyy/4+UksePDHP9xu6qZIfLodh7V0rJj46PcvOWXlNHDM8nmz4wZ5F1Qm/32myKp116ROejBwdZ3eG1tJbJDJTRzuGxw0NMzMX4wBus9bKhdi+iZqPtN1d//LnTk7gdtqo8iV1rOnh5MFB1Jf2/HxlibVdLzUbCyL8OTFU+xq9WCVPQRDJaXfaK/yf7XtM6Lm8xwWjfvFlrsfzpyer1ER7eP8B0OMH9N2+s6Vz8Fg4NSaUlA1ORiqIzQghanJXPjD46NMvZiRC/vGtVpadZNkZ4PF8F+dGhWUZnY6a8P4phhMdDFS5qCnFyJMiqdo8p09BWl/C0h2ciPNs/wTt291lWqGcWymjrSCn50s/OsnW5nzdUGOKthlrC4/FkmnA8ZZKnbZbRnmDP+s6qJjpd09dBMJqsyrubnIvx3OlJfunqlTWHhves78Rlt/HDY6MVH3sxEKW71VXV35/L3g1d/OzUREXH7Ds+xpUr22paMBis7W5hfXcLz1R4DgbxZJr/99Mz7N3QVbOISDanbb6nPTIbJZ5KV1wH4XU5iCQqW0Q8fmQYh01wxw7rxDqKhcd/dHwMIeDNW6012nabYHWHl369KLJWTozO1dSfnUu710lbgV5tKSX/8ztHcdhsvHORh8ZBGe0Mz5+e5PhIkPe/Yb1lvdm51OJFGBKmHTWE2swcGjIxF+P4SJAbN1W32DFkFatp/fr+0RFSacmdV9de7en3OHnj5h6eeHWk4pxurT3aBrdsW8a5yTBnyrzxBcJxDp6f5lYTvaibt/Ty/OnJqqq2v3PoIsMzUT705k01n4eVOW2jcrzSOogWV2WethYaH+KmK3osLWz1ux247La84fF9x8fY1ddRU5FiuVyzpsOUFs5kKs3psbmaNMfn09fZkjeK9rUXLvDD10b5w7dtW7S92bkoo63z4E/P0N3q4q5r6iOJ7rDbaHXZq/K0axkWYuBzO3A7bKbktH9+ZhKgbFGV+Wzq9dHqsldVjPb4y8Ns7G3N9BjXytuuWsnFQIQjg5Xl2MsdklIKwxv60fGxEntq/OTkOKm05JYrTTTam3uJJFKZWoFySaUl//iT0+xY1cabttSuZOivsVizGMYc7Up7+ztbXQzPlN+e+PLgDIPTEVMWlcUQQtDV6loQHh8PxjgyGDB1UVeMa9Z0MDgdqdkZODcZJp5Km1I5bpBvrnb/WJD/8x/HeOPmHt5343rTXstKlNEGXr04w9Mnxnn/GzbUHN6sBL/HWVXor5ZhIQZCCK1Xu8oWp1x+dHwMv9vBVaur02m32wRX9bVXvEIfm43y87OT3Hn1KtOiI7dduRyHTfC9V4fLPmZiLsa5yRAbahA2MVjT1cLmZT5+fKI8o/2j42N0t7pMHWxz/aZuHDZRcYj8yaMjnJkI8cE3bzLl/2G0DlrhaZ+bDOO0i4qjI6/f0MXhC4Gye7Uff3kIl93GWy0MjRvkkzJ9+sQYUmLqoq4Yu/QWzlq97WzluLme9vmpMI+8OMBMJEEsmeK/fuMwLS4Hf/3OXYs+l22gjDbw+af78XscvPeGdXV93dWd3qp0t2ci2gezFqMNuipajS1fc7EkT7wywi9dvRJHDUpP16zp5NjwbEWV0997ZRgp4ZdN9GLaW5zcdEUPT7xSfoj8kQMDJFKSd1xnTpTmlm3L2H92irl5annh+KXiGclUmqdPjPPmrcuwm3jD8bkd7F7XyTMVFKNJKfn80/1s6GnlbTvN+X/YbQK/u7YRtoW4MBViTWdLxdftxk3dxFNpXjw3f1rxQtJpyX+8MszNW3pq/qyWQ7dvoZTpj46PsaLNY1okqhQ7V7dht4madBdAqxy3CUzt4vmlq1ewst3Dx755hNf9nx/y9r9/lmPDs/zfd1zNMotmTFhB0xvt/rEgT7w6wn03rDelSrES3rK1lyODMxnJx3LJ5rRry5GZ4Wl/78gwkUSKd+6prYDjmjXtJFKS1yoQ9nj8yDBbl/tNaQnJ5RevWsGFqXBZYypTacnXX7jADRu7uWKZOefxlm3LSKQkPzuVNZqJVJpf+dxz3PY3z2TysS9dCDATSVhSFXzzll6ODc9mlM1K8fTJcV69OMvvvGmjqQsIq8Zznp8MVzVcZu+GLpx2wbOnS0chDl6YZngmyp1XW1c1nkt366WedjyZ5pmT47xl27K61OmAJkCzZbm/ZqN9ciTI+u5WUyOfu9d18fRH38x3PnwT77l+HXOxJL/1hg3ctn25aa9RD5reaD/w9Bk8Djvvu2l93V/7lm3am+Xp45W118zUOOHLwIyhIf92cICNva01i0ZUGlYbCkQ4cH7aklzhbdtXYLcJnigjRP6Tk2MMTkd4z/XmRWl2r+vE73Fcktd+6LlznBgNMjEX4x0PPM/RoRn2HR/FYRO8cYv53Q5GTvpn/aXfm+m05K+ePMGaLi+/cq251bfVppCKIaXk/GSY9VUo6LW4HFy7ppPn+idL7vv4y0O4HTZ+oU5GYb7++P6zU4Tiqbrlsw2MYrRqWzhBU/gzM59tIITgmjUdfPKXt/Psx2/hj+/cbvprWE1TG+2BqTDfOXyRd+9dS3eNU2Sq4cqVfla2e9h3vPwWo0A4zsMvDuB3O2jz1Ca119fpZXwulinKqZSzEyFePDfNr+3uq3klv6LNwzK/mxfPTXNiJMi+10Z5eP+FgouK772iGdQ7Leh97Wp1ccPGbr5XRoj8X35+gV6/m7fuMO/G7LTbuHlLLz8+MU46LZmYi/F3+05x85ZeHvvdN+CyC+75ws/5zqGL7N3QZUmEaPvKNrpbXTxzsrRH+b1Xhzk6NMvv/8IW08U7/B5HJrJkFlOhOHOxZNUDZm68optXh2YIhAunlhKpNN97dYRbti2rSNa3Frp9biKJFOF4kuMjs3zlubO4HbaKVQpr5Zo17cxGk5ydrE6g5+xEiHOTIdPavS43mtpoP/jMGWwCfvtm68VU8iGE4JZty/jpqYmy2mtmIgne+6X9nBkP8bnfuK6mHDLAO/eswWmz8cBP+qs6/tGDA9iEOXNnjRXwf7wyzO2feYYPPHSAj3/rFd7zxRcW5HaD0QTf2H+BnavbMjO5zeZtV63g7EQoo8qUj4GpMD8+Mca7X7fG9MlNb9m6jPFgjGPDs/z1D04Siaf45J1XcsUyH9/80I0sb/dYKphhswnesLmHn54aL+oxJVNp/uYHJ9my3GdJ58WOVW28dCHAtElyuwD/emAg87ur4aYrepAy2zWRj4dfHGA8GONde9ZU9RrV0K23lL3pL5/mjs/8lB++NsZvvH5dRqmsXlyzRou6VVOMFk2k+PDXXqLN6+Se19Xv2i0lmtJoT8zF+NyP+/nXAwP82u4+VrbX3qpTLbdeuYxwPMULZ4oXtszFkvzml/dzfGSWf3zvddxsQkvNinYP/+l1a3j04CAXA+UNiDdIpSXfPHiRN23pZblJRRwfvX0rH3/bNj777mv59odu5B/fs5tTY3N85OHDGcMRTaT4rYcOcH4yzEffutWU183HW7evwCa0nH0hvr7/AgK4Z+9a01//zVt7EQI+9+N+Hn7xAvfesD6TM1/Z7uXR37mBP7h9K//JwhvbzZt7mZiLFx0g8s2XBjkzEeK/v3Wrqblsg3v2riWeTPOtQxdN+X0nR4N85qlT/OJVK3j9xupaFHf1ddDisvNsgRB5KJbk7354ir3ru3jz1to/p+Vy9Zp2VrV7uHZNB5/+1at44X/cyid/uf7h3yuWaS2c1eS1//x7r3FseJa/+rVdpugeXI40ldF+8dwU//Ubh7jhz/fxl0+eYM+6Tj7yC1saek43burB47QV7cudCsV535f3c2Rwhr9/93WZXLgZ/I4ugvGPT5+u6Lif9U8wMhvlnSZ6EluW+/mdN23i7btWce3aTu7YuYJP3rmdH742yl/+4ASJVJoPfe0l9p+b4q/ftctShadev5sbN/XwwE9O8z++/coC+cNYMsUjLw5w65XLLbm59PjcXN3XwROvjtDZ4uL3fmHzJc93tLj48FuuyPQyW8Ebt/Rgtwn+6DuvZgai5BJNpPjMD09xzZoO3mpR3vbKlW1cu7aDb+y/UPMQk2QqzUf/7WV8Hgf/+66dVf8el8PG3g1dBYvRvvSzs0zMxfjDt22rWwEYwLYVbTz3iVt58N493LN3rWmL6UqptoXz+68O89Dz5/mtN2yoWx3AUqSpjPaff+81fnxijPdcv44f/rc38fXfvr5hb2wDj9POTZt62Hd8NO9N6cC5KX7x737KywMz/N0913DHTnP7PVd3eHnHdX3864EBRiuYaf1vBwboaHFyq8X9n/fesI5ff/1aHnj6NL/2wHP86PgYf3bXzrqI4PzNf9rFPa9by6MHBnnLXz3NHz56hL996iR/8thR/vM/H2QyFDe1AG0+t+iLko++dWtdWobms8zv4R/efS39o0Hu/Puf8Wx/1kil0prs7/BMlI/dvtVS4/TuvWvpH5vjwPnKxF7m84VnznBkcIY/u2tnRsa3Wt5wRQ9nxkOMzBNamZyL8YWfnOaOHSvYvc662dWLnV1rOipq4RyYCvMHjx5h15oOPnbHNovPbmmzeIeGWsDfvOsalrW5M3NxFwu3XrmcfcfH6B+by7QvpdOSB396hr988gR9nV6+9aEb2VmleEkpPvTmK/i3g4N84SdnygqnDQUi/ODYKL++dy1uh7X5MiEEf/r2HZwdD/H8mUk+dsdWSw1lLsv8Hv7s7p188M2beODp0/zriwPEU2naPA46Wly8becK3mhhkc99N66jy+eyNAReirddtZItK/z8zj8f5L1feoFfva6PkZkohwcCzMWSvHFzDzdaXOh059Ur+bN/P8bXX7jA66qcB31iJMhnfniSX7pqJb9kQseBIdlrDJkw+Psf9RNNpvmDO6xL3SwFrl3TkWnhvLZEZ8nITJT7vrwfgH9497WWTiK7HFhc1sti1ltUtFQrRjHRvuNjbF7uZ2Qmyse+eYRnTo7zS1et5M/fcZWlPeRru1u4+5rVfO2F83zwzZuKahQHowne/5UXcdlt3Fcn2T+n3cb/u28PrwzOVDz60wxWdXj5s7t38sd3XonDZrMkd5uPjhYX763TAqUYm3p9fOfDN/FH336F7xy6yJblfu6+dhXXrunkdpMjP/locTn4letW8/CLA/yvX95esT7BVCjO7379Jdo8Tv73XTtMOadtK/x0tbp49nTWaF+YDPO1F87zrj1r2NRr7WjfxY7Rwnl4IFDUaA9Mhfn1L/6c6VCCL7/vdUtC+7vRNJXRXqysaPewY1Ub+14bZUWbh09+91USKcmf3b2T97x+bV3yYh9+yya+fWiQT3zrCJ9997V5oxFGTrl/bI4vv+91llVu58PndjTEYOdidVRhMdPqdvCZe67lr991Td0WLbm8e+9avvr8eb710kXeX8Ho3EA4zm988QUuTIX58vteZ1prp80muGFTN8/2T3B2IsSjBwd49OAgdpvgI/PqD5qRle1elre5i+a1T4/P8Z4vvkA4nuJrv/X6jKFXFEcZ7UXCrduW8dkf9fPiuWl2r+vkr965q65GcWOvj0/euZ0/ffwY7/rC83zx3texoj2b75dS8kfffoWfnprgL95xNW/cXL+qWMXioREGGy4tSHvfTeVN4psJJ3jPl17g9PgcX7x3T9VT6Apx06Ye/uPIMG/5q6exCU2Q5v6bNzW8TmaxsKuvg8MDARKpNDYhEMD5qTAvnZ/m0MA033tlBJuAh++/nivrJLN6OVB3oy2EuAP4O8AOfFFK+el6n8Ni5O3XrOKbL13kvTes47ffaK4UZLn85k0bWNPVwn/9xiHu+tzP+Idfvw67TXB8OMhzpyd4/Mgw/+WWK3iX6p9UNIB3713Lxx49wscePcKtVy7jhk1ZTe94Ms10OM5UKM50KM5kKM4Xf3qGkyNzfOG9u01pkZzP7TuWs++1UXav7+RXr+27ZJGrgGvXdvKDY6Ns/qMnFjxn6Nv/zzu3m6ov3gyIWtsoKnoxIezASeA2YBB4EXi3lPJYoWP27NkjDxw4UKczVAC8NjzLbz104JLebb/bwTt29/G/fnl7XdtYFAqDaCLFH37zCD88NkoonsImtDDsbCRBMLZQn9zlsPG5X79uyWlLXy4EwnEeOTBAPJkmLSEtJcvbPFy3tpMrlvkaFrVZCgghDkop9+R9rs5G+wbgT6SUt+s/fwJASvnnhY5RRrsxTMzFePzlIVZ1eLlyZRt9nV5lrBWLgkQqzaELAX7WP8GFyRAdLS66W110trrobHHR1eqi2+diRbun7kOAFAozKGa06x0eXw0M5Pw8CLx+/k5CiPuB+wHWrjVfbUpRmh6fm9+8qTHyrgpFMZx2Tdxk74bq2r8UiqVMvRvi8rlqC1x9KeWDUso9Uso9vb2q4EmhUCgUCqi/0R4EcquY+oChOp+DQqFQKBRLknob7ReBzUKIDUIIF3AP8Fidz0GhUCgUiiVJXXPaUsqkEOJ3gSfRWr7+SUp5tJ7noFAoFArFUqXufdpSyu8B36v36yoUCoVCsdRRyuwKhUKhUCwRlNFWKBQKhWKJoIy2QqFQKBRLhLoqolWDEGIcOF/Dr+gBJkw6ncsRdX2Ko65PcdT1KY66PsVR1yc/66SUeUVKFr3RrhUhxIFCcnAKdX1Koa5PcdT1KY66PsVR16dyVHhcoVAoFIolgjLaCoVCoVAsEZrBaD/Y6BNY5KjrUxx1fYqjrk9x1PUpjro+FXLZ57QVCoVCobhcaAZPW6FQKBSKywJltBUKhUKhWCIsSaMthPgnIcSYEOLVnG3XCCF+LoQ4LIQ4IITYq2+/TQhxUAjxiv79lpxjduvb+4UQnxVC5Jv3veSo5PrkPL9WCDEnhPhozrbL7vpUem2EEFcLIZ4XQhzVr4VH337ZXRuo+LPlFEI8pF+H14QQn8g5ppmuzy79PfKKEOLfhRBtOc99Qr8GJ4QQt+dsb/rr04z3ZlOQUi65L+Bm4Drg1ZxtPwDepj/+ReBp/fG1wCr98U7gYs4x+4EbAAE8YRy/1L8quT45z38T+Dfgo5fz9anwveMAjgC79J+7Afvlem2quD6/DjysP24BzgHrm/D6vAi8SX/8fuDP9MfbgZcBN7ABON2k759C16fp7s1mfC1JT1tK+QwwNX8zYKxw24Ehfd9DUsohfftRwCOEcAshVgJtUsrnpfYu+Spwt+UnXwcquT4AQoi7gTNo18fYdllenwqvzVuBI1LKl/VjJ6WUqcv12kDF10cCrUIIB+AF4sBsE16frcAz+uOngHfoj+9CW9TEpJRngX5gr7o+2vVpxnuzGdR9NKeFfAR4UgjxV2hh/xvz7PMO4JCUMiaEWA0M5jw3CKy2/Cwbx0fIc32EEK3AHwK3AR/N2b+Zrs9HyP/e2QJIIcSTQC/aDfgvaK5rA4Wvz6NohmkYzdP+fSnllBBiD811fV4F3g58F3gnsEbfvhr4ec5+xnVIoK7PfJr53lwRS9LTLsAH0W4aa4DfB76U+6QQYgfwf4H/bGzK8zsu5/63QtfnT4G/lVLOzdu/ma5PoWvjAN4A/Ib+/VeEELfSXNcGCl+fvUAKWIUW/v3vQoiNNN/1eT/wYSHEQcCPFnGAwtdBXZ8c1L25Mi4no30f8C398b+h3VAAEEL0Ad8G7pVSntY3DwJ9Ocf3kRMyvgwpdH1eD/yFEOIcmkf1P4QQv0tzXZ9C12YQ+ImUckJKGQa+h5ava6ZrA4Wvz68D35dSJqSUY8CzgOFlN831kVIel1K+VUq5G/gGWu4atOuQ61Ua10FdHx11b66cy8loDwFv0h/fApwCEEJ0AP8BfEJK+ayxs5RyGAgKIa7XKxPvRQvfXK7kvT5SyjdKKddLKdcDnwH+PynlPzTZ9cl7bYAngauFEC163vZNwLEmuzZQ+PpcAG4RGq3A9cDxZrs+Qohl+ncb8MfAP+pPPQbco+dpNwCbgf3q+mjXR92bq6TRlXDVfKGt1obJ5oY+gBa+PIhWrfkCsFvf94+BEHA452uZ/twetHzLaeAf0BXilvpXJddn3nF/wqXV45fd9an02gDvQSuSeRX4i8v52lR6fQAfmud9FDgG/EGTXp/fA07qX5/O/VuBP9KvwQlyKqDV9WnOe7MZX0rGVKFQKBSKJcLlFB5XKBQKheKyRhlthUKhUCiWCMpoKxQKhUKxRFBGW6FQKBSKJYIy2gqFQqFQLBGU0VYoFAqFYomgjLZCoVAoFEuE/x/to34sWUCRNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "URL = \"https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/lynx.csv\"\n",
    "lynx = pd.read_csv(URL, index_col=0)\n",
    "data = lynx[\"value\"].values\n",
    "print(\"Length of time series:\", data.shape[0])\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(lynx[\"time\"], data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time series has a length of 114 (a data point for each year), and by looking at the plot, we can observe [seasonality](https://en.wikipedia.org/wiki/Seasonality) in this dataset, which is the recurrence of similar patterns at specific time periods. e.g. in this dataset, we observe a cyclical pattern every 10 years, but there is also a less obvious but clear spike in the number of trappings every 40 years. Let us see if we can model this effect in NumPyro.\n",
    "\n",
    "In this tutorial, we will use the first 80 values for training and the last 34 values for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train, y_test = jnp.array(data[:80], dtype=jnp.float32), data[80:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model we are going to use is called **Seasonal, Global Trend**, which when tested on 3003 time series of the [M-3 competition](https://forecasters.org/resources/time-series-data/m3-competition/), has been known to outperform other models originally participating in the competition:\n",
    "\n",
    "\\begin{align}\n",
    "\\text{exp-val}_{t} &= \\text{level}_{t-1} + \\text{coef-trend} \\times \\text{level}_{t-1}^{\\text{pow-trend}} + \\text{s}_t \\times \\text{level}_{t-1}^{\\text{pow-season}}, \\\\\n",
    "\\sigma_{t} &= \\sigma \\times \\text{exp-val}_{t}^{\\text{powx}} + \\text{offset}, \\\\\n",
    "y_{t} &\\sim \\text{StudentT}(\\nu, \\text{exp-val}_{t}, \\sigma_{t})\n",
    "\\end{align}\n",
    "\n",
    ", where `level` and `s` follows the following recursion rules:\n",
    "\n",
    "\\begin{align}\n",
    "\\text{level-p} &=\n",
    "\\begin{cases}\n",
    "y_t - \\text{s}_t \\times \\text{level}_{t-1}^{\\text{pow-season}} & \\text{if } t \\le \\text{seasonality}, \\\\ \n",
    "\\text{Average} \\left[y(t - \\text{seasonality} + 1), \\ldots, y(t)\\right] & \\text{otherwise},\n",
    "\\end{cases} \\\\\n",
    "\\text{level}_{t} &= \\text{level-sm} \\times \\text{level-p} + (1 - \\text{level-sm}) \\times \\text{level}_{t-1}, \\\\\n",
    "\\text{s}_{t + \\text{seasonality}} &= \\text{s-sm} \\times \\frac{y_{t} - \\text{level}_{t}}{\\text{level}_{t-1}^{\\text{pow-trend}}}\n",
    "+ (1 - \\text{s-sm}) \\times \\text{s}_{t}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more detailed explanation for SGT model can be found in [this vignette](https://cran.r-project.org/web/packages/Rlgt/vignettes/GT_models.html) from the authors of the Rlgt package. Here we summarize the core ideas of this model:\n",
    "\n",
    "+ [Student's t-distribution](https://en.wikipedia.org/wiki/Student%27s_t-distribution), which has heavier tails than normal distribution, is used for the likelihood.\n",
    "+ The expected value `exp_val` consists of a trending component and a seasonal component:\n",
    "  - The trend is governed by the map $x \\mapsto x + ax^b$, where $x$ is `level`, $a$ is `coef_trend`, and $b$ is `pow_trend`. Note that when $b \\sim 0$, the trend is linear with $a$ is the slope, and when $b \\sim 1$, the trend is exponential with $a$ is the rate. So that function can cover a large family of trend.\n",
    "  - When time changes, `level` and `s` are updated to new values. Coefficients `level_sm` and `s_sm` are used to make the transition smoothly.\n",
    "+ When `powx` is near $0$, the error $\\sigma_t$ will be nearly constant while when `powx` is near $1$, the error will be propotional to the expected value.\n",
    "+ There are several varieties of SGT. In this tutorial, we use generalized seasonality and seasonal average method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are ready to specify the model using *NumPyro* primitives. In NumPyro, we use the primitive `sample(name, prior)` to declare a latent random variable with a corresponding `prior`. These primitives can have custom interpretations depending on the effect handlers that are used by NumPyro inference algorithms in the backend. e.g. we can condition on specific values using the `condition` handler, or record values at these sample sites in the execution trace using the `trace` handler. Note that these details are not important for specifying the model, or running inference, but curious readers are encouraged to read the [tutorial on effect handlers](http://pyro.ai/examples/effect_handlers.html) in Pyro."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sgt(y, seasonality, future=0):\n",
    "    # heuristically, standard derivation of Cauchy prior depends on\n",
    "    # the max value of data\n",
    "    cauchy_sd = jnp.max(y) / 150\n",
    "\n",
    "    # NB: priors' parameters are taken from\n",
    "    # https://github.com/cbergmeir/Rlgt/blob/master/Rlgt/R/rlgtcontrol.R\n",
    "    nu = numpyro.sample(\"nu\", dist.Uniform(2, 20))\n",
    "    powx = numpyro.sample(\"powx\", dist.Uniform(0, 1))\n",
    "    sigma = numpyro.sample(\"sigma\", dist.HalfCauchy(cauchy_sd))\n",
    "    offset_sigma = numpyro.sample(\n",
    "        \"offset_sigma\", dist.TruncatedCauchy(low=1e-10, loc=1e-10, scale=cauchy_sd)\n",
    "    )\n",
    "\n",
    "    coef_trend = numpyro.sample(\"coef_trend\", dist.Cauchy(0, cauchy_sd))\n",
    "    pow_trend_beta = numpyro.sample(\"pow_trend_beta\", dist.Beta(1, 1))\n",
    "    # pow_trend takes values from -0.5 to 1\n",
    "    pow_trend = 1.5 * pow_trend_beta - 0.5\n",
    "    pow_season = numpyro.sample(\"pow_season\", dist.Beta(1, 1))\n",
    "\n",
    "    level_sm = numpyro.sample(\"level_sm\", dist.Beta(1, 2))\n",
    "    s_sm = numpyro.sample(\"s_sm\", dist.Uniform(0, 1))\n",
    "    init_s = numpyro.sample(\"init_s\", dist.Cauchy(0, y[:seasonality] * 0.3))\n",
    "\n",
    "    def transition_fn(carry, t):\n",
    "        level, s, moving_sum = carry\n",
    "        season = s[0] * level**pow_season\n",
    "        exp_val = level + coef_trend * level**pow_trend + season\n",
    "        exp_val = jnp.clip(exp_val, a_min=0)\n",
    "        # use expected vale when forecasting\n",
    "        y_t = jnp.where(t >= N, exp_val, y[t])\n",
    "\n",
    "        moving_sum = (\n",
    "            moving_sum + y[t] - jnp.where(t >= seasonality, y[t - seasonality], 0.0)\n",
    "        )\n",
    "        level_p = jnp.where(t >= seasonality, moving_sum / seasonality, y_t - season)\n",
    "        level = level_sm * level_p + (1 - level_sm) * level\n",
    "        level = jnp.clip(level, a_min=0)\n",
    "\n",
    "        new_s = (s_sm * (y_t - level) / season + (1 - s_sm)) * s[0]\n",
    "        # repeat s when forecasting\n",
    "        new_s = jnp.where(t >= N, s[0], new_s)\n",
    "        s = jnp.concatenate([s[1:], new_s[None]], axis=0)\n",
    "\n",
    "        omega = sigma * exp_val**powx + offset_sigma\n",
    "        y_ = numpyro.sample(\"y\", dist.StudentT(nu, exp_val, omega))\n",
    "\n",
    "        return (level, s, moving_sum), y_\n",
    "\n",
    "    N = y.shape[0]\n",
    "    level_init = y[0]\n",
    "    s_init = jnp.concatenate([init_s[1:], init_s[:1]], axis=0)\n",
    "    moving_sum = level_init\n",
    "    with numpyro.handlers.condition(data={\"y\": y[1:]}):\n",
    "        _, ys = scan(\n",
    "            transition_fn, (level_init, s_init, moving_sum), jnp.arange(1, N + future)\n",
    "        )\n",
    "    if future > 0:\n",
    "        numpyro.deterministic(\"y_forecast\", ys[-future:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `level` and `s` are updated recursively while we collect the expected value at each time step. NumPyro uses [JAX](https://github.com/google/jax) in the backend to JIT compile many critical parts of the NUTS algorithm, including the verlet integrator and the tree building process. However, doing so using Python's `for` loop in the model will result in a long compilation time for the model, so we use `scan` - which is a wrapper of [lax.scan](https://jax.readthedocs.io/en/latest/_autosummary/jax.lax.scan.html#jax.lax.scan) with supports for NumPyro primitives and handlers. A detailed explanation for using this utility can be found in [NumPyro documentation](http://num.pyro.ai/en/latest/primitives.html#scan). Here we use it to collect `y` values while the triple `(level, s, moving_sum)` plays the role of carrying state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another note is that instead of declaring the observation site `y` in `transition_fn`\n",
    "\n",
    "```python\n",
    "numpyro.sample(\"y\", dist.StudentT(nu, exp_val, omega), obs=y[t])\n",
    "```\n",
    "\n",
    ", we have used [condition](http://num.pyro.ai/en/stable/handlers.html#numpyro.handlers.condition) handler here. The reason is we also want to use this model for forecasting. In forecasting, future values of `y` are non-observable, so `obs=y[t]` does not make sense when `t >= len(y)` (caution: index out-of-bound errors do not get raised in JAX, e.g. `jnp.arange(3)[10] == 2`). Using `condition`, when the length of `scan` is larger than the length of the conditioned/observed site, unobserved values will be sampled from the distribution of that site."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we want to choose a good value for `seasonality`. Following [the demo in Rlgt](https://github.com/cbergmeir/Rlgt/blob/master/Rlgt/demo/lynx.R), we will set `seasonality=38`. Indeed, this value can be guessed by looking at the plot of the training data, where the second order seasonality effect has a periodicity around $40$ years. Note that $38$ is also one of the highest-autocorrelation lags."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lag values sorted according to their autocorrelation values:\n",
      "\n",
      "[ 0 67 57 38 68  1 29 58 37 56 28 10 19 39 66 78 47 77  9 79 48 76 30 18\n",
      " 20 11 46 59 69 27 55 36  2  8 40 49 17 21 75 12 65 45 31 26  7 54 35 41\n",
      " 50  3 22 60 70 16 44 13  6 25 74 53 42 32 23 43 51  4 15 14 34 24  5 52\n",
      " 73 64 33 71 72 61 63 62]\n"
     ]
    }
   ],
   "source": [
    "print(\"Lag values sorted according to their autocorrelation values:\\n\")\n",
    "print(jnp.argsort(autocorrelation(y_train))[::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us run $4$ MCMC chains (using the No-U-Turn Sampler algorithm) with $5000$ warmup steps and $5000$ sampling steps per each chain. The returned value will be a collection of $20000$ samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                      mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "      coef_trend     32.33    123.99     12.07    -91.43    157.37   1307.74      1.00\n",
      "       init_s[0]     84.35    105.16     61.08    -59.71    232.37   4053.35      1.00\n",
      "       init_s[1]    -21.48     72.05    -26.11   -130.13     94.34   1038.51      1.01\n",
      "       init_s[2]     26.08     92.13     13.57   -114.83    156.16   1559.02      1.00\n",
      "       init_s[3]    122.52    123.56    102.67    -59.39    305.43   4317.17      1.00\n",
      "       init_s[4]    443.91    254.12    395.89     69.08    789.27   3090.34      1.00\n",
      "       init_s[5]   1163.56    491.37   1079.23    481.92   1861.90   1562.40      1.00\n",
      "       init_s[6]   1968.70    649.68   1860.04    902.00   2910.49   1974.42      1.00\n",
      "       init_s[7]   3652.34   1107.27   3505.37   1967.67   5383.26   1669.91      1.00\n",
      "       init_s[8]   2593.04    831.42   2452.27   1317.67   3858.55   1805.87      1.00\n",
      "       init_s[9]    947.28    422.29    885.72    311.39   1589.56   3355.27      1.00\n",
      "      init_s[10]     44.09    102.92     28.38   -105.25    203.73   1367.99      1.00\n",
      "      init_s[11]     -2.25     52.92     -2.71    -86.51     72.90    611.35      1.01\n",
      "      init_s[12]    -12.22     64.98    -13.67   -110.07     85.65    892.13      1.01\n",
      "      init_s[13]     74.43    106.48     53.13    -79.73    225.92    658.08      1.01\n",
      "      init_s[14]    332.98    255.28    281.72    -11.18    697.96   3685.55      1.00\n",
      "      init_s[15]    965.80    389.00    893.29    373.98   1521.59   2575.80      1.00\n",
      "      init_s[16]   1261.12    469.99   1191.83    557.98   1937.38   2300.48      1.00\n",
      "      init_s[17]   1372.34    559.14   1274.21    483.96   2151.94   2007.79      1.00\n",
      "      init_s[18]    611.20    313.13    546.56    167.97   1087.74   2854.06      1.00\n",
      "      init_s[19]     17.81     87.79      8.93   -118.64    143.96   5689.95      1.00\n",
      "      init_s[20]    -31.84     66.70    -25.15   -146.89     58.97   3083.09      1.00\n",
      "      init_s[21]    -14.01     44.74     -5.80    -86.03     42.99   2118.09      1.00\n",
      "      init_s[22]     -2.26     42.99     -2.39    -61.40     66.34   3022.51      1.00\n",
      "      init_s[23]     43.53     90.60     29.14    -82.56    167.89   3230.17      1.00\n",
      "      init_s[24]    509.69    326.73    453.22     44.04    975.15   2087.02      1.00\n",
      "      init_s[25]    919.23    431.15    837.03    284.54   1563.05   3257.27      1.00\n",
      "      init_s[26]   1783.39    697.15   1660.09    720.83   2811.83   1730.70      1.00\n",
      "      init_s[27]   1247.60    461.26   1172.88    544.44   1922.68   1573.09      1.00\n",
      "      init_s[28]    217.92    169.08    191.38    -29.78    456.65   4899.06      1.00\n",
      "      init_s[29]     -7.43     82.23    -12.99   -133.20    118.31   7588.25      1.00\n",
      "      init_s[30]     -6.69     86.99    -17.03   -130.99    125.43   1687.37      1.00\n",
      "      init_s[31]    -35.24     71.31    -35.75   -148.09     76.96   5462.22      1.00\n",
      "      init_s[32]     -8.63     80.39    -14.95   -138.34    113.89   6626.25      1.00\n",
      "      init_s[33]    117.38    148.71     91.69    -78.12    316.69   2424.57      1.00\n",
      "      init_s[34]    502.79    297.08    448.55     87.88    909.45   1863.99      1.00\n",
      "      init_s[35]   1064.57    445.88    984.10    391.61   1710.35   2584.45      1.00\n",
      "      init_s[36]   1849.48    632.44   1763.31    861.63   2800.25   1866.47      1.00\n",
      "      init_s[37]   1452.62    546.57   1382.62    635.28   2257.04   2343.09      1.00\n",
      "        level_sm      0.00      0.00      0.00      0.00      0.00   7829.05      1.00\n",
      "              nu     12.17      4.73     12.31      5.49     19.99   4979.84      1.00\n",
      "    offset_sigma     31.82     31.84     22.43      0.01     73.13   1442.32      1.00\n",
      "      pow_season      0.09      0.04      0.09      0.01      0.15   1091.99      1.00\n",
      "  pow_trend_beta      0.26      0.18      0.24      0.00      0.52    199.20      1.01\n",
      "            powx      0.62      0.13      0.62      0.40      0.84   2476.16      1.00\n",
      "            s_sm      0.08      0.09      0.05      0.00      0.18   5866.57      1.00\n",
      "           sigma      9.67      9.87      6.61      0.35     20.60   2376.07      1.00\n",
      "\n",
      "Number of divergences: 4568\n",
      "CPU times: user 1min 17s, sys: 108 ms, total: 1min 18s\n",
      "Wall time: 41.2 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "kernel = NUTS(sgt)\n",
    "mcmc = MCMC(kernel, num_warmup=5000, num_samples=5000, num_chains=4)\n",
    "mcmc.run(random.PRNGKey(0), y_train, seasonality=38)\n",
    "mcmc.print_summary()\n",
    "samples = mcmc.get_samples()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forecasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given `samples` from `mcmc`, we want to do forecasting for the testing dataset `y_test`. NumPyro provides a convenient utility [Predictive](http://num.pyro.ai/en/stable/utilities.html#numpyro.infer.util.Predictive) to get predictive distribution. Let's see how to use it to get forecasting values.\n",
    "\n",
    "Notice that in the `sgt` model defined above, there is a keyword `future` which controls the execution of the model - depending on whether `future > 0` or `future == 0`. The following code predicts the last 34 values from the original time-series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictive = Predictive(sgt, samples, return_sites=[\"y_forecast\"])\n",
    "forecast_marginal = predictive(random.PRNGKey(1), y_train, seasonality=38, future=34)[\n",
    "    \"y_forecast\"\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get sMAPE, root mean square error of the prediction, and visualize the result with the mean prediction and the 90% highest posterior density interval (HPDI)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sMAPE: 63.93, rmse: 1249.29\n"
     ]
    }
   ],
   "source": [
    "y_pred = jnp.mean(forecast_marginal, axis=0)\n",
    "sMAPE = jnp.mean(jnp.abs(y_pred - y_test) / (y_pred + y_test)) * 200\n",
    "msqrt = jnp.sqrt(jnp.mean((y_pred - y_test) ** 2))\n",
    "print(\"sMAPE: {:.2f}, rmse: {:.2f}\".format(sMAPE, msqrt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot the result to verify that we get the expected one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAEICAYAAAByPazKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAB3QUlEQVR4nO29d3hkV324/57poxn1stJK25vL2uuy7mDABkzH1JjegvMlkIQkhJpfKhAgpEAIECAEEwhgik21Ke7d7Hp37W1eaYtWq96naeo9vz/uvaM7oxlpurSr8z6PHs3cMvfMmZn7OZ8upJQoFAqFQqFY+diWewAKhUKhUCgKQwlthUKhUCjOEpTQVigUCoXiLEEJbYVCoVAozhKU0FYoFAqF4ixBCW2FQqFQKM4SlNBW1AwhREgIsbkKr/t3QojvVPp1ixzDRiGEFEI4lnMc1UIIcZcQ4h2L7P+WEOKTtRzTSkAI8XwhxJkCj130eyqEcAshDgshOis3wrMDIcTFQohHl3scZwNKaK9whBCnhBBzhsAz/9Yu97iWQghxvxDiD63bpJR+KeWJ5RrTSqGYG/1KuY6U8qVSytuM132nEOLhMsf2HiHEUSFEUAgxKoT4pRCi3rJ/txDiF0KIaSHEjCHMPiWEaBZCHLL8FlJCiKjl+cfLfa/LyK3Ag1LKEQAhRJMQ4jYhxJjx93fWg42F4n1CiIgxly+07NtlzNOEEOLPLdudQognhBDr8g0i3wLUujCzHGPO+ykhxEctx0ohRNjYNymEuEcI8QdZr5e+R0gpnwZmhBCvLGHeVhVKaJ8dvNIQeObfUDEnn6van+LsRAjxPODTwJuklPXA+cDtlv3XAvcDjwDnSSmbgJcASWCXlPJC87cAPAR8wPLb+HRt301F+SPgfy3P/w2oAzYCVwJvE0K8y7L/e8A+oBX4BPAjIUS7se+fgA8Bu4C/tmjvfwH8WEo5UKExNxmfw5uAvxFCvMSyb5exbwfwLeBLQoi/XeS1vos+B4rFkFKqvxX8B5wCXphjuxv4d2DI+Pt3wG3sez5wBvgIMIJ+I7ABHwWOA5PoN8kWy+s9B3gUmAEGgHca21+OfmMIGNv/znKOB/iO8XozwO+BNcCngBQQBULAl4zjJbDVePwt4D+BXwJB4Algi+W1Xww8C8wCXwYeAP4wzxz9HfAd4/EvgT/J2v80cLNlDP8P6AWmjTEIY99XgB9ZzvsscI+5P+s17cDngQngBPB+47Udxv53AUeM93YC+CNjuw+YAzRjbkLAWvSb8mPGPA4DXwJcxjkC/QY+ZszH08BOy/fg88BpYBT4KuDNd52s97DJuJ7NeP4NYMyy/zvAB43H9wN/iC5go8bnGwJmCvk8s677IeDORb7zDwP/UeDv4/5834us78cPjfcTBJ4BtgMfM+Z0AHix5fi1wM+AKaAPeK9ln9d4r9PAYeCvgDNZ5/4YGAdOAn+a63uaY4zrjc/LYdk2AVxhef5x4CHj8XYgBtRb9j8E/D/j8RHm7wePG9+v9cCTgHOJ+dqI5bts2f4t4JP5jkH//X8o+7du2f9647vTmuuzA7qNOXCXe988l/+Upn328gngauAS9NX0lcBfW/Z3Ai3ABnSz258CNwPPQ7+xmAILIcR64C7gP4B24zX3G68TBt4ONKEL8PcJIW429r0DaATWoa/2/x8wJ6X8BJka0AfyvIc3AX8PNKPfHD9ljKcN+BH6TbUVXXhfW+C83Aa81XwihNiFfjP4leWYVwBXoM/bG4GbjO1/CVxsmH+fC7wHeIc07ihZvNd4nUuB3eg3JCtjxv4GdAH+b0KIy6SUYeClwJDMtJykgD8H2oBrgBuBPzZe68XA9eg36ibgD9AXSqAvLLajf2Zbjff6N4tcJ42U8iT6YuxSY9NzgZAQ4nzj+fXoiyXrOUfQP+fHjNdssuzO+Xnm4AngJiHE3wshrhNCuM0dQgif8f5/nOfcUnkl+uK1GX0R+mv0hWw38A/Af1mO/R76onct+uf6aSHEjca+vwW2GH83of8GzLHbgJ8DB4zXvRH4oBDiJpbmIuCElDKZtV1kPd5pPL7QOD5o2X/A2A5wEHixEKIHXcAeB74IfFhKmShgPAUjdK4zrr1vkUN/CjjQ71ULkFIOAgl0zVyRByW0zw7uNPx6M0KIO41tbwH+QUo5JqUcR79Zvs1yjgb8rZQyJqWcQzc7fUJKeUZKGUNf9b/eMJ2/BfidlPJ7UsqElHJSSrkfQEp5v5TyGSmlJnW/0/fQBT/oP7BW9BV1Skq5V0oZKOJ9/URK+aRxo/ouuuABeBlwSEr5E2PfF9EtBoXwU2CbEGKb8fxtwA+klHHLMZ+RUs5IKU8D95nXlVJG0AX+v6JrZX8ipcznE34j8O9SygEp5RS6OTKNlPKXUsrjUucB4DfoQjEnxtw9LqVMSilPoQsR6zzXA+eha/1HpJTDQgiBvnj4cynllHED/zRwy5KzNM8DwPMs5tMfGc83oS84DhTxWvk+z+z3+hDwWuAydM18Ugjxr0IIO7pQtWH5vIUQnzO++2EhxF/nes0CeEhK+WtjbD9EX5x+xhBg3wc2Gj7kdehWp49IKaPG7+AbzP+23gh8ypjvAfTvpskVQLuU8h+klHGpx298ncI+jyZ0K4CVu4GPCiHqhRBbgXejm8sB/OhWFyuz6N8T0K0Z70O3GPw5cJ3x+ieEED8VQjwghHjDEmOasNx3ZoA35zoG3SLxDeCjUsp78r2YMdcT6MpEPoLoc6HIg/J1nh3cLKX8Xda2tUC/5Xm/sc1kXEoZtTzfANwhhNAs21Lo5ux16CvxBQghrgI+g77Cd6GbY39o7P5f49zvCyGa0AXdJ4pYyVsFcQT9RmS+t7TPTUopCw2oklLGhBC3A28VQvw9uvaXrQXnuy5SyieFECeADix+1hxkjJHMzwIhxEvRtbLt6EKoDt0smxMhxHb0xcJu41gHsNcY071CiC+hW0bWCyHuQL8pe4xj9+ryW38pdNN9oTwAvApds3wQ3WT5NnQz5kNSSi3/qQvIO6/ZSCnvAu4ytNMXoH+nnkX/DmlAF3DUOPbDwIeFHnld6j1r1PJ4DpiQUqYszzHGuxYwF0Am/eifCyz+uW8A1hoCzsSObnVaimnmBa7Jn6Jbv3rRLSvfQ/8+g+6aaMg6vgFD8Esp+9EXvwgh6tBdXzcZr/cD9MXSQSHEPcaiMxdtVs1fCPGtpY5ZDCGEE32xlO96oM/BTCGvt1pRmvbZyxD6TcJkvbHNJNukOwC8VErZZPnzGCapAXRzXy7+D321vk5K2YjuMxWgr5yllH8vpbwA3Xz9CnRTeq7rF8Mw0GM+MTTKnvyHL+A2dOvBjUBESvlYoScKId6PvjAZAj68xBitEbjrLa/hRjfvfh5YY5iQf8W8qTPX3HwFXUhtk1I2oPsv05JYSvlFKeXl6CbI7ei+1Al0gXOh5TNtlHrwT77rZPMAugXg+cbjh9G1sueRZRq3ULHWgIYF5x7gXnQ/fRjdfP7aSl2jSIaAFmGJZEf/bAeNx3k/d/Tf0cms31i9lPJlBVz3aWCzNWjU0ObfIqXslFJeiH6/ftLYfcg43jrOXcb2bP4G+IaUchTdDL9HSjmLvlDbWsDYKsWr0YMJn8y1U+hZMS70xZsiD0pon718Dz0qtN3wAf8NupaSj68CnxJCbAAwznu1se+7wAuFEG8UQjiEEK1CiEuMffXomkdUCHElFhOZEOIFQoiLDLNmAN2Ma2ovo0CpOdm/BC4SQtxs3MTej+6jLwhDSGvAv5AZjbsohrb7SXQT+dvQtbtL8hx+O/CnQogeIUQzepCfiWmRGAeShtb9Ysv+UaBVCNFo2VaPPochIcR56KZNc1xXCCGuMjSVMEYgmKEFfx3dX95hHNtt8aHmuk4GUspedMH/VvR0o4Bx3uvIL7RHgR4hhCvf6y6GEOLVQohbhJ6+JYzv1fPQA6ZAXyy9WwjxUcv76kEPnKsqhsn7UeCfhBAeIcTF6LEN3zUOuR34mDH2HuBPLKc/CQSEEB8RQniFEHYhxE4hxBUFXPcMukad9vcKIbYYv0W78R26Ff37iZTyGHrcyd8a43wNcDFZsQBCiAvQF2RfMTadBG4QQqwBtqEHMFYVIUSLEOIt6Jaiz0opJ/Mc+nzgXsN9p8iDEtpnL58E9qCv0J8BnjK25eML6Brzb4QQQfQb5FUAhm/3ZeiBWFPoN4Ndxnl/DPyDcc7fkGky7kT3gQbQo1UfYH7h8AV0n/m0EMLq91sSKeUE8Abgc+hmwQuM91rMj/nb6FpFQUVXjMXBd9BvKgcMYfZx4H+FJVDKwtfRg5kOoM/9TyzjD6KbNm9HN3u+GX3uzf1H0RddJwx/4Vp0c/eb0c2bX0c3YZo0GNum0c2xk+haPOgZAn3A40KIAPA7jECePNfJxQPApPE9MJ8L8gcV3Yuu0Y0IISbyHLMY0+i++F707853gH+WUn7XGPfDwA3ogXDHDHPz3eim+/8o4XrF8ib04K0h4A702JDfGvv+Hv0zOIkep5BeFBrm9lei+/JPoltCvoEerFkI/0VmXMrl6L/tIHrMxFuklFZN+hZ0s/00ugvr9VKPb7Hyn8CfWVwBH0P/bh4CPi2NnPAqcUAIEUL/fv4heuzF3yxy/FvQlQvFIpipLgrFisXwe55Bv2ndV+A5bwdulVI+p6qDUygqhLE43AfcKKUcXu7x1BIhxEXA16SU1yz3WFY6SmgrViSGifcJdNPtX6GbyDdLPRJ+qXPr0LXBL0spv13VgSoUCkUNUeZxxUrlGvSI9gl0k+PNBQrsm9B9yaPoQXQKhUJxzqA0bYVCoVAozhKUpq1QKBQKxVnCii+u0tbWJjdu3Ljcw1AoFAqFoibs3bt3QkrZnmvfihfaGzduZM+ePcs9DIVCoVAoaoIQoj/fPmUeVygUCoXiLEEJbYVCoVAozhKU0FYoFAqF4ixBCW2FQqFQKM4SlNBWKBQKheIsQQlthUKhUCjOEpTQVigUCoXiLEEJbYVCoVAozhKU0FYoFIoVRkpTPSEUuVFCW6FQKFYY0URquYegWKEooa1QKBQriJQmiSe15R6GYoWihLZCoVCsIJKaRiKlhLYiN0poKxQKxQoipUniSmgr8qCEtkKhUKwgUpokkVKBaIrcKKGtUCgUKwhdaCtNW5EbJbQVCoViBZFUgWiKRVBCW6FQKFYQStNWLIYS2gqFQrGCSCqftmIRlNBWKBSKFYSmNG3FIiihrVAoFCuIpBLaikVQQluhUChWEClNI5mSSKlM5IqFKKGtUCgUKwhTyVYFVhS5UEJboVAoVhBJTRfWKhhNkQsltBUKhWIFYbblTKhcbUUOlNBWKBSKFUTSFNrKPK7IgRLaCoVCsYLQDKGtfNqKXCihrVAoFCuIeU1b+bQVC1FCW6FQKFYQKWUeVyyCEtoKhUKxgjA1bdU0RJELJbQVCoViBZFKp3wpoa1YiBLaCoVCsUKQUmLIbOXTVuRECW2FQqFYIZimcVCatiI3SmgrFArFCiFlEdoq5UuRCyW0FQqFYoVgFdqplEznbCsUJkpoKxQKxQohmSWklbatyEYJbYVCoVghpLKEtvJrK7IpSGgLIf5cCHFICHFQCPE9IYRHCNEihPitEKLX+N9sOf5jQog+IcSzQoibLNsvF0I8Y+z7ohBCVONNKRQKxdmI2eHLREWQK7JZUmgLIbqBPwV2Syl3AnbgFuCjwD1Sym3APcZzhBAXGPsvBF4CfFkIYTde7ivArcA24+8lFX03CoVCcRaTJbOVpq1YQKHmcQfgFUI4gDpgCHg1cJux/zbgZuPxq4HvSyljUsqTQB9wpRCiC2iQUj4mpZTAty3nKBQKxaonW9NWVdEU2SwptKWUg8DngdPAMDArpfwNsEZKOWwcMwx0GKd0AwOWlzhjbOs2HmdvX4AQ4lYhxB4hxJ7x8fHi3pFCoVCcpSiftmIpCjGPN6Nrz5uAtYBPCPHWxU7JsU0usn3hRim/JqXcLaXc3d7evtQQFQqF4pwgO3pc+bQV2RRiHn8hcFJKOS6lTAA/Aa4FRg2TN8b/MeP4M8A6y/k96Ob0M8bj7O0KhUKhgAV52UrTVmRTiNA+DVwthKgzor1vBI4APwPeYRzzDuCnxuOfAbcIIdxCiE3oAWdPGib0oBDiauN13m45R6FQKFY9Kk9bsRSOpQ6QUj4hhPgR8BSQBPYBXwP8wO1CiPegC/Y3GMcfEkLcDhw2jn+/lDJlvNz7gG8BXuAu40+hUCgU5PBpq0A0RRZCD+ReuezevVvu2bNnuYehUCgUVWf/wAwTwVj6ucdp5znb2pZxRIrlQAixV0q5O9c+VRFNoVAoVggqelyxFEpoKxQKxQohW2inNNU0RJGJEtoKhUKxQsgurgIqGE2RiRLaCoVCsULI1rRBmcgVmSihrVAoFCuE3EJbmccV8yihrVAoFCsEpWkrlkIJbYVCoVgBpDRJrgxc1TREYUUJbYVCoVgB5NKyQWnaikyU0FYoFIoVQH6hrXzainmU0FYoFIoVQK50L1CatiITJbQVCoViBZBP085uIqJY3SihrVAoFCuAfEI7lUcDV6xOlNBWKBSKFUB+oV3jgShWNEpoKxQKxQrANIOPBqKcnopYtiuprZhHCW2FQqFYAZia9k+eGuQ/7u3FbJucTwNXrE6U0FYoFIoVgKlph2JJpiMJRgN6X20ltBVWlNBWKBSKFYApnGPJFABHhgMAaLnKpClWLUpoKxQKxQrAFNrRhO7DPjoSBEDTUD21FWmU0FYoFIoVgBlwFk3omvbRkUBay1a52goTJbQVihoSjCZ4xX88lDZ9KhQmaU07maLe4yAcTzFgRJErv7bCRAlthaKG9E9GODgYYN/pmeUeimKFoXf5ksQSGrt6moB5E3lK+bUVBkpoKxQ1JBxLAjAdiS/zSBQrjZQmiSc1JLCmwU1noydtkUmppiEKAyW0FYoaEo7rQnt2LrHMI1GsNJKaJGr0zvY47ZzfWU/vWIhkSlMFVhRplNBWKGpIKKYHGc0oTVuRRUqT6SA0j9PO+V0NxJIaJyfCyjyuSONY7gEoFKuJUNQ0jytNW5FJhtB22Ni2ph4BHBkJqkA0RRqlaSsUNcT0ac8qoa3IQhfa8+Zxv9vB+tY6jo4ElNBWpFFCW6GoISFDaM/MKfO4Yh4pJSlNEkumaCJIszYBwPmdDRwfD6ctNAqFEtoKRQ2Zjx5XmrZiHms1tB+4/pE3PPEGbMko53XWk9Ik+8/MLO8AFSsGJbQVihoSspjHpQouUhiYFc+i8SRbxBDuZBBPZIiNbT4ATk2El3N4ihWEEtoKRQ0xhXY8pTFnBB0pFKamLeMhHEL3a7ujY9S57AhUiqBiHiW0FYoaYprHQZnIFfOYmraIzZe3dc+NYROCOpedwJzyaSt0ChLaQogmIcSPhBBHhRBHhBDXCCFahBC/FUL0Gv+bLcd/TAjRJ4R4Vghxk2X75UKIZ4x9XxRCiGq8KYVipRKOzWvXKldbYWJ28bLHZtLb3HNjAPjcDoJRtcBT6BSqaX8BuFtKeR6wCzgCfBS4R0q5DbjHeI4Q4gLgFuBC4CXAl4UQduN1vgLcCmwz/l5SofehUJwVBGNJ6j16eQSV9qUwMTVtZyJT0wZdaAdU9LjCYEmhLYRoAK4H/htAShmXUs4ArwZuMw67DbjZePxq4PtSypiU8iTQB1wphOgCGqSUj0k9AufblnMUilVBOJaku8kLKPO4Yp5UWmgH09tc0VEAfC670rQVaQrRtDcD48D/CCH2CSG+IYTwAWuklMMAxv8O4/huYMBy/hljW7fxOHv7AoQQtwoh9ggh9oyPjxf1hpYbKSXJlKoTrMhNOJakp1kX2ipXW2Fi1hZ3JWbT26yadiimNG2FTiFC2wFcBnxFSnkpEMYwhechl59aLrJ94UYpvyal3C2l3N3e3l7AEFcO//zrZ3ndVx5d7mEoVighi6Y9ozRthYHZD8SdDKW3pYW2y6GKqyjSFCK0zwBnpJRPGM9/hC7ERw2TN8b/Mcvx6yzn9wBDxvaeHNvPKZ4dCXLgzKxaGSsWkEhpxJIabX43HqdNpfEo0piatic1bx53R8dASnxuO5F4SpUyVQAFCG0p5QgwIITYYWy6ETgM/Ax4h7HtHcBPjcc/A24RQriFEJvQA86eNEzoQSHE1UbU+Nst55wzTIRiABwbDS5xpGK1YaZ7+dwOmrwupsPKPK7QMQPRvKl5TdueiuFIBPC5HUhQfm0FUHiXrz8BviuEcAEngHehC/zbhRDvAU4DbwCQUh4SQtyOLtiTwPullGaey/uAbwFe4C7j75xiIqTfiI8OB7lsffMSRytWE6b1xe920FTnZEZp2gqDmNEopE7LXOy758bwuVoAmArHaapz1XxsipVFQUJbSrkf2J1j1415jv8U8Kkc2/cAO4sY31mFlDKtaR8dCSxxtGK1YeZo+wyhrVK+FCbxlPHdkOGM6B/33Bg+tx7XMxWOs/nsCvFRVAFVEa2ChOMpYkl9xXx0RJnHFZmEYrqQ9nsM87gqrqIwiCU1pJTUS908Hq3rAsAdHaXOpetW6vuiACW0K8pEUNey/W4HR4cDqiGEIoOQoWn73XZlHldkEEvqQYr1RAAINWwFTE1br001pWIgVgyJZUzrVUK7gkyGdaF99eYWAtEkI4HoMo9IsZLICESrc6lOXwpAL6ySSkliSY1GoXfzCluFtqFpqxTBlcNyWj2U0K4g40H9g7xuaxugB6MpFCZmIJrPpfu04ymNSFx1+lrtxA2XWjSRooGFQrvO0LRXs9CeW2G/k+W0eiihXUFMTTsttJVfW2HBLJBR73HQ5HUCKBO5Ii20Y7EoPhFDw064fhMArugYDpsNj9O2qivoBVZQuls8qRFcxmI3SmhXkEkj3Wtjq4+1jR4VQa7IINs8DqrTlwJiSV2LlHP6/SLm8BPzrgEyq6Kt5myDlVSIKJZMLWvxLCW0K8hEKEaj14nLYeO8rgaeVZq2wkIonsTlsOG022iqMzTtVXwjVuiYGSciNg1A3NlA3Ku3cnBHx0Fq+NyOFSW4aomUksAKeu/xpEYqJZfNZK+EdgWZDMVp9esa1I7OevrGQmnTl0IRjiXxu/WgIiW0FSZxIxJZRPVmIQlnA5rdTdzVhE1L4IzN4HPZV63QTqTkior9MD+v5dK2ldCuIOOhGG1+NwDnddaT1CQnJkJLnKVYLYSi80K72TSPr2I/pULHrIZmj+nm8aSrHsCibY+u6p7a8ZSma7crpPa6qYgpoX0OMBmK0WZo2ud1NgAqglwxTyiWwmcI7Uav0rQVOqbmZjfacqZcjQDEPIbQnhszhPbq/K4kDCE5l1gZ2rbpzliuzmtKaFeQiVCcVp+uaW9u9+G0CxVBrkijm8f19B2P065HBKtAtFWPqbk54rqmLd2G0PZahLbLTjCaXJV5/eaiJhJfGZYGpWmfI8STGrNzibR53Gm3saXdryLIFWnC8WRa0wbdRK40bYUZPe6Mm+bxHELb7SClScIryLdbK9J57PGVER9katqReBJtGUz2SmhXCLNCjhmIBnB+V4MyjyvSWH3aoJvIp5XQXvWYJTFdSf1eYfq0583jo5aqaKvPMmNq2ivHPG6k6El9IV5rlNCuEONG3XFT0wY9gnwkEF2VPzTFQkKxTKHdVOdkVgWirWriSQ3NUCA9htBOOLM07eh8/fHVaJmJrzCftjUjaDlM5EpoV4hJo6xdm0XTPq9TXzErv7YCdJ+2Mo8rrMQtjSc8KVPT1oNYrQVWzO/Nakz7Sqwgn7amSZKpeZN4WAnts5eJHJq2GUGuiqwoNMMf6cvStJV5fHUTs2iPdZqeHppIC+15n3adS2naZmrcchLP6u61HOVMldCuEGbdcatPe02DG4dNMKq6fa16TN9XfYZP28XsXHxVRgQrdKxCwBTaSacutOMevYeBKzaBX88QXJV5/eYcpTRJdJlN5LGsYlnKPH4WMxGK43LYMnyWQgilTSkACBu9tLM17ZVW7UlRW6z+UZ8ptA1NW9qcxNytCKnRLPUc7tWsaQMrQGhnXj+W0GreW1sJ7QoxEYrR7ncjhMjY3lTnUoFoivm2nEZAEUCzUcp0OXvzKpYXq+bmN9pymoFoMF8VzZ+YwGkXK6oGdy2QMtOHvNzBaLnKUtfar62EdoWw1h230lznVDdlRfqH7c8yj8Pq1J4UOmkhICX1UhfaZsoXZBdYcay6e0m2D3m5rVK5hHat/dpKaFeICUvdcStNKkJYgVXTzjSPw+qMCFbomJq2LRHCITRiwoO0Oef3Z5UyXW33kmwhuVydtUyyfdpQe7+2EtoVYjIUp9WnNG1FbkI5NG2zaYj6fqxeTB+pNjcDwJy9PnN/Vq72avuuJFKZQZrL7dNW5vFzBCklk+EYrTk07eY6F9ORhIoQXuXkMo+r9pwKUwhIQ2hHHdlC28zV1quirTarzAJNe7mFdo6gs2SNS5kqoV0BAnNJEimZUVjFpKnORTypLfuXTbG8hHOYx81OX6vtRqzQsRbqkEYv7ZijIeOYuLsFAFdsWpnH0aO1l6Pet0kuTbvWKKFdAcZDCwurmMxHCK+uH5sik2AOTdvjtON12pkOry6Tp0LHqrWJ6AwACWempm02D3EkAvhc9lUXPZ5Ls11OBUgJ7XOEyUWEdpPpt1Q35lVNOJbEbhN4nJk/uaY6JzOr7EacTd9YcFX+PqwVvmyGpp1wZmraZnU0RzyAz+0gmtRq7tetdR6ylVxCcrkiyBMpjdQyavkmSmhXALPueL6UL1B+y9VOOJbC57IvyONv9DpX/XfjTV9/gi/c07vcw6g5sdS88LEZbTlTrkyhbVZHc8Znl63++Inx8LJpmLkWDMsVjLYStGxQQrsiTIQWljA1afapCGHFwg5fJs2rvPhOJJ5kPBhjcGZuuYdSc6xCwBHXNe2UuzHjmETaPB7Etwz1xyPxJIMzkZxm6lqQS2gvl3lcCe1ziIlQHCGgpS5XIJqpaa/eG7NC76XtyyG0V7t5fGRWr8tvtrZdTcQyhLauacssoZ1y1CGFDUcyTL1Zf7yG95Le0RCaltnYpJbkEpTLlau9XAuXbJTQrgAToRjNdS4c9oXT2eQ1Ne3Ve2NW6A1D/J48QnsVL+hMoW1aq1YTVoHkMnppS0+m0EbY0n7uJptujajVIm86HE8vppZLYOW67nL5tFdClzFQQrsiTIZiOdO9gHQTEWUeX93kM483eJwElqG930ph2KJpr7ZaBlah7U7omjaepgXHJY2I8kabXuZ0tkYKwLHR+ZbCy2Eazu5dnR7Lsi0gVkbarhLaFWAiFKfVtzBy3ETXppSmvZoJx5L4XAuFts/tIJ6sfaeglcKI0bY2ltSWpc3hcmI1j7tTuoDM9mnDfNpXAxGgNu05h2fnMmpq5yrfWW3yCefkMv1WlmMOclGw0BZC2IUQ+4QQvzCetwghfiuE6DX+N1uO/ZgQok8I8awQ4ibL9suFEM8Y+74oskNpz1ImQzHa6vML7RafS2naq5x8Pu06I7goElsZq/haMzw7H4C22vzaVu3Va5jHs1O+YD7ty6cFsYnaBKINTGUGBi6Hpp1PaEu5PIL7rBPawJ8BRyzPPwrcI6XcBtxjPEcIcQFwC3Ah8BLgy0IIsx/hV4BbgW3G30vKGv0KIV/dcZMmo5SpYvUSiiWpz+HTNk3mofjq0jJNTJ826BarlcB4MFYToWA1t3o1s8PXQqFtpn25EgHqPc6qp3xFE6kFRVyWQ2AlFrlmrUuHwlkWPS6E6AFeDnzDsvnVwG3G49uAmy3bvy+ljEkpTwJ9wJVCiC6gQUr5mNSdV9+2nHPWEk2kCMaSeX3aoOdqr+Zgo9WOlJJwPJXRS9vE1L4jq8w0bDI8G2VtowdYGZp2SpPc+C/386X7+qp6nURKQ7PIAJ+ma9rJRTRtRyKI3+2oeiBarqBAs7FJLVnMd70c7qSzSmgD/w58GLCOeo2UchjA+N9hbO8GBizHnTG2dRuPs7cvQAhxqxBijxBiz/j4eIFDXB7Mwiq5qqGZNNe5VmXFJ4VOLKlXUsplHjcF+Wrz55qMzEa5qEf32a6ECPJQNEkgmuTOfYNVDYyzaq4iFcdDjBQ2Uo66BceagtwRn8Xntlc9EC3X4mk5BFYimX/+s7t/VRsp5YqJO1lSaAshXgGMSSn3FviaufzUcpHtCzdK+TUp5W4p5e729vYCL7s8jBmBNLk6fJk01ekRwssVQKFYXsyAnlzR42ZwWngV+rRjyRST4TjndzVgt4kVoWmbpudTkxGODAeXOLp0MgqrJPTrhIQfcoT5mCZzZ0IvZVrNQLSUJnPG3yRTsuaNOhaL1k5qtb2XxlMaKyW5oRBN+zrgVUKIU8D3gRuEEN8BRg2TN8b/MeP4M8A6y/k9wJCxvSfH9rOakxO6L2pT28IVsonZN1l1c1qd5GrLaWJq3+FV6NMeC+hCem2Tl1afa0Vo2oHo/G/0roPDVbuOVWg7jXSviM2f89hEWtMO4HNVt9PXZChGPnlYa792fAVp2ivFNA4FCG0p5ceklD1Syo3oAWb3SinfCvwMeIdx2DuAnxqPfwbcIoRwCyE2oQecPWmY0INCiKuNqPG3W845a+kbC2G3Cda3+PIe06Q6fa1qQjnacpqkhfYqNI+bOdpdjR7a/O4VpWnXux388pnhqpnIrT5isxpaxFaf89hk2qcdoM5VXfP4eChGNJHiod7xBe89l+AKRBP0T4arMpbFfNq1tlqulMhxKC9P+zPAi4QQvcCLjOdIKQ8BtwOHgbuB90spzW/o+9CD2fqA48BdZVx/RdA3FmJDax0uR/6pNDVtFYy2Ollc07ZnHLOaMNO9uho9tNe70y1ulxMzavrVl67lxHiYY6Ohqlwn0zyuC+05Rx5N28jTdhqdvoKxZFU0PyklE6E4j52Y5LbH+ukbz3zvsRzm6s/cdZQ3f/2Jio8FFg82U5p2gUgp75dSvsJ4PCmlvFFKuc34P2U57lNSyi1Syh1Syrss2/dIKXca+z4gz4ESSMfHQ2xtz/1jMzGFttK0VyeLadr+tHl89fm0zXSvzkYvbX43EytI037D5esQAn71THVM5Bl1x2N6s5CYPY+m7ZzXtM3vy2S48nM1E0mQSGoMz+ify7MjmT79XGU8nzw5xVgwWhWLxGKCstY+7eXs4Z2NqohWBomURv9khK0diwvttHlcRZCvSkKLaNpepx0hVqumHcXvduB3O2ivdzMRii97KVPTp72lw88VG1uq5tfOMLcaQjueI90LLIFo8QANRq7/RLDy9xIzpmDIsIBkWxmyzdWzcwn6xkIkUrIq5uPFzeO1/Z4sVzvQXCihXQb9k2GSmlxSaKv2nKsbMzI8l9AWQuBzOVZlytfIbJROI0e7ze8intIIzC3vPMzOJbDbBD6XnZdf1MWx0RB9Y5WPIrdqkSI6A0DCkVvTTlhSvhq8ugJQjaA9M6bAjDXoGw9l+I6zNe0DAzPpx8EK189PaZLUIoK51ulXSmifI/SN6SvRpYS2z2XHaRfKPL5KCafN4wuLq5jbV2MZ0+FAlC5DaLcbZYDHQ9HFTqk6gbkkDR4HQghesrPTMJGPVPw6GYFoc5MAxN1NOY81G4Y4EkEajIVfpf3/kXiSSDxFOJZkdi7B5jYf8aTGqclI+phszXff6Zn040ovOpcSymZFtIGpCK/4j4c4Mx1Z9Phyia6QDl+ghHZZHB/XoyY3L+HTFkLQVOdSgWirlKAptHM0DDG3r8YypiOzc3Q2ZAntKph9i2F2LkGjoc2uafCwe0Nzxf3a2d2rvGG95lSkLmetKaTdRcruxSZTNDv0+am0pm1ag0wt+3nb9foYi3X62jcwnX4cjFZWIVnK3G4K9cdPTHJwMMBdVVhYmUgplaZ9rtA3FqKr0ZPT7JlNc51TmcdXKXqHLzs2W+7+OD63Y9X5tJMpjfFgbF7T9pua9vIGowWiibQJGuCF56/h6EiQyQqOK1tj9c8NAjBX15PrcMDSNIQwXqe94j5tU/M3I/q3r6mnu8mbEYxmtQ5ommTf6Rk2t+mprqEKm8eXitY2Fz2m4vRgb/UqZ8aSWYVVpESklu9eroR2GfSNhZY0jZssV9OQPaemeLh3oubXVcwTjuXu8GWyGs3j46EYmtQjx2Fe017uCHKrpg2kfe6VrPedoUVKSX1UF9rxhvV5zzEjyJ3xWVqqUIjGNP8OzUZx2gWtPhc71tTTOx5KR2pbTdYnJ8PMziW43tDIK90TfqlFrDmm40Za2hMnp6qmDc9lZXact/dvuP7n1+GO6BaYp05Pc3qyuuZ5K0pol4imSY6Ph9iyhGncZLmahvztzw7x9z8/VPPrKuYJxZKLWmNWYyCatbAKQKPXidMull/TnkvQ4JkX2ubjSgZaWTVWZ3wGdypCQHpx+1vznmNtGlIdoT2vaXc2eLDZBNs7/cSTeoYMgKbNa8CmP/u529qAyvu0l3o9TdOD1U6Mh2iqcxJPajxxcmrRc0olmtUspXXkYZzxWTrO/JakpvGF3/XyzUdOVuXauVBCu0SGA1Ei8VTBmnbzMmjagWiCw8OBjPaHitoTWlLTdqy6MqbzOdq60BZCrIhc7dm5ZIZ5vMGrf27ZrSrLwWr69Rj+7DOyg3rLdbNJpguszNJU56y40Da1/+GZKF2G9WN7hx4Al8tEvu/0NPVuB5esawIq79O2LpI6zvyaSx94F87YdMYxkXiS/skIr7m0G5fDxoPHqmMit2raIhXHM6dr2K2jD3F6MkIsqXHlppaqXDsXSmiXyPECI8dNzEC0Wuah7jk1hZR6IFSlf1SKwpmJZJpcs9F92qvLPJ6taYPeKW/ZNe1oIi2oAeqroGlbhbYZhDZIB057/ttxOoI8HqDR66x47/FYIkUsoTdwWdukfyYNXidrmzw8myMYbd/pGS5Z35Re4FTSp61pkohlEbuu99u0jj5Cx5lfZxzXPxkhqUkuXNvIlRtbeKhKfm1r5LgnMoyQ+vPmsSc5MaJr91dsVEJ7xWOmexVjHk+kZE0rX1nNRUrbXj5GA1HWNHjy7ve57KsuEG1kdg6P05axmNELrCyf0I4mUsSTWsaYTPN4oIKLXqtP2xvSuxiP2dcseo5ZytSR0IX2dCRe0VzlWFJjOGAupLzp7TvW1NM3Nu/Xjqc0IvEkR0cCXLquCafdhsdpS2dIVIJQPJkR+OUJ6z7/+umDGceZ+fNb2n1cv72NY6OhdCBdJbFWQ/OG57tO21NzOAefpNMow1srlNAukb7xEI1eJ21+V0HHp0uZ1rAq2hMnpvA69dzgISW0lwVNk4wFY6xpyP+j9rkdzCVSpGrc+nA5GZ7VzbDC0oqyze9a1qYhpgnc6tOu91TXPG5q2pPOrkXPMTVtZ3yWBo8TKWGqQveSRErv926WL7VaP3asqSdm8WvHEhpPn5lFk3Dp+mYA/G5nRS0RVq1daAk8c3o6V8NUltC2pNyaAXEPHat80G3MKrRDAxn7Ns4+znmduYviVAsltEvEjBwXOfrf5sKsilbNtnpWwrEkBwdnuelCfQU/UoUVqGJpJsIxUppM+25zYQapRVaRX3tkNprO0TZpr3czGYrXvG+ziVl33Kpp17ns2G2isubx1EKhPeNeu+g585p2ML2QqNQCxxqEZheCDssCc/saXSCZ+drxlJYOQjP92Q0eR0Xdb9YgNHdkJG2O9s8ew5aaf8+nJsK0+d00ep3sWFNPR727Kqlf1kA0U9Oear8KgKvlASW0zxaOjy3dKMRKc7o9Z2007adOT5PUJK/ctRYhYGhGadrLweisfpNZ1Dyebs+5evzauqadOSdtfjdJTVY0vaoYTKFtDUQTQlDvcVTWPJ5YGIgW9OYurGJiTfkyvy+VciWY5vqh2SgdDW4ctnmxYPq17zkyxv6BGWIJjadOT7OpzZdWRPyeymY/WBdI3shg+rFNJvHPHE0/PzURZku7nicuhOC529p5uG+ioharaCKV0V/cFNrDG19DXLi50NbPrubaZgUpoV0C0+E4k+F4wUFooAeiQe2E9pMnp7DbBFdtbqXN71Y+7WVixPATZmuVVtLtOVeJpq1pUvfzNy7UtKE6dbULwRTM2UGDDZ7Kmn/jZotLqaWF0pxvCaFtSfmaF9qVuZdYNe21Fn+2ybuu3USd286X7uvj735+kL3901xqaNmguxAqah63LABMf7ZJg8WvfXoqklGN8vrtbcxEEjwzOFuxsWTnfpvm8XDDFg45LwJga/D3FbteISihXQJmQn8xQtvUtGtlHn/ixBQ71zbgdztY2+hJd+5R1Ja00F7EPG6WN10twWiT4ThJTebUtKFyZt9iSWvansz0vHqPo2I+7URKS2tu7rkxbFqCcdmIx7u4idXaNKTOWORVUtNOpDTGLBXqrGxq8/E3r7iAN+7u4dBQgKlwnEvXN6X3+92OikWPx5IpEjl8/nGXfj3Trx2MJghEk2lNG+A5W9sQAh6qYOpXRs1xKakLnwYg4uvhnsSFADQNP1Sx6xWCEtolUGzkOMyv3muhaUcTKfYPzHDVZr1YQ2ejR2nay8RYIIpNQKsvf8CiqTmtlgIr6RztHD5tWEZN2+gwVk1NOyNy3DC1Dsj2tJ/aJDtUxtqe0ykEHqetYjnt0USK0UAUKaGrKffi0mGz8eILOvn0ay/iIy85j9dcNl9ytd7jrJhPO1v4m0J7vPtF+rUMTdv8Dm2xKE6tfjc71zbyUAUrQFojxx3xWRyJEEmHj5GEj7ujutBuHHoIatjfWwntEugbC+F22OhuXmhKyofDbqPB46iJpr1/YIZ4SuNKI3ewq9Gbzotd7ewfmOG2R0/VLF9+ZDZKe70bxyI5uKZ5fLWUMjXTcrqyTLErRtPOEtqV9GnHcwrtjgyhbbMtbONqTfkSQtDqq1x6XCyppe8PuczjVhrcTm69fnPG+PxuR8VSvrIXrqZ5fLz7hWjCjj/Qiy05l7ZgbWnLVJx2djfSN57ZB7wcrIVVzM9rzr+OY2Mh+mQ3YXcHrugkjB7M9xIVRwntEugbD7G53Y89TwOIfDT7XDXRtJ88OYUQ8wn/XY0eQqrACj/ae4Y3fvUx/vZnhzgzXRt3wUhgYZR0NulAtFXi087nMmjwOHA5bMtWYCUwl6DOZV9Q5KTBWzlNO0Noh3QtckC247cIbbfDvqCC3nzKVwCAVr+roj7toZk5BIsHTJpkN/NoMALRKhH1nz3Pps8/3LCFcMNWhNSonznKiFEjPVtx6mn2MhWOV8zVlCtyfM63jt7REHUuB9Ndz9V3Hr+nItcrBCW0S+DURDjd3aYYatU05ImTk5zX2UCj4UfvatK/2KtV205pkk//6ggf+uGBtDntySrVKc5mqcIqMK9VrRbz+MmJMB6nbYHLQAhBu9+9rJq2NUfbpMHjrJhP21p33IwcH5Ad1Lvnr+t22KhzZfZeTzl8aMKOPTWHSMVpqatc/XFT026rd+NyLC0SYlm1uP0eB1JCpAINO6y/AZGK454bRQobUW8ngRY98Kth+pn0YjhbcVrXUgdQsUV5NL4wR3vOt45eI3toyhTafUpor1hSmmRwZo71rXVFn1uLpiHxpMbe/mmustTCNYNLVqPQllLyvu/s5WsPnuBtV2/gzvdfS73HwZ7+2gjtkdnookFoQPoGvVrM40/1T7Orpylnq9K2+uUT2oFo7nKz9R6933klNMlchVVOZ5nHPc6FmjZCzJcyTQRp9lWm/ngipZFKyZwpePnI1rTnS72Wt7DJLl/qmdNztGPeNUi7i0DzTv16Uwfz/q56DM37zHRlum7l0rSn3WsZCUTZtsbPdMc1pBxecPmgRi43JbSLZCQQJZGSrGsuRWhX3zz+9JkZogktt9CeWX0R5P2TEX5zeJT3v2AL/3jzTtwOO7s3NPP7U9NLn1wmc/EUgWhySU3bjB5fDZr2XDzFoaEAuzc259zf7ndXvK52oczOZdYdN2nwOtM1/MsllkNoj4gO3BYNN5emDZamIYkAjV4XU+F42TnJsaSGJvUUvK4CTOPmOVbSlqIyXQjheDIjnsv0Z8/59KC3YMu80B4PxXL+ruaFdvn3ulgyK0fb0LSPxfXOZts66km4m9nzB/vgzT9YGD1YJZTQLhKzb+r6luKFdlOdk5lwdc3j//fkaepcdq7d2pbetqbBgxCrU9PuNSL9bzx/vrbz7o0t9I2FKlYGMh+jhu92KaFtswnqVkn98f0DMyQ1ye4NuRsstNcvXynTwFwyr6at7y//t2sKPNP0m8JGyNWZUVnR7bDjczkWyIBEVtMQrQKlTKOJFDORBElN0lGi0E7PT7lCO8vSZC5q5ur0HPZg43loNif+4HE8Unc7ZVs/2v1u3A4bA1Pla9rReOb7NDXtZ8JNOGyCjYa1VdoLK2VdKZTQLpKB6dKFdnOdi2AsWdFC/1bGglF+fmCIN1zek3HzcdpttPvdVSmmv9Lpy9GNzQzQ29tfXW27kMIqJnUuR02bySwXew23xGXr82vaU0bp11qzmE8bKtPpyzQteyODCCQTtna83sy69B6nDZtN4HFmatvJjKYhlSllGktq6dfI1UfBG+xn9z1/wK6H37fgPZiYQrtcS1EolrkoMjXtqFF4RtpdhBq3I5BcKE7R1eAhkZVqJYSgp9lbEU3bahoXWgJPZBiJ4JlwI12NnkUzQqqJEtpFMjAVwSby5zMuRouvuk1DvvNYP0lN8q7rNi3Y19XoWRGa9n/c08v3njxds+v1jgVZ0+DOuBlf3NOIy25jz6nq+rVH01HSS3cA8rtXh6a9t3+abR3+dJBkNm317opokKWgt+XMJbQd6f3lYgZxmUFog6KD+iz/tduhC+tsE7m1lKnfCFwr168dTaTS0frZnapahx/gyt+9jqbJfbQP3YPDiFyPZgWiVcqnna2peyOGpu2bzwk3/doX206wpsFDMrVwcdfTXMeZmfI1bWu6lycyjE2miNV1cjqQYm1T4em+lUYJ7SIZmIqwtsm7aO/bfJgr2ckq3JCiiRTfeeI0N563ho05IttXQq724Mwc/35PLz/cM7D0wRXiuNHYxYrHaefinkZ+X2WhbRaAKMTsqPfUPreFtqZJ9vZP5/Vnw3yudq0LrKQ0STCazC20vZXRtDVNpoVMOt1Ly0z3AnA79XtLdjBawixlGg+kzyl3nmIJjYlgDCHmlQqkZOPhr3DJQ7fiTATSx9aFTqXPsVIxn3aeHO2opcSrKbQvd57C67LntFqua6mMpp2rJWe4roepcLzgoL1qoIR2kZyeipQUhAZ6xR6AySoE2ty5b5CpcJz3PGehlg0royratx87RUqTNRuHlJK+sRDbOhaWiNy9sYVnBmcX1BauJKOBGHUu+wJNKhc+l+Ocz9PuGw8RiCa5PI8/G+Yrx9Va0zYFTjV92rm6e51MtmVEjgtBOihtQa62WRUtEcBfoVKm0aSuabf6XOlGIRuPfJWtB/8NgeT4hX/KmFGNzBvSLWTZKV/m+MtZ1GiaXLAY8JqBaHXzQjtopH1dZDsJQCKPpj0TSZSt+UdztOScMFqoditN++zh9NRcSf5smL8hTYYrq0VIKfnmIye5oKuBqzfnviGubdILrFSyW1ExROJJvv+k/sUfDdbGZzk8GyUcT+WsEX/FxmYSKcmBgZmqXX/UyCUtpH2rz20/57t87TEi9i/fkF/TbjWsUbXWtPPVHde3Vcb8myty/ESqPW1eBnA5bOnviy/LPD5ffzyAw2bD5bCVHWkfS+g+bdPCgZai5/j/AfDM1f/KyQs/QKR+IwB1wVMAJFMy4/drBs2VE12fbXK35mjH6jrT28PGWNZoYyAlyRzlQ02lqlxtO5emPYge0NqlhPbZwVw8xUQoxrqW0j6wVp9p+qusFvFQ7wTHRkO8+zmb8goIs2Tkcmnbd+wbZHYuwWsu7SalScaC1R9Hb44gNBNTcOypYjDaSAGFVUxWg3l8b/80rT5XOuo2Fy3Gb6TWmna+Dl9A2hRdbnR0LEfO7xnZnmGJMf3ZoAcnWrEGoiU1qafHlRmIZmra7YbQbh5/As/cKBFfD6PrXg5AxL9RH49hHs9+LzabwO8qr6f2XFYQpmduBIEk6u1E2uY/k5DmIiDrcJHAGZ/J49OuTNqXVfM33RknUu04bCI9X8uBEtpFYEaOrytR027wOnDYBJMV1iL+55GTtPndvHJXV95jTB/M0DLkaksp+dYjp9jZ3cArLtbHWAv/eu9oEIBtOYR2U52L7Wv8VfVrF1JYxWQ1mMf39k9x+YbmRS0PTV4nNlF7oZ2v7jjo2Rd1LnvZmrY16jqjGlpGYZX5W7LLYcNpyd+2ljJNahptfldZJV+TKY1wNEkwmkwHoXWduhOAkQ2vTucdz2va/elzo9l+bU95nb6iydzpVVFLEBro0fKjUl9wu+bGcvq0TaFdTtpXNJHKsCaY4zkabaWzcWEltlqihHYRmF+CUs3jQgha/a6K+7T39E/zkp1rMlbp2ZjmnOXQtB/um6B3LMS7rt2U1viHZ6o/juPjIZrrnOlYgmx2b2xhb/90VUz1UurWhI6GwlbkuqZdG/P4H/3vHj5399GaXMtkPBjj1GRk0SA00LW2Fp+rKsGai2H6q3Np2mC25yxvUWUKbXsihCs+Q8LmZpzGBXXHrVhN5NamIZqmx8iUY7WLJrW0G6K93o0tGaFj8DcADG+4OX1cxL8BAG9oXmjn8muX49NeoGmnC6tk9hnXhXYTAO7oGMkcv90Wn4s6l70sTTs7fc0U2gfCTUs2Vak2SmgXwemp8jRt0E3klfRpRxMpgtHkgo5J2XTUuxEChpZBaP/PI6do87t5xa4u1jaZJVWrr/H3juYOQjO5YmMzwWiSY4ZGXkmmwnESKVlQjjYYKV/xZNW7j40Govz60Ch3Hxyp6nWyMXPiF/Nnm7T4XBW3Ri3FYpo2GPXHK+TTNn3D065uQGT4tN1Ztb+tJvIFTUN85dUfjyVS6Rztdr+bjsHf4khGmGm9lLn6Denj4p52ko46XPEZHLEZIIem7XaUlaedHRBqBqFF6zKF9lgwxih63I4nMppT057P1S5d07a6qhzxWZzxWZL2OnrDnpLSfSuJEtpFMDA1R53Lvmhv5KWoZHceIONHtxhOu42OejcjNS6wcnIizL1Hx3jLVetxO+w0ep14nfaqm8ellHpR/zX5e56bVbmqYSIvprAKQJ1bb7owV8VodoDfHRkF4MREmNkaNK8x2ds/hcthY2d345LHtvrcK8qnDZXp9BVPC2098nnUpZt+rT7t7IIqZttWsGja8VlAr7A4FY6XXBM9mtTS5vW2enfaNG7VsgEQgjn/en3shra9UNMur6d2ttDOLmFqMhGKMW1vBQxNO4dPG/QI8oEKadpm0GDAoy+yVrymLYRYJ4S4TwhxRAhxSAjxZ8b2FiHEb4UQvcb/Zss5HxNC9AkhnhVC3GTZfrkQ4hlj3xdFIWG1Kwgz3aucYbf5K6tp5yuMkIvOZcjVvvfoGAC3XLkO0FfBeqGX6i4eJkJxZucSbG3PL7R7mr2017vZX4UI8nQJ00J92jXq9PXbw6Npf9z+MzNVvZaVvf3TXNzduKgLx6TFX3vz+OxcArtNLIjYNqlET21TvPiCJwAYtPdgFyKjiMqimraZp53QLUNNdU5Smiy5n4GpaXuddlpSE7SMPYZmczK67qULjp0PRjOEdg6fdjnR49mLVbOwStS3UNOOuDoAcM+N5oweB1hXtqZtTffSU93GHHoU+9qzQNNOAn8ppTwfuBp4vxDiAuCjwD1Sym3APcZzjH23ABcCLwG+LIQwv5VfAW4Fthl/L6nge6k6A1ORskzjoJu0KunTTmvaBQjttctQFa1vLEiLz5Vhvu9qqv44eseMILRFNG0hBJtafVXprT0yq38uxZjHYWH95UoSiiV5tG+S11/WgxCw//RM1a5lJZpIcXAwwOVL+LNNKv0bKYTAXJIGjyPvgrzBU7me2qamfVqsxZ91zcU07bR5PBEAqdHkNdPjSpuraELXtNvr3XSd/gVCakx0vYCku2nBsaZf2zTtZ9cfbyjDp50rR3sxn3aizhTaYznztEHXtIPRZNrtUSzWoFD/7DEATolu7DZR0L22miwptKWUw1LKp4zHQeAI0A28GrjNOOw24Gbj8auB70spY1LKk0AfcKUQogtokFI+JnXH3bct56x4pJQMTEdKTvcyafW7icRTGS3oyqEYod3Z6GF4Zq7qflMrfUbf2YxxNHirHoh2fJF0LyvdzV4GqyC0RwNRhCjsc4F5jaqaaV8PHRsnntK4+dJutrb7OVAjTbt/MkI8pXHh2qVN46D7tGfnElWr0Z8LvcNXbtM4mIFolXEn+Ayh3aetTVcTM8nWtL1OO0bNE6TNSdLhQ0gNezKSjjov1a8dS6aYCMZp97no6r8TgOGNN+c81owg9+Yxj/vdpUeP587RHkMTdmLe+RztZEpjKhJH1usZKO653D5tKK9FZzSRImVZDNRPHwbgQGojnQ2edBGa5aKoqwshNgKXAk8Aa6SUw6ALdqDDOKwbsNapPGNs6zYeZ2/PdZ1bhRB7hBB7xsfHixli1ZgMx4nEUyVHjpuYxSMqpUmMZ5cgXIS1jV7C8VRFWgwWgulX3pIlONc2eRgLRklW8abcOxbC73Ysqen2NHsZCVR+LKOBKK0+d8Hlbs2bdzWF9m8Pj9JU5+SKjc1csq6J/QMzNVnA9U+GAdhQ4G/HjPavdhtbK/l6aZuYPu2y50tqaW312WRnRrqX02Fb0GNcCIHXmSsYbTY93lKbhoRjSSZCMXZ6RvDPHiPuamKi8/qcx6Y1bSNXO57UMuai3uNkLpEqaaG1IHI8MoRAEvN2Im3z730yHEdKsDetBXRNO1f0OMwHCw9MFb8gz3ZR1c/oQvuJSPeyli81KVhoCyH8wI+BD0opA4sdmmObXGT7wo1Sfk1KuVtKubu9vb3QIVaVctO9TCpdf3w8FKOlzlWQcOhM99WujYl8MhxnJpJYoO12NXrRpO6fqhZ9Rs3xpeIPupu8pDRZcXP9SCBaUKMQE9OnXa1c7WRK495nx7hhRwcOu41L1jcxFY6XdFMrFjPrYsMiRVWspCsH1tBEnq/Dl0m9x0E8pS0wCxeLe24Ue2qOmLuVkbg7Q2hna9kmOYPREuXVH48lUwzNRElqkp3iFAAz7VfkbTM5lzaP94OUSJlpIi9n0Zmdo236zbP92ebipK6pC4nAFZtASyRyLqTK0bSt78EZncIzN0rC4WN/uHlZy5eaFCS0hRBOdIH9XSnlT4zNo4bJG+P/mLH9DLDOcnoPMGRs78mx/aygEuleMF8VrVIpLWOBWMEm2FqmW0HutpgwX+ilmuPozdEoJBfdxo97sMJFZ0ZmowX7s2E+H7daPu09/dPMRBK88AK9DOOuniYA9g1Utz0p6L+dBo+DprrCsi5alqH+eGBuCU3bEOjlBqP5AnoQWqR+E8Foknr3/DWz/dkmuYLRnLEZXHYbLrutpAIrs3OJtLDfpJ0CINS4Pe/xcU8bSYcPZyKAM65/Z6x+6HLqj2dr2s1jT+pjbLk4Y7u5yG9t9BP3tCKkhis6kdOv3eh14nc7SopXsWrappY96duOhm3Z072gsOhxAfw3cERK+a+WXT8D3mE8fgfwU8v2W4QQbiHEJvSAsycNE3pQCHG18Zpvt5yz4jE17VKbhZhU3DweKlxom8FgQzXStE2hnV2RrCu9eKjOOGYjCcaDsZyV0LLpqVCd4mxGiyhhChZNu0rm8d8dHsVlt3H9dt1ydV5nPR6njQMDs1W5npX+yQjrC9Syobrd8PIxO5ekwZu/sct805DyPp86I3I8VL+JSDxVtKYdNXy8nsgwKanHTIwFShDaxm8EoCum+9hDjTvynyCEJRhtoV+7HKGdne7VMvYoAFNrrsvYPh6M4bLbaPQ6iXmMYLToWM4I8nJyta0LZ9Offdq9FWDZ072gME37OuBtwA1CiP3G38uAzwAvEkL0Ai8yniOlPATcDhwG7gbeL6U0Z+F9wDfQg9OOA3dV8s1Uk4GpOdr8brx5UkIKJV1/vEJpXxPBWMF1cNc0eHDYBIMV6DVbCH1jIXwu+wI/ULWrovWN65HjhWja5tgqGYwWS6aYjiRKEtrVSPmSUvLbI6Ncs6U1bcZ02G1c1N3I/hpp2htaFraLzUe6/ngNC6zk66VtMt+es0xN2whCm/bqAtCfUcI0973FY0mTM03GnsggiZRGd5O3JCvR7FxiPh4m1AcsrmkDROozK6NFMzTt0ufHKrSdsSnqpw+TsrmYabs847jxUIy2ehdCCGLe+bSvRDJ/BHkpi3Gri6p+5hAAR9mIXYiCKxxWkyV7BkopHya3PxrgxjznfAr4VI7te4CdxQxwpXB6KsL6MiPHAbwuOz6XvSKatpSS8WCM9gK/SHaboKvJU5UUp1z0GUFo2X7lBo+DOlf1Cqz0jpoafv5qaCYep52OendZOZ3ZmJpPKebxSLzy5vFjoyH6JyO897mbM7Zfsq6J2x7rJ57UcOXR8solpUnOTEe46cLOpQ82MOuP10rTjiZSxJPaoj7thgo1DTHTvcZcugexEE3bKszNNpXe8CCJpMbaJk/RTW80TRKIJhgPxdhQl8A7N0zK7k5r0vmYy87Vtmja/jIWndYc7eaxxxFIZtouR3Nk/n7GgzE6/Pq2mFd387jnxkjkydXuafby+IlJpJQF19ZYGDl+BIB98fWsaXAve+Q4qIpoBaOne5VnGjdp9bsr4tMOzCWJp7SiOs70NJW2+iyFXOleUP0CK31jITxOW9pfvRTdzaVpK/kYKbKwCuiar9thq6h5fCwQ5TN3HeX1X3kUl8PGiwx/tskl65qJJzWOjiwWV1oew7NzJFKy4CA00OuPN9fVrsDKUnXHweLTLjPtyxTaw3azGlr+EqbW7abMMRtoeMKDxFMaa5u8jMxGi6qfH4gm0DRdCF7uHQYg3LAVbItbEU1NO1eudqnm8ewc7daRRwCYWnNtxnFSygxX4Lymnb8q2rqWOkKx4nK1rYsOeyKEL3QKzebkiVD7srbjtKKEdgEkUhpDM6X30c6mtUIVn8ZDunAoJtm/3Jq8hRKMJhgJRBeke5l0NXqrVgf92dEgWzv8BXfiKdWMlg+zk1oxmjaUX7/ZRErJZ+46ynM+ex9fe/A41+9o544/vnaBuX7XOj0SuZo9xU+XmHXR6ncxVaPo8aXqjlv3lVNgRSTn8EaG0IQj3Ze5EPO4zSbSlpCoT0938kQGSWmStU0ekposKu3LfL/joRgXOfQiJkuZxsGa9pVD0zaFdpHf34wcbSlpGc3tz56dSxBPahahbWraS+dqF5MhYV00+2f0pjrB+m0MhzXWroB0L1BCuyCGZ6JosvzIcZNWX3ndeUzGiiisYtLTXMdoILagOEKlyReEZtLV6KlKHXQpJQcHZ7mwq7BCHqCnfQ3PzlWs29fxsRA2UXiKk4nP7aiIeXx4NspXHzjO83e0c+9fPp//fPNlOQubdDd5afO72VdNoT1ZmtDWO33VxqdtRoQ3eAoIRCvDp+2e1bXsOf96AsbLZPbSzn87NgV6tM4U2iOQStLZYGY/FL4Qn4kk0o2GtqGX6CxMaG8EjFxtKTN82g0l+rStkePe0Gm8kUHiriaCTednHJddrtmqaecT2uvSQaaFz4110dxgRI4PerYiZfnpvpVCCe0COF2hyHGTNn9luhiZq+uOIjVtqH4Eeb50L5OuJi9jwVjFq14NzUaZjiTY2d1Q8Dk9zV4SKb2VZiXoHQuxodWXV3PKR53LXhFNe59RnvT9L9jKxrb8AWBCiHSRlWrRPxXBYROsLdK0qHfDq62mvZh53Ou047CJsgLRPLN65Hi4fhPBaALBvC/Ybhc4Fqm1YAp0ze4m5unAJpO45sbS8SyDRfyerele65OngCUixw0S7haSTj+ORAhnbCpj4e922HDYRNFV0aw52qaWPd1xzQJTfXblx3T0+NwoM3nM32b1ypNGcZ9CyIwc1/3ZB5IbcNgE53cVfk+pJkpoF4BZx7pYzSkfrX5XWd15TOY7fBVutimn6EAx9I2HcNlteVenXY0eZBUKrBwa1FOYLiygm5RJOle7Qiby3rFQQelm2fjdjor4tPednsblsBV0k7l0fRMnxqvX8ev0VISeZm/BrgoT8zdSC8w0rsXM40KIsntqu2ePAxCp30wwmsTndqQroHmWaKSSEYxmmMi9kUHajUj7oQJjMiLxJPGkZtw7JB1RIwWtAE1bT/vaCOjatqbNdy4z56dY94FV024dNf3Z1yw4zox0bzNy+M3ypu7oGLORRM57ab3HSXeTl2dHCm+9mxk5rmva9852sqOzvuhFeLVQQrsAHjw2zsbWuqK1hXy0+twkjQjOchgP6XmLi+WXZtPTUp285GyOj4XY2FaXV3tIF1ipcFGTg0MBbALO7yx8VbwuvZApfyzxpMapifCijUry4XM7CFfAPL5/YIaLuhsLigi/ZF0TUL0iK6cnI6xvLTzdy6TF52ImUpv644Vo2mCWMi3HPG4K7U0EY8kMf7bbufhnlZH2ZUSQe8KDuJx63nKhC06rP7uTKdzJIHFXE3FPYZUnl/JrF2spSqd7aSmaxx4HYDLLnw364r6lzpW+nyTcTWg2J874LDI+lzfYbEdnfcFC2xo5LlJxfIE+JIJHQl3pYkQrASW0l2AunuLR45O84LyOpQ8uELPASrl+7fGgHk1ZTKvQNfVuHDZRdU17qYpk6UIvFQ5GOzQ4y5Z2f1H59OZirBIR5P2TYZKaLCjdLBuf2162ph1PajwzOMulhjBeisvWN+Oy23js+GRZ181H/2S4pFRJs5RpLeqPmxHhi6V8gdmes/TPJ20eb9DN44X6swE8FqFudr7yhvVc7bVN3oI17RnDotI/GeESt16QMtS4Awq8hyzW7aveXfyixhTaDTOHcCYCRHzriPrXLThuIruIlLBlFFjJ50rZvqae4+OhghZ/1gWHP9CHTUsw4eohgoddPYVb7qqNEtpL8NiJCWJJjRsqKLTb/JUpZToejNFWZJs4h91W9VztaCLFwFSErYsILrMqWqWD0Q4OzbKzCNM46GUiW32uiixkegvsLpYLn8tBpEyhfXQkQCypcen6wtpgel12LtvQxMN9E2VdNxezkQSBaLKowiom6QIrNTCRz84l8DrtS1om9PacJWraUuKezSphWkDkuInbmUPTjgySSEq6mzwFLzhnIglSmuSZwVmub9ArT4cbtxX8NsINWwDwzfbqY0lkatpFm8eN81vSqV4LtWzQNe3s1FZrgZV8i7vzOutJpCQnJ5b2a1sXzGYltINyE91N3nQTm5WAEtpLcO/RMepcdq7c1FKx12ytUJnG8WCsqCA0k2rnap+cCKPJxQVXg0evDVzJgLixYJTRQIwL1xYfMNLd7K3InBwbDSIEbMmRn74UvgqkfJlBaJeubyr4nOu2tHF4OFBxAdk/pd8oiylhamL+RmqR9jUciNJWv3Rd9LJ82sER7Imw3klLq2dkNprRfKJ+kch1yNS001XRwoMktMKroiVSGuFYkr6xEJF4iotdhad7mYSMqO56Ix3KqmkX21PbmqPdMvYYsDA/G0hHumdnyVjTvgJziZyd+nZ06orD0QJM5Nb7cbqz11z3itKyQQntRZFSct/Rca7b2oZ7iUCRYqhU0xDTPF4s1c7VTkeOLyG4Ohs9jFTQPH5oSC8SUqymDfqcVCIQrXcsxLrmupLK3frcdsLxVFntH/ednqaj3l1UC8HrtrUhJRU3kZeaow3z5vGJGmjah4cCBcVAlKVpT+qaaaR+k94SFbhsw7w1pHmJZipux3xf7bR5PDJvHg9Gk0vGyJh+3/1nZnDYBOuKiBw3idRvJGV3440M4ogHMgqjFFtnwMzRFqk4jZP7AJjuuGrBcfmyZKxpX1LCVA5te0u7H4dN8OwSBYRCsWTGAtEU2ge1jezK4WoqNrCykiihvQjHRkMMzsxV1DQO0FznRIjyfNoJoyF8MdXQTKqdq907FkII2Ny+uFm00lXRDhtC+4JSNG1DWym3X3LfaIjtJQShga5ppzRZVvvH/QMzXLq+qag4h4u7G6l3O3jkeGVN5P0l5miDpdNXleuPh2JJTk6EC1ro1Xucpfu0J0yhvZmnjIVVj6Fp+9yOglrrmsFo87nawyQSqXRMxlK1/Aen9e/3gYEZLlzjpT6oB8aFijCPS5uDUIN+vH/maEZxlPoiFzVm5Hj9zBHsqRjh+s0k3AvdOmaOdrYr0Jr2BTAdXnhtl8PG5nbfksFoZj0B0BcRZmGVfucWNmUFUgpRuZodpaCE9iLc96zu83nBjsoKbYfdZpRpLP2GNGU0hC9F0+6ucq72cUPbXMpP19XoqWgg2sHBWTa01i0ZUJSLnuY6YkmtpDaHJsmUxomJ0KK+/MUopycx6N+JU5ORgv3ZJg67jas2t/BIhf3apycjtPld6WYoxdBU58Imqu/TPjKsL/QKcak0eHVNsqQiPJN6U47pug0cHQ5y2frm9MKq2VfY99WMMNccXmLuVmxaAkIjlt9z/gXwdDjOeDDGSCDKWDDGCzqC2LQEc74eUs7iFplWE3mGpm1Ejxe68DVztE0te6btspzHmda4BZp23Xz9ccj/XdnR2cCzo/mFdjypMRKYn7umyadwJCP0yh7W9axLp+WZmKlky4US2otw79ExLuhqoLMK5etafa6ymoZkFxsohmrnavcVmKfc1ehlIhRL53qWy8GhWXbmqPxVCOaPsBwTef9UhERKlpSjDfN9k0vtqW127Co0ctzKdVvb6J+MpFvQVgK9yU5pGondqD9ebfP4QSOvv1BNGyi6gAiQ1rQPxTpISclllpiDJm9hfcbdObp92WZPp7+7Z/IIbSklxwyhZRbSuapuBCjOn20SbDoP0DXk7PaciVRhlqJESuOM8V0zhfZs6yU5j+0dC9HZ6MnoKw5kRI+DvtjNbvMJsGONn4Gpubym+zPTEaw9R1qHHwDgntQlXJzDn93odeJ12WmsK145qARKaOdhNpJgb/80LzivsPzFYmn1rwShXflgtGA0wYmJENvWLK1trm3SC6yMBsrXtmcjCQam5riwiEpoVnpayp+TdHexEs3jfqNvsrXAQzHsOz2D3Sa4qITAmeu2tgHwaAVN5KenImwoIUfbpMVX/frjh4YCtPldBQV0NpRTynTn65g47638bnoNzXXOjEp1TQXe/DMLrOhCW06fJpHScNpFXk17eDaaDhA7MDDL+pY6OospqpKFWWLUP3OEZEqmLQ9mCttSwWjxpMbe/un0cU0TptBeqGlrmqRvLMSOHPcTayCaSa4o8h1GvMKxHNq2pskFv/m24QcBeEhekrMcsvl5FdtboFIooZ2HB3vHSWmy4v5sk1a/u6ye2vPV0IoX2p0NHuxVytW+6+AIiZRc0FEq5zjMvtoVMJEfGjI0pnI17TJytXuNm0IpkeMw31O7VPP4vtMznNdZv0AjKYRtHX7a69080leZYLR4UmNodq4s318tqqIdGgpwwdrGgmIAzIppJQntS95E7xX/wH2jHi5d34zNuF6dy15wpa2MCHJLgZWjw0Ha/e6cQjulSY6P64vJYDTB8fEQu3oa8c8eA0oT2iFD0/YHehFaIm15KKSndiyZYm//dPocd2QYz9wICWcD4YbNC44/MzPHXCKV00VpDUTDMMnn+r6cZ0SQ5/JrjwajGZY+T3gQf6CXkPQy17k7Z0CpWYSno8FdaHp7RVFCOw/3HR2juc7JJeuK8w8WSlu55vFQ6Zq2w26jq9FTsbKdVu54apCNrXUZ5r98mAFb+ytQjcuMHC8l3Qv0G06j11nWQqZ3LERPs7ckHy7Mm8dLSfvSND3AqJhULytCCK7b0sqjxyfKDsYD3eQoJWwoR2j7ylvYLkUsmaJ3NMjOAr8z6aYhJaZ97Ts9TSKVaRovxsRqFe7ptK/IoDE2Zzprw0r/ZDjtd376zCwSuKLLQbORYhVoubjYt0HK6SfiW4dNS1AXOMFTA9OMzEaX7Kk9O6dbL62L0nnT+C4QC8WRqR2/6MI1CwLRUk4/SYcPeyqKI6H//nMJ7e4mLz6XPafQtgagAbSaWra2k5fuWthf3OOcX2S5HXaafYW5NiqJEto5iMST3PfsGM/b3l610P5Wvzvdbq4UxoMx6j2Okuvh9lQoL9nK0Mwcj5+c5OZLuwvSXLoavZzXWc89R8bKvvbBoVm6Gj1lFUHobiov7avUmuMm5k2vlE5fx8dDBGPJshaZ125tYyIUXzRop1D6zXSvMur1t/iqq2kfGwmR1GTODmi5KLWTlcljJybxux0Z1fKWSvWykss87g3rQrvF52J4JsqTJ6fYPzDD02dmODg4m47gBz3Vq7nOye7gvTiSEabbr2DOv76k95IORps9Siqld9YzU66s5nFNkwzPzvHkySl+f3KKSFa8xmKmcdB/Ux31brqbvGzr8C/QbDO0bSCW0Pj9qSn29k+zf0Cfg7FgjO1r6jP6xgeiCQ4PBRaY8psG7wfgeNO1bMrRbCe71O1ymMiV0M7BF+7pZTqS4C1XL1xpVQqzeESpZRpLzdE2qXQPaYCf7h9CSnjNpd0Fn3Pj+R3s6Z8uqlF9Lg4OzhZ8881HOQsZ0wxZiC8/Hz7Dp12Kpl1KUZVsTL/2w73l+7XNgLZyNG2z/niuohmVIO1SKTAOwhTapaR9xZIp9pya5tJ1TRmKQKH+bACPI4d5PKKXIm31u5ieizMdjjMRjDEWiDEyG037mxMpjUNDAXZ1N9Jz4gcAnNl8S9Hvw8T0a5udsGA+QO+ZwVn29k/x+IlJHuwd59BgIF0qNpvGyacAmG27ZME+M4DuciOf3ed2LOj/MO/Xnl/4z0YS6XkYmY1ycHCWBq+Tw0MB+ifCPHFikidPTC1wJ4hUPF3kZc2lL8853uzPq6PeXfOcbSW0szg6EuC/HzrJG3f3cMXGylVBy8YssDJRYopRqdXQTHqavYwGoxXL1ZZScse+M1y+obmo4KMbzltDSpM8eGy85GuHY0lOTISLaseZi+7m0nO1B6YixJNaSeVLTXyu0n3ajx6foNHrXJBTWgzdTV42tfkqkvrVPxnB47SVtbBsSy9sq9OB7ODQLPVuR8Etd83GPKVo2g/3TjCXSGUUVHE7bUXFHzjsNux2XUCkc7XDgyA1WnwupISZPErAnv5p4kmNl7aO0DBzmLirifGeFxf9PkzMCHL/7NH0Nq9hCRidjTId1n3dyVT+35ItGaV++ghS2Jht2bVg/2ggRjCa5Jotreltm9t96TmAzFKmi9HV6CEQTbLHEgCXjfPMo7hllFOOzTR3bsx5THYnOIfdli5LXSuU0LagaZK/vuMg9R4HH33p+UufUAbmDalUv/Z4KEZ7femmmZ7mOqRcuiBDoRwaCnBsNFSUlg16l6kWn4t7j5ZuIt/bP42UpQehmfQ01xGJp0oSEmbN8XLM474SzePT4Ti/OjjCq3atXZBTWizP39HOI8cn8978CyGlSR4/McnGVl9RRV6yMeuPl1PPYDH0ILSGgufMdF+U4tP+9mP9NHqd6aAoKDzVy4pZYCXl9BF3N2PX4riiE+kKcvlKI9//7BidDR6um/05AMMbX4NmL13YpDXtmSPpIDBTaM/lSLvKRcP0M9hkklDj9py54seMlsimBQh0P7LVemNq2v7ZZxe9ViFprnOH7wYgvP4FOffbbSKdQWBlTYMS2svGD/cOsKd/mo+97Px0RaZqYfpeS70hjecooF8MlU77umPfIE674OUXdRV1nt0meP72du57dqy0ohXA1x86QZvfzXO2tS198CKYLTqPDi9e8jAXZsBMOZq2y2HD67QXHcH+w70DxJMab62AO+f1l/cQT2rcuW+w5Nf474dPcGgowPuev6WsscxXRau8XzulSY4MB4pyqTjsNnwue9Ga9sHBWR44Ns4rd63NqHxWjGncJFcEuTc8mLbc5YoBOD0Z4fh4mJu21tE58EsABjf/QdHXthKr6yLhbMAVm8Zl5El7XPrYCrUUNab92Zfm3N87GqLJ62Rzlm95XUtd2rc91v0iANae/BH2xMJAPJOlskOGZubYFtBN45ENN+Y8psHryLkILSVboxyU0DaYDMX4p7uOcuXGFl5/WU/Vr5duzxks/oYUiScJxRYW0C+GShZYSaY0fnZgiBfs6CgpmvKG8zuYiSTYd7r4KPKnz8zwUO8E73nOprKb1F+3tY1Gr5NvPXqq6HP7xkKsbfSk015K5UUXrOEXB4bSJR6XQtMk33n8NFduakk3RyiHC9c2srO7gdv3nCnp/BPjIf7lN8d40QVreNWutWWNpa1CjXVycWI8RDShFZ1toJcyLU5of/n+PurdDl66szNjeym/lVzBaJ7IYLqqWi6hff+xMVx2G691PmoEoF1JJEd6lRW300ZHg5stHX4uWd+0sAOaEJnaNuCw2Vjb5MkZxZ6LJrMSWh6hfWw0mLMkr9NuS6diBVovYbr9CpyJIN3Hf5D3WvPZIQuF9mggyu2/eZDNYpi4s4FAniIvjSVYRqqBEtoGn/zlEULRJJ98zc6yTYyFUO924LLbSkppKaewisl8rnb5mvYjxycZD8Z47WXFmcZNnrutHYdNcE8JJvIv33eceo+Dt15dWhSsFZ/bwduv2cBvj4ymc1sLpXcsyNYygtBM3nLVegLRJD9/eqig4x/sHef0VIS3VTBo8g92r+PwcCBdLaxQUprkwz96Go/Tzqdu3lmWaRzmNe1yG+vkotTmMk11zqLqCvSNhbjr4Ahvv3ZDRiqg02FLm9uLwZOrRWd4ELfDjt/tWLDAicSTPH5yiqs2NrO5/4cAnNmyeABaW72b67a0cXFPE5vafLT53TnvNaEcwWgXdzdxbDS09KJTSku610KhPRmKMRmOZ/izrVijuE/t+EMA1vf+DyKVf4GXqxvaaCDKP//6WV4gHwdgqvM6pC3355IdOb5cKKEN/PbwKHfsG+SPX7CV7RW48RaCEIJWv4uxwPIIbYfdRmeDpyKa9m2PnqKpzskLSixE0+h1csXGFu4rUmj3jQX59eER3nntxrI1XJO3X7MRp93GNx46UfA5vaNBjgwHubiE7mLZXLmphW0dfr77eH9Bx3/n8X7a/G5uurBz6YML5FWXdON22PjB7weKOu9bj55iT/80f/vKC+ioQCpMU50LUaX64wcHZ3E7bGxZoqlNNtdvb+ex45MFj+mrDxzH7bDxrus2ZWxvLdH9ltmiU7dkmGlfrf6F/QwePT5JPKnxlvY+6meOEHc1MdadPwCtq8nDrp7GBYpLrqDXXMFou3oaSUmZjszPhzd0Gldsiri7JWfamRkjYvVnW7H2F5jsej6hxu145sboPP3zvNfsbvYyMBXhP+/v43dHRjk0NMvnf/MsF2lH+KDtdgBG170s7/mluDOqwaoX2jOROB+/4xnO72rgAy/YWtNrX7a+mXuPjuWsl7sY5VRDs1KJXO3fn5ri3qNj/NH1W8pqX3rj+R0cHQkWtYj4yv0n8DjsC26I5dBe7+b1l/fw46cG0/O8GFJK/uEXh/G57Lz7OeWPQwjBW65az4EzszxzZv7GJ6XkH35+mI/f8Uw6SGxgKsI9R8d405XrFpovy6DR6+SlOzu5c/9gwd/NkxNh/vnXR7nhvI6igxHzYdYfr4Z5/NBQgPO6GnAU0F3Lymsv6yapSX5+YGlLyODMHHfuG+SWK9YviDA23WPF4nFYzeO6G88bPg1SLshrl1Jy/7PjvKnpCC965i8AOLP1LUh77muvb63jwjzV4Vp8Lhz2zO3zNcjnhfbmdj8+l50DZ/ILbVsywvpj/wMYpvEc1zs2GqTOZee8PC1TM6K4hUhr2xuPfh1k7hTBG3Z0cM3mVs5MzfH93w/wb7/rpSM1ytdc/4ZdJhjY+hbGe27KeW6d215QJ7ZasDJGsYz83c8OMR2O8/k3XFzRG18hvPXqDczOJQq6AVgxO2OVo2mDHi09UIamLaXkc3cfpaPezTuv3VjWWEwtvVBte2Aqwp37B7nlynUVDxp873M3k0hp3FaAb/veo2M81DvBn71we8XG8drLe/A67Xz3iXlt+ysPHOebj5zk/544zQv/9UHuPjjM9548jQDedGX5roFs3rh7HcFokrsPjix5bCyZ4k++9xQep51Pv+aiss3iVlqrUGBFGppgKdXzzuts4IKuBn7y1NI+/68/qFtr3nv9Qv9xqd+VTJ+2/rm3jj7Klb97HS/XHiAQCqdTFo+OBNkduo9PRj+NPRXjzOY3cuKCD+R83Va/a1EroxBiwf0m3LAVTTioC57CG9S/q3abYGd3I88MzqJlB5ZqKdae+BHX3vVi1h3/PwBG178i5/V6x0JckpXTbqXe7Uj3F9df5+VE67rwBU/QNnRfznPa692867pN/NNrL+Kzr72IP722g9vr/w1vYprJNc/h2CWfyPv+qx2YXAy1DXtbYfzm0Ah37h/iz27cVnZhjlK4erNuCv3fx/t5w+51Sx4vpeS7T5zms3cfZX1LXdlfpIt7GvnxU2e465lhXlpk1DfA/c+O8/tT03zy5p05a/QWw+Y2Hxtb6/jvh0/ySN8kE6EYM3MJ3nTlet593cYMQSCl5D/v68MmdAFbaTa1+bjpgk7+9/F+3vf8LXnLksaTGp/85RE2t/t4+zWV8yk3eJy8atdafrp/iI+//Hz2nprmn3/9LK/atZZbr9/MR378NP/vO0/hsAleeP6aBQUnKsHVm1tZ1+Ll9j0D3LyE5vyZu45ycDDA19++u+Id8drr3fSOhdA0WbFYkwNnZglEkyWnCL72sm4++csj9I2F8mYL9I2F+N6Tp7n50u4FbRwbvM6SrVJui2IRbtjC8Qv/jHW9t9EwfZA/5CCvs/uZ+fl6RmQTG2IO3u98CBuSUzveS9/FH8qp1QIFuTM66j0ZKaLS7mKm7XJaxp/gml+/jIGtb+Hk+X/Mrp4mnjg5yfDQaS50jdIw/QwNU8/QOLkfr1EMJtC8k95dH2G646oF1/nN4RGGZ6O8YxFFwGYT+N3OdNEWaXPSv/1d7Nj/abYf+Cc6zvwaIVMIJIHmCxlf+0Lm6vXfqEjF2Rb6PS8e+AqNoeOEGrby9DVfyOvLBl3BWSmsOqEtpeT4eJj7jo7x1QeOc35XA++vsVncRAjB267ZwN/89BAHBmbYtUhLxdm5BB/7ydP86pkRrt/ezr+8YVfZlXjefNV6frT3DB+/4xl2b2wpSnPXNMln7z7KhtY6/uCKpRccSyGE4E1Xrue/HjzB8fEQbX43DR4H//iLw5wYD/H3r7oQh91GOJbk43c8w0/3D/HOazdWRWAB3Pq8zdx9aITvPN7PHz0vd+rStx87xcmJMP/zzisqbjp769Ub+MGeAf71N8f48VNnuKCrgc++7mK8Ljt3vv86vv7QCb758EluzaHFVQKbTfDGy9fxL789Rv9kOG/BnN8eHuV/HjnFO6/dWFCTmGJ5w+4e/vwHB/jpgUFec2n5WR0nJ8L84W176GzwlDzeV12ylk//6gh37DvDX9103oL9c/EU7//uU/jcDj704h0L9pdqGgf9c3E5bHr5YyE4eeH76d/xHjpP/4L2w9+kPdJHc/QwaUeNgN6L/pL+8/8o72sKUZirrdXnwm4XpCwFUw5e/S9sffrzdPX/lA3HvsXakz/hMu8aPug+Td2jC91Lc3Xd9F30F4yuf3nOWuO/eHqIO/cPcdWmFt6zhLupwevIqLQ2tOkNbDr8ZepCp6kLnU5v7zz9C7Yf+Cyhhq1E/BtoGXscRzIMQNzdzP7n/BcpV34rQ7PPVVLQYLUQlWgOUE12794t9+zZU5HX+pffPMvPDgyl6/Ge39XAF265pGbBZ7kIRhNc9el7eNlFXXz+DQurAmma5OdPD/G5u59lNBDlr27awXufu7liWkfvaJCX/8fDXL+tna+//fKCTZs/3T/In31/P1+45RJefUllfJjZaJrkc79+lq8+cJzn72jngy/czod+eIAT4yH+4kXb+ePnb61qpP+bvvY4j52Y5LzOem6+tJtX7lpLg8dBJJ5iIhTjlq89zuUbmvnWu66syvVf/aWHOXBmlhafi5994Lqar/aHZ+d43ufup7HOyV+8aDtv3L0uY6E4NDPHy774EN1NXn7yx9eWFdOQD02TvPo/H2EyFOPeDz2/rLS+4dk5Xv+Vx5hLpLj9j65ma0fpv/t3/s+T9I6GeOjDL1jwHfyrHx7gR0+d4bZ3Xcn12+db++4fmGEiGOOKjS1l9WJ+8uRUzrKgmqYxNdhLt22aRm0K99w4ocbtTK+5ZtHXa/Y5uXxDYdUfDw7OMpIjer5++jBbD3yWVqMMKEAAP7J5I8HmnQRadhJouYhww7acGq2Ukjv2D/KrZ0a4ZnMrn3rNTjYv0S1vaGaOw0OZNRXqAn00TexH2mxI4cCWitM89hhtw/fjTMzX1A827mBi7Q0Mbno9Uf/iSsfF6xrpKKOQVSkIIfZKKXfn3LeahPZf/GA/05E4N5y/hhvO61hgtlouPnHHM/xo7xke/9iN6dxNKSUP9U7w2buPcmgowPldDXzqNTu5bH3lu45946ETfPKXR/jn119ckJk+HEvysi8+RJ3LwS//5DlVT5H73pOn+es7D5LSJG1+F1+45dK8UaWVJBBNcMdTg9y5fzBd29uKwya4+4PXl1VQZTF+8fQQH/nR0/z3O6/g6s25U1+qzb7T03zyl0fY2z/N9jV+3nHtRsaDMfrGQjxl1Iz/xZ8+N2dzhUrx+IlJbvna4/zVTTtKtopNheO88b8eY2Q2yvfee3VJPcetmIvW77336oy0pB/tPcOHfniAP7lhK3+ZpWXvH5hhdi7B9dvayvL7Hx8PcXI8XPL52ezorC+4hepYIMrT+YLMpMQ/cwQB/GLAyf/un+Wzr71oySY+0USK7/9+gIf7Jrh+Wxvve/4Wrt7cuuQchWJJHj9eWCtZkYrTPPF7POFhptZcTdRXmNXG67Jz7Zalx1JplNA2kFLWfPIL4ehIgJf8+0N8/GXn8d7nbubR45N88Z5enjg5RU+zlw+9eEdFSlTmQ9Mkt3z9cY4MBbjrg89dVKMbmIrw3m/v4dhokNvefSXP3dae99hK8nDvBD87MMhfvGhHxf2mhXBqIszvjoyiSYnP7cDvdnB+V0PVrTTxpFbzAMlspJTcfXCEz9x9lP7JCELomQdb2/2867pNGdpktXjvt/fwaN8E9//VC4oOwHx2JMgH/u8pTk9FuO3dV1ZkATQXT3HFp37HS3d28s9v2IWUkt+fmubt33yCS9Y18d0/vHqB+2r/wAwOI1Cr3GtXoj68yXO2tRVswTD7BCxVvXBkNspf//Qgb75yPTcskgp6bDTINx85yWQozssu6uI1l67lik2tBeVESym5/9h4hrm+0mxfU19Wp7pSWVFCWwjxEuALgB34hpTyM4sdX0mhvZJ541cf48x0hM5GD0+dnqGj3s37nr+FN1+1vipmx2wGpiK89AsP4bQLPvKS83jj7nULFglPnJjkfd99imRK40tvvqwmN2vFyiGe1Dg9FaGn2Vt29bliOTEe4sX/9iBvvGIdn37NRQWdYwZu/uMvDlPvcfDFWy7l2gpaaD78owP88ulhbrqwk4f6JtKd9375J8/JGdi1f2CGNQ1uuhrLt/Dt7Z9iOlx+I5UGr5MrNxXXGOnAwMyS6ZBSSj5x50E66t188IXbM7YHokmGZ+fYPzDDPUfGaPO7edd1G9m+Rtf4i6nsV6l5yIXdJnjOtrZlSfVaTGjX1LsuhLAD/wm8CDgD/F4I8TMp5eFajmMl8vZrN/CB/9uHEIJ/vHknb7i8p6Y3xnUtdfz4fdfy/915kI/+5Bm+//sBPnzTDlJSMjA1R+9YkP99rJ/1rXV84+27l/Q3Kc49XA5b1VwBS7G53c9br97Atx87hd/t4KLuRi5c20BPcx3hWJJgNEkwliAYTRIyHv/64Ch3Hxrhudva+Nc3XlJ2imQ2f3DFOm7fc4b7j41z3dY2nrutjRvP61jUHFyp1KGuRm9FhFUpc9LV5FlSaAsh2NXTxD1HR/nEnc8gECAgMJfIaIjzvO3t6Xud21l8sZsGj7NqQrurybNicrOt1FTTFkJcA/ydlPIm4/nHAKSU/5TvnNWiaUspOXBmlgu6GpbVHCql5M79g3zql0cz2oY67YIbzuvgn9+wK6MakUJRK6bDcd733b3s7Z8mUYBJ1GETfOimHdxawcDNbMaDMVp9roJe/8R4qGKL3ZQmebC3fNPwNVta86Y0LsaeU1PMLNENbzQQ5RdPD5PSJJohZ3xuB12NHroaPXQ3eWmqm1/EXNzTWHQlvdFANKMIUaUQQk97LGVuKnP9FaJpA92AtTbiGWBhot4qRAjBJYukfNVyHK+5tIcbz1/DQ8cmaPW7WN9SxxqjVrlCsVw0+1x8/9ZriCc1jo0GOTwUYDQQxe/RYwzqPU7qPQ7jz0mr31X1BWYxmupSAVnFYLcJOurdZbXWrXPbSxZKWzv87Dm1eIOfNQ2eJdO2THpavCWVvq1WPfAt7f5lE9hLUetR5brrL1gqCiFuBW4FWL++8tWeFEvT4HHy8ouLL7iiUFQbl8PGzu7GsgO6ak2uXszl0N3kLUtol5PG1FTnor3eXVCp36Vo9jnZXmL6ncdpn89brxBrGjxsrGI2RLnU2g57BrDmFPUAC2p4Sim/JqXcLaXc3d6ugp0UCsXZT6UzV5rqXNSVWImwqc7J2qbysjC2dvjzFVgrGI/Tzs7uhQ1KiqGhgtq2z+3gghLK29aSWgvt3wPbhBCbhBAu4BbgZzUeg0KhUJwTdBVRa8LpsLGupY6rt7Sye2MLda7yNH+f21FWRUKbTS9cUm52TE+zl/oKWDEcdsGudY0r3g1YU/O4lDIphPgA8Gv0lK9vSikP1XIMCoVCca7Q3eTFaRckUpJkSiOpSYQAmxDYBLjsdnyG77oa2Sib233MRBK4nTa8Tjtepx2HXWC3CexC/++w23Aa2wCSKUkipWG3iYq01G3zu2nzuwlEEwzNzDETSaBpEk2CJiUOm8DpsOG023AYY5ASJBKHzYbbacPjtNPgcZS9kKkFNR+hlPJXwK9qfV2FQqE413A5bMvazMLtsGdUhCvsnOqMpcHjpKHz3M9sWXlJaAqFQqFQKHKihLZCoVAoFGcJSmgrFAqFQnGWoIS2QqFQKBRnCUpoKxQKhUJxlqCEtkKhUCgUZwlKaCsUCoVCcZaghLZCoVAoFGcJSmgrFAqFQnGWUNN+2qUghBgH+st4iTZgokLDORdR87M4an4WR83P4qj5WRw1P7nZIKXM2S1rxQvtchFC7MnXTFyh5mcp1PwsjpqfxVHzszhqfopHmccVCoVCoThLUEJboVAoFIqzhNUgtL+23ANY4aj5WRw1P4uj5mdx1PwsjpqfIjnnfdoKhUKhUJwrrAZNW6FQKBSKcwIltBUKhUKhOEs4K4W2EOKbQogxIcRBy7ZLhBCPCyH2CyH2CCGuNLa/SAixVwjxjPH/Bss5lxvb+4QQXxRCiOV4P5WmmPmx7F8vhAgJIT5k2XbOzU+xcyOEuFgI8ZgQ4pAxFx5j+zk3N1D0b8sphLjNmIcjQoiPWc5ZTfOzy/iOPCOE+LkQosGy72PGHDwrhLjJsn3Vz89qvDdXBCnlWfcHXA9cBhy0bPsN8FLj8cuA+43HlwJrjcc7gUHLOU8C1wACuMs8/2z/K2Z+LPt/DPwQ+NC5PD9FfnccwNPALuN5K2A/V+emhPl5M/B943EdcArYuArn5/fA84zH7wb+0Xh8AXAAcAObgOOr9PuTb35W3b25En9npaYtpXwQmMreDJgr3EZgyDh2n5RyyNh+CPAIIdxCiC6gQUr5mNS/Jd8Gbq764GtAMfMDIIS4GTiBPj/mtnNyfoqcmxcDT0spDxjnTkopU+fq3EDR8yMBnxDCAXiBOBBYhfOzA3jQePxb4HXG41ejL2piUsqTQB9wpZoffX5W4725EjiWewAV5IPAr4UQn0c3+1+b45jXAfuklDEhRDdwxrLvDNBd9VEuHx8kx/wIIXzAR4AXAR+yHL+a5ueD5P7ubAekEOLXQDv6DfhzrK65gfzz8yN0wTSMrmn/uZRySgixm9U1PweBVwE/Bd4ArDO2dwOPW44z5yGBmp9sVvO9uSjOSk07D+9Dv2msA/4c+G/rTiHEhcBngT8yN+V4jXM5/y3f/Pw98G9SylDW8atpfvLNjQN4DvAW4/9rhBA3srrmBvLPz5VACliLbv79SyHEZlbf/LwbeL8QYi9Qj25xgPzzoObHgro3F8e5JLTfAfzEePxD9BsKAEKIHuAO4O1SyuPG5jNAj+X8Hiwm43OQfPNzFfA5IcQpdI3q40KID7C65iff3JwBHpBSTkgpI8Cv0P11q2luIP/8vBm4W0qZkFKOAY8Appa9auZHSnlUSvliKeXlwPfQfdegz4NVqzTnQc2Pgbo3F8+5JLSHgOcZj28AegGEEE3AL4GPSSkfMQ+WUg4DQSHE1UZk4tvRzTfnKjnnR0r5XCnlRinlRuDfgU9LKb+0yuYn59wAvwYuFkLUGX7b5wGHV9ncQP75OQ3cIHR8wNXA0dU2P0KIDuO/Dfhr4KvGrp8Btxh+2k3ANuBJNT/6/Kh7c4ksdyRcKX/oq7Vh5n1D70E3X+5Fj9Z8ArjcOPavgTCw3/LXYezbje5vOQ58CaNC3Nn+V8z8ZJ33d2RGj59z81Ps3ABvRQ+SOQh87lyem2LnB/Cja96HgMPAX63S+fkz4Jjx9xnrewU+YczBs1gioNX8rM57cyX+VBlThUKhUCjOEs4l87hCoVAoFOc0SmgrFAqFQnGWoIS2QqFQKBRnCUpoKxQKhUJxlqCEtkKhUCgUZwlKaCsUCoVCcZaghLZCoVAoFGcJ/z9z3elmakTjFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(lynx[\"time\"], data)\n",
    "t_future = lynx[\"time\"][80:]\n",
    "hpd_low, hpd_high = hpdi(forecast_marginal)\n",
    "plt.plot(t_future, y_pred, lw=2)\n",
    "plt.fill_between(t_future, hpd_low, hpd_high, alpha=0.3)\n",
    "plt.title(\"Forecasting lynx dataset with SGT model (90% HPDI)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can observe, the model has been able to learn both the first and second order seasonality effects, i.e. a cyclical pattern with a periodicity of around 10, as well as spikes that can be seen once every 40 or so years. Moreover, we not only have point estimates for the forecast but can also use the uncertainty estimates from the model to bound our forecasts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Acknowledgements\n",
    "\n",
    "We would like to thank Slawek Smyl for many helpful resources and suggestions. Fast inference would not have been possible without the support of JAX and the XLA teams, so we would like to thank them for providing such a great open-source platform for us to build on, and for their responsiveness in dealing with our feature requests and bug reports."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] `Rlgt: Bayesian Exponential Smoothing Models with Trend Modifications`,<br>&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Slawek Smyl, Christoph Bergmeir, Erwin Wibowo, To Wang Ng, Trustees of Columbia University"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit 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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}