Skip to content

Instantly share code, notes, and snippets.

@yannabraham
Last active March 13, 2026 14:16
Show Gist options
  • Select an option

  • Save yannabraham/5f210fed773785d8b638 to your computer and use it in GitHub Desktop.

Select an option

Save yannabraham/5f210fed773785d8b638 to your computer and use it in GitHub Desktop.
How to do Dose/Response curve fitting in Python for Drug Discovery
Display the source blob
Display the rendered blob
Raw
{
"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>6.087886</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.956018</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.712320</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>54.164809</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>78.656826</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" compound dose logDose response\n",
"0 A 10.000000 5.000000 6.087886\n",
"1 A 3.333333 5.477121 12.956018\n",
"2 A 1.111111 5.954243 28.712320\n",
"3 A 0.370370 6.431364 54.164809\n",
"4 A 0.123457 6.908485 78.656826"
]
},
"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": [
"/var/folders/2p/m_28_h056155fhzfq4bryffr0000gp/T/ipykernel_38712/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[0]\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>6.383218</td>\n",
" <td>106.383218</td>\n",
" <td>0.40</td>\n",
" <td>367.625170</td>\n",
" </tr>\n",
" <tr>\n",
" <th>B</th>\n",
" <td>0.7</td>\n",
" <td>5.613085</td>\n",
" <td>91.613085</td>\n",
" <td>1.30</td>\n",
" <td>240.023681</td>\n",
" </tr>\n",
" <tr>\n",
" <th>C</th>\n",
" <td>2.0</td>\n",
" <td>30.172700</td>\n",
" <td>158.172700</td>\n",
" <td>0.02</td>\n",
" <td>354.063751</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" b c d e residuals\n",
"compound \n",
"A 1.0 6.383218 106.383218 0.40 367.625170\n",
"B 0.7 5.613085 91.613085 1.30 240.023681\n",
"C 2.0 30.172700 158.172700 0.02 354.063751"
]
},
"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>6.383218</td>\n",
" <td>106.383218</td>\n",
" <td>0.40</td>\n",
" <td>367.625170</td>\n",
" <td>1.0</td>\n",
" <td>0</td>\n",
" <td>100</td>\n",
" <td>0.40</td>\n",
" </tr>\n",
" <tr>\n",
" <th>B</th>\n",
" <td>0.7</td>\n",
" <td>5.613085</td>\n",
" <td>91.613085</td>\n",
" <td>1.30</td>\n",
" <td>240.023681</td>\n",
" <td>0.7</td>\n",
" <td>0</td>\n",
" <td>86</td>\n",
" <td>1.30</td>\n",
" </tr>\n",
" <tr>\n",
" <th>C</th>\n",
" <td>2.0</td>\n",
" <td>30.172700</td>\n",
" <td>158.172700</td>\n",
" <td>0.02</td>\n",
" <td>354.063751</td>\n",
" <td>2.0</td>\n",
" <td>24</td>\n",
" <td>152</td>\n",
" <td>0.02</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 6.383218 106.383218 0.40 367.625170 1.0 0 100 \n",
"B 0.7 5.613085 91.613085 1.30 240.023681 0.7 0 86 \n",
"C 2.0 30.172700 158.172700 0.02 354.063751 2.0 24 152 \n",
"\n",
" ref_e \n",
"compound \n",
"A 0.40 \n",
"B 1.30 \n",
"C 0.02 "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHpCAYAAACLCpbcAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAlfhJREFUeJztnQd8VGXWxp/pk15JAULvVQVEEEEFRUVX7LqoqKxtFeu6q679s+uurg17l7WsHRVERBBBqvQQOgRIIb1P/37nvZlhEgIkIZPMTJ6/e/fe+743MzcTMvPkvM85R+fxeDwghBBCCAlD9G19A4QQQgghgYJChxBCCCFhC4UOIYQQQsIWCh1CCCGEhC0UOoQQQggJWyh0CCGEEBK2UOgQQgghJGyh0AEgpYTKysrUnhBCCCHhA4UOgPLycsTFxak9IYQQQsIHCh1CCCGEhC0UOoQQQggJWyh0CCGEEBK2UOgQQgghJGyh0CGEEEJI2EKhQwghhJCwhUKHEEIIIWELhQ4hhBBCwhYKHUIIIYSELRQ6hBBCCAlbKHQIIYQQErZQ6BBCCCEkbKHQIYQQQkjYYmzrGyCEEBJauD1uZBZloqSmBPHWePRP7A+9jn83k+CEQocQQkijWZqzFG+tews7ynbA6XbCqDeie2x3TBs8DSPTR/KVJEEHJTghhJBGi5xHljyCzcWbEWmMRHJEstrLuYzLPCHBBoUOIYSQRi1XSSSn0lGJlMgUWI1WtVwlezmXcZmX6wgJJih0CCGEHBHx5MhyVZwlDjqdrs6cnMu4zMt1hAQTFDqEEEKOiBiPxZNjNpgbnJdxmZfrCAkmKHQIIYQcEcmuEuOx3WVX5zXOGlTYK9RekHGZl+sICSbaVOgsXLgQ55xzDjp27KhCn1999dVB12RmZuJPf/oT4uLiEBUVhREjRmD37t2++ZqaGtx0001ISkpCdHQ0LrjgAuTl5bXyd0IIIeGNpJBLdlVBdQF2lu7E7vLd2FexT+3lXMZlXq4jJJhoU6FTWVmJoUOH4uWXX25wftu2bRgzZgz69euHX375BWvXrsX9998Pq9Xqu+b222/Ht99+i88++wwLFizAvn37cP7557fid0EIIeGPGI9HdxyNKkcVqp3VasygM6i9nMu4zLOeDgk2dB6Px4MgQCI6X375JSZPnuwbu/TSS2EymfDBBx80+DWlpaXo0KEDZs6ciQsvvFCNbdq0Cf3798eSJUtwwgknNPh1NptNbV7KysqQkZGhHi82NrbFvzdCCAl1JJvqhrk3YEPBBjg9TtjddkA+PXSAWW+GUWfEwOSBePW0Vyl2SFARtB4dt9uN7777Dn369MHEiRORkpKCkSNH1lneWrlyJRwOByZMmOAbk+hPly5dlNA5FE888YRaCvNuInIIIYQcOesqOTIZXWO7IiMmAx2jO6q9nMs4s65IMBK0Qic/Px8VFRV48sknccYZZ+DHH3/Eeeedp5alZIlKyM3NhdlsRnx8XfNbamqqmjsU99xzj4reeLfs7OyAfz+EEBIuWVcSgY8wRiDaHK32cs6sKxKsGIM5oiOce+65yocjHHPMMVi8eDFeffVVjBs3rtmPbbFY1EYIIaFKa/eb8s+6kiKBkm3lbQEh58y6IsFK0Aqd5ORkGI1GDBgwoM64+G8WLVqkjtPS0mC321FSUlInqiNZVzJHCCHhSFv0m/JmXW0o3KCe8yCPjt6IgUkDm511JXZReVyH2wGXxwWX26W8QDLmf672tWPevbcas+w98p/Hc8S9G251/3Lu+zptoM6518baM74nBiTV/TwioUHQCh1ZkpJU8qysrDrjmzdvRteuXdXxsGHDlFl53rx5Kq1ckOsl/XzUqFFtct+EENIa/aak5YJUI5YlI4mmePtNPTDqgYCIHW/W1fLc5UpcGPQGlXUl4qPKWQUDDEiPSsecnXNQbi8/eHOUqyiQzWVDjasGNmft3mVTx7JXQiNIuXrQ1RQ6IUqbCh3x4GzdutV3vmPHDqxevRqJiYnKUHzXXXfhkksuwdixY3HKKadg9uzZKpVcUs0FMRJPmzYNd9xxh/oayZiaPn26EjmHyrgihJBw6TflbcUgS0cWgwX5VflqfkTaiBZbxpLU8d1lu7GrbBc+2/yZek75T6Ip8p8XF1z4YusXamtJJFIkGV2y94qr+udyT3rofa+H/7n6z7v3Hvude18n75icy17734Gv6Rqj/YFNQo82TS8XwSICpj5Tp07Fu+++q47ffvttlSW1Z88e9O3bFw8//LDy7fgXDLzzzjvx3//+V6WMS4bWK6+80qSlK0kvF9HE9HJCSDAjy0a3zb9NdQwXcVMfiZhIdOX5U55Xy0hNRSIvq/NXI6s4C5uKNiGrKEsJnMZEWryCo09CH6RFpiHGHFNnizZpxmURZBajBVaD9cDeYFHfj0SnvCJGNtbkIWFVR6ctodAhhIQCv+39Dff9dh+SI5KVCKhvCJaIj1QofvTER3FipxOP+HgSGVqZt1ItRy3LXabETUPdx2WJLN4Sj5zKHMSYYpQg8X50yHGUKUqJoaY8NyFo7x4dQgghDWc+ldpKUWYvU74WryFYoiKx5tgj9puSr5u/ez5+2vUTFu9brJmK/ZC6OIOSB6FvQl/0S+yHvol9lbDyRpNkKae5z01IW0ChQwghIYJkNCVaElXkRTAatOUdySCqdlSrTcRJ/cwnib6syl+Fzzd/jh93/aiJlFo6RXfCCeknYHjacAxPHY60qLQWfW5C2hoKHUIICTVqjbLe9Gg1pNMd5KWRVOy5u+bizXVvKt+Nl17xvXB6t9MxocsEdew18bbkcxMSLFDoEEJIiCAFAotsRUiNTD1o+Ug8OrJ8JPMbCzciuzwbL69+WZmJBTECn9X9LFzQ+wK1NNUkcdPAc4s/SMSNiB7/55brmmOEJiRQUOgQQkiItWEQz4yYg6UOjURtJM1aspdEeORW5uK+RfdhW+k2n5F4Sv8p+HO/P6vjo31uEUy+HJZakSXn4s+RVHS5jpBggkKHEEJChPptGNQSktIa2tLRvop9KtpSYitRguSaQdfgygFXItIU2SLPLRlZORU56rlEXEmtGvHoSGRJxmMtsTQjk6CDQocQQkKEQ7Vh8G9jIMgS1e3Dbj+ksbg5SBaWRI+kKrK0fPCKLG9hPrkXmZfrCAkmgrZ7OSGEkIbbMFQ5qtQykRcRH16Rc16v8/DU2KdaVOQIYmbW6/UqkiNtH5Swkp5RHrc6VxEevb6O6ZmQYIBChxBCQgQRFVL7RlVGNlh9jS69lYmjjdHKo9NQ0b+jRbw38hzSz0oVJ4Rba6gJtzqXcZmnR4cEG1y6IoSQEEEymqRjeWJEoqpC7HFpURwxJkvbBVk+kvlAZD55/UEmvQldYrocZIQWn45D76BHhwQdjOgQQkiIINESh8uhmndWOCqUTybJmoREa6JaNpJ2DOLdCURUxesPkqrMgpido81a/ypBxmWeBQNJsEGhQwghIYJkT5U7ylXjTq9nRwTG7nKtu7gcB6oNgzzXtMHTVF8rEVpSR0eWyGQv5zIu82zESYINCh1CCAkBJJLzxto3VGq5INEcWTbydvmWFgx5lXmqTUOgoioj00figVEPqA7lIrZk+Uz2ci7jMk9IsEGPDiGEBDmS3fTo0kfx695fD6R16/Rt0oZBxMyItBHKByRLZBI9EmHFSA4JVih0CCEkyPk462N8seULJXLEeCxZTodrARHoNgwiatjmgYQKFDqEEBLELM9djqeWPaWOz+t9HhbuWahaOYjgEbEjS1liQhaRIxEdWU5iijchB6DQIYSQIEXaKtzxyx2qXs2kHpNwUZ+LVB0dMR3Xj+hI2wcRO4EyIxMSqlDoEEJIECIZTff/dr8SMOKBeWjUQypyI2bjTUWb1DVGg2ZElqJ9YkaWrV9iP6Z4E+IHs64IISQI+STrEyzNXaqK8T0z7hnlwfGh08zHyowsncQ9teeaT5kQ4geFDiGEBBlSE+e5lc+pY2nO2TW2qzoWk7GYjVMjUxtswyDjXjMyIUSDS1eEEBJESFuF+xbdp5p2jkwbiUv7XeqbE5OxVD5OjkhWZuT6bRhoRibkYCh0CCEkiPh8y+dYvX+1qjT8yImP1KlP4+03JZlWEsHxtl/wYnPaaEYmpB5cuiKEkCCh3F6Ol/54SR1PP3Y6OkZ3PGS/KeXN8UPO2W+KkIOh0CGEkCDhtTWvodhWjB5xPXBx34sPmme/KUKaDoUOIYQEiQH5o00fqeO7RtwFk97U4HXsN0VI06BHhxBCgoB/rfiXMhqP6TRGbYeD/aYIaTwUOoQQ0sasyluF+dnzYdAZcNfwuxr1New3RUjj4NIVIYS0Ma+vfV3tJ/eajB7xPdr6dggJKyh0CCGkDVlfsB6/7ftNRXOmDZ7GnwUhLQyFDiGEtCFvrH1D7aVpZ0ZMBn8WhLQwFDqEENJGbC7ejJ+zf4YOOkZzCAkQFDqEENJGvLn2TbU/vdvpqnYOIaTlodAhhJA2ILs8G7N3zlbH1w6+lj8DQgIEhQ4hhLQBn2V9pppwntjxRPRN7MufASEBgkKHEEJaGZvLhi+3fqmOL+l7CV9/QgIIhQ4hhLQyP+78ESW2EqRHpWNs57F8/QkJIBQ6hBDSynyc9bHaX9jnQhj0Br7+hAQQCh1CCGlFMgszsXb/Whj1Rpzf+3y+9oSEs9BZuHAhzjnnHHTs2BE6nQ5fffXVIa+94YYb1DXPP/98nfGioiJMmTIFsbGxiI+Px7Rp01BRUdEKd08IIU3nk6xP1P60LqchOSKZLyEh4Sx0KisrMXToULz88suHve7LL7/E77//rgRRfUTkbNiwAXPnzsWsWbOUeLruuusCeNeEENI8KuwV+H7H9+r44r4X82UkJNy7l5955plqOxx79+7F9OnTMWfOHEyaNKnOXGZmJmbPno3ly5dj+PDhauzFF1/EWWedhWeffbZBYUQIIW3FvN3zUO2sRrfYbhiWOow/CELau0fH7XbjiiuuwF133YWBAwceNL9kyRK1XOUVOcKECROg1+uxdOnSQz6uzWZDWVlZnY0QQgLNd9u/8/W1kqV4Qkg7FzpPPfUUjEYjbrnllgbnc3NzkZKSUmdMrk9MTFRzh+KJJ55AXFycb8vIYCM9QkhgKaguwNJc7Q+wSd3rRqcJIe1Q6KxcuRL/+c9/8O6777b4Xz733HMPSktLfVt2dnaLPj4hhNTnhx0/wO1xY0iHIciI5R9XhKC9C51ff/0V+fn56NKli4rSyLZr1y7ceeed6Natm7omLS1NXeOP0+lUmVgydygsFovK0vLfCCGkVZatGM0hpP2YkQ+HeHPEb+PPxIkT1fjVV1+tzkeNGoWSkhIV/Rk2TDP2/fzzz8rbM3LkyDa5b0IIqc/O0p3YULgBBp0BE7tN5AtESHsROlLvZuvWrb7zHTt2YPXq1cpjI5GcpKSkOtebTCYVqenbV2uA179/f5xxxhm49tpr8eqrr8LhcODmm2/GpZdeyowrQkjQ8N0OLZozquMoJEXUfV8jhITx0tWKFStw7LHHqk2444471PEDDzzQ6Mf46KOP0K9fP4wfP16llY8ZMwavv/56AO+aEEIaj8fjwffbtdo5Z/c4my8dIa2MziO/he0cSS+X7CsxJtOvQwhpSbKKsnDhtxfCarBiwSULEGmK5AtMSCsStGZkQggJB37O/lntT+h4AkUOIW0AhQ4hhASQ+bvnq/2pGafydSakDaDQIYSQAJFTkYPMokzodXqMyxjH15mQNoBChxBCAsT8bC2ac0yHY5BoTeTrTEgbQKFDCCEBFjqnZJzC15iQNoJChxBCAkCZvQwrcleo41O6UOgQ0lZQ6BBCSAD4dc+vcHqc6BnXE11ju/I1JqSNoNAhhJBALlsxmkNIm0KhQwghLYzD7cBve39Tx/TnENK2UOgQQkgLs27/OlQ4KhBviceg5EF8fQlpQyh0CCGkhVm8b7Haj0ofpWroEELaDv4GEkJIC7Nk3xJft3JCSNtCoUMIIS1Iqa0U6wvXq2MKHULaHgodQghpQZblLoPb41Zp5WlRaXxtCWljKHQIISQQ/hwuWxESFFDoEEJIC+HxeLB4L4UOIcEEhQ4hhLQQu8t3Y1/lPhj1RgxPHc7XlZAggEKHEEJaeNnquJTjEGmK5OtKSBBAoUMIIS0E/TmEBB8UOoQQ0gK43C5ft3IakQkJHoxtfQOEEBKKSAp5ZlEmSmpKEG+NBzxQbR+iTdHol9CvrW+PEFILhQ4hhDSRpTlL8da6t7CjbAecbqcyH1sNVjV3bMqxMOgNfE0JCRIodAghpIki55Elj6DSUYk4SxzMBjPsLjuyy7PVfEpkCl9PQoIICh1CCGnCcpVEckTkiKDR6XRq3GKwqDkhqyhLHbOZJyHBAc3IhBDSSMSTI8tVEsnxihzB5rLB5XFBBx32V+9X1xFCggNGdAghpJGI8Vg8ObJcJdQ4a9S5RHiECGOEEjxyHSEkOKDQIYSQRiLZVWI8lg7lZfYyFcmRbCsRN4IsV8m8ysIihAQFXLoihJBG0j+xPxIticirzEO1o1oJG4POAI+oHUBFdmReriOEBAeM6BBCSFPRyf90KprjNSETQoITRnQIIaSRiMm4yFaE1MhUWI1WuOGus2yVFpWm5mlGJiR4oNAhhJAmmpEl66pLTBdkxGQowSMkWBLUuMzTjExI8EChQwghTTQjS4FASS+X/5QhGUCUKUqN04xMSHBBjw4hhDQSMRl3j+2ODYUbVOTGWz9H2F+1HyaDCQOTBtKMTEgQwYgOIYQ09g1Tp8fojqNR5ahCtbPal20l1Lhq1LjMsyoyIcEDhQ4hhDQSybBavG8xIo2Rqjigx6MJHVnCknMZl3lmYhESPFDoEEJIE1tAJEcmo2tsV9XjSkiKSFLnMi7zzLoiJHig0CGEkGa2gPAaka0GqzInyzizrggJLmhGJoSQZrSAKLYV+zw6OZU56jzWHMusK0KCDEZ0CCGkGS0gbE4tmqOHHga9QbWEkHG2gCAkuGhTobNw4UKcc8456Nixowr7fvXVV745h8OBf/zjHxg8eDCioqLUNVdeeSX27dtX5zGKioowZcoUxMbGIj4+HtOmTUNFRUUbfDeEkHaDzu9Qp7WCUHu/cUJIcNCmQqeyshJDhw7Fyy+/fNBcVVUVVq1ahfvvv1/tv/jiC2RlZeFPf/pTnetE5GzYsAFz587FrFmzlHi67rrrWvG7IIS0xxYQXiTzSlpBSIVkGWcLCEKCizb16Jx55plqa4i4uDglXvx56aWXcPzxx2P37t3o0qULMjMzMXv2bCxfvhzDhw9X17z44os466yz8Oyzz6ooECGEtLQZOd4S7/PnJFgTlMgRf46MFVQXsAUEIUFESHl0SktLVXhYlqiEJUuWqGOvyBEmTJgAvV6PpUuXHvJxbDYbysrK6myEENJYM3JhdeGB9yVbKfKq8rC7fLc6ZgsIQoKLkBE6NTU1yrNz2WWXKT+OkJubi5SUlDrXGY1GJCYmqrlD8cQTT6iIkXfLyMgI+P0TQsLHjCwZVoJUQBZhI3uakQkJTkJC6Igx+eKLL1Zr4TNmzDjqx7vnnntUdMi7ZWdnt8h9EkLaF1IRWVVHphmZkKDFGCoiZ9euXfj555990RwhLS0N+fn5da53Op0qE0vmDoXFYlEbIYQ01YxcWFOoCZxaj45q6qmDz6fjNSNLc09CSNtjDAWRs2XLFsyfPx9JSUl15keNGoWSkhKsXLkSw4YNU2MihtxuN0aOHNlGd00ICWczslRDFpEjYicjJkP1tZI6OlIdmWZkQoKPNhU6Uu9m69atvvMdO3Zg9erVymOTnp6OCy+8UKWWS9q4y+Xy+W5k3mw2o3///jjjjDNw7bXX4tVXX1XC6Oabb8all17KjCtCSEDMyN5IjkRwIk2RdealiCDNyIQEFzqPt/1uG/DLL7/glFNOOWh86tSpeOihh9C9e/cGv06iOyeffLI6lmUqETfffvutyra64IIL8MILLyA6OrrR9yFZV2JKFr+O/9IYIYT4I9Gbif+biNyqXCRYEpAene6bk7fS/Kp89Enog1dPe1UZlAkh7TyiI2LlcDqrMRpMojszZ85s4TsjhJCDEfESbY4GqgCH24EaZ41q5Gl32VVqeZQpCtMGT6PIISSICGqPDiGEBBNSLHBvxV513DOuJ/Kq81BmL1PLVRLJEZEzMp3+QEKCCQodQghpJNtLt6PaWY1IYyTePfNdZBVnKYOyeHekxg6XqwgJPih0CCGkkWwo2KD2A5IGqCgOU8gJCX7oliOEkEayoVATOhQ4hIQOFDqEENJINhZuVPuBySwGSEioQKFDCCGNwOFyYFPRJnU8KGkQXzNCQgQKHUIIaQRbSraolPIYcww6x3Tma0ZIiEChQwghTfTn6HQ6vmaEhAgUOoQQ0oSMq0HJXLYiJJSg0CGEkEYgHckFqZdDCAkdKHQIIaQRFZG3FmsNiPsm9uXrRUgIQaFDCCFHYFfZLtjddkQYI5ARk8HXi5AQgkKHEEKOQFZRltr3TujNNg+EhBgUOoQQcgSkp5XQN4HLVoSEGux1RQghISJ03G4PNuwrQ1GVHYmRZgzsGAu9nqnuhBwOCh1CCDkCm4s2t7kRefHWAsxYsA3b8ivgcHlgMujQMyUaN47ridG9ktvsvggJdrh0RQghh6Gopgj7q/f7PDptJXLu/XIdMnPKEGUxIiXGovaZOeVqXOYJIQ1DoUMIIY0wIneJ6YIoU1SbLFdJJKfC5kRarBVWk0EtV8k+LdaCCptLzct1hJCDodAhhJDDsLm4bZetxJMjy1UJkeaDWk/IeXykSc3LdYSQg6FHhxBCGhHR6ZPQp03MwPJc4skxG/TwwIMauxtOtxtGvR5Wsx4Wgx6lbo+6jhByMBQ6hBDSiIwrT006pr6zrNXNwCKo5LlKqh0orXbA5nTB45FoDmAxGhAXYYJJr1PXEUIOhktXhBByCBwuB7aXblfHMxc528QMLFGjpGgzckqrUW13Qq/TwWjQqb2cy7jMy3WEkIOh0CGEkEMgIkf6XBk8kaiqjml7M7B3lcz7VCyhQ8gRodAhhJAjLFvBno7ESEubmIHlcQsr7EiPsyLCZIDb44HT7VF7OU+Ls6p5mpEJaRh6dAgh5AhGZNg7tpkZ2GtGluUyybyqcfg9v0mv/Dr5FTaakQk5BBQ6hBByhIiO0dWpzczAXjOy3eWGxVQvCK8DbE53WJmRZQnQG7GSvcvlgUsdu+GS89pN5jweuQ7qWo/f3ncM75h33v8a7WvlKu846j2Wthqp7bsnR6F/On1QoQiFDiGENIB8EHpbPyQau2JrabWyxBgNeuj02oepmIGr7E5lBA6UGVgeV7K71u4pgdPlUYLHK7QkyiTG5CGd45v9/A6XG9UOF6rtLlTZtX21w4lqu1t9b9452cu8CCu7062+zn9vc7nhkLl643a5Z6dLRaU0gSKCBXC5JTLl8QkbJWBqRUYwcv24HhQ6IQqFDiGENIC0fSi2FUOv08PkTgdQ07AZOMAfzGJ8Hts7GUu2FSoxIMLGoNeiEyI8DHodRvVIRE5ZDYoq7CistKGo0q6iT2XVTpTXOFBe40S5re55mdocSoyEAlKuSL5Xtem0vb72WESf+KXkxyHZaHrvudp7x7R577VyjYxpj61d6/1aeB/D93g6dE6IbOuXgDQTCh1CCDmMPyc9oguK90KZgbWlqwMRFTEDx0aYfGbgwZ3jWjSiJGJkX3E1vl69D0a99mEs0Q9/cSXnT8/ZrLajQYRDpMkAq9mASLNBfW8RtXs5lywzOZa92aiHyaBXe7NBV7vXw1S79577Xyd7WYLzpsfL96OO9XoYDAfEixqv3fsLG3ZpJ82FQocQQg7jz0mL6I7CAJmBZXlnb3E1dhZWYndRFXYWVGFXYaU6zymtURGbRr+Z63Wqnk5ilAVJUWbERZoQa5XNiBirUQky2cdYTAeOa89F0IgIqZ9VRkg4QKFDCCGHiej0ju+DzUdpBpbloR0FldiUW4bNeeXIyi3HlvwK7CmuVstRhyPKbECN062iKUqMyH86qAhJpFmvoiLF1Q7866KhOLlvCn+WhNSDQocQQg4T0Tmp6xBsyrI02gwsURoRMquzS7AmuwTr9pZi236tbURDWIx6dEuKQpekSHRLikTXpCh0TYpUnhApULg1vwLXf7BCRVu0pTOn7/krbFrWl6S4J0VZ+HMkpAEodAghpB41zhrsKtuljvsl9cPY3uWHNAOLaVUqJD/+fSb+yC7B+r2lKspTnxiLEX3SYtAnNQb90mLQOzUaPZKj1ZLY4fwn3hYQ4gFqi6wvQkIdCh1CCKnH1pKtcHvcSLAkIMmSjIVbdihDrggdqaPjqNUxElWRQM1nK/fW+XrxvhyTEY+hneMxNCMe/dNj0Ck+4ug9MG2Q9UVIqEOhQwgJKaTuikQ3xPwrvhiJZLR0Ro7Xn9MnsQ825pRjc265yjaSCI7UgKmvLSQzaEL/VEwcmKqETfekqBa7J/8WEK2Z9UVIuEChQwgJGaRLuDTQlN5S4nkRc64U07txXE+M7pXc4v6cmspU3PLxH8grt9WZl+eV7uXRFqNKyRYz8J9HdsG4Ph0QyBYQ0lertMqpfEBiRo6LNAIeHVtAEHIYKHQIISEjcu79ch0qbE6V5i1GYDEGZ+aUq/HHzxt81GJHCu19ty4HX21cDhiA3zOtcJZVqjmrUYSFGH8N0OlEZEnNGb1KN2+NFhANtaAoqQ58CwpCQh0KHUJISCxXSSRHRI5kInm9Lla9AWmxeuSW2dT8CT2SmrxkJIX5ft9ehJnLdmPO+lzYXS5E99mj7C+DOvTDhacMwrdr96nU8IoaJ4pc9hZtwXAkaEYm5Oig0CGEBD3iP5HlKonk1Df0yrks6ch8U3wq0v7g0+XZSuBs369FbYQ+nZzIMdTAqDPis2smw2Qwqd5PK3cVH7IFg7RoaJXKvTQjE9Jk6lW/al0WLlyIc845Bx07dlRvVl999dVBf2k98MADSE9PR0REBCZMmIAtW7bUuaaoqAhTpkxBbGws4uPjMW3aNFRUVLTyd0IICSRen4pEUKTbtDSZlJ5NspdzqSPjcHsaVZ04v6wGT/yQiROf+BmPfpepRI5kVF12fAa+vXkM7j5Xi8z0jO+pRI5EkxZuKVDXyCaIIVnQxoxqXq4LBP5mZDEfe7t6y17O0+KsPjMyISTIIjqVlZUYOnQorrnmGpx//vkHzT/99NN44YUX8N5776F79+64//77MXHiRGzcuBFWq1VdIyInJycHc+fOhcPhwNVXX43rrrsOM2fObIPviBDS2j4Vi7FxPhWpTPzqL9vw5R97lbdH6JUSjatGd8O5x3REjNWkxmas0YzIfRP71okmpcZaYTbqDjID25yeJkeTmmtGbukWFIS0B9pU6Jx55plqawiJ5jz//PO47777cO6556qx999/H6mpqSryc+mllyIzMxOzZ8/G8uXLMXz4cHXNiy++iLPOOgvPPvusihQ1hM1mU5uXsjL+JURIMHM0PpXc0hr8Z94WfLoi29duYXjXBNwwridO7Zdy0JLT5iKtOWafhD51hIaIo5xS+0FmYLmvxkaTjkbkyfOrxpoqqqRFloQap4tmZEKCdenqcOzYsQO5ublqucpLXFwcRo4ciSVLlqhz2ctylVfkCHK9Xq/H0qVLD/nYTzzxhHos75aRkRHg74YQElCfSgOUVNnxxPeZGPfMfPx32W4lck7p2wH/u2EU/nfjaEwYkNqgr8abWu6N6IjQkGUiab5Z43D5um/LXs5lXJatApX1JOJNUuiLqxzqD0B/5LykyqHmWRmZkBATOiJyBIng+CPn3jnZp6TUbWJnNBqRmJjou6Yh7rnnHpSWlvq27OzsgHwPhJDW96mIoHl/yU6MfXo+Xlu4XRXYG9EtAZ/dMArvXH08hndLPOTzVDoqkV2uvR/0TdCETv+0GLg8HvW4YkIWgSP/yV7OZVzm5bpAIGJM6gRFWwwqu6za4VLCSvZyLuMy3ypmaEJCkHaZdWWxWNRGCAkNGutTWbazCHd/sdZnzJWeUv84ox9O7tuhUe0XthRryQ4pESlIsCao48zcchh0slymg7SwMuo9vs4L6lyysHTadYGqTCz1gaROkLdYYqnbo5arpLVESxdLJCTcCFqhk5aWpvZ5eXkq68qLnB9zzDG+a/Lz8+t8ndPpVJlY3q8nhIQ+/j4Vi6leIFqnpXlX1jjxf7M2qqFYqxF3TeyLP4/sqtK/m9P6wV9k6XV6dIyLQGGlrV4LBq1reJXDFXAzsIgZqRMU6PYXhIQbQSt0JMtKxMq8efN8wkZMw+K9ufHGG9X5qFGjUFJSgpUrV2LYsGFq7Oeff4bb7VZeHkJIeOD1qazdUwJnrTHYKzZkCcku4qP22guO64x7zuqH5OimR219/pzaZSt/kWU26tEtKeqgaFKN0w2Ty90qlYlF1LCfFSEhJHSk3s3WrVvrGJBXr16tPDZdunTBbbfdhkcffRS9e/f2pZdLJtXkyZPV9f3798cZZ5yBa6+9Fq+++qpKL7/55ptVRtahMq4IIaGHvrYo35Jthb6ifXqd+HQAR61BNynKjFevGIYRh/HgHIn6RmR/kSWtJtJiLXWynrxmYFlCohmYkOCkTYXOihUrcMopp/jO77jjDrWfOnUq3n33Xfz9739XtXakLo5EbsaMGaPSyb01dISPPvpIiZvx48erbKsLLrhA1d4hhIQP/kX7ROhISrW3aJ8ghQT7psVgWJeE5j+Hx+3z6PhHdLxmYOmnJeZfqcIsBQptLrcSOTQDExLc6Dz18xXbIbIkJmnmkoElFZYJIcHFuj2luP6DFUroSL8rERyCUa9D54QIVVenyubEa1cMb/bSzq6yXTj7y7NhMVjw+59/h1FvPHTn9FozcCA6pxNC2olHhxBCvIj5VkzAInLKapxqTKohd4qPUGZjifhIJtLRGIK9RuRe8b0OEjkCzcCEhCYUOoSQoEcyqqT1g9TOkRwjqUYcZTYoE7LVrC0jHakFRHP8OfWhGZiQ0INChxAS1MzdmIe//W+NT+RYjHolesQfI1lX4s8Rc/KQzvFHZQj2pZbXtn4ghIQHQVsZmRBC3l60A9d9sELVyemaFAkpGSNLWIJUJRZUDR2bS2VlHU1NmYZSywkhoQ+FDiEk6JDMqoe+2YBHZm1U9XIuG5GhTMdRFqMyJGvXaNfKeaTZqLKyxKvTHEptpcitzD2oWCAhpJ0LHbvdjqysLFWNmBBCWgLx3dzy3z/w7uKd6vzuM/vhsuO7YMf+SqTGWtEtORKpMVZVN0f2cp4Sa1HZUN7WD01lc7HWsbxTdCfEmpl5SQjau9CpqqrCtGnTEBkZiYEDB2L37t1qfPr06XjyySdb+h4JIe2EKrsTf3l/Bb5bl6OqEb/052Nxw7ieKK52qF5XUhF5V2E18sprUFhpV3s5d7jcKuW7uVlX9OcQEr40S+hI9+81a9bgl19+qVO8b8KECfjkk09a8v4IIe2EshoHrnhrGRZu3q86kr81dQTOHqJVOJdsKulUvre4GjUOl2r7oFVH1qlzGZdlq+ZmXTUm44oQ0o6yrr766islaE444YQ6HYElurNt27aWvD9CSDugvMaBK99ahtXZJaoh5ztXj8CwrgdaOfRPi4HL41HeHbNREziC2uk9sDs9al6uO5qIzhGNyG43kLsGqCoEIpOAtKGSc96s5ySEBLHQ2b9/P1JSUg4al3YN/sKHEEIaJXLe1kSOFAH86C8jMahT3erGmbnlMOigojiSdGXUa6nmYj1W5wadmpfrmloZ2el2YlvJtiMLne0LgEXPAQVbALcD0JuA5N7AmNuBHuP4gyYkSGnWnyLDhw/Hd9995zv3ips333xTdRQnhJDGUGlz4qp3luOP3YcWOYJ4b/Q6PTrGRSDCpFfLWFJXR/ZyLuPS6645Hp2dpTthd9sRaYxEp5hOhxY5s24D8jYA5iggOlXby7mMyzwhJHwiOo8//jjOPPNMbNy4UWVc/ec//1HHixcvxoIF/IUnhBwZm9OFGz5ciZW7itVy1aFEjiDeGzEnm416dEuKQo3DDafbDaNeD6tJjxqnGyaXu1keHa8/RwoFiphqcLlKIjm2CiAmvXa9DIApAjBagfIcbb7bSVzGIiRchI50EV+9erXKsBo8eDB+/PFHHHfccViyZIk6J4SQwyFemzs+WYNfazuSv3fN8YcUOYJUPJYGmpk55UiLtSBC1dLR6ulIX2Kpktw/PaZZlZGPaEQWT44sV0UkHBA5XuRcxmVerut4bJOfnwQYKcTkcQNuZ+3mOrD3uOqN1577xv3GYjsCid3542pPLSB69uyJN954o2XvhhAS9ogwuf/r9b4U8teuGIZjuyQc9muk4rF0Cb/3y3Wqc3l8pAkWg9bjSkROtMWg5ptTGXlz0ebDt34Q47F4cowW7dxRpX346Q2AKVIbrynRrmvPYsJlBxzVgLNGe40cNbXn1dq+/pxcrzaH33FDY7XHTpvfeL15r1hRAsVfpNSOtQQn3gac9nDLPBYJfqGzatUqmEwmX/Tm66+/xjvvvIMBAwbgoYcegtnc/MZ6hJDw5pVftmHm0t0qGPL8JcfipN4dGvV10j388fMGY8aCbao4oHQrl0aeEskRkSPzzeGIER3JrhLjcXUxUF2ifVgrG7S4o61ARLw2L9eFEiI25HuqKQVs5YC9XFueU8eyL9PO1bF3rkw7t1ceEC9eMSNRk1BDliqlU71sOoMmXtVWf8wIRB7IAiTtQOhcf/31uPvuu5XQ2b59Oy655BKcf/75+Oyzz1Qxweeff77l75QQEvJ8u2YfnpmjCYuH/zQQk4akN+nrRcyc0CNJVUAW47F4cmS5qrk9rgqrC1FQXQAddOgd37vhiySFPCoZyF2nnRvkbVO8PO7a6EQVkDZYu66toikiVir3AxX52l6JsvpbSd1zESeBEg8q0mXVfExeL5OMmawHomAGC2A0AwbvZqp3bGlgrPbY/+tEhMiczitQ/MTJkcaYJdwuaJbQ2bx5M4455hh1LOJm3LhxmDlzJn777TdceumlFDqEkINYuasId362Rh1PG9MdV47q1qxXSURNU1PIjxTN6RLbBZHyAXwkfB+MngPnIjQCgZigRbSU7QFK9wLluUBFHlCZD1Tsrz2uFTcuW/OeQz70rbGAJQawxALmaMAiW0ztcYzfcfSBayTjzCyCJcJPvNTulehgmRES4kJH1tjd8ksI4KeffsLZZ5+tjjMyMlBQUNCyd0gICXn2FFfh2vdXqj5Wpw1Ixb1n9W/+g7Vg0b4j+nMEea7KAi3jSiInapnGo32Yywe9NU6bb6oZ2eUESncDRduB4l1A6R6gbK8makTclO3TPCiNxRwDRKcAUR20ZRYxSast3u/Yf0vURAxFCQlypAvDKaecguLiYsTHx7eO0JE6Oo8++qhq+SDp5DNmzFDjO3bsQGpqanMekhASplTbXbj+g5UoqrRjUKdY/OfSY2Bo5lJTSxft8/lzDlco0GtGliiGyuCRQdnXRnNkGUU8Kw2ZkcUoW7RDEzPerbj2vGS3ZpY9LDogJg2I7QTEpmv1e6JSNEGjRI3fsSwREUJaRuiIB2fKlCmqFcQ///lP9OrVS43/73//w+jRo5vzkISQMESiv3d/sVZ5aqTb+GtXDEekuZnJnt6ifWKKlYiE+DwkE8dbtO/s55ssdhrV40qiRhJFEmEiGCSt3aR5dFw12rglXhM7G78B9m8C8jO1vVeQHQrxoUjKckI3IC4DiBNB07l2L1tHbSmIENJsmhXvHTJkCNatW4fS0lI8+OCDvvFnnnkG7733XvPvhhAS9EjzzHV7SrFg8361l/ND8daiHfh69T4VwXl5ynHoFN/MqEP9on0SvVCm1wjtXMZlvnZJvTHYXXbsKNlx5IhO6mDAU5uurBORpquty+LWIjoSlakuAD69Evj0CmD+Y8CGL4D8jZrIMUVpZuUB52qRp3NeAKbOAm7fCPwzF7hpKfDnT4BJz2rzQy4Cuo4GErpS5LRTxBry9NNPqyCCxWJBly5d8Nhjj6k5+ew99dRTERERgaSkJFx33XWoqKjwfe1VV12FyZMnq8K+ssISHx+PRx55RBX3veuuu5CYmIjOnTurTGkvO3fuVB0OPv74YxWskGbdgwYNOqgAsJwff/zx6p7S09NVUpI8rpdu3bod5NEVP69kY3uR55EuCueddx4iIyPRu3dvfPPNN3W+5vvvv0efPn3U9yhLVnJ/bVJHR7Db7cjPz/f5dbzID4UQEn4s3lrgS+92uDyqDo4U8msovXvFziI88cMmdXz/pP4qW6rZBKBo3/bS7XB6nIgxxyAtKq3hi+S9LfMbbQlKkAjOoRAzbod+QEr/2v0AIKWfFqFh40/SBO655x5Vp+65555TBXpzcnKwadMm1U9y4sSJqtXS8uXL1efvX/7yF9x888149913fV//888/KzGzcOFClSQ0bdo01blg7NixWLp0qWrKLdnTp512mrrOiwghESpSKubf//43zjnnHGVJEUG1d+9enHXWWUpIvf/+++p+rr32WiWK/IVMY3j44YeVkJPgyIsvvqhWiHbt2qVEWHZ2tsrivummm5SIW7FiBe688862ybryvnD1w9Si1lyuFirQRAgJKpEjBfsqbE4kRJphNuhhd7lVtWIZlxo3XrFTUmXHLf/9Q1VAPveYjpg6unkZVoEs2uffsdzXjFiWwvb9AexeAuz+XdvkcRtEp3mExOwrRenOexXoffrRfZ+k3VNeXq7aKr300kuYOnWqr0CvCB4RPzU1NUpoREVFqTm5TgTJU0895fPIimB44YUXVP+3vn37KlEhpV/uvfden5CSzgaLFi1SmdJeRDBdcMEF6li8t7Nnz8Zbb72Fv//973jllVdUwpE8n/y+9OvXD/v27cM//vEPPPDAA+q5GouIpcsuu0wdS+RJ7nXZsmU444wz1PPK9/uvf/1Lzcv9SxRLvr9WFTpXX301jEYjZs2apcJX7FhOSHgjy1MSyRGRkxZr9f3OW/UGpMXqVbVimZeojUz97bO12Fdag+7JUXjsvMFH/x4RgKJ9Pn+OMRaY/ziwYyGwd9XBqdry+F7TsdNRG91xa6nZqqaLRZuXbCdCjpLMzEzYbDaMHz++wbmhQ4f6RI5w4oknqlWVrKwsn9AZOHBgHeGRmpqqlqK8GAwGFaWRiJA//k255TNeEo/kOb3PLfP+v8vy3LJstmfPniat5Ij9xYt8L7Gxsb57kecZOXLkIe+r1YSO9LlauXKlUnSEkPBHzMSyXCWRnPqiRc6lJYPMy3UrdhXhp8w8FfF58bJjEW05qhXyli/aV7AV2PYzsrK+UKd9V38GVFQemJdMpi4nAF1GAV1GAimDgLdPP/DcRr/nllTzsuq2LRhIwgrxpRwt0rmg/u+oqYGx+raTo0XElazs+ONwHGzGb4178adZ70Cyfsd6OYS0H6QKsXhyRLx44EGN3a97uFmv+k5JS4aVu4vx+PfaX4D3ntXvsI06m01Ti/bJElf2MiDrOyDrB6Bwq/rKrC6dVAZVX10EMHAi0ONkoNsYILFHXR+Q/xuweq5aE7LOE9iCgaRdIuZcETvz5s1T/ht/+vfvr7w44tXxRnXEg+Ndojpafv/9d+XjEcRkLAENWc7yPvfnn3/us6h4nzsmJsbn8+nQoYPyE3kpKytTHp+mIM9T35ws99XqQkfWymTNTtbWpA1EfXUmYShCSPggrRbEeFxS7UBptQM2p8tXM89iNCAuwgSDDnhtwTYliE4fkHr0vpyjKdonXputPwGbvgM2z67r3dGbkNdlOEp1e2HQ6dHzlvWHr0HjfW4pwldVVLf2jSxfyXhzCgYS0gBi7hXfi3zGSt9IWR7av38/NmzYoEy7kuks3h0xAMv49OnTccUVV7RIDbuXX35ZCS0RG2KElgJ911xzjZr761//qozK8nwifmSpTO7ljjvu8C2TSTaYCDHxDEm2l3h3ZJmsKdxwww3KnyPGaBF6Irb8jdatJnSkUKBQfw2RZmRCwhPpJ5UUbVZLU/K3nNGgV9ndojWq7U5U2Z1K7BRXOVQK+TMXDm1Z715jivZJivm2n4FlbwKZ3wK20gNfL0Ko90Sg75lArwnYXLAamHcTusf1gOVIhfbkuaVGjjSzFKSujTyvRHQk2iOeIWmP0J67l5MW5f7771ceGREKYvgVL6wIAEnHnjNnDm699VaMGDFCnYt5WDKkWoInn3xSbWJPkdR2iawkJ2sJBp06dVJp3yJAxCckhmdJSrrvvvt8Xy8mZ4ngSLeEuLg4/N///V+TIzri9ZHI0e23364ysiSdXYIqXsHVHHSe+gtqjaB+bn19pPdVKCHhNfmhSF0gRqMIadiMfO7LizShI0EUvV4JHnnzkCUsVVKmdmXnk+tG4fjuLdzpWTKhPrwIqCnWziXbSnwyku0kS1O+ztl+b2cxHbXaNf3O0vw2foX33lj7Bl744wWc1f0sPDX2CNkce1YB75yhRXLEfOyv3+TppE2DNIm8ejbQ+biW/K4JaRV27tyJ7t27448//vD1sQwnmhXRCTUhQwg5OkTgFFbYkR5nrV26cvtWjqxGPaod2vmfhnRseZFTv2ifVBMWRHioZSQ/cSP9mwZOBgZdqImbQ6S8NqoisheJ3PiozfSqc97QdYSQYKHZ6RAlJSUqv96beibpbBJaksgIISQ8zcgpMRaVeVXjOGBG3l9hg9vuVtWPzx7aMTA3kCcZT3qtGrJK/64vKsQsFANM+RToPKJJNXSOiKS0S60eR6UmrFQ0qTae5V/LR64jhIRHCwipVCgFfcSsVFRUpDZZI5SxVatWtfxdEkKCwowsBQL9AxpSV0ciPEJChEnV2GlxxB/zx4da4T5ZqqoTRZE2ENFAQnfNUCxG5SNQ5ajCrrJdjY/oSG0e8QZJ2rkYn2WdTkWT3Nq5jMt8E2r4EBJMdOvWTXlsw3HZqtkRHTEJ/elPf1JVGsUw5U1FE4f0bbfdpspOE0LCy4wsrR7W7imB0+VRgkd8O65azWEx6tC/Y6y6rsUo3Qssex1Y+a5fdWKpRmwApL+WWiqq/VvN28m8EWJjS8kWlSKfHJGstiMi9XGkQ7o0DxVBJW0gVGTHCBisQEUukDqQdXQICbeIjqS/eUWOIMeSDidzhJDwQq/XYWzvZFTaXKiyay1e/BePRPzIvFx31OSsAf43DfjPEOC35zWRIwJDunnDm2UlIsOidRJXRfv2aQUFG1G0b1Oh1n+rX2IjC56Kz0eabUpmlYgauQdztLaXc1kyk3n2syIkfISOZCbt3r37oHFpxiXFgwgh4YVEbxZuKUCk2YAosw4ZnhwVVBG5E2EEoiwmNX+4TuZHRNovzLwUeG0ssP5/WtSk6xjg0pnATcu1FgvifvYW7fNmWzUxjX1TcROFjtBjHHD281rkRpbSKvK0vZyf/Zw2TwgJn6WrSy65ROXPP/vss6qlu7dCouTXext1EULCrwXExMgsnG/7Gnc4tGaDfzbOx4XWFZhpvgCL8weq6wZ3bmJCwp6VwIIngS0/HvDdDDwfGD0d6HjMgfTyFira1+SIjhcRM91O0p5DaubIMplEkBjJIST8hI4IHCkGduWVVypvjiDVkW+88UZVbIgQEn5ZV0Mda3Gb43W8YT8NOZ4kdNCV4HLjfKS5CnF7zQxU6q5DUdWBZn2NiuDMf0yrYOwVOIMvBsb+TfPEBKBon9PtVB6dZgkdQUQNqx8TEv5CR8pSSxv5J554Atu2bVNjknElVRoJIeFHYoQRV3m+RLYrDh84tcro1xu+hcvtRo4uAam6Ylyl+xKREVcd+cGKdgA//x+w/vMDEZkhl2gCJ6lnw19jTdQad8pSVf2ifbraon32Ku26w7CzdCdsLhsijZHIiMlo/AtACAlZjqqtsAgb6WfhPSaEhCcDdTuRjxxc47gNLhhwgn4jTjBILRod3B4Pij1R6G7ch2TdTkk0b/hBZMlp4TPAsje0LClRKCJwTv6H1kizFYr2ef05klaulwgSISTsadZvuixXSS8OKQ4o+feyybH0vGioJXtzcblc6nmkNLV0c5WokfTO8O9aIcfSD0R6gcg10odryxYtNE0IaSGqCvGt83hs9HRDBGy4xfgFolEFK6R4H2CDEUapXNzQ0pHLCSx9HXjhGOD3VzSR0/NU4PqFwPmvHVnk+BftU6nlzgMdxJUp2dnoon1NKhRISBjidnuwbk8pFmzer/ZHlUDQBJYsWaIafE6aNAkhEdGR7qVffPEFnn76aYwaNcr3TUg31cLCQsyYMaNFbk66pMtjvffee6rysqSuX3311UpU3XLLLeoauYcXXnhBXSOCSITRxIkTsXHjRtUFlhBy9KwoMuEFx2R1fIvhCwzXZfl6Xdl1ZpTromH3GLCtyoI67ppdi4Hv7wLy1h9o5XDaw0Cvug2BG120TzbVvbzmQGTH273ce91hyCzSKrn3T+rftOcnJAxYvLUAMxZsU4kFUulcioBKfawbx/XE6F6NqCl1FEgnBdEOspdGpR07BqiKeks19RSh8fHHH+PMM8+sMy6dTSXrSppjtgTSAVVaz8sL40U6tUrk5sMPP1TRHHmx7rzzTvztb39T8/Lc8jXS1v3SSy9t8HFtNpva/Jt6ZmRksKknIYdg6ttLsWBzAQboduIL8wMqGcENHfTwwACXEj1Z6Ib9l83GuL6pQHke8ON9wLpPtQewxgPj7weGXV3bQqGJiOH4w/O1on3RaYcu2nf5F4fMgpL3i7GfjEWJrQQfn/0xBiYN5M+btCuRc++X61Q1c2njYjboVeHP4ioHoi0GPH7e4ICJnYqKCrXqIsGKBx98EEOGDMG9996LoF66slgsarmqPhJREaNySyGp6/PmzcPmzZvV+Zo1a7Bo0SKfwJL277m5uWq5yl+EjRw5UkWYDoWYqOU67yYihxDSMBv2lWLh5gJ1/KDpfRjgUb4cEQ7aXudzzSRGmIBVHwAvj6gVOTpN3ExfBYz4S/NETgsV7curylMix6gzold8L/64SbvB7faoSI6IHGnTYjUZVHFP2afFWlBhc6n5QC1jffrpp+jXrx/69u2Lyy+/HG+//XYdC0pQCp2bb75ZeWX8oyJy/Nhjj6m5luLuu+9WURl5gSR9/dhjj1UtJqZMmaLmReQIEsHxR869cw1xzz33qOiNd5NCh4SQg5E3o4e/3aiEzOmGleiIAtTABCM8MMOl9nKe605Eiq4Eg344H/jmZm15Kf0Y4Lr5wDnPA1Et0AfqKIv2bSrSjMjd47vD4u2ATkg7qoOVEGlW0Vh/5Dw+0qTm5bpAIKsyInCEM844Q33uLliwAEHt0fnjjz9UpKVz584YOnSoL9pit9sxfvx4nH/++b5rxctzNCrwo48+wsyZM5VHZ/Xq1UroyHLV1KlawbLmRqRkI4Qcnu/X5WLZjiJYDcCtxi/h9HgjMh6/vQcRqEGipwK6nCLNM3PqP4GRNwKGo0rsbNGifT5/TiL9OaT91cFyuDxquaohLAY9St0edV1Lk5WVhWXLluHLL7/0tYuSosMifk4++WS0Bs16F5KUcvHK+BOI5R+ptOyN6giDBw/Grl271NKTCJ20tDQ1npeXp9b/vMh5uHZhJaS1qHG48Pj3mji4uL8F6VsLEaOrVNlOkmLuUR4dN6JgQ4xei+5WdTgGkZe9CyR2D9yNNbNoHzOuSHslMdKsjMfiybE2sHxsc7lh0uvUdS2NCBrJ1PY3H0ukWIINL730krKPBKXQeeedd9AaVFVVQV/vLzVJT3OLMbHWEyRiR6JLXmEjxuKlS5eqKs2EkObzxsLt2FtSjY5xVow5ph/0W0XeuGDXmdS8EW4Y4fRlX1XDguWn/hfjElsvm6IpeJeumHFF2hsDO8aq7KrMnHKkxerrLF+J6CipcqB/eoy6riURgfP+++/jX//6F04//fQ6c5MnT8Z///tf3HDDDQhKoVNdXa1eHG+RQImySFhqwIABB30zR8M555yjfD9dunRRS1eyZPbvf/8b11xzjZqXH5YsZT366KPo3bu3L71clKO8iISQ5pFTWo1XftGqnt99Vn90d26DG3oVydHEjWRbactX8v9OGFADCzra5GuCT+iU2cuwt2KvOu6T0Ketb4eQVkWv16kUcsm6yi2zKU+OLFdJJKekNutK5uW6lmTWrFkoLi5WvTHrR25kVUiiPa0hdJplRj733HOVShNKSkpw/PHHK8Um4y1VQ0d48cUXceGFF+Kvf/0r+vfvr1LIr7/+emWE9vL3v/9d5eZfd911GDFihEpjmz17NmvoEHIUPDMnC9UOF4Z3TcA5Q9LRM9IGvd6AQsSq5SoROSJwJLFcygbuQwfo9UZ1XTDiXbbqFN0JcZbAh8oJCTZG90pWKeQSuamyOZFfYVN7OQ9UarkIGcmKbmh5SoSOpJuvXbsWQVlHJzk5WTmmJcry5ptvKkEi0ZbPP/9cVSnOzNTW9UMFWe6SH4Q4wWNjWzZ0R0iosXFfGSa9+KsqPPzNzSdiSOd4YM8KON+bDIOj3LdUJXV0JMpjgwkViEGC/JV4+cygbHr5wcYP8PTyp3Fqxqn4z6n/aevbIaTNcLs9KrtKjMfiyZHlqpaO5AQbxuZ6Z2JiYtTxjz/+qLKsxEtzwgknqGUsQkjo8tTsTUrknD0kXRM5FfnAvEdgdJSreXHIOZRDR68KBlphRxQKoYsbomVABbE/p1kdywkJI/R6HQZ3bl9RzWYtXfXq1QtfffWVqj8zZ84cny8nPz+fERFCQrx6qvTAMep1uGtiXyB7OfDaWGDHAkA1wdRBpzeqDA2zzgOTHjDo65obgxEKHULaL80SOrI8JX4ZqY4s/hxvvyuJ7khRP0JIaIa0n5ytRT6mjOyCrjs/A949CyjPAeK7ApEdgKhkWe+G3u2AweNQezlHZCJQWaDVtwky7C47tpdsV8eM6BDS/mjW0pUYhMeMGYOcnBxfwUBBigWed955LXl/hJBW4rt1OVi7pxTxZjfudr0KfPuBNtH/HGDIJcDX0wF7hTZmMAEenbj8tD5U1SVae4aGupe3MVtLtsLpcSoTclqUVnuLENJ+aFZER5D6NeLTmTt3rko3FyTrSdo1EEJCC7vTjWd/zEIySvFD3NOIWCsiRweMfwC4+AMgphPgqAI8bsBgrm2maajdm7VxexVgTUSwkVmY6YvmBPsSGyEkSIROYWGhit706dMHZ511lorsCJIrL53ECSGhxX+X7YalKAvfWB9AevlawBoHTPkfcNKdUrBKi9z4qJ+o6Xde57rgYEPhBrVnt3JC2ifNEjq33367arK5e/duX9FAQfpXSA0bQkjoUGlzYvlPn+Fz80PoiP1AYk/g2vlA7wkHLqouBkyRWvdxt1OL4Ig3R/ZyLuMyL9cFGesL1qs9hQ4h7ZNmeXTEdCzZVtLU0x+pTsz0ckJCixWfPY3nXc/AqHPD3fVE6C/5UDMX+yPNM81R2iadyZ01tZEcndbEUyJA3uuCCJvLhi0lW9TxwOSBbX07hJBQETqVlZV1IjleioqK2BWckFDB7YZt9n0Yt/VlpVd2ZUxG1yveAIwNNPaT+jjJvYG8DUBCd8BVUxvJEY+OFajIBVIHBl0dnS3FW+B0O5FgSUDHqOBrTUEICdKlq5NOOsnXAkIQg5802nz66adxyimntOT9EUICgcsBfHUDLMteVqdvWa5A56veaVjkCNJcd8ztWmaViBpRRuZobS/nlhhtvl4T3rZmQ4HmzxmQPIBGZELagKuuukr97nm3pKQknHHGGa3S+sFLs96VnnnmGbz++us488wzYbfbVb+pQYMGYeHChXjqqada/i4JIS2HrQKYeQmw9hM4PXrcYb8BaZP+CYPhCG8HPcYBZz+vRW7slUBFnraX87Of0+aDDBqRCamHlIPY9wew9SdtL+cBRoSNJC3JNm/ePBiNRpx99tkI2qUrh8OBW265Bd9++61KLZcUc2mkKW0gbrrpJqSnpwfmTgkhR48U9fvoImDfKtj1VlxXMx25KWPx7KBG1pcRMdPtJK0woNTMEU+OLFcFWSTHC4UOIX5sXwAseg4o2AK4HYDepC1JSzQ2gH+oWCwWVZJGkP3dd9+tVob279+PDh06IOiEjmRbScgpISEB//znPwNzV4SQlqdsH/Den4DCLXBbE3Bl1R343d0Tr5/Wp2lN/UTUBGHjzvpUO6uxrWSbOmbGFWn3bF8AzLpNi+hGJABGC+C0ab47GZdobStEZSUw8uGHH6pWUrKM1Ro068+wyy+/XLVfJ4SESPi5eCfw9hlK5CC2M17r+Qp+t/fEkM5xOG1AKsKRrKIsuDwuJEckIyUypa1vh5C2w+3WIjkicmLSAVOE1rtO9nIu4zIfoGWsWbNmITo6Wm2yCvTNN9/gk08+Uc3Agzbryul04u2338ZPP/2EYcOGISoqqs78v//975a6P0LI0YafC7YC7/8JKNurMqYKLvgMz7+6VaWH33Fan7A16fovW4Xr90hIo8hdo71fSCSn/u+CnMu4zMt1AYjWSpLSjBkz1HFxcTFeeeUV5fFdtmwZunbtGpxCZ/369TjuuOPU8ebNm+vM8Q2FkCAKP+dtBN4/F6jMB5L7Ald+jZcXFMHmdGNY1wSM6xP49fG2zrjishVp91QVan8UyftFQ8h4TUnAetVJMESWqry8+eabiIuLwxtvvIFHH300OIXO/PnzW/5OCCFNCz97/zKT8LPRqnUZl3kxC0tIOGetFsmRasWpg4Erv8J+dwxmLl2tvuzW8b3D+g8TX0SHhQJJeycySYv8yh9F8n5RHxmX+VYq+CnvO7Js5e2TGZRChxASBOFnabLpdh1ov+AffpZGmxLJEZHTaRhw+edq/s0fMlU0Z2hGPE7qnRy2P8ZKRyV2lO5QxwOSBrT17RDStqT5FfyUP4r8/8CRVi7qj6HAFfy02WzIzc31LV299NJLypR8zjnnoDWg0CEk1MLPshXuq9eGwQpEd9Dm9q0Gfn4UqC4COh4HXPGlatFQXGnHB0t2qYe65dReYR3N2Vi4ER54kBqZqszIhLRr9LUFP2V5WyK//sveInICXPBTemB6S8+IGblfv3747LPPcPLJJ6M1oNAhJFSQsLIsX5Xs1s4lkqMSJ92As1obN8cA8x7W3rzShwJXfOHrQ/X2bztQZXdhYMdYnNovvLOQ1u7Xqq4O6TCkrW+FkOCgR23BT28ig3hyZLlKIjkBrKPz7rvvqq0todAhJFQQn43HqS1XGSx+4We9diz9p+TNS6I8cu0VX2l/uQEorXbg3d92quPpYR7N8Rc6QzsEV+8tQtqUHqFV8LOloNAhJFTIWwfoDFojTSV4dIBHB+g8fvUvPEBCN2U89u9A/t7inSi3OdE3NQanD2hkFeQQxePxYM3+NeqYER1CQrPgZ0sS3jKOkHBC/gKT5aqoZM1AKI053XZt73Fp14gQOvV+7ZpaKmxOtWwl3HRqr6ZVQQ5B9lXuQ2FNIYw6I/on9m/r2yGEtDGM6BASah6dmtpaFwaTsufA4zhwjfhxknrW+bIPf9+FkioHeiRHYdLg8O9F51226pfYD1YxaRNC2jWM6BASih4dMRFK9AbOehfV+nNqqba78Oav29XxX0/pBUOYR3MEGpEJIf5Q6BASih4dSSN32bQlLEGN14ofua6Wmct2o6DCjozECJx7TEe0B7z+HBqRCSEChQ4hoebRie2sZVn5ixwpGBiXoc3XlnGvcbjw2gKte/dfT+4FkyH8f91tLhsyizLVMY3IhBAh/N/5CAknj47OCFTlA27vkpUsRdWKHm+Dz9oy7p+v2oP8chvS46y44LjOaA9kFmbC6XYi0ZqITtGd2vp2CCFBAIUOIaGC1LsQj469UjuXJSwx2xoMWsHAsn1atlXaULjcHryxUPPmXHtSD5iN7eNX3b9+TrjXCiKENA5mXRESKqx8F6jI0451+trKyN7lK7+lLABzNuRiZ2EV4iNNuPT4DIQrbo9bLVWV1JQg3hrP+jmEkIOg0CEkFNgyF/j+Tu1Y2jwIEsURcSMixxihpZZXFsCTsxqvLqhSl1w5qhsizeH5a740ZyneWvcWdpTtUMtVRr0RRdLfi0ZkQoKO3NxcPPbYY/juu++wd+9epKSk4JhjjsFtt92G8ePHB/S5w/MdkJBwImcN8OlUwOPWlqqikuGp2K8FcNT/6aDzeKAzmtWy1pJtBVi7xwGrSY+po7oiXEXOI0seUV3K4yxxMBvMqHJUwS4FFCWt3lHd1rdISEhEQfsn9odeIsQBZOfOnTjxxBMRHx+PZ555BoMHD4bD4cCcOXNw0003YdOmTQF9fgodQoKZkmzgo4sBRyXQ6TigaDfcJdlwezxwQQ83jNB7PDA4qqEv3g19RAJmrNfetC4enoGkaAvC8Y1aIjkiclIiU3xeHFdtdWiDzoAPMz/EmM5jAv4GTkgo0VAUtHtsd0wbPA0j00cG7Hn/+te/qt/TZcuWISoqyjc+cOBAXHPNNQg0fBcgJFipKQVmXgxU5AIpA4A/fw6HywmPxwUH9PBAL7EctVfnHhfW2FLx626bKgwoJuRwRP4alTdqieT4G44loiNEmaLUvDfNnBACXxR0c/FmRBojkRyRrPZyLuMyHwiKioowe/ZsFbnxFzleJMoTaCh0CAlGnHbgkyuA/I1AdBow5TO4i3ehyqlFckxwQ6f1f1B7OZfxV2tOU18urR4yEiMRjkjIXf4aleUqocZZgwp7BSocFeo8xhyj5uU6QggOioJKaxSJdspezmVc5uW6lmbr1q2q0W6/fv3a7EfBpStCgg3x3cy6DdixADBHA1M+BeI6Y1fmKsS4dShHChJQAjMcSuSI1LHBjPXuHpjjGqYe4vpx4RnNEcRXICH3Ulspyuxlqkigxy0CUFu6crgdal6uI4TgkFFQQc5l3BsFHZg0sEVfMhE5bQ2FDiHBxm//AVZ/pFU8vuhdIH2oGi5CDKwwwqU3Yh86wuSxwaA+3g1w6Cz42HkK3NBjaIoeAzvGIVwR82SiJRGbijQDo9FghEfvgUt6gAEoqCpQDT3ZuZwQjfpR0PrIuPzREIgoaO/evZWYCrTh+HBw6YqQYGLT98BPD2nHZzwJ9NaWogRzp2OwW9cJUe5y2F1ulLtMKHFZ1X6/04rvncPVdReN6o92gU77a1RKCbmlq7uWf6YViyaEHBQFtbu0rMT6yHigoqCJiYmYOHEiXn75ZVRW1hY79aOkJPBLzBQ6hAQLeRuAL67VigAOvwY4Xo4PMLBTPL6JuRjlbgs6oAhWnR16nUftf3Ydo5avEo02XHZ8F4QzEl4vshUhNTJVeQzcyp/k8v1lKuMyTzMyIRoS3ZTsKlnurb+UJOcyLvOBioKKyHG5XDj++OPx+eefY8uWLcjMzMQLL7yAUaNGBfzHFPRCRwoLXX755UhKSkJERITKv1+xYkWdH9IDDzyA9PR0NT9hwgT1IhISUlTsB2ZeCtgrgO5jgTOf1goB1mOtaSjuc/4FWZ4uiEINOqAEeo8bH7tOVfPRsfFh3/rAG4YXX0GXmC51elplxGSocZqRCTmAXqdXKeSSkZhfla8M/GI8lr2cy7jMB6ocQ48ePbBq1SqccsopuPPOOzFo0CCcdtppmDdvHmbMmNG+PTrFxcWqyJC8OD/88AM6dOigRExCQoLvmqefflqpwvfeew/du3fH/fffr8JkGzduhNVqbdP7J6RROG3AJ5cDpbuBxB7ARe8BBtNBl23YV4bCCjuKYofhr9WD0NO5HQkox2Z3J1QgQqWUO5xudd3gznHtIgwvER2bvH7yZqY3qoiOvHnTjExIXUamj8QDox7w1dERT478nvRJ6BPwOjqCBCNeeukltbU2QS10nnrqKWRkZOCdd97xjYmY8Y/mPP/887jvvvtw7rnnqrH3338fqamp+Oqrr3DppZc2+Lg2m01tXsrKygL6fRBySCSM/O1tQPbvgCUOuOwTIDKxwUuLquxwuDxIibEgIdKMCsdgFLtcyCupUY+TEm2B0+NR17WHMPyGwg0qclMtrTCkYKDbhZ2lO9Wbt2SO0IxMSF1EzIxIG9HqlZHbmqD+7r755hsMHz4cF110keqLceyxx+KNN97wze/YsUP1z5DlKi9xcXEYOXIklixZcsjHfeKJJ9R13k3EFCFtwuIXgTUztSadF70DdOhzyEsTI80wGXTKiOw13FbaXXC6PSqaE2ExwKTXqevCGXlTHt1xtCoQKCLHU9vYVJbs5FzGZT7c37wJaQ56nV79IXBipxPVvj38ngT1d7h9+3a1fifpadIT48Ybb8Qtt9yilqkEETmCRHD8kXPvXEPcc889KC0t9W3Z2dkB/k4IaYBtPwM/PagdT3wC6HX4xnYDO8aiZ0o08strsGN/JXYWVqCgQove6HVAYYVNzct14Yx4CxbvW6yquloNB5anJeMqwhihxmU+EMXPCCGhR1AvXUnKqER0Hn/8cXUuEZ3169fj1VdfxdSpU5v9uBaLRW2EtBlFO4DPrtYadR5zOTDy+iN+iV6vw9jeyViyrRAutwd6vz9TZEmr0uZS83Jdeyh+lhyZrIoDZpdnq/5WYkQWoSMFBANV/IwQEnoEdURHzEsDBgyoM9a/f3/s3r1bHaelpal9Xl5enWvk3DtHSNBhr9TMx1Kcq9MwYNK/Gsywqo/b7cHCLQWINBvUVls6RkVzotSYUc3Lde2l+Jm3v5VEdiQEL8tXMs6sK0JISAgdybjKysqqM7Z582Z07drVZ0wWQSMpav7G4qVLl7ZKbj4hzTIff30zkLceiOoAXPwBYGpcdqBkU23Lr0BqrBWpsZZaZwrQJSES3ZOjkBJrUfNyXTjj3wKiuKZYjYk3Z3f5buwq26XGmXVFCAkJoXP77bfj999/V0tX0hhs5syZeP3111UXVEH+ervtttvw6KOPKuPyunXrcOWVV6Jjx46YPHlyW98+IQ2bjzd8AeiNwMXvA3EHasAcCW/Wldmg93lzoi0GGA16ZU62GPRwuNtH1pW0gMitzFXFAgWT3qQiOtWOauRV5ql5Zl0RQoLeozNixAh8+eWXyjz8yCOPqAiOpJNPmTLFd83f//53VVb6uuuuU6Wkx4wZo1rCs4YOCWrzsbR36Dq6SV/uzbraX2FDWY1TjVXbXdhVVAmL0YC4CFO7yLqqj2r7IHi0P368WViEEKLeIzzB0Fq0jZHlLkkzlwys2NjwzlghbUTxTuD1k4HqYs18fO5LjfLl+CPem3NfXoR1e8t83hyzUa9Ww5wu6WKuZWZ9fdOYsDYkS/2c2+bfppaoZMlKhI4IHPlPCgjGmmOV2Hn+lOdpRiaEBPfSFSFhYz7+eIomcppgPm4Ih9TQqUVq5/iCF+Graxo0IztcDpVd5Yvo1L4O8neb+HNoRiaEeKHQIaS5SNrTvj+ArT9pe28alD8SbvnmlmaZj+sjJuM9RVoVYG/ARooFuj0eRJgMSIuzqhYR7cGMLGnl3jo5Br1BZVrJXsRPTkWOmgtEJ2ZCSOgR1B4dQoKW7QuARc8BBVsAtwPQm4Dk3sCY24Ee4w5ct/Q1YP3/mmU+rs+ekipVCVnonBChdJVEeEwGPeIijYBHh/wKW9ibkfsm9FV9rrzRHKmhI+ihpZfb3XbVDkKuI4QQRnQIaY7ImXUbkLcBMEcB0anaXs5lXOaF7GXAj//Ujk9/rMnm4/r8vq1IrdAY9ToVuckrr0FhpbbfVViNkmpHuzAjZxVnwenRzNiCRG9kyUr2Mi6RHb1er64jhBAKHUKagoRRJJJjqwBi0gFThNanSvZyLuMyX54PfHYV4HYCA89rVOXjw2FzuvD9+hx1LFWRaxwu6HU6GA06ta+2O5FTWo2kaHPYt4DYX7lfeXAEMR9LirnL41J7OU+PSldRHvHyEEIIl64IaQq5a7TlqogEzVAslXndLkBvAEyR2vj+zcDHlwJle4HkPsCfXmy2+djLl6v2Yn+5TaWXO12SR1074W9Gbif5k3sq9qi9iJmuMV1hc9vUUpVEcqRCsvh0HHoHPTqEEAWFDiFNoapQ8+TIVrgPcNbUKgwdYLQC0R201g7lezXhI74cS8xRvcYSwXl94XZ1LKbjmBgjSqsdsDllyUbTUDIeG2HymZEHd44L25/r1pKtai/ZVeLJkf5WXmQJS9LO+yT0YcFAQoiCQoeQphCZpC1flWj91lQkR60AuwFntVYvx9s1+5wXgJT+R/36zt2Yi+0Flaq/lRiP4yPMSIg0o8bhhtPthlGvh9Wk1dMJdzOy+HAW7lmojqVeTn5VPuIscSrrSgzKInKiTFGYNniaqpRMCCF8JyCkKaQOBsQIK8tVOqPmz5FlI/WhajggcoZdDQy56KhfW4lQvLpAi+ZMGpwOi1EPu8utRTLMBsRYTWov5zbJwApzM/KGgg3YX71fiZlHTnxERW6qnFUoqC5Qezl/YNQDGJk+sq1vlRASJDCiQ0hTyFsHSDqzpIsrwSPF6kTpuDXjsSBzx1zWIq/r79uLsDq7RFVA/tvpfZFbVoPMnHKkxWqp1P6CqKTKgf7pMWFtRp6fPV/tT+x4IsZ0GoPRHUcjsyhTGY+lbo70t2IkhxDiD4UOIU316MhyVVQyUJF/QNz40AFSqM5W3iKv64wF29T+4uGdkRpnxY3jeuLeL9cht8yG+EiTauQpkRwROdLgU+bDuf2DV+icnHGy2ouoGZg0sI3vihASzHDpipDmeHQq92vnBpMW4fGioiw67bqjZMO+UizcvF9VQb7upJ5qbHSvZDx+3mAVuamyOZUnR/ZyLuMyH65kl2crI7JkW43tPLatb4cQEiIwokNIcz06Bou0xQU8Dm1OV+vRkXm57ijxenPOHtIRXZIifeMiZk7okaSyq8R4LJ4cWa4K50iO8Ev2L2p/XOpxyoBMCCGNgUKHkOZ6dCTF3Gs+9vplZFzm5bqOxzb7td1VWInv1u5Tx9eP63HQvIiacE4hP9yy1SkZp7T1rRBCQgguXRHSHI9ObOdaceOt0ifVkSOBuAxtXq47Ct74dTvcHmBcnw4Y2LF9CZqGkKyqlXkr6/hzCCGkMVDoENIUxHsjDTxriv2MyLW+HClk423weRQeHamA/OkKrfrvjSdr3pz2zpydc1QNncHJg5ERk9HWt0MICSEodAhpCmlDtb5W3j5KEr2RisgGg1YwsGyflpEl1zWTd37bAbvTjWO7xGNk90T+fAB8v+N79Tqc1f0svh6EkCZBoUNIU6gs0KofK3S1GVe1y1dH2c9KKK9x4IPfd6ljSRX3r5XTXskuy8ba/WtVKvkZ3c9o69shhIQYFDqENBbJtJJmnWp5ygiYojQzskt6X7kAY20HcxFD0vyzGcxcuhvlNU70SonGhP6p/Nn4RXNGpo1EckT4ps8TQgIDhQ4hjeWXJ4G9miEW0WnaXgVzJMVclScGjGZNCDXDjFzjcOHNRTvU8fVje4R9unhjkIrP3+34Th2f1YPLVoSQpkOhQ0hj2DIXWPi0dmyOASpyAVeN5s2Rejqyl3Np9ikFBZthRv7yj73KiJweZ8W5x3TizwVAVnEWdpTugFlvxvgu4/maEEKaDIUOIUdCxMsX1x5o1mkwNtzUU1c73oyCgS63B6/Vtnv4y0k9VG8rAny3XYvmjMsYhxgRmIQQ0kT4bkrI4XDagE+vBKqLgY7Hac06/Zt6qkrIngMVkf0LBjaBORtysbOwCnERJlw6gunT6qV3O/H9ds2fM6n7JP47JYQ0CwodQg7HnHuBfX8AEQnAxe9pzTolpVwKA4r5WJappJ6O7OW8GQUDxYcy4xctmjN1dDdEWViwXFiQvQD51flItCbipM4n8d8pIaRZ8B2VkEOx9lNg+ZtaGvn5bwDxXWorI5u0Zp6JPbTaOe7aSI4IHWcN4GpawcBFWwuwbm8prCY9rhrdjT+PWj7d/KnaT+41GWaDma8LIaRZMKJDSEPkbQS+vVU7HnsX0Ps07VgKASb31payBGn7YInV9oKMy3wTCga+OG+r2l92fBckRvEDXdhdthuL9y2GDjpc1Oci/hslhDQbCh1C6lNTBnx6BeCoAnqcApx8t99vjB4YcztgiQbKcwBHtebPkb2cW2K0ebmuEfy+vRDLdhbBbNDj+rFs9+Dls82fqf2JnU5E55jO/DdKCGk2FDqE+CPG4q9vAgq3ao07L3hL89z402MccPbzQOpAwF4JVORpezk/+zltvpG8+PMWtb94RGekxVn5swBgc9nw1dav1GtxSd9L+JoQQo4KenQI8WfJy0DmN5oPR8zHUYfw2oiY6XaSVgFZfDviyZHlqkZGcoSVu4rx29ZCGPU63DCO0RwvP+78ESW2EqRFpeGkTjQhE0KODgodQrzsWgzMfUA7PuMJoPPww782Imo6Htvs1++l2mjO+cd1QueEWo9PO0cy0D7e9LE6vrD3hTDUj6YRQkgT4dIVIUJ5HvDZVYDHBQy+CBjxl4C+Luv2lGJ+1n5Il4e/ntyLP4Nalucux9qCtaoS8gV9LuDrQgg5aih0CHE5gf9drXltOvQHzvlPi3Qib4w3R1o9dEuO4s+gljfWvaH25/c+nw08CSEtAoUOIfMeBnb9BpijgUs+AMyBFR6ZOWX4cWOe0lI3ncJojpd1+9fh95zfYdQZcfWgq/nvkhDSIlDokPbNxm+AxS9ox+e+rNXACTAvzdfq5pw1OB29UqID/nyhFs2Z1GMSOkZ3bOvbIYSECRQ6pP1SsBX46q/a8aibgYGTA/6UW/PL8f26HHU8/VRGc7xsLt6M+dnzVYHAaYOnBfznQAhpP1DokPaJvUpr1mkvB7qMAiY81CpP+/L8bapUz+kDUtEvLbZVnjMUeGOtFs05retp6B7Xva1vhxASRlDokPaHKI1ZtwP5G4CoFODCd7TeVQFmR0Elvl69Vx1PPzXwS2ShwobCDZi9c7Y6vm7IdW19O4SQMINCh7Q/Vr4DrP0Y0BmAi94BYtNb5Wn/89NmuD3Aqf1SMLhzXKs8Zyjw/Mrnfd6cvol92/p2CCFhRkgJnSeffBI6nQ633Xabb6ympgY33XQTkpKSEB0djQsuuAB5eXltep8kiMleBnz/d+14woNAtzGt8rRb8srx9Zp96viO0/q0ynOGAgv3LFSZVia9CdOPnd7Wt0MICUNCRugsX74cr732GoYMGVJn/Pbbb8e3336Lzz77DAsWLMC+fftw/vnnt9l9kiCmLAf45ArA7QD6/wkYfUurPfXzP21RK2YTB6ZiUCdGcwS7y46nlj2lji/vfzk6RXdqtZ8HIaT9EBJCp6KiAlOmTMEbb7yBhIQE33hpaSneeust/Pvf/8app56KYcOG4Z133sHixYvx+++/t+k9kyDDadPMxxW5WlHAyTMCXhTQy4Z9pfhuXY56utsZzfHx3ob3sLt8NzpEdMD1Q69vlZ8FIaT9ERJCR5amJk2ahAkTJtQZX7lyJRwOR53xfv36oUuXLliyZMkhH89ms6GsrKzORkIQtxvY9wew9SdtL+eH4od/AHuWAZY44NKPAEvr1a95bq5WBfnsIR3DLtPK7XErM/Fve39TezlvDNtKtmHGmhnq+I7hdyDKxOrQhJB22tTz448/xqpVq9TSVX1yc3NhNpsRHx9fZzw1NVXNHYonnngCDz/8cEDul7QS2xcAi54DCrZoS1HSbVyK/Y25Xess7s/KdzUDMnTAhW8BSa3XKXxNdgl+ysxTPa1umxBemVZLc5birXVvYUfZDjjdThj1RnSP7a7q4IxMH3nIr5Nr/7non3C4HRjbeSwmdZ/UqvdNCGlfBHVEJzs7G7feeis++ugjWK3WFnvce+65Ry17eTd5HhJiImfWbUDeBq1dQ3SqtpdzGZd5f/Pxd3/Tjk+9D+h9Wqve6r/mblb7ycd2Qs8O0WElch5Z8ogq9BdpjFR9qWQv5zIu84finfXvqOhPjDkGD456UCUYEEJIuxQ6sjSVn5+P4447DkajUW1iOH7hhRfUsURu7HY7SkpK6nydZF2lpaUd8nEtFgtiY2PrbCREkOUpieTYKoCYdMAUAej02l7OZVzm5br65uOT7mzhW/GoLuQLNu9Xezn3Z8XOIizcvB9GvQ63jg+faI4sT0kkp9JRiZTIFFiNVuh1erWXcxmX+YaWsUQIvbLmFXV8z/H3qOsJIaTdLl2NHz8e69atqzN29dVXKx/OP/7xD2RkZMBkMmHevHkqrVzIysrC7t27MWrUqDa6axJQctdoy1URCZqZ2FEFuF2A3gCYIrVxmd+7HJhzX8DMx4u3FmDGgm3Yll8Bh8sDk0GHninRuHFcT4zulQyPx4OnZm9S1140vDO6JoWPByWzKFMtV8VZ4g6Kxsi5jMu8XDcwaaBvrtxejjt/uVMtXZ2ccTLO7nF2G9w9IaS9EdRCJyYmBoMGDaozFhUVpWrmeMenTZuGO+64A4mJiSoyM336dCVyTjjhhDa6axJQqgq1CI1shfsAZ42UOtb8N0YrEN1Bm1vwjGY+tra8+VhEzr1frkOFzYmESDPMBj3sLjcyc8rV+OPnDUaV3YXlO4thMepxSxhFc4SSmhIlVswGszqvcdb4PDoS1ZHxMnuZus6LRHfu/fVe7CzbidTIVDw06iEuWRFCWoWgFjqN4bnnnoNer1cRHcmmmjhxIl55RQuNkzAkMklblirZrZ1LJEetwLoBZ7U2brACW+dq4ueCt1vUfCzLUxLJEZGTFmv1fVhb9QakxeqRW2bDK79sRV6ZTY1ffWJ3pMdFIJyIt8YrUVNqK1WCxuay+bSmxWBBrDlWzct1XiTD6pc9v8CsN+M/p/wHSRFJbfo9EELaDyEndH755Zc652JSfvnll9VG2gGpgwGPU1uuMlj8lqP02rGrBnBXaEPjHwB61y1JcLRs2FemlqskktPQsk18pAnr95aipNqpjm88ufUyvFqL/on9kWhJxKYibWnOaDAqj44bblQ7qtXWL7Gfuk6Ys3MOXl3zqjp+cPSDGJh8YDmLEEICTcgJHdLOyVun9ajSG2sFjw7w6ACdp24dnd6na6nmLUxRlV15cmS5ygMPauxuON1uGPV6WM16mPQ6lNU41bU3n9ILcRGBbxYqy0Lih5GlIomiiMAQ4RFw5GWX//PI/zw+sec9FhZkL8DdC+9Wx1P6T8Gfev4p8PdFCCF+UOiQ0PPoyHJVVDJQkQ+4NVFRB6mpM/zqgFQ+Tow0K+NxSbUDpdUO2Jwu1dpBp3PDEpkHuC1wexLQIdqMy0/oimCtZXM0iKgqshUpr039pSvx6MjSlcx/kvUJnln+DJweJ87sdibuGn5XQO6HEEJCNr2ckEN6dCr3a+cGE6AzwSOfslL1WqfDpsg4uKW2TgAY2DEWSdFm5JRWo9ruhF6ngyl6G8yd3oY7+WNUVGt+nJSMJVhTsCJoa9m0hBlZsqu6xHRBRkwGOkZ3VHs5l/EqR5USOVIUcHyX8XjspMdgUH4qQghpXSh0SOh6dCRyozci1yjBBA+qdTrclJKM6QlW3LD+1YB90NdZuonYCkOHz6Gz5MBeNAZwW2Gw7EeF5eeAio2jqWXTUmZkacopS1Vq+Uq9HDp1XlBdoCI9InLGdBqDp8c+rbqTE0JIW0ChQ0Lao7ND50aaw6Gmnk+IR6nBgEiPB5uLMgMiNMSMXFhhR3qcFREmHfQJ8wG9Da7qDDhLj1XXWBJ/R6wpKaBioym1bFoa8QDJ8pgImp2lO1Vjzn0V+7CrbBc2F21W48JZ3c9SGVbeNHRCCGkLKHRIaHp04jLg1pvQxa6lcf83JhqrIyKgN5hhhQ4pxuiACA2vGTk+woz0DsUwWQtg8ETBWXySrKPBFLkbhohdcHkQULFRv5ZNfWRc5v1r2bQUEjka3XG0Wp6qlpR+9UaiV6+z+HGEE9JPwJMnPUmRQwhpcyh0SOh5dGQZxFEJj7Ma4vqYHxGBr6K9lYeVMxg6gykgQsNrRpYCgW69pLG74KjqDkdlD1XLJ6LDIuXbkbYPgRQb/stH3qJ9FfYKtRdkvH4tm5ZCBM3ifYuVHyjCGKHO7R67L9vKahCpWTf7ihBC2gpmXZHQIm0oYI0F8jcqkfOHxYzXE+KUuKmBR0UUjAYLrKYImOE5qEJvS5iRpdXD2j0lcFYZ4Ek1oHr/GDVnilsNj2kfzEYTrCYDbK6agIkN7/KRNMcUMWV3232ZT1KUT55X2i94a9kEZNnMGqdeW2/ETJ63c0xnddxQCwhCCGkLKHRIaGErA4q2q8ManQ6PJCehUqdDmUGvMq403LCU726wQu/RotfrMLZ3MpZsK4TLnQqj6SR4HMmAvhLmDt/BY6yGQ6fHrrLSgIoN7/LR8tzlcHlcKqPJoDPABZdaTpJjmQ9EPZ2cihxVFVmiR96ojbzWSdYkZYYW4dPSApMQQpoLhQ4JHVwO4MMLtP5WOj3MxgiYPR5sN2n/jI0eEQAiczwNVuhtqRYQC7cUINJsgNOlQ1XJCDVu6fAjdIYa6DxSIRgBFxv+y0cSxZKIjggeiejIcpJRZ1TzVw68ssWev7imGB9mfogPN354wJuj06tlKvHrVDmrDtkCghBC2goKHRIaSFW+WbcBe2tr08SkA9WlvmlvLEfV0/F4VK1AKZjc0nhbQKTGWlFYYUOVwwS9OQ+muLW19QlF5uhgMUTArG95sVF/+Sg5MlmJixpXDVxuLbIjHhkp4tcSy0fShX3N/jX43+b/qVYO8jyCPKf4gGRetYDAoVtAEEJIW0KhQ0KDhc8Cf3yoIjkwRwMVecg0m1BkMCDV5UKZXg+bkjtiRhZDrBmx1gRVobclvSLerCuJ7MixPFlE8q8wupOgg2QcueB065AcFQezyRUwr4p/1pWkk0sUx5+GOog3FhEvUnRw7q65+HHXj9hRusM3J+LlL4P/gjfXvomskqwjtoAghJC2hkKHBD8r3wXmP6odn/Ek8MsTqmBgiSESTh2QLOnebhck1uCSJG8pnJfYTUV0pKZLS3pFJOvKqAeyi2XpRgdD1CbAshcOlx46txEGnUmlMpoMepgNhoB5VfyzrsQXI34ZbwsIOW9q1lVeZR6W5S5Tnh/Z763Y65uTCNEZ3c/AhX0uxJDkIdhYtBHF9mJfCwivV0dEj38LCJqRCSHBAIUOCW42fAXMqm3OKU06M473FQyMdzmUL8cuERwtkKOllsvmrIFNr29xr4hkXYkhWdLLlek5eR6gcwEeafJpg1PGDMGddSXiL6soS3Ufl708hhT980ce48ROJ+K0rqfh5IyTEWOOOSiaJFEkif4oap9bLWXpjcrDQzMyISQYoNAhwcv2X4AvrgUkffm4qcD4B4Ft83wFA/tX7kd3hxMbzEY4dTrYlUlGW74yV+6D0Wht8ayn/eU25JRqPhV5Op3HCp1RKgGLEVjiSR44dDrsKitrk6wrMSaLKdgAAzpGdcQ3277BnvI9yC7P9m0ltpIGH0/u8/i04zEibQSGpQ5DpCmywecW4SZmaMm+kkiOPLfXoyPeIBmPtcTSjEwICQoodEhwsncl8PEUQAri9f8TcPZzmrLwFgw0mKBP7IHRujIs10tStdQl1jY5rnY7YHC4Wzzr6R+fr4XbI0tTUhDQCFdlbxgjtgE6N+AxQCcLVx5Po7OuVDVht1OZfGW5SZaBRCx4NzlX4y5tvNxerjZZMpq9Y7Z6DIlgyWPIf17kFfl86+dqq48sMXWN7aoMw30T+6r9kA5D1JJTY+ib0FcZn0VgSeTH2+tKxI7ci8oAc7vUdYQQ0tboPL7Yc/ulrKwMcXFxKC0tRWxs497sG+KLLV/gl+xf6owdZMxs4NVuyLxZf6yhH1Njvq4xzxfIx27W8zmqgdx18LgdgDUO6DBAEzneJ83PBOxV8BhM2Kp3oar2Vtx+lZG9H7qSft0jvkej7vPgb6/uQHmNA9v2SzVkwGLUq1R2uz4f0NkOfgDV7BIqqiPLPiJo5MNfIi7eYxEKgTLuiriSTSJKvRN6+7qMe7dDRWsagyx13fjTjUpwyc9SBJ23ErJ8T/K6y/c8Y8IMFgwkhLQ5jOi0IFuKt2B+9vyWfMj2i9kbn7EB+/9oYF5fG7upjyaIZBlFNES5o1ylR7cUhtrkpgOxk0OhSRjp4F1UU9Sox/aaeSVjStK3xQQsx7K3GC1qTASEbFKw77e9v6ljEVMinOQZ5ZooU5Q6Fi/OjUNvVF6blkS8NyJu0qPS1fcmkSavR0fuP9GaSI8OISRooNBpQSZ2m3hQ9MAb1j/cWP3u04f6uoOuqd+1uqHnasQ1R/qahr6u2c9Vb6jONdUlwC9PQleeA8R2BE75J2CJafh+8jZiU+b/MNNVjFhVN0cHhywlmaJgMEXApDepD3tZ4pnSbwr6JfU7qp/FV3/swTdrc5AQaUKMxYQ9xdXQW/fClLgIHpe1tm2cU2KkMOmNyEiIVc8vgmT6sdMxPHW4irCIn0WK+XmPRTB4a9/IPTf0/IeKqqwrWKfuTTw3/mIj0EX7vBlfcr8SKWqoho9D76BHhxASFFDotCDHpByjNtIMKvKBd88GCrYrozEu+waI0/omNUjX05DefTy+/vlmuN1ulLmqYZPO2c4qwHWgQq9EN07ucvJRLaGs31uKb3+vhMsdj4fOORavLdwGV2UZPK5omOKXQ6eTkswVgE7iPB64dDrsr7ap55cIx8DkgQ0unx0NYhxOtCSqzClBFe3TtU7RPm/Gl9TaSYlMqVPDR5ayRNz1SejDgoGEkKCA3ctJ21NZCLx/LlCQBcR2AqZ+e3iRU0v/5AFIjExFnr0U1S67+qCXSIPs5YNeasOIGDiaD3ub04U7Pl0Np9uDSYPT0TUpCoUVdqTHWWF1Z8DjjAKMpYBeKx6ok+U2j67Fnv+ISOaXKgNd631SFiUZDNxTyus7bfA0JSLzq/KVYVp8R7KXcxmX+UC0viCEkKbCdyLStlQVaSInfyMQnaaJnMTuQfNh//xPW7A5rwLJ0Wb83+RBvsrI8RFmJXrMRun1pEkcX3a776YQMKQYnxTlk6J9EjUSE7DD41B7OZdxb9G+QDAyfSQeGPWAitxIOrv4gWQv5zIu84QQEgxw6Yq0HdXFwAeTgbx1QFQKcNUsIKlnsz7sA1Ghd9XuYry2YJs6fuy8wUiMMqvKyJJaLgUDdZY98BgqoHfHw6Ovkgo62hdK0T6DFQnWuIBVCA6Gon0iZqTmjnx/8jzi3ZHoFSM5hJBggkKHHB1uN5C7Bqgq1GrcpA0F9I0IFNaUAh+cD+SsASKTtUhOcu+g+bCvcbjwt8/WqJo55x3bCRMHpvkqI/dMiUZmTjnirJJq7pIWnrWp7QduQFLPAyk2gqVonzeFnRBCghUKHdJ8ti8AFj0HFGwBpOaNFPITsSKtGnqMO/xy1UcXAvtWARGJwNRvgJR+QfVh/8ycLGzfX4nUWAseOufAB7m0f7hxXE/c++U6FJeb4Y5ww6OX9HFZMhOBp4PRoAu42GDRPkIIaRz06JDmi5xZtwF5GwBzFBCdqu3lXMZlviEq9gPvnaNVPo5IAK78GkgdeNQf9iplu/afs+zlXMabU6H3t60FePs3rWP3kxcMQVykqc786F7JePy8weif3AceaU8h9Xw8UqBQrzw7Rv3RPX9jyCrOgl6vpah7ixBKFEtVWvY4NdGn16vrCCGkPUOhQ5q3XCWRHFsFEJMOmCIAybCRvZzLuMzLdf6U7QPePQvIW1/ryfkeSB8SVB/2hRU23P7JaunigMuO74JT+qY0eJ2InXv+FI9YqwUGvREGvQcmg+Y/bg2x4V+0T/xIEsUSYSV7OZdxmWdjTUJIe4dLV6TpiCdHlqskIlO/wJ2cy7jMy3Udj9XGi3cB7/8JKN4JxHbWlquaYDxujQq9IpLEl5NfbkOvlGg8cPaAw15fai+B2WhAnLX1KwSzaB8hhDQOCh3SdMR4LJ4co0U7d1QBbpfWVVx6KMm4fLjLdcL+LOCD84CyvUBCd03kxHcJug/7t3/biflZ+9Xy00t/PhYRqg1F6z1/U2DRPkIIaRxcuiJNx9tBXNLDC7cBRTuAkl3aXs5lXOblut1LgbdO10ROcl/g6h9aROT4f9hLJV5Bsq+izdG+Sr0yLvONKdi3OrsET/6g1Zy5f1J/9EuLbdXnbyos2kcIIY2DQoc0HUkhj0rWPDcSzZF0chE2spdzGZf50n1aMUCJ7nQeAVwzG4hND7oP+6JKO/764UpVCPCMgWm4/ISurfr8zYVF+wgh5MjoPL4CJO2XsrIyxMXFobS0FLGxR/5Lvt0jJuM3TgFy12meHFmyUjZcj7aEJf+kpCmnCB6PC+g9EbjoHS0rKwAszVmKt9a9hR1lO1RdHVlOkkiKiIwjVeh1uT2Y+vYyLNpagO7JUfjm5hMRYzW12vO3BCKuWLSPEEIahkKHQqfp7PsD+PjyA4X/nNWauBHRY7BqkR1bmTZ/zBTgnP8AhqaJh9b6sH92ThZemr8VESYDvrrpRPRNi2nV5yeEEBJYaEYmzTcjS4SmtreUFs2RXQ0gXcSFgecB5758cGZWkFTo/W5tjhI5wpMXDG62yGnu8xNCCAk8FDqk6YjJWJavSnZr5wbp2G0EVHp17UqoORo48daAihy324MN+8pUo03pQSXtGaRycWNYt6cUd362Wh1PG9Md5x7TKWD3SQghpO2g0CFNJ3WwFrURP45BUsw9msjRQjsaslQl1wWIxVsLMGPBNmzLr1AmYmm0KT2opD2DFPM7HPllNbj2/RWocbgxrk8H3HtWy2dFEUIICQ5oIiBNR7qN6wyA3gi47YDTT+SocZO2l+sCJHKk11RmThmiLEakxFjUXhptyrjMH4oqu1OJnNyyGlUU8MU/HwtDI6NAhBBCQg8KHdI8j44YbZVHx6/Ng662YGBchpaJ5S0Y2MLLVRLJqbA5kRZrhdUkbRZ0ap8Wa0GFzaXm5br6OF1u3DzzD6zZU4r4SBPevHI4YpuYYUUIISS0oNAhTUcqH0tWla+1ge7AJh4dbydz8fK0MOLJkeWqhEgzdPX8P3IuAkbm5Tp/pIrCfV+tx8+b8mEx6vHW1BHolhyYdHdCCCHBQ9ALnSeeeAIjRoxATEwMUlJSMHnyZGRl1W2SWFNTg5tuuglJSUmIjo7GBRdcgLy8vDa757Amdz3w7W2As6Z2QAcYzYDRqpmSJdXcWzBQCgu2MGI8Fk+O2aCHBx5U210or3GovZxbDHo43B51nT/P/bQFHy/PhqxSvXjZsRjWNaHF740QQkjwEfRCZ8GCBUrE/P7775g7dy4cDgdOP/10VFZW+q65/fbb8e233+Kzzz5T1+/btw/nn39+m9532CGRmhXvAG+OBwq3ahEb+ecjS1SyfKUKBboDnkou2VViPC6pdmBnQRV2FVViT3G12su5jJv0OnWdlzcWbscL87ao4/+bPAinD0wL6D0SQggJHkKuYOD+/ftVZEcEzdixY1U14w4dOmDmzJm48MIL1TWbNm1C//79sWTJEpxwwglHfExWRj4CNWXAt7cCG77QzjNGAsW7tcyrqiKt+rG/TycyUcvGuvTDA93LWwjx3pz78iK1NCWSymjQK20l/4rFgyP/mCXN/Oubxijvzge/78L9X61XX3vnaX0wfXzvFr0fQgghwU3QR3TqI8JGSExMVPuVK1eqKM+ECRN81/Tr1w9dunRRQqchbDabEjf+GzkEe1cBr43VRI5kWZ32f8DYv2lLV9UlB1LJ9eYD1Y9l3FEZEDNyHbzBI69UrxdM+mxFtk/k/PXknrj51F6BvR9CCCFBR0jV0XG73bjttttw4oknYtCgQWosNzcXZrMZ8fHxda5NTU1Vc4fy/Tz88MMIC6RwX+4aTVSI+Vd8MdKC4WhxOYBf/w0sfBpwO4G4LsCFbwMZI4A9q7TmnbJUZTDXFRhy7LID9irAqonRlkQiOYUVdqTHWVFa7YDN6fZ1n5A2DrERJjX/77mb8fIvWtXjq0Z3w10T+x5kXiaEEBL+hJTQEa/O+vXrsWjRoqN6nHvuuQd33HGH71wiOhkZGQg5ti8AFj0HFGw5kOmU3BsYczvQY1zzH1ce74vrgH2rtPMBk4Fzngciag28Ov/VTjn2FxB+c3Wua1kzstTOkcwrKfrndLth1OthNemV6NlVVOVr7XDlqK548JwBFDmEENJOCRmhc/PNN2PWrFlYuHAhOnfu7BtPS0uD3W5HSUlJnaiOZF3JXENYLBa1hTQicmbdBtgqNAEiKd9SuC9vgzZ+9vNNFzsSHVr+BjD3QS17yhoHTPo3MOiCuibj6mKtXo4sT0m0p373cjmXebkuQGZku8sNi6lu5EqyrnLLbKrGjnDtSd1V1WNGcgghpP0S9EJHvNLTp0/Hl19+iV9++QXdu3evMz9s2DCYTCbMmzdPpZULkn6+e/dujBo1CmGJCBKJ5IjIiUk/IEJMEVqad3mONt/tpMYvY+VvAr69Bcheqp33PFVryBnb8eBrZYlMigXKprqX1xyI7BgjNIHkva6FEaOxtHpYu6cETpdHCR7NTq89v7O2UOBNJ/fE37hcRQgh7R5jKCxXSUbV119/rWrpeH03cXFxiIiIUPtp06appSgxKMfGxiphJCKnMRlXIYl4cmR5SSI59X0nyqySoM3LdUfKepIokIiiX/+leWukGeeEh4ARfzl0qrj4gGSJTKJHCd0BV01tZMcIGKxARS6QOjAgdXQkk2ps72Qs2VYIl9sDo0EHvc4Dh1sTxcL4fh1w1xn9Wvy5CSGEhB5BL3RmzJih9ieffHKd8XfeeQdXXXWVOn7uueeg1+tVREcyqiZOnIhXXnkFYYsYj8WTI8tVghiD/ZeMZFyqFh8p62nXYmDW7cD+Tdp5nzOASf8C4g4sDTaIRInEByRLZCJqRFiJQBLRJOeWGG2+JUzRDaSXL9xSgEizQQkdm9OlRI4XMSTbXR51XWM7mRNCCAlfQq6OTiAIuTo6+/4APr5ci7hIKnedpSMrEBGvFZY5VB2bshxg7v3Aus+088hk4KyngYHnN63gX6DM0Idh3Z5SXP/BCtXEU/7p7iysUstVomm6JEaqujpVNideu2I4BneuXUIjhBDSbgn6iA5pAFkSkhYLubXdwQ3yY5ToiVuL7siWNvjgpSOnHVg6A1jwNGCv0ITRsKnA+Ae1In9NRcSM+IACkd5+hKyrGrsLe0urIZYcaQfRNSlSNfaUSE5pAy0gCCGEtE8odEIdXwSmNjDnLRPsj5xnfgP89DBQtE0b6zwCOOuZo65c7IYOG9w9UOTujES3GQOhC2gVyhiLEVV2J/ZXuHxLVcnRZvUtStaVzeU+qAUEIYSQ9guFTigiEZTKAi3jSmU9VWtiRueX9STzcp2jRlum2rNc+9qoDsBpjwBDLj3qyMvirQWYsWCb6hYuURZJ+5aMqBvH9cToXsloaXYVVuKBb9aj0q6JHBE0DpcLOaU16luXyI6Yk4d0jlfZWYQQQgiFTiibkSW9WwtlaBEdT200R7qJi3fn+7sOCBwxKY+erm1iFj5KROTc++U6VbNGCveJyJBU78yccjX++HmDW1TszFq7D/d8vg7lNqeK4ogJ2enLuoJawqqyu2CozcqiEZkQQohAoROKiBdGaumU7NbODVKwT/pMubXoTtEOTfjsKdKabB53JXDy3UBMy3TtFh+MRHJE5KTFWn0F+ax6A9Ji9apon8yf0CPpqAVHQYUND3y9Ht+v08oKDOuaAIMO2JRbrrKuRFy5apumSyaWQa9XWVnTxvSg2CGEEEKhE5KkDtY6h0tKuXQJV74ct9afSvZeJItKBE6Hvi3eb0qWqySSU7/qsJzHR5rUvFzX3Mwnyaj6Zs0+PPTNBhRXOVSkRooAntovFX/9aCVSY60wG3UorXLCIb4cgx5xkUbYnJ6jfm5CCCHhAyM6oUjeOi1SI5vbrgI5dQSOZFOJT+fEW1pc5PhnPslylRiAa+x+/abMelgM+qPKfNpbUo0Hv96AnzLz1Hn/9Fg8c+EQDOoUhwWb96vnlkhOTqldLWF57Ukl1QYkRZvhYNYVIYSQWih0QpHyXK2KsZiJXVpfJw2d5sWJTtX6UB2pYOBR9psqqXbUdhA/IDYsRgPiIkzNynyqtrvw6oJtapOu5PIc00/tjRvG9YTZqBmn5THdHg/2Flerc4n06PSaNanG4VLj8vzMuiKEECJQ6IQSkkm1/C2tFo6tzG+idvlIIjxy7C3eF4BeU4JkNEnkRJaH5JmlSJ9XbFTbnSr9W65pbOaTLFN9uzYHT3yfqTKohOO7J+KRcweiX1rdx+ifFgOXx6P8ObJ0pa9dOlM7vQd2p0fNy3WEEEIIhU6wI6bjHb8Aq94HMmdpIkYQISP9pVQ1ZNOBgoFiRi6rbrhgYCCoV8bH28S8sQJHjMPPzd2M1dklaqxTfAT+Oak/zhyU1mDX8czccmVGlmwrpxsw6j2+p1TnBp2al+vo0SGEEEKhE6yU7QP++Aj44/0D2VVCx+OAkTcAS14G8tYfMCKrtaPa9aMAd/WQSE5hhR3pcdbapSutg7jqJ2oyIDbCpOYPZQgWgfPb1kI899NmrNxVrMbk6248uSeuG9tDVTg+FOL70ev06BgXgcJKW73n1iMpyoIqh4uVkQkhhCgodIIJWzmw6Ttg3f+AbfMOGIwtccCQi7Q08fShWq8r1XIhEagqqo3s4MDylYx7CwYeZeXjw5mRU2IsKvOqxuFnRjbplfDIr7AdJDZkuWleZh7e+HU7lu/UBI7FqMflJ3TF9eN6ICXG2mh/kHh2uiVFHfTcNeLtcbnp0SGEEKKg0GlrpHLxlh+B9f8DNs+pbdBZS5fRmrgZcC5gjjwwLiLHXlnbr0ocuSatWKBEdGSpS4oFWqIDbkaWzCeLqV51ZR00I7GfGVnq7fxvRTbeWbwTuwqr1JgIlT8f3wV/PbknUmKPLHC8iO9Hqi9LYcK0WAsizBL9MfgiRSVVDvRPj2FlZEIIIQoKnbagpgzY+hOQ9QOweXZdY3FSL2DQhcDgC7VO4A1hTdQad0rEx2A+4JMR5FgysuxV2nUBwCs21u4pgbM21du7fOTfhkFu5cGv1+OLP/aivEaLOklG1GXHd8FVo7shLa7xAseLFCCUFhNSfVkKE0rNHklnlx5XInKiLQY1z8rIhBBCBAqd1qJ0L5D1vbbt+PWAqViI7QwMOh8YdIG2NNWACbcOErnxIcf+13sOcV3Loa9ts7BkW6FajlIGYL3WhkH6UEkx5M155Tj7pUW+r+meHIVrTuyGC4Z1RqT56P7ZSWsJaTHh7bMlNXskgiSRnED12SKEEBKaUOgEClla2rUY2DYf2P4LkL+h7nxSb6DvmUC/SUDn45vWYLO6WKuXI7VyxJ+jr00rF5Ej1ZLlXOblugAgLSAkW0paLojQkTo6jlo/tJr3AHllNrW8dfqANFw0vDPG9u7QolEWETPSYkIMz+IFkmUyiTQxkkMIIcQfCp2WJHed5rMRYZO9VFtC8qEDMkYC/c4C+p516GWpxiD1caShp2yqe3nNgciOt3u597oAIOIiK7ccFpNBFfmTtG5/xH9jNerxypRhGNM7cNEVETVMISeEEHI4KHRaksUvAWs/PnAe1wXoeTLQ4xSg+zggqoWEh9THEaGUtwFI6A64amojO0bAYAUqcoHUgS1aR0eKAEoq+KItBZi1Ngf55bY685I9Jf4b2cSnI1lXUriPEEIIaUsodFoSWYqSFPGepwA9TwUSexzZb9McZJlrzO3ArNs0URORAJijAadNO7fEaPNNWQ6rR1mNAyt3FuP3HYVYur0I6/eWwilrUn5IOrfUzLEaDOqpJMVbsrAk5bs5LSAIIYSQlkbnkZzcdk5ZWRni4uJQWlqK2NjGtS0ICrYvABY9BxRsOdD2QSI9InJ6jGv0w8jy08acMiVm1u4pxbq9JdiaX6G8Nv5I1eKRPRJxcp8OmLlsNzJzyg6bdfXe1cfTM0MIIaRNYUQnlBEx0+0krTCgKiCYpC1XHSKSI8bh3UVVKiNqi2z5FcprI3uZq0/XpEiM7J6Ikd2TlMDpnHCglk9eWQ1W7Cw+KOuqyu5SjTYlK4vGYEIIIW0NhU6oI6LGr/qxZEDtLaxQgia7qErtZZNCfdsLKmGv7xyuJTnagiGd4zCoUxyGyNY57pCF/OpnXUlEx+XWIjoyZtDr1fy0MT0odgghhLQpFDohgqwwltU4kV9Wo1K3JaKSV16D/NpjMQfnlFQjp6zmsK2uxFfTKyUafVJi0Cs1Gr1TYjCoUyzSYq0NNtE8VNaV1K9JjbWqDuKlVU44XG6YDHrERRphc3rU/KF6XRFCCCGtBYVOKyMRlUqbU7VFqLQ7UVHjRHGVA8VVdpTW7uW8tNqO4sra8WoHiirtqrVCY5CoSpfESGQkRqq9d+vZIRqdEyKOOsri7XUlkZycUrkvl8+jU1JtQFK0GQ63h401CSGEtDkUOi3I16v34tctBQeEjNpcStB4j0UcHA2Svp0aa1HRFGmCmSLHMdp5apxVCZqkKHOjozPNQbKp3B4P9hZXq3Px5Oj0WsHAGodLjct9MuuKEEJIW0Oh04Kszi7B/1buadS1Uncm2mJEtNWI+EgzEiJNiI8w1R6bVQ8n2bzHsu8QY4HVpDWwbEv6p8WoGjniz5GlK32tqFI7vQd2p0fNy3WEEEJIW0Kh04KM75eqoizSWDLKYtQ2s+wNStQcGDPAKGlKIUpmbjkMUoTZoFNVkY16j7cBhXYuWVg67Tp6dAghhLQlFDotiLQ7CGTLg2BBPDp6nR4d4yJQWGlT3iGvRyfCpEdSlAVVDhc9OoQQQtocCh3SZMR7Iw07padVt6QoVQnZ6XarysiS1VXjdMPkctOjQwghpM0J3fUT0mZIl/CeKdEqO0yIMBsQYzWpvVBS5VDzch0hhBDSllDokKb/o9HrcOO4nsqLlFtmQ7XDpYoIyl7OZVzmWRmZEEJIW0OhQ5rF6F7JePy8weifHoMqm1N1K5e9nMu4zBNCCCFtDZt6hnJTzyBAIjlSAVkMyuLdkeUqRnIIIYQECzQjk6NCRA1TyAkhhAQrXLoihBBCSNhCoUMIIYSQsIVChxBCCCFhC4UOIYQQQsIWCh1CCCGEhC1hI3RefvlldOvWDVarFSNHjsSyZcva+pYIIYQQ0saEhdD55JNPcMcdd+DBBx/EqlWrMHToUEycOBH5+fltfWuEEEIIaUPComCgRHBGjBiBl156SZ273W5kZGRg+vTpuPvuuw+63mazqc2/YKBcz4KBhBBCSHgR8hEdu92OlStXYsKECb4xvV6vzpcsWdLg1zzxxBOqErJ3E5FDCCGEkPAj5IVOQUEBXC4XUlNT64zLeW5uboNfc88996jojXfLzs5upbslhBBCSGvSLltAWCwWtRFCCCEkvAn5iE5ycjIMBgPy8vLqjMt5Wlpam90XIYQQQtqekI/omM1mDBs2DPPmzcPkyZN9ZmQ5v/nmmxv1GF4/tpiSCSGEkEASExMDnU7HF7mVCHmhI0hq+dSpUzF8+HAcf/zxeP7551FZWYmrr766UV9fXl6u9jQlE0IICTTM8G1dwkLoXHLJJdi/fz8eeOABZUA+5phjMHv27IMMyoeiY8eOypDcVJXtTUuXr42NjUV7ha8DXwf+W+DvBN8bGv8eKZ81pPUIizo6bfmPV9LT27s65+vA14H/Fvg7wfcGvkcGKyFvRiaEEEIIORQUOoQQQggJWyh0jgKpxSP9tdp7TR6+Dnwd+G+BvxN8b+B7ZLBCjw4hhBBCwhZGdAghhBAStlDoEEIIISRsodAhhBBCSNhCoUMIIYSQsIVC5yh58sknVTXl2267De2Jhx56SH3f/lu/fv3QHtm7dy8uv/xyJCUlISIiAoMHD8aKFSvQnujWrdtB/x5ku+mmm9CecLlcuP/++9G9e3f1b6Fnz574v//7P18/vfaEtNaR98WuXbuq12L06NFYvnw5wpmFCxfinHPOUdX25d//V199VWde/h1IBf/09HT1mkyYMAFbtmxps/ttL1DoHAXyS/vaa69hyJAhaI8MHDgQOTk5vm3RokVobxQXF+PEE0+EyWTCDz/8gI0bN+Jf//oXEhIS0N5+F/z/LcydO1eNX3TRRWhPPPXUU5gxYwZeeuklZGZmqvOnn34aL774Itobf/nLX9S/gw8++ADr1q3D6aefrj7Y5Q+DcEV6LA4dOhQvv/xyg/Pyb+GFF17Aq6++iqVLlyIqKgoTJ05ETU1Nq99ru0JaQJCmU15e7undu7dn7ty5nnHjxnluvfXWdvUyPvjgg56hQ4d62jv/+Mc/PGPGjGnr2wg65PehZ8+eHrfb7WlPTJo0yXPNNdfUGTv//PM9U6ZM8bQnqqqqPAaDwTNr1qw648cdd5znn//8p6c9IB+vX375pe9cfhfS0tI8zzzzjG+spKTEY7FYPP/973/b6C7bB4zoNBMJyU+aNEn9hdJekZCrhGh79OiBKVOmYPfu3WhvfPPNNxg+fLiKXKSkpODYY4/FG2+8gfaM3W7Hhx9+iGuuuaZJTXLDAVmemTdvHjZv3qzO16xZoyKdZ555JtoTTqdTLeNZrdY647Jc0x4jv8KOHTtU02n/zwzplThy5EgsWbKkTe8t3AmL7uWtzccff4xVq1aF/Xrz4ZBfznfffRd9+/ZVSxUPP/wwTjrpJKxfv75ddebdvn27Wqq44447cO+996p/E7fccgvMZjOmTp2K9oj4EkpKSnDVVVehvXH33XerJrfiVzMYDOrD/rHHHlN/CLQn5D1g1KhRyp/Uv39/pKam4r///a/6QO/VqxfaIyJyBHkt/JFz7xwJDBQ6TSQ7Oxu33nqrWnuu/9dKe8L/L1TxKInwEdPhp59+imnTpqG94Ha7VUTn8ccfV+cS0RGxJ2vw7VXovPXWW+rfh0T72hvy7/+jjz7CzJkzlYdt9erVypArr0V7+/cg3hyJ6nXq1EmJvuOOOw6XXXYZVq5c2da3RtoZXLpqIvJLmp+fr35pjUaj2hYsWKAMZnIsf8G1R+Lj49GnTx9s3boV7QnJnhgwYECdMfkLtj0u4wm7du3CTz/9pIyo7ZG77rpLRXUuvfRSlX13xRVX4Pbbb8cTTzyB9oZknMl7Y0VFhfoDcdmyZXA4HGqpuz2Slpam9nl5eXXG5dw7RwIDhU4TGT9+vMogkL/UvJv8RS+haTmWv1zaI/Jmtm3bNvXB356QjKusrKw6Y+LPkOhWe+Sdd95RXiXxr7VHqqqqoNfXfVuV9wSJ/LVXJLNI3hckQ3HOnDk499xz0R6RkgMiaMTD5UWWOSX7Spb5SODg0lUz1p4HDRp00C+y1FCpPx7O/O1vf1P1IuQDfd++faqLu7yhS2i6PSF/rYsBVZauLr74YvVX6+uvv6629oZ8mIvQkSUaiW62R+R3Qjw5Xbp0UUtXf/zxB/7973+rJZz2hogaST4SH59EeiXaJd6lq6++GuH8B59/VFsMyPIHcGJiovo3IcuYjz76KHr37q2Ej9RckmXNyZMnt+l9hz1tnfYVDrTH9PJLLrnEk56e7jGbzZ5OnTqp861bt3raI99++61n0KBBKk20X79+ntdff93THpkzZ45Kqc3KyvK0V8rKytR7QZcuXTxWq9XTo0cPlU5ts9k87Y1PPvlEff/yHiFp1TfddJNKpw5n5s+fr34H6m9Tp071pZjff//9ntTUVPV+MX78+Hb9+9Ja6OT/2lpsEUIIIYQEAnp0CCGEEBK2UOgQQgghJGyh0CGEEEJI2EKhQwghhJCwhUKHEEIIIWELhQ4hhBBCwhYKHUIIIYSELRQ6hBBCCAlbKHQICTNOPvlkVWqeEEIIhQ4h5Ai8++670Ol0apN+ZgkJCRg5ciQeeeQRlJaW8vUjhAQ1jOgQQo5IbGwscnJysGfPHixevBjXXXcd3n//fRxzzDGqqSshhAQrFDqEhDHFxcW48sorVRQmMjISZ555JrZs2VLnmjfeeAMZGRlq/rzzzlPdtuPj4+tcI9GctLQ0pKeno3///pg2bZoSPNKt+e9//7vvOpvNhltuuQUpKSmwWq0YM2YMli9fXud+pkyZgg4dOiAiIkJ1cZaO516ys7NVF3h5fun4fO6552Lnzp0BfY0IIeENhQ4hYcxVV12FFStW4JtvvsGSJUsgPXzPOussOBwONf/bb7/hhhtuwK233orVq1fjtNNOw2OPPdaoxxYxI6JFHtvlcqkxET2ff/453nvvPaxatQq9evXCxIkTUVRUpObvv/9+bNy4ET/88AMyMzMxY8YMJCcnqzm5J7k2JiYGv/76q7q36OhonHHGGbDb7QF7jQghYU6r9UknhLQK48aN89x6662ezZs3e+RX/LfffvPNFRQUeCIiIjyffvqpOr/kkks8kyZNqvP1U6ZM8cTFxfnO33nnnTrn/syYMUM9R15enqeiosJjMpk8H330kW/ebrd7Onbs6Hn66afV+TnnnOO5+uqrG3ysDz74wNO3b1+P2+32jdlsNnW/c+bMafbrQQhp3zCiQ0iYIhETo9GojMNekpKS0LdvXzUnZGVl4fjjj6/zdfXPD4dEiLxLW9u2bVNRmRNPPNE3bzKZ1ON5n+/GG2/Exx9/rLw9Ev2R5S8va9aswdatW1VERyI5ssnyVU1NjXpsQghpDsZmfRUhhNSKKTEqi4ASs/KREI/Qrl278P3332Pu3LkYP348brrpJjz77LPK7zNs2DB89NFHB32deHoIIaQ5MKJDSJgipmGn04mlS5f6xgoLC1UUZ8CAAepcojv+ZmGh/vmhyM/Px8yZMzF58mTo9Xr07NkTZrNZeWu8SIRHHs/7fF7RMnXqVHz44Yd4/vnn8frrr6vx4447Thmlxfsj3h7/LS4u7qhfD0JI+4RCh5AwRTKaJGvp2muvxaJFi9TS0OWXX45OnTqpcWH69OkquiKZViIyXnvtNWUUlqWo+ktUubm5KmojUZy3334bo0ePVgLkySefVNdERUWppam77roLs2fPVqZjee6qqiqVpSU88MAD+Prrr9US1YYNGzBr1iwlyAQxNosxWe5NzMg7duzAL7/8orK4JK2dEEKaA4UOIWGMpG7LctDZZ5+NUaNGKcEiwka8M4L4aV599VUldIYOHaoEyu23365Sw/0pKytTqeUikuRxRBBJVOaPP/5Q415E9FxwwQW44oorVIRGBM2cOXNUersgEZ977rkHQ4YMwdixY1UBQvHsCJLevnDhQnTp0gXnn3++L41dPDqyPEYIIc1BJ47kZn0lISQskSjMpk2bVFSFEEJCHZqRCWnniBFY6ufI0pMsW0kNnFdeeaWtb4sQQloERnQIaedIJWLxwpSXl6NHjx7KtyNFBAkhJByg0CGEEEJI2EIzMiGEEELCFgodQgghhIQtFDqEEEIICVsodAghhBAStlDoEEIIISRsodAhhBBCSNhCoUMIIYSQsIVChxBCCCEIV/4fg6LEc5nIHzMAAAAASUVORK5CYII=",
"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": "dr_curve_python",
"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.13.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
@jfhc
Copy link
Copy Markdown

jfhc commented Aug 6, 2015

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?

@yannabraham
Copy link
Copy Markdown
Author

@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.

@akashbahai
Copy link
Copy Markdown

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'

@cbmII
Copy link
Copy Markdown

cbmII commented Dec 13, 2020

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])

@yannabraham
Copy link
Copy Markdown
Author

Thanks! Fixed the code as you suggested, should work now @akashbahai @cbmII

@Fae14
Copy link
Copy Markdown

Fae14 commented Mar 15, 2021

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?

@yannabraham
Copy link
Copy Markdown
Author

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!

@amnfar
Copy link
Copy Markdown

amnfar commented Sep 29, 2022

Nice attempt. May be code should be more explanatory and have more comments,

@resminer
Copy link
Copy Markdown

resminer commented May 14, 2024

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'

@meeuosh
Copy link
Copy Markdown

meeuosh commented Dec 10, 2024

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?

@yannabraham
Copy link
Copy Markdown
Author

@resminer @meeuosh thanks for your patience, the last commit should solve the issue (mixing positional and named arguments is never a good idea!)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment