-
-
Save yannabraham/5f210fed773785d8b638 to your computer and use it in GitHub Desktop.
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Dose Response Curve Fitting in Python" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Setting up the stage" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"While I would naturally fit dose-response curves in R using *drc*, recently I have started using iPython more and more and was wondering : \"Can I do curve fitting in Python directly?\"\n", | |
"\n", | |
"Turns out this is quite easy, although not necessarily well documented. Let's create a toy example and see how things are done using *scipy.optimize* & *pandas*" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import pandas as pd\n", | |
"import numpy as np\n", | |
"import scipy.optimize as opt\n", | |
"import matplotlib.pyplot as plt\n", | |
"import seaborn as sns\n", | |
"%matplotlib inline" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Let's use a 4-parameter sigmoidal function as a starting point:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def ll4(x,b,c,d,e):\n", | |
" '''This function is basically a copy of the LL.4 function from the R drc package with\n", | |
" - b: hill slope\n", | |
" - c: min response\n", | |
" - d: max response\n", | |
" - e: EC50'''\n", | |
" return(c+(d-c)/(1+np.exp(b*(np.log(x)-np.log(e)))))\n", | |
"\n", | |
"def pDose(x):\n", | |
" '''This is just a helper function, to compute easily log transformed concentrations used in drug discovery'''\n", | |
" return(-np.log10(1e-6*x))" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Generating data" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Lets generate a set of basic curves using a range of parameters" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"params = [{'compound':'A', 'b':1, 'c':0, 'd':100, 'e':0.4,'startDose':10, 'nDose':8, 'dilution':3},\n", | |
" {'compound':'B', 'b':0.7, 'c':0, 'd':86, 'e':1.3,'startDose':30, 'nDose':8, 'dilution':3},\n", | |
" {'compound':'C', 'b':2, 'c':24, 'd':152, 'e':0.02,'startDose':3, 'nDose':8, 'dilution':3}]\n", | |
"paramsCompound = [item['compound'] for item in params]" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"replicates are generated by adding random noise to the reference curve:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"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>compound</th>\n", | |
" <th>dose</th>\n", | |
" <th>logDose</th>\n", | |
" <th>response</th>\n", | |
" </tr>\n", | |
" </thead>\n", | |
" <tbody>\n", | |
" <tr>\n", | |
" <th>0</th>\n", | |
" <td>A</td>\n", | |
" <td>10.000000</td>\n", | |
" <td>5.000000</td>\n", | |
" <td>5.393205</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>1</th>\n", | |
" <td>A</td>\n", | |
" <td>3.333333</td>\n", | |
" <td>5.477121</td>\n", | |
" <td>12.261337</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>2</th>\n", | |
" <td>A</td>\n", | |
" <td>1.111111</td>\n", | |
" <td>5.954243</td>\n", | |
" <td>28.017640</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>3</th>\n", | |
" <td>A</td>\n", | |
" <td>0.370370</td>\n", | |
" <td>6.431364</td>\n", | |
" <td>53.470128</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>4</th>\n", | |
" <td>A</td>\n", | |
" <td>0.123457</td>\n", | |
" <td>6.908485</td>\n", | |
" <td>77.962146</td>\n", | |
" </tr>\n", | |
" </tbody>\n", | |
"</table>\n", | |
"</div>" | |
], | |
"text/plain": [ | |
" compound dose logDose response\n", | |
"0 A 10.000000 5.000000 5.393205\n", | |
"1 A 3.333333 5.477121 12.261337\n", | |
"2 A 1.111111 5.954243 28.017640\n", | |
"3 A 0.370370 6.431364 53.470128\n", | |
"4 A 0.123457 6.908485 77.962146" | |
] | |
}, | |
"execution_count": 4, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"drData=[]\n", | |
"for curve in params:\n", | |
" # generate base curve\n", | |
" curData = pd.DataFrame(data={'compound':curve['compound'],\n", | |
" 'dose':curve['startDose']/np.power(curve['dilution'],range(curve['nDose']))})\n", | |
" curData['logDose'] = pDose(curData.dose)\n", | |
" curData['response'] = curData.dose.apply(lambda x: ll4(x,*[curve[i] for i in ['b','c','d','e']]))\n", | |
" # generate replicates\n", | |
" repData = []\n", | |
" for i in range(5):\n", | |
" rep = curData\n", | |
" rep.response += 0.25*np.random.normal(len(rep.response))\n", | |
" repData.append(rep.copy())\n", | |
" repData = pd.concat(repData)\n", | |
" drData.append(repData)\n", | |
"# assemble data\n", | |
"drData = pd.concat(drData)\n", | |
"drData.head()" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Fitting the curve" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"is as simple as" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 5, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stderr", | |
"output_type": "stream", | |
"text": [ | |
"/tmp/ipykernel_1567/3180035174.py:7: RuntimeWarning: invalid value encountered in log\n", | |
" return(c+(d-c)/(1+np.exp(b*(np.log(x)-np.log(e)))))\n" | |
] | |
} | |
], | |
"source": [ | |
"compoundData = drData.groupby(['compound'])\n", | |
"fitData = []\n", | |
"for name,group in compoundData:\n", | |
" fitCoefs, covMatrix = opt.curve_fit(ll4, group.dose, group.response)\n", | |
" resids = group.response-group.dose.apply(lambda x: ll4(x,*fitCoefs))\n", | |
" curFit = dict(zip(['b','c','d','e'],fitCoefs))\n", | |
" curFit['compound']=name\n", | |
" curFit['residuals']=sum(resids**2)\n", | |
" fitData.append(curFit)\n", | |
"fitCompound = [ item['compound'] for item in fitData]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 6, | |
"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>b</th>\n", | |
" <th>c</th>\n", | |
" <th>d</th>\n", | |
" <th>e</th>\n", | |
" <th>residuals</th>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>compound</th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" </tr>\n", | |
" </thead>\n", | |
" <tbody>\n", | |
" <tr>\n", | |
" <th>(A,)</th>\n", | |
" <td>1.0</td>\n", | |
" <td>5.756833</td>\n", | |
" <td>105.756833</td>\n", | |
" <td>0.40</td>\n", | |
" <td>354.207787</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>(B,)</th>\n", | |
" <td>0.7</td>\n", | |
" <td>5.653138</td>\n", | |
" <td>91.653138</td>\n", | |
" <td>1.30</td>\n", | |
" <td>279.283932</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>(C,)</th>\n", | |
" <td>2.0</td>\n", | |
" <td>29.441756</td>\n", | |
" <td>157.441756</td>\n", | |
" <td>0.02</td>\n", | |
" <td>283.328028</td>\n", | |
" </tr>\n", | |
" </tbody>\n", | |
"</table>\n", | |
"</div>" | |
], | |
"text/plain": [ | |
" b c d e residuals\n", | |
"compound \n", | |
"(A,) 1.0 5.756833 105.756833 0.40 354.207787\n", | |
"(B,) 0.7 5.653138 91.653138 1.30 279.283932\n", | |
"(C,) 2.0 29.441756 157.441756 0.02 283.328028" | |
] | |
}, | |
"execution_count": 6, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"fitTable = pd.DataFrame(fitData).set_index('compound')\n", | |
"fitTable" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 7, | |
"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>b</th>\n", | |
" <th>c</th>\n", | |
" <th>d</th>\n", | |
" <th>e</th>\n", | |
" <th>residuals</th>\n", | |
" <th>ref_b</th>\n", | |
" <th>ref_c</th>\n", | |
" <th>ref_d</th>\n", | |
" <th>ref_e</th>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>compound</th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" <th></th>\n", | |
" </tr>\n", | |
" </thead>\n", | |
" <tbody>\n", | |
" <tr>\n", | |
" <th>(A,)</th>\n", | |
" <td>1.0</td>\n", | |
" <td>5.756833</td>\n", | |
" <td>105.756833</td>\n", | |
" <td>0.40</td>\n", | |
" <td>354.207787</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>(B,)</th>\n", | |
" <td>0.7</td>\n", | |
" <td>5.653138</td>\n", | |
" <td>91.653138</td>\n", | |
" <td>1.30</td>\n", | |
" <td>279.283932</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>(C,)</th>\n", | |
" <td>2.0</td>\n", | |
" <td>29.441756</td>\n", | |
" <td>157.441756</td>\n", | |
" <td>0.02</td>\n", | |
" <td>283.328028</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" <td>NaN</td>\n", | |
" </tr>\n", | |
" </tbody>\n", | |
"</table>\n", | |
"</div>" | |
], | |
"text/plain": [ | |
" b c d e residuals ref_b ref_c ref_d \\\n", | |
"compound \n", | |
"(A,) 1.0 5.756833 105.756833 0.40 354.207787 NaN NaN NaN \n", | |
"(B,) 0.7 5.653138 91.653138 1.30 279.283932 NaN NaN NaN \n", | |
"(C,) 2.0 29.441756 157.441756 0.02 283.328028 NaN NaN NaN \n", | |
"\n", | |
" ref_e \n", | |
"compound \n", | |
"(A,) NaN \n", | |
"(B,) NaN \n", | |
"(C,) NaN " | |
] | |
}, | |
"execution_count": 7, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"fitTable = pd.DataFrame(fitData).set_index('compound')\n", | |
"paramTable = pd.DataFrame(params).set_index('compound')[['b','c','d','e']]\n", | |
"paramTable.columns = ['ref_'+i for i in paramTable.columns]\n", | |
"fitTable.join(paramTable)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The model can be plotted as follows:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 8, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHpCAYAAACLCpbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACgkElEQVR4nOzdeXxU5dXA8d9dZslkmWQSsgnIDjEiboggLiitK1WhrbZq1VK3qq3aV6u20tZat7ZqrRZapVpbbau1WlfcBZVNUBAxsgZZQhKyTpbZ7vL+cTNDwpqEJJOE833f+czMvTNzn6SSnDzPec5RbNu2EUIIIYToh9RkD0AIIYQQortIoCOEEEKIfksCHSGEEEL0WxLoCCGEEKLfkkBHCCGEEP2WBDpCCCGE6Lck0BFCCCFEvyWBDmDbNsFgECkpJIQQQvQvEugADQ0N+P1+Ghoakj0UIYQQQnQhCXSEEEII0W9JoCOEEEKIfksCHSGEEEL0WxLoCCGEEKLfkkBHCCGEEP2WBDpCCCGE6Lck0BFCCCFEvyWBjhBCCCH6LQl0hBBCCNFvSaAjhBBCiH5LAh0hhBBC9FsS6AghhBCi35JARwghhBD9lp7sAQghhOhbLNuipKaEunAdmd5MigJFqIr83Sx6Jwl0hBBCtNuS7UuYu2oupcFSDMtAV3WGZgxl5tiZTCiYkOzhCbEbCcGFEEK0y5LtS7hz0Z2srV2LT/eRk5KDT/extnYtdy66kyXblyR7iELsJqmBzoIFC5g2bRqFhYUoisKLL76422tKSkr4xje+gd/vJzU1lfHjx7N58+bE+XA4zLXXXkt2djZpaWnMmDGDioqKHvwqhBCi/7Nsi7mr5tIUayLXl4tX96IqKl7dS64vl6ZYE3NXzcWyrWQPVYg2khroNDU1MW7cOB599NE9nt+wYQOTJ09mzJgxvP/++3z22WfccccdeL3exGtuvPFGXn75ZZ577jnmz59PWVkZ06dP76kvQQghDgolNSWUBkvxe/woitLmnKIo+D1+SoOllNSUJGmEQuxZUnN0zjzzTM4888y9nv/Zz37GWWedxf333584Nnz48MTj+vp65s6dyzPPPMOpp54KwBNPPEFRURGLFy/m+OOP777BCyHEQaQuXIdhGbg1NwBhI5zI0fHqXtyam2A0SF24LrkDFWIXvTZHx7IsXn31VUaNGsXpp59Obm4uEyZMaLO8tXz5cmKxGFOnTk0cGzNmDIMHD2bRokV7/exIJEIwGGxzE0IIsXeZ3kx0Vac+Us9Xwa/Y3LCZssYyNjds5qvgV9RH6tFVnUxvZrKHKkQbvTbQqayspLGxkXvvvZczzjiDN998k/PPP5/p06czf/58AMrLy3G73WRmZrZ5b15eHuXl5Xv97HvuuQe/35+4DRo0qDu/FCGE6POKAkUEPAEqmioIxUKoioqu6qiKSigWoqKpgoAnQFGgKNlDFaKNXhvoWJaT0Hbuuedy4403cuSRR3LrrbdyzjnnMGfOnAP67Ntuu436+vrEbcuWLV0xZCGE6P8UJycHG2zbBrvlubL/twqRDL020MnJyUHXdQ477LA2x4uKihK7rvLz84lGo9TV1bV5TUVFBfn5+Xv9bI/HQ0ZGRpubEEKIvSupKaEmUkOeLw+v7sW0TWJ2DNM28epe8nx51ERqJBlZ9Dq9NtBxu92MHz+eNWvWtDm+du1aDj30UACOOeYYXC4X77zzTuL8mjVr2Lx5MxMnTuzR8QohRH8WT0bWVd2ZyQGI39k2uqpjWIYkI4teJ6m7rhobG1m/fn3ieWlpKStWrCAQCDB48GBuvvlmLrjgAk466SSmTJnCvHnzePnll3n//fcB8Pv9zJw5k5tuuolAIEBGRgbXX389EydOlB1XQgjRhTK9mVi2xfbG7djYaKqGioqFRcSMsL1xOxmeDElGFr2OYidC8573/vvvM2XKlN2OX3rppTz55JMA/PWvf+Wee+5h69atjB49ml/96lece+65ideGw2F+8pOf8M9//pNIJMLpp5/On/70p30uXe0qGAzi9/upr6+XZSwhRJ/Q0/2mDMvglH+fQjAaxK2629TSsW2bqBUlw53B+xe8j67qWLZFY6yRhmhD4hY2wkTMCGEzTMRouTcjicdRM4phGZi2iWmbGJaReG5YBqZlYtg7j8Wf27aNbdtYWInZJhsby7awsRPnE/9nt72PFzls/dzGyT+yW6atLhh9AVeNu6rbvr+i+yQ10OktJNARQvQlyeg3tbp6NVe/dTUN0QZs204EOvHgAEBVVLK8WYSNME2xpsTx/uDywy/npmNuSvYwRCdIoIMEOkKIviPeb6op1oTf48etuYmaUeoj9aS6Upk1cdYBBzsxM8bG+o18WfMlX9Z8yabgJtbVrqOiuePtdTyah3R3OmmuNFL0FDyaB4/uwat58WgevLpzH7/pqo6marhUF5qioSkauqo7x1sea6qGruiJ7e3xG4CCgqIoqKjODjEUVEVNHAfaPG9z3/La+Ptan8vyZpHryz2g76tIDuleLoQQfcSu/abiv7jjwUJlcyVzV81lfP74Di1j1YXrWF6xnKXlS/mk8hPW163HsIy9vt6tOtWRFUVBV3VS9VRsbAzL4CfH/oRxA8aR7k4n3Z2eqKQsRLJIoCOEEH1ER/pNFWcX7/VzbNvm86rPeWvzWyzctpC1tWt3W2ZKd6UzOjCaMYExDMscxsC0gcxZMYd1deswbZOoFQULTNtM7Loqzi5m2vBp3ZorJERHSaAjhBB9xIH2myqtL+WFdS/wWulruy1DDfcPZ3z+eI7NP5axOWMpSC3YLZhaU7OGz6o+w7RNNNVZVjIxCRkhNEVjUuEkCXJEryOBjhBC9BGt+00Fo0EiZsSpZaM4uTAZ7ozd+k1ZtsU7m9/h6ZKnWV6xPHHcp/s4aeBJTBk0heMKjiMnJWef17Zsi4VlC/HpPgzbIGpFMW0TFEjRU9AVnYVlC/le8fck2BG9igQ6QgjRR8T7TX1Z8yUAuuYk41pYhGIhQrEQYwJjKAoUYVgGr5W+xuOrHqe0vhRwknBPPOREzh95PpMPmYxH87T72vFlsxxfDh7NQzAaJGpGcWtuMtwZRMxIu5bNhOhpEugIIURf07IryLKdujG24mz3jufZLCpbxG8//i0b6jcAkO5O5ztjvsO3Rn2L/NT21xhrLb5sFrNiVDRXtJlNqovUEfAGpDKy6JUk0BFCiD4i3m8q05NJXaQO0zKdE7YzW5PuTmdD/QaufvtqALI8WVxafCkXjL6ANHfaAV17j5WRW2aTwkZYKiOLXksCHSGE6CPqwnU0x5oJxUJgg67qKDgzOaZlUh+pB5yg57tjvsvV467G7/F3ybVHZ43GtJyKxW7VjdLSrlxFRVEUJ2fHMhmdNbpLridEV5FARwgh+gi/x0/YCDvBRsvOK9u2idmxxLKVgsJ9J97HGUPP6NJrr6ldg6qqaKqGYRtoaDuDrJZdWKqqsqZ2jeToiF5FUuOFEKKPSAQziuL0drItolY00atJQUFXdA5JP6TLr10XrkNTNApSC/DqXiwsTNvEwsKreylILUBTNMnREb2OzOgIIUQfEYwE8epeQkaImBXDwkqcU3FmW1L0FIKRYJdfO7613aW6GJw+mLAZxrScmRyv5iViRoipMcnREb2OzOgIIUQfkenNxOfykeZK2y3ISXGlkJOSg8/l65ZgoyhQxNCMoYk8oBQ9hTS3078KoD5Sz9CMoRQFirr82kIcCAl0hBCijygKFJHlyaI2Ups4Fm88aVkWTbGmbgs2VEVl5tiZpLpSqWyuJGyEsWxnx1VlcyWprlRmjp0pxQJFryP/RQohRB+xvGI562rXJZ5riubsgFIUwmaY5lhzt7ZhmFAwgVkTZzEqaxTNRjNVoSqajWZGZY3qkq7pQnQHxbZte/8v69+CwSB+v5/6+noyMjKSPRwhhNjN51Wf8/03vk/ICKErTq5MzI4liva5VTe6olOcU8ycr83p1pkVy7YoqSmhLlxHpjeTokCRzOSIXkuSkYUQopfb0byDH7/7Y0JGCLfqpjCt0Ek6TlIbBlVRZQu56DMk0BFCiF4saka58f0bqQxVkp+aT9R0mmlubtgsbRiEaAeZaxRCiF7Ktm3uXnI3K3esJN2Vzi3jbwFge+N2wkYYVVHRVaexZ7wNg2VbssVbiFYk0BFCiF7qlY2v8Py651FQuP/k+5kyaEqiDYOu6KgtP8JVVHRFx7RNacMgxC4k0BFCiF6ovKmce5bcA8APj/whkw+ZvFsbhnj3csu2nLYMrdowCCEcEugIIUQvY9s2v1z4SxpiDYzNGcsPxv4AkDYMQnSGJCMLIUQv8/y65/mo7CPcqpu7TrgLXXV+VEsbBiE6TmZ0hBCiFylvKue3H/8WgB8d/SOGZQ5LnJM2DEJ0nAQ6QgjRizyw/AGajWaOHHAkFxdd3OactGEQouPkX4MQQvQSKypX8Hrp6ygo3DbhNjRV2+010oZBiI6RHB0hhOgFLNtKLFmdN+I8Dss+bK+vnVAwgfH546UNgxDtIIGOEEL0Aq+Xvs5nVZ+Roqdw/VHX7/f10oZBiPaR8F8IIZIsbIR5cPmDAFwx9goG+AYkeURC9B8S6AghRJI9v+55KporKEgt4JLDLkn2cIToVyTQEUKIJIqaUf76+V8B+MHYH+DVvUkekRD9iwQ6QgiRRC9teInK5kpyU3I5b8R5yR6OEP2OBDpCCJEkhmUwd9VcAC47/DLcmjvJIxKi/5FARwghkuT10tfZ2riVLE8WM0bOSPZwhOiXJNARQogksGwrMZtzyWGX4HP5kjwiIfonCXSEECIJPtr2ERvqN5DmSuPCMRcmezhC9FsS6AghRBI8u+ZZwKmCnO5OT/JohOi/JNARQogetr1xOwu2LQDgW6O/leTRCNG/SaAjhBA97Lm1z2HZFhPyJzDMPyzZwxGiX0tqoLNgwQKmTZtGYWEhiqLw4osv7vW1V199NYqi8NBDD7U5XlNTw0UXXURGRgaZmZnMnDmTxsbG7h24EEJ0UsyM8fy65wH49uhvJ3k0QvR/SQ10mpqaGDduHI8++ug+X/fCCy+wePFiCgsLdzt30UUXsXr1at566y1eeeUVFixYwJVXXtldQxZCiAPyzuZ3qAnXMCBlAFMGT0n2cITo95LavfzMM8/kzDPP3Odrtm3bxvXXX88bb7zB2Wef3eZcSUkJ8+bN4+OPP+bYY48F4I9//CNnnXUWv/vd7/YYGAkhRDL9e82/AZgxagYu1ZXk0QjR//XqHB3Lsrjkkku4+eabKS4u3u38okWLyMzMTAQ5AFOnTkVVVZYsWbLXz41EIgSDwTY3IYTobpuDm1lWsQxVUaVAoBA9pFcHOvfddx+6rvOjH/1oj+fLy8vJzc1tc0zXdQKBAOXl5Xv93HvuuQe/35+4DRo0qEvHLYQQe/LqxlcBmFgwkfzU/CSPRoiDQ68NdJYvX84f/vAHnnzySRRF6dLPvu2226ivr0/ctmzZ0qWfL4QQu7Jtm9dKXwPg7GFn7+fVQoiu0msDnQ8++IDKykoGDx6Mruvous5XX33FT37yE4YMGQJAfn4+lZWVbd5nGAY1NTXk5+/9ryWPx0NGRkabmxBCdKcvqr9gU3ATXs3LqYNPTfZwhDhoJDUZeV8uueQSpk6d2ubY6aefziWXXMLll18OwMSJE6mrq2P58uUcc8wxALz77rtYlsWECRN6fMxCCLE3r5Y6y1anDDqFVFdqkkcjxMEjqYFOY2Mj69evTzwvLS1lxYoVBAIBBg8eTHZ2dpvXu1wu8vPzGT16NABFRUWcccYZXHHFFcyZM4dYLMZ1113HhRdeKDuuhBC9hmmZzCudB8iylRA9LalLV8uWLeOoo47iqKOOAuCmm27iqKOOYtasWe3+jKeffpoxY8Zw2mmncdZZZzF58mT+8pe/dNeQhRCiw5aWL2VHaAd+j58TCk9I9nCEOKgkdUbnlFNOwbbtdr9+06ZNux0LBAI888wzXTgqIYToWvEk5K8f+nVcmtTOEaIn9dpkZCGE6A+iZpS3v3obkGUrIZJBAh0hhOhGS8uX0hhrZEDKAI7KPSrZwxHioCOBjhBCdKP3Nr8HOLutVEV+5ArR0+RfnRBCdBPLtnh/y/sATBkkDTyFSAYJdIQQopt8Uf0FlaFKfLqPCQVS20uIZJBARwghusl7W5xlqxMOOQG35k7yaIQ4OEmgI4QQ3SQe6MiylRDJI4GOEEJ0g60NW1lXuw5N0Thp4EnJHo4QBy0JdIQQohvEZ3OOzjsav8ef5NEIcfCSQEcIIbqBLFsJ0TtIoCOEEF2sMdrIJxWfAE79HCFE8kigI4QQXWxp+VJM2+TQjEMZlD4o2cMR4qAmgY4QQnSxhWULAZhYMDHJIxFCSKAjhBBdbPH2xQBMLJRAR4hkk0BHCCG60LbGbXwV/ApN0Tgu/7hkD0eIg54EOkII0YUWlS0C4IgBR5DmTkvyaIQQEugIIUQXSuTnyLKVEL2CBDpCCNFFTMtkyfYlgCQiC9FbSKAjhBBd5IvqLwhGg6S70jk85/BkD0cIgQQ6QgjRZRZtd/JzJhRMQFf1JI9GCAES6AghRJeR/Bwheh8JdIQQoguEjTCf7fgMgOMLjk/yaIQQcRLoCCFEF1hVtYqYFSM3JVfaPgjRi8gishBCdIJlW5TUlFAXriPTm8nH5R8DcEzeMSiKkuTRCSHiJNARQogOWrJ9CXNXzaU0WIphGeiqTigWApxARwjRe0igI4QQHbBk+xLuXHQnTbEm/B4/bs1NxIhQ3lQOgEtzJXmEQojWJEdHCCHaybIt5q6aS1OsiVxfLl7di6qo2NgAKCi8vvF1LNtK8kiFEHES6AghRDuV1JRQGizF7/G3ycNpNpoBSNFT2NSwiZKakmQNUQixC1m6EkKIdqoL12FYBm7NDThbyg3LoCHaAECaOw3DMqgL1yVxlEKI1iTQEUKIdsr0ZqKrOvWReoLRIBEzgm3ZmJgAWJaFrulkejOTO1AhRIIsXQkhRDsVBYoIeAJUNFUQioVQFRVN1RLnq8PVBDwBigJFSRylEKI1mdERQoiOUpzEY2wwbXPnYamfI0SvIzM6QgjRTiU1JdREasjz5eHVvZi2mQh0XKqLPF8eNZEaSUYWoheRQEcIIdopnoysqzq2bWPbduKcqqjoqi7JyEL0MhLoCCFEO2V6M7Fsi+2N24mYkTb5OTEzxvbG7Vi2JcnIQvQiEugIIUQ7jc4ajWk5y1W6otNSJxAFBZfqcpayLJPRWaOTO1AhRIIEOkII0U5rategqs5OK8M2Evk5CgqGbaCpGqqqsqZ2TZJHKoSIk0BHCCHaqS5ch6ZoFKQW4NW9WLS0elDAq3spSC1AUzTJ0RGiF5FARwgh2ileMNCluhiYNjBxfEDKAAanD8alutBVKRgoRG+S1EBnwYIFTJs2jcLCQhRF4cUXX0yci8Vi/PSnP2Xs2LGkpqZSWFjI9773PcrKytp8Rk1NDRdddBEZGRlkZmYyc+ZMGhsbe/grEUIcDIoCRQzNGEpVqIpNwU2J41WhKr4KfkVVqIqhGUOlYKAQvUhSA52mpibGjRvHo48+utu55uZmPvnkE+644w4++eQT/vvf/7JmzRq+8Y1vtHndRRddxOrVq3nrrbd45ZVXWLBgAVdeeWVPfQlCiIOIqqhMKpxEc6yZiBlxjqGiKAohI0RzrJlJhZNQFZksF6K3UOzWhSCSSFEUXnjhBc4777y9vubjjz/muOOO46uvvmLw4MGUlJRw2GGH8fHHH3PssccCMG/ePM466yy2bt1KYWFhu64dDAbx+/3U19eTkZHRFV+OEKIfsmyLq9+6mtVVq2kymjBtExUVVVVxq250Rac4p5g5X5sjwY4QvUSf+pdYX1+PoihkZmYCsGjRIjIzMxNBDsDUqVNRVZUlS5bs9XMikQjBYLDNTQgh9qekpoTSYCk5vhynBQSQ7k4nz5fHoRmHkuPLoTRYKpWRhehF+kygEw6H+elPf8p3vvOdxKxLeXk5ubm5bV6n6zqBQIDy8vK9ftY999yD3+9P3AYNGtStYxdC9A/xysghI4RhGwA0RhupaK5gc8NmYlZMKiML0cv0iUAnFovx7W9/G9u2mT179gF/3m233UZ9fX3itmXLli4YpRCiv4tXRq5oqgBaCgVqLlRFJWyEpTKyEL1Qr+9eHg9yvvrqK9599902OTT5+flUVla2eb1hGNTU1JCfn7/Xz/R4PHg8nm4bsxCif4pXRo7Xz4nn4cQTkqNWVCojC9HL9OoZnXiQs27dOt5++22ys7PbnJ84cSJ1dXUsX748cezdd9/FsiwmTJjQ08MVQvRz8crIrdm2jWVbUhlZiF4qqTM6jY2NrF+/PvG8tLSUFStWEAgEKCgo4Jvf/CaffPIJr7zyCqZpJvJuAoEAbreboqIizjjjDK644grmzJlDLBbjuuuu48ILL2z3jishhGivunCdM3uDgt3S6Mq0zURl5IA3QMgISY6OEL1IUgOdZcuWMWXKlMTzm266CYBLL72UX/7yl7z00ksAHHnkkW3e995773HKKacA8PTTT3Pddddx2mmnoaoqM2bM4OGHH+6R8QshDi6Z3kxQwMZGRWVQ+iAs20JTNbyal4gZIabGJEdHiF4kqYHOKaecwr7K+LSnxE8gEOCZZ57pymEJIcQeFQWKyHRnUh2qxqt78bl8iXO2bVMfqWdU1iipjCxEL9Krc3SEEKI3URWVwRmDAad4YNgIJ+4rmytJdaUyc+xMKRYoRC/S63ddCSFEb1IVqgJgUPogGmINBKNBdFVnVNYoZo6dyYQC2QghRG8igY4QQrRTzIyxpsbZUfXglAdpjDVSF64j05tJUaBIZnKE6IUk0BFCiHZaX7eeqBUl3Z3O4PTBKIqS7CEJIfZD/vwQQoh2Wl29GoDi7GIJcoToIyTQEUKIdmod6Agh+gYJdIQQop1WV7UEOjkS6AjRV0igI4QQ7RAxI6yrXQfIjI4QfYkEOkII0Q5ra9Zi2AZZniwKUguSPRwhRDtJoCOEEO0Qz885LOcwSUQWog+RQEcIIdpBEpGF6Jsk0BFCiHaIFwo8LHBYkkcihOgICXSEEGI/YlaM9XXrARgVGJXk0QghOkICHSGE2I9N9ZuIWTFSXakcknZIsocjhOgACXSEEGI/1tQ6y1YjM0dKPysh+hj5FyuEEPuxtmYtAKMDo5M8EiFER0lTTyGE2I/4jM6orOTm51iWzeqyIDXNUQI+N8WFGaiqbHUXYl8k0BFCiP1YW5v8GZ2F66uYPX8DGyobiZk2Lk1heG4a15w8nEkjcpI2LiF6O1m6EkKIfagOVVMVqkJBYWTmyKSMYeH6Km5/YRUl24OkenRy0z2kenRKtjdw+wurWLi+KinjEqIvkEBHCCH2Ib5sNThjMD6Xr8evb1k2s+dvoDFikJ/hxevSUFUFr0sjP8NDY8Rk9vwNWJbd42MToi+QQEcIIfYhnoicrPyc1WVBNlQ2kuVz79Z6QlEUMn0uNlQ2srosmJTxCdHbSY6OEELsQ+tE5GQkA9c0R4mZNm5NxcYmHLUwLAtdVfG6VTyaSr1lU9Mc7dZxCNFXSaAjhBD7EA907HABlz6xtMeTgQM+Ny5NoS4Uoz4UI2KY2DYoCnh0DX+KC5eqEPC5u20MQvRlsnQlhBB7ETNjlNaVAvDMh0ZSkoGLCzPITnOzvT5EKGqgKgq6pqAqCqGowfb6ENlpzuySEGJ3EugIIcRebKzfiGEbaLaP5lB68pOB46tk9i7PhRB7JYGOEELsRXzZimgBAZ8nKcnAq8uCVDdGKfB7SXFpWLaNYdlYtk2KSyPf76W6MSrJyELsheToCCHEXqypiQc6hUlLBo4nI+eme8j0uahvNoiZFi5Nxe/TwVaobIxIMrIQeyGBjhBC7EV8Rkc3C5OWDLyvZOS6UO9NRrYsm6hpETMtooZFzLQxLAvTapmRark3W25G4t7Cski8tvV5y7YxzJZjtnNv2zaWTeLesp11Pcu2sW3aHrNaXsvO99gt5+PH7V0+K/78hBE5fL04P5nfUtFJEugIIcQe2LadqKET0A9lfX0IBdA1FUUF24ZQ1KA5alBcmNFtycDxZOTVZcEuuX44ZtIQNmgIxwiGDZoiBqGoSXPMJBw1aY4ahGIWoahBKGbSHDUJxUxCLfeRmNUmgIkaFhGj5XkiqHECm/7E69Yk0OmjJNARQog9qApVURupRVVUXFYhEN5zMnBP/j5vdX3btrFxAp6GsMH/VmyjpjlGTVOEmqYo9aEYwZAT0DSEDYItgU3UsHpwwDvpqrNbTFMUNFVB11RnB5nqPNfUto/jz9X4fctuM01V0RTQVBVVAU1VUBQnX0pVFBRAbXmuKLQ6pqCqzjdRjR9vuSd+XmHnexLnnefHDgkk5fsmDpwEOkIIsQfxZav8lEHUboMCv7dl6chKLB2luDQyUlyJZOCxA/1dcu365hjb6kKUB0Ms21RLaVUTXpdGJGYS2UOgsqm6mRufXdmha6R7dNK9OmlenRS3TopLxefWSXFppLg1UlwaPreGt+U+peWx16Xh1lTcuoJb03DrKi5Nwa2reHQVl6bi1tWW17Q811Tpsi6SRgIdIYTYg3gicn7KMGpakoGzfG7CsVbJyC4V26bDycC2bVPdFOWr6iY2VTU799U77+tDsXZ/lq4q2MCovDSGD0gjO9VNINVJXM5I0Un3uEj36mSkOPfpXhdpHh1NAg9xkJBARwgh9iA+ozMycxTrNIWoaeFx7VKRQ4GIYe0zGbgpYrCmooE15c7ty/Iga8obqG3edzCTneom3+8l1a3zeVk9KS4Nr1vFbpnQ8egaPo9KxLBpjhjcP2Ncl80oCdGfSKAjhBB7EE9EPvHQI1izxsNnW+swTGcnUXzpyq2p6JrCEQMzKS7MIGZafLm9gRVb61i5pY4VW+rYsKMRew95PIoChf4UDs32cWh2KkPi9zk+Bgd8+NzOj2fLsrn0iaV8trWOuiZzn9cXQuxOAh0hhNhFxIywKbgJgNGBUZw0solFG6oxLbslIdbZttwUNVEVJxj51p8X8fm2+j3m0ORleBiVl86Y/HRG52cwJj+dEblpeF3afseiqgonjczZ4/WboyZay3nJgRFizyTQEUKIXayvW49pm/g9fgZ4c1mw7mN8bg3DtIiYFqZFYpbGsuGjDdWJ9/pTXIwblMmRA/0cOTiTIwZmkpPm6fRYLMtmwboqfG4Ns6U2jWk5Mzo+t4amqixYV8XMycMk2BFiDyTQEUL0KZZls7osSE1zlIDPaWbZ1b/g48tWo7NG88nmOlZvqwcgYtqYu0zYaKqCW1O44sRhnHfUIQzNSd2tVcSBWF0WZENlI3kZXjy6ulsydNiwEi0oJEdHiN1JoCOE6DMWrq9i9vwNbKhsJGbauDSF4blpXHPycCaNyOmy65RUfwnAtopMLlqwhGir6EZVwOvSSHXrZPp0XKrKjqYoxwwJMGxAWpeNIS7eAsKtqbs38VTo9hYUQvR1EugIIfqEheuruP2FVTRGDLJ8btyaStS0KNnewO0vrOLu88ceULBj2zZLSmt44ZNtvFa1BFJg/bYMDNNCawluTNvGSFT/jRKKmUltAdETLSiE6OuS2r18wYIFTJs2jcLCQhRF4cUXX2xz3rZtZs2aRUFBASkpKUydOpV169a1eU1NTQ0XXXQRGRkZZGZmMnPmTBobG3vwqxBCdDfLspk9fwONEYP8DC9el4aqKnhdGvkZHhojJrPnb8CyOl6muK45yuMfbOS0B+Zz4V8W8+9lm7HdZQCcd9h4Xrr2BIoK0gnFTGKGhaY6O51URSEUNdheHyI7zd3tLSC214cIRY1EheCeur4QfV1SA52mpibGjRvHo48+usfz999/Pw8//DBz5sxhyZIlpKamcvrppxMOhxOvueiii1i9ejVvvfUWr7zyCgsWLODKK6/sqS9BCNED4nkqWT73bvkviqKQ6XMl8lTa64uyIDc9u4IJd7/DXa+WsHFHE6lujXOPSUXRQmiKxt3nTOXwQ/w7r7mnFhA9KdnXF6IPSurS1ZlnnsmZZ565x3O2bfPQQw/x85//nHPPPReAp556iry8PF588UUuvPBCSkpKmDdvHh9//DHHHnssAH/84x8566yz+N3vfkdhYWGPfS1CiO7TOk/FxiYcbZWQ61bbnadi2zaLNlYzZ/5GFqzdkTheVJDBxccP5twjD2F55Ue8+y4M9Q/FrblZtbWe6sZoj7aAaG11WTCp1xeir+u1OTqlpaWUl5czderUxDG/38+ECRNYtGgRF154IYsWLSIzMzMR5ABMnToVVVVZsmQJ559//h4/OxKJEIlEEs+Dwfb/FSiE6HkHmqdi2zZvl1TyyLvrWLnV2UGlKnD2EYVcfsIQjhqUmZi1WVvr7LgalTUK2Blk5aY7bRXqmw1ipoVLU/H7dLCVDreA6IjW1++qFhRCHEx6baBTXl4OQF5eXpvjeXl5iXPl5eXk5ua2Oa/rOoFAIPGaPbnnnnv41a9+1cUjFkJ0l3ieyuqyIAqgayqK6tSyCUUNmqMGxYUZe8xTWbyxmvvmfcmnm+sA8Ogq3z52EFecOIzB2b7dXh9v/TA6MBrYd5BVF+r+ZOD49aOmhbel4SbsLDQYNkxJRhZiH3ptoNOdbrvtNm666abE82AwyKBBg5I4IiFEu+0pT2UPOciry+q5f94a5rcsUXldKpdNGsoPThy6zwJ+8Waeo7OcQOdAgqyuUFyYwfDcNEq2N5CfobbJUbJtm7rmGEUF6ZKMLMRe9NpAJz8/H4CKigoKCgoSxysqKjjyyCMTr6msrGzzPsMwqKmpSbx/TzweDx5P5yuVCiF6VkfyVAYHfNz/xpc8s3Qztu10977wuEH86NSR5GZ493mdkBFic8NmYOeMThvtDLK6kqoqXHPycG5/YRXlwQiZPhceTSViWtQ1x0jzaFxz8nCpiizEXvTaQGfo0KHk5+fzzjvvJAKbYDDIkiVLuOaaawCYOHEidXV1LF++nGOOOQaAd999F8uymDBhQrKGLoToYu3JU6loCPO/ldt44ZNtVDc5+SrnHFHA/319NENyUtt1nfW167Fsi4A3QE6KU5OnNyQDTxqRw93nj00US6y3bFyqQlFBepcXSxSiv0lqoNPY2Mj69esTz0tLS1mxYgWBQIDBgwdzww03cNdddzFy5EiGDh3KHXfcQWFhIeeddx4ARUVFnHHGGVxxxRXMmTOHWCzGddddx4UXXig7roToR1rnqXhcu1TFUKAhbBAMxXj8g1IARuamcee5hzNxeHaHrhPPz4knIkPvSQaeNCKH44dld3v7CyH6m6QGOsuWLWPKlCmJ5/G8mUsvvZQnn3ySW265haamJq688krq6uqYPHky8+bNw+vdOf389NNPc91113HaaaehqiozZszg4Ycf7vGvRQjRfeJ5Kp9trcMwncaWTlNNG1VRiJrO+lGKS+PHU0fy/ROG4tY7XiYsvuMqnp8DvSsZWFUV2UIuRAclNdA55ZRTsO29L3ArisKdd97JnXfeudfXBAIBnnnmme4YnhCil1BVhZNG5rBoQzWmZbdUBraJmWC0JMkMG5DKP2ZOoDAzpdPXSSQit8rPkWRgIfq2pFZGFkKI9rAsmwXrqvC5NXxuDcuyiZo784C9ukqh30v+fpKN98W2bdbVOi1mWi9dxZOB0zwa5cEIoZiJZdmEYiblwYgkAwvRy0mgI4To9eItIAake3DrKi0rVXh0lZG5qQwM+Ni4o6lDLSB2VdZURkOsAV3VGeYf1uZcPBm4qCCd5ohBZWOE5ohBUUH6ATcTFUJ0r16760oIIeJqmqOEYxa1zVFCMQsAf4pOTqoHj0vDtmhXC4h9iS9bDfcPx6W5djsvycBC9E0S6Agher0tNc3UhaLEm5NrKjRFTJqjze1qAdEeu1ZE3hNJBhai75GlKyFEr/a3hZuY9b/PE0GOCuiq2pKQrBCKGmyvD5Gd5j6ghOC1NW17XAkh+geZ0RFC9EqmZfObV0v460dObZzMFBf1oRhK/M+zLq5OvGszTyFE/yAzOkKIXiccM/nh08sTQc73Jh5KikulMNNLikvDsm0My8aybVJcGvl+b6I6cWc0x5rZ0rAF2PfSlRCi75EZHSFErxIMx/jBk8tYuqkGt67y+2+NIyPFxeuryslNdZPpc1HfbBAzLVyait+ng60cUHXitbVrsbEZkDKAgDew9xdaFpSvhOZq8GVD/jhQ5e9FIXqzAwp0otEopaWlDB8+HF2XmEkIcWCqGyNc+sRSPt8WJN2jM/ey8Rw3NMCqrfW4NIW6UKyl35SZ6DdVFzrwZOTEslVgH8tWG+fDhw9C1TqwYqC6IGckTL4Rhp3cqesKIbpfp/4UaW5uZubMmfh8PoqLi9m82en2e/3113Pvvfd26QCFEAeHimCYb/95EZ9vC5Kd6uafVx7PcUOd2ZXiwgyy09xsrw8RihqoitKlyciJishZe1m22jgfXrkBKlaDOxXS8pz7itXO8Y3zO3VdIUT361Sgc9ttt7Fy5Uref//9Nn2npk6dyr///e8uG5wQ4uBQGQzznccWs2FHEwV+L89ePZHDD9nLNu542Rp7l+cHILG1fE+BjmU5MzmRRkgvAFcKKKpzn17gHP/wQed1Qohep1OBzosvvsgjjzzC5MmT2/R9KS4uZsOGDV02OCFE/7ejIcJ3HlvMxh1NHJKZwrNXTWT4gLQ2r1ldFqS6MUqBv+uTkS3b2mPrh4Tylc5yVUqWs1bWmqI4x6vWOa8TQvQ6nUqs2bFjB7m5ubsdb2pqahP4CCHEvtQ0Rbno8Z0zOf+84ngGBXy7v645Ssy0yU33kOVzE45ZGJaFrqp4XSq2TaeTkbc1bKPZaMatuhniH7L7C5qrnZwc3bPnD9A9EK5zXif6r3hSmOhzOhXoHHvssbz66qtcf/31AIng5vHHH2fixIldNzohRL/VGDG4/ImlrK1oJD/DCXIGZ+8e5AAEfG5cmkLUtPC6NFLcGqAlzocNs9PJyPFlq+GZw9HVPfxI9GU7icdGxFmuijWDZYKqgcvnHFddzuvE7mzb+R4ZIYiFne+fEQYzCmbMuTciOx+3Pt7mto/XWoZzs03nfxvLaHtvmztfY7V6jW3ufmxvn3PCj+Frv0r2d1N0QqcCnbvvvpszzzyTL774AsMw+MMf/sAXX3zBwoULmT9fkvKEEPsWMUyu/vtyVm6tJ8vn4h8/mMCQnNS9vr64MIPhuWmUbG8gP0NtM3Ns2zZ1zTGKCtI7lYy839YP+eOc3VVln4JpgBnBSRBSQPOApkPhUc7r+jLTgGiDk3MUaYBoI0SCzvNoy7FIy7FoE8RCLcFLq1simNnlXFdUdEw220z2CEQndSrQmTx5MitWrODee+9l7NixvPnmmxx99NEsWrSIsWPHdvUYhRD9iGXZ3PTsSj5cX4XPrfHE5ccxIjdtn+9RVYVrTh7O7S+sojwYIdPnwqOpREyLuuYYaR6Na04e3qkGm/vdcaWqMPxU2PSB89e9puP86DQh1gSm5pzvDfV0bNtZRgvVOrfm2p2P93YL1ztBjBHq/vEpLbNguse5ae5WN5dzr+/hmOZp9Tj+upZjqsu5VzRnlk3VW93rTuJ4/LHa6jVKq9eorV6jaHt+rXvvgbjo3Tpd/Gb48OE89thjXTkWIcRB4N55X/LqZ9txaQp/ueRYjhyU2a73TRqRw93nj2X2/A1sqGyk3rJxqQpFBelcc/JwJo3I6dR44jV09jqjY1mw4V1wpzlLGEYEMADF+eWn6s75idd1X7ATDkJDOTRWQFMlNLbcWj9urISmHU4+0YHQPOBJA086uNNbPW6596Q7X7fL5yzl6d6Wx96W5ynO/W7nfE5AIkQP61Sg88knn+ByuRKzN//73/944oknOOyww/jlL3+J2935DsJCiP7rmSWb+cuCjQD87lvjmDyyY8HJpBE5HD80iw2rFhKqqyQlM5fhY49F1bT9v3kPGqINbGvcBuyjx1V811VanvOL2wg5AY+qO7/UjfDOXVeFR3V8EEYUgluhfivUb3Pugy2Pgy3PIx3cTeZOc3aDpWS23O/j5vW3BDAZzvt0+fkt+pdOBTpXXXUVt956K2PHjmXjxo1ccMEFTJ8+neeee47m5mYeeuihLh6mEKKvW7B2B3f873MAbpw6inOPPKTjH7JxPuqHDzKydXXizzpfnTi+rTzPl4ffs5e6Pa13XSVyg1ruFaV9u65iYajdBDUbd7/VbwG7HTV4PH5Iy3UCrrQBkJrb8rzlWOqAlvucve8QE+Ig1KlAZ+3atRx55JEAPPfcc5x88sk888wzfPTRR1x44YUS6Agh2tiwo5Frn/4E07KZftQh/Oi0ER3/kHh14kijMxOhe5xlpHh14nMe6nCws99EZNi56ypUC6E6ZwYnnoyse51Zk/iuq0gj7FgDO0qgsuW2Y40zM7OvhFw9BfwDwX8IZAxs9fgQ8A+CjEJnCUkI0WGdCnRs28ZqqQL69ttvc8455wAwaNAgqqqqum50Qog+ryEc48qnltEQMRg/JIt7ZozteL2tXasTx98fzwNp2O6cH3Jih/Jk9puIDM5uqtQcKF/lPNecJqJgObuPoo3O0s+z34O6zXv/HHc6ZA+DwB5uaXlSo0WIvXj//feZMmUKtbW1ZGZmdvj9na6jc9dddzF16lTmz5/P7NmzASgtLSUvL68zHymE6Icsy+Ynz65kw44m8jO8/OmiY/Doncin6Uh14g7kybSrmSc49VpsG7BbPW4lEtyZR5OaC7ljIPcwGDDGueWMdGZ8JJgRosd1KtB56KGHuOiii3jxxRf52c9+xogRzjT0f/7zHyZNmtSlAxRC9C6WZbO6LEhNc5SAz2mkubdt3Y++t543v6jAranMueQYBqR3Mndk1+rEuxbt60R1YtMyEzk6bWZ0LMtZetq8CDYvdpbMmip3nm8d4yiqs+ykuZxickXfAF+gc1+jEKJbdGov5BFHHMGqVauor6/nF7/4ReL4b3/7W/72t7912eCEEL3LwvVVXPrEUq76+zL+79mVXPX3ZVz6xFIWrt99yfqj9VU88LYzY3LXeYe3exv5HrXOk6neADWlUPeVc1+9wTnewerEWxq2EDbDeDUvgxuq4IPfw9PfgvuHwOxJ8OpPYNVzO4MczQ2KjvNjUwE00H2QnucEWv6BEuSILmFZFvfffz8jRozA4/EwePBgfvOb3wCwatUqTj31VFJSUsjOzubKK6+ksbEx8d7LLruM8847j7vvvpu8vDwyMzO58847MQyDm2++mUAgwMCBA3niiScS79m0aROKovCvf/2LSZMm4fV6Ofzww3crADx//nyOO+44PB4PBQUF3HrrrRiGkTg/ZMiQ3XJ0jzzySH75y18mniuKwuOPP87555+Pz+dj5MiRvPTSS23e89prrzFq1ChSUlKYMmUKmzZtOqDv5wEVfYhGo2zdupXNmzezefNmKisr2b59+wENSAjROy1cX8XtL6yiZHuQVI9ObrqHVI9OyfYGbn9hVZtgp6oxwg3/XoFtw4XjB/Ht8YMO7OLxPJlgmTObo6pOYKOqzvNgmXO+vdWJa79izfI/AzAi1IT2+FR4505Y96ZTQM+VCsNOgVNug7MfgJQcZ2eUAuguJy9I18EMO3k5liUtIESXue2227j33nu54447+OKLL3jmmWfIy8ujqamJ008/naysLD7++GOee+453n77ba677ro273/33XcpKytjwYIFPPDAA/ziF7/gnHPOISsriyVLlnD11Vdz1VVXsXXr1jbvu/nmm/nJT37Cp59+ysSJE5k2bRrV1c4s6bZt2zjrrLMYP348K1euZPbs2cydO5e77rqrw1/fr371K7797W/z2WefcdZZZ3HRRRdRU1MDwJYtW5g+fTrTpk1jxYoV/OAHP+DWW2/t5HfS0eldVzNnzmThwoVtjtu2jaIomKaUyhaiP7Esm9nzN9AYMcjP8CaSib2qRn6GSnkwwuz5Gzh+mPPL/qZnV7KjIcLI3DR+Ma24aweTyHOxdz7fNWdmty/AhC1L4MtXYe08qF7Pmiw/ZPoZHQk5NWSGngRDJsPg4yFvbEsFZJzWCO/+uqUqcust5qrz2IyAbTjvEeIANTQ08Ic//IFHHnmESy+9FHAK9E6ePJnHHnuMcDjMU089RWqqU6n5kUceYdq0adx3332JHNlAIMDDDz+MqqqMHj2a+++/n+bmZm6//XZgZyD14YcfcuGFFyaufd111zFjxgwAZs+ezbx585g7dy633HILf/rTnxg0aBCPPPIIiqIwZswYysrK+OlPf8qsWbNQO7AJ4LLLLuM73/kO4LSUevjhh1m6dClnnHEGs2fPZvjw4fz+978HYPTo0axatYr77ruv09/TTgU6l19+Obqu88orr1BQUCAdy4Xo51aXBdlQ2UiWz73bv3dFUcj0udhQ2cjqsiALN1SxYO0OPLrKI989uqUB5wEqXwlNVc6Oq3D9Llu8U5yid01VbZORo02w4T1Y85oT3LTO31E01mbkAhFGHX0lTLh5Z2Czq4pVO9sC2AZYirPrSrGdACveIqBiVecKBgrRSklJCZFIhNNOO22P58aNG5cIcgBOOOEELMtizZo1iUCnuLi4TeCRl5fH4YcfnniuaRrZ2dlUVrbKPYM2Tbl1XefYY4+lpKQkce2JEye2+fd/wgkn0NjYyNatWxk8eHC7v8Yjjjgi8Tg1NZWMjIzEWEpKSpgwYcJex9UZnQp0VqxYwfLlyxkzZswBXVwI0TfUNEeJmTZuTcXGJhy1MCwLXVXxulU8mkq9ZbN0Uw2/fcPZsv3LbxQzOj+9awYQT0ZOy3OWiMK1YMScZSRvFmA77REayuHL15zcmjWvt+3f5PXDyNNhzFkw/FTWvPJNaCpn9PCv7z3IiV9b1ZylscZKpypynKI5xftsq0OJ0ELsTUpKygF/hsvVttWGoih7PBYvE9NVVFXF3mV2NRbbvSVJT4yltU4FOocddpjUyxHiIBLwuXFpCnWhGPWhGBHDxLadlRuPruFPcaEp8Of5GzAsm3OOKODCA83LaW1fRfuaasCd4szgPP8Dp65NXOZgGH22E9wMnpjotVQfqae8qRzYR+uH1te2LAi3BDKaa+eMjmU5/aW8WZKjI7rEyJEjSUlJ4Z133uEHP/hBm3NFRUU8+eSTNDU1JWZ1Pvroo8QS1YFavHgxJ510EgCGYbB8+fJE/k9RURHPP/98IkUlfu309HQGDhwIwIABA9rk6QaDQUpLSzs0hqKiot2SkxcvXtzprwk6mYx83333ccstt/D+++9TXV1NMBhscxNC9C/FhRlkp7nZXh8iFDVQFQVdU1AVhVDUYHt9iHDMorIhwuCAj7und6Io4L7sKRlZ0ZyZlFijE2zEmp0gJ70Ajr8WrngXfvwZnHmvk3/TqqFkvH7OIWmHkO7ez6xT3tiWJSuzJQFaBy3e4drlHJccHdFFvF4vP/3pT7nlllt46qmn2LBhA4sXL2bu3LlcdNFFeL1eLr30Uj7//HPee+89rr/+ei655JIuqWH36KOP8sILL/Dll19y7bXXUltby/e//30AfvjDH7Jlyxauv/56vvzyS/73v//xi1/8gptuuimxTHbqqafy97//nQ8++IBVq1Zx6aWXonWwD93VV1/NunXruPnmm1mzZg3PPPMMTz755AF9XZ2a0Zk6dSrAbmuIkowsxEFgl1xgFCfeqAvFUBV48IIjyfB2U5dq23auv6eifSnZ8M0nYOhkZ6lpHxKFAvc3mwO75+jYGs4gbLBNydERXe6OO+5A13VmzZpFWVkZBQUFXH311fh8Pt544w1+/OMfM378eHw+HzNmzOCBBx7okuvee++93HvvvaxYsYIRI0bw0ksvkZPjNN495JBDeO2117j55psZN24cgUCAmTNn8vOf/zzx/ttuu43S0lLOOecc/H4/v/71rzs8ozN48GCef/55brzxRv74xz9y3HHHcffddycCrs5Q7F0X1Nph1731uzr55I4310umYDCI3++nvr6ejIyMZA9HiF5n1dZ6rvr7MhSFlqUrK7F05dYUQjELy4YZRx/C7799ZNcPYO0b8J+ZEGvapQGm4uyEilcdvvAf7Qo2Zn00ixfWv8BVR1zFdUddt+8Xr38bXvyhs+W8acfuva5SBzjjOu9PMGLqgXyVQiTFpk2bGDp0KJ9++mmij2V/0qkZnb4WyAghDkw8GTk33UOmz0V9s0HMtHBpKg2RGFbUQlMVzji8oGsvvO0TWPwn+Py/zuwJsHNKiZ0zLS6Pk6PTzoTgdjXzjIvnB2kupy+VEXISklXd2fFlhMHsWLFCIUTP6VSgA1BXV8fcuXMTW8+Ki4v5/ve/j9/v77LBCSF6h70lI9vYmC0TLFkpLvIzvAd+Mct0dkwtehQ2t6rVpbpagh2lZZeUClg7i/a1MyHYsAzW164H9tPMMy5/nNOrqmK1k//j8u08Z9tOgnRecfuLFQohelSnkpGXLVvG8OHDefDBB6mpqaGmpoYHHniA4cOH88knn3T1GIUQSdY6GTkcjTJWKWWSuhrF2pmPV5DppbjwAJZ+LRM+exYenQD/vsgJclQdjrgAfvAueNJa6ta4nB5TCi33eocSgkvrS4laUVJdqQxMH7j/cakqTL7RuX7DdoiFWpKgQ85zT7pzvgMF04ToTYYMGYJt2/1y2Qo6OaNz44038o1vfIPHHnsMXXc+wjAMfvCDH3DDDTewYMGCLh2kEKJ3mKis5mr9JYYrZcyKXYqBxiClkoFU0cjxnftQ04DPn4cF90O1M9OC1w/Hfh+OuxIyCqHs0y4r2vdlzZeAM5ujKu0MToadDOc8BB8+6HRJD9c5AVdesRPkDJPlfCF6q04FOsuWLWsT5IBTRfGWW27h2GOP7bLBCSF6h9VlQYYEl3OL56+kWM3MM4/lbetYVEzu1h9niF7N/UEXq8uOYOzAdi5fmwZ8/h+Yfz/UbHCOpWTBxOucAMfbanaoC4v2ralx8nPGBDpY8HTYyTDkRKf6cnO1s0yWP05mcoTo5ToV6GRkZLB58+bdKiNv2bKF9PQuqoQqhOg1aprCXBT7L6lKiG1qAY9EzwfgfH0x+a4QqXaIi4z/UtP0HWA/gY5tO7uo3v4F7HBmV0gJwKSWAMezh58hXVi078ta55odDnTACWpkC7kQfUqnAp0LLriAmTNn8rvf/Y5JkyYBToXEm2++OdGoSwjRfxSG1pLONoKk87RxCuV2gBylnktc74GiELTTGWJvoyG0FthH4bJtn8Cbd8BXHzrPU7Jg0o/guCv2HODEtS7ap3lAbV2M0G53Y03btncuXbVnx5UQos/r1Jzr7373O6ZPn873vvc9hgwZwpAhQ7jsssv45je/eUAdRndlmiZ33HEHQ4cOJSUlheHDh/PrX/+6TS8N27aZNWsWBQUFpKSkMHXqVNatW9dlYxBCwHBfBK9qUmIW8FzsBAB+rD1PmllH1DBptjW8qsVwX2TPH1C/zamD89gUJ8jRPHDCDfCjFXDiTfsOcmAPRfssZ2bItpznrXN09vUxzRXUR+rRFZ0RmSM6/o0QQvQ5nZrRcbvd/OEPf+Cee+5hwwZnbX348OH4fL79vLNj7rvvPmbPns3f/vY3iouLWbZsGZdffjl+v58f/ehHANx///08/PDD/O1vf2Po0KHccccdnH766XzxxRd4vV2w1VUIgZqWjap7mBM6BxON09TlXKi9hw1EcFNv+dA8btS0XZaOjCgsfhTm/9YpqocC4y6EKT+DzA70worn6PgH7aFoX8rOon37ydGJz+YMyxyGW3N35FsghOijOl1HB8Dn85GZmZl43NUWLlzIueeey9lnnw04W+D++c9/snTpUsCZzXnooYf4+c9/zrnnngvAU089RV5eHi+++CIXXnhhl49JiIORlXcEf4+eyEp7BCmEuV1/GgMNBRsvEVLUCJut4aTmHbFzmnj9O/D6LTt3Ug2aAGfeD4VHdnwAXVS0Lx7odCo/RwjRJ3Vq6cowDO644w78fn9i6crv9/Pzn/98jy3ZO2vSpEm88847rF3r9KVZuXIlH374IWeeeSYApaWllJeXJ3pvAfj9fiZMmMCiRYv2+rmRSEQakQrRAYs2VvPn8NcB+LH+AgPVahSFlpuCAkRiJqvLgtC4A567HP4x3QlyUnPhvDnw/Tc6F+TAzqJ9oVrnucsHnoydxftCtc75/RTta721XIiDkWXZrNpaz/y1O1i1tR7L6nAXqE5ZtGgRmqYlJi56UqdmdK6//nr++9//cv/99zNx4kTA+SJ++ctfUl1dzezZs7tkcLfeeivBYJAxY8agaRqmafKb3/yGiy66CIDy8nKA3bq25uXlJc7tyT333MOvfvWrLhmjEAeDv7/7KUFSGaGUcaH2Hio2ChY2ChHcNCppZNr1hJc+AOv/7gQeigoTroZTbnXq4hyIeNG+V25wivSlZIHuASPiXKudRftkRkcczBaur2L2/A1sqGwkZtq4NIXhuWlcc/JwJo3I6dZrz507l+uvv565c+dSVlZGYWFht16vtU4FOs888wz/+te/EjMrAEcccQSDBg3iO9/5TpcFOs8++yxPP/00zzzzDMXFxaxYsYIbbriBwsJCLr300k5/7m233cZNN92UeB4MBhk0qAP5AkIcRNaUN/DmphigcIX+GtsYQIbVgAuTGBpBJR0vMQZRjnvVw86b8sbCuX/s2q3YB1i0LxgNsq1xGyA7rsTBZ+H6Km5/YRWNEYMsnxu3phI1LUq2N3D7C6u4+/yx3RbsNDY28u9//5tly5ZRXl7Ok08+ye23394t19qTTgU6Ho+HIUOG7HZ86NChuN1dl+B38803c+uttyZybcaOHctXX33FPffcw6WXXkp+fj4AFRUVFBTsbCZYUVGxz1LWHo8Hj8fTZeMUor+ybZtfv/IFlq0wRVvJJGUlmTTjUaIoOL2ucqjHY8dQFbA1F8optzlbxjVX1w/oAIr2ra1xlsALUwvxe6Qnnzh4WJbN7PkbaIwY5Gd4URSnPINX1cjPUCkPRpg9fwPHD8tGbVO6oWs8++yzjBkzhtGjR3PxxRdzww03cNtttyXG0d06laNz3XXX8etf/5pIZOdW0kgkwm9+8xuuu+66Lhtcc3Mz6i4/wDRNw7KcLoJDhw4lPz+fd955J3E+GAyyZMmSxJKaEKLz3i6p5MP1Vbg1hcvd71Kg1OAlgoWKgYqGTYriBDkhvNhXfggn/qR7gpy4eNG+EVOd+3ZWJpb6OeJgtbosyIbKRrJ87t2CC0VRyPS52FDZ6OTYdYO5c+dy8cUXA3DGGWdQX1/P/Pnzu+Vae9KpGZ1PP/2Ud955h4EDBzJunJP8t3LlSqLRKKeddhrTp09PvPa///1vpwc3bdo0fvOb3zB48GCKi4v59NNPeeCBB/j+978POP8D3XDDDdx1112MHDkysb28sLCQ8847r9PXFUJAxDD5zatfAPCNcYUM+KIaRVGwbQUFEw8WiuKUs7EUhU0UYETz2H9bzeSQ/BxxsKppjhIzbdzanv8o8Ggq9ZZNTXO0y6+9Zs0ali5dygsvvAA47aIuuOAC5s6dyymnnNLl19uTTgU6mZmZzJgxo82x7shx+eMf/8gdd9zBD3/4QyorKyksLOSqq65i1qxZidfccsstNDU1ceWVV1JXV8fkyZOZN2+e1NAR4gA9+dEmNlU3MyDdw3kFNWR+EaSKAAGlFg/OrKoNhBUPQSWDTLuBsm0rYNCUpI57b9bUdrLHlRB9XMDnxqUpRE0Lr6rtdj5iWrhUhYCv62tLzZ07F8Mw2iQf27aNx+PhkUcewe/v/mVkxW5dZvggFQwG8fv91NfXk5GRsf83CNHP7WiIMOV379MYMfjtN4/gWONT/PN+SBoh3DgNNS3ARCWGmxolE68dpuGMPzL0+G8kd/B7EDNjHPfMcRiWwRsz3qAwred2fAiRbJZlc+kTSynZ3kB+hqfN8pVt25QHIxQVpPO3y4/r0hwdwzAYOHAgt9xyC1//+tfbnDvvvPP4v//7P66++uouu97edCpHJxQK0dzcnHj+1Vdf8dBDD/Hmm2922cCEEMnzwFtraYwYHDHQz4yjB3JodD1ZNODGwAZiaERxY6HhJkauvQOXCocO7J27FzfUb8CwDNLd6RSkFuz/DUL0I6qqcM3Jw0nzaJQHI4RiJpZlE4qZlAcjpHk0rjl5eJcnIr/yyivU1tYyc+ZMDj/88Da3GTNmMHfu3C693t50KtA599xzeeqppwCoq6vjuOOO4/e//z3nnntul20tF0Ikx/rKRp5dtgWAWWeNRJ13C+q7v2rZZQVhdEx0QMFGJYaKhoVPB7Wgd2botM7P6amdHkL0JpNG5HD3+WMpKkinOWJQ2RihOWJQVJDebVvL586dy9SpU/e4PDVjxgyWLVvGZ5991uXX3VWncnQ++eQTHnzwQQD+85//kJ+fz6effsrzzz/PrFmzuOaaa7p0kEKInvO7N9ZgWjbTR7k49v3LYPNC54Segm3EcGFjYcY7TTnFAxUdl647TTW7snZOF1lTI/k5QkwakcPxw7JZXRakpjlKwOemuDCjW7aUA7z88st7PXfcccfRU5kznQp0mpubSU93ug2/+eabTJ8+HVVVOf744/nqq6+6dIBCiJ7zyeZa5q0u50h1PfdV/wmayp1WCxN/CMueQFVUlMZKNNvY+SZFQ0nLcTqJ76epZrLIjishHKqqMHbgwVVHqlNLVyNGjODFF19ky5YtvPHGG4kko8rKSknmFaKPsm2be1//km9q83nOcxeupnLIGQVXvAujzgDLgqYdKICiuVBUt3MPTkdxy9pvU81ksG07MaMjPa6EOPh0KtCZNWsW//d//8eQIUM47rjjEsX53nzzTY46qvdNWwsh9u/9LyuYsuVRfuf6My47CqPPgh+84zTLzBsLtgGW6bRdUHXQNOdedTnHbcN5XS+zrXEbDbEGXKqLYZnDkj0cIUQP69TS1Te/+U0mT57M9u3bEwUDAU477TTOP//8LhucEKJnmJEm9P9ezjV6Sz7OSbfAKbftrDpcsQqUlsDGNsDWIJ6ebJvOcUXrlTk68WWrEZkjcKndWLFZCNErdWpGByA/P5/09HTeeustQqEQAOPHj2fMGFkDF6JPaaigfvbXOTG2kCg6zWfPhlN/1ra1QnM1qBr4B4Ge4ixTWYZzr6c4x1WtV+borK5eDcBh2YcleSRCiGToVKBTXV3NaaedxqhRozjrrLPYvn07ADNnzuQnP/lJlw5QCNGNKkuwHzuVQN3n1NhpvHbUHHzjv7v763zZzhKV5oLAMAgMhczBzn1gmHNcdfXKHJ3VVU6gU5xTnOSRCCGSoVOBzo033ojL5WLz5s34fL7E8QsuuIB58+Z12eCEEN1o82L46+kowa1ssAq4wn0vZ5w1fc+vzR/n5OqEap3nLp+zG8vV8u8/VOuczx+35/cniW3biRmd4mwJdIQ4GHUq0HnzzTe57777GDhwYJvjI0eOlO3lQvQFa16Hp86FcD2fKaOZEf0l5516Il7X7n1wAGcZa/KN4EmDhu0QCznbyWMh57kn3Tnfzk7iPWVr41aC0SAu1cXIzJHJHo4QIgk69VOpqampzUxOXE1NDR6P54AHJYToRp/+A/51ERhhNuecyLdDt5Kamcu3x++nfcOwk+GchyCvGKJN0Fjh3OcVwzkPOud7mfhszqisUbg0SUQW4mDUqUDnxBNPTLSAAFAUBcuyuP/++5kypXd2LhaiX7EsKPsU1r/t3FtW+9734UPwv2vBNjGO+A7fqr2WMB6uP3UEHn0vszmtDTsZLv4vXPgPOO9Pzv3F/+2VQQ7AF9VfALJsJUSyXHbZZSiKkrhlZ2dzxhln9Ejrh7hObS//7W9/y6mnnsqyZcuIRqPccsstrF69mpqaGj766KOuHqMQorWN8+HDB6FqHVgxJwk4Z6SzdLS3gMO24b27YcH9zvMTbuBx1yVUNK1hcMDHjGMG7vl9e6KqvW4L+d58UdUS6EgishAOy4Lylc4OSV+2k1fXzUvOZ5xxBk888QQA5eXl/PznP+ecc85h8+bN3XrduA4HOrFYjB/96Ee8/PLLvPXWW6Snp9PY2Mj06dO59tprKSiQzsBCdJuN8+GVGyDSCClZoHvAiEDFauf4OQ/tHuzYNrx1Byz8o/N86i9pHH89f77vXQB+dNpIXFrvyq3pCpZtyYyOEK115o+kLuDxeMjPzwec0jS33norJ554Ijt27GDAgAHddt24Dgc6LpeLzz77jKysLH72s591x5iEEHtiWc4PqUgjpBdAvAu3KwV0r5MU/OGDMOTEnX+hWRa8fjN8/Ljz/Mz7YcJVPPnuOmqbYwzLSeW8IwuT8/V0sy0NW2iINeDRPFIRWYjO/JHUDRobG/nHP/7BiBEjyM7umXIUnfoz7uKLL2bu3LldPRYhxL6Ur3T+EkvJ2hnkxCmKc7xqnfM6cNoyvHx9S5CjwLQ/wISrqA/F+MuCjQD8eOpI9H44mwM76+eMDoyWisji4LbrH0muFFBU5z69wDn+4YPtz/XroFdeeYW0tDTS0tJIT0/npZde4t///jdqD+3S7FSOjmEY/PWvf+Xtt9/mmGOOITU1tc35Bx54oEsGJ4RopbnamW7WW3Y2xppbek9pTj0b3QPhupbXWfDSj2DFP5wfaOfNgXEXAPDkR5sIhg1G5qZxzhH9czYHWlVEDkhFZHGQ68gfSd2QfzdlyhRmz54NQG1tLX/6058488wzWbp0KYceemiXX29XnQp0Pv/8c44++mgA1q5d2+acsus3UQjRNeLViUO1EKoDIwzYgOIsXaVkOudTAvDqjS1BjgYzHofDnUKAjRGDv35UCji5OZraf/+9JgoFSiKyONjt+kfSrlr/kdQNUlNTGTFiROL5448/jt/v57HHHuOuu+7qlmu21qlA57333uvqcQgh9id/HKTmQPkq57mm46w+W87sTqwZ8g6HFc/A8iedmZzz/5wIcgD+sfgr6kMxhg1I5ayx/XfjgGVblFSXAJKILETijyQj4ixX7cqI9GgLF0VRUFU10Sezu3Uq0BFCJFli5tTe+dyyoKEMPn4MUODcR+GIbyXeEoqaPP6Bk5tz7Skj+vVszqb6TTQbzaToKQz1D032cIRIrngLl4rVzuxv65UX23ZmifOKu62FSyQSoby8HHCWrh555BEaGxuZNm1at1xvVxLoCNFXlK+EpioneTBc33bpSvOCrkDTDue10/4AR7ZtzvnPpZupaowyKJDCN/rpTqu4VVXOrFdRoAhdlR9z4iAXb+Hyyg3O7szWu65Ctd3ewmXevHmJ0jPp6emMGTOG5557jlNOOaVbrrcr+QkgRF8RX2dPy3OmmMO1YMRAd0EsAk0VzuuOuwKOubTNWyOGyZ8XbADgmpNH9Mu6Oa2t3OHsPDtiwBFJHokQvUS8hUu8jk64zlmuyivu1jo6Tz75JE8++WS3fHZ7SaAjRF+xt2Rk23JuAB4/HHnRbm/9z/KtVAQjFPi9zDjmkB4ddjJ8tsMpLz9uQO/qpi5EUg072amz1cOVkZNNAh0h+oo9JSNbrYIcFAgM3W2dPWZazH7fmc256qRh7etp1Yc1x5pZV7cOkBkdIXbTh1q4dBUJdIToixTFqaFjGfEDLbfdvfjpNrbWhshJc3PhcYN7bIjJsrp6NZZtkZ+aT64vN9nDEUIkWf+erxKiP2mdjKy62gY5rlTIKHTOxysjA6Zl86eW2ZwrThyG19W/Z3OgVX5OjszmCCFkRkeIviOejOzJADPiHNM8kJbrFAnEhsaKNkW/XvmsjNKqJjJ9Li4+vvsrkPYG8fwcWbYSQoAEOkL0Hb5sQIGaDS15OYozq9NQ7iQnxysjtxT9siybR99bD8DME4aS6un//9xt25ZEZCFEG7J0JURfkTEQQjU7k481F2huJ7kw1gzBMidZuSUZ+a2SCtZWNJLu0fnepCHJG3cPKmsqozpcja7qjAmMSfZwhBC9gAQ6QvQFsRD86ztgRp3nitZS3bRVZeRWbNtmznwnN+eSiYfiTzk4unfHZ3PGZI3Bq3uTPBohRG/Q/+eyhejrLAteuBq2fgwokJrrzOAYIad8u6KAngJefyIZeVn0UD7dXIdbV7nshCHJ/gp6jOTnCCF2JYGOEL3d27+AL14EVQd3OnhSIdrUMpljg604AY/udo43V/Pnj5wdWTOOHkhu+sEzsyGBjhBiV7J0JURv9vHjsPBh5/HJtzjBTt1mbDOMpWqYqhtL1bDNMNRtBstibSiDt0sqURS44sT+3dDSsi1WV6/mo20fsaJyBSU1TsdySUQWoncpLy/n+uuvZ9iwYXg8HgYNGsS0adN45513uv3aMqMjRG+19k147Wbn8ZSfweSbYPFsbMskigsnJdnJ0VFRcRNDsQ3+UuLk45x+WD7DBqQlZeg9Ycn2JcxdNZfSYCmGZWDZFjErRro7nUPS+n+bCyE6w7ItSmpKqAvXkenNpChQhKp075zHpk2bOOGEE8jMzOS3v/0tY8eOJRaL8cYbb3Dttdfy5Zdfduv1JdARojfavhKeu8zZYXXkxXDSzbB9BVFLRUFFwwRUbBQUbDQsDFS2GwH+t6IMgKtOHpbUL6E7Ldm+hDsX3UlTrAm/x49bc1PZVAmAYRksLV/KhIIJSR6lEL3Lrn8c6KrO0IyhzBw7s1v/vfzwhz9EURSWLl1Kampq4nhxcTHf//73u+26cbJ0JURvE9wOz1wAsSYYdgpMewgUBauxmuaYRRk5RHCjYqNjomITwU0ZA3g8dDIxCyYMDXDU4KxkfyXdwrIt5q6aS1OsiVxfLl7di6qoRC1nR5qCwtxVc7ESPcCEEPE/DtbWrsWn+8hJycGn+1hbu5Y7F93Jku1LuuW6NTU1zJs3j2uvvbZNkBOXmZnZLddtTQIdIXqTWAj+9V1o2A45o+HbTzn1coANzR7CloatuClTC9mu5lOh5rJdzadMLSRIOv8xJgNw9cnDk/lVdKuSmhJKg6X4PX4URSFshGmINNAUawIgy5tFabA0ka8jxMFub38ceHUvub5cmmJN3fbHwfr167FtmzFjklfXSpauhOgtbBv+dx2UfQIpWfDdfzlbxluUpYzCwyGM4iuqySaEFydHR0HF5k3jKJrxMjAzhVNGD0jal9Hd6sJ1GJZBzIpREawgYkawLRu7JV9JQcGwDOrCdckdqBC9xK5/HLSmKAp+jz/xx0FxdnGXXtu27S79vM7o9TM627Zt4+KLLyY7O5uUlBTGjh3LsmXLEudt22bWrFkUFBSQkpLC1KlTWbduXRJHLEQnffB7+Pw/zs6qb/8dAm1zbAKpXp52TaeJFDLNKlQjRMwwUI0QKUYdTxunAvCtYwft9sOsP8n0ZmLZFtsbtxM2wqiK2ubrLW8qx7ItMr2ZyRukEL1I/I8Dt+be43m35u62Pw5GjhyJoijdnnC8L7060KmtreWEE07A5XLx+uuv88UXX/D73/+erKyduQf3338/Dz/8MHPmzGHJkiWkpqZy+umnEw6HkzhyITqo5GV499fO47N+B0NP3O0lxYUZbMo4hlsj3+dLazCpSphcpZ5UJcxfzTOpIQOXpnB1P05CBhidNRrTMjFtE13RUVETU+6aomHaJqZlMjprdJJHKkTvkOnNRFd1ovHK6ruImlF0Ve+WPw4CgQCnn346jz76KE1NTbudr6ur6/Jr7qpXL13dd999DBo0iCeeeCJxbOjQnXVBbNvmoYce4uc//znnnnsuAE899RR5eXm8+OKLXHjhhT0+ZiE6rHwV/Pcq5/FxV8Gxl+/z5YvsYpaYRRxhbyagNFBtpfOp6XQmz0lz49J69d8vB2xN7RpUVUVTNQzbABus+GZ720ZTNVRVZU3tmi6fhheiLyoKFDE0Yyhra9fi0TxtZkBt26Y+Us+orFEUBYq65fqPPvooJ5xwAscddxx33nknRxxxBIZh8NZbbzF79mxKSro3n65X/0R86aWXOPbYY/nWt75Fbm4uRx11FI899ljifGlpKeXl5UydOjVxzO/3M2HCBBYtWrTXz41EIgSDwTY3IZKicQf88zs7d1idfvdeX7q6LEh1Y5QCvxevy8Uqeyjvm0fwmT0EUFBaCiSvLuvf/z3XhevQFI0sTxY2NqZtJs7Z2GR5stAUTXJ0hGihKiozx84k1ZVKZXMlYSOMZVuEjTCVzZWkulKZOXZmt9XTGTZsGJ988glTpkzhJz/5CYcffjhf+9rXeOedd5g9e3a3XLO1Xh3obNy4kdmzZzNy5EjeeOMNrrnmGn70ox/xt7/9DXAqLQLk5eW1eV9eXl7i3J7cc889+P3+xG3QoEHd90UIsTdGFJ69BOq3QGA4fOtJ0PY+yVrTHCVm2mSmuDk020deupeAz4XW8tdZdqob03Ze15/Fc3Rqw7Vgk/jhrLT8X224VnJ0hNjFhIIJzJo4i1FZo2g2mqkKVdFsNDMqaxSzJs7q9rpTBQUFPPLII2zatIlIJMLWrVv53//+xymnnNKt14VevnRlWRbHHnssd9/t/JV71FFH8fnnnzNnzhwuvfTSTn/ubbfdxk033ZR4HgwGJdgRPe+N22HzIvBkwHf/7ey02oeAz41LU6gLxagPxYgYJpZlY7ZsatBVFVVxXteftc7RcatuYnYMcAIeXdGJWlHJ0RFiDyYUTGB8/vger4ycbL36qysoKOCwww5rc6yoqIjNmzcDkJ+fD0BFRUWb11RUVCTO7YnH4yEjI6PNTYge9enT8HHLMuz0xyBn5H7fUlyYQXaam+31IUJRA1VRiG/cVIDKhjDZaW6KC/v3f8+tc3RiViyRiKygYNhGmxwdIURbqqJSnF3MCYecQHF2cb8PcqCXBzonnHACa9a0/WG1du1aDj3USbwcOnQo+fn5bZqCBYNBlixZwsSJE3t0rEK027ZP4JUbncen3A6jz+j4Zyhg2TZWS6TTj3eT7yaeo1OQWoBH9ySO29h4dS8FqQWSoyOESOjVgc6NN97I4sWLufvuu1m/fj3PPPMMf/nLX7j22msBp9DRDTfcwF133cVLL73EqlWr+N73vkdhYSHnnXdecgcvxJ407oB/XwJmBEaf5fSwaqfWycgpLg3DihfIA59bI9/vpbox2u+TkeNbZV2qC7/bKajoUl3kp+YzOH0wLtXVbVtlhRB9T6/O0Rk/fjwvvPACt912G3feeSdDhw7loYce4qKLLkq85pZbbqGpqYkrr7ySuro6Jk+ezLx58/B6vUkcuRB7YMacRp3BrZA9As6fA2r7/9aIJyPnpntI9eisrWgE4JDMFDJ9LmwbKhsj/T4ZOb5VdnX16kTbB9MyqWiuoDZci67qFGcXd9tWWSFE39KrAx2Ac845h3POOWev5xVF4c477+TOO+/swVEJ0QlvzYKvPgR3Glz4TJv2Du0RT0aOmlYimPHqKl6XBgpEDAuXqvT7ZGRVUZlUOIml25di4mwt1xQNG5uQEUJTNCYVTjoocg+EEPsnPwmE6Akr/w2L/+Q8Pn8ODOj4jqDiwgyG56ZREQxR3egEOlHT4quaJkp3NLGjIczw3LR+n4xs2RYLyxbi1bxtjgGk6Cn4dB8LyxZK93IhBCCBjhCdZ1lQ9imsf9u5t/byi3X7Snj5R87jk26GommdupyqKpw0MofGyM4CeXrLv+DmqElTxOSkkTmoav/OTI43KHTrzsyVV/MSSAmQ58vj0IxDyfHlSPdyIURCr1+6EqJX2jgfPnwQqtaBFQPV5WwRn3wjDDt55+uaa+BfF4MRhhFfg1Nu6/QlLcvm/TU7iDcD1hSwbKciss+toakqC9ZVMXPysH4d7MQbFDYbzQBErahTPFCBukgdAW9AupcLIRJkRkeIjto4H165ASpWgzsV0vKc+4rVzvGN853XWRa8cBXUb4asoTDjMVC1Tl92dVmQz8vqsQGXpjAkO5WBWSkcGkhlaE4quRkeNlQ2HhS7rkzLJGJGACc/R1d1VEUlbITZ3rhdKiMLIRIk0BGiIyzLmcmJNEJ6AbhSQFGd+/QC5/iHD7a87gFY9yboXvj2U/utfLw/VY0RGsPOslV2mrtNYz4U8GgqMcvu97uuRmeNTgQ5AGrLjzEVpzKydC8XQrQmS1dCdET5Sme5KiVr9yp9iuIcr1oHy5+A937jHD/rd1BwxAFfek15A6ZtoyjQEDKoboxi285lPbqGP8V1UOy6WlO7JtHIU0HBxsb5f6fBp3QvF0K0JjM6QnREc7WTkxOvyBtrhkiDcw/OcSMM7/wKbAuOvBiOvuSAL2vbNq99vr3lMYRjJqqioGsKqqIQihpsrw8dFC0g4jk6AG7NjYWFaZtYWFIZWQixG5nREaIjfNlO4nGoFkJ1TlCDDSjOEpXXD5EgmFHIOxzO+m2XXHZJaQ2fba0nMYcUf9C62ZW929v6JcM2EjM6g9MHO88tZybHq3mJmBFiakxydIQQgAQ6QnRM/jhIzYHyVc5zTceZGLWcWZ1oI2CDO93Jy3H7uuSyf56/AYAUt0aWz9XSvdxKLF2luDQyUlyJFhBjB3asGGFfsrVhKwC6ojutIBRX4pxt29RH6hmVNUoqIwshAAl0hOg8ZddpFXvn43MfgezhXXKZL8uDvLdmB4oCXl0jM8VNls9NOGZhWBa6quJ1qQdNC4j5W5xdbWnuNCqbK/F7/Lg1N1EzSn2knlRXKjPHzpTKyEIIQHJ0hOiY8pXQVOXssNJTnN1VlgGm6eTkgNPiIevQLrvkX+ZvBGDS8By8LpWoaaEoCilujXSvixS3hqIoRMz+3wKiIdrAxxUfA/B/x/4fo7JG0Ww0UxWqotloZlTWKGZNnMWEgglJHqkQoreQGR0hOiKejJyW5+TrhGshFoVwNZiAy+fcmqu75HLb6kK8tLIMgP/7+igeeGstJdsbyM9Q22wvt22buuYYRQXp/ToZ+aOyjzAsgyEZQzh3xLlMGz6NkpoS6sJ1ZHozKQoUyUyOEKINCXSE6Ig9JSNbBomEZK8fFM15XRd4/IONGJbNpOHZHDU4i2tOHs7tL6yiPBgh0+fCo6lETIu65hhpHo1rTh7er6siv7f5PQCmDJoCOA0+ZQu5EGJf5E8fIToinowcLGvZUt4qLwcFGiud8/njDvhStU1R/rV0CwBXn+zk+0wakcPd54+lqCCd5ohBZWOE5ohBUUE6d58/lkkjcg74ur1VzIrxwbYPAJgyeEqSRyOE6CtkRkeIA9FSz8Vp7aCQaETVBf62aBOhmElxYQYnjtwZwEwakcPxw7JZXRakpjlKwOfUzunPMzkAn1Z8SkO0gSxPFkfkHHgBRiHEwUECHSE6Ip6MnJYHDdt3Hrdx2kB4/c758pVQeFSnL9McNfjbwk0AXHXy8LbtHnA6mffnLeR78t4WZ9nqpIEnoR1AzzAhxMFFlq6E6IjmaqcYYLiettX6WmZzdLeTrHyAycjPfryF2uYYgwM+zjo8/0BH3efZtp0IdOL5OUII0R4S6AjREb5siDZBrMl5rrmcisiaBmYY6jY7W84PIBk5Zlo89kEpAFecNAxdk3+mn1d9zrbGbaToKUwsnJjs4Qgh+hD5CSpER8TCO4McRQdVb5nQUZ3nlgm2AXljO32JVz/bzra6EDlpbr51zMCuGXcf91rpawCcMugUfK6uqTYthDg4SKAjRHs1VcGz8QadLc2l4sUCLcMJcFTd2V5esapTl7Btmzkt7R4uP2EoXpfkopiWybxN8wA4e+jZSR6NEKKvkUBHiPawTHj+B9C0wwlk0nKd42YMrKhzb9vO1nJV63SOzvtrdvBleQOpbo2LJ3RddeW+7OOKj6kKVeH3+JlUOCnZwxFC9DES6AjRHvPvh43vOfk4noydgYzmAtXt3IMTCB1Ajs7sltmc704YjN/n2s+rDw6vbXSWrb5+6NdxafI9EUJ0jAQ6QuzP+rdh/n3O47MfcFatLNOpkKzqTiKyqjvPDyBH55PNtSwtrcGlKcycPKxrv4Y+KmJGePurtwE4a+hZSR6NEKIvkkBHiH2p3wrPXwHYcMzlkHeYs3Sl6k5AY1vOkpVtHXCOzpz3ndmc8486hHy/t4u/kL7pw60f0hBrIM+Xx9F5Ryd7OEKIPkgCHSH2xojCc5dBqAYKxsEZ9zpLVqoG/kFtu5dblvPcP6hTOTrrKxt484sKFAWuPGl493w9fdCrpa8CcObQM6VZpxCiU6QyshB789YdsPVjp9rxt58Cl3dnU0/NBYFhYIScQEfVnUDHCIPp6nCOzp/nbwTga0V5jMhN646vps+pCdfw/pb3AVm2EkJ0nvyJJMSefP48LJnjPD5vDmQNcR7nj4OckU73cgCXz0lOjtd2CdU65zvQ1HN7fYgXV2wD4OpTZDYn7n/r/0fMilGcXUxRdlGyhyOE6KMk0BFiVzvWwks/ch5PvhHGtJpNUFXnmCfN6XUVCzn5ObGQ89yT7pxX2/9P6y8LNhIzbSYMDXD04Kwu/mL6Jsu2eG7tcwB8e/S3kzwaIURfJoGOEK1FGp2igNFGGHIiTPn57q8ZdjKc8xDkFTvtIBornPu8YjjnQed8O+1oiPDPpZsBuP7UkV30RfR9i7cvZkvDFtJcaZwx5IxkD0cI0YdJjo4QcbYNr9wAO76EtHz45l9B28s/kWEnO4FQ+Uon8diX7SxXdWAmB+DxDzcSjlkcOSiTE0Z0vj9Wf/PcGmc2Z9rwadLyQQhxQCTQESLu48dh1XPO9vBvPbGz+vHeqCoUHtXpy9U2RfnHoq8AuP7UESiK0unP6k8qmioSncq/PUqWrYQQB0aWroQA2Loc5t3mPP7ar+DQ7m818MTCTTRFTQ4ryODUMfsJqg4i/13/X0zb5OjcoxmRNSLZwxFC9HES6AjRVA3Pfg+sGBRNg4nXdfslg+EYT3xUCshsTmsRM5JYtvrW6G8leTRCiP5AAh1xcLNM+O8VENwKgeFw7qPQA0HH3xd9RUPYYGRuGqcX53f79fqK/63/HztCO8jz5XH6oacnezhCiH5AAh1xcFvwW9jwjlPs74K/O8UBu1lz1ODxD5wCgdedOgJVldkcAMMy+OvnfwXg8sMvlwaeQoguIYGOOHitfxvev9d5fM6DzvbwHvD04s3UNscYku3j7LEFPXLNvmDepnlsa9xGwBtg+sjpyR6OEKKfkEBHHJzqtrRq1nkZHPmdHrlsOGbyl5bZnB9OGYGuyT9BcAoEPv7Z4wBcctglpOgpSR6REKK/kJ+y4uBjROC5S1s167yvxy797LIt7GiIcEhmCucfdUiPXbe3e2/Le2yo30CaK40LRl+Q7OEIIfoRCXTEweeNn8G25eDN3NmsswdEDYs5728AnJ5WLpnNAZzZnD+v/DMA3xnzHdLd6UkekRCiP+lTP2nvvfdeFEXhhhtuSBwLh8Nce+21ZGdnk5aWxowZM6ioqEjeIEXvtuKf8PFjzuPpf9nZrLMH/HvZFsrqw+RlePjWMQN77Lq93bzSeZTUlJDqSuWSwy5J9nCEEP1Mnwl0Pv74Y/785z9zxBFHtDl+44038vLLL/Pcc88xf/58ysrKmD5dEhnFHpR96rR4ADj5pzCq57Yvh2Mmj767HoBrp4zA69J67Nq9WcyM8fCnDwPw/cO/T5ZXmpoKIbpWnwh0Ghsbueiii3jsscfIytr5g7C+vp65c+fywAMPcOqpp3LMMcfwxBNPsHDhQhYvXpzEEYseYVlO8LL+befesvb+2qYq+PclYIRh1Blw8q09N07gn0s3Ux4MU+j3csH4QT167d7sHyX/YFvjNnJScri46OJkD0cI0Q/1iUDn2muv5eyzz2bq1Kltji9fvpxYLNbm+JgxYxg8eDCLFi3a6+dFIhGCwWCbm+hjNs6Hf0yHf10ML/7Quf/HdOf4rkwDnrsM6rc4RQHP/3OHm28eiFDU5NH3nNyc604diUeX2RyAHc07mLNyDgA3HH2DNO8UQnSLXh/o/Otf/+KTTz7hnnvu2e1ceXk5brebzMzMNsfz8vIoLy/f62fec889+P3+xG3QIPkLu0/ZON9ZgqpYDe5USMtz7itWO8d3DXbe/gVs+gDcaXDhM5CS2aPD/fviTVQ1RhgUSOFbx0puTtxvl/2WZqOZI3KOYNrwackejhCin+rVgc6WLVv48Y9/zNNPP43X23U7Y2677Tbq6+sTty1btnTZZ4tuZlnw4YMQaYT0AnClgKI69+kFzvEPH9y5jPXZc7DoEefxebMhd0yPDrcpYjBnvlM35/pTR8pOqxbvfPUOr5e+jqqo3D7hdlRFvi9CiO6hJ3sA+7J8+XIqKys5+uijE8dM02TBggU88sgjvPHGG0SjUerq6trM6lRUVJCfv/f+QR6PB4/H051DF92lfCVUrYOUrN17UimKc7xqnfM6VYeXrnfOnfgTOOwbPT7cJxduoqYpytCcVKb3w7o5lm1RUlNCXbiOTG8mRYGi/QYtteFa7lx8J+AkIBfn9ExFaiHEwalXBzqnnXYaq1atanPs8ssvZ8yYMfz0pz9l0KBBuFwu3nnnHWbMmAHAmjVr2Lx5MxMnTkzGkEV3a652uozrLYFqrNlpzKlq4PI5x8N1ULsJ3voFGCEYfhpM+VmPDzUYjvGXBc5szo9PG9nvqiAv2b6EuavmUhosxbAMdFVnaMZQZo6dyYSCCXt93z1L7qEmXMOIzBFcM+6aHhyxEOJg1KsDnfT0dA4//PA2x1JTU8nOzk4cnzlzJjfddBOBQICMjAyuv/56Jk6cyPHHH5+MIYvu5ssG1QWhWgjVObuosAEFdK+Tf6PosOhRqPvKqZMz43EnEOpilmWzuixITXOUgM9NcWFGmwadcz8opT4UY0RuGtPGFXb59ZNpyfYl3LnoTppiTfg9ftyam6gZZW3tWu5cdCezJs7aY7DzxqY3eH3T62iKxl0n3IVbcydh9EKIg0mvDnTa48EHH0RVVWbMmEEkEuH000/nT3/6U7KHJbpL/jhIzYHylpk+TcdJNbOc2Z1YM/gCsPVjZ4bngqed511s4foqZs/fwIbKRmKmjUtTGJ6bxjUnD2fSiBx2NEQSHcpvnDoKrR91KLdsi7mr5tIUayLXl4vSsoTo1b14NA+VzZXMXTWX8fnj2yxjbazbyC8W/gKQJSshRM9RbNu2kz2IZAsGg/j9furr68nIyEj2cMS+WBY8NsUJdBSlZaZGAWxnCcsynccA33oSis/v8iEsXF/F7S+sojFikOVz49ZUoqZFbXOMNI/G3eeP5Y3V5fxt0VeMG+jnxWtPSAQD/cHq6tXc8N4N+HQfXn33TQJhI0yz0cxDUx6iONsJZhqiDXz31e+yKbiJY/OO5S9f/wsu1dXTQxdCHIT6/IyOOMiUr3SK/6UXQLi+7dKV6gLLcF53zGXdEuRYls3s+RtojBjkZ3h3zmaoGvkZKuXBCA+8tZYVW+oA+OmZY/pVkANQF67DsIzEslPYCCdydLy6F7fmJhgNUheuA5wZoNs/uJ1NwU3k+fL43cm/kyBHCNFjJNARfUs8GTktz8nXCdeCEXNmdxpbepzpHhhzdrdcfnVZkA2VjWT53LsFMIqikOlzsbqsHsOyOXnUACYNz+mWcSRTpjcTXdWpj9QTjAaJmJFErOnRPGS4M9BVnUxvJrZt88CyB3h/6/u4VTd/mPIHslOyk/0lCCEOIhLoiL5lT8nItgW26ZxXXZCSDakDuuXyNc1RYqaNW1OxsQlHLQzLQldVvG4Vy7IJxSwU4Kdn9GzNnp5SFCgi4AnwZc2XAOiajqqoWFiEYiFCsRBjAmMoChQxe+Vs/vbF3wD4xaRfSF6OEKLHSaAj+pZdk5FVDVqnmVkmpOU6r+sGAZ8bl6ZQF4pRH4oRMUxs25lQ8ugaMdMpVHjy6AEcVtjP870UUFCwbAvbtrEVG0VRsFtypJ74/Almr5wNwE/H/5RvDO/5OkZCCCGBjui7FMUJbOyWKshK9/eQKi7MIDvNzeqyIAqgayqK6sRazREDCyc1+s5v9NzMRWeK9h2IkpoSaiI1ZHoyqYvUYVots2k2qIqK3+1nc8NmHvrkIQB+fPSPufgwadgphEgOCXRE39I6GbmpCqwwAJaiUuJLp86bRma4iqLtn6Ieckz3jiWeomODbVtYShRsN/60MIWZXdeyZF86W7TvQNSF62iONROKhcAGXdVRcGZyTMukLlKXmNW56oir+MHYH3TLOIQQoj36V6lW0f/Fk5FtA0wnyNmm61ydl8MNgVR+nqZyQzpc/fHdLNm+pMsvv7osSHVjlAK/lxSXhmXbWN514P8AbDeoYYzcv/C9167oluu3Fi/at7Z2LT7dR05KDj7dlyja113X93v8hI0wpm3i1tzoio6maGhoiYAH4PLiy7nuqOu6ZQxCCNFeEuiIvsWXDUYUGpzu9FWazpV5A1jrcuGzLHJiUXwWrA1VdMsv+3gycmaKm0OzfWRnb8aT+yqx+qMASAksQ1FNNjWs79ZgY9eifV7di6qoeHUvub5cmmJNzF01Fyu+rNeF4oGMoijEy3BZtkXEimDhXE9XdL425Gtdfm0hhOgoCXRE3+LLgUgdADYKswZk06Sp5FoWXtv5D9qLTW5qQbf8sm+djPxVdRNB75uEa4/FNtNQXLW4Mlah2m6yPQO6NdgoqSmhNFiK3+NHURTCRpjGaCNhI4yiKPg9fkqDpZTUlHT5tYORYCKwilkxomaUqBVNnNcVHZ/LRzAS7PJrCyFER0mOjug7Ig3w9/MSu6zWeDysc2n4TQtllwLfihFu88s+XqH3QLVORlY9W3ETI1br5MK4s9/DVBrxaC5S3OmoWtdfPy5etC9mxagIVuxWyybgDWBYRqJoX1fK9Gbic/nQFK1NPo6CQoqegt/jx8Ym05vZ5dcWQoiOkhkd0TeYBvzn+1C9HhQVUgdQraoYgNu2AJuwotCouQgrCrRU7u2uX/YAit5MpHoKoKH51qFnLEfRa4mpO9jc8BUxK9atwYZlW2xv3E7YCKMqKrrq1LMJG2G2N27Hsq3uCTZspxpybaS2TZCjKiq2bdMUa2JoxlCKAkVdf20hhOggmdERvZ9tw+s3w7o3QfOAngKhWjJdOjpQr2oENZVIolKxgidcRQZmokJvV2mdjFwdySYWygZMPHmvoCgKCs5e83iwkeHJ6JZgY3TWaEzLdBKCVTdKyxYwFRVFUYhaUUzLZHTW6C65nm3brNixgsc+e4wPtn3Q5pymaOiKjolJ2AwTs2JMKpzUrVvchRCiveQnkej95t8Hy/4KKDD9MVCd+jlFBgQsmwpdI6QoqLaNbtuoNoTMKBVNFQQ8gS6dWYgnI/u9LuyY097BlbUYzV2NqqgtbSEUFDRM2+zSYKO1NbVrUFUVTdUwbAPDNjAtM/FYUzVUVWVN7ZoDuk5DtIF/fvlPvvnyN/ne69/jg20foCoqOSk5pOqp+HQfiqJgtlSmTtFT8Ok+FpYt7JbcJCGE6CiZ0RG929LH4P17nMdn/RayBjuFAVXd2WJux5dOwFnAUrCVln7m3dBLM56MvKMxQsSwUdRmPDnvgWIllnHAGVbrYKM7cnQ0RSPLk0V1uHq3on0BTwAbu1PLZvWReuZvnc9bX73Fwm0LE4nGHs3DOcPO4cSBJ3Lv0nvJS83Do3kIm2FMy0RTNbyal4gZ6bbcJCGE6CgJdETvteo/8NrNzuOTb4XjroD1bzttH/yDKInWUKOpZJomdZqG2TKbAjYqNpmeLGoiNV2ejDwwkMLS0loA3NnzsW0VJxPY6XGFAvYBBhv7E8/RCUaCuxXtsyyL2nBtu5fNmmJNLK9YzsflH7O0fClf1nzZZjZmROYIvjnqm5wz7Bz8Hj8fbfso0b3cWa5zIkoFBUVRduteLoQQySSBjuid1r8DL1wN2DD+B3DKrc7xeFNPzUWdv5BmtZFQy8yNTjzMcUKOhkgDKa6ULv2Fq6oKwZDhPFEi6P5lKFrIeW5r2ICmKmB3LNjoqN1ydFp1UrdVe485OrZtsyO0gw11G1hTs4Yva7/ky+ovKQ2W7rbMNCJzBF8/9Ot87dCvMTxzeJvP70j3ciGESDYJdETvs3UZ/PsSpwJy8XQ4836nrxU4zTpzRkLFavyuHMIKmIA7/l7bBlUFzU3UjBJu2WbeVd4pqeDL8gYAUjQ3imKDYoHthFkKYNsKbk3r8oTg1nbN0VFtdWcbBttMJCT//KOfE4qF2NK4ha0NWwkZoT1+3sC0gRxXcBzH5R/H+Pzx5Ppy93rtjnQvF0KIZJNAR/QuO9bA09+EWBMMmwLn/9lZqopTVZh8I7xyA3ZzNXg8iVkcJd5GXNWxbTsxC9E6d+ZAhKImt/3X6Zqe6XMxILCDCkXFRnOCHZwlLNuGmGU7u44U+LjiY0ZmjkzMwMSTlC3bcpKHLYOIGSFshJ3gzAwnnkfMyM6bEaEh1kBDtIEtwS00RJyAy7RNTMy2g7UhZIR4deOrbQ6risrAtIGMDoymKFDE6MBoxgTG7DOw2auW7uXYu1RL7qLvtxBCdAXFtu2D/qdSMBjE7/dTX19PRkZGpz/nL5/9hf+t/1+bY7v+0N/Tt7s9vxja+75OXW8Pl++yz96Dvb7PMiESdLqRqzq2J33nTM6u77MMDCNEc+slF2XP2cdezYuuto3p9/j17TKuXV8TMy1ipnM9VVGwsQCj7Yd0QwJ0RygouFRXYmv3aYeexhE5RzA4YzCD0gdRmFqIS3Md0DVWV6/mhvduQEHZ69KVjc1DUx6SZGQhRNLJjE4Xqg3Xsrlhc7KH0bepCqABNkTb0UJgL8FNa2EzzK4THp0VLw3Tmb8OdEVHUzWnAaaiObuyWgr9eTUvbs2NV/Pi0T27Pfdozi3dnU66O500Vxp//+LvbGnYgmVbxKwY2KCoCrqqo6s6xdnF3D357i6vZxOvypyTkkOmJ3O3XVc2NlWhKklGFkL0ChLodKHvFn2Xrw/5eptjSjv+xFf28Mt61/ft6XN2fd8er6Xs+rQd12rHePakPePZ47Uayp3E44YyyDwU5bw5kJrT+k17/BzbtvnFwl+wrm4dhum0Q7CxE7MauqYzMnMkd55wp1Pjph3j2dO1bAuu++enrNpWz+QROexoiLCusgEFG0/Bcyie7YAFigk4S2Ypugdd1Tks+zDmfG3ObjNKXaE2XMsfPvkDpm0mAigTk5ARQlO0bivaF09GjppRvLqXFD2lzfmIEZFkZCFEryGBThcalD6IQemDkj2MviVYBv+9Fmo3Q9ZQuOQVyChs99tPG3wan1d9vvOXPc4v+6gVxbRNTht8GkP9Qw9oiH9ZsIHPNmmkeXK5Z9qJXPP0cuyYCxSwmorRvJtbEpI14nk6YTOMZmlMPmRytwQ5lm2xsGwhPt2HYRuJrxfFKdqnKzoLyxbyveLvdXmwUxQoYmjGUNbWrsWjedru+LJt6iP1jMoaJcnIQoheQSojiwNjWVD2qVPfpuxT53l7NVTA374BtaWQORgufblDQU7rX/bxWYWurtC7vrKB3725FoA7zimitjmWaAGR4lJQfGvBdoPl7PtSFOdaHq17KwTHu5fn+HI4NONQ8nx5ZHmzyPPlcWjGoeT4crqte7mqqMwcO5NUVyqVzZWEjTCWbRE2wlQ2V5LqSmXm2JnSAkII0SvIjI7ovI3z4cMHoWqdsxVcdTlbvyffCMNO3vd7G3fAU9+A6nWQMRAufQUyOzYb1vqXfXdU6DVMi588u5KoYXHK6AF8+9hBLFhXRcy0yU334Esrp9JbBZYfBReoIWzbwLI0clKzcbvMnule3ty2e3ldpK5bu5cDTCiYwKyJs5i7ai6lwVKC0SC6qjMqaxQzx85kQsGEbrmuEEJ0lAQ6onM2zodXboBII6Rkge4BIwIVq53j5zy092CnoQKeOhd2fAnphXDZy5B1aIeHEP9l310Vev+8YCMrt9aT7tW5d/oRKIqSaAERNS1sVxOKYoJqYalV2MRo6cFAdSRMjpbdI93LbexEYrOF1e0NReMmFExgfP54SmpKqAvXkenNpChQJDM5QoheRQId0XGW5czkRBohvWDnzidXCuheaNjunB9yolP3prX6rc5yVc0G572XvgyBYZ0aRndW6P2yPMhDbztLVr+cVky+3ws4LSCG56ZRsr2BLHcqNhaWWoNzYRVsFUWBqBXpV93L90ZVVNlCLoTo1eRPL9Fx5Sud5aqUrN23dyuKc7xqnfO61mo3wRNnOkGOfxBc/hrkjOj0MOIVeiuaKgjFQomt2qqiEoqFOt29PGpY/N9zK4mZNlOLcpl+9CGJc6qqcM3Jw0nzaFTXZbfU0jFxghzne+HSVHRF79Hu5ZZtYdt2oghhV3UvF0KIvk4CHdFxzdVOTo7u2fN53eOcb67eeaxqPfz1TKhr2V11+eudnsnZjdKyTd1uKfpntzzvZPG+B95ay+fbgmT6XNw9fexu2+Ynjcjh7vPHMrSgvqWhp4azvdzGpTl1mrs72Ih3Ly9ILcCre7GwMG0TCwuv7qUgtQBN0aSWjRDioCdLV6Lj4o01jYizXBVrdqoaqxq4fM5x1eW8DmD7Z/CPGdBUCTmj4Xv/g4yCAx5GSU0JNZEa8nx5uy1deXUvGe6MDncv/2DdDubM3wDAvdPHkpvu3ePrJo3IwfQewm0fuHEpGTQadRh2FNu2sHGuH/AGCBmhbsvR0VUdl+picPrgPSZix9SY1LIRQhz0JNARHRdvrFn2KZgGmBESEYbmAU2HwqOc15UugH9+F6INkDcWvvdi22KAB6CrK/RWNUa46Vlnue27EwZzxuH7DsYCKVmkuNz4dC8DtEN7NNhoXcsm15fbpmif1LIRQoidZOlKdJyqwvBTIdroNN9UAEV37mNNzvHhp0LJS85MTrQBDp0Ml7/aZUEOtK3Q61QjTiHNnUaKnuIk5JrRdicjW5bNzc+tZEdDhJG5adxx9mH7fU882KiP1AO0uT5AfaSeoRlDuyXYkFo2QgjRPvJTUHScZcGGd8GdBu5UZzLHNpx7d6pz/NO/w3OXgRmFomlw8fPg9XfpMFoHGrs15GyZ1WhvoDF7/gbeW7MDt67yx+8eRYpb2+97kh1sxGvZjMoaRbPRTFWoimajmVFZo5g1cZbUshFCCKR7OdB13csPGmWfwr8udoIa3QtGCCwDVB00LwS3QqjGee0xl8PZv3fyd7rBku1LuHPRnTTFmvB7/Lg1N1EzSn2knlRXart+4X+0vopL5i7Bsp28nAuPG9zhMcQL5xmWga7qDM0Y2mOF8yzbklo2QgixFxLoIIFOh61/G178IaTlOe2848nIigJNVRDPiRn7bZj+l3Z1GD8QBxJobK8Pcc7DH1LdFOVbxwzk/m8escempvsjwYYQQvROkowsOi6+6ypUC6E6MMJOi2/bwlm/ArxZMPGH3R7kQOcr9EYMk2uf/oTqpihFBRn8+rzDOxXkgBTOE0KI3koCHdFx+eOcpOLyVc5zVW1p5hmfHFSdlg7543psSB0NNGzb5mcvfM4nm+tI9+rMufhovK7uWV4TQgiRPBLoiANkgxlreaw4S1k9xLJsVpcFqWmOEvC5KS7MQFXbNyMz98NS/rN8K6oCj3z3aA7NTu3m0QohhEgGCXREx5WvdLqPu3zO1vE4RXWOef1Ork75SqeeTjdYuL6K2fM3sKGykZjpVCQenpvGNScPZ9KIfW9hf+/LSu5+rQSAn599GCePGtAtYxRCCJF8ki0pOq5uMzTvaBXkKDtvtg26e/cWEF1o4foqbn9hFSXbg6R6dHLTPaR6dEq2N3D7C6tYuL5qr+9dXVbP9f/8FMuGC8cP4vIThnTLGIUQQvQOvT7Queeeexg/fjzp6enk5uZy3nnnsWZN295B4XCYa6+9luzsbNLS0pgxYwYVFRVJGnE/V7Ea3viZUx8HnFkc3e1sM9c0MMNOIGRZO1tAdCHLspk9fwONEYP8DC9el4aqKnhdGvkZHhojJrPnb8Cydt9MuLW2mcuf+JjGiMHxwwLceW7nk4+FEEL0Db0+0Jk/fz7XXnstixcv5q233iIWi/H1r3+dpqamxGtuvPFGXn75ZZ577jnmz59PWVkZ06dPT+Ko+yHbhuVPwmOnQv2WnbupVLcT7Ci03OvOVnPbcFo+dLHVZUE2VDaS5XPvFqQoikKmz8WGykZWlwXbnKtrjnLZEx9T2RBhdF46f77kWNx6r//PXwghxAHq9Tk68+bNa/P8ySefJDc3l+XLl3PSSSdRX1/P3LlzeeaZZzj11FMBeOKJJygqKmLx4sUcf/zxyRh2/xJpgJdvgM//4zwfdJzTjTzS4AQ0lgK2AortBESqDooGFau6PEenpjlKzLRxayo2NuGohWFZ6KqK163i0VTqLZua5mjiPc3R/2/vzsOjLM9Hj39nz2SZbJBJwp6wR0ABSQMotmApIhcuVeqJFgHt0YNI4FTrcpCqFZVa9FgRxcuidUN/vx7UUpAfpSrrDyIIsqQYAkqALLIkQ5hktvc5f7zJhLCoQCZvMrk/1zVX8i4zua/3CnDzPPfz3EHuevML9lXWkO6KYcmUK0l02po1LiGEEK1Tq090zlRdrfcVSklJAWDr1q0EAgHGjBkTvqdv37507dqVTZs2nTPR8fl8+Hy+8LHH4znrnjZD0/SiX+8xfaoofZC+3Lu5lO2A/5gCx0v05GX0HHDnwEf36UvMayr1XZEbmCwQ31HfUycCNTopsXZsFhNVtQGqawP4giGU0geYHFYLiU4bNrOJlFg7AHWBEL/561a++PYECTFW3ph6JZlJzh/4KUIIIaJFm0p0NE2joKCAESNGcNlllwFQXl6O3W4nKSmpyb1ut5vy8vJzfs7TTz/N448/HulwI2//57D+eTharBf/mm16V/GRsyBr1KV9thaCDS/Ap0/rn+3qDL/8C3TN1VtAaBrU1ScyFlvjiI6mwanv9A0DI1Cjk5PpIjXezu4jHkyA1WLGZNYHkmr9Qbz+IDmZLnIyXQRCGve9u431+44Sa7fwxpRh9E2Xna+FEKI9aVNFCtOnT2fXrl0sXbr0kj7n4Ycfprq6OvwqLS1tpghb0P7PYXmBXhxsj9PbMdjj9OPlBfr1i3V8PywZB2ue0JOcvtfDPev0JAf02hsV1JMhs62+x5VF/2q2RbRGp4mGEh11xjEQCGkULN3OP4sqcVjNvD75SoZ0S45sPEIIIVqdNjOic99997F8+XLWrl1L586dw+fT09Px+/1UVVU1GdWpqKggPT39nJ/lcDhwOByRDjlyNE0fyfHVQEJGY2GwzamvfjpZpl/vftWFTWM1FByvehQCp8CeANfNh0G3NW3lULFTn6IyW+u7llvQswwFKhTRGp3dRzwcq/GTkRhTP3WlhaeunDYLLqeNoyd9TF6yhf/efxybxcQrtw8hL7v5R5eEEEK0fq0+0VFKMWPGDJYtW8Znn31Gjx49mlwfMmQINpuNNWvWcPPNNwOwd+9eDh48SF5enhEhR175Dn26ypl8di8pk0k/f7T4wjbsO1YCf58J36zTj7uNgBsW6a0czuQ9pncjT+yiT1MF69CHVUxgdUJcRz1RikCNTkMxclqCg+RYO3WB04qRbWZCmqLku1OUeY5jt5p59fYh/LRvWrPHIYQQom1o9YnO9OnTeffdd/noo49ISEgI190kJibidDpJTExk2rRpzJ49m5SUFFwuFzNmzCAvLy96V1x5j+lTStbzjEpZHXoH8R+TaISC8N8L4dN5esJidcLPHoWf/C89mTmXhqaeFhukZEGwVi9INlv19wfrIGSLSI1OQzGyP6ThsDUdrQopxTfHvPhDGnarmdcnD+WqXrLrsRBCtGetPtFZtGgRANdcc02T80uWLOHOO+8E4Pnnn8dsNnPzzTfj8/kYO3YsL7/8cgtH2oIaEo2gT5+uCnjr62UseguGoE+//kOJxpEv9VGcsh36cY+rYcL/1ZOX75M+SC96rtitT53ZYhuvKaV3NXfnRKSpZ06mi+y0eL46VEUwpPCH9KkrUGgKNAVmE7xx55U/2ApCCCFE9Gv1iY5SZ+9we6aYmBgWLlzIwoULWyCiVqAh0TjypT4iE/IRnjqyOMBi1aeszpdonDoG/3oCtr6pvy8mEcbOg8vzz54KOxezWV/ZtbxArwdyJuujSEGfnuQ4EvTrzbnMPfyjTVzdqwObSo4R0hRWiwmTSREINd5z5/DukuQIIYQA2tiqK1HPbIbsn4G/Rq+FMaHvSGxCP/bX6NfPTDS0EBS+Di8N0YuOUTDgVpheCFfc/uOSnAZZo+D6F/SRG/8pqKnQv7pz4PrnL315+3lommJtsb5cPNZuQVONSY4JiLdbKa6sOWcLCCGEEO1Pqx/REeegaVDyL7DH67UxQR8QBEz6EnOzVb+ed19jslPyL/ivx/SVUADuy+C6P0K34RcfR9YofWVXJDcsPENDCwi3K4YaX5BT1XqW47Ca6ZEaS1ARbgExoHNixOIQQgjRNkii0xY1rLqKd+vLyc9VDNyw6spshdWP6YkOgCMRfvZ/YOhUfYrrUpnNzb6E/Psc9/rxBzW8fh9VtQEA4uwW3K4YrFYzFo2zWkAIIYRovyTRaYtOX3UVnm6q/2oy6ee9R+Gfv6/fOFDpxcnD7oarfgtxzbcaStMUu494OO71kxJrJyfThdkcuY7gXl+QqtoAwfqpKbNJb/NQesJ7zhYQQggh2jdJdNqihlVXtSegtqrpPjaW+n/gA6dg/2f69zk36T2qfmg11QXauO8oiz4voaSyhkBIYbOYyE6L595R2REpBv5o+2Ee/n87G5McwGYxYzKduwWEEEIIIYlOW5Q+SG+oWV5fb2Ox6nmOFtATnAbZo+Gnj0LnIc0ewsZ9R3lk2U5qfEGSY+3YLWb8IY2ispM8smwn824c0GzJzsm6AH9YXsT7X+itOuLsFrz+EKaGUqDTW0BIDbIQQojTyKqraBAK6EvMldZ4LrUX5P9nRJIcTVMs+ryEGl+QdFcMMTYLZrOJGJuFdJeDGl+IRZ+XNMvKp7Vff8cvXljH+1+UYjLBr67sgivGSmZSDE6bvuoqqCk0pXDaLKQnxnCsxs/uI224I70QQohmIyM6bVHZl1BVqhcah3ynXahvwRCbAoHaC2sBcQEaVj4lx9oxnbEk3WQykRRru+SVT566AE+dNorTJcXJszcPJBBSrCmqJC3OTlKsjWpvkEBIw2YxkxhrBWWissYnxchCCCEASXTaFr8XdrwH656D2tPbOzQUIlvAZAarXd/TJgK9pqCx35TdYkahqPOf1m/KbsZhMV/0yielFCt2lvOHf+yhrLoO0DcAfGBsH+IcVnYeqsZmMVFVG6hv6hkKN/WsqpViZCGEEE1JotMWVOzWdzH+6n29hxXQmNyY65eJmwENQnVQdRBikiPSawoa+02dK9m4lJVPe454ePzvu9l84DgA3VNjmf/LQQzrkRK+JyfTRWq8nd1HPJgAq8WMySzFyEIIIc5NEp3WylcDu/4G296Ew1sbzyd1g2H/E9bOh7pqffVVePrIrH8f8oEKgntAREJr7mTj+Ck/f/qvvby35SCa0jf/u2dUNveMysZpP09jUQjnelKMLIQQ4nwk0WlNtBB8sw52/gfs/lBv5QB6LU6f62DIZMj6KZR/Beut+nkVBM0EygQmpWcbZqs+jVWxM/Kb+V1CslHl9fPauv28seEbTvn1HY7HD8zgkev60SnJec737D7i4ViNn4zEmPrRJC08muS0WXA5beFiZNkZWQghhCQ6RlMKDn0Bu/4Tdi/Te0Y1SMnWk5tB/wPiOzae9x7TO5XHdYCaSn1X5AYmi36v0iJWo3OpyUZ1bYDX1x9gyfoDnPTpsedkuphzfX9+kvX9020N9UFpCQ6SY+3UBU6rD7KZUQopRhZCCBEmiY4RtBAcKoR//wP2fARV3zZecyZD/4l6s81uw8/daDM2Ve93VVefyFhsjSM6mganvotojc7FJhuHq2p5c+M3vLf5YDjB6ZuewKxre/Pz/u6zVnCdS0N9kD+kEWOz1E9tNU5v1QVDUowshBAiTBKdluL36jsV7/0H7P1Eb9HQwBYHfcfDgF/qU1PWH/hH2j2gfsoqBBaH3gchTEW8Ruf0ZMNhO2MrJhP4glqTZOPLgyd4ff0BVu4qJ1S/t05vdzwFY3rzi5z0C2oZkZPpIjstnqKyk6S7zE2SI6UUVd4A/TISpBhZCCEEIIlO5CgF3/0bSj6F/Z/CgXV6880GjkToda2e4PQeq3cd/7EqdupTVA01OspCuDhGhSJeo9OQbHx1qIpgSOEPNU5d2S1mrBYT/TNcbDt4goeXfcWuw42b943omcq0kT24pnfaRfXEMptN3Dsqm0eW7aTc4yMp1obDYsYX0qjyBoh3WLh3VHZE+20JIYRoOyTRaU413+mjNiX/0pObk2VNryd20YuK+14H3UboU04Xo6FGJ7GLPk11eq8rqxPiOuqtICJUo2M2m7i6Vwc2lRwjpCmsFhMWM4Q0xSl/CBNQ+M0JNu3Xl4nbLWYmDMpk2sge9G+GkZbhPTsw78YB4T5b1ZrCZjbRLyMhYn22hBBCtE2S6DSnFf9br7lpYI3R62yyfgrZPwN3zrlrbi5UQ1NPi01v1Bms1QuSzVY90QnWQcgWsRodTVOsLT5KrN1CMKThC2kEQ42LrRQQ1BS93fFMurIrN17RiZS45q2ZGd6zAz/JSm3RzulCCCHaHkl0mlP2aDh+ALJ/qic3XfPAFtP8Pyd9EHTopW8kmJABttjGa0rpXc3dOfp9EbDlwHF2Ha5GAb6QInRaiy2L2USc3YLdYua5Xw5iYJekiMQA+siSLCEXQgjxfSTRaU6Df60vB480sxlGzoLlBfr0mDMZrA4I+vQkx5GgXzc3T89Wf1Djq0NVbD5wnHXF31F44AQh1bhZjtkEsXYLrhgbyXG2cL+pE7WBZvn5QgghxMWSRKc5Nce01I+VNQqufwHWPw9Hi/XWEGabPpIzcpZ+/SJ5/UF2lFaz+cAxthw4zraDJ6gLaE3usdZ3Kw9qGsGQRl1Awxf04akLSr8pIYQQrYYkOm1Z1ijofpXepdx7TK/JSR90QSM5dYEQRWUedh6u5qtD1ew8VE1x5Um0M3Y3To2zM6xHCj/JSmVU747MeG9b0xYQJuk3JYQQovWRRKetM5t/1BJyXzDEgaOnKK6oobjiJF9X1FBceZJvjnnDe9uczu1ycGX3FHKzUvlJjxR6psWH96zRTr9f+k0JIYRoxSTRiRJKKaprA3x7zMvB4/qrtP7rt8e8lFXXnjVK06BDvJ0BnRIZ0DmJgZ0SGdA5Ebfr/EXU0m9KCCFEWyGJTisX0hSe2gDHvX4qPT4qT9ZR4amjwuOjwlNHpcdHRf25M+tozpTgsNLLHU9vdwK93An0StO/d7scP6r9QoPTW0Akxdqo9gYJhDRsFjOJsdZwMbL0mxJCCGE0SXQiTNMUp/xBvP4QNb4gp3zB+q8hTvmCnKwLUOUNcMIboMrr54TXH/6+qjZAdW0AdQFTQW6Xg64psXRJiaXrGa+OCReW0JxPQwuIhvh8wVB4RKeq1iLFyEIIIVoNSXSa0eK1JXy84winfI1JjdcfapbPjndYSXM5cCfE4HY5cLtiSHM1fu9OiCHN5SDGZvnhD7tEOZkuUuPtTYuRzVKMLIQQovWRRKcZVXh8Tfo6na5hI714h5VYh5U4h5V4h36cHGsnKdZOcqyNpFhb/ff6cWKsjSSnHbu1efbEaXZSjCyEEKIVk0SnGU26sgsje3Ug3mElzm6tT2r0ZMZhNTfLtFFrIMXIQggh2gpJdJpRb3cCvd0JRocRcacXIyfH2qkLaAQ1DavZTIzNjFJIMbIQQohWQRIdccEaipH9IY0YmwWn3QI01gbVBUNSjCyEEKJVaKWFH6I1y8l0kZ0WzwlvAHXGkjClFFXeANlp8VKMLIQQwnCS6IgLZjabuHdUNvEOC+UeH7WBEJqmqA2EKPf4iHdYuHdUNmZzdNQkCSGEaLsk0REXZXjPDsy7cQD9MhLw+oJU1vjw+oL0y0hg3o0DGN6zg9EhCiGEEJjUmXMP7ZDH4yExMZHq6mpcLpluuRCapth9xMNxr5+UWDs5mS4ZyRFCCNFqSDGyuCRms0mWkAshhGi1ZOpKCCGEEFFLEh0hhBBCRC1JdIQQQggRtaIm0Vm4cCHdu3cnJiaG3NxctmzZYnRIQgghhDBYVCQ677//PrNnz2bu3Lls27aNQYMGMXbsWCorK40OTQghhBAGiorl5bm5uVx55ZW89NJLAGiaRpcuXZgxYwYPPfTQD75flpcLIYQQ0anNj+j4/X62bt3KmDFjwufMZjNjxoxh06ZN53yPz+fD4/E0eQkhhBAi+rT5ROfo0aOEQiHcbneT8263m/Ly8nO+5+mnnyYxMTH86tKlS0uEKoQQQogW1uYTnYvx8MMPU11dHX6VlpYaHZIQQgghIqDN74zcoUMHLBYLFRUVTc5XVFSQnp5+zvc4HA4cDkdLhCeEEEIIA7X5ER273c6QIUNYs2ZN+JymaaxZs4a8vDwDIxNCCCGE0dr8iA7A7NmzmTx5MkOHDmXYsGG88MILnDp1iilTphgdmhBCCCEMFBWJzqRJk/juu+947LHHKC8v5/LLL+eTTz45q0D5fBpW2MvqKyGEEJGWkJCAyWQyOox2Iyr20blUhw4dkpVXQgghWoTs2dayJNFBr+k5cuTIBWfZHo+HLl26UFpa2q5/aeU56OQ5yDNoIM9BJ89Bd+ZzkBGdlhUVU1eXymw207lz54t+v8vlatd/iBvIc9DJc5Bn0ECeg06eg06egzHa/KorIYQQQojzkURHCCGEEFFLEp1L4HA4mDt3brvffFCeg06egzyDBvIcdPIcdPIcjCXFyEIIIYSIWjKiI4QQQoioJYmOEEIIIaKWJDpCCCGEiFqS6AghhBAiakmic4meeeYZTCYTBQUFRofSon7/+99jMpmavPr27Wt0WIY4fPgwt99+O6mpqTidTgYMGMAXX3xhdFgtqnv37mf9PphMJqZPn250aC0qFAoxZ84cevTogdPpJDs7myeffJL2uObj5MmTFBQU0K1bN5xOJ8OHD6ewsNDosCJq7dq1TJgwgczMTEwmEx9++GGT60opHnvsMTIyMnA6nYwZM4bi4mJjgm1HJNG5BIWFhbz66qsMHDjQ6FAMkZOTQ1lZWfi1fv16o0NqcSdOnGDEiBHYbDZWrlzJnj17+NOf/kRycrLRobWowsLCJr8Lq1evBuCWW24xOLKW9eyzz7Jo0SJeeuklioqKePbZZ5k/fz5//vOfjQ6txd11112sXr2at956i507d/Lzn/+cMWPGcPjwYaNDi5hTp04xaNAgFi5ceM7r8+fP58UXX+SVV15h8+bNxMXFMXbsWOrq6lo40nZGiYty8uRJ1atXL7V69Wo1atQoNXPmTKNDalFz585VgwYNMjoMw/3ud79TI0eONDqMVmfmzJkqOztbaZpmdCgtavz48Wrq1KlNzt10000qPz/foIiM4fV6lcViUcuXL29yfvDgwerRRx81KKqWBahly5aFjzVNU+np6eqPf/xj+FxVVZVyOBzqvffeMyDC9kNGdC7S9OnTGT9+PGPGjDE6FMMUFxeTmZlJVlYW+fn5HDx40OiQWtzHH3/M0KFDueWWW0hLS+OKK67gtddeMzosQ/n9ft5++22mTp3a7hoXDh8+nDVr1vD1118DsGPHDtavX8+4ceMMjqxlBYNBQqEQMTExTc47nc52OfILcODAAcrLy5v8m5GYmEhubi6bNm0yMLLoJ009L8LSpUvZtm1b1M83f5/c3FzeeOMN+vTpQ1lZGY8//jhXXXUVu3btIiEhwejwWsz+/ftZtGgRs2fP5pFHHqGwsJD7778fu93O5MmTjQ7PEB9++CFVVVXceeedRofS4h566CE8Hg99+/bFYrEQCoV46qmnyM/PNzq0FpWQkEBeXh5PPvkk/fr1w+12895777Fp0yZ69uxpdHiGKC8vB8Dtdjc573a7w9dEZEiic4FKS0uZOXMmq1evPut/K+3J6f9DHThwILm5uXTr1o0PPviAadOmGRhZy9I0jaFDhzJv3jwArrjiCnbt2sUrr7zSbhOd119/nXHjxpGZmWl0KC3ugw8+4J133uHdd98lJyeH7du3U1BQQGZmZrv7fXjrrbeYOnUqnTp1wmKxMHjwYG677Ta2bt1qdGiinZGpqwu0detWKisrGTx4MFarFavVyueff86LL76I1WolFAoZHaIhkpKS6N27N/v27TM6lBaVkZFB//79m5zr169fu5zGA/j222/55z//yV133WV0KIZ44IEHeOihh/jVr37FgAEDuOOOO5g1axZPP/200aG1uOzsbD7//HNqamooLS1ly5YtBAIBsrKyjA7NEOnp6QBUVFQ0OV9RURG+JiJDEp0LNHr0aHbu3Mn27dvDr6FDh5Kfn8/27duxWCxGh2iImpoaSkpKyMjIMDqUFjVixAj27t3b5NzXX39Nt27dDIrIWEuWLCEtLY3x48cbHYohvF4vZnPTv1YtFguaphkUkfHi4uLIyMjgxIkTrFq1iokTJxodkiF69OhBeno6a9asCZ/zeDxs3ryZvLw8AyOLfjJ1dYESEhK47LLLmpyLi4sjNTX1rPPR7Le//S0TJkygW7duHDlyhLlz52KxWLjtttuMDq1FzZo1i+HDhzNv3jxuvfVWtmzZwuLFi1m8eLHRobU4TdNYsmQJkydPxmptn3+1TJgwgaeeeoquXbuSk5PDl19+yYIFC5g6darRobW4VatWoZSiT58+7Nu3jwceeIC+ffsyZcoUo0OLmJqamiaj2gcOHGD79u2kpKTQtWtXCgoK+MMf/kCvXr3o0aMHc+bMITMzkxtuuMG4oNsDo5d9RYP2uLx80qRJKiMjQ9ntdtWpUyc1adIktW/fPqPDMsTf//53ddlllymHw6H69u2rFi9ebHRIhli1apUC1N69e40OxTAej0fNnDlTde3aVcXExKisrCz16KOPKp/PZ3RoLe79999XWVlZym63q/T0dDV9+nRVVVVldFgR9emnnyrgrNfkyZOVUvoS8zlz5ii3260cDocaPXp0u/7z0lJMSrXDLTuFEEII0S5IjY4QQgghopYkOkIIIYSIWpLoCCGEECJqSaIjhBBCiKgliY4QQgghopYkOkIIIYSIWpLoCCGEECJqSaIjhBBCiKgliY4QUeaaa66hoKDA6DCEEKJVkERHCPG93njjDUwmEyaTCYvFQnJyMrm5uTzxxBNUV1cbHZ4QQnwvSXSEED/I5XJRVlbGoUOH2LhxI7/5zW/461//yuWXX86RI0eMDk8IIc5LEh0hotiJEyf49a9/TXJyMrGxsYwbN47i4uIm97z22mt06dKF2NhYbrzxRhYsWEBSUlKTe0wmE+np6WRkZNCvXz+mTZvGxo0bqamp4cEHHwzf5/P5uP/++0lLSyMmJoaRI0dSWFjYJJ78/Hw6duyI0+mkV69eLFmyJHy9tLSUW2+9laSkJFJSUpg4cSLffPNNRJ6NEKJ9kERHiCh255138sUXX/Dxxx+zadMmlFJcd911BAIBADZs2MA999zDzJkz2b59O9deey1PPfXUj/rstLQ08vPz+fjjjwmFQgA8+OCD/O1vf+PNN99k27Zt9OzZk7Fjx3L8+HEA5syZw549e1i5ciVFRUUsWrSIDh06ABAIBBg7diwJCQmsW7eODRs2EB8fzy9+8Qv8fn8Eno4Qol0wuHu6EKKZjRo1Ss2cOVN9/fXXClAbNmwIXzt69KhyOp3qgw8+UEopNWnSJDV+/Pgm78/Pz1eJiYnh4yVLljQ5Pt2iRYsUoCoqKlRNTY2y2WzqnXfeCV/3+/0qMzNTzZ8/Xyml1IQJE9SUKVPO+VlvvfWW6tOnj9I0LXzO5/Mpp9OpVq1adUHPQAghGsiIjhBRqqioCKvVSm5ubvhcamoqffr0oaioCIC9e/cybNiwJu878/j7KKUAfWqrpKSEQCDAiBEjwtdtNhvDhg0L/7x7772XpUuXcvnll/Pggw+ycePG8L07duxg3759JCQkEB8fT3x8PCkpKdTV1VFSUnLhD0AIIQCr0QEIIdquoqIiXC4XqamplJWV/eD948aN49tvv2XFihWsXr2a0aNHM336dJ577jlqamoYMmQI77zzzlnv69ixYyTCF0K0AzKiI0SU6tevH8FgkM2bN4fPHTt2jL1799K/f38A+vTp06RYGDjr+HwqKyt59913ueGGGzCbzWRnZ2O329mwYUP4nkAgQGFhYfjngZ60TJ48mbfffpsXXniBxYsXAzB48GCKi4tJS0ujZ8+eTV6JiYkX/RyEEO2bJDpCRKlevXoxceJE7r77btavX8+OHTu4/fbb6dSpExMnTgRgxowZrFixggULFlBcXMyrr77KypUrMZlMTT5LKUV5eTllZWUUFRXxl7/8heHDh5OYmMgzzzwDQFxcHPfeey8PPPAAn3zyCXv27OHuu+/G6/Uybdo0AB577DE++ugj9u3bx+7du1m+fDn9+vUDID8/nw4dOjBx4kTWrVvHgQMH+Oyzz7j//vs5dOhQCz45IUQ0kURHiCi2ZMkShgwZwvXXX09eXh5KKVasWIHNZgNgxIgRvPLKKyxYsIBBgwbxySefMGvWLGJiYpp8jsfjISMjg06dOpGXl8err77K5MmT+fLLL8nIyAjf98wzz3DzzTdzxx13MHjwYPbt28eqVatITk4GwG638/DDDzNw4ECuvvpqLBYLS5cuBSA2Npa1a9fStWtXbrrppvAy9rq6OlwuVws9MSFEtDGphmpCIYQA7r77bv7973+zbt06o0MRQohLJsXIQrRzzz33HNdeey1xcXGsXLmSN998k5dfftnosIQQolnIiI4Q7dytt97KZ599xsmTJ8nKymLGjBncc889RoclhBDNQhIdIYQQQkQtKUYWQgghRNSSREcIIYQQUUsSHSGEEEJELUl0hBBCCBG1JNERQgghRNSSREcIIYQQUUsSHSGEEEJELUl0hBBCCBG1/j+wSE6sGV089wAAAABJRU5ErkJggg==", | |
"text/plain": [ | |
"<Figure size 584.486x500 with 1 Axes>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"refDose = np.linspace(min(drData.logDose)*0.9,max(drData.logDose)*1.1,256)\n", | |
"refDose = (10**-refDose)*1e6\n", | |
"sns.lmplot(data=drData,x='logDose',y='response',hue='compound',fit_reg=False)\n", | |
"for fit in fitData:\n", | |
" plt.plot([pDose(i) for i in refDose],[ll4(i,*[fit[i] for i in ['b','c','d','e']]) for i in refDose])" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Python 3 (ipykernel)", | |
"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.10.12" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 4 | |
} |
@jhfc Thanks for your feedback, and sorry about the late answer: those functions belong to the numpy package and were incorrectly called. The source is now fixed.
I am getting this error. I checked and found that there is no plt.plot in seaborn.
sns.plt.plot([pDose(i) for i in refDose],[ll4(i,*[fit[i] for i in ['b','c','d','e']]) for i in refDose])
AttributeError: 'module' object has no attribute 'plt'
import matplotlib.pyplot as plt is missing from the imports
and
sns.plt.plot needs to be modified to plt.plot
for fit in fitData:
plt.plot([pDose(i) for i in refDose],[ll4(i,*[fit[i] for i in ['b','c','d','e']]) for i in refDose])
Thanks! Fixed the code as you suggested, should work now @akashbahai @cbmII
How would you go about implementing a 1/Y^2 weighting and returning their parameter values as well as visualising the line on the graph?
Hey @Fae14, this link should give you a few pointers:
https://stackoverflow.com/questions/27696324/using-scipy-optimize-curve-fit-with-weights
Hope this helps!
Nice attempt. May be code should be more explanatory and have more comments,
Copied as-is into a new jupyter session, I get the following after assigning fitCompound
:
[/var/folders/2s/9l3t4cfj0c5g1s4_7shzc5900000gp/T/ipykernel_19131/1118045928.py:7](http://localhost:8888/var/folders/2s/9l3t4cfj0c5g1s4_7shzc5900000gp/T/ipykernel_19131/1118045928.py#line=6): RuntimeWarning: invalid value encountered in log
return(c+(d-c)[/](http://localhost:8888/)(1+np.exp(b*(np.log(x)-np.log(e)))))
Then a type error making the lmplot:
TypeError Traceback (most recent call last)
Cell In[81], line 3
1 refDose = np.linspace(min(drData.logDose)*0.9,max(drData.logDose)*1.1,256)
2 refDose = (10**-refDose)*1e6
----> 3 sns.lmplot('logDose','response',data=drData,hue='compound',fit_reg=False)
4 for fit in fitData:
5 plt.plot([pDose(i) for i in refDose],[ll4(i,*[fit[i] for i in ['b','c','d','e']]) for i in refDose])
TypeError: lmplot() got multiple values for argument 'data'
Copied as-is into a new jupyter session, I get the following after assigning
fitCompound
:[/var/folders/2s/9l3t4cfj0c5g1s4_7shzc5900000gp/T/ipykernel_19131/1118045928.py:7](http://localhost:8888/var/folders/2s/9l3t4cfj0c5g1s4_7shzc5900000gp/T/ipykernel_19131/1118045928.py#line=6): RuntimeWarning: invalid value encountered in log return(c+(d-c)[/](http://localhost:8888/)(1+np.exp(b*(np.log(x)-np.log(e)))))
Then a type error making the lmplot:
TypeError Traceback (most recent call last) Cell In[81], line 3 1 refDose = np.linspace(min(drData.logDose)*0.9,max(drData.logDose)*1.1,256) 2 refDose = (10**-refDose)*1e6 ----> 3 sns.lmplot('logDose','response',data=drData,hue='compound',fit_reg=False) 4 for fit in fitData: 5 plt.plot([pDose(i) for i in refDose],[ll4(i,*[fit[i] for i in ['b','c','d','e']]) for i in refDose]) TypeError: lmplot() got multiple values for argument 'data'
I have the same problem, anyone found a solution?
When I try running it, in the fifth block I get
NameError Traceback (most recent call last)
in ()
3 # generate base curve
4 curData = pd.DataFrame(data={'compound':curve['compound'],
----> 5 'dose':curve['startDose']/power(curve['dilution'],range(curve['nDose']))})
6 curData['logDose'] = pDose(curData.dose)
7 curData['response'] = curData.dose.apply(lambda x: ll4(x,*[curve[i] for i in ['b','c','d','e']]))
NameError: name 'power' is not defined
any idea what's going on here?