wb_model_D3_best_fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 3: Fitting functional connectivity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goals:\n",
    "\n",
    "- Load the BOLD data, calculate the experimental functional connectivity and visualize it.\n",
    "   - What can you see in the structure?\n",
    "- Load the structural connectivity, and calculate the predictive power relative to the experimental functional connectivity.\n",
    "  - Note: The experimental and structural connectivity are on different brain regions which are only partially overlapping. To calculate the predictive power, one has to find the common subset of regions and use only those. Use the function `predictive_power()` defined below.\n",
    "  - Don't forget to remove the diagonal and normalize as before.\n",
    "- Run a single stochastic simulation with $G$ in multistable regime and model parameters as before.\n",
    "  - Use noise with $\\sigma$ = 0.015.\n",
    "  - Use both `Monitors.Bold(period=2000)` and `Monitors.TemporalAverage(period=1)`.\n",
    "  - Run for 60 seconds.\n",
    "  - Look at the results:\n",
    "    - Plot the source states of all regions across the time. Can you see some structure?\n",
    "    - Plot the BOLD signal and compare with the source signal.\n",
    "    - Calculate the functional connectivity and visualize. When calculating the FC, ignore the first 20 seconds.\n",
    "  - Tip: write a function `def simulate(conn, G, simlength, dt, Io, w, sigma)` that will setup and run a stochastic simulation with specified parameters and return the simulation results; it will help you structure your code.\n",
    "- And finally, find the best fit! \n",
    "  - Run the stochastic simulations for $G$ in the identified range.\n",
    "  - For each $G$ calculate the predictive power of simulated FC relative to empirical FC.\n",
    "  - Plot the relation of $G$ and predictive power. What's the best value?\n",
    "  - Redo the previous step for the best fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-09-03 15:52:18,459 - WARNING - tvb.simulator.common - psutil module not available: no warnings will be issued when a\n",
      "    simulation may require more memory than available\n",
      "   INFO  log level set to INFO\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mariolavanga/Documents/GitHub/eitn_school_2021_tvb_project/env/lib/python3.7/site-packages/tvb/datatypes/surfaces.py:61: UserWarning: Geodesic distance module is unavailable; some functionality for surfaces will be unavailable.\n",
      "  warnings.warn(msg)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "\n",
    "from tvb.simulator.lab import *\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last part we will fit the network model to the empirical connectivity. That is, we will find an optimal $G$ that produces the simulation output most similar to what was recorded in live mouse. To this end, we will simulate the network model again, this time with noise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experimental functional connectivity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us load the data and calculate and visualize the functional connectivity. We quantify the functional connectivity between two region as a correlation coefficient of their time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = os.path.abspath(\"./dataset\")\n",
    "subj = \"m143\"\n",
    "fmri = np.load(f\"{datadir}/fmri/{subj}/session0.npz\")\n",
    "fc_emp = np.corrcoef(fmri['bold'])\n",
    "reg_emp = fmri['region_labels']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Regions')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAEoCAYAAADG9UOjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABON0lEQVR4nO29eXyc5Xn3+7tm04xGq63F8iobY8xiVkMC2QhkoSSFtOUlpGkLKackbThNT3LakDYnPU1z3jd9z5ut2WlCQtIkhCalkISwhEDICrYBg21svMmbZMmSrV0azXK9f8wIze+SNVrQWBrr+n4+87Guebb7uZ9nbt/P77kWUVU4juMsVAJz3QDHcZy5xAdBx3EWND4IOo6zoPFB0HGcBY0Pgo7jLGh8EHQcZ0Hjg6DjOCWDiNwlIh0ism2C5SIi/yoie0TkeRG5eLJ9+iDoOE4p8U0A1xRY/nsAzsx9bgPw5cl26IOg4zglg6o+CeB4gVWuB/AtzfI7ADUi0lRon6HZbKDjOA4AvPWNce06np72dlueT2wHMJz31Z2qeuc0drEMwKE8+3Duu7aJNvBB0HGcWafzeBpPPbx82tuFm/YOq+rGIjRpQnwQdBynCCjSmpmLAx8BsCLPXp77bkJcE3QcZ9ZRABnotD+zwAMA/iz3lvjVAHpUdcJHYcBngo7jlBAi8j0AVwKoE5HDAP4RQBgAVPUrAB4EcC2APQAGAbxnsn36IOg4TlHIYPYfh1X1XZMsVwDvn84+fRB0HGfWUSjSJZKr1AdBx3GKwixpfEXHB0HHcWYdBZAukUHQ3w47M0ZE+kVkTYHlXxGR/+cVHuPKnADulBhz9HZ42vhM8DRDRFoANALId9f/pqrePtvHUtWKSZa/b7aPaRERRfYt4OgvKKWqNbllVQA+DuAPASwC0A7gRwA+oaqdxW7bQkaBktEEfSZ4evL7qlqR95n1AXAyRCR4Cg93Qd651uSOHwHwGIBzkQ24rwJwOYAuAJedwrYtWDIz+MwFPgguIETkFhH5tYh8RkS6RWSfiFyR+/5QLkXRzXnrfzP3SPuoiPSJyC9EZFXechWRtXnrfllEHhSRAQBvzH33ibz1rxeR50SkV0T2isg1ue/fIyIv5o6xT0TeOwun+2cAVgL4A1XdoaoZVe1Q1X9W1QdnYf9OARSK9Aw+c4EPgguPVwF4HsBiAN8FcA+ASwGsBfAnAL4gIvmPue8G8M8A6gA8B+A7Bfb9xwD+PwCVAH6Vv0BELgPwLQB/C6AGwOsBtOQWdwB4O7KztfcA+MxU8sBNwpsAPKSq/a9wP85MUCA9g89c4IPg6cl/5WZ6o5+/yFu2X1W/oappAN9HNs7y46qaUNVHAIwgOyCO8hNVfVJVEwD+AcDlIpIfm5nP/ar669ysa9gsuxXAXar6aG75EVXdCQCq+hNV3ZtLf/QLAI8AeN00zveZvHP919x3i1Egc4hTXLJhc6XxOOwvRk5P3qGqP5tgWXve30MAoKr2u/yZ4MtpiVS1X0SOA1gKTlc0bt2TsALZkKZxiMjvIRv+tA7Z/5jLAbxQYF+Wi1V1j/muC0DBPHJOMRGkIXPdiCnhM0FnMl6e9eUekxcBaJ1g3UIPNIcAnGG/FJEyAD8E8L8ANOZebDwIvOJf0M8AvFVE4q9wP84MUAAZnf5nLvBB0JmMa0Xktbm3rf8M4HeqWmjGNxFfB/AeEblaRAIiskxE1gOIACgDcAxAKjcrfMsstPvbyA68PxSR9bljLhaRvxeRa2dh/84kpHOzwel85gIfBE9PfpRzZB793PcK9vVdZB9VjwO4BNmXJ9NGVZ9G7qUHgB4AvwCwSlX7APw1gHsBnED25coDr6C9o8dLIPtyZCeARwH0Anga2Rc8T73S/TuFyUaMlMYgKFoiDo3OqUdEvgngsKp+dK7b4pQW554f0e/+uHHa21246vAWzyztOE7JMzoTLAX8cdhxnAWNzwSdCVHVW+a6DU5pohCkS2SONSetFJFrRGRXrkr8HXPRBsdxiktGZdqfueCUzwRzgfVfBPBmZGuCbhKRB1R1x6lui+M4xaGUNMG5eBy+DMAeVd0HACJyD7JV4yccBOsWBbV5Rfhle2+ikpYP9UbJlklqPmfiJkAnwxcr3Mt2xuRDyVSZ7Yd4Qq1l5o17mvdn26f2KpjdBxNmfdOewIhdn3eQinP7AilzOLu/JLdfA1JweSZSeHmyrrAHggyY/jP9Yc83NMB2qoL3LylzPft4+Ugtbx8YNte7zBzPRB+LOR2ZJOg1GTftGTDtqTT9Z69PxOzQ3B8hE6A42fULJvgGzITMA2GeOTx4AsmRgRmMZoK0lsbj8FwMgierEP+qQhs0rwjj6YfHwlVv3Hc1Ld/28Flkh/vMDsw92r9xiOzMCP/Kmh7hbklU8z3QfzX/CgM7OK3ecDOPWtLP+4sc5+ON1PFNGTSDQkULH3+kmkxUtfCvorKFz69jIwdNlB/j9YcW8fEqj/CvMGkG0fL2JNn9S/lXGuvk7Y++h3+lah57Ipu5/xK1fMFS5j+dxVu4PZ2v5uNFOrm/l/2C27v/nWSiYie3f6CZr0fjrwsP8uF+bl/ADIpHL+X9NzzD7TlyJbc31mHut1Xm/hjk8699kUzEj/L+U+V8v1Xs7SE70cD3RyYytv9nf/WvmAnZ2GEfBF8RInIbgNsAYOWyedtMx3EmwB+HJ2ZKFeJV9U4AdwJA7fp6zZ/93bvmMVr3gj6eCZZ38P/MkV62B1by807VQf4fKzTMM4tUjP8nvXTFQbI37zyH7GBXmOxMhGcGiSX8PzUi5vHV5CNdtJNnAh2XcPtrnz9BtvTyTLVyCcsFlTuPc/suXky2fXxKlrM9bmbRys/j0f1dZKdTNWTX1fJUPd3BM5FAko83nObrEx7i/rQzv4y5q0MD5vlSeAWxj5d9fLwRkz87WcHLa/aZ62v6J2jkCg3Zx1NeHuSJPGJtvL/QIC9PR+3zuem/WvNkYZZHW/l6aGSsfwIjM8vtouqPw4XYBOBMEVmN7OB3E7KhUo7jnEZkfCZ4clQ1JSK3A3gYQBDZHHPbT3U7HMcpHtm3wz4TnJBcenNPce44py3+ODyrDPVG6Q2w1QC3/t2XyL7gX/6K7JDRkAIj5u2kcaEIDbAGd+z3WRcZvJ81wDLrwmI0H6sJVu1kzXCowbxt7OPty46whhdfWkd27/oa3n6giuxElbkZ1y8iM9rN51v+Emt6wbWsGZbv4eU9FzWQnQnXk131BN9mwS7WNMuGWLMLpFkDq3ueNdSyQ6yBQvn4lQdYNAv28tvpyu3cf5Zlj5vjdbBIp2VGozPXR8tZg80mrhkj+qOnyY6tuoLs2pf4huq8wPjsmKfMyoPcf6F+bn/dZtb8Ms/vZPt1F5GtwbEDqHWfmSL+dthxnAVPeo4iQKaLD4KO48w6pRQ77IOg4zhFIeOa4OwhaY4CsX6AVgPc+mHWCM/89l+SXbOL99+7mu3yDtZ8xIRVjQt7M2FdqXITthQ1foBR3mC8Xxvbx17DGpsNu2q9itePHWHNMbmBd5gc4uUNj7NdVst+ewNN3MDQUA3ZPav5Zq/fyh00Um394nj9vuXcH+P6M8oRF5F61gCPXWIiOtIxsit3sqNe2oShDTYnzXJeIXqc+6fiCJ9fYu1ysq3f4dBi46f32gvJtmPFQJNpoHEDHB9WZ/wCG1hDjBgNM7qKiwX215jz23Tg5b9l2Dg5ThF/O+w4zoJGISWjCZbGUO04jlMkfCboOE5RcBeZ2SZPF7GxwNYP0GqAu//0y2S/+V3vIbvnTPbrOnwTa0RLfswai01FFD3Gfl3tKCc7kDSa2bPsd3b01axh1e7m4x8/mzWb6n0mS02SNaSy49wfQx28fxurHBkoHB9qj2djh21qqfJtXJa4Os6aWXpc6i2juSZsKijef7SL+yd6jK9PJmj2d4z9+CI9xk9yC/8M4u0susVb2M9ODh4lu+zslWQHktyfJ87i1G8h47e4eBv3Z2DEiM7C16/iCOt0wUEjEhrCh9mvU43GWtbJ11eH89o3w2LAqnBnacdxFjLiscOO4yxcFD4TdBxngeMuMrNIJp6hbNA2H6CNBbZ+gFYDfPR73yD77K+yn+GaT7Hmcuga1mRScdZJyk6wBjjUaDMjs8aTivH+Aia9YPda1gDrtrJm03Uea5irHjT5/Fo4tjY0zH518aOTpHc3mtSJczihXryNG7xoB2tSw+uWkH3sYv4xJCtZM6vmUFYkTXuinSZ9vpGpyrpNfsFebn9y3TKyuzdweyv28s+g9wwTuxznVN6pi2t4+1buj3Af21Um36AkuX2H3sIaXfygySy9krfvPovXr9vK90vVjm6yBzY08f63s6aZquTtE68bi83P/MLELU8RxdwVTpouJTEIOo5TevhM0HGcBYvCw+Ycx1nQiNcYmVUyQhXhbE0Qmw/QxgJbP0CrAb74Xo413niE/QxtLLDFxrpaKcRWjys3mtzAMhNb283bd57Puky8jTW1YxewRlS5iGONu9fz/vqaWQOKH+b2DDeyxmn9MC2ZsNHwdreTLVdyrGrlfhM7bWQn259Djbz/eDv3Z/c6Xj+Q5Nu6bh9rpOFu7p8MdwdCpsJkaMhomLs5FnuokTXe4CBrgr1rjGa9m/s31sbLyzv5eCPVfL42NtkGF6dq+H6P9BjR2dQYie3i65U6MFYMUtQUNJkiPhN0HGfB4zNBx3EWLKriM0HHcRY27iw9i4R7BU2PjDXV1gW2NUFsPkAbC2z9AK0GuPnjHGt8zXV/QvZQk9F0DrNG1L+G/erCJja3rJ11Fg1yTRAbu3rsfBatwoO8v2in0ZxeYA0s3sbt7V3FmlHVPo5lDvVz/8QOcXtkiP0WR1ZyLC6UNarGTXx9hhbz9Ykf5f0n4/zj6VljNESjQdbu4MMPLOXli1oOkV2zkzXB6hY+/0yQjx/d28HL24+RHT/I1w9pPt9Yew3Z4Va+PhVtrCmWt3JssQb4epW38/WJvHAAhUh3cuxwptzEtldxbHOwaux8pP+V1BiZ/cdhEbkGwOeQrVT5NVX9pFm+EsDdAGpy69yRK+w2ISUxCDqOU2rMfrU5EQkC+CKANwM4DGCTiDygqvn/DX4UwL2q+mUROQfZqpbNhfbrg6DjOLNO9u3wrM8ELwOwR1X3AYCI3APgegD5g6ACGJ3KVgPglEYnwQdBx3HmE3UisjnPvlNV78z9vQxAvrZxGMCrzPb/L4BHROT/BBAH8KbJDlgSg2AmCCTy6lSkYqwR2brAtiaIzQc4LhbY+AFaDfChB/6d7PM+x36GcjbHlsoVrPnUxFjj6biPY2sH2I0O/StZA1zxKGuIg018PsP1fL4dV7BGl6g1+fVM+rlkOffH4m28fKiRNaPKPb28fdzUIFnMGlmy3Ph19pt8e+v4fGPHeHnNHtbYKnZ3kz14OddFXvljzh+IM9lxtH8l98dIFZ9/pNf43ZVz7G1FkO+/4WY+ftDUUQ6OmPyIff28fY3xE2xnO2FqtPStYI2wbC07Si5+wWjUV51Jds3TPDnqu4Dvx8qnD44ZgzN/pJ1h2Fynqm6c8UGBdwH4pqp+SkQuB/BtETlPVSdMmlkSg6DjOKVFkRIoHAGQP2VYnvsun1sBXAMAqvpbEYkCqAPQgQkojXfYjuOUHBkEpv2ZhE0AzhSR1SISAXATgAfMOgcBXA0AInI2gCiAYyiAzwQdx5l1sun1Z3cmqKopEbkdwMPIur/cparbReTjADar6gMAPgTg30Tk/0L2Jcktqlow7rMkBsFMVQb9V4/pHJeuOEjLB+8/h2xbF9jWBLH5AC3WD9BqgNs+wLHG1+66luwXW1hD6jvKmlrEaIDJGpYrAkbTbLuc22Pz59l8frE21oRar2TNsuoA98fAEta4EnUcizxUx/9Dp2K8P7WxqAe5fSdu4PaMJPi2q9jCmlzX+ebHY/bft7yO7J4N7Afac1Yt2at+yn6Nw2tZox3uYU0yOMTnO9DN7R1YwtfX1gGO9PP6XRdzf9dvZg2ubw1vb/0CB7hEC9XgBoAhown3NVu/UFP3uIWvXyDJ1yvZ3DjWFtM306EY+QRzPn8Pmu8+lvf3DgCvmc4+S2IQdByntMhqgqWhtvkg6DhOUSiVBApFG6pF5C4R6RCRbXnfLRKRR0Vkd+7f2kL7cBynNBl1lp7uZy4o5kzwmwC+AOBbed/dAeAxVf2kiNyRsz886Z6GAgjsGIvH3byTNcAylnzG5aOzdYFtTRC7vo0Ftn6AVgN88CwOTVz3W45FtnV1K4xmJhn+vyjDkhyix20dWpvvjpcH27vJDvez3571W4seN5qkWW7r+Ea7WARTs3xoBWugmZf4Nqs+DF5uZKeqfWynYrz/cD+3L9DPF7DmRd4+0sp+jfFtXHMlzRLcuJom4QFTU6aH+3ukgq9fzMT2Vu41NWU6e/iA4Psr3sHXI1FrY6d56zLjFlm1h/0Qk+V8PQIJvn4ZUwearueMx6XSeRwuWitV9UkA5vLgemSDm5H79x3FOr7jOHNLJld7eDqfueBUa4KNqtqW+/sogMaJVhSR2wDcBgChan9qdpxSohguMsVizl6MqKqK2OKJtPxOAHcCQHTpisI+LY7jzDtK5XH4VA+C7SLSpKptItKEAqEs+WiZYrh5TNcLdrEooiH+H8fGArejcF1g+x+WzQdoY4GtH6DVAF+6hfMRXrLlRj5+gv3cBleyRhMwfmoNT3MDe5uN356pMRLdsIrsvtWsMQ0u4cueLuP+qDhkbgt7L5v+OnYhr9/wLPvtWT/IwTTvMMjp/DDIbnQImeXR40YTNHWnR0wsLkKsqfWv4/aVtfL91HsG77/uWd5dWTdvn45w//et4tju/tWsIWbqTKy5kXyHa7h/Rqqthszrh01NlHTM+D0ajdfmexys5/6JFoyvmBqlVHf4VA/VDwC4Off3zQDuP8XHdxznFLHgNUER+R6AK5FNjXMYwD8C+CSAe0XkVgAHANw48R4cxylVipRPsCgUbRBU1XdNsOjqYh3TcZz5g2uCs0laIHnxmJmI8WMztkaN31vSaGhVxq/O1AW2NUFsPkAbC2z9AK0GuOWSe8k+Y/f7yF60rJv3389+ZRWHeP9Hr+bLVv9L1oAipubKUKOpmcElJ8b54VUdMLG4q83+T7AfXJWpI2xrvgQHuL1Ww4sZv7xIDy8fVxd4mNcPmZx3YZMPMB1nzU6GTR1f64aZ4uMPNLEdSLHmlzZ+dvaprmoXHy9Zy46Jts6xzbcYSNnjmzrVw/b3YDTXhNEEjaZY0Wpq9nSPibCSmjANX2Hm0Pl5upTGUO04jlMkSmMm6DhOSVGsanPFwAdBx3GKQqk8DpfEIChpIHJ8THdKLGHNqmoni0apKGtU9c+yo1kqZuq8HjWxoaYusK0JYvMB2lhg6wdoNcC9N32F7HVP/hnZ0RhrbmETqxvbz3Vz4+3cHzYWuno3qx6VBzmWuus81rgs0ROsC43UsMZWtZ/7d6SWl5cf4R9D0pTprd9kYmmNH9uJDexXZzWwMnbjRHmnEfmefoHM6Fuv4OVG9ipvK6zxWTt23MRSB4xml+D+D3Vz/1fv4/7KmF9l3CSQDw2Z2OZBthM1vINxsd5RXh47YBIUpkz/zQB/O+w4zoLHB0HHcRYspRQx4oOg4zhFwV+MzCIaAkbq8nSKCIs4Qw3GT8qc1dFXm3xuLKFhYJnNp8aila0LbGNhbT5AGwts/QCtBvjS679F9t+0cdnVrSsvJDvIkiH6VphY4LCJpeZQaACsQdl8eTa/XMbEZocH+fwGlpu6vT28PMFlkJGs4v4bWMUNTMaNhmY0QDGa4XC9zffHoii3Dhius/kTC9c0sZpjaMhopBUm35+5/0Yqzf21vrLgchsrbfM1SsbUBDF+nvGjpu7xsIlN72O/15GlfL+H+sd+IBqYoRed+uOw4zgLGH8x4jjOgscHQcdxFiz+YmS2yXB8b0pYgwn3GU2FS4SgdjeLgN1rTb61bl4/3s4aSv9KXt/WBbY1QWw+QBsLbP0ArQb42abNZF/TeRa35xzWhBb/wOSfq+T+UaNpRU+wH1hPs4ntNbGmVrMa51fXamKrV3Fs7Pg60Cb/n6nRUdbDB0hHzfrm/CRtNLeAiZU1Gl95m4mt5eYjZUREG3tr8wmilu+PQdO+WJeNBebN+0xdYBy3/c+LrUZr/QaDCZO/cSlfj1jQ1pHm9o5U58XpvzTzgUx9EHQcZyHjb4cdx1mwqL8ddhxnoeOPw7NIMAFUtIx16KKdLDKVHeHKnsdew7G1x89mzaZuK8dudp7PsbPHzuf1VzzKscRtl3PNElsX2NYEsfkAbSyw9QO0GuBDP/oO2deufz2394/OJdtqmqkYa0DDpo5txVFuv61jvPQn7WRrP4uumRUcW73413w90mXLyI52s2bVvYbbM67OsPFjTFaY/HpGorOxtHjVBjKrWvj4kT4+33AP7zDcYWJrj3IRjnDQaIA9XOd4/3+/jOwzP3+AbMksJzu+l7cfWsWioBqNs2K7uT6mpkqk3WjE+w+RnbmC75/YL196+e9gvxFMp4y/GHEcZ4FTKjNBT6rqOM6CxmeCjuPMOh4xMstoEBjJSynXcQlrePGlnL/PakjV+2z+PPabirexRhQeNH5WTXy8su7Cfmi2LrCtCWLzAdpYYOsHaDXAB3c+SfY113Od4cD+Vm7PlWvJrtrLml7f6jjZlS2dZB/6Q9aslv28m+xMjM8vuYjPr38F989QA2tUjZv4+iSNn91wjdHcOvn6qHB/29jl4B5OyDe0kTXXQMrkL1xrNV92HFy01dQpPovzHdp8gnXP8f6PvZmvl61HNFhfS3bSxhaP00i5Dnb1jm6yEw18faNGgwwNsgaaObt5rG3bCueanBAdlxZy3lISg6DjOKWH+wk6jrNgUZTOixEfBB3HKQLuIjOrBEbYt6v2eQ6G7V1fQ3brVby9JDm4d9WDLMIdu4CXRzv54g3Xs71oh6nTavzqUmZ/ti6wrQli8wHaWGDrB2g1wIfu/zbZl2/9I7KHR7iGx3vWPUH2V/7nH5Ddv4H9/qoO8Pn1r2a/teNns2a3/GfsV7fyp9zfwX7WAAebWVPrW879sWgHrx/pYE1zuIY1tFhLN9lSxe2t3sv9bzXIuhf4eLaOrxq/wLITpoaHyf83sISv/6Kt7AfYY/ILlvVyf3fHeHvrJ2nrRFsipqYJTOxwXzNroDXPjPlBzrjuMFwTdBxngeOPw47jLFhUfRB0HGeB45rgLBJMZFDZMlZ4QXpZEwoPcI2E2BHWUMpMfrZoC2uKlcavreYFXt5xBRfJiLXx8YPt3bz/DazZRQZY47F1gW1NEJsvz8YCWz9AqwH+9oIfkr3+V39K9r/tf23B/YcGub2RtK3pwZqnjd2VtKnBkjAJBdvZDzFSy5pUrNNobm2soSUbuCZJvMMk6Gvn2F5tauD9dXE8bCbCfoDWby7QyZpmoIPvj2iCNc1MzN5/RpMzfoGLnjzIX5jY39gB9vNL1nN/BUz/2vZlBjj2PXMW35/BEb6+2taRd7DCemMhXBN0HGdB44/DjuMsWBRSMoNg0RIoiMgKEXlcRHaIyHYR+UDu+0Ui8qiI7M79WzvZvhzHKT10Bp+5oJgzwRSAD6nqMyJSCWCLiDwK4BYAj6nqJ0XkDgB3APhwwR3FA+jYOKaLVC7h2N9EFY/lyQ2s2Q11GM1nmDWi7vV8vHgbay6JWv4frfVK1oDC/axJ9q02dZEbuX3Vu02NB1MX2NYEsfkAbSyw9QO0GuDO17If4V8ceg3Z2ys59jp6jDWzo5ezH1vU1Myw0VGJOu7v4+ewhljeXkO2jZ3t2mBrkHD7Kg+xTtW+kTW4isZzyG54lPPnHbmONeCwSRfYeZ7RHNv451neYa73ctbwrEZqY38Xb2ONMLWCY3mjHbz8yBtYEwz38/7LO/h6lJv8jfl1hAEg0MuFjcuOmyI5a/Las8ssmyol9Ha4aDNBVW1T1Wdyf/cBeBHAMgDXA7g7t9rdAN5RrDY4jjOHlMhU8JTkExSRZgAXAXgKQKOqtuUWHQXQOME2t4nIZhHZnBoaONkqjuPMY1Rl2p/JEJFrRGSXiOzJPUmebJ0b82S47062z6K/GBGRCgA/BPA3qtoreWmnVFVF5KTjv6reCeBOAChvXFEiL9sdxykWIhIE8EUAbwZwGMAmEXlAVXfkrXMmgI8AeI2qnhCRhpPvbYyiDoIiEkZ2APyOqv5n7ut2EWlS1TYRaQLQMfEesgRSQPmxMd2jcifXsMB69uNLDrFGhAhrJvGjPKb2NZv8cLZurnFDs7G01s9qcAl3a7SLt688aPzGYGKXTV1gWxPE5gO0scDWD9BqgP+24tdkvzZ9Htnd61iDqjjC7bF1gofrjB9mO/ulNRi/w1Ava442H5+tI237I3aI/QYrli0mu/ZFPr6OsJ9f7UvGb9FMQCpbTd1fc76xdtbUMhHWkG3d3wE198MRbn//OpM/sIrvB5u/Ml1m6ih3GD9G6zdo64QYB77EIuPXuCvPDzVlbv5pUAQ/wcsA7FHVfQAgIvcgK6/tyFvnLwB8UTWrrKvq5OPLrDczh2SnfF8H8KKqfjpv0QMAbs79fTOA+4vVBsdx5obRVFozeByuG5XBcp/b8na7DED+W67Due/yWQdgnYj8WkR+JyLXTNbWYs4EXwPgTwG8ICLP5b77ewCfBHCviNwK4ACAG4vYBsdx5gIFMLO3w52quvEVHDkE4EwAVwJYDuBJEdmgqt2FNigKqvorjHvQeJmri3Vcx3HmB0V4HD4CYEWevTz3XT6HATylqkkA+0XkJWQHxU0T7bQkIkYyQWBo0diTe+Zi1oCi3abu7+OscUQGWKMZMX5b8cN8tar2seaTLGe/t4ElrNFFj/P+02W8P1tHt+s8rttga6L0NPNlsXWBbU0Qmw/QxgJbP0CrAf7qX79K9oX/46/ITvHhUHmAz3fJ77j9gW52ZBs2+QdhYo/DJra6+YcsonZdyu0fXsr7y5i7OFHH/RvazppY/zKjsbZwfw2aGiiLX+DzkaSpU9zL2ycruUFiJMjec1nDHq5mVSraw/0b4tsRGSN59y3n8637FceWH3v9Ut5fgq9XpMfEtqfy7FcykM3+ILgJwJkishrZwe8mAH9s1vkvAO8C8A0RqUP28XhfoZ2WxCDoOE6pMfthc6qaEpHbATwMIAjgLlXdLiIfB7BZVR/ILXuLiOwAkAbwt6raNfFefRB0HKdYFMGxTVUfBPCg+e5jeX8rgA/mPlPCB0HHcWafEgqbK4lBMJBUVB4Z011sXdfyl3i2W1ZrRCy7vxHjh9fIfl6hfva7WryNt0/UmXx6xk+w4hB362Q1IDIRPp+A0WxsDRNbF9jWBLH5AG0ssPUDtBrgcx/5EtnXnvMGsnU517lN1bBfZaqF8+OVrWZNL9Rr8+uxJja0qoZsWwc6tvMo2Zkwa16BpKlLvJKXNz3G+QYz+7i9FUvYvzZTw7HEma0vkh2t51jkWJQ1upFmE6t8lGO9Y0tNPkJTA6RqBwc3y5DpvxG+v1IruL9rd5n8l8ZPc2gFx0Knj431j+rM/QTnLCPCNCmJQdBxnFLEZ4KO4yxkfCboOM6C5nQaBHMJUb8BoA/A15DNCHOHqj5SxLa9jAYEyfiYTpIs52l2cC37DQ408WlV72MN5cQ5rPGEhvhqxQ6xDjLUyH5pQ3Ws2WRMnVkbjNizmh27oidYs8qEbI0Rtpf+pJ3sQ3/I+edsLLOtCWLzAdpYYOsHaDXAB3f8guzXvf+9ZAeSpv8u4jrJ/ctYQ5UlbMfbWIO1/REaMrG4G1jj61vJ17tmt9nftp1kd9x+BdnBK1hDS8VN7LLJnxhrupSPb+ok26fAila+n0ZqzPGi5njHef2uS/j+tvkXrR9s5cM7yO5+O1+P6t3mfEwsduryC8aMrb/BjJh5xMgpZ6qxw3+uqr0A3gKgFtlwuE8WrVWO45Q82bKb0/vMBVN9HB4d0q8F8O2cg2JpDPOO48wNp9PjMLKp8R8BsBrAR3Lp8jOTbOM4zkKmRB6HpzoI3grgQgD7VHVQRBYDeE/RWmUIJBXl7WO+UKlyju0s38N+gqGhGrLt+vG2wn571g+rco/RTGLs1xXtMr5U5tpHTrBGNVJjYmcHzfbmvxftZz+vZT/vJrvfxObausBW07L58WwssPUDtBrgL7/IscZrn7iF7DP+eDvZNaENfPxa9qOz/RUe4P6wfqHRl1gjLd9n8j8OmxoeAVPH2GiyZSZWN19/BsZr0Pa3bTXe8jb2w0tHC//MyjpN3eduzocY7ubrObicRdzutSbWfH0z2ZUHuT0Nn2e/yOd+wLHkTZ/K0wHVBC6fhkxpEFTVjIi0AzhHRPyNsuM4k3LynPHzj6m+Hf4XAO9ENoPr6KtFBfBkkdrlOE4pM5c1NKfJVGd17wBwlqravPCO4zgnQU47TXAfgDCAORkEMxFB/9IxXaSilTW2nos41rNnNWs6dlq+aAdrMJmw8dNbaWqWxLmb1LwYV+MneOxCEzu83+Sv2886y8ByzldY3soaTmYFxwZnYrz/42cXrntrNTdbE8TmA7SxwNYP0GqAe678JtlvW/JW3p/pn5DR/IYaWCMcl+/vQOH41aHVXKMjtp8LN4eaTEFD0x8DTXy8mPULPMgdKsYPs6yLr5dkePmQqZNdtYnzgA6cZ2Kxzf021MB+pnZsqTzE7Q10s4YcMLHZez7LdZmHrmcNctY4zWaCgwCeE5HHkDcQqupfF6VVjuOUPqfZIPhA7uM4jjM1TqdBUFXvFpEIsqmqAWBXLoe/4zjOeEoobG6qb4evBHA3gBZkFZUVInKzqp6St8OBpCKW50sV3c9+gZkw52ur38qxseXbuObC8DrW2KK72e/Mxu+EFnO+tdhBo/msYD+9hmf5/4eQqaExUst+X5Ee1rz6TN3jxb/mOsvJRXy+y39m8s2lWSNK1LHmaOsC25ogNh+gjQW2foBWA/zJMw+Tfe3ruQYKwkYzbetmeyuvnmk3+f/WrmTbxBrb/ICaNPkhn+D10401vH2E2xc6zhobutlvFGlTRKSBY30rf2bas4SvX/luzg+JPj5e4g1ryLYabfl9T3FzUJjqw21k1/yc82k2/m7sfo7cPPOqvKeViwyATwF4i6ruAgARWQfgewAuKVbDHMcpcUpkEJzqMB8eHQABQFVfQvZtseM4Tkkz1ZngZhH5GoB/z9nvBrC5OE1yHOd04HR7HP5LAO8HMOoS80sAX5p49dklWac4+p4xX6x0qoaWVz3BpzFSzZpPdZzz7x272PgRXrmC7MZNrKoky3n9EzewZpN5iY+frGFNLjjAy8uPcPsS7JY4rk5tumwZ2f0rePuVP2XNSxK8g+PnsAbZYGqQ2LrAtiaIzQdoY4GtH6DVAB988j6y37j9erJPPMDnN7CMfz3pOPv5xdrYry+5ga+Hvo5VmlU/ZffW3bdweyu3sp/iwHK+fvWbWVMNJm1+P1vHmv1AD93AsbnLfsl+hYeu5uNHeozfapXRANv4fhx87+VkL9rB++9Zwxpz3UN7yU6tYT/FjneOadDJVu7raXE6vRjJRYp8OvdxHMcpzOkSNici96rqjSLyAk5ySqp6ftFa5jhOaXM6DIIAPpD79+3FbojjOKcXp4UmqKptuX8PnJrmFGrLmL5QV8t+ccEu1lSCCdZM0qaub7LS1GQwsb1Di9mO9Jt8fAlTw+Qwt3UwzccPjJjjs9shklVGQ0yY/HndvHyogdsX7Dch3e3sd1beXkN2yNSdhck/aOsC25ogNh+gjQW2foBWA3z83PvJvuJr7+P9V5n+S7Ed5suPRCtrdpF+7r9IK9f5veZcjt39+YGLydYGPv+RqqixTXv6TY2VDhPLbbonUcuOFaFBbu9QI1/vcK/tDz6ezWcYSPABq1vM9a5mDVgjM/cFLMjpMAiOIiJ9GH9KPci+If6Qqu6b7YY5juOcCqb6dvizAA4D+C6yESM3ATgDwDMA7gJwZRHa5jhOKVMiM8GpzoOvU9Wvqmqfqvaq6p0A3qqq30e2+tw4RCQqIk+LyFYR2S4i/5T7frWIPCUie0Tk+7mYZMdxTiNEZ/aZC6acSktEbgTwg5x9A4BRoWGipicAXKWq/SISBvArEfkpgA8C+Iyq3iMiX0G2fsmXCx1cBgKIbB6rFZzu4BoLZUMm9na5za9n/Aa5DC0ypuRF/Cjv78Q61nAqtrAGlTGxM0FTliHWyV1Uv4k1qoFVXAfZ1gDpXsPn07iJNavBZq55EqnlWFBbp7b3LF4/bGKbYfLP2brANh+fzQdoY4GtH6DVAH/zma+QvfqB28iu3mH9QPn4oQFuUNU+1sgGz2RHzO3/xH5+dWV8/okDJv9fC59/WRvHDo+YutThTo7NDg2Z/T3HsepDteynt+xxvoH6mvl+G1eTZdBogiN8Pv2r+H6oMjVLZMTUmInk3dCvpKhkifgJTnUm+G5kaw13AGjP/f0nIhIDcPvJNtAso5H54dxHAVyFscH0bmSzVjuOc7qhM/jMAVN1lt4H4PcnWPyribYTkSCALQDWAvgigL0AulV1dKp1GMCyCba9DcBtABCuOukTt+M485hScZGZ0kxQRNaJyGMisi1nny8iH51sO1VNq+qFAJYDuAzA+qk2TFXvVNWNqroxGItPvoHjOPOL02kmCODfAPwtgK8CgKo+LyLfBfCJqWysqt0i8jiAywHUiEgoNxtcDuBI4a0BDQGJ2rEeshpfIM2amZpwR0lw7yYrbY0QXt/WnY0dY82k63yrQfH2g5yucFwsqM1XaI9n6+CmYtbPkRvct5wvY6yTl3dt4O3DfWw3/5DzMw6tqiHb5uuzdYHH1QQx+QBtLLD1s7Ma4P7r7iT7wh1/RXbUaKzd53F/DS/i9oip4zy4xGierbxCxvwqbA2awTU1ZJfv6+YNkqaGjdmfDrBmmKwwfoKmJom9fvFWtmMdnL8ycIIdKeOtrAFmtr7I9hvZT5L8PGeqCc7hi47pMlVNsFxVnzbfFax+IyL1IlKT+zsG4M0AXgTwOLIvVgDgZgD3n3QHjuOUNqfZTLBTRM5ArpkicgOAtsKboAnA3TldMADgXlX9sYjsAHCPiHwCwLMAvj6zpjuOM68pkZngVAfB9wO4E8B6ETkCYD+yb4wnRFWfB3DRSb7fh6w+6DjOaUypPA5P5+3wm0QkjuysbhDZqJFTElOsQSCVF187bGJz655nTSQVZQ3E+slZTWmokXWPHuOXV7PH+NEZncRqdiHjJ2j9CE9sYEe30LDJnxe1GhwvH67h9i3awX6D4/zYKjg/YPQEn0/Xpbzc+p2FhowfmcmfZ+sC25ogNh+gjQW2foBWA3zuDk5dufFjf0l2xXI+34q13B+RT7OfYO9qU+Olj/tjqN5cMHO9g6Y/Btay90J8l9FYG/j6ScT46ZnBYrDeapZGlzMap6XvIvY7jLZzf4TWNJMdOMGxxekdL738dzaL3ulNQU1QRKpE5CMi8gUReTOyg9/NAPYAuPFUNNBxnBLlNNEEvw3gBIDfAvgLAP+ArL/6H6jqc8VtmuM4JUsJvR2ebBBco6obACBXY6QNwEpVHS68meM4C57TZBB8WWxT1bSIHJ6LATA0ACzeMvbkHh7i3i07dILsSH0D2dEu1gzt/1DxdlYFrF9Yxe5usvuWGw3N5JOLHjd1i43mFzS2GL/BEeMHaP3IYp0sCkU6uMZGsoFjkSsP8fnHDrGGNryUY19jO4+SPbBhKdnRl0ydZoOtC2xrgth8gDYW2Gq2VgPc/HEONV/zn+8lu2yL0Rx3ccLH2ko+H/tjDZr7K5Dk/o62sOY33MyxyMjw+k2/NZqq0awrD7MmGTMaXt9K9hu0lJnrHz1kNNoqjj3WExy7LiHjV7kxrybK9l8XPHZBTpNB8AIRGf3FCIBYzhZkw4OrJt7UcZyFiuA0eRxWtbEUjuM4U6REBsEi5dV2HMcpDabqLD2npCoUna8e0zkinTYYkzXAY5eYGh3HON9dWbeJPV3Hu6vdwfbg5az59GwwsZr9Jn+hqSkSGjSxwSxhYrjeaIRpExvNh4MK72+4xvipdbAm1L7R5ENcxuczPlaWNbO+laZu8j7WqIZW8/FtrLGtC2xrgth8gDYW2PoBWg1w3x9+ley1te8hu+5Z1jy7zrP5JgvHGsdO8PkPXcbB4baGTaq8nuy21xhNt51VpI5LTKx6O+f/62s2muQxXj+YNH6nQ7bmiNkeTGJ5DdmSybsfgzOcJxXp7bCIXAPgcwCCAL6mqp+cYL0/QjZl36WqurnQPn0m6DhOcZhlP8FcCO4XAfwegHMAvEtEzjnJepXIVsp8airN9EHQcZziMPvO0pcB2KOq+1R1BMA9AK4/yXr/DOBfMJb9viA+CDqOUxRmWGOkTkQ2533y86wtA3Aozx6XlFlELgawQlV/MtV2loQmKCkhHdBqWJUHOD9bIG1qgARNHdpe9ssKJHmHA0t5/ZU/Pk52z1msgdVwejaM1JjY317+L668k49f1mPyIQZ4/fAg25Ee1vxiLd3cABO7W9HITwy1L3J/JepYM7V+cTW7ucaGDLMfW2w/i5yZfQfJ1tddQratC2xrgth8gDYW2PoBWg1wzxu/QfY1/4tzfUTYTQ5Nv+Qv2i9njS1+oJ/sTJQ11pFqtss6h43NmmTwud1sX3UB2SETK167nfurvIuvv40NHllkYqNNbHCmnzXa0BPPcXtqx85fhl5B7PDMNMFOVd04kw1FJADg0wBumc52PhN0HGf2mcmj8OSD5hEAK/Jsm5S5EsB5AJ4QkRYArwbwgIgUHFRLYiboOE7pUYS3w5sAnCkiq5Ed/G4C8MejC1W1B8DL4Vwi8gSA/9vfDjuOMzfM8kwwV5LjdgAPI5ul/l5V3S4iHxeR62bazJKYCYb7FMt+MeYsFzI1LoK9rHlU7jQa1jHW9JLruMBd3T7WtBa1HCIbZ64mc9VPTb66VvZjg4nFTMdNffmnXyDTVJUdX9fhVRvIDO7hsixSxZqTNrHfZMOjfD46wv0T2s79pyvZTzCzjQs1pwJ8fqEmzheoSd7/+P5iDc7WBbZ+ejYfoI0Ftn6AVgN86CffIfttl3PhxNQB7p/GbayRaoLbb6tumLLVCETZE6++jq9fZoA1uabfmP4xdYtlhO/3dLXxszzQQXb4oImF7+b+tmNNaAVf7/z+UDW5NKdBMfwEVfVBAA+a7z42wbpXTmWfJTEIOo5TgpRI2JwPgo7jzD5zmCR1uvgg6DjOrCMYLxvMV0piEBypBfa/M+8L4WZXbuf8fmkjwUV6WFPq3sAaS7ibYz1rdrLdv5Iv5/Ba1tDi21iD619n8hcOs0YTfesVvL86FsHK23j9qhZePrTxLLKr9/Lxyrq4fUeu4/OpfYl1Hls3uOkx9jPsuJ3bW3bCiHbmbl/8BH+x+xa2rzmXNc3t/8SxzLYusK0JYvMB2lhg6wdoNcCf/PZHZJ/xc/YzbKrjHfQ9zLHCKSPiBk1cQtUBU8PlXG7fip4LyT52IauK0U4+X1sDJ2Bc9+Qyzh9Z/6wpchNgzTbUxcsT9RyrHImOtUdafoUZ4zNBx3EWMqdFPkHHcZwZ44Og4zgLGh8EZ4/AsKBi55hOYv3ILIPNrJFFt/BpVuxl29YFrm5hzWTE1GgY7uEN0iZBW1krLxframXab/MPWo3J1sUNpPjuSpqaJJkIt9fW9LAaXlWLqUlhYn+DV7DmWtZj6u42Gb/IxhqyK7ey5vXzAxeTXVfG5xdvNTVUzPnbH5fNB2hjga0foNUA917FscZnfZ1rmoSNI6CtgWL709YIiaxgzS0VN/kZ203+P5M/sozdXMedv83/riZW3uZ3DCfYjzMwYm7g43n9l5q5n2CpUBKDoOM4JcZpVHLTcRxnZvgg6DjOQsZngrNIpgwYaB7TJkJ97Ee27HEWUdIR9rOKt7Pm1XuGiX01NS4ypq5CxOQDDA4Zvz+jEfWeYWqGpEw+uDYjIplYYeuHFu7h8zuxljWmuhdYgwoN8vqd55k6xK3cvsEGk79vCfs9puKmZkjc1MToYk0rEzH5GZebursN3N7EAdakbL7IoXrWWG1dYKsR23yANhbY+gFaDXDXraau8c/+nOxUi6lZc3432S2ra8iu2sft7Wk258MS3bgaIcOLTX93FB5dxmmA/Xw/9J3HfqPhAdb99Kyx2Hp91gjm08EHQcdxFjI+E3QcZ+FSQrHDRc8nKCJBEXlWRH6cs1eLyFMiskdEvi8ikcn24ThOCTL7maWLwqmYCX4A2QSIo8VW/wXAZ1T1HhH5CoBbAXx5oo0BINQPNP56TOcYYYkLZR3s1xc9zjpGvIUd5eriheu0RvdyfrZUeRPZA93cbWFTE6LuWW7fQJPVANm0dYiDCVNjpIPbHz1uaqhE+P+yQCeLTPE23t9IBa+/+AVTQ6OGOzhqNL9kuamje5A1p9BxzpdXv5nbO1LFGmBVC7c3Ey6smdoaKLYusK0JYvMB2lhg6wdoNcB9b7qL7Is+8Vdk96xkjba82/TPcdbchqu5/ysP8fmHe9kOJsZlnOTlw7z/ZBX3R7SFHQ0D9dz/0UMm2Dow1j5JTeKUOwGC0nkcLupMUESWA3gbgK/lbAFwFbJFkQHgbgDvKGYbHMeZI3wmCAD4LIC/Q7YACgAsBtCdS5MNnKRk3ii5Unu3AUCkvPZkqziOM48RLY2pYNFmgiLydgAdqrplJtur6p2qulFVN4ai8VluneM4RaU41eaKQjFngq8BcJ2IXAsgiqwm+DkANSISys0Gbcm8kyIKBJJjPZQ0mpaWGT+3I6yRyMGjZKcuriG7ejdrWBlbtzfI+x9YwhphWY+pI9zNGlkgZUQnI3lZTdJuj6PcnkVbWfNU075AB4uM5R1VZMfaWUOVJLc/s5ULKceaLuXjWckube7ebq65EkxyvsCRKr5+ZW28/uCaGt7earYtXWQPXcYan60LbJN7Wj9M6+dp/QCtBvjsR79E9vqvsZ9h0KTzCw1w+2vbWfOzdYuH6vg/fes3GRq2sct8/WNtps5wJWuWFS/w7wEmlhj591N6Zpog4JogVPUjqrpcVZuRLY33c1V9N4DHAdyQW+1mAPcXqw2O48whJTITnIuSmx8G8EER2YOsRvj1OWiD4zhFRnT6n7nglDhLq+oTAJ7I/b0PwGWn4riO48whJfI4XBIRI5JWhPvHtIkaE4sZOsJ+UIm1y8kuO3sl2RWtrLkNNbJIFD/IGtpwM2taAQ5FHud3Z2OX0xHrN8Y7GKkwCeFqWSMKG82v9yz2cyw7wfuLJnh533Kbb5A1onCv2b6eY0v7lhu/M1NjxNY0QdrEogZsLKvxW2zkusHl+7rJHljL3gH2etj+tRqbrQts8zWO89s0scDWD9BqgDv/D3ZzvWTLjdyedq5xk7H3Q7vR8MJ8P9r8jeEe1vB6zuT2pcuMH+bznWQPnM01R8pbusnGkfa8xswwn2AJpdKai8dhx3GceUNJzAQdxylBSmQm6IOg4zizTimFzZXMIBjI80VLlLPGpeWsgdj8cjbWNNzHmmDQ5N+zmlZwiDWzSD93W8z4ffWtKuwXaDUy6wc2aGqGxHrYj85uP66mRIw1MVuzIpgwscCV5nyihdtfbvzQJGPu9gbW7OKH2XEu1sHnF+4c5O2T3N/xXewXiAy3P1XOGmZZJ7cvEDUamakLbGuC2HyANhbY+gFaDXDLJfeS/aZPcyxy53ncnmCC+7tiH1/vvrWsUZeZWG7rZ2r9EjXGGnXU9I8M8flnmvPqOr/0CvKblEjESMkMgo7jlBY+E3QcZ+FSQvkEfRB0HKcoTFYad75QEoNgMi44eumYNmFrMgBcF3fI1GQ4cRb7odmaD71rrN9WDdnBEV6/62LWlCr3sl9X/2peXrWL2xNMsD1SaY5vNJ/9/519y+ue4/YMLDF+ccdZ40ma/Q8oX3ZbF3mkmTW2ilbW6NJR3n5oCWtclT/jusWHbjiPbOtnGRqapMZIA59v02+5f9pew+dX1snXu75uA9ld57ImaesC2/vD5gO0mpv1A7Qa4M++w/kI33b575PdcTX7tWZC3H5b13mwgfM9RvpMzRU1nm8mH2OwlTXWdAP7YQbb8vwKU+ZiTQefCTqOs5BxTdBxnIWLwt8OO46zsPGZ4CwSHlA0PDPm7Kamrmr0R0+TXfHaC8kO9Rq/KJM/r2a3iaVt5Xx82sc1K+o3c/66QCfXaMjUcexuspY1r1A3a3a6njUgq5md+fkDZB978yqyF21lvzIbDLl4Gx8veoTX7z2XNa3wUT6fkRrWXC1VmzglpC5hTXHZL7n/EyY2uuq5drJ1gP0GxcRia5TtWDv70QWf2012ZoDzRa7ouZDsVJx/BrYusK0JYvMB2lhg6wdoNcCf/PZHZL/lhpvJDh9izS7WzppdJmw05WrTH0e4Jo2Gjaa4gZO5x36zi2wsqhn7W15BZK0Pgo7jLFQ8YsRxnIWNqmuCjuMsbHwmOIuMVAqOXDnW1CBLXIituoJs6ya1eBtrIofeYjSUNtZ0KtrY72+4hpf3rbEtZA3Q+t2FBnj76n18fOsn2LeKT0Ay7Edmz6/HaIqLnmQ/vdQK3r5/ncnPZzSv2FI+n1TU+uGxaDlwHtdcKd/N+esOXc2xsaFB3t9QLW+frODl9sdUeZg7uOMS44d51QVkN/2Gb5hjF3J7ytvZ78/6odq6wDZfoc0HaGOBrR+g1QAf+cHdZK/+r9vIji9hTbOmnO3UvzeQXXaC76+BZXw/V29pI3vkgjPIJg39qGuCjuM4M6JUZoKeVNVxnAWNzwQdx5l9FIBNsTZPKYlBMJACYh1jOpHN51b7Ems+A02siQRGWEOKH2TNqbyTNaHyVtZ4yttN/r4A+4HFO3j74RqeYEf6efm4OrLmfHCcb574XlOXt541vbJeK0KyBhrtMLHEVSa/nKlhkQly+6OmJkqom/34rJ8d+lizivQYDbCRj7fsce4AG4s8WG80S5P/L9bOfp6hAe6/iMlXGO0094fJt2jz84V7WRMcVxfY1ASx+QBtLLD1A7Qa4P533En2uV/guseJo+wXWdXK7Qsav9iKlMkv2M/XJ9THschyoHXMGDGdMx1KYwwsjUHQcZzSo1Q0QR8EHccpDiXiJ+gvRhzHKQrFKL4uIteIyC4R2SMid5xk+QdFZIeIPC8ij4nIqpPtJ5+SmAlmIkD/qjHdK9bGmlfnBaYmhu1MMfn+Vpq6t8ZPzmp+iWqTj4/dvpCo5faMVPP+AynePs6htuNqhIywhIShVfyFzQ/YbWqKxA6wZnXkDWyXdZu6zUaTrNrBsaddl3DNkHA3a2pDDXz8xBvYkXKkio8X7uX+7mvm69O1wfRXq/GjXMnXp6+ZNa/a7cbPcIQ1zaFG4/fIZasxvNjmf+T2WU3X1gW2NUFsPkAbC2z9AK0GuP32L5H9obaLyX7qf1zK7R3m9o7U8vUq32P6o4nvj/Tq9WN///xhzIgiZJYWkSCALwJ4M4DDADaJyAOquiNvtWcBbFTVQRH5SwD/E8A7C+3XZ4KO48w62dhhnfZnEi4DsEdV96nqCIB7AFyfv4KqPq6qo2/CfgfATFnGUxIzQcdxSpCZpdevE5HNefadqjr6unwZgEN5yw4DeFWBfd0K4KeTHdAHQcdxisIUZnYno1NVN77iY4v8CYCNAN4w2bqlMQhmgODg2JN7yJSptXVxbT6+iiPsR9V9FmsktiBMua0jvII1qDBLZsiwJDa+7rHRBENDpiaEcSrNmHyJampEWD84G2ubrDf5ETkdItJlpk6xbb+pQ2tjlQeXs4aktv+T3L7yNt5BIGV+HGZ7qwFONqOIHuP9l3eZmijVrJEFTOy51a5iHYV/vKFho3H28P1i6wLbmiA2H6CNBbZ+gFYD/FTTM2RvrOXJUNUO43fZaGq4rF1B9kiV8SvtGvMNlPQMhb3iVJs7AiC/8ctz3xEi8iYA/wDgDapqr/Y4XBN0HKcI6Fg6rel8CrMJwJkislpEIgBuAvBA/goichGArwK4TlU7ptLS0pgJOo5Tcsy2s7SqpkTkdgAPAwgCuEtVt4vIxwFsVtUHAPz/ACoA/Idkn6AOqup1hfZb1EFQRFoA9AFIA0ip6kYRWQTg+wCaAbQAuFFVT0y0D8dxSpQiOEur6oMAHjTffSzv7zdNd5+nYib4RlXNTzB3B4DHVPWTOWfHOwB8uNAOQsNA7YtjdjrKnVt5kDUgW/MhOMjL67YaEcyIF5EXuKZH2dp1ZA/VF/YzC5v8gUGrIQ0aTS9WWDOs2M41OJIVnH+v6gDHdwYSHEtcbmKbyzuM5rnc+FmaeNGIqbPbvZZvm8pDZv/3PUX24HsvJzt6wsTmDrId6ygcr1rWwRpaMMn5D6Mmtjh0gJ+K5DLW6JQlsXEEh01stvAGPWea2GUTe2zrAtuaIDYfoI0Ftn6AVgPc/E9fJnvdN/+S7OQibn/TE3z+FQdYQwwfPPby34HhGcYOa+kUX58LTfB6AKNZJO8G8I45aIPjOMVm9jXBolDsQVABPCIiW0RkNFVGo6qOprY9CqDxZBuKyG0isllENqeGB062iuM48xmdwWcOKPbj8GtV9YiINAB4VER25i9UVRU5uXyac5C8EwDidStKIxLbcZySo6iDoKoeyf3bISL3IRv20i4iTaraJiJNACZ9jR1IKuJH87QJ6zfXz7rFcIPRuAxVO7rJTtVET75ijsUv8Ey0r5k1oKo97IiXNrG8mQhPuBM13O3xo6xZBhMm/5vJD1ht2m8JdPB7pvIyrjNrNcO6X7WSnVrBdYYrH95BdsX6Zt5fN/ePUdCwaAfntwsk+HxtvsfACXbE7LuINdDoIZPf0GhwI4tYcwsf5P6vf5Y1MBu7bf00k1V8vWJtfD7pMlNX2mioYhwtbV1gWxPE5gO0scDWD9BqgC/dwhrhm97953y8Fs5nqF3mvWQ8//6e+fxjhs7Sp5yiPQ6LSFxEKkf/BvAWANuQ9esZrTRzM4D7i9UGx3HmkBLRBIs5E2wEcF/OVycE4Luq+pCIbAJwr4jcCuAAgBuL2AbHceYCxUxjh085RRsEVXUfgAtO8n0XgKuLdVzHceYewZSywswLSiJiRAOCVPmYLjZcy0/xdZtZY4mUsYYWPswayMAG1pgiPawppjt5/f6rziS719QFTpZzvr/gCF/8YILtqIltDQ6zPbiUNaZIO59PooFjdyPd7BeXGeDgaquZBvpZczr2+qVk1+5ija/77eeSXXnQaHyBwqpKzxo+n+oW3r5/FWus8VZTA8X4/WWqbCyw8cs7wfvPdPdwgwLskGA1wLDpr2gLO4JmKo0m/DzXWdYYt99q2Brm62nrAtuaIDYfoI0Ftn6AVgP82XfuIvutf/hnZAfDPAykXto71lY1gfjTwQdBx3EWND4IOo6zYHFN0HGchY5rgrNIMJFGxd4xXafCaCyZ58kHG9FVnC9No6ypxLcf5QPY/ZWz5lPzNPvRVbRwrKr1exv3GGD+R9Qod3ugjzWsmKn7q/sPkR3NmIzhZv3MWVxbJtBrioiY9oWMZmn91Kp38wk0fP4g2Xs+ew6vf7iN7LqH9oJXYA21ytQsyWx9kezQmmay9QRrfNbLM2Pq6tqfYqjLxMomTCz1efVkB+r5CBUv8P0zcDZrjNFO4+fXyhrz4Ab226zewv1l6wLbmiA2H6CNBbZ+gFYDfPg/v0X2OV/imiYrPmGu10zxQdBxnIXL3Pn9TRcfBB3HmX0UPgg6jrPA8Rcjs0cmFCDfuGgr+wVmXncR2f01HLtb1sl+ZqlKU6d3F+frC1SxZtV3wRJebmpo2PyFg/XsB1bRyppO7AC3f2Qp15RIxXj7zBXspxcaZD82G8ts/RTLjrPmlljE5x/pYT+zoRXcntihXrKf+8F5vP717JdY83NuT2oN+2WqiaWWEf61ZN7INTUCxu9PTCx1YnkN2aEnnmN7BftBJkwNlsCIqSEzwP0RPWT8DI2GWN7Sze0zNVrSDVxnOPabXWSPXHAGt7ePNT5bF9jWBLH5AG0ssPUDtBrgjr/iusZv/cSFmA38xYjjOAsbHwQdx1mwKICMD4KO4yxY/O3w7BLgnHwa4WbbfHAVm7hGiA6zppR43Vlkpw6wH16wijWxyqfZLy7ZzH5h9vjRY2Qi1G389FKsOdnY3pFqk7/uly+RnTm7meyaZ/iA2mZSNK5hv8KyXez3qKY96WO8v9TlnAej6VO/QSEaf8eaasc7j0+wZu74EVPzxWhY6R18/rKRNUlbtzlYy36c9vpGoibf5HHW/PQs9uODjY0OmqIkR1hTzjSzBhls49hiLKohM2T8MuUAX5/06vVk59cFBrgmCACTD5BjgYHxfoBWA3y49bmX/77srbbI9zQokUHQ6w47jrOgKY2ZoOM4pUeJzAR9EHQcZ/YpoRcjoiUwWovIMWSzUNcB6Jxk9bnE2zdz5nPbgIXbvlWqWj/5akx1WaNesfTd0z7YQy2f2aKqG6e94SugJGaCoxdBRDaf6g6aDt6+mTOf2wZ4+2ZECUywgBIZBB3HKTFK6HHYB0HHcYqDzwSLwp1z3YBJ8PbNnPncNsDbN31KZBAsiRcjjuOUFtWRBr2i/p3T3u6h1i/4ixHHcU4DFECmNHJp+SDoOE5xKJGnzJIImxORa0Rkl4jsEZE75kF77hKRDhHZlvfdIhF5VER25/6tLbSPIrdvhYg8LiI7RGS7iHxgPrVRRKIi8rSIbM21759y368Wkady1/n7IhKZbF9FbGNQRJ4VkR/Pt7bl2tMiIi+IyHMisjn33by4vi+jOv3PHDDvB0ERCQL4IoDfA3AOgHeJyDmFtyo63wRwjfnuDgCPqeqZAB7L2XNFCsCHVPUcAK8G8P5cn82XNiYAXKWqFwC4EMA1IvJqAP8C4DOquhbACQC3zlH7AOADAPIrPs2nto3yRlW9ME9Dmy/XF4BmXWSm+5kD5v0gCOAyAHtUdZ+qjgC4B8D1c9kgVX0SgE2Ncj2Au3N/3w3gHaeyTfmoapuqPpP7uw/ZH/MyzJM2apb+nBnOfRTAVQB+kPt+ztonIssBvA3A13K2zJe2TcK8uL4Acpm0MtP+zAWlMAguA5CfC+lw7rv5RqOqjtZOPAqgsdDKpwoRaQZwEYCnMI/amHvcfA5AB4BHAewF0K2qo7UI5vI6fxbA32GsSsZizJ+2jaIAHhGRLSJyW+67eXN9AZTMTNBfjBQBVVURmXNVWEQqAPwQwN+oaq/k1Vee6zaqahrAhSJSA+A+AOsLb3FqEJG3A+hQ1S0icuUcN6cQr1XVIyLSAOBREaHi23N9fXONmNPDT5VSGASPAMivNr089918o11EmlS1TUSakJ3hzBkiEkZ2APyOqv5n7ut51UYAUNVuEXkcwOUAakQklJtxzdV1fg2A60TkWmTrulcB+Nw8advLqOqR3L8dInIfsrLR/Lm+qiXjIlMKj8ObAJyZezsXAXATgAfmuE0n4wEAN+f+vhnA/XPVkJyG9XUAL6rqp/MWzYs2ikh9bgYIEYkBeDOyuuXjAG6Yy/ap6kdUdbmqNiN7r/1cVd89H9o2iojERaRy9G8AbwGwDfPk+pYa834mqKopEbkdwMMAggDuUtXtc9kmEfkegCsB1InIYQD/COCTAO4VkVuRTft149y1EK8B8KcAXsjpbgDw95g/bWwCcHfuzX8AwL2q+mMR2QHgHhH5BIBnkR3I5wsfxvxpWyOA+3LyRgjAd1X1IRHZhPlxfbOUyOOwh805jjPrVAfr9NWxt017u0cGvuVhc47jnA54tTnHcRYynk/QcZwFzxw5P08XHwQdx5l1FICWyEywFFxkHMcpNVSzM8HpfiZhsmQqIlKWS3CxJ5fwonmyffog6JwUEUnnMpRsE5Efjfr1zWA/S0XkB5Ov6ZxuaEan/SnEFJOp3ArgRC7RxWeQTXxREB8EnYkYymUoOQ/ZZBHvn8lOVLVVVW+YfE3ntGP2Z4JTSaaSn0TiBwCulvx40ZPgmqAzFX4L4HwAEJEzkP3fuB7AIIC/UNWdue+/AyCObKTC36hqRe5x5Meqep6IRAF8GcBGZNN9fVBVHxeRWwBcB6AcwBkA7lPVv8v9z//13PqKrKP8Z07VSTszpw8nHv6Z/qBuBptGR/Mj5rhTVUfrp5wsmcqrzPYvr5MLtOhBNgHGhDWZfRB0CpIbiK7GWITEnQDep6q7ReRVAL6EbJqpzwH4nKp+T0TeN8Hu3o9sbP8GEVmPbBaUdbllFyKb7SYBYJeIfB5AA4BludkoZvpI7px6VNXm25y3+OOwMxGxXMjdaEqmR3NZaa4A8B+5ZV9FNgQOyCZA+I/c39+dYJ+vBfDvAKCqO5EN7RodBB9T1R5VHQawA8AqAPsArBGRz4vINQB6Z+/0nBJkKslUXl5HREIAqgF0FdqpD4LORAyp6oXIDkaC7CwugGxevQvzPmfP0vESeX+nAYRU9QSACwA8AeB9yCU5dRYsU0mmkp9E4gZkE2AUfOPig6BTEFUdBPDXAD6ErAa4X0T+G5DNViMiF+RW/R2AP8r9fdMEu/slgHfntl0HYCWAXRMdW0TqAARU9YcAPgrg4ld2Nk4pk0tjNppM5UVkE29sF5GPi8h1udW+DmCxiOwB8EFMocSAa4LOpKjqsyLyPIB3ITuIfVlEPopsWvx7AGwF8DcA/l1E/gHAQwB6TrKrL+W2fQHZFyO3qGqiwMu7ZQC+ISKj/1l/ZJZOySlRVPVBAA+a7z6W9/cwgP82nX16FhlnVhCRcmQfoVVEbgLwLlWd01owjjMVfCbozBaXAPhCzierG8Cfz21zHGdq+EzQcZwFjb8YcRxnQeODoOM4CxofBB3HWdD4IOg4zoLGB0HHcRY0/xv8z5jkXZ5XyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "plt.imshow(fc_emp, vmin=0, vmax=1)\n",
    "plt.title(\"Empirical FC\")\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"Regions\")\n",
    "plt.ylabel(\"Regions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Structural connectivity"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "conn = connectivity.Connectivity.from_file(f\"{datadir}/conn/allen84.zip\")\n",
    "\n",
    "np.fill_diagonal(conn.weights, 0.)\n",
    "conn.weights /= np.max(conn.weights)\n",
    "nreg = conn.weights.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_conn(name):\n",
    "    conn = connectivity.Connectivity.from_file(f\"{datadir}/conn/{name}.zip\")\n",
    "    np.fill_diagonal(conn.weights, 0.)\n",
    "    conn.weights /= np.max(conn.weights)\n",
    "    nreg = conn.weights.shape[0]\n",
    "    return conn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING  File 'average_orientations' not found in ZIP.\n",
      "WARNING  File 'cortical' not found in ZIP.\n",
      "WARNING  File 'hemispheres' not found in ZIP.\n",
      "WARNING  File 'areas' not found in ZIP.\n"
     ]
    }
   ],
   "source": [
    "separator   = ''\n",
    "str_subj_sc = separator.join([subj,\"_det\"])\n",
    "conn        = load_conn(str_subj_sc)\n",
    "np.fill_diagonal(conn.weights, 0.)\n",
    "conn.weights /= np.max(conn.weights)\n",
    "nreg = conn.weights.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stochastic simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us run one simulation of the whole-brain network with noise. We also add the BOLD monitor, so that we have simulated time series that can be compared to the empirical fMRI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate(conn, G, simlength, dt, Io, w, sigma):\n",
    "    ...\n",
    "    return [(t1, tavg), (t2, bold)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU TIME = 4141.17589\n"
     ]
    }
   ],
   "source": [
    "T = 120000\n",
    "t0 = time.time()\n",
    "(t1, tavg), (t2, bold) = simulate(conn, ..., T, ...)\n",
    "CPU_TIME = time.time() - t0\n",
    "print(f\"CPU TIME = {CPU_TIME:.5f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results for all of the regions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the timeseries for one region only:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can calculate the simulated functional connectivity. To avoid the interval affected by the initial conditions, we use only the time series starting at 20 seconds after the simulation start."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The comparison between the simulated and empirical FC is bit difficult: due to the technical constraint of the imaging methods, the structural and functional connectome are defined on different brain regions, which are only partially overlapping. We can plot the simulated and empirical FC on these overlapping regions only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def common_regions(regions1, regions2):\n",
    "    common_regs = [reg for reg in regions1 if reg in regions2]\n",
    "    ncommon = len(common_regs)\n",
    "    inds1 = [list(regions1).index(reg) for reg in common_regs]\n",
    "    inds2 = [list(regions2).index(reg) for reg in common_regs]\n",
    "    \n",
    "    return inds1, inds2\n",
    "\n",
    "inds_sc, inds_fc = common_regions(conn.region_labels, reg_emp)\n",
    "\n",
    "# ncommon = len(inds_sc)\n",
    "# offdiag_mask = ~np.eye(ncommon, dtype=bool)\n",
    "# fc_sim_r  = fc_sim[inds_sc, :][:, inds_sc]\n",
    "# fc_emp_r  = fc_emp[inds_fc, :][:, inds_fc]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let us calculate the predictive power, or the correlation coefficient of the simulated FC with the empirical FC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predictive_power(conn1, conn2, regions1=None, regions2=None):\n",
    "    \"\"\"Calculates the predictive power of the connectome `conn1` to `conn2`. \n",
    "    \n",
    "    If they are defined on different regions, list of region names in both\n",
    "    connectomes (`regions1` and `regions2`) must be provided.\n",
    "    \"\"\"\n",
    "    \n",
    "    if regions1 is None and regions2 is None:\n",
    "        assert conn1.shape == conn2.shape\n",
    "        ncommon = conn1.shape[0]\n",
    "        inds1 = inds2 = np.r_[:ncommon]\n",
    "    else:\n",
    "        inds1, inds2 = common_regions(regions1, regions2)\n",
    "        ncommon = len(inds1)\n",
    "    \n",
    "    conn1 = conn1[inds1, :][:, inds1]\n",
    "    conn2 = conn2[inds2, :][:, inds2]\n",
    "    \n",
    "    offdiag_mask = ~np.eye(ncommon, dtype=bool)\n",
    "    \n",
    "    pp = np.corrcoef(conn1[offdiag_mask],\n",
    "                     conn2[offdiag_mask])[1, 0]\n",
    "    \n",
    "    return pp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding the best fit"
   ]
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}