Using Pytorch for group based trajectory models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using pytorch to estimate group based trajectory models\n",
    "\n",
    "[Andy Wheeler](mailto:apwheele@gmail.com)\n",
    "\n",
    "This code illustrates using pytorch to estimate group based trajectory models. The idea is basically we have some outcome that is a function of time, so say a regression equation of the form:\n",
    "\n",
    "$$\n",
    "\\begin{matrix}\n",
    " y = 1.0 + 0.5(t) + 0.0(t^2) + \\epsilon_1 & \\text{: Group 1 Model} \\\\ \n",
    " y = 0.0 + 0.0(t) + 1.5(t^2) + \\epsilon_2 & \\text{: Group 2 Model} \n",
    "\\end{matrix}\n",
    "$$\n",
    "\n",
    "But we don't know whether an observation belongs to group 1 or group 2 (the latent classes or clusters), and we also do not know the actual equation for the time trajectory (so we need to estimate the coefficients that we multiply by $t$).\n",
    "\n",
    "Here I am going to illustrate using [*PyTorch*](https://pytorch.org/) to estimate these models using a simple, simulated example. First the upfront libraries I need, and setting the seed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "####################################\n",
    "#Set up needed to run this code\n",
    "\n",
    "import numpy as np\n",
    "np.random.seed(10)\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "torch.manual_seed(0)\n",
    "#if you don't have a GPU change this to cpu\n",
    "#it runs pretty fast either way with this small of data\n",
    "#device = torch.device(\"cpu\")\n",
    "device = torch.device(\"cuda:0\")\n",
    "####################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next I simulate some data. Here I do two groups (same as the initial equations I showed earlier) with 50 individuals in each group and 11 time points. To note I am doing the simulations in long data format (many gbtm models make you have the data in wide format). This makes it easier though to pipe in new variables, and it also allows you to do unbalanced panel data in this format as well. I include `t0` (a constant term), as this allows me to estimate the intercept in the subsequent pytorch model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>ID</th>\n",
       "      <th>group</th>\n",
       "      <th>t</th>\n",
       "      <th>t0</th>\n",
       "      <th>t1</th>\n",
       "      <th>t2</th>\n",
       "      <th>outcome</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.527132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.818182</td>\n",
       "      <td>0.669421</td>\n",
       "      <td>0.542984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.636364</td>\n",
       "      <td>0.404959</td>\n",
       "      <td>0.695183</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.454545</td>\n",
       "      <td>0.206612</td>\n",
       "      <td>0.797608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.272727</td>\n",
       "      <td>0.074380</td>\n",
       "      <td>0.863487</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1095</td>\n",
       "      <td>99</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>0.090909</td>\n",
       "      <td>0.008264</td>\n",
       "      <td>0.061558</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1096</td>\n",
       "      <td>99</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>0.272727</td>\n",
       "      <td>0.074380</td>\n",
       "      <td>0.077348</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1097</td>\n",
       "      <td>99</td>\n",
       "      <td>1</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>0.454545</td>\n",
       "      <td>0.206612</td>\n",
       "      <td>0.301391</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1098</td>\n",
       "      <td>99</td>\n",
       "      <td>1</td>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>0.636364</td>\n",
       "      <td>0.404959</td>\n",
       "      <td>0.626088</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1099</td>\n",
       "      <td>99</td>\n",
       "      <td>1</td>\n",
       "      <td>10</td>\n",
       "      <td>1</td>\n",
       "      <td>0.818182</td>\n",
       "      <td>0.669421</td>\n",
       "      <td>0.989794</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1100 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      ID  group   t  t0        t1        t2   outcome\n",
       "0      0      0   0   1 -1.000000  1.000000  0.527132\n",
       "1      0      0   1   1 -0.818182  0.669421  0.542984\n",
       "2      0      0   2   1 -0.636364  0.404959  0.695183\n",
       "3      0      0   3   1 -0.454545  0.206612  0.797608\n",
       "4      0      0   4   1 -0.272727  0.074380  0.863487\n",
       "...   ..    ...  ..  ..       ...       ...       ...\n",
       "1095  99      1   6   1  0.090909  0.008264  0.061558\n",
       "1096  99      1   7   1  0.272727  0.074380  0.077348\n",
       "1097  99      1   8   1  0.454545  0.206612  0.301391\n",
       "1098  99      1   9   1  0.636364  0.404959  0.626088\n",
       "1099  99      1  10   1  0.818182  0.669421  0.989794\n",
       "\n",
       "[1100 rows x 7 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "##############################################\n",
    "#Simulate data\n",
    "n_obs = 50 #per group\n",
    "periods = 11\n",
    "\n",
    "#Group coefficients\n",
    "gc = [(1.0, 0.5, 0.0),\n",
    "      (0.0, 0.0, 1.5)]\n",
    "gc_np = np.array(gc).transpose()\n",
    "gps = len(gc)\n",
    "\n",
    "#Creating panel data over time in long format\n",
    "t = np.tile(range(periods), n_obs*gps)\n",
    "id_col = np.repeat(range(n_obs*gps), periods)\n",
    "group_col = np.repeat(range(gps), n_obs*periods)\n",
    "t_dat = pd.DataFrame( zip(id_col, group_col, t), columns=['ID','group','t'])\n",
    "\n",
    "#Latent classes (long format)\n",
    "lat_np = pd.get_dummies(t_dat.group).to_numpy()\n",
    "\n",
    "#Prepping the time design matrix\n",
    "center = periods/2\n",
    "t_dat['t0'] = 1\n",
    "center = periods/2\n",
    "t_dat['t1'] = (t_dat['t']-center)/center\n",
    "t_dat['t2'] = t_dat['t1']**2\n",
    "vars = ['t0','t1','t2']\n",
    "t_np = t_dat[vars].to_numpy() #this is our design matrix\n",
    "\n",
    "#Now creatings the outcome data\n",
    "pred_full = t_np.dot(gc_np) #This gives us predictions for both groups\n",
    "y = (pred_full * lat_np).sum(axis=1) #Element wise multiplication here\n",
    "\n",
    "#Add some error to the true underlying curves\n",
    "y_werr = y + np.random.uniform(low=-0.05, high=0.05, size=y.shape)\n",
    "t_dat['outcome'] = y_werr\n",
    "\n",
    "#Show what the data looks like in the end\n",
    "t_dat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows the two groups over time (you can clearly pick them out from the plot) and the typical error. So this should be a fairly easy problem to solve (much easier than real data). This also prints out the mean square error for the true data, so I should not be able to fit a model that has any less error than this (if I am, it is over-fitting). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFzCAYAAADxBEqxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1Tc550v/vczjaGKJiSKABVA9A6SkIS6UJctK5Zjx47jXSeb+G7J3bv33rO/k83Ze37Zm92z+SW7KY7s2HKJI8u2HMm2ZDUEqIAomqGIKlTpVSDKMO35/YFhQQLmOzPfafB5neMTwzzPfD+REZ952udhnHMQQgghxDlJHB0AIYQQQmZHiZoQQghxYpSoCSGEECdGiZoQQghxYpSoCSGEECdGiZoQQghxYjJHPTgwMJBHRkY66vGEEEKI3VVUVPRwzheb08dkomaMvQNgL4AuznnCLG02AfglADmAHs55rqn3jYyMRHl5uTmxEkIIIS6NMXbf3D5Cpr6PAcib46G+AH4LYD/nPB7AYXODIIQQQsjMTCZqznkRgL45mnwbwEnO+YNv2neJFBshhBCy4ImxmSwagB9jrIAxVsEYe1mE9ySEEEIIxNlMJgOQDmArAHcAxYyxEs5545MNGWOvA3gdAMLDw0V4NCGEEDK/iTGibgHwNed8mHPeA6AIQPJMDTnnRznnGZzzjMWLzdr0RgghhCxIYiTqUwA2MMZkjDEPANkA6kR4X0IIIWTBE3I8608ANgEIZIy1APgnjB/DAuf8Tc55HWPsawBVAIwA3uac19guZEIIIWThMJmoOecvCGjzbwD+TZSICCGEEDKJSogSQgghTowSNSGEEOLEKFETQgghTsxhiVqj0Tjq0YQQQojLcFii7uvrw9DQkKMeTwghhLgEh059nzx50pGPJ4QQQpyeQxP13bt3YTQaHRkCIYQQ4tQcvpmsoKDA0SEQQgghTsvhifrq1auODoEQQghxWg5P1JxzNDU1OToMQgghxCk5PFEDwGeffeboEAghhBCn5BSJemxsDIODg44OgxBCCHE6TpGoAeD99993dAiEEEKI03GaRN3b2wuDweDoMAghhBCn4jSJGqC1akIIIeRJTpWo6+rqHB0CIYQQ4lScKlEDdK6aEEIImcrpEvXly5cdHQIhhBDiNJwuURuNRjQ3Nzs6DEIIIcQpOF2iBoATJ044OgRCCCHEKThlotZqtejt7XV0GIQQQogoDAYDTp06ZVFfp0zUABVAIYQQMj88fPgQ//7v/w61Wm1Rf5nI8YhmcHAQQ0ND8PLycnQohBBCiNnGxsbw1Vdfobq62qr3cdoRNQD86U9/cnQIhBBCiFk456iursYvfvELq5M04MARtY+Pj8k2bW1tGB4ehqenpx0iIoQQQqzT1dWF06dPo7W1VbT3dFiiFjqlfebMGRw+fNjG0RBCCCGWGx4exqlTp9DU1CT6ezt0jXrLli3Iz8+fs01tbS1GRkbg4eFhp6gIIYQQ0zjnaG1tRX5+Pu7evWuz5zg0Ua9Zs8ZkogaAoqIi5OXl2SEiQgghZG6PHz9GVVUVKioq0N/fb/PnOTRRy+VybNiwAVeuXJmz3Y0bN7Bp0yYolUo7RUYIIYT8F71ej4aGBlRWVuL27dvgnNvt2SYTNWPsHQB7AXRxzhPmaJcJoATA85zzT4UGkJOTYzJRA0BxcTE2b94s9G0JIYQQq3DO0d7eDrVajerqamg0GigUCrsmaUDYiPoYgF8DmLUCCWNMCuDnAM6ZG4CbmxvWrFmDkpKSOdsVFRVh3bp1cHNzM/cRhBBCiGBDQ0OoqqqCWq1Gd3c3ZDIZli5ditbWVmi1WrvHYzJRc86LGGORJpr9NwCfAci0JIgNGzaYTNQAUFpaig0bNljyCEIIIWRWer0ejY2NUKvVk1PboaGhiImJQVNTE1paWhwWm9Vr1IyxUADPANgCE4maMfY6gNcBIDw8fPL7Hh4eSE1NhUqlmvNZly9fRnZ2NhQKhbVhE0IIWeBmmtr29vZGdHQ0hoeHHZqcpxJjM9kvAfxPzrmBMTZnQ875UQBHASAjI2PaJP+mTZtMJmrOOcrKypCTk2NdxIQQQhasmaa2g4ODAQAtLS1oaGhwcITTiZGoMwAc/yZJBwLYzRjTc87/bM6b+Pj4IC4uDrW1tXO2KygoQFZWFuRyucUBE0IIWVhmmtoOCAhAUFAQenp68PDhQ0eHOCurEzXnfPnEvzPGjgH40twkPWHr1q0mE7Ver0dFRQXWrFljySMIIYQsEFOntmtqajA6OgqlUgkvLy8MDQ25zHXKQo5n/QnAJgCBjLEWAP8EQA4AnPM3xQzG398fK1aswJ07d+Zsd/nyZaSnp9OomhBCyFMmprYrKyvR1dUFxhhksvF0p9FooNFo7B6TUqmEVCq1qK+QXd8vCH0zzvl3LYpiip07d+J3v/vdnG20Wi1UKhWysrKsfRwhhJB5YGJqW6VSobm5edpZZ845dDqd3WNijE3GYc2HA6e7jzooKAihoaEmbx4pLCxEWlra5KckQgghCwvnHC0tLSgoKMC9e/dgNBqfaiORSMAYg8FgcEh8YnDKLLdr1y68/fbbc7YZGRmBWq1GRkaGnaIihBDiDHp6evDZZ5+hs7NzxmQok8mg1+sBYMbk7WqcMlGHhoYiKCgIXV1dc7YrKipCamqqxfP+hBBCXAPnHCqVCpcuXcLIyMhTr09Nzo4YPc/Gz88PcXFxiIqKgr+/P37605+a/R5OmagBIC8vD++/P2vVUgDjN5hUVlYiLS3NTlERQgixp8HBQVy4cAG1tbVzjo4nkjQg3pSzpUJDQ5GVlYXY2FhRNj07baKOjIzEokWLMDAwMGe7q1evIiUlBRKJxE6REUIIsSWDwYBbt26hqKjIZY5Q+fn5ISMjA0lJSfDy8hL1vZ02UTPGsGvXLhw/fnzOdv39/aiurkZycrKdIiOEEGILXV1duH79Ompqapxq+nrq7u2p3N3dkZKSguTkZCxZssRmz3faRA0A0dHRkwfT51JUVITExEQaVRNCiIvRaDSorq5GSUkJ+vr6HB3ONBMJemqSlkqliI2NRXJyMlasWGGXvOPUiZoxhh07duDkyZNztuvr68OtW7eQmJhop8gIIYRYinOOu3fvoqKiAvX19SZ3Zsvlcoecg56aoCMiIpCcnIy4uDi7X7fs1IkaAOLj43HmzBmTh8WLioqQkJAAUxeDEEIIcYxHjx5BrVajoqLC5EwpYwxKpRKjo6MOSdLAeLXM5ORkJCUlwdfX1yExAC6QqCUSCbZu3YqvvvpqznY9PT2ora1FfHy8nSIjhBBiik6nQ319PVQqFe7evSu4H+cco6OjNoxsuolpbnd3dyQkJCApKQmhoaFOMfhz+kQNAKmpqTh//rzJT1VFRUWIi4tzij9YQghZqCYuw1CpVKiqqoJWq3XKPUQTyVkikSA6OhrJycmIiopyutocLpGopVIpcnNzcfHixTnbdXV1ob6+HrGxsXaKjBBCyITh4eHJe54nLsOY4IwVwkJCQpCcnIz4+Hh4eHg4OpxZuUSiBoCsrCzk5+eb/I9dVFSE1atX06iaEELswGg04vbt21Cr1aivr58coQKOLzwyEx8fn8l158DAQEeHI4jLJGq5XI5169bh6tWrc7br6OhAU1MToqOj7RQZIYQsPD09PVCr1VCr1RgeHp521tjZRs9SqRQJCQlISUlBRESEyw3kXCZRA0BOTg6uX79u8oegsLAQUVFRLvcfgxBCnNnY2BgqKytRVlaGnp6eaa854+jZz88PGzduRHx8vCilPB3FpRK1UqlEeno6ysrK5mzX1taG5uZmrFq1yk6REULI/MI5x+PHj9Ha2gq1Wo2HDx/adRe2NcLDw7F//34EBAQ4OhRRuFSiBoDc3FyTiRoYH1WvXLmSRtWEEGJCQ0MDSkpKMDIygpGREWg0mmmXXLiKhIQE7NmzB0ql0tGhiMrlErWnpycSExNRXV09Z7uWlhbcvXsXK1assFNkhBDiWvr6+vDpp5+ivb3d0aFYJSoqCgcPHnTqndvWcLlEDQDbtm0zmaiB8R3glKgJIWQ6rVaLU6dOoba21tGhWMXDwwOvvfYa/P39HR2KTblkovbx8UF0dDQaGxvnbHf//n3cu3cPkZGR9gmMEEKcmE6nw9mzZ6FSqRwdisUCAwORnp6O5ORkuLu7Ozocu3DJRA0AO3fuNJmoGWMoKiqiRE0IWdB6enpQUFCAW7duifJ+EonErhdpLFu2DOnp6YiJiZl3689CuGyi9vf3R3h4OB48eDBrm4kbWh48eIDw8HA7RkcIIY6l1+tRV1eHa9euobOzU1CfqWehZ7uDGRB2TtraJL1kyRJkZmYiPj5+QSbnqVw2UQPArl278Pvf/37ONhOj6pdeeslOURFCiOP09PSgrKwMKpXK7GTp6+uL4eFhaLVah5yL9vb2RlpaGrKzsxfMtLYQLp2oly5diuDg4Dl3LHLO0dzcjJaWFoSFhdkxOkIIsQ+9Xo/a2lrcuHEDbW1tZvVdtWoVFAoFGhoa0N/fb/azpVIpwsLC0NraatGRLplMhujoaGzdunXebwqzlEsnamB8VP3OO+/M2WZiVP3tb3/bTlERQojtdXd3o7y83OzRs1wux759+9DX14eCggKT7aVSKWQyGcbGxia/5+3tjYyMDJSUlOD+/ftmxx4UFIQNGzYgLi7OKW/WciYun6iXLVuGgIAA9Pb2ztqGc46mpia0t7cjODjYjtERQoi4dDod6urqUFJSYtH5Zz8/PwwNDeHkyZMzvi6VSrFs2TKsWLECOp0OTU1N6OjogMFgAAC4ubnB398fHR0duHz5slnPlsvlSElJQXZ29rypGmYPLp+ogfFR9YcffmiyXVFREZ5//nk7REQIIeLq6upCeXk51Gq1VRu1Zpre9vDwQHJyMtauXQulUgm1Wo3i4uKn2kokEoyNjZn9AWHJkiXIyspCYmKiS9fcdpR5kahXrFgBHx8fDA4Oztmuvr4eHR0dWLp0qZ0iI4QQy+l0OtTW1qKkpAQdHR2iv39kZCS+9a1vwd3dHUNDQygrK0NZWdmsNb3NuRVLKpUiPj4emZmZCA0NpXLOVjCZqBlj7wDYC6CLc54ww+svAvif33w5BOCvOOeVokZpAmMMO3bswKeffmqy7ZUrV3D48GE7REUIIZbp6upCWVkZ1Gq1TWpuBwcH4/Dhw9DpdKioqEBtbS06Ojpm3ek911GtJy1atAgZGRlITU2Fp6enmGEvWEJG1McA/BrA+7O8fhdALue8nzG2C8BRANnihCdcXFwcPDw8MDIyMme72tpadHV1ISgoyE6REUKIaTqdDjU1Nbh27dqce26soVQqkZqais7OTrz55pvQarWztvXw8IBWq4VerxeUpFeuXInMzExERUXR5jCRmUzUnPMixljkHK9fn/JlCQCHnIFijGHbtm04ffq0ybZXrlzBoUOH7BAVIYTMrb29HQUFBbh9+7ZZU8tzmWkELJfLodFoUFxcPGffwMBA9PX1mRz0AOOJPyUlBRkZGbQ5zIbEXqN+DcDZ2V5kjL0O4HUANqkUlpycjK+//nrOT4kAUFNTg9zcXAQGBooeAyGEmDI4OIizZ8+iqalpcje1uaRS6ax9ZxoBz7QBjTEGhUIBDw8PGAwGDA4Ooqenx+Szly5diqysLCQkJNDmMDsQLVEzxjZjPFGvn60N5/woxqfGkZGRIXrZG4lEgk2bNuH8+fMm2165cgXPPPOM2CEQQsg0nHN0d3dDrVajtrYWAwMDVr3fokWLMDAwYFGCl0qlkMvlkEgk0Ol00Ol0GBsbm3Y+ejaMMSQkJCArK4s2h9mZKImaMZYE4G0AuzjntllcESgrKwv5+fkmN2BUVVUhNzeXKuEQQkTFOcfNmzfR0NCAtrY2DA8Pz9leoVBAIpFAo9HM2W5iOtuaRG8wGMxO8O7u7li7di3S0tJoc5iDWJ2oGWPhAE4C+A7nfO7rrOxAKpUiJycHhYWFJtteuXIFBw4csENUhJD5bmhoCCqVCteuXZtzhMoYm1zP7evrM7lUN8Hc2tuMscnRs9FoFPycCYsXL8aWLVsQHR1Nm8McTMjxrD8B2AQgkDHWAuCfAMgBgHP+JoCfAAgA8NtvpkL0nPMMWwUsRE5ODq5cuWJyY0ZlZSU2btwIPz8/O0VGCJlPjEYjmpqaoFKp0NjYOGcylcvlMBqNMBgMgtaBzcEYw9KlSyc3gvX29kKj0ZidnAFg+fLl2LNnD20OcyLMETekAONr1OXl5TZ7/6+//ho3btww2S4tLQ379u2zWRyEkPmnt7cXKpUKlZWVGBoagkwms8l5Z3vx9PRETEwMdu7cCYVC4ehw5jXGWIW5g9l5UZlsJps2bUJpaanJ6SKVSoWNGzdi0aJFdoqMEOKKtFotamtroVKp8ODBAwBAQEAARkZGXC5JSyQSrFq1Cnl5eTSj6ALmbaJWKpVITExEVVXVnO0457h27Rp2795tp8gIIa6Cc47W1laoVCrU1NRAq9XCx8cHS5cuRUdHh80KkwjBGIOHhwc454LOPANAaGgo8vLy6MpfFzNvEzUA7Nixw2SiBoCbN29i/fr18PHxsUNUhBBnNzw8jKqqKqhUKnR3d0MmkyEkJARDQ0Po6+vD48ePHRJXUFAQfH190dnZiYGBAZM7ygHA19cXGzduRFJSEqRSqR2iJGKb14l6Yt2loaFhznYGgwHXr19HXl6enSIjhDgbo9GI5uZmqFQqNDQ0wGg0YsmSJVixYgVaWlomp7sB83dgW8vT0xNjY2Po6upCV1eXyfZyuRxJSUnYtGkTvLy87BAhsaV5naiB8SswTSVqAKioqMD69evph5qQBaa/v39yY9jg4CDc3d2xcuVK9Pf3o7Ozc7LdxDEnaygUCixZsgRKpXJyd7YQQkbOwHjFx82bNyMiIoIKkswj8z5RL1q0CJGRkbh3796c7fR6Pa5fv44dO3bYJzBCiMPodDrU1dVBpVLh3r17YIwhLCwM3t7e6OzsRFNTE4Dx5CyVSqHT6SxO0owxbNy4EXq9Hnfu3MHDhw8teg9g5pG8j48PMjIysGbNGirnOU/N+0QNAHv37sWvf/1rk+3Ky8uRk5ND1XcImafa29tx8+ZN1NTUQKPRwNvbG6Ghoejt7Z0xgXLOZ6yRLZSfnx8GBwdRWFg4uflLKKlUCqPRCM75UwlaKpUiOjoa27Zto+qKC8CCSNQBAQEICQlBW1vbnO10Oh2Ki4uxbds2O0VGCLG10dFRVFVVQa1Wo6OjA1KpFL6+vuCc4/Hjx09tDAsICIBEIkF3d7fVa9EKhQJpaWkYGBhAU1OTySnsJy/aePL5QUFByMnJQWJiIk1tLyALIlEDwJ49e/DWW2+ZbFdWVoacnBy4u7vbISpCiC1wznH37l2oVCrU1tbCaDRCqVROJsKpa8MymQxeXl7QaDTQaDQWH7maWMNWKpVITk6Gm5sb6uvrUVZWJvg9PDw8MDIyMq0mt1KpnNwYRr+XFqYFk6hDQkIQGBhosnSfVqtFSUkJNm/ebKfICCFiGRgYgFqtxs2bNzE4ODht1Dn10oupa756vR6PHj2y+JkKhQJarRYymQzx8fHo6uoSVBVxQnBwMLRaLXp7e6eN7sPDw7F161abXAlMXMuCSdTA+Kj6vffeM9mupKQEa9euhVKptENUhBBr6PV6NDQ0oLS0dNoRKmD2Y1RiHK+aGP3qdDoEBQXh0aNHUKlUgvoyxhAREYGuri60t7dPft/LywvZ2dlYu3YtnXkmkxZUoo6MjISPjw8GBwfnbKfVanHjxg3k5ubaKTJCiLk6Oztx+fJlNDU1WX1syhwT101OVAPjnAs62zzBw8MDo6OjkydRJBIJVq5ciZ07d9JFGGRGCypRA+Pnqj/++GOT7YqLi7FmzRq4ubnZISpCiBD9/f04d+4cmpubHVZf29rR+ESC9/HxwZYtW5CUlEQbw8icFlyijomJmZyymsvY2BhKS0uxYcMGO0VGCJlJWVkZrly54rCyneaaGHHPRCKRIDo6Gnv27KHiSkSwBZeoGWPYvn07Tp06ZbLtlStXkJ2dTde+EeIAHR0deP/99zE6Oiq4j4eHB5RKJR4/fmzV+WcAcHd3h0ajMXsEPVN7b29v7NixA/Hx8TR6JmZbcIkaAJKTk/H1119jbGxsznY6nQ6lpaVYv369nSIjhLS1teHUqVOC1n2lUilCQkLAGENnZydGRkYE3yQ1m4ld3OZ8QJhNZGQkDh06RKNnYpUFmagZY8jNzcX58+dNti0qKkJ2djaV5iPEhjjnuHPnDs6ePWvyHLNEIkFAQAAYY+jp6bGoJOdctFqtVf09PT2RnJyMzZs3QyZbkL9iicgW7E9RdnY2Ll++bHJ6TKfToaysDOvWrbNTZIQsHAaDAVVVVbh06ZLJql0KhWKypGd3d7edIhQmNDQU2dnZiIuLo2NVRHQLNlFLJBKsW7cOhYWFJttevnwZmZmZNKomRCQjIyMoLi5GSUnJrLu3n7ytSuhId2Lq2pYUCgWioqKQnZ2NsLAwWncmNrVgEzUArF+/HleuXDF5BlOv16O4uBgbN260U2SEzE8dHR24cOEC7ty5M2sbmUwGvV4v6Gz01B3WE/9uqyS9aNEiJCQkIDU1lc47E7ta0IlaJpMhIyMDpaWlJtsWFRVh3bp1tOZEiJk451CpVLhy5YqgUp2mzke7ublNbgSdusNajGpjTwoODkZycjLi4uLg7e0t+vuThWVqGVtzLPiss3XrVpSVlZn8S24wGHD58mVs377dTpER4tqGhoZw4cIF1NbWilqcZCJJe3h4QKfTWX0MayqpVIrly5cjKSkJUVFRVEaYiKa9vR1Hjx61qO+CT9QKhQJJSUmorKw02fb69evYuHEjVSsjZBacc9TX16OgoGDyeJW167dyuRx6vX7yw/TEGrS1x7AmuLm5YfXq1YiPj8fy5ctp1ozYhJCKmLOhn0gAO3fuFJSoAeDzzz/HkSNHbBwRIa7l8ePHKCoqQlVV1eQa8cRas7VT0k+OmMVYg/b09ERSUhJiY2MRGhoKiURi9XsSMpu6ujoMDAxY3J8SNcYrEEVHR6OxsdFk24aGBvT29tJmErLgGQwGNDQ04OrVq9NugJpI0OZMd89VdnOu18zh5uaGrKwsJCQkYPHixbRTm9gF5xx//vOfrXoPStTf2Lt3L37xi18IanvixAn81V/9lY0jIsQ5dXV1oaysDJWVlZOj3anJ1JL16LkSsTVJmjGGyMhI7Nmzhz5cE4coKiqyehaIEvU3vL29ERERgfv375ts29XVhYaGBsTExNghMkIcT6PRoKamBmVlZTOW9rQkmYo1Up6JVCrFmjVrkJubS/UPiMOMjY0JqtVhCiXqKfbt24df//rXgtqePn0af//3f0/TZ2Te4pzj3r17qKioQF1dneh3PtsiScvlcuTk5NAVtcQpnD59WpSfc5OJmjH2DoC9ALo45wkzvM4A/ArAbgAjAL7LOb9pdWQOEBAQgODg4GnrbbMZGRlBQUEBNm/ebIfICLGf/v5+FBcXo7q6+qlzn7YcBVtDLpdjzZo1WLt2Ldzd3R0dDiHo7e1FbW2tKO8lZER9DMCvAbw/y+u7AER98082gN99878uad++fYLPul29ehXZ2dnw8PCwcVSE2NbIyAiuXLmCmpoaDA0NzdrO2ZK0TCZDZmYmcnJy4Onp6ehwCJn0ySefiPZeJhM157yIMRY5R5MDAN7n43+DSxhjvoyxYM656WGpEwoODoa/vz/6+vpMtjUajTh16hReeOEFO0RGiPX0ej2qqqrQ2dmJjo4O9Pb2mrwMwxkxxpCRkYENGzZQxTDidJqamtDZ2fnU9y1dKhVjjToUwNR75lq++Z5LJmoA2L9/P44dOyaobWNjI1pbWxEaGmrboAix0oMHD/DHP/7R5hdWCCV0Gv3JyzmSkpKwefNm+Pr62jI8QixiMBjw+eefz/iapTNSYiTqmT4izBgNY+x1AK8DQHh4uAiPto2IiAi4u7sLvjj+5MmTeOONN2hjGXFKBoMBn332Gerq6hwdyjRCf2lNJGk/Pz+89NJL8Pf3t2VYhFilpKRkxtxhzf4OMRJ1C4BlU74OA9A2U0PO+VEARwEgIyPDuRa7nnDgwAEcP35cUNu+vj7cvHkT6enpNo6KEPPU1tbi888/F7XWtj1JJBIsXrwYmzdvpuOQxOkNDQ0hPz9/xtes2d8hRqI+DeANxthxjG8iG3DV9empYmJiIJVKYTAYBLU/d+4c4uLiaMcpcQq9vb04fvw4enp6HB2KIEqlEpxzjI2NQSqVIiYmBsnJyVi1ahWV9yQu49y5c6IfYwSEHc/6E4BNAAIZYy0A/gmAHAA4528COIPxo1m3MX4861XRo3SQHTt24OzZs4La6nQ6XLhwAfv377dxVITMbnR0FH/+858FlcN1FKVSifDwcHh6eqK3txcPHz6ERqNBSEgIUlJSkJCQQB94ictpbW1FTU3NU98X40ijkF3fc25p/ma394+sisJJZWZmCk7UAKBSqZCRkYGQkBAbRkXI07RaLc6fP4+KigpHhzKjgIAA5OTkwMfHBw0NDaipqcHo6Ci8vLywdu1apKSkYPHixY4OkxCLcM5x+vTpWV+zFlUmmwNjDGlpabh5U3j9li+//BJ/+Zd/SRvLiF1otVqcOnVKtMIKYlIoFMjNzUVcXBxu3bqF4uJidHd3QyqVYvXq1UhJScGKFStoapu4vKqqqhlL6z55YsFSlKhN2Lt3r1mJur29HSqVCmlpaTaMiix0LS0t+OKLL2b85eBogYGByMvLw+joKCorK3Hx4kVwzhEWFoY9e/YgISEBSqXS0WESIoqxsbFZZ17FWq+mRG0CYwxhYWFoaWkR3OfixYuIjY2ldTYiqtHRURQUFECtVjvNWegJjDFER0cjJSUFt2/fxqeffgqNRgMfHx/k5OQgOTkZgYGBjg6TENEVFRVhbGzMps+gRC3Aiy++iJ///OeC24+OjiI/Px979uyxYVRkITAajaisrJycNnY2crkcqamp8PDwQHV1NT7++GPIZDLExsYiOTkZy5cvp6ltMv6iTO4AACAASURBVG/19vaiuLhYcHtL/y5QohZAqVQiLi7OrHXA8vJypKWlITg42IaRkfmqvb0d+fn5uHPnjk2Oe1jL29sb0dHR6O/vR2lpKYDxIkbr1q1DfHw83VxFFoSzZ8+atVns0KFD+MlPfmL2cyhRC/Tcc8/hZz/7meDCEYwxnDlzBt/73vdoYxkRZHR0FIWFhaiqqhJcFc/eAgICEBgYOHn95aJFi7Bx40YkJydTxTCyoDQ1NaG5uVlwe39/f8TFxVn0LErUAjHG8N3vfhdvv/22oPacc7S0tECtViM1NdXG0RFXxTmHWq3G9evXnbo4SWBgIPR6PXp7ezE4OIi4uDgkJycjMjKSPoiSBcdgMOCrr74yq893vvMdi59HidoMoaGhiI6ONquYxIULF7B69WraWEamaWtrw8WLF3Hv3j2nuzpyglQqhaenJwYHB9HT04OIiAjk5uYiNjaWprbJglZSUoKBgQHB7dPS0qy6RIYStZkOHTqEf/3XfxVcWnR0dBSXL1/G7t27bRwZcXb9/f34+uuv0dzc/NTPjxjVi+ZizvvLZDJwzmEwGCCRSJCbm4vk5GT4+fnZLD5CXMXQ0BAKCgoEt5fJZFb//qdEbSaFQoGDBw/is88+E9ynrKwMqamptLFsgdFoNLh37x7Kyspw//79p5KzTCaDr68vent7bT6q5pxDIpGAcz7rsyaKM0gkEsTFxSElJQXh4eE0tU3IFBcuXDDrkpuDBw9CKpVa9UxK1BaIj49HSUkJWltbBbWnjWULg0ajwf3793Hv3j00NTWht7f3qTaMMcTExEAmk6G2ttau69Kmdo9HREQgOTkZsbGxUCgUdoqKENfR2tqKqqoqwe0DAwMRHx9v9XMpUVuAMYZDhw7hP/7jPwS1n9hYVllZiZSUFBtHR+zFaDTizp07aG5uxv3799HePvulcYsXL8aSJUvQ09OD+vp6O0Y5N09PT2RlZSEpKcmqNTRC5jvOOb788kuz+nz7298W5dmUqC3k5+eH3NxcFBYWCu5z/vx5rF69msonujiNRgO1Wo3S0lL09/eDMTZjIQOpVAp3d3cMDw+ju7vbqQqWKBQKHDhwALGxsTTLQ4gAlZWV6OjoENw+JSVFtH0dlKitsH79eqhUKgwODgpqP1GxjDaWuabe3l7cuHEDlZWV0Gq1kx+4JjZePclgMGBoaMjeYU4jkUgglUqh0+kAjJ+DTkpKwrp16yCT0V9/QoQYGxvD+fPnBbeXyWTYu3evaM+nv6lWkMlkOHjwIN5//33BfcrKypCWloalS5faMDIiFs45mpubUVJSgubm5mmjT41G48DInqZUKifraff392N4eBhGoxFhYWFYvXo1Vq9eTTu3CbFAYWGhWUWI9u/fb/UGsqkoUVtp+fLliI+Px61btwT3OXPmDF599VWacnRiWq0WarUa165dmzZjIpVKYTQap+2clsvlk6Nqe52Jlkgk8Pf3R0xMDJRKJXp6etDc3IyWlhZIJBKsWLECq1evRkxMDLy8vOwSEyHzUW9vL0pKSgS39/f3R2JioqgxUKIWQV5eHhobGyenF015+PAhbSxzUn19fcjPz0d9ff3kdLZUKkVERAS6u7vx+PHjp/oI/e9uLU9Pz8la2m1tbairq0N5eTnGxsYgl8sRFRWF1atXIyoqivZBECISc+t5v/jii6LHQIlaBF5eXti2bdusd5LOZKJiGf1CdTyj0YiKigoUFxejv78fwPjO/oiICISEhKC5uRl37txxSGwymQzbtm1DYmIiGhsbUV9fj8uXL0Ov18Pd3R2xsbFYvXo1VqxYAblc7pAYCZmvGhsbzarnnZSUZJOa95SoRZKRkQGVSiV4V+DIyAguX76MXbt22TgyMhPOOVpbW1FQUIC7d+9OnjH29/dHamoq9Ho9SktLcf/+fYfE5+Pjgy1btmBsbAz19fU4d+4cOOfw9vZGamoqYmNjERERQVdIEmIjBoPBrMGXVCrF/v37bRILJWqRSCQS7Nu3D2+99ZbgPqWlpUhNTaWNZXbU3d2NioqKaTdUKRQKxMXFTRayyc/Pd1j9bV9fXyxfvhxdXV3485//DGB8p/a6desQGxuLkJAQ2ttAiB2UlJTg0aNHgtvv2bNH1A1kU1GiFlFISAgyMjJQXl4uuA9tLLO9vr4+1NTUQKVSTfuLt2TJEuTk5KCzsxNqtRpqtdphMSqVSri5ueHRo0dQqVQIDg7G5s2bERsbi8WLFzssLkIWosePH5tVz9vPz8+mtyRSohbZ1q1bUVtbi5GREUHtHz58iKqqKiQnJ9s4soVlYGAAt27dQk1NzbSKYTKZDAkJCQgMDER1dTVOnjzpwCj/67KMsbExLFmyBGvXrsXq1auxaNEih8ZFyEJ26dIls+p5i1WBbDaUqEWmVCqxc+dOfP7554L7nDt3bvKYDbHc48ePUVtbi1u3buHhw4cAMLmG6+vri4iICAwMDKCystIprpaUSCRYuXIlYmNjER0dDU9PT0eHRMiCN1HuWai4uLjJ+gW2QonaBhITE6FSqXDv3j1B7UdHR1FQUIC8vDzbBjYPjYyMTCbniT9vpVI5OVINDAyE0WhET0+PWetNtrRs2TJkZWUhKiqK7nUmxIlwzvHVV18Jbi+RSPDMM8/YMKJxlKhtgDGG3bt347e//a3gPjdu3EBqaiqWLFliw8jmB41Gg7q6Oty6dQt37twB5xxeXl7w8fHB4OAgtFot5HI5tFoturq6HB0ugPGfiYyMDGzfvp2OURHipNRqtVn1vPPy8uxSipcStY0sXrwY69evx9WrVwX3+fLLL+kqzFmMjY2hsbERNTU1aG5uhsFggI+PD4KDg9Hd3T2tprbRaIRWq3VgtNPl5ORg69at9N+VECem0Whw7tw5we19fHyQmZlpw4j+CyVqG9q4cSOqq6sxMDAgqH1LSwuqq6uRlJRk48hcg06nQ1NTE27duoXGxkbo9Xp4eHggICAAg4ODk//Yg1wuh16vN2ttOzY2FocOHbLZkQ1CiHgKCwsxNjYmuL0tKpDNhhK1DcnlcuzevRt/+tOfBPc5e/YsoqOjF+TGMs45BgcHUVtbi4aGBrS0tMBgMEwbiY6MjAjeUW8OiUQyWfRkJuaUCQ0KCsJLL70Eb29vMUIjhNhYT0+PWfW8o6OjERQUZMOIpqNEbWPR0dGIiYlBQ0ODoPYajWbBVSwbHR2FWq1GUVHRjDdS2WOH9lxJWig3NzccOXIEkZGR1gdECLELzjm++OILwe0lEgmee+45G0b0NEGJmjGWB+BXAKQA3uac/98nXg8H8B4A32/a/C/O+RmRY3VZu3btQnNzs+BzeaWlpUhLS5v3G8taW1tRWlqKmpoaURKlozDGsHnzZuTk5FBJT0JcTFNTEx48eCC4/bZt2+y+IdRkomaMSQH8BsB2AC0AyhhjpznntVOa/T8ATnDOf8cYiwNwBkCkDeJ1SYsWLcKmTZtw8eJFwX2++OILvPbaa/NuA5JOp0NNTQ2Ki4vR3d3t6HAsIpPJsGTJEoSHhyMwMBDR0dF0lSQhLkiv10+W6hXC09MTa9eutWFEMxMyos4CcJtzfgcAGGPHARwAMDVRcwA+3/z7IgBtYgY5H6xZswZqtRo9PT2C2re2ts6rjWUTa0CVlZVmVfxxFosXL0ZWVhbi4uLg4eHh6HAIISK4fv36ZM1/Iey5gWwqIYk6FMDDKV+3AMh+os1PAZxnjP03AJ4Ats30Royx1wG8DgDh4eHmxurSpFIp9u7di2PHjgnuc+bMGcTExLhsUQyDwYC6ujpcvXoVnZ2djg5nGsYYZDLZnJvEgoKCkJGRgfj4eErOhMwz5tbzXr58OYKDg20X0ByEJOqZ5l6f3N3zAoBjnPN/Z4ytBfABYyyBcz5t4ZFzfhTAUQDIyMhwfA1HO4uIiEBSUhKqqqoEtR8bG0N+fr7LbSwbGBjAtWvXUFlZ6VTnmYHxD0wGgwGc8xmTdEBAANLT05GUlEQlPQmZx06dOiV4oypjDEeOHLFxRLMTkqhbACyb8nUYnp7afg1AHgBwzosZY0oAgQCcoyyUE9mxYwfq6+sFJ7DS0lKkp6fb9SiAJTjnaGhoQGFhoVmVfezNYDA89T1fX1+kpqYiPT2dkjMhC8DDhw/R3NwsuP3mzZuhUChsGNHchCTqMgBRjLHlAFoBHAHw5FUhDwBsBXCMMRYLQAnANXcK2Zinpye2b99uVj3Zzz//HK+//rpTbiwbHh5GYWGhXUfPCoUCOp3OqmNbXl5eSEpKwpo1a+i8MyELCOfcrNoW7u7u2LBhgw0jMs1kouac6xljbwA4h/GjV+9wzm8xxv4ZQDnn/DSA/w7gLcbY32F8Wvy73BmuJ3JS6enpqKioEDzy7OjocKqNZZxzNDY2orCwcNoVkrYWGhqKvr4+szZ/TKVUKhEXF4eNGzfSNZKELFA3btww63fICy+8YMNohBF0jvqbM9FnnvjeT6b8ey2AHHFDm78YY9i/fz+OHj0quM+XX37p8I1lo6OjyM/PR1VV1ayj54lbq8Qik8mQmJiI9vZ2tLa2WtQ/KioKW7ZssflVdIQQ52ZuPe+wsDAsW7bMdEMbo8pkDhIcHIzMzEyUlZUJaq/T6XD+/Hns27fPxpE9raGhAZcuXZr13PPU5CxWkvby8sKaNWvQ1tYGlUplVl/GGEJDQ7Fjxw6n+EtGCHEOx48fN6u9o45jPYkStQNt3boV1dXVM5bNnMnNmzeRnZ1tl41lw8PDOHv2LOrq6kxWDZNKpaKcjZZIJFi6dCk2bNiA+vp6swrEAOOFZTZv3oykpCSnXM8nhDhOR0cH7t+/L7j9+vXrnebOBUrUDuTm5oY9e/bgs88+E9znk08+wQ9/+EObJCLOOUpKSnD16lWTF1/I5XIYDAYYjUarkrRSqURsbCzS0tIwNDSEwsJCfPzxx4L7KxQKJCcn0z3PhJBZcc7xzjvvCG6vUCiwdetWG0ZkHkrUDhYfH4/S0lI8fPjQdGOMV/gSc2PZxLGqS5cumayaxhiDm5sbNBqNWbdJzSQsLAyZmZmT/3/Mnd5etmwZ8vLyEBISYlUchJD5r6CgwKzfWc8//7wNozEfJWoHY4zhwIED+M1vfiN4fff06dNWbSwzGAy4c+cOrl69OmMxerlcPu2HWqlUQqvVwmg0Cp6mn83EqLelpQUtLS1m9VUoFNiwYQOysrIceqaREOI6dDodioqKBLdfsmQJVqxYYcOIzEeJ2gkEBARg7dq1uH79uqD2BoMBX331FZ599lnBz9Bqtbh9+zZUKhWam5uf+lDg5uYGvV4Pg8EAnU4HiUQCzjk451Yn56ksGYkHBARg7969iIiIoLVnQohZfve735nV/uWXX7ZRJJajRO0kNm3aBLVabXJteEJ1dTXWr18/58aykZERNDQ04NatW7hz585TyVkqlQIYT/xjY2PTXnOGayfDwsLwzDPPwN/f39GhEEJc0O3bt9Hf3y+4fUZGhlPW9adE7STkcjkOHjyIjz76SHCfDz/8EH/3d3/31ChzZGQEBQUFJo9+zVRO0xkEBgbixRdfhK+vr6NDIYS4KM65Wb9P5XI5du/ebcOILEeJ2olERUVh+fLluHv3rqD2jx8/Rnl5OTIzMwGMJ96SkhLk5+cLGhFLJBKnGDkD4+vgiYmJyMrKosIkhBCr/fGPfzSrrsOhQ4ecdmmNErWTOXjwIH75y18K/gE7c+YMEhMTcePGDVy9etWso1KOTtJyuRwJCQnYsmULvLy8HBoLIWT+6OzsNOvSjYCAAMTExNgwIutQonYyPj4+2LhxIwoLCwX3+fnPf26yzWylPcUu+SlEaGgocnNzsWrVKqf9BEsIcU1GoxFvv/22WX1eeeUVG0UjDkrUTmjjxo0oLy/H8PCwaO85WzK2V5L29PRESkoKcnNzqTAJIcRmzp49a9bMYlJSktPfoEeJ2glJJBI899xzeO+99xwdilUkEgkiIyOxdetWKkxCCLG57u5ulJeXC24vlUpx8OBBG0YkDkrUTioyMhJhYWFmFwWZja03jk2dQvf19UV2djaysrIgkUhs9kxCCJlgNBpx7Ngxs/rs37/fJZbfKFE7mdHRUdTU1KCystKiax2nkkqlk0ewbL1xTCqVIjo6Gtu2bYOfn59Nn0UIIU+6fPmy4DoUwPh+ILFKMdsaJWonYDQa0dzcjMrKysnbqsQokSn2OWmZTPbU2k9QUBDWr1+PhIQEl/hkSgiZfzo7O3H16lWz+rz66qs2ikZ8lKgdqKurC5WVlaiqqsLQ0BDkcjkUCgU0Gg0YY9NGxI4ikUgglUqh0+kmk7RCoUBiYiI2b94MT09Ph8ZHCFnYDAaDWYVNACA6OtqlCipRorazkZGRyanttrY2AICfn9/kRRiLFi2CRCIxawrHFmQyGYxG4+Q/ABASEoItW7Zg5cqVDo2NEEImFBUVYXBwUHB7xpjT3Y5lCiVqOzAYDGhuboZarUZDQwOMRiN8fX2xePFidHd3Y2BgAKtWrcLAwAA6OzsdGuvEprOpo+fMzExs2LDB4tu6CCHEFtrb2826GQsA8vLyXG6TKyVqG+rs7IRarUZ1dTWGh4fh7u6OiIgIPHr0CP39/fDw8MC6devQ2dmJxsZGR4cL4L82nQUFBWHv3r1YtmyZgyMihJCnGQwGfPLJJ2b18fDwQFZWlo0ish1K1CIbGRlBdXU1Kisr0d7eDolEguXLl0Mul+Pu3bu4e/cufHx8sHz5cnR2dgq+2tIepFIpEhMTsWvXLrrvmRDi1AoKCsy6GQtw/gpks6FELQKDwYCmpiZUVlaisbFxcmp72bJlGBgYeKrm7ODgoFlrKrbm6emJ7du3IykpiXZuE0KcXkdHh9m7vCMiIua8FtiZUaK2kF6vR2NjI1QqFe7fvw+dTjetqMijR4/w6NEjB0c5t5CQEOzevRuhoaGODoUQQgQxGo04efKk2f1eeuklG0RjH5SoBdDr9WhoaEBbWxs6OjrQ0dHx1K5sqVTqoOjMFxsbi23btsHf39/RoRBCiFmKi4vR3d1tVp9NmzZBJnPddOe6kdtBT08PKioqoFarodFopr32ZElOR593FiI1NRWbN292+gL0hBAyk/7+fly8eNGsPm5ubsjNzbVRRPZBifoJer0edXV1qKiowP3798EYm/GTmKPvcraERqOhJE0IcUmcc4suKnrxxRdtEI19UaL+Rnd3N27evInKykqMjo5i0aJFCAoKQldXF3Q6naPDm2aiOIq56urqoFKpkJqaaoOoCCHEdvLz8zEwMGBWn6VLl86LI6YLOlHr9XrU1taioqICDx48gEQiQUxMDACgvr7eJnc1T+yqNve9/f39MTw8jLGxMas+OJw+fRqhoaEuu/uRELLw9PT0mL3LGwC++93vih+MAwhK1IyxPAC/AiAF8Dbn/P/O0OZbAH4KgAOo5Jx/W8Q4RdXd3Y2KigpUVlZCo9HA398fW7ZswePHj1FRUSHqtPbU6x8nCE3SE/W+9Xo9+vr6RIvp7bffxo9//GMolUrR3pMQQmzBaDTizTffNLtfdnb2vKmmaDJRM8akAH4DYDuAFgBljLHTnPPaKW2iAPxvADmc837GmNMN13Q63eTa88ToOTY2FomJiWhqasLly5dFHUHL5XJ4eHhgYGAASqUSnp6e6O3tNesZnPOnbqsCxhP4ypUrkZOTA19fX/zqV78yKzadTodjx47h+9//Pp2bJoQ4taNHj5q9WVcmk2Hnzp02isj+hIyoswDc5pzfAQDG2HEABwDUTmnzlwB+wznvBwDOeZfYgVqqq6trcu15YvS8fft2hIaGoqioCMePHxf1eeHh4UhNTcWVK1fQ19cHT09PDA8PP7Vr3BJ+fn5Ys2YNkpKSpo2Go6OjzS5B2tnZibNnz2L37t1Wx0UIIbZw+fJli+4/OHz48LwahAhJ1KEAHk75ugVA9hNtogGAMXYN49PjP+Wcfy1KhBbQ6XSTa88PHz6EVCpFbGwsUlNTMTQ0hMLCQly4cEGUZ3l7e+Px48fw9fXF2rVrUVtbi1OnTk2+Pjw8bNX7y2QyJCYmIjMzE8HBwTO2OXz4MH72s5+ZPSNQVlaGyMhIxMXFWRUjIYSI7f79+2ZfuAGM7+eJjo62QUSOIyRRz/Sx5MmMIAMQBWATgDAAVxhjCZzzaaW5GGOvA3gdGB95iq2rqwsVFRWoqqqCRqNBQEAAtm/fjpiYGNy4cQPHjx8XZQc3YwwhISHo6urC8PAwgoKC8OjRI5w9e1aE/xfjQkNDkZWVhdjYWMjl8jnbymQybNmyBZcuXTL7OZ988gl+9KMfITAw0NJQCSFEVIODgzh27JhFfV999VVxg3ECQhJ1C4Cp+9vDALTN0KaEc64DcJcx1oDxxF02tRHn/CiAowCQkZEhyoLwbKPn9PR0SKVSXLhwQbTRMzCeQPv7+9Ha2jr5va4u82f6pVLpU+suSqUS6enpSEtLM7tq2Pr161FUVGTRB5E//OEP+Nu//dt5s/GCEOK6dDodjh49alHfxMREeHl5iRyR4wlJ1GUAohhjywG0AjgC4Mkd3X8G8AKAY4yxQIxPhd8RM9AnzTR63rFjBxITE1FZWYkTJ05gdHRUlGdJpVJ4e3tjYGBgWoI2F2MMjDEYjcZpSToqKgrp6emIioqy6p7UI0eO4IMPPjC7n0ajwQcffIDXXnttXq3rEEJcC+ccH330kcVLhgcPHhQ5IudgMlFzzvWMsTcAnMP4+vM7nPNbjLF/BlDOOT/9zWs7GGO1AAwA/gfnvFfsYHU6HW7duoWKigq0tLRAKpUiLi4OaWlpkMvl+Oqrr3DhwgXRzz8bDIbJCzYYY3BzczNrc9jE6JlzPhmbj48PMjMzkZycLFq1sBUrVsDHx8eim7laW1tx7tw55OXliRILIYSY6+LFi7h3755FfQ8ePGjVQMeZMVsU9RAiIyODl5eXC2rb2dk5OXoeGxtDYGAgUlJS4O/vD7Vajebm5mkj1JnOLpvryVreSqUSOp3OqpreEokE8fHxSEtLQ0REhE1Gr48ePTL7uNZUhw8fps1lhBC7q62txSeffGJRX3d3d/zDP/yDyBHZBmOsgnOeYU4fp61MptVqcevWLdy8eXNy9Lxq1Sr4+fmhu7sb+fn5sxYmsTRJT03ORqMRMpkMXl5eePToETQajcWJNSgoCBkZGUhMTLR5kRFfX1+EhoZaPEX/6aef4kc/+hECAgJEjowQQmbW2dlpcZIGgO9973siRuN8nC5RPzl69vHxQXh4OAYHB9HQ0CDKM6auFU+QSqWT31MoFOCcQ6fTTbtT2pwPAHK5HMnJyUhLS5v1WJWtvPzyy/iXf/kXi/pyzic3lykUCpEjI4SQ6UZGRvDuu+9a3D80NHTen1pxikQ9MXouLy9HW1sbGGPw9PSETqfD4OCgRWuuM5mYEp+6Vjxh6pS2Vqu1+BnLli1Deno64uLiTB6rshWFQoG4uDjU1taabjyD0dFRvPfee/iLv/gL2lxGCLEZg8GADz74AGNjYxa/x3yp5z0XhyVqzjlu376N8vJy3L59GwaDYdqFFaOjozNObVuz/jzR78n1Z2u5u7sjLS0NqampTjNlfOjQITx8+BCPHz+2qH9bWxu+/vpr7Nq1S+TICCFk3OnTp9HR0WFx/yNHjsx4DfF847DNZCEhIfz73/++4PZibBAT28qVK5GRkYGoqChIpVJHh/MUnU6Hf/u3f7OqyMuhQ4eQkJAgYlSEEAKUlpZaVSTqueeeQ3x8vIgR2ce82kz2JGdJ0lKpFGvXrkVmZiZ8fHwcHc6c5HI53njjDfzyl78E59yiDzsnT57E0qVL5/0aECHEfu7du2dVkt67d69LJmlLzc9DZyKRSCRwc3ODr68vgoODkZeXh3/8x3/E1q1bnT5JT/Dx8cF3vvMdAOMfdsw9ZzixucyadXtCCJnw6NEjfPTRRxb337p1K9LT00WMyPm5zIja1mQyGZYuXYrY2FjExsbC19d33mykWr58ObZv344LFy7AaDRCLpebNR2u0Wjw7rvv4vXXX583fyaEEPvTarV4//33LV6Oy8rKwvr160WOyvkt2EStUCgQGRmJ9PR0rFy50inXmMW0bt06tLa2ora2Fjqdzuxk3dHRga+++gp79+61YZSEkPmKc44TJ06gv7/fov4xMTELtnLigknUUqkUwcHBWLt2LWJiYuZ9Yp7Js88+i66uLvT09ECn0814MchcKioqEBYWhpSUFBtGSQiZjy5duoTm5maL+gYHB+Nb3/rWgp3Rm9eJ2tfXF1lZWUhKSoKnp6ejw3E4qVSKV155Bf/5n/8JrVZrUTnU06dPIyQkBEFBQTaIkBAyH9XU1ODatWsW9fXx8cGrr746b+t4CzHvEnVQUBC2b9+OlStXLthPX3Px8vLCd77zHfzhD38AYP6xN8453nnnHfz4xz+mymWEEJPa29tx8uRJi/rK5XL84Ac/cFjxKGcxLxK1XC7HmjVrsHbtWri7uzs6HKcXFhaGffv24YsvvgDn3Owp8LGxMbz11lv44Q9/SB+GCCGzGhoawgcffGDR8VrGGH74wx/S73S4YKJWKpWQyWQYGhqCl5cX1q1bh/T0dBrdmSktLQ2tra24efMmDAaD2cm6p6cHn3/+OZ599lkbRkkIcVV6vR5//OMfMTo6alH/1157Db6+viJH5ZpcJlEHBASAc46+vj74+flh7969SE5OXhDl42xl9+7d6OzsRGtr62QJV3M++VZXV2PZsmXIzMy0YZSEEFfDOccXX3xhcXnQZ599FqGhoSJH5bqcfnVeKpXCw8MDvb29kMlkePbZZ/HGG28gPT2dkrSVpFIpjhw5MrnRzpLpqTNnzlh8pSYhZH4qLi5GVVWVRX3XrVuHxMREkSNybU6dqCUSCQwGA/z9/fHCCy/gBz/4ARITExf07j+xeXl54YUXFbo02AAAHzNJREFUXpj8M7Xk2NqxY8csnt4ihMwvt2/fxoULFyzqGxkZie3bt4scketz6owXERGBV155Bd/73vcQHR1NG5dsJDQ0FPv27QMwfu2cuev9er0eb731ltPUYyeEOEZvby8+/vhji/p6eXnhlVdeETmi+cGhc8cTa6JPro1OfKoKCQlxYHQLS0pKCtrb21FaWgqtVmt25bL+/n58/PHHOHLkiA2jJIQ4K41Ggw8++AB6vd7svhKJBH/9139tg6jmB4cl6qmj44kkHRAQgGeffZYStIPs2LEDHR0dePDggUW1eBsaGnD16tUFWYuXkIXMaDTixIkTGBgYADCeeI1Go+D+dFZ6bg6b+uacTyZomUyGffv24Y033qAk7UBSqRTf+ta34O3tbfF7XLp0CXfv3hUxKkKIszt//vzk33tzk/SBAwewePFiW4U2Lzh8jToyMhJ/8zd/g7S0NEeHQgB4enpO21xmiQ8//BCDg4MiRkUIcVZqtRo3btyY/NqcJJ2cnEx3BwjgsETNGMOuXbvw8ssvw8vLy1FhkBkEBwfjwIEDFvc3Go146623LKolTghxHQ8fPsSpU6cs6hsYGGjV75mFxGGJevHixcjKyqKd3E4qKSkJa9assbj/ROlAQsj8NDAwgI8++siivnK5HK+++ir9/hfIYYmaipU4v+3bt2P58uUW979//77F5ykJIc5Lp9Phww8/hEajsaj/yy+/DA8PD5Gjmr8cvkZNnJdEIsFzzz2HRYsWWfzJ9/r166irqxM5MkKIo3DO8dlnn6Gnp8ei/jt27EBYWJjIUc1vlKjJnDw8PHDkyBFIJBKLk/WJEyfQ29srcmSEEEcoLCxEQ0ODRX2jo6OtWlJbqChRE5OWLl2KgwcPWlV57Pe//z3GxsZEjIoQYm91dXUoLCy0qK+Pjw8OHTpE69IWoERNBElISMC6dess7q/T6ajMKCEurLOzEydOnLCoL2MML730El1HbCFBiZoxlscYa2CM3WaM/a852j3HGOOMsQzxQiTOYuvWrVixYoXF/Xt7e/HJJ5+IGBEhxB6Gh4fx7rvvWtz/mWeeoaImVjCZqBljUgC/AbALQByAFxhjcTO08wbw1wBuPPkamR8mNpf5+flZXBClrq4O169fFzkyQoitGAwGvPvuuxYvXWVkZNC1lVYS8ts2C8BtzvkdzrkWwHEAM51S/z8A/hWAZfv1iUtwd3fH888/D6lUanGyvnDhAu7duyduYIQQ0XHOrdoMumTJEuzcuVPkqBYeIb9pQwE8nPJ1yzffm8QYSwWwjHP+5VxvxBh7nTFWzhgr7+7uNjtY4hyWLFmCAwcOwGg0Wrwx5L333pss4E8IcU75+flobGy0qK9CocCRI0eoZoYIhCTqmX4TT+4IYoxJAPx/AP67qTfinB/lnGdwzjNovcK1xcfHIycnx6rNYb/97W8tuqWLEGJ71dXVuHr1qsX9n3vuOfj6+ooY0cIlJFG3AFg25eswAG1TvvYGkACggDF2D8AaAKdpQ9n8t2XLFqxatcri/lqtFm+++SbtBCfEybS2tuLkyZMW99+wYQOioqJEjGhhE5KoywBEMcaWM8YUAI4AOD3xIud8gHMeyDmP5JxHAigBsJ9zXm6TiInTkEgkePbZZ63aXNbX12fxkQ9CiPiGhoas2uEdERGBTZs2iRcQMZ2oOed6AG8AOAegDsAJzvktxtg/M8b22zpA4tzc3d1x5MgRSKVSSKVSi96jvr7e4iIKhBDx6PV6vPnmmxbffOfh4YHDhw9bdU0ueZqgP03O+RnOeTTnfCXn/P/95ns/4ZyfnqHtJhpNLyxBQUF45plnYDAYLN5cVlBQQDXBCXEgzjnefvttDA8PW/wezz//PDw9PUWMigBUmYyIJDY2Fhs2bLBqvfnEiRNob28XMSpCiFAnTpxAZ2enxf137NiB8PBwESMiEyhRE9Fs2rTJ6g0kR48eRVdXl0gREUKEuHjxIurr6y3uHxMTQ5dt2BAlaiKa/7+9ew2K6kzzAP5/uhuaO8hFjQJKEFFHjRq8ojHqrEST8RKI4q0sx6pMkjXXnUrN7n7Yrf002WztZKoSU0bMTcWIZjJRK4nGEUWRGG4aFaNx1DFEUbwgjaB94d0PgIsI0tCn+xya/68qFeg+/ZynurD/fc55z/u2DC6Ljo7u9vVqAHj//fdx7tw5DTsjoo7k5+ejsLCw26+PiIjAggULuNiGFzGoSVNBQUH3Bpd58g9348aNKCsr07AzImpNKYVvvvkGBQUF3a4hIli6dCmCgoI07IzaYlCT5uLi4rBw4UKP74/euXMn9uzZo1FXRNTC4XAgLy8PR448fGmGzr5sz5s3D/369dOyNWoHg5q8YtiwYZg+fbrHdYqKirBp0yY0NjZq0BUR2Ww2fPzxx51ekzabzQ/9sj169GiMGTNG6/aoHQxq8prp06cjNTXV4zp///vf8d5776GhoUGDroh6r6qqKqxfvx6XLl166HaBgYEPvZc6OjoazzzzjNbtUQcY1OQ1IoKFCxciNjbW44EmN27cwDvvvAMu5kLUPadPn8aGDRtgs9keup3ZbIbdbu/weZPJhGXLliEgIEDrFqkDDGryKqvVisWLFyMgIMDjsLbb7Vi7dq1Ht5EQ9TZKKRw+fBifffZZpzOOiUin22RlZSE6OlrLFqkTDGryutjYWGRmZmq2+MbWrVtRUFDAxTyIOuFyubBz5058++23AODxv5m0tDQMHz5ci9aoCxjU5BNDhw7VdKL+/Px8bNu2DU6nU7OaRP6koaEBmzZtQnl5uVvbBwYGPjTI4+LiMGfOHK3aoy5gUJPPPPHEExg2bJhm9U6dOoWcnBzU1dVpVpPIH1y/fh3r16/HhQsXOt1WRBAYGPjQ69IWiwUrVqzgYhs64btOPiMiWLBgAWJjYzWreeXKFaxduxZVVVWa1STqyc6fP48PPvgANTU1nW5rMpmglHpoSANNi22Eh4dr1SJ1EYOafMpqtSI7OxsWi0Wzmg0NDVi/fj1X36Jer6ysDJ9++ikcDodb16PdmZ9g0qRJGDJkiBbtUTcxqMnnYmJisGjRIk1rNjY2Ii8vDwcOHOAgM+p1GhsbsWfPHuzcuROA54PGWvTr1w+zZ8/WpBZ1H4OadJGSkoIZM2Z0+Lwn61pv374dDoeju60R9Sh2ux15eXkoKirStK7FYsHKlSu52IYBMKhJN9OmTetw5jJPjggqKiqwYcMG1NbWdrsGUU9w69Yt5OTk4PTp026/xt2V7ZYtW4bg4ODutkYaYlCTbkQEzz77LKKiojSvfeXKFaxbtw6VlZWa1yYygl9++QXr1q1ze7a+liPjziY0AYD09HQMHjzYk/ZIQwxq0lVgYCBWrFjh1rf8rp6Cq6+vx0cffYQffvihu+0RGdLJkyfx4Ycfuj3/vdVqdfssVb9+/TBr1ixP2iONMahJd9HR0Vi8eHGn23XndHhjYyO++OIL7N27lytwUY+nlMKBAwewffv2B/6eO/oia7VacffuXbfqm81mrFq1itelDYZBTYaQkpKCjIwMr9UvLCzE1q1b3f7AIjIap9OJzz//HPv373/guZCQkA6/yHblb37FihWwWq3dbZG8hEFNhjFp0iS89NJL6Nu3r1fqnzlzBhs2bMDNmze9Up/IW27fvo2cnBycPHnygeeio6NRX1/vUf2IiAi8+OKLGDRokEd1yDtEr3tO09LSVElJiS77JmNTSqG4uBi7d+/W/HS1iCAoKAiLFi3iYBnqEaqrq/HRRx89cD1aRBAeHu7x3Q3jx4/H7NmzNZ2EiDomIqVKqbQuvYZBTUbV0NCA3NxczUdut8xXPGfOHKSldenfC5FPnTlzBlu3bn3gC2tQUBAcDodbI7g7YrVasWDBAk3n36fOMajJL5WXl2PXrl2aHl2LCJRSGD9+PDIyMty+t5TIVwoKCpCfn//A45GRkbh165ZHtePj45GVlYXIyEiP6lDXdSeoea6DDG/s2LFITU3Fpk2bcPnyZU1qtnxBLS4uxrVr1/Dcc89xcgcyhJbpcNubxCQ0NNTjkJ46dSpmzJjBlbB6EB5RU4+hlMKRI0ewe/fudp9vOUrujujoaGRnZyMuLs6TFok80tDQgJycHNy4ceOB50wmk0dnlYKDg5GVlYVHH33UkxbJQzz1Tb3CtWvXsHHjRk2nCDWbzbBYLMjMzERKSopmdYncVVVVhQ0bNsDpdGpeOykpCc8++yzCwsI0r01d052gduvch4g8JSKnReSsiPyhneffEJEKEflBRP4mIhzjT14TGxuLl19+GePHj9espsvlgtPpRG5uLoqKirgCF/nUsWPHsG7dOq+E9MyZM7FixQqGdA/WaVCLiBnAewDmABgBYImIjGizWTmANKXUaADbAfy31o0StWaxWDB37lwsWbJEs9tKWkbQ7tmzBzt27PDKhyZRWzt37sRf//pXzeuGhIRg1apVmDZtGmca6+HcOaKeAOCsUuqcUsoO4DMA81tvoJTKV0q13HH/HYB4bdskat/QoUPxyiuvaD5JytGjR/HJJ5+grq5O07pELRwOB9auXYuysjKPa7UdGDZkyBCsWbMGiYmJHtcm/bkT1AMB/Nzq98rmxzqyGsDXnjRF1BXh4eF44YUXMHHiRE3rVlZWYt26daiqqtK0LtHVq1fx1ltvub3yVWdaDzKbPXs2li5dyrsY/Ig7Qd3eOZN2L+CJyHIAaQDe7uD550WkRERKtPoDJQKaRnw/9dRTWLp0qaa3ndTV1SEnJwenTp3SrCb1bqWlpXj//fc9mqykPYGBgXj++ecxefJknur2M+58olUCSGj1ezyAS203EpFfA/h3APOUUu3OAq+U+kAplaaUSuNtMOQNKSkpeO211zQdOONyuZCXl4f8/HwOMqNuU0ohNzcXu3bt8rhW2y+jCQkJeOONN/DII494XJuMx52gLgaQIiJJIhIIIBvAjtYbiMhYAOvQFNJXtW+TyH3h4eF47bXXkJSUpGndgoICbNmyBQ6HQ9O65P9qa2vx9ttv46effvK4Vtv7qWfOnIlVq1Zx1Ss/1mlQK6WcANYA2A3gFIA8pdRJEfkvEZnXvNnbAMIAbBORoyKyo4NyRD5hNpuxYsUKpKena1r3p59+wrvvvqvpPdzk306cOIE//elPDyyq0V0tIS0iWL16NUd19wKc8IT83vHjx/HFF19oetpaRLB8+XLO8kQdamxsxLZt2/Djjz9qXjssLAxr1qzhUXQPxLm+idoxatQoREZGYvPmzbDb7ZrUVEph48aNGDlyJDIzMzWpSf7j5s2bWL9+vWZH0a0NHToUS5Ys0bwuGRePqKnXuHbtGjZt2uTxogZtmUwmLFq0CKmpqZrWpZ7HbrcjPz8f3333nVfqz507V9MZ+cj3ONc3USfq6uqwZcsWXLr0wI0LHrNarZg3bx6GDx/Oa4a9jMPhQHFxMfbv3++1wYYvvvii5hP7kO8xqIncYLfb8fnnn+PMmTMer0jUHqvVihkzZmDs2LEIDAzUtDYZi8PhQElJCQoKCnDnzp0Hng8ICPA4uAMDA/H6668jKCjIozpkDAxqIjc1Njbi66+/RklJCYKCgtr9kPWU2WzGhAkTMHHiRERGRmpen/TTEtCHDh1CfX39A89bLBZN1o7u27cvfve733HtaD/CoCbqAqUUDh8+jL179yIqKgo1NTVe29eIESMwZcoUDBz4sNl3yegcDgdKS0tx6NAh3L59u91t4uPjUVVV5fGiLmPGjMH8+fM735B6FI76JuoCEUF6ejoiIiLw5Zdfok+fPqipqfHK7GMVFRWoqKjAwIEDMWXKFAwbNoxHST1IS0AXFhairq4OZrP5gW0sFgucTicqKys93t/TTz+NtLQufZaTH+MRNRGACxcuYOvWrTCbzRCRbq2a1ZXr3WFhYZgyZQrGjh3La48G1jagrVYr7t5td4ZkzaxevRrx8VyA0F/x1DeRB6qrq7F582bU19cjISEB586d61aduLg4XL9+3a3QNpvNGDduHCZPnow+ffp0a3+kPafTee8Ud11dHUJDQzs81a2VgIAAvPrqqwgNDfXqfkhfDGoiD9lsNuTm5uLKlSsYN24cSktLu1XHZDJh0KBBOH/+vNuvSUpKwvTp05GYmMjbu3TSEtCFhYWw2WwIDw9HXV2dZpdDWk6PtxUVFYWXX36Zl0N6AQY1kQbsdju2bduGs2fPYty4cSgrK+t2LRFBSEhIl47GIiIiMGPGDIwaNarda6GkPafTibKyMhw6dAg2mw2RkZGw2Wya3rrXUUgnJydj+fLlmu2HjI1BTaSRxsZG7Nq1C+Xl5Rg+fLgu61G33N41depUhISE+Hz/vUHbgI6IiIDNZtN8QKGItFtz6tSpmDVrlqb7ImNjUBNpSCmFgwcPIj8/H/Hx8ZqM5u2uhIQE/OY3vwHXcdeG0+lEeXk5Dh48CJvNhtDQUNTX1/t0vfHnnnsOI0aM8Nn+yBgY1ERecOzYMezYsQORkZG4efOmrr0EBwcjIyMDo0eP5nXsbmgJ6EOHDqG2trZLk910dFTcHS+88AL69eunSS3qWRjURF5y7tw55OXlwWQytbsikpYf4u4QEYwYMQLz589HQECAz/bbUzmdThw9ehQFBQWw2WwdXi/2hd///vcc2d2LMaiJvOjKlSvIzc1FfX29bh/y7YmIiEBWVhYSEhL0bsVwXC4XysvLceDAAdTV1fn8C1Vbb775JoKDg3XbP+mPQU3kZbW1tcjNzcXVq1d1/cBvT8s92RkZGb1+tLjL5UJZWRn279/f7lzcenjppZc4xoAY1ES+cPfuXWzdurVL90j7Wp8+fZCVlYUBAwbo3YpPuVwuHDx4EEVFRbDb7Xq3c092djbXKycADGoin3G5XNi5cyeOHTumdysPZTKZMGbMGGRkZPj1kpsNDQ3YsWMHTp8+bbgzHU8++SSmT5+udxtkEAxqIh9SSuH48eOorq5GTU0Nbty4gdraWty+fdtwYQE0rWs8ZswYjBo1ChbL/6/H07ZXo/0ONM25fefOHTQ0NKCurg63bt2CzWbDtWvXvD61pydSU1ORnZ2tdxtkIFw9i8iHRASjR49+4HGn04mff/4ZZWVlOHXqFFwulw7dPchut+P777/H999/r3crvUJMTAxDmjTBoCbSmMViQVJSEpKSkqCUwokTJ7Bnz55urcjlS63nmVZKGfKsgJFZrVZER0dj8ODBSEpKwpAhQ/RuifwEg5rIi0QEo0aNwsiRI1FRUYF9+/bhxo0berfVLk/mtRaRe/+11RL4/hb+kZGRGDx4MEaPHo2kpCROQENew6Am8gERwa9+9SsMHz4cx48fx/79+1FTU6N3Wx5ruS+5sxBuL8hbB3hPYLVakZKSgqFDhyI5OZnzr5PPMKiJfMhkMuGxxx7DyJEj782UVVtbC5PJpOlKTb7ibsj21KPp/v37IzU1FSkpKRgwYACPmkkXDGoiHZjNZjz++ON47LHHUFpaioMHD+L27du6Tm3ZEwUFBcFut2v2JUdEkJycjFGjRiE5OZlTfZIhMKiJdGSxWDBx4kSMHTsWxcXFKCwshNPp7NJiEb1R//79UVNTo8l7JCJ49NFHMW3aNCQkJNw3qI7ICBjURAYQGBiI9PR0pKWl4bvvvkNRUREAIDQ01ND3Cfua2WyGy+VCVVWVx7UiIyORmZnJOdLJ8BjURAZitVoxffp0TJgwAYcPH8aRI0cANC28UVtbq3N33uXOghla3JOenJyMrKwsBAUFeVyLyBfcCmoReQrAnwGYAeQopf7Y5nkrgE8BPA7gOoDFSqkL2rZK1HsEBwdj1qxZmDhxIgoLC1FcXAwRQXh4eJcDW+8Vo9zl7R4HDRqEBQsWICoqyqv7IdJap0EtImYA7wH4JwCVAIpFZIdSqqLVZqsB3FRKDRGRbABvAVjsjYaJepOwsDBkZGRg8uTJOHjwIMrKymAymRAUFOT2qlBKKZjNZpjNZkMtVOErUVFRyMzMRHx8vN6tEHWLO0fUEwCcVUqdAwAR+QzAfACtg3o+gP9s/nk7gHdFRFRP+BpP1ANERETg6aefxpQpU1BQUICjR4/CbDbDZDLB4XB0+nqXywWXywWz2YyQkBA4nU7Y7XbDTG/qDRaLBfPmzcPIkSN5WxX1aO4E9UAAP7f6vRLAxI62UUo5ReQWgBgA11pvJCLPA3geABITE7vZMlHv1adPH8yfPx/p6ek4cOAATpw4cW+AlTtcLhdsNpuXu9TfpEmTMGvWrPsWHyHqqdz5K27vq2jbI2V3toFS6gMAHwBNq2e5sW8iakdsbCwyMzMxdepU7N+/Hz/++KPuk6a0zDzmq8lNRARWq/W+W7QSExOxePFizhpGfsWdoK4E0Pr+hXgAlzrYplJELAAiARhzQmMiP9KvXz8sXrwYly5dQn5+Ps6ePavb4DFfzj6WmJiIixcv3gvp0NBQLF++HP379/fJ/ol8yZ2gLgaQIiJJAH4BkA1gaZttdgBYCaAIQBaAfbw+TeQ7AwYMwLJly3Dx4kXs27cP//jHP/RuySvCwsLgdDpx8eJFAE1Tss6dOxePP/64zp0ReU+nQd18zXkNgN1ouj3rQ6XUSRH5LwAlSqkdADYA2CgiZ9F0JM1FWIl0kJiYiJUrV+L8+fPYu3cvLl++7FE9EYHZbIZS6r7r4CaTCTExMYiLi0NcXBz69u2LuLg4REdHw2w2w+Fw4MaNG6isrMSVK1dQXV2Nmpoa2Gy2++p05eg/PDz8vuvrI0eOxMKFCzmTGPk90evANy0tTZWUlOiyb6LeQCmFs2fP4ptvvnFraU2z2YzGxsb7gjMgIAAxMTGIjY1FdHQ0YmNjERMTg8jIyHsB2XpEdWc/K6VQW1uL6upqXL16FVevXkV1dTWqq6s77Cs4OBh37ty511efPn3w29/+FmFhYW6+E0TGISKlSqm0rryGQyKJ/JSIICUlBUOGDEFFRQW++uqrh9573d7IcYfDgaqqKk2m7OyuhoYGAE1fJLKzszFkyBDdeiHSA4OayM+1Xgu7tLQUxcXFCAsLQ3h4OMLCwhAWFobAwMD7jqSN8rPT6UR1dTX69u2LmTNn8n5o6pUY1ES9hMlkwvjx4zF+/Hi9WyGiLuAoDCIiIgNjUBMRERkYg5qIiMjAGNREREQGxqAmIiIyMAY1ERGRgTGoiYiIDIxBTUREZGAMaiIiIgNjUBMRERkYg5qIiMjAGNREREQGxqAmIiIyMGm9rJxPdyxiA3Bal533HrEArundRC/A99n7+B57H99j30hVSoV35QV6LnN5WimVpuP+/Z6IlPA99j6+z97H99j7+B77hoiUdPU1PPVNRERkYAxqIiIiA9MzqD/Qcd+9Bd9j3+D77H18j72P77FvdPl91m0wGREREXWOp76JiIgMTJegFpGnROS0iJwVkT/o0YM/E5EEEckXkVMiclJEXtW7J38lImYRKReRXXr34q9EJEpEtovIj81/05P17snfiMjrzZ8VJ0Rki4gE6d1TTyciH4rIVRE50eqxaBH5VkR+av5/H3dq+TyoRcQM4D0AcwCMALBEREb4ug8/5wTwL0qp4QAmAfhnvsde8yqAU3o34ef+DOAbpdQwAI+B77emRGQggFcApCmlRgIwA8jWtyu/8DGAp9o89gcAf1NKpQD4W/PvndLjiHoCgLNKqXNKKTuAzwDM16EPv6WUuqyUKmv+2YamD7aB+nblf0QkHsDTAHL07sVfiUgEgCcAbAAApZRdKVWjb1d+yQIgWEQsAEIAXNK5nx5PKVUA4Eabh+cD+KT5508ALHCnlh5BPRDAz61+rwRDxGtEZDCAsQCO6NuJX3oHwJsAGvVuxI89CqAawEfNlxhyRCRU76b8iVLqFwD/A+AigMsAbiml9ujbld/qp5S6DDQdUAHo686L9AhqaecxDj33AhEJA/A5gNeUUrV69+NPROQZAFeVUqV69+LnLADGAXhfKTUWwG24ebqQ3NN8nXQ+gCQAAwCEishyfbui1vQI6koACa1+jwdPs2hORALQFNKblVJ/0bsfP5QOYJ6IXEDT5ZuZIrJJ35b8UiWASqVUyxmh7WgKbtLOrwGcV0pVK6UcAP4CYIrOPfmrKyLyCAA0//+qOy/SI6iLAaSISJKIBKJp0MIOHfrwWyIiaLqmd0op9b969+OPlFL/qpSKV0oNRtPf8D6lFI9CNKaUqgLws4ikNj80C0CFji35o4sAJolISPNnxyxwwJ637ACwsvnnlQC+dOdFPl+UQynlFJE1AHajaXThh0qpk77uw8+lA1gB4LiIHG1+7N+UUl/p2BNRd70MYHPzF/tzAFbp3I9fUUodEZHtAMrQdMdIOThLmcdEZAuAJwHEikglgP8A8EcAeSKyGk1fkJ5zqxZnJiMiIjIuzkxGRERkYAxqIiIiA2NQExERGRiDmoiIyMAY1ERERAbGoCbqRZpXonpJ7z6IyH0MaqLeJQoAg5qoB2FQE/UufwSQLCJHReRtvZshos5xwhOiXqR5NbVdzesOE1EPwCNqIiIiA2NQExERGRiDmqh3sQEI17sJInIfg5qoF1FKXQdQKCInOJiMqGfgYDIiIiID4xE1ERGRgTGoiYiIDIxBTUREZGAMaiIiIgNjUBMRERkYg5qIiMjAGNREREQGxqAmIiIysP8DI0UBAyYpfb0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean Square Error compared to true trajectory\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9370665441703778"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Plot to show the functions\n",
    "fig, ax = plt.subplots(figsize=(8,6)) \n",
    "t_dat.groupby('ID').plot(kind='line', x='t', y='outcome', ax=ax, c='grey', legend=False)\n",
    "plt.show()\n",
    "\n",
    "#Should add in what the best case loss is (mean square error)\n",
    "print(\"Mean Square Error compared to true trajectory\")\n",
    "((y - y_werr)**2).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next is the actual good stuff. Basically I create a module class to host my parameters, and do the subsequent manipulations to go from my time series data, coefficients, and the latent classes to the end predicted result. The example is borrowed pretty heavily from the pytorch [example documentation](https://pytorch.org/tutorials/beginner/examples_nn/two_layer_net_module.html).\n",
    "\n",
    "To do this in practice you need to choose the number of groups beforehand (and well as how to estimate the function over time). \n",
    "\n",
    "It takes quite a few iterations to converge in this example using the Adam optimizer. (For this simulation it did not converge to the correct estimates when using stochastic gradient descent.) This prints out the loss in a loop, and so you can see that here it took around 25k iterations to get where it needed to be! But we do see the results calm down to the error that we pretty much expected based on how I simulated the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99 638.4146728515625\n",
      "1099 458.09051513671875\n",
      "2099 338.63336181640625\n",
      "3099 262.8341064453125\n",
      "4099 215.86953735351562\n",
      "5099 184.10838317871094\n",
      "6099 158.9413299560547\n",
      "7099 136.25038146972656\n",
      "8099 114.6995849609375\n",
      "9099 92.67115020751953\n",
      "10099 71.51577758789062\n",
      "11099 53.634376525878906\n",
      "12099 40.28346252441406\n",
      "13099 30.74210548400879\n",
      "14099 23.511497497558594\n",
      "15099 17.468568801879883\n",
      "16099 12.363271713256836\n",
      "17099 8.303993225097656\n",
      "18099 5.275670051574707\n",
      "19099 3.1743240356445312\n",
      "20099 1.8652191162109375\n",
      "21099 1.1816942691802979\n",
      "22099 0.9176636934280396\n",
      "23099 0.8551831245422363\n",
      "24099 0.8466705083847046\n",
      "25099 0.8454619646072388\n",
      "26099 0.8452392220497131\n",
      "27099 0.8452112674713135\n",
      "28099 0.8452098965644836\n",
      "29099 0.8452096581459045\n"
     ]
    }
   ],
   "source": [
    "##############################################\n",
    "#Adapted heavily from\n",
    "#https://pytorch.org/tutorials/beginner/examples_nn/two_layer_net_module.html\n",
    "\n",
    "class LatGBTM(torch.nn.Module):\n",
    "    def __init__(self, groups, poly, obs, periods):\n",
    "        \"\"\"\n",
    "        Construct parameters for the coefficients and the \n",
    "        latent groups \n",
    "        Note that poly includes the intercept, so a quadratic\n",
    "        function should be poly=3\n",
    "        \"\"\"\n",
    "        super(LatGBTM, self).__init__()\n",
    "        self.coef = torch.nn.Linear(groups, poly, bias=False)\n",
    "        self.latent_flat = torch.nn.Linear(groups, obs, bias=False)\n",
    "        #setting the default weights for latent to equal prob\n",
    "        #torch.nn.init.constant_(self.latent_flat, 1/float(groups))\n",
    "\n",
    "    def forward(self, x):\n",
    "        \"\"\"\n",
    "        In the forward function we accept a Tensor of input data and we must return\n",
    "        a Tensor of output data. We can use Modules defined in the constructor as\n",
    "        well as arbitrary operators on Tensors.\n",
    "        \"\"\"\n",
    "        #here periods could also be long tensor with varying N\n",
    "        latent_expand = self.latent_flat.weight.repeat_interleave(periods, -2)\n",
    "        #softmax to scale probabilities to sum to 1\n",
    "        pred_groups = torch.mm(x, self.coef.weight) * torch.nn.functional.softmax(latent_expand, 1)\n",
    "        #sum up to get the final predictions!\n",
    "        y_pred = torch.sum(pred_groups, 1)\n",
    "        return y_pred\n",
    "        \n",
    "\n",
    "#Set up a basic object\n",
    "model = LatGBTM(groups=2, poly=3, obs=100, periods=11)\n",
    "#Initialize latent weights to equal probability (maybe do 0 here instead?)\n",
    "torch.nn.init.constant_(model.latent_flat.weight, 1/float(2))\n",
    "criterion = torch.nn.MSELoss(reduction='sum')\n",
    "#optimizer = torch.optim.SGD(model.parameters(), lr=1e-4)\n",
    "#SGD converged faster, but not even close right answers!\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=1e-4)\n",
    "\n",
    "#Now create the data tensors for my input data (time points)\n",
    "#And outcome data\n",
    "x_ten = torch.Tensor(t_np)\n",
    "y_ten = torch.Tensor(y_werr)\n",
    "\n",
    "#Needs alot of iterations to converge!\n",
    "for t in range(30000):\n",
    "    # Forward pass: Compute predicted y by passing x to the model\n",
    "    y_pred = model(x_ten)\n",
    "\n",
    "    # Compute and print loss\n",
    "    loss = criterion(y_pred, y_ten)\n",
    "    if t % 1000 == 99:\n",
    "        print(t, loss.item())\n",
    "\n",
    "    # Zero gradients, perform a backward pass, and update the weights.\n",
    "    optimizer.zero_grad()\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    \n",
    "##############################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can check out how close our actual coefficients are, as well as our latent trajectories.\n",
    "\n",
    "The order of the groups is not meaningful (here they line up, but if you run with a different set of starting parameters they may be swapped). But you can see we are kind of close in terms of coefficients. We have some coefficients that should be 0 off by ~0.1. \n",
    "\n",
    "For the latent classes, we need to apply the softmax function to get our probabilities. You can see that they did not converge to 0/1 (most are around 95% in the correct group), but are stacked up like you would expect. So even this simple problem did not perfectly recover the correct groupings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter containing:\n",
      "tensor([[-0.0722,  1.0735],\n",
      "        [-0.0400,  0.5390],\n",
      "        [ 1.6049, -0.1103]], requires_grad=True)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor([[0.0752, 0.9248],\n",
       "        [0.0583, 0.9417],\n",
       "        [0.0534, 0.9466],\n",
       "        [0.0596, 0.9404],\n",
       "        [0.0752, 0.9248],\n",
       "        [0.0708, 0.9292],\n",
       "        [0.0669, 0.9331],\n",
       "        [0.0837, 0.9163],\n",
       "        [0.0635, 0.9365],\n",
       "        [0.0672, 0.9328],\n",
       "        [0.0786, 0.9214],\n",
       "        [0.0796, 0.9204],\n",
       "        [0.0524, 0.9476],\n",
       "        [0.0792, 0.9208],\n",
       "        [0.0783, 0.9217],\n",
       "        [0.0504, 0.9496],\n",
       "        [0.0669, 0.9331],\n",
       "        [0.0698, 0.9302],\n",
       "        [0.0535, 0.9465],\n",
       "        [0.0793, 0.9207],\n",
       "        [0.0675, 0.9325],\n",
       "        [0.0741, 0.9259],\n",
       "        [0.0664, 0.9336],\n",
       "        [0.0662, 0.9338],\n",
       "        [0.0478, 0.9522],\n",
       "        [0.0613, 0.9387],\n",
       "        [0.0633, 0.9367],\n",
       "        [0.0669, 0.9331],\n",
       "        [0.0598, 0.9402],\n",
       "        [0.0669, 0.9331],\n",
       "        [0.0559, 0.9441],\n",
       "        [0.0628, 0.9372],\n",
       "        [0.0644, 0.9356],\n",
       "        [0.0660, 0.9340],\n",
       "        [0.0656, 0.9344],\n",
       "        [0.0828, 0.9172],\n",
       "        [0.0561, 0.9439],\n",
       "        [0.0616, 0.9384],\n",
       "        [0.0612, 0.9388],\n",
       "        [0.0551, 0.9449],\n",
       "        [0.0601, 0.9399],\n",
       "        [0.0723, 0.9277],\n",
       "        [0.0512, 0.9488],\n",
       "        [0.0725, 0.9275],\n",
       "        [0.0823, 0.9177],\n",
       "        [0.0838, 0.9162],\n",
       "        [0.0559, 0.9441],\n",
       "        [0.0623, 0.9377],\n",
       "        [0.0751, 0.9249],\n",
       "        [0.0872, 0.9128],\n",
       "        [0.9354, 0.0646],\n",
       "        [0.9514, 0.0486],\n",
       "        [0.9287, 0.0713],\n",
       "        [0.9467, 0.0533],\n",
       "        [0.9273, 0.0727],\n",
       "        [0.9348, 0.0652],\n",
       "        [0.9335, 0.0665],\n",
       "        [0.9215, 0.0785],\n",
       "        [0.9351, 0.0649],\n",
       "        [0.9454, 0.0546],\n",
       "        [0.9188, 0.0812],\n",
       "        [0.9417, 0.0583],\n",
       "        [0.9214, 0.0786],\n",
       "        [0.9531, 0.0469],\n",
       "        [0.9184, 0.0816],\n",
       "        [0.9431, 0.0569],\n",
       "        [0.9392, 0.0608],\n",
       "        [0.9396, 0.0604],\n",
       "        [0.9304, 0.0696],\n",
       "        [0.9326, 0.0674],\n",
       "        [0.9279, 0.0721],\n",
       "        [0.9455, 0.0545],\n",
       "        [0.9224, 0.0776],\n",
       "        [0.9405, 0.0595],\n",
       "        [0.9154, 0.0846],\n",
       "        [0.9319, 0.0681],\n",
       "        [0.9386, 0.0614],\n",
       "        [0.9252, 0.0748],\n",
       "        [0.9397, 0.0603],\n",
       "        [0.9562, 0.0438],\n",
       "        [0.9341, 0.0659],\n",
       "        [0.9253, 0.0747],\n",
       "        [0.9343, 0.0657],\n",
       "        [0.9377, 0.0623],\n",
       "        [0.9310, 0.0690],\n",
       "        [0.9333, 0.0667],\n",
       "        [0.9448, 0.0552],\n",
       "        [0.9569, 0.0431],\n",
       "        [0.9528, 0.0472],\n",
       "        [0.9332, 0.0668],\n",
       "        [0.9467, 0.0533],\n",
       "        [0.9386, 0.0614],\n",
       "        [0.9391, 0.0609],\n",
       "        [0.9454, 0.0546],\n",
       "        [0.9435, 0.0565],\n",
       "        [0.9249, 0.0751],\n",
       "        [0.9236, 0.0764],\n",
       "        [0.9194, 0.0806],\n",
       "        [0.9250, 0.0750],\n",
       "        [0.9258, 0.0742]], grad_fn=<SoftmaxBackward>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "##############################################\n",
    "#Model coefficients\n",
    "print( model.coef.weight )\n",
    "#Probabilities need to use softmax function\n",
    "torch.nn.functional.softmax(model.latent_flat.weight, 1)\n",
    "##############################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I make a plot to illustrate how close the estimated trajectories are compared to the true trajectories. And you can see they line up pretty well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAFzCAYAAADSXxtkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeZxN5R/H38c29p3RjGWQJcY2hskSlTWJUkirZPsliZBUVJKkBUWWLK1EEpWkZCdhsocxxjIjY2eGmTHL8/vjO3ONzFwzc889d5nn/XrNi3vunXOe+dzn3u95nu9mKKXQaDQajUbjXeRx9QA0Go1Go9GYjzbwGo1Go9F4IdrAazQajUbjhWgDr9FoNBqNF6INvEaj0Wg0Xog28BqNRqPReCH5XD0AMylbtqwKCAhw9TA0Go1Go7GEHTt2nFVKlcvoOa8y8AEBAWzfvt208+3Zs4d69eqZdr7ciNbQcbSGjqM1NAeto+OYraFhGMcye05v0dvh3Llzrh6Cx6M1dBytoeNoDc1B6+g4VmqoDbxGo9FoNF6INvB2aNCggauH4PFoDR1Ha+g4WkNz0Do6jpUaepUP3mxiYmIoVaqUq4fh0WgNHUdr6DjequG1a9cIDw/n6tWrrh6KxskULlyY6tWrU6BAgSz/jjbwdjhy5AiVK1d29TA8Gq2h42gNHcdbNQwPD6dkyZLUqlWLPHn0hqy3kpKSQnR0NOHh4dxxxx1Z/j09IzQajcZDuXr1Kr6+vtq4ezl58uTB19c32zs1elbYQefUO47W0HG0ho7jzRpq4547yMn7rLfo7VC6dGlXD8Hj0Ro6jtbQcbSGzuHcuXO0adMGgFOnTpE3b17KlZOaK3/99Ve2/MW34s8//2T48OGcOXMGwzBo1aoVU6ZMoVChQqZdw9vQt352CA0NdfUQPB6toeNoDR1Ha+gcypQpw86dO9m5cycDBw5k6NChtsdpxl0pRUpKikPX+ffff+nZsycffvghBw8eZP/+/bRp04bY2Ngs/b4ZY/BEtIHXaDQajakcPnyYwMBABg4cSFBQECdOnKBkyZK25xcuXEjfvn0BiI6Oplu3bgQHB9O0aVP+/PPPm8738ccf8+yzz9K0aVNAtqt79uxJuXLleO2115g8ebLttbVr1yYyMvKmMYwbN47Ro0fbXvfZZ58xdOhQAD7//HOaNm1Kw4YNee6550hJSSEpKYknn3ySevXqERgYyNSpU52ilTPRW/R28Ma0GqvRGjqO1tBxcoOGL74IO3eae86GDSGd7cwW+/fvZ968ecyYMYOkpKRMX/fCCy8wcuRI7rzzTo4ePUrnzp3Zu3fvDa/Zu3cvAwYMcGgMp06dolWrVrzzzjsAfPvtt4wbN469e/eydOlSNm/eTL58+ejfvz8LFy6kevXqnD17lj179gBw8eLFbF/f1WgDbwdd1MFxtIaOozV0HK2h9VSvXp0mTZrc8nW///47Bw8etD2+cOECcXFxpvjW04+hQoUKVKxYke3bt1O5cmUiIiIICQlhypQpbNu2jeDgYADi4uKoVKkSHTp04ODBgwwZMoROnTrRvn17h8djNdrAZ8aCBRxbsYIqX37p6pF4NOvWraN169auHoZHozV0nNygYU5X2s6iSJEitv/nyZMHpZTtcXx8vO3/SqlbBuTVrVuXHTt2cP/999/0XL58+W7wr6c/d/oxAPTs2ZNFixYREBDAww8/jGEYKKXo06cP48aNu+ncu3fv5pdffmHq1KksWbKEWbNm3eKvdi+0Dz4zQkOptGABREe7eiQeTfoPtSZnaA0dR2voWvLkyUOpUqUICwsjJSWFpUuX2p5r27Yt06ZNsz3emYGfYfDgwcyZM8fWLVQpxeeff86ZM2cICAhgx44dgETunzhxItNxPPLII3z//fcsXLiQnj172q6/aNEizp49C0hmwPHjxzlz5gxKKbp3786bb77pkYGa2sBnRp8+5ElOBr2CdwjDMFw9BI9Ha+g4WkPXM3HiRDp27EibNm2oWLGi7fi0adPYtGkT9evXp06dOsyePfum3/Xz8+Obb75hyJAh1K5dmzp16vDnn39StGhRunfvTnR0NI0aNWLOnDlUq1Yt0zGUKVOG22+/nX///ZegoCAA6tWrx9ixY2nbti3169enffv2REdHc+LECVq1akXDhg3p16+fzXfvSRjedGcbHByszOwHT/PmcPEi7NsH+gtCo9G4GTt27KBx48auHobGIjJ6vw3D2KGUCs7o9XoFb4cT7drBP/9ABmkbmqyxa9cuVw/B49EaOo7WUJMb0QbeDkdDQuDuuyEx0dVD8VguXLjg6iF4PFpDx9EaanIjOoreDsmFC8OaNa4ehkaj0Wg02Uav4O2QFoTBhQtw4IBrB+Oh2DTU5BitoeNoDTW5EW3g7XD+/Hn5zz33QP/+rh2Mh2LTUJNjtIaOozXU5Ea0gbfD0aNH5T+9esGGDXDokEvH44nYNNTkGK2h42gNNbkRpxl4wzDmGoZx2jCMvZk8f7dhGJcMw9iZ+jMm3XMdDcM4aBjGYcMwRjlrjFnmqacgb16YN8/VI9FoNBq3Im/evDRs2ND28+6772b62h9++IH9+/fbHo8ZM4bff//d4TFcvHiR6dOnZ/v33njjDd5///0Mn/viiy8IDAykbt261KlTJ9PXuTPOXMHPBzre4jUblFINU3/eAjAMIy8wDbgPqAP0MgyjjhPHmSm2ggm33QadOsH8+WCnaYLmZuwVndBkDa2h42gNnUehQoVsLWJ37tzJqFGZr8n+a+Dfeust2rZt6/AYcmrgM+OXX35h8uTJrFq1in379hEaGkqJEiWy/Pv2mutYidMMvFJqPZATx1dT4LBS6ohS6hqwEOhq6uCywPHjcPZsuje0b184dQo2bbJ6KB5NsWLFXD0Ej0dr6DhaQ+sZNWoUderUoX79+gwfPpzNmzezfPlyRowYQcOGDQkPD6d379589913AAQEBDB69GiaNWtGcHAwoaGhdOjQgerVqzNjxgwAYmNjadOmDUFBQdSrV49ly5bZrhUeHk7Dhg0ZMWIEAJMmTaJJkybUr1+fsWPH2sY1fvx4atWqRdu2bW9ocpOeCRMm8P777+Pn5wdAwYIF6devHwB33323rWTu2bNnCQgIAGD+/Pl0796dBx54gPbt29OzZ09WrFhhO2fv3r1ZsmQJycnJjBgxwja2mTNnAtLzPq1yXmBgIBs2bHD4PXB1mlwzwzB2ASeB4UqpfYA/kL6YcCQQYuWgEhMhKAgaNIhn9epUI9+pk0TS16pl5VA8nl27dnH33Xe7ehgejdbQcXKNhhn9jT16wHPPwdWr8j32X3r3lp+zZ+GRR258bu3aW14yLi6Ohg0b2h6/8sortGvXjqVLl3LgwAEMw+DixYuULFmSLl260LlzZx7573VSqVSpElu2bGHo0KH07t2bTZs2ER8fT926dRk4cCAFCxZk6dKlFC9enLNnz3LnnXfSpUsX3n33Xfbu3WurY79q1SrCwsL466+/UErRpUsX1q9fT5EiRVi4cCF///03SUlJBAUFZVgJcO/evTmqELhlyxZ2795N6dKlWbp0Kd9++y2dOnXi2rVrrF69mk8//ZQ5c+ZQokQJtm3bRkJCAi1atKB9+/Z8//33dOjQgVdffZXk5GSuXr2a7ev/F1ca+FCgilIq1jCMTsAPQA0go5qwmdbTNQyjP9AfpF7x2tQJWa1aNYoVK2arYFWmTBnq1q3L+vXrAelA1LJlS0JDQ7l8+TIAwcHBthrE99xTnR9+8GffvnOcOSP9gMuXL0/NpCQ2btwIgI+PD82aNWP79u3ExsYCEBISQmRkJFFRUQDUqlWLvHnz2ralKlSoQNWqVdmyZQsg21shISFs3bqVuLg4AJo1a0ZERASnTp0CoE6dOiQnJ9vuNv39/alYsSJbt24FoGjRogQHB7NlyxYSEhIAaNmyJYcOHeL06dMABAYGkpCQQFhYGCAfJF9fX9udaPHixQkKCmLjxo227aVWrVqxb98+zp07B0jLzZiYGI4cOQLIHXfp0qVtTRhKlSpFgwYNWLduHUopW/3vXbt22QqNBAUFcf78eVvQk6PvE0CNGjXw8fGx9ZAuX748NWvW9Jr3KTk52TavnfU+tW7d2qvfp/j4eJuGnvx5+u/75A47E2lb9OlJSkqiYMGC9O3bl/vvv5/OnTtn6VxdunQBpD58bGwsxYoVo1ixYhQsWJCLFy9SpEgRRo8ezfr168mTJw9RUVFEZ9AQbNWqVaxatYpGjRoBsvIPCwsjJiaGhx56iMKFC99wPbNo164dpUuXBuC+++7jhRdeICEhgZUrV9KqVSsKFSrEqlWr2L17t23n4tKlS4SFhdGkSRP69OlDYmIiDz744A03TelJm8dpnye7KKWc9gMEAHuz+NqjQFmgGfBruuOvAK9k5RyNGzdWZnHokFKg1Lhx6Q4mJirVo4dSEyeadh1vZ/fu3a4egsejNXQcb9Vw+/btrh6CKlKkSIbH4+Pj1c8//6yefPJJdc899yillHr66afV4sWLba9J/7hKlSrqzJkzSiml5s2bpwYNGmR7Xdpz8+bNUz169FDXrl2zHY+IiFARERGqbt26ttcPGzZMzZgx46YxffTRR2rMmDG2x0OHDlWTJk266XUtW7ZUq1evzvDvatOmjdq6datSSqkTJ06oKlWqZDhmpZR64okn1LJly1SvXr3U8uXLlVJKdevWTa1cuTLDc0dFRalZs2apwMBA9fnnn9/0fEbvN7BdZWITXZYmZxhGBSN1iWcYRlMkHuAcsA2oYRhGVcMwCgCPAsutHl+NGtCunWLmzHRxdfnywblzMGMGpOs/rMmcunXrunoIHo/W0HG0htYSGxvLpUuX6NSpE5MnT7at8IsVK0ZMTEyOz3vp0iXKly9P/vz5WbNmDceOHcvwvB06dGDu3Lm2naCoqChOnz5Nq1atWLp0KXFxccTExPDjjz9meJ1XXnmFkSNH2nbnEhISmDp1KsAN7WnTVuGZ8eijjzJv3jw2bNhAhw4dbGP79NNPSUwtgX7o0CGuXLnCsWPHKF++PP369ePZZ581pT2tM9PkFgBbgFqGYUQahvGsYRgDDcMYmPqSR4C9qT74qcCjqTckScDzwK/AP8AiJb55y2nVai+RkbA8/e1F374QEQF//OGKIXkcaVu4mpyjNXQcraHzSPPBp/2MGjWKmJgYOnfuTP369WndujUfffQRIAZv0qRJNGrUiPDw8Gxf6/HHH2f79u0EBwfz9ddfU7t2bUBcRi1atCAwMJARI0bQvn17HnvsMZo1a0a9evV45JFHiImJISgoiJ49e9KwYUMefvhh7rrrrgyv06lTJwYNGkTbtm2pW7cujRs3trlahg8fzqeffkrz5s1tPeQzo3379qxfv562bdtSoEABAPr27UudOnUICgoiMDCQAQMGkJSUxNq1a2nYsCGNGjViyZIlDBkyJNv6/BfdLtYOq1evpU+fu6lRA2ypmvHx4O8P7drBwoWmXctbWbt2be4IbnIiWkPH8VYNdbvY3IVuF2siPj75GDAAVq9OV4q+YEF48klYulSiTjV2yZfP1Ykano/W0HG0hprciDbwdmjZsiV9+0KBAnBDDYX+/WHkSDAyCvjXpKdly5auHoLHozV0HK2hJjeiDbwdQkNDKV8euneHzz+H1HgNqFMHxo2DMmVcOj5PwIxAkdyO1tBxtIaa3Ig28HZIy+d97jm4fBm+/jrdk0lJEn23e7drBuchpGmoyTlaQ8fxZg1TdEZPriAn77M28FmgWTNo2BCmTQNbTGJCgvjiPbABgUaj8Q4KFy5MdHS0NvJeTkpKCtHR0bYCPVlFR9HbITY2lqJFiwIwZ45kyK1dC61bp75g0CB54uRJSK1epLmR9BpqcobW0HG8VcNr164RHh5uSllTjXtTuHBhqlevbku3S8NeFL0OLbVDdHS07Uvhscckru7jj9MZ+P79Jfruq6/ghRdcN1A3Jr2GmpyhNXQcb9WwQIEC3HHHHZZdLzw8nOrVq1t2PY/j8GH4+Wf46SdYt04am5w6Bb6+kopVqhThsbGWaai36O2QVkMboFAhePZZ+OEHsB1u0ACaNoVZs9Lt3WvSk15DTc7QGjqO1tActI7/ITFR3LUAX3whJVBffBGiomDIEFiz5vrubu3a4OtrqYbawGeD554TO57auVDo31+6NKWWNNRoNBqNF3PhAixYINu65ctfj75u0wamTIHwcNi/HyZNku5++fO7bKjawNuhRo0aNzwOCIAHHpAFe3x86sGnnpJtmdtus3x8nsB/NdRkH62h42gNzSFX6xgfD/feC+XKiXFfvRoeegjSXCT+/uKqrVbN7mluv906DbWBt0NGrfgGD5YCdosWpR7Inx/y5JFtmtT2lJrr3LKdoeaWaA0dR2toDrlGx+Rk2LwZRo2CYcPkWMGCYtxffhm2bIF//4W5cyXN6hZcuQKffQZNmsAff1jX4lcbeDuk9a1Oz733yg3bxx+nc7tHR0PFivIOam4gIw012UNr6DhaQ3Pweh3XrZNgKz8/aNECPvgAjhy5/mX/7bcwfjzceacs7G7Bnj2SbOXnB/36yRowMvKYk/+I62gDn00MA55/HrZvh61bUw/6+kLVqjBzpg6202g0Gk/h7FkpU5rmc/39d/juO1nJLVgAZ85IZHU2ypLHxUm8XfPmUL++ZFJ36QIbN4rBb978nJP+mJvRBt4O5cuXz/D4U09B8eKyircxYADs2yfbOhobmWmoyTpaQ8fRGpqDV+h49ChMniwBcL6+0Lu3rNwBRowQo75gATz6KJQsmeXT/vOPBND7+8PTT8O5c7IBEBUFX34pGwKGYa2GutCNHZKSkjLtQvXii5ICf+xYanzdlSuyD/PAA5IXrwHsa6jJGlpDx9EamoNH6qiUrNALFYJdu6QsKUBgIHTtKoFyQUE5ah6WkABLlsjm7fr1EpLVrRsMHCj1UjI6pdka6naxOWTjxo2ZPvf881KO3pYyV6SILO0XL5ZbNw1gX0NN1tAaOo7W0Bw8RseUFAmEGzFCctNfekmO16sHU6dCWJjsl7/9NjRunG3jHhYmp/b3h8cfh8hImDhR/l24UDYHMjullRp62K2Y+3D77XD//WLgR48GHx8k2rJHD122VqPRaFzF66+L4/vff2VJ3aYNtGolz+XJI6lQOeDaNVi2TFbrq1dDvnyyATBggFwiCzF3lqMNvB1ulRIyZAi0ayd3bE8/jQTaVa1qzeA8hFyTVuNEtIaOozU0B7fTMSEBfvsN/vhDHN6GATExEuHWrRt06pQtP3pGRETA7NmSERcdDVWqyMK/T5+clT+xUkPtg3cApcSN4+MDO3akbsmcPy93kN27yz6NRqPRaMzj6lVYuVKi3X/6SQx6iRLiX69SxZRLJCXJqWfOhF9/le/2zp1ltd6hA+TNa8plTEH74HPIrW4WDEMKF/39t6RAAOKLX7ToPyH2uRcrb7i8Fa2h42gNzcFlOsbGwsWL8v9Vq+Dhh+XfHj1gxQo4fdoU437iBIwdK6d66CFx048ZI4H3y5bJhoCjxt1KDbWBt0NsbOwtX/Pkk1CqlMRtALKcf+YZmQ0nTzp3gB5AVjTU2Edr6DhaQ3OwVMeYGElX69ZNKsilLZo6dBAn+KlTUlzsvvvgPy1Us0NysjSA69JFypGPGyd9xH74QQz7G29ApUpm/EGClRpqA+8ghQtLhaKlS+H48dSD/fvLrJkzx6Vj02g0Go9DKejZ83rN961b5Uv2vvvk+UKFpBCNg6lmJ0+KMa9WTbbf//pLKtMeOSKbAl27OnwJl6N98HaIi4ujUKFCt3zd8eMSWzd8uKRKANC+vVQ+iIjw/FniAFnVUJM5WkPH0Rqag1N0jIkRh/euXfDuu3Ksf38x5D16SK13k0LUU1IkJm/mTFi+XNZhbduKb71rV2sav5mtofbB55DIyMgsva5yZfHXzJ4t9W4ACbG/7z7xHeVisqqhJnO0ho6jNTQH03S8ckXquj/8sLRcfewxKRAWEyPPz5olrVdbtDDFuEdHy71DjRrQsSNs2CBZzWFhYvAfecS6rq5WzkVt4O0QFRWV5de++KK0Cf7ii9QD998vk9TBFA1PJzsaajJGa+g4WkNzcEjHuLjrHTe/+EJKwW7eLNvvGzbIVmgx8zqtKQVr1shuf6VK8Mor8u+CBVKQ5r33pJ6J1Vg5F7WBN4kWLSA4WEocp6SkHlQK/vwznXNeo9FochHXrsn2+5NPyko9rYx3jx5ifSMjJUK5ZUvTtuHPnpWU+Nq1xVX/229SefSff2DtWrmvcLd0fmehDbwdatWqleXXGgYMHQqHDsEvv6QePHsW7rpLtppyKdnRUJMxWkPH0RqaQ5Z1TEqCvn2hQgXpz/Hzz7KUbtRIni9TRuqEmJRQrpRsAjzxhJSPHT4cypaVRnFRUfDhh2Lw3QEr56I28HbIm83J1727TK6PPko9UK6c+JjmzpXiDLmQ7GqouRmtoeNoDc0hUx2VktrvaZlD+fJJjtn998sKPi2lLTjDWLAcc+GCbADUrSvVaH/6SeLzdu+GTZukPYi7xVZaORe1gbfD/v37s/X6/PllK2j1aplgADz3nBRoWLjQ/AF6ANnVUHMzWkPH0Rqaw0067tkjzTiqVZPysC+9JOVjQfbGv/xSjLwDeer/Jc3z+cwz0sBzyBBx3c+ZI6v1jz+WnjLuipVzURt4k+nfX3LjJ09OPXDXXXJ7OW2azEyNRqPxBiZPhvr1JVqtdm3ZDz9+/LqDOwftV+1x+TJ8+qns8jdrJpVqn34aQkMlVb5PHykkqrmONvB2qFChQrZ/p3RpmXRffy2pGRgGDBok21W5sLJdTjTU3IjW0HG0hg5y5gxMm0bTYcNkHxykOsy0adK17ZdfZD+8eHHTL71jhyyc/PxkQzRPHuniefKk/Jvm1vcUrJyLutCNHRISEnLU+efgQbmhHTtWyhwSHy+rd3dzBllATjXUXEdr6DhawxyQlASLF0vk+6+/QnIyKXXrkuedd6SuqxOJjRWv5syZsH27fHX26iUFaZo0MX1zwFLMnou60E0O2bJlS45+r1YtcTtNny62nYIFZYampFz3T+UScqqh5jpaQ8fRGmaRpCRZoYAslUeNEj/78OGwezfrP/nEqcZ91y5Zpfv5SXp8fLz41E+eFB9706aebdzB2rmoDbyTGDZMdrXS0j6JiZFlvc05r9FoNG6AUuLIHjZMKsHcdRckJoqBX79e3Ivvvuu0yLWrV2H+fPGrN2woSUcPPigdOnfvlsDlXF4vLMdoA28HR+oF33OPTNYPP0wtfFOsGFSsKMv65GTzBunm6PrfjqM1dBytYSasWCFBwI0bwyefiJWdMeP6MrlKlRsK0Jip4/79EgHv7y8R8RcuSIrxyZNS6K5FC89frWeElXNR++CdyFdfSQGnn3+WPsIsXSqtD5culVtUjUajsZLLl2HJEggJgTp1pLTb669LhZju3SVK2IkkJMjlZ8yQwjT580upkIEDJY/dGw26s3GJD94wjLmGYZw2DGNvJs8/bhjG7tSfzYZhNEj33FHDMPYYhrHTMAyXWeytW7c69Ps9e8rd6QcfpB544AHZAvvkE8cH5yE4qqFGa2gGuVrDpCQJknv8caks16ePBM+BVJPbsEGi17Jg3HOqY1iYuPH9/WUYJ09K582oKKkN37p17jHuVs5FZ/YxnQ98AnyRyfMRQGul1AXDMO4DZgEh6Z6/Ryl11onjuyVxaY0Rckj+/PDCC/Dyy7BzJzRsmE8iSF55Rfan6tQxaaTui6MaarSGZpBrNUxJgcBACZwrVQp695Z0tpCQW/5qRmRHx2vX4IcfJBL+jz+kuF3XrrJav/de00rPexxWzkWnGXil1HrDMALsPL853cM/gYrOGosr6d8fxo2TVfyXXyL1mWvUgJo1XT00jUbjbZw5A998I2VjFywQKzp4MNx2m6T2WJAqeOSItM6eOxdOn4aAABg/Xvzst93m9Mtr0uFUH3yqgf9JKRV4i9cNB2orpfqmPo4ALgAKmKmUmpWV67lLHvx/GTpUduUjIiTOLjeh848dR2voOF6t4bVrEiw3f74E/CQlQVCQlIo12aeemY6JiVL/ZsYMWLVKesh07iyr9XbtTOsp4xVYmQfvcgNvGMY9wHSgpVLqXOoxP6XUScMwygO/AYOVUusz+f3+QH8APz+/xl9//TUA1apVo1ixYuzatQuAMmXKULduXdavl9Pky5ePli1bEhoayuXLlwEIDg4mOjqaEydOAFC4cGGqVavG3r0SRlC+fHlq1qzJxo0bAfDx8aFZs2Zs376d2NhYAEJCQoiMjLT1/K1VqxZRUfkJCSlDjx4nGDv2KlUrV+bkoEHE+/py+aGHCAkJYevWrbatm2bNmhEREcGpU6cAqFOnDsnJyRxMzU/19/enYsWKNl9O0aJFCQ4OZsuWLSSk5tm3bNmSQ4cOcfr0aQACAwNJSEggLCwMgEqVKuHr60vaDVHx4sUJCgpi48aNJCUlAdCqVSv27dvHuXPnAGjQoAExMTEcOXIEgICAAEqXLk1oaCgApUqVokGDBqxbtw6lFIZh4OvrS0JCAhcuXAAgKCiI8+fPc/ToUdPepxo1auDj4+Pw+5Q3b15bnegKFSpQtWpVW85qoUKFXPY+XblyhfPnzzv1fWrdujW7du3y2vdp27ZtNr08+fN0w/uUnExQkyYkzJ9PuSFDuFaqFPHdu6Oeeoq/ExOd8j7Fx8dTokQJ2/sUHp7Izz/fxm+/VeHUKYNy5eK5//5/GTy4MBUrGm75ebL8feLGz1N8fDx16tQx7fNUrly5TA08Simn/QABwF47z9cHwoGadl7zBjA8K9dr3LixMpM1a9aYdq4ePZQqUUKpy5dTD9x5p1I1aiiVnGzaNdwRMzXMrWgNHcdrNDx9WqnJk5Vq0ECp996TY1euKPXTT0olJjr98mvWrFFJSUr9+KNSnTsrlSePUoahVKdOSi1bZskQPB6z5yKwXWViE10W5mAYRmXge+BJpdShdMeLGIZRLO3/QHsgw0h8T2L4cLh0SXxTgFRvCAuT/SyNRqOxx4oVkk/m7w8vvmJPnuIAACAASURBVCgRvFWqyHOFC4t/PZ8zY6bT8tOrULWqJARt3y7xwkeOiGegSxenD0GTTZy2RW8YxgLgbqAsEA2MBfIDKKVmGIbxGfAwcCz1V5KUUsGGYVQDlqYeywd8o5Qan5Vrmu2DP336NOXLlzftfPfcA4cPywcif0qCfECDguTD66WYrWFuRGvoOB6pYVSUGHSADh3g77+lsMYzz0hkvAWkpIgrf+ZMWL5canS1aydZdV26yH2GJnuYPRft+eCdGUXf6xbP9wX6ZnD8CNDg5t+wnmSTK86NGCE32gsXwpNP+sD//ifdaA4elAL2XojZGuZGtIaO4zEaxsbCokUSgv7nn3DsmBj5efOgXDnLLGp0tAxh9mwJDi5bVlq9d+16mubNPexGyc2wci7m0kzErJEW3GEW990nN96TJqW2hh84EB56yKtL15qtYW5Ea+g4bq/h8ePw7LNSiObZZ+HcOZgwQbbfQbqvONm4p6TA6tXQo4dk+4weLZuMCxdCZKQUprl2bb9Tx5AbsHIuao+JhRiG+OJ795bCUh07+sL337t6WBqNxhWcPi2G/I47xHm9ZAk8+qhUmmvWzLLSbmfPwuefyzZ8WJhk1g0eLNvwXrqxmGvQK3g7+Kf5v0ykVy/ZcZs0Kd3Bw4elJrQX4gwNcxtaQ8dxGw2Tk2HlSnjkEfkiGDxYjvv5yb74Z59B8+ZON+5KSYXaxx+XYQwfDr6+UowrKkqaZGVk3N1GRw/GSg31Ct4OFZ1QlaZAAQmCHTFColCDg5E79qgoOHTI6ypCOEPD3IbW0HHcQsNZs+Dtt+HECXFqv/CCbMenYUEhngsXpFPbrFlSLbtECVmp9++ftbg9t9DRw7FSQ72Ct4OzmgL07w/Fi6dbxQ8efD3XxMvI1U0+TEJr6Dgu0TAxUVxwV6/K46tXZTt+0SK5of/gA0v6USgllWt795aNghdfhKJFJYju5EmYOjXrQfl6LjqOlRpqA+8CiheX+LrvvoPwcKR1bMWK8knTaDSeTXi4JIhXqiS560tTs36HDJHgm+7dZSvPyVy6BNOnQ4MGsuu/ZIkY+b//hq1bJdsuLYZP451oA2+HokWLOu3cQ4ZIXM0HHyDRsYMGSQjrXo+v6XMDztQwt6A1dBxLNLx8Gdq2hdtvl+25kBD48UfpGw2WBc1t3w79+slqfdAg+XqZOVNW659+Cg0b5vzcei46jpUaOrUWvdWYXejG2fTvL/6wY8fAN985+WL46CO5zdZoNO7P4cPSC/qRR2Qv/MEHJbCmT5/rRWosIDZWmsfNmAGhobIy79VLdgqDM65SrvES7BW60St4O6Q1RnAWI0ZII6gpU4AyZcQv52XG3dka5ga0ho5jqobXrsHixbJar1FDWkDHx8sKfdkyeP11y4z7rl3w3HOyWu/fX9z+n3wiq/XPPjPfuOu56DhWaqgNvB3SOhQ5ixo15MZ/2jTxl9kcYmfOOPW6VuJsDXMDWkPHMU3DH3+EypWlGkxYGIwbJ+HoBQuac/4scPWqdIZt1ky22+fNk3pZmzaJwR80SKLjnYGei45jpYbawLuYl18W193MmekO1KsH+oOk0biepCQpwr5zpzyuWhWaNpWMlyNH4LXXZPlsAfv2SWadv78EyF28KB69qCgpVGNB+rzGw9A+eDskJSWRz4L2SO3bw549UvO54Ibf5MD8+fD0006/trOxSkNvRmvoONnWMG2Pe/ZsqdM6YIA4uC0mPl6ybWbOhI0bJfj+4YdlOK1aWW/Q9Vx0HLM11D74HHLo0KFbv8gERo2CU6ck4I62baFuXbk194KbL6s09Ga0ho6TLQ3/9z/Zhh87VvLUv/9eHNsWcvCgNHfx95cGcqdOSWB+ZCR88w20bu2a1bqei45jpYbawNvh9OnTllznnnugSRN47z1ITjGkEsWuXV5RvtYqDb0ZraHj2NXwwgVZIqc1fQoIgGHDJEL+11/FwW3BqvXaNfj2W7j3XqhdW8pitGkDv/8uBn/4cGko50r0XHQcKzXUBt4NMAxZxYeHS3Aujz8upSxtjnmNRmM6O3ZIqVh/f8kn27xZjr/8stxtV69uyTCOHJHPf6VK0msmIgLGj5eKtosWiZHPo7+pNTlA++DtcPbsWcqWLWva+eyRkiI78wUKSDyP8ddWOeDhhSWs1NBb0Ro6zg0aRkVBt27w11+SufL445Jr5kgFmGySmCgB+TNnwqpV0oLigQfEt96+vfsadD0XHcdsDe354HW0hB0sTWfII9Utn35aAnQ7dw6x7NrORKfVOI7W0HGSDh+WCjDt20vP9RIlZA/8qaecl1OWAceOSdzenDniV69YEd588/pGgruj56Lj6DQ5NyEsLMzS6/XqJe6/8eNT4+vWr5dKFefOWToOM7FaQ29Ea5hDUlLEh96lC77Nm0tuWXKyLJdXrZImTxYY9+RkWa3ff79k2b3zDjRuLNl3EREwZoxnGHfQc9EMrNRQG3g3In9+GDkS/vwzNb6uTBnxE376qauHptF4FitWSKRax46wdSvHH39cPlgWtmOOioK33pKb9i5dpMnLa6/B0aPw00+yJa8zzjTORBt4O1SqVMnyaz7zjOwgjh+P+OA7dpQUnfh4y8diBq7Q0NvQGmaRffskMg2gUCEJVP36azh+nKQ33pAoNieTkgIrV0pJ+ipVbsy0O3ZMDH7lyk4fhtPQc9FxrNRQG3g7+Pr6Wn7NggUl/3X1amnpyEsvQXS0JL96IK7Q0NvQGtohKUnasd57rzQ1/+ADOX733RIV/9hj4OPjdA1PnZKt9+rV4b775NIjRkhmTFqmXf78Th2CJei56DhWaqgNvB1c1Zlu4EAoVUq+MGjTBurXly8uD8x48KTufu6K1jATpkwRi9qtm+SsT5gge+BwUxUYZ2iYkiI34j16yObAq6+Kj33hQilIM2ECVKtm+mVdip6LjmOlhtoD5IYULSr94t94A/bsNaj37rtSeDolxVIfokbjdoSFSZcmgN27xYJ+9JE4uS1yaJ85I5WkZ82S+4rSpaVGfP/+UKuWJUPQaLKENvB2KF68uMuuPXgwvP+++OIXLrzPZeNwFFdq6C3keg3TwtCnTJHo07/+ktKPM2Zked/bUQ2VkqSWmTNhyRKpOteypdyEP/ywpc3kXEqun4smYKWGutCNG/PKKzBxonSjrO0fI31lO3WSLXuNxtu5ckWM+CefSOh55crSC7VfP/FhWcD589IjYuZMOHAASpaU1Pn+/SUGVqNxNbrZTA7ZuHGjS68/bJgEA7/zDrKKGT9eLL4H4WoNvYFcp+HVq/KvUvD222LYlyyRiLWRI3Nk3LOjoVISJPf005KfPnSoGPZ58yT1bcqU3Gvcc91cdAJWaqgNvB2SkpJcev1y5STg7ptv4PDZkrJy+fZbybfxEFytoTeQKzRUSorP3H+/bL+npEgwyoEDsG6dBNI54GPPioaXLskmWYMG0KKFBOc/84yUjt6yBXr3lsq2uZlcMRedjJUaagPv5gwfLt9rEyYgXeYMAyZPdvWwNBpzuHpVotUCA6FDByns1LOnFGsHcHJKkVLi0n/2WfDzg+efF7f+rFnSEn76dDH4Go0non3wdkhJSSGPG3R9eOEFKWYXFgYBY56SqhnHj0v4rpvjLhp6Ml6t4bffSgu1Ro3kBrZnT/DxMf0y/9UwJkZ2xmbOlApzhQtLyvyAAVIdWpMxXj0XLcJsDbUPPofs27fP1UMAxO2YJ0+q+33ECAnfvXDB1cPKEu6ioSfjVRqGhsKTT8KkSfK4WzcJT9+xQ6LXnGDc4bqGO3eK28vPT/5NTpZt+X//lSYw2rjbx6vmoouwUkOdJmeHc27S5KViRfEFzp0Lr71WD/+VK109pCzjLhp6Mh6vYUqKFF//8EPxpxctej1hPH9+uOsup17+6lX46qv8tgy7ggVl02DAAAgJuakmjsYOHj8X3QArNdQreA9h1Cj5nrQF0R87Bps2uXRMGk2W6N8funaV1mmTJkm9+LSKc05k3z5xb/n5wXvv1SYmRsJXTp6UiPg779TGXePdaB+8HS5cuEApi/Jts0K/fvDll5It5N+9uZTUOnDAravbuZuGnojHaXjypOx79+snrdS2bpUbUgcj4bNCfDx8952kz2/aBAUKwCOPwOOPx3DffcW0QXcQj5uLbojZGmoffA6JiYlx9RBuYPRo8RlOnIgk5x4+LLk8boy7aeiJeIyGu3dL8nhAgKR9/PGHHA8JkYLtTjTuBw9K3Qh/f3Hxnz4tlSCjoqShXGDgBW3cTcBj5qIbY6WG2sDb4ciRI64ewg1UrSq5uLNmQVRIN2m08d57bt2Ext009ETcXsOUFMlfb9BACtIMHCgpH336OPWyCQnS2OWee6T1+8cfS2+m1atlY+ull6RjLHiAhh6C1tFxrNRQG3gPw7aKfz+vRNRv23Z9paTRWMW1a7Bihfw/Tx5pev7OO5K+OXWq3Hw6ifBwePllCT7t1Ut2/ydMkA5uixZJ51idyaXRONnAG4Yx1zCM04Zh7M3kecMwjKmGYRw2DGO3YRhB6Z572jCMsNSfp505zswICAhwxWXtcsMqvu3TUsZz1y5XDytT3FFDT8OtNLx4UXxEVavKqn3PHjk+aZI0T3BSbYbERNkcaN8ebr9duiffdResXCmeqlGj7NfEcSsNPRito+NYqaGz73PnAx3tPH8fUCP1pz/wKYBhGKWBsUAI0BQYaxiG5ZEdpd20kIxtFT+l4HXno5virhp6Em6h4blzsuddqZJY0zp1xLoGBjr1sseOScB95coSLHfgALz1lmwUfP+9FL/LymrdLTT0ArSOjpGSYq2GTjXwSqn1wHk7L+kKfKGEP4GShmHcBnQAflNKnVdKXQB+w/6NglMIDQ21+pJZ4oZV/LnUPpVu6htzVw09CZdqeOWK/JsvnzRB79pVitX89ptYVydEriUlwfLlskFQtapsvwcHS8fYiAh4/XVJfcsOeh6ag9Yx51y9Kk2KJk6MtOyarvZU+QMn0j2OTD2W2XFNKmmr+HffRSrg3H67rOY1GkdRSuI6OnaE5s3lcYkSspz+6ispK+sEIiOlv3pAgNxH/P23rN4jIsS4d+7s1hmhGo1dPv1UdqB8fRMsu6arK9lldPuv7By/+QSG0R/Z3sfPz4+1a9cCUK1aNYoVK8auVP90mTJlqFu3LuvXrwcgX758tGzZktDQUC5fvgxAcHAw0dHRnDgh9xYFChTg7Nmz7N0rIQTly5enZs2atnZ/Pj4+NGvWjO3btxMbGwtASEgIkZGRREVFAVCrVi3y5s3L/v37AahQoQJVq1Zly5YtABQqVIiQkBC2bt1KXFwcAM2aNSMiIoJTp04BUKdOHZKTkzmYasD9/f2pWLEiHTpcYObMCjwwtyrtfXw4PXw4+196CYCWLVty6NAhTp8+DUBgYCAJCQmEhYUBUKlSJXx9fUmrG1C8eHGCgoLYuHGjrdtRq1at2Ldvn63yUoMGDYiJibFFgQYEBFC6dGnbXX2pUqVo0KAB69atQymFYRiUKlWKXbt2cSG1tG5QUBDnz5/n6NGjpr1PNWrUwMfHx23fp61btwJQtGhRgoOD2bJlCwkJCVl+n4oWLWqb1856n1q3bs2uHTvIt2wZlRcupFhYGMlly3LswQc58fvvVK1VS96n1HGY9T5Vr16TyZP3s3y5H3/+WQalDEJCLjFgwAmaNTtHixZNiYyMZO1ax94nwzBsGjrrfbLi89S6dWuXfp7i4uLYsmWLR3+eXPE+1awZxIQJRQgOvkSNGqc4fvy4ae+TXZRSTv0BAoC9mTw3E+iV7vFB4DagFzAzs9dl9tO4cWOVmzh6VKn8+ZUaMEAp9fzz8uDECVcPS+OpfPONUqBUzZpKzZ6tVFyc0y518qRSb7+tVJUqcsny5ZUaNUqp8HCnXVKjcRkTJsg837LF/HMD21UmNtHVW/TLgadSo+nvBC4ppf4FfgXaG4ZRKjW4rn3qMUtZt26d1ZfMFlWqSLGwOXPgRI+XJILjww9dPawbcHcNPQGnaXj+PLz9tnRZAXj4YVi2DP75B/r2laLtJpKSAr//LsFylSvL9nv16pLaduKE+NqrVTP1kjb0PDQHrWP2uXxZkkw6dZLyyFZq6Ow0uQXAFqCWYRiRhmE8axjGQMMwBqa+ZAVwBDgMzAaeA1BKnQfGAdtSf95KPWYpyo0LyKQxerT4JcfOC5B+l99+KznKboInaOjumK5hZKRkXlSuLBFrf/0lxwsUgC5dTE8iP3NGvuBq1oR27WDtWukMe+iQFKXp3l0u7Uz0PDQHrWP2mTJF7qXfekseW6mhU33wSqlet3heAYMyeW4uMNcZ48oqhgfUtvT3h//9T6p4vbZuItWmFnL+t2U28AQN3R1TNfzgA8lXT0mRlmovvwz16pl3/lSUki6wM2dK/vq1a5K3/uabslFg8ubALdHz0By0jtnjwgX5yHXtCo0bpcC1JEs11M1mvIDoaEkn6tZNgpxRSiqDuJGh17iQv/6S1XqFCrJHvmyZ5LQ7oeDG+fPw+eeSwnngAJQsKeXp+/eX1HmNJjcxZgyMGye1yOofXQ6DBslnMK1dsgnoZjM5ZJcbV4hLj68vPP88fPMNHNhxBRo2lD1RN8BTNHRncqShUvJF0qaNNHuZOlWOt20r2z0mGnelpHPbU09JfvqwYWLY58+XZi+TJ7veuOt5aA5ax6xz7pzM/e7doX49JfEu+fOzKzXzwAq0gbdDWoqDJzByJBQpAmMmFpEi3R99dL1IiQvxJA3dlWxruGwZNG0qDu9//rleRtZkLl6ETz6B+vWhZUv44Qd49llZrWzZIiv3woVNv2yO0PPQHLSOWWfiRIiNldoOGIZY++nTuaC7yWmyS9my0kF28WI42P01uX2cNcvVw9JYRXLy9f8vXCjWd9YsqRIzfDgUK2bKZZSS9u59+shqffBg8afPnn29DXz9+qZcSqPxWE6elI2yJ59Mt3vVvLkUj7IQ7YO3w+XLlylevLhp53M2ly6JL75ZM/g57l6pbBcebn1EUzo8TUN3xK6GcXFSyfD996XcW2CgRPYUL25q2beYGHEBzZgBO3fKbtFjj8GAAdC4sWmXcRp6HpqD1jFrDBok99cHD0K1qA3w9dfSbbF0adM11D74HHL+vOWZeQ5RooT0AVmxAvZ0eVVuI7/5xqVj8jQN3ZEMNbx8WfYAAwIkAMPPTxqkA5QqZZpx37lT2rv7+cm/SsH06TK1Zs3yDOMOeh6ahdbx1hw5Ip+Nfv1S6zqMGyf+q0KFAGs11AbeDmnlHz2J55+XYOlBS+5Ffb9U9ohciCdq6G7cpGFionStGDVKAirXroWNG02ztleuyKZASIiUnf/iCylOs2WL1If/3/9kg8CT0PPQHLSOt+bNN6U302uvIf6s334TN1mqgbdSQ23gvYzChaV2yYaNBqsKPwj587t6SBozOHVKAieVkvf07bcl/e3XX6F1a1O6uu3dKz51Pz8JlouNlSIdUVEwb55U4dJp0BpN5uzfL6nKaZtqjB8PpUvL9pcL0AbeDtWcVTfTyfTtKzu3o0eD+nYRtGolPThdgKdq6DYcP07jefPkDR0+HPbtk+NPPw1Nmjh8+rg4WaG3aCH1bmbPhgcekCI1e/fCCy/Ijr+no+ehOWgd7TNmjMSnvPwy4t/68UcYMgSKFrW9xkoNtYG3QzGTIo+tpkAB2SYKDYXNoQVhwwaX+eI9VUOXc+6c3KlVr07RBQvE1XLwoATRmcCBA5J14e8v9wpnz0qcXlSUrEDuusu7Vut6HpqD1jFzduyQqo0vvSRZTZQuDc89J3fJ6bBSQ23g7eDJRR0ef1zSM5794QFUg4aypZs+lcoiPFlDl5BWu6BIEenJPnAgf375pSytb7/doVMnJEgG3d13wx13SEpb+/ZymQMH5IupTBnH/wR3RM9Dc9A6Zs6oUdfTlQGpHjltmlR9SoeVGmoD76XkzSvun4OHDH5vPgbCwqQRjcY92b0bevSQJPLEREltPHgQPv6YBF9fh059+LAUQqpYEXr1guPHpXPbiRNi8O+5x7tW6xqN1fz+u/y8+mpqAOrEibBtm6uHpQ28Pcp4+HKma1fJie+zrCspgfUkXcPiVbyna+h0QkPhoYegQQNYuVIawKR1A0wNkMyJhomJ8N13UsyuRg3pInzXXRKTd/iwrDYcvG/wKPQ8NAet480oJYUiK1eWDBP27JEP2PLlGb7eSg11oRs7pKSkkMfk1plWs369BFkv6PcHjz6SJN/4Fi7XvEFDp7F5s0S3lSwp/VMziWjLjoZHj8pu/pw50oSocmXJx02rPJdb0fPQHLSON7N4sWy+zZ8v8Sz06CE360ePih/+P5itoS50k0PWr1/v6iE4TKtW0KkT/G/xvVxo0t7yvVhv0NBUNm+WqlYgeWfTp8OxYzB2bKbh6rfSMClJys936iSFNd59V0rR//STFN147bXcbdxBz0Oz0DreSGKibMsHBsITTyBZLt99J/mmGRh3sFZDbeBzARMmSBnbD8ZdlVSr77939ZByHxs2SCe3Fi0kxSE5GfLkcahqTGSkNLIICIAHH5QmL6+/LguH5cvh/vtNrVar0Wj+w9y5Et70zjupn7Vx4yRAdtgwVw8NgHyuHoA7ky+fd8hTv75E1X803YexVX4h/4oV4qC34NvfWzTMMTt3yod9zRpxen/wgRRwz4b26TVMThY/+syZskJXCjp0kK5unTtLBS3NzeT6eWgSWsfrXL0q9+otWshnD6Ukdal+fbvpKFZqqH3wuYSICKhVC6a2XMTANT0lL75XL1cPy3tJSAAfH6k29+CDEsbev3+O+6f++6/41WfPlih4X1/xq/frJw2GNBqNtUyYIMXENmyQdsmuQvvgc0hoaKirh2AaVatKzYXn1z5CfI3A69vETsabNMwSa9dKovlzz8njpk1lz/zFF7Nt3FNSYNUqaNPmApUqyfZ7zZoS1HP8uGwLauOeNXLdPHQSWkfhzBkx8F27phr3f/4R12cWFsxWaqgNvB0uX77s6iGYymuvQZFiefio+BuSY71ggdOv6W0aZkqaYb/nHjh06MbGLwUKZOtUp09LGm2NGrL9vmNHEYYNk9P+9ps0fsnmKXM9uWYeOhmtozBunGzRv/tu6oExY6B3b7h48Za/a6WG2sDnIsqWlXzNV3c8xPGewz2n16e788EH1w37lCkQHn59BZ9FlJJ7hEcflYI0o0bJv998A4sWbeG998TgazQa13L4MHz6qVSSrl0biW797jspYedmjRu0D94OsbGxFE3XJMAbiIuTbd4KFaSTobNTWr1RQ0Acb8WLS4Ga48dh6VLxsae2hMwq58/D559L0NzBg5IS37u3nOqOO+Q1XquhhWgNzUHrKGnuK1aIoa9QASlUtWaNuOL+U5Y2I8zWUPvgc0h0dLSrh2A6hQpJWfrt2+GXyQeluIoTO815nYabN0uxoFatxAkHUk1myJAsG3elYNMm6R/j5ydB9qVLi6E/eVK6wqYZd/BCDV2A1tAccruOW7dKDMzw4anGfccO+OEHaeSQBeMO1mqoDbwdTpw44eohOIUnnpCF5w8TpdY5X3zhtGt5jYbbtsF990lOzK5d0npt7txsneLiRZG7Xj0JzFm+XLb5du+W+4annsr4HsFrNHQhWkNzyM06KgUjRkgGy0svpR48f15cnUOGZPk8VmqoDXwuJG9eeO89+Oz0A/xbualE1CckuHpY7s2vv4qRnzhRcg5feilLUfFKyV1/WqnYF14QIz57tqzWP/lEDL5Go3Fvli8Xz9wbb4Ct42u7dvK9kMNiVc5GG3g71PDiqKb27aF9e4NB598WH/JnnznlOh6r4e7d4ltbvFgeDx0qhn3kSKlUdQsuX5ZAnEaNpCLtokWyJb9jh3wf9O2bpdMAHqyhG6E1NIfcqmNionz0a9eGZ59NPbh4McTHZ7v8t5UaagNvBx8fH1cPwalMmgQ/xLYl3L+VOOavXjX9Gh6n4T//QM+e4sP44w+4cEGOFymS7rY9c0JDpVCdn58E0huGGPqTJyWQLigo+0PyOA3dEK2hOeRWHWfMkCSZ999PbfK4dq1E282Zk+1zWamhNvB22Lt3r6uH4FTq14dn+xo8e2o8Fzo95pRgO4/SMK1rxIoV8v+jRyWc/RZcuSKf8yZNxB335Zfy2d+6VQz+wIGO7eB5lIZuitbQHHKjjhcuyLZ8mzbS0Aml5PvBz098b9nESg11YeFczrhxcPuClvS92JIl7ulGci7Hj0uBgMKFxRk+bJjsxZUrd8tf3b1bVuVffSVb8nXrwtSpshWfxYBajUbj5rz9thj5Dz5I3Y1f8YtExc6Yke20WKvRK3g7lC9f3tVDcDoVKkjxm++/h90f/m56RL3banjqlES81aghLVtBqsxMmmTXuMfFSTpb8+ayiz9nDnTpIsE3e/ZIl0izjbvbauhBaA3NIbfpePiwZL706SOfd1JSZPVerVqOVu9grYa60I0dkpKSckX3pLg4aUTzxeUHaa3WYEREZNrLOLu4nYbnzkkKwccfw7Vr8iF97TXJZbfDP//Iav3zzyXdrVYt8bU/9ZTdxlGm4HYaeiBaQ3PIbTo+8gisXCktYW+7DSlC36WLBNg8+WSOzmm2hrrQTQ7ZuHGjq4dgCYUKSc2WwZfGQUyMpIKZhNtp+MQTskrv1g0OHIBZszI17gkJUiq2dWvpAjl9OnTsKEWr/vlHAuudbdzBDTX0QLSG5pCbdNywAZYsgZdfTjXuILt7mzdL/+0cYqWG2sBrAOkcW6hJPb4v+Dhq6lSIinL1kMzh6lUJff33X3n87ruyl/7VV3D77Rn+SliYFLSoWFE+x5GR8muRkdKf5+67s50Zo9FoPIiUFLmB9/dPV9Rm61aIjpYPv7NrfJuEZ4zSReSmlJA8eaRE6oi4N0m+lizRdybgMg2vXZP8wCPnYwAAIABJREFUtNtvF2v9/fdyvEEDiYbL4OWLF0PbtlKr/6OPZOW+apUY/JdfBle5H3PTPHQWWkNzyC06fv651KyYODG1nlVCgqTG9Ozp8Lmt1FD74DU30KsXBC4ey8Cx5Snz+iBXDydnfP21NE+PiJCasOPHS+34DIiIkKpyc+fKzXnlytCv3/XKcxqNJndx+bLc5FerJj0jDAOYPFmW9L//LvlyboT2weeQ3HizMHEijC/wJoP2mWPcXaLhihUSyr5iBaxff5NxT0qS/hAdO0L16vI3N20KP/8MR45IzJ07GffcOA/NRmtoDrlBx3fekZv9KVNSjfvly7JIaNvWFONupYZONfCGYXQ0DOOgYRiHDcMYlcHzHxmGsTP155BhGBfTPZec7rnlzhxnZsTGxrrisi6lcmXZjv7u2yQOvTxHGqs4gCUa/vGH5K3t3i2PZ8yQdnn33XeDs/z4cRgzBqpUkSq0e/fK46NHpc50p05Sp9/dyI3z0Gy0hubg7TqGh4t77qmnpHAVIAnwZ89e7x7pIFZq6LR8B8Mw8gLTgHZAJLDNMIzlSqn9aa9RSg1N9/rBQKN0p4hTSjV01vg0mTNiBCyafYXyH4xE7WqKsfIXVw8pY7Zvh9Gj4bffJCLu1Ckpz5eupGxysqS5zJghC3qlZOU+fTrcfz/koowfjUZzC4YPl1K0N9jykyfF/x6c4S64e6OUsvsDPA+UutXrMvi9ZsCv6R6/Arxi5/WbgXbpHsdm95qNGzdWZnL16lVTz+dJLFig1EtMUgqUWrMmx+dxmoZPPy1jK1tWqQ8/VCou7oano6KUGjdOqcqV5WW+vkqNHq1URIRzhuNMcvM8NAutoTl4s46rV8t3xfjxGTyZmGjadczWENiuMrGJWdmir4CsvhelbrlnNUHIH0jf+DYy9dhNGIZRBagK/JHucEHDMLYbhvGnYRgPZvGaphIZGemKy7oFPXtC6J2DiMxTiaSXRsrSNweYqmF09PVx1KgBY8fKntrQoVCwICkp0tW1WzdxNbz+ugTLLF4MJ06IGy0gwLzhWEVunodmoTU0B2/VMTFRWroHBEi1akC+NNLqxpu41WelhrcctVLqNcMwXgfaA88AnxiGsQiYo5QKt/OrGd0IZGYlHgW+U0olpztWWSl10jCMasAfhmHsyeh6hmH0B/oD+Pn5sXbtWgCqVatGsWLF2JXqQy5Tpgx169Zl/fr1AOTLl4+WLVsSGhrK5cuXAQgODiY6OpoTJ+S+5Nq1a5QqVcrWHKB8+fLUrFnTVqjAx8eHZs2asX37dptfJSQkhMjISKJS88hr1apF3rx52b9fPBMVKlSgatWqbNmyBYBChQoREhLC1q1biYuLA6BZs2ZERERw6tQpAOrUqUNycjIHDx4EwN/fn4oVK7J161YAihYtSnBwMFu2bCEhta97y5YtOXToEKdPnwYgMDCQhIQEwsLCAKhUqRK+vr62gI/ixYsTFBTExo0bSUptOjNxaitea/oW80OfYd8bb+D34ovExMRw5MgRAAICAihdujShoaEAlCpVigYNGrBu3TqUUhiGgVKKq1evciG1K1tQUBDnz5/n6NGjWX+f6tQheuhQyi5cyL433qDa4MFEP/qovE+hoZQsWYslS4ozd24eTp4sRKlSSbz4Yh4aNdqGv38cPj4+5M/vue/T8ePHbePM6H1q1aoV+/bt49y5cwA0aNAg2+9T69at2bVrl2Pv0y0+TzVq1MDHx8cln6fw8HDbOVz1efKG9yk2NpazZ8963ffe8OFH2Lu3Gm+9tZcCBeqwZ88+yg8dStlNm7h84AAxhmHa+xQbG4uPj49p75NdMlva//cHaABMBg4AnwJ/A+/ZeX2Wt+hTz9XczrnmA4/caoxmb9GvcWBr2lsY0DdJfWv0UBGf/Z6j33dIwytXlHrnHaVKlFDKMJR66imljh5VSimVkiJbaj16KJU/v2yttW4troX4+Jxf0h3R89BxtIbm4I06RkfLV0z79vK9opRSavt2+VIZPdr065mtIXa26LNi2F8AdgC/At2B/KnH8wDhdn4vH3AE2XovAOwC6mbwulrAUVJz8lOPlQJ8Uv9fFggD6txqrGYb+JMnT5p6Pk/kzBmlSpVS6p570k3+bJBjDVNSlAoOlinaubNSu3fbxvP++0rVrClPlSql1NChSv3zT84u4wnoeeg4WkNz8EYd+/RRKl++dN8hKSnyhVe2rFIXL5p+PbM1tGfgs+KDLwt0U0p1UEotVkolpq78U4DOdnYGkpAAvV+Bf4BFSql9hmG8ZRhGl3Qv7QUsTB1oGncA2w3D2AWsAd5V6aLvrSKvO+ZMWUzZstIucfuay+x7dJzkhGaDbGmoFPz4o5SVMwzJYVu/HrX8RzZcrMcTT0jpyOHDpST0F19IRd0PP4TatbP5h3kQeh46jtbQHLxNx7/+kiJXL76Y7jtk5UppODFmDJQoYfo1LdUwM8vviT96i945JCUp9VjNbUqBShj5WrZ+N8sarl2rVEiILMvnz1dKKXX+vFJTpihVp44cLl5cqeefV2rPnmz+AR6OnoeOozU0B2/SMTlZqSZNlLrtNqUuXUr3xLRpSgUGKpWQ4JTrWrlFryvZaW5J3rzw3NxgFtITPvpQ8kLNYs8e6NxZOrhERqI+m8PW2x/nmWdktT5kCBQpIn3XT56ULq+BgeZdXqPR5E7mz4dt26R7dPHi6Z547jn4+28oUMBVQzMNbeDtUKFCBVcPwW1o0QK2PzgeEhO59OKYLP+eXQ2VkqLvmzYR/+ZEZo0Io9HHfbizZT6++06qSYWGyjZanz5i6HMjeh46jtbQHLxFx/PnpWJnixbpOr9evSrb80o5tQKWlRpqA2+HqlWrunoIbsVL06szs8Bgii2ei9q1O0u/c5OGFy9Ksffz58Ew2PfyFwzpHE7Z90Yy4MVC5MkjVedOnpR/GzXK+Ly5CT0PHUdraA7eouOrr8KFC1LR0lbZZfJkKW/tYHnuW2GlhtrA2yEtZ1Mj3HYb5H/jNZbwMKvWZa3loU3DhASJhqteHfXOO/zx8q8EB0Ng9zv47PvS9Owp7ZZ37IABA26oNpvr0fPQcbSG5uANOm7bBjNnwuDBUtkakCJaEyZA167Q0LkV0q3UUBt4TbboO6IU4xsspu+kWmS5Z8I330iI6ksvsb9IMC0LhdLms14kJMAnn8hqfc4c6eiW5TqJGo1Gk02Sk+F//4MKFeDNN9M9MWYMxMeLQ96L0AbeDoUKFXL1ENyOfPng008hT+Qxtt87Unqv2sEwCnP8/UUcOl2Sdqwi6PSvVH+4IZs2SfO3QYOckoniVeh56DhaQ3PwdB1nzZJdwg8+SBdYt2cPfPaZfBnVrOn0MVipoaFUzmqMuyPBwcEqN/Qrdgemt/2e51Y/zMkxM/B7c8CNT+7dS8zg0Uyp9D4f/FgTdfEit9UszoD/5eGpp6B0adeMWaPR5F5On4ZatSSuZ/XqdLuFK1dKAfqNGz3yy8kwjB1KqQxb3ekVvB3Sah5rbqbHgofYnO8uCk4Yg7okxW/ij5zk8N19Sa7XgOS16wn95gBNmpzlhzUl2X8gDy++6JGfH5ej56HjaA3NwZN1HDkSrlyBadP+4wrs2FGaylj05WSlhtrA2yGtCYLmZsqWM4ge+QGlE0/z54MTWHHXOyRXr0HldV/wRYkX+HJsODP/7cLo0Xu5+27tW3cEPQ8dR2toDp6q49q18Pnn8NJLcMcdqQeTkuRgYiLksc4UWqmh85L9NF7NtWuQULcxC/I/Sbe1H/INj7OrYmfU+Hd4+onqVn5eNBqNJlMSEmDgQKhaVVpI2/jsM4m4K1cOOnVy2ficifbB2yEhIeHW7fhyGUeOwOxZiqMzVvLqpRFMLTmGRhfXsrPrWGb+4HvT67WGjqM1dBytoTl4oo5vvQVjx8Ivv8huPCD1OGrUgLp1pe68hVuMZmuoffA5JCIiwtVDcAuSkmDpUvlwPFR9F/dO7MCCS52o6pfAp9+W4fio6cxa5svatTf/rtbQcbSGjqM1NAdP0/HQIXjnHejZM51xB7H6585JcRuL/YdWaqgNvB1OnTrl6iG4lOPHJT20ShXo1g0e3DSCnTTi3pI7YPJkikTsI2/7Nrz+OnTx2865rn1IuJp8wzlyu4ZmoDV0HK2hOXiSjkrJDnzBgvDRR+meOHBAmlr06+f0ojYZYaWG2sBrbiA5GX76CR54QHxWH467QlCDZJYtg34TqmEMG0reI4elC0xqM4bCheHNJw/z8OV5rHp0rov/Ao1Go4GvvoI//pACdbfdlu6Jq1eheXMYN85lY7MK7YO3w+nTpylfvrxp53NnoqKkmtxnn8GJE3CbbwofN/mCrtteJd+Et+GZZ+yfQCkOlL+L0mcPcWlbGDWCpXpNbtLQWWgNHUdraA6eouO5cxItX60abN5saZD8LTFbQ+2DzyHJycm3fpEHk5ICv/4KDz0k2/Bjx0pF2XVvriXqtmAe/ukZ8lWpmC6vxA6GQdkvJ1OeM+x46G1SUuSwt2toBVpDx9EamoOn6DhsmDSTmT07nXG/dg3eflsC7FyIlRpqA2+HgwcPunoITiE6Gt59F26/XQJPNm2C4cPh8GFYFTiMVmPvwTh3Fr7+GrZsgTvvzNJ5y3YM5mDLZ3k4cjKL3z4AeK+GVqI1dBytoTl4go6//gpffAGjRkG9eume+OQTyZPbtMllYwNrNdR58LmElBQp9jBzpkTEJybC3XdLhOlDd1/Ap1gBabjerh2ULQtDh0IOaibX/O4dZjWvyWvvV+Wuvqb/GRqNRpMpsbHSjbJWLWkJayM6WrrLdOzotTnvGaFX8Hbw9/d39RAc5uxZeP992Xpv0wZ++w2efx7++QfWrErk0bOf4FP3dlnSg/RDHj06R8YdwPAtT5tfRxKb6MOgQeDn5/kauhpvmIeuRmtoDu6u45gxcOyYbM0XLJjuiZdfhrg4l6TF/RcrNdQG3g4VK1Z09RByhFL8v73zDo+i6v7450IwhCZFA1KkKEWKEYgiGoHXV7EC/sQCiggWFAv2gu1VbCi+dnztFQG7YkGsNAmh10gAqaEKQQgtkOT+/jibECDZlJ2d2dmcz/PkSXZndubkO3f2zL333HOYPBmuuAIaNIB77oH4eBm2WrdOyrK3WjleiiHfeqssFbnkEsfOf/zx8GG/n7jt627MSz7aseOWV/zaDiMJ1dAZIlnHGTPgpZcka90ZZxTYMG2apKS9807p2nuMmxqqgw+C3worbNsmDbxNG+jaFX74QYarFi2SQklXXRXomD/8sAxTZWfD11/DL79AQoKjtlx8eSW6MYm/Bj/Ptm2OHrrc4bd2GImohs4QqTru2wfXXSd13vMGI/OpWxeuvhoeesgT2w7FTQ11Dt7nWAvTp8vc+iefwN690KkTvPceXHaZrFEHZN1ITo505S+9FGrVkrH6wFp2p6l41r/Ydval3P7zCB654VpGfNo4LOdRFEV5+mkp6/7NN3DkkYdsPO44eP99L8zyHO3BB6FatWpem1AkO3bAa6/J6Pppp8EXX8CAATB3rjj8AQMCzn3/fnj5Zcm7fNdd8uETT5ThqjA59zxqvf0cFWIMp352Jz/+GNZTRTWR3A79gmroDJGo44IFsvqtb1/o2bPAhq1boV8/mZSPINzUUBPd+IxZs6S3PmaM1DZu317mnPr2herVD9l5/Hhx5EuWwFlnSb7Gtm1dtTd72FPE/OdBLjr6Dz5YdtrhT9eKoihlZP9+WcWbng6LF8sCoHwGD5Zou3nzXP/ecxNNdFNGkpOTvTYBkKUfb70FiYlw8skwejT06QMzZ8KcOTBoUCHO/cUXZZ49JwfGjYOffvKkkc/qejrLhn/Bt1s6c889rp8+KoiUduhnVENniDQdR4yQ78CRIw9x7jNnSk/ollsizrm7qaHOwQchKyvL0/PPmydt9OOPITNT2ukrr0iwXKE94X/+gYwMyc/Yp48sfg/jPHtJ2Gstze+7mLu2wisj9nDZZXGcdZZn5vgSr9thNKAaOkMk6bh4sSxtv+SSQxYB5eRI771ePakaF2G4qaH24COM3bslQO7UU2X4/f33JZXsH3/IXNMttxTi3HNy5EmgeXOZfAdp3C7Ms5eUx0/9nnUVj+Wx/n+Rmem1NYqi+JmcHLjmGhm5HDnykI1vvw2zZ8uUZI0antgXMVhro+anY8eO1kn279/v6PGCsXixtUOGWHvkkdaCta1aWfvii9Zu3VrMBydOtDYhQT7UpYu1c+e6Ym9JydcwPd1mx1Wz33G+vWFQrrdG+Qw322G0oho6Q6ToOHy4fOWNHl3IxowMa196ydrcyPyecVpDYJYtwidqDz4IS5cuDevx9+6VkoZnnCFr119/XabNJ02C1FSpyFq7dpADfP655Jvdtg0+/VRy0XpQ3zgY+Ro2aEDFJx7jAn5g05tfM2GCt3b5iXC3w/KAaugMkaDjwoWSse6SS2Qm8iByc2UJ8JAhnmesKwo3NVQHH4TNmzeH5bhpabJirWFDmU/fuBGefVYiQUePhi5dgrTN3btl8gngggskyuTPP2VtewQ26IM0vPVWctu247WYIQwZmKkJcEpIuNpheUI1dAavddy3D/r3h5o1ZZnwQV95v/0GHTvCypWe2VcS3NRQHbxL7NsniWjOPFPywr/8MvzrX5JELi1N0skeHSyrq7XSSz/hBOjRQ9aHxMVJGbj8bDYRTqVKVHjzDerlrqfDxh8YMsRrgxRF8RNPPHEg+Pig78usLLj5ZkkQUq+eZ/ZFGurgg9DWgeUVf/0lZQsbNpThpJUrpYLb2rXw2WdSAKZCcVdhwQJ5Grj8chl+ev99qFQpZNvc4DANO3fGLFtGi4cvZ9Qo+PJLb+zyE060w/KOaugMXuo4c6Z8d/bvDxdddMjGZ5+VfB+vvlrmQllu4aaGukwuCGVdzrB/P3z7rcyp//wzVKwone4bboDu3Uvg0Asyc6aE1NeqBf/7H1x/vRzQJxSqYbNmPPQQpH06n5tuaEtSUkXi4923zS9E0tIkv6IaOoNXOu7ZI+nkjzlG6m0cRFqadO0vv1yqYUY4ukwuQli2bFmp9l+9WuoZHHss9O4tU+OPPSbvf/WVlCIukXPPyZFeO0h2m2eegaVLJWWdj5w7FK1hpUVzGbO0A1dsG8mgQTIDoRROaduhcjiqoTN4pePQofJ9+s47Mv9+EM8/L9OUL77oiW2lxU0N1cGHSHa2JIq74AJo2lSGkDp2lB78qlUS7Vmq8r9//CHp6pKSpJi7MTLPHjSc3oecdBKme3eGxzzIrG/Sefddrw1SFCUS+fln6bXfequMgB7Gq69KgJ3OvR+GOvggNGrUqMht69ZJ77xpU+jVS4q8PPigzLF/9x1ceGEpO9sbNkhIfVIS/P235KatUyf0f8JjitTQGHjtNSpVyGHsUbdy222wfLm7tvmFYO1QKRmqoTO4rWNGhuTuOuEEGcg8bOP27RKP1L69q3aFgpsahtXBG2PONcakGWOWG2PuL2T7AGPM38aYeYGf6wpsu9oYsyzwc3U47SyKunXrHvQ6Nxd+/FECPBo3hkcfhdatJVBs9Wp4/HF5v9Rs3AgtW0qU/AMPSLDI5ZdH5LK30nKohgfRtCnm0UdJ2vI1F9sv6NdPRkSUgwmqoVIiVENncFNHa2VW8u+/JV/IYbFzt90meT/27nXNJidwU8OwOXhjTEVgJHAe0Broa4xpXciun1hrTwr8vB34bG3gP0An4BTgP8aYWuGytSjyKtNt3Cj1ho87TmI4pk2TUfPly2HCBEklW6ag9ryEB3k5kxcvhiefhKpVnfsnPKbY6n533gmnncZtfTaRkiL/vnIw0V4h0Q1UQ2dwU8dRo2Sl0bBh0KHDIRsnTJAd+vWDypVds8kJ3NQwnFH0pwDLrbUrAIwxY4FeQGoJPnsO8LO1NiPw2Z+Bc4ExYbL1MHJzYfbsmowcCV9/LT3Lbt1g+HBx6CGleF+9Whzb11/Los527eD2250y3V/ExMCUKXSsUIF++2QU5JxzZOGAoijlk1WrZFn7GWdweBXKzEwpodmqlcyLKkUSTgffAFhb4HU60iM/lN7GmC7AUuAOa+3aIj5bmlC1kNi/X0Z+UlNPonZtyXo4aJCMoofE3r3w3HMSiQfyaNq8ecj2RjI1SlLsIbC04PXuX5L1Ux369evK3LmFlMAtp5RIQyUoqqEzuKFjdjZccYXMUH74YSGxTA8+KIlEpkzxXe8d3G2L4XTwhU0gH7oY6ltgjLU2yxhzI/ABcGYJPysnMWYQMAigfv36TJw4EYBmzZpRvXp15s+fD0CdOnVo06YNkydPBiAmJoakpCTmzJnDjh07AEhMTGTTpk2sXbuWU05pzA03HMnFF+eyfPkiNmyAnJx4WrRowdSpUwGIjY2lc+fOzJo1i507dwLQqVMn0tPTWbduHQAtW7akYsWK/LlwIYnXXUfVVavIufhiZlx2GVl16xI3fz6dOnUiJSWFPXv2ANC5c2dWrlzJxo0bAWjdujU5OTmkpaUB0KBBAxo2bEhKSgoA1apVIzExkeTk5Pw1lklJSSxdujQ/LWLbtm3JysrKX6LRqFEj6tatmz9cVKNGDTp06MDUqVPJDkyEd+nShcWLF7N161YAEhISyMzMZMWKFQA0adKE2rVrM2fOHABq1apFQkICkyZNwlqLMYauXbsyf/58tgXy0nbo0IGMjAxWrVp10HVaMGsWJz9wG29bQ/0Vf9K7dyYPPLCk2OsE0Lx5c2JjY1m0aBEA8fFlv06pqTLAVK9ePZo2bZpfuzkuLs6z69SiRYv8du31dSrr/eT1dapSpUq+huXhfgrndUpOTg7r/fTUU5VJTm7C669vJzZ2NxMnFrhOxxzD7hkz2HPRRWyMjSUxYI/frtOaNWscu05BKaoKTag/QGdgQoHXQ4GhQfavCGwP/N0XeKPAtjeAvsWd0+lqclOmTAn9IBs2HPj7tdes/emn0I/pI0ql4cSJ1oL947S7LVj74Yfhs8tPONIOyzmqoTOEW8eJE62tUMHaAQOC7JSba21WVljtCCdOa4hH1eRmAs2NMU2NMUcAfYBxBXcwxhxT4GVP4M/A3xOA7saYWoHguu6B91wlO5SQ7r17D6yj+/57eW/wYDj7bGeM8wml0rBrVxg0iM7Tn2dQQgo33aRL5yDEdqgAqqFThFPHjAyJmTvuOHjllUJ2GDVK4peMCTEIylvcbIthc/DW2mzgFsQx/wl8aq1dbIwZZozpGdhtiDFmsTFmPjAEGBD4bAbwOPKQMBMYFnjPH3z/vdR/ffRRWSSfkOC1Rf5hxAhMgwa8umsgVWL20bevFOpRFCV6sRauuw42bYIxY6BatUN2mDcPBg7UZTalxNgoyhGamJhonVyCkJubS4VSJY5HWuk770iE58iRUj6uHFMmDX/5BTZv5svYvvS+xHDPPVJLorxSJg2Vg1ANnSFcOv7vf3DTTRKDfNddh2zct0+ye27eLEuJfZ7V02kNjTGzrbWJhW3TFh+ExXl114sjK0vyx4NUfXvmGZg/v9w7dyiFhgU56yy44gou7m0YfEMuI0bA+PHO2+YXyqShchCqoTOEQ8d58+COO6RWxx13FLLDE09IbY433/S9cwd326I6+CDkRVEGZcIEaNtWHkEBrrwS7r3X13NETlIiDYvigw94ddapdGybxVVXQXq6c3b5iZA0VADV0Cmc1nHHDrjsMsnK/eGHhRTjmjv3QI3YHj0cPbdXuNkW1cGXlbVr4ZJL5LHTGEmWrDhLfDwVZs9kwunDyMqCPn00la2iRAvWSgntv/6CsWPh6KML2alZM0kC5pNKcZGGOvggJBQVHPfhh+LQf/hBgj4WLoR//9td43xCkRqWhPPOg4EDqfPWcL64N4U//oCHH3bONr8QkoYKoBo6hZM6vvWWOPbHH5eMdYdhLRx5pEzM13I9U3nYcLMtqoMPQmZm5sFv5AUkNmok8+upqVIcprhkA+WYwzQsLS+8AA0a0H1Uf265ZjfDh8tzVXkiZA0V1dAhnNJx/nzJENq9O9x/WBkyYOZMqbvtUf35cOJmW1QHH4S8zEVs3gxXXy0VZkAC6caNgyZNPLPNL+RrWFaOPBLefx+WL+e5s34kIUGm49asccQ8XxCyhopq6BBO6Lh9u8xu1q4NH31UyLz77t1SOvvvv4sYt/c3brZFdfDByMmB11+XJPRjxkRVlTdfceaZkJZGbN+L+ewzqRVw6aWyeEFRFP9grdR3X7VKKsXFxxey0/33Q1qaPNjXrOmugVGGOviiSE2l8113Sfa59u1lTGnYMK+t8h1NnBrlOP54AJpvSWbUyO3MmFHEkpooxDENyzGqoTOEquNzz0kRzREj4PTTC9nh558ljd2QIVEb1+RmW1QHXxQVK1IpI0PSI/76q0bJl5HaTq5bXb8eunWjx89DuOceWZn40UfOHT5ScVTDcopq6Ayh6DhxonTOL70UbrutiJ1efVWShA0fXubzRDputkV18EXRsiVT3ntP1rWbworbKSUhr+KSI9SvLzkGPvyQpxO/oGtXWWazYIFzp4hEHNWwnKIaOkNZdVy/Xpa5Nm8uiT6L/Er97DP48UeIiyu7kRGOm21RHXwQ7GGFiBXPeeQRSEyk4o3X8+kL66hZE3r3hn/+8dowRVEKY98+SWazcyd8+SVUr17ITn/8Adu2SYKwxo1dtzFaUQcfhFpRtPbSKxzXsFIl+PhjyMoi/p6r+eyTXFatkoGWvGzB0Ya2w9BRDZ2hLDredpv473ffhdatC9lhzRq44AKp41EOcLMtarEZxZ+8954so7n7bl5/swKDB8ODD0raakVRIoO334brr4f77itiWj0nR5Ydz50rSemPO851G/2OFpspI5MmTfLaBN8TNg0HDpT5+AoVuGGQ5frrJang55+H53Reou0wdFRDZyiNjsnJcPPNcM4l7UPiAAAgAElEQVQ5Qaq8Dh8OU6ZI5c1y4tzdbIvq4IMQTaMbXhF2DcePx/yrG688s5vOnWWN7cKF4T2l22g7DB3V0BlKquP69RIb07AhjB4NhYYzpaTAf/4j0XdXXeWsoRGMm21RHXwQjEbPh0zYNYyJgcmTiR16J198ATVqwEUXQTQVD9N2GDqqoTOURMe9e8W579gha96LXBV27LHQr5+sdy1H18fNtqhz8Ir/uf9+eOYZ+PRTpje6lK5dJYnGhAkSk6coijtYK1m9P/pIpst69y5iJ2sLyVGrlAWdgy8j8+fP99oE3+OKho8/Dp06wfXXc2rdlbz9Nvz+O9x664H6QH5G22HoqIbOUJyOzz4rzn3YsCKcO8AHH0iWum3bnDfQB7jZFtXBB2FbOW2ATuKKhpUqSa0AgFGjuOoq6dS/8YYkxvI72g5DRzV0hmA6jhsHQ4fKlPpDDxWx059/SuQdyHxaOcTNthjj2pkUJZw0bSrLbAJJMp58Ur5Lbr8dWrSQSF5FUcLDggVwxRVS4fXdd4uYUt+zRzLeVK0quSw0kVjY0Tn4IOzYsYMa5fQp0yk80XDJEti0iZ0du5KUJJWrkpP9W05A22HoqIbOUJiOGzfCqadKlceZMyWjdKHceKMMq40fD+eeG35jIxSn26LOwZeRjIwMr03wPa5raC1cey307k21bWsZNw4qV4bzz4dNm9w1xSm0HYaOaugMh+q4axf06CE5p8aNC+Lc//lHol7vu69cO3dwty2qgw/CqlWrvDbB97iuoTEyRpiVBX36cOwx+/n2W3HuPXvC7t3umuME2g5DRzV0hoI65uRIiug5c2DsWBmeL5KaNSVb3eOPh93GSMfNtqgOXok+WraUHJnTpsH993PyyRKDN3OmLLuN1pz1iuIm99wD33wDL74ovfhC2bVLAmKyssTJ67pVV1EHH4RmzZp5bYLv8UzDyy+XaN3nn4effqJXL3jhBfjqK/li8hPaDkNHNXSGPB1HjpT76bbbZDlqoVgLgwfDww/L07UCuNsWNYo+CNULrWuolAZPNfzvfyVXZpcugHwZrVghX0xNmsCQId6ZVhq0HYaOaugM1atX56uv5N7p2VNusSJ56y1ZFP/oo5CU5JaJEY+bbVF78EHQ5Bih46mGsbGyIL5yZQny2b2b55+XVLa33w6ffuqdaaVB22HoqIbOMGrUKvr2hVNOkWmvIle6zZ4tXftzzpEevJKPJrpRFCfZsyc/013FCpbRoyWV7VVXwW+/eW2coviDxYvhgQfa0qQJfPstVKlSxI65uVL1KT4eRo3SlLQeosoHoU6dOl6b4HsiQsO4OOjfX8pavfIKcXGypKd5c+nNz5vntYHBiQgNfY5qGBpr18rqtsqV4ccf4aijguxcoYKE1X/9dTE7lk/cbIua6CYIubm5VNCnz5CIGA1zc+Hii+H77+HXX6FLF9LT4bTTJEHHH39ApMZhRYyGPkY1LDsZGRLGsnYtTJyYS/v2QXRMTpasN+WoOlxpcbotaqKbMjJ58mSvTfA9EaNhhQpS5KJZM7j0Uli3joYNJfdGVhZ07w4bNnhtZOFEjIY+RjUsG5mZcN55sGyZrEDZvj2Ijp98Ik/Mo0e7Z6APcbMtqoNXyg9HHinfUu3bSx15JH3tDz9Ius3u3aW3oiiK1HW/6CKJl/v0UzjzzCA7L1gA11wjwS2XXuqajUpw1MEHISZGVxGGSsRp2Lq1TCLWrQvZ2WAtp54qCTuWLpXeSmam10YeTMRp6ENUw9KRnQ19+0oQ6rvvQq9e8n6hOm7ZIk8CNWtKEfgjjnDXWJ/hZlvUOXilfLJrl3wpnXMO3H03IDFBl1wCXbvKVH3lyh7bqCgekJsLAwfChx/Cyy8HSWQDkszmrLMkiGXSJFmtoriKzsGXkTlz5nhtgu+JWA2rVIFateDee2WMHvH3770nvZbLLoN9+zy2MUDEaugjVMOSYS3cdJM498ceO9y5H6ajMZJB6r331LmXEDfbojr4IOzYscNrE3xPxGpojHwpJSTIWOSffwKyNv6112Sdb58+EmHvNRGroY9QDYvHWvHVb7whRd8Ky09zkI55Uak9e8o9pJQIN9tiWB28MeZcY0yaMWa5Meb+QrbfaYxJNcYsMMb8aoxpXGBbjjFmXuBnXDjtVMopVavK5HvlylItY8sWQNJnv/SSxONdeaXMRypKNGOtDGa98grccQc8/XQxK91+/hmaNpW5LCVysdaG5QeoCPwFNAOOAOYDrQ/Z519AlcDfg4FPCmzbWdpzduzY0TpJZmamo8crj/hCw2nTrD3+eGsXLz7o7eeesxasveIKa7OzPbLN+kTDCEc1LJrcXGsfeEDa+s03y+uiyMzMtDY11dqaNa1t187aHTvcMzRKcLotArNsET4xnD34U4Dl1toV1tp9wFig1yEPF79ba/MqdE8HGobRnlKzadMmr03wPb7QsHNnGaJv3VpeBwJP77oLnnpKlvUOHOhdmVlfaBjhqIaFYy088oi08+uvl6C6YD33LYsXw/nnS52Hb78FLeJTatxsi+F08A2AtQVepwfeK4prgfEFXlc2xswyxkw3xlwUDgOLY+3atcXvpATFNxrGxIgHHzIEHn88/+2hQ+XlRx/J/LwXw/W+0TCCUQ0Px1qpxfTEE3DddfD668Wkjc/KovbAgbBpkzj3xo2D7KwUhZttMZwL8gp7Dix0TZ4xph+QCHQt8Pax1tr1xphmwG/GmIXW2r8K+ewgYBBA/fr1mThxIiA1d6tXr55fuadOnTq0adMmP4tQTEwMSUlJzJkzJz/oITExkU2bNuVfgH379rFlyxYWLVoEQHx8PC1atGDq1KkAxMbG0rlzZ2bNmsXOnTsB6NSpE+np6axbtw6Ali1bUrFiRVJTUwGoV68eTZs2JTk5GYC4uDg6depESkoKe/bsAaBz586sXLmSjRs3AtC6dWtycnJIS0sDoEGDBjRs2JCUlBQAqlWrRmJiIsnJyWRlZQGQlJTE0qVL2bx5MwBt27YlKyuLZcuWAdCoUSPq1q1L3rLCGjVq0KFDB6ZOnUp2wIt16dKFxYsXs3XrVgASEhLIzMxkxYoVADRp0oTatWvnR4XWqlWLhIQEJk2ahLUWE+gKzJ8/n23btgHQoUMHMjIyWLVqlWPXqXnz5sTGxoZ+nSpUIHfpUuq98grr4+KoM2QIycnJJCXBTTc147XXjmXTpq0MHbqImBjr2nXKycnJb9fhuk5du3b1z3Uqw/20d+/efA39fD85dZ1mz57D8OFH8/nnjbjuuv3ce+8aJk8u5jpNm0a9jh3Z0L8/LU8+Wb/3yniddu7cyZo1axy7n4JS1Nh9qD9AZ2BCgddDgaGF7HcW8CcQH+RY7wOXFHdOp+fg09PTHT1eecR3GmZlWdu1q7WVKln7668Hbfrvf2Weslcva/fudc8k32kYgaiGB8jNtfbWW6UtDxkSfM49n7//ttaqjk7gtIZ4NAc/E2hujGlqjDkC6AMcFA1vjGkPvAH0tNZuLvB+LWNMbODvo4DTgdQw2looxT4dKcXiOw2POELC51u0gP/7P0nBGeDOOyXK+JtvpG5NoOMRdnynYQSiGgo5OTBokLTju+6CF18sQV2YV16R+2HJEtXRAdzUMGwO3lqbDdwCTEB66J9aaxcbY4YZY3oGdhsBVAM+O2Q53AnALGPMfOB3YLi11nUHnzdEpZQdX2pYqxaMHy/1rFevPmjTLbfIOuHx46V85vbt4TfHlxpGGKqhJG7q2xfefhseeghGjCiBc//sM1kc37UrNG+uOjqAmxqGNSmutfYH4IdD3nukwN9nFfG5aUC7cNqmKEFp1AhSU6FSJXmdm5sfgTRoENSoIUF3Z54pqe2PPtpDWxWlGHbvht69pa0+95z03otl4kTo1+9AhbiKFcNtpuIwmskuCPHx8V6b4Ht8rWGec3//fcm3vXt3/qY+fWSoPjUVzjhDamWHC19rGCGUZw23b5eSCz/9JL33Ejn3JUukwszxx8O4cRAXB5RvHZ3CTQ3VwQehRYsWXpvge6JCwypVpDdz2WUH5a49/3z50tywQapkpoZpEikqNPSY8qrhunXyAJqSAmPHwrXXlvCDTZpI7/3HH6F27fy3y6uOTuKmhurgg5C3fEcpO1Gh4WWXySLh77+HAQNkuD7AGWdIEa39+8XJT5ni/OmjQkOPKY8apqZKDqeVK6XplqhM+4YNkJEh6ZtHjpSpqgKURx2dxk0N1cErSkkYNEgSdI8eLclwCpRZPukkSE6WEvNnny0lsRXFS6ZOlQfO/fth8mRpl8WydavseOGFB7Vvxb+ogw+CLgkJnajS8L774J57JLr+EJo0kZLYHTtKh//FF507bVRp6BHlScMvvpCQkfh4mDYN2rcvwYcyM+G882D5cnjyySLD68uTjuHCTQ2NjaIntcTERJuXoUhRwoK1B778/v77sPD5PXukAt1XX8mSuhdekCy4ihJurIXhw+GBB2Roftw4OOqoEnxwzx4JKJkyRRpujx5ht1VxDmPMbGttYmHbtAcfBH1YCJ2o0zDPuS9fDq1aiQcvQFycLB2+6y549VUZ7Qx1rXzUaegB0a7hvn1wzTXi3Pv2hd9+K6FzB6kPO2kSfPBBsc492nV0Azc1VAcfhLw8y0rZiVoNmzSBf/9b0tuNHHnQpooVZa3xm2/Cr7/KMuKVK8t+qqjV0EWiWcOtW6F7d1nN+eij8PHHEiNXYh5+WELsr7yy2F2jWUe3cFNDdfCKUhZiYuSbtFcvGYt/663Ddrn+epgwAdavh1NOkU6SojjJwoXStqZPl/jP//ynBNnpQLr8L70k5REbNJDAESXq0Dn4IOzZs4e4QIIHpWxEvYZZWZKzfvx4+OEHCVQ6hKVLoWdP+OsveP55eR4o0ZdwgKjX0AWiUcPPPpNVm0ceCV9+CaeeWsIP7t8Pl18u8+0//ihZcEpINOroNk5rqHPwZSQ9Pd1rE3xP1GsYGyvfrsOGSd7aQmjRAmbMkDimIUPkS7k0hWqiXkMXiCYNc3Jg6FDpdCckwOzZpXDuWVmyIP6rr6QHXwrnDtGlo1e4qaE6+CDk1TZWyk650LByZZnHjI2VJCEff3zYLjVqyHfqY4/Bhx9KgpxACe9iKRcahplo0XDLFrjgAomWHzQIfv8djjmmhB/es0dGm775Bl5+WZ42S0m06OglbmqoDl5RnGTECEnx+fzzh22qUAEeeUSWLy1bJuuTx40r5BiKUgjTpklSpd9/l4qGb7whz5QlZskSWQr3xhtw661hs1OJHNTBB6Fly5Zem+B7yp2Gw4bBJZfIOrlhwwrNCNajB8yZA82aSYze3XcflOL+MMqdhmHAzxpaC//9r1RsjY2VrImDBpXiANnZ8rt9ewkEKdWHD8bPOkYKbmqoDj4IFbU8YsiUOw0rVYIxY6B/fwlpvvPOg3LX53HccZL57qab5Mu7S5fDSs/nU+40DAN+1XDLFrjoInkI7NFD5ts7dCjlAU4/Hf73P3kdYiUzv+oYSbipoTr4IKSGqzxYOaJcahgTA++9J3Oc48bBtm2F7pZXz2PsWFi8WAKmxow5fL9yqaHD+FHDX36BE0+UBRrPPy8paGvWLMUB1q6VYI/582UpnAP4UcdIw00N1cErSjioUEES0s+cCXXqyBj83r2F7nr55TBvHrRuDVdcIVP4oWa/U/xLVpb02M8+W5bAzZghyeZKs7SSJUuk575+vSRj6NkzbPYqkYs6+CDUq1fPaxN8T7nW0JgDtbQHD5YlSUX05ps1k6pfjz4qPfqEhAOJccq1hg7hFw0XLJAlb//9rzSZ2bMlsK5UZGTInE9WFkycKJP3DuEXHSMZNzVUBx+Epk2bem2C71ENA5x5pkRHJSXJ0GkhxMTItP2UKZLutls3GeWPj1cNQyXS2+H+/fDEE5CYKJ3ub76B116DKlXKcLDateVgU6eWsJRcyYl0Hf2Amxqqgw9CcnKy1yb4HtUwwBVXyFBperp00RYsKHLXzp1l8623wiuvQLt2uZrmNkQiuR0uWiTX/OGHoXdvicco04j6668fGPYZNAiaN3fUTohsHf2Cmxqqg1cUt/jXv6RXZYxkK8nKKnLXqlUlF8nEifK6WzeJuP/nH1csVVxg715x6h06wJo18PnnEmRZ4ipweeTmwn33yZj+22+HxVbFn6iDD4LmXA4d1fAQ2rWTyiAffiiLmoupBdG1K3z88SJuv13yk5xwgszRR1EJCVeItHb466/SFJ54Avr0kV57795lONDu3VIf9tlnxcG/957jthYk0nT0I25qqMVmFMVLXnpJqtG89JJMwgdh9my48UaYNUvKg44cCccf75KdiiNs3Aj33gsffSTX7vXXpepwmfjnH/nw3Lni4O+6q5Sh9ko0oMVmykhKSorXJvge1bAYNmyQaKrzzy9y/D1Pw44dpfP/8ssSr9emjYzM7tjhpsH+xOt2mJUlPrh5c/jkExmaX7gwBOcOUuDgpJMk18Ldd7vi3L3WMRpwU0N18EHYU5qSX0qhqIbFMHw4vPOOTLYnJkrE1SEU1LBiRQm+S0uTuL1nn5Vqde++W2jCPCWAV+3QWvj2W2jbVh7GzjxThuOHDZNER2Vi9GhYsUJyLbzzDlx4oaM2B0Pv59BxU0N18IriNddcIw5+1y447TT4++9iP3LMMTLdOmOGrKG/9lp5PpgwQefnI4WUFImr7NlTMhhPmCDL38o8rbJvH9x8M1x5pTzZKUox6Bx8ELKysogtVbkm5VBUw1KwYYPkJ73qKnltLRhTrIbWSvT1gw9KCdpu3eDpp0tRI7wc4GY7XLJErsWXX0rq94cfhhtuECdfZtLTpY779OkyHP/008XGbIQDvZ9Dx2kNdQ6+jKxcudJrE3yPalgKjjnmgHP/+WeZoF23rlgNjZHh+rQ0WTefmirrqnv1ksA8xZ12mJYGV18tw/E//QSPPSbF2265JUTnvnChrKVbtAg++0xKEnvg3EHvZydwU0N18EHYuHGj1yb4HtWwjGzbJuPvJ53Evm++KdFHjjhCnMlff8kc7+TJMmx//vkSlFeeCWc7TE2VB6zWrcX/Dhki1+CRR6BaNQdOcNxxMtY/c6aUIvYQvZ9Dx00N1cErSiRy2WWyHq5+fU68/36J0ApWNL4A1arJsPDq1fDUU+IXTjtNAry+/16D8ZzAWkkp/H//Jz32vED2Vauk8luIVVkl803//pCZKflqP/kEWrVywnSlHKEOPgitW7f22gTfoxqGQKtWMH06e/r3l6CqL74o1cdr1IChQ8XpPPecLLe/8ELpab7+uuRIKS841Q7375dEQ6ecIvVcJk8+EPvwzDMOOHaAr7+W5W9ffRU0pbEX6P0cOm5qqA4+CDk5OV6b4HtUwxCJi+Of4cPh99+lriyIpy5FcGzVqpIDZeVK+PhjqF5dkp41bAi33y5DzNFOqO1w9Wp46CE49lhJHLd9u6QvWLsWHn+8DOllC2P7dhg4UIYFmjWTBDann+7AgZ1D7+fQcVNDdfBBSEtL89oE36Mahk5aWpqExhsjkfYnnyy57DdsKNVxKlWSueIZM2R4uXt3cVJt2sAZZ0j23J07w/M/eE1Z2mFWlkTCX3ABNG0q0x0dO8pw/JIl8pBUpmpvRTF4sFyEBx+EadMiMk2h3s+h46aG6uAVxU/Uqyee5vffZaz9/fdLvfDdGKlaO3YsrFsnQdmbNkkEeN268hDw3XclnvKPKnJzpSDb9deL1L17S0f6oYdkGP6776BHD8kx4wg7dx7Ie/DEE+LYn3hCIiYVJUTUwQehQYMGXpvge1TD0DlIQ2Mk2cm8eRLdNXAgnHtu0Mp0wTj6aAkOS0uTXn3//pKQpUcPcXADBsiU8K5dzvwvXhGsHe7bJ8vabroJGjWSwZIxY0SDCRMk3m3YMBmed5Tx42X4ZNAged2sGXTq5PBJnEXv59BxU0NNdBOEPXv2aPWkEFENQ6dIDXNzJVouNRVefVXeCyTHCYU8hzd2rPRYt2+XtKpnny3PEmedJTnV/VTX5FAN16yRVAM//ww//ij/Y5Uq8v/17i05BKpWDZMxmzdL8MOYMVIe8K23Im6uvSj0fg4dpzUMlugGa23YfoBzgTRgOXB/IdtjgU8C21OAJgW2DQ28nwacU5LzdezY0TrJ77//7ujxyiOqYeiUWMO5c61t397aKVMcO/e+fdb+8ou1t95qbePG1soThLUNG1o7YIC177xjbWqqtTk5jp3ScXJzrf3442Q7apS1gwdb26LFgf+jXj1rBw60dtw4a3fvdsGYSZOsrV3b2iOOsPaxx6zdu9eFkzqH3s+h47SGwCxbhE8MWzokY0xFYCRwNpAOzDTGjLPWFozZvRbYZq093hjTB3gGuNwY0xroA7QB6gO/GGNaWGs1hFNRiuKff2DrVomYy1taV7duSIesVEkS6v3731LR9q+/pJb5L79IsNn778t+NWvK6HKHDlLnvF07aNkyxAxuZSAnR+qwLFggP/PnS3bXTZskb2+1aiLPjTfKiESbNi6NROzeLUMEbdpA167w5JPSe1eUMBLOfIenAMuttSsAjDFjgV5AQQffC3g08PfnwKvGGBN4f6y1NgtYaYxZHjieq/m4qjmShqp8oxqGTok17NZNhuufekoi5776SvKl3nGHI3YYI4Hdxx8vudVzc2XF3vTpB35GjIDsbNm/UiVJwta0qUwvN20qc9zx8fJTty7UqlXygLXcXMn7snmz1FXftEkWEqxYIQ8ef/0lf+/dK/tXqCCV9s45B+rXX02fPo1p21Yq8rnGmjWSjGDJElm+UKeOhOb7FL2fQ8dNDcPp4BsAawu8TgcOjSDJ38dam22M2Q7UCbw//ZDPuh7dkZhY+LSGUnJUw9AplYZVq0rvsH9/yX6XV2M+b1TasfBvOVSrVvIzYIC8t2+fBOwtXCg/y5aJ0502Tea5C6NyZTG7atUDJVTzQoOysyXAb+fOohPzVKkiDxItWsgcetu2cOKJssjgwFRnY4f+6xKSkSEPWnmxEffcI/+Mq08XzqP3c+i4qWE4HXxhA1+HRvQVtU9JPisHMGYQMAigfv36TJw4EYBmzZpRvXp15s+fD0CdOnVo06YNkydPBiAmJoakpCTmzJnDjh07ABF+06ZNrF0rzyXWWtq1a8eiQI3u+Ph4WrRowdSpUwGIjY2lc+fOzJo1i52BBcSdOnUiPT2ddevWAdCyZUsqVqxIaiCbSL169WjatCnJgeTgcXFxdOrUiZSUlPw6wZ07d2blypX5OYtbt25NTk5O/vrJBg0a0LBhQ1JSUgB5IkxMTCQ5OZmsQDR1UlISS5cuZfPmzQC0bduWrKwsli1bBkCjRo2oW7cueUGJNWrUoEOHDkydOpXsQBesS5cuLF68mK1btwKQkJBAZmYmK1asAKBJkybUrl2bOXPmAFCrVi0SEhKYNGkS1lqMMRxxxBFUqVKFbdu2AdChQwcyMjJYtWqVY9epefPmxMbGRu11Wr9+fX5yjBJfp7g4Vtx+O+Tm0mTVKuKTk8l9+GFWDhhA7oUXknDSSQddp65duzJ//vyQr9OiRXKd6teHnj0Pvk5HH92CjIw4Zs5czbZtR5CdXZuYmHiWLUtnz56KZGdXonr1eLZt25r/v8XH1wF2kZu7g7i4HBo3rk18POzZs5JatfbRqlVNOnY8lunTD79OKSkHrtOUKVOICRRoCff9VG3ZMk66805idu1i07nnsuLqq8mqW5culSqxeOHCkO4np65TWe+nXbt2Ubt2bV/fT2587wW7Trt27aJdu3aOXaegFDU5H+oP0BmYUOD1UGDoIftMADoH/o4BtiDO/aB9C+4X7EeD7CIP1TB0HNHwu++sPe446cd37Civc3NDP65PCHs73LnT2tmz5e+sLGuvv97aBQvCe04P0Ps5dNwMsgvnOviZQHNjTFNjzBFI0Ny4Q/YZB1wd+PsS4LeAweOAPsaYWGNMU6A5MCOMtipKdHPBBTIP/O67Mnx84YXQp4/XVvmfHTskCX2TJrJwft8+SVLz5psSaagoHhLWdfDGmPOBF4GKwLvW2ieNMcOQJ45xxpjKwEdAeyAD6GMPBOU9CFwDZAO3W2vHF3c+p9fBZ2dn5w/rKWVDNQwdxzXcvx8++EBCyvv0gT17YNQo6Nev4KR1VOG4hunp8MILsoY9M1Mm/x9+WMr2RTF6P4eO0xoGWweviW6CkJqaqtWTQkQ1DJ2wa/jJJ+Lojz4arrtOMqs1aRK+83mAIxpaKxkDK1eWTEDnny9lfe+8E8pJ8Jnez6HjtIbBHLymqg1CXqCGUnZUw9AJu4aXXQYTJ0rv85lnZE3bhRf6Pz9tAULSMCNDkgC0bSuFYEAW0a9aBaNHlxvnDno/O4GbGqqDV5TyjjGSfOXrr8VpPfSQrIHLy9U6apRUXImi0b4SMWGC1IatX19Sy1arJpX8QDRr2NBb+xSlGHSIPghbtmzhKEcKPZdfVMPQ8VTDrCzJSLN9u2Rhu/RSqVferp2vktGXSMP9+yVbT1KS/G8DBkgy/ssvl/JyJ53kiq2RjN7PoeO0hjpEX0ayylihSzmAahg6nmoYGysp4l57TfLRPvYYJCTAc8/J9v37y1zJzk2K1DAjQ4bZ+/aVGIQuXSS/Lcj/uGEDjBypzj2A3s+h46aG6uCDkJccQSk7qmHoeK5hnToweDBMnQrr18sSsJ49ZdtPP0m+2XPOEYc4Y4YsFYsw8jXcsQMCyUeYOBGOOgquvBJ++w0uuUSmKVq0kO1HHeV+Mv0Ix/O2GAW4qaGud1AUpeTUqyfD1Xk0biyvf/lF0rGC9Pr//FOSzy9dKsnhW7WS9eFukpMDM2fC7Nm0+vZbKfi+ZAk8/rgEy7VvD48+Ct27wymnOJrGV1EiAXXwQWjUqJHXJvge1TB0IlrDtm0lwhykdz9tGsyaJY4f4Pnn4Y03ICZGeusJ3H0AAAieSURBVMatWkni+Gefle1r1sh891FHlW0N/tq1sHq1/F67Vo7XogUMGSLbu3WDrCyZ8zz1VFkOeN55su3II+GRR0L698sbEd0WfYKbGmqQXRB27typ1ZNCRDUMHV9ruHy59KIXLZLqM8uXy/uBHOWcfz6MD+SwiouT4f62bSWCHWRN/uzZUqglb76/XTsZSgepSbt06YHzHXmkOPHXX5fXv/0GzZuzs2ZNqlWvHv7/N8rxdVuMEJzWMFiQnfbggzBr1iy6devmtRm+RjUMHV9rmFdftijuvVei8rduhS1bZH68Xr0D2+Pi4JhjZASgUiUZ5m/b9sD2556T9xo1kp9DnfiZZwIwa+JE/2oYQfi6LUYIbmqoDl5RFO/o1k1+iiJv+L8oevRw0hpFiSo0qiQINWrU8NoE36Maho5qGDqqoTOojqHjpoY6B68oiqIoPkUT3ZSRqVOnem2C71ENQ0c1DB3V0BlUx9BxU0N18EHIzs722gTfoxqGjmoYOqqhM6iOoeOmhurgFUVRFCUK0Tn4IOTm5lJBs1uFhGoYOqph6KiGzqA6ho7TGuocfBlZvHix1yb4HtUwdFTD0FENnUF1DB03NVQHH4StW7d6bYLvUQ1DRzUMHdXQGVTH0HFTQ3XwiqIoihKFqIMPQkJCgtcm+B7VMHRUw9BRDZ1BdQwdNzVUBx+EzMxMr03wPaph6KiGoaMaOoPqGDpuaqgOPggrVqzw2gTfoxqGjmoYOqqhM6iOoeOmhurgFUVRFCUKiap18MaYv4HVDh7yKGCLg8crj6iGoaMaho5q6AyqY+g4rWFja+3RhW2IKgfvNMaYWUUlEFBKhmoYOqph6KiGzqA6ho6bGuoQvaIoiqJEIergFUVRFCUKUQcfnDe9NiAKUA1DRzUMHdXQGVTH0HFNQ52DVxRFUZQoRHvwiqIoihKFqIMvBGPMucaYNGPMcmPM/V7b40eMMY2MMb8bY/40xiw2xtzmtU1+xRhT0Rgz1xjznde2+BFjTE1jzOfGmCWB9tjZa5v8hjHmjsB9vMgYM8YYU9lrmyIdY8y7xpjNxphFBd6rbYz52RizLPC7VjhtUAd/CMaYisBI4DygNdDXGNPaW6t8STZwl7X2BOBU4GbVsczcBvzptRE+5iXgR2ttKyAB1bJUGGMaAEOARGttW6Ai0Mdbq3zB+8C5h7x3P/CrtbY58GvgddhQB384pwDLrbUrrLX7gLFAL49t8h3W2g3W2jmBvzORL9UG3lrlP4wxDYELgLe9tsWPGGNqAF2AdwCstfustf94a5UviQHijDExQBVgvcf2RDzW2slAxiFv9wI+CPz9AXBROG1QB384DYC1BV6no44pJIwxTYD2QIq3lviSF4F7gVyvDfEpzYC/gfcC0xxvG2Oqem2Un7DWrgOeA9YAG4Dt1tqfvLXKt9S11m4A6QQB8eE8mTr4wzGFvKdLDcqIMaYa8AVwu7V2h9f2+AljzIXAZmvtbK9t8TExQAfgf9ba9sAuwjwsGm0E5ol7AU2B+kBVY0w/b61SSoI6+MNJBxoVeN0QHY4qE8aYSohz/9ha+6XX9viQ04GexphVyFTRmcaYUd6a5DvSgXRrbd7o0eeIw1dKzlnASmvt39ba/cCXwGke2+RXNhljjgEI/N4czpOpgz+cmUBzY0xTY8wRSDDJOI9t8h3GGIPMe/5prX3ea3v8iLV2qLW2obW2CdIOf7PWas+pFFhrNwJrjTEtA2/9G0j10CQ/sgY41RhTJXBf/xsNVCwr44CrA39fDXwTzpPFhPPgfsRam22MuQWYgESLvmutXeyxWX7kdOAqYKExZl7gvQestT94aJNSPrkV+DjwwL4CGOixPb7CWptijPkcmIOsjpmLZrQrFmPMGKAbcJQxJh34DzAc+NQYcy3y4HRpWG3QTHaKoiiKEn3oEL2iKIqiRCHq4BVFURQlClEHryiKoihRiDp4RVEURYlC1MEriqIoShSiDl5RooxA9bSbCryuH1jmFHEYY34wxtQsxf6PGmPuDqdNihItqINXlOijJpDv4K216621l3hoz2EYoYK19nwt/qIo4UEdvKJEH8OB44wx84wxI4wxTfJqUhtjBhhjvjbGfGuMWWmMucUYc2egEMt0Y0ztwH7HGWN+NMbMNsZMMca0OvQkgd70R8aY3wL1ra8vsO0eY8xMY8wCY8xjgfeaBOqxv4YkTWlkjFlljDkqsP3OQL3xRcaY2wsc60FjTJox5hegJYqilAjNZKco0cf9QFtr7UmQX82vIG2R6n6VgeXAfdba9saYF4D+SAW7N4EbrbXLjDGdgNeAMws514nAqUBVYK4x5vvA8ZsjpZcNMM4Y0wXJ3NUSGGitvSlgG4HfHZEMc50Cn0kxxkxCOiF9AvbGIA8GWnxHUUqAOnhFKX/8bq3NBDKNMduBbwPvLwRODFQAPA34LM8BA7FFHOsba+0eYI8x5nfEqScB3ZGUpgDVEIe/BlhtrZ1eyHGSgK+stbsAjDFfAmcgDv4ra+3uwPtaF0JRSog6eEUpf2QV+Du3wOtc5DuhAvBP3ghAMRya69oiPfCnrbVvFNwQGEnYVcRxCivTXNQ5FEUpAToHryjRRyZQvawfttbuAFYaYy6F/IC4hCJ272WMqWyMqYMU1piJFGq6JjASgDGmgTEmvpjTTgYuClQsqwr8HzAl8P7/GWPijDHVgR5l/b8UpbyhPXhFiTKstVuNMX8EAuvGAyPLcJgrgf8ZYx4CKiH16OcXst8M4HvgWOBxa+16YL0x5gQgOTDEvxPoB+QEsXmOMeb9wPEA3rbWzgUwxnwCzANWI05fUZQSoNXkFEUpE8aYR4Gd1trnvLZFUZTD0SF6RVEURYlCtAevKIqiKFGI9uAVRVEUJQpRB68oiqIoUYg6eEVRFEWJQtTBK4qiKEoUog5eURRFUaIQdfCKoiiKEoX8P+5lyhXBMj9VAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "##############################################\n",
    "#Making graphs to show actual trajectories vs the estimated ones\n",
    "\n",
    "t_orig = np.linspace(0,periods-1,100)\n",
    "t_cent = (t_orig - center)/center\n",
    "t_design = pd.DataFrame( zip(t_cent**0, t_cent, t_cent**2)).to_numpy()\n",
    "\n",
    "#Actual\n",
    "true_curve = t_design.dot(gc_np)\n",
    "\n",
    "#Estimated\n",
    "est_coef_np = model.coef.weight.detach().numpy()\n",
    "est_curve = t_design.dot(est_coef_np)\n",
    "\n",
    "#Now superimposing on a plot\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(t_orig, true_curve[:,0], color='b', label='True Curves')\n",
    "ax.plot(t_orig, true_curve[:,1], color='b')\n",
    "ax.plot(t_orig, est_curve[:,0], linestyle='--', color='r', label='Estimated Curves')\n",
    "ax.plot(t_orig, est_curve[:,1], linestyle='--', color='r')\n",
    "ax.set_xlabel('time period')\n",
    "ax.set_ylabel('y')\n",
    "ax.grid(True,linestyle='--')\n",
    "ax.legend(facecolor='white', framealpha=1)\n",
    "plt.show()\n",
    "##############################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this is a nice illustration, but there are a bunch of things I would need to update to make it viable in practice. Some of those things are:\n",
    "\n",
    " 1. A function to estimate the latent classes for out of sample data. \n",
    " 2. Different loss functions for count or binomial data.\n",
    " 3. The ability to estimate an equation predicting the latent class.\n",
    " 4. Standard errors/covariances for model coefficients.\n",
    " \n",
    "For 1, you want to do [cross-validation](https://www.tandfonline.com/doi/abs/10.1080/03610926.2012.719986) to figure out the correct number of classes. But that isn't directly possible in this code. I have an idea how to do it (use the coefficients to solve for the latent classes, then estimate the error), but I need to confer with others to see how they solve the problem.\n",
    "\n",
    "For 2 my work often uses counts over time, or it is common in criminology to be looking at binary 0/1 data as well. That you should just need to update the loss function and a link function for the `forward()` prediction. So I don't think that will be much work.\n",
    "\n",
    "For 3 we often want to estimate an equation that predicts the latent classes based on time-invariant characteristics. I'm not sure how to do that in this set up though.\n",
    "\n",
    "This is basically set up out of the box to model however complicated temporal functions you wanted though, and can also incorporate inconsistently measured panel data. So if you had one observation measured at time periods `[1,3,4,8]`, and another measured at `[2,3.4,6,9.1]` this code would just work out of the box. If you did not want to do polynomial time functions, but something [like splines](https://link.springer.com/article/10.1007/s40865-016-0025-6), you just feed it those functions. Also can incorporate whatever other time varying attributes as well.\n",
    "\n",
    "You might also want to have a model for balanced vs unbalanced data. For balanced data with array broadcasting, you really only need to fit the outcome vector into memory. So can potentially do this for giant datasets. \n",
    "\n",
    "For 4 I don't know if you can get standard errors for the coefficients using the gradients. I know you can do Bayesian MCMC [modelling in Tensorflow](https://blogs.rstudio.com/tensorflow/posts/2019-05-06-tadpoles-on-tensorflow/), so maybe that is a way to get posterior distributions for the coefficients.\n",
    "\n",
    "Also don't take this as the best way to fit these models. Not sure if there are potential optimizations internally into how I created the `forward()` function to do the prediction, nor of the optimizer (this is basically just me doing the predictions and seeing if the machine could figure it out!). So I have no clue if those are written in the best way possible. If you follow the link before about cross-validation, these models tend to be notorious to fit in practice, and so maybe doing a better job with the initial starting parameters would go along way. Pytorch ends up being a pretty simple & fast tool through to fit these equations.\n",
    "\n",
    "If you have suggestions always feel free to get in touch!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}