{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finishing the encoder\n",
    "\n",
    "__Contents:__\n",
    "\n",
    "- <a href=\"#raw2feat\">Create features from raw stimulus</a>\n",
    "- <a href=\"#algorun\">Initialize model and execute algorithm</a>\n",
    "- <a href=\"#eval\">Evaluate performance</a>\n",
    "\n",
    "___\n",
    "\n",
    "To the end of putting all the pieces together to get a fully-functional learning machine, which takes advantage of the various functions and classes we prepared in previous lessons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic workflow of the final procedure is as follows:\n",
    "\n",
    "    Load raw stimulus --> Generate and save visual features --> Fit sparse linear model.\n",
    "\n",
    "We shall take this one piece at a time. After going through a simple and fast prototypical implementation of this encoder, we will put forward a series of exercises that constitute the bulk of the work to be done here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"raw2feat\"></a>\n",
    "## Create features using filter bank\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import tables\n",
    "\n",
    "import gaborfil\n",
    "import models\n",
    "import dataclass\n",
    "import algorithms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/vim-2/stimulus_ds.h5 (File) 'vim-2: stimulus'\n",
      "Last modif.: 'Tue Mar 27 21:14:47 2018'\n",
      "Object Tree: \n",
      "/ (RootGroup) 'vim-2: stimulus'\n",
      "/test (Array(64, 64, 3, 8100)) 'Testing data'\n",
      "/train (Array(64, 64, 3, 108000)) 'Training data'\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Establish connection with the file objects.\n",
    "h5_X = tables.open_file(\"data/vim-2/stimulus_ds.h5\", mode=\"r\")\n",
    "print(h5_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the parameters that specify the first filter bank.\n",
    "PIX_W = 64\n",
    "PIX_H = 64\n",
    "max_cycles = 32 # the maximum cycles per image.\n",
    "myparas = {\"freqs\": max_cycles/max(PIX_W,PIX_H),\n",
    "           \"dir\": 0,\n",
    "           \"amp\": 0.1,\n",
    "           \"sdev\": max(PIX_W,PIX_H)/20,\n",
    "           \"phase\": 0}\n",
    "mygrid_h = 4\n",
    "mygrid_w = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Images processed so far: 0\n",
      "Images processed so far: 2160\n",
      "Images processed so far: 4320\n",
      "Images processed so far: 6480\n",
      "Images processed so far: 8640\n",
      "Images processed so far: 10800\n",
      "Images processed so far: 12960\n",
      "Images processed so far: 15120\n",
      "Images processed so far: 17280\n",
      "Images processed so far: 19440\n",
      "Images processed so far: 21600\n",
      "Images processed so far: 23760\n",
      "Images processed so far: 25920\n",
      "Images processed so far: 28080\n",
      "Images processed so far: 30240\n",
      "Images processed so far: 32400\n",
      "Images processed so far: 34560\n",
      "Images processed so far: 36720\n",
      "Images processed so far: 38880\n",
      "Images processed so far: 41040\n",
      "Images processed so far: 43200\n",
      "Images processed so far: 45360\n",
      "Images processed so far: 47520\n",
      "Images processed so far: 49680\n",
      "Images processed so far: 51840\n",
      "Images processed so far: 54000\n",
      "Images processed so far: 56160\n",
      "Images processed so far: 58320\n",
      "Images processed so far: 60480\n",
      "Images processed so far: 62640\n",
      "Images processed so far: 64800\n",
      "Images processed so far: 66960\n",
      "Images processed so far: 69120\n",
      "Images processed so far: 71280\n",
      "Images processed so far: 73440\n",
      "Images processed so far: 75600\n",
      "Images processed so far: 77760\n",
      "Images processed so far: 79920\n",
      "Images processed so far: 82080\n",
      "Images processed so far: 84240\n",
      "Images processed so far: 86400\n",
      "Images processed so far: 88560\n",
      "Images processed so far: 90720\n",
      "Images processed so far: 92880\n",
      "Images processed so far: 95040\n",
      "Images processed so far: 97200\n",
      "Images processed so far: 99360\n",
      "Images processed so far: 101520\n",
      "Images processed so far: 103680\n",
      "Images processed so far: 105840\n",
      "(108000, 16)\n"
     ]
    }
   ],
   "source": [
    "# Construct features using the specified filter bank (TRAINING).\n",
    "X_tr = gaborfil.G2_getfeatures(ims=h5_X.root.train.read(),\n",
    "                          fil_paras=myparas,\n",
    "                          gridshape=(mygrid_h, mygrid_w),\n",
    "                          mode=\"reflect\", cval=0, verbose=True)\n",
    "print(X_tr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Images processed so far: 0\n",
      "Images processed so far: 162\n",
      "Images processed so far: 324\n",
      "Images processed so far: 486\n",
      "Images processed so far: 648\n",
      "Images processed so far: 810\n",
      "Images processed so far: 972\n",
      "Images processed so far: 1134\n",
      "Images processed so far: 1296\n",
      "Images processed so far: 1458\n",
      "Images processed so far: 1620\n",
      "Images processed so far: 1782\n",
      "Images processed so far: 1944\n",
      "Images processed so far: 2106\n",
      "Images processed so far: 2268\n",
      "Images processed so far: 2430\n",
      "Images processed so far: 2592\n",
      "Images processed so far: 2754\n",
      "Images processed so far: 2916\n",
      "Images processed so far: 3078\n",
      "Images processed so far: 3240\n",
      "Images processed so far: 3402\n",
      "Images processed so far: 3564\n",
      "Images processed so far: 3726\n",
      "Images processed so far: 3888\n",
      "Images processed so far: 4050\n",
      "Images processed so far: 4212\n",
      "Images processed so far: 4374\n",
      "Images processed so far: 4536\n",
      "Images processed so far: 4698\n",
      "Images processed so far: 4860\n",
      "Images processed so far: 5022\n",
      "Images processed so far: 5184\n",
      "Images processed so far: 5346\n",
      "Images processed so far: 5508\n",
      "Images processed so far: 5670\n",
      "Images processed so far: 5832\n",
      "Images processed so far: 5994\n",
      "Images processed so far: 6156\n",
      "Images processed so far: 6318\n",
      "Images processed so far: 6480\n",
      "Images processed so far: 6642\n",
      "Images processed so far: 6804\n",
      "Images processed so far: 6966\n",
      "Images processed so far: 7128\n",
      "Images processed so far: 7290\n",
      "Images processed so far: 7452\n",
      "Images processed so far: 7614\n",
      "Images processed so far: 7776\n",
      "Images processed so far: 7938\n",
      "(8100, 16)\n"
     ]
    }
   ],
   "source": [
    "# Construct features using the specified filter bank (TESTING).\n",
    "X_te = gaborfil.G2_getfeatures(ims=h5_X.root.test.read(),\n",
    "                          fil_paras=myparas,\n",
    "                          gridshape=(mygrid_h, mygrid_w),\n",
    "                          mode=\"reflect\", cval=0, verbose=True)\n",
    "print(X_te.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the above features only cover one specific filter setting, namely that with spatial frequency set to zero. Once again linking up with Nishimoto *et al.* (2011), the filter bank that we are using here corresponds to their \"static\" model (i.e., motion information is not considered). The filter orientations they used were \"0, 45, 90 and 135 degrees\" (quoting their appendix). We've covered the 0 degree case, so let's handle the rest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding features using dir = 45.0 degrees\n",
      "Adding features using dir = 90.0 degrees\n",
      "Adding features using dir = 135.0 degrees\n"
     ]
    }
   ],
   "source": [
    "todo_dir = math.pi * np.array([1,2,3]) / 4\n",
    "\n",
    "for mydir in todo_dir:\n",
    "    print(\"Adding features using dir =\", mydir*(360/(2*math.pi)), \"degrees\")\n",
    "    myparas[\"dir\"] = mydir\n",
    "    \n",
    "    tmp_X = gaborfil.G2_getfeatures(ims=h5_X.root.train.read(),\n",
    "                               fil_paras=myparas,\n",
    "                               gridshape=(mygrid_h, mygrid_w),\n",
    "                               mode=\"reflect\", cval=0, verbose=False)\n",
    "    X_tr = np.concatenate((X_tr, tmp_X), axis=1)\n",
    "    \n",
    "    tmp_X = gaborfil.G2_getfeatures(ims=h5_X.root.test.read(),\n",
    "                               fil_paras=myparas,\n",
    "                               gridshape=(mygrid_h, mygrid_w),\n",
    "                               mode=\"reflect\", cval=0, verbose=False)\n",
    "    X_te = np.concatenate((X_te, tmp_X), axis=1)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<closed File>\n"
     ]
    }
   ],
   "source": [
    "# All finished with this file, so close it.\n",
    "h5_X.close()\n",
    "print(h5_X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's do some temporal down-scaling. The responses only have one observation per second, while there are 15 frames of stimulus per second. Let's take averages over disjoint 15-frame windows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7200\n",
      "540\n"
     ]
    }
   ],
   "source": [
    "framerate = 15\n",
    "n_tr = X_tr.shape[0]//framerate\n",
    "n_te = X_te.shape[0]//framerate\n",
    "print(n_tr)\n",
    "print(n_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_tr shape after down-sampling: (7200, 64)\n",
      "X_te shape after down-sampling: (540, 64)\n"
     ]
    }
   ],
   "source": [
    "# Training data.\n",
    "tmp_X = np.zeros((n_tr,X_tr.shape[1]), dtype=X_tr.dtype)\n",
    "idx = np.arange(framerate)\n",
    "for i in range(n_tr):\n",
    "    tmp_X[i,:] = np.mean(X_tr[idx,:], axis=0)\n",
    "    idx += framerate\n",
    "X_tr = tmp_X\n",
    "print(\"X_tr shape after down-sampling:\", X_tr.shape)\n",
    "\n",
    "# Testing data.\n",
    "tmp_X = np.zeros((n_te,X_te.shape[1]), dtype=X_te.dtype)\n",
    "idx = np.arange(framerate)\n",
    "for i in range(n_te):\n",
    "    tmp_X[i,:] = np.mean(X_te[idx,:], axis=0)\n",
    "    idx += framerate\n",
    "X_te = tmp_X\n",
    "print(\"X_te shape after down-sampling:\", X_te.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again following Nishimoto et al. (2011), let us compute the so-called Z-score of the feature values. This is simply\n",
    "\n",
    "\\begin{align}\n",
    "z = \\frac{x - \\bar{x} }{\\sqrt{\\widehat{v}}},\n",
    "\\end{align}\n",
    "\n",
    "where $\\bar{x}$ is the empirical mean, and $\\widehat{v}$ is the empirical variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Z-scores\n",
    "Z_tr = X_tr - np.mean(X_tr, axis=0)\n",
    "Z_tr = Z_tr / np.std(Z_tr, axis=0)\n",
    "#print(\"Mean =\", np.mean(Z_tr, axis=0), \"StdDev =\", np.std(Z_tr, axis=0))\n",
    "Z_te = X_te - np.mean(X_te, axis=0)\n",
    "Z_te = Z_te / np.std(Z_te, axis=0)\n",
    "#print(\"Mean =\", np.mean(Z_te, axis=0), \"StdDev =\", np.std(Z_te, axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, a hard truncation of outliers is carried out (anything beyond three standard deviations from the mean is truncated). This is done separately for training and testing data, to ensure the learner does not gain unfair oracle information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Truncation of outliers.\n",
    "thres = 3\n",
    "for j in range(X_tr.shape[1]):\n",
    "    stdval = np.std(Z_tr[:,j])\n",
    "    Z_tr[:,j] = np.clip(Z_tr[:,j], a_min=(-thres*stdval), a_max=thres*stdval)\n",
    "    stdval = np.std(Z_te[:,j])\n",
    "    Z_te[:,j] = np.clip(Z_te[:,j], a_min=(-thres*stdval), a_max=thres*stdval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a new hierarchical data file, `features.h5`, to store the features that shall be used for both training and evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "Last modif.: 'Sat Apr  7 17:14:23 2018'\n",
      "Object Tree: \n",
      "/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Open file connection, writing new file to disk.\n",
    "myh5 = tables.open_file(\"data/vim-2/features.h5\",\n",
    "                        mode=\"w\",\n",
    "                        title=\"Features from vim-2 stimulus, via 2D Gabor filter bank\")\n",
    "print(myh5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "Last modif.: 'Sat Apr  7 17:14:34 2018'\n",
      "Object Tree: \n",
      "/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "/train (Array(7200, 64)) 'Training data'\n",
      "\n",
      "data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "Last modif.: 'Sat Apr  7 17:14:34 2018'\n",
      "Object Tree: \n",
      "/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "/test (Array(540, 64)) 'Testing data'\n",
      "/train (Array(7200, 64)) 'Training data'\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Add arrays.\n",
    "myh5.create_array(where=myh5.root, name=\"train\", obj=Z_tr, title=\"Training data\")\n",
    "print(myh5)\n",
    "myh5.create_array(where=myh5.root, name=\"test\", obj=Z_te, title=\"Testing data\")\n",
    "print(myh5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<closed File>\n"
     ]
    }
   ],
   "source": [
    "# Close the file connection.\n",
    "myh5.close()\n",
    "print(myh5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"algorun\"></a>\n",
    "## Initialize model and execute algorithm\n",
    "\n",
    "Assuming the above has all been run once, then we can feel free to restart the kernel and run from here.\n",
    "\n",
    "We shall run the `Algo_LASSO_CD` routine implemented in a previous lesson here, over a grid of $\\lambda$ parameters controlling the impact of the $\\ell_{1}$ norm constraint. As for performance metrics, citing the Nishimoto et al. (2011) work, from which this data set was born:\n",
    "\n",
    "> *\"Prediction accuracy was defined as the correlation between predicted and observed BOLD signals. The averaged accuracy across subjects and voxels in early visual areas (V1, V2, V3, V3A, and V3B) was 0.24, 0.39, and 0.40 for the static, nondirectional, and directional encoding models, respectively.\"*\n",
    "\n",
    "Every element of our implementation is less sophisticated than theirs, from the filter bank used to create inputs, to the learning procedure used to set parameter values, and thus this performance should be considered an upper bound for the performance we can achieve in our pedagogical exercise here. In particular, since we are only using two-dimensional Gabor filters, our encoding model corresponds to a simplified version of their \"static\" encoding model (which achieved average accuracy of 0.24).\n",
    "\n",
    "After running our algorithm, as an output we get an estimate $\\widehat{w}$ called `w_est`. Given a new collection of features $X$ (this will be `X_te`) and a response $y$ (this will be `y_te`), the goal is to estimate as $\\widehat{y} \\approx y$ with $\\widehat{y} = X\\widehat{w}$. Evaluation of performance then can be done with the correlation\n",
    "\n",
    "\\begin{align}\n",
    "\\text{corr}\\,(\\widehat{y},y) = \\frac{\\text{cov}\\,(\\widehat{y},y)}{\\sqrt{\\text{var}\\,(\\widehat{y})\\text{var}\\,(y)}},\n",
    "\\end{align}\n",
    "\n",
    "implemented in `scipy.stats.pearsonr`, and the *root mean squared error* (RMSE), defined\n",
    "\n",
    "\\begin{align}\n",
    "\\text{RMSE}\\,(\\widehat{y},y) = \\left( \\frac{1}{m} \\sum_{i=1}^{m} (\\widehat{y}_{i}-y_{i})^2 \\right)^{1/2},\n",
    "\\end{align}\n",
    "\n",
    "implemented in `mod.eval` as a method of the model object, where $m$ represents the number of samples in the test data (here $m=540$).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import tables\n",
    "\n",
    "import gaborfil\n",
    "import dataclass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "Last modif.: 'Sat Apr  7 17:14:42 2018'\n",
      "Object Tree: \n",
      "/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
      "/test (Array(540, 64)) 'Testing data'\n",
      "/train (Array(7200, 64)) 'Training data'\n",
      "\n",
      "data/vim-2/response.h5 (File) 'vim-2: BOLD responses'\n",
      "Last modif.: 'Mon Apr  9 09:14:08 2018'\n",
      "Object Tree: \n",
      "/ (RootGroup) 'vim-2: BOLD responses'\n",
      "/sub1 (Group) 'Data for subject 1'\n",
      "/sub2 (Group) 'Data for subject 2'\n",
      "/sub3 (Group) 'Data for subject 3'\n",
      "/sub3/idx (Group) 'ROI-specific voxel indices'\n",
      "/sub3/idx/v1lh (Array(653,)) ''\n",
      "/sub3/idx/v1rh (Array(713,)) ''\n",
      "/sub3/idx/v2lh (Array(735,)) ''\n",
      "/sub3/idx/v2rh (Array(642,)) ''\n",
      "/sub3/idx/v3alh (Array(164,)) ''\n",
      "/sub3/idx/v3arh (Array(118,)) ''\n",
      "/sub3/idx/v3blh (Array(88,)) ''\n",
      "/sub3/idx/v3brh (Array(138,)) ''\n",
      "/sub3/idx/v3lh (Array(504,)) ''\n",
      "/sub3/idx/v3rh (Array(627,)) ''\n",
      "/sub3/resp (Group) 'Response arrays'\n",
      "/sub3/resp/test (Array(4381, 540)) 'Testing data'\n",
      "/sub3/resp/train (Array(4381, 7200)) 'Training data'\n",
      "/sub2/idx (Group) 'ROI-specific voxel indices'\n",
      "/sub2/idx/v1lh (Array(470,)) ''\n",
      "/sub2/idx/v1rh (Array(573,)) ''\n",
      "/sub2/idx/v2lh (Array(733,)) ''\n",
      "/sub2/idx/v2rh (Array(926,)) ''\n",
      "/sub2/idx/v3alh (Array(135,)) ''\n",
      "/sub2/idx/v3arh (Array(202,)) ''\n",
      "/sub2/idx/v3blh (Array(83,)) ''\n",
      "/sub2/idx/v3brh (Array(140,)) ''\n",
      "/sub2/idx/v3lh (Array(714,)) ''\n",
      "/sub2/idx/v3rh (Array(646,)) ''\n",
      "/sub2/resp (Group) 'Response arrays'\n",
      "/sub2/resp/test (Array(4622, 540)) 'Testing data'\n",
      "/sub2/resp/train (Array(4622, 7200)) 'Training data'\n",
      "/sub1/idx (Group) 'ROI-specific voxel indices'\n",
      "/sub1/idx/v1lh (Array(490,)) ''\n",
      "/sub1/idx/v1rh (Array(504,)) ''\n",
      "/sub1/idx/v2lh (Array(715,)) ''\n",
      "/sub1/idx/v2rh (Array(762,)) ''\n",
      "/sub1/idx/v3alh (Array(92,)) ''\n",
      "/sub1/idx/v3arh (Array(160,)) ''\n",
      "/sub1/idx/v3blh (Array(104,)) ''\n",
      "/sub1/idx/v3brh (Array(152,)) ''\n",
      "/sub1/idx/v3lh (Array(581,)) ''\n",
      "/sub1/idx/v3rh (Array(560,)) ''\n",
      "/sub1/resp (Group) 'Response arrays'\n",
      "/sub1/resp/test (Array(4120, 540)) 'Testing data'\n",
      "/sub1/resp/train (Array(4120, 7200)) 'Training data'\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Open file connections with data to be used in learning and evaluation.\n",
    "h5_X = tables.open_file(\"data/vim-2/features.h5\", mode=\"r\")\n",
    "print(h5_X)\n",
    "h5_y = tables.open_file(\"data/vim-2/response.h5\", mode=\"r\")\n",
    "print(h5_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Subject ID.\n",
    "subid = 1\n",
    "y_node = h5_y.get_node(h5_y.root, \"sub\"+str(subid))\n",
    "\n",
    "# Number of voxels.\n",
    "num_voxels = y_node.resp.train.nrows\n",
    "\n",
    "# Set up the model and data objects.\n",
    "mod = models.LinearL1()\n",
    "data = dataclass.DataSet()\n",
    "data.X_tr = h5_X.root.train.read()\n",
    "data.X_te = h5_X.root.test.read()\n",
    "\n",
    "# Basic info.\n",
    "n = data.X_tr.shape[0]\n",
    "d = data.X_tr.shape[1]\n",
    "\n",
    "# Dictionaries of performance over all voxels.\n",
    "dict_corr_tr = {}\n",
    "dict_corr_te = {}\n",
    "dict_l0norm = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the long routine: run the full learning procedure for each individual voxel (using random initial values each time)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Voxel: 1 of 4120\n",
      "Voxel: 101 of 4120\n",
      "Voxel: 201 of 4120\n",
      "Voxel: 301 of 4120\n",
      "Voxel: 401 of 4120\n",
      "Voxel: 501 of 4120\n",
      "Voxel: 601 of 4120\n",
      "Voxel: 701 of 4120\n",
      "Voxel: 801 of 4120\n",
      "Voxel: 901 of 4120\n",
      "Voxel: 1001 of 4120\n",
      "Voxel: 1101 of 4120\n",
      "Voxel: 1201 of 4120\n",
      "Voxel: 1301 of 4120\n",
      "Voxel: 1401 of 4120\n",
      "Voxel: 1501 of 4120\n",
      "Voxel: 1601 of 4120\n",
      "Voxel: 1701 of 4120\n",
      "Voxel: 1801 of 4120\n",
      "Voxel: 1901 of 4120\n",
      "Voxel: 2001 of 4120\n",
      "Voxel: 2101 of 4120\n",
      "Voxel: 2201 of 4120\n",
      "Voxel: 2301 of 4120\n",
      "Voxel: 2401 of 4120\n",
      "Voxel: 2501 of 4120\n",
      "Voxel: 2601 of 4120\n",
      "Voxel: 2701 of 4120\n",
      "Voxel: 2801 of 4120\n",
      "Voxel: 2901 of 4120\n",
      "Voxel: 3001 of 4120\n",
      "Voxel: 3101 of 4120\n",
      "Voxel: 3201 of 4120\n",
      "Voxel: 3301 of 4120\n",
      "Voxel: 3401 of 4120\n",
      "Voxel: 3501 of 4120\n",
      "Voxel: 3601 of 4120\n",
      "Voxel: 3701 of 4120\n",
      "Voxel: 3801 of 4120\n",
      "Voxel: 3901 of 4120\n",
      "Voxel: 4001 of 4120\n",
      "Voxel: 4101 of 4120\n"
     ]
    }
   ],
   "source": [
    "for voxidx in range(num_voxels):\n",
    "    \n",
    "    # Set up the responses.\n",
    "    data.y_tr = np.transpose(np.take(a=y_node.resp.train.read(),\n",
    "                                     indices=[voxidx],\n",
    "                                     axis=0))\n",
    "    data.y_te = np.transpose(np.take(a=y_node.resp.test.read(),\n",
    "                                     indices=[voxidx],\n",
    "                                     axis=0))\n",
    "\n",
    "    # Set up for a loop over trials and lambda values.\n",
    "    todo_lambda = np.logspace(start=math.log10(1/n), stop=math.log10(2.5), num=50)\n",
    "    num_loops = 15\n",
    "    t_max = num_loops * d\n",
    "    \n",
    "    # Storage for performance metrics.\n",
    "    corr_tr = np.zeros(todo_lambda.size, dtype=np.float32)\n",
    "    corr_te = np.zeros(todo_lambda.size, dtype=np.float32)\n",
    "    l0norm = np.zeros(todo_lambda.size, dtype=np.uint32)\n",
    "    \n",
    "    # Initialize and run learning algorithm.\n",
    "    w_init = 1*np.random.uniform(size=(d,1))\n",
    "\n",
    "    for l in range(todo_lambda.size):\n",
    "\n",
    "        lamval = todo_lambda[l]\n",
    "        \n",
    "        if (voxidx % 100 == 0) and (l == 0):\n",
    "            print(\"Voxel:\", voxidx+1, \"of\", num_voxels)\n",
    "            #print(\"Lambda value =\", lamval, \"(\", l, \"of\", todo_lambda.size, \")\")\n",
    "\n",
    "        # Use warm starts when available.\n",
    "        if l > 0:\n",
    "            w_init = al.w\n",
    "        \n",
    "        al = algorithms.Algo_CDL1(w_init=w_init, t_max=t_max, lamreg=lamval)\n",
    "\n",
    "        # Iterate the learning algorithm.\n",
    "        for onestep in al:\n",
    "            al.update(model=mod, data=data)\n",
    "\n",
    "        # Record performance based on final output.\n",
    "        corr_tr[l] = gaborfil.corr(w=al.w, X=data.X_tr, y=data.y_tr)\n",
    "        corr_te[l] = gaborfil.corr(w=al.w, X=data.X_te, y=data.y_te)\n",
    "        l0norm[l] = np.nonzero(al.w)[0].size\n",
    "    \n",
    "    # Save the performance for this voxel.\n",
    "    dict_corr_tr[voxidx] = corr_tr\n",
    "    dict_corr_te[voxidx] = corr_te\n",
    "    dict_l0norm[voxidx] = l0norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save results to disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "# Method name\n",
    "mthname = \"DefaultMth\"\n",
    "\n",
    "# Lambda values used.\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".lam\"\n",
    "with open(fname, mode=\"bw\") as fbin:\n",
    "    pickle.dump(todo_lambda, fbin)\n",
    "\n",
    "# Correlation over lambda values.\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrtr\"\n",
    "with open(fname, mode=\"bw\") as fbin:\n",
    "    pickle.dump(dict_corr_tr, fbin)\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrte\"\n",
    "with open(fname, mode=\"bw\") as fbin:\n",
    "    pickle.dump(dict_corr_te, fbin)\n",
    "\n",
    "# Sparsity over lambda values.\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".l0norm\"\n",
    "with open(fname, mode=\"bw\") as fbin:\n",
    "    pickle.dump(dict_l0norm, fbin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<a id=\"eval\"></a>\n",
    "## Evaluate performance\n",
    "\n",
    "Feel free to restart the kernel here for evaluation, since the core results should all be saved to disk at this point. First let's load the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preparation.\n",
    "import math\n",
    "import numpy as np\n",
    "import pickle\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Method name\n",
    "mthname = \"DefaultMth\"\n",
    "\n",
    "# Subject ID.\n",
    "subid = 1\n",
    "\n",
    "# Lambda values used.\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".lam\"\n",
    "with open(fname, mode=\"br\") as fbin:\n",
    "    todo_lambda = pickle.load(fbin)\n",
    "\n",
    "# Correlation over lambda values, on the training and test data.\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrtr\"\n",
    "with open(fname, mode=\"br\") as fbin:\n",
    "    dict_corr_tr = pickle.load(fbin)\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrte\"\n",
    "with open(fname, mode=\"br\") as fbin:\n",
    "    dict_corr_te = pickle.load(fbin)\n",
    "\n",
    "# Sparsity over lambda values.\n",
    "fname = \"results/\"+mthname+\"sub\"+str(subid)+\".l0norm\"\n",
    "with open(fname, mode=\"br\") as fbin:\n",
    "    dict_l0norm = pickle.load(fbin)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, some visualization of performance for each voxel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results: subject 1 , voxel id 3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzwAAAG9CAYAAADQhXO6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl4lNXZx/HvnRDIyr6ooIKCW1ERAXdwF7TiUndt1aqgrbZvUau0al1bbKtSFTdcW+v2ulSqaFFfqRsiYKmCqICiIAqRBEhCWBLO+8d5RsYwgUmYmeeZye9zXblmMvPMzJ0oM/k955z7mHMOERERERGRXJQXdgEiIiIiIiLposAjIiIiIiI5S4FHRERERERylgKPiIiIiIjkLAUeERERERHJWQo8IiIiIiKSsxR4JOeY2QIzO7yZjz3IzD5JdU3pYmYHmNlcM6s2s+PNrJuZvWFmVWZ2i5n9xszuT+J57jGzqzNRs4iINM7MZpvZwRl4nWa/72/J56xIGBR4JOXM7Awzmx78Ef61mb1kZgeGXVciZubMrHfse+fcm865ncOsqYmuB+50zpU65/4BjAC+Bdo65y51zv3eOXf+5p7EOXehc+6GLS3GzA42s0Vb+jwiIplmZgea2TtmtsLMKszsbTMbmOk6nHM/cM5NDmq61sweTdPrpOR9H8DMfmVm3wS/uwfNrE0qnlckVRR4JKXMbBQwFvg90A3YDrgLOK4Zz9UqmdtauO2B2Q2+/8hpR2ERkaSZWVvgBeAOoCPQHbgOWJOG18qpzzEzOwq4EjgM6AnsgP/dpft1zcz0d6wkRf+jSMqYWTv8iMPPnXPPOudqnHPrnHP/dM5dHhzTxszGmtni4Gts7ExQbHTAzK4ws2+AhxLdFhz7QzObaWbLgzNyezRS0yAzmxIc97WZ3WlmrYP73ggO+28wGnVqwxEKM9vVzCYHj59tZsPj7nvYzMaZ2YvBFLKpZrbjJn4/sbOHy81soZmdE/u9mdlfzazczL4ws6vi38TN7KdmNsfMKs3sX2a2fXD7fPwHyz+D+h8HzgZ+HXx/eMOzg5uo4WEzuzHuuEZ/v8FUhsvM7IPgbN6TZlZoZiXAS8A2wetXm9k2jf0+REQiZCcA59zjzrl651ytc26Sc+4DADM7JxjxuSN43/vYzA6LPdjMzg3ep6vM7DMzGxl3X6LPts5m9kLwHlthZm/G3veD99jDzWwo8Bvg1OD99L9mdrKZzYgv3MwuNbN/NPyBzOw0M5ve4LZfmdmE4Pp37/tm1iGopzz4rHnBzHok+bs7G3jAOTfbOVcJ3ACc09jBm/vsNLP9zWxa8HueZmb7x9032cxuMrO3gVXADsFtNwafVdVm9k8z62RmfzezlcFz9EzyZ5EcpcAjqbQfUAg8t4ljfgvsC/QD9gQGAVfF3b8V/uza9vjpWRvdZmb9gQeBkUAn4F5ggiUeQq8HfgV0Duo7DPgZgHNucHDMnsGUsCfjH2hmBcA/gUlAV+AS4O9mFj/l7XT8mawOwDzgpkQ/tJlthw8DdwBdgp9/ZnD3HUA7fHgZAvwEODd43PH4D7wTg8e9CTwe1L8j8CVwbFD/6cDfgT8G37/ahBrij0vm93sKMBToBewBnOOcqwGGAYuD1y91zi1O9PsQEYmYT4F6M3vEzIaZWYcEx+wDfIb/PPkd8KyZdQzuWwr8EGiLf/++LXgvjWn42XYpsAj/XtwN/z7/vZF559zL+NkSTwbvp3sCE4BeZrZr3KFnAX9LUO8EYGcz6xN32xnAYwmOzcOfUNwePzOjFrgzwXGJ/AD4b9z3/wW6mVmnTTwm4Wdn8Pt8Ebgd//lzK/Big+f6Mf53WAZ8Edx2WnB7d2BHYErw83QE5uD/e0kLpsAjqdQJ+NY5V7eJY84ErnfOLXXOlePf8H4cd/964HfOuTXOudpGbrsAuNc5NzU4E/cIftrBvg1fzDk3wzn3rnOuzjm3AP/H+5Akf559gVJgjHNurXPu//BTHk6PO+ZZ59x7wc/8d3yIaOznfjU4e7jOObfMOTfTzPKBU4HRzrmqoMZb4n4nI4E/OOfmBK/xe6BfbJSniRLWkOC4ZH6/tzvnFjvnKvChsLGfW0Qk8pxzK4ED8aFjPFBuZhPMrFvcYUuBscH755PAJ8AxweNfdM7Nd96/8SfKDop7bMPPsXXA1sD2wfO9mcxUZOfcGuBJfMjBzH6An0b2QoJjVwHPE3xmBcFnF3wQanjsMufcM865Vc65KnwASfazshRYEfd97HrZJh7T2GfnMcBc59zfgs/tx4GPgWPjHvtwMJpU55xbF9z2UPD7X4E/sTffOfdq8Pz/C+yV5M8iOUqBR1JpGdDZNj0/eRs2nJEhuB4/7ancObe6wWMa3rY9cGkwFWC5mS0Htm3wPACY2U7B0Pw3ZrYSHxg6J/nzbAMsdM6tb1Bv97jvv4m7vgr/xp/ItsD8BLd3Blqz8e8k9hrbA3+J+zkrAGtQQ7Iaq6GhZH6/yf7cIiJZITixdI5zrgfQF/+eNzbukK8ahJLvPr+CUaF3g+lpy4Gj+f5nTcPPsT/hRzYmBVPgrmxCqY8AZ5iZ4U+OPRUEoUQeY8NJujOAfwRB6HvMrNjM7jU/rXol8AbQPjgptznV+JGtmNj1KvOdQmNTnO+JO6axz5CGfyPAxp+7CxPUsCTuem2C7/UZ1cIp8EgqTQFWA8dv4pjF+D+oY7YLbotJdIar4W0LgZucc+3jvoqDM0EN3Y0/O9THOdcWP23ANvNzxNe6rX1/UeR2wFdJPr5hzYnW93yLP9PX8HcSe42FwMgGP2uRc+6dFNaQ6Lhkf78NqVmCiGQ959zHwMP44BPTPQgZMdsBi4Ppvs8Afwa6OefaAxP5/mdNw+lqVc530twBP3oxKn5NUGOPCx77LrAWP4J0Bomns8VMwp+I7IcPPomms4GfYrczsE/wWRmb8p3M5+Vs/BT1mD2BJcGo0e/jpjhfmMRzNfwbATb+3NXnjDSZAo+kTDCUfA0wzvyeMMVmVhCc+fpjcNjjwFVm1sXMOgfHN7Xl5njgQjPbx7wSMzvGzBINn5cBK4FqM9sFuKjB/Uvwa2cSmQrU4JsAFJjfF+FY4Ikm1gt+yP5wMzvFzFoFCyr7OefqgaeAm8ysLJiqNooNv5N7gNHBtIVYg4OTm/H6jdaQ4Lim/H4bWgJ0Mt/AQkQkK5jZLuYX//cIvt8WHxDejTusK/CL4PPgZGBXfLBpDbQByoE6MxsGHLmZ1/uhmfUOAtRK/HrT+gSHLgF62sbdyP6KX2NT55x7q7HXCaZ0PY0fUeoIvNLIoWX4kZDlwTqapqx5+StwnpntFqx9ugofFptjIrCT+e0tWpnZqcBuJJiyJ9IUCjySUs65W/F/sF+Ff/NfCFwMxDrI3AhMBz4APgTeD25rymtMx68zuROoxE8LOKeRwy/DnwGrwv8h/2SD+68FHgmmbp3S4HXWAsPxC/G/xbfX/klw5q9JnHNf4qc4XIqfljaTDWfELsEHq8+At/Bn4B4MHvcccDPwRDDNYFZQT5Ntpob445ry+2342I/xofaz4HeqLm0ikg2q8E0JpppZDT7ozMK/X8ZMBfrgPw9uAk4KRjGqgF/gT15V4j9zNlon00Af4FX8dLApwF0u2Hungf8NLpeZ2ftxt/8NP/q0qdGdmMeAw4H/dY2vsR0LFOF/tneBl5N4XuC75gp/BF7HTz/7gmY2CXDOLcM3f7gUP03+18APnXPfNuf5RGIsiTVyIiIiIi2W+Rb+5zvnIrGJtpkV4Zso9HfOzQ27HpGo0wiPiIiISHa5CJimsCOSnJza7VdEREQkl5nZAnwzgU01CBKROJrSJiIiIiIiOUtT2kREREREJGclNaXNzIYCfwHygfudc2Ma3D8Y3+FjD+A059zTcfdtB9yP37jQAUcHu8kn1LlzZ9ezZ8+m/RQiIpJyM2bM+NY51yXsOqJIn1UiIuFL9nNqs4En2GV3HHAEsAiYZmYTnHMfxR32Jb5t7WUJnuKv+E0MXzGzUmB9gmO+07NnT6ZPn765skREJM3MrOGO5xLQZ5WISPiS/ZxKZoRnEDDPOfdZ8MRPAMcB3wWe2IiNmX0vzJjZbkAr59wrwXHVyRQlIiIiIiKSCsms4emO3zwyZlFwWzJ2wu/a+6yZ/cfM/hSMGImIiIiIiKRdMoHHEtyWbGu3VsBB+KluA4EdSLBju5mNMLPpZja9vLw8yacWERERERHZtGSmtC3CNxyI6QEsTvL5FwH/iZsO9w9gX+CB+IOcc/cB9wEMGDBAfbJFJHTr1q1j0aJFrF69OuxS0q6wsJAePXpQUFAQdikiIiIpl0zgmQb0MbNewFfAacAZST7/NKCDmXVxzpUDhwJa5Skikbdo0SLKysro2bMnZokGunODc45ly5axaNEievXqFXY5IiIiKbfZKW3OuTrgYuBfwBzgKefcbDO73syGA5jZQDNbBJwM3Gtms4PH1uOns71mZh/ip8eNT8+PIiKSOqtXr6ZTp045HXYAzIxOnTq1iJEsERFpmZLah8c5NxGY2OC2a+KuT8NPdUv02Ffw+/OIiGSVXA87MS3l5xQRkZYpmaYFIiIiIiIiWUmBR0QkgpYvX85dd93V5McdffTRLF++PA0ViYiIZCcFHhGRCGos8NTX12/ycRMnTqR9+/bpKktERCTrJLWGR0REMuvKK69k/vz59OvXj4KCAkpLS9l6662ZOXMmH330EccffzwLFy5k9erV/PKXv2TEiBEA9OzZk+nTp1NdXc2wYcM48MADeeedd+jevTvPP/88RUVFIf9kIiIimaXAIyKyOf/zPzBzZmqfs18/GDu20bvHjBnDrFmzmDlzJpMnT+aYY45h1qxZ37WOfvDBB+nYsSO1tbUMHDiQH/3oR3Tq1Ol7zzF37lwef/xxxo8fzymnnMIzzzzDWWedldqfQ0REJOIUeEREssCgQYO+t0/O7bffznPPPQfAwoULmTt37kaBp1evXvTr1w+AvffemwULFmSsXhERkahQ4BER2ZxNjMRkSklJyXfXJ0+ezKuvvsqUKVMoLi7m4IMPTriPTps2bb67np+fT21tbUZqFRERiRI1LRARiaCysjKqqqoS3rdixQo6dOhAcXExH3/8Me+++26GqxMREckeGuERSacvv4Yly6BvbygqDLsaySKdOnXigAMOoG/fvhQVFdGtW7fv7hs6dCj33HMPe+yxBzvvvDP77rtviJVKU1RVQefOzXvsCSfA/fdDaWlqaxIRyXUKPCLpsvAb+Pwrf/3DudBvF2hdEG5NklUee+yxhLe3adOGl156KeF9sXU6nTt3ZtasWd/dftlll6W8Pmm61q1h1KimP66yEsaPhzlz4PnnoWfPlJcmIpKzFHhE0uGrJfDZIujSAbbpCh9+CrPmwZ47QX5+2NWJSEjatIE//KF5jz3hBDj1VBg4EJ55BgYPTm1tIiK5Smt4RFLt63KYtxA6tYddekH7Mth1B6iqgTmfg3NhVygiWeioo+C996BTJzjsMLjvvrArEhHJDgo8Iqm0ZBl8+gV0aAu77QB5wT+xzh2g93awbDnM+1KhR0SaZaedYOpUOOIIGDkSLr4Y1q0LuyoRkWhT4BFJlfIK+PhzP6Lzg94bwk5M967QoxssLvfre0REmqFdO/jnP+Hyy2HcOD/yU1MTdlUiItGlwCOSCt8u99PV2pb6jmz5jfzT2qEHdO3omxksWZbZGkUkZ+Tnwx//CI88Aq+/DtdfH3ZFIiLRpcAjsqXW1cGcz6C0CHbvvemmBGawc09oVwafLIDKlZmqUkRy0E9+Aj/9Kdx6K8yeHXY1IiLRpMAjsqUqV8D69bDjdtAqicaHeXnQd0coagMffQar16S/Rsk6y5cv56677mrWY8eOHcuqVatSXJFE1c03Q9u28LOfaXmgiEgiCjwiW2rZCihoBW1Lkn9Mq1Z+nY9b70PP+vXpq0+ykgKPJKtzZxgzBt54Ax59NOxqRESiR/vwiGwJ56BiJXRs66erNUVxIezcCz6aD/MXQp/t01OjZKUrr7yS+fPn069fP4444gi6du3KU089xZo1azjhhBO47rrrqKmp4ZRTTmHRokXU19dz9dVXs2TJEhYvXswhhxxC586def3118P+USQDzjsPHngALrsMfvhD6NBh08evWOGbH4iItAQKPCJbYmUN1NX5PXeao0sH37lt0RLf8KBbp9TWJ6kx70uoTvGISWmxb1XeiDFjxjBr1ixmzpzJpEmTePrpp3nvvfdwzjF8+HDeeOMNysvL2WabbXjxxRcBWLFiBe3atePWW2/l9ddfp3PnzqmtWSIrLw/uvhsGDICrrvLd2xJxzq/3+fWvYfx4v/5HRCTXaUqbyJaoWO4vO7Zt/nP06u7DzqdfQE1tauqSnDJp0iQmTZrEXnvtRf/+/fn444+ZO3cuu+++O6+++ipXXHEFb775Ju10yr5F22svvy/P3XfD9Okb3796NZx9th8FatUKxo7Vmh8RaRk0wiOyJZatgHalyTUraExent+kdMZHMHs+9N8VWm2i05tk3iZGYjLBOcfo0aMZOXLkRvfNmDGDiRMnMnr0aI488kiuueaaECqUqLj+enjqKbjoInj33Q1NIxcvhhNOgPfe88d06+Y3Lp0yBfbfP9yaRUTSTSM8Is21eq0fkemYgrPqbVr70FO7Gj5doNOuQllZGVVVVQAcddRRPPjgg1RXVwPw1VdfsXTpUhYvXkxxcTFnnXUWl112Ge+///5Gj5WWpV07P2Vt+nS47z5/27RpMHCgb1v97LNw9dVwxhlQVgb33BNuvSIimaARHpHmqljhL5u7fqeh9m399LbPv4K2S/3aHmmxOnXqxAEHHEDfvn0ZNmwYZ5xxBvvttx8ApaWlPProo8ybN4/LL7+cvLw8CgoKuPvuuwEYMWIEw4YNY+utt1bTghbotNPg/vth9Gior/dT2LbaCt55B/bYwx9TWgo//rFvdHDbbdBJywdFJIeZi9iZ5AEDBrjpiSYfi0TNrLlQXQv77N70Dm2Ncc5Pa6tYAQN/AEWFqXleabI5c+aw6667hl1GxiT6ec1shnNuQEglRVrUP6s+/tiHm3XrYPBgePpp6NLl+8d88AHsuSfccguMGhVOnSIiWyLZzylNaRNpjvXrobIKOrVLXdgB/1w79PDBZ3l16p5XRFqUXXbxzQt++1t45ZWNww74QLT//n5aW8TOfYqIpJQCj0hzLK/yoScV63caKmrjVxpX1aT+uUWkxTjvPLjxRmjduvFjLrwQ5s4FzXwUkVymwCPSHMtW+O5q7begHXVjzKCsWIEnAqI25TddWsrPCWBm7c3saTP72MzmmNl+ZtbRzF4xs7nB5Wa27cwdJ50EHTuqeYGI5DYFHpGmcs7vv9O+DPLT9E+orMR3gFu/Pj3PL5tVWFjIsmXLcj4MOOdYtmwZhYUtZr3YX4CXnXO7AHsCc4Argdecc32A14LvW4SiIjjnHHjuOViyJOxqRETSQ13aRJpq1WrfknrbrdL3GmUlPlhVr/KbkkrG9ejRg0WLFlFeXh52KWlXWFhIjx49wi4j7cysLTAYOAfAObcWWGtmxwEHB4c9AkwGrsh8heEYMcK3sn7wQd/ZTUQk1yjwiDRVrB11OtbvxJQV+8sqBZ6wFBQU0KtXr7DLkNTaASgHHjKzPYEZwC+Bbs65rwGcc1+bWddEDzazEcAIgO22C3cz2lTaeWc49FC/b8+vf71hs1IRkVyhKW0iTbVsBZQUQWGb9L1Gm9ZQ0ErreERSqxXQH7jbObcXUEMTpq855+5zzg1wzg3okqjtWRa78EJYsAAmTQq7EhGR1FPgEWmKujpYWZ3e0R0IGheUKPCIpNYiYJFzbmrw/dP4ALTEzLYGCC6XhlRfaI47Drp1U/MCEclNCjwiTVGx0q+t6ZTmwAM+8KxaDXX16X8tkRbAOfcNsNDMdg5uOgz4CJgAnB3cdjbwfAjlhap1a9/G+oUXYOHCsKsREUktBR6RpqhYAa3yM7OupqzEX1ZrlEckhS4B/m5mHwD9gN8DY4AjzGwucETwfYtzwQX+fM7994ddiYhIaqlpgUiynPOBp2M7P+Us3drGNS5Ix34/Ii2Qc24mMCDBXYdlupao6dkThg3b0LygpCTsikREUiOpER4zG2pmn5jZPDPbaIGnmQ02s/fNrM7MTkpwf1sz+8rM7kxF0SKhqKqBdXXpX78TU1AAha21jkdEMuY3v4FvvoEbbwy7EhGR1Nls4DGzfGAcMAzYDTjdzHZrcNiX+H0NHmvkaW4A/t38MkUiYFkG2lE3pMYFIpJBBxwA554Lf/4zfPRR2NWIiKRGMiM8g4B5zrnPgk3angCOiz/AObfAOfcBsNG28Ga2N9ANULNLyW6VK6FtiW8XnSllJX6T07XrMveaItKi3XwzlJXBz3/uZ/KKiGS7ZAJPdyC+Z8ui4LbNMrM84Bbg8qaXJhIh6+r8SEuHDK+lid+AVEQkA7p0gTFjYPJkeKyxeRsiIlkkmcCTaHV2sud8fgZMdM5tssmlmY0ws+lmNr28vDzJpxbJoOUr/WWHDE5nAygNVg1rWpuIZND558OgQXDppbB8edjViIhsmWQCzyJg27jvewCLk3z+/YCLzWwB8GfgJ2a2UbvPXN69WnJE5UrIz/dT2jKpVT4UFyrwiEhG5eXB3XdDeTlcdVXY1YiIbJlkAs80oI+Z9TKz1sBp+E3aNss5d6ZzbjvnXE/gMuCvzrmNuryJRJpzfsPR9mWZaUfdUKxxgSbTi0gG9e/v1/HcdRfMmBF2NSIizbfZwOOcqwMuBv4FzAGecs7NNrPrzWw4gJkNNLNFwMnAvWY2O51Fi2RU7RpYsxY6hrQXTlmJX0O0Zm04ry8iLdYNN0DXrnDRRVBfH3Y1IiLNk1S7KefcRGBig9uuibs+DT/VbVPP8TDwcJMrFAlbZWz9TliBJ65xQWGbcGoQkRapXTu49VY480wYPx4uvDDsikREmi6pjUdFWrTKlX4D0LDCRmmxn0qndTwiEoLTT4dDDoHRo2Hp0rCrERFpOgUekU1Zv953aOvQNpz1O+BXD5cUKfCISCjMYNw4qKmB446Db74JuyIRkaZR4BHZlKoaqF+f+XbUDZWV+CltalwgIiHYdVd4/HH44AMYOBDefz/sikREkqfAI7IpsfU77cvCraOsxK8Yrl0Tbh0i0mL96Efw9tt+0PnAA+HJJ8OuSEQkOQo8IptSudKHjYKk+nukz3eNCzStTUTC068fTJvmW1afdprfo2f9+rCrEhHZNAUekcbU1cHKmvC6s8UrKfKnVRV4RCRkXbvCa6/BeefBTTfBiSdCVVXYVYmINE6BR6Qxy4NP8CgEHjPfrU2BR0QioE0b36b69tvhhRdgv/3gs8/CrkpEJDEFHpHGVKyE/DxoWxJ2JV5ZCVSv0vwREYkEM7jkEnj5ZVi82DczeP31sKsSEdmYAo9IYypX+mYFeRH5Z9K2GNY7qFkddiUiIt85/HCYOhW6dYMjjoC77gq7IhGR74vIX3IiEVO7GlavCb8ddbyyYKRJ09pEJGL69IEpU2DoUPj5z+Gii2Dt2rCrEhHxFHhEEom1o47C+p2YwjbQKl+BR0QiqV07eP55uOIKuOceOPJIKC8PuyoREQi5165IRFWuhDatoahN2JVsEGtcUFMbdiUiIgnl58OYMbD77r6L28CBcOONUFi48bGFhX4KXJsIvc2KSG5S4BFpyDmorIIuHXzIiJJW+bCuLuwqREQ26cwzYaed4Pjj4cc/bvy4/feHZ5/1639ERNJFgUekoaoaqK+P1nS2GDN1aRORrDBwIHzyCSxYkPj+99+HCy+EAQP8VLj+/TNanoi0IAo8Ig1VRHD9Tkxenh+BEhHJAqWl0Ldv4vv69vVT3447Dg48EB5+GE45JaPliUgLoaYFIg1VroSyYiiI4PkAM9+aWkQkB+y1F0yb5kd3Tj0Vrr5ag9giknoKPCLx6uphZXU0R3fABx6N8IhIDunWDV57zTc5uPFGOPFEqKoKuyoRySUKPCLxVgSfsu0jGnjyFHhEJPe0aQPjx8Ptt8MLL/hmBp99FnZVIpIrFHhE4i2v8qMo7UrDriQxTWkTkRxlBpdcAi+/DF99BYMGweTJYVclIrlAgUck3ooqaFvimwNEkaa0iUiOO/xwmDoVunb1+/TcfXfYFYlItovoX3UiIairh6pV0L4s7EoaF5vSptAjIjmsTx+YMgWOOgp+9jO46CJYty7sqkQkWynwiMTE1u+0i3DgiW2EqsAjIjmuXTu/P88VV8A99/jRnm+/DbsqEclGCjwiMSuqfaBoWxJ2JY2LTbVT4BGRFiA/H8aMgUcfhXff9ZuZzp0bdlUikm0UeERilldBWYn/hI2q2AiPGheISAty5pnw5puwbBn8/vdhVyMi2UaBRwSC9Ts10V6/A5rSJiIt1sCBcOyx8OKLUF8fdjUikk0UeETAbzYK0Q88eRrhEZGWa/hwKC/3XdxERJLVKuwCslZtLSxYAKtW+a/a2g3X16zxu6gVF0NRkb+MfbVr57eVjvK0qZYotv9OlNfvgEZ4RKRFGzoUWrWCCRP85qQiIslQ4GmO+fPhsMPgiy+a9/j8fNhqK+je3X9ts42/3G472H57/7XNNgpFmbSiCsqKo/87j43wuPXh1iEiEoJ27eDgg33gGTMm7GpEJFso8DTV3LlwyCGwejU89BB07Pj9EZyiIj+6s3Zt4tGfigq/hXTs65NP4P/+D1as+P7r5OdDjx4bAlDPnv6rVy9/2aMHFBSE8AvIQfXB/js9uoVdyeZZMAtVU9pEpIUaPhx+8Qv/cdynT9jViEg2UOBpik8+8WFn3TofUvbYI3XPXV0NCxf6UaMvv/SXseuTJ/twtD7urH5eng89u+0Ggwb51ZwDB/rpctI0K6r9FLGor98BTWkTkRbv2GN94PnnP2HUqLCrEZFsoMCTrI8+gkMP9X9ovv469O2b2ucvLYVdd/Vfiaxb5wPRggUbvj7/HP77X7jxxg1haNttffAZNMhPcB4wwI86SeNiG462LQ23jmSoaYGItHA9e/rzjRMmKPCISHIUeJIxa5Zfs5OX58OZQAmNAAAgAElEQVROY6EknQoKYIcd/FdDNTXw/vswbdqGr2ef3fC4vfeGAw7Y8NW1a2Zrj7rl1X7/nVYRX78DGuEREcFPa/vDH/y+PJ06hV2NiESd2lJvzgcf+GlsrVr5qWVhhJ3NKSmBgw7yp7oefxzmzfN9OydMgF/9yq8HuuMOOPFEP+Vtjz3g5pth0aKwKw9ffZbsvxOjwCMiwvDh/u37pZfCrkREsoECz6bMn+/DTmEh/PvfsPPOYVeUvM6d/UTnm2+Gt96ClSvh7bf996WlcOWVvivcoYf65gsrV4ZdcThW1vjw0C5LAs93U9rUpU1EWq6994att/bn9URENkeBZ1NefNF3VXv5ZejdO+xqtkybNn5Nz69/De+849vb/O53fl3QT3/qR35OO82vAl27NuxqM2d5sH6nXRas34ENXdo0wiMiLVhenj+n9/LLfus7EZFNUeDZlIoKf7nLLuHWkQ69e/vA8+mn8O67cP758Nprfp7AVlvByJF+VCvXRxJi++9kw/odiNuHR4FHRFq24cOhqsp/VImIbEpSgcfMhprZJ2Y2z8yuTHD/YDN738zqzOykuNv7mdkUM5ttZh+Y2ampLD7tKiv9LmdR34xyS5jBPvv4NT6LF/tRraOPhr//3e/utt12cPnl8OGHYVeaevXr/ZS2bJnOBhvW8KhLm0izmNkCM/vQzGaa2fTgto5m9oqZzQ0uO4Rdp2zeoYf67e80rU1ENmezgcfM8oFxwDBgN+B0M9utwWFfAucAjzW4fRXwE+fcD4ChwFgza7+lRWdMRYXfWLSlKCjwYefRR2HJEt8AYa+9YOxY3+jgpz/1t+eKlVm0/06MRnhEUuEQ51w/59yA4Psrgdecc32A14LvJeKKiuDII33g0VuiiGxKMiM8g4B5zrnPnHNrgSeA4+IPcM4tcM59AKxvcPunzrm5wfXFwFKgS0oqz4SWFnjilZRsWNPz9dd+lOfRR2GnneDWW/2+QNluRZat3wGN8Iikx3HAI8H1R4DjQ6xFmmD4cL8U9b//DbsSEYmyZAJPd2Bh3PeLgtuaxMwGAa2B+U19bGhacuCJ17kz/PGPfj+iAw6ASy/1Iz6TJoVd2ZZZXgWlxb7leLZQW2qRLeWASWY2w8xGBLd1c859DRBcJtyszMxGmNl0M5teXl6eoXJlU445xr8talqbiGxKMoHHEtzWpL+2zGxr4G/Auc65jVbBR/ZDRIHn+3baCSZOhBdegLo6OOooOP54WLAg7Mqabn0Wrt8BtaUW2XIHOOf646dp/9zMBif7QOfcfc65Ac65AV26ZM9khVzWtSvst58Cj4hsWjKBZxGwbdz3PYDFyb6AmbUFXgSucs69m+iYyH6IKPAkdswxfrRnzBjf2W3PPeHJJ8Ouqmli++9k0/od0AiPyBYKplfjnFsKPIeftr0kODEXO0G3NLwKpamGD4cZM7SXtog0LpnAMw3oY2a9zKw1cBqQ1LmU4PjngL865/63+WWGYP16BZ5NadMGrrjCd2/bbTe/3ueCC6CmJuzKkpNt++/E5GkfHpHmMrMSMyuLXQeOBGbhP9PODg47G3g+nAqlOYYP95cvvBBuHSISXZsNPM65OuBi4F/AHOAp59xsM7vezIYDmNlAM1sEnAzca2azg4efAgwGzglagM40s35p+UlSrarKhx4Fnk3r2RPeeANGj4YHHoABA+CDD8KuavNWVEFpERRk0fqdeGpaINIc3YC3zOy/wHvAi865l4ExwBFmNhc4IvhessQuu/it5TStTUQak9Rfe865icDEBrddE3d9Gn6qW8PHPQo8uoU1hiO26agCz+YVFMDvfw+HHQZnnQWDBvlObhddtGEKVpQ4B1WroGsW/rc1818a4RFpMufcZ8CeCW5fBhyW+YokFcz8KM+dd0J5OURpZryIRENSG4+2SJWV/rKD9p9L2mGH+d6ghx4KP/85/OhHfqQsalavgfp6KCsOu5LmyTON8IiIxDn/fD8p40rtoCQiCSjwNEYjPM3TtaufSH3LLX5+wSGHwNKIrf+tXuUvS7M08GiER0Tke3bd1e+Y8OCD8PbbYVcjIlGjwNMYBZ7my8uDUaPg+efho4/83j2ffx52VRtUBYGnpCjcOprLDDbu7i4i0qJdfTVsu62fTV1XF3Y1IhIlCjyNUeDZcscc49tWV1TA/vvDzJlhV+RVr/JhJy9L//fPy9OUNhGRBkpK4PbbffPQO+4IuxoRiZIs/YsvA2KBR2t4tsx++8Fbb/nGBkOGwOTJYVfkA0+2TmcDTWkTEWnEccf5c23XXANffRV2NSISFQo8jamogOJiKCwMu5Lst+uuflJ1jx5w1FHwzDPh1bJmLayry+7Ak6fAIyKSiJkf5amr8zOrRURAgadx2nQ0tbbdFt580+/Tc/LJcP/94dSR7Q0LwH+ia0qbiEhCO+wAv/0tPPUUTJoUdjUiEgUKPI1R4Em9jh3hlVf8KM/IkfDqq5mvIRcCj0Z4REQ26fLLoU8fv0PC6tVhVyMiYVPgaYwCT3oUF/vTbrvuCqecAvPnZ/b1q2uhsA20ys/s66aSmd9wQkREEmrTBsaNg3nz4E9/CrsaEQmbAk9jKivVsCBdysp8y2rwK0wzuTlptjcsADUtEBFJwhFHwKmnwk03Zf7cmohEiwJPYzTCk1477uhHeubMgbPPzsyIRV0drF6jwCMi0kLceiu0bg2jR4ddiYiESYGnMQo86Xf44XDLLfDcc3DDDel/vVxYvwPah0dEJEnbbAMXXgjPPguLF4ddjYiERYEnkdpav8pRgSf9fvlLP8Jz7bXwj3+k97VigacsywOPRnhERJI2YgTU18ODD4ZdiYiERYEnkdimowo86WcG99wDAwfCj38Ms2al77Wqa6F1gf/KZurSJiKStN69/Xqe++7zwUdEWh4FnkQUeDKrsNBPayst9U0MKivT8zq50LAAtA+PiEgTXXghLFwIL70UdiUiEgYFnkQUeDKve3c/yfqLL+Cqq1L//PXroaY2NwJPnoFTW2oRkWQdeyxsvbWfUCAiLY8CTyIKPOHYbz9/Gu6ee1I/ta2m1l/mQuDRCI+ISJMUFMD558PEif68moi0LAo8iSjwhOfaa6FtWxg1KrXrVHKlQxuA5WkNj4hIE51/vj9fNH582JWISKYp8CQSW0OijUczr3Nn+N3v4JVX/Km4VKleBfn5UNg6dc8ZljyN8IiINNV228Exx8D998O6dWFXIyKZpMCTSEUFtGrlF9FL5v3sZ7DTTnDppan7VKpeBaVF/vRetlNbahGRZrnwQliyBJ5/PuxKRCSTFHgSiW06mgt/HGej1q39hqSffAJ3373lz+cc1ORIhzbYEHgUekREmuSoo2D77dW8QKSlUeBJJBZ4JDzHHOM3Trj2Wli2bMuea9VqPwWsrCQlpYUuLwjiCjwiIk2Sn+83In3tNfj007CrEZFMUeBJRIEnfGZw662wYgVcd92WPVcuNSyADSOPCjwiIk3205/6Wev33Rd2JSKSKQo8iSjwREPfvjByJNx1F8yZ0/znqV7lR0WKC1NXW5g0wiMi0mxbbQUnnAAPPQSrV4ddjYhkggJPIgo80XHddb55xKWXNv85qldBSY40LADflhrUqU1EpJkuvNB/1D/9dNiViEgmKPAkosATHV26wDXXwEsv+a+mci7o0JYj09lAU9pERLbQIYf4ZqBqXiDSMijwNLRuHVRVKfBEycUXQ+/efpSnrq5pj129FurqcyvwxKa0aYRHRKRZzPyM6bffhlmzwq5GRNJNgaeh5cv9pTYdjY7WrWHMGL+O57HHmvbYXGtYABrhERFJgbPO8m+nzz4bdiUikm4KPA1VVPhLjfBEy4knQv/+vk312rXJPy4WeEpyKPCoaYGIyBbr2hX23RcmTAi7EhFJNwWehhR4oskMbrwRPv8cHnww+cdVr/Ld2fJz6H/12AjP+vXh1iEikuWGD4cZM2DRorArEZF0yqG/AlNEgSe6hg6F/feHG26A2trkHpNrDQtgQ5c2jfCIiGyR4cP95QsvhFuHiKSXAk9DCjzRZQY33QSLFyfXWmftOv+Va4FHTQtERFJi111hxx01rU0k1ynwNKTAE20HHwyHHw5/+ANUV2/62FxsWABqWiAikiJmfpTntdc2/5EiItlLgaehigr/DtiuXdiVSGNuvBHKy+Evf9n0cVU5GnjUtEBEJGWGD/e9cCZNCrsSEUkXBZ6GKiqgfXvIzw+7EmnMPvvAscfCn/4ElZWNH1dVA0VtoKBV5mrLBNOUNhGRVDngAL8Thaa1ieQuBZ6GKis1nS0b3HADrFgBt9zS+DHVNVBWkrmaMuW7KW3q0iYisqUKCuDoo33jgvr6sKsRkXRIKvCY2VAz+8TM5pnZlQnuH2xm75tZnZmd1OC+s81sbvB1dqoKT5uKCm06mg323BNOOQXGjoWlSze+f+06WLMOynJsOhtAnrq0iYik0vDhsGwZTJkSdiUikg6bDTxmlg+MA4YBuwGnm9luDQ77EjgHeKzBYzsCvwP2AQYBvzOzaKeJigqN8GSL667z7alvvnnj+1bW+MtcHuHRlDYRkZQ46ig/0qNpbSK5KZkRnkHAPOfcZ865tcATwHHxBzjnFjjnPgAazrE5CnjFOVfhnKsEXgGGpqDu9FHgyR677AI/+QmMGwdfffX9+6qDwJNrDQtATQtERFKsXTvfBFSBRyQ3JRN4ugML475fFNyWjC15bDgUeLLL737nJ103XMtTVQMlRbnZfEIjPCIiKTd8OHzyif8SkdySTOCxBLcl+5dWUo81sxFmNt3MppeXlyf51Gmwfr2aFmSbnj3hxBPh4Yf99DbwIx9Vq3Jz/Q5oHx4RkTQ49lh/+c9/hluHiKReMoFnEbBt3Pc9gMVJPn9Sj3XO3eecG+CcG9ClS5cknzoNVq70oUeBJ7uMHOmD6tNP++/XrIV1dbm5fgc2TGlbry5tIiKpsv32vh+OprWJ5J5kAs80oI+Z9TKz1sBpQLJvB/8CjjSzDkGzgiOD26KposJfKvBkl0MOgT594N57/fdVsfU7ORp4NMIjIpIWw4fD22/Dt9+GXYmIpNJmA49zrg64GB9U5gBPOedmm9n1ZjYcwMwGmtki4GTgXjObHTy2ArgBH5qmAdcHt0WTAk92MoMRI/yn1OzZfjqbGZQWhV1ZeijwiIikxfDhfvB84sSwKxGRVEpqHx7n3ETn3E7OuR2dczcFt13jnJsQXJ/mnOvhnCtxznVyzv0g7rEPOud6B18PpefHSJHKSn+pwJN9zjkHWrf2ozxVNT7s5OXovrpm/ktNC0REUqp/f9hmG01rE8k1OfoXYTPFRni08Wj26dwZTjoJ/va3IPDk6HS2GDON8Ig0g5nlm9l/zOyF4PteZjY12Bz7yWDqtrRQeXm+ecHLL8Pq1WFXIyKposATT1PastvIkVDWDurX527Dgpg8BR6RZvolfnp2zM3Abc65PkAlcF4oVUlkDB8ONTUweXLYlYhIqijwxNMIT3Y76CA49HB/PVdbUsdoSptIk5lZD+AY4P7gewMOBYIWjzwCHB9OdRIVhx4KxcWa1iaSSxR44lVUQEkJtGkTdiXSHGZwzHC/H8+8T8OuJr00pU2kOcYCvwZiPd07AcuD5jywmc2xI7NnnKRVYSEMGwZPPbXhPKiIZDcFnngVFZrOlu169oL5c+G++8KuJL3yTPvwiDSBmf0QWOqcmxF/c4JDGz2TEJk94yTtrrkGli+H3/wm7EpEJBUUeOIp8GQ352D1OqAeHn0UqqvDrih9LE8jPCJNcwAw3MwWAE/gp7KNBdqbWavgmKZsrC05bI894Be/8OfOpk4NuxoR2VIKPPEUeLJbTa0f9dhzd6iqgieeCLui9FHTApEmcc6NDrZP6InfQPv/nHNnAq8DJwWHnQ08H1KJEjHXXgtbbw0XXQT19WFXIyJbQoEnngJPdquq8Zf994S+ff2ePLlKTQtEUuUKYJSZzcOv6Xkg5HokItq2hdtug//8B+6+O+xqRGRLKPDEq6xU4MlmVasgPx+KCn2L6unTYcaMzT8uG6lpgUizOecmO+d+GFz/zDk3KNgc+2Tn3Jqw65PoOPlkOOII+O1v4Ztvwq5GRJpLgSfGOY3wZLuqGt+O2gzOOguKinJ3lCdPIzwiIulmBnfe6TchveyysKsRkeZS4ImprYU1a7QHT7Zav96v4YltONq+PZx2Gjz2GKxYEW5t6WAGTl3aRETSbaed4Ior4O9/h9dfD7saEWkOBZ6YWLN9jfBkp+pVfpQuFngAfv5zv132gw+GV1e6qGmBiEjGjB4NvXrBz34Ga9eGXY2INJUCT4wCT3arWuUvy4o33Lb33nDQQfCXv0BdXeLHZSvL05Q2EZEMKSryU9s+/hhuvTXsakSkqRR4YhR4sltVDRS0gjatv3/7qFHwxRfwj3+EU1e6qGmBiEhGHX00nHACXH89LFgQdjUi0hQKPDEKPNmtqsZPZ7MGG6cfeyzssEPunZLTlDYRkYwbO9ZPaXtAzctFsooCT4wCT/aqr4dVq78/nS0mPx/+539gyhR4993M15Yu2odHRCTjttsO9toL3ngj7EpEpCkUeGIUeLLXd+t3ShLff+650K6d30EuV2iER0QkFEOGwNSpvlW1iGQHBZ6Yykpo3RqKE4wSSLRV1fjLxgJPaSlccAE884xfz5MLNMIjIhKKwYP9LhbvvRd2JSKSLAWemNimow3XgEj0Va3yzQpaFzR+zCWX+Ms77shMTemmfXhEREJx0EH+LVjT2kSyhwJPTEWFNh3NVlU1idfvxNtuOzjpJBg/HqqqMlNXOuWpLbWISBg6dIDdd4d//zvsSkQkWQo8MbERHskudXWweg2UNjKdLd6oUbByJTz0UPrrSje1pRYRCc2QIfDOO7BuXdiViEgyFHhiFHiyU3Wtvywt2vyxgwbB/vv7vqL19emtK91iUy8VekREMm7wYFi1CmbMCLsSEUmGAk+MAk92qokFniSbTYwaBZ9/DhMmpK+mTMgLAo+mtYmIZNxBB/lLreMRyQ4KPDEKPNmpZhW0yt90w4J4xx8PPXtm/0akGuEREQlNt26wyy4KPCLZQoEH/LbJ1dUKPNmoutaP7iTbXS8/H37xC3jrLZg2Lb21pVNshEed2kREQjF4MLz5ZvbPkBZpCRR4wO/BAwo82cY5P6WtJIn1O/HOOw/KyuDmm9NTVyZY8E9XU9pEREIxZIjvg/PBB2FXIiKbo8ADCjzZqnYNrF+f/PqdmLZt4Ve/8huRTpmSntrSTVPaRERCNXiwv1R7apHoU+ABv34HFHiyTc0qf9nUER6Ayy+HrbaCSy/NztCgpgUiIqHq0QN22EHreESygQIPbAg82ng0u8RaUjcn8JSWwg03+BGep59ObV2ZoBEeEZHQDR7sA896LacUiTQFHtAIT7aqqYXiQshr5v/G557rt8u+4gpYsya1taVbngKPiEjYBg+GZctgzpywKxGRTVHgAQWebFW9CkqauH4nXn4+/PnPfl+eO+9MXV2ZEBvh0WlFEZHQDBniLzWtTSTaFHjABx4zaNcu7EokWXV1sGYtlDZjOlu8I4+EoUPhxhv9abpsoSltIiKh69ULundX4wKRqFPgAR94OnRo/tQoybwtWb/T0J/+5HuL3nDDlj9XpqhpgYhI6Mz8KM8bb+j8k0iU6S988IFH09myS00QeJrakjqRvn393jzjxsHcuVv+fJkQ24dHn7AiIqEaPBi+/hrmzQu7EhFpjAIPKPBko5pV0CofWhek5vmuvx4KC30Dg2ygER4RkUjQOh6R6FPgAb/xqAJPdqmu9aM7sbUsW2qrrXzYee657PjU0hoeEZFI2Hln6NJF63hEoiypwGNmQ83sEzObZ2ZXJri/jZk9Gdw/1cx6BrcXmNkjZvahmc0xs9GpLT9FNMKTXZzzU9pSsX4n3qhRfvXppZdGv/uZAo+ISCSYbdiPR0SiabOBx8zygXHAMGA34HQz263BYecBlc653sBtwM3B7ScDbZxzuwN7AyNjYShSYk0LJDvUrvGBJBXrd+IVF8NNN8H06TB+fGqfO9Xy1JZaRCQqhgyBL77wXyISPcmM8AwC5jnnPnPOrQWeAI5rcMxxwCPB9aeBw8zMAAeUmFkroAhYC6xMSeWp4hwsX67Ak01qVvnLVI/wAPz4x3D44X60J8orUDXCIyISGYMH+0uN8ohEUzKBpzuwMO77RcFtCY9xztUBK4BO+PBTA3wNfAn82TlXsYU1p1ZtrT9LXlYWdiWSrFS2pG4oLw8eeghat/bhp64u9a+RCnnq0iYiEhW77w7t2yvwiERVMoEn0arwhn9lNXbMIKAe2AboBVxqZjts9AJmI8xsuplNLy8vT6KkFKqu9pclJZl9XWm+mlVQXJi+fZN69IC77oJ334UxY9LzGlvK1KVNRCQq8vLgoIPUuEAkqpL5i3ERsG3c9z2AxY0dE0xfawdUAGcALzvn1jnnlgJvAwMavoBz7j7n3ADn3IAuXbo0/afYEjU1/rK0NLOvK81XXQslKV6/09Dpp8Npp8F118GMGel9rebQlDYRkUgZMsRv5fb112FXIiINJRN4pgF9zKyXmbUGTgMmNDhmAnB2cP0k4P+ccw4/je1Q80qAfYGPU1N6isRGeBR4skNdHaxZC6VpmM7W0Lhx0K0bnHWWn/oYJdqHR0QkUrSORyS6Nht4gjU5FwP/AuYATznnZpvZ9WY2PDjsAaCTmc0DRgGx1tXjgFJgFj44PeSc+yDFP8OWUeDJLt+t30nzCA/4VuUPPQQffwxXbtSNPVwa4RERiZS99vJ/SijwiERPq2QOcs5NBCY2uO2auOur8S2oGz6uOtHtkaLAk11qgsCTiREegCOOgF/8Am6/HX74Q/99FHwXeNSWWkQkClq1ggMOUOARiaI0rfrOIgo82aVmlf9UaV2QudccMwZ22QXOOcfv2RQVeXma0iYiEiGDB8OsWfDtt2FXIiLxFHgUeLJLda0f3bFEjQHTpKgIHn0Uli6Fiy/O3OtujpmmtImIRMiQIf7yrbfCrUNEvk+BR22ps4dzfkpbOvbf2Zy994arroLHH4fXXsv86yeSp8AjIhIlAwZAYaGmtYlEjQKP2lJnj9o1fpPY0gw0LEjkiiugVy+/pmfdunBqiGemKW0iIhHSpg3su6/24xGJGgWe2AhPcUh/REvyalb5yzBGeMCfths7Fj76CO68M5wa4mmER0QkcoYMgZkzYcWKsCsRkRgFnupqP50tT7+KyPuuJXVIgQfg2GNh6FC49lpYsiS8OkAjPCIiETR4sJ+M8PbbYVciIjH6K7+6WtPZskXNKiguDDecmsFf/uI3Ig17bx4ztaUWEYmYffeFggKt4xGJEgUeBZ7sUV0b3vqdeDvtBKNGwcMPw5Qp4dWRl6cpbSIiEVNcDAMHKvCIRIkCjwJPdqirgzVrw53OFu+qq2CbbeCSS6C+PpwaNKVNRCSShgyBadM29EUSkXAp8MTW8Ei0fbd+JwIjPOBD8p/+BDNmwIMPhlOD9uERaRIzKzSz98zsv2Y228yuC27vZWZTzWyumT1pZq3DrlWy2+DB/jzdu++GXYmIgAKPP/2iEZ7oqwkCT2lERngATj8dDjoIRo+GiorMv36eRnhEmmgNcKhzbk+gHzDUzPYFbgZuc871ASqB80KsUXLA/vv7WcdqTy0SDQo8mtKWHapqoKAVtC4Iu5INzOCOO6CyEq65JpzX1wiPSNKcF+xFQEHw5YBDgaeD2x8Bjg+hPMkhbdtC//5axyMSFQo8CjzR5xxUroT2Zf6P/CjZc0+46CK4+26/8UIm5alLm0hTmVm+mc0ElgKvAPOB5c65uuCQRUD3sOqT3DF4sJ/StmZN2JWIiAKPAk/01dTC2nXQoV3YlSR2/fXQqRNccIGftJ0palog0mTOuXrnXD+gBzAI2DXRYYkea2YjzGy6mU0vLy9PZ5mSA4YM8WHnvffCrkREFHgUeKKvcqW/7Ng23Doa07Gjn9o2fbrfoydTTG2pRZrLObccmAzsC7Q3s1bBXT2AxY085j7n3ADn3IAuXbpkplDJWgce6C81rU0kfC078NTVwerVCjxRV7nSbzjaJsKNk045BYYPh6uvhnnzMvOaalog0iRm1sXM2gfXi4DDgTnA68BJwWFnA8+HU6Hkko4dYffd1bhAJApaduCJNchXW+roqq+H5VXRnc4WYwZ33eW3177ggsyMvKhpgUhTbQ28bmYfANOAV5xzLwBXAKPMbB7QCXggxBolhwwZAu+8A+vWhV2JSMumwAMa4YmyFdX+j/qoTmeL1707/PnPMHkyjB+f/tfLU+ARaQrn3AfOub2cc3s45/o6564Pbv/MOTfIOdfbOXeyc07LzCUlBg/2f2r85z9hVyLSsrXswFMddCdV4ImuypV+JKNdlvw3Ov98OOQQuPxyWLQova+lpgUiIpE2eLC/1LQ2kXAp8IACT5RVrPBhJz8/7EqSY+ZHd9at8+2q0zkCoyltIiKR1q0b7LyzGheIhE2BBxR4omrNWli1GjpGfP1OQzvuCDfeCC+8AE8+mb7XyTNYr314RESibMgQePNNvyRVRMKhwAMKPFEVa0fdIQvW7zT0y1/CoEFwySWQrv06LPjnq1EeEZHIGjwYVqyADz8MuxKRlkuBB9SlLaoqVkLrAigpCruSpsvPhwce8J9yF14Iy5al/jXM/KUCj4hIZMXW8Wham0h4WnbgUZe26HLOj/B0aLvhD/ts07cvXHstPPssbLUVHH00/PWvPgSlQl7we1HjAhGRyNp2W+jZE956K+xKRFqulh14NKUtuqpX+Y1hs3E6W7zf/AZmzIBRo+Cjj+Dss/0q1hNOgCee2BC6m0MjPCIiWWHgQP9RICLhUOABBZ4oqghGQbI98AD07w833wyffw5Tpvjube+9B6efDl26wIym5JQAACAASURBVCmnwDPPQG1t0543T4FHRCQbDBgAn30GFRVhVyLSMinwtGoFrVuHXYk0VLkSSov9Gp5cYQb77gu33QYLF/qNGc4911+edJIPP2ecAc8/D6tXJ/d8oE5tIiIRN2CAv9Qoj0g4FHhKS7N3jUiuqquHlTW5MbrTmLw8v5J13Dj46it49VU480yYNAmOPx623homTNj8c4BGeEREIq5/f385fXq4dYi0VAo8paUwax58uzzsaiRmeZX/Iz6XA0+8Vq3gsMPg3nvh66/h5Zehd28ffG65pfFAY2paICKSDdq3hz59FHhEwqLAs3V3WLZ8w5oRCV/lCj960a4Frq0qKICjjvLT3H70I7jsMhg5Etat2/hYNS0QEckaAwYo8IiEpWUHnpoa2HFHf33N2nBrkQ0qV0L7sg1Ttlqi4mJ48knf5W38eBg2DCorv3+MmhaIiGSNAQPgyy9h6dKwKxFpeVrwX5T4EZ5te/rrCjzRULvGf3VsIdPZNiUvD266CR5+2O9Yt99+MH/+hvs1pU1EJGuocYFIeBR4tunuryvwRENlrB11u3DriJKzz/ZNDcrLYZ994J13/O3fTWlTlzYRkajbay//tq1pbSKZp8DTdSt/va4e6uvDrUf8dLY2raGoTdiVRMvgwTB1KpSUwKWX+tvyNMIjIpItyspgl10UeETCoMDTsfOGtSJrEiwMl8xZvx4qq/x0NrUK31jv3nDAAbBsmf9eTQtERLKKGheIhKNlBx4Mior9AnnQtLawrazxo2yazta40lIf1EH78IiIZJkBA2DxYv8lIpmTVOAxs6Fm9omZzTOzKxPc38bMngzun2pmPePu28PMppjZbDP70MwKU1f+FnAOOnXx12ML5BV4wvVtpR+1UMOCxsUHHjUtEBHJKmpcIBKOzQYeM8sHxgHDgN2A081stwaHnQdUOud6A7cBNwePbQU8ClzonPsBcDAQjXlja9fCttv56+0VeELnnN/8tWM7yM8Pu5roigUe5zSlTUQky/Tr5wfnNa1NJLOSGeEZBP/f3t3HyV3W9/5/fZK9S3YgQAgRwq2AlCBSNYB3hVZBaK1BbC14ePykRw+Qtng8x2qLbfV4rPbGPqrWFq147IGHpxUUbY2W4q/FG9TDragIghICSABJAiGwudvd7HX+uGbYyWbDzmxm5zvzndfz8djHd+Y739n97LWTfOc91/W9LtaklNamlEaBq4FzphxzDnBV9fa1wGsiIoDXAnemlH4IkFJ6IqXUGTMDjIzAEUfBxE5YOAT9fQaeIo1sze1/4H5FV9LZKpUccLZtc9ICSeoyCxfCCScYeKR2ayTwLAMerru/rrpv2mNSSuPAZmAx8AIgRcTXIuKOiPiDvS+5RWqBZ2w0f1I+OGDgKdKG6qKaiw08z6lSyduREXt4JKkL1SYu8L9uqX0aCTzTTZc19Z/pno7pA14FXFDdnhsRr9ntB0RcHBG3R8TtGzZsaKCkFhgZgSOOhFqH0+CAs7QVaeNTefKI/r6iK+ls9YFnnuvwSFK3WbEC1q+Hhx+e+VhJrdFI4FkHHFZ3/1Bg6vwizx5TvW5nEfBkdf+3UkobU0pbgeuAl0z9ASmlK1JKK1JKK5YsWdL8bzEbz4zA8w6ZfNM42G8PT1G2bINt2+HA/YuupPPt0sNT/efrkDZJ6hq1iQsc1ia1TyOB5zbg2Ig4KiIGgPOB1VOOWQ1cWL39m8DXU0oJ+BrwoohYWA1CpwM/bk3pe2nL1rwdrPYoDA64+GhRNlaHs3n9zsym7eEx8EhSt3jRi6Cvz8AjtdOM44dSSuMRcSk5vMwH/iGldHdEfAC4PaW0GvgM8NmIWEPu2Tm/+txNEfERcmhKwHUppX+do9+lOTvG8m+/oDpL9uBAdf8oLFxQWFk9aeNTsM/w5N9AezbdNTz28EhS1xgaghNPNPBI7dTQBRMppevIw9Hq972v7vZ24E17eO7/IU9N3VkmEoyPw77VYVS1N9vbDTxttX1HnqHt+YcWXUl3qAWeLVuctECSutSKFXDttbuuMCBp7jS08Gg5zYdHHoZ99sl3a4Fn1IkL2mrjU3nrcLbG1PfwQD5TGngkqausWAGbNsEDDxRdidQbejfw9A/Agw9MvoEc7M/b7U5c0FYbN8HwgsmhhXpuw8N5Wws88wImnKVNkrqJExdI7dWbgWdiAgaH4GcPTL6BnDfPxUfbbXQMNo/Yu9OM3Xp45tnDI0ld5oUvhIEBA4/ULr0ZeLbtyG8UH3sU5s+f3D/k4qNt9URtOJvTUTdswYI8jK2+h8fAI0ldZWAATjrJwCO1S28Gnq3b8nbD47vuHzDwtNWGTTA0mIe0qTHz5uVeyfpreJylTZK6zooV8L3vOSpZaofeDDxbtudPxTdv2nX/0ECerlpzb3wcnnomD2dziprmVCr28EhSl1uxAp5+GtasKboSqfx6M/Bs3Q5PP5V7dOoN9OeFR8ddfHTOPbE5v1F3OFvz6gOPPTyS1JWcuEBqnx4NPNtgw/rJC8BrhuoWH9Xc2vhUDpj7DhddSfeZGnjs4ZGkrrN8eV6E1MAjzb3eCzwp5R6exx7ZPfAMGnjaYudOeHKzw9lma7chbQ4Al6Ru09cHL36xgUdqh94LPNt35NDz8EMGnqJsejpfpelwttlxSJsklcKKFXDHHflzQElzp/cCz5btefvA/bsHnoHq4qMGnrm1YRP0zYdFlZmP1e52CTyuwyNJ3WrFCtiyBe66q+hKpHLrvcBTm5L6vp/uHnjmzcuhx5na5s6O0Rx4lhyQ21vNmzqkzR4eSepKZ50F/f1w1VVFVyKVW++949y6PYea9Y/n9UymGuy3h2curXs890gc9ryiK+leTlogSaWwdCm88Y1w5ZWwbVvR1Ujl1ZuBZ8EQbN26ew8P5Ot4DDxzY2wcHt0ABx0ACwaLrqZ7uQ6PJJXGqlWwaRN84QtFVyKVV28FnpTykLaB+fm+gae9Hnk8T1Zw+MFFV9LdKhUYHYWxMSctkKQud/rpcNxx8Pd/X3QlUnn1VuAZHYOdE1CbCXlPgWfnBIyPt7W00hvfCY+sh8X7wfCCoqvpbrXX7ZYt1SFtTkstSd0qIvfy3HQT/PCHRVcjlVNvBZ4t1QGyE9VJCfYUeMCJC1rtsQ059BzutTt7rXbt2chInvjBHh5J6mpveUtehPRTnyq6EqmceivwbK1OSb2jup028Dg1dctNTOTJCvbbB/Z1Kuq9Vnvdjow4aYEklcABB8B558FnPwvPPFN0NVL59F7g6ZufhwLBDD08Bp6W+fnGPJzQa3daoz7wOGmB1LCIOCwivhER90TE3RHxjur+AyLi3yPivurWVZHVdqtW5f/WP/e5oiuRyqfHAs82WDgEW6uBZ7ppqV18tLVSgod/DvsM5x4e7b2pPTwOaZMaNQ78fkrpeOBlwO9FxHLgMuCGlNKxwA3V+1JbnXoqnHRSnrzAz7Gk1uqxwLMdFi6YnNJ3uh6eZxcfNfC0xPonYfto7t2JmPl4zcwhbdKspJQeSyndUb39DHAPsAw4B6gt/XgV8IZiKlQvq01e8P3vw223FV2NVC69E3jGxvI6MAuHnjvwQHVqaict2Gspwc8ey7OyLV5UdDXlMd2QNkOP1JSIOBJ4MXALsDSl9BjkUAQcVFxl6mUXXJD/i3eKaqm1eifwbKlOVDA8Qw8PuBZPqzzxVO5VO+x59u600tQeHjDwSE2IiArwReC/pZSebuJ5F0fE7RFx+4YNG+auQPWsffbJoefqq/NipJJao3cCT22GtoZ6ePrzMCzfRM5erXdnaBAOOqDoaspllx6e6j9hX6tSQyKinxx2/jGl9KXq7scj4uDq4wcD66d7bkrpipTSipTSiiVLlrSnYPWcSy6BbdvyjG2SWqO3As+8ebn3pvZGcWho+mMHB/JUyjt3trfGMnniKXhmq707c6F+HZ5a2zpxgTSjiAjgM8A9KaWP1D20GriwevtC4Mvtrk2qefGL4ZRTnLxAaqXeCTxjY3kygoj8RrFS2fMb8drU1Nsd1jYr23fATx7Mwweft7joasqnvx8GBx3SJjXvlcD/B7w6In5Q/fo14C+AMyPiPuDM6n2pMKtWwT33wLe/XXQlUjn0FV1A24yNQ3/1192yZfopqWvq1+KpLJz72spkYgJ+vDa/AV9+9OSQK7VWpTI5aQEYeKQGpJS+A+ypy/k17axFei7nnQf//b/DJz8Jp51WdDVS9+udd6NjYzBQDTy1Hp49eTbwOFNb09aug2e2wHFH5uulNDdqgcchbZJUOgsXwoUXwhe/COunvaJMUjN6J/CM1vXwzBh4XHx0VjZsgkfWw7KDYIkTFcyp3Xp4JoqtR5LUUpdckj+rvfLKoiuRul9vBJ6UqkPaqkFmpsAT4eKjzdq2PV+3s88wPP/QoqspP3t4JKnUli+H00+HT30qjxaXNHu9EXh2TuTQ02gPD7gWTzN2TsCP788j45c/3+t22mF4OF+LFk5LLUlltWoVrF0L//EfRVcidbfeeGc6Vr0Wp5nAM2Tgadj9P4ORbfALR+V1dzT3pg5ps4dHkkrn3HNhyZI8RbWk2euRwDOet80EnoGBPGmBn5w/t8efgMc25vV2Fu9XdDW9Y+qQNl+nklQ6g4Pw1rfC6tXwyCNFVyN1rx4LPNVreGaalhpgqD8Pmh138dE9enoEfvoQLKrAUcuKrqa3GHgkqSdcfHFeB/0znym6Eql79VjgafIaHnBY255s2wF3rclTfS8/es+LuGpu7DakzStaJamMnv98OOss+PSnYXy86Gqk7tQbgWe0eg3PQB+MjuYvA8/sjY3DXfflXoUTj80z2qm9KpXcU1ljD48kldaqVbBuHVx3XdGVSN2pocATEWdHxE8iYk1EXDbN44MRcU318Vsi4sgpjx8eESMR8a7WlN2ksfH8Sfi8eZNvEg08szMxAXffn3t4TjgGFi4ouqLeVKnkkDO6I9838EhSaf36r8Mhhzh5gTRbMwaeiJgPXA78KrAceHNELJ9y2NuATSmlY4CPAn855fGPAv+29+XOUm0Nnog8DAgamLSgtvjo2NzW1k1SytfsbH4GjjsS9tun6Ip6V+31u21b3jpLmySVVl8fXHQRXH89PPBA0dVI3aeRHp5TgDUppbUppVHgauCcKcecA1xVvX0t8JqIfFFHRLwBWAvc3ZqSZ2FsfNfrd2DmwBPhWjxTPfRYnpXtiENg6eKiq+lttdfv1q15aw+PJJXaf/kv+a3Jpz9ddCVS92kk8CwDHq67v666b9pjUkrjwGZgcUQMA38I/M+9L3UvzCbwAAz2G3hqfr4RHno0B50jDi66Gk0NPPbwSFKpHXoovP71eba2Ud+aSE1pJPBMN/3W1HdXezrmfwIfTSmNPOcPiLg4Im6PiNs3bNjQQElNGhubDDy1a3hmmpYa7OGp2bApD2Xbbx94wRHOyNYJaoGn9nq2h0eSSm/VKli/Hv7lX4quROoujQSedcBhdfcPBR7d0zER0QcsAp4ETgU+HBEPAv8N+KOIuHTqD0gpXZFSWpFSWrFkyZKmf4kZzbqHpxp4evXN5OgY/Pj+/DW8AE44Ok/8oOLtFnicllqSyu61r4Ujj3TyAqlZjbx7vQ04NiKOiogB4Hxg9ZRjVgMXVm//JvD1lP1SSunIlNKRwMeAP0sp/V2Lam/MxATsnJhcdLTZwDORJqe17hUp5Wt1brsbNj4FRx4CL/6FfNWkOkPt9Vt7PTukTZJKb948uOQS+MY34M47i65G6h4zBp7qNTmXAl8D7gE+n1K6OyI+EBErq4d9hnzNzhrgncBuU1cXZrS6StfALHp4arOQbXiy9XV1qh2jeUHRex+ABYPw0uV5kgJ7djrLs4Hnmbzt1V5ISeoxb30r7L8/vOY1cOONRVcjdYeGPrJPKV0HXDdl3/vqbm8H3jTD93j/LOrbe2PV3pnZDGmrLIR9huGxjbBsabmvXZmYyBMTrF2Xr746+jBYdlC5f+duZg+PJPWkgw6Cm2+GlStz6Ln8crj44qKrkjpb+T+2H6v28EwNPI1MWgBwyBLYuh02P+e8C91rfCc8/HO49S6472c54K04AQ4tecDrdrXX75Yt+e9kD48k9YwXvABuuQXOOCMPcbv00snPdyXtrocCT901PIODk/dnsmR/mD8fHpuD2eOKtGM09+bcfGfeDg3AC4+BF70gD2VTZ1uwYHIh3XkGHknqNYsWwVe/Cu96V+7lOesseOKJoquSOlP5r0Kf2sOzZUvjvTuQw87SA/KwtmPqZnvrVlu2wbqfw+NP5jfJS/aHQ58H+zbRJipeRB7WNjKSbzukTZJ6zvz58Fd/BSeeCBddBCefDKtXwwtfWHRlUmfp8nfvDaj18fbNz9uRkcau36l38BJ4dEOeuezQpa2trx1Sgs3PwMOPw5Ob8wQEBy/Jv4u9Od2rPvA4LbUk9ay3vCUPczv3XHj1q2HdOhgYKLoqqXP0xpC2/r7J61FmE3ienbxgQ3cNHUoJ1j8Jd9wDP/wpPLMlTzH9shfBsYcbdrpdLfDMm9ddr0tJUsu97GVwxRWwYQN861tFVyN1lvL38IyOw0Dd9TqzCTyQJy/4yYN58oLadNWdZGIih7ux8bxu0JZt8Mj6fK3OgkE49ghYuhjmlz/j9gyHtEmS6pxxRr7Ec/VqOPPMoquROkf5A8/YlOtuZht4luwPax7OvTxFBp6dO+HpLXmI2uYR2L4j/447pxnStG8FjjkcFi9yxrUyeraHx0kLJEk57Lz2tTnwfPzjnvqlmt4IPJUFk/dHRuDww5v/PkVMXpBS7qEZ2ZrDzeaRfLv25rayMIea/r7qV//k7cF+WDA09zWqOJUKbNxoD48k6VkrV8KXvwx33gknnVR0NVJn6IHAMwb9dT0ys+3hgdlPXpBSXu9mdCz30KSUF/esPVbbbh+Fbdth247Jbe3xiHwd0aFLcw/TvpXJiRjUmyoVePBB1+GRJD3rda/Lp4XVqw08Uk25A08taNT3xjQ7LXW9+skLlh00fV/x9h2w7vEcXkbHJr8afUMaka+5WTAEByyChUP5a5/hfHG6VFM/pG3CWdokSbB0aZ7AYPVqeO97i65G6gzlDjxTFx2FvevhgcnJC54egUV1PUcTE3mSgAcfBVIOLAP9OawM9MNAX97O74NgMiw9uwUGB/KXg27ViF2mpbaHR5KUrVwJ73kPPPIILFtWdDVS8Xok8FR/zZT2PvDUJi94dMNk4Hl6BH76UJ4ZbfF+eaKAISfA1xzbZZY2e3gkSVkt8Hz1q3DJJUVXIxWv3GOkRquLjtYCz7ZtOfTsTeCpTV6wYVO+xua+h+D79+ZwdcLR8MJjDDtqj0plcmFdJy2QJFUdfzwcfXQe1iap7IGn1sNTW4dnZCRv9ybwQJ68ICW47a7c07PsIDj5hXDg/nv3faVm1F7HEzsd0iZJelZE7uW54YbJtz5SL+uNwFPr4WlV4KkshAP2zdNdv+T4PITNGdPUbrXJN3YaeCRJu1q5EnbsgH//96IrkYrXG4GnFka2bMnbvQ08ACe+AF6yPM+eJhWh9joeH3dImyRpF698Jey/v8PaJCh94BnLYac2nXOth2e201JLnaQ+8NjDI0mq098Pv/ZreeKCnTuLrkYqVskDz/iua/C0akib1AmeDTxjztImSdrNypWwcSPcfHPRlUjF6oHAM2UNHjDwqBxqr+OxJha2lST1jLPOym+DHNamXtcDgcceHpVU7XU8OmrgkSTtZtEi+OVfNvBIBh6pW9UHHictkCRNY+VKuPde+OlPi65EKk55A09KBh6VW+11vGNHfr3byyNJmuL1r8/br3yl2DqkIpU38IxX1yYZqAs8W7bk1bgWLCiuLqlVarMN7thebB2SpI51xBFw0kkOa1NvK2/geXbR0SmTFixcODlNtdTN+vpgaAi2VwOPM7VJkqaxciV85zvwxBNFVyIVo7zv/MfG8nbqkDaHs6lMKhXYvi3fdkibJGkaK1fmz8S+/OWiK5GKUeLAU+vhMfCoxCoV2Fbr4THwSJJ295KXwPHHw6WXwjXXFF2N1H4GHqmbVSqwdWu+bQ+PJGka8+bBN7+Zg8/558Of/ImjoNVbeiDwTLmGx8CjMqlUYOuWfNvAIz2niPiHiFgfEXfV7TsgIv49Iu6rbvcvskZprhx0EHz96/C2t8GHPgRvfCM880zRVUntUd7AMzqeP9KYX/crGnhUNvU9PA5pk2ZyJXD2lH2XATeklI4Fbqjel0ppYAA+/Wn4+Mfhq1+Fl78c1q4tuipp7pU38IyN7TolNeRpqQ08KpPh4fy6Bnt4pBmklG4Enpyy+xzgqurtq4A3tLUoqc0i4O1vh+uvh0cfhZNPhm98o+iqpLlV4sAzZdFRyD08tbVLpDKoVCYX1LWHR5qNpSmlxwCq24MKrkdqizPOgFtvhaVL4cwz4ROfKLoiae6UPPD077rPIW0qm0qlrofHK1CluRQRF0fE7RFx+4YNG4ouR9prxxwDN98MZ58Nv/d7sGoVjI4WXZXUeiUPPNP08Bh4VCaVCoxUrzq1h0eajccj4mCA6nb9ng5MKV2RUlqRUlqxZMmSthUozaV9983r81x2GXzqU7m3xzyvsumdwDM+nlekN/CoTCqVyWl2vIZHmo3VwIXV2xcCLs2onjN/Pvz5n8M//mMe5nbyyXDnnUVXJbVOOQPPzp15gvn6wFMb9mPgUZlUKjnMg4FHmkFEfA64CTguItZFxNuAvwDOjIj7gDOr96We9J/+E9x4Y5736RWvgC99qeiKpNbom/mQLrSnNXjAwKNyqVRgfCzfdkib9JxSSm/ew0OvaWshUgc7+WS4/XY491z4jd+Aj30M3vGOoquS9k5DPTwRcXZE/CQi1kTEbmsURMRgRFxTffyWiDiyuv/MiPheRPyoun11a8vfg9Fa4LGHRyVXqeSP4sAeHklSSxx8MHzzm7ByJbz73XDvvUVXJO2dGQNPRMwHLgd+FVgOvDkilk857G3AppTSMcBHgb+s7t8IvD6ldCJ5bPRnW1X4c6q9Aaxfh6fWw+O01CqTSgV21oa0OUubJKk1hobyIqXDw3kGNz9TUzdrpIfnFGBNSmltSmkUuJq8UFu9+oXbrgVeExGRUvp+SunR6v67gaGIGGxF4c9pbJoeHoe0qYwqlcnXu0PaJEktdNBBeTKDr38dPve5oquRZq+RwLMMeLju/rrqvmmPSSmNA5uBxVOO+Q3g+ymlHbMrtQkGHvWK+mt4/PhNktRiF12Ur+t55zth8+aiq5Fmp5HAE9Psm/rO6jmPiYgTyMPcLpn2B7R6MbexcYjI8yzWGHhURvWztNnDI0lqsfnz4ZOfhPXr4b3vLboaaXYaCTzrgMPq7h8KPLqnYyKiD1gEPFm9fyjwz8BbUkr3T/cDWr6YW20NnqjLYQYelZGTFkiS5thLXwq/+7tw+eVwxx1FVyM1r5HAcxtwbEQcFREDwPnkhdrq1S/c9pvA11NKKSL2A/4VeE9K6butKnpGY2O7DmcDA4/KySFtkqQ2+OAH4cAD4Xd+Jy91KHWTGQNP9ZqcS4GvAfcAn08p3R0RH4iIldXDPgMsjog1wDuB2tTVlwLHAO+NiB9Uvw5q+W8x1ej47oGnNi21s7SpTIaH64a0eQaSJM2N/faDv/5ruPVW+F//q+hqpOY0tPBoSuk64Lop+95Xd3s78KZpnvdB4IN7WWPzxsZhwZRgMzICfX0wMND2cqQ5MzQ0GXTs4ZEkzaELLshh57LL8sKkrbgKQWqHhhYe7Tpj0/Tw/OxnsO++u17XI3W7iDysbWKnkxZIkuZUBHziE/DMM/CHf1h0NVLjyhd4JiZg585dA89DD8E118Bv/VZxdUlzpVLJr3t7eCRJc2z5cvj934f//b/hO98puhqpMeULPNOtwfOhD+WPJf7oj4qpSZpLlUoO+QYeSVIbvPe9cPjheQKD2kShUicrceDpz9sHHsgfQ1x0ERx22J6fJ3WrWuBxSJskqQ2Gh+Fv/gbuugv+9m+LrkaaWYkDT7WH54MfzKtmvec9xdUkzaXa1NT28EiS2uScc+B1r4P/8T9g3bqiq5GeWwkDT7Vvtb8P7r8frroKVq2CZcuKrUuaK5VKDvoGHklSm0TAxz+eV0Z45zuLrkZ6buULPKN1PTx/+qd5aJtTiajMKpUc9F2HR5LURs9/PvzxH8MXvgBf+1rR1Uh7Vr7AUxvS9sBa+Oxn4Xd/Fw4+uNiapLlUqcDYqD08kqS2e/e74dhj4dJLYfv2oquRplfOwNPfl6/dGRyEP/iDoiuS5lalAqM7nLRAktR2g4Nw+eWwZg18+MNFVyNNr5yBZ2In/NM/5Y8bli4tuiJpblUqMGoPjySpGGeeCeedB3/2Z/nyaanTlDDwjMHDP4MFC3I/q1R2tcCzc2fRlUiSetRHPgIDA/mzZj9/U6cpX+AZ2QL3/RTe/nZYsqToaqS5V6nkaXLGx4uuRJLUow45BD7wAbj+evjnfy66GmlX5Qw8W7fAu95VdCVSewwPVwOPPTySpOJceimcdFIeYGMvjzpJuQLPnXfC4BAcdxwsXlx0NVJ71KaldkibJKlAfX1w8cWwdi088EDR1UiTyhV4Dj8C5s2DV76i6Eqk9qlUYNx1eCRJxTv99Lz91reKrUOqV67A0z+Qt/vtW2wdUjvVruEx8EiSCnb88XmQzY03Fl2JNKmv6AJaanAATjga9hkuuhKpfWqBxwHTkqSCzZsHp51m4FFnKVcPT998OHD/HHykXlG7hse8I0nqAKedlq/jWbeu6EqkrFyBR+pFz/bwFF2IJEmT1/HYy6NOYeCRul1t0oIouhBJkuBFL4J993XiAnUOA4/U7YaH85A2E48kqQPMnw+vepU9POocBh6p282fDxEwz8AjSeoMp58O994L69cXXYlk4JHKYV5AzHOmNklSRzjttLy1l0edPefFFwAADFtJREFUwMAjlcG86j9lA48kqQO89KWwcKGBR53BwCOVQV91Sa0JA48kqXj9/fCKVzhxgTqDgUcqg/nz8zZNFFuHJElVp58OP/oRPPlk0ZWo1xl4pDKwh0eS1GFOOy2PtP7ud4uuRL3OwCOVQS3weA2PJKlDnHIKDA46rE3FM/BIZTBgD48kqbMMDcGppzpxgYpn4JHKoH8gb+3hkSR1kNNOgzvugGeeKboS9TIDj1QGgwYeSVLnOf102LkT/u//LboS9TIDj1QGA4N5Oz5ebB2SJNV5+cvzZaZex6MiGXikMhiqBp6t24qtQ5KkOsPDeRFSr+NRkQw8UhnUAs82A48kqbOcfjrcequnKBXHwCOVwdBQ3m7dWmwdkiRNcdppMDYGN99cdCXqVQYeqQwWLMjbbduLrUOSpCle9SqIcFibitNQ4ImIsyPiJxGxJiIum+bxwYi4pvr4LRFxZN1j76nu/0lEnNW60iU9qxZ4tjteQJqNmc5zkmZv0SL4xV904gIVZ8bAExHzgcuBXwWWA2+OiOVTDnsbsCmldAzwUeAvq89dDpwPnACcDXyi+v0ktdLChXm7fUexdUhdqMHznKS9cPrpcNNNMDpadCXqRX0NHHMKsCaltBYgIq4GzgF+XHfMOcD7q7evBf4uIqK6/+qU0g7ggYhYU/1+N7WmfEkADC+E0R0wMgKbny66GnWKyjDM9zOmBjRynpO0F047DT72Mfj2t+HUU4uuRp1iaChPWz7XGvkRy4CH6+6vA6a+VJ89JqU0HhGbgcXV/TdPee6yWVcraXr77gubNsDiQ+AHPy26GnWKA4fhhOOLrqIbNHKek7QXfumX8nU8Z5xRdCXqJDfcAK9+9dz/nEYCT0yzb+py7ns6ppHnEhEXAxcDHH744Q2UJGkXRx4Bd/zQdXi0q+PbcBYpB89V0hw78ED44hfh/vuLrkSd5Jhj2vNzGgk864DD6u4fCjy6h2PWRUQfsAh4ssHnklK6ArgCYMWKFbudZCQ14I0ri65A6laeq6Q2OPfcoitQr2pklrbbgGMj4qiIGCBPQrB6yjGrgQurt38T+HpKKVX3n1+dxe0o4Fjg1taULklSSzRynpMkdakZe3iq1+RcCnwNmA/8Q0rp7oj4AHB7Smk18Bngs9VJCZ4knyyoHvd58oWf48DvpZR2ztHvIklS0/Z0niu4LElSizQ0L0JK6Trguin73ld3ezvwpj0890PAh/aiRkmS5tR05zlJUjk0tPCoJEmSJHUjA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0oqUUtE17CIiNgAPAYuAzXUPHQhsnKMfO/Vnteo5Mx2zp8en2z9130z356q9ZtNWjT6vqPbqtNdWo8+zvZp73nMdszdtNd2+Tm6vRp+zCNgvpbSk6ap6QPVc9RTdf56a6bhuPE9N97Na9RzP6809z/NUc89r1b/F6faXtb0aO0+llDryC7hiyv3b2/WzWvWcmY7Z0+PT7Z+mPWa6PyftNZu26vT26rTXlu3V/vbam7bqtvZq9Dmz/Vv00lcZzlMzHdeN56m5bC/P693RXp322mpFezXTVrbX7l+dPKTtKx3+sxp5zkzH7Onx6fZP3TfT/bky259je7X+ebZXc897rmP2pq2m29fJ7dXoc9r5f3C3KsN5aqbjuvH/kdn+LP/fbf3zbK/mnteqf4vT7e+19tpFxw1p25OIuD2ltKLoOrqF7dU426o5tldzbK/e4d+6ObZXc2yvxtlWzemF9urkHp6prii6gC5jezXOtmqO7dUc26t3+Lduju3VHNurcbZVc0rfXl3TwyNJkiRJzeqmHh5JkiRJaoqBR5IkSVJpGXgkSZIklZaBR5IkSVJplSLwRMRwRHwvIn696Fo6XUQcHxF/HxHXRsTvFF1Pp4uIN0TEpyPiyxHx2qLr6XQR8fyI+ExEXFt0LZ2q+v/VVdXX1QVF16P28VzVGM9TzfE81RzPUzMr43mq0MATEf8QEesj4q4p+8+OiJ9ExJqIuKyBb/WHwOfnpsrO0Yr2Sindk1JaBfwWUOo511vUXv+SUroI+G3gvDkst3Ataq+1KaW3zW2lnafJtnsjcG31dbWy7cWqaZ6rGud5qjmep5rjeWr2ev08VXQPz5XA2fU7ImI+cDnwq8By4M0RsTwiToyIr075OigizgB+DDze7uILcCV72V7V56wEvgPc0N7y2+5KWtBeVX9SfV6ZXUnr2qvXXEmDbQccCjxcPWxnG2vU7F2J56pGXYnnqWZcieepZlyJ56nZupIePk/1FfnDU0o3RsSRU3afAqxJKa0FiIirgXNSSn8O7DYMICJ+BRgm/6G2RcR1KaWJOS28IK1or+r3WQ2sjoh/Bf5p7iouVoteXwH8BfBvKaU75rbiYrXq9dWLmmk7YB35ZPIDiv/QSQ3wXNU4z1PN8TzVHM9Ts9fr56lCA88eLGMyVUJu9FP3dHBK6Y8BIuK3gY1lPIHMoKn2iohfJndVDgLXzWllnamp9gLeDpwBLIqIY1JKfz+XxXWgZl9fi4EPAS+OiPdUTzi9ak9t93Hg7yLidcBXiihMLeG5qnGep5rjeao5nqdmr2fOU50YeGKafWmmJ6WUrmx9KV2hqfZKKX0T+OZcFdMFmm2vj5P/4feqZtvrCWDV3JXTVaZtu5TSFuA/t7sYtZznqsZ5nmqO56nmeJ6avZ45T3ViN9U64LC6+4cCjxZUSzewvZpjezXH9po9267c/Ps2zrZqju3VHNtr9nqm7Tox8NwGHBsRR0XEAHA+sLrgmjqZ7dUc26s5ttfs2Xbl5t+3cbZVc2yv5thes9czbVf0tNSfA24CjouIdRHxtpTSOHAp8DXgHuDzKaW7i6yzU9hezbG9mmN7zZ5tV27+fRtnWzXH9mqO7TV7vd52kdKMQ44lSZIkqSt14pA2SZIkSWoJA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwqOtFxMgcfM8HI+LAIn72bOqQJHUuz1NSsQw8kiRJkkrLwKNSiojXR8QtEfH9iPiPiFha3f/+iLgqIv7/6qdSb4yID0fEjyLi+ojor/s2746IW6tfx1Sff1RE3BQRt0XEn9b9vEpE3BARd1S/1znT1PQ7EfHhuvu/HRF/W739LxHxvYi4OyIunua5R0bEXXX33xUR76/ePrpa+/ci4tsR8QvV/W+KiLsi4ocRcePetqkkqXU8T3meUvsYeFRW3wFellJ6MXA18Ad1jx0NvA44B/g/wDdSSicC26r7a55OKZ0C/B3wseq+vwE+mVI6Gfh53bHbgXNTSi8BfgX464iIKTVdC7yx7v55wDXV229NKb0UWAH814hY3MTvegXw9urz3wV8orr/fcBZKaWTgJVNfD9J0tzzPOV5Sm3SV3QB0hw5FLgmIg4GBoAH6h77t5TSWET8CJgPXF/d/yPgyLrjPle3/Wj19iuB36je/izwl9XbAfxZRJwGTADLgKXUnWxSShsiYm1EvAy4DzgO+G714f8aEedWbx8GHAs8MdMvGREV4BXAF+rOW4PV7XeBKyPi88CXZvpekqS28jzleUptYuBRWf0t8JGU0uqI+GXg/XWP7QBIKU1ExFhKKVX3T7Drv4nUwO2aC4AlwEurJ6kHgaFpjrsG+C3gXuCfU0qpWt8ZwMtTSlsj4pvTPHecXXtka4/PA55KKf3i1B+UUloVEaeSPw38QUT8YkppxpOTJKktPE95nlKbOKRNZbUIeKR6+8JZfo/z6rY3VW9/Fzi/evuCKT9vffUk8ivAEXv4nl8C3gC8mclhAouATdWTyC8AL5vmeY8DB0XE4ogYBH4dIKX0NPBARLwJILKTqrePTindklJ6H7CR/ImcJKkzeJ7yPKU2sYdHZbAwItbV3f8I+ZOyL0TEI8DNwFGz+L6DEXEL+YOBN1f3vQP4p4h4B/DFumP/EfhKRNwO/ID8ydhuUkqbIuLHwPKU0q3V3dcDqyLiTuAn1XqnPm8sIj4A3EIe9lD//S8APhkRfwL0k8eC/xD4q4g4ljyM4YbqPklS+3me8jylAsVkL6kkSZIklYtD2iRJkiSVloFHkiRJUmkZeCRJkiSVloFHkiRJUmkZeCRJkiSVloFHkiRJUmkZeCRJkiSV1v8DOIW5DGn2jTgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7d56f0c978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Single-voxel performance evaluation.\n",
    "voxidx = 3\n",
    "\n",
    "print(\"Results: subject\", subid, \", voxel id\", voxidx)\n",
    "\n",
    "myfig = plt.figure(figsize=(14,7))\n",
    "ax_corr = myfig.add_subplot(1, 2, 1)\n",
    "plt.title(\"Correlation coefficient\")\n",
    "plt.xlabel(\"Lambda values\")\n",
    "ax_corr.set_xscale('log')\n",
    "ax_corr.plot(todo_lambda, dict_corr_tr[voxidx], label=\"train\", color=\"red\")\n",
    "ax_corr.plot(todo_lambda, dict_corr_te[voxidx], label=\"test\", color=\"pink\")\n",
    "ax_corr.legend(loc=1,ncol=1)\n",
    "\n",
    "ax_spar = myfig.add_subplot(1, 2, 2)\n",
    "plt.title(\"Sparsity via l0-norm\")\n",
    "plt.xlabel(\"Lambda values\")\n",
    "ax_spar.set_xscale('log')\n",
    "ax_spar.plot(todo_lambda, dict_l0norm[voxidx], color=\"blue\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we look over all voxels, and then look at average performance over all voxels in ROIs of interest. First is to collect the \"best\" performance (over all $\\lambda$ values) achieved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_voxels = len(dict_corr_tr)\n",
    "\n",
    "best_corr_tr = np.zeros(num_voxels, dtype=np.float32)\n",
    "best_corr_te = np.zeros(num_voxels, dtype=np.float32)\n",
    "\n",
    "for v in range(num_voxels):\n",
    "    \n",
    "    # Best absolute correlation value.\n",
    "    best_corr_tr[v] = np.max(np.abs(dict_corr_tr[v]))\n",
    "    best_corr_te[v] = np.max(np.abs(dict_corr_te[v]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we iterate over ROIs, collecting the relevant indices each time. Fortunately, our hierarchical data set will come in very handy here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_roi_corr_tr = {}\n",
    "dict_roi_corr_te = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = tables.open_file(\"data/vim-2/response.h5\", mode=\"r\")\n",
    "tocheck = f.get_node((\"/sub\"+str(subid)), \"idx\")\n",
    "for idxnode in tocheck._f_iter_nodes():\n",
    "    idx = idxnode.read()\n",
    "    roi_name = idxnode._v_name\n",
    "    dict_roi_corr_tr[roi_name] = np.mean(best_corr_tr[idx])\n",
    "    dict_roi_corr_te[roi_name] = np.mean(best_corr_te[idx])\n",
    "\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAAGrCAYAAADtiicXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xu4ZXV95/n3RyqAXCyIXOQmpQ0OrSmwnVI7KkK3bcRG28wEATUKTlpCMzqxWxJN1IjOaIOBJ6bbRIKxhygxhZSaQbFRYgcUIzpVgkCpKJKSsqq8oFBy0ULLb/+x1sHN8Zw6+5y9z6Xq9349z35qXX5rre/ea+9T67N/a62dqkKSJEmSWvKIxS5AkiRJkhaaQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSQOSrEhSSZbNcfk/SvJX465rrpL89yRn7GD+pUn+nx3Mvy/J4+enOs2nJCcm+fYO5l+c5M0LWdMwkuyR5CtJHjPGdb4syafG3XaIdX0kyUnjWJek8TMISZpWkg1JftwfDN+d5KokR4xpvf9mHDUupqkONKvqHVX17xerpsmq6vlV9dcASc5Mcv0sl9+nqu6Yy7b7kPVg//75YZJrkhwzqc3hSf4myQ+S3J/ki0leMKlNJTlqLjVoelV1dlX936OsY5j3VJJTk/xjkgeSXDvEas8CPlNV3+mX32FYH0ZV/U1V/ca42w7hfODtY1qXpDEzCEmayQurah/gEOC7wH9d5HqGNlWvzlx7ejRn7+zfP4cBm4D3TcxI8qvA9cCDwJOAA4A/BT6Y5JRFqHVaSXZb7Bp2Yj8E3kUXCobxu8AHhl35Uv5MV9UXgUclWbXYtUj6ZQYhSUOpqp8Aa4AnTkzrT2G5MMmdSb7bn2rzyH7eAUk+nuSevjfgs0kekeQDwGOBj/U9BX8w1faSvCjJTUl+lOSbE6eXJDk0yZX9Om9P8qqBZc5LsibJZUl+BJw5zbRHJHlDv94fJPlQf1A+VR2vTPLVJPcmuSPJ7/bT9wb+O3Bo/zzu62s7L8llA8v/uyTr+9fh2iT/fGDehiTnJrk5ydYklyfZc6Z9keRx/foe0Y//VZLvDcy/LMlr++Frk/z7frsXA7/e13rPwCr373v77k3yhST/bGBdD/XG9N/M//l0bXekqn4MfAh48sDk/wjcB/xOVX2nqn5cVX9L9w36RUkyxGtxcpIb+/fJxiTnDcy7OsmrJ7X/cpL/vR8+pu+l+mGS25KcOtDu0iTvSfKJJPcD/2pH2+qXeUWSb/XvqTdnoOdzpvdc/x546TTPccrPUj/vYb1lmaL3JN3pmnf19bxsurZJXpDuM3dPuh6cYwfmHZHuNK/v9/W/e4b31EOq6u+r6kPA5qnmT6r1scA/A77Qj58FvAz4g34bH+unb0jy+iQ3A/cnWTbw+t6b7tS6/21gvQ/ruepft7OTfCNdb/efT7zfZtl2tyQX9a/vPyV5dX751NprgZNneu6SFp5BSNJQkuwFnAbcMDD5AuAJdAe3R9F96//H/bzXAd8GDgQOBv4IqKp6OXAnfU9TVb1zim09DXg/8PvAfsCzgQ397L/t13socArwjiTPGVj8RXSBbT/gb6aZ9n8Bvwmc0K/nbuDPp3nq3wNeADwKeCXwp0meUlX3A88HNvfPY5+qetiBXpIn9PW+tn8dPkEXAHcfaHYqcBLwOOBY4Mxp6nhIVf0T8CPgX/STjgfuyy9C1rOB6yYt81XgbODzfa37Dcx+CfBWYH/gdnZ8Ks9s2j4kXXB8Sb/MhOcCH66qn09q/iG6sPyEIVZ9P/AKun17MvAfkvxmP++D/TYnangicCRwVV/PNX2bg/p2f5HkSQPrfmn//Pal67madlv9uv+C7qD9EGA53edhwg7fc1V1bFV9cJrnOOVnaYjXBuAxdD1thwFnAJck+V8mN0ryFOC/0fXGPBr4S+DKdF927AZ8HPgWsKJf1+oZ3lNztRK4o6p+BlBVl9B9Zt/Zb+OFA21fQrcf9uvbf5Pus7Cc7j16WZJDdrCtFwBPBY6j+xw+bw5tX0X3d+DJwFPo9vFkX+2Xk7TEGIQkzeTv+m96f0R34PonAP03oq8C/mNV/bCq7gXeAZzeL/dTugPCI6vqp1X12aoa9uDtd4D/VlXXVNXPq2pTVX0t3fVJzwJeX1U/qaqbgL8CXj6w7Oer6u/65X48zbTfBd5YVd+uqm3AecApmeIUm6q6qqq+WZ3rgE/RHWwN4zTgqv55/BS4EHgk8IyBNv+lqjZX1Q+Bj/HwHpMduQ44Ib+4oHxNP/44utD25SHXA/CRqvpifzD5NzPUMJu2AOf275976fbd4L46ANgyxTJbBubvUFVdW1W39Pv2ZrrgeUI/+6PAk5Mc2Y+/rK9/G92B7Yaq+n+r6mdV9SXgw3ThesL/V1Wf69f9kxm2dQrwsaq6vqoepPtCYPD9PvR7bgqjfJYA3lxV2/r371V0B/KTvQr4y6r6QlVt768r2wb8S+BpdOHt96vq/v61mNW1ZrOwH917ZRj/pao2TnzOq+qK/rP086q6HPgGXe3TOb+q7qmqO4F/YMfv5enangr8Wb9f72bq0//u7Z+XpCXGICRpJr/Zf9O7B/Bq4Lr+4PtAYC9gXX8qzT3A1f106ALT7cCn0p1S9oZZbPMIum93JzsUmAhdE77Fw7953zjFcpOnHQl8dKDurwLb6b5tf5gkz09yQ39K0j3Av2WIA/SBer81MdL3fGycVO93BoYfAPYZct3XASfS9f58hu70mxP6x2en6GXZkdnUMGXb/vSriVMELx5oc2H//lkB/BgY7I24i+4Af7JDBubvUJKnJ/mH/pStrXQ9FAcA9O+Tq/hFOD+dX/QSHgk8feI90O/bl9H1oEx42PtmR9ui29cPta+qB4AfDCw+9HtuCqN8lu7uey8nfKuvdbIjgddNej2O6NseAXxropdmnt1N1wM3jMn75xUDp/bdA/waO/6sjvy+Z9J+n1xTb19gytMGJS0ug5CkofTfEn+E7uDtWXQHqT8GnlRV+/WP5f2F8VTVvVX1uqp6PPBC4D8NnMI207fZG+muE5hsM/CrSQYPlB5LdxH+Q6VOVf4U63/+QN37VdWeVTW4HpLsQddLcCFwcH9A/wlg4tqVmZ7HZroDzIn1he6gctO0SwzvOrqeqRP74euBZ9IFoeumWWY2vQizUt3d8iZOETx7ivl3Ar8H/Fn668iAvwd+K/31LgNOpdtHXx9i0x8ErgSOqKrldNesDF5b9LfAS5L8Ol1v3D/00zcC1016D+xTVf9hsOxZbGsLcPhEw/45Pnpg2aHec1OZ4bP0AN0XEhMm33J6//40wAmPZeprdTYCb59U317VXbO1EXjsNL1X435P3Qw8ftK2ptvGQ9P7Xr/30n1Z8+j+s3orD38vzIeH7Xe6z/dk/5zZ9dBKWiAGIUlDSedFdNeGfLXvcXgv3TUzB/VtDkvyvH74BUmO6g/+f0QXoLb3q/susKPfpnkf8Mokz0l3kflhSY6pqo3APwL/Ocme6S7m/h1+8S3/sC4G3j5xylSSA/vnNtnudD1h3wd+luT5wOBtdb8LPDrJ8mm28yHg5P55/ArdtR7b+ucwo/6i6xOnmldV36ALor9Nd6vhH/X1/BbTB6HvAodPukZpwVTVNXQH4Wf1k/6U7jS+9yV5TL9PXwK8ke40rGEOsvel6yX8SX9t2eQbDnyCLoy+Dbh8oKfs48ATkrw8ya/0j6cOXGc1222tAV6Y5Bn96/tWHn4QvsP3XLqL/8+caqMzfJZuAl6a7qL9k/jFqXqD3ppk9yTH050SeMUUbd4LnN33eiXJ3uluDrEv8EW6A/7z++l7Jnlmv9yM76m+tj2BZcAj+uV/Zaq2VfVtfvmUtpn+XgDsTReMvt9v85V0PULz7UPA7/V/o/YDXj9FmxPobqwiaYkxCEmayceS3Ed3APZ24IyqWt/Pez3dKTs3pLsj29/zi1Ofju7H7wM+D/xFVV3bz/vPwJv6U1jOnbzB6m45+0q6A+WtdAf2Ez0rL6E7zWoz3TUgb+kPsGfjz+i+2f9UknvpbgDx9CnquJfuIvcP0Z2y89J+uYn5X6Prcbijfy6HTlr+Nrqg8l/petBeSHeTiAdnKjDJ4XSv3S07aHYd8IO+t2ViPMCN07T/H8B64DtJZjztbJ78Cd0dwPaoqh/Q9S7uCXyF7lSy/wS8vL/GYxjnAG/r9+Mf0+2rh/TX43wE+Dd0PToT0++lC7Wn072XvkN384895rKt/jPxGmA1XWi4l+5GG9v6JtO+5/oQ8WgefiOSQTv6LP0e3ftq4tS+v5u07Hfo3rub6b4wOLt/3z5MVa2lu07o3X372+lv3FFV2/ttHEV3o5Nv013/BsO9p15OF9rfQ9eL+WO64DWdv+Th15K9D3hi/xmb/Pwm6v8KcBHd6/NdupsufG4H2xiX99JdN3gz3efuE8DP6INqkqcC9/d/0yQtMZnd9ZaSpIWQ5LfpTjv8w8WuRbOXZB+6cHJ0dXf521HbZwH/Z1W9ZEftxi3J+4Hbq+ptC7ndmfSnpN4IPKeqprqZxpLV9xpfXFUTPX8fBt5XVZ9Y3MokTcUgJEnSGCR5IfBpul65i+h6fJ4yyzu8LYj+GpzPAe+uqqF/vFQP118L9q/oeoUOprum8Iaqeu2iFiZpKJ4aJ0nSeLyI7hS0zXSns52+FENQ7zt0PVYfXuxCdnKhux7sbrperK/yi99Sk7TE2SMkSZIkqTn2CEmSJElqzjC/aL1kHHDAAbVixYrFLkOSJEm7qB/8oPst5Ec/+tEztNRStW7duruq6sCZ2u1UQWjFihWsXbt2scuQJEnSLurSSy8F4Mwzz1zUOjR3Sb41TDtPjZMkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKas2yxC5iNWzZtZcUbrlrsMiRJkjSkDeefvNglSFOyR0iSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmjPWIJTkyCTrktyUZH2Ss3fQdkOSA8a5fUmSJEkaxrIxr28L8Iyq2pZkH+DWJFdW1ebBRkl2G/N2JUmSJGloc+4RSnJBknMGxs8DXlNV2/pJewyuP8l9Sd6W5AvAr/eTX5PkS0luSXLMXGuRJEmSpNkY5dS41cBpA+OnAlckOSLJzcBG4IKB3qC9gVur6ulVdX0/7a6qegrwHuDcqTaS5Kwka5Os3f7A1hHKlSRJkqTOnINQVd0IHJTk0CTHAXdX1Z1VtbGqjgWOAs5IcnC/yHbgw5NW85H+33XAimm2c0lVraqqVbvttXyu5UqSJEnSQ0a9WcIa4BS6nqHVgzP6nqD1wPH9pJ9U1fZJy0+cRred8V+vJEmSJElTGjUIrQZOpwtDa5IcnuSRAEn2B54J3DbiNiRJkiRprEbqhamq9Un2BTZV1ZYkzwUuSlJAgAur6pZxFCpJkiRJ4zLy6WhVtXJg+Brg2Gna7TNpfMXA8FrgxFFrkSRJkqRhjPUHVSVJkiRpZ2AQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnNG/kHVhbTysOWsPf/kxS5DkiRJ0k7OHiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpObsVLfPvmXTVla84arFLkOSJEmLYIM/o6IxskdIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkpoz9iCU5Mgk65LclGR9krMH5m1IckCSFUluHfe2JUmSJGkYy+ZhnVuAZ1TVtiT7ALcmubKqNs/DtiRJkiRp1kbqEUpyQZJzBsbPA15TVdv6SXvsYBu7JXlv32v0qSSPHKUWSZIkSRrWqKfGrQZOGxg/FbgiyRFJbgY2AhdM0xt0NPDnVfUk4B7gt6baQJKzkqxNsnb7A1tHLFeSJEmSRgxCVXUjcFCSQ5McB9xdVXdW1caqOhY4CjgjycFTLP5PVXVTP7wOWDHNNi6pqlVVtWq3vZaPUq4kSZIkAeO5WcIa4BS6nqHVgzP6nqD1wPFTLLdtYHg783O9kiRJkiT9knEEodXA6XRhaE2Swyeu90myP/BM4LYxbEeSJEmSxmLkXpiqWp9kX2BTVW1J8lzgoiQFBLiwqm4ZdTuSJEmSNC5jOR2tqlYODF8DHDtNuxX94F3Arw1Mv3AcdUiSJEnSMMb+g6qSJEmStNQZhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktScsfyg6kJZedhy1p5/8mKXIUmSJGknZ4+QJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzdqrbZ9+yaSsr3nDVYpchSZKkOdjgz6BoCbFHSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaM+sglOTIJOuS3JRkfZKzB+bdN+Q6Lk1yymy3LUmSJEnjsGwOy2wBnlFV25LsA9ya5Mqq2jzMwknmsk1JkiRJGpsd9ggluSDJOQPj5wGvqapt/aQ9Jq8jyUVJvpTk00kO7Kddm+QdSa4Dfq9v+uwk/5jkDnuHJEmSJC2kmU6NWw2cNjB+KnBFkiOS3AxsBC4Y6A3aG/hSVT0FuA54y8Cy+1XVCVV1UT9+CPAs4AXA+dMVkOSsJGuTrN3+wNahn5gkSZIkTWeHQaiqbgQOSnJokuOAu6vqzqraWFXHAkcBZyQ5uF/k58Dl/fBldEFnwuU83N9V1c+r6ivAwUyjqi6pqlVVtWq3vZbP4qlJkiRJ0tSGuVnCGuAUup6h1YMz+p6g9cDx0yxbA8P3T5q3bWA4Q9QhSZIkSWMxTBBaDZxOF4bWJDk8ySMBkuwPPBO4bWB9E9f7vBS4frzlSpIkSdLoZryDW1WtT7IvsKmqtiR5LnBRkqLrybmwqm7pm98PPCnJOmArD7++SJIkSZKWhKFuZV1VKweGrwGOnabdPv3gmydNP3HS+JnTLCdJkiRJ827WP6gqSZIkSTs7g5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaM9QPqi4VKw9bztrzT17sMiRJkiTt5OwRkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqTqpqsWsY2h6HHF2HnPGuxS5DkiRJu6iTdv8aAFc/eMyU8zf4Uy5LXpJ1VbVqpnb2CEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkpoz5yCU5Mgk65LclGR9krMH5t03zTKXJjmlH96Q5IC5bl+SJEmS5mrZCMtuAZ5RVduS7APcmuTKqto8ptokSZIkaV4M1SOU5IIk5wyMnwe8pqq29ZP2mLyuJBcl+VKSTyc5cJpVv6Zvc0uSY+ZQvyRJkiTN2rCnxq0GThsYPxW4IskRSW4GNgIXDPQG7Q18qaqeAlwHvGWa9d7Vt3kPcO5UDZKclWRtkrXbH9g6ZLmSJEmSNL2hglBV3QgclOTQJMcBd1fVnVW1saqOBY4CzkhycL/Iz4HL++HLgGdNs+qP9P+uA1ZMs+1LqmpVVa3aba/lw5QrSZIkSTs0m5slrAFOoesZWj04o+8JWg8cP82yNc30iVPrtjPa9UqSJEmSNLTZBKHVwOl0YWhNksOTPBIgyf7AM4HbBtZ7Sj/8UuD68ZQrSZIkSaMbuhemqtYn2RfYVFVbkjwXuChJAQEurKpb+ub3A09Ksg7YysOvL5IkSZKkRTWr09GqauXA8DXAsdO026cffPOk6WcODK8YGF4LnDibWiRJkiRprub8g6qSJEmStLMyCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5s/pB1cW28rDlrD3/5MUuQ5IkSbuoSy/9PgAXn+kx567OHiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOakqha7hqHtccjRdcgZ71rsMiRJkrSLOmn3rwFw9YPHALDBn27Z6SRZV1WrZmpnj5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNWdOQSjJkUnWJbkpyfokZ89hHdcmWTWX7UuSJEnSKJbNcbktwDOqaluSfYBbk1xZVZuHWTjJbnPcriRJkiSNbMYeoSQXJDlnYPw84DVVta2ftMfgepK8J8navqforQPTNyT54yTXAy/uJ784yReTfD3J8eN4QpIkSZI0k2FOjVsNnDYwfipwRZIjktwMbAQuGOgNemNVrQKOBU5IcuzAsj+pqmdV1ep+fFlVPQ14LfCWqTae5Kw+WK3d/sDWWTw1SZIkSZrajEGoqm4EDkpyaJLjgLur6s6q2lhVxwJHAWckObhf5NQkXwJuBJ4EPHFgdZdPWv1H+n/XASum2f4lVbWqqlbtttfyoZ+YJEmSJE1n2GuE1gCnAI+h6yF6SFVtTrIeOD7JOuBc4KlVdXeSS4E9B5rfP2m9E6fXbZ9FLZIkSZI0kmHvGrcaOJ0uDK1JcniSRwIk2R94JnAb8Ci6sLO17yF6/vhLliRJkqTRDNULU1Xrk+wLbKqqLUmeC1yUpIAAF1bVLQBJbgTWA3cAn5unuiVJkiRpzoY+Ha2qVg4MX0N3M4Sp2p05zfQVk8ZPHBi+i2muEZIkSZKkcZvTD6pKkiRJ0s7MICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmDP2DqkvBysOWs/b8kxe7DEmSJO2iLr30+wBcfKbHnLs6e4QkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkpqTqlrsGoa2xyFH1yFnvGuxy5AkSdIu6qTdvwbA1Q8es8iV7Bw2LMGftkmyrqpWzdTOHiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJas6cg1CSI5OsS3JTkvVJzh5imWuTrOqH75vrtiVJkiRpFMtGWHYL8Iyq2pZkH+DWJFdW1eYx1SZJkiRJ82KoHqEkFyQ5Z2D8POA1VbWtn7TH4LqSvCfJ2r6n6K07WO/bk3w5yQ1JDp7bU5AkSZKk2Rn21LjVwGkD46cCVyQ5IsnNwEbggoHeoDdW1SrgWOCEJMdOsc69gRuq6jjgM8CrptpwkrP6ULV2+wNbhyxXkiRJkqY3VBCqqhuBg5IcmuQ44O6qurOqNlbVscBRwBkDvTqnJvkScCPwJOCJU6z2QeDj/fA6YMU0276kqlZV1ard9lo+9BOTJEmSpOnM5hqhNcApwGPoeogeUlWbk6wHjk+yDjgXeGpV3Z3kUmDPKdb306qqfnj7LGuRJEmSpDmbzV3jVgOn04WhNUkOT/JIgCT7A88EbgMeBdwPbO17iJ4/3pIlSZIkaTRD98JU1fok+wKbqmpLkucCFyUpIMCFVXULQJIbgfXAHcDn5qFuSZIkSZqzWZ2OVlUrB4avobsZwlTtzpxm+okDw/sMDK+hO/VOkiRJkubdnH9QVZIkSZJ2VgYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNWdWP6i62FYetpy155+82GVIkiRpF3Xppd8H4OIzPebc1dkjJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUnJ3q9tm3bNrKijdctdhlSJIkaQFs8GdTNI/sEZIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5ow1CCV5cpLPJ1mf5OYkp+2g7bVJVo1z+5IkSZI0jGVjXt8DwCuq6htJDgXWJflkVd0z2CjJbmPeriRJkiQNbc49QkkuSHLOwPh5wAur6hsAVbUZ+B5wYD9/Q5I/TnI98OJ+sRcn+WKSryc5fq61SJIkSdJsjHJq3Gpg8NS3U4ErJkaSPA3YHfjmQJufVNWzqmp1P76sqp4GvBZ4y1QbSXJWkrVJ1m5/YOsI5UqSJElSZ85BqKpuBA5KcmiS44C7q+pOgCSHAB8AXllVPx9Y7PJJq/lI/+86YMU027mkqlZV1ard9lo+13IlSZIk6SGjXiO0BjgFeAxdDxFJHgVcBbypqm6Y1P7+SePb+n+3j6EWSZIkSRrKqOFjNfBe4ADghCS7Ax8F3l9VV+xwSUmSJElaJCPdPruq1gP7ApuqagvddULPBs5MclP/ePIY6pQkSZKksRn5dLSqWjkwfBlw2TTtVkwaP3Fg+C6muUZIkiRJksZtrD+oKkmSJEk7A4OQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktScZYtdwGysPGw5a88/ebHLkCRJkrSTs0dIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5O9Xts2/ZtJUVb7hqscuQJEnSErfBn1zRDOwRkiRJktQcg5BNDeQTAAAMoElEQVQkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmjD0IJXlyks8nWZ/k5iSnDcy7Nsmqfvi+cW9bkiRJkoaxbB7W+QDwiqr6RpJDgXVJPllV98zDtiRJkiRp1kbqEUpyQZJzBsbPA15YVd8AqKrNwPeAA6dZ/u1JvpzkhiQHj1KLJEmSJA1r1FPjVgOnDYyfClwxMZLkacDuwDenWHZv4IaqOg74DPCqqTaQ5Kwka5Os3f7A1hHLlSRJkqQRg1BV3QgclOTQJMcBd1fVnQBJDgE+ALyyqn4+xeIPAh/vh9cBK6bZxiVVtaqqVu221/JRypUkSZIkYDzXCK0BTgEeQ9dDRJJHAVcBb6qqG6ZZ7qdVVf3w9jHVIkmSJEkzGkf4WA28FzgAOCHJ7sBHgfdX1RU7XFKSJEmSFsHIt8+uqvXAvsCmqtpCd53Qs4Ezk9zUP5486nYkSZIkaVzGcjpaVa0cGL4MuGyadicODO8zMLyG7hQ7SZIkSZp3Y/9BVUmSJEla6gxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJas5YflB1oaw8bDlrzz95scuQJEmStJOzR0iSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTmpqsWuYWh7HHJ0HXLGuxa7DEmSJO2iTtr9awBc/eAxY1/3Bn8GZkEkWVdVq2ZqZ4+QJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5hiEJEmSJDVn7EEoydVJ7kny8RnaXZtk1bi3L0mSJEkzmY8eoT8BXr6jBkl2m4ftSpIkSdJQ5hyEklyQ5JyB8fOSvK6qPg3cO0X7DUn+OMn1wIv7yS9O8sUkX09y/FxrkSRJkqTZGKVHaDVw2sD4qcAVMyzzk6p6VlWt7seXVdXTgNcCb5lqgSRnJVmbZO32B7aOUK4kSZIkdeYchKrqRuCgJIcmOQ64u6runGGxyyeNf6T/dx2wYprtXFJVq6pq1W57LZ9ruZIkSZL0kGUjLr8GOAV4DF0P0UzunzS+rf93+xhqkSRJkqShjBo+VgPvBQ4AThi9HEmSJEmafyPdNa6q1gP7ApuqagtAks/SXSv0nCTfTvK80cuUJEmSpPEZ+XS0qlo5aXzKu79V1YpJ4ycODN/FNNcISZIkSdK4zcfvCEmSJEnSkmYQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnNG/kHVhbTysOWsPf/kxS5DkiRJu6hLL/0+ABef6THnrs4eIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5qSqFruGoe1xyNF1yBnvWuwyJEmStIs6afevAXD1g8csciXjt6GRn6FJsq6qVs3Uzh4hSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWrOvAShJFcnuSfJxydNvzbJqn74vvnYtiRJkiTNZL56hP4EePk8rVuSJEmSRjJSEEpyQZJzBsbPS/K6qvo0cO8Qy789yZeT3JDk4FFqkSRJkqRhjdojtBo4bWD8VOCKIZfdG7ihqo4DPgO8aqpGSc5KsjbJ2u0PbB2pWEmSJEmCEYNQVd0IHJTk0CTHAXdX1Z1DLv4gMHEN0TpgxTTbuKSqVlXVqt32Wj5KuZIkSZIEwLIxrGMNcArwGLoeomH9tKqqH94+plokSZIkaUbjCB+rgfcCBwAnjGF9kiRJkjSvRr5rXFWtB/YFNlXVFoAkn6W7Vug5Sb6d5HmjbkeSJEmSxmUsp6NV1cpJ48dP0+7EgeF9BobX0J1iJ0mSJEnzbr5+R0iSJEmSliyDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUnGWLXcBsrDxsOWvPP3mxy5AkSdIu6tJLvw/AxWd6zLmrs0dIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJak6qarFrGFqSe4HbFrsOzegA4K7FLkI75D5a+txHOwf309LnPto5uJ+Wvp1pHx1ZVQfO1GjZQlQyRrdV1arFLkI7lmSt+2lpcx8tfe6jnYP7aelzH+0c3E9L3664jzw1TpIkSVJzDEKSJEmSmrOzBaFLFrsADcX9tPS5j5Y+99HOwf209LmPdg7up6Vvl9tHO9XNEiRJkiRpHHa2HiFJkiRJGplBSJIkSVJzlkwQSnJSktuS3J7kDVPM3yPJ5f38LyRZMTDvD/vptyV53kLW3ZK57qMkz02yLskt/b//eqFrb8kon6V+/mOT3Jfk3IWquTUj/r07Nsnnk6zvP1N7LmTtLRnhb96vJPnrfv98NckfLnTtrRhiHz07yZeS/CzJKZPmnZHkG/3jjIWrui1z3UdJnjzwt+7mJKctbOVtGeWz1M9/VJJNSd69MBWPSVUt+gPYDfgm8Hhgd+DLwBMntTkHuLgfPh24vB9+Yt9+D+Bx/Xp2W+zntKs9RtxH/wI4tB/+NWDTYj+fXfUxyn4amP9h4Arg3MV+PrviY8TP0jLgZuC4fvzR/r1bkvvppcDqfngvYAOwYrGf0672GHIfrQCOBd4PnDIw/VeBO/p/9++H91/s57SrPUbcR08Aju6HDwW2APst9nPaFR+j7KeB+X8GfBB492I/n9k8lkqP0NOA26vqjqp6EFgNvGhSmxcBf90PrwGekyT99NVVta2q/gm4vV+fxmvO+6iqbqyqzf309cCeSfZYkKrbM8pniSS/SXdAsH6B6m3RKPvoN4Cbq+rLAFX1g6ravkB1t2aU/VTA3kmWAY8EHgR+tDBlN2XGfVRVG6rqZuDnk5Z9HnBNVf2wqu4GrgFOWoiiGzPnfVRVX6+qb/TDm4HvAQcuTNnNGeWzRJL/FTgY+NRCFDtOSyUIHQZsHBj/dj9tyjZV9TNgK923ocMsq9GNso8G/RZwY1Vtm6c6Wzfn/ZRkb+D1wFsXoM6WjfJZegJQST7Zn6LwBwtQb6tG2U9rgPvpvsG+E7iwqn443wU3aJT//z12WBhjeZ2TPI2up+KbY6pLDzfn/ZTkEcBFwO/PQ13zbtliF9DLFNMm39d7ujbDLKvRjbKPupnJk4AL6L7V1vwYZT+9FfjTqrqv7yDS/BhlHy0DngU8FXgA+HSSdVX16fGWKEbbT08DttOdzrM/8Nkkf19Vd4y3xOaN8v+/xw4LY+TXOckhwAeAM6rql3ojNBaj7KdzgE9U1cad8dhhqfQIfRs4YmD8cGDzdG360w2WAz8cclmNbpR9RJLDgY8Cr6gqv9GZP6Psp6cD70yyAXgt8EdJXj3fBTdo1L9311XVXVX1APAJ4CnzXnGbRtlPLwWurqqfVtX3gM8Bq+a94vaM8v+/xw4LY6TXOcmjgKuAN1XVDWOuTb8wyn76deDV/bHDhcArkpw/3vLmz1IJQv8/cHSSxyXZne6i0ysntbkSmLiryynA/6ju6qwrgdP7u/c8Djga+OIC1d2SOe+jJPvR/SH7w6r63IJV3KY576eqOr6qVlTVCuBdwDuqaue6+8vOYZS/d58Ejk2yV3/gfQLwlQWquzWj7Kc7gX+dzt7AvwS+tkB1t2SYfTSdTwK/kWT/JPvTnanwyXmqs2Vz3kd9+48C76+qK+axRo2wn6rqZVX12P7Y4Vy6/fVLd51bshb7bg0TD+DfAl+nO//zjf20twH/rh/ek+5OVrfTBZ3HDyz7xn6524DnL/Zz2VUfc91HwJvozpe/aeBx0GI/n131McpnaWAd5+Fd45bkPgJ+m+5mFrcC71zs57IrP0b4m7dPP309XVD9/cV+LrvqY4h99FS6b7vvB34ArB9Y9v/o993twCsX+7nsqo+57qP+b91PJx07PHmxn8+u+hjlszSwjjPZye4al75wSZIkSWrGUjk1TpIkSZIWjEFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5/xNarj7m2bdAUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7d56f31cf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Training\n",
    "xvals = list(dict_roi_corr_tr.keys())\n",
    "yvals = list(dict_roi_corr_tr.values())\n",
    "myfig = plt.figure(figsize=(14,7))\n",
    "plt.barh(range(len(dict_roi_corr_tr)), yvals)\n",
    "plt.yticks(range(len(dict_roi_corr_tr)), xvals)\n",
    "plt.title(\"Best correlation, within-ROI average; subject \"+ str(subid)+\" (training)\" )\n",
    "plt.axvline(x=np.mean(np.array(yvals)), color=\"gray\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0cAAAGrCAYAAAALo+xTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+cXXV95/HXW1L5lQCpBDQFiRZdV0qg7ojrDwy71lpbXfvYRkCtgttKWR6ydbe0ttUq+litsbBdtz90se5SxTZIxC7+WJW6JYordScEgalFFCiBRAUMMRANGj77xzkjh3EmuXPvnbmZ5PV8PO5jzu/zOfneOznv+Z5zbqoKSZIkSdrfPWbUBUiSJEnS3sBwJEmSJEkYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkvZKSVYkqSSL+lz/95P8xbDr6leS/53krN3MvzTJf97N/AeSPHluqtNMkhyY5B+SPH4e9vWqJJ8d0rauTPILw9iWpP2L4UjSXi/JHUm+154gb03yySTHDmm7PzeMGkcpyWlJ7upOq6p3VtWvj6qmqarqxVX1lwBJzk5y7SzXX1xVt/Wz7zZ4PdS+f76T5OokT5uyzDFJPpzkviQPJvlykpdMWaaSHN9PDQvYOcDnq+qbsOcQ26vpwn9Vfbiqfn7QbbfeBbxjSNuStB8xHElaKF5aVYuBJwDfAv5kxPX0bLren357hNS3d7fvn58C7gY+MDkjyU8C1wIPAScARwJ/DPxVktUjqHVGSQ6Y513+BvChed7nwKrqy8BhScZGXYukhcVwJGlBqarvA+uAp09Oay/9uSjJnUm+leR9SQ5u5x2Z5BNJ7m97Db6Q5DFJPgQ8Efh426PwO9PtL8nLktyQ5LtJvjF5qU6S5Umuarf59SSv66xzYZJ1SS5L8l3g7BmmPSbJ77bbvS/JR9oT9enqeG2SrybZnuS2JL/RTj8U+N/A8vY4HmhruzDJZZ31/02Sifbf4Zok/7wz744kFyS5Mcm2JJcnOWhPbZHkSe32HtOO/0WSb3fmX5bkDe3wNUl+vd3v+4Bnt7Xe39nk0rZXcHuSv0/y051t/ajXpu29+LOZlt2dqvoe8BHg5M7k/wg8APxaVX2zqr5XVX9N0/NwcZL08G/xS0k2tu+TTUku7Mz7dJLXT1n+K0n+bTv8tLY36ztJbklyeme5S5O8N8mnkjwI/Kvd7atd5zVJ/ql9T/1BOj2ks3zPPRH4aeDv2/FzgFcBv9O23cfb6cuTfDTJPUluT/IfOts4Jcl4W+u3kvyXdtbn25/3t9t6dqb0KLZtfm6SW9P0GP/ZZFskOSDJxUnubff5+vz4ZajXAL+0+5aTpEczHElaUJIcApwBXNeZvAZ4Ks0J7/E0vQNvaef9FnAXsAw4Gvh9oKrq1cCdtD1SVfXuafZ1CvBB4LeBI4DnA3e0s/+63e5yYDXwziQv6Kz+MpoQdwTw4Rmm/Qfgl4FV7Xa2An82w6F/G3gJcBjwWuCPkzyjqh4EXgxsbo9jcVVtnnIcT23rfUP77/ApmlD42M5ipwO/ADwJWAmcPUMdP1JVtwPfBX62nXQq8EAneD0fWD9lna8C5wJfams9ojP7FcDbgKXA19n9ZVGzWfZH2jD5inadSS8EPlpVD09Z/CM0AfqpPWz6QeA1NG37S8C/T/LL7by/avc5WcPTgeOAT7b1XN0uc1S73J8nOaGz7Ve2x7eEpodrxn212/5zmhDzBOBwms/DpNm8504EbquqHwJU1SU079t3t2330jYYfxz4SrufFwBvSPKidhvvAd5TVYfRBK2PtNOf3/48ot3Wl2ao4SXAM4GTaN6jk9t9Hc37/mTgGe0xTfXVdj1J6pnhSNJC8TdtL8N3aU5m/wig/Uvy64D/WFXfqartwDuBM9v1fkBzknhcVf2gqr5QVdXjPn8N+B9VdXVVPVxVd1fVP6a53+l5wBur6vtVdQPwF8CrO+t+qar+pl3vezNM+w3gTVV1V1XtBC4EVmeaS+6q6pNV9Y1qrAc+SxNGenEG8Mn2OH4AXAQcDDyns8x/q6rNVfUdmpPdk6fZznTWA6vyyA3769rxJ9EEua/0uB2AK6vqy+3J+If3UMNslgW4oH3/bKdpu25bHQlsmWadLZ35u1VV11TVTW3b3kgTRle1sz8GnJzkuHb8VW39O2lO/u+oqv9ZVT+squuBj9IE7kn/q6q+2G77+3vY12rg41V1bVU9RPNHgu77vef3HE342r6HQ38msKyq3l5VD7X3hb2fR3/+jk9yZFU9UFXXzbil6b2rqu6vqjuBv+ORdj6dJnTdVVVbae4xmmp7ewyS1DPDkaSF4pfbXoYDgdcD69sT8mXAIcCG9hKv+4FPt9OhCVFfBz6b5nK0353FPo8FvjHN9OXAZBCb9E88+i/0m6ZZb+q044CPder+KrCLpofrUZK8OMl17aVX9wO/SA8n7Z16/2lypO0h2TSl3m92hncAi3vc9nrgNJqegM/TXMq0qn19YZremN2ZTQ3TLpvmKX2Tlxe+r7PMRe37ZwXwPeCfdebdSxOgp3pCZ/5uJXlWkr9rLy3bRtM7diRA+z75JI8EhjN5pDfxOOBZk++Btm1fBXSfDveo983u9kXT1j9avqp2APd1Vu/5PUfTq7RkD4d+HM0lnd36f7+zvV+j6Xn7xyT/L1MectGDmd4TjzpOpv+8LQHun2a6JM3IcCRpQamqXVV1Jc0J3fNoTly/B5xQVUe0r8Pbm++pqu1V9VtV9WTgpcB/6lz+tqcepE00lwJNtRn4ySTdE8cn0tzo/6NSpyt/mu2/uFP3EVV1UFV1t0OSA2l6Ey4Cjm5P8j8FTN4Ls6fj2ExzEju5vdAEv7tnXKN362l6sE5rh68FnksTjtbPsE6vPXez1j6lb/LywnOnmX8n8JvAe9Lelwb8LfAr7SViXafTtNHXetj1XwFXAcdW1eE091V171X6a+AVSZ5N02v3d+30TcD6Ke+BxVX177tlz2JfW4BjJhdsj/FxnXV7es+1bgSePKVXabr38O1Ttrekqn4RoKpurapX0FwyuAZY115KOOh74FHHSfN+nuqfM7ueS0kyHElaWNJ4Gc29Jl9teybeT3MPzlHtMj81ec9DkpckOb4NBN+lCVW72s19C9jdd+d8AHhtkhe0N7L/VJKnVdUm4P8Cf5jkoCQraf5C/uHdbGs67wPeMXm5VZJl7bFN9ViaHrN7gB8meTHQfeTxt4DHJTl8hv18BPil9jh+guY+rJ3tMexRe6P7adPNq6pbacLpr9I88vm7bT2/wszh6FvAMVPueZo3VXU1TWA8p530xzSXAH4gyePbNn0F8Cbgt3u8DHMJTW/i99t71V45Zf6naALq24HLOz1qnwCemuTVSX6ifT2zc9/WbPe1Dnhpkue0/75v49Ehrdf3HFV1F3ArcEpn8tTPzJeB7yZ5Y5KD2wcl/EySZ7bb/9Uky9rjnezF2UXzXn6Y3X/+ducjwG+2n8kjgDdOs8wqmoeVSFLPDEeSFoqPJ3mAJuC8AzirqibaeW+kuXTuujRPgvtbHrls6int+APAl4A/r6pr2nl/CLy5vRzogqk7rOZxwK+lOXneRnOyP9kD8wqaS7Q209xT8tb2pHs23kPTA/DZJNtpHjLxrGnq2E5zI/1HaC51emW73uT8f6TpmbitPZblU9a/hSa8/AlNT9tLaR5E8dCeCkxyDM2/3U27WWw9cF/bKzM5HmDjDMv/H2AC+GaSPV6yNkf+iOapawdW1X00vZAHAf9AcxnafwJeXVWX97i984C3t+34Fh558AAA7f09VwI/R9PzMzl9O03QPZPmvfRNmh6WA/vZV/uZOB9YS9O7sp3mYR4720V6es91/HcefX/WB4Cnt++zv6mqXTTvp5OB22neX39B8yAIaB7yMdF+dt8DnNneN7WD5nP8xXZb/3I3NUzn/TT33d1I8z77FPBD2j98tOHswfYzLEk9S29/EJMk7Y+S/CrNJYu/N+paNHtJFtP02DylmqcLznb9A2nCxwuqarqHVuwV2t7U91XVZI/YR4EPVNWnRluZpIXGcCRJ0j4kyUuBz9H03l1M0zP0jB4vD1wQ2nup/hVN79HRNPfkXVdVbxhpYZIWPC+rkyRp3/Iymkv0NtNcVnrmvhSMWqG5n2orTc/WV3nku80kqW/2HEmSJEkS9hxJkiRJEgDTfSP2XuvII4+sFStWjLoMSZIkzbP77mu+z/hxj3vcHpbU/m7Dhg33VtWyPS/54xZUOFqxYgXj4+OjLkOSJEnz7NJLLwXg7LPPHmkd2vsl+ad+1/WyOkmSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEgCLRl3ArGzfAevHR12FJEmS+rVqbNQVSDOy50iSJEmSMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEjDkcJTkuCQbktyQZCLJubtZ9o4kRw5z/5IkSZLUr0VD3t4W4DlVtTPJYuDmJFdV1ebuQkkOGPJ+JUmSJGkgffccJVmT5LzO+IXA+VW1s510YHf7SR5I8vYkfw88u518fpLrk9yU5Gn91iJJkiRJgxrksrq1wBmd8dOBK5Icm+RGYBOwptNrdChwc1U9q6qubafdW1XPAN4LXDDdTpKck2Q8yfg927YOUK4kSZIkzazvcFRVG4GjkixPchKwtarurKpNVbUSOB44K8nR7Sq7gI9O2cyV7c8NwIoZ9nNJVY1V1diyw5f2W64kSZIk7dagD2RYB6ym6UFa253R9hhNAKe2k75fVbumrD95Cd4uhn//kyRJkiT1bNBwtBY4kyYgrUtyTJKDAZIsBZ4L3DLgPiRJkiRpzg3UW1NVE0mWAHdX1ZYkLwQuTlJAgIuq6qZhFCpJkiRJc2ngS9mq6sTO8NXAyhmWWzxlfEVneBw4bdBaJEmSJKlfQ/0SWEmSJElaqAxHkiRJkoThSJIkSZIAw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRIwhC+BnVdLDoFVY6OuQpIkSdI+yJ4jSZIkScJwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBCy0R3lv3wHrx0ddhSRJ0v7Br1DRfsaeI0mSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJwByEoyTHJdmQ5IYkE0nO7cy7I8mRSVYkuXnY+5YkSZKkfi2ag21uAZ5TVTuTLAZuTnJVVW2eg31JkiRJ0lAM1HOUZE2S8zrjFwLnV9XOdtKBu9nHAUne3/YufTbJwYPUIkmSJEmDGPSyurXAGZ3x04Erkhyb5EZgE7Bmhl6jpwB/VlUnAPcDvzLdDpKck2Q8yfg927YOWK4kSZIkTW+gcFRVG4GjkixPchKwtarurKpNVbUSOB44K8nR06x+e1Xd0A5vAFbMsI9LqmqsqsaWHb50kHIlSZIkaUbDeCDDOmA1TQ/S2u6MtsdoAjh1mvV2doZ3MTf3P0mSJElST4YRjtYCZ9IEpHVJjpm8fyjJUuC5wC1D2I8kSZIkzZmBe2uqaiLJEuDuqtqS5IXAxUkKCHBRVd006H4kSZIkaS4N5VK2qjqxM3w1sHKG5Va0g/cCP9OZftEw6pAkSZKkfg39S2AlSZIkaSEyHEmSJEkShiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJwJC+BHbeLDkEVo2NugpJkiRJ+yB7jiRJkiQJw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRKw0B7lvX0HrB8fdRWSJEmab/dvb37uTeeCfsXMPseeI0mSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJQB/hKMlxSTYkuSHJRJJzO/Me6HEblyZZPdt9S5IkSdJcWdTHOluA51TVziSLgZuTXFVVm3tZOUk/+5QkSZKkObXbnqMka5Kc1xm/EDi/qna2kw6cuo0kFye5Psnnkixrp12T5J1J1gO/2S76/CT/N8lt9iJJkiRJGrU9XVa3FjijM346cEWSY5PcCGwC1nR6jQ4Frq+qZwDrgbd21j2iqlZV1cXt+BOA5wEvAd41UwFJzkkynmT8nm1bez4wSZIkSZqN3YajqtoIHJVkeZKTgK1VdWdVbaqqlcDxwFlJjm5XeRi4vB2+jCb8TLqcR/ubqnq4qv4BOJoZVNUlVTVWVWPLDl86i0OTJEmSpN718kCGdcBqmh6ktd0ZbY/RBHDqDOtWZ/jBKfN2dobTQx2SJEmSNGd6CUdrgTNpAtK6JMckORggyVLgucAtne1N3j/0SuDa4ZYrSZIkSXNjj0+Oq6qJJEuAu6tqS5IXAhcnKZoen4uq6qZ28QeBE5JsALbx6PuVJEmSJGmv1dNjtavqxM7w1cDKGZZb3A7+wZTpp00ZP3uG9SRJkiRpJGb9JbCSJEmStC8yHEmSJEkShiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJQI9fArvXWHIIrBobdRWSJEmab7ff3Pz0XFBzyJ4jSZIkScJwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBCy0R3lv3wHrx0ddhSRJ0uz5CGppr2fPkSRJkiRhOJIkSZIkwHAkSZIkSYDhSJIkSZIAw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRJgOJIkSZIkwHAkSZIkScAA4SjJcUk2JLkhyUSSczvzHphhnUuTrG6H70hyZL/7lyRJkqRhWjTAuluA51TVziSLgZuTXFVVm4dUmyRJkiTNm556jpKsSXJeZ/xC4Pyq2tlOOnDqtpJcnOT6JJ9LsmyGTZ/fLnNTkqf1Ub8kSZIkDUWvl9WtBc7ojJ8OXJHk2CQ3ApuANZ1eo0OB66vqGcB64K0zbPfedpn3AhdMt0CSc5KMJxm/Z9vWHsuVJEmSpNnpKRxV1UbgqCTLk5wEbK2qO6tqU1WtBI4HzkpydLvKw8Dl7fBlwPNm2PSV7c8NwIoZ9n1JVY1V1diyw5f2Uq4kSZIkzdpsHsiwDlhN04O0tjuj7TGaAE6dYd2aYfrkZXm7GOz+J0mSJEkayGzC0VrgTJqAtC7JMUkOBkiyFHgucEtnu6vb4VcC1w6nXEmSJEmaGz331lTVRJIlwN1VtSXJC4GLkxQQ4KKquqld/EHghCQbgG08+n4lSZIkSdrrzOpStqo6sTN8NbByhuUWt4N/MGX62Z3hFZ3hceC02dQiSZIkScPU95fASpIkSdK+xHAkSZIkSRiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGz/BLYkVtyCKwaG3UVkiRJkvZB9hxJkiRJEoYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRJgOJIkSZIkYKE9ynv7Dlg/PuoqJEmSNEx+VYv2EvYcSZIkSRKGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkC+gxHSY5LsiHJDUkmkpzbxzauSTLWz/4lSZIkadgW9bneFuA5VbUzyWLg5iRXVdXmXlZOckCf+5UkSZKkObHHnqMka5Kc1xm/EDi/qna2kw7sbifJe5OMtz1Kb+tMvyPJW5JcC7y8nfzyJF9O8rUkpw7jgCRJkiSpH71cVrcWOKMzfjpwRZJjk9wIbALWdHqN3lRVY8BKYFWSlZ11v19Vz6uqte34oqo6BXgD8Nbpdp7knDZsjd+zbessDk2SJEmSerfHcFRVG4GjkixPchKwtarurKpNVbUSOB44K8nR7SqnJ7ke2AicADy9s7nLp2z+yvbnBmDFDPu/pKrGqmps2eFLez4wSZIkSZqNXu85WgesBh5P05P0I1W1OckEcGqSDcAFwDOramuSS4GDOos/OGW7k5fm7ZpFLZIkSZI0dL0+rW4tcCZNQFqX5JgkBwMkWQo8F7gFOIwmAG1re5JePPySJUmSJGn4euqtqaqJJEuAu6tqS5IXAhcnKSDARVV1E0CSjcAEcBvwxTmqW5IkSZKGqudL2arqxM7w1TQPXJhuubNnmL5iyvhpneF7meGeI0mSJEmaD319CawkSZIk7WsMR5IkSZKE4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSMIsvgd0rLDkEVo2NugpJkiRJ+yB7jiRJkiQJw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRKw0B7lvX0HrB8fdRWSJEmaLb+ORQuAPUeSJEmShOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkoABwlGS45JsSHJDkokk5/awzjVJxtrhB/rdtyRJkiQN26IB1t0CPKeqdiZZDNyc5Kqq2jyk2iRJkiRp3vTUc5RkTZLzOuMXAudX1c520oHdbSV5b5LxtkfpbbvZ7juSfCXJdUmO7u8QJEmSJGlwvV5WtxY4ozN+OnBFkmOT3AhsAtZ0eo3eVFVjwEpgVZKV02zzUOC6qjoJ+Dzwuul2nOScNmiN37Nta4/lSpIkSdLs9BSOqmojcFSS5UlOArZW1Z1VtamqVgLHA2d1en9OT3I9sBE4AXj6NJt9CPhEO7wBWDHDvi+pqrGqGlt2+NKeD0ySJEmSZmM29xytA1YDj6fpSfqRqtqcZAI4NckG4ALgmVW1NcmlwEHTbO8HVVXt8K5Z1iJJkiRJQzWbp9WtBc6kCUjrkhyT5GCAJEuB5wK3AIcBDwLb2p6kFw+3ZEmSJEkavp57a6pqIskS4O6q2pLkhcDFSQoIcFFV3QSQZCMwAdwGfHEO6pYkSZKkoZrVpWxVdWJn+GqaBy5Mt9zZM0w/rTO8uDO8juayPUmSJEkaib6/BFaSJEmS9iWGI0mSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJmOWXwI7ckkNg1dioq5AkSZK0D7LnSJIkSZIwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAEL7VHe23fA+vFRVyFJkiTwK1a0z7HnSJIkSZIwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSMORwlOTkJF9KMpHkxiRn7GbZa5KMDXP/kiRJktSvRUPe3g7gNVV1a5LlwIYkn6mq+7sLJTlgyPuVJEmSpIH03XOUZE2S8zrjFwIvrapbAapqM/BtYFk7/44kb0lyLfDydrWXJ/lykq8lObXfWiRJkiRpUINcVrcW6F42dzpwxeRIklOAxwLf6Czz/ap6XlWtbccXVdUpwBuAt063kyTnJBlPMn7Ptq0DlCtJkiRJM+s7HFXVRuCoJMuTnARsrao7AZI8AfgQ8Nqqeriz2uVTNnNl+3MDsGKG/VxSVWNVNbbs8KX9litJkiRJuzXoPUfrgNXA42l6kkhyGPBJ4M1Vdd2U5R+cMr6z/blrCLVIkiRJUt8GDSRrgfcDRwKrkjwW+Bjwwaq6YrdrSpIkSdJeZKBHeVfVBLAEuLuqttDcd/R84OwkN7Svk4dQpyRJkiTNqYEvZauqEzvDlwGXzbDciinjp3WG72WGe44kSZIkaT4M9UtgJUmSJGmhMhxJkiRJEoYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRJgOJIkSZIkwHAkSZIkSYDhSJIkSZIAWDTqAmZlySGwamzUVUiSJEnaB9lzJEmSJEkYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkoCF9ijv7Ttg/fioq5AkSdIg/GoW7aXsOZIkSZIkDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBMxBOEpycpIvJZlIcmOSMzrzrkky1g4/MOx9S5IkSVK/Fs3BNncAr6mqW5MsBzYk+UxV3T8H+5IkSZKkoRio5yjJmiTndcYvBF5aVbcCVNVm4NvAshnWf0eSryS5LsnRg9QiSZIkSYMY9LK6tcAZnfHTgSsmR5KcAjwW+MY06x4KXFdVJwGfB1433Q6SnJNkPMn4Pdu2DliuJEmSJE1voHBUVRuBo5IsT3ISsLWq7gRI8gTgQ8Brq+rhaVZ/CPhEO7wBWDHDPi6pqrGqGlt2+NJBypUkSZKkGQ3jnqN1wGrg8TQ9SSQ5DPgk8Oaqum6G9X5QVdUO7xpSLZIkSZLUl2EEkrXA+4EjgVVJHgt8DPhgVV2x2zUlSZIkaS8x8KO8q2oCWALcXVVbaO47ej5wdpIb2tfJg+5HkiRJkubSUC5lq6oTO8OXAZfNsNxpneHFneF1NJfnSZIkSdJIDP1LYCVJkiRpITIcSZIkSRKGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEnAkL4Edt4sOQRWjY26CkmSJEn7IHuOJEmSJAnDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJErDQHuW9fQesHx91FZIkSZpv929vfu6v54J+nc28sOdIkiRJkjAcSZIkSRJgOJIkSZIkwHAkSZIkSYDhSJIkSZIAw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRIwB+EoyaeT3J9FpY3sAAAJtklEQVTkE3tY7pokY8PevyRJkiT1Yy56jv4IePXuFkhywBzsV5IkSZL61nc4SrImyXmd8QuT/FZVfQ7YPs3ydyR5S5JrgZe3k1+e5MtJvpbk1H5rkSRJkqRBDdJztBY4ozN+OnDFHtb5flU9r6rWtuOLquoU4A3AW6dbIck5ScaTjN+zbesA5UqSJEnSzPoOR1W1ETgqyfIkJwFbq+rOPax2+ZTxK9ufG4AVM+znkqoaq6qxZYcv7bdcSZIkSdqtRQOuvw5YDTyepidpTx6cMr6z/blrCLVIkiRJUt8GDSRrgfcDRwKrBi9HkiRJkkZjoKfVVdUEsAS4u6q2ACT5As29Ry9IcleSFw1epiRJkiTNrYEvZauqE6eMT/vUuapaMWX8tM7wvcxwz5EkSZIkzYe5+J4jSZIkSVpwDEeSJEmShOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEjCEL4GdV0sOgVVjo65CkiRJ8+32m5ufngtqDtlzJEmSJEkYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkoCF9ijv7Ttg/fioq5AkSdJ8u39789Nzwb3LPvZodXuOJEmSJAnDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBcxSOknw6yf1JPjFl+jVJxtrhB+Zi35IkSZLUj7nqOfoj4NVztG1JkiRJGrqBwlGSNUnO64xfmOS3qupzwPYe1n9Hkq8kuS7J0YPUIkmSJEmDGLTnaC1wRmf8dOCKHtc9FLiuqk4CPg+8brqFkpyTZDzJ+D3btg5UrCRJkiTNZKBwVFUbgaOSLE9yErC1qu7scfWHgMl7kjYAK2bYxyVVNVZVY8sOXzpIuZIkSZI0o0VD2MY6YDXweJqepF79oKqqHd41pFokSZIkqS/DCCRrgfcDRwKrhrA9SZIkSZp3Az+trqomgCXA3VW1BSDJF2juPXpBkruSvGjQ/UiSJEnSXBrKpWxVdeKU8VNnWO60zvDizvA6msvzJEmSJGkk5up7jiRJkiRpQTEcSZIkSRKGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAFg06gJmZckhsGps1FVIkiRpvt1+c/PTc0HNIXuOJEmSJAnDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJAKSqRl1Dz5JsB24ZdR3qy5HAvaMuQn2x7RY222/hsu0WLttu4bLtFq5u2x1XVcv62cii4dUzL26pqrFRF6HZSzJu2y1Mtt3CZvstXLbdwmXbLVy23cI1rLbzsjpJkiRJwnAkSZIkScDCC0eXjLoA9c22W7hsu4XN9lu4bLuFy7ZbuGy7hWsobbegHsggSZIkSXNlofUcSZIkSdKcMBxJkiRJEntROEryC0luSfL1JL87zfwDk1zezv/7JCs6836vnX5LkhfNZ93qv+2SvDDJhiQ3tT//9XzXvr8b5HPXzn9ikgeSXDBfNasx4O/MlUm+lGSi/fwdNJ+17+8G+J35E0n+sm2zryb5vfmuXT213/OTXJ/kh0lWT5l3VpJb29dZ81e1oP+2S3Jy53fmjUnOmN/KNcjnrp1/WJK7k/zpHndWVSN/AQcA3wCeDDwW+Arw9CnLnAe8rx0+E7i8HX56u/yBwJPa7Rww6mPaX14Dtt3PAsvb4Z8B7h718exPr0HarjP/o8AVwAWjPp796TXg524RcCNwUjv+OH9nLpi2eyWwth0+BLgDWDHqY9qfXj223wpgJfBBYHVn+k8Ct7U/l7bDS0d9TPvLa8C2eyrwlHZ4ObAFOGLUx7S/vAZpu8789wB/Bfzpnva3t/QcnQJ8vapuq6qHgLXAy6Ys8zLgL9vhdcALkqSdvraqdlbV7cDX2+1pfvTddlW1sao2t9MngIOSHDgvVQsG+9yR5Jdp/nOfmKd69YhB2u7ngRur6isAVXVfVe2ap7o1WNsVcGiSRcDBwEPAd+enbLX22H5VdUdV3Qg8PGXdFwFXV9V3qmorcDXwC/NRtIAB2q6qvlZVt7bDm4FvA8vmp2wx2OeOJP8COBr4bC8721vC0U8Bmzrjd7XTpl2mqn4IbKP5i2cv62ruDNJ2Xb8CbKyqnXNUp35c322X5FDgjcDb5qFO/bhBPndPBSrJZ9pLEH5nHurVIwZpu3XAgzR/tb4TuKiqvjPXBetRBjnn8HxltIby75/kFJrei28MqS7tWd9tl+QxwMXAb/e6s0WzKm3uZJppU58xPtMyvayruTNI2zUzkxOANTR/0db8GaTt3gb8cVU90HYkaX4N0naLgOcBzwR2AJ9LsqGqPjfcEjWDQdruFGAXzWU9S4EvJPnbqrptuCVqNwY55/B8ZbQG/vdP8gTgQ8BZVfVjPRSaM4O03XnAp6pqU6/nK3tLz9FdwLGd8WOAzTMt015ScDjwnR7X1dwZpO1IcgzwMeA1VeVfYebXIG33LODdSe4A3gD8fpLXz3XB+pFBf2eur6p7q2oH8CngGXNesSYN0navBD5dVT+oqm8DXwTG5rxidQ1yzuH5ymgN9O+f5DDgk8Cbq+q6Idem3Ruk7Z4NvL49X7kIeE2Sd+1uhb0lHP0/4ClJnpTksTQ3oF41ZZmrgMknu6wG/k81d1hdBZzZPt3nScBTgC/PU90aoO2SHEHzi+b3quqL81axJvXddlV1alWtqKoVwH8F3llVe34CjIZlkN+ZnwFWJjmkPfFeBfzDPNWtwdruTuBfp3Eo8C+Bf5ynutXopf1m8hng55MsTbKU5mqJz8xRnfpxfbddu/zHgA9W1RVzWKOm13fbVdWrquqJ7fnKBTRt+GNPu5u60l7xAn4R+BrNNZxvaqe9Hfg37fBBNE/F+jpN+HlyZ903tevdArx41Meyv736bTvgzTTXz9/QeR016uPZn16DfO4627gQn1a3oNoO+FWaB2ncDLx71Meyv70G+J25uJ0+QRNof3vUx7I/vnpov2fS/KX7QeA+YKKz7r9r2/XrwGtHfSz726vftmt/Z/5gyvnKyaM+nv3pNcjnrrONs+nhaXVpF5YkSZKk/dreclmdJEmSJI2U4UiSJEmSMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQB8P8BNMk/WF38GsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7d5363b160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Testing\n",
    "xvals = list(dict_roi_corr_te.keys())\n",
    "yvals = list(dict_roi_corr_te.values())\n",
    "myfig = plt.figure(figsize=(14,7))\n",
    "plt.barh(range(len(dict_roi_corr_te)), yvals, color=\"pink\")\n",
    "plt.yticks(range(len(dict_roi_corr_te)), xvals)\n",
    "plt.title(\"Best correlation, within-ROI average (testing)\")\n",
    "plt.axvline(x=np.mean(np.array(yvals)), color=\"gray\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is clear that such a simple toy example is wildly inefficient---the output is basically junk. Some serious modifications to the following factors will be required. Keep them in mind when reading the major tasks below.\n",
    "\n",
    "- Setting of `w_init`.\n",
    "- The size and range of the $\\lambda$ grid.\n",
    "- Number of iterations `t_max`.\n",
    "- All of the filter bank parameters (in `myparas`), especially `freqs`, `dir`, `sdev`.\n",
    "- The quantity and variety of filters used to generate features; both fine and coarse grids, and a wide variety of frequencies and orientations is likely necessary (see the \"major tasks\" below for more ideas)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "0. Focus on \"the early visual areas\" looked at by Nishimoto et al. (V1, V2, V3, V3A, and V3B) in both hemispheres. Train a model for each voxel in these regions, and compute the performance (on test data) for the best lambda value (determined on the *training* set). Average the error/correlation over all the voxels in these regions.\n",
    "\n",
    "0. Complete the above exercise for each subject. Is there much difference in performance between subjects? How does your best model perform against the cited work?\n",
    "\n",
    "0. Another approach is to capture temporal delays in the feature vectors. For example, if our original feature vectors are $x_{i}$ for $i=1,\\ldots,n$, then re-christen the feature vectors as $\\widetilde{x}_{i} = (x_{i},x_{i-1})$ for a delay of *one* step (here, one second) for all $i>1$ (we lose one data point, now $n-1$ total). The dimension grows from $d$ to $2d$. Analogously, for a delay of $k$ steps, this would be $\\widetilde{x}_{i} = (x_{i},x_{i-1},\\ldots,x_{i-k})$ for all $i>k$, and lose $k$ data points for a total of $n-k$ now. Try several temporal delays; which seem to work best?\n",
    "\n",
    "0. How does performance depend on ROI looked at? Which ROI saw comparatively good/bad performance? Provide visuals to highlight performance in each region.\n",
    "\n",
    "0. Be sure to experiment with the learning algorithm parameters (number of iterations, size and range of $\\lambda$ grid, etc.). What strategies did you find particularly effective? If you made any modifications to the algorithm, describe them.\n",
    "\n",
    "0. Make a large data set of stimulus from any two subjects, and use the third subject's data as a evaluation of *inter-subject generalization* ability. How does performance compare with the more standard by-subject approach? Are there subjects that are particularly difficult to predict for? In contrast to this, considering the by-subject training approach we have considered thus far, how does our interpretation of *generalization* change?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## References:\n",
    "\n",
    " - Nishimoto, Shinji, et al. \"Reconstructing visual experiences from brain activity evoked by natural movies.\" Current Biology 21.19 (2011): 1641-1646.\n",
    " - Description of dataset vim-2 (visual imaging 2), at CRCNS - Collaborative Research in Computational Neuroscience. https://crcns.org/data-sets/vc/vim-2/about-vim-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}