{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Getting Started: RAM Report Data\n", "\n", "In the process of generating reports for both record-only and stimulation sessions, the RAM reporting pipeline produces a number of summary data files. These files contain all of the information requried to reproduce the plots, summary statistics, and any analyses contained in the reports. However, another use for these summary files is for conducting alternative post-hoc analyses using the RAM data. In particular, some common use cases for this data are:\n", "\n", "1. Post-hoc assessment of classifier performance and/or re-assessment using a different classifier\n", "2. Behavioral analyses comparing cohorts of subjects or comparing stim/no-stim experimental conditions\n", "3. High-level summary information (either behavioral or electrophysiological) across a broad group of subjects\n", "\n", "**Data Types**\n", "\n", "Each of the summary files saved as part of generating a report has a corresponding python class with various properties and helper methods defined. These classes are documented in the [Seraializeable Data Structures](https://pennmem.github.io/ram_utils/html/data.html) section of the ramutils documentation.\n", "\n", "1. Session Summary\n", "2. Math Summary\n", "3. Classifier Summary\n", "4. Target Selection Table \n", "5. Retrained Classifier\n", "\n", "\n", "**Common Patterns**\n", "\n", "The remainder of this guide will go through some common usage patterns for interacting with the saved data including:\n", "\n", "1. Finding data\n", "2. Loading multiple files of a single type\n", "3. Extracting metadata from a file path\n", "4. Building an item-level dataset\n", "5. Reloading a trained classifier\n", "6. Pandas magic\n", "7. Combining across summary file types" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finding Data" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Could not import MorletWaveletFilterCppLegacy (single-core C++ version of MorletWaveletFilter): cannot import name 'MorletWaveletTransform'\n", "You can still use MorletWaveletFilter\n" ] } ], "source": [ "from ramutils.tasks.misc import load_existing_results" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "rhino_root = '/Volumes/RHINO/'\n", "report_db_location = \"/Volumes/RHINO/data10/RAM/report_database/\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For loading all saved data associated with a single subject/experiment combination (either a single session or multiple sessions), you can use the load_existing_results helper function. Internally, this function is used to load data for producing aggregate reports. A dictionary is returned whose keys indicate the type of underlying data. The values are a list of the summary objects, one per session. For this example, we are loading a single session." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "[I 2018-05-16 15:11:47,221 _wrapper.py:21] calling load_existing_results\n" ] }, { "data": { "text/plain": [ "['target_selection_table',\n", " 'classifier_evaluation_results',\n", " 'session_summaries',\n", " 'math_summaries',\n", " 'hmm_results']" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# .compute() is required because load_existing_results is actual a task, which\n", "# is lazily evaluated by default\n", "results_dict = load_existing_results(subject='R1389J', experiment='catFR5', \n", " sessions=[1], stim_report=True, \n", " db_loc=report_db_location, rootdir=rhino_root).compute()\n", "list(results_dict.keys())" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Data is saved as a object.\n", "AUC for this session: 0.6512345679012346\n" ] } ], "source": [ "classifier_data = results_dict['classifier_evaluation_results'][0]\n", "print(\"Data is saved as a {} object.\\nAUC for this session: {}\".format(type(classifier_data), classifier_data.auc))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Data is saved as a object.\n", "Session Completed On: 2018-02-28 16:47:02.190000+00:00\n" ] } ], "source": [ "session_data = results_dict['session_summaries'][0]\n", "print(\"Data is saved as a {} object.\\nSession Completed On: {}\".format(type(session_data), \n", " session_data.session_datetime))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Data is saved as a object.\n", "Subject answer 79.76190476190476% math problems correctly\n" ] } ], "source": [ "math_data = results_dict['math_summaries'][0]\n", "print(\"Data is saved as a {} object.\\nSubject answer {}% math problems correctly\".format(type(math_data),\n", " math_data.percent_correct))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading Multiple Files of a Single Type" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Found 58 FR5 or catFR5 session summaries\n" ] } ], "source": [ "import os\n", "import glob\n", "\n", "from ramutils.reports.summary import FRStimSessionSummary\n", "\n", "fr5_session_summary_locations = glob.glob(os.path.join(report_db_location, '*FR5*session_summary*'))\n", "print(\"Found {} FR5 or catFR5 session summaries\".format(len(fr5_session_summary_locations)))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "session_summaries = []\n", "for summary_loc in fr5_session_summary_locations:\n", " summary = FRStimSessionSummary.from_hdf(summary_loc)\n", " session_summaries.append(summary)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extracting Metadata from the File Path\n", "\n", "As noted in the documentation, underlying files are stored with a specific naming convention so that they can be easily differentiated despite residing in the same directory. ramutils has a utility function for extracting this metadata from the file path, which can be handy when combining multiple files of a single type or across types" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'subject': 'R1275D',\n", " 'experiment': 'FR5',\n", " 'montage': 0,\n", " 'sessions': [1],\n", " 'file_name': 'session_summary',\n", " 'file_type': 'h5'}" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ramutils.utils import extract_report_info_from_path\n", "\n", "metadata = extract_report_info_from_path(fr5_session_summary_locations[0])\n", "metadata" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Building an Item-Level Dataset\n", "One of the properties of the session summary objects is a dataframe containing all of the encoding events from the session that has been enriched with stimulation parameter information and other metadata in the case of a stimulation session. One common use case is to perform some sort of behavioral analysis that requires item-level information. Oftentimes, this analysis is done across subjects and experiments. Luckily, there is a helper method to make it rather trivial to combine this item-level information and return a Pandas Dataframe" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
serialpossessionsubjectexperimentmstimetyperecalledlistis_stim_listphase...is_stim_itemis_post_stim_itemthreshclassifier_outputlocationamplitudepulse_freqstim_durationstimAnodeTagstimCathodeTag
011R1275DFR51497112517408WORD010BASELINE...000.50.565596nannannannannannan
121R1275DFR51497112520131WORD010BASELINE...000.50.647678nannannannannannan
231R1275DFR51497112522816WORD010BASELINE...000.50.358846nannannannannannan
341R1275DFR51497112525335WORD010BASELINE...000.50.459033nannannannannannan
451R1275DFR51497112527751WORD010BASELINE...000.50.446794nannannannannannan
\n", "

5 rows × 21 columns

\n", "
" ], "text/plain": [ " serialpos session subject experiment mstime type recalled list \\\n", "0 1 1 R1275D FR5 1497112517408 WORD 0 1 \n", "1 2 1 R1275D FR5 1497112520131 WORD 0 1 \n", "2 3 1 R1275D FR5 1497112522816 WORD 0 1 \n", "3 4 1 R1275D FR5 1497112525335 WORD 0 1 \n", "4 5 1 R1275D FR5 1497112527751 WORD 0 1 \n", "\n", " is_stim_list phase ... is_stim_item is_post_stim_item \\\n", "0 0 BASELINE ... 0 0 \n", "1 0 BASELINE ... 0 0 \n", "2 0 BASELINE ... 0 0 \n", "3 0 BASELINE ... 0 0 \n", "4 0 BASELINE ... 0 0 \n", "\n", " thresh classifier_output location amplitude pulse_freq stim_duration \\\n", "0 0.5 0.565596 nan nan nan nan \n", "1 0.5 0.647678 nan nan nan nan \n", "2 0.5 0.358846 nan nan nan nan \n", "3 0.5 0.459033 nan nan nan nan \n", "4 0.5 0.446794 nan nan nan nan \n", "\n", " stimAnodeTag stimCathodeTag \n", "0 nan nan \n", "1 nan nan \n", "2 nan nan \n", "3 nan nan \n", "4 nan nan \n", "\n", "[5 rows x 21 columns]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "combined_item_df = FRStimSessionSummary.combine_sessions(session_summaries)\n", "combined_item_df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We previously found 58 sessions worth of data. At 300 encoding events per session, we should have around 17,400 rows in our dataset, unless some sessions were incomplete" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "16548" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(combined_item_df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reloading a Trained Classifier\n", "\n", "Another useful set of information is the classifier itself, rather than just the classifier evaluation results. There are two ways to get this data. The best way is to load the classifier that was actually used during the session. Ramutils has a helpfer function to do this loading for you. The benefit of using this classifier is that it should allow you to exactly match what occured during the session, since this is the classifier that Ramulator was using. If channels were excluded because of some artifact detection routine(s), there will be no need to deal with this post-hoc. Unfortunately, having a consistent storage format is a relatively recent development, so for older subjects it will not be possible to load the actual classifier used. Instead, you can use the classifier trained with the default parameters (frequencies, timings, etc.) that is saved to the report database directory. \n", "\n", "You can find more information about the storage format of the classifier by looking at the [documentation](https://pennmem.github.io/classiflib/html/index.html)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(997, 872)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Reloading the used classifier\n", "from ramutils.classifier.utils import reload_classifier\n", "\n", "classifier_container = reload_classifier('R1387E', 'catFR5', 1, mount_point='/Volumes/RHINO/')\n", "classifier_container.features.shape # n_events x n_features power matrix" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(997, 968)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Loading the classifier from the reports database\n", "from classiflib.container import ClassifierContainer\n", "\n", "retrained_classifier_loc = glob.glob(os.path.join(report_db_location, 'R1387E_catFR5_all_retrained_classifier*'))[0]\n", "classifier_container = ClassifierContainer.load(retrained_classifier_loc)\n", "classifier_container.features.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, we can see that the classifier used during the session differs from the one that does not take into account artifactual channels. In fact, we can deduce that 12 channels were exclude due to artifact detection (12 * 8 = 968 - 972). \n", "\n", "Once a classifier has been loaded, the main use cases are to make out of sample predictions and check on the AUC from the training session." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Training AUC: 0.61\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/zduey/miniconda3/envs/ramutils_dev/envs/ramutils_test/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n", " warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4m9WBLvD3SLLkRfIqy/tux44TEoc4JCQEAiSUAoWyFmh7YbpAO9AyLXfa6S3TmWk7t9PbjbnTdi60ZaAdwlIKlLIvSUjI7uyLkzi25X2R5V3etJz7h500gBPLtqSjT3p/z5Nn7Ea1329kvz053znfEVJKEBGRduhUByAiotlhcRMRaQyLm4hIY1jcREQaw+ImItIYFjcRkcawuImINMbgz4uEEHYAQwC8ADxSyupghiIiovPzq7inXCml7AlaEiIi8stsittvVqtVFhYWBuNLExFFpH379vVIKdP9ea2/xS0BvC2EkAAek1I+fqEXFxYWoqamxs8vTUREQogmf1/rb3GvkVK2CyFsAN4RQpyQUm79yDe9D8B9AJCfn+93WCIimh2/VpVIKdun/m83gJcAXDLNax6XUlZLKavT0/0a7RMR0RzMWNxCiAQhhOXMxwCuAXA02MGIiGh6/kyVZAB4SQhx5vUbpZRvBjUVERGd14zFLaVsALA0BFmIiMgP3DlJRKQxLG4iIo1hcRMRaUxQdk4SBdK4x4uTnUM43j6Ihh4XRiY8GHf7kJJgREFaPCoyE7EsLxk6nVAdlSgkWNzkl427m0P6/TxeH453DOJo+yBOdQ1hwuMDABh0AiaDDga9Dq5xDzy+ycOuE2MNuCgnCatLrUiJN5736969kpvDSPtY3BRWhsbc2NXgxJ7GXrgmvDCbDFiam4wymxnZyXFIjo+BbnJpKnxSYnDUDbtzBEfaBrCrsRe7G3uxusSKdeXpiI3RK74aouBgcVNYGHN7sbXOge2ne+DxSpRnWnBpSRpK0s1ni/qjdEIgOd6IqngjqvKS0T8ygXeOd2FrnQOHW/vx2ZUFyEmJC/GVEAUfi5uUklLiUGs/Xj3cgZEJL5bkJmF9RQasFtOsv1ZyvBG3V+dhZXEant3TjMe21uPGpdmoLkwNQnIidVjcpMzAqBsv7m9FXfcw8lLi8DdrcpCTPP8Rcn5qPB64shTP7W3BiwfaMOr2Ym0Zn59DkYPFTUqc6hrC8zUt8HglbliShVXFaeedEpmLBJMB96wuxPM1LXjjaCeEELis1Bqwr0+kEoubQkpKiU0nu7Gpthu2RBPuuiQfNktsUL6XXidwR3UepJR4/UgHYvSCq0ooInADDoWM1yfx0oE2vFfbjaq8ZHz1itKglfYZep3AZ1bkozzDgr8caseuBmdQvx9RKLC4KSTcXh+e3t2EmqY+XFmejtuW58JoCM2P32R55yE1wYS/fXo/WvtGQvJ9iYKFxU1B5/VJPLe3BSc6h3Dj0mxsqMyECOB8tj9iY/T4/KoCuL0+3P+HfRj3eEP6/YkCicVNQeWTk9MjxzsGz96EVCXdYsIv7qjCsfZB/PztU8pyEM0Xi5uC6r3aLuxv7sPVFTasLlG/qmN9ZQbuXpmPx7c1YGc957tJm1jcFDS1HYPYfNKB6oIUXFVhUx3nrEeuX4jCtAQ8/PxBDIy6VcchmjUWNwWFc3gcf9zXguzkWHxqaXbI57QvJN5owC8+U4XOwTH8+M0TquMQzRqLmwLO4/Nh455mCAjcfUkBYvTh92NWlZeML6wpwsbdzaix96qOQzQr4fcbRZq39ZQDHQNjuPXiXKQmnP8Rq6p9Y8MC5CTH4X+9dOTsY2OJtIDFTQHVNTiGzSccWJKbhMrsRNVxLijBZMD3b1qEU13D+M22BtVxiPzG4qaA8UmJF/e3whSjww1LslXH8cvVCzPwiUUZ+OWm0+gaHFMdh8gvLG4KmBp7H1r6RnH9RVkwm7TzGJzvXlcJr0/iJ2+dVB2FyC/a+e2isDbu8eK92i4UpMajKi9ZdZzzOt8RbCuLU/Gnfa3ISIyd96Nl+SArCjaOuCkgPjjdg6FxDz65OPTb2QPhynIb4ox6vHa4A1JK1XGILojFTfM2NObGtlM9WJydiPy0BNVx5iQ2Ro/1CzNgd7pwrH1QdRyiC2Jx07xtOtENj8+HaxZlqo4yLysKU2GzmPDmsU54vFweSOGLxU3zMjTmxr6mPiwvSIHVPPtzIsOJXidw3UVZ6HVNYCef201hjMVN87Kj3gmvT0bMmY4LMixYkGHGphPdGB73qI5DNC0WN83ZmNuL3Y1OLMpO1Pxo+1zXLc6C2+vDphPdqqMQTYvFTXO2196LMbcPly+IjNH2GbbEWCwvSMXexl70uiZUxyH6GBY3zYnH58P20z0oTk9Abkq86jgBd1WFDUJMPk+cKNywuGlOjrcPYnDMg7Wl6g9HCIakuBhcWpKGgy396ORWeAozLG6akz32XiTHx6Asw6I6StBcsSAdphgd3jnWqToK0Yf4XdxCCL0Q4oAQ4tVgBqLw1zM8jgaHCysKU6HT4C5Jf8UbDVhblo7aziE0OV2q4xCdNZsR90MAaoMVhLRjr70XOgEsz09RHSXoVpekwWwy4K1jXdwKT2HDr+IWQuQCuB7Ab4Mbh8Kdx+vDvqY+VGQmIjEuRnWcoDMZ9Liywga704W67mHVcYgA+D/ifhTAtwBwH3CUO94xiJEJLy4pSlUdJWRWFKYgJT4Gbx3rhI+jbgoDMxa3EOIGAN1Syn0zvO4+IUSNEKLG4XAELCCFlwPN/UiKi0Gpzaw6SsgYdDqsX5iBjoExHG0bUB2HyK8R9xoANwoh7ACeBXCVEOK/P/oiKeXjUspqKWV1enpkbcigScPjHtR1D2FpblJE35ScztK8ZGQmxuKd413w+jjqJrVmLG4p5XeklLlSykIAdwLYJKX8XNCTUdg52jYAn5wssWijEwIbKjPgdE2gpomnwpNaXMdNfjvc2g+bxYTMxFjVUZSoyLQgPzUem09081R4UmpWxS2l3CKlvCFYYSh89Y9MwO4cwdK8ZE2ecBMIQgh8YlEmBsc82MXHvpJCHHGTXw63Tt6UW5obfdMk5yqyJmBBhhnvn3JgdMKrOg5FKRY3+eVQaz/yUuKQmmBUHUW5ayozMer2YlsdV0+RGixumlGT04WOgTFclJOkOkpYyE6Ow5LcJGyv78HQmFt1HIpCLG6a0TvHJx9tWpnN4j5jw8IMeH0Sm0/ysAUKPRY3zejtY13ITIzlNMk50swmVBemYg8PWyAFWNx0QT3D46hp6kVldqLqKGHnqnIb9DqBd3nYAoUYi5suaFNtN3wSqMxicX9UYlwMLi224lBLPzoGRlXHoSjC4qYLevt4J3KS45CVFJ2bbmZy9rCF4xx1U+iwuOm8XOMebK3rwTWLMqJ2081M4ox6XF6WjhOdQ7D38LAFCg0WN53XtroeTHh82FCZoTpKWFtdYoXFZMDbxzt52AKFBIubzmvLyW5YTAasKIyeZ2/PhdGgw7oKG+zOER62QCHB4qZpSSmx5aQDa0qtiNHzx2QmZw5bePt4J3x87CsFGX8jaVqnuobROTiGdeV8tro/DDodrl6Ygfb+MbzJU+EpyFjcNK0tUzsCr2Bx+60qLxnpFhN+9vZJeLx87CsFD4ubprXlpAMVmRZkJcWpjqIZOiGwYWEG6h0uvHSgTXUcimAsbvqY4XEPapp6Odqeg0XZibgoJwmPvluHcQ8f+0rBweKmj9l+ugdur8QVC1jcsyWEwN9/ohxt/aN4Znez6jgUoVjc9DFbTjqQYNSjuoDLAOdibZkVq4pT8cvNpzEy4VEdhyIQi5s+5oPTDlxaYoXRwB+PuTgz6u4ZnsB/bberjkMRiL+Z9CEtvSNo6R3FZaVpqqNo2vKCVFxdYcNj79djYISHLVBgGVQHoPCy/XQPAGBNqVVxEu3aODW3XZmdiPdOdOOhZw/gmkWZAf0ed6/MD+jXI23hiJs+5IPTPbBZTCi1mVVH0bysJB5xRsHB4qazfD6JnfVOrCm18mmAAbJ+6oizLSd5sDAFDoubzjrZNQSna4LTJAFkNZuwLD8Fe+29GOSomwKExU1n/XV+mzcmA2ndgnT4pMS2Uxx1U2CwuOms7ad7UGxN4Db3AEszm1CVl4w99l7OdVNAsLgJAOD2+rCnsZfTJEGyboENHq/EB3U9qqNQBGBxEwDgcOsAXBNerC7hNEkwWC0mLM1Lxq5GJ4bHuZuS5ofFTQCAXQ1OAMAlRdzmHizrFqTD45Vn7yUQzRWLmwBMFnd5hgVpZpPqKBHLlhiLxTlJ2NngxAhH3TQPLG6C2+vDvqY+rCzmaDvYrqywYcLjw/Z6jrpp7ljchCNtAxiZ8GJVMee3gy0zMRaLsxOxo96J0Qk+r5vmhsVN2N3QC4Dz26FyZYUN4x4fdjRw1E1zw+Im7GpwosxmhpXz2yGRlRSH8gwLdtU74ebZlDQHLO4o5/H6UGPv5TRJiK0ts8I14cX+5j7VUUiDZixuIUSsEGKPEOKQEOKYEOJfQhGMQuNo+yBcE17emAyxImsCcpLj8EFdD3xSqo5DGuPPiHscwFVSyqUAqgBcK4RYFdxYFCpn1m+vLOKIO5SEELh8QTqcrgkcbx9UHYc0ZsbilpOGpz6NmfrDIUKEqLH3otiagHQL57dDbVF2IlITjNhW54DkqJtmwa85biGEXghxEEA3gHeklLunec19QogaIUSNw8GnoGmBzydR09SHFYWcJlFBJwTWlFrR0jeKJueI6jikIX4Vt5TSK6WsApAL4BIhxOJpXvO4lLJaSlmdnp4e6JwUBKcdw+gfcaO6MEV1lKi1PD8F8UY9ttVxsEP+m9WqEillP4AtAK4NShoKqb32yfXbHHGrYzTosKo4DbWdQ+geGlMdhzTCn1Ul6UKI5KmP4wCsB3Ai2MEo+GrsfbCaTShIi1cdJaqtKk6DQSf4yFfymz8j7iwAm4UQhwHsxeQc96vBjUWhsNfeixWFKTxfUjGzyYCLC1JwoKWfBy2QX/xZVXJYSrlMSrlESrlYSvn9UASj4OoYGEVr3yinScLEZaVW+HwSuxt7VUchDeDOySi11z65Y4/FHR6sZhMWZFiwp7EXHm6DpxmwuKNUjb0X8UY9FmZZVEehKatL0jA87sGRtgHVUSjMsbij1F57Hy7OT4FBzx+BcFFqMyPdbMKOeic35NAF8bc2Cg2OuXGic5DTJGFGCIFLS9LQ1j+Klr5R1XEojLG4o9C+pj5ICazgxpuwsyw/GSaDDjt4Qg5dgEF1AJq/jbubZ/X6t491QieAU13DsHOrdVgxGfSoLkjBzgYnBhe7kRgXozoShSGOuKOQ3TmC7OQ4GA18+8PRquI0SAkuDaTz4m9ulPF4fWjtG0FBKndLhqs0swnlmRbsaXRyaSBNi8UdZdr6R+HxSRRaE1RHoQtYXTJ5Qs5hLg2kabC4o8yZx4cWpLG4w1lJ+uQz0nfU93BpIH0MizvK2J0uWM1GmE28Lx3OhBBYXZKG9v4xtPTyBjJ9GIs7ivikRJNzhKNtjajKS0ZsjA47p46XIzqDxR1FHEPjGHV7Ucji1gSTQY/l+Sk42jbIpwbSh7C4o4jd6QIAFPL525qxsjgNXimxx86lgfRXLO4o0uQcgdlkQGqCUXUU8tPkUwPN2NPYC6+PNylpEos7ijQ5XShIi+fBCRqzqjgNQ2MeHGvn0kCaxOKOEgOjbvSNuDm/rUELMixITTDyJiWdxeKOEn+d32Zxa41OCKwqSkWTcwTt/XxqILG4o0aT0wWjXofMpFjVUWgOlhekIkYvsIujbgKLO2o0OUeQnxoPvY7z21oUZ9SjKi8Zh1r7MTLhUR2HFGNxR4ExtxedA2Mo4DJATVtVnAa3V2JfU5/qKKQYizsKNDlHIAE+WErjspLiUJiWgF0NTi4NjHIs7ijQ5HRBJ4C8FI64te7SkjT0jbix5WS36iikEIs7CtidLh6cECEqsxKRGGvAkzvsqqOQQvxNjnBurw8tfaNcBhgh9DqBS4rSsK2uB/WOYdVxSBEWd4Rr7RuF1ydRxPntiLGiMAUxeoE/7GxSHYUUYXFHuDMbb3hUWeSwxMbg+ouy8MK+VgyPc2lgNGJxRzh7jwsZiSbE8+CEiHLP6kIMj3vwx5oW1VFIARZ3BPP6Jg9O4Px25FmWn4LlBSl4YnsjlwZGIRZ3BOsYGMWE18f57Qj15bVFaOkdxVvHOlVHoRBjcUcwew8fLBXJNlRmoiAtHr/Z1qA6CoUYizuCNTpHkJpgRGJcjOooFAR6ncAX1hThQHM/9jXxhJxowuKOUD4pYe9xoYij7Yh2e3UukuJi8JutjaqjUAixuCNU95mDgTm/HdHijQZ8dmU+3jreiaappZ8U+WYsbiFEnhBisxCiVghxTAjxUCiC0fycmd/mjcnId8/qQhh0Ak98wFF3tPBnxO0B8LCUciGAVQAeEEJUBjcWzZfd6UJirAEp8ZzfjnQZibG4cWkOnq9pRf/IhOo4FAIzFreUskNKuX/q4yEAtQBygh2M5k5OzW8XWhN4MHCU+NLaIoy6vXh6d7PqKBQCs5rjFkIUAlgGYPc0f3efEKJGCFHjcDgCk47mpNc1gcExD5cBRpGFWYlYW2bFUzvsGPd4VcehIPO7uIUQZgB/AvB3UsrBj/69lPJxKWW1lLI6PT09kBlpls48n4Tz29HlvsuL0T00jj/ta1MdhYLMr+IWQsRgsrSfllK+GNxINF/2nhHEG/VIt5hUR6EQuqzUiqV5yfjP90/D4/WpjkNB5M+qEgHgdwBqpZQ/D34kmq9GpwuFaQnQcX47qggh8LUrS9HSO4pXDrWrjkNB5M+Iew2AzwO4SghxcOrPdUHORXM0MOpGr2sChTwYOCpdvdCGikwLfrX5NB8+FcH8WVXygZRSSCmXSCmrpv68HopwNHtn5re58SY6CSHw4FWlqHe48MbRDtVxKEi4czLCNDpcMBl0yEqKUx2FFPnk4iyU2sx49N06jrojFIs7wtQ7hlFkTYBex/ntaKXXCXxj/QKc7h7GXzjXHZFY3BGkf2QCTtcEitPNqqOQYp9cnImKTAseffcUV5hEIBZ3BGmYej5JSTrnt6OdTifw8DXlsDtH8OJ+ruuONCzuCNLgGEa8UY+MxFjVUSgMrF9ow9LcJPz7e3UYc3M3ZSRhcUcIKSXqHS4Up5u5fpsATK4w+fa1FWjrH8Xvd9pVx6EAYnFHCKdrAgOjbk6T0IesLrViXXk6frnpNJ8cGEFY3BGi3jEMACix8sYkfdg/fLICw+Me/GrzadVRKEBY3BGiwTH5/O00s1F1FAozFZmJuG15Lp7a0YSW3hHVcSgAWNwRwCcl6h3DKEk38/nbNK1vbiiHXifwr6/Vqo5CAcDijgDt/aMYmfCi1MZpEppeZlIsHryqFG8e68QHdT2q49A8sbgjQF335Px2WYZFcRIKZ1+8rAj5qfH4l78cg5ubcjSNxR0B6rqGkJ0UC7PJoDoKhbHYGD0euX4h6rqH8YedTarj0DywuDVuaMyN5t4RjrbJLxsqM3D5gnT8/J1T6BwYUx2H5ojFrXE76p3wSaCM89vkByEEfnDTIri9PvzLX46pjkNzxOLWuK2nHDAadMjnwQnkp4K0BHz96jK8cbQTm050qY5Dc8Di1jApJbbWOVBsTYBBx7eS/PfltcUos5nxjy8fw/C4R3UcmiX+tmuY3TmClt5Rzm/TrBkNOvzbrRehfWAUP3qda7u1hsWtYZtOdAMAylncNAfLC1LxpcuK8PTuZmyrc6iOQ7PA4taw92q7UGYzIzWB29xpbh6+phwl6Qn41guHMTjmVh2H/MTi1qiBUTf2NPZifWWG6iikYbExevzsjip0DY7hh68eVx2H/MTi1qj3Tzng8UmsX2hTHYU0riovGV+5ogTP17RylYlGcKudRr1X24XUBCOq8lJwsnNYdRwKsY27mwP69TITY5GZGIuHnj2Ih64uQ7zRgLtX5gf0e1DgcMStQW6vD5tPdOOqChtPc6eAMOh1uHV5LlzjHp4MrwEsbg2qsfdhcMzDaRIKqJzkOFxVYcOh1gEcaO5THYcugMWtQe/VdsGo12FtWbrqKBRh1pXbUJiWgD8faoe9x6U6Dp0Hi1tjpJR442gn1pSmIYFPA6QA0wmBO6pzoRcCX3/2ACY8fPxrOGJxa8yh1gG09Y/i+iXZqqNQhEqON+KWi3NwuHUAP3v7pOo4NA0Wt8a8drgdMXqBDVy/TUG0KDsJn1uVj8e2NmDrKe6qDDcsbg2RUuL1I51YW5aOpLgY1XEowj1yfSUWZJjxzecPwTE0rjoOnYPFrSEHW/rR1j+K6y7KUh2FokBsjB7/cdfFGBpz45vPH4TXJ1VHoiksbg15/UgHp0kopMozLfjnGxdhW10PfrnptOo4NIXFrRGcJiFV7lyRh1uW5eDR907xhPgwweLWiL32PrT1j+KGJZwmodASQuCHNy9GaboZDz17AF2DPKtStRmLWwjxhBCiWwhxNBSBaHp/rGmB2WTAtYszVUehKBRvNOA/P3cxRt1efG3jAXi8XN+tkj8j7icBXBvkHHQBrnEPXjvSgesvykK8kZtuSI1SmwU/uuUi7LH34qdvn1IdJ6rNWNxSyq0AekOQhc7jzaOdGJnw4rbqXNVRKMrdVJWDu1fm4/+9X4/3avkIWFU4x60BL+xrRWFaPKoLUlRHIcL3bqjEouxEfPP5Q2jpHVEdJyoFrLiFEPcJIWqEEDUOB3daBUpL7wh2Njhx2/JcCMFHuJJ6sTF6/PqzF8MnJR7cuB/jHq/qSFEnYMUtpXxcSlktpaxOT+dT6wLl+ZoWCAHcfDGnSSh8FKQl4Ce3LcWh1gH879d4SnyocaokjI17vHhmTzOurrAhJzlOdRyiD7l2cSa+eFkRntrZxMMXQmzGJQpCiGcArANgFUK0AvgnKeXvgh2MgNcOd6BneAL3rC5UHYWikD/HoxWmJSA/NR4P//EQGh0uWC0mv78+j0abO39WldwlpcySUsZIKXNZ2qEhpcSTO+wotZlxWalVdRyiael1AneuyINBJ7BxTzOf3x0inCoJUwda+nG4dQD3XFrAm5IU1pLjjbijOg9dg2N4+WAbpOTDqIKNxR2mntxuh8VkwC28KUkasCDDgqsX2nCwpR87G5yq40Q8FncYanK68NqRDnxmRR6PJyPNWFduQ0WmBa8f6eB5lUHG4g5Dv95cD71O4L7Li1VHIfLb5HmVeUiJN+KZPc0YHHWrjhSxWNxhpqV3BH/a34q7L8mHLTFWdRyiWYmN0eOzqwow7vFh455meHy8WRkMLO4w85/v10MnBO6/gqNt0qbMxFjccnEOmntH8PqRDtVxIhKLO4y09Y/ijzUtuGNFLrKSuOGGtGtJbjIuK7ViV0Mv9jf1qY4TcVjcYeTHb5yATgh8dV2p6ihE8/aJRZkosibg5YNtaOsfVR0norC4w8S+pl68cqgd911ezO3tFBH0OoG7LslHgsmAp3c3YWTcozpSxGBxhwGfT+L7r9bCZjHhK1eUqI5DFDBmkwF3X5KPoTEPnubNyoBhcYeBlw+24VBLP759bQXXbVPEyUuNxy3LctDY48IrB9u5szIA2BKK9QyP44ev1WJpXjJuXpajOg5RUCzLT4FjaBxbTjlgs5hwWRkf/TwfLG7Fvvfnoxge8+Anty2BTsdnklDkWl+ZAcfwON442gmr2f+nCNLHcapEoVcPt+P1I514aH0ZFmRYVMchCiqdELh9eR6ykmPxbE0LajsGVUfSLBa3Ih0Do/jen49hSW4S7ufWdooSRoMOn19VCJNBhy89VYPuwTHVkTSJxa3AhMeHv316P8bdXvz8jioY9HwbKHokxcXg86sK0Dcygf/xxB4MjPCZJrPFxlDgh68dx4Hmfvzk9qUotZlVxyEKudyUeDz2+eWodwzji0/txegEDxyeDRZ3iD23txm/39mE+y4vxnUXZamOQ6TM2rJ0PPqZZdjX3IcHNu6H28s13v7iqpIQOHN2X23HIP57VxPKbGbkpcT7daYfUSS7fkkW+kcX47svHcW3XjiMn92+lKur/MDiDpHGHhee2dOMnJQ43L0yH3r+cBIBAD67sgB9rgn89O1TSI6PwfduqORxfTNgcYdAQ88wfr+zCcnxRtxzaSFMBr3qSERh5YErS9HrcuOJ7Y3QC4HvXr+Q5X0BLO4g23yiG09utyMlwYgvrCnilnaiaQgh8I83LIRPSvz2g0Z4peTI+wLYIkH03N5mPPLyUWQkxuLe1YUsbaILEELgnz5VCZ0QeGJ7I8bcXvzgpsVcLjsNNkkQeLw+/PC1Wjy5w461ZVZcWW5DbAynR4hmcmbkHW/U45ebT8M5PIH/e9cy/v58BP+nLMDa+kdx929348kddnzpsiL8170r+ENHNAtCCPzPT5Tjnz5VibePd+Hzv9sN5/C46lhhhcUdQK8d7sAnH92KY20D+MVnluKRGyr5zzyiOfqbNUX4j7uW4VDrAG761Xac7BxSHSlssFUCoHNgDPf/oQYPbNyPonQzXn9oLW5elqs6FpHmfWppNp6//1KMe3y45dfb8ZdD7aojhQUW9zxMeHz47bYGbPj5+9hy0oFvX1uBF75yKQrSElRHI4oYVXnJeOXBNSjPtOBrzxzAd186gjF3dG+R583JOZBS4q1jXfjxmyfQ2OPC2jIrfnDTYhRaWdhEwZCVFIfn7r8UP33rJB7b2oA9jb34ye1LUZWXrDqaEizuWZBSYtOJbvzi3VM42jaIkvQE/Ne9K7CuPJ3rTYmCLEavw3euW4hLS9LwnReP4JZfb8eXLy/G168qi7qlttF1tXM04fHhlUPt+M3WBpzsGkJ+ajx+evtSfLoqmzcfiUJsXbkNb33jcvzrq7V47P0GvHygDd/55ELcVJUdNQMoFvcFnOoawnN7W/DSgTb0uiZQkWnBT29fipuqshHDwiZSJjE2Bj++bQnuWJGHf37lGP7uuYN4fGsDHlpfhmsqMyK+wFncH9HnmsCbxzrx3N4WHGzpR4xeYENlBu5ckY+1ZdaI/4Eg0pLlBSn48wNr8NKBNvzHpjrc/4d9KM+w4N41hfh0VQ7lEQgvAAAGk0lEQVTijJG5hyLqi1tKiXrHMN6t7cZ7tV3Y19QHnwQWZJjxyPULcfOyHKTxYFOisKXTCdy6PBc3VWXj5YPt+N0HjfjOi0fwo9drccPSbNy8LAfL81Mi6nGxfhW3EOJaAP8OQA/gt1LKfwtqqiDy+iROdg5hf3Mf9jf3ocbeh+beEQBAZVYiHryyFBsqM7E4J5GjayINMeh1uG15Lm69OAd7GnuxcU8zXtrfho27m5FuMeHK8nSsK7fhkqJUzZ8yP2NxCyH0AH4FYAOAVgB7hRCvSCmPBzvcfHh9Em19o6jvGUaDw4V6xzDqu4dxtG0ArqljktISjFiWn4IvX16MqytsyE6OU5yaiOZLCIGVxWlYWZwG17gH7xzvwju1XXjjaCeer2kFABSnJ2BRdhIW2MxYkGlBeYYFeanxmnlOvj8j7ksAnJZSNgCAEOJZADcBCHhxD4y44fH54PVJuH0SXq88+7nHJzHm9sI17sXwuBvD4164xj0YnvrTP+KGY2gMXYPj6B4ag2NoHD7516+dHB+DYmsCbl2ei4vzU3BxfgryUuM4qiaKYAkmAz69LAefXpYDt9eHw60D2GvvRY29Dwdb+j60E9Oo18GWaEJGYiwyEk2wWWKRbjHBbDLAEmuA2WSAOdaAuBg9YvQ6GPQCMXodYnR//dio1yEpPibo1+VPcecAaDnn81YAK4MRZuWP3sWYe/bnzhl0AklxMbBN/T98YZYFGYmxyEmOQ4nNjJJ0M1ITjEFITERaEaPXYXlBCpYXpABXTP5nrnEP6rqHcaprCPWOYXQPjqNrcAwnO4ew7VQPhsY9s/oeVrMRNY9sCEL6D/OnuKcbksqPvUiI+wDcN/XpsBDi5HyCBYgVQI/qEEHGa4wMUXeNn1UYJEisTUCP+Mc5//cL/H2hP8XdCiDvnM9zAXzsSS9SyscBPO7vNw4FIUSNlLJadY5g4jVGBl6j9oXy+vzZRbIXQJkQokgIYQRwJ4BXghuLiIjOZ8YRt5TSI4R4EMBbmFwO+ISU8ljQkxER0bT8WsctpXwdwOtBzhIMYTV1EyS8xsjAa9S+kF2fkPJj9xmJiCiM8UlJREQaExHFLYS4VghxUghxWgjxD9P8/VeEEEeEEAeFEB8IISpV5JyPma7xnNfdJoSQQgjN3b334328VwjhmHofDwohvqQi51z58x4KIe4QQhwXQhwTQmwMdcb58uM9/MU5798pIUS/ipzz4cc15gshNgshDgghDgshrgt4CCmlpv9g8oZpPYBiAEYAhwBUfuQ1ied8fCOAN1XnDvQ1Tr3OAmArgF0AqlXnDsL7eC+AX6rOGsTrKwNwAEDK1Oc21bkDfY0fef3XMLnYQXn2AL+PjwP46tTHlQDsgc4RCSPus1vypZQTAM5syT9LSjl4zqcJmGYDUZib8Rqn/ADA/wEwFspwAeLvNWqVP9f3ZQC/klL2AYCUsjvEGedrtu/hXQCeCUmywPHnGiWAxKmPkzDNvpf5ioTinm5Lfs5HXySEeEAIUY/JYvt6iLIFyozXKIRYBiBPSvlqKIMFkF/vI4Bbp/75+YIQIm+avw9X/lzfAgALhBDbhRC7pp7KqSX+vocQQhQAKAKwKQS5Asmfa/xnAJ8TQrRicjXe1wIdIhKK268t+VLKX0kpSwB8G8AjQU8VWBe8RiGEDsAvADwcskSB58/7+BcAhVLKJQDeBfBU0FMFjj/XZ8DkdMk6TI5GfyuE0NJpuH79Lk65E8ALUkqtHdfuzzXeBeBJKWUugOsA/GHqdzRgIqG4/dqSf45nAXw6qIkCb6ZrtABYDGCLEMIOYBWAVzR2g3LG91FK6ZRSjk99+hsAy0OULRD8+TltBfBnKaVbStkI4CQmi1wrZvO7eCe0N00C+HeNXwTwPABIKXcCiMXkc1oCR/VkfwBuFhgANGDyn11nbhYs+shrys75+FMAalTnDvQ1fuT1W6C9m5P+vI9Z53x8M4BdqnMH+PquBfDU1MdWTP6TPE119kBe49TrygHYMbWPREt//Hwf3wBw79THCzFZ7AG9Vs0fXSbPsyVfCPF9TBb0KwAeFEKsB+AG0AfgHnWJZ8/Pa9Q0P6/x60KIGwF4APRicpWJJvh5fW8BuEYIcRyAF8DfSymd6lLPzix+Tu8C8KycajYt8fMaHwbwGyHENzA5jXJvoK+VOyeJiDQmEua4iYiiCoubiEhjWNxERBrD4iYi0hgWNxGRxrC4iYg0hsVNRKQxLG4iIo35/0ucWOJnIMJvAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import seaborn as sns\n", "\n", "%matplotlib inline\n", "\n", "classifier = classifier_container.classifier\n", "\n", "print(\"Training AUC: {:.2f}\".format(classifier_container.classifier_info['auc']))\n", "\n", "# Create 100 random new events and calculated predicted probabilities\n", "out_of_sample_powers = np.random.normal(size=(100, classifier_container.features.shape[1]))\n", "predicted_probs = classifier.predict_proba(out_of_sample_powers)[:, 1]\n", "sns.distplot(predicted_probs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pandas Magic\n", "\n", "One pattern that has come in handy for almost every cross-subject analysis, is the ability to embed these summary objects in a pandas dataframe. Bear with me as this sounds strange, but what it allows you to do is to write your own functions that operate on a summary object and return some information. These functions can then make use of the pandas API to quickly apply the transformation to all session summaries in a series. In this way, you can easily extend the built-in set of methods that operate on these summary objects. Let's go through a quick example.\n", "\n", "First, we will define some custom function that takes a summary object as its input. In this case, we have some behavioral measure of recall performance that we want to calculate that is not part of the standard set of methods." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def normalized_item_delta_recall(summary):\n", " \"\"\" Calculate normalized delta recall at the item level \n", " \n", " Parameters\n", " ----------\n", " summary: ramutils.reports.summary.FRStimSessionSummary\n", " Helper class containing events, powers, montage, and other information needed to analyze a single\n", " stimulation session\n", " \n", " Returns\n", " -------\n", " normalized_delta_recall: float\n", " Item-level measure of the change in recall performance due to stimulation. Calculated by\n", " comparing the recall rates for stimulated items with items for non-stim lists that would\n", " have been stimulated\n", " \n", " \"\"\"\n", " if type(summary) != list:\n", " summary = [summary]\n", " df = FRStimSessionSummary.combine_sessions(summary)\n", " nonstim_low_bio_recall = df[(df.classifier_output < df.thresh) &\n", " (df.is_stim_list == False)].recalled.mean()\n", " stim_recall = df[df.is_stim_item == True].recalled.mean()\n", " overall_recall = df.recalled.mean()\n", " normalized_delta_recall = 100 * ((stim_recall - nonstim_low_bio_recall) / overall_recall)\n", " return normalized_delta_recall" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we will calculate normalized item-level delta recall using two different methods. This first method uses list comprehension, and the second uses the session summary objects embedded in a pandas dataframe." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "normalized_delta_recall = [normalized_item_delta_recall(summary) for summary in session_summaries]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, so great. Now we have this list of values that we wanted to extract. We can do some simple things with it:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Average item-level delta recall: 4.24%\n", "Biggest item-level delta recall: 87.69%\n", "Biggest adverse item-level delta recall: -60.29%\n" ] } ], "source": [ "print(\"Average item-level delta recall: {:.2f}%\".format(np.nanmean(np.array(normalized_delta_recall))))\n", "print(\"Biggest item-level delta recall: {:.2f}%\".format(max(normalized_delta_recall)))\n", "print(\"Biggest adverse item-level delta recall: {:.2f}%\".format(min(normalized_delta_recall)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This works pretty well for doing some basic manipulations and even plotting, but what if we want to ask some slightly less trivial questions? For example, what is the mean, min, max delta recall split by FR5/catFR5? What if we want to split based on stimulation region or some other value? This is no longer straightforward to do with the list-based method. Instead, let's get these summary objects into a pandas dataframe and really let the analyses run wild." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "\n", "session_data = {\n", " 'subject': [],\n", " 'experiment': [],\n", " 'session': [],\n", " 'session_summary': [],\n", " 'montage': [],\n", " 'anode': [],\n", " 'cathode': [],\n", " 'electrode_type': [],\n", " 'location': [],\n", " 'amplitude': [],\n", " 'pulse_freq': [],\n", " 'stim_duration': [],\n", " 'date': []\n", "}\n", "\n", "for session_summary_loc in fr5_session_summary_locations:\n", " metadata = extract_report_info_from_path(session_summary_loc)\n", " session_summary = FRStimSessionSummary.from_hdf(session_summary_loc)\n", " session_df = pd.DataFrame.from_records(session_summary._events)\n", " stim_params = FRStimSessionSummary.stim_parameters([session_summary])[0]\n", " \n", " session_data['subject'].append(metadata['subject'])\n", " session_data['experiment'].append(metadata['experiment'])\n", " session_data['session'].append(metadata['sessions'][0])\n", " session_data['session_summary'].append(session_summary)\n", " session_data['montage'].append(metadata['montage'])\n", " session_data['anode'].append(stim_params['stimAnodeTag'])\n", " session_data['cathode'].append(stim_params['stimCathodeTag'])\n", " session_data['electrode_type'].append('')\n", " session_data['location'].append(stim_params['location'])\n", " session_data['amplitude'].append(float(stim_params['amplitude']))\n", " session_data['pulse_freq'].append(float(stim_params['pulse_freq']))\n", " session_data['stim_duration'].append(float(stim_params['stim_duration']))\n", " session_data['date'].append(session_summary.session_datetime)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentmontagesessiondatesession_summaryelectrode_typelocationpulse_freqstim_duration
0R1275DFR5012017-06-10 17:15:09.603000+00:00<FRStimSessionSummary(_events=[ ( 1, 1, 'R1275...Right Cerebral White Matter200.0500.0
1R1275DFR5022017-06-11 17:29:07.115000+00:00<FRStimSessionSummary(_events=[ ( 1, 2, 'R1275...Right Cerebral White Matter200.0500.0
2R1292EFR5002017-07-03 16:12:42.221000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1292...Left Middle Temporal Gyrus100.0500.0
3R1304NFR5002017-05-29 15:35:00.879000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1304...Right Cerebral White Matter200.0500.0
4R1308TFR5002017-06-19 20:33:04.856000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1308...Left Middle Temporal Gyrus200.0500.0
\n", "
" ], "text/plain": [ " subject experiment montage session date \\\n", "0 R1275D FR5 0 1 2017-06-10 17:15:09.603000+00:00 \n", "1 R1275D FR5 0 2 2017-06-11 17:29:07.115000+00:00 \n", "2 R1292E FR5 0 0 2017-07-03 16:12:42.221000+00:00 \n", "3 R1304N FR5 0 0 2017-05-29 15:35:00.879000+00:00 \n", "4 R1308T FR5 0 0 2017-06-19 20:33:04.856000+00:00 \n", "\n", " session_summary electrode_type \\\n", "0 \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentmontagesessiondatesession_summaryelectrode_typelocationpulse_freqstim_durationdelta_recall
0R1275DFR5012017-06-10 17:15:09.603000+00:00<FRStimSessionSummary(_events=[ ( 1, 1, 'R1275...Right Cerebral White Matter200.0500.035.411540
1R1275DFR5022017-06-11 17:29:07.115000+00:00<FRStimSessionSummary(_events=[ ( 1, 2, 'R1275...Right Cerebral White Matter200.0500.0-7.532247
2R1292EFR5002017-07-03 16:12:42.221000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1292...Left Middle Temporal Gyrus100.0500.05.376644
3R1304NFR5002017-05-29 15:35:00.879000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1304...Right Cerebral White Matter200.0500.055.067482
4R1308TFR5002017-06-19 20:33:04.856000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1308...Left Middle Temporal Gyrus200.0500.0NaN
\n", "" ], "text/plain": [ " subject experiment montage session date \\\n", "0 R1275D FR5 0 1 2017-06-10 17:15:09.603000+00:00 \n", "1 R1275D FR5 0 2 2017-06-11 17:29:07.115000+00:00 \n", "2 R1292E FR5 0 0 2017-07-03 16:12:42.221000+00:00 \n", "3 R1304N FR5 0 0 2017-05-29 15:35:00.879000+00:00 \n", "4 R1308T FR5 0 0 2017-06-19 20:33:04.856000+00:00 \n", "\n", " session_summary electrode_type \\\n", "0 \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
delta_recall
meanminmax
experiment
FR53.739698-60.29362787.689632
catFR54.622367-39.46265661.873683
\n", "" ], "text/plain": [ " delta_recall \n", " mean min max\n", "experiment \n", "FR5 3.739698 -60.293627 87.689632\n", "catFR5 4.622367 -39.462656 61.873683" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "session_df.groupby(by='experiment').agg({'delta_recall': ['mean', 'min', 'max']})" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
delta_recall
meanminmax
location
---12.920880-23.027612-0.691512
Left Cerebral White Matter-28.169914-28.169914-28.169914
Left Middle Temporal Gyrus10.327724-36.40040487.689632
Left caudalmiddlefrontal-26.389477-35.246487-17.532468
Left inferiortemporal29.32419229.32419229.324192
Left middletemporal6.031860-39.46265661.873683
Left precentral6.6398696.6398696.639869
Left rostralmiddlefrontal1.6685921.6685921.668592
Right Cerebral White Matter5.663287-60.29362755.067482
Right MCg2.332751-4.1323048.797806
Right Middle Temporal Gyrus16.338604-45.99394778.671155
Right caudalmiddlefrontal11.73222911.73222911.732229
Right middletemporal-3.304355-28.24295921.634249
Right parstriangularis12.677112-3.59195428.946178
Right rostralmiddlefrontal13.4619965.02173222.924525
Right superiorfrontal-4.072398-4.072398-4.072398
Right superiortemporal-8.582085-12.624835-1.666667
\n", "
" ], "text/plain": [ " delta_recall \n", " mean min max\n", "location \n", "-- -12.920880 -23.027612 -0.691512\n", "Left Cerebral White Matter -28.169914 -28.169914 -28.169914\n", "Left Middle Temporal Gyrus 10.327724 -36.400404 87.689632\n", "Left caudalmiddlefrontal -26.389477 -35.246487 -17.532468\n", "Left inferiortemporal 29.324192 29.324192 29.324192\n", "Left middletemporal 6.031860 -39.462656 61.873683\n", "Left precentral 6.639869 6.639869 6.639869\n", "Left rostralmiddlefrontal 1.668592 1.668592 1.668592\n", "Right Cerebral White Matter 5.663287 -60.293627 55.067482\n", "Right MCg 2.332751 -4.132304 8.797806\n", "Right Middle Temporal Gyrus 16.338604 -45.993947 78.671155\n", "Right caudalmiddlefrontal 11.732229 11.732229 11.732229\n", "Right middletemporal -3.304355 -28.242959 21.634249\n", "Right parstriangularis 12.677112 -3.591954 28.946178\n", "Right rostralmiddlefrontal 13.461996 5.021732 22.924525\n", "Right superiorfrontal -4.072398 -4.072398 -4.072398\n", "Right superiortemporal -8.582085 -12.624835 -1.666667" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "session_df.groupby(by='location').agg({'delta_recall': ['mean', 'min', 'max']})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We just answered the two questions we posed previously in one line of code of each. Try doing that with a recarray all you doubters out there." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Combining Across Summary File Types\n", "\n", "So far, we have been taking a look at the data contained in each of these summary objects in isolation. However, for many analyses, you may want to combine the information contained in each. For example, the most commone usage pattern is to combine information about the trained classifier, the classifier results from the session, and the session summary data. The three of those together open up a world of possibilities for your analyses. To do this, we are going to generalize the \"Pandas magic\" approach by creating a dataframe containing the summary object of interest and then rely on the merging mechanics available through pandas to get all of this data into a single data structure." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "from ramutils.reports.summary import ClassifierSummary\n", "\n", "# Repeat the process for getting summary objects into a dataframe that we used for the session summaries\n", "classifier_summaries = glob.glob(os.path.join(report_db_location, '*FR5*classifier_session*'))\n", "\n", "classifier_summary_data = {\n", " 'subject': [],\n", " 'experiment': [],\n", " 'session': [],\n", " 'classifier_summary': [],\n", " 'montage': [],\n", " 'auc': [],\n", " 'pvalue': []\n", "}\n", "\n", "for summary_dir in classifier_summaries:\n", " summary = ClassifierSummary.from_hdf(summary_dir)\n", " metadata = extract_report_info_from_path(summary_dir)\n", " classifier_summary_data['subject'].append(metadata['subject'])\n", " classifier_summary_data['experiment'].append(metadata['experiment'])\n", " classifier_summary_data['session'].append(metadata['sessions'][0])\n", " classifier_summary_data['classifier_summary'].append(summary)\n", " classifier_summary_data['montage'].append(metadata['montage'])\n", " classifier_summary_data['auc'].append(summary.auc)\n", " classifier_summary_data['pvalue'].append(summary.pvalue)\n", " \n", "classifier_df = pd.DataFrame.from_dict(classifier_summary_data)\n", "classifier_df = classifier_df[['subject', 'experiment', 'montage', 'session', \n", " 'classifier_summary', 'pvalue', 'auc']]" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentmontagesessionclassifier_summarypvalueauc
0R1275DFR501<ClassifierSummary(_predicted_probabilities=[ ...0.0150.616795
1R1275DFR502<ClassifierSummary(_predicted_probabilities=[ ...0.4100.515080
2R1292EFR500<ClassifierSummary(_predicted_probabilities=[ ...0.3900.515818
3R1304NFR500<ClassifierSummary(_predicted_probabilities=[ ...0.3150.526984
4R1308TFR500<ClassifierSummary(_predicted_probabilities=[ ...0.0000.705267
\n", "
" ], "text/plain": [ " subject experiment montage session \\\n", "0 R1275D FR5 0 1 \n", "1 R1275D FR5 0 2 \n", "2 R1292E FR5 0 0 \n", "3 R1304N FR5 0 0 \n", "4 R1308T FR5 0 0 \n", "\n", " classifier_summary pvalue auc \n", "0 \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentmontagesessiondatesession_summaryelectrode_typelocationpulse_freqstim_durationdelta_recallclassifier_summarypvalueauc
0R1275DFR5012017-06-10 17:15:09.603000+00:00<FRStimSessionSummary(_events=[ ( 1, 1, 'R1275...Right Cerebral White Matter200.0500.035.411540<ClassifierSummary(_predicted_probabilities=[ ...0.0150.616795
1R1275DFR5022017-06-11 17:29:07.115000+00:00<FRStimSessionSummary(_events=[ ( 1, 2, 'R1275...Right Cerebral White Matter200.0500.0-7.532247<ClassifierSummary(_predicted_probabilities=[ ...0.4100.515080
2R1292EFR5002017-07-03 16:12:42.221000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1292...Left Middle Temporal Gyrus100.0500.05.376644<ClassifierSummary(_predicted_probabilities=[ ...0.3900.515818
3R1304NFR5002017-05-29 15:35:00.879000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1304...Right Cerebral White Matter200.0500.055.067482<ClassifierSummary(_predicted_probabilities=[ ...0.3150.526984
4R1308TFR5002017-06-19 20:33:04.856000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1308...Left Middle Temporal Gyrus200.0500.0NaN<ClassifierSummary(_predicted_probabilities=[ ...0.0000.705267
\n", "" ], "text/plain": [ " subject experiment montage session date \\\n", "0 R1275D FR5 0 1 2017-06-10 17:15:09.603000+00:00 \n", "1 R1275D FR5 0 2 2017-06-11 17:29:07.115000+00:00 \n", "2 R1292E FR5 0 0 2017-07-03 16:12:42.221000+00:00 \n", "3 R1304N FR5 0 0 2017-05-29 15:35:00.879000+00:00 \n", "4 R1308T FR5 0 0 2017-06-19 20:33:04.856000+00:00 \n", "\n", " session_summary electrode_type \\\n", "0 \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentmontagesessiondatesession_summaryelectrode_typelocationpulse_freqstim_durationdelta_recallclassifier_summarypvalueauc
\n", "" ], "text/plain": [ "Empty DataFrame\n", "Columns: [subject, experiment, montage, session, date, session_summary, electrode_type, location, pulse_freq, stim_duration, delta_recall, classifier_summary, pvalue, auc]\n", "Index: []" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check if we are missing any classifier summaries\n", "combined_df[combined_df['classifier_summary'].isnull() == True]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the classifier information joined with the session information, we can now ask a wider set of question. For example, how does classifier performance correlated with delta recall?" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/zduey/miniconda3/envs/ramutils_dev/envs/ramutils_test/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n", " warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAGoCAYAAAAuIBCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl83VWd+P/XuXv2fWvTLTTd0xboAt0oSCltERRRQVQQUVQGnR+PceQxivMVdVxm3GEEBEdQWRSXolTEAqWydgG6l26kTdLs+3L3z/n9cZOQpFnuTe6e9/PxyOMmN3c5N03vO+ec9/t9lNYaIYQQItGZYj0AIYQQIhwkoAkhhEgKEtCEEEIkBQloQgghkoIENCGEEElBApoQQoikIAFNCCFEUpCAJoQQIilIQBNCCJEULLEeQBRIKxQhRKJTsR5AIpAZmhBCiKQwGWZoQsTcY2+cCcvjfGzl9LA8jhDJSGZoQgghkoIENCGEEElBApoQQoikIAFNCCFEUpCAJoQQIilIQBNCCJEUJKAJIYRIChLQhBBCJAUJaEIIIZKCBDQhhBBJQVpfCRFHnB4/Z9udtPV48PgDfbWzHFZy02wUZtpjPDoh4psENCFirMPpZV91G29XtVHb7hrxdg6riZ3HGtlUUczGhcWk2uS/rxADKa2T/nSVpH+BIv4N15y4rcfDjnca2XO6BUPDtJwU5pVkUpqdQn66HavFhKE1HU4vjZ1u3m3qprbdRU2bkzSbmQ+cP5VPr5lFWUF6DF6RiDI5PiYIEtCEiIKBAc3nN3jpWCMvHWtEa1g2M4fVs/PJTx97SfH65dPYXdnC7/dW8/TbZ/EaBpsWFfP/XT6H8qKMSL4EEVsS0IIgAU2IKOgLaGdaevjD3moau9wsLs1i48JiclJtQT/OwONjGjpdPPrqaX71aiU9Hh8fPL+Uf728nGm5qWEfv4g5CWhBkIAmRBT85vXT7HingReONpCZYuUDS6cyZxwzquHOQ2vp9vDzHSd45LXTaK25YcV0vvS+cvKCmPGJhCEBLQgS0ISIsIYOFx998HXebepm6bRsrl4yBYfVPK7HGu2Az9p2Jz974QRP7q4i1WbmS+8r55MXz8RmkeqcJCABLQgS0ISIoNdONnPH42/R7vRwzdKpXDA9Z0KPF8yJ1cfrO/nmM0fYeayRsvw0vnbVfC6dW4hS8p6YwOQfLwgS0ISIAMPQPLDzFP/996PMzE/jqsVTKM50TPhxgwloAFprXnyngW/99QinmrpZN6eAu7fMl8SRxCUBLQgS0IQIs9ZuD19+ah/bjzSwZXEJ3/vQYp5++2xYHjvYgNbH4zN49LVKfvL8cXo8fm66eCZ3XjGHdLvUsCUYCWhBkIAmRBjtrmzhi4+/RVOXm69uns9Nq2ailBq2Dm08Qg1ofZq73PzgH8d4fNcZijIc/L+rF7BxYbEsQyYO+YcKguwWCxEGfkPzs+eP89EHXsNuMfHHz6/m5tWz4iZg5KXb+a8PVvDHz68iJ83G537zJrc+soeqlp5YD02IsJEZmhATVNfu4s7fvc2rJ5u5ZukUvvWBRWQ4rINuE64ZWjj4Dc1rJ5vYfqQBjWbTohJWzsodd/Ad76xRhCQ+/jKKc7KQLsQ4GYbmid1VfGfbEbyGwfc/tJgPLyuNm1nZSMwmxZryAhZNzeJPb9Xw9L6zHKnt4NoLSslKsY79AELEKZmhCTEOlU3d3PXH/bx+qoWLy/L47ocqmJGXNuLt42mGNpDWml2VLWw7UIvZpLh6yRSWlGaHFJRlhhYV8f1XUpyQGZoQIXD7/Pzy5Up+8vwxrCYT3722go8unxb3s7KRKKVYOSuP2QXp/H5vNb/bU83Ruk4+sHTquIu/hYgVCWhCBEFrzbYDdXz32SNUtTjZsKCIb16ziOKsideWxYO8dDufXVfGzmONbD9ST1VLD9cvny59IUVCkYAmxBjermrjW389zJ7TrcwrzuDXn17B2vKCWA8r7ExKsX5uIWUF6Ty5+wwP7DzJhvlFrJ1TgClBZ6BicpE9NCFGcLbNyRd++yZvV7WRZrdwxfwiLpyZMyne3J0eP39+u4YDNe2cV5DGh5dNI9MxfMKI7KFFRfL/0oWBBDQhhuh2+7j/pZM8uPMUfkOzenY+l8wpmHR7Slpr9p5u5S/7z2I1m/jwhaXMLc4853YS0KJCAloQJKAJ0cswNE+9Wc3//P0dGjrdXLW4hPnFmeSkBX9eWTJq6HDxxO4q6jpcrD4vj40Li7GY3+vJIAEtKiSgBUECmhDA66ea+eZfD3PobAdLp2Vz91ULuHBGTtym20eb12/w7ME6XjvVzJQsB9cvn05+RuC8tXAFtFi3B4tzEtCCIEkhYlI73dzNf207wt8P1TMly8FPrl/K+xdPwWSS94+BrGYT718yhdmF6fzhzWruffEE718yhQumZ8d6aEL0k4AmJiWPz+CBl07ysxdPYDEp/u2KOdy6tmzS7ZOFan5JJndcVs7v9lTxhzerOXS2nXVzCiS9X8QFCWhi0tld2cJ//PEAxxu62LK4hK9ftYCiMJxVNllkpVj59JpZvHKiieePNHD5D1/iM2vL+My6MmmdJWJK9tBEUhpuP8bp8fPsoTp2V7aQnWLl6qVTmDdM1p4IXrvTy+HaDv6y7yyZDgs3r57Fx1ZMH1fBueyhjUrWwIMgMzSR9LTWHKhp56/7a+l2+1gzO5/3zS/EbpHlxYnKSrHysxvO53OXlPHj7cf52QvHue/FE1w6t5ArFhZx6dxCCnqTR4SINAloIqm1dnvYuq+GY/VdTM1O4aZVM5manRLrYSWdhVOy+MUnl3GmuYff7jrN1rfOsv1IPQClOSksnJLJtJxUirMcpNosWEwKs0lhMStMSmFozVtnWgPLKRosZoXNbCLNbiHDYSEzxTopCtrFxMiSo0hKv37tNK+ebGL7kXoUissXFHFxWR5myV4Mu+GW+LTWHK7t4J/HmzhY087h2g7OtjlxeY1xPYfFpCjMsFOSncKs/DTK8tPITh2+PlCWHCcvmaGJpLOnsoV7XzxOfYebecUZXL1kyohvfiIylFIsnJLFwilZ/ddprelw+XB5/Xj9Bn5D4zM0hqExmRTb9tf2n1rg9Rt4fAbdHh8dTh/NXW7qOlwcPtvB3tOtAEzNDsz8zp+eI8koApCAJpJIa7eH7/7tKE/uqSIrxcrHV05nfklmwh7tkmyUUmSlWEcMPnnpY++1GVpT3+HieH0Xh86289zherYfqWdecSYrZuUyuzA93MMWCUQCmkh4Lq+fR1+r5H93nKTL5eO2dWUUZzkk6SMJmZSiJCuFkqwU1s0poKXbw653W9h7uoXDtR3kptlwevxcv2IaGSM0UxbJS/bQRMLy+Q2e2lvNj7cfp67DxSVzCviPzfOZW5whLauiKB5aX/n8BodqO3jjVDOVzT1k2C3csHI6N6+ayZTkSAKSZYYgyAxNJBy3z8/Tb5/l5y+d5FRjN+dPz+ZHH13KxeflxXpoIkYsZhNLSrNZUppNxdQsfvHPUzz88rv88uV3uWpxCbeuLWPR1KyxH0gkNJmhiYRR2dTNH96s5vFdVTR1uZlblMGdV8zhigVF5+yTyQxt8uqbMVa39vB/r1TyxK4zdHv8rDovj8+sK2P9nIJE3FdNuAHHggQ0Ebe01hxv6GL7kXqeO1TP21VtmBRcMqeAW9bMYs3s/BHfmCSgTV5Dl0DbnV6e2HWG/3ulkroOF3OK0rlhxXSuXjIlqESUOCEBLQgS0ERc8PoNzrY5ebepm0NnOzhQ3c7uyhaauz0AVEzNYlNFMdeeXxpUWyUJaJPXSHt6Hp/BMwfO8suXKzlQ047FpLh0XiEfuqCUS+cVxHsSkQS0IEhAS1Dx1vduuPEYWuP2GnS7fXR7fHS7/XR7fPS4fXR7/HS7fbQ7vbT2eGjr8Q76h5qem8qymTlcNCuPtXPyKckKbWNfAtrkFczv9NG6Dv74Zg1/equGxk43aTYzF5+Xx9ryAtaW5zMrPy3eliXjajDxSgJaDBiGptPlo80ZeCNvd3ppcwYuO3ov23uvb3d66XL7AoWmfgOv38Dn1zg9fkymQNsgsylwXpXdYsJuMeOwmrBbzTh6v7ZbTTiGfm0x47CauWnVjEGnD480XqfXT4fLS4fT13vppcMVGGdLj5c3TjX3B6me3uDV4/FhjPDTt5pVoK2R3UJeup2cVCu5aTY+smwa80oyJ1woKwFt8grljzSf3+CfJ5p4/kg9O481caalB4ApWQ4WTc1iXkkmC0oymFecSWlOypj/VyJIAloQJKANo6HTRUOHG7+h8etAJ4O+z/29nxta4zfo7XZg4PT4cXr9OD1+ejx+XN7AZafLS2tPb8Dq8fQHrtF+7A6rqb8ANSvFSrrdgs1iwmo2YTObMJsUJxu7MDT9Y/H6DVxeA7fP33/p9hpBvXi7xYTFpDD19tczK4XP0Hh8vQF0pKg0QIrVTJrdTJrNQprdQqrNTJo98Hla3+c2C2l2M6m2wOsZTjykgIvENpHfodPN3ew83sTrp5o5UttBZVN3/x9lJgUFGXaKMx0UZTooyLCT4bCS4Qj0m0y3W7BbzP19KK1mExazClz2/vGpVOBMuXG0YJOAFgRJ2x/GY2+c4cfbj0/oMewWE6k2M5kpVrJTrGSl2piRm0p26ntfZ/cGrOzU94JXZoo1qEMmg3nDNrTG6zNw+QzcXj9un4GrN9D1Bb45RRn0eHz4BgVqjcWk3guilsAMLyvFSqbDSmaKpfcy8J85O8XK7/ZUT+jnJUQ8mJGXxify0vjERTOAwJFDxxs6OVrbSXWbk7p2J7XtLiqbu9lzupVOlxevP7S/mQ99YyNpdnnrjYSkn6EppZ4F8mPw1PlAUwyedzQypuDImIIjYwpOOMbUpLW+MhyDSWZJH9BiRSm1R2u9LNbjGEjGFBwZU3BkTMGJxzElq5jtcAohhBDhJAFNCCFEUpCAFjkPxnoAw5AxBUfGFBwZU3DicUxJSfbQhBBCJAWZoQkhhEgKEtCEEEIkBQloQgghkoIENCGEEEkh6QPalVdeqQn0c5QP+ZAP+UjUj6Al4Xte0JI+oDU1xVsXHCGEiJzJ/J6X9AFNCCHE5CABTQghRFKQgCaEECIpyKE8Iil4vV6qq6txuVyxHooQ4+ZwOCgtLcVqndiJ7ZOVBDSRFKqrq8nIyGDmzJkoJYf7isSjtaa5uZnq6mpmzZoV6+EkJFlyFEnB5XKRl5cnwUwkLKUUeXl5ssowARLQRNKQYCYSnfwOT4wENCGEEElBApoQImTPPvssc+fOZfbs2Xz3u98d9jb3338/FRUVLF26lDVr1nD48GEAPB4Pn/rUp6ioqGDJkiXs2LEjiiN/T0tLCxs2bKC8vJwNGzbQ2to67O0eeeQRysvLKS8v55FHHum/fu/evVRUVDB79my++MUv0ncU1913383ixYtZunQpV1xxBWfPngXgt7/9LYsXL2bx4sWsWrWKffv2Rf5FTjZa66T+uPDCC7VIfocPH471ECLC6/VG5HF9Pt+E7ltWVqZPnjyp3W63Xrx4sT506NA5t2tvb+//fOvWrXrjxo1aa63vvfdeffPNN2utta6vr9cXXHCB9vv94x7PeH35y1/W3/nOd7TWWn/nO9/R//7v/37ObZqbm/WsWbN0c3Ozbmlp0bNmzdItLS1aa62XL1+uX331VW0Yhr7yyiv1tm3btNaDX/dPfvITfdttt2mttX7llVf677tt2za9YsWKYcc1wu/yZH7PC/q1ywxNiDCorKxk3rx53HTTTSxevJjrrruOnp4eIPCX/CWXXMKFF17Ixo0bqa2tBeAXv/gFy5cvZ8mSJXzoQx/qv/3NN9/MnXfeyaWXXspXvvIVXnrpJZYuXcrSpUs5//zz6ezsRGvNl7/8ZRYtWkRFRQVPPvkkADt27GD9+vVcd911zJs3jxtvvLF/5jBz5kzuuece1qxZw+9///txv9Zdu3Yxe/ZsysrKsNlsXH/99WzduvWc22VmZvZ/3t3d3b8/dPjwYd73vvcBUFhYSHZ2Nnv27AHg1ltv7f98oJtvvpnPfe5zrF27ljlz5vDXv/513OPvs3XrVm666SYAbrrpJv785z+fc5u///3vbNiwgdzcXHJyctiwYQPPPvsstbW1dHR0cPHFF6OU4pOf/GT//Ud63atWrSInJweAiy66iOrq6gm/BjGYpO0LESbvvPMODz/8MKtXr+aWW27hf//3f/nSl77EHXfcwdatWykoKODJJ5/kq1/9Kr/85S+59tpr+cxnPgPA1772NR5++GHuuOMOAI4dO8b27dsxm828//3v57777mP16tV0dXXhcDj44x//yNtvv82+fftoampi+fLlrFu3DoC33nqLQ4cOMWXKFFavXs0rr7zCmjVrgECd08svv3zO2H/729/y3//93+dcP3v2bJ566qlB19XU1DBt2rT+r0tLS3njjTeG/Zncd999/PCHP8Tj8fDCCy8AsGTJErZu3cr1119PVVUVe/fupaqqihUrVvDQQw+N+POtrKzkpZde4uTJk1x66aWcOHECh8PR//3Ozk7Wrl077H0fe+wxFixYMOi6+vp6SkpKACgpKaGhoeGc+w33WmtqaqipqaG0tPSc6/t89atf5dFHHyUrK4sXX3zxnMd9+OGH2bRp04ivdSI8PiMij5sIJKAJESbTpk1j9erVAHz84x/npz/9KVdeeSUHDx5kw4YNAPj9/v430YMHD/K1r32NtrY2urq62LhxY/9jffjDH8ZsNgOwevVq7rzzTm688UauvfZaSktLefnll7nhhhswm80UFRVxySWXsHv3bjIzM1mxYkX/m+3SpUuprKzsD2gf/ehHhx37jTfeyI033hjU6+yb8Q00Unbe7bffzu23385jjz3Gt771LR555BFuueUWjhw5wrJly5gxYwarVq3CYhn7regjH/kIJpOJ8vJyysrKOHr0KEuXLu3/fkZGBm+//XZQryFYI73WsX4G3/72t/n2t7/Nd77zHe69916+8Y1v9H/vxRdf5OGHHx72D4twaOnxRORxE4EENCHCZOibet8b38KFC3nttdfOuf3NN9/Mn//8Z5YsWcKvfvWrQckRaWlp/Z/fddddbNmyhW3btnHRRRexffv2Yd9Q+9jt9v7PzWYzPp9v2McdKJQZWmlpKVVVVf1fV1dXM2XKlBHHA3D99dfz+c9/HgCLxcKPfvSj/u+tWrWK8vLyUe8Pw/98Bwp1hlZUVERtbS0lJSXU1tZSWFh4zv1KS0sH/btUV1ezfv16SktLBy0ZjvQz+NjHPsaWLVv6A9r+/fu59dZb+dvf/kZeXt7oL3i8QjpwJbnIHpoQYXLmzJn+wPX444+zZs0a5s6dS2NjY//1Xq+XQ4cOAYE34JKSErxeL7/97W9HfNyTJ09SUVHBV77yFZYtW8bRo0dZt24dTz75JH6/n8bGRnbu3MmKFSvGPfYbb7yRt99++5yPocEMYPny5Rw/fpx3330Xj8fDE088wdVXX33O7Y4fP97/+TPPPNMftHp6euju7gbgH//4BxaLpT/YfPKTn2TXrl3DjvH3v/89hmFw8uRJTp06xdy5cwd9v2+GNtzH0GAGcPXVV/dnLT7yyCNcc80159xm48aNPPfcc7S2ttLa2spzzz3Hxo0bKSkpISMjg9dffx2tNY8++mj//Qe+7qeffpp58+YBgd+Pa6+9ll//+tfMmTNn2NcoJkZmaEKEyfz583nkkUe47bbbKC8v5/Of/zw2m42nnnqKL37xi7S3t+Pz+fjXf/1XFi5cyDe/+U1WrlzJjBkzqKiooLOzc9jH/fGPf8yLL76I2WxmwYIFbNq0CZvNxmuvvcaSJUtQSvH973+f4uJijh49GvHXabFYuPfee9m4cSN+v59bbrmFhQsXAvD1r3+dZcuWcfXVV3Pvvfeyfft2rFYrOTk5/cGjoaGBjRs3YjKZmDp1Kr/+9a/7H3v//v39S7JDzZ07l0suuYT6+nruv//+Qftn43HXXXfxkY98hIcffpjp06f3J8rs2bOH+++/n4ceeojc3Fzuvvtuli9f3v/6cnNzAfj5z3/OzTffjNPpZNOmTf17YnfddRfvvPMOJpOJGTNmcP/99wNwzz330NzczBe+8IX+n+NwCTATNYknaKjRli6SwbJly3QkfmlEfDly5Ajz58+P2fNXVlZy1VVXcfDgwZiNIdF1dHTw6U9/etgMzJtvvpmrrrqK6667LgYji64RfpeDbiEytXyRrjmeVL+HQb92WXIUQsSFzMzMCZUTCCFLjkKEwcyZM2V2FkG/+tWvYj0EkQBkhiaSRrIvn4vkF47f4cn8/0ACmkgKDoeD5ubmSf2fWSQ23Xse2kSTXSbz/wBZchRJoa8uqLGxMdZDEWLc+k6snpBJHNEkoImkYLVa5ZRfIZjU8UyWHIUQIplM5mV3CWhCCJFEJm84k4AmhBBJZRJP0CSgCSFEMtGTeI4mAU0IIZKIzNCEEEIkBQloQgghkoIxiSOaBDQhhEgiEtCEEEIkBWPyxjMJaEIIkUxkhiaEECIpSEATQgiRFLQG/yRdd5SAJoQQSabL5Yv1EGJCApoQQiSZ1h5PrIcQExLQhBAiybQ5vbEeQkxIQBNCiCTTJjM0IYQQyaCtR2ZoQgghkoDM0IQQQiSFFpmhCSGESHQWk6KhwxXrYcREzAOaUuqXSqkGpdTBAdflKqX+oZQ63nuZ03u9Ukr9VCl1Qim1Xyl1QexGLoQQ8cdqNlHbLgEtVn4FXDnkuruA57XW5cDzvV8DbALKez8+C/w8SmMUQoiEYDWbqJcZWmxorXcCLUOuvgZ4pPfzR4APDLj+UR3wOpCtlCqJzkiFECL+WcyKOglocaVIa10L0HtZ2Hv9VKBqwO2qe68bRCn1WaXUHqXUnsbGxogPVgghYmnge57X5aStx4vL64/1sKIuXgPaSNQw153ThVNr/aDWepnWellBQUEUhiWEELEz8D0vKzMdYFIuO8ZrQKvvW0rsvWzovb4amDbgdqXA2SiPTQgh4pbVHHhbr5uEiSHxGtCeBm7q/fwmYOuA6z/Zm+14EdDetzQphBACrKbegDYJZ2iWWA9AKfU4sB7IV0pVA/8JfBf4nVLq08AZ4MO9N98GbAZOAD3Ap6I+YCGEiGMWS2BnZjLO0GIe0LTWN4zwrfcNc1sN3B7ZEQkhROIyK0WqzTwpZ2jxuuQohBBinIozHZIUIoQQIvEVZTom5ZKjBDQhhEgyxVkO6jvcsR5G1ElAE0KIJFOYaaeh00Ug7WDykIAmhBBJpiDdjtevaXdOrmNkJKAJIUSSKciwA9DUNbmWHSWgCSFEkukLaA2dkyugxbwOTQgRn3YcbeCBnaeoau1hWk4qt60rY/28wrHvKGKusDegNU6ygCYzNCHEOXYcbeDrTx+iodNFdoqVhk4XX3/6EDuONox9ZxFzBekOQAKaEELwwM5TWM2KVJsFpQKXVrPigZ2nYj00EYTMFAs2s4lG2UMTQkx2Va09pFjNg65LsZqpbu2J0YhEsFq6PTy+q4oUm5ldp1p47I0zsR5S1EhAE0KcY1pOKs4hB0Q6vX5Kc1JjNCIRqgyHhS63L9bDiCoJaEKIc9y2rgyvX9Pj8aF14NLr19y2rizWQxNBSrdb6HRJQBNCTHLr5xVyz9ULKcxw0O70Upjh4J6rF0qWYwLJcFjonGQzNEnbF0IMa/28QglgCSzdbqXH7cNvTJ72VzJDE0KIJJThsKCBbs/kmaVJQBNCiCSUbg8swHVNon00CWhCCJGEMhy9AW0S7aPJHpoQIu5JG67QZTisAJMq01FmaEKIuCZtuMbnvSXHyXOEjAQ0IURckzZc42OzmLBZTJMqdV8CmhAirkkbrvHLmGTF1RLQhBBxTdpwjV/6JGt/JQFNCBHXpA3X+MkMTQgh4oi04Rq/dIeVLvfkSQqRtH0hRNyTNlzjk+mw4PIa9Hh8pNqS/+1eZmhCCJGkMlMCtWh17a4YjyQ6JKAJIUSSypKAJoQQIhn0BbSzEtCEEELEktvnp6bNOe77vzdDG/9jJJLk3yUUQogEo7WmrcdLm9OL1uM/z8xqNpFqM1M7SWZoEtCEECKOuLx+GjvdeP1GWB4vK8XK2QnM8hKJBDQhhIgDWmtauj20O8NbN5abZuN08+RoEyZ7aEIIEWNOj5/qVmfYgxlAQbqd0y09eHzhmfHFMwloQggRI4ahaex0U9vuDNsS41AFGXb8huZMS/LP0iSgCSFEDPR4fFS3OumM8HllBRl2AE42dkX0eeKB7KEJIUQU+Q1Nc7ebrig1Dc5Pl4AmhBAizLrcPpq73PiN8afih8phNVOYYedkQ3fUnjNWJKAJIUSE+fwGTV0eejyxOcrlvIJ0TjR0xuS5o0kCmhBJbsfRBh7YeYqq1h6m5aRy27oy6VwfRe1OL63dHowJFEhP1MIpmTz6+mm8fgOrOXlTJ5L3lQkh2HG0ga8/fYiGThfZKVYaOl18/elD7DjaEOuhJT2Pz+Bsm5PmLndMgxlARWkWHp/B8frk3keTgCZEEntg5ymsZkWqzYJSgUurWfHAzlOxHlrSCrSt8lDT5sTl9cd6OABUTM0C4GBNe4xHElkS0IRIYlWtPaRYzYOuS7GaqW5N/pqkWHB5A82EW7o9E+rBGG4z89JIt1s4IAFNCJGopuWk4hwyS3B6/ZTmpMZoRMnJMDTNXW7OtjnjsiOHyaRYOCVTApoQInHdtq4Mr1/T4/GhdeDS69fctq4s1kNLGn0F0pFoWxVOFVOzOFLbEbGOJPFAApoQSWz9vELuuXohhRkO2p1eCjMc3HP1QslyDAOf36C+w0VduwufEf9BYvG0bNw+g6O1yZu+L2n7QiS59fMKJYCFWTyk4odq+cwcAHZXtlBRmhXj0USGBDQhQiR1XZOXy+unuduDO06yF0NRkpXC1OwU9pxu4ZY1s2I9nIiQJUchQiB1XZOTYWiaepM+EjGY9VkxK5fdla1xlYEZThLQhAiB1HVNPl3uQNJHR5wnfQRj2cwcGjvdSXvgpwQ0IUIgdV2Th9dvUNfuoqEjMZI+grF8Zi4Q2EdLRhLQhAiB1HUlv75OH9UGr3dGAAAgAElEQVStzpg1E46U2QXpZKdaJaAJIaSuK9k5PX6qW+Ov00e4mEyKlbNyeeVEc1K+PslyFCIE6+cVcg+BvbTq1h5KJcsRSPzMz2gfuhlLa8sL+Puhek41dXNeQXqshxNWcR3QlFKVQCfgB3xa62VKqVzgSWAmUAl8RGvdGqsxislH6roG68v8tJrVoMzPeyAhfk4drkBNWTQP3YylS+YUALDzWGPSBbREWHK8VGu9VGu9rPfru4DntdblwPO9XwuRNHYcbeCGB19nzfde4IYHX4/7koBEzfx0ef1Ut/bQ1BndE6RjbVpuKrPy09h5rDHWQwm7uJ6hjeAaYH3v548AO4CvxGowQoRTIs52qlp7yE6xDrounjM/vX6D1m4PXe7kX17s89gbZwZ9XZRp5+UTTbh9fuwW8wj3SjzxPkPTwHNKqb1Kqc/2Xlekta4F6L0853+5UuqzSqk9Sqk9jY3J91eISF6JONtJlMxPf29H/OpWZ9IFs4HveZ1tY2cwlhdm4PVr9lYm125NvAe01VrrC4BNwO1KqXXB3Elr/aDWepnWellBQUFkRyhEGCVinVu8Z372peFXtfTQ7vQmZXbfwPe8jOzcMW9flp+GWSl2JNmyY1wHNK312d7LBuBPwAqgXilVAtB7Gd8bDEKEIFFmOwPFc0f/LrePqpZAGn4iNRKONLvVzHmFaTx7sC6pAnzc7qEppdIAk9a6s/fzK4B7gKeBm4Dv9l5ujd0ohQiv29aV8fWnD9Hj8ZFiNeP0+uNqtjOSeMv89PkNmrs9dCfZ0mI4LZySxZ/equFwbQcLpyRH9/14nqEVAS8rpfYBu4BntNbPEghkG5RSx4ENvV8LkRTiebaTKDpdXmranBLMxjC/JBOTgmcP1sV6KGETtzM0rfUpYMkw1zcD74v+iISIjnib7SQKj8+gqcuNK4G74UdTut3Cyll5bDtQy50b5qCUivWQJiyeZ2hCCDEmw9C0dHuoaXMmVTDTWnOktoMfPHcsYs9x1ZISTjZ2c7CmI2LPEU1xO0MTQoixdLq8tHZ7k6YbPkCH08v2I/VsO1DHqaZuAO6L0HNdtXgK3/jLYZ7aW5UUp1hLQBNCxMx4e0Am8snRwzG0Zl9VG9sO1LHzeCNe/3uZh2X5aRF73qwUKxsXFrN131n+Y8v8hC+yloAmhIiJ8XRF8fkNWno8SdNEuLnLzd8P1bPtYC1n21z916dYzVw2r5Ati4uZW5QR0TF8+MJS/rLvLM8faWBzRUlEnyvSJKAJIWJiYFcUgFSbhR6Pjwd2njonoGmtaXd6aevxJnw9md/Q7Hq3hWcO1PL6qWYGtpFcUJLJlsUlrJ9TQIotOrOl1bPzKcly8PiuMxLQhBBiPILtAdnt9tHS7cHrT+x9srNtTv52sI5nD9XR3OXpvz7TYeGKhUVsWlTCrAguL47EbFLcuHI6//PcMU40dDK7MLIzwkiSgCaEiIlpOak0dLr6Z2gwuCuK2+enpduD05O4+2Qen8ErJ5p45kAtb55pG/S9C2fksKWimFXn5WOzxDbh/IYV0/npCyf4v1cq+fYHK2I6lomQgCaEiImRuqJ8Zs0smrrcdDi9sR7iuL3b1M22A7X843A9HQP2+/LSbWxaVMymRcWUZKXEcISD5aXb+cDSKfzxzRq+vHEu2am2WA9pXCSgCSFiYrjTvz950QzOK0pPyGDm9PjZ8U4Dzxyo5XBtZ//1JgUXl+WxuaKEFbNyMZvis4D5U6tn8bs91Ty+q4rPrz8v1sMZFwloQiSg8aa7x5u+rigDu3wk0mGbWmveqe/kmf11vHC0YVBj6SnZDjYvKmHjwiLy0u0xHGVw5pdksmZ2Pg+//C6fWj0ThzXxUvgloAmRYBLxENCRaK1p7fEm3LEugeLnBrYdrOVUY3f/9VazYl15AVsWl7CkNCvh2kn9y2Wzuf7B13li1xluXj0r1sMJmQQ0IRJMKOnu8czp8dPU5U6Y7MWxip83V5Rw+fxCModkbiaSlbNyWT4zh/tfOsUNK6cnXKG1BDQhEkyw6e7xyt/be7HTlRj7ZKMVP79vfiGbFhUzrzgj4WZjw1FKccdl5Xzyl7v4w94aPrZyeqyHFJIxA5pSatTjT7XWY5/3LYQIm7HS3eNZl9tHc5c77vfJ+oqftx2o5bVhip83VxRz6dzCqBU/R9Pa8nyWTsvmvhdPcO0FUxNqLy2YGdpeQAPD/fmhgfg+eVCIJJOIh4B6fAbN3e64rymrbXey7cDwxc8bFhSxuSI2xc/RpJTi3zfO5WMPvcFvXj/NrWvj9/dqqDEDmtY68XYGhUhiw6W7x2uWo9aath4vbXGc9NFX/LztQC17hxQ/XzA9my0VJayeHfvi52haNTufteX53PviCT68bBpZCbIvGMyS4wWjfV9r/Wb4hiOECEYiHALq9vlp7HTj8cVn0keiFT9HwmNvnBnxe0tKs/nn8Sa++PhbbFxYPOh78bq3FsyS4w9G+Z4GLgvTWIRIaolaOxbquOM5Ff+94uc6Dte+d6hlX/HzlsUlLJ8Zv8XP0TQlO4UlpVm8erKJi8ryEmKWFsyS46XRGIgQyWCkN/9ErR0LddzxmIqfTMXP0bZhQTEHazp44Wg9Hzy/NNbDGVNIaftKqUXAAsDRd53W+tFwD0qIRDTam3+8146NFIiDHbff0DR3u+PqnLJOl5d/HB6++PmSOQVsWlTMkmnZmJIg3T5SctNsrCjL5fWTzaw+L5/CTMfYd4qhoAOaUuo/gfUEAto2YBPwMiABTQhGD1rxXDs2WiAOZtztTi+t3Z64OKdMa82+6na2HajlpWPDFT8Xc/n8ooQufo62S+cW8ubpVv52sI6bVs2M9XBGFcoM7TpgCfCW1vpTSqki4KHIDEuIxDPam3+wtWOx2GcbLRCPNm6nx09zd3wkffQVP//tYB01bc7+64ee/JwMxc/Rlm63cNm8Qv52sI536jqZWxy/56WFEtCcWmtDKeVTSmUCDUgNmhD9RnvzD6Z2LFb7bKMF4m9es+iccXt8Bh9dVkptu3OER4wOv6HZXRk4+fm1k0OLnzPYUlHC+iQtfo62i8/L6z9le3ZheqyHM6JQAtoepVQ28AsCxdZdwK6IjEqIBDRa0AqmdixW+2yjBeKB465q6aY4K4XrLixl8bTsiI1nLLXtvSc/H6yjKY5Ofk5mFpOJLRUlPPr6aV4/1cwnLp4R6yENK+iAprX+Qu+n9yulngUytdb7IzMsIRLPWEFrrNqxWO2zjTV7XD+vkAtn5tDa7cVnxGZ50eMzePVkE8/sP7f4+cLp2WyehMXP0Ta3OIPywnSeP1pPc5c7LrNCQ0kK+SDwgta6XWtdqZTKVkp9QGv95wiOT4iEEmrB88A9sw6nF5/foCDjvUyyaPRoHC0Qu31+mro8uL2xaVlV2Rwofn7u0LnFz1cuDBQ/T8lO7uLneKGUYnNFCT974Tg//Mcxvv3BilgP6RyhLDn+p9b6T31faK3bejMfJaAJEaSBASzDbqGxy01WipXsFCt+w6ChM7CElp9uj2qPxqGB2G9omrrcMTk52un1s+OdRrYdqOXQ2cHFzxeV5bElzk9+TmZFmQ5WluXx+K4zfPyiGcwvyYz1kAYJJaANN5eX42eECNLQpI8TDV34DE2azYKyKfLTAzOzbrefdqc3Zj0a251e2no8Ue2I31f8vO1AoPi5xyPFz/Hq8nlFHK3t4J6/HOaxz6yMq8zRUJNCfgjcR6Dl1R0EkkOEEEEYmvTh1xqTgqYud39dVF6aHYvJyz+/Ev2Oci5voMtHNNPwO129Jz8fqOXkMCc/b66Q4ud4k2Izc+eGOdy99RB/P1THlYtKYj2kfqEEtDuAu4Ene79+Dvha2EckRJIamvRhM5vw+g08A9pExeJcs2gfuDl28XPin/yc7G5YMZ3fvH6Gb287wvq5hXFzZlooWY7dwF1KqXStdVcExyREUhqaHl+QYae61YnFrNBaT3jPbDxF2QMP3Nx1qoUndldR2+GkJDOF65dPY0XZqOf7hqSl28PfD9Wx7cDwxc+bK5Ln5OeJsJpN2CwmbGYTdqsJuyU+gsVAFrOJr79/ATc+9AYPv/wut186O9ZDAkLLclxFoDNIOjBdKbUEuG1AOr8QYhRD0+PNJkVOqpW8NNuE98xCLcr2G5rmLjdd7kDm4K5TLfzkheNYTIpMh4Xmbjc/eeE4X6J8QkFt9OLn5D75eTRKKaxmhc1swmo2YbWY+r9OlIC+enY+GxcWcd+LJ/jQBaUUZ8W+z2MoS44/AjYCTwNorfcppdZFZFRCJKHh0uPv3rIgLEkfwRZlG4am0+WjzTk46eOJ3VVYTIqU3qWjvnq0J3ZXjSug1bW72Hawdtji58sXFLFlEpz83MekFHZr34zLjK13BpYMvrp5AZf/6CW+/+xRfvjRpbEeTmhZilrrqiF/PcT3eepCxJlIHcw5VlG239C0O710OL3DNhGu7XCS6Rj8duCwmqjrCL69VX/x84E63jzdysBnmSwnP5uUCiwX9n7YLfG5ZBgu0/NS+czaWdz34kk+fvEMLpieE9PxhBLQqnqXHbVSygZ8ETgSmWEJIUIxUvuqqdkptPd4ae0ZvRt+SWYKzd3u/hkagMtrUJw5dtHyZC5+NimFw2omxWrGbjXFTXJENH1h/Wx+v6eabzx9iD99YTWmGNYHhhLQPgf8BJgKVBPIcrw9EoMSQoRmuPZVbp/Bhy4opbnbPeb9r18+jZ+8cByn14/DasLlNfAZmuuXTxv29qMVP19clsfmJC1+NpsCAcxhMeOwJffsK1hpdgt3bZrHnb/bxx/fquG6C2N3EGhQAU0pZQY+obW+McLjEUKMw8D9uTMt3RRlOPjIsmlcODO4JaAVZbl8iXKe2F1FXYeT4mGyHCdj8bPFZMJhDex9OeI04zAefGDpVB597TTfe/YoVy4qJt0em54bQT2r1tqvlLqGQGKIECIOLZuVy+yi9HEXRq8oyx02AWS04ue1vcXPS5Og+LkvecNuMffufZmwmJN3vy+cTCbFf75/AR/831e578UTfOXKeTEZRyhh9BWl1L0ECqv7f6u11m+GfVRCiKBorel0+2jv8eL1h6/Dx2jFz7Py09hSUcz75heRFeHi50jVximl+oOWrTdxI5mTVaLh/Ok5fGDpFH758rvcdPHMmKTxhxLQVvVe3jPgOg1Ev0ePEJOc1poOVyCQhfNIl3gpft51qoUH/3mKyuZuLGZFfppt3LVxfZmHErzC57E3zgx7/ezCDHz+Wu54/C0+eP7UCT3Hx1ZOD/k+oXQKuXS07yulbtJaPxLyCIQQQTMMTYfLS7vTy2snmsMyexnr5OfNFSVRLX7uK/Ju7nJjVqANaOh0U5jhwGJSI9bGmZTC2tthw2Y2YbWoQNGyLBtGTW6ajRWzcnnj3WbWzs4nPyO6+6nh3Ln7EiABTYgQBdOyamgdWTg6e4xW/LxhQRGbY1T83FfkbWiNyaRQKAwDWns8lOakDKqNc1jNpNrMpNosMuuKE+vnFrD3dCv/OFLPDStCn2VNRDgDWmLvCAsRA30tqzw+P50uH3XtLt4808rt68/ji5fPwec3AoHM5UPriXf2GKv4eXNFCWtiXPzcV+RtNZvwGRqlQCl6GzlrpuWkUpKVgs1iSrqygGSQ4bCyenY+L77TwGUdLooyo7eXFs6AFr3Dk4RIEg/sPIXH56e524MJhcWk8GvNvS+eYGZeGhXTsgcFsj6hdvaobO7mbwfqeO5wPe0DDu2Mx+LnkqwUWrvd5KfbqWt3oXvfWixmE1rD7ZfOnnS9HxPN6vPyePlEI/883hTVujSZoQkRQ1WtPXS6fJhQmEyBrvsmBV6/5pevVPLDjy4Z9n7BdPYYrfh55aw8tiwuZuWsvJjOcgY26bVZAp02/vV95Xz96UPYLIop2Q7qO9x4DYOy3DTu2jQ/6geeitCl2i1cOCOX3e+2sGFB5LNh+4QzoL0SxscSYlKYlpNKXbsrkPygA8schgabWY3aR/H8aVn8ZtcZDENjs5hIs5mxWsx8dFkp79R1su1ALc8PKX4uyXKwpaKEKxYWkR+j4melAsukjt42UcO1ihraxPn86TkxOblbTMya2fm8caqZ1042Re0Q0JACmlJqC7AQ6F8U1Vrf03v5L+EdmhDxYTznjAXDMDQfXzmdvadb+k+v1hrQkJFiGbGP4q5TLTx7uJ7sFAudLh9un4HXb3Dh9HQeeuXduCp+7qv3SrGaSbEFCpaDSfmPVBNnET25aTYWTc3ijXdbuHReYVS6rIRyHtr9QCpwKYFz0a4DdkVoXELEhZ9uP8Z9O07iNzR2iwmf3xj1nLFguLx+utw+ulw+5k/J5MYV0/nNrjP4/BqbWZGRYsFiNo/YR7EvISTdZsNusdDu9NLl9vFGZWv/baJZ/DyUbUAAc1jMMW1WK2Jr1Xl5HKhp52BNBxfOiHwn/pAKq7XWi5VS+7XW31BK/QD4Y6QGJkSs7TjawH07TmJojcWk8Pk1zd0e8tJs55wzNhav36Db7aPT5Tuno8cnVs1kbnHmOX0UAe58ct85dWbVbT1oranr8A3q4KGATRXFbF5UwvyS6J38LAFMjGR6bip5aTb2nm6Nu4DWt6Dfo5SaAjQDs8I/JCHiwwM7T+EzDKxmEwqFUoABnS5f/zljo/H6DXo8frrdPlze0Y8OHNpHcWidWVOXi+/9/SglWSmDasYAHBYTKTYzpdmp/NsVc8f1WkNhNQf2v1Js5v6Tt+NBpJaGxfgppVg6PZsXjjTQ6fKS4YjsakEoAe2vSqls4L+BNwnsXz8UkVEJEQeqWnuwm034daAOCgKXbp9BaU7qsPdx+/z0uP10e3zjbhIM7y0rWkyK5i4PHS4fPkPT2uPtH0eazUxumg2twWdobhxHq6BgmE2qP3ilWM1x2bC3r57PalZkp1hp6HRNeGlYhMeCkkyeP9LA0bpOls+ceB/O0YQS0L6vtXYDf1BK/ZVAYogrMsMSIvam5aTiNwyau7wYBAp8/VpjNiluW1cGBBI7XD4/Lm9gSTEcDYI9PoN3m7vw+gx6vIMfz2ZW/PuV87CbTfzhzZoRj3oJxXANgNfMySfVZiHVNnwmYrx5YOcprGbVf8Bpqs1Cj8cX8tKwCL/iTAdZKVaO18dXQHsNuACgN7C5lVJv9l0nRLLpOzQzLx3ae7y4/QYWk4nPrSvj/Ok5nG1z4vYZwxY+j8fg4uf3Tn42mxRZDgt2i4mizBQu632DXl2eP+Hn7Fva7JvZtDk93LfjBMVZjoQKBFWtPWQPSX5JsZqDWhoWkaWUYkZeKu82daO1juje7pgBTSlVTOCU6hSl1Pm8V0CdSSDrMeqUUlcSOD3bDDyktf5uLMYhktvAeqgq1U1xVmD2csGMnKBOgQ7GaMXPVrOJDLuF7FQLbp8e9QTpUFnNge7zf3izmhSriTR7IBjYLOaEnNlMy0mlodPVP0ODwM92pKVhEV0zclPZX91Ou9NLdqotYs8TzAxtI3AzUAr8cMD1ncB/RGBMo+o9Pfs+YANQDexWSj2ttT4c7bGI5OX1Gzi9fhZMyeR7H1o84SNaBi7rFWc4WFuez+mWnmGLnzdXFLNxYTGnGrpHPUE6FObevo+O3r2wvg70tR2uEWc2iZRk0Teb7vH4+vtaev26f2lYxFZhbz/Hpi5PbANa75EwjyilPqS1/kPERhK8FcAJrfUpAKXUE8A1gAQ0MW5aa9y+wD5Yj8cf1sMy+5b1TAT23A7XdbCvpr3/+1azYs3sfLYsLukvfg7HwZZWs4k0++j7YCPNbNLtloRKshjaXaQ0zgPwZNPXmaa5281s0iP2PMEsOd453Od9tNY/HHpdhE0FqgZ8XQ2sHHgDpdRngc8CTJ8e3eMLRGLQWuPyGvR4Ap02wrkXNvR5Hnr5XdqdXpwe/6AO3g6LiVvXzjqn+Hm8R8MopXBYTaRaLaTazUGdAzbSzMZq0gmXZDGZu4sMfM/LL57YwZqRkGYP/EE1cDUiEoJZcsyI6AhCN9yO4qB3Iq31g8CDAMuWLZNTAASGEZiBeXxGILXe48eIQADr09Lt4blDdWw7WEd163s9GZWCDLuFrBQLHp/BtRec24k8lKNhTEqRag+cB5ZqDb2oeaSZzde2HpQkiwQy8D2vbP7iuHvPs5hMWExqzHrMCT/PWDfQWn8joiMIXTUwcGe8FDgbo7GIOKO1xuMPBC6Pz2DnO408+vppzrZP7FTnYPSd/LztQB2vnWrGP+DoZ5tZkZNqI8NuwWRSOL1+SrKGT1gY62iYviCWbreQYjWHnDW242gD33v2KKeaAj0fZ+WlntPFftpOSbIQ4WVSigj+DQmE1stxDvBzoEhrvUgptRi4Wmv9rYiNbni7gXKl1CygBrge+FiUxyDiRN/el9Pjx+n1D1o6DMepziMZuMeVm2qjJCuFAzXtNHa9l/2Y6bBw+YIipuek8uSewKxLmQKBYbSMxZGOhpmSnUJBhp00m2Xc7aV2HG3gy0/to7XHS99DnGjs5t+e2sf/XLekP6hJkoUIJ601Xr+BxRzZrjKh1KH9Avgy8ACA1nq/UuoxIKoBTWvtU0r9C/B3Amn7v9RaH4rmGETsuHuLmN0+Px6fgdevR9z7Gu+pzmPZdaqFHz9/LJAJ6fFT3+HmSF1n//eHO/m5ONMRdMbi9cun8ZMXjuPy+Umxmvs7jnzxsvKQWwcNzVRs6/HQ6fJhNqn+zvvK0HS5B++PSZKFCCe3z0ADjgh33A8loKVqrXcNWd7wjXTjSNJabwO2xeK5RfT0zb7cXqO3G4d/0DLeWEI91TkYp5u7+cE/jtHc7WbgUEwKCtLt/OAjS4Y9+Xlor8aRKKW4dH4hWalWHnm1kpo257iDyXDtoCqbA42N+wJt4DkDy6VD98cmc5KFCK+23pZt2anx08uxSSl1Hr0JGEqp64DaiIxKTCp9+15ev8bbe7aX22fgM0aefQUjmFOdg+H0+nmpt/j54IDiZwj0U8xKsZJqM9Hl9g8bzILhsJpJd1hIs1kwmxSbKkrYVDGxQxGHawdlNSvcXo0e0J9S60CdmuyPiYn42Ci9RJ89GAgVN6yYzpJp2REbQygB7XYCWTTzlFI1wLvAjREZlUhKfTMub2/ShtcfWFcPZ83XQH1Ld06vH4fVhMtrBN1tQ2vNsfouth2o5bnD9bgHNBq2mBRpNjN5abb+Rr1Orz/kQGk2KdLtFjIc1kEzpnAZrh1UUYadqlYnfkOjVeCPBUNDtt0q+2MiYvZXt2MxKeYWRzZpPqQ6NALLfC8CJqAb+BCDu4cI0c/jCywVugfseUXTirJcvkR5SN02Ol1enj/SwLYDdZxo7Br0vXR7oMuG12+glMJraMxmHVKgVCqwpxeYjYWeoRiK4YqmLWYTcwrTUUr1ZznOzj83y1GIcHrzTCvzSzIj3ug6lDq0ucByYCuBWrBPADsjNC6RYDw+Y1C6vMsb2TqvYAWzd6W1Zn9NO8/sr2Xn8aZBgddhCZz9lZdm6z/3y+n1YzUpMlNsQQdKk1JkpljJSrFG7fywkTIV796yIKGDVyK15BLQ5fax93Qrt6yJ/PGZQdehKaWeAy7QWnf2fv3/gN9HdHQirvQtD/oMja93udAThr2uWBmp+NlhNXHZ3EK2LC7hG385RFaKFTWgnt9hNdHp8vHQzUvGfA6lAiUD2am2qB+EmYyZinLuWeJ5+XgjXr/m0rmR//cJZQ9tOjDwqFwPMDOsoxFxwTACSRp9nTU8fgOvz4iLGddEjVb8PL8kg82LSrh0XkH/Mt2UrNRxJZaYlCIjRoFsoIGZin0zm69tPZiwMxs59yzxPL3vLHlpNi6ckRPx5woloP0a2KWU+hOBTMcPAo9EZFQiaoYuFQZmXNHd64qGunYXzx6s428H684tfp5fxOaKYsoKzm2aGmpiidmkyHRYyYzi0mIwkmVmI+eeJZb2Hi/bDzdw40XTg+otOlFBBzSt9beVUn8D1vZe9Smt9VuRGZYIN7+h+wOW2z92UXK4haN7fKi8foNXTjSz7UAte0+3Dmr4ef70bLYMKX4eTrCJJTaLiawUK+l2y6iJHrHa/0mWmY2ce5ZYnth9Bo/f4LoLz+1ZGgmhzNDQWr8JvBmhsYgwGNrLMLBcqGM664pkC6rhnG7uZlv/yc/e/uvz0mxcuaiYKxcVMzWEerGREkuUCqTvZ6ZYg8reiuUsKVlmNtKSK3F4fAb/90olq2fnsXBKVlSeM6SAJuLLoA7yMZh1BStSLagGGqn42aRg5aw8NlcUc1FZXliWAS0mExkOCxkOS38dWjBiOUtKlplNMia6JKun952lrsPF965bHLXnlICWIDwDCpI9/YXJibHXFYkWVDC4+PmFow10j3Dyc9/hghNlNZvISbONu35stFlSpJcik2lmIy254p/Wmgd3nmRecQbryvOj9rwS0OJMX2r84MAVf7OuUISrBVWfvuLnZw7UcrKxu/96q1mxtryAzYuKWTo9u7/57kRZTCay06xkjLE/NpaRZklpNnPElyJlZiOi6bnD9Ryr7+KHH1kS0eYBQ0lAiyG/oXH3dtLo66iRDKnxQ02kBVUfrTX7q9t55sC5xc8z81LZXFHChgWDT36eqL2VrfxuTxVn251Mz02bcAAYaZZkM5uishQpMxsRDT6/wfeePUpZQRrvXzIlqs8tAS1KBtZ29QWxRFkynKjxtKDqM1bx8+aKEuaXZIz5V2AoWZZWs4mD1e389IXj2CwmclJtYZk1yenQI5PuH8njid1VnGrs5sFPXBiVVP2BJKBFQN/Mqz9NPoH2uyIl2ONTIPTi57EEm2VpNimyU2xkplj496f2Y7OYwj5rGm6WNNlPh06WGjkRaHP14+3HWDEzlw0LiqL+/BLQJqhvv8ud5IXJ0TBS8XOGw8KGBUVsXjR88fNYxsqyNJsUWSlWMh3W/pOgo5nmnlye/Y4AABeKSURBVEwJG+ORLDVyAh586SRNXR4euml+VPfO+khAC9Jw9V3JuucVTWMVP29eVMza8oIJHa8yUpZlfYeTvHQ7mY5zkz2imeY+2RM2kqVGbrI73dzN/TtP8f4lU1gawTPPRiMBbRh953b1Ba2+QJbImYbxpq/4+R+H62kbUvy8cWERmypKQip+Hs05WZYKvD7NjLy0EZNIoj1rmswJG8lSIzeZaa35f08fwmY28bUt82M2jkkf0Aa1hPL5Zb8rgsYqft6yuJiVs8JT/DxQX5aly+cn1WbG4zMwNHzukvNGvM9knzVF02Rfck0Gh2s7ePGdRu6+agFFmY6YjWNSBTTfkA7ybq/sd0Wa1prjDV08s7+W54820DOk+HnTokDxc0FGeIqfh7PyvDzuspl5fNcZatqcQQenyTxriib54yGxuX1+/rq/lnnFGdx08YyYjiXpA5qhNbXtTjw+Y1C2nIisLpeP7Ufqzzn52WpWrJmdz5aKkrAWPw9HKUWa3UxOqo1Z+YHlxb7U8Ad2ngIkiy5eyB8PiWt7b8/Uh29aFlIruEhI/oBmgHPArEBETt/Jz9sO1PHSscZzip+3LC7h8vnhLX4eSarNQl66rb8ORlLDhQi/My09vHqymRWzclk2M7KnZwQj6QOaiLxwFT+Hg0kpctJs5wRNSQ0fTAqZxUT5/AZ/eLOazBQrVy4sjvVwAAloYpzGKn7etKiEy0Iofg4Hh9VMQYZ92O4EiZQaHulgI7NVMR4fWzl90Nc/eO4dGjvd/N+nlnPp3Pj4vZGAJkIy3pOfI2mkWdlAiZIaHo1gI7NVMVGHz3bw8x0nufb8qXETzEACmgjCaMXPS6cFTn5eWz76yc+RkmqzkJ9uG3MzOlFSw6MRbBJptirij9dv8OWn9pGdauXuqxbEejiDSEATIwr3yc/hZDWbyEu3Bb2kuX5eIddVt/HQy+/S7Qkc2XLrmlkjBolY7TFFI9gkymxVxKcHd57i0NkOfn7jBeSk2WI9nEEkoIlBonny83iYlCI71UpWijWkJJMdRxt46s0aCjLsTO+doT31Zg2LS7PPCVSx3GOKRrBJlNmqiD8nGrr4yfPH2bSomE0VJbEezjkkoAkAjtV38syBWl44cu7Jz9Eofg5GZoqVnFTbuIJpKEt5sdxjikawkUJmMR5+Q/OVP+wnxWrmG9csjPVwhiUBbRLrcvl4/mg9zxyo40RDbIqfg5Fut5CdapvQHl0oS3mx2GMauMSZYbegtabd6Y1YsJFCZhGqR1+rZO/pVn7w4SUUZsSuvdVoJKBNMmMVP0fi5OfxSrVZyEmzYreYJ/xYoSzlRXuPaegSZ9+s7JvXLJKgI+JCS7eH+148wSVzCrj2gqmxHs6IJKBNEi3dHp47XM+2A7UxL34ei81iIi/NTopt4oGsTyhLedHeY5I0ehHPtNb8+a0azCbFf11bERfvESORgJbE/IZm7+lWnjlQy6snBxc/zyvOYEtFaCc/R5rFZCInzUqGI/yzw1D2jaK9xyRp9CKe7T3dyonGLr75gUUxy2oOVny8k4mwqusIFD8/e7COhs4hJz/PL2JTRTHnRbn4eTRKKTIdFnJSbf0nRkdCKPtG0dxjkjR6Ea86nF62HaxlZl4aN66YPvYdYkwCWpLw+g1ePRkoft5TGZmTnyPBYTWTn26Pu3FFk6TRi3iktWbrvrP4/JprL5ga0T82w0UCWoI709zDMwdqzzn5OTfNxpULi9i0qISpOfG3TGBSitx0G5kRWF5MNONZ4pTmwiLSDtS0c6S2gysXFpOfHtuSnWBJQEtATq+fnccCxc8HauKv+HksaXYLeWljt6uaTEJZ4pTmwiLSut0+/rLvLFOzU1g9Oz/WwwmaBLQEkgjFz6OxmALtqtLs8ms3XjuONvDFJ96i2+PDYQmcLpDhsEpWpAirZw7U4vT6ueWCqXH7h/Fw5J0lziVK8fNYMlOs5EY46SPZ9c3Mejx+LCaFz9CcbXMxJTtQfC5ZkSIc3qnr4O2qNi6bV0hJVvxtV4xGAlocGq34eUZeKlsqStgwv4is1Pjff7JZTOSn23FYw1dTNln11avZLSZ8fo3JpDDQNHa6MZuUZEWKCXN5/fz57bMUZthZP6cg1sMJmQS0ODJi8bPFxPq5hWxZXMyCksy4LmzsY1KKnFRbQgTdWAs2waOvXi0/3c7ZdicYgNK4fIZkRYqweO5wPR1OL5+75LyE3OOWgBZjfkOz53Tg5Oehxc9z+4qf5xYk1L5Tqs1CXrpt2JOjxWChJHj01atl9hZhN3W5cfs0aTYL91y9UPbPxITUtDp541QzK8vymJabmLP9xHmXTDIjFT+n2y1sWBA4+Tmeip+DYTYp8tLtpCdQ8I21UNpeDaxXy3BYsJgVXr+WYCYmzNCarftqSLMHmi8kKnnniaKxi59jd/LzRGU4rOSmje9ol8kslLZXcuyLiJTdlS1Utzr5yLLSsPZQjTYJaFFwprmHbQdree5QYhU/B8NuNZOXZpOkj3EKte2VHPsiwq3L7eO5Q/XMyk9jSWl2rIczIRLQIsTl9fPSCMXPK2blsnlRCReV5SbkxisElhdz0qTTx0RJ2ysRa88fqcft83P1kikJkXA2Gglo/3979x4jV3necfz7zMzOZW/emxcv9uILmBgTI9wYU1IVTEKD47ShDqQyJKU0NFBECm1TKaS0CIFoaNoGkZZGkChKGind0lQ0LhhSYXCsplgYhMHYMbZjTNZgfAff7b08/WOO08Ee747tnTlnzvl9pBFnzpzd/el4Oc++73nP+46x9dv2sWT1uyz9+bYTHn6e/+EJzI/4w8+VaM5n6GzKqXtxDKgbUWrhhkvLTyy8cfs+/vrHr3PDpefwZ791fo1TjT0VtDEQl4efR9OQLj5TVs997FGkbkQJy9eWrKOxIc2fXlX/xQxU0E6bu7P67fd5KgYPP1eitdBAZ1O27rskRKTofzfuZOm67Xxl/oy6mXx4NCpop2ikh5+vnNHNp6q88vOLm3bTt7KfrXsP0dNaYNElvcyd1lGVnwVqlYnEkbvz4DPrmNhW4A9/Y0rYccaMCloFRlv5eUGNHn5+cdNuHn5uA5lUcUHMXQeO8PBzG7iT6WNe1MyKD/q2NTaoVSYSM8ve2MFrW97nb6+dFasRyipoI4jaw899K/vJpIxC8At4bFRc38r+MS1oTbkMHU2a6UMkjtydh5duYGJbgYWzJ4UdZ0ypoB1npIefL+5tY8GsCfzmeV3kQvirZuveQ7TmP/hPlm9I8e7eQyf5ilOj7kWR+Fu+YSer+t/jbxbOqstJHEYSuYJmZvcCXwR2BLv+0t2XBJ99FbgZGALucPefjNXPrYeHn3taC+w6cORXLTSAwwPDTGg9s1zHJhJuLWTUvSgSc99cuoGzx+W59iMTw44y5iJX0AIPufvfl+4ws5nAIuBC4GzgWTM7392Hyn2DShwOVn5+avW7rH77/V/tP/bw86dm9XDp1Og8/Lzokl4efm4DhwaGyDekODwwzOCws+iS3tP+ns35DB2NWj1aJAlW9b/Hy2/t4Z7fnkkuE7+emKgWtHKuAfrc/QjwppltBOYCL5zqNzrZw88TWvN8clZ0H36eO62DO5lO38p+3t17iAlnMMqxIZ1ifIvWKRNJkn95YTNN2TSfnROve2fHRLWgfcnMbgReAr7s7nuAicCKkmO2BPsqMtrDzwtm9TC7Dh5+njut44wGgGj0okgy7dp/hCdf3cqiub20xHTKulAKmpk9C0wo89HdwLeA+wEP/vsPwBeAcldfL7MPM7sFuAWge+JkHnx6HT9dv4MjpQ8/dzSy4KIePhGjh59Hk80UW2Vx7GoQSbLSa17XhPJ/5/et7Ofo0DA3Xja5ltFqKpSC5u5XVXKcmX0beDJ4uwUovVk0CXjnJN//MeAxgFzPdP/vtduA/3/4ecGs+ln5eay0NWZpV6tMJJZKr3nTLrjohD/03Z2+lb/ko+d2cl53S83z1UrkuhzNrMfdtwZvFwKvB9uLgR+a2TcoDgqZDrxYyfes15Wfx4LulYnIK/3v0b/7EHd+PB5zNp5MFK/uXzeziyl2J24GbgVw9zVm9jiwFhgEbq9khOPUria+9blfq2Lc6GrJF+dfTGlWfJFEW7zqHbKZFFdfWL+rUVcicgXN3X9/hM8eAB44le+XT+D9onTK6GrOJa41KiInGhp2nnxtKx+f0R3bwSDH6IoXM1qrTERKvbnzADv3H+Gai88OO0rVqaDFhKatEpFy1r27l1wmxbwPxX/NPRW0OpdOGW2NWVrzmrZKRE60fts+Lju3MxEDw1TQ6pRZcQmZ9kYN+hCR8nYfOMrO/Ue54vzxYUepCRW0OpRrSNPVnNUD0jGybN12Hl2+if49B+ltb+TWy6cxb0b8u4ikutZv2weQiO5GAM1IW0fMjM6mHBPbCipmMbJs3XbuWbyG7fsO01ZoYPu+w9yzeA3L1m0PO5rUuY3b99Pe2MCUzsawo9SEClqdaEin6BmXT8w0XUny6PJNNKSNxmzxPmhjNkND2nh0+aawo0kdc3fe2n2QKZ1Nibm/ri7HOtCcz9DVlNO9spjq33OQtsIH/1ApNKTZsudgSImknnU0Zbnh0nPYvPMAB44Mcl1MZ9YvRy20CCu2ygp0t+RVzGKst72RQwMfnPTm0MAQk9qT0U0k1fHyW3sAmDP59FfnqDcqaBGUMqOjKcuk9oKeK0uAWy+fxsCQc/DoIO7F/w4MObdePi3saFLHXv7lHlryGaZ3N4cdpWbU5RgxzbkMHU1aQTpJ5s3o5j6K99K27DnIJI1ylDGw5u33mTVxXKJ6d1TQIkKz4ifbvBndKmAyZoaGnTe27eNzl8Z37bNyVNBCZmaMKzRorTIRGTNv7jzA4YFhLuhpDTtKTamghUgPSItINazduheAC3riu5hnOSpoIUinjPamLK0xX8pBRMKxYds+UgbnJWhACKig1VxLvoGOpqyWdxGRqtm8qzi4KGm9PypoNZJrSNPZlNWgDxGpus07DzClqynsGDWnglZlKTPaG7OaskpEambzrgPMPqct7Bg1p4JWRY3ZDJ3NWRr0TJmI1MjgsLPv8CCTO9VCkzGQThmdzTmaczq9IlJbRweHAZjalbyp03TFHWPN+QydTTkN+hCRUBwdLM4LqhaanLaGdIqu5pzmXhSRUB0dclLAxLZC2FFqTgVtDIwrFIfia6YPEQnb4NAw4wsNiRxRrYJ2BjT/oohEzcCQc1ZrLuwYoVBBO01qlYlIFA0MDXNWaz7sGKFQQTtFapWJSJQNDjvjW9RCk1G0FhroaMwman0hEakvg2qhyUgyqWKrTCMYRSTqHOhWC03Kac5n6GrKqVUmInWjvTEbdoRQqKCdRCaVoqslS2NWp0hE6su4QjLnjtXVuoziCtK6VyYi9alVBU20grSIxIFaaAmmFaRFJE5U0BLIzGjJZ2hv1ArSIhIfKmgJU8im6WzKkc1orTIRiQ+DxF7XElfQ0imjoylLi7oXRSSGkjyYLVEFrTmXobNZa5WJSHylSO71LREFTWuViUhSJHm+9NgXtHTKmNRe0Kz4IpIIqQRf62J/59AMFTMRSYwkX+5iX9BERJJELTQREYmFJI95U0ETEYmRJN9iUUETEYmRBNczFTQRkThJcD1TQRMRiZUEVzQVNBGRGLEEVzQVNBERiQUVNBGRGNGgEBERkTqngiYiEiMJbqCpoImISDyooImISCyooImIxIimvgqBmX3WzNaY2bCZzTnus6+a2UYze8PMri7ZPz/Yt9HM7qp9ahERiaowW2ivA58BlpfuNLOZwCLgQmA+8M9mljazNPAI8ElgJnB9cKyIiASS2z4LccVqd/85lG0eXwP0ufsR4E0z2wjMDT7b6O6bgq/rC45dW5vEIiJ1IMEVLYr30CYC/SXvtwT7Trb/BGZ2i5m9ZGYv7dixo2pBRUSioPSaN3x4f9hxQlPVgmZmz5rZ62Ve14z0ZWX2+Qj7T9zp/pi7z3H3OePHjz+d6CIidaP0mnd2d1fYcUJT1S5Hd7/qNL5sC9Bb8n4S8E6wfbL9IiKScFHsclwMLDKznJlNBaYDLwIrgelmNtXMshQHjiwOMaeIiERIaINCzGwh8I/AeOApM1vl7le7+xoze5ziYI9B4HZ3Hwq+5kvAT4A08F13XxNSfBERiZgwRzk+ATxxks8eAB4os38JsKTK0UREpA5FsctRRETklKmgiYhILKigiYhILKigiYhILKigiYhILKigiYhILKigiYhILKigiYhILJh72fl9Y8PMdgBvhfCju4CdIfzckShTZZSpMspUmbHItNPd51dyoJk9U+mxcRP7ghYWM3vJ3eeMfmTtKFNllKkyylSZKGaKK3U5iohILKigiYhILKigVc9jYQcoQ5kqo0yVUabKRDFTLOkemoiIxIJaaCIiEgsqaCIiEgsqaGfIzOab2RtmttHM7irz+Z+b2Voze83MlprZ5Ahk+mMzW21mq8zsf8xsZtiZSo67zszczKo+zLmC83STme0IztMqM/ujsDMFx/xe8Du1xsx+GHYmM3uo5BytN7P3IpDpHDN73sxeCf7fWxCBTJODa8BrZrbMzCZVO1PiuLtep/kC0sAvgGlAFngVmHncMVcCjcH2bcC/RSBTa8n2p4Fnws4UHNcCLAdWAHPCzgTcBPxTxH6fpgOvAO3B++6wMx13/J8A3w07E8WBGLcF2zOBzRHI9O/AHwTbHwN+UKvfraS81EI7M3OBje6+yd2PAn3ANaUHuPvz7n4weLsCqPZfZZVk2lvytgmo9sigUTMF7ge+Dhyucp5TyVRLlWT6IvCIu+8BcPftEchU6nrgXyOQyYHWYHsc8E4EMs0Elgbbz5f5XM6QCtqZmQj0l7zfEuw7mZuBp6uaqMJMZna7mf2CYgG5I+xMZjYb6HX3J6ucpeJMgWuDLqIfmVlvBDKdD5xvZj8zsxVmVu0pjir+HQ+606cCz0Ug073A581sC7CEYssx7EyvAtcG2wuBFjPrrHKuRFFBOzNWZl/Z1o6ZfR6YA/xdVRNVmMndH3H3c4GvAH8VZiYzSwEPAV+uco5SlZyn/wKmuPtFwLPA9yOQKUOx23EexdbQd8ysLeRMxywCfuTuQ1XMA5Vluh74nrtPAhYAPwh+z8LM9BfAFWb2CnAF8DYwWMVMiaOCdma2AKV/tU+iTNeGmV0F3A182t2PRCFTiT7gd6uaaPRMLcCHgWVmthn4dWBxlQeGjHqe3H1Xyb/Xt4GPVDFPRZmCY37s7gPu/ibwBsUCF2amYxZR/e5GqCzTzcDjAO7+ApCnOElwaJnc/R13/4y7z6Z4PcDd369ipuQJ+yZePb8o/rW8iWI3y7EbwRced8xsijeLp0co0/SS7d8BXgo703HHL6P6g0IqOU89JdsLgRURyDQf+H6w3UWxm6sz7H874EPAZoLJGiJwnp4Gbgq2L6BYXKqWrcJMXUAq2H4AuK/a5yppr9AD1PuLYnfG+qBo3R3su49iawyKXVXbgFXBa3EEMj0MrAnyPD9ScalVpuOOrXpBq/A8fS04T68G52lGBDIZ8A1gLbAaWBR2puD9vcCD1c5yCudpJvCz4N9uFfCJCGS6DtgQHPMdIFer85WUl6a+EhGRWNA9NBERiQUVNBERiQUVNBERiQUVNBERiQUVNBERiQUVNBERiQUVNBERiQUVNJEKmNl/mtnLwRpktwT79pd8fp2ZfS/YPsvMnjCzV4PXR0OKLZIombADiNSJL7j7bjMrACvN7D9GOPabwE/dfaGZpYHm2kQUSTYVNJHK3GFmC4PtXkaeEPhjwI0AXpx5XhPQitSACprIKMxsHnAVcJm7HzSzZRRnby+dNy4fQjQRKaF7aCKjGwfsCYrZDIrL2wBsM7MLgnW2FpYcvxS4DcDM0mbWiohUnQqayOieATJm9hpwP7Ai2H8X8CTFFZq3lhx/J3Clma0GXgYurGFWkcTSbPsiIhILaqGJiEgsqKCJiEgsqKCJiEgsqKCJiEgsqKCJiEgsqKCJiEgsqKCJiEgs/B/bxDE4Qqj0kAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.jointplot(x='auc',\n", " y='delta_recall',\n", " kind=\"reg\",\n", " data=combined_df)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "# Trained Classifier Dataset\n", "import classiflib\n", "\n", "retrained_classifiers = glob.glob(os.path.join(report_db_location, '*retrained_classifier*'))\n", "training_classifier_summary_data = {\n", " 'subject': [],\n", " 'montage': [],\n", " 'training_auc': [],\n", " 'classifier': []\n", "}\n", "for classifier_loc in retrained_classifiers:\n", " metadata = extract_report_info_from_path(classifier_loc)\n", " subject = metadata['subject']\n", " montage = metadata['montage']\n", " \n", " # Only add the retrained classifier if one has not already been added\n", " # for the given subject\n", " if subject not in training_classifier_summary_data['subject']:\n", " classifier = classiflib.container.ClassifierContainer.load(classifier_loc)\n", " training_classifier_summary_data['training_auc'].append(classifier.classifier_info['auc'])\n", " training_classifier_summary_data['subject'].append(subject)\n", " training_classifier_summary_data['montage'].append(montage)\n", " training_classifier_summary_data['classifier'].append(classifier)\n", "\n", "trained_classifier_df = pd.DataFrame.from_dict(training_classifier_summary_data)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "combined_df = combined_df.merge(trained_classifier_df, how='left')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's again check that the merge has worked as expected" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "assert len(combined_df) == len(session_df)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentmontagesessiondatesession_summaryelectrode_typelocationpulse_freqstim_durationdelta_recallclassifier_summarypvalueaucclassifiertraining_auc
\n", "
" ], "text/plain": [ "Empty DataFrame\n", "Columns: [subject, experiment, montage, session, date, session_summary, electrode_type, location, pulse_freq, stim_duration, delta_recall, classifier_summary, pvalue, auc, classifier, training_auc]\n", "Index: []" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "combined_df[combined_df['classifier'].isnull() == True]" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentmontagesessiondatesession_summaryelectrode_typelocationpulse_freqstim_durationdelta_recallclassifier_summarypvalueaucclassifiertraining_auc
0R1275DFR5012017-06-10 17:15:09.603000+00:00<FRStimSessionSummary(_events=[ ( 1, 1, 'R1275...Right Cerebral White Matter200.0500.035.411540<ClassifierSummary(_predicted_probabilities=[ ...0.0150.616795<classiflib.container.ClassifierContainer obje...0.716635
1R1275DFR5022017-06-11 17:29:07.115000+00:00<FRStimSessionSummary(_events=[ ( 1, 2, 'R1275...Right Cerebral White Matter200.0500.0-7.532247<ClassifierSummary(_predicted_probabilities=[ ...0.4100.515080<classiflib.container.ClassifierContainer obje...0.716635
2R1292EFR5002017-07-03 16:12:42.221000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1292...Left Middle Temporal Gyrus100.0500.05.376644<ClassifierSummary(_predicted_probabilities=[ ...0.3900.515818<classiflib.container.ClassifierContainer obje...0.558089
3R1304NFR5002017-05-29 15:35:00.879000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1304...Right Cerebral White Matter200.0500.055.067482<ClassifierSummary(_predicted_probabilities=[ ...0.3150.526984<classiflib.container.ClassifierContainer obje...0.503983
4R1308TFR5002017-06-19 20:33:04.856000+00:00<FRStimSessionSummary(_events=[ ( 1, 0, 'R1308...Left Middle Temporal Gyrus200.0500.0NaN<ClassifierSummary(_predicted_probabilities=[ ...0.0000.705267<classiflib.container.ClassifierContainer obje...0.599638
\n", "
" ], "text/plain": [ " subject experiment montage session date \\\n", "0 R1275D FR5 0 1 2017-06-10 17:15:09.603000+00:00 \n", "1 R1275D FR5 0 2 2017-06-11 17:29:07.115000+00:00 \n", "2 R1292E FR5 0 0 2017-07-03 16:12:42.221000+00:00 \n", "3 R1304N FR5 0 0 2017-05-29 15:35:00.879000+00:00 \n", "4 R1308T FR5 0 0 2017-06-19 20:33:04.856000+00:00 \n", "\n", " session_summary electrode_type \\\n", "0 \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subjectexperimentsession_summarysessiondelta_recall
0R1275DFR5{<FRStimSessionSummary(_events=[ ( 1, 1, 'R127...213.939646
1R1292EFR5{<FRStimSessionSummary(_events=[ ( 1, 0, 'R129...15.376644
2R1304NFR5{<FRStimSessionSummary(_events=[ ( 1, 0, 'R130...155.067482
3R1308TFR5{<FRStimSessionSummary(_events=[ ( 1, 2, 'R130...310.134115
4R1315TFR5{<FRStimSessionSummary(_events=[ ( 1, 0, 'R131...1-4.448695
\n", "" ], "text/plain": [ " subject experiment session_summary \\\n", "0 R1275D FR5 {