fc_pgm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functional Connectivity, PGM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from nilearn.connectome import ConnectivityMeasure\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Average functional connectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "connectome_measure = ConnectivityMeasure(kind='correlation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_data():\n",
    "    data = []\n",
    "    for f in os.listdir(\"PCA_results/\"):\n",
    "        if not f.startswith(\".\") and f.startswith(\"sub\"):\n",
    "            data.append(pd.read_csv(os.path.join(\"PCA_results\", f)).values)\n",
    "    return np.stack(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = extract_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove procedual \n",
    "data = np.delete(data, 4, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reorder to [\"Visual\", \"Goal\", \"Retrieval\", \"Imaginal\", \"Manual\"] [3, 1, 4, 20]\n",
    "reorder = [3, 1, 4, 2, 0]\n",
    "data = data[:, :, reorder]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation = connectome_measure.fit_transform(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_correlation = np.mean(correlation, axis=0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heatmap(data, row_labels, col_labels, ax=None,\n",
    "            cbar_kw={}, cbarlabel=\"\", **kwargs):\n",
    "    \"\"\"\n",
    "    Create a heatmap from a numpy array and two lists of labels.\n",
    "\n",
    "    Arguments:\n",
    "        data       : A 2D numpy array of shape (N,M)\n",
    "        row_labels : A list or array of length N with the labels\n",
    "                     for the rows\n",
    "        col_labels : A list or array of length M with the labels\n",
    "                     for the columns\n",
    "    Optional arguments:\n",
    "        ax         : A matplotlib.axes.Axes instance to which the heatmap\n",
    "                     is plotted. If not provided, use current axes or\n",
    "                     create a new one.\n",
    "        cbar_kw    : A dictionary with arguments to\n",
    "                     :meth:`matplotlib.Figure.colorbar`.\n",
    "        cbarlabel  : The label for the colorbar\n",
    "    All other arguments are directly passed on to the imshow call.\n",
    "    \"\"\"\n",
    "    # bigger fig\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(8, 8)\n",
    "\n",
    "    if not ax:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    # Plot the heatmap\n",
    "    im = ax.imshow(data, **kwargs)\n",
    "\n",
    "    # Create colorbar\n",
    "    cbar = ax.figure.colorbar(im, ax=ax, **cbar_kw)\n",
    "    cbar.ax.set_ylabel(cbarlabel, rotation=-90, va=\"bottom\")\n",
    "\n",
    "    # We want to show all ticks...\n",
    "    ax.set_xticks(np.arange(data.shape[1]))\n",
    "    ax.set_yticks(np.arange(data.shape[0]))\n",
    "    # ... and label them with the respective list entries.\n",
    "    ax.set_xticklabels(col_labels)\n",
    "    ax.set_yticklabels(row_labels)\n",
    "\n",
    "    # Let the horizontal axes labeling appear on bottom.\n",
    "    ax.tick_params(top=False, bottom=True,\n",
    "                   labeltop=False, labelbottom=True)\n",
    "\n",
    "    # Rotate the tick labels and set their alignment.\n",
    "    plt.setp(ax.get_xticklabels(), rotation=-30, ha=\"left\",\n",
    "             rotation_mode=\"anchor\")\n",
    "\n",
    "    # Turn spines off and create white grid.\n",
    "    for edge, spine in ax.spines.items():\n",
    "        spine.set_visible(False)\n",
    "\n",
    "    ax.set_xticks(np.arange(data.shape[1]+1)-.5, minor=True)\n",
    "    ax.set_yticks(np.arange(data.shape[0]+1)-.5, minor=True)\n",
    "    ax.grid(which=\"minor\", color=\"w\", linestyle='-', linewidth=3)\n",
    "    ax.tick_params(which=\"minor\", bottom=False, left=False)\n",
    "\n",
    "    return im, cbar\n",
    "\n",
    "\n",
    "def annotate_heatmap(im, data=None, valfmt=\"{x:.2f}\",\n",
    "                     textcolors=[\"white\", \"black\"],\n",
    "                     threshold=None, **textkw):\n",
    "    \"\"\"\n",
    "    A function to annotate a heatmap.\n",
    "\n",
    "    Arguments:\n",
    "        im         : The AxesImage to be labeled.\n",
    "    Optional arguments:\n",
    "        data       : Data used to annotate. If None, the image's data is used.\n",
    "        valfmt     : The format of the annotations inside the heatmap.\n",
    "                     This should either use the string format method, e.g.\n",
    "                     \"$ {x:.2f}\", or be a :class:`matplotlib.ticker.Formatter`.\n",
    "        textcolors : A list or array of two color specifications. The first is\n",
    "                     used for values below a threshold, the second for those\n",
    "                     above.\n",
    "        threshold  : Value in data units according to which the colors from\n",
    "                     textcolors are applied. If None (the default) uses the\n",
    "                     middle of the colormap as separation.\n",
    "\n",
    "    Further arguments are passed on to the created text labels.\n",
    "    \"\"\"\n",
    "\n",
    "    if not isinstance(data, (list, np.ndarray)):\n",
    "        data = im.get_array()\n",
    "\n",
    "    # Normalize the threshold to the images color range.\n",
    "    if threshold is not None:\n",
    "        threshold = im.norm(threshold)\n",
    "    else:\n",
    "        threshold = im.norm(data.max())/2.\n",
    "\n",
    "    # Set default alignment to center, but allow it to be\n",
    "    # overwritten by textkw.\n",
    "    kw = dict(horizontalalignment=\"center\",\n",
    "              verticalalignment=\"center\")\n",
    "    kw.update(textkw)\n",
    "\n",
    "    # Get the formatter in case a string is supplied\n",
    "    if isinstance(valfmt, str):\n",
    "        valfmt = matplotlib.ticker.StrMethodFormatter(valfmt)\n",
    "\n",
    "    # Loop over the data and create a `Text` for each \"pixel\".\n",
    "    # Change the text's color depending on the data.\n",
    "    texts = []\n",
    "    for i in range(data.shape[0]):\n",
    "        for j in range(data.shape[1]):\n",
    "            kw.update(color=textcolors[im.norm(data[i, j]) > 0.01])\n",
    "            text = im.axes.text(j, i, valfmt(data[i, j], None), **kw)\n",
    "            texts.append(text)\n",
    "\n",
    "    return texts    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:109: DeprecationWarning: In future, it will be an error for 'np.bool_' scalars to be interpreted as an index\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcUAAAGKCAYAAAB9z1o9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl8VNX5x/HPMxMSCGFJ2AQCymYFERcWqW2tG6JWsBtWrFpqf+5Iq9WqbbXS1tpaf12s2qotiLWKVvsTVAQUxdq6AKKCLAqCQNiXsG/JzPP7417CZCMBcjOMfN+v131l7jnn3jl3Mskzz7nn3jF3R0RERCCW7g6IiIgcKhQURUREQgqKIiIiIQVFERGRkIKiiIhISEFRREQkpKAon3lm9hczu70enudTMzsr6udJFzPraGZbzSxei7b18pqL1DUFRYlUGCh2hP9M9yztIny+YWb2n9Qyd7/a3X8R1XPWlpn1M7MJZrbRzDaY2TQz+266+1WdikHe3Ze6e567J2raNvU1N7PTzKwoyr6K1BUFRakPg8J/pnuWFenuUH0zs88DrwKvA12BFsA1wLnp7JeIlKegKGlRVfaQmpmY2Z1m9rSZPWZmW8xsjpn1SWnbwcz+ZWZrzWy9md1vZt2BvwCfDzPSjWHbR83slynbXmFmC8NsbXxq5mpmbmZXm9kCMys2swfMzMK6Lmb2avh868zsH2bWvJaH/FtgjLv/xt3XeeBdd7+wDvo1zMz+Y2b3hnWLzezclG2bmdnfzGylmS03s1+mDoGGzzsvfJ3nmtlJZvZ3oCPwfPha/sjMjgr7kWVmF5nZjAq/vxvMbHzqa25mjYGXgHapIwVmtt3MWqRs2zv8XTao5espEgkFRTmUDQbGAs2B8cD9AOE/9BeAJcBRQHtgrLvPA64G3goz0koBy8zOAO4GLgTahvsYW6HZ+UBf4Piw3cA9m4fbtgO6Ax2AO2s6CDPLBT4PPLOPNgfTL4CTgY+AlsA9wN/2BE1gDFBKkKGeCJwN/E/4vEPCY7gMaErwmq9390uBpezN8u+p0JfxwOfMrFtK2cXAE6mN3H0bQTa8osJIwdTwGPa4hOB3WFLdayRSHxQUpT48F55H22hmz+3Hdv9x9wnhOay/EwQDgH4Egelmd9/m7jvd/T/V7qW8bwOj3H2mu+8CbiPILI9KafNrd9/o7kuB14ATANx9obu/7O673H0t8Dvgy7V4znyCv7WVUfQrtMTdHwlfqzEEgbWNmbUhCEo/CF+rNcDvgYvC7f4HuMfdp4fZ60J3X1LTAbn7dmAcMBQgDI7HEATL2hhDEAj3fMgZSvA7FkkrBUWpD1919+bh8tX92G5VyuPtQEMzyyLI0Ja4e+kB9KUdQRYGgLtvBdYTZJvVPW8egJm1NrOx4RDkZuBxgsysJsVAkiBQ1Xm/KtaFAYuw/kigAbByzwcT4CGgddimA/BJLY6hKk8QBkWCLPG5lOeuyTigh5l1BgYAm9x92gH2Q6TOKChKumwDcveshNlCq1puuwzoGAbIimr62pcVBIFiz/M2Jpj0srwWz3t3uP9e7t6UINOxfW9SFqTeAr4RUb/2ZRmwC2iZ8sGkqbsfm1Lfpbqu17DvyUBLMzuBIDg+UU27Svtx953A0wQZ8qUoS5RDhIKipMvHBJnfV8LJFT8Fcmq57TSCochfm1ljM2toZl8I61YDhWaWXc22TwDfNbMTzCwH+BXwjrt/WovnbQJsBTaaWXvg5lr2F+BHwDAzu3nPBBMzO97M9pw3PJh+VcvdVxIEr/81s6ZmFgsnDO0Z9v0rcFM40cXMrKuZ7QnOq4HO+9h3KcF50t8CBcDL1TRdDbQws2YVyh8DhhGcx3z8AA5PpM4pKEpauPsm4FqCf8rLCTLHWl3LFp43G0QwcWRpuN23wupXgTnAKjNbV8W2U4DbgWcJAmsX9p5fq8lI4CRgE/Ai8K9aboe7vwmcES6LzGwD8DAwoQ76VZPLgGxgLsFQ7jOEQ7nu/k/gLoKgvAV4jiDAQZAZ/zQcdr2pmn0/AZwF/LO64Wx3nw88SXDcG/fMqnX3/xIMK8882OAvUldMXzIsIuliZq8CT7j7X9PdFxFQUBSRNDGzvgRDrh3cfUu6+yMCGj4VkTQwszHAKwSXiiggygExs1FmtsbMPqym3szsvvCmGLPM7KQa96lMUUREMpGZnUow+e0xd+9ZRf15wPXAeQQ3uPiju5+8r30qUxQRkYzk7v8GNuyjyQUEAdPd/W2guZnt63phBUUREfnMak9wLe4eRZS/IUYlVV38/FmmsWIRkRRmhrvXeBOKA3HOOebrKl0YVXvvvsscYGdK0cPu/vB+7KKq49pnHDjcgqKIiNSTdetg+jsHvn0si53u3qfmltUqIriV4R6FBHePqtZhGhQj+VCUAfZ+QLrIqv3Chs+8sf7NsseX2+F5u81R3q/s8dNtLk5jT9LnwtV770q39wtFDi+HwUTL8cDw8M5RJxPcY3dfN+Y/XIOiiIhEzoFEdLs3syeB0wjuwVsE/IzgBvi4+18I7hh1HrCQ4Cb6361pnwqKIiISDQePMCi6+9Aa6h24bn/2qaAoIiLRiTAoRkFBUUREohFxphgFXacoIiISUqYoIiLRiHiiTRQUFEVEJBIOeDLdvdg/CooiIhKNDMwUdU5RREQkpExRREQik2mzTxUURUQkGhk4fKqgKCIikVGmKCIiAkGmmGGzTzXRRkREJKRMUUREIqPhUxEREdBEGxERkVS6o42IiAhkZKaoiTYiIiIhZYoiIhIJz8DvU1RQFBGR6OicooiICJCBmaLOKYqIiISUKYqISHQ0fCoiIoKGTw9nl18OrVtDz55V17vDiBHQtSv06gUzZ+6tGzMGunULljFj6qe/UVjJe0xgBC8ynHn8X7XtlvEWTzGEDXxSVraRJbzCj3mJG5jIjSTYXR9djkQRb/EsQ3iGbzCL6n+hnzKF0ZzMOuYB8AkTGcclZcto+rOej+ur23Xq/V0r+MH6Fxix/nme2za3yjZv7VzKjetf5IfrX+S+TW+Wla9LbOOu4te4Yf2L3Lj+RdYkttZXt+vUwIEDmT9/PgsWLOCWW26pVH/VVVcxa9Ys3nvvPd544w26d+8OQFZWFo8++iizZs1i7ty53HrrrfXd9bqVPIglDSLLFM1sKnC3u09KKfsB0Ato6u7frMPnuhPY6u731tU+99ewYTB8OFx2WdX1L70ECxYEyzvvwDXXBD83bICRI2HGDDCD3r1h8GDIz6/X7h+0JAne5W+cxu00ooCXuY129KEZHcq1K2EHC3iJArqV2/Zt7uNkriefo9jFFox4fR9CnQiO5bcM5E/k0prnGUZHvkRzOpdrV8I25vI0rTi2rKwL59CFcwDYwEKmcDMtOLpe+18Xkp5k1JZ3+Un+6bSINeK24sn0yWlPYVazsjYrS7fw3PY5/Dx/AHmxbDYld5bVPbD5bb7WuAe9stuyM1mCmaXjMA5KLBbjgQceYMCAARQVFTF9+nTGjx/PvHnzyto88cQTPPTQQwAMGjSI3/3ud5x77rkMGTKEnJwcevXqRaNGjZg7dy5PPvkkS5YsSdfhHDhliuU8CVxUoewiYHRdBsRDxamnQkFB9fXjxgUB0wz694eNG2HlSpg0CQYMCLbNzw8eT5xYf/2uKxtYSBOOII82xGlAR77AcmZUajebsRzDBcRpUFa2ig9ozpHkcxQAOTQhlqFBcR1zaUIhTWhPnAZ0ZgBL+XeldjN5iOO4lDg5Ve5nMZPpzNlRdzcSC0s30CYrjzbxPLIszik5HZm+q6hcmyk7F3J2o6PJi2UD0CzWEICi0k0kSNIruy0ADWMNyLHMO8vTr18/Fi5cyOLFiykpKWHs2LFccMEF5dps2bKl7HHjxo1xdwDcncaNGxOPx2nUqBG7d+9m8+bN9dr/w1mUQfEZ4HwzywEws6OAdkCRmX0Ylh1rZtPM7H0zm2Vm3czsqD31YZubwkwQM7vCzKab2Qdm9qyZ5UbY/zq1fDl0SEmaCguDsurKM80ONtCIFmXruRSwg/Xl2hSzmB2spx29y5VvYSUAr/NLJvEj5jEu+g5HZDtraEybsvVcWrONteXarOcjtrGaDnyx2v0s5pWMDYobEttpEdv7p9kilktxcke5NitLt7AysZnbi1/mJxsm8/6uFUF5YguNLZt7N73BLRte4vGt75HMtJtnAu3bt2fZsmVl60VFRbRv375Su2uvvZaFCxdyzz33MGLECACeeeYZtm3bxsqVK1m6dCn33nsvxcXF9db3Opdhw6eRBUV3Xw9Mg3A8KMgSnyK4G94eVwN/dPcTgD5A+Y+Tlf3L3fu6+/HAPOB7NfXDzK40sxlmNuPhhx/e38OoM+6Vy8yqL/9s2HsgTpL3eJQTqDy+7CRYx3z6M4Iz+QXLeYfVzK7PjtaZKn6dWIXXYRp/oC/fr3Yfa/mQOA3Jp0sEPYxeVa9BRUmcVaVb+VnzM/l+s1N4aMs0tiV3k/Ak80rWcmneifwqfyCrE1uZunNx5H2ua1UN+XoVf+wPPvggXbt25ZZbbuGnP/0pEGSZiUSCdu3a0alTJ374wx/SqVOnyPsciXD49ECXdIh6XGLPEOq48OflFerfAn5iZoUEAW9BDecPeprZL4HmQB4waV+NAdz9YWBPNKzN32skCgsh5YMjRUXQrl1QPnVq+fLTTqvv3h28RhUyw+1soBF7x5NL2MEmlvEqdwKwk428wW/4EreQSwta0YMcmgLQlpMoZhFtOK5ej6EuNKY121hdtr6dNeTSsmy9hO0U8wkTuRaAHaznFW7iLO6lJcFEi0W8nLFZIkCLeC7rk9vL1tcnt5Mfa1SuTUEsl24NWpBlMVrH82gXb8rKxBYK4rl0ysqnTTwPgL7ZhSwoLT/ikAmKiorokDIEVFhYyIoVK6ptP3bsWP785z8DcPHFFzNx4kRKS0tZu3Yt//3vf+nTpw+LF2fehwMg4y7JiHr26XPAmWZ2EtDI3WemVrr7E8BgYAcwyczOAEor9KthyuNHgeHufhwwskLdIW3wYHjssSAzfPttaNYM2raFgQNh8mQoLg6WyZODskxTQFe2sJKtrCZBCUv5L+3pU1afTWO+xigG8SCDeJAWdONL3EIBXTiC49nIUkrZRZIEa5lLUwrTeDQHriXd2cwytrCCBCUs4mU6cGpZfTZ5XMxkhvAcQ3iOVvQsFxCdJJ8yhc4MSNchHLQuWQWsKt3CmsRWSj3Bm7uW0ien/O+zb0575pQEHx42J3exMrGFNvE8umYVsNV3szmcePNhyWoK403r/RgO1vTp0+nWrRtHHXUUDRo04KKLLmL8+PHl2nTt2rXs8Ve+8hUWLFgAwNKlSznjjDMAyM3NpX///syfP7/+Ol+HXJliee6+NZyFOoogayzHzDoDi9z9vvBxL+ANoLWZtQC2AucDe6aeNAFWmlkD4NvAIXP2bejQIONbty7I/kaOhJKSoO7qq+G882DChOCSjNxcGD06qCsogNtvh759g/U77tj3hJ1DVYw4J/E9XucunCSdOZ1mdGA2YymgC+3pW+222eTxOc7nZW4FjHacWOm8Y6aIkUV/bmIyI3CSdGMQ+XRmJg/Rku50TAmQVVnFezSmNU2ofP4pU8QtxuVN+vCrjVNJunNao850yGrG01tn0blBAX1yCjk+uy2zdq/ixvUvEsP4dt4JNIkFk44uzTuBXxS/igOdGxRwZqPMG0ZOJBIMHz6cSZMmEY/HGTVqFHPnzmXkyJHMmDGD559/nuHDh3PWWWdRUlJCcXEx3/nOdwB44IEHGD16NB9++CFmxujRo5k9OzNPJ2Qiq2qcu06fwOxrwL+A7u4+P5xw84K79zSz24BLgBJgFXCxu28wsxHACGAxQeD71N3vNLNrgB8BS4DZQBN3H7Yfl2SEB/uZOWm3n/b+ri+yZ9LYj/QamzL5+XKblsaepM8o71f2+Ok2F6exJ+lz4eonyh5n4mUfdcHdMTPcPZIX4MRu5q//4cC3b3Y+77p7n5pb1p3I5zq7+/+REoXc/VOgZ/j4buDuKra5D7ivivI/A3+uovzOOuuwiIjUmUybPJx5FwCJiEhmcEAX74uIiGQmZYoiIhIZDZ+KiIhAcPG+gqKIiEgow84pKiiKiEgkPAMzRU20ERERCSlTFBGR6GRYpqigKCIi0cjA4VMFRRERiY4m2oiIiJCRmaIm2oiISEYys3PM7CMzW2hmt1ZRf6SZTTGzWWY2Nfzu3n1SUBQRkegkD2LZBzOLAw8A5wI9gKFm1qNCs3uBx9y9F/BzqvgCiooUFEVEJBrh8OmBLjXoByx090XuvhsYC1xQoU0PYEr4+LUq6itRUBQRkegcXKbY0sxmpCxXpuy5PbAsZb0oLEv1AfCN8PHXgCbhF9hXSxNtRETkULVuH18yXNUXI3uF9ZuA+81sGPBvgi+tL93XEyooiohIZCKcfVoEdEhZLwRWlHtu9xXA1wHMLA/4hrtv2tdONXwqIiLRcCKbaANMB7qZWSczywYuAsanNjCzlma2J87dBoyqaacKiiIiEgkHPHHgyz737V4KDAcmAfOAp919jpn93MwGh81OAz4ys4+BNsBdNfVZw6ciIhKNiC/ed/cJwIQKZXekPH4GeGZ/9qlMUUREJKRMUUREopNht3lTUBQRkWhk4L1PFRRFRCQ6CooiIiKBTMsUNdFGREQkdJhmihXvBHT4GevfTHcXDgmjvF+6u5B2F65+It1dSDv3w/d/wvW2Orqd77l4P4McpkFRRETqQ6Z93lBQFBGR6ChTPPRdZPt1g4PPjPJDplXdYP5wsfej6yU2OY39SJ/H/eyyxz/JyrD/WnXkrtK9UypGxNaksSfpc1+ydbq7cMg5LIOiiIhEz3WdooiIyF4KiiIiIqDZpyIiIqkyLVPUxfsiIiIhZYoiIhKdDMsUFRRFRCQarov3RURE9lKmKCIiEtBEGxERkQylTFFERCKhO9qIiIik0kQbERGRQKZlijqnKCIiElKmKCIi0dC9T0VERPbKtOFTBUUREYmOJtqIiIgQ3OYtwzJFTbQREREJKVMUEZHI6IbgIiIie2TY8KmCooiIREK3eRMREUmVYcOnmmgjIiISUqYoIiKRybThU2WKdWQl7zGBEbzIcObxf9W2W8ZbPMUQNvBJWdlGlvAKP+YlbmAiN5Jgd310uc5dfjm0bg09e1Zd7w4jRkDXrtCrF8ycubduzBjo1i1Yxoypn/5GZQXTGc/3GMcw5vBUte2W8gb/YCDr+bisrJhFTOIHvMAVvMBVGfte+MQn8pdEd/6cOJo3k7+pVD8r+Si/T7Thr4mT+GviJN5P/hWATb6EUYm+/DVxEg8njmNm8i/13fU6s8Rf5e9+Co/5yczw+yrVz/OxPOI9eNLP4Ek/gzn+eErdUzzm/XnM+zPPq38PHfI8PK94gEs6pDVTNLM2wO+B/kAxsBu4x92rjypV7+co4AV3r+bfcbSSJHiXv3Eat9OIAl7mNtrRh2Z0KNeuhB0s4CUK6FZu27e5j5O5nnyOYhdbMOL1fQh1YtgwGD4cLrus6vqXXoIFC4LlnXfgmmuCnxs2wMiRMGMGmEHv3jB4MOTn12v360SSBNN5gDO4m1xaMpHrKaQ/zTiyXLsStjOf52jBMeW2fZN7OIWbyacLu9icke+FpCeYlLyeobFJNKWQ0cmT6eaDaGU9yrXrYRcyMPancmV5tOWy2H/Ishx2+1YeSfaimw+mibWrz0M4aElPMJVb+SpPk0c7nmIgnX0gBfa5cu26cQGn2d3lynZ6MdO4l28xGTCeYgCdfCANrXk9HkEdUqZYO2ZmwHPAv929s7v3Bi4CCtPVpwO1gYU04QjyaEOcBnTkCyxnRqV2sxnLMVxAnAZlZav4gOYcST5HAZBDE2IZ+I8Q4NRToaCg+vpx44KAaQb9+8PGjbByJUyaBAMGBNvm5wePJ06sv37XpfV8RBPa0YS2xGnAkZzGMt6q1O4DxtCDIcTJLitbybs0pxP5dAEgh6YZ+V5YwTTy6UK+dSZu2fSwb7HAx9dq27hlk2U5AJSyC8+0/6ih1cykOZ1oZkcRt2yO5qssonZv6qW8Rge+TEPLp6E1pwNfZimvRtthKZPO4dMzgN3uXjY+4u5L3P1PZtbQzEab2Wwze8/MTocgIzSzN8xsZrickrbep9jBBhrRomw9lwJ2sL5cm2IWs4P1tKN3ufItrATgdX7JJH7EPMZF3+E0Wb4cOqQkz4WFQVl15ZloB+vJpVXZei4t2cG6cm02sJDtrKWQ/uXKt1AEGK/yYyZwHXN4uj66XOe2sJymtvcX2oT2bKHyL3S+/4tHEifwbGIIm31ZWflmX8YjiRO4P3kk/e1HGZclAmxjFXns7Xce7djKqkrtPuEFnvDTmODfY4sHr9FWVtGkFttmikwbPk1nUDwWmFlN3XUA7n4cMBQYY2YNgTXAAHc/CfgWUHmgvgIzu9LMZpjZjIcffrhuel4rVvbISfIej3IClccVnQTrmE9/RnAmv2A577Ca2fXYz/pT1ZvcrPryTORVzj8v/154l4c4iSsrtUqSYC0fcgq3cDb/SxFvsor3IuxtVPb9GgB0tUFcF1vEFfH36WRn8nzyu2V1Ta0DV8Tf55rYx8z2x9jqqyPub92r6n1Q8S19FGczjBlcbFPpwKm8wvVlW1eWsX8QwfDpgS5pcMhMtDGzB8zsAzObDnwR+DuAu88HlgBHAw2AR8xsNvBPoEd1+9vD3R929z7u3ufKKyv/I6oLjSpkhtvZQCP2jiOWsINNLONV7uR5rmU9C3iD37CBT8ilBa3oQQ5NySKHtpxEMYsi6We6FRbCsr0JAUVF0K5d9eWZKJeWbGdt2fp21pUbRQjeC5/yCj/iOS5jHfN4nZ+xno/JpRVt6EVDmpFFQ9rRlw0sTMdhHJQmFJbL/LawvFzmA5BrLcqGSU+wK1jFu5X3Y+1oZT1YxhvRdjgCebRlKyvK1reygsYcUa5NIysgHr4Gx3IJa5hVtu2WCtvm0aYeeh0NZYq1Nwc4ac+Ku18HnAm0ovqPRTcAq4HjgT6QckImjQroyhZWspXVJChhKf+lPX3K6rNpzNcYxSAeZBAP0oJufIlbKKALR3A8G1lKKbvCTGEuTTPvtGqtDB4Mjz0WvNnffhuaNYO2bWHgQJg8GYqLg2Xy5KAsE7Xgc2xhOVtZRYISljC13DBpNo35Jv/kqzzGV3mMlnTny4ykBUfTlt4Us5hSdpIkwRpm0YyOaTyaA9OOvhSzkI2+mITvZq4/RTcbVK7NVl9Z9ngB42lBdwA2exElvgOAHV5Mkb9JC8pPTskEbTiRjSxiky8h4bv5mOfoRPk39baUDHgxk8gPJ+B15HSWMZWdvpGdvpFlTKUjp9dr/+uKE1yScaBLOqRz9umrwK/M7Bp3/3NYlhv+/DfwbeBVMzsa6Ah8BDQDitw9aWbfgUNjFkKMOCfxPV7nLpwknTmdZnRgNmMpoAvt6Vvtttnk8TnO52VuBYx2nFjpvGOmGDoUpk6FdeuC7G/kSCgpCequvhrOOw8mTAguycjNhdGjg7qCArj9dugbvkx33LHvCTuHshhx+nAdr/JjnCRdOJvmHMUHjKEFR1PI56vdNocmdOfrTOR6gvdCP9pzcv11vo7ELIuzY/cxNnkuSRIcb9+llR3L68mf0dZ6c7QNZrr/iQXJ54mRRUMKOD82CoD1zOOV5M0YhuOcbDfS2o5L8xHtv5hl8WW/m/FcRJIEPRhKCzuGt/03tOZ4Ots5fMAjLPbJGHEa0pyzwrNBDS2fvn4jT4dBtC8/pKFl4FTsDGWexluYm1lbgksyTgbWAtuAvwDjwp+9gVLgRnd/zcy6Ac8C24HXgOvdPW8/LslwgIvsmbo/mAww1r+Zspah5yjqxN73/CU2OY39SJ/H/eyyxz/JyswZngfrrtK9A2UjYmvS2JP0uS/ZmhGxNdyXbB3JP4RjC8yfPvPAt+/5DO+6e5+aW9adtF6n6O4rCS7DqMqwKtovAHqlFN0Wln8KpOUaRRERqUbENwQ3s3OAPxKMGv7V3X9dob4jMAZoHra51d0n7Gufh8xEGxER+Qzyg1j2wcziwAPAuQSTLoeaWcXJlz8Fnnb3EwkSsAdr6q6CooiIZKJ+wEJ3X+Tuu4GxwAUV2jjQNHzcDFKm9VZDNwQXEZHIHOS0lZZmlnp7sIfdfc8F5+2BlIu5KIJKM9PuBCab2fVAY+Csmp5QQVFERKJx8OcU1+1jok1Vk4MqhuChwKPu/r9m9nng72bW0736XikoiohIdKK7wKEIyn3rQiGVh0e/B5wD4O5vhXdGa0lwd7Qq6ZyiiIhEJsI72kwHuplZJzPLJphIU/HO80sJbgqDmXUHGkLKLaeqoKAoIiIZx91LgeHAJGAewSzTOWb2czMbHDb7IXCFmX0APAkM8xouztfwqYiIRGPPDcGj2n1wzeGECmV3pDyeC3xhf/apoCgiIpFw0ndj7wOloCgiIpFRUBQREdkjw26tq4k2IiIiIWWKIiISjTR+WfCBUlAUEZHIKCiKiIjskWFBUecURUREQsoURUQkMlF+yXAUFBRFRCQatfiy4EONgqKIiERCd7QRERFJkWlBURNtREREQsoURUQkOhmWKSooiohIZDJt+FRBUUREouG6JENERGSvDMsUNdFGREQkdFhmimP9m+nuwiEgwz6+ReRxPzvdXUi7u0r12fi+ZOt0d+EzK9POKeqvQUTkMLbE/xHtE/hBLGlwWGaKIiISPdf3KWaGy21auruQFqO8X9njS2xyGnuSXuWHTC1t/Uivvf+p7s/9JI39SJ/h27uUPb7Afp/GnqTPOL8h3V045ByWQVFEROqHMkUREZE9FBRFREQCyhRFRERCmRYUdUmGiIhISJmiiIhEI43XGx4oBUUREYlMpg2fKiiKiEgkHAVFERGRvTIsKGqijYiISEiZooiIREbDpyIiIgC6IbiIiEiKDAuKOqcoIiISUqYoIiKR0fCpiIhISEFRRERkDwVFERGRjLz1qSbaiIiI7KFMUUREoqHrFEVERFIoKIqIiAQyLVPUOUUREYmM+4EvNTGzc8zsIzNbaGa3VlH/ezN9lt1TAAAgAElEQVR7P1w+NrONNe1TmaKIiGQcM4sDDwADgCJgupmNd/e5e9q4+w0p7a8HTqxpv8oURUQkMhFmiv2Ahe6+yN13A2OBC/bRfijwZE07VVCsI0W8xbMM4Rm+wSzGVNvuU6YwmpNZxzwAPmEi47ikbBlNf9bzcX11u06tYDrj+R7jGMYcnqq23VLe4B8MLHecxSxiEj/gBa7gBa4iwe766HKdu/xyaN0aevasut4dRoyArl2hVy+YOXNv3Zgx0K1bsIyp/i2UEeYmXucXOwcwcucZTC75S6X6t0uf5bYdffn1zkH8eucg3izd+34ZV3IPv9p5Lr/aeS7vlr5Yn92uU2uYxxTu5hXuYgFTqm23gg8Yz41sZFm58u0U8yK3spDXou7qoaylmc1IWa5MqWsP5V60orCsEjM7EugEvFrTE9Y4fGpmCWB22HYxcKm7Vzsua2bNgYvd/cF9tHnT3U+p6bn3h5ltdfe8utxnbSVJ8Da/ZSB/IpfWPM8wOvIlmtO5XLsStjGXp2nFsWVlXTiHLpwDwAYWMoWbacHR9dr/upAkwXQe4AzuJpeWTOR6CulPM44s166E7cznOVpwTLlt3+QeTuFm8unCLjZjxOv7EOrEsGEwfDhcdlnV9S+9BAsWBMs778A11wQ/N2yAkSNhxgwwg969YfBgyM+v1+7XiaQn+GfJnVyXPYbmdgS/3fV1joufSdtYt3LtTox/hQuz7yxX9mHiNZYl53BLzvOUsps/7rqYHvFTaWRN6vEIDp6TZBb/4vNcTSOa8W9+zxEcSxOOKNeulJ0s4g3y6VhpH3N4jtZ0r68uR6K25wb3YZ2796mmzqp6ymraXgQ84+6Jmp6wNpniDnc/wd17AhuA62po3xy4tqqKcAyYug6I6baOuTShkCa0J04DOjOApfy7UruZPMRxXEqcnCr3s5jJdObsqLsbifV8RBPa0YS2xGnAkZzGMt6q1O4DxtCDIcTJLitbybs0pxP5dAEgh6bEMjQonnoqFBRUXz9uXBAwzaB/f9i4EVauhEmTYMCAYNv8/ODxxIn11++6tCT5AS3tSFrGOpJl2fSOf4XZiVdqte2q5EK6xvoRtyxyLJf2sWOYl6j8t3SoK2YpjWlJY1oQI4v2nMgqPqzUbj4v0ZXTidGgXPlKZpNLC5rQpr66nImKgA4p64XAimraXkQthk5h/4dP3yIlPTWzm81supnNMrORYfGvgS7hbJ/fmtlpZvaamT1BkHFiZlv3tQ8z+42ZXZvS5k4z+6GZ5ZnZFDObaWazzWxf48f1ZjtraJzy5s2lNdtYW67Nej5iG6vpwBer3c9iXsnYoLiD9eTSqmw9l5bsYF25NhtYyHbWUkj/cuVbKAKMV/kxE7iOOTxdH11Oi+XLoUPKn3FhYVBWXXkm2shq8q1t2XpzO4KNvrpSuw8Sk7h751f4267rKE4G/8vax45hbvJ1dvsOtvoGFiTfodhX1lvf68pONtGI5mXrDWnODjaVa7OJInawkSNSRo4AStnFQl7lcwysl75GLcJzitOBbmbWycyyCQLf+IqNzOxzQD5U8Sm9CrWefRpmeWcCfwvXzwa6EZzsNGC8mZ0K3Ar0dPcTwnanhW16uvviCvusbh9jgT8Ae4ZgLwTOAXYCX3P3zWbWEng7nG2U1ithqnpyS8nsnSTT+ANf5PZq97GWD4nTsCxbyjRezauwtz7JuzzE5/lhpVZJEqzlQ87hT2SRwxRupQXdOKLmiWIZp6p3qln15Zmp8sFYhZGu4+Jn0Dt+Pg0sh/+UPsHfS37EiJzH6R7/EkuTs/ndrgvJswI6xU4kbpk4SX7fr4GT5EPGcSJDK7X7iEl05stkVTOilHEi+u/s7qVmNhyYBMSBUe4+x8x+Dsxw9z0BcigwtrZxojbvtkZm9j5wFPAu8HJYfna4vBeu5xEEuKVV7GNaxYC4r324+9/MrLWZtQNaAcXuvtTMGgC/CgNnkiBrbQOsqq7z4YnZKwEeeughrrzyyuqaHrDGtGYbez8Jb2cNubQsWy9hO8V8wsRwVHkH63mFmziLe2kZnjNYxMsZmyVCkBluT8mOt7OORrQoWy9hB5v4lFf4EQA72MDr/IwvM5JcWtGGXjSkGQDt6MsGFn4mg2JhISxLmRpQVATt2gXlU6eWLz/ttPruXd1ozhHlsruNvopm1rpcm8a292TpKfFvMa7knrL1gQ2uZWCD4G/l0d030MqOirS/UQgyw71TL3aykYY0LVsvZRdbWMV/eQCAXWzhHf7GyXyPYpawgg+Yy/OUsAPDiJNFJ75U78dRF/ygToXs+xSgu08AJlQou6PC+p3784y1CYo73P0EM2sGvEBwTvE+gjTgbnd/KLWxmR1VxT62VbPvKvcRegb4JnAEQeYI8G2CINnb3UvM7FOg4b467+4PAw/vWd1X2wPVku5sZhlbWEEurVjEy3yZX5TVZ5PHxUwuW3+Ja+jLiLKA6CT5lCmcR1UvQ2ZowefYwnK2sopGtGAJU/kCe6+lzaYx3+SfZesvczMncQUtOJo82jKXf1LKTmI0YA2zOIavp+MwIjd4MNx/P1x0UTDBplkzaNsWBg6EH/8YiouDdpMnw913p7evB6pjrBdrfQnrkstobm14N/Eiw7J/V67NJl9TFihnJ6dwhAUjJElPsIPNNLZ8lifnsyI5n2Ma/Lbej+FgNacD21jLNtbTiGYs5z1O4tKy+gY04pyU/xH/5QGOZTDN6cAXub6sfD4TySInYwMiGO7ZNTer1o4660lt1Xpcwt03mdkIYJyZ/ZkgZf2Fmf3D3beaWXugBNgC1HaqWJX7cPc1BIHwEaAl8OWwfTNgTRgQT4cKUxvTJEYW/bmJyYzASdKNQeTTmZk8REu605FT97n9Kt6jMa1pUvVs4owQI04fruNVfoyTpAtn05yj+IAxtOBoCvl8tdvm0ITufJ2JXA8Y7ehHe06uv87XoaFDg4xv3bog+xs5EkpKgrqrr4bzzoMJE4JLMnJzYfTooK6gAG6/Hfr2DdbvuGPfE3YOZXHLYkiDn/Hg7u/iJOgfH0Lb2NG8WPIHOsZ6clz8LF4vHcPsxBRiZNHYmvHt7CBTTFDKH3ZdBEBD8rgs+38zcvg0Rpzj+Dpv8zBOko70oylHMJ+XaE4HjqCaa3Yk7aymYdaKlzqY2fPA0+7+dzP7PvA/YdVW4BJ3/yScVNMLeAl4EbjJ3c+vap/V7SOsm00wJff0cL0l8DzQAHgf+AJwrrt/WstLMhzgcptWQ7PPplHer+zxJTZ5Hy0/2x731GHqjD1xd5D2/t3fn/tJGvuRPsO37z1/f4H9Po09SZ9xfgMX2O8Z5zdE8ofQLSfLf3/EgV9OM2jpxnf3cUlGJGr8CFYx0Lj7oJTHfwT+WMU2F1comlrdPqvbR1h3XIX1dVB1ypGuaxRFRKQafrDDp/Uv88YlREQkIziGo6AoIiLCwU+0qX+696mIiEhImaKIiEQkBp5ZNyFQUBQRkYjonKKIiEhI5xRFREQyljJFERGJSOZligqKIiISkRjonKKIiAi47mgjIiKyR+YFRU20ERERCSlTFBGRiOg6RRERkVAMMmz4VEFRREQiknmZos4pioiIhJQpiohIRDJv9qmCooiIRMIVFEVERPYwoEG6O7FfFBRFRCQaGXhHG020ERERCSlTFBGRiBiu4VMRERGAWMYNnyooiohIRDTRRkREJKSJNiIiIhnrsMwUR3m/dHch7R73s9PdhUOEp7sDaTd8e5d0dyHtxvkN6e7CZ5Jroo2IiGSSFnwhwr0b7gqKIiIi6JKMDPF0m4vT3YW0uHD1E2WPf5KVTGNP0uuu0r2n0u/P/SSNPUmf8kOmlrZ+pNfeofPLbVoa+5E+OpVU2WEZFEVEpD4YaPhUREQENHwqIiIScjfcMyvMZFZvRUQkg2RepqiL90VERELKFEVEJCKZlykqKIqISEQMdE5RREQk4BkWZjKrtyIikkEy7zZvmmgjIiISUlAUEZFIBN+SkXXAS03M7Bwz+8jMFprZrdW0udDM5prZHDN7oqo2qTR8KiIiEbHIzimaWRx4ABgAFAHTzWy8u89NadMNuA34grsXm1nrmvaroCgiIhEx3ONR7bwfsNDdFwGY2VjgAmBuSpsrgAfcvRjA3dfUtFMNn4qIyKGqpZnNSFmuTKlrDyxLWS8Ky1IdDRxtZv81s7fN7JyanlCZooiIROigwsw6d+9TTV1V33nmFdazgG7AaUAh8IaZ9XT3jdU9oYKiiIhEJNIbghcBHVLWC4EVVbR5291LgMVm9hFBkJxe3U41fCoiIhExnPgBLzWYDnQzs05mlg1cBIyv0OY54HQAM2tJMJy6aF87VaYoIiKR8Ahnn7p7qZkNByYBcWCUu88xs58DM9x9fFh3tpnNBRLAze6+fl/7VVAUEZGM5O4TgAkVyu5IeezAjeFSKwqKIiISDSfKSzIioaAoIiIRsdqcGzykKCiKiEhkFBRFRESA4PsUMyso6pIMERGRkDJFERGJjGdY7qWgKCIikXBNtDl8vb9rBY9unUkS54yGXfhq4x6V2ry1cyn/3DYbA47MymdEs1MAWJfYxkObp7EuuR0Dbm3+ZVrH8+r3AOrAJz6Rl5M34CQ43r7HKbFbytXPSj7KFL+FJuE9e/vYtZwQ+x82+RKeTX6TJAmSlNDHruOk2NXpOIQ6MTfxOs+W/JIkCT4fv5CzG5Q/lrdLn2Vcya9pZkcAcGrWJZyS9S0AxpXcw5zEawAMzBpO76yv1G/n68Dll8MLL0Dr1vDhh5Xr3eH734cJEyA3Fx59FE46KagbMwZ++cvg8U9/Ct/5Tr11u84V8Rbv8DucJEczmF5UfTCfMoXX+DGDeJSWdOcTJvIhj5fVb2Ahg3mMFhxdX12vU5/ZTNHMtrp7pP+pzexNdz/lALd9FHjB3Z+p217VLOlJRm15l5/kn06LWCNuK55Mn5z2FGY1K2uzsnQLz22fw8/zB5AXy2ZTcmdZ3QOb3+ZrjXvQK7stO5MlmFV1n9tDW9ITTEpez9DYJJpSyOjkyXTzQbSy8h8OetiFDIz9qVxZHm25LPYfsiyH3b6VR5K96OaDaWLt6vMQ6kTSE/yz5E6uyx5DczuC3+76OsfFz6RtrFu5difGv8KF2XeWK/sw8RrLknO4Jed5StnNH3ddTI/4qTSyJvV4BAdv2DAYPhwuu6zq+pdeggULguWdd+Caa4KfGzbAyJEwYwaYQe/eMHgw5OfXa/frRJIEb/NbBvIncmnN8wyjI1+iOZ3LtSthG3N5mlYcW1bWhXPoQvBlDhtYyBRuztiAmIkOqRB+oAEx3RaWbqBNVh5t4nlkWZxTcjoyfVdRuTZTdi7k7EZHkxfLBqBZrCEARaWbSJCkV3ZbABrGGpBjmZfAr2Aa+XQh3zoTt2x62LdY4BVvQ1i1uGWTZTkAlLILJxllVyO1JPkBLe1IWsY6kmXZ9I5/hdmJV2q17arkQrrG+hG3LHIsl/axY5iX+HfEPa57p54KBQXV148bFwRMM+jfHzZuhJUrYdIkGDAg2DY/P3g8cWL99bsurWMuTSikCe2J04DODGAplX+XM3mI47iUODlV7mcxk+nM2VF3N0KGe+yAl3TY72c1s9PM7HUze9rMPjazX5vZt81smpnNNrMuYbtBZvaOmb1nZq+YWZuwvJWZvWxmM83sITNbEt6oFTPbmvIcU83sGTObb2b/sDB9MrM7zGy6mX1oZg/vKU+nDYnttIjllq23iOVSnNxRrs3K0i2sTGzm9uKX+cmGyby/K7iZ+8rEFhpbNvdueoNbNrzE41vfI+mZFxS2sJymtveG9U1ozxaWV2o33//FI4kTeDYxhM2+96vQNvsyHkmcwP3JI+lvP8rILBFgI6vJt7Zl683tCDb66krtPkhM4u6dX+Fvu66jOBm8F9rHjmFu8nV2+w62+gYWJN+h2FfWW9/ry/Ll0CHluw0KC4Oy6soz0XbW0Jg2Zeu5tGYba8u1Wc9HbGM1HfhitftZzCuZHxSjuyF4JA40FB8PfB84DrgUONrd+wF/Ba4P2/wH6O/uJwJjgR+F5T8DXnX3k4D/AzpW8xwnAj8AegCdgS+E5fe7e1937wk0As7fV0fN7Mo9X1D58MMP7/+R1kLFL/CqShJnVelWftb8TL7f7BQe2jKNbcndJDzJvJK1XJp3Ir/KH8jqxFam7lwcST+jVdWrUP7zSlcbxHWxRVwRf59OdibPJ79bVtfUOnBF/H2uiX3MbH+MrVUEksxQ+XWwCq/DcfEzuLPhVG5r+CKfi3+Bv5cEfxrd41/i2Nhp/G7XhTy6+wY6xU4knoGjBjXxKt4qZtWXZ6Kq/xospT7JNP5AX75f7T7W8iFxGpJPlwh6WH+CyTYHtqTDgQbF6e6+0t13AZ8Ak8Py2cBR4eNCYJKZzQZuhrJB8y8SBEncfSJQXM1zTHP3IndPAu+n7Pf0MAOdDZyRst8qufvD7t7H3ftceeWV+2p6wFrEc1mf3F62vj65nfxYo3JtCmK59MlpT5bFaB3Po128KSsTWyiI59IpK5828TziFqNvdiGLS6t7SQ5dTSgsl/ltYTlNKJ/t5VqLsmHSE+wKVvFu5f1YO1pZD5bxRrQdjkhzjiiX3W30VTSz1uXaNLZ8GoSvwynxb7EsuXc2ysAG13Jrw+cZnjMGx2llR9VHt+tVYSEsS/m+9KIiaNeu+vJM1JjWbGPvB7vtrCGXlmXrJWynmE+YyLX8k6+ylg95hZtYx7yyNot4OcOzxODDwWd++DS0K+VxMmU9yd7JO38iyOqOA64CGobltQ3/qc+RALLMrCHwIPDNcL+PpOw3bbpkFbCqdAtrElsp9QRv7lpKn5zCcm365rRnTknwR7I5uYuViS20iefRNauArb6bzeHEmw9LVlMYb1rvx3Cw2tGXYhay0ReT8N3M9afoZoPKtdmaEiwWMJ4WdAdgsxdR4sFw8w4vpsjfpAWfq7/O16GOsV6s9SWsSy6j1HfzbuJFjoufWa7NJl9T9nh2cgpHBGccSHqCbR58IFqenM+K5HyOiVU/tJapBg+Gxx4LMsO334ZmzaBtWxg4ECZPhuLiYJk8OSjLRC3pzmaWsYUVJChhES/TgVPL6rPJ42ImM4TnGMJztKInZ3EvLcO/CSfJp0yhMwPSdQiHrSjHZppB2Uml1LnI/wEuBH5jZmcD+zO3bE8AXGdmecA3gXqfbVpR3GJc3qQPv9o4laQ7pzXqTIesZjy9dRadGxTQJ6eQ47PbMmv3Km5c/yIxjG/nnUCTWJAtXJp3Ar8ofhUHOjco4MxGmTdcErMszo7dx9jkuSRJcLx9l1Z2LK8nf0Zb683RNpjp/icWJJ8nRhYNKeD82CgA1jOPV5I3YxiOc7LdSGs7Ls1HdGDilsWQBj/jwd3fxUnQPz6EtrGjebHkD3SM9eS4+Fm8XjqG2YkpxMiisTXj29n3AJCglD/sugiAhuRxWfb/ZuTw6dChMHUqrFsXZH8jR0JJSVB39dVw3nnB5RhduwaXZIweHdQVFMDtt0PfvsH6HXfse8LOoSxGFv25icmMwEnSjUHk05mZPERLutMxJUBWZRXv0ZjWZZcvZS7LuEsyzKsayK+qYXhJhpmdBtzk7ueH5VPD9RmpdWZ2AfB7gsD4NtDX3U8zs9bAkwTB8HXgW0And9+1j+e4n+BLIx81s18SfMPyp8AyYIm731nLSzIc4Ok2F9fy5flsuXD1E2WPf5KVeZN56spdpXv/SO/P/SSNPUmf4dtTP3hl6Im7g7b3f9/lNi2N/UifUd6Py20ao7xfJG+Cwlhvv77Bmwe8/a27G77r7n3qsEs1qvXH0D3XKLr7VGBqSvlpKY/L6tx9HDCuil1tAgaG35r8eeD08Nzkvp5jeMrjnwI/raJ/w2p7LCIiUj8yLVNMx9hMR+BpM4sBu4Er0tAHERGpF5k1ElHvQdHdFxBcbiEiInJIybyz+CIikiEMd2WKIiIiwXWKGj4VEREJZFpQzKxpQSIiIhFSpigiIpGp3ZXwhw4FRRERiYajiTYiIiJ7KFMUERFhz+zTzKKJNiIiIiFliiIiEplMyxQVFEVEJDIKiiIiIqHafj3hoUJBUUREIpNZIVETbURERMooUxQRkUgEl2RkVq6ooCgiIpFRUBQREQGCkJhZQVHnFEVERELKFEVEJDLJDMsUFRRFRCQSmmgjIiKSQhfvi4iIAEGemEx3J/aLJtqIiIiEDstM8cLVT6S7C2l3V6k+DwEM394l3V04BGTW8FYURnm/dHfhMyvTMsXDMiiKiEhglJ9MlB+MMm2ijdIFERGJhIfnFA90qYmZnWNmH5nZQjO7tYr6YWa21szeD5f/qWmfh2WmaGbp7kJapM4CGxFbk8aepNd9ydZljy+w36exJ+kzzm8oe3y5TUtjT9Kn/JDp4fk/oT6Gzp1EJPs1szjwADAAKAKmm9l4d59boelT7j68tvtVpigiIpmoH7DQ3Re5+25gLHDBwe5UQVFERCJy0MOnLc1sRspyZcrO2wPLUtaLwrKKvmFms8zsGTPrUFOPD8vhUxERqR/Jg5t9us7d+1RTV9WYd8Xx4OeBJ919l5ldDYwBztjXEyooiohIJBzHPZpzigSZYWrmVwisKPf87utTVh8BflPTTjV8KiIimWg60M3MOplZNnARMD61gZm1TVkdDMyraafKFEVEJDJRzT5191IzGw5MAuLAKHefY2Y/B2a4+3hghJkNBkqBDcCwmvaroCgiIhHxyIIigLtPACZUKLsj5fFtwG37s08FRRERiUi0QTEKCooiIhIJB5IZFhQ10UZERCSkTFFERCKi4VMREZEyCooiIiJAkCmWprsT+0XnFEVERELKFEVEJBLB7cAzK1NUUBQRkYhooo2IiEiZTDunqKAoIiIRybzhU020ERERCSlTFBGRSHgGXpKhoCgiIhFxnJJ0d2K/KCiKiEhEMu+cooKiiIhExHHPrExRE21ERERCyhRFRCQSuqONiIhImcybaKPh0zoycOBA5s+fz4IFC7jlllsq1V911VXMmjWL9957jzfeeIPu3bsDkJWVxaOPPsqsWbOYO3cut956a313vc4s8Vf5u5/CY34yM/y+SvXzfCyPeA+e9DN40s9gjj+eUvcUj3l/HvP+zPOn6rPbdW4N85jC3bzCXSxgSrXtVvAB47mRjSwrV76dYl7kVhbyWtRdjUwRb/EsQ3iGbzCLMdW2+5QpjOZk1jEPgE+YyDguKVtG05/1fFxf3a5Tl18OrVtDz55V17vDiBHQtSv06gUzZ+6tGzMGunULljHVv3wZwElScsBLOtRJUDQzN7O/p6xnmdlaM3uhLvZfyz4MM7P76+v5UsViMR544AHOPfdcevTowdChQ8uC3h5PPPEEvXr14sQTT+See+7hd7/7HQBDhgwhJyeHXr160bt3b6666iqOPPLIdBzG/7d35nFWVFce//5YlDWAC6KAgoqJ45I4UTTGuAXUuBGNGy6RiVskaMZo1HFLUMeFEM0kcSMTUccoxiUCBhA3EpcoEpfgMCoIM4AsogICCkLzmz/ufc2jWUTp169f9/l+Pu9D1a1bxanqqjr3nHvOqY1ipasYx6UczX2cwrO8zZ/40G+t0a8Hfeirp+mrp9lFpwKw1PMZz2BOYDQnMIbxDGapF9T1KdQKZiX/4BH24WwO5hLe5RUWMWeNfitYylSepQPbrrHtv3mUjuy8RnulsJIqXuQXHMKvOIZhTGUsC5i6Rr/lLGESf2RLdqlu24HD6MO99OFevsXPacPWbM5OdSl+rdGvH4wZs+7to0fD5MnpN2QInHtuav/wQxg4EF56CcaPT8vz59eJyAG1ZykuAXaV1DKv9wberaVj13t69uzJlClTmDZtGsuXL2fYsGH06dNntT6LFi2qXm7dujW2AbBN69atadq0KS1btuTTTz/lo48+qlP5a4O5vEJ7utNO3WiqTdiJ7zKV9bwRipjOM3TlAFqoAy3Unq4cwHSeLq3AJWI+02nNFrRmc5rQjM7swRzeWKPfm4xmRw6iCc1Xa5/NRFqxOW3Zqq5ErnXeZxJt6UJbOtOU5mxPb6bz1zX6vcId7MZpNGXTtR5nGmPZnkNKLW7J2H9/2GyzdW8fPhy+/32QYJ99YMECmD0bHn8cevdO+3bokJbXp1zrM87u0y/6Kwe16T4dDRyRl/sC9xc2SOop6QVJr+Z/v5zb+0l6RNIYSZMlDSraZ3HR8nGS7srLR0l6KR/rSUllf3t07tyZGTNWucBmzpxJ586d1+jXv39/pkyZwqBBgzj//PMBeOihh1iyZAmzZ89m+vTpDB48mPkVOCxcwhzasE31ehu2YfFaLKR3eIz7fCCjfAaLnMZNi5lD2w3YtxJYykJa0r56vQXt+YSFq/VZyEw+YQGdiiwkgBUsYwpP82UOrRNZS8XHvEfrIqXeio4sYd5qfT7gLZYwl67st87jTOPJilaKn8W770LXrqvWu3RJbetqr0waqfs0Mww4SVILYHfgpaJtbwL7294DuAq4rmjb14ATgd2AEyUV3Q5r5Tlgn3ysYcDF6+ss6WxJEyRNGDJkyOc6oQ1F0hptBUuwmFtvvZUdd9yRSy65hCuuuAJIVmZVVRXbbLMN3bt358ILL6R79+4lkbOUmDXPt+ZV6cYh9GMCJ2scXdmfJzmveu81WfOaVgZruw4q2rqSNxjOLvRZo99bPM72HECzdVhOlcLa/5qrX4Px/Iq9+PE6jzGPN2hKCzqwQwkkrB+s5RWBtO72yqTyLMVaiz61/Q9J3UhW4qgam9sBd0vqQXpmin1GT9leCCBpErAd1Ig8WJ0uwAOStgY2AaZ9hlxDgII2XNvzutHMnDmTrkVDuy5dujBr1qx19h82bBi33XYbACeffDJjxoxhxYoVzJs3j+eff54999yTadPWe1r1jjZszWJWnfNiZtGaTvCqy/EAABSbSURBVKv1aalVvqRdOJUXuKZ635m8sNq+Xdi3xBKXhmQZrpoPXcoCWvCl6vUVLGMRc3ieWwBYxiJe4vfszRnM5/+YxetMYiTL+QQhmtKM7nyrzs9jY2hNR5Ywt3r9Y96jFVtUry/nY+bzDmPoD8AnfMCTXEQvBrNFnkudyhMN2kqEZAEWOZiYORO22Sa1jxu3evuBB9a1dLWDs6VYSdR29OkIYDBFrtPMNcAztncFjgJaFG1bVrRcxSpFXazAivv/Bvit7d2Ac2psKwsvv/wyPXr0oFu3bjRv3pyTTjqJESNGrNZnxx13rF4+4ogjmDx5MgDTp0/n4IMPBqBVq1bss88+vPnmm3UnfC2xFXuwgKks9P9R5U95m0fpXsMNuMSrXpTTeJwO9ABgWw5iBuNY6gUs9QJmMI5tOahO5a8t2tOVJcxjCR+wkhW8y6tsxarww+a05DCuoTdX0psr6cB27M0ZtKcr+3Fedfv27E8PelWcQgTYgp35iBksYhZVLGcqT9CV/au3b0IbTmYsx/Mox/MoW7LragrRrOR/eYrt6V2uU6gTjj4a7rknWYYvvgjt2sHWW8Ohh8LYsSm4Zv78tHxoZXvUK4razlO8E1hoe6KkA4va27Eq8KbfBh5rrqSdgbeAY4BCpErxsU7fKGlriaqqKgYMGMDjjz9O06ZNufPOO5k0aRIDBw5kwoQJjBw5kgEDBtCrVy+WL1/O/PnzOf30JPott9zC0KFDeeONN5DE0KFDmThxYpnP6PPTRM04wNczgpNYSRX/RF8211d40TfSka+yvQ7jdX7HNI9FNKUF7elFSttooQ7s5Z/wx6xE9+JCWqhDOU/nC9OEpuzGsbzIEMxKtqUnX6ITbzKa9nSlE+uIz29ANKEZ+3ARYzkfs5IeHEUHtucV7mALdmbbIgW5NubwKq3pSFvWnJevJPr2TRbf++8n62/gQFiejaYf/hAOPxxGjUopGa1awdChadtmm8GVV8Jee6X1q65af8BO/aby8hS1trmvz30QabHtNjXaDgQusn2kpG8AdwPzgKeB02x3k9QP2NP2gLzPY8Bg2+MkHQfcSHKlvgG0sd1PUh/gZpJifBHYy/aBNY+1Dpz/n40+50qk+G99fpP3yihJefn1yo7Vy310cxklKR/DfUH18g80voySlI873bNorXG+E9IrUYBLcgHaq7P3y27yL8KfueLvtvesRZE+k1qxFGsqxNw2DhiXl/8GqyUbXZnb7wLuKtrnyKLlh4CH1nLc4cDwtbSvdqwgCIKgvKQ5xU/LLcbnIsq8BUEQBCUiAm2CIAiCoGIJSzEIgiAoEcbhPg2CIAiCmFMMgiAIgiJiTjEIgiAIKpawFIMgCIKSYFaG+zQIgiAIEjGnGARBEASZUIpBEARBAFRm9GkE2gRBEARBJizFIAiCoESspKrCLMVQikEQBEFJqET3aSjFIAiCoEQ4LMUgCIIggMq0FCPQJgiCIKhIJB0m6S1JUyRdup5+x0mypM/8YHFYikEQBEGJKJ37VFJT4BagNzATeFnSCNuTavRrC5wPvLQhxw1LMQiCICgJztGnX/T3GfQEptieavtTYBjQZy39rgEGAUs3ROZQikEQBEFJcLYUN0IpbiFpQtHv7KLDdwZmFK3PzG3VSNoD6Gr7sQ2VOdynQRAEQX3lfdvrmgfUWtpcvVFqAtwM9Ps8/2EoxSAIgqBEmCqWlergM4GuRetdgFlF622BXYFxkgA6ASMkHW17wroOGkoxCIIgKAnGrChdSsbLQA9J3YF3gZOAk6v/b3shsEVhXdI44KL1KUQIpRgEQRCUCLOyZErR9gpJA4DHgabAnbb/W9LVwATbI77IcUMpBkEQBCXBmOUlTN63PQoYVaPtqnX0PXBDjhnRp0EQBEGQCUsxCIIgKAnGLCtdoE1JkO3P7tVwaFQnGwRB8NkI8NrSGzb+yNLGvnP/vp6UjJLQ2CzFkvzhP5cA0tm2h5RbjnIS1yCuAcQ1KFD+61BSW+FxiiJAvwDv15YgG0pjsxTLjqQJdT3yqW/ENYhrAHENCsR1qF9EoE0QBEEQZEIpBkEQBEEmlGLd0+jnUIhrAHENIK5BgbgO9YiYUwyCIAiCTFiKQRAEQZAJpVjHSOqTC9g2WpRL1jd2JJ0q6VvllqOcSNpZUpPGfk9I6icpIlDrAaEU6xBJTYG+wG7llqXMtMgvw36SOpVbmHIgaXPgHuBGSXfm9UaFpObAK8CWbsTzOJLaAXcCl0naJLc16kFCOQmlWIfYriJ946vmhzAbDZK2A+4DbgC+CfxV0iXllaos3A0cZ3tfYBPg0jLLUw4uA/7D9lxofM9CEb8DZgMjbBeqZzcvozyNmsZ6E5YFSXsDX7Y9sqj5FEmXSTqkXHLVFZI6AD8E5pEUwlnAfsDukq6X1KKsAtYRknYHtrL9SG56Nbc3aywWQraMz7F9qaStJbW3vTJva1pm8eoMSV8DdgYuBmbktu2A+yVNkbRDOeVrjIRSrFt2AEZCGhVLOgL4KTAHOF7STuUUrg7oDbQBfmN7uaTmtt8DzgM+JJfha8iKIZ/bn4DlkraQdDjp6+ELba8ouBEbwTzbPSRFAPBzYLikM6Hao9Kg7wNY7V44Dvgq8Imkb5LeCVuQPpy7oHwSNk5CKdYBRW6h1kDBGvox0B+4Hphk+yzbb+f+DXWkvCfwoe2Jeb05gO0Pgf+0/UledwO+Bn2Bl0ju48nAAOAZ4FYASVtJamV7Zb4ODU4xSPoq0NH2vZIOAg4F/kDyGDwkqSNU3wcN7vyLOAn4KzAN6AUsBK4DJpIU4hDbHxQ6N+Bnol4ReYp1iKQnSPNpW5LchjcAi0gvxT1sz6jRv0nBpdQQkPQccJPtR7IreSDwJ9t35Bfl94DZtm8rq6AlQtKmpL/3lrYXShoCTLd9raT2wIFAT2BfYKzt64r2bWt7UTnkrk2ykpsNfIf08v8NMMf2wBxk8jppgLAdcL/tv5dN2BJSdC90BDYDXiQl8b8LzAJOBU4geRH2s31fmURtdISlWEfkOZRmwNHAEUB/2y8AlwPXAoslnSfpqexWpSEpxMy9wD8B2H4JuBC4SNIjwL8CTYGDJD0mqWXxjg3EYmgLnG57YV7/O+mlB/AvwJXAivyF8J0knVa07wBJu9aZpKVje2Ck7VdJltJWwG/ztmtICuFlYDowRNKXinduQNZS4V4ouEdvB35FUo6nAncXReQeJOk5Sb3LIGejI5Ri3bEMaE966M+yPVPSiaSXxH8ARwKdSQ/GVZKubSCKoJjngKMlPZDnT3cnzZ0sBW63fbntE4AqoDtAQTk2BJeq7fdt31/U9A1gSV7eAXgQOFjSzaQBVHNJe0j6C9DX9ht1K3HtY/sd22dJagscAzxt+wNJewB7ARfZfhEYTgpAaiepTSGHrzDfWOkU3wu2pwI/I1mO3wX+1/ZjknYjuZZ/QvKqHAKrDxAbccRuyYgLWkfYXmx7D9s/sv12dpedD/wyjwjnkB6Aibb3Bu7LiuDEHJXYBCrbYrL9hu2ewDiSdXAYcAswleQ2QtJXgF2AOdnFdKOkX0lq2VBeiEX8ERifFcRuJHfhfsBMkvtwIcli6kJSDj/K1wRY9UKspHuiIGt2Bfe3XbASrwKeBCbmczyYNAc/B9iUNJgaX/CiNFD2IT0TN0m6kJSyshfwLHAK8JU8SNxX0pE5UK0QsVsx90C9x3b8yvAjjfrG5OVO+d9TgF8CzfL6t0nzL+3LLW8Jzn9TkjXUA/hzUfudwOV5uU3+/RvwGtClqN+PgO+U+zxq4ToU5vXvAH5U4/q0I70Y78ltp5Gsy/1qHKNTuc9jI867CWlebQiwU27bhRSVeVBe/zrQihSh+RTQtcaxmpb7fGrpmmwKHJ6XbwaOysvNgf8BLsjr1wJvAo8W3zPxq51fWIplwvZY0ui3B/Dr7B76CykA4et5tHwBcL3tBZK2k3RuzmECKjts3/Yy2yuAD4D5kl6Q9BuSS/VeSf9CcivfDTxMci3uBNV5XOcDi8sifC3i/JYDbgKOyNbQkcBy0nzjCcAvcp8HgTOAuyT9p6ReuX2ApO/VpdwbS+G8naJs37N9tpMHpTVJ+c8AnpP0NHAuKWK3Hyk4Z+vCcbK1VEjhqHT3+jLbo/LqfOD6bBkPAlYCt+Ro3Z1Jg+d+wMmSTig+TqW+E+oLzcotQGPGqXrFZEkPkNyI80kupLmk3KU2tn+du3+TlNLwmqR5QAunVIaKjlLN53CqpKNJcyrXkOYTe5OU4SfAUOCfSe4lSHlcz9h+VpKKFEvFYvst4PCsEOfYXinpAmCU7YnZbbYfsDdJSbQHmuXgm+tsfwyVeS/UkPkTkkJ8ljSd8L7tM3PAzW+BbsAV+X7ZH+gl6Re2/2C7qhLPf23YvlrSq8COpNStfnlTL+B1278DyPqvdV4+LG+bXecCNyBCKdYDbD8MPCzpO8BTtj+V9BjJbUiecP8q8Jrtv0nqCVwt6XXgUq+aV6jYF4LtEYVlpXqo7Uju5SWSvk8KPnhdqYD2/iSl2eCw/RiAUnWf9qTqP5DC9o8luVKfyH1ECko6HBiT96+uCuMKmYMtvmfz8i0AklqRonD3Bb5EyvH9A2m+9XKSdf0AqfDFsaQgnWnZk3Aryf1Ykc8DgHPlK0nzbN8v6XSgE/Dr3H4M6f4Ynu+FUaSBRCjFjSDcp/UI26OzQmxHchO1y5v6kV4Ij+b1ubYPI1WAuUepsDLAUTlqtUtdyl0C3iFV8rhN0jCgD8ldCqkKyn225+ZBQMVbiWvD9lLgHGBTSf1I88tbAcU5nLcAD9sek93rZxVcaZWiEGtS7PpzikIdRLr//xXYnFQn9CpgqO0HbL9M8qz0IUV4A7xHivBeKaniB/627ysaAL2ZB4fNSYUgRmdvy00kr8ITkjaRdJhSdZxqIlJ1w6j4G6Yh4pTYPZhkDX6P5E78Cak02BP5349I1TCWO5VM60gKSFjOqjD/isT2ElJN2K4kJTDYKWz/NKCD7Rtyv4q1AjYE2+8DF0hqQ1IKo2wvBpDUjZTP1l2pbu6ZwFjgGEknk5TCvNy3YlzMBTkLMjslrd8nqS/QJCu6KuCFot1OAq6xPUtSf+Bi293y8Vbk41WsFwWqr8u9RUr+p6QPC9yfn5MfAIU81h6kfOh3geclbWl7Xr52FX0d6oJQivUUp0oeR+QAnEOB50kJvlNtnyPph6SaoU/lEeDBpITgh23PL5fctYlThZ8jAbJi+DlwSV5vNA93VoTX1mi+i5TTWkV6Ie4GvGO7r6Tbgc2y+/GDgiKtJIqUY5McjFOc39mE9NL/h6TLgLa2B+Ztl5JK5yHpp8Crtp/MCqFiBgfroqDkgb8BE2x/JOl64L9sz1Aqur8nybN0e16/OgcwXVgYKEEqQF90vCAT5nQ9x/YE2/9uexkpZ61QDWUO6Vt0I0ll4/YGppAelgZFfpktBn5m+yFo+Fbi+lCKWN7V9lWkCNVlpCpJKyW9DRxECr44DPhxVo4VSUGZQfVXRDqTvCbdJf2SNFAamLf/Aphse4Sko0gJ8edIGiHpa5WuEIux/YxTBDvAJFa9yw8nDZBHOtVNbW77XFLhjD9L2haqI7gvkXRAHYte7wmlWFncRgo8GEaK0pydgy6+S7IY/pIj8BpUSHbhZWb73nLLUh+wPZkUbAIph21r21NtX06aexwE/IOU49nR9sdK6TvbFh+nUuaYipRZW1KpwPNICq878McchLI1yWI+K3sVjiN9jeV4kpfliprPhSo8haOIB4FOku4jzb3Ptv0npRSnmyQ9SYpsHwUUyid+nRTZurwcAtdnKuKhCBK259j+LikqdTRprmVXUj3RV5zLgDWkEXGwdvK8K6SgpPmSRko6yPYM278nBWF9DxiW+10KXCfpgMK8VLbCKkYx5GmBb5Pm0oYCb5Hq50Iq+jDGqWTa8SQFenXe1oyU5mJJnSX1ycdrEJ+ocsrzPJZUO/cvpEj2o0gR6/2BG0kDiUKZuG6k9KaJTvWXgyJCKVYgtqfZvtj2ayQXmUkPQ8U/4MHnw/ZC2ycCjwHXSOqZFV1X4AmgKofunwr8O6ng9I8l3SqpXaVFqeb5xcuAo23/m+25edOZwGmSdiRFZd5v+xOlDzp3IZUWhJQEv7ek17LiKJ6/rOiC6051ZS/OEbkfkSK4F5Ei2V8BnnXKh/02KbL9obIJW4+JT0dVOHnUv5PtSeWWJSgvRRbgCknbkKynlaTC2uNtD1eqlNSKlOZwDsnN+Gbe53hSwYihZTmBL4iKcjKzi/jMPN+KpCtJLsNBJHfybqSBQUfSp9vOsv1+HjicAgz0qu99Viz5XriVdJ6bkJTgcaTc1/OAF23fUz4J6y8RfVrh5OixUIhBcWQipApIh5CKjt9ge1EOqjgW2IaU77kC2Murvr5xO6vyQSuGYmvX9nRSHiOSDiS5EIeQ3nUvk6yjG0hu58Wkr923IXlcniRV06l48r1wtqQtSdG4i2zPzuk6SwkrcZ2EUgyCBojtByVNB6Zlhbg7yTqcTJqPvgv4FilqFUk3kOaY/lAeiWuHGmkXz5G+0fk3knv1Bdvn5Zzeu0jpTUsk/YA0QHjBq75v2CDIKRg/A1CqhPUN4PfOZQGDNQmlGAQNFKcPOReYRZpbu8b2VKXSaYtsj85Koj8praeiKQ4yy9bSUwCS5pLSOL5BKrq/ABglaStSCsNjpC9PNFhsj5f0S9sNLm2rNgmlGASNg4XA26QvskwlFRUvfGXjDuBB2/9TLuFKje0HJC0DziJ9ZeJJ209Lupj0/cqXnAr0N2hCIX42EWgTBI0ISduTvtEo22dki3E46VuVy9a/d2VSs5JNnmt8j1T15SLgrqJE+KCRE0oxCBohklrmlIXZwM22B5VbplJTszRgjrb9CvA723PKJ1lQnwilGASNkIL1JKm/7VvLLU+5kNQqgk6CYkIpBkHQ6GgIxcGD0hBKMQiCIAgyUeYtCIIgCDKhFIMgCIIgE0oxCIIgCDKhFIMgCIIgE0oxCIIgCDKhFIMgCIIgE0oxCIIgCDL/D4fvccyL9z0GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 576x576 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "regions = [\"Visual\", \"Goal\", \"Retrieval\", \"Imaginal\", \"Manual\"] \n",
    "fig, ax = plt.subplots()\n",
    "ax.set_title(\"Functional Connectivity\")\n",
    "im, cbar = heatmap(mean_correlation, regions, regions, ax= ax, cmap=\"gnuplot\")\n",
    "texts = annotate_heatmap(im, valfmt=\"{x:.2f}\")\n",
    "fig.set_size_inches(8, 6)\n",
    "fig.savefig(\"func_conn.png\", dpi=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PGM Inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import combinations\n",
    "from pgmpy.models import BayesianModel\n",
    "from pgmpy.estimators import MaximumLikelihoodEstimator, BayesianEstimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"PCA_results/sub-001.csv\")\n",
    "df.columns = [i[:-7] for i in df.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = df.columns\n",
    "X = df.values\n",
    "scaler = MinMaxScaler()\n",
    "X = scaler.fit_transform(X)\n",
    "df = pd.DataFrame(X, columns=labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add nodes and edges\n",
    "model = BayesianModel(list(combinations(labels, r=2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.fit(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}