copula_path_forecast.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copula Model Simulation Study in Path Forecasting\n",
    "Compare computational time of the copula model to standard path forecasting (via sequential simulation) as a function of: Samples drawn from the forecast distribution (nsamps)\n",
    "Forecast horizon (k) State vector size (p)\n",
    "Test standard path forecasting both with and without the novel interpolation to speed up the variational Bayes step. Compare resulting forecast distributions to ensure that they match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import matplotlib\n",
    "from pybats.conjugates import pois_conjugate_params\n",
    "from pybats.forecast import forecast_aR\n",
    "from pybats.plot import ax_style\n",
    "matplotlib.use('Agg')\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import cProfile, pstats, io\n",
    "import time\n",
    "from pybats.dglm import pois_dglm, bern_dglm\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define timing function and simulation from Poisson DGLM\n",
    "def time_func(f, nruns, **kwargs):\n",
    "    start = time.time()\n",
    "    for i in range(nruns):\n",
    "        f(**kwargs)\n",
    "    end = time.time()    \n",
    "    return (end-start)/nruns\n",
    "\n",
    "def time_path_forecast(p, k, nsamps, repetitions = 1, copula=True, interpolate=True):\n",
    "    a0 = np.ones(p)\n",
    "    R0 = np.diag(0.1 * np.ones(p))\n",
    "    \n",
    "    pmod = pois_dglm(a0 = a0, R0 = R0,\n",
    "                 ntrend = 1,\n",
    "                 nregn = p - 1,\n",
    "                 deltrend = .99,\n",
    "                 delregn = .99,\n",
    "                 interpolate=interpolate)    \n",
    "    if p > 1:\n",
    "        X = np.random.normal(0, 0.1, [k, p-1])\n",
    "    else:\n",
    "        X = None\n",
    "                    \n",
    "    time = time_func(pmod.forecast_path, nruns=repetitions, k = k,\n",
    "                     X = X, nsamps = nsamps, copula=copula)    \n",
    "    return time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k array: [5, 10, 20]\n",
      "p array: [5, 10, 20]\n",
      "nsamps array: [   1   10   50  100  500 1000]\n"
     ]
    }
   ],
   "source": [
    "#Define empty dataset to hold results\n",
    "kps = [5, 10, 20]\n",
    "nsamps = np.array([1, 10, 50, 100, 500, 1000])\n",
    "print('k array: ' + str(kps))\n",
    "print('p array: ' + str(kps))\n",
    "print('nsamps array: ' + str(nsamps))\n",
    "\n",
    "rows = []\n",
    "i = 0\n",
    "for k in kps:\n",
    "    p = k\n",
    "    for n in nsamps:\n",
    "        rows.append(pd.DataFrame({\n",
    "            'k':k,\n",
    "            'p':p,\n",
    "            'nsamps':n,\n",
    "            'standard_time':0.0,\n",
    "            'standard_time_interp':0.0,\n",
    "            'copula_time_interp':0.0}, index=[i]))\n",
    "        i += 1            \n",
    "res = pd.concat(rows, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n"
     ]
    }
   ],
   "source": [
    "# Run Simulation\n",
    "i = 0\n",
    "for k in kps:\n",
    "    p = k\n",
    "    for n in nsamps:\n",
    "        if 0 < n < 10:\n",
    "            repetitions = 32\n",
    "        if 10 <= n < 50:\n",
    "            repetitions = 32\n",
    "        if 50 <= n < 100:\n",
    "            repetitions = 16\n",
    "        if 100 <= n < 1000:\n",
    "            repetitions = 10\n",
    "        if 1000 <= n < 10000:\n",
    "            repetitions = 5\n",
    "        res.loc[i,'standard_time'] = time_path_forecast(p, k, n, repetitions=repetitions, copula=False, interpolate=False)\n",
    "        res.loc[i,'standard_time_interp'] = time_path_forecast(p, k, n, repetitions=repetitions, copula=False, interpolate=True)\n",
    "        res.loc[i,'copula_time_interp'] = time_path_forecast(p, k, n, repetitions=repetitions, copula=True, interpolate=True)       \n",
    "        i += 1\n",
    "        print(i)        \n",
    "res['speedup_interp'] = res['standard_time_interp'] / res['copula_time_interp']\n",
    "res['speedup_total'] = res['standard_time'] / res['copula_time_interp']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>k</th>\n",
       "      <th>p</th>\n",
       "      <th>nsamps</th>\n",
       "      <th>standard_time</th>\n",
       "      <th>standard_time_interp</th>\n",
       "      <th>copula_time_interp</th>\n",
       "      <th>speedup_interp</th>\n",
       "      <th>speedup_total</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>1000</td>\n",
       "      <td>3.61</td>\n",
       "      <td>0.74</td>\n",
       "      <td>0.01</td>\n",
       "      <td>89.65</td>\n",
       "      <td>435.29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>7.66</td>\n",
       "      <td>1.63</td>\n",
       "      <td>0.02</td>\n",
       "      <td>99.79</td>\n",
       "      <td>470.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>20</td>\n",
       "      <td>20</td>\n",
       "      <td>1000</td>\n",
       "      <td>15.71</td>\n",
       "      <td>3.30</td>\n",
       "      <td>0.07</td>\n",
       "      <td>45.68</td>\n",
       "      <td>217.45</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     k   p  nsamps  standard_time  standard_time_interp  copula_time_interp  \\\n",
       "5    5   5    1000           3.61                  0.74                0.01   \n",
       "11  10  10    1000           7.66                  1.63                0.02   \n",
       "17  20  20    1000          15.71                  3.30                0.07   \n",
       "\n",
       "    speedup_interp  speedup_total  \n",
       "5            89.65         435.29  \n",
       "11           99.79         470.11  \n",
       "17           45.68         217.45  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Save / load simulation results\n",
    "pd.to_pickle(res, 'sim_study_data')\n",
    "res = pd.read_pickle('sim_study_data')\n",
    "res[res.nsamps == 1000].round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f93b2891f70>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAEGCAYAAACJsIcWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABjOklEQVR4nO3dd3gUVRfA4d9JJ6QQeg2ht9BDFwSpAoqAdJGmKKBiwd7bJyj2CgKCdAQEBKQqTem99xZqAiQB0pP7/TGTGCCETbKb3Wzu+zx5slN25mTh7J1y51xRSqFpmqZpmnNxsXcAmqZpmqZZn27gNU3TNM0J6QZe0zRN05yQbuA1TdM0zQnpBl7TNE3TnJCbvQPIjsKFC6ugoCB7h6FpDm/79u3hSqki9o4jIzqfNe3eMpPLubqBDwoKYtu2bfYOQ9McnoictncM96LzWdPuLTO5rC/Ra5qmaZoT0g28pmmapjkh3cBrmqZpmhPK1ffg05OQkEBoaCixsbH2DsVpeXl5Ubp0adzd3e0diubkdD7bls5l5+Z0DXxoaCi+vr4EBQUhIvYOx+kopbhy5QqhoaGUK1fO3uFoTk7ns+3oXHZ+TneJPjY2lkKFCukvAxsREQoVKqTPqPIYEZkkIpdFZF+aeQVFZKWIHDV/B5jzRUS+EZFjIrJHROpldb86n21H57Lzc7oGHtBfBjamP988aTLQ4bZ5rwGrlVKVgNXmNMCDQCXzZyjwY3Z2rP+/2Y7+bJ2bUzbwmpYnKAXXL+XQrtQ64Opts7sAU8zXU4BH0sz/VRk2AQVEpESOBKppudX1i0ZOW5Fu4K3s1KlTBAcH2zuMexo4cCDlypWjTp061KlTh127dtk7JC0zkhJh4Qj4uRVE397u5phiSqkL5uuLQDHzdSngbJr1Qs15dxCRoSKyTUS2hYWF2S7SLNL5rOUIpWDOANj4vVU363Sd7DTLffbZZzz66KP2DkPLin1zYdd0aPk65AuwdzQopZSIZPr0Qyk1HhgPEBISYt3TlzxG53MuJgLdJ0ByolU3q8/gbejEiRPUrVuXrVu3Zms7LVu2ZOTIkdSpU4fg4GC2bNlipQi1XKtWLxjwB7R8zfhysI9LKZfezd+XzfnngDJp1ittzsvVdD5rVrd7Nsx/CpKToUAZKGjdpxmc+gz+/T/2c+B8lFW3Wb2kH+8+VOOe6x0+fJjevXszefJkateufcuy69ev07x583TfN2PGDKpXr37H/OjoaHbt2sW6desYPHgw+/btu2X54cOH6dWrV7rbXLNmDQUKFLhj/ptvvskHH3xA69atGT16NJ6envf8uzQ7un4Rfn8KOn0BhSpAuRb2jmgRMAAYbf5emGb+MyIyC2gERKa5lJ9lOp8NOp+dxL/fwYo3Iag5JMaAR36r78KpG3h7CQsLo0uXLsyfPz/d5Pb19c30PbI+ffoA0KJFC6KiooiIiLglyatUqZKpbX7yyScUL16c+Ph4hg4dypgxY3jnnXcyFZOWg64ch6ld4WY4RJ0zGvgcJCIzgZZAYREJBd7FaNjniMgQ4DTQ01x9KdAROAZEA4NyNFgr0/msWZVSsOo9+OcrqPYwdPsZ3L1ssiunbuAtOTK3BX9/fwIDA9mwYUO6XwhZOeK//XGW26cze8RfooTRqdnT05NBgwYxduzYu/49mp1d2A3TukNyEgz8A0rVz/EQlFJ97rKodTrrKmCEtWPQ+WzQ+ZzLLXsdNv8I9QdBp8/BxdVmu7pnAy/G/7x+QHml1AciEggUV0rpG0d34eHhwe+//0779u3x8fGhb9++tyzPyhH/7NmzadWqFRs2bMDf3x9/f/9blmf2iP/ChQuUKFECpRQLFizIFT2F86Tzu2ByZ8hXAB6bD0UqZ2tzOp8zT+ezZlXVHjI6xt7/is37z1hyBv8DkAw8AHwAXAfmAQ1sGFeulz9/fhYvXkzbtm3x8fHh4Ycfztb2vLy8qFu3LgkJCUyaNCnb8fXr14+wsDCUUtSpU4effvop29vUbKBwZajeBVq9Af7pPmmWWTqfs0Dns5YtMRFwbBXUfBSCmhk/OcCSBr6RUqqeiOwEUEpdExEPG8eVawUFBaV2mClQoEC2e9ymeOyxx/jqq6+ssi2Av/76y2rb0mxg/+9QoTV4+cEjVn02VudzJuh81rLt+kWY2g2uHIUyDaFAYI7t2pLH5BJExBVQACJSBOMMQNM0a1MK1o2F3wbCv9/YYg86nzUtp1w5DhPbwrVT0Hd2jjbuYNkZ/DfA70BREfkYeBR4y6ZRabdYs2aNvUPQckJyMix/w+iAU7Mn3P+qLfai89nOdD7nEed3GZ1jUTBwMZTK8phLWXbPBl4pNV1EtmP0lhXgEaXUQZtHpml5SVICLBgOe+dAo2HQ/n/gYv06VDqfNS2HnN8J7t7Q/3coXNEuIVj6mNwlYL25fj4RqaeU2mG7sDQtj7kZDqf/gdbvwH0v2rp3rc5nTbOVm1cgfyEIGQQ1e4Cnj91CseQxuQ+BgcBxzPt25u8HbBeWpuURsZHg4Qt+JWD4RvDyv/d70ohPTMbdVSwe9lPns6bZ0NaJsPJdGLQUStSya+MOlp3B9wQqKKXibR2MpuUpkedgWjeo2Abaf5zpxv1k+E2GT9/B403K0qehxZ13dD5rmrUpBWvHwJpPoHIHKGSfS/K3s+Qm3z6ggI3jcBq5ZXjJ7777jooVKyIihIeHp85XSvHcc89RsWJFatWqxY4d+sqtTYQdgYntjEa+codMv33p3gs89O0GLkTGUNwvU2UudT5ngs5n7Z6Sk2DpKKNxr90Xek0DD297RwVYdgb/CbBTRPYBcSkzlVLZq/Sg2VWzZs3o3LkzLVu2vGX+n3/+ydGjRzl69CibN29m2LBhbN682T5BOqvQ7TD9UaNE5aAlUKL2vd9jik9M5n9LDzL531PUKVOA7/vVo1SBfJnZu85nJ6Tz2Y52ToWtE6Dpc9D2A3uO7ngHSxr4KcAYYC/6edlMOXHiBN27d2f8+PE0aJD1QmEtW7akdu3arF27lsTERCZNmkTDhg2zFVvdunXTnb9w4UIef/xxRITGjRsTERGRWgZTs4K4GzCjB3j6Gr1rMzFoTOi1aEbM2MnusxEMblaO1x6siodbpnva63zOIp3PWrrqPAbehYwStA7GkgY+Willk4obOeKXTnfOq/EINHwS4qNheo87l9fpC3X7Gb0h5zx+67JBSyzabW4YXjI9586do0yZ/4byLl26NOfOndNfCNbi6QPdJ0DR6uBb3OK3/XXoEi/M3k1ysuLHfvV4sGaW/z10Pqel81nLihthsOQFePAzo4OsAzbuYFkDv15EPsEY4zntJT19M+cucsPwkloO2/Kz8Uxs3X5QwfIO64lJyXy+8gg/rjlO9RJ+/NCvHkGFszVutM7nTNL5rN3i2mlj6Oao89DgiNHAW8H5iBhWH7pM/8ZlrbI9sKyBT7n20zjNvNzzWE1GR+ge3hkvz1/I4iP8tHLD8JJ3U6pUKc6ePZs6HRoaSqlSVhnkJG9SCtaMhrWjoWpn42zSwnt0l6JieXbmTracvEqfhoG8+1B1vNyzPbSkzudM0vmspbq036grnxgLjy+EwEZW2eyK/Rd5ee4ekpIV7asXo2jmOs7elSWV7FpZZU95SG4YXvJuHn74Yb777jt69+7N5s2b8ff315fzsio5Cf58xeiAU6cfPPSNxY37P8fCGTlrJzfjkviyV2261i1tlZB0PmeezmcNMCrT/doF3PPD4GVQtFq2NxmbkMT/lh7k142nCS7lx7d96lmtcYcMGngReUwpNU1EXkxvuVLqC6tF4YQcfXjJb775hk8//ZSLFy9Sq1YtOnbsyIQJE+jYsSNLly6lYsWKeHt788svv2R7X3lScjLMG2KMCpeJ3rXJyYrv/j7Gl6uOUKGIDzOfrEelYr7WiKgggM7nrNH5rFGgLJRrYZSRtsKgMccuX+eZGTs5dPE6Q+4rxysdquDplu0rdLcQpVT6C0SGKqXGi8i76S1XSr1v1UiyICQkRG3btu2WeQcPHqRatewfWTmSli1bMnbsWEJCQuwdSipn/Jytbs0YcPeCZiMtWv3KjTien72L9UfDeaROST7uWpP8npZWk86YiJxWSgXpfLY/R8tnZ/yMrerIcih3v5HLVqCUYs62s7y36AD5PFz5vEdtWlUtavH7RWS7Usqi/zwZfXs0BcY7QuJrWq5xMxwiQ6FkHWhp+Whw205d5ZkZO7kaHc//utakT8MyFpeftZAPOEZDrmm5xj9fw8p34IG3ocWobG8uKjaBN+bvZfGeCzStUIgve9WhmBUvyd8uowa+ls32qmWKHl4yl4g4Y/SuTYiB53aCm+c936KUYsL6k4xedojSAfmYP6wpwaUyV7LWQpmqhqPZjs7nXCA5GVa+DRu/g+Duxm22bNp1NoJnZ+7gfEQsL7evwtP3V8DVxbZFcTJq4L1FpC7GkJJ30I/VaFoalw8ajXt8NPSdbVHjHhmdwKi5u1l54BIdahTn0x618PNyt1WELjqfNc0CSQmw8BnYMwsaDoUOY7I1dHNysmL8+hOMXX6YYn5ezHmqMfXLFrRiwHeXUQNfCvic9L8QLH6sRkRcgW3AOaVUZxEpB8wCCgHbgf5KqXgR8QR+BeoDV4BeSqlTlv4hmmY3Z7cYBVbcvIxRpIrfu3b53tBIhs/YzoWIWN7pXJ1BzYKsfUn+dh5YIZ81zelFhsLR5dDqLeOyfDby8vL1WF6as5v1R8PpWLM4n3SrhX8+mx3E3yGjBv6YUsoaST8SOAj4mdNjgC+VUrNE5CdgCPCj+fuaUqqiiPQ210v/QVBNcyRbxoN3QaP0bEBQhqsqpZi26TQfLj5IYR8P5jzdhHqBATkRZZyV8lnTnFP8TaMYVcFy8Mw2yF84W5tbdySMF+fs4npsoq361dxT1q87WEBESgOdgAnmtGCcKcw1V5kCPGK+7mJOYy5vLTn9aWhaZiSaI64+/B0MWXnPxv1GXCLPzdrF2wv307RiIZY81zynGndN0zISeQ5+fgDWf25MZ6Nxj09M5pOlB3l80hYK5vdg0TP30bdRYI437pBxA295F+C7+wp4hf8GtSgERCilEs3pUIxbAZi/zwKYyyPN9W8hIkNFZJuIbAsLC7NCiNaVW4aX7NevH1WqVCE4OJjBgweTkJAA6OElLbbxB+MLISbCeHzmHl8Ihy5G8fB3G1iy5zwvt6/CpAENCMjvkTOxGkJzcmfOQudzHpB26OYy2Rv058yVaHr89C/j1p2gb6NAFo64jyrFrVLHIkvu2sArpVZkZ8Mi0hm4rJTanp3t3E4pNV4pFaKUCilSpIg1N52n9OvXj0OHDrF3715iYmKYMGECcOvwkuPHj2fYsGF2jtTBKAWr3oflrxuX8tzv3Tn9t21neeT7f7gem8j0JxozolVFXGzcezYdUTm9Qy3n6HzOotDtMKk9JMUbZYzLtcjyphbuOkfHb9ZzMvwmP/arx/+61iSfh3UL12SWLS/RNwMeFpFTGJ3qHgC+BgqISMq9/9LAOfP1OaAMgLncH6OzXa514sQJ6taty9atW7O1nZYtWzJy5Ejq1KlDcHAwW7ZsyXZsHTt2REQQERo2bEhoqHGCd7fhJTUgKRH+eA42fAH1B0GPyRn2lo+JT+KVubt5ee4e6pYJYMlz99Gkwh0XpbRcQuezk4m+ClMfAS8/GLIcStS+51vS3Ux8Iq/M3c3IWbuoUtyXpSObZ2e0R6uyWQOvlHpdKVVaKRUE9Ab+Ukr1A/4GHjVXGwAsNF8vMqcxl/+l7lZmLxMGLRvEgmMLAEhITmDQskH8cfwPAGISYxi0bBDLTi4D4Hr8dQYtG8Sq06sAuBZ7jUHLBrHm7BoAwmPCLd7v4cOH6d69O5MnT75j7Ojr169Tp06ddH8OHDiQ7vZShpf84YcfGDx4cLr7u9s2IyIi7hpnQkICU6dOpUOHDsDdh5fUgNXvwY5focUr0PlLcLn70fmJsBt0/eEf5mwL5dkHKjLtiUYU9bVdQYu8Quezzmer8S4IXb6HwSugYPksbWL/+Ug6f7uB37aH8kyriswe2pjSAd5WDjTr7lkHU0QqAy8DZdOun40eua8Cs0TkI2AnMNGcPxGYKiLHgKsYBwW5Um4aXnL48OG0aNHirqNhaWk0Hg6FK0O9xzNcbfGe87w6dw8ebi5MHtSAllUsL0NpazbIZ6en89nJbB5n1JKv8iBUz9qYAkopft14mo+XHKSAtzvThzSiacXs9bq3BUsKXf8G/AT8DCRlZSdKqTXAGvP1CeCOngxKqVigR1a2n5FfOvw3uIK7i/st0/nc8t0y7evhe8t0gFfALdOF81n2D5hbhpd8//33CQsLY9y4canz9PCSt7lxGTb9CA+8BX4lM2zc4xKT+N+Sg0zZeJp6gQX4rm89ShZwuAJy2c5ne9L5rPM5y5SCvz6C9WMh+FGjgc+CazfjeWXeHlYeuESrKkUY26M2hXzuXdjKHixp4BOVUj/aPBInkhuGl5wwYQLLly9n9erVuKSp0qSHl0zj6kmjOt2NS0a5ygwK2Jy9Gs0zM3awOzSSJ+4rx6sPVsXd1aZPoWaV1fNZRF4AnsAomLMXGASUIJ2CVtbcb07R+ewEkhJhyYuwY4pxkN7pyyxtZvOJKzw/exfhN+J4q1M1htxXzi6Pv1nKkgb+DxEZDvwOxKXMVEpdtVlUTsDRh5d8+umnKVu2LE2aNAGgW7duvPPOO3p4yRQX98G0bkbv2gF/ZNi4rzpwiRfn7EIBPz1Wnw7BxXMuzsyzaj6LSCngOaC6UipGROZg3F7rSPoFrXIlnc+5WFIC/DYQDi2G5qOMq3GZbJQTk5L59q9jfPvXUQILejN/WDNqlrbJmBFWZUkDn9Lx7eU08xSQtV4JTi4oKIh9+/YBUKBAgWz3uE3x2GOP8dVXX1llWwCJiYnpzhcRvv/+e6vtJ1c6/S/M6A2ePkbjXqRKuqslJCUzdsVhxq09QY2SfvzQrx5lC+XP4WAzzRb57AbkE5EEwBu4gPHUTMqp7hTgPXJhA6/z2Qm4uBm31zqMgcZPZ/rt5yNieH7WLracukq3uqX44JFgfKw0jLOt3TNKpVS5nAhE0xyGuELBIOg1HQqUSXeVi5GxPDtzB1tPXaNfo0De7lwdL3f7PvNqCWvns1LqnIiMBc4AMcAKjEvydytodQsRGQoMBQgMDLRmaFped/0SxEVB4UrQ8bMsbWLF/ou8Mm8P8YnJfNGzNt3qlbZykLZlSS/6dHsVKaV+tX44Wnr08JI55OI+41J8YCN4cs1dR5BafzSM52ftIiYhia9716FLndzTacna+SwiARhlpssBERid+DpY+n6l1HhgPEBISEi2H4vNDXQ+54CrJ4z+M64eMHxTho+0pic2IYlPlhodZoNL+fFtn3qUK+zwV+fuYMl1hrQPfXoBrYEdGCO/OSSllEN3fMjtrFCewPFs+ApWvQt9Zhm9a9Np3JOSFd+sPso3fx2lUlEffuhXj4pF7VeGMousnc9tgJNKqTAAEZmPUeSqgIi4mWfxaQtaZZrOZ9txyly+sBumPQrJCdBvbqYb92OXr/PMjJ0cunidIfeV45UOVfB0c/yrc+mx5BL9s2mnRaQARu9Yh+Tl5cWVK1coVKiQ/lKwAaUUV65cwcvLSYq2JCfDqnfg32+hRjeo0Drd1cJvxPH8rF1sOBZOt3ql+OiRYLw9csd9uLRskM9ngMYi4o1xib41xvDQKQWtZnFrQatM0flsO06XywAn18OsvuDpBwMX37X/THqUUvy2LZR3F+0nn4crvwxsQKuqjlPDIiuy8g11E+NynEMqXbo0oaGhOOJANM7Cy8uL0qVz172odCUlwKLnYPcMaPAkPPhpumfuW05e5dmZO4iITmBM95r0DMn5YR9tKFv5rJTaLCJzMa4CJGIUrxoPLCH9glaZovPZtpwml8F4zn3Dl+BbAvrPB3/L/66o2ATe/H0ff+w+T5Pyhfiqdx2K+eX+Ax9L7sH/gdHLFozSttUx7rM5JHd3d8qVc9jjD82RnFhrNO4t34D7X7nj0ZnkZMX49Sf4bPlhygTk45fhDale0s9OwVqHLfJZKfUu8O5ts9MtaJVZOp81iyQlgKs79PgFkpOMMrQW2nU2gmdn7uB8RCwvt6/C0/dXwDXnB4OyCUvO4MemeZ0InFZK6aEntdxLKaMxr9QGnlqX7iATEdHxjPptN6sOXqZjzeKM6V4LXy93OwRrdTqfNeehFKwbC8dXQ//fwcvyZ9NTDuDHLj9MMT8v5jzVmPplLT8wyA0ybOBFxBXYr5QKN6c9gIEi8oJSqlpOBKhpVhV1AWb3g3YfQ9km6Tbuu89GMHz6Di5fj+W9h6ozoGmQU1yS1/msOZXkZFj2KmwZD7V6Gc+7WyjsehwvztnF+qPhPBhcnNHdauHv7RQH8Le46yciIr2BccBNETkKfAxMArYC/XImPE2zoivH4ddHIOaqUaHuNikDSHy05ABFfb2Y81QT6gYG5HycthGAMYiTzmct90uMhwVPw7550OQZaPvhXR9rvd26I2G8OGcX12MT+bhrMH0bBjrFAXx6MjrkeQuor5Q6JiL1gI3Ao0qpP3ImNE2zovO7YFp34/XAxVCy7i2Lr8cm8Nr8vSzZc4EHqhbli561KeDtkfNx2k5JoKbOZ80pLHnBaNzbvA/NRlpUejY+MZnPVxxm3LoTVC7mw/QnGlOleK57zDVTMmrg45VSxwCUUjtE5Kj+MtBypbDDMLkz5Asw7tMVrnjL4oMXohg+fQdnrkbzaoeqPNWiPC5O0skmjWSdz5rTaPY8BDWH2paNKn7mSjTPztrJ7rMR9G0UyNudqpPPI3c+254ZGTXwRUXkxTTTBdJOK6W+sF1YmmZFhSpCwyeg4VCjJnUac7ae5e2F+/DP586MJxrRqHwhOwVpc+46n7VcLeIs7J4JLV42ys8WrmTR2xbtPs+b8/eCwA/96tGxZt4ZTS+jBv5nwDeDaU1zbLtmQrnmxvOwbd67ZVF0fCJvL9jPvB2hNKtYiK961aWIr2OO6WwlYeh81nKrywdhajeIv2l0qAsoe8+3RMcn8t6i/czZFkr9sgF83bsOpQO8cyBYx3HXBl4p9X5OBqJpVpPy6MzfH0HDp6Djp7csPnb5BsOnb+fo5Rs817oSI1tXcprnXjNwQee0liud3QLTe4CbFwxaalHjfuB8FM/O3MGJ8Js806oiz7ephJurZZ3wnEnuq7WpaRlJToZlr8GWcVCrN7T/+JbFC3ed4/X5e/Fyd2XKoIa0qFzEToFqmnZPR1fC7P7gV8LoPxMQlOHqKU/CfLz0IAXyuTN9SCOaViycM7E6IN3Aa84jMR4WDIN9c+94dCY2IYmPlhxg2qYzhJQN4Nu+dSnhn8/OAWualiGloFgNYxAon4wPxiOi43l57h5WHrhEqypFGNujNoV8nPq22z3pBl5zHomxEH74jkdnzlyJZviM7ew7F8XQFuV5uX0V3PPg5TpNyzUu7Tca9srtoGKbez7jvvnEFZ6fvYvwG3G81akag5uVc8YnYTLNogZeRDoBNTCGlwRAKfWBrYLStEyJvgru+cDLD55YDW7/HbUv33+RUb/tRoDx/evTrkZx+8XpIHQ+aw5LKWPY5n+/hUHLILBRho17UrLi27+O8s3qowQW9Gb+sGbULG15uVpnZ8lgMz8B3kArYALGEJBbbByXplkmMtToXVu0GvScktq4JyQlM+bPQ0zYcJJapf35vm89yhTMWz1o06PzWXNYSYnwx0jYNQ1ChkDpkAxXvxAZw8hZu9hy8ird6pbig0eC8fHUF6XTsuTTaKqUqiUie5RS74vI58Cftg5M0+4p7DBM7Qpx16Hzl6mzz0fE8MyMHew4E8HjTcryZqdqeLo5f1ELC+l81hxPfDTMHQxH/oSWr8P9r2ZYnW7F/ou8Mm8P8YnJfNGzNt3qOcmQt1ZmSQMfY/6OFpGSwBUg71QK0BxT6DaY/ii4uMPAJVCiFgBrj4Tx/KydxCcm822fujxUu+Q9NpTn6HzWHM+BhXBkGXT6HBo8cdfVYhOS+GTpQaZsPE1wKT++7VOPcoXz52CguYslDfxiESkAfAbswBhLeoItg9K0DCXGw9xBxtCQ/X+HguVJSlZ8teoI3/19jCrFfPm+Xz0qFPGxd6SOSOez5jhShm6u3dvoVGceqKfn2OUbPDtzJwcvRDG4WTlefbCKvjJ3D/ds4JVSH5ov54nIYsBLKRVp27A0LQNuHtBrOvgUA99iXL4ey8iZu9h44go96pfmgy7BeaLOdFbofNYcRvgxmDsQuo7LsHFXSvHb9lDeXbiffB6uTBoYwgNVi+VsrLlURsPFdstgGUqp+bYJSdPuYvN4iA6HVm+kfhlsOnGFZ2fu5HpsAp8+WoueIWXsHKTDKnC3nNb5rOW4czuMW2xIukM3p4iKTeCt3/exaPd5mpQvxFe961DMz+uu62u3yugM/qEMlilAfyFoOUMpWPMJrB0DVTtDchLJuPDj2uN8vuIwQYXyM3VIQ6oW97N3pI6sAHfPaZ3PWs45/jfMfgy8C0L/BVCoQrqr7TobwbMzd3A+IpZR7SozrGXFvFBS2qoyqkU/KCcD0bR0JSfB0lGwbRLU7Q+dv+JaTBIvztnO34fD6FyrBKO719KPx9zbKZ3Tmt2d2WTUlS9cGR6bZ5SgvU1ysuLn9Sf4bPlhivl5MeepxtQvW9AOweZ+ljwH/05683VhDM3mlIL5Q43Ss/e9AK3fZefZCJ6ZsZOw63F82KUGjzUui2TwOI12K53Pml2VrAuNh0HzlyBfgTsWh12P46XfdrPuSBgPBhdndLda+Hu753ycTsKS056baV57AZ2Bg7YJR9PSEIFK7aBkXVSTEfzyzyk++fMgxfy8mDusCbVKF7B3hLmRzmctZykFW8ZDzR7GZfl2H6a72rojYbw4ZzfXYxP4uGswfRsG6oP3bLKkF/3naadFZCyw3GYRadqNMLi4Byq2htq9iIpN4NXpO/hz30XaVCvG5z1q66P6LNL5rOWo5CRY+jJsm2iM5d78xTtWSUhKZuyKw4xbe4JKRX2Y/kQjqhT3tUOwzicrNy69AV02SLONa6dh6iMQfQWe38v+qzBi+g7OXovhjY5VebJ5eX1Ub106nzXbSIyD+U8aRWyaPW/cZrvNmSvRPDtrJ7vPRtCnYSDvdK6uH3G1Ikvuwe/F6GUL4AoUAe55v05EvIB1gKe5n7lKqXdFpBwwCygEbAf6K6XiRcQT+BWoj1Fdq5dS6lSm/yLN7mITY9lxeQeNSzTGRTIxatul/UZd+cQYVJ/ZzNoTybuL9lPQ24PZQxsTEqQ72mRXVvNZ0zIlNgpm94OT66Ddx9D0mTtWWbT7PG/O3wsC3/etR6dauqCitVlyBt85zetE4JJSKtGC98UBDyilboiIO7BBRP4EXgS+VErNMge+GAL8aP6+ppSqKCK9gTFAr8z8MZr9RMRG8OPuHxleZzjrQtfxxoY3mNxhMvWL1bdsA2c2wYye4O5NzGNLePOfJObv3EvzSoX5qledPD+usxVlNZ81zXIJMRB1wShiU7v3LYui4xN5b9F+5mwLpV5gAb7uXVcPBGUjltyDPy0i9YD7MI78NwA7LXifAm6Yk+7mjwIeAPqa86cA72E08F3M1wBzge9ERMztaA5u66WtzD86n26VutEuqB353PJRr2g9ACbuncj+K/sZ02IM7i53uXd+eCl4F+bkg9MYOieMY2E3eKFNZZ55QD/7ak1ZzWdNs0jUechfFHyLwbB/jaqTaRy8EMUzM3ZwIvwmI1pV4Pk2lXF3zcRVPi1T7vnJmo/VTMG4pF4YmCwib1mycRFxFZFdwGVgJXAciEhzxhAKlDJflwLOApjLI8193r7NoSKyTUS2hYWFWRKGlgPalm3L8keXU6VgFTxdPWlTtk3qvXIXccFVXFMb9/Wh6wmPCTfeGBtl/G79HksaTaPjr2e4Fh3PtCGNGNmmkm7crSw7+axpGbq4D8a3hBVvGtNpGnelFL9uPEWX7//hemwi04Y04uX2VXXjbmOWXKLvB9RWSsUCiMhoYBfw0b3eqJRKAuqYg1v8DlTNcqT/bXM8MB4gJCREn93bWWRcJOdunKN6oeoU9Er/Hvmg4P/qq8QmxjJq7SjaB7XnAykKG78jZsAyPlh3nZlbztAwqCDf9q2ry1HaTpbzWdPu6tQ/MLMPeOSH+gNvWRQRHc8rc/ew4sAlWlUpwtgetfUttxxiyeHTeYznZVN4AucysxOlVATwN9AEoyZ2yoFF6TTbOgeUATCX+2N0ttMc2E+7f6L/0v7/nZHfg5ebF7M6zeTJqGhY8SZHi9ak6aJBzNm7jqfvr8CMJxvpxt22sp3PmnaLQ0tgalfjsvyQFVC0WuqiLSev8uDX6/n78GXe6lSNiQMa6MY9B1lyBh8J7BeRlRj37NoCW0TkGwCl1HPpvUlEigAJSqkIEclnvm8MRkP/KEZP+gHAQvMti8zpjebyv/T9d8f3dO2nqV+sPoXzFbbsDUmJlFv7Oeycxpnyvel+piGq0GLGdGnBo3WqcuTaEcKjw2lcMpM98DVLZSmfNS1dMdfg92FQPBj6/gb5jbuqScmKb/86yjerjxJY0Jv5w5pRs7S/nYPNeyxp4H83f1KssXDbJYApIuKKcaVgjlJqsYgcAGaJyEcYnXsmmutPBKaKyDHgKtA7vY1qjiEuKQ4PFw/8Pf1pU7aN5W/c8AXsnMbaEkMYcOABapcuwHc9+6T2op1xcAYrTq1gdc/V5HPLR1JyEq4u+rlYK8pqPt+VeQtuAhCMcdAwGDgMzAaCgFNAT6XUtezuS3Mw+QKg/+9QpAp4+gBwITKG52ftYvPJq3StW4oPHwnWY0XYiVhykiwiHkBlc/KwUirBplFZKCQkRG3bts3eYeRJ7298n/M3zvND6x8y1QCfvxzGzGk/8+3l2gxsGsQbHavh4fbfmXp8UjxHI45So1ANAJ5Y/gRVClbh5QYvW/1vyEtEZLtSKsR8bdV8FpEpwHql1ARz297AG8BVpdRoEXkNCFBKvZrRdnQ+5xLJyUZHuoAgaPTULYtWHrjEy3N3E5+YzIddguleX9dQsra0uXwvlhS6aYnR6/YUIEAZERmglFqXjRi1XK56oeoUzVfUssb9+iX460PWVhzFyHlHSEyqz/d9a6Vb2MLD1SO1cU9KTqJaoWqU8TXGeE9Wycw4OIMO5TpYfktAu4W181lE/IEWwEAApVQ8EC8iXYCW5mpTMK4UZNjAa7lAUgIsGA5750Dj4amzYxOSGP3nISb/e4oaJf34tk9dyhfxsWOgGlh2if5zoJ1S6jCAiFQGZmJUnNPyqB6Ve1i24tUTqKldSYi8xJebKlG8WF1+6FfPouR3dXHlpZCXUqf3hu9lzNYxFM5XmA7lOpCQnICruOp79Zlj7XwuB4QBv4hIbYzqlCOBYkqpC+Y6F4Fi6b1ZRIYCQwECAwOzGIKWI+JvwpwBcGwlPPC2MSIccDzsBs/M2MnBC1EMblaOVx+sgqebvq3mCCz5ZnRP+TIAUEodwShao+VBP+7+kbVn11q28oU9JE1ox/XIq/SMeZ2qIa1YMKJZlo/saxepzR+P/MEDgQ8AsPDYQjrO72hxD34NsH4+uwH1gB+VUnUxRqt7Le0KZmfZdO8FKqXGK6VClFIhRYoUyUYYmi38uv9XVp1eBUmJRk/546vhoa+hxSgUMGfbWTp/s4FLUbFMGhjCOw9V1427A7HkDH6biEwAppnTjwH6RlkeFJcUx6rTq4iKi+L+MvdnvPKZTSRO7U54gheDk95l8KMdeNQK9+OC/INSX5f0KUmjEo0o5GX03P3z5J/4e/rTtGTTbO/HiVk7n0OBUKXUZnN6LkYDf0lESiilLohICYxiV1ouopRi3tF51Clax+hIW6Mb08vXo2H5ZhSPTeDN3/exaPd5mpQvxFe96+jHWx2QJQ38MGAEkPL4zDqM0rJaHuPp6snMTjNJVskZrpecrJi2N4aKcYF84/siX/bvYJPhH5uWbJramCulmLB3AiXyl0idF50Qjbe7rnF9G6vms1LqooicFZEq5pWB1sAB82cAMJpbH4fVHNz+8P2cjDpJ5/KdWdB0DNFRRpmEq3V68tmcB+jpUZZl/1TlXMR1ujS7zHtt7iMgn27cHdFdG3jzOfYiSqkDwBfmDyJSA/DDuO+m5RHrQ9fTqEQjPFw9Mlwv6tDfPPuPN2uPRvNw7R+Y2K0m+XPgERkRYWanmUTERQBGhb3289rzUshLlvcXcG5uIlLdRvn8LDDd7EF/AhiE+WisiAwBTgM9sxO8lnOmHZzG9kvbecA1AO9Z/cifryCM2EIBjwAGl5nMN38eo5iP4r0e3ny652seCatAq8BWXIu9xvkb56lWqJruF+MgMvpX+BajVvXtCgJf2yYczRGdiTrDiNUjmLhvYobrXVz3C36zHqHkyXl89EgwX/eukyONewoPVw+KehcFIEkl0aNyD2oXqQ3AqchT/LznZyLjInMsHgcTiI3yWSm1y7yPXksp9YhS6ppS6opSqrVSqpJSqo1S6mp29qHlnA+bfcikhu/iPbMPePnDY3MJi05i4OStfLH8HG2rlmPpc83pV+cBZnaaSeOSjQFYcWoFvZf05nTUaQDCY8Lzcr45hIy+fSum9+iMUmq9iOhL9HlIoF8gP7b5kTpF69x1nfCjWwn46xW2SQ36PfkywWXt22GqoFfBW3rgb7ywkR93/0jXSl0B4wzf18M3L51peOp81jKy4dwG6hWth3dSAqUXjgQXN3h8EevDvXlh9nquxybw0SPB9GsUmDqQVHDh4NT3twtqRwGvAgT5BQHGKJLzj85nQ+8NuLu6ExUfha+7b+p7NdvLqIHP6Kap7kWfR9xMuEl+9/w0K9XsrutEXbtMwsx+XFW+5H9sKtXs3Linp0/VPrQt2zb1+fn3N77PhRsXmNFpRl75wsmoa7PO5zwuPCackX+NpEeVHrwW6wbXTpL42ALGborhp7X7qFTUh+lPNMqwL02AVwDtg9qnTj9U4SGqFKyCu6vx3+vVda8SmxjLLx1+ASAxORE3F13hzpYy+nSPiUhHpdTStDNF5EGM+2yak9sXvo+hK4fydauvaVC8QbrrxCcmc3xcP2okXWF/+1nUrVghh6O0XNriOB2COhAVH5XauH+x7QtaBbaibtG69grP1mJ1Pmt3UzhfYX5q+xOVAyqDuy9xgc0ZsCKZTSeO06dhIO90rk4+j8w9/la9UHWqF6qeOt2pfCeSkpMAo1Ns14VdaVu2Lc/V08Mf2EpGDfzzwBIR6YlRvAIgBGNEuM42jktzAAW9CnJ/6fupWjD9UX6TkxWvzN3N5ag2PN/wERo2bZvDEWZdu6B2qa/DY8JZcGwBxfIXo27RuiQmJ3It9hpFvB3vSkQ2nAW+0vmspZWskjkddZpy/uVoEHUV8l0n0c2XEWtd2XwynM971LZaudnO5f/7b5aYnEirwFap3y3RCdH0WtyL5+o9R9uyued7xNHd9QakUuooUBNYizFgRJD5upZZHENzciV9SvJJ80/w9Uj/stw3S7awYNd5mrbpRsNHRuRwdNZTOF9hVvVYRfdK3QFYF7qOtnPbsidsj50js6o4dD5rt5l1aBbdF3XnyOEFMLs/avmbvPH7XlYdvMQHD9ewWS15d1d3Xqz/YuqBdmRcJOX8yxHgGQDAiYgTPLXyKY5eO2qT/ecVGd4AUUrFAb/kUCyagzh45SAzDs3g5QYv4+fhl+46v69ay8BtfShT8Vm6teqYwxFaX9rH/6oUrMKTtZ5Mvby4+MRizt84z5DgIbl6ZDudz9rtHiz3ILE3w6i0+DXwLcZ3+Z5mzj+hPNe6Ev2bBOVYHCV8SvDNA9+kTl+JvcKFmxdSTy42nt/ImrNrGF5nOP6eethZS+WZLsSa5fZd2cemC5u420iDK3edoNq6Ebi5utClax+n66RWyqcUI+qMSO0AtOPSDtacXZPauJ+KPHXPYj+a5siiE6JRShHgmo/BOxcisVH8XmUsn/9zlb6NAnmhTSW7xtegeAMWPbKI4vmLA3A84jgrTq9ILVy18vRKZh+afdfvKM2guzBqd+hRuQedy3cmn1u+O5ZtP3WFuPkjqOwSSkLP33ArFJTzAeawd5q8Q2xiLGAMZ9v/z/60LduWd5q8Y+fINC3zklUyL619CW83b8Ym+iHntrOl4Te8sDaRDjWK82GXYIc7aH+s+mP0qdon9SB7xakVnI46Ta+qvQCjTHVR76LUL6bHQEsro0p244E/gVVKqes5F5JmL8cjjhObFEuNQjXSbdyPh93grykf8LLLv0Tf9wbeVfNOZxgvN6MUp4jwZuM3Ke1j3JsMjwnnf5v/x4g6I6hQwHGfIADKikhXdD7neYJwX6n78HD1QALbc0CVpd+aIjQuH8BXvevg6uJYjXuKtLfHPm3xKdcTjP/GSim+2v4VtYvUTm3gV59eTXDhYIrlT3cQwzwjo0v0E4HawFIRWS0ir5rDQWpO6qvtX/Hs6meJT4q/Y9nl67EMmLSFZHHlZsXOeD/wsh0itD93F3c6BHVILfBxMvIk2y9tTz3juXTzEpejHXJclXB0Pud5SilEhH7+1ekR1IldYck8ur44FYv6Mv7xELzcc0cfExFJ7R8kIvze5ffUwlaRcZG8uPZF5h6dCxhXLLZf2k5CcoLd4rUXseQehogUAtoBD2L0xN0JLFNKzbFteBkLCQlR27bpge2sJSo+ilORp6hVpNYt82/EJdJ73L8cD4tm1tDG1C7tDw52Cc+e0hbs+GjTR/xx/A/W9FqT7lUQexGR7UqpEPO1zuc86ErMFUasHsFrlftRZ94Iosq14/4jPfHxcmPesKYU9XWOAWOUUhyPOI6Phw/F8xdnb9he+i7ty5jmY+hYviPRCdFEJ0bfUhcjN0mby/diUSc7s670TKXU4+aYz98D9u2FoVlNZFwkSin8PPzuaNwTkpJ5dtpmRoW/zeyW16hdpoBu3G+TthrX49Uf572m76U27h9s/IBZh2bZK7R06XzOmyLjI0lMiiPfyndIcnFj4Mk2uLoIUwc3cprGHYwz+ooBFVM76FUoUIEvW36ZWo1z9ZnVtJrTimPXjgFGh8OUAjzOJku96JVS25VSH1s7GC3nJSUnMXz1cEatHXXHMqUUr8/fS8uTX9LSZRe1ijvPl4CtBPoF8mC5BwHjzP7cjXOExRgDtSml2Hpxq8N9meh8zhvK+5ZlTpRQ+cpZXnF5mSOxAUwe1JCgwvntHZpNebt706Zsm9TH6+oUqcOokFGUL1AegAl7J9D6t9aptyYdLT+zQ/eiz+NcxIVHKj6S7vPuX648QtLOmQzwWAlNnoEaXe0QYe7l5uLGuLbjUh+p2xO+h8HLB/NRs4/oUrGLnaPT8ootF7aw+eJmhl+NwPX4an7ye44/rgQxeVB9gkvlvWfKy/iVYUCNAanTIcVD8HT1TK2F8dr614hPiufrB3L/oKm6gc/jRCTd8dJnbD7Dqr9XsdBrIqrsfUib9+0QnXNIGbGuesHqfHb/ZzQv1RwwHvVZcmIJ7zd9nwJeBewYoebMNpzfwN9n/mZA409YsS+CMRcb812fOjStmDvvQVtb05JNaVqyaep0cOHgWzrkDV0xlMYlGzM4eLA9wsuWe16iFxFvEXlbRH42pyuJiK5dncslq2ReXvsy60LvGEGU1Qcv8daCvTxVZD9uPoWRRyeDqz4WzC53V6MHfn5345LojYQbhMeE4+dpXD3ZdXmXzXvg63zOe16s2ItpD07lo3+SGHXhAT54uAadapWwd1gOa0CNATxR8wkAEpITCPAKwMfdx5hOSmDoiqFsOLfBniFazJJv7V8wBqdoYk6fA34DFtsqKM32rsZe5Xjkca7GXr1l/q6zETwzYyc1SvrT5smvkcQI8HGqQVccRrdK3ehasSsiglKKt/55i6LeRZnUfpItd6vzOY/4/ejvNPYtR4lpPTns24Y5JzvleAna3M7dxZ0xLcakTl+OuUxEXASJyYkAnL9xnnF7xjE4eDBl/craK8y7sqSTXQWl1KdAAoBSKhrQ3ahzucL5CjO702y6VPjvXvCp8JsMnryVJ71WMaWzD/m93HXjbmMpz8+LCD+0/oFRIUZnx+iEaB5d9Ch/nfnL2rvU+ZwHRMRG8NnWT5m09AkSYq7z1qlaDlGCNrcr5VOKOQ/NoWWZloBRB2PlqZWp/Wz2he/jl32/cCP+hh2j/I8lDXy8iOQDFICIVMAYmUrLhZJVMnOPzCU+KR53V/fUBib8RhwDftlC8+StvJjwMwX3/2rnSPOeQL/A1AFursZeJcArILXn7+Xoy6wLXWeNHr46n/OAAp7+zHKvwAtnjzE85mnKV2/gkCVoc7tmpZqxtvdagvyCANh0YRPf7/o+tere1otbWXN2jd1q5lvSwL8LLAPKiMh0YDXwik2j0mxm68WtvL/xfVaeXpk6Lzo+kSGTt+IReYrP3X+AErWhvX5qyp5K+5bm53Y/p5beXHhsIc+sfoaL0RcBsvOFofPZiSUmJ7Lpwib452sC9//Bjwk9uB7UzqFL0OZ27i7/nSg9UfMJVj26KrUOxtQDUxm7bWzq8q0Xt3I66nSOxXbPe/BKqZUisgNojHEpb6RSKtzmkWk20ahEI6Y+OJXaRYwqpYlJyTwzYyfHzl1iU5EfcEtwg55Twd1xqrBpMLDGQOoUrUMpn1IAfLDpA2ITY/mk+SeZ2o7OZ+c298hcPt78MZ+WHU6UasOqwv2ZnYtK0DqDtE/EjL1/LBduXkidfuefd6hQoALftf4OgL1he6kUUCl1rAtrs7RrdCnA1Vy/hdkpaL5NItJsQilFeEw4RbyLUKdondR5by/cx1+HLrMgeBu+x47CY/MgwPE6i+R17q7uNCjeIHW6qHdR4hKzfGVd57OT6h7Ukbg4L95Y6I1v/qrMHdIQPy93e4eVZ3m4etzS+e7ndj8TkxgDGP1sBiwbQL9q/Xgp5CWUUlyJvWLVErr3bOBFZBJQC9gPpAyCrQD9hZCLLDi2gNFbRjOt4zQqBRgdbb776xgzt5xleMsK1GnTFk61hoqt7RypZolhtYdl6X06n53TlZgreCfG4z6pM1evt8TVpQO/Dm7oVCVonUFp39Kpr91d3fn2gW8p4WM8sngi8gSXbl6icKkcbOCBxkqp6lbbo2YXjUo0onfV3qlDmv627SyfrzzCs9Vu8mLzIuDmARXb2DlKLQfofHYySileXPMCSZcOMOnaKXaoQCYPdf4StLmdu4t7an18gACvAMr4lrHqPizpZLdRRPQXQi5X0qckL9R/ARdxYe2RMF6fv5fO5YQXL7+JzH/C3uFpOUfns5MREQbHufH45VDeSxzMsP798mQJ2tyuoFfB1HK51mJJA/8rxpfCYRHZIyJ7RWSPVaPQbOaP43/wxvo3iE6IBmDfuUiGTdtO9aJefO3yBRJ/E9r/z85RajlI57MTiUuKI2n7VO7fNY9zUS24r9eLugStlsqSS/QTgf7AXv67Z3dPIlIG48ukGMY9vvFKqa9FpCAwGwgCTgE9lVLXxHiO4GugIxANDFRK7bD8T9HSEx4TzoWbF/B09eTs1WgG/rKVAG8PZpX9A9fdW+HRX6BoNXuHqeWcLOWz5ngu3LhA/z/780BEOZom1ca94yd0rKlL0Gr/saSBD1NKLcrCthOBl5RSO0TEF9guIiuBgcBqpdRoEXkNeA14FXgQY0zqSkAj4Efzt5YNg4IH8Xj1x4mKSWLAL1uIT0xicauLeK+cZIwQF9zN3iFqOSur+aw5GHcXdzySAvn5zH14NX+JF5tWsHdImoOxpIHfKSIzgD9IU/HqXo/VKKUuABfM19dF5CDG4zldgJbmalOANRgNfBfgV2VU8NgkIgVEpIS5HS2TNpzbQIBnADUK1yAhCYZM2UrotRimDWlE8SKJcPM5aP2uvcPUcl6W8llzMPHRxP/ci+KXWlKnQT1eaFvF3hFpDsiSBj4fxhdBuzTzMvVYjYgEAXWBzUCxNI32RYxL+GA0/mfTvC3UnHdLAy8iQ4GhAIGBgZaGkKcopfh6x9d4uHowpf1Unpu5k51nI/ipR2UaBvoZI8O1+9DeYWr2ke181uzr7zN/MXfV24y+doD6ZbvytC5Bq92FJZXsBmVnByLiA8wDnldKRaX9j6iUUiKSqZqbSqnxwHiAkJAQ+xT4dXAiwsT2E4mMjeSDxQdYceAS73SqSvsDr8L+ZKOYjf5CyJOym8+a/e3bOJmIhDB+9+3PkMHDdAla7a7u2sCLyCtKqU9F5FvMgSnSUko9d6+Ni4g7RuM+Pc0lwEspl95FpASQMgD2OSDtQ4ClzXlaJoReD6WkT0n8PPyY/m8Yv248zZPNyzE4cTYcWwmdvtCNe95UDCA7+ZwREXEFtgHnlFKdRaQcMAsohDE8bX+lVHx29qHBiX/mMuLAIqp5NKfhi2N1CVotQxmdwR8wf2/LyobNXvETgYNKqS/SLFoEDABGm78Xppn/jIjMwuhcF6nvv2dOVHwUjy19jHZB7ajuMYAxyw7xUO2SvF7hNMwaA3X6Qchge4ep2Ues+TtL+WyBkcBBwM+cHgN8qZSaJSI/AUMwOs5qWTR643fIulV0cS1PrWFT8Mtn3WemNeeTUQP/HLBYKTUli9tuhvk4jojsMue9gdGwzxGRIcBpoKe5bCnGI3LHMB6T05cSM8nX3Zdn6j5DYnRpXv5tN43LF2Rsax9cJj5ljBDX6XN99p53FQXIRj7flYiUBjoBHwMvmgf3DwB9zVWmAO+hG/gsOx4ezowDc3HxqUGfR8dRtGBBe4ek5QIZNfDZqpaglNqAMVpVeu4oeG72nh+RnX3mZUopRITqPu3oOX0j5Qv7MK5/CJ6Rh6FgOej5qx4hLm+zdGCprPgKY8hZX3O6EBChlEo0p1M6zN5Bd5q9t8gbMRwZN4zAxB58NLgXgSV0IRvNMhklfQERuetD0vqxGscRFR/FUyue4rEqw/jgt3h8PN34ZVAD/PO5Q75gGLpGn7lrbrbIZxHpDFxWSm0XkZaZfb/uNJuxyJgYPv+lO+8l/EPFJq2pVEY37prlMmrg/YHOpH8Wrh+rcSARsRHEJyXy2dLTRMcV5bdhTSh5aApcOQodRoOrHi5SwxXb5HMz4GER6Qh4YdyD/xrjBMHNPIvXHWazIDEpmY+mjGCZ31kaFuhC5/ZP2zskLZfJqIE/rZTSPbJygaL5SuFy8QXOX45k8uD6VI3bByvehErtQHQvWw2AOFvks1LqdeB1APMMfpRSqp+I/AY8itGTPm1nWs0CSinGTZvBRxcW0L5YHdo8NcneIWm5UEYNvL6m6+BuxN9g5sGZ7NhXi60nr/FNn7o0LZIA4wZAgbLQ9SdwsWQ8IS0PyOl8fhWYJSIfATsxnqjRLPT2knUEn5xIdL5StBk0zyhOpWmZlNH/mv45FoWWJX+d/Ytvd33LjZMjeKNjGx6uURimdIb4G/D4QvDSQ0ZqqU7YegdKqTUYpadRSp0AGtp6n85owvoTzDv7GX+V9+DRhxZBvgB7h6TlUndt4JVS+3IyEC3zLp8P5vrxFxkQ0oAnm5eH8zvh4j7o8h0U00N+a7eIvfcqmr0t2BHKlWVjuL9yXwa0Lo9bQDl7h6TlYvq6Ty4UnRDNvF2H+GjJBTrUqMnbnasbtahL1YORu8CnqL1D1DQtk9YeCWPfwjd5y30BCVVCcC/V2N4habmcvkGbC7299nPG7B1KnbJufNW7Dq6X98H2yaCUbtw1LRfadTaCX2d9xuKgbcys3AL3Rk/YOyTNCdzzDF5EmmFUoSprri8YdWnK2zY0LT1HLl1nxcZK+AV0ZfLgVnglRsHsxyAxDmp01ffdtQzpfHY8xy7f4KNJc/mZKfycFEirTt/quhWaVVhyiX4i8ALGgBFJtg1Hy8j5iGgGTNqCl0thZvd9mAJebjBzKESeg0F/6sZds4TOZwdyMTKWAZO28AUT8fX255WeC8GvpL3D0pyEJQ18pFLqT5tHomUoKjaB7rPf4LrvNWZ2/ZIyBb3h70/g6AqjxnyZBvYOUcsddD47iMjoBAZM2kKk6yb+V6MYk0I+pZBu3DUrsqSB/1tEPsOodBWXMlMptcNmUWm3iE9M5ump27l2w4W2wWWpWaoghB+DdZ9C7b4QMsTeIWq5h85nBxATn8SQKVupfGU1jzxUh0M3L+Ef2MTeYWlOxpIGvpH5OyTNPIUxWpRmY8nJipfn7ubf41f4vMeLdKtnjtlRuCL0mwtlm+r7dVpm6Hy2s8SkZJ6ZsYPKoXP5n/tEcP0SHvja3mFpTuieDbxSqlVOBKKlb8zyQyw9PY8BLZvRvX5piL8Jlw9C6RCoeMegfJqWIZ3P9qWU4vX5e4k6vI5yJRewxL8hneoNsHdYmpOypBf9O+nNV0p9YP1wtLSm/HuKcesOUbTaPyhfF1Ad4I+RcGAhPLcL/NMdgVPT7krns319uvww/27fyR9+3/KCtx8Fqz0ILnq8CM02LLlEfzPNay+MEakO2iYcLcWyfRd574/9tKlahrE9F+LqIrB5HOz9DVq9pRt3Lat0PtvJxA0n+XHNMdYWnEBAUhITH5qDFK5k77A0J2bJJfrP006LyFhguc0i0th++iojZ+2katkIvu7djvye7nB6ozFCXOUHoflL9g5Ry6V0PtvHgp3n+HDxAdrVCGBemeY8XfItvItWtXdYmpPLSiU7b4zxnTUbOB52gyFTtlGs0HUueH/GzCNT4GY4/DYACgTqEeI0a9P5bGNrDl9m1G+76VXmGt3vi2Hq2WUc9C9i77C0PMCSe/B7MXrZArgCRYAPbRlUXnX5ulH0ws1FmNq/EzsjPGhVphW4+0LDoVClI+QrYO8wtVxM53PO2nnmGsOm7WBAwF7eDvsfXP+Z2t3+pHj+4vYOTcsDLLkH3znN60TgEuBpm3DyrhtxiQz6ZStXbsQz48kQyhbOT9nCj0BMhHHG3mKUvUPUnIPO5xxy7PINBk/eSqP8F+iZ/C37S9emRrWHKe7uZe/QtDwiw2u9IlIK4wj/glLqNJAAfAAczYHY8oyEpGSGT9/BoYvXeatrAK9t6ceesD2wezZ8Ww/CDts7RM0J6HzOORciYxgwaQsF5ToTPD7ns4L+vOTvSYKr7jGv5Zy7nsGLyPPAm8AxwFNEfgDGAL8C9XMkujxAKcVr8/ay7kgYY7rXpF65ONaEB1I65rrxSFypelBQjwOiZVtRYBc6n20uIjqeAZO2EBUTx/KSP+N2+TKfdZ3DRb9iuLu42zs8LQ/J6BL9UKCKUuqqiAQCR4BmSqntORNa3vDFyiPM2xHKyNaV6NUgEICfm38K41sa99t7TAZX/aWgZVsRoJTOZ9uKiU/iiSnbOBUezeTBDTh35WEqew0hoFxLAuwdnJbnZHSJPlYpdRVAKXUGOKy/DKxr+ubTfPvXMXqFlKFLAze+2fENcQkxMN8cIa7nVD2+u2YtyTqfbSulBO32M9f4oWsgpYtE0+fwBH7mmr1D0/KojM7gS4vIN2mmS6SdVko9Z7uwnN+qA5d4e8E+WlYpwkddg5l2cApzjsyhb8WueLp5wYOj9QhxmjV56Hy2nZQStKsPXWZ8i1jaLG+D8vuV1xq+Rtuybe0dnpZHZdTAv3zbtD7at5KdZ67xzMwd1Cjpz/d96+Hu6sKg4EE8XOFhCuUrBD1/tXeImvMJ5dYc1vlsRWOWHea37aG83cybtnufIdKvJP6BDenp5W/v0LQ87K4NvFJqSk4GklecDL/JkCnbKOrrxaSBDbgWf5Fr8VA6IZ5Ci56DLt9BQJC9w9SczxWd07YxYf0Jflp7nEEhhRl89nlm5nNlQpH8TE+MpgS6gdfsx5Ln4DUrCb8Rx8BftqCUYvKgBhTx9WTE6k84dOUgy67E4h55zt4hapqWCQt2nuOjJQfpFFyUdxK/QsIOUf+Rrzkde45i+YvZOzwtj9MNfA6Jjk9kyOStXIqKZcaTjSlfxAeAtxq+ycmlI3G/tNIY312fvWtarpBSgrZJ+UJ80asOal01pGxTqtR+nNftHZymkbVa9FomJSYlM2L6Dvaei+TbPvWoFxhAXFIcACUO/EHTgyug1ZtQqY2dI9U0zRIpJWgrF/NlfL9aeLi58ZLLFX701l+pmuPIqNDNt/xXs/oOutetZZRSvL1wH38fDuOjR4JpW70YSilGrB5BkG8gb+1ZrUeI03JCmdt60d9C57PlUkrQFvH1ZFrnfPj+3JiE7hPw8fAhv1t+e4enaakyukS/zfzdDKgOzDanewAHbBmUM/n2r2PM3HKWEa0q8FjjsgAkqSTqFa1HMe9iMGgUJCfqEeI0W4vG6Dmv8zkbUkrQurq4ML1PeQrO7QQqGfcCZfmwzIcodddzIk3LcffsRS8iw4D7lFKJ5vRPwPqcCS93m7PtLF+sPEK3eqUY1a5K6nw3pRgeEQXV+oOHtx0j1PKQK0qpKTqfsy6lBG1kTAKzh9SjzMrHiYwJ58N6nRjloigOiIi9w9S0VJacNgYAfmmmfcx5GRKRSSJyWUT2pZlXUERWishR83eAOV9E5BsROSYie0SkXmb/EEez5vBlXp+/l/sqFmZ0t1qICEopPtn8CXuXPANrPoGT6+wdppb3ZCmf87qY+CSGmCVox/evR43dH8GZfzna8mW2RBwlPCbc3iFq2h0saeBHAztFZLKITAF2AP+z4H2TgQ63zXsNWK2UqgSsNqcBHgQqmT9DgR8t2L7D2hsayfDpO6hSzJcfH6uHh5vxMV+Ovszq44vZe3QxNB4OVTvZOVItD8pqPudZKSVod5y5xle969A0yB9uhMF9LxLS9CWWdV9GcOFge4epaXfI8DE5EXEBDgONzB+AV5VSF++1YaXUOhEJum12F6Cl+XoKsAZ41Zz/qzJuYG0SkQIiUkIpdcHCv8NhnL0azaDJWwnw9uCXQQ3w9fpvoJhi1y+z4OQJvEvWgbYf2C9ILU/KTj5nsM0yGCPSFcPolDteKfW1iBTEuM8fBJwCeiqlcl1R9rQlaD96JJiONUsAsKvVKC7HhNEO8HbXt9k0x5ThGbxSKhn4Xil1USm10PzJ8pcBUCxNo30R40sBoBRwNs16oea8O4jIUBHZJiLbwsLCshGK9V29adyjS0hKZsrgBhTz8wKML4k1Z/4mefHz+Hj54/LoZD1CnJbjbJDPAInAS0qp6kBjYISIVOfuV+tylZQStM+3qcRjlZPh1y4QcZZfD03jq51fpz7uqmmOyJJL9KtFpLtYufeIebae6S6nSqnxSqkQpVRIkSJFrBlSthjDRG4lNCKGCQNCqFjUN3XZxvMbefbv51jWeCD0mQm+usKVZjdWzWel1AWl1A7z9XXgIMbBeReMq3SYvx+xxv5yUkoJ2n6NAhl5X3GY1RfO74KkeMY0H8PP7X7G09XT3mFq2l1ZUsnuKeBFIFFEYgHBaJ/9Mn5bui6lXHoXkRLAZXP+OaBMmvVKm/NyhaRkxchZO9l5NoIf+tajQVDBW5Y3iU/my/s/p1Vga3BxtVOUmgZYN59vYd6Sqwts5u5X63KFlBK0HWsW54OHqyNz+kPYYdZ3/JAG/qXwcnWnlE+6Fxk1zWHc8wxeKeWrlHJRSnkopfzM6ax+GSwCBpivBwAL08x/3OxN3xiIzC3335VSvLdoPysOXOKdztV50LxHlyL64B/IlM60uXgcV924a3Zm5XxOJSI+wDzgeaVU1G37vOvVOke85Za2BO2XverguvYTOLyECw+8znMHJ/Dj7lzdB1jLQyyqRW8+zlYJ8EqZp5TK8BkvEZmJ0aGusIiEAu9i9OCdIyJDgNNAT3P1pUBH4BhGQY5Bmfor7OjHtceZuuk0Q1uUZ1CzcrcsW75vKqO3jGZSiWqUqz/QPgFq2m2yks/32J47RuM+XSk135x9t6t1t1BKjQfGA4SEhNi9SswtJWgfr4+niodDS6Buf0rcN4ofKtxHzcI17R2mplnkng28iDwBjMS4bL4LoyPNRuCBjN6nlOpzl0Wt01lXASPuFYujmb8jlE+XHeah2iV5rUPVWxfGR1N2ww80U4mU6TkN3PPZJ0hNSyOr+ZzB9gSYCBxUSn2RZlHK1brR3Hq1zmGlLUE7eXDKEzDuJA9extnoy5QVoUnJJvYOU9MsZkknu5FAA+C0UqoVxj22CFsGlRtsOBrOK3P30KR8Icb2qIWLy219lpaOourFg3zU5nvcClWwT5Cadidr53MzoD/wgIjsMn86YjTsbUXkKNDGnHZYFyJjeHziZlxdXJg6pCFFJQqWvgLx0Uw9voDuS/twKvKUvcPUtEyx5BJ9rFIqVkQQEU+l1CERqXLvtzmv/ecjeXradioU8eGn/vXxdLv13vpfZ/5ij39+hrf7GA89QpzmWKyaz0qpDRgd9dJzx9U6RxQRHc/jE7cQFZvIrKGNKevvDlP6w4XdUO9xOpXvRLJKpqxfWXuHqmmZYkkDHyoiBYAFwEoRuYZx/zxPCr0WzaBftuLr5cbkwQ3wz3fb8+yxkewO283GG6cZcb8uEKY5HJ3PaaSUoD19JZrJgxsQXNIPFj0LZzcR3XUc+YrVoLAIg4JzTbcgTUt1zwZeKdXVfPmeiPwN+APLbBqVg4qMTmDgL1uJSUhi7tNNKeF/23316xdhfEteaPY8Tz84BXcXXcxGcyw6n/+TkKYE7fd969G0QmHYPA52TiXpvhd57uIqil3fy8f3fWzvUDUtSzIaD75gOrP3mr99gKs2ichBxSYk8eSv2zhjHulXKe576wpJCez6rQ+FEq5TplwL8rnpTnWaQ3FNJ6fzbD6nW4I27jqsGwtVOiGt3qTpwV8p6JXe16Cm5Q4ZncFvx3h2VYBA4Jr5ugBwBih313c6meRkxUtzdrPl1FW+6VPXONK/fZ3lb/Je0nnyV6jJtKLV7npTUtPspDqwDZ3PgFGCdm5KCdrG5r11T194YiUqX0FcXN0YHDzYvkFqWjbdtRe9UqqcUqo8sAp4SClVWClVCOgMrMipAB3BxA0nWbL3Aq8/WJWHa5e8c4U9v+GyZRw/lXyQ99t8p8eE1hzRXp3Phin/nuKntcd5rHEgI1tXgtgo2PIzKMVlT2/6rBrK/vD99g5T07LNksfkGiullqZMKKX+BJraLiTHsv98JJ8uP0S76sUY2qJ8uutExkVA+VYU7/AZFQMq5myAmpY5eTqfj4fd4OOlB2ldtSjvPxyMqGSY/yT8+Spc2k9EXARxSXF4uXnde2Oa5uAs6UV/XkTeAqaZ0/2A87YLyXHExCcxctYuArw9GN29Vrpn5jcTbtLj9G90Du7Ec3qEOM3x5dl8Tk5WvDp3D/ncXfmke01cXQRWfQRHlkHHsVA8mMrAvIfn4SKWnPtommOz5H9xH6AI8Lv5U9Sc5/T+t/Qgxy7f4POetSmY3+PWhcnJ8NsgPA4upmvFrrQofb99gtS0zMmz+Tx102m2nb7G252rU9TXC/bOhQ1fQP2B/FOyGuP3jCdZJevGXXMaljwmdxWj+lWesvrgJaZuOs2Q+8rRvFI6w9Ku+wz2z8e9bFOGNRyW8wFqWhbk1Xw+ezWaMcsOcX/lInSvVwpirsEfz0NgU3jwM9ZuH8vWi1t5vPrj+vK85jQsqUVfGRgFBKVdXymVpdrVuUHY9ThembuHqsV9ebl9OkW+jq4keu1oXqlYm2HlmlAj50PUtCzJi/mslOKN3/ciwMddg41bbfkCoM9MKFIV3Dx4veHrRMVH6cZdcyqW3IP/DfgJmAAk2TYc+1NK8crc3dyIS2Tm0MZ4ud82xOvVkzBvCOeKVeWIpycxSbH2CVTTsiZP5TPAb9tDWX80nA+61KC0awTsXw41ukK55sw/Op/7St1HUe+i+Hv62ztUTbMqSxr4RKVUnhkAeeqm0/x9OIz3H65B5WK+d65w8A9AqNRzBkv8S+OuO9ZpuUueyufLUbF8tPgADYMK8lj5GJjQCeKvQ7n7CZdkRm8ZTa8qvXgp5CV7h6ppVmdJA/+HiAzH6JATlzLTvJfnVI5cus7HSw7SqkoRHm+S/sASsY2Gsszfn4cDyuKuO+NouU+eyWelFG8v3EdcYjJfNYnG5Zde4OYFA/4A74IUBmZ3nk1Jn3RqW2iaE7CkgR9g/n45zTwFpP9QeC4Vl5jEczN34uPpxqeP1r7zkbjds6BoNRZFHebD7WOpUKweNYvUtE+wmpZ1eSKfAZbuvcjy/ZcYV/8cJRe9BQXKwGPzSPAvxa6LW2lQvAHl/PNUAT8tj7GkF32eyIDPlh3m0MXrTBwQQhFfT2NmbCQcWgL75sPx1VDtIXr0mEKlgEq6cddypbySz9duxvPuon3ULOVPm6IHIbI29JkF+Qux7sxqnv/7eeY+NJcqBfP0yNeak7PkDB4RCcaoZZ3axVQp9autgsppG46GM2HDSfo3LkvrasWMmevGwtoxkBQP/oHENRlObKNh+ItQt2hd+wasadng7PkM8MEf+/GNDmXM4O7cKFiT1+OP0fHSRjqX70yj4o0YFTJKN+6a07vnTWQReRf41vxpBXwKPGzjuHLMtZvxvDF7E08E7OS9mE8g4oyxoFgNaPAkPLEant/DuMLFeGR5fyJiI+war6Zlh7PnM8Ca/aF4nX2e4YU/pLrPDfw8/YhOiiUhKQEAHw8fBtQYcI+taFruZ0kvsUeB1sBFpdQgoDbGGNJOYeG8X/kzYTBvxXzGqOs7mXfAqOCpKneg7Y2tTLi2C0RoF9SOvlX7UsCrgF3j1bRscup8vh55lXxz+xLme541pSuDbwlEhMkdJtO1Uld7h6dpOcqSS/QxSqlkEUkUET/gMlDGxnHliBMXr+B54TMivYqRv/cPXDs0kegCxp8mIrQo1YJyfsYty6oFq1K1YFV7hqtp1uC0+bz13CbiZo2kcfIR/Cp8SLW2z9g7JE2zK0sa+G0iUgD4GWOM+BvARlsGlVOmLvqSucU9iCvTjYHlmjO5XPNblr/d5G07RaZpNuOU+ayU4sM1r5Hf9xqXC79NN924a5pFveiHmy9/EpFlgJ9Sao9tw7K9XWcj+OVEYwbXL0Tf+0fYOxxNyxHOms/xScnEXHqe8gkbaD/yWXuHo2kOwZJOdqtTXiulTiml9qSdlxsppfhu8T8U9vHkxYefxcPV495v0jQn4Iz5nJCcwIS/D3H0sjtdurxCfi9dXVLTIIMGXkS8RKQgUFhEAkSkoPkTBJTKsQht4O9NK7ngNorHK6/Cx9OiJwU1LbcTZ83nd5c+w+6jjzK0yg1aVS1q73A0zWFk1Lo9BTwPlMS4V5dS2i0K+M62YdlOUrIifs1HlCiQTKvGj9g7HE3LKUUw8tip8lklxlPi2L/kF8UTvdvaOxxNcyh3beCVUl8DX4vIs0qpb3MwJpvatGIWHeJ2U6bIq9QIbGLvcDQtp1xWSoU4Wz7vmTeaZyPOsLb+txQtWNDe4WiaQ8noEn0DESme8mUgIo+LyEIR+ca81JfrxMbFsXf/xxxyK0m1Lnr0KC1P8Xa2fF5/8E+unBrPNs+GNO/U397haJrDyaiT3TggHkBEWgCjgV+BSGC87UOzvhnLfuXHQrAwuC0u7p72DkfTclJZnCyfv/93LJ8U8sW/+1hcXOTeb9C0PCaje/CuaYaQ7AWMV0rNA+aJyC6bR2ZlUbEJfLO7LA1KPM/znR+zdzialtPEmfJ5/dEwNh17lmcbRFGpcm17h6NpDimjM3hXEUk5AGgN/JVmWa7rej7tz5VERCfwQqceeLp53fsNmuZcxFnyOfJGJF/NW035QgUY8VA/e4ejaQ4ro8SeCawVkXAgBlgPICIVMS7r5Rqnju9g3tVR9C5fhxolO9k7HE2zhys4ST5/MmcwSYX282HIFLzcXe0djqY5rIx60X9sFsAoAaxQSilzkQuQq0pFXVj0Nk3d4niglR5BSsuzLgIvkcvz+fjRQ7QK20yyfyDN6jaydzia5tAyvDSnlNqUzrwjtgpGRDoAXwOuwASl1OjsbvPwluU0idxAXMknaRncPtsxalpuldvzWSlFxPznaRETT+P+U60So6Y5M0uGi80RIuIKfA88CFQH+ohI9exsMyEhgVmbXmGnWyFC+r5jjTA1TbOALfJ57tJPOe62lz0Vn8a/ZAVrhKlpTs2ROtc0BI4ppU4AiMgsoAtwIKsbnLjid5b6JuFdsgN1ffysFKamaRawaj4nJCWz8NRKzhcoyIJuuoaFplnCkRr4UsDZNNOhwB032URkKDAUIDAwMMMNtm/UiZubi/Bc+6ZWDFPTNAtYNZ/jE5MJKvkl3QPj8fP2t3KomuacHKmBt4hSajxmYY6QkBCV0brlCufnpU6tciQuTdMyz9J8zu/pxkdda+ZYXJrmDBzmHjxwDiiTZrq0OU/TtNxH57Om2ZkjNfBbgUoiUk5EPIDewCI7x6RpWtbofNY0O3OYS/RKqUQReQZYjvFYzSSl1H47h6VpWhbofNY0+3OYBh5AKbUUWGrvODRNyz6dz5pmX450iV7TNE3TNCvRDbymaZqmOSHdwGuapmmaE9INvKZpmqY5IflvUKncR0TCgNP3WK0wEJ4D4WSGjskyOibLWBJTWaVUkZwIJqt0PluNo8UDOiZLWTWXc3UDbwkR2aaUCrF3HGnpmCyjY7KMI8ZkK474tzpaTI4WD+iYLGXtmPQlek3TNE1zQrqB1zRN0zQnlBca+PH2DiAdOibL6Jgs44gx2Yoj/q2OFpOjxQM6JktZNSanvwevaZqmaXlRXjiD1zRN07Q8RzfwmqZpmuaEnLqBF5EOInJYRI6JyGs5tM8yIvK3iBwQkf0iMtKcX1BEVorIUfN3gDlfROQbM8Y9IlLPhrG5ishOEVlsTpcTkc3mvmebw3oiIp7m9DFzeZCN4ikgInNF5JCIHBSRJvb+nETkBfPfbZ+IzBQRr5z+nERkkohcFpF9aeZl+nMRkQHm+kdFZIA1YrMXe+SyuV+dz5bHo/M5/Rjsl89KKaf8wRii8jhQHvAAdgPVc2C/JYB65mtf4AhQHfgUeM2c/xowxnzdEfgTEKAxsNmGsb0IzAAWm9NzgN7m65+AYebr4cBP5uvewGwbxTMFeMJ87QEUsOfnBJQCTgL50nw+A3P6cwJaAPWAfWnmZepzAQoCJ8zfAebrAFv//7fR/xO75LK5b53Plsej8zn9OOyWzzZPEHv9AE2A5WmmXwdet0McC4G2wGGghDmvBHDYfD0O6JNm/dT1rBxHaWA18ACw2PwPFA643f55YYzh3cR87WauJ1aOx99MPrltvt0+J/ML4ayZRG7m59TeHp8TEHTbF0KmPhegDzAuzfxb1stNP46Sy+a+dT6nH4/O54xjsUs+O/Ml+pR/3BSh5rwcY17iqQtsBooppS6Yiy4CxczXORXnV8ArQLI5XQiIUEolprPf1JjM5ZHm+tZUDggDfjEvM04QkfzY8XNSSp0DxgJngAsYf/d27Ps5pcjs52L3//9W5BB/i87nDOl8zpwcyWdnbuDtSkR8gHnA80qpqLTLlHEIlmPPJ4pIZ+CyUmp7Tu3TAm4Yl61+VErVBW5iXKpKZYfPKQDogvFlVRLID3TIqf1bKqc/F03nswV0PmeRLT8XZ27gzwFl0kyXNufZnIi4Y3wZTFdKzTdnXxKREubyEsDlHIyzGfCwiJwCZmFc1vsaKCAibunsNzUmc7k/cMXKMYUCoUqpzeb0XIwvCHt+Tm2Ak0qpMKVUAjAf47Oz5+eUIrOfi93+/9uAXf8Wnc8W0fmcOTmSz87cwG8FKpk9Jj0wOk0ssvVORUSAicBBpdQXaRYtAlJ6Pg7AuJeXMv9xs/dkYyAyzaUbq1BKva6UKq2UCsL4HP5SSvUD/gYevUtMKbE+aq5v1SNMpdRF4KyIVDFntQYOYMfPCeNSXmMR8Tb/HVNistvnlEZmP5flQDsRCTDPZNqZ83Iju+Qy6HzOREw6nzMnZ/LZmp0aHO0Ho0fiEYweuG/m0D7vw7jcsgfYZf50xLiXsxo4CqwCCprrC/C9GeNeIMTG8bXkv1635YEtwDHgN8DTnO9lTh8zl5e3USx1gG3mZ7UAo3eoXT8n4H3gELAPmAp45vTnBMzEuGeYgHFmNCQrnwsw2IztGDAoJ/7/2/D/bY7nsrlfnc+Wx6LzOf0Y7JbPulStpmmapjkhZ75Er2mapml5lm7gNU3TNM0J6QZe0zRN05yQbuA1TdM0zQnpBl7TNE3TnJBu4G1MRJSITEsz7SYiYWKOAJWF7RUQkeFZeJ+PiIwTkeMisl1E1ohIo0xuY42IhGRi/c5m2crdYozG9VRm485kfO+JyKhMrF9MRBaniW+pOb+kiMy1XaRabqRzWedybuN271W0bLoJBItIPqVUDMZAFdmp2FQAY9SjHzL5vgkYg0FUUkoli0g5jFGxLCIirpnZmRjVv8YDDZVSoSLiiTHggiP5AFiplPoaQERqASilzvNfIQxNS6FzWedyrqLP4HPGUqCT+boPRuEDIHVc4AVijP27KeU/pnkEO8k80j4hIs+ZbxkNVBCRXSLymbnuyyKy1dzG+7fvXEQqAI2At5RSyQBKqZNKqSXm8gXmmcB+ERma5n03RORzEdmNMepS2m32EZG9YoyzPCadv9kX4wDyirm/OKXUYfO9D4kx3vJOEVklIsXS/M1TRGS9iJwWkW4i8qm5n2XmFw0icirN/C0iUjG9v9l8z3Zze1XTibEERuEJzBj3mO8NEnPsZjEGzdhl/oSJyLuWfOaa09K5rHM597BllSX9owBuALUwajN7YVTCasl/1ae+Bd41Xz8A7DJfvwf8i1F5qTBGcrlz57CD7TCOrgXjgG0x0OK2GB4Gfs8gxpQqSvkwKj4VMqcV0DPNemuAEIyBG84ARTAS/y/gkXS2OwGjxvJMoB/gYs4PgNQiS08An6f5mzeYf2dtIBp40Fz2e8o+gFOY1cyAx9N8lu8Bo8zXqzHOcMD4QvwrnfjaAxEYpSvfBEqa82/5jM15ZYGD5u97fub6x/l+dC7rXM5tP/oSfQ5QSu0RY6jJPhhnAGndB3Q31/tLRAqJiJ+5bIlSKg6IE5HL/DekYFrtzJ+d5rQPUAlYl4kQnxORrubrMub7rwBJGINs3K4BsEYpFQYgItOBFhjlKVMppZ4QkZoYgz6MwrikORBjoITZYgyy4IFxuTHFn0qpBBHZC7gCy8z5e7n1suDMNL+/TLtfMUb+agr8JiIpsz1v/yOUUstFpDzGCFMPAjtFJPj29UQkpYTls0qp0yLyLNn/zLVcSOeyzuXcRDfwOWcRxtjELbF8jOG4NK+TSP/fS4BPlFLjMtjOfqC2iLgqpZJuebNIS4ykbaKUihaRNRhnJwCxt6+fWUqpvcBeEZmKkfwDMc50vlBKLTL3/16at8SZ70sWkQRlHnJjjHud9u9Xd3kNxpF4hFKqjgXxXQVmADPE6CzVAmPM6LR+AuYrpVaZ05Z85prz0rmsczlX0Pfgc84k4H0zSdJaj3HJKyVBw9Vt403f5jrGPbEUy4HB5pEuIlJKRIqmfYNS6jjGIBDvi3kYbN6b6oQxJOI18wuhKtDYgr9lC3C/iBQWo8NOH2Bt2hXE6OnbMs2sOsBp87U//3VOGkDW9Erze2PaBebnd1JEepixiIjUvn0DIvKAiHibr32BChiXK9OuMwLwVUqNTjP7np+55tR0LutczhX0GXwOUUqFAt+ks+g9YJKI7MG4T5VhkiilrojIP2bHkT+VUi+LSDVgo5nvN4DH+G984RRPAJ8Dx0QkBggHXsYY+elpETkIHAY2WfC3XBCR1zDudwnG5ceFt60mwCsiMg6IweiBPDDN3/ybiFzDuOdX7l77TEeA+ZnFYXwp3a4f8KOIvIVxH3AWsPu2deoD34lIIsbB7gSl1FbzEmyKUUCCiOwyp39SSv1k4WeuOSGdyzqXcws9mpyW64jIKYxhFMPtHYumaVmnc9m29CV6TdM0TXNC+gxe0zRN05yQPoPXNE3TNCekG3hN0zRNc0K6gdc0TdM0J6QbeE3TNE1zQrqB1zRN0zQn9H+36LmKUPEorgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Simulation Visualizations\n",
    "#Copula Model Speedup\n",
    "## Plot the computational speedup results\n",
    "fig, axes = plt.subplots(ncols=2, figsize=(8, 4))\n",
    "\n",
    "f1 = (res.k == 5) & (res.p == 5)\n",
    "f2 = (res.k == 10) & (res.p == 10)\n",
    "f3 = (res.k == 20) & (res.p == 20)\n",
    "axes[0].plot(res['nsamps'][f1], res['speedup_total'][f1], color='tab:blue', label='k = p = 5')\n",
    "axes[0].plot(res['nsamps'][f2], res['speedup_total'][f2], color='tab:orange', label='k = p = 10', linestyle='dashed')\n",
    "axes[0].plot(res['nsamps'][f3], res['speedup_total'][f3], color='tab:green', label='k = p = 20', linestyle='dotted')\n",
    "axes[0].set_ylabel('Standard Run Time / Copula Run Time')\n",
    "axes[0].set_xlabel('Monte Carlo Sample Size')\n",
    "\n",
    "axes[0].legend()\n",
    "\n",
    "axes[1].plot(res['nsamps'][f1], res['speedup_interp'][f1], color='tab:blue', label='k = p = 5')\n",
    "axes[1].plot(res['nsamps'][f2], res['speedup_interp'][f2], color='tab:orange', label='k = p = 10', linestyle='dashed')\n",
    "axes[1].plot(res['nsamps'][f3], res['speedup_interp'][f3], color='tab:green', label='k = p = 20', linestyle='dotted')\n",
    "axes[1].set_ylabel('Standard Run Time / Copula Run Time')\n",
    "axes[1].set_xlabel('Monte Carlo Sample Size')\n",
    "\n",
    "axes[1].legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABA00lEQVR4nO3de5zWc/7/8cdLihBTSmqSkrBOiw2R3XXYXexa5dhMOS5a+7WWPaADCp3wc9hdh90cVhYzSjWiyKmwDlGKio2i1IgKYQnJ6/fH5zPXXDNzXddczVzn63m/3eY21/W5Tq8u9fL8vN/vz+dj7o6IiIhIIdks2wWIiIiIpJoCjoiIiBQcBRwREREpOAo4IiIiUnAUcERERKTgbJ7tApJxzDHH+OOPP57tMkSkcZbtAlJFfUckb8TsO3kxgrN27dpslyAiRUZ9RyS/5UXAEREREdkUCjgiIiJScBRwREREpOAo4IiIiEjBUcARERGRgqOAIyIiIgVHAUdEYqqaV02fsc/QffA0+ox9hqp51dkuKcLM7jaz1Wa2sN72C83sv2a2yMyui9o+xMyWmNliMzs68xWLSKblxYn+RCSzquZVM2TyAtZv2AhA9br1DJm8AIB++5dms7Qa9wC3APfWbDCzI4C+wA/d/Rsz2yHcvidQBuwFdAaeMrPd3H1jxqsWkYSq5lVz/YzFfLBuPZ1LWnPJ0bs3uedoBEdEGrh+xmLWb9jI8W8+y4gn/4H596zfsJHrZyzOdmkAuPtzwCf1Nv8OGOvu34TPWR1u7wtUuvs37v4esAQ4KGPFikhSquZVM2TSG5w+5VZOef2JyI5VU0ePFXBEpIEP1q3n+Def5aZHb2CPNcto9d2GyPYcthvwYzObbWbPmtmB4fZSYEXU81aG2xows0FmNsfM5qxZsybN5YpItKumLuTiJ+/k/Fcm84M17wE0a8dKAUdEGjhj2Uvc9OgNvNplT84+eQTftNwCgM4lrbNcWUKbA+2A3sAlwAQz26RrY7n7OHfv5e69OnTokI4aRSSGqtdWcv70cfz2lcmMP+BXXHXUoMhjTd2xUsARkboqKxkxcQxzdtqLs08ewfpWWwLQsoVxydG7Z7m4hFYCkz3wCvA90B6oBnaKel6XcJuI5AJ3vvzjnyPhZvjPzoeofZOm7lgp4IhIrcpKGDiQj/c7kPNOqQ03AHj2ykpSFXAEgJntBrQC1gJTgTIz28LMugM9gVeyVaSIRHGHyy5j4HMTuHf/huEGaPKOlQKOiAQefBAGDoTDDqN/3yv4fPMt6jy84XvPmUXGZlYBvATsbmYrzewc4G5gl/DQ8UrgzHA0ZxEwAXgTeBy4QEdQieQAdxg8GK6/nkm9j+fKnzcMNyWtWzb5KCodJi4iQbgZMAAOOwymTePdkc/GfFquLDJ29/I4D50W5/mjgFHpq0hENklNuLnuOvjd72hx7lBaT1kYOTUFQOuWLRhx/F5N/ggFHJFiVzNy06cPTJsG22xD55LWVMcIMzm+yFhE8oE7DBkSCTfceiv9zMAsZefAAQUckeI2YQIMHMjafXtRdtRlLB35LJ1LWnPEHh2YNLe6wd5Uji8yFpFcVxNurr02CDe33BKZluq3f2lKTySqNTgixWrCBBgwgLX79uLon1/GkvXBOuLqdeuZNLeak35USmlJawwoLWnNmBP3yZWzGItIPnKHoUPrhpvN0hdDNIIjUozCcMMhh1D2s8F8XG82av2Gjcz87xpeGHxkduoTkcJSE27GjoXzz4dbbqHq9VVcP2Mx1evW08KMje6UpmBqqoZGcESKzcSJkXDDY4+xNM664VxZUCwiea5+uLn1VqpeX8WQyQsia/02enAeiuZeniGaAo5IMZk4EcrLI+GmZkFxLFpQLCLN5g7DhgXh5re/hVtvhc02i1zvLpZUXfdOAUekWESHm+nTYZttgOAkWq1btqjzVC0oFpFmqwk3Y8YE4ea226h6fRV9xj4T8yjNaKkYQdYaHJFi8NBDQbjp3TsIN23aRB6qmetO5eGZIlLk3OHyy4NwM2hQJNwMmbwg7shNtFSMICvgiBS6hx6CsrIg3Dz2WJ1wUyPVh2eKSBGrCTejRwfh5vbbYbPNuOqRRUmFm1SNICvgiBSySZMaDTciIinjDldcEYSb886LhJuqedV8+tWGuC9Lx1FUCjgihWrSJOjfn+U99+HXB13M56Oeo4UZ5QfvxMh++9R5atW8ak1RiUjz1ISbUaOCcPOPf0TOc5No0XBpSeu0nJJCAUekEIUjN8t77sMvjxnKl62C+eyN7tz38vsAkZBTNa+6zrx4zWGagEKOiCTHHa68Mma4gcSLhtN1QIOOohIpNJMnB9NSBx3E8ccO48sttmrwlIrZKyK3Yx2umarDNEWkCNSEm5Ej4dxzG4QbiL9ouDlXC2+MAo5IIZk8Gfr3h4MOgsce47NWsZtKzUm1IP6elU70JyKNcofhw2vDzT//GfPyC/FOR9Gcq4U3RlNUIoWiJtwceGCwoHjbbSML9+prEV7cDtCVw0WkaWrCzTXXRMJNzeUX6q/ny8bpKNIecMysBTAHqHb348ysO1AJbA/MBU5392/TXYdIQZsyBfr35+M9f0jZL4awZPTzdC5pTe9d2vLC0k8aPL384J0ity85evcG56bQif5EJKHocHPOOZFwk2g9X6ZPR5GJKaqLgLei7l8L3OTuuwKfAudkoAaRwjVlCpx6Kh/v+UOOPnoI73y9WeSq4K+9/xl9erSLjNi0MOO03l3rHEXVb/9Sxpy4j64cLiLJcYcRI2rDzbhxcS+/kM31fGkdwTGzLsCvgFHAn8zMgCOBAeFTxgMjgNvTWYdIwQrDDb16UXb0UNZ+XXefZf2GjSz7eD1Lx/wy4dvoRH8ikrQRI+Dqq+E3v4mEG8i99XzpHsG5GbgU+D68vz2wzt2/C++vBGJ2VTMbZGZzzGzOmjVr0lymSB6qqoqEG2bMYMnXsf85F+JiYTO728xWm9nCGI/92czczNqH983M/mZmS8zsDTM7IPMVixSI6HBzxx2Ra0t1HzyNzaLW9kXL1nq+tAUcMzsOWO3uc5vyencf5+693L1Xhw4dUlydSJ6rqoJTTgnCzeOPw7bbFttVwe8Bjqm/0cx2An4BvB+1+VigZ/gzCI0YizTNiBFw1VVw9tmRcDNk8gKq163HIeYBDdlcz5fOEZw+wPFmtoxgUfGRwF+BEjOrmRrrAlSnsQaRwvPww0G4+dGPgnCz3XZAsFi4ZYu6e1AtW1hBLhZ29+eAhqun4SaCUePoTtsXuNcDLxP0oE4ZKFOkcESHmzvvjLvmBoK1frmwni9ta3DcfQgwBMDMDgf+4u4DzWwicDJB6DkTeDhdNYgUnIcfhpNPDsLNjBmRcBNRfweq4Q5VwTKzvgRHa75udYfKS4EVUfdrpsZXxXiPQQSjPHTt2jV9xYrkk6uuioSbqguu4vrrZvFBOGoTy/fuvDf2VxktMZZsnOjvMoIFx0sI1uTclYUaRPJP9MhNjHBz/YzFbPi+bsvZ8L0XxRmJzWwrYChwZXPeR1PjIrWq5lVz58/OghEjePRHx3DFcRcxpGpRZEoqnlyZFs/Iif7cfRYwK7z9LnBQJj5XpGDUhJv99489ckPuHcGQYT2A7kDN6E0X4DUzO4hgGnynqOdqalykEVXzqllx8WAufO4+Ju79My496v/glZWNDgrn0jm0dCZjkVw3dWptuHniiUi4qX8F8JKtWvLpVxsavDxX9qbSyd0XADvU3A/X/vVy97VmNhX4vZlVAgcDn7l7g+kpEam19tJhkXBz2bEX4pZ4wscgI2cn3hQKOCK5bOrUYM1NjHBT/4yh8WzVqvAuOWdmFcDhQHszWwkMd/d4093TgV8CS4CvgLMzUqRIHqi/o3TEHh3odusNnPvUeB7a+yguO/ZCvt+sRcL3KC1pzQuDj8xQxclTwBHJVY88Uhtu6k1LxTt6IZZ3Vn+Zrgqzxt3LG3m8W9RtBy5Id00i+WbgHS/VuZRL9br1tL3hWs79z/08tPdRXHrsHxqEG6PusQu5NCVVX+Ht2okUgkcegZNOgv32C8JNSUmdh4tkXY2IpMnlVQsaXKfuwhcq+PN/7mfS3kfGDDetW7ZgYO+ueXNZF43giOSa6HDzxBMNwg3EvwK4iEgyKmavqHP/9y9WRsLNJcde1CDclObY+ppkKOCIZNnlVQuomL2Cje78bOmr/KNqFMu79OSEQ//M52NfAKBPj3Z077BN5HlGw6HiePr0aJfO8kUkj1TNq2bE1EV1zjr8+xcr+cvz9zFpryPihptcXGPTGAUckSy6vGoB970cXFXgyCWvcNuU0SzaoTun9x3O51tuE3neC0s/qTOcHC/YdGzTio+++DZyv0+Pdtx/3iFpqV1E8kvVvGr+NGE+0afLuuDFB2vDzS8vjjktlatrbBqjgCOSRTXDxEcueYXbq0bz1g7dOb3/NXXCzaZY+78NLMuBM4iKSO656pFFDcLNJc//u0642QzoVNI6clRVvk1LRVPAEcmije4csfRVbq8azX87NC/c1LyfiEgs0efJ+r+XJjQINxCMDufjdFQsCjgiWXTU0le5bcqolIQbCC5yJyISrf7h4P/30gQufe5eJseYliqkE4Mq4Ihky7Rp/LNqNG926JaScANQfvBOjT9JRIrGz2+cVedcWNHh5i/1wk0+r7eJRefBEcmG6dPhxBN5d8funNZ/ZINws2WLuiMxfXq047TeXSMjNC3M6LnD1nXun9a7KyP77ZOZ+kUk511etSBmuJmy5+ENwk2un9OmKTSCI5Jp06fDCSfAPvtw8mGXxBy5+Xqjx1wsrAAjIsmomlfN/eERmlA33Pz5V39scLRUoay7iaaAI5JJNeFm773hySf5/NoXs12RiBSQmmtLRZ8I9HcvT+TS5+6las+fxgw3pQW07iaaAo5Ipjz2WG24eeopaNs22xWJSAGJPq9Wjd+9PJHLnh1P1Z4/5U+/+lODcNOyhRXUuptoCjgiaVRz9MLhS+fwzykjebv9zpx22CV81sjITc8dts5QhSJSCKrmVTcIN+e//FAk3Pw5Rrhpu1VLhv96r4JadxNNAUckTWrCzU/fncs/p4zinfY7c1r/kXzWuk2jr92hzRYZqFBECsWwKQvq3D//5YcY/Ow9kXCzMSrc3Nx/v4INNdF0FJVImtSEm3GTR/JO+64MTDLc1LxWRCQZA+94iS+/3Ri5/9vZQbh5+AcNw01pSeuiCDeggCOSNk0NNyIiyap/Er/fzn6IIbOCcPOn4+qGm0I7z01jNEUlkiLRJ9RSuBGRTIgVbqb+4CcNwk1pnl9XqikUcERSIDrc/CQF4aZPj3apLrGgmNndwHHAanffO9x2PfBr4FtgKXC2u68LHxsCnANsBP7g7jOyUbdIKl1eVbvuZtDsSZFw88fj/hwJNz132Jon/3R4lirMLk1RiaRAdLi5Y/JIlmy/U1LhZtstWjQIM316tOP+8w5JW60F4h7gmHrbngT2dvd9gbeBIQBmtidQBuwVvuY2M2uBSJ6qmldNn7HPRI6aGjR7EkNn/atBuDmtd9eiDTegERyRlKkTbsoSh5tYZymW5Ln7c2bWrd62J6LuvgycHN7uC1S6+zfAe2a2BDgIeCkTtYqkStW8ai5+cH6dbefNnszQWf/ikT1+XCfcgM58rhEckRT48Xuv1Qk361pvm+2Sit1vgMfC26XAiqjHVobbRPJGvHAzbNbdPLLHj7n413+pE25O6901wxXmHo3giDTBvsMf5/NvgsMyf/zea9w56RqWbt8lqXCz7RaaHUknMxsGfAfc34TXDgIGAXTtqv9BSO64bNIbde6f+0oQbh6NEW5AozegERyRTVY/3NwxeSRLt+/CgLJRSYWbN66qv3REUsXMziJYfDzQ3T3cXA3sFPW0LuG2Btx9nLv3cvdeHTp0SGutIsnaY9h0vvnu+8j9c1+ZzOUzg3BzUYxwI4G4IzhmlvAwDnfXmcikKNWEm8Pem8cdk0fybrtSBvaPP3Kj9TYNpaO/mNkxwKXAT939q6iHpgIPmNmNQGegJ/DKpr6/SDZ0Gzytzv1zXpnSaLjRUZiBRFNUcwEHLMZjDuySlopE8sBh783jzsnX8F7bzgzsP5JPt9ou2yXlm2b1FzOrAA4H2pvZSmA4wVFTWwBPmhnAy+5+vrsvMrMJwJsEU1cXuPvG2O8skjtihZsrZt7Fo7sfljDc6CjMQNyA4+7dM1mISL7os2x+JNwMKBulcNMEze0v7l4eY/NdCZ4/ChjVnM8UyaRdh9QLN69W1Yab4y8pymtLbaqkFhmbWVuCYd0ta7a5+3PpKkokl0TvRfVZNp+7Jl2tcJNC6i8ide0xbDrfee39c16t4opn7uTR3Q9rsKC4mK4ttakaDThmdi5wEcHCvPlAb4LzRxyZ1spEckBzw43W3ySm/iJSV4NpqTDcTNu9Dxf/+i9816L2f9vFdm2pTZXMCM5FwIEE89lHmNkewOj0liWSW5INNwo0m0z9RQQ4eNSTfPTFt3W2RYebi359SZ1wU4zXltpUyQScr939azPDzLZw9/+amSKjFI1DNS2VTuovUvRihZvfvPowVzxzJ9N3O7RBuNFC4uQkE3BWmlkJUEVwdMKnwPJ0FiWSKw5dNp+7w3AzUOEmHdRfpOjFCjdXPnMH03c7lD8cf6nCTRM1GnDc/YTw5ggzmwlsR+0p0EUKTs0ceDByc00k3HyicJNy6i9S7LrXW3Nz9pz44UZT4Jum0TMZm9m/a267+7PuPhW4O61ViWRJTbg5ZPnr3DXpGpaX7Jh0uFHz2XTqL1LMLq9aQNTBUpw952GGPx073OjaUpsumSmqvaLvmFkL4EfpKUck+w5Z/jp3P3Q1y0t2ZED56JjhprSkNS8M1oE+KaD+IkXrvpffj9w+a85Uhj99B4/FCDega0s1RaJLNQwBhgKtzexzas84+i0wLgO1iWRcdLhJNHLzwbr1Ga6ssKi/SLGKdVXws+ZMZcTT43hst0O5sF642dxgyRiNDjdFojMZjwHGmNkYdx+yqW9sZlsCzxGcOn1z4CF3H25m3YFKYHuC07Wf7u7fxn8nkdS4vGoBFbNXsNGdFmaUH7xTnT2oQ5a/wd0PXc37JR0ZWDaKj7cuiftenUtaZ6DiwtXc/iKSb6Iv0hvtzLmPMOLpcTy+2yENwo0WFDdPMouMh5jZ8cBPwk2z3P3RJN77G+BId/+fmbUE/mNmjwF/Am5y90oz+wdwDnB7E+sXScrlVQvqhJmN7jHCzVW8X9KRAWWjE4YbnVwrdZrRX0TyRvfB0+qstalx5txHuOqpf/L4bofw++MvazAtpXDTPMksMh5DcDKuN8Ofi8ys0RNxeeB/4d2W4Y8TnKH0oXD7eKDfppctsmkqZq+I+1jv95MLN0aw9mbMifvo5Fop0tT+IpIvBt7xUsxwc0ZUuKk/cgM6aCEVkllk/CtgP3f/HsDMxgPzCObPEwoXDM4FdgVuBZYC69z9u/ApK4GY/6cws0HAIICuXbV6XJpno8dqMUG4+dfE5EZu3lPDSYcm9xeRfPDC0k8abDtj7iNc/dQ/mdGzNxcefykbWrSs83gLswavkU2X1MU2gRKg5r9S0icDcfeNwH7hibymAHtswmvHES427NWrV+z/O4kkqYVZg5BTE25WbNd4uFHDSasSmtBfRHJZ/WtK1YgON7/ve1mDcANQfvBO6S6vKCQ6iuoJd/8FMAaYF56Eywjmygdvyoe4+7rw9YcAJWa2eTiK0wWobnL1InHEay41aqalVmzXkfLyxOEG1HBSLZX9RSTXxOs/p7/2aMJwU3Pwgw4JT41EIzgdANy9wsxmEVwQD+Ayd/+wsTc2sw7AhjDctAZ+DlwLzAROJjiS6kzg4aaXL9JQY+Hm4PcXcPdDV7Fy244MKG94tNRpvbs2ONpKDSflmtVfRHLVwDteirn99Nce5Zon/8ETMcKN1tukR6KAs52ZnRhj+6FmhrtPbuS9OwHjw3U4mwET3P1RM3sTqDSzkQRz7Xc1qXKRJjj4/QX866ERkXCzduu2MZuLAk3aNbe/iOSkWGtuosPNBVHhRsEmvRIGHOA4ak/AFc2BhA3I3d8A9o+x/V3goE2oUSQlasJN9bY7RMKNZE2z+otIrol1RXCA016bFjPcaFVf+iUKOMvd/TcZq0QkjQ5asTASbsrLRyvcZF+z+ouZ3U0QkFa7+97htnbAg0A3YBlwqrt/amYG/BX4JfAVcJa7v9a88kVq7TFsOl9vbHgszGmvTWPkk7fz5K4H1wk3oKMyMyFRwFHAlLyRaN3NQSsWcs/E4cHITZnCTY5obn+5B7gFuDdq22DgaXcfa2aDw/uXAccCPcOfgwlOLHpwMz9fBAhOIhoz3MybHgk3/9dvcJ2RG4WbzEgUcE7PWBUizdBYuPnXxBF80KYDA8pGs2abuuFGc+BZ06z+4u7PmVm3epv7AoeHt8cDswgCTl/gXnd34GUzKzGzTu6+qjk1iEDsk4ieNm86I5+4jSd3PahOuFG/yaxE16JamMlCRFKtJtysatOe8vIxrNkm9oJiybw09ZeOUaHlQ6BjeLsUiP6/UM0JRhsEHJ1gVDZV/fNr1Q03Q2Ke50YyI9kT/YnklQNjhBspHu7uZrbJJwjVCUYlWbFGjgdGhZsL+g7RoeBZlsy1qC5KZptIrjhwxULuUbjJCynuLx+ZWafwPToBq8Pt1UD0mRp1glFplnjhZtQTt/FUjwO5oO8Qvt1c4SbbkhnBOZPgCIRoZ8XYJpIRidbc1A03DdfcSM5JZX+ZGr7fWOqeRHQq8HszqyRYXPyZ1t9IU1XNa5iNB8x/LBJu/q/f0Ei4UbDJrkSXaigHBgDdzWxq1EPbUnvdGJGMShRueq1cxD0TR/BhJNy0q/O4mk3uaG5/MbMKggXF7c1sJTCcINhMMLNzgOXAqeHTpxMcIr6E4DDxs1P0x5AiUzWvmosfnF9n24D5jzF6xq0Nwo1kX6IRnBcJFuG1B26I2v4F8EY6ixLZVL1WLmL8hOF82KY9ZVHhRqEmZzWrv7h7eZyHjorxXAcuaEKNIhGxdq7K5z/O6Bm38nSMcHNz//0yWJ3EkugoquXAcjP7GbDe3b83s90Irgi+IFMFijQmeuSmLMbIjeQe9RfJJ/HCzZgZt/B0jwP5Xb1w03OHrem3f2kmS5QYklmD8xzwYzNrCzwBvAr0BwamszARaPzCmT9a+Sb3TBzBR9u0U7jJT+ovkncShZs+Pdpx/3mHZLE6qdHoUVSAuftXwInAbe5+CrBXessSSS7cjJ84PAw3YxRu8pP6i+S0+n2oLAw3z+zSq0G4ubn/fgo3OSSZERwzs0MI9qjOCbe1SF9JIo2LDjflZbFHbrT+Ji+ov0hOirWDVTb/ccaG4eb8E4Y1CDealsotyQSci4AhwBR3X2RmuwAz01uWSHwHrHyrTrhZ3WZ7QIEmT6m/SM6JFW76vz4jEm5+d0LdkZvTendVuMlBjQYcd3+OYJ685v67wB/SWZRIPAesfIt7J17ZINxIflJ/kXzQ//UZXPv435m5y4/43QlD+WbzVpHHTuvdlZH99slidRJPowHHzDoAlxLMi29Zs93dj0xjXVKEGltzUxNuVm/dVuGmQKi/SK6p34eiw835JwyrE246tmmlcJPDkpmiuh94EDgOOJ/gDKFr0lmUFJ9kws34MNyUlY9pEG40PZW31F8kJ8TqQae+/kTccAMwe9jPM1WeNEEyAWd7d7/LzC5y92eBZ83s1XQXJlLjgOog3KzZui3l5cHIjQJNwVB/kayLF26ue/xvzOoeO9yoB+W+ZALOhvD3KjP7FfABoONxJSMOqH6L8RNqw81HbdpnuyRJLfUXyTmnvPEEYx//O7O6/4jfnqhwk6+SCTgjzWw74M/A3wmuFfPHtFYlQm24Wbt1icJN4VJ/kayqP3pzyhtPcO1jf+e57gco3OS5ZI6iejS8+RlwRHrLkWKSaN3N/tX/jYSbsvIxCjcFSv1FsiVW/6kJN89331/hpgA0eiZjMxtvZiVR99ua2d1prUoKXmPh5t4JVwQjN2UNw42aTOFQf5FsiB1unoyEm0EnXq5wUwCSmaLa193X1dxx90/NbP/0lSTFrGbk5uOtgnDz4bbt1VgKm/qLZF0Qbv6mcFNgkrkW1WbhhfAAMLN2JBeMRDZJTbj5ZKvtKCsPwo0UPPUXyaiGa26CcPOfbvsp3BSYZBrJDcBLZjYxvH8KMCp9JUkx2u+DxQo3xUn9RTIi1rTUyQueioSb8xRuCk4yi4zvNbM5QM2ZRU909zfTW5YUmkRrbvb7YDH3PniFwk0RUn+RTIgXbq6b/tfacNNyizqPK9zkv2SmqCA4L8WX7n4LsMbMuqexJikwyYab8vLRDcKNmkxRUH+RjDppwdMKN0UgmWtRDQd6AbsD/wJaAvcBfdJbmhS6H4bh5tOttqW8fDSrtu0AqLkUk3T0FzP7I3Au4MAC4GygE1AJbA/MBU5392+bVbzkhfo7WCcteJrrp9+scFMEkhnBOQE4HvgSwN0/ANqksygpfD/8YDH/DsNNWfmYSLiRopPS/mJmpQRXI+/l7nsDLYAy4FrgJnffFfgUOKeZdUuO6zZ4msJNkUtmkfG37u5m5gBmtnWaa5I819iFM3/4wWL+PeFKhRuB9PSXzYHWZrYB2ApYRbDGZ0D4+HhgBHB7Cj5LclCsHnTiwiDcvLDzDxVuikQyIzgTzOyfQImZnQc8BdyR3rIkXzUWbvZd9XYQblq3qTMtJUUrpf3F3auB/we8TxBsPiOYklrn7t+FT1sJlMZ6vZkNMrM5ZjZnzRpd1LxQnLjwaf7ftCDcnHvSFQo3RSLhCI6ZGfAgsAfwOcE8+ZXu/mQGapMCs++qt7nvwStYt+U2lJeP5oNtd2jwHDWa4pGO/hKeU6cv0B1YB0wEjkn29e4+DhgH0KtXL29qHZI99XeyTlj4TCTcnHeSRm6KScKAEw4dT3f3fQCFGmmy6HBTNmBMJNyouRSvNPWXnwHvufsaADObTLBgucTMNg9HcboA1Sn6PMkRVfOqufjB+XW2nbDwGW6YdhMv7rwv5510OV+33LLO4+o/hS2ZKarXzOzAtFciBWufVe/EDDcipL6/vA/0NrOtwhGio4A3gZnAyeFzzgQeTuFnSpZVzatmyOQFdbb1WzQzEm7OPekKhZsilMwi44OBgWa2nOBIByPY+do3rZVJXmhszc0+q97h/gcvD6elFG6kgZT2F3efbWYPAa8B3wHzCKacpgGVZjYy3HZXKoqX3FB/5Kbfopnc+OiNCjdFLpmAc3Taq5C8tKnhpnq7uuFGTUZIQ39x9+HA8Hqb3wUOSvVnSfbV70M14ealnfdRuClyyVyqYbmZ/RD4cbjpeXd/Pb1lSb4LpqUu57OocKPGIvWpv0hTxdrBqpmWemnnfTjnpCsVbopco2twzOwi4H5gh/DnPjO7MN2FSf7a+8Ml3Pfg5Xy+5TaUxRi5Eamh/iJNESvc9A3Dzctd91a4ESC5KapzgIPd/UsAM7sWeAn4ezoLk/y094dLuL9ymMKNJEv9RZqt76KZ3DjtJmbvpHAjtZIJOAZsjLq/MdyW+EVmOwH3Ah0Jrgkzzt3/ambtCM590Q1YBpzq7p9uWtmSLYnW3USHm/Ly0Qo3kowm9RcpXvV70PFvzoqEm9+crHAjtZIJOP8CZpvZlPB+P5I7AuE74M/u/pqZtQHmmtmTwFnA0+4+1swGA4OByza5csm4ROFmr3Ba6osttqa8fDQrt+tY53E1GYmjqf1Fikys/nP8m7O46dEbNXIjMSWzyPhGM5sFHBZuOtvd5yXxulUEp0rH3b8ws7cITo/eFzg8fNp4YBYKOHltrw+XcP+Dl/O/VltRNmAMK7frqMYiSWlqf5HikijcvLLTXpxz0pWsb6VwI3XFDThmdo+7nxXe3cfd/9bUDzGzbsD+wGygYxh+AD4kmMKK9ZpBwCCArl27NvWjJc1ihRuRxqSyv0jxOf7NZyPh5jcnDVe4kZgSHUX1w6jbFzX1A8xsG2AScLG7fx79mLs7wfqcBtx9nLv3cvdeHTrogoy5SOFGmiEl/UUKX8M1N89y06M3KNxIoxJNUTX7QnNm1pIg3Nzv7pPDzR+ZWSd3X2VmnYDVzf0cSY+Ea24+WqpwI82hC1lKQrH6z68VbmQTJAo4XczsbwRHNNTcjnD3PyR64/A6MHcBb7n7jVEPTSW4FsxYdE2YnNVouKkcxv9atY4ZbtRkJAnN6i9S2OKFm5sfvYFXu+ypcCNJSRRwLom6PacJ790HOB1YYGbzw21DCYLNBDM7B1gOnNqE95YsqRNuymvDjZqLbKLm9hcpItHh5uyTRyjcSFLiBhx3H9+cN3b3/xD/fBZHNee9JTv2+mgp91VeXhtuSnbMdkmSp5rbX6Rw1R+9Oe6t57j50RuYo3Ajm6jRSzWIAOz50bvcV3k5X7baknKFGxFJg1jh5q+P/D/mdNmTsxRuZBMlc6I/KQKJ1tzs+dG73F85jK9aBuFmhcKNiKRYzJGbMNycfbLW3Mim0wiOJB1uygbEDjdqNCLSHPV70K/eep6bH/l/zC39AWefPJyvWrWu87h6jiSj0REcM9uS4IJ4ewGRCO3uv0ljXZIDfrA6drhRc5FUUX+RWOHmr49cH4SbU0Yo3EiTJTOC829gR+Bo4FmgC/BFOouS7PvB6nd5oKJmWmq0pqUkXdRfilj9cPPL//5H4UZSJpk1OLu6+ylm1tfdx5vZA8Dz6S5M0ifRlBTUjNxczvqWW1BePpr323bKUGVShNRfilSscPO3qdcp3EjKJDOCsyH8vc7M9ga2A3ZIX0mSTsmGm683b0VZ+RiFG0m3lPcXMysxs4fM7L9m9paZHWJm7czsSTN7J/zdttmVS5PFCzevle4RM9yINEUyAWdc2AyuIDgL8ZvAdWmtSrIimXCjvShJsXT0l78Cj7v7HgTXvHoLGAw87e49gafD+5IFicLNWadcFTPcqO9IUzQ6ReXud4Y3nwV2SW85ki17rH4vZrhRY5F0SnV/MbPtgJ8AZ4Xv/y3wrZn1BQ4PnzYemAVc1tzPk01TP9wcGz1yc7KmpSS14gYcM/tTohfWu76U5LE9Vr/HA5XD+HrzVlpzIxmRxv7SHVgD/MvMfgjMJbhaeUd3XxU+50Mg5tVhzWwQMAiga9euTSxBYokVbv4+9TrmdQ7CzZdbbFXncYUbaa5EIzhtMlaFpFWidTfByM0wvmnRkvLy0Sxv2zmDlUkRS1d/2Rw4ALjQ3Web2V+pNx3l7m5mMa9m7u7jgHEAvXr10hXPU6R+Dzpm8QuRcHPWKQ3DjUgqJLoW1VWZLETSI1G42X3NMu6vHMa3LVpSNmBMg3CjPShJlzT2l5XASnefHd5/iCDgfGRmndx9lZl1Alan6fOlnljh5paHr00YbtR7JBWSOdHfv4AGezI6EVd+233NMh6oGFon3KipSKalur+4+4dmtsLMdnf3xQQX9n0z/DkTGBv+frjpVUuy4o3czO+8e8xwox4kqZTMeXAejbq9JXAC8EF6ypFMiA43mpaSLEtHf7kQuN/MWgHvAmcTHDE6wczOAZYDpzbzM6QRuw6pG26OXvwif596Ha932o0zT7lK4UbSLpmjqCZF3zezCuA/aatI0mq3MNxsaLE55eWjWdauNNslSRFLR39x9/lArxgPHdWc95Xk1R+5OXrxi9wy9Vpe77QbZ8UINyLp0JSrifdEJ/rLWYnW3Oy2ZhkVYbgpKx+jcCO5SP0lzzUIN2/XDTf/05obyZBk1uB8QTBHbuHvD9H5I3JSY+HmgcphccONGoxkg/pLYYkZbh6+ljd27Bkz3KjvSDolM0Wlw8XzXE24+W6zFpRHhRs1F8k29ZfCkSjcnHnq1Qo3knFJTVGZ2YnAYQR7WM+7e1U6i5LUqR9u3tO0lOQY9Zf8t6nhRiQTGr0WlZndBpwPLAAWAueb2a3pLkyar+ea5TxQOYyNCjeSo9Rf8l/9cPOLt1/iloevZcGOu8YNNxq9kUxIZgTnSOAH7u4AZjYeWJTWqiQpidbc9FyznIrKoWzcrAVlCjeSu9Rf8liscHPrw2NZsOOunHHqNZqWkqxK5mriS4Doi7LsFG6TLGos3NSM3MQLN2o0kiPUX/JU/R7083dejoQbrbmRXJDMCE4b4C0zeyW8fyAwx8ymArj78ekqTjbdrmvf54HKYXy/2WaUl42OhBs1F8lR6i95KFa4ua1qDAs7BuHmiy22zlJlIrWSCThXpr0KSYld175PRcXQSLh5d/su2S5JpDHqL3km5shN1VgWdtyVM/rHDjfawZJsSOYw8WfNrCPBnhXAK+6uC9VlWKIpKagNN26mcCN5Q/0lv9TvQz97Zza3Vo1lUcceMcONgo1kUzJHUZ0KvAKcQnD9ltlmdnK6C5NaSYWbyiDclJWPUbiRvKH+kj9ihZvbqsYo3EjOSmaKahhwYM1elZl1AJ4CHkpnYZKcSLjBKC+PPXKjRiM5TP0lD8QLN2923CXutJRItiUTcDarN2T8MckdfSVp1mPtijrhZun2OwEKNJJX1F9yXP1wc9SS2nBzev9rtOZGclYyAedxM5sBVIT3+wOPpa8kSUaPtSuorBzSINyI5Bn1lxwWK9zcPiUcuYlxtJSCjeSSZBYZXxJ1KnWAce4+Jb1lSaJ1NzXhBlC4kbym/pK74oeb7pxx6tV8vuU2dR5XuJFcEzfgmNmuQEd3f8HdJwOTw+2HmVkPd1+aqSKLTcJw83FtuCkrH9Mg3KjJSD5Qf8lt9XvQkUteiQo31zQINyK5KNEIzs3AkBjbPwsf+3Ua6pEEeny8gsqKMNyUBeFGgUby1M2ov+SkWOHmH1NGJww36kOSixIFnI7uvqD+RndfYGbd0leSxNLj4xVUVAwFwnDTXtNSktfS2l/MrAUwB6h29+PMrDtQCWwPzAVOd/dvm/s5he7IJa9we9Vo3tohdrhRsJFcluhohZIEj7VOcR2SQE24MZzystEKN1IIShI8lor+chHwVtT9a4Gb3H1X4FPgnBR8RkE7Yumr3F41mv926M7p/RVuJP8kGsGZY2bnufsd0RvN7FyCPSBJkURrbnb5eGUQbtwpLx/NkvZd4z5XJI+krb+YWRfgV8Ao4E9mZgRXLR8QPmU8MAK4vTmfU0h2HTKN77z2/hFLX+UfU0bFDTci+SBRwLkYmGJmA6ltOL2AVsAJaa6raDQWbiorhsQNN9qDkjx2MenrLzcDlxJcyBOCaal17v5deH8lUNrMzygY9cPN4UmEG/UeyQdxA467fwQcamZHAHuHm6e5+zMZqazI7fLxSioqG47cqLFIIUhXfzGz44DV7j7XzA5vwusHAYMAunYtjtHS+uHmn1NGsbhDN01LSd5L5jw4M4GZGahFQjXhZrPvv2dA2ShNS0nBSkN/6QMcb2a/BLYEtgX+CpSY2ebhKE4XoDpOPeOAcQC9evXyWM8pVNHh5rT+IzUtJXkvbadEN7O7zWy1mS2M2tbOzJ40s3fC323T9fn5qn64eafDztkuSSRvuPsQd+/i7t2AMuAZdx9IEKJqLuJ5JvBwlkrMGd0GT4tMkScTbjq2aZXpEkWaJZlLNTTVPcAtwL1R2wYDT7v7WDMbHN6/LI015JxEa266f1JNReVQWny/kfKy0Qo3IqlzGVBpZiOBecBdWa4nq6L70OFL5/DPKaN4u/3OCcPN7GE/z2SJIs2WtoDj7s/FOJ9FX+Dw8PZ4YBZFFHAaCzeVFUMShhvNf4skz91nEfQY3P1d4KBs1pMrGoabkbzdfmcGlo3SmhspKOkcwYmlo7uvCm9/CHTM8OfnpO6fVFMRI9youYhIuvz03bmNjtyI5LNMB5wId3czi7uIr1iOZugWhpvNv9+oNTcikjbRIzc/fXcu4yaP5O32XTmt/0g+a90mwStF8lOmA85HZtbJ3VeZWSdgdbwnFsLRDImmpCAIN5VR4ebtDt0yU5iIFJVY4eadRsKNRpAl36XtKKo4phIcwQAFfiRDsuGm5cbvFG5EJG2ie9FPosLNQIUbKXBpG8ExswqCBcXtzWwlMBwYC0wws3OA5cCp6fr8XBZMSw2l5cbvKC8fHTPcqMGISCr95N253NFIuFHfkUKSzqOoyuM8dFS6PjMf7PzpB1RUDKXVxg0MiAo3aiwikkr1R24aCzcihSbTU1RFbedPP6DygSGRcLNY01IikgbR4ebH773GHZNHsmT7nTQtJUUla0dRFaJE624UbkQkE+qHmzsnXROEmzKFGykuCjgpkmy4GVg2qkG4UXMRkVSrH27Wtd62wXPUe6SQKeCkWc2amy02bmBA2Sj+u0N3NRURSblY01JLt+8SN9yIFDqtwUmjrp+uoqJiKFt+920k3IiIpFp0uDnsvXncMXkk77YrZUDZqLjhRjtaUugUcNKk66erqKwYonAjImlVP9zcOfkahRsRNEXVZInW3ESHm4FlIxVuRCTt+iybHwk3A/trzY2IAk4TNBZuaqalBpaN5K0ddqnzuBqMiKRCdB/qs2w+d026mvfadmZg/5F8utV2WaxMJDco4KTQTus+pKJiKFtt+JoB5aMi4UahRkRSKV64GVA2Km64UR+SYqM1OCmy07oPqXxgSINwIyKSStHh5lCFG5G4NIKTAtHhZmCZwo2IpN7AO17ihaWfRO4fumw+dzcSbhRspJgp4CQh0Zqb+uHmzY4KNyKSWnsMm87XGz1yP5lwI1LsNEXViMbCTUVF4nCjPSiRzDGzncxsppm9aWaLzOyicHs7M3vSzN4Jf7fNdq3J6jZ4WoNwc9eka4IFxZqWEolLIzhN1CUMN1t/+zWnlY2MhBs1FZGs+g74s7u/ZmZtgLlm9iRwFvC0u481s8HAYOCyLNaZlH2HP17n/iHLX+euSdewrG0nBpaN4hOFG5G4FHCaoMu6D6msGMI2365nYNkoFnXske2SRARw91XAqvD2F2b2FlAK9AUOD582HphFHgScz7/ZGLl9yPLXufuhqxOGGwUbkVoKOPUkmpIChRuRfGFm3YD9gdlAxzD8AHwIdIzzmkHAIICuXbtmoMr4ukf1oppws7xkx7jhpmObVpksTyTnaQ1OlEbDzWcfUVkxVOFGJMeZ2TbAJOBid/88+jF3d8Bjvc7dx7l7L3fv1aFDhwxUGlu3wdMiBR6y/I1IuBlQPjpuuJk97OeZLVIkx2kEJ0ldPvuIygeGsM23XzGw/8iY4UbDwyLZZ2YtCcLN/e4+Odz8kZl1cvdVZtYJWJ29ChP7+Y2zIreDcHNV3JEb9RyR+BRwktAg3Oy4K6DmIpJrzMyAu4C33P3GqIemAmcCY8PfD2ehvIQur1pAxewVbPRg7Kb3+0G4eb+kIwPLRvHx1iV1nq/+I5KYAk4j4oUbEclJfYDTgQVmNj/cNpQg2Ewws3OA5cCp2SkvtsurFnDfy+9H7vd+/w3+NTEINwPKRjcINyLSuKIPOInW3XT57CMqKobS5psvGVA2SuFGJMe5+38Ai/PwUZmsZVPUDzc1Izfxwk2fHu0yWJ1IfirqRcbJhJttv/5fsKC4XrjR8LCIpELVvOrI7Zpws3LbxOHm/vMOyWCFIvmp6EdwYin9bHWdcLNwx10VaEQkLa6fsRiAg99fEAk35eUNw40B76kPiSStqEdwYin9bDWVFUPqhBsRkXSpXreeg99fwL8eGhGM3JQ3XFAMMLB3ds/LI5JvFHCiRIeb0/qPVLgRkbSrH27Wbt3wMlmn9e7KyH77ZKE6kfxVVFNUidbcBNNSteFmQaeeGaxMRIrR83dOShhuNDUu0nRFM4KTKNx0/jwINyVxwo2ajIik3HPP0ev/TqN62x3ijtyISNMV1QhOLJ0/X03lA0G4GRgVbhRqRCSVquZVc/2MxXywbj3HfvoOf/v3MFZu24Hy8tExw03rlkWz/ymSFkUdcKLDzWn9r9G0lIikRdW8aoZMXsD6DRs5cMVCrp84gmXbdeB3Z13H2hZtYr5mzIn7ZrhKkcJStAEnmJYaGgk3b3TaLdsliUiBun7G4ki4uWfiCFa1aU95/9Fs2KYdrb/7nvUbNkaeawRHTPXbvzR7BYsUgIINOD+/cRbvrP4y5mM14abt+i84/dSrFW5EJOWip6Qc6oab8jGs2aYttn4DN/XfL/K8ziWtueTo3RVuRFKgIANOonDT6fM1dcLN6513b/Acrb8RkeaInpKC2OEGoHNJa/rtX6pAI5IGBRlwEoWbyoohtP3qc07vf00k3CjQiEgq1UxJAfRauSgq3IyOhJvWLVtwydENd7BEJDUKMuDEEi/ciIik2gfr1gNBuBk/YXhUuGmHgaaiRDKgYAJOUtNSCjcikgGdS1rTaeEc7pk4gg+jwk1pSWteGHxktssTKQoFEXAShZsdP19LRcVQ2n31GWfECDeWiQJFpKiM6fg5Pxo+gg+32Z6yMNxoSkokswoi4CQKN5UVQyLhZn6McKOr84pISr3wAj/5w+l80bkzfywfw9rvt6JUU1IiGVcQASeW6HBz5qlXM7/z7go0IpJeL7wAxxwDnTvTZtYspnbqlO2KRIpWQQacHT9fS0XlELb/ah1nnHoN80r3AIJ5cRGRtHjxxUi4YdYsULgRyaqsXOzEzI4xs8VmtsTMBjf3/XrusHXkdscvgnDT/su64Ubz3yLFLdV9p44XX4Sjjw7CzcyZCjciOSDjAcfMWgC3AscCewLlZrZnc97zgiOCa0h1/CKYlqoJN8t67oMBpSWtGXPiPpr/FilS6eg7EdEjNzNnBr9FJOuyMUV1ELDE3d8FMLNKoC/wZlPf8PoZi+uEmzNPvZp5pXtQ2mpz5l35ixSVLSJ5LOV9B6gNNzvuqHAjkmOyMUVVCqyIur8y3FaHmQ0yszlmNmfNmjUJ3/CDdevp+L9P2OK7DZx56tW8VvqDyHYREdLQdwBYsgRKS4M1Nwo3IjklK2twkuHu49y9l7v36tChQ8Lndi5pzRudduPwQXdEwk3NdhGRZG1K3wHgjDNg/nyFG5EclI2AUw3sFHW/S7ityS45endat2zBt5u3jGzTomIRiZLyvhOxxRYpeRsRSa1sBJxXgZ5m1t3MWgFlwNTmvGG//UsZc+I+lJa01qJiEYkl5X1HRHJbxhcZu/t3ZvZ7YAbQArjb3Rc193377V+qQCMiMaWr74hI7srKif7cfTowPRufLSLFSX1HpLjk7CJjERERkaZSwBEREZGCo4AjIiIiBUcBR0RERAqOAo6IiIgUHHP3bNfQKDNbAyxP8untgbVpLCedVHv25HP9uVT7Wnc/JttFpMIm9J1c+v6zTd9FLX0XtdL9XcTsO3kRcDaFmc1x917ZrqMpVHv25HP9+Vx7IdD3X0vfRS19F7Wy9V1oikpEREQKjgKOiIiIFJxCDDjjsl1AM6j27Mnn+vO59kKg77+Wvota+i5qZeW7KLg1OCIiIiKFOIIjIiIiRU4BR0RERApOwQQcMzvGzBab2RIzG5ztehpjZneb2WozWxi1rZ2ZPWlm74S/22azxnjMbCczm2lmb5rZIjO7KNye8/Wb2ZZm9oqZvR7WflW4vbuZzQ7//jxoZq2yXWsiZtbCzOaZ2aPh/byqv1DkW99JNTNbZmYLzGy+mc0Jt+V8H2iuTenfFvhb+HfkDTM7IHuVp16c72KEmVWHfy/mm9kvox4bEn4Xi83s6HTWVhABx8xaALcCxwJ7AuVmtmd2q2rUPUD9ExMNBp52957A0+H9XPQd8Gd33xPoDVwQft/5UP83wJHu/kNgP+AYM+sNXAvc5O67Ap8C52SvxKRcBLwVdT/f6s97edp30uEId98v6jwn+dAHmuseku/fxwI9w59BwO0ZqjFT7qHhdwFBP9ov/JkOEP77KAP2Cl9zW/jvKC0KIuAABwFL3P1dd/8WqAT6ZrmmhNz9OeCTepv7AuPD2+OBfpmsKVnuvsrdXwtvf0HwP9pS8qB+D/wvvNsy/HHgSOChcHtO1l7DzLoAvwLuDO8beVR/Acm7vpMhOd8HmmsT+3df4N6w97wMlJhZp4wUmgFxvot4+gKV7v6Nu78HLCH4d5QWhRJwSoEVUfdXhtvyTUd3XxXe/hDomM1ikmFm3YD9gdnkSf3h9M58YDXwJLAUWOfu34VPyfW/PzcDlwLfh/e3J7/qLxSF0neaw4EnzGyumQ0Kt+VFH0iDeH/uYv178vtwSu7uqGnKjH4XhRJwCo4Hx+/n9DH8ZrYNMAm42N0/j34sl+t3943uvh/QhWDvYY/sVpQ8MzsOWO3uc7NdiwhwmLsfQDANc4GZ/ST6wVzuA+lUrH/uKLcDPQiWAawCbshGEYUScKqBnaLudwm35ZuPaoYuw9+rs1xPXGbWkiDc3O/uk8PNeVM/gLuvA2YChxAMG28ePpTLf3/6AMeb2TKCKZEjgb+SP/UXkkLpO03m7tXh79XAFIIdhrzqAykU789ddH9P3P2jcEfye+AOaqehMvpdFErAeRXoGR5J0opgEdPULNfUFFOBM8PbZwIPZ7GWuMI1H3cBb7n7jVEP5Xz9ZtbBzErC262BnxOsIZoJnBw+LSdrB3D3Ie7exd27Efw9f8bdB5In9ReYQuk7TWJmW5tZm5rbwC+AheRBH0iTeH/uqcAZ4dFUvYHPoqayClK9NUYnEPy9gOC7KDOzLcysO8HC61fSVoi7F8QP8EvgbYL1FMOyXU8S9VYQDN1tIJiHPIdgLcXTwDvAU0C7bNcZp/bDCIZf3wDmhz+/zIf6gX2BeWHtC4Erw+27hP/QlgATgS2yXWsSf5bDgUfztf5C+Mm3vpPiP/suwOvhz6KaP38+9IEU/NmT7t+AERxttxRYAPTKdv0Z+C7+Hf5Z3yAINZ2inj8s/C4WA8emszZdqkFEREQKTqFMUYmIiIhEKOCIiIhIwVHAERERkYKjgCMiIiIFRwFHRERECo4CTgEws43hFVsXmtlEM9sqwXMPN7NDo+7fY2Ynx3t+jM+o+emWovKbxcwujvXnNbMpYZ1LzOyzqLoPjfHcfslcJDG8Qu5fUlW7SL5Sz1HPyQcKOIVhvQdXbN0b+BY4P8FzDwca/IPbhM+o+VmWzIuizq6bLhcDDZqNu5/gweUYzgWej6r7xRjv0Y/gatAikhz1nHrUc3KPAk7heR7Y1cx+bWazzWyemT1lZh3DPaDzgT+GexY/Dl/zEzN70czeTWbPqoaZ7WdmL4cXVJtSc0E1M5tlZjeb2RzgIjP7kZk9G16Qb0bU6cx3DWt73cxeM7MeZraNmT0d3l9gZn3D525tZtPC5y40s/5m9gegMzDTzGYmUW83M3smrPdpM+sa7l0dD1wffic9zOw8M3s1/KxJifZORUQ9J0G96jnZlO2zIOqn+T/A/8LfmxOcHvx3QFuInMjxXOCG8PYI4C9Rr72H4My3mxHsUSyJ8xkbqT1r8ZRw2xvAT8PbVwM3h7dnAbeFt1sCLwIdwvv9gbvD27OBE8LbWxLsFW0ObBtua09wZl4DTgLuiKpnu/D3MqB9gu/mcGrP9vsIcGZ4+zdAVdR3cHLUa7aPuj0SuDDWd6cf/RTrj3qOek4+/KR7KE8yo7WZzQ9vP09wnajdgQfDPZdWwHsJXl/lwUXR3jSzjnGes96D4VcAzGw7oMTdnw03jSdoWjUeDH/vDuwNPGlmAC2AVRZcw6bU3acAuPvX4fu2BEZbcFXi74FSoCPBab9vMLNrCZrH8wn+PPEcApwY3v43cF2c5+1tZiOBEmAbYEYTPkukkKnnJEc9J4sUcApDnUYAYGZ/B25096lmdjjBnkA830S/NEU1fRn1fovc/ZB69bWJ87qBQAfgR+6+wYKrZm/p7m+b2QEE1/4ZaWZPu/vVKaq1vnuAfu7+upmdRbBHJiK11HNS6x7Uc1JOa3AK13bUXob+zKjtXwDx/qEnzd0/Az6NmlM/HXg2xlMXAx3M7BAI9pbMbC93/wJYaWb9wu1bhPPO2wGrw0ZzBLBz+Hhn4Ct3vw+4HjigCX+eFwmu+AxBU6vZI6v/Hm0I9vhahs8Tkcap5zSknpNFCjiFawQw0czmAmujtj8CnFBvwV9TnUmwUO4NYD+COfE63P1b4GTgWjN7nWA+veaIitOBP4SvfxHYEbgf6GVmC4AzgP+Gz90HeCUcFh9OME8NMA54PJkFf8CFwNnh550OXBRurwQuCRdH9gCuIJirfyHq80UksRGo59SnnpNFupq4iIiIFByN4IiIiEjBUcARERGRgqOAIyIiIgVHAUdEREQKjgKOiIiIFBwFHBERESk4CjgiIiJScP4/lzFR2lCzM+QAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Distribution Comparison between Copula-based Forecasting and Standard Path Forecasting\n",
    "#Sample from joint (path) forecast distribution using copula-based forecasting and standard sequential simulation\n",
    "#Take two example runs p=k=20 and p=k=5, both with nsamps=1000\n",
    "\n",
    "def sim_study_path_forecasting(p, k, n=10000):\n",
    "    a0 = np.ones(p)\n",
    "    R0 = np.diag(0.1 * np.ones(p))\n",
    "    pmod = pois_dglm(a0 = a0, R0 = R0,\n",
    "                 ntrend = 1,\n",
    "                 nregn = p - 1,\n",
    "                 deltrend = .99,\n",
    "                 delregn = .99,\n",
    "                 interpolate=True)\n",
    "    if p > 1:\n",
    "        X = np.random.normal(0, 0.1, [k, p-1])\n",
    "    else:\n",
    "        X = None\n",
    "    path_forecast = np.sort(pmod.forecast_path(k, X=X[:k].copy(), nsamps=n, copula=False).sum(axis=1))\n",
    "    path_forecast_copula = np.sort(pmod.forecast_path(k, X=X[:k].copy(), nsamps=n, copula=True).sum(axis=1))\n",
    "    \n",
    "    corr = pmod.forecast_path(k, X=X[:k].copy(), nsamps=n, copula=True, return_cov=True)\n",
    "\n",
    "    tmp = np.sqrt(corr.diagonal()).reshape(-1,1) @ np.sqrt(corr.diagonal()).reshape(-1,1).T\n",
    "\n",
    "    corr = corr / tmp\n",
    "    \n",
    "    return path_forecast, path_forecast_copula, corr\n",
    "\n",
    "\n",
    "path_forecast_20, path_forecast_copula_20, corr_20 = sim_study_path_forecasting(20, 20, 10000)\n",
    "path_forecast_5, path_forecast_copula_5, corr_5 = sim_study_path_forecasting(5, 5, 10000)\n",
    "\n",
    "# Plot comparison of forecast distribution from the two methods\n",
    "fig, axes = plt.subplots(ncols=2, figsize=(8, 4))\n",
    "axes[0].scatter(path_forecast_5, path_forecast_copula_5)\n",
    "axes[0].plot(np.array([0, np.max(path_forecast_5)]), np.array([0, np.max(path_forecast_5)]), color='r')\n",
    "axes[1].scatter(path_forecast_20, path_forecast_copula_20)\n",
    "axes[1].plot(np.array([0, np.max(path_forecast_20)]), np.array([0, np.max(path_forecast_20)]), color='r')\n",
    "axes[0] = ax_style(axes[0], xlabel='Path Forecast Total', ylabel='Copula Forecast Total')\n",
    "axes[1] = ax_style(axes[1], xlabel='Path Forecast Total', ylabel='Copula Forecast Total')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARcAAAEYCAYAAAB2hbIBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAmGklEQVR4nO2debwcVbW2n/ckOeckIQmYGIEEFAVERAZB5klBDTigggiIXBSNfF6c9Yqz4vAJzhdwCMJVQSZRuBEiIEhAGWRQRGYDDkkAgQQyz1n3j10dK50eVk/npPqsJ7/+5XT36rVXVVev2nvXW2vLzAiCIGg3PYMdQBAE3UkklyAIOkIklyAIOkIklyAIOkIklyAIOkIklyAIOkIklyAY4kg6T9KTku6t8r4k/bekWZLukfRyj99ILkEQ/BiYUuP9w4DtssdU4Psep5FcgmCIY2Y3AfNrmBwB/NQStwGbStqint/h7QowCILO86rPXtKwpP6GLx/zXlKPo8Q0M5vWgItJwOzc8znZa4/X+lAklyDocrJE0kgyaQuRXIKgQPRoUJqdC2yVez45e60mMecSBEE9pgMnZFeN9gYWmFnNIRFEzyUICoXU/q6LpIuAg4EJkuYAnwdGAJjZD4AZwOHALGAp8E6P30guQTDEMbNj67xvwH826jeSSxAUiEGac2mKmHMJgqAjRM8lCApEJ+ZcOkX0XIIg6AjRcwmCAlGgjkv0XIIg6AzRcwmCAtFToK5L9FyCIOgI0XMJggIRV4uCIBjyRM8lCApEgTou0XMJgqAzRM8lCApED8XpukTPJQiCjhA9lyAoEDHnEgTBkCd6LkFQIELnEgTBkCd6LkFQIOLeoiAIhjzRcwmCAlGgjkv0XIIg6AzRcwmCAhFzLkEQDHkGrOfyqs9eYh67LZ8z2uVvmCODr1qz1uVrdJ9vN3g1BitWrXHZjR3V21Z/3pPa6jWur4Llq1a77HqH1T9HrVjt+y62nrCJy272vMUuu/4Rvu92lPMYWLJ8lctuxWrfd/azj7yhoa5IgTouzfdcJF0j6SlJn2lnQEEQdAet9FxOAg4lrXgfBMEAMCQUumY2p52BBEHQXXR0QlfSVEl3SrrzsT9e18mmgmBI0CM1/Bi0WDvp3MymmdkeZrbHli8/tJNNBUGwkRE6lyAoEAWacmk+uUg6B9gX6JO0h5m9qW1RBUFQeJpOLmb2nnYGEgRBfYp0tWjAhkVecdxj85e47Eb31w99pFNA9cSCpS47rz+v6G3ZSp9Ibc1an+jNy1rz+esd7puSW+0QKy5cttLlq8/Z5qhe33cxf/EKl92CpT47565jWE9xkkCniDmXICgQRcpZrSh0d5N0s6SbJP1W0gvbGVgQBMWmlZ7L48AUM1sk6XDgi8A72hNWEASVUIHWLWplQveJ3NMVgG8CIQiCIUHLIjpJo4EvA1+v8N46he5fb7261aaCYMjTo8YfgxZrKx+WNAK4BDjdzO4vfz+v0N1unymtNBUEQcFoRUTXA1wAXGFmV7QtoiAIqlIknUsrPZe3AK8Djpc0U9KZbYopCIIuoJUJ3cuAy9oYSxAEdehUz0XSFOC7wDDgR2b2tbL3twZ+Amya2ZxqZjNq+RwwEZ2nLCX4lLcAS5bXvzg1dqS3jKQvNq+ytb93mMtuxHCfnbe0ovf2eq+d90D2qFGHO0phAqxx7uNVThV03wjfPvYqgxcu830XI50K4o0BScOAs4FXA3OAOyRNL5tH/QxwqZl9X9KOwAzgBbX8RoHuICgQHbpatCcwy8weNbOVwMXAEWU2BozN/h4HPFY3Vv9mrY+ksZJuyeZbbpd0SLO+giAYVCYBs3PP52Sv5fkCaX51DqnX8v56TlvpuSwGDjSzg4FjgK/VNg+CoFUkNfNYpzfLHlObaPpY4MdmNhk4HDg/u2JclVYmdNcCpdthxwL3NOsrCILOYWbTgGk1TOYCW+WeT85ey3MSMCXzd6ukfmAC8GQ1p62K6CZJ+j1wLXB5hffXZcyHbwmFbhC0itT4w8EdwHaStpHUSxqJTC+z+SdwSIpBLwH6gadqOW0puZjZXDPbnzQhdFaF99cpdLffNxS6QbAxYmargVOAa4AHSFeF7pN0mqQ3ZmYfBd4j6c/ARcCJZrUv7bWi0O0zs1KFnYXAomZ9BUHgo1PV/DPNyoyy1z6X+/t+YL9GfLZyMX4nSd8G1mR+PtSCryAIuoxWJnTvAg5sYyxBENShSPcWDZiM0LsovLdOrUd9+/gzvtq4E8b0u+y8KtPhw5zKVueBMmbkCF+7Pb74vIuke2rjAtS5IgnA6D7fNnjLBXtLCfQ6FbojnApdnApdr9K4mymORjkIgqFRQ7eEpO0lrZK0fzsCCoKgO2hHz+WzwI1t8BMEQR2GzJyLpL2AJ0hXjIIgCNbR6rDo09S4pyiv0J11Wyh0g6BVOqTQ7Qit3BX9OuBOM5tXzSav0N1271DoBsFQopVh0a7AwZL2BV4G7CDpbWb2j7ZEFgRBoWlFRPcV4CsAkn5MKo0XiSUIOkin5P+doC06FzM7sR1+giDoHgZMRDe6z9fUEwt8qlpP3Vuv8vbpRctddhPHjXTZLVnuu3jmqT0LsNKpqO111uRd45TBjnR+Z6sc8S1Z4VuQs8+5Dd7vYuGylS47r0LXWx+5UxRpOddW67ksy8pczpR0UruCCoKg+LTac5mblbkMgmAAKNCUS8s6l80l3Sjpl5Je0I6AgiDoDlrtubzAzJ6W9FrgXLIyeEEQdIYiXS1qtczl09n/1wDPL38/r9B98OZft9JUEAQFoxWF7ibZSm1I2hl4utwmr9DdYb/DWggzCAIolvy/lWHRjsAPJS0ircb23vaEFARBN9CKQvd2YLc2xhIEQR2KVHIh1ooOgqAjDJhC15txvTV01zpqlHpr3nrVnk8uWOay22KzUS67Ub2+bfXuO49SFmDZKp/dSGf92TrL1wDQ7/S12eg+l12fNzaf4Nvtb/FyXw1d775rlCFT5lLS7pKulXSDpDPaFVQQBMWnlUXRekmFot5iZrEgWhAMAENlzmUfYDFwoaTfSjqgTTEFQdAFtDLnsiWwC6lo1Bjgekkvqbd+bBAEzTNUFLrzgVvMbKGZzSWJ6J6bN8grdB/4/YyKToIg6E5aSS5/ALaXNFzSGGAisF493bxC9yX7H95KnEEQMEQUumb2rKQzgZnACOATZhZLjARBALSoczGz84Hz2xRLEAR1GCpXi4IgCKoyYArdFU5VqNfOU8t0+DBflvfWvPUqbx9/xicLnTx+tMvOW/PWW3/Wq1z2trtq9dq22AA4m2TRMp9S1qPkBljmrPHb69x3C5b6avc2SpEUuq2I6HYEvpc97QO2N7PxbYkqCILC08qE7v3AwQCSjgZe1aaYgiCowlCcczkeuKBNvoIg6AJannORNB7YAbi59XCCIKjFkFm3KONtwM8ryf7zCt2Hb7m6DU0FQVAU2pFc3k6VIVFeobv9vlPa0FQQDG161Phj0GJt5cOSXgj0mdkDbYonCIIuoVWF7qPAHm2KJQiCOgzFq0VBEATrMWAK3bGjel12y1b6lJIjHGrUYc4sP8w5MPXWvPUqb+fMW+Ky89b4XeGsoetV3nrr2Y7qq79fvIraJSt8dhPG9LfVn+d4Ali9xqc0Hu+Mr1GKpNBtZVE0STpL0q2S7pB0bDsDC4Jg4JA0RdJDkmZJOrWKzdGS7pd0n6QL6/lspefyUuClZrZPVs/lbuCiFvwFQVCHTsy5ZCunng28GpgD3CFpeqbCL9lsB3wS2M/MnpE0sZ7fVuZcHgNWShpBKnM5vwVfQRAMHnsCs8zsUTNbCVwMHFFm8x7gbDN7BsDMnqzntJXk8gzwV+BhUq/lyy34CoLAgaRmHuvErNljapnbScDs3PM52Wt5tidVnrxZ0m2S6grXWkkur84C2JYk//+qpPVmAPMbde9NV7XQVBAEzZIXs2aPaU24GQ5sR7pZ+VjgHEmb1vpAK8lFwDNZactFQC+w3pR7fqN2OvB1LTQVBAF0TKE7F9gq93xy9lqeOcB0M1tlZn8jjVi2qxmrf7M24DqgR9LvgVuAM828i2cGQbARcQewnaRtssUOjwGml9lcwb9LrEwgDZMereW0lXoua4ATm/18EASN04mrRWa2WtIpwDWk0cd5ZnafpNOAO81sevbeayTdD6wBPm5m86p7HUARXRAEGy9mNgOYUfba53J/G/CR7OFio6uh61WPLlleX3k5ZuQIl6+VTmWr96zh3Qav8vbJBctcduOcKuje4b7R8L+e9Y1yRzqUy4sd3xf46x4/66xR610AdPFynz93LWhnzeBGKdCtRS3fFf1FSbdImilp53YFFQRB8WmlQPeuwJ5mtq+krYCfAq9sV2BBEGzIUFkrenvgLgAzmw1sU65zCYJg6NJKcrkXOFhSr6RdSNfGN2tPWEEQVKJIa0U3nVyym5ouBH4DfBC4D3gqb5NX6D7w+xkVvARB0K20NKFrZt8zs4OAbwF/KV+IPq/Qfcn+h7fSVBAEpDmXRh+DRUuXoiVdm/mYB/xnWyIKgqAraLWG7mvaFUgQBPUZausWBUEQbMCAKXTbPfTzjCWH9/hyZ6+zfuoqp5K3z+nPW/PWq7xd4FSt9o3w7ZetJ4xx2S1ZUb/usVd5u0m/b1u9LHfWZPYeKyvlU95u6qw/3CgFkrnU77lIukbSU5I+kz2XpDMl/U7SlZKe0/kwgyAoGp50fRLw8dzz1wKjzOwA4FLgvzoRWBAEG1Kkq0V1k4uZzSl76SDgyuzvX2XPgyAI1qOZOZfxpPq5AM8SqtwgGDC6as6lAvOBTbO/x/HvRLMBodANgqFLM8nlRqAktz08e16RUOgGQXtppvr/YFF3WCTpHGBfoE/SHsBbgNdL+h2wEDihsyEGQVBE6iYXM3tPhZdD6h8Eg8BQqecSBEFQlQFT6K5e46tlutZZ89STwb0KWG/N22XO+qnDh/lytrddb81br/J2xSqfynT5St/2Ll5WXxns/FoZ1uNT1PY4F+RZ5lTojun31Vte6/zOvDV5G8W52RsFzSh0XyTpLkmLJe3f+RCDICginp7LScChpEpzAI+TlnL9VqeCCoKgMoN59adRPBO6c/IblK2quLRIGxkEwcATi6IFQYGIq0UZeYXugzf/upNNBUGwkdHR5JJX6O6w32GdbCoIhgRFqv7fjEL3BOCXwI7ASyXNMLPPdzbMIAiKRrMK3UM7EEsQBHUo0oWUUOgGQdARBuxq0fJVPqWkV43qyeCr1/iUqCP7fLth5AhfbVyv8nYzZ53Vfz271GXnrXnrVd7OnrfYZTdx3EhHm+1V3nrrGU8cWz82aL+qet6i5S67Rul2he4Jkm6XdJOki2N96CAIKtFMDd3fA/uY2YHAP4HjOxFYEAQb0lX1XCoodB/Nvb0C8PV3gyAYUjQ9oStpB2AKcEkNm3Uiulm3Xt1sU0EQZBSp59JUcpE0GfgJcIyZVZ25yovott1nSrMxBkFQQBq+WiRpAvAL4GQze6T9IQVBUI0iXS1qRqE7B5gEfDvrcp1vZud2NMogCApHswrdUzoQSxAEdSiSQnfARHS9TpGSV/g2zNE/lHxtegVZ5qzVuGq1bxtGecV7vT47z4Lw4CtLCT5xHMCTC5bVtRnZ6xMger8LL88u9W3raGeZy9H9vu+iQDmgY0Q9lyAoEEWac2lGoXuQpJsl3SjpBklbdT7MIAiKRjMK3VvNbD8zOwg4H/hARyILgmAD1MQ/l19piqSHJM2SdGoNuyMlWXZxpyZ1k4uZzSl7nh/EjgXuqecjCIKNF0nDgLOBw0h1mo6VtGMFuzHAB4E/ePw2K6J7naQ7gfcBt9awW6fQfeiWUOgGQat0SKG7JzDLzB7NOg8XA0dUsPsScDrguuW7qeRiZleZ2R7AZ4Cv1rBbp9B98b6h0A2CwSB/ks8eU8tMJgGzc89LWra8j5cDW5nZVd52m1Ho9uck/88CvmIjQRC0TDNXi8xsGjCt2TaVNB3fAk5s5HPNKHSvlPQOYC3pruj3NhxtEAQbE3OB/FXfydlrJcYAOwEzs2HW5sB0SW80szurOW1WofsjT8RBELSXDil07wC2k7QNKakcAxxXetPMFgATcjHMBD5WK7HAAIroVjhVqwud6lFPWcLRfT7VpVfZ2u8sc+lV6C5atsplt3i5z274MN+B510U3lua0qO+XeYsrbmJUyk7Yrjvu1ixyteuVxnuXYh+RZuVxp3EzFZLOgW4BhgGnGdm90k6DbjTzKY34zcUukFQIDp1W4GZzQBmlL32uSq2B3t8NqzQzb3+Tkm+U2oQBEMOT8/lJNI6RZNLL0jqB44k1dANgmCA6Kq1ossVuhkfAH4AOEfvQRAMNRoW0UnaDDjQzK502P67hu5todANglYp0lrRzSh0Pwmc4TFcr4bu3qHQDYKhRDPJZXvgU5KuBraQVLX6fxAE7aVHavgxWDSs0DWzN+Xem2Vmb+tgfEEQFJRmFbql97ZtbzhBENQiauhWYOsJm7js+pwL0a9xyEydYkr6nGpP78Lx3naXrGiv8naT/l6X3bCegV8U3qu8fWqhbwH3Fzx3jMuu31l/2PvdehW/3rrH3UwzZS4PlvS4pJnZY/fOhxkEARTralFTIjrgKjN7d2dCCoKgG2hWRPdaSb+TdKYk3/oTQRC0TJGuFjVzKfouYDszOwBYCHysvSEFQdANNJxczGxRrhLdz4CqVcDzCt27Z9YV9AZBUIcizbk0I/8fl3v6KuCharZ5he6uB7++mfiCICgozZS5vFbSu0i1c58G3tXZEIMgKOFdh2hjoFkR3fc6EEsQBF1EKH2CoEAUqZ7LgCWX2fMWu+xGOZWNqxxKSe8yDBPH+a6m9zlr6Hpr404Y0++ye3apr66wl3Yqb714a956lbd/f2qRy27L54x22a0chBrP3U5TZS4lnSDp+mwh+uNqfT4IgvZRpKtFDSt0Jb00e36ombeOfBAEQ41mFLpHka4UXSvpckmTK3wsCIIO0KG1ojtCMwPDLYHxwGuAc4FvtDWiIAi6gmaSy3zgmmxIdA3wsmqG69XQvTVq6AZBq/So8cegxdrEZ2byb8n/7sAj1QzXq6G7T9TQDYKhRDMK3TcDU7L1YnuAqR2NMAiCdXRVJboqCt0PdyCWIAi6iFDoBkGBCIVuBfpH+Jqav3iFy86jlu11Kmq9qktb6jJjrVP+462h65UTLV/pq427zGk3caxPuexREHtrz3pr3nqVt4/NX+Kymzze529Uny8+bx3lbqaZGrrH5ernPiDpF50PMwgC6HKFrpldCFwIIOl7wE0diy4IgsLSbA1dJI0ADgP+t91BBUFQmW5X6JY4DLjJzJa1K5ggCLqHVpLL8cAFtQzyCt0Hb/51C00FQQDdr9BF0liSOvf6WnZ5he4O+x3WTFNBEBSUZheiPwq4wsx8FXaCIGgLXa/QNbPzOhNOEATdQih0g6BAFKjjMnDJxatsXLDUqdAdXn+6aITDphE7bw3dZSt8ClhvXdnFy511W3t82zGmf4TPn7MO7GiHv9VrfCPozUb3uey8NW+9yts583xK3hc9b6zLbk1IdJtS6G4m6VpJN0q6WdLOnQ8zCALovrWiTwI+nnv+duBmMzsI+HT2CIIgWI9mFLoPAKW+4WbAk+0OKgiCyqiJf4NFMzqXu4C9Jd0LnAl8s70hBUEw0EiaIukhSbMknVrh/Y9Iul/SPdmyQs+v57OZ5PJfwC/MbCfgrcDZNQJep9C996armmgqCII8nVDoShpG+h0fBuwIHCtpxzKzPwF7mNnOwGXAGXVjbXTjAJEWoIc0JHpONcO8QnenA1/XRFNBEAwAewKzzOxRM1sJXAwckTcwsxvM1lU0uo2sSkItmqmh+z7gfEnvAkYCn2hoM4IgaJoOKXQnAbNzz+cAe9WwPwmoe7NgszV0D6n3uSAINg4kTWX9QvrTzGxak76OJ63+cVA921DoBkGBaKbjkiWSWslkLrBV7vnk7LWytnUoSXpykJnVVbsOWHJZstxbL9bnb+Eyhz+PDdDf61XK+vz1OpWtXtWqt/7sSvn8rXWqR70q09H99Q8jb5vebfXWPfYqw73K20f+tdBl99yx/S67jYQ7gO0kbUNKKscAx+UNJO0G/BCYYmYu+UkzCt1Rki7LauheLmnTBjckCIIm6YRC18xWA6eQVlB9ALjUzO6TdJqkN2ZmXwc2AX4u6W5J0+v5bbiGLmnsdqeZfU3S20jq3VDpBkGBMbMZwIyy1z6X+/vQRn02o9DdHrgz+/t24JWNNhoEQXN0ew3dvwClhZ8Pp4bOJQiCoUszyeVcoF/SDaTr449VM8wrdB+6JWroBkGrFGndooaTi5mtNLNTzOyVwN9JUuBqtusUui/eN2roBsFQohmF7qeA7wFrgHtYvxxDEAQdpKvWiq6i0D24/aEEQdBNhEI3CApEgTouA5dcVqz2KS+HOVdxGtlbP/Q1Xrmvk5HOGroLlvrUo+PH+FScK5z1Yjd11p/11uSdt2i5y85zwHu/f8/3Cv76vt5Stl41sld5+9RC377rZjxzLrsBZ5HmWFYD7yZdIToX2Br4J3CSmcXeDIIOU6Q5F0/6f5x0P8GBwDeALwInAg+a2QHAQ9nzIAiCdXgUuk+Y2aLs6QpS7+Ug4MrstV/huP06CILW6UqFrqTRwJdJNzCNB57J3nqWUOkGQVCGK7lIGgFcApxuZvcD84FNs7fHZc8rfW6dQnfWrVe3IdwgGNp0lUJXUg9wAWnh+Suyl28k3VdE9v+NlT6bV+huu8+USiZBEHQpnut+bwFeBzwvK3H3F9IKAOdJ+h2p3uY7OxdiEAQlegZxHaJG8Sh0L6Py/UPHtj+cIAi6hVDoBkGBKJDMBcxs0B7A1HbZtdNX2MV3Nph2tR6PzHnSGn202mazj2bqubSTqfVN3Hbt9BV2A2O3Mcc2mHZV6UqdSxAEQSPEnEsQFAjnfb0bBYPdc/Gu+uaxa6evsBsYu405tsG06wqUTTQFQVAA/vnEvIZ/sFtvPn5Q+juD3XMJgqBLiTmXICgQg3n1p1EGpeeiOntIiZ0k7VDHrkfSoZJe47B7u6QxDrsPSnq+w+7rknZx2N0i6at1bI6UdFQdX5K0u6S9atl1kuw+s3b685X26zLqHf/dwoD1XLIduoeZ3WFmJklWYcInO4AvItWO6ZV0o5l9v4rdxaSSD9tI2sHM/rtK80cA/wOcIOl6M3uqir/zgTvM7B91tuNy4Doz+7OkEWa2wQr1mb9pwD+y7djSzB6r4OtK4BFgL0n7mdmHq7R5FfBX0oLhj5jZ+6vFWP7ZSvt5MCnFZGZrJA0zM18NzI2UWvs4Ow4+C9wPPGNm19U6/usRV4vKyHbwjcCnJB0NUNrBFczPBf5mZicA3yHVjqnEucAcM5sKvIu09MlLq9jeB9xBWvP6sKw2TTn/C4w1s+9IOjPrmXypgt044HYzO1PSRcBPJZ1V4aw+DXiKdFPneGA32OCsdSIwP0sUBwCrJL2kQptvAh4zsw8CrwdeLemnVbZ1HVnPoL9Cu26yntWPJJ0F/D9Jm9ewlaRRksbVcXu+pD8ClBJMDX871evhZHZvkLSJY1vOkPRih7+JkibW2Y7Sfu3L/V3OT0jHzBbAUZI+C+n4r+e76AzUsOiFpFq7PwJeJumtsC7BrIshO4jON7NTs5eeBQ6ukgxOM7OPZH9/FTiGdKf2eeWGZvYw8HPgJlLVvLco1ZrJV1v+PPBiSY8Cs0i9mP0lfavMnYBXSrqM1Os4CZgInJ3bjh7gLDP7pKXawlcBn5e0VdlB9RdgsqRDss+/CjhL0gVlbc4CJkrazczWAl8C9pD00Qr7pfQjuhQ4E/iGpAPLk3lm835JW1fykdmUemlzScl3K+C51dok1fw5E/iCpEk17O4HJkn6paThWYIZXsHul8A+tXo2ud7flma2uOz1cn8XAK8Btsle2yBpZXaXk46pc5RqSFfcjux7+hZwhqQDKuzjTUi1p0/PetX/n9Tz/Eglnx4UCt31MbNZZnY8cAPpwNqt1IMhO7NmdmuA38G6AlWPAQvMbImkKZIOz9n+LbMbC9xtZrub2V7AMEmbluyys5BI1fLuJh00nyD1Gobl/P0ROAo428y+a2b3AG8GxuQPfDN7hlSwfGtgqZktNbOjgS1LZ3UzW2tmd2ftD7d0Z/m1wPOz10rt3gV8gVTWYmtLtW8OAcZL2jG3Cx8CrgY+LekTwHHAh6n+/X0SeJS0YN11wNmSXluW2K7L2n5rjQSzE3CzmX3ezH4DLAWOLzfK9u9lwMPAaaTvbYyk3tz7pf23lpSEPgDMBmZkb21e5u83pB7iOZLeKuk1qtyrO4I0RP2hpI9Kep+k/fM/9Oz/C0n78WvAR7JYKiWtHwF/NbN3A9cAY0vfV1nC+i7phPk54PfA9/NJXNLOwItICwfukrX3D+B00gllywptdxUDerXIzJZKmgGsBV4o6fvAOEknmtnKzGZV7v9Vkh6RdDLph//eCj4XAt8EUJoU3YJU57f0vmXvXQEcCWxJSjLzSGeVvK97SAdDicNI62EPz/sEfgvsTPphLgOM1PVdWiG+0ucWkFarPLx0UGexzZQ0l3T2ezHprDoGeDLnY6WkH5OGd9sCJwMvAHbKEt+assTxL2C4mS0BLpe0EPiypKfM7I/Zj+VbwG3Ap4G3S7rIzP6e7SsBu5rZnyT9U1JPlhRuB16R2bwUeMLM5pESw+Vmdn723muBlwMm6UzgNkmvMLM7svjGAK83s3dI+q6kUq/2tOz955DmlyZJug54AFgC7CDpC6XEnbGU1Bu4CLgXWET6oX/QzH6b2bwMuN7MzsnimyJpqplVErX9nVSjCOAQ4JXAhyT90Mzy5RT/AfzdUn3pSyW9CfhidqyeRpozXAlMBvaW9Axwl5ndl538NqXGOuvVKNKcy4BfijazBcAlSmP4fYFjSoklT3aA9wFHk4ZHR5nZI9X8SnoL6Ww4Nd89zjELmADMNrPTpNoTapKOAz4EnGBly6aY2UKl4dL+pB7QIuD9WaLbYDuydr5JGnZtbmZPlJnNBWaSzoITgPea2dNlbS4lzVvdKOlQ0tDonaXkle2vPc3sD6T5pRdJ2hV4wMyul7Rt1v6w7Ed+Zfa500kTjsdK+i7pLPsd4AlJ52e9rhJPZ6+/HvgP0hzM7mZ2F2kYiaRXkpLbV0nf716ks/W8LIFdamZ3S/pTFtMqYBkpUQp4uZndpXSF7ZPAldk8WC/wPmBHSX8GXmFmtwN/JJ0ElprZV7IYniYNaW8gS5JkJ42sjRmkkwO510r77i7gFZKuJZ00jiIlmDdniW73zG52ZrcSWEg6WT0IfBu4x8xOVbqytxup3vRHgZuzHvk22We6mkHRuShdEt4CONrMHqpkk/0gl2dJ6Eozm1XD33DSQfrubH6lkr+Fkj5rZqV6vyIdPJX8jQCWkxLLg1X8LQauljQzC3dFte3IEsxaSSfnejJ5m6WSziZd/cLM/lVtWzMeBt5R2idK8wQzgfmSfmJml0t6NWke6irSUHMTUs/pT5IuNLNLs7aekPRFUiK9gNRrmAOcRzrjysx+nrW7GXAG6Qd9CmlOZH7eH3Cvmb09i+tVwI6k4cNFmb8eM7sYWEwaqswwsxdLei4peT4j6adm9gtJ/5XZlXpvm2XbcSMpWf3MzC6T9DBpHu1IM/sFqRc5MbObX7a9ln1np2afuyC3784zs+nAVZJOAZaY2SKlOb9RpB7rfEnnmtklSnMqe5PmFD8PHEiWWLJ9sRg40sxenbW1VWZ/spmVekcNoUGcQ2mUQUku2Rd2XLUfZBnfrjI2zvtbLenKWj2RzG4+rOtNVF3G0MxWSbq8nr/Mtu5icCU/lRJLzmY16YxfFzP7Z9lLpQnzi4B9JS0nDXveQzrjvo80R/Rn4BzSj9xKScPMnpQ0B3gr8AYze1DSKGAsaX5M2Y9zAekHdiopMeeTBlmv5CkASceQ5kNOKfO3i6TFZjZN0ioz+59sGzbN+dsna3Ndr0nSsaQ1yk8Ddsht6zLgB6STy96S3kzq/X0n83kR6TJ/+fZ+lTR026Gs3f7ctu6vNM+zO2nVi+Mzu/0krTWzc7PYhpN6N2/P9nnpBDWbLDmSkvazZvapWt9tNxH3FnUR2Q/4CGBX4HdmdqWkycBI0rBhbs7mZaSJ8Euz1/YjXdp/IOdvHDCF1I2fTLpS9DEzm12hzby/nUnDiYuq+NuOdHn+eWZ2nP49p1PN3ytIw+Mfm9l9NbZ1EmkovdwyTVEVf/3AaKAnlwxLdjsDt5jZryQdQUpas0o94pzdLpm/i7Ne12uBP5Xiy23zGcDfSPN9J9fqgXt4cv6Chn+wE58zblC6O5FcuozcD/hFpIQwjjS8W1PBZhvSxPAY4ESrIAbM7M8i9RqOrDSMrZCEJpDmoDYQK+b8HUCab3ugwvt5f5NISW0Df55trbG9/1Hek8zZbUsaVk1w+Nua1EN6dzYvVrIpzRk+yL/nDFtKLFCs5BI3LnYZZrbAzC4hXRXbD/hy+Y8jZzMZ2CezqZZYSvNjFRNLBX8HAV+okVjy820bJJYK/g6u5s+zrTW2t9LcV8lui6xdj7/9gS/lE0tmY9mQ+SxSEm05sUC6WtToY7CIGxe7EM+EuccG/PNjOX9HDaQ/73Z00O6ttexwzBl2KzEs6lIk9Tl+wHVt2t1mJ/xt7HbtZN4zCxv+wY7fbGwMi4L24Tno2/3DGCx/G7vdUCWGRUFQIAokc4meSxAEnSF6LkFQIHoK1HWJnksQBKWbOR+SNEvSqRXe75N0Sfb+HyS9oJ7PSC5BUCA6oXNRukv+bNINoDuSbmLdsczsJFIlvW1JN2eeXjfWRjcuCIKuY0/SLQ6PWqpQcDHpFoc8R5Cq6kGq3XOI6txFGXMuQVAgxo4d0/Cki6SprL9O9TRbv5bNJNJNliXmkEplUMnG0o3CC0j3hz1NFSK5BEGXkyWSAV/tMYZFQRDMJdWaKTE5e62iTVZiYhypQFZVIrkEQXAHqVToNkoV/44BppfZTCdVH4RUTuO3VufeoRgWBcEQJ5tDOYVUkHwYcJ6lWr+nAXdaqs53LmlJmFnAfFICqkncuBgEQUeIYVEQBB0hkksQBB0hkksQBB0hkksQBB0hkksQBB0hkksQBB0hkksQBB3h/wChQkVYyUTg/wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 288x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot correlation matrix from simulation of k=20 forecast horizon\n",
    "labels = [str(i) for i in range(1, 21)]\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(4, 4))\n",
    "axes = sns.heatmap(corr_20, vmin=0, vmax=1, center = 0,\n",
    "                   cmap=sns.diverging_palette(10, 240, as_cmap=True),\n",
    "                   cbar=True,\n",
    "                   square=True, ax=axes,\n",
    "                   xticklabels=labels,\n",
    "                   yticklabels=labels)\n",
    "axes.set_xticklabels(labels = labels, rotation=45, size=9)\n",
    "axes.set_yticklabels(labels = labels, rotation=0, size=9)\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e70966c555988d806fbf33aba289056abaf362856512e47686ded46e4756b9ad"
  },
  "kernelspec": {
   "display_name": "Python 3.8.9 ('venv': venv)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.9"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}