neural_encoding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><ul class=\"toc-item\"><li><ul class=\"toc-item\"><li><span><a href=\"#there-are-3-types-of-data:-test-data,-train-data,-validation-data\" data-toc-modified-id=\"there-are-3-types-of-data:-test-data,-train-data,-validation-data-0.0.1\"><span class=\"toc-item-num\">0.0.1&nbsp;&nbsp;</span>there are 3 types of data: test data, train data, validation data</a></span></li><li><span><a href=\"#Linear-regression-model:-B-=-W-*-r-+-c\" data-toc-modified-id=\"Linear-regression-model:-B-=-W-*-r-+-c-0.0.2\"><span class=\"toc-item-num\">0.0.2&nbsp;&nbsp;</span>Linear regression model: B = W * r + c</a></span></li><li><span><a href=\"#evaluate-model-with-test-data\" data-toc-modified-id=\"evaluate-model-with-test-data-0.0.3\"><span class=\"toc-item-num\">0.0.3&nbsp;&nbsp;</span>evaluate model with test data</a></span></li><li><span><a href=\"#build-neural-belief-data-using-the-learned-model\" data-toc-modified-id=\"build-neural-belief-data-using-the-learned-model-0.0.4\"><span class=\"toc-item-num\">0.0.4&nbsp;&nbsp;</span>build neural belief data using the learned model</a></span></li></ul></li></ul></li><li><span><a href=\"#\" data-toc-modified-id=\"-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span></a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "written by Minhae Kwon, modified by Zhengwei Wu \n",
    "Linear regression for \"Encoding\" Process\n",
    "This code works very well!\n",
    "\n",
    "ref: https://chrisalbon.com/machine_learning/linear_regression/linear_regression_using_scikit-learn/\n",
    "     https://datatofish.com/multiple-linear-regression-python/\n",
    "\n",
    "cross validate score: Coefficient of determination\n",
    "https://en.wikipedia.org/wiki/Coefficient_of_determination\n",
    "\"\"\"\n",
    "\n",
    "from pandas import DataFrame, read_csv\n",
    "#from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import cross_val_score, StratifiedKFold\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "from sklearn.metrics import r2_score\n",
    "#from sklearn.pipeline import Pipeline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "# hyperparameters\n",
    "TEST_SIZE = 0.2  # ratio of test data set \n",
    "N_SPLITS = 10 # \"K\" in K-fold cross validation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "# pdf.fontype and ps.fontype must be set to 42 in order for text in exported figure to be editable in Adobe Illustrator\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "\n",
    "plt.rcParams[\"font.family\"] = \"CMU Serif\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read data\n",
    "bb_df = read_csv('./data/bb_df.csv') #behavior belief =[belief for box1, beleif for box2]\n",
    "r_df = read_csv('./data/r_df.csv') # 100 neuron response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### there are 3 types of data: test data, train data, validation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# separate test data set\n",
    "r_df_train, r_df_test, bb_df_train, bb_df_test = train_test_split(r_df, bb_df, test_size=TEST_SIZE, random_state = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(16000, 100)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_df_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# linear regression is used for encoding process\n",
    "regr = LinearRegression(n_jobs = -1) # n_jobs = -1: parallelize taks to use all cpu cores in your machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we use only train data \n",
    "# why use StratifiedKFold?: The folds are made by preserving the percentage of samples for each class.\n",
    "X = r_df_train.to_numpy()\n",
    "y = bb_df_train.to_numpy()\n",
    "k_fold = StratifiedKFold(n_splits=N_SPLITS) # seperate train data and validation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "score [0.94607343 0.96460421]\n",
      "score [0.94210557 0.96719142]\n",
      "score [0.94606846 0.96564456]\n",
      "score [0.94840315 0.96333097]\n",
      "score [0.94583238 0.96749751]\n",
      "score [0.9484271  0.96227921]\n",
      "score [0.94490721 0.9629764 ]\n",
      "score [0.94899425 0.96383938]\n",
      "score [0.94256708 0.9634059 ]\n",
      "score [0.95023253 0.96487268]\n"
     ]
    }
   ],
   "source": [
    "# why argmax(1)? split cannot simply handle multidimension y. \n",
    "# ref: https://stackoverflow.com/questions/48508036/sklearn-stratifiedkfold-valueerror-supported-target-types-are-binary-mul\n",
    "for i, (train_index, val_index) in enumerate(k_fold.split(X, y.argmax(1))):\n",
    "    x_train_kf, x_val_kf = X[train_index], X[val_index]\n",
    "    y_train_kf, y_val_kf = y[train_index], y[val_index]\n",
    "    regr.fit(x_train_kf, y_train_kf) # fit the model\n",
    "    nb_val = regr.predict(x_val_kf) # predict based on current model -> use validation data for evaluation\n",
    "    print('score', r2_score(y_val_kf, nb_val, multioutput='raw_values')) # get r^2 score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear regression model: B = W * r + c "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  [0.04252137 0.37150612]\n",
      "Coef:  [[-7.10895602e-01  1.80012789e-01  8.67413097e-02  8.24642875e-02\n",
      "  -1.97557005e-02  1.34251865e-01  8.08297418e-02  3.60320742e-02\n",
      "   6.80391504e-02 -2.99511450e-01  3.57181318e-01  5.91749170e-01\n",
      "   3.51609437e-02 -1.88462981e-01 -1.04955899e-03  1.24486007e-01\n",
      "   1.64516468e-01 -3.27198410e-04  5.42315746e-01 -2.49562887e-01\n",
      "  -7.92409965e-03 -2.56738233e-01 -3.66067593e-01  9.23172359e-02\n",
      "  -5.41251761e-02  6.19467212e-01 -6.33234462e-02 -1.57810612e-02\n",
      "  -2.17461262e-02 -1.31455540e-02 -3.38482364e-01  2.29909202e-01\n",
      "  -2.19702129e-01  2.72408100e-02 -9.43849154e-02 -1.05292371e-02\n",
      "  -4.35298505e-01  1.05205708e-01  5.96381708e-01 -7.97842009e-01\n",
      "   8.53594591e-02  3.84171245e-02 -1.00261657e-01 -2.17585388e-01\n",
      "  -3.25154538e-01  1.45847030e-01  7.06878811e-01  2.97623408e-01\n",
      "  -2.41096452e-01 -3.05675927e-01  7.64043078e-02 -5.55676532e-01\n",
      "  -2.75102239e-02 -4.51175585e-03 -2.47226463e-01  6.43720703e-02\n",
      "  -2.05233965e-03  4.12377548e-01  2.95578363e-01 -1.85409906e-01\n",
      "   6.39528511e-03 -4.45355653e-02 -5.91337324e-02  5.25003347e-01\n",
      "  -2.20393503e-01 -2.63974012e-01 -3.24397012e-01 -7.77883047e-02\n",
      "   6.11642048e-02 -6.30355854e-02 -1.15741607e-01 -1.26742263e-01\n",
      "  -2.17058181e-01 -4.67817884e-01 -2.06638885e-01  1.02526522e-01\n",
      "  -3.23392476e-01 -3.40447925e-01  1.10629164e-01 -3.61873600e-01\n",
      "  -1.93771843e-01 -2.11482863e-02  6.58836300e-02  3.57302972e-01\n",
      "  -1.02695143e-01 -2.02651421e-02  2.68827914e-01  3.11656621e-01\n",
      "   3.29219483e-02 -3.45631682e-02 -6.94652397e-02 -1.50683249e-01\n",
      "   4.88804768e-02 -2.88700717e-01 -1.36327865e-01  1.49680476e-01\n",
      "   3.38042263e-02  1.84933428e-01  2.86284665e-01  1.44476046e-01]\n",
      " [ 2.01012313e-01  3.23023440e-01 -3.63142471e-01 -3.71824932e-01\n",
      "  -2.36265449e-01  8.82779922e-02  9.53130550e-02  4.25061492e-01\n",
      "  -2.01038416e-02  1.12033408e-01  1.06505575e-01 -2.39864538e-01\n",
      "  -2.33291823e-02 -3.54558484e-01  2.25749327e-01 -2.42462777e-01\n",
      "   5.41959249e-02 -1.13483329e-01  1.47925291e-01 -4.29210842e-01\n",
      "   5.37344272e-01 -1.00155418e-01 -3.81766461e-02  1.12139798e-01\n",
      "   6.75502013e-01  8.32131822e-02  2.74398436e-01 -3.44490337e-02\n",
      "   9.65179787e-02  2.20931985e-01 -2.21074752e-01 -5.76647505e-01\n",
      "   2.59360049e-01  5.06462046e-01  6.08578876e-02 -1.25776534e-01\n",
      "  -2.36498340e-01  1.01078434e-01 -1.88467397e-01  2.47354355e-01\n",
      "   1.60582576e-01  2.98119260e-01  4.75426231e-01  3.59610313e-01\n",
      "  -3.09261412e-01 -1.03282361e+00  1.12699158e-01 -1.20989563e-01\n",
      "  -1.05526191e-01 -2.71148564e-01 -4.28882064e-01  9.30240179e-02\n",
      "  -4.60317270e-01  1.32759566e-01  2.49995523e-01 -1.18031603e-01\n",
      "   2.33520269e-01 -1.36543320e-01  2.47358275e-01 -2.82622593e-01\n",
      "  -2.08579646e-02  7.26436446e-02 -4.15322827e-02  2.82457401e-01\n",
      "  -2.74729612e-01  2.65146659e-01 -5.24205298e-01  4.42842961e-01\n",
      "   8.82659026e-02 -1.25430435e-01  2.11672864e-01 -1.47994072e-02\n",
      "   2.60219536e-01  5.63249087e-02 -1.96580291e-01  7.91424609e-02\n",
      "  -3.25997505e-01 -3.87674253e-02 -1.28740992e-01 -2.87552917e-01\n",
      "   4.01663611e-02  3.97111536e-01 -4.17540764e-01 -1.98212043e-01\n",
      "   9.19569440e-02 -5.32693844e-02 -8.28320028e-02  5.17819811e-02\n",
      "   7.28739502e-02 -3.08474845e-01 -1.85193022e-02  2.89926579e-01\n",
      "  -2.12275817e-01 -1.34857318e-01  3.96986412e-03  4.66553427e-01\n",
      "  -4.08274893e-01 -2.84862033e-01  4.67153218e-01  3.26658944e-01]]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "c = Intercept\n",
    "B (behavior belief) N x 2 dimension (N: # of data)\n",
    "W (Coef): 2 x 300 dimension\n",
    "r(neural response - 300 neurons) N x 300 dimension\n",
    "\"\"\"\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coef: ', regr.coef_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pickle\n",
    "path = os.getcwd()\n",
    "encoding_params = [regr.coef_,regr.intercept_]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoding_params = {'encoding_params': encoding_params}\n",
    "\n",
    "# create a file that saves the parameter dictionary using pickle\n",
    "para_output = open(path + '/Results/encoding_params_twoboxCol' + '.pkl', 'wb')\n",
    "pickle.dump(encoding_params, para_output)\n",
    "para_output.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### evaluate model with test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "encoding mean squared error: \n",
      " behavior_belief1    0.002735\n",
      "behavior_belief2    0.002674\n",
      "dtype: float64\n",
      "encoding error std:\n",
      " behavior_belief1    0.052291\n",
      "behavior_belief2    0.051710\n",
      "dtype: float64\n",
      "score [0.94549446 0.96596505]\n"
     ]
    }
   ],
   "source": [
    "nb_test = regr.predict(r_df_test)\n",
    "\n",
    "encoding_error = bb_df_test - nb_test # true - estimate\n",
    "print('encoding mean squared error: \\n', np.mean(encoding_error**2, axis=0))\n",
    "print('encoding error std:\\n', np.std(encoding_error, axis=0))\n",
    "print('score', r2_score(bb_df_test, nb_test, multioutput='raw_values'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_df_test = DataFrame(nb_test, columns = bb_df.columns) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nq = 10\n",
    "belief_array = (np.arange(nq) + 0.5) / nq\n",
    "belief_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # plot: compare raw data vs prediction\n",
    "# S = 5 # dot size for scatter plot\n",
    "# plt.scatter(bb_df_test['behavior_belief1'], nb_df_test['behavior_belief1'], s=5, c='blue', marker='o')\n",
    "# plt.xlabel('Beliefs from behavior model (true)')\n",
    "# plt.ylabel('Beliefs from neural code (estimated)')\n",
    "# plt.plot(belief_array, belief_array, c='b')\n",
    "# plt.xlim([0,1])\n",
    "# plt.ylim([0,1])\n",
    "# plt.title('Encoding: Box 1 (test data)')\n",
    "# plt.show()\n",
    "\n",
    "# plt.scatter(bb_df_test['behavior_belief2'], nb_df_test['behavior_belief2'], s=5, c='red', marker='x')\n",
    "# plt.xlabel('Beliefs from behavior model (true)')\n",
    "# plt.ylabel('Beliefs from neural code (estimated)')\n",
    "# plt.plot(belief_array, belief_array, c='red')\n",
    "# plt.xlim([0,1])\n",
    "# plt.ylim([0,1])\n",
    "# plt.title('Encoding: Box 2 (test data)')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "neural_color = 'blueviolet'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'cmunrm.otf' can not be subsetted into a Type 3 font. The entire font will be embedded in the output.\n",
      "'cmunrm.otf' can not be subsetted into a Type 3 font. The entire font will be embedded in the output.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAAKrCAYAAAApj5ziAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4FNX6wPHvEAIBC80GuqIoF4UFG/yE0GtCL0pHRZoFUCxrvaAgILjIBQSuFKWoqCAiICWA1CRYuAK6ioiKMCBKE6Smnt8fs8EkJiE7O7uzu3k/z7PPSWZm57xreXP2zCmaUgohhBChoZjdAQghhPibJGUhhAghkpSFECKESFIWQogQIklZCCFCiCRlIYQIIcXtDkCIcKVpWgzQArgCuBJYp5Tabm9UItxJUhaW0zStOjAfuAko6z18DvjeWwJEARUBh/dngF5KqfeDGKq/YoD2QDPgZqAvcCEpa5p2JfAZsEcpdY8dAYrwI0lZWE4p9T1QG0DTtJ1ALeAFpdSk3NdqmhYNPAa8AlwfzDj9pZQ6ATykaVpb4NM8LqkF1ASqappWSil1Lo9rhMhBkrIItOPe8s+8Tiql0oDXNU0rjdFqDken8zm+ERgC7JeELApLkrIIFdOAt+wOwkpKqQyMzyVEocnoCxESlFLHgVS74xDCbpKUhW00TZuT69AvtgQSAJqhgqZpt2qa1kjTtFJ2xyTCg3RfCFtomlYSY/TFBUqp573nFgE1MIaZTQRmAk8BJYDy3vc9p5Tamc+9ewKdMFre6UAGMC33cDVN0y4HXgBKea8tDfwBjFdKpeRx32LAoxgP8P703vsokFccdYElwNXe328A9vn7+TRNqwa8BJwBLgc07/vrAZdgjAR5QCm1K69/NiIMKKXkJa+AvYANgMJIFFnHigEjgA35vKci0N37vlHAG8CV2c6/DewFonK9r5j33G7gpmzHnd7rq2Y7VgvYBbTMdY/ngB3A9bmOlwU2AR8AxbMdrwck5v6M2c7P8p6rbMHnq4HxULFztmP3e+/TEiNBvwRcZ/e/d3mZf9kegLwi+5UtKe/y/pzsTSwqv6TsfV8x7zW/ZE+m3nP9veea5jr+iPd4XK7jH3qPP+r9PRrwANPzqXstsCbXsXnAESA6j+ufLCApv5w7Kfvx+T7CaKFr2Y6Vwvgm8IXd/67lZc1L+pRFsIxTSjVVSsVifO1+rqCLlVKZ3h+PK6X25Dq911uWyzrgHe88FvhNKZWQ6/qpGC3WZd7fH8JodeY1thiMboeW3vHHaJp2J0aLdLEyhvDltq+gj5LnQR8/n9dNQIryZmPvfc4BvwPVCohBhBFJyiLovAnpdeBs7nOapt2Y69CPedwiw1uWznbsVowuhp/yqG+LUmqQUuqA91Ajb7k397VeP3vLBt4y1lt+l8/1/ijs5wOjm+Rqb184cKFvviywNQCxCRtIUha2UEqlY0xBvkDTtErAgFyX/uOBWzZR2X7+l7c8VIjqb73Ivc97yxreMqsV+o8/IhYo7OcDox/+d+A5TdOKeR88PovRHTQ4ALEJG8joC2EbpdTEXIfqYPQ9m5E1nO7aQly7F+PhX37D1Epluw7+bjmXNBeaZcphdMX8hdHHnQ78ClRTxpRvEQEkKYtQ0g540+R7PRgt2X/ld4GmadWVsS7HFxgLCVUn7y4Jp7f8wlt+6S2r5Hdrn6M1pymwQym1Ikj1CRtI94UICZqmXQfcB/xg5v1KqVSMr/dXaZrWMY/71wL6eH99A9gPdMjjOg1jjPM2YKH33skYDwk7e8/ndpe3DPT/T6X4u59bRChpKYtAyxpBUCavk5qmRQF3AzOAP5RSZ7zHL/Neclkeb8t60BWT6/gk770ma5r2o/JOoNA0rTzgAoYCKKX+8k4w+UjTtG5KqYXe6zRgAnAN0M7b753lEWA18B9gWLb4WwL1vb820jRtC7BXGeteZP/c15BtlIbJz7cFSPY+3DsKpGFMIjkGfKaUOprHvUSY0bKNrhHCEt71lOdhrDGctZ7yGYz1lLMebGkYiep6jHHDAAlKqXhN0z7CSK7XYYxE+B6jL3UVxljdqhjJ/pT33EMq2+w3TdPuA7oBJzH6X9OAl5VSOVaq0zStAsY44qsxuj7KA98CY5RSeY0MiQFe9MZ8HCPhbvR+trcwkuMxYKD3s43DeKiYlUTnKaWeNvv5NE0rAazAWFg/t0xgJXB/7s8pwoskZSHCgHfB/A0YXS/zlVLnvKMvYjD+SMQBw4FlSql+9kUq/CXdF0KEh8eBQ0qpGVkHvOO9z2L0w/+gadqfGK1zEcbkQZ8Q4eE4/5zhl9tlwG9BiEUEkCRlIcLDW8BpTdOmZp/Rl0XTtHYYK809GfTIhKWkT1mIMOFd36MrxlC+YxiTRzSMB5S7gbeVUgfti1BYQZKyEEKEEOm+EEKIECJJWQghQogkZSGECCGSlIUQIoRIUhZCiBAiSVkIIUKIJGUhhAghkpSFECKESFIWQogQIklZCCFCiCRlIYQIIZKUhRAihEhSFkKIECJJWQghQogkZSGECCGSlIUQIoRIUhZCiBAiSVkIIUKIJGUhhAghkpSFECKESFIWQogQIklZCCFCiCRlIYQIIUUiKbudeiO3U2/k/bmu26m3sDsmIYTIS5FIysD3wAy3U78beNv7uxBChJwikZRdHsdR4EUgCXjV5XH8ZnNIQgiRJ00pZXcMQgghvIpES1kIIcJFxCVlt1N/1O3Uldup97M7FiGE8FXEJWXgTm/5P1ujEEIIEyIxKd8FnAe+szsQIYTwVUQlZbdTLwnUAL5xeRzpdscjhBC+Km53ABarCUQD/3M79duB4UATIAZjONyTLo/DY194QghRsIhqKWN0XQBUxkjCCngL+BpoCax3O/WyNsUmhBAXFWkt5aykfDfQwOVxbM864Xbq84H7gEeBsTbEJoQQFxVpLeWskRfPZE/IXtO8Zc0gxiOEED6JmKTsdurRGAl3HzA/j0t+95YxQQtKCCF8FDFJGSMhlwCW5zPyorK33Be8kIQQwjeRlJSzui5+zed8Z2+5JuuA26lPcDv11YEMSgghfBFJSTnrIV/53CfcTr0i8BCwh2xJGagDfBn40IQQonAiZpU4t1P/EiPJ7gVqujyOM97jlwJLgGZAW5fHsdrb/3wGY0xzlu9dHkeNIIcthBA5RMSQOLdTL47Rp7wdKA3scDv1T4CSQBegEjDM5XFkdVVkAPWAbRjD5/YDKcGOWwghcouIpIwxtToGI8mOBN7A6K4A2Ar0cXkcG7Mudnkcmd4ujVPAVy6PIzK+Lgghwl5IJWVN08pjJNV1SqmlhX2fy+PYCWjZDnUpxNvuAHZKQhZChJKQSMqapt0B9AIOAEMwWryBdjtGd4cQQoSMkEjKSqnteBOkpmmTglTtbcCqINUlhBCFEklD4nxVHLjF7dQrySJFQohQUZST8otAD4wuk1dtjkUIIYAQHKesaZoC+iql5hVwzSBgEED16tXv+u472WRECGF4Jdsj/+G+pzft4pcEVli2lJVSM5VStZVStUuVKmV3OEIIYZmwTMpCiPCyqCu8UswoRcEkKQshAu6HxYDylqJAkpSFEAF3bd2cpchfSCVlTdOyFggqYWsgQghLpZ01yvSz9sYRDkJi8oimaZWBh/l7q6anNE1zAhuUUp/YF5kQwgpV28Lx3XBzW7sjCX0hkZQBHWPNi3SMFdyKYcQWWuP1hBCmfPkGpJ83ymZj7I4mtIVEUlZKZQLnsx3K8L6EEBEg7VTOUuQvpPqUhRCR6ZZ7Ac1bigJJUhZCBFyFf0FUSaMUBZOkLIQIuKRXIeO8UQbaKX7nR1YGvqIACYk+ZSFEhFO5ygD5+eefmU0rzvMnT7AXKBPYCgNAWspCiIiwfft2YmNjSeEk97GamDBMyCBJWQgRAdavX0/jxo2JiYmhP4lcx//ZHZJpkpSFEGFt0aJFtG7dmsqVK5OcnMyV3GJ3SH6RpCyECFvTp0+ne/fu1KlTh82bN3PttdfaHZLfJCkLIcKOUooRI0YwePBg2rdvz9q1aylXrpzdYVlCRl8IIcJKRkYGjz76KDNnzqRfv37MmDGD4sUjJ5VJS1kIETbOnz9P165dmTlzJi+88AKzZ8+OqIQM0lIWQoSJEydO0LFjRzZv3szkyZN57LHH7A4pICQpCyFC3qFDh4iPj2fXrl28//779OjRw+6QAkaSshAipO3Zs4dWrVpx5MgRVqxYQcuWLe0OKaAkKQshQta2bdto06YNSik2btxI7dq17Q4p4ORBnxAiJK1du5amTZtSunRpkpKSikRCBknKQogQ9MEHH9C2bVtuvPFGkpOT+de/is6an5KUhRAhZcqUKfTs2ZO6deuyefNmKlWqZHdIQSVJWQgREpRSvPjiizz++ON06tSJhIQEypYta3dYQScP+rKZNGkSJ06csDuMgCtbtizDhg2zOwwhLkhPT+fhhx/mrbfeYuDAgUyfPj3iJoUUVtH81Pk4ceIEL7/8st1hBFxR+IwifJw7d44ePXqwbNkyhg8fzsiRI9E0ze6wbCNJWQhhmz///JMOHTqQlJTEG2+8wZAhQ+wOyXaSlIUQtvjtt9+Ii4tj9+7dfPDBB3Tr1s3ukEKCJGUhRNDt3r2buLg4jh07xqpVq2jevLndIYUMScpCiKD68ssvadOmDVFRUWzatIk777zT0vtHXwZpp4wyHMmQOCGKoN3LYMbtRhlMe0igWbNmlClThqSkJMsTMkDnd+Hq24wyHElLWYgiaOMIOLwTNo2Aah2CU+dO3mMJfal5cw1Wr17NNddcE5B6qnUI3mcKBEnKPjp58iS//PILx44d49prr+XWW2+1OyQhfFa1LRzbDTe3DU59yUxiNU9wA43ZtGkpZcqUCU7FYUi6L3zk8XgYO3YsLVu2ZOHChXaHQ1pamt0hiDDkWQAZ540ykJRSPPfcc6zmCarThftYLQn5IiQp+6h+/fosXLiQ8uXL2xrHDz/8QIcOHRg/frytcYjwpIJQR3p6Ov369WP8+PHU5iG6sZBoYoJQc3iTpGyCpmlcfvnlttQ9a9YsHnjgAdasWcPy5cvJyMiwJQ4R3uInGw/D4iYH5v5nz56lc+fOzJ07l5deeon2/JdiRAWmsggjSTnMDBw4kHnz5kXs/mQiOKp1gEE7AvNA7Pjx47Rs2ZIVK1Ywffp0Xn75ZTSK7rRpX0lStll6ejrp6el5Hk9LS0MpRWpqKkoF4wunKCoCNSTuwIEDNGzYkG3btrFw4UIeeeQRaysoAmT0hR8yMjKYNm0ap06dIjU1FV3XGTx4MLfffvuFaw4cOMDEiROpVKkSSikOHjzIoEGDqF69OgDLli1jypQpbNq0CU3TmDVrFv3792flypV06tSJ0qVLExsby4cffki5cuXs+qgiwiQ8Did/NUqrWsu7du0iLi6OEydOsHr1apo2bWrNjYsYScp+WLhwIWvXrsXhcACg6zqNGzdmxowZtGzZkp07d9KzZ0+WLl1K1apVATh8+DAtW7bE7XbTqlUrunTpQpcuXRgxYgRut5v4+HgAWrZsyV133cWCBQsuvFcIq5w9lrP01+eff07btm2Jjo5m06ZN3HHHHdbcuAiS7gs/NGvW7EJCBnA4HPTu3ZtBgwaRkpJCv3796NixY46ketVVVzFw4EAGDBjA6dOnLxwfNWoUzZs3p0+fPmRmZjJu3Djee+89ScgiINLP5iz9sXLlSpo1a0a5cuVISkqShOwnScp+KFmy5D+O3X333fz66698/vnnfP3119SpU+cf19StWxdd10lKSspxfM6cOezevZtOnTpRs2bNIrUvmQiumPI5S7Pmz59Phw4duOWWW0hKSuKmm27yP7giTpKyxUqVKgUYfckAMTH/HJeZdc13332X4/iVV17J+PHjWblyZY4WuBBWaz/bGBLXfrb5e0yYMIEHHniAxo0bs3HjRq6++mrrAizCJClbLKtL4pZbbsnxe17X5G5VpKSkkJSUxNNPP02PHj04efJkgKMVRVXSOPhjp1H6KjMzk6effhqXy0XXrl1ZuXKlbeP2I5EkZT+cOXPmH8c2btzIzTffTK1atahSpco/uigAtm7dSvny5WnQoEGO46+++iovvPACo0eP5pprrmHAgAEBi10UbQe35iwLKy0tjb59+/L6668zePBg3n///Ty78YR5kpT9sGvXLvbt23fh9x9//JHFixczd+5coqOjeeedd1i8eDE7d+68cM2BAweYMmUKM2bMoEKFCheOT5kyBaUU119/PcWLF2f+/PksXbqUqVOn5ll31kw+mdEnguXMmTN07NiRd955h1deeYU33niDqCiZpWc1GRJnUuXKlVm5ciXTp08nLS2N8+fP8/vvv7N+/XpuvvlmAGJjY9m4cSPjxo3jiiuuAOCPP/5g0aJF3HXXXYAxbXrq1Kl888031KxZk5EjRwKwatUq0tLSGDp0KCtWrGDIkCG0bduWhIQEFi1axK5duwB48803+eGHH6hYsSKTJwdozqwo8o4dO0bbtm356quvmDlzJgMHDrQ7pIilhftMsdq1a6tt27ZZcq+XX3456Ds9Z2ZmomlaoXfvTUlJIT09nejoaKKjo8nIyCAtLY20tLRC9+vZ8TlFaHkl239uwy+SAvbv309cXBx79+7l/fffp3PnzgGtz2a2zweXlrLNihXzrQepZMmSOfrwihcvTvHixS+M6BDCSt999x3x8fGcOnWKNWvW0KhRI7tDinjSpyyEyFNycjINGzYkPT2dzZs3S0IOEknKQoh/+PTTT2nRogUVKlQgOTmZWrVq2R1SkSFJWQiRw9y5c+nUqRPVq1cnKSmJG2+80e6QihRJykIIwNi6afz48Tz44IM0bdqUDRs2cNVVV9kdVpEjSVkIQWZmJk899RTPPfccPXr0YMWKFVx22WV2h1UkyegLIYq41NRUHnzwQRYsWMDQoUOZNGmSz6OChHUkKQtRhKVwmg4d7iUhIYGxY8fy3HPPFXrMvAgMScpCFFFnOMK7tOXQ2v8xe/Zs+vfvb3dIAknKQhRJJ9jHPFpxkv0sWbKEDh0CsIOqMEWScjZly5YtEtOPy5Yta3cIIpf1L8LnE6Huk9BsTGDr+vbbb5lFPGmc5X7W0KFDw8BWKHwiSTmbYcOG2R2CKKI+nwgZ5+GLiYFNylu2bKF9+/bAJfRnM1dTM3CVCVPkEasQIeAa77Z2Vwdwe7ulS5fSqlUrrrnmGgaSLAk5RElSFiIEHNqWs7Ta7Nmz6dKlC7Vq1SIxMZGyVA5MRfkoWSZnKfInSVmIEFCybM7SKkopxo4dy8CBA2nZsiWfffbZhbW9g+nGloDmLUWBJCkLEQIqNwY0b2mRzMxMHn/8cV588UV69+7NsmXLuPTSS62rwAc/LgWUtxQFkqQsRAj44WNAeUsLpKSk0KtXL9544w2eeOIJ5s+fT4kSJay5uQlRMTlLkT8ZfSFEKFC5Sj+cOnWKLl26sG7dOsaPH4/L5bJ9ll50KUg7ZZSiYJKUhQgFGkZC9jN3Hj58mDZt2rBjxw7mzJlD3759LQjOf1m7zoX57nNBIUlZiBBQvioc322UZu3du5dWrVpx8OBBPvnkE9q1a2ddgH7KaqjLshoXJ0lZiBBwfHfO0lc7d+4kPj6elJQU1q1bR2xsrHXBWSDtXM5S5E8e9AkR5jZt2kSjRo0oXrw4iYmJIZeQAdLP5ixF/iQpCxHGlixZQlxcHNdeey3JyclUr17d7pDyFFM+ZynyJ0lZiDA1c+ZM7r33Xu644w62bNmCw+GwO6R8tZ8NV99mlKJg0qcsRJhRSjF69GhGjBhBmzZtWLhwIZdccondYRWoWgfjJS5OWspChJGMjAyGDh3KiBEjuP/++/nkk09CPiGDsTTp2FJGKQomSVmIMJGSkkLPnj2ZNm0aLpeLuXPnEh0dbXdYhZJ9aVJRMEnKQoSBv/76izZt2rBo0SImTJjAa6+9ZvssPV/UfRKKx8DdT9odSeiTPmUhQtwff/xB69at+fbbb3nnnXfo06eP3SH5rNmYwO+oEikkKQsRwn7++Wfi4uI4dOgQy5Yto3Xr1naHJAJMkrIQIWr79u20bt2a9PR01q9fz9133213SCIIpE9ZhI3dy2DG7UYZ6TZs2EDjxo0pUaIEiYmJkpCLEEnKImxsHAGHd8KmEcGpz64/Ah999BHx8fFcf/31JCcnc8sttwQ3AGErScoibDQZZcwKazwqOPUF+48AwJf8l27dulGnTh22bNnCddddF7zKRUiQpCzCRrUOMGhH8GaGBfOPgEKxnpf4lEdp164da9asoVy5coGrsFiuUoQMedAnRD6CNTU4IyOD5QxmGzO4k358/PEMiheX/zWLKvk7KYSNzp8/T7du3djGDBryHB2ZHZyEnJmrFCFD/hwLYZOTJ0/SsWNHNm3aRGsmUY/H7Q5JhABpKQuRj0COvjh06BCNGzcmOTmZBQsWSEIWF0hSFiIfywcYoy+WD7D2vnv27KF+/fr89NNPfPrpp/Ts2dPaCgqhZJmcpQgdkpSFyMe5YzlLK/zvf/+jfv36nDp1ig0bNtCqVSvrbu6D2oMhKsYoRWiRpCxEPqKic5b+WrduHU2aNKF06dIkJiZSp04da25sgmeBsZSmZ4FtIYh8SFIWIh8ZKTlLfyxcuJA2bdpw4403kpycTLVq1fy/qR+UrbWLgkhSFiIfWomcpVlTp06lR48e1K1bl82bN1OpUiX/g/NT/GRjYkzcZLsjEbnJkDgh8qFSc5Y+v18pRowYwejRo+nYsSPvv/8+pUqVsi5AP8ieeaFLkrIQAZCens4jjzzC7NmzGTBgAP/9739llp4oFOm+EMJi586do2vXrsyePZt///vfzJw5UxKyKDT5L0UIC504cYIOHTqQmJjIlClTGDp0qN0hiTAjSVkIi/z222/Ex8fzww8/8P7779O9e3e7QxJhSLovhLDAjz/+SGxsLHv37mXlypUhn5DXvwhjSxmlCC2SlIXw01dffUX9+vU5e/YsGzdupEWLFnaHdFGfTzQmj3wx0e5IRG6SlIXww5o1a2jatCmXXnopSUlJ3HXXXXaHVCh1n4TiMXD3k3ZHInKTPmUhTFqwYAEPPPAA1atXZ/Xq1VSsWNHukAqt2RjjJUKPtJSFMGHy5Mn07t2b+vXrs3nz5rBKyCK0SVIWwgdKKZ5//nmGDRtGly5dWL16NWXKyPqXwjqSlIUopPT0dPr378+4ceMYNGgQCxcuJCYmxu6wTAnkAv7CP5KUhSiEs2fP0rlzZ+bMmcOIESN48803iYqKsjss0zaOMBbw3zTC7khEbpKUhbiIsxynVatWrFixgmnTpjFy5Eg0TbM7LL80GWWsEtd4lN2RiNxk9IUIG7uXGS28JqOCt8LZSQ7wDvGc+GoPCxcu5N577w1OxQEmq8SFLmkpi7DxcU/jK/fHQdrS7gi7mE0sJ9nPqlWrApqQi5fKWYqiS5KyCBvpZ3OWgfTFF1/wFg3JIJV+bKJZs2YBrS/60pylKLokKQvTIvUJ/urVq2nWrBklKcMAkqjIHQGvM/18zlIUXZKUhWmR+AT/3XffpX379lSrVo2BJFOem4JSb8b5nKUouiQpC9Mi7Qn+xIkTue+++2jUqBEbN27kUq4OWt31XMZaFPVcwakvUr/lRAJJysK0ah1g0I7wf4qvlOKZZ57hqaee4t5772XlypVcfvnlwY8jiHVF4recSCFJWRRpaWlpPPjgg7jdbh555BE++OADSpYsGfQ4trqNrout7uDUF2nfciKJjFMWRdaZM2fo1q0bK1euZOTIkQwfPty2SSFRMZCZZpTBIOOUQ5e0lEWRdOzYMVq0aMHq1at58803GTFihK2z9EpVMMrSFYJTn/Qphy5JysK0cN1SSNd1GjZsyNdff82iRYt46KGH7A6JSrUBDSrWDk590qccuiQpC9PCcUuh77//ntjYWA4ePEhCQgJdunSxOyQA9nwKKPjp0+DUJ33KoUuSsjAt3LYU2rp1Kw0aNCA9PZ3NmzfTpEkTu0O6oGo7QIOb2wWnvkgZOROJJCkL0669G8pXM8pQt2LFCpo3b06FChVISkritttuszukHI7vART8ucfuSITdJCkL01Y/bvRLJjxudyQFmzdvHh07dqR69eokJSVRpUoVu0P6B+lOEFlkSJwwLWthoLQgLBBkhlIKt9vNs88+S/PmzVmyZAmXXXaZ3WHlSYaoiSzSUhamFS9tlNGl7Y0jL5mZmTz99NM8++yzdO/enRUrVoRsQobwHckirCdJWZhWs5fxoM/Zy+5IckpNTeX+++9n4sSJDB06lAULFtgyS88X4TiSRQSGJGVh2p4VxlKTP62wO5K/nT59mg4dOvDee+8xZswYJk+eTLFiof+febBHX4jQFfr/tYqQFWoPp44ePUrz5s1Zu3Yts2bN4oUXXgibvfRk9IXIIklZmHbwCzi62yjttm/fPho0aMA333zDxx9/zIABA+wOySeh9gdO2EdGXwjTksb+XTYbY18cHo+HuLg4zpw5w5o1a2jYsKF9wZgkoy9EFmkpi7CWmJhIw4YNUUqxZcuWsEzIQmTnc0vZ7dTrAtEYa3J/5fI4UtxOvSnQBygJvOvyOFZbG6YISVFAhre0wbJly+jevTvXX389CQkJ3HDDDZbeP6oUZJwzSiGCxUxL+XFgJFALKOl26vHAWqAMsAiIdzv1e6wLUYSszFxlEL399tt07tyZmjVrkpiYaHlCBih5Wc4ykGQpTZHFTFLeDrRweRzTXR7HX8B4QAe6uTyOpS6PYxhQx8ogRWjSiuUsg0Gh2Myr9O/fnxYtWrB+/XquvPLKgNQVzMkxspSmyGLmfyfl8jgyAdxO/RqgJrA465jXaSuCE6Gt3M05y0DLJJNVPME6XqBnz54sX76cSy+9NGD1BXNyjIy+EFnMJOW/sv3cEqNveWOua8qYDUj4J5hfg4/vzlkGUmpqKovpw+dMpi6P8+6771KiRImA1hnMyTGylKbIYiYpV8z2cw+MJL0p64DbqdcE9vkZlzApEr8Gnzp1inbt2vEt79OCV2nNf4IyS69qW6OlfHPbgFclxAVmxikvcTv1dUAJIBaCBk4vAAAgAElEQVQY5PI4Trmd+v8B3YGHgacsjFH4oMkoIyFHytfgI0eO0KZNG7Zv304n3uZOHgxa3Z4FRkvZs8DecdiiaPG5ueHyOHYCPYH/ADVcHsfb3lMOYCfwCHDcsgiFTyLpa/Cvv/5K/fr18Xg8LFmyJKgJGSD1TM5SiGAwNaPP5XEcAZbkOrY462e3U+/kZ1yiiPvmm2+Ij4/n3LlzrFu3jvr167PD7qCECIJAdcwNCdB9RRGwefNmGjVqRLFixUhMTKR+/fp2hyRE0JiZ0VcM4wFfJ6A0kHu+UzRQ1//QRFH0ySef0KNHD2644QbWrFnD9ddfb1ssJS6Bc0eMUohgMdNSHgfMBWoAlwBarldxk/cVFgjnmWGzZs3innvu4fbbbycxMdHWhAwQN9kYOxw32dYwRBFjpk/5eqCcy+PI9/GH26lvyu+cCKzsQ+LC5WGfUooxY8YwfPhwWrduzaJFi7jkEvubp7Jym7CDmRbtdwUlZK9JZoIR/gu3mWEZGRkMHTqU4cOH06dPH5YuXZp/Qi6eqxQiAplJyoWZRnXUxH2FRZTdARRSSkoKvXr1Ytq0aTz11FPMmzeP6Ojo/N+QnqsUIgKZScrvuJ36v91OvaAFDV82GY/wU7jM6Pvrr79o06YNCxcuxO12M2HChJDbSy+c++dF+DLzRbArcAWwy+3UfwB+I2fjLAaoZ0FswoRwmNH3xx9/0KZNG3bu3Mm8efO4//77C/W+YiUgM9UogyEc++dF+DOTlB/DGHVxBPiX95VdDMZi98IG2ffNC8VE8ssvv9CqVSt+++03li1bRps2bQr9XhXk9ZvD4Q+ciDxmkvIRoJ7L4/glvwvcTv1z8yEJf3w+ETLOwxcTQ2+9hh07dhAfH09aWhrr16+nbl3fhrPHPmN8rrufDFCAucjoC2EHM514bxeUkL3mmrivsEDdJ42VzYKVuAprw4YNNGrUiBIlSpCYmOhzQga49m4oX80ohYhUmlL+Pat3O/VoINrlcZy1JiTf1K5dW23bts2Oqou8V7S/fx5ewH9GH330Eb179+bmm29m9erVOBwOU/XNuN3o4736NmPRJSECQLv4JYFl+nG326l3djv1r4BTwEm3U/e4nfog60ITkeDNN9+kW7du1K5dmy1btphOyGCsbxwl6xuLCGcqKbud+nSMTVLB2HVkPcYIjDfdTn2BNaGJcKaU4uWXX+aRRx6hTZs2rF27lvLly/t1zz0rjP7yYOwEIoRdfE7Kbqd+L1ABqOjyOOq4PI54l8cR5/I4amKsqVzJ7dT7WB2oKJxQGFubkZHB4MGDGTlyJH379mXJkiWULu3/7qPhNltRCDPMtJQbAr28ayrn4PI4DgJxQEh/wQyFxBUodk8eOX/+PN27d+e///0vzz77LG+//XbBs/R8EEkL+AuRHzNJ+aTL48jI76TL40ghxHcesTtxBZKdrcmTJ0/SunVrFi9ezMSJExk3bhyaZvtzEyHCipmkXJj3WNM0CpBI/hqcNA7+2GmUwfT777/TpEkTEhMTeffdd3niiSeCG4AQEcJMUi7tduoN8jvpduptCPE1cYL5NTjYXSUHt+Ysg+EYP1G/fn327NnDp59+Su/evYNXuRARxsyMvrFAotupfwp8BhzDGNtXEYgHOgB3WhZhmIv09RN+42veoTUxJzNYv349//d//2d3SEKENTO7WR8FugPNgRXAViAZ+BhoAdzj8jh+tzJIqwWz9RrJXSW/sJ45NKE4MSQlJUlCFsICZnez3ul26rWBdhit4hjgf8CKQiyAb7tgtl4jdf2EhQsX8g73UYF/cT+rqVbtWrtDEiIimErKbqde0jvKYqn3FVZk9S//TJs2jaFDh+KgPr1ZRinK2R2SEBHDzOSRCRjTqsN6zeSQfhIZopRSDB8+nCFDhtC+fXseYI0kZCEsZmb0xY3AfOCn/C5wO/WypiMKgkgepxwo6enpPPTQQ4wePZr+/fuzePFioilo8xkhhBlmkvKPwFN5zejLZr7bqYfsrIFIfvgWVSpnaYXz58/TtWtXZs2axYsvvsisWbMoXlx2LxUiEMz8nzUO6Od26n8Cm4DfXB5HWtZJt1MvB1R3eRwh20MQqQ/fADLScpb+OnHiBB07dmTLli1MmTKFoUOHWnNjIUSezCTlLzFGW5QHSgO4nbqVMQl/WLjj86FDh4iPj2fXrl0sWLCAHj16+H9TIUSBzCTlGIxxyT+Q9/OyUsDD/gQVaLuXGf3KTUZFbovZX3v27KFVq1YcPXqUlStX0qJFC7tDEqJIMLtHX/+CdhrxjmEOWZE+y85f27Ztu7Ch6YYNG6hdO6T/dQoRUcw86OtViK2fHjUTTLBE8oM+f61du5amTZtyySWXkJSUVGBCLl4qZymE8J+ZadY/FnTe7dRjCPHuC5Bxynn54IMPaNu2LVWqVCEpKYmqVasWeH3pq43ykquDEJwQRYTpPfoKUAdoHYD7WiaY45TXvwhjSxllKJsyZQo9e/akXr16bNq0iUqVKl30PfGTjW8ccZODEKAQRYTPfcpup34AY0W4goT0IvfBnGad/BqodNj6GjQbE/j6fKWU4t///jdjx46lc+fOLFiwgJiYmEK9N5KHFgphFzMP+o4DXwDfkLMXoAxwF3AUYzPVkBXMZFKqPJw9DDH+7RkaEOnp6Tz88MO89dZbDBo0iOnTpxMVFWV3WEIUaWaS8jFgiMvjOJTXSe8i94VratkkmEPi2s0KzcWPzp07R48ePVi2bBnDhw9n5MiRPm/dJEMLhbCemT7lh/JLyAAuj2MlIb7IfTD7lENxs88///yTVq1asXz5cqZOncqoUaNM7aUna4gIYT3LR194WTCfLHCK8pC4gwcP0qhRI7788ks+/PBDBg8ebPpeRfmfoxCBYvmqMt4V4m6w+r5WKqoPqHbv3k2rVq34888/WbVqFc2aNfP7njK0UAhrmRl98TlwZT6no4CrAVm1JsR8+eWXtGnThqioKDZu3Midd/rfwyQzI4Wwnpk+5UsxRmAkAZsxVorbBGwAFgMdXR7HbMsiDHOhME45ISGBZs2aUaZMGZKSkixJyCDdF0IEgtkhce1cHsdfVgcTiT6fCBnn4YuJ9oxTfu+99+jbty9Op5NVq1ZxzTXXWHbvotoNJEQgmWkpvygJufDqPgnFY+DuJ4Nf96RJk+jTpw8NGjRg48aNliZkCO6u4EIUFT63lF0exxYAt1MvA9wBpLg8jq3enUZiXR5HksUxhr1gPwxTKNbyPIlPjOeee+7h3XffLfQsPV9In7IQ1jO19oXbqT8PHAA+A14G8O40Utnt1N9wO/VArKlhmWC28LJ3XwRDBul8Qj8SGc/DDz/Mhx9+GJCEDFC1rfEt4Oa2Abm9EEWSmd2s+wP9ABdQG9iadc7lcSwA3gBs+LJeeMGc9HDNHUZ59R2Br+vs2bO8T2e2M5emvBzwadN7VkD6efhpRcCqEKLIMdOivRu4zeVxvOnyOLYDZ7Kf9E4uCendrMtXBTQoV/DKlJb4fbtR/rE9sPUcP36cli1bsocVtOe/NOUlU7P0fCGjL4SwnpnRF3ohFrnPbxxzSNjzKaDgp08DX1fVdvDDYri5XeDqOHDgAHFxcfz00090YxE1uCdwlWUjoy+EsJ6ZlvLluX7P0Rzzzui7ynREQRDMERHH9wAK/twTmPvv2rWL2NhYDhw4QEJCAjU0b0IObCNZCBEgZpLyYbdTn+V26hW8v18YXOB26rcDa4FFVgQXKM3GwPPngjNuOJBf8T///HMaNGhAamoqmzZtokmTJn//25D5z0KEJTNJeQLG2haH3U79EPCw26mvdTv1n4H/Ad95H/iFrGCOvjj4BRzdbZRWWrlyJc2aNaNcuXIkJydz++23W1uBEMIWZlaJU0Ac8DRwCiNBNwFOY+xy3de68AJj9ePG6IuExwNfV/JrxpC4ra9Zd8/58+fToUMHbr31VpKSkqhSpYp1NxdC2MrUKnEujyMT+A/wH7dTLw1kujyO85ZGFkDB7G6Niob0dCgWbc39JkyYgMvlonnz5nz88cdcfnnOLv5r68HBrUYphAg/fk/ycHkcZ3MnZLdTn+bvfQPJ2QuiYowy0DJSc5ZmZWZm8vTTT+NyuejWrRsrVqz4R0IGSPOOi0m/2PgYIURIMr2esreFXBooletUNNDYn6ACbc8Ko0vhpxWBf9inRYHKMEqz0tLS6N+/P++88w5Dhgxh8uTJFCuW99/Tqm3h+G6ZZSdEuDKznnIVYAbg/wrpNgnmbtaZqTlLX505c4auXbuyatUqRo8ezQsvvFDgpJDss+xCcfdsIUTBzLSUpwJHgOeAw/xz8FUpYKSfcQVUMCc9RF8GaaeM0lfHjh2jbdu2fPXVV8ycOZOBAwde9D3B/IMjhLCemaT8g8vjKHDahdupdzYZT1AEcxfmm+KMGX03xfn2vv379xMXF8fevXtZvHgxnTp1KvR7ZYiyEOHLzIO+44W4pruJ+wZNMBck2v0JoODHTwr/nu+++47Y2FgOHTrEmjVrfErICUEc7ieEsJ6ZpHze7dQvtlr6U2aCCZZgLqQTVcIoi5Uo3PVJSUk0aNCAjIwMNm/eTKNGjXyqL/VMzlIIEV7MJOVJwANup35PXsnZ7dSjgZBepqZaBxi0Izj9yhk+POhbvnw5LVq04MorryQ5OZlatWr5XF+JS3KWQojwYqZPeScQA1wDxLidurURBUEw+5Qr1TEmc1SsU/B1c+bMYeDAgdxxxx2sXLmSK680t9Ces5exoH4wxmALIaxnJilfirGw/S7yfqZUCnjYn6ACLZjbGF1sPWWlFOPHj+f555+nZcuWLF68mMsuMzFUw0uGxAkR3swk5aPAoII2T3U79drmQwq8YA4bK2g95czMTJ566ikmTZpEz549mTt3LiVKFLLzOR8yJE6I8GYmKfcqxG7Wj5oJJliyr9wW6Jby/s2A8pbZpKam8uCDD7JgwQIee+wx/vOf/+Q7S88XsvC8EOHNzCpxuwtxTYCWdLdGMDczPXskZwlw+vRp2rdvz4IFC3j11VeZNGmSJQkZgrssqRDCeiG963SgBHPnkdyLzh85coRmzZqxbt063nrrLZ577jlL99IL5hhsIYT1imRSvvZuKF/NKAOtxOV/l7/++isNGjTg22+/ZcmSJfTr18/y+mQzUyHCW5FMysFsTRaPMcrDxb8hNjaWw4cPs3btWjp0CFzHr0yzFiJ8FcmkXLWtsZ5yMJa3TD8Pv7KF6ccboWkaW7ZsoUGDBgGrT7ovhAhvRTIpf7vAeNDnCcJOgjtPLWU+LbmUa0hOTsbpdAa0vmB2X8hDRSGsZ3qR+3B23ruk0rnCLK3kh9mzZ/OBeohK1GbQFSuoXPmKwFZIcIfEBXMSjhBFRUBaym6nHtI7xKWezllaTSnFmDFjGDhwILdd3YoHWc+tTQKfkINNHioKYb1AtZRfAuIDdG//5RqmZqXMzEwef/xxpk6dSp8+faj6wdtoRPPjUuvrsptMVBHCeqaSstupxwKdyH+PvoZ+xhVQpa6Ac0eM0kopKSncf//9LFy4kCeffBK3283EhGKcOwIly1pblxAiMpnZo+8ZYBxwBmMdjNztzRjvK2SVvxkOHjFKq5w6dYrOnTvz2Wef8dprr+FyuQCo3NhY+6JykLaSDeYKeEII65lpKccCtxY03drt1L8wH1LgHdyas/TX4cOHadOmDTt27GDu3Lk88MADF86Z2XnEH/LwTYjwZuZB39eFWP9ijplgwtHevXupX78+33//PUuXLs2RkAGKRRmlFhWceOThmxDhzUxLuTCJfL2J+wZNiXKQ+qdR+mPnzp3Ex8eTkpLCZ599Rr16/xx0kpGSsww0efgmRHgz01Je4Xbqgy5yzVQzwQRL6omcpRmbNm2iUaNGFC9enMTExDwTMuRc+yIYZEKHEOHN55ayy+P4yu3Uy7qd+jbgK+AQOR/2xQAX2fzIXrfcYzx8u+Uec+//+OOP6dWrF1WqVCEhIQGHw5HvtZdWhON/GWUwSJ+yEOHNzOiLB/i7z/jOfC4L6TVxTh0ElLf00YwZM3j00Ue5++67+fTTTylfvnyB1x/fnbMMNNl5RIjwZqZPuS/GxJDtLo/jSO6TbqdeEvjSz7gCyszoC6UUr7zyCi+99BJt27Zl4cKFlC5d+qLvK1/NSMjlq5kM1oSQ/osohCiQmT7lJJfHsSavhAzg8jhSgCX+hRVaMjIyGDJkCC+99BIPPPAAS5YsKVRChuA/6JNV4oQIb2aSckYhrplk4r5BUywmZ1mQlJQUevTowfTp03nmmWeYM2cO0dHRha4r2K1WGRInRHgzk5S3up36vRe55n0zwQRLdMmcZX7++usvWrduzUcffcTrr7/O+PHjfd66qWYvY6F7Zy+TwfqoWgcYtEMe8gkRrsz0KZ8B7nI79UXA/4Df+OfoizssiC1gUk7mLPPy+++/07p1azweD++88w59+vQxVdf22cZC99tnQ7Mxpm4hhChCzCTlpUAZQAPyG1QW1s+afv75Z1q1asXvv//O8uXLiY83v+Bd1prN5wO8drMQIjKY6b44CnQDbgRuyON1K/CTBbEFzLX1cpbZbd++ndjYWE6ePMn69ev9SsgA1ToBGvyrk1+3EUIUEWZaytuAT1weR74P/NxOfYP5kALvsCdnmWX9+vV06tSJcuXKkZCQwC233OJ3Xcf3AAr+3OP3rYQQRYDPLWWXx9GroITsveZh8yEFXtqpnCXAokWLaN26NZUrVyY5OdmShAwyGkII4RvTO4+4nboDeAyojbHQ/VfADJfH4SnwjSFo+vTpDBkyhNjYWJYvX065cn6uVJRLWHewCyGCytQefW6n3hv4HngKqIcx2mIwsM3t1B+1LrzAUihGjBjB4MGDad++PWvXrrU8IctkDiGEL3xOym6nXgsjGT8AXObyOGJcHkdJ4EqgP/C826nXtjZMixWHTDJYrj3MK6+8Qr9+/Vi8eDGlSuXe2cp/0n0hhPCFme6LQUBbl8dxKPtBl8dxDHjP7dSTgdFAbwviC4j0zPMsohe71BJeeOEFRo8e7fOkEF9I94UQorDMdF+czJ2Qs3N5HHuBU/mdt9vJkyeZr+LZxRLaFp/MmDFjApqQpftCCOELM0k5vRDXZJq4b1B8//33HFJf05X3uTvjsYDXJ90XQghfmEnKV7id+jX5nXQ79RuBIO2z4bt69erxBHupSY+g9StI94UQorDMJOWZwEq3U2/nduoX1q/07kbSA1gHTLEqwEAoTYWg1SXdF0IIX5iZPLIT+A/GSnCn3E79rNupnweOAXOBiS6PI6QXuQ8m6b4QQvjC1OQRl8fxjtupbwEex9gSKgZjxbiZLo9jh4XxhT3ZXVoI4QvTM/pcHsevwBPWhRI8l1SCM78ZpRBChBIzG6f2Ah4C+rs8jpBeDS4/Z37LWQohRKgw86CvP1ARKGFxLMETlasMoN3LYMbtRimEEBdjpvtiu8vjaF7QBW6nPhUY6vI4Cj0aTNO0rhhrMR8BKgNfKqU+NhHfxWXkKgMo++gL6VsWQlyMmaQ81+3UXwQ+LKD7op6PCbkd8DQQq5TK0DQtCkjWNO2cUmqViRhDRpNRRkKW0RdCiMIwk5SfAEoCX7id+l8Ye/Sdz3a+AuD08Z7/Bj5USmUAeBPzh8DzQFgnZRl9IYTwhZk+5dYYy3V+A+wFUjD268t6nfXlZpqmlcVYkzn3ehpHgXqapoXs7EAhhLCamZbyUaCZy+M4mt8Fbqe+yYf7VcF45HYm1/FTGPHdCOz0NUghhAhHBSZlt1N/GngEuNXlcaR6D48pKCF7TfIhhqxFjPNKygClcx1H07RBGEuIApzWNG23D/VRkTvv1NA0hVIjtK+/9uW9QoiItlop5d9uyX66WEv5NmAfOccpFGZrju98iCHrgWB0ruNZv/+ji0UpNRNjDQ4hhIgoF+tTPuzyOJrl2ii1dSHuO8nt1Au7SHFWC/myXMezfvepj1oIIcLZxVrKFdxOfTCwwOVx/Ok9Vtrt1K/DeKiXlyuAWj4MifsVSOOfSflyjHWZ9xXyPkIIEfYulpSHA58BU9xOPfvxghKlhg9LCCulTmqatgm4Otepa4DtSqnjhb2XEEKEO02pgvOn26mXAhpiJEkNGIaxdGd+SgMvuzyOfBfC/0cQxuSRV4A6Sql0TdOKA9uAkUqpJYW9jxBChLuLJuXc3E59mcvjKHA6hNupr3Z5HD49wdQ0rRtQF/gBqATsVkq971NwQggR5syMU+5biGu6+3pTpdRCYKHP0QghRATxuaUshBAicMxMsxZCCBEgkpSFECKESFIWQogQIklZCCFCiCRlIYQIIZKUhRAihAQkKbuduixML4QQJlielN1OvSSQYPV9hRCiKMh3Rp/bqXfD2IvPFxrGjtRX+ROUEEIUVQVNsx6LsRVTXkt0qnyOZ507lc85IYQQBSgoKR8FRgGbMdY1znIl8CwwFTiR6z2VMVaRe97CGIUQosgoKClvBz52eRynsx90O/X7gQddHkfuPfUAvnE79Z3Ag8CX1oUphBBFQ75J2eVxPJLPqbR8EnLW+/bL6AshhDDHzOiL3DuE5CXz4pcIIYTIzUxSvtTt1Gvld9Lt1Ksgoy+EEMIUM4vcjwc2uJ36W8BS4HeMERcVMHa6Hgp0sSxCIYQoQkwtcu926m2BtzF2rs7uJDDY5XHINk5CCGGC6Z1H3E79SuB+4E6MbpDtwDsuj+OQ6WA0rTwwElinlFpq9j5CCBGuQmI7KE3T7gB6AQeASUBfpdQ8e6MSQojgM9t9UQ54CmgAfO/yOB51O3UNmIzRWv7KdECappCkLIQoonwefeF26hWB74D+GH3IlwC4PA4FPAm0dzv126wMUgghigozQ+JeBUa6PI6KLo+jI0aCBsDlcaQDLwG9LYpPCCGKFDNJeb/L45iR7fcc/R/eFnN+ixUJIYQogJlxyudz/Z5XAr7OxH0LTdO0QcAggOrVq9/13XffXeQdQgg7ndRh23So/SiUcQS2rleyZaThvj8ys71BaaalfLPbqd+c7fccH9u7YNFpAkgpNVMpVVspVbtUqVKBrEoIYYEyDmj+auATciQw01KeAHzudurvAx6MJF0fqALcA9QD8p2GLYQQIn8+t5RdHsf3wH1AJ+C/wACMNZfnAdWAti6P4w8rgxRCiKLCTEsZl8exyu3Ub8FIzLWANOB/wCqXx5G7z7nQNE2L9v5Ywuw9hBAinJlKygDeNZXf8778omlaZeBhoKb30FOapjmBDUqpT/y9vxBChAufk7LbqVd2eRz7ch27BegG/AlscXkcO3y8rY6x5kU6kIHRrVKcXA8RhRAi0pkZfTEl9wGXx/GDy+MYBcwHbnI79cd9uaFSKlMpdV4pla4MGUqpFKVUqon4hBAibJlJyvmO43N5HCddHsdijJEYQgghfGSmTznfLgW3Uy8D1Afqmo5ICBFxTuqwbRrUHixjlS/moknZ7dQ/Av4F1Mh2LKOAt6RjDJkTQgjASMjJ4wHNmEQSaAqFZv/kPFMumpRdHse9AG6nfjuwGEgF8ttZ5ACQ5PI4dlsWoRAi7FXrDPs2QbVOga9rB++wm+XcywL8GGBmm0JH7PI4driden/gSe9DPSGEKJTk1+Dg50bZbXHg6nn99df5mKe5kaakcx64NHCVBYivf0YSMfdwUAhRhO3fnLO0WmZmJs8++ywTJkygBvdyD+9SnJKBqSzAfEqwLo8j3eVxrHc79evcTr1F7vNup/6M26mXtS48IUQkuLJGztJKaWlpPPjgg0yYMIFHH32UrnwQtgkZzO08UhdjrYsEt1Nvmuv0IuBVt1O/xorghBCR4fShnKVVzpw5Q6dOnZg/fz6jRo1i6tSpFCPK2kqCzExXRHegLTAe+DL7CZfHsRd4HBjmf2hCiEiRkeYt062757Fjx2jRogWrV69mxowZDB8+HE0LzxEX2Zkap+zyOHYBL+R10uVxpLqdukyPFkJckHrGKNMsWml9//79xMXFsXfvXhYtWkSXLl2suXEIMJOUoy9+CeVN3FcIEaHOHTbKs4f9v9f3339PXFwcf/31FwkJCTRu3Nj/m4YQM90XFdxO/f/yO+l26s2BCuZDEkKIvG3dupUGDRqQnp7O5s2bIy4hg7mW8jhgk9upvw0kAD9jrIdRFaOvuS/G7iNCCGGZFStW0LVrV6699lrWrFnDjTfeaHdIAWFm55FvMHYbeQBYBfwI7AY+BXoD/V0ex7dWBimEKNrmzZtHx44dqV69OklJSRGbkMHkRBDvSnDVgOeBD4CPgBeBGi6P40PrwhNCRIKo0jnLwlJK8dprr9G3b1+aNGnChg0buOqqq6wPMIT4s/PIMeA1C2MRQkQolZGzLIzMzExcLhcTJ06ke/fuzJs3j5Ilw3dSSGEFZMq026lPD8R9hRDhKTMlZ3kxqamp3H///UycOJGhQ4eyYMGCIpGQ4SItZbdTrwhc7/I4vsh2rC4FD4srC7QDHrUkQiFEkXL69GnuvfdeEhISGDNmDM8//3xETAoprIt1X2wGqridemy2xDwfuKmA92jI3npCCBOOHj1K27Zt2bZtG7NmzWLAgAF2hxR0F0vKnwGZwC/Zjh0FPsYYDpdXD1Ep4E1LohNCFBn79u0jLi6Offv28fHHH9OxY0e7Q7JFgUnZ5XE8nMfhr4HXXB7H8fze53bqn/sbmB0mTZrEiRMn7A4j4MqWLcuwYbI8iQgdHo+HuLg4zpw5w5o1a2jYsKHdIdnG59EXLo9jSCGu6WkuHHudOHGCl19+2e4wAq4ofEYRPhITE2nfvj2lSpViy5Yt1KxZ06/7adGg0owyHAVq9IUMlRNCXNSyZcto2bIlV111FcnJyX4nZIBi3qZmlCRlg9up34Yx+kIIEaJO6vDZc0Zpl7feeovOnTtTs2ZNEhMTueGGGyy5b4b3SWdTy/AAACAASURBVJeVy4QGU77dF26nfgCoaPK+R02+TwgRBMHeXTo7pRTjxo3jhRdeoFWrVixevJhLL7VuL73r68H+TXB9rGW3DKqC+pSPAl8BOzFGYGTRgPswhsvlfipWGbgdY3SGECJE1R4MaFA7yLMJMslk2LAnmDJlCr169WLOnDmUKFHC0jqurG4k5SuqW3rboLlYUh7s8jh+y37Q7dQHAve5PI6teb3J7dTvB05ZF6IQwmplHMFvIaeTyhIe4NspHzBs2DBef/11ihWz/rFWlXjY/QlUibP81kFR0D+Rh3InZK8y+SVkAJfHMR+IvEVOvU6ePMn27dtZt24du3btsjscIcJCCqd4l7Z8yweMGzeOiRMnBiQhA2x43tgLcEOeeyOFvnz/qbg8jp/zOVWY3arPmAsn9Hk8HsaOHUvLli1ZuHChLTHMmjWLYcOGce+999KwYUOmTp1KZmbmxd8ohFcwH/QdPnyYOTTjVzbQibd59tlnAzptOuvW4Toz28wqcZXcTj3a5XGk5XXS7dQ14Gr/wgpd9evXJzY2liuuuMKW+kePHk3VqlWZNGkSANu3b6dJkyZs3bqV9957z5aYRPgJ1oO+vXv3EhcXx2F0erCEW2gfuMq82r0Fa4ZBq0kBryogzHx/SACWup36DblPuJ36VcACjIeDEUvTNC6//PKg1/vHH3/w448/0r179wvH7rjjDgYMGMCCBQv47LPPgh6TCE+1B0Psc4F90PfNN99Qv359jh49Sl8+C0pCBrjubui31SjDkZmdRz4E/gB+djv1HW6nvtrt1Fe5nfqXwEGglMvjeMPqQAV89dVXfPDBB0yYMCHH8dhYY+zP5s2b7QhLhKsALhu2efNmGjVqRLFixdiyZQvXE7zxaaEwBtsfZnceeRBjS6gooDnQErgUeAzobFl0RUB6ejrp6f8c5Z6enk5aWhpKKVJTU1FKUalSJZRS7Nu3L8e1UVFRAJw6JYNeROFkdV9sC8DK55988gmtWrWiYsWKJCcnU6NGDesrKUAgP1sw+LPzyBxgjtupxwCay+M4Z11Y4SEjI4Np06Zx6tQpUlNT0XWdwYMHc/vtt1+45sCBA0ycOPFCQj148CCDBg2ienVjEOWyZcuYMmUKmzZtQtM0Zs2aRf/+/Vm5ciWdOnWidOnSxMbG8uGHH3LnnXdy8uRJSpfOuafO119/DUCdOnWC9+FFWIupAFpxiClv7X1nzZrFww8/TJ06dVixYgUVKgR/Y3u7xmBbRill6vVajf3/eq3G/r6v1djfyfu79lqN/T1fq7E/yuw9zbzuuusuZZWXXnqp0NfecMMNqlq1amr//v0Xju3fv1/deOONas2aNUoppXbs2KFuvfVW9eOPP1645o8//lC1atVSCQkJOe43fPhwFRMTow4cOKCUUurs2bOqdu3aOd6bl/Pnz6vKlSurGjVqqPT09ELF7svnFJFpTCmlRqHU2NLW3C8zM1O98sorClCtW7dWp0+fznF+FH+/QlzQcld+L1PdF26nPhP4HngbeATA5XEoQAdmeFvPEa9Zs2Y4HI4LvzscDnr37s2gQYNISUmhX79+dOzYkapVq1645qqrrmLgwIEMGDCA06dPXzg+atQomjdvTp8+fcjMzGTcuHG89957Od6bl5deeon09HSWLl16oRtDiIu54lZveYv/98rIyGDo0KEMHz6cPn36sHTpUi655BL/b2xSketTdjv1Z4HqGIsOVQDWZ51zeRyJwL+BZ6wKMJTltWfY/7d332FSFOkDx79DTpJUFLEPBUxnk5SfIiBJZZclSZAT8NADRQkKp7QIBsBDAdsACggqIIKKBOEQkCWzsIsJUa6RoGJokgGRKLDs9u+P6tHZYePsxJ338zw8xfb0VNWkd2qqK9x44418//33fPTRR3z++efZdik0btwY27ZJTU3NcnzmzJns2rWL22+/nbp163LllVfmWv5HH33EO++8w7p166hdO7fNYITI6revs6aBOn36ND179mTy5Mk88sgjzJo1i5IlI7s827rHVZ/yuscjWo2ABdJSvhRobljaCsPSDpN1XQwMSzsIBHcyewwpW7YsoPqSAcqUOfdHg/ec7du3Zzl+4YUXMn78eJYvX56lBZ6dffv2MXDgQFauXJlna1oIfxlnVHr2TOB5HD16lKSkJObNm4dpmjz//PMhm6VXED//T6W/WJGtR6ACeQYPGpaW1/SxSwKpTFHg7ZK4+uqrs/yd3Tn+rdvTp0+TmprK0KFDufPOOzly5Ei2ZZw4cYI+ffowe/bsP8sBddFQiPxo/AjgcdMA/PTTT7Rq1YoNGzYwa9Yshg4dGtT6FUbSVKjRGNq+GumaBCaQoHyRO2vPK8tkRlO3NfI3FTvmnThx7mzy9evXU6dOHerVq0etWrXO6aIA2Lx5M1WrVqVZs2ZZjo8dO5YRI0YwZswYLr744mw3jczMzOTBBx/kueee+3MEB8A333zDN998E4RHJeLBjvmAAzsWFPy+e/bsoWnTpuzYsYMlS5bQu3fvoNevMGJ98kggQ+K+AD40dftxYDvuEHRTt0sDbYHngVgdjFIgO3bs4IcffqBmzZoA7N69m4ULF/Luu+9SsmRJZs+eTffu3enTpw/169cHVLfGyy+/zLRp07IMF3r55ZdxHIe//e1vALz11ltce+21TJo0iUGD/tqB69FHH8W2bebMmcOcOXNwHIejR4+yfft2nnjiiTA+ehHLjrtLjR3fV7D7ffHFFyQmJpKens7atWtp3Lhx8CsX5wLZo2+GqdvNUGstZwJHTd3uh1oQvwzwrGFpK4NbzehTs2ZNli9fzpQpU0hPT+fUqVMcPHiQtWvXUqdOHUDNtFu/fj3jxo37c62Mn376ifnz53P99dcDalznpEmT2LZtG3Xr1mX06NEAfPjhh6Snp/Pggw+ybNkyBg0axPnnn88LL7wAwOrVq8+p0zXXXBOOhy6KgPSTWdP8WL9+PZ06daJSpUqsW7dO3m8h4nGcwOZamrrdGdUirgekA1uACYalrQte9fLWqFEj57PPPgtKXqNGjQr7pqKZmZl4PJ6QrprlLxKPU0SX//i83Z7MRwhYuHAhPXv2pE6dOqxYsSLPC9GFLS+CIr62XGFm9C0CFgWxLnEpGq5WC5GbqVOnMmDAAG666SY++OADqlYN8jRAkUUg45Rbm7o9y9TtQPfvE0LEAMdxGD16NP379ycpKYlVq1ZJQA6DQJppo4FuqPHKQogiKCMjg4EDBzJq1CjuueceFi1adM6aKyI0AgnKm4HzDEv7NKcTTN1+MvAqCSEi6dSpU/zjH//g1VdfZdiwYcyYMSPis/TiSSBBeSkw2NTt3PYEbxdgfYQQEXTkyBHatm3LwoULefHFFxk3blxYL0KLwC70tQIqAbtM3d4N7AdO+dx+PtAwCHUTQoTRwYMHadu2LZZlMWfOHHr16hXpKsWlQILyAKAc8CtQ0/3nq0yA+QohIuSbb74hISGBn376iaVLl5KQkBDpKsWtQILnL0CSYWk/5nSCqdsfBV4lIUQ4ff7557Rt25aMjAzWrl3LDTfcEOkqxbVA+pSn5xaQXW8GkK8QIsz2sJaWLVtSpkwZUlNTJSBHgUA2Tn0pH+dMDaw6QsSnSCzMbjGf2bSlZs2apKWlcdVVV4WsrLLVVFquWsiKKDIC7vs1dTsRGIS6qOeg1sJ4ybA02VJZiALybvaJB24ZG/ryPmEKyxiERlNSUpZQpUqV0BcKODKQI08BBWVTt18BBqKC8QFUi7sT0MHU7VGGpY0JXhXDp3LlynGxJkTlynGxsmpMuaozfL8Brro9tOU4jsPIkSNZyn+4mo7cwVyqVCkb2kKBPw656a8hLyrmFTgom7rdA+gCPADMMSztpHv8PKAPMMLU7U2Gpa0PZkXDYciQIZGugohTuxbB/o9g1+LQrQOckZHBgAEDeO2117iOPnRgGsXDNFCqyTBIG6tSkbtAXpHuQAvD0rKsqG5Y2jFgoqnb64HHgPWFrp0QcaLRQMADjUK0EvmpU6fo2bMnixYtYsSIEZR+dgyeMC6Idvow4MDp38NWZMwKdDuoHLe4MCztS9SEEiFEPh3bD9+vU2mw/f777yQkJLBo0SImTpzIM888E9aALAomkKB87hbO5/Kd4Yep2xUDKEeIuJE8GPZ/DCuD3IN24MABWrRowebNm3n33Xd56KGHgltAPtW7W+2bVy+6do6KSoEE5a9M3b4+pxtN3a4D+O8WOjeAcoSIGw3vg5LlocG52zIG7Ouvv6ZJkyZ8++23LFu2jDvvvDN4mRfQZ5Nh30cqFbkLpE/5G2C6qdvTODf4XoJajGimqdv/dI+VRtbCECJXn7wE6SfgkwlwXd/C5/fZZ5+RlJSE4zisX7+eRo0aFT7TQtj3cdZU5CyQoDwTtSBRbt95zfz+ju4NYISIMO+ubAHuzpbFqlWr6NKlCxdccAHJyclceeWVhc+0kIq5kaa4rACap0C6Lw4BicDlwGX5+Hc18HUh6ihEkddqLFSoDq2eLVw+c+fOpV27dtSqVYvU1NSoCMjw1+NrWcjHFw8CaSl/WtDdqk3dDutmqkLEmj0r4PgB2JMMV3cMLI+XX36ZwYMH07x5c/773/9G1SShfWnq8e3bHPjjixeBrH3RI4D79C/ofYSIJ4e/c9M9Bb+v4zg8/vjjDB48mM6dO5OcnBxVARnUjMVLGod+xmJRIFspCxEF9m7OmubX2bNnue+++3j22Wfp168f8+fPp0yZMsGvYCFtm6VmLG57K9I1iX4SlIWIAhlnsqb58ccff9C1a1emT5/Ok08+ydSpUylevHhoKijCRnYIESIKZKa76dn8nX/48GE6duxIamoqkyZNYuDAgaGrXBA0HQ6lK4VuGnlRIkFZiCjgHQmXnyFx+/btIzExkd27d/Pee+9xxx13hLRuQSMDY/NFgrIQ0eCsX5qDnTt3kpCQwOHDh/nwww9p3bp1yKsWDKljYcurcPoYJMmsvlzl2qds6nZJU7cvD1dlhBA5++STT2jWrBmnTp1i/fr1MROQRcHkdaFvGrDL1O0K4aiMENEkEls05SQ5OZnWrVtTqVIlUlNTue666yJdpQKpd7caEicLEuUtr6B8AhhjWNqfa1yYup3njw9TtzsXtmJCRJp3i6bPpkS2Hm+//Tbt27fniiuuIDU1lTp16kS2QgHwXcRf5C6vPuVyhqU97Xfs0nzk+xCwKLAqCREdQr3wfH5MmDCBf//737Rq1YrFixdTsWJsroIbDc9lrMgrKB8wdXsB8Dbg3V1LM3X7LshxlewLkFXhRBFQSQvPJqbZcRyH4cOHM378eLp27cqcOXOiclJIfkXyuYw1eQXl/6BWhZsL+K7vlNe8HBn8IkSAvLP03nzzTR544AEmTZokk0LiSK5B2bC000BPU7fLoVrAHmA6aoPUnJQDlgSthkJEyBFb9Ss3GqhaeuFwhpN07vwPli5dyqhRo3jqqafweGTrpniSr3HK7o7VPwKYuv2tYWk/5na+qdtrg1A3ISLKe6EPT3h+ep/kN96mA3uXbebVV1/lgQceCH2hYRKJL7hYVeDJI4al3Z+Pc4rOu0nErXBenDrCXt4igd/4hvnz59O1a9fQFxpGMnkk/wKa0WfqdgngHtQOI98Zljba1G0P8AwwNa+WtBCxIFwXp3bs2MEbJHCKI/Qmma5dW4a8TE9pcE6rVESXAq8S5+5M/QUwBTXK4u8AhqU5wHPAo6Zu1wxmJYUoqj766COaNWtGBmfowwYup2VYyi1RKmsaak2HQ5PHoOlj4SkvlgWydOd44AOgqmFp9YEt3hsMS/sdeBiQ0Ygi5oV6Rt/y5ctp3bo1VapU4V7SqE6D0BSUjfRjWdNQ8/7qkP7kvAUSlH83LG24zyy/LMPfDEs7A2QWumZCRFgoZ/S99dZbdOzYkWuuuYbU1FSqUiv4heTC4w5wLSYbmUadQIKy/3drduN15PtQxLyrOkONEGxh9Pzzz3P33XfTsmVL1q1bx0UXXRTcAvKhVDmVliwX9qJFHgIJyrVN3a7q83eWlrKp27cWrkpCRIddi2BfENdryMzMZOjQoRiGQffu3Vm2bFnEpk3XuClrKqJHIKMvZgJfmro9AbCA803drgHUAroCvYDGwauiEJFR5QooWR6qBGH9n/T0dPr27cvs2bMZNGgQEydOpFixyO3GVsVdkLdKeHtNRD4Espv1JmCk+285YKAmlqwHegD3GJb2bRDrKERErBsB6Sdg/eOFy+fEiRN06tSJ2bNnM2bMGF5++eWIBmRQS2nWkKU0o1JA7wzD0maghsI9BSwA3gWGArphacuCVz0hIkdrptJLmwaex6FDh7jllltITk7mtdde4/HHH4+KadPB7poRwRPwdlCGpe1FTRYRokiqkwR7klUaiB9//JGEhAS+++47Fi5cyO23B/mKYSHIUprRKyS/oUzdlrWURcxbO0x1X6wNYMLD9u3badKkCQcOHGDlypVRFZBBxg1Hs1xbyqZu1weuNiztPZ9j3YCyudytMtAkONUTInKKl8qa5ldqairt27enTJkypKSkUK9eveBXThRZeXVfLAEuNXV7v2FpG91jY4AryHmRe5D1lEURcOakStNP5P8+H3zwAd27d0fTNJKTk7n8ctl3WBRMXt0XrwOrUWtdeP0CDAFqA5dl8+8a4OtgVlIICP9Gpk6GSjPzOT915syZdO7cGV3XSU1NlYAsApLXIvdjsjmcDLxpWFqOs+ZN3f6wsBUTwl+41zdOdxcSyGt9CMdxGD9+PMOHD+e2225j4cKFnHfeeaGvoCiSAllPObtA7X/OvwOrjhA5i8YRA5mZmTzyyCNMmDCBHj168Oabb1KqVJiWXhNFUqhGX7wfinyFiKarFWfOnOGf//wnEyZM4KGHHmLOnDkSkEWh5dhSNnW7OhDIBNPSwHUB10iIHIS7+yI3x48fp2vXrqxcuZKxY8cybNiwqJgUImJfbt0XbYAZAeTpIaraM6KoiJbui19++YV27dqxZcsWpk+fTp8+ue0jLETB5BaUfwW2A48ApwuQZ1lgWmEqJUR2wrU9U26+//57EhIS+PHHH1m0aBEdO3aMbIUCJBuZRq/cgvIvwCrD0lYVNFNTt7cGXiUhshfpQLJt2zYSExP5448/WLVqFc2aNQt/JYIkmrqCRFY5BmXD0j4BPgkkU8PSOgdcIyFyEMlAsnHjRjp06ED58uVJSUmhbt264a1AkEVLV5A4V8ALErmL2TcDDhqWNtXdzfrfwOu5jWEWIlA1mkD56uFfmH0Hi3nmtju57LLLSE5OpmbN4O8LfFED+OkLlYZDNHQFiewFspt1SVO3VwIrUTP7kuDP3az/C7zktzOJEEGRMhJOHICNo8JX5hbeYC5dqV+/Pps2bQpJQAao7C42X7l2SLIXMSSQccpjgEPA3w1Lqwxs8t7gLm7/MOrioBBBVf7irGkoOY7DBp7hv9xHbW5jzZo1XHDBBSErb/+nKj0QUIehKEoC6b4oblhaD5+//XezPmrqYVqcQMSVoz9mTUMlMzOTwYMHs4ZJ1KMntzOTChVCOynk2AGVHj0Q0mJEDAikpXw4H+dcGkC+QuSq+vUqvTiEU5NOnz5Njx49mDRpEjfxb7owmxKEfpaed95JMZl/EvcCCcqaqdulff7O8jZy12CWjctF0F3dDSpUV2koHDt2jHbt2jFv3jzGjx9PIi9QLDQrEZxDc1cgv1RWIo97gbzjlgOfmrrd1dTtq3C7QEzd1kzdHgKsAMYFsY5CAOpC3/EDkDIq+Hn//PPPtGzZkvXr1zNz5kweffRRPLkuGR5cF/7dTa8NW5EiSgWyStwSU7cbAu+hWsmZpm6PRgX4U8AQw9K2BLeaQkBxd7+bEmWCm++ePXtISEhg3759LF68mPbt2we3gHyodzcc2Cq7S4vAd7MeDdyIWhvjc+BjYBJwvWFprwevekL85YD7Vb8/iF/5X375JU2bNuXQoUOsXr06IgEZ1O7S+2V3aUHhdrPeAtwXxLoIkaviJeDsGZUGw4YNG+jYsSMVK1ZkzZo1/P3vfw9OxgG4qjP8sAGuiq79VUUEBOUqhqnb1UzdviIYeQmRkz83Mi2d+3n58f7775OQkECNGjVIS0uLaEAGNYV830cqFfEtX20OU7ebAO2AY8Bsw9L2ucfLA28A3d2/vwG6GZb2v9BUV8Qz/S7YMgn0XoXLZ9q0aQwYMIAbbriBpUuXcv755wengoVw8POsqYhfebaUTd1+EkgBHgOeAb40ddsdMcok4B+oC34e1C7XG0zdDt3UJxG3Pn8ta1pQjuPw9NNP88ADD5CYmMjq1aujIiDDX2OvQzkGW8SGXIOyqdvNgNGodZUnAm+6N71j6vbVQCfgIaA6cCGqjzkTGBGi+oo45pzJmhZERkYGgwYNYuTIkfTu3ZvFixdTvnz54FawEBoNhBqN3dXbRFzLq/vCAMYYlvaU94Cp2wZq4aFngX/4rbc83dTtH4DxQa+pEAE6ffo0d911FwsWLMAwDMaPHx91Wzdtm6X6lLe9BZfeGOnaiEjKq/uimm9ABjAs7TegJ2r42zkL4BuWtppCjOoQIpiOHj1K27ZtWbBgAc8//zzPPfdc1AVkIXzlFZR/ye6gYWk2auZeTn4IuEYiZhyxYc1jKo1GBw8epEWLFmzcuJHZs2fzyCPRu3hhvbvhksYyeUTkHZTTc7ltb4D3E0WEdyeQz6ZEuibn+vbbb2natCm7d+9myZIl3HXXXZGuUq5k8ojwKkw3g/wGjHPRuqXQ1q1bSUxMJCMjg7Vr13LjjdHfSRutz6UIv7xayhfmcpuTy21VAqiLiDHH9sP361QaLdauXUuLFi0oXbo0mzZtiomADH9tzyQ7S4u8WspNTd3+GjiZzW3VTN3ObhHFMoBsahMHkgfD/o9h5RDosznStYEFCxbQq1cvrrjiClasWMGllxZyWe/iQIabChEmeQVlD7kH2ItyOJ5bK1oUEQkTVUBuMyHSNYFXX32VgQMH0qRJEz744AOqVCn8j7XSFeD0EZUKES55BeW9QH3D0vKz2wgApm6fj1o5ThRx510CNVuoNFIcx2HUqFE8/fTTdOjQgblz51KuXHD2WKhxE+xZIQvPi/DKq095R0ECMoBhaYeAHYFXScSK1LFq9EVqhLY0yMjIoH///jz99NP06dOH999/P2gB2Zf87BPhlFdQHhpgvkaA9xMx5MxxNz0W/rJPnTpF9+7dmTZtGsOHD+eNN96gRIngzlnam5Y1FSIccn0XB7ram6wSFx9KuX2tpc4Lb7mnOEJiYic2bNjAhAkTGDx4cGgKkkGfIgJkOrQIWNPhULpSeMfWHuMAs2nLobSveOedd+jRo0fIyrq4Afy4Aao3DFkRQpxDgrIonDB2uB7ia94igRP8zNKlS2nTpk1Iyzt9JGsqRDiEZ/90USSte1xd6Fv3eOjL2rJlC2/QlNMc41+sC3lABmg0CEqWh+tlOU0RRhKURcD+3C1ja2jLWb16NS1btqQk5biXVGrwf6Et0LX1dUg/AV+8EZbihAAkKItCOM+dEnxeISfO5ea9994jKSmJyy+/nHtJ4wKuDF1hfpqOgArVocnwsBUphARlEbhDO910V2jyf+WVV+jRoweNGzcmJSWFioR3lsrOBXD8AOxcGNZiRZyToCwCdnSfm+a2iGsAHMfhiSee4KGHHqJTp04kJydTuXLl4BaSDz//L2sqRDhIUBYBc9KzpsFw9uxZ+vXrxzPPPMN9993H/PnzKVu2bPAKKIDmo6HCJdB8VESKF3FKgnIRE87dQEpXyZoW1h9//MEdd9zBG2+8wZNPPsm0adOCPkuvIPasgOP7YU9yxKog4pCMUy5ivLuB4FHr84bS6cNZ08L4/fff6dixI5s2beKVV15h0KBBhc9UiBgkQbmIuaoz/LABrro90jXJv/3795OYmMjOnTuZO3cu3bt3j3SVgMjMWBRCgnIRs2uR2qp+1+LY2Kp+9+7dtGnThkOHDrF8+XJuvfXWSFfpT97dQIQIJwnKRUws7fX26aefkpSUhMfjYf369Vx//fWRrpIQEScX+oqYWNnrbeXKlbRq1YrzzjuP1NRUCchCuCQoi7B75513aNeuHXXq1CE1NZUrrrgi0lUSImpIUC5i9n4M0xurNBpNnDiRXr160bRpUzZs2ED16tXzfd/ipbOmQhRFEpSLGN8dpqOJ4zgMHz6cIUOG0KVLF1asWEGlSpUKlEdGRtZUiKJIgnIR0/A+tdxkg3sjXZO/nD17lr59+zJu3Djuv/9+5s2bR5kyZQqcT0l3Yl+pyEzwEyIsJCgXMZ9NUstNbpkc6ZooJ0+epHPnzsycOZORI0fy6quvUrx48YDyuvDvbnptECsoRJSRoFzEVK6l0kq1IlsPgN9++402bdqwbNkypkyZwqhRo/B4At/47rfdKg3VqnT+or1/XhRNEpSLmFLl3bRcZOuxd+9emjdvzqeffsq8efPo379/ofMse4GbXljorPIlWvvnRdEmQbmIidQO07527NhBkyZN+PHHH1mxYgXdunULSr4N7wNPCWgYpv7yhIlQozG0mRCe8oQACcpFTr274ZLGUK93ZMr/+OOPadasGWfOnGHDhg20atUqaHmnjQfnLGx+LmhZ5urSG6HP5tiYri6KDgnKRcy2WbD/I9j2VvjL/vDDD2ndujVVqlQhNTWVhg0bBjX/cuertGzVoGYrRFSRoFzEHP7OTfeEt9zZs2fTsWNHrrrqKlJTU6ldu3bQy7jFVIvO32IGPWshooYE5SJmvztSYP8n4SszlRfo3bs3zZs3Z/369Vx00UUhKWfnArXo/M4FIcleiKggQbmIKe7OyQjHVGQHh2QeJZmh3HHHHSxfvpyKFSuGrDzvXnm/WCErQoiIk6BcxKSfcNOTIS4nPZ1F/ItUTG5gIO+++y6lJM6BpgAAIABJREFUS4f2m8C7Z97No0JajBARJUE5xMK5Zx5AiTJZ01A4ceIEt99+O18wi9Y8TTteCXiWXkHsS1PdF/s2h7woISJGFrkPsXDumQdw8hc3/Tk0+R86dIj27dvzySef0JFpNKJfaArKRrgX8D9iq9ev0cDoX59aFB0SlEMs7DuBOH5pENm2TUJCAnv27GH+/Pls79ol+IXkItzbM4X7C1UIkKAcHiEIkOH21VdfkZCQwNGjR0lOTqZFixZsj3SlQiyWttYSRYf0KYdY6ljV2kodF+maBG7z5s00a9aMs2fPkpKSQosWLSJSj3D3z8fK1lqiaJGgLHK1bNkybrnlFs4//3zS0tKoX79+xOri7U74bErEqiBEyEn3RYjVuxsObI3cWhSFMWvWLPr27UuDBg1Yvnw51apVi2h9pDtBxANpKYdYJNeiCJTjODz33HPcc889tGrVinXr1kU8IIN0J4j4IEFZZJGZmcnQoUMZNmwYd955J8uWLeO88yK4DqgQcUa6L0Ks6XAoXSk2fnKfOXOGPn368Pbbb/Pggw8yYcIEihWT720hwkmCcoiFe2xtoI4fP063bt1ITk7m2Wef5bHHHst766ZiQCbye0uIIJKPU4jFwj5vv/76K7fccgurVq3ijTfeYPjw4fnaS6/JY4DHTYUQQSFBOcSWP6CW01xe+C3qQuKHH36gWbNmbNu2jUWLFtG3b9983/f7NYADP6wNXf2EiDcSlEOsWt2saTSxLIsmTZpw8OBBVq5cSceOHQt0/6YjoEJ1aDI8RBX0E+7JI0JEggTlEGs0UO2ZF20X+jZt2sTNN9+M4zhs3LiRm2++ucB57FkBxw/AnuQQVDAbMnlExIO4DMrhbHHtWqTGKe9aHPqy8mvJkiXcdtttVKtWjbS0NOrWDawZf/wnNz0YxMrlotFA1X8dbV9wQgRTXAblcLa4ajRRP/Fr3BT6svJj+vTpdO7cmbp167Jp0yYuu+yygPPy7gASrp1AZPKIiAdxGZTD2eJKGal+4m8cFfqycuM4DmPHjuXee+/l1ltvZe3atVx44YWRrZQQ4hxxGZTD2eLyXuC7UA99WTnJzMxkyJAhjBgxgp49e/LBBx9QoUKFQufrfUyRfGxCFDUyeSTEGg2EQ1+7i+lEwJkzZ7j77ruZO3cuQ4YM4YUXXgjaLL0K7qbVFS4OSnZCCOK0pRxOkVyQ6NixY7Rr1465c+cybtw4XnzxxaBOm653N9RoHJsr4AkRrSQoF1En+IXWrVuzbt06ZsyYwbBhw/I1S68gdi2CfWEcWSLjlEU8kO6LEKuVqIJWrYTwlXmY73mLNpzcvpfFixfTvn37kJQT7vWNZc88EQ8kKIfYGkONvlhjwNUFmzAXkINsYzaJnOUUq1evpkmTJiErK9yLLcki9yIeSFAOsaN73XRf6MtKSUlhBh0pRQX6spEmTa4NfaFhFCsr7glRGNKnHGIZZ9z0dGjLWbx4MW3atKECF3MvaVSjaAVkIeJFXAblcF4wKl7KTUuHrozXX3+drl270qBBA/qyicr8LXSFCSFCKi6DcupYdcEodVzoyypdSaVlKgU/b8dxGDNmDP369SMhIYE1a9ZQnguCX5AQImykTznE0k+o9MyJ4OabkZHB4MGDmTx5MnfddRczZsygZMmSwS1ECBF2cRmUw7lvXoWL4bejwZ31dvr0aXr37s28efN45JFHeO6552QvPSGKiLj8JIdz7YvL27jpbcHJ7+jRoyQlJTFv3jxM0+T555+PWECWyRxCBF9cBuVwBpMtk7KmhfHTTz/RqlUrNmzYwKxZsxg6dGjhMy2EcC86L18CIh7EZfdFLM4M27NnD23atGH//v0sWbKEpKSkSFdJZvQJEQJxGZRjbWbYF198QWJiIunp6axdu5bGjRtHukqAzOgTIhTisvsirDtYePzSAlq3bh3NmzenVKlSbNq0Kc+AfFEDN20YWHnRTHYeEfEgLoNyWDl+aQEsWLCAxMRENE0jNTWVa665Js/71E4CPFC7bcHLE0JEngTlKDV16lS6d+9Oo0aN2LhxI5qWv+Zh2njAcVMhRMyJy6AczVfxHcdh1KhR9O/fn6SkJFatWkXVqlXzff/y1VTq3RVECBFb4jIoh3OadUFkZGQwYMAARo8ezT333MOiRYsoV65cgfKodatKL78lBBUUQoRcXAblM8ezptHg1KlTdO/enalTpzJs2LCAp01f3Q0qVFepECL2xGVQ/vl/WdNIO3LkCG3btuX999/nxRdfZNy4cQFv3ZQyUi2qv3FUcOsohAiPuByn/NMXbro1svUAOHjwIImJiWzfvp05c+bQq1evQuVXuZZ6fJUuD1IFhRBhFZct5UtudNMIz8H45ptvaNq0Kd988w1Lly4tdECGvy7wBXMBJCFE+MRlS7npCFj+gFotLlI+//xz2rZtS0ZGBmvXruWGG24ISr7hXAFPCBF8cdlSThkJJyLY77pmzRpatGhBmTJlSE1NDVpABpn1JkSsi8ugXLmWSiPR7zpv3jySkpK47LLLSEtL46qrrgp/JYQQUSsug/LZk276R3jLnTx5MnfeeSc33HADKSkp1KhRI7wVEEJEvbgMyns/ypqGmoPDGp5k0KBBdOjQgZUrV1KlSpXwFC6EiClxGZRLlFVpybKhLyuDsyzhfjYwhr59+7Jw4ULKlg1DwUKImBSXoy8yTqk01N0Xp06dYh492cEimvM4r7/+n4AnhQgh4kNctpQr1cyahsLvv/9OQkICO1lMEhO5lTFhCcjRvNiSECJvcdlS/nVn1jTY9u/fT2JiIjt37qQb71CXO0NTUDZkyyQhYltcBuWLGsL+zXBxCHbn2L17NwkJCfz6668sX76czbfdGvxCciFbJgkR2+IyKJ854qZHg5vvp59+SlJSEh6Ph3Xr1tGoUSM2B7eIPIV73zwhRHDFZZ/y+VertGoQ522sWrWKVq1aUaFCBVJTU2nUqFHwMhdCxI24DMpewbru9u6779KuXTtq165NWloaV1xxRXAyFkLEnbgMyvs/Vem+Twqf18SJE+nZsydNmjQhJSWF6tWrFz5TIUTcisugfOxg1jQQjuMwYsQIhgwZQpcuXVixYgWVKlU65zxPyaypEELkJi4v9BXzQKabBuLs2bPcf//9zJgxg379+jFlyhSKFy+e7blOetZUCCFyE5ct5cwzWdOCOHnyJF26dGHGjBk89dRTTJ06NceALIQQBRWXLeVAHT58mA4dOpCWlsbkyZMZMEAGAwshgkuCcj7t27ePhIQEvv76a+bNm0e3bvncLroYbl9JKGsnhCgq4jIol6qkJpCUOve6XLZ27txJQkIChw8f5sMPP6R169b5LyzTLxVCiFzEZVA+cyJrmptPPvmEpKQkSpQowYYNG2jYMARzs4UQwhWfP6rz2XpdsWIFrVq1olKlSqSmpkpAFkJky9Ttf5u6/e9g5CVBOQdvv/02HTp04MorryQtLY3atWuHpWpCiNhi6vb9wFhgrKnb/QqbX3wG5Ty89NJL3HXXXdx8881s2LCBiy66KPDMPH6pEKLIMHW7FvA88CAwGHjB1O1CbcksQdmH4zgMGzaMhx9+mG7durF8+XIqVqxYqDy9W0+VKBeECgohooap28WAN4EFhqW9bljaNGAR8KZ7W0Di8kJfdtLT07nvvvuYNWsW/fv355VXXgnKpJDSldTu2WUKF9uFEFHGsLRMoLnfsd6FzTcuW8reoXDe9OTJk3Tu3JlZs2YxevRoJk+eHLRZeqfdNZtPBXntZiFE0RSXLeX0P/5Kf/vtN9q3b8/HH3/M1KlTuf/++4NaVsbprKkQQuQmLoNyseKQARz12DRrlsC3337L/Pnz6dKlS/DLKgEZZ1UqhBB5ictQkfEH/MIOZp1uA/uOkpycTMuWLUNSVmZm1lQIEftM3R4ATAb6GpY2I5h5x2Wf8o+ezbxBMzI5S0pKSsgCMkB1d77JJdeFrAghRPh5P9Fbgp1x3AXl1NRUZjm3UJaq3Esq9evXD2l50qcsRJF0PXAK2B7sjOMuKNetWxedO7mXVKpSK+TlJU2FGo2h7ashL0oIEQambpcGrgW2GZZ2Ntj5x12fcsWKFelcfIa60heGtenPuwRqtlCpEKJIqAuUBLaYut0AeBJoCZQBUoGHDUuzAs087lrKgArIvmkIfTYZ0sbDZ1NCX5YQIiyud9OaqCDsANOBz4HbgLWmblcONPO4aymHW40mUKE61Lgp0jURQgSJNyjfCDQzLG2r9wZTt98C/gkMAJ4NJPP4bCmH0RoDjh+ANY9GuiZCiCDxjrx41Dcguya7ad1AM5egHGLH9mdNhRCxy9TtkqiA+wPwVjanHHTTMj73ed7U7RX5LSMug3Kx0m5aJvfzgqF4yaypECKm1QVKAR/kMPKippv+4HPs/4BP8ltA1PQpezyeO4BrgF9QD+wTx3HeD0VZme6Y4cxTocg9qw5vwof9ZUicEEWEt+vi+xxu7+ymK91W9QnUSI3mpm4/CXxlWNq1uRUQFUHZ4/G0B4YCTRzHyfB4PMWBNI/H84fjOB8GvcDihG1I3NUd1T8hRJHgvchX1f8GU7erA/cDXwMrUXsb3QR8hroo+COQ5zSyaOm+eAJ4z3GcDAA3fQ8YHpLSZIdpIURgvEG5h6nb5b0HTd2ugOpjLg08ZFjaWXe95erAMeBTw9IOGpZ2OK8CIt5S9ng8lYFGwES/m34FbvJ4PBUdxwnuasSOXxpCO5fA8gfUzD5pMQsRu0zdLoHqU94KlAO+MHV7MSoQdwEuAYYYluZ7Ua8h8KVhafmONtHQUq6F6kg44Xf8GOpLo1D7XUXa8gfgxAHVryyEiGnXokZVfAbcAvwP1V1xD2oNjNaGpb3sd58GqCCebxFvKQPuLnbZBmVQ30hZeDyefoB319jTHo+nQFMaL+F6708QnvJsCfoqT77KUqVSRf5W8+j+H3942HP4SCjLEkIUmuU4jp7dDYalfUnWLZDzswB7faBA18WiISh7m/X+g8a8f5/Tmncc5zXgNQCPx/OZ4ziNQlc9IUS88Hg8nwU5yxLA1aZuXwKcNCzt97zuEA3dF94W8nl+x71/nwxjXYQQIpgeB+4E9gJj83OHaGgpfw+kc25QrogaH/GD/x2EECIWGJb2NvB2Qe4T8Zay4zhHgA3ARX43XQxsdRzntzyyeC0kFRNCxKOIxxOP44RhXFhelVCTR/4D/J/jOGc9Hk8J1BXO0Y7jLIps7YQQInyiIigDeDye7kBjYCdqvN8ux3HejWythBAivKImKAshhIiOC30BCecCRkKIos/j8VQFRgOrHcf5b6TqEZNBOewLGAkhiiyPx9MQ6IkatjYIdT0rYiI++iJA4V3ASAhRZDmOs9VxHMNxHP/1dyIi5oKyzwJGB/xu+nMBo/DXSgghgiPmgjJFfAEjIUR8i8WgXOAFjIQQIlbEYlAu8AJGQggRK2IxgMkCRkKIIisWg/L3yAJGQogiKuaCchAWMBJCiKgVc0HZNRHo7i5chJt2A56JaK2EEDHL4/F4r0uVimg9YnXtC1nASAgRDB6PpybwAGpT1HbALiAZWOc4zuKw1ydWg7IQQgSDx+MphmodnwUyUD0IJQDHcZwzYa+PBGURi0zdvhG4ELgAWGVY2r4IV0mIoIjJBYlE9DJ1ux7wPGoFv0tR48rTUKv5gXrPXQj8DJiGpW0MsKh/ALcBOtAS+DMom7r9BGq3806GpRVoe3efPMoAT6K6xk4A5YEHDEs7HWB9hciXWL3QJ6KUYWnbDEtrAzRwD31nWFozw9I6u/86GJbWGPgSSDF1e2SA5TwMPJLDzc0BDbg6kLxdM4BrDEv7F+oL5Z+oLxkhQkqCsggJw9IOuf/NaTLPU8D/gJGmbjcLsJicWq13oFrP7wWSqanb5VGjed4GMCxtNFDVsLRvA8lPiIKQ7gsREYalOaZur+GvK96bgpi3dyx7oC5ATdv/1SfPo4WtlxD5IS1lEUmV3PRYrmcJEUekpSwiwtTtqkAn4Hdgdja3dwKSULtBVEYNVXrCsLRchyiZuv0kkIC6mPikYWnz/G5vBvQB9gAV3LyHua1rTN1+Dajtnj7R1G0b2GtYWv8AH6oQBSJBWYSdqdsa8C5wGOhpWJrtd/tIVL/w/xmW9od7bC6wDDXiIjcvA0dQsz7L+t5g6nYf4D9AA8PSfnGPjQPSTN1uaFjaGcPS+pm6XRO1xspLhqXNKtSDFaKAJCiLUKvmN8KiBpAIfA70MyztK9+TTd2+HnURsIc3ILtGA1+Zut3GsLSVORVmWNoRU7eXooKyb76ae+w/3oDsegY1iuMe4LWCPjghgk2Csgi1n93RC1mYut0GSDV1+7/AIMPSjrs3PYC61vE/3/MNS9th6vZRoBmQY1B2ZWRzrDequ8I/32Ombu9w85WgLCJOgrKICMPSVpq6/RSqu6EsajIIQH3UEqxdTd1O97vbRuBggEXWd9NW7gQXX7uA7wLMV4igkqAsImk+Kih3N3V7sGFpB1FrZWNY2phszh9fiLK8Af4V/z5sIaKJDIkTkXTW5//e4XFbgGLuxbZzmLod6Ht2i5vWCnK+QgSVvBFFJHVy028NS9vl/v9V4BTQ1f9kU7frAA8HWNYc1GSQ7PKtgqzFLaKEBGUREqZun+/+t0wOt98KPAccRU1pBtQFPWAwMMzU7YY+51cCRgCv+2Tjzbu0X/be45V98v0Z+BfQy73I6M23NDAGmOZzf+9QugtyfoRChIYs3SmCyr2IZqJWidNQq8SloMYkg7qOURM1QiIFeMYNmP75/B8wDPgDNYEE4AXD0n51b18MXOeWYQOfG5Z2u6nbE4C2wJWolnGaYWmdfPK9AhgFeNz7FQOmGJb2nXv7NKA1UAc103Ar8LJhaQsL+9wIkR8SlIXIganbHlTwdgxLkw+KCAsJykIIEUWkT1kIIaKIBGUhhIgiEpSFECKKSFAWQogoIkFZCCGiiKx9EeVM3W6Cmm12AfCHYWk3RLhK+eaOJb4UVfdehqWlRrhKYWfq9s3A0wT59XNnIdZ388WwtAXByFdEngTl6PcRaqeM91C7aWRh6vZy4BLgJr/1h6NBH6A/asZcgd9rpm7rwJOoBedrAUtjcNH5VNTzMBeoFsR8NdQa0F2Ab4GIBmVTt/+Fep3vMyxteSTrEuuk+yLKGZaW6c422+F/m6nbZYCmqM1Hz/e/Pb9M3Q7JbDXD0n4D0gK5r6nblwDrgLcNSxsGNEIFoZiS2+tXyHy3GZZ2D7A0mPkWwk2oxkH9vE4UuZOgHDsy/Q8YlnYKaABcZ1ja3nPvkjdTt0uQw8ppQXJOvfOpPVDasLQl7t/XkvdWUNEsVLO0Toco34IajJqe/kKkKxLrpPsixnnXbCiE9vgs3BNFqgOHvH8YlnYygnUReXC7ztZFuh5FgQTlOOV2fdwKvIFaeEeIAjF1uxSQLuuCBJcEZT+mbqegLqKcj/q53B7VanNQF9qm+G7caeq2CXRAXQV/FLUaWhPgZuB+w9K+dM8rBgwFLkLttqwBHxuW9oZf+R5gANAQ2A+cQF3oyq6uC9x8LgRaZLMrdEvUnne/8tfOG6+5y2N2BnqiVmG7yNTtD3zueqdhaSd88mmGuli1B7XPXWVgmGFpR/zKuw54EPgN+BkoDnyRXd1zYur231E7jFwJXOZTr1nA34C+7uN9FViLen1uAkYZlrbaJ59bUBfBfkO9z88Dnjcs7XufczagRodciOqXb+qWW9bNc5xhactN3U4AGqO6+24A5gZ6wdF9fP9ErRldGbXM6Ej/lfLy+37xOb8S8BBq9b2/od6/A7PJtzlwvfsYK7r5jjAs7Qefc/oBg9znJB34wLC0ru5tbYF5qPfBZ0AX91gv1PM43bC0LF0Ypm5fiFp29Q+3fpWBFMPS5vuc8zD5fG2LOlmQyI+p28VRAeg11EWqfxqWtse9rTGQjFpu8jn3WDFUAP0MtX39d6iNPfcCYwxLe9LNcwXwo2Fpfd37lQe+BGYalvaMe6wksBy1pGRfbwvE1O3OwBTglGFpl/vUtTpqmcxewGV+H6yHgMeB5t4F5N0P5FvANd6RGqZuzwRa+ubr93z0cR9XA+8u0KZuj0N9ETU0LO2Me+w+1JfSrd56uPV7B2jplrGhAK/DSNSH0eN3vCawHdXCt1E7VB8FFhqW9k/3nMdQy3d2MCztqHvsGuB94CHD0la5x8qgRoe86Nb9Y8PSUtzb7gcmAO2AKt6lO03dTkRdXKvjG+Dz8XhmotaNfsl9XJnu8U6o17aNYWnb3WP5er/45NsQtbXWJO8Xpanbm4HthqXd63Pu31E7sDxsWNqr7rHuwHSgvvd97h4vCWwDThuW1sDvsQwALjcszfCpWwKw0H1so33OrY8K4v28r7+b99vAIcPS+vucm+drGw/kQp8fw9IygK/dP8f6vlENS/sI9QF6xrv5pnt13bvVUGdglmFp+4BWqA8gqBbMrcBwn7xOoC6KDHfHnAI8gWqtPeT7k9CwtEXAzmzqegDVosjC1O2rUcF6rM+OHqACUE1U6yhPpm5rqA/GRG9Adj2DalHe4553DWqvvbG+Xwxu/d7MT1n55eb/K2qj1dcNSzuLeq4fc+vS2K3fQ96A7N5vB+rD/qbbqvReKN3qnnKxNyC7tqJasYP91lL+HPULoH0A1S+Feo7+vPhpWNp/UcH2dZ8tqfL7fvGqC6z2++WyFbjd77yyqBb6KZ9856G6rx71PdGwtHRgKlDfDea+rsRnpxa3bov9H6x7EXkWsNz3C9nNexDQz9Ttrj7Hc31t44UE5ex5A2J2fa2zUT+HH8jmts+9HzjD0ta7Q8JAdUf8nM1i7h8D5YGGbrfFENSH63g2ef+eQ10zsjn2ECoArMjmeAvD0nbnkJe/3qifqf/zPWhY2jHUEK9m7qEBqAB2zgeTnOtdGA6w2xt0DUv72P0iBLUw/m/ebiM/yahhW7398oK/9vDz8r5WH/se9HkNywVQ75M5jCX/EPUz/Xr37zzfL37Hz6B+qfn62b+OhqVtMSytimFpM/3O/RbI7pfSbNTojr7eA6ZulwOKGZaW5XX1/aLx0Q41RO6crgf3sW1FdWv4yu21jQvSp1xw36KGefl/MED1AWfh/rSrDfxg6vYwv5vLod6wp4ErUH18B4NQxwaoN3eWkRlua/eXbO+RPe+Y01beXwY+dvnkfz3qgs9vhM85z7VPXQ7lcJv3eKN85OfkcNyrbA7HA/GrmzY0dfsr8vd+8fWT+wvPl5NTHU3dvgpog3p9DwE1gB/8zzMs7TdTt98Hepu6PdztquqG6o7ID++XTG6vRytTt0sblub7mHJ6zuOCBOWC8+5GkV3L4FQ2xzLdf3sMSxufU6ZulwMEZzxrCTcfT14nZlOPcsCF7k9J78XBV/wvIvopRujG4eYku+caVIs9p5mN3ucju2CVU/3D8bi8v1gd8vl+8ZOvOroNhBmoL6WhqD5ox9Tt3KZ+vw70ADqiZg02NSztrXzWy7tXYk7vQw+qK6gkWb9ocnpt44J0XxTclag3k//P3Wy5P1e/IpcJGm5f4m7gJIWYmedjK+q1rZ1DebkF62qoC3Pw12PMtu4+faBbgVKmblcocE2DbwtQzaduvi5xU/+f+pFW3U23FOD9EohpqH7mdoalLfK5bvHn+8HUbf/3zHrgG6Cve9vX5J/3/VM9h9svAXbl0F0XtyQo5+6SbI79C/VNPqUA+byEGt51nf8Npm73QrU+MoHJwC3u1Wl/fyP/Ld+pwFlUn55/ec2Au3wO7QV8LxxV5q9+4Dmon9Zd8eNebPJe7PGW1zaHelOAuhfWJNRjaJzNbQmoxzYnTHXxV97U7YrZHO+ECn7ei455vl8CLL8TakSG/0XjOj7/zzLKwQ3c01HdHSNQ/cz5tQzVxZXof4Op2zVQszRfKUB+cUGCcu56mbp9kfcPd+xrX2CQ7xvb3aYecmjluhdWZqOusFfxud9VQBPD0ja6h0ahWiJjfe/vDvG6GDjfHRHhy/sT8c9ZeYalbUVdsR7hDkny5lMZ9aXyns/91wOV3MV/QLWSP3Lz+dk9v5ep2238Hu8YVMsL96La48Aod0yq97yGwN3un3oOXzY5Kevmkd1zWoacn+tlqJEnU0zdvsCnLv+HupB6r2Fpvn2WZfzSXI/7vNYFeSxen6NeE9+W6b9QY6X7eFuuBXi/eOvnX/c/6+2OgPDajAr2fz537pC8FP5aLCm7C5hvorpUKhuW9lN2D8znefF9Hx5HNQA6mrrd3u/cd1B901OzqXcwfi3GLBmnnA1Tt1ugglUHVIvLg7rqfTEw3g163nNHo4bC1UX1ZW4F5njHgfrlexfqQon3gtvvwIvu0B/vOSVRQ+MuRk18qIAaXzsAaI5quU5HtdQXoy44VkNN7JhtWNoon7xaoQLRAdSklgzgOd+hYu55g1EfnhTgf4alvel3+xWoLwwPavxoMdQkmu/8zktAjZk+iBr9sRd1YXSmW2/bsLRW/s+LXx5Xuo/tBtSEjx8ACzVi4h+oYXg3oFrmW4FlvuNiffJJch/T725dSqHGl+/yOee/7vOnuY9rm2Fp7U3dnoYainUF6nXaiupXHcxfr/XvqC+vhwxLy/MnvTue+BXUc9gbNSnofNQIn9HuiBb/++T4fnF/8TyNmqhUGvgUFciT3efv/1AXjj8HlhiWNtr9Un4G0IFN7nOY5j6+VaiLt68alnbOdGlTt99BjZFelc1tj7nPTz338XwK3O59TKZaXGoU6n1zxn3cy3z7pk3d7k8+X9uiToJyNnyCckEnPHgvAuY0REgUkNt/6vEfXSDPddEV76+tBOVsBBqUhRCisKRPOXul/VIhhAgLaSn7MdUWRg1RowZ+BLYaluY/XVUIIUJCgrIQQkQR6b4QQogoIkFZCCGiiARlIYSIIhKUhRAiikhQFkIvzD8fAAAAD0lEQVSIKCJBWQghosj/A6Bfct5egCvAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x720 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_encoding, ax = plt.subplots(2,1,figsize=(5,10))\n",
    "plt.suptitle('Encoding', fontsize = 24, fontweight='bold')\n",
    "fig_encoding.tight_layout(rect=[0, 0.1, 0.9, 0.9])\n",
    "\n",
    "ax[0].scatter(bb_df_test['behavior_belief1'], nb_df_test['behavior_belief1'],\n",
    "              s=10, c=[(130/255,0,1)], marker='.')\n",
    "ax[0].plot(belief_array, belief_array, c='black')\n",
    "ax[0].set(xlim = [0,1], ylim = [0,1],xticks=[], yticks = [0, 1])\n",
    "ax[0].tick_params(labelsize=14)\n",
    "#ax[0].set_title('box1', fontsize = 20.0)\n",
    "ax[0].spines['top'].set_visible(False)\n",
    "ax[0].spines['right'].set_visible(False)\n",
    "ax[0].text(0, 1.02, '$b^{\\check}_t$',  fontsize = 20.0, color = neural_color)\n",
    "ax[0].text(.1, .8, 'box1',  fontsize = 20.0,\n",
    "        bbox={'facecolor': 'none', 'alpha': 0.5, 'pad': 8})\n",
    "for label in (ax[0].get_xticklabels() + ax[0].get_yticklabels()):\n",
    "    #label.set_fontname('Arial')\n",
    "    label.set_fontsize(18)\n",
    "    \n",
    "\n",
    "ax[1].scatter(bb_df_test['behavior_belief2'], nb_df_test['behavior_belief2'], \n",
    "              s=5, c=[(130/255,0,1)], marker='.')\n",
    "ax[1].plot(belief_array, belief_array, c='black')\n",
    "ax[1].set(xlim = [0,1], ylim = [0,1],xticks=[0, 1], yticks = [0, 1])\n",
    "ax[1].tick_params(labelsize=14)\n",
    "#ax[1].set_title('box2', fontsize = 20.0)\n",
    "ax[1].spines['top'].set_visible(False)\n",
    "ax[1].spines['right'].set_visible(False)\n",
    "ax[1].text(1.02, 0, '$\\hat{b}_t$',  fontsize = 20.0, color = neural_color)\n",
    "ax[1].text(.1, .8, 'box2',  fontsize = 20.0,\n",
    "        bbox={'facecolor': 'none', 'alpha': 0.5, 'pad': 8})\n",
    "for label in (ax[1].get_xticklabels() + ax[1].get_yticklabels()):\n",
    "    #label.set_fontname('Arial')\n",
    "    label.set_fontsize(18)\n",
    "\n",
    "fig=ax[0].figure\n",
    "fig.text(0.45, 0.08, 'Belief \\n predicted from behavior', color=neural_color,\n",
    "         va='center', ha='center', fontsize = 22)\n",
    "fig.text(0, 0.5, 'Belief predicted from neurons', color=neural_color,\n",
    "         va='center', ha='center', rotation='vertical', fontsize = 22)\n",
    "fig_encoding.savefig('encoding.pdf', format = 'pdf', \n",
    "                     bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### build neural belief data using the learned model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#neural estimate of belief using current linear regression model\n",
    "nb= regr.predict(r_df)  # input is ALL data\n",
    "nb_df = DataFrame(nb, columns = bb_df.columns)\n",
    "#save neural estimate belief\n",
    "nb_df.to_csv(path_or_buf='./data/nb_df.csv',index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # plot: compare raw data vs prediction\n",
    "# S = 5 # dot size for scatter plot\n",
    "# plt.scatter(bb_df['behavior_belief1'], nb_df['behavior_belief1'], s=S, c='blue', marker='.')\n",
    "# plt.xlabel('Beliefs from behavior model (true)')\n",
    "# plt.ylabel('Beliefs from neural code (estimated)')\n",
    "# plt.plot(belief_array, belief_array, c='b')\n",
    "# plt.xlim([0,1])\n",
    "# plt.ylim([0,1])\n",
    "# plt.title('Encoding: Box 1 (ALL data)')\n",
    "# plt.show()\n",
    "\n",
    "# plt.scatter(bb_df['behavior_belief2'], nb_df['behavior_belief2'], s=S, c='red', marker='.')\n",
    "# plt.xlabel('Beliefs from behavior model (true)')\n",
    "# plt.ylabel('Beliefs from neural code (estimated)')\n",
    "# plt.plot(belief_array, belief_array, c='r')\n",
    "# plt.xlim([0,1])\n",
    "# plt.ylim([0,1])\n",
    "# plt.title('Encoding: Box 2 (ALL data)')\n",
    "# plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig_encoding, ax = plt.subplots(2,1,figsize=(7,15))\n",
    "# plt.suptitle('Encoding', fontsize = 24)\n",
    "# fig_encoding.tight_layout(rect=[0, 0.1, 0.9, 0.9])\n",
    "\n",
    "# ax[0].scatter(bb_df['behavior_belief1'], nb_df['behavior_belief1'], s=S, c='blue', marker='.')\n",
    "# ax[0].plot(belief_array, belief_array, c='b')\n",
    "# ax[0].set(xlim = [0,1], ylim = [0,1],xticks=[0, 1], yticks = [0, 1])\n",
    "# ax[0].tick_params(labelsize=14)\n",
    "# ax[0].set_title('box1', fontsize = 20.0)\n",
    "# ax[0].spines['top'].set_visible(False)\n",
    "# ax[0].spines['right'].set_visible(False)\n",
    "\n",
    "# ax[1].scatter(bb_df['behavior_belief2'], nb_df['behavior_belief2'], s=S, c='red', marker='.')\n",
    "# ax[1].plot(belief_array, belief_array, c='r')\n",
    "# ax[1].set(xlim = [0,1], ylim = [0,1],xticks=[0, 1], yticks = [0, 1])\n",
    "# ax[1].tick_params(labelsize=14)\n",
    "# ax[1].set_title('box2', fontsize = 20.0)\n",
    "\n",
    "# fig=ax[0].figure\n",
    "# fig.text(0.45, 0.08, 'Beliefs from behavior model $\\hat{b}_t$', va='center', ha='center', fontsize = 18)\n",
    "# fig.text(0, 0.5, 'Beliefs from neural encoding $\\check{b}_t$', va='center', ha='center', rotation='vertical', fontsize = 18)\n",
    "# #fig_encoding.savefig('figures/encoding.pdf', format = 'pdf', bbox_inches=\"tight\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig_encoding, ax = plt.subplots(2,1,figsize=(5,10))\n",
    "# plt.suptitle('Encoding', fontsize = 24, fontweight='bold')\n",
    "# fig_encoding.tight_layout(rect=[0, 0.1, 0.9, 0.9])\n",
    "\n",
    "# ax[0].scatter(bb_df['behavior_belief1'], nb_df['behavior_belief1'],\n",
    "#               s=10, c=neural_color, marker='.')\n",
    "# ax[0].plot(belief_array, belief_array, c='black')\n",
    "# ax[0].set(xlim = [0,1], ylim = [0,1],xticks=[], yticks = [0, 1])\n",
    "# ax[0].tick_params(labelsize=14)\n",
    "# #ax[0].set_title('box1', fontsize = 20.0)\n",
    "# ax[0].spines['top'].set_visible(False)\n",
    "# ax[0].spines['right'].set_visible(False)\n",
    "# ax[0].text(0, 1.02, '$b^{\\check}_t$',  fontsize = 20.0, color = neural_color)\n",
    "# ax[0].text(.1, .8, 'box1',  fontsize = 20.0,\n",
    "#         bbox={'facecolor': 'none', 'alpha': 0.5, 'pad': 8})\n",
    "# for label in (ax[0].get_xticklabels() + ax[0].get_yticklabels()):\n",
    "#     #label.set_fontname('Arial')\n",
    "#     label.set_fontsize(18)\n",
    "    \n",
    "\n",
    "# ax[1].scatter(bb_df['behavior_belief2'], nb_df['behavior_belief2'], \n",
    "#               s=5, c=neural_color, marker='.')\n",
    "# ax[1].plot(belief_array, belief_array, c='black')\n",
    "# ax[1].set(xlim = [0,1], ylim = [0,1],xticks=[0, 1], yticks = [0, 1])\n",
    "# ax[1].tick_params(labelsize=14)\n",
    "# #ax[1].set_title('box2', fontsize = 20.0)\n",
    "# ax[1].spines['top'].set_visible(False)\n",
    "# ax[1].spines['right'].set_visible(False)\n",
    "# ax[1].text(1.02, 0, '$\\hat{b}_t$',  fontsize = 20.0, color = neural_color)\n",
    "# ax[1].text(.1, .8, 'box2',  fontsize = 20.0,\n",
    "#         bbox={'facecolor': 'none', 'alpha': 0.5, 'pad': 8})\n",
    "# for label in (ax[1].get_xticklabels() + ax[1].get_yticklabels()):\n",
    "#     #label.set_fontname('Arial')\n",
    "#     label.set_fontsize(18)\n",
    "\n",
    "# fig=ax[0].figure\n",
    "# fig.text(0.45, 0.08, 'Belief \\n predicted from behavior', color=neural_color,\n",
    "#          va='center', ha='center', fontsize = 22)\n",
    "# fig.text(0, 0.5, 'Belief predicted from neurons', color=neural_color,\n",
    "#          va='center', ha='center', rotation='vertical', fontsize = 22)\n",
    "# fig_encoding.savefig('encoding.pdf', format = 'pdf', \n",
    "#                      bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig_encoding_jitter, ax = plt.subplots(1,2,figsize=(15,7))\n",
    "# plt.suptitle('encoding of neural activity', fontsize = 24)\n",
    "# fig_encoding_jitter.tight_layout(rect=[0, 0.1, 0.9, 0.9])\n",
    "\n",
    "# ax[0].scatter(bb_df['behavior_belief1'] + 0.06 * (2 * np.random.random(bb_df['behavior_belief1'].shape)-1), nb_df['behavior_belief1'], s=S, c='blue', marker='o')\n",
    "# ax[0].plot(belief_array, belief_array, c='b')\n",
    "# ax[0].set(xlim = [0,1], ylim = [0,1],xticks=[0, 1], yticks = [0, 1])\n",
    "# ax[0].tick_params(labelsize=14)\n",
    "# ax[0].set_title('box1', fontsize = 20.0)\n",
    "\n",
    "# ax[1].scatter(bb_df['behavior_belief2']+ 0.06 * (2 * np.random.random(bb_df['behavior_belief1'].shape)-1), nb_df['behavior_belief2'], s=S, c='red', marker='x')\n",
    "# ax[1].plot(belief_array, belief_array, c='r')\n",
    "# ax[1].set(xlim = [0,1], ylim = [0,1],xticks=[0, 1], yticks = [0, 1])\n",
    "# ax[1].tick_params(labelsize=14)\n",
    "# ax[1].set_title('box2', fontsize = 20.0)\n",
    "\n",
    "# fig=ax[0].figure\n",
    "# fig.text(0.45, 0.08, 'Beliefs from behavior model $\\hat{b}_t$', va='center', ha='center', fontsize = 18)\n",
    "# fig.text(0, 0.5, 'Beliefs from neural encoding $\\check{b}_t$', va='center', ha='center', rotation='vertical', fontsize = 18)\n",
    "# fig_encoding_jitter.savefig('figures/encoding_jitter.pdf', format = 'pdf', bbox_inches=\"tight\")\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}