{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Author: Julia Jurkowska"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "ed8a8166dd34c963"
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Source localization and signal reconstruction - case study for oddball data"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "83bba7a26693beba"
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Introduction\n",
    "In this tutorial, we will learn how to **localize sources** from EEG data and **reconstruct signals** at those sources using **MVPURE-py**, an extension to MNE-Python.\n",
    "Source localization allows us to move beyond sensor-level analysis to estimate where in the brain the measured activity originates. Once sources are identified, we can reconstruct time series from vertices of interest for further analysis.\n",
    "\n",
    "We will cover the following steps:\n",
    "1. Reading all necessary data for the ``sample_subject``. You can download this dataset [here](https://figshare.com/articles/dataset/Sample_subject_data_/30102451?file=57853861).\n",
    "2. Computing data and noise covariance (R and N, respectively).\n",
    "3. Analysis of $RN^{-1}$ eigenvalues to guide the number of sources to localize and select an appropriate optimization parameter.\n",
    "4. Localizing the specified number of sources.\n",
    "5. Reconstructing source signals for vertices of interest and plotting the results.\n",
    "\n",
    "All steps will be repeated for two time frames: \"sensory\" (50-200 ms after stimuli) and \"cognitive\" (350-600 ms after stimuli).\n",
    "\n",
    "By the end of this tutorial, you will understand the basic workflow of source localization and signal reconstruction using the MVPURE-py package."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "e1aa86bb8a404230"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import mne\n",
    "import os\n",
    "\n",
    "mne.viz.set_3d_backend('pyvistaqt')\n",
    "\n",
    "from mvpure_py import localizer, beamformer, viz, utils"
   ],
   "metadata": {
    "collapsed": false,
    "is_executing": true
   },
   "id": "cccde4627289188e"
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will use data the ``sample_subject`` dataset [provided on Figshare](https://figshare.com/articles/dataset/Sample_subject_data_/30102451?file=57853861). If you wish to start from the beginning, please complete tutorial [Preprocessing data from oddball paradigm] first."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "1e3d629de1266e9f"
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading /Volumes/UMK/oddball/subjects/sample_subject/_eeg/_pre/sample_subject_oddball-epo.fif ...\n",
      "    Found the data of interest:\n",
      "        t =    -199.22 ...     800.78 ms\n",
      "        0 CTF compensation matrices available\n",
      "Not setting metadata\n",
      "621 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "EEG channel type selected for re-referencing\n",
      "Adding average EEG reference projection.\n",
      "1 projection items deactivated\n",
      "Average reference projection was added, but has not been applied yet. Use the apply_proj method to apply it.\n",
      "Created an SSP operator (subspace dimension = 1)\n",
      "1 projection items activated\n",
      "SSP projectors applied...\n"
     ]
    }
   ],
   "source": [
    "subject = \"sample_subject\"\n",
    "subjects_dir = \"subjects\"\n",
    "\n",
    "# Reading mne.Epochs\n",
    "epoched = mne.read_epochs(os.path.join(subjects_dir, subject, \"_eeg\", \"_pre\", f\"{subject}_oddball-epo.fif\"))\n",
    "forward_path = os.path.join(subjects_dir, subject, \"forward\", f\"{subject}_ico4-fwd.fif\")\n",
    "trans_path = os.path.join(subjects_dir, subject, \"_eeg\", \"trans\", f\"{subject}-fit_trans.fif\")\n",
    "\n",
    "# We will be using only data for 'target' stimuli\n",
    "target = epoched['target']\n",
    "sel_epoched = target.copy()\n",
    "sel_epoched = sel_epoched.set_eeg_reference('average', projection=True)\n",
    "sel_epoched.apply_proj()\n",
    "sel_evoked = sel_epoched.average()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:10:03.008966Z",
     "start_time": "2025-09-11T21:10:02.807604Z"
    }
   },
   "id": "4ea0b442e5c104c3"
  },
  {
   "cell_type": "markdown",
   "source": [
    "To perform source localization, we need a **forward model** that links activity at source locations to the sensors (in this case EEG channels).\n",
    "Here, we load the forward solution and convert it to a fixed-orientation representation."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "8b6a2f6dcde0c66e"
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading forward solution from /Volumes/UMK/oddball/subjects/sample_subject/forward/sample_subject_ico4-fwd.fif...\n",
      "    Reading a source space...\n",
      "    [done]\n",
      "    Reading a source space...\n",
      "    [done]\n",
      "    2 source spaces read\n",
      "    Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available\n",
      "    Read EEG forward solution (5124 sources, 128 channels, free orientations)\n",
      "    Source spaces transformed to the forward solution coordinate frame\n",
      "    No patch info available. The standard source space normals will be employed in the rotation to the local surface coordinates....\n",
      "    Changing to fixed-orientation forward solution with surface-based source orientations...\n",
      "    [done]\n"
     ]
    }
   ],
   "source": [
    "# Reading mne.Forward \n",
    "fwd_vector = mne.read_forward_solution(forward_path)\n",
    "\n",
    "# Using fixed orientation in forward solution\n",
    "fwd = mne.convert_forward_solution(\n",
    "    fwd_vector,\n",
    "    surf_ori=True,\n",
    "    force_fixed=True,\n",
    "    use_cps=True\n",
    ")\n",
    "   \n",
    "# Leadfield matrix\n",
    "leadfield = fwd[\"sol\"][\"data\"]\n",
    "\n",
    "# Source positions extracted from forward model\n",
    "src = fwd[\"src\"]"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:10:03.113362Z",
     "start_time": "2025-09-11T21:10:03.009744Z"
    }
   },
   "id": "924bcc8322c22de1"
  },
  {
   "cell_type": "markdown",
   "source": [
    "### \"Sensory\" processing\n",
    "We will start with analysing processes in \"sensory\" time window.\n",
    "\n",
    "In an oddball paradigm, participants are presented with a sequence of frequent (standard) and infrequent (target) stimuli. The early neural responses to these target stimuli reflect **sensory processing** - the brain's initial registration of the incoming stimulus before higher-level cognitive mechanisms are engaged.\n",
    "We assume that sensory processing for given oddball paradigm occurs within the **50-200 ms** window after the stimuli. We will therefore compute the data covariance in this time range.\n",
    "To estimate the noise covariance, we use a baseline period **-200ms to 0 ms**, i.e., the interval before stimulus onset. This baseline is assumed to be free of stimulus-locked activity and provides reference for separating signal from noise."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "c60a4a5cd016cc63"
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Created an SSP operator (subspace dimension = 1)\n",
      "    Setting small EEG eigenvalues to zero (without PCA)\n",
      "Reducing data rank from 128 -> 127\n",
      "Estimating covariance using EMPIRICAL\n",
      "Done.\n",
      "Number of samples used : 4056\n",
      "[done]\n",
      "    Created an SSP operator (subspace dimension = 1)\n",
      "    Setting small EEG eigenvalues to zero (without PCA)\n",
      "Reducing data rank from 128 -> 127\n",
      "Estimating covariance using EMPIRICAL\n",
      "Done.\n",
      "Number of samples used : 3042\n",
      "[done]\n"
     ]
    }
   ],
   "source": [
    "# Compute noise covariance\n",
    "noise_cov = mne.compute_covariance(\n",
    "    sel_epoched,\n",
    "    tmin=-0.2,\n",
    "    tmax=0,\n",
    "    method=\"empirical\"\n",
    ")\n",
    "\n",
    "# Compute data covariance for range corresponding to sensory processing\n",
    "data_cov_sen = mne.compute_covariance(\n",
    "    sel_epoched,\n",
    "    tmin=0.05,\n",
    "    tmax=0.2,\n",
    "    method=\"empirical\"\n",
    ")\n",
    "\n",
    "# Subset signal for given time range\n",
    "signal_sen = sel_evoked.crop(\n",
    "    tmin=0.05,\n",
    "    tmax=0.2\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:10:03.257319Z",
     "start_time": "2025-09-11T21:10:03.121816Z"
    }
   },
   "id": "d1c9366013218b82"
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### $RN^{-1}$ eigenvalues analysis\n",
    "Before attempting source localization, we need to decide **how many sources** to model and with what **rank**. Our proposition is to analyze the eigenvalues of the product of data covariance matrix $R$ and the inverse of the noise covariance matrix $N$. \n",
    "For a detailed theoretical background, see [PAPER]."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "79c997845cc13d8c"
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1200x600 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8gAAAITCAYAAADb3DVaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAATVpJREFUeJzt3Ql83HWdN/Bf0vQA2qaFgi09oICCpVBE0WVBDkGRxSJQUBGXIj4CS0FYHkHBVawulOPZ9SiIoAjrIoKgsNBdVFApWznkvspZkCKklCptWo6WNPO8vv84cZJm0iRNMjOZ9/v1Gqbzn+s3M9+E+eR31eRyuVwCAACAKldb6gYAAABAORCQAQAAQEAGAACAFgIyAAAACMgAAADQQkAGAAAAARkAAABaCMgAAAAgIAMAAEALARkAAAAEZAAAAGghIAMAFe+SSy5Ju+66axo8eHD6+te/XurmAFChBGQAoOKNGzcuC8YzZswodVMAqGACMkA3xBfwmpqatGzZsk5vd+WVV2a3++Mf/9hvbasWlfLedrVWKuk1lbNDDjkkHXzwwWnUqFG9/tgD/fMZ6K8PoDsEZGBAe/TRR9Phhx+ettpqqzRs2LA0fvz49OEPfzjNnTs3VbI777wzC2DLly8vdVOo8nrJh6v8qa6uLvs5O+aYY9JLL71U9Pbx89jR9fvss0+aOnVqr7cTALpCQAYGdCh43/velx5++OH0+c9/Pl100UXp//yf/5Nqa2vTd77znT597n/8x39Mb775ZhbM++q1zZ49W0AeIAZCvXzjG99I//mf/5m+//3vpwMPPDBdddVVae+9905vvfVWh7dfvXp1Ou+881Il6OvPB4DyUVfqBgD0lXPOOSfV19ene++9d51hl0uXLu3T5x40aFB2gmqplwjF8QepEH+IGjNmTDr//PPTTTfdlD7xiU+sc/tddtkl/eAHP0hnnnlm2nLLLYs+7p577pl+//vfd3jdV77ylfSv//qvqa8NhM8HgK7RgwwMWIsWLUo77rhjh3MSt9hii9Z/x1DQrbfeuugc0o7EvNL40j9y5Mi02WabpVNOOaVNT1mxOX0xpPTYY49N73jHO9LQoUOz9v3oRz9a5/Hjdp/73Oey4BC3mzx5cvqnf/qntGbNmqxdp59+ena7OJ4f2lr4XE8++WRavHjxet+jlStXplNPPTV7/fE88b7EEPQHHnggu/6FF15IJ554Ytp+++3TRhttlL3WI444Yp3XlX+vnn766fSZz3wm+8PE5ptvnr761a+mXC6XXnzxxfTxj388e7/Gjh2b/u3f/q3o+x1t7+y9Laar7+2GvB89qZf11Uo51Etf+OAHP9j6c9iRs846K61du3a9vcgLFizIaqijU2+E4668x8U+n9tvvz37o0AMF992223TpZdeWrQOuvI8+fs+++yzWZ3F7674WfrsZz+b3njjjew2119/fXab+fPnr/Mc8fxx3WOPPdatn9+OdLfOu/L6uvLzBVBqepCBASuGQ951113Zl8XentMYgSe+5M2ZMyfdfffd6bvf/W567bXX0o9//OOi93nllVfS3/3d32VfLk866aQsQN5yyy1ZsGlsbMy+OIaXX345vf/978+Gwx533HFphx12yL58xhfj+JJ82GGHZUH0pz/9afrWt76V9dSFeLy8d7/73dnw1vgC35kTTjghe9xoz5QpU9Kf//znLJA88cQT2ZY50fsew3M/9alPpQkTJmRfrGM7nZgnunDhwrTxxhu3ebxPfvKT2XNH6Pnv//7vLMBsuumm2Rf3D33oQ1mP4k9+8pP0xS9+Me22225pr7326tf3dn3W9370RE9eT3/XS1/Ih7DRo0d3eH2E9aOPPjrrRf7yl7/caS9yVzQ1NWWnCN1xHn+EiC2fOuv53ZC6efDBB9NHP/rRbPXsGL4ezxvDzDt6X7v7PFEz8f5EzUR4/OEPf5iFyfj5Oeigg9Lw4cPTz372s+xnvNC1116bBdP877vu/vz2VFdfX1/8fAH0uhzAAPXrX/86N2jQoOy0++67584444zcr371q9yaNWva3G7mzJm5rbbaap37n3322bn2vybzxw4++OA2x0888cTs+MMPP5xdvuKKK7LLzz//fOttPve5z+XGjRuXW7ZsWZv7fupTn8rV19fn3njjjezy0Ucfnautrc3de++967Spubk5O7/wwgvXefxCcd3ee++93vconnfWrFlFr8+3qdBdd92VPf6Pf/zjdd6X4447rvVYU1NTbsKECbmamprceeed13r8tddey2200UbZ+97f7+2Gvh/dqZeuvp5yqJcNkW/7bbfdlnv11VdzL774Yu7666/Pbb755rmhQ4dmlzu6fbR30aJFubq6utwXvvCF1uujbnfcccdutyP/fhee4rk609X3uKPPZ/r06bmNN94499JLL7Uee+aZZ7LX0/73RlefJ/8ajj322Da3O/TQQ3ObbbZZ6+Ujjzwyt8UWW2Q/Y3kNDQ1ZHXzjG9/o9s9vR6+vO78Xu/r6uvLzBVBqhlgDA1YM3Yse5Nj6JRbquuCCC9IBBxyQrbAb8yI3xKxZs9pcPvnkk7Pz//mf/+nw9pFZf/7zn6fp06dn/45ht/lTtGnFihVZT1Fzc3O68cYbs9vl53MWKjbku6PnW1/vcYghnPfcc0/WC9mRGJaZ9/bbb2c9Ptttt112v46GRcbc07zouYvXEG2JnqTC54whn88991y/vrddsb73oye6+3pKUS+9Yf/99896DidOnJitHL/JJptkP2fRc1nMNttsky2Addlll6WGhoYNev4Y+tt+CHYMEy5mQ+omeotvu+22bGupwp7v+NmIudgb+jzR09p+uHr87EVvbH6kRqyjUPgzHj2zUQ9xXU9/fnuiO6+vL36+AHqbgAwMaDGM9xe/+EU2pPUPf/hDtiBQzIOLL/AxxLCn3vnOd7a5HPMPY3XsYnP7Xn311WwIbASBCBGFp5hfGOILb9wuvgT31zY38UeDGIIeoSaG6UbIKAyusXLv1772tez6mDMYw3OjzfFa4otve5MmTWpzOeZPxvzM/LDewuPxmfTne9sb70dPdPf1lLpeYt7ykiVL2pwiEK7PxRdfnG699dYsqP3DP/xDFpCiZtbnX/7lX7Ih0f29ovWG1E0cj5+NCJvttT/Wk+dp/3OUH6ae/5mJod3xMxRDqvPi37Hw2bve9a4e//z2RHdeX1/8fAH0NnOQgaowZMiQLCzHKb5Axhe36667Lp199tlFe9m6Egq62lMXPTshFrCaOXNmh7fZeeedsx6Y/hRzHaN36oYbbki//vWv04UXXpjNc4w/KkRPWPR2XnHFFdkcwt133z37Uh6vNeY05l9ToY7mexabA9rV19pb721vvB+dtaer9dKVXt1S1kvMWd13333bHHv++ec7XLCpUASefC929KzG6tOf/vSn01NPPZXNme2sFzleZwSsmIvcX3qzbnr7edb3MxNhN97jqNPvfe972RzgWOn73HPPbXP77v78FupqnXfn9XXl5wug1ARkoOrkv8Tnh3RG70xH+8PGCrDFPPPMM9kiOnmx6mx8USwWIqI3ZcSIEdmXyxiKWkw8Rqx2nF+FtpjeHDobiwzFSrdxip6eWCwntsiKL6zRGxhfegtXnY7Fj/pyP92+em974/3oSb109/WUul6mTZuW9QQXipXHuyMCXiwwFUE79h9fX/CNXuTYNznCUn/ZkLqJBbNiZER8lu21P9bb9ZkXQ6n/4z/+I/3mN7/JFrmK8Fw4vDpsyM9vV+u8u69vfT9fAKVmiDUwYP3ud7/rsIctP/cz5sHmh7zGcMNHHnmk9TYRnqOXo7PhpIXmzp2bnRf7kheBYcaMGdlcvY7CTAxTDDH0NnqGbr755nTfffetc7v864n5naHYF92ubPMUX2jbD7OML/4xp3L16tWt7W7/HsZr7U7venf11Xu7Pl15P3pSL919PaWol/bBKIJO4SnCYHfFSsnRq/ztb397vdt0xXsaPZCx2nkM6e4PG1I3cd94X2L+d+F82gjHsXpzbz1PZ+L5Y4X4GFodp3ivC/8Qs6E/v12t866+vq7+fAGUmh5kYMCK4YWxzc2hhx6abX0Tcytj+Gh8mYzeu/z8uBhu+KUvfSm73Re+8IXsPrEVSgzFLraQTQw5jcW/Yi5gLAQWvV8xnDR634qJOZYR2j/wgQ+kz3/+89k2J3/5y1+y54gFf+LfIYZJxvDD2MIltu2JbZPii2kMCY8tUWKhm/e+973Zbb/yla9k7Y/tbGKRnHwQ6so2TzEXOxZQivnY0e4YBhvtiK1h8j1OH/vYx9J//ud/ZkMzo73xWuM2sZ9qX+nL97YzXXk/elIvPXk9/V0vfSX2X459d2Mf4fYLT7UXbYtaiyHZsVVRf9iQuon5s/G+77HHHtme0xEAo7c85oM/9NBDvfY8xcRnGFt4XXPNNen1119P/+///b91brMhP7/dqfOuvL6u/nwBlFypl9EG6Cu33HJLtl3KDjvskBs+fHhuyJAhue222y538skn51555ZV1toSaOnVqdpvtt98+d9VVV3W6zdPChQtzhx9+eG7EiBG50aNH50466aTcm2++2em2KSGeN7Y5mThxYm7w4MG5sWPH5vbbb7/cZZdd1uZ2L7zwQrZ9T36rnG222Sa73+rVq1tv881vfjM3fvz4bGuX9s/VlW2e4rFOP/303LRp07LXsckmm2T//t73vtdmS6bPfvazuTFjxmTv4QEHHJB78skns+1fCrdpyr8vsc1PobhNPG57HW3l01/v7Ya8H92pl66+nnKolw1RuG1Te2vXrs1tu+222Sm/JVFnt496iet6ss1TT3XlPS72+fzmN7/Jvec978nqIF7jD3/4w9z//b//Nzds2LAePU+xn6Niz3/rrbdmx2MrtfbbaXXn57fY43f192JXXl93fr4ASqkm/lPqkA4w0Fx++eXZlkcvvvhip9vc0LZHbvbs2dmQzParXg906mXgiCHvjz/+eDb3HIDKYw4yQB+IIa6xMFLMEYT1US+VKbZRKhShONY4iPnXAFQmc5ABelFstxIrx37/+9/PtlXZeOONS90k/rpA0PoWQ4o5kZ1tR9QX1Etliy2qjjnmmOw8VneOObqxpdwZZ5xR6qYB0EMCMkAviu1WYmGiWFH2Bz/4Qambw1/F0OX2K/y2F3tixzDv/qReKlssvPbTn/40W3k79iaOP3LEomnvfOc7S900AHrIHGQABrzYZihWdO5M9ALGCQCoXgIyAAAAlGKIdXNzc3r55ZfTiBEjsgVJAAAAoC9Fv3Dsyb7lllum2tra8gnIEY4nTpzY308LAABAlXtxPVsq9ntAjp7jfMNGjhyZSv1XhBUrVqT6+nq92agH2lAPtKcmKKQeKKQeKKQeylNjY2PWUZvPo2UTkPNFEuG4HAJynKIdihf1QCH1QHtqgkLqgULqgULqobyt7zMpPvgaAAAAqoiADAAAAAIyAAAAlGgOMgAAQDlYu3Ztevvtt3v1MWP+8Zo1a9Jbb71lDnI/Gjx4cBo0aNAGP46ADAAAVJUIsUuWLEnLly/vk8dvbm5Of/7zn/vksSlu1KhRaezYsRv0hwkBGQAAqCr5cLzFFlukjTfeuFd7eiN8R8909GbqQe4f8Z6/8cYbaenSpdnlcePG9fixBGQAAKBqRHjNh+PNNtus1x9fQC6NjTbaKDuPkByfbU+HW1ukCwAAqBr5OcfRc8zAsvFfP9MNmVcuIAMAAFVH7+7AU9MLn6mADAAAAAIyAAAAtBCQAQAAQEAGAACgI4ceemgaPXp0Ovzww1O1EJABAABYxymnnJJ+/OMfp2oiIAMAALCOffbZJ40YMSJVEwG5A4teXZXOuP7h9LG5/5udx2UAAIBSOuaYY7KtjOI0ePDgNHny5HTGGWekt956a53bnHfeeW3ue+ONN5Zsa6s77rgjTZ8+PW255ZZZG6It5aqu1A0oNxGGp89dkFY3Nae1zbn0RMPKNO+RhnTzyXumbTcfXurmAQAAVeyjH/1ouuKKK9Lbb7+d7r///jRz5swsdJ5//vmttxk2bFh2+fjjj8/mEBezyy67pKampnWO//rXv87CbG95/fXX07Rp09Kxxx6bDjvssFTOBOR2Lp2/qDUchziPy3H8gsOnlbp5AABAFRs6dGgaO3Zs9u+JEyem/fffP916661tAnIce/bZZ9OcOXPSBRdcUPSxHnrooX5p84EHHpidKoEh1u0sbGhsDcd5cXlhw8qStQkAAKC9xx57LN15551pyJAhbY4PGjQonXvuuWnu3LnpT3/6U68/77nnnpuGDx/e6Wnx4sWpEulBbmfKuJHZsOrCkDyotiZNGVddk9MBAKCaNL3xRq88Ti6XS2vXrk25QYM6nfNbt/HGPXr8efPmZQE0hkavXr061dbWposuuqjDLZpiCPXZZ5+dLr/88h491/77758efvjhbIj0hAkT0nXXXZd23333dMIJJ6RPfOITnd63N4do9ycBuZ3j9942m3OcH2Yd4XhoXW12HAAAGJhu2Wmnfn2+6YsW9eh+++67b7rkkkuy0Pqtb30r1dXVpRkzZnR42xh2/aEPfSh98Ytf7NFz3XbbbR0e33TTTbPTQGSIdTuxEFcsyDVj1/Fp6vj67NwCXQAAQDnYZJNN0nbbbZctevWjH/0o3XPPPUV7iPfaa690wAEHpDPPPLNX23CuIdbVJcKwBbkAAKB6HPjoo706xDrmAff1tkoxvPqss85Kp512Wvr0pz+dNtpoo3VuE9s9xVDr7bffvtee9wRDrAEAAAauns4J7igg1/RTQA5HHHFEOv3009PFF1/c4VDqnXbaKR111FHpu9/9bq8956bdHGK9atWqbFXtvOeffz5bQTseY9KkSamcGGINAABQoWIO8kknnZRt5xTzkjvyjW98IzU3N6dSue+++9J73vOe7BSixzv+/bWvfS2Vm5pc/ImjHzU2Nqb6+vq0YsWKNHLkyFRK8dKjHdGe/vjrDuVNPVBIPdCemqCQeqCQeqgsb731VtaDOXny5DRs2LBef/z+HGJN1z/bruZQPcgAAAAgIAMAAEAPAvLWW2+dDRNof5o1a1Z3HgYAAAAqexXre++9NxtPn/fYY4+lD3/4w9nKaQAAAFA1AXnzzTdfZ0+tbbfdNu2999693S4AAACojH2Q16xZk6666qpsie7OVmdbvXp1dipcPSy/uls/L6C9jnwbSt0OyoN6oJB6oD01QSH1QCH1UFnyn1N/fGZqonw+265+Fj0OyDfeeGNavnx5OuaYYzq93Zw5c9Ls2bPXOR7La5e6YOL5Y9PqYAl21AOF1APtqQkKqQcKqYfKEh19sSdwTB0tnD7am0q553A1W7t2bfber1y5sk0nbWFHbZ/tg3zAAQekIUOGpJtvvrnT23XUgzxx4sQsXNsHmXKiHiikHmhPTVBIPVBIPVTeXrl//OMfswWI+2If5JDfB5ny+Wwjh44aNWq9+yD3qAf5hRdeSLfddlv6xS9+sd7bDh06NDu1l18Bu9QKV+MG9UAh9UB7aoJC6oFC6qFy5D+jvvq8Cvsf1UP5fLZd/Sx6tA/yFVdckbbYYot00EEH9eTuAAAAUHa6HZBjTHcE5JkzZ6a6uh5PYQYAAIDKDsgxtHrx4sXp2GOP7ZsWAQAAUHKHHnpoGj16dDr88MNTteh2QP7IRz6Sjat/17ve1TctAgAAoOROOeWU9OMf/zhVkx7NQQYAAGBg22effdKIESNSNRGQAQAAKsAxxxzTukLz4MGD0+TJk9MZZ5yRbW/U/jbnnXdem/veeOONJVtV+5JLLkk777xztr1SnHbfffd0yy23dHjbz372s+lf/uVf0pw5c9Juu+2WBfRYIPqQQw5JTz31VJ+3VUAGAACoEB/96EdTQ0NDeu6559K3vvWtdOmll6azzz67zW1iD+Dzzz8/vfbaa50+1i677JKmTp26zunll1/u1TZPmDAhC+z3339/uu+++9KHPvSh9PGPfzw9/vjj6+wfPW/evHTwwQen+fPnp1mzZqW777473Xrrrentt9/Opvu+/vrrqS9ZhhoAAKBCDB06NI0dOzb798SJE9P++++fBcgIxHlx7Nlnn816YS+44IKij/XQQw/1S5unT5/e5vI555yT9SpH+N1xxx1bj995551Zz3j0HP/yl79sc58rr7wy60mOkL3XXnv1WVv1IAMAAFSgxx57LAuVQ4YMaXN80KBB6dxzz01z585Nf/rTn3r9ec8999w0fPjwTk+x81FHopf4mmuuyXqCY6h1oZtuuikL0x0NBV+xYkV2vummm6a+pAcZAACoem+vfbtXHid2/IkQ2JyaO53zO3jQ4B49fgxBjgDa1NSUVq9enWpra9NFF13U4RZNMYQ6hl9ffvnlPXqu/fffPz388MNZmI1h0tddd10Wak844YT0iU98otP7brnllm0uP/roo9l9Y750tP+GG25IU6ZMaXOb//qv/8qGjbfX3NycTj311LTHHntkQ8D7koAMAABUvYtvv7hfn+/U/U7t0f323XffbHhyhNYIk3V1dWnGjBkd3jaGXcd83y9+8Ys9eq7bbrutw+PRi9vdntztt98+G9IdPcHXX399mjlzZjbPOB+Sn3jiiWzu83777bfOfWMucvSWL1iwIPU1Q6wBAAAqxCabbJK22267NG3atPSjH/0o3XPPPUV7iGOu7gEHHJDOPPPMkg+xHjJkSNbu9773vdnc6Gj/d77znTbDqz/84Q9nC4wVOumkk7Je89/97ndZL3Zf04MMAABUvVn7zOrVIdYxD7ivt1WK4dVnnXVWOu2009KnP/3ptNFGG61zm1g9OoZaRw9ubzmhB0OsOxo2HUPEC4dXH3fccW3ex5NPPjkbin377bdnW1r1BwEZAACoej2dE9xeBLvaVNsvATkcccQR6fTTT08XX3xxh0Opd9ppp3TUUUel7373u732nJt2c4h19GAfeOCBadKkSWnlypXp6quvzkLvr371q+z6pUuXZts/RS9y4bDquF0E59gLecmSJdnx+vr6Dv8Q0FsMsQYAAKhQMQc5hiHHdk7F9gj+xje+kfXYlsrSpUvT0UcfnfVixxzje++9NwvHMaQ63Hzzzen9739/GjNmTOt9Yp51zFfeZ5990rhx41pP1157bZ+2VQ8yAABABYi9gDvy5S9/OTsVu83WW2/dZjhzf7t8PatoRy/xwQcfvE5PfCnoQQYAAKBk9txzz3TkkUemcqAHGQAAgJI544wzUrnQgwwAAAACMgAAALQQkAEAAEBABgAAgBYCMgAAAAjIAABANWpubi51EyjDz9Q2TwAAQNUYMmRIqq2tTS+//HLafPPNs8s1NTW99vi5XC6tXbs2DRo0qFcfl87f8zVr1qRXX301+2zjM+0pARkAAKgaEaAmT56cGhoaspDcVz2Z8Tz0r4033jhNmjRpg957ARkAAKgq0cMYQaqpqSnr7e3t3syVK1emESNG6EHuR9FjX1dXt8HvuYAMAABUnQhSgwcPzk69HZBXr16dhg0bJiBXIP3+AAAAICADAABACwEZAAAABGQAAABoISADAACAgAwAAAAtBGQAAAAQkAEAAKCFgAwAAAACMgAAALQQkAEAAEBABgAAgBYCMgAAAAjIAAAA0EJABgAAAAEZAAAAWgjIAAAAICADAABACwEZAAAABGQAAABoISADAACAgAwAAAAtBGQAAADoSUB+6aWX0mc+85m02WabpY022ijttNNO6b777uub1gEAAEA/qevOjV977bW0xx57pH333TfdcsstafPNN0/PPPNMGj16dN+1EAAAAMotIJ9//vlp4sSJ6Yorrmg9Nnny5L5oFwAAAJRvQL7pppvSAQcckI444og0f/78NH78+HTiiSemz3/+80Xvs3r16uyU19jYmJ3ncrnsVEr5NpS6HZQH9UAh9UB7aoJC6oFC6oFC6qE8dfXz6FZAfu6559Ill1ySTjvttHTWWWele++9N33hC19IQ4YMSTNnzuzwPnPmzEmzZ89e5/iKFStKXjTx/KtWrcr+XVNTU9K2UHrqgULqgfbUBIXUA4XUA4XUQ3nKd9SuT02uGyk1gvD73ve+dOedd7Yei4AcQfmuu+7qcg9yDNNevnx5GjlyZCqleOkR1Ovr6xUv6oE21APtqQkKqQcKqQcKqYfyFDl01KhR2WfTWQ7tVg/yuHHj0pQpU9oce/e7351+/vOfF73P0KFDs1N7USzlUDD5dpRDWyg99UAh9UB7aoJC6oFC6oFC6qH8dPWz6NY2T7GC9VNPPdXm2NNPP5222mqr7rUOAAAAyky3AvI///M/p7vvvjude+656dlnn01XX311uuyyy9KsWbP6roUAAABQbgF5t912SzfccEP66U9/mqZOnZq++c1vpm9/+9vpqKOO6rsWAgAAQD/o1hzk8LGPfSw7AQAAQNX2IAMAAMBAJSADAACAgAwAAAAtBGQAAAAQkAEAAKCFgAwAAAACMgAAALQQkAEAAEBABgAAgBYCMgAAAAjIAAAA0EJABgAAAAEZAAAAWgjIAAAAICADAABAi7q/ntMFi15dlS6dvygtbGhMU8aNTMfvvW3advPhpW4WAAAAvUBA7kY4nj53QVrd1JzWNufSEw0r07xHGtLNJ+8pJAMAAAwAhlh3UfQc58NxiPO4HMcBAACofHqQuyiGVefDcV5cXtiwMvu34dcAAACVTUDuogi9May6MCQPqq1JU8aNMPwaAABgADDEuouiR3hoXW0WikOcx+U4bvg1AABA5dOD3EXRExw9wi3DqFdmPcf5YdTrG34NAABA+ROQuyHC8AWHT+vW8GsAAAAqgyHWfTz8GgAAgMqgB7mPh18DAABQGQTkPh5+DQAAQGUwxBoAAAAEZAAAAGghIAMAAICADAAAAC0EZAAAABCQAQAAoIWADAAAAAIyAAAAtBCQAQAAQEAGAACAFgIyAAAACMgAAADQQkAGAAAAARkAAABaCMgAAAAgIAMAAEALARkAAAAEZAAAAGghIAMAAICADAAAAC0EZAAAAOhuQP7617+eampq2px22GGHvmsdAAAA9JO67t5hxx13TLfddtvfHqCu2w8BAAAAZafb6TYC8dixY/umNQAAAFApAfmZZ55JW265ZRo2bFjafffd05w5c9KkSZOK3n716tXZKa+xsTE7z+Vy2amU8m0odTsoD+qBQuqB9tQEhdQDhdQDhdRDeerq59GtgPyBD3wgXXnllWn77bdPDQ0Nafbs2emDH/xgeuyxx9KIESM6vE8E6LhdeytWrCh50cTzr1q1Kvt3zKemuqkHCqkH2lMTFFIPFFIPFFIP5SnfUbs+NbkNSKnLly9PW221Vfr3f//39LnPfa7LPcgTJ07M7jty5MhUSvHSI6jX19crXtQDbagH2lMTFFIPFFIPFFIP5Sly6KhRo7LPprMcukErbMUTvOtd70rPPvts0dsMHTo0O7WXXwW71ApX5Ab1QCH1QHtqgkLqgULqgULqofx09bPYoH2QY+jAokWL0rhx4zbkYQAAAKDkuhWQv/jFL6b58+enP/7xj+nOO+9Mhx56aBo0aFA68sgj+66FAAAA0A+6NcT6T3/6UxaG//znP6fNN9887bnnnunuu+/O/g0AAABVE5CvueaavmsJAAAAlNAGzUEGAACAgWKDVrGmaxa9uipdOn9RWtjQmKaMG5mO33vbtO3mw0vdLAAAAAoIyP0QjqfPXZBWNzWntc259ETDyjTvkYZ088l7CskAAABlxBDrPhY9x/lwHOI8LsdxAAAAyoeA3MdiWHU+HOfF5YUNK0vWJgAAANYlIPexmHM8qLamzbG4PGXciJK1CQAAgHUJyH0sFuQaWlfbGpLjPC7HcQAAAMqHRbr6WCzEFQtytaxivTLrObaKNQAAQPkRkPtBhOELDp9W6mYAAADQCUOsAQAAQEAGAACAFgIyAAAACMgAAADQQkAGAAAAARkAAABaCMgAAAAgIAMAAEALARkAAAAEZAAAAGghIAMAAICADAAAAC0EZAAAAEgp1ZW6AdVs0aur0qXzF6WFDY1pyriR6fi9t03bbj681M0CAACoSgJyCcPx9LkL0uqm5rS2OZeeaFiZ5j3SkG4+eU8hGQAAoAQMsS6R6DnOh+MQ53E5jgMAAND/BOQSiWHV+XCcF5cXNqwsWZsAAACqmYBcIjHneFBtTZtjcXnKuBElaxMAAEA1E5BLJBbkGlpX2xqS4zwux3EAAAD6n0W6SiQW4ooFuVpWsV6Z9RxbxRoAAKB0BOQSijB8weHTSt0MAAAADLEGAACAFgIyAAAACMgAAADQQkAGAAAAARkAAABaCMgAAAAgIAMAAEALARkAAAAEZAAAAGghIAMAAICADAAAAC0EZAAAAEgp1ZW6AXRs0aur0qXzF6WFDY1pyriR6fi9t03bbj681M0CAAAYsATkMg3H0+cuSKubmtPa5lx6omFlmvdIQ7r55D2FZAAAgD5iiHUZip7jfDgOcR6X4zgAAAB9Q0AuQzGsOh+O8+LywoaVJWsTAADAQCcgl6GYczyotqbNsbg8ZdyIkrUJAABgoNuggHzeeeelmpqadOqpp/Zei8gW5BpaV9sakuM8LsdxAAAAymyRrnvvvTddeumlaeedd+7dFpEtxBULcrWsYr0y6zm2ijUAAEAZBuRVq1alo446Kv3gBz9I//qv/9r7rSILwxccPq3UzQAAAKgaPQrIs2bNSgcddFDaf//91xuQV69enZ3yGhsbs/NcLpedSinfhlK3g/KgHiikHmhPTVBIPVBIPVBIPZSnrn4e3Q7I11xzTXrggQeyIdZdMWfOnDR79ux1jq9YsaLkRRPPH73hIeZSU93UA4XUA+2pCQqpBwqpBwqph/KU76jt1YD84osvplNOOSXdeuutadiwYV26z5lnnplOO+20Ng2bOHFiqq+vTyNHjkyllA/o0ZZKKt5Fr65Kl92xKD3R0JjePW5kOm4v85OruR7oG+qB9tQEhdQDhdQDhdRDeerqZ9GtgHz//fenpUuXpl133bX12Nq1a9Mdd9yRLrroomwo9aBBg9rcZ+jQodmpowaWQ8Hk21EObelqOD74ot+n1U3Nf90beVWa98iSbFEvIbn66oG+pR5oT01QSD1QSD1QSD2Un65+Ft3a5mm//fZLjz76aHrooYdaT+973/uyBbvi3+3DMb0vVrbOh+MQ53E5jgMAANBz3epBHjFiRJo6dWqbY5tssknabLPN1jlO31jY0NgajvNaepJXlqxNAAAAA0G3epApvSnjRqZBtW2HB8Tl2CsZAACAft7mqdDtt9++oQ9BNxy/97Zp3iMNrcOsIxwPravNjsf85BhqHb3MEaTjmHnJAAAA/RSQ6V8ReGNBrpYgvDLrOY4gHKbPXdAanJ9oWJkFaYt3AQAAdI2AXIEi8F5w+LQ2x864/uGii3e1vy0AAADrEpCrZPEuw68BAAA6JyAPEBF6Y1h1YUjOL94V4djwawAAgM5ZxXqAiB7hWKwrv8J14eJd9k4GAABYPz3IA3zxrjhu72QAAID1E5AH+OJd6xt+DQAAQAtDrKt8+DUAAAAt9CBX+fBrAAAAWgjIVT78GgAAgBaGWAMAAICADAAAAC0Msa5yi15d9de5yY3ZatfmJgMAANVKQK7ycDx97oK0uqk52wIqtoKa90hDtqCXkAwAAFQbQ6yrWPQc58NxiPO4HMcBAACqjYBcxWJYdT4c58Xl2AoKAACg2gjIVSzmHA+qrWlzLC7HPskAAADVRkCuYrEg19C62taQHOdxOY4DAABUG4t0VbFYiCsW5GpZxXpl1nNcuIq1Fa4BAIBqIiBXuQi8Fxw+bZ3jVrgGAACqjSHWdMgK1wAAQLURkOmQFa4BAIBqIyDTIStcAwAA1UZApkNWuAYAAKqNRbrokBWuAQCAaiMgU5QVrgEAgGoiINOrK1xHT7KeZQAAoBIJyPTaCtcPLl6uZxkAAKhYFumi11a4blrbbO9kAACgYgnI9NoK13WDau2dDAAAVCwBmR6vcD1j1/Fp6vj67Dwuv2fSKHsnAwAAFcscZHpthevoWY45x/lh1vZOBgAAKomATL/tnQwAAFDOBGT6Ze9kAACAcmcOMgAAAAjIAAAA0EJABgAAAAEZAAAAWliki36z6NVVf13hujFNGTfSCtcAAEBZEZDpt3A8fe6C1j2Sn2hYme2ZHNtCCckAAEA5MMSafhE9x/lwHOI8LsdxAACAciAg0y9iWHU+HOfF5YUNK0vWJgAAgEICMv0i5hwPqq1pcywuTxk3omRtAgAAKCQg0y9iQa6hdbWtITnO43IcBwAAKAcW6aJfxEJcsSBXyyrWK7OeY6tYAwAA5URApt9EGL7g8GnrHLf9EwAAUA4EZErK9k8AAEBFzkG+5JJL0s4775xGjhyZnXbfffd0yy239F3rGPBs/wQAAFRkD/KECRPSeeedl975znemXC6X/uM//iN9/OMfTw8++GDacccd+66VVO32T4ZfAwAAZRmQp0+f3ubyOeeck/Uq33333QIyPRKhN4ZVF4bk/PZPhl8DAAAVMQd57dq16brrrkuvv/56NtS6mNWrV2envMbGxuw8eqDjVEr5NpS6HdXsuL22Sf/9yMtpdVMuC8Et2z/VZMcvnf9sWtO0NjU351JsDhXna5py2fHzZ0zLAvRldyxKTzQ0pnePG5mO22vDepfVA4XUA+2pCQqpBwqpBwqph/LU1c+j2wH50UcfzQLxW2+9lYYPH55uuOGGNGXKlKK3nzNnTpo9e/Y6x1esWFHyoonnX7VqVfbvmpqW/XnpX2OGpHTtZ6eln9//p/T8sjfS5DEbpxnvnZDGDFmbli57LY3dqH2N5NLSZcvT439sSKde+1B6O4J1LpfuWbEyPfDsy+nbn9wlTRi9cY/aoh4opB5oT01QSD1QSD1QSD2Up3xH7frU5LqZUtesWZMWL16cBdzrr78+/fCHP0zz588vGpI76kGeOHFiWr58ebbQVynFS4/XUV9fr3jL0Jd+/nD6xQMvrTP8+rBdx2f/LnZd9C73hHqgkHqgPTVBIfVAIfVAIfVQniKHjho1KvtsOsuh3e5BHjJkSNpuu+2yf7/3ve9N9957b/rOd76TLr300g5vP3To0OzUXhRLORRMvh3l0BbaOn7v7dK8R5a0zkGOADykrjY7fso1D6am5rjV3z63uLywYdUGfZbqgULqgfbUBIXUA4XUA4XUQ/np6mexwfsgNzc3t+khht4S84ljQa6WVaxXZgt35Vex7mxxLwAAgJ7oVkA+88wz04EHHpgmTZqUVq5cma6++up0++23p1/96lc9enJYnwjDFxy+7pDpCMqxonVh7/LQrHd525K0EwAAqLKAvHTp0nT00UenhoaGbEz9zjvvnIXjD3/4w33XQuhm7zIAAECfB+TLL7+8R08C/dm7DAAA0BO1PboXAAAADDAbvEgXlJtFr67669DrxmwxL0OvAQCArhCQGXDhePrcBa2Ld8VK17GYV8xXjpAsPAMAAMUIyAwoEX7z4TjEeVyO4xGGOwvPAABAdTMHmQEleoYL90YOcTlWuu4sPAMAAAjIDCgxbDr2RC4Ul2MbqM7CMwAAgIDMgBLDqIfW1baG5DiPy3G8s/Acc5O/9POH06nXPJidx2UAAKC6mIPMgBJziWNOcctCXCuz8JtfiCvOY85xfph1Pjz/w07jsrnJa5rWprEb5dKCxW+keY8sMTcZAACqjIDMgBOh9oLDp3U5POfnJjd3MDe5o8cBAAAGJgGZVO3hOT83uXDwtbnJAABQfQRkql7MTY4tn/I9yIVzk4O9kwEAoDoIyFS9/NzkNU0RkFvmJg/568JeEY472ztZeAYAgIFDQKbq/W1u8rNp6bLl6e/HjErH771ddvyM6x8uundyhOHOwjMAAFBZBGT4a0g+f8a0tGLFilRfX59qalpmJHe2d3J+ca+OwrPFvQAAoPLYBxk60dneyZ2FZwAAoPIIyNCJGEYdeyXnQ3J+7+Q43ll4BgAAKo+ADF2Ynzxj1/Fp6vj67Dw/x7iz8AwAAFQec5ChB3snt13cK1axXpn1HFvFGgAAKpeADH0QngEAgMojIEMfsD8yAABUHgEZ+iAc2x8ZAAAqj0W6oJd1tj8yAABQvvQgQy9b3/7Ihl8DAEB5EpChl0XojWHVhSE5vz+y4dcAAFC+DLGGXtbZ/siGXwMAQPnSgwy9rLP9kdc3/BoAACgdARn6cX/kzoZfAwAApSUgQz+KnuSYc5wfZl04/NriXQAAUFoCMpTB8Otg8S4AACgtARnKYPj1Gdc/XHTxro6GagMAAL3PKtZQBizeBQAApScgQxmIOcf5baHyLN4FAAD9S0CGMt87GQAA6B/mIEOZ750crHANAAB9T0CGMt87OcKxFa4BAKDvGWINZS56joutcA0AAPQePchQ4StcG34NAAC9Q0CGMhehN4ZVF4bk/ArXhl8DAEDvMcQaKniFa8OvAQCg9+hBhgpe4bqz4deGXgMAQPcIyFDBK1wXG349YfQwQ68BAKCbDLGGATj8uibVdDr0OnqXz7j+4fSxuf+bncdlAACodnqQYQAOvz7lmgc7HXqtdxkAANYlIMMAHH7d2crXnS3s1dEwbgAAqBaGWEOVrXy9voW9DL0GAKBa6UGGKlv52sJeAADQMQEZqmzl6wjKEXzzQbgrC3sZeg0AQDXo1hDrOXPmpN122y2NGDEibbHFFumQQw5JTz31VN+1Duiz3uUZu45PU8fXZ+dx+cXX3ig69BoAAKpBt3qQ58+fn2bNmpWF5KampnTWWWelj3zkI2nhwoVpk0026btWAiVd2AsAAKpBtwLyL3/5yzaXr7zyyqwn+f7770977bVXb7cN6EfFhl7HcQAAqAYbNAd5xYoV2fmmm25a9DarV6/OTnmNjY3ZeS6Xy06llG9DqdtBeaj2ethmzCbpppP2SJfdsSjrSX73uBHpuL22zY7HexIrWrdc15jePW5kdt1AXryr2uuBdakJCqkHCqkHCqmH8tTVz6Mm18NPrrm5OR188MFp+fLlacGCBUVv9/Wvfz3Nnj17neMvvPBCGjlyZCqleOmrVq1Kw4cPTzU1LdvhUL3UQ3F/eu2NdOq1D6W3m3JpbS6XBtXUpMF1Nenbn9wlTRi9cRqI1APtqQkKqQcKqQcKqYfyFB21W221VdbJ21kO7XEPcsxFfuyxxzoNx+HMM89Mp512WpuGTZw4MdXX15dFQA7RFsWLeiju3Nv+mP7YmCuYnxxDsFP68f1L0/kzpg3I3mX1QHtqgkLqgULqgULqoTx19bPoUUA+6aST0rx589Idd9yRJkyY0Olthw4dmp06amA5FEy+HeXQFkpPPXQsVrJuao5//e19icsLG1al55a9ng6+6Petc5fj2LxHlgyI/ZPVA+2pCQqpBwqpBwqph/LT1c+itrt/DYlwfMMNN6Tf/va3afLkyT1tH1BBYoXrWLSrUH6F69gnudj+yQAAUElquzus+qqrrkpXX311thfykiVLstObb77Zdy0ESi5Wso4VrfMhuXCF64UNjfZPBgBgQOjWEOtLLrkkO99nn33aHL/iiivSMccc07stA8pGDJWOIdPRKxzBN3qOIxzH8c72T465yS33acxul79P6Ow6AAAo+4BsqXKoXhFeLzh8Wpf3T/6Hncal6XMXtB6PEB23i6Adil0nJAMAUJH7IAMU611e39zkYtfl76tnGQCA/iYgA33Su9z53ORch9c9uHi5nmUAACpjkS6A3lj5uth1TWubrYgNAEDJCMhAv698Xey6ukG1VsQGAKBkDLEG+n3l61Bs3vJzy17vcEXsYOVrAAD6koAM9PvK18WuK7YidhyPcGx+MgAAfckQa6Dsep1n7Do+TR1fn53nA/D6VsUGAIANpQcZqIhe585XxQYAgA2nBxmo+FWxAQCgN+hBBipCZ/OTO1vAy8JeAAB0lYAMVPyq2MUW8PreUbumE3/ygIW9AADoEgEZqPj5ycUW8Jp90+NFF/Yqtro2AADVS0AGKl6xBbwaGt/qdGEvw68BACgkIAMVL8JtDJ8uDMMxR3ncyGFp8WtvrnM8hmfbVxkAgPasYg1UvOj5jQW78qtc5xfwOvvgHTs8HrfvbF/lCM9f+vnD6dRrHszO4zIAAAOfHmRgQC/gVex4sWHZDy5envUsr2lam8ZulEsLFr+R5j2ypLVn2bBsAICBS0AGBvQCXsWOFxuW3bS2OetJbu6gZznCcGfDsoVnAIDKZog1UJWKDcuuG1RbdGGv9Q3LjvD88wdeSo+91Jidx2XDswEAKoeADFSl/PDrGbuOT1PH12fncfk9k0a1hub2C3sVG5a9vvAMAEBlMMQaqFodDb+OnuUYNr2mKYJuLgvHQwoW9upoWPb6wjMAAJVBQO7E22vfLnUT6Ee5XC77zONUU9O2B5HqMWnToemGWR9IP/rf59Krf16e9thsVDr2g9tkxz/3wUnpl4+9lFY35bLw2zIsuyY7/qP/fT4988qKdcLzjuM2Tk+98lp2/RNLGtO7x45Mx35wctpmjLnJlcbvCAqpBwqpB6q5HgYPGpwGkppcfIL9qLGxMdXX16cVK1akkSNHplKKlx7tiPZ0VLzf/s23S9IuAACASnDqfqemStDVHGoOMgAAABhi3blZ+8wqdRMooxEFVJfeqocZl/w+LXx53XnIU7YcmS48Yud0xCV3tS7ulV9J+7p/2j27jWHZ5cXvCAqpBwqpBwqph8omIFfReHrW/8ssPvM4+WVGb9XDDmNHp8dffmOduck7jB2VLv/fxen1NbGYVzx+TWpqTtnpW7cuSnc8vaw1OMf9b35kaet+y5SG3xEUUg8UUg8UUg+VzRBrgBLstxzHi618ffdzf7FlFABACehBBuiH/ZYj3MaWT7ElVITjOD5l3MgOt42KKN3ZllGLXl3118drzB4j/3gAAGwYARmgBPstF+653H4O8t9ts1n69cJXOtxvOcLx9LkLWu8TATsew/BrAIANJyADlFnvcpj/9KvrBOe4Lm5bbPh1hHC9ywAAPScgA5Rh73KxYdnF5i3H7TrrXQ6CMwBA5wRkgAoKzsXmLUeILta7fOGvnmyzKrZh2QAAHbOKNUCVr4odPc9nXP9w+tjc/83O4zIAQDXSgwxQxatiW/QLAOBvBGSAKl4V26JfAAB/IyADVPGq2Kdc86BFvwAA/kpABqjiVbH7YtEvvc4AQKUSkAGqQHeHZXfWu9zZol9xP3OaAYBKJSADVLHeXvSrsznN8bjFepb1OgMA5UBABqhyvbnoV7Gtph5cvLzT+cx6nQGAcmAfZAA67V2esev4NHV8fXYel794wPZF92KO3t/88by43LS2uWjPcme9zgAA/UkPMgC9tuhXsV7nukG1RYdlp5Tr5DrDrwGA/iMgA9BrwbnYnOa4/Nyy1zsclh2KraTd062mhGoAoCcEZAD6PDx3tlp2KHZdT7aaCuY0AwA9ISADUNLVskOx64ot+tXZVlOhJytpAwAIyACUdFh2Z9f1ZKupYnOaO1tJW0gGAIJVrAEoW9HD29GK2bHVVEerZUfvc09W0s7PWz7j+ofTx+b+b3YelwGA6qIHGYCKG5od5j/9arfmNHe2knZni4FFG1oW/Xo2LV32WtpizOh0/N7b6XUGgAGo2wH5jjvuSBdeeGG6//77U0NDQ7rhhhvSIYcc0jetA6DqdXerqWLXdbaSdmd7Mcd9IzyvaVqbxm6USwsWv5HmPbJkvStpAwBVEJBff/31NG3atHTsscemww47rG9aBQC9PKe5s5W0T7nmwaK9y/nw3NyNlbT/1ussPAPAgA7IBx54YHaqBk1vvFHqJtCPcrlcWvvmm6lp8OBUU9N2/iLVRz0MPFttUpv+6//smi5f8Hx6csmqtMPY4elze07Ojk/dbEha9OKadXqXp242OD25+NU0aM1bqS7l0uC3UxrSFMuA1aT7n2pIudVNaVCE7b/eJ9dUk35w6+PZ4x7+/bvSmr+G50UvLku/euCP6foTds9u19KGlWmHsSOy204eIzhXGr8jKKQeqOZ6qNt44zSQ1OTiE+zpnWtq1jvEevXq1dkpr7GxMU2cODEtX748jRw5MpVSvPQVK1ak+vr6Dot33nbblaRdAAAAleBjzz6bKkHk0FGjRmX5r7Mc2ueLdM2ZMyfNnj17nePRsA3I5r0inn/VqpZVSqvhrzsAAAC9acWKFalSAnJX6EHupAfZEOvqEvUQ9Rl16Q8mqAcKPb/s9fSj3y9Kry5bkTYfU5+O3aNltexPXHpn6zDqGJI9pK42/ez4v08/+v1z6aaHXl5nyPaIoXVp+Ztvr/P4U7asT9uPHd7hfQ7eZcvs38WuO3aPbbLne2rJyrT92BHZ5cljNunjdwS/IyikHqjmeqirkCHWZdODPHTo0OzUXhRLORRMvh0dtWXwJr5gVNsvs7qmpuxzL4fapLTUA4Xetckmac6kzdf5o+rP/3n/DlfS/vzGG6ebnlyemgrC86C62rTruzZPv174yjpB952TNk+PNTSmN2qHplTb9rkf+3NTVGSH192/ZHW66fIHWhcKe2TZX7LnXd9CYRYQ23B+R1BIPVBIPZSnrn4W9kEGgF5eSbsn+zfHbWMl7I62oQodXde0tnm921O1X2X7e0ftmk78yQNFV98GgGrW7YAcc3afLZiI/fzzz6eHHnoobbrppmnSpEm93T4AqIr9mzvbhip0dF3doNr1bk/VPjzPvunxoqE62qt3GYBq1u2AfN9996V999239fJpp52Wnc+cOTNdeeWVvds6AKjyXud8OO3ourj83LLXO+x1joDbUXhuaHyraKiOcNxRr7O9nQGoFt0OyPvss0/JV58GgGoKz8Wu66zXudiQ7XEjh6XFr73ZYagu1uvc2ZDtCM9BcAZgIGi33AcAUCnyvc4zdh2fpo6vz87zvb0RUiMsR/gN+fB89sE7dng8bl+s17mzIdsX/urJLDj//IGX0mMvNWbncTl6mwGg0likCwCqbMh2sePR+1tsobBi4fnu5/5iTjMAA4aADABVGJ47Ot6TIdvRD92TOc1BcAag3AjIAECms97lYuH577bZrMO9nTub0xzDsu94epnFwAAoOwIyANAnezufcs2D3R6W3dliYMIzAH1NQAYA+mRv52Jzmjsblt3TlbSFZAB6g4AMAGyQ7s5p7mxYdk9W0s6HZz3LAGwoARkA6BM9GZZdbDGwzsLzg4uX61kGoFcIyABA2QzL7slK2k1rm3u01ZT5zAC0JyADABWxf3Ox8Fw3qLbbW01976hd04k/eUCvMwBtCMgAQEWH57j83LLXu7XV1OybHu9Rr3Noue7ZtHTZa2mLMaPT8XtvJ1QDDBACMgBQ0eG5s2HZxbaaamh8q9u9zhHOQ1y3pmltGrtRLi1Y/Eaa98iS1usM2QaobAIyAFDROhuWXWyrqXEjh6XFr73ZrV7nOB7i383trrvwV0+mO55eZv9mgAonIAMAVbfV1NkH79hmDnJXep0jfKeUy/5d0+66u5/7i/2bAQYAARkAqMre5e72OsdtQlzX3O66CMz2bwaofAIyAFC1i351p9c5v4dzXLemKcJuy3VD6mrT322zWfr1wld6df9mw7IB+p+ADABQoLPe5dByXaxivTz9/ZhR2SrWYf7Tr/ba/s2GZQOUhoAMANDF3uX8defPmJZWrFiR6uvrU01Ny4zk3ty/ubNh2evbhgqAnhOQAQDKbP/mYsOy17cNlZAMsGFqN/D+AAB0MTzPO3nP7Dzfsxw9yRGKQ+Gw7OgVzh/v7jZUAPScgAwAUAL5nuUZu45PU8fXZ+f5XuDOwnNnvcsAbBhDrAEAKmRY9vq2oTI3GWDDCMgAAGWou9tQ/cNO48xNBthAhlgDAAyAodn/8+jfQnMwNxmg+/QgAwAMgN5lc5MBNpweZACAAaCzla8B6BoBGQBgAOhs5WsAusYQawCAAaCzla8B6BoBGQBggK98DUDXGGINAAAAAjIAAAC0EJABAABAQAYAAIAWAjIAAAAIyAAAANBCQAYAAAABGQAAAFoIyAAAACAgAwAAQAsBGQAAAARkAAAAaCEgAwAAgIAMAAAALQRkAAAAEJABAACghYAMAAAAAjIAAABsQEC++OKL09Zbb52GDRuWPvCBD6Q//OEPPXkYAAAAKBt13b3Dtddem0477bT0/e9/PwvH3/72t9MBBxyQnnrqqbTFFlv0TSsBANggi15dlS6dvygtbGhMU8aNTMfvvW3advPhnV7XX/ep/DY8m5Yuey1tMWZ0On7v7ar4fdCGcq6Hvni8gagml8vlunOHCMW77bZbuuiii7LLzc3NaeLEienkk09OX/7yl9d7/8bGxlRfX59WrFiRRo4cmUopXnq0I9pTU1NT0rZQeuqBQuqB9tQElVwP8QV3+twFaXVTc1rbnEuDamvS0LradPPJe2bXd3Td947aNZ34kwf6/D4DoQ1rmtamsRvl0pI3a9KQukFV+z5oQ/nWQ1883rYVFpK7mkO71YO8Zs2adP/996czzzyz9VhtbW3af//901133dXhfVavXp2dChuW/x9LN7N5r8u3odTtoDyoBwqpB9pTE1RyPUSPVnxpb27OpYjzcb6mKZcdDx1d942bHuuX+wyMNjQXXJeq+H3QhnKth754vPNnTEuVpKu/r7sVkJctW5bWrl2b3vGOd7Q5HpeffPLJDu8zZ86cNHv27HWOR3Iv9f9U4vlXrVqV/bsS/vpL31IPFFIPtKcmqOR6iOGe0aPVVi4tXbY8O+/oupqmN/vlPgOhDVEBW2yUv6Z63wdtKN966IvHW7FiRaok+Y7aXp+D3F3R2xxzlgsbFkOyo3u7HIZYh0oZHkXfUg8UUg+0pyao5HqIuZALFr+RDZHMi6GSfz9mVPbvjq6bNHqjtGTlm31+n4HQhuhVCy+9HqMrq/d90IbyrYe+eLz6+vpUSbr6u7pbq1iPGTMmDRo0KL3yyittjsflsWPHdnifoUOHZkG48JRvoJOTk5OTk5OTU9+fYqGgmAsZU+OifyvO43IcL3bd1w6e2i/30QZt0IbKfLyaMvjd1t1Tny3S9f73vz/NnTs3uxxj7CdNmpROOukki3RR0dQDhdQD7akJKr0e/rYS7co0ZdyIIqvUtr2uv+5T+W2IVYuXpy3GjOpg1eJqeh+0oZzroS8er5J0NYd2OyDHNk8zZ85Ml156aRaUY5unn/3sZ9kc5PZzkzekYf2hEv/nRt9RDxRSD7SnJiikHiikHiikHqpoFevwyU9+Mr366qvpa1/7WlqyZEnaZZdd0i9/+csuhWMAAAAoVz1apCuGU8cJAAAABopuLdIFAAAAA5WADAAAAAIyAAAAtBCQAQAAQEAGAACAFgIyAAAACMgAAADQQkAGAAAAARkAAABaCMgAAAAgIAMAAEALARkAAAAEZAAAAGhRl/pZLpfLzhsbG/v7qTtsS7SjpqYmO1Hd1AOF1APtqQkKqQcKqQcKqYfylM+f+TxaNgF55cqV2fnEiRP7+6kBAACoYitXrkz19fVFr6/JrS9C97Lm5ub08ssvpxEjRpT8LyrxV4QI6i+++GIaOXJkSdtC6akHCqkH2lMTFFIPFFIPFFIP5Slib4TjLbfcMtXW1pZPD3I0ZsKECamcROEqXvLUA4XUA+2pCQqpBwqpBwqph/LTWc9xnkW6AAAAQEAGAACAFlUdkIcOHZrOPvvs7BzUA4XUA+2pCQqpBwqpBwqph8rW74t0AQAAQDmq6h5kAAAAyBOQAQAAQEAGAACAFgIyAAAAVHtAvvjii9PWW2+dhg0blj7wgQ+kP/zhD6VuEv1gzpw5abfddksjRoxIW2yxRTrkkEPSU0891eY2b731Vpo1a1babLPN0vDhw9OMGTPSK6+8UrI20z/OO++8VFNTk0499dTWY2qh+rz00kvpM5/5TPaZb7TRRmmnnXZK9913X+v1sbbl1772tTRu3Ljs+v333z8988wzJW0zfWPt2rXpq1/9apo8eXL2WW+77bbpm9/8ZlYDeeph4LrjjjvS9OnT05Zbbpn9v+HGG29sc31XPvu//OUv6aijjkojR45Mo0aNSp/73OfSqlWr+vmV0Nf18Pbbb6cvfelL2f8vNtlkk+w2Rx99dHr55ZfbPIZ6qAxVG5CvvfbadNppp2VLsD/wwANp2rRp6YADDkhLly4tddPoY/Pnz88Cz913351uvfXW7JfaRz7ykfT666+33uaf//mf080335yuu+667PbxC+6www4rabvpW/fee2+69NJL084779zmuFqoLq+99lraY4890uDBg9Mtt9ySFi5cmP7t3/4tjR49uvU2F1xwQfrud7+bvv/976d77rkn+zIU//+IP6YwsJx//vnpkksuSRdddFF64oknssvx+c+dO7f1Nuph4IrvBfH9MDpUOtKVzz7C0OOPP55935g3b14Wso477rh+fBX0Rz288cYbWZ6IP6jF+S9+8Yus8+Xggw9uczv1UCFyVer9739/btasWa2X165dm9tyyy1zc+bMKWm76H9Lly6NroDc/Pnzs8vLly/PDR48OHfddde13uaJJ57IbnPXXXeVsKX0lZUrV+be+c535m699dbc3nvvnTvllFOy42qh+nzpS1/K7bnnnkWvb25uzo0dOzZ34YUXth6LOhk6dGjupz/9aT+1kv5y0EEH5Y499tg2xw477LDcUUcdlf1bPVSP+L1/ww03tF7uyme/cOHC7H733ntv621uueWWXE1NTe6ll17q51dAX9ZDR/7whz9kt3vhhReyy+qhclRlD/KaNWvS/fffnw2Fyautrc0u33XXXSVtG/1vxYoV2fmmm26anUdtRK9yYX3ssMMOadKkSepjgIoRBQcddFCbzzyohepz0003pfe9733piCOOyKZgvOc970k/+MEPWq9//vnn05IlS9rURH19fTZNR00MPH//93+ffvOb36Snn346u/zwww+nBQsWpAMPPDC7rB6qV1c++ziPYbTxOyUvbh/fOaPHmYH//TKGYkcNBPVQOepSFVq2bFk2r+gd73hHm+Nx+cknnyxZu+h/zc3N2XzTGFI5derU7Fj8D2/IkCGtv9AK6yOuY2C55pprsuFQMcS6PbVQfZ577rlsSG1MwTnrrLOyuvjCF76Q1cHMmTNbP/eO/v+hJgaeL3/5y6mxsTH7w9igQYOy7w7nnHNONkwyqIfq1ZXPPs7jD22F6urqsj/Iq4+BLYbZx5zkI488MptvHNRD5ajKgAyFPYePPfZY1iNA9XnxxRfTKaecks0FisX6IP5oFn/dP/fcc7PL0YMcvyNijmEEZKrLz372s/STn/wkXX311WnHHXdMDz30UPZH1ViARz0AHYmRZ5/4xCeyRdziD65UnqocYj1mzJjsL8HtV6KNy2PHji1Zu+hfJ510UrZAwu9+97s0YcKE1uNRAzEMf/ny5W1urz4GnhhCHQvz7brrrtlfceMUC3HFoivx7+gJUAvVJVajnTJlSptj7373u9PixYuzf+c/d///qA6nn3561ov8qU99Klud9h//8R+zhftiN4SgHqpXVz77OG+/+GtTU1O2krH6GNjh+IUXXsj++J7vPQ7qoXJUZUCOoXLvfe97s3lFhb0GcXn33Xcvadvoe/EXvQjHN9xwQ/rtb3+bbd9RKGojVrAtrI9YiTC+IKuPgWW//fZLjz76aNYrlD9F72EMn8z/Wy1Ul5hu0X7bt5h/utVWW2X/jt8X8UWmsCZiCG7MH1MTA0+sTBvzAwvFH9jjO0NQD9WrK599nMcfWOOPsXnxvSPqJ+YqMzDDcWz1ddttt2VbBRZSDxWk1KuElco111yTrTR45ZVXZqvKHXfccblRo0bllixZUuqm0cf+6Z/+KVdfX5+7/fbbcw0NDa2nN954o/U2J5xwQm7SpEm53/72t7n77rsvt/vuu2cnBr7CVayDWqgusepoXV1d7pxzzsk988wzuZ/85Ce5jTfeOHfVVVe13ua8887L/n/xX//1X7lHHnkk9/GPfzw3efLk3JtvvlnSttP7Zs6cmRs/fnxu3rx5ueeffz73i1/8IjdmzJjcGWec0Xob9TCwdzh48MEHs1N8Zf73f//37N/5VYm78tl/9KMfzb3nPe/J3XPPPbkFCxZkOyYceeSRJXxV9EU9rFmzJnfwwQfnJkyYkHvooYfafL9cvXp162Ooh8pQtQE5zJ07N/viO2TIkGzbp7vvvrvUTaIfxC+1jk5XXHFF623if24nnnhibvTo0dmX40MPPTT7JUf1BWS1UH1uvvnm3NSpU7M/ou6www65yy67rM31sb3LV7/61dw73vGO7Db77bdf7qmnnipZe+k7jY2N2e+D+K4wbNiw3DbbbJP7yle+0uYLr3oYuH73u991+H0h/nDS1c/+z3/+cxaAhg8fnhs5cmTus5/9bBa0GFj1EH9AK/b9Mu6Xpx4qQ038p9S92AAAAFBqVTkHGQAAANoTkAEAAEBABgAAgBYCMgAAAAjIAAAA0EJABgAAAAEZAAAAWgjIAAAAICADAABACwEZAAAABGQAAABoISADAACQSOn/A7VxAc38qXmbAAAAAElFTkSuQmCC"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Suggested number of sources to localize: 62\n",
      "Suggested rank is: 42\n"
     ]
    }
   ],
   "source": [
    "sugg_n_sources, sugg_rank = localizer.suggest_n_sources_and_rank(\n",
    "    R=data_cov_sen.data,\n",
    "    N=noise_cov.data,\n",
    "    show_plot=True,\n",
    "    subject=subject,\n",
    "    s=14\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:10:03.407036Z",
     "start_time": "2025-09-11T21:10:03.257192Z"
    }
   },
   "id": "37881de806459c46"
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Localize\n",
    "Based on the eigenvalue spectrium above, we will localize **62 sources** using **rank of 42**.\n",
    "We will use function ```mvpure_py.localizer.localize```, which performs the actual source localization.\n",
    "The main parameters are: \n",
    "- ``subject``: the subject ID (here: ``\"sample_subject\"``)\n",
    "- ``subjects_dir``: directory containing the subject folders.\n",
    "- ``localizer_to_use``: the algorithm variant. Here we choose ``\"mpz_mvp\"`` because it provides the highest spacial resolution. Other possible options include: ``\"mai\"``, ``\"mpz\"``, and ``\"mai_mvp\"``. For details, see [PAPER] or the function documentation.\n",
    "- ``n_sources_to_localize``: number of sources to localize. We will use the suggested number of sources from $RN^{-1}$ analysis.\n",
    "- ``R``: data covariance matrix\n",
    "- ``N``: noise covariance matrix\n",
    "- ```forward```: the ``mne.Forward`` object for this subject\n",
    "- ``r``: optimization rank parameter. We use the suggested value from the eigenvalues analysis, but it can be any integer smaller than ``number_of_sources_to_localize``."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "70553011409ac322"
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001B[36mCalculating activity index for localizer: mpz_mvp\u001B[0m\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5124/5124 [00:16<00:00, 302.37it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 311.68it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 311.65it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 306.47it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 309.97it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 307.57it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 303.06it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 304.61it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 303.29it/s]\n",
      "100%|██████████| 5124/5124 [00:16<00:00, 302.00it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 301.35it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 295.38it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 294.21it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 295.95it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 293.16it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 291.34it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 282.57it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 287.68it/s]\n",
      "100%|██████████| 5124/5124 [00:17<00:00, 285.50it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 269.37it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 268.82it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 272.12it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 275.31it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 278.55it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 274.03it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 277.67it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 272.22it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 272.04it/s]\n",
      "100%|██████████| 5124/5124 [00:18<00:00, 271.20it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 266.81it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 266.88it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 264.43it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 264.33it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 264.21it/s]\n",
      "100%|██████████| 5124/5124 [00:19<00:00, 260.33it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 231.19it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 226.51it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 230.34it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 231.64it/s]\n",
      "100%|██████████| 5124/5124 [00:23<00:00, 215.33it/s]\n",
      "100%|██████████| 5124/5124 [00:30<00:00, 168.64it/s]\n",
      "100%|██████████| 5124/5124 [00:28<00:00, 181.66it/s]\n",
      "100%|██████████| 5124/5124 [00:33<00:00, 152.17it/s]\n",
      "100%|██████████| 5124/5124 [00:34<00:00, 146.65it/s]\n",
      "100%|██████████| 5124/5124 [00:36<00:00, 141.51it/s]\n",
      "100%|██████████| 5124/5124 [00:38<00:00, 131.63it/s]\n",
      "100%|██████████| 5124/5124 [00:39<00:00, 129.57it/s]\n",
      "100%|██████████| 5124/5124 [00:36<00:00, 139.15it/s]\n",
      "100%|██████████| 5124/5124 [00:38<00:00, 133.81it/s]\n",
      "100%|██████████| 5124/5124 [00:38<00:00, 134.63it/s]\n",
      "100%|██████████| 5124/5124 [00:38<00:00, 133.06it/s]\n",
      "100%|██████████| 5124/5124 [00:39<00:00, 128.37it/s]\n",
      "100%|██████████| 5124/5124 [00:37<00:00, 136.43it/s]\n",
      "100%|██████████| 5124/5124 [00:37<00:00, 137.27it/s]\n",
      "100%|██████████| 5124/5124 [00:37<00:00, 136.45it/s]\n",
      "100%|██████████| 5124/5124 [00:38<00:00, 132.86it/s]\n",
      "100%|██████████| 5124/5124 [00:39<00:00, 129.94it/s]\n",
      "100%|██████████| 5124/5124 [00:39<00:00, 128.66it/s]\n",
      "100%|██████████| 5124/5124 [00:40<00:00, 126.19it/s]\n",
      "100%|██████████| 5124/5124 [00:41<00:00, 123.15it/s]\n",
      "100%|██████████| 5124/5124 [00:42<00:00, 121.05it/s]\n",
      "100%|██████████| 5124/5124 [00:45<00:00, 113.27it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Leadfield indices corresponding to localized sources: [31, 3557, 795, 1213, 1690, 2506, 2697, 2602, 1225, 83, 1966, 2085, 994, 2850, 2212, 4404, 4304, 4882, 608, 2522, 2325, 1876, 6, 3255, 1860, 3714, 4804, 2690, 84, 371, 2454, 3624, 5108, 1971, 265, 650, 992, 2727, 42, 4698, 329, 2574, 1258, 4002, 3368, 4920, 3, 33, 1405, 2887, 2333, 2771, 4604, 3333, 302, 306, 4422, 2597, 2580, 16, 4266, 3815]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "locs_sen = localizer.localize(\n",
    "    subject=subject,\n",
    "    subjects_dir=subjects_dir,\n",
    "    localizer_to_use=[\"mpz_mvp\"],\n",
    "    n_sources_to_localize=sugg_n_sources,\n",
    "    R=data_cov_sen.data,\n",
    "    N=noise_cov.data,\n",
    "    forward=fwd,\n",
    "    r=sugg_rank\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.706360Z",
     "start_time": "2025-09-11T21:10:03.407450Z"
    }
   },
   "id": "a7ac386cf184320c"
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [],
   "source": [
    "# Transform leadfield indices to vertices\n",
    "lh_vert_sen, lh_idx_sen, rh_vert_sen, rh_idx_sen = utils.transform_leadfield_indices_to_vertices(\n",
    "    lf_idx=locs_sen[\"sources\"],\n",
    "    src=src,\n",
    "    hemi=\"both\",\n",
    "    include_mapping=True\n",
    ")\n",
    "\n",
    "locs_sen.add_vertices_info(\n",
    "    lh_vertices=lh_vert_sen,\n",
    "    lh_indices=lh_idx_sen,\n",
    "    rh_vertices=rh_vert_sen,\n",
    "    rh_indices=rh_idx_sen\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.713803Z",
     "start_time": "2025-09-11T21:36:19.712555Z"
    }
   },
   "id": "40a1929b9cc56c6b"
  },
  {
   "cell_type": "markdown",
   "source": [
    "Optionally, we can plot the localized sources on the brain surface:"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "ec69bb6603ecf240"
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [],
   "source": [
    "# locs_sen.plot_localized_sources()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.718550Z",
     "start_time": "2025-09-11T21:36:19.713504Z"
    }
   },
   "id": "5cc5f1d4eda54f64"
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, the size and color of the markers indicate the order of localization:\n",
    "- large, red foci: sources localized earlier\n",
    "- small, white foci: sources localized later"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "7109589ea45acf66"
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Reconstruct\n",
    "Now that we have localized sources of interest, the next sgtep is to **reconstruct their activity**.\n",
    "First, we restrict the original forward model to only include the localized sources. This reduces the forward solution to the relevant subspace:|"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "f782df22e0978f93"
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "# Subset mne.Forward\n",
    "new_fwd_sen = utils.subset_forward(\n",
    "    old_fwd=fwd,\n",
    "    localized=locs_sen,\n",
    "    hemi=\"both\"\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.795935Z",
     "start_time": "2025-09-11T21:36:19.716622Z"
    }
   },
   "id": "eb99cd2d30ffe436"
  },
  {
   "cell_type": "markdown",
   "source": [
    "To compute the filters, we will use ``beamformer.make_filter``. \n",
    "\n",
    "This function works similarly to ``mne.beamformer.make_lcmv``, but with additional parameters specific to MVPURE.\n",
    "\n",
    "We provide these in a dictionary called ``mvpure_params``:\n",
    "- ``filter_type``: type of the beamformer to use. Options are: ``MVP_R`` and ``MVP_N``. In this case, we will use ``MVP_R`` as it is generalization of commonly used LCMV filter.\n",
    "- ``filter_rank``: optimization rank parameter. For best performance, we use the same rank as in the localization step.\n",
    "\n",
    "Note: setting ``filter_rank=\"full\"`` reduces the method to a standard LCMV filter.\n",
    "For theoretical details see [PAPER]."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "7b72fd38ee20c861"
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [],
   "source": [
    "# MVPURE filter parameters\n",
    "mvpure_params = {\n",
    "    'filter_type': 'MVP_R',\n",
    "    'filter_rank': sugg_rank\n",
    "}"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.808188Z",
     "start_time": "2025-09-11T21:36:19.737964Z"
    }
   },
   "id": "f92cd1accedf6dc0"
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing rank from covariance with rank=None\n",
      "    Using tolerance 4.4e-13 (2.2e-16 eps * 128 dim * 15  max singular value)\n",
      "    Estimated rank (eeg): 86\n",
      "    EEG: rank 86 computed from 128 data channels with 1 projector\n",
      "Computing rank from covariance with rank=None\n",
      "    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)\n",
      "    Estimated rank (eeg): 86\n",
      "    EEG: rank 86 computed from 128 data channels with 1 projector\n",
      "Making MVP_R beamformer with rank {'eeg': 86} (note: MNE-Python rank)\n",
      "Computing inverse operator with 128 channels.\n",
      "    128 out of 128 channels remain after picking\n",
      "Selected 128 channels\n",
      "Whitening the forward solution.\n",
      "    Created an SSP operator (subspace dimension = 1)\n",
      "Computing rank from covariance with rank={'eeg': 86}\n",
      "    Setting small EEG eigenvalues to zero (without PCA)\n",
      "Creating the source covariance matrix\n",
      "Adjusting source covariance matrix.\n",
      "Computing beamformer filters for 62 sources\n",
      "MVP_R computation - in progress...\n",
      "Filter rank: 42\n",
      "Filter computation complete\n"
     ]
    }
   ],
   "source": [
    "# Build beamformer filter (similar to LCMV but with MVPURE options)\n",
    "filter_sen = beamformer.make_filter(\n",
    "    signal_sen.info,\n",
    "    new_fwd_sen,\n",
    "    data_cov_sen,\n",
    "    reg=0.05,\n",
    "    noise_cov=noise_cov,\n",
    "    pick_ori=None,  # not needed with fixed orientation forward\n",
    "    weight_norm=None,\n",
    "    rank=None,\n",
    "    mvpure_params=mvpure_params\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.868241Z",
     "start_time": "2025-09-11T21:36:19.745153Z"
    }
   },
   "id": "faa8b0c9ebbd7671"
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [],
   "source": [
    "# Apply filter to cropped evoked response\n",
    "stc_sen = beamformer.apply_filter(signal_sen, filter_sen)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.887427Z",
     "start_time": "2025-09-11T21:36:19.866069Z"
    }
   },
   "id": "cdb609b79a691207"
  },
  {
   "cell_type": "markdown",
   "source": [
    "We then attach the resulting ``mne.SourceEstimate`` to the localized sources object, making it easier to visualize:"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "d5d1a62c564e44a"
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [],
   "source": [
    "# Add reconstructed source time course\n",
    "locs_sen.add_stc(stc_sen)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:19.888613Z",
     "start_time": "2025-09-11T21:36:19.886324Z"
    }
   },
   "id": "7553c0bd5dcc339c"
  },
  {
   "cell_type": "markdown",
   "source": [
    "Finally, let's plot the localized sources with their reconstructed activity:"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "664b440fc132f9f5"
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using control points [2.91332758e-09 3.43983961e-09 7.89233742e-09]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<mne.viz._brain._brain.Brain at 0x1583aa270>"
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "viz.plot_sources_with_activity(\n",
    "    subject=subject,\n",
    "    stc=stc_sen,\n",
    "    background=\"white\"\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:23.240780Z",
     "start_time": "2025-09-11T21:36:19.890610Z"
    }
   },
   "id": "c0dd032cd6283f1e"
  },
  {
   "cell_type": "markdown",
   "source": [
    "### \"Cognitive\" task\n",
    "After examing the early sensory responses, we now turn to the later cognitive stahe of processing in the oddball paradigm.\n",
    "In EEG, target stimuli typically evoke a P300 component — a positive deflection peaking around 300–600 ms after stimulus onset.\n",
    "This response is thought to reflect higher-level cognitive processes, such as attention allocation and stimulus evaluation, in contrast to the earlier sensory responses.\n",
    "\n",
    "For this dataset, we will therefore define the cognitive time window as **350–600 ms**.\n",
    "The pipeline remains the same as before:\n",
    "\n",
    "- Compute noise covariance (always from −200 to 0 ms).\n",
    "- Compute data covariance in the cognitive window (350–600 ms).\n",
    "- Subset the evoked signal to this time range."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "ad22f685ecd15c6f"
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Created an SSP operator (subspace dimension = 1)\n",
      "    Setting small EEG eigenvalues to zero (without PCA)\n",
      "Reducing data rank from 128 -> 127\n",
      "Estimating covariance using EMPIRICAL\n",
      "Done.\n",
      "Number of samples used : 5070\n",
      "[done]\n"
     ]
    }
   ],
   "source": [
    "# Compute data covariance for range corresponding to sensory processing\n",
    "data_cov_task = mne.compute_covariance(\n",
    "    sel_epoched,\n",
    "    tmin=0.35,\n",
    "    tmax=0.6,\n",
    "    method=\"empirical\"\n",
    ")\n",
    "# There's no need to compute `noise_covariance` again as it is the same time interval\n",
    "\n",
    "sel_evoked = sel_epoched.average()\n",
    "\n",
    "# Subset signal for given time range\n",
    "signal_task = sel_evoked.crop(\n",
    "    tmin=0.35,\n",
    "    tmax=0.6\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:23.360636Z",
     "start_time": "2025-09-11T21:36:23.246567Z"
    }
   },
   "id": "237a6b2afb47240e"
  },
  {
   "cell_type": "markdown",
   "source": [
    "From here, we can repeat the same steps as in the sensory section:\n",
    "- analyze eigenvalues of $RN^{-1}$,\n",
    "- localize sources,\n",
    "- reconstruct signals with MVPURE filters,\n",
    "- and finally visualize the results."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "141d93859715c998"
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1200x600 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8gAAAITCAYAAADb3DVaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAASshJREFUeJzt3Qt43FWdN/CT9A40oVBqSylawAVKsYqii7CAwoosFsSCN1y5uAJSbvIKCiosrFx9XVkLIrCK64qioCh0X1RwpdiV+x3LtbAWIaVUadNS2tJk3uf3jxMn6UyapEnm9vk8zzCd/9zOzJyE+eac8zsNuVwulwAAAKDONZa7AQAAAFAJBGQAAAAQkAEAAKCDgAwAAAACMgAAAHQQkAEAAEBABgAAgA4CMgAAAAjIAAAA0EFABgAAAAEZAAAAOgjIAEDVu+KKK9Juu+2WRowYkf75n/+53M0BoEoJyABA1Zs0aVIWjGfNmlXupgBQxQRkgD6IL+ANDQ1p6dKlPd7uu9/9bna7//3f/x2yttWLanlve9tXquk1VbIPfvCD6eCDD06bb775gD92rX8+tf76APpCQAZq2qOPPpoOO+yw9MY3vjGNHj06TZ48Of393/99mjNnTqpmv/vd77IAtmzZsnI3hTrvL/lwlT8NHz48+zk76qij0gsvvFDy9vHzWOz6fffdN02fPn3A2wkAvSEgAzUdCt7xjnekhx9+OH36059Ol112Wfqnf/qn1NjYmP7t3/5tUJ/7H//xH9Nrr72WBfPBem3nnnuugFwjaqG/nHfeeek///M/07e+9a104IEHpu9///tpn332SatXry56+zVr1qSLLrooVYPB/nwAqBzDy90AgMFy/vnnp+bm5nTvvfeuN+1yyZIlg/rcw4YNy05QL/0lQnH8QSrEH6LGjx+fLr744nTTTTelD3/4w+vd/q1vfWu6+uqr05lnnpm23nrrko+71157pf/5n/8pet0Xv/jF9JWvfCUNtlr4fADoHSPIQM1auHBh2mWXXYquSZwwYULnv2Mq6Jve9KaSa0iLiXWl8aW/qakpbbnllumUU07pMlJWak1fTCk95phj0hve8IY0atSorH3f+c531nv8uN2nPvWpLDjE7aZOnZo+85nPpLVr12btOv3007PbxfH81NbC53riiSfSokWLNvgerVixIp166qnZ64/nifclpqA/8MAD2fV/+MMf0gknnJB23HHHNGbMmOy1Hn744eu9rvx79dRTT6VPfOIT2R8mttpqq/TlL3855XK59Pzzz6dDDjkke78mTpyYvva1r5V8v6PtPb23pfT2vd2Y96M//WVDfaUS+stg+Lu/+7vOn8NizjrrrNTW1rbBUeT58+dnfajYaSDCcW/e41Kfz+233579USCmi2+//fbpyiuvLNkPevM8+fs+88wzWT+L313xs3T00UenVatWZbe54YYbstvMmzdvveeI54/rHnvssT79/BbT137em9fXm58vgHIzggzUrJgOeeedd2ZfFgd6TWMEnviSd+GFF6a77rorfeMb30ivvPJK+t73vlfyPi+99FL627/92+zL5YknnpgFyFtuuSULNq2trdkXx/Diiy+md77zndl02GOPPTbttNNO2ZfP+GIcX5I/9KEPZUH0hz/8Yfr617+ejdSFeLy8nXfeOZveGl/ge3L88cdnjxvtmTZtWvrTn/6UBZLHH3882zInRt9jeu5HP/rRtM0222RfrGM7nVgnumDBgrTJJpt0ebyPfOQj2XNH6Pmv//qvLMBsscUW2Rf39773vdmI4rXXXps+97nPpd133z3tvffeQ/rebsiG3o/+6M/rGer+MhjyIWzcuHFFr4+w/slPfjIbRf7CF77Q4yhyb6xbty47ReiO8/gjRGz51NPI78b0mwcffDC9//3vz6pnx/T1eN6YZl7sfe3r80Sfifcn+kyEx3//93/PwmT8/Bx00EFps802Sz/+8Y+zn/FCP/rRj7Jgmv9919ef3/7q7esbjJ8vgAGXA6hRv/rVr3LDhg3LTnvssUfujDPOyP3yl7/MrV27tsvtjjzyyNwb3/jG9e5/zjnn5Lr/mswfO/jgg7scP+GEE7LjDz/8cHb5mmuuyS4/99xznbf51Kc+lZs0aVJu6dKlXe770Y9+NNfc3JxbtWpVdvmTn/xkrrGxMXfvvfeu16b29vbs/Ktf/ep6j18orttnn302+B7F886ePbvk9fk2Fbrzzjuzx//e97633vty7LHHdh5bt25dbptttsk1NDTkLrroos7jr7zySm7MmDHZ+z7U7+3Gvh996S+9fT2V0F82Rr7tt912W+7ll1/OPf/887kbbrght9VWW+VGjRqVXS52+2jvwoULc8OHD8+dfPLJnddHv91ll1363I78+114iufqSW/f42Kfz8yZM3ObbLJJ7oUXXug89vTTT2evp/vvjd4+T/41HHPMMV1ud+ihh+a23HLLzssf+9jHchMmTMh+xvJaWlqyfnDeeef1+ee32Ovry+/F3r6+3vx8AZSbKdZAzYqpezGCHFu/RKGuSy65JB1wwAFZhd1YF7kxZs+e3eXySSedlJ3/v//3/4rePjLrT37ykzRz5szs3zHtNn+KNi1fvjwbKWpvb08/+9nPstvl13MWKjXlu9jzbWj0OMQUzrvvvjsbhSwmpmXmvf7669mIzw477JDdr9i0yFh7mhcjd/Eaoi0xklT4nDHl89lnnx3S97Y3NvR+9EdfX085+stA2H///bORwylTpmSV4zfddNPs5yxGLkvZbrvtsgJYV111VWppadmo54+pv92nYMc04VI2pt/EaPFtt92WbS1VOPIdPxuxFntjnydGWrtPV4+fvRiNzc/UiDoKhT/jMTIb/SGu6+/Pb3/05fUNxs8XwEATkIGaFtN4f/rTn2ZTWu+5556sIFCsg4sv8DHFsL/e/OY3d7kc6w+jOnaptX0vv/xyNgU2gkCEiMJTrC8M8YU3bhdfgodqm5v4o0FMQY9QE9N0I2QUBteo3Hv22Wdn18eawZieG22O1xJffLvbdtttu1yO9ZOxPjM/rbfweHwmQ/neDsT70R99fT3l7i+xbnnx4sVdThEIN+Tyyy9Pt956axbU/uEf/iELSNFnNuRLX/pSNiV6qCtab0y/iePxsxFhs7vux/rzPN1/jvLT1PM/MzG1O36GYkp1Xvw7Cp/9zd/8Tb9/fvujL69vMH6+AAaaNchAXRg5cmQWluMUXyDji9v111+fzjnnnJKjbL0JBb0dqYuRnRAFrI488siit3nLW96SjcAMpVjrGKNTN954Y/rVr36VvvrVr2brHOOPCjESFqOd11xzTbaGcI899si+lMdrjTWN+ddUqNh6z1JrQHv7WgfqvR2I96On9vS2v/RmVLec/SXWrL7nPe/pcuy5554rWrCpUASe/Ch2jKxG9emPf/zj6cknn8zWzPY0ihyvMwJWrEUeKgPZbwb6eTb0MxNhN97j6Kff/OY3szXAUen7ggsu6HL7vv78FuptP+/L6+vNzxdAuQnIQN3Jf4nPT+mM0Zli+8NGBdhSnn766ayITl5UnY0viqVCRIymjB07NvtyGVNRS4nHiGrH+Sq0pQzk1NkoMhSVbuMUIz1RLCe2yIovrDEaGF96C6tOR/GjwdxPd7De24F4P/rTX/r6esrdX2bMmJGNBBeKyuN9EQEvCkxF0I79xzcUfGMUOfZNjrA0VDam30TBrJgZEZ9ld92PDXT/zIup1P/xH/+Rfv3rX2dFriI8F06vDhvz89vbft7X17ehny+AcjPFGqhZv/nNb4qOsOXXfsY62PyU15hu+Mgjj3TeJsJzjHL0NJ200Jw5c7LzUl/yIjDMmjUrW6tXLMzENMUQU29jZOjmm29O991333q3y7+eWN8ZSn3R7c02T/GFtvs0y/jiH2sq16xZ09nu7u9hvNa+jK731WC9txvSm/ejP/2lr6+nHP2lezCKoFN4ijDYV1EpOUaVL7300g1u0xXvaYxARrXzmNI9FDam38R9432J9d+F62kjHEf15oF6np7E80eF+JhaHad4rwv/ELOxP7+97ee9fX29/fkCKDcjyEDNiumFsc3NoYcemm19E2srY/pofJmM0bv8+riYbvj5z38+u93JJ5+c3Se2Qomp2KUK2cSU0yj+FWsBoxBYjH7FdNIYfSsl1lhGaH/Xu96VPv3pT2fbnPz5z3/OniMK/sS/Q0yTjOmHsYVLbNsT2ybFF9OYEh5bokShm7e//e3Zbb/4xS9m7Y/tbKJITj4I9Wabp1iLHQWUYj12tDumwUY7YmuY/IjTBz7wgfSf//mf2dTMaG+81rhN7Kc6WAbzve1Jb96P/vSX/ryeoe4vgyX2X459d2Mf4e6Fp7qLtkVfiynZsVXRUNiYfhPrZ+N933PPPbM9pyMAxmh5rAd/6KGHBux5SonPMLbwuu6669Krr76a/u///b/r3WZjfn770s978/p6+/MFUHblLqMNMFhuueWWbLuUnXbaKbfZZpvlRo4cmdthhx1yJ510Uu6ll15ab0uo6dOnZ7fZcccdc9///vd73OZpwYIFucMOOyw3duzY3Lhx43Innnhi7rXXXutx25QQzxvbnEyZMiU3YsSI3MSJE3P77bdf7qqrrupyuz/84Q/Z9j35rXK222677H5r1qzpvM2//Mu/5CZPnpxt7dL9uXqzzVM81umnn56bMWNG9jo23XTT7N/f/OY3u2zJdPTRR+fGjx+fvYcHHHBA7oknnsi2fyncpin/vsQ2P4XiNvG43RXbymeo3tuNeT/60l96+3oqob9sjMJtm7pra2vLbb/99tkpvyVRT7eP/hLX9Webp/7qzXtc6vP59a9/nXvb296W9YN4jf/+7/+e+z//5//kRo8e3a/nKfVzVOr5b7311ux4bKXWfTutvvz8lnr83v5e7M3r68vPF0A5NcR/yh3SAWrNt7/97WzLo+eff77HbW7oOiJ37rnnZlMyu1e9rnX6S+2IKe+///3vs7XnAFQfa5ABBkFMcY3CSLFGEDZEf6lOsY1SoQjFUeMg1l8DUJ2sQQYYQLHdSlSO/da3vpVtq7LJJpuUu0n8pUDQhoohxZrInrYjGgz6S3WLLaqOOuqo7DyqO8ca3dhS7owzzih30wDoJwEZYADFditRmCgqyl599dXlbg5/EVOXu1f47S72xI5p3kNJf6luUXjthz/8YVZ5O/Ymjj9yRNG0N7/5zeVuGgD9ZA0yADUvthmKis49iVHAOAEA9UtABgAAgHJMsW5vb08vvvhiGjt2bFaQBAAAAAZTjAvHnuxbb711amxsrJyAHOF4ypQpQ/20AAAA1LnnN7Cl4pAH5Bg5zjesqakplfuvCMuXL0/Nzc1Gs9Ef6EJ/oDt9gkL6A4X0BwrpD5WptbU1G6jN59GKCcj5ThLhuBICcpyiHTov+gOF9Ae60ycopD9QSH+gkP5Q2Tb0mZSefA0AAAB1REAGAAAAARkAAADKtAYZAACgErS1taXXX399QB8z1h+vXbs2rV692hrkITRixIg0bNiwjX4cARkAAKgrEWIXL16cli1bNiiP397env70pz8NymNT2uabb54mTpy4UX+YEJABAIC6kg/HEyZMSJtsssmAjvRG+I6R6RjNNII8NOI9X7VqVVqyZEl2edKkSf1+LAEZAACoGxFe8+F4yy23HPDHF5DLY8yYMdl5hOT4bPs73VqRLgAAoG7k1xzHyDG1ZZO/fKYbs65cQAYAAOqO0d3a0zAAn6mADAAAAAIyAAAAdBCQAQAAQEAGAACgmEMPPTSNGzcuHXbYYaleCMgAAACs55RTTknf+973Uj0RkAEAAFjPvvvum8aOHZvqiYBcxMKXV6Yzbng4fWDOb7PzuAwAAFBORx11VLaVUZxGjBiRpk6dms4444y0evXq9W5z0UUXdbnvz372s7JtbXXHHXekmTNnpq233jprQ7SlUg0vdwMqTYThmXPmpzXr2lNbey493rIizX2kJd180l5p+602K3fzAACAOvb+978/XXPNNen1119P999/fzryyCOz0HnxxRd33mb06NHZ5eOOOy5bQ1zKW9/61rRu3br1jv/qV7/KwuxAefXVV9OMGTPSMccckz70oQ+lSiYgd3PlvIWd4TjEeVyO45ccNqPczQMAAOrYqFGj0sSJE7N/T5kyJe2///7p1ltv7RKQ49gzzzyTLrzwwnTJJZeUfKyHHnpoSNp84IEHZqdqYIp1NwtaWjvDcV5cXtCyomxtAgAA6O6xxx5Lv/vd79LIkSO7HB82bFi64IIL0pw5c9If//jHAX/eCy64IG222WY9nhYtWpSqkRHkbqZNasqmVReG5GGNDWnapPpanA4AAPVk3apVA/I4uVwutbW1pdywYT2u+R2+ySb9evy5c+dmATSmRq9ZsyY1Njamyy67rOgWTTGF+pxzzknf/va3+/Vc+++/f3r44YezKdLbbLNNuv7669Mee+yRjj/++PThD3+4x/sO5BTtoSQgd3PcPttna47z06wjHI8a3pgdBwAAatMtu+46pM83c+HCft3vPe95T7riiiuy0Pr1r389DR8+PM2aNavobWPa9Xvf+970uc99rl/PddtttxU9vsUWW2SnWmSKdTdRiCsKcs3abXKaPrk5O1egCwAAqASbbrpp2mGHHbKiV9/5znfS3XffXXKEeO+9904HHHBAOvPMMwe0DReYYl1fIgwryAUAAPXjwEcfHdAp1rEOeLC3VYrp1WeddVY67bTT0sc//vE0ZsyY9W4T2z3FVOsdd9xxwJ73eFOsAQAAald/1wQXC8gNQxSQw+GHH55OP/30dPnllxedSr3rrrumI444In3jG98YsOfcoo9TrFeuXJlV1c577rnnsgra8RjbbrttqiSmWAMAAFSpWIN84oknZts5xbrkYs4777zU3t6eyuW+++5Lb3vb27JTiBHv+PfZZ5+dKk1DLv7EMYRaW1tTc3NzWr58eWpqakrlFC892hHtGYq/7lDZ9AcK6Q90p09QSH+gkP5QXVavXp2NYE6dOjWNHj16wB9/KKdY0/vPtrc51AgyAAAACMgAAADQQUAGAAAAARkAAAA6CMgAAAAgIAMAAEAHARkAAAAEZAAAAOggIAMAAICADAAAAB0EZAAAABCQAQAAKObQQw9N48aNS4cddliqFwIyAAAA6znllFPS9773vVRPBGQAAADWs++++6axY8emeiIgAwAAVIGjjjoqNTQ0ZKcRI0akqVOnpjPOOCOtXr16vdtcdNFFXe77s5/9LDteDldccUV6y1vekpqamrLTHnvskW655Zaitz366KPTl770pXThhRem3XffPQvoEyZMSB/84AfTk08+OehtFZABAACqxPvf//7U0tKSnn322fT1r389XXnllemcc87pcpvRo0eniy++OL3yyis9PtZb3/rWNH369PVOL7744oC2eZtttskC+/3335/uu+++9N73vjcdcsgh6fe//32X27W1taW5c+emgw8+OM2bNy/Nnj073XXXXenWW29Nr7/+enrf+96XXn311TSYhg/qowMAADBgRo0alSZOnJj9e8qUKWn//ffPAmQE4rw49swzz2SjsJdccknJx3rooYeGpM0zZ87scvn888/PRpUj/O6yyy6dx3/3u99lI+MxcvyLX/yiy32++93vZiPJEbL33nvvQWurEWQAAIAq9Nhjj2WhcuTIkV2ODxs2LF1wwQVpzpw56Y9//OOAP+8FF1yQNttssx5PixYtKnrfGCW+7rrrspHgmGpd6KabbsrCdLGp4MuXL8/Ot9hiizSYjCADAAB17/W21wfkcXK5XBYC21N7j2t+Rwwb0a/HjynIEUDXrVuX1qxZkxobG9Nll11WdIummEId06+//e1v9+u59t9///Twww9nYTamSV9//fVZqD3++OPThz/84R7vu/XWW3e5/Oijj2b3jfXS0f4bb7wxTZs2rcttfv7zn2fTxrtrb29Pp556atpzzz2zKeCDSUAGAADq3uW3Xz6kz3fqfqf2637vec97sunJEVojTA4fPjzNmjWr6G1j2nWs9/3c5z7Xr+e67bbbih6PUdy+juTuuOOO2ZTuGAm+4YYb0pFHHpmtM86H5Mcffzxb+7zffvutd99Yixyj5fPnz0+DzRRrAACAKrHpppumHXbYIc2YMSN95zvfSXfffXfJEeJYq3vAAQekM888s+xTrEeOHJm1++1vf3u2Njra/2//9m9dplf//d//fVZgrNCJJ56YjZr/5je/yUaxB5sRZAAAoO7N3nf2gE6xjnXAg72tUkyvPuuss9Jpp52WPv7xj6cxY8asd5uoHh1TrWMEd6Ac348p1sWmTccU8cLp1ccee2yX9/Gkk07KpmLffvvt2ZZWQ0FABgAA6l5/1wR3F8GuMTUOSUAOhx9+eDr99NPT5ZdfXnQq9a677pqOOOKI9I1vfGPAnnOLPk6xjhHsAw88MG277bZpxYoV6Qc/+EEWen/5y19m1y9ZsiTb/ilGkQunVcftIjjHXsiLFy/Ojjc3Nxf9Q8BAMcUaAACgSsUa5JiGHNs5ldoj+LzzzstGbMtlyZIl6ZOf/GQ2ih1rjO+9994sHMeU6nDzzTend77znWn8+PGd94l11rFeed99902TJk3qPP3oRz+qnBHkmCrwz//8z+n73/9+luBj2Pyoo45KX/rSl4bkryMAAAD1KvYCLuYLX/hCdip1mze96U1dpjMPtW9voIp2jBIffPDB643El0OfAnJUQYsk/x//8R/Zhs4xDH700Udnw9wnn3zy4LUSAACAmrTXXnulj33sY6kS9CkgxybUhxxySDrooIM6/xLxwx/+MN1zzz2D1T4AAABq2BlnnJEqRZ/WIL/73e9Ov/71r9NTTz2VXY5No2MvqlhwDQAAANWsTyPIMa+9tbU17bTTTllVtliTfP7552dV0UqJue6F893j/vk55eWaV56Xb0O520Fl0B8opD/QnT5BIf2BQvpDdcl/TkPxmekTlfPZ9vaz6FNA/vGPf5yuvfbarNx2rEF+6KGH0qmnnpoV6zryyCOL3ic2gT733HPXOx4VycrdYeL5V65cmf1bkTH0BwrpD3SnT1BIf6CQ/lBd1q5dm1V0jsG+OA2GclaMrmdtbW3Zex9bSXUvSpYfqN2QhlwfUuqUKVOyUeTYkyrvK1/5SlbV+oknnuj1CHI8zrJly1JTU1Mqp3jpEdSjyJhfZugPFNIf6E6foJD+QCH9obqsXr06/e///m9WT2n06NGDFtRixi2V89lGDt18882zn9WecmifRpBXrVqVGhu7LluOD76nv5CMGjUqO3UXvzwq4RdIvh2V0BbKT3+gkP5Ad/oEhfQHCukP1SP/GQ3W51U4/qg/VM5n29vPok8BeebMmdma42233TabYv3ggw+mf/3Xf03HHHNMXx4GAACgrEyDrj3tA/CZ9ikgz5kzJ335y19OJ5xwQlqyZEm29vi4445LZ5999kY3BAAAYLCNHDkymxX74osvpq222iq7PJAjvTGCnJ9ibQR5aMR7HmvLX3755eyzjc90SALy2LFj06WXXpqdAAAAqk0EqKlTp6aWlpYsJA/WSGb3pakMvk022SSb7bwx732fAjIAAEC1ixHGCFLr1q0b8ErWMZoZVZRjcNEI8tCJEfvhw4dv9HsuIAMAAHUngtSIESOy00AH5NjFJ6ooC8jVx7g/AAAACMgAAADQQUAGAAAAARkAAAA6CMgAAAAgIAMAAEAHARkAAAAEZAAAAOggIAMAAICADAAAAB0EZAAAABCQAQAAoIOADAAAAAIyAAAAdBCQAQAAQEAGAACADgIyAAAACMgAAADQQUAGAAAAARkAAAA6CMgAAAAgIAMAAEAHARkAAAAEZAAAAOggIAMAAICADAAAAB0EZAAAABCQAQAAoIOADAAAAAIyAAAAdBCQAQAAIKU0vNwNqCYLX16Zrpy3MC1oaU3TJjWl4/bZPm2/1WblbhYAAAADQEDuQzieOWd+WrOuPbW159LjLSvS3Eda0s0n7SUkAwAA1ABTrHspRo7z4TjEeVyO4wAAAFQ/AbmXYlp1PhznxeUFLSvK1iYAAAAGjinWvRRrjmNadWFIHtbYkKZNGpv92/pkAACA6iYg91IE3lhznJ9mHeF41PDG7Lj1yQAAANXPFOteiqAbgXfWbpPT9MnN2Xk+AFufDAAAUP2MIPdBhOFLDpux3nHrkwEAAKqfEeQBEGuOY8p1ocL1yQAAAFQ+AXkAxDrkWI+cD8mF65MBAACoDqZYD+D65I4q1iuykWNVrAEAAKqLgDzI65MBAACoDqZYAwAAgIAMAAAAHQRkAAAAEJABAACgg4AMAAAAAjIAAAB0EJABAABAQAYAAIAOAjIAAAAIyAAAANBBQAYAAAABGQAAADoIyAAAACAgAwAAQIfhfzlnEC18eWW6ct7CtKClNU2b1JSO22f7tP1Wm5W7WQAAABQQkIcgHM+cMz+tWdee2tpz6fGWFWnuIy3p5pP2EpIBAAAqiCnWgyxGjvPhOMR5XI7jAAAAVA4BeZDFtOp8OM6LywtaVpStTQAAAKxPQB5kseZ4WGNDl2NxedqksWVrEwAAAOsTkAdZFOQaNbyxMyTHeVyO4wAAAFQORboGWRTiioJcHVWsV2Qjx6pYAwAAVB4BeQhEGL7ksBnlbgYAAAA9MMUaAAAABGQAAADoICADAACAgAwAAAAdBGQAAABQxbq8Fr688i/bP7WmaZOabP8EAABQRgJyGcPxzDnz05p17amtPZceb1mR5j7Sku2ZLCQDAAAMPVOsyyRGjvPhOMR5XI7jAAAADD0BuUxiWnU+HOfF5QUtK8rWJgAAgHomIJdJrDke1tjQ5VhcnjZpbNnaBAAAUM8E5DKJglyjhjd2huQ4j8txHAAAgKGnSFeZRCGuKMjVUcV6RTZyrIo1AABA+QjIZRRh+JLDZpS7GQAAAJhiDQAAAB0EZAAAABCQAQAAoIOADAAAAAIyAAAA9DMgv/DCC+kTn/hE2nLLLdOYMWPSrrvumu67776+PgwAAABU7zZPr7zyStpzzz3Te97znnTLLbekrbbaKj399NNp3Lhxg9dCAAAAqLSAfPHFF6cpU6aka665pvPY1KlTB6NdAAAAULkB+aabbkoHHHBAOvzww9O8efPS5MmT0wknnJA+/elPl7zPmjVrslNea2trdp7L5bJTOeXbUO52FLPw5ZXpqjsWpsdbWtPOk5rSsXtvn7bfarNyN6umVXJ/YOjpD3SnT1BIf6CQ/kAh/aEy9fbz6FNAfvbZZ9MVV1yRTjvttHTWWWele++9N5188slp5MiR6cgjjyx6nwsvvDCde+656x1fvnx52TtNPP/KlSuzfzc0NKRK8cdXVqVTf/RQen1dLrXlcunu5SvSA8+8mC79yFvTNuM2KXfzalal9gfKQ3+gO32CQvoDhfQHCukPlSk/ULshDbk+pNQIwu94xzvS7373u85jEZAjKN955529HkGOadrLli1LTU1NqZzipUdQb25urqjO+/mfPJx++sALqa39rx/NsMaG9KHdJqeLZ80oa9tqWaX2B8pDf6A7fYJC+gOF9AcK6Q+VKXLo5ptvnn02PeXQPo0gT5o0KU2bNq3LsZ133jn95Cc/KXmfUaNGZafuorNUQofJt6MS2pK3oGVFWtce//prm+LygpaVFdXOWlSJ/YHy0R/oTp+gkP5AIf2BQvpD5entZ9GnbZ6igvWTTz7Z5dhTTz2V3vjGN/atdfRo2qSmbMS4UFyeNmls2doEAABQ6/oUkD/72c+mu+66K11wwQXpmWeeST/4wQ/SVVddlWbPnj14LaxDx+2zfRo1vLEzJMd5XI7jAAAADI4+TbHefffd04033pjOPPPMdN5552VbPF166aXpiCOOGKTm1aeoVn3zSXulK+ctzKZbx8hxhGNVrAEAACokIIcPfOAD2YnBFWH4ksMU5AIAAKjIKdYAAABQqwRkAAAAEJABAACgg4AMAAAAAjIAAAB0EJABAABAQAYAAIAOAjIAAAAIyAAAANBBQAYAAAABGQAAADoIyAAAACAgAwAAQAcBGQAAAARkAAAA6CAgAwAAgIAMAAAAHQRkAAAAEJABAACgg4AMAAAAAjIAAAB0GP6Xc6rIwpdXpivnLUwLWlrTtElN6bh9tk/bb7VZuZsFAABQ1QTkKgzHM+fMT2vWtae29lx6vGVFmvtIS7r5pL2EZAAAgI1ginWViZHjfDgOcR6X4zgAAAD9JyBXmZhWnQ/HeXF5QcuKsrUJAACgFgjIVSbWHA9rbOhyLC5PmzS2bG0CAACoBQJylYmCXKOGN3aG5DiPy3EcAACA/lOkq8pEIa4oyNVRxXpFNnKsijUAAMDGE5CrUIThSw6bsd5x2z8BAAD0n4BcI2z/BAAAsHGsQa4Rtn8CAADYOAJyjbD9EwAAwMYRkGuE7Z8AAAA2jjXINSIKcsWa4/w06+7bPyngBQAA0DMBuQ62f1LACwAAYMME5DrY/qmnAl7Fbg8AAFCPrEGuAwp4AQAAbJiAXAcU8AIAANgwU6zrvICX4l0AAAAdBOQ6LuAVFO8CAADoICDXcQGvM254WPEuAACAv7AGuY4p3gUAAPBXAnIdU7wLAADgrwTkOhbrkKNYVz4kFxbvAgAAqDfWINexUsW7FOgCAADqkYBc54oV7wIAAKhHplgDAACAgAwAAAAdBGQAAAAQkAEAAKCDgAwAAAACMgAAAHQQkAEAAEBABgAAgA7D/3IO61n48sp05byFaUFLa5o2qSkdt8/2afutNit3swAAAAaFgEzJcDxzzvy0Zl17amvPpcdbVqS5j7Skm0/aS0gGAABqkinWFBUjx/lwHOI8LsdxAACAWiQgU1RMq86H47y4vKBlRdnaBAAAMJgEZIqKNcfDGhu6HIvL0yaNLVubAAAABpOATFFRkGvU8MbOkBzncTmOAwAA1CJFuigqCnFFQa6OKtYrspFjVawBAIBaJiBTUoThSw6bUe5mAAAADAlTrAEAAEBABgAAgA4CMgAAAAjIAAAA0EGRLvpl4csr/1LhujXbM1mFawAAoNoJyPQrHM+cMz+tWdee2tpz6fGWFWnuIy3ZtlBCMgAAUK1MsabPYuQ4H45DnMflOA4AAFCtBGT6LKZV58NxXlxe0LKibG0CAADYWAIyfRZrjoc1NnQ5FpenTRpbtjYBAABsLAGZPouCXKOGN3aG5DiPy3EcAACgWinSRZ9FIa4oyNVRxXpFNnKcr2KtujUAAFCtBGT6JULvJYfN6HJMdWsAAKCamWLNgFHdGgAAqGYCMgNGdWsAAKCaCcgMGNWtAQCAaiYgM2BUtwYAAKqZIl0MSXVrAACASicgM+jVrfNsAQUAAFQyAZkhYQsoAACg0lmDzJCwBRQAAFDpBGSGhC2gAACASicgMyRsAQUAAFQ6AZkhYQsoAACgpgPyRRddlBoaGtKpp546cC2ipreAmrXb5DR9cnN2rkAXAABQE1Ws77333nTllVemt7zlLQPbIupyCygAAICqHEFeuXJlOuKII9LVV1+dxo0bN/Ctou62gDrjhofTB+b8NjuPywAAAFUxgjx79ux00EEHpf333z995Stf6fG2a9asyU55ra2t2Xkul8tO5ZRvQ7nbUc8iDB9y2V/3R36ipTX91yMvpp+fOPTTr/UHCukPdKdPUEh/oJD+QCH9oTL19vPoc0C+7rrr0gMPPJBNse6NCy+8MJ177rnrHV++fHnZO008f4yGh1hLzdC79o6n0viR7altRL4v5NKwhvZ07R2Pp1P2/5shbYv+QCH9ge70CQrpDxTSHyikP1Sm/EDtgAbk559/Pp1yyinp1ltvTaNHj+7Vfc4888x02mmndWnYlClTUnNzc2pqakrllA/o0Radtzzua1mdFq3s/oeSXLqvZU32uQwl/YFC+gPd6RMU0h8opD9QSH+oTL39LPoUkO+///60ZMmStNtuu3Uea2trS3fccUe67LLLsqnUw4YN63KfUaNGZadiDayEDpNvRyW0pR7tPKk5LWhZmU2vzostoHae1FSWz0R/oJD+QHf6BIX0BwrpDxTSH+okIO+3337p0Ucf7XLs6KOPTjvttFP6/Oc/v144hg2JfZDnPtLSuQbZ/sgAAEC59Ckgjx07Nk2fPr3LsU033TRtueWW6x2HvuyPfOW8hWlBy4o0bdLYLBzbHxkAAKiafZBhoNgfGQAAqImAfPvttw9MSwAAAKCMjCBT8fskd0y/bk3TJjWZfg0AAAwaAZmKDscz58zvLOD1eMuKrKBXrFkWkgEAgIHWOOCPCAMkRo7z4TjEeVyO4wAAAANNQKZixbTqwv2RQ1yOatcAAAADTUCmYsWa49gXuVBcjq2gAAAABpqATMWKglyjhjd2huQ4j8txHAAAYKAp0kXFikJcUZCro4r1imzkWBVrAABgsAjIVLQIw5ccNqPczQAAAOqAgEzVskcyAAAwkARkqpI9kgEAgIEmIFNzeyTHSLKRZQAAoK8EZGpqj+QHFy0zsgwAAPSLbZ6oqT2S17W1lxxZBgAA6IkRZKpSTJuOkeF8GM7vkTx8WGPRkeXYJioo7AUAAJQiIFNTeyTH5WeXvtolJEd4jusV9gIAAHoiIFNTeySXGlnOh+dS06/ttQwAAAjI1MXIchwvVdgrP/0aAACobwIydTGyHGLNcUyrLjb9GgAAQECmbvQ0/bqjeNczacnSV9KE8ePScfvsYF0yAADUGQGZVO/Tr0MU71q7ri1NHJNL8xetSnMfWax4FwAA1BkBmVTv06/PuOHhbFS5XfEuAACoawIydS9fvKuh4Ji9kwEAoP4IyNS9fPGu/AhysHcyAADUn8ZyNwDKLUaEo1hXhOLQ272TAQCA2mIEmbr31+JdUcV6WXr3+M07q1jbOxkAAOqHgAx/CckXz5qRli9fnpqbm1NDQ8dosr2TAQCgfphiDf2cfg0AANQWI8jQj72T8wW6VLgGAIDaISBDP/ZODipcAwBAbTHFGvpJhWsAAKgtAjL0kwrXAABQWwRk6KdYc5wv3pWnwjUAAFQvARn6SYVrAACoLYp0wSBVuAYAAKqLgAyDUOEaAACoPqZYAwAAgBFkGByxR3LH1OvWrJiXqdcAAFD5BGQYhHA8c878zj2SH29ZkeY+0pKtVxaSAQCgcpliDQMsRo7z4TjEeVyO4/kAfcYND6cPzPltdh6XAQCA8jOCDAMsplXnw3FeXI5K10aXAQCgchlBhgEWa47zeyPnxeXYBqqn0WUjywAAUF4CMgywKMg1anhjZ0iO87gcx0uNLj+4aFk2svyTB15Ij73Qmp3HZSEZAACGjoAMAyymSseU6Vm7TU7TJzdn5/kp1KVGl9e1tfe4bhkAABh81iDDIIgwfMlhM9Y7HqPIseY4H4bzo8vDhzWWXLcMAAAMDSPIUAGjy2/bdvOS65YBAIChYQQZKmB0udTIchwHAACGhoAMFTSyHGuOY1p1jBxHOM5v/RTFujqua83WMRdeBwAADAwBGSp83bK9kwEAYGhYgwwVrqe9kwEAgIEjIEOFK7V3sgrXAAAwsARkqHCl9k5W4RoAAAaWgAwVLgpyRUXrfEhW4RoAAAaHIl1Q5RWuAQCAgSEgQxVXuA62gAIAgIEhIEMVswUUAAAMHGuQoYrZAgoAAAaOgAxVzBZQAAAwcARkqGK2gAIAgIFjDTJUsSjIFWuO89OsC7eAUrwLAAD6RkCGGtwCKijeBQAAfSMgQw1uAXXGDQ+XLN5VarsoAACod9YgQw1SvAsAAPpOQIYapHgXAAD0nSnWUGfFu0KpAl4KewEAUM8EZKij4l35EFysgNc3j9gtnXDtAwp7AQBQtwRkqKPiXSFCc7ECXufe9PseC3sZXQYAoNYJyFBnShXwamldXbKwV6lRZ6PLAADUEkW6oM6UKuA1qWl0ycJepUad4zgAANQKARnqTEyNjoJd+TCcL+B1zsG7FD0et7dtFAAA9cAUa6gzPRXwKnU8Rp1jWnVhSLZtFAAAtUZAhjpUqoBXqeMb2jYKAABqgYAMbFBPo8sAAFArBGSgV0qNLgdbQAEAUAsEZGCj2AIKAIBaISADG6WnLaBiJNnIMgAA1UJABjZKqS2gHly0zMgyAABVxT7IwEaJkeH83sl5cXldW3vJkWUAAKhERpCBjVJqC6jhwxqLjixHFeygsBcAAJVGQAYGZQuouPzs0le7hOQIz3F9T4W9guAMAEA5CMjAoGwBVWpkOR+ei02//uovn0h3PLW05Lplo84AAAwmARkY0pHlOF6qsNddz/65x4rYin4BADCYBGRgSEeWQ4z+RsDtPv06Sn2VWrfc03ZSxZ4DAAD6ShVrYMjFaHBMt85Xv85Pv/7b7bYsWhE7Rp9LjTrni34BAMDGMoIMVMz06zDvqZdLrlsuNuoc9w3WJwMAsLEEZKCipl+XWrfcU9GvnqpiC8kAAAzKFOsLL7ww7b777mns2LFpwoQJ6YMf/GB68skn+/IQAL0KznNP2is7zwfc/KjzrN0mp+mTm7PzfADuaX0yAAAMygjyvHnz0uzZs7OQvG7dunTWWWel973vfWnBggVp00037ctDAQzYqLP1yQAADHlA/sUvftHl8ne/+91sJPn+++9Pe++994A0CKCvSlXFjinaPa1Ntm4ZAIABW4O8fPny7HyLLbYoeZs1a9Zkp7zW1tbsPJfLZadyyreh3O2gMugP1evYvbdL//XIi2nNulzB+uSGdOD0iengOb/tnH79REtrdrufn7hXdr9DLptf9LoIyfoD3ekTFNIfKKQ/UEh/qEy9/Tz6HZDb29vTqaeemvbcc880ffr0Htctn3vuuUXDdbk7TTz/ypUrs383NHTdWob6oz9Ur/EjU/rR0TPST+7/Y3pu6ao0dfwmadbbt0k/uf8PafzI9tQ2Iv+7JpeGNbSna+94/C/3K37dKfv/jf7AevQJCukPFNIfKKQ/VKb8QO2gBeRYi/zYY4+l+fPn93i7M888M5122mldGjZlypTU3NycmpqaUjnlA3q0RedFf6hu8bnt8qZJXY6dOfeZtGhl9z/E5dJ9LTGrJVfyungs/YHu9AkK6Q8U0h8opD9Upt5+Fv0KyCeeeGKaO3duuuOOO9I222zT421HjRqVnYo1sBI6TL4dldAWyk9/qC07T2pOC1pWrrc2eedJHX+cK3Xds0tfTVfOeyYtWfpKmjB+XDpunx2sTSbjdwSF9AcK6Q8U0h/qJCDHX0NOOumkdOONN6bbb789TZ06tb/tAxh0Pe2dHIpd9w+7Tsr2VF67ri1NHJNL8xetSnMfWdy5pZTCXgAAtWt4X6dV/+AHP0g///nPs72QFy9e3Dl9YMyYMYPVRoB+ye+d3BFoV2RVrQsDbbHr8nsqtxfZUzmuj/CcD9VROTtCdj48AwBQRwH5iiuuyM733XffLsevueaadNRRRw1sywAGce/kUtfl91RuKLKncj48txUJz/E4RpcBAKpbn6dYA9TDnsr5EeTCPZXz4TkVCc8Rjo0uAwBUt8ZyNwCgksSob6xFjlAcCtctR3jOH+8ennsaXQYAoDoIyABF1i1/aLfJaYetNsvO86PAPYXnnkaXAQCoDv3eBxmgVkUYvnjWjLR8+fIuexj2VPQrPzW7+7ZRcRtrkwEAqoOADDAARb9KbSmV3zbK2mQAgMonIAMMgFKjy/2tfG3UGQBg6AnIAAOkp22j+lL5+ptH7JZOuPYBo84AAENMkS6AQdSfytfn3vR7FbEBAMpAQAYYRP2pfN3SurrHitgx8nzGDQ+nD8z5bXYelwEA2HimWAMMov5Uvp7UNDoteuW1khWxeyr6Ze0yAED/CcgAFVb5+pyDd+myBrlw1Lmnol9xvYrZAAD9Z4o1QJlHl2ftNjlNn9ycncflfXecUPR43L6nol89hWcAADbMCDJABY4ulzpealp2TL/uKTwDALBhRpABaqToV08VsxX2AgDYMCPIADVS9KvUmuZ/2HWSwl4AAL0gIAPU0LTsYuF5Ywp7Cc8AQD0RkHvwetvr5W4CQyiXy2WfeZwaGrpOU6X+VGt/2HaLUen8Q6d1OfbE4ldSQ2pLw7stqnli8bJ01bynUlv7utSQcp3Xt7W3Z8eP+bup6fAr7uwMz0+/tDz94rEX0vWf2SNtN77+QnK19gkGh/5AIf2Beu4PI4aNSLWkIRef4BBqbW1Nzc3Nafny5ampqSmVU7z0aEe0p1jnvfTXl5alXQAAANXg1P1OTdWgtzlUkS4AAAAwxbpns/edXe4mUEEzCqgvtdYfnl26Mn3nt8+lxxevSDtPHJtNn45p0nG8cBp1vrBXTKOO2//soRfX21LqvTtvlf7n6T8VvU8tT72utT7BxtEfKKQ/UEh/qG4Cch3Np2fDv8ziM4+TX2bUWn/Y8Q3j0sWHjSt6/MbZexetin3sPiPSzY8sSeva/xqEhzU2ppQbnl5dG2uV431pSOvaU3b69m8XZcXDShX2qvaCX7XWJ9g4+gOF9AcK6Q/VTUAGqHN9rYp9ynUPdhlVDnE5bhMhuFhV7G8esVs64doHSlbLBgCoBAIyAH0KzzH6GwG3+9TrCNCltpQ696bfl9xqqqdRZwCAoSQgA9AnEV5j9Lf7GuSeRpdbWlf3edR5Q3sxC9UAwEATkAHok1JTr+N4qdHlSU2j06JXXuvTqHMcj8c1ZRsAGCq2eQKg31Ov5560V3aeD6URaGM0OcJvyI8un3PwLkWPx+1jBLjU6HJ/pmyHGF0+44aH0wfm/DY7j8sAABtiBBmAIRld7uuoc9ymVHju75TtYFo2AFCKgAzAkFXFLna8pzXNEWYHasr2V3/5RLrjqaV9XusMANQPU6wBKKv86PKs3San6ZObs/N8aB3IKdt3PfvnktOy86POP3nghfTYC63ZeVw2NRsA6osRZACqbi/m/kzZjijd17XOXbeheiYtWfpKmjB+XDpunx2MLgNADRKQAaiLKdt/u92W6VcLXurTWufCNc1r17WliWNyaf6iVWnuI4s3uKbZlG0AqD4CMgA1pdTocpj31Mt9WutcuKa5vQ9rmkNPezsDAJVJQAag5pQaXS41LbunQmGnXPdgdqxjtXPv1jSHnvZ2NrIMAJVJQAagbvRnrXN+TXN+BLk3a5pTyhW97sFFy3ocWTYtGwDKS0AGgF6saV67LgJvx+jyyA2saQ7Fpmyva2vvcWS5P+FZqAaAgSMgA0AP/jq6HFWsl6V3j988q2Ld05rmUGzK9vBhjf2qpF0qPH/ziN3SCdc+YK0zAAwQARkANiDC5sWzZqTly5en5ubm1NDQsSK51LTsUtfF5WeXvtrnStqlwvO5N/2+F9tTGV0GgN4SkAFggKdll7qup2JgPVXSLhWeW1pXb3B7qlJVtgVnAFifgAwAQ6SnYmD9Cc+TmkanRa+81uP2VG192J7KmmYA6p2ADAAVXkm7VHg+5+BduqxBLrY9VerD9lQbKhQGALVOQO7BulWryt0EhlAul0ttr72W1o0Y0bm+kPqlP1COPvHGTRvTBf/w5vX+XxTHf/5Pu6Vvz38uPbF4Zdpp4mbpU3tNTVPHb1b0eNx++pYj08Ln1643ujzy9fY0bO3raVi3535q0dJ09a2rU271a2lYhO38617XkK6+9ffpKx/cdVBec7XyO4JC+gP13B+Gb7JJqiUNufgEh1Bra2tW4CQKnTQ1NaVyipfeveBKoZu376hECgAAwPpmLlyYqkFvc2jjkLYKAAAAKpQp1j048NFHy90EhnhGQfxlKf6iVA/TYeiZ/kAt94nnlq4sOl07jh/2rTvT2oI1zSOHN6Ybjt8ju/3PH3pxvSnb22w+Jv1x2fqFwg5569bZ45Z6vNDRhhVpp4ljO9tQLWqpP7Dx9AcK6Q/VTUCuo/n0bPiX2bDXX88+d7/M0B+o5T7x5m03SRd9fELR4z/97H5FC4U99qcn02uNI9ebe/b8a7m0usjxx/70err6npa0Mo1IbY25zuvXpob0tTsWdamk/ejStemmJ17pVSXtSqmyXUv9gY2nP1BIf6huAjIAsMEq2xFG+7rVVKn9m/tbSTvY2xmAwSQgAwAb1J+tpkrt3xzjKcWCc4xal9q/OY4HezsDMJgEZABgg3rap7mv+zf/7XZbpl8teKlPo87x2CnlBnxEWnAGoJCADABs1PTrno4XC89h3lMv92nUOe4bBmpEekOjzgDUJwEZABg0pcJzX0ed88F6oEakexp1jvaalg1QnwRkAGDI9XXUOR9OB2pEuqdR5wjHpaZl/3VN8zNpydJX0oTx49Jx++wgPAPUiIZc1CEfQrEnWHNzc1q+fHm2N1g5xUuPdkR7lGBHf6CQ/kB3+kRl++uIb9dQ3T3s5sPzPn+zVdFR51m7Tc7+/ZMHXih6XX5N89p1bWnimFxa/Frs7TzM1Ow65/cDhfSHytTbHGoEGQCoegO5DvqU6x7c4Jrm9hIFwUzLBqhuAjIAUNP6ug661J7PhWuaC8eE4vKDi5b1Ylq28AxQ6QRkAKAulQrOPRUKy69pzo8gh7h+XVu7raYAaoCADABQoDd7O69dF0G4IzyPHN6Yhg9rtNUUQA0QkAEA+rymOapYL0vvHr95VsU6QvCzS18d0K2mABh6AjIAQB9ESL541owuVWp7My27L1tNAVAejWV6XgCAmpEfWY6toKZPbs7O81OlIyRHWI5QHPLh+W+327LzWPdRZwDKwwgyAECFbTUFQHkIyAAAFbbVFADlISADAFRYcAagPARkAIAKtPDllfZIBhhiAjIAQAWG45lz5pfcI1l4BhgcAjIAQIWJ8Ftqj+QIwz2FZwD6zzZPAAAVJkaGS+2R3FN4BmDjCMgAABUmpk2X2iO5p/AMwMYRkAEAKkxMo449kfMhuXCP5J7CMwAbR0AGAKgwsZY41hTP2m1ymj65OTvPrzHuKTwDsHEU6QIAqKI9kvPhuaOK9Yps5FgVa4CBISADANRIeAZg45hiDQAAAAIyAAAAdDDFGgCgRix8eeVf1ia3ZtWurU0G6BsBGQCgRsLxzDnz05p17dm+yI+3rEhzH2nprH4NwIaZYg0AUANi5DgfjkOcx+U4DkDvCMgAADUgplXnw3FeXI6toADoHQEZAKAGxJrjYY0NXY7F5dgnGYDeEZABAGpAFOQaNbyxMyTHeVyO4wD0jiJdAAA1IApxRUGujirWK7KRY1WsAfpGQAYAqBERhi85bEa5mwFQtUyxBgAAAAEZAAAAOgjIAAAAICADAABABwEZAAAA+huQL7/88vSmN70pjR49Or3rXe9K99xzz8C3DAAAACp5m6cf/ehH6bTTTkvf+ta3snB86aWXpgMOOCA9+eSTacKECYPTSgAANsrCl1f+ZY/k1jRtUlOXPZJLXTdU96n+NjyTlix9JU0YPy4dt88Odfw+aEMl94fBeLxa1JDL5XJ9uUOE4t133z1ddtll2eX29vY0ZcqUdNJJJ6UvfOELG7x/a2tram5uTsuXL09NTU2pnOKlRzuiPQ0NDWVtC+WnP1BIf6A7fYJq7g/xBXfmnPlpzbr21NaeS8MaG9Ko4Y3p5pP2yq4vdt03j9gtnXDtA4N+n1pow9p1bWnimFxa/FpDGjl8WN2+D9pQuf1hMB5v+yoLyb3NoX0aQV67dm26//7705lnntl5rLGxMe2///7pzjvvLHqfNWvWZKfChuX/x9LHbD7g8m0odzuoDPoDhfQHutMnqOb+ECNa8aW9vT2XIs7H+dp1uex4KHbdeTc9NiT3qY02tBdcl+r4fdCGSu0Pg/F4F8+akapJb39f9ykgL126NLW1taU3vOENXY7H5SeeeKLofS688MJ07rnnrnc8knu5/6cSz79y5crs39Xw118Gl/5AIf2B7vQJqrk/xHTPGNHqKpeWLF2WnRe7rmHda0Nyn1poQ/SACWPy19Tv+6ANldsfBuPxli9fnqpJfqB2wNcg91WMNsea5cKGxZTsGN6uhCnWoVqmRzG49AcK6Q90p09Qzf0h1kLOX7QqmyKZF1Ml3z1+8+zfxa7bdtyYtHjFa4N+n1poQ4yqhRdejdmV9fs+aEPl9ofBeLzm5uZUTXr7u7pPVazHjx+fhg0bll566aUux+PyxIkTi95n1KhRWRAuPOUb6OTk5OTk5OTkNPinKBQUayFjaVyMb8V5XI7jpa47++DpQ3IfbdAGbajOx2uogN9tfT0NWpGud77znWnOnDnZ5Zhjv+2226YTTzxRkS6qmv5AIf2B7vQJqr0//LUS7Yo0bdLYElVqu143VPep/jZE1eJlacL4zYtULa6n90EbKrk/DMbjVZPe5tA+B+TY5unII49MV155ZRaUY5unH//4x9ka5O5rkzemYUOhGv/nxuDRHyikP9CdPkEh/YFC+gOF9Ic6qmIdPvKRj6SXX345nX322Wnx4sXprW99a/rFL37Rq3AMAAAAlapfRbpiOnWcAAAAoFb0qUgXAAAA1CoBGQAAAARkAAAA6CAgAwAAgIAMAAAAHQRkAAAAEJABAACgg4AMAAAAAjIAAAB0EJABAABAQAYAAIAOAjIAAAAIyAAAANBheBpiuVwuO29tbR3qpy7almhHQ0NDdqK+6Q8U0h/oTp+gkP5AIf2BQvpDZcrnz3werZiAvGLFiux8ypQpQ/3UAAAA1LEVK1ak5ubmktc35DYUoQdYe3t7evHFF9PYsWPL/heV+CtCBPXnn38+NTU1lbUtlJ/+QCH9ge70CQrpDxTSHyikP1SmiL0RjrfeeuvU2NhYOSPI0ZhtttkmVZLouDovefoDhfQHutMnKKQ/UEh/oJD+UHl6GjnOU6QLAAAABGQAAADoUNcBedSoUemcc87JzkF/oJD+QHf6BIX0BwrpDxTSH6rbkBfpAgAAgEpU1yPIAAAAkCcgAwAAgIAMAAAAHQRkAAAAqPeAfPnll6c3velNafTo0eld73pXuueee8rdJIbAhRdemHbfffc0duzYNGHChPTBD34wPfnkk11us3r16jR79uy05ZZbps022yzNmjUrvfTSS2VrM0PjoosuSg0NDenUU0/tPKYv1J8XXnghfeITn8g+8zFjxqRdd9013XfffZ3XR23Ls88+O02aNCm7fv/9909PP/10WdvM4Ghra0tf/vKX09SpU7PPevvtt0//8i//kvWBPP2hdt1xxx1p5syZaeutt87+3/Czn/2sy/W9+ez//Oc/pyOOOCI1NTWlzTffPH3qU59KK1euHOJXwmD3h9dffz19/vOfz/5/semmm2a3+eQnP5lefPHFLo+hP1SHug3IP/rRj9Jpp52WlWB/4IEH0owZM9IBBxyQlixZUu6mMcjmzZuXBZ677ror3Xrrrdkvtfe9733p1Vdf7bzNZz/72XTzzTen66+/Prt9/IL70Ic+VNZ2M7juvffedOWVV6a3vOUtXY7rC/XllVdeSXvuuWcaMWJEuuWWW9KCBQvS1772tTRu3LjO21xyySXpG9/4RvrWt76V7r777uzLUPz/I/6YQm25+OKL0xVXXJEuu+yy9Pjjj2eX4/OfM2dO5230h9oV3wvi+2EMqBTTm88+wtDvf//77PvG3Llzs5B17LHHDuGrYCj6w6pVq7I8EX9Qi/Of/vSn2eDLwQcf3OV2+kOVyNWpd77znbnZs2d3Xm5ra8ttvfXWuQsvvLCs7WLoLVmyJIYCcvPmzcsuL1u2LDdixIjc9ddf33mbxx9/PLvNnXfeWcaWMlhWrFiRe/Ob35y79dZbc/vss0/ulFNOyY7rC/Xn85//fG6vvfYqeX17e3tu4sSJua9+9audx6KfjBo1KvfDH/5wiFrJUDnooINyxxxzTJdjH/rQh3JHHHFE9m/9oX7E7/0bb7yx83JvPvsFCxZk97v33ns7b3PLLbfkGhoaci+88MIQvwIGsz8Uc88992S3+8Mf/pBd1h+qR12OIK9duzbdf//92VSYvMbGxuzynXfeWda2MfSWL1+enW+xxRbZefSNGFUu7B877bRT2nbbbfWPGhUzCg466KAun3nQF+rPTTfdlN7xjnekww8/PFuC8ba3vS1dffXVndc/99xzafHixV36RHNzc7ZMR5+oPe9+97vTr3/96/TUU09llx9++OE0f/78dOCBB2aX9Yf61ZvPPs5jGm38TsmL28d3zhhxpva/X8ZU7OgDQX+oHsNTHVq6dGm2rugNb3hDl+Nx+Yknnihbuxh67e3t2XrTmFI5ffr07Fj8D2/kyJGdv9AK+0dcR2257rrrsulQMcW6O32h/jz77LPZlNpYgnPWWWdl/eLkk0/O+sGRRx7Z+bkX+/+HPlF7vvCFL6TW1tbsD2PDhg3Lvjucf/752TTJoD/Ur9589nEef2grNHz48OwP8vpHbYtp9rEm+WMf+1i23jjoD9WjLgMyFI4cPvbYY9mIAPXn+eefT6ecckq2FiiK9UH80Sz+un/BBRdkl2MEOX5HxBrDCMjUlx//+Mfp2muvTT/4wQ/SLrvskh566KHsj6pRgEd/AIqJmWcf/vCHsyJu8QdXqk9dTrEeP3589pfg7pVo4/LEiRPL1i6G1oknnpgVSPjNb36Tttlmm87j0QdiGv6yZcu63F7/qD0xhToK8+22227ZX3HjFIW4ouhK/DtGAvSF+hLVaKdNm9bl2M4775wWLVqU/Tv/ufv/R304/fTTs1Hkj370o1l12n/8x3/MCvfFbghBf6hfvfns47x78dd169ZllYz1j9oOx3/4wx+yP77nR4+D/lA96jIgx1S5t7/97dm6osJRg7i8xx57lLVtDL74i16E4xtvvDH993//d7Z9R6HoG1HBtrB/RCXC+IKsf9SW/fbbLz366KPZqFD+FKOHMX0y/299ob7Ecovu277F+tM3vvGN2b/j90V8kSnsEzEFN9aP6RO1JyrTxvrAQvEH9vjOEPSH+tWbzz7O4w+s8cfYvPjeEf0n1ipTm+E4tvq67bbbsq0CC+kPVaTcVcLK5brrrssqDX73u9/Nqsode+yxuc033zy3ePHicjeNQfaZz3wm19zcnLv99ttzLS0tnadVq1Z13ub444/Pbbvttrn//u//zt133325PfbYIztR+wqrWAd9ob5E1dHhw4fnzj///NzTTz+du/baa3ObbLJJ7vvf/37nbS666KLs/xc///nPc4888kjukEMOyU2dOjX32muvlbXtDLwjjzwyN3ny5NzcuXNzzz33XO6nP/1pbvz48bkzzjij8zb6Q23vcPDggw9mp/jK/K//+q/Zv/NViXvz2b///e/Pve1tb8vdfffdufnz52c7JnzsYx8r46tiMPrD2rVrcwcffHBum222yT300ENdvl+uWbOm8zH0h+pQtwE5zJkzJ/viO3LkyGzbp7vuuqvcTWIIxC+1Yqdrrrmm8zbxP7cTTjghN27cuOzL8aGHHpr9kqP+ArK+UH9uvvnm3PTp07M/ou600065q666qsv1sb3Ll7/85dwb3vCG7Db77bdf7sknnyxbexk8ra2t2e+D+K4wevTo3HbbbZf74he/2OULr/5Qu37zm98U/b4Qfzjp7Wf/pz/9KQtAm222Wa6pqSl39NFHZ0GL2uoP8Qe0Ut8v4355+kN1aIj/lHsUGwAAAMqtLtcgAwAAQHcCMgAAAAjIAAAA0EFABgAAAAEZAAAAOgjIAAAAICADAABABwEZAAAABGQAAADoICADAACAgAwAAAAdBGQAAAASKf1/PpHxBtG8sykAAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Suggested number of sources to localize: 69\n",
      "Suggested rank is: 50\n"
     ]
    }
   ],
   "source": [
    "# Suggest number of sources to localize\n",
    "# and optimization parameter to use for both localization and reconstruction\n",
    "sugg_n_sources, sugg_rank = localizer.suggest_n_sources_and_rank(\n",
    "    R=data_cov_task.data,\n",
    "    N=noise_cov.data,\n",
    "    show_plot=True,\n",
    "    subject=subject,\n",
    "    s=14\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T21:36:23.456621Z",
     "start_time": "2025-09-11T21:36:23.350760Z"
    }
   },
   "id": "9fc08f372dc5e18d"
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001B[36mCalculating activity index for localizer: mpz_mvp\u001B[0m\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5124/5124 [00:22<00:00, 227.29it/s]\n",
      "100%|██████████| 5124/5124 [00:23<00:00, 219.92it/s]\n",
      "100%|██████████| 5124/5124 [00:24<00:00, 206.87it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 230.55it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 239.85it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 232.09it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 241.26it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 234.91it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 242.70it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 243.91it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 234.48it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 241.44it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 239.30it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 248.39it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 246.96it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 243.96it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 250.55it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 252.83it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 250.98it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 250.03it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 249.71it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 252.17it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 245.33it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 252.36it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 247.02it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 251.69it/s]\n",
      "100%|██████████| 5124/5124 [00:23<00:00, 215.61it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 234.60it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 243.44it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 245.63it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 240.35it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 247.86it/s]\n",
      "100%|██████████| 5124/5124 [00:20<00:00, 246.37it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 243.88it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 243.70it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 243.27it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 238.34it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 239.47it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 238.29it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 242.33it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 235.48it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 228.16it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 225.14it/s]\n",
      "100%|██████████| 5124/5124 [00:21<00:00, 233.82it/s]\n",
      "100%|██████████| 5124/5124 [00:24<00:00, 206.37it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 227.36it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 224.77it/s]\n",
      "100%|██████████| 5124/5124 [00:27<00:00, 185.77it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 226.67it/s]\n",
      "100%|██████████| 5124/5124 [00:22<00:00, 231.82it/s]\n",
      "100%|██████████| 5124/5124 [00:30<00:00, 165.56it/s]\n",
      "100%|██████████| 5124/5124 [00:33<00:00, 154.77it/s]\n",
      "100%|██████████| 5124/5124 [00:31<00:00, 163.60it/s]\n",
      "100%|██████████| 5124/5124 [00:31<00:00, 165.03it/s]\n",
      "100%|██████████| 5124/5124 [00:31<00:00, 162.19it/s]\n",
      "100%|██████████| 5124/5124 [00:32<00:00, 159.90it/s]\n",
      "100%|██████████| 5124/5124 [00:33<00:00, 154.72it/s]\n",
      "100%|██████████| 5124/5124 [00:33<00:00, 151.46it/s]\n",
      "100%|██████████| 5124/5124 [00:34<00:00, 148.82it/s]\n",
      "100%|██████████| 5124/5124 [00:34<00:00, 148.85it/s]\n",
      "100%|██████████| 5124/5124 [00:35<00:00, 142.92it/s]\n",
      "100%|██████████| 5124/5124 [00:37<00:00, 137.03it/s]\n",
      "100%|██████████| 5124/5124 [00:39<00:00, 128.70it/s]\n",
      "100%|██████████| 5124/5124 [00:40<00:00, 127.15it/s]\n",
      "100%|██████████| 5124/5124 [00:41<00:00, 123.43it/s]\n",
      "100%|██████████| 5124/5124 [00:42<00:00, 120.62it/s]\n",
      "100%|██████████| 5124/5124 [00:42<00:00, 120.70it/s]\n",
      "100%|██████████| 5124/5124 [00:42<00:00, 119.88it/s]\n",
      "100%|██████████| 5124/5124 [00:43<00:00, 117.34it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Leadfield indices corresponding to localized sources: [1689, 1291, 1489, 2159, 4405, 5051, 2039, 800, 2020, 3224, 1792, 4998, 821, 1865, 3762, 4165, 2372, 2471, 4817, 495, 1755, 2548, 3334, 1572, 2322, 2008, 4379, 3733, 4426, 2527, 2230, 2362, 4774, 356, 2545, 902, 108, 3720, 211, 481, 64, 2740, 2303, 4641, 148, 170, 1460, 3454, 68, 4850, 1732, 641, 2310, 4678, 764, 669, 1376, 3730, 2856, 4542, 4680, 2932, 3896, 600, 3172, 2564, 2631, 4982, 1342]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Localize\n",
    "locs_task = localizer.localize(\n",
    "    subject=subject,\n",
    "    subjects_dir=subjects_dir,\n",
    "    localizer_to_use=[\"mpz_mvp\"],\n",
    "    n_sources_to_localize=sugg_n_sources,\n",
    "    R=data_cov_task.data,\n",
    "    N=noise_cov.data,\n",
    "    forward=fwd,\n",
    "    r=sugg_rank,\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T22:05:58.982307Z",
     "start_time": "2025-09-11T21:36:23.455793Z"
    }
   },
   "id": "8000a464fcae6679"
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [],
   "source": [
    "# Transform leadfield indices to vertices\n",
    "lh_vert_task, lh_idx_task, rh_vert_task, rh_idx_task = utils.transform_leadfield_indices_to_vertices(\n",
    "    lf_idx=locs_task[\"sources\"],\n",
    "    src=src,\n",
    "    hemi=\"both\",\n",
    "    include_mapping=True\n",
    ")\n",
    "\n",
    "locs_task.add_vertices_info(\n",
    "    lh_vertices=lh_vert_task,\n",
    "    lh_indices=lh_idx_task,\n",
    "    rh_vertices=rh_vert_task,\n",
    "    rh_indices=rh_idx_task\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T22:05:59.014700Z",
     "start_time": "2025-09-11T22:05:58.985058Z"
    }
   },
   "id": "cd9cd7db0a043160"
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [],
   "source": [
    "new_fwd_task = utils.subset_forward(\n",
    "    old_fwd=fwd,\n",
    "    localized=locs_task,\n",
    "    hemi=\"both\"\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T22:05:59.022372Z",
     "start_time": "2025-09-11T22:05:58.986535Z"
    }
   },
   "id": "f068189124c7cd57"
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing rank from covariance with rank=None\n",
      "    Using tolerance 5.2e-13 (2.2e-16 eps * 128 dim * 18  max singular value)\n",
      "    Estimated rank (eeg): 86\n",
      "    EEG: rank 86 computed from 128 data channels with 1 projector\n",
      "Computing rank from covariance with rank=None\n",
      "    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)\n",
      "    Estimated rank (eeg): 86\n",
      "    EEG: rank 86 computed from 128 data channels with 1 projector\n",
      "Making MVP_R beamformer with rank {'eeg': 86} (note: MNE-Python rank)\n",
      "Computing inverse operator with 128 channels.\n",
      "    128 out of 128 channels remain after picking\n",
      "Selected 128 channels\n",
      "Whitening the forward solution.\n",
      "    Created an SSP operator (subspace dimension = 1)\n",
      "Computing rank from covariance with rank={'eeg': 86}\n",
      "    Setting small EEG eigenvalues to zero (without PCA)\n",
      "Creating the source covariance matrix\n",
      "Adjusting source covariance matrix.\n",
      "Computing beamformer filters for 69 sources\n",
      "MVP_R computation - in progress...\n",
      "Filter rank: 50\n",
      "Filter computation complete\n"
     ]
    }
   ],
   "source": [
    "mcmv_params = {\n",
    "    'filter_rank': sugg_rank,\n",
    "    \"filter_type\": \"MVP_R\"\n",
    "}\n",
    "\n",
    "filter_task = beamformer.make_filter(\n",
    "    signal_task.info,\n",
    "    new_fwd_task,\n",
    "    data_cov_task,\n",
    "    reg=0.05,\n",
    "    noise_cov=noise_cov,\n",
    "    pick_ori=None,  # because scalar forward\n",
    "    weight_norm=None,\n",
    "    rank=None,\n",
    "    mvpure_params=mcmv_params\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T22:05:59.086012Z",
     "start_time": "2025-09-11T22:05:59.002464Z"
    }
   },
   "id": "de8ad11b531964"
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using control points [2.84206626e-09 3.11743662e-09 5.08812233e-09]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<mne.viz._brain._brain.Brain at 0x16c554cd0>"
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stc_task = beamformer.apply_filter(signal_task, filter_task)\n",
    "\n",
    "# Add source estimate to mvpure_py.Localized object\n",
    "locs_task.add_stc(stc_task)\n",
    "\n",
    "# Plot\n",
    "viz.plot_sources_with_activity(\n",
    "    subject=subject,\n",
    "    stc=stc_task,\n",
    ")"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2025-09-11T22:06:01.164271Z",
     "start_time": "2025-09-11T22:05:59.087331Z"
    }
   },
   "id": "80e9e875ebb67275"
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   },
   "id": "a4f53f92586fdefe"
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "mvpure-tools-venv",
   "language": "python",
   "display_name": "mvpure-tools-venv"
  },
  "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.10.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}