Laminar information flow - solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "16da5cb6",
   "metadata": {},
   "source": [
    "# MaDeLaNE: Workshop on effective connectivity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dd9e21c",
   "metadata": {},
   "source": [
    "Connectivity can be characterized in two ways, functional connectivity and effective connectivity. While functional connectivity describes statistical dependencies between data, effective connectivity deals with the causal influences between the data. \n",
    "\n",
    "<b>Mutual information</b>[1] and <b>directed information</b>[2] are two metrics that can be used to study the connectivity between two neurons. \n",
    "\n",
    "<b>Mutual information</b> quantifies the amount of information that one random variable (or neural signal, in our case) contains about another. It is a symmetric measure, and better suited to study the functional connectivity between neural signals. \n",
    "    \n",
    "<b>Directed information</b>, on the other hand, is a causal metric. It measures the amount of information that one time series carries about another when predicting its future values, thus indicating the direction of influence. Directed information provides an estimate of the effective connectivity between the signals. For two random variables, directed information is described as, \n",
    "    $$I(X^{n} \\rightarrow Y^{n}) = H(Y^{n}) - \\displaystyle \\sum_{i=1}^{n} \\left[H(Y_{i}|Y^{i-1}, X^{i})\\right] $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20ee6e24",
   "metadata": {},
   "source": [
    "Lets first import some basic dependencies and simulate single unit activity to understand how DI is estimated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fcb4ed1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib as mpl\n",
    "from random import seed, random\n",
    "import pandas as pd\n",
    "from scipy import signal\n",
    "from scipy.stats import f_oneway\n",
    "\n",
    "import utils\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e915dee",
   "metadata": {},
   "source": [
    "We will use the <b>Izhikevich neuron model</b> [3] to simulate spiking activity in two neurons that are connected such that neuron 1 causally influences neuron 2. Both of these neurons are provided with input current simulating the presense of a stimulus during the entire trial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4d2c23c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initializations\n",
    "duration = 1 # in seconds\n",
    "num_trials = 10  \n",
    "num_units = 2\n",
    "step_size = 1e-3 # in seconds\n",
    "Fs = 1/step_size # sampling frequency\n",
    "T =int(duration / step_size) # duration in samples\n",
    "\n",
    "# weights - 2D matrix summarizing the strength of connectivity between the units\n",
    "#         - rows corresponds to effect variables, columns are the causal variables\n",
    "weights = np.array([[0, 0], [0.7, 0]]) # neuron 1 -> neuron 2, with a strength of 0.7\n",
    "\n",
    "input_current = [100, 100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "690c0a50",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>spike_times</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[[59, 138, 205, 278, 319, 397, 452, 528, 596, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[[48, 80, 151, 216, 283, 300, 354, 415, 470, 5...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                         spike_times\n",
       "0  [[59, 138, 205, 278, 319, 397, 452, 528, 596, ...\n",
       "1  [[48, 80, 151, 216, 283, 300, 354, 415, 470, 5..."
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simulate spiking\n",
    "\n",
    "simulated_spikes = utils.simulate(num_units=num_units, num_trials=num_trials,\n",
    "                                  duration=duration, step_size= step_size, weights = weights, inp_curr=input_current)\n",
    "\n",
    "simulated_spikes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8053fc14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAEUCAYAAAAhnVB2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZF0lEQVR4nO3dfbRsdX3f8fcnXBVBUTBHIw9GbKxd1GiEiZXYmCloi2BDlnVZk/qU2tB2tY0aU4PVVaPLtI1xWapJTYgPodFiE0MTYhW1VmK6kkUcEQW8Kvh85RKOWnwuYPj2j9nU4/Gee8/DzN4z83u/1jprZvaZ8/t9f3tmvvdz99lnJlWFJEmS1JrvG7oASZIkaQgGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQVjqWZKnJvmzJN9McuXQ9UjSMkryqiQ3JPlako8leebQNWn57Bu6AGnRJdlXVd+e4ZBfBi4C/gZw1gzHlaSFNYde+g3g7wOfAH4UuCLJjVX1ZzOcQyvOI8JaWkk+k+QXk3wkyVeS/LckR2/4/pOSXJPk1u4I7CM2fK+S/NCG27+T5BXd9XGSA0l+KcnNwJuS3CPJRUlu6r4uSnKPTfd/QZJbkhxM8rNb1V1V/7Oqfg+4aR77RZJ2Yol76Uur6mNVdWdVXQX8KXDmHHaRVphBWMvuqcA5wKnAI4BnAyQ5HXgj8E+B+wG/BVx+V8Pdhh8ATgB+ELgAeDHwGOBHgEcCjwZesun+9wFOAp4D/EaS43e/LEnq1VL30iT3ZHpU+Ppt1iUBBmEtv9dU1U1V9WXgj5k2V4CfA36rqq6qqr+qqkuA25g24O24E3hpVd1WVd8C/hHw8qq6parWgZcBz9hw/zu6799RVe8Avg48bM+rk6R+LHsv/U3gw8C7tlmXBBiEtfxu3nD9m8C9uus/CLyg+1XerUluBU4BTtzmuOtV9X833D4R+OyG25/dNNaXNp37trEWSVp0S9tLk/wa8HDgqVVV26xLAgzCWl2fB36lqu674euYqrq0+/43gWM23P8HNv385mZ6E9N/EO7yIDzHV9LqW+hemuRlwBOBv1tVX93tOGqXQVir6reBf5bkb2Xq2CTnJbl39/1rgJ9JclSSc4CfOMJ4lwIvSbKW5PuBfwu8eTeFdXMezfRdW74vydFJ7rabsSRpzha5l74I+BngCVX1pd2MIRmEtZKqasL03LZfB/4PcCPdH390nsv0bXduZXrO2h8eYchXABPgI8C1wNXdtt14BvAt4HXAj3fXf3uXY0nS3Cx4L/13TI8o35Dk693Xv9nlWGpUPJ1GkiRJLfKIsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUn7hi5gO84555y64oorhi5DkvqWWQ5mL5XUqC176VIcEf7iF784dAmStPTspZL03ZYiCEuSJEmzZhCWJElSkwzCkiRJapJBWJIkSU2aWxBO8sYktyS5bsO2E5K8J8kN3eXx85pfkiRJOpx5HhH+HeCcTdsuBN5bVQ8F3tvdliRJkno3tyBcVe8Hvrxp8/nAJd31S4Cfmtf8kiRJ0uH0fY7wA6rqIEB3ef+t7pjkgiSTJJP19fXeCpSkVWIvlaStLewfy1XVxVU1qqrR2tra0OVI0lKyl0rS1voOwn+Z5IEA3eUtPc8vSZIkAf0H4cuBZ3XXnwX8Uc/zS5IkScB83z7tUuDPgYclOZDkOcB/AJ6Q5AbgCd1tSZIkqXf75jVwVf30Ft86e15zSpIkSdu1sH8sJ0mSJM2TQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCY1FYTH4zHj8Xjhx+xr/HnXvtt5+qpr3vPudry91DGrNSzKvhhq3FnOMdTzedaGeh3Pcv9tNdYy/tuw13kWvb55jbtK/XXR9u1exh6yTzYVhCVJkqS7GIQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJalKqaugajmg0GtVkMhm6DEnqW2Y5mL1UUqO27KUeEZYkSVKTDMKSJElqkkFYkiRJTWoiCI/HY8bj8dzuP89aln3+Wc3XZ93zmmse4w69f2e5plmMtUyP3TLbvD/2un8W7Xm0l7GXtVf1UXcfr/FZr2PIvrQoPXGR/+2ahSaCsCRJkrSZQViSJElNMghLkiSpSYME4STPT3J9kuuSXJrk6CHqkCRJUrt6D8JJTgJ+HhhV1cOBo4Cn9V2HJEmS2jbUqRH7gHsm2QccA9w0UB2SJElqVO9BuKq+ALwK+BxwEPhKVb178/2SXJBkkmSyvr7ed5mStBLspZK0tSFOjTgeOB84FTgRODbJ0zffr6ourqpRVY3W1tb6LlOSVoK9VJK2NsSpEY8HPl1V61V1B3AZ8GMD1CFJkqSGDRGEPwc8JskxSQKcDewfoA5JkiQ1bIhzhK8C3gZcDVzb1XBx33VIkiSpbfuGmLSqXgq8dIi5JUmSJPCT5SRJktSoVNXQNRzRaDSqyWQydBmS1LfMcjB7qaRGbdlLPSIsSZKkJhmEJUmS1CSDsCRJkprUZBAej8eMx+OlnGOW4+51rFnVsmjj9D32LOfpq855zrfIz4e+9+8y2GqfLFufXcSxFq0fLPNrarfz9Pma7+u1NIvxFvH1sltNBmFJkiTJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSk1JVQ9dwRKPRqCaTydBlSFLfMsvB7KWSGrVlL/WIsCRJkppkEJYkSVKTmgvC4/GY8Xi8kD/Tx1jzGn/eNfYxz6zG3us4s15jX4/NLOdbxNdP3/tx0W21P3a7n4bus/MYb1bj9t335jHfIr0OF7G/bHfsRdiPy/K6267mgrAkSZIEBmFJkiQ1yiAsSZKkJg0ShJPcN8nbknwsyf4kZw5RhyRJktq1b6B5/xNwRVU9JcndgWMGqkOSJEmN6j0IJzkOeBzwbICquh24ve86JEmS1LYhTo14CLAOvCnJh5K8PsmxA9QhSZKkhg0RhPcBpwOvq6pHAd8ALtx8pyQXJJkkmayvr/ddoyStBHupJG1tiCB8ADhQVVd1t9/GNBh/l6q6uKpGVTVaW1vrtUBJWhX2UknaWu9BuKpuBj6f5GHdprOBj/ZdhyRJkto21LtG/CvgLd07RnwK+NmB6pAkSVKjBgnCVXUNMBpibkmSJAn8ZDlJkiQ1yiAsSZKkJqWqhq7hiEajUU0mk6HLkKS+ZZaD2UslNWrLXuoRYUmSJDXJICxJkqQmrWwQHo/HjMfjhR9zGcedx3yL9Hj1td/mPc8i7dNZjzGvcWddW9+vwT7sZk07+ZlFed5u52eWqVcs037diz4ek1nNsdU4i/Rv6bzXughWNghLkiRJh2MQliRJUpMMwpIkSWqSQViSJElN2lEQTnJ8kkfMqxhJkiSpL0cMwkmuTHJckhOADwNvSvLq+ZcmSZIkzc92jgjfp6q+CjwZeFNVnQE8fr5lSZIkSfO1nSC8L8kDgacCb59zPZIkSVIvthOEXw68C7ixqj6Q5CHADfMtS5IkSZqvfUe6Q1X9PvD7G25/CvgH8yxKkiRJmrctg3CS1wK11fer6ufnUpEkSZLUg8MdEZ70VoUkSZLUsy2DcFVd0mchkiRJUp9SteXZD9M7JGvALwGnAUfftb2qzppvad8xGo1qMvEAtaTmZJaD2UslNWrLXrqdd414C7AfOBV4GfAZ4AMzKUuSJEkayHaC8P2q6g3AHVX1J1X1j4HHzLkuSZIkaa6O+PZpwB3d5cEk5wE3ASfPryRJkiRp/rYThF+R5D7AC4DXAscBz59rVTMwHo8BuPLKKxd2nnnUOOsxF3W8eT2+837e7HT8ZV1nX3Mt2uPVCvvr3sdcxLFmMc6Rxli0nrZI/8Ytwhpa7Knb+UCNuz5W+SvA35lvOZIkSVI/DveBGi+sqldu9cEafqCGJEmSltnhjgjv7y59rx1JkiStnMN9oMYfJzkKeHhV/etZT9yNPQG+UFVPmvX4kiRJ0uFs+fZpSfZV1V8BZ8xp7ufynaPOkiRJUq8O9z7Cf9FdfijJ5UmekeTJd33tZdIkJwPnAa/fyziSJEnSbm3n7dNOAL4EnMX0j+bSXV62h3kvAl4I3HurOyS5ALgA4EEPetAeppKkdtlLJWlrhzsifP8kvwBcB1zbXV7fXV632wmTPAm4pao+eLj7VdXFVTWqqtHa2tpup5OkptlLJWlrhzsifBRwL6ZHgDf7nrdT24HHAj+Z5FzgaOC4JG+uqqfvYUxJkiRpRw4XhA9W1ctnPWFVvQh4EUCSMfCLhmBJkiT17XCnRhzqSLAkSZK0Eg53RPjseU9eVVcCV857HkmSJGmzLY8IV9WX+yxEkiRJ6lOq9vJ3b/0YjUY1mfhJz5KaM9NT1Oylkhq1ZS893DnCkiRJ0soyCEuSJKlJBmFJkiQ1aeWD8Hg8Zjwe9/6zfY8/j1pnNeYy7cdZjbcs+27IOVf9eTHEYzcPO1nH0Ptoq/EW8fW417H66lNDv44WeYydjt9HT+jrObqIz4vdWvkgLEmSJB2KQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJqWqhq7hiEajUU0mk6HLkKS+ZZaD2UslNWrLXuoRYUmSJDXJICxJkqQmGYQlSZLUpGaC8Hg8Zjwer8xc85hjL2POa817HbfFx2LeYy/ba2nZ6l0Gh1vnvPfBIr82jjTWor32Fvn1tZNxd1rDoj6HlunfqyEe93nN2UwQliRJkjYyCEuSJKlJvQfhJKckeV+S/UmuT/LcvmuQJEmS9g0w57eBF1TV1UnuDXwwyXuq6qMD1CJJkqRG9X5EuKoOVtXV3fWvAfuBk/quQ5IkSW0b9BzhJA8GHgVcNWQdkiRJas9gQTjJvYA/AJ5XVV89xPcvSDJJMllfX++/QElaAfZSSdraIEE4yd2YhuC3VNVlh7pPVV1cVaOqGq2trfVboCStCHupJG1tiHeNCPAGYH9Vvbrv+SVJkiQY5ojwY4FnAGcluab7OneAOiRJktSw3t8+rar+N5C+55UkSZI28pPlJEmS1CSDsCRJkpqUqhq6hiMajUY1mUyGLkOS+jbT08jspZIatWUv9YiwJEmSmmQQliRJUpOaC8Lj8ZjxeLxU8/RR87zm6Gt/z2u+3Y63aI//Mj6H9jrevNfc93N7SDtZ6zL12L7HWLTXyG7HXJR9P6uxFqmWeY83qzFXJTNAg0FYkiRJAoOwJEmSGmUQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJqaqhazii0WhUk8lk6DIkqW+Z5WD2UkmN2rKXekRYkiRJTTIIS5IkqUnNBeHxeMx4PF6q8edV827Gnef+m8XYi15fH2POe95F3g+zrm2ox2eRbXefLPLzpM/x9jrHLGtcxB47717U52t4lfvYoq5tFpoLwpIkSRIYhCVJktQog7AkSZKaNEgQTnJOko8nuTHJhUPUIEmSpLb1HoSTHAX8BvBE4DTgp5Oc1ncdkiRJatsQR4QfDdxYVZ+qqtuBtwLnD1CHJEmSGjZEED4J+PyG2we6bd8lyQVJJkkm6+vrvRUnSavEXipJWxsiCB/qY+6+53Oeq+riqhpV1Whtba2HsiRp9dhLJWlrQwThA8ApG26fDNw0QB2SJElq2BBB+APAQ5OcmuTuwNOAyweoQ5IkSQ3b1/eEVfXtJP8SeBdwFPDGqrq+7zokSZLUtt6DMEBVvQN4xxBzS5IkSeAny0mSJKlRqfqeN2xYOKPRqCaTydBlSFLfDvUuO7tmL5XUqC17qUeEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJalKqaugajijJOvBZ4PuBLw5cTt9ccxtccxt2uuYvVtU5s5p8Qy/dTS2roLU1t7ZeaG/Nra0XdrfmLXvpUgThuySZVNVo6Dr65Jrb4JrbsEhrXqRa+tLamltbL7S35tbWC7Nfs6dGSJIkqUkGYUmSJDVp2YLwxUMXMADX3AbX3IZFWvMi1dKX1tbc2nqhvTW3tl6Y8ZqX6hxhSZIkaVaW7YiwJEmSNBNLE4STnJPk40luTHLh0PXMQpJTkrwvyf4k1yd5brf9hCTvSXJDd3n8hp95UbcPPp7k7w1X/d4kOSrJh5K8vbu90mtOct8kb0vyse7xPrOBNT+/e15fl+TSJEev4pqTvDHJLUmu27Btx+tMckaSa7vvvSZJ5liz/ZTlfc5tZC9d7V7aQh8dvIdW1cJ/AUcBnwQeAtwd+DBw2tB1zWBdDwRO767fG/gEcBrwSuDCbvuFwK9210/r1n4P4NRunxw19Dp2ufZfAP4r8Pbu9kqvGbgE+Cfd9bsD913lNQMnAZ8G7tnd/j3g2au4ZuBxwOnAdRu27XidwF8AZwIB3gk8cU712k+X/Dm3ad320hVdcyt9dOgeuixHhB8N3FhVn6qq24G3AucPXNOeVdXBqrq6u/41YD/TJ/75TF/sdJc/1V0/H3hrVd1WVZ8GbmS6b5ZKkpOB84DXb9i8smtOchzTF/obAKrq9qq6lRVec2cfcM8k+4BjgJtYwTVX1fuBL2/avKN1JnkgcFxV/XlNO/p/2fAzs2Y/nVra59xd7KVN9NKV76ND99BlCcInAZ/fcPtAt21lJHkw8CjgKuABVXUQps0duH93t1XZDxcBLwTu3LBtldf8EGAdeFP3K8zXJzmWFV5zVX0BeBXwOeAg8JWqejcrvOZNdrrOk7rrm7fPw6rt6+/RUD+9CHvpyvbSxvtobz10WYLwoc7zWJm3u0hyL+APgOdV1VcPd9dDbFuq/ZDkScAtVfXB7f7IIbYt1ZqZ/o/+dOB1VfUo4BtMf9WzlaVfc3c+1/lMf3V1InBskqcf7kcOsW2p1rxNW62zz/Wv9L5upZ/aS1e/l9pHD2nmPXRZgvAB4JQNt09m+uuBpZfkbkyb9luq6rJu8192h/npLm/ptq/Cfngs8JNJPsP0V7JnJXkzq73mA8CBqrqqu/02ps18ldf8eODTVbVeVXcAlwE/xmqveaOdrvNAd33z9nlYtX39/zXWT+2lq99LW+6jvfXQZQnCHwAemuTUJHcHngZcPnBNe9b9ReMbgP1V9eoN37oceFZ3/VnAH23Y/rQk90hyKvBQpieHL42qelFVnVxVD2b6OP6vqno6q73mm4HPJ3lYt+ls4KOs8JqZ/irvMUmO6Z7nZzM9Z3OV17zRjtbZ/erva0ke0+2vZ274mVmzn35n+9I+5+ylwOr30pb7aH89dFZ/9TfvL+Bcpn8F/EngxUPXM6M1/W2mh+4/AlzTfZ0L3A94L3BDd3nChp95cbcPPs6c/qq8x/WP+c5fOq/0moEfASbdY/2HwPENrPllwMeA64DfZfpXviu3ZuBSpufv3cH0qMRzdrNOYNTtq08Cv073gUdzqtl+usTPuUOs3V66omtuoY8O3UP9ZDlJkiQ1aVlOjZAkSZJmyiAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYSyvJ/ZJc033dnOQL3fWvJ/nPc5rzeUmeOYNx3prkobOoSZL2wl6qlvk+wloJSX4Z+HpVvWqOc+wDrgZOr6pv73GsnwCeXlU/N5PiJGkG7KVqjUeEtXKSjJO8vbv+y0kuSfLuJJ9J8uQkr0xybZIrktytu98ZSf4kyQeTvOuuzzjf5Czg6rsad5Irk/zHJO9Psj/Jjya5LMkNSV7R3efYJP8jyYeTXJfkH3Zj/Snw+O4fBElaOPZStcAgrBb8NeA84HzgzcD7quqHgW8B53UN/LXAU6rqDOCNwK8cYpzHAh/ctO32qnoc8JtMP9f8XwAPB56d5H7AOcBNVfXIqno4cAVAVd0J3Ag8cqYrlaT5sZdq5RiE1YJ3VtUdwLXAUXQNtLv9YOBhTBvue5JcA7wEOPkQ4zwQWN+07fINY11fVQer6jbgU8Ap3fbHJ/nVJD9eVV/Z8LO3ACfucW2S1Bd7qVaOv0pQC26D6ZGDJHfUd06Mv5PpayBMG++ZRxjnW8DRhxq7G+u2DdvvBPZV1SeSnAGcC/z7JO+uqpd39zm6G1OSloG9VCvHI8ISfBxYS3ImQJK7Jfmbh7jffuCHdjJwkhOBb1bVm4FXAadv+PZfB67fXcmStHDspVo6HhFW86rq9iRPAV6T5D5MXxcX8b2N9Z3A7+5w+B8Gfi3JncAdwD8HSPIA4FtVdXAvtUvSorCXahn59mnSDiT578ALq+qGPY7zfOCrVfWG2VQmScvDXqpF4akR0s5cyPQPPfbqVuCSGYwjScvIXqqF4BFhSZIkNckjwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCb9P1P/zhffrZSyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make the raster plot\n",
    "fig, ax = plt.subplots(1, num_units, sharey = True, figsize =(10, 4))\n",
    "spike_times = simulated_spikes.spike_times.values\n",
    "for unit in range(num_units):\n",
    "    ax[unit].eventplot(spike_times[unit], linelength = 0.8, color = 'k')\n",
    "    \n",
    "    # plot clean up \n",
    "    ax[unit].set_title('neuron ' + str(unit+1))\n",
    "    ax[unit].spines['top'].set_visible(False)\n",
    "    ax[unit].spines['right'].set_visible(False)\n",
    "    ax[unit].set_xlabel('Time (ms)')\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "ax[0].set_ylabel('Trials')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c09746b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to evaluate spike rates\n",
    "def calculate_rate(T, win, spikes = None, spike_times = None):\n",
    "    if spikes is None:\n",
    "        spikes = np.zeros(T)\n",
    "        spikes[spike_times] = 1\n",
    "    S = np.zeros(T)\n",
    "    for t in range(0, T):\n",
    "        S[t] = np.sum(spikes[t - win:t + 1])\n",
    "    return S/win\n",
    "\n",
    "# window for boxcar filter\n",
    "rate_w = int(0.01*Fs)\n",
    "\n",
    "# evaluate spike rate for each trial\n",
    "rates = np.array([[calculate_rate(T, win = rate_w, spike_times = spike_times) \n",
    "                   for spike_times in row.spike_times] for _, row in simulated_spikes.iterrows()])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81243fb1",
   "metadata": {},
   "source": [
    "How does the correlation changes if we were to change the strength of connection between these units?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bbf0b3b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to evaluate the avg. cross correlation between rates\n",
    "def avg_cross_correlation(rates1, rates2, maxlags): \n",
    "    num_trials, T = rates1.shape\n",
    "    cross_correlations = [np.correlate(rates1[i], rates2[i], mode='full') for i in range(num_trials)]\n",
    "    \n",
    "    if maxlags is not None:\n",
    "        lags = np.arange(-T + 1, T)\n",
    "        trim_idx = (lags >= -maxlags) & (lags <= maxlags)\n",
    "        cross_correlations = [cross_corr[trim_idx] for cross_corr in cross_correlations]\n",
    "        lags = lags[trim_idx]\n",
    "    else:\n",
    "        lags = np.arange(-T + 1, T)\n",
    "\n",
    "    return lags, np.mean(cross_correlations, axis=0), np.std(cross_correlations, axis=0)/np.sqrt(num_trials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f47b7573",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAE4CAYAAAD2NJx0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABc7UlEQVR4nO3dd3ib5dX48e/Rsrwdj0zi7BASCCEJGwIBwihtoVB4Wza0L3TTvi20vy6g5aW0lK63LS0FyiiUQguFloYRIISwMskkZG873kOWZK3798cjGcXxkmMt+3yuS5etR8+4Jcs6uu/nfs4RYwxKKaVUOtnS3QCllFJKg5FSSqm002CklFIq7TQYKaWUSjsNRkoppdJOg5FSSqm002CkVAqJyMMismgA9jNeRIyInDYQ7VIq3TQYZRERGSkifhGpFhFnutujUkNEtorI7Z0W7wFGAe+lvkVKDTwNRtnlBuAFoB64KM1t6TMRcaW7DQOlu+eS6udojAkbY6qNMcFUHlepZNFglCVExAb8N/BI9HZj3GP/LSLNIpLbaZtvi8i+6LaIyDkisi7au1orImdEh3quSrAtBSLyKxHZIyLtIrJTRL4bfSw2fHSliPxHRNqAu8TyLRHZLiIBEdkmIl/vtN+LRGS1iHhFpElElonIcdHHnCLyCxHZGz1mlYg82Ye2zhGRF0WkRUQ80X2eGPf4tSKyMbrPvSJyp4g44h5fLCIPisiPRaQK2Nfdc4yu/xkReT/6Gu+Mtjm/h/bNFpGFIlITbd9yETk//vjAJOC26DFN9PiHDNOJyJEi8kJ0Px4R+ZeITI57/DoRCYnIqSKyKvo6LxeROb29jkolnTFGb1lwA84HagAH1vBMAJgYfawY8AGf7bTNeuCn0d/HAF7gAWA6cDawEjDAVQm0Q4DFwHbgYmAiMA/47+jj46P73AtcFX18AvDlaBtvBKYAXwD8wOei242MPqdbo+sfBVwBHBN9/H+i+zwTqASOB77eS1tnAG3AX4G50eN+Fjg5+viFQBj4f8BU4L+ARuDHcftYDLQCf4i+bsf08Byvi25/ddzrshZ4LG5/DwOL4u6fCVwb3fdU4M7o6zA1+ngpsAP4efQ1GgnY49pwWnS9XGAX8CowJ3p7HdgKuKLrXAdEgCXA6cA04GVgG+BI93tcb0P7lvYG6K2Pfyh4Fvhl3P3/AHfF3X8SWBh3f3b0w2pG9P7/AjsBe9w655N4MDo7us3cbh6PfUj+oNPyPcDPOi37JbA9+vtx0e3Gd7PfXwOvAZJAWx8D1gC2bh5/E3iq07KbsYJm7AN8MbA5fh89PMedwBc6LZsXXXdY9P5Bwaibdq0Bvhd3fytwezevcywYfQ7ry0Z53Dojos/lmuj966LbzI5b56TosiPT/R7X29C+6TBdFhCRUcDHsYbnYh4Gro8bUnoUWCAiI6P3rwZWGmM2RO9PB5YbY8Jx+3inH82ZAzQaY1b0st6yuPYXAUdgfSOP9wYwXkTysHoQLwHrReRZEblZRMbGrftnrF7JVhH5g4hcGjtPIyKVcUNTHhH5Q1xbXzXGRLpp44xu2uTGGhqLWdnNPuKfYwUwDvhFfFuAhdFVJnexPSJSISK/F5FN0aFJT7Rd47ppc3dmABuNMXWxBcaYA8CH0cc6FmMFu5h90Z8jEjyeUgPK0fsqKgN8DutvtUJE4pfbgU8Cz2B9kNcCV4rIr7GGo+7qtJ/OKdr7m7K9L9u19WG7jidjjAmLyAVYw2/nAJcCd4vIZcaYfxtj3heRCcACYD5WT+nHInISsB+YFbfflgTa2l2b4pd39Vw6L499sbsZa3iss73d7ONhrGHHW7GG43xYvdz+TIjo6rlKp+WRTl9IYo/pF1OVVvoGzHDRyQefxwosszrd/kJ0IkP0A+YJ4BrgXKxzDX+N29VG4HgRscctO7kfTVoJlIrI3L5uYIxpwfowPqPTQ/OAHcYYb3Q9Y4xZZoy5yxgzD6uXcn3cfjzGmGeNMV/DOgd0FHCGMSZkjNkad6uJa+s5sQkcXdjQTZt8WOfE+izaC9mDNdy1tYubv5tN5wG/N8Y8b4xZB1RhnW+KF8D64tGTDcAMESmPLRCREVjnoTZ0u5VSGUKDUeY7H+ub8x+NMevjb1hDVwtEZHx03UeAmVjnhxYaY2rj9vN7rKGY+0TkKBGZH10P4r45i8irIvKTHtrzGta5lr9FZ79NiM7O+nwvz+MnwFfFmvk3RURuAr7IR7PQThGRH4jIidFht7Ojz2Vj9PFborPXZkR7SDdgTT7Y3MMxf4Y1aeFxEZkrIpNE5DIRiQXhnwCXish3RGSqiFwO3A7ca4wJ9PJ8uvI94Gsi8n0ROTo6u+1iEfljD9t8iNWbPUZEZmF9gegceHYAp0Zfl/JugusTWD3jv0Vn6M3B6mHtA/7Wj+eiVGql+6SV3nq+Ac8B73TzmB2oBu6MW7YaK7hc2sX652DNsGvHOkdzfud1sU7CP9xLmwqB/8P6Fh/A+rD8TvSx8cSdWI/bRoBbousGsXoeX497fAbWpIzqaPt2Affw0USCm7B6Oi2AB1gOXNSH1+8EYBHWkFor1kWiJ8Q9fi3wQfR57MMK0I64xxcDD3TaZ5fPMfrYxVjn4rzRtr4P/DDu8Yc5eDbdMcDbWL2xncCXou19OG6dudHn7osed3xXbQCOjL6Gnujt38DkuMevA0Kd2ntEdD9npvu9rrehfRNjtNLrUCUisaGwmcYaIlJKqbTQYDSEiMgXsWZS7ceaXfdLrJlxJ6W1YUqpIU9n0w0t47Au8ByBNRz2CvDttLZIKaXQnpFSSqkMoLPplFJKpZ0GozjRpJgPpLsdauiRfiSsVWow0WCUBNJ1/RmVAUTEISJ3iZUdvFVE6kTkJYnL5K1SS0TcIvLn6N8kICJb090mlXoajFTGkeTWBsrByjxxL1aS0DOxLhZdJCKTethuyBJLMos52rGu87of60JdNQRpMDqUTUTujn5jbhGRB+TQOkFfjSa29IvIFhH5XixhqXRff2apiNwZt487oo+dE7fsDRH5Wdz9BSLyloj4xKpL9GcRKevUlh7r58SGHqPZDapFpEGs0tfd1tiJbmdE5Esi8li0B7FHRG7ttI5DRG4XkR3R42+IZlbovJ+rOi1bJCIPx93fKVYdod+LSD3wVnT5x0RkpVi1hmqij8c/t4ej+7pRRHZF/17PiZW0tEvGmDZjzHxjzF+MMRuMlcniWqxsDh/r5TX5XxH5QKw6QHvESthaHPd4n+oFich8sepJxepKze/puNFtbo/2uC+KvvfaROT1zgFUrPpNL4uVqLVWRJ4RkXGd99Npm9Ni79NOz2O+iKzGugj5PBEpFJE/RvfrF5EVInJu3H5iNZYuF6uWkles+lVX9/Tcon+Tm4wx95FgGiY1eGgwOtSngTKsei9XYiUi/WnsQbGG376FNUX6KKzEmDcBt0VXuQTrSvp7seoOjcLKWfYaVvmFmLOwvpGfHd1vLtY39dei98/Cyr7wJFZanIuxrrp/VsTKlioi1wH3RY81HSsv3TlYtXc6P6dSrF7AFdF93UrvbsPKaj0LKxvCTzt9cD4Qfb43RV+LH0XX+Vwf9t3Z17DqNZ0MXCsiM4Hn445/LVbm8s7P7XisxKkXYmWUmIVV+ycRuYATqOtlvVg9pulY2QzOBH7TaR0bVpqhm7HKeDQCT8lHX1ZGY2VGWBl9/JtYSV/7YhRWCqUrgVOAEuCh2IMiMh3rIuZ3sLI2nIUVZF8REXcfjxH/PH4Wbd80rMwVDwHnYdVwOg7rS8O/RWRap23vxirfMRN4CviziExJ8PhqqEl3CohMumGlftnJwTV/bsT6ZpgP5GGleTm/03bXAE1x97uqP3MmEAKKovtpx/pHXxZ9fAHWUEV+XFvu7rSPSqzULbOi93fSe/2cxcDaTuv8gW5SDMWtY4DfdFq2CfhJ9PcJWIXapnVa54fA+532c1WndTqnu9mJVeohfp3HYq9N3LKLosccF73/MFZAz4lb5ztAVYJ/9weibShIcLtPRf+Otuj96+ilXhBW8bxdHJxy6ONdvU6djnV79P1TEbfsM9HXwx33ejzZabuc6Hv24rj9bO20zmnE1ZKKex6nx60zObrsY522XQU8FP19fHSd/4l73IGVmuimPr6mh7RPb0Pjphe9HmqZOTjF/ltY6fwnYf1j5wL/EJH4C7TsgFtEKszByUnjvYMVbOZh5WbbhVWD6O7oUM9Z0WPHyhIcD5wkIl/pYl9TRGQfH9XPie8JxEogTMbK3wZWfrR4+7Aye/emq+1idW/mRo/VuayFA+vbeKKWdbo/g2gvMc4b0WNOx3r9AD4wxrR308ZeicjdWD3Fs4wxnl7WvQT4OtZrW4TVe3BhVV/dH12tp3pBH0bbvswYE4pbZ2kfm7u/0/trH9brMRzYjfWemSxWTaR4bqyEsYlaHvf79OjPzvWflnBo9vf3Y78YY0IicgCtl6R6ocGod/GftLFhzcvoOlt0Q3c7Mca0i8jbWMNyAeA1Y0ytiGzC6jWdhVWTKP5YP8XqIXRWjdVTg77Vz+mcgdrQtyHanraL/TwF65t35/Xif5dOj3d1Mrwv9Y+6Wt5VGzsf7xDRoc5Y3aezjTFre1n/ROBprCG4W7CG307CypQeP+Git3pBnesL0cX97nT1XOP3bcN6v9zdxbb1sfbRt79H2HRf9iJeV8+nv+83NYRpMDrU8SJij/tAORnrn2sb1j+eH5hojPlPD/vorv7Ma8B/RR//WdyyT2FVJY1PzbMCq2R4d9NcPSISq5/zp96f1oBbGf1ZaYz5dw/r1QCjY3dEJAfrW/aOXvbfVa2hM7A+2DYm1tSDiVXT6QGsc0xnmo+q4fbkNKDOGPP9uP18uh+H3wBc3ek9dlo/9tOVFVjnabYZY7oLcDXA8E7Hn92Hfcdeo3lYmcFjTsfKFK/UYdFvK4cqA34nVs2fC4EfA38y1owfD1b9nbtE5Cti1auZIdaMtp/G7aO7+jOvYZUMmMVHvZnXsE5IBzm4DPgPgYtE5JciMkusWjzni8iD8tHsvv7UzxkQ0SD5EPAnEblaRCaLyLEicoOIxAfVRcAXRORkETka67xGX6Zu3wPMFmt24DQROR+rbMXjxpjd/W13dCLBk1gTUy4H6kVkZPRW0MOmHwIVIvI5EZkoItdglXtI1H1ABXB/9D12Nh/VlTpcd2FNJPmLiJwgVq2p+SLyaxGJFex7Heuc5Y+j76nLgC/3tmNjzDasnuHvReS86N/k18DRWH+rwyIi08Wq5zQScEXf87MkudP8VQbRYHSov2PVvVmK9aH1H+Jmnhljfgx8A6v66proet/AOgEecxtQjPUBVos18QCsMfhWYKMxJjZzK3Ye5K34cx/GmNexhu6OwSpmtxYry3YrVuDCGPMY1gfqhVjnXJZjnQCOnadIthujbfoeVm/lVaxZb/HTc7+FVUPpJWAh1jmG5fQiOmz2Saze0Bqs4acXgC8cZpuP4KPZhUuwajLFbt/qoT3/xgoadwHrsCYP3JLowY0x+4BPYNVZeh9rqPB/Et1PN/v+AGvYtADr9d4I/AnrPGdTdJ0Pgf/Gav96rCKF3+3jIT4f3e9fsP4mpwIfN8ZsGoDm/werh3UTMDb6+2rietVqcNNEqUoppdJOe0ZKKaXSToORUkqptNNgpJRSKu00GCmllEo7DUZKKaXSbtBc9Hr++eebF198Md3NUEplpl6zcqj0GjQ9o7q63hIuK6WUylSDJhgppZTKXhqMlFJKpZ0GI6WUUmmnwUgppVTaDZrZdEopdThWrlw53OFwPICViVy/qA8sIyLN4XD4z5FI5L45c+Z0rnmlwUgppQAcDscDI0eOPKqioqLRZrNpBukBZIwhEAg49+/f/9WWlpbZWNn9D6LRXymlLEdXVFS0aCAaeCJCTk5OcNy4cc10U0xSg5FSSllsGoiSK/r6dlUFW4ORUgpuueBRbrng0XQ3Qw1hGoyUUmoAVFe1HFld1XJkutuRrTQYKaW6pL2loeu5554rnDBhwozc3NzjTjzxxKmbN292dbfugQMH7AsWLJiUm5t73OjRo4/5wx/+UNqfY2owUkop1aGqqspx1VVXTfrBD36wv76+/v1Zs2Z5L7vssondrf/5z3++0uVymerq6jV//vOfd9xyyy2VK1ascCd6XA1GSg0x2uPJLr/+9a/LzjrrrMmx+5WVlUd/7GMf6wgOI0eOnPn222/nDtTxHn/88ZLJkyf7b7jhhsa8vDxzzz337P/www/zVq9efUiAaWlpsb344ovDfvKTn+wrLi6OnHfeeZ6zzz67+aGHHipL9LgajJRStPuCtPuC6W6G6sKCBQtaly9fXhAOh9m1a5czFArJihUrCgA2btzo8nq9thNPPNHX1baFhYWzurt997vfHdnVNhs2bMidPn26N3a/qKgoMnbs2PY1a9YcEozWrVuXY7fbmTlzZnts2cyZM72bNm1KODjqRa9KDWFtzX4e/vFiPlxZBcDtn3mKy24+mRknj01vwzJAfV3b2GAwnNfX9YPBcC5YExn6uo3TafeWlefv6Wmd6dOnB/Lz8yPvvPNO3saNG3POOOOMlvXr1+etXr3avWTJkvy5c+d67PYuZ0vT2tr6fl/bEtPW1mYrLy8PxS8rLCwMt7S0HHKQ1tZWe0FBQTh+WXFxcdjj8XTdoB5oMFJqCDLG8OY/P+D3t7xE4wEPZaMKsDtsrH97D++8sJmjTjgCb4uforIBG/1Rh+HEE09sXbRoUeHWrVtz5s2b11pcXBxetGhRwbvvvltw+umntw7ksfLz8yOdA4/H47EVFRWFO69bWFgYbmtrO2iEraWl5ZAA1RcajJQaYgLtIfZtbeDOq//B5GNHcvuTl/PA9xcB8OO/f4aXHlvDP377Lgd2NTN8bFGaW5s+vfVYOov1iEaOKvpwoNsyb9681hdeeKFk7969rjvuuKNq2LBh4ccff7x01apVBTfffHNNd9vl5eUd191jX/va16ruvvvu6s7LZ8yY4XviiSfKY/dbWlpse/bsyTn22GP9ndc95phj2kOhkKxbty7nmGOOaQdYu3Zt7rRp07ocNuyJnjNSaggJ+EPs3FBLW3M7N/1kAb9ZfANHzhnd8bg738VFXzieP7//ZUoq8qjd28K+bQ1pbLEC67zRe++9V+j3+22TJk0Knnvuua1Lliwpbm5udpxyyine7rbzer2ru7t1FYgArrjiiqYtW7a4H3744RKv1yvf/va3R02dOtV33HHHHRKMioqKIuedd17Td7/73dEtLS22l19+OX/RokUlN9xwQ32iz1GDkVKDWOeZcw/e9hp+b5DKaeVc8pUTsTu6/giwO2yMmjAMsQn3/79Xut2fSo2ZM2e25+XlhU844QQPQGlpaWTs2LHts2fP9jgcAzvANXr06NBjjz227Y477hhTWlp63MqVKwueeuqp7bHHv/Od74ycN2/elNj9Bx98cJfP57ONGDHi2GuvvXbiPffcs3vu3LmHBK7e6DCdUkPEyle38c/fL6NsVAFFpb2fC3K67IwYW8y7C7ew/OWtHH/u5F63UclTW1u7Nv7++vXrP0jWsS6++OLWiy++eENXj0V7VB29qhEjRoQXLVq07XCPqT0jpYaA5jovP7/pX1ROK2fU+JI+b1c2upAjJpdy360vEwwkfE5aqT7TYKTUIGeM4VdffYGWBi/fefBibPZD/+3vWXgN9yy85pDlNpvwhZ+dx75tDfzz98tS0dysNXJU0YfJmLwwVKQsGInIV0RkhYi0i8jDvaz7DRGpFpFmEXlIRHJS1EylBp2GA228/e8Pue62+Uya2eV1jj06fsEkTrpgCo//9E3tHamkSeU5o/3AncB5QLcD1iJyHvAd4KzoNs8Cd0SXKaUS0O4Lsn97I7POGM+lXz0JoMseUFfi17vp7gXcePwfqdrRSOWR5T1spVT/pKxnZIx5xhjzT6C3KX/XAg8aYzYYYxqBHwPXJbl5Sg1K1buaEYFv/fGT2GzS7/2MnljKpV87iaZaL15PYABbqJQlE88ZzQDWxN1fA4wQkUMS74nIjdGhvxW1tbUpa6BS2WDnxhqa67yUjy6kYszhX7x6+ddPRmxC4wHPALROqYNlYjAqAJrj7sd+L+y8ojHmfmPMXGPM3IqKipQ0Tqls8Ze738RmF8pHH/Kv0y/5xW6Ky3JpqvUSaA/1voFSCcjEYOQB4r/GxX4f0PxLSg1GsYtSd2yo4c1nP6B8dCEOZ8I5K7s1bHg+4VCE917cMmD7VAoyMxhtAI6Nu38scMAYk3B6CaWGqsfvfpO8QteA9YpiCkrcOFx2Fj2xbkD3Oxg89vDKIx97eKWWHe+nVE7tdoiIG7ADdhFxi0hXs/keBT4nItNFZBjwfeDhVLVTqWznawvw5j8/4OIvnjCgvSIAEWFYRR7LX95KU23bgO5bZY5Eyo7HrFu3LicnJ2f2RRddNKE/x0xlz+j7gA9rivZV0d+/LyKVIuIRkUoAY8yLwM+A14Fd0dttKWynUlmtZnczeUU5fOrLJ3Z7MevhiA3Vvf70Bs1VNwglWnY85gtf+ELl0Ucf3e9vKKmc2n27MUY63W43xuw2xhQYY3bHrfsLY8wIY0yRMeZ6Y0x7T/tWSll8bQGa631c/MXj+5R/rj/c+S6mHDeKRU+s7X1lddgyuex4zP333z+suLg4fMYZZ/T73H4mnjNSSvVTze5mbHbhU186MSn7j/W0zvnsMWxdU42vTa85SrZMLjsO0NDQYLvrrrvG/OY3v0mo/lNnmrVbqUFi5we1NNf7GD62KGm9opgzPz2D+7+7iMaaNnIn9Ho6ISstfGHT2Lratj6XHa+v9+aCNZGhr9uUV+R7L7hwWtaWHQf45je/OebKK6+smzx5cjDRY8XTYKTUIPGv+1cgwoDPoOtKSUU+J5w3meUvb00oC7jqn0wtO/7222/nvvnmm0Xr16/feLjH1WCk1CDQ1tLOq0+uo6Qif8Bn0HVnwRUzeeeFzbQ2JlxHLSv01mPpLNYjuvq6OUOm7PiiRYsK9+3b5zriiCNmAni9XlskEpHp06e7N27cmFC9JQ1GSg0Cr/1tHT5PgDGThqXsmMefNxm7w0ZjjU7xTrYFCxa0/vCHPxxbVlYWnDRpUnDYsGGtN91004RwOCy9lR1P9FhXXHFF0x133HHEww8/XHL55Zc391R2/Bvf+Ebddddd11GX/s477xy5e/du10MPPbS787q90QkMSmU5Ywz/+tNKphw3ityC1J2/ceU4KKnIo6XBR1uLTnhNpkwtO15YWBiprKwMxW4FBQWRnJwcM3r06ITzRWnPSKksFbu+55rvncGuD2r5n999nEV/Te1065KKfOqrPCx/eStnfnpGSo891GRq2fF4v/jFL/b395jaM1Iqyz3/p5UUlLg5Iw3BIK/QhcNpY+nzm1J+bDW4aM9IqSwWDIR56/lNXHTT8bjznCk/vohQVGalB2r3BcnJTX0bMkUyJi4MJRqMlMpiDdUewqEIH//8HKDvVVwHUnFZLg3VHla9toOTL5ya8uOrwUGH6ZTKUsYYGqo9zDl7ImMml6atHQXFbgpK3LylQ3XqMGgwUipLtdT7CAbCHb2idBGbcNIFU3hn4WZCwUOui1SqTzQYKZWl6qs9OHPsnHj+lHQ3hVM/OQ1Po581b+7STN6qXzQYKZWFqnc14WnyUzqiALsj/f/Gc86eiDvfqUN1qt90AoNSWeiNv1uXgAwbnp/WdsRPmDj+3Mm8/a8POWJKKSKSxlapbJT+r1RKqYS9/vcN5BW6cLkz5/vkaZ+cRmNNG97WoVlW4qtnPnTkV898SMuO95MGI6WyzK5NtexYX0NJRXp7RZ0df+5knC47zXXdpkpTWSKRsuMffvih64wzzphcVFQ0q7y8/NhrrrmmMhhMvJqEBiOlsszipzdgswnF5X0utZMS+UU5HDd/As31Xowx6W6O6qdEy47feOONleXl5aHq6uo1q1at2vDuu+8W/PSnPx2e6HE1GCmVRYwxvP70eo49YzxOV2pKRSTitE9OI9gext92WHXWVJxMLzu+Z8+enMsvv7wxLy/PVFZWhubPn9+ycePGhNujwUipLPLhyv1U7WhifoYmJT3pY1YGhuZ6HaobKJledvwLX/hCzZNPPjmstbXVtmPHDudrr71WfP755zcn+jwz5+ynUqpbset2Jh07EqfLzqmfnMZ518xKb6O6UFyeR35xDs31XX42ZpV7v/ivsTs/qO3zWOieD+tywZrI0Ndtxh9V4f3mfZ/I6rLj55xzTusjjzxSPmzYsOPC4TCXXHJJ/VVXXdWU6HG1Z6RUljDG8MY/NnL8uZMpKOnyS2pGKCrNpd0bpGZPwl+OVTdiZceXLFlSOG/evNZTTz21ddGiRQVLliwpTGfZ8XA4zIUXXjjl4x//eGNra+uqqqqq95uamuxf+tKXjkj0uNozUipLtDW301Dt4czLMnOILqZwWC5VO5pYsWgbH7t+drqb02+99Vg6i/WI/m/xDUOm7HhNTY2jurradeutt9bm5uaa3Nzc8HXXXVf/ox/9aAywN5HnqD0jpbJEU20buQWujEj/05OcXAdOl50Vi7b3vrLqkwULFrS+9957hX6/3zZp0qTgueee27pkyZLi5uZmR29lx7u7dRWIwCo7vmXLFvfDDz9c4vV6paey46NGjQqNGTMmcO+991YEg0Hq6ursjz76aNlRRx2V8ElDDUZKZYFIxNBc7+PkC6empW5RIkSEwmFuVi/eoYlTB0imlh0HeOqpp7YuWrSouKysbNaUKVOOdjgc5r777kuoVwk6TKdUVvA0+giHIsy/7Oh0N6VPCofl0nCgjk3L93H0KZXpbs6gkKllx0855RTfsmXLDntoUntGSmWBpjovdoeN2WdNSHdT+qSgxI3NLqx8VYfqVN9oz0ipDOf3Bmmp91EyPB+HM/MudO2K3WHjqOPHsGLRNq79wZnpbk5KJGPiwlCiPSOlMtzyl7cSiRhKMiz9T2/mnDOJLauraKptS3dTVBbQnpFSGe6Nf2zE4bSRX5yT7qb0SaysxIcr9/PonW+w6vUdLPzzqoMeU6qzlPWMRKRURJ4VkTYR2SUiV3SznojInSKyT0SaRWSxiGT2hRVKJYnPE2DZS1soLs/LuhpBk2eNpKg0lxWLtqW7KX1lNMFrckUiEQEiXT2WymG63wEBYARwJXBfN0HmMuAG4HSgFHgHeCxVjVQqk7y7cDPtvlDWDdEB2O02Zp89kVWvbs+KLN4i0hwIBDJ73nyW8/l8bhHp8vqmlAQjEckHLgV+YIzxGGOWAs8DV3ex+gRgqTFmuzEmDPwFmJ6KdiqVaZY88wFlowrJK8qOIbrO5p4zicaatqzI4h0Oh/+8f//+/Oi3dzWAIpGItLW15e7cudMVCoXu6GqdVJ0zmgqEjTGb45atAc7oYt0ngf8SkanADuBa4MWudioiNwI3AlRW6rUManBpa2ln+StbufBzc9i+tssvkxlvztlWpYPWRh+5Bd3WZ8sIkUjkvpaWltnr1q07DciOaYvZIyIi1aFQ6I7Zs2e/1NUKqQpGBUDnrInNQGEX61YBbwIfAmFgD3BWVzs1xtwP3A8wd+7czB8HUCoB77zwIcH2MGdcMp0v/vTcdDenX0pHFDBp5giqdzYxfGxxupvTozlz5gSwvvyqNEhVMPIARZ2WFQFdZZu9DTgeGIt1le9VwGsiMsMYo0VS1KAXKxfhzncxfGwRR50wJs0tOjxzzp7EU796m3Coy/PWSgGpm8CwGXCISHyGx2OBrtJNHAv8zRiz1xgTMsY8DAxDzxupISQUirDqte3M+9T0rJtF19ncBRPBgKf5kDybSnVISTAyxrQBzwA/EpF8ETkVuIiuZ8ktBy4TkREiYhORqwEnsDUVbVUqE7TUewkFI8y7JPu/g00/cSw2u9DaqMFIdS+VF71+CXgIqAHqgS8aYzaISCWwEZhujNkN/BQYDrwP5GMFoUuNMU0pbKtSadVU62XUhBKmzh6V7qYcNqfLTn6xG0+jD2NM1vf0VHKkLBgZYxqAi7tYvhtrgkPsvh/4cvSm1JATCobxNPm58IbZg+aDu7DETWuDj/3bGxkzqTTdzVEZSHPTKZVhmut9AJxxafYP0cUUDrPKpGsWb9UdDUZKZZjmOi8ut4OJx4xId1MGTE6uE5fbwcrsSQ2kUkyDkVIZpLXRh6fZn5W56HpTWOJmzZu7CAa0+qs6lAYjpTLIuwu3gIHistx0N2XAFQxz4/ME+GDZ3nQ3RWUgDUZKZZC3nt+E02XP+NQ5/VFQ7MbusGVTFm+VQhqMlMoQPk+Ala9up6gsd9AN0UG0+usJY3QSg+qSBiOlMsSKRdsI+EMUZ1G5iCcfX82Tj6/u8/pzzpnE1veraapt45YLHu1IfaSUVnpVKkMsfX4TxeV55GdwuYhY4PnMlccltF2swuvmVft55EeLWfWa9o7UwTQYKZUBAu0hlr24hXmXTOcbv/14uptzWHoKWJOOtaq/6lCd6kyH6ZTKAO8v3om3NcBpn5yW7qYkRWw4z263MfusCax6bUdWVH9VqaM9I6XSKHbOZPTEUvKKcph15oQ0tyj55pw9icV/30jZqIJBOWtQ9U9CwUhEzgVmEZdLDsAY88MBbJNSQ4oxhnde2MxJF0zB6cqcAqP9PT/Um47qr01+DUaqQ5+DkYj8FrgceB2IL3KnfW2lDkNbczvN9V5O/cTgHKLrrGxUIRNmDKd2bzPDj+hcc1MNVYn0jD4LzDLG7ElWY5QaiprrveTkOph7zsR0NyVlZp81kX/89l2qq7oq9qyGokQmMNQDTUlqh1JDkjGGlnofc8+ZhDt/6AxZzTnHqv4a9AXT3RSVIRLpGd0LPC4iPwEOxD9gjNF5mkr1g88TIBgIc2qaZ9H15fzQrp2N1NV5sdmEd9/eRVGxm/b2EA5H4pNyjzmlEgTavaF+t1kNLokEo/uiPztfBGGAzDnrqlQWaa7zIgInnj8l3U3p1oHqVpYs3s7OHY3YbFaaojff2HHQOq+8tJlTThtPfh97dy63A6fbQUB7Riqqz8HIGKPXJCk1gIwxNNf7KChxU1DiTndzDhEKRWhu8vPon1fidjuYf/YktmyuRUS45LKZtLT4+dc/N+DzhVizej8b1lVz/IljiURMR9DqiSvXQVuDn7ZmP/nFmff8VWolfJ2RiFQCY4C9OplBqf7bvu4AAX+IigycUbZlc13H5IITT67kxJMqyXE72LqlDgCXy055eT65uU5yc5185opZvLlkB28v3YXNJhQV5xAOR7Dbu/8O68p10oafDe/u5YTzJgPJm06uMl+fezsiMkpE3gC2As8A20RkiYiMTlrrlBrElj63Cci82kUrl+/ln/9Yj9NpZ+SoQuadOZEcd8/fW0vL8rjoUzO48prZOJ02mhr9PHT/MjZ9UNNtpgVXdJ9rl+4a8Oegsk+i54zWAB8zxrSJSD5wF/AH4JPJaJxSg9mb//yA/OIcHM7MOOVqjOG1RVtZuXwvk6eU4fMF+zTcFm/0mCIqhufj94cA4V//3MjykYWEIxHcnQKa2ASn28H6t3YP4LNQ2SqR80CnAd80xrQBRH/eCpySjIYpNZjt2lTLns31GVMuIhIx1Nd5Wbl8L7PnjuGiS45OOBDFiAi5uU6uvWEuF3x8Gl5vgNqaNhrqvUQiB/eSXLkONq+uwt8WGIinobJYIsGoEZjeadmR6LVHSiVs6XObEIHisvQHI2MMdXVefL4Q88+exNkLpvQ7EMWz2YSjjxnJ5246gcKiHNragrzy4uaDhu1cuQ7CoQgb3tNS5ENdIsN0PwMWiciDwC5gHHA98INkNEypwWzpc5uYfuIR/OKV61J+7M6TBDasq6bdH6JkmJu5J4wd8OM5HHZKStwIsHZNFQ6HjbMWWBMWnLlObHZh3dJdzDlr6GSgUIdKZGr3n0RkG3AFMBPYD3zWGPNashqn1GBzywWP0u4Lsn3dAW76yYJ0NwefN8ji17bhyrFTkOSkpUXFOUydVsGKZXtxOG2c/z8nISK8/cBq1i7V80ZDXUJTu6OBR4OPUoehud4HwKmfPDLNLYE3Fm+jvT1MxfB8RHofmjucKdciwplnTSIUjLDs3T0UFedQXOxm5qnj+Od9y2jv4gJYneo9dPQYjETke8aY/43+/qPu1tMSEkr1XXOdl6lzRjOisiSt7di7p4l1a6wLVaurWlJyTBHhnPOmEApFWL+uGhHh2NMqefrX7/DB8n0svPcdQIPPUNRbz+iIuN8HfjBZqSEm4A/h86S/oqsxhlde3ExRUQ6nnDaeZ55ee8g6fQkI/QkaIsJ5HzuSbdvqaW7yE3LbEEGneA9xPQYjY8wX436/PvnNUWpwa663SoGddlF6g1Fra4DmJj+f+vTRuNJQ0M9mE8rKcqk3hrff2cPoyWV68esQl0gGhoZultf0cftSEXlWRNpEZJeIXNHDuhNF5N8i0ioidSLys762U6lM1lzvw53nZMyk0rS1IRSK0NLsZ/KUMiZPKU9bO0SEsvI8JkwsxVaSw8Z392IiWqtzqErkOiNn5wUi4qTvGbt/BwSAEcCVwH0iMqOLfbqAV7AmSozEGir8SwLtVCoj1Ve34m1pT+mFrk8+vrpjEgBYw3ONDdYEirMXpD9TuIhw0SUzmDp3NMFAGFfh0KnppA7W62w6EXkTq0yEW0SWdHr4CODtPuwjH7gUONoY4wGWisjzwNXAdzqtfh2w3xjzi7hlhw5oK5Vl3v73ZiC9ueg2b6rF7w9RXOKmKC5T9kBOGEh0X06nnRv+51RevW8FzqIcXEOoyKD6SF+mdj8ACHA88GDccoNVZK8vU72nAmFjzOa4ZWuAM7pY9yRgp4gsjB5zPfBVY8y6PhxHqYz11vObyMl1kJN3yCBDSvh9QRa9sgWn00ZhhvVAykcWMm56BQ07mpj1XzNobPAyrDT92SlU6vQajIwxjwCIyLvGmE39PE4B0NxpWTNQ2MW6RwDzsZKvvgrcDDwnItOMMQclsBKRG4EbASorK/vZNKWSr6XBx5olOykbVdin63mS4Y3F2/F5gwwfUZC2NvRk5mnj+PcDK4mEwjzz9DquvHZ2l+vptUeDU5/PGRljNonICBH5hIhcLyI3xG592NwDdC7aUgS0drGuD1hqjFkYDT4/B8qAo7po0/3GmLnGmLkVFRV9fSpKpdx7CzcTCZu0DdHt3dPE2vermHP8EWmZPdcXM0+rxEQMWxftoKnJz7/+ubHb8hNq8ElkNt3FwDbgR8Afga9Gf17dh803Aw4RiT9jeiywoYt112INASo1aCx9/kMqjigiN8kpd7pijOGlhZspKs7h1NMnpPz4YPVieuvJzDxtHAAtVR7OOW8KO3c00tTkT0XzVAZIZDbdncD1xpjjgLbozxuBlb1tGC038QzwIxHJF5FTgYuAx7pY/S/ASSJyjojYga8DdcAHCbRVqYzh8wRY+eo2Tv3ktLQMj7W0tNNQ72XBeVMztlcEUFKRjyPHTsAb5NhZo5k9dwye1gAej5aXGAoSCUaVxpinOy17BLimj9t/CcgFaoC/Al80xmwQkUoR8UTLmWOM+RC4CqtoXyNW0Ppk5/NFSmWL5a9sJdge5rRPpD4XXTAYpqW5naOmD2fipLKUHz9ROXlOAr4Q/rYA88+ehNvtoLHBR22tJ91NU0mWSKLUGhEZYYw5gDXb7WSsHkufvmoZYxqAi7tYvhtrgkP8smewelJKZb23nt9EcXke008eyz0L+/rdbWA0N/kRgfnnTO5Ylskn/l25Ttoa/ax7ew/HL5hEaVku1VUeXl64mSuuPi4jJ16ogZFIz+hPWNVeAX4JvI41Pfv3A90opQaLgD/Eey9u5ZSPH4ndnsi/W//EX+Ratb8Fny9EYVEO+Vly7Y4r1wECq17bDoDdbqOkxM3+fS28v3p/mlunkimRekY/jfv9URFZDOQbY/RcjlJduOWCR2lp8FmJUdOQi+7NN3ZgswmFhTkpP3ZfdNVDE5vgcjtYvXhHx7K8fCcFRTksWbw9remLVHL1+6uaMWa3BiKletZc5yW/OIdj541P6XF37Wxk185GiopyBqSEeCq58pzsWF9DY411nkhEOPe8qUTChldf2ZLm1qlk6a2e0R76MM3aGKNXnCrViTGGlgYfZ356Os4UzmIzxvDmGzsoLMyhIMMyLfRFTp4TT72P1Yt3diwbVprHKaeNZ8ni7ZSV55GXpiwWKnl6G6a7KiWtUGoQamtuJxyKcOonUjtE5/eFqKvzcu4FU/lgw4GUHnsgOHLsFAxzs/r1HYw57aOSanNPOIIPNh6god6L251QkWqVBXqrZ/RGqhqi1GDTXO9FbMKccyal7JjGGJqb/ZQMy+XoY0ZmXTC64JsnA7Bl4TZWvb6dNe/uRkT4zJXHYbfbOO+CI/nLI6to1othB51EMjDkiMj/ish2EWmOLjtXRL6SvOYplZ0iEUNzvY/CYW7cKRxS8nqDBIMRTj19fEpm7w20WKaG2WdNoG5fK+Fg5KDHR40uoqDARSAQJhyOdLMXlY0S6ev+EhiDVYtoYXTZhujy3w5wu5TKahve3k0oEE5p7aJwOEJLcztOp42jpg8HMvuaop4cd6aVtqjdG8TR6XxbcYkbEbIy2KruJRKMPgVMNsa0iUgEwBizT0TGJKdpSmWv1/++AbEJxaWpS4y6fm01oVCE8vK8rL84dPTEUkaOL6Gxto38EvdBj2Xb7EDVN4l8tQjQKXiJSAVQP6AtUirLhYJh3nz2A4pKc7El8dt7/AWuwWCYt5buxOWy484dHCf3jztzAgFfSDN3DxGJ/Kc8DTwiIhMARGQU1vDck8lomFLZatVrO2hp8DGsInVDdKtW7KXNE6CkxJ31vaKY4+ZPwEQMQX8o3U1RKZBIMPousBNYB5QAW4D9wB0D3iqlstjiv6+nYJibgmGpGaLz+YK8985uJk4qJWcQTXmeFb1Q2GnTc0NDQZ/+ytFSDt8Hvm2MKQBGAIXGmG9oNm2lPuL3Bnn735s57ZPTUnZu4713dtPeHmbemRNTcrxUKS7PIzffiUencQ8JffoaZYwJi8iXgduj92uT2SilstV7L27B5wkw/7IZzDoj+YXsQqEIq1fuY/rRI6gYXtD7BlmmoMRN7f5WfJ5AWgoTqtRJpP/7CPCFZDVEqcFg8dMbKB1ZwDHRqqXJ1tLsJxIxnHb6+JQcL9UKStxgYP07e9LdFJVkiQwwnwB8VURuBQ7KWWeMmTfQDVMqm9xywaOEQxE2r6ri4/89JyXXwASDYdragsyeO4biktRNIU+l/KIcRGDNGzs5fkHqMlmo1EskGP0pelNKdaG5zkswEObMT89IzfGihfNOOiU1vbB0sNlt5BXmsPqNHb2vrLJan4JRdALD9cB5xpj25DZJqezUVOdl1MRhHDln9IDvO3Y9USyjwt49Tfh8IYqKDy6cl60ZF3pSUOJm25pqWhp8FKXwImKVWn0aSzDGhIEJwOC4gEGpARYMhPE0+Tnz0ulJv84nEjEsenkLdnvmFs4bSAXFORgD65buSndTVBIlMrB9B/AHERknInYRscVuyWqcUtmiuc4LwJmXHZ30Y615fz+1NW2UDMsdEqlxcgtzyMlz8v6SncBHyVTV4JLIOaMHoj+vjlsmWBMZUlc5TKkM1FjThjvPyfijKpJ6HJ83yNI3djC2sgRjhkbWaptNOPqUsbz/xs50N0UlUSK9mgnR28S4W+y+UkPWphX78HkClI5M/nU+S9/cQXt7iLMXTB40aX/6Yta88ezeVEfDAU+6m6KSpM89I2PMLoDosNwI4IAZKl/NlOrB839cjs0uDBuen9TjBAJh1qzez6zZYwblBa5duWfhNQBsWV0FwPtv7OSsy5M/FKpSL5HiekUi8ijgB/YBPhF5RESKk9Y61WcHqls5UN2a7mYMOQ0HPLzxj40MG56P3ZG806fGGJoafbjdjkF7gWtPJs4cQUGJmzU6VDdoJfLf8xsgHzgayAWOAfKiy1UK9TXwaIBKvv88tIpQMEL56MKkHsfnDdLeHub0Mybizk1d5dhMYbfbmHnauI5JDGrwSWQCw/nARGOMN3p/s4hcD2wb+GapmFgwGTFyYD7sBnp/Q1kwEObfD6zk+HMn0e4NJu044XCEpiY/TqeNY44dlbTjZLpZZ4zn7X9/yM3zH8LldnQM4anBIZFg5AcqgPjJ/uWAXgSrhqQ3n91IY00bF33hhKSmqtmw/gDhsGFY6cFTuYfa9OZZZ44HwNPkT8lkEZVaiU7tfkVEfoEVkMYB30BTBKkh5pYLHgWg3R/iiMmlzDl7YCeUxmdbiEQivPv2LpwuO+5BVKuoPyqPLGfY8Hw8zRqMBqNE3t3/izVx4UpgNFZhvZ8ZYx5MRsNUauiwXf94W9vZuuYAX/r5eUm98HTjhhqam/yUl+cNqancXRERjj1jPEv/+YGWIh+EEpnA8GvgQ2PMOcaY6caYc4APRORXyWna0KSTDrJD3f5W8gpdLLhiZtKOEYkY3n1rFxXD83HnDu1eUcyseeMJBSO0+7QU+WCTSDD6LLCi07KVwBV92VhESkXkWRFpE5FdItLrdiLymogYEdH/RJUxgu0hmuq8nHv1LPKSmBtu0wc1NDb6OPnU8UO+VxRzXNx5IzW4JBKMukr7Y09gH78DAlgXzF4J3Cci3ebaF5ErSWwYUamUqK/2gIFP3jg3accwxuoVlZfnMfXI8qQdJ9uMHD8MZ44dT7MGo8EmkWD0JvDjWGLU6M/bo8t7JCL5wKXAD4wxHmPMUuB5Ds5zF79+MXAbcGsC7VMDRIcKu+dtbae+ykNRaS5jJpUm7Tg+X4j6eq/2irpQUOLG5wnoeaNBJpGex83Av4EqEdkFVAJVwCf6sO1UIGyM2Ry3bA1wRjfr3wXcB1Qn0D6lku7fD6wkHIowfGxR0o5hjKGl2U9paS5TpyU38Wo2GjW+BJvdpkF6kEkkN91eEZmNVX58LFbp8WV9zE9XADR3WtYMHDKFS0TmAqdiBb8jetqpiNwI3AhQWVnZh2Zklr7MZDPGYIzB6w1gDJiIIRSyXvKGBi+RcIRw2BAIhAChpcVPbq4TRxJT0wxV7b4gz/z2PQpK3Ek9V+TzhQgGI5x06riOmXpD7ZqinjicWiRgMEronEw08LwbvSXCA3T+KlkEHDQWFB36+z1wszEm1Ns3H2PM/cD9AHPnzh1UfXZjDB5PgEAgDEBdbdsh67S1BbDbBJtdEJGO/GVNjT4cDhuRiMFut2GMOexvkToFHF589H0aa9qYeMzwpB0jErF6RQ6HjaOmJ+84SmWaVE0Q2Aw4RGSKMWZLdNmxwIZO6xUBc4G/RT88Y1+B9orIZcaYXs9PZTtjDJGIYf++ZsJhg4iVl6u8PB+xWUGnrtZKoz9y1EfxPRYsysrz8PlC+H1BfL4gkUiY6qpWCgtzyIsrT60SEwyEefpX7zDjpCOw2ZM3PPTBhgMEgxHKynKx2bR3q4aOlAQjY0ybiDwD/EhEPg/MAi4CTum0ajPWBbUxY4FlwBygNgVNTSuv96OeUE6OnbKyXJqafIgIrpyP/lQ99XIcDjuFhXYKC3OormohErE6jA0NXhobfdhs6IdcAmLZFs7+zDHU7m3h5l9/jOPPnTxg+4/PthAKRVj65g6cThu5eUMvGaoa2lI5dfpLwENADVAPfNEYs0FEKoGNwHRjzG7iJi2IiDv66wFjTFZf5dbTMFc4HKGx0Ye3LYAIOBw2ho8oRESQw5jCKiLY7cLwEQW0t4fxtPrxeoOEw2FqDrRSUJgzIEN4g50xhr/98m0mzxrJ3CTmoFuzej8tze1UVGi2BTX0pCwYGWMagIu7WL4ba4JDV9vsxCptPmiFwxGqqlqIhA1FxW78vqAVhAbww0hEcLsduN0FVFe1EA4bgsFwx3mogTqvNFg11XnZv62RH/zl00l7jQLtId59exeV40qIRLRmpRp69KLSNAmFIgSDYSIRg9NpZ/jwPFwuBwf8ye0AiggOh9Vb8vmC1Nd5rYC4v4Wy8nxycvr+lhgKkxqMMdTuaaFyWjmnfOLIpB1n+bI9eL1B5p05kTde16osPdHSEYOTBqMk6G1IrqXFT2uLVXnDbhdGjipMea9ERMjLc9HqaicSiRCJGA5Ut1JYpEN38Vrqffi9QT7zzVOTlhA1HI6wfNlepk6rYNTo5F2/pFQm02CUIsYYwmHD/n0tGGPIz3cRDIYHfEiuP2w2GyNGFtDU6KW1pT163kqv5fB5Auzf0Yg7z8mZn+42c9Vha2luJxQMc/q8CUk7hlKZToPRYeptqCoSMbR52jtmyblznZSU5OJy2TMq5Y7NJpSW5ZOb56K2xhM9p+ShuCQX5xC9yPDxu98k2B6mcmY59iRdRBwKRfB4AsycNYrSsjxAL3BVQ5MGoyQxxtDa4qelxd9xvZDDYWf48MwuCpab68TlshMOR/D5gni9QQoKht7Q3Y4NNfzjt+9SOiKf/KLkZVtobvIjAqecNj5px1AqG2gwGmCBQJhQKBxN0eMjJ8dBWbmb5ixKeW9NcrBTXpFPc5Mfjyd2fqv3WXfZPqnhlgsetYZUQ4bCklxGji9J2rEOVLfi9QYpLMqhMInphZTKBnr1Y4K6ymgdiRiamnzs39/cMXVaxJqxNmJkIW53/y9gHDGyMG0f7Ha7jdKyPEaNLsJmE8LhCNXVrdE8eINXw4E2Nr63l8/feXbS8qAZY1j82jZsNqEoiT0vpbKF9owOQyAQIhAIY4wh2BwmJ8dBYWkObZ5A9NqexINQXwJPV+skM2A5nXacTmvoLhyKUF3VSnGxe1AO3YWCYap3NnHMqZUsuHIm51517IDuP5Zx4cSTKtm9q4mSYe6kli1XKltoMOpBd0NOxhjq6trwtgUAqwcxclQhdrvV0fS2Bfu0/2wbyrLbbYwYmU9jg5fmZv+gnHVXtaOJcDjCV391QdICrTGGxa9vp6TETUGB5gtUCnSYrkNfCsqFQhFCoTCBQBifN0BRUQ4ulx2Hw9YRiLqTjuG2ZBzTbrdRXlFAeUU+xkAwGKYtGpSz3arXt9NY00bFmCLGJbGOkLctSF1tG/POnDjoepZK9deQDEaJVjINBsLU17V1ZNK22YRRo4spGTZ0c4jl5blwueyIQH1dG42N3qyuvOlp8nPvF/9FTq6DEUksnBeJGJqb/YwaVaiF85SKo8N0PYhETEfuOBEoKMwh0B6KzjbrPo5n8vDbQLZNRHA67bhcDlpb2jvOn8UH6EyeXRfLyH3Pwmv4/a0v0VDtYeLRw7H10ss9HJ7WdsJhwxlnTUJE9JoipaKGZM+oN4FAiJoDrR2544qL3YweU0xp6aE9oXTOdssEIkJpWR6lZXm0+60JHbGyFdli6XObePWv67ji1tOSWsHV6w3Q0tKOO9fB2MqSpB1HqWykwSiOMYb6ujaqq1ppD4Sx2224XHaKS3J7PSeUrQYqmBYU5DAyup9gMEygPTumfwcDYX5z83+YctwoPnvLaUk91puLd2AMlJS4e19ZqSFmcH7CJigYDFsfoAHrZHxhUQ6jRxfhcNiG7Dmh/nDlOHC5rNl1Bw600p7kDOSHyxjD3q0NeD3t3Hr/RUm7pghg/75m1q6porDQNWTTKynVkyF7zsgYg9cboLW1veND02YTRo4qytrzQQPlcJ6jiOBy2YlEDDU1rTgctoytLNtY00Zrg4+b7l5A5bRyIDnlCSIRw6KXt1BQ6KKoWHtFSnUlMz8lkiwcjhAIWMXlQsEIxSVuXC7rws6eApHqGxFhxMhCHA4bwWCEcPjgYnGJzmYcCLdc8GjHhAWAbWur2b+9kfziHC7+4glJPfaa1fs5UO1h/tmT9QJXpboxZHtGIkJZeR65uU5EBL/v0CGlodAL6qtEXwu73Sqdvn9fs5WZurWdggzJv1a9q4nvX/IkdoeNyqllSQsQTz6+mnA4Qn29j8rxJRw5rYL3V+1LyrGUynZDshtgs1lDSXl5Lj0nlER2uw2n046I0NDgpaEh/dcitdR7+d6n/kqgPcSEGcNxJlDZtj+am/wEA2HOOXeKvteU6sGQ7Bnph8LA6EtvyboWyUZOjoPWVquIXLpy2kXCEX54+VMc2N3ET567kkfvXJzU47W3h2hrC3LCSWMpK8sHtFaRUt0ZksFIpZaIMKw0D6fTTkODt8ucdsm+ONYYw+4P62lt9PG9xy7lmFMrk3KcmHA4QmODD7tdOPnUcUk9llKDgQajKD0/NDB6eh0LCnNwOG3UHLAqyba2+CkozEl6L8kYw/7tjbQ0+PjSz8/j9IuOSurxAN54fTvBYITy8jxcLv03U6o3+l+iUsrttirJBoMRGht9tLUFOsptJ8vz96+gvspD+ZhCLrrp+I7lyZjGDbBlcx0rl++loMBFbl7/a1kpNZRoMFIpFzuPVFjkprHBS3VVK3a7DFiWi/iccytf284fvv0yhaW5jEpi1daYpiYfC1/YxMiRhTicem5Sqb4aksFIh+RSp7vXWkTIz3eRm+ugqdGHxxMgHA7T5mknL9+a5Xi455H2bqnnrmufofLIcvKLkjsc+OTjqzHGEAoZMIZPXDydF/+zKWnHU2qwGZJTu1XmsNlslJblR6eAQ329lwMHWrvNbdfXC2ZDoQi3Xf437HYbd/ztcuwpuJi5qclPdVUr5184jZJhuXzmyuN09pxSfTQke0Yq89hsVjmKgsIcmhp9VFe3YrP1XKoDup6FZ4xh96Y6/G0B7v7XVYwcPyypbQfweoN4WgPMnjuGqUdqnSKlEqXBSKVcT0N3BQU55OW5aG72ddRIOlDdSn6Bi7y83kt0RyKG3Zvr8TT5+cZvL0z6FG6A/ftaaKj34nLZOfOsSUk/nlKDkQYjlXFsNmHYMKs+UjhsFThsqPfS2OBFxJroELtw9q4r/w7Ar1+9nkjE8H9f/w/NtV7KRhdy/rUfDZF968+f6tOxO/e0uup5xS+rq23jH0+txWYTiopzBm2pEaWSTYORyghd9ZasirrC8BEFVjYDT4C2tgCRSJiq/S3k57soKsujpd5LJGL4zc3/YeHDqykbVUD5mN4nPfR3gsRLCz8E4KJPzeDpJ9dgt9vIL3AdEogyucqtUpkmZV/jRKRURJ4VkTYR2SUiV3Sz3rUislJEWkRkr4j8TEQ0aA5hIoLb7aSsPB+Xy8qsbnfYaG7289/3nMtl3zqVe7/0LxY+vJrP3nIq5UcUJf1C2nA4wlNPriUYinDZZ2b2Odt7OjKWK5UNUvkh/zsgAIwAZgEviMgaY8yGTuvlAV8H3gMqgOeBbwF3p6ylKiN011uy24URIwoJhcL87dfvsGnZPpb9Zwsfu3EO5//3HCI2YcvK/UlrVyRiaGr0A3D5Z4+lYngB511wZL/21dswoFJDRUqCkYjkA5cCRxtjPMBSEXkeuBr4Tvy6xpj74u7uE5HHgfmpaKfKfN+KXtD62OovEwkb/vmb92hp8HH5/5zCJV8/CZ83yDlXH8tZVxxDTY2HgvyBzYJw4EArDQ0+wqEIn778GMYcUTxg+1ZqKEtVz2gqEDbGbI5btgY4ow/bzgM6955Uhnry8dXAR9mpO98/3GW26PBbS4OPH13xNC0NPvLLcrnh9vmICCUlhuqqViJ2IRgIUecLdtQrstsl4YzhsfNDV183h+XL9vDm4h0AlAxzM2FSWZ/3c7i0t6QGu1QFowKgudOyZqDH/ywRuR6YC3y+m8dvBG4EqKxM/hRedbCugkUqhIJh/mfBw1TvbKJoRD7uwo/qUokIr7xkfee5+ro5+P0h2jzteFrbiURs7NndRE6OA7fbQSQSQaT3ABUOR3jqr2vYvauJKVPLaWryZUTFVg1QajBJ1QQGD1DUaVkR0O2ZXBG5GOs80QXGmLqu1jHG3G+MmWuMmVtRoRcaDgVBf4iGPS001rbxk+euwF3Y/bVHIkJurpPyigLeWrqLdWuqKCzMIWIMzc1+gkGr/PzePU1U7W+mpsbD9m31bPqghtbWdlpb/JiIoaHeR9X+Fs67YCoXXTIjIwKRUoNNqnpGmwGHiEwxxmyJLjuWbobfROR84E/AhcaYdSlqo+pBOnpBC+9956BjLn1uEw17W7DZbfzyleuoPLK8Y9iuN1b5by/DSq0M4ZFwhP+8sIm8XCfHzRlDKBQhFIowfEQBDoeNte/vZ+2aKjytASqG53P8iWMpLMxh395m5p9tXdhac6AVsQk2EYLBMCJCoD2E02XXAo5KJSglwcgY0yYizwA/EpHPY82muwg4pfO6InIW8DjwKWPMslS0Tw0cYwzGQGOjD583gM8bJBwxvPfObrzeAF5vkNqaNiIRw58fWE4wGCYYCOPzBQH45T1LOvY1/ZJpAPziZ2+wfeketry+k6JRBUxdMJHNOxtoaQ9xztdO6PO06ng2u42Gei8NwFkLphCJGHy+IP9+fiNtngDBYISysjxKjnBTUprL2MoSImFDOBLB5w1ijCESMUTCEUzEmmEHhurqVkQgx+0kFIpgswmRiElZb0qH7lS2SuXU7i8BDwE1QD3wRWPMBhGpBDYC040xu4EfAMXAf+K+Xb5pjLkghW1VvYh9eAcDYV5+8UOam/y0NPtpbPRhDDzwh/cOWn/J4u04HDby8p2EIwa7TRg2LBen047DaWP7tnoEOGrGiI5tXnh4NZFwhHBjO1uW7GLGGeOwl7lxuOysX3+A1aus6dsi8NCfllFc7Kao2I3fH0IEXn91Ky0t7bQ0+ak54AHgd795C6fTjtNpp6nJh4kYfvurtzqCIVgZIM6/8EhmHD2SvzyyktqaNoqK3B2PHwh2PR3bGENhkZt2fxC/P0Q4HCEchr17mnA67bhy7ITDEUSs816pSN6qVLZIWTAyxjQAF3exfDfWBIfYfZ3GnWY9DckFg2HWr6tmxbK9NDX6ANi8qY7iYjflwwsIBMPY7TZOOXUcuXkuli7Zjt1u4/LPzsLlsh+0/4svPfqQY54x/6Pcbn+6eSFN+1sJtoe5/rb5/Nc3T+G6OdbM/9tWfJGGBi/P/WM9gWCY0tI8mpv97N/Xgt9vZfx+f9V+iopyKCp243Y7EBEmTyknFIwQDIZpafYjDmHqtAry8pzk5TlZsXwvOTl2jpk5CiCh64diZTHy863zWNVVLUQihrx8F4H2ED5vMNqDgv37W+K2g7q6NhwOGw6HrWNihVJDiWY2UH0SDkfweAL88Xfv4vMFGTW6ELs9D7fbwRVXz+5YLxZUjo5+mC9/bzdARyDqTufzQ8te2kr97mZMxHDbE5dxyicODgo2m1Benk9+gYt8Dg5sTzy2CmMMV1w9u+ND/dE/rwAODi6xZeeeP7Vj2fp11X18RXoXu0C3pCQXsIYw43tQoVAET2s7xhja/SG84chB21ftbyE310lurjPhKekxOmynsoUGI9WjcDjCqpX7qNrfijEwaXIZJ5w0ljFHFPO3J94f8OO1Nvr443de4ZUn1uJw2Skek39QILrgmyf3ug/r/Iwc9OHdWzBMBRHpuBUU5ADQHu3FjRhZiDGGcChCba2HSMR6Hi0tflpa/B3b19V6sNlt2G0SHfLTHpQaHDQYDXE9Dcnt2FbPa4u20tDgw+12UFLi5pLLjklaW/xtAW468X4aazxccetpvPqPDUiKp1H3N63PQBARHE47NpsNm80KUJFwBJ8/RGODF2MMgUCYcDiIMR9td6C6hfz8HPLyey+xYa2vvSWVeTQYqUMEAmGam/38/al1lAzL5ZJPH82y93YP2LfwlX+xZuvHAmD1rib8TX6aar1MmDGcO/52OVOOG8X6t3cfsu1ATi2/5vq5A7avZLHZbeTnu/C0tgMfBZDYkF84bAhHDA0NXhobPyqxoVS20WA0xHTXEwqHI2zZXMfqlfs4UO1BBObNn8icuUfgcNhYvmxPn/bfeb89BY/mOi9P/vwt/vWnFdY1PmOL+L8ln8OZwJBaX4NTX9fLlt5CbLgvVmIjEAjT5gng8bQTiYSprfFQXOLufUdR2dRbyqa2qr7TYDTEhcMRPK3WxIS2tgDFJW6KS9zk5zs58aSeUyz1t5cSCUeo29/Kdcf+Dr8nwLlXHcvOD2pw5TgOCkT3LLymX/sfyN5TNnzgiQg5OQ5ychwEAlZBQn97CF+VVbpde0oqG2gwGqLq69pYvmwP+/dZ3zInTirluNljGD+xlKf++v6AHeeWaJbtexZeQzgU4aXH3mfTyipCgTAnXziV62+fz7hpFR3rZYNMDlCx3lLF8AJaW/w0N/utYoRVLRTku/p8XilTaC9o6NBgNIh1HpKLTSFuaW3noT8tx+GwKpQWFrq49PKZPe6rr72N+OATY4zh7X99yEO3v8aezfXkFboYd2QZtz95eX+eVkbKtA9Lm00oLsnF5wsSDluzHRobfTQ2Wkle+5IZQmstqVQaksEoXdmm0yUSMWzZXMfy93ZTU9OGzSacfOo4Zs8dw/PPJqc6hzGGvVsaqK9qpbGmjTuueJqxU8u47YnLeOZ37x4yGaK/Q3KZLBM+sGM9pREjCwkErNLtra3tRCKGvXuayHE7yM11EokY+jM/5XCCkwY7FW9IBqPBqKsAG4kYvG0BHrx/GU2NPopL3JQMc5Of7+K0eRO63Vd/e0GeJj/11R48TX4+O/lXNNa0AeDMsXPzbz7GeVfPwu6w8ezv3+t2n4Nd5w/ZVH7oulwOXKUO2ttDGAPuXAd+X7Ajk4YINDX6yMt34nT2/7qsZAaZWEJbML2uq7LLkAxG/f0WmCl66tlZU349rF9XTdX+ViIRw8hRhXzy4ulMObJiQM8HxWxbW83z96/g9afW0+6zslaf/qmjmHnaOBY+vAqX28HHrv8oS8Ng7AUNpGQHKGsmHgwblgfDIBQKU3PAQzhsOi6ydThsHVkfAoEQDoc9baUzjLHa5fUGCLSHo8+BfmelUJlpSAaj1pZ2PJ4A77y1i1mzR5Ob68zYobu+tiscitDmDfLwgyuoq23Dbhfcbgf5BS6uunZ2t/+0/e0FtbW003jAQ321hy+d+gA5uQ7mX340W9dUk5vv5NsPXAzAa3/TCiADoasANVBBy+GwY7fbsNuhvCIfrzeI1xuIZoewKucC1gXI0QDQ2ODF7rB1ZIEYqMAQy/re0uIn0B4mELB6cYGAD6fTTnGJG2+bVb1XA9HgMiSDkdvtIBAIs3TJDt57ZxfHHDuKUChySCmCdASovh7TGMO+vc1s31bP9m0NHVmpR40uYsH5U5k2rYJ/PrMeYMD+acOhCK8+uY4lz37Ayle3EWwP43I7uOnuBSy4YiaFw3KzalbcYHS4Acput1FYmENhYQ7VVS0YA8XFVh69UDiCty2AMQaPp/2gLBB7djfhcNhwOu2EQlZtp7a2ACLW+y8SMRhjaG1tJxQKEw5FCASsVEj79jZhzEeBCKzhQrvdCjg2m3UtVWzo0O8LHdZzVJlpSAajHLeDCreDsxdMYfl7e3h/1X4iEYPb7eD91fuZOKn0oJIB6RaJGKr2t1Bf10ZdnZfamjYCgTBPPLYaERhzRDHFxTnk5jm56trZPe6rc/aDrma/xS8LBcMse2krOz+opbXBx4Z391I+ppCPf34Oa9/cRV6hi0u+fGK3x9MhufQ6nOAUG86Lnw4eDFjDZMNHFBCJGGoOeDDGykweDIYJBsPR2XuG+rq2Q/bZ2OAFiH7xs/bvdjutY9nA2xZEBIaPKMThsHWcazqcc1gqOwzJYBRTMbyAj33iKE4/YwKPP7oKrzfIKy9uBqzhivb2ELluB5FIBJvN6jUNZG+pq329+Mt3cZe4mTK1guqqVqr2txIKRfjLI6sAsNutb4ruXAdnnzOF8ROG4c51cu3s3wNw9XUfpbjpnAk7ET5PgPu+/TKvP7We5jovDqeNslGFfPeRSzjqhCOw2aTLXpAGn8w3EMN7sYzksaS0sczkYJXOACgrz8dEDBFjokHIuv4p1uOJBZqy8vyObQPt1rL+FExU2W1IBqPOH86FRW5KhuVSXOJmwflHsn2rNfRVV9tGa0s7v/3120yYMIyJk8oIhyMHXdEeCwKPrPpSwssW3vsOYhdOOmUce3Y3sWd3E0d9aho2m/Daoq3k57to2NWMr9HHjT+cT3l5PiXD3Fw/9w8AfPlrp/b4PIePKOj2MU+TnwO7m2iu9xLwh7jv2y9Ts7uZA7ub2LmhhnDYsHNjLSd9bCoLrpzJ0796GxFhxkljezymyk4DOWkiNiwc35uJfZnTIKO6MySDUXdErBo55eX5nHBSJU88tgq/P8So0UVs31bPpg9qAQi0Bfj739ZSXpFPyfhi2lsCNDX5yMvrfUpsQ72XqqoWqqtamTh/PO7SXJ5+ci0iMHJUIfWb6/HW+/jJE5dRUJjTUUxu6pEVPe43FniCgTA1e5qp2dNMQ7WHYDDMfd9+mea6NppqvWxeVUWgPcSlY39+0PYvPrKaEZUljKgspqm2DXeek1+9dgNFpdY33r//+p1Djqm9oMEtU6/1ydR2qcMjxgyO+fpz5841K1asGNB9xvdmYlOm7/nmf3AXu5kydzQN9d6Oq9tjHA4b3mY/kWCEsVPLOpbv2VyPM8+JPZp7zem00bTfg7fOyxdvm8/oMUV4W9r58pkPYYDbHrkUgDuu/QcYwy2//wT+tiDtviCP3bWESCjC6Z+a3hFk1r21i2B7mFAwTOc/aV6hi+LyPErK89m7rQFXjp1LvnISI8YV8+TP38KVY+cXr1zX8Y22q/NISh2Ovl57lMSLXnXqXYYbkj2j3k7ax8QPc4kII0cVUrepHmMM37nrXKp2NnLHtc+ACCecO4m2lnaaPQF2bKwBoHZdTcf29VWtmLBh4lEV5OY5EbeD1c9tIhI2/O+GWmr2NNMeN0vo5rP+fFCb/98nnzjkeTzz23cpqcinuDwPm91GQYmTj39+DsPHFjNiXDEP/fA1HC4797507SHP89NfOwmA5/+wvOP5KZUsfZ2arr2eoWtIBqPetDX72f1hHQ0HPATbw9x360s01bbRVOelblcz4VCET1fee9A2//nDCkQgJ89JwB8CEcQf7ni8odqaer0z+FFp6YAvhM0ujDuqghPOm8KIccW88OBKRITP//hsAB74wauIwJfvvQB3vhN3nouff+E57HYb97587SG9mav+37yO/bvc/fvzao9IKZVqQzIYNda24Wn08/Obnu9YtmdzPcFAmCun/Zq6aCbrmFeeWEtxeR7FZfnkFbhw5ti59KsnMaKyhL/+fCkul527/30VLrcDEenzBIZYAPnBXz7dsWzpPz8A4PhzJwPw1C/fBuCYUz8q5+DKsf5s/enNdA40GniUUplgSAajgC+Ep9nP2qW7OpZ5mv3YHTaOPX0y46ZVUDmtgr/e8ybOHEeXw1yXftUa5nruD8sAyMl1dqzT1Sy2nma2JYsGGqVUthiSwWhEZTEjKou7PGd0658u7lj2zG/fTXXTDtFVQOnrMqWUyhZDMhilgwYVpZTqnk7tVkoNBTpdNMPp5dBKKaXSToORUkqptNNgpJRSKu00GCmllEo7DUZKKaXSLmXBSERKReRZEWkTkV0ickUP635DRKpFpFlEHhKRnFS1UymlVOqlsmf0OyAAjACuBO4TkRmdVxKR84DvAGcD44GJwB2pa6ZSSqlUS0kwEpF84FLgB8YYjzFmKfA8cHUXq18LPGiM2WCMaQR+DFyXinYqpZRKj1T1jKYCYWPM5rhla4BDekbRZWs6rTdCRMq6WFcppdQgkKpgVAA0d1rWDHRVvKTzurHfD1lXRG4UkRUisqK2tnZAGqqUUir1UpWbzgMUdVpWBLT2Yd3Y74esa4y5H7gfQERqRWRX53XSqByoS3cjDoO2P720/QPrRWPM+eluhOpeqoLRZsAhIlOMMVuiy44FNnSx7oboY0/FrXfAGFPf0wGMMRUD1diBICIrjDFz092O/tL2p5e2Xw01KRmmM8a0Ac8APxKRfBE5FbgIeKyL1R8FPici00VkGPB94OFUtFMppVR6pHJq95eAXKAG+CvwRWPMBhGpFBGPiFQCGGNeBH4GvA7sit5uS2E7lVJKpVjK6hkZYxqAi7tYvhtr0kL8sl8Av0hNy5Lm/nQ34DBp+9NL26+GlEFTz0gppVT20tx0Siml0k6DkVJKqbTTYDSARCRHRB6MJoJtFZHVInJBp3XOFpFNIuIVkddFZFy62tsVEflK9ELidhF5uIvHM739fU7Imwl6er2z4LXu8f2e6e1XmUWD0cByAHuAM4Bi4AfAUyIyHkBEyrGmuP8AKAVWAH9LS0u7tx+4E3io8wNZ0v4+JeTNIF2+3lnyWnf7fs+S9qsMohMYkkxE1gJ3GGP+ISI3AtcZY06JPpaPdZX6ccaYTelsZ2cicidwhDHmurhlGd3+aHsagaNjeRBF5DFgnzHmO2ltXC86v96Z/lp3J/Z+B8rIwvar9NGeURKJyAisJLGxTBMHJYGNXgy8ja4TxmaiTG9/Igl5M12mv9aH6PR+z7r2q/TSYJQkIuIEHgceifsmmEjC2EyU6e3P9PYlIqueSxfv96xqv0o/DUYJEJHFImK6uS2NW8+GleooAHwlbheJJIwdcH1tfw/S2v4+yPT2JSJrnks37/esab/KDBqMEmCMOdMYI93cTgMQEQEexDqBfqkxJhi3i1gSWKLr5gOT6DphbFra34u0tr8POhLyxi3rLiFvpsv01xro8f2eFe1XmUOD0cC7DzgK+IQxxtfpsWeBo0XkUhFxAz8E1mbSCV0RcUTbZgfsIuIWkVjaqIxuf4IJeTNCD693Rr/Wcbp7v2dL+1WmMMbobYBuwDjAAH6sYYrY7cq4dc4BNgE+YDEwPt3t7vQcbo8+h/jb7VnU/lLgn0AbsBu4It1t6u/rnQWvdY/v90xvv94y66ZTu5VSSqWdDtMppZRKOw1GSiml0k6DkVJKqbTTYKSUUirtNBgppZRKOw1GSiml0k6DkUorEdkgImem6FjTRWRFEvb7jIicP9D7VWoo0euMVFKJiCfubh7QDoSj928yxjyewrb8A3jaGPPkAO/3BOA+Y8ycgdyvUkOJBiOVMiKyE/i8MWZRGo49Cisv2mhjjD8J+98CfNYYM+A9L6WGAh2mU2klIjtF5Jzo77eLyNMi8pdoGet1IjJVRP6fiNSIyB4ROTdu2+Jo2esqEdknIneKiL2bQy0AVsUHouixbxGRtdEy5Q+KyAgRWRg9/iIRGRZd1x1tV72INInI8mj9npjFwIUD/gIpNURoMFKZ5hNYiU2HAauBl7Dep2OAHwF/jFv3ESAETAaOA84FPt/Nfo8BPuxi+aVYgWpq9NgLge8C5dHjfi263rVYpbXHYlUx/QJWzrWYD4jLUq2USowGI5Vp3jTGvGSMCQFPAxXA3cYqTfAkMF5ESqK9kguArxtj2owxNcAvgc90s98Suq6l83/GmAPGmH3Am8B7xpjVxph2rMzTx0XXC2IFocnGmLAxZqUxpiVuP63RYyil+sHR+ypKpdSBuN99QJ0xJhx3H6wqoqMBJ1BlldQBrC9Xe7rZbyNdVxntfLzO9wuivz+G1St6UkRKgL8A3zMf1e8pBJq6e1JKqZ5pz0hlqz1YM/PKjTEl0VuRMWZGN+uvxRqK6xdjTNAYc4cxZjpwCvBx4Jq4VY4C1vR3/0oNdRqMVFYyxlQBLwP3ikiRiNhEZJKInNHNJq8As6OF3hImIvNF5JjoBIkWrGG7cNwqZ2Cdb1JK9YMGI5XNrgFcwEasYbi/A6O6WtEYcwB4Davya3+MjO6/BWuywhtYQ3WIyPFAmzFmWT/3rdSQp9cZqSFDRKZjzcA7wQzgGz96Me2Dxpj/DNQ+lRpqNBgppZRKOx2mU0oplXYajJRSSqWdBiOllFJpp8FIKaVU2mkwUkoplXYajJRSSqWdBiOllFJpp8FIKaVU2v1/Zv91x63eVjoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x324 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set maxlag\n",
    "maxlags = int(0.025*Fs)\n",
    "\n",
    "fig, ax = plt.subplots(figsize =(6, 4.5))\n",
    "\n",
    "# colors for plotting \n",
    "colors = [cm.Purples(i) for i in np.linspace(0, 1, 25)][5::9]\n",
    "\n",
    "for n, w in enumerate(np.arange(0, 1.1, 0.4)):\n",
    "    # adjust the connectivity strength\n",
    "    weights_n = np.array([[0, 0], [w, 0]])\n",
    "    \n",
    "    # simulate the network\n",
    "    simulation_n = utils.simulate(num_units=num_units, num_trials=num_trials,\n",
    "                                      duration=duration, step_size= step_size, weights = weights_n, inp_curr=input_current)\n",
    "    \n",
    "    \n",
    "    # calculate spike rates\n",
    "    rates_n = np.array([[calculate_rate(T, win = rate_w, spike_times = spike_times) for spike_times in row.spike_times] for _, row in simulation_n.iterrows()])\n",
    "    \n",
    "    # calculate the trial averaged cross correlation \n",
    "    lags_n, mean_corr_n, sem_corr_n = avg_cross_correlation(rates_n[1], rates_n[0], maxlags=maxlags)\n",
    "    \n",
    "    # plot it\n",
    "    ax.errorbar(lags_n, mean_corr_n, yerr = sem_corr_n, linestyle = '-',\n",
    "                color = colors[n], label = 'w = ' + str(np.round(w, 2)));\n",
    "\n",
    "# plot clean up  \n",
    "ax.legend(bbox_to_anchor = (1.01, 1))\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.set_title('Avg. cross-correlation\\nbetween neuron 2 and neuron 1')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('correlation')\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88dd9582",
   "metadata": {},
   "source": [
    "So the simulations are working! \n",
    "\n",
    "Important: Cross-correlation serves as an excellent initial step to assess functional connectivity, particularly for capturing linear relationships between neural signals. However, interactions among single neurons in the brain are often non-linear. Consequently, relying solely on cross-correlation analyses might overlook significant interactions in such cases. To capture non-linear interactions between units, <b> directed information </b> serves as a better metric. \n",
    "\n",
    "Now, let's see how to calculate directed information between these units."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b6147c3",
   "metadata": {},
   "source": [
    "## Directed Information"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d74f306",
   "metadata": {},
   "source": [
    "To calculate directed information on point process data[4], we need to consider several key components:\n",
    "\n",
    "1. <b>Rate function</b>:  Determine the intensity function, representing the spike rate. Recommended types of intensity functions: boxcar or Gaussian filter. \n",
    "<br> </br>\n",
    "2. <b>Rate function filter width</b>: The width of the filter used is crucial in determining the temporal resolution at which interactions can be detected. \n",
    "<br></br>\n",
    "3. <b>Memory</b>: Define the length of past history to consider when evaluating directed information. This parameter can be determined based on experimental requirements, the window of interest, or through empirical testing with various values.\n",
    "<br></br>\n",
    "4. <b>Kernel width of the Gram matrix</b>: This parameter influences how similar closely related data points are considered when assessing temporal dependencies and interactions within the point process data.\n",
    "<br></br>\n",
    "Advanced hyperparameters include: \n",
    "1. Kernel to evaluate the Gram matrix\n",
    "2. Order of alpha entropy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c42f125d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from calculate_DI import directed_information"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b45dbd02",
   "metadata": {},
   "source": [
    "Since we modeled our data such that neuron 1 $(n_1)$ causally influences neuron 2 $(n_2)$, we should see that \n",
    "\n",
    "$$ DI (n_1 \\rightarrow n_2) > DI (n_2 \\rightarrow n_1)$$   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "057df112",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 4min 51s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Initialize\n",
    "di = {(n1, n2): np.zeros(num_trials) + np.nan for n1 in range(num_units) for n2 in range(num_units) if n1 != n2}\n",
    "memory = 15\n",
    "\n",
    "# Evaluate DI for each trial and in both directions\n",
    "for trial in range(num_trials):\n",
    "    for n1 in range(num_units):\n",
    "        for n2 in range(num_units):\n",
    "            if n1 != n2:\n",
    "                model = directed_information(memory=memory)\n",
    "                di[(n1, n2)][trial] = model.DI(rates[n1, trial], rates[n2, trial])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c457683b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEUCAYAAABnHY4tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgtklEQVR4nO3de7hVdb3v8feHi4oCooDsUJdoFxM18JEsb1tKfVRU6qT78a54LLzUycums7EtSmianh6z9i6NDopQmWYWXjna0dSMw9YyQHZKbRFBAREVgbwAfs8fYyydTOZcjLXWnL81mX5ezzMfxhjzN+f4Mlh+/a3f+I3fVxGBmZnVX7euDsDM7MPCCdfMLBEnXDOzRJxwzcwSccI1M0vECdfMLJEeXR1ArRx99NExc+bMrg7DzEzV3miaHu6rr77a1SGYmbWpaRKumVmjc8I1M0vECdfMLBEnXDOzRJpmloJZM1i3bh1Llizh7bff7upQrA3bbLMNu+yyCz179mzX55InXEkfB+YBd0bE6VXaXAz8C9AL+BVwfkS8ky5Ks66xZMkS+vTpw5AhQ5Cqzi6yLhQRrFy5kiVLlrD77ru367NdMaTwQ+DJam9KOgoYDxwODAH2AL6VJDKzLvb222/Tv39/J9sGJon+/ft36LeQpAlX0snAG8D/baPZWcCUiJgfEa8DVwJj6h+dWWNwsm18Hf03SjakIKkvMIms53pOG033BmaU7M8BBknqHxEry75zLDAWoKWlpbYBf8g9+uijSc932GGHJT2fWVdI2cO9kqznungz7XoDq0r2W7f7lDeMiMkRMSIiRgwcOLBGYZqZ1UeShCtpOHAE8L0CzdcAfUv2W7dX1zgsM0vg0ksv5YYbbujqMNrlgAMOYP78+TX/3lQ93JFkN8BelLQMGAecIOlPFdrOB4aV7A8DlpcPJ5hZekOGDKFXr1706dOHfv36cdBBB3HTTTfx3nvvbdTmt7/9LQArVqxg2rRpnHvuuV0VcoeMGzeOyy+/vObfmyrhTgY+CgzPXzcB9wFHVWg7DThH0lBJOwCXAVOTRGnWYCTV9dUR99xzD6tXr2bRokWMHz+ea6+9lnPOqXxbZurUqYwaNYpevXp15jIkN3r0aB555BGWLl1a0+9NknAj4u8Rsaz1RTZs8HZErJDUImmNpJa87UzgOuARYFH+uiJFnGZW3Pbbb8/o0aO5/fbbufXWW3nmmWc2afPAAw9sckP0pptu4thjj+WrX/0qAwYMYPDgwTz00EOdjqdaBfKOnG+bbbZh//3358EHH+x0XKW65NHeiJjY+tBDRLwYEb0j4sWS96+PiEER0TcizvZDD2aN64ADDmCXXXbh8ccf3+S9efPmseeee250bO7cucyaNYvRo0fzyiuvcO6553Lttdd2KobHH3+cUaNGVZwb29Hz7bXXXsyZM6dTcZXzWgpm1mmDBw/mtdde2+T4G2+8QZ8+G08wmjt3LuPHj+eoo46iW7duDB069P33Vq1axQEHHEDv3r0r9pirOfjgg9lpp50YPXr0Jkm3rfPNmjWLAw88kMMOO4xTTjmFdevWvf9enz59eOONNwrHUIQTrpl12ksvvcSOO+64yfEddtiB1as3nmA0b948jj/++Pf3n3nmmfeT4Lbbbst9993HiSeeWPE8M2fOrDgW3b17d6ZNm8ZDDz3EjTfeWPh8u+22Gw8//DCPPvooe+yxBzNmfPAIwOrVq+nXr1/7LsRmOOGaWac8+eSTvPTSSxxyyCGbvPepT32KBQsWvL+/cOFC1q9fv9Eww9NPP83w4cMB6NmzJ23NqT/66KOJiE1eGzZs4Mwzz+TII4/kvPPOK3y+wYMHv39Dr0ePHnTr9kFK/Mtf/sKwYaUTpjrPCdfMOuTNN9/k3nvv5eSTT+b0009n33333aTNqFGjNnpqce7cuey7774bJbann36604ntiSeeYPny5cyYMWOjGRFFz7dw4UIeeOABjjvuOADeeecd/vjHP3LkkUd2Kq5yXp7RzNrl+OOPf783OHToUC655JKNepWlzjzzTIYPH85bb71Fr169mDt37vu9S4CVK1eybNky9tlnn07FdOihh/LAAw9sMtWtyPnefPNNzjrrLKZPn85WW20FwN13383IkSMZPHhwp+Iq54RrZoW98MIL7WozYMAAzjzzTH784x9z0UUXMWHChI3a9u/ff6MbVZ1RaV7x5s63fv16TjnlFCZOnLjRsMN3v/tdpkyZUpO4SnlIwayBVRqvrOUrhauvvpqLLrqocPtRo0bx4IMP8pWvfIWpU6fWLS6A2267jdmzZzNp0iRGjhzJ7bffDsDs2bM73euuxD1cM2so999/f7JznXHGGZxxxhnJzucerplZIk64ZmaJOOGamSXihGtmlogTrplZIk64ZmaJOOGamSXihGtmlogTrplZIskSrqSfSloq6U1JCyR9uUq7MZI25GV3Wl8jU8VpZrXlqr0fSNnDvQYYEhF9gdHAVZL2r9J2Vl52p/X1u2RRmllVrtrbOcnWUoiI0v9dRP76KPDHVDGYbWlK15Kth/ICj0Xcc889HHHEEaxatYpHH32UCy+8kNmzZ3PLLbds0nZLrtp73nnnsXTpUj7ykY/U7HuTjuFK+pGkvwPPAkuBaqtU7Cfp1XzoYYKkiv9jkDRW0lOSnlqxYkW9wjazCly1t/2SJtyIuADoAxwK3AVUqsb7GLAPsBNwAnAK8I0q3zc5IkZExIi2ynKYWf24am9xyWcpRMSGiPg9sAtwfoX3n4+IhRHxXkTMAyYBlSvKmVlDqFXV3raq6Lalo1V726oS3GxVe3uQjeFuTgCbLuVuZg2jVlV726qiC7Wv2ttWleAttmqvpJ0knSypt6Tuko4iGyp4uELbYyQNyrc/CUwAZpS3M7PGUMuqvW1V0YXaV+1tq0rwlly1N8iGD5YArwPfBS6KiBmSWvK5ti1528OBuZLWkt1Uuwu4OlGcZlZQPav2llfR3ZzOVu0tt0VX7Y2IFUDF+ScR8SLQu2R/HDAuRVxm1n71rtpbqYru5nSmam8lrtprZl2u3lV7q1XRLaIjVXurqVfVXidcswbWkQcTGs3VVxcfESytojtp0iTOP/98TjrppDpGlw17/PnPf+a5557j3HPPZcyYMcyePbsu53LCNbOGkbqKLqStEuzVwszMEnHCNTNLxAnXzCwRJ1wzs0SccM0aTLVVr6xxdPTfyAnXrIF079698IIt1nXWrVtHjx7tn+TlhGvWQPr168fy5cs3qqBgjeW9995j+fLlbL/99u3+rOfhmjWQAQMGsGTJEp577rmuDsXasN122zFgwIB2f84J16yBdOvWjZaWls03tC2ShxTMzBJxwjUzS8QJ18wsESdcM7NEnHDNzBJJlnAl/VTSUklvSlog6ctttL1Y0jJJqyTdLGnrVHGamdVLyh7uNcCQiOgLjAaukrR/eaO8wOR4stpmQ4A9gG8ljNPMrC6SJdyImB8R77Tu5q9KZdLPAqbk7V8HrgTGpInSzKx+kj74IOlHZMmzF/A0WVXecnuzcVn0OcAgSf0jYmXZ940FxgJNMVm8tLqpmTWfpDfNIuICoA9wKFn583cqNOsNrCrZb93uU+H7JkfEiIgYUa22vJlZoyiUcCVdIml4vv1ZSS9Kel7Sge09YURsiIjfA7sA51dosgboW7Lfur26vecyM2skRXu4FwML8+1rgOuBbwM3dOLcPag8hjsfGFayPwxYXj6cYGa2pSmacLePiFWS+pAlwH+LiClAocLxknaSdLKk3pK65zMRTgEertB8GnCOpKGSdgAuA6YWjNPMrGEVTbiLJR0EnAw8FhEbJPUFNhT8fJANHywBXge+C1wUETMktUhaI6kFICJmAtcBjwCL8tcVhf9GZmYNqugshW8AdwLvAifkx44D/qPIhyNiBXBYlfdeJLtRVnrserJhCzOzplEo4UbE/cDgssO/zF9mZlZA1YQraY+C3/F8jWIxM2tqbfVw/0Y29qo22gTQvaYRmZk1qaoJNyK8kpiZWQ05qZqZJVLoppmkHsAFZDMNBlAyzBAR/1if0MzMmkvRHu73gHOBx4D9gV8BO1H5wQUzM6ugaML9EnBMRHwfWJ//+UXgc/UKzMys2RRNuNsCi/PttyRtGxHPAvvVJywzs+ZT9EmzvwCfJnuy7ClgoqQ3gZfqFZiZWbMpmnAv5IN1Ey4BbiRbn3ZsPYIyM2tGRR/tfbJk+6/AEXWLyMysSbX1aO8/RsRj+fbnq7WLCM9UMDMroK0e7o+AffLtKVXaBFlVXTMz24y2Hu3dp2T3YxFRdO1bMzOrYLPTwiR1B9ZI2jpBPGZmTWuzCTfv2S4A+nf0JJK2ljRF0iJJqyU9LemYKm3HSNqQV4FofY3s6LnNzBpF0WlhPwPulfR9sjI50fpGwZtmPcgenDgMeBEYBdwhad+IeKFC+1kRcUjB2MzMtghFE25rOfOJZccL3TSLiLVln71X0kKydRleKBiDmdkWreg83N1reVJJg4BPkJVEr2Q/Sa8CrwHTgWsiYn0tYzAzS63QWgqSZlQ5fld7TyipJ9kQxa35egzlHiObjrYTWcHKU8iKWFb6rrGSnpL01IoVK9obiplZUkUXr6m2KtjI9pxMUjeyHuu7wNcqtYmI5yNiYUS8FxHzgEnAiVXaTo6IERExYuDAge0JxcwsuTaHFCRNyje3KtlutQewqOiJJInsAYpBwKiIWFfwo5urq2ZmtkXY3Bjurvmf3Uq2IUuCi9n0JlpbbgT2Ao6IiLeqNcqni/0pIpZL+iQwAZdjN7Mm0GbCjYizAST9ISJ+0tGTSNqNrGLEO8CyrLML+bHHgf8EhkbEi8DhwFRJvYHlwE+Bqzt6bjOzRlF0lkKHk23++UW0PSzQu6TtOGBcZ85nZtaIXLXXzCwRJ1wzs0SqJlxJ/6tku+p6uGZmVkxbPdzS8jm/qXMcZmZNr62bZnMk3Uk2g2DrCvNwAYiIy+sSmZlZk2kr4Z5I1svdjWyGwa4V2kSFY2ZmVkFbFR9eAa4CkNSjdU6umZl1TNF5uGdL2gE4HtgZeAm4NyJeq2dwZmbNpOhqYQcC/wWcB3yK7Amxv+XHzcysgKILkN8AXBARv2g9IOkk4AfAp+sQl5lZ0yn64MMngDvKjt0JfKy24ZiZNa+iCfevwMllx/6JbJjBzMwKKDqkcBFZHbKvk62BOwT4OHBcfcIyM2s+RWcp/EHSR4FjgcHAPcD9nqVgZlZc0R4uEfE62dq0ZmbWAV4tzMwsESdcM7NEij740KnELGlrSVMkLZK0WtLTee2yau0vlrRM0ipJN0vaujPnNzNrBJtNpJK6A2s7mfR6kBWdPAzYnqww5B2ShlQ431HAeLLaZkPIqgN/qxPnNjNrCJtNuBGxAVgA9O/oSSJibURMjIgXIuK9iLgXWAjsX6H5WcCUiJif36i7EhjT0XObmTWKorMUfkY2D/f7wBJKlmWMiIfbe1JJg8ieXptf4e29gRkl+3OAQZL6R8TK9p7LzKxRFE245+d/Tiw7HmS/8hcmqSdZAr81Ip6t0KQ3sKpkv3W7D7BRwpU0lrwyRUtLS3vCMLOc1FZB7fqIqL6Udup42oql1oo++LB7LU6W33ybDrwLfK1KszVA35L91u3VFeKaDEwGGDFihBdDN7OGVnj2gaSekg7NVwlD0naStmvH5wVMAQYBJ0TEuipN5wPDSvaHAcs9nGBmW7qi08L2Jbtx9hOypAnZjIOb23GuG4G9gOMj4q022k0DzpE0NF/0/DJgajvOY2bWkIr2cG8ELo+ITwKtPdNHgUOKfFjSbmSLlg8Hlklak79Ok9SSb7cARMRM4DrgEbKFchYBVxT9C5mZNaqiN8325oN1FAKyqV6SehX5cEQsIitEWU3vsvbXA9cXjM3MbItQtIf7AmVzZiUdAPyt1gGZmTWroj3cCcB9km4CtpJ0KVl9s6/ULTIzsyZTqIebPxl2DDCQbOx2N+BLEfFgHWMzM2sq7VkP90/ABXWMxcysqRWdFraVpEmS/ippbf7nlZK2qXeAZmbNomgP90ZgT6C1ptluwKXAzsB/r09oZmbNpWjC/SLw0Yh4I9//T0mzyWYpOOGamRVQdFrYMmDbsmO9gKW1DcfMrHlV7eFK+nzJ7nRgpqR/I1uecVfgq2SP4ZqZWQFtDSlMqXDsm2X75wLX1i4cM7PmVTXh1mpJRjMzy7hqr5lZIkXn4Q6T9LCk1yS9m7/WSXq33gGamTWLotPCbgN+RTYPt621bM3MrIqiCfcfyNbDdRkbM7MOKjqGeytwaj0DMTNrdkV7uN8BZkn6JrC89I2I+Hzlj5iZWamiCfdOYCHwazo4hivpa8AYYF/gtogYU6XdGLI5wKXnOS4ifteR85qZNYqiCXc40D8iOjMr4WXgKuAosseC2zIrIgrVSzMz21IUHcN9HBjamRNFxF0R8RvA5c7N7EOpaA93IfCgpF+z6Rju5TWPCvaT9CrwGtk6DtdExPryRpLGAmMBWlpa6hCGmVntFE242wL3AVuRLVxTT48B+5Ctu7s3cDuwHrimvGFETAYmA4wYMcJT1sysoRVKuBFxdr0DKTnX8yW78yRNAr5BhYRrZrYlKZRwJe1R7b2yBFkPAajO5zAzq7uiQwp/Y9PE1/orfPciXyCpR36+7kD3vB7a+vKxWUnHAH+KiOWSPklWov2XBeM0M2tYRcukd4uI7vmf3YDBZGOnZ7TjXJeRza0dD5yeb18mqUXSGkmtd70OB+ZKWgvcD9wFXN2O85iZNaTCZdJLRcQySRcBC4CfF/zMRGBilbd7l7QbB4zrSFxmZo2sM+vh7smmdc7MzKyKojfNHueDMVvIEu3ewKR6BGVm1oyKDin877L9tcCciPhrjeMxM2taRefh3lrvQMzMml3RIYWtyFb6Gk7JDS6AiDiz5lGZmTWhokMKtwLDgHsoW0vBzMyKKZpwjwZ2j4g36hiLmVlTKzot7EVg63oGYmbW7Ir2cKcBMyR9n02XZ3y45lGZmTWhogn3a/mf5Y/YBlB1YRszM/tA0Wlhu9c7EDOzZteZR3vNzKwdnHDNzBJxwjUzS8QJ18wsESdcM7NEnHDNzBJJlnAlfU3SU5LekTR1M20vlrRM0ipJN0vyU25mtsVL2cN9GbgKuLmtRpKOIqt7djgwhOzBim/VOzgzs3pLlnAj4q6I+A2wcjNNzwKmRMT8iHgduJJsaUgzsy1ah4pI1tnewIyS/TnAIEn9I2KjZC1pLDAWoKWlBbO2SEp6voio+l7qWBpNI/39axVLW//erRrxpllvYFXJfut2n/KGETE5IkZExIiBAwcmCc7MrKMaMeGuAfqW7Ldur+6CWMzMaqYRE+58suoSrYYBy8uHE8zMtjQpp4X1kLQN0B3oLmkbSZXGkKcB50gaKmkH4DJgaqo4zczqJWUP9zLgLbIpX6fn25dJapG0RlILQETMBK4DHgEW5a8rEsZpZlYXyWYpRMREYGKVt8srAV8PXF/nkMzMkmrEMVwzs6bkhGtmlogTrplZIk64ZmaJOOGamSXihGtmlogTrplZIk64ZmaJOOGamSXihGtmlogTrplZIk64ZmaJOOGamSXihGtmlogTrplZIk64ZmaJpCyxs6OkX0taK2mRpFOrtBsjaUNeBaL1NTJVnGZm9ZKs4gPwQ+BdYBAwHLhP0pyImF+h7ayIOCRhbGZmdZekhytpO+AEYEJErImI3wN3A2ekOL+ZWSNINaTwCWBDRCwoOTYH2LtK+/0kvSppgaQJVar7mpltUVIlst7AqrJjq4A+Fdo+BuxDVq13b+B2YD1wTXlDSWOBsQAtLS01DNfMrPZS9XDXAH3LjvUFVpc3jIjnI2JhRLwXEfOAScCJlb40IiZHxIiIGDFw4MCaB21mVkupEu4CoIekj5ccGwZUumFWLgDVJSozs4SSJNyIWAvcBUyStJ2kg4EvANPL20o6RtKgfPuTwARgRoo4zczqKeWDDxcAvYBXgNuA8yNivqSWfK5t6yDs4cBcSWuB+8kS9dUJ4zQzq4tkd/8j4jXgixWOv0h2U611fxwwLlVcZmap+NFeM7NEnHDNzBJxwjUzS8QJ18wsESdcM7NEnHDNzBJxwjUzS8QJ18wsESdcM7NEnHDNzBJxwjUzS8QJ18wsESdcM7NEnHDNzBJxwjUzS8QJ18wsESdcM7NEkiVcSTtK+rWktZIWSTq1jbYXS1omaZWkmyVtnSpOM7N6SdnD/SHwLjAIOA24UdLe5Y0kHQWMJ6ttNgTYA/hWujDNzOojScKVtB1wAjAhItZExO+Bu4EzKjQ/C5gSEfMj4nXgSmBMijjNzOopVQ/3E8CGiFhQcmwOsEkPNz82p6zdIEn96xifmVndpara2xtYVXZsFdCnQNvW7T7AytKGksYCY/PdNZKe63yoXWoA8GpXB9GAtsjrIinFabbIa5NA8utS8u89MyKOrtQmVcJdA/QtO9YXWF2gbev2Jm0jYjIwuRYBNgJJT0XEiK6Oo9H4ulTna1NZo16XVEMKC4Aekj5ecmwYML9C2/n5e6XtlkfEygptzcy2GEkSbkSsBe4CJknaTtLBwBeA6RWaTwPOkTRU0g7AZcDUFHGamdVTymlhFwC9gFeA24DzI2K+pBZJayS1AETETOA64BFgUf66ImGcXalphkdqzNelOl+byhryuigiujoGM7MPBT/aa2aWiBOumVkiTrgJFV1PQtI+kv6PpFclfSjGfNpxbc6S9EdJb0paIuk6SammNybXjutysqTn8vVHXpF0q6TyqZhNpT3rs5R85mFJ0VU/M064aRVaTwJYB9wBnJMwtq5W9NpsC1xENrH9M2RrboxLFGNXKHpdngAOjojtydYf6QFclSzKrlH02gAg6TTSPXtQOQbfNEsjX0/idWCf1kecJU0HXoqI8VU+8zHgrxGR5JGlrtKRa1Py2UuAz0XE8fWPNK2OXhdJvYEfAQMiYlSSYBNr77WRtD3wJHAmMAvoGRHrE4YMuIebUnvWk/iw6cy1+UcqP0DTDNp1XSQdImkV2VOZJwA31D3CrtPen5mrgRuBZfUOrC1NO/bVgNqznsSHTYeujaSzgRHAl+sUV1dr13XJV+HbXtLOwFeAF+oaXdcqfG0kjQAOBi4Edql/aNW5h5tOe9aT+LBp97WR9EXgO8AxEdGsi7d06GcmIl4CZgK/qFNcjaDQtZHUjWx45cKuGEIo54SbTnvWk/iwade1kXQ08BPg+IiYlyC+rtKZn5kewEfrElVjKHpt+pL9FnS7pGVk47gASyQdWv8wN+abZglJ+gUQZL8CDwfuBw6KiPll7QRsTXa3eT7ZI9EREe8kDTihdlybzwO/BP5bRDyWOs7U2nFdTgMeBxYDLWRrkqyMiC8lDTihItcm/29pUMnHdgX+g2xoYUVEvJssYICI8CvRC9gR+A2wFngRODU/3kL2K1JLvj+E7Aep9PVCV8ffINfmEWB9fqz19UBXx98A1+XbwJK83RKytQT6d3X8jXBtyj7T+t9Wj66I2T1cM7NEPIZrZpaIE66ZWSJOuGZmiTjhmpkl4oRrZpaIE66ZWSJOuFaYpKmSrsq3D5X0XFfHVImkFyQdUeW9myRNKPg9vSTdk68x+8vaRtl5kk6T9GBXx2HFefEa65CIeBzYs9bfK2kk8NOIqMsiIxFxXjuan0j2lFL/6OLn8CUNARZSsqxgRPwM+FlXxmXt4x6u1VwTVWDYDVjQkWTbRNfAasgJ16qStJ+kP0laLel2YJuS90ZKWlKy/4Kkf5E0F1grqYekz0r6g6Q3JM3Je6+t7XeUdIuklyW9Luk3+aLSDwCDJa3JX4MldZM0XtJ/SVop6Q5JO5Z81xl5iZWVkv51M3+n0mGRkXmZnn/Oy9IszZd8RNK3gMuBk/I4zsnjuCw/1yuSpuULWyNpSF665RxJLwIPSxoj6QlJ38uvwfOSDsqPL86/46yS2I6V9LSy8kGLJU0sCb113Yg38ngOzL/n9yWfP0jSk/kQyJOSDip573eSrszjWS3pQUkDNvczYLXlhGsVSdqK7Dn16WTPrP+SbFHrtpwCHAv0I/tV/D6yMi87kpXB+ZWkgXnb6WTlcvYGdgK+FxFrgWOAlyOid/56Gfg68EXgMGAw2Ur/P8zjHEq2sPQZ+Xv9ad+ap/8AbA/sTFbS6IeSdoiIK8gWrb49j2MKMCZ/fY5sYaHewL+Xfd9hwF7AUfn+Z4C5eVw/J1sy8dPAx4DTgX9XVqEBsjUBzsyv37HA+cqWoYRsoXWAfnk8s0pPmv8P6D7gB/m5rgfuk9S/pNmpwNlk13srmrs0UUNywrVqPgv0BG6IiHURcScfLG1XzQ8iYnFEvEWWTO6PiPsj4r2IeAh4Chgl6SNkifW8iHg9//5H2/jec4F/jYglka2YNhE4Mf+1/UTg3oh4LH9vAvBeO/6e64BJeQz3ky16Um1s+jTg+oh4PiLWAJcCJ5cNH0yMiLX5NQBYGBG3RMQG4Hay1aomRcQ7EfEgWU2ujwFExO8iYl5+veYCt5El8CKOJSvHND0i1kfEbcCzQGnpoVsiYkEe2x1kK2xZQh5nsmoGk9WHKl3daNFmPrO4ZHs34J8klf4H35Nsta9dgdci4vWCsewG/FpSaSLdQNaLHlx63ohYK2llwe+FbAnD0jHav5P1XCsZzMbXYBHZf0Oly/8tZmPLS7bfymMsP9YbQNJnyBZV34esB7o12W8WRZTH1hrfziX7peVl2vp7Wp24h2vVLAV2ztcTbdWymc+UJufFwPSI6Ffy2i4ivpO/t6Okfpv5jtLvOqbsu7aJrLLBUrIEDoCkbcl+pa6Hl8mSf6sWsqUiSxNoZ5bf+zlwN7BrZNV3bwJar//mvrc8ttb4XupEPFZjTrhWzSyyZPL1/AbYl4AD2vH5nwLHSzpKUndJ2+Q3qXaJiKVkN8d+JGkHST0ltY5RLgf6t96Myt0EfFvSbgCSBkr6Qv7encBxygoobgVMon4/17cBF0vaPR93bR3jrdWUsT5kPf+3JR1ANubaagXZUMkeVT57P/AJSafm/14nAUOBe2sUm9WAE65VFNlK+F8iu0n0OnAScFc7Pr8Y+ALwTbJksRj4Bh/8zJ1BNn76LPAKcFH+uWfJEtvz+Z39wcD3yXp+D0paDfw/sptRRLa6/1fJeodL81jfnz1RYzeT3ex7jGxO7NvA/6jh918ATMr/jpeTjbMCEBF/J1tk/In8uny29IMRsRI4DvhnYCXwP4HjonnrvW2RvAC5mVki7uGamSXihGtmlogTrplZIk64ZmaJOOGamSXihGtmlogTrplZIk64ZmaJOOGamSXy/wFwrDSEmun7tAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the histogram of DI in either directions for comparison \n",
    "fig,ax = plt.subplots(figsize = (5, 4))\n",
    "lim = np.max(np.concatenate((di[(0, 1)], di[(1, 0)])))\n",
    "ax.hist(di[(0, 1)], color = 'k', bins = np.arange(0.05, lim, 0.025), label = r\"DI($n_1\\rightarrow n_2$)\", width = 0.04)\n",
    "ax.hist(di[(1, 0)], color = 'silver', bins = np.arange(0.05, lim, 0.025), label = r\"DI($n_2\\rightarrow n_1$)\", width = 0.04)\n",
    "\n",
    "# plot clean up\n",
    "ax.set_xlabel('directed information')\n",
    "ax.set_ylabel('number of trials')\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "ax.legend()\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b91f8a0",
   "metadata": {},
   "source": [
    "### Hyperparameter tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e4e75aa",
   "metadata": {},
   "source": [
    "How do memory and Kernel width of the Gram matrix influence the calculation of directed information?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0d9b02a",
   "metadata": {},
   "source": [
    "#### Memory size "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "95260b44",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2min 31s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "trial = 8 # picking a random trial\n",
    "memory_range = [2] + list(np.arange(5, 30, 5)) # memory range to sweep over\n",
    "\n",
    "# Initialize\n",
    "di_memory = {(n1, n2): np.zeros(len(memory_range)) + np.nan for n1 in range(num_units) for n2 in range(num_units) if n1 != n2}\n",
    "\n",
    "# loop through the memory values\n",
    "for n, memory in enumerate(memory_range):\n",
    "    for n1 in range(num_units):\n",
    "        for n2 in range(num_units):\n",
    "            if n1 != n2:\n",
    "                model = directed_information(memory=memory) \n",
    "                di_memory[(n1, n2)][n] = model.DI(rates[n1, trial], rates[n2, trial])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3e4095ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEUCAYAAABnHY4tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4VElEQVR4nO3deXgUVdb48e8hBMKSBBM2QQgwLiPqqD9RGUVRQWULOCyDCMiIGkXGXZBFEBF01HdG5hXHDRfQuOAIA8giyOqGiq+C4gIqJIgIJJCVLCQ5vz+6k+nsDenuSnefz/P0Q+d2ddUpmhxu37p1j6gqxhhj/K+B0wEYY0y4sIRrjDEBYgnXGGMCxBKuMcYEiCVcY4wJEEu4xhgTIA2dDsBX+vTpo6tWrXI6DGOMkepeCJkeblpamtMhGGNMjUIm4RpjTH1nCdcYYwLEEq4xxgSIJVxjjAmQkJmlUJ2SkhLS0tLIyMiguLjY6XBMNSIiImjRogUtW7akQQPrB5jQFPIJ95dffkFE6NSpE5GRkYhUO2PDOERVOXr0KPv37+eXX36hY8eOTodkwlxycjJTp04lNTWVjh07Mnv2bEaOHFnn/YZ8VyI3N5f27dvTqFEjS7b1lIjQqFEj2rdvT25urtPhmDCXnJxMUlISKSkpqCopKSkkJSWRnJxc532HfMIF7CtqkLDPydQHU6dO5ciRI+Xajhw5wtSpU+u8b/sXbowxHlJTU4+p/VhYwjXGGA/x8fFVtvvi2oIl3CA1efJk5syZ43QYx+SCCy5g+/btTodhTLU2bNjA4cOHKw1vNW3alNmzZ9d5/5Zw64FOnTrRpEkToqOjadGiBRdddBHPPvssJSUl5bZ5//33ATh48CALFizglltucSrk43Lfffcxffp0p8MwpkrffPMN11xzDaeeeirPPPMMCQkJiAgJCQk8//zzPpmlEPLTwoLFsmXL6N27N5mZmWzcuJE777yTTz/9lJdffrnStq+88gr9+vWjSZMmDkR6/AYOHMitt97Kvn37OPHEE50Ox5gyv/zyC3379qVp06asXLmShIQEkpKSfH4c6+HWIDk5mU6dOtGgQQM6derkk2khtYmNjWXgwIG89dZbzJ8/n2+++abSNitXrqRnz57l2p599ln69+/P+PHjadmyJe3atWPNmjV1jqe6qs7Hc7yoqCjOO+88Vq9eXee4jPGVzMxM+vXrR0ZGBitWrCAhIcFvx7KEWw1/zsXzxgUXXMBJJ53EBx98UOm1r7/+mtNOO61c27Zt2/jkk08YOHAgBw4c4JZbbuGxxx6rUwwffPAB/fr1Iz8/v9Jrx3u8008/na1bt9YpLmN8pbCwkMGDB/Pdd9+xaNEizjnnHL8eL+yGFO666y6++uqrWrfbvHkzBQUF5dqOHDnCjTfeyAsvvFDje8855xyfXNBq164dhw4dqtSekZFBdHR0ubZt27YxadIkrr76agC6du1alqwzMzO58sor+fbbb9m8eTNnnnmmV8e/+OKLad26NQMHDmTp0qVERUV5dbxPPvmEe+65h0aNGtGuXTsWLFhAZGQkANHR0ezbt+8Y/yaM8b2SkhJuuOEG1q1bx/z587nyyiv9fkzr4VajYrKtrd0f9u7dS1xcXKX2E044gezs7HJtX3/9NYmJiWU/f/PNN3Tt2hVwXWFdvnw5Q4cOrfI4q1atQkQqPSIiIliwYAFr1qzhmWee8fp4CQkJrFu3jo0bN9KlSxeWLFlStl12djYtWrQ4tr8IY/xg8uTJvP7668yePZvrr78+IMcMux6utz3PTp06kZKSUqk9ISGBDRs2+DaoKnz++efs3buXHj16VHrtD3/4Azt27OD8888HYNeuXRQVFZUbZvjyyy+55pprAIiMjKRVq1bVHqtPnz5VjtWW9gD27dvHrbfeWtZe2/HatWtX1t6wYcNyU2y+++47Ro0aVcvZG+Nfc+fO5fHHH+fWW29l8uTJATuu9XCrMXv2bJo2bVquzVdz8WqSlZXFu+++y7XXXsuoUaM466yzKm3Tr18/Nm7cWPbztm3bOOuss8olti+//JKzzz67TrF89NFH7N+/nyVLlpSbEeHt8Xbt2sXKlSsZMGAA4Pp28MUXXwTkq5sx1Vm0aBF33HEHAwcOZO7cuQFdYyXserjeKp1z548Vg6qSmJhY1hvs2rUr99xzT7lepafrr7+ec845h7y8PJo0acK2bdvKDfanp6fz22+/eT1WW51LLrmElStXVvoH6c3xsrKyGDNmDK+++iqNGjUCYOnSpVx22WXlesDGBNJHH33EyJEjufDCC3njjTeIiIgIbACqGhKP8847T6vy7bffVtke7CZPnqxPPvnkMb1nzJgx+vXXX/snIA9Hjx7Vfv366dq1a8u1X3DBBbUeP1Q/L+O87777Tk844QQ95ZRT9ODBg/48VLV5SrSaeZbBplu3brply5ZK7d999x2nn366AxHVL/369eOrr74iISGBW265hb/85S9+O9arr77K3XffXdbjHTduHMOHD/fqvfZ5GX/47bff6N69O3l5eXzyySd06dLFn4erdozChhTCxIoVKwJ2rNGjRzN69OiAHc+YmmRnZ9OvXz8OHjxYNnPGKZZwjTEh6+jRowwbNoxt27axdOlSunXr5mg8lnCNMSFJVbn55pt57733mDdvHv369XM6JJsWZowJTdOnT2f+/Pk8+OCD3HjjjU6HA1jCNcaEoOeee45Zs2Zx44038uCDDzodTpmAJVwRiRORxSKSKyIpInJdNduJiMwSkb0ikikiG0TkjEDFaYwJbsuWLeO2226jb9++PPPMM/WqeGwge7hPA4VAG2Ak8Ew1iXQYMBa4BIgDPgFeDVSQxpjg9emnnzJ8+HDOPfdcFi5cWLZoUn0RkIQrIs2AIcA0Vc1R1Q+BpUBVc4c6Ax+q6s+qWgy8BnQNRJzGmOC1c+dOBgwYwIknnsjy5ctp3ry50yFVEqge7qlAsaru8GjbClTVw30TOFlEThWRSGAMsCoAMRpjgtSBAwfo27cvqsrKlStp06aN0yFVKVDTwpoDmRXaMoHoKrbdB3wA/AAUA3uAK6raqYgkAUngm4qaxpjgk5uby4ABA9i7dy/r1q3j1FNPdTqkagWqh5sDxFRoiwGyq9j2QeB8oAMQBTwErBORphU3VNXnVbWbqnarafnBUGRVe42BoqIihg8fzhdffMGbb77JH//4R6dDqlGgEu4OoKGInOLRdjZQ1W/f2cBbqvqLqhap6ivACYTwOK5V7TXm2Kkqt912G8uXL2fu3LkMGjTI6ZBqFZCEq6q5wCJgpog0E5GLgUFUPfvgc2CYiLQRkQYiMhqIBH4MRKxOWbZsGdnZ2aSkpDBp0iQee+yxaidrB3PV3vXr11uJHeMTs2bN4oUXXmDy5MmMGzfO6XC8EshpYbcBTYADwBvAOFXdLiIdRSRHREoHYR/DdUHtKyADuBsYoqoZAYwVgP3797N582Y2btzI5s2b2b9/v9+PaVV7jandyy+/zPTp0xk9erTfiwL4UsASrqoeUtVrVLWZqnZU1dfd7amq2lxVU90/56vqeFU9UVVjVPX/qWrAZyns37+fHTt2lNUwKygoYMeOHQFJumBVe42pzqpVq7j55pvp3bs38+bNq1c3NtQm7Bav+fHHH8nJyal1u6ysrEq9u5KSEn744YdavxI3b96ck08+uU5xgu+q9tZURbcmx1u1t6YqwVa119TFF198wdChQznrrLN45513yqqJBAtbS6Ea1X2VDuSC7b6q2ltTFV3wfdXemqoEW9Vec7x27dpF//79iY+PZ/ny5cTEVJz4VP+FXQ/X257n5s2bqyyJ3rhx43L1vPzFl1V7a6qiC76v2ltTlWCr2muOR3p6On379qWgoID169cHbV086+FWo3PnzpUSU4MGDejcubNfj+vPqr0Vq+jWpq5Veyuyqr3meOTl5ZGYmMju3btZunRpUJdgCrserrdKbw3ctWsXBQUFNG7cmM6dO/vtlkF/V+2tqopubepStbcqVrXXHKvi4mKuu+46Nm/ezMKFC7nkkkucDqluaqowGUwPq9pbveqq6PpbxSrBVrXXHIuSkhIdP368AjpnzhynwzkW1eYpxxOlrx7hlnCPxYIFCzQ+Pl579uypPXv21DfffNPvx+zbt6+eeOKJ2r17d3355Ze9fp99XqbU3/72NwX03nvvdTqUY2Vl0k1wsM/LACQnJzNq1CiuvfZakpOTK11PqeeqnRgcVGdhjAl9a9eu5YYbbuCyyy7jlVdeCbZkW6PQORNjTNDbunUrf/rTnzjttNNYvHgxjRs3djokn7KEa4ypF1JTU+nXrx8xMTGsWLEiJG+QsWlhxhjHHT58mL59+5KTk8OHH35Ihw4dnA7JL8Ii4apqUC1wEa5C5QKuOTb5+flcc8017Ny5k/fee6/Km31CRcgPKURGRpKXl+d0GMYLeXl59a7KqvGvkpISxowZw6ZNm5g/fz6XX3650yH5lVc9XBGJA+4DzsFVn6yMql7q+7B8p3Xr1uzdu5f27dvTpEkT6+nWQ6pKXl4ee/furbfF/4x/3HfffSxcuJDHH3+cESNGOB2O33k7pPA60BhYCBzxXzi+V7qi0K+//srRo0cdjsZUJzIykjZt2gTlClDm+Dz55JM8+eST3H777dx3331OhxMQ3ibci4BWqlp5+awgEBMTY7/IxtQjCxcu5J577mHw4ME8+eSTYfPN09sx3G3ASf4MxBgTHjZt2sTo0aO5+OKLee2114iIiHA6pIDxtoe7DlglIi8Dv3m+oKov+TwqY0xI2r59O4MGDaJLly4sXbo06Aqh1pW3CfcS4Beg4kKmCljCNcbUau/evfTt25eoqChWrlxZZTWTUOdVwlXV0J6rYYzxq8zMTPr168fhw4fZtGkTnTp1cjokR3h944OInAAkAu2BvcAyVT3sr8CMMaGhsLCQIUOG8O2337J8+XLOPfdcp0NyjFcXzUTkj8BPwK3AH4BbgJ/c7cYYU6WSkhLGjh3L2rVrmTdvHldddZXTITnK2x7uHOA2VX2ztEFEhgP/C5zvh7iMMSFgypQpJCcnM2vWLMaMGeN0OI7zdlrYqbhuevD0b8C7ErjGmLDz9NNP89hjj3HLLbcwZcoUp8OpF7xNuDuBayu0DcM1zGCMMeUsXryY22+/ncTERObOnRs2NzbUxtshhbuAd0XkDiAF6AScAnhXb9sYEzY+/vhjrrvuOi644ALefPNNGjYMi0UJveLttLCPReR3QH+gHbAMWKGqh/wZnDEmuPzwww8kJiZy0kknsWzZMpo2bep0SPWK1//1uKeAvebHWIwxQey3336jT58+REREsGrVKlq1auV0SPVOtQlXRFapah/38w9w3VVWSX1fntEY43/Z2dn079+fAwcOsGHDBn73u985HVK9VFMPd4HH83n+DsQYE5yOHj3KsGHD2Lp1K0uWLOH8822maHWqTbiq+rrHj9+r6qcVtxGRC/wSlTEmKKgqSUlJvPfee7zwwgv079/f6ZDqNW+nha2ppn2VrwIxxgSfBx98kFdeeYXp06dz0003OR1OvVfjRTMRaQCI66mI+3mp3wFFfozNGFMPJScnM3XqVFJSUgDo2bMnM2bMcDaoIFFbD7cIKASaup8f9Xh8C/zLr9EZY+qV5ORkkpKSypItwOeff87rr79ew7tMKampNLWIJODq1W4EPGcjKHBQVetNOdxu3brpli1bnA7DmJDWqVOncsm2VEJCArt37w58QPVTtbfV1djDVdUUVd2tqgnu56WP1GNNtiISJyKLRSRXRFJE5Loatu0iIu+KSLaIpInI48dyLGOMf6Smph5TuynvWNbDHQj0BFrikcFV9Xovd/E0ruGJNrjKrS8Xka2qur3CcRrhukj3NDAcKMa1eI4xxkHLli2jum/EHTt2DHA0wcnb9XAfBJ5zbz8MSAeuBjK8fH8zYAgwTVVzVPVDYCkwuorN/wL8qqr/UNVcVc1X1W3eHMcY4x8LFy5k8ODBdO7cuVIdsqZNmzJ79myHIgsu3k4LGwtcqap3A4XuPxNxLWLjjVOBYlXd4dG2FTijim27A7tFZKV7OGGDiJxV1U5FJElEtojIloMHD3oZijHmWCxYsIARI0bQvXt3vvrqK1544QUSEhIQERISEnj++ecZOXKk02EGhRovmpVtJJKpqrHu5weA9qp61LO9lvdfArytqm092m4GRqrqZRW2XQ1cDgwE1gJ3AuOA36tqYXXHsItmxvjec889x6233kqvXr1YsmQJzZo1czqkYHB8F808/CQipb3Rb4BxIjIa8LamWQ4QU6EtBsiuYts84ENVXelOsP8DxAOne3ksY4wPzJkzh1tvvZX+/fvz7rvvWrL1AW8T7gO4kh7AZOAO4AngHi/fvwNoKCKneLSdDWyvYtttVLNQjjEmMB555BHuvvtuhgwZwqJFi4iKinI6pJDgVcJV1RWqusn9/FNVPVlV26rqIi/fnwssAmaKSDMRuRgYBLxaxeavAd1FpLeIROBa/DwN+M6bYxljjp+q8sADDzB16lRGjhzJm2++SaNGjZwOK2Qcy7SwprhqmDX3bFfVj73cxW3AS8ABXLMcxqnqdhHpiOuuta7u+b0/iMgo4FmgNfB/wMCaxm+NMXWnqtx77708+eST3HTTTTz77LNEREQ4HVZI8Srhisj1wFxc82g9b3hQwKsJeO7qENdU0Z5K5SS+CFeP2BgTACUlJYwfP55nn32W22+/nTlz5tCggbcjjsZb3vZwHweGqGp1q4YZY4JUUVERN954IwsWLOD+++/n0UcftaKPfuJtwi0ENvgxDmOMA44ePcqoUaNYuHAhM2fO5IEHHrBk60fefmeYBvxDRFr6MxhjTODk5+czZMgQFi5cyBNPPMG0adMs2fqZtwl3B64bEfaLSLH7USIixX6MzRjjJ0eOHGHQoEEsW7aMp59+mvvuu8/pkMKCt0MKr+KqcfYW5S+aGWOCTHZ2NgMGDOCDDz7gpZde4oYbbnA6pLDhbcKNB6arN/cBG2PqrYyMDPr06cOWLVtITk5mxIgRTocUVrwdUniZqlf2MsYEibS0NK644gr+7//+j7ffftuSrQO87eFeAPxVRKYC+z1fUNVLq36LMaa+2LdvH7179+bnn39m6dKl9OnTx+mQwpK3CfcF98MYE2T27NlDr169+PXXX1mxYgWXX3650yGFrVoTrns9gxuAq1W1wP8hGWN85eeff+aKK67g8OHDrF69mosuusjpkMJarQlXVYtFpDM1rPFojKl/vv/+e3r16kV+fj7r1q3jvPPOczqksOftRbOHgGdFJEFEIkSkQenDn8EZY47Ptm3b6NmzJ0VFRWzYsMGSbT3h7RjuPPefnjMVBNfiNbackDH1yJYtW7j66qtp0qQJa9eu5bTTTnM6JOPmbcLt7NcojDE+8dFHH9GvXz/i4uJYu3YtXbp0cTok48GrhKuqKQDuIYQ2wH5VLfFnYMaYY7Nu3ToSExNp3749a9eupUOHDk6HZCrwtkx6jIgsAPKBvUCeiMwXkVoLSBpj/G/lypX079+fzp07s2nTJku29ZS3F73+F2gGnAk0Ac4CmrrbjTEOWrx4MYMGDeL0009nw4YNtG3btvY3GUd4O4bbB+iiqkfcP+8QkRuAn/wTljHGG2+88QajR4/m/PPPZ+XKlbRo0cLpkEwNvO3h5gOtKrS1BOxGCGMc8tJLLzFy5Eh69OjB6tWrLdkGgWOZFrZGRP4BpAAJwN3A8/4KzBhTvaeffpq//vWvXHXVVSxevJimTZs6HZLxgrcJdzbwK3Ad0M79/HFcVXiNMQH0xBNPMHHiRAYOHMjChQtp3Lix0yEZL1WbcEXkLVUd7v7xL6r6EpZgjXGMqjJz5kxmzJjBn//8Z1577TUiIyOdDsscg5rGcK+W/xY4+mcggjHGVE1VmTRpEjNmzGDMmDG8/vrrlmyDUE1DCh8An4jIDiDKPQ+3ElW93i+RGWMAKCkp4c4772Tu3LmMGzeOuXPn0qCBLWMSjGpKuMOAobgukCk2BcyYgCsuLubWW29l3rx53HPPPfzP//yPVdYNYtUmXFXNB14DEJFIVX0oYFEZYygqKiobPnjggQeYOXOmJdsg5+1aCjPct/GeBjSv8No6fwRmTDgrLCxkxIgRLFq0iEceeYTJkyc7HZLxAa8SroiMAf4F5ABHPF5SwJYjMsaH8vLyGDp0KCtWrGDOnDnceeedTodkfMTbebiPAENVdaU/gzEm3OXm5jJw4EDWr1/Pc889R1JSktMhGR/yNuE2BFb7MxBjwl1mZib9+/fnk08+Yf78+YwePbr2N5mg4u3ckseAB6ykjjH+cejQIXr37s2nn37Km2++ack2RHnbw70baAtMFJF0zxdUtaPPozImjBw4cIArr7yS77//nkWLFpGYmOh0SMZPvE24o/wahTFhau/evfTu3ZuUlBTeffddrrzySqdDMn7k7bSwjf4OxJhwk5KSwhVXXMGBAwdYtWoVl156qdMhGT+rafGaqao62/18ZnXbqep0fwRmTCjbuXMnvXr1Ijs7m/fff58LL7zQ6ZBMANR0Eewkj+cdanh4RUTiRGSxiOSKSIqIXOfFe9aJiIqIt0MfxtR73377LZdeeil5eXmsX7/ekm0YqenW3nEez2/wwbGeBgpxVf09B1guIltVdXtVG4vIyJriMyYYffnll1x11VU0bNiQjRs30rVrV6dDMgEUkGleItIMGAJMU9UcVf0QWApUOffFfRvxg8DEQMRnTCB8+umnXHHFFTRp0oRNmzZZsg1DgZpXeypQrKo7PNq2AmdUs/0jwDPAb/4OzJhA2LRpE7179yYuLo4PPviAU045xemQjAMClXCbA5kV2jKB6Iobikg34GLgqdp2KiJJIrJFRLYcPHjQJ4Ea42tr1qyhT58+nHTSSWzatImEhASnQzIOCVTCzQFiKrTFANmeDe472f4F3KmqRbXtVFWfV9VuqtqtVauKRYWNcd6yZcsYMGAAp5xyChs3bqR9+/ZOh2QcVNO0sCu82YGXyzPuABqKyCmqutPddjZQ8YJZDNANeMu97meEu/0XERmmqh94E5Mx9cHbb7/NddddxznnnMN7771HXFyc0yEZh9U0C+DFCj+3x7UcYzoQDwjwC14sz6iquSKyCJgpIjfhmqUwCLiowqaZuKoCl+oAfAacB9iYganXkpOTmTp1KqmpqcTFxZGens7FF1/M8uXLiY2NdTo8Uw/UNC2sc+lzEZmCK8lOU9UjItIUmIkr+XrrNlxVfw+43zdOVbeLSEfgW6CrqqbicaFMRKLcT/d7M8RgjFOSk5NJSkriyBHXctHp6ek0aNCAsWPHWrI1ZURVa99I5CDQTlWPerRFAr+qar0YPO3WrZtu2bLF6TBMmOrUqRMpKSmV2hMSEti9e3fgAzJOqrYOkrcXzXKBCyq0nU/56g/GhK3U1NRjajfhyds7uaYBq0RkGbAH19jqAGC8vwIzJlhkZGTQuHFj8vPzK73WsaOtXmr+y6serqq+ClwIfIdrJsH3QHd3uzFh66effuKPf/wjhYWFNGrUqNxrTZs2Zfbs2Q5FZuojr9cqUNVvReR7oI2q7vNjTMYEhU2bNjF48GBUlXXr1vHLL7+UzVLo2LEjs2fPZuTIkU6HaeoRb6v2tsB1Q8JQ4CjQTEQGAheo6gP+C8+Y+umVV14hKSmJLl268O6773LyyScDWII1NfL2otmzuObIJuBa8QvgE2C4P4Iypr4qKSlh0qRJ3HDDDVx66aV88sknZcnWmNp4O6TQC/e0MBFRAFU9KCKt/ReaMfVLbm4uo0ePZvHixdxyyy089dRTREZGOh2WCSLeJtxMoCVQNnbrvmHBxnJNWNi7dy+JiYls3bqVOXPmcMcdd+C+/dwYr3mbcOcB74jIVKCBiPwR1xKKz/otMmPqiS+++IKBAweSlZXF0qVL6d+/v9MhmSDl7RjuY8BCXFUbInHdorsE+Kef4jKmXli0aBGXXHIJDRs25OOPP7Zka+rE24TbRlXnqGpXVW2mqqer6hxc5XKMCTmqyqOPPsqQIUM4++yz+eyzzzjrrLOcDssEOW8T7o5q2r/1VSDG1BcFBQX85S9/YcqUKYwYMYL169fTpo31LUzdeZtwK10dEJEYoMS34RjjrLS0NHr37s2CBQt46KGHSE5OJioqqvY3GuOFGi+aicgeXGvgNhGRiqtwxANv+CswYwLtu+++Y8CAAezdu5c33niDa6+91umQTIipbZbCKFy92xWUr7CruNao/cFfgRkTSGvWrGHYsGE0btyYDRs20L17d6dDMiGoxoSrqhsBRKSlqtpSjCYkPfPMM9x+++107dqVZcuWWZFH4zfejuG+JiKXeDaIyCUi8m8/xGRMQBQVFXHHHXdw22230adPHz766CNLtsavvE24PYGPK7R9Alzu23CMCYysrCwGDhzIU089xd13382SJUuIjo52OiwT4ry90ywfaAZkebQ1x7VymDFBZdeuXSQmJvLDDz/w3HPPkZSU5HRIJkx428N9D3jOPRWsdErYXGCVvwIzxh8+/vhjLrzwQvbu3cuqVass2ZqA8jbh3our0sNhETkAHAJigbv8FJcxPpecnMzll19ObGwsmzdvplevXk6HZMKMV0MKqnoY6C8ibXHVM9ujqr/V8jZj6oWSkhJmzJjBww8/TM+ePXnnnXeIj493OixTj+3fv59du3ZRUFBA48aN6dy5s0/uNvS2h4uIxANXAper6m8i0k5ETqpzBMb4UV5eHiNGjODhhx9m7NixrF692pKtqdH+/fvZsWMHBQUFgOtW7x07drB///4679urhCsiPYEfgJG4KvgCnAI8U+cIjPGTffv2cdlll/H222/z+OOPM2/evEqFHo2paNeuXZSUlF+1oKSkhF27dtV5397OUpgDDFfVtSJy2N32KXBBnSMwxg+2bt1KYmIi6enpLFq0iGuuucbpkEw9V1RUxOHDh8t6thVV134svE24nVR1rfu5uv8sPIb3GxMwy5YtY8SIEbRo0YIPP/yQc8891+mQTD2Vn59Peno66enpZGRkoKrVbtu4ceM6H8/bhPmtiFytqu95tPUGvq5zBMb4iKryj3/8gwkTJnDeeeexZMkS2rVr53RYph5RVbKzs8uSbG5uLgBNmjShffv2xMfHk5+fz86dO8sNKzRo0IDOnTvX+fjeJtx7gXdFZDmulcOeAxKBQXWOwBgfKCwsZPz48cybN4+hQ4cyf/58mjZt6nRYph4oLi7m8OHDpKenc+jQIQoLXYXHY2Nj6dKlC/Hx8ZX+rYiIX2YpeJtwPwP+gGv1sJeAPcAFqvpLnSMwpo4OHTrE0KFDWb9+PVOnTmXmzJk0aOD1BBwTggoKCsoNFZSUlBAREUFcXBzx8fHExcXVWHG5TZs2fll0vtaEKyIRQA7QQlUf93kExtTBjh07GDBgACkpKSxYsIDRo0fX/iYTclSVnJycsiSbk5MDQFRUFCeeeCLx8fHExsY6/h9xrQlXVYtFZAeuBcd/9X9Ixnhn/fr1DBkyhIiICNauXUuPHj2cDskEUElJSbmhgtJZBDExMXTu3LlsqKA+lbP3dkghGdcY7j+BX/jvTAVUdZ0/AjOmJvPmzWPcuHGceuqpLFu2jC5dujgdkgmAwsLCsl7s4cOHKSkpoUGDBsTFxdGpUyfi4uLq9VxrbxPuOPefMyq0K2D/0k3AFBcXc//99/P3v/+dq666ioULFxIbG+t0WMZPVJXc3NyyJJudnQ24pmi1bduWuLg4TjjhBMeHCrzl7VoKdZ8PYUwd5eTkMHLkSJYuXcr48eOZM2cODRvaVPBQU1JSQkZGBocOHSI9PZ38/HwAoqOjSUhIoGXLljRr1qxeDRV4y/61mqCwZ88eEhMT+frrr3nqqaf461//6nRIxoeOHj3KoUOHSEtL4/DhwxQXF9OgQQNatGhBhw4diI+P98mNB06rNuGKyHeqerr7eWn13kpUtaOfYjMGgM8++4xBgwZx5MgRli9fTp8+fZwOydSRqpKXl0d6ejppaWlkZblqGzRq1IhWrVoRHx/PCSecQEREhMOR+lZNPdybPZ6PquuBRCQOeBG4CkgDJqvq61VsNwa4A9fiOFnA68AUVS2qawwm+CxcuJAxY8bQtm1b3n//fc444wynQzLHqaSkhKysLNLS0jh06BB5eXkANGvWjI4dOxIfH090dHRQDhV4q9qEq6ofejzf6INjPY1r/YU2wDnAchHZqqrbK2zXFNfC5p8CrYClwH3A33wQgwkSqsrs2bOZNm0aF198MYsXL6ZVq1ZOh2WOUVFRUdlY7KFDhygqKkJEaNGiRdmttFFRUU6HGTA1DSnM9GYHqjq9tm1EpBkwBDhTVXOAD0VkKTAamFRhf55LPu4VkWSsWGVYyc/P56abbiI5OZlRo0Yxb968kBi/CxelQwXp6elkZmaiqkRGRhIfH182VBCuFztrOusOHs+jcCXMz4EUoCOupRnf8fI4pwLFqrrDo20rrmrAtbkUqNgLBkBEkoAkgI4dbSg5FBw4cIA//elPfPzxx8yaNYspU6aE9FfMYFSxGkKnTp1o0qRJWZI9cuQIAE2bNuWkk04iPj6emJgY+xypeUjhhtLnIvImMEJV3/FoGwwM8/I4zYHMCm2ZQI11qUXkBqAbcFM1MT4PPA/QrVu36tdVM0Hhm2++ITExkf379/P2228zdOhQp0MyFZRWQyhdSaugoIAffvgBcC34EhsbW3YrbZMmTZwMtV7ytl/fF1e1B09LgJe9fH8OriKUnmKA7OreICLX4Bq37a2qaV4exwSplStXMnz4cJo1a8bGjRs5//zznQ7JeCguLiYrK6vSsoWlGjZsyIUXXhi2QwXe8vZv50dgPPC/Hm23AT95+f4dQEMROUVVd7rbzqb6oYI+wAtAf1W1NXdDmKoyd+5c7rrrLv7whz+wdOlSOnToUPsbjV+VJtiMjAwyMjLIzs6ucXHuoqIiS7Ze8PZv6CZgsYhMBPYC7YEiYLA3b1bVXBFZBMwUkZtwzVIYBFxUcVsRuQLX2g1/UtXPvIzPBKGioiLuvPNO/vWvfzFw4ECSk5Np3ry502GFpeLiYjIzM8nMzKyUYKOjoznppJNo0aJFueKKnuyipne8vbX3SxE5BegOtAP2AZ+o6tFjONZtuNbSPQCkA+NUdbuIdAS+BbqqaiquIpWxwAqPQfYPVLXvMRzL1HMZGRkMHz6c1atXM2HCBB599NGQm+Ren3mbYGNiYsr1XDt37lxuDBd8Vw0hHHj9HcCdXD843gOp6iHgmiraU3FdVCv92aaAhbiff/6ZAQMGsHPnTl588UXGjh3rdEgh73gTbEWli3L7oxpCOLBBF+N3ycnJTJ06ldTUVFq3bk1OTg6NGzdmzZo1XHbZZU6HF5J8lWCr4q9qCOHAEq7xq+TkZJKSksrmZu7fvx8R4aGHHrJk60OlCTYjI4PMzEyfJljjO1LTlcdg0q1bN92yZYvTYZgKOnXqREpKSqX2hIQEdu/eHfiAQkRtCbZFixaWYJ1T7R0e9kkYv/ntt9+qTLYAqampAY4muFWXYEWkXA82NjbWLj7WY5Zwjc/l5OTw97//nSeeeKLabexW7JpZgg1NlnCNzxQVFfHiiy/y4IMPsn//foYNG0b37t2ZNm1a2RguuO6xnz17toOR1j+eCTYjI4OcnBxLsCHIEq6pM1Vl6dKlTJo0ie+//54ePXrwn//8h+7duwOuq9qlsxQ6duzI7NmzGTmy4p3i4cUSbHiyi2amTjZv3syECRP48MMP+f3vf89jjz1GYmJiWK4MVXEVLc/5qRUTbGkxxNIEW5pcLcGGBLtoZnxr586dTJkyhX//+9+0adOG5557jrFjx4btFfHqVtE6cOAAR48erZRgO3bsaAk2DIXnb4c5bgcPHuThhx/mmWeeoXHjxjz00EPcc889Yb8Gws8//1xpFS1V5dChQ8TExFiCNYAlXOOlI0eOMGfOHP72t79x5MgRkpKSmD59Om3btnU6tIArLi4mJyeHrKwssrKyyM7OprCwsNrtzz333ABGZ+ozS7imRsXFxcyfP59p06bx66+/cs011/Doo4/y+9//3unQAkJVyc/PL0uuWVlZ5Obmlt1kEBUVRWxsbFm9ropsFS3jyRKuqZKqsnLlSu6//36++eYbunfvzltvvUWPHj2cDs2vioqKynqtpQm2NJFGREQQHR1Nhw4diI6OJiYmhkaNGgGVx3DBVtEylVnCNZVs2bKFiRMnsn79ek4++WT+/e9/M3jw4JCbeaCq5ObmlkuwFecLt2zZkpiYGKKjo2nWrFm1fwe2ipbxhiVcU2bXrl1MnTqVN954g1atWjF37lySkpKIjIx0OjSfKCwsLDfump2dTXFxMQCRkZFER0fTunXrsgRrq2gZX7OEa0hPT2f27Nk8/fTTRERE8MADDzBhwgRiYiqWoQseJSUllS5s5efnA66pWc2bN6dNmzbExMQQExNDVFRUyPXgTf1jCTeM5eXl8dRTT/HII4+QnZ3N2LFjmTFjBu3bt3c6tGOiqhQUFJS7sFV65xa4LlzFxMTQrl07YmJiaN68uU3NMo6whBuGSkpKyhYF37NnDwMGDOBvf/sbZ5xxhtOheaWoqKhsSKA0wR496qr21KBBg7JbY0uHBmymgKkvLOGGmTVr1jBx4kS++uorunXrxvz587n88vpb1UhVOXLkSLkLW7m5uWWvN2nShLi4uLLk2rx5cxsaMPWWJdwwsXXrViZOnMjq1avp3Lkzb7zxBn/+859p0KCB06GVc/To0XLjrllZWWUXtho2bEh0dHS5mQOhckHPhAdLuCEuNTWVadOm8eqrr3LCCSfw5JNPMm7cuIB+za5uUZeSkpKyaVmlCTYvL6/sfaUXtkrnvDZp0sR6ryaoWcINURkZGTz66KP885//BGDixIlMmjSJFi1aBDSO6hZ12b17N4WFhWXtjRo1IiYmhrZt25b1Xu3Clgk1lnBDTEFBAf/617+YNWsWhw8f5vrrr2fmzJkBr7BQelPBjz/+WOWiLgUFBbRv377chS3rvZpQZwk3RJSUlPDWW28xZcoUdu/ezdVXX81jjz3G2WefHZDjFxcXlw0NZGZmlht7rYqq8rvf/S4gsRlTX1jCDQHr169nwoQJfPHFF5xzzjmsXr2aK6+80q/HLCwsJDMzsyy5es57bdasGa1btyY2Npaff/65ypW0bKqWCUeWcIPYN998w6RJk1i+fDkdOnRgwYIFjBw50uczDzynZpUm2NKLW6XzXjt06EBMTAyxsbGVbom1RV2McbGEG4T27t3Lgw8+yMsvv0x0dDSPP/44t99+O1FRUT7Zf0lJCdnZ2WXJNTMzs2zFrMjISGJiYjjxxBOJjY2lefPmNSZ4W9TFmP+yhBtEsrKyePzxx/nHP/5BcXExd911F1OmTCE+Pr5O+y2d+1o6RFBakhtcNxaUznuNjY09rqlZtqiLMS6WcINAYWEhzz//PA899BBpaWlcd911zJo167i+lpcuqO3Zey1dkrB0UZf27dsTGxtbbr1XY0zdWcKtx1SVd955h8mTJ/Pjjz9y+eWX88QTT3Deeed5vY/Smws8L3CVXsSKiIggNja27AKXzX01xr8s4dZTH374IRMmTGDz5s2ceeaZrFixgj59+tT6db60YkFpcs3Kyiq7YBUVFVVWjjsmJqbGBbWNMb5nCddhpat2paam0rFjR8aPH89HH33EkiVLaNeuHS+99BLXX399tT3PgoKCStOzSjVv3pwTTzyxbPzVpmIZ4ywpvTgS7Lp166ZbtmxxOoxjkpycTFJSUrmyLuCaozp9+nTuuusumjZtWtZeeveW5/hrQUEB4JpqVZpYSxfVPtaKBcYYn6j2a6Ml3ABTVQ4fPsyuXbvo06cPZ599NjfffDOtW7fmwIEDvPDCC/zwww/s2bOH4uLisulZFe/eatSoUVlyLZ2eZcMDxtQLlnAD6ciRI+zatavaR1ZWFgC9evViwoQJ5ebPFhYW8sUXX3DVVVeVu3uradOmxMbGliVZKwljTL3lfMIVkTjgReAqIA2YrKqvV7Pt3cD9QBPgHWCcqhbUtP9AJtyjR4+SmppabUI9cOBA2bYRERG0bt2a008/nZNPPpmEhATat29P69atAcoNGZRS1bLkWppgbd1XY4JGtQk3kIN8TwOFQBvgHGC5iGxV1e2eG4nI1cAk4ArgV2Ax8JC7zScWLVqEqhIXF8ehQ4cQEQYPHlz2eklJCb/99lulRLpv3z7S09PJzc0tW6M1JiaGFi1acMopp9CjR4+y6gNRUVG1jqHW9J/dueee66vTNcbUEwHp4YpIM+AwcKaq7nC3vQrsVdVJFbZ9HditqlPcP/cCklW1bU3H8LaHu2jRIpo1a1bpa/z27dvJzs6msLAQVaV58+ZlSweWzlGt6RbWhg0bEhkZWenPmtq2bNlSdtHLU+PGjenevXut52KMqZcc7+GeChSXJlu3rUDPKrY9A1hSYbs2IhKvqul1DURVK6050KhRo7IeZWFhIUVFRYgIjRo1olmzZmU91poS6fGMp3bu3NkWdjEmjAQq4TYHMiu0ZQLRXmxb+jwaKJdwRSQJSAK8XmA7Li6uyvaSkhJ69uwZ0BpftrCLMeElUAk3B4ip0BYDZHuxbenzStuq6vPA8+AaUvAmkEOHDtGyZcsq250oqGgLuxgTPgKVYXYADUXkFI+2s4HtVWy73f2a53b7fTGcAK4FWvLz88u15efn2xQrY4zfBSThqmousAiYKSLNRORiYBDwahWbLwBuFJGuInIC8ADwiq9iGTx4MLm5uaSlpVFSUkJaWhq5ubnlZikYY4w/BHJa2G3AS8ABXGOx41R1u4h0BL4FuqpqqqquEpHHgfX8dx7ug74MxJKrMcYJdqeZMcb4VrXjk4G/SmSMMWHKEq4xxgSIJVxjjAkQS7jGGBMgIXPRTEQOAinH+LaWuFYuqw/qUyz+Eg7nCHaeoeR4zjFNVftU9ULIJNzjISJbVLWb03FA/YrFX8LhHMHOM5T4+hxtSMEYYwLEEq4xxgRIuCfc550OwEN9isVfwuEcwc4zlPj0HMN6DNcYYwIp3Hu4xhgTMJZwjTEmQMIu4YpIYxF5UURSRCRbRL4Ukb4OxrNBRPJFJMf9+MGpWHxFRP4qIltEpEBEXqnwWi8R+V5EjojIehFJcCjMOqvuPEWkk4iox2eaIyLTHAz1uNX2+xIKn2dN5+jrzzKQyzPWFw2BPbjqqaUC/YCFInKWqu52KKa/quo8h47tD78Cs4CrcS2xCYCItMS1LvJNwDLgYeAtIFgrZlZ5nh5aqGpRYEPyuWp/X3BVZwmFz7Omcyzlk88y7BKuezH0GR5N74rILuA8YLcTMYUaVV0EICLdgJM8XhoMbFfVt92vzwDSROT3qvp9wAOtoxrOM2TU8vsSTwh8nrWc4xe+PFbYDSlUJCJtcFUVrqrcT6A8KiJpIvKRiFzmYBz+dgauKsxA2T/0n9ztoShFRH4RkZfdvfugV+H3JSQ/z2pygk8+y7BOuCISCSQD8x38H/l+oAvQHtecv2Ui8juHYvG3Y6neHMzSgPOBBFy9pGhc/86CWhW/LyH3eVZxjj79LMM24YpIA1w11QqBvzoVh6p+qqrZqlqgqvOBj3CNIYWiY6neHLRUNUdVt6hqkarux/Xv6yoRqXjuQaOa35eQ+jyrOkdff5ZhmXDFVaL3RaANMERVjzockielhhIdQa5cRWYRaQb8DmeHcwKh9O6ioPxca/h9CZnP8xhyQp0+y7BMuMAzwOlAoqrmORWEiLQQkatFJEpEGorISOBS4D2nYvIF97lEARFAROn5AYuBM0VkiPv16cC2YLrA4qm68xSRC0XkNBFpICLxwP8CG1S14tfvYFHd70sofZ5VnqPPP0tVDasHrrEYBfJxfSUqfYx0IJZWwOe4voJlAJuBK53+O/LBec1w/x17Pma4X+sNfA/kARuATk7H6+vzBEYAu4BcYB+wAGjrdLzHeY41/r6EwudZ0zn6+rO0tRSMMSZAwnVIwRhjAs4SrjHGBIglXGOMCRBLuMYYEyCWcI0xJkAs4RpjTIBYwjXGR0TkFhGZ4+N9NnavN9val/s1zrB5uMb4gIg0wrVSVndV3evjfU8E2qjqvb7crwk86+EaU4G4HOvvxiDge18nW7fXgTEi0tgP+zYBZAnXBJSI7BaRCSKyTURy3aVN2ojISnd5k/dF5ASP7buLyMcikiEiWz3XCxZXeaJZ7tdzRGSZiMSLSLKIZInI5yLSyWP7i9xtme4/L6qwr9ki8hFwBLhXRMotPi0i94rIf6o5tb7ARo9tS0uz3CAie0TksIjcKiLnu889Q0Tmemx/sohsdMeWJiJvlb6mqr8Ahwm+SgqmIqfvY7ZHeD1wVdXYjGtVpvbAAeD/gHOBxsA64EH3tu2BdFzLVTYArnT/3Mr9+gbgR1wrVMUC3wI7cN3f3xDXfe8vu7eNw5W0RrtfG+H+Od5jX6m4Fs9u6I7lEHC6R+xf4lpJqqrz+hwY5vFzJ1z35z8LRAFX4bpX/z9Aa49z7+ne/g1gqvs8o4AeFfa/FLjD6c/PHnV7WA/XOOEpVd2vrq/fHwCfquqXqlqAawWqc93bjQJWqOoKVS1R1TXAFsqvF/yyqv6krtWbVgI/qer76qo/9bbHvvoDO1X1VXWtbfoGrkVXEj329Yqqbne/XoCrPtcoABE5A1cSfbeac2pB1evAPqyq+aq6GtcCKG+o6gGPcy+N7yiuRVTaubf/sMJ+st3HMEHMEq5xwn6P53lV/Nzc/TwBGOb++p0hIhlAD+DE49hXOyClQhwpuHqapfZUeH0+cJ17rdTRwEJ3Iq7KYaqudOBtfBNxrbH6mYhsF5GxFfYTjWtFORPEwq6IpAkqe4BXVfVmH+zrV1wJ3FNHYJXHz+Wm7KjqZhEpBC4BrnM/qrMNVx2s46KqvwE3A4hID+B9Edmkqj+6Nzkd+Pvx7t/UD9bDNfXZa0Cie5H20gW+LxOR46mQuwI4VUSucy8SPhzoSvVDBKUWAHOBoiq+5lfcf8/jiAsAERnmcV6HcSX/Yvdr7XGNQW8+3v2b+sESrqm3VHUPrulWU4CDuHq8EziOf7eqmg4MAO7FdeFtIjBAVdNqeeurwJnuP2uyDPi9iLQ71tjczgc+FZEcXBfI7lTVXe7XrsNV1LC64QwTJOzGB2NqICJNcM0m+H+qurOWbZOArqp6lw+P3xhXKfJLVfWAr/ZrnGEJ15gaiMg9uHrCVzgdiwl+dtHMmGqIyG5cMweucTYSEyqsh2uMMQFiF82MMSZALOEaY0yAWMI1xpgAsYRrjDEBYgnXGGMCxBKuMcYEyP8HGQO2vCKloZ4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the DI values for the memory range\n",
    "fig, ax = plt.subplots(figsize = (5, 4))\n",
    "ax.plot(memory_range, di_memory[(0, 1)], color = 'k', marker = 'o', label = r\"DI($n_1\\rightarrow n_2$)\")\n",
    "ax.plot(memory_range, di_memory[(1, 0)], color = 'silver', marker = 'o', label = r\"DI($n_2\\rightarrow n_1$)\")\n",
    " \n",
    "# plot clean up \n",
    "ax.set_xlabel('memory (ms)')\n",
    "ax.set_ylabel('directed information')\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.set_xticks(memory_range)\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "ax.legend()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15032bd7",
   "metadata": {},
   "source": [
    "As we increase the memory size, the directed information also grows. It's essential to choose an appropriate memory size to effectively capture the interactions within spiking data. For instance, at a memory size of 2 ms, we observe no discernible difference between the causal and non-causal directions. This observation aligns with our simulation, where we modeled a time delay of ~10 ms for the affected neuron to receive the input. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee114abe",
   "metadata": {},
   "source": [
    "#### Kernel width of the Gram matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee91d76b",
   "metadata": {},
   "source": [
    "Next, lets look at the effect of the kernel width of the Gram matrix on DI. The Gram matrix captures the variance in the data and is an important component for entropy estimation. \n",
    "\n",
    "\n",
    "For a variable $X = \\{\\textbf{x}_i\\}_{i = 1}^{N} $, the Gram matrix, $K$, is an $N \\times N$ matrix whose entries are given by $\\kappa(\\textbf{x}_i\\ ,  \\textbf{x}_j), i, j = 1, ..., N$, where $\\kappa$ is a kernel function. This function dictates the extent of similarity between various realizations of $X$. \n",
    "\n",
    "For our puposes, we use a Gaussian kernel to evaluate the Gram matrix of each of the input variables and its kernel width dictates the extent of similarity between data points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d9399054",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2min 51s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "trial = 8 # picking a random trial\n",
    "kernel_widths = [1e-1, 2, 5, 10, 20, 30] # range of widths to sweep over\n",
    "memory = 15\n",
    "\n",
    "# Initialize\n",
    "di_kw = {(n1, n2): np.zeros(len(kernel_widths )) + np.nan for n1 in range(num_units) for n2 in range(num_units) if n1 != n2}\n",
    "\n",
    "# Evaluate DI for different memory values\n",
    "for n, width in enumerate(kernel_widths):\n",
    "    for n1 in range(num_units):\n",
    "        for n2 in range(num_units):\n",
    "            if n1 != n2:\n",
    "                model = directed_information(memory = memory, multiplier_X = width, multiplier_Y = width)\n",
    "                di_kw[(n1, n2)][n] = model.DI(rates[n1, trial], rates[n2, trial])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "3df31528",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAAEUCAYAAACvXaa+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABAgUlEQVR4nO3deXxU1fn48c8z2QNJIIBohIRN+IqiKLiC4lIUERCXfrEGQb4qaLUIhP6Ka1sLVSyouFRFrRJEwaIoKKCIlSqKIqWyKXGBgIghrGFJQkie3x8ziZPJTDJJZskkz/v1ui9mzpw595m5CU/OPfecK6qKMcYYYypzhDsAY4wxpiGyBGmMMcZ4YQnSGGOM8cISpDHGGOOFJUhjjDHGC0uQxhhjjBfR4Q4gVAYMGKBLly4NdxjGGGMaFvH1QpPpQe7evTvcIRhjjIkgTSZBGmOMMbVhCdIYY4zxwhKkMcYY44UlSGOMMcaLJnMVqzHGBFpJSQk//vgjRUVF4Q7FVCM+Pp527doRExNTuzeqapPYevXqpfXxyiuvaEZGhoqIZmRk6CuvvFKv9oLdbmOJx5iG7IcfftD8/HwtKysLdyjGh7KyMs3Pz9cffvjBVxWfeSPsiStUW30S5CuvvKKJiYkKVGyJiYn1Th7BarexxGNMQ7dp0yZLjhGgrKxMN23a5Otln3lDVJvG/SB79+6tX375ZZ3e26FDB3Jzc6uUR0dH07lz5zrH9P3333Ps2DGv7Xbp0gWRX+avlj+urqy+9detW8fRo0erxJORkcHWrVtr9+GMaQK+/vprTj755HCHYfxQzbHyuVCAjUH6Ydu2bV7Ljx07Rs+ePevc7ubNm32226NHj4rn5X/EuP8x41lW19fcH3tLjgC5ubl8++23nHTSSf58LGOMaRQsQfohPT3daw8yIyODuXPn1rndVatW+Wz39ddfr3O7deWrpwzQtWtXzjvvPEaOHMmwYcNo0aJFaIMzxpgQs2kefpgyZQqJiYmVyhITE5kyZUqDbDfQ8Tz55JNMnTqVAwcOcNttt3H88cczbNgwFi9e7PUUsTEmct199908/vjj4Q6jVs4++2w2btwY+IarG6BsTJtdxVr/eMrKynT16tV65513ampqqgJ6/PHHa1ZWlq5bty6MURsTHtVc+NEgZGRkaHx8vDZv3lxTUlL0vPPO02eeeUZLS0sr1Vm2bJmqqu7atUvT0tL0yJEj4Qq5TubNm6fXXHNNtXXqcpFO2BNXqLb6JkhTWXFxsb755pt61VVXaXR0tAJ65pln6owZM3TXrl3hDs+YkIiEBFme/Pbv369vv/22dujQQW+66SavdR555BG95ZZbwhJrfRQWFmrLli31p59+8lmnLgnSTrGaOomNjeXqq6/mrbfe4qeffmLGjBkA3HXXXaSlpTF06FAWLFjg88IfY5qyOXPm0KFDBxwOBx06dGDOnDlB32dKSgpDhgxh3rx5zJo1iw0bNlSps2TJEvr161ep7Nlnn+XKK6/kjjvuoHXr1qSlpbFs2bJ6x6PqfQZFXfYXHx9Pr169eP/99+sdlztLkKbe2rRpw9ixY1mzZg3r1q3jrrvuYtWqVVxzzTWkpaXxu9/9ji+//NLnL4QxTcmcOXMYPXo0ubm5qCq5ubmMHj06JEkSnON17dq14+OPP67y2vr16+nWrVulsnXr1vHZZ58xZMgQdu3axZgxY5g6dWq9Yvj4448ZOHCg1xWI6rq/k08+ma+++qpecXmyq1hNQPXo0YNp06bx8MMP8/777zNr1iyef/55nnrqKbp3787IkSMZPnw4aWlp4Q7VmIAaN24c//3vf2ust2rVKoqLiyuVHTlyhJtvvpnnn3++2vf27NkzIBfQpKWlsXfv3irl+/fvJykpqVLZunXrmDRpEpdffjkA3bt3r0iuBw4coH///mzatIlVq1Zx6qmn+rX/Pn36cNxxxzFkyBAWLlxIfHy8X/v77LPPmDBhArGxsaSlpZGdnV2xfFxSUhI7d+6s5TdRPetBmqCIjo5m4MCBzJs3j507d/Lss8+SkpLCH/7wB9q3b8+AAQN47bXXKCwsDHeoxoSUZ3KsqTwYduzYQWpqapXyli1bcvDgwUpl69evZ/DgwRXPN2zYQPfu3QHnVe7vvvsu1113ndf9LF26FBGpskVFRZGdnc2yZct45pln/N5fRkYGH374IStWrKBTp068/fbbFfUOHjwY8OlnIetBikgq8CJwGbAbuFtVX/VS73rgz8DxQDGwBPidqha4Xv8IOBcon1+wQ1W7ebZjGo6WLVsyZswYxowZQ05ODtnZ2WRnZ3PDDTeQnJzMsGHDGDlyJOeff36l1X2MiST+9ux8zTfOyMjgo48+CmxQXqxevZodO3bQt2/fKq+ddtpp5OTkcNZZZwGwZcsWjh07Vum069q1axk6dCgAMTExtGnTxue+BgwY4HVopaysjFGjRrFz505uu+22ivKa9ud+5ik6OhqH45c+3tdff83w4cNr+PS1E8oe5NPAUaAtkAk8IyKneKm3EuijqilAJ5xJfLJHnTtVtblrs+QYQbp27crkyZPZunUry5cvZ+jQocyZM4e+ffvStWtX/vKXv9iydqZRC9f854KCAt555x2uv/56hg8fXmm1rnIDBw5kxYoVFc/XrVtHjx49KiWitWvXcvrpp9crlpUrV5KXl8fbb79NQkJCrfe3ZcsWlixZwqBBgwBn73vNmjX079+/XnF5CkmCFJFmwLXA/ap6SFU/ARYCN3rWVdXtqrrbragU6BKKOE3oOBwOLrnkEmbNmkVeXh4vv/wy7dq144EHHqBjx45cfPHFvPzyyxw6dCjcoRoTUJmZmcycOZOMjAxEhIyMDGbOnElmZmZQ9jd48GCSkpJo3749U6ZMYcKECbz00kte644YMYLFixdXDH2sW7eu0nKae/bs4eeff/Z7rNGXCy64gCVLllRKjv7ur6CggJEjRzJ79mxiY2MBWLhwIRdddFHgr22obg5IoDbgDKDQo2wisMhH/b7AAZx3lDgMXOb22kdAPs7TtCuBi6rZ72jgS+DL9PR0n/NjTMOxZcsWffDBB7Vz584VdxO58cYb9YMPPqg0udmYhqChz4Osi7vvvlsfe+yxWr1n5MiRun79+uAE5KakpEQHDhyoy5cvr1R+9tln17j/BrtQAHAB8LNH2a3ARzW870TgT0BXt7JzgCQgDhgJHAQ61xSDLRQQWcrKyvSTTz7RW2+9VZOTkxXQ9u3b6z333KObN28Od3jGqGrjTJC1dcUVV+gJJ5yg5557rr700ktB3Vd2dra2atVK+/Xrp/369dO5c+f6/d4Ge7srETkDWKmqiW5lWa7e32Df7wQRORf4u6qe6eP1pcC7qvpkde3U53ZXJrwKCwt5++23mTVrFu+//z5lZWWce+653HTTTbZwugkru91V5KjL7a5CdZFODhAtIu73Szod8Gd12WigupsuKtV8QBP5EhISuP7661myZAnbt2/nkUce4eDBg7ZwujEmqEKSIFX1MPAm8KCINBORPsBVwGzPuiKSKSLp4pQBTAGWu15rISKXi0i8iESLSCZwIfBeKD6HCb+0tDR+//vfs379elavXs2tt97K8uXLufLKK2nfvj0TJ05k/fr14Q7TGNMIhHKax2+BBGAX8Bpwu6pudCXDQyKS7qrXHfgUOITzIpzNOMcrAWJwTvkov0jnd8BQVfV+52HTaIkIvXv35sknn+Snn37izTff5JxzzmHGjBmcdtppnHnmmcyYMYP8/Pxwh2qMiVAhGYNsCGwMsmnIz8/ntddeY9asWfznP/+pWNFn5MiRDBo0qOKycGMCwcYgI0dDHoM0JiTcF05fv34948aN44svvuDaa6/lhBNO4M4772T16tWoaljuqGCMiRzWgzSN3rFjx1i2bBmzZs3irbfeori4mLS0NPLz8ykpKamol5iYGNQJ26bxsR5k5LAepDFeREdHc8UVVzB37lx+/vlnnnvuOfbs2VMpOYLzjgr33ntvmKI0xjQ0liBNk9KiRQtGjx7t80bO27ZtC3FExpiGyhKkaZLS09O9lqekpNiNnY0xgCVI00R5u6NCVFQU+/fv54YbbuDw4cNhisyY8Lr77rsDclPmUDr77LPZuNGfdWdqxxKkaZK83VFh1qxZPPTQQ8ybN48+ffqwZcuWcIdpTL106NCBhIQEkpKSaNGiBeeffz7PPvssZWVllep88MEHgHOaVHZ2NmPGjAlXyHUyceJEHnjggYC3awnSNFmZmZls3bqVsrIytm7dSmZmJpMmTWLx4sXk5ubSu3fviv84jIlUixYt4uDBg+Tm5jJp0iSmTp3KzTff7LXuyy+/zMCBA6vchqqhGzJkCP/617/YuXNnQNu1BGmMhwEDBrB69WrS0tK4/PLLmTZtmo1LmoDKy8tj1apVrFixglWrVpGXlxf0faakpDBkyBDmzZvHrFmz2LBhQ5U6S5YsoV+/fpXKnn32Wa688kruuOMOWrduTVpaGsuWLat3PL5+p+qyv/j4eHr16sX7779f77jcWYI0xosuXbrw2Wefcc011/D73/+ezMxMjhw5Eu6wTCOQl5dHTk4OxcXFABQXF5OTkxOSJAnO8bp27drx8ccfV3lt/fr1dOvWrVLZunXr+OyzzxgyZAi7du1izJgxTJ06tV4xfPzxxwwcOJCioqIqr9V1fyeffDJfffVVveLyFB3Q1oxpRJo3b87rr7/O1KlTueeee9i0aRMLFiygY8eO4Q7NNEDfffcdhw4dqrFeQUFBld5TWVkZmzdvrvEUYfPmzenSpUu94gTnov979+6tUr5//36SkpIqla1bt45JkyZx+eWXA9C9e/eK5PrZZ58xYcIEYmNjSUtLIzs7m5iYmBr336dPH4477jiGDBnCwoULiY+P92t/Bw4coH///mzatIlVq1Zx6qmnVrwvKSnJTrEaE0oiwqRJk3j33XdtXNIEhK9Ti6E8jb9jxw5SU1OrlLds2ZKDBw9WKlu/fj2DB/9y294NGzbQvXt3ADIyMvjwww9ZsWIFnTp14u2336703qVLlyIiVbaoqCiys7NZtmwZzzzzjN/7S0xM5N133+W6666rEvvBgwcDfm9Y60Ea44crrriC1atXM3ToUC6//HKmTp1KVlYWInYrUuPkb89u1apVFadX3cXFxdGzZ88AR1XV6tWr2bFjB3379q3y2mmnnUZOTg5nnXUWAFu2bOHYsWOVTruuXbuWoUOHAs6eaLno6Ggcjsp9rgEDBnhN/GVlZYwaNYqdO3dy2223VZTXtL+YmBjatGnj9XN9/fXXDB8+vIZPXzvWgzTGT126dGHVqlU2LmnqpWPHjlUSicPhCPqp+4KCAt555x2uv/56hg8fTo8eParUGThwICtWrKh4vm7dOnr06FEp3rVr13L66adXet+WLVtYsmQJgwYN8iuWlStXkpeXx9tvv13pill/9+epuLiYNWvW0L9/f7/27y/rQRpTC+Xjkg8//DD33nuvjUuaWmvbti3gTCrFxcXExcXRsWPHivJAGzx4cEXvrnv37kyYMKFSr83diBEj6NmzJ4WFhSQkJLBu3bpKvdo9e/bw888/Vxr7KygoYOTIkcyePdvv28ldcMEFLFmypMoZGH/2583ChQu56KKLKvVoA0JVm8TWq1cvNSaQFi9erC1atNDU1FRdtmxZuMMxYbBp06ZwhxBwd999tz722GN+1S0pKdGBAwfq8uXLgxuUh5EjR+r69esrnp999tmVnntTzbHymTfsFKsxdVQ+LnnCCSfYfEnTaPz1r39l3LhxftV97bXX+Pzzz3nwwQe56KKLmDdvXnCDw3ka+P333+fWW2/l5ZdfBuDzzz+vsZdZF3Y/SGPq6dChQ4waNYr58+fzm9/8hhdeeKHKOq+mcbL7QUYOux+kMWFQPi7517/+lblz53L++efbOq7GNAIhS5AikioiC0TksIjkisgNPupdLyKbReSAiOwSkVkiklzbdowJJRHh7rvvtvmSxjQioexBPg0cBdoCmcAzInKKl3orgT6qmgJ0wnml7eQ6tGNMyHmOS06fPt3GJY2JUCFJkCLSDLgWuF9VD6nqJ8BC4EbPuqq6XVV3uxWVAl1q244x4eI+X3LixIk2X9KYCBWqHmRXoFRVc9zKvgK89vxEpK+IHAAO4kyIj9exndEi8qWIfJmfn1/Pj2CM/2xcsumwMwQNX12PUagSZHPggEfZASDJS11U9RPXKdZ2wN+ArXVsZ6aq9lbV3r6WJzImWGxcsvGLioqipKQk3GGYGpSUlBAdXft1cfxKkK4LY/4qIotF5N/um5/7OQQke5Ql4+wh+qSqO4ClwNz6tGNMONm4ZOPVokUL8vLyKCsrC3coxoeysjLy8vJISUmp9Xv9TamvAnHA60BdBlNygGgROUlVv3WVnQ5s9OO90UDnALRjTNiU319y1KhRTJw4kTVr1th8yUagdevW/Pjjj2zevDncoZhqNGvWjNatW9f6fX4tFCAiBUAbVa26BL2/OxKZCyhwC9ATWAycr6obPeplAh8D24F0IBvYo6rX1KYdT7ZQgGkIVJWHHnqI++67j9NOO83WcTUm/Oq9UMA6nOOB9fFbIAHYBbwG3K6qG0UkXUQOiUi6q1534FOcp1NXApuBW2tqp56xGRMSIsI999xj45LGRAB/e5APAr8BXgJ+dn9NVf8RnNACy3qQpqH57rvvGDp0KF9//TWPPPIIEyZMsPtLGhN6Pn/p/E2Q//LxkqrqJXWNKpQsQZqG6ODBg4waNYo33njD1nE1JjzqlyAbA0uQpqGycUljwqr+i5WLSEsRGSEid7v+bRmY2Ixp2ryNSy5fvjzcYRnT5Pk7D/I84HvgNuA0YAzwvavcGBMA7vMlL7vsMh599FGbL2lMGPnbg3wc+K2qnq+qv1HVPsDtwBNBi8yYJqh8vuTVV19NVlYWw4cPt3VcjQkTfxNkV5yLBLibj2sRcWNM4CQlJfHPf/6TKVOm8Nprr9GnTx+2bt0a7rCMaXL8TZDfAtd7lP0a52lXY0yAuY9LbtmyxcYljQkDfxPkOOApEVklIvNE5HPg78DYoEVmjKkYl2zbtq2NSxoTYn4lSFX9FOd6qE8Ba4AngS6ucmNMEJ100kmsWrWKoUOH2rikMSFk8yCNiRDu8yVPP/10FixYQIcOHcIdljGRrvYLBYjIUlUd4Hr8Mc4FwqtQ1QsDEWGwWYI0jcXixYu54YYbiI6OZt68eVx66aXhDsmYSOYzQVZ3u6tst8cvBC4WY0x9DBw4kNWrVzN06FAuu+wy/va3vzF+/Hhbx9WYAPOZIFX1Vben36jq5551ROTsoERljKlW+bjkTTfdRFZWFmvWrOH555+3dVyNCSB/r2Jd5qN8aaACMcbUTlJSEvPnz7f5ksYESbUJUkQcIhLlfCjiel6+nQQcC02YxhhvyudLvvPOOzZf0pgAq6kHeQw4CiS6Hpe4bZtwzoU0xoRZ+bikzZc0JnBqSpAdcc5//BHo5LZ1BJJV9U9Bjc4Y4zfP+ZI33nijzZc0ph6qTZCqmquqW1U1w/W4fNumqoWhCtIY4x/3cclXX33VxiWNqQe/FwoQkSFAP6A1bvNGVHVEcEILLJsHaZoa9/mSr7/+Opdcckm4QzKmIarfDZNF5I/Ac676vwb2AJcD+/2OQCRVRBaIyGERyRWRG3zUGykia0SkQER+FJFHRCTa7fWPRKRIRA65ts3+xmBMU+I+Ltm/f38blzSmlvyd5vF/QH9VHQ8cdf07GOhQi309jfOCn7ZAJvCMiJzipV4izsXRWwPnAJcCEz3q3KmqzV1bt1rEYEyTUj4uedVVV9m4pDG15G+CbKGqG1yPj4pIjKp+gfOUa41EpBlwLXC/qh5S1U+AhcCNnnVV9RlV/VhVj6rqDmAO0MfPOI0xHsrHJSdPnmzjksbUgr8J8nu33t4G4HYRuRHY5+f7uwKlqprjVvYV4K0H6elCYKNH2UMisltEVorIRX7GYEyT5XA4uPfee1m0aFHFfMkPP/ww3GEZ06D5myDvA1q5Ht+N8z6QfwMm+Pn+5sABj7IDQFJ1bxKRUUBvYJpb8R9wTjU5EZgJLBKRzj7eP1pEvhSRL/Pz8/0M1ZjG68orr7RxSWP8FJLbXYnIGcBKVU10K8sCLlLVwT7eMxTnhUG/UtX11bS9FHhXVZ+sLga7itWYXxw8eJCRI0eyYMECMjMzmTlzpq3japqq+l3FCiAiiSJymoic7775+fYcINq1PF2506l66rR8XwOA54HB1SVHF6WaD2iMqcrGJY2pmV89SBEZATyF8ypU9wUCVFXT/dqRyFycyewWoCewGDhfVTd61LsE+Cdwtar+2+O1FjivbF2Bc+m7YThPs56pqtVO97AepDHevfvuu2RmZtp8SdNU1bsH+Qhwraq2VtX2bptfydHlt0ACsAt4DbhdVTeKSLprPmN5W/cDKcBit7mOS1yvxQCTgXxgN/A7YGhNydEY45vnuORjjz1m45LG4H8PchvQWVVLgh9ScFgP0pjq2bikaaLq3YO8H3hURFoHJh5jTEPjOS7Zt29fcnNzwx2WMWHjb4LMAYYAeSJS6trKRKQ0iLEZY0LMfb7kDz/8QK9evWy+pGmy/E2Qs4FsnFeednVtJ7n+NcY0MuXjkscddxyXXXaZjUuaJsnfBNkKeEBVN6jq9+5bMIMzxoTPSSedxOeff86QIUOYMGGCreNqmhx/E+RLeFk31RjTuJWPS/7lL3+xcUnT5PibIM8GXhCRzSLyb/ctmMEZY8LP4XBw3333sWjRIr7//nsblzRNhr/TPEb6ek1VZwU0oiCxaR7G1F9OTg5Dhw4lJyeHv/3tb4wbNw4RW8jKRDSfP8DRvl6oeKdIFDAKuFxViwMZlTEmsnTt2pXPP/+ckSNHMmHCBNasWWPzJU2jVeMpVlUtBTpi650aY7BxSdN0+DsG+WfgWRHJEJEoEXGUb8EMzhjTMNm4pGkK/E1wLwAjgB9wLlhegnOx8Ihdes4YU382X9I0Zv4myI6urZPbVv7cGNOElY9Lus+XLCwsrPmNxjRwfiVIVc1V1VxgO84e5Ha3MmNME+c5LtmnTx8blzQRz68EKSLJIpINFAE7gEIRmSUiKUGNzhgTMTzHJXv37s2//vWvcIdlTJ35e4r1CaAZcCrOezr2ABJd5cYYU6F8XLJNmzb079+fxx9/3MYlTUTyN0EOAG5U1RxVLVbVHJxzIwcELzRjTKQqH5ccPHgw48ePZ8SIETYuaSKOvwmyCGjjUdYasIUDjDFeJSUl8cYbb/Dggw8yZ84cG5c0Eac20zyWichtInKFiNwGvAfMDF5oxphI53A4uP/++1m4cKGNS5qI42+CnAI8DFwHTHf9+4ir3BhjqjVo0CAblzQRx2eCFJF5bk9vUtV/qOqvVLW7698XtRY/4SKSKiILROSwiOSKyA0+6o0UkTUiUiAiP4rIIyISXdt2jDENS9euXVm1apWNS5qIUV0P8nL5ZZn+GQHY19M451C2BTKBZ0TkFC/1EoFxOMc4zwEuBSbWoR1jTAOTnJxs45ImYvi83ZWILMJ5YU4OcD0w11s9VR1R405EmgH7gFNdV8AiIrOBHao6qYb3TgAuVtXB9WnHbndlTMPyzjvvkJmZSWxsLK+//joXX3xxuEMyTZPPG3FU14P8NfAUsBlQ4Hsfmz+6AqXlSc3lK8Cfnt+FwMYAtGOMaUBsXNI0dD7vB6mqRcArACISo6p/rsd+mgMHPMoOAEnVvUlERgG9gVvq0o6IjAZGA6Snp9cuYmNM0JWPS44cOZLx48dX3F8yISEh3KEZ4/darH8SkRQROVtELnHf/NzPISDZoywZOOjrDSIyFOeVs1eo6u66tKOqM1W1t6r2btPGcxqnMaYh8ByX7Nu3L9u2bQt3WMb4vRbrSOAnYBHwotv2gp/7yQGiReQkt7LT+eXUqef+BgDPA4NVdX1d2zHGRAb3+ZLfffcdvXr14qOPPgp3WKaJ83ce5F+B61S1rap2dNv8ut2Vqh4G3gQeFJFmItIHuAqY7VnX1SudA1yrql/UtR1jTORxH5f81a9+xYwZM2xc0oSNvwkyGni/nvv6Lc6FzncBrwG3q+pGEUkXkUMiUj5IeD+QAix2lR8SkSU1tVPP2IwxDYT7fMlx48YxcuRImy9pwsLnNI9KlZxTLZKAv6hqWdCjCgKb5mFMZCkrK2PKlCk88MADnHnmmSxYsMAutjPB4HOah78JcjtwPM4J+nvcX1PViPiJtQRpTGRyny/5z3/+k4suuijcIZnGpd4Jsp+v11R1RR2DCilLkMZErs2bNzN06FC+/fZbpk+fztixY/lloS9j6sXnD5LPeZDuIiUJGmMap27duvH5558zYsQIxo0bx5o1a3juuedsvqQJKp8JUkTuVdUprscP+qqnqg8EIzBjjHGXnJzMm2++yeTJk/njH//Ixo0bbVzSBFV1V7G2c3vcvprNGGNCwuFw8MADD7Bo0SKbL2mCzq8xyMbAxiCNaVxsXNIESJ0WKzfGmAarfFxy0KBBNl/SBIUlSGNMxCofl/zzn//M7NmzbR1XE1CWII0xEc1zXLJ3797cd999dOjQAYfDQYcOHZgzZ064wzQRyMYgjTGNxubNm7n44ovZuXNnpfLExERmzpxJZmZmmCIzDVjtFwrw91ZWqvphHYMKKUuQxjQN6enpbN++vUp5RkYGW7duDX1ApqGr00IBL3o8PxFQnEvNtXI1+iPg1x09jDEmFH788Uev5TY2aWrL5xik+22tcN6b8UmgpaqmAS2BJ1zlxhjTYPhaOEBVGTRoEB999JHdQsv4xd+LdMYDk1T1CIDr37uBCcEKzBhj6mLKlCkkJiZWKktISOCaa67hiy++4OKLL6Z379689tprlJSUhClKEwn8TZCHgbM9ys4CjgQ2HGOMqZ/MzExmzpxJRkYGIkJGRgbPP/88b7zxBrm5uTz33HMcPnyYG264gc6dO/Poo49SUFAQ7rBNA+Tv3TxuBP4OLAK241xibhBwh6rODmqEAWIX6RhjypWVlbF48WKmTZvGihUrSE5OZvTo0YwdO5b27W0FzSamfivpuJLgOcDXQDLwDXBupCRHY4xx53A4KsYjV69ezcCBA3nsscfo1KkTw4cPZ+3ateEO0TQAtZoHKSIOoK2q7qyxcgNjPUhjTHVyc3OZMWMGzz//PIcOHeKSSy4hKyuLAQMG4HDYmiqNWP16kCLSQkReBYqA71xlQ0RkcmDiM8aY8MrIyODRRx9l+/btPPLII2zevJkrr7ySHj168OKLL1JUVBTuEE2I+ftn0bPAASADOOoq+wwYFoygjDEmXFq0aMHvf/97fvjhB2bPnk1sbCy33HILGRkZTJ48mT179oQ7RBMi/ibIS4GxrlOrCqCq+cBx/u5IRFJFZIGIHBaRXBG5wUe9U0XkPRHZLSJVzv+KyEciUiQih1zbZn9jMMYYf8XGxjJ8+HD+85//8MEHH9CrVy/uv/9+2rdvzx133MF3330X7hBNkPmbIA8Ard0LRCQdqM1Y5NM4e59tgUzgGRE5xUu9EuB14OZq2rpTVZu7tm61iMEYY2pFRLj00ktZvHgxGzZs4De/+Q0vvPACXbt25ZprrmHlypW28EAj5W+CfAF4Q0QuBhwich4wC+ep1xqJSDPgWuB+VT2kqp8AC4EbPeuq6mZVfRHY6GdsxhgTEqeccgovvvgiubm53HPPPXz00Uf07duX888/n/nz51NaWhruEE0A+Zsgp+Ls1T0NxAD/AN4GZvj5/q5AqarmuJV9BXjrQfrjIdcp2JUicpGvSiIyWkS+FJEv8/Pz67grY4yp7Pjjj2fy5Mls376dp556ivz8fH7961/TtWtXnnzySQ4dOhTuEE0A+Jsg26rq46raXVWbqerJqvo4ztOl/miO8zStuwNAkp/vd/cHnAuknwjMBBaJSGdvFVV1pqr2VtXebdq0qcOujDHGt2bNmnHHHXewefNm3njjDY4//njGjh1Leno699xzT5XbbpnI4m+CzPFRvsnP9x/CucCAu2TgoJ/vr6Cqn6vqQVUtVtVZwEpgYG3bMcaYQImKiqoYj/z000+55JJLePjhh8nIyGDUqFFs2LAh3CGaOvA3QVaZSCkiyUCZn+/PAaJF5CS3stMJzDijUs1ET2OMCaXzzjuP+fPn8+233zJmzBhef/11evTowYABA/jggw/sgp4IUm2CFJHtIrINSBCRbe4bzitY3/JnJ6p6GHgTeFBEmolIH+AqoMpSdeIUD8S6nseLSJzrcQsRudxVFi0imcCFwHt+f2JjjAmBzp078+STT7J9+3amTJnCf//7X/r370/Pnj3Jzs7m6NGjNTdiwqrapeZEpB/O3tli4Aq3lxTIU1W/5yCKSCrOi3v647zp8iRVfdU1XWQT0F1Vt4lIB2CLx9tzVbWDiLRxxfI/QCnONWHvV9VlNe3flpozxoRTcXExr776KtOnT2fjxo2kpaUxduxYxowZQ4sWLcIdXlPm8wykv3fzSCy/F2SksgRpjGkIVJX33nuPadOmsXz5cpo3b87NN9/MuHHj6NChQ7jDa4rqtxYr8IqIXFCpRZELRGR+vcIyxpgmRkQqxiPXrl3L0KFDefrpp+ncuTPDhg3jiy++CHeIxsXfBNkP+NSj7DPg4sCGY4wxTUfPnj2ZPXs2W7ZsYeLEiSxdupRzzjmHCy+8kIULF1JW5u91kCYY/E2QRUAzj7LmOJeFM8YYUw/t2rVj6tSp/Pjjjzz22GPk5uZy1VVXcfLJJ/Pss89SWFgY7hCbJH8T5HvAc66pHeVTPJ4ClgYrMGOMaWqSkpIYN24c33//PXPnziU5OZnbb7+d9PR0/vSnP7Fr165wh9ik+Jsgs3BO7N8nIruAvUAKMC5IcRljTJMVHR1dMR65YsUKzjvvPP785z+Tnp7O6NGj+eabb8IdYpPg11WsFZVFjgfaA9tV9eegRRUE9b2KNS8vjy1btlBcXExcXBwdO3akbVt/V9oLfbuNJR5jjNM333zDY489RnZ2NkVFRQwaNIiJEydy4YUXImJrpdRD/aZ5AIhIK5xLup2gqo+ISBrgUNUfAxNjcNUnQebl5ZGTk1NpwNzhcNCxY0dat25dzTurt3v3brZs2eKz3UD90PvbTn5+Pj/88EOVeLp27WpJ0pgGIj8/n7///e889dRT7N69m169ejFx4kSuu+46oqOjwx1eJKr3PMh+wBvAl0AfVU1ylU1U1cEBCzOI6pMgV61aRXFxcYAjihxRUVF06tSJhIQEEhMTiY2Ntb9YjQmzwsJCsrOzefTRR8nJySE9PZ1x48Zx8803k5zsufS1qUa9E+RanMlwuYjsU9WWruXgclU1IroW9UmQK1as8Pla165d6xoSOTm+1oD3r91ArelY3o6/d0h3OBwkJiZWJEz3f+0vWGNCq6ysjHfffZdp06bx73//m+TkZMaMGcPYsWNp165duMOLBPVOkPtUtaXr8V5VTRURB5Cvqq0CF2fwBKMHGRcXx7nnnlvnmILVbjDi6dmzJ4WFhRw5cqTSv0VFRZXqxsbGVkqY5Y/j4+NxOPy9JswYUxerV69m+vTpzJ8/HxFh2LBhZGVlccYZZ4Q7tIas3glyJfCgqr7nliAvA+5R1YsCF2fwBGMMsr5jc8FqN5TxlJWVUVhY6DV5lpRUnibrnjDd/7VTtsYE1tatW5kxYwYvvPAChw4d4pJLLmHixIkMGDDAfteqqneCPBd4B3gX+F8gGxgMXKWqqwMUZFDZVayhj6ekpMRr4iwsLKyUhKOiorwmTjtla0z97N+/n5kzZ/LEE0+wY8cOunfvTlZWFpmZmcTFxYU7vIai3gnSARwPDAcygO3AK5FyBSvYYuUNiapSXFxcKWH6e8rWPXnaX8LG+Ofo0aPMmzeP6dOn89VXX9G2bVt+97vfcdttt9GqVUSMkgVT3ROkiEQBh4AWqhqxl3JagowMpaWlFBUVVel1HjlyhGPHjlXUExHi4+O9XiwUExNjydMYL1SVDz/8kGnTprF06VISEhL4v//7P8aPH0/nzp3DHV641LsH+RVwhar+FMioQskSZOQrKSnxmjgLCwsrXdEbFRXl8yrbqKioinoN7fS2MaG0YcMGHn30UebMmUNJSQlXX301WVlZnH/++eEOLdTqnSD/H3A9MAP4EecNkwFQ1Q8DEGDQWYJsvFSVoqIir4nT86rcuLg4EhISADhw4EClxGqLIpimaOfOnTz11FM888wz7Nu3j/POO4+srCyGDh1a6Q/KRqzeCXKLj5dUVTvVNapQsgTZNJWWlnq9yvbgwYNe68fExHDOOec0lf8YjKlw+PBhXnrpJR577DF++OEHOnXqxPjx4xk1ahTNmnnezKlRqf9Sc5HOEqRxV93iDw6HgxYtWpCamkqrVq2Ij48PYWTGhFdpaSlvvfUW06ZNY9WqVbRs2ZLbb7+dO++8kxNOOCHc4QWDJUhLkMadr0URYmJiaNOmDXv37q24ojYxMZHU1FRSU1NJSUmxBQ9Mk/Hpp58yffp0FixYQExMDJmZmUyYMIFTTz013KEFUu0TpIh8raonux5vx23c0Z2qpvsVgUgq8CJwGbAbuFtVX/VS71RgOtALaKWqUpd2PFmCNO5qWhRBVSksLGTPnj3s3bu3YrwyKiqKli1bViRMm0tmmoLvvvuOxx9/nJdeeokjR44wYMAAsrKyuPTSSxvDFeN1SpB9VfUT1+N+vhpQVd/nqiq39xrO+0/eDPTEuejA+aq60aNeN6AvzuT3lpcE6Vc7nixBGk+1uYr12LFj7N+/vyJhHj16FIDmzZtXJMvk5OTG8J+FMT7t2bOHZ599lieffJK8vDxOP/10srKyGDZsGLGxseEOr67Ce4pVRJoB+4BTVTXHVTYb2KGqk3y8pwvwrXuCrEs75SxBmkBRVQ4fPszevXvZs2cPBQUFgPMmt+XJMjU1lZiYmDBHakxwFBUV8eqrrzJ9+nQ2bdrEiSeeyNixYxk9ejQtWrQId3i1Vace5IP+tKyqD9S4d5EzgE9VNcGtbCLQz9ftsnwkyFq1IyKjgdEA6enpvXJzc/35SMbUSklJCfv27WPv3r3s3bu3Yg3a5OTkimTZvHlz612aRkdVWbp0KdOnT2f58uU0b96cW265hbvuuosOHTqEOzx/1SlBvuT2NB64FlgN5ALpwNnAG6r6mxr3LnIB8E9VPd6t7FYg09di5z4SZK3bKWc9SBMKqsrBgwcrkmX5dJLY2NiKZNmyZUtbY9Y0OmvXruXRRx9l7ty5qCrXXXcdWVlZnHXWWeEOrSb1ngc5F2diesOt7Brg134myDOAlaqa6FaWBVxUhx5krdopZwnShMPRo0crkuXevXspLS1FREhJSalImImJida7NI3Gjz/+yBNPPMFzzz1HQUEBF154IRMnTuTKK69sqFeA1ztBHgBSVbXUrSwK2KuqKX68v3zs8BRV/dZVlg38VMcxSL/bKWcJ0oRbWVkZBQUFFcny8OHDAMTHx1ckyxYtWtgiBaZRKCgo4MUXX+Txxx9n27ZtdOvWjfHjxzNixIiK1awaiHonyDXALFV9wq3sd8AoVT3TrwicvVAFbsF59elivF/FKkAc0AnYCCTgXLGnuDbteLIEaRqaoqKiimS5b98+ysrKKi1SkJqa2tD+IzGm1o4dO8b8+fOZNm0aa9asoXXr1txxxx3ccccdtGnTJtzhQQAS5BnAAiAa2AGcCBwDrlHV//gVgXP+4j+A/sAeYJKqvioi6cAmoLuqbhORDoDn0na5qtqhunZq2r8lSNOQlZWVsX///oorY22RAtPYqCr//ve/mTZtGu+88w7x8fGMGDGCCRMm0K1bt3CGVv9pHiISA5wLpAE7gc9UtaT6dzUcliBNJDly5EhF73L//v22SIFpVL755hseffRRsrOzKS4uZvDgwWRlZXHhhReGYzzelpqzBGkiVWlpaaVpJOVL5NkiBSbS7dq1i7///e88/fTT7N69m969ezNx4kSuvfbaUF7pbQnSEqRpDNwXKShfAg9skQIT2QoLC8nOzmb69Ol8++23ZGRkcNddd3HLLbeQlJQU7N1bgrQEaRojX4sUJCUl0apVK1ukwESUsrIy3nnnHaZNm8bHH39MSkoKo0ePZuzYsbRr1y5Yu7UEaQnSNHa2SIFpTL744gumT5/O/PnzcTgcXH/99WRlZbFx40buvfdetm3bRnp6OlOmTCEzM7M+u7IEaQnSNDXuixTs27ePY8eOISIkJydX9C5tkQLT0G3ZsoUZM2bwwgsvcPjwYRwOR6W78CQmJjJz5sz6JElLkJYgTVOmqhw4cKDKIgVxcXEVydIWKTAN2b59++jUqRP79++v8lpGRgZbt26ta9OWIC1BGvMLb4sUiAgtWrSoSJi2SIFpaBwOB95ylohU6lXWks8EaYMRxjRB8fHxpKWlkZaWVmmRgr179/Ldd98BkJCQQGpqKq1atbJFCkyDkJ6ejre7MqWnpwdlf5YgjWniHA5HxUU84LzkvvzG0D/99BM7duzA4XDQsmXLit6lLVJgwmHKlCmMHj2aI0eOVJQlJiYyZcqUoOzPEqQxppKEhATatWtHu3btKC0tZf/+/RUJc8+ePQA0a9asondpixSYUCm/ECfAV7H6ZGOQxhi/qGrFEnh79uyptEiB+xJ4sbGxYY7UmFqxMUhjTP2ICM2aNaNZs2a0b9+eY8eOVSxSsGfPHvLz8wHnIgXlyTIpKcl6lyZiWYI0xtRJdHQ0bdq0oU2bNqgqhw4dqkiWubm55ObmEhMTU2kJPFukwEQS+2k1xtSbiJCUlERSUhIZGRkcPXq0Uu8yLy8PgJSUlIpk2axZM+tdmgbNxiCNMUGlqhQUFFRMIzl06BDgXKTAfQk8W6TAhImNQRpjwkNESElJISUlhY4dO1JcXFyRLHft2sXOnTsrFikoT5iJiYnhDtsY60EaY8KnrKys0hJ45fPbyhcpKF8CzxYpMEFkPUhjTMNTvgBBy5Yt6dy5M4WFhRXJcufOnZUWKShPmPHx8eEO2zQRliCNMQ1GQkICJ554IieeeGLFIgXlF/p4LlKQmppKcnKy9S5N0ITsFKuIpAIvApcBu4G7VfVVH3XHA38AEoA3gNtVtdj12kfAucAxV/Udqtqtpv3bKVZjIpf7IgV79+7lwIEDqCpRUVGVppGUL1KQl5fHli1bKC4uJi4ujo4dO9K2bdswfwrTQDWIU6xPA0eBtkBP4F0R+UpVN7pXEpHLgUnAJcBPwALgz66ycneq6guhCNoYE37VLVKwd+/eSosUxMbGsnfv3oq7PhQXF5OTkwNgSdLUSkgSpIg0A64FTlXVQ8AnIrIQuJHKiQ9gJPBieeIUkb8Ac7zUM8Y0Ub4WKXBfL9ZdWVkZ33//Pa1atbLFCozfQvWT0hUoVdUct7KvgH5e6p4CvO1Rr62ItFLV8p/8h0TkYWAzcK+qfhSEmI0xEcBzkYIVK1Z4rVdSUsLKlSuJjo4mISGBhIQE4uPjKx4nJCQQExNjixeYCqFKkM2BAx5lB4AkP+qWP04C9uAcm9yE83Tt9cAiEempqt97NiQio4HRELz7hRljGpa4uDiKi4urlEdHR5Oenk5hYSGFhYUcOHCAXbt2VarjcDh8Js+4uDhLnk1MqBLkISDZoywZOOhH3fLHBwFU9XO312aJyG+AgcCTng2p6kxgJjgv0qlT5MaYiNKxY0dycnIq3WHe4XDQpUuXKmOQZWVlFBUVUVhYWPFvYWEhhw8fZs+ePZXuXi8ilZKm52O7mrbxCVWCzAGiReQkVf3WVXY6sNFL3Y2u1153q5fndnrVk1LNVUjGmKalPAn6cxWrw+EgMTHR68o9qkpxcXGV5Fne+ywtLa1UPy4uzmfv05bRi0whSZCqelhE3gQeFJFbcF7FehVwvpfq2cDLIjIH2AncB7wMICItgHOAFTineQwDLgTGBfUDGGMiStu2bet9xWp5j9HbwgSqSklJidfkuXv3bkpKSirVj4mJ8Zk8o6Oj7dRtAxXKy7l+C/wD2IVzLPF2Vd0oIuk4xxS7q+o2VV0qIo8A/+KXeZB/dLURA0wG/gcoBb4Bhqrq5hB+DmNMEycixMbGEhsbS0pKSpXXjx075jV57t+/v8r4aFRUlM/kGRsba8kzjGwtVmOMCaHS0lKv455FRUUUFRVVGvd0OBwVSdPbRUM27hkQDWKhAGOMafKioqIqFj3wpKo+k+e+ffsqXXgEVHvRkI171p8lSGOMaSBEpCLJeVJVjh496jV57tq1i2PHjlWqHxsb6zN5xsTEhOojRTRLkMYYEwFEhLi4OOLi4ry+7uuiob1793L06NFKdX0tlhAfH2/jnm4sQRpjTCMQExNDTEwMycmeU86d457ekmdBQYHfiyWUX9HblJKnJUhjjGnkoqKiaN68Oc2bN6/ymq/FEo4cOVLtYgnexj8b20VDliCNMaYJq8tiCUVFRdUuluAteUbiIvGRF7ExxpiQqOtiCXv27PG5WIK35NlQF4m3BGmMMabW6rpYgrdF4ssXS/CWPD0XiQ/lzbAtQRpjjAm46OjoituQeSorK/OaPH0tEl+eMMvKyjhw4EDIboZtCdIYY0xIORyOGhdL8EyeRUVFHD58uEr9srIytmzZYgnSGGNM4+a+WELLli0rvebrZtje7v8ZCI3rmlxjjDGNlq9FEnyV15clSGOMMRGhY8eOVeZaOhwOOnbsGJT92SlWY4wxEaE2N8MOBEuQxhhjIkYgbobtLzvFaowxxnhhCdIYY4zxwhKkMcYY44UlSGOMMcYLS5DGGGOMF+K+5l1jJiL5QG4AmmoN7A5AO6Fqt64aWjzB1JQ+a1Njx7bxCtSx3a2qA7y90GQSZKCIyJeq2jtS2q2rhhZPMDWlz9rU2LFtvEJxbO0UqzHGGOOFJUhjjDHGC0uQtTczwtqtq4YWTzA1pc/a1NixbbyCfmxtDNIYY4zxwnqQxhhjjBeWII0xxhgvLEF6EJFUEVkgIodFJFdEbvBR71QReU9EdotIjeepRSRORF50tXlQRNaKyBXV1B8vIj+LyAER+YeIBPSOoCLykYgUicgh17Y5XLEEmojcKSJfikixiLzs8dqlIvKNiBwRkX+JSEY17fj1s2BCp6bfIzu+kU1EXhGRnSJSICI5InKL22shP7aWIKt6GjgKtAUygWdE5BQv9UqA14Gb/Ww3GtgO9ANSgPuB10Wkg2dFEbkcmARcCnQAOgF/rs2H8NOdqtrctXXzViGEsQTST8Bk4B/uhSLSGngT53efCnwJzKumHX9/Fkzo+Pw9suPbKDwEdFDVZGAIMFlEeoXt2Kqqba4NaOb6Uru6lc0GHq7mPV2cX2Od9rcOuNZL+avAX92eXwr8HODP+hFwix/1gh5LEI/nZOBlt+ejgU89jnch8D+B+FmwLWzHeR1wrR3fxrUB3YCdwP+G69haD7KyrkCpqua4lX0FBPyvShFp69rfRi8vn+Lar3sMbUWkVYDDeMh1iniliFzko06oYgmFSp9FVQ8D3+P9+IbsZ8HUncfvkR3fRkBE/i4iR4BvcCbIxYTp2FqCrKw5cMCj7ACQFMidiEgMMAeYparf+BFH+eNAxvEHnKdLT8Q5n2iRiHQOUyyhUpvjG5KfBVN3Xn6P7Pg2Aqr6W5zH4QKcp1WLCdOxtQRZ2SEg2aMsGTgYqB2IiANnd/8ocKefcZQ/Dlgcqvq5qh5U1WJVnQWsBAaGI5YQqs3xDfrPgqk7H79HdnwbCVUtVdVPgHbA7YTp2FqCrCwHiBaRk9zKTsf7adBaExEBXsQ5cHytqpb4qLrRtV/3GPJUdU8g4vBBAWkgsQRLpc8iIs2Azng/vkH9WTB1V83vkR3fxieaX45h6I9tuAdiG9oGzAVewznQ2wdn1/wUL/UEiAe640wu8UBcDW0/C6wCmtdQbwDws6vtlsCHBPDiAaAFcLkr5micV3kdBrqFOpYgHcNo12d7CGcvo/xztnEdz2tdZVOBVfX9WbAt5MfX6++RHd/I3oDjgOtxniKNcv0fdRi4KlzHNuxfSkPbcF5C/JbrwGwDbnCVp+Psuqe7nndwJUb3bWs17Wa46hS52infMj3bdtWfAOQBBcBLNSXfWn7GNsBqnKcc9rv+s+nv7XMGO5YgHcM/eTk2f3K99iucg/+FOK/k7eD2vnuAJTX9LNgW1mPr8/fIjm9kb67/l1a4/k8qANYDt7q9HvJja2uxGmOMMV7YGKQxxhjjhSVIY4wxxgtLkMYYY4wXliCNMcYYLyxBGmOMMV5YgjTGGGO8sARpjDHGeGEJ0oSdiGwVkV81hjhEZKOvO6OIyMsiMjmY+68P99jDHUukse+rcbIEaUwAqeopqvqRP3Ub2n+qtYk9Uri+4w7hjsNEJkuQplEQkehwxxBq4nSGl/LTRSQqHDFVJxKPUSTGbALHEqRpUETkf0Rki4hcLyJpIvKGiOS7ysZ61N0qIn8QkXXAYRGJdpVNFJF1InJAROaJSLyrfrXtVRPTKBFZ5Pb8OxF53e35dhHp6RbTr1yPzxCR/4jIQRGZh3OR5fL3zMa57u0iETkkIv/P9VJPb7H70AF4X0QGuLV7AbAcOLmaz/MHEdnhimuziFzqGbtH/UkiMt+jbIaIPOF67PN7reYY/d71OQ+LyIsi0lZElrhi+kBEWnrs/3vXa5tE5GqP9r0eb38/t5d63mJWEeniVsfn6fK6/pyZBijcC9TaZhuwFedCxGfiXFh4EM4/3tYADwCxOG/u/ANwucf7/gu0BxLcyr4A0nAuWPw1cFst2vuVl/g64VxA2QGcAOQCO9xe2wc4PD5LrKveeCAGuA4oASb72p+v2Gv47i4E9gAXA2cD+cCAaup3A7YDaa7nHYDOnvF4PM4AjgDJrudROO/0fm5N32s1x2gVzttVnQjsAv4DnAHE4bxjzB/dYv616ztxAMNwLkB9Qm2+s+o+t4+fR8+YFejiVudlb8eypu/DtsjarAdpGooLgIXASFV9BzgLaKOqD6rqUVX9AXge5+1w3D2hqttVtdCj7CdV3QssAnrWor0qXHUPutrpB7wH7BCR/3E9/1hVyzzedi7OxPi4qpao6nycd1CpibfYq4vt38ANwHzgHWC0qi6t5i2lOJNQdxGJUdWtqvp9DfvIxZnAhrqKLgGOqOoq/PtevR2jJ1U1T1V3AB8Dn6vqWlUtBhbgTJbl+/+n6zspU9V5wLc4/xhwb7+m76y2n9tbzP6o88+ZaXjs/LppKG4DVqjqv1zPM4A0EdnvVicK53+m7rZ7aetnt8dHcPYu/G3PlxXARUAXfrklTz/gPNdzT2k4e5nut8vJ9WM/3mKvyTbgGM7ey9bqKqrqdyIyDuctwU4RkfeACar6Uw37eBX4DZCNMyG/6ir353v1dozy3B4XennevPyJiIzAecu1Dq6i5kBrt/o1fmd1+NzeYvZHfX/OTANiPUjTUNwGpIvIY67n24EtqtrCbUtS1YEe7/P3fm3+tudLeYK8wPV4Bc4E2Q/vCXIncKKIiFtZeh1j90lEOgPLgD/g/A4Xi8gp1b1HVV9V1b78cm/FqX7s6p/ARSLSDriaXxKkP99rnT+niGTg7IHdCbRS1RbABpw3LK+VWn5uz5iPAIluz4/38b76/pyZBsQSpGkoDgIDgAtF5GGc40oFroslEkQkSkROFZGz6th+fdtbgXOcL0FVf8TZIxgAtALWeqn/Gc5e3VjXRR7XUPm0IDh7TZ3q8mHAeTEIzgtypqjqy6r6BjAR54U7XtsVkW4icomIxOG86XAhztOP1VLVfJw3qX0JZwL42vVSoI+Tp2Y4k1W+K/5RwKm1baSun9vNf4EbXJ9vAM4/jLwJ9vdhQsgSpGkwVHU/0B+4AuepsME4x5O2ALuBF4CUOrZdWp/2VDUH553rP3Y9L8B58cVKV9ue9Y8C1wA34byIZxjwpke1h4D7RGS/iEys9YdyXpyTparPuO13DvBbnBe+eBMHPIzz8/8MHIfzbuz+eBXnhSjlvcd6f681UdVNwHScf3DkAT2AlXVoqj6fG+AunJ9zP5CJ82713uIN6vdhQksqD5EYY4wxBqwHaYwxxnhlCdIYY4zxwhKkMcYY44UlSGOMMcYLS5DGGGOMF5YgjTHGGC8sQRpjjDFeWII0xhhjvLAEaYwxxnjx/wHncOjSVJbYbwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 468x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the DI values for different kernel widths\n",
    "fig, ax = plt.subplots(figsize = (6.5, 4))\n",
    "ax.plot(kernel_widths, di_kw[(0, 1)], color = 'k', marker = 'o', label = r\"DI($n_1\\rightarrow n_2$)\")\n",
    "ax.plot(kernel_widths, di_kw[(1, 0)], color = 'silver', marker = 'o', label = r\"DI($n_2\\rightarrow n_1$)\")\n",
    "  \n",
    "# plot clean up \n",
    "ax.set_xlabel(r\"kernel width $\\times$ silverman's rule\")\n",
    "ax.set_ylabel('directed information')\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.set_xticks(kernel_widths)\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "ax.legend()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4295f5d",
   "metadata": {},
   "source": [
    "Very high values would group almost all points as being similar, meaning that no additional information is needed to explain the variance in the effect variable since all information can be obtained from its own values. No knowledge from the past values of the causal variable is needed and we see a drop in directed information. \n",
    "\n",
    "Very low values, on the other hand, would group all points as being dissimilar, and so we see an increase in directed information. In real scenarios, in the absence of ground truth, this could lead to falsely detected causal influences. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26971918",
   "metadata": {},
   "source": [
    "#### Questions: \n",
    "Directed information is unbounded. So, in the absence of ground, what are some ways by which you can find a significant causal effect between units? \n",
    "\n",
    "How about filtering out indirect causal effects like $n_1 \\rightarrow n_2$ and $n_2 \\rightarrow n_3$ gives $n_1 \\rightarrow n_3$?   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ece893e",
   "metadata": {},
   "source": [
    "### References \n",
    "\n",
    "[1] Cover, Thomas M. Elements of information theory. John Wiley & Sons, 1999.\n",
    "\n",
    "[2] Massey, James. \"Causality, feedback and directed information.\" Proc. Int. Symp. Inf. Theory Applic.(ISITA-90). Vol. 2. 1990.\n",
    "\n",
    "[3] Izhikevich, Eugene M. Dynamical systems in neuroscience. MIT press, 2007.\n",
    "\n",
    "[4] Akella, Shailaja, et al. \"Measurable fields-to-spike causality and its dependence on cortical layer and area.\" bioRxiv (2023)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59ae1aae",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}