Multiple_Event_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook will now have \n",
    "\n",
    "1. Regression of Behavior + GCamp \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy\n",
    "import _pickle as cPickle\n",
    "\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import ListedColormap\n",
    "import matplotlib.ticker as ticker\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "import statsmodels.api as sm \n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "\n",
    "\n",
    "import seaborn as sns\n",
    "import os\n",
    "\n",
    "import rpy2.interactive as r\n",
    "import datetime\n",
    "\n",
    "# to render plots with white background\n",
    "import seaborn as sns\n",
    "sns.set_style('white')\n",
    "sns.set_context('notebook', font_scale=1.75)\n",
    "\n",
    "%load_ext jupyternotify\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython\n",
    "%R library(rms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -o f\n",
    "test_data <- data.frame(time = c(1,1,1,2,2,2,3,3,3,4,4,4),\n",
    "status = rep(c(1,0,2),4),\n",
    "n = c(12,3,2,6,2,4,2,0,2,3,3,5))\n",
    "# # fit a Kaplan-Meier and plot it\n",
    "f <- ols(time ~ ., test_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# multiple event regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'pd' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-b11c50bd80cf>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdf_corr_all\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"data/int_pc_qvals_gcamp_dms.csv\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex_col\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# first pull out the original gcamp. start with one mice (mice 9) for simplicity\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m data = pd.read_csv('./data/rachel_csv_6_19.csv', names = ['gcamp', 'TrialStart', 'NosePokeEnter', 'NosePokeExit', 'LeverPresentation',\n\u001b[1;32m      5\u001b[0m                                                        \u001b[0;34m'IpsLeverPress'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'ConLeverPress'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'CSplus'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'CSneg'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'RewardEnter'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'pd' is not defined"
     ]
    }
   ],
   "source": [
    "df_corr_all = pd.read_csv(\"data/int_pc_qvals_gcamp_dms.csv\", index_col = 0)\n",
    "\n",
    "# first pull out the original gcamp. start with one mice (mice 9) for simplicity\n",
    "data = pd.read_csv('./data/rachel_csv_6_19.csv', names = ['gcamp', 'TrialStart', 'NosePokeEnter', 'NosePokeExit', 'LeverPresentation',\n",
    "                                                       'IpsLeverPress', 'ConLeverPress', 'CSplus', 'CSneg', 'RewardEnter',\n",
    "                                                       'RewardExit', 'MouseID', 'RecordLoc', 'Session', 'HighProb'])\n",
    "#         1 = VTA/SNc -> Dorsomedial Striatum\n",
    "#         2 = VTA/SNc -> Nucleus Accumbens\n",
    "#         3 = VTA/ SNc -> Dorsolateral Striatum (forgot i had these)\n",
    "#         4 = VTA/ SNc -> Dorsomedial Striatum Cell Bodies\n",
    "#         5 = GFP controls (DMS and NAc together)\n",
    "data_g = pd.read_csv('./data/dan_table.csv')\n",
    "\n",
    "data['gcamp'] = data_g['gcamp']\n",
    "\n",
    "data['LeverPress'] = data['IpsLeverPress'] + data['ConLeverPress'] # add Lever Press in general \n",
    "\n",
    "# for gcamp (dt = dt of gcamp recordings)\n",
    "dt = 15\n",
    "LEN_BEFORE = dt\n",
    "LEN_AFTER = dt * 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalize\n",
    "\n",
    "def normalize(x, mu, std):\n",
    "    return (x - mu) / std\n",
    "\n",
    "\n",
    "# normalizes values (gcamp or q values across mice)\n",
    "def normalize_across_mice(df_src, values):\n",
    "\n",
    "    for mouse_id in np.unique(df_src['MouseID']):\n",
    "        mice_gcamp = df_src.loc[df_src['MouseID'] == mouse_id, values].values\n",
    "        mice_mu = (np.mean(mice_gcamp))\n",
    "        mice_std = (np.std(mice_gcamp))\n",
    "\n",
    "        df_src.loc[df_src['MouseID'] == mouse_id, values] = normalize(mice_gcamp, mice_mu, mice_std)\n",
    "    return df_src\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rpy2.robjects import IntVector, Formula\n",
    "from rpy2.robjects.packages import importr\n",
    "from rpy2.robjects import pandas2ri\n",
    "pandas2ri.activate()\n",
    "from rpy2 import robjects\n",
    "\n",
    "base = importr('base')\n",
    "rms = importr('rms')\n",
    "stats = importr('stats')\n",
    "\n",
    "import rpy2.robjects.numpy2ri\n",
    "rpy2.robjects.numpy2ri.activate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/rpy/lib/python3.6/site-packages/pandas/core/indexing.py:543: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  self.obj[item] = s\n"
     ]
    }
   ],
   "source": [
    "data_src = normalize_across_mice(data[data[\"RecordLoc\"] == 1], 'gcamp').reset_index()\n",
    "\n",
    "recordloc = \"DMS\"\n",
    "df_trial_src = normalize_across_mice(df_corr_all[df_corr_all['RecordLoc'] == recordloc].reset_index(), 'Q_ch_diff')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# building design matrix for different timelocks\n",
    "\n",
    "event_names =  ['NosePokeEnter' , 'LeverPresentation', 'LeverPress']\n",
    "design_matrix = np.zeros((len(data_src), len(event_names) * (LEN_BEFORE + LEN_AFTER)))\n",
    "for event_i, event in enumerate(event_names):\n",
    "    # shift the event indicator column accordingly for each event. \n",
    "    for shift_i, shift in enumerate(np.arange(-LEN_BEFORE,LEN_AFTER)):\n",
    "        design_matrix[:, (event_i * (LEN_BEFORE + LEN_AFTER)) + shift_i] = data_src[event].shift(shift).fillna(0).values\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# con/ips/reward previous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting cond: 0\n",
      "starting cond: 1\n",
      "starting cond: 2\n",
      "starting cond: 3\n",
      "CPU times: user 17.1 s, sys: 28.6 s, total: 45.7 s\n",
      "Wall time: 42.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# buliding design matrix for different conditions\n",
    "\n",
    "\n",
    "full_design_matrix = np.zeros((np.size(design_matrix, 0), np.size(design_matrix, 1) * 4))\n",
    "\n",
    "# adding dummy 'trial_starts' for fencepost problem for final trial \n",
    "trial_starts = np.append(np.where(data_src['TrialStart'] == 1)[0], len(data_src))\n",
    "\n",
    "for cond_i, (action, reward_last) in enumerate([('Con', 0.0), ('Con', 1.0), ('Ips', 0.0), ('Ips', 1.0)]):\n",
    "    print(\"starting cond:\", cond_i)\n",
    "    cond = df_trial_src.index[(df_trial_src['Action'] == action) & (df_trial_src['Reward_last'] == reward_last)].values\n",
    "    trial_starts_cond = trial_starts[cond]\n",
    "    trial_ends_cond = trial_starts[cond + 1]\n",
    "    \n",
    "    cond_design_mask = np.zeros(len(data_src), dtype = bool)\n",
    "    for (trial_starts_i, trial_ends_i) in zip(trial_starts_cond, trial_ends_cond):\n",
    "        cond_design_mask[trial_starts_i:trial_ends_i] = True\n",
    "    cond_design_matrix =  np.copy(design_matrix)\n",
    "    cond_design_matrix[~cond_design_mask,:] = 0\n",
    "    full_design_matrix[:, cond_i * np.size(design_matrix, 1) : (np.size(design_matrix, 1) * cond_i + np.size(design_matrix, 1))] = cond_design_matrix\n",
    "    \n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.55 s, sys: 14 s, total: 20.6 s\n",
      "Wall time: 25.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# add intercept term\n",
    "A = np.hstack((np.ones(len(data_src))[:, None], full_design_matrix))\n",
    "df_reg = pd.DataFrame(full_design_matrix) # fit automatically adds an intercept here "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_reg['gcamp'] = data_src['gcamp']\n",
    "# df_reg['mouseInfo'] = data_src['MouseID'] # can use cluster on a vector "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "print('performing least squares')\n",
    "\n",
    "# perform least squares\n",
    "w_hat = np.linalg.inv(A.T @ A) @ A.T @ data_src['gcamp'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-15 17:49:04.124245\n",
      "2019-01-15 18:03:07.280759\n"
     ]
    },
    {
     "data": {
      "application/javascript": [
       "$(document).ready(\n",
       "    function() {\n",
       "        function appendUniqueDiv(){\n",
       "            // append a div with our uuid so we can check that it's already\n",
       "            // been sent and avoid duplicates on page reload\n",
       "            var notifiedDiv = document.createElement(\"div\")\n",
       "            notifiedDiv.id = \"8485cd32-20e5-4b0a-826b-55eb1cd1923a\"\n",
       "            element.append(notifiedDiv)\n",
       "        }\n",
       "\n",
       "        // only send notifications if the pageload is complete; this will\n",
       "        // help stop extra notifications when a saved notebook is loaded,\n",
       "        // which during testing gives us state \"interactive\", not \"complete\"\n",
       "        if (document.readyState === 'complete') {\n",
       "            // check for the div that signifies that the notification\n",
       "            // was already sent\n",
       "            if (document.getElementById(\"8485cd32-20e5-4b0a-826b-55eb1cd1923a\") === null) {\n",
       "                var notificationPayload = {\"requireInteraction\": false, \"icon\": \"/static/base/images/favicon.ico\", \"body\": \"Cell execution has finished!\"};\n",
       "                if (Notification.permission !== 'denied') {\n",
       "                    if (Notification.permission !== 'granted') { \n",
       "                        Notification.requestPermission(function (permission) {\n",
       "                            if(!('permission' in Notification)) {\n",
       "                                Notification.permission = permission\n",
       "                            }\n",
       "                        })\n",
       "                    }\n",
       "                    if (Notification.permission === 'granted') {\n",
       "                    var notification = new Notification(\"Jupyter Notebook\", notificationPayload)\n",
       "                    appendUniqueDiv()\n",
       "                    notification.onclick = function () {\n",
       "                        window.focus();\n",
       "                        this.close();\n",
       "                        };\n",
       "                    } \n",
       "                }     \n",
       "            }\n",
       "        }\n",
       "    }\n",
       ")\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%notify\n",
    "print(datetime.datetime.now())\n",
    "rdf_reg = pandas2ri.py2ri(df_reg)\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "# took about 15 minutes for DMS "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-15 18:03:10.709603\n",
      "2019-01-15 18:16:07.431538\n"
     ]
    },
    {
     "data": {
      "application/javascript": [
       "$(document).ready(\n",
       "    function() {\n",
       "        function appendUniqueDiv(){\n",
       "            // append a div with our uuid so we can check that it's already\n",
       "            // been sent and avoid duplicates on page reload\n",
       "            var notifiedDiv = document.createElement(\"div\")\n",
       "            notifiedDiv.id = \"bf65ad27-ac35-4d3d-b97b-cbe263d609db\"\n",
       "            element.append(notifiedDiv)\n",
       "        }\n",
       "\n",
       "        // only send notifications if the pageload is complete; this will\n",
       "        // help stop extra notifications when a saved notebook is loaded,\n",
       "        // which during testing gives us state \"interactive\", not \"complete\"\n",
       "        if (document.readyState === 'complete') {\n",
       "            // check for the div that signifies that the notification\n",
       "            // was already sent\n",
       "            if (document.getElementById(\"bf65ad27-ac35-4d3d-b97b-cbe263d609db\") === null) {\n",
       "                var notificationPayload = {\"requireInteraction\": false, \"icon\": \"/static/base/images/favicon.ico\", \"body\": \"Cell execution has finished!\"};\n",
       "                if (Notification.permission !== 'denied') {\n",
       "                    if (Notification.permission !== 'granted') { \n",
       "                        Notification.requestPermission(function (permission) {\n",
       "                            if(!('permission' in Notification)) {\n",
       "                                Notification.permission = permission\n",
       "                            }\n",
       "                        })\n",
       "                    }\n",
       "                    if (Notification.permission === 'granted') {\n",
       "                    var notification = new Notification(\"Jupyter Notebook\", notificationPayload)\n",
       "                    appendUniqueDiv()\n",
       "                    notification.onclick = function () {\n",
       "                        window.focus();\n",
       "                        this.close();\n",
       "                        };\n",
       "                    } \n",
       "                }     \n",
       "            }\n",
       "        }\n",
       "    }\n",
       ")\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%notify\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "\n",
    "fmla = Formula('gcamp ~ .')\n",
    "\n",
    "\n",
    "fit = rms.ols(fmla, rdf_reg, x = True)\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "# took 7 minutes for DMS_CB \n",
    "# took around 13 minutes for DMS "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c20a46b00>"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEJCAYAAAB/pOvWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VMX6wPHv7qb3XkhCQgIcIKFIkSJNqggolyglAoYLXrAEBL2Iwo0FFS4dBNQLqD8iSlSKiDRFUUCKCBikHJBAIIE00uum7O+PzS6E3SQbSIPM53nyLJw5c86sD+Y9c2bmHYVGo0EQBEEQbqes7wYIgiAIDY8IDoIgCIIBERwEQRAEAyI4CIIgCAZEcBAEQRAMiOAgCIIgGBDBQRAEQTAggoMgCIJgQAQHQRAEwYAIDoIgCIIBERwEQRAEA2b13QBTSZJkCXQBbgAl9dwcQRCE+4UK8AZ+l2W50NRK901wQBsYDtR3IwRBEO5TvYCDpp58PwWHGwAbN27Ey8urvtsiCIJwX0hMTOSZZ56Bst+hprqfgkMJgJeXF76+vvXdFkEQhPtNtV7HiwFpQRAEwYAIDoIgCIIBERwEQRAEAyI4CIIgCAZEcBAEQRAMiOAgCIIgGLifprIKQv0pLYG1/aDDM9D1X3V228zMTFJTU1Gr1XV2T6HhU6lU2Nvb4+LigqWlZa3cQwQHQTDFpZ/gximwcqyz4FBQUEBSUhK+vr5YW1ujUCjq5L5Cw6bRaCgqKiIrK4urV6/StGnTWgkQ4rWSIJjiZJT2M+GEthdRB1JSUnB3d8fGxkYEBkFPoVBgYWGBm5sbzs7OpKWl1cp9RHAQhKrkpsL5neDoB+psSL1QJ7ctKCjAzs6uTu4l3J8cHBzIzs6ulWuL4CAIVYn5CkqLYPD72r/H/14nty0uLsbMTLz5FSpmbm5OSUnt9GRFcBCEymg0cPJz8OkMrYaBlRPEH6+z24vXSUJlavPfR7UfSyRJ6gr8BvSRZbnK9K+SJH0DhFZQ3MuUawhCvbl+EpLPwLDloFSCT6c6DQ6CUF+qFRwkSQoCtlC9Hkc74E9gsZEyuTr3F4Q6dzIKzKwhZKT2775d4NeFUJgDlmI8QHhwmRwcJEkaBnwGuFajjg0QBKyWZfnzardOEOqTOg9OfwNtntROYQXw7QyaUm2Polmv+m2fINQik3oAkiR9AXwHJAFfVuP6bcvu8Vf1myYI9ez8DijMgofG3Trm00n7WUeD0o3djz/+yAsvvEDfvn0JCQmhe/fuPP/88xw+fLi+mwZAWloakZGR9OnTh/bt2xMaGsrevXtNrn/t2jVmzpzJI488wkMPPcS4ceM4duxYLbbYdKa+HmoNzAE6AtWZx9e+7PMMgCRJtpIkiUFw4f5wMgqcm0FAz1vHbFzAJQgS/qi/djUC2dnZPP/887z44oskJyczatQoIiMjGTVqFKdPnyY8PJx169bVaxtzc3MJDw9n69atDB06lNmzZ6NSqYiIiGDbtm1V1k9KSmLcuHEcPHiQMWPGMHPmTG7evMnEiRMbRPAz9bVSV1mW1QCSJFXn+rrgMEaSpK/RbnKdJ0nSFmCmLMsp1bmYINSZtMtw+VfoNxfunBHi2xli92tnMonZRLXijTfe4KeffmLOnDlMmDChXNlzzz3H+PHjWbRoEc2bN6dv37710saNGzciyzKrVq1i4MCBAISGhjJ69Gjmz5/PoEGDsLGxqbD+6tWrSU5O5uuvvyYkJASAJ554guHDhzNv3jx27txZJ9+jIiY9xesCw11oV/bZBXgdeArYAIQBByVJcrzL6wpC7Tr1BaCA9mGGZb5dICcJMuPrvFmNwS+//MLevXsZNmyYQWAAsLOz45133gEgKiqqrpunt23bNnx8fPSBAcDCwoIJEyaQkZHBr7/+WmHdkpISvvvuOzp37qwPDACOjo6MGjWKS5cucfr06Vptf1Vqe4XNZ8AeYIEsy8VlxzZLkiQDy4BXgMhaboMgVE9piTY4NO8Pjj6G5bePOzj51W3bGoHt27cDEBZmJDCXadu2LTt27CAoKKjc8fPnz7NmzRqOHj1Kbm4uPj4+PP7440yZMgUrKyv9eZIkMWHCBDp06MDatWu5dOkS9vb2DBgwgFdeeQVHx8qfW7Ozs4mNjWXIkCFG2wYQExPDY489ZrT+xYsXycvLo3379gZlt9fX/bk+1GpwkGV5fQVFa4BFwCBEcBAamtj9kBUPg98zXu4ZAmZW2nEH3RTXOrT5j3i+On6tzu9rqlGd/Qjt5HvX9WNiYjAzM6vyF2OLFi3K/f3IkSM899xzWFlZMXbsWDw9PTl69Chr1qzh0KFDbNiwoVyA2LdvH1u2bGHs2LGEhYWxb98+oqOjyc7OZtmyZZXeOykpCY1Gg5eXl0GZp6cnAAkJCZXWB/D29r6r+nWhXtbmy7KsliQpA7Cvj/sLQqVORoG1C0iGT4UAmFmAd3uxGK6WpKSk4OzsjIWFhcl1SktLmTt3Lkqlks2bN9O0aVMAnnnmGVavXs3KlStZv349L774or7O9evX2bJlC23atAHg6aefZujQoezdu5f8/Hysra0rvJ8un5GxMQVdAMrPz6+yvrF7mFK/LtRacJAkyQ/YCfwhy3L4HWXugBvaldaC0HDkpcH576HzJDCrJA2ybxf4fR2UFIHKvO7aB4R28r2nJ/OGTqVSUVxcXPWJtzlz5gzXrl1jzJgx+sCg89xzz7Fu3Tp2795dLjgEBQXpAwNoU1EEBwdz6dIlMjIyKg0OGo1GX+dOumOVpba41/p1oTanlSYADsAoSZJa3lGm669/Vov3F4TqO/01lKjhoWcqP8+nExQXQJJYwlPTPDw8yMzMrNYGR9euaV+z3fmqCbSDxP7+/sTHl59A4OpquJ5X11vRJbNLS0sjJSWl3E9BQYG+x2Ds6b6goACg0oy6ldXXHbO3r98XKzXWc5AkqR3a2Um/ybIcK8tyqSRJLwDfAgckSVoNpAFPAAOBDbIsb62p+wtCjYiJ1r4y8qpiINC3i/Yz/jg0eaj229WIdOnShdjYWE6cOEG3bt0qPG/atGnY2tryn//8R/8kXpGSkhKD11RKZdXPxk899ZTBu3/dNFW4NXZwu8TERACj4xE6vr6+FdbXHdONPdSXmnytNBJ4E5gIxALIsvy9JEl9gP+gnZlkiTaf0kvAhzV4b0GoGelXoM2Iqs9z9AU7T21wePi5Wm9WYzJkyBCio6P54osvKgwOFy5cYM+ePQQFBWFjY4Ofn3bW2MWLFw3OLSwsJD4+Hn9//2q3ZdGiRRQWFpY71rx5c+zs7GjWrJnR6aYxMTEARmci6QQGBmJra8tffxn2PE2pXxeqHRxkWX4LeKsaxw8BxudzCUJDUloK+elg7Vz1uQqFNo23SKNR47p3786jjz7Knj17+OyzzwgPDy9XnpaWxsyZMwF4+eWXAQgODsbX15dt27YxadKkcuMOa9euJS8vr9x6BFN16tSpwrJhw4bxwQcf8NNPP9GvXz8A1Go1UVFRuLi40Lt37wrrmpubM3jwYL799lvOnj2rH/vIysrim2++oWXLluXGQ+qD2ElEEHQKs7RJ9WxcTDvftzPI32sHsU2tI5hkwYIFTJ06lfnz57Njxw4GDRqEo6MjsbGxbNmyhezsbCIiIvSvd1QqFfPmzWPKlCmEhoYSFhaGh4cHx44dY/fu3QQHBzNp0qQabePEiRPZvn07M2bMIDw8HG9vb7Zu3crZs2dZtGhRuWmzJ0+e5OrVqzzyyCO4ubkBMH36dH7++WcmTpzIxIkTsbOz44svvuDmzZssXLiw3gekRXAQBJ38dO2nKT0H0AYH0O4r3WJA7bSpkXJycmLDhg1s376db7/9lqioKNLT03FwcKBbt248++yzdO7cuVydHj16EB0dzZo1a9i0aRN5eXn4+voybdo0Jk+ejKVlJbPP7oKtrS0bN25kyZIlbNq0CbVaTVBQEKtXr2bAgPL/HqKjo9m6dSsbNmzQBwcvLy82bdrE4sWLWbt2LQBt2rThnXfeMfhu9UFR1UBOQyFJUgBwed++ffrBHEGoUQl/wNp+MHZTxWscbleYDQuaQu9Z8OjrNd6cc+fO0bp16xq/rvBgqerfSXx8PP379wdoJsvyFVOvKzKkCoKOvudg4isiS3twbw0JYjGc8OARwUEQdPKq+VoJtK+W4o9rM7QKwgNEBAdB0NH1HKozuOzbGQoy4Oal2mmTINQTERwEQUcXHKycTK/jUzZwKKa0Cg8YERwEQSc/DSwdQVWNSXzuEljYi3EH4YEjgoMg6OSng3U1eg0AShX4PCQytAoPHBEcBEEnL616g9E6vl20CfiK6jfFsiDUJBEcBEEnP/3uVjr7dIbSYrjxZ823SRDqiQgOgqCTf7c9B92gtHi1JDw4RHAQBJ38dNMXwN3OzgNsPSDlXM23SRDqiQgOggBlGVkz7q7nANpZSylyzbZJEOqRCA6CANqFbGjuITi00gYHsVJaeECI4CAIcHero2/nLmlTfmffqLk2Cfz444+88MIL9O3bl5CQELp3787zzz/P4cOH67tpBv78809at27N8eOmjz2lpaURGRlJnz59aN++PaGhoezdu7cWW2k6ERwEAaqfrvtO7q20nynna6Y9jVx2djbPP/88L774IsnJyYwaNYrIyEhGjRrF6dOnCQ8PZ926dfXdTL2rV6/y0ksvUVpaanKd3NxcwsPD2bp1K0OHDmX27NmoVCoiIiLYtm1bLbbWNGI/B0GA6mdkvZM+OMgQ1K9m2tSIvfHGG/z000/MmTOHCRMmlCt77rnnGD9+PIsWLaJ58+b07du3fhpZ5ueff2b27NlkZGRUq97GjRuRZZlVq1bpd6kLDQ1l9OjR+n2qbWxsaqPJJhE9B0EA7QI4uPueg62bNrCInsM9++WXX9i7dy/Dhg0zCAwAdnZ2vPPOOwBERUXVdfPKeeWVV5g6dSpubm4MGzasWnW3bduGj49Pue1LLSwsmDBhAhkZGfz666813dxqEcFBEODexxwUiluD0sI92b59OwBhYWEVntO2bVt27Nih30FN5/z580ybNo2uXbsSEhLC4MGDWbFiBQUFBeXOkySJ9957j++//54RI0bQtm1bevToQWRkJJmZmSa39dKlS8yYMYOtW7cSEBBgcr3s7GxiY2Np37690e8GEBMTY/L1aoN4rSQIcFtGVse7v4a7BGe2amcs1fP+v/ezmJgYzMzM9L8kK9KiRYtyfz9y5AjPPfccVlZWjB07Fk9PT44ePcqaNWs4dOgQGzZsKLev8759+9iyZQtjx44lLCyMffv2ER0dTXZ2NsuWLTOprV999RUWFhbV/o5JSUloNBq8vLwMyjw9PQFISEio9nVrkggOggDa1dFWjtpEehVIzi7AwcocK/MKznGXtFNic1O0C+Nqy6kv4eTntXf9e/XQOOgw9q6rp6Sk4OzsXK1fuqWlpcydOxelUsnmzZtp2rQpAM888wyrV69m5cqVrF+/nhdffFFf5/r162zZsoU2bdoA8PTTTzN06FD27t1Lfn4+1tbWVd73bgIDaHsOgNExBV0Ay8+v31xd4rWSIECVq6OTswvov/gXlv1woeJruEvaTzHucE9UKhXFxcXVqnPmzBmuXbvGiBEj9IFB57nnnsPGxobdu3eXOx4UFKQPDAAKhYLg4GCKi4urPbhcXZqy9TAKIz1M3TFjZXVJ9BwEAarMyLrshwtkFxZzOPZmxde4fcZSs9413MDbdBh7T0/mDZ2HhwdXrlxBrVab/GR+7do1wPBVE2if7v39/YmLiyt33NXV1ei5ACUlJYB2HYLuzzr29vblXk/dDV2PwVjvQDc+Ymdnd0/3uFei5yAIUGlG1vOJWUT/fg17SzPOXs8iX11i9DzsvcHSQQxK36MuXbpQWlrKiRMnKj1v2rRpvP766+Tl5emfxCtSUlJiEGiUyqp//T311FP07Nmz3M/OnTur/hJV8PX1BbRjD3dKTEwEMDoeUZdEz0EQQBscXIOMFr33/TnsrcyJHNaGV77+k5j4DLoGGj51amcsSeK10j0aMmQI0dHRfPHFF3Tr1s3oORcuXGDPnj0EBQVhY2ODn58fABcvXjQ4t7CwkPj4ePz9/avdlkWLFlFYWFjuWPPmzat9nTvZ2dnRrFkzTp8+bVCmm6VkbCZTXRI9B0GAsnTdhj2H/XIyBy6mEtGvOY+20g4y/3E1veLriAR896x79+48+uij7Nmzh88++8ygPC0tjZkzZwLw8ssvAxAcHIyvry/btm3j6tWr5c5fu3YteXl55dYTmKpTp0706NGj3I+HR81MNhg2bBhXrlzhp59+0h9Tq9VERUXh4uJC7961+GrSBKLnIAilJVCQaTDmUFxSyvs7z+HvasOE7gFYmCkJdLPlRFwlg5XurbQzifLS7n7NhMCCBQuYOnUq8+fPZ8eOHQwaNAhHR0diY2PZsmUL2dnZREREMGjQIEA7iD1v3jymTJlCaGgoYWFheHh4cOzYMXbv3k1wcDCTJk2qt+9z/vx5ZFmmY8eO+l7OxIkT2b59OzNmzCA8PBxvb2+2bt3K2bNnWbRo0T2Pa9wrERwEIb/sl/0dwSH6+DUuJOXw0biOWJhpO9kd/Z356XwyGo3G+GyS2wel/bvXZqsfaE5OTmzYsIHt27fz7bffEhUVRXp6Og4ODnTr1o1nn32Wzp07l6vTo0cPoqOjWbNmDZs2bSIvLw9fX1+mTZvG5MmTsbS0rKdvAz/88AOrVq1i/vz5+uBga2vLxo0bWbJkCZs2bUKtVhMUFMTq1asZMGBAvbVVR1HVQE5DIUlSAHB53759+sEcQagRqX/Dqk4wci20GwVAdkERjy7eT6CbHdFTuukDwRdHr/LG1tPsf7UvAW62htfKuArL28Kw5dB54j0169y5c7Ru3fqeriE8+Kr6dxIfH0///v0BmsmyfMXU64oxB0HIN8yr9OH+S6TmqJkztHW5HkInf+05f8RVMO7g4AvmtmLcQbjvieAgCHdkZE3IyGf9wcuM6NCE9n5O5U5t4WGHvaVZxYPSSiW4txQzloT7nggOgqAPDtpAsGi39hf7vx9rZXCqUqmgQ1MnTlTUcwBwEzOWhPufCA6CoEvXbePCqWsZbDt1nUk9m+HjZDy3Tid/Z+SkbLILioxfz12C7OvaGVCCcJ8SwUEQ8tMBBVg6smSvjJudBc/3Nb4gDqBjU2c0GvjzWgW//PUzlirJwyQIDZwIDoKQn6Z9paRUcuZ6FoOCvbC3Mq/w9A5NnVAoKhmUFgn4hAeACA6CUJaRtbiklPQ8Ne52lc+Hd7Ayp6WHPScqGpR2DgCVpQgOwn1NBAdByE8Ha2fSctVoNOBmX/ViqY7+zpy4mk5pqZF1QkoVuLUUg9LCfU0EB0EoS3WRkqNNsOZmW3Wa6I5NncguKOZSSo7xE0SOJeE+V+30GZIkdQV+A/rIsnzQhPPdgXnAUMAVOAPMl2V5S3XvLQi1Ij8d3CVu5qgB03oOty+Ga+Fpb3iCeyv46xtQ54KFkZXUgtDAVavnIElSELDF1HqSJNkBPwLhwJfATKAE2CxJ0oRqtVQQakvZmEOqrudQxZgDQDM3W5xtzKselE4VM5aE+5PJwUGSpGHAUaBJNa7/EtAOGCvL8ixZlj8C+gAngaWSJIlHKqF+lRRBYRZYO+uDg6td1a+VFAoFHZs6VzwofXsCPkG4D5naA/gC+A5IQtsDMNUE4Iosy1t1B2RZLgRWoH3FNKQa1xKEmndbRtbUHDUWZkrsLU1729rR35lLKblk5KkNC12agdJczFgS7lum9hxaA3OAjoBJ/WRJkhyBVmh7G3f6vezzYRPvLwi1Q5c6w0b7WsndztLkjd07NtWOO5y8amR/B5W5dmc50XO4K0ePHkWSJGbPnl0v98/NzWXx4sX079+fdu3aMWzYMKKjo02un5aWRmRkJH369KF9+/aEhoayd+/eWmxxzTM1OHSVZfn9sqd+U/kACuCakbKEss+AalxPEGrebXmVUnPUuJnwSkmnvZ8jKqWi8nEH0XO475SUlPDSSy+xfv16evXqxZw5c/D09CQyMpJVq1ZVWT83N5fw8HC2bt3K0KFDmT17NiqVioiICLZt21YH36BmmNR/lmXZSL+5So5ln7lGyvLKPsWYg1C/9Om6XUjNzsXL0fTdt2wszGjtXcliOPdWcO47KMoHc+N5moSGZ/fu3fz222+8/vrrhIeHAzBq1CimTp3KRx99RGhoKN7e3hXW37hxI7Iss2rVKv3WpKGhoYwePZr58+czaNAgbGxs6uKr3JPaXOeg65sb201IU0mZINQdfc/BmZu5hdXqOQB0aurMqWsZFJeUGha6S6AphZt/10BDhbqybds2LC0tGTt2rP6YQqFg0qRJFBUVsXPnzirr+/j4lNuz2sLCggkTJpCRkcGvv/5aa22vSbUZHHSrg4yFSN0xkbZSqF9lGVlLrZy5maM2aRrr7Tr6O5OnLkFOyjYsFDOWasyWLVuQJIkjR44wd+5cHn74YTp16sSkSZOIiYkpd25ubi5vv/02AwYMICQkhO7duxMREcG5c+dMuldMTAytWrUy2Fa0bdu2+vKKZGdnExsbS/v27Q3KTKnfkNTmHtJXyj59jJTp9vmMr8X7C0LV8tNBoSKz1JriUk31g0PZoPSJuHSCmziWL3RtDgqlGHeoQW+88QZmZmZMmTKFvLw8NmzYwPjx44mKiqJdu3YATJ8+nWPHjjFu3DiaNWtGUlISUVFRHDx4kN27d+Pp6Vnh9fPz88nIyKBr164GZdbW1jg4OJCQkGCkplZSUhIajQYvLy+DMt19K6vfkNRacJBlOUuSJBnoYqRYN0vJ2EwmQag7ZRlZU3O1w2qmrHG4na+zNR72lvwRl8747gHlC80swSWwxoPD9kvb2Xpxa9Un1pN/tPgHTwQ9USvXLikpYdu2bTg4OAAwcOBARo4cycKFC/n8889JS0vjwIEDhIWFMWvWLH29Nm3asHjxYs6cOVNpcMjO1vYAKxoTsLKyIj8//67qW1lpx7Mqq9+Q1HZupS+AlpIkDdcdkCTJApgGpAC7avn+glA5/epobXCoKiPrnW4thjMynRW0r5bEa6UaEx4erg8MAK1ataJ3794cP36ctLQ07OzssLe3Z9euXURHR5OSkgJAv3792LlzJ/369av0+hqNdhi0ounMCoWi0qnOldXXHTN1qnR9q7GegyRJ7dCuhv5NluXYssNLgfFAtCRJy4CraFNpPASMk2X5/gihwoOrLCOrPnWGCXmV7tTJ35ndZxJJyS7E/c767hLIu6BYDWbV65VU5ImgJ2rtybyha9GihcGxwMBAfv75Z+Lj42nXrh3z58/n9ddfJzIyksjISFq2bEmvXr0IDQ0lKEi7iVNGRgZFReV38rOxscHWVjuBsqKn+4KCApo0qThJhK7HYKx+QUEBAHZ2diZ80/pXkz2HkUAU0Ft3QJblHKAX8BUwBVgCqIARsixvrMF7C8LdKcvIWp28Snfq6K/de9rolFb3VqApgbRL99RMQcvc3HATptJS7UwxMzPts+7AgQP55ZdfWLZsGaGhoeTm5rJ+/XqGDx/Orl3alxURERH07Nmz3M8nn3yCnZ0dTk5OJCUlGdwnPz+fzMxMo+MJOr6+2uFUY/UTExMBKq3fkFS75yDL8lvAW9U4noi2tyAIDU9+BngGk5pTiEqpwMm64h3gKhLcxBFzlYJT1zIYHHzH//i6BHzJ58CjdQ00uHGLi4szGCyOjY3FzMwMPz8/cnNzkWUZHx8fHn/8cR5//HEAjh8/Tnh4OGvXrmXIkCG89tprZGVllbuOn58fACEhIZw4cYKioqJywUg3y8jYTCQdOzs7mjVrxunTpw3KTKnfkIj9HITGLT8NrF24maPGxdYCpbL674OtzFU097Dn3I0sw0K3lqBQQfLZGmisEBUVhVp9a03umTNnOHDgAL169cLe3p5Lly4xduxYPvzww3L1goODsbCw0PcuQkJC6NGjR7kfXXAYPnw4eXl5bNq0SV9fo9Hw6aefYm5urg84FRk2bBhXrlzhp59+0h9Tq9VERUXh4uJC7969K6ndcNTmVFZBaNiK1aDO0Y45JBXe1SslnTbeDhy4mGJYYG6tDRA37o+57Q1dXFwcY8aM4cknnyQ9PZ2oqCgcHBz0OZjatWtH7969+fLLL8nOzqZLly4UFhayfft28vLymDhxYpX3eOKJJ4iOjmb+/PnExcUhSRJ79uzhwIEDzJgxo9xsp/PnzyPLMh07dtQHl4kTJ7J9+3ZmzJhBeHg43t7ebN26lbNnz7Jo0SL9rKWGTgQHofEq0GVkdSKlmnmV7tSmiQObT8QbH5T2agtXqtwXSzDB7NmziYmJYcWKFZibm9OnTx9mzJih/8UMsHz5ctatW8euXbvYt28fKpWK4OBgPv74Y/r06VPlPZRKJWvXrmX58uXs3r2br776Cn9/f959912efvrpcuf+8MMPrFq1ivnz5+vbYGtry8aNG1myZAmbNm1CrVYTFBTE6tWrGTBgQM3+B6lFIjgIjVfZ6mhsXLiZU0ig292n+mrjrZ1eee5GFu727uULvdrC6a8gNxVs3e76Ho1J165dkWXDKcC2trYsWLCABQsWVFjX1taW6dOnM3369Lu+v52dHXPnzmXu3LmVnhcREUFERITBcXd390rbeD8QYw5C41WWV0ljpZ3KWlXPQX31KiVZRsYVuBUczhobd/DWrtwlUbxaEu4fIjgIjVdZRtZ8c0cKikorHXMoun6d2Cee5NJjQ8jatUu/2EnH0cYcHydrzl43Ehy8dMHBcAaLIDRUIjgIjVdZzyGtRPs6ybWS4JC8bDmUlmLu5UXCjJnER0RQlJxc7pw2TRyM9xxsXMDBVwxKC/cVERyExqtszCGlVLuqtaLXSvmnT5P13Xe4hIcT8FU0Hq++Qu6Bg8QOG07G5i36XkQbbwdiU3LIV5cYXsSrreg53IORI0ciyzJPPvlkfTel0RDBQWi88tNBaUZSvnahk7HXShqNhqQF/0Xl6orrv55DYWaG6+TJNNu2FcuWLbgxZw7XJj+HOj6BNk0cKNVgPH23dztooG+8AAAgAElEQVS4eRHUeYZlgtAAieAgNF66vEplGVkNpqAC2T/8QP4ff+AeEYHqtpw4ls2a4b9hA56R/yH/5Elin3iCln+fBKhg3KGtduMfsRhOuE+I4CA0XmWro3V5lVxsy79W0qjVJC9egmWL5jg9FWpQXaFU4hIWRuCO77AMDKTwzTdol3eDszeM7GGlG5S+8WeNfw1BqA0iOAiNV1nP4WaOGicbc8xV5f93SPviC4quXsVj1iwUZhUvCTJv0gS/jz7EzMWFOYfWE38hzvAkp6Zg5SjGHYT7hggOQuOVl67PyHrneENxejqpaz7EtmdP7Hr1qvJSZm5u+H70IdbFhTy9ZTlFObnlT1AotL0HsdZBuE+I4CA0Xrft5eB6xyul1DUfUpqTg8esf5t8OauWLUmYNgf/jOvEzngVTVkqaT2vtpB0BkqKa6L1glCrRHAQGq/8NP1rpds3+Sm8fJn0L7/E6amnsGrZslqXbDpkAP9r+ySlB/aTsnRp+UKvdlBcIPZ2EO4LIjgIjVNRARTlgbUzKTmF5bYHTV68BKWFBe7TDHPmVKWFhz27WvTkcs/HuLluPRnffHOr0Kut9lMshhPuAyI4CI1TWUbWIgsnsguK9Qvgco8eI2ffPlz/9S/M3G4lycsoyGD8zvGsjVlb6WUtzJQ093Tg6+6jsH3kEW689Ta5R49pC90lUFmIcQcTHD16FEmS9Km461N8fDwPPfQQ3377rcl1cnNzWbx4Mf3796ddu3YMGzaM6OjoWmxlzRPBQWicylZHZyntgVupM1JWrsTM2xuX8Gf1pxaVFPHy/pc5lXKKlSdXsu/qvkov3cbbgTOJufgsX4ZFgD/x06ahjosDlbl2NzgRHO4baWlpTJ06lbw80xcvlpSU8NJLL7F+/Xp69erFnDlz8PT0JDIyklWrVtVia2uWCA5C41SWVylTo82r5GZniaa0lIK//sJh8GCUZRuyaDQa3jr8Fn8k/cG8R+YR7BrMnINzuJx5ucJLt2niQHJ2IWlY4PfRRwAkzJqFprj4VhqNOxL3CQ3PqVOneOqpp7h48WK16u3evZvffvuN1157jbfeeovRo0ezbt06+vbty0cffcSNGzdqqcU1SwQHoXEqy8iaWqoLDhYUJyWhKSzEIsBff9ra02vZfmk7L3Z4kRHNR7Cs7zIslBbM+HkGeUXGnyZv39vBwtcX7zcjKfgzhpvr1oFXe8i7CVnXa/kLCvdiyZIljBkzhpKSEkaPHl2tutu2bcPS0pKxY8fqjykUCiZNmkRRURE7d+6s6ebWChEchMaprOeQXHSr56COuwqAhb82OOy+vJsPTn7A8MDhTGk3BQBvO28W9lnI5azL/OfQfwxSd4Ph3g4Ojz+Ow9ChpKxaTX62tkwshqueLVu2IEkSR44cYe7cuTz88MN06tSJSZMmERNT/jVdbm4ub7/9NgMGDCAkJITu3bsTERHBuXPnTL7f+fPnmThxIt9//z0dOnSoVltjYmJo1aoVlpbl1860bdtWX34/EMFBaJzKgkOC2hrQBQftymYLf39OJZ9izsE5dPToyFs93kKhUOirdvPuxrSHprE3bi8bzm4wuLSxvR28Iv+DmYsL15duoLQEMe5wl9544w2OHTvGlClTCA8PJyYmhvHjx5f7hTt9+nQ2b97MoEGDePPNN3nmmWc4duwYYWFhJCUlmXSf1atX89prr2F3Wz4tU+Tn55ORkYGXl5dBmbW1NQ4ODiQkJFTrmvVFbBMqNE55aaCyIClPia2FCmsLFVlX41BYWHDDtojpu6fjZevFikdXYKEyTOX9z5B/cjr1NMv+WEYb1zZ08epSrvzOvR1Ujo54v/8+1yZPJsXeD8+Quw8OGdu2kbl5y13Xr22OoSNxGjGiVq5dUlLCtm3bcHDQ9sAGDhzIyJEjWbhwIZ9//jlpaWkcOHCAsLAwZs2apa/Xpk0bFi9ezJkzZ/D09KzyPhYWd7efeHa2NiOvjY2N0XIrKyvy8/Pv6tp1TfQchMbptoysugVw6rg4VH4+RPw8jeLSYlb3X42TlZPR6gqFgncfeRc/ez9e/eVVEnMTy5Ub29vBrucjOIeFkfZnCbl/iAR8dyM8PFwfGABatWpF7969OX78OGlpadjZ2WFvb8+uXbuIjo4mJSUFgH79+rFz50769etXq+3TvWa8vad5O4VCUWFZQyN6DkLjdFtGVl3qjKK4OP62zeFqdjb/G/g/AhwDKr2EnYUdyx9dztjvx/LKL6/w6eBP9b2M2/d26OB3K8B4/PtVcvd9z/UfbxKYkoDK3afaTXcaMaLWnswbuhYtWhgcCwwM5OeffyY+Pp527doxf/58Xn/9dSIjI4mMjKRly5b06tWL0NBQgoKCAMjIyKCoqKjcdWxsbLC1tb2n9unqV9Q7KCgooEmTJvd0j7oieg5C45SfcSt1Rtk0VvXVa5yxTme0NNrgNVFFgpyCmPfIPGJSYlj4+0L9cf2g9B17OyitrWkyazLF+SqS3o6sue/TSJibmxscKy3LYWVWljl34MCB/PLLLyxbtozQ0FByc3NZv349w4cPZ9euXQBERETQs2fPcj+ffPLJPbfPzs4OJycno2Mb+fn5ZGZmGh2PaIhEz0FonPLSwKUZqQmFdApw1k9jve6kpI9L62pdanDAYP5K/YvPznxGO/d2PBH0BL7O1thbmRnd28G69zDc2rxF6o+/YffDDzgMHFhT3+qBFxcXR9euXcsdi42NxczMDD8/P3Jzc5FlGR8fHx5//HEef/xxAI4fP054eDhr165lyJAhvPbaa2RllQ/cfn5+NdLGkJAQTpw4QVFRUblgphs0b9++fY3cp7aJnoPQOOWnU2rlRFqeutw01hsu0MLZ8NVFVaZ3nE4Xry68c/gdzqedR6FQ0MbbgXM3jGwZau+F28PWWPnYkxj5JsWpqff6bRqNqKgo1Gq1/u9nzpzhwIED9OrVC3t7ey5dusTYsWP58MMPy9ULDg7GwsJC37sICQmhR48e5X5qKjgMHz6cvLw8Nm3apD+m0Wj49NNPMTc31weshk70HITGKT+dfDNHNBrtAjjdNNZkFxWBjoHVvpyZ0oyFvRcy+rvRzPh5BpuGbaK1twNfHb9GaakGpbL8IKSiSVua9LtO7KYckpcspcn892vkaz3o4uLiGDNmDE8++STp6elERUXh4OCgz8HUrl07evfuzZdffkl2djZdunShsLCQ7du3k5eXx8SJE2u0PdeuXePEiRNIkkSrVq0AeOKJJ4iOjmb+/PnExcUhSRJ79uzhwIEDzJgxw6TZUg2B6DkIjU9RPhTnk1uWV8nNzhL11TiKzZTYNmmKlZnVXV3WzdqNJX2XkJiXyBsH36C1tx156hLi0oyspPZqi2XJ37hOGEfm1q3k/ylmL5li9uzZtGzZkhUrVvDll1/Sp08fvvnmGwICAvTnLF++nBdeeIEzZ86wYMECVq5cia2tLR9//DFDhgyp0fb8/vvvzJo1ix9++EF/TKlUsnbtWsLCwti9ezfz5s0jKSmJd999l6lTp9bo/WuT6DkIjY8urxLaBU66BXCpziqau1Zv/4Y7dfDowKwus3j/6Pt4BrYEmnP2ehbN3O6YBePdDkqLcB3Zm8ztO0h89z0CojehUIrnNYCuXbsiy7LBcVtbWxYsWMCCBQsqrGtra8v06dOZPn16jbRl5MiRjBw5slpldnZ2zJ07l7lz59ZIG+qD+JcoND5lGVnTNLfyKhVeucI1p+K7Gm+40xhpDMMCh/FN7Hos7C4YHZTGqx0Aqqy/8fj3qxScPk3m1m33fG9BqCkiOAiNT1nPIbVYGxxcbMxRX71KojO0dLq3ngNoFzpFdo+khXMLrHw2cfK6kZ3fXALB3BYST+MwfDjWHTqQvHQpJdlGBrAFoR6I4CA0PmUZWW8UW2OhUmKTlYZCXUSis6JGeg4A1mbWLO+7HJUS/ipeRUFxQfkTlCrwDIYbMSgUCjznzqUkLY3U1Wtq5P6CcK9EcBAan7Kew/UCa9zsLCgqm8aa5m6Jr71vjd3Gz8GPIZ4z0VjE85+D7xhmcPVup83OWlqKdUgwTk+Fkvb55xReEntM32nkyJHIssyTTz5Z301pNERwEBqfsuBwNd8K17KZSgDW/oEoFTX7v8SwFv0pTOnH7rjv2Hxxc/lCr7agzoaMKwC4v/wySmtrkt6fbzQVuCDUJREchMYnLw3MrEjMv7XGocgMPJq1qfFbBTdxRJ06AF/LDrx/9H3+Sv3rVqF32UrZhBMAmLm64h7xErmHDpHz00813hZBqA4RHITGR5eRNVu7Ojrn8t8kOUILl3sfjL6To7U5kqcjLnkTcbd2Z8b+GaQXaHsueLYFS0e4/Iv+fOexY7FoHkTS/AWi9yDUK5PXOUiS5A7MA4YCrsAZYL4sy1Umlpck6TjQqYJiP1mW401thyDcs/x0NNbO3EwrxM3ekvzLl7jhokCqocHoO3UOcGb7qet88dJSwndPYNavs/howEeoVGbQrBdc2q/dU1qhQGFujtcbb3D1n5MozclBo9HcNymehbpXmw8QJvUcJEmyA34EwoEvgZlACbBZkqQJVdRVAW2AfcB4Iz9pd9l2Qbg7+emUWDpRVKLB1doMRUISic53l1PJFA83cyG7sBhlkS9zu83lyI0jrD61WlsY2Bcyr0L6Zf35tj16YD9wIJqkJNQFBUavKQgARUVFqFSqWrm2qT2Hl4B2wEhZlrcCSJL0KXAYWCpJ0mZZlnMrqNsSsAZ2yrL8+b02WBDuWW4K+Y7aV0jexTko1cVke9jjYuVSK7frHKC97vEr6Tzb4x/8mfIna0+vJcQthH6Bj2pPit2vXftQxuO1WaRv3UpWfDzuRvYwEASArKws7O3ta+Xapo45TACu6AIDgCzLhcAKtK+YKktY0q7s869KzhGEupOTRI65KwBumcmAdt/o2uLjZE0TRyt+v6LtJL/e9XWCXYOZc3AOceZm4OCrDQ63sfD1xRkFSVeukJOWJsYfBD2NRoNarSY1NZX09HRcXGrnoabKnoMkSY5AK+ArI8W/l30+DHxTwSV0ycvPlF3PDsiVZVn8axfqXlEBFGSSqXQGwP7mDQCcmldvD4fq6hzgwtHLN9FoNFiqLFnadymjd4zm5f0z2NisFzbyLigt0S6OK+M9YTxZL8/gqlqNIrD6mWKFB5dKpcLe3p6mTZtiaWlZK/cw5bWSD6AArhkpSyj7DKikfntADbwmSdIzgAuQIUlSFPB6Ja+jBKHm5Wp7CjcV2q07S+JlNCrwbVa7G7B0CXBm+5/XiU/Px8/FhiZ2Tfhv7//y/I/P84o9rCzIwPzGn+DTUV9H5eiI16N9SXrlVfz+9zF2vXvXahsF4XamvFZyLPs09ktcl4u4so1X2wEWaMceXgTGALuBCGCPJEmG+/4JQm3J0QaHpBIHlApQX7tEkhO0dGtVq7fVjTvoXi0B9GjSg/90+w8Hs/7mLTcXNJd+NqjnPGYM5k2bkrxoEZri4lptoyDczpTgoJtHZ+w1kKaSMiRJUgDLgBmyLD8my/ImWZajZVkei3a84hG0M5YEoW7kaPf2TSh2wMXWktJrCdxwURDoVLuvbSRPe+ytzPj9Snq540+1fIoXOrzAdns7VlzealBPYWGBx8yZFF78m4ythuWCUFtMCQ45ZZ82Rsp0x4zkJAZZljWyLC+VZXm5keKVZZ+DTGiDINSMsp5DXKEd7rZmWCWmk+fpgLWZda3eVqlU0NnfmeNXDGduT203lVHW/qwnk41/fWZQbj94ENYdOpCyciWlueItrFA3TAkOV8o+fYyU6bKU3c0itqSyz9qZhyUIxpQFhysFNgSQj1lRKUo/Y/+0a17nABcuJueQnqsud1yhUPBGhwgG5Obx3z+WsvvyboNyj9dmUZKSys1PP6uTtgpClcFBluUsQAa6GCl+uOzzqLG6kiR1lSTprCRJ7xop1k0P+duUhgpCjchJAhtXknJLaFqgfT6xD6ybdQRdysYd/ohLNyhTBfRkQWomHS1ceP3g6xy5caRcuc1DD2H/2GPcXL+eouTkOmmv0LiZus7hC6ClJEnDdQckSbIApgEpwK4K6p0HmgKTy9Jv6OqaoU3FoQE23EW7BeHu5CSBnSep2WpcsrXPJV5Sxyoq1Yx2vo5YqJTlBqX1LGyx9OvKyqxiAhwCePnnlzl381y5UzxmzkBTXEzqBx/USXuFxs3U4LAU7RN+tCRJ70mSNAX4BXgI7WBzPoAkSd0lSRonSZIngCzLmcCrgCdwVJKkWZIkvQwcAh4D5smy/EfNfiVBqEROMiU27uQXlWCfFkeRCpq16Fwnt7YyV9HW19F4cAAI7IvDjb/46JH5OFg48PyPz3M957q+2KJpU1zCxpKxeQsFFy7USZuFxsuk4CDLcg7QC+1CuCnAEkAFjJBleeNtp04Borj1yghZlj8ChqMdl3gTeA9tjyFMluU3a+A7CILpcpIosHQDwC71OsnOCvwca2919J26BLhwOiGTgqISw8LAvoAGz6RzfDTwIwpLCnll/yuoS26NUbhOnYrSzo7kxYvrqslCI2VyVlZZlhPRJt6r7JxwY+fIsrwD2FG9pglCDdNoICeZbF9t6gz71HQyPexRKWsncZkxXQKc+egXDX9ey6BroGv5wiYPaVN4X/qZwOB/8G7Pd3n555dZcGwBkd0jATBzdsZt6lSSFy4k87sdOA4fVmdtFxoXsZ+D0HgUZkNxPplKJxSaUpxT8yn19azTJnTy16btOG5kUBpdCu+yPEv9m/bnnyH/5OsLX/Pt39/qT3OZMB7rjh1JfPNN1Feu1EGrhcZIBAeh8SibxnoTJ1yLErEs1m4NWpecbCxo6WnHscsVjzuQEQdp2hTeEQ9F8LDXw8w7Mg85TQZAYWaGz5LFKMzNiZ85k1K12vi1BOEeiOAgNB5leZWSSh3xLbwIgFvLkDpvRucAF07EpVNSaiSxQGBf7WdZ78FMacZ/e/8XRwtHZuyfQZY6CwBzb2+858+n8Ow5khcuqpN2C42LCA5C43Fb6oyAIm3OyKZtutV5M7oEOJNdWIycmG1Y6NocHHzKpfB2s3ZjSd8l3Mi5wZyDcyjVlAJg3+9RXJ59lvTPPyfrhx/qqPVCYyGCg9B46FJnqO3wyUuiSAUeAbWbqtsY3WK443FGXi0pFBD4qHZf6dJS/eEOHh14tcur7L+2n0/++kR/3OOVmViFhHBjzlyKEhIMrycId0kEB6HxyEkCpRlXcy3xysogy80aRS1tsVgZHydrvB2tDJLw6QX2hfx0SPyz3OGwVmEMCRjCByc/0K+gVlhY4LN0CZSWkjDzFTRFRbXbeKHREMFBaDxyksDWgxs5Bbin51HUxK1emqFQKOgc4MLvlyvY4S2wj/bzjt3hFAoFb/V4i2YOzZj1yyx+jPsRjUaDRdOmeM97h/w//yRlxYra/wJCoyCCg9B45CSjsfMgKTcBr3QN5v5N660pXQKcScwqID4937DQzgM8gg2CA4CNuQ3LH12Oq7UrM/bPYNLeSchpMg5DhuA0ejQ3160n59dfa/8LCA88ERyExiMniSIrNxyL/sayGJwC6368QaezfyXjDgBBj0LcYSgyDB4BjgF8Pfxr5nSdw8X0i4zaMYq3D7+NxYwpWLZsyfVZr1F48WJtNl9oBERwEBqPnGRyzV3xzb8CgHerukm4Z4zkZXzzH72gR6GkEC7sMVpspjRjTKsx7PjHDsJahbHt4jaG7RzJsYi+YG5O3IRnKTh3zmhdQTCFCA5C41BaCrkppCudaZKfCIB9YMt6a45KqaBTBZv/ANoZSy6B8NtKbdqPCjhaOvLaw6+x+cnNdPDowHvXP+HNcSqyFAVcHBfGmYPfUlQqBqmF6hPBQWgc8tOhtJhUHPHJuUmJmQJzb696bVKXABcuJBlu/gOAUgXdX4KEPyDutyqvFegYyIcDPmRN/zW4tQjh/Qm23DQvoOCF2Uz478OM3zmehb8vZO+VvWQUZNTCtxEeNCYn3hOE+1rZArhzRSq8s/JRe7rUyzTW2+nWOxz4O5Un2jcxPKFDGPz8PhxaAQGPmHTNXr696OXbC82jGhIGxpD2rwje2JROtF02X6V9RdTZKBQoaO3amu7e3enepDsdPDpgqbKsya8mPABEz0FoHMqCw4niDHxTNFgF1m1OJWM6+Tvj52JN1OErxk8wt4auU+DiHkiu3viBQqHAN6g9bTZtxsbPn3HrrvBj08VEDYnihQ4vYKWy4v/O/B+T906m55c9mfLDFNadXsdv138jvaCCcRChURE9B6Fx0O0dnXsdnzRweahrPTdIO+7wbPcA3v3+HH8lZBLi42h4UpfJcHAZ/PYBjFhT7XuYubvjv2EDVydNJvHFaQStWE6HflOZ2n4quUW5HE88zpEbRzh8/TArTtxaI+Fl60Vrl9a0dm1NG5c2tHdvj5OV0718XeE+I4KD0DiU9RzskuMAsG3bvj5bozeqix/LfrjAJ4cus3RUB8MTbFyg4wT4fT08Ogccfap9DzMXF/w//YSrz/2L+BdexH7wYNynT8c2sBl9/PrQx0+76C6zMJPzaec5d/McZ9POcu7mOfZf248GDeZKc4Y0G0JY6zCCXYPv9WsL9wERHITGISeJHAtb/K+mAGAV0jB+wTlYmfNUJ1++OHaV2UNa4WFvZXhStxfg2Fo4+iEMeveu7qNycsL/s0+5+cmnpH36Kdk//ojTyJG4vfQi5p7aPS0cLR3p6t2Vrt63elV5RXmcSzvHrsu72H5pO9svbae9e3ueaf0MA/wHYK40v6v2CA2fGHMQGoecZE7Zu9I8UUOOqzNmzs713SK98EeaUVSiYeORq8ZPcPaH4H/A8c+gIPOu76O0tcU94iWCftiLc1gYGdu2cWnQYJIXL6Ykw/gMJhtzGzp5dmJut7nse3ofs7rMIq0gjVm/zmLwN4P56M+PSMpNuus2CQ2XwmhulwZIkqQA4PK+ffvw9fWt7+YI95v/e4Ll6iQeWpeFdZtePPzp/+q7ReX887PfiYnP4NDsfliaGZlFdeNP+Lg3DHgber5cI/dUx8eT+sEHZG7/DqW9PW4vPI/Ls8+iUCgqrVeqKeVgwkG+OPcFh64fAqC1S2t6+/amt29vQtxCUCqq/9yZmp/KX6l/ce7mOUo0JVibWWNjbqP9NLPBxtwGO3M7gl2DMVeJHoup4uPj6d+/P0AzWZavmFpPvFYSGofcFM5pShmcCSVtjbzbr2cTHwlg/Ppj7PjzBqGdjDz8eLfXLow78iF0ex7M7n3qqYWvL03++19c/jmJ5CWLSV7wXwpiYvB+7z2U1tYV1lMqlPpAcCXzCj9e/ZED8QdYe3otH8d8jIuVCz19etLbtzcBDgGoFCqUSqX2U3Hr80rWFf5K/Uv/k5Sn7YEoUKCh4odWN2s3xkhjGCWNwtmq4fQAHzQiOAiNQlFOEuocWwCcO7ar59YY6tncjRYednxy6DIjO/oYf3p/ZDpEjYCYr6Dj+Bq7t5XUEr+PPyZt/XqSlyxFffUavqtX6cciKhPgGMDktpOZ3HYyGQUZHLp+iF/jf2X/tf1sv7TdpPs3tW9KR8+OhLiGEOIWQiuXVlibWVNQUkBeUR75xfnkFeeRV5RHSn4Kmy9sZtWpVaw9vZZhgcMY13oczZ2b3+t/BuEOIjgID76SIi4UZ+N/w5pSwLNTw+s5KBQKJj7SjDe2nub3K+k83MzF8KTAvuDVVjuttcMzoKy5IUOFQoHr5MlYBAZx/dVXufLU0/iuWY1127YmX8PJyomhgUMZGjiU4tJi/kr9i5T8FEo0JZSWlmo/NaWUarR/bmLbhGC3YBwtjUzhBazNrLE2M+zBDPQfyKWMS3x+7nO+u/Qdmy9upkeTHoxvM54eTXrc1SstwZAYcxAefFnXiVr3MPzkQJNUN/ofaZgprfPVJXRfsI/uga58OK6T8ZNOfwObJ8HYTSANqZV2FMgXiH/hBYpTU/F+7z0chw2tlfvUhPSCdL6+8DWbzm8iJT8FP3s/QluE8mTzJ3Gzrp/9Ohqaux1zECFWePDlJHHSypKgG5Ds3XBfP1hbqBj7cFP2nEnkWlqe8ZPajADHptqUGrXESmpJwNdfYdU2hOuvvkryihVobtuytCFxtnLmX+3+xZ7QPSzotQBPG0+Wn1jOwK8HMnP/TA5fP6zfc1uoHhEchAeeJjuJWLUFzrkacgLqLxOrKcZ380ehUBB1JM74CSoz6BEBVw/Dty8a3e+hJpi5uOD/ySc4PhXKzQ8/ImH6dErzKghYDYC5ypyhgUP59LFP+XbEt4S1DuNY4jH+9cO/GLZ1GGtj1vJ74u9kFt79VODGRow5CA+8+DQZlxQlUEppy/rb4McUTZyseSzEi03HrvLygBbYWBj5X7TLJMhNgV8XwfVT8PT/gVvN94gUFhZ4z5uHVYsWJP13IVeeGYffmtWYe3vX+L1qUqBjIP/u8m+mdZzGj3E/8vWFr1l5cqW+3MPagxbOLfQ/krNES+eWVU7hbWxEcBAeeCfSzhF4Q0OxQolNm1b13Zwq/fORAL6PucHmEwmM7+ZveIJSBf3mgF9X2PIc/K8PPPEBhIys8bYoFApcnn0Wi4AAEma+wuVRo/BbvRrrdg1vxtedLFWW+gHylLwU5HSZi+kXtT8ZFzl27ph+rwsPaw96+/Wmj28funp3NToQ3tiI4CA88E7mxBF8XUOcgzeebsZnxjQkHZs6097XkU8PXeaZh5uiVFbwRNtiAEw9AN/8E76ZqH3VNOjdGlkDcSe7Pn0I2PQl155/gbjxE/B+/z0chzbcgeo7udu4427jTk+fnvpjxaXFXM26SkxqDL/G/8quy7v45sI3WKos6ebdjd6+2mDhaVv1lN4HkRhzEB54JwtSCEqEC05+NHFq+E+ECoWCf/ZsRmxKLusOxlZ+sqMvhH+v3Rjo2P/gk8GQft7fzdgAACAASURBVKVW2mXZosWtgepXXiVl5QcNdqDaFGZKMwKdAhnRfARL+y7lwOgDfDzwY0JbhPJ3xt/MOzKPAd8M4B/f/oOFvy/kQPwB8ooa7rhLTRM9B+GBll6QTk56AdYFcNHZjwmORhLbNUDD2zVhz5lE3t95Hl9nGx5vW8l7fpU5DH4P/HvA1udhTXdtLqb2Y8H/kRpdD2Hm7EzTTz4h8a23SV2zhsLYWJrMf7/SFdX3C3OVOT2a9KBHkx7Mfng2lzIucSDhAIevHyb6fDRRZ6MwV5rzkMdDdG/SnR5NetDKpdUDu65CBAfhgXYq+RRBN7RreWJdm+Jqa1HPLTKNUqlg6agOJGYe4eXoU3g6WNLJ38jCuNu1GgpTf4VfF8OZbXBqIzg11QaJ9mO0e1LXRNssLPB+710sg4JIXryYS6dO4TjiSZz+8Q8s/I2MkZiqpBhSZUg8rX015uCrTVFu56kdZ6lDCoWC5s7Nae7cnJGB4zhzI5Wf445yIuUYMTdOcixxBStOrEBZ6oCLoh1Bdp3p5N6VFu7u+Lva4Odig5V5/e40eK/EIjjhgbb0+FJKV67j8T/gxQkr+fm1AfXdpGpJy1Uzcs0hMvOL2PLCIzRzszWtojoPzn+vDRCx+wENNO0O0uNg5QjmNtqd5sytb/3Z0h7sPMDCDkycuZN7+DA3P/mU3EOHoLQU6/bBOPXtgEOwM8q8BO1Jdu5g66G9tq2b9s+2bpAZDwkn4PoJuH5Sm1zQ2GsbpRnYN9EGCkc/8O0MzfqAu2RyO6sjM6+I/ReS2XcumcOxN0nJLtSXWZkrae5hR1P3EgrMzhGbe5y00tNo/r+9M4+O66jz/ef2vqm7ta+WLdlSyba8x1s2h3hJ7DgQSBggAQLvwTADA28Ih4E3M4Q1hIEhwAzDMGceCQwQSEhiZw9eSBzHTmzL8b5cWZblRdbW2rvVe9/3R7W8yo7iyJIl1eecOtXqutVd/VP3/Vb96ldVpjCGYSLZN4lEUGCKVDCroIIbJxezeHI2s0r82CwjM8JQG+8pFAOwq6WG+06nCGX5yfMP8sZ6DZHltvHrTy/gg7/Ywqcf284zn7+BrMGMfmwumPlhmbobYe8TsOcPsP4b71zX6gJ3ruyxe/LOCkYqCUYSUol0SuJOJXEv6SBeZqN7T4Du+j007TlAsyWFtzSOtyyOK6vn8h1/ixMKZ8Lc+6F4LhTMlK/ffQp6Tsn29zTK/PhW2PekrOcpgLKboXyJFAv/hEHZdCAaAiE2HGph46FWtjd0kEwZ5Hhs3FSRy9TCDCryMpiS56HY7zwnQOAWQE5s72ndw7qG19h8ajMn3S8DL3PI0DhQn8l/HsrDlChgkq+cBcVTWTZ5FtWFOfhc1/bOsmrkoBizRBIRrv/dIn7z4wgN5VP4y73f4ycfufb2VRoMO4938LH/3kZ1kZfHP7voylwWhgHhTrlwLh6WvfRz80g3hFrlkarBVnl6XrAVgs1yJGK2gmaWLh6TJZ3M4MyULquscozMMsLNKbo27aN3w2uk+vrQHA7c82binluFZ3oxNk8SQgHIyIeiuZBbJRf3DZbOBqjfBMc2wbHX5ZoPgMwyKJoNedPSaap8Lj3nEk+mON0V5mRHmBMdfZzo6ONkRx+HmnuobwsBIPIzWDo1j2XT8pld4r90pNhlaA41s6dtD/Xd9Rxur+NgWy0tkZMYJOW/Iekk1nE93vgtiNyCM8JTkeehLNdNjtt+Re97KdTIQaG4gAPtB8hrT2CJa9Rm5FMwSiajB2LexCx++pHZfOHxt/nyE7v5j3vnvvsbiKbJY0evIhrgAlx3QCocpm/7doKvbyb4+usEt2ynBbCVl+O56SY8S6tw5U5FM79LocucBPMmwbz7peC1HiRR9xqpY2/AiZ3YDqw5c2lUc3DSPIHaVAnHYxkEDC/thpd2vPRoXmy+AspzCvjEomksm5rPhCzXe7ZBgbuAAnfBec8lUglO9p5kZ9Mh1tQ+x17zRhK8QUPsBnbvWkwonHHmWpvZRKHfQZHPSZHfSbHfwQSviWXVpWR6hj5M+VIocVCMWXa17jozGb3fN4lbR7E4AKyaUcg/rpzKQy8d4nsvHuLrK6tGzI89GExOJ54lS/AsWYJhGMQaGgi9/jrB1zfT+fjjdPzmN5izs8lYtgzvbStwLViAZrHQGYpxujtMW2+UQDBGIBhNP5YpGE3SF03QF0sSiiXoiyaJJScBkwBwEaFCO8VU8ynm2JuYajrF9ezjNmsnZiN5fiP7gBMa9EyA+grIqYDsKem8ArxFQzKvYTFZKPOVUeYr456qVRztOsqj+x/lxfoXsZW/zu0lK1mQeTexcDanusL0BJrIDGymtGU30+MHmKY1sG3fp7n+sz95z20ZdJuH7Z0UimHm7Za3mddiw2Tp44BnEvf6Rn+45WduKuNkZx+PbjnGUztPctv0Au6YWcgNU3Kwmq9dodA0DXtZGfayMrLuv59UKERw82Y6X3qFzmefo+uJJwg7PdSUzGBdzjR25VaSPGeiwmk1k5thJ9tjw++0Uux34LJZcNvMuOzp3GahwOeg0OegyO8kx2PHfO7oyjDSrrMA9AWkOyoUgN5maK+D9iPw9lsQD52tY/fJOY2K5TB5qZwUHwIm+yfz0I0P8YXZX+DX+x/jmSPPsP7E88yz+JkV6mVBVzMzojHyNCupsusIFtzFwuv/ekjee7AocVCMSVJGit1tu7m3CeyZCQKan8JRPnIAeZP91p3TWVKZy4t7m3hlfzN/2nkKv8vKbdOkUCyenH1NCoVhGBxv72NHQwc1DZ3sPuniiHs5luW3MK9FZ0XgAAuP7+KmI2+S9PpJLV1BxgfuIm/ODNz2IbhVaRo4/TJxib2oDAN6myBwRIrF6d1QtxEOpQ8uypsOU5ZKsZiwCCxXEBodbIPGGjhVQ1FjDf/YuIvPxXt53JvBFneC39gtJPJz5ds585iZW8mM3FI+YLOTfWWf/IoYtMWFELnAd4E7gGzgAPCwruvPDKJuGfB94H2AB9gJPKjr+qYrabRC8U4c7TpKX7iHzOYU1koTKUxjQhxAroFYOjWfpVPziSaSbK4N8OK+Jl7c18QTNSdx28xMzvNQnuOmLMdDea6bshw35bnugTfyew8YhkEyJQ/1NAwwMOiPcUmmDI62BdnR0ElNQwc7GjoJBGVYqM9pZU6pn9urC5g9wc+MkpXkeOykYjFCb2yh+9lnCb6whsiaJ2mpqsL/wbvwrl6NJfsq3x41TbqSvEVyxED6g7Uegrr1ULdBHtW69d/AZIWsMumGyp6cztPJ6oLOY9B+FDrqoeNYOj8qJ/pBTu7nT4PqD5Fdch1fLL6OL+ZUEjXiHGo/xP7AfvYG9rKvbR8bTmwgGAvypblfurqf/1xTDCZaSQjhAbYAAvg3oB74FLAQuF/X9f+5TN0iYBvgTtdtBz6PlO7bdF3/y2AaOiqjlVIpGYrXVisX97TpEKiFRET2QPL7UzW4h7NPMPZ5Un+S3679Dj98LIn7FjvLcx7m8HdvH9M7b0biSV6vbWNLXYD6QIj6thCnu8Oc+xMv8DqYnOemPMfD5Fw35blSPIp8zstOcEcTSY6393G0NUh9IMTR1iBH24LUt4XojSbesW0TspzMn5jFdZOymD8pk8m5nnecUE90dtLz0kt0r32WyL59YLHguflmfO9/P55blmByjJDYR3vh2GY4uU26ozrqpQgko5eu4ymQApJVBjlCrtUonC1DjkmLbHs7RjSK2e9Hc7nO+672xHrIsGZc0ff3akcr/R0wE/iQrutrAIQQjwFvAo8IIZ7WdT10iboPAkXAAl3Xd6br/g7YB/wcmDbYxl7z9H9pjm6Ek9vlF+fcRT2ubPnFcPjgyDrY/buzZRmFUijypp3Nc8VV2URtrNId7WZny052NO9gw4kNLGx3Az0Ec6RLaSwLA4DDambF9AJWTD8bKROJJ2lol0JxLBA6c0Nfu7uR3kjinLom8jIcJFMGqfRo4GwOvZE4qXNEpsjnoDzXw4fmFpPltmPSZKe738aaBhoaxZlO5k/KpPAK5nssmZlk3XcfWffdR/TIEbrWrqX7uecI/uUvmNxuMpYvx7t6Ne5FC9Esw+ght2dA1SqZ+kkl5VqM9jopFLGQFIKstCDYzq6xSXR2EqurI7JlLbG6OqJH6ogeOUKyq+vMNZrNhjkzM538WDIzsXz847jmzh22jzlYi34SaOgXBgBd16NCiJ8BvwZWAk9dWEkIYQbuA17vF4Z03U4hxH8D3xJCXKfres17+AwjRyoFLfukT/LoX+DEW5CKyyHlhIUw6UbIqZQ3+RyB4fATra0lFQxi/9BUzPRBy35oOSBT834Zt52MydfXzLK30R+3nT9dniHsL70qK0OvVQzDIJFKEE1Gz0uxZIxTvafY0bKDmuYaajtrMTCwm+3Myp3F+2MpzPa3aHRnjeow1veCw2qmqsBLVYH3vOcNwyAQjJ0Ri6NtQdqDUUwmDbOmYTZpZx6bNPC5bEzOdTM510NZjnto5gDeBfaKCvK/+lXyHniAvu3b6X7hBXr/vI7utWsx5+TgXbkS3x2rcMyY8e5DY4cCk1n+Lv2lMPnWi4oNwyD42msE/vOXRPbuPVvN48FeUUHG8uXYK6ZgcrtJdnWR7Owk0dlJslM+jhw6TOz4iWtLHIQQPqAKeHKA4h3pfAEDiAMwHTnHsO0d6o6cOKSSEOkm3tNIoOMILV3HaOk9SWuomZZIO53xEBgpMAxMGGiGgclIAQZGPExfMkqfyUSf3UNoylTCFht9GFjNERZaUlzfEKBicwuWfb8ivGcPqdDZAZatrAxHdTXO6uk4qj+LY0UVJodN9jxaD8rUclBuK3DwWSDddbP7zgpFQbXM86Zf2eTYNYRhGDSHmjnccRi9U6e2s5bDHYdpDDZe9qhHh9nBrLxZfGH2F7iu4Dpm5MzAZrZR/8sPYsmK0ZjwUjQGIpWGEk3TyM2wk5thZ1H56HFpamYz7sWLcS9eTOrBBwlu2kTPCy/S9cQTdP72t5jcbpyzZuKcPRvnnDk4Z83C7PW+8wtfJQzDIPjqqwT+4xdEDhzAWlJC7gMP4Jhahb2iAkt+/jU7oh2M/Bcj17acHKAsvXlKOsD4YvonB66k7pBhRIIc+fH9NBltNGXHabVECKRiBEgQ0KDNbKbdbMK44J9kNyALEyY0UhoYmkb/LSoFmKwuXLY8nM5scuIuKtpSFDTHyGsK420IkHXyOcwGhIGmAivB+SV45y0gO38i0QOH4HAdti2vYn/+efmaGnTnuYhMzMc8pRzftBkUzvs6uZOnY0qEoe0wNO+VI4zmfbDrd2fD7hw+qFwJU++UPRfb4BbzGMkkiUA7iZZm4i0tJFpaSbS0kAgEwDDQrNazydaf27FNmoRdVGIrLR1UT60p2ERNSw29sV76En30xfsIJ8KE4iH6En20h9up7aylJ9YDSJdEqbeUqqwqbpt0G06LE7vZflHKdmYzPXs6VvP5WxGkIhGitUfwiCjHo55xO3IYy5jsdrwrVuBdsYJkby/B1zYR3vU2fbt3E/jlf8mRPWCbMhnXnDlSMGbPxlZWhjaEO9UOhGEYBDdupO0XvyB68BDWCRMofOghfO+/E816bW+b0c9gxKH/dJSB5hT6HeqX2rTmvdQdMgLN9YT+eJC8OOQBzX44XWAmt9BNqMQDRZkUJn3kx71kRV14IzacvUlM7d0ku7rQ7HZMTicmpwPN4cTkcKA5HRjxuPQdHjlCsi1w5v3MPh/2qdNwrppFpyhkT26Ibb37qGmpoSv6NHQiZ2GKgFuhIOymOuCiskUjr7GPnNoG8rYfAzbSwU85bYX2Qhe9xX56i/2EijKJlMwlOWM5tkQER18HZV2nWXjkFbx7/4hhcZEqeR/xnOtJOCtI9IRJdrSTaO8g2R4g0d5BoqOdZKCdRHs7JC9YGGSxYMnJAZOGEY9DLI4Rj2MkEhix2HmXag4H9ooK7KISR6XAXlmBbcIELAUFNIdbWXd8HeuOr2Nv297z6pk0E26LG7fJQWHIRlHYwWJLNRNNuRRqPnJTHswdUVL7g2iOGBnvuwHnnDmDdhlEDx+GZBJnVpyWlI/5o+AcB8WVY87IwHfnanx3rgYgFQoR3reP8O7d9O3aRc+69XT9STo3TF4vzlmzcM6eJQVj1izMHs+QtCMVidC7fj3tjz5G9NAhrKWlFH7/+/juXD1qRKGfwYhDf3d6oLAm4zJl77XukJFdOoPn/vXHFLc2Ijp7mVLXwISDh4i/egrZr2+7qGHRzEwseXmY/X6MaJR4dzdGOEwqEjmTo2nYJ0/Gc+NN2Csr5U2ysgJLbu6ZoWIeMsTro8jY+yOdR2gLt5FpzyTTIdOFRxImU0ka245yet9bdB7YQ6LuKLaGZor3tZDx5ukz10UtcDpbozEbjgE9QTcFIR++ngTm+C5g13mvq9ntWLKzMWdnY83Nw1E1FUteLtb8fCz5+Vjy8rEW5GPOyrpkz8owDIxIhOjReqK6TrS2lkitTnDjX+h+6umzn8EEAS84fRqr8vx8onwBZaWzcHb1QVOAVFMz8dOnSbS0nOnhgX6mfheA2YzJ48Ho66PjV49izskhY9lSMpYvx71gwWV/bOF9+wFwZMVow0ehV40cxhMmtxv3okW4Fy0C5Pc2dqyB8K5dhHfvJrx7N4GfvyHDVE0m7ELgmjcP17y5OOfNw5qXN+j3MgyDyIGDdD39FD0vvEiqtxfbpEkU/uBhfKtXD+9k+RAymFYH0/lAfor+57qvQt0hozsc5ydbzQSjRbhsZhbOuZEbPpzDDfk2JrSfIn7iBGa/X94o8/Iw5+Zisg29/96kmRBZAoG47HVmk5nS/EpK8yvhgh2mk11dROvridbVETtaT1b9USqOHiVmJOj2mmnMirDV1kN7hkbU76DYl8HUeCtzoy0UWMxoldPlQTCVt8moi3eJpmloTifO6uk4q6djGAYHOw6y6cRr1BzcQORILXndMDWazdRoNrN6NLSGNpI7tpJiKyGTCUtBPraiYtwLFmAtLsJaXIwlvwCzNwOTx4PJ48Hs8aA5nWiaRjIYIvT6JnrWraf7uefp+uMTmH0+PLfeimNGNaneIMnubpLdXSS7u0l1dROtr8eSmYHVdZrWqJ9CvxKH8YymadjLy7CXl+G/W561neztJbxnL+Fdu+jbuZOup5+m83cygtBaWopr3jycM2dgzsnB4vdj9vtl9JDPh2a1kuzqovv5F+h6+mmihw+j2e1krFiB/+67cS2Yf9VdV1ebwYhDQzofaN14/5zCqUvUPfYe6g4ZmW4bW752K2/Wt7OlLsCWowFeffEQANluG4vKpzAJF0UxJ0VBKyWWKEV+07BHZAwGs9+Pa+7cy0YtdEQ6ePP0m2w9vZUXGrfwm4gZKKLM7GZx1y4Wv7SR654Dz5RlMP2DxCcv5XgkQF13HXWdddR11XG85zhem5eSjBKKPcVnUklGCV6bl5qWGl47+RqbTm2ita8VDY1ZubO49e4vs2LiCkq9pee1KRWJkOzqwpKd/a6H12aPG++qVXhXrSIViRB64w161q2jd8MGutfIADrN4cDs851Jrnlz8ZQCIZ02w3dFoZSKsY05IwPPjTfgufEGAIx4nMihQ/TV7KRv506Cr7565vt1ISaPByMaxYjHcUybRv6D38C3evWITn4PNe9499N1vUcIoQPzByhekM4HikYCOAz0XmHdIcXnsnJ7dQG3V8sY8KbuMFvqpFjUHO/glQPNJFPne7h8TiuTctysqi7gA7OLR82kZpYjizvK7+CO8jswDIPazlreanqLN5ve5JnmnTzuNGNGY3rvLvre3EFDjZVE2g1m0kyUZpQy0TuR3lgvbzW9RVtfG8YA3j+nxckNRTewZMISbiq+iWznpaNeTA4HpoKCS5YPFpPDQcayZWQsW4YRi5Ho7MLs8w68IGr9N0lstRK1ZJB5je+drxh5NKsV58yZOGfOJPt/fRojlSLR0kKys1OGl3Z1kUiHmSa7utFsVnyrV+OYOnWkm35VGGzX+HHg20KIO3Vdfx5ACGEDvoR02L88UCVd1+NCiKeATwoh5ui6vitd1w/8b2A/FzrGh4lCn5N75pVwzzw5gEmmDFp7I5zuCnOqM8zpLvl4X2M3D798mB+8cpgbJufwobnF3Da94JocVQyEpmnSlZUluH/6/cSSMXa37ubNpjd5u2UnJckES/qCTGmpZUpvgLKUGfuUasi7BSpvB7uHWDJGU6iJxt5GTgVP0R5uZ2buTOYXzMdmHrnwWc1mw5p/Gd9wsJUecxZFbtc1Gy6ouHbRTCashYVYCy9zfvcYZrB3uEeATwBPCCF+ApxAbp8xB/i4ruthACHEYmAysF7X9fQGIjwI3AmsF0I8AvQgt8/IAz6h6/o1cdqQ2aRR6HNS6HMy74JjcI8FQqzZ1ciaXad44Mk9OK37z4xCCn0Ostw2st12nLZr/8xYm9nGgsIFLChccH5BKikX8R1cCwefg8MvyBO6Km/DVn03EyuWM9H7Hs4HHgmCLQTwU6AmoxWKd82gxEHX9aAQ4ibgB8DnAAdwELhL1/Vnz7n0c8D9yA32WtJ1Twkhrgf+Bfha+rpdwOd0Xd88JJ/iKlOW4+aB5ZV8eVkFNcc7eebtRl7Ye5o1uxrPu85hNZHttpPltpHptpHhsOB1WPA6rGQ4LGQ4rHidFnxOK5kuKShZHhtum3nke7YmM0y6Qabb/wVOvAkHnpEH1R9cC7YMeYB99d1QfsvoWHAXbKUllaEmoxWKK0AdE3qFROJJDpzuoSMUoyMUpT0UoyMYo6MvRkcoRmdfnN5InJ5wgt5InGji0it8bRYTWS4bWW4bU/I8rKwu4BaRd22MRJIJaHgd9j8jty2OdMtjIa//Eiz6PFiv3Ruv8aMKnuip5uSND/PV26pGujkKxYigjgkdZhxWM/MmZg76+mgiSW8kQW8kQVdfjM6+GO1BKSQdfWlhCcV4oy7Ac3tO47SaeV9VLiurC7m1Km/k5jjMFrnievKtcMcjcg+pmkdh47dlvvSbcjRxrYXtpZLQF6BFRSopFFeEEodhwm4xY/eYyfHYudyi8EQyxfaGDl7e18zL+5t5aV8zdouJJZW5LJ+Wz82VueSPlA/dYgNxu0z1m2DdP8Ezn4G3fgG3fR8mLh6Zdg1EKIBmpGgz/FSPkigzheJaQonDNYbFbOL6yTlcPzmHb71/OjUNHby8v5mX9zex7qCc46/M93BTRS43VeSwsCx7ZNxP5Uvgr1+HvX+Ejd+Fx26X+zot+7bcSXakCbUCqDUOCsUVosThGsZs0lhYns3C8mweXD2Nw829bD7SxuYjAX771nF+9cYxbGYT88syuX5yDgvLsphZ4h++Q+dNJph9L0y7C978ObzxUzj8kjxGccZfyf3ubVd966yBSZ+21WaMjeNBFYrhRonDKMFk0phW5GVakZfPLZlMOJZkR0PHGbH40Z/lvkR2i4m5pZksLM9iYVk2c0r9OKxXeWRhc8GSf4C5n5RHKO77Exz5DFjdMHW1FIryW+T8xXARlCOHHksmfrUATqF41yhxGKU4bWZurszl5kp5EHlHKMb2Yx1sP9bBtmPt/GzjEQzjCDaziSUil3sXlHJzZS7mdzia8T2RUQDLvy0nqU9shb1PyjDYvU+AOxemvh+K58mzKHKrrm6kU3rkYPMWjHyYsEIxClHiMEbIctvO2x6kOxxn5/EOttS18+zuRtYfbKHI5+Aj80v5q/klV9cPbzLJU/Am3QirfgRH1sO+J2HPH6DmV/IazQw5FfL87PzpUDATiueCK2to2hBsJaw58fv9Q/N6CsU4Q4nDGMXntHJrVT63VuXztdur2Hiohce3n+AnG2r52cZabq3K42MLSllYno3naobJWuzStTR1tQwv7aiXR6M2p49HPbkN9p9ziGD2FCiZLw9gL5kvT7h7t+6oaBBO7yaAX01GKxRXiBKHcYDNYmLljEJWzijkRHsff9xxgidrTrHhkDydNdttY0KWi4nZLiZmuSjNdlOa5aLI7yDf68BqHqIJblN6tJBTIbcN7yfcJY9CbayBUzVQt0GOMkBu4TFxMSz8W6hYfvmzsw0DDqyBdf8MPY08mfiwmoxWKK4QJQ7jjNJsF/9wexVfXl7J5iNt6M1BTnT0caIjxM7jnTy/5zTnbk6raZDrsVPoc1Dgc1Doc1Lsd3JjRQ5VBRlD4893+mVobPkS+bdhQNcJOLVDpkPPw+MfhrxpcmV29d0Xb9/Reghe+io0bIaCmQRW/pJ//00vD6mtMxSKK0KJwzjFajadcTudSyyRorErzImOPpq6wjR1R2jujnC6O0x9W4itde30RhMATMp2cXt1ISurC5hZ4hu6iV9Ng8yJMs24B5Z/F/Y/DVv/Ddb+DWz8Diz+PMy9HzDgtX+Bbb+Uhxfd8QjM+xTHT/YAW9XIQaG4QpQ4KM7DZjFRluOmLOfS6xNaeyOsP9jCK/ub+X+b6/nlpqMU+RzcVl3AqhmFXDcxc2gjhCw2mP0xmPVR6XLa8jPpOtr0I1kWCsC8++HWB8Etz5Ro7o4AqDkHheIKUeKgeNfkZTi4b+FE7ls4ka6+GBsOtfLK/iZ+v+0Ej21poDzHzb0LS7l7bgmZ7iHcvVXT5LxDxXJo3Albfw7RHnjfP8lIp3No6g4DqJGDQnGFKHFQvCf8LtuZQ5OC0QSv7G/mD9tP8L0XD/HDP+usqi7gvkUTh340UTwPPvzYJYubuiM4rWZ8TrUATqG4EpQ4KIYMj91yRigON/fw+LYTrHm7kbW7T1OR5+H9s4owkGswusNxetJ5dzhOXyx55ijSC3eRz3LbWD2zkA/MLn7HTQd7InFe3NvEnw80U+h3qAVwCsUVosRBcVWoKvDynQ9U8/WVVTy/5zSPbzvBj9fXAuC2mfE6rficVrxOKyWZejGoEwAAB6lJREFULjx2M6Zzb+TnPKxvC/H9lw7zg5cPc2NFLnfPLWbFtIIzGw4mUwZbjwZ4aucpXtnfTDSRoiLPw5eWVgznR1YoxhRKHBRXFZfNwkfml/KR+aV0h+O4bOYrWjdR3xZkza5Gnnm7kf/zx924bWZWzSgkJ8PO2l2NNHVH8Dmt/NV1E7hnXsnQRk8pFOMQJQ6KYeO9+P/Lcz18ZYXgy8sq2d7QwTNvn+Klfc30xRLcIvL4xuppLJ2ah91yDZyep1CMAZQ4KEYVJpPGovJsFpVn850PVBONp/CpXVcViiFHiYNi1OKwmq/+duQKxTjlGjv4V6FQKBTXAkocFAqFQnERShwUCoVCcRFKHBQKhUJxEUocFAqFQnERShwUCoVCcRGjKZTVDNDc3DzS7VAoFIpRwzn3zHcV9z2axKEQ4L777hvpdigUCsVopBA4OtiLR5M47ABuApqA5Ai3RaFQKEYLZqQw7Hg3lTTjwv2RFQqFQjHuURPSCoVCobgIJQ4KhUKhuAglDgqFQqG4CCUOCoVCobgIJQ4KhUKhuAglDgqFQqG4CCUOCoVCobiI0bQI7ooQQuQC3wXuALKBA8DDuq4/M6ING2aEEAuBrcASXdffuKBsFvAQsAiwAW8AX9N1fd+wN3QYEEJUAN8GlgKZyIWVa4EHdV3vPue6cWMXIcQ04HvA9YALqAEe0nV94wXX3Qx8B5iDXIy6DmmT48Pb4uFFCGFG/v/n67puuaBsTNpkTI8chBAeYAPwKeAPwAPIf97TQohPjmDThhUhxGTgGQb4fwshqoHXgZnAvyK/5DOBrUKI6cPZzuFACFEMvAmsBh4DvghsBP4O2CyEcKWvGzd2EUJUAW8hheHfgX9CiuZ6IcQ951y3FPl7ygG+BfwEuA1pk4JhbvZw8w1kJ+E8xrJNxvrI4e+QP+gP6bq+BkAI8Rjy5vCIEOJpXddDI9nAq40QYjXwa+SoaSD+FSkaN+i6fjJd5ylgf7ps5TA0czj5IfLGt1jX9e3p5/5LCLEH+CnwN8AjjC+7/BCwAwt0XT8MZ34nh5Gf9SkhhIYUjiakTbrT161D/p6+CfztCLT9qiOEWAT8MxDlnHvmWLfJmB45AJ8EGvqFAUDX9SjwM+TNciz9wC9CCPE48DzQghw5XViej+zlPNV/AwTQdb0B+BOwQgiRNzytHTaWArvOEYZ+/ied3zye7JK+wcWAP/ULA4Cu60HkaGKiECILmA9MBR491/Wm6/o2YBPwUSHEmOtspr0PvwdeRtrjXMa0TcasOAghfEAVsG2A4v4NqBYMX4tGhKlIF8FcoHaA8oXp/FI2MgHXXZ2mjRgLgXsHeD43nScYR3bRdd3Qdf0eXdc/fu7zQggrctTdBXRz1qVyKZv4gcqr2dYR4t+BDOAzA5SNaZuMWlUbBMWABpwcoKwxnU8attaMDAt1XY8BCCEGKi9J5+PGRpeZJPxKOn+VcWiXftKjhGnITkUF8EVd15NCiMHa5OBVb+QwkZ5v+RRwl67rLQP8hsa0TcbsyAHwpfOB5hT60rl7mNoyIvQLw2UY9zYCEEJ8FPgs0AA8yvi2ywvAZuB2ZATXr9PPjyubpAMX/gv4la7rz17isjFtk7EsDlo6H2hPcuMyZeOJcW8jIcTHgN8if+D36LoeZnzb5RHgg8jJ+dXANiFEJuPIJul5mN8gXWp/f5lLx7RNxrJbKZjOXQOU9T/XPUDZeGJc20gI8TXgYaQd7tB1fWe6aNzaRdf1p9IP1wohaoFfAF9ifNnkAeBW4C7AIYRwpJ+3AgghcoA4Y9wmY3nk0JDOiwco6/cVnhqeplyzHEvn48pGQgiTEOIXwA+ANuB9uq5vPueScWmXAfh9Or+O8WWTO5GjgmeR34/+dD3yVLW2dNmYtsmYFQdd13sAHRludiH9UUoDRRmMJ2qQw95L2cgALgz5HNWkXQa/Qsaf68Cic0YM/YwbuwghJgkhjggh/nOA4ox0HuZshN+lbNKNXBcxFvgKsHyAtBdIpR9/hTFukzErDmkeByqFEHf2PyGEsCGHyW3I2OVxi67rTcBryHjs/p4OQohJwD3AS7qut49M664aX0FGoOxFLlw6duEF48wuJ5EL4O5Nf75z+cd0vgbZkaoH/loI4e2/QAixALgZ+J2u66PWv34uuq7v1HV9w4UJ6ASM9N87GeM2GctzDiAn1z4BPCGE+AlwAnljmAN8PD35ON75CrAF2CKE+Cmyw/D3yHj/fxjJhg016YnVb6X/XAOsHCA8sSm9n9C4sEs6TPVvkW6SLUKI/wB6kK6VFchFf3/Udd0QQnwRuahyixDil8iV5l9BCsz3RuQDjCBj3SZjeuSQXuV5E/Ak8Dngx0if4V26rv/+cnXHC7qu7wKWIF0s30X2FvcAN+m6Pirjsy/DYs6GFn4TGaV0Yfq/ML7souv6i8jPug/4OnLLjALkCPuj/b1fXddfAlYBvelrvgi8Atys63rzCDR9xBnLNtEMY9SOehQKhUJxlRjTIweFQqFQXBlKHBQKhUJxEUocFAqFQnERShwUCoVCcRFKHBQKhUJxEUocFAqFQnERShwUCoVCcRFKHBQKhUJxEUocFAqFQnERShwUCoVCcRH/H0hbr5LJQ/RFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "w_hat = fit[0]\n",
    "_w_hat = w_hat[1:] # take out intercept\n",
    "\n",
    "nosepoke_i = np.arange(0, len(_w_hat), 45 * 3)\n",
    "leverpress_i = np.arange(45, len(_w_hat), 45 * 3)\n",
    "choice_i = np.arange(90, len(_w_hat), 45 * 3)\n",
    "\n",
    "event_i = leverpress_i \n",
    "for cond, cond_name in zip(event_i, ['Con-0.0', 'Con-1.0', 'Ips-0.0', 'Ips-1.0']):\n",
    "    plt.plot(_w_hat[cond:(cond + 45)], label = cond_name)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-15 18:16:08.409580\n",
      "2019-01-15 18:18:44.530928\n"
     ]
    },
    {
     "data": {
      "application/javascript": [
       "$(document).ready(\n",
       "    function() {\n",
       "        function appendUniqueDiv(){\n",
       "            // append a div with our uuid so we can check that it's already\n",
       "            // been sent and avoid duplicates on page reload\n",
       "            var notifiedDiv = document.createElement(\"div\")\n",
       "            notifiedDiv.id = \"dbc88340-3d72-4131-b0a1-24e60dec61ba\"\n",
       "            element.append(notifiedDiv)\n",
       "        }\n",
       "\n",
       "        // only send notifications if the pageload is complete; this will\n",
       "        // help stop extra notifications when a saved notebook is loaded,\n",
       "        // which during testing gives us state \"interactive\", not \"complete\"\n",
       "        if (document.readyState === 'complete') {\n",
       "            // check for the div that signifies that the notification\n",
       "            // was already sent\n",
       "            if (document.getElementById(\"dbc88340-3d72-4131-b0a1-24e60dec61ba\") === null) {\n",
       "                var notificationPayload = {\"requireInteraction\": false, \"icon\": \"/static/base/images/favicon.ico\", \"body\": \"Cell execution has finished!\"};\n",
       "                if (Notification.permission !== 'denied') {\n",
       "                    if (Notification.permission !== 'granted') { \n",
       "                        Notification.requestPermission(function (permission) {\n",
       "                            if(!('permission' in Notification)) {\n",
       "                                Notification.permission = permission\n",
       "                            }\n",
       "                        })\n",
       "                    }\n",
       "                    if (Notification.permission === 'granted') {\n",
       "                    var notification = new Notification(\"Jupyter Notebook\", notificationPayload)\n",
       "                    appendUniqueDiv()\n",
       "                    notification.onclick = function () {\n",
       "                        window.focus();\n",
       "                        this.close();\n",
       "                        };\n",
       "                    } \n",
       "                }     \n",
       "            }\n",
       "        }\n",
       "    }\n",
       ")\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%notify\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "fit_adjusted2 = rms.robcov(fit, cluster=data_src['MouseID'].values) #, method = 'efron')\n",
    "# took about 3 minutes for CB, DMS \n",
    "\n",
    "print(datetime.datetime.now())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, ..., 7, 7, 7])"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_src['MouseID'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [1] \"coefficients\"      \"residuals\"         \"effects\"          \n",
      " [4] \"rank\"              \"fitted.values\"     \"assign\"           \n",
      " [7] \"qr\"                \"df.residual\"       \"var\"              \n",
      "[10] \"stats\"             \"linear.predictors\" \"x\"                \n",
      "[13] \"call\"              \"terms\"             \"Design\"           \n",
      "[16] \"non.slopes\"        \"na.action\"         \"scale.pred\"       \n",
      "[19] \"fail\"              \"sformula\"          \"orig.var\"         \n",
      "[22] \"clusterInfo\"      \n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(fit_adjusted2.names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "coef = np.asarray(fit_adjusted2[0])\n",
    "var = np.asarray(fit_adjusted2[8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pull the variables to put them in a dataframe \n",
    "varnames = list(np.repeat('NosePoke', 45))\n",
    "varnames.extend(np.repeat('LeverPrez', 45))\n",
    "varnames.extend(np.repeat('Choice', 45))\n",
    "varnames = np.tile(varnames, 4)\n",
    "varnames = np.insert(varnames, 0, 'Intercept')\n",
    "\n",
    "time = list(np.tile(np.arange(-1, 2, 1/15), 3))\n",
    "time = np.tile(time, 4)\n",
    "time = np.insert(time, 0, -1)\n",
    "\n",
    "cond = ['Intercept']\n",
    "cond.extend(np.repeat(['Con-0.0', 'Con-1.0', 'Ips-0.0', 'Ips-1.0'], 45 * 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_fit = pd.DataFrame(coef, columns = ['Estimate'])\n",
    "df_fit['StdError'] = np.sqrt(np.diag(var))\n",
    "df_fit['Variable'] = varnames\n",
    "df_fit['Time'] = time\n",
    "df_fit['Cond'] = cond"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_fit.to_csv(\"data/\" + recordloc + \"/multiple_event_reg_gcamp.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_hat = np.asarray(fit_adjusted[0])\n",
    "variance = np.asarray(fit_adjusted[])\n",
    "\n",
    "import pickle\n",
    "pickle.dump( w_hat, open( \"data/\" + recordloc + \"/multiple_event_reg.p\", \"wb\" ) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# gcamp regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-15 20:23:14.328191\n",
      "2019-01-15 20:23:55.731669\n"
     ]
    },
    {
     "data": {
      "application/javascript": [
       "$(document).ready(\n",
       "    function() {\n",
       "        function appendUniqueDiv(){\n",
       "            // append a div with our uuid so we can check that it's already\n",
       "            // been sent and avoid duplicates on page reload\n",
       "            var notifiedDiv = document.createElement(\"div\")\n",
       "            notifiedDiv.id = \"40ebfde8-bd81-4747-9068-1fea6ba1c178\"\n",
       "            element.append(notifiedDiv)\n",
       "        }\n",
       "\n",
       "        // only send notifications if the pageload is complete; this will\n",
       "        // help stop extra notifications when a saved notebook is loaded,\n",
       "        // which during testing gives us state \"interactive\", not \"complete\"\n",
       "        if (document.readyState === 'complete') {\n",
       "            // check for the div that signifies that the notification\n",
       "            // was already sent\n",
       "            if (document.getElementById(\"40ebfde8-bd81-4747-9068-1fea6ba1c178\") === null) {\n",
       "                var notificationPayload = {\"requireInteraction\": false, \"icon\": \"/static/base/images/favicon.ico\", \"body\": \"Cell execution has finished!\"};\n",
       "                if (Notification.permission !== 'denied') {\n",
       "                    if (Notification.permission !== 'granted') { \n",
       "                        Notification.requestPermission(function (permission) {\n",
       "                            if(!('permission' in Notification)) {\n",
       "                                Notification.permission = permission\n",
       "                            }\n",
       "                        })\n",
       "                    }\n",
       "                    if (Notification.permission === 'granted') {\n",
       "                    var notification = new Notification(\"Jupyter Notebook\", notificationPayload)\n",
       "                    appendUniqueDiv()\n",
       "                    notification.onclick = function () {\n",
       "                        window.focus();\n",
       "                        this.close();\n",
       "                        };\n",
       "                    } \n",
       "                }     \n",
       "            }\n",
       "        }\n",
       "    }\n",
       ")\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%notify\n",
    "print(datetime.datetime.now())\n",
    "Q_vals = 'Q_ch_diff'\n",
    "# create design matrix with q vals\n",
    "\n",
    "full_design_matrix = np.zeros((np.size(design_matrix, 0), np.size(design_matrix, 1) * 3))\n",
    "trial_starts = np.append(np.where(data_src['TrialStart'] == 1)[0], len(data_src))\n",
    "\n",
    "# q values, con, con * q values, ips, ips * q values \n",
    "\n",
    "# q values first \n",
    "q_val_cond = df_trial_src[Q_vals].values\n",
    "trial_len = trial_starts[1:] - trial_starts[:-1]    \n",
    "\n",
    "# have to add 0's to before the trial starts \n",
    "q_val_cond = np.hstack([np.zeros(trial_starts[0]), np.repeat(q_val_cond, trial_len)]) \n",
    "\n",
    "\n",
    "\n",
    "# block 1: Contra = 1, Ipsi = 0\n",
    "\n",
    "cond = df_trial_src.index[(df_trial_src['Action'] == \"Con\")].values\n",
    "trial_starts_cond = trial_starts[cond]\n",
    "trial_ends_cond = trial_starts[cond + 1]\n",
    "\n",
    "cond_design_mask = np.zeros(len(data_src), dtype = bool)\n",
    "for (trial_starts_i, trial_ends_i) in zip(trial_starts_cond, trial_ends_cond):\n",
    "    cond_design_mask[trial_starts_i:trial_ends_i] = True\n",
    "cond_design_matrix =  np.copy(design_matrix)\n",
    "cond_design_matrix[~cond_design_mask,:] = 0 # 0 out all ipsi trials \n",
    "\n",
    "full_design_matrix[:,  : np.size(design_matrix, 1)] = cond_design_matrix\n",
    "\n",
    "# block 2: Q values * design_matrix \n",
    "\n",
    "full_design_matrix[:, np.size(design_matrix, 1): np.size(design_matrix, 1) * 2] = design_matrix * q_val_cond[:, None]\n",
    "\n",
    "# block 3: Q values * cond_design_matrix \n",
    "\n",
    "full_design_matrix[:, np.size(design_matrix, 1) * 2: np.size(design_matrix, 1) * 3] = cond_design_matrix * q_val_cond[:, None]\n",
    "\n",
    "\n",
    "# add intercept term\n",
    "A = np.hstack((design_matrix, full_design_matrix))\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "# takes 30 seconds for DMS_CB "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_reg = pd.DataFrame(A) \n",
    "df_reg['gcamp'] = data_src['gcamp']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-15 20:23:56.153880\n",
      "2019-01-15 20:41:07.733835\n"
     ]
    },
    {
     "data": {
      "application/javascript": [
       "$(document).ready(\n",
       "    function() {\n",
       "        function appendUniqueDiv(){\n",
       "            // append a div with our uuid so we can check that it's already\n",
       "            // been sent and avoid duplicates on page reload\n",
       "            var notifiedDiv = document.createElement(\"div\")\n",
       "            notifiedDiv.id = \"b0a31916-f3ef-437d-92ab-008f36f885c6\"\n",
       "            element.append(notifiedDiv)\n",
       "        }\n",
       "\n",
       "        // only send notifications if the pageload is complete; this will\n",
       "        // help stop extra notifications when a saved notebook is loaded,\n",
       "        // which during testing gives us state \"interactive\", not \"complete\"\n",
       "        if (document.readyState === 'complete') {\n",
       "            // check for the div that signifies that the notification\n",
       "            // was already sent\n",
       "            if (document.getElementById(\"b0a31916-f3ef-437d-92ab-008f36f885c6\") === null) {\n",
       "                var notificationPayload = {\"requireInteraction\": false, \"icon\": \"/static/base/images/favicon.ico\", \"body\": \"Cell execution has finished!\"};\n",
       "                if (Notification.permission !== 'denied') {\n",
       "                    if (Notification.permission !== 'granted') { \n",
       "                        Notification.requestPermission(function (permission) {\n",
       "                            if(!('permission' in Notification)) {\n",
       "                                Notification.permission = permission\n",
       "                            }\n",
       "                        })\n",
       "                    }\n",
       "                    if (Notification.permission === 'granted') {\n",
       "                    var notification = new Notification(\"Jupyter Notebook\", notificationPayload)\n",
       "                    appendUniqueDiv()\n",
       "                    notification.onclick = function () {\n",
       "                        window.focus();\n",
       "                        this.close();\n",
       "                        };\n",
       "                    } \n",
       "                }     \n",
       "            }\n",
       "        }\n",
       "    }\n",
       ")\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%notify\n",
    "print(datetime.datetime.now())\n",
    "rdf_reg = pandas2ri.py2ri(df_reg)\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "# takes about 9 minutes for CB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-15 20:41:07.965685\n",
      "2019-01-15 20:53:34.550119\n"
     ]
    },
    {
     "data": {
      "application/javascript": [
       "$(document).ready(\n",
       "    function() {\n",
       "        function appendUniqueDiv(){\n",
       "            // append a div with our uuid so we can check that it's already\n",
       "            // been sent and avoid duplicates on page reload\n",
       "            var notifiedDiv = document.createElement(\"div\")\n",
       "            notifiedDiv.id = \"37a2f119-4a70-4c2a-82c1-6358aa26915b\"\n",
       "            element.append(notifiedDiv)\n",
       "        }\n",
       "\n",
       "        // only send notifications if the pageload is complete; this will\n",
       "        // help stop extra notifications when a saved notebook is loaded,\n",
       "        // which during testing gives us state \"interactive\", not \"complete\"\n",
       "        if (document.readyState === 'complete') {\n",
       "            // check for the div that signifies that the notification\n",
       "            // was already sent\n",
       "            if (document.getElementById(\"37a2f119-4a70-4c2a-82c1-6358aa26915b\") === null) {\n",
       "                var notificationPayload = {\"requireInteraction\": false, \"icon\": \"/static/base/images/favicon.ico\", \"body\": \"Cell execution has finished!\"};\n",
       "                if (Notification.permission !== 'denied') {\n",
       "                    if (Notification.permission !== 'granted') { \n",
       "                        Notification.requestPermission(function (permission) {\n",
       "                            if(!('permission' in Notification)) {\n",
       "                                Notification.permission = permission\n",
       "                            }\n",
       "                        })\n",
       "                    }\n",
       "                    if (Notification.permission === 'granted') {\n",
       "                    var notification = new Notification(\"Jupyter Notebook\", notificationPayload)\n",
       "                    appendUniqueDiv()\n",
       "                    notification.onclick = function () {\n",
       "                        window.focus();\n",
       "                        this.close();\n",
       "                        };\n",
       "                    } \n",
       "                }     \n",
       "            }\n",
       "        }\n",
       "    }\n",
       ")\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%notify\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "\n",
    "fmla = Formula('gcamp ~ .')\n",
    "\n",
    "\n",
    "fit = rms.ols(fmla, rdf_reg, x = True)\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "#  takes about 8 minutes for CB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-15 20:53:35.130860\n",
      "2019-01-15 20:56:42.811845\n"
     ]
    },
    {
     "data": {
      "application/javascript": [
       "$(document).ready(\n",
       "    function() {\n",
       "        function appendUniqueDiv(){\n",
       "            // append a div with our uuid so we can check that it's already\n",
       "            // been sent and avoid duplicates on page reload\n",
       "            var notifiedDiv = document.createElement(\"div\")\n",
       "            notifiedDiv.id = \"178d4ce5-1002-4957-a2a9-a8a768b2ca10\"\n",
       "            element.append(notifiedDiv)\n",
       "        }\n",
       "\n",
       "        // only send notifications if the pageload is complete; this will\n",
       "        // help stop extra notifications when a saved notebook is loaded,\n",
       "        // which during testing gives us state \"interactive\", not \"complete\"\n",
       "        if (document.readyState === 'complete') {\n",
       "            // check for the div that signifies that the notification\n",
       "            // was already sent\n",
       "            if (document.getElementById(\"178d4ce5-1002-4957-a2a9-a8a768b2ca10\") === null) {\n",
       "                var notificationPayload = {\"requireInteraction\": false, \"icon\": \"/static/base/images/favicon.ico\", \"body\": \"Cell execution has finished!\"};\n",
       "                if (Notification.permission !== 'denied') {\n",
       "                    if (Notification.permission !== 'granted') { \n",
       "                        Notification.requestPermission(function (permission) {\n",
       "                            if(!('permission' in Notification)) {\n",
       "                                Notification.permission = permission\n",
       "                            }\n",
       "                        })\n",
       "                    }\n",
       "                    if (Notification.permission === 'granted') {\n",
       "                    var notification = new Notification(\"Jupyter Notebook\", notificationPayload)\n",
       "                    appendUniqueDiv()\n",
       "                    notification.onclick = function () {\n",
       "                        window.focus();\n",
       "                        this.close();\n",
       "                        };\n",
       "                    } \n",
       "                }     \n",
       "            }\n",
       "        }\n",
       "    }\n",
       ")\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%notify\n",
    "print(datetime.datetime.now())\n",
    "\n",
    "fit_adjusted = rms.robcov(fit, cluster=data_src['MouseID'].values) #, method = 'efron')\n",
    "# took about 3 minutes for CB \n",
    "\n",
    "print(datetime.datetime.now())\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take off added intercept \n",
    "coef = np.asarray(fit_adjusted[0])[1:]\n",
    "var = np.asarray(fit_adjusted[8])[1:, 1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "varnames = list(np.repeat('NosePoke', 45))\n",
    "varnames.extend(np.repeat('LeverPrez', 45))\n",
    "varnames.extend(np.repeat('Choice', 45))\n",
    "varnames = np.tile(varnames, 4)\n",
    "\n",
    "time = list(np.tile(np.arange(-1, 2, 1/15), 3))\n",
    "time = np.tile(time, 4)\n",
    "\n",
    "\n",
    "cond = (np.repeat(['Intercept', 'Action:Con', 'Q values', 'Interaction'], 45 * 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_fit = pd.DataFrame(coef, columns = ['Estimate'])\n",
    "df_fit['StdError'] = np.sqrt(np.diag(var))\n",
    "df_fit['Variable'] = varnames\n",
    "df_fit['Time'] = time\n",
    "df_fit['Cond'] = cond"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_fit.to_csv(\"data/\" + recordloc + \"/multiple_event_reg_qvals.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'DMS'"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "recordloc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:biomath]",
   "language": "python",
   "name": "conda-env-biomath-py"
  },
  "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.5.4"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}