Skip to content

Instantly share code, notes, and snippets.

@ravila4
Created April 10, 2018 03:32
Show Gist options
  • Save ravila4/26b6ceb21c7e87af80be01f4620a7a58 to your computer and use it in GitHub Desktop.
Save ravila4/26b6ceb21c7e87af80be01f4620a7a58 to your computer and use it in GitHub Desktop.
Notebook for ROC/AUC and enrichment factor analysis
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# DUD-E screening results for receptor ID: *comt*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### **Note:** The scores in this analysis we obtained using only the top docking pose from Smina."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reading DLScore output"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>nnscore_01</th>\n",
" <th>nnscore_02</th>\n",
" <th>nnscore_03</th>\n",
" <th>nnscore_04</th>\n",
" <th>nnscore_05</th>\n",
" <th>nnscore_06</th>\n",
" <th>nnscore_07</th>\n",
" <th>nnscore_08</th>\n",
" <th>nnscore_09</th>\n",
" <th>nnscore_10</th>\n",
" <th>...</th>\n",
" <th>dlscore_01</th>\n",
" <th>dlscore_02</th>\n",
" <th>dlscore_03</th>\n",
" <th>dlscore_04</th>\n",
" <th>dlscore_05</th>\n",
" <th>dlscore_06</th>\n",
" <th>dlscore_07</th>\n",
" <th>dlscore_08</th>\n",
" <th>dlscore_09</th>\n",
" <th>dlscore_10</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>4.536008</td>\n",
" <td>5.396756</td>\n",
" <td>9.277304</td>\n",
" <td>7.494246</td>\n",
" <td>2.945559</td>\n",
" <td>1.765776</td>\n",
" <td>17.833014</td>\n",
" <td>3.187081</td>\n",
" <td>7.635727</td>\n",
" <td>1.553746</td>\n",
" <td>...</td>\n",
" <td>6.280250</td>\n",
" <td>6.189641</td>\n",
" <td>5.945060</td>\n",
" <td>5.594811</td>\n",
" <td>5.756034</td>\n",
" <td>6.155606</td>\n",
" <td>5.419800</td>\n",
" <td>6.898730</td>\n",
" <td>5.973953</td>\n",
" <td>6.864763</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>7.450654</td>\n",
" <td>11.730812</td>\n",
" <td>6.265531</td>\n",
" <td>1.070353</td>\n",
" <td>3.418288</td>\n",
" <td>1.639919</td>\n",
" <td>9.641711</td>\n",
" <td>7.422944</td>\n",
" <td>5.057000</td>\n",
" <td>5.830049</td>\n",
" <td>...</td>\n",
" <td>4.938908</td>\n",
" <td>5.334916</td>\n",
" <td>4.293091</td>\n",
" <td>5.267609</td>\n",
" <td>5.085522</td>\n",
" <td>6.309859</td>\n",
" <td>4.435493</td>\n",
" <td>5.776924</td>\n",
" <td>4.851845</td>\n",
" <td>5.440871</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>5.279146</td>\n",
" <td>12.310115</td>\n",
" <td>3.841030</td>\n",
" <td>0.835435</td>\n",
" <td>2.924070</td>\n",
" <td>3.057468</td>\n",
" <td>5.750599</td>\n",
" <td>2.682566</td>\n",
" <td>9.104843</td>\n",
" <td>-0.774676</td>\n",
" <td>...</td>\n",
" <td>4.803917</td>\n",
" <td>5.444674</td>\n",
" <td>3.924729</td>\n",
" <td>5.376449</td>\n",
" <td>5.361448</td>\n",
" <td>6.503062</td>\n",
" <td>4.328010</td>\n",
" <td>6.359634</td>\n",
" <td>4.692744</td>\n",
" <td>5.696967</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>7.545798</td>\n",
" <td>6.963721</td>\n",
" <td>6.307883</td>\n",
" <td>7.976962</td>\n",
" <td>1.827056</td>\n",
" <td>1.603488</td>\n",
" <td>13.498619</td>\n",
" <td>6.406145</td>\n",
" <td>5.750445</td>\n",
" <td>18.367493</td>\n",
" <td>...</td>\n",
" <td>4.375763</td>\n",
" <td>4.890980</td>\n",
" <td>3.990071</td>\n",
" <td>5.115571</td>\n",
" <td>4.450336</td>\n",
" <td>5.754834</td>\n",
" <td>3.881885</td>\n",
" <td>5.338937</td>\n",
" <td>4.852479</td>\n",
" <td>5.326446</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>4.421091</td>\n",
" <td>12.009138</td>\n",
" <td>1.918143</td>\n",
" <td>8.824357</td>\n",
" <td>2.579612</td>\n",
" <td>2.455463</td>\n",
" <td>5.053625</td>\n",
" <td>7.954544</td>\n",
" <td>7.999264</td>\n",
" <td>1.751066</td>\n",
" <td>...</td>\n",
" <td>6.104701</td>\n",
" <td>5.777875</td>\n",
" <td>5.372001</td>\n",
" <td>5.696952</td>\n",
" <td>5.494438</td>\n",
" <td>6.711897</td>\n",
" <td>4.867440</td>\n",
" <td>6.407499</td>\n",
" <td>5.344415</td>\n",
" <td>6.223027</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>5 rows × 30 columns</p>\n",
"</div>"
],
"text/plain": [
" nnscore_01 nnscore_02 nnscore_03 nnscore_04 nnscore_05 nnscore_06 \\\n",
"comt 4.536008 5.396756 9.277304 7.494246 2.945559 1.765776 \n",
"comt 7.450654 11.730812 6.265531 1.070353 3.418288 1.639919 \n",
"comt 5.279146 12.310115 3.841030 0.835435 2.924070 3.057468 \n",
"comt 7.545798 6.963721 6.307883 7.976962 1.827056 1.603488 \n",
"comt 4.421091 12.009138 1.918143 8.824357 2.579612 2.455463 \n",
"\n",
" nnscore_07 nnscore_08 nnscore_09 nnscore_10 ... dlscore_01 \\\n",
"comt 17.833014 3.187081 7.635727 1.553746 ... 6.280250 \n",
"comt 9.641711 7.422944 5.057000 5.830049 ... 4.938908 \n",
"comt 5.750599 2.682566 9.104843 -0.774676 ... 4.803917 \n",
"comt 13.498619 6.406145 5.750445 18.367493 ... 4.375763 \n",
"comt 5.053625 7.954544 7.999264 1.751066 ... 6.104701 \n",
"\n",
" dlscore_02 dlscore_03 dlscore_04 dlscore_05 dlscore_06 dlscore_07 \\\n",
"comt 6.189641 5.945060 5.594811 5.756034 6.155606 5.419800 \n",
"comt 5.334916 4.293091 5.267609 5.085522 6.309859 4.435493 \n",
"comt 5.444674 3.924729 5.376449 5.361448 6.503062 4.328010 \n",
"comt 4.890980 3.990071 5.115571 4.450336 5.754834 3.881885 \n",
"comt 5.777875 5.372001 5.696952 5.494438 6.711897 4.867440 \n",
"\n",
" dlscore_08 dlscore_09 dlscore_10 \n",
"comt 6.898730 5.973953 6.864763 \n",
"comt 5.776924 4.851845 5.440871 \n",
"comt 6.359634 4.692744 5.696967 \n",
"comt 5.338937 4.852479 5.326446 \n",
"comt 6.407499 5.344415 6.223027 \n",
"\n",
"[5 rows x 30 columns]"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"\n",
"# Generate column labels for the csv\n",
"nn_labels = [\"nnscore_\" + str(i).zfill(2) for i in range(1, 21)]\n",
"dl_labels = [\"dlscore_\" + str(j).zfill(2) for j in range(1, 11)]\n",
"header = nn_labels + dl_labels\n",
"\n",
"# Read csv files\n",
"all_actives = pd.read_csv(\"output_actives_comt/comt_predicted_values.csv\", names=header)\n",
"all_decoys = pd.read_csv(\"output_decoys_comt/comt_predicted_values.csv\", names=header)\n",
"\n",
"# Display scores\n",
"all_actives.head(5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculate the average scores for each tool"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>nnscore_avg</th>\n",
" <th>dlscore_avg</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>6.851305</td>\n",
" <td>6.107865</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>5.433614</td>\n",
" <td>5.173504</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>3.999191</td>\n",
" <td>5.249163</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>5.908439</td>\n",
" <td>4.797730</td>\n",
" </tr>\n",
" <tr>\n",
" <th>comt</th>\n",
" <td>4.857370</td>\n",
" <td>5.800025</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" nnscore_avg dlscore_avg\n",
"comt 6.851305 6.107865\n",
"comt 5.433614 5.173504\n",
"comt 3.999191 5.249163\n",
"comt 5.908439 4.797730\n",
"comt 4.857370 5.800025"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Calculate row averages for each subset\n",
"all_actives['nnscore_avg'] = all_actives[nn_labels].mean(axis=1)\n",
"all_actives['dlscore_avg'] = all_actives[dl_labels].mean(axis=1)\n",
"all_decoys['nnscore_avg'] = all_decoys[nn_labels].mean(axis=1)\n",
"all_decoys['dlscore_avg'] = all_decoys[dl_labels].mean(axis=1)\n",
"\n",
"# Display results\n",
"all_actives[['nnscore_avg', 'dlscore_avg']].head(5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generating ROCs"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.metrics import roc_curve, auc\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def plot_roc (actives_list, score_list, tool, receptor):\n",
" \"\"\" Calculates and plots and ROC and AUC.\n",
" Parameters:\n",
" actives_list - binary array of active/decoy status.\n",
" score_list - array of experimental scores.\n",
" tool- a string name of the tool used.\n",
" receptor - a string name of the protein used.\n",
" \"\"\"\n",
" # Get arrays for true positive rates and false positive rates\n",
" fpr, tpr, _ = roc_curve(actives_list, score_list)\n",
" # Area under curve\n",
" roc_auc = auc(fpr, tpr)\n",
" # Plot figure\n",
" plt.figure()\n",
" lw = 2\n",
" plt.plot(fpr, tpr, color='darkorange',\n",
" lw=lw, label='ROC curve (area = %0.3f)' % roc_auc)\n",
" plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\n",
" plt.xlim([0.0, 1.0])\n",
" plt.ylim([0.0, 1.05])\n",
" plt.xlabel('False Positive Rate')\n",
" plt.ylabel('True Positive Rate')\n",
" plt.title('Receiver operating characteristic for ' + tool + \" in \" + receptor)\n",
" plt.legend(loc=\"lower right\")\n",
" plt.savefig(\"figures/\" + tool + \"_\" + receptor + \"_ROC.svg\", format=\"svg\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate a binary array of true values."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1 1 1 ..., 0 0 0]\n"
]
}
],
"source": [
"num_actives = len(all_actives)\n",
"num_decoys = len(all_decoys)\n",
"true_values = np.repeat([1, 0], [num_actives, num_decoys])\n",
"print(true_values)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate arrays of combined active and decoy scores."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(4011,)\n",
"(4011,)\n",
"(4011,)\n"
]
}
],
"source": [
"nnscore_values = np.concatenate((all_actives['nnscore_avg'], all_decoys['nnscore_avg']))\n",
"dlscore_values = np.concatenate((all_actives['dlscore_avg'], all_decoys['dlscore_avg']))\n",
"\n",
"# Make sure that the arrays are the same size\n",
"print(true_values.shape)\n",
"print(nnscore_values.shape)\n",
"print(dlscore_values.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the ROC curves:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mkdir: cannot create directory ‘figures’: File exists\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3Xd4VNXWwOHfSgJJ6ASEiyBSpQookSJKlSJgL2DBcrEgigVFVLChXBURFGk25FOuomLjAoKAAhYQghSlCAgIAaQTakLK+v44J2EIKZOYyZSs93l4yJy65syZWWfvfc7eoqoYY4wx2QnzdwDGGGMCmyUKY4wxObJEYYwxJkeWKIwxxuTIEoUxxpgcWaIwxhiTI0sUBUBEbhGRb/0dh7+JSHUROSoi4YW4zxoioiISUVj79CURWSMi7fOxXr7OQRGJFpH/iUiCiHyW1/VDnYhMFJGn/R2Hv0moPUchIluBykAqcBSYDTygqkf9GVcoco/1Xao6z48x1AC2AMVUNcVfcbixKFBXVTf5eD81KKD3LCJ9gAHAxf4+fibv3IuKKapazZf7CdUSxRWqWgpoBlwAPOnnePLFn1fJoXKFnhdF9HifC2zIT5IItnMk2OINKKoaUv+ArcBlHq9HADM9XkcCI4FtwG5gIhDtMf8qYCVwGPgT6OZOLwu8B+wCdgAvAuHuvDuAH92/JwIjM8X0NTDQ/fts4HNgL85V4YMeyz0HTAOmuPu/K4v3Vxb4wF3/L2AoEOYRx0/Am0ACsB7olGndnN7DT8Bo4IA7rzbwHbAf2Af8FyjnLv8hkAacwCm5PQ7UABSIcJdZALzgbvcI8C1Q0SOe29z3sB94OvNnl+l9RwOvucsnAD+609L3ebv7me4Dhnis1wJYDBxy3/dYoLjHfAXuBzYCW9xpbwDb3c9gOXCpx/LhwFPuuXHEnX8OsMjd1jH3ePRyl++Jcz4dAn4GmmQ6VwcDq4EkIMLzGLixx7lx7AZGudO3ufs66v5rjcc56C7TCJjrfpa7gaeyOKbPAyeBZHc7fXEuHoe6x3kPzrlW1l0+/Vj3dWNYlMU22wPxwKPu+ruAOz3mTwbGATPd4/cLUNudJzjn3x73M14NNM7p83fnXQmscY/xAqBBLsc42+9gFu9nMvCiN+8ti3VjgPeBncBB4CuPeXcDm9zPZzpwdqZzsj/OOXkE5ztUG+c8Pgx8ChQHSuJ8/9I8zoWzs4vnH/2uFuaPeGH84/QvWjXgN+ANj/mvux9MDFAa+B/wkscXMwHo7H5hqgL13XlfAW+5H04lYClwrzvvDk4lirY4PzLp1Xrl3Q/zbHeby4Fn3A+6FrAZ6Oou+xzOl/Zqd9noLN7fBziJpzTOF3cD0NcjjhTgEaAY0Mt9PzFevocUnGqICJwvZh33WEQCZ+H8GL6e1bHO9EPimSj+BM5zt7cAeNmd19A9sS9xj8VI971nlyjGuetXxfmxvtiNK32f77j7aIrzg9DAXa850Mp9TzWAdcDDmb6Uc3HOh/QfnluBCu46jwJ/A1HuvEE451Q9nB+2pkAFj23V8dj2hTg/KC3dmG93j1mkx/FbiZNoojMfU5wfhj7u36WAVlkd5yzOwdI4P2KPAlHu65bZHNfncKou0l//G+cHrJa7zy+ADzPt9wOccyir87M9znk0DOcc7A4cB8q78yfj/Di2cI/vf4Gp7ryuON+Pcu6xbQBUyeXzPw8nOXd29/e4G3/xrI4xuXwHs3g/kzk9UWT73rJYdybwCc5vQDGgnTu9I84FzYXue3gTj6TrHuPpQBmchJ8EzHdjLQusBW73iCne57+rhf1D7vM35JwYR3EysboHOP0qWNyTqrbH8q05dSX5FjA6i21Wdj8sz5LHTcD3WXxJBedqq637+m7gO/fvlsC2TNt+Enjf40t7xlWax7LhbhwNPabdCyzwiGMnbpJypy0F+nj5HrZlt293mauBFZmOdW6JYqjH/P7AbPfvZ4CPPeaVwLm6PSNR4Hy5TwBNs5iXvs9qmd5z72zew8PAlx6vFeiYy/s+mL5v4A/gqmyWy5woJgAvZFrmD079YGwF/p3F+ZueKBbhXPVXzLTMacc5i3PwJs/PKZf39hynJ4r5QH+P1/VwEniEx35r5bC99u5n5RnbHk4lucnAux7zugPr3b874lz4tMItJXvx+T8NfJpp2R1A+6yOMbl8B7PY/mROTxTZvrdM61XBudI/I4nglOpHeLwu5R7jGh7nURuP+cuBwR6vX8O9YKOQEkWo1tldrarzRKQd8BFQEadYehbOD9JyEUlfVnB+gMG56piVxfbOxbki2OWxXhhOyeE0qqoiMhXny7oIuBmnKil9O2eLyCGPVcKBHzxen7FNDxVxroL+8pj2F85VVrod6p5BHvPP9vI9nLZvEakEjAEuxbkqDcP50cyLvz3+Po7zpcCNKWN/qnpcRPZns42KOFfGf+Z1PyJyHjAKiMX57CNwvnieMr/vR4G73BgV58quojv7nFzi8HQucLuIDPCYVtzdbpb7zqQvztXrehHZAjyvqjO82G9eYszsbM48vyJwLjTS5RQzwH49vc3D83OHbD4rVf1ORMbilB6qi8iXwGM4n312n/9p8apqmohs5/TvhGe83nwH/8l7S3cOcEBVs/q+nA386hHzUffcr4qT2MCpLkx3IovX//Iy3gIRqo3ZAKjqQpwrgpHupH04B7mRqpZz/5VVp+EbnBOqdhab2o5zNV7RY70yqtoom11/DFwvIufiXMF87rGdLR7bKKeqpVW1u2fYObylfThXHud6TKuOcwWVrqp4ZAJ3/k4v30Pmfb/kTmuiqmVwqmQkh+XzYhdO1SDg3KaJU92TlX1AIll/NrmZgNNWU9d9D09x+nsAj/chIpfi1GnfiHM1WA6n+i59nezOkaxsB4Zn+rxLqOrHWe07M1XdqKo34VQTvgJME5GSOa2Tjxgz28mZ51cKp/9Q/ZPPPUeqOkZVm+NUuZyHU9WX0+d/WrzuuX8Op38nPOP15jtYELYDMSJSzouYS+Kc+zuyWDY3PvssPIV0onC9DnQWkWaqmoZTlz3avVpGRKqKSFd32feAO0Wkk4iEufPqq+ounIbY10SkjDuvtltiOYOqrsBpKHsXmKOq6VcvS4HDIjLYvX89XEQai8hF3rwRVU3FacgaLiKl3UQ0kFMlFnB+VB4UkWIicgNOPe+svL4HV2mcarxDIlIV50vraTdOvWl+TAOuEJGLRaQ4ThVL5h9wwLlKBCYBo0TkbPe4tRaRSC/2UxqnAfCoiNQH7vNi+RSczy9CRJ7BKVGkexd4QUTqiqOJiKQnuMzH4x2gn4i0dJctKSI9RKS0F3EjIreKyFnu+08/h1Ld2NLI/tjPAP4lIg+LSKR7rrT0Zp84FzmPiEhNESkF/Af4RAvh1lkRucg9VsVwqogTgdRcPv9PgR7ud7YYTrtMEs6NA1n5R99Bb7nft2+A8SJS3v0+tnVnf4TzO9PMfQ//AX5R1a352NVuoIKIlC2QwLMR8olCVffiNL6lPzQzGKexa4mIHAbm4dTDoqpLgTtx7rxIABZyKvPfhlNtsBan+mUaTj1kdj4GLsM5KdJjSQWuwLltdwvOldK7OA1U3hqA8yXajHPnx0c4X6J0vwB13W0PB65X1fQqnby+h+dxGtwScBrmvsg0/yVgqIgcEpHH8vAeUNU17nuZilO6OIJT35uUzSqP4TQiL8NpDH0F787fx3Cq/47g/HB/ksvyc3C+4BtwqjQSOb3qYhTOj9O3OAnoPZxGUnDq+//PPR43qmocThvVWJzjvQmnLcFb3YA1InIU506s3qqaqKrHcT7bn9x9tfJcSVWP4DTuXoFTzbMR6ODlPifh3NG2COccTcT5nApDGZzP6CCn7oZLrw3I8vNX1T9wSrpv4pzzV+DcHn8yqx0U0HfQW31wagDW45zbD7sxzMf5Pfoc59yvDfTOzw5UdT3Ob81m91w4O7d18iPkHrgrykTkDpxbai/xdyx55V69HsKpItri73iMMaeEfInCBC4RuUJESrh1tCNxrhi3+jcqY0xmliiMP12F07C3E6e6rLdaEdeYgGNVT8YYY3JkJQpjjDE5CroH7ipWrKg1atTwdxjGGBNUli9fvk9Vz8rPukGXKGrUqEFcXJy/wzDGmKAiIn/lvlTWrOrJGGNMjixRGGOMyZElCmOMMTmyRGGMMSZHliiMMcbkyBKFMcaYHPksUYjIJBHZIyK/ZzNfRGSMiGwSkdUicqGvYjHGGJN/vixRTMbpJjk7l+P071MXuAdngBljjDEF7OTJ1H+0vs8euFPVRSJSI4dFrgI+cDuBWyIi5USkijvghzHGmH/iix6wZRaD/teZFTtzGnYmd/5so6jK6QPCxHP6OLcZROQeEYkTkbi9e/cWSnDGGBPUtswCoPG/9vDD5ur/aFP+7MIjq2Evs+zKVlXfBt4GiI2Nte5ujcmOexVpiq61f5/FrzuqcGtz5/VtH6+g3csJ1Kz5Yr636c9EEY8zCHq6ajjjEhhj8suSRJF1/GQxXpzXllcXXEx4mNLq3HjqXNQKEaFGjXL/aNv+TBTTgQdEZCrQEkiw9gljCsijVvAuSr75ZiP33z+LLVsOAdD37uZUePwpKB+dy5re8VmiEJGPgfZARRGJB54FigGo6kRgFtAdZ8D548CdvorFmIBhVUOmAO3YcZiHH57DtGlrAWjSpDITJ/agdetzclkzb3x519NNucxX4H5f7d+YgFQYSaJmd9/vwwSE+++fxddf/0GJEsUYNqw9Dz3UioiIgr9HKejGozCmUPj6yt+qhkw+paSkZSSDV165jGLFwnnttS5Ur17WZ/u0LjyMyYovk4Rd8Zt8SEhIZMCAWfTo8RFOhQzUq1eRzz67wadJAqxEYcyZvuhx6m+78jd+pqp89tlaHn54Nrt2HSU8XFi58m8uuOCfPUSXF5YojMksvTRhV/7Gz/788wAPPPANs2dvAqB162pMnNiTJk0qF2ocliiMyc61M/0dgSnCRo78maef/p7ExBTKlYvilVcu4667LiQsLKtnlX3LEoUpWuz2VBMkjh9PJjExhT59mjByZBcqVSrpt1gsUZiixdskYdVOppDt3XuMP/7YzyWXOP0yDR7chvbta9C27bl+jswShQlF3pQarJHaBIi0NGXSpBU8/vhcIiLCWL/+AWJioomMjAiIJAGWKEwoyi1JWGnBBIjff99Dv34z+OknpyPtzp1rcfx4MjExBdP1RkGxRGFCl5UaTIA6duwkw4YtZNSoJaSkpFG5cklef70bvXo1QqTwG6tzY4nCBB9rkDZB7vrrP2P27E2IQP/+sQwf3oly5aL8HVa2LFGY4ONNkrDqJRPABg9uw+7dR5kwoQctW1bzdzi5skRhAl92JQirWjJBICUljTff/IWtWw/xxhuXA9C+fQ3i4u7xyzMR+WGJwgS+rJKElRhMEFi6dAf33juDlSv/BuCee5rTqFElgKBJEmCJwgQTK0GYIHHoUCJPPTWfiRPjUIVzzy3L2LHdM5JEsLFEYQKXNVqbIDR16u88/PBsdu8+RkREGI8+2pqnn25LyZLF/R1avlmiMIHLM0lYVZMJEt9++ye7dx+jTZtzmDChB+efX7gd+PmCJQoTOKzR2gShpKQUduw4Qq1a5QEYMaIzl15andtvbxZU7RA5sYGLTOCwRmsTZL77bgtNmkykR4+POHkyFYCKFUtw550XhEySACtRmEBkJQgT4HbvPspjj81lypTVANSvX5H4+MMZpYpQY4nC+Jc1WJsgkpamvPPOcp54Yj6HDiUSFRXB0KGXMmhQG4oXD/d3eD5jicL4V+YkYVVNJoBdc80nTJ/+BwBdu9Zm3Lju1K4d4+eofM8ShfEfG5vaBJlrr63P0qU7eOONbtxwQ8OA7MDPFyxRGP+xsalNgJs+/Q/i4w/Tv/9FANx2W1OuvbYBpUtH+jmywmWJwhSOnNoibGxqE2C2bUvgwQe/4euv/yAyMpxu3epQq1Z5RKTIJQmwRGEKS3ZJwkoTJoAkJ6cyZswvPPvsAo4dS6Z06eK8+GJHzj23rL9D8ytLFKZwWVuECVBLlsRz770zWL16NwA33NCQ0aO7UrVqGT9H5n+WKEzBsVtdTRB7+unvWb16NzVrlmPs2O50717X3yEFDEsUpuDYWNUmiKgqR46cpEwZp81h7NjL+eCDVQwZ0pYSJYr5ObrAYonC/DNZlSKseskEuD/+2Ef//rMQgblz+yAi1KtXkeHDO/k7tIBkicL8M/bAnAkiiYkpvPTSD7z88k+cPJlKhQrRbN16iJo1Q7PrjYJiicIUDCtFmAA3d+6f9O8/i02bDgDw7383Y8SIzlSoUMLPkQU+n/YeKyLdROQPEdkkIk9kMb+6iHwvIitEZLWI2OVosPiiB7xWNJ5KNcFNVfn3v7+mS5cpbNp0gIYNz2LRojt4772rLEl4yWclChEJB8YBnYF4YJmITFfVtR6LDQU+VdUJItIQmAXU8FVMpgDZoEImSIgINWqUIzo6gmeeacfAga1DugM/X/Bl1VMLYJOqbgYQkanAVYBnolAg/SblssBOH8Zj/ilruDZBYuXKv9m16wiXX+7c4jp4cBv69GlibRH55Muqp6rAdo/X8e40T88Bt4pIPE5pYkBWGxKRe0QkTkTi9u7d64tYjTes4doEuCNHkhg4cA7Nm7/N7bd/xYEDJwCIjIywJPEP+LJEkVUFdubLz5uAyar6moi0Bj4UkcaqmnbaSqpvA28DxMbG2iWsv1kpwgQYVeWrr9bz4IOziY8/TFiYcPPN51OsmA3iWRB8mSjigXM8XlfjzKqlvkA3AFVdLCJRQEVgjw/jMt6yJ61NEPjrr0M88MA3zJixAYDY2LN5662eXHhhFT9HFjp8mW6XAXVFpKaIFAd6A9MzLbMN6AQgIg2AKMDqlgKFjWFtApyqct11nzJjxgbKlIlk7NjLWbKkryWJAuazEoWqpojIA8AcIByYpKprRGQYEKeq04FHgXdE5BGcaqk7VNXqNfwluxKEVTWZAJOWpoSFCSLCyJFdmDgxjtGju1KlSml/hxaSJNh+l2NjYzUuLs7fYYSmrJ6LqNndxoswAWP//uM88cQ8AN5550o/RxNcRGS5qsbmZ117Mts4bFhSE8BUlQ8+WMVjj81l377jFC8ezrPPtqdaNesCvDBYojAOG5bUBKh16/Zy330zWbjwLwDat6/BhAk9LEkUIksU5nRWzWQChKryzDPf88orP5GcnEbFiiV47bUu9OnTBBHrPqYwWaIwp1c7GRMgRIQdO46QnJzG3XdfyMsvX0ZMTLS/wyqSLFEYq3YyAWPnziPs23ecJk0qAzBiRGf69r2ANm2q+zmyos0eWyzKMvcAa9VOxk9SU9MYO3YpDRqMo3fvaZw8mQpAxYolLEkEACtRFGXWA6wJAL/+uot7751BXJzTcUPbtudy+HASFStaF+CBwqtE4T5ZXV1VN/k4HuMPdjus8YPDh5N4+unvGDt2GWlpSrVqZRgzphtXX13fGqsDTK6JQkR6AKOA4kBNEWkGPKuq1/g6OOND1oBt/EhVadv2fVat2k14uDBwYCuee649pUtH+js0kwVv2iiGAS2BQwCquhKo48ugTCGwBmzjRyLCI4+0okWLqsTF3cNrr3W1JBHAvKl6SlbVQ5mKglZXEYyy6svJGrBNITh5MpVRoxYTHi4MGtQGgNtua8qttzYhPNzuqQl03iSKdSJyIxAmIjWBh4Alvg3L+IQNPGT84Icf/qJfv5msXbuXyMhwbrutKZUrl0JECA+3tohg4E2ieAB4BkgDvsDpDfZJXwZlfMwar00h2LfvOI8/Ppf3318JQN26MYwf34PKlUv5OTKTV94kiq6qOhgYnD5BRK7FSRomGNgARKYQqSqTJ69k0KC57N9/guLFw3nyyUt44olLiIqyO/KDkTeVg0OzmDakoAMxPmTPS5hCNmXKb+zff4KOHWuyenU/nnuuvSWJIJbtJyciXXGGKa0qIqM8ZpXBqYYygS5zScKqnIyPHD+eTEJCIlWqlEZEGD++O8uW7eSWW863ZyJCQE4pfg/wO5AIrPGYfgR4wpdBmQJiJQlTCL75ZiP33z+LWrXKM3duH0SEevUqUq9eRX+HZgpItolCVVcAK0Tkv6qaWIgxmYJmJQnjAzt2HObhh+cwbdpaAEqXjmT//hPW9UYI8qbSsKqIDAcaAlHpE1X1PJ9FZYwJWKmpaYwbt4yhQ7/jyJGTlCxZjGHDOvDggy2JiLBnIkKRN4liMvAiMBK4HLgTa6MwpkhKS1PatZvMTz9tB+Dqq+vzxhvdqF69rJ8jM77kTaIooapzRGSkqv4JDBWRH3wdmMknuxXW+FBYmNClS222bUtg7NjuXHllPX+HZAqBN4kiSZzbFv4UkX7ADqCSb8My+WZPX5sCpKp8+ukaIiLCuO66hgAMHtyGgQNbU6pUcT9HZwqLN4niEaAU8CAwHCgL/NuXQZl8sFthTQH7888D9O8/i2+//ZOzzipBx441KV8+msjICCKt/74iJddEoaq/uH8eAfoAiEg1XwZl8sFuhTUFJCkphVdf/Znhw38gMTGF8uWjGD68I2XLRuW+sglJOSYKEbkIqAr8qKr7RKQRTlceHQFLFoHIShLmH1iwYCv33TeT9ev3AdCnTxNGjuxCpUol/RyZ8ads72UTkZeA/wK3ALNFZAjwPbAKsFtjA4kNQmQKQGpqGv37O0miXr0KfPfdbXzwwTWWJEyOJYqrgKaqekJEYoCd7us/Cic04zUbhMjkU1qakpiYQokSxQgPD2PChB4sWvQXjz/ehshI65vJOHI6ExJV9QSAqh4QkfWWJAJM5gZsG4TI5MFvv+2mX7+Z1K9fgffeuwqAdu1q0K5dDf8GZgJOTomiloikdyUuQA2P16jqtT6NzOTOGrBNPhw7dpJhwxYyatQSUlLS2LLlIAcPnqB8+Wh/h2YCVE6J4rpMr8f6MhDzD1gDtvHS//73Bw888A3btiUgAv37xzJ8eCfKlbM7mkz2cuoUcH5hBmLyyBqwTR6kpKTRq9c0vvhiHQDNmv2Lt97qSYsWVf0cmQkG1loVrKwB2+RBREQYZctGUqpUcV54oQMPPNDCOvAzXhNV31VbiEg34A0gHHhXVV/OYpkbgecABVap6s05bTM2Nlbj4uJ8EG0A8qbfJqt2Mtn45Zd4AFq2dB552r//OCdOpFCtWhl/hmX8RESWq2psftb1ukQhIpGqmpSH5cOBcUBnIB5YJiLTVXWtxzJ1gSeBNqp6UESsDylPuSUJK02YLBw6lMiTT87jrbeWU79+RVau7Efx4uFUqGDjRJj8yTVRiEgL4D2cPp6qi0hT4C5VHZDLqi2ATaq62d3OVJxnM9Z6LHM3ME5VDwKo6p68v4UiwEoNxguqyscf/87AgXPYvfsYERFhXHllPVJT03AK9cbkjzclijFAT+ArAFVdJSIdvFivKrDd43U80DLTMucBiMhPOGfyc6o624tthz5rrDZ5sHHjfvr3n8W8eZsBaNPmHCZO7EnjxlZIN/+cN4kiTFX/yjRAeqoX62U1onrmS+MIoC7QHqfvqB9EpLGqHjptQyL3APcAVK9e3YtdhwBrrDZeSk5OpWPHD4iPP0xMTDQjRlzGnXdeQFhYVl9BY/LOm0Sx3a1+UrfdYQCwwYv14oFzPF5Xw+kGJPMyS1Q1GdgiIn/gJI5lngup6tvA2+A0Znux79BhT1ubbKgqIkKxYuEMH96R77/fyogRl3HWWdY3kylY3twfdx8wEKgO7AZaudNyswyoKyI1RaQ40BuYnmmZr4AOACJSEacqarN3oRtTNO3efZQ+fb7kxRcXZUy77bamvP/+VZYkjE94U6JIUdXeed2wqqaIyAPAHJz2h0mqukZEhgFxqjrdnddFRNbiVGcNUtX9ed2XMUVBWpryzjvLeeKJ+Rw6lEi5clE8/HArSpe2UYSMb3mTKJa5VUKfAF+o6hFvN66qs4BZmaY94/G34pRWBnq7TWOKolWr/qZfv5ksWeI8G9GtWx3GjetuScIUCm9GuKstIhfjVB09LyIrgamqOtXn0RlTxCUnp/Lkk/N5/fUlpKYqVaqU4o03unH99Q3JdIOJMT7j1TP8qvqzqj4IXAgcxhnQyBS0L3rAa+L8Mwan640VK/4mLU0ZMKAF69bdzw03NLIkYQqVNw/clcJ5UK430AD4GrjYx3EVTZmfxLZbY4ukbdsSSE1No2bN8ogIEyf2ICEhidjYs/0dmimivGmj+B34HzBCVX/wcTxFU+Y+nexJ7CIpOTmVN974hWefXUDr1tWYO7cPIkLduhX8HZop4rxJFLVUNc3nkRRlNgBRkbd48Xb69ZvJ6tW7AYiJieb48WRKlizu58iMySFRiMhrqvoo8LmInHGJayPc+YCVJIqcgwdP8MQT83j77V8BqFmzHOPGdefyy+v6OTJjTsmpRPGJ+7+NbOdL1qdTkZWUlEKzZm+xbVsCxYqFMWjQxQwZ0pYSJYr5OzRjTpPTCHdL3T8bqOppycJ9kM5GwCsI1qdTkRUZGUHfvhcwf/4WJkzoQcOGZ/k7JGOylOvARSLyq6pemGnaClW9wKeRZSOkBi7ybMS2aqeQl5iYwksv/UC9ehW5+ebzAWeI0vBwsdtdjc/5ZOAiEemFc0tsTRH5wmNWaeBQ1muZPLHSRJExd+6f9O8/i02bDlCpUkmuuaY+0dHFbDhSExRyaqNYCuzH6fV1nMf0I8AKXwZV5FgPsSHr77+PMnDgHD7++HcAGjU6i4kTexIdbe0QJnjk1EaxBdgCzCu8cIoQa8QOaampabz11nKeemo+CQlJREdH8Oyz7XjkkdYUL26jzZngklPV00JVbSciBzl9wCHB6c8vxufRhTKrdgppqanKm28uJSEhie7d6zJ27OXUrFne32EZky85VT2lD3dasTACCWmZn7z2ZNVOIePIkSRSU5Vy5aIoXjycd965gt27j3LttQ2ssdoEtWxb0jyexj4HCFfVVKA1cC9go6PkRXZJwkoTIUFV+eKLdTRoMI5HH52TMf2SS6pz3XXWy6sJft504fEVcJGI1AY+AGYCHwFiMnTlAAAgAElEQVQ9fRlYSLA+nELe1q2HGDDgG2bMcEYH/v33vSQmphAV5c1Xy5jg4M29eWnumNbXAq+r6gCgqm/DChHWh1PISk5O5ZVXfqRhw3HMmLGBMmUiGTv2cn7++d+WJEzI8WooVBG5AegDXO1Os3v78sJKEiHl+PFkWrV6l99+2wNA796NGTWqC1WqlPZzZMb4hjeJ4t9Af5xuxjeLSE3gY9+GFeRyarw2Qa9EiWLExp7N8ePJjB/fgy5davs7JGN8ypuhUH8XkQeBOiJSH9ikqsN9H1oQsyqnkKKqfPDBKmrXjuGSS6oDMHp0V4oXD7cH50yR4M0Id5cCHwI7cJ6h+JeI9FHVn3wdXNCzKqegt27dXu67byYLF/5FgwYVWbmyH8WLh1O2bJS/QzOm0HhT9TQa6K6qawFEpAFO4shX51LGBIMTJ5IZPvwHRoz4ieTkNM46qwRPPnkJxYpZ30ym6PEmURRPTxIAqrpORGzYLROyZs/exP33z2Lz5oMA3H33hbz88mXExET7OTJj/MObRPGriLyFU4oAuAXrFDBr1ogd9I4ePUmfPl+yb99xGjeuxMSJPWjTprq/wzLGr7xJFP2AB4HHcdooFgFv+jKooGWN2EEpNTWNtDSlWLFwSpUqzhtvdCM+/jCPPNKKYsWsAz9jckwUInI+UBv4UlVHFE5IQSC3koM1YgeN5ct3cu+9M7jqqno8/XQ7gIxBhYwxjmxb5kTkKZzuO24B5orIvwstqkCXU5KwkkRQOHw4iYce+oYWLd5l+fJdfPjhapKTU/0dljEBKacSxS1AE1U9JiJnAbOASYUTVpCwkkPQUVWmTVvLQw/NZteuo4SHCwMHtuL55ztYNZMx2cgpUSSp6jEAVd0rInZfoAlqR44k0avXNL75ZhMALVtWZeLEnjRr9i8/R2ZMYMspUdTyGCtbgNqeY2er6rU+jcyYAlaqVHGSklIpWzaSl1++jHvuaU5YmHUBbkxuckoU12V6PdaXgRjjC4sW/UWVKqWoW7cCIsKkSVcSFRVB5cql/B2aMUEjpzGz5xdmIMYUpH37jvP443N5//2VdOpUk7lz+yAinHtuOX+HZkzQsY7zTUhJS1MmT17JoEFzOXDgBMWLh3PppdVJTVUiIqyayZj88GkDtYh0E5E/RGSTiDyRw3LXi4iKiPUfZfJtzZo9tG8/mb59p3PgwAk6darJb7/dx7PPticiwu7FMCa/vC5RiEikqiblYflwYBzQGYgHlonIdM9+o9zlSuM8+f2Lt9s2JrOEhERatXqPo0dPUqlSSUaN6sLNN59v41UbUwByvcwSkRYi8huw0X3dVES86cKjBc7YFZtV9SQwFbgqi+VeAEYAid6HbYxD1XmWpWzZKAYPbkO/fs1Zv/5+brmliSUJYwqIN+XxMUBPYD+Aqq4COnixXlVgu8freDKNtS0iFwDnqOqMnDYkIveISJyIxO3du9eLXZtQt2PHYa6//lOmTFmdMW3IkEuZMKEn5ctbL6/GFCRvEkWYqv6VaZo3fR1kdTmX8Siz+wDfaODR3Dakqm+raqyqxp511lle7NpHvugBr9lVqj+lpKTxxhtLqF9/HJ9/vo5nn11AamoagJUgjPERb9oototIC0DddocBwAYv1osHzvF4XQ3Y6fG6NNAYWOB+wf8FTBeRK1U1zpvgC531DutXy5btoF+/mfz66y4Arr66PmPGdCM83BqqjfElbxLFfTjVT9WB3cA8d1pulgF1RaQmzjCqvYGb02eqagJQMf21iCwAHgvYJOHJ+ngqVMeOnWTw4HmMH78MVahevSxvvnk5V15Zz9+hGVMk5JooVHUPzo98nqhqiog8AMwBwoFJqrpGRIYBcao6Pc/R+tMXPfwdQZEVERHGvHmbCQsTBg5szbPPtqNkSRtk0ZjCkmuiEJF38GhbSKeq9+S2rqrOwul11nPaM9ks2z637flVerWTVTkVij//PEC5clFUqFCCyMgIPvzwGqKiIjj//Mr+Ds2YIsebyt15wHz3309AJcDr5ylCgmdp4tqZ/oujCEhKSuHFFxfRuPEEBg+elzH9oouqWpIwxk+8qXr6xPO1iHwIzPVZRIHIShOFYsGCrdx330zWr98HOHc4paamWWO1MX6Wn76eagLnFnQgQcFKEz6xZ88xBg2aywcfrAKgXr0KTJjQgw4davo5MmMMeNdGcZBTbRRhwAEg236bgl5u42GbArVv33EaNBjHgQMniIwMZ8iQS3n88TZERlp/lcYEihy/jeI84NAU5/ZWgDRN7zMhVGWXJKzayScqVizBVVfVIz7+MOPH96BOnRh/h2SMySTHRKGqKiJfqmrzwgooYNizEj5x7NhJhg1bSI8e59G2rVODOX58DyIjw+3JamMClDethEtF5EKfR2JC3v/+9wcNG45nxIif6d9/JmlpTjKOioqwJGFMAMu2RCEiEaqaAlwC3C0ifwLHcPpwUlW15GG8sn17Ag89NJsvv1wPwAUX/Iu33upp41UbEyRyqnpaClwIXF1IsZgQk5KSxpgxv/DMM99z7FgypUoV58UXO3D//S1sICFjgkhOiUIAVPXPQorFhJjDh5N46aUfOXYsmeuua8Drr3ejWrUy/g7LGJNHOSWKs0RkYHYzVXWUD+LxL+vP6R87dCiR6OgIIiMjiImJ5q23ehIZGU6PHuf5OzRjTD7lVP4PB0rhdAee1b/QY09g55uq8tFHv1Gv3lhGjPgpY/q11zawJGFMkMupRLFLVYcVWiT+ktUDdvYEdp5s2LCf/v1nMn/+FgAWLdqGqtqdTMaEiFzbKEJe5iRhpQmvJSam8MorP/Kf//zIyZOpxMRE8+qrnbnjjmaWJIwJITklik6FFkUgsAfs8uTvv4/Stu37bNx4AIA77mjGq692pmLFEn6OzBhT0LJNFKp6oDAD8QtrvM63ypVLcs45ZYmICGPChB60a1fD3yEZY3ykaPe8Zo3XXktLU955ZzkdOtTkvPMqICJ89NG1lC8fTfHi4f4OzxjjQ/bUE1jjdS5WrfqbNm0m0a/fTPr3n0l6v5CVK5eyJGFMEVC0SxQmR0ePnuS55xbw+utLSE1Vzj67NP36xfo7LGNMIbNEYbL01VfrGTDgG+LjDxMWJgwY0IIXX+xImTKR/g7NGFPIim6isIbsbO3YcZjevaeRlJRK8+ZVmDixJ7GxZ/s7LGOMnxTdRGEN2adJTk4lIiIMEaFq1TIMH96R4sXD6d//Ihuz2pgiruj9AnzRA17zeBjMGrL5+eftNG/+NlOmrM6Y9uijFzNgQEtLEsaYIpgoPJ/ELuKliQMHTnDvvf+jTZtJ/PbbHsaPjyPUR7o1xuRd0a16KsJPYqsqU6as5tFHv2Xv3uMUKxbG44+3YciQS63rDWPMGYpOosiq878iaPfuo9x00+d8//1WANq1O5cJE3rQoMFZ/g3MGBOwik6isConAMqVi2LXrqNUrFiCkSM7c9ttTa0UYYzJUdFJFOmKYJXT3Ll/cuGFVahQoQSRkRF89tkNVKlSigoVrAM/Y0zuil5jdhGya9cRbrrpc7p0mcLgwfMypjduXMmShDHGa0WvRFEEpKam8dZby3nyyfkcPpxEdHQE9epVsMGEjDH5UjQSRRF6CvvXX3fRr98Mli3bCUCPHnUZO7Y7NWqU83NkxphgVTQSRRF5Cnvr1kO0aPEOqalK1aqlGTPmcq65pr6VIowx/4hPE4WIdAPeAMKBd1X15UzzBwJ3ASnAXuDfqvqXzwIK8aewa9Qox513NqN06Uief749pUtbB37GmH/OZ43ZIhIOjAMuBxoCN4lIw0yLrQBiVbUJMA0Y4at4QtHWrYe44oqPWbhwa8a0t9++glGjulqSMMYUGF+WKFoAm1R1M4CITAWuAtamL6Cq33ssvwS41YfxhIzk5FRGjVrM888v5MSJFPbtO87ixX0BrJrJGFPgfJkoqgLbPV7HAy1zWL4v8E1WM0TkHuAegOrVqxdUfEHpxx+30a/fDNas2QtA796NGTWqi5+jMsaEMl8miqwubbN82k1EbgVigXZZzVfVt4G3AWJjY4veE3PAwYMnGDRoLu+9twKA2rXLM358D7p0qe3nyIwxoc6XiSIeOMfjdTVgZ+aFROQyYAjQTlWTCjyKELk1Ni1N+frrPyhWLIwnnriEJ5+8hOjoYv4OyxhTBPgyUSwD6opITWAH0Bu42XMBEbkAeAvopqp7fBJFEN8au379PmrWLEdkZAQVKpTgv/+9lurVy1K/fkV/h2aMKUJ8dteTqqYADwBzgHXAp6q6RkSGiciV7mKvAqWAz0RkpYhM91U8wXRr7PHjyQwZMp8mTSYwYsRPGdO7dKltScIYU+h8+hyFqs4CZmWa9ozH35f5cv/BaPbsTfTvP5MtWw4BsG/fcT9HZIwp6orGk9lBYOfOIzz88Gw++8y5e/j88ysxcWJPLr74nFzWNMYY37JEEQA2bNhPbOzbHDlykhIlivHcc+14+OFWFCsW7u/QjDHGEkUgqFs3hosuqkrJksV4883LOfdc68DPGBM4LFH4weHDSTzzzPf0738R551XARFh+vTelCxZ3N+hGWPMGSxRFCJVZdq0tTz00Gx27TrK+vX7mD3b6bXEkoQxJlBZoigkmzcf5IEHZvHNN5sAaNWqGq+8Yjd9GWMCX+gmii96nHrYzo9Onkxl5MifeeGFRSQmplCuXBQvv9yJu+9uTliYdeBnjAl8oZsoPJOEH5/K3r49gWHDFpKUlMott5zPa691oXLlUn6Lxxhj8ip0E0W6Rwu/D8GDB09QrlwUIkLt2jG88UY36tSJoVOnWoUeizHG/FM+68KjKEpLUyZNWkGdOm8yZcrqjOn33htrScIYE7RCM1H4ocfYNWv20L79ZPr2nc6BAycyGq2NMSbYhWbVUyH2GHv8eDIvvLCQkSMXk5KSRqVKJRk9uis33dTY5/s2xpjCEJqJIp2Pe4zdsGE/XbtOYevWQ4hAv37N+c9/OlG+fLRP92uMMYUptBJFId8Se+65ZYmKiqBp08pMnNiTVq2qFdq+TeBITk4mPj6exMREf4diDFFRUVSrVo1ixQpuYLPQShQ+viU2JSWNiRPjuOmmxlSoUILIyAhmz76FqlXLEBERms09Jnfx8fGULl2aGjVqIGLPxhj/UVX2799PfHw8NWvWLLDthlaiSOeDW2KXLt1Bv34zWLHib1au/Jt333XGXrIO/ExiYqIlCRMQRIQKFSqwd+/eAt1uaCaKApSQkMiQId8xfvwyVKF69bJcdVU9f4dlAowlCRMofHEuWqLIhqryySdreOSROfz991EiIsIYOLAVzzzTzjrwM8YUKVaxno1Vq3Zz002f8/ffR7n44nP49dd7eOWVzpYkTEAKDw+nWbNmNG7cmCuuuIJDhw5lzFuzZg0dO3bkvPPOo27durzwwguonqqe/eabb4iNjaVBgwbUr1+fxx57zB9vIUcrVqzgrrvu8ncYOXrppZeoU6cO9erVY86cOVkuM3/+fC688EKaNWvGJZdcwqZNzvNWo0aNomHDhjRp0oROnTrx119/ZazTrVs3ypUrR8+ePU/bVu/evdm4caPv3pAnVQ2qf82bN9dsjcT5l08pKamnvX7kkdn6zjvLNTU1Ld/bNKFv7dq1/g5BS5YsmfH3bbfdpi+++KKqqh4/flxr1aqlc+bMUVXVY8eOabdu3XTs2LGqqvrbb79prVq1dN26daqqmpycrOPGjSvQ2JKTk//xNq6//npduXJloe4zL9asWaNNmjTRxMRE3bx5s9aqVUtTUlLOWK5u3boZ58u4ceP09ttvV1XV7777To8dO6aqquPHj9cbb7wxY5158+bp9OnTtUePHqdta8GCBXrXXXdlGU9W5yQQp/n83bWqJ9f332+hf/9ZvPVWT9q2PReAUaO6+jkqE3Re81FbRR5u0GjdujWrVztdyHz00Ue0adOGLl26AFCiRAnGjh1L+/btuf/++xkxYgRDhgyhfv36AERERNC/f/8ztnn06FEGDBhAXFwcIsKzzz7LddddR6lSpTh69CgA06ZNY8aMGUyePJk77riDmJgYVqxYQbNmzfjyyy9ZuXIl5co5N3/UqVOHn376ibCwMPr168e2bdsAeP3112nTps1p+z5y5AirV6+madOmACxdupSHH36YEydOEB0dzfvvv0+9evWYPHkyM2fOJDExkWPHjvHdd9/x6quv8umnn5KUlMQ111zD888/D8DVV1/N9u3bSUxM5KGHHuKee+7x+vhm5euvv6Z3795ERkZSs2ZN6tSpw9KlS2nduvVpy4kIhw8fBiAhIYGzzz4bgA4dOmQs06pVK6ZMmZLxulOnTixYsOCMfV566aXccccdpKSkEBHh25/yIp8o9uw5xqBBc/ngg1UAjBq1OCNRGBNsUlNTmT9/Pn379gWcaqfmzZuftkzt2rU5evQohw8f5vfff+fRRx/NdbsvvPACZcuW5bfffgPg4MGDua6zYcMG5s2bR3h4OGlpaXz55Zfceeed/PLLL9SoUYPKlStz880388gjj3DJJZewbds2unbtyrp1607bTlxcHI0bn+rpoH79+ixatIiIiAjmzZvHU089xeeffw7A4sWLWb16NTExMXz77bds3LiRpUuXoqpceeWVLFq0iLZt2zJp0iRiYmI4ceIEF110Eddddx0VKlQ4bb+PPPII33///Rnvq3fv3jzxxBOnTduxYwetWrXKeF2tWjV27Nhxxrrvvvsu3bt3Jzo6mjJlyrBkyZIzlnnvvfe4/PLLcz2+YWFh1KlTh1WrVp3xGRe0Ipso0tKU9977lcGD53HwYCKRkeEMHdqWQYMu9ndoJpj5obdigBMnTtCsWTO2bt1K8+bN6dy5M+BULWd3F0xe7o6ZN28eU6dOzXhdvnz5XNe54YYbCA8PB6BXr14MGzaMO++8k6lTp9KrV6+M7a5duzZjncOHD3PkyBFKly6dMW3Xrl2cddZZGa8TEhK4/fbb2bhxIyJCcnJyxrzOnTsTExMDwLfffsu3337LBRdcADiloo0bN9K2bVvGjBnDl19+CcD27dvZuHHjGYli9OjR3h0cOK3NJ11Wx3f06NHMmjWLli1b8uqrrzJw4EDefffdjPlTpkwhLi6OhQsXerXfSpUqsXPnTksUXsnjE9lbthzk1lu/5OeftwPQpUttxo3rTp06Mb6K0Bifio6OZuXKlSQkJNCzZ0/GjRvHgw8+SKNGjVi0aNFpy27evJlSpUpRunRpGjVqxPLlyzOqdbKTXcLxnJb5yfSSJUtm/N26dWs2bdrE3r17+eqrrxg6dCgAaWlpLF68mOjo7Lu9iY6OPm3bTz/9NB06dODLL79k69attG/fPst9qipPPvkk995772nbW7BgAfPmzWPx4sWUKFGC9u3bZ/lUfV5KFNWqVWP79u0Zr+Pj4zOqldLt3buXVatW0bJlS8BJnt26dcuYP2/ePIYPH87ChQuJjIzM9nh4SkxMzPHYFZTQuOspj09klykTyYYN+/nXv0oxdep1zJ59iyUJExLKli3LmDFjGDlyJMnJydxyyy38+OOPzJs3D3BKHg8++CCPP/44AIMGDeI///kPGzZsAJwf7lGjRp2x3S5dujB27NiM1+lVT5UrV2bdunUZVUvZERGuueYaBg4cSIMGDTKu3jNvd+XKlWes26BBg4y7g8ApUVStWhWAyZMnZ7vPrl27MmnSpIw2lB07drBnzx4SEhIoX748JUqUYP369VlW/4Bz9b9y5coz/mVOEgBXXnklU6dOJSkpiS1btrBx40ZatGhx2jLly5cnISEh41jPnTuXBg0aAM5dXffeey/Tp0+nUqVK2b6nzDZs2ECjRo28Xj6/QiNRpHtUs+0IcM6cTSQlpQBQoUIJpk/vzfr199OrV2N7WMqElAsuuICmTZsydepUoqOj+frrr3nxxRepV68e559/PhdddBEPPPAAAE2aNOH111/npptuokGDBjRu3Jhdu3adsc2hQ4dy8OBBGjduTNOmTTOutF9++WV69uxJx44dqVKlSo5x9erViylTpmRUOwGMGTOGuLg4mjRpQsOGDZk4ceIZ69WvX5+EhASOHDkCwOOPP86TTz5JmzZtSE1NzXZ/Xbp04eabb6Z169acf/75XH/99Rw5coRu3bqRkpJCkyZNePrpp09rW8ivRo0aceONN9KwYUO6devGuHHjMqrdunfvzs6dO4mIiOCdd97huuuuo2nTpnz44Ye8+uqrgJOwjx49yg033ECzZs248sorM7Z96aWXcsMNNzB//nyqVauWcevt7t27iY6OzvW4FwTJqm4tkMXGxmpcXNzpE9PvNMmifnj79gQefHA2X321nhde6MDQoW0LIUpTlKxbty7jytD4xujRoyldunTAP0tRmEaPHk2ZMmUyblzwlNU5KSLLVTU2P/sK/hJFNoMUpaSkMWrUYho0GMdXX62nVKnixMRY99/GBKP77rvP63r7oqJcuXLcfvvthbKv4G/MzmKQoiVL4unXbwarVu0G4LrrGvDGG92oWrWMPyI0xvxDUVFR9OnTx99hBJQ777yz0PYV/Ikinds28csv8Vx88XuoQo0a5Rg79nJ69DjPz8GZUJfTbajGFCZfNCeETqJwtWhRla5d63DBBf9i6NC2lChRcIN3GJOVqKgo9u/fT4UKFSxZGL9SdzyKqKioAt1u0CeKjXtjeGR6N0ZdsZ/zznO+qDNn3kxYmH1hTeGoVq0a8fHxBT4GgDH5kT7CXUEKzkTxRQ+SNs7h5e8u4aXv+pOUEkHUU/OZNu1GAEsSplAVK1asQEcTMybQ+PSuJxHpJiJ/iMgmETnjKRURiRSRT9z5v4hIDW+2O//bdTR57T6e+7YDSSkR3NlxLxMn9sx9RWOMMXnmsxKFiIQD44DOQDywTESmq+paj8X6AgdVtY6I9AZeAXqdubVTtmw5xGXLnVvCGjSoyMSJPa0TP2OM8SFflihaAJtUdbOqngSmAldlWuYq4P/cv6cBnSSX1sCDB44TFZHMfy6fx8qV/SxJGGOMj/myjaIqsN3jdTzQMrtlVDVFRBKACsA+z4VE5B4gvcP4pMSU4b8/9Q08FRmcTSwFqCKZjlURZsfiFDsWp9ixOKVeflf05S9tViWDzDf4erMMqvo28DaAiMTl9zH0UGPH4hQ7FqfYsTjFjsUpIhKX+1JZ82XVUzxwjsfrasDO7JYRkQigLHDAhzEZY4zJI18mimVAXRGpKSLFgd7A9EzLTAfSOyu5HvhOg62XQmOMCXE+q3py2xweAOYA4cAkVV0jIsNwBvmeDrwHfCgim3BKEr292PTbvoo5CNmxOMWOxSl2LE6xY3FKvo9F0HUzbowxpnAFfzfjxhhjfMoShTHGmBwFbKLwVfcfwciLYzFQRNaKyGoRmS8iIfsUYm7HwmO560VERSRkb4305liIyI3uubFGRD4q7BgLixffkeoi8r2IrHC/J92z2k6wE5FJIrJHRH7PZr6IyBj3OK0WkQu92rCqBtw/nMbvP4FaQHFgFdAw0zL9gYnu372BT/wdtx+PRQeghPv3fUX5WLjLlQYWAUuAWH/H7cfzoi6wAijvvq7k77j9eCzeBu5z/24IbPV33D46Fm2BC4Hfs5nfHfgG5xm2VsAv3mw3UEsUPun+I0jleixU9XtVPe6+XILzzEoo8ua8AHgBGAEkFmZwhcybY3E3ME5VDwKo6p5CjrGweHMsFEgf4rIsZz7TFRJUdRE5P4t2FfCBOpYA5USkSm7bDdREkVX3H1WzW0ZVU4D07j9CjTfHwlNfnCuGUJTrsRCRC4BzVHVGYQbmB96cF+cB54nITyKyRES6FVp0hcubY/EccKuIxAOzgAGFE1rAyevvCRC441EUWPcfIcDr9ykitwKxQDufRuQ/OR4LEQkDRgN3FFZAfuTNeRGBU/3UHqeU+YOINFbVQz6OrbB5cyxuAiar6msi0hrn+a3Gqprm+/ACSr5+NwO1RGHdf5zizbFARC4DhgBXqmpSIcVW2HI7FqWBxsACEdmKUwc7PUQbtL39jnytqsmqugX4AydxhBpvjkVf4FMAVV0MROF0GFjUePV7klmgJgrr/uOUXI+FW93yFk6SCNV6aMjlWKhqgqpWVNUaqloDp73mSlXNd2doAcyb78hXODc6ICIVcaqiNhdqlIXDm2OxDegEICINcBJFURy7djpwm3v3UysgQVV35bZSQFY9qe+6/wg6Xh6LV4FSwGdue/42Vb3Sb0H7iJfHokjw8ljMAbqIyFogFRikqvv9F7VveHksHgXeEZFHcKpa7gjFC0sR+RinqrGi2x7zLFAMQFUn4rTPdAc2AceBO73abggeK2OMMQUoUKuejDHGBAhLFMYYY3JkicIYY0yOLFEYY4zJkSUKY4wxObJEYQKOiKSKyEqPfzVyWLZGdj1l5nGfC9zeR1e5XV7Uy8c2+onIbe7fd4jI2R7z3hWRhgUc5zIRaebFOg+LSIl/um9TdFmiMIHohKo28/i3tZD2e4uqNsXpbPLVvK6sqhNV9QP35R3A2R7z7lLVtQUS5ak4x+NdnA8DlihMvlmiMEHBLTn8ICK/uv8uzmKZRiKy1C2FrBaRuu70Wz2mvyUi4bnsbhFQx123kzuGwW9uX/+R7vSX5dQYICPdac+JyGMicj1On1v/dfcZ7ZYEYkXkPhEZ4RHzHSLyZj7jXIxHh24iMkFE4sQZe+J5d9qDOAnrexH53p3WRUQWu8fxMxEplct+TBFnicIEomiPaqcv3Wl7gM6qeiHQCxiTxXr9gDdUtRnOD3W8211DL6CNOz0VuCWX/V8B/CYiUcBkoJeqno/Tk8F9IhIDXAM0UtUmwIueK6vqNCAO58q/maqe8Jg9DbjW43Uv4JN8xtkNp5uOdENUNRZoArQTkSaqOganL58OqtrB7cpjKHCZeyzjgIG57HS3W9EAAAItSURBVMcUcQHZhYcp8k64P5aeigFj3Tr5VJx+izJbDAwRkWrAF6q6UUQ6Ac2BZW73JtE4SScr/xWRE8BWnG6o6wFbVHWDO///gPuBsThjXbwrIjMBr7s0V9W9IrLZ7Wdno7uPn9zt5iXOkjjdVXiOUHajiNyD872ugjNAz+pM67Zyp//k7qc4znEzJluWKEyweATYDTTFKQmfMSiRqn4kIr8APYA5InIXTrfK/6eqT3qxj1s8OxAUkSzHN3H7FmqB08lcb+ABoGMe3ssnwI3AeuBLVVVxfrW9jhNnFLeXgXHAtSJSE3gMuEhVD4rIZJyO7zITYK6q3pSHeE0RZ1VPJliUBXa54wf0wbmaPo2I1AI2u9Ut03GqYOYD14tIJXeZGPF+TPH1QA0RqeO+7gMsdOv0y6rqLJyG4qzuPDqC0+15Vr4ArsYZI+ETd1qe4lTVZJwqpFZutVUZ4BiQICKVgcuziWUJ0Cb9PYlICRHJqnRmTAZLFCZYjAduF5ElONVOx7JYphfwu4isBOrjDPm4FucH9VsRWQ3MxamWyZWqJuL0rvmZiPwGpAETcX50Z7jbW4hT2slsMjAxvTE703YPAmuBc1V1qTstz3G6bR+vAY+p6iqc8bHXAJNwqrPSvQ18IyLfq+penDuyPnb3swTnWBmTLes91hhjTI6sRGGMMSZHliiMMcbkyBKFMcaYHFmiMMYYkyNLFMYYY3JkicIYY0yOLFEYY8z/bxTgBQDVZMeT07/OWQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fdb64306550>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%mkdir figures # The function needs a file called figures\n",
"\n",
"plot_roc(true_values, nnscore_values, 'nnscore', 'comt')"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3XmcjXX7wPHPNTPMDMa+PEX2fS+TSKHFEpKiaNFT6ZGkkpJEm/JUEiVri/zK06NSyoOIUlqIUShLEWKQ3VhnzHL9/rjvGceYOXNmzDlnluv9enk5516vc8997ut8v9/7/n5FVTHGGGMyExLsAIwxxuRtliiMMcZ4ZYnCGGOMV5YojDHGeGWJwhhjjFeWKIwxxnhliSIXicjtIvJlsOMINhGpKiLHRSQ0gPusLiIqImGB2qc/ich6EWmfg/VydA6KSKSI/E9E4kTk4+yu7+M+VERqu69niMgL/thPbsjp8S+oCmyiEJHtInLKvWD97Z6YJfy5T1X9j6p29Oc+8iL3WF+b+l5Vd6hqCVVNDmZcweJ5QcwpVW2kqt9ksZ9zkuN5nIO9gEpAOVW9OQfrFyi+HP9gEpFnRWRmoPZXYBOF63pVLQE0By4Ghgc5nhwJ5q/kgvILPTsK6fGuBvyhqknZXTGvnyN5Pb58QVUL5D9gO3Ctx/sxwHyP9+HAWGAHsBeYCkR6zL8BWAMcBf4EOrvTSwHvAHuAXcALQKg77y7ge/f1VGBsupg+B4a4ry8EPgH2A9uAhzyWexaYDcx0939vBp+vFPCeu/5fwEggxCOOH4A3gDhgE3BNunW9fYYfgPHAIXdeLeBr4CBwAPgPUNpd/n0gBTgFHAceB6oDCoS5y3wDPO9u9xjwJVDeI5473c9wEHgq/d8u3eeOBF51l48Dvnenpe7zn+7f9AAwwmO9lsBy4Ij7uScCRT3mK/AAsBnY5k57Hdjp/g1WA1d6LB8KPOmeG8fc+RcBy9xtnXCPR293+W4459MR4EegabpzdRiwDkgAwjyPgRt7jBvHXmCcO32Hu6/j7r/WeJyD7jKNgMXu33Iv8GQGx/Q54DSQ6G6nH86PyJHucd6Hc66VcpdPPdb93BiWZfK3Guoe693APe46td15M4AX3NflgXnusTkEfMeZc/ki4FOc8/wgMNGdnu34gFbusT8CrAXa+3L9wPk+fuTu4xiwHoj2sm6GxxznmvOaezx2u6/D3XntgVic788+97j1ALoAf7jbSt1O53R/r7V+v54G4qIdjH/p/tBVgF+B1z3mvwbMBcoCUcD/gBc9vphxQAf3hKwM1HfnfQZMA4oDFYGVwH3uvLs4kyja4lxkxH1fBudieqG7zdXA00BRoCawFejkcWImuidKCB4JzCP+93AST5T7xfgD6OcRRxLwCFAE6O1+nrI+foYk4EGcC1YkUNs9FuFABZyL4WsZHet0X1TPRPEnUNfd3jfAS+68hu7JfoV7LMa6nz2zRDHJXb8yzsX6cjeu1H2+5e6jGc5Ft4G7XgucC0WYu+xGYLDHdhXny1029XgDdwDl3HUeBf4GItx5Q3HOqXqAuPsr57Gt2h7bvgTny3+ZG/M/3WMW7nH81uBcFCPTH1OcBNfXfV0CaJXRcc7gHIzCueA8CkS47y/L5Lg+C8z0eH8PsAXn3CyBc7F+P91+38M5hzI6PzvjXCQbu8t8QOaJ4kWcH1ZF3H9Xusc0FOeCPt7dRgRwRU7iwzlfDuJceENwzueDQAUfrh/PAvHuuqFuvCsyWS/TYw6MAlbgfOcq4CSt59157XG+d0+7x+BfOMnxA3cbjdwYamb09/L79TQYF/GAfDDnD30c5xeAAl9x5lew4Pziq+WxfGvO/JKcBozPYJuVcC4+niWPW4GlGXxJBefXTFv3/b+Ar93XlwE70m17OPCux0mQ4a80d36oG0dDj2n3Ad94xLEbN0m501YCfX38DDsy27e7TA/gl3THOqtEMdJj/kBgofv6aeC/HvOK4fxaOidR4HzBTwHNMpiXus8q6T5zn0w+w2Bgjsd7Ba7O4nMfTt038DtwQybLpU8UU3AvCB7TfgfaeRy/ezI4f1MvVMtwfvWXT7fMWcc5g3PwVs+/Uxaf7VnOThRfAQM93tfDSeBhHvut6WV703F/DLjv65J5ohiF86OndrpttMa5WIZlsP1sxYdTYns/3TYWAf/MJH7P4/8ssMRjXkPgVCbrZXrMcX4sdfF43wnY7r5uj3Nup5bso9zPcJnH8quBHhn9vfz9r6C3UfRQ1SicP0J9nCIuONm8GLBaRI6IyBFgoTsdnF92f2awvWo42X6Px3rTcH4hnEWdv+YsnBMH4DacKpvU7VyYug13O0/iXMRT7fTyucrj/Pr+y2PaXzi/mlLtcmPwnH+hj5/hrH2LSEURmSUiu0TkKE6VWHmy52+P1ydxfgXixpS2P1U9ifNLLyPlcX6lZfS38bofEakrIvPcGxuOAv/m3M+Q/nM/KiIb3TuBjuBU2aWuk9k5kpFqwKPp/t4X4Xz2DPedTj+cC+0mEVklIt183G92YkzvQs49v8Lw/Ry9MN38vzJbEHgFp3TwpYhsFZEn3OkXAX9pxu0m2Y2vGnBzur/BFcAFXuLylP68isik7cPbMc8oZs9z4KCeuQHklPv/Xo/5pzjzvQmogp4oAFDVb3F+wYx1Jx3AOeiNVLW0+6+UOg3f4JxgtTLY1E6cX+PlPdYrqaqNMtn1f4FeIlINpxTxicd2tnlso7SqRqlqF8+wvXykAzi/nqp5TKuK096QqrKISLr5u338DOn3/aI7ramqlsSpkhEvy2fHHpyqQcC5TROnuicjB3CK3xn9bbIyBaetpo77GZ7k7M8AHp9DRK7E+RV6C1BGVUvjVN+lrpPZOZKRncDodH/vYqr634z2nZ6qblbVW3GS+cvAbBEp7m2dHMSY3m7OPb+SOPvC5W3/e3Aump7rZ0hVj6nqo6paE7geGCIi1+DEXzWTC3J249uJU6Lw/BsUV9WXvHyGnPB2zDOKeXcO93M+37lsKxSJwvUa0EFEmqtqCk5d9ngRqQggIpVFpJO77DvA3SJyjYiEuPPqq+oenIbYV0WkpDuvloi0y2iHqvoLTtH5bWCRqh5xZ60EjorIMPf+9VARaSwil/ryQdxfHR8Bo0Ukyk1EQ3B+6aeqCDwkIkVE5GagAbAgu5/BFYVTjXdERCrj1M972otTV5wTs4HrReRyESmKU8WS/gIOgPt3mw6ME5EL3ePWWkTCfdhPFE5j8HERqQ/c78PySbhVHyLyNFDSY/7bwPMiUkccTUUkNcGlPx5vAQNE5DJ32eIi0lVEonyIGxG5Q0QquJ8/9RxKdmNLIfNjPw/4h4gMFpFw91y5zJd94vzIeUREari3lf8b+DCTX/cZ+Qi4S0Qaikgx4JnMFhSRbiJS2/1hc9T9bMk435M9wEvuMYsQkTY5jG8mznnWyT1vIkSkvYhUyWT5nPJ2zP8LjBSRCiJSHqfaNae3uO4FqotIQK7hhSZRqOp+nMatp9xJw3CKuyvcqoglOPWcqOpK4G6cRrQ44FvO/BK4E6faZwNOnfVsvBdf/wtci9MolRpLMs4vp+Y4dzwdwLnwlMrGR3oQp51lK86dPx/gXERT/QTUcbc9GuilqqlVOtn9DM/hNMjGAfNxGg49vYjzBTgiIo9l4zOgquvdzzIL56JwDKfhNyGTVR7DaURehXMnyMv4dh4/hlP9dwznwv1hFssvAr7AuUngL5ySjGdVxjici+GXOBe3d3AaTcGpP/4/93jcoqoxOG1UE3GO9xactgRfdQbWi8hxnDux+qhqvFtNNxr4wd1XK8+VVPUYTqPt9ThVJ5uBq3zc53ScO9qW4Zyj8Th/J5+o6hc4P86+xvm8X3tZvA7O9+84TsP9ZFX9xuN7UhunvS8W58aMbMenqjtx7mR8EifB7sT5wZOr18AsjvkLOHevrcM5h392p+VE6kORB0Xk5xwH7KPUO3JMASIid+HcUntFsGPJLvfX4RGcKqJtwY7HGFOIShQm7xKR60WkmFvvPhbn19b24EZljEllicLkBTdw5iGkOjhVK1bUNSaPsKonY4wxXlmJwhhjjFf5rrOs8uXLa/Xq1YMdhjHG5CurV68+oKoVsl7yXPkuUVSvXp2YmJhgh2GMMfmKiHh7Ot4rq3oyxhjjlSUKY4wxXlmiMMYY45UlCmOMMV5ZojDGGOOVJQpjjDFe+S1RiMh0EdknIr9lMl9EZIKIbBGRdSJyib9iMcYYk3P+LFHMwOkeOTPX4fTrUwfojzOwjDHGmFx2+nRy1gt54bcH7lR1mYhU97LIDcB7budvK0SktIhc4A6sY4wxJr1Pu8K2BdlaZej/OvDLbl9HfM1YMNsoKnP2QDCxnD3mcxoR6S8iMSISs3///oAEZ4wxeU42kwRA43/s47utmY5E65NgduGR0XCXGXZlq6pvAm8CREdHW3e3xpjC59OuZ14/mvllcMOG/fz88x7uuKMpAHeq0u6lOGrUyOlgesFNFLGcPfh6FXI+0LgxxhRsqaWJGl0ynH3yZCIvvLCMV175kdBQoVWrKtSuXRYRoXr10ue162AmirnAIBGZBVwGxFn7hDHGZOGm+edM+uKLzTzwwAK2bTsCQL9+LShXLvKc5XLKb4lCRP4LtAfKi0gs8AxQBEBVpwILgC44A6+fBO72VyzGGJNveWnA3rXrKIMHL2L27A0ANG1aialTu9K69UUZLp9T/rzr6dYs5ivwgL/2b4wxBYJnkkhX7fTAAwv4/PPfKVasCKNGtefhh1sRFpb79yjlu/EojDGmwPHltle3ATspKSUtGbz88rUUKRLKq692pGrVUn4Lz7rwMMaYYMsqSdToQlxcPA8+uICuXT/AqZCBevXK8/HHN/s1SYCVKIwxJvfl4ME4IMPbXlWVjz/ewOAGk9iz5zihocKaNX9z8cXn9xBddliiMMaY3JaTJJHBba9//nmIQYO+YOHCLQC0bl2FqVO70bRppfONMFssURhjjL94eTAuK2PH/shTTy0lPj6J0qUjePnla7n33ksICcnoWWX/skRhjDG5IafVTZk4eTKR+Pgk+vZtytixHalYsXiubTu7LFEYY0xuSJ8kMnmCOjP795/g998PcsUVTr9Mw4a1oX376rRtWy23IswxSxTGmMIrl0sBQLarm1JSlOnTf+HxxxcTFhbCpk2DKFs2kvDwsDyRJMAShTGmMMvtJJHNUsRvv+1jwIB5/PCD05F2hw41OXkykbJlc6/7jdxgicIYY86j0TknTpw4zahR3zJu3AqSklKoVKk4r73Wmd69GyES+MbqrFiiMMYUTp7ddgdYr14fs3DhFkRg4MBoRo++htKlI4IWT1YsURhjCqcsuu32p2HD2rB373GmTOnKZZdVCfj+s8sShTGm4PPWaJ1Bt925KSkphTfe+Int24/w+uvXAdC+fXViYvoH5ZmInLBEYYwp+DJLEn4uTaxcuYv77pvHmjV/A9C/fwsaNaoIkG+SBFiiMMYUJgFqtD5yJJ4nn/yKqVNjUIVq1UoxcWKXtCSR31iiMMYULP54NiIbZs36jcGDF7J37wnCwkJ49NHWPPVUW4oXLxq0mM6XJQpjTMESpGqmVF9++Sd7956gTZuLmDKlK02aBLYDP3+wRGGMKRjSlyQCVM2UkJDErl3HqFmzDABjxnTgyiur8s9/Ns9X7RDe2MBFxpiCwcuQof7y9dfbaNp0Kl27fsDp08kAlC9fjLvvvrjAJAmwEoUxpqAJQEli797jPPbYYmbOXAdA/frliY09mlaqKGgsURhj8q4gN0ynl5KivPXWap544iuOHIknIiKMkSOvZOjQNhQtGhrs8PzGEoUxJu/KbpLwc5XTjTd+yNy5vwPQqVMtJk3qQq1aZf26z7zAEoUxJm/JqBQR4E77MnPTTfVZuXIXr7/emZtvbpgnO/DzB0sUxpi85TwHAMpNc+f+TmzsUQYOvBSAO+9sxk03NSAqKjxoMQWDJQpjTPB4a4MIYilix444HnroCz7//HfCw0Pp3Lk2NWuWQUQKXZIASxTGmGAK8sNx6SUmJjNhwk8888w3nDiRSFRUUV544WqqVSsVlHjyCksUxpjgywNtECtWxHLfffNYt24vADff3JDx4ztRuXLJIEcWfJYojDH+kcdubc3KU08tZd26vdSoUZqJE7vQpUudYIeUZ1iiMMb4h69JIkjVTKrKsWOnKVnSaXOYOPE63ntvLSNGtKVYsSJBiSmvskRhjDlbbpcE8kC1Unq//36AgQMXIAKLF/dFRKhXrzyjR18T7NDyJEsUxpiz5WaSCOKtrRmJj0/ixRe/46WXfuD06WTKlYtk+/Yj1KhRMLveyC2WKIwxGcuDJYHzsXjxnwwcuIAtWw4BcM89zRkzpgPlyhULcmR5n18ThYh0Bl4HQoG3VfWldPOrAv8HlHaXeUJV80/rlzF5WT5rTPYXVaVfv7m8++4aABo2rMDUqV258spqQY4s//BbohCRUGAS0AGIBVaJyFxV3eCx2EjgI1WdIiINgQVAdX/FZEyhcj5JIo9VGZ0PEaF69dJERobx9NPtGDKkdYHuwM8f/FmiaAlsUdWtACIyC7gB8EwUCqTepFwK2O3HeIwpnApYFZIv1qz5mz17jnHddc4trsOGtaFv36bWFpFD/hy4qDKw0+N9rDvN07PAHSISi1OaeDCjDYlIfxGJEZGY/fv3+yNWY0wBcOxYAkOGLKJFizf55z8/49ChUwCEh4dZkjgP/kwUGXWrmP6nza3ADFWtAnQB3heRc2JS1TdVNVpVoytUqOCHUI0x+ZmqMmfORho2nMz48SsAuO22JhQpYoN45gZ/Vj3FAhd5vK/CuVVL/YDOAKq6XEQigPLAPj/GZYwpQP766wiDBn3BvHl/ABAdfSHTpnXjkksuCHJkBYc/0+0qoI6I1BCRokAfYG66ZXYA1wCISAMgArC6JWOMT1SVnj0/Yt68PyhZMpyJE69jxYp+liRymd9KFKqaJCKDgEU4t75OV9X1IjIKiFHVucCjwFsi8ghOtdRdqlr4Wt6MOR+F8DbYlBQlJEQQEcaO7cjUqTGMH9+JCy6ICnZoBZLkt+tydHS0xsTEBDsMY/KOV72MslajC9w0P3Cx+NnBgyd54oklALz1VvcgR5O/iMhqVY3Oybr2ZLYx+VX6kkQBvg1WVXnvvbU89thiDhw4SdGioTzzTHuqVLEuwAPBEoUx+ZVnkihAD8ilt3Hjfu6/fz7ffvsXAO3bV2fKlK6WJALIEoUx+V0BLUmoKk8/vZSXX/6BxMQUypcvxquvdqRv36aIeKluM7nOEoUxgVIIG53Ph4iwa9cxEhNT+Ne/LuGll66lbNnIYIdVKFmiMCZQ/JEkCliV0+7dxzhw4CRNm1YCYMyYDvTrdzFt2lQNcmSFmyUKY3JbViWHAlpVdD6Sk1OYMiWGESO+pnLlKNasGUDRoqGUL1+M8uUtSQSbJQpjcpu3JFHASgC54eef93DfffOIiXE6bmjbthpHjyZQvryNE5FX+JQo3Cerq6rqFj/HY0zBYSUHr44eTeCpp75m4sRVpKQoVaqUZMKEzvToUd8aq/OYLBOFiHQFxgFFgRoi0hx4RlVv9HdwxgSFNTr7narStu27rF27l9BQYciQVjz7bHuiosKDHZrJgC99PY0CLgOOAKjqGqC2P4MyJqhyI0lYFZNXIsIjj7SiZcvKxMT059VXO1mSyMN8qXpKVNUj6YqCVqY2+ZsvpQarOso1p08nM27cckJDhaFD2wBw553NuOOOpoSGWlfgeZ0viWKjiNwChIhIDeBhYIV/wzLGz7JKElYiyDXfffcXAwbMZ8OG/YSHh3Lnnc2oVKkEIkJoqLVF5Ae+JIpBwNNACvApTm+ww/0ZlDEBY6UGvzlw4CSPP76Yd99dA0CdOmWZPLkrlSqVCHJkJrt8SRSdVHUYMCx1gojchJM0jMl/Pu0a7AgKNFVlxow1DB26mIMHT1G0aCjDh1/BE09cQUSE3ZGfH/lSOTgyg2kjcjsQYwImtdrJqpf8ZubMXzl48BRXX12DdesG8Oyz7S1J5GOZ/uVEpBPOMKWVRWScx6ySONVQxuQ92bm1tQCN0xBsJ08mEhcXzwUXRCEiTJ7chVWrdnP77U3smYgCwFuK3wf8BsQD6z2mHwOe8GdQxuSYr0nCShO55osvNvPAAwuoWbMMixf3RUSoV6889eqVD3ZoJpdkmihU9RfgFxH5j6rGBzAmY7KvEA3ik1fs2nWUwYMXMXv2BgCiosI5ePCUdb1RAPlSaVhZREYDDYGI1ImqWtdvURmTXYVkEJ+8IDk5hUmTVjFy5NccO3aa4sWLMGrUVTz00GWEhdkzEQWRL4liBvACMBa4Drgba6MweZWVJPwqJUVp124GP/ywE4AePerz+uudqVq1VJAjM/7kS6IopqqLRGSsqv4JjBSR7/wdmDFprO+lPCMkROjYsRY7dsQxcWIXunevF+yQTAD4kigSxLlt4U8RGQDsAir6NyxjPFgDddCoKh99tJ6wsBB69mwIwLBhbRgypDUlShQNcnQmUHxJFI8AJYCHgNFAKeAefwZlTBrPh+OsWimg/vzzEAMHLuDLL/+kQoViXH11DcqUiSQ8PIxw67+vUMkyUajqT+7LY0BfABGp4s+gjEljD8cFXEJCEq+88iOjR39HfHwSZcpEMHr01ZQqFZH1yqZA8pooRORSoDLwvaoeEJFGOF15XA1YsjCBYw/HBcQ332zn/vvns2nTAQD69m3K2LEdqVixeJAjM8GU6b1sIvIi8B/gdmChiIwAlgJrAbs11vjXp13hVXuiN5CSk1MYONBJEvXqlePrr+/kvfdutCRhvJYobgCaqeopESkL7Hbf/x6Y0EyhZs9FBERKihIfn0SxYkUIDQ1hypSuLFv2F48/3obwcOubyTi8nQnxqnoKQFUPicgmSxLGL7zd/moN2H7z6697GTBgPvXrl+Odd24AoF276rRrVz24gZk8x1uiqCkiqV2JC1Dd4z2qepNfIzOFR2ZJwkoSfnHixGlGjfqWceNWkJSUwrZthzl8+BRlykQGOzSTR3lLFD3TvZ/oz0CMsdKD//3vf78zaNAX7NgRhwgMHBjN6NHXULq03dFkMuetU8CvAhmIMcZ/kpJS6N17Np9+uhGA5s3/wbRp3WjZsnKQIzP5gbVWGVMIhIWFUKpUOCVKFOX5569i0KCW1oGf8ZlfzxQR6Swiv4vIFhHJcAwLEblFRDaIyHoR+cCf8RhTmPz0Uyw//RSb9v6VVzqwceMDDB7cypKEyRafSxQiEq6qCdlYPhSYBHQAYoFVIjJXVTd4LFMHGA60UdXDImJ9SBlzno4ciWf48CVMm7aa+vXLs2bNAIoWDaVcORsnwuRMlj8rRKSliPwKbHbfNxORN3zYdktgi6puVdXTwCycZzM8/QuYpKqHAVR1X7aiN8akUVU++OBX6tefyNSpqwkNDaF793okJ9uoAOb8+FKimAB0Az4DUNW1InKVD+tVBnZ6vI8FLku3TF0AEfkBCAWeVdWFPmzb5EfWXbjfbN58kIEDF7BkyVYA2rS5iKlTu9G4sRXSzfnzJVGEqOpf6QZIT/ZhvYz6X0h//2MYUAdoj9N31Hci0lhVj5y1IZH+QH+AqlWr+rBrkyd5SxL2zESOJSYmc/XV7xEbe5SyZSMZM+Za7r77YkJCrAsUkzt8SRQ7RaQloG67w4PAHz6sFwtc5PG+Ck43IOmXWaGqicA2EfkdJ3Gs8lxIVd8E3gSIjo62m+3zGxvP2i9UFRGhSJFQRo++mqVLtzNmzLVUqGB9M5nc5cutD/cDQ4CqwF6glTstK6uAOiJSQ0SKAn2AuemW+Qy4CkBEyuNURW31LXSTb1i/Tblq797j9O07hxdeWJY27c47m/HuuzdYkjB+4UuJIklV+2R3w6qaJCKDgEU47Q/TVXW9iIwCYlR1rjuvo4hswKnOGqqqB7O7L5NHWUkiV6WkKG+9tZonnviKI0fiKV06gsGDWxEVZaMIGf/yJVGscquEPgQ+VdVjvm5cVRcAC9JNe9rjteKUVob4uk2Tj1hJItesXfs3AwbMZ8UK57mIzp1rM2lSF0sSJiB8GeGulohcjlN19JyIrAFmqeosv0dnCgYrSeRYYmIyw4d/xWuvrSA5WbngghK8/npnevVqSLobTIzxG58ez1TVH1X1IeAS4CjOgEbGZM5zrGuTY2FhIfzyy9+kpCgPPtiSjRsf4OabG1mSMAGVZYlCRErgPCjXB2gAfA5c7ue4TH5nY13n2I4dcSQnp1CjRhlEhKlTuxIXl0B09IXBDs0UUr60UfwG/A8Yo6rf+Tkek1/4+vCcjXXts8TEZF5//SeeeeYbWreuwuLFfRER6tQpF+zQTCHnS6KoqarWB4A5my9JwkoTPlu+fCcDBsxn3bq9AJQtG8nJk4kUL140yJEZ4yVRiMirqvoo8ImInNMaaSPcGcAaqs/T4cOneOKJJbz55s8A1KhRmkmTunDddXWCHJkxZ3grUXzo/m8j2xVG1i+T3yUkJNG8+TR27IijSJEQhg69nBEj2lKsWJFgh2bMWbyNcLfSfdlAVc9KFu6DdDYCXkFmVUt+Fx4eRr9+F/PVV9uYMqUrDRtWCHZIxmRInGfevCwg8rOqXpJu2i+qerFfI8tEdHS0xsTEBGPXBYuvJQarWso18fFJvPjid9SrV57bbmsCOEOUhoaK3e5q/E5EVqtqdE7W9dZG0RvnltgaIvKpx6wo4EjGa5l8w0oMAbV48Z8MHLiALVsOUbFicW68sT6RkUVspDmTL3hro1gJHMTp9XWSx/RjwC/+DMoEkJUY/Orvv48zZMgi/vvf3wBo1KgCU6d2IzLS2iFM/uGtjWIbsA1YErhwjF9ZA3XAJCenMG3aap588ivi4hKIjAzjmWfa8cgjrSlaNDTY4RmTLd6qnr5V1XYicpizBxwSnP78yvo9OpO70icJq1rym+Rk5Y03VhIXl0CXLnWYOPE6atQoE+ywjMkRb1VPqcOdlg9EIMbPPPtesuomvzh2LIHkZKV06QiKFg3lrbeuZ+/e49x0UwNrrDb5WqYtaR5PY18EhKr8s/lwAAAgAElEQVRqMtAauA+w0VHyG+t7yW9UlU8/3UiDBpN49NFFadOvuKIqPXtaL68m//PllovPcIZBrQW8h9Mx4Ad+jcr4j/W9lKu2bz9C9+6z6NnzI3btOsZvv+0nPj4p2GEZk6t8SRQp7pjWNwGvqeqDQGX/hmVM3paYmMzLL39Pw4aTmDfvD0qWDGfixOv48cd7iIjwpQs1Y/IPn4ZCFZGbgb5AD3ea3dtnCq2TJxNp1eptfv11HwB9+jRm3LiOXHBBVJAjM8Y/fEkU9wADcboZ3yoiNYD/+jcskyvsdli/KFasCNHRF3LyZCKTJ3elY8dawQ7JGL/yZSjU30TkIaC2iNQHtqjqaP+HZs6b3Q6bK1SV995bS61aZbniiqoAjB/fiaJFQ+3BOVMo+DLC3ZXA+8AunGco/iEifVX1B38HZ86D3Q6bKzZu3M/998/n22//okGD8qxZM4CiRUMpVSoi2KEZEzC+VD2NB7qo6gYAEWmAkzhy1LmUCRC7Hfa8nDqVyOjR3zFmzA8kJqZQoUIxhg+/giJFrG8mU/j4kiiKpiYJAFXdKCI27FZ+YbfDZtvChVt44IEFbN16GIB//esSXnrpWsqWjQxyZMYEhy+J4mcRmYZTigC4HesUMG/zrHYy2XL8+Gn69p3DgQMnady4IlOndqVNm6rBDsuYoPIlUQwAHgIex2mjWAa84c+gzHmyaqdsSU5OISVFKVIklBIlivL6652JjT3KI4+0okgR68DPGK+JQkSaALWAOao6JjAhmWzL7DZYq3bK0urVu7nvvnnccEM9nnqqHUDaoELGGEemLXMi8iRO9x23A4tF5J6ARWWyJ6MkYaUJr44eTeDhh7+gZcu3Wb16D++/v47ExORgh2VMnuStRHE70FRVT4hIBWABMD0wYZkcsdtgs6SqzJ69gYcfXsiePccJDRWGDGnFc89dZdVMxmTCW6JIUNUTAKq6X0TsvsC8xJ66zrZjxxLo3Xs2X3yxBYDLLqvM1KndaN78H0GOzJi8zVuiqOkxVrYAtTzHzlbVm/wamfHOnrrOthIlipKQkEypUuG89NK19O/fgpAQ6wLcmKx4SxQ9072f6M9AjI/SlySsusmrZcv+4oILSlCnTjlEhOnTuxMREUalSiWCHZox+Ya3MbO/CmQgxkeeScJKEZk6cOAkjz++mHffXcM119Rg8eK+iAjVqpUOdmjG5DvWcX5+ZSWJDKWkKDNmrGHo0MUcOnSKokVDufLKqiQnK2FhVs1kTE74tYFaRDqLyO8iskVEnvCyXC8RURGx/qNMjq1fv4/27WfQr99cDh06xTXX1ODXX+/nmWfaExZm92IYk1M+lyhEJFxVE7KxfCgwCegAxAKrRGSuZ79R7nJROE9+/+Trto1JLy4unlat3uH48dNUrFicceM6ctttTWy8amNyQZY/s0SkpYj8Cmx23zcTEV+68GiJM3bFVlU9DcwCbshgueeBMUC872Eb41B1quBKlYpg2LA2DBjQgk2bHuD225takjAml/hSHp8AdAMOAqjqWuAqH9arDOz0eB9LurG2ReRi4CJVnedtQyLSX0RiRCRm//79PuzaFHS7dh2lV6+PmDlzXdq0ESOuZMqUbpQpY728GpObfKl6ClHVv9L9OvOlr4OMfs6ltcC6D/CNB+7KakOq+ibwJkB0dHTBb8W1h+kylZSUwqRJKxk5cinHj5/m55/3cNttTQgNDbEShDF+4kui2CkiLQF12x0eBP7wYb1Y4CKP91WA3R7vo4DGwDfuF/wfwFwR6a6qMb4EX2BllSQK6W2xq1btYsCA+fz88x4AevSoz4QJnQkNtYZqY/zJl0RxP071U1VgL7DEnZaVVUAdEamBM4xqH+C21JmqGgeUT30vIt8AjxX6JOHJboEF4MSJ0wwbtoTJk1ehClWrluKNN66je/d6wQ7NmEIhy0ShqvtwLvLZoqpJIjIIWASEAtNVdb2IjAJiVHVutqMtyKy6KVNhYSEsWbKVkBBhyJDWPPNMO4oXt0EWjQmULBOFiLyFR9tCKlXtn9W6qroAp9dZz2lPZ7Js+6y2V6BZ301n+fPPQ5QuHUG5csUIDw/j/fdvJCIijCZNKgU7NGMKHV+qnpZ4vI4AbuTsu5lMbirk1U0JCUm88sqPjB79Hbff3oS33+4OwKWXVs5iTWOMv/hS9fSh53sReR9Y7LeITKH1zTfbuf/++WzadABw7nBKTk6xxmpjgiwnfT3VAKrldiCm8Nq37wRDhy7mvffWAlCvXjmmTOnKVVfVCHJkxhjwrY3iMGfaKEKAQ0Cm/TaZHPi0a7AjCJoDB07SoMEkDh06RXh4KCNGXMnjj7chPNz6qzQmr/D6bRTnAYdmOLe3AqRoap8JJvekNmQXwgbs8uWLccMN9YiNPcrkyV2pXbtssEMyxqTjNVGoqorIHFVtEaiACpX0t8TeND94sQTIiROnGTXqW7p2rUvbtk4N5uTJXQkPD7Unq43Jo3xpJVwpIpf4PZLCqJANQvS///1Ow4aTGTPmRwYOnE9KilM4jYgIsyRhTB6WaYlCRMJUNQm4AviXiPwJnMDpw0lV1ZJHbingt8Tu3BnHww8vZM6cTQBcfPE/mDatm41XbUw+4a3qaSVwCdAjQLEULoWgATspKYUJE37i6aeXcuJEIiVKFOWFF67igQda2kBCxuQj3hKFAKjqnwGKpXApBA3YR48m8OKL33PiRCI9ezbgtdc6U6VKyWCHZYzJJm+JooKIDMlspqqO80M8hYNnaaKANWAfORJPZGQY4eFhlC0bybRp3QgPD6Vr17rBDs0Yk0Peyv+hQAmc7sAz+mdyqgCWJlSVDz74lXr1JjJmzA9p02+6qYElCWPyOW8lij2qOipgkRQWBbA08ccfBxk4cD5ffbUNgGXLdqCqdieTMQVElm0UJpcVoNJEfHwSL7/8Pf/+9/ecPp1M2bKRvPJKB+66q7klCWMKEG+J4pqARVEY5fPSxN9/H6dt23fZvPkQAHfd1ZxXXulA+fLFghyZMSa3ZZooVPVQIAMx+UulSsW56KJShIWFMGVKV9q1qx7skIwxfmI9rxmfpKQob721mquuqkHduuUQET744CbKlImkaNHQYIdnjPEje+opkPLpQ3Zr1/5NmzbTGTBgPgMHzie1X8hKlUpYkjCmELASRSDls4bs48dP8+yz3/DaaytITlYuvDCKAQOigx2WMSbALFEEQz5oyP7ss008+OAXxMYeJSREePDBlrzwwtWULBke7NCMMQFmiSJQ8lG1065dR+nTZzYJCcm0aHEBU6d2Izr6wmCHZYwJEksUgZLHq50SE5MJCwtBRKhcuSSjR19N0aKhDBx4qY1ZbUwhZ1eAQMjjT2P/+ONOWrR4k5kz16VNe/TRy3nwwcssSRhjLFEERB4tTRw6dIr77vsfbdpM59df9zF5cgw20q0xJj2regqkPFKaUFVmzlzHo49+yf79JylSJITHH2/DiBFXWtcbxphzWKLwtzzWiL1373FuvfUTli7dDkC7dtWYMqUrDRpUCG5gxpg8yxKFv+WxaqfSpSPYs+c45csXY+zYDtx5ZzMrRRhjvLJE4Q+fdj2TIFIFsdpp8eI/ueSSCyhXrhjh4WF8/PHNXHBBCcqVsw78jDFZs8Zsf0ifJIJUmtiz5xi33voJHTvOZNiwJWnTGzeuaEnCGOMzK1H406PBuYMoOTmFadNWM3z4Vxw9mkBkZBj16pWzwYSMMTliiSK3Bbnx+uef9zBgwDxWrdoNQNeudZg4sQvVq5cOalzGmPzLEkVuC2Lj9fbtR2jZ8i2Sk5XKlaOYMOE6bryxvpUijDHnxa+JQkQ6A68DocDbqvpSuvlDgHuBJGA/cI+q/uXPmAImCI3X1auX5u67mxMVFc5zz7UnKso68DPGnD+/NWaLSCgwCbgOaAjcKiIN0y32CxCtqk2B2cAYf8VTEG3ffoTrr/8v3367PW3am29ez7hxnSxJGGNyjT9LFC2BLaq6FUBEZgE3ABtSF1DVpR7LrwDu8GM8/pXRLbF+kpiYzLhxy3nuuW85dSqJAwdOsnx5PwCrZjLG5Dp/JorKwE6P97HAZV6W7wd8kdEMEekP9AeoWrVqbsWXuzyThB/bJ77/fgcDBsxj/fr9APTp05hx4zr6bX/GGOPPRJHRT9sM7xcVkTuAaKBdRvNV9U3gTYDo6Oi83Wudn26JPXz4FEOHLuadd34BoFatMkye3JWOHWv5ZX/GGJPKn4kiFrjI430VYHf6hUTkWmAE0E5VE/wYT76WkqJ8/vnvFCkSwhNPXMHw4VcQGVkk2GEZYwoBfyaKVUAdEakB7AL6ALd5LiAiFwPTgM6qus+PseRLmzYdoEaN0oSHh1GuXDH+85+bqFq1FPXrlw92aMaYQsRvdz2pahIwCFgEbAQ+UtX1IjJKRLq7i70ClAA+FpE1IjLXX/H4zadd4dXcbUA+eTKRESO+omnTKYwZ80Pa9I4da1mSMMYEnF+fo1DVBcCCdNOe9nh9rT/3HxC53Ii9cOEWBg6cz7ZtRwA4cODkeW/TGGPOhz2ZnVvOsxF79+5jDB68kI8/du4ebtKkIlOnduPyyy/KYk1jjPEvSxQ5kcvPTPzxx0Gio9/k2LHTFCtWhGefbcfgwa0oUiQ01/ZhjDE5ZYkiJ3K5G/E6dcpy6aWVKV68CG+8cR3VqlkHfsaYvMMSxfnIYXXT0aMJPP30UgYOvJS6dcshIsyd24fixYvmcoDGGHP+LFEEkKoye/YGHn54IXv2HGfTpgMsXOj0WmJJwhiTV1miCJCtWw8zaNACvvhiCwCtWlXh5Zfz/01fxpiCzxJFdmVzYKLTp5MZO/ZHnn9+GfHxSZQuHcFLL13Dv/7VgpAQ68DPGJP3WaLIrmwOTLRzZxyjRn1LQkIyt9/ehFdf7UilSiX8GKAxxuQuSxQ55WVgosOHT1G6dAQiQq1aZXn99c7Url2Wa66pGcAAjTEmd/itC4/CKCVFmT79F2rXfoOZM9elTb/vvmhLEsaYfMsSha+y6NNp/fp9tG8/g3795nLo0Km0RmtjjMnvrOrJV5n06XTyZCLPP/8tY8cuJykphYoVizN+fCduvbVxEII0xpjcZ4kiuzwesvvjj4N06jST7duPIAIDBrTg3/++hjJlIoMYoDHG5C5LFFnx0q9TtWqliIgIo1mzSkyd2o1WraoEODiTFyQmJhIbG0t8fHywQzGGiIgIqlSpQpEiuTewmSWKrHgkiaSqXZg6cSW33tqYcuWKER4exsKFt1O5cknCwqy5p7CKjY0lKiqK6tWrI2LPxpjgUVUOHjxIbGwsNWrUyLXt2tXNRyuvjKXli9158MEvGDZsSdr0atVKW5Io5OLj4ylXrpwlCRN0IkK5cuVyvXRrJYosxJ0KZ8QX1zB56NuoQtWqpbjhhnrBDsvkMZYkTF7hj3PREkUmVJUPP1zPI2MG8fexKMLCQhgypBVPP93OOvAzxhQqVmeSibVr93LrrZ/w97EoLq++g59/7s/LL3ewJGHypNDQUJo3b07jxo25/vrrOXLkSNq89evXc/XVV1O3bl3q1KnD888/j+qZu/e++OILoqOjadCgAfXr1+exxx4Lxkfw6pdffuHee+8Ndhhevfjii9SuXZt69eqxaNGiDJdRVUaMGEHdunVp0KABEyZMOGv+qlWrCA0NZfbs2WnThg0bRuPGjWncuDEffvhh2vQ+ffqwefNm/3yYjALPT/9atGih/pKUlHzW+0fadtK3br5Ekz/u4rd9mvxvw4YNwQ5Bixcvnvb6zjvv1BdeeEFVVU+ePKk1a9bURYsWqarqiRMntHPnzjpx4kRVVf3111+1Zs2aunHjRlVVTUxM1EmTJuVqbImJiee9jV69eumaNWsCus/sWL9+vTZt2lTj4+N169atWrNmTU1KSjpnuenTp2vfvn01Odm51uzduzdtXlJSkl511VV63XXX6ccff6yqqvPmzdNrr71WExMT9fjx49qiRQuNi4tTVdVvvvlG77333gzjyeicBGI0h9ddq3pyLV26jYEDFzBtWjfaHhgI2xYwrrs7s9f5jYdtChEvT++fl2wMktW6dWvWrXO6kPnggw9o06YNHTt2BKBYsWJMnDiR9u3b88ADDzBmzBhGjBhB/fr1AQgLC2PgwIHnbPP48eM8+OCDxMTEICI888wz9OzZkxIlSnD8+HEAZs+ezbx585gxYwZ33XUXZcuW5ZdffqF58+bMmTOHNWvWULq0M3pj7dq1+eGHHwgJCWHAgAHs2LEDgNdee402bdqcte9jx46xbt06mjVrBsDKlSsZPHgwp06dIjIyknfffZd69eoxY8YM5s+fT3x8PCdOnODrr7/mlVde4aOPPiIhIYEbb7yR5557DoAePXqwc+dO4uPjefjhh+nfv7/Pxzcjn3/+OX369CE8PJwaNWpQu3ZtVq5cSevWrc9absqUKXzwwQeEhDiVORUrVkyb98Ybb9CzZ09WrVqVNm3Dhg20a9eOsLAwwsLCaNasGQsXLuSWW27hyiuv5K677iIpKYmwMP9eygt9oti37wRDhy7mvffWAjBu3HLaXpnxU9jG5HXJycl89dVX9OvXD3CqnVq0aHHWMrVq1eL48eMcPXqU3377jUcffTTL7T7//POUKlWKX3/9FYDDhw9nuc4ff/zBkiVLCA0NJSUlhTlz5nD33Xfz008/Ub16dSpVqsRtt93GI488whVXXMGOHTvo1KkTGzduPGs7MTExNG58pqeD+vXrs2zZMsLCwliyZAlPPvkkn3zyCQDLly9n3bp1lC1bli+//JLNmzezcuVKVJXu3buzbNky2rZty/Tp0ylbtiynTp3i0ksvpWfPnpQrV+6s/T7yyCMsXbr0nM/Vp08fnnjiibOm7dq1i1atWqW9r1KlCrt27Tpn3T///JMPP/yQOXPmUKFCBSZMmECdOnXYtWsXc+bM4euvvz4rUTRr1oznnnuOIUOGcPLkSZYuXUrDhg0BCAkJoXbt2qxdu/acv3FuK7SJIiVFeeednxk2bAmHD8cTHh7KyJFtGTr0cph4q7NQDoc6NYVYkM6ZU6dO0bx5c7Zv306LFi3o0KED4FQtZ3YXTHbujlmyZAmzZs1Ke1+mTJks17n55psJDQ0FoHfv3owaNYq7776bWbNm0bt377TtbtiwIW2do0ePcuzYMaKiotKm7dmzhwoVKqS9j4uL45///CebN29GREhMTEyb16FDB8qWLQvAl19+yZdffsnFF18MOKWizZs307ZtWyZMmMCcOXMA2LlzJ5s3bz4nUYwfP963gwNntfmkyuj4JiQkEBERQUxMDJ9++in33HMP3333HYMHD+bll19OO16pOnbsyKpVq7j88supUKECrVu3Pqv0ULFiRXbv3m2Jwh+2bTvMHXfM4ccfdwLQsWMtJk3qQu11fWFixk9hG5OXRUZGsmbNGuLi4ujWrRuTJk3ioYceolGjRixbtuysZbdu3UqJEiWIioqiUaNGrF69Oq1aJzOZJRzPaenv3S9evHja69atW7Nlyxb279/PZ599xsiRIwFISUlh+fLlREZm3u1NZGTkWdt+6qmnuOqqq5gzZw7bt2+nffv2Ge5TVRk+fDj33XffWdv75ptvWLJkCcuXL6dYsWK0b98+w+cOslOiqFKlCjt37kx7Hxsby4UXXnjOulWqVKFnz54A3Hjjjdx9992AU2rq06cPAAcOHGDBggWEhYXRo0cPRowYwYgRIwC47bbbqFOnTtr24uPjvR673FIo73oqWTKcP/44yD/+UYJZs3qycOHt1K5dNtOO/4zJL0qVKsWECRMYO3YsiYmJ3H777Xz//fcsWeI8JHrq1CkeeughHn/8cQCGDh3Kv//9b/744w/AuXCPGzfunO127NiRiRMnpr1PrXqqVKkSGzduTKtayoyIcOONNzJkyBAaNGiQ9us9/XbXrFlzzroNGjRgy5YzvTHHxcVRuXJlAGbMmJHpPjt16sT06dPT2lB27drFvn37iIuLo0yZMhQrVoxNmzaxYsWKDNcfP348a9asOedf+iQB0L17d2bNmkVCQgLbtm1j8+bNtGzZ8pzlevTowddffw3At99+S926dQHYtm0b27dvZ/v27fTq1YvJkyfTo0cPkpOTOXjwIADr1q1j3bp1ae1N4FTvNWrUKNNjkFsKTaJYtGgLCQlJAJQrV4y5c/uwacy39I5tgowLObsR8lH1OjCRMXnZxRdfTLNmzZg1axaRkZF8/vnnvPDCC9SrV48mTZpw6aWXMmjQIACaNm3Ka6+9xq233kqDBg1o3Lgxe/bsOWebI0eO5PDhwzRu3JhmzZql/dJ+6aWX6NatG1dffTUXXHCB17h69+7NzJkz06qdACZMmEBMTAxNmzalYcOGTJ069Zz16tevT1xcHMeOHQPg8ccfZ/jw4bRp04bk5ORM99exY0duu+02WrduTZMmTejVqxfHjh2jc+fOJCUl0bRpU5566qmz2hZyqlGjRtxyyy00bNiQzp07M2nSpLRqpC5durB7924AnnjiCT755BOaNGnC8OHDefvtt71uNzExkSuvvJKGDRvSv39/Zs6cmVb1tHfvXiIjI7M87rlBMqpby8uio6M1JibG5+V37ozjoYcW8tlnm3j++asYObLtmZkZ3aFSo4slCZMtGzdupEGDBsEOo0AbP348UVFRef5ZikAaP348JUuWTLtxwVNG56SIrFbV6Jzsq8C2USQlpTBhwk88/fRSTpxIpESJopTd8S682u7cha3R2pg87f777+fjjz8Odhh5SunSpenbt29A9lUgE8WKFbEMGDCPtWv3AtCzZwNef70zlWeVOndha4swJs+LiIgI2EUxv0htCA+EApcofvoplssvfwdVqF69NBMnXkfXrnXPXshKECaXebsN1ZhA8kdzQoFLFC1bVqZTp9pcfPE/GDmyLcWK5d7gHcZkJCIigoMHD1pX4ybo1B2PIiIiIle3m+8TxebNB3nkkUWMG9eJunWdL+r8+bcREmJfWBMYVapUITY2lv379wc7FGPSRrjLTfk2USQkJPHSS9/z4ovfk5CQTEREGLNn3wJwJkl4GcbUmNxSpEiRXB1NzJi8xq/PUYhIZxH5XUS2iMg5T6mISLiIfOjO/0lEqvuy3a++2krTplN59tlvSUhI5u67mzN1ardzF8woSVjjtTHGZIvfShQiEgpMAjoAscAqEZmrqhs8FusHHFbV2iLSB3gZ6H3u1s7Ytu0I1177PgANGpRn6tRutG1bzXsw1nhtjDE55s+qp5bAFlXdCiAis4AbAM9EcQPwrPt6NjBRRES9NNsfPnSSiLBEnu7wLY+2W07RVYNgVWZLG2OMOV9+ezJbRHoBnVX1Xvd9X+AyVR3kscxv7jKx7vs/3WUOpNtWfyC1w/jGwG9+CTr/KQ8cyHKpwsGOxRl2LM6wY3FGPVWNynqxc/mzRJHRbUfps5Ivy6CqbwJvAohITE4fQy9o7FicYcfiDDsWZ9ixOENEfO/7KB1/NmbHAhd5vK8C7M5sGREJA0oBh/wYkzHGmGzyZ6JYBdQRkRoiUhToA8xNt8xc4J/u617A197aJ4wxxgSe36qeVDVJRAYBi4BQYLqqrheRUTiDfM8F3gHeF5EtOCWJPj5s+k1/xZwP2bE4w47FGXYszrBjcUaOj0W+62bcGGNMYBWagYuMMcbkjCUKY4wxXuXZROGv7j/yIx+OxRAR2SAi60TkKxHJ4lH1/CurY+GxXC8RUREpsLdG+nIsROQW99xYLyIfBDrGQPHhO1JVRJaKyC/u96RA9uUjItNFZJ/7jFpG80VEJrjHaZ2IXOLThlU1z/3Dafz+E6gJFAXWAg3TLTMQmOq+7gN8GOy4g3gsrgKKua/vL8zHwl0uClgGrACigx13EM+LOsAvQBn3fcVgxx3EY/EmcL/7uiGwPdhx++lYtAUuAX7LZH4X4AucZ9haAT/5st28WqJI6/5DVU8Dqd1/eLoB+D/39WzgGimYgwFkeSxUdamqnnTfrsB5ZqUg8uW8AHgeGAPEBzK4APPlWPwLmKSqhwFUdV+AYwwUX46FAiXd16U495muAkFVl+H9WbQbgPfUsQIoLSIXZLXdvJooKgM7Pd7HutMyXEZVk4A4oFxAogssX46Fp344vxgKoiyPhYhcDFykqvMCGVgQ+HJe1AXqisgPIrJCRDoHLLrA8uVYPAvcISKxwALgwcCEludk93oC5N3xKHKt+48CwOfPKSJ3ANFAO79GFDxej4WIhADjgbsCFVAQ+XJehOFUP7XHKWV+JyKNVfWIn2MLNF+Oxa3ADFV9VURa4zy/1VhVU/wfXp6So+tmXi1RWPcfZ/hyLBCRa4ERQHdVTQhQbIGW1bGIwuk08hsR2Y5TBzu3gDZo+/od+VxVE1V1G/A7TuIoaHw5Fv2AjwBUdTkQgdNhYGHj0/UkvbyaKKz7jzOyPBZudcs0nCRRUOuhIYtjoapxqlpeVauranWc9pruqprjztDyMF++I5/h3OiAiJTHqYraGtAoA8OXY7EDuAZARBrgJIrCOHbtXOBO9+6nVkCcqu7JaqU8WfWk/uv+I9/x8Vi8ApQAPnbb83eoavegBe0nPh6LQsHHY7EI6CgiG4BkYKiqHgxe1P7h47F4FHhLRB7BqWq5qyD+sBSR/+JUNZZ322OeAYoAqOpUnPaZLsAW4CRwt0/bLYDHyhhjTC7Kq1VPxhhj8ghLFMYYY7yyRGGMMcYrSxTGGGO8skRhjDHGK0sUJs8RkWQRWePxr7qXZatn1lNmNvf5jdv76Fq3y4t6OdjGABG50319l4hc6DHvbRFpmMtxrhKR5j6sM1hEip3vvk3hZYnC5EWnVLW5x7/tAdrv7araDKezyVeyu7KqTlXV99y3dwEXesy7V1U35EqUZ+KcjG9xDgYsUZgcs0Rh8gW35PCdiPzs/jnS5CAAAANVSURBVLs8g2UaichKtxSyTkTquNPv8Jg+TURCs9jdMqC2u+417hgGv7p9/Ye701+SM2OAjHWnPSsij4lIL5w+t/7j7jPSLQlEi8j9IjLGI+a7ROSNHMa5HI8O3URkiojEiDP2xHPutIdwEtZSEVnqTusoIsvd4/ixiJTIYj+mkLNEYfKiSI9qpznutH1AB1W9BOgNTMhgvQHA66raHOdCHet219AbaONOTwZuz2L/1wO/ikgEMAPorapNcHoyuF9EygI3Ao1UtSnwgufKqjobiMH55d9cVU95zJ4N3OTxvjfwYQ7j7IzTTUeqEaoaDTQF2olIU1WdgNOXz1WqepXblcdI4Fr3WMYAQ7LYjynk8mQXHqbQO+VeLD0VASa6dfLJOP0WpbccGCEiVYBPVXWziFwDtABWud2bROIknYz8R0ROAdtxuqGuB2xT1T/c+f8HPABMxBnr4m0RmQ/43KW5qu4Xka1uPzub3X384G43O3EWx+muwnOEsltEpD/O9/oCnAF61qVbt5U7/Qd3P0VxjpsxmbJEYfKLR4C9QDOckvA5gxKp6gci8hPQFVgkIvfidKv8f6o63Id93O7ZgaCIZDi+idu3UEucTub6AIOAq7PxWT4EbgE2AXNUVcW5avscJ84obi8Bk4CbRKQG8BhwqaoeFpEZOB3fpSfAYlW9NRvxmkLOqp5MflEK2OOOH9AX59f0WUSkJrDVrW6Zi1MF8xXQS0QqusuUFd/HFN8EVBeR2u77vv/f3t2jRBAEARR+lRpoaGjqEQRPYGZm6CW8hqksBmImIhiJggay2WLk717CQEQWBMEyqFmDZbZlQ+F9YTP09EzQRVc13cCwy+mvZOYlVSju23n0QR173ucc2KbuSDjt2hYaZ2Z+USmkjS5ttQxMgPeIWAW25oxlBGxOvykiliKib3Um/TJQ6L84AHYjYkSlnSY9z+wAzxFxD6xTVz6OqQn1OiIegRsqLfOnzPykTtc8i4gn4BsYUJPuRdffkFrtzDoGBtNi9ky/b8AYWMvMu65t4XF2tY99YC8zH6j7sV+AIyqdNXUIXEXEbWa+UjuyTrr3jKh/Jc3l6bGSpCZXFJKkJgOFJKnJQCFJajJQSJKaDBSSpCYDhSSpyUAhSWr6Ac+aLxDrkWCuAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fdb64204668>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot_roc(true_values, dlscore_values, 'dlscore', 'comt')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Enrichment factor analysis"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def EF(actives_list, score_list, n_percent):\n",
" \"\"\" Calculates enrichment factor.\n",
" Parameters:\n",
" actives_list - binary array of active/decoy status.\n",
" score_list - array of experimental scores.\n",
" n_percent - a decimal percentage.\n",
" \"\"\"\n",
" total_actives = len(actives_list[actives_list == 1])\n",
" total_compounds = len(actives_list)\n",
" # Sort scores, while keeping track of active/decoy status\n",
" # NOTE: This will be inefficient for large arrays\n",
" labeled_hits = sorted(zip(score_list, actives_list), reverse=True)\n",
" # Get top n percent of hits\n",
" num_top = int(total_compounds * n_percent)\n",
" top_hits = labeled_hits[0:num_top] \n",
" num_actives_top = len([value for score, value in top_hits if value == 1])\n",
" # Calculate enrichment factor\n",
" return num_actives_top / (total_actives * n_percent)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"NNScore enrichment factor at 0.1: 3.765\n",
"DLScore enrichment factor at 0.1: 1.412\n",
"\n",
"NNScore enrichment factor at 0.01: 1.176\n",
"DLScore enrichment factor at 0.01: 0.000\n",
"\n"
]
}
],
"source": [
"for n in [0.1, 0.01]:\n",
" ef_nnscore = EF(true_values, nnscore_values, n)\n",
" ef_dlscore = EF(true_values, dlscore_values, n)\n",
" print(\"NNScore enrichment factor at {}: {:.3f}\".format(n, ef_nnscore))\n",
" print(\"DLScore enrichment factor at {}: {:.3f}\".format(n, ef_dlscore))\n",
" print(\"\")"
]
}
],
"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.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment