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": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "from IPython.display import set_matplotlib_formats\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "import jax.numpy as jnp\n",
    "from jax import lax, random, vmap\n",
    "from jax.nn import softmax\n",
    "\n",
    "import numpyro; numpyro.set_host_device_count(4)\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.diagnostics import autocorrelation, hpdi\n",
    "from numpyro import handlers\n",
    "from numpyro.infer import MCMC, NUTS, Predictive\n",
    "\n",
    "if \"NUMPYRO_SPHINXBUILD\" in os.environ:\n",
    "    set_matplotlib_formats('svg')\n",
    "\n",
    "assert numpyro.__version__.startswith('0.2.4')"
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAewAAAD4CAYAAAAuGtVZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eXRk133f+bm1r9jRK5psdrPZ3MRNrcWirMikRUmOLcqxnZEnYzGOEk5sjRMfe04iz2RGx3KccU5yxmONjzUjR3SoeJEd2Y5kWRuHoqxdZHNrbt3sbrK7gQbQqAaqgELty50/3rtVBaCW9+q9V4Vm3c85OAAuHlCvCq/e7/62709IKdFoNBqNRrO78Q37BDQajUaj0fRGG2yNRqPRaK4BtMHWaDQajeYaQBtsjUaj0WiuAbTB1mg0Go3mGiAw7BPoxszMjDx8+PCwT0Oj0Wg0moHx9NNPX5VSzm5f39UG+/Dhw5w8eXLYp6HRaDQazcAQQlxst65D4hqNRqPRXANog63RaDQazTWANtgajUaj0VwDaIOt0Wg0Gs01gDbYGo1Go9FcA/Q02EKI40KI51o+NoQQvyqEmBJCPCaEOGt+njSPF0KITwohzgkhTgkh7mn5Ww+Zx58VQjzk5RPTaDQajeaNRE+DLaU8I6W8S0p5F/BmIA/8NfAx4HEp5THgcfN7gPcDx8yPh4FPAQghpoCPA28D3gp8XBl5jUaj0Wg03bEbEr8fOC+lvAg8CDxqrj8KfND8+kHgs9LgB8CEEGI/8F7gMSnlmpQyDTwGvM/xM9BoNJo3EH/z/CJrufKwT0OzC7FrsD8E/Jn59V4p5RKA+XmPuX4QmG/5nQVzrdP6FoQQDwshTgohTqZSKZunp9FoNNcuG8UKv/Jnz/Kfv3dh2Kei2YVYNthCiBDwAeC/9jq0zZrssr51QcpPSylPSClPzM7uUGbTaDSaNyz5Ug2Alxc3hnwmmt2IHQ/7/cAzUsor5vdXzFA35ucVc30BONTye3PAYpd1jUaj0QDFimGwX1nSBluzEzsG++dphsMBvgioSu+HgC+0rH/YrBZ/O7Buhsy/BjwghJg0i80eMNc0Go1GAxSrhsG+nCmwnq8M+Ww0uw1LBlsIEQPeA/xVy/LvAO8RQpw1f/Y75vqXgdeAc8AfAr8MIKVcA34LeMr8+IS5ptFoNBqgWKk3vn5Ze9mabVia1iWlzAPT29ZWMarGtx8rgY92+DuPAI/YP02NRqN546NC4mAY7B85Ot3laM2ooZXONBqNZpdQaDXYuvBMsw1tsDWaAfKds1f5jb96YdinodmllEyDPZMI68IzzQ60wdZoBsg3z6zwZ09eol7f0dGo0TRy2PdcN8HZlSzlar3Hb2hGCW2wNZoBoqqAc+XqkM9EsxtROex7rp+kUpOcW9kc8hlpdhPaYGs0A0R5ULlSrceRmlFEGew3X2+MWdBhcU0r2mBrNANEFRVtlrSHrdlJ0QyB37J/jHDAp1u7NFvQBlujGSCqqCinDbamDYWycX3Egn5u3pfUleKaLWiDrdEMEO1ha7pRrNYIBXz4fIJbD4zxyvIGhrSFRqMNtkYzUFQOWxtsTTtKlTqRgHFbvnX/GJl8haX14pDPSrNb0AZboxkgKuSpQ+KadhQrNSJBP2DksUELqGiaaIOt0QyQRluXNtiaNhQrNaIhw2DfbBpsXSmuUWiDrdEMkGJZ5bB1W5dmJ8VKnUjAMNiJcIDD0zFdKa5poA22RjNAVNuO9rA17ShUakSCzdvyLfvHtMHWNNAGW6MZIIWyrhLXdKZYqRE2c9gAx/clubia3zLFSzO6aIOt0QwIKaXOYWu6UqzWG0VnAFPxEKA3eBoDbbA1mgFRqtZRLbX6BqxpR6lSa7R1AcRCAUBv8DQG2mBrNAOiVGlOXtIGW9OO1ipxgETY+FpfLxrQBlujGRiFljyk9pg07WitEgeIhw0PO1/WOWyNNtgazcAobjHY+gas2cn2KnFlsLWHrQGLBlsIMSGE+LwQ4rQQ4hUhxI8IIaaEEI8JIc6anyfNY4UQ4pNCiHNCiFNCiHta/s5D5vFnhRAPefWkNJrdiPKwYyG/vgFr2tKqdAZGLzboiIzGwKqH/XvAV6WUNwN3Aq8AHwMel1IeAx43vwd4P3DM/HgY+BSAEGIK+DjwNuCtwMeVkddoRgHlYc8kwuTK+gas2YqUklK1vqWtK64NtqaFngZbCDEGvAv4DICUsiylzAAPAo+ahz0KfND8+kHgs9LgB8CEEGI/8F7gMSnlmpQyDTwGvM/VZ6PR7GIKDYMd0jdgzQ5KpqjOlpB4SBWd6RSKxpqHfQRIAX8khHhWCPGfhBBxYK+UcgnA/LzHPP4gMN/y+wvmWqf1LQghHhZCnBRCnEylUrafkEazW1FV4jOJMJWapFTVN2FNExWBiWoPW9MBKwY7ANwDfEpKeTeQoxn+bodosya7rG9dkPLTUsoTUsoTs7OzFk5vd/G7j73KCwvrwz4NzS5EedjTiTCgC880W1GjV1tz2EG/j1DAp1MoGsCawV4AFqSUPzS//zyGAb9ihroxP6+0HH+o5ffngMUu628YipUav/f4Wf7ymYVhn4pmF1JsCYkDbBb1TVjTRG3oWkPiYBSeaQ9bAxYMtpRyGZgXQhw3l+4HXga+CKhK74eAL5hffxH4sFkt/nZg3QyZfw14QAgxaRabPWCuvWHImjfghXR+yGei2Y0UWorOQLfqaLaiNnStfdgA8bBfR2M0gBHutsKvAH8ihAgBrwG/iGHs/0II8RHgEvBz5rFfBn4COAfkzWORUq4JIX4LeMo87hNSyjVXnsUuYaNYAWB+rTDkM9HsRlTIc9r0sHWYU9NKw2AHtxnsUEBv7jSARYMtpXwOONHmR/e3OVYCH+3wdx4BHrFzgtcSGwXDYC+k80gpEaJd2l4zqhS1h63pQrscNhiFZzokrgGtdOYqG2ZIPFeukc5Xhnw2mt1GsVLDJ2AyZnrY+iasaUFNctuew9YGW6PQBttFlIcNOo+t2UmhbKhYJSK6VUezk1KHkHgi7CentcQ1aIPtKiqHDTqPrdlJsVojGvSTCCl9aH0T1jQpdMlh682dBrTBdpWNQvNNNa89bM02CuU6kaCfuDkyUd+ENa00c9g7Q+K63kED1qvENRbYKFYI+X1EQ34dEtfsoFg1JjEF/D7CAZ822JotdGrrUn3YupBVoz1sF9koVBiLBjg0FdUhcc0OiuXmJKZEOEBWG2xNC8rDjoa2GuxY2E9dNn+uGV20wXaRjWKVsUiQuYmY9rA1O1A5bNCVv5qdKA87HNipdAa6DVCjDbarbBQqJKNBDk1FWUgXMFrSNRqDQouHrQ22ZjvFao1wwLcj7B03ixTzWmhn5NEG20U2ihXGIgHmJmOUqnVSm6Vhn5JmF1Gs1FtC4n7tMWm20JoyaSWuPWyNiTbYLmLksA0PG3Rrl2YrxUqtUQFseNi6rUvTxNjQ7bwlJxojNvX1Mupog+0ijRz2ZAzQ4imarRQrzRy2nsCk2Y7RRdDOw9ZtgBoDbbBdRFWJz00aHvZCWnvYmiaFytYqcR3i1LTSuqFrRYfENQptsF2iWKlRqtYZiwSJhQLMJELMr2kPW9OkWKk3WnZ00ZlmO8VKnXAXg62vF4022C6hZmGPRYMAHJyMaQ9b00BKaXjYgZYcdrlGva47CTQGxZbroxUlZav1xDXaYLuE0hEfMwc7HJqManlSTYNS1ZSdDDWrxEHPxNY0KVZ0DlvTHW2wXUJN6lIe9txkjMVMgZr2oDTslJ2M68pfzTY6VYlrKVuNQhtsl1hXBjtiGOxDU1EqNcmVjeIwT0uzS2gOdmgWnYEuJNI06VQlDrpIUWOgDbZLbJg57PGocSNutnbpPLamOToxGjJz2CFdSKTZSqcqcdBFihoDbbBdYmO7hz2pxFN0HlvTLSSub8Iag1YlvO3EQn49P12jDbZbNIrOGlXiuhdb00R52KroLBnRIXHNVoqVGuE2OWwwQuJaS1xjyWALIS4IIV4QQjwnhDhprk0JIR4TQpw1P0+a60II8UkhxDkhxCkhxD0tf+ch8/izQoiHvHlKw2GjUCUU8DV2yOGAn71jYV0prgG6eNj6JqwB6nVJqVrfMQtboUPiGrDnYf+YlPIuKeUJ8/uPAY9LKY8Bj5vfA7wfOGZ+PAx8CgwDD3wceBvwVuDjysi/ETAGfwS3rB2ajOmQuAZoGuymcIrxWYc5NdDS9qeLzjRdcBISfxB41Pz6UeCDLeuflQY/ACaEEPuB9wKPSSnXpJRp4DHgfQ4ef1ehZElbmZuM6pC4BmitEjfeco0q8aK+CWtaIjAdQuLxsF+3AGosG2wJfF0I8bQQ4mFzba+UcgnA/LzHXD8IzLf87oK51ml9C0KIh4UQJ4UQJ1OplPVnMmTU4I9WDk3FWFovUKnVh3RWmt1CwVSpUlXA0aAfn9BFZxqDYnXr9bEdHRLXgHWDfa+U8h6McPdHhRDv6nKsaLMmu6xvXZDy01LKE1LKE7OzsxZPb/io0Zqt7BuPUJdwVc/FHnnUDVmFPIUQxEM6zKkx2N6nv514KECuXEVKLcQ0ylgy2FLKRfPzCvDXGDnoK2aoG/Pzinn4AnCo5dfngMUu628IjBz21pB40vS49c5Y0+6GrL0mjaJ3SDxAXTavI81o0tNgCyHiQoik+hp4AHgR+CKgKr0fAr5gfv1F4MNmtfjbgXUzZP414AEhxKRZbPaAufaGYKNQ3eFhJ808ZVbnKUeedjfkeNivq8Q1QLPtr920Lmhqz+uIzGgT6H0Ie4G/FkKo4/9USvlVIcRTwF8IIT4CXAJ+zjz+y8BPAOeAPPCLAFLKNSHEbwFPmcd9Qkq55tozGTLtqsT1HFuNolip4RMQ8jcNtlH5qwuJNDvb/rbTKrQzmwwP7Lw0u4ueBltK+RpwZ5v1VeD+NusS+GiHv/UI8Ij909zdFCs1ytX6jirxhFaz0pgUyoZOtLnxBSAR0SFxjUFpWxfBdvTmXwNa6cwVmqM1t4XEIzokrjEoVnfqRMdD2mBrDLb36W9Hb/41oA22K2wUjDfR9hy21ovWKArlnTrRWgxDo2h0EXQIicdCen66RhtsV2h62FtD4nFdKKIxMUYnbn27xbXB1pj0autK6PnpGrTBdoXGpK5tHnY44Cfk95HVN+WRp1jeOetYt3VpFEpYp1cOW18vo4022C6gZmFvz2GDLizSGLTLYSfCfio1SamqvaZRZ7uwznZ00ZkGtMF2haaHvbPoPhEOaL1oTaNKvJW4DnNqTFRIPBzo4GGrHLa+VkYabbBdoFOVOOg8pcagWNlZdKbDnBpFqWLUOLS2/bUS8PuIBH266GzE0QbbBbbPwm4lqQ22BqNtZ3t+MqHDnBoT4/poHw5XaO15jTbYLtBO5UyRiOg3mca4Ie/MYWsPW2NQrNQ7tnQp4uEAeX2tjDTaYLvAeptZ2AqjEljnnUadQhsPShcSQTpX5n/56xcaVdKjSqFNBGY7cS1lO/Jog+0CG4UuHnY4oJXONGYOu31IfJQ3dN87v8qf/vASpxYywz6VoWIlJJ4I+3U0ZsTRBtsFNoo7J3UpEmE/m6XKgM9Is5uQUlJoExJvCuuM7vWxli8DzdbIUaVYrXec1KWIhwO66GzE0QbbBbKFnbOwFYlwkGKlTrWm59iOKqWq2bLTIYc9ymHOdM4w2OuF0d20gKpxsBIS1wZ7lNEG2wU2ihXGO3nYER32HHUagx10W9cO1rTBBlRbV68qcR0SH3W0wXaIlJKNQveQOEB2hMOeo04nneig30co4Bvpm3A6rw022KkS1xv/UUYbbIeUqnXKtXqXojNjXXvYo0uhMTpx59stHvKTH+EKaeVhb4y4wbZSJZ4wc9hSygGdlWa3oQ22Q7rJkkIzJD7KhUWjjgqJt/OgIkF/4+ejiPawDSwJp4QD1GVzA6gZPbTBdkg3WVJohsRHubBo1FE32Eiog8Gujm5BYjpnvH+0wbZmsGG0+/ZHHW2wHbJeMCd1dcxhG+t6AMjo0s3DDgd8I+1h65C4QbG6U2t+O2rzr9Nro4tlgy2E8AshnhVCfMn8/gYhxA+FEGeFEH8uhAiZ62Hz+3Pmzw+3/I3fMNfPCCHe6/aTGQZND7uT0pnutR11GlXinTzsETXYxUqtEX0YZQ+7XpeUqzuFdbYTC+muglHHjof9L4FXWr7/98DvSimPAWngI+b6R4C0lPJG4HfN4xBC3Ap8CLgNeB/wB0KI7lvKa4BmDru9h51UHrbeFY8szSrxnW+3SNBHqTKaIXGVv/b7xK4x2OVqfeCqa6pPv7eHrQ32qGPJYAsh5oC/D/wn83sB3Ad83jzkUeCD5tcPmt9j/vx+8/gHgc9JKUtSyteBc8Bb3XgSw0QpNHXKYTc8bB0SH1mUTvb2PmwwbtKjWkSkwuGHJqO7xmB/5cUlPvD73+XcyubAHrNR49BhFrai0bev1c5GFqse9v8F/CtAuQLTQEZKqa6cBeCg+fVBYB7A/Pm6eXxjvc3vNBBCPCyEOCmEOJlKpWw8leGgPOxkh5C4mmOrQ+KjS7Fq3pDbGezA6IbEVcHZ4Zk4pWp9V7wOahPx5OtrA3vMRo2DxRy2jtaNLj0NthDiJ4EVKeXTrcttDpU9ftbtd5oLUn5aSnlCSnlidna21+kNnY1ihXCHWdiKRDio32QjjPKw2xrsoK9h0EcNpSN+eDoO7I7CM9UTf/LC7jPYWhlP094t3Mq9wAeEED8BRIAxDI97QggRML3oOWDRPH4BOAQsCCECwDiw1rKuaP2da5Z0rsxErH04XJHUM7FHmmaOsl0O29/IcY8aSkf88HQMMArP9oxFhnlK5M1w85MDNdjWctjaYGt6ethSyt+QUs5JKQ9jFI19Q0r5j4AngJ81D3sI+IL59RfN7zF//g1pSPN8EfiQWUV+A3AMeNK1ZzIkltaL7BuPdj0mHvazWRy+96AZDoVyDZ+AkL+TwR5RDztXRgi4zjTYG7vgPaJaphbSBZbWCwN5zGbKpEcOO6TnEow6Tvqw/zXwa0KIcxg56s+Y658Bps31XwM+BiClfAn4C+Bl4KvAR6WU1/yVt5gpcGC8u1eQCAf0m2yEUaIYRu3lVsIjXiU+Hg0yGQsBu6O1K99S0PXUhfRAHtNqSNzvE0SCPl10NsLYMthSym9KKX/S/Po1KeVbpZQ3Sil/TkpZMteL5vc3mj9/reX3f1tKeVRKeVxK+RV3n8rgkVKytF5kfw8POxEOkNVhrJGl3SxsRSTgp1yrU6uPnj70Wq7MVCzUmHS3Gwx2rlzj8HSMWMg/sDy2VYMNxr3E6/Ta5568xGe/f8HSsd88s8JiZjCRCI1WOnPERqFKvlzjwIQVD1sb7FGlWOmsYqXWSyNYeJbOl5mMtxjs/PANdr5UJRkJcs91kwP0sDvXOGwnEQ6Q9bhF9I9/eJFPf+u1nseVq3X+6aMnefT7Fzw9H00TbbAdcNncWR6Y6OFh66KzkabYZRJT1FwfxcKztVyFyVioITqkZH6HSa5cIxby85bDU5xe3hiI199NunY749Gg59X0i5kiC+kC2R41BRdXc1TrUmtMDBBtsB2gilL298hhx8MBfVGPMN0GO6j1USw8S+fKTMWDBP0+4iH/rgiJF8o14uEAbzk8iZTwzCXvvWy1WWsnXbudsWjQ0+K8YqXW6EV/9Uq267HnU4a4zKgK/wwDbbAdsLheBHp72MlwgHKtPpJhT42adawNditSStbMkDh4b4iskitXiYX83HXdBAGfGEge246HPRYNerqxac1Hn17uZbBzQFNnQOM92mA7YClTIOATzCTCXY9ragDrC3sUKXYrOhvRkHi+XKNcrTNlVoiPe2yIrJIv1YiHAsRCAW47OM5Trw/AwzY38mELOeyxSJAND1MHi5li4+vTS90NtpJv1R724NAG2wFGD3YEv6+diFuTxhxbHRYfSQqVzpOYwsrDHrHoiwq7tnrYu8Fg58rVRmj6LddP8txCxvPIWLFcQwhj1GovVA7bkLZwH+Vh7xuLcHp5o+uxjZD4ED3sVLbEa6nB6b4PG22wHXA5U+BAj5YuaOqM68Kz0aTULSQeGM2QuJrU1ephD1uaVEpJvlxrDOx5yw1TlKt1Xry87unjFqt1wgFf2z797YxFVXrNm4jM4noBIeDv3TTL6eVsx42BlJLzu8DD/j++8gq/8JlrXn/LMtpgO2BpvcD+Hi1d0OJha4M9knTPYRtvwVETT9nuYe+GkLjqh1dzp09cPwnAkx6HxbsVJW7H6571xUyBPckwb5obJ1usNup0tnNlo0TO9KyH6WFfWs1zOVNgJdv+PN9oaIPdJ/W6ZNmCaAroObajTvcc9mh62Bmz53pqFxnsvFljEjdD4tOJMEdm4jzrcaV4t+tjO8pgexWNWMwUOTAR5eZ9SQDOdAiLq3D4/vFIY2DKMFgyNxQvXe4evn+joA12n1zNlajUZE/RFGiGxLXa2WhS6NKHHRnxHLYKiY9FguTLNSq14UUalORnLNyciTSTDJPxeCNR6CKss52xiPce9oHxKDeZBvuVDoVnymDffnB8aJvNWl1yZcMw2C94nLbYLWiD3SeqmtJKDntYRWf1uuR3H3uV+bX8QB9X00RKSbFS71klXiiPVkg8nS/j94nGZnY8anweZh5beYqxln7oeMjvecg3W6w0onC9aHjYHrTASSlZXC9wYCLCWCTI3GS0Y2vXuZVNkuEAh6djQ8thr26WqJqSvl7XGXTi6maJJ06v9BSZcQttsPtkyaymtJLDHlZI/Fxqk997/CxfeXFpoI+raaKKg8K66GwLa7kyk7EgPrPDYjw2fD1x9f5UU7EAYqHAloEgXpAtVhsbl16MeZjDTucrFCv1hq7EzfuSXUPiR/YkiIYCFCo1z6rWu6HC4WORwNAM9g9eW+UX//NTXBqQU6QNdp80RFOseNih4YTEz5i747Xc8NtlRhVliHvmsEcsJJ7OlxtTusD7YiorFNp42NEBedhWDXYzh+3+vUS1dKm6nJv3jXE+lWvb1nZ+JcfR2TjRoB8p8axqvRtKafK+m/ewuF5kdbM08HN4dTmLT8DR2cRAHk8b7D5ZyhSIBH1MmJ5BN3w+QTzkH3hIXBnstJkv1AyeQo9JTKr3dtSEUwwPe3cZbFX1HA+3eth+8h5HPwwPu/d9BJr1MF68Tmo2wkHTwz6+L0mtLjm/ktty3GapyvJGkaOziYYW/jAqxZWH/Z5b9wHw4uLgC8/OXMlyeDpuuQbBKdpg98nSepED41FLvZNgDAAZdEhc5Z/W8tpgDwvlCXXyoHw+QSjgozRiIfF0rsJkvGmkdoPBVqHv7R6211XQmzZC4kp33Ytc/1JjmJGR5rtlv1F4tl1ARQmVHJ1NNFrgvN7UtGN5vUgo4OOdx2aA4eSxz17Z5Ka9yYE9njbYfXI5Y60HWxEfwBzb7Sjx/jXtYQ8NJRDS6k1uJxLwjV4OO19utHRBs/p5mEVnSjo41prDDgYoV+tUPaper9clm2XrHjZ4pwq3uF4kHPA1/i+Hp+OEAr5GpE6hJElv3BMnYm5uhuFhL64X2T8eYTwa5Prp2MANdrFS48JqrlFRPwi0we6TpXVrKmeKZDgw0Bx2rlRtFELokPjwyJgGu1vqJBL0j1RIXEpJeltI3MtiKqs0POxw08NW3rZXHuRmuYqURuGUVcY9GpRyOVPgwEQzahjw+zi2J8Er2wz2+dQmAZ/g+ul4ozZjGAZ7eb3QmJR4+8Hxgbd2nVvZpC7huPawdzeVWp2VbIn9PaZ0tTLokLjyrvePR1jVBntopLcJhLQjEvSPVNFZtlSlWpdbXpNI0E844GNjiHr7jbau4NaQOHhnkLLF7imTdoxFPPKwM4UduhI37xvj9NLWkPj5lRzXTccI+n2NDc0wWruWWoSr3nRwnIV0obFBHgRnV4x77E17B1NwBtpg98WVjSJSwoEec7BbiYcGOxNbGewfOTLNeqHiWUhP0x0rIfFo0D9SIXEV8dn+moxHg6znh1l0ViUc8BHwN2+LDQ/bM4NtPF+7IXEvqsSXMsUdUcNb9idZyZa2pNXOpzYbVdGq2GrQBrtuiqbsUx72gXEAXhyg4tmZ5U2CfsHhmfjAHlMb7D5Q1Yl2PexB5rBPL2eJBv3cMWdcyOkh3ghHmUy+Qjjga3hq7YgEfSMVEm+onMXbGOxhhsRLtS0V4tDMZ3vVi92Xhx0NuP46VWp1rmSLjR5sxXEzP/v8fAaAaq3OhdVcw2A3PGyPe9W301CabITEx4DBKp69eiXL0dkEQf/gzGjPRxJCRIQQTwohnhdCvCSE+E1z/QYhxA+FEGeFEH8uhAiZ62Hz+3Pmzw+3/K3fMNfPCCHe69WT8hrVr2jHw04OuOjszHKWm/YmmEkas7rTulJ8KGzP1bYjPGoedn7r4A/FsA12rlzd0S8f8zwkbt/D9iKHvbxuRg23hcTvPDTBbDLMv/jcs3z7bIr5dIFKTXJ01vAqo0PysJdNp2mfGRGYiIWYm4zy4uLgDLZxjx1c/hqsedgl4D4p5Z3AXcD7hBBvB/498LtSymNAGviIefxHgLSU8kbgd83jEELcCnwIuA14H/AHQojBNK+5jJIlteNhqyrxQSkCvXoly/F9yYZWs64UHw7pfLlnr76Rwx4lD9vM67cLiQ/dw25vsL0LifeXw84Wq9Tq7t1LGk7ItnvaWCTIX//yOzg4EeUf/9FT/IevnQbg6B7Dw27m+Ad7/TbuwS1O05sOjg+sUnyzVOVypjDQ/DVYMNjSQE0ID5ofErgP+Ly5/ijwQfPrB83vMX9+vzDKDh8EPielLEkpXwfOAW915VkMmKX1AmORgGX9XzBC4rW6HIgi0NXNElc3yxzfN9bwYrTBHg7pfKWnhx0ZsT7sRg47vnUj41W7klXyldqWli5oGiSvDPZGHwZb9ay7qV+t0nzbDTbA3GSMz//SO/ix47N8+YVloKns1Xx9BhsSX15XqmxNg337wXEuruYHcg2dvaIKznafh40Qwi+EeA5YAR4DzgMZKaX6Ly0AB82vDwLzAObP14Hp1vU2v9P6WA8LIU4KIU6mUin7z2gAqBF0dlDGPTuAwiQ2Ty8AACAASURBVLNXzTaM43uTTGuDPVTS2/qN2xEZsZD4Wr5M0C92bHi9aleySr5UbeNhG+dYqHiVwzae75jNojNwV570ciPN1/6+lggH+H9/4QS/ct+NfODOA41NQ3RI42GXNoqE/L4t763bDxr1Oi8NICyuinqPD7AHGywabCllTUp5FzCH4RXf0u4w83M76S/ZZX37Y31aSnlCSnlidnbWyukNnKWW/j+rqJvTIPLYSuHs+L4kEzokPlQy+YqFkPhoFZ2pvP52lcCxqPuhXjvkyjs97EGExIN+0ZCotYIXqnCLmQKTsWDX4ki/T/DrDxznkz9/d2Mt6PcR8ImBz8ReyhgV4q3X0O0HjMKzTiNB3eTM8iaRoI9DkzHPH6sVW+VtUsoM8E3g7cCEEEJd3XPAovn1AnAIwPz5OLDWut7md64pFjMFW/lrGOzErjPLWabiIWYSIUIBH8lIQBtsk9PLG42CFa+p1yWZfO+is1Hrw17LtY86eBHqtUO+XN0iSwotId+Sd0VnyUjQssQxNEVW3IxGLJqiKf0QDfmHUnS2b5vTNBU37ncrG96/v8+uZDm2J9mYNjcorFSJzwohJsyvo8CPA68ATwA/ax72EPAF8+svmt9j/vwb0qi0+iLwIbOK/AbgGPCkW09kUBTKNdL5SkMg3yqDDImfuZLl+N5k4yYwFQ/pKnHg6YtpHvz97/KJL700kMfLFqvUZXeVMxi9kHimQ15/2HriuVIbDzvovYdtJ38N3owi7SfNpxiGjsDSRmFHl44QgtlEmNQApnYNo0IcrHnY+4EnhBCngKeAx6SUXwL+NfBrQohzGDnqz5jHfwaYNtd/DfgYgJTyJeAvgJeBrwIflVJec3eppTbFDlZIRAbjYdfrslEhrpiMhUbew764muOfffYkpWp9ICEzsCaaAkpLvD6UmcLDYLuOuGLYBrtQrhLf5mEH/D5Cfh95z3LY9g22F7rri+s7DaBVYgMYkNJKvS5ND3vnBmMmEeLqprf3unSuzEq2xPF9g60QB+h5pUgpTwF3t1l/jTZV3lLKIvBzHf7WbwO/bf80dw8rWWP3tndsd+awF9IF8uXaFoM9HQ+xPIAw0W4lky/zi3/0FFJKPnDnAb50apFipeb5SLy1fPtq6O2EzfMoVesDG9M3LOp1yWKmwI+aE5ZaGfNwdGQv6nVpVIm36fyIhb2biZ0tVkiGrRecgfsbm41ihWyx2reHHQl6PzO8ldVcmUpNtnWaZhLhRsW7V7w6pApx0EpntkmZBnvWFCSxSiMk7rHBPtOmenEyProedqla4+HPPs1CusCnP3yC99y6l7o05BW9pjn4o3cOGwZfaTsMFteNDWW7m50XoV6rFKs1pGSHhw1GWHw3hcRjIT9+n3Ath72U6dzSZYVB57BVDUo7gz2bDHPV45D4q+a0skFXiIM22LZpGOyETYM9oJD4GXN2besNcdo02KMScm3lv55c4MkLa/yHn7uDtxyearwuZ694b7DTHQRCthMJGm/DUagUP2ve7I7t2RlOHPegXckqzdGaOw12NOSlh21vtCYYuVo3RWYWt83BtkvMw9enHYuNtGS7kHiY1VyZuoedBq8uZ0mGA+yzGWV1A22wbZLaLBH0i8bNxSrRoB+fwPMBIC8tbjA3Gd3S4zoZD1Gq1gfeerEbeC2VIxr084E7DwBww0ycgE80wlpeYj2HPToethKcuLGLwR6Gh90YrRlqExIPBTwTBtkoVmx72GD2rLu0sblipsvspvkU0eBwPOztVeJg5LBrdelpke2ZK1lu2pe0VdnvFtpg2ySVLTGTCNsu5xdCNORJvaJUrfHts1d5541b84OjLE+6tG6MDFRvrlDAx+GZOK8OwMPO5Cv4RG8Vq0ZIfARau85e2WQ2GW6bJogG/QT9YkgG23jttwungOFhe7HZrdclm6WqrVnYirGIewNA+k3zKQadw15aN0RTptsULqrZCV4Wnp1f2WwbIRoE2mDbJJUt9X1hez0A5HvnV9ksVXnvbfu2rKuK3FFs7WrXX3rT3kRjlq2XGDrioZ6bu1ELiXe62bkd6rVDdw/bGw8yV64ipb3BH4oxF1XhVrIlxqNBwoH+Ch69en06sbReYO94e6dpJqEMtjd57EK5xmquzKGpwQqmKLTBtkkqW7Kdv1bEwwFPRSG+/tIyiXCAd9w4vWVd6YmvjqCHfTlT3NEzf2xPkktrec+9AisqZzA8ecdBI6XkXA/vZCwSdLVdySoqh93Ow46F/J7UnvQz+EPhpu56KltiT59OCAw+JL60XmT/WPsCOa8Ntsqf95vvd4o22DZJbfbvYY9Fg5552LW65LGXr/Du47M7dsoqdJQessHO5Mv84z96kgtXcwN5vFK1xtXN0o7ilJv2JpEDqBRPW1A5g2Zb1xvdYC9vFNksVbmxSzvMsAaAdPOwo8GAJ5u7psG272G7mcN2ck8DiIYCA62PWV4vsr+DwVTOlArzu83ltGGwD05oD3vXU6tLVh1c3MlIwLMK2Gcupbm6Wd4RDgd2zcSur764zDfPpHjywtpAHm+5MYFo65tbjcTzuvBszcIsbBidkLiqzO/mYQ8rJN6tSjwW8pP3YDPVnIXdTw7biES40fmxki06M9hBP+VqfSAa8E3RlPYGeywaIOT3eZbDdlpR7xRtsG2wlitTl/QdPhqLeDeN6GsvLhPy+3j38Z0DU8YiAQI+MXSD/djLV4DBefpqAtH2kPjhmThBv/C88MyQ4OztPUUawilvbA+7W0uXwhgAsrs87FjYm6IzJyHx8WiQcq3ueFyvlNJ5SDxkmJFBhMXX8mXKtTr7O1S0CyFMtTOPQuKZAj7Rf0W9U7TBtoHTasqxaMATLXEpJV97eZl33DjdNrwmhGByyHriuVKVb5+7CgzO01/sIAgR9Pu4YSbeaDHyinS+3IhudGNUhFPOrRhDaaa71IAkI968R3rRrUo8Fgx44kFuNDzsforO3FGFy5aqFCt1xyFxYCCV4s2Wrs4iLzPJsGch8YVMgX1jEYL+4ZhObbBtoETl+zbYLoaxWnllKcv8WqFtOFwxFQux6rHGbje+fTZF2fQGBmewDQ+7Xfjs2N4kr3pYKV4o1yhV65aKziKB0QmJt+u/bkUZ7EGL/OTKNYRo9sS30hyx6e5GQm1M+mnraorMODPYyrDtSfbvMQ6yaFL1jHcKiYNReOalh92vIpwbaINtg6bKWX8XdzISpFqXroeOvvbSMkLAj9+yt+Mxw57Y9fWXrzAeDXJ8b3JgBntpvcBMItRWn/umPUnm1wqeCWJYFU2B0fCwpZRdW7oUYxF3Qr12yZeqxIL+tq1CasSm2x6kk6IzNQDEqYftNGoITYM9iMKzlcYGo/P5ehsS73+qmRtog20DdXHPJHvfhNuhwlhuF5597aVlTlw/2fVNNxUPDa2tq1qr843TK9x/8x72jIUbQzG85nKXN5cqPDu/4k3FetNgW89hv5E97NRmifVCxYLBdn/WsxVy5VojtLudpofttsGuEPCJRtGhHRoetsPXacUFg61en0HksFc2zHtwl7TKTCLM6qb78qT1umRpvcDBSW2wrwlS2RKJcKBtYYoV1K7YzaKa+bU8p5ezXcPhYEyMGlZb11MX0mTyFd5z696BjvpczBQ40CHXdcxsLfKqUjyTN/7HVjxsv08Q9Is3tNLZObPAr9eEo2TjPTLYPHa+XG2bvwYvDbYx+KMficsxl2RcUxY81l6oDecgctgr2SKTsSChQGfTNZMIU61L17sNUpslKjWpPexrBaf9imMu7YpbUQbnxOGprsdNxcNkCpWBtF5s5+svLxMK+HjXTbNMDWhymJTGGMdO/ZqHp2OE/D7P8tjqOVopOgMjdzqIkPhmqcq/+vzznhXldEJdpzfu7Z3DhmEY7FrHjbjyvN3PYVf6CoeDe4NSVrKGzKfd2QitND1s7/9nK9lSz3x7U57U3Wu82XUynApx0AbbFqlssW+VM2jejNwMiTc9ue5vuKlYECmbIx8HhZSGoMs7b5whHg4wFQ+RLVYbBWhesVGoki/XdrR0KQJ+H0dm455N7WqO1rR2IwwHB2Ow/+5Mir84ucATZ1Y8f6xWzq5sMh4N9nz/JD2IQlkhX662Ha0J3nvY/ZB0aXa4klp2MsiimeP3PqWzki2xZ6z7NdQQT3HZYC9mhiuaAtpg28KJjjg0Q+JuetjqDdtrhzxlXsSDLjx7ZSnLQrrAA7caBXFK19zrjcPlhsBB5/DVsb1Jz0LiaXMjNRG16GEHfQPJYT97KQ0YU8wGiSo462UYvKrz6EWuVCMW7uBhe1RU5cRgB/0+4iG/K1XiMw7uadD6+nj/P7tq4R48a9YYuS2eolTOhiWaAtpg22LFqcFWNyMXw33KYPcKramJXYNu7Xrs5SsIAfffstVge10At2jBYN+0J8FCuuCJTnQ6XyYRDnTNtbUSGZCH/ex8BvBelnU751Y2OdYjHA7D9bBjbboJwJgBAO6HfDcchMTBHRlXJ7MRFMrD9vr6bYq89AiJeyRPupgpkIwEHP3PnKINtkWKlRrZYtUdD9vFYoj1gjFP199jItRk3HjsQXvY3z13lTvmJhqv29SAdM0bIv1d+jVV4ZkXxsvq4A+F4WF7e8MrV+u8cHkdgNcGaLBXN0us5crcuKd7wRkML4dteNjXTkgcTD1xhxublIUQcy+Uh+11lXgmX6Fcq/cskBuPBgn6hQc57J2DhAZNT4MthDgkhHhCCPGKEOIlIcS/NNenhBCPCSHOmp8nzXUhhPikEOKcEOKUEOKelr/1kHn8WSHEQ949Lfdp9mA7q6YM+X2uh8StGIbpuHHeg27tWkjnOTobb3w/OA+7SNAvurZ/qDefUk9yE6uDPxRG0Zm3IfFXljYoV+scno5xaS1PtTaYNjIrkqSKRCiAEIP3sAuVGvGORWde9WFXSHYIw1thLOLMw67U6qzmyo497MiA+rAbPdg9NhhCCKbjYa667GFfzhR2v8EGqsCvSylvAd4OfFQIcSvwMeBxKeUx4HHze4D3A8fMj4eBT4Fh4IGPA28D3gp8XBn5awGnKmcKt+VJ1wsVSxWeDQ97gAa7WquzvLF1Vzqo2dyLmQL7x6NdZ1HPeJTrAiOHbbVCHMyQuMdtXc+Y+et/cM8clZpk3szJeU3DYFsIift8gkQ44GrayAq5UrWzh+2BQZJSslmqOg6JO8n1q/SYUw/b7xOEAz7PPeyVrLGxtrLBmEm6L54ybJUzsGCwpZRLUspnzK+zwCvAQeBB4FHzsEeBD5pfPwh8Vhr8AJgQQuwH3gs8JqVck1KmgceA97n6bDzEDUUgcH/er1WDHQ74SYQDrOUG57msZEvUJVvGW06Y5+p1Lt14c3XPdTWiDh6oIqVzZUuiKQojh+2tx/vspQz7xiLce+MMMLiw+EI6T8jvY5/FgQleDslpR9VUVuvkYQf8PkJ+n6sGO1euUZf9Df5QjEUDjjxsOwawF9GQn6LXHvaG8rB7X0eGPKl795jNUpX1QmX3G+xWhBCHgbuBHwJ7pZRLYBh1YI952EFgvuXXFsy1TuvbH+NhIcRJIcTJVCpl5/Q8xS2DnYwGXS86s9pDORkPspYbXP9tu1F0Ab+PiVhwIB52J9EURSjgYywS8ETG0HZIPOij5LGH8ux8mnuun2ikKAZVeGa3dWjQA0DU6Mx2ozUV0ZDf1SrorIPBHwqnOWy37mlg5LEHFhK3cL6zLuuJL6mWriGqnIENgy2ESAB/CfyqlHKj26Ft1mSX9a0LUn5aSnlCSnlidnbnqMhhkcqWEKIZ0u2XsUjAVQ87k68wbrF1aCoWYi0/OM9lcb39tKypmLcyqdVanSvZkqXd8EwyzFWXz6Vaq5MtVm0WnXlbJZ7KlphfK3D3oUkmYiGm4qGBtXbZbYc0DPbgrtN8YxZ2Z283FnLXIDkZrakYiwTJFqt9iyE1VM5cGBUZDfk9D4mnsiXiIX+jar8bM0nDYLs1RGZhF4imgEWDLYQIYhjrP5FS/pW5fMUMdWN+VkoMC8Chll+fAxa7rF8TpDZLTMVCjseqGW8yd25GUko2bHjYhsrY4D3s/dsqtafiIU9z6SvZErW6NQnBGQ+KUzIF67KkikjQR9FDMZnnzHauu6+bAODITHxgBntlw57BNtJGA/SwTc+5kzQpmAbJVYOtPGxnVeIAm31GI5THOpNw5oSA4WF7LU26ki1a3lzMJMJUau7Jk1ppEx0EVqrEBfAZ4BUp5f/Z8qMvAqrS+yHgCy3rHzarxd8OrJsh868BDwghJs1iswfMtWsCp6IpirGoewU1xUqdcq1uIyQeIj3AHPZSh77FSY/lSduF4jsxk3Tf27ercgbeS5M+cylNwCe4/eA4AEdm47x2dUAh8c2SLa3qZCRAtjRAD7ts1cN2bxOx4WBSl2LcoZ54KltiIhYk3GakqF1iA/Cw7ehgqE2IW2HxxUyBgE84GkPqBlbcxXuBXwDuE0I8Z378BPA7wHuEEGeB95jfA3wZeA04B/wh8MsAUso14LeAp8yPT5hr1wSuGWwXi86sqpwpjFD04DzsTn2L014b7A6h+HZMx92fnatUzuykT1RI3Ks50M9eSnPrgbFGC87R2QRXN8use5wiKVfrrOXKtm50STPUOyiUcE4naVIwjLkXIfF+ZmErnA4AWXEotdxKJDiYkLjVjV9DnjTrzn3mcrrAvvFIT70Lr+l5tUgpv0P7/DPA/W2Ol8BHO/ytR4BH7JzgbiGVLXFkJt77wB4kIwFK1Tqlas3xzjZTMC5GqwZ7OhGmWKkbqk59Thyzw9J6YUc4HExPP19GSulIw7gTnULx7ZhJhMnkK1RqdcfpDoUK99sNidclVGqSUMDd16Raq3NqYZ1/eKKZkToya44XvbrJPdd51115tY92SFV05tX1sR1liKNdDbbf1U2mG0VnU2arZr/jat1yQsAIiXs9UGZlo8iPHd/T+0DcHwAy7DnYCq10ZgEppeNJXQq1K3bDg1DekdXQ67QKE7m06+xFp77F6XiISk2S9UASVD3umEUJQfWauHkzTvcTElczsT3oxX71yib5cq2RvwYjJA7ea4rbqexVjEWD1OpyIPOVAXKNHPa1VXSm2hL7rUtZseGx9sLrkHiuVCVXrtkIibtrsC9nCsxpg31tsGFOl3IrJA7uyJPaDYl7NcWmHYVyjXS+fd+i8jzXPOrFtiNw4IXucNrGLGxFOOidHvOz84Zgyt2Hmp70dVMxAj7heS92yqI6VSteTLXrRjOH3aXoLBhwvejM7xNdH7MXU4n+5wMoXW7XPGyXNzTbsbvxm4gGCfjckSdVAlDaw75GcLNf0c0BIHYNttu7zm40tLzbFH6pG02/obyej20jfKWKU9wsPEvny4T8Pls344g5JKTkgXjK0xfTTMdDHJpqviZBv4/rpmID8LCNegK7OWwYnDxpvpHDHlzRWbZYJREOOAr5J8MBQn5fXwIh2VKVUrXuWhFVJOitcIrdjZ/PJ5hOhFyJJtrpOvEabbAt4IaOuCLpgYc9ZtVgJ92tnOzGUsa4Ue9vI14y5bWHvd5b5UzR2MS46GFncoa+u52bccQDD3uzVOV/+28v8lfPXOadx2Z2nM+R2YTnleIrG4Z+wbSN1qGGhz2gwrOc8rC7tHXFQv7GcW7gdPAHmJrZif5aNZVqmFsettch8X42fjMuiacs7hLRFLBQdKZphpCdau5CMyTuSg67UMEnsDxAoCnF6X0OuznsvY3BjnvnYedKVTIdQvHtUMUpblbP21U5g6bBduum98SZFf7Xv3qBpY0i/+TeG/if33vTjmOOzsb51qspanXpWfXrSta+fsHYoD3schW/TxDqco7RkJ9yte7aa5V1OFpTMRUP9fV+djNqCEbRWbUuKVfrlkfK2qEhS2rjfGcSYVfSf5d3iWgKaA/bEk0P2/k/rBkSd8fDHosGuw64aCUU8DEeDQ4sJC4E7G0jdNAw2B60dl1czQNw/ZS1iv54yE840F9YsRNXNoq2N3eRoPFWdENP/O9eTfGLf/QU8XCAv/yld/C//9StbbsCjszGKdfqLKTzjh+zE/3kSccG7WGXasRC/q4RkeaITXfOacMFDxuMzo9+lPqaHqtbOWw1M9wbL3slWyJkyhpbZSbhjijSQlp1nQzfw9YG2wIp82JRxtYJbhedWc1fK6YT7k+xacdipsBsItx2tx0L+QkFfJ6onV1cNXKy10/HLB0vhHAtdKboZ26umyHxJ19fJeAT/M2vvLNry5Zq7fIyj52yoU6lGHQOu1DuPFpT0TBILoXFs8Wqox5sxUyf6oVeeNjgTdEkmD3jNvToQU3sKjvWNnh5cYO5yaglSVSv0QbbAnaHF3QjFvLj9wlXPGxDR9yewTZ2nd6HxJfWi+zvYLSMebXe6IlfUB62RYMNRuGZWx52sVLj6qY1HfNWIgH3bnjnV3JcNx1rbAI6oXQFvBwCspIt2a79UJ7noMRTcuXOozUV8ZC7IzaHHhLfNJwQu/ePTkRDhinxqlK8n0jNnmSEcq1OxqE40HPzGe46NNH7wAGgDbYFUpulRq7TKUII16YR9eNhuz3FphPGsPfOntVkzBs98UtrOabjIVs3w5lE2LURm8umypp9D9sMibugJ34+tcnR2d6zp6fiISZiQV676o2HXa8brUN20wNqUzu4HHZvDzvmusF2LySeL9dsh+pTG+45IWC0vYF7EYjt2NWjB5gzi8RUDrqvx80WuZwpaIN9LXF+ZbPxz3cDt+RJ7Qz+UMwkQp73YUspWcoUu+Z8phMeedhX87a8a3UubmoOg/0hAW6FxKu1OhdWc5YMthCCIzNxzq9442FnChWqdWk7T6o2tYPqw86Vql1VzqA1R+v8nKSUbJbcMtj99WK76YRAUyXOjdenHXb16KFpsJ3UaDx3yRiaow32NcLVzZKxw5pz7x/m1gCQfjzsmUSYbLFKyQNFLUUmX6FQqXU1WpOxkCczsS+u5jg8bU9C1vCwnee6oHUMX38G2+lM7Pl0gUpNNmZe9+LobMKzkHg/rTiKQY7YNDzs7gbbTQ87X65Rq0tXQuLTfRZwrmy4p3IGzRx2oey+jkA/evQAc5PGxl0VjfXDc/OZLUNzho022D04tWDssO6Yc+8flgw797CllP0Z7KT3rV0N0ZQuWt5T8ZDrfdjFSo2ljSLX2faww1Tr7oziW8yY1fHjw6kSV97y0T29PWyAW/aPcXWz3DCubuKk1zcZHtwAECOH3aPoLOiewXZDllQxnbDfliilZHG94KrBjoXcbUts5WqfbbXj0SDJSMCRwX5+IcPN+5M960EGhTbYPXhufh2fwNUd1ljUeQ47V65RrUtbbQ4wGLWzxUzvaVlT8RDZkiH56hYL6TxS0oeHrQRlnG8gVHW83cEuboXElbd8dMaawb71wBhgVMK6TT864go33iNWKdjysJ2fkxuDPxTKw7Zz7V5YzZMtVl29p0WC7ra9teLkOpqbjPUdEq/XJafm13dNOBy0we7JqYUMx/YkXS3pH4sEHVeJ25UlVUy7PCe2HUumh72/S9GZ6sV2Myx+4ar9CnFwdxOzmCn2pYgU9Pvw+4Tj4R/nU5vMJMKMW9zI3bLfMNgveWCwnbQOJV14j1glV+o9vU793A0Pe6NhsN3LYdsJiT+ntOWvc88QqQ2NF21dKxv9p1bmJqN9e9jnU5tkS1XudDEd6hRtsLsgpeT5+Qx3HnI3fzEWdR4SV5O6+qkSB28ndl3OFAj5fczEO9+ovRBPubimDLb9HDa4kyawM3hkO5GAz3lIPJWznL8G4/qZm4zy8pIXHnaReMjf12bXrU6KXpSrdTZL1Z7yvo2iKhcM9stLWcB6FKQbsVCAaNBvq8vh2UsZ4iE/x/YkHT++opnD9sBg9zFARjE3GWV+Ld9Xfcqz80Y61M2NjVO0we7CQrpAOl/hDpd3WMlIwAhp1/q/OdvVEVfMDGBi11KmyL7xSFcFNk8M9mqOZCTApM00gVtRByml2c7Wp8EO+h15KFJKzq1sWs5fK247MMYrHoXE7YqmKNyIQlnhcqZAXRrTy7rhZtHZ0xfWmE2GtwxjccJ0wl4v9rOXMtwxN+GqHK3a0OQ98LBTWVOPPm5P7heMkHiuXOurF/u5+QzJSIAjLmys3EIb7C48N+9NSb9SO9t0MA96vWDOXI7au4ijIT/xkN/jHHbv4RteGOwLq3kOT8dt95ZOxkL4BI57sVdzZUrVetdiu24YBrv/Tdxarsx6oWKppauVW/eP8/pqjpzL88lTffTOKsYiATZLVep155X73VDKeId7pFGCfh8hv88Vg33yYpo3XzfpWg/0dDxkWZ60WKnxytKG615jOOBDCDyZ2LWSLTEdDxGwoUevaLZ22Q+LP3cpw51zE5alnweBNthdOLWQIRTwcXyfe6EjaHrFTvpMGzlsm94kGJXiXlaJL60XOdBDd9cLg31pNWe7QhzA7xNMxcOkHL4m/fZgK8JBn6McthJAsRMSB6PwTEo4vZzt+7Hb0U/vrCIZCSKlUcHtJUp73sp1Ew35KTg8nysbRRbSBU4c7iwZa5fpRNiyPOmLl9ep1iV3d5Gs7QchBNGgNzOxU9kis32OAe23F7tQrnHmSnZXFZyBNthdeX5+ndsOjNmaNGSF5nCD/kN+/RadgXtj59pRq0tLw94nzPN2y2BXanUW0oWenlInZlwQT3E6hi8S8Dvqw260dNn0sG9rVIqv9/3Y7VjZKPbtYQ9qxObF1TyxkN+SfKoxE9uZQXr6olHw9ebrXTTYNuRJn/VQCCQa9GbE5kq2/41fv73YLy6uU6tL7rzWDLYQ4hEhxIoQ4sWWtSkhxGNCiLPm50lzXQghPimEOCeEOCWEuKfldx4yjz8rhHjIm6fjHtVanRcur3tSIejGTOz1QgW/T/RsR2mHG8apEyvZIrW67FohDhAwJ++4ZbAXMwWqtp/KjwAAIABJREFUdWm74Ezhhjzp5Ux/sqSKSNBZ0dn51CbhgM+2h79/PMJELOhq4VmuVCVXrvVV2QuDGwByaS3HdVMxS+HpaMjvOEd78kKacMDHbQfcK2SdMnPYVgqrnp1PMzcZdW3oRytRD2ZiSyl5LZXreyPe7MW252HvNoUzhRXX8T8D79u29jHgcSnlMeBx83uA9wPHzI+HgU+BYeCBjwNvA94KfFwZ+d3KudQmhUrN9QpxaB2x2b/3kMlXmIgG+8qDGR62NyFxO2HhqVjItZnYzbGa/b2xp10YAHI5XSAW8vc9UMFp0dn5VI4bZuK2i4mEENy6f8zVXuyUg95ZaL5HvK4Uv7BqXco2FvI7roJ++uIadx6acHVm9Ew8TLlWt1QT89yljOvhcEU06Pz12c7iepHNUpWbHKQljV5sex72c/MZDk54s7FxQs+rRkr5LWBt2/KDwKPm148CH2xZ/6w0+AEwIYTYD7wXeExKuSalTAOPsXMTsKs4NW+EB92uEIeWEZsOQ+L9GobpRJh0vuyoSr0Tl8zWKitepptqZ43ioZnhediqpavfYqKIw5Di+ZT9CnHFrfvHOL2cde2acNKKA4PxsOt1yaW1vOWoTCwYcFSYVyjXeGlxgxMuhsPBup748nqRxfUid3vkNcY88LBfNesqju91YrDt92I/N5/hrl3UzqXod5u3V0q5BGB+3mOuHwTmW45bMNc6re9anlswSvpv6DPE2o1m0Zkzg223pUsxmwghpbsFX4onX0+TDAcaoxu7MRl3T0/8wmqeSNDXt0c3nQiRK9cceQiL6/33YIMKiff3+MVKjfm1vO38teK2g2OUqnXXJncpqVOnOWwvPewr2SLlat2yh+005PvcfIZqXbpacAbNAs5e8qRKMMUrQxTxoOhMFUIec2ywrfdiz6/luZwpeLaxcYLbRWftXAvZZX3nHxDiYSHESSHEyVQq5erJ2eHUQoY75sY9KelPhJ3fjPqZ1KVoKnu5b7C/d/4qbzsyZakFYzru3hzqi6t5rp+y39KlcEPtbNFBDzYYRWf95rAvruapS/sV4opb9xupH7fC4kpHvP8ctpk2ckHfvRMNZbwpix62w6Kzpy8agcp7XA5JWxX+efZShpDf1ygydJtoyFlKpx2vXslyYDziaG636sVOW+zF/m/PXgbgvbft6/sxvaJfg33FDHVjfl4x1xeAQy3HzQGLXdZ3IKX8tJTyhJTyxOzsbJ+n54xipcbppaxnknR+nyAZDjgKiWecGOykN3riC+k8F1fzvOPojKXjj8zGubpZcuU8Lq7mbEuStjLjUDylWKlxdbPcdQZ4L8JBf99T1Boa4n162Edm44QCPtcKz1KbJYJ+YVvERtFMG3nnYV9aM6IJtjxsRwY7zY17EkzE7AuAdKMREu8RMXt2PsOtB8Zs69xbxY0c/3bOLGcd5a/BXmuXlJK/fGaBtx+Z4lCf9TBe0q/B/iKgKr0fAr7Qsv5hs1r87cC6GTL/GvCAEGLSLDZ7wFzblby0uEG1Lj3JXysMeVJnfdh2B38ovBoA8v3zqwC848ZpS8erDdHzpkBNv9Trkotr+b7z1+BcnnRpvffAk144qRJXLV1H+vSwg34fN+9LuuphzybCjvL5Ib/P05D4xdU8Qb9gv0Whm3go0Pdwi3pd8vTFtOv5a2gJiXd5P1drdU4tZDyV2XQ7JF6t1TmX2nSUvwZ74ilPX0xzYTXPz775UM9jh4GVtq4/A74PHBdCLAghPgL8DvAeIcRZ4D3m9wBfBl4DzgF/CPwygJRyDfgt4Cnz4xPm2q7k6y8t4/cJ7rneu4s7Genfw67XpcOQuDcDQL53fpXpeMjyG+xNc+P4hHODrXKRveQluzHtcBNzOe1MNAWMKtt+Q4rnU5scnIj2HGLRjVv3j/HS4rorc8FXsv33YCu8nol9cTXP3GTMsoKWk5D4udQmG8Wqq/3XinDATzIc6JpeOr2cpVipe9qm5OT6bcfFtTzlap2bHBts1Yvd28P+/NMLxEJ+3n/77guHA/R8d0spf77Dj+5vc6wEPtrh7zwCPGLr7IZArlTlT5+8xPtu39d3/s0KY5Fg3zejzXKVuuxPNAWMHHo44HM1hy2l5Hvnr/IjR6cte1WxUICb9iZ5bsGZYIfKRdodq9nKdNxaWLETDdEURx62n2pdUq3Vbcswnk/l+vauFbceGONzT82zvFFkfw+lul6ksqXGjbJfjE2thx622YNtlWjIT6lap1aXtlvnTl4wCr5OHJ6y9XtWmU6EuhaRqkEWbufPW3G7SrxRIe4wJD4eDTJmYS52oVzjb08t8f7b97s6ndFNtNLZNj7/9ALZYpWPvPMGTx9nLBroOySuJnX1WyUuhDB6sbPuedjnUzmubJQs568Vdx2a4Pn5jCOvTrV0OclhR4KGl5Lq8zW5nCkgBOzrU0fcOAfj7Vi0OSO8WqsbLV195q8Vt+53bzZ2Klvqu6VLMRbtf1PbCymlUaho45pRA0D6MUonL64xHQ/1LQDSi6l4qGuV+HfPXmXfWKQRHvYCpXTmRoQG4MyVLELAjX22KrZipRf76y8vky1V+dk3zzl+PK/QBruFWl3yyHdf5+7rJjzdiYKzaURKlnTCQeXkTCLk6sSu75+/CsC9FvPXirsOTbBeqHBhtb8h82CEzoJ+4SgcDabGugMPe28y4kjGNhLsb6bwM5cy5Ms13naDM+/tlv1jrqQoKrU6a/ly3y12Ci9HbKbzFbLFqi1lvGhjJra9c6rXJT84v8qJw+4N/NjOdKLzfIBStca3z6a475Y9nj0+GK+PlFCyueHsxJnlLIen4433hRNUa1c3Pv/0AnOTUcfvIy/RBruFx1+5wsXVvOfeNTibie1ER1wx0+UN3g/fPbfKwYmo7Tyy0up1YiSevpDm6GzC8bjAmUSo76iD0YPtLIUSCfRnsB8/fYWgX/DOY/aiG9uJhwPcMTfBt89ddfR3vvbSMlLiWH4zGfbOw25EZWxcrzHTcORL9v4/z86nWVwvetomNJMIddxsPvn6Grlyjftv3tP2524RNSNEbhWenbmS5aa97oy2VB52J+9/ab3Ad85d5R/cM7erpnNtRxvsFj7zndc5OBHlfQPov0s6GB/oZFKXws0BIPW65PuvrdrKXyuO7UkQDfobo0ztsrxe5KmLa7z/9v19/X4re5IR5m0ILLRyOe1MNAWavcd2awu+8coKb71hqqEO5oR3HZvh+flMI+1iFyklf/DEeY7Mxh0biGSk/7RRL5QiXz8hcbsG6W+eXyIU8PGeW/fa+j07TMWNHHa7+8njr6wQDvhsp6vsomZiuzGmtVipceFqznGFuGJuMkq+Sy/2Xz1zGSnhZ+7Z1Xpe2mArXry8zg9fX+Ohd1zf19xVu4xFgtT7HB+ohrE78rCTxo7cjXnDLy9tsF6o2A6HgzEE5E0Hx/s22H/7whJSwk/e6dxgv+PGaRbSBc5csTdmsl6XLK4XHRWcQbMg6TtnrQsGXVrNc3Zlk/tudscYvOumWeoSvnu+Py/7m6+meHlpg3/+94469lS8zGFfuJpHCGz12kYbOWzr79laXfLlF5b4seOzrmyoOjEdD1Ory8ZmXiGl5PHTV7j3xpnG+XvFdaYAzTlTE8AJ51Ob1CWOe7AV3Xqx59fy/D/fPM+PHpvpe3jQoNAG2+SR77xOLOTnv3vLdQN5PCcDQJo57P4FGGYSxhs844KS1PfMm3u/O/i7rpvg5cUNyn3kvr50apFb9o85LrgCeODWffgEfPmFZVu/t5orU67WHXvYs8kwd86N843TK70PNvnG6SsAroU77zo0QTIS4Fuv9qcy+KknzrN/PMIH73LuqSQjAXLlGjUXNpXbubiWY99YxFZ+NNbIYVv3sJ+6sMZKtsRP3nHA9jnaoZN4yvnUJvNrBe7zOBwORpumcKEGAgyFM3CmId6K6liYX9taeFat1fnVP38OgH/3029y5bG8RBtsjIKhvzm1yD88cciR12qHMQcjNtcLFUJ+X6OquB/cFE/53vlVjs7G2TvWXw73zrkJyrU6p5ftVScvpPM8eynDT7ngXYNhMN96wxRfeWHJ1u+50dKlePfxPTw7n7Gs8/746RWOzMYdica0EvD7uPfoDN96NWU7NXDywhpPXljjn/3oEVemUSmPdNODwrNLq3nb9RZTceN8lEiOFb50apFo0M/9t3hrMKfjSvhn6/tZbf4GYbAT4QA37Un2HS1r5czyJkG/cO26PtjBw/79J87x9MU0//anb9+Vymbb0QYb+PS3XkNK+Kc/6n2xmaI5jag/D3usz9GaCrUjd9ralStV+eFra9x7Y//5MTXC1O7O/G9PGYb1J9/knvfyE2/az9mVTc7aCIs/fdHosb3BYR80GDdWKeHvXu3tZW+ar73bxUQ/etMMi+tFzqfsDQL5g2+eZzIW5ENvdUclqqEn7kFY/MJq3nbf/pGZBNPxUEPRrxfVWp2vvLDMfbfscSRoY4VOHvbjr6xwy/4xx9Efq9x5aNxxmybAmeUNjs4mHHVdtKJ6sb/60jLfO3+VWl3y9MU1Pvn4WX767oM86EJEaBCMvMG+ulnic09d4oN3H3Qs9GAHteOz61UCrBfKjEed3QBmlYftcGLX376wRKFS4wN39m80D05EmUmEG+IOVvmbU4vcOTfOdS72tr73tn0IAV950VpYXErJn/zwIncdmnAlLP+mg+PMJMJ84/TOkPR6vrLlRvids1cp1+qu5a8V7zpmaPjbCYu/srTBN06v8E/uvcE14+TGGNp25EpVrm6WbF83Pp/gR45O891zVy0ZpB+8tsZqrsxP3eFOBKgbzRGbzQ34er7CyYtp7rt5cDMZ7jo0STpfaRT19curVzYdC6Zs56F3HObMcpb//g9/yDt+53H++R8/w4GJKL/54G2uPo6XjLzBfuQ7r1Oq1vmldx8d6OMeno5xeDpmK1+pMHTEnQ0QaITEHXrYnz+5wJGZuCPJRSEEd5k7c6tcuJrjxcsbrucG945FOHH9JF+2GBb//murnE/l+IW3X+/K4/t8gncfn+XvzqxsmU391IU17vm3j/Gxv3yhsf6N01dIRgKuj2s8NBXjyEycb9kofvv9J84RD/n58I8cdu08xjwasdlPhbji3htnWMmWGsNWuvE3zy8SD/l593Hvw9GTsZ0e9t+dTVGrS9c3dN1Q0TInYfFsscLlTMGxJOl2fv2B45z8Nz/O//3zd3PH3ARSwu996K7GxvBaYKQN9nqhwn/5/kV+4vb9rnhHdhBCcN/Ne/ne+VXbQgzrDnTEFePRIAGfcJTDvnA1x5MX1viZN885FmS4c26C86mcZW/qS6eMYW9/3wPv5f237+f0cpbXLNyU//gHF5mIBV09j/tu3sNGscozl4ybXq0u+fgXXiIc8PHnJ+f5pT95hkK5xjdOp3jXTbOuhQ1beddNs/zgtVVL08NevLzO355a4hfvvcFRq+F2nKSNuqF6sPuRsr3XLKz87rnuYfFytc5XX1rmgdv2uSL80Yug38dELLhFW+Ebr1xhKh7yVD98O8f3JokEfY4M9qtXNht/y21ioQA/decB/vDDJzj5b36cN1+/e0VS2jHSBvu/fP8C2VKVX/6xwXrXivtv2UO5Wu/55m/liTMrnFnOOpYY9PkEByejPHMp3fff+PzTC/gE/Mw9zqX8lIDKt15N8cylNF98fpGvvLDUMfT4pVNLnLh+0pPc3PtM4f9eYfGVjSJff+kKP/fmOVdvyu88NkPAJxrRl784Oc/LSxv8zs/cwW9+4Db+v1eu8FO//x2ubpY8E8P40WMzFCv1hgZ2N/7j188wHg3yz951xNVzUDns7a1KTrloqur1k0q5bjrG3GSU7/YQl/n22RTrhQo/OYBwuGLalCctVmo8cWaFJ86kePfxWceCQnZw2qYJ8HdnjOve7ZD4G4GRNdj5cpVHvnuBHzs+61iRqV/ecniKRDjQaM3pxXfOXuV//C9Pc9PeJL/+nuOOH/9/eNv1/OC1NZ6+aH9wWq1uzI19102zjvSzFWrU5v/0p8/yD/7ge/yLP3uWX/qTZ/jk4+d2HPvtsylOL2c9uxkemIhy93UTfOXF7mHxzz01T7Uu+UdvcyccrhiLBDlxeJJvnllhvVDhP37tDG85PMlP3bGfh95xmN/70N1cXM0hBJ6FW99+ZJqgX/TMYz/5+hrfPJP6/9u797CqynyB498fV0UQREBCIBRRUgsvCGmZZmmXqXTSTmYX83JqavKpps4pp6Y5HZ+ZpuyZbp4eTxfrmDWmk9NlctLM1ElNBcUQFURSREVFTE1TBN7zx1o4W9sbBfeFzf59nmc/LN691mbt3/Oy3rXeKw8MTXf7CIukmLbERITydTOajVypqa3nbxt2kxLbttlVoVekx/Ft2UGXw82MMby2tJTE9m0YnOG99uOO7cJZUVJFn/9ezIR31nGqrp5xOd4ZpuqoT0oMRc0cppm/s5r/WbadX1x6kV/02va2gCuwjTFsKD/EI3MLqD5Ww0PDuvnsXMJCgriqexxfbdl/zk4sa8oOMnn2OrrGtWPOpFy3VD3eeXkqse3CnBaK57KytIq9h09wm5vWjY2OCOWVsX149pZevD0+m8WPXsWt/Trz0pKSM9qTC3b9wK/ey6d7p0hu9eAk/Tf2vohNu49Q7mKO89q6ej5YU87gjDi3DT1xNCwzga2VR/ntgkKqj9fw+5t7nW52uCUriQ/+/XJeHJN1ei1kd2sXHkL2xbEsb6TANsYwfdFWEqLCGe/GtusGYSFBjO6XzKKiymYvynK215ZuY2vlUZ69pfkdjQZ168iRE7Vs2u18lbkvNlVSsOsHfjO8u1uGt52vgekdiY8K5/bsFN6dMID1vxvusdXBGpOVEkNNbdOHaf5wvIYpH2ygc0xbnhvd8sdE+0LAFNg/1dQxd205N732Db98fRUrS6uYMqybz9swhmV2Yv/RkxQ1skLSytIqJr67juQOEcyZnEsHN12kI8JCmDy4C8tLDjS5Cmt+fgUxEaFc29N9T3gj+3Rm/KA0rrmkE907RfHcrZfSLzWG38wrYNPuwxRXHuXed9bSMTKc9yblerSzyPW9EwkOEsa+sZo53+78WVvuV1v3U3nkBHe5qbPZ2RrGzX5euJfbs1Po3fnMWqABabGM9vCqQldnxrO18igvLip2+jS5rPgA63YcYso1GR6bReuOnFRq7dqcC1VYcZjXl21ndL/kC+qI1TBBkLPZ4E7V1TN9UTEZCZHc6uVpLh8d3p2vHx/KsyN7M7RHglfazp3p04z1AYwxPD7/Ow78eJIZ4/r6VUcwbwqYAnvP4Z94ckEhdfWGaaN6s+apa3lsxIVXK1+ooT3iEbHGS56trt7w8pIS7np7DUkxbflgcu7p3t3ucs/ANKLbhjJj6bbzPubw8VMsKqpkZFYS4SGeuyiEhwTzv3dn07FdOJP/L4+7315DWHAQcyblNnuSlvOVEhvBe5NySIxuw9Mfb2Lo9GX8eXExz3+xlakLCnlu4RYuim7jsTbk9PhIUmLbEhUewuPX+Saf3jMwjX/LTmbG16Xc+87aMyZz+fFkLdMXFZMaG8Ht2e6pZXGmW0IkOV1imbu2/IKm0T1ZW8dj8wuIiwzjmZt7XtA5xUeFk5kYxSonfU/m5e2irOoYT1yf6ZUpjlsia5hmWJOGac5auYMlW/Yx9YZLuCzZe53k/E3LXKXbA9LjI/nHw4PJTIzy6BJzTRUXGU7flBiWbt3Hw9dmnE4/cPQkj3y4gZWlB7m1b2emjertkUXVI8NDmHRlF/78ZQmbdh/+2ZOcMx/mlVNTW89tHrxQN4iPCufNe7IZM3MVYSFBzLt/oFvHXTdmUHocHz3QkW9Kq3h5yTZeXVpKSJAQExFGTEQoDw3L8NhFWUR4cUwW9Qa336SdrzahwbwwJot+qR145tMibnr1nwzpkcCG8kOU7DtKvT0sxtPVvuNyUnnkwwJWlx1s9gQ9r361jZJ9P/LOvQPc0tY+KD2O99fs5MSputNPssdranl5yTYGpHXw+MxmLZk1TDPmvJ+wP924hz8u3MLwnp2YcEWaZ0/OzwVMgQ3Wer8t0TWXdGL6omL2Hz1BQlQblm7dx3/+tZCjJ07xwujLuC37wodNNWb8oDTeXFHGjKWlzLy7f6P7frOtihe+KGZI93h6JXknnj2T2vPxr6+gTUiw1wrrBiLC4Ix4BmfEc+JUHeEhQV674cvt2vTFVDxhbE4qvZKimfKX9fx94x76pMYwolcil3eN9fgKUGA1T8R8FsoHa8ubVWAv2byPmcvLGNM/mavdVCNyRbeOzFr5PevLD52OwaxvvufA0ZPMvKtfi3oo8IWs5BiWbNl/ziGo8/J28cRH3zEgLZaXbu8T8HE7l4AqsFuqYZkJTF9UzGcb91JceYR5eRVkJkYxZ3IOmYmeLxSj24Yy4Yo0Xl1aykf5FS7bRrdWHuGBOfl0S4jktXF9vfrP5e5JFJrDV22CLcGlydF8/fhQjMHr6wW3CQ1mdL9kZq/eQdWPJ5tU4/B18X4efH89vZPa8/sLrAp3lNMlluAgYUVJFTW19czPr2BxUSUjenbyeb+YlqBPqlWtXVhx2OU67bNX7+CZT4oYnBHHG3dne3w1sdZAC+wWIDMxiqToNkz7+2aCBB4cms7D12Z4tH34bPcNSWfdjkM8Nn8j2w/8yOMjepxxYa48fIIJ76wjIjyYWfcO0E4hAUhE8NUD0B05qbz9zfd8lF/B/UPOb96EFSUHrGGQiZHMnpjr1uUto9qEkpUczczl25m53JpD/c7ci5niw1EnLUlDO3T+zkPkdo0lSIRTdfUU7TnM+p0/sHZHNV9u3sfwnp2YMa6vV691/szrBbaIXA+8AgQDbxlj/uTtc2hpRIRxuaksLKxk2qjeFzTNZ3NFhocwe1IOz3yyideXbafswDGm3phJWdUxtuw9woL1uzny0ynm/Wqg1xYSUKpBQ+ezN/9pTSV8ZUYcl3WOJiQ4CGMMx2rqOHSshoPHajh0rIby6uP8ceEW0uMjeW+ie4ZBnm3y4K58tnEPN2clcc0lCVroOIhuG0p6fDteWlLCS0tKfvZ+SmxbJl3ZhSdvyPTITH2tlVzoqipN+mMiwUAJMByoANYBdxhjNjvbPzs72+Tl5Xnt/JQ1vGLWyh384fPNOHbKTYltyx9GXcpV3b03EYRSjgp2/cDTHxdStOcIxlgzobULC6H6eI3TSToyE6N4f3IuHX3UaS/Q5e+s5tuyaurrDQYQoHtiFH1TY0iI8uwoD38nIvnGmOyfpXu5wB4I/Jcx5jr796kAxpjnnO2vBbbvrNtRzZa9R+jRKYrMi9p7bZ1wpc6l+lgNq7ZXsWr7QWrr6unQLozYiDA6tAujY7t//UyKaatPb8ovuSqwvV0l3hnY5fB7BZDruIOI3AfcB5Ca6v1p9ZRlQFosA3wwS5JS5xLbLoybLkty+0ptSrV03r79dNZl5YxHfGPMG8aYbGNMdny8Vr8qpZRS4P0CuwJwnG0jGdjj5XNQSiml/I63C+x1QIaIdBGRMGAs8KmXz0EppZTyO15twzbG1IrIQ8AirGFds4wxRd48B6WUUsofeX0ctjFmIbDQ239XKaWU8mc65kEppZTyA1pgK6WUUn5AC2yllFLKD3h1prOmEpEDwM4L/Jg4oMoNp9NaaXwap/FpnMbHNY1N4zQ+rl1sjPnZRCQtusB2BxHJczbFm7JofBqn8Wmcxsc1jU3jND5Np1XiSimllB/QAlsppZTyA4FQYL/h6xNo4TQ+jdP4NE7j45rGpnEanyZq9W3YSimlVGsQCE/YSimllN/TAlsppZTyA35ZYIvILBHZLyKbHNL6iMi3IlIgInkikmOn3yki39mvVSKS5XDM9SJSLCKlIvKkL76LuzUlNg7vDxCROhEZ45A2XkS22a/x3vwOntTU+IjIUDu9SESWO6S3urwDTf7fihaRz0Rkox2fCQ7HBFL+yRKR1SJSaMejvcN7U+08Uiwi1zmkB1L+cRofERkuIvl2er6IDHM4pr+dXioir4qI+OL7tDjGGL97AVcB/YBNDmmLgRvs7RuBZfb2IKCDvX0DsMbeDga2A12BMGAj0NPX382bsXGIw1KsBVnG2GmxQJn9s4O93cHX380HeScG2Ayk2r8ntOa804z4/BZ43t6OB6rteARa/lkHDLG3JwLT7O2edt4IB7rYeSY4APOPq/j0BZLs7d7Abodj1gIDAQH+0ZD/Av3ll0/YxpgVWBeHM5KBhjvbaGCPve8qY8whO/1bINnezgFKjTFlxpgaYC4w0qMn7gVNiY1tCvARsN8h7TrgS2NMtR27L4HrPXPG3tXE+IwDFhhjyu1jG2LUKvMONDk+Boiyn34i7eNqCbz80wNYYW9/CYy2t0cCc40xJ40x3wOlWHkn0PKP0/gYYzYYYxryUhHQRkTCReQioL0xZrWxSu/ZwCjPn33L5/XlNT3oEWCRiLyIVdU/yMk+k7Du1gA6A7sc3qsAcj16hr7jNDYi0hn4JTAMGOCwv7PYdPbOqfqEq7zTHQgVkWVAFPCKMWY2gZV3wHV8ZgCfYhXgUcDtxph6O18FUv7ZBNwCfALcBqTY6Z2xHhIaOMYhkPKPq/g4Gg1sMMactPNPhcN7rT3/nDe/fMJ24QHgUWNMCvAo8LbjmyJyNVaB/URDkpPPaK1j3FzF5mXgCWNM3Vn7B1JswHV8QoD+wC+wnhp/JyLd0fg0xOc6oABIAvoAM+z2yUCLz0Tg1yKSj3XjUmOnu4qDxseBiPQCngfub0hy8hmtOT7nrTUV2OOBBfb2fKxqJwBE5DLgLWCkMeagnVzBmXd6yZxZVdyauIpNNjBXRHYAY4DXRWQUgRUbcB2fCuALY8wxY0wVVrVeFhqfhvhMwGoyMMaYUuB7IJMAi48xZqsxZoQxpj/wF6z2aXAdB42Fk+euAAABlElEQVSPTUSSgb8B9xhjHOOW7PARrTo+TdGaCuw9wBB7exiwDUBEUrEuNncbY0oc9l8HZIhIFxEJA8ZiVe+1Rk5jY4zpYoxJM8akAX8FHjTGfAwsAkaISAcR6QCMsNNaK6fxwarCGywiISISgVVtuYXAyjvgOj7lwDUAItIJq62yjADLPyKSYP8MAp4GZtpvfQqMtdtluwAZWJ2pAir/uIqPiMQAnwNTjTErG/Y3xuwFjorI5Xb/iHuw/heVr3u9NeeFdZe2FziFdTc2CbgSyMfqcbkG6G/v+xZwCKvqrgDIc/icG4ESrDu+p3z9vbwdm7OOexe7l7j9+0SsTjKlwARffy9fxQf4D6ye4puAR1pz3mlqfLCqwhcDhXZ87grQ/POwnRdKgD9hzyBp7/+UnUeKcejpHGD5x2l8sArvYw7X5gL+NRIj285T27H6SoivvlNLeunUpEoppZQfaE1V4koppVSrpQW2Ukop5Qe0wFZKKaX8gBbYSimllB/QAlsppZTyA1pgK6WUUn5AC2yllFLKD/w/YvhA4mTeACsAAAAASUVORK5CYII=\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",
    "$$\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",
    "\n",
    ", where `level` and `s` follows the following recursion rules:\n",
    "\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}\n",
    "$$"
   ]
  },
  {
   "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": [
    "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 `jax.lax.scan` instead. A detailed explanation for using this utility can be found in [lax.scan documentation](https://jax.readthedocs.io/en/latest/_autosummary/jax.lax.scan.html#jax.lax.scan). Here we use it to collect expected values while the pair `(level, s)` plays the role of carrying state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scan_exp_val(\n",
    "        y, init_s, level_sm, s_sm, coef_trend, pow_trend, pow_season, future=0):\n",
    "    N = y.shape[0]\n",
    "    duration = N + future\n",
    "    seasonality = init_s.shape[0]\n",
    "\n",
    "    def scan_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 = moving_sum + y[t] - \\\n",
    "            jnp.where(t >= seasonality, y[t - seasonality], 0.)\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",
    "        return (level, s, moving_sum), exp_val\n",
    "\n",
    "    level_init = y[0]\n",
    "    s_init = jnp.concatenate([init_s[1:], init_s[:1]], axis=0)\n",
    "    moving_sum = level_init\n",
    "    (last_level, last_s, moving_sum), exp_vals = lax.scan(\n",
    "        scan_fn, (level_init, s_init, moving_sum), jnp.arange(1, duration))\n",
    "    return exp_vals, last_level, last_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our utility function defined above, 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 `substitute` 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": 5,
   "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",
    "    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(\"offset_sigma\", dist.TruncatedCauchy(\n",
    "        low=1e-10, loc=1e-10, scale=cauchy_sd))\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",
    "    exp_val, last_level, last_s = scan_exp_val(\n",
    "        y, init_s, level_sm, s_sm, coef_trend, pow_trend, pow_season, future=future)\n",
    "    if future == 0:  # training\n",
    "        omega = sigma * exp_val ** powx + offset_sigma\n",
    "        numpyro.sample(\"y\", dist.StudentT(nu, exp_val, omega), obs=y[1:])\n",
    "        # we return last `level` and last `s` for custom forecasting\n",
    "        return last_level, last_s\n",
    "    else:  # forecasting\n",
    "        exp_val = exp_val[y.shape[0] - 1:]\n",
    "        assert exp_val.shape[0] == future\n",
    "        omega = sigma * exp_val ** powx + offset_sigma\n",
    "        numpyro.sample(\"y\", dist.StudentT(nu, exp_val, omega))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that all prior parameters are retrieved from [this file](https://github.com/cbergmeir/Rlgt/blob/master/Rlgt/R/rlgtcontrol.R) in the original source."
   ]
  },
  {
   "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": 6,
   "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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                      mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "      coef_trend     31.71    109.28     13.07    -94.03    157.59   2223.27      1.00\n",
      "       init_s[0]     87.94    111.74     63.44    -68.05    235.47   3378.28      1.00\n",
      "       init_s[1]    -19.40     71.06    -25.47   -126.41     91.89   5181.19      1.00\n",
      "       init_s[2]     31.82     91.29     21.08   -106.81    176.86   5820.17      1.00\n",
      "       init_s[3]    126.90    126.38    110.88    -57.80    314.41   3946.57      1.00\n",
      "       init_s[4]    450.66    250.66    402.09     80.56    809.16   3572.03      1.00\n",
      "       init_s[5]   1178.64    477.78   1086.83    469.48   1872.74   2261.22      1.00\n",
      "       init_s[6]   1980.18    667.78   1869.90    948.43   2969.04   1685.41      1.00\n",
      "       init_s[7]   3683.29   1110.94   3576.69   1991.56   5381.30   1788.64      1.00\n",
      "       init_s[8]   2594.08    830.16   2458.22   1340.25   3902.45   1903.57      1.00\n",
      "       init_s[9]    950.58    442.95    872.55    297.35   1586.94   3819.57      1.00\n",
      "      init_s[10]     48.05    115.09     29.45   -104.54    203.32   1944.86      1.00\n",
      "      init_s[11]      0.74     51.49     -1.21    -78.70     71.10   1051.24      1.01\n",
      "      init_s[12]     -9.68     64.99    -13.42   -110.83     86.89   6341.63      1.00\n",
      "      init_s[13]     67.36    108.94     47.79    -82.89    217.10   1968.48      1.00\n",
      "      init_s[14]    328.46    250.74    277.57    -27.40    685.78   3291.15      1.00\n",
      "      init_s[15]    969.34    382.30    900.11    380.28   1547.19   2617.55      1.00\n",
      "      init_s[16]   1249.08    466.92   1180.00    571.64   1962.99   1434.37      1.01\n",
      "      init_s[17]   1393.10    561.44   1287.51    512.45   2198.97   2642.67      1.00\n",
      "      init_s[18]    631.61    309.09    578.08    168.21   1064.59   2431.41      1.00\n",
      "      init_s[19]     17.99     91.00      5.65   -125.10    142.38   4486.76      1.00\n",
      "      init_s[20]    -30.45     65.18    -24.85   -144.74     59.57   6396.10      1.00\n",
      "      init_s[21]    -15.73     46.61     -5.95    -91.76     44.06   2125.94      1.00\n",
      "      init_s[22]     -2.63     39.24     -2.19    -66.78     54.17   8246.59      1.00\n",
      "      init_s[23]     41.04     86.69     27.38    -85.88    162.45   4360.63      1.00\n",
      "      init_s[24]    536.16    335.72    477.67     49.87   1006.04   1714.93      1.00\n",
      "      init_s[25]    929.64    438.76    846.44    250.63   1568.07   2613.02      1.00\n",
      "      init_s[26]   1809.15    700.33   1706.28    743.22   2817.44   1873.94      1.00\n",
      "      init_s[27]   1267.05    477.93   1186.64    577.52   1981.32   2563.01      1.00\n",
      "      init_s[28]    216.87    169.62    188.16    -32.93    455.16   5057.82      1.00\n",
      "      init_s[29]     -7.31     85.30    -17.81   -130.93    124.72   2685.79      1.00\n",
      "      init_s[30]     -2.51     88.73    -14.41   -140.48    119.15   4614.12      1.00\n",
      "      init_s[31]    -33.77     73.45    -34.27   -153.91     67.47   1450.09      1.00\n",
      "      init_s[32]     -8.69     84.26    -16.32   -137.97    113.76   2231.66      1.00\n",
      "      init_s[33]    117.30    139.20     89.81    -81.13    316.03   5162.83      1.00\n",
      "      init_s[34]    517.04    290.37    465.07     85.27    930.02   3820.63      1.00\n",
      "      init_s[35]   1065.83    452.55    986.97    385.24   1706.04   2613.49      1.00\n",
      "      init_s[36]   1842.99    643.31   1743.51    862.68   2815.47   2146.56      1.00\n",
      "      init_s[37]   1440.90    556.17   1350.66    611.01   2260.78   1734.89      1.00\n",
      "        level_sm      0.00      0.00      0.00      0.00      0.00    696.21      1.01\n",
      "              nu     12.30      4.88     12.66      5.29     19.97    915.91      1.01\n",
      "    offset_sigma     33.43     31.75     23.41      0.00     78.19    508.38      1.01\n",
      "      pow_season      0.09      0.04      0.09      0.02      0.15   1249.64      1.00\n",
      "  pow_trend_beta      0.26      0.17      0.23      0.00      0.50    290.54      1.01\n",
      "            powx      0.62      0.13      0.61      0.41      0.84   3712.07      1.00\n",
      "            s_sm      0.09      0.10      0.06      0.00      0.21    242.35      1.02\n",
      "           sigma      9.64      9.39      6.94      0.40     20.33   5204.77      1.00\n",
      "\n",
      "Number of divergences: 4411\n",
      "CPU times: user 1min 20s, sys: 247 ms, total: 1min 20s\n",
      "Wall time: 48.7 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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictive = Predictive(sgt, samples, return_sites=[\"y\"])\n",
    "forecast_marginal = predictive(random.PRNGKey(1), y_train, seasonality=38, future=34)[\"y\"]"
   ]
  },
  {
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sMAPE: 64.12, rmse: 1252.74\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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAEICAYAAABLWh2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9d7wkV3nn/T2d7+2b7+SkkTSjTBAaCdnEBZNtwr6AYUm2sVlYnNZ+1xh71zgb/NqAWb9gA8IWYQnGAYwJkglCYCQ0kkBhRtLMaPLMzbH7dqqqs3+cU93Vfbtv6njVz/fzuZ/bXVXdVV3VXb/zPOcJSmuNIAiCIAibl1C7D0AQBEEQhPoQMRcEQRCETY6IuSAIgiBsckTMBUEQBGGTI2IuCIIgCJscEXNBEARB2OSImAttRyn1BqXUbU1671NKqZ9qxnuv4xh+Xyn16XYeQ7NQSu1TSqWUUuEVttFKqQOtPK5OQCn1XKXUuTVuu+J3RCkVV0odUUrtaNwRbg6UUi9XSn2u3cfR6YiYb1KsSGXsjdT/29Xu41oNpdR+e3OP+Mu01p/RWr+wncfVKSil/l4p9cebZT9a6zNa6z6ttWvf9ztKqV+s47iGlFKfUEqNKaUWlVKPKaXeFVivlFK/rJR6QCm1ZLf7jlLqdXb9w4Hfg6uUygae/069n7eNvA34rtZ6DIrn6Val1IT9+/3gxvZ39m17jh4JDmiVUs9XSp1USl1USv1sYPmQUuo+pVR/rYOoNUAJXne7jWfP+aJS6lGl1M8HjksHrsm4UuorSqkXVLxfcRCutf4ycJ1S6skbOG9dg4j55uZn7I3U/7uwnhcHBVUQOoQPAH3A1cAg8HLgRGD9h4BfB34TGAV2A/8TeDGA1vpa//cA3An8cuD38aet+xgN578Cnwo8/wDQC+wHbgLe5Aum5bPA/Zhz9LvAF5VSW+26DwI/gzlnHwl4Vf4MeK/WerEBx3vBXoMB4F3Ax5RS1wTWD9n1TwFuB/5ZKfVzK7zfZzEDGqEWWmv524R/wCngp2qseznwMDAHfAe4uuJ17wIeAHJABNgF/CMwCZwEfjWwfRj4HcwNdRG4F9hr1/0VcBZYsMufFXjdTcBhu24ceL9dfgbQQMr+/QTwc8D3Aq/VwNuBY8As8P8DKnA8fwlM2WP9Zbt9ZKXzBOwAloDRwLob7GeO+scA/IXd50ngJXa7EeAcZvAERmyOA2+usc9LgTvs+bod+Gvg04H1/wCMAfPAd4Fr7fK3AQUgb8/Nv9rlvx04/0eAVwXe64Dd17w9J58PrLvK7n8GeBR47Ur7qfgMfwD8b/s4CqSBP7fPe4AsMIwRE435Hv0J4Np1KeCvV7ueVfb7EPDKGuuusO9/aI2/ke8Av7jKNr9vr8en7fl90O7n3cAE5vv9wsD2u4Av23N6HPilwLoe4O/tZzwC/A/gXMVra/3Ofj/4Hak4xn1AhsB33F7rGwPPfwe4M3CeckB/YP2dwNvt48cDy8eAbZjf69fXcE6fG/xM1c51tW3sZ3518PtSsf7/xdwnQtXub8AzgJNrvT9241/bD0D+Nnjhaoi5/SGngRdgbsK/ZW86scDrfgTstTefEEaIfw+IAZcBjwMvstv/D3uDuxJQmJH0qF33RszIP4KxlMaAhF33A+BN9nEfcLN9vOzHTHUx/wowZG9kk8CL7bq3Y26UezBi8u/Vbg7VzhPwVeAdgXUfoCRYP4cRuF/CDBjeAVygNIh4YeDG9zHgiytcmx8A7wfiwLMxIhEU818A+u36DwI/Cqz7e+CPK97vNRghCAE/a6/vTrvusxjLKwQkgGfa5UmMEP28vT5PwwjAtbX2U7HP5wEP2sc/iRlM3B1Y9+Nq15MqArrS9ayy349jBqI/DxysWPd24NQ6fiPLjqXKNr+PGXy8yJ6nT2KE9ncxv59fIiAimIHTh+25fqr9LM+3696LEc0RzO/rIayosfrv7PepLeYvAx6uWDYF3BR4/rvArH38KuBoxfZ/Tem7fhfmd/wUzHc8CvwHcMUazulzWYeY28/9Ksxv68rK70vg9ZfZ5VdXu7/Zc6qBgbVe/277Ezf75uZflFJz9u9f7LKfBf5Na3271rqAsTR7MDdknw9prc9qrTPAjcBWrfUfaq3zWuvHMWL1OrvtLwL/U2v9qDb8WGs9DaC1/rTWelpr7Wit/xIjTlfa1xWAA0qpLVrrlNb6rnV+tvdqree01meAb2NunACvBf5Ka31Oaz2LuYGulVsxAxCsa/H1lLsuT2utP6bN/O+twE5gu/2st2EsuG9ibq7/tdoOlFL7MOf0f2mtc1rr7wL/GtxGa/0JrfWi1jqHuYk/RSk1WOugtdb/oLW+oLX2tNafx1i4N9nVBeASYJfWOqu1/p5d/tMY4fs7e33uw1iFr179NAFmQHJQKTWKGZDcAuxWSvUBz8GI2nqodT0r+RXgMxiPyxGl1HGl1Evsui2YAVURpdQ5+/3PKqUuWecx+dyptf6G1trBXOOt9ngLwOeA/XY+eS/wTOBd9lz/CDP4eJN9n9cCf6K1ntFan8VMCfis9jtbiSHMgDDI14HfVkr12+DCX8C43cEMnucrtp/HDCDBDIr+CvioPfZ3YL7XCaXUN+xc+3NWOJ5dgfvOnFJqDnNelm2DGXS8BzOwf3SF9/SnCEdqrPc//9AK79HViJhvbl6ptR6yf6+0y3YBp/0NtNYexkLbHXjd2cDjS6j4cWJcdtvt+r2Uz1kWUUr9plLqqFJq3r5uEHPDBXgrxkvwiFLqHqXUT6/zswVv2kuYG5T/+YLHH3y8Gl8CrlFKXYbxXMxrrX9YbZ9a6yX7sC+w/qPAdcDf+QOaKuzCWEjpwLLi9VBKhZVS71VKnVBKLWAsECidt2Uopd6slPpR4PpcF9j+tzAekx/a4K9fsMsvAZ5ecV3fgJluWBU70DuMEe5nY8T7PzDuzo2Iea3ruWy/Wus/1VrfgPH6fAH4B6XUCDCNGWAFt9+DORdxzHnYCOOBxxlgyg7o/OfY490FzOjyOeXTlH5bld/N04HHq/3OVmKWkhD7/Ko9tmOY7/VnMVNBYKY4Biq2H8AKotb6R1rr52qtn47xcv0C8KeYgckfYLwin1JK1TqfFwL3nSGt9RBmiqraNiNa66dqrVeLRvfP4UyN9f7nn1vlfboWEfMnHhcwNw7ARP9iBPl8YBsdeHwW40YM/jj7tdYvDay/vHInSqlnYebeXwsM2x/0PPaGqrU+prV+PcYt/T5MAE6yYt8b4SLGxe6zd60v1FpnMeLwBoxF8qmVX1HCWvJ/i3HDvkPVTrW6CAzbz+qzL/D4vwCvwMzjD2LcjlASorLzY63Nj2Es1VF7nh+idJ7HtNa/pLXehfEWfNge21ngjorr2qe1fke1/dTgDoxL/XrgHvv8RRivwHdrvKbe61t6I60XMCKTxMQhfAvYo5Q61Kh9rJMLwEhFtPc+Sr+ti5R/H4PXfbXf2Uo8AFymyjNAZrTWb9Ba79BaX4u5l/sD04ft9sHjfIpdXskHMJ63DPAk4LDW+hTG9b61yvbN4lWYGIVa1vvVGE/TQusOaXMhYv7E4wvAy2z6SRQzl53DWFXV+CGwoJR6l1Kqx1qO1ymlbrTrPw78kVLqoDI82bpe+wEHM2cYUUr9HgFrQCn1RqXUVusZ8EfTrt3ew8yRbfTz/ZpSardSaggzoFgPn8TMj78cE/S0Vvy0pl/ATF18UlXJrdZan8ZYtH+glIoppZ6JiRz26cdcj2mMW7Qywnqc8nPjD4AmAWzE8nX+SqXUa5RS/uBm1m7rYuaor1BKvUkpFbV/Nyqlrq6xn2rcAbwZOKK1zmPnRTGiNFnjNWt535oopf6XPc6YUioB/Brm+/OoddP+LfA5pdQL/O8r5VNITcO6zv8D+DOlVMKmSr0VMy0A5rv5bqXUsL0mvxJ4+Wq/s5X2e47yqRWUUpcrpUbt+7wEE9T4x3b7xzBxMe+xx/kq4MmYaRYC7/ECTIzLV+yik8DzlFLXYjwdtbxPDUMptV0p9csYV/y77f2iGs8Bvtbs49nMiJg/wbA3vDcC/xszX/UzmCjsfI3tXbvNUzE/5imMgPtzuO/H3KRuw0Sm34KZg/8G5sf1GMadmKXcxfhi4GGlVAozP/c6O8+4hIl6/r51N968zo/4MXssD2BSb76KGVS4K70o8Hm/jxlM3GctkFVRSt0A/AYmet3FeBo0Jsq8Gv8FeDrGZfgezADC55OY83Ue4+KsjCW4BTMVMKeU+het9RFM9P4PMEL5JOD7ge1vBO625/nLwK9prU9aV/ALMXOyFzBu7vdhbtLL9lPjc/wH5lr7VvgRzHWuZZWDudavVkrNKqU+tMJ2tdDA32G+hxcw0yEv01qn7Pp3Yuai3485v+eAP8LEipzZwP7Wy+sx3pQLwD8D79Fa327X/QHm2p7EfEeLnp81/M5W428pzc2DycR4EOM6/zPgDVrroOX9OuAQZoD3XuDVwQGYUioO/H+YwZLPrwB/gwkq/W+BqYZmMKeUStvP8FLgNVrrT6yw/esx50CogR+pKwibEmuV/I3Wes3BT0qpbwH/R2v98eYdmSA0Diu+92Mi5y+2+3haiVLqZzABdK9t97F0MiLmwqZCKdUD/CeM5bMd4zq8S2v962t8/Y2Y3Ou9ujHFMQRBENqOuNmFzYbCuDNnMZbKUUzu7uovVOpWjAvx10XIBUF4IiGWuSAIgiBscsQyFwRBEIRNzqqNNpRSn8BUk5rQWl9nl40An8dEdZ7C1HyetTnNf4WJTlwCfs5WnkIp9RZMQwQwZSRvtctvwJSW7MFEJv+aXoO7YMuWLXr//v1r/ZyCIAiCsKm59957p7TWVfP/V3WzK6Wejako9MmAmP85phLSe5VSv40pGvIupdRLMekNL8Wk5vyV1vrpVvwPY1IlNKZG8Q12APBDTHrEXRgx/5DWetV8wkOHDunDhw+v5fMLgiAIwqZHKXWv1rpq0aRV3eza1JauLLH3Ckztauz/VwaWf1Ib7gKGlFI7MVWjbrdVi2Yx0cQvtusGtNY/sNb4JwPvJQiCIAjCGtjonPl2P9fR/t9ml++mvHDIObtspeXnqiwXBEEQBGGNNDoArlphfr2B5dXfXKm3KaUOK6UOT07WqiYpCIIgCN3FRsV83LrIsf8n7PJzlDca2IMpe7jS8j1VlldFa/1RrfUhrfWhrVtb2QNAEARBEDqXjYr5l4G32MdvwbTg85e/2TbkuBnTYvIipo73C20DgmFMzehv2HWLSqmbbST8mwPvJQiCIAjCGlhLatpngecCW5RS5zCNI94LfEEp9VZMc4PX2M2/iolkP45JTft5MO36lFJ/hGmjCPCHWms/qO4dlFLTvoZ0xhEEQRCEdbFpK8BJapogCILQTdSVmiYIgiAIQmcjYi4IgrAJ8LzN6UUVWoOIuSAIwiYg73q4IuhCDUTMBUEQNgGupym4XrsPQ+hQRMwFQRA2AY6IubACq6amCYIgCO3H9TRuuw9C6FhEzAVBEDYBrqfZrKnEQvMRMRcEQdgEuJ7G8cTNLlRH5swFQRA2Aa7WFFyxzIXqiJgLgiBsAlxX40gAnFADEXNBEIRNgKs1eRFzoQYi5oIgCJsA1/NwxM0u1EDEXBAEYRMgeebCSoiYC4IgbAJcT9zsQm1EzAVBEDYBrqfFzS7URMRcEARhEyC12YWVEDEXBEHYBJgKcEh6mlAVEXNBEIRNgGPbn0rhGKEaIuaCIAibAM+KuQTBCdUQMRcEQdgEuLbJirjZhWqImAuCIGwCxM0urISIuSAIwibAK4q5WObCckTMBUEQOhw/kh1EzIXqiJgLgiB0OME+5uJmF6ohYi4IgtDhBLRcLHOhKiLmgiAIHU65ZS5iLixHxFwQBKHDKbfMxc0uLEfEXBAEocMJWuaSZy5UQ8RcEAShw/ELxoBUgBOqI2IuCILQ4bheScwdV6O1uNqFckTMBUEQOpzKPuYyby5UUpeYK6X+u1LqYaXUQ0qpzyqlEkqpS5VSdyuljimlPq+Uitlt4/b5cbt+f+B93m2XP6qUelF9H0kQBOGJhacrxVxc7UI5GxZzpdRu4FeBQ1rr64Aw8DrgfcAHtNYHgVngrfYlbwVmtdYHgA/Y7VBKXWNfdy3wYuDDSqnwRo9LEAThiYbjiZgLK1Ovmz0C9CilIkAvcBF4HvBFu/5W4JX28Svsc+z65yullF3+Oa11Tmt9EjgO3FTncQmCIDxh8JaJubjZhXI2LOZa6/PAXwBnMCI+D9wLzGmtHbvZOWC3fbwbOGtf69jtR4PLq7xGEASh6xHLXFiNetzswxir+lJgF5AEXlJlU/9bqGqsq7W82j7fppQ6rJQ6PDk5uf6DFgRB2IS4FWJeGRAnCPW42X8KOKm1ntRaF4B/An4SGLJud4A9wAX7+BywF8CuHwRmgsurvKYMrfVHtdaHtNaHtm7dWsehC4IgbB4qxVxyzYVK6hHzM8DNSqleO/f9fOAI8G3g1XabtwBfso+/bJ9j139Lm2TJLwOvs9HulwIHgR/WcVyCIAhPKMTNLqxGZPVNqqO1vlsp9UXgPsAB7gc+Cvwb8Dml1B/bZbfYl9wCfEopdRxjkb/Ovs/DSqkvYAYCDvBOrbW70eMSBEF4olGZmiZudqEStVkrCR06dEgfPny43YchCILQdH5wYpp0zik+H07GuOGS4TYekdAOlFL3aq0PVVsnFeAEQRA6HCkaI6yGiLkgCEKHUzlnLm52oRIRc0EQhA7H9cotcbHMhUpEzAVBEDoYrTUVWo7r6WVV4YTuRsRcEAShg6nMMfeRXHMhiIi5IAhCB1M5X+4jrnYhiIi5IAhCB1MZye4jQXBCEBFzQRCEDkYsc2EtiJgLgiB0MG4NC7wgAXBCABFzQRCEDsat4WYvOGKZCyVEzAVBEDqYWtHsTmW+mtDViJgLgiB0MLXFXNzsQgkRc0EQhA6mpphLNLsQQMRcEAShg6llgdcSeaE7ETEXBEHoYHzR/ps7TvCZu0+Xlm/S9tVCcxAxFwRB6GB8MT87s8Q9p2aLRWTEMheCiJgLgiB0ML5oZx2PVM7h3EwGkDlzoRwRc0EQhA6mKOYFF4AjFxeA2mVehe5ExFwQBKGDcbXG05qcLRJz1Iq5pKYJQUTMBaEDSOcc3vf1R4rWlyD4uJ5HrmCEPBJSHJtIUXA9XCkaIwQQMReEDuCux6f5yHdOcN+Z2XYfitBhOK4m65hB3rW7Bsi7HicmU3geaHG1CxYRc0HoAFI5B4C5pUKbj0ToNFytix6bJ+8ZIqRK8+biahd8RMwFoQNYzBoxn13Kt/lIhE7D9TRZ62Yf6o1y6ZYkRy8uFtcJAoiYC0JHkBbLXKiBEXNjmSciYa7ZOcCp6TRLeUfEXCgiYi4IHUDJzS6WuVCOF3CzJ6Ihrt45gNbw6NiiuNmFIpF2H4AgCDJnLlTH8zSeZwrGACSiYXYP9RCPhDhycQFPxFywiGUuCB1AqjhnLmIulHDKCsZoEpEQkXCIg9v7OHpRLHOhhIi5IHQA6by42YXl+FXecvk8/xr7XZ55+JcBuGbnAGMLWS7MZdp5eEIHIWIuCB2AH80+lxHLXCjhW96x7AxPCp1i29gdKM/hwNY+AI6OLbTz8IQOQsRcEDqAtATACVXwo9XD+XkAFJpYdor+RBSAeZmWESwi5oLQAaRzJlp5bqkgVb2EIr6YR/IlCzyeGScZDwPiyRFK1CXmSqkhpdQXlVKPKKWOKqV+Qik1opS6XSl1zP4fttsqpdSHlFLHlVIPKKWeFnift9jtjyml3lLvhxKEzYYfze54uvhYEBxbfz1aCIh5doKeaJiQgnnx5AiWei3zvwK+rrW+CngKcBT4beCbWuuDwDftc4CXAAft39uAjwAopUaA9wBPB24C3uMPAAShW0jlHJIxa22J61Sw+L1UYkExz0yglKI3FmFeLHPBsmExV0oNAM8GbgHQWue11nPAK4Bb7Wa3Aq+0j18BfFIb7gKGlFI7gRcBt2utZ7TWs8DtwIs3elyCsNnQ2ljje4Z7ASnpKpTwLfOYs1hcFs9MANAXj7CQES+OYKjHMr8MmAT+Til1v1Lq40qpJLBda30RwP7fZrffDZwNvP6cXVZr+RMKrbUUeBCqknM8XE+zZ7gHEMtcKOFb5j3ucjFPxsMs5OS7IhjqEfMI8DTgI1rr64E0JZd6NVSVZXqF5cvfQKm3KaUOK6UOT05Orvd428pvfOHH/Orn7m/3YQgdiJ+WtndELHOhHN8y7w2IeSxrxTwmlrlQoh4xPwec01rfbZ9/ESPu49Z9jv0/Edh+b+D1e4ALKyxfhtb6o1rrQ1rrQ1u3bq3j0FvPwxfm+Y8T0xKpLCzDT0sTy1yoxC8a0+ulissSmXEAkvFIsXKgIGxYzLXWY8BZpdSVdtHzgSPAlwE/Iv0twJfs4y8Db7ZR7TcD89YN/w3ghUqpYRv49kK77AnFdCrPTDrP5GKu3YcidBh+9PruIRFzoRy/aExSl8Q8FnCzS+aD4FNvo5VfAT6jlIoBjwM/jxkgfEEp9VbgDPAau+1XgZcCx4Eluy1a6xml1B8B99jt/lBrPVPncXUUjusxY12nR8cW2TaQaPMRCZ2Ef0Me7I3SH4+Im10o4noax/UYIl1cFsvPEXJzJOMRMgWXgusRDUvJkG6nLjHXWv8IOFRl1fOrbKuBd9Z4n08An6jnWDqZmaU8vnf96MUFnnPF5poiEJqL7yrti0cY7I1KFTihSMHVZAsegwExB2OdJ2NxwHhytvbH23F4Qgchw7kWMLVYujk/clFqKQvl+E1W+uIRhntjUtVLKJIruGQdlwFlxDyX2AKYwjF+Fbj5jAz+BBHzljCVMvPkgz1RHhlbXGVrodtYDFjmQ71RaYMqFMm7HtmCywBLAKQHDgAmPa0vbhyrs2kRc0HEvCX4Yv6Tl49yfCJFznHbfERCJ+FHs/clIgz1xqREpwCY2hR5xyOXy9OvMmgU6f7LAFufPWbEfKaLxdyV2h1FRMxbwHTK/NiedXArjqc5MZFe5RVCN5HKOYQU9ETDDItlLlhyjofWoLOmY1ou0k+udwdgLPOktcynu1jMF7Od9VtxXK9t+xYxbwFTqRyxSIgb95uS849ID2IhQCrnkIxHUEox1BtjIVsQi0Mg7wtDbs48jw6Q6zEFNY2Yd3ctf611cYqqE/A8Xcxaagci5i1gMpVjSzLGpVuSxCIhjkoQnBAglXWK859DPVG0hgUJgut6cgUj5iFrmTvRfnKJ7UB557RuTWUsuJqlfOdMWeZdr62/WxHzFjCVyrOlP04kHOKK7X0SBCeUkc47RZfpcDIKdO8NWijhW+bhnBXz2GDAMh8vdk7rVss873pkCp0j5rmC19YudiLmLWBqMceWPpMHevWOAY5eFDEXSiyWWeYxAJk3F8hZoYrkjZh78cEyNzuYKnDd+l0pOB6ZDrLMc47LQtZpW8luEfMWMJ3OsaXP3KSv2jnAVConZV2FIulcQMx7jWUuucNCzjGWeSRvBv9OfBAnOoAbThBx0oQLKdtspTvFPO96ZDsoMyjneLiuJt2mAYaIeZPxPM10Kh+wzPsBCYITSqQCYj7cay3zdHfeoIUSvpjHHHOvcKIDoFSZdd4Xj7TVtdtO8lY88077IsiD+CnH7boeIuZNZj5TwPE0o1bMr9o5AMAj4moXLOmcW5wz9y1zmTMXfJGK+2IeM/eOXMKKedakp813WHpWq/BjCjpl3jxrAxbn2zTtIWLeZPyCMb6bfSQZY/tAXCLahSKL2QL9CSPmA4koIdW+0b3QOfiWXsIxA/9CdNAsr0hP61Y3e8GKebZDxNz3pCy0aXAlYt5kJq2Yb+0rNUK4ascARyWiXcDkyqbzbjFnOBRSDPZExTLvcvzqbwAJ186Zx8wUXZmYxyIs5d2isHUT/vnplCA4f/CVzjltqRMhYt5kpmz1ty2BrkZX7ezn+MRiV/4AhXKyBQ/X0/TFo8VlQ72xro1QFgx51yt2Wuz1TC/zQsy3zG2ueaAKXDd6coqWeYcEwfmWebvqRIiYN5npopu9JObX7Byg4GpOTKbadVhCh+D3Mu+zljmYefN2zbsJnUEuENSVtGLuRO2cuRXzWHa8q6vA5TrIMnc9jeuWrPF2uNpFzJvMVCpHOKQY6ilZXldsN+6yY+Mi5t1OsMmKz3BvTNzsXU4wQjupKyxzGwCXCDRbmevC70vBimcnBMBVNs9qh6dExLzJTC3mGUnGCIVUcdlW63Lv5m5HgsG3zP2bMpiSrt1oaQklgpb5gDaNmUqWuRHzWMDN3m3fF611salJJwTA+aV3fRYyra8ZL2LeZKZSuTIXO1C00sX6ElJVLPOh3lhXWlpCCb/6m/Zc+rFibgPg8omtgM0zj1k3e5fNmQdjCjyv/YKeq8h1zxbclre6FjFvMkbMY2XLIuEQA4kIs2KZdz2prD9nHhTzKOm82zHFMITW44uDl1skrDTZUA86ZIwAN5rEifYR9vIMhYwLvtsGfwW3PFq80jJuNdWEu9WudhHzJjOVypelpfmMJCViWTBNVqBczIdt4Zg5KenatfgDOW/J1GXPhvvL1mdtENygM41S3edmrxzotnvevNIyh9a72kXMm4jW2ljm/cvFfEiCnAQo9mMut8yNJ6fbbtBCiWKaU2bWPI+Ui3kxCC47STIW6bqBX2Vab9vFvIpnYCkvYv6EIZVzyDkeo8nYsnXDvVIYRChFsycr3OyATMN0MUW3bWYOgHxkoGx9PtAKtS/efW1Ql1nmbU5Pq+Zmb3XzNBHzJlIsGFPFzT6cjEkzDYFUzkEp6I2V8sz9ZivdFtQklChanraXeSFWLuZ+RHsiM05vLNx9Yt5plnkHxLeImDeRYl32Km52ySUWwHZMi0VQqpS66Fvm3RbUJBhyjotntSHsi3m0wjKPbwEgmpshGY+05V4yvpBt+T59Ki3zXNvFvP3pcSLmTWRqsbzJSpDh3ihL+danLwidRSrrlKWlQWnOXAIku5OgUIXzRsxdWzDGxy8gE8kvkIy33jLXWnNiIoVutS/ZUjlnnnXcth1L3vGKg692ImLeRKbsnGe1aPbhpAQ5CSaaPThfDpCMhY4cXNMAACAASURBVImGVVd/N1xP845P38sPT860+1BaTtBlG8mb7opeolzM/Xao0cJCW3qaTy7mrDHSHhWrtMw9r32u7k4xyETMm4hvmQ9XDYAzy6QKXHezmHXKItkBlFJdXzhmfCHL1x4a447HJtp9KC0nKErRghFzHa9lmc+TjEVI5ZyWNm46PbMELJ+7bhXV9tuuILhOqQchYt5EplI5hnujRMPLT/Nw0ZXavTdswUSzV4o5mCqB3fzduDhv5mMn7YC4mwiKQ8yKOYmhsm2cqElVi+YXWt45bW4pX2wE1C4hq7bfdgXBZUXMn/hUK+XqM5z0g5y615Uq2AC4KmI+3Bvr6u/GxfkMUMoI6SaCbtuEY3qZu7Us88I8yVhrO6edsVY5tEfMTV325fPj7SrpKpZ5FzCVytcWc3GzC0A65y6bMwcY7I12ZY9qn7EutsyDBUh6XCPmlalpftOVcsu8+feSpbxTdk3aIWS1XPvtssxlzrwLqFX9DST9SDAsZgv0J5aLeV88UmzC0o1cmOteMQ+KVcL1e5mXW+ZupBdPRQi7WfqjZvtWWOZnZzJlxVDaMWdeawDRLsu83XXhfeoWc6VUWCl1v1LqK/b5pUqpu5VSx5RSn1dKxezyuH1+3K7fH3iPd9vljyqlXlTvMXUK06l81epvAPFImGQsLOlHXYzWmnTeJRkPL1uXjIeL1eG6kbEF382ew/Pak3LULoLikPSqW+YoVYxoHw6ZrmqtuJeMVeSWt8Myr2yy4pN32vM96YSCMdAYy/zXgKOB5+8DPqC1PgjMAm+1y98KzGqtDwAfsNuhlLoGeB1wLfBi4MNKqeV3t01GtuCSyjnF3uXVGOqNScnOLiZb8HA9TV88umxdMh4hnesM91078APgHE93XSW8vGuvu9YktW+ZDyzbzp83H7QtUpvt5ZvPFChUFmtph5u9xj6dNiV7PyHc7EqpPcDLgI/b5wp4HvBFu8mtwCvt41fY59j1z7fbvwL4nNY6p7U+CRwHbqrnuDqByRUKxviYzmki5t1KsZd5Fcu8LxYh73odE1zTai7OZYuBXZ3gas8WXN50y908dH6+qfsJFiAJuTliOOSJ4kUSy7b1Bb5Ppwip5kezT6eWX4d2CFmtFLxWpub5aK075jdar2X+QeC3AP/TjAJzWmvfP3gO2G0f7wbOAtj183b74vIqr9m0FEu51giAAzNvPiNu9q6lKOZV5sz9oKZudLU7rsfEYpbrdhvLsxPE/PT0Encem+LW/zjV1P0ExTFSMAOHdKiv6ra+6z3mLLak2cp0FS9iO4SsljfA82j5lEzO8VreUKUWGxZzpdRPAxNa63uDi6tsqldZt9JrKvf5NqXUYaXU4cnJyXUdb6tZqcmKz0iyuwuDdDvFjmmx6gFwQFcGwU2mcnganrJ3yD5vXw1wH9/qve3IeFMFLPjeUVv9bSnUX3XbYER7X6K59dnzjsdCFcvfcXXLBXQlC7zVAXmdMl8O9VnmzwBerpQ6BXwO417/IDCklPLvTnuAC/bxOWAvgF0/CMwEl1d5TRla649qrQ9prQ9t3bq1jkNvPn4u5q6hnprbDMuceVdT7GW+gmW+1ObWju3Aj2R/UgdZ5r6QzWcKfP/EVNP2U17K1Vjm2XAty7xUBa7ZJV1n0vmaFmirBXSlwZTTcsu8c36fGxZzrfW7tdZ7tNb7MQFs39JavwH4NvBqu9lbgC/Zx1+2z7Hrv6VNZfwvA6+z0e6XAgeBH270uDqF4xOLDPVGV5wzH+qNspB1cNpUElFoL+ninPlyMe+18+jdaJn7OeYHtvXREw13hJj7QhlS8G8PXGzafoJiHs4ZyzwbWR78BoH67PkF+pvsZp9K5XBcr6r132rrdCXLvNX30k5JS4Pm5Jm/C/gNpdRxzJz4LXb5LcCoXf4bwG8DaK0fBr4AHAG+DrxTa905w50NcnwixcFtfWWtLSsZSUrf6m7GF+pqRWP6unjO3K/+tmuwh6398Y4S8/905TZue3isaa72svfNzgGQi1R3sxeigSpwTW6DOpPO881HJvif//LQshro1c7F733pIX7j8z9qyrGsdO7FzV4nWuvvaK1/2j5+XGt9k9b6gNb6NVrrnF2etc8P2PWPB17/J1rry7XWV2qtv9aIY2onWmuOTaQ4sK36j9DHb3Up8+bdiS/m/VXE3J9H704xz9ITDTPQE2FLX4zJKlHUrcYX89fdtI+FrMP3jzfH1R502yor5vkqaWlQYZknokyn8k1pAzqfKZB3PE5Np8k5HscmFsvWVwqo1pqvPjjGA02K/F9JsKuVeW0m7SpUUw2pANcEptN55pYKHNhWfa7LZ9hWgZtJi2XejaTXYJl3q5t951ACpVRHWeb98QjPuWIr/YkIX2mSqz1o6amcEcNqOeZQimaP5OcZ6ImQKbikmxBj4aek+bEMj45ViHmFdXpuNsNUKsditvH3Nc+rXpfdR8RcaCjHxk2hh4Orirl0TutmUjkHpaA3Vr0CHHSnZX5hPsPOQZNX3SlivpAtMNATJRYJ8aJrd3DbkbGmBD8FhTFsA+CceC3L3LjZo4UF+m3hoakmnKvpdB7X08Xqb4+Oryzm952ZBUoBno1kNTd6q93s7aoHXw0R8yZwfNKI+aqWeVLc7N1MKufQF4tUjaso5pl3YTT72HyWHQMmC2RrX4LZpUJbCoIEWcgUGOwxgvmyJ+1kMevwvWONd7WX5ZnbADgvPlR120LUt8wX6O8x35dGT0kUXJOSNrmYw/U0W/vjnJ5ZYilfEuplYn7aiPlS3sVtcHT5at8Dvwrc6ek0v/elh5qaRuh5nVMwBkTMm8Lx8UWSsXDRuqjFSLFzmrjZu5FU1qmalgYQj4SIhFTXudlNwZgcu4ZKljmYPgftZD4g5s84sIWBRIR/e7CxrvaCW6r+BtCbHTPLe7ZX3T5omffZGItGW+azSyYl7fycCUp87hVb0Roes95HWJ6edd+ZueLjVIOt89XE03ezf+WBi3zyB6e56/Hphu4/SNZxO6ZgDIiYN4XjkykOrBLJDtATCxOPhMQy71LSeafqfDmAUsrWZ+8uMZ9MGQtwx2C5mLfb1T6fKTBgrd9YJMQLrtnBtx6ZaOg+KiOjB7PnzPL+vdU2L7PMB+xAo9GWuR+57mcYPOPyLURCqmzePCiwS3mHIxcX2D5grttCg+fNazVZKa03x3LCeke/+1jziotlOygtDUTMm8Kx8dUj2X2Ge2Nt6Wl+5MICZ21hG6E9LGadqjnmPn1d2GzFb7Cya9C62X0xb3MVuKBlDrB/tJe5Brv/c4H5V+XmGchP4GqFN7in6vZuNIlWISJOmoGoJqQab5n7gnVhLsuWvhh9iQgHtvXxyNhC6bgD5+CBc/O4nubZB01Rr0bPm682R+0XjTkxYcT8jiaKeeWxxLJTJFKlyuSPT6Za6oYXMW8w85kCE4u5VefLfYaTsZa3QdVa8wt/fw9/8m9HV99YaBrp3Mpi3o1tUC/aiOlOs8wXMk6ZmPuWcCPFqqyP+dIFQnhc0FvoTfRWf4EKFa3zmJNiONn4ND4/WtsEJZoB1pU7+jk3mylOAbmuLs6N+8Fvz77CiHmjp4lWc9sXHA+tNScm0/REwxybSBWnCBpNZb799d99Kzff9nKi2RlSWYf//oUf87E7H6/x6sYjYt5gjk+sLZLdZ7g32nI3+6npJcYWspydFcu8naRWEfPeWIR0vsvE3Lpz/XgTv4JiO8U873hkCm6FmJvrVq1e+UYJVhPrTZ0B4HxoO6FQ7em6Yq55YZ7RZIzJxcbeS7IFE8Q2Np9ll70mV23vRwOPjS93td93eo7LtiTZN2IGII1OTwu+3/4jH+GGb7+RkFO6jxU8zcRijlTO4dU3GI9Gs1ztwbS0kJOhf+4oESfN6PidxVz8my4dacq+qyFi3mB89856LPOZFou5HxTiuzSF9hCch61GXzzSdQFwY7ZgjC+c8Yh53E4x9wvGDATFPGEeN3JOODhn3mPFfDy0c8XXBOfNh3ubYJk7ninl6ulin4n9W5LEwqFl8+Zaa+4/M8v1+4bpt4GdjfRcOK5X1qtgz4n/w/DkDxmeuLtsG/8e/KJrd7BrMMEdjzZfzHvSJff66MXv8thEimhY8eQ9g03ZdzVEzBvMsYlFYpEQe0dquMYqMJZ5a93sd1sxn0nnO6roQTfhuB6Tizm2D9TOeOhKN/t8lp2DibLg0a398bZWgfPFPGiZ9/tinmnc9QlGhfviMBnbteJrihHt+XmGe2MNnTN3PU3B8bhg3dQ7bYZBNBzi8m1JHgmIec51OTOzxHQ6z9MuGSpmaTTSMg8ObJWbJ54xAYiDMw8Wl2tdyoM/sK2P51y5le8fn2pKamNwztz3pACMjt3JibE5rtzeTzyyvIZEsxAxbzDHJ1JctiVJeAXXWJDhXtMGtVVtBLXW3PX4DLGIufRinbcHv83njhXSF5NdGQCXKYqGz9a+9haOqWqZ9zRerPJVLPP5+O4VX1OsAldYYKg3ymQq17CSrv7g4kJFUCLAVTsGOD+XKX7+vOMV58uftm+46LlYbOBgNGjlJ5YuoGyn7IHZB8u2Oz6Roi8eYftAnOdcsZXFnMP9gXS5RlCZY96TOl18HMvPMTD3MNfuap1VDiLmDefYRIqD29cWyQ5GzD3d+BSOWpyZMfPlL7ja5K5ebFJwiLAyfmewHStY5t3qZvcLxvi0uwqc/9scbKWbPW3S0hZ7qkey+zjRkmU+2BMl73gNE1A/wOvCXIaR3hiJaMnKvGqHucf5VnDe8bjv9Bx98QhXbO8nHgkRDauGutmD7+WfH4CBmQcIJnyfmExx+dYkSil+8sAWwiHFHY81No2wMse8x0axu2Hz3X22eoBrd1Wv3NcsRMwbyFLe4fxchgNb1zZfDjCcNDeFVkW0+/Plr7rejPgviGXeFsZtacyV3ewmz7wZzTM6EdfTjC/mlhVbaruYV3Gz+1Z6U9zsWhfFIZ2snmPuU6rPvlB0/TfqXGXt4OLifLZYxMfnktFe4pEQ//rjiyYFy/W49/QsT907RDikUErRn4g21HMRfK+gmMdysySWzhefn5xKc7m9Bw8kojxt31DDU9Qqc8x7rWV+8ZKXA/Dc8I+5csfajbpGIGLeQB6fTKM1HNy+djEfanF99rsen2FLX4xnHtwCiGXeLoqW+Qpu9r54BMfTLa833S78kqHL3Oz9cdJ5t23xA0U3e6Ik5slYmJBqnGUerP4Wy04RcZeY1X3EksMrvs5vwhIpLBSDzho1b54tuHietlMf5d6SSCjELz3rMlI5hz/92iO872uP8MjYAk/bVyo92xePNMwy9zxdltkRDDgDGLDz5tmCy/hCjssDAcjPuWIrD51faOiAsDLH3J8WOX/ZaykQ4SmhEyTdhWovbRoi5g3k+Doj2aFU0nW2BYVjtNbc/fg0T790lEQ0zEgyJpZ5mxhbyBELh4rXvxrJmN9spTvmzS9UpKX5bOkzueZTbQqCm19abpn7lmejUtPKXexGqE7rbUWBBqrG4ZTaoM7TZ5vzNCpYMFfwmErnKLi6mJYW5Kl7h/iTV17Hi6/dwZ3HpvA0XH9JafDRn4g0rJxrOu+UlbpNWMt8KWmmIXwx9wfJl28Nivk2AO481jjrPJhjrrwCiaULaBSzyYPc411JGI/+83c2bH9rQcS8gRybWCQcUuwfTa75NaXOac13s5+dyXBhPsvNl5ncx52DiWJer9BaxuYzbBuIr5hDXGy20iXz5v6NeOfg8jlzaF+u+UK2QE80XAwa9RnoaZzlWS347azeVla7PxENE4+WH0OhGM2+UPy+NMwyd9xi29NdFZZ58JhefcMe/uw/P4nffenVPOvAluK6/kTjzk/l+/SkjVt9Yu9LgJKYX7TTVwe2le7B1+4aYDQZ484GNsYJZgEl0hcIaZdc7w6Oz+b5tvsUAPrP39Gw/a0FEfMGcnwixf7R3mU/+pUozpm3wDK/66SZL3/6ZaOAuWn6Fbe6na8/NMZn7j69+oYNYmwhu2LwG3RfT/OLRTFfHs0O7RPzylKuPgOJaMPc7MG0ND/N6bTeXubaj0dD9MbK6xI4gWj2eDRMOKQaZpln824pLW2VplFb+uL80rMvIxIu3fv6G3h+Kn8D/pz5+B4r5rMPgfa4OJ8hpGDfSEnMQyHFtbsHi4VcGkF5jrm5Xkt9l/DYeIrvelbMz91BmTuhyYiYN5BjE6l1udjB3LAjIdWSOfO7Hp9mJBkrVqfbNZQouja7FdfTvPdrj/D2T9/Le770cFlrx2YyvpBj+yo3yN4us8wvzmVIREPLhLNUn719Yl6tuM9AItqwALhg9beSm317mZs9HgnRGyvPWy4Eotk9TzOajDHVoCpwOcfj4nyW4d7oskFEJa6ncSpiO/obOGceDH4LF9LEcjN4oSiLw9eQ7dlBxEnTu/i4qVQ31LPMoNo73MO52cbd68pyzBeNEZBJ7uWx8UUyQ1eQ7dlBNDMJ4w/WeouGI2LeIFxPc2Z6icvWEckOZu5tqLc19dnvfnyGp186UizIsXOwh8Ws0zWWXyXzmQJvvfUe/uaOE9y4fxjH0/yowfmo1dDalMfcuaplbm7c3XJ9Hp9Ks2+kd1m3wZFkrClNRNZKTcu8J9IwyzMY5OhHsp/V20jGgmIeLnsOAcs8v4DWMNrXmCpwecfD9XRZTfZVX1Mp5onGpVaW55gbF3umdzeoEAsjTwJM8Zix+Sy7q0wJ7Bk2jXEaEV2/LMfcWuap5D4en0xzxY5+pnY+26w8dnvd+1srIuYN4uJ8BsfTXLLGym9BRpLRprvZz84scX4uw9MDtYL9dJNujGhP5Rz+84e/z/eOTfHHr7yOW37uRpSCu0/ONH3fCxmHTMFdMZIdgnPmT/wAOM/T3Ht6lqftWx69HQ4pRvvaVwVuvqLJik9DA+AKy8VhKrarLKYiHgnRU2mZF3uazwMwmow3JFAw67h4WldNS6tFZYew/kSUVANSKzN5t9inHEou9qwNfvPFvG/6QZPaWEXM946YZWdn6r/XVeaYB6dF8q7HFdv7md7xLLPy+Dfr3t9aETFvEGemTbH/fRsQc2OZN1fMv3/cBH/48+VQCjTqxoj2w6dmODGZ5oOveypvvPkSBhJRrt4xwD2nmi/mY2vIMQeKVlg3uNmPT6aYzxQ4tL96Y4p2VoFbyBTKqr/5mDnzBrnZ7Zx5yFkinp2iQIRMbFvZNvFoiGS8XMydaD8aRaSQAs9lJBlryHnKFlxm03nyjrd2y3yZmEdwPV1WT30jVFrTvphnKsS8d+rHuJ5mh+2lHmTvsLkvN6K5VGWOuR+w+MCS+e4e3NbHzPZncOHGd8PL/rLu/a0VEfMGccb2Bl9rTfYgza7PrrXmkz84zYFtfVwZqE7XzZb5sXGTRviMy0vRtzddOsL9Z+aaUsc5iC/mq1nmfgBcN3RO8wdRN+6vnlfdzsIxCyu42VM5p9j+sx58IfRd7BdD20n2lKctxiNheqJhQsG7tgrhRM1vOmpLuk41oKRrruAVPSF+AGIlfXOPEMuU0r1yFWLe16BmK5UV7RKVYj5sxHx44VGiOMsqCELpvnx2pn4xL8sx115RzA8vDLFjIEF/Ioob7WPyyW+H7dfUvb+1ImLeIM7MLBEJqVWjPqsx0uTOaT84Mc2Riwv84jMvLXPbbR9IoFR3WubHJhbZ0hdjOFm6Yd64f4RMweWh8/NN3ff4Gkq5Qnelph0+NcvW/nhNz1a7xNz1NIs5pyyq3Mdf1ohcal8Ie62L/SzbixXdfOKREEopEhXNO0oR7YsM9UQpuLpY6GajZAsuE/Z8+wGIQfYc+zQ33/ZynvL9dyz7DD7+8ady9R3L8rQ0K+Z9Rsyd2ADpvv1EdJ4r1Fm29i+v3TDcGyUZCzckCC6YYx7PTBD28uTio5xYUOwZXpsXoxmImDeIMzNL7BnuKUvNWCsjyRiz6eY1W/nYnY+zpS/GK68vb9oQDYfY1h9vu2XueZo3fPyulqaGHZtIcXBbebnFGy81VmGzXe1+CtZqbvZYJEQsHCLVBXPmh0/PcOiS4WXBbz5+57RWl7atVsrVp1jStc6gKsf1ita9b5mfcrcxEIhkV8qIOZSyHHxKbVDnGew1x1TvvHm24DG1mCOsFCPJcnG85OhHuer+PwSgf+4oynPsa8q/p34kfr1TEbXc7Nne0v3Md7XfHD9NPBJe9j1RSrF3pLchlnm2rFuauWctJfcytZirmY/fCkTMG8TZmaUNudjBBK04nm5Ks5XjE4t8+9FJ3nTz/rJGCT47B3va3jntm49M8P3j08V5/Wajteb4eGpZ2d1t/Qn2j/byw5OzTd3/2EKW0WRsTfUIuqEN6th8lrMzmZrz5WDymBthca6Xau1PfXyxrfeYqvUxP+FsLU6zgBnY+QOdyvS0UhW4heIxTdaZnpZ1jGU+2hcrVZ7Tmsse/CAHH/wLNAo3HCfkFUjYAi6VlvlAA9zsjuuVBQeidWDOvFS33hfzG2MnASi4ywd9e4Z7G2KZl+WYp/xgxd1oWHOwYDMQMW8QZ2aWNhT8BiadBJpTrvKW750kHgnxxpv3VV3fCbnmt3zvcYBitalmM7aQZTHnFPPtg9y4f4TDp2ea2pJ2fCG7qlXu4zdbeSJz+PTK8+UAW+xvZLoFxZWCVOuY5tPfoM5ptUu5lvYZC3j8luWa24j2SGGe/rhttlK3Ze4ymcqVzZfvf+Rvuezoh/FUmIef/ufMjT7NHE/KCGiuwjLvs8dSTzpYZQ30SGGBSCGFE+mlEC99XxaGrgXgCoy4OlWKtewd6eHs7FLd3p1gAJzf+vSC2gmUt4ltNSLmDWAhW2B2qbBhMd9arD3d2BvVVCrHP953nv/nhj2M1ghi8avAtasz18MX5rnr8RnikVDLSsv6wW/VWtXedOkIc0sFjk+mmrb/sfnsmmMruqEN6uFTs/TGwlyzs3bLyNGk+f5ON/g3shrVepn7lHqa13d9yku5lgrGDFSUcvVZVgUuWrLMkw1otuLnUU8u5orz5SEnwyWPfAyAh25+P2OXvIKl/kvN8SyeAqrNmdsKhnWcn2UNTYJpaYEpmemIifwf8czAsOAsv5/tHe5lKe8yU8eA0PN0RbU+c71OuNsIhxTbqkTStwoR8wZQT1oaUBTaRt+oPvWD0+Qdj7c+89Ka2+wcTJApuC13X/rc8r2T9MbCvP6mfUws5poeSQ5mvhyoapnfZPPwm5lvPr6QXbX6m09vLPyEj2a/59QMT907tGK8iT9vO93iXPOV3ex+G9R6LXMrDp5Ljy2IUlmXPTglU9Myz88X+4jXY5nnHI9U1mEp7xbFfPu5rxMtLDI/8uRiPfSl/v3meKyYuxXFVPob4GbP5ivSwCoi2X3OFcyAZsCZAc+lUNUy99PTNm40LO9jbizzo9lRdgwkiITaJ6ki5g3gbB1paVByITbSze55mk/fdZqfunpbWQehSvyAjVa5uINMLGT51x9f4DU37OGqHf1oXerz3UyOTywykoxV9VbsG+llW3+ce5ok5jnHZTqdXzWS3ScZj7QsAO7CXKYsUrcVpHIORy8urDhfDu1zs68o5sUAuPoGW75Fm8iMEfIKLEa3kiNW5maPB8Q8EQ0TDpes0uCcueNpUzimDss8W3CLmQO+mO9+/POAafHpUxLzk4HPUvr+JGMRlKrPzZ51yr+PiZQV897yYN7xtMeUHiCESyw3VVZkxqdUOGbjQXBlXgati9kH96WHq1aeayUi5g3AzzHfN7oxMR/qNeUqG2l1zCzlmU7nedbBrStu57t729E97VN3ncbxND//jEuLVZtaMah4bLx2DX2lFDdeOsI9p2aaMvUwsWCu8VrFvC8eYakFbvZsweVFH/wuH/z3x5q+ryD3n5nF0yvPlwPFFMJ2udmrzpnHjVjVa5mXcsyNlTcV2wVQ5maPVwSv9gae+/XZI4V5Cq4uRv5vFD/4DYyYJ+ePMzR1H04kyfjelxW3W+ozHr9k6lTptYH55FBI0ReL1DXYqRxc+jEF2QrLfGIxx7g236F4ZqKqh29PAwrHBKe8ovlZIoUUhUgfJ9Lxtga/gYh5Qzgzs8Rwb7RqLupaCIdM+sdkA29Uvmhsq5IjGqRombc4oj1bcPnM3Wd4/lXb2b8lWeyX3OxBhdaaY+OLXLG9trfipv0jXJzPNrQxg89aC8b4tCoA7s5jUyxmHe5qQTnbIPecmiWk4PoqZVyDRMOmAct0urVu9oWMQywcIhFdfqssiVVj3OxJa+FeDO9GQUVd9vL9B+fNnZgtGpNfwHE9tvTF6vLyZQte8fVb++LsOvkFAMb2/TRutNSNLJvcjaciJJYuEnIyZZ/Fp9767LXmzP0cc5+JxRwzIePdSWTGq4p5XzzCcG+0rpKuwdLKfiT7XM9eQLU1LQ3qEHOl1F6l1LeVUkeVUg8rpX7NLh9RSt2ulDpm/w/b5Uop9SGl1HGl1ANKqacF3ustdvtjSqm31P+xWks9kew+W/riDbXMi9WbVhHzLX1xIiHV8lzzL//4AjPpfHE+v1WW+eRijoWssyzHPMiN1uXbjHxzv2f3WsW8VQFwtz08BsDRCwvLbsjN5PCpGa7eOVCWhlWL0b5YW9zsAz2RmvnvAz31d07L22Ct3gWT1XEmtIdkPLKsLnuQ3kBZ11I0u99spb4CO37BmP5EhN6Qw85T/wLA+ct+tmw7HYqQ6TPpYX6+dWWp0/5EtD43e4WYF5usLLPMsyxGjRcynhnHqZGNsnekl3MNssz9muwTYRvJvlnFHHCA39RaXw3cDLxTKXUN8NvAN7XWB4Fv2ucALwEO2r+3AR8BI/7Ae4CnAzcB7/EHAJuFenLMfRp9o5qwFuC2/pVFIxxSbB9ItDzX/PCpGbb0xbj5MiOcffEI/YlI0y3zlYLffK7a0U88EuLoxYWG739sjQVjfJLxMOm829RsA8f1+PejSpzhtAAAIABJREFU44wmY+Rdj4cvNP5zV6Pgetx/Zq44eFqNLcnGDnjXQq267D79iUhDOnFByTI/qXeWudjBlHINEgyCC0azg6l2Np3aeBEqf858W3+credvJ5afY2HoGhZHrlu27fKI9uWW+UYD4AquVz73rTU9aV/M95ZtO7mYI9OzHTBu9so68T5768g19zxNphAQc3u9zrCdSEixrUbGUKvYsJhrrS9qre+zjxeBo8Bu4BXArXazW4FX2sevAD6pDXcBQ0qpncCLgNu11jNa61ngduDFGz2uVuO4HudmM3Vb5o3qduSzVsscbK55iy1zvwJb0OLZNdjTdMv8sfFFoHpamk8opNg93MP5JpyTsYUsPdHwspt1LZJx06yiMu2nkdx7epbZpQK//LwDANzfgjawAI+OLZIpuNxwydrG7iPJWFvmzKvNl/sM9EQbVuzJF4fH3B1lkeyhEMsKDAXd7MFodoDh3hiOp5nb4Fx+tlBKS9t9wg98+9mq2y717S879krLvK8OMa+0ymPZKcJulnxsCDfaV7bdQtbB7fXFvLZlvmekh/OzmQ0NdJYKLsEg+f65owA86Oxh52CizJPSDhoyZ66U2g9cD9wNbNdaXwQj+IDf+mc3cDbwsnN2Wa3lm4KL81nT+nSDwW8+xs3e2Dnz/nhkWcvEarS6ClytCmw7hxItscyHeqPF6OhaNKpaVCVjCybHvJbbthLf/dxMV/ttR8aJRUK85tBedg0muP9Mcyvg+ZywufxXrDCwCjLaF6srR3gjrCrmifrd7GDyuHuWzuOpCI/lRysi2Zf/hssscz+avWAs83pKumqtWcwWmF3Kc3VsipHJu3HDPYxd8jNVty9GtNsguOWWeXTD393aOebl8lCcUhgwgYPxzAROjRTXvcO95F2P8cX13+8qY1f6Zx8C4PtL+9ruYocGiLlSqg/4R+DXtdYr+eeq3b30Csur7ettSqnDSqnDk5OT1TZpOfWmpfmM9sVI5Zxlo9GNMpnKrckqByOiY/PZplY9CzK+kKtaga0Vg4rj4ykObutbVUz3DPc0RczH59de/Q2a3wZVa81tR8Z45oEt9MUjXL9vmB+dbY1lvt76DKN9cWaW8g3pUrZWFrKFFQNbBxL1B8BBSQwzffuYz1VEslcp+xsNh4ja5U6xNvsCaI8hO/jYyLx5KucwsZBDAzd4D5j32fW8Mks4yDI3+7I5841PQ1TmmFd2S/PxI++jQ76Yj1ct5wrB7mnr/20HByWxzCSJzASFSJIHloY3v5grpaIYIf+M1vqf7OJx6z7H/p+wy88BwYmOPcCFFZYvQ2v9Ua31Ia31oa1bV065ahXFtLQ6xbxUBa4xrvbJhbWL+a7BHvKu17LgomMTxtV9oCIIbddggpl0vmEDmkq01jw2sbiii91nz3APM+l8w0X04nx2zcFvQLF/dbMs80fGFjk7k+GF1xgX5fX7hjg3m2FiA5bLejk9s8T2gfiavEcAo8kYWsNsEzsMVrImN3sDCi75bup0/6Wkc05ZQGA1yxwgac+bDkVwIkkUmkghVSzWspF7yXymULz2lzinAFgYXj5X7pOuyDV3PV0WSd6f2HhqWqVlPjj9I7PPgQNly/1BS8+oEfl4jWh2gL3DG881D+aY9889bPbddxWaUDEbp53UE82ugFuAo1rr9wdWfRnwI9LfAnwpsPzNNqr9ZmDeuuG/AbxQKTVsA99eaJdtCkqtT+sbmfn12Rvlal+XZd7iXPNSOdVKN3uPPY7mCMlUKs/cUmHF4DcfPye1kda552kmFtdpmdub+lKTirnc9vA4SsHzrzZi/tS9QwD8qAXz5meml7hkJLn6hhb/N9IqV7vn6Zq9zH0GbOpVvV6t5IIRw9meS9BQ7mavkhYHlA2CCjFz3aL5uaInwU9PXQ9zS4ViWekdWRNdnxq6oub2+cQ2nEgvsfwc0ZyZngkOxvvjEfKOt6EMicpB/ej49wCY3vHMsuUTi1n64hEifVvxQlGihQW8fLrqe+4e7kGpjeWaBwf2AzNGzE9FzcBis1vmzwDeBDxPKfUj+/dS4L3AC5RSx4AX2OcAXwUeB44DHwP+G4DWegb4I+Ae+/eHdtmmwG99Gq4z+GG0wZb5xEJ21Uh2n912tHq+CW7lahybSDHcG2W0orWiX3ShWWlyvkdgpbQ0H78vcT1pLJXMLOUpuHpdPe+TTZ4zv+3IGE/bN1wc+F23e5BISLXE1X5qOr2uQkt+ffZmNCSqRirv4OnqBWN8BnqieJq6S+72LhrhnIqbhkj9q7jZoTwPPdu7A4BE+jzRcIieaHhDg+IFa5nHIoqhRVNAKDVYW8xRqhgE5xe9yZWVdN14z/egZZ5InyO5eJJCtJ+FkaeUbedH3qMUuYQJ0YqmJ6pOx8QjYbb3J9btZjeR7IFBirXMj3AZ0bAqa0jTLtYWUlsFrfX3qD7fDfD8Kttr4J013usTwCc2eizt5MzMEvtG125d1GJLAy3zdM4hHairvBrNsEJX4vjE4rJIdih1HGpGFLnZb3WPQDVKYt64Y1lvWhqUAuCaMWd+bnaJhy8s8O6XXFVcloiGuWbXQNMj2jN5k8t8yTqmpxrtvVqN+SW/yUrt22SxPnvWKbOm10uxYEzUzDiWi3l1N3uw+Uo2uQem7qUnfY6Cp202xvoGogXXYylv0tKuSaaI5hbIx4bIJ7at+Lql/v0MzB0huXiShdGnllvmgfrstZo91SIonqNjxiqf3fYT6FD59ZhM5Yolq3M92+lZOl90tYdDy8/d3pGedQ/S03mnrCb7wKwR88O5fewc7Gl7JDtIBbi6MQVj6nexFK2OBlS48ueQVqv+5jPYE2UgEamrzOFa0VqbcqpVBHVH0d3fHDf7sfEUA4nIms7L1r448UiooZb5+Dqrv0HJMm+0mE8u5nj/7cbyeuG1O8rWXb93iAfOzTU10MyPNblkyzrc7MnWutlXan/q4wt9XfPmWhfnnM+FzLxvWfvTGpZ5sCqdHxTWkz5HwfXYPbT+1Eq/dO1kKsfT4iZsKTV4ZVl3smqs1D3N/xzrTU8ruB5uIIhtdOxOYLmL3bGxPr7hkvNzzbPVS7rCxnLNy8q4ZmdILF3EifRyz+Jw28u4+oiY18F8psBcHa1Pg/TEwiRjYaYW679RrSfH3KdZqViVTKXyzGeqz1snomFGk7Gmzd0/Nr7IgTVEsoOp0b67wRHt/iBlrXXZAfpivpu9MXPmJyZTvPufHuAZ7/sW/3z/eV5/0z4urRDUp+4bIp13i9MSzeD0tJnTXI9l3oweBiuxUvtTn/4GdE6LZiaIOGnysSEmXHMt+gMBcNVKyZrlJavTF/NE+jx51zOW+Tq/u/OZAlprphbzXBM2keMrutgtld3Tgpa571lab0R70CpXXoGRiR8AML3jWWXbTaXzaF0yXHI9xosQXxqv2mwFjNft4nxmXR0ag4Np38U+P3gV0xmvrT3Mg2zYzS6UIiIbIeYAW/rjDak9XazLvo7eunuGezg5VT1opJGUItmru7p3DiWaUjjG8zRHLizwiut3rfk1jR7gnJ5OE4+E1jXI8qPZG2GZj81nedmH7sTT8Oob9vCLz7yUy6p01Lt+rynicv+ZOa7aUbvHeD0ULfN1zJmHQ4rh3hhTrbLMV2iy4jOwQcszSM/8CQCW+i9jMVtAQVk0e6xGa9h4JEQoBJ5XajzSkz6H42p2D/Uwu1RgKe8s639ei/lMgflMgbzrcbk2899rEnO/cEzKeBdy1dqgrvP7mw0EfA5M/5hIIUW6/9KaOeaVlnkiM161DSrAnpFePG26BF6yxinS4GB6YMbkl59PmHPTCcFvIJZ5XTQqx9xnNFlfgwSfSZtasp6gDFOzONPUsqEQmLeuEYRmcs0bb5mfmVliMefwpN2Da36NyTVvnJv9sfEUl2/tW1ewZCQcIh4JNUTM7z45Tbbg8fm33cyfvupJVYUcjMAO90abWjzm1HSagUSEod6Vi/dUMtoXY6ZVc+ZrEXPfzV5HrnmPrcmeHriMxaxTVpc9HFY1+7wrpYrz6ZVudj/mY62VHbXWNi3N3H/25I0wp4auXPW1Jcv8NGivzDLf6GCn2nx5pVUOy8tWZ62Yx7Ljy/LUffbaGKEz60hPK7fMjwBwf2E/YaXWlB3TCkTM6+B0oy3zBlWBm1jMEbFWzFrZM9xDpuA2Pdf8+ESK/niE7TW8BrsGE1xsgmX+4HlT6vK6dYr57FKhYZHkq3Vrq0Wjmq3cf2aO3lh41QGNUoqn7h1qahDc6emlNVtFQUaTjfFerQW/sttKbvaBBrjZfTFf6r+Uxayzpkh2H9/VnuvZjqfCxLMTKCdXDLJcq2cpnXdxXc1kKkcEh5GMzXuvyOmuhhMbJB8fIexmiGcmqlvm6xzsBMvCjo5Xny8HM6UYj4SKRXby1s2eyIzXjD/yvYJ+iuxquJ4ua8X6f9s77/g26zuPv3/akiVPOZ6JEychk2zIZAXKammhQOG6xxXoddHeXTl617ve9a7X3R6lPbpuAFegLS1QGlYzgEAI2cRJvOLEM96yLcvW/t0f0iNLtmRrOc543q9XXlHkR3rkXx493993fb5K8dsrAyUsLLFG6lpmGtWYZ8C+U/3MLjRnVMUaTZHVGOnxzIQepwe71ZhSheXss1TR3hAufkuUty7LN+P0+LM2vEKhpn0Qg1aTVFuaglLln42WPafbR8egOynBmvHkGHVZ6TM/2OJgRWVeQk8vmtVzCmjsGc6a7vh4WvpH0pJALrSePX32wVEfGjFWtxAPxVhlMrPbFO4xH7HNw+nxJVXJHnltOJ8uNTrcltD0LpOrnWJbaCOfbBFcpPjN6aFa04lW+hnJmZ1Q+W080RrtgcCYcIw1qpo9FRTPXO/pJ7e/hqBGj6P48gnHdYfFsZT7ScQzH+1mIIFaYLHNiN1qoLYzuYFC0RtpnWcglMrQmNgzZI/oMpwLqMY8TUa9AXY39nLt4pKsvafdaqDf5cm4irjb6UkpXw6hAQSQnjJSKoQGrCS+QZRNU0V7Tccgi8tsCSuD45HNXnMlvZCsDnk0FoM2Y8/c7QtwvGOINVPMDVdYN7cAKeHtpuxLPvgCQdodo2kZc3vO2RuDOhiemDbZplin1WAxaDPzzMM5c1fEM4/WZU/OM4fYvHmBOTTaONmNqNKG1zXkZm2kkn3qfLnCyDglOMU714dnwad6/SqecGHXmwgkDvs6grqJ10v3sCemOyVSADfaTTAgE6oFLiq1UdeZXIFnjFhMuPit3biAAFpWVqrG/LxnT1MvHn+QrYsn78FMBbvVSFDCQIZylT1OT8oiBmej19zh8tI77JnUO67ITy3XlwxSSmrah1IKsUN2oxVKSC/dMHumOfN32gbxB2XSxnxtVQEmvYbdjb0ZnTceHQOjoeFEKai/KRTmGENFWtM4RU5hKilXhVxTBpPTfKMYXaEBK6PWOSFjnkQlu4I5QUV7QEpK80wpeeaBcJHo+pwzQLgtLUmUcLx1MNTuGNtrnvpMc3dYMU7Jl/fHyZcHg5JeZ6zSZVBnwae3oQ160XsdCaM4i0tzqetyJuU4xeTLHaF8+ZFAFeX5ppSKWacb1ZinyfYT3eQYtKyvTm4WczIoohiZhtrT8cytRh0FFv209po3hqdkxesxV5gOSdfW/lEGR30sL0/NmNuthqz1mtd3OTHqNJFNUyrkZMGYK8Vsq+Yk50kYdVoun1c0Lca8WRmwkoZnrnxHzoY+u2PEm5wxN6c/5pP+JgSSUetsBj2hkG5BlDLi1GH2+J650muezKbYFwji8vhp7B7G5Q2wXBduS5tExnU8zvyQ8JBtoBaYWNGeShrC6w/3mEtJYUTCdfOE4xwjXvxBOUHp0hM11zxRG+PiUhtuXzDSIjkZjpGxjYgyKW23q4JV55BXDqoxTwspJTtqu7liYfGUX7ZUUIRjMumjDQQl/a7UPXMYq2ifLiKa7JOE2UtsRjQiu555TUeo+C2VSnYIFYJla3paffcwC2alVsmukI0CuIMtDqqKLNhTuC6uWGCnsXs4690FzWm0pSlkUylxKuo6nSxIUPEfTUaeeW/Ikx2xzeNQWEJ3RdR1Gm2s4zGpcEySveZKiuBw6wA6jaBcqWRPIcw+HDbm1sG6CRXtNmNqmx0lX252tWIa7cZrLGA4b/GE4xQ9jfEiUGOh9i5GvIGY4jUFpeVyqlD74IgvJoWiFL8dDc5j5TmULwfVmKfFiTNOzgy62bokeyF2IFK0kkkfbd+wh6CE4hSESRQqC8y0TWPOvKHbicWgnVRkQafVMMuW3V7zo+2D6LWCS0pTD3Fnq9c8VMmeer4cQr3mrgxEY6SUHGwZSDrErrB5gR2ANxr70j53PFrC/fYlSc4OiKZQ2fBOc0V795CbbqeHZUlsAEOT09LcbPU2AqF8+cFmB7NsxkitBkC+ZfLIgEkXHWYP9WCbRsLtaflmOofcU4qjOEZCYjGHWh2sKdFhGekgoDEwGi5qSwavyY7HZEfnd2F2tcWMQrWZ9AynsNnxhI25MiVtsGh1XBW6bmd8cawxz7wLiK/lv7DEikbAiSmMeXT7ms4zgGW4BS96Oo1VE8SWZhrVmKfB9hOhi+SaRdk15hFJ1zTmECtELvB0PPMCC20Do9M217wx7J1OVWVflm/KqjdY0z7IJSW2tKIo2eg1H3L7ODPoTkoTPh6ZhtnbHKP0OD2sSTLErrC41IbdamB3Q0/a547H6b4R5hRa0tKzPlv67Eo0Z3n51KI5tkxmmvc1AOAwVXGic4i1VQWRymyLUYt+is4DjUZEpqpFe+Zev6Q830xQjs0EiIfbF6C1f4T2gVF6h71cZw9t3EZy50/QQJ8KZ/4SAKwDJyI5b1Bmmqfumef1HQJgsGhV3OO6hzxoNYLCcS240WF2IG7BpEmvZa49h9oziSva3b5AzCjgwu43ATgkF7Ks0o4mCSXJs4lqzNNge203K2fnZ734Ic+sR6cRGXkdEV32FHPmEDJcXn8wEr7KNg1dwwmV36IpzzNnLWceKn4bTDnErlBZYMm41zxS/JZCW1w0VqMuPOghvU2WEr5dnaJnrtEINs23s7uxL6tiQi196bWlAdgjnvk0G/P20E0+Kc/clMFM896QMT80WkxQhgoPFfLNyelEKKF4r6mYgMaAweMg4HaOTUOcJGVVHy4CU6bkrTO1A+BMofhNQQmF2wZqx3nmqRlz5dgxzzy+MW/qHaayYOKQE2VymuKZO0a8cR2UJeEiuES09I/EDFexn3kVgO3+lQlb0vQpdMtkG9WYp0iP08ORtgGuzWIVu4JGIyjMMWSkz96TgWdeWahUb2c/1D4w4qVzyJ2UMS/LM9ExkB01uvaBURwjvpQr2RUqszAetiF8w0g/zK4jKGOFNFLhYLMDs17L4tLUz79loZ3eYc+kN71UkFKGhxOlF6LMNetCG95p1mevaR+k2p4TI6s62WdyutPcbOUU4zUVsb0nD7vVEKNVP1WIXSESaheaiNzpSHdTxHNMdO32Dnsi0s+HWweotudQPBJqk0ulkl1BKYKzDtTGeOZWY/LV7G2OEToH3Wj8o1gHapFCw1DBpROO8wWCNPW44kabPJbYMHsgIBmIs9laXGqjuW8kbtTLHwjGboJkMDLsZTerWVI28buUY9SxKM3veDZQjXmK7KrrRkqy2pIWTZE1M4UrJSyUTtRgdoHSa579IrgXazoB2DzfPuWxZflmPP5gTBVputS0p1f8ppCNXvP6rmHMem1MLjQVMp1pfigFsZjxbAnnzXc3ZKeqvcfpYdQXSNszFyK04Z3uMPuxjqGkvHIIeeb+cfOuk+ZDv+G19+7hrS7BmjkFMWJKyRpzs2FiEZxxuC0SLYhX8xEMSurD+eJ+l5fTfSOsmp0faS1LpZJdQSmCsw3WEQjIyHQ7m0kXUpibIn3X7XRHCtJyHTVoZIDhvEUE9BM3fqd6XfiDkq1xdD7Gh9khflHxovDmtj7ORrVjwB0ztc3mOI7R3csZitCXLo2bsrNbDZTmmSjIyY6IWKqoxjxFdtR2U5JrZFkSubR0sFsNGbWm9Tg95Jp0U1bBxmOs1zz7nvnTB9tYMMvKisqpb5CKd5KsqMNkHG0fRKcRkS9uqmSj/76h25lUrUAirBkMW3H7AhzrGGJNVWohdoXyfDPVxTlZa1HLpJJdIbThnT5j3u/y0j4wmlS+HMbkXtMtgjvY4iAQlDEhdoNOk/SAlGjDEmlPG25Drw3JnB5uHcDhig01n+5zRVQFjygh9jJ9RHc8Hc98xDaPgMaA2dWG1uvkcKuDHqcnomg32WZ0YMTLsfahSFg7r3fyfHl9lxMBXLHQHtPKBxPD7BBKy4yPnCwpC/3/1o67z0gpJ7To2jtfA2CnfyU3Li+L+5mUee2LSnOnmho7LajGPAXcvgCv1fewdXFJUmM008FuNWY0bCXUY57efF2TXovdasy6Z97c52LfaQfvX1OR1LptmF+EQathR23XlMdORU37EAtLbGltbiA7veb1Xc60i9+AyE09Hc/8aHtqYjHxuGKBnb1N/Xj8mUvKng5P5ktHl12hKMcwrdXsqUZzIvrsaRbBvX2qnwKLPqY6Opn+doWYUaiWsYp2CBmYriE3B5od7KrvZm9TH8c7hjgd1V99uHWAEpuRVQOvoPOP4LCvi2icp4LU6HDlLQTANlhLMAhH2wciBjpeqN3tC9A15OZw60CM5z5V8Vtdl5P5s6zkWwyRiKKC12RHCg1GTx8iGDrnsNvP9hPd7KzrZndDL/tO96PXCnIM2pgiOK8/SEP38IR2toL2XQCcKtgUV5JZqxXkh//PrEZd1oZvpYJqzFPgpzsbcXkD3Loq+TGaqVKUYQgxHfW3aGYXmmkbyK5n/vuD7QgBt62umPpgQl+G9dWFbD/RPfXBkzBW/JZ+FCXTXvPBUR9dQ5608+UwNg4zHc9cEYtZnWIlezSbF9gZ9QWyMnilpX8EjRhT+kuHomnWZ1cq2ZclKTIU0WdPowhu2OPnSNsAa6sKYqqjkw2xA5gNUcbcOlbRDrEbn2AwVFzWMTCKMh10xOuntsvJqso8KpueAqBt/l+k/HtEfp+8sby5cs7+8PlP9bo43euirtPJoRYHr9X3sLuhl6Ntg7Gzx6Ukr/8IAANFqyecwx8McrLHxfp5IcGuYpsxUtEPoU2F1xhKDxlGYzsxAgGJ2xdgcMRHTfsQFQVmjrYPMuL1c+LMEG809tLSF3v/03kGyHe8g1dqmbPupri/d6HFEBN5q7bnxHyms4FqzJOkocvJf756kttWV7C+umjazmO3GRn1BdJuRUpH/S2abM/wllLy+0NtbJ5vp2yS/vLxXLekhKZeF009yU02iseZQTd9Lm/axW8KlQWWtJXxlOK3TMYkKjlzlzf1a+LtUw7mFKYmFjOeDfOL0GpEVvLmzX0jlOebU9LIH09RjjGSj50OjrUPMbvQTF6SBjUSZk/DM99R240vIFk7LnKSbCU7gClqLaNV4AAKcwz0xwkxKxxoDoX4byjoIHfgOF5DPt2VN6T6a0QYrwQHY5GDmrZBGruHae0foW/Ym1CS1+Rqw+juxWvIZ9RaNeHnzX0jeP3BSD2HEGLC5tAdNT1tMsryzNR2OnmjoZd2x2jcvL6pZRcagtQbL6WsuDju+9jH1SjptJqz3oeuGvMkCAYlX/3DUXKMOv7+3Uum9VxFOen30UopM/bMKwtCEpCZDntR2N/soLV/lPevSc4rV7g2LMiTiXeutNtkbszT98zru9IfsKKg5MyHUxSO6R5ys6uum3ctzWwYUK5Jz5o5+bx8vDOjDoNAUHKsYzCjfDmEPPNhjz9GZSyb1HSk1sqYm+ZkMIAn326hMMfA/KjNnkZDzOS0qdBpNei0Ia8wWp8dQvcTX0DGlVNVlCwr8s2s738WgI55tyO1qc2Yj2ZMCW7MmCuRg2QLBGNa0uKk5ZSCtcvmjUlpl+eb0URZM6UIzhTe1CSiIt/MiDcwabFtoO5lAHzV1yU8pihn4pplUx00GVRjngS/2d/KvtMOvnrTkow8nGRQ3j/RLN7JGPb4GfUFMvLMZxdY8AUkXUPZ6fN++kAbFoOWG5aVpvS6ygILi0tt/PlE+nnz/33zNKW5ppQ12eN9loERX1reYH2XE7Nem1FY2WoMeX6OFM//67db8AclH9kw0btJlfevqaS+azjSs54O/7X7FCd7XNy5dnZGnyWy4Z0G73xw1Edz30jSIXaILoBLzTM/0NzPmyf7uHl5WUyIPdc0+aS2eCjer89QgF9nQe9zovMORoqy4l27J3tctDpGuXmBmdLWPwHQXn1XUufLMepYUj6x0CvSnjbYAMGQ8bboQxuTZMf4xii/xaGuy0nVuGiTSa+l2DpWKzRgXwdARdNvJj1X5RS9+J0DIyxyvQ2Au2pr3GOsaRYcZxvVmE9Bj9PDN7edYP28Qu5cVznt54sY8zRU4HoSyBumQmVB9kahun0B/vTOGW5aXhYJFafCdUtK2N/siIxnTIW3T/Wz91Q/915VnVFIF+CqS0KhtSfebkn5tQ3doeK3dCvZIaQ9XW3P4bkjHUm/xusP8n97W7h6UTFzsxDuu2VlORaDlifTWAMIqf999+U63rW0hPdlWHMSMVDTkDc/pii/peCZpzvT/OEdjRRY9JEolEIq+XKFyPQ0IWJC7WMbn4n3kx213Zj1Wm7VvI424KavZBOj4VGmibCZdKyYncfG+UVU5JsnFOr5DXm4LWVoA24sw6cBKMgJiWE19yV3T5ms+C0YlDR2D8d45Qrl+WPGvL36A/j1Vgp79pIb3hzEQ9lkxytw7Rv28PL2F7GLIVzmMly58+O+hzIvYKZRjfkU/Mvzx3H7gvzbbZdOWwV7NBG5yjS8DkXKdfwUoVSYXZi9sZ8vH+/C6fFze4ohdoVrl8wiEJTsqk891P7wzkaKcgzcfdmctM4dzdLyXK5TqjiOAAAZZklEQVS6pJj/fuNUyqHdhq7hSUe+JoNGI/jQhioONDsixmYqXjrWSY/Tw8c2zs3o3ApWo45bVpTzxyNnUh5n6Q8E+evfHsFi0PJvty3P+HtUmKPMMMh+RfuxsPJbsm1pEAqnmvSalLpQatoH2VnXw19eUT3Bq8tLIV+uYIo7CrUtofztwIiXA80Otswvoup0yHttr7570nMUWQ1cPq8w5v5SFCdS6cyLzZsbdVoWl9k43DYwZZpG43djU8RiCieKxbQ4RnD7gpG5AdEUWAyRUHvAYKN1/gcBmFv784TnyzHqKLQYJtzvepwevvNSHev9+wFwlF8dN+QPYzLcM41qzCfhT++c4Y9HOvj81gVJKZdlA+XL1zNDnnl5vgkhyHgUqpSSx/c0U55nYkOaBYMrK/OxW428cjy1UPuR1gFeqw/dKKMrfTPhvqvm0zvs5XcHJs/BRXO0bZBup4elWdAkuGNNJSa9hsffak7q+Ef3nGZOoSUSVcgGd18+m1FfgD8eOZPS637x+imOtA7wL+9bntFGU2E6J6fVdAxSnmeKa6Qm47K5hbxyvCvpuQY/2dmIzaTjIxsnpkDS8cxjpqeF29PMrjYsBh1mvXaCc/B6Qy8BKbmr+DTWoQY8pmJ6Kq5N+P55Fj0rKvMnbMTieaXD+bEV7QCrKvPpcXqmlGnOdRxFI/0M5y4koJ94z1Xy5ZvmT7ynaDQCm2ls7VoXfoyAxsCs9j+TM9iY8JyVBWaOdQzxzKF2jnUM0uYY4bsv1THP18C92ucB6CmPH2LXakVa/1/TgWrME9DtdPMPzxxlZWUen7k6fnhlOjDqtCyYZWVHbere6Jhnnr4xN+q0lNhMGXvmu+p6ePt0P/ddPT/tELNGI9i6uJhX63umnPwUzcM7G8kz6/nwhsy9coUN1YWsnJ3Pz19rwp/EZ5FS8vU/HsNuNWQlPZNn0fO+lRU8c6iDwajcbDAo+f7LdTz2VnPEkBzrGGTfaQcf3ViVUXh/PKtm57OoxMaT+5IPtdd1OvnhK/XctLyUW1bEF9tIlbE8cPY986Ptg0krv0Vz+5pK2hyj7D3VP+Wx9V1OXqjp5OOb5kZ61BWSGa4SD3O8uebDrUC4oj1q4+MPBtlV38Pdxc1cdfjLQCgsLTXxjZLFqGXV7Py443ttJv2EFqxIRXtUEdyK8OzvI22Jay5E0EdJSyh3nyhfXt81THm+iZIEWhr5UWF/r7mYM/NuB6Cq7hcJz7t18SzsVgPbas7wwz838PU/Hic/0Mv/mn+ALuimY+7t9JVeGfe1RTmGsxKxTQbVmMdBSslXf1+Dyxvg+x9YmZYMZiZ8ZEMVh1sHIspMydLtdGPQalISnIhHZYE5o5x5ICj59ou1VBVZMg5zX7ukBKfbz74kbpIAJ84M8crxLj6xeW7MLj1ThBB85qr5tPSP8EJYmnYynj3cwYFmB1+5cfGEG3a6fGRjFaO+AE9HRQd+tL2BH+9o5GvP1HDHI2/S0OXksT3NmPSajAvNxiOE4O7LZ/NO22BS4f4Rr5/P/foguWYd37g18/C6Qo5Bi0GnybpnPuzxc6rXlVbB5A3LSrEadTx9cOrIzU93NmIxaPnE5nkTfpbuRtwYHWa3hv7fK5ueZMXuz3Ctvob+4bHN+aGWAa7y7OJfh/8RvW+I7orrOL343gTvq2HNnIJJNxjjw8zK9DSb4wQaf+i8hTkG5hRaIh0mMcggs1q3sfHFm5l98tcA9JZfM+GwoJQ0dDljlPLGM76dsHnRp5BCQ2nzHzG64tecLK/I4x/evZSH7l7N/dcu5K6Vhfw290fkeHpwFF/GibX/nDDEno1IU7ZQjXkcnj7Yzp9PdPGVGxaxIMN8Zzq8f00FOQYtj+5JLqQqpeSxt5r5nzdOs6jUlvFNc3GZjUOtA3E1i5PhmUPt1HY6+ZvrF2VcfHbFQjsGnYYf/bmBB39/lE8/up8P/fItXo8zllNKyY93NGA16vj4prkZnTce1y8tobo4h//cdXLS3N+wx883t51gZWUed6zJXtHk8oo8Vs/J5/G3mpFS8vw7HTy0vYE711byw7tWcqrXxc0Pvc7vD7Zz2+qKpPukU+G21RUYdBqe2tc65bH/+OwxGnuG+dFdq7PaBSKEoNhqpKnXNfXBKfD2qT6khEsrU0+LmA1a3n1pGduOnplUI6KmfZDnjnTw4Q1Vkdx/NMVpGofoMHtfyWbOVL0XKbTM6tjOvzi/xlOj91C97W7Mz9/H8re/wkOGn6CTPloWfox3Nv6YoC7+eecUWqas1B4fah+1zsFjLMLo7mbztndRfvIpRNDPyso8mnpcON0+jCOdzGp9gUsOfZMNL93Cij33YxluxmWbx5FNP45rzLcdPYPLG+DKhYlTR+MdmVHrHLpm34xG+lly8J+oqv0FVSd+RlXtLyjo3osIjG0ITToNl+U6+Bvndyly1jJincM7mx5O2Kqn04qMoqDZJvUS4wuY1v4RdtV1850X67h8bmHcnfPZwGbSc/vaSp7c18pXb148af7O4fLywNPv8PLxLq66pJjv3bky4/N/8dpLeOFoJ/c/eZhnPrs5JYPs9gX4wSv1XFqRx7svzTysajHouHl5Kc8d6aCp14XdasDp9vPx/97H19+7LNJ2NeT28eDTR9l2tJMvbF1AviX7FaYajeDeK6t54OmjvFrfw9UJ5tn/ZGcj3U4Pj3xkbVbD3AAf3VjFl546ws9ea+JHf65nbVUB/3rbcow6LVcuLOYbzx8Ph3Cn59rNtxi4aXkpfzjUzoM3LUlYk/C7A2387kAbX9i6gC0Lpx6ukyrvWVnGz15t4mjbIJcmofc/FSfODPGlp44wu9DMZXMnVkonwx3rKnlqfysv1nRy+9qJm7jBER+f+b8DlOSauO+qiak7s0GbdlTNqNOi1QgCQUlQZ+bY+u/RsPLvKG/6LUV1v6bM1wXDfVSHj5cI6lc9SOslH5/0fROFs6MpzDEgBGPjQoWGI1seYfHBfybXUcPSA1+jqu5XLNEXc7+hjcrnBzEGYyN/HtMsmpZ9jo55d0yYoy6l5OmD7bx4rJNN84smVZI06rRYDNqYNrjTiz9Nacvz2M+8GhljquDX5eCYtR6foYCC7j2YR0Leu0+fy+EtP8NnTBwFKMubOH51JhHZnFN8Nlm3bp3cv39/xu8TDEr+/YUT7KzrobE7JPCxcJaVX33sMuZkKG6RCY3dTq77wWv87Q2L+Ow1Cyb8PBCUPHOone+8VEu/y8sDNy7mk5vnZe3ievlYJ/c8doC/uno+X7lxcdKv++XrTfzrn07w679cz6Y4FafpIKUkKInk7IY9fr74xCG213bzic1zuW11BZ9/4hBtjlH+9oZF3HNF9bR9yTz+AFd/dxfdTg9bFth578pyrl9WglmvxeUJcLrPxZ2P7OGWleV8/wOZb6zG4/YF2PStHfS7vJTnmXj2c1smFDz6AsG08q7Jsrepj7t+/hZLynJ58KbFXDmuyK6hy8l7H36DlbPz+L+/3BA315opTrePa763i3n2HH5z78aMolFNPcN84Gd70Gk0/Pa+jWnrakspufp7uyjPM/PEPRtifhYMSu55bD+76nr4zX0bY7Ty95zsw+XxM9duySgSuO90f9w2TrfXS+/pGuboh5hFP0Z3LwP2NQwWr5v0/fItetYlubE50NyPwzXu3DJISesLzK/5IZbh2DoLn97GYNEqBotWM2hfzUDR2rjRgaCU/HpvC7vqe7jqkmK+/K5LWJlglrjCsY5BzgzEFtqVNj+HzXEcKTQgBFr/KAXde7EONcQc5zXk45i1kdOLP42zcPmk51lfXZjVVF4yCCEOSCnj/sdd9MYc4NafvIHNpOOaRbO4ZvGssy7Dl4gP/fItTvW4eO0r10Ty9lJKXjnexfderqO+a5jlFbl86/0rMlY5i8cDv3uH3x5o5al7NyblrfS7vGz9/i5WVObz6Ccvz/rniSYQlHxz2wl+tfsUAKW5Jn78wdVpe1Wp0No/whNvt/Ds4Q7aB0ZjvRJCOd2df3N12gNvpuKh7Q387NWT/Oa+jSmJm2ST59/p4Nsv1tLaP8qWBXZuX1tBu2OUkz0u3jzZiz8g2fbFK5Ly7NLl13tb+OofjvLTD63h5jSjQG2OET7wyB48/iBP3bsx466Vh7Y38INX6tn9wDWRiXsA/7nrJN9+sZZ/umXphIifYswvry7MqL6ipW8k7dRYPBaV2pLe2JzudUWcofGIgJfC7jeRQs9vGwK80ir4xgc2oZ9CIW3Y4+ext5o50OzgxmWl3LmugvXV9ilnzLcPjHKiY2jSYxSMI50Udb6OzufEMWt9KN8vpt4I20y6aZX1ToRqzKcgGJTnVLhE4aVjndz72AEe+fBarl9awvbabh7e2ciR1gGq7Tn89fWLuGl56bR99mGPn5v/43Ukkm1fuGLSXeiJM0Pc89h+ugY9PPPZzVlpx0qGp/a1sO+0gwdvmjwdMR1IKTnY4uDV+l60QpBj1GI16lg3t2Baay2klIx4A2kJ8WQTjz/A42+18OMdDQyEPcKK8MjU+69byNqq6d1YBYKSdz/0OsMeP3/+8lUpq3AdaO7nC08cxun28cQ9G7KyMWrtH+GK7+zky++6hC9cu5BAULL9RBf3PX6Amy4t4+G/WD0hirDnZB9BKeP2TqeCxx9gd0Mv2bilCwFXLCxOOsXmdPvY2zR1keo7bQM8tKOR+69dOKkDcqRtgEf3NDPs9vP+NRVcv7SEhSW2pASQXB4/e072JfW50yWVjU42OS+MuRDiRuA/AC3wSynltyY7PpvG/FwlEJRc+Z2dkZ1oXZeT2YVmPnv1Au5YW3lWquwPNPdz5yN7KMsz87X3LOWGZRPHv247eoa//s0Rcs06HvnwWlZnMG5T5fzD6fbR2j/KXLsl6Rnc2eKNxl4+9Mu9PHDj4qRbSP2BID/ZeZKHdjRQnm/ipx9cm5W8u8LdP99Dc98Ia+YU8MbJXgZGfFQX5/Dc57bE9Sr3nOyj2GbIygbwYIsjK8p4RVZDyt/j3Q29U4oq+QJBvvjUYTZVF/HhKJnhoJQ4XF46h9y81dTPnqY+KvLNfGrzPOYUWci36FlbVZB0OuXV+h58CQa5ZIpWI9iy0D6tqaxETGbMz4kCOCGEFvgJ8C6gDdgnhHhOSnl8Zj/ZzKLVCD66sYp/f6GWBbOs/PCuldyyovystsqtrSrkiU9v4B+fPcZ9jx/gioV2/urqBaEbuGOU4x1DPH2wjTVz8nnkw2unLbSscu5iM+lZWj4zwhmbF9i5bkkJD+9oYMjtY1l5LkvLcim2GXF5Agx7fDjdfoY9fpxuP063j98daGPfaQe3rirnG7cuz3re84Prq/jCE4c40OzguiUlbFlgZ+uSWZOGh7P1vSnNNWXFmKeTHim2GadsadVrNSwry2V3Yy/HzgwhAAE4Rn2RKWoaAe+5tIz3rChDp9Wg1QiWluemVBeRb9anJbyVDMU244wY8qk4JzxzIcRG4OtSyhvC/34QQEr574leczF45hDyzo+2D7KiIm9GUwH+QJDH32rm+6/Ux0yHshi0vH9NBV97z9KzPiVIRQVCueLPP3mIY+2D+JNQYLMadXzj1mXctnp6Zi1IKekd9mK3JicocqR1YMqirmTxB4K81tATmVeeDhoNXLmwOGWnweMPsOdkX+xs8jic7nWxvbaboJRICRJJnllPWZ6Z0lwT5fmmmA3W4jJbTP1BMjT3uWjoSn988mSsrSqgIE5b4dngnPfMgQogunG1DVg/Q5/lnEKrEazK0hc9E3RaDR/fPI9bVpazv9lBaa6J2YUWCiz6c0YBSeXiZE6RhWc/uxmPP0BD1zDHO4ZwjHixmfRYTTpsRh02kw6rSYfVqMNuNU7rlCshREqSytksEtRpNRRbTRlNPbRbjWlF/4w6LZeU2Dg+RfHZXHsOn9qSXOtkaZ4pZUMOE/vNs0Wh1TBjhnwqzhVjHs8aTNjeCSHuAe4BmDMne1KdKslTZDWmPM5UReVsYNRpWV6RNy2dHdNJtqduleQZMzLmpRlsLsrzzZwZdKc8rjceBTl6lpalV0gbGiNLRhGK8RRaDaysnHnHKhHnSuC/DYjWnqwEJmjvSSl/LqVcJ6VcV1ycvQESKioqKjNFtmtg7DlGdNrUo2V6nYb5s6wZq/UtKbNFppelS45Rx4rK/IzmOmRLRhlCBYGrKuPr058rnCvGfB+wUAgxTwhhAO4Gnpvhz6SioqJy3qHRiJQ0wy1GLYtKbWxZYGeePSfj2hyLQcc8e/r9+gadhtVz8jMuMltekUd1cU7GKZVim5GVGWwszhbnRJhdSukXQnwOeIlQa9p/SSmPzfDHUlFRUTkvmWu3YLcZMOq0GHUadBoRqu73+nF5/GjDnmuuWTcthatziyxYDFq8/iDeQDAy9VAgIiJLHn8Ajz+IxxfEqNdgM+mwmfQU5RiyUtNg0mupLrYyz55Dn8uL0+0PF91JAsFQm5w3EMTnDxKQEp1Gg04r0GkEOUYdeWY9uSZ9xvMlzhbnhDEHkFJuA7bN9OdQUVFROd+xGHQTev7zLJppGb4TDyHEtKr/pYIQArvVmNVhP+ci58eWQ0VFRUVFRSUhqjFXUVFRUVE5z1GNuYqKioqKynmOasxVVFRUVFTOc1RjrqKioqKicp6jGnMVFRUVFZXzHNWYq6ioqKionOeoxlxFRUVFReU8RzXmKioqKioq5znnxDzzdBBC9ADNGbyFHejN0se5EFHXZ3LU9ZkcdX0mR12fxKhrk5gqKWXcKWPnrTHPFCHE/kRD3lXU9ZkKdX0mR12fyVHXJzHq2qSHGmZXUVFRUVE5z1GNuYqKioqKynnOxWzMfz7TH+AcR12fyVHXZ3LU9ZkcdX0So65NGly0OXMVFRUVFZULhYvZM1dRUVFRUbkgUI25ioqKiorKec4FZcyFEP8lhOgWQtREPbdKCPGWEOKwEGK/EOLy8PMfEkK8E/7zphBiZdRrbhRC1AkhGoUQfzcTv8t0kMr6RP38MiFEQAhxR9RzHxNCNIT/fOxs/g7TRaprI4S4Ovz8MSHEq1HPX/TXjhAiTwjxRyHEkfD6fCLqNRfctQMJ12elEGKPEOJoeD1yo372YPgaqRNC3BD1/MV0/cRdHyHEu4QQB8LPHxBCbI16zdrw841CiIeEEGImfp9zEinlBfMHuBJYA9REPfcycFP48c3ArvDjTUBB+PFNwN7wYy1wEqgGDMARYOlM/25ne32i1mIHsA24I/xcIdAU/rsg/Lhgpn+3s3zt5APHgTnhf89Sr52Y9fkq8O3w42KgP7weF+S1M8n67AOuCj/+JPCN8OOl4WvDCMwLXzPai/D6SbQ+q4Hy8OPlQHvUa94GNgICeEG5/tQ/8sLyzKWUrxG6ccQ8DSg74jygI3zsm1JKR/j5t4DK8OPLgUYpZZOU0gs8CbxvWj/4WSKV9QnzeeBpoDvquRuAV6SU/eH1ewW4cXo+8dkjxbX5IPB7KWVL+LXK+qjXztjztrDXZA2/zs8Feu1AwvVZBLwWfvwKcHv48fuAJ6WUHinlKaCR0LVzsV0/cddHSnlISqlcS8cAkxDCKIQoA3KllHtkyLI/Ctw6/Z/+/EA30x/gLHA/8JIQ4nuE0gqb4hzzKUK7PIAKoDXqZ23A+mn9hDNL3PURQlQAtwFbgcuijo+3PhVn56OedRJdO5cAeiHELsAG/IeU8lHUa0dZn4eB5wgZdxtwl5QyGL6mLpZrB6AGeC/wLHAnMDv8fAUhB0Iheh0upusn0fpEcztwSErpCV8/bVE/u9Cvn5S4oDzzBHwG+JKUcjbwJeBX0T8UQlxDyJg/oDwV5z0u5P69ROvzI+ABKWVg3PEX0/okWhsdsBZ4NyFv82tCiEu4uNYGEq/PDcBhoBxYBTwczodebOvzSeCzQogDhDY13vDzidZBXZ8ohBDLgG8D9ypPxXmPC3l9UuJiMOYfA34ffvxbQqEsAIQQK4BfAu+TUvaFn24jdodYSWzo+UIj0fqsA54UQpwG7gB+KoS4lYtrfRKtTRvwopTSJaXsJRQqXMnFtTaQeH0+QSgNIaWUjcApYDEX2fpIKWullNdLKdcCTxDKh0PidVDXJ4wQohL4A/BRKWX0ulVGvcUFvT6pcjEY8w7gqvDjrUADgBBiDqEb0UeklPVRx+8DFgoh5gkhDMDdhEKGFypx10dKOU9KOVdKORf4HfBXUspngJeA64UQBUKIAuD68HMXInHXhlBY8AohhE4IYSEUCj2Beu0o69MCXAsghCghlBtt4uK6dhBCzAr/rQH+AXgk/KPngLvDeeB5wEJChV0X1fWTaH2EEPnAn4AHpZRvKMdLKc8ATiHEhnA9xkcJfRdV4IKrZn8COAP4CO3iPgVsAQ4QqgzdC6wNH/tLwEEoHHgY2B/1PjcD9YR2in8/07/XTKzPuNf9D+Fq9vC/P0moaKcR+MRM/14zsTbA3xKqaK8B7levnZjvVjmhSvej4fX58IV87UyyPl8MXwv1wLcIK26Gj//78DVSR1RF9kV2/cRdH0KG3RV1bz7MWMfIuvA1dZJQbYaYqd/pXPujyrmqqKioqKic51wMYXYVFRUVFZULGtWYq6ioqKionOeoxlxFRUVFReU8RzXmKioqKioq5zmqMVdRUVFRUTnPUY25ioqKiorKeY5qzFVUVFRUVM5z/h+9S/RVfTO8ggAAAABJRU5ErkJggg==\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 (pydata)",
   "language": "python",
   "name": "pydata"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}