2_MOU-EC_ClassificationTasks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Classification to identify tasks\n",
    "\n",
    "**Created by: _Matthieu Gilson & Gorka Zamora-López_**\n",
    "\n",
    "The Human Brain Project SGA0 - SGA2 (SP4); SGA3 (WP1)  /  Pompeu Fabra University, Barcelona.\n",
    "\n",
    "<img src=\"Logo_HBP_small.png\" style=\"float:left; width:254px; heigth:50px \" alt=\"HBP logo\" />\n",
    "<img src=\"Logo_UPF_small.png\" style=\"float:right; width:149px; height:50px\"  alt=\"UPF logo\" />\n",
    "\n",
    "\n",
    "<br/>\n",
    "\n",
    "\n",
    "-----------------\n",
    "The purpose of this notebook is to showcase, test and validate the integration of effective connectivity tools in EBRAINS and for use in broader workflows. This notebook should:\n",
    "\n",
    "- Retrieve and load the results of whole-brain EC, estimated previously along the workflow.\n",
    "- Use machine learning tools that are integrated in EBRAINS (Scikit-learn) to perform classification of the data.\n",
    "- Save / export the results and use EBRAINS visualization tools to create fancier outputs, e.g. a brain visualiser.\n",
    "\n",
    "-----------------\n",
    "\n",
    "Here, we will load whole-brain effective previously estimated (see for example the notebook '1_Mou-EC_estimation.ipynb'. The goal is to compare the connectivity measures (in particular, effective versus functional connectivity) in identifying tasks performed by subjects in the scanner. \n",
    "\n",
    "\n",
    "#### References and documentation\n",
    "\n",
    "1. M Gilson, G. Deco, K. Friston et al. (2018) \"Effective connectivity inferred from fMRI transition dynamics during movie viewing points to a balanced reconfiguration of cortical interactions.\" NeuroImage 180, 534-546.\n",
    "DOI: https://doi.org/10.1016/j.neuroimage.2017.09.061\n",
    "2. Gilson M, Zamora-López G, Pallarés V, et al. (2020) \"Model-based whole-brain effective connectivity to study distributed cognition in health and disease\", Network Neurosci. 4(2) 338-373.\n",
    "DOI: https://doi.org/10.1162/netn_a_00117\n",
    "\n",
    "\n",
    "#### Software\n",
    "This Notebook employs Scikit-learn for the classification via machine-learning methods\n",
    "Web: https://scikit-learn.org/stable/\n",
    "GitHub: https://github.com/scikit-learn/scikit-learn\n",
    "\n",
    "\n",
    "#### Data\n",
    "We employ structural connectivity and resting-state fMRI from published data. It consists of recordings for a cohort of 22 healthy participants udergoing five fMRI sessions each: two sessions in resting-state (R1, R2) and three sessions while watching three distinct segments of a movie (M1, M2, M3).\n",
    "\n",
    "1. J. Hlinka, M. Palus, M. Vejmelka, et al. (2011) \"Functional connectivity in resting-state fMRI: is linear correlation sufficient?\" NeuroImage 54, 2218–2225.\n",
    "DOI: https://doi.org/10.1016/j.neuroimage.2010.08.042\n",
    "2. D. Mantini, U. Hasson, V. Betti, et al. (2012) \"Interspecies activity correlations reveal functional correspondence between monkey and human brain areas.\" Nat. Methods 9, 277–282.\n",
    "DOI: https://doi.org/10.1038/nmeth.1868\n",
    "\n",
    "\n",
    "_**TO-DO**: Reproduce this notebook using a dataset accessible in EBRAINS. This \"movie\" dataset is not part of the EBRAINS catalog. At least not yet. See file \"0\\_README.md\" for suggestions._\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------\n",
    "### TO-DO list towards full EBRAINS integration \n",
    "Missing technical points, for full integration of the workflows in this notebook.\n",
    "\n",
    "- Input data is the result of whole-brain EC estimation, which can simply be stored in Drive and loaded here. But… any fancier interface shall be built?\n",
    "- Interface to visualization services. Visualise results in a \"standard brain map\".\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----------\n",
    "### Import the dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import standard library modules\n",
    "import os\n",
    "# Import 3rd party packages\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import scipy.stats as stt\n",
    "\n",
    "import sklearn.linear_model as skllm\n",
    "import sklearn.neighbors as sklnn\n",
    "import sklearn.discriminant_analysis as skda\n",
    "import sklearn.preprocessing as skppc\n",
    "import sklearn.pipeline as skppl\n",
    "import sklearn.metrics as skm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and prepare the data\n",
    "\n",
    "The data required in this notebook is the output of the MOU-EC estimation method previously calculated and saved in, for example, a notebook like '1_MOU-EC_Estimation.ipynb'\n",
    "\n",
    "So far, we will load datafiles stored in the local Drive of this collab: _'.../Effective Connectivity estimation/Datasets/Movie/'_ . The required data files are:\n",
    "\n",
    "- 'ECdata_Movie.npz': zipped file containing various results from the MOU-EC estimation.\n",
    "- 'BOLD_Movie_filtered.npy': BOLD-fMRI (all sessions and subjects) needed to compute the 0-lag covariances\n",
    "\n",
    "We load the data and calculate the features: 'EC' stands for effective connectivity (previously in the workflow, see for example the *1_MOU-EC_Estimation.ipynb* notebook), 'FC' stands for the BOLD correlations, 'FC + mask' for the BOLD correlations only concerning the ROI connections that exist according to the structural data (diffusion tensor imaging), 'PC' for the partial correlations.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/mnt/user/drive/Shared with groups/Effective Connectivity estimation/Jupyter_Notebooks'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of subjects: 22\n",
      "Number of sessions: 5\n",
      "Number of ROIs: 66\n"
     ]
    }
   ],
   "source": [
    "## Load the data needed for the fMRI data from the Drive\n",
    "datapath = '/mnt/user/drive/Shared with groups/Effective Connectivity estimation/Datasets/Movie/'\n",
    "\n",
    "# 1) Load data previously generated, e.g., in a notebook like '1_MOU-EC_Estimation.ipynb'\n",
    "ec_data = np.load( datapath + 'ECdata_Movie.npz' )\n",
    "\n",
    "# Effective connectivity matrices\n",
    "EC = ec_data['EC_mod']\n",
    "# Mask of existing connections in SC (structural connectivity)\n",
    "mask_EC = ec_data['mask_EC']\n",
    "# Extract the number of subjects, sessions and nodes\n",
    "n_sub, n_run, N, N = np.shape(EC)\n",
    "\n",
    "\n",
    "# 2) Load the BOLD signals and calculate the 0-lag covariance matrices\n",
    "ts_emp = np.load( datapath + 'BOLD_movie_filtered.npy' )\n",
    "\n",
    "FC0 = np.zeros([n_sub,n_run,N,N], np.float)\n",
    "for i_sub in range(n_sub):\n",
    "    for i_run in range(n_run):\n",
    "        FC0[i_sub,i_run] = np.cov(ts_emp[i_sub,i_run])\n",
    "\n",
    "\n",
    "# Generate a triangular mask to retain half of the matrix elemenets \n",
    "mask_tri = np.tri(N,N, k=-1, dtype=np.bool)\n",
    "\n",
    "# Print some feedback\n",
    "print( 'Number of subjects:', n_sub )\n",
    "print( 'Number of sessions:', n_run )\n",
    "print( 'Number of ROIs:', N )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Toggle to True to create directory and store results there\n",
    "save_outputs = False\n",
    "if save_outputs:\n",
    "    import os\n",
    "    res_dir = datapath + 'classif/'\n",
    "    if not os.path.exists(res_dir):\n",
    "        os.mkdir(res_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Perform the task-based classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done\n"
     ]
    }
   ],
   "source": [
    "# Calculate features based on EC (effective connectivity)\n",
    "for i_sub in range(n_sub):\n",
    "    for i_run in range(n_run):\n",
    "        # apply z-scoring to EC\n",
    "        EC[i_sub,i_run,mask_EC] = stt.zscore(EC[i_sub,i_run,mask_EC])\n",
    "\n",
    "# Calculate the features based on FC (functional connectivity)\n",
    "corr = np.copy(FC0)\n",
    "for i_sub in range(n_sub):\n",
    "    for i_run in range(n_run):\n",
    "        corr[i_sub,i_run,:,:] /= np.sqrt( np.outer(corr[i_sub,i_run,:,:].diagonal(), \n",
    "                                                   corr[i_sub,i_run,:,:].diagonal()) )\n",
    "\n",
    "# Calculate the features based on PC (partial correlations)\n",
    "PC = np.copy(FC0)\n",
    "for i_sub in range(n_sub):\n",
    "    for i_run in range(n_run):\n",
    "        PC[i_sub,i_run,:,:] = -np.linalg.pinv(PC[i_sub,i_run,:,:])\n",
    "        PC[i_sub,i_run,:,:] /= np.sqrt( np.outer(PC[i_sub,i_run,:,:].diagonal(), \n",
    "                                                 PC[i_sub,i_run,:,:].diagonal()) )\n",
    "\n",
    "print( 'Done' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code performs the task classification by splitting the fMRI sessions into a train set and a test set (`n_rep` times). Here the 22 subjects are split in 80% for training and 20% for testing. We compare 3 classifiers:\n",
    "\n",
    "- MLR is the multinomial logistic regression;\n",
    "- 1NN is the 1-nearest-neighbor;\n",
    "- LDA performs the linear discriminant analysis.\n",
    "\n",
    "See the [scikit-learn](https://scikit-learn.org/stable/) webpage for details.\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running classification: tasks ...\n",
      "  Repetition 0\n",
      "  Repetition 1\n",
      "  Repetition 2\n",
      "  Repetition 3\n",
      "  Repetition 4\n",
      "  Repetition 5\n",
      "  Repetition 6\n",
      "  Repetition 7\n",
      "  Repetition 8\n",
      "  Repetition 9\n",
      "  Repetition 10\n",
      "  Repetition 11\n",
      "  Repetition 12\n",
      "  Repetition 13\n",
      "  Repetition 14\n",
      "  Repetition 15\n",
      "  Repetition 16\n",
      "  Repetition 17\n",
      "  Repetition 18\n",
      "  Repetition 19\n",
      "  Repetition 20\n",
      "  Repetition 21\n",
      "  Repetition 22\n",
      "  Repetition 23\n",
      "  Repetition 24\n",
      "  Repetition 25\n",
      "  Repetition 26\n",
      "  Repetition 27\n",
      "  Repetition 28\n",
      "  Repetition 29\n",
      "  Repetition 30\n",
      "  Repetition 31\n",
      "  Repetition 32\n",
      "  Repetition 33\n",
      "  Repetition 34\n",
      "  Repetition 35\n",
      "  Repetition 36\n",
      "  Repetition 37\n",
      "  Repetition 38\n",
      "  Repetition 39\n",
      "\n",
      "Running classification: rest vs. movie ...\n",
      "  Repetition 0\n",
      "  Repetition 1\n",
      "  Repetition 2\n",
      "  Repetition 3\n",
      "  Repetition 4\n",
      "  Repetition 5\n",
      "  Repetition 6\n",
      "  Repetition 7\n",
      "  Repetition 8\n",
      "  Repetition 9\n",
      "  Repetition 10\n",
      "  Repetition 11\n",
      "  Repetition 12\n",
      "  Repetition 13\n",
      "  Repetition 14\n",
      "  Repetition 15\n",
      "  Repetition 16\n",
      "  Repetition 17\n",
      "  Repetition 18\n",
      "  Repetition 19\n",
      "  Repetition 20\n",
      "  Repetition 21\n",
      "  Repetition 22\n",
      "  Repetition 23\n",
      "  Repetition 24\n",
      "  Repetition 25\n",
      "  Repetition 26\n",
      "  Repetition 27\n",
      "  Repetition 28\n",
      "  Repetition 29\n",
      "  Repetition 30\n",
      "  Repetition 31\n",
      "  Repetition 32\n",
      "  Repetition 33\n",
      "  Repetition 34\n",
      "  Repetition 35\n",
      "  Repetition 36\n",
      "  Repetition 37\n",
      "  Repetition 38\n",
      "  Repetition 39\n",
      "\n",
      "Finished.\n"
     ]
    }
   ],
   "source": [
    "# WARNING: Running this cell will take a few minutes\n",
    "# Labels of sessions for classification\n",
    "n_conn = 4\n",
    "\n",
    "# distinct labels for movie runs\n",
    "task4_labels = np.repeat( np.array([0,0,1,2,3], \n",
    "                            dtype=np.int).reshape([1,-1]), \n",
    "                            n_sub, axis=0 )\n",
    "# rest versus movie, same labels for movie runs\n",
    "RM_labels = np.repeat( np.array([0,0,1,1,1], \n",
    "                            dtype=np.int).reshape([1,-1]), \n",
    "                            n_sub, axis=0 ) \n",
    "\n",
    "# number of repetitions and storage of results\n",
    "n_rep = 40\n",
    "# first index: task4_labels/RM_labels; last index: MRL/1NN/LDA\n",
    "perf = np.zeros([2,n_rep,n_conn,3])\n",
    "\n",
    "# Loop over classification type\n",
    "for i_classif in range(2):\n",
    "\n",
    "    if i_classif==0:\n",
    "        # discriminate rest and 3 movie sessions individually\n",
    "        lab_tmp = task4_labels\n",
    "        print( 'Running classification: tasks ...' )\n",
    "    else:\n",
    "        # discriminate rest vs movie\n",
    "        lab_tmp = RM_labels\n",
    "        print( '\\nRunning classification: rest vs. movie ...' )\n",
    "\n",
    "    # Repeat classification\n",
    "    for i_rep in range(n_rep):\n",
    "        print( '  Repetition %d' %i_rep )\n",
    "        \n",
    "        # Split run indices in train and test sets \n",
    "        # (choose some subjects for train and rest for test)\n",
    "        train_ind = np.ones([n_sub,n_run], dtype=bool)\n",
    "        while train_ind[:,0].sum() >= 0.8 * n_sub:\n",
    "            train_ind[np.random.randint(n_sub),:] = False\n",
    "        test_ind = np.logical_not(train_ind)        \n",
    "\n",
    "        # Loop over connectivity measures\n",
    "        for i_conn in range(n_conn):\n",
    "\n",
    "            if i_conn==0:\n",
    "                # vectorized EC matrices (only retaining existing connections)\n",
    "                vect_features = EC[:,:,mask_EC]\n",
    "            elif i_conn==1:\n",
    "                # vectorized FC matrices (only retaining low triangle)\n",
    "                vect_features = corr[:,:,mask_tri]\n",
    "            elif i_conn==2:\n",
    "                # vectorized FC matrices (only retaining SC existing connections)\n",
    "                vect_features = corr[:,:,mask_EC]\n",
    "            else:\n",
    "                # vectorized PCmatrices (only retaining low triangle)\n",
    "                vect_features = PC[:,:,mask_tri]\n",
    "\n",
    "            # dimension of vectorized EC\n",
    "            dim_feature = vect_features.shape[2]\n",
    "\n",
    "            # Number of components for LDA\n",
    "            if i_classif==0:\n",
    "                n_comp_LDA = 3\n",
    "            else:\n",
    "                n_comp_LDA = 1\n",
    "\n",
    "            # Classifier and learning parameters\n",
    "            c_MLR = skppl.make_pipeline( skppc.StandardScaler(), \n",
    "                                        skllm.LogisticRegression(C=0.1, \n",
    "                                                                 penalty='l2', \n",
    "                                                                 multi_class='multinomial', \n",
    "                                                                 solver='lbfgs', \n",
    "                                                                 max_iter=500) )\n",
    "            c_1NN = sklnn.KNeighborsClassifier(n_neighbors=1, \n",
    "                                               algorithm='brute', \n",
    "                                               metric='correlation')\n",
    "            c_LDA = skda.LinearDiscriminantAnalysis(n_components=n_comp_LDA, \n",
    "                                                    solver='eigen', \n",
    "                                                    shrinkage='auto')\n",
    "            \n",
    "            # Train and test classifiers with subject labels\n",
    "            c_MLR.fit(vect_features[train_ind,:], lab_tmp[train_ind])\n",
    "            perf[i_classif, i_rep,i_conn,0] = \\\n",
    "                c_MLR.score( vect_features[test_ind,:], lab_tmp[test_ind] )\n",
    "        \n",
    "            c_1NN.fit(vect_features[train_ind,:], lab_tmp[train_ind])\n",
    "            perf[i_classif, i_rep,i_conn,1] = \\\n",
    "                c_1NN.score(vect_features[test_ind,:], lab_tmp[test_ind])\n",
    "        \n",
    "            c_LDA.fit(vect_features[train_ind,:], lab_tmp[train_ind])\n",
    "            perf[i_classif, i_rep,i_conn,2] = \\\n",
    "                c_LDA.score(vect_features[test_ind,:], lab_tmp[test_ind])\n",
    "\n",
    "\n",
    "print( '\\nFinished.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Saving the results into files ... \n"
     ]
    }
   ],
   "source": [
    "# Save the results and the metadata into a file for later use\n",
    "res_dir = datapath + 'classif/'\n",
    "if not os.path.exists(res_dir):\n",
    "    os.mkdir(res_dir)\n",
    "\n",
    "np.save( res_dir + 'perf_tasks.npy', perf )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The performance is interpreted by comparing the classification accuracy with the chance level (the dashed line), in relation with the gap towards 100% correct. Here we can calculate the chance level straight from the session properties.\n",
    "\n",
    "<br>\n",
    "\n",
    "#### Plot the task-classification results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEZCAYAAACU3p4jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZxcVZ338c+veu900umskIUkAgkBlIwEfBAVNDoCwjjKJggSHhWQcWEUUVAfxQVmECQgiPCghGVUNkdkCTCAgAiKiUQWNcgSSBDJ1mbrves3f5zbSXX1Wt23lu77fb9e9UrXufd2nTqp/tapc889Ze6OiIiMfqliV0BERApDgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgmhwBdJODNbamZ3Fbsekn+mefhSiszMgWPd/bZi12W0M7N6Qhb8o9h1kfwqL3YFZHQxs0p3byt2PYbDzCrcvb3Y9SgUd99c7DpIYWhIR4bFzB42s6vM7GIzWw/8JiqvN7NrzGydmW01s0fMbGHGcfVmdmO0vcXMXjKzs6Jtq6PdbjUzz7if/dg/MbPbs8pSZrbGzD4f3TczO8fMXjSzZjN7xsxOyth/dvQYJ5jZQ2bWDJzeX/2i49zMjsl67NVmdnbG/dPN7Pno+A1mdp+Z9drJyqjHR6K2ajazp8zsLWa2r5k9bmbbzewxM5uTdezpZvaCmbVF/34yxzbqNqQzUJvJCObuuuk25BvwMLAVuATYC5gPGPAYcDdwILAH8C1gC7BrdNz3gZXR9lnAoYQhHIDJgAOfAHYBJvfx2EcALUB9Rtm7gY6Mx/kOsAo4DJgDnAhsBz4QbZ8dPdZq4Jhonxn91S86zoFjsuqzGjg7+nlhVI+PRsfvB/w7UN7Hc+mqx6roee0F/Ap4Lvr33cA+wHLgzozjPgS0A58G5gKfie4flUMbLQXuytjeb5vpNnJvRa+AbiP7FgX+01ll7wG2ATVZ5SuBc6Kffwn8uJ/f2yNQe9mnHHgD+HhG2bXA/dHPY4Bm4J1Zxy0B7ol+7graL2Ttk3P9sgL/w8BmYOwg27GrHqdnlB0ZlX04o2wxsC3j/m+y6xkF+GODaaOM/e8abJvpNnJvGtKROKzIur8/UAusN7NtXTdgX2D3aJ+rgOPN7I/RcNAhuT6ou3cANxN60ZhZFXA0cFO0y95ANXBvVj0+lVGPLsuz7g+3fv8DvAK8bGb/ZWanmNnYQRz3dMbPb0T/PpNVNsbMaqP784mG0TI8Rnjug2mjbLm0mYwwOmkrcdiedT9FCKZ39rLvFgB3X2Zms4DDgUXA3WZ2q7ufmuNj3wQ8YWbTgbcBlcDPM+oBcBTwatZx2Sdluz2HQdTPCUNXmSoyjt9qZm8F3gW8DzgXuMDMDnD3v/XzfDLr5f2UDdRZy5x+118bZculzWSEUeBLPvwBmAqk3f2lvnZy9w3AjcCNZrYM+KmZneHurYRwKRvogdz9STN7ATgBOAi4w923RZv/BLQCs9z9oVyfxAD1Ww/s2rWvmU3NvB8d3wE8BDxkZl8H1hGGaa7JtS79+DNwMPCjjLJ3EJ57Vz36a6Nsw2ozKW0KfMmHBwjDDHeY2TnAXwgnXw8DHnD3X5vZNwlvDM8RXocfBl6KwhTCePgiM3sEaHX3xn4e778IJ3hnR78H2NHLvhi42MwMeBSoA/4P4c2oz+AdRP0eAv7NzB4HOoELCCdHu44/kjAE8iiwiXCidCwhoOP0XcJsphXA/YQ2/igZ7RDptY2yDafNZAQo9kkE3Ub2jXDS9opeyscClwFrgTZgDfAzYPdo+1cIYdpECMR7gPkZxx8F/JXQ0189QB3eRBjCeIOsWTCEYZfPsLPnup4wvv6+aPvs6NiFWccNVL9pwDLCyekXCePiq9l50vYdhNk1GwknQZ8FTu3nOfSoB2GmjwOzM8oOi8rqMsrOAF6I2uoF4JM5ttFSus/S6bfNdBu5N11pKyKSEJqlIyKSEAp8EZGEUOCLiCREwQLfzH4crUvybKEeU0REdipkD38pYYaBiIgUQcHm4bv7o2Y2O5djJk2a5LNn53SIiEiirVixYoO7T+5tW0lfeDV79myWL89e4mT43vyN+wa13zPfeH/PwvZm6GjpXvaT4wf3wCfe3LOssg7KKroVHfSTgwb165448YnBPW7MhtV+LVvAO7uXDaP9msoqafeObmWffvDTg/p1Vyy6otv9MiujrrJucHUZhiG3XzoNrb0sXT/E9uvE2JbKXh1i6O1XkaqgtqK2j71z19aRprmts0f5x6///aCO/9EpB3S7X1meoqZywIu3B62ts43mjuZuZUNtO4Ca8hoqyyqHXS8ze6WvbSUX+GZ2GnAawG677Vbk2vTiH6/CpqzVAlq3dL/fGLV3w6zu5a/+tufv23U/GDctvvqVutatsO5PkM4I6SG1n0HDLGom7wUdzbR0ttDa0UprZytt6TbSnt5xW7t1LQAzxs4gZakdt00tm6guq6aqvIrqsmqqy6vjf75xMoPG1bD1793Lh9J+liI1ZR9aKipY17SOxtZG0p4GYFt791UXMtsv01PrngKgtryWiTUTmVYX7+vYcf70+ha2t3Z/Q9/S0n1JnzWbmgCYOaH7m83vV2/a8XN5mfFPMxtiDfwyK6Olo4XG1ka2tG6hqaMp57YbUzGGcZXjaKhuYGzlYNbWG56CXngVDenc5e77Dmb/hQsXej56+MOyfQOsHaCH8eD54d9FX+9/P0vBnEOgosSDJm4dbbB5TXjzzP60BP23X6ocxk0PYVY5ZsCHak+3s3jZYtKkueGwG6jI+jQ1Im3fAJtehqYNvW/vr/2sDOqnw4Q3QUXNjuK0p9natpXNrZtpbG1kc+vmHW8Al/3hMgA+99bPAVBdXk1DVQPjq8Yzvno8VWVV8T23LOm089o/mnllYxMt7T17+wAX3fcXAM55/149tpWVGTPG1zBr4hgqy/N7yrI93c62tm1sadtCY0sjW9q2cOmKS4HQdilLUV9Vz/iq8dRX1lNXWUd5Kv4+t5mtcPeFvW0ruR5+yRszCcZMhu3rh/+7GuYkL+wByith4u4hdLavD73Wpo39H1M5JrTXuGmQGnwvrSJVseNj8qgIe4heg5OgdRv84xXY/FrPYbJsFTUwfhbUz+gxhAjsCKP6qnp2Yzc60h1sbN7Iq1vDgpmGMa1uGtPGTCvIsNeOeqWMmRNqmT6+hnVbW3ll43a2tnQMeFxVRYrdJtQybXwNFWWFmZtSkaqgobqBhuoGZo2bRXu6natWXkVLZwtz6ucwrW4aFanivgYLFvhm9lPCtwZNMrO1wNfd/Uf9H1WidnkzrH4MOofx1a3V42HiHvHVaSQyg7op4bZtHfxtZe/BNWlueHOwnuPNiVZVB1P3gYl7wqYXo6GcrE/sZZUweV74VJRD+5Wnypk6ZiqTaydTniqnuqyauQ1z461/DlIpY5f6anapr2bd1hae//u2Xnv8ZSlj98l1zGioIdXL+YlCqkhVUF1eTVmqjFnjZg18QAEUcpbOCYV6rLwrr4JdF8DaJ4d2fKoCpi2AlK5726FuCtROhO3rupenKhT2AymvhCnzwyfPtRlDoJV1MPPA8HodopSlqCqrysvQw1BNGVtNQ20ly1c3dhvfL0sZC2c3MLa6tD7JlVLbKXGGasxEqJ85tGMnz+02fiqEYZ3ssAdIt8Pfn4H0AEMWEoZ5xu6y8/6kucMK+y6d6U48+5NDkVWUpZg5ofvf0OSxVSUX9sCOcyGloHTeekaiSXvCltcgl//Qitqhv1GMRq3b4I3noHlT3/tseS1sn7pvCDXpW+b4fAznLFo7w8ynUgt8CLNUM3WmS6+OAC0dLWFp4hL4lKrAH47yKqibCltfH/wx42dqeKLL1r/D638c3Btme3OYHTVh9/AJqR+n3tv9WxL/sukvvZZfd9h1udV3JGjPmPXU3tz3fv3IbKft7dtZs3UNACffc3K34Ylit9/m5u7TM7OnaxZD9mvszxv/TKd3cvKyk7udsC1W2ynwh2vctNwCf2yC5twPpPGV3D4dATS+HE5253D+Y17DvJLsoeZF5knvgWbuDEJFqmLHydpS6KFmas/q4nd0lt7/8Z4Ne5L2NNbj64+LQ4E/XLUTw9zmwfxxVY1L5jTMvkzdG/72FLRlfwd6H1IVYYbUAGFf7J5nUWVemzDE6ZMjpf3qayrYtG3nTLlxNcUfvy/1ttNJ2+FKlUFNw+D21fhzd1VjYdY7wrTBfucnG4zfDea8C8ZOLVj1RqS6rvZJham/o9j08TXd3vuzT+JKT+rhx6FmfN9XPWYa5X+AQ5JKhWmX9TNh/apwBW6m2okwZe8w51wGVjMBSIUTtqN82m91RRnTxoeQL0sZU8bq0/NAFPhxGMQl/jntl0RlFbDLvlA7YWdZw2yYvJdOcucilQqfOkto7nc+7TIuhHxlga6mHenUSnEY7KJbMcyJHvXGTYOy6nBeRGE/RJaYdqutDG9sZUW+qnakSEY3IN8Gu7ZLQnpdw1ZWAUZiQit2ns599tMI1TX3vvTm55Qm9fDjYIMIfEspwHKithoS9zBjLCFXJjdH6+mkS/Siq1KjLmccBvOlBTF8scGodd0Hut9f92fAe5afenfBqjSiZLaTp8MyFdiobb/jr975xT9tnWnWbGoiZdatHODm0wf3RUJJosCPQ3llCPT+Vs8s4JKyI96u+6EP6UNkKZj2T8WuRcFUlqV4ywzNfhssBX5caif2f8Vt5uwT6W6U9DyLJmHtp5770GkMPy5jdx3edhGRPFPgx2XM5L5XJ6yZAJXxfbmziMhQKPDjkkr1vTBa/fTC1kVEpBcK/DhlfvlEF0tlrG8iIlI8Cvw41TT0XASspiGWL6IQERkuBX6czKA2a+XMGs3OEZHSoGmZcci8wKWzDf7xavj54f8I6+x0LamQsOlzIlJaFPhxK6uEXZNz4YuIjBwK/Dio5y4iI4DG8EVEEkKBLyKSEAp8EZGEUOCLiCSEAl9EJCEU+CIiCaHAFxFJCAW+iEhCKPBFRBJCgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgmhwBcRSQgFvohIQijwRUQSQoEvIpIQCnwRkYRQ4IuIJIQCX0QkIRT4IiIJocAXEUkIBb6ISEIo8EVEEkKBLyKSEAp8EZGEUOCLiCSEAl9EJCEU+CIiCaHAFxFJCAW+iEhCKPBFRBJCgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgmhwBcRSQgFvohIQpQXuwIiIkX148Mg3dn3dktBqhxOvbtwdcoTBb6IJFtnG7Q37bzf+Er4t2FW+Le8OgT+KDA6noWIyFAdez288ezO+w+eH/5d9PXw76S5MHH3wtcrDzSGLyKSEAp8EUm2besG2P4GuBemLnmmwBeR5Nr0EmwfIPBbNsP6VYWpT55pDF9EkscdNjwfAn8wGl+GdDtM2QdSI7efrMAXkWRJp+H1pwYeysm2eS20NcH0/aFsZEbnyH2rEhEZivV/yT3suzRvgjeeibc+BTTowDezJWa2bz4rIyKSdwON2Q9kqG8WJSCXHv4BwB/N7EkzO83MxuarUiIieTNmyvCOrxvm8UU06MB394OBvYFfAV8HXjezG8zskHxVTkQkdpPmQkXN0I4tq4DJ8+OtTwHlNIbv7qvc/UvATOAjQB1wv5n91cy+bGYT8lFJEZHYlJXDxD2HduyEN0FFdbz1KaChnrStAMYB9UAZ8CpwMvCqmZ0YU91ERPKjfIihPdTjSkROgW9mC83sB8DrwEXAb4E93X2Ru+8DfBG4NP5qiojEJJ2GjX8d2rEbX4TOjnjrU0C5zNJ5BnicMJyzGJjl7l9x95czdrsVmBxrDUVE4tLZAa+tgObGoR3ftg3WPgkdbfHWq0By6eHfAsxx96Pc/Zfu3mMBaXff4O6a2y8ipSedhrW/h6YNw/s9LZthze+gsz2eehWQ+SAXBTKzSiDl7i1Z5dVA2t1jf8sbO3as77///t3KvvrVr/Le976XlStXctZZZ/U45oILLuDtb387jz/+OOedd16P7UuWLGHBggU88MADfPvb3+6x/eqrr2bevHnceeedXHLJJT2233jjjcycOZObb76Zq666qsf22267jUmTJrF06VKWLl3aY/s999xDbW0tP/jBD7jlllt6bH/44YcBuPjii7nrrru6baupqWHZsmUAfOtb3+LBBx/stn3ixIncfvvtAJx77rk88cQT3bbPmDGDm266CYCzzjqLlStXdts+d+5crrnmGgBOO+00nn/++W7bFyxYwJIlSwA46aSTWLt2bbftBx10EBdeeCEARx99NBs3buy2fdGiRXzta18D4PDDD6e5ubnb9iOPPJKzzz4bgEMPPZRsxx13HGeeeSZNTU0cccQRPbYvXryYxYsXs2HDBo455pge2z/1qU9x/PHHs2bNGk4++eQe27/whS9w1FFHsWrVKk4//fQe2/XaG+GvvbM+Cev+xOH/98s0t2TEVePLHLlvA2e/bxoAh160PJRnzOQ57q0TOfM719LU3MIRHz9v5/byKqC0XnsHH3zwCndf2GMnclta4VbgEeB7WeVnAIcC/5rD7xIRKayxu4Rb7USwjM7G5vIQ3tUN4X5VL5cYVdTC7u+BpiaoGV+Y+uZBLj38DcCh7v5sVvk+wK/cPfarERYuXOjLly+P+9eKiIxaZtZnDz+X8fZaoLfT02lAV92KiJS4XAL/aeCEXspPBJ7tpVxEREpILmP43wTuMLM9gIeiskXAscCH4q6YiIjEK5e1dO4BjgJmAZdHt92Af3H3u/o7VkREii+nVfzd/V7g3jzVRURE8kgXSYmIJEQuSytUmtn5Zva8mbWYWWfmLZ+VFBGR4culh/8t4BTgEsJUzC8CVwIbgTPjr5qIiMQpl8A/DjjD3a8GOoE73P2zhC9DeV8+KiciIvHJJfCnAn+Kft4GdF1ffC/wz3FWSkRE4pdL4L8KTIt+fgF4f/TzQUBzr0eIiEjJyCXw/5twoRXAZcD5ZvYysBS4NuZ6iYhIzAY9D9/dz834+TYzWwMcDDyvC69ERErfoALfzCqAm4Dz3P1FAHf/HfC7PNZNRERiNKghHXdvJ5yYHdxayiIiUnJyGcP/OfDhfFVERETyK5e1dF4Fvmpm7wSWA9szN7p79jdhiYhICckl8BcDjcBbolsmp+dXH4qISAnJZZbOnHxWRERE8kurZYqIJMSge/hmdnl/26N1dUREpETlMob/5qz7FcBeQBnwVGw1EhGRvMhlDP/d2WVmVg38CPh1nJUSEZH4DWsM391bgAuAr8RTHRERyZc4TtpOAupi+D0iIpJHuZy0/Xx2EbAr8FHgnjgrJSIi8cvlpO1nsu6ngfXAdcCFsdVIRETyQhdeiYgkxKDH8M2sMpqVk11ebWaV8VZLRETilstJ21uBM3spPwO4JZ7qiIhIvuQS+AcD9/dS/j/A2+OpjoiI5EsugV8LdPRSngbGxlMdERHJl1wC/2nghF7KTwSejac6IiKSL7lMy/wmcIeZ7QE8FJUtAo4FPhR3xUREJF6D7uG7+z3AUcAs4PLothvwL+5+V36qJyIiccmlh4+73wvcm6e6iIhIHuUyD/8QMzukj/J3xVstERGJWy4nbS8FGnopHxdtExGREpZL4M8D/thL+bPRNhERKWG5BH4zYXXMbNOBtniqIyIi+ZJL4N8H/KeZ7RjWMbMJhJUy74u7YiIiEq9cZumcDTwKrDazp6OytxCWSD4+7oqJiEi8clke+XUz24/whScLouLrgZ+4e1M+KiciIvHJaR4+Yaz+OWAr0LUk8jFmhrvfEGvNREQkVrl8xeFewJ3AHMLXG3ZGx7cDrYACX0SkhOVy0nYJsAKoB5qA+cBCYCVwdPxVExGROOUypHMAcIi7bzezNFDu7n8ws3OA7xNO4IqISInKpYdvhJ49hJk506Of1wJ7xFkpERGJXy49/GeB/YCXgCeBL5lZJ/BJ4IU81E1ERGKUS+B/BxgT/fxV4G7gV8AG4LiY6yUiIjHLZR7+fRk/vwTMj660bXR3z0flREQkPrnOw+/G3TfFVREREcmvXE7aiojICKbAFxFJCAW+iEhCKPBFRBJCgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgmhwBcRSQgFvohIQijwRUQSQoEvIpIQCnwRkYRQ4IuIJIQCX0QkIRT4IiIJocAXEUkIBb6ISEIo8EVEEkKBLyKSEAp8EZGEUOCLiCSEAl9EJCEU+CIiCaHAFxFJCAW+iEhCKPBFRBJCgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgmhwBcRSQgFvohIQijwRUQSQoEvIpIQCnwRkYRQ4IuIJIQCX0QkIRT4IiIJocAXEUkIBb6ISEIo8EVEEkKBLyKSEAp8EZGEUOCLiCSEAl9EJCEU+CIiCaHAFxFJCAW+iEhCKPBFRBJCgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgmhwBcRSQgFvohIQijwRUQSQoEvIpIQCnwRkYRQ4IuIJIQCX0QkIRT4IiIJocAXEUkIBb6ISEIo8EVEEkKBLyKSEAp8EZGEUOCLiCSEAl9EJCEKGvhmdpiZrTKzF8zsy4V8bBGRpCtY4JtZGXAlcDiwN3CCme1dqMcXEUm6QvbwDwRecPeX3L0N+BnwwQI+vohIopUX8LGmA2sy7q8F3pa9k5mdBpwW3d1mZqsKULd8mARsKHYlRjC13/Co/YZnJLffrL42FDLwB8XdrwGuKXY9hsvMlrv7wmLXY6RS+w2P2m94Rmv7FXJI5zVgZsb9GVGZiIgUQCED//fAnmY2x8wqgY8Avyzg44uIJFrBhnTcvcPMPg3cB5QBP3b35wr1+EUw4oelikztNzxqv+EZle1n7l7sOoiISAHoSlsRkYRQ4EvJi875iMgwKfALwMymRVcayyCZ2YFmdoWZ3Q98uNj1SRozq47+Tfzr1szqi12HuCjwY2ZmFv27n5ldaWZ/AL4I7F/cmo0cZnYsYRmOV4B/B17rCiDJPzP7f8D3ANy908zmm9mEIlerIDL+fhdGHY7HgU+a2cwBDh0RdNI2RmZW5e6tZjYLuBp4CPgpsAmY5u5/LWoFS5CZ1QInEd4QPw848AiweJTP4ioZZpYiZEFndP9Q4CJ3PzAaTnsM+IS7P13EahaMmR0AXArcAfwc6ADa3f1vRa1YDEruStuRJgqsEwiLwt0JXA98AHjR3S/K2FVhnyHqSf2CcDHes8DP3X171JNaDbRG+5VFvUxz9U7ywt3TsONcSQPwONBqZgcBBwO/HK1hH/39HgtUAtdGr7GzgR+4+0+KWrk80JDOEFlwB6H3cwhwg7tfH20+HPhjtF9l1/5FqWiJiv6wGoEfuvvH3P0X0abpQDtQFe3XmbG/DEP0mu0xJm9mbzGzG4DngQ9GixsuAy4g/H9cFu03avIi6+/3PYSFHL9rZpMJr72N0X4VXfsXq65xGjX/gYXWS2BlXjX8IPDu6OeOjP2luwuBj5jZzOh8x+WEK7LrgHkZ46lmZlVmNmu0/OEVgwddwzYNGZuOAP4M7BmtZQXh09d84ETgFDPbo+uTwGiQ8ff7/939FOAbhF7+8cBTwO7Rfu0Z+494Cvzh6QqsGVFgPWJm7wF+Q1hGYpy7pzOCa28z26OoNS4h7r4KWADcDrwMXBoF0q+BRcBh0X5OCKUPjpY/vHzrrTdvZnVmdq6Z/Qb4YXRyHMLqinsCnzKzY6Mx7BcJV8X/EDDgMjNbYmazC/Yk8u9CwvkjgOeA3QjPezvh77fr03kq+neemU0tRkXjopO2w2Rmmwgvkp8B/+3uL0XldwMvAd929zei8dCjCUM/o3I8dCjM7Gyg090vzSgbAxwHfBZ4GHgn0Ab8R9YnKckQdSxSXb34qKyc0L5uZu8gDF9cB0wDriVMed0KLCb0eA8jdAT/E9gLONHd3xsNdRwNXOfurYV7VvllZm8QOhPviW7HEHr3XwMe73pdmtkiwvmM77t7Y5GqO2wK/GHqLbCi8nnAJ4C5hB7UNuA2wsmgtoJXtESZ2RTg1+4+z8zK3b0jY9tM4Ejgt+7+VNEqOQJFJyO/Qwjpy4Ebo47HDOBLwNsJwf6Eu5+Zcdz06LjzgRZCD3gJ4TU+aoZ0upjZZwnP9yLCxIHnovJ3A18hfG/HAkKv/0fA9ZlvqCONAn+YMgMrul8HHAA8Gs0uORBodvdnilnPUmZmzwHHu/uzZpYajcESl/5mK5nZbsCZhK8QfZgwHfhRwnTXbe7+ZTM7A5jv7p8zs0OAK91936gHewahg/Ir4Dx3b8r/Myqu6JPLo+4+v5dt5YQOx2p3X1nwyuWBxvCHyd3XAR1m9lkzuxb4LfA+YFy0/UmF/YAuJLQZCvv+RUMz5dknr83sy4RpwS8CdwP/BmyNhhh/RjgnAvDPhHNMEDom88zsvYSpsDcBB7j7WUkIewB3Xw+kzWw/6DZeb+7e4e6/GC1hD+rhx8LMjiZ8/LsEuN3dW4pcpRFFc+z7ZmaVmUOAZnYEYTjmSnd/PaP8OOBqd2+I7l8LPEC4cChNmDn2OcLJ2aOj33Er4ZPAI+6+pSBPqASZ2UnA5Oxh2dFIgR8DBZbEKRoGPIUQzj92959lbNsHOAc4P2OCQIowd/x3wKnuvsLMPg9MAa5w97VmdiEw192Pjn5HxWjquQ5Hkv5+NaQTg6S8WCT/bIB1hKKTiusIc+OvNLMXgXe5ezNwP/CxaNeHCePxXWPTVwDf7fodCvudkvT3q6UVRIrEel9H6Bz6WEcomrV0CmHpjjrCrJG3ufuGaJefE6Zc4u5/MLPHCGPzuPtr6DukE089fJECy7qs/x3AMnffDkwkBHRbtF9ZxjG1hMv/xwOnA/cCP3X3DRkncJcDv+u6uM/dv+dasE8yqIcvUmDRTJtG4O6MpQxg5zpCldF+nRnHNBGGZQAwsw8A/2pmP+w64Rqd3P0YIn1QD1+kOPpbR2ivvtYRsmgxL8IMnI1EnwZEBkOzdESKJFqW4wXgFsJ03pfN7IvAbOAud18W7fchYKa7X56kGSUSPwW+SJFoHSEpNAW+SJFoHSEpNAW+SBFpHSEpJJ20FSkurUXWZLsAAAA+SURBVCMkBaMevkgR6SSsFJICX0QkITSkIyKSEAp8EZGEUOCLiCSEAl9EJCEU+CIiCaHAFxFJCAW+iEhC/C9b4oa0K/fPCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEjCAYAAAAxP7roAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2dd5hkV3mn369Cd1d1nJmepMmSUAKBjEcEkWQLmyDELiYIMCwSxiKsjbTGYAOSAcsIgwlCy4IlgrSAjYXAgEGAyAavRBiCAkIJMdKMNKOZ7ulc1dUVzv5xbs1UVVe6Vafi/d7nqaerbjj31O1bv3vud74gxhgURVGU/ifU6Q4oiqIo7UEFX1EUJSCo4CuKogQEFXxFUZSAoIKvKIoSEFTwFUVRAoIKvtJTiMg7ReQOh+3tFBEjIrtdtVnjeD8QkY8UfI6LyBdEZM7rx87SbVrYFyMiL2r1cZTuQQVf8YWIvNUTipqCpIJSlj8B3lrw+dXA04GnApuBfWW2aQoRuU5EvlZm1Wbgq66Oo3Q/kU53QOkdRORJwEXAbZ3uS69ijDlSsuhE4DfGmNsLlpVu06q+HGzHcZTuQUf4Sl2IyDjwL9gR6Uwd2+/13t7gjfT3estPEJGviMhBEVkSkV+IyPNK9v0TEblNRJIickRE/lNENlY4znYRuUtE/q+IlB3AiMiAiFwhIg+ISEpE7heRN1bYNiwinxSR33nHv1dE3iIioYJtTheR74rIvIgsisitIvIH3rqoiFwlIg97x9onIv9YsO9Rc42I/AC4GHi6d45+ULpNrf7X6q+IvBN4FXCudwwjImd764qewLzv9Z2C836d93/Pr79ORL4mIheLyEMiMiMi14pIvNy5VLoPHeEr9XIN8AVjzPdF5B11bH8mcAj4c+BrQNZbPgJ8A7gUSALnA/8uIo81xtwlIpuAf8OaNL7obf+kcgcQkVOBbwE3AG8ylfOE/F/gaVhx/SWwA9hWYdsQ8BDwEuAw8ATvu08Dn/S2+VfgVm9dBjgdWPbWvRF4AfBSYC+wFTi5wrH+BHg/cIr3fqWB/tfq7/uBU4G1wCu9fVY9QYjIMHAT8FOvjbXAx4FPAS8s2PRpwAHgmV4fPg/cA7ynQt+VLkIFX6mJiPw51vTwinr3McYcFhGA2ULTgTHmVqxY5nm3iJwHvAj4B+A4IIq9uTzgbbNqklZEngjcCHzIGPPuKn1/FFZ8n2OM+aa3+P4q/U4Df1ewaK+IPB54GccEfwfwfmPMXd7n+wq234EVwB95N6AHgZsrHOuIiCSAlUrmlVr9r9VfY8yiiCSBVA0TzsuBYeCVxpgF79gXAd8XkRONMfnvOA+8zhiTBX4jIjcA56CC3xOoSUepioicDFwBvNwTl3LbvM0zbeRf26u0Nywi7xOROz2TwCKwG8jvcyvwHeAOEfmiiLxeRNaXNLPF2+a91cTe4/eAHPD9ml/2WB9fJyJ7ROSw17//VdA/gA8CnxCR74nI20XklIJ11wFnAPeIyP8RkXMLzUENULP/dfS3Hk4FbsuLvcfN3rFPK1h2pyf2eR4GNvg8ltIhVPCVWjwZmAR+LSIZEckAzwDe4H0eBP4ZK3L518NV2ns/8GLgMq+dM7BmhAEAT0z+2HvdBvwZcK+IPK6gjSngx8BLRWSNqy8KICLnA1dihftZXv8+mu+f18d3YkXwy8BZwG0i8mpv3S+AnViTVAhrjvl2k6LfVH8dUGgqK73pG1RHegb9Rym1+DLWRl0o6HuwdvYzsOaII8aY+wpeGW/fNBAuae+pwKeNMV80xtwG7AdOKNzAWG4xxrwLOxfwMNbWnycFPB87efxtEZmo0v9fYa/zP6jz+z4V+Ikx5iPGmF94powTSjcyxtxrjLnKGHMu1tTzmoJ1C8aYLxhjXg+cC/wh1iTWCLX6X09/V1j9fyjlN8DpIjJasOws79i/8d9tpRtRwVeqYoyZNcbcUfgCloAj3udqBRX2AueIyKaCkfg9wAtE5PEicjrwWWAov4OIPElELhWRMz3T0POxk4N3lvQrCZwHzFFF9I0x92AnFj8hIi8UkV0i8jQReWW57b3+PV5EniMijxKR/JNIvn8xz1RzttggqSdiRfdOb/1ficjLRORUETkRaxufx97YfFNH/6v212Mv8BgROVlEJkUkWuZQ/wIkgE973jpPB64G/r3Afq/0OCr4Sit5E3Zkug/rXQLwV1jvnR9hvXV+7L3PMwc8BevZcy/wAeByY8xnSxv3RP95WEGtNtL/H1jPmquAu7Dmj/EK216NFdh/BX6GNc98oGB9FljjtXE38CXgFu97ASwAb8aaqX6BfQp6jjEmUeF49VCt/7X6C9bb5jfYJ7PD2PNbhNe/ZwFjXt+/4n2vVzfRb6XLEK14pSiKEgx0hK8oihIQVPAVRVECggq+oihKQGib4IvIp0TkkDhMbasoiqLUT9smbT03r0WsD/Zj6tlncnLS7Ny5s6X9UhRF6Sd+/vOfTxljSqPTgTbm0jHG/FBEdvrZZ+fOnezZs6c1HVIURelDROSBSuu6zoYvIhd5eUH2HD58uNPdURRF6Ru6TvCNMdcYY3YbY3avX1/2qURRFEVpgK4TfEVRFKU1qOAriqIEhHa6ZX4Om5vjZBHZLyJ/1q5jK4qiKO310nlZu46lKIqirEZNOoqiKAFBBV9RFCUgqOAriqIEBBV8RVGUgKCCryiKEhBU8BVFUQKCCr6iKEpAUMFXFEUJCCr4iqIoAUEFX1EUJSCo4CuKogQEFXxFUZSAoIKvKIoSEFTwFUVRAoIKvqIoSkBQwVcURQkIKviKoigBQQVfURQlIKjgK4qiBAQVfEVRlICggq8oihIQVPAVRVECggq+oihKQFDBVxRFCQgq+IqiKAFBBV9RFCUgqOAriqIEBBV8RVGUgBDpdAc6wflX33L0vTGGnDm2LiSCiH1//Wuf3OaeWS785oVH32dNtmidIITE3qevffa1be1XnsLzB5AzBmMgHJKi5Z06f4qilCeQgl9IMp3jvkMLAGxbG2coGiYWDXe4V8dYXFlk38I+ALaObiUaijIcHe5wr4q58+F5csbw2K0Tne6KoihVEGNM7a06xO7du82ePXtaeoyb75vinV/9NQBvedYpjAxFeNLx61p6zHpZXFlkzyN7+PAvPgzAxY+/mHgkzhM2P6HDPSvmeVf9iJVsjm9c/PRVo3xFUdqLiPzcGLO73LpA2/CX01kSK8Umk8XlDCuZXId6VMxSemnVsmQmSc50R//y5IcM6Wx39UtRlGICLfhHllbKLp9NlF/ebpYyqwXfYEikEx3oTWVy3lOiCr6idDeBFvyZCsI+k0i3uSflqSTs5Ub+nSRvFUx1yZORoijlCbbgL5UX9koj/3ZTSdjLjfw7RS5njo7wl9PZGlsritJJAiv4y+lsRYFaSmU6bp7ImRzLmeWy65LpZJt7U5lkwTlUwVeU7iawbpmzNcw2s4k060cH29Sb1STSCQzlPai6yaSzlMocfb+Y6g7BL4xjqEan4hiU/qXbr73ACv5csrrgzyU7K/iL6cWK65KZJNlclnCo8/EC88vHBH9huTvmPgym6Onod3O/A2DX+K6jywbDnfvfKsHhriN3kTM5Tlt3Wqe7AgRY8Gt54swlO2vHX1ypLPgGw2J6kfHB8Tb2qDzzBSKfSudYTmcZ6nDg2hVPvYJ7Zu45+jkfx/C6x73u6LKdYzvb3S0lAJSO3F/1jVcxvzLPx//440RCnZfbQNrwsznDYoEpohzzyQydDEqbX5lvan07yOUMcyWmsVqmsnbwSOIRJ9soSr8RSMGfS6appeXZnCkyV7STbC5b1aQDMJeaa1NvKjO/nCabKz6R00upDvXGspxZruvcJDPJrrhpKko7CaTg1+t2OdMh98yZ1EzNaNqZ5drbtJqpxdXnp9MurYcSh+rfdqn+bRWlH1DBr8J0h8RrKjlVc5usyTKzPNOG3lTm8MLq0Xwqnas5Id5KfAl+8lBHzXZKcKjkcdduAif4y+ks83UK0mxipe3++NlclsOJw3Vt20k7dGIlU+SSWUi5G0E7SKQTNU1hhaxkV9Sso7SH7tD74An+ofn6xcgYONRm8TqUOLQqB34lppJTpLOdGU1XE/VOCf7hZH03ykL8PBEoSqPk6I60I533E2ozD8/5i1I9MJtky0SsRb1Zzf7F/XVvmzM5Hlp8iJ3jO1vXoQpUE/WlVIbESob4QJ2X17Xn1rfdhTdW71OdT0ZF+yQPc+LEiYi0L61zaQGZSmgBmd4nb8rJ5rLQ+bCZYAn+XCLNok/Pm9lEmsVUhpHB1p+qqeSU7yjahxYfYuvo1rb6+C6nszXdLw/Np9g52UCfTBYO3GrfbzgNIkN17ZbMJH2Zc/LkzTrtimlYTmeZTaYpdBN78IhNkrd9bfzoslBIyOUMIa0v0NPk54jqfWpvNYES/H0zjaUV3nckwambxxz3ZjUPzD/ge590Ls2BxQNsG9vWgh6Vpx6TzaGFFDsn66zMVThyXzwEn36+ff/fPwrrT66rialE7Ynuivsmp9om+PcdWuTNf1z8nd53010AvPlZpxQt3zeTYMe67qpupvgjP8JfyXZHQsbA2PCX01kemS+fjKwWB+aSLS+KMp2cZmFloaF99y3ss4+MbeLQQu3zOJ9Mk1xpoE8rBU84VaKNS2nEfu9iXz/MJlY4OFf/NXj/1JImpOtx8q7TKvhtZv9MomawVSVyOXhotrUZKhsZ3edZya3w8NLDDntTmeV0tmJa6VIausEWivxKfeatZr1tljPLVVNZuMAYw10H/d3Qs1nDbw+3tl9Ka8kLfjLTHRluAyH4uZxh/0xzJ9zeMFrjWzWXmmvaPXD/wv62+JT7EfEDPkazRyka4Sfs3bYG08lp/8cpbWO5+TaqcWBu2ff8EcCB2eWifEVKb5G33XdLhttACP6hhRSZbHNimErnykaWuuDhxeZH56lsquWiBf6edJZSmVW5dqpiDKQKR8GmLrPOkeUj9R+jhW1UIpcz3H+48R/8bw/pKL8XMcYcNbU2aq51TSAEv1HbfavaKSRncnVF1tZDq33KZ5ZWSPjMeb9/1sdEeToBuZJRcKr6k48xhplU8xHH86l5MqXHdsSB+eWmbPHTiysdjV5WGmMxvXhs0ja30hW1qPte8HM54yy/y9RiyrnZZC4158xlq9WpFvLug354ZH6ZVKbO77dcJulZuWUFLKWXnAi1wbQs6vbB6eZ/6PsaOPdKZyn9PboYmDRL37tlLqQyqzI6Nkoma0isZBl26JPvcrIwnUuznFlmqE7fdT8spTINRdDmclasTtwwWnvjZJkfRLllBbgU6fnUPGuH1jprD2xm1kopKPxwaGGZTHaUSLi1Y7RUNsWrv/nqmtuJCJ997mdb2pdKFAauZbyaygJES85NpwPXSk2s08lptoxs6VBvLH0v+IkVt4/pSysZp4K/nHVrJkplUy0R/L3Tjdug980k2bFueNUPchWJMnb01AJk0xCOlt2lkWCrSrhsK4+rNBO5nE36t2HM/f+2kAfmHyCRKX6a2L9go7+3jm4tWp7MJIlF2heFXo7ldJb7PU+mx26d6GhfCkllU6vSdM+mZknn0kRD5a/ldtD3gu/af951e67txumce1vvcjrry3+8lGzWekntqhaIlUlVnqBNTMPopvKrHNpFW+E6V6uymh9mEumWCn46l+bg0kEufvzFRcvzFcNKlz+08BAnrjmxZf2pROHI/T/vOcy7b7wTgE9dcKbTwVgzlEvzkTM5ppPTbBoufy23g7634bsy5+Spw0vQF65zuAjuQ/EfPNJ4DEOefUcS5Kr9LxarTDhXWecyoCWVdZ/0rVZltU61VY4Diwd81Vg4mDjY1oC/UpbTWdIFA7CFDhUsKkelTLYHlw62uSfF9L3gtzMpViOExW1GJdc5dbI54yTobCWT45FqEbrVBH/pEJXuOC6fkFw/beVypml34EJaGe1tjPEdvJfJZTqabbQ0zXm3xCsk0omKbpizqdmWDCzqpe8F33XuqZDjM+banue6vanFFFlHolUxECuzAtVy4WTTsFR+vevCEi69sLKOPbpyLQysm0nNsJzxb7ZrV4R3OWaT3VdPGeyTTzUeWepcHYu+F/ywY8WPOFb8aIXJyIbbcyz40w6DzWYTK+UFdeEA1DIlzJdPGx0Sd/+PkIScPhGGXJvrWviw2qipYWFloWNRpKXu1gvL6bYXLCpHraeeThYuUsH3SbeP8F2bdJJpd2aOXA6W02V+kHP7au+8eMg+CZQwGB500DPLQHjAWVtgrz2X199Ai1wys7lsU8F/nRCw5XR2VaoKYzpfU3kuNVfzSWkpvdSxm2R3TGm3ENejLNftubThux6hQkXTeePtlZpgEkdK0ilU2jFnbwzrTihaHI/Gnfnit8LFcDAa8h2dXImhaIVrpckCMkeWj/iarC3lcOIwx48f3/D+jTC1WN4OfnghxcYWu65Wo94iPIcSh9g1vqvFvVlN34/wux23JgT3/86KItMAIjAUKWlv5nf1NzD7wCo3qZHoiIOeWUYH6ggO88lwvVW/6iA+UOV/kUvbNBSpeXj4l/aVmvduptXv2s2mh05mkm1PG1Cp9OjUYqq6N1iLmVqu70nJVToVv/T9CN/5CNV5e+4abEW2zHUjA0354BcyER8oruC0slTdO6eUTMra+8ePRSuuGVrjpG8AE4PuA3eGByPOgq8qVl278EbY/3PrzQTw3XfZv3/wNvt302OLzlkhxhgnKTmmk9PEo/HaGzoglckyU8F0k8kappZSbBht/yh/Kb1U98R3fttWBElWo+9H+N3uKdHMo3Qr28qzYXSISNjNU8iq2sBHfIzu85Q8EQxHh53Y8cMSbonguyyNWTGoaHnumNiX48hvK45UFtOLToL12pkn5tB8qurA65G5zrg9+s242urcV+Xoe8F3/XjnOpDLZa1Lg3Eu+uGQOCmzFx8Ms3GsQJgzKZh/yH9DqQVYLDZBbIhvaLJ3MBmbbIlJbHjQjUksFKpi0pm6r/rOK0swX959cjY122TPLHOpubbUY4DadRYOLy53xFtndtnfuexEMrX+F/wuH+G7Lm7cisjH7WvjTdvyT9o4WjxfMbuvtitmJWb2Fn3cGN/YeMfybQw330Y5XNnwY9FI+fme5Gz10X2e6fvKjvJd5WnPmuyqHDytILGSWRVwVUou5y6HUb0YY5hbqZ7ZtRRXN1s/BEDw3bbnehDjekSew/3IJhwSTtrU+OTo+tFBJkcKRve5nJ2AbZTEVJFnz8jASFMTroPhQdYMupsLKCQUEgajzf/MKo7up2uM7vOkE2WfqFy6B7a6TCTAI/P1CXkraldUI5FJ+I7UXsmuNBTs1gx9L/juI227O/dNK3LpgLXlrxvx76ceCtnRfRFLh6DZHDhzxYFYxw0f13BTm4c3tzQFhwtPp7I3jeQsLPnwsJkutuUbY5wmjGtH3dZ6R+4ziRUybTTrzNco1FOJ0oyarabvBb9mSl6/7TmP3HXrKOW6vUKsWcbfPtvXDhMrHZ02YrsvZf6hIvFaH1/fUEyDIGwe2dx8f6rgImCq7HXsx6UV7Ch/8ViQVDqXdvqE2eocMSuZXE1zTp5cDo44zFRaC7/mnDytKrpTib53y6zqu9wAq8QLCgJfcuVtPhICpGzgi8tgn8HwYEsmHvMMD0bYODZUt5tmOCzsWFfiqpfLVsyL44ts2hZHiduCJZFQhA3xDRxYOuCrmXWxdU6jdcvhIto2WhrinU7CQgMRrjN7j6aadp1KuxWpuQuZ8SngM0vptrlnNjrCb7fg160OInKliDymlZ1pBaNDUWc5SMIhqT4Jt5KA5Vk4cKt9Lc/aVxU73ejAqDORHhscc9JONbaXCngVjhuPrR6ZJo40PllbSsmNo5E84+3ITe4iOntVE3MPUSugqizJGUhZW7vr+SPXDgil+K3r2646wOlsuuEJ68WVxbammPYzwj8T+EsR+TnwCeBzxpjuKMVehXBImIgPVAzU8MPa4YHyNvwLb7Qj199+zxbhzge+nPMO+3doAnaUL7cWkhBrh9Y6ibybHJpsuo1ajA1FiQ+G60oXsKlciHuNGrW+KGlrfHCcWCRWty05Eoo4L2lYDhf5l1Zdd82YxeYfhvUnNdehMrRq/ijPgs/0x4upNMaYlqdIb9ScA8dqKbsMIKxG3ZeiMeYpwGnA94F3AAdE5NMi8oxWdc4Vm8fdPNZVbWfxkBX7cizPWl/oCrgYZUZCESZjrRd8oNjjpgKRsDAWKzOecOQGWKktPz7562PrW2oCy+NCCItaWJ639vhGWbSZMV3XYgiH3LZXyqLPnES5HCRWWj969ut/v2r/Nrpn+rrajTF3G2P+BtgGvBQYAb4lIveKyN+KSOuHSw2wcWyIaKS5H/ZgNFRd6GqNuOYqr183tI54pLmw9ONGjmv5Dy7P2FDtDJ9jsWj5kVXaoRtaJrVqzsTPTa9dN0gXA8yiNvx45pRjZQlWEs6zgw6E3LZXyEomV1Tdql6WHNe0Lkezgt1OT51GVTAKjAHjQBh4EHgl8KCIvNxR35wRDgnb1jQ3Obpj7XBll8x0svZE5Pz+ivURRYSd4zsb7lskFGHb6LaG9/dLPdGjFVMKOK4qVdre6MBoXUIWklD7HqMdKH64sI1yxd79kpgmEoo4FelWFjRPNjhSb3S/eknn0iymm4s/mF+Zb5sd35fgi8huEfkocAB4H/Bj4FHGmHOMMY8G3gx8yH03m2fb2njDOWEGoyG2VLthzDxAzQm0TOroo3Q51sfWNxw8tGN0h/O8+tWox698VVbMo7Q+Eq6eIKrxwfG2mHMAIg68dI56+hhjTYTNkrRh/cMDzafNyOMyc2kpiw2O1FtdB9jF6Dxncm3z1vHjpXM7cDPWnHMBsMMY83ZjTKEz8A3Aeqc9dEQ0HGLXZGMX967J4cquddl0fQU8oGqyMBHhxIkTffctFomxZbR8JsRWEQ2Hat48hwYqXFquRbZMe/UkQWtForRKuIi0HczfZFcW3TwleTeN8YHx5tvCPjG1Ir10nqUGhXvJUS2CSjRrvz/aTpvs+H6uxM8Du4wx5xlj/sOY1T5YxpgpY0zXBnNtW+M/J0x8ILw6y2MhMw/U/wNMzVf1nR4fHPedCOyEiRPaNlItpNZ5rLjepd1YQhBebTqqxz21leJUSjzafLhLLH8+lx2NBFcSkMs6M2uNDYy1dA6p3oCrUhZT6Zbmx3c1Mu+6ET7wXmC6dKGIDIlI62ZrHBIKCbvW+xvln7BhpLJbV2bFf7Tj9L1VE/LsGt9Vt4CPD463beKxlFoBbfF2CH6FtuKReM1z2E7BHx4MNzVxGx8MH3vCrKc6WF0YSC0wNjDmZPK2lddhLmeY9+mSeWxfWFhujVknZ3JN2+/zuEpiVws/gn8D8IYyy1+HHf33BJvHhspHy5ZhZCjChtEqnjnT9/p/vE4trMoDU0gsEmPzcH2h/u0uK1dIxdzs2NF9pFI6AZdFMiq0JSJVvZ4GwgNtnfOIhEOM1uHZVImJWIEgV3Hv9c3KEiLChlhz6aUFcZKiuhJzyXQlf4e6aFWKhUQ64Sx4LZPLtCWRmh/BfwrwrTLLvw2c5aY7rScUEnbWacvfuW648ug+tWBT/DbC1D2QrXyj2D62veYIdWJwgvFBN/bXRhiPVRawsv73eRxOElZrazja2LpWsXa48VH05Gih4DvMSOndPJrNJbQuts65i2ch000GTU5XqH/bLK6TxbUjvbQfwY8D5VQqB7Tv+dgBm+vwyx+KlhTsKOXw3TTscZJdgSP3V1w9GB5kfaz63PfW0a2NHdsR47HKKSvWxKv8+AcdXipV2qpWbq/ZmIdGqHotVSEcFtYNFxaOcRnHYAVrODrc1OCh3ifSRmk2t/1cMs1KAz78tXCdLG6l2QyydeBH8G8DXlZm+cuBO9x0pz2EQsKWierRt1vWxCqP7pemmw9+mdlbNQjpuJHK6X4Hw4OsG1rX3PGbJBoOMVZhlF91NDsw4s5Tp5rgVxH1dtVeLWR0KMrIkP/J2w2jg8fs99m0uzxEYOegPLaMNObpFYvEWpqeIrmSbdhDJ48xML3kfpTvOneQ33z6jeDnCvx74CsiciLwPW/ZOcCLgRe47lirOW4ixt6pyo9QVdMoHPlt8x0wWSv6G04pu3p8cJyhyFBZu97G+MaW5weph3XDA8wliifTYgPhqvZ9QiEr+g1mFyyiijdONVFvZYBQNY4bj3HPsr/JuePGC/rawqC19bH1DIYHfY9at4xsaem1OOXIHDO1sMLm8c7837sJP7l0vg6cB+wArvJe24HnG2O+1prutY74QKTiCHXNcLSyW+FKAhKrnJUaY35/VY+dSsnQJuOd8cwpZbLMhHY9eXYYcjD3MDBS1iUzTywSq5jDphM2fIBN40O+EqnFB8KsKXxacl4z9lh7IuJ7lB8JRVqebbRZ+/2xdlLOa+72Ui2LPH5z6XzTGPNUY8yw93qqMeYbrepcq1lfwQOnag7thINc7nmy6apRk+tiq8020VCU0Wh3TJmMDUVXBRVN1lMVK+Yg6KlGGyEJMRRZ/X+MhCItz39fiYFIiPUj9SfyO640/sO1n3tJ8rTNI5t9xXRsim9qqUgZY5h15GGTyRqWHKdZGAq7zbXfjifPvi+AUo7zr74FgKwxzCfT7DtiTTvvu+kuwE5IhkS4/rVlUhqvOJ5JTychVhz8cuE3LwS81KmpefYvWDfOD//iw0RDUT5x+ycAuPbZ17rtS53kzx/YbIT58/dPN93NRPzYU1PZ8wervm9DDFUW/Pz5W0ovkc6li85fJBThk7d/EujM+ds8MVRXvVUR+0RQRHgAmzfT0Ug1UnxzjoaibIhv4OBS5RQghVSbZ3LBcjpHJutuVL6wnK6c46kBXMZyCNKWJ8+6v70XXPV27MTtdmwCtaMYY9qTqtEhYREmYgNMbCm+8KuaJJ3bKyu3JwhjA2M8et2j69m8I8QHwpy+xRPfevs2MGzFqxmvhHjticJ4NA4GxtcVmJA6fP7WDQ8QjYRqZn6ciA+sNiuKQDTWXGrkQgZW577ZPLy5LsGfGJxo+eR3wnGmS9epkgfCAwxHh8d0kFkAABxPSURBVJ0Ugh8dGG2LScfPES4Hzgfeg02Q9mZgJzZN8mXOe9ZCKo4866HMj6QpBle316mRe700df7yxNfCQn0jyVVEBqv64Hfz+RMR1o8M8vBsdR/uigF/sQl3gl9mLmV8cJx4JF7TJ7wdlcIyjlMiuHxayDMZm3Qi+LXcsF3hx4b/EuB1xpirgSzwFWPMG7HFUP6oFZ3rSoYncTZMjAy59UvvJeJNuJU2s28XUM88R8XJ72FHwiAhiJV/Slofr36MkIQ6ltKj23Bx4wtJiI3DGx30po5j+dh2I3Cn934RyBtRvwn8sctOdTWRQRhx9KMbb18O+66jGU+jLvFSapSJaoFp2KC/iuk/hte7iWMYnqzo5VRrtLlmcE17zA8O0koX4qKYfCmxSKzpm9/62PqWRioX4ufKeRDIz9LcBzzLe/9kwG2McbezZlfzbUgYJrY3306vMhBv0Dwm7ka5HWIgEqqaz6lqaopwFEYdmFOqDDZGBkaqCtDaCk8GrvGb2bYW9ebQ8sv2seZ+x83u7wc/gv8lbKAVwIeBd4nI74DrsEXNg0N8bVUvkboY37LKSyJwNCJc8bV9cd6qeYvU9CRZs7O5g0fjNW+a1YrI1FNgxgXxgbCTAvB5RgZa81QyNjDWcLTx+vj6tsaF1H0GjDFvLXj/BRHZh02odk8vBl41zZodcKCJogUTO9z1pVcZ3QzT9/nfpw+oll46XkuYhsatWavRmJC1x9f0NhsfHOeRxOraDdFQtG2pKUSE8ViUmaXGUiMXEg4Jow2ktqiXXeO7OLLsr/SkIOwac2At8EFd908RiYrI9SJyQn6ZMeYnxpgPBlLsAUY2rgpcqZvBsbLeOYFjcMRf1K2E3JgzuoBq5opYPaaMdSfU3qYc0RiM1Y6oreRjPuLaS60Ga4fdBMmNx6OVa1I7YHRgtOZkdymbhje1Pa9TXYJvjEljJ2ZbVzqm1wiFGw8g6nEbtFPGfWT9HN1kbdh9QLWyh3WVRIyvbew6Wnci9dhJhqPDZVNTtLJubTnqityug/X1pPxokl1juyqm8yglJCF2jLX/Kd+PhezfgT9pVUd6kkZdKnV0f4zR46Bej48+muQerFjkHQYqFY8pZfIkfwcdGKlrdA9WkAYjq0Wy3amlR4eiNaur1UKkchoVl8Sj8brdKzcPby6b+qPV+DFqPQhcKiJPA/YARdEGxpgPuuxYT9DoaLNNLlg9QThiR/kze6tvNzTuJiVDlzBYoR5DNBKq3/QwNGbnNBYO1Lf95Em+IsVjkdiqbK2dyDS6aXyI+w83Hty0ZrhM1HKL2D66nUeWHsFUMYaEJNRWz5xC/Aj+BcAM8FjvVYgBgif4jeI862GPs2anLQZfzWLowhW2i6g0iq90I6jI5KO8iOUa19TQOIz6C+4plxysE6PSzeOxpgS/aqpzx8SjcSZjkxxOVq6XsSG+oWMJ/Px46fTXL84F2Qa9B3LNex30FdGYFaNKqRYiQ30zWZsnFJKyOXV8C/7AsHXxrVInGfBv/mG1uAvStgChQmJemuiZBlIlR8JSPfttC9g6urWq4G8d6Vy1OodergGk0XJzmdbU2OxpqvmWr9nRgqR1nWeojLhXs+1XZO3xVE33MTTupQTxR+koNBqO+kqf7JKtaxozJW0ej7UkwrYa+XxE5RgdGG27p1MhfrJlXlVtvZdXJ1g0WoGo0SeDfia2xgrT8lzxcgn3bQqKoWiYheXia6ihaNCB4epPSGuPb6B3qwW/U2YIsF42A5GQ79q0Wxq8UTTLxuGN/G7ud6uXx9uTM6cSfmz4p5d8jgKnAGHgl8561EvkGky36rgWZt8wsR0O3l68bGxz37hillLO+6QuH/xyTOwoL/iRIRsz0gClAu+64IcfQiHxypLWb8tfMxx1mv/eD+tj68sKvl9ffdf4seH/QekyERkCPgn8yGWneoZGE0i1IfFUTzK6GR65s3hZnW6EvUS+gEwqkyOxkikqwDM2FD1qgvCVhjq+1o70V0oEcXyrb3NYYQGeudTc0QIyl//48iIvnXanod66JsYD00t1+zxsW9P+YvV54tH4qjmQ4ehwR5+SoMmKV8aYZRG5Apsx85/ddKmHGBgucU71sZ+ymlAYRjYc+xwZqqvQSa8yEBZCgxFO3XysGHtT9uaxLTB1T8myxqtS5aswnbz2ZNu3RiPLHTEUDbNhtL6KYUPRcFt876sxMThR9XMncDHUnASCGUkUXwczqx/batKmbIM9SaHgF77vI5wUkCnHyIZiwR8YaWhw0c0FZLatjdUl+FvXxJAOT/SPDYxV/dwJ/Eza/lXpImAz8KfA1112qmeIr/Nfqi8+CdHO2UK7nsJc9z2e977tDI7aeg15GvDM6XYm4gOMDEVYXK7sMBEKlSkA3wFKvXFc1sBtFD8j/L8s+ZwDDgPXYsseBo9QyLoTlj5GV2OthjNUJTJgPXNM1pbzU/xRmLa72RTeXcq2tXF+8/B8xfUbx4YY8BvP0AIKXTMF6UiUcikaeNUsEztsWoB6RvmxtX056nJOKAxZUzxaVeqjMPvoUOdNCK1g09gQ9z6yULFG7daJzk3WFhIJRQgRIkeOochQx01M4CPwSkQGPK+c0uVDIhLc5DDhSP1RjBtOaW1f+gUJ15XRUSnD0XS7UvC+vwiHhE0V0iWMDEUYj3ePG28+UK0TKSnK4edXdQPwhjLLXwd83k13epTxrTbHfTXGtvjL/R5kRNzUbQ0i+acikb6MTs6zeby8eeS4Css7RX5U32l3zDx+flVPAb5VZvm3gbPcdKdHEYENp1ZZH24ol0lgkRBVUwUolTma66a/z994LEp8sNhNVAQ2jneHsObJj/A7kYOoHH4mbeNAuanxHND56edOcu259u/KAmQzXuZH4LuX27+RQfh/H4YLb+xM/3qRPh6dtpSQJ4J9fP7ygWvL6SypTI59R5IAfPDb93DND+8/ul3L3F+rkA9ay7N3bi85k+Oy/7qsSPQ75frqR/BvA14GvKNk+cuBO5z1qJfJu131Ud72ztC/YtV6gnPuhqJhhqJhHru1e02lp66r8uTfAfwI/t8DXxGRE4HvecvOAV4MvMB1x3oKHbkr3YLxl1ysF+nEyL1eujloDXzY8I0xXwfOA3YAV3mv7cDzA1vIXGkRWiCmYfKCr6dQKYOv1ArGmG9i8+YoSmtRwWqMoym79QQqq/Hjh/8MEXlGheVPd9stJdAYgwpWg+RTdmsZTaUMftwyPwSUm40c89YpitJpdISvVMGP4J8M3Fpm+R3eOkVxiApWQxytl2waL9Cj9C1+BD+JzY5ZyhbAf3VhRamIin3DFIq8Cr5Sgp9J25uA94rI840xMwAishabKfOmVnTu7rvv5uyzzy5adumll/LMZz6TX/3qV1xyySWr9rniiis466yzuPnmm3nb2962av2VV17JGWecwXe+8x3+4R/+YdX6q6++mpNPPpmvfvWrfOADH1i1/jOf+Qzbtm3j+uuv52Mf+9iq9V/4wheYnJzkuuuu47rrrlu1/utf/zrxeJyPfvSjfP7zqzNS/OAHPwDg/e9/P1/7WrHzUywW4xvf+AYAl19+Od/97neL1q9bt44vfvGLALz1rW/llltuKVq/detWPvvZzwJwySWX8Ktf/apo/UknncQ111wDwEUXXcQ99xRnAT3jjDO48sorAXjFK17B/v37i9Y/+clP5j3vsYlTX/jCFzI9PV20/pxzzuGyyy4D4DnPeQ7JpA2YyZc1fN7jNvDXzz4epu7h7I8fhg+tK9r/JW+8nDe84Q0kEgme+9znUsoFF1zABRdcwNTUFC960YtWrX/961/P+eefz759+3jlK1+5av2b3vQmzjvvPO6++25e+9rXrlrfE9deKMt1txzmup8chuufWZSiQq+9MtcewMHbj117GM7+x1u8wLVjMQ0vOXMzb/jMHT1x7VXDj+D/NfBDYK+I3OYteyw2RfL5PtpRlOpsehwM/LjTvehNiiZr9UnJNyZny0SK2AIyfYYYH7P5IhLHFjw5w1v0S+BfjTGJFvSN3bt3mz179rSiaUXpT6Z/C9f/qX3/mu9qOU2/LM/Bp55lo+Zf8+1O96YhROTnxpjd5db5LXG4AvwaWADyiSFeJCIYYz7dRB8VRVE6T7a/vZz8lDg8BfgqsAtr3Mp6+6eBFKCCryidpjCtdIeLjvck+UJGfRrH4MdL50rg58A4kABOBXYDvwJe6L5riqL4JhQu/16pj4w3mdunOYn8mHTOBJ5hjFkSkRwQMcb8QkTeAvxv7ASuoiidJFTwk9YRvn9WluzfPnVp9TPCF+zIHqxnzhbv/X7gRJedUhSlQUL58n6iZSIbITlr/+bKlf7offyM8O8AHgfcD/wU+BsRyQJ/DtzXgr4pilIP+QI8YCNt8wV4CpeDpvGuRXoZVha9DzlYnu+7QvB+BP/dQN7H61LgRuD7wBTwEsf9UhSlEUIR2Py4TveiN1k8WPx54WBwBd8Yc1PB+/uBU71I2xnjx5lfURS36MjdDXPF0bvM74d1J/aVaaypb2KMOaJiryhKz7M0BamF4mWZ1OpRf4/TP7cuRVGURpn+bYXl9/WVT74KvqIowWZpCpJHyq9bWYL5h9rbnxaigq8oSnAxBg7fVX2bqXsh1x+BWCr4iqIEl/mHV9vuS8ksw+zetnSn1ajgK4oSTHI5mL63vm2n74dsuvZ2XY4KvqIowWR+P6STtbeD4oC2HkYFX1GU4GEMHLnf3z4zewvSJ/cmKviKogSPhYP1j+7z5NL2qaCHUcFXFCV4zD7Y3v26BBV8RVGCxUqist99zX2XIDnjtj9tRAVfUZRgsXCguf3nm9y/g6jgK4oSLBYf6ez+HUQFX1GU4JBJwfJck20s21z5PYgKvqIowSEx3V3ttBkVfEVRgoOrCddGJ307jAq+oijBIV+ztlmaNQt1CBV8RVGCQS5XULO2STIp++oxVPAVRQkG6QQYh2mOV5bctdUmVPAVRQkGflMp1Gwv4ba9NqCCryhKMMg6NsGoSUdRFKVLyWXdtufSPNQmVPAVRQkGIp3uQcdRwVcUJRiEIt3dXhtQwVcUJRhEY47bi7ttrw2o4CuKEgwGRt22Nzjitr02oIKvKEowCEdg0JHohwdgYNhNW21EBV9RlOAwvN5RO5Nu2mkzKviKogSHkY1u2hnd7KadNqOCryhKcIhNND/ZGo5CXEf4iqIo3c/41ub2H9sCod6Uzt7staIoSqOMbwVpQvrGt7nrS5tRwVcUJVhEBmF0U2P7xid70h0zjwq+oijBY2JHY/utaXC/LkEFX1GU4BGbgKEJf/tEY+7cOjuECr6iKMHE72h9YkfPJ2BTwVcUJZiMbLIulvUgoea9e7oAFXxFUYJJKARjdYr4qI+bQxejgq8oSnCpd9Tew66YhajgK4oSXAZHYHCs+jaRIYitaU9/WowKvqIowaaWT/7opp6frM2jgq8oSrCplVDNVcK1LkAFX1GUYDM4Ys025QhF/PvrdzEq+IqiKPF15ZfH1vRsorRy9M83URRFaZT42vLL+2SyNo8KvqIoSiWzjQq+oihKnzEwbO31pdRy2ewxVPAVRVFEVhc4j8Zt4fM+QgVfURQFVo/mh/prdA/QX7cvRVEUv1x7rv2bXYF0EmYetJ9v/Otid80Lb2x/3xyjgq8oigIQHrCv436v0z1pGSr4iqIEmz4YudeL2vAVRVECggq+oihKQFDBVxRFCQgq+IqiKAFBBV9RFCUgqOAriqIEBBV8RVGUgKCCryiKEhBU8BVFUQKCCr6iKEpAUMFXFEUJCCr4iqIoAUEFX1EUJSCo4CuKogQEFXxFUZSAoIKvKIoSEFTwFUVRAoIKvqIoSkBQwVcURQkIKviKoigBQQVfURQlILRV8EXk2SJyt4jcJyJ/285jK4qiBJ22Cb6IhIH/AzwHOA14mYic1q7jK4qiBJ12jvCfANxnjLnfGLMC/Bvw39p4fEVRlEATaeOxtgD7Cj7vB55YupGIXARc5H1cFJG729C3VjAJTHW6Ez2Mnr/m0PPXHL18/nZUWtFOwa8LY8w1wDWd7keziMgeY8zuTvejV9Hz1xx6/pqjX89fO006DwHbCj5v9ZYpiqIobaCdgv8z4FEisktEBoCXAv/RxuMriqIEmraZdIwxGRH5C+AmIAx8yhjz63YdvwP0vFmqw+j5aw49f83Rl+dPjDGd7oOiKIrSBjTSVlEUJSCo4CuKogQEFXxFUZSAoILfBkRkq4gcJyLifZZO96nXEJGhTvdBKU+/X88iskFERjvdDxfopG2LEJEJ4C+AlwFJ4HfAPcaYt3e0Yz2EiEwCrwXOxkZm32mM+ScRCRljch3tXAARkXHg5cCvjDG3iEjYGJPtdL9agYisBf4n8EJgFrgHuMkY88WOdqxJui7Sth8QkUcBnwAGgSuBm4HtwHwn+9UreIn23gJcjL1RXod9Gj2zn0WmBzgN+Evg68At+YUisgt4PPAfxph0h/rmDBH5Q+BDQAr438BdwCnAfZ3slwt0hO8AEXk28ERjzLtEJA58ERgDXmWM6fmLpB2IyDOwI/jDIvJ84Crg3cAndTTfWfJPVN6N+GPYXC0vNcbMeOt/jL0B/BVWU3rq/yUiawq+yxbsDe0XwFuNMQc72jnHqA2/QUTkdBH5pIjMA/8CDIpIBPg94CnAB4wx9xXaN/vd1umXgnM4C3wL2OWtuhT4pTHm44XioeevfYglApD/H3hPVj8F1gJ/6G33Z9hEY39vLD0h9t61d42I3A7cICJ/7dnp/wib6PGDhWLvnY+ev/5U8H0iIo/2RjS3AiPA+cB6Y8zbjDEZ4HhgCLjdGxkdfYQy+jgFlD2HLwXixpifisjpwEl4JoPCH5mev/bhiXcGQESeICLneKtuBuaAc0QkClwCXJ4fIXc7JdfeGPBe4CDwTuAdwHogbYy5PX/Dg6Pno+evP7Xh+ycBjAPvNMb8fYX1EWCkdLQjItIPF40Dqp3Dh4BRIFHOXq/n0D3lzqk3Qfs3wIVAFvhP4LvGmDtF5NfA6cCbgN8AX/P2GfBqXXQzq649Efl37MTshcDnASMik8aYqcJz0w/Xno7w/bMX2AMcTZ0qIgMi8noR+SesWD0E/Km3LpzfruDCOb6dHe5C9lL+HL4UOA74L6zJYMRbp6N8x4hIKH9tVjinfw48AzuCPxV4j4jEvHU/xIrmFcDvAy8WkSFjzIqIdLum7MVee2eKyLA3qEhgbwTT2KfzFeA8b/tofseC3+/WtvbYId3+z+k6vH/6l4GNIvI2EbkBeAQ72rkP+AHwPeCVInJymRHqc4FPeZ48gaTCOZzBnsMJ4BvAs4AzC7YHQETWich3RaTvcpW3E2NMLn9tisjTRWR7Xqw9U82fAjcbY67HmjjuMMYkvd1/DDyIHfX/M/BmYI+IXAzESo/VTRRce5uARxtjst6NbzfwI6x3zhL2hkfhE4t3k3wj8DbPZbj3MMboy+cLO2n1U6x//QeAJwHhgvWnAYeB24D/gX38PR74X9ibwV92+jt0+lVyDt/vncOQt24E+5T0G+AMrK01gq2ncBn2B7u909+hV17581qy7ARsRsgjWPG+HXivt24AW3/6AeCzWBfjTwHXAyd527zL+/9txY6C/xbr2XIAeH2nv3ON87EOOzD7F+BqYAH4NXCWt/4iIOetO9v7jjuxXkjfAf57p79Dw9+90x3o1RfwEU+8R7zPIe8V9j4/HfgKdvR/K9ZG+EvsyCncqX5306vgHA57n6Xg/P0RdsSVA/4fdjS5BHwbOLvTfe+FVzmh95bHsfb5L2A9yuLAc71ByssLtrkS64d+KfA+bwDzbWAYeBrwE+yEbb7dGPDkTn9vH9dewjsHTy+z/rXAvdiAv58Bi97fF1U6r73w0knbxrkBOzJ9LvB5U+w+eJox5oci8jPsSPaxwK3GmIc709WuJX8On4P94YWwAo8x5tsi8hPsuTsLyACvMMbsq9BW4PFMMmI8U42xvvNjWE+yOeBLxgZGZbFi/UljJyZPAB6DHfmeKyI/8s7zJV67g8aYlIjch7Xbg30ieAgIFUxmLhtjbqE3+DzwZOAz3m81r4VZY7naMzWego07+KExpucr9GngVYN4AVZfAu4yxlwsIicDf4YdwceAU4wxhzrZx27HO4dfBn5rjHl9p/vTi3gT2mWDnUTkDOyNNIn1lf8W8EZjzJy3fiPwT9j5kp9hR/hPA/7CGPNNERnGmjMewZrc3oS9UVxmjDEiEjPH7Po9Ra1rrx88csqhk7YNYuzM/i3AS0Tkd1gb4BOBi40xa1Xsa+Odw5uBs4I8id0M3mg0ByAivy8inxOR74jIecAbsOJ8OtbH/PHAq7xtQ8AFwMnAucaY5wFvxNqqn+DdSM7EC4LD2rt/A1zlib30qthD7WuvH8Ue1A+/Wb4IPBprX/6UdxEp/vgS1hWz56MYW4VUSRbnedT8LVaUX4KdgMxPuh7E2urBms92Y33Nr/K2uQA7yv2lt82rsS6Jfwj8G3ZAMwF8whjzn4XH7RNBDNy1pyYdRekBKpkYvHQAN2FNLv9ojHmb2EyPlwIXGmPWFGx7PtaN8hxjzC9E5CZgA/CP3ibPAR7Guiy+yxjzQEu/lNJ21KSjKB2knkAlz/f7I57duXS/RWxCsyx2xIox5oj3flhsUro8t2I9T17lfb4M63nzYezk+R7g74wxr1ax709U8BWlA+QFu8D+PlSwTgq3waaIPgfYUrqfN+q/BSv8jyk4xF3YAKn/WbDsAeD7WNMNxpifAq8BzjTGbDPGfMR4+XOU/kQFX1E6gOcyOSkibxeR7wAfE5G3eKslv43396tYd9UzCm4QZ4vI073tD2DNOhcWHOII1m7/7IJjJrGC/zkRWeMtS6ura3BQwVeUFiGWcJnlYRF5K3AH8DysMP8UOMXL7ZIrs++XgAtE5EMi8lusi+UzAIwxS9hJ1rNEZLO3LIsNWBsRkVfmGzHGfNMYc5HpkeyWilt00lZRHOOZXUzhJKuInA382hQXeLkcuLaSB4633zjwYmy+mhOwN4aPAN8oFG1P6G/GetS821s2gjUF/RCY7RPPGqUJVPAVpUWIze1/CbYuagx4mrE5/38G7DfGvKBk+8JUvBHgg1gzzb3YfDYXYd1/rypzrDDwSeBkY8yTW/i1lB5GTTqK4hipXuDlMcCjsKPxsqmfPeHPYO3yT8GWz/wo8DnguSKyrsy+WeASFXulGir4iuKewiIb53t283ya7IexNROSnr2+tPDI0VG+MeZGY8xt2DxCYDNX7gYe560v2tcYM9uyb6T0BSr4iuKevTgs8OKlMgh7ybu+DKRb2nulb1HBVxTHeKLdbIGX3/c+528GeRfN1xhjftS2L6P0FSr4itIavo/1p78MG/B0DtYW/1/YHPOz2OjZM0RkXEQiIrIdm/BsHpu5kgLzjnpXKE2jydMUpQUYY46IyE+xEbDvMMYs5X3rvfcXAH+HrRJ1C9ZOvxs7mftuY8yDneq70r+oW6aitAgvj837saUDv+C5TuYKvHHGKC7wcoNGvSqtRAVfUVqEFnhRug214StKi9ACL0q3oTZ8RWktgSuyoXQvatJRFEUJCGrSURRFCQgq+IqiKAFBBV9RFCUgqOAriqIEBBV8RVGUgKCCryiKEhBU8BVFUQKCCr6iKEpA+P/seGJkhssFnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_conn = 4\n",
    "\n",
    "# Plot rest versus movie\n",
    "# chance level: predicting always movie gives 60% (as compared to predicting rest, 40%)\n",
    "chance_level = 3. / n_run \n",
    "\n",
    "plt.figure()\n",
    "plt.violinplot(perf[1,:,:,0], positions=np.arange(n_conn)-0.25, widths=[0.25]*n_conn)\n",
    "plt.violinplot(perf[1,:,:,1], positions=np.arange(n_conn), widths=[0.25]*n_conn)\n",
    "plt.violinplot(perf[1,:,:,2], positions=np.arange(n_conn)+0.25, widths=[0.25]*n_conn)\n",
    "plt.plot([-1,n_conn], [chance_level]*2, '--k')\n",
    "plt.axis(xmin=-0.6, xmax=n_conn-0.4, ymin=0, ymax=1.01)\n",
    "plt.xticks(range(n_conn), ['EC','FC','FC+mask','PC'], rotation=20)\n",
    "plt.yticks([0,1])\n",
    "plt.ylabel('accuracy', fontsize=14)\n",
    "plt.title('rest versus movie', fontsize=14)\n",
    "\n",
    "#plt.savefig(res_dir + 'perf_RM.png', format='png')\n",
    "\n",
    "# Plot 4 classes (rest, movie1, movie2, movie3)\n",
    "# Chance level: predicting always rest gives 40% (as compared to predicting \n",
    "# a movie session, 20%)\n",
    "chance_level = 2. / n_run \n",
    "\n",
    "plt.figure()\n",
    "plt.violinplot(perf[0,:,:,0], positions=np.arange(n_conn)-0.25, widths=[0.25]*n_conn)\n",
    "plt.violinplot(perf[0,:,:,1], positions=np.arange(n_conn), widths=[0.25]*n_conn)\n",
    "plt.violinplot(perf[0,:,:,2], positions=np.arange(n_conn)+0.25, widths=[0.25]*n_conn)\n",
    "plt.plot([-1,n_conn], [chance_level]*2, '--k')\n",
    "plt.axis(xmin=-0.6,xmax=n_conn-0.4,ymin=0,ymax=1.01)\n",
    "plt.xticks(range(n_conn), ['EC','FC','FC+mask','PC'], rotation=20, fontsize=14)\n",
    "plt.yticks([0,1])\n",
    "plt.ylabel('accuracy', fontsize=14)\n",
    "plt.title('4-task classification', fontsize=14)\n",
    "\n",
    "#plt.savefig(res_dir + 'perf_4tasks.png', format='png')\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Confusion matrices ...\n",
    "\n",
    "... are useful to understand the errors performed by the train classifier on the test set. Each row represents the classification of a test sample, with the diagonal element corresponding to the correct classification and off-diagonal elements to errors. Here we separate the 5 sessions to show that the two rest sessions are not separable by the MRL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use vectorized EC as features\n",
    "vect_features = EC[:,:,mask_EC]\n",
    "                   \n",
    "# Individual labels for each session for classification (separating the 2 rest sessions)\n",
    "# Distinct labels for each run/session\n",
    "run5_labels = np.repeat( np.array([0,1,2,3,4],dtype=np.int).reshape([1,-1]), \n",
    "                         n_sub, \n",
    "                         axis=0 ) \n",
    "\n",
    "# Confusion matrices (1 dimension per condition)\n",
    "conf_mat_MLR = np.zeros([n_rep,n_run,n_run])\n",
    "conf_mat_1NN = np.zeros([n_rep,n_run,n_run])\n",
    "\n",
    "# Perform classification\n",
    "for i_rep in range(n_rep):\n",
    "\n",
    "    # Split run indices in train and test sets \n",
    "    # (choose some subjects for train and rest for test)\n",
    "    train_ind = np.ones([n_sub,n_run],dtype=bool)\n",
    "    while train_ind[:,0].sum() >= 0.8*n_sub:\n",
    "        train_ind[np.random.randint(n_sub),:] = False\n",
    "    test_ind = np.logical_not(train_ind)\n",
    "    \n",
    "    # Train and test classifiers with subject labels\n",
    "    c_MLR.fit( vect_features[train_ind,:], run5_labels[train_ind] )\n",
    "    conf_mat_MLR[i_rep,:,:] = \\\n",
    "            skm.confusion_matrix( y_true=run5_labels[test_ind], \n",
    "                                  y_pred=c_MLR.predict(vect_features[test_ind,:]) )\n",
    "\n",
    "    c_1NN.fit( vect_features[train_ind,:], run5_labels[train_ind] )\n",
    "    conf_mat_1NN[i_rep,:,:] = \\\n",
    "            skm.confusion_matrix( y_true=run5_labels[test_ind], \n",
    "                                  y_pred=c_1NN.predict(vect_features[test_ind,:]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the confusion matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAAEbCAYAAABOR5jyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAepUlEQVR4nO3de7QcZZnv8e9vJ4SEIOIhICqQCMhNhKjhooCiHoYA44jCKCgwETFzwQEcXIiKHJQliHIcZSEO+3BTuYkoCDHiIMoKSACDxAQIcpGA3Ey4DSEJJJDn/FHVpNN0967eqerqrv37ZPVK16Wrnt1hP7yXet9XEYGZWRUNlB2AmVlRnODMrLKc4MysspzgzKyynODMrLKc4MysspzgzKyynOBGCEkXSQpJ5zc5dkZ6bEbdvlMk3dXmejemnwlJKyQ9KOl0SesW9TOk99277r6Nr+3qznudpFMl3SNpuaS/pTEfKsn/3Y8Qo8sOwLrqr8DHJR0TEUsBJI0GjgAeGcb1LgS+DIwBdkm3Ab6U9QKSTgEmRcS0Du/9duCZhn2L02tuCNwMvAE4CbgdWAHsCXwVmA0s7PB+1oec4EaWecCbgY+zOhkdALwIzAI26vB6yyLiyfT9I5I+CfwdHSS4tbAoIp5qcew04K3AthHxaN3++yVdVnxo1itcVB95zgeOrNs+kiTZrdWYPUk7A3sAK9fmOmsrrX4eAlzSkNwAiIgXI+LF7kdmZXCCG3kuBaZIepukTYGpwEXDvNZ0SS9IegmYC2wMfDufMIe0ML137VVLZhNIqqYLuhSH9TBXUUeYiHhW0lUkJbfngBsj4hFJw7ncT4CvARsAXwSejYiftfuApE8B59btGpPs1sF1+/45Ii4Z4t4fAJ6t236ldossgdvI4AQ3Ml0A/BB4ATh5La7zPxHxAICkw4C7JU2LiIvafOYa4La67WOAt5AkyJq/Zbj3Qy3a4BaTJO7tM1zDKs5V1JHpBpJexQnA1XlcMCJWkjTuny5pvTbnLYmIB2ovkp7QNfZFxJK1iGMVcDnwKUmbNR6XNFbS2OFe3/qLE9wIFMkkgDsBb42Il9qcOlbS5IbXNm3Ov5Sks+JzecbbwiaSNm14jUmPfYXksZfbJH1a0tslbS3pcOAOYNMuxGc9wFXUESpjKWkr4M6GfXcAU1pcc4Wks4ETJP1gbUpiGdzdZN8+wG8i4hlJuwMnkFR9JwHPA/cApzK8Z/6sD8kz+ppZVbmKamaV5SqqmfUVSQuBJSSPBr0cEU2bTMAJzsz60wfaDNV7lauoZlZZfd3JMGHChJg4aYuyw6ikJSuK7ADN3+iBUWWHkNlzLz1fdggdeeLeJ56KiI2H+3lNGBusWJXt5CUr7yaZ/KFmMCIG17ie9BDJKJYAzm08Xq+vq6gTJ23B72+7uewwKmnWEzeUHUJHNhrb6UQo5fnFg9eVHUJHTt391IfX6gIrVsFum2Q79zePvdiuTS21Z0Q8JmkT4HpJ90bErGYnuopqZsWTsr0yiIjH0r8XAVcBu7Y61wnOzIolYJSyvYa6lDRe0utq70nmH2w583RfV1HNrE/kN8fLG4Gr0tlvRgOXRkTLOr8TnJkVLHv1cygR8Rdg56znO8GZWbFEaY1hTnBmVrycSnCdcoIzs+KVNM+yE5yZFavWi1oCJzgzK56rqGZWWa6imlklCRhwCc7MqsolODOrJAlGlfMgnBOcmRXPJTgzqyz3oppZZbkEZ2aVVGIvatda/iSFpIvrtkdLWixpRrr9EUnzJM2VNEfSnt2KzcwKpoyvnHWzBLcU2FHSuIhYTrIK+WN1x28AromIkLQTcAWwXRfjM7OilDRUq9t9tzOBA9L3hwKX1Q5ExAuxegWc8SQLSphZv8s6XXkBHRHdTnCXA4dIGgvsBNxWf1DSRyXdC/wSOLLZBSRNT6uwcxYvHnJZRDPrBSVVUbua4CJiHjCJpPQ2s8nxqyJiO+BA4NQW1xiMiCkRMWXjjScUGa6Z5WWElOAArgHOpK562ihdAmxLSc5gZlUwkPGVszIeE7kAeC4i5kvau7ZT0tbAg2knw7uAdYGnS4jPzPI0kgbbR8SjwFlNDh0EHCFpJbAc+ERdp4OZ9bOqJ7iIWL/JvhuBG9P3ZwBndCseM+siD9Uys0oqqIc0Cyc4MyuYUMYSXN5tUk5wZlY4JzgzqyQBozJ2MqzK+d5OcGZWLGUvweXNCc7MCucEZ2YVlb2TIW9OcGZWuJLymxOcmRVLuIpqZlUlGJCXDTSzinIJzswqy21wZlZJQgyUlOHKqRib2YgiKdOrg+uNknRnbVW+VlyCM7NiCQbynw/uWGABsEG7k1yCM7NC1R4TyasEJ2kzktX5zhvqXJfgrKnN19+i7BA6Mjj/4qFP6hF/N3GPskPoug6qnxMkzanbHoyIwYZzvgucALxuqIs5wZlZwTpqX3sqIqa0vJL098CiiLijfk2XVpzgzKxY+c4msgfwD5L2B8YCG0i6OCIOa3ay2+DMrHB5LYsaEV+KiM0iYhJwCPDbVskNXIIzs4IJGBjwUC0zq6giHvStX5WvFSc4MytWxupnEZzgzKxQ8oSXZlZlKmlhVCc4MyucS3BmVlkFjEXNxAnOzAolLxtoZtXlTgYzqzAnODOrLD8HZ2aVJHmolplVWN9UUSWNpWEWkohYlltEZlY5ZVVRM5UbJW0m6SpJ/wMsBZY0vLJcIyRdXLc9WtLi2qIRkj4laZ6k+ZJukbRzxz+NmfWgbNOVF1HKy1qCuxDYFPg88DgQw7jXUmBHSeMiYjmwD/BY3fGHgPdHxLOS9gMGgd2GcR8z6zG9XkXdHdgzIv60lvebSbJYxJXAocBlwF4AEXFL3Xm3Aput5b3MrAeU+aBv1q6NRzo4t53LgUPSdrydgNtanPcZ4Fc53M/MesDAgDK9cr9vxvP+AzgtXa5r2CJiHjCJpPQ2s9k5kj5AkuC+2OL4dElzJM1ZvPiptQnHzLolrznLO5S1inoxyRJdD0t6HlhZfzAiNungntcAZwJ7AxvVH5C0E8lah/tFxNPNPpwuITYI8O4p7xpOW6CZdVXvD9U6Kcd7XgA8FxHz65f9krQF8HPg8Ii4L8f7mVmZen1G34g4N68bRsSjwFlNDp1MUqI7J832L7dbH9HM+kNtZfsyZH7QV9I6wD8CO5A8JnI3cGVEvJzl8xGxfpN9N5IuGhERRwFHZY3HzPpHTyc4SduQ9Gq+EViQ7v48cKqk/SPi/oLiM7MKKGvCy6y9qGcB9wFbRMQuEbELMBF4APheUcGZWQVkHMVQ5kiGvYDdI+KZ2o6IeFrSicDvc4/KzCqjH9rgXgJe04YGjAdW5BeOmVVRr49k+BVwrqR313ZImgL8APhlEYGZWXWUVUXNmuCOAZ4A/iBpuaTlJMOsHgeOzT0qM6sOlTdUK+tzcE8D+0p6B7BduntBRNyVe0RmVil9s7J9RMwH5hcUi5lVVM8lOEnfAr4WEUvT9y1FxAm5R2ZmldGLQ7X2Atape9+KB7ybWWu9uPBzRLyn2Xszs4718mMikk6QNK7J/rGSXD01s5YEjBpQplfesj4mcjrJfHCNxqfHzMxa6P2hWqJ5W9vbgWfzC8fMKkcwkFPySpc7mAWsS5K/royI/9Pq/LYJTtJiksQWwD2S6pPcKOD1JCtumZk1lfNY1JeAD0bEC+kUbjdL+lVE3Nrs5KFKcCel8Z0DfAt4vu7YCmBhRPwuh6DNrMLyWLEKICICeCHdXCd9tXySo22Cq83kK+kh4LcRsbLd+VYdV9x3VdkhdOSoHT9ZdgiZLXh2wdAnVUjSyZA5xU2QNKduezBdh2X19aRRwB3A1sD3I6LV6nxtH/RdLyKWpZs3AeukRcLXqDvPzKyBOmmDe2qopQoi4hVgsqQNgask7dhq2Gi7EtwSSW+KiEUkRcJ2D/SOGipqMxuhCnrQNyKek/Q7YCrQcYLbH6hNcLlfzrGZ2Qgh8muDk7QxsDJNbuOAfYAzWp3fbiTDr5u9NzPrVF6PiQBvAn6YtsMNAFdExIxWJ7drg8u8mHNajTUzayqvKmpEzAPemfX8dlXUJ8k+kN5tcGbWlIBRvTbYnjVnENkS+CbJqvSz033vAT4NfKmY0MysGjrqRc1Vuza4V1fLkvR14PiIuLzulJmS7gb+BfhxcSGaWT9TjkO1OpW1c+O9JA/WNfojsGt+4ZhZFfX6ojN/BT7bZP9RwKP5hWNmVTQgZXrlLetsIscDV0ral9VtcLsD2wP/mHtUZlYZSl9lyLqq1rWSdgCOZvWqWjcAB0XEg0UFZ2ZVIEZnH4uaq8yraqWJ7D8KjMXMKki9uCZDI0ljgB2BTWhou4uImTnHZWYV0nOPidST9H7gMmDTJocDP+hrZm2U1QaXtWJ8NnAjyfxL6wHj6l7rFRKZmVWC6P1e1C2Bj0bEX3KPwMwqTp1MeJmrrHe9DdiqyEDMrJpq0yVleeUt6zW/C5wp6TBJO0vaof6V5QKSQtLFddujJS2WNCPd3k7SbEkvSfpC5z+KmfUklTeSIWsV9er07x+lf9dmGaktJ5ilk2EpsKOkcRGxnGSiusfqjj8DHAMcmDEmM+sTPd2LSjJiIQ8zgQOAK4FDSXpm94JX55RbJOmAnO5lZj2g1slQhqwjGf6c0/0uB05Oq6U7kUy/tFf7j6xJ0nRgOsDmW2yeU1hmVqSyHvTN3K4naVtJg5JulnSTpHMlbdPJzdLZOCeRlN6G9XBwRAxGxJSImLLxxhOGcwkz6yoxSgOZXnnLdEVJ+wHzgXeQDLa/laQENl/S1A7veQ1wJkn11MwqrjYfXC8/B3ca8J2IOHHNwPVN4HTgug7ueQHwXETMl7R3B58zsz6lksYydNLJ8PEm+88HjuvkhhHxKHBW435JmwJzgA2AVZKOA3aIiOc7ub6Z9Z5eH2y/GNgZuL9h/87psSFFxPpN9t1IMgSMiHgS2CxjPGbWJ9SLazI0uBAYlDQRuCXdtwfwZZqUxszM6qmQcQpDy5rgTgGWk6yg9b/SfU+TrCj97fzDMrMqKWssatbn4FaRdCacLmnjdF+mqqmZjWxK/5Qh63xw2wCjI+Ke+sSWjkNdGRGNbXNmZok+WDbwfODdTfZPBs7LLxwzq6JeXzZwZ1avplXv9vSYmVlTyXRJ2f7kLWsnQwCva7J/AzxduZm1JQZ6fMLLm4AvSqsHi6XvTwRuLiIwM6uOAZTplbesJbgTgVnAAkmz0n3vI1lh6325R2VmlSF6fDaRiLiLpENhBsn6DFsC1wKTI2J+ceGZWd/rg8H2RMQjwPG5R2BmFdfjz8GZmQ1XMqNvD49kMDNbG2UluHLuamYjSLb2tyxtcJI2l/Q7SfdIulvSse3OdwnOzAolcp3w8mXg+Ij4o6TXAXdIuj4i7ml2ckclOEnrp+uirpNHpGY2MuRVgouIJyLij+n7JcAC4C2tzs862H488APgMGAVsA3wF0lnA09ExDeyXCdvf35mIR+87Mgybt2xz7/n78sOoSNffPcJZYdQWVttsG3ZIXSXQNnb4CZImlO3PRgRg00vK00C3gnc1upiWauopwPbAu8FflO3/7+BrwOlJDgz6wcdPSbyVERMGfKK0vrAz4Dj2i1rkDXBfQT4eETcJinq9t9D8tCvmVlTIt8JL9Mmsp8Bl0TEz9udmzXBbQwsarJ/fIexmdkIlNc4UyVjvs4HFkTEd4a+bzZ3APvXbddKcUfSfBolMzNg9VjUnOaD2wM4HPigpLnpa/9WJ2ctwX0FmClpu/QzR0t6O7A38P6M1zCzEUmddDK0FRE3Q/biYNbB9rNIEtkmwGPAx4ClwB4Rcfsw4jSzEaTXp0siIu4APpF7BGZWaVKPj0WVtF674xGxLJ9wzKx6illvIYusJbgXWN2x0IynLTezloqofmaRNcHt17C9DskTxEcBX801IjOrlKQXtYerqBHx6ya7Z0i6j2T41o9yjcrMKqS8CS/XNq3OAT6YRyBmVl1lrYs67OmSJI0BjiZ5bMTMrKVe70VdzJqdDAI2BFYARxQQl5lVRLLwc293MpzUsL0KWAzcEhHNxqiamSUKqn5mMWSCkzQaWAnMjIgniw/JzKpGJa2OMORdI+Jl4Gxg3eLDMbMqKquTIWtavR3YOY8bSnolnQHgLknXStow3T9Z0ux0IYl5kjwszKwChBilgUyvvGVtgzsb+L+S3kwyddLS+oOtFnxoYXlETAaQ9EOSnthvAMuAIyLi/tp9JP06Ip7r4Npm1oN6feHnK9K/z0n/rvWoKn0/3KFas4GdACLivtrOiHhc0iKSiTad4Mz6XM92MqS2z/vGkkYBHyKZnbPx2K7AGODBvO9rZt2VLBvYg8/BSboAODYi/pzjPcdJmkuy1NcC4PqGe74J+DHwTxGxqklM04HpAGM2ajvJiZn1hPIeExkqrf4TMC7ne9ba4CaSJPejawckbQD8EvhKRNza7MMRMRgRUyJiyjobuGPXrB+UNeHlUAmusLSbziF3DHC8pNHp0K+rgB9FxJVF3dfMuqs24WWWV96ytMG1mwdurUTEnZLmAYem93kfsJGkaekp0yJiblH3N7Pu6OVOhieHCi4iMveiRsT6Ddsfrtu8OOt1zKxfqDc7GVLT8aMaZrYWBnq4BHetB9Sb2XAlj4n0ZoIrrP3NzEaOXm2DKycqM6sQ9eaElxFRTlRmVhnJhJc9mODMzNaaereKama2lspbVcsJzswK5xKcmVWS2+DMrNpcgjOzanIbnJlVmNvgzKyyyirB+UFeMyucMv4Z8jrSBZIWSbory32d4MysUEqHauU04eVFwNSs93YV1cwKl1cVNSJmSZqU9fy+TnDjx4xht4mblx1GJnf87e6yQ+jIhyd9rOwQrCo6G6o1QdKcuu3BiBgc7q37OsGZWX/ooAT3VERMyeu+TnBmVijhx0TMrLLKe9DXvahmVri8elElXQbMBraV9Kikz7Q73yU4Mytcjr2oh3ZyvhOcmRWqlxedMTNbS3Ing5lVmROcmVWR6M1VtczM8uA2ODOrJLkNzsyqzCU4M6ssJzgzqyxXUc2skmoTXpbBCc7MCucqqplVmBOcmVVUOemthOmSJL0iaa6kuyRdK2nDumPXSXpO0oxux2VmxZGU6ZW3Mlr+lkfE5IjYEXgGOLru2LeBw0uIycwKpYyvfJU94eVs4C21jYi4AVhSXjhmlr+sq6JWKMFJGgV8CLimw89NlzRH0pzlzy0rJjgzy400sqqo4yTNBZ4E3ghc38mHI2IwIqZExJRxG65XSIBmVg2ltcEBE0kq3UcPcb6Z9bkRV0WNiGXAMcDxkvy4ilmFjbgEBxARdwLzgEMBJN0E/BT4ULpizr5lxmdm+SirDa7rJaeIWL9h+8N17/fqdjxmVl2uGppZwcpb+NkJzsy6wAnOzCqomDEK2TjBmVnhPOGlmVWW2+DMrMKc4MyskspbNrDs2UTMzArjEpyZFSrpRXUV1cwqywnOzCpqwI+JmFk1lfeorxOcmRVuxKyqZWYjUX6LzkiaKunPkh6QdGK7c53gzKxYOa7JkK7l8n1gP2AH4FBJO7Q63wnOzApVe0wkpxl9dwUeiIi/RMQK4HLgIy3vHRH5/BQlkLQYeLiAS08AnirgukXop1ihv+Ltp1ihuHgnRsTGw/2wpOtIYstiLPBi3fZgRAzWXetgYGpEHJVuHw7sFhGfa3axvu5kWJsvvR1JcyJiShHXzls/xQr9FW8/xQq9G29ETC3r3q6imlk/eQzYvG57s3RfU05wZtZP/gC8TdJbJY0BDqHN4vF9XUUt0ODQp/SMfooV+ivefooV+i/ejkXEy5I+B/waGAVcEBF3tzq/rzsZzMzacRXVzCrLCc7MKmtEJzhJr0iaK+kuSddK2rDu2HWSnpM0o8wY67WKV9JkSbMl3S1pnqRPlBxnSLq4bnu0pMW171LSdmm8L0n6QnmRvhrfUPF+Kv1e50u6RdLOPRzrR9JY50qaI2nPsmLtBSM6wQHLI2JyROwIPAMcXXfs28Dh5YTVUqt4lwFHRMTbganAd+uTdQmWAjtKGpdu78OaXfnPAMcAZ3Y7sBaGivch4P0R8Q7gVMptzB8q1huAnSNiMnAkcF6X4+spIz3B1ZsNvKW2ERE3AEvKC2dIr8YbEfdFxP3p+8eBRUAhD0F3YCZwQPr+UOCy2oGIWBQRfwBWlhFYC+3ivSUink03byV59qpM7WJ9IVb3HI4HRnQvohMcrw7g/RBtnqfpJe3ilbQrMAZ4sNtxNbgcOETSWGAn4LaS4xlK1ng/A/yqa1E11zZWSR+VdC/wS5JS3Ig10hPcOElzgSeBNwLXlxzPUNrGK+lNwI+BT0fEqhLie1VEzAMmkZQwZpYZSxZZ4pX0AZIE98XuRfZaQ8UaEVdFxHbAgSRV6hFrpCe45WlbxUSSSQ+OHuL8srWMV9IGJP/H/kpE3FpSfI2uIWlnu2yoE3tEy3gl7UTSnvWRiHi624E1MeR3GxGzgC0lZR3oXjkjPcEBEBHLSBq9j5fU86M7GuNNh6xcBfwoIq4sN7o1XAB8LSLmlx1IRk3jlbQF8HPg8Ii4r5TIXqtVrFsrnVhN0ruAdYFeSMil6Plf5m6JiDslzSMp9v9Y0k3AdsD6kh4FPhMRvy41yDoN8QbwPmAjSdPSU6ZFxNyy4gOIiEeBsxr3S9oUmANsAKySdBywQ0Q83+UQ19AqXuBkYCPgnDR3vFz2rB1tYj0IOELSSmA58Im6TocRx0O1zKyyXEU1s8pygjOzynKCM7PKcoIzs8pygjOzynKCq6B0tpFT6rYXljFrh6Qp6ewXk7p976EM5zuSNE3SC4UHZ7lxgusCSRelv+ghaaWkv0g6U9L4LoWwC3BOlhNH8C/xGt9R+m91cMM5PwG27GpUtlb8oG/3/IZk+qV1gL1Ihv2MB/612cmS1omIXGbbiIjFeVynyrJ8RxGxnOThWesTLsF1z0sR8WRE/DUiLgUuIRkMjaS90xLD/pJul7QC2Dc99mFJd0h6UdJDkr6RDs0iPb6JpF9IWi7pYUmvmT2isfol6fWSfiDpifS6CyR9QtLewIXA+LoS5ynpZ8ZIOkPSo5KWSfqDpH0b7jNV0r3pNW8CthnqS0mve1oa+0tp6faYuuPvk3Rbes2/SfrPhp//RknnpNd4StKitHQ8UHdOR9+RpIXp7p+m38HCdP9rSreS/lnSA5JWpH9/tuF4SJou6aeSlqY/32FDfS+WD5fgyrOcpDRX7wzgeOABYEmaQC4BjgVmAVsA/0UyvrCWsC4iGXz/v0kmvvxPkpkmmkrHKc4E3gB8GrgP2JZkRfFbgOOA04Ct0o/UfqEvTPd9EngU2B+4VtIuEfEnSZsDVwP/D/g+yTQ+38nwPfyQpER7LHBn+rNsnsb6FpKpiX4MTEvvfx6wKv2eaj4FfA94LzAZuBS4g9UD0Tv6jkiqq4uAzwIzgFeanSTpo8DZwOeB/yb5n9I5kp6MiGvrTj0ZOBH4EslsJBdImhURj7SJwfIQEX4V/CL5BZtRt70r8BTwk3R7b5LxpAc1fG4W8NWGfQeSJB2RlJAC2KPu+ESSX8hT6vYtBL6Qvt+HJEFs3yLWacALDfu2Sj+zRcP+q4Fz0venkSRL1R0/KY1vUot7vS09PrXF8W8A9wMDDfG9BKyXbt8IzG743PXAeen7jr+jdDuAg9t9N8DvSZata/y3vrnhOqfXbY8mSbKHlf3f5Uh4uQTXPVPT6s1okpLbL4B/bzhnTsP2u4FdJdXPPzYAjAM2BbYnSTy31w5GxMOSHm8TxzuBJyJiQQexv4skod6TDjavWRf4bfp+e+DWSH+LU7OHuO47SeL/XYvjtWvWz213M8mEnlsD89J98xo+9ziwSd01Ov2OstqeZFaPejcD/9Cw79X4IlnXc3FdfFYgJ7jumQVMJ5mm+/Fo3oGwtGF7APga8NMm59Y3ihc9Y8JAeo9deO0042U1utf/zI0xBa9tX+7mrBKN98oSnxXAX3L3LIuIByLi4RbJrZk/Atuln2t8vQzcS/JvuGvtA0rmLntzm2veCbxJ0vYtjq8gWTG88TMCNm0SR23BkwXAblqziLf7ED/f3DT+D7Q4vgDYvb7DANgzjTHrlOzD+Y4gSUqN30Oz+PZo2LcncE/G2KxgTnC97evAJyV9XdKOSpbbO1jStwAi4s/AdcC5kt4jaTJJG1C7UtUNJHP4/0zSvpLeKmkfSQemxxcCY9N9EyStF8kkj5cAF6X331LJQ7xfkPSx9HP/RdJw/11J2yp5huxf2v1w6XWvAM6TdFAay16SaquZnUOSiM6RtL2kA4BvAmdHMunnkIb5HdW+hw9J2lTSG1qc823gcElHS3qbpH8n6fD4VpbYrHhOcD0skgk2DyAp4dyevk4E6nvfppEsa/db4FqSHsSFba65CtiPpIH8YpJSyPdI2rWIiFtIktVlJNXgE9KPfpqkJ/VbJKWiGSSTbD6cfu4R4GMkyxb+iaRn8cQMP+YRacxnpde9CHh9es3H0ljfSVLauyCN68sZrltvGh18R6njSb73v5KUYF8jIq4maUf9PEmp7Vjg32LNHlQrkSe8NLPKcgnOzCrLCc7MKssJzswqywnOzCrLCc7MKssJzswqywnOzCrLCc7MKuv/A9OzZ5E2/IePAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAAEbCAYAAABOR5jyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAehklEQVR4nO3de7gcVZnv8e9v50JCEHEICCOQDCo3EaJEvHAR5XC4jcPNURmFiYo5nhMFFQbjoB7ERxDlOMog6h5EFBAEFIZLwEGUCUgAA8QECCJyk0smAeRAQiAJeeePqoZO05fqpKqru/bvk6efdFWtrnp3h/2yaq1aaykiMDOroqGyAzAzK4oTnJlVlhOcmVWWE5yZVZYTnJlVlhOcmVWWE5yZVZYT3AghaU9Jl0t6VFJImtakzPXpsSMa9k+TtLRue6+03D2SRjeUfVDScYX9IMk1osXrk3VlJOkoSXMkPSvpGUm3Szpe0oZFxmf9wwlu5NgAuBM4BljeptzzwFclrZfhnJOAj69LUGmyfHAtPvoJYPOG14/rjp8L/CswC9gb2An4EvAe4NB1CNkGyOjORawKImIWyS87ks5pU/RnwP7ADOBbHU57OnCipPMiYlkecXbh6YhY1OyApA8AHwYOi4hf1B16ELhK0kY9iM/6gGtw1mgp8BXghAyJ4F+BlcDnCo+qOx8G7m1Ibi+JiKd7HI+VxAnOmhkGngRmdij3PMlt3z9J2qTwqNZ0rqSlDa83p8feCNzT43isDznB2StExCrgBOBoSa/rUPxcklu/L2U5t6St6pMScDWwVUOi+n6GU/0TMKXh9YfaZbLEYtXnNjhrKiIuTntDTwJuaFNutaSZwGWSvpPh1I+RJKOatwOnAnvV7Xsmw3kWRcR9LY7dC2yf4RxWcU5w1s7xwHXAU+0KRcQsSb8FvtbphGnt8KXEJGkLYFWbZLU2fgpcKOnQZu1wkjZyO9zI4AQ3QkjaAHhDujlEcls4BXgqIh5u9pmI+E9J1wCfAl7scInjgZtJOh16YSNJmzXsWxoRS4GLgEOA8yWdDFwDLCap1X0auBg4p0dxWoncBjdyTAXuSF/jSXpK7yC5BW1nJjC208kj4nfAJUCW5+fy8G/A4w2vmWksARxO8szf3wK/ARYApwD/Cfy8RzFayeQZfc2sqlyDM7PKchucmQ2UdGjfsyTtwqsiYmqrsk5wZjaI3hMRT3Qq5FtUM6usge5kmDhxYkyavFXZYVTSihdXlB1CV0YPDc7NyJ2P/6nsELqy6vGlT0TEWg/F08RxwYrV2Qo/u/IukiGANcMRMbzG+aQHgL8AAfyg8Xi9wfmvoolJk7fit7fcWHYYlfTYsqaPxvWtjcdtWnYImW17ykFlh9CV/zrxxofW6QQrVsPbM/77/OrR59u1qaV2j4hHJW0KXCvpnoiY3aygb1HNrHhStlcGEfFo+vdi4FJg11ZlneDMrFgCRinbq9OppAmSXlV7D/xPkolcmxroW1QzGxD5ze/yWuBSJbW90cBPI+KaVoWd4MysYNlvPzuJiPuBnbOWd4Izs2KJ0hrDnODMrHg51eC65QRnZsUraY5lJzgzK1atF7UETnBmVjzfoppZZfkW1cwqScCQa3BmVlWuwZlZJUkwqpwH4ZzgzKx4rsGZWWW5F9XMKss1ODOrpBJ7UXvW8icpJJ1Xtz1a0hJJV6bbB0maL2mepLmSdu9VbGZWMGV85ayXNbhlwI6SxkfEcmAf4NG649cBl0dESNoJuAjYrofxmVlRShqq1eu+21nAgen7w4ELagciYmm8vALOBJIFJcxs0GWdrryAjoheJ7gLgQ9JGgfsBNxSf1DSIZLuAa4CPtbsBJKmp7ewc5cs6bgsopn1g5JuUXua4CJiPjCZpPY2q8nxSyNiO+Bg4KstzjEcEVMjYuomm0wsMlwzy8sIqcEBXA6cRt3taaN0CbCtJTmDmVXBUMZXzsp4TORs4OmIWCBpr9pOSW8A/pR2MrwVWA94soT4zCxPI2mwfUQ8Apze5NBhwJGSVgLLgQ/WdTqY2SCreoKLiA2a7LseuD59fypwaq/iMbMe8lAtM6ukgnpIs3CCM7OCCWWsweXdJuUEZ2aFc4Izs0oSMCpjJ8PqnK/tBGdmxVL2GlzenODMrHBOcGZWUdk7GfLmBGdmhSspvznBmVmxhG9RzayqBEPysoFmVlGuwZlZZbkNzswqSYihkjJcOTfGZjaiSMr06uJ8oyTdUVuVrxXX4MysWIKh/OeDOwZYCGzYrpBrcGZWqNpjInnV4CRtQbI631mdyroGZ02NHTW27BC6snL1irJDyGxoaOTVK7q4/ZwoaW7d9nBEDDeU+TZwPPCqTidzgjOzgnXVvvZERExteSbpb4HFEXFb/ZourTjBmVmx8p1NZDfg7yQdAIwDNpR0XkR8pFnhkVdXNrOey2tZ1Ij4QkRsERGTgQ8Bv26V3MA1ODMrmCiv3dEJzswKV8SDvvWr8rXiBGdmxcp4+1kEJzgzK5Q84aWZVZlKWhjVCc7MCucanJlVVgFjUTNxgjOzQsnLBppZdbmTwcwqzAnOzCrLz8GZWSVJHqplZhU2MLeoksbRMAtJRDyXW0RmVjll3aJmqjdK2kLSpZL+P7AMeLbhleUcIem8uu3RkpbUFo2Q9GFJ8yUtkHSTpJ27/mnMrA9lm668iFpe1hrcj4DNgM8CjwGxFtdaBuwoaXxELAf2AR6tO/4A8O6I+Iuk/YFh4O1rcR0z6zP9fov6DmD3iPj9Ol5vFsliEZcAhwMXAHsARMRNdeVuBrZYx2uZWR8o80HfrF0bD3dRtp0LgQ+l7Xg7Abe0KPdx4OocrmdmfWBoSJleuV83Y7nPASeny3WttYiYD0wmqb3NalZG0ntIEtznWxyfLmmupLlLljyxLuGYWa/kNWd5l7Leop5HskTXQ5KeAVbWH4yITbu45uXAacBewMb1ByTtRLLW4f4R8WSzD6dLiA0D7DL1rWvTFmhmPdX/Q7W+mOM1zwaejogF9ct+SdoK+AVwRETcm+P1zKxM/T6jb0T8IK8LRsQjwOlNDn2ZpEZ3ZprtV7VbH9HMBkNtZfsyZH7QV9IY4O+BHUgeE7kLuCQiVmX5fERs0GTf9aSLRkTEUcBRWeMxs8HR1wlO0jYkvZqvBRamuz8LfFXSARHxx4LiM7MKKGvCy6y9qKcD9wJbRcTbIuJtwCTgPuA7RQVnZhWQcRRDmSMZ9gDeERFP1XZExJOSZgK/zT0qM6uMQWiDewF4RRsaMAFYkV84ZlZF/T6S4WrgB5J2qe2QNBX4HnBVEYGZWXWUdYuaNcEdDTwO/E7ScknLSYZZPQYck3tUZlYdKm+oVtbn4J4E9pX0ZmC7dPfCiLgz94jMrFIGZmX7iFgALCgoFjOrqL5LcJK+AXwlIpal71uKiONzj8zMKqMfh2rtAYype9+KB7ybWWv9uPBzRLyz2Xszs67182Miko6XNL7J/nGSfHtqZi0JGDWkTK+8ZX1M5BSS+eAaTUiPmZm10P9DtUTztrY3AX/JLxwzqxzBUE7JK13uYDawHkn+uiQi/m+r8m0TnKQlJIktgLsl1Se5UcCrSVbcMjNrKuexqC8A742IpekUbjdKujoibm5WuFMN7otpfGcC3wCeqTu2AngwIn6TQ9BmVmF5rFgFEBEBLE03x6Svlk9ytE1wtZl8JT0A/DoiVrYrb9Xx6rF/VXYIXXnhxefLDiGzxx8fWYslJZ0MmVPcRElz67aH03VYXj6fNAq4DXgD8N2IaLU6X9sHfdePiOfSzRuAMWmV8BXqypmZNVA3bXBPdFqqICJeBKZI2gi4VNKOrYaNtqvBPStp84hYTFIlbPdA76hOUZvZCFXQg74R8bSk3wD7AV0nuAOA2gSX++ccm5mNECK/NjhJmwAr0+Q2HtgHOLVV+XYjGX7Z7L2ZWbfyekwE2Bz4cdoONwRcFBFXtircrg0u82LO6W2smVlTed2iRsR84C1Zy7e7RV1E9oH0boMzs6YEjOq3wfasOYPI1sDXSValn5PueyfwUeALxYRmZtXQVS9qrtq1wb20Wpakk4BjI+LCuiKzJN0FfBI4t7gQzWyQKcehWt3K2rnxLpIH6xrdDuyaXzhmVkX9vujMn4FPNNl/FPBIfuGYWRUNSZleecs6m8ixwCWS9uXlNrh3ANsDf597VGZWGUpfZci6qtYVknYAZvDyqlrXAYdFxJ+KCs7MqkCMzj4WNVeZV9VKE9nnCozFzCpI/bgmQyNJY4EdgU1paLuLiFk5x2VmFdJ3j4nUk/Ru4AJgsyaHAz/oa2ZtlNUGl/XG+AzgepL5l9YHxte91i8kMjOrBNH/vahbA4dExP25R2BmFaduJrzMVdar3gK8vshAzKyaatMlZXnlLes5vw2cJukjknaWtEP9K8sJJIWk8+q2R0taIunKdHs7SXMkvSDpuO5/FDPrSypvJEPWW9TL0r9/kv5dm2Wktpxglk6GZcCOksZHxHKSieoerTv+FHA0cHDGmMxsQPR1LyrJiIU8zAIOBC4BDifpmd0DXppTbrGkA3O6lpn1gVonQxmyjmT4Q07XuxD4cnpbuhPJ9Et7tP/ImiRNB6YDbLnVljmFZWZFKutB38ztepK2lTQs6UZJN0j6gaRturlYOhvnZJLa21o9HBwRwxExNSKmbrLJxLU5hZn1lBiloUyvvGU6o6T9gQXAm0kG299MUgNbIGm/Lq95OXAaye2pmVVcbT64fn4O7mTgWxExc83A9XXgFOCaLq55NvB0RCyQtFcXnzOzAaWSxjJ008nwgSb7fwh8ppsLRsQjwOmN+yVtBswFNgRWS/oMsENEPNPN+c2s//T7YPslwM7AHxv275we6ygiNmiy73qSIWBExCJgi4zxmNmAUD+uydDgR8CwpEnATem+3YB/pkltzMysngoZp9BZ1gR3IrCcZAWtv0r3PUmyovQ38w/LzKqkrLGoWZ+DW03SmXCKpE3SfZluTc1sZFP6pwxZ54PbBhgdEXfXJ7Z0HOrKiGhsmzMzSwzAsoE/BHZpsn8KcFZ+4ZhZFfX7soE78/JqWvVuTY+ZmTWVTJeU7U/esnYyBPCqJvs3xNOVm1lbYqjPJ7y8Afi89PJgsfT9TODGIgIzs+oYQpleectag5sJzAYWSpqd7tuTZIWtPXOPyswqQ/T5bCIRcSdJh8KVJOszbA1cAUyJiAXFhWdmA28ABtsTEQ8Dx+YegZlVXJ8/B2dmtraSGX37eCSDmdm6KCvBlXNVMxtBsrW/ZWmDk7SlpN9IulvSXZKOaVfeNTgzK5TIdcLLVcCxEXG7pFcBt0m6NiLubla4qxqcpA3SdVHH5BGpmY0MedXgIuLxiLg9ff8ssBB4XavyWQfbTwC+B3wEWA1sA9wv6Qzg8Yj4Wpbz5O3ORffxpm8dVMalu3bjjB+WHUJXXrPeYC3oM2ZobNkhZPbnb15Rdghd2fJ7r1+3EwiUvQ1uoqS5ddvDETHc9LTSZOAtwC2tTpb1FvUUYFvgXcCv6vb/B3ASUEqCM7NB0NVjIk9ExNSOZ5Q2AH4OfKbdsgZZE9xBwAci4hZJUbf/bpKHfs3MmhL5TniZNpH9HDg/In7RrmzWBLcJsLjJ/gldxmZmI1Be40yVjPn6IbAwIr7V+brZ3AYcULddq8V9jObTKJmZAS+PRc1pPrjdgCOA90qal74OaFU4aw3uBGCWpO3Sz8yQ9CZgL+DdGc9hZiOSuulkaCsiboTs1cGsg+1nkySyTYFHgUOBZcBuEXHrWsRpZiNIv0+XRETcBnww9wjMrNKkPh+LKmn9dscj4rl8wjGz6ilmvYUsstbglvJyx0IznrbczFoq4vYzi6wJbv+G7TEkTxAfBXwp14jMrFKSXtQ+vkWNiF822X2lpHtJhm/9JNeozKxCypvwcl3T6lzgvXkEYmbVVda6qGs9XZKkscAMksdGzMxa6vde1CWs2ckgYCNgBXBkAXGZWUUkCz/3dyfDFxu2VwNLgJsiotkYVTOzREG3n1l0THCSRgMrgVkRsaj4kMysalTS6ggdrxoRq4AzgPWKD8fMqqisToasafVWYOc8LijpxXQGgDslXSFpo3T/FElz0oUk5kvysDCzChBilIYyvfKWtQ3uDOD/SfprkqmTltUfbLXgQwvLI2IKgKQfk/TEfg14DjgyIv5Yu46kX0bE012c28z6UL8v/HxR+veZ6d+1HlWl79d2qNYcYCeAiLi3tjMiHpO0mGSiTSc4swHXt50Mqe3zvrCkUcDeJLNzNh7bFRgL/Cnv65pZbyXLBvbhc3CSzgaOiYg/5HjN8ZLmkSz1tRC4tuGamwPnAv8YEaubxDQdmA4weqNxOYZlZsUo7zGRTmn1H4HxOV+z1gY3iSS5z6gdkLQhcBVwQkTc3OzDETEcEVMjYurQBC/PajYIyprwslOCKyztpnPIHQ0cK2l0OvTrUuAnEXFJUdc1s96qTXiZ5ZW3LG1w7eaBWycRcYek+cDh6XX2BDaWNC0tMi0i5hV1fTPrjX7uZFjUKbiIyNyLGhEbNGy/r27zvKznMbNBof7sZEhNx49qmNk6GOrjGtwVHlBvZmsreUykPxNcYe1vZjZy9GsbXDlRmVmFqD8nvIyIcqIys8pIJrzswwRnZrbO1L+3qGZm66i8VbWc4MyscK7BmVkluQ3OzKrNNTgzqya3wZlZhbkNzswqq6wanB/kNbPCKeOfjueRzpa0WNKdWa7rBGdmhVI6VCunCS/PAfbLem3foppZ4fK6RY2I2ZImZy0/0Alu7Hpj2HLSa8sOI5P7n7mv7BC6sssmE8sOobJGa4StJdLdUK2JkubWbQ9HxPDaXnqgE5yZDYYuanBPRMTUvK7rBGdmhRJ+TMTMKqu8B33di2pmhcurF1XSBcAcYFtJj0j6eLvyrsGZWeFy7EU9vJvyTnBmVqh+XnTGzGwdyZ0MZlZlTnBmVkWiP1fVMjPLg9vgzKyS5DY4M6sy1+DMrLKc4MyssnyLamaVVJvwsgxOcGZWON+imlmFOcGZWUWVk95KmC5J0ouS5km6U9IVkjaqO3aNpKclXdnruMysOJIyvfJWRsvf8oiYEhE7Ak8BM+qOfRM4ooSYzKxQyvjKV9kTXs4BXlfbiIjrgGfLC8fM8pd1VdT8E1xpbXCSRgF7Az/s8nPTgekA601cv4DIzCxP6m5VrVyVUYMbL2kesAh4LXBtNx+OiOGImBoRU8dsOK6QAM2sGkprgwMmkdx0z+hQ3swGXFm3qKW1wUXEc8DRwLGS/LiKWYWNuAQHEBF3APOBwwEk3QBcDOydrpizb5nxmVk+ynpMpOc1p4jYoGH7fXXv9+h1PGZWXb41NLOClbfwsxOcmfWAE5yZVVAxYxSycYIzs8J5wkszqyy3wZlZhTnBmVkllbdsYNmziZiZFcY1ODMrVNKL6ltUM6ssJzgzq6ghPyZiZtVU3qO+TnBmVrgRs6qWmY1E+S06I2k/SX+QdJ+kme3KOsGZWbGU33xw6Vou3wX2B3YADpe0Q6vyTnBmVqjaYyI5zei7K3BfRNwfESuAC4GDWl47IvL5KUogaQnwUAGnngg8UcB5izBIscJgxTtIsUJx8U6KiE3W9sOSriGJLYtxwPN128MRMVx3rvcD+0XEUen2EcDbI+JTzU420J0M6/KltyNpbkRMLeLceRukWGGw4h2kWKF/442I/cq6tm9RzWyQPApsWbe9RbqvKSc4MxskvwPeKOlvJI0FPgRc3qrwQN+iFmi4c5G+MUixwmDFO0ixwuDF27WIWCXpU8AvgVHA2RFxV6vyA93JYGbWjm9RzayynODMrLJGdIKT9KKkeZLulHSFpI3qjl0j6WlJV5YZY71W8UqaImmOpLskzZf0wZLjDEnn1W2PlrSk9l1K2i6N9wVJx5UX6UvxdYr3w+n3ukDSTZJ27uNYD0pjnSdprqTdy4q1H4zoBAcsj4gpEbEj8BQwo+7YN4EjygmrpVbxPgccGRFvAvYDvl2frEuwDNhR0vh0ex/W7Mp/CjgaOK3XgbXQKd4HgHdHxJuBr1JuY36nWK8Ddo6IKcDHgLN6HF9fGekJrt4c4HW1jYi4Dni2vHA6eineiLg3Iv6Yvn8MWAwU8hB0F2YBB6bvDwcuqB2IiMUR8TtgZRmBtdAu3psi4i/p5s0kz16VqV2sS+PlnsMJwIjuRXSC46UBvHvT5nmaftIuXkm7AmOBP/U6rgYXAh+SNA7YCbil5Hg6yRrvx4GrexZVc21jlXSIpHuAq0hqcSPWSE9w4yXNAxYBrwWuLTmeTtrGK2lz4FzgoxGxuoT4XhIR84HJJDWMWWXGkkWWeCW9hyTBfb53kb1Sp1gj4tKI2A44mOSWesQa6QluedpWMYlk0oMZHcqXrWW8kjYk+T/2CRFxc0nxNbqcpJ3tgk4F+0TLeCXtRNKedVBEPNnrwJro+N1GxGxga0lZB7pXzkhPcABExHMkjd7HSur70R2N8aZDVi4FfhIRl5Qb3RrOBr4SEQvKDiSjpvFK2gr4BXBERNxbSmSv1CrWNyidWE3SW4H1gH5IyKXo+1/mXomIOyTNJ6n2nyvpBmA7YANJjwAfj4hflhpknYZ4A9gT2FjStLTItIiYV1Z8ABHxCHB6435JmwFzgQ2B1ZI+A+wQEc/0OMQ1tIoX+DKwMXBmmjtWlT1rR5tYDwOOlLQSWA58sK7TYcTxUC0zqyzfoppZZTnBmVllOcGZWWU5wZlZZTnBmVllOcFVUDrbyIl12w+WMWuHpKnp7BeTe33tTtbmO5I0TdLSwoOz3DjB9YCkc9Jf9JC0UtL9kk6TNKFHIbwNODNLwRH8S7zGd5T+W72/oczPgK17GpWtEz/o2zu/Ipl+aQywB8mwnwnA/25WWNKYiMhlto2IWJLHeaosy3cUEctJHp61AeEaXO+8EBGLIuLPEfFT4HySwdBI2iutMRwg6VZJK4B902Pvk3SbpOclPSDpa+nQLNLjm0r6d0nLJT0k6RWzRzTefkl6taTvSXo8Pe9CSR+UtBfwI2BCXY3zxPQzYyWdKukRSc9J+p2kfRuus5+ke9Jz3gBs0+lLSc97chr7C2nt9ui643tKuiU9539J+peGn/96SWem53hC0uK0djxUV6ar70jSg+nui9Pv4MF0/ytqt5L+l6T7JK1I//5Ew/GQNF3SxZKWpT/fRzp9L5YP1+DKs5ykNlfvVOBY4D7g2TSBnA8cA8wGtgK+TzK+sJawziEZfP8/SCa+/BeSmSaaSscpzgJeA3wUuBfYlmRF8ZuAzwAnA69PP1L7hf5Ruu8fgEeAA4ArJL0tIn4vaUvgMuDfgO+STOPzrQzfw49JarTHAHekP8uWaayvI5ma6FxgWnr9s4DV6fdU82HgO8C7gCnAT4HbeHkgelffEcnt6mLgE8CVwIvNCkk6BDgD+CzwHyT/UzpT0qKIuKKu6JeBmcAXSGYjOVvS7Ih4uE0MloeI8KvgF8kv2JV127sCTwA/S7f3IhlPeljD52YDX2rYdzBJ0hFJDSmA3eqOTyL5hTyxbt+DwHHp+31IEsT2LWKdBixt2Pf69DNbNey/DDgzfX8ySbJU3fEvpvFNbnGtN6bH92tx/GvAH4GhhvheANZPt68H5jR87lrgrPR9199Ruh3A+9t9N8BvSZata/y3vrHhPKfUbY8mSbIfKfu/y5Hwcg2ud/ZLb29Gk9Tc/h34dEOZuQ3buwC7Sqqff2wIGA9sBmxPknhurR2MiIckPdYmjrcAj0fEwi5ifytJQr07HWxesx7w6/T99sDNkf4Wp+Z0OO9bSOL/TYvjtXPWz213I8mEnm8A5qf75jd87jFg07pzdPsdZbU9yawe9W4E/q5h30vxRbKu55K6+KxATnC9MxuYTjJN92PRvANhWcP2EPAV4OImZesbxYueMWEovcbbeOU042U1utf/zI0xBa9sX+7lrBKN18oSnxXAX3LvPBcR90XEQy2SWzO3A9uln2t8rQLuIfk33LX2ASVzl/11m3PeAWwuafsWx1eQrBje+BkBmzWJo7bgyULg7VqziveODj/fvDT+97Q4vhB4R32HAbB7GmPWKdnX5juCJCk1fg/N4tutYd/uwN0ZY7OCOcH1t5OAf5B0kqQdlSy3935J3wCIiD8A1wA/kPROSVNI2oDa1aquI5nD/+eS9pX0N5L2kXRwevxBYFy6b6Kk9SOZ5PF84Jz0+lsreYj3OEmHpp/7PknD/bclbavkGbJPtvvh0vNeBJwl6bA0lj0k1VYzO5MkEZ0paXtJBwJfB86IZNLPjtbyO6p9D3tL2kzSa1qU+SZwhKQZkt4o6dMkHR7fyBKbFc8Jro9FMsHmgSQ1nFvT10ygvvdtGsmydr8GriDpQXywzTlXA/uTNJCfR1IL+Q5JuxYRcRNJsrqA5Db4+PSjHyXpSf0GSa3oSpJJNh9KP/cwcCjJsoW/J+lZnJnhxzwyjfn09LznAK9Oz/loGutbSGp7Z6dx/XOG89abRhffUepYku/9zyQ12FeIiMtI2lE/S1JrOwb4P7FmD6qVyBNemllluQZnZpXlBGdmleUEZ2aV5QRnZpXlBGdmleUEZ2aV5QRnZpXlBGdmlfXfPSM74JsQf74AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.imshow(conf_mat_MLR.mean(0), origin='bottom', cmap='Greens', vmin=0, vmax=5)\n",
    "plt.xticks(range(n_run), ['R1','R2','M1','M2','M3'])\n",
    "plt.yticks(range(n_run), ['R1','R2','M1','M2','M3'])\n",
    "plt.colorbar(ticks=range(6))\n",
    "plt.xlabel('Predicted condition', fontsize=14)\n",
    "plt.ylabel('True condition', fontsize=14)\n",
    "plt.title('MLR + EC', fontsize=14)\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(conf_mat_1NN.mean(0), origin='bottom', cmap='Greens', vmin=0, vmax=5)\n",
    "plt.xticks(range(n_run), ['R1','R2','M1','M2','M3'])\n",
    "plt.yticks(range(n_run), ['R1','R2','M1','M2','M3'])\n",
    "plt.colorbar(ticks=range(6))\n",
    "plt.xlabel('Predicted condition', fontsize=14)\n",
    "plt.ylabel('True condition', fontsize=14)\n",
    "plt.title('1NN + EC', fontsize=14)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}