{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shared response model\n",
    "\n",
    "Authors: Javier Turek ([javier.turek@intel.com](mailto:javier.turek@intel.com)), Samuel A. Nastase ([sam.nastase@gmail.com](mailto:sam.nastase@gmail.com)), Hugo Richard ([hugo.richard@ens-lyon.fr](mailto:hugo.richard@ens-lyon.fr))\n",
    "\n",
    "This notebook provides interactive examples of functional alignment using the shared response model (SRM; [Chen et al., 2015](https://papers.nips.cc/paper/5855-a-reduced-dimension-fmri-shared-response-model)). BrainIAK includes several variations on the SRM algorithm, but here we focus on the core probabilistic [`SRM`](https://brainiak.org/docs/brainiak.funcalign.html#brainiak.funcalign.srm.SRM) implementation. The goal of the SRM is to capture shared responses across participants performing the same task in a way that accommodates individual variability in response topographies ([Haxby et al., 2020](https://doi.org/10.7554/eLife.56601)). Given data that is synchronized in the temporal dimension across a group of subjects, SRM computes a low dimensional *shared* feature subspace common to all subjects. The method also constructs orthogonal weights to map between the shared subspace and each subject's idiosyncratic voxel space. This notebook accompanies the manuscript \"BrainIAK: The Brain Imaging Analysis Kit\" by Kumar and colleagues (2020).\n",
    "\n",
    "The functional alignment ([`funcalign`](https://brainiak.org/docs/brainiak.funcalign.html)) module includes the following variations of SRM:\n",
    "* [`SRM`](https://brainiak.org/docs/brainiak.funcalign.html#brainiak.funcalign.srm.SRM): A probabilistic version of SRM\n",
    "* [`DetSRM`](https://brainiak.org/docs/brainiak.funcalign.html#brainiak.funcalign.srm.DetSRM): A deterministic version of SRM\n",
    "* [`RSRM`](https://brainiak.org/docs/brainiak.funcalign.html#brainiak.funcalign.rsrm.RSRM): Robust SRM for better filtering idiosyncratic components and outliers in data\n",
    "* [`SSSRM`](https://brainiak.org/docs/brainiak.funcalign.html#brainiak.funcalign.sssrm.SSSRM): Semi-supervised SRM for   labeled data  \n",
    "* [`FastSRM`](https://brainiak.org/docs/brainiak.funcalign.html#brainiak.funcalign.fastsrm.FastSRM): A faster version of SRM with reduced memory demands"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annotated bibliography\n",
    "1. Chen, P. H. C., Chen, J., Yeshurun, Y., Hasson, U., Haxby, J., & Ramadge, P. J. (2015). A reduced-dimension fMRI shared response model. In C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, R. Garnett (Eds.), *Advances in Neural Information Processing Systems, vol. 28* (pp. 460-468). [`link`](https://papers.nips.cc/paper/5855-a-reduced-dimension-fmri-shared-response-model) *Introduces the SRM method of functional alignment with several performance benchmarks.*\n",
    "\n",
    "2. Haxby, J. V., Guntupalli, J. S., Nastase, S. A., & Feilong, M. (2020). Hyperalignment: modeling shared information encoded in idiosyncratic cortical topographies. *eLife*, *9*, e56601. [`link`](https://doi.org/10.7554/eLife.56601) *Recent review of hyperalignment and related functional alignment methods.*\n",
    "\n",
    "3. Chen, J., Leong, Y. C., Honey, C. J., Yong, C. H., Norman, K. A., & Hasson, U. (2017). Shared memories reveal shared structure in neural activity across individuals. *Nature Neuroscience*, *20*(1), 115-125. [`link`](https://doi.org/10.1038/nn.4450) *SRM is used to discover the dimensionality of shared representations across subjects.*\n",
    "\n",
    "4. Nastase, S. A., Liu, Y. F., Hillman, H., Norman, K. A., & Hasson, U. (2020). Leveraging shared connectivity to aggregate heterogeneous datasets into a common response space. *NeuroImage*, *217*, 116865. [`link`](https://doi.org/10.1016/j.neuroimage.2020.116865) *This paper demonstrates that applying SRM to functional connectivity data can yield a shared response space across disjoint datasets with different subjects and stimuli.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents\n",
    "- [](#example-fmri-data-and-atlas)\n",
    "- [](#estimating-the-srm)\n",
    "- [](#between-subject-time-segment-classification)\n",
    "- [](#summary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary python modules\n",
    "from glob import glob\n",
    "import nibabel as nib\n",
    "import numpy as np\n",
    "from nilearn.plotting import plot_stat_map\n",
    "from scipy.stats import zscore\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import brainiak.funcalign.srm\n",
    "from brainiak.fcma.util import compute_correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example fMRI data and atlas\n",
    "To work through the SRM functionality, we use an fMRI dataset collected while participants listened to a spoken story called \"[I Knew You Were Black](https://themoth.org/stories/i-knew-you-were-black)\" by Carol Daniel. These data are available as part of the publicly available [Narratives](https://github.com/snastase/narratives) collection ([Nastase et al., 2019](https://openneuro.org/datasets/ds002345)). Here, we download a pre-packaged subset of the data from Zenodo. These data have been preprocessed using fMRIPrep with confound regression in AFNI. We apply the SRM to a region of interest (ROI) comprising the \"temporal parietal\" network according to a cortical parcellation containing 400 parcels from Schaefer and colleagues ([2018](https://doi.org/10.1093/cercor/bhx179))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-12-07 02:06:15--  https://zenodo.org/record/4300825/files/brainiak-aperture-srm-data.tgz\n",
      "Resolving zenodo.org (zenodo.org)... 137.138.76.77\n",
      "Connecting to zenodo.org (zenodo.org)|137.138.76.77|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 9927319659 (9.2G) [application/octet-stream]\n",
      "Saving to: ‘brainiak-aperture-srm-data.tgz’\n",
      "\n",
      "brainiak-aperture-s 100%[===================>]   9.25G  4.37MB/s    in 13m 37s \n",
      "\n",
      "2020-12-07 02:20:30 (11.6 MB/s) - ‘brainiak-aperture-srm-data.tgz’ saved [9927319659/9927319659]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Download and extract example data from Zenodo\n",
    "!wget -nc https://zenodo.org/record/4300825/files/brainiak-aperture-srm-data.tgz\n",
    "!tar --skip-old-files -xzf brainiak-aperture-srm-data.tgz\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get filenames for example data and atlas\n",
    "data_fns = sorted(glob('brainiak-aperture-srm-data/sub-*_task-black_*bold.nii.gz'))\n",
    "atlas_fn = 'brainiak-aperture-srm-data/Schaefer2018_400Parcels_17Networks_order_FSLMNI152_2.5mm.nii.gz'\n",
    "\n",
    "# Load in the Schaefer 400-parcel atlas\n",
    "atlas_nii = nib.load(atlas_fn)\n",
    "atlas_img = atlas_nii.get_fdata()\n",
    "\n",
    "# Left temporal parietal ROI labels\n",
    "parcel_labels = [195, 196, 197, 198, 199, 200]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in functional data and mask with \"temporal parietal\" ROI\n",
    "data = []\n",
    "for data_fn in data_fns:\n",
    "    voxel_data = nib.load(data_fn).get_fdata()\n",
    "    \n",
    "    # Take union of all parcels (brain areas) comprising the full ROI \n",
    "    roi_data = np.column_stack([voxel_data[atlas_img == parcel, :].T\n",
    "                                for parcel in parcel_labels])\n",
    "    data.append(roi_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAADJCAYAAAAHFcoVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deZgU5bn2754ZRmSRTVBkURYVBAWjJhfGHNEIHENEIUD8XFgURRARVHKMghvESDYkcADRsGpEERFcQFCDmhijqBg1JqACAgqyKTvDMPP9Mefurr6rnumefQae33VxNd1Ty1tVb1W99/tssfz8/Hw4juM4jlPpyKjoBjiO4ziOE42/pB3HcRynkuIvacdxHMeppGQV9seOHTti69at5dUWxylVGjZsiFWrVlV0MxzHcYpNoUraX9BOVcb7r+M4VR2f7nYcx3GcSoq/pB3HcRynkuIvacdxHMeppBTqOGaxadMm3HLLLXj22WfTWr5379646667cNxxx2HYsGH4/PPP0bx5c7z66quRy9etWxfdunXDU089VZzmVWrmz5+PdevWYdSoUeW2zwkTJqBx48a48sor01r+pJNOwrnnnovFixentfxtt92GXr164YILLoj8+9tvv41mzZol/bZ//35s3LgRTzzxBB599NGkv9WsWRNDhgzBZZddhqZNm2LXrl1YuXIlpk2bhvfeey++XNOmTfGPf/wDV1xxBd5999202uo4TnpMmDABAFC9evWkz2OOOQYAkJmZCQA4fPgwAGDXrl0AgJtuuintfUyaNAkAUKdOHQBAtWrVkrZ54MABAMCgQYOKeRRVn3JR0vfccw9WrFiBCy+8ECtWrMDMmTPRoUMHc/m77roLvXv3Lo+mHRXcc889GDx4cNrL//73v8dFF11Uqm2YPHkyOnbsGP/3k5/8BG+++Sbuu+8+9OjRI75c3bp1sXjxYlx++eX4/e9/jwsvvBADBgzAzp07sWDBAvTt27dU2+U4jlOZKZaSLirHHXcc3n77bWzatCmt5WOxWBm36Ohi9+7dRVq+LM7/vn37krytt27dijFjxuCiiy7CZZddFlftY8eORa1atdCtWzd8++23AICNGzfiww8/xNatW/Hggw/i3Xffxdq1a0u9jY5zNDJz5kwACYWclVXwWmjZsmXS79nZ2QASapffyb59+wAAf/7znwEAO3bsCC23f/9+AECDBg0AAO3atQMA1KpVC0BCnRMq6jfeeAMAsHfvXgAJ1b59+3YAwNChQ4tyyFWKUlHSV111Fd544w189tlneOWVV9CnTx8ABdORmzZtQrVq1TBhwgS8/fbbmD9/Plq0aIHbb78db7/9dmhbt912G6666iqcf/752LRpE5o2bVroPgCgU6dO+Pzzz9G5c2e8+eab+OyzzzBv3jw0btwY48aNw6effopVq1bh5ptvjq8zYcIEPPzwwxg/fjxWr16N9957DyNGjEhqy+mnn445c+bgk08+wccff4yJEyeiXr168b9v2rQJd9xxB1auXImVK1eiYcOGOOOMMzBnzhz861//wtq1a/HGG2+kPSvQqVMnrF+/Ht27d8fbb7+N1atXY/bs2TjppJPiyzRt2hSPPPIIPvroI6xbtw5vv/120vTShAkTMHXqVMyfPx+ffvop+vXrhwkTJmDevHlJx/XEE0/gs88+w8qVKzF+/Hgcd9xx8fV/9KMfoW/fvvFBVd26dfH73/8e77//PtatW4f33nsPo0ePLpWX+aFDh5CbmwsAqFevHnr06IFHH300/oIOMnHiRBw6dAhXXXVViffrOI5TFSixku7Xrx9uv/123HXXXfj4449xzjnnYNy4cQCABQsWoGPHjli5ciUeeOABLFy4EHl5eViyZAleeukl/O///m9oe9OmTUOLFi3QvHlzDBo0CNu3by90H/PnzwdQMFq74447cPPNN6NatWqYPXs2li9fjieeeALdu3dHr169cNddd2HZsmVYs2YNAKBHjx5YsmQJunfvjjPOOAO/+c1vkJubi8mTJ6Np06Z47rnnsHz5cvTq1Qt16tTBr371K8ybNw+XXnop8vLyABQMHq699lpUq1YNe/bswfLly7F8+XJcdtllAIDBgwfjN7/5DVasWIFt27alviBZWfjlL3+JUaNGYceOHXjwwQfx+OOPo0uXLjh8+DBmzZqFjRs3ok+fPjhw4AB69+6NMWPG4M0338Qnn3wSP64xY8bgzjvvxK5du3D22WfHt3/iiSdiwYIFmDdvHu655x7UqVMHo0ePxmOPPYa+ffvinnvuQfPmzfHNN9/gnnvuAVDwcqxfv3582vniiy/GuHHj8O677+Lll18uVr+pXr06BgwYgNNOOw0PPvgggILkOVlZWVi5cmXkOjk5OXjvvfdw7rnnFmufjuMknpn169cHALRq1QoA4s80QkVN5azfqbA5WKfNmmr40KFDABIKGwBq164NAGjevHlSG7gut8W28JNKm2KiYcOGAIAmTZoAAJ5//nkAwNdffw0AuPHGG9M6F1WBEr+khw8fjj/84Q948cUXAQDr169H06ZNccstt2D+/PnxKc7du3fHpz8OHz6MvXv3xr8H2bdvHw4cOICcnJz4uqn2AQAZGRkYP348/vnPfwIA/va3v6FDhw749a9/DaDAJjpy5Eicfvrp8Zf0jh07MHLkSOTk5GDNmjU49dRTMXDgQEyePBn9+/fHrl27cNttt8WV3pAhQ/D666/joosuiju9Pf300/GXY4MGDTB9+nTMmDEj7vAwadIkXH311WjZsmVaL2kAuP/++/Hmm2/Gj/2tt97CBRdcgH/84x+YP38+Fi1ahM2bNwMA/vCHP2D48OFo06ZNvB1btmzBjBkzIrfdr18/rF+/Pj7IAQqmit577z2cc845eO+993Do0CEcOHAgfv7/8pe/4K233sLq1asBALNnz8bQoUPRtm3btF/St956a3wmIxaL4ZhjjsGnn36KIUOGYPny5QASziM7d+40t7Nz5874De44jnOkU6KXdP369dG4cWOMGTMGd999d/z3zMxMZGVloVq1avHRVFnug6xbty7+/3379uHLL7+Mf+dLM2gf+eCDD5CTkxP//v7772PkyJGoV68eTj/9dKxatSr+ggaAzz77DNu3b8fpp58ef0kH97F9+3bMmTMHffr0Qfv27dGiRQucccYZ8famy9///vf4/9evX49t27ahTZs2eP311zFz5kz06NEDHTt2RIsWLdCuXTtkZmYmbT/YJqV9+/Zo3759/IUbpHXr1kne02TOnDno1q0brrrqKrRs2RJt27bFSSedhIyM9K0ls2fPxuzZs5GZmYmuXbtixIgReOqpp5I8yPly5qg5ijp16sTtUI7jpOaPf/wjAOCUU04BkFCxNWvWBID4fXzw4EEACQVMVatK+thjjwWQUNJ8RlL11qhRA0Bi0E37MVAwkwckFDW3xbYQqy36HFXPc+579uzZAID+/ftHn5QqRIle0jyBY8aMSXqxkOALrjz2ofvT6ZtU7WMHyM/Pj3cSJSMjI2k9vvwBoFGjRnj++eexefNmLF++HK+88go2b96MpUuXFtqOdNqVn5+PY489Fs899xwyMzPx4osv4q233sIHH3yAf/zjH0nLB9uk5OTk4PXXX49PZQexXn5z585Fq1atsHDhQjzzzDNYtWpVkcPjvv322/ggaurUqcjLy8PYsWOxfft2LFq0CACwatUqHDx4EN///vfjswJBqlWrhrPPPvuIDM1zHMeJokQv6d27d+Prr79Gs2bN8OSTT8Z/v/baa9GuXTvceeedxdpufn5+me8DKFCVsVgsvr/vfe972LhxI7799lusXr0avXv3RlZWVvyleeqpp6JevXqRKhQALr30UtSsWRM9e/aMDxAuvPBCAEXzmD7zzDPjcb8tW7ZEvXr18PHHH+P8889H+/bt0a5du7hjVatWrYqk0levXo2ePXti48aN8eNq1qwZxo0bhwcffBD/+c9/ks7/aaedhs6dO6Nbt274+OOPARQo3UaNGpXIceyRRx5Bt27d8OCDD+Ktt97C1q1b8d1332H+/PkYMmQInnvuudC095AhQ1CjRo2496jjHEnQdBd0ii0OEydOBJCYNaSirVu3LoCEauUnlTKFCW3IfIbx71SrXI/PCZ0t5fNI46qBhIKmKtdPtpnPFlXp3DbbxE/OBnB5zsYxBwO3VxXjrUvs3T1x4kTceOONuPrqq3HyySfjiiuuwL333otvvvnGXGfPnj1o2bIlTjjhhMi/7927F40bN0azZs2QmZlZrH2kQ8uWLXH//fejVatW6NWrF66//npMnToVQEFYQu3atfGHP/wBp512Gs477zxMnjwZn3zyCf76179Gbm/Hjh2oXbs2unfvjiZNmqBLly4YP348gHC4QmH8+te/xrnnnouzzjoLEydOxAcffIC///3vcaXbs2dPNGnSBBdccAGmTZsGIPlGKIyZM2eiTp06mDBhAtq0aYOzzjoLU6ZMQYsWLfDFF18AKLg+zZo1Q5MmTfDdd9/h0KFD8cQi55xzDmbMmIHq1asX6ZiiGDVqFKpXr46xY8fGfxs7diy2bNmCxYsXx89ju3bt8MADD+C2227DXXfdFW+n41R1NmzYgFNPPRVAwcu5pC9op2xp0qQJfvzjH6NLly7o2rUrfvSjH+EnP/kJPvzwwzLbZ4kdx+bOnYvs7GwMGTIEY8eOxebNm/Hwww9j8uTJ5jrTp0/H2LFjceGFF+Kss85KUm4A8NRTT6Fbt25YsWIFevXqVax9pMM777yDGjVqYOnSpdi+fTseeughzJo1CwCwbds2/L//9/8wevRovPTSS9i/fz+WLVuGcePGmdP4zz//PDp27IixY8eiZs2aWLduHSZMmIBbbrkFHTt2xIoVK9Jq1/z58zFt2jTUqlULr7zyCsaMGYP8/HysWrUKDzzwAG6++Wbcfffd2LRpE5588kns3LkTHTp0wNy5c1Nue+vWrbjyyitx991344UXXsCBAwfwt7/9DYMHD46PiGfPno1Jkybh9ddfR6dOnXDbbbfh9ttvx6BBg7B58+b4lH7Hjh3TOh6Lzz//HJMmTcKoUaOwYMECLF++HHv27EHv3r1x3XXXYeTIkTjllFOwZ88evPPOO+jVqxfef//9Eu3TcY5UONNIb20+p/ipPiRUlxxs85Oqls8DqlVV3sy/wO3zOc5Pbp/bC/6ff9PZOLV/U0FT3WuMNtU/l9Njbdy4MYCEV/iyZcsAIO54269fPxSV+fPnx73SgYKIpNGjR8c9zEubWL6+IQPQvf1IpKipMsuDTp064ZlnnsG5554bDyVwSka6CXQcpyLYsGEDLr744njESUngS5pTyvqSptMWp705JUxnK8KEI+m+pLk898PvjA4J3oNMkHL88ccDSLw82QadSv/uu+8ApH5Jsy00A3LfXJ77ob9OcV/STZo0wUcffRR/Sefm5uL+++/HunXr0hJJxaFcMo45juM4ZQNrKFA1UoXypblnzx4AYaXLbF6McFGvbC6n8dFcnutzP/yuDrtBBc//85NKWttEezYHFDTncT3ugxnIODDgsfJY+DKlFzm3z7/TafXyyy9HuvTp0wexWAw7duzAMcccg0suuSSe57ws8Je04ziO46QJp7s/+ugjXHvttTj//PPjMwNlwVH7kh45cmRFNyHE3//+9yPaxOA4TulBFdioUSMAYedRzfyltmJVsVyOKpNTyqpeuRwVtSpptaAGo09UQROuw33wOxU1VTz/rlPt/OS+qJzVC5zHxr9zyn/JkiUACiJ00uXMM8/Efffdh5EjR6J9+/ahSn+lhdeTdhzHcZxicMUVV6Bjx4647777ymwfR62SdhzHqQzs27cvHoZFFi9ejLZt20Yu//TTTwNAvPAOnaI0VphqVXNsE1XSqoDVG5xKWR3LVEmroo7Kp2DtS9tE1M5ORc19cT0qZFX//M6/k2DGSqB4Nupx48ahS5cuWLFiBTp37pz2euniL2nHcZwKolmzZh6BUIWIulatW7cu09K5hb6kGzZsmFQD2HGqEqyU4zhHAsyexRzcajOmSqR6ZHiS2pxV6Voe14TbY/gSw6E0v7aGfKkXePBvXEdTGGtsNW3M/NS28nfOFmi8tHqo08atSp3LM8FWZcr9XehLetWqVeXVjhLzpz/9CUDiYrED8eQzPIGOWawLzYvGjsROQ1d+xit/9dVXABKdoLBSaMwCpgnnOT3FjsB4RnY0TiOxWhZjJ5mgXp0nhg0bZp8Qx3Ecp8rj092O4ziVgN69ewMAnnnmmci/c1BPFUilysG7VpRSpa2qVZW05tFWBU6PagoYKmpVzBRIFBfBSoNch23j3/hdE6MwTppCRu3sFFlcnttXBc3lrVhvwuUppB5++GEAwIgRI1BRlNtLurSSx5Pp06cDSIQftGjRAkDiYrI4AzPQWInceZE1qbw6IvAGYQ1sTofQaQNIKGCdftJi6ERLv1lhB2wLt8u/z5w5E0Di5hsyZEjEmXIqgtLu747jHJ2U6Us6mPLOH1aVnyZNmqBNmzbIyMhALBbD/v37Ubt2bfz6179Ghw4dKrp5VQqW0/R+75QEDsSBRFpPCgw10alqpIDQylEUIppJzIqL5vaZopOKmtm+qEpV4VP4BG3cqrap0vldY6/ZdraVbWNbOWvANmrlLgogtderfZ2fhObQ1q1bo6KpctPdfPgxiTyD0bUik14MDZK3wgPYafjJTqGdgRc1WGyD6/DmsDq9tsUK7Od6tGlzu1TSWs6NuXtZLas4NuvyTh7vOI7j2FS5l7RTfuTm5mLTpk1xu5DjOOUDK/w1b948/htNburJrAqWdl0KC+biViWtdaK5Xc3+RcWsXt2E61Ptqp03WG+aokrjlYmKK6p2rsc2cl/8ncuprZrChr+rctZKXro+mTNnTvz/xamcVRIq/Ut6xowZABJTPfTO1pNP2CH4u/7dUtJUreoswU5DVazVVHhDAOFpI07VaPJ4q9Slhk5wOd6cDRo0AJC4+fSmYMfjDUCP9+uvvz5yf1GUd/J4x3Ecx6bSv6Sd8qW8k8c7jhOmadOmABIDdCAcO0wxoFWn1GSnzq9cXmOQVcBQJFCQqGOrVRNaRUmUkNGqVlbWM7Vd85i5HLdDG7Lm7qao4jFrtSzdrjr58hgoEiuCSvuSpn2V9Uc5ZcOTyIvAk6vKV50lrM6g9l92OPXy1qkg7jcY1qCqXVPXaYC9hab1o3LmTcYpHu1IWrSd+6fzycCBAwvdb5DySh7vOI7j2FTal7RT8VxxxRWYN28e7rvvvvjUueM4ZcfWrVuxYMGCeJxucOCvQkKVMFFlTDRPNVHbNrGEjOb01u2rzTsoZFT8qFOt1SZLaXM7VMwUNNpW2p75SYWudnmiJsmgDf3xxx8HAFxzzTUoD8r8JV3c5PG0PdNpScunsfOq/ZeKWqdVtBMo2mmsT7V1B7dreXFrMXSdJtKk8uwgVMTseBp3rfCYubxOVzG2nFM+t956q31C/o+yTh7vOI7j2JTpS9qTx1ctKiJ5vOM4wJQpU7Bp0ybUqFED9erVC4kPIJyv2sprrdm5LNSLW5fndigKGCetlaVUuFAEqOgAwjZmbbu1b03spI6/Gk9t2cXpAa/nLuiBHjzGYLY0Ut7RLpVmupsu7jTQ0/6qtme9mGoz1sTtXJ+/W4pZp1FUeattW78DiQ6idm1ryiaVbVrt79ph1Q6vqfBou+ZUEKd62DEZ5uE5wB3HcSon0cYFx3Ecx3EqnHJT0lHJ46dOnRr/v7q+a7iBTl+oMwOVNG3XWhKN0x26fZ0+sbavqfiilLTaybUIunp5a8Fyna5SBa3niFNHuh0eG+GshFb8YnsZiw4A1113HRzHKV/q16+PatWqIRaLITMzM/TMABL3tzqApQqF0ulp/bu1vmZb5CdhO6xnblQRi1QOYPo819Sm1uyjtT0rLwa/azlPfQ9oYpZg28qrnKUracdxHMeppJSrTXrz5s2YM2dOfCQSHJmp4tVC4pYHtI6QNPG6JmpXGzfRguXcnuW4oEo/uE8N9mfgPEeWtJfrKE7Pi84O8LsmtNcKX0RHgTwG5uamzTo4Ml+4cCGAhBOZ26vLh0mTJgEIqxmdsaEq4v3h16dqQ7+Qli1bJjl86fMISNzvXE5TWGqCEP1UPxvta1YEDLfLpEacmdNnMisE6nMu+Ixkf1bHMJ1F5DOQzz7LYYykynuhznSWf5BG4kQpd26bz8+yxpW04ziO41RSykVJz549G5s3b0ZmZiZq1aoVHx0FRylaBo2jMrWL8LumviOqqLUilebX5shJPaXVDqxx0ZrRLNhGqnW1e+sIU2cNuC/m6ObvHLGpTckKxOd6moecyzF7G0fEwe2q/fqRRx4BAAwePBhO+nz11VcAEokPrBzz/OQ9kSq8RW13ixYtStoO66cDwKBBg0rteJyy4ZRTTgFQcM/TJl2rVq3I8CWdGbPUpSpsK9rEslErXI7PNT532DYW4NDZTc3OCIRTmvKZqdkSo0LQorZt2bj1GLWKIX9P5S9knSsgcV7LOrmJK2nHcRzHqaRUmjhpx3Gco4XHHnsMANC+fXsABTN+mZmZyMjIiM/+KZYHs86oqeK11rewvLyJZlSkzw19bNRvJyq1qW5DbdVW9E26bdbZBZ4b2tWt8puaDjTK8z1qhqAsKdOX9KxZswAkOmBWVhaOO+64kDMAkLhInObmJ+GJ0WlBPbncpk7bajIULqfTHTotzvWtaZWgs5Y1VcNPTm9r2zkNzpuLFV3YJs2iw7apM5s6lHEKlNvhzc/pcz03QDglKZkyZQoAYOjQoXDC8KFLUwIfTExvqw9MXnt+au1efuoUJ/s1ryX7jK4PAM8++yyAAodNwK+d41RFXEk7juOUMw0bNgSQGGTVqlUrSUlHpfZUVafe21Y6UCtGWLG8u1VJazQJB6D8TnUaHDASDjLVfq0K2mqTNZuQKi6acPDM/VCoaLpSy8YdbJOqbuYAYU6Q0qJMXtIsjciphZo1ayIrKwuZmZmoXbt2yAkGCKtFOpDt2rUr6XeqTHYMqkeqVE0+QuWtF5/whGsxC1XeqlK13cHj0eIeWhSEHYJox1ClTFTtE+0s3B5vFrZRw8/YnmCFFw194/mjIpw4cSKA9IpzHMnMnTsXQOLcMZ+vlhHl9+A5BsJOhLxG/L5161YAif7Pv7P/U0mrQ1pwP+yP7AdMWsPrz2vM6UoP53KcyocracdxnHKGgywt75ifn4/c3NyQOgYSAzU12Slq71WsegVEPZtV4KgXeaNGjQAkVDIFlpocgcRglr9RBKgXt3qiR2X+CrZNl9djt2YdNKujHntUvDSX0YxqFKWlTam+pKdNmwYAaNGiBYDkjpiVlYVYLIZjjz3WrGsKJC6sVieh+ty5cyeAxMnlyeMJ40XgxddkJERDXbg9nTZRZWlNtwT/b4XPsC1US6mmoawQC/1dE1xw+5r+U50z1P4e3IY6RzABSqqSn0cqTGHbuHFjAAUV3oDEtWYIIVGHFX0YaSiI+jOwv7Ffc30qc02cE/Uw0SlE9nFeS00vu2DBAgAJf5AbbrghxVlxHKescSXtOI5TTvzxj38EAHzve98DkBhI7969O/7/HTt2hDJyAWFzlTVYt5S05ojQbFyE7aAZhNtX9a85Kfg7B5InnHACgGSbttrNVdXrsdCMyU/NxkhUcatDsXqaW/Z7/q71EILnUkt7ag7z0s4tUaovadotOVJnp6JDRCwWQ3Z2dugiB9HgcVXK7DicUtHpIstrluh3nULixaHS4adVFrIw9CZQZW3ZmFNhJYHXWQcrEF/VW/A66MyBOocw0QpLi/br169Iba9q0Gv75JNPBhD2quY5tmZHNOkNr7Xl26Ce/poYR7fH39WTP7gPPjz0Iar3Dn/nbNa8efMAJJS1e4c7TvnjStpxHKec4KCKypDk5+fHB1rbt2+PD6SjlLSqSCtLlgoSDsLU+dYSRtu3bwcArFu3DkDCI51Z0jSHN9HQ0yCp8oXzWOg4uWbNmqS/n3XWWQDCjpjqcU2HSwoWtTmnsu+r3ZmfwX3oeebx8ryUFqXykmY8dOvWrQGEO1N+fn6SslAlx2WCn6r+1COcipon05o2sZwriP5dvW6thOs6XROFqiUuqzMJqUInrOksTWWnHSpVqjx1ughui9vWG47X9sQTTzSP+0jg0UcfBZCwPdPhJTg7BCQeAoX1A8AuR6ohIXqNtF9rH9BCDMFrybYydlujCvRBpWpd28xZBYbWHO0e/o5THriSdhzHKSc0ljiYhCkvLy/+qaYIIKG+tY581AANCDvTEiq9qEF6EA4M6RRJdUvzB0UZFXZRzHYqQDjw+/LLLwEA//73vwEA27ZtA5AQBZbqV0dfijieZ4ounhOtLqfniOeU5zo48xFU1cF1uY/SplRe0nQQsAo95OXlJZ3EqCmGoOoGoqdKgPDJ5AnTVHR6ESzUpqhFL9T131K3Ub9Zn5oMwCqjpk4LqnKtjqSew1Y4ArcfLFKiykyVNNvMG/2pp54CAPz85z/HkQBLRjZv3hxAQoVqtID6SqSasbHSIaqnPc+7Pjz4O68tl9d+GuzvdLLhvthWrsO/q28C/652dPUkpxPU8OHDCz12x3GKjytpx3GccoKOl3SuDQ6QGKZau3btSKdOK0RTTSXqQKoDQ8sMxk+ux8FZ06ZNASTSy27ZsgVAouIaHYbVNBQVaqsOvlTnX3zxBYBEHXsORjk4Zuhj1AxD1DFoQicVMBqOqmLOMiUG264mQSuuvKSUaGu0UbVp0wZAorFRMcWcysnJyQm5rgOJDqElJC2nBMKLwIumHtSWmtQOSW/xVB1blWZh6eJUmah9m8ok3XJs2mH4XcMUuD+eS0tJF1YOz1LQapsuqwD+8mby5MkAEsdDJxvtT/zUEqlEZ0X0vGmJSn2QWbHNvC6cFtTsdPwedCrSjHLcthVCom22UidaU4SO45Q+rqQdx3HKCapMmouiai3XrFkzsvoS0UGSigEO5DiY10GZZZrR8FNC1c9BH23SVNT//Oc/AQBr164FkPD+ZthicCDKdekxzrrrbDP3QRs0VTrt3mybDjyJFf9tWiQAACAASURBVPvNc6HK2aq1bVXRCu7bmpEo7YRPJXpJa5J4NaCz02RmZsaVdNDDNNhBrRhi7VA6DcRPzTWtLvp60oNJBICEy76luNUlXytYBdtgVeZSpwVVrTpNYk3BqBOETulYN6XuJ0pJa6fUPOKqxnieZ8+eDQDo378/qiKcTtPsbJb/gF5jnerSTHY6o5OqwADPPx/mGgvPzHuaqSz4QNSHuS5r+W9oKI/Oqug9x2xsPHc8xkGDBsFxnJLhStpxHKeMoTnlzDPPBJAYIFFUZGdnx23S9evXT+mICIQFB9WolT5Y1yOpcn1zcKapbflJe/KGDRsAAC+88AKAxICdswdAwmv75ZdfTtrHhRdeCCChvqmkuU+i5iJVzFY4oaZq1gQ/Viniwkw7en45cOe2p0+fDgC48cYbQ+sWhRK9pJl9SdWVehtnZWVF2qS1ehMQtqPqSbPst5o1KdXJ50WmgmZb9cSrR7U6ZQRd8zV2O5WS5TY0CN7y8ta2q4LmflX1alUutYMGO6DWwraUtNpAmSChqvHkk08CSBQJSJWlzfIP0BSJmoxCHVGsrHFcn57Xeq+oD4UWM4hCZ3RU7Wu/THXvqd8IQ3TUX6Gqz644TmXAlbTjOE45oYN7Doiys7Nx+PBhxGIx7Nu3L2T+A8IDRQ7WNYmROqgSy5PZSlGbytGQCpkD3Ndeey1pOxykBcMzVUETeoar2UkT6lj5yomVBEgH36qoVXxYNbuDy1hJpLRcckkp1kua1a46dOiQ9Lt6PmseYQbqR5Ufs+xfUWobCNvu1Bs2VTk2blerPlkKx4prDToUqJesqifCtvGGZYfhLIA11aU3rqp/q341OyJnPjSePfgg0BvZcrTQGFq2varAJPicXlOv7VSOJETPkyppwqlIXjPt72oH1n7MPqTb1+1EPUzUj4L9wSrAoOtb2dLoUMQZHfWx4H6YkXDAgAFwHKdouJJ2HMcpYzhwocOfFj/Jzs5Gbm4uYrEYdu7cGXIcBMLheGryUGdWDuys6leW4uZy6hypzo8URkz4Y1FYJq5LLrkEQEKNW5nA1NZshdaqScc6Vk1SZSWWUodhIFzEST9p9tTMZMWlWC9pqiY9EelUiIrFYoU6RagdV+1kOqpXJa2kq4R0GsVSkiQqFR/tiGqz05kF7SCqoqwpGp1m0tkDrVzF37W8m16voCe8ZcPXv1sew6Vdpq2sYDy0OpDozax9WqfL2P+shwMfUOo/wGulsfJU8vyd106nOK3+HjUrorNDGrevKp3nwjo2HhMfzgyrYXIL7X/czowZMwAA1113XWTbHccJ40racRynjKGtlQMXDpQ4wMnOzsahQ4eQmZmJPXv2hKo1AeHwVDVjUMFx2zpI17BULWmrKlVTz1qhnalsr4XFDdMEp8mnLBOLmntU+GjdaA441eFXnXwpPNVhk8vRwTi4DX6qxzlNW6WVeaxYW1G3eLXvqhoJ2s8OHz4c6fKvisVSsIruw1JA+qlTQDolo/GxXN6aQgq2WT161aFD28DzaVU80ptCbc9UNHou9KbTcxYV+K+d3grYJ6qk2ckrK6xu1aJFCwDhIvBqAyZ6zrgejzeYiCKIJpTQKmsaT63w/Fqxy+xrvPeCM0Xqbc1+pg939dzX7H9a61pt3dbLQ/0XuP8//elPAIDrr78+8pgdx0ngStpxHKeMoIdzy5YtASQGVRwEUnVx0JaXl4f9+/ebg0UgbE7Qwim6nA7quZwOFC2FrOa6VDWhlcLMm2qDtryqdeBoebir0ywVsgpJHaiqEud2aKYKKmmtYKYDZg5Gue+SxksX6yWtDghqjI+qHBV0W48qpF1YNZ+o79oxLHu45bLPzkEPVb0pVAXw4nH6I3hs+n9V35q2j/viRdQY11TJCHS5VNNN6SYzABLnXzuv5aCh3y1FWFngrIPly2BVHuO1U7u/lb2L66sNmvBGVj8Gy8mH14MxydoebWfwGPVBxSlWOhxZ2dU0Z4EqaZ25YQZCtaOrAxIfoIwSuemmm+A4TjSupB3HccoItblq+VAOkA4cOICsrCzk5+ejRo0akSGelne1DvCiQmCBcHpiyyyp6lYVuCbkCaZ6jqIw26xuSwf/mohJ19NZBHWG5QBVB4qaojlVhbHgMfCa6PnSRFD8XtIiRMV6SWuQu3qNkqCtK2gnVbUSXFdjLS1Xe5KqVrJlm+ZyVLE88ZwC0k5jeWhH2dW5DrfNG9NyblAP9qh628E2pToHih6LlVUNCNfV5ned6rJuvNIu01baqIOKVbWKx692W37XKUe9F3hO1ZmHy7FP8JPKWs+flTtew2K0LwX/ry8Eqzqb1b9032pX5zlkKkc+GLk9Tg9yFoDnhG1++OGHAQAjRoyI3L/jHM1U7ieq4zhOFUZVqSZh4iDt4MGDqFatGvLz83H88ceHEhEBiUGR9akmD923FSOsv6tAUjWrA8WgvTYKVcFBNKZYzUc6qFaRZJnjrBTJOqhVG7bawjkTEkSFnIYsEv69pImeSvSStqZVdPolKysLGRkZiMViqFmzZmiaBgirc1WVRBWxtVy6Xt5qe9bOoB04nRhvPSaNT1b1lqrkWSplnUopK3pugtdPlaA6Wmjn57WMqq9dmZg4cSKARBk9y8Ndb2ZVwPpw0LJ1Wu9Za3zzd97QWkKPSp3t0rSFmnWO+4+altNsZmwLla318NR+ZiVt4PbUi1y3azkB8QFdWfuM41QGXEk7juOUEWoa0GQ3QbMH/9agQYPQoA9IJIvRUsBqZ9UBpCZXsoQLsepOc/DFNnF76kxLrrjiCgCFF39hLu/WrVsDSAxKdTBqzQZYpYXVHEQs0yDPZbBgEpA4d8ECQjpwV6GiA+2SmgCLtbYVx6s1iIMer5mZmYjFYqhdu3ahSppY3tqpYoitsmuq1IuKFRoQ7MhWwfGoePGobVh29FSzClYblVTl2ICwl7I+CDQpQLp2zYqGN7/VPlWiqqT1U51lNHOYXgteQ65Pz2r13rbi/jllpseh6Qujrr0qal5LyxPferCp74lmP9PZCCu8Ru35fBDPnDkTADBw4MDIdjnO0YgracdxnFKG5hUmzaHK5eBMBzq1atVCXl5eyNzG5YGEkqYNWIuxaNIYDuw4CLIESiqzmaLKUENAWRu6adOmAJKV9KWXXgoAWLJkSdI6WhVMB346uE4VcquDVh30ql2Zy/N3HhPPZTA/uTqLsu28VuojQKZOnQoAGDJkCIpCsV7SaoxXFaJqjEoaKFAc6lmr/wdsb+xUv1u2Y+2AqXJTW3m3SWGF0vWCqy1Oj9W6GSw7unbcVL+nKiAfbE8qO6VVC1wdMiobvLH0Jrbyo/OccCqPdtzgQ/TyJzea+8v95iDOOvHYUEayVDXANQOZOhDplKOmNww+EDWMh1N2mgtA7fDWjI3Gjmvbicbv6+wCsSIxqkr+d8cpD1xJO47jlDJMksRPDeFTRZ2Tk4Pc3FxkZmZi7969oXA5IByXrAlxtOiLOh1qulbdLlEho/Zf7ofHRjHSo0ePpP0yPjiYKpelYS+66KKktrHNOrBLN7w0VfpiHfSq6Ubt+ZoIKMo0y4EyZzg4ANdBalRFs6JQooxj+qkqJDhCDhr0Vc0E17Hq5So66lfbqGLZay3bNtVCMF1fcPtRVbBUafAiRiUmiDo2rqdZpXhxrST0SqpZgsI6tCpia7ZEvZV57JrXvSJh/C2Q8OrmNdDrqV7W2re1H6h6jCIzMzN+7TQEhw9a3ti80fV8q1MQt0MFzfXYnuADUZWy5njXhA76YrCUszULofeYVVZQ28XtsO2s8uY4TjFe0qfc+WL8/wv6lCz+y3GqItcvPwDArpHrODQtcOChzrQamhesfLV79+74QCU4sKcaZOyumg1V/aknuYU1qNLBmA6y2J5mzZoltUvrWwcVJAfxJ510UtI6/M62q1iynA9TmS8J28rrwnOlYYXbt29P2l5UjDMH7BwgB/0GgutoqdeomOt0KNJL+o9//COAVvHveoJ0OiV44oIG+yjXdKvYN5exOoqqQ7aBN4GlYq2YYnXh1+pYVhxx8P+qejQPs6IxtZ9//jkA4KuvvgKQKHPXtm1bAOGOrFj2ez2HUc4XasfmA0BtopqOT7NPVQaCqlJrG6vTB79rf1HVqVn1CiMWi4UeqFTO1gOR7dH98SGidae1FnRwpkhj13nt+AJQnwkr25kVI87fdeZMbdx6LDodqJ7vfKh7tSzHKYaSvu2kz9GxY8e0lu3x5w3x/+d+U3h+V8epjFy1eGex131/4250mvg+/n7r90qxRU5VgAMNVdCEvwdtobFYLF7OlwOWoC2UA00Ovjjg0wGcFvAhlupUE481qFc7rnpia8ldKsygCZLrNmnSBAAiZwyCx6gDQB0IahlV61g1DTSPlZ7yaovWcOLg4JfnXbOkcR1NMMRjKW7msSK9pHXUr8pRl4vi0KFDkcHdqVShtbzaFNetWwcgoUgYBpBKzZJUqkGD3oPxsrx47Kxq47OUre6bnXvx4sVJf7/lllsAJBS1hj7oLAPRm03DP6KcM9S+aNkbNQ80v0+ePBkAMGzYsNC2y4vgjW/VBadtl6qS1640q3llZ2eHKkNp/9Wyg/Rn0GT+fDhZhQmCfYr9iNukTdp6iOg2La/sVEkwtI1W6kd9ganjU3EdbRznSMK9ux3HcUqJKVOmAEB8tjGqiA0QnY44OBjSlKtAeDDNgSXFgVWOV9MSK5bDoOUsqwLGCsnl3zdt2hTfF8uZavIfPT8q5KxEUOrISlR5c3k9Z2qK0dzqUSGlPC71RNciN2rupABgHxk6dGho21EU6SWtnqkDX94X/9uMrseGMg1F0evpr0O/Lbvu1NBvqTyYtUOxI6xatQpA2FM0aJ+M2r6GG6iCZifizRNVJ5gXXhPcW7Zg65it2YP//Oc/ABJOFuzw6j1rxUdr/LZmuQquY2VYS7WcqqKKJCpsQo9Dvad1+aUDCnww/nvW58VuR25urnk+eE20lKE6p7BPsV06k8OHCvsgEE7cQLTmNZWubottsTKKWbW41enHyvWtSlqnDflQ87hp52jGlbTjOE4pwcGWlcfZyj2dl5cX///evXsjB5UcqKlyVuVqOTlaBYB00KafVnhrqjBZdfoFwlWu1IRKrFrYlm0/VU1rHZDSFMSBI23TaqpU58dg26xCNGpPJzz2onp5F+klzY0nYmETrufVq1dPS0lHEcxIpp0zVWwvO+ynn34KAFi2bBmARLA8lQgD7y1FaKldtR1qcvnghaJi4QWlIlCbn9W5UwXu89hOPbVg5oHhBGpv17rVun1N+h9sh2WLTuUzoDmnC/NLKC+CbVB1p7ZgwuPgDacKuzgcOHAg5F2uNmbNeKbXiP1OZ3R0yjH4cNMXAvets0Cq1tUbXJW0frdQG7feB2rD1mQZPGel6R/gOFWNIivpu/5ZB/jn7tDvJfGCveiRfwEA3hh6ZrG34TiVlQunfhz///u/6FSBLXHKCx1UWc53QSfO3Nxc5OfnY8eOHZEmKzV50BFQbaMqlnTfqkpVoJBUNZstgaN1qYMOnCpQrJBajdHWwayKAi0IpKhA0QEkz4XOHlAEBk1GOtC3ChapOUcHoelSpJd0SYtXp+KYY44JdTS9mPr71q1bAQDPPPNM0rZ40fWiqJODpWatjGa8+aIymKnnN216bINOfakSUfXWuXNnAMCKFSuSltu4sSBndMuWLZOOSY/BSkKvN1PUdJVViUtRpc3ObGV+qyhUxfHG4TXS68xrkKqgfVHh+eGNqg9a9WvQ9ljt09SEQPie0SxyRGP+rRAb9ZC37qF0H/ok1exNZZiVcZyKotCXdO/evZO+b9u2Dbm7Cg9hssjfURAznbvsd+YyN/yrtjm9ajlD8YW4fv36pN83b94MIDxlZoWV6H4sxzF1top6wWmIiTW1aTl2aV5fPoCJDkysJBS6fT2WwqbXU+XCJdY++LLgFH1F8M0338T/bxWBsIo86DlsExg8fvj1PqTC6u/fjzgd5zQ7zqz1a9X+tWybUf1aPYm1v1iqRZWY2iZVJVrKQZfXaXNLRalT33PPPRc6trLk5z//Ofr06VOkddSLWM+tFTLH1Mn5+fnIzc2Nrxe0tXKbDRo0ABD2GtbY3lRKmm211L410LbKqVqD/qADp15j9aq2wlQtO28qM5t+V9OJvh80HzfPf/CYdFBrqXq9Hqkchy2KESddePq1kpAqhjmIJptXUnmHF7Ut6Zxg68Gqf7fiSy0FTNuzvmR1+ivVsaYT75oqlluXszK3Ffe8lybBNlh5o6364OnMNpQW+fn5pi+G9dK2KkpFqU7tN+leI+sap/JTsPpQqnNoDUQIM/AxusFxjgYKfUnrFPLLL7+MwX/JNZYuHCqKrK53mMvM/EUncypYpyk5un7rrbcAADNmzEhannVLW/xfPVd+algJURuRTidSlX355ZeR6wXbyJcqs+poSk0r2Qnzxr777rtJ+9YiD3wQn3/++QAS096anEJfPDy2HTt2AEh4OgbDyHRaOJXC09Ej19uwoUBJDho0CBXFggUL4v/XbEgMzdNrobYojqx5buvVq4fLnkjuA1Gk09/JoyPPDTl+abpV9kc1nWiikKD9j9ddU5/ymNmvVDExiYhmg2J/5OfOnTuT9kOFR7MY28Q2snAMp/i5fQ0r4zlnP9U6ypU5FIvOtfxUs5uatIIqlvWkMzMzQ9mrALtesRW3nK6zp9ZeJtYgjeh+rWdF1DYt04eiToXqsGuZZCwxpXm1eS65PV4fPu91lijYJmtgb9nm1fs+XYqkpGOxGKZfXC1+M7IDZWdn45LH/lOkHVvbV1Vo2a94w3Jam9COqzGX1vYs9aoPSe6P36PCC3jy1a6o091Wp+dLnA9JPpzYNvXO5YtQk9OrGrP2pw9+wK6fbSlOjXHlMXCgUhFMnDgRQCLbHJA4Lp5DbbdOg1kl88pihuAHE1YCAP46rEPo4czzyYeHRheofTjoN8Jt8fqyP/FlqVEA7HeW846Gz+ig0MpEqL4YVviNlQeekRlRmQod50jHe73jOE4JmTRpEgCgXbt2ABIDHiuxCwciHHAHvYdjsVh8xi84o8jfONjiTIYKC+6DMxapMo5Z5XqtT53NtBS1zhYE96V+B5pJzEo2ZYXGWjOwRGcJrCxqaj7V0Mngvoj6algzEurzkS5FVtJAOM4yMzMTz/Q+IaReu85YU6TGBA/emlbl9y1btgBITKUpnEpj4XH1PiZW0nj15qaSfuWVV5LWv/DCC+P/f/311wEAF198cdK6vLlU1WtHpArijc7f6c2tDmmcRdi2bRuA8JSuptqzMj0Fp4YsxwvLw1adVCpDXWl9aAV/03aqXV/rMus0flmquQsmf1jo3z/4n/PNOH/ek8GZHXVgUYctPuQ1E5gur2kmdTv6kNLyf5xZUiWuy2tNb04ZqwNZUdMqOk5VxpW04zhOCbF8HFThacUpLr9//35kZWUhPz8fWVlZ8fWjErnwNw6yqKy1nrSmmtVBuYbopUr9m8o5VBWjmv2AxAyDmjpUpavaTDc9tLUeB3oUTio02VZ+8nryXEblUOe2OHjUQbEOSlVJp1uEyAMQHcdxHKeSUiQlbRVz5wjByhmbLhkZGWboCffNqbGvvy4o1KFTuhy90GmII1VNs6hezJbHJEerWjaScIo7yGuvvQYA6NmzJ4DwNLd6PeoUKkdxZ555ZlJbOL2tiTdYnrNRo0YAwg5kajLgyFqPPbgvkirGUhNcRFXvKW9GjBgBAHj88cfjv2mfVa9Ongs1KWgazooMLdu/f39cHannddR0tx6Lmjm4Lr2G+ck+r/1DVZ32W55bjZXX5D6E/VdjYK3CNMVNO+w4VRmf7nYcxykhnMbV6VorOY4OnLKzs5GVlYW8vDzUqFEj0glJIyj4Xb32uQ/17lc/EcvfRAfqVoInK9EOB18c9AWFFKM+Ug3iU+Vn0IgZRSNpNAe+DiAp/jj45fWM2r4mquH1UDOD+nBo27ivVBQjmUmic+jF0b+/ekMbAMCPH/13odv96O7/Cv1mhQ/x5LJD8qbo2rUrAKBVq4LSgizjyA7CDsOLoSkQdVSvpQKLw8KFCwEAP/vZz5L2zW2zI6ijDr/TO/Pss88GAKxduxZAIp6aHY6OY+pAxu3pzAavDztJMARLHeiIpaCttKlFzapTFgTboPly2R80DElr8JKgc9PsS2vFrxEdE3mDMpwqyAvXnBwq16ixy/pQv/G18GzUtm3b4vth3+F3XvOoYilK8MUQPHZ+8t5inDKX5z51loH75LVnn+B29IHLvqIzS/owszxwK0OiHMcpL1xJO47jlBAOeFKF3yhBZR2LxZCRkYEaNWqEar5zGSCs1HTQxUgUVYsa/aFYyVFUyKii1uU16Urw2Pk3HoOqck3EZH0SjZDR5dRhjOeEUUFcj7H4HIgGU5kGjz24LQ4qacZR85JVNU4TKqWiREpaE/Zb5SZTcfjwYbN8o3VRqJj1gBm2wRPDjqLeezzpPKHswBp/pyFbxYEd09q2egVqh6FaatOmYGZCQ1s0OQU/tYQl4Tnj/oKqMVVpUPU/0DhPjdWsSILTdlbOX0XLNWpyEb3h0kmxWbdu3XgfsB7WmpUI+C60zLfffhvKu27lJAfsBxz3pVN66nmsD3mdVdKynvow53rsv9onNHe3blf9AbSPOc7RQLko6beGn43Bn9RCLBbDY7edF3ooOM6RAktRDvy0YLD4v8M6FGs7sy8teLGpKSkd/mvKR2kt985t56W9zaFv5APQgUjB/fvbc9LezBGLlU3QquQV5ZDJf9nZ2ZG1CXRbHDhaaUBTZQ9Uc4NV4EQHtJaNmm2lmSTqWDnQ0+pz2pZUNmdLOet558BO0+dSzHEAqYNwonUSgISStkKt1CSopjSSbi6JIr2krWo5mp9Wp0VisRgOHTqEzMxM7NmzJ+R9HEQ9ktU5gSeTubh1SkVH2TyR7BSaTlHtbJwq4kXlCe7WrRuAgvzlqbjkkkuStqmqnW1im61Uo2wDlYRm+OH5p/1dS4laatG6MaLWVW9onYbjp2ZUonKsSIIPOD1negOpDVrjNtUzWfsn+416JuuDVFWqPpTUJ0Jt2DobYxXmSAfuS3N5a5algvO41dwOH1o6W6B+JPSl4P50VkJjXnW2RmNbHedowG3SjnOUcv4fPwj99tbws4u8nfv+3fD//tcQtzf5ooStqproYEkHIKkKSWRkZMQHKfT0Dq4PhAuW8G8cVKmpj4MZ/m5Vy7NUv37qgFIHjAzd42AtSoRx4MxyuxwIpioxqepdz7OVcEVTk3KAqMln1HFTB4g0KQKJAb0OwFX9qwDQTHvphhIW6SWtNif15FUlzRN1zDHHIC8vLzTNHZUGVEfhar/lydb8tWovV1uzXhRVmexYnKphZ2LxCrare/fuAKLtaXpTcF9sKy+4uvnr9BW3w5uRSkQ7nDouWErd6sAk+F3PN9uQqvKOVQGmIrnpppvi/58/fz6AxDnWB5IqXXXSUS92tZNqNR5NZakFM9gHdEpRK5PpA1LtsvpyKOl5z8vLMyoF2Uo6SO3ateP91UrbOX36dADhc58qXzOP+ZZbbkn/gByniuNK2nEcp4ToYEkHHul4eatDZ3A7QNjUxMGQJqfhAFCVneU8aTnrWrMAOnjigJUhoEzmVBgXXXQRgESYKeOn9bxZFQStght6vnlOND6bA1HNP6+OxjRjBc+dmod031yWA3Id6HPgnm4oYZFe0lSbVAWWS77a86pXr46MjAxkZGTg2GOPDdU3BhIdjweiB6bKzarLqlmJeNJ5UfjJfVOJ017G7/SwZifidqiSC7P98SLqdAineLgPjVfmcpoBijMWREMwdNaB50I7vBX7HOwsOiui+2SbVTnr7ElliJMOonWhiR6HEqz3G/xOdMpKH4R8gOny2s+J7idVZZ1UPPrj7Pg2ua9rX9xlLp+qyEcqGjVqhJ/85CeFLnPjjTdG/s4So+qQQ4pah9dxjgRcSTuO45QQnapX80k63vlBlRc12E125LPLNKqqT+XNrVipmdWLm3Zaig+qz06dOkWuH2wDRc8333wDIOwhnarkpB6TVW5TY5gVDvzYdo2nZnuCKXHVi1tt0FayoOJSLCVNNak5qdkolft16tRBVlYWMjIykJ2dHT8xnEoAwieF6lIruPCCUxGzo3J9jdlUj2h1KNDsSrQTa5o/ndLQTqH/D54HrsPsVKqouRzbpJVatKNYdlNLYWsHVI/lYAfW86MqXjuixgvzUxVkRTNgwAAAwJ///Oek363EDFa+aCucgudL+xuvtd4T6j+g5UzZVzQ3t84cRU1FzrsiOR0k7wley8cuKWjDoFdKP5Y9qmpTutx6661J36ms2Y+HDx9e/IY5ThXFlbTjOE4JCSrenz6+PnKZd2//fui3qEF+ZmZmpL1SQ+S0cI6KgihzFpAY2KkYUJOj1l7XNMR0cmS7TjvtNADhELvgQFK9pblNmhs1URZRM1OqHBs6u6AhjhQuFJ48Foo77p+D5KAa1vOqmd/0U/ddlDBJoJje3argrIxFwUDxjIwM5OXlIScnJ36BgkqaF5IdUL1i1eVeHQ14slVxaLwpYYekmlUHA1VKPMFqKw9OY1mezbot2rnZQXk+eN60o6hnsPWpMxs8BzxG7odtpt09mJ5OlTHPG9ukN7DeyNyn5dlb0Wj4iipSy6vacg5hP+J2+V3Ph3pna5527l9nNzSTnk7P8fqoz0awzTrFWJa23VgsFvek79OnT4m2pcracY5GXEk7zlHK5AuSkxJVq1atxFPgA5buBVADvz1nX0mbV6XIzc3FDa/mAIhW0QBw3u/fAVBQeCgqPJIDxF27dkWW5VSTS1ELjXDgyME6y/3yO/dJJaxFW7g+TTEnn3wygPDAsrB2abwzhRrNQjpo1mIumgKZAkZTKKtZ0iryooVu1OwUlfnMDF7fuQAAIABJREFUMoFZyXn0/FltsijSS/qGG24AALz66qtJjVAlx89gxq38/HwcPnwYe/bsCamO4LKqYNSbW5WOlgez8lozVpm/qzc4VaWecKoO7k+D36NQ9am2evVe1SkX7XjsOBoPbcWda0pBtWMS7jd4LNwn28BPXgeNhdcprZJUDSsPNJ80v+sDxjqnVvy+xklrf1GbvX6qR78mQuB10IeI9vdgrnlr2lNnP5LDT0qnwlRR0pg6jmPjStpxHKeEFAyW0hvgfPfdd5GFWTiI3r59eyiMFUgMzHQgbdWu1kEaB6SfffYZAGDTpk1J++AnB4o0tVA5N2nSJGl/xSkZqselJjeNBdeUs2re4TnjOWnevHnSdq1MZnouec4sVRyc0bCqWimqnIub8KlYL2mO2mknowrROtPB7Eh5eXnIy8vDvn37IoO51U1dPZ41axdPmtVRUyWFp5JSdaHt0GkV9V6NCi/gRdMKWtwHj0nPG1W/2jP5qXZU9TDXnMtaG9qqNsTrGTz+oGd+EC6r3vK8LnQAqaxY8d16rfRm1YxjOuWos0capUB4vjRxvz6c9AHLPqLRB9ovg7NTWpZPZz10lioWi2H6xeH+VqNGDfRfsgdFoTJknCtvJl+Qj2bNmuHyJzcWulyvp7+O/D13q8eBO2FcSTuO45QQiogTTjihxNvavXt3ZPUsLcqiOaGtkrQckDIm+T//+Q+AxEBRTVrqHEoFrWGwxDIJRWGlJOY6NA+p0GCIrDq/ss2bN29O+k6zkDoe81yo17a2o7Cyy6nyiquI0oF5uShpqiV6QmslKr1IOTk5cSV98ODBeGOj3NrVUK9qgKN+/p371ixmGiNKJUJlqEpKlZPuhxdVvZ+DWAnpeTH0ptJ84xqjrPmydRZBO5SGZBArL7omAgDCSplt0yT9WjeabRo8eHDovFQm1Duf7ddzxOO2Evrr9qxShXqtrHKCOq1mOZdof9dqW8GZHraNDz71CNf71ipZWJTa4NM6Z2Lnzp04CoW045QJrqQdx3FKCMMcg2GlxWXgywWDqcd/WifS3qke0NYgnHBA/dVXXwFIDNbUdKdOshQ0wVziwe1bCZ4Kw1LdegzcJsUBv2tCJu6T+cJZUviUU05JOkarzVZaXj2HUUWIVBGrklalXK42aWZvWrx4ccFGRF3oCaRNOj8/P8nDOCqQnx1C7bJWzKc11aAx22q75vrqKa0XR2OPC6sgpReB+1aPdD1fan/UKRq2lR1U7ZTsaKqsdXYgSjkHlwfClc60mpiqMnXwqOxY8dtz584FkDhu9aK20HPK9YPhNEGsHAPqwGLd0FpvWmuRBx867B9aeIH9i22zHkiafjIdcnJyPMe245QirqQdx0nJ1AszkJmZiRtfO5R64aMQls984YUXUFphbNe8kKhhvPDnJ4XS+aYbL61Vqoia7vidKlRt0Yqa53SAT6IchC0BoXCA2bRpUwAJ2zPhoPeHP/whgLA3eIMGDQCEhRCxPLNV/AWPyUohnEop69/THcyW6CXdo0cPAMBLL70EIBwbSqxSaMGLrwnQebI1F7Lm/ebvvJhagkyrGFll29SBgIpE42Gt9YPraq1fTbKv27YUtnp/82ajfd46R+ococde2PSYdkrdtoZzsA29evUKnY+qBNUlj4fXm32ZSlfPj2aT02vAa6XX2KoHrbkGeL6pajUuWh+IUaEi7G+cOuRDnspaH5TaX5Lj6wt/Se/Zs6fS5W13nKqMK2nHcdJmRtdj4wMQTbmrA4ajkYKBTemXaa1Xr15o0GUlVFLTnKrLt956CwDQuXPnpO21aNECQCIHt2bd4nXnoE2deTUPd5Q65rY46KXS1eJBqnxPPPFEAAmVv3bt2qTtqsMvM5g1a9Ys6e+KlddcE0EFVXEqBW2V69UytOmaCEvlJf3FF18ACKthtU3r/4NTDTxwutirclCPVB4w1QEvuqaR031peAFPqHrh8uJwv2oH1pzMweMNZloLovmuuQ31ylbVxb9zPSohfVhydoEqS8MMrJCA4LSLxu9qx+Inb/j16wvSINJho6qinvda9F3rietDRR+YaqO2HGasG1tndNQ3QHMR6HUL/o3Hog8czjrxmDSZhap1nb7jMfMe4Xo33XQTnNKl87RPAAB/HdahglvilDeupB3HcYpJ7969k75v27YNubuKZ5PO37EBAJC77HfmMkP/VTulB7L+zkEUB9SE9l0N30uVUUzttaqwtchMVLiqlttNVcmLcJ8cEHIAqemIqdDVpKPHYCV40pjmKApzIg6SKonU/Pnz8cwzz5j7KZWX9LBhwwAAU6dOBQC0bNkSQLKL/+HDh5Gfn5+k2qKmEIJZyoBwujbGZqs6VPVpXWRi1b7WwuaaSYpTPDzBUY4H7BjcpnYI7aDaqbXDWuEJar/U6UfNYKZx2eoxHGyjdlK1serDqapz2WWXAQAWLFgAINF32Z84NcVpNPWB0ClIneXQ2HmdNdG/a/lAvdY69anTc8Hf1COc19CKsdai9TqdrTNEbINPd5c972/cHf//uc3rFLKkc6TgStpxHKeYqAJ6/PHHMfrjesXaFhV0Vtc70lp+xi86FSwv6V51wLd69WoAwBtvvJH09zPOOAMA0Lp1awBhQaJCRxW0mgQ5gF23bh2AhG06KMY48KNtmbZmmvDUvGOZkbivlStXAgBefPHFqFOEsWPHAgAaN26cdExskwohHYxTjAUVtZqmVOxYibn4O7f505/+NLLNSqm+pIcMGQIAmDRpEoBEirxq1aolpXeLclFXexenK9TWbGVsUrWpy2smKc08pXZadnwuzxP7yScFtqGlS5ea54Fx5OyI2rkJt622bE1SoLHfmuBdy6zRbsmbSFPqaZL7YCywVUebx3/ttdeax30kQJWpObh5TtXOyweUVqvSqmC8Rjpzo/3R8vCnKqZy1wIBqnqD+9Ic8HwQ6cNF7ew8ZvXo19CboiSzcBynaLiSdhzHKSUKYpGLp6SLSqr4aBUszCBG8USzpKY8VpOfDuo1uZKGpWpyJSp4APjv//7vpGUtE6CaLfWYGzZsCABo06YNAFtJa+EcVbkq4qww2XS8u1Ml1LLalIoyeUkzsJ9MnToVOTk5iMVi2LFjR7zRVAVAwo6rgfWqMnX6Qz2ptVi4Bv5r59ATrDGi/H3NmjUAClfQZNasWQCA4cOHA0iUerM8eq3aw1aMsraRWInceb5Zvo3LUQ0GnS50yobnlX4HRzrquGL5AfA7rxHPpcau9+nTp1jt4GwUfTA4Lcjt8t5h/9YHQRDLsUdnoTTHvBXio7NWbBP7u+M4pYcracdxnFJixIgRePjOaFVXGrxz23kpkxJZlaVofmvUqBGAcPphK4OYChZ6VtM0pLZsbpemGCa9Cm5LlbI6z6pjr9rJCTORDRw4EEDCHKmFanTWwfKMZ7s0b3lQ/Wo4qqp/NQOpGVPDg1NRLi/pIUOG4NVXX8WWLVtQo0aNuCoI2kJVAWucKjua2s009phw21phyrIPE714TKVXmIu8xeeffw4goYas7FRWej4eM49B60NrInetTMWbiMvzpmFnohrjckBCVfNzxIgRRT7uqsj06dMBJGY9VJnqjIsVpsJzz5SKxUVnowi9zzllyYeH3hdAeJbJqmBG9N7Q7GhcXvOG86HtOE7p40racRynisCSv0EsU6AWLaIdl4Nz/l29tdXMptuNCvcDwgVdKDKCzpBqUlGnRA0X1dhjdazksXOWgF7ZWkAp3ZhvVdRaqhgIJw7Scr0aimgVSkqXcn1JZ2VloW7duiGPbSCsCjndoKN2/l2Vg8Z6qqu9FWhu5RUnLO+mMMNWcPklS5YkLcMQBCoN3iRRBd2j2qQ1iTVTmWZl0/R8XF+VPGcpeE6DHZBe9UV1bqjq6LSY2mV1tkJvRPUHKCt+9rOfAQBmzJgBIDHdx6xhQT8PnU2y7gUei4abaO5uVdKc5erfv38pHJlTGG+POKeim+BUEK6kHcdxKjl79+41c0RTmHCgRnSQTydDNbtZ6WN14Ko2ay0XTCg+oga0WgCJ29JyvmpTVvuvHgOFX/PmzQGEzaSWkubAlUpewwyjih+p4646k1pmSf5+/fXXR7bFolxf0rFYDNnZ2fETEZwu0WB2zS/M39X7mFgu+/rdUtRWh2Q7LrnkEgCJzsFOEdw+s1ZpbCyPVz2HraxolkOIdmCtvc1P3qw61aN5z3kMQbgPLQh/pMN80/PmzQMQzhanqHJWf4Ky5rrrrgMATJ48GQDQqlWrpHYAif6gCSL0Ic/+yYe05SvB/sN4ey0b6BSw7qHu8f+fkqYT2V8Gn4FhH9VALBbD5JvamUUanKMPV9KO4zgVzJ49e5K8p3XgDaTOma02YFV8VoInLX+qRVvUaZJ/p8iiwFLFqKltgYQZU9Pmphrkqte3ChGuR3MdxR3PiSbc0fTHPHbNCU6CStqqha0x1TrzoPnF06VCXtJRuX61cpM6FmhubQ2c13rQqZLQW8uzTVQLun0rz3ZUW9lG2oi5Lju1FbNtVUQiXE9jZNXhQx1C1CNebdtAeNrpaIOe7jy3OpWljivqPDNo0CAA6cXTlwaMX582bRqA6MpyVrIK7VecItRpVH7nA4zOObfeemtpHorjOBG4knYcxykjOPX9wgsvYNhf7QxhlpK2zC2AreC0DK+VyEm3rdtR8xzRQZ/WYC6s7aqgtZAS/66DX1X9qqA1IVRUSGLw2PW7FpXRdgQFpc5Q6PlUEwWP1XJCTkW5vqTz8vKwf//+eHax4AlU5auxmqqcLe9btedan7pfnshNmzYBSKgKKioup/Wug968/L96nPNCM/E848Q5NWMF6lul0KxUeprAXUMqNAhfnSOC/7fiyI90qITnzJmT9HvQAx4IP9D69etXDq2zoU390Ucfjf/GWG2NkyY61cc+T3u34zgVjytpx3GcMubLL78EcHLo95ndCpRf0Ht7//79kSpWBYsWAbJMfJaDquaWVlOjqlAraRI/tc008zEUFQjbfLVyl4ozbTPbSOFhhauq3dhCY8mJmjeDyp5tt0yD6iuwc+dOAMVPm1uuL+nDhw9j9+7dcaN+sBOoQV9rH+sUC1FPaXVysJQz4Qn/+uuvASSmJOj5rLZwtkftzMG2WUXUqcL52bZtWwAJRZ2q2LrlFR6liIPngG3UusAkuL8oZ4+jESrjxx57DEDCUUXLAla2GYcbbrgh9Nvjjz8OIDwrxGtMBR21ruM4FYsracdxnDJm6NCh+E1EOFZQAcZiMcRiMRxzzDEh1QyE1aEKGA1DVSdbzYjFfVvhlmpKpJpUJ1srNJRiLOjVrG2nCLLimVWtp7LVW7nBiZpT9Xc1HWp7g23RxEXqMc7lrrzyysi2pku5vqTz8/PjahpIZMICwgpa6/PyJPEEUMny71ofWj2m9eJyO1QR7EiqVhUt5B1UnGyj1onWddm5OQ3EtmsmtlQZybQDE66vub+ppLi8Tl8BCe/mwYMHR+77aIM2akIvap7TqmC/veaaayq6CQ4KnMg4q8EiF46TClfSjuM45QQHTMuWLQOQrIKDhXVU/Qb/r57I6umstmn18rbK+mqInsY7a4pmDvY1xlmdfINKnU7DPG7+zTJnsg2aiElDdNVcapkO9ZyofVnPDQkKJq3Qpd7c3BbNmiWlXF/SJ554Ivr37x8fTdIWC4Rjh/ViUX1bNW/VkUAdCIhO3fC72usUS9VG2ST1ImoRcA20V9uxzipo2IKi1a40O5uWTuPfo5wiNPTBSYZe1I7jOOWBK2nHcZxyhh6/NHFlZ2cjFoshIyMDNWrUiLSFWopZK0apxzMH72rTthIzUU3SpEhTIJdXGzS3pyZHdaQNrqNqXW3wam+niFKhQdQr3PLq1qxtVgUwLfhUWMy3toGC8uc//3lkG4rK0ZlWynEcx3GqABWipGmXeeaZZ+K/McEHR04cydAeooUfONJSm4WmQtRYP6JTx8wpyxGYVkPR0auOvIBwmBa3xd/ZJu6Lo0qdliYaKsM281xoykompeD2+HeeE06La+lCjuqDx+s4TtlBlfXaa68BSK7xXK1atUjnUXX0VFMb1+G2NG5ZTYEa76zb5/OEaBUtte9apsKgXZ1tUjs328DnO02BGrut8dNacMnKrmZlWbParEo/6MltJcjisRTExJcerqQdx3Ecp5JSoTbp4Eht27ZtAIDjjz8+aRlV0FYyCVWZOrrhiEjj8XQ7lk1DR1zqqBZcR5OYRFW0idon285Rpo7+NC5Pnex0lKqxf1o6k96HtD0BwMCBA+E4TvkQTJ6Ul5eHjIwM5OXlhbJWBf9vpflVD2m142rctM5G0paqjq7cjnpEa1ireqQXprDVWdbKZsZZAW7LKiGrnu7pJn7SutSWsg6i6p1tZeGZ0g5fdSXtOI7jOJWUClXSwaIEDMviiIrxdPR+VGXMEZAqYQ2V0vg3zZKjIVwab0eVq+UdOXILjuiYOpKjP8I268hTyx9qqJWmJFUlbWX80RKU3L7WNaWN2tWz41QM9M9Zvnw5cnNzkZmZiYMHD4Y8s4Fw8iL1PFZ1avmu6PPCKtRj+Qfx+aHPYg1j1QpXQSybsebw5rHo81a9vNUPyMpjrrW4Nf20PpOjQoJ1BoNt5GxwaeNK2nEcx3EqKZUmTpojyqeffhpAYrTG0YrapKNGmsHfuZ7agzm604QjausmamfmSItKP2hXZ3IW9aq2lLMei46QrRg9HTlrLl4eE0e6TH2qRUFcQTtO5WDNmjXx6le8XxV9PmhJXK3kxGcTl1dfGc0MxhlAK4aYPizcrj6/dHmtjhVsszWTqcupgtYsaZbXt/omcT22XWcZNeZb2xe0XbOtvE47duwAAIwYMQJlgStpx3Ecx6mkVBolTfr27QsAeOGFFwCER4+W97XGrHF5VdyWh7TmtbXqlHJk17BhQwAJuy6QiH9Wr24d+epojftUuzjbqkpabdc6W8D1NAaRxxL0BXAcp+IZOnQoXnvtNWzZsgVffPEFWrduDSA6W5dm4dLsW5rDQb2xqZh1Zo/b5fJ8tnE73P8JJ5yQtB6xsnmxncHj4bOI0Tbcl6p/baNmIrOysBGNaOEMLY+J22E71BatnupA4rnKtvbp0wdliStpx3Ecx6mkVDolTX76058W+vdHHnkEQHiUyJEPR2xUtxzNqd2GqMcjR0mWPbl+/fpJ+wHC9haNb9bKLTri1HyxbIN+quLWtnlpQsdxnCODSvuSdhzHOdo44YQTcM0112DhwoUAgObNm8f/ptPbmqRIzWIcvGtJSQ7mOcWrg3xNHML9qrlOE4dYiaWCYbHqkKWmPw2h1Sl+dSizprnViY1Os1puk+h21HzKcxjcVmmn/7Sosi/p4mZ1mTdvHoBErnDNQEbY8a182uwkwfWsGG3aQ9jpNVaPswH0Fty+fTsAoH///kU9PMdxHOcIosq+pB3HcY5UevbsCQBYtmxZ/DeqSi3Mo4k71HGVUCxoWlB1vtJCQVSRVJAUGRQdXE9Ng1GJQFSh6u/qGKzJSiznWzUtso2qoPXvPDZ1DOa55HoMPwOADRs2AAAGDBiA8qDcXtLBilcVyZVXXgkAmDVrFgCgWbNmAMIdn3lsOc1i1WoNqmdV49wWFbJWr2Hn5d9Lq/6o4ziOc2RQbi/ppUuXYvjw4Vi9enX8txdffBGTJk1CTk4OmjRpgokTJ8YdshwnyP33348XXnghnkSmVatWmDZtGnJycjB69Gi88847AICLLroIo0ePThq9O05VpWvXrvH/v/LKKwAS4Z9aOtJK4KHlfLmcrm8paYZscTtUpzQFapIUtS8H70UNcVLTnzr+aqiVKmiuz1kCTQClxUis5Fb6qeV9v/766/gx9O7dG+VJubykv/jiC4wdOzZpiuPDDz/E6NGjsXjxYjRr1gz33nsvxo8fj/Hjx5dHk0JTFXPmzAEQvjj16tUDYDtnBFGnA63LSjv65MmTAQDDhg0r2UEcRaxcuRJTpkzBeeedl/T7zJkzsWPHDrz22mvIy8tDz5498fzzz+OKK66ooJY6juOUHoW+pEeNGoUGDRrgzjvvBAAsWLAAY8aMQZMmTULL3n333ejcuXPo9/3792P48OG49957cfPNN8d/X7BgAa688sr4dPPtt9+OnTt3luRYnErK008/jYcffhjLly9HLBbDpZdeimHDhmH69OmhZaP60cGDB/HJJ59g6tSpuPPOO9GyZUvcd999aNKkCQYPHozrrrsOGRkZ2L59O3bt2hVX2xWJmy6c0uaSSy4BUFCMA0gkFWERIi1soWV6VURYoZzq6awOsZpCWW3h1vLBZTSNJxWzpmFWL25um8qYYorqXo9FFbOVJIXfNRHMpk2bAAA9evRARVHoS7p///7o168f7rjjDmRlZeGJJ57AlClTIl/GFr/4xS9wzTXXoG3btkm/r127Fm3btsXAgQOxYcMGtG3bFvfdd19xjqFUsLJwzZgxA0C4OhbRDhmEHWnkyJFJvx9tCrpv37544403MG7cOOTk5OAHP/gB+vTpk3amni1btuCHP/whfvGLX+D000/HtGnTMHDgQLz88suIxWKoVq0aHnzwQcycORMdOnTAD37wgzI+otSUdRYix3GODgp9Sbdv3x7NmjXDq6++ipYtW2LLli3IyMhAly5dQstGKaBZs2YhKysLV155Zdwjjhw6dAjLly/HU089heOPPx7jxo3DqFGj4i9F58jioYceQpcuXVC9enUsWbIEb7zxBsaOHRtaLqofNW/eHHPnzo1/v+mmm/Dwww9jw4YN8TjSu+66C6NGjcKoUaNw5513YuLEiWV6PI5TUfD5u3TpUgBA48aNASRUp5bnVXuuxkNzea6vUIhouKqWxCSanCnKJq0KVj3X1U7OT0tBa4pkTVmqXt0a681zwe3RBl0ZzGYpbdIDBgzAvHnz0LJlS1x99dX4r//6r/h0i/Lb3/42HjLQtWtXrFixAvv370eXLl1w6NAhHDhwAF26dMHcuXNxwgknoG3btmjUqBGAgulB5u2uTLCjE6vKSvA3K0/40czWrVtx8OBB5OTkYMuWLUXqR927d8e//vWvJIeN/Px8ZGVl4d1330X9+vXRqlUrVKtWDX379sWYMWPK5Zgcx3HKmlh+YZ5QKJib/+EPf4jMzEy89NJLxfa+3rBhAy6++GKsWbMGQEH839ixY7Fo0SLUr18fEyZMwAcffBB34KosMP3o97//fQCJ0aSm9AQSMX18Sa9fvx4A0K1bt/JpbCXl0KFDuPzyy3HNNdcgLy8PTz75JJ577rmk8m+F8e9//xu9evXC0qVL0bx5c8yaNQsLFy7EokWLMGHCBLz//vuYOXMmMjIy8D//8z849thj8cADD5TxUTlO6cJBaFHDVRcsWAAAccFD+65+Wt7XfG7RDqxKWO28RNUpBQ1LN65btw4AcMopp8TXoSOuZbdWT3VdjgqaIbJqhrTarMmpNJSWswIbN24EULlSK6dU0tnZ2ejevTu2bdtWquFRXbt2xddff43evXsjLy8PTZs2xe9+97tS235pQY/slStXJv2uCQSA8IVXFX608tBDD+H444/HVVddBQB4+eWXMX78eIwePTqt9du0aYOxY8diwIABOHz4MBo3bowpU6YAAG6++Wbce++96NKlCzIyMnDeeefhl7/8ZZkdi+M4TnmSUknv27cPvXr1wq9+9Succ8455dWuSgdf0vSk1DKUAELKkLMGqYqFOI7jFFdJE/pttGjRAkAinpqzfxqbTNVJFct4aEJ1qeV71TNal9+yZQuAgtBbAPGym0AiNpvboMqnmufvFDrchyaGou1YM5HpsVmKmdth6UrOet5yyy2obBSqpFesWIGbb74Z/fr1O6pf0ABw7rnnFnmd008/vQxa4jiO4xwtFPqS7ty5Mz755JPyaovjOI5TTK699loAicRMWi1L7bZUlVSfVKNcj2ihICpyqmCNMaa9mN+DFae4LhW1ZjnTXNr6d81QRqWtnulqL9e20Hu7KuQzsIN8HcdxHMepULwKluM4zhEEEzNNnToVANCuXTsAwPHHHw8gHHNMqDZVxRJNc0y4HDNGUklT3VKBA0CdOnWS9kFFzG1S+dKOzuW4b7aZbeXvWpGL6p7KmeV/P/vsMwDlV8GqNHAl7TiO4ziVlJTe3Y7jOE7pk5+fjxEjRqBt27a46aabABQoyfvvvx8rVqzA4cOHMXjwYDNlcVF59tlnAQAnnXQSgIRNmYqXcdL08qYapRLmJ9Wp1lxmnmutLx0MU2XMNG3StEHzk+uyDWwj/06VTu9u/p3bszKHXXrppalOT6XFlbTjOE45s2bNGvTt2xcvvvhi0u9z587F2rVr8dprr+HFF1/EY489hg8++KCCWulUBtwm7TiOUwJmzZqFJ554Iv59zZo1OPPMM5O8mglz08+aNQtXXXVVqKLg0qVLcfXVVyMrKwt169bF5ZdfjmeffRZnn312idvZq1evpO+LFi0CAJx44okAEiqW7abCtipTUVl/9913ABKKWfNCBOO+aXumrZk2Zs3ZrZnCqJDV5kyFTRs145+Z7ez66683z0dVwV/SjuM4JWDAgAFxR6TZs2dj3rx5eOqpp+LOT1H86le/AgC8/vrrSb9/9dVX8elooKB4xqefflr6jXaqDP6SdhzHKQWWLFmCRx55BM899xxWrlyZdpW3IHl5eaFUw4WVwy0Jl19+edL3hQsXAkgoa2ZXVEVN9btt2zYAiQxjTPurBCsgalUrqnON5dbaCFTMVOJUzMwTvnXrVgDA0KFDUx53VcNf0o7jOCXk3Xffxd1334158+ahUaNGaNSokVnlrTCaNGkSf+kBBS9AlqJ0jk78Je04jlMC1qxZg8GDB2Py5Mk47bTTSrStbt26Yd68eejSpQv27t2LRYsW4aGHHiqllhZOz549I3+fPXs2gLAt+rrrrktruyNHjgz99thjjwF4vXqUAAAAtElEQVRA3CZP1a61rqmgmWOb6n3QoEFp7ftIwF/SjuM4JeDee+9FTk4Oxo4dG58K7tChQ7Gq+vXr1w/r1q1Dly5dkJOTg2uvvRadOnUq7SY7VQiPk3Ycx3GcSorHSTuO4zhOJcVf0o7jOI5TSfGXtOM4juNUUvwl7TiO4ziVFH9JO47jOE4lxV/SjuM4jlNJ8Ze04ziO41RS/CXtOI7jOJUUf0k7juM4TiXFX9KO4ziOU0n5/0kGqXaDq2HXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 475.2x187.2 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The left \"temporal parietal\" ROI comprises auditory association\n",
      "cortex extending from anterior superior temporal cortex to the\n",
      "temporoparietal junction.\n"
     ]
    }
   ],
   "source": [
    "# Visualize the left temporal parietal ROI\n",
    "sns.set(palette='colorblind')\n",
    "roi_img = np.zeros(atlas_img.shape)\n",
    "for parcel in parcel_labels:\n",
    "    roi_img[atlas_img == parcel] = 1\n",
    "\n",
    "# Convert to a NIfTI image for visualization with Nilearn\n",
    "roi_nii = nib.Nifti1Image(roi_img, atlas_nii.affine, atlas_nii.header)\n",
    "\n",
    "# Plot plot left temporal parietal ROI\n",
    "plot_stat_map(roi_nii, cmap='tab10_r', cut_coords=(-53, -46, 10),\n",
    "              colorbar=False, title='left temporal parietal ROI');\n",
    "plt.show()\n",
    "\n",
    "# Print short \"figure caption\" describing visualization\n",
    "print('The left \"temporal parietal\" ROI comprises auditory '\n",
    "      \"association\\ncortex extending from anterior superior \"\n",
    "      \"temporal cortex to the\\ntemporoparietal junction.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once data is loaded, we divide the data into two halves for a two fold validation.\n",
    "We will use one half for training SRM and the other for testing its performance.\n",
    "Then, we normalize the data each half."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the number of subjects and TRs\n",
    "n_subjects = len(data)\n",
    "n_trs = data[0].shape[0]\n",
    "\n",
    "# Set a train/test split ratio\n",
    "train_test_ratio = .5\n",
    "test_size = int(n_trs * train_test_ratio)\n",
    "\n",
    "# Split/compile data into training and test halves\n",
    "train_data = []\n",
    "test_data = []\n",
    "for subject in np.arange(n_subjects):\n",
    "    \n",
    "    # Take the first chunk of TRs as training\n",
    "    train_data.append(zscore(data[subject][:-test_size, :], axis=0).T)\n",
    "    \n",
    "    # Take the second chunk of TRs as testing\n",
    "    test_data.append(zscore(data[subject][-test_size:, :], axis=0).T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating the SRM\n",
    "Next, we train the SRM on the training data. We need to specify desired dimension of the shared feature space. Although we simply use 50 features, the optimal number of dimensions can be found using grid search with cross-validation. We also need to specify a number of iterations to ensure the SRM algorithm converges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting SRM\n",
      "SRM has been fit\n",
      "Share response shape: 50 Features x 275 Time-points\n"
     ]
    }
   ],
   "source": [
    "# Set the number of features of shared space and number of iterations\n",
    "features = 50\n",
    "n_iter = 10\n",
    "\n",
    "# Create an SRM object\n",
    "srm = brainiak.funcalign.srm.SRM(n_iter=n_iter, features=features)\n",
    "\n",
    "# Fit the SRM data\n",
    "print('Fitting SRM')\n",
    "srm.fit(train_data)\n",
    "print('SRM has been fit')\n",
    "print(f'Share response shape: {srm.s_.shape[0]} '\n",
    "      f'Features x {srm.s_.shape[1]} Time-points')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After training SRM, we obtain a shared response $S$ that contains the values of the features for each TR, and a set of weight matrices $W_i$ that project from the shared subspace to each subject's idiosyncratic voxel space. Let us check the orthogonal property of the weight matrix $W_i$ for a subject. We visualize $W_i^TW_i$, which should be the identity $I$ matrix with shape equal to the number of features we selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASoAAAESCAYAAABKE5jvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVxU9d7A8c/AiBsqoaDeLCrXrmU+ZiFuuIUkDBipYaV2vZaay5O2uKGWV9PMIjWza49pFpWK+4b2aKkB5i520zITzEIYJVRin/N7/vBxEgVmZIA5g9+3r/N6ceYsv9+Mw5fffgxKKYUQQuiYm7MzIIQQtkigEkLongQqIYTuSaASQuieBCohhO5JoBJC6J6uA9U///lPPvnkE+v+mTNnaNmyJe+++671tYsXL/LAAw9w5cqVEu+TlpZGZGSkzfR69OjB8ePHiz02dOhQMjIybiH39lu9ejUxMTHFHvviiy9YsmRJhae7cOFCZsyYUSHplIdBgwYRFxdX5P/y119/ZcyYMbd0nxMnTtCrVy8iIiI4d+5cRWS1WBMnTmTp0qXFHgsPD+fy5ctlum9SUhLTpk0r9lhGRgbDhg2jT58+hIaGcvjw4TKloQe6DlRdu3blu+++s+5//fXXdO/enZ07d1pf27dvH+3ataNOnTol3qdhw4Z8+eWXDuUlPj7eoetLc+jQIXJzc4s9NnDgQF544YVKT1evrv+//P333zlz5swtXb9z5078/f1Zu3YtTZo0qYgs3rINGzZQt27dMl37888/k5aWVuyxN954g/bt27N161befvtt/vu//5ucnBxHsuo0RmdnoDRdu3Zl0aJFaJqGm5sbX3/9NePGjWP8+PGcPXuWu+++m8TERLp16wZcLTnNmDGD1NRUCgoKCAkJYcSIEZw7dw6TycSRI0fIyclh+vTpHDt2jDp16tCsWTMA5syZA8DKlSuZPn06GRkZhIeHM27cOCZNmgTAkCFDWLJkCY0bN7bmceHChZw9e5a0tDTMZjOtW7fG39+f9evXc+7cOV599VVCQ0O5cOEC06ZN4+LFi5jNZu68807ee+89Dh8+zK5du4iPj6dGjRpkZGRw9OhR0tPTadmyJX5+fvzxxx+8+OKL9O3bl1mzZhEYGMh7773HsWPHWLp0KW5uf/29KSgoYM6cOSQmJuLu7k6bNm2YNGkSnp6e9OjRgzZt2vDjjz8yfvz4IukC/PLLLwwaNAiz2UyDBg1499138fX15dSpU8yYMYPMzEwMBgNDhw6lb9++ACxZsoTY2Fhq165N+/bt2blzJ7t27eLKlSu88cYbnDx5EoPBQJcuXRg/fjxGo5EHH3yQF154gfj4eNLT0xk2bBhPP/002dnZvP7666SkpJCZmUnt2rWZN28e9913n/X9Xfu/PHjwIFFRUaSlpfHPf/6T9u3b8/PPP/POO+8AcPDgQWbOnMn69eut127cuJEvvvgCi8VCbm4u77zzDosWLWLLli24u7tz7733MnXqVHx8fBg0aBD16tXjl19+YeDAgQwaNMh6H7PZzIQJE/jjjz8ACAwM5KWXXmLt2rVs376df//73wA37R86dIjt27eTlZVFp06dmDBhAkajkZYtW5KYmIi3tzerV6/miy++QNM0vLy8mDp1Kk2bNuXPP/9k5syZHD58GHd3d3r16sXAgQNZsGABV65cYdKkScyePduax8LCQr755humT58OwP33388999zD3r17CQoKKtsvpBPpukR17733UrduXX788UcuXbrEmTNnaNu2LV27dmXXrl0AJCYmEhgYCMCrr77Kk08+ydq1a4mNjSUhIYGtW7cWuecHH3yAxWJh27ZtLF++nB9++KHI8erVq7N27VpWr17Nxx9/TGpqqvUL8MknnxQJUtccOnSIRYsWsW7dOvbs2cPp06eJiYlh6tSpLFy4EIAtW7bQtm1bVq5cyc6dO6lRowYbNmzgscceo0ePHjz33HM888wzAPz222+sW7eOefPmWdNo0KABc+bMYerUqXz11VesX7+ed955p0iQAli8eDHp6els2LCBDRs2oGkac+fOtR5v3rw527ZtKzbdX3/9lfnz5xMXF0fdunVZvXo1hYWFjBw5kkGDBrFp0yY++ugj3n33XY4cOcLevXutn/XatWv5888/renMnDkTLy8vNm3axJo1a/jxxx/5+OOPAcjPz+eOO+7gyy+/ZMGCBcyePZu8vDz27NlD3bp1WblyJdu3b+eBBx4osUrs7u7OzJkzufvuu1m6dCkDBgzgm2++ITMzE4BVq1bdVN0PCwsjMjKSPn368M4777BmzRr27t1LbGwsmzZtonnz5kycONF6ft26ddm6dWuRIHXt3k2aNGHdunXExMSQkpJSatPDNefPn2f58uWsX7+ekydPsmrVqiLH9+/fz/r164mJiWH9+vUMGzaM0aNHA7BgwQLy8vLYunUr69ev5/Dhw5w9e5axY8fSvn37IkEK4I8//kDTNLy9va2vNWzYkPPnz9vMpx7pOlDBX9W/PXv20LFjR9zc3OjevTvffvst586dw2Aw0LRpU7Kzszlw4ADz588nPDycAQMGkJqaysmTJ4vcb/fu3fTr1w83Nzc8PT154oknihwPDQ0FwMfHhwYNGnDx4kWbeezYsSN16tShRo0a+Pr60qVLFwDuvvtu6y/OkCFDaNeuHcuWLeP111/n1KlTZGdnF3u/tm3bYjTeXNjt3Lkzffr0YcyYMcybN6/Il/CaPXv2EBkZSbVq1XBzc2PQoEHs3bvXerx9+/Ylvo9OnTpZ79mqVSsyMjJITk4mLy/P+le4YcOGBAUFsXfvXnbv3k1wcDB169bFYDBYA961fDz77LMYDAY8PDyIjIxkz5491uM9e/YEoHXr1uTn55OdnU1wcDBPPPEEn376KTNnzmT//v0lfkY3ql+/Pt26dWPDhg1cunSJb7/9FpPJVOo1e/bsISIiglq1agEwePBg9u3bR35+fqmfVZcuXdixYwfPP/88K1eu5OWXXy616eGa8PBwatWqhYeHB2FhYSQkJBQ5/s0335CSkkJkZCTh4eG8/fbbXL58mczMTBISEujXrx/u7u54eHjw2Wef4e/vX2JamqZhMBiKvKaUwt3d3WY+9UjXVT+4GqhiY2OpXr269csdEBBAVFRUkWqfpmkopfjyyy+pWbMmcLUxsXr16tYiOoDRaOT66Y03lkiuDxAGgwF7pkJ6eHiUeI9r3n77bZKSknjyySfx9/ensLCwxHtf+8W5kVKK06dP06BBA44ePVrsL9KNX1BN0ygoKLB57xvzfe29WyyWYr/whYWFN32W1/8SFJePwsJC63716tWt6Vy75+eff86qVat45plnMJlMeHl53VKD9zPPPMPrr7+O0WgkKCiI2rVrl3q+rTyW9Fm1adOGnTt3kpiYyL59++jfvz8fffTRTd+X6z93KPr5KKVu+p5omkZ4eDivvvqqdT89PZ169ephNBqL5DU1NdVaZS9O/fr1UUqRmZmJl5cXAOnp6TRs2LDEa/RM9yUqf39/Tpw4wf79+60llRo1atC6dWs+++wza7XP09OTtm3bsmzZMgAuX77MwIEDizS8w9X2hDVr1qBpGjk5OWzevPmmX8TiuLu7F/kS36pvv/2WIUOG0LdvX+rXr09CQgIWi+WW7r18+XKys7NZs2YNy5cvJykp6aZzunTpwhdffEFBQQGaphETE0OnTp3K/J7uu+8+jEYjO3bsAK62A27fvp2OHTsSGBjIjh07rNWe2NhY63WdO3fms88+QylFfn4+q1atomPHjqWm9e233/LEE0/Qv39/7r33Xnbt2mX9jErK//XBoF27dri5ubF06VK7enm7dOnCmjVrrKW2Tz/9lEceeeSmPzw3mjdvHh988AG9evViypQpNGvWjFOnTuHt7c2pU6fIy8ujoKCA7du3F7luy5Yt5Ofnk5eXx7p16+jatWuR4507d2bLli2kp6cDV3t8hwwZAlz947xu3To0TSM/P5+xY8dy4MCBEv8PjUYj3bp1s1YvT548yenTp0sthemZ7ktUNWvW5J577qGgoKBI8TowMJC33367yAc/b948/vWvf2EymcjPzyc0NJSwsLAif5WHDx/OjBkzMJlM1KlTh/r165f6l+ma4OBgBg0axMKFC2nRosUtv49Ro0Yxd+5c5s+fT7Vq1WjXrh1nz54FrpYarzXml+SHH37gww8/JDY2loYNGzJ58mRefvll1q1bh6enp/W8kSNH8tZbb9G3b18KCwtp06YNU6dOLfae9qRbrVo1PvjgA2bOnMnChQuxWCyMGjWKDh06ADBgwACeeuopatSoQfPmza2l2aioKGbOnInJZKKgoIAuXbowYsSIUtMaOnQo06ZNswa8tm3b8tNPP5V4frNmzahevTr9+vVj9erVGAwGIiIi2Lp1K61atSo1LYB+/fqRmppK//790TQNPz+/Iu2CJRkyZAgTJ04kNDQUDw8PWrZsSUhICG5ubjzyyCM8/vjj+Pj44O/vz48//mi9rkmTJjz99NP8+eefPPbYYzc1O3Tu3Jnnn3+eoUOHYjAY8PT05P3338dgMDB69GhmzZpFeHg4FouFPn36EBQUREpKCosWLWL06NG8//77Re43ffp0oqKiCA0NxWAwMHfuXLuqqHpkuN2WedmyZQuenp4EBgaiaRpjxoyhU6dOPP30087Omss5fvw4R44cYfDgwQAsW7aMY8eO8d577zklP4WFhYwePZqwsDD69OnjlDzcqqysLB5++GGOHTtm1x/M25Xuq37lrXnz5ixevJjw8HBCQ0Px9fWlf//+zs6WS7r33ns5ePAgoaGhmEwmEhMTrUM5KtvPP/9MQEAAd9xxB8HBwU7Jw61KSkoiODiYiIgICVI23HYlKiGE67ntSlRCCNcjgUoIoXsSqIQQuieBSgihexKohBC6p9tAtWnTJuugtpImpjpbVlYWoaGh1gGlCQkJmEwmgoKCiI6OdnLu/vL+++8TEhJCSEiIdYKyXvMKMH/+fPr06UNISIh1poGe8wvw1ltvWSc0nzhxgoiICHr37s2UKVMcmtEg/p/SofPnz6vu3burP/74Q/3555/KZDKpU6dOOTtbRRw9elSFhoaq1q1bq19//VXl5OSowMBAdfbsWVVQUKCGDh2qvvnmG2dnU8XHx6unnnpK5eXlqfz8fDV48GC1adMmXeZVKaW+++47FRkZqQoKClROTo7q3r27OnHihG7zq5RSCQkJyt/fX02YMEEppVRISIg6cuSIUkqpSZMmqZiYGGdmr0rQZYkqISGBDh064OXlRa1atejduzdxcXHOzlYRq1atYvr06fj6+gJXB+/5+flx1113YTQaMZlMusizj48PEydOxMPDg2rVqtG0aVOSk5N1mVeARx99lBUrVmA0Grl48SIWi4XLly/rNr+ZmZlER0dbpwf99ttv5Obm0rZtWwAiIiJ0k1dXpstAlZ6ejo+Pj3Xf19e3xFUMnWXWrFlFVi/Qa56bN29u/aVJTk5m27ZtGAwGXeb1mmrVqrFgwQJCQkIICAjQ7WcLMG3aNMaNG2ddofPGvPr4+Ogmr65Ml4HqxuU3lFJ2rXDgTHrP86lTpxg6dCivvfYad911l67zCjB27FgSExNJTU0lOTlZl/ldvXo1jRs3JiAgwPqa3r8HrkqXqyc0atSIgwcPWvfNZrO1iqVXjRo1wmw2W/f1lOdDhw4xduxYJk+eTEhICPv379dtXk+fPk1+fj73338/NWvWJCgoiLi4uCJrOeklv1u3bsVsNhMeHs6lS5fIzs7GYDAU+WwvXLigi7y6Ol2WqDp27EhiYiIZGRnk5OSwY8eOm9bu0ZuHHnqIM2fOkJKSgsViYfPmzbrIc2pqKqNGjWLevHmEhIQA+s0rXF0TPSoqivz8fPLz89m5cyeRkZG6zO+yZcvYvHkzGzZsYOzYsfTo0YPZs2dTvXp1Dh06BFx9cIMe8urqdFmiatiwIePGjWPw4MEUFBTQr18/2rRp4+xslap69erMmTOHMWPGkJeXR2BgoC5m8S9dupS8vLwi605FRkbqMq9wdZ2xpKQk+vbti7u7O0FBQYSEhODt7a3L/BZn3rx5REVFkZWVRevWra3L4Iiyk9UThBC6p8uqnxDi9nDjoOmSSKASQjjFsWPHGDhwIMnJyTbPdUob1aZNm1i8eDGFhYUMGTKkyGOWcnNz+f777/Hx8XHZR/sIoXcWiwWz2cwDDzxQ5tVFMzMzycrKsutcT09P69Nwrrk2aPq1116zeX2lB6q0tDSio6NZu3at9Xlv/v7+1icWf//990UClxCi4sTExJT6rMeSZGZm8livnly+Yl+gqlevHjt27CgSrGbNmmV3epUeqK6fHgNYp8dceyLstVG9nyx6m0a+DQBo6S9rmgtRnozuBprcWbvIKPpbkZWVxeUrWaz44G0a+pZ+j7R0M4NffJWsrKybSlX2qvRAVdx0iOufT3etutfItwF3Nr76sMTCQumYFKIiONq80rCBN3c2bFD6SVrJz2a0V6UHKpliIEQVYim8utk6x0GVHqjsnR7T0r+/tSSV+/vem47X+FuXisukEMIuCoVSms1zHFXpwxNccXqMEKIEmmbfVopdu3bRpEmTUs+p9BKVK06PEUKUQGlXN1vnOMgp46hMJhMmk8kZSQshypOm2W4st1GisocuJyULIVxEVS5R3ariGs5vbGCXxnUhKp/SClE2evWU5oK9fkKIKkRTtqt2muO9fhKohBBlJ1U/IYTuaRY7GtNdcGR6ebmxTUrarIRwAilRCSF0T7PYniJzO5eohBA6YMfIcxlHJYRwKqU0lCq9xGRrLqA9qkygstVmVdw5QggHSRuVEEL3lB1VPwlUQginkhKVEEL3LIVgKbB9joMkUAkhyk56/RwjE5mFqARS9RNC6J6yY1KykknJQghnkqqfEELvlKUQZaMx3dZ6Vfa4rQKVDAoVopxJG5UQQvdkwKcQQvekRCWE0D1ZilgIoXtSoqp4MihUCAdZCqHQRq+e9PoJIZxKKTtKVFL1E0I4kwz4FELonrRROYc83UaIWyAlKiGE/tnRRoW0UQkhnKnQjl4/W8ftIIFKCFF2Stnu1ZNev4onE5mFKIW0UQkhdE8WzhNC6J4MTxBC6J7FcnWzdY6D3By+gw1ZWVmEhoZy7tw5ABISEjCZTAQFBREdHV3RyQshKtK11RNK3XRe9Tt27BhRUVEkJycDkJuby+TJk/n0009p3Lgxw4cPZ/fu3QQGBlZkNsqVTGQW4jqVtHBehZaoVq1axfTp0/H19QUgKSkJPz8/7rrrLoxGIyaTibi4uIrMghCiIl1ro7K1OahCS1SzZs0qsp+eno6Pj49139fXl7S0tIrMghCiImkKZatqp/eq3400TcNgMFj3lVJF9oUQLqYqjqNq1KgRZrPZum82m63VQlcmE5nFbUvTbPfqlUOgqvBev+s99NBDnDlzhpSUFCwWC5s3b6Zr166VmQUhRHmy2eNnR4nLDpVaoqpevTpz5sxhzJgx5OXlERgYSHBwcGVmQQhRnqpS1W/Xrl3WnwMCAti4cWNlJCuEqGgyKdl1yURmcduQuX5CCN3TlO3hB642PEEIUcVU0lw/CVRCiDJTmoayUfWzddweEqiEEGWnsF21c7zmJ4GqMshEZlFlyXpUQgjdk8Z0IYTuWSxQKI3pQgg9U3Y810/GUbkumcgsqgSp+gkh9E6GJwgh9E9KVEII3VN2BCppo6o6ZCKzcEkWOxbOs0jVTwjhRMqONdNtrqluBwlUQoiyk6qfEEL3qtIKn0KIKkp6/W5vMpFZuAQJVEIIvVOaQtno1ZPGdCGEc0mJSgihd7ocnpCWlsavv/5K+/btHU5Y3DqZyCx0Ry8lqs8//5xDhw4xZcoUIiIi8PT0JCgoiJdfftnhxIUQLk4BtkYflMNSxDYf6R4bG8ukSZOIi4ujZ8+ebNmyhfj4eMdTFkK4PFWo2bU5ymagMhgMNGjQgMTERDp06IDRaEQrhwFcQogqQLNzc5DNqp+HhwcfffQR+/fvZ+bMmXz++efUrFnT8ZSFw2Qis3A2pexoTC+HKTQ2S1SzZs0iOTmZt956i3r16nHo0CFmzpzpcMJCiCpALyWq++67j6lTp5KSkoJSipkzZ0qJSghxlR3DE8qj189miero0aP06tWL4cOHk5aWRrdu3Th8+LDDCQshqoBKKlHZDFRz585l+fLleHl50ahRI+bOncusWbMcT1kI4fKUBVShjc3xp2XZrvrl5ubSrFkz635gYCDR0dGOpyzKnUxkFpWtkh6UbDtQGY1GLl26hMFgAOCXX35xPFUhRNVgT9WuMgLViBEjePbZZ7lw4QLjx48nPj6eGTNmOJ6yEMLlVdLzR20Hqi5dutC0aVPi4+PRNI1Ro0bRtGlTu27+/vvvs23bNuBqlfG1114jISGB2bNnk5eXx+OPP864ceMcewdCCKfRTdWvX79+bNiwAT8/v1u6cUJCAt9++y3r1q3DYDAwbNgwNm/ezLx58/j0009p3Lgxw4cPZ/fu3QQGBpb5DYjSyaBQUZGUxYD6/2ah0s5xlM1ev5o1a3L+/PlbvrGPjw8TJ07Ew8ODatWq0bRpU5KTk/Hz8+Ouu+7CaDRiMpmIi4srU8aFEM53repX6lYZVb+cnBx69uxJo0aNqFWrlvX1TZs2lXpd8+bNrT8nJyezbds2nn32WXx8fKyv+/r6kpaWVpZ8CyF0QGl2lKg0x0tUNgPVlClTHErg1KlTDB8+nNdeew13d3eSk5Otx5RS1t5EIYQL0kDZ+hWujDaqFi1alPnmhw4dYuzYsUyePJmQkBD279+P2Wy2Hjebzfj6+pb5/kII51LKgLIRqWwdt4fNQNWhQwcMBkOR0o+Pjw979uwp9brU1FRGjRpFdHQ0AQEBADz00EOcOXOGlJQUmjRpwubNm3nyyScdfhPCfjIoVJQnZUeJqlJ6/U6ePGn9OT8/n82bN3PmzBmbN166dCl5eXnMmTPH+lpkZCRz5sxhzJgx5OXlERgYSHBwcBmzLoRwNk0zoFF6pNIqo43qeh4eHkRERBAREWFzKeKoqCiioqKKPbZx48ZbSVYIoVNKM6BsBKpKaUzPzMz8K0Gl+P7777l8+bLDCQshXJ9uAtX1bVQA9evXd7gnUOiLPN1GlJVStsdJVco4qh9++AE3t6LjQi9duuR4ykII16cMtktM5dDrZ3NkenG9cs8884zDCQshXN/VEpXBxuZ4OiWWqIYMGcLx48fJzc2lXbt21tc1TePBBx90PGUhhMuzWAxYbAzatlTkOKpFixaRmZnJ5MmTmT179l8XGI1FpsGIqkcmMgt7KWVHY3pFBipPT088PT1ZsWLFDYkqUlJSuOeeexxOXAjh2uwamY4BG7HMJpuN6V9++SVz584lJyfH+pq3t7c8LVkIcbWNytY5UPGBasmSJSxbtozFixfz0ksv8fXXX5dp2RchRBVkxzgqMNjRbVc6m5d7eXnx0EMPcf/993Px4kVGjhzJgQMHHEtVCFElWDQ3uzZH2f1wBz8/P5KSkujUqRMWSzk8/0a4DJnILEpid9XPQTZD3YABAxg+fDjdunVj5cqVRERE2L1muhCiatMwoCkbm6MNVNi5ZnqfPn2oVasWK1eu5Pjx43Tu3NnhhIUQru/agM9SzymHdGyWqDRN44svvmDChAnUrl2b06dPU61atXJIWgjh8tRf8/1K2sojUtksUc2dO5eMjAyOHz+OUoq9e/diNptLXMJF3B5kIrMArNW7Us8ph6qfzRJVYmIic+bMoXr16tSpU4ePP/5YxlAJIQCd9fpdv3qCh4cHRuMtrbcnhKii7KnZlUcblV0Pd4iJicFisfDLL7+wfPlyWrVqVQ5JCyFcnbKj6md7QKhtdj0u68033+TixYsMHDiQLl26yMJ54iYykfn2ZPdcPweVGKimTJnCrFmz2Lt3L2+++abDCQkhqh4N24/tK4eH0JQcqBISEjh8+DALFizAz8/PuhTxNa1bty6H5IUQrkxhxzIvFVmiGjBgAK+99hrnz59n9OjRRY4ZDAZ27tzpcOJCCNdmwUChjaqfrYX17FFioBo5ciQjR45k3LhxREdHO5yQEKLqudrrV/Ej0202pkuQEmUhE5lvD05voxJCCFuc3kYlhBC2SIlKCKF7GgYsNkpMFbrMi8lkKvXCTZs2OZy4uL3IROaqRwNsPX+0QktU2dnZ5OXlERYWRpcuXXB3dy+H5IQQVYmG7YXxKrREtXPnTg4ePMi6det444036NGjBxERETRr1szhRIUQVYMuJiW3b9+e9u3bk5uby1dffcXs2bPJysoiPDycp59+uhySF0K4MoXtql2ljKMCqFGjBo8//ji1atVi2bJlREdHS6ASDpOJzK5PM9hR9avIkenXHD16lPXr1/PVV1/RunVrBg4cSK9evRxOWAjh+iz/v9k6x1ElBqr333+fjRs3UqtWLfr27cuGDRto0KBBOSQphKgqNIMd46gcL1CVHqj+9re/0ahRI/bt28e+ffuKHP/www8dT10I4dKUHb1+FToyffbs2Q7ffP78+Wzfvh2DwUC/fv34xz/+QUJCArNnzyYvL4/HH3+ccePGOZyOEMI5nN7r98QTT5R4kT0Pd9i/fz/79u1j48aNFBYW0qdPHwICApg8eTKffvopjRs3Zvjw4ezevZvAwMCy5V5UKTKR2fVUVtWvxMdD/Oc//yEyMpIRI0aQkZEBwO+//87o0aMZOXKkzRs/+uijrFixAqPRyMWLF7FYLFy+fBk/Pz/uuusujEYjJpOJuLg4x9+FEMIpNDs3R5UYqF5//XWCgoJo0qQJixcv5n//938JCwsjJyeHDRs22HXzatWqsWDBAkJCQggICCA9PR0fHx/rcV9fX9LS0hx/F0IIp9AMYLGxVWhj+pUrVxg6dCgWi4XevXuzbds23njjDUJCQm4pgbFjx/L8888zYsQIkpOTMVw3pkIpVWRfCOFanL56Qs2aNQFwd3cnLy+PJUuW8Pe//93uG58+fZr8/Hzuv/9+atasSVBQEHFxcUXmDJrNZnx9fR3IvqjqZFCovlVWoAZpiAAAAAyxSURBVCqx6nf9wxzuuOOOWwpSAOfOnSMqKor8/Hzy8/PZuXMnkZGRnDlzhpSUFCwWC5s3b6Zr165lz70QwqmUwb7NUSWWqDRN49KlS9aAdf3PAF5eXqXeODAwkKSkJPr27Yu7uztBQUGEhITg7e3NmDFjyMvLIzAwkODgYMffhRDCKZw+1++nn36iQ4cO1uDk7+9vPWYwGDhx4oTNm48ZM4YxY8YUeS0gIICNGzeWNb9CCB1x+hSakydPlsPthRBVmdOn0MDVdiqLxYLRaCQrK4uEhARatmyJn5+f4ykLUQYyKFRfnN6Y/vPPP9OzZ0/27t1Lbm4u/fv3Jzo6mmeffdaukelCiKrvWhtVaVt5tFGVGKjmzp3LSy+9RPfu3dmyZQsAW7ZsYdWqVSxcuLAckhZCuDpl5+aoEqt+qamphIWFAfDdd9/Rs2dP3NzcaNy4MVlZWeWQtBDC1Tm9jcrN7a/C1pEjR4iKirLu5+XlOZ6yEOVEnm7jPE7v9atXrx4nT54kKysLs9nMI488AsDhw4dp2LBhOSQthHB1CoVmo3KnyqHyV2KgGj9+PM899xxZWVm88sor1KpVi6VLl/Lhhx+yaNEihxMWQrg+p8/1a9u2LXv27CE3N5e6desC8F//9V+sXr2ae+65pxySFkK4OqcvnAfg4eGBh4eHdb9du3blkKQQFUsmMlcep5eohBDCFguKQkPpZSZLRbZRCSGELbqo+gkhRGmk6ieE0D3NjuEJto7bQwKVqPJkInPFKo+qnS0SqIQQZSZVPyGE7llQNnv1pNdPCOFUUqISogLJRObyouyYyyclKiGEE0mJSgihe1cDla3hCY6TQCWEKDMZmS5EJZKJzGVjQVEovX5CCD1TdjSmV+jCeUIIYYs0pgshdE9KVEII3ZMSlRBOJBOZ7aOhsCgbwxNsLKxnDwlUQogyk2VehBC6J21UQgjdkzYqIXRGJjLfzOkPIBVCCFuuTqGxFagcJ4FKCFFmFmW71688ptC4OXwHG9566y0mTpwIwIkTJ4iIiKB3795MmTKFwsLCik5eCFGBrvX62docVaGBKjExkXXr1ln3X331VaZNm8b27dtRSrFq1aqKTF6IClXjb12KbLm/771pq+oUfzWol7SVR9WvwgJVZmYm0dHRjBgxAoDffvuN3Nxc2rZtC0BERARxcXEVlbwQohIoO/85qsIC1bRp0xg3bhx169YFID09HR8fH+txHx8f0tLSKip5IUQlcOmq3+rVq2ncuDEBAQHW1zRNw2AwWPeVUkX2hRCuRyll1+aoCun127p1K2azmfDwcC5dukR2djYGgwGz2Ww958KFC/j6+lZE8kKISqLZ8bgs3U6hWbZsmfXntWvXsn//fmbPnk1oaCiHDh3i4YcfZsOGDXTt2rUikhfCKW7HicxVcq7fvHnziIqKIisri9atWzN48ODKTF4IUc6UwmbVziUGfEZERBAREQFAq1atiI2NregkhRCVpEqWqIQQVYusniBEFVDVJzJrdkyhkRKVEMKpNBQGqfoJIfRMApUQQvfsGdApbVRCuJiq9kRmWThPCKF7snCeEEL3LEpDqdLn7GrlsGq6BCohRJlJG5UQQve0/6/82T7HMRKohHAi15/IbM/CeBKohBBOpKEw2Kr6ySPdhRDOpOyo+kkblRDCqSxKYVCl9+rpdoVPIUTZudKgULuWGlYKRxcdl0AlhCgze6p+IIFKCOFE2tUlPks/SSmHnyIjgUoIUWb2lqgcJYFKCFFmmtJstlHZGr5gDwlUQuicngeFXm1Mt3mSw+lIoBJClJlmxygpWwvr2UMClRCizKREJYTQPc2OQCVtVELcpvTydBt7piRLr58Qwqk01NWxVKVwdAwVSKASQjjAnjaqq8cdG5sugUoIUWaaUmh2jfeUQCXEbc9ZE5ntL1E5RgKVEKLMrj4uq3SOTkgGCVRCCAdIiUoIoXsWpWGp+DnJEqiEEGWnKWwOTzDYeO6fPSRQCVEF2ZrI/FtqGr37/8PhdOx6rl85lKjKYyyWEOK2pWz+K6nut2nTJvr06UNQUBAxMTGlpiIlKiFEmSnseVLyzdLS0oiOjmbt2rV4eHgQGRmJv78/zZo1K/YeugtUFosFAKN7eXRqCiGu+S01zfrz+fQLwF+/b2Xl5gZu7rbPuVFCQgIdOnTAy8sLgN69exMXF8fo0aOLvYfuApXZbAagyZ21nZwTIaqW4tqkzGYzfn5+t3wvT09P6tWrZ/f59erVw9PT07qfnp6Oj4+Pdd/X15ekpKQSr9ddoHrggQeIiYnBx8cHd3cboVoIUSYWiwWz2cwDDzxQpuu9vLzYsWMHWVlZdp3v6elpLT0BaJqGwfBXrUkpVWT/RroLVDVq1KB9+/bOzoYQVV5ZSlLX8/LyKhJ8bkWjRo04ePCgdd9sNuPr61vi+dLrJ4SodB07diQxMZGMjAxycnLYsWMHXbt2LfF83ZWohBBVX8OGDRk3bhyDBw+moKCAfv360aZNmxLPN6jyeDC8qHRHjx7lnXfeITMzE6UUjRo1YsKECTRv3hyAli1b0qJFC9zc3DAYDOTk5ODp6cnrr7/Ogw8+yNq1a5k0aRKjRo1i7Nix1vsqpejVqxc1a9Zk8+bNN6W7ePFiVq5cSUBAALNnz77lfF+5coVRo0axYsWKsr95cftRwuXk5eWpRx99VH3//ffW19avX68CAwNVYWGhUkqpFi1aqIsXLxa57n/+53/UgAEDlFJKrVmzRnXr1k317NmzyDn79+9XHTt2VCEhIcWm3aNHD3XgwIEy5/3XX39Vbdu2LfP14vYkbVQuKCcnhytXrpCdnW19LSwsjKlTp5Y4LqawsJDU1NQiXcotWrSgVq1aHD582PraunXrCAsLK/YeL730EmlpaUyZMoWtW7dy5coVJk6cSEREBCaTiTfffJPCwkIAYmNj6d+/P3379qV79+58/vnnAEyaNInc3FzCw8OxWCy0bNmSjIwMaxrX9r/77jvCwsKIjIzEZDKRn5/Prl27rPeMjIzkyJEjAJw+fZrIyEgiIiJ44oknbI5yFi7I2ZFSlM3HH3+s2rRpo3r06KFeeeUVtXr1apWdnW093qJFCxUaGqpCQ0NVp06dVI8ePdS//vUvdeHCBaXU1RLVCy+8oJYuXaqmTZumlFIqOztbBQUFqfj4+BJLVN27d1dJSUlKKaUmTpyoVqxYoZRSqrCwUL3yyitqyZIlKisrSw0YMEBlZGQopZQ6cuSItRR1Y4nqxpLftf19+/apVq1aqXPnzimllDpz5owKDQ213vOnn35SnTp1Un/++aeaNGmS+ve//62UUio9PV299NJLymKxOPgJCz2RxnQX9Y9//IP+/ftz4MABDhw4wEcffcRHH31EbGwsderUAeCTTz7B29ub//znP7zwwgv4+/tTv379IvcxmUyEh4czZcoUvvrqK3r06GH3+LVvvvmG48ePExsbC0Bubi4AtWvX5sMPP2T37t0kJydz8uTJIqU/ezVu3Jg777wTgPj4eNLT03nuueesxw0GA2fPnuWxxx5jwoQJJCUlERAQQFRUFG7FDYcWLksClQs6dOgQR44cYdiwYXTv3p3u3bszfvx4QkNDiY+PJzg4uMj5rVu3ZtKkSUycOJH777+fJk2aWI/5+Pjw97//nT179rB+/XomTpzIH3/8YVc+NE1j/vz5NG3aFIDLly9jMBg4f/48Tz31FAMGDODhhx8mODiYr7/+2ub98vPzi+zXqlWrSFoBAQG899571tdSU1Px9fWlVatWbN++nYSEBBITE1m0aBFr166lUaNGdr0PoX/yZ8cFeXt7s3jx4psGzGVlZdGiRYtirwkNDaVNmzbF9tT17duXZcuWceXKlRKvL07nzp1Zvnw5Siny8/MZOXIkn332Gd9//z3e3t68+OKLdO7c2RqkLBYLRqMRi8Vincjq7e3N8ePHAYrtZbwmICCA+Ph4Tp8+DcDu3bsJCwsjNzeXl19+ma1btxISEsL06dPx9PTk7Nmzdr8PoX9SonJB9957L4sWLSI6Oprz589TvXp16tSpw5tvvsl9991X4nVTp04lLCyMvXuLLvzfq1cvpk+fzrhx424pH1OmTGHWrFmYTCYKCgro2LEjw4YNo7CwkNjYWIKDgzEYDDz66KN4e3uTkpKCn58fbdq0ISQkhJiYGKKiopgxYwZ169alY8eOReZ/Xa9Zs2bMmDGD8ePHo5TCaDSyePFiateuzYsvvsiUKVNYuXIl7u7u9OrVi0ceeeSW3ovQNxlHJYTQPan6CSF0TwKVEEL3JFAJIXRPApUQQvckUAkhdE8ClRBC9yRQCSF07/8AaDcDrOKr14MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weight matrix shape: 935 Voxels x 50 Features\n",
      "\n",
      "This test confirms that the weight matrix for subject 0 is orthogonal.\n"
     ]
    }
   ],
   "source": [
    "# Use the first subject as an example\n",
    "subject = 0\n",
    "\n",
    "sns.set_style('white')\n",
    "fig, ax = plt.subplots(1)\n",
    "m = ax.matshow(srm.w_[subject].T.dot(srm.w_[subject]))\n",
    "ax.set_title(f'Weight matrix orthogonality for subject {subject}', pad=10)\n",
    "ax.set_xlabel('SRM features')\n",
    "ax.set_ylabel('SRM features')\n",
    "ax.tick_params(length=0)\n",
    "cbar = fig.colorbar(m, ax=ax, ticks=[0, 1])\n",
    "cbar.ax.tick_params(length=0)\n",
    "plt.show()\n",
    "\n",
    "print(f'Weight matrix shape: {srm.w_[subject].shape[0]} '\n",
    "      f'Voxels x {srm.w_[subject].shape[1]} Features\\n')\n",
    "\n",
    "# Check against identity matrix\n",
    "if np.allclose(np.identity(features), srm.w_[subject].T.dot(srm.w_[subject])):\n",
    "    print(\"This test confirms that the weight matrix for \"\n",
    "          f\"subject {subject} is orthogonal.\")\n",
    "else:\n",
    "    print(\"Weight matrix is not orthogonal.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Between-subject time-segment classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we trained SRM above, we learned the weight matrices $W_i$ and the shared response $S$ for the training data. The weight matrices further allow us to convert new data to the shared feature space. We call the `transform()` function to transform test data for each subject into the shred space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform subject-space test data into shared space\n",
    "test_shared = srm.transform(test_data)\n",
    "\n",
    "# z-score the transformed test data\n",
    "test_shared = [zscore(ts, axis=1) for ts in test_shared]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We evaluate the performance of the SRM using a between-subject time-segment classification (or \"time-segment matching\") analysis with leave-one-subject-out cross-validation (e.g. [Haxby et al., 2011](https://doi.org/10.1016/j.neuron.2011.08.026); [Chen et al., 2015](https://papers.nips.cc/paper/5855-a-reduced-dimension-fmri-shared-response-model). The function receives the data from `N` subjects with a specified window size `win_size` for the time segments. A segment is the concatenation of `win_size` TRs. Then, using the averaged data from `N-1` subjects it tries to match the segments from the left-out subject to the right position. The function returns the average accuracy across segments for each subject."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_segment_classification(data, win_size=10): \n",
    "    n_subjects = len(data)\n",
    "    (n_features, n_trs) = data[0].shape\n",
    "    accuracy = np.zeros(shape=n_subjects)\n",
    "    n_segments = n_trs - win_size + 1\n",
    "    \n",
    "    # Set up container for training data\n",
    "    train_data = np.zeros((n_features * win_size, n_segments), order='f')\n",
    "    \n",
    "    # Training data (includes test data, but will be removed)\n",
    "    for m in range(n_subjects):\n",
    "        for w in range(win_size):\n",
    "            train_data[w * n_features:(w + 1) * n_features, :] += \\\n",
    "                data[m][:, w:(w + n_segments)]\n",
    "            \n",
    "    # Analyze each subject (leave-one-out)\n",
    "    print(\"Between-subject time-segment classification accuracy \"\n",
    "          \"for each subject:\", end=' ')\n",
    "    for test_subject in range(n_subjects):\n",
    "        test_data = np.zeros((n_features * win_size, n_segments), order='f')\n",
    "        for w in range(win_size):\n",
    "            test_data[w * n_features:(w + 1) * n_features, :] = \\\n",
    "                data[test_subject][:, w:(w + n_segments)]\n",
    "\n",
    "        A = np.nan_to_num(zscore((train_data - test_data), axis=0))\n",
    "        B = np.nan_to_num(zscore(test_data, axis=0))\n",
    "\n",
    "        # Compute correlation matrix\n",
    "        correlations = compute_correlation(B.T, A.T)\n",
    "\n",
    "        # Correlation-based classification\n",
    "        for i in range(n_segments):\n",
    "            for j in range(n_segments):\n",
    "                \n",
    "                # Exclude segments overlapping with the testing segment\n",
    "                if abs(i - j) < win_size and i != j:\n",
    "                    correlations[i, j] = -np.inf\n",
    "\n",
    "        max_idx =  np.argmax(correlations, axis=1)\n",
    "        accuracy[test_subject] = sum(max_idx == range(n_segments)) / n_segments\n",
    "\n",
    "        # Print accuracy for each subject as we go\n",
    "        print(f\"{accuracy[test_subject]:.3f}\",\n",
    "              end=', ', flush=True)\n",
    "        \n",
    "    # Get a rough estimate of chance (accounting for excluded segments)\n",
    "    chance = 1 / np.sum(~np.isinf(correlations[n_trs // 2]))\n",
    "        \n",
    "    print(\"\\nThe average accuracy among all subjects is \"\n",
    "          f\"{np.mean(accuracy):.3f} +/- {np.std(accuracy):.3f}\")\n",
    "    return accuracy, chance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute time segment matching accuracy for the anatomically-aligned data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Between-subject time-segment classification accuracy for each subject: 0.094, 0.173, 0.090, 0.004, 0.158, 0.064, 0.064, 0.075, 0.143, 0.135, 0.132, 0.162, 0.083, 0.098, 0.053, 0.083, 0.169, 0.132, 0.030, 0.188, 0.075, 0.124, 0.083, 0.094, 0.135, 0.180, 0.195, 0.135, 0.102, 0.128, 0.109, 0.120, 0.060, 0.132, 0.041, 0.113, 0.132, 0.135, 0.165, 0.113, \n",
      "The average accuracy among all subjects is 0.113 +/- 0.044\n"
     ]
    }
   ],
   "source": [
    "# Time-segment classification on anatomically-aligned data\n",
    "win_size = 10\n",
    "acc_anat_test, chance = time_segment_classification(test_data, win_size=win_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we compute it after transforming the subjects data with SRM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Between-subject time-segment classification accuracy for each subject: 0.267, 0.462, 0.331, 0.038, 0.808, 0.188, 0.278, 0.462, 0.654, 0.876, 0.688, 0.553, 0.120, 0.327, 0.282, 0.440, 0.459, 0.312, 0.079, 0.523, 0.109, 0.583, 0.391, 0.436, 0.575, 0.511, 0.665, 0.308, 0.429, 0.320, 0.538, 0.571, 0.481, 0.372, 0.395, 0.229, 0.410, 0.451, 0.613, 0.526, \n",
      "The average accuracy among all subjects is 0.427 +/- 0.186\n"
     ]
    }
   ],
   "source": [
    "# Time-segment classification on SRM data\n",
    "acc_shared_test, chance = time_segment_classification(test_shared, win_size=win_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we plot the classification accuracies to compare methods. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEjCAYAAAAypHaFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXiNd/7/8efJIkpU0CRKMKi1EVpjSWhUEWsSO6VUFy2qGb52UhRRglHaqU7U1FRjqY6dohiUxFJqHWvVWlkUkY2c5Ny/P8T5NSWO7ZyQvh7XlevKfX/uc9/v+87JeZ17+9wmwzAMRETkT88prwsQEZHHgwJBREQABYKIiGRTIIiICKBAEBGRbAoEEREBFAh56vz581SrVo3Q0FBCQ0MJDg6mU6dO7Nmzx+Zrz507x/vvv++AKu1n+PDhzJkz545toaGhXLt27YHme+DAAUaPHn3HtsWLFxMdHQ3AggULiIqKeqBl5CfJycn07Nnzvl5zt7/dg9q4cSMTJkwA4MiRIzRt2pT27dvz1VdfWcc/iPDwcA4dOgTAqFGjiImJeST15kcueV3An13BggVZvny5dXjNmjWMGDGC9evX3/V1v/76K7/88ou9y8szv98m9+vkyZPEx8ffsW3Pnj1UqlQJgFdfffWBl5GfJCUlcfDgwbwugyZNmtCkSRPgZjjUq1ePiIiIh55vTEwMXbp0AXgk88vPFAiPmatXr+Lp6Wkd3rRpE7NmzcJsNlOwYEGGDRuGn58f4eHhxMfH89Zbb+Hm5kbjxo3p1KkTP/30E127dmXDhg2UKVOGzz77jNTUVIYMGcKsWbNYv349FouF0qVLM2bMGLy9vUlOTiYiIoLjx49jNpvx9/dn6NChuLi4UKNGDd555x22b99OQkICb7/9Nt26dbut7sTERIYNG8aVK1cAaNSoEQMGDGDJkiWsW7eOf/7znwC3De/Zs4d169aRkpJCgwYNGDZsGC4uLlSpUoXY2FiKFy/O4sWLWbBgARaLBQ8PDz744AMqVqxIamoqEyZMYO/evTg7O9O0aVNeffVVZs6cSXJyMiNGjOCjjz6y1vj999+zadMmtm/fTsGCBbl8+TJXrlxh9OjRvPLKK7Rp04YdO3aQlJTE22+/zd69ezl8+DAuLi7MmjULb29v4uPjGTduHBcvXsRsNtO6dWv69Olz2/bIzMxk/Pjx7N27F1dXV3x8fPjoo48oXLgwe/fuZerUqaSnp+Pk5ET//v1p3LgxWVlZREZGsmnTJooUKYKfnx8///wz8+bNo0ePHjz//PPs27ePy5cv07lzZy5dusSuXbtIT0/n448/pkqVKg/0txwxYgTXr18nNDSUJUuW4OzsbF2PO23jgQMH5ljXb7/9lkWLFmE2m0lKSqJ3795069Yt1/eErfdK69atWbBgAVlZWVy/fp0GDRpY3zOJiYmMGTOGU6dO4eTkRNeuXenZsyf79u1jypQpZGRkkJiYSEBAABMnTmT69OkkJCQwePBgIiMjmTp1Kt27d6dFixZs2LCBTz/9FIvFQuHChRkxYgR+fn588sknXLhwgcTERC5cuIC3tzdTpkzBy8vrgf6nnyiG5Jlz584ZVatWNUJCQoyQkBDj5ZdfNp5//nlj8+bNhmEYxi+//GK0adPGuHz5smEYhnH8+HGjQYMGRmpqqrFjxw6jdevWhmEYxtKlS43333/fMAzDmDFjhtGgQQNj4cKFhmEYRocOHYz9+/cbS5cuNQYMGGCYzWbDMAxj4cKFxttvv20YhmEMHz7c+OqrrwzDMIzMzExj8ODBRlRUlGEYhlG5cmVj3rx5hmEYxsGDBw1fX1/j+vXrt63Lp59+anzwwQeGYRhGamqqMWDAAOPatWvGf/7zH+Odd96xTvf74WHDhhnt2rUzUlNTjRs3bhivvfaaER0dbV3ub7/9ZuzcudPo1q2bkZaWZhiGYfzwww9GixYtDMMwjIkTJxoDBw40MjMzjRs3bhjdu3c3duzYcdsyf2/YsGHGF198YRiGYcycOdP48MMPDcMwjMaNGxsTJ040DMMwVq9ebVStWtU4cuSIYRiG0a9fP2PWrFmGYRhGjx49jI0bNxqGYRjXr183evToYaxevfq25ezevdto0aKFYbFYDMMwjMjISGPPnj3G1atXjaCgIOPcuXOGYRhGXFycERgYaFy4cMFYsGCB0b17d+P69evGjRs3jDfffNN47bXXDMMwjNdee83o37+/YRiGsW/fPqNy5crWOiIiIozw8PAH/lueO3fOqFWr1h23V27b+NZ2TElJMTp37mx9j/7000/WeeX2nriX98rv/za/H//ee+8ZkydPNgzDMK5du2a0bt3aOH36tDFw4EBjx44dhmEYRkpKilGvXj3j4MGD1r/tgQMHrNvxu+++M06ePGkEBAQYZ8+eNQzDMGJiYowGDRoYycnJxsyZM40mTZoYycnJhmEYxrvvvmvMmDHjjtsnv9EeQh774yGjmJgY3nvvPVasWGH9JterVy9ru8lk4uzZsznm0bhxYz766CMyMzPZtm0bffv2Zfv27bz88stcvnyZGjVqMGfOHA4ePEiHDh0AsFgspKenA7B582YOHjzIt99+C8D169dzzP/Wbvzzzz9PRkYGaWlpuLm55ZjmpZde4p133uHixYsEBAQwaNAgihQpYnP9Q0NDKVSoEAAhISFs2bIlxx7I5s2bOXPmDF27drWOu3btGlevXiUmJoYRI0bg7OyMs7MzX3/9NXBzL+RBBAUFAVCmTBmeeeYZqlatCkDZsmVJSkoiLS2N3bt3k5SUxIwZMwBIS0vj6NGjtGrVKse8KleujLOzM506daJhw4Y0b94cPz8/tmzZQmJiIu+99551WpPJxLFjx9iyZQuhoaHWbdulSxfmzZtnna5Zs2bW+uDmNr9V365du6zb637/lneT2zZeunQpAIULF+bzzz9ny5YtnD59mqNHj1rnmdt74kHfK7fqGTJkCABFihRh1apVAEyaNImtW7fy+eefc+rUKW7cuHHXdduxYwf169e3bkt/f3+KFy9uPddQt25d3N3dAahevTpJSUn3VN+TToHwmAkICKBs2bIcPHgQi8WCv78/H3/8sbX94sWLeHl58eOPP1rHFS1alGrVqvHf//6XlJQUQkND+eyzz9iwYQNNmzbFZDJhsVhyHO7JyMiwvsktFgszZsygYsWKwM0PXJPJZJ3/rQ+oW+MMw2DUqFHWf56uXbvy6quvsnHjRmJjY9mxYwedOnVi9uzZmEwmjN91l2U2m3Os7+8PTxiGgYtLzrekxWIhNDTU+iFgsVhISEigaNGiuLi45Kjz4sWLFCxY8L629+8VKFDA+rurq+tt7RaLBcMwWLhwIU899RQAly9fxs3NjQULFrBw4UIAfH19iYiIYPny5ezdu5cdO3YwYMAA3nrrLZ599lkqVqzI4sWLrfONj4+nePHitwWZk1POaz5+X9/darzfv+Xd2NrGcXFxdOnShc6dO1O7dm1atGjBf//7XwD8/Pzu+J7Ibfy9+GM9586do1ixYrz55ptUqVKFl156iZYtW7J///67rpvFYskxn1vbIjMzEyDHOv7xPZyf6Sqjx8wvv/zChQsXqFatGv7+/mzfvp2ff/4ZgC1bthASEsL169dxdnbO8eHarFkz/v73v+Pv74+7uzt/+ctfmD17tvVbb8OGDfn2229JSUkBYMaMGQwdOtTaNnfuXAzDICMjg759+1q/Cebm1gfe8uXLefXVV5k6dSqfffYZTZs2ZdSoUTz33HOcOHGC4sWLc+LECW7cuIHZbGbdunU55rN69WoyMjK4ceMGS5cuJTAwMEd7w4YNWb16NQkJCcDNK4Nef/114Oa3uqVLl2KxWMjIyCAsLIzdu3fj7Oxs/cf+o7u12eLu7k6tWrX48ssvgZsftreC8NVXX7Vuj4iICP773//Sq1cvXnjhBd5//33atm3LoUOHqFWrFmfOnGH37t3AzatpmjdvTnx8PI0aNWLFihVkZGSQmZlp/RZ+Px7kb+ni4kJWVtYdP/Ry28a3HDp0iOLFi9OvXz8aNmxoDYOsrKxc3xO5jb8X/v7+/Oc//wFuXh31+uuvc/r0aQ4ePMjgwYMJCgoiLi6Os2fPYrFYgDv/zf39/dm2bRvnzp0DIDY2losXL1KzZs17qiO/0h5CHrt1Mu8Wi8XCuHHjKF++PADjxo3j//7v/6zfnmfNmkXhwoV57rnncHNzo2PHjixevJimTZsyfvx4Bg8eDNz8YIiOjubFF18EoFOnTsTHx9O5c2dMJhPPPvsskyZNAm5eihcREUFwcDBms5mAgADefvvt+1qP119/neHDh9OmTRsKFChAlSpVaN26NU5OTtSpU4eWLVvi6elJvXr1OHbsmPV1Pj4+dOvWjdTUVJo1a0a7du1yzLdhw4b07t2bN998E5PJhLu7O59++ikmk4n+/fsTERFBaGgoWVlZtGrViqCgIM6cOcM//vEP+vfvz6effppjfoGBgdb1fhBTp05l/PjxBAcHk5GRQZs2bQgJCbltusDAQLZu3UqbNm0oVKgQRYsWZfz48RQvXpyZM2cSGRnJjRs3MAyDyMhIfHx8aN++Pb/88gtt27alUKFC+Pj4WPdE7tWD/C09PT3x8/OjdevWREdHU6xYMWtbbtt406ZNADRo0IBvv/2WFi1aYDKZqFu3LsWLF+fMmTO5vieSkpLuOP7W4Z+7GT16NGPHjiU4OBjDMHj33Xfx9fXlnXfeoV27dhQqVAhvb29efPFFzpw5g7+/P82aNWPIkCGMHTvWOp/nnnuOMWPG0L9/f7KysihYsCCff/75PR+6yq9Mxp9lX0ieGCkpKdSuXZv9+/c/1CGgJ822bdv47bffrF8QJkyYgJubm/VwmYi96ZCRPFYOHDhAixYtaN++/Z8qDAAqVarEsmXLCA4OpnXr1ly5cuWOl7SK2Iv2EEREBNAegoiIZFMgiIgI8IReZXT9+nUOHTqEp6dnjuvYRUQkd1lZWSQmJuLr63vHc3RPZCAcOnSI7t2753UZIiJPpOjoaP7617/eNv6JDIRbnb9FR0dTsmTJPK5GROTJEBcXR/fu3XN0oPl7T2Qg3DpMVLJkSXx8fPK4GhGRJ0tuh9p1UllERAAFgoiIZFMgiIgIoEAQEZFsCgQREQEUCCIikk2BICIigAJBRB4Tvr6+mEym+/rx9fXN67LzlSfyxjQRyX9uPaP7j/5MzzTOa9pDEBERQIEgIiLZFAgiIgIoEEREJJsCQUREAAWCiIhkUyCIiAigQBARkWwKBBERARQIIiKSTYEgIiKAnQNh5cqVtGrViqCgIKKjo29rP3z4MB06dCAkJIR3332Xa9eu2bMcERG5C7sFQnx8PNOnT2f+/PksW7aMRYsWcfLkyRzTREREEBYWxooVKyhfvjxz5syxVzkiImKD3QIhJiaG+vXr4+HhQaFChWjevDlr167NMY3FYiE1NRWA9PR0ChYsaK9yRETEBrsFQkJCAp6entZhLy8v4uPjc0wzfPhwwsPDadiwITExMXTt2tVe5YiIiA12CwSLxYLJZLIOG4aRY/j69euMGjWKuXPnsm3bNrp168awYcPsVY6IiNhgt0AoWbIkiYmJ1uHExES8vLysw8ePH8fNzQ0/Pz8AunTpwq5du+xVjoiI2GC3QAgICCA2NpbLly+Tnp7O+vXrCQwMtLaXK1eOuLg4Tp06BcDGjRupUaOGvcoREREb7PYITW9vbwYOHEjPnj0xm8107NgRPz8/evfuTVhYGDVq1OCjjz5iwIABGIZBiRIlmDhxor3KERERG+z6TOXg4GCCg4NzjJs9e7b190aNGtGoUSN7liAiIvdIdyqLiAigQBARkWwKBBERARQIIiKSTYEgIiKAAkFERLIpEEREBFAgiIhINgWCiIgACgQREcmmQBAREUCBICIi2RQIIiICKBBERCSbAkFERAAFgoiIZFMgiIgIoEAQEZFsCgQREQHuIRDatWvH4sWLSU9Pd0Q9IiKSR2wGQnh4OD/++CPNmjVj3LhxHD9+3BF1iYiIg7nYmqB27drUrl2ba9eusXLlSvr164eXlxc9evSgZcuWjqhRREQc4J7OIVy7do3ly5fzzTffUKRIEVq2bMny5csJDw+3d30iIuIgNvcQBg8ezObNm2ncuDFjx47lhRdeAODVV18lICCACRMm2L1IERGxP5uB8NxzzzFy5EiKFy+e84UuLixYsMBuhYmIiGPZPGT08ssvM3z4cACOHTtGaGgop06dAqBixYr2rU5ERBzGZiCMHTuWTp06AVClShXef/99xowZY/fCRETEsWwGQnp6Os2aNbMON23alJSUFLsWJSIijmczEEwmE0ePHrUO//zzzzg56QZnEZH8xuZJ5b/97W/06NGDypUrA3Dq1CmmTp1q98JERMSxbAZC48aNWbt2LXv37sXZ2ZmaNWtSokQJR9QmIiIOdE/HfuLj4ylWrBhFihThxIkTfPPNN/auS0REHMzmHkJ4eDgbN27kxo0beHl5cfbsWWrXrk3nzp0dUZ+IiDiIzT2EmJgYNm7cSLNmzYiKiuLLL7+kYMGCjqhNREQcyGYgeHp6UqhQISpUqMDx48epV68ecXFxjqhNREQcyGYguLq6snv3bipWrMjWrVtJTk4mLS3NEbWJiIgD2QyEIUOGsHDhQho1asTRo0epX78+ISEhjqhNREQcyOZJ5UOHDjFt2jQAvvnmG5KTkylSpMg9zXzlypXMmjWLzMxMXn/9dbp3756j/dSpU4wZM4akpCQ8PT35+9//TtGiRR9gNURE5GHZ3EP4Y4+m9xoG8fHxTJ8+nfnz57Ns2TIWLVrEyZMnre2GYdC3b1969+7NihUrqFatGlFRUfdZvoiIPCo29xDKly9PeHg4f/3rXylUqJB1fFBQ0F1fFxMTQ/369fHw8ACgefPmrF27lv79+wNw+PBhChUqRGBgIAB9+vTh2rVrD7wiIiLycGwGwtWrV7l69SpnzpyxjjOZTDYDISEhAU9PT+uwl5cXBw4csA6fPXuWZ555hpEjR3LkyBEqVKjABx988CDrICIij4DNQJg3b94DzdhisWAymazDhmHkGM7MzGTXrl18/fXX1KhRg48//phJkyYxadKkB1qeiIg8HJuBkNsjMm09T7lkyZL8+OOP1uHExES8vLysw56enpQrV44aNWoA0KZNG8LCwu6paBERefRsnlT28PCw/hQuXJhdu3bd04wDAgKIjY3l8uXLpKens379euv5AoAXXniBy5cvW7vW3rRpE88///wDroaIiDwsm3sIt04C39K7d2/69u1rc8be3t4MHDiQnj17Yjab6dixI35+fvTu3ZuwsDBq1KjBP/7xD8LDw0lPT6dkyZJERkY++JqIiMhDsRkIf+Tu7k5CQsI9TRscHExwcHCOcbNnz7b+XrNmTb799tv7LUFEROzgvs4hGIbB4cOHqVChgl2LEhERx7MZCLfuI7glJCREXVeIiORDNk8qv/vuu5QtW5b+/fvTpUsXMjIyctygJiIi+YPNQBg3bhybN2++ObGTE3v27GHixIn2rktERBzM5iGjn376iVWrVgFQokQJZsyYQWhoqN0LExERx7K5h2A2m8nIyLAOZ2Zm2rUgERHJGzb3EF5++WXeeustQkNDMZlMrFq1ikaNGjmiNhERcSCbgTB06FDmz5/Pxo0bcXFxISgoiC5dujiiNhERcSCbgWAYBkWLFmXWrFkkJiayevVqR9QlIvlUSZ+yxF84d1+v+X3HmPfCu3QZ4s6fva/XyD0EwtixY0lLSyMkJMR6ldH58+dtdm4nInIn8RfOYRq0wr7LmKZ7pR6EzUDYt2+frjISEfkT0FVGIiIC6CojERHJdk9XGUVHR1uvMmrWrBldu3Z1RG0iIuJANgPB2dmZnj170rNnT+u4tLQ09WckIpLP2AyEDRs2MHPmTNLS0jAMA4vFwtWrV/npp58cUZ+IiDiIzUCIjIxkwIABLFiwgN69e7NhwwYKFy7siNpERMSBbF5l9NRTT9GqVStq1aqFm5sbY8eOtfZ+KiIi+YfNQHBzcyMjI4OyZcty5MgRnJyc7vuuQRERefzZPGT0yiuv8M477zB58mS6dOnCnj17KFasmCNqExERB7IZCH369CEkJARvb28+++wzdu/eTZs2bRxRm4iIOJDNQAAoVaoUANWrV6d69ep2LUhERPKGzXMIIiLy56BAEBERQIEgIiLZ7ulO5YkTJ5KUlIRhGBiGgclkYu/evY6oT0REHMRmIEyZMoXhw4dTvXp13X8gIpKP2QyEp59+mqCgIEfUIiIiecjmOYSaNWuyZcsWR9QiIiJ5yOYewpYtW/j6669xdXXF1dVV5xBERPIpm4Ewd+5cB5QhIiJ5zWYglC5dmjVr1vDDDz9gNptp2LAhbdu2dURtIpIPfTfUhSrPdrDrMo4NvadOGOQPbG61OXPmsGLFCtq1a4dhGHz55Zf8+uuv9OvXzxH1iUg+0zIyE9OgFXZdhjEtBGOyXReRL9kMhGXLlrFgwQLc3d0B6NixI507d1YgiIjkM/d0p/KtMAAoUqQILi7aHRMRyW9sBkLp0qX597//jdlsxmw2M3fuXGvvpyIikn/YDIQPP/yQDRs2UKtWLWrVqsX69esZPXr0Pc185cqVtGrViqCgIKKjo3OdbvPmzbzyyiv3XrWIiDxyNo/9eHt7M2/ePNLT07FYLBQuXPieZhwfH8/06dNZsmQJBQoUoGvXrtSrV4/nnnsux3SXLl1i8mSd/RERyWu5BkJERASjRo2iT58+d2z//PPP7zrjmJgY6tevj4eHBwDNmzdn7dq19O/fP8d04eHh9O/fn2nTpt1v7SIi8gjlGgj+/v7AzQ/yB5GQkICnp6d12MvLiwMHDuSY5quvvqJ69erUrFnzgZYhIiKPTq6BcOuY/pkzZxgwYECOtgkTJtCuXbu7zthiseToHfVWlxe3HD9+nPXr1zN37lzi4uIeqHgREXl0cg2EmTNncu3aNdasWUNKSop1vNlsZtu2bYSHh991xiVLluTHH3+0DicmJuLl5WUdXrt2LYmJiXTo0AGz2UxCQgLdunVj/vz5D7M+IiLygHK9yqhmzZp4eHjg5OSEh4eH9adkyZJMnTrV5owDAgKIjY3l8uXLpKens379egIDA63tYWFhrFu3juXLlxMVFYWXl5fCQEQkD+W6h9CoUSMaNWpEYGAgfn5+9z1jb29vBg4cSM+ePTGbzXTs2BE/Pz969+5NWFgYNWrUeKjCRUTk0bqnB+RMmDCBtLQ0DMPAYrFw5swZFi5caHPmwcHBBAcH5xg3e/bs26bz8fFh06ZN91G2iIg8ajZvTBs0aBBms5mffvqJ0qVLc/LkSSpXruyI2kRExIFsBkJqaioffvghDRs2JDAwkC+//JJ9+/Y5ojYREXEgm4Fw68aycuXKceLECZ5++ukcl4+KiEj+YPMcQrly5YiIiKBdu3aMGjWKtLQ0MjMzHVGbiIg4kM09hLFjx/LXv/6V6tWr06lTJ3bu3Mm4ceMcUZuIiDiQzUDIysri9OnTwM1LUZ999lmqVatm77pERMTBbAbCiBEjuHr1KoD1/MEHH3xg98JERMSxbAbC6dOnGTZsGHDzaWkjR47kxIkTdi9MREQcy2YgZGZm5ujLKDU1FcMw7FqUiIg4ns2rjNq2bUunTp1o0aIFJpOJ77//nvbt2zuiNhERcSCbgfDuu+/y3HPPERsbi4uLC4MHD6ZRo0aOqE1ERBwo10BISUnB3d2dq1evUrt2bWrXrm1tu3r1qvWGNRERyR9yDYQePXqwdOlS6tevf8cH3Rw5csQhBYqIiGPkGggdO3YE4Pvvv6dMmTIOK0hERPJGrlcZzZs3D8MwCAsLc2Q9IiKSR3LdQyhfvjy1atUiMzOTF1980Tr+1iGjvXv3OqRAEclfvEuXIX5aiN2XIfcv10D4xz/+QVxcHL179yYqKsqRNYlIPhZ3/ux9TW8ymXTvk4PkGghpaWmUKlWKf/3rX7i5uTmyJhERyQM2rzJq1KjRbQmtq4xERPKfXANh6dKlABw9etRhxYiISN6x2ZfRpUuX2LhxIwBTp07l9ddfV0iIiORDNgNh+PDhnDt3jtjYWLZu3UpoaCgTJkxwRG0iIuJANgPh6tWr9OrVi61bt9KmTRvat29Penq6I2oTEREHshkIZrMZs9nMDz/8QEBAAOnp6aSlpTmiNhERcSCbgdCkSRP8/f0pVqwYvr6+dOrUiTZt2jiiNhERcSCb3V+HhYXRuXNnvL29gZsnlqtWrWr3wkRExLHu6Sqjw4cPYzKZmDJlCh999JGuMhIRyYfu6yqjH374QVcZiYjkU7rKSEREAF1lJCIi2XSVkYiIALrKSEREstkMhIyMDA4ePMiOHTsAyMrK4rvvvmPgwIF2L05ERBzHZiAMHDiQc+fOkZiYSPXq1dm/fz9169Z1RG0iIuJANs8hHDlyhCVLltCkSRNGjhzJggULSEpKckRtIiLiQDYDwcvLCxcXF/7yl79w/PhxKlWqRHJysiNqExERB7IZCIUKFWLlypVUrVqV7777jmPHjt3zZacrV66kVatWBAUFER0dfVv7hg0bCA0NJSQkhH79+mnPQ0QkD9kMhA8++IAjR47QoEEDnJyceO2113jzzTdtzjg+Pp7p06czf/58li1bxqJFizh58qS1PSUlhbFjxxIVFcWKFSuoUqUKn3zyycOtjYiIPDCbgVC+fHmGDh2KyWTi448/Zvfu3XTr1s3mjGNiYqhfvz4eHh4UKlSI5s2bs3btWmu72WxmzJgx1stZq1SpwsWLFx9iVURE5GHkepVRcHDwXV+4cuXKu7YnJCTg6elpHfby8uLAgQPW4WLFitGsWTMArl+/TlRUFD169LinokVE5NHLNRA++OCDh5qxxWLBZDJZhw3DyDF8S3JyMu+99x5Vq1alXbt2D7VMERF5cLkeMqpbty5169albNmyrFmzhrp16/LMM88wd+5cypcvb3PGJUuWJDEx0TqcmJiIl5dXjmkSEhLo1q0bVapUISIi4iFWQ+6Vr68vJpPpnn98fX3zumQRcZB76nBWm+sAABISSURBVP66QoUKAJQuXZq6desycuRImzMOCAggNjaWy5cvk56ezvr16wkMDLS2Z2Vl0adPH1q2bMmoUaPuuPcgj96hQ4cwDOO2H+CO4w8dOpTHFYuIo9i8U/nKlSv07NkTADc3N3r16sWyZctsztjb25uBAwfSs2dPzGYzHTt2xM/Pj969exMWFkZcXBz/+9//yMrKYt26dcDNb6/aUxARyRs2AyErK4v4+Hjr1UCXLl2yfqO0JTg4+LaT07NnzwagRo0aevKaiMhjxGYg9OrVi7Zt2/LSSy9hMpmIiYlh6NChjqhNREQcyGYgdOzYEV9fX3bs2IGzszNvvfUWlStXdkRtIiLiQDYDAaBq1ap6BoKISD5n8yojERH5c1AgiIgIoEAQEZFsCgQREQEUCCIikk2BICIigAJBRESyKRDyoZI+Ze+rR9NbHQve72tK+pTN4zUVkUfpnm5MkydL/IVzmAatsP9ypoXYfRki4jjaQxAREUCBICIi2XTIKB/6bqgLVZ7tYPflHBuqt49IfqL/6HyoZWSmQ84hGNNCMCbbfTEi4iA6ZCQiIoACQUREsikQREQEUCCIiEg2BYKIiAAKBBERyabLTvMh79JlHNKthHfpMnZfhog4jgIhH4o7f/a+X2MymTAMww7ViMiTQoeMREQEUCCIiEg2BYKIiAAKBBERyaZAEBERQIEgIiLZFAgiIgIoEEREJJsCQUQeC76+vphMptt+gDuON5lM+Pr65nHV+YvuVBaRx8KhQ4fyuoQ/Pe0hiIgIoEAQEZFsCgQREQHsHAgrV66kVatWBAUFER0dfVv7kSNHaN++Pc2bN2fUqFFkZmbasxwREbkLuwVCfHw806dPZ/78+SxbtoxFixZx8uTJHNMMGTKE0aNHs27dOgzD4JtvvrFXOSIiYoPdrjKKiYmhfv36eHh4ANC8eXPWrl1L//79Abhw4QLXr1+nVq1aALRv356ZM2fSrVu3e17GF198QZEiRazDtWvX5uWXXyYjI4NPPvnktun9/f0JCAggJSWFf/7zn7e1BwYGUqdOHS5fvsyXX355W3vTpk2pWbMmcXFxd9zjadWqFdWqVePcuXN3DLe2bdtSsWJFfv75Z5YtW3Zbe+fOnSlTpgxHjhxhzZo1t7V3796dkiVLsn//fjZs2HBb+xtvvEHx4sXZvXs3W7duva393Xffxd3dnZiYGGJjY3O0tWnThoyMDAoUKMDmzZvZs2fPba8fNGgQAOvXr+fgwYM52lxdXQkLCwNg9erVHD16NEd74cKF6dOnDwBLly7l1KlTOdo9PDx46623AFi0aBHnz5/P0e7l5UWPHj0AmDdvHgkJCTnafXx86NKlCwBz5szh6tWrOdorVKhAu3btAPj8889JTU3N0V61alVat24NwMyZMzGbzTnaa9SoQVBQEADTpk3jj/Tee/D3HsD777+v9x72f+998cUXt7X/nt0CISEhAU9PT+uwl5cXBw4cyLXd09OT+Ph4e5Uj2Xx9fTGbzVSuXPm2tkKFCpGVlUX16tWpUKECcPOfqU6dOo4uU0TygMmw02OyZs2axY0bNxgwYAAA33zzDYcOHWLcuHEA7Nmzh2nTpjF//nwATp8+TZ8+fVi7dq3NeZ8/f54mTZqwceNGfHx87FG+iEi+Y+uz027nEEqWLEliYqJ1ODExES8vr1zbL126lKNdREQcy26BEBAQQGxsLJcvXyY9PZ3169cTGBhobS9dujRubm7W44XLly/P0S4iIo5lt0Dw9vZm4MCB9OzZk7Zt29KmTRv8/Pzo3bu39aTQ1KlT+eijj2jRogVpaWn07NnTXuWIiIgNdu3LKDg4mODg4BzjZs+ebf29atWqfPvtt/YsQURE7pHuVBYREUCBICIi2RQIIiICPKHPQ8jKygIgLi4ujysREXly3PrMvPUZ+kdPZCDcun+he/fueVyJiMiTJzExkXLlyt023m53KtvT9evXOXToEJ6enjg7O+d1OSIiT4SsrCwSExPx9fWlYMGCt7U/kYEgIiKPnk4qi4gIoEAQEZFsCgQREQEUCCIikk2BICIigAJBRESyKRD+JGbOnMmPP/740PPZuHEjM2bMuO/X7dy50/pcWpG1a9fSvn17QkJCCA4Otj7rt0ePHjRr1ozQ0FBCQ0Np0qQJvXr14tKlS9b22rVrk5GRkWN+oaGhen89Ak/kncpy/3bv3k29evUeej5NmjShSZMmj6Ai+bOKj49n8uTJLFmyhGLFipGamkqPHj0oX748ABMmTLC+Vy0WC2FhYXz55ZcMGTIEAHd3d7Zt28Yrr7wCwKlTp0hISODpp5/OmxXKRxQIj7nMzEzGjh3LiRMnuHTpElWqVGHQoEEMGjSISpUqceTIEUqUKMGMGTPw8PDg66+/Zvny5aSnp+Pq6sq0adM4cOAAhw4dIjw8nE8//ZQCBQowevRorl69SqFChRg1ahR+fn4MHz6cp556iv/9739cu3aN//u//2P58uUcPXqUpk2bMnz4cJYsWcKuXbuYNGkSMTExTJo0CcMwKFWqFNOmTQNg5MiRxMfHk5CQgL+/PxEREXm8FeVxcuXKFcxmM9evXwegcOHCTJo0CTc3t9umTUtL48qVK/j5+VnHBQUFsW7dOmsgrFmzhubNm/Pzzz87ZgXyMR0yesz99NNPuLq6smjRIr7//nuSk5PZsmULR48e5Y033mDVqlU8/fTTrFy5kpSUFDZs2MC8efNYtWoVL7/8MtHR0bRt2xZfX18mTJhAlSpVGDJkCD169GDlypWMGDGCv/3tb9Zd8ISEBBYtWsQ777zDiBEj+PDDD1m2bBnffPMNycnJ1royMjIYPHgwkydPZuXKlVSuXJmlS5eyefNmqlWrxqJFi1i3bh27d+/m8OHDebX55DFUtWpVmjRpQtOmTenYsSNTpkzBYrFY+9YJDw8nJCSEhg0b0qVLFwICAujVq5f19YGBgezatQuz2QzA5s2bady4cV6sSr6jPYTHXJ06dfDw8CA6OppTp05x+vRp0tLSKFGiBNWrVwegUqVKJCUl4e7uzrRp01i9ejWnT5/mhx9+oFq1ajnml5qaytmzZwkKCgKgVq1aFC1alFOnTgFYn2tdqlQpKlWqRIkSJQDw8PAgKSnJOp9jx47h7e1tnf+gQYOsbQcOHGDu3LmcOnWKq1evkpaWZqetI0+qDz/8kH79+rFt2za2bdtG586dmTp1KvD/Dxnt3buXsLAwmjVrRoECBayvLVCgALVr1yYmJoZnn32WMmXK3LFfHrl/2kN4zG3cuJHBgwdTsGBB2rdvT506dShVqlSO3WuTyYRhGFy8eJEuXbqQnJxMYGAg7dq1449dVd2p6yrDMKzd4bq6ulrHu7jk/n3B1dUVk8lkHU5OTiYuLo558+YRGRlJ8eLFee2116hYseIdlyl/Xps3b2bNmjV4e3vToUMHpk+fTnh4+G2P033xxRfp0aMHgwYNIjMzM0dbixYtWLduHd999x2tWrVyZPn5mgLhMRcbG0vLli3p0KEDTz/9NDt37sy1L/ODBw9Srlw5evXqRY0aNdiwYYN1WmdnZ7KysnB3d8fHx4f169cDsG/fPi5dukSlSpXuq67y5cvz22+/cfLkSQC++OILFixYwPbt2+nSpQshISHcuHGDo0ePYrFYHmILSH5TsGBBpk2bxvnz54GbX0iOHDly294swBtvvEFqaiqLFi3KMT4wMJCdO3eydetW616tPDwdMnrMderUicGDB7N69WpcXV158cUX2blz5x2nbdCgAQsWLKBVq1YYhkGdOnU4ceIEAC+99BJjxoxh8uTJTJkyhbFjx/LJJ5/g6urKJ598kmOX/F64ubkxZcoUhg4ditlspmzZskRGRnLgwAHGjh1LVFQU7u7uvPDCC5w/f56yZcs+9LaQ/KF+/fr079+fPn36WM8DvPTSS7z33nu89dZbOaYtUKAAAwYMYOLEiYSEhOQY/+KLLwLc8WS0PBh1fy0iIoAOGYmISDYFgoiIAAoEERHJpkB4gp0/f956t+aMGTPYuHFjHlf0/507d46RI0fmdRkich90lVE+8be//S2vS8jh119/5dy5c3ldhojcBwXCEyK3Po1uGT58OHXr1qV9+/Z89dVXfP311xQpUoQKFSpQtmxZ3n//fRo2bEjz5s3Zs2cPzs7OfPzxx5QpU4ZXXnmF1q1bs337dlxcXOjXrx//+te/OHPmDMOGDaNVq1ZcunSJ0aNHExcXh8lkYtCgQQQEBPDJJ58QHx/PmTNnuHDhAp06daJv375MmDCB8+fP8+GHHzJmzJg83HIicq90yOgJkVufRn909OhRoqOjWbJkCfPnz+fMmTPWtsTERPz9/Vm2bBl16tQhOjra2vbMM8+wZMkSKlasSFRUFP/617+YMmUKUVFRAERERNChQweWLFnCrFmzGD16NCkpKcDNbizmzJnD4sWLiYqK4tq1a4SHh+Pr66swEHmCaA/hCZFbn0Z/FBsbS+PGjXF3dwegdevWXLt2zdr+0ksvATf7P/r98xF+34eRl5cXLi4ulCpVyvramJgYTp06xcyZM4Gbeyy3DgnVq1ePAgUKUKJECTw8PHJ0giciTw4FwhNi48aNzJw5k549e9K+fXuuXLlCqVKlbpvOycnprl1F3Lqr81b/R7fY6sPIYrHw73//Gw8PD+Bmr6glSpRgw4YNd+xXSUSePDpk9IS41z6N/P392bJlCykpKWRkZLB+/focndA9qPr16zN//nwATp48SXBwMOnp6blO7+zsfFuHZCLyeNMewhPiXvs0qly5Mj179qRLly4UKlSIYsWKPZK+XsLDwxk9ejTBwcEAREZGWg9L3UnFihVJTk5myJAhTJky5aGXLyL2p76M8plffvmFLVu2WB8o0rdvXzp16mS9X0FEJDfaQ8hnSpcuzcGDB2nTpg0mk4mGDRvqaVIick+0hyAiIoBOKouISDYFgoiIAAoEERHJpkAQ4WZfUHPmzAEgNDQ0x93dj5PNmzczY8aMvC5D8ildZSTyB8uXL8/rEnJ18OBBkpKS8roMyacUCPKnYrFYmDhxIvv37yc1NRXDMJgwYUKOaapUqUJsbCxFixYlMjKSTZs2UaRIEfz8/Pj555+ZN28ePXr0oFatWuzdu5eLFy/i7+/P+PHj+fXXX3n99ddp0KABhw4dIisri7CwMBYtWsSpU6fw9fXl73//O05OTuzdu5epU6eSnp6Ok5MT/fv3p3HjxixZsoTvv/8eJycnzpw5Q8GCBZk8eTIpKSksXLiQrKwsihQpwsCBA/NoK0p+pUCQP5X9+/eTkJDAokWLcHJyIioqitmzZ1v7aPq9xYsXc/jwYVatWoXJZKJv37452s+ePcu8efNIS0ujZcuW7Nq1Cx8fH86fP0+jRo0YN24cY8aMISIighUrVuDq6kqTJk3Yt28fFStWZMSIEcyZMwcfHx/i4+Pp3LkzVapUAWD37t2sWrWKkiVLMn78eKKiopg8eTJdu3blypUrCgOxCwWC/Km88MILFC1alIULF3Lu3Dl27txJ4cKF7xgIW7ZsITQ01Nr1R5cuXZg3b561vXHjxjg5OeHu7k65cuVISkrCx8cHV1dX653hZcuW5YUXXrB28+Hl5UVSUhL79u0jMTGR9957zzo/k8nEsWPHAHj++ecpWbIkANWrV+f777+3zwYR+R0FgvypbN68mYiICN544w2aNGlChQoVWLFixR2n/WOvr05OOa/BKFiwoPX33/fy6urqmqNDwd/3JHtLVlYWFStWZPHixdZx8fHxFC9enJUrV+Y6bxF70lVG8qeyfft2GjduTLdu3fD19WXDhg137DUWoFGjRqxYsYKMjAwyMzNZunTpI6ujVq1anDlzht27dwNw5MgRmjdvTnx8/F1fp15kxZ4UCPKn0rVrV3bt2kVwcDDt2rWjTJkynD9//o7PkGjfvj1+fn60bduWrl274urqylNPPfVI6ihevDgzZ84kMjKSkJAQhg4dSmRkJD4+Pnd9Xf369dm2bRvjx49/JHWI/J76MhLJxbZt2/jtt98IDQ0FYMKECbi5uTFkyJA8rkzEPhQIIrmIj49n+PDhXLp0CYvFQtWqVRk7dixFihTJ69JE7EKBICIigM4hiIhINgWCiIgACgQREcmmQBAREUCBICIi2RQIIiICwP8D1+m5bqhK7E8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SRM functional alignment provides a marked improvement in between-\n",
      "subject time-segment classification over anatomical alignment.\n",
      "The dotted line indicates chance performance (chance = 0.004)\n"
     ]
    }
   ],
   "source": [
    "# Box plot for the classification results\n",
    "labels = ['anatomical\\nalignment', 'SRM']\n",
    "\n",
    "plt.figure()\n",
    "plt.boxplot([acc_anat_test, acc_shared_test], vert=True,\n",
    "            patch_artist=True, labels=labels)\n",
    "plt.axhline(chance, linestyle='--', color='.4')\n",
    "plt.xlabel('alignment')\n",
    "plt.ylabel('classification accuracy')\n",
    "plt.title('Between-subject time-segment classification')\n",
    "plt.show()\n",
    "\n",
    "print(\"SRM functional alignment provides a marked improvement in \"\n",
    "      \"between-\\nsubject time-segment classification over \"\n",
    "      \"anatomical alignment.\\nThe dotted line indicates chance \"\n",
    "      f\"performance (chance = {chance:.3f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "The SRM allows us to find a reduced-dimension shared response spaces that resolves functional–topographical idiosyncrasies across subjects. We can use the resulting transformation matrices to project test data from any given subject into the shared space. The plot above shows the time segment matching accuracy for the training data, the test data without any transformation, and the test data when SRM is applied. The average performance without SRM is 11%, whereas with SRM is boosted to 40%. Projecting data into the shared space dramatically improves between-subject classification."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "* Chen, P. H. C., Chen, J., Yeshurun, Y., Hasson, U., Haxby, J., & Ramadge, P. J. (2015). A reduced-dimension fMRI shared response model. In C. Cortes, N.D. Lawrence, D.D. Lee, M. Sugiyama, R. Garnett (Eds.), *Advances in Neural Information Processing Systems, vol. 28* (pp. 460-468). https://papers.nips.cc/paper/5855-a-reduced-dimension-fmri-shared-response-model\n",
    "\n",
    "* Haxby, J. V., Guntupalli, J. S., Connolly, A. C., Halchenko, Y. O., Conroy, B. R., Gobbini, M. I., Hanke, M., & Ramadge, P. J. (2011). A common, high-dimensional model of the representational space in human ventral temporal cortex. *Neuron*, *72*(2), 404-416. https://doi.org/10.1016/j.neuron.2011.08.026\n",
    "\n",
    "* Haxby, J. V., Guntupalli, J. S., Nastase, S. A., & Feilong, M. (2020). Hyperalignment: Modeling shared information encoded in idiosyncratic cortical topographies. *eLife*, *9*, e56601. https://doi.org/10.7554/eLife.56601\n",
    "\n",
    "* Nastase, S. A., Liu, Y. F., Hillman, H., Norman, K. A., & Hasson, U. (2020). Leveraging shared connectivity to aggregate heterogeneous datasets into a common response space. *NeuroImage*, *217*, 116865. https://doi.org/10.1016/j.neuroimage.2020.116865\n",
    "\n",
    "* Anderson, M. J., Capota, M., Turek, J. S., Zhu, X., Willke, T. L., Wang, Y., Chen P.-H., Manning, J. R., Ramadge, P. J., & Norman, K. A. (2016). Enabling factor analysis on thousand-subject neuroimaging datasets.  *2016 IEEE International Conference on Big Data, pages 1151–1160*. http://ieeexplore.ieee.org/document/7840719/\n",
    " \n",
    "* Turek, J. S., Willke, T. L., Chen, P.-H., & Ramadge, P. J. (2017). A semi-supervised method for multi-subject fMRI functional alignment. *2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1098–1102*. https://ieeexplore.ieee.org/document/7952326\n",
    " \n",
    "* Turek, J. S., Ellis, C. T., Skalaban, L. J., Willke, T. L., & Turk-Browne, N. B. (2018). Capturing Shared and Individual Information in fMRI Data, *2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2018), pages 826-830*. https://ieeexplore.ieee.org/document/8462175\n",
    "\n",
    "* Richard, H., Martin, L., Pinho, A. L., Pillow, J., & Thirion, B. (2019). Fast shared response model for fMRI data. *arXiv:1909.12537*. https://arxiv.org/abs/1909.12537"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}