Skip to content

Instantly share code, notes, and snippets.

@jessegrabowski

jessegrabowski/nonlinear_regression.ipynb

Created May 7, 2024 04:50
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save jessegrabowski/969df4dc1daf3683d177bc1cef0b999a to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save jessegrabowski/969df4dc1daf3683d177bc1cef0b999a to your computer and use it in GitHub Desktop.
Download ZIP
Non-Linear Regression by hand in python
Display the source blob
Display the rendered blob
Raw
nonlinear_regression.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"id": "35df5f55",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from scipy import optimize, stats"
]
},
{
"cell_type": "markdown",
"id": "84fe9624",
"metadata": {},
"source": [
"# Basic way -- no gradients"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "14c030e4",
"metadata": {},
"outputs": [],
"source": [
"def expected_value(params, t):\n",
" δ0, δ1, δ2 = params\n",
" return δ0 / (1 + np.exp(-(δ1 + δ2 * t)))\n",
"\n",
"def logp_data(params, t, y):\n",
" *mean_params, sigma = params\n",
" mu = expected_value(mean_params, t)\n",
" \n",
" # IMPORTANT -- return the NEGATIVE logp! Scipy only knows how to minimize!\n",
" return -stats.norm(loc=mu, scale=sigma).logpdf(y).sum()\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "718f50c8",
"metadata": {},
"outputs": [],
"source": [
"t = np.linspace(0, 1, 100)\n",
"true_params = [3, -5, 10]\n",
"true_sigma = 0.3\n",
"true_mu = expected_value(true_params, t)\n",
"data = stats.norm(loc=true_mu, scale=true_sigma).rvs()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "4c8c95c7",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1520e78d0>]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3fElEQVR4nO3debxbdZ0//tdJcpPcJbn72nu7t7cbpaWl0ELZioUWUWdwxMFhUXGmI4jaYfhZnO+4ff3WGRmmMioMCFQEhdEC4rBI1S5AW6ClhUJX2tu771tyl+zn98fJ52TPTXKz3r6ej8d92Jub5J4bS/O+7+0jybIsg4iIiChHaDJ9AURERETxYPBCREREOYXBCxEREeUUBi9ERESUUxi8EBERUU5h8EJEREQ5hcELERER5RQGL0RERJRTdJm+gGTzeDzo6OiAyWSCJEmZvhwiIiKKgSzLsFqtqKurg0YTPbcy5YKXjo4ONDQ0ZPoyiIiIKAGtra2or6+Pep8pF7yYTCYAyg9vNpszfDVEREQUC4vFgoaGBvV9PJopF7yIUpHZbGbwQkRElGNiaflgwy4RERHlFAYvRERElFMYvBAREVFOSWnw8vDDD2Pp0qVq/8nq1avx6quvRrz/7t27IUlSyMeJEydSeZlERESUQ1LasFtfX48f/ehHmDt3LgDgl7/8JT796U/j8OHDWLx4ccTHnTx5MqDZtrKyMpWXSURERDkkpcHLjTfeGPD5D3/4Qzz88MM4cOBA1OClqqoKJSUlqbw0IiIiylFp63lxu9149tlnMTo6itWrV0e97/Lly1FbW4t169Zh165dUe9rt9thsVgCPoiIiGjqSnnwcvToURQVFcFgMGDTpk144YUXsGjRorD3ra2txaOPPoodO3bg+eefR2NjI9atW4e9e/dGfP6tW7eiuLhY/eB2XSIioqlNkmVZTuU3cDgcaGlpwdDQEHbs2IFf/OIX2LNnT8QAJtiNN94ISZLw0ksvhf263W6H3W5XPxcb+oaHh7mkjoiIKEdYLBYUFxfH9P6d8g27er1ebdhduXIl3n33XfzkJz/Bf//3f8f0+EsvvRRPP/10xK8bDAYYDIakXCsRERFlv7TveZFlOSBTMpHDhw+jtrY2hVdEREREuSSlmZf7778fGzZsQENDA6xWK5599lns3r0br732GgBgy5YtaG9vx1NPPQUA2LZtG2bOnInFixfD4XDg6aefxo4dO7Bjx45UXiYRERHlkJQGL93d3bj11lvR2dmJ4uJiLF26FK+99ho+8YlPAAA6OzvR0tKi3t/hcODee+9Fe3s78vPzsXjxYrz88svYuHFjKi+TiIgIAPDn492w2lz4zPJpmb4UiiLlDbvpFk/DDxERkeBweXDBd/8Ih9uDd+6/FpUm9lOmUzzv3zzbiIiICEDr4BjsLg9kGWjuH8305VAUDF6IiIgAnOvzBSxtg+MZvBKaCIMXIiIiAOf6x9Q/tw2ORbknZRqDFyIiIjDzkksYvBAREQE418/gJVcweCEiIkJg8NLKslFWY/BCRETnPYfLg3a/bEvH0Djcnim1SWRKYfBCRETnvZaBMXhkoECvhU4jwemW0WO1ZfqyKAIGL0REdN4TzbozywtRW2IEwL6XbMbghYiIznui32VmRQEaSgsAcFw6m6X0bCMiIqJcoAYv5YXoM9gBAG0DzLxkK2ZeiIgo6ewuN976uA82pzvTlxKTc31KlmVmRSHqvZkXThxlLwYvRESUdNvfOocv/OJtPP5mU6YvJSYi8zKrohD1pfkA2POSzVg2IiKipPu4ZwQA8GH7cIavZGJ2lxsdQ0qgMqO8ALJ3QprBS/Zi5oWIiOKy70wfvvjkO2gdiFxW6RtR+kaa+7O/9NLqHZMu1GtRWWRQMy/c9ZK9GLwQEVFcfv12C3ad7MX/ftAZ8T69avAyClnO7gDAv99FkiRUm43I00pweWR0W7jrJRsxeCEiorhYbS4AiLrErc/qAACMOtzoH3Wk5boS5T9pBABajYS6Eva9ZDMGL0REFJdRuwhe7GG/7vHIatkIULIv2aypz7fjRRClo2ilMcocBi9ERBSXEW/w0msJH7wMjzvh8usVyfa+l+DMCwDUl4hFdcy8ZCMGL0REFJcRe/SykX/WBQDOZXvw4u15mVXhF7yo49LZfe3nKwYvREQUF1E26o1QNgq+PZvLRjanGx3DYkzaL3gpY89LNmPwQkREcRm1K1tzRx1uNZDx1zsSHLxkb/aidWAMsgwUGXSoKNKrt6vnGw1l77Wfzxi8EBFRzOwuNxxuj/p5uKZdkXmZU6lkMrI58yJKWjMrCiBJknq7OCKgc8gGl9/PS9mBwQsREcVMZF2EcKWjvhFlNHrljDIAwOCYE8PjztRfXALOeSeN/EtGAFBlMqi7XrqCdr1sef4ort+2F0Nj2T0CPpUxeCEiopgFl4nCNe2Kht0ZFQVqKaYlS0tHTeJMo6DgRaORMC3MrpemvlH85p0WnOiyRl3SR6nF4IWIiGI2Ehy8hBmXFtmYiiKDmtFoHsjO0tE5dcdLYcjXROnIP3j5n4Ot6p9f+7BrUt9bluWQySyKDYMXIiKKWXDmJbg5F/BlXiqLDJhRrgQA2dq0K65rlt+COqGhLHBc2uX24HeH2tSv7z/bj8FJbA/+zz+dxsr/+yfsPdWb8HMAgNXmxMO7z+Bs78iknieXMHghIqKYxZN5qTQZMKMse5t2I41JC8GZl10ne9FrtaOiSI/51UVwe2TsPN6d8Pc/3DIIADjUPJjwcwDA7w614d9eO4G/fnhfTpzinQwMXoiIKGbBDbvBPS8ej6yeZVRRZFBX7mfjoroW75i0yaBDeaE+5OvBi+qee7cFAHDTRfW44YI6AJMrHfV7G5s7hia3S+aMN+MyNObE3z52QA2KpjIGL0REFLMRuzI1pNcqbx/B00ZD4064vUcDlBfpMb1MCV6ysWHXv9/Ff0xa8J1vNI6uYRv+cqIHAPC5ixuw4YIaAMCbp/tgtSU2SdU/qrx2IvsTzgdtQ7j18bfxUUfkjErLgPJ4k0EHq82Fv/vF23inaSCha8oVDF6IiChmI97Mi+hlCQ5exOelBXnI02rU84K6LDbYnIFZm0xTzzQK06wL+MpGXRYbnnu3FR4ZWDWzDHMqizCvqgizKwvhcHvUoCYesiyrmZfOocincz/7biveON2H595tjXifNu/hkf958zKsnl2OUYcbtz/xDt76uC/u68oVDF6IiChmomFXvOH3jzrg9FviJpp1K4oMAICSgjyYjDoASpnG36HmQSz97h/xzNvNKb/ucMT1TPc25garLDJAr9PA7ZHx+JtnAShZFwCQJAkbl9QCAF49Gn/pyDLuUg+vbB8ahyzLYe8nTrWOVHbzeGS1J6exxoQnv3gxrpxfiXGnG1/+5bvon6LTTAxeiIgoZiJ4aSgtgFajlFpEBgEIbNYFlDd5kX0RZRrhibeaYLG58Ofj8WcukkFkPOpKwgcvGo2Eeu/XLDYXTAYdNnrLRQBw/RLlz7tP9WDMEXpMQjR9o76gwu7yYCDC1JIITCI1PHdbbXC4PdBpJNQWG2HM0+LR21ZgWkk+bE4PTnVPzQmklAYvDz/8MJYuXQqz2Qyz2YzVq1fj1VdfjfqYPXv2YMWKFTAajZg9ezYeeeSRVF4iERHFQUwbmfN9ZwH5N+0GZ14AYHqYcekxhwt/8QYtmcoOdAx7g5fi8MELAEwr9X3tU8vqUKDXqZ8vrjOjoUwJEvacDBx37rHasO9MX8SMin/ABwAdYUpHHo+Mdm/w0jY4HpDhElq9/S51JfnQefuQDDot6kqMABAxKMp1KQ1e6uvr8aMf/QgHDx7EwYMHcc011+DTn/40Pvroo7D3b2pqwsaNG7F27VocPnwY999/P+655x7s2LEjlZdJREQxEsFLkUGHKpPyBuk/Li32vojMCwDMFMGL36K6v5zowbi3B6ZvJDNvsJ3eRtla7xt9OKLvBQA+f/H0gK9JkoQNonTkN3X0+kdduPY/9uCWx97Gu+fCT/4EB2ztYSaOeqx29Rwpt195yJ8oKzUElb7KvNNTA6MsG8XtxhtvxMaNGzF//nzMnz8fP/zhD1FUVIQDBw6Evf8jjzyC6dOnY9u2bVi4cCHuvPNOfOlLX8IDDzyQysskIqIYibJRoUGnBij+i+r8t+sK6pZdv8zL/77vW63fP2qPmKEAkJKDEccdbgyNKVNCtVEyLyIoWFhrxpJp5pCvi9LRn493w2pz4nt/+Ah//6tDsNiU1+l0jzXs8/YFZUQ6w0wctQ4G9rmcC1M68vXtBC7ZK/e+/v0pyLy09I9F/f8rHdLW8+J2u/Hss89idHQUq1evDnuf/fv3Y/369QG3XXfddTh48CCczuw81IuI6Hwy4he8VHmDF//Mi8iiiJISAMwoCywbjdhd2HXS1+dic3ow5gg/iXSubxTLvr8TW185HvW62ofGQ7b/RiPGk4sMOpiNuoj3+5sVDfjk0lr88K+WhB2nXlZfgtpiI0Ydbqz/z7148q1zAIBqc+hr4y848xJu10tbUPDS3BcavIgAxz9DBEDdWxNcnpqsMYcL1/7nHly69c9hD+VMl5QHL0ePHkVRUREMBgM2bdqEF154AYsWLQp7366uLlRXVwfcVl1dDZfLhb6+8CNfdrsdFosl4IOIiFJDLKkrMmh9wYtfz0twwy7gm0xqH1L6Nv58vBt2lwezKgphzFPehiK9yb57bgAjdlfUTbYdQ+O46se7cO2De9AU5g0+HNGsW1tsDBuUCJUmA356y0W4aHpp2K9rNBKuW6xkXzqHbSgpyMMvbluJv12llJjCHVwJ+H5eMYkVrudF9LMI4SaOfGWjwODFVzZKbvCy/0w/HC4PdBpNQICabikPXhobG3HkyBEcOHAA//iP/4jbb78dx44di3j/4L9EIjUV6S/X1q1bUVxcrH40NDQk7+KJiCiAWjbS+5WNrP6Zl9CyUZXJAGOeMnLcPjiunsb8yaW1KC9U7tcXoTdDlKTaBsbh8YQvVXzQNgynW0bnsA03//d+fByhVOOvQ+13iVwyitXfrKyHXqvByhmlePmetbh2UTWqzUofTXekzIv3571gWjGA8D0vIjCpK1aeK1zZSAQ4wWUjEbz0J7nnRey0uXpBZdSgL9VSHrzo9XrMnTsXK1euxNatW3HhhRfiJz/5Sdj71tTUoKsrcF6+p6cHOp0O5eXlYR+zZcsWDA8Pqx+trZEX+RAR0eSoDbtGHSpFw643eHF7ZPU3/Sq/zIskSeoZRx92DKuTOTcsrVV/e4+UeRGBkcPtQXeELEaLXyNwj9WOzz96ACe7ogcw6ph0ceRm3VgtrivG+99Zj99uWo1p3mBIlI26LeGvWZTXLqhXgpdwPS+iQffyeRUAQg+3tDnd6mvSUBoYhIngMZmZF1mWsdv7/901C6qS9ryJSPueF1mWYbeHjwRXr16NnTt3Btz2+uuvY+XKlcjLywv7GIPBoI5iiw8iolx1qHkQX3nqYMQ3vUwLmDYyB2ZeBscccHtkSJLvN39BjEv/4o0mONwezK0qQmO1yddYGmFcuscvqxPpiAHxpn7LJdOxqNaMvhEHPv/o/qiHFIoek2jNuvHI12sDMhFiEiti5mUkMPPSY7XD4QpsTBb9LJfPq1Q+HxgLaF5WltsBBXptyOudirLR6Z4RtA+NQ6/TYPXsiqQ9byJSGrzcf//9eOONN3Du3DkcPXoU3/72t7F792584QtfAKBkTW677Tb1/ps2bUJzczM2b96M48eP44knnsDjjz+Oe++9N5WXSUSUNR7440nsPNaN599rz/SlhPB4ZLWx1r9ht9eqTAuJklFpgV7dOSKIcekjrUMAgBsuqIUkSb7MS4Q3Wf+SVPCG3uDblzWU4NdfuQRL64sxOObELY8dCFmMJ3TEMCY9GaJs1D9qD7ufRfy886tN0Os0kOXALI3L7UGndw/NxTNLYdBp4PLIAeWlVr9Jo+ASTrlf8BKp3BavXd6S0erZ5cjXa5PynIlKafDS3d2NW2+9FY2NjVi3bh3efvttvPbaa/jEJz4BAOjs7ERLS4t6/1mzZuGVV17B7t27sWzZMvzgBz/AQw89hJtuuimVl0lElBVsTjcOeU8Ezsa17qN+W2SLDDq1NOFwezA87vQ16/r1uwjTywPPD/rkUmU/isi89EX4efv8gpfWCMGLyLzMKCtASYEeT995CRbXmWGxufCH9zvCPqYzhgV1k1FeqIdWI0GWQ382p9ujjmlXFBnUUpN/YNI5bIPbI0Ov06DaZFTPkvJv2m31lpWCJ40AoNQbvHhk5bDMZBATYlc3Vibl+SYj8nxYEjz++ONRv759+/aQ26688kq89957KboiIqLsdfDcoFo6SMV+jskSk0ZajQSDTgNJklCcn4fhcSd6rHZfs64pdApFZF4AoLHahHnVJgATj/T2TJB5cbk96pu+KE2ZjXm4YWktPuqw4Exv6Hp8WZbROZTazItGI6HKZEDnsA09FntAeUqUcjQSUJKfh7oSI5r6RgPGpdUR6JJ8aDQSZpQX4lT3iPeYAF8ZCQht1gWAPK0GZqMOFpsLA6P2kLJSvCw2Jw56F+5dneF+F4BnGxERZY03/U4Bzsbgxb/fRZQp/He99FmVaw6XeRENu4DSqCtUqMvUQjMvYw6X+j2B8MFLx1BghkKYU1kEADjTG1o2sthcGPWWv1KVeQGAKnXiKLB/SQR5ZYUGaDSSGtiIbBDga9YVxxOI4O9cn1/mJcJ2XUF9bZOw6+Wt031weWTMrixUlw5mEoMXIqIsse+ML3jJxrXu/sGL4Nuya1PHmivCBC91JUaYDDpIkq9kBADlUaaNRDAktAyETuSIIwemlxVAo/H1ffiCl5GQbbBisqekIC+lvRvV3temO2iZW3/QIr+6MGWjtqD9LSJg8B+XFtmZhjBlIyC5Tbu+klHmsy5AistGREQUm+ExJ476TccMZOi8n2h8RwP43vADMy+hC+oEnVaDJ754MUbsLsz2BhYAfHtewvy8vSNKJqKkIA9DY070jdgx5nAFHI7o3+/ib0Z5AXQaCWMON7ostoCyjW9BXeqyLgDUaayeoMyLyDKJwG2at3TlXzZqU/tZROYlNHgR01fTy6MHL8FHEcTL45GxyzsinS3BCzMvRERZYP/Zfsiyb+Nq/6gj4+fHBPM/GkAQpZEeqz1q5gUALp5ZFvLmJ7IPA6P2kKkYsVp/dkUhivOVdRnBW2fVs32C3sDztBr1tjM9gaUjMWmUjB0v0VSbwpeNROZFBG4i8xKu50VkVWZWKP/bOjAGt0fG8JhTPT+pvjR8ECaCo8kGwsc6Lei12lGg1+LiWeE3DacbgxcioiwgSkbXe1fN210etS8jW4yGKxsV+cal1UMZw2ReIok2FSOCoSqTUW1KDe57aYmQeQECS0f+1MxLipp1hUhbdkWWSQQXas+L3xEBIkgTgUltcT70Wg2cbhkdQ+NqcFNRpA/IRPkTwdFkS5BiRPryuRUw6DI7Ii0weCEiygJveZt1r11UrZ73k22lI/+jAQS1NGK1qW/K4Rp2I8nTalBSoGRVgsfD/c9JihS8NEfIvACRgxd1x0uaykahmZfADFWdN4iy2l2w2Jywu/w253p/bq1GUhtzz/WPRjzTyJ/viIDJ/T1S+12yYMpIYPBCRJRhXcM2nOkdhUYCLp1drv7GnOxzaSZrRBzKaAxt2O0atqm/4YcblY5GjEsH9734By/iTdp/14ssy2jpFw27oRMwcyqV2yJlXqYl4VyjaETmJfj0ZRFMiJJZgV6HUm8A1zE0jo4hG2QZyM/Tqq8N4N/3MjZhsy4QvRk6VgOjDhz2Lha8Kgv2uwgMXoiIMkxkXS6YVozi/LyUnQg8WSN2pazjXzYSDbstA2PwyMruEhF8xao8wri02PFSFSHz0j/qwKjDDUkKPy48p8qbeQnqeelUMy/pKRv1jzoCVv+LzIv/6+Tf99ImdryU5gdszhWnczf3jaqvQ6QxaSA500ZvnO6FLAMLakwpz1TFg8ELEVGGveXtd1kzVzkvJlnp/mQTS+r8p43E4Yyi17bMu1k2HpEOZ5yobCQmjWrNxrC9GKJs1GWxqc3GsiyjQ2zXTXHmpbQgD3la5bXo9SuJBfe8AL4SVvuQTe13CS4JzfTbshvpNGl/k/17tP9MP/79tZMAsqtkBHBUmogoo2RZxr6P+wEAl81RgpfyrM28hE4bmY06GHQa2L2ZhUiTRtGoZbIoPS8l+cpr0jowBo9HhkYj+TbMRhgVLs7PQ6XJgF6rHWd7R7C0vkTNgkiSLzOSKpIkocpkRPvQOLotNkwryYcsy2qGyf+1EuPSnX4TR8FTRP67XjzeSbRoZSPx/INjDvU1i8W4w41/e+0Etu87p17HbatnxPTYdGHmhYgog872jaLLYoNep8HKmcoYqq9XIbt6XsJNG0mSFLDXJdyOl4mIn9d/H4nH4zvosdJkQG2JEVqNBLvLo2YxROYlWvYhuO9F9LtUFBmg16X+LTB418uYww2bUwn0/DMv/mUjcWZRcGAiel5a+sfUPTDRGnZLC5Tnd3tkWGyxnW90qHkQGx96Qw1cbrlkOl77xhVZVTICGLwQEWXUPm+/y4rppTDmKaWPMrVhN0szL0GjuVV+AUtCmZei0MzL4JgDLm8tqqLIgDytRp3KEaUjsV032rp6deLI2/eSrh0vgm/Xi/KzidJYfp42YMTZF7zY1IxScOalrsSIPK0Eh9sDh8sDrUaK2rej12nUvUHhlgAG+7B9GJ9/dD+a+kZRYzbil19ahf/3VxcEBKvZgsELEVEGvSVKRnPL1duyvWzkP20EKHtYhEQyLxVhDmfsVc//0SNPq7xVqX0v3oxLS0yZl8BxafVAxjRlEqqDxqX7grbrCiIwax8aj5hV0Wk1AdmYuhIjdNrob+Px/F369z+ehNMtY+28Cvzxm1fgyvnZM10UjMELEVGGuD0y9p8VwUuFenu2ThuFKxsBgQFLRdCbcix800Z+wYvod/HL5AQ37YodLzMi9LwAfhNHIngZTs+COsF/AzHgt103KEMlMi9dFptaLgu3Odf/Z40WtAni+0y0qG7fmT7sPdULnUbCDz9zgbrROFsxeCEiypDjnRYMjzthMuhwwbRi9fayJOznSAXftFGyy0Ziz4vvDVYcDSB6RgAE7HoZd7jVAGdGmB0vguh5Odc3Bpfb45s0SlvmJfCIAHVBXWFgkFdlUnp63N5SmcmgCxtA+JfIojXrCrFMHMmyjH/zThXdcsn0iA3Q2YTBCxFRhrzfNgQAWDa9JCD9X5G1S+pE5iVwLNk/wEisbKQ8xmpzwe5SAiRRNoqUeRHZF7NRh+KCyFmCuuJ8GPM0cLg9aBsc95WN0pR5qVYbdpWfRwRowWUjrUZCjd/0U31ZQcCOF2GmX2ARrVlXUMtGUQLhP37Ujfdbh5Cfp8Xd18yd8DmzAYMXIqIM+aBVOUV6aX1xwO0i82JzejDmcKX9usKRZdnvVOloZaP4gxdzvg467xivKJX1hjmh2j94ae6fuFkXADQaCbMrfKUjtWyU7syLd91/X4SyERC48TfSYYtiUR0QW/AyUebF5fbggdeVrMuda2cF9C9lMwYvREQZIjIvS+tLAm4v1GvVMd5sKR3ZXR51+ie452WyDbuSJIWsso8WvPRY7TjVbVVui6HEIfpeTnWPoMsiFtSld9poaMwJm9OtBhHlhaG9Qf7ZoEgloZkBZaOJA7CJgpfn32vHxz0jKC3Iw1eumD3h82ULBi9ERBkw7nDjdI/SRHphUPAiSVLWTRyJkhEQOipdV5IPjaSM/4rdIvESi+pEWaXHm6nwD16K8/PU0d83vSPm4U6TDib6Xg6c7YfbI0OrkdKWYTDn69RAtNdqDzmU0V9dDJmXaaX5KDLokKeVMKsietbJ//uEa9i1Od34zz+dAgDcdfVcmI3Z3aTrL/uGt4mIzgMfdQzD7ZFRZTKgJsyujrJCPTqHbVkTvIiSUYFeG7KptaxQj//624tQaNDGfTSAEEvmRZIkTC8rwEcdFrzXPAQg+qSRIMal325SJrtqzMaErzNekiSh2mxA64CyZbc/zNEAgn/wEqkklKfV4Kkvr4LN4UZJDIFiWZgxdOHpA83oHLahttiIv7s0uzboToTBCxFRBrzfJvpdSsJ+PdvONwp3NIC/G5bWTur5K4IOZ+z1O5TRnwheHG5lS20sfR8ieBGbbVN9IGOwapMRrQPj6LHa1Z8v3OGV0/zKRpEyLwBw0fTSmL93tLH73x5sAwDcs26euiAxV7BsRESUAR94+10uDGrWFSrCbJ3NJDEmnaptq+V+GQKb0w2LTQmWKoPKO8G7TSZq2AWAWRWF8B/cqU3xgYzBRNNux9C4GkSE24fj30QcS1AWC5HhGRh1QPaehwQAYw4XTvcofUPXZNmhi7Fg8EJElAEfiMxLQ0nYr2fbojrfpFFqfkMX0zd9Iw4166LXaWAO2ubr/6au12oCxosjyddrAyZ50nU0gCBGyU91W9XTt0vDNOzOrizE7MpCrJlTnrQgUfw9cnlkWMZ9fUvHOizwyMood6oPqEwFlo2IiNJseNyJpj5l1HfptPCZl2wrG1kjbNdNFrXnZdQesOMleNeJf+alviw/5t6VOZVF6tr9tJeNvMHBsU4LAKCkIE898sCfQafFzm9eiWS24xh0WpgMOljtLvSP2tWdOCJ4vmBaSfK+WRox80JElGZHvW8c08sKwv4GDmTf+UaRjgZIlgq/hl2138Uc2hfiH7zEMmkkiL4XIBNlI2/mpUuZLou2C0erkcIup5uMsqLQv0tH28PvGMoVDF6IiNLMt98l8htHtmVeIi2oSxbRwNo/YlfPAaqMME4sMhOxnO0jzKny9cak62gAQex6EU3G4Xa8pJL4u+R/srToubogQuYv2zF4ISJKs/dbhwCE7nfx52u0TF7D7rm+Uew705fQYyeaNpos9XyjUQd6LaE7XgS9TqM2tk6PoVlXCMy8ZKbnRUhkC/FkBGfxrDYnznrLlksYvBARUSzUZt0omRdfJiJ5mZdNTx/CLY+9jQ+9JYN4pLpsJH5eh8uDpn7l3KJI23rF6xZPyaOx2gS9ToOKIn3aMx9VQQ2x4Xa8pJJ4bUUg/FGHBbKsNC4nshE5G7Bhl4gojXosNnRZbNBI0X/rFX0KYw43bE73pPdwWG1OnOhSRmPf/Lgv7t+4R8SJ0vrUvG3k67Uo1Gsx6nDjhLexNdIW3H//7FL841VzIu7ICae0UI/fbVqNAr026T0lEzEZdMjP02LcqbyG4Xa8pJJ6Srk38yJ6ri7I0X4XgJkXIqK0Esvp5lYVRS3BmLwr4IHk9L2c9AYuAHDw3EDcj1dPlDam7ndeMS4tShqRsgImY15cgYuwtL4Ec6tMCV9fosSWXSH9mZfAstEH7dEXJOYCBi9ERCng8ch45u1mnOiyBNzuW05XEvXxkiT5dr0koXR0vNN3He+eG4THI0e5dyhf2Sh1m1jFm7rbe225WtIIx790FG5BXSoFHxFwNMebdQEGL0REKbH/bD++/cKH+OzD+/FRh6/H5P0JltP5KysMXJk/Gcf9Mi/D4071UMhYpbphFwgtpwQfDZDL/BfBlae7YVc9esGB4TEnznl7ihi8EBFRgM5hZWJmxO7CHU++i9aBMciyPOGxAP4qgg4rnAyRedF554zfjbN0lOpRaSA0I5Hu8koqVfsFYuluGPaVjez40BtIN5TlR9wxlAsYvBARpcDQmC/g6LXacdsT7+D9tmEMjTmh12qwoMY84XMk64gAj0dWe16uX1IDIPHgJVXTRkBgsFJSkAeDLrcOC4wmk5kX/79H6qRbjm7WFVIavGzduhUXX3wxTCYTqqqq8JnPfAYnT56M+pjdu3dDkqSQjxMnTqTyUomIkmp43AkAuG5xNaaV5KOpbxS3Pf42AGBhrTK2O5FkLaprGRjDmMMNvU6Dz61sAAC82xRf8DKSjuDFr2wUbkFdLhO7XvK0Ush5Takm/h453bK65yeXJ42AFAcve/bswV133YUDBw5g586dcLlcWL9+PUZHRyd87MmTJ9HZ2al+zJs3L5WXSkSUVCJ4mV9twlNfXoXSgjz1pORYpzz80/2TIUpGjdUmrJxZCp1GQsewDe1D4zE/R1qCF7/MS7ijAXKZyLxUhDmvKdWMecoYOgC8fVYJWiOdqZUrUhr+vfbaawGfP/nkk6iqqsKhQ4dwxRVXRH1sVVUVSkpKUnh1RESpMzSmBC/F+XmYU1mEJ+64GLc89jbGnW4si6FZF/A17E62bCSCl4W1JhTodVg8rRjvtw7h3aYBTFs+bcLHu9we2JzKavvU9rxM3czLihml+MyyOlwyuzwj37+8yIDRgTH1iILFOR68pLXnZXhYqbWVlZVNeN/ly5ejtrYW69atw65duyLez263w2KxBHwQEWXa0LgveAGA5dNL8fSdl+Cea+bikxfWxvQcySobHetU+l1En83FM0oBAO/E2Pcy6nCrfy5Mw6g0MLXGpAEgT6vBts8vx9+ump6R71/m15w7s7xA/XuZq9IWvMiyjM2bN+Pyyy/HkiVLIt6vtrYWjz76KHbs2IHnn38ejY2NWLduHfbu3Rv2/lu3bkVxcbH60dDQkKofgYgoZqJsVFLge9NYMaMUm9c3xtyImqxpI7FrZmGtN3iZpfwCGeuyOtGsm6eVUtpE69/zEmm7LiXGf8LpghxeTiekrWvo7rvvxgcffIA333wz6v0aGxvR2Niofr569Wq0trbigQceCFtq2rJlCzZv3qx+brFYGMAQUcYNe6eNSgoS/w03GdNGFpsTbYNKb8vCWmW77Epv5uVU9wgGRx0TjsymY0waAEoL8iBJgCxPvcxLpvlnXnK93wVIU+bla1/7Gl566SXs2rUL9fX1cT/+0ksvxenTp8N+zWAwwGw2B3wQEWVacNkoESITMWJ3we5yT3Dv8E54S0a1xUY1C1ReZMCcSuVE5oPNgxM+hzUNzboAoNNqUOq9RgYvyeU/np3rk0ZAioMXWZZx99134/nnn8df/vIXzJo1K6HnOXz4MGprY6sRExFlmscjwyLKRpMIXsz5OnWpXKLZl+CSkbAqjtJROna8CNcvqUF9aX7cB0dSdKJsJEnA4rrc/yU/pX8T77rrLvz617/G73//e5hMJnR1dQEAiouLkZ+fD0Ap+7S3t+Opp54CAGzbtg0zZ87E4sWL4XA48PTTT2PHjh3YsWNHKi+ViChprHYXxNFB5kkEL5IkobRQj16rHf0jDtQW58f9HP6TRv4unlmG37zTGlPTbrrKRgDw//7qAsiynPZx4qlOlI1mVxTCZMztZl0gxcHLww8/DAC46qqrAm5/8skncccddwAAOjs70dLSon7N4XDg3nvvRXt7O/Lz87F48WK8/PLL2LhxYyovlYgoaYa9Y9LGPA2MeZNrcC0XwUuCmRcxaRScebl4ppJ5Odo2jHGHG/n6yNc5YldKVukIXgAwcEmBtfMrcNH0EnVJYa5L6d9EWZ741NLt27cHfH7ffffhvvvuS9EVERGlnjpplD/5s2PE+HAii+rcHhmnugLHpIX60nzUmI3osthwpHUIq+dE3j+SjhOlKbWqTEY8/9XLMn0ZScOzjYiIkmxofPKTRoJ6snQC49LN/aMYd7phzNNgVkVhwNckScLKmcrU0UTnHKVjuy5RPBi8EBElmdiuO5l+F6F8EuPSx70lo8ZqE7Sa0FKMaNp9+kAzdp3sifg8I2nseSGKBYMXIqIkG07CpJEwmV0vYtIo0gnWNy6tw/SyAvRY7fjik+/iG88eDvt90jltRBQLBi9EREnm266bvOAlkYbdSJNGQmmhHq9+fS2+dNksSBLw4pEOXPvgHvz+SHvA/Zh5oWzD4IWIKMmGvNt1k3F+jO+IgPgbdo9HmDTyV2jQ4V9vXITn/3ENGqtNGBh14OvPHsGfjnWr90nnqDRRLBi8EBElWbhzjRKV6MnSw+NOtA8pxwIsiBK8CMunl+IPX7scf7NC2YL+s90fqxOjo95RaU4bUbZg8EJElGSiYTcZmZdEy0YnvCWjaSX5MV+HXqfBP1/fCL1Og8MtQ+rRAb7jAXJ/uRlNDQxeiIiSLBnnGgnVZgMkCbDaXGgbHIv5cSe7xX6X8P0ukVSZjLjpomkAgP/ecwaAf9mImRfKDgxeiIiSzJLEhl2TMQ+rvNtwXz3aFfPjzvSMAADmVhfF/T3vXDsbkgT86XgPPu6xctqIsg6DFyKiJBNlo2Rs2AWAG5YqB9O+fLQz5sec7RsFoJxlE685lUX4xMJqAMBje5s4bURZh8ELEVGSiQ27ySgbAcpJy5IEHGkdUptwJ3K21xu8VMafeQGAf7hyNgDghcPt3LBLWYfBCxFREtmcbticHgBAcRLKRoDSh3KxWjqaOPtic7rRMawEOYlkXgBgxYwyrJxRCofbA3FMHTMvlC0YvBARJZHod9FIgCmJb/Ybl9QAAF6JIXg51z8KWQbMRp06rZSIv79itvpnSQIKJnlCNlGyMHghIkoiMWlkzs+DJsx5QonacEEtJAl4r2UIHROUjvxLRpKU+DVcu7AasyuVzE2hXpfUn4doMhi8EBElUTLPNfJXbTZi5QzlFOhXP4w+dXS2V5k0EoFHojQaCX+/Vsm+JKt/hygZGLwQESWRuqAuCdt1g228QJk6mqjvZTKTRsH++qJ6bLpyDu7fuHDSz0WULAxeiIiSKJnnGgW73tv3crB5EF3Dtoj3m+ykkT+9ToNvbVigjmsTZQMGL0RESZSqshEA1BbnY4VaOgqffZFlOWllI6JsxeCFiCiJhpO4XTccUTqKNHXUP+qAxeaCJAEzyxm80NTE4IWIKImSeShjOBsv8JWOui2hpaMmb79LXXE+jBxtpimKwQsRURINJ/FQxnBqi/Nx0fQSyDLwWpipI5aM6HzA4IWIaAI//uMJ3PXr99Q1+dEMqWWj5E8bCaJxd9fJnpCviWbdOUlo1iXKVgxeiIiiGLW78LNdZ/DyB52465n34HR7ot5/OIXTRsIV8ysBAAfO9sPucgd8TYxJz0rCmDRRtmLwQkQUxcluq/rnPad68e0XjkIWh/2EkeqGXQBorDah0mSAzenBoXODAV9j2YjOBwxeiIiiON5pAQDUFhuhkYD/OdiGh/78ccT7D6VwVFqQJAlr51UAAPae7lNvd7k9aBkYA5CcHS9E2YrBCxFRFCc6lczLpy6sw/c/vQQA8J9/OoXfHmwNua/HI6e8YVe4Yp5SOnrjdK96W9vgOJxuGcY8DWrNxpR+f6JMYvBCRBTFiS4l87Kw1oy/u3QG/vGqOQCALc8fxb4zfQH3tdpdEBUlc4qDl8vmKpmXjzos6BuxAwDO9iklo5nlhTxEkaY0Bi9ERBHIsqxmXhbUmgAA/7y+EZ+6sA4uj4zH9p4NuP+wd8dLfp425TtWKk0GLKo1AwDe+lgJojhpROcLBi9ERBG0DY7DanchTythdoUSEGg0Em5fMxMAcLzTGnD/ofHUTxr5Wzvf2/dySglezvRy0ojODwxeiIgiONGlBCdzq0zQ63z/XDbWKFmYLosNg6MO9fZ0TBr58+97kWUZTX2cNKLzA4MXIqIITngnjRZ6gxWhyKDD9LICAMBxb08MkPqjAYKtmFEKY54GPVY7TnWPJPU0aaJsxuCFiCgCEZiIfhd/C723nfArHQ2ladJIMOZpccmscgDKQY09VqVxl2UjmuoYvBARRSACk4Xexlh/C2qU28QeGACwpLlsBEDd9/Kbd1oAABVF+rQFT0SZktLgZevWrbj44othMplQVVWFz3zmMzh58uSEj9uzZw9WrFgBo9GI2bNn45FHHknlZRIRhRh3uNHUr5RhRKDiT828dPllXrxHA6TyXKNg4qgAkXURjcVEU1lKg5c9e/bgrrvuwoEDB7Bz5064XC6sX78eo6OjER/T1NSEjRs3Yu3atTh8+DDuv/9+3HPPPdixY0cqL5WIKMCpbitkWclkVJoMIV8X2ZiT3Va4vOcdpbvnBQDmVRWh2uy7PpaM6HygS+WTv/baawGfP/nkk6iqqsKhQ4dwxRVXhH3MI488gunTp2Pbtm0AgIULF+LgwYN44IEHcNNNN6XycomIVKIcFC7rAgANpQUo0Gsx5nDjXP8o5laZ0rZd159yVEAlfneoDQAnjej8kNael+HhYQBAWVlZxPvs378f69evD7jtuuuuw8GDB+F0OkPub7fbYbFYAj6IiCZLlIMWhmnWBZR9L2JkWux7GcpAzwvg63sBOGlE54e0BS+yLGPz5s24/PLLsWTJkoj36+rqQnV1dcBt1dXVcLlc6OvrC7n/1q1bUVxcrH40NDQk/dqJ6PwzUebF/2vivsMZKBsBwOVz/YMXZl5o6ktp2cjf3XffjQ8++ABvvvnmhPeVpMAzOcTx88G3A8CWLVuwefNm9XOLxcIAhogmRZZlNfMSbkxaWBTUtKsuqctPX8MuAJQXGfCvn1yEwTEHZrPnhc4DaQlevva1r+Gll17C3r17UV9fH/W+NTU16OrqCritp6cHOp0O5eXlIfc3GAwwGEKb6YiIEtU5bMPwuBM6jYS5VZHLMAu8TbtimZ04HiDdZSMA+NLls9L+PYkyJaVlI1mWcffdd+P555/HX/7yF8yaNfF/XKtXr8bOnTsDbnv99dexcuVK5OVxdwERpZ44SXpOZREMusgHLIqel45hG3osNticytRRqk+UJjrfpTR4ueuuu/D000/j17/+NUwmE7q6utDV1YXx8XH1Plu2bMFtt92mfr5p0yY0Nzdj8+bNOH78OJ544gk8/vjjuPfee1N5qUREquOdE5eMAMBszEN9aT4A4EDTAABAIwEmQ9oq8kTnpZQGLw8//DCGh4dx1VVXoba2Vv147rnn1Pt0dnaipaVF/XzWrFl45ZVXsHv3bixbtgw/+MEP8NBDD3FMmojSRu13idKsK4j7HDjbD0Bp1tVoQvvziCh5UvrrgWi0jWb79u0ht1155ZV47733UnBFREQTE9NDkcak/S2sNeFPx7sDghciSi2ebURE5MfmdONs7wiA8GcaBROZF3Gic3EajwYgOl8xeCEi8vNxzwg8MlBakIeqMMcCBAvOzpQw80KUcgxeiIj8fNShbAJfWGsOu1sq2IzyQhjzfP+UsmxElHoMXoiIvN5pGsAPXz4OALiwoSSmx2g1EhqrfdmXTOx4ITrfMHghIgLw2oed+LvH34bF5sJF00uw6Yo5MT/WvzeGZSOi1OMyAiI67z21/xy+89JHkGXgE4uq8V9/uxzGvMjL6YItqPFlXrigjij1GLwQ0XntwddP4qG/fAwAuOWS6fj+pxZDp40vKb3AP/PCaSOilGPwQkTnrcMtg2rg8k+fmI+7r5kbU5NusIU1LBsRpRN7XojovPWzXUrgctNF9fjaunkJBS4AUFyQh4Yy5ZiAyhjGq4locph5IaLz0rEOC/50vAeSBNx1dezNuZH8x98swwdtQ1haX5yEqyOiaBi8ENF56We7lazLDRfUYnZl0aSfb9WsMqyaVTbp5yGiibFsRETnnY97RvDK0U4AwF1Xz83w1RBRvBi8ENF55+HdZ9Sx6FjOLyKi7MLghYjOK60DY3jxSDsA4G5mXYhyEoMXIjqvPLznDNweGWvnVcR8BAARZRcGL0SUc1xuD378xxM4cLY/rsd1Ddvwu4NtAICvXTMvFZdGRGnA4IWIcs7uk7342a4zuOuZ92BzumN+3K8OnIPD7cGqmZwMIsplDF6IKOe0DIwBAPpHHfi9t38lFqe7RwAAn7ywNiXXRUTpweCFiHJO+9C4+udfvNEEWZZjetzgmAMAUFHELbhEuYzBCxHlnA6/4OV0zwj2nOqN6XH9o0rwUsrDE4lyGoMXIso5IniZXVEIAHj8zaaYHjfoDV7KChm8EOUyBi9ElFFvn+3Hyx90xvUYUTa67/oF0EjAG6f7cKLLEvUxbo+MoXEnAAYvRLmOwQsRZYzbI+POpw7irl+/hx6rLabH2Jxu9I0oGZRLZpVhwxKl+fbxN6JnX4bGHBCtMSUFeYlfNBFlHIMXIsqYpr5RWG0uAECv1R7TYzqHlSCnQK9FSUEe7lw7CwDw+yMdUQOgAW/JqDg/D3la/tNHlMv4XzARZczxTl+pRwQxE2kfVEpGdSX5kCQJy6eXYsWMUjjcHvxqf3PExw2w34VoymDwQkQZk0jwIpp160ry1dvuvFzJvjx9oBnjjvBL6xi8EE0dDF6IKGNOdFnVP1u8zbQTEc260/yCl/WLa1Bfmo/BMSf2nOoJ+7iBMY5JE00VDF6IKGMCMy+xBS8davBiVG/TaiQs8x6y2D4Uvu9FjEmXM/NClPMYvBBRRgyNOdTmWyCOnpcwZSMAqDIpwUykxl91QR2DF6Kcx+CFiDLiWGfgXharPb6el2lBwUulSVn5H2niyLegjmPSRLmOwQsRZcTxTmvA57H0vHg8Mjq82ZrQzIsSvEyUeSkr5LlGRLmOwQsRZYTodxHTP7GUjfpG7XC4PJAkoKbYGPC1ygmCF3EoIzMvRLmPwQsRZYQIXi6eWQoAsMTQsNvhbcatNhlDFs1VmUXZKHzwMjDCzAvRVJHS4GXv3r248cYbUVdXB0mS8OKLL0a9/+7duyFJUsjHiRMnUnmZRJRmTrcHp7tHAACrZpUDACwxZF7UfpfS/JCviYbdgVEHHC5PyNfFqHQZR6WJcl5Kg5fR0VFceOGF+OlPfxrX406ePInOzk71Y968eSm6QiLKhLO9o3C4PSgy6LCo1gwgtlFp/+26wUry86DTSACA/tHA7MuYwwWbUwloyooYvBDlOl0qn3zDhg3YsGFD3I+rqqpCSUlJ8i+IiLKCKBktqDHBnK/8MxRLz4tvTNoY8jWNRkKlyYDOYRt6LHbUFvsCHLFdV6/VoFCvnfT1E1FmZWXPy/Lly1FbW4t169Zh165dmb4cIkrQX050Y++p3pDbRfCysNYMs1FpoI0l8yLKRvVhMi9A5ImjwVHlucsK9ZAkKcarJ6JsldLMS7xqa2vx6KOPYsWKFbDb7fjVr36FdevWYffu3bjiiivCPsZut8Nu9/1DZbFYwt6PiNKr22LDV546BFmWsfveqzG9vED92rEwwYvN6YHD5YFeF/l3qkgL6gTfrpfA4EWUkbigjmhqyKrgpbGxEY2Njernq1evRmtrKx544IGIwcvWrVvxve99L12XSEQx2nemD26PDAD4xZtn8f1PL1G/Jna8LKw1ocjo+2fIanOivCjyNFC4Qxn9VXqbdoMX1XFMmmhqycqykb9LL70Up0+fjvj1LVu2YHh4WP1obW1N49URUST7Pu5X//w/B1vRP6JkP3qtdvSN2CFJQGONCVqNpPahROt7GXO4MDimlH/CTRsBkXe99HNMmmhKyfrg5fDhw6itrY34dYPBALPZHPBBRJklyzL2nVGCF5NBB5vTg6f2NwPw9bvMKi9EgV7JupjzlYxItF0vYseLyaBTS03BqiKUjdTMSwEzL0RTQUqDl5GRERw5cgRHjhwBADQ1NeHIkSNoaWkBoGRNbrvtNvX+27Ztw4svvojTp0/jo48+wpYtW7Bjxw7cfffdqbxMIkqyloExtA+NI08r4V9vXAQAeGr/OYw5XL5Jo1qTen+TceKJo4n6XYDIwcsAjwYgmlJS2vNy8OBBXH311ernmzdvBgDcfvvt2L59Ozo7O9VABgAcDgfuvfdetLe3Iz8/H4sXL8bLL7+MjRs3pvIyiSjJRNZleUMp/vqievx018do7h/Dbw+2+SaNanxZUlMME0fRFtQJomzUFzF4YeaFaCpIafBy1VVXQZbliF/fvn17wOf33Xcf7rvvvlReEhGlwVsf9wEAVs8ph1Yj4StrZ+NfXvwQj71xFsY8pb9lYa0veDF7My/Rtux2RNnxIlSZla/1Wu2QZVkdi2bmhWhqyfqeFyLKLbIsY78383LZ3AoAwGdX1KO8UI+2wXF83KMcC7CwLjTzEu1k6WjbdYUK7/Zch9uDoTHfc4ngpZSZF6IpgcELESXVqe4R9I86YMzTYFlDCQDAmKfFHWtmqvcxG3Wo8zsVOp6el2lRgheDTosSb1Nu74ivdCSmlMqZeSGaEhi8EFFSiZLRxTPLAhbO3bp6BvL9Skb+m27FtFG04KVjeOLgBfBr2rUowYvbI6vTRsy8EE0NDF6IKIDN6caW54/iz8e7E3r8vqCSkVBSoMcXLpkOAFg5szTgaya15yV82cjtkdHpHZWOVjYC/LfsKvcfHndCtN6V8kRpoikhqzbsElHm7TvTh9+804L3mgexbmF1XI91uT14+6wSvKyZUx7y9f9vwwKsnFmGtfMCA5uJpo16rXa4PDK0GgnV5sgNuwBQZfI17QLAgPdoALNRhzwtf18jmgoYvBBRgD7vNtqWgbGAiZ1YfNhhgdXugtmow+K64pCv52k1uH5JTcjt5gl6XkS/S43ZCK0m+vUE73oZ8DuUkYimBv4aQkQBBr2TOeNOtxrIxGrfGaXf5dLZ5RMGGf58J0tHD14m6ncBQo8IEJkXBi9EUweDFyIKMDDmC1haB8ci3u9/3m3F6x91BexyEiPS4UpG0UzU8xLLgjohuOeFmReiqYdlIyIKMDTqCyBaB8Zw0fTSkPuc6rbivh0fAABWzijFv3xyERbWmvDuuQEAwJqgZt2JmCbIvMSyoE6oDCkbMfNCNNUweCGiAP6ZlzbvYrhgp7qt6p8PNg/iMz97C5fMKoPN6UFFkQHzqori+p7mfNHz4gzbZ9MRw7lGQmjDrhKMlTJ4IZoyWDYiogBD/mWjgfBlo3N9owCAqxsr8dcXTQMAvN3kzbrMKY+ryRfwZV6cbhk2pyfk620xbNcVqsxK5sVqc8HmdKs7XsoZvBBNGQxeiCiAWKUPRO55OesNXlbOLMODn1uGl+6+DBfPLIUkKUcBxKtQr4Xo7w03Lt1t8e54KZ44eDEZdDB4l+P1WOzoF0cDcMcL0ZTBshERBRj0OxOoJULmpckbvMwsLwQALK0vwW83rYHT7Ulol4okSSgy6GCxuWCxuVDlO/bImz1RrqnaPPF6f0mSUGU2oHVgHL0jNnV6ij0vRFMHMy9EpPJ45ICyUceQDS53aBlHlI1mVRQG3D6ZJXC+IwICMy+id8Wg06A4P7b1/pVFviMCBhi8EE05DF6ISGWxOeHxTj7naSVlLf+wLeA+Q2MONRMys6Igad9bPVk6aOKoy1syqjYbY+6lEU27PVYGL0RTEYMXIlKJoKTIoEN9qRKYBPe9iJJRjdmIAn3yKs++k6UDMy+i36VmgmMB/Imm3ZaBMYw73QAYvBBNJQxeiEglshSlhXloKFOCl7aBwHHppgglo8mKtGW3y5v5qYqh30UQZaOTXcpId55W6akhoqmBwQsRqQb9JnMavNtsgzMvot9lZtKDF++W3fHkZV5OeIOXskJ93OPbRJS9GLwQkUrsRCkt0KuZl+CJIzEmPTvJwYspwuGM3RalYXei06T9iZ6XvhHlsRyTJppaGLwQkUoEL2WFejSInpeB8D0vyc68+I4ICMy8qA27xbEHL+KIAKG8iMEL0VTC4IWIVKJht6QgDw1lomzk63mRZTnimPRk+Y4ICMy89IjgxRR7z0tV0H2ZeSGaWhi8EJFKXehW4Mu89FrtsHkndnqtdow63NBIwPSy5I1JA/6j0r7MiyzLaualJo7Mi9LjEvg5EU0dDF6ISCWmjUoK9SgpyIPJO6HT5m3aFSWj+tIC6HXJ/edD9Lz473mx2FzqWUfx9LzotBqUF/qyLwxeiKYWBi9EpBrylo3KCpTpnPoy0feilI5SNSYNhB+VFpNGxfl5MOZp43o+/9IRgxeiqYXBCxGpBsZ8e14AhIxLpzJ4MYUZlU5kTFqoZPBCNGUxeCEi1dBY4AnM6rh0fzqCl9Bpo0QW1AkBmRc27BJNKQxeiAiAciijmDYSmYpImZdkj0kDvmmjEbsLsqwcsJS0zAtHpYmmFAYvRARA6TVxe09lLCnwlo38el7cHhnN3p0vyV5QB/h6XjwyMOpQppsSWVAnMPNCNHUxeCFK0L2/fR+f/umbsLvcmb6UpBAL6gr1Whh0SnOsGrwMjqFjaBwOlwd6rQZ1JflJ//4GnQZ5WmW+WfS9JLKgTqjyC3hK2fNCNKUweCFKgN3lxvPvteH9tmGc7h7J9OUkha9Z1/dGL3a9WG0uvN82BACYXl4ArSb55wRJkuTX96JMHCWyoE4QZSOTUYc8Lf+pI5pK+F80UQLaBsfhrbCg13t+Tq4LbtYFgHy9FhXeE5r3nuoFkJpmXcGsnm8UmHmJZ0GdsKSuGEummXHTRfXJu0Aiygo8I54oAc39o+qfey1TI3gZGFUChuASS0NZPvpG7Nh7qg9AaoMX/8yL2yOj15p4z0u+Xov//drapF4fEWUHZl6IEnCuz3dY4VTJvIijAUq9zbqCKB2JLEhqgxexZdeJvhE7PDKgkaBmf4iIAAYvRAnxz7yIvoxcNximbARAPaBRSG3ZSJxv5FLHpCtNhpT02BBR7kpp8LJ3717ceOONqKurgyRJePHFFyd8zJ49e7BixQoYjUbMnj0bjzzySCovkSgh5/p9mZce6xTJvHiDl+BttCLzIqQj82K1OdUFdYnseCGiqS2lwcvo6CguvPBC/PSnP43p/k1NTdi4cSPWrl2Lw4cP4/7778c999yDHTt2pPIyieIW0PMyVYIX0fMSVDbyPz26QK8N2J+SbOrJ0uMudHtf1yoGL0QUJKUNuxs2bMCGDRtivv8jjzyC6dOnY9u2bQCAhQsX4uDBg3jggQdw0003pegqieLjcnvQNjiufj5VMi/hRqUB364XAJhZXghJSl0Jxz/z0j2sfB9mXogoWFb1vOzfvx/r168PuO26667DwYMH4XQ6wz7GbrfDYrEEfBClUseQDS4xJw0l8yLW2ecyX8NuYPBSW2xUe05mVaauZAQA5nzftNFkxqSJaGrLquClq6sL1dXVAbdVV1fD5XKhr68v7GO2bt2K4uJi9aOhoSEdl0rnsXPektE075bZcacbI3ZXJi8pKcS5RsHBi06rQa03gEjFsQD+AjIv3uAllWUqIspNWRW8AAhJSYvfaCOlqrds2YLh4WH1o7W1NeXXSOc30e+yqM6MQr2yRj/X+15kWVaX1AU37ALAnMoiAMC8alNKr8Osjkr7po2YeSGiYFm1pK6mpgZdXV0Bt/X09ECn06G8vDzsYwwGAwwG/mZG6SMmjWaWF+DjHiOa+kbRY7VjtvcNPhdZ7S61FFYS1LALAP9yw0JcPrcC1y+uSel1mNUldc5JHcpIRFNbVmVeVq9ejZ07dwbc9vrrr2PlypXIywv9B5UoE0TmZUZ5ISq9y9Mmk3l59p0WXL9tL9oGxya+c4qIfpcCvRbGPG3I1+dVm/CVK2ZDr0vtPxli2qjXasew93BGBi9EFCyl/xKNjIzgyJEjOHLkCABlFPrIkSNoaWkBoJR8brvtNvX+mzZtQnNzMzZv3ozjx4/jiSeewOOPP4577703lZdJFBdf5qUQlWYleJnMxNEzb7fgRJcV+z7uT8r1JWIgQrNuuomeF9F/Y8zTqKUkIiIhpf8qHDx4EFdffbX6+ebNmwEAt99+O7Zv347Ozk41kAGAWbNm4ZVXXsE3v/lN/OxnP0NdXR0eeughjklT1nB7ZLR4g5cZ5QWTzrzIsowzvcqp1BZb+Im6dBgSzbqFmc1wimkjocZsTOloNhHlppQGL1dddVXUEdLt27eH3HbllVfivffeS+FVESWuy2KDw+1BnlZCbbERVWrmJbEjArosNow53ACQ0YmlbMu8CFxQR0ThZFXPC1G2a+5T+l0aSgug02pQZVLeXBPNvHzcM6L+2WrLXPAS6VyjdMvTamDM8/2zxAV1RBQOgxeiOJzzKxkByqGBQOLByxm/4GUkC4KXcGPS6SaadgGg2sxJQiIKxeCFKA7NA75JI8C3QC3Rht0zvb4zkqz2zPW8DIyGX1CXCf4Nupw0IqJwGLwQxaG5z7fjBfBlXgZGHXC6PXE/X7aUjYbUc40yv5IgMPPC4IWIQjF4IYqDOBpghndNflmBXj33p28k/uyLmDQCMhu8ZEvDLhDYtMvtukQUDoMXohjJsoxmvx0vAKDRSKgoUt7w4+17sdicAeUmazaMSmdB8OI/Ls2GXSIKh8ELUYx6rXaMO93QSL5DGQGoE0c9lviCF/9mXSDDo9JZVDby73mp5KGMRBQGgxeiGIlJo2ml+QFr8kXTbm+cZSPRrCsCoUyVjWRZVo8HyKZpo9KCvLBHFRARMXghipHodxElI0FkB+LOvHj7XZZNLwEAjDnccHsiL3VMlRG/QxmzoWxkMiiZFzbrElEkDF6IYuQ7kLEg4HbfuHR8W3bFpNGy+hL1tkzsehn0jknn54U/lDHdir2nWjN4IaJIGLwQxehcULOukOiiOpF5WVhrVstQmdj14tuum/l+FwBYt7Aaa+dV4PY1MzJ9KUSUpXhcK1GMfJmX4ODF27AbR/DidHvUAx7nVBXCbNShb8SRkb4XX7Nu5ktGgNID9KsvX5LpyyCiLMbMC1EMAsekA8tGiWRemvtH4fLIKNRrUWM2qk2qmZg4yqZmXSKiWDB4IYrB4JgTVpsLkgQ0lIXveem12qOeou7v4x4lizOnqgiSJKHI26SaiV0vg94dLyVZ0KxLRBQLBi9EMRCTRjVmY0hTq8i8ONweWMZjy5yIfpc5lUUAfFtlM1E2UjMvWdLzQkQ0EQYvRDGINGkEAMY8rbpYLdaJI7Ggbk6l0j/jy7xkIHjJsp4XIqKJMHghikHrwDgAYHpZaPACAFXesd5Y+15E5mVulci8KFmPjAYvLBsRUY5g8EIUA9FIG6kvpLJI7HqZOHiRZVndrhtcNhrJwKi0eigjMy9ElCMYvBDFYNQbvBTowy9xqzLHPnHUbbFjxO6CViOpY9eZ7XkRhzKy54WIcgODF6IYjDncAIBCffjVSL7My8Q9L6JkNKOsQF1OJ3pe0r1hV5ZltA8pJbHa4vwJ7k1ElB0YvBDFQM28GKJnXmIpG4ngZba3ZAT4el4saQ5eBkYdGLErI+D1pQxeiCg3MHih807n8Dgsce5TmSjzUmWKvWFXnGk0p8q3qbcoQz0vzQPK4r1wI+BERNmKwQudVz5sH8ZVP96Nv/vF23E9btQRvedFPVk6jszLnIDMS2Z6XsQRBZGmqIiIshGDFzpvON0e/PPvPoDd5cHJLmtcjx2zezMvhkiZl9gbds94t+uKMWkA6p6YdB8PII48CLe/hogoWzF4ofPGf+85g+OdFgCA3eWBw+WJ+bGxZl6Gx52wOd0Rn8dqc6LLojT1zqnwBS9FhszseWkeCH/YJBFRNmPwQueFj3useOjPHwfcFk+WQ+15iZB5Kc7Pg16r/OfUNxI5+yL2u1QUGVDsN5qs7nlJc/DSOsCyERHlHgYvNOW5PTLu+90HcLg9uLqxUs2exHMI4kR7XiRJiqnv5d2mAQDA4jpzwO2iYdfh9kTN3CQby0ZElIsYvNCU99T+c3ivZQhFBh1++FcXxN0c63J7YPeWmCJNGwG+0lG0vpc9p3oBAFfOrwy4vcjvedNVOhp3uNVAi5kXIsolDF5oSmsdGMO/v3YSAPCtDQtQV5If9yGIY36ZkEh7XoCJJ47GHC684828XNkYGLxoNJJvUV2amnZbvCUjs1EX8dgDIqJsxOCFprT//NMpjDvdWDWrDLesmg7A/xDE2MpGYtJIp5HUvpZwJpo4OnC2Hw63B/Wl+ZhdEdog68sIpWfXi++kbDbrElFuiZwDJ5oCxEj0V9bOhkYjAfA/BDG2DIf/pJEkSRHv51tUF/6IgD0nfSWjcM8z2SMC2gbH8ODrpzBid0EGIMsAIGNhrRmbPzE/5HuKzMt09rsQUY5h8EJTWof33B7/1ffx9rxMtONFqClWMi/HO8PvkInU7xJ8XYkeEfDg66fw/OH2kNv/dLwH1y2uwZJpxQG3i+BlBvtdiCjHMHiJk9sjQ6uJ/Ns3ZY9xhxuDY0oJps7v0EGTd6dKIpmXaK5qrIJOI+FI6xA+bB8OCBbO9Y3iXP8YdBoJa+ZWhH18kTG+6/I3PO7Ey0c7AQD/9In5KC8yQJKAX+1vxrFOCw63DoUEL83crktEOSotPS8///nPMWvWLBiNRqxYsQJvvPFGxPvu3r0bkiSFfJw4cSIdlxrR4KgDq374Jyz4P6/C5Y59uRllTuewknUp1GthzvfF6UVqhiPGnhdv8DJR5qXabMSGC2oBAL/cdy7ga3tPK1mXlTNL1fJQsMn0vPz+SDvsLg8W1Jhw9zVzccsl0/G3q6bj2oVVAIDDLYMhj2HZiIhyVcqDl+eeew7f+MY38O1vfxuHDx/G2rVrsWHDBrS0tER93MmTJ9HZ2al+zJs3L9WXGpU5Pw/9ow443TL6Rx0ZvRaKTceQ0ntSW5If0O8Rb9lo1Fs2mijzAgB3rJkJAPj9+x0Y8Pt74ut3qYr4WHOCi+pkWcZv3mkFANx8cUPAz7p8eikA4EjrUMBj3B4ZbYNixwsbdokot6Q8eHnwwQfx5S9/GXfeeScWLlyIbdu2oaGhAQ8//HDUx1VVVaGmpkb90Goze+KtViOp0yRdw+EbMim7dHgzL7XFxoDb422MVTMvUXa8CBdNL8EF04rhcHnwm3eUAN3ucmPfmX4Akftd/K/LGmfZ6MN2C453WqDXafBXy6cFfO3ChhIAwNneUQyN+YKpjqFxON0y9FoNasyBrw8RUbZLafDicDhw6NAhrF+/PuD29evXY9++fVEfu3z5ctTW1mLdunXYtWtXxPvZ7XZYLJaAj1Sp8v4j321h8JILRLPutJL8gNvNcY5Kq5mXCcpGgLJpV2RfnjnQDJfbg4PnBjHudKPSZMDCWlPEx/pGuOMLXn7zrhIkbVhSE7KvpaxQj5nespB/9kWUjOrL8tnDRUQ5J6XBS19fH9xuN6qrqwNur66uRldXV9jH1NbW4tFHH8WOHTvw/PPPo7GxEevWrcPevXvD3n/r1q0oLi5WPxoaGpL+c6jX7c28MHjJDZ2ibFQcGLzEOyrty7zElv375IW1KC/Uo2PYhp3HugOmjKKNWvuW58Xe8zLmcOGlIx0AlJJROOFKR+qxAGzWJaIclJZpo+B/sGVZjviPeGNjIxobG9XPV69ejdbWVjzwwAO44oorQu6/ZcsWbN68Wf3cYrGkLICpKRaZl8jr3yl7iLJRXUlQ2SjenheH6HmJ7T8Xg06Lv101HT/d9TG27zuHIe/EU7SSERB/Lw4AvPxBJ0bsLswsL8Dq2eVh77N8egleONyOwy1D6m0tPJCRiHJYSjMvFRUV0Gq1IVmWnp6ekGxMNJdeeilOnz4d9msGgwFmszngI1WqWTbKKaJsVFcSnHmJrzwzZhfTRrH3Xf3dpTOg1Uh4u2kAJ7ut0EjA5RFGpH3XFf/xAM++qzTqfi6oUdffMm/fy5HWIcjK5jq0DCjbdaezWZeIclBKgxe9Xo8VK1Zg586dAbfv3LkTa9asifl5Dh8+jNra2mRfXtxE8NLF4CXrybKMzmFRNgrfsBtzz0ucmRdAydJdv6RG/fzChhKUFkY/PyjeYwtOd1txqHkQWo2Ez15UH/F+C2rMMOg0GB53oqlPCVpYNiKiXJbystHmzZtx6623YuXKlVi9ejUeffRRtLS0YNOmTQCUsk97ezueeuopAMC2bdswc+ZMLF68GA6HA08//TR27NiBHTt2pPpSJ1Rt9h68x7JR1hsed2LMG3QEZ17MfhmOaCVMwbfnJb6Jty+umYmXP1AWx01UMgLin4J6zpt1WbegSm0mD0ev0+CCacU42DyIwy1DmFVRiBYRvHDHCxHloJQHLzfffDP6+/vx/e9/H52dnViyZAleeeUVzJgxAwDQ2dkZsPPF4XDg3nvvRXt7O/Lz87F48WK8/PLL2LhxY6ovdUJq2SjC2TWUPcSOl7JCPYx5gUGHyHB4ZGDM4Z5w+Zxvz0t8/7msmFGKi2eW4nDLEDYsmThzGE/Pi8vtUY8C+PyqiXu8ljWUKMFL6yCuXlCljmM3MPNCRDkoLQ27X/3qV/HVr3417Ne2b98e8Pl9992H++67Lw1XFT8RvAyNOWFzukPeFCl7dEZo1gUAY54GWo0Et0eG1eaaMHiJd9pIkCQJT35xFQZHHTEFCaKReMThgscjqwdJhtNjtWNg1AGdRsIV8ybO6igTR0040jqkniZdbTbw7zAR5aS0HA8wVZiNOhjzlJeMpaPsJpp1g8ekASWo8DXHTtxfMhLHnpdgRQZdzNkNsX9Gln3nKUXSN6L8/Ssv0kOnnfg/4+XTSwAoh0ae6lYOjpxRxmZdIspNDF7iIEkSm3ZzRIe3WbeuOHwviOgvieUE50QzL/Ey6DTI0yrZlokmjkTwUlFkiOm5a4uNqDIZ4PbIePmoMv3HM42IKFcxeIkTx6VzQ6QxaUH0vcTSHJtoz0u8JEnym4SaIHixKqv+Yw1eJElSsy9veg+J5KQREeUqBi9xyqbgpWNoHD/9y2kM8qDIEJ1+hzKGY4oxSAB8mZdIp0EnU6zj0r1xZl4A36Zdj7LqhZkXIspZaWnYnUqy6YiAn+/+GE8faIHTLeObn5if6cvJKmK77rQwDbuA/2RP9CDB45HVkeuCOEelExFz5kUEL6bou2P8iWV1Ak+TJqJcxcxLnLLpiIDT3SMAgGOdqTuMMhe5PbJ68ne4hl0g9m224063+udYTpWerFjHpftGlGxbZRyZl6X1xfAfYOLRAESUqxi8xClZJ0t7RO5+Es55R17F9Agp+kbscHlkaCSgyhT+zV2MJU/UsCumfiQJ6qRZKsUaVPVZ4y8bFeh1WFCjHJ9hMuhQWpCX4FUSEWUWg5c4JaNs9P0/HMPFP/yTmh1IxJjDpWZ/WgbG1L4MAtq9zbo1ZmPEMeJYG3bHvM26hXrdhJt4kyHWnpd4p42EZd6m3enlBWn5eYiIUoHBS5z8y0bikLt4yLKMHe+1oX/UgYPNAwlfhzibRnlO4OOekYSfa6qZqFkXiP18I5F5KUjxmLQQ6xEBifS8AMA1jVUAlO2/RES5ig27caoyKcHLuNMNi82F4vz4Uu9tg+MYHlfeMCeTeTnnPWBPONllxdL6koSfbyrxbdeNHLyYY+wtEc26E23hTRZTDOUsp9uDwTHl71C8mZdrF1Xjf792OWZXslmXiHIXMy9xytdr1Te+ngRKRx+0Dat/nkzp6Zxf5gVg34s/UTaKtKAO8FvFP0Fvyag9zZmXGK5rwDsar5GA0oL4Mi8AsGRaccp31hARpRKDlwRMZuLoaLsveOmaxMSSyLzUeBuIT3azbCSoZaMowYvJEFtviZp5SdObfSw9L73eZt2yQgO0Uc4/IiKaqhi8JGAyRwR86Be8dE+ibNTknTRav7gaAHCqi5kXoSOGspE6khxr5iUNO14AXzkrWubF16wbf9aFiGgqYPCSgES37MqyHJR5STx4EScDr19Uoz7X8NjEhwyeDzq8mZdowUtRvD0vacq8xLKkTt3xEmEMnIhoqmPwkoBqs/KmEW/PS+uAr1kXUAKORCaW/Mekl0wzY5r3TfpUD7MvdpdbzUxEb9iNbVQ63dNGvrJRLJkXBi9EdH5i8JKARMtGIuuyoMYEAHC4PBhKIFsixqSL8/NQUqBHo/f5TrB0pE5wGXSaqEvYRIZj3OmG0+2JeD91z0uapo1iyrxYWTYiovMbg5cE+MpG8TXcftA+BAC4aEap+saaSOlINOvOrFDGXedXK8EL+158JaNpJflRl7CJshEQPfuS/szLxPtnROaFZSMiOl8xeElAtJ6XUbsL/SPhgxrRrLt0WvGkmn5Fs+4s76nAjTVFAICTHJdGh3dMujbCgYxCnlajrvuP1hyb7syLCF7sLg8crvAZIdHzwrIREZ2vGLwkQO15sdoDziiSZRm3PHYAV/14N9oGA/ewyLKMo94dL0umFfvGrROYOGruU547JPPSbU2oh2YqEQvqIh3I6E/0l1iiZDkytWEXiBxUseeFiM53DF4SUFlkgCQppxf3exeGAcCZ3hG83zYMq92F3x5sC3hMy8AYLDYX9FoN5leb1P0sieyKEZmXmeVK8DKnsggaCRgac6o7QOJxrMOCax7YjacPNMf92GzTMTzxpJGgHoIYpWyU7mkjnVaD/Dxt1Oti8EJE5zsGLwnQaTXqG4d/6egvJ3rUP//uUFtAVkY06y6sNUGv00yqbBTc82LM06p/jrd0JMsy/s/vP8TZvlH8/kh73Ncykdc+7MKqH/4JB88lfo5TPDpi2K4rmGJojk33nhfA/4iA0IyQ2yOrG3bjPdeIiGiqYPCSIFE6ihS8tA+NY9+ZfvVz/5IR4L+lN77gZczhQo83uzKr3Hc+TaO3dHQyzqbdl4924lDzIACoz5tMvzvUih6rHc+925r05w6nM4YdL4J6snS0npc0Z16A6EcEDIw64JEBSQLKEjgagIhoKmDwkqDgss/wuBMHzylBwFWNlQCA/znoe8MWmZcLRPAiMi9x9ryc8/a7lBTkodhvFNi/7yVWNqcbP3r1hPp5T4InZUdzvFO5nvdaBpP6vJH4tutOnHmJ5WTpdPe8ANF3vYiSUVmBHjot//MlovMT//VLUFVQ2eeN071weWTMqSzEP32iEQDw2kddGB5zBmzWvaBeCV4S3dJ7LqjfRRC7XuI54+iJt5rQNjiujtyOO90THlQYD4vNqR6SeKZ3FENjjgkeMTlWm1N9w4+tYXfiE5zTPW0E+B8REBpUsd+FiIjBS8KqTUrwIbbsipLRuoXVWDLNjAU1JjhcHvz+/XY094/BanNBr9OoGRJRNuofdcDucsf8fUXwMqsiMHgRz3u62xrQaxNJr9WOn+86AwDYsmGBmoVIpOE3kuAS1uGWoaQ9dzhix0txfl5MwUYsJzhnIvMSbVGdGryw34WIzmMMXhJUU6z85ttlscHjkbHnZC8A4OrGKkiShM+tbACglI7UZt0aE/K8qf7SgjzovX/uiWPiSG3WDcq8zCwvgF6rwZjDrWY7onlw5ymM2F1YWl+MzyybhiqTb/w7WU50WgI+T2XpqMdqwz//7n0AwOzKwgnurZjoBGdZln09L2nMvJiinLvUZ+WOFyIiBi8JqvLreXm/bQj9ow6YjDqsnFkKAPjM8mnI00r4sN2i9r6IkhEASJKEqjBNvxM5p+54KQi4XafVYE6Vd1ndBE27J7oseO7dFgDAv9ywCBqNpJaOkhm8HPdeh3ijTVXwcrLLir/62T580DaM0oI8/OsnF8X0OPMEo9J2lwdubxYrvZmXiXteGLwQ0fmMwUuC/MtGomR0xbxKNbNSVqhXT3x+43QfAF+zrlCTwLh0pJ4XAGisjm3T7o9ePQGPDGy8oAarZpUB8AVj8R42GY3IvHxuZT0A4EjLkBoMJMveU7347MP70D40jlkVhXjhq5dh+fTSmB470TlCo37lpII0ThtFOyKgl8ELERGDl0T596z88aMuAMDVC6oC7vM33jdtYUlQ8FJdHN/E0ajdNyYdLniZXzPxuLTd5cab3mBqs7exGFAW7wHJ63nxeGT1Oj61rA5FBh1GHe64pqEm8tL7Hfji9ndhtbuwalYZnv/HNeq+m1hMdIKzKBkZ8zTQaiKfk5Rspii9OL6jAdjzQkTnLwYvCfLvWTnVPQJJ8o1IC2vnVarZFf9mXaEmzokjkXUpDRqTFhpjGJdu6huFyyPDZNRhjl9vSJU5uWWj9qFxjDrc0Gs1mFtZhAsblMAtWaUjWZbxo1eOw+2R8ZlldfjVl1ehtDC+N3TRsGuN0LArmnXTueMFmKjnRTTsMvNCROcvBi8J8u9ZAYAL60tCUvlajYTPrlCyLwtrzWpJSfCVjWILGJr7A880CiaCozO9I3C6wx/qd8o7St1YbQo4dVk07CYr83LcWzKaW1UEnVaDi7ylnPeah5Ly/G2D4+gYtkGnkbD1r5fCoIu/J2WiE5xHvWPS6dyuC/h6XsL14qgnSrNsRETnsfT+SjnFVJuNaBtUJnuuCSoZCV9ZOxvdFhv+avm00MfHeThjU584TTp88FJfmg+TQQer3YWPe0awsNYccp9T3lKOKDEJVaKHx5qcnpcT3u+zoFb5PiJ4OZykzMu73uMGLqgvRn6CzbTieIBIo9JjGc68BB8P4PE7S4s9L0R0PmPmZRKq/TIvkYKX4oI8/PhvLsSauRUhX1PLRjEGDMFnGgWTJAkL65SA5aMOS9j7iGbe+d7JJCHZZaMTXcr3X1ijXM/y6SUAgLN9oxgcnfyyunealOBl1cyyhJ/Dv+cl3GZhNfOSxkkjIPL+mcExh9rwXM6eFyI6j6UlePn5z3+OWbNmwWg0YsWKFXjjjTei3n/Pnj1YsWIFjEYjZs+ejUceeSQdlxk3sSW3ymTA4rrQLMdE/I8IiGUtvygbzSgviHifxWrwMhz266IfJjjzIsoQQ2POuJbmRXKiMzDzUlKgV/evHG6dfPblHW/m5eJJBS9KkOD2yLA5Q8tsauYljTteAN//Fz0WO8Ydvv8vRLNuSUFeSAmSiOh8kvJ/AZ977jl84xvfwLe//W0cPnwYa9euxYYNG9DS0hL2/k1NTdi4cSPWrl2Lw4cP4/7778c999yDHTt2pPpS4ybKMjcsrQ3oH4mVyHbYXR4Mj0c+X0doirBd19/iOqUxNlzmZdzhRsuAEgA1BjUPl/g1IMfa9+Jye9QejODvI651QY0vqEtW30vfiB1ne5XnF3t1ElGg10IMEYXrexl1ZCbzUl+aj2kl+XC4PTjQ5DvckzteiIgUKQ9eHnzwQXz5y1/GnXfeiYULF2Lbtm1oaGjAww8/HPb+jzzyCKZPn45t27Zh4cKFuPPOO/GlL30JDzzwQKovNW43XVSPZ+68BN/asCChxxvztCj1Tg0F73p54s0mLPg/r+KC7/4Rl/6/P+OaB3arQcWMCD0vALDIG1Ad77CEHBPwcc8IZFkZsy0PegOUpPgX1W199QRW/fBP2HOqN+D2U91W9ftU+k3FqMHLJPte3vWWjBbUmFAyiZOVJUlSd72EO99ozJ6ZnhdJknDFfKXMKDY3A/7BC0tGRHR+S2nw4nA4cOjQIaxfvz7g9vXr12Pfvn1hH7N///6Q+1933XU4ePAgnM6JsxPppNVIuGxuRUKTLkJ1mNOlZVnG4282web0wGpzoctiw1lvv8vsikIU54eOSQvzqoug12pgtbvQOjgW8DW13yUo6yJUxjlx9MePuuCRgZ/v+jjgdtHv4p91AYCLZpQAAN5vndyyumSUjATR9xKuaVfNvKR52ggArpyvjN3vPe0LXsT/L8y8ENH5LqW/Uvb19cHtdqO6ujrg9urqanR1dYV9TFdXV9j7u1wu9PX1oba2NuBrdrsddrvvzdZiCd+omq1qio040WUN2PVyosuK9qFxGHQavHT35XC6PRi1uzDmdE/YW5On1WB+TRE+bLfgow5LQJbm1ATBSzznG/VYbeqk1dtNAzjWYcEi77UdF/0uQX0186pMKDLoMGJ34WSXVb1/vMSk0cWzkhG8RB6XzlTmBQDWzK2AViPhbO8oWgfG0FBW4LegjsELEZ3f0tL1F9wPIsty1B6RcPcPdzsAbN26FcXFxepHQ0NDEq44fXxNu76AQRw3cPncCjTWmLBkWjEumV2Oqxur1JHmaBbXKn0vx4L6XsTG24jBi7cHpzeGpXlHgk6I/uW+cyHfZ0HQqLZWI2FZQwkA4FCCpSOrzan+XJOZNBJE2SjcThVfz0v6gxezMQ8XeSe0RPZF3fHCBXVEdJ5LafBSUVEBrVYbkmXp6ekJya4INTU1Ye+v0+lQXl4ecv8tW7ZgeHhY/WhtbU3eD5AG1WHON/rT8W4AwDULw49fT2TxtPATR6e9mZfGmqKQxwBAZZHY9TJx5uVI6xAAX3blxSPtGBh1QJZlv7JRaJAk3pAPNycWvBxqHoRHBqaXFahHNExGtG22vmmj9JeNAOWsLMDX98KeFyIiRUqDF71ejxUrVmDnzp0Bt+/cuRNr1qwJ+5jVq1eH3P/111/HypUrkZcX2uthMBhgNpsDPnKJeAMWZaO+EbsaGKxbED7Am8jiMLteLDYnOrx9NfMmyLzEErwc9mZe7lgzExdMK4bd5cFv3mlBj9WOwTEntBoJc6tCg6TlMybXtPtuEvtdAL9dL+F6XuyZy7wAwJXe4yb2nemH02+yi2UjIjrfpbxstHnzZvziF7/AE088gePHj+Ob3/wmWlpasGnTJgBK5uS2225T779p0yY0Nzdj8+bNOH78OJ544gk8/vjjuPfee1N9qRlRE9Swu+tED2QZWDLNnHBmYUGNGZKkBCGiyVNkXWqLjTAbwzf8xnpEgNsj4/22IQDA8umluGPNTADA0weacbRNyfbMqiiEMS80Y3FRQykkCTjXP4bWgbGQr0/k3SYl6Fk1K/ERaX9F0XpeMpx5WVJXjLJCPUbsLrzXPIg+K3teiIiANAQvN998M7Zt24bvf//7WLZsGfbu3YtXXnkFM2bMAAB0dnYG7HyZNWsWXnnlFezevRvLli3DD37wAzz00EO46aabUn2pGVEddDjjn48r/S7XJJh1AZSlamIXjCgdnexSzjSK1O8CxH5EwKluK8YcbhQZdJhbVYRPXliLiiI9Oodt+C/v5FG4khGgbBxePVsp//3hg444firA5nTjiDdoSl7mJXLZKJM9LwCg0UhYO08Zmd59qhf9ozyUkYgISFPD7le/+lWcO3cOdrsdhw4dwhVXXKF+bfv27di9e3fA/a+88kq89957sNvtaGpqUrM0U5E4YqB/1IERuwtveJszr02w30UIXlZ3Su13iRK8eK+lb8QRdZRZlLWW1hdDq5Fg0Glxy6rpAJQxaABhz1USPr2sDgDw0pH4gpcP2obhcHlQUWSIuqgvHqYoDbvqtFGGMi+Ar+/lfz/ogNPtPRogztOziYimGu4Yz7CyQr262fYP73dg1OFGpcmAJd7gI1Gi7+VYUPASLfNSXqiHJClloYEo5w+JwxXFeUUA8IVLZ0Cn8U2DRcq8AMD1i2uRp5VwosuqTibFQvS7rJpVmtBG43B8PS/hykZK5iUTo9LCWu+yutYBZSzdZNSFLccREZ1PGLxkmCRJasbjmbebAQDrFlRBo5ncm3PwGUdq5iVK8KLTatTf6qOVjkSz7vIGX99JtdmIjRf4dvAEj0n7Ky7Iw1WNSmbppffbo/0YAcRhjMkqGQG+UenwZaPMZ16qTEZ1azLgO/eIiOh8xuAlC4im3Q/blSxJpBOq4yHKRuf6x9DcP4q+EQckCWEngPxVmqKPS1tsTnzcq/TPLPPLvADAFy+b6X0OA+omaDZWS0fvd8R0KKXbI+OQd7w6mcFL1FHpDE8bCWLqCGCzLhERwOAlK1T7vdHrdRpc7m3SnIyyQj1qvc/74mGlt2R6WQHyJzhkcKKJow9ahyHLQENZfsgb6fLppfjll1bhyTsunrCss25BNQr1WrQOjOOwt08mHKfbg6Ntw/j5ro8xYnfBZNBF7aeJl1o2Cpo2crg8cLiVk6YzWTYCfH0vAFBhYr8LEVFm/1UmAL7MCwBcNqc8ab/pL64zo3PYhhePKKWZaP0uwkTBi9rv0hB+VFmcyTORfL0W6xfX4IXD7XjpSId6aKPw/HttePbdVnzQNgSb06PefsnscmgnWVLzJzIvwWcbjXv7XcS1ZtKKGaUo1Gsx6nAz80JEBGZesoJ/8HLNwsRHpIOJXokm76GO0fpdBHVRXYQjAsSk0fKgklEiPuUtHf3vB51wuX0BynPvtmDz/7yPd5oGYHN6YDbqcMX8Snzj2nn44V8tmfT39RepbCT6XfRaDfS6zP5notdpsGauko1LxlZhIqJcx8xLFvAvG61LQr+LsChoYmledfR+F8DXEBqu50WWZbXEI84omozL51agtCAPfSN27D/bj7XzKvGnY93Y8vxRAMDtq2fg1tUzMLuiaNINzJGIht0xhxtuj6xmdcSCukycKB3O/7lhEWaUFagj6URE5zNmXrLAolozNJJy0GBdSX7Snjf4BOpoO16EKnPkht2WgTEMjDqg12oSPhHaX55WgxuWKhNKvz/SgUPNg7j7N+/BIwN/s6Ie3/3UYsytMqUscAF8G3aBwF0v4miATPe7CNPLC/Avn1yEkgL2vBARZce/zOe5uVVF+OM3rkj6acH1pfkozs/D8LgTOo2E2RUTZ15Ez0u4UWlRMlo8zQyDLjkZiU8vm4anD7TgtQ+78Kfj3bA5PbhmQRW2/vUFSdvlEo1Bp4Vep4HD5YHV7kRxgdLAK8pGBRnudyEiolDMvGSJedWmpP9WLUmS2vcyq6Iwpt4NcURAr9UeMsIs9rsko2QkrJheirpiI0bsLgyNObGsoQQ/vWU5dNr0/dU0h+l7UcekDYzviYiyDYOXKU6UjmKZNAKgZn9sTk/IScu+zbrJORQRUM7v+fTyaQCA2ZWFeOKOi9O+V8U3Lu1XNhIL6ph5ISLKOvy1coq747KZaB0cwz9eNSem++frtTAZdLDaXeix2NUTqG1ON451Kkv0licx8wIAd189F7XFRly/pAZlGTi3RzTtjvgdETCaJQvqiIgoFP9lnuLqSwvw37eujOsxlWYDrL0u9Fht6kbe95oH4XTLqCjSo740eU3FgHIK9m2rZyb1OeMRblx6LAuOBiAiovBYNqIQ4RbV/XL/OQDAJxZVp6WRNp3CnW/EzAsRUfZi8EIh/Jt2AaC5fxSvH+sGAHzpslkZu65UET0v/sHaGHteiIiyFoMXCuEbl1bezJ986xxkGbiqsRLzYmz8zSUrZigNyC8eaYfHo0xYqaPSnDYiIso6DF4ohP8RAcPjTvzPwVYAwJcvn3pZFwD4zPI6mI06NPePYfepHgC+UWlmXoiIsg+DFwpR6Zd5efadFow53GisNuHyuZM/7TobFeh1uPniBgDA9n3NAJh5ISLKZgxeKIToeekctmH7vnMAgC+vnTXlGnX93XrpTEgSsPdUL870jmDMwcwLEVG2YvBCIUTPS1PfKDqHbagoMuDT3hOgp6rp5QXqoZhP7TuHUbs4HoCZFyKibMPghUKIzItw66UzknaWUTa7Y43S0/O7Q21qszL3vBARZR8GLxTCnK9Tz0HS6zT4u0unZ/iK0uOyueWYW1WEUYcbbYPjAJh5ISLKRgxeKIQkSagsUkpHN100DeVFyT3tOltJkoTbV88IuI2ZFyKi7MPghcJav7galSYD/uGK2M5Emir++qJ6mPwmjAqZeSEiyjoMXiis79y4GO/cvw4zKwozfSlpVWjQ4W9WNqifF3DaiIgo6zB4oYim8mh0NLetngGNBOTnaVFkZOaFiCjb8F9moiAzKwrxqy9fAknCeTFlRUSUaxi8EIVx2RTdJkxENBWwbEREREQ5hcELERER5RQGL0RERJRTGLwQERFRTmHwQkRERDmFwQsRERHllJQGL4ODg7j11ltRXFyM4uJi3HrrrRgaGor6mDvuuAOSJAV8XHrppam8TCIiIsohKd3zcsstt6CtrQ2vvfYaAODv//7vceutt+IPf/hD1Mddf/31ePLJJ9XP9Xp9Ki+TiIiIckjKgpfjx4/jtddew4EDB3DJJZcAAB577DGsXr0aJ0+eRGNjY8THGgwG1NTUpOrSiIiIKIelrGy0f/9+FBcXq4ELAFx66aUoLi7Gvn37oj529+7dqKqqwvz58/GVr3wFPT09Ee9rt9thsVgCPoiIiGjqSlnw0tXVhaqqqpDbq6qq0NXVFfFxGzZswDPPPIO//OUv+I//+A+8++67uOaaa2C328Pef+vWrWpPTXFxMRoaGsLej4iIiKaGuIOX7373uyENtcEfBw8eBBD+VGJZlqOeVnzzzTfjhhtuwJIlS3DjjTfi1VdfxalTp/Dyyy+Hvf+WLVswPDysfrS2tsb7IxEREVEOibvn5e6778bnP//5qPeZOXMmPvjgA3R3d4d8rbe3F9XV1TF/v9raWsyYMQOnT58O+3WDwQCDwRDz8xEREVFuizt4qaioQEXFxCfurl69GsPDw3jnnXewatUqAMDbb7+N4eFhrFmzJubv19/fj9bWVtTW1sZ0f1mWAYC9L0RERDlEvG+L9/Go5BS6/vrr5aVLl8r79++X9+/fL19wwQXyJz/5yYD7NDY2ys8//7wsy7JstVrlf/qnf5L37dsnNzU1ybt27ZJXr14tT5s2TbZYLDF9z9bWVhkAP/jBD37wgx/8yMGP1tbWCd/rU7rn5ZlnnsE999yD9evXAwA+9alP4ac//WnAfU6ePInh4WEAgFarxdGjR/HUU09haGgItbW1uPrqq/Hcc8/BZDLF9D3r6urQ2toKk8kUtbcmERaLBQ0NDWhtbYXZbE7qc5MPX+f04OucPnyt04Ovc3qk6nWWZRlWqxV1dXUT3leS5VjyMwQo/4cVFxdjeHiY/2GkEF/n9ODrnD58rdODr3N6ZMPrzLONiIiIKKcweCEiIqKcwuAlDgaDAd/5znc4mp1ifJ3Tg69z+vC1Tg++zumRDa8ze16IiIgopzDzQkRERDmFwQsRERHlFAYvRERElFMYvBAREVFOYfAS5Oc//zlmzZoFo9GIFStW4I033oh6/z179mDFihUwGo2YPXs2HnnkkTRdaW6L53V+/vnn8YlPfAKVlZUwm81YvXo1/vjHP6bxanNXvH+fhbfeegs6nQ7Lli1L7QVOEfG+zna7Hd/+9rcxY8YMGAwGzJkzB0888USarja3xftaP/PMM7jwwgtRUFCA2tpafPGLX0R/f3+arjb37N27FzfeeCPq6uogSRJefPHFCR+TkffBuA4rmuKeffZZOS8vT37sscfkY8eOyV//+tflwsJCubm5Oez9z549KxcUFMhf//rX5WPHjsmPPfaYnJeXJ//ud79L85Xnlnhf569//evyv/3bv8nvvPOOfOrUKXnLli1yXl6e/N5776X5ynNLvK+zMDQ0JM+ePVtev369fOGFF6bnYnNYIq/zpz71KfmSSy6Rd+7cKTc1Nclvv/22/NZbb6XxqnNTvK/1G2+8IWs0GvknP/mJfPbsWfmNN96QFy9eLH/mM59J85XnjldeeUX+9re/Le/YsUMGIL/wwgtR75+p90EGL35WrVolb9q0KeC2BQsWyN/61rfC3v++++6TFyxYEHDbP/zDP8iXXnppyq5xKoj3dQ5n0aJF8ve+971kX9qUkujrfPPNN8v/8i//In/nO99h8BKDeF/nV199VS4uLpb7+/vTcXlTSryv9Y9//GN59uzZAbc99NBDcn19fcqucSqJJXjJ1Psgy0ZeDocDhw4dUg+RFNavX499+/aFfcz+/ftD7n/dddfh4MGDcDqdKbvWXJbI6xzM4/HAarWirKwsFZc4JST6Oj/55JM4c+YMvvOd76T6EqeERF7nl156CStXrsS///u/Y9q0aZg/fz7uvfdejI+Pp+OSc1Yir/WaNWvQ1taGV155BbIso7u7G7/73e9www03pOOSzwuZeh9M6anSuaSvrw9utxvV1dUBt1dXV6OrqyvsY7q6usLe3+Vyoa+vD7W1tSm73lyVyOsc7D/+4z8wOjqKz33uc6m4xCkhkdf59OnT+Na3voU33ngDOh3/aYhFIq/z2bNn8eabb8JoNOKFF15AX18fvvrVr2JgYIB9L1Ek8lqvWbMGzzzzDG6++WbYbDa4XC586lOfwn/913+l45LPC5l6H2TmJYgkSQGfy7IccttE9w93OwWK93UWfvOb3+C73/0unnvuOVRVVaXq8qaMWF9nt9uNW265Bd/73vcwf/78dF3elBHP32ePxwNJkvDMM89g1apV2LhxIx588EFs376d2ZcYxPNaHzt2DPfccw/+9V//FYcOHcJrr72GpqYmbNq0KR2Xet7IxPsgf73yqqiogFarDYnge3p6QqJKoaamJuz9dTodysvLU3atuSyR11l47rnn8OUvfxm//e1vce2116byMnNevK+z1WrFwYMHcfjwYdx9990AlDdZWZah0+nw+uuv45prrknLteeSRP4+19bWYtq0aSguLlZvW7hwIWRZRltbG+bNm5fSa85VibzWW7duxWWXXYZ//ud/BgAsXboUhYWFWLt2Lf7v//2/zI4nQabeB5l58dLr9VixYgV27twZcPvOnTuxZs2asI9ZvXp1yP1ff/11rFy5Enl5eSm71lyWyOsMKBmXO+64A7/+9a9Zr45BvK+z2WzG0aNHceTIEfVj06ZNaGxsxJEjR3DJJZek69JzSiJ/ny+77DJ0dHRgZGREve3UqVPQaDSor69P6fXmskRe67GxMWg0gW9zWq0WgC87QJOTsffBlLYD5xgxhvf444/Lx44dk7/xjW/IhYWF8rlz52RZluVvfetb8q233qreX4yIffOb35SPHTsmP/744xyVjkG8r/Ovf/1rWafTyT/72c/kzs5O9WNoaChTP0JOiPd1DsZpo9jE+zpbrVa5vr5e/uxnPyt/9NFH8p49e+R58+bJd955Z6Z+hJwR72v95JNPyjqdTv75z38unzlzRn7zzTfllStXyqtWrcrUj5D1rFarfPjwYfnw4cMyAPnBBx+UDx8+rI6jZ8v7IIOXID/72c/kGTNmyHq9Xr7ooovkPXv2qF+7/fbb5SuvvDLg/rt375aXL18u6/V6eebMmfLDDz+c5ivOTfG8zldeeaUMIOTj9ttvT/+F55h4/z77Y/ASu3hf5+PHj8vXXnutnJ+fL9fX18ubN2+Wx8bG0nzVuSne1/qhhx6SFy1aJOfn58u1tbXyF77wBbmtrS3NV507du3aFfXf22x5H5RkmbkzIiIiyh3seSEiIqKcwuCFiIiIcgqDFyIiIsopDF6IiIgopzB4ISIiopzC4IWIiIhyCoMXIiIiyikMXoiIiCinMHghIiKinMLghYiIiHIKgxciIiLKKQxeiIiIKKf8/54RKyHC4N1qAAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(t, data)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "e06eb458",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" message: Optimization terminated successfully.\n",
" success: True\n",
" status: 0\n",
" fun: 23.329378869713718\n",
" x: [ 2.960e+00 -5.312e+00 1.076e+01 3.056e-01]\n",
" nit: 215\n",
" nfev: 362\n",
" final_simplex: (array([[ 2.960e+00, -5.312e+00, 1.076e+01, 3.056e-01],\n",
" [ 2.960e+00, -5.312e+00, 1.076e+01, 3.055e-01],\n",
" ...,\n",
" [ 2.960e+00, -5.312e+00, 1.076e+01, 3.055e-01],\n",
" [ 2.960e+00, -5.312e+00, 1.076e+01, 3.055e-01]]), array([ 2.333e+01, 2.333e+01, 2.333e+01, 2.333e+01,\n",
" 2.333e+01]))"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"initial_guess = np.full(4, 0.8) #number of parameters\n",
"res_nograd = optimize.minimize(logp_data, x0=initial_guess,\n",
" args=(t, data),\n",
" method='nelder-mead' # Good choice for low-dimensional problems with no gradients\n",
" )\n",
"res_nograd"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "87af776a",
"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>True</th>\n",
" <th>Nograd</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>δ0</th>\n",
" <td>3.0</td>\n",
" <td>2.960324</td>\n",
" </tr>\n",
" <tr>\n",
" <th>δ1</th>\n",
" <td>-5.0</td>\n",
" <td>-5.311618</td>\n",
" </tr>\n",
" <tr>\n",
" <th>δ2</th>\n",
" <td>10.0</td>\n",
" <td>10.759441</td>\n",
" </tr>\n",
" <tr>\n",
" <th>σ</th>\n",
" <td>0.3</td>\n",
" <td>0.305551</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" True Nograd\n",
"δ0 3.0 2.960324\n",
"δ1 -5.0 -5.311618\n",
"δ2 10.0 10.759441\n",
"σ 0.3 0.305551"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.DataFrame(np.c_[np.r_[true_params, true_sigma], res_nograd.x],\n",
" columns=['True', 'Nograd'],\n",
" index=['δ0', 'δ1', 'δ2', 'σ'])"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "ec1d95c8",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAACIZElEQVR4nO3dd3hTZfsH8G+aNOlOd5sOOiiU2TLKKAKylzJEfVERxfUTcaOi4OtAfMEBiqCAg6ECilhQBCy77FUou3u3dO+ZeX5/1MSm2WlmuT/XlUtzznPOeRqR3H3GfbMYhmFACCGEEGInHKzdAUIIIYQQQ1DwQgghhBC7QsELIYQQQuwKBS+EEEIIsSsUvBBCCCHErlDwQgghhBC7QsELIYQQQuwKBS+EEEIIsSsca3fA1GQyGe7cuQN3d3ewWCxrd4cQQgghemAYBg0NDQgKCoKDg/axlS4XvNy5cwehoaHW7gYhhBBCjFBYWIiQkBCtbbpc8OLu7g6g7Yf38PCwcm8IIYQQoo/6+nqEhoYqvse16XLBi3yqyMPDg4IXQgghxM7os+SDFuwSQgghxK5Q8EIIIYQQu0LBCyGEEELsilnXvGzYsAEbNmxAXl4eAKBv3754//33MXXqVLXtk5KSMHbsWJXjqamp6NWrl8n6xTAMJBIJpFKpye5JiL7YbDY4HA5t5SeEECOZNXgJCQnBJ598gqioKADAjz/+iJkzZyIlJQV9+/bVeF16errSYls/Pz+T9UkkEqGkpATNzc0muychhnJxcYFAIACXy7V2VwghxO6YNXiZPn260vv//e9/2LBhA86fP681ePH394enp6fJ+yOTyZCbmws2m42goCBwuVz67ZdYFMMwEIlEqKioQG5uLnr06KEzGRMhhBBlFtsqLZVKsWvXLjQ1NSE+Pl5r24EDB6K1tRV9+vTBf//7X7VTSXJCoRBCoVDxvr6+XmNbkUgEmUyG0NBQuLi4GP5DEGICzs7OcHR0RH5+PkQiEZycnKzdJUIIsStm/5Xvxo0bcHNzA4/Hw4IFC7Bnzx706dNHbVuBQIDvvvsOCQkJ2L17N6KjozF+/HicPHlS4/1XrlwJPp+veOmTXZd+0yXWRn8GCSHEeCyGYRhzPkAkEqGgoAC1tbVISEjADz/8gBMnTmgMYDqaPn06WCwW9u7dq/a8upGX0NBQ1NXVqSSpa21tRW5uLiIiIui3XWJV9GeREEKU1dfXg8/nq/3+7sjs00ZcLlexYDcuLg6XLl3CV199hW+//Vav64cPH45t27ZpPM/j8cDj8UzSV0IIIYTYPouPXTMMozRSoktKSgoEAoEZe3T32rp1q1kWRpuTPfaZEEKIaZk1eFm6dClOnTqFvLw83LhxA++++y6SkpIwd+5cAMCSJUvwxBNPKNqvWbMGf/zxBzIzM3Hr1i0sWbIECQkJeOmll8zZTbswf/58sFgsldeUKVP0uj48PBxr1qxROjZnzhxkZGSYobfKLB1wyD+b8+fPKx0XCoXw8fEBi8VCUlKSxfpDCCHEtMw6bVRWVoZ58+ahpKQEfD4fMTExSExMxMSJEwEAJSUlKCgoULQXiUR48803UVxcDGdnZ/Tt2xf79+/HtGnTzNlNuzFlyhRs2bJF6VhnpsycnZ3h7Ozc2W7ZpNDQUGzZsgXDhw9XHNuzZw/c3NxQXV1txZ4RQgjpNKaLqaurYwAwdXV1KudaWlqY27dvMy0tLYpjMpmMaRKKrfKSyWR6/1xPPvkkM3PmTK1tPvjgAyY0NJThcrmMQCBgXn75ZYZhGObee+9lACi9GIZhtmzZwvD5fKXrY2NjmU2bNjGhoaGMq6srs2DBAkYikTCffvopExAQwPj5+TEff/yx0nNXr17N9OvXj3FxcWFCQkKYF154gWloaGAYhmGOHz+u8uwPPviAYRiGEQqFzFtvvcUEBQUxLi4uzNChQ5njx48r3XvLli1MaGgo4+zszMyaNYtZtWqVUp/VAcD897//ZTw8PJjm5mbF8YkTJzLvvfceA0DpOUVFRcx//vMfxtPTk/H29mZmzJjB5ObmKs5fvHiRmTBhAuPj48N4eHgwo0ePZi5fvqzyzO+//56ZNWsW4+zszERFRTF//vmnxj6q+7NICLG+06dPM7t377Z2N+5K2r6/O7JYnhdb1SKWos/7B63y7NsfTYYL1zT/CX7//Xd8+eWX+PXXX9G3b1+Ulpbi2rVrAIDdu3cjNjYW//d//4fnnntO632ys7Px999/IzExEdnZ2XjooYeQm5uLnj174sSJEzh79iyefvppjB8/XjGq4eDggLVr1yI8PBy5ublYuHAhFi9ejPXr12PEiBFYs2YN3n//faSnpwMA3NzcAABPPfUU8vLy8OuvvyIoKAh79uzBlClTcOPGDfTo0QMXLlzA008/jRUrVmD27NlITEzEBx98oNfnMXjwYERERCAhIQGPP/44CgsLcfLkSXzzzTdYvny5ol1zczPGjh2LUaNG4eTJk+BwOPj4448xZcoUXL9+HVwuFw0NDXjyySexdu1aAMDq1asxbdo0ZGZmwt3dXXGvZcuW4bPPPsPnn3+OdevWYe7cucjPz4e3t7ee/xUJIdZ2/vx5VFRU4IEHHrB2V4gWd33wYk/27dun+OKXe/vtt/Hee++hoKAAgYGBmDBhAhwdHdGtWzcMHToUAODt7Q02mw13d3cEBgZqfYZMJsPmzZvh7u6OPn36YOzYsUhPT8eBAwfg4OCA6OhofPrpp0hKSlIEL6+99pri+oiICCxfvhwvvPAC1q9fDy6XCz6fDxaLpfTs7Oxs/PLLLygqKkJQUBAA4M0330RiYiK2bNmCFStW4KuvvsLkyZPxzjvvAAB69uyJs2fPIjExUa/P66mnnsLmzZvx+OOPY8uWLZg2bZpKqYlff/0VDg4O+OGHHxTZlrds2QJPT08kJSVh0qRJGDdunNI13377Lby8vHDixAncf//9iuPz58/Ho48+CgBYsWIF1q1bh4sXL+q9LokQYn2ZmZmora21djeIDnd98OLsyMbtjyZb7dmGGDt2LDZs2KB0TP5b/cMPP4w1a9YgMjISU6ZMwbRp0zB9+nRwOIb9Jw4PD1caTQgICACbzVZKqhYQEIDy8nLF++PHj2PFihW4ffs26uvrIZFI0NraiqamJri6uqp9zpUrV8AwDHr27Kl0XL6oFmgryNnxt5/4+Hi9g5fHH38c77zzDnJycrB161bFyEl7ly9fRlZWltLPDLTlYcnOzgYAlJeX4/3338exY8dQVlYGqVSK5uZmpfVaABATE6P4d1dXV7i7uyt9ToQQ21dUVATGvOnPiAnc9cELi8Uy2dSNubm6uipy5nQUGhqK9PR0HD58GEeOHMHChQvx+eef48SJE3B0dNT7GR3bslgstcdkMhkAID8/H9OmTcOCBQuwfPlyeHt74/Tp03jmmWcgFos1Pkcmk4HNZuPy5ctgs5WDOPnoUmf/AvHx8cH999+PZ555Bq2trZg6dSoaGhpU+jF48GBs375d5Xr5KM38+fNRUVGBNWvWICwsDDweD/Hx8RCJRErttX1OhBD7wDAMZcC2A/bxrU304uzsjBkzZmDGjBl48cUX0atXL9y4cQODBg0Cl8uFVCo1+TOTk5MhkUiwevVqxf/wv/32m1Ibdc8eOHAgpFIpysvLMWrUKLX37tOnj8p2547vdXn66acxbdo0vP322ypBEgAMGjQIO3fuhL+/v8aMjqdOncL69esVu94KCwtRWVlpUD8IIbavtbUVTk5OcHBwQHNzM9XAs2EUXtoRoVCI0tJSpZf8S3Tr1q3YtGkTbt68iZycHPz8889wdnZGWFgYgLbpoJMnT6K4uNikX7zdu3eHRCLBunXrFM/duHGjUpvw8HA0Njbi6NGjqKysRHNzM3r27Im5c+fiiSeewO7du5Gbm4tLly7h008/xYEDBwAAr7zyChITE/HZZ58hIyMDX3/9td5TRnJTpkxBRUUFPvroI7Xn586dC19fX8ycOROnTp1Cbm4uTpw4gVdffRVFRUUAgKioKPz8889ITU3FhQsXMHfu3C67xZyQu1l2dja6d++O8PBw5OfnW7s7RAsKXuxIYmIiBAKB0mvkyJEAAE9PT3z//fe45557EBMTg6NHj+Kvv/5SrB/56KOPkJeXh+7du6ssWu2MAQMG4IsvvsCnn36Kfv36Yfv27Vi5cqVSmxEjRmDBggWYM2cO/Pz88NlnnwFoWxj7xBNP4I033kB0dDRmzJiBCxcuKIprDh8+HD/88APWrVuHAQMG4NChQ/jvf/9rUP9YLBZ8fX3B5XLVnndxccHJkyfRrVs3zJ49G71798bTTz+NlpYWxUjM5s2bUVNTg4EDB2LevHl45ZVX4O/vb+hHRQixcVlZWYiKikJ4eDjy8vKs3R2ihdkLM1qatsJOVAyP2Ar6s0iI7Vm1ahUGDx6MlpYW5Ofn44UXXrB2l+4qhhRmpJEXQgghBG3bpHv06EEjL3aAFuwSQgghAIqLixEUFAQvLy8KXmwcjbwQQggh+HebtKurK5qbm63dHaIFBS+EEELuei0tLUq7CLvYctAuh4IXQgghZvHTTz9Zuwt6k2+TlnN1dUVTU5MVe0S0oeCFEEKIWbz77rt2k2U6MzNTKYM55XqxbRS8EEIIMZhQKNR6XiwWo6ioCHfu3LFQjzonKysLPXr0ULynHUe2jYIXQgghButYNLWjiooKsFgsRYFTWyffJi1HwYtto+CFEEKIQcRiMdLS0rS2KS0tRd++fZGVlWWhXnXOnTt3EBQUpHhPwYtto+DFTsyfPx8sFktR5TkgIAATJ07E5s2bDZpT3rp1Kzw9Pc3XUUJIl9fQ0IDS0lKtO3LKysowcuRIuxl5YRgGLBZL8T4sLIyCFxtGwYsdmTJlCkpKSpCXl4e///4bY8eOxauvvor7778fEonE2t0jhNwlGhoa0NLSgvr6eo1tSktLMWLECLsYeVFXQdrFxYVyvdgwCl7sCI/HQ2BgIIKDgzFo0CAsXboUf/75J/7++29s3boVAPDFF1+gf//+cHV1RWhoKBYuXIjGxkYAQFJSEp566inU1dUpRnE+/PBDAMC2bdsQFxcHd3d3BAYG4rHHHkN5ebmVflJCiC1raGgAAJSUlGhsU1ZWhsjISLsIADpukya2j4IXOzdu3DjExsZi9+7dAAAHBwesXbsWN2/exI8//ohjx45h8eLFANqqO69ZswYeHh4oKSlBSUkJ3nzzTQCASCTC8uXLce3aNfzxxx/Izc3F/PnzrfVjEUJsWH19PXx9fVFaWqqxTVlZGQICAgDYfsK3jot15dzc3BS//BHbQrWNuoBevXrh+vXrAIDXXntNcTwiIgLLly/HCy+8gPXr14PL5YLP54PFYiEwMFDpHk8//bTi3yMjI7F27VoMHToUjY2NcHNzs8jPQQixDw0NDejZs6fWkZfS0lIEBATAz88PlZWV8PPzs2APDZOVlYVhw4apHJfneunbt68VekW0oeDlHy+88AKKi4st9rzg4GBs2LDBJPdqv9Ds+PHjWLFiBW7fvo36+npIJBK0traiqakJrq6uGu+RkpKCDz/8EFevXkV1dbViEXBBQQH69Oljkn4SQroGfYIX+S8+3bt3R3Z2tk0HL5mZmXj88cdVjst3HFHwYnsoePmHqQIJa0hNTUVERATy8/Mxbdo0LFiwAMuXL4e3tzdOnz6NZ555BmKxWOP1TU1NmDRpEiZNmoRt27bBz88PBQUFmDx5MkQikQV/EkKIPWhoaECPHj20ThsBAIvFQlRUFLKysjB8+HAL9c5wJSUlEAgEKsfDw8PVLjgWiUSoqKhAcHCwJbpH1KA1L3bu2LFjuHHjBh588EEkJydDIpFg9erVGD58OHr27KmS3ZLL5UIqlSodS0tLQ2VlJT755BOMGjUKvXr1osW6hBCN9Bl5kZOPvNiyjtuk5TTlevnpp5/wyiuvWKBnRBMKXuyIUChEaWkpiouLceXKFaxYsQIzZ87E/fffjyeeeALdu3eHRCLBunXrkJOTg59//hkbN25Uukd4eDgaGxtx9OhRVFZWorm5Gd26dQOXy1Vct3fvXixfvtxKPyUhxNbV19cjJCQEdXV1as+LRCI4OjoCgGLkxVap2yYtFxYWpra+UUJCAhobG6lwoxVR8GJHEhMTIRAIEB4ejilTpuD48eNYu3Yt/vzzT7DZbAwYMABffPEFPv30U/Tr1w/bt2/HypUrle4xYsQILFiwAHPmzIGfnx8+++wz+Pn5YevWrdi1axf69OmDTz75BKtWrbLST0kIsXUNDQ3w8PDQeL68vFyx08jLywu1tbUW6pnhsrKylAoytufs7IzW1lalYzdv3kRUVBQeeOABJCYmdvr5Z86c6fQ97ka05sVObN26VZHLRZvXX38dr7/+utKxefPmKb3fsGGDyhqfRx99FI8++qjSMVvf3kgIsY6Ghga4u7trPN9+mzRg23+XaNomrcmmTZvw7LPPIiAgAG+99RYefPBBo58tFArxyiuv4PLly0bfQ+7FF1/EM888g0GDBnX6XvaARl4IIYQYRB68cLlctdWl5duk5Tw8PLRm47UmbSMvQFuuF3lSPqFQiNTUVMTGxiIwMBA1NTU6q2trU1FRgYKCAqOvby81NRVLly7FhQsXTHI/W0fBCyGEEIPIt0EHBgairKxM5XxZWZlSLilbXrSra+Sl/aLdP/74A7NmzVKcmzhxIg4fPmz0s8vLyxVrD7XRtlsUaAuqPD098dtvv2HZsmU4deqU0X2yFxS8EEIIMYhEIgGHw0FgYKDaHUcdR15sedFuaWmpStLO9toHLzt27FCaXp89ezYSEhKMfnZ5eTlcXV1RWFiotd3kyZO1ns/Pz0d4eDg8PDzw22+/4bPPPsPRo0eN7pc9oOCFEEKIUQQCgdrgxZCRl9TUVLP1T1/qtknLyYOX3NxceHt7g8/nK86FhoaitLRU58iIJuXl5Rg0aJDWqaOWlhacPHlS6/RUbm4uIiIiALRNc/3222/4+uuvceLECaP6ZQ8oeCGEEGIQ+Ze9QCBQm6iu44JdTSMvEokEY8aMgUQiMV9ntZCPIGkjD142b96MZ555RuX82LFjcfz4caOeX1ZWhsGDB2sNXgoKCiCVSpGTk6OxTU5OjiJ4Adp2SX366afYu3evUf2yB2YNXjZs2ICYmBh4eHjAw8MD8fHx+Pvvv7Vec+LECQwePBhOTk6IjIxUyVNiCra88p3cHejPIOkKNE0bdSxHEhAQoHZtTFJSEqqrq1FZWWnWfmqiKbNue2FhYcjJycG5c+dwzz33qJx/8MEHjZ46Ki8vR1xcnNZpo9zcXPTr10/rtFtubi4iIyOVjvn5+XXpZKNmDV5CQkLwySefIDk5GcnJyRg3bhxmzpyJW7duqW2fm5uLadOmYdSoUUhJScHSpUvxyiuvdGpOsT150iR7KNFOujb5n0H5n0lC7JGmaaOONE3L/Pbbb3jooYes9iVbWFiI0NBQrW2cnJyQkpKCKVOmqP05unfvjry8PJXM5RcvXsRPP/2k9d7y4EXbyEteXh4mTpyoM3gJDw9XOubp6WnT+XU6y6x5XqZPn670/n//+x82bNiA8+fPqy10tXHjRnTr1g1r1qwBAPTu3RvJyclYtWpVp/bSy7HZbHh6eir+R3FxcdE610mIqTEMg+bmZpSXl8PT0xNsNtvaXSLEIPKirQDg7++PiooKva7j8XhoaWmBs7MzgLYdNAUFBZg6dapNBy9A23fRE088ofH8yJEjcfr0adx7771gGAZffvklLly4gNbWVq3X1dTUoEePHlo/w9zcXEycOBH79u3T2KalpUUlS3BX/26zWJI6qVSKXbt2oampCfHx8WrbnDt3DpMmTVI6NnnyZGzatAlisdgkv6XKF5F15eE0Yvs8PT217nAgxFY1NTXBzc0NQNsvhB1HHIRCIbhcrsp1kZGRyM3NVVSpP3bsGCZMmAB/f3+dfx9nZGQgKioKDg6mnSwoKCjAkCFDdLbbv3+/1mc/+OCD2LhxI/r374/nn38eo0aNwq+//ooZM2ZovS/DMDp/pry8PIwcOdKmigdnZ2fDzc1NaV2TpZk9eLlx4wbi4+PR2toKNzc37NmzR/GHt6OO2+uAtrlSiUSCyspKtXOTQqFQaRW2rkRILBYLAoEA/v7+Rq8QJ6QzHB0dacSF2C1d2XXblwZoT75oV/73/2+//Yb3338fWVlZuHHjhtZnLlmyBJ988onWfCynTp1CYGCgQdlyCwsLMXv2bJ3tdAUYffr0wfXr1/HQQw/hs88+Q1xcnN59ANq+lzQVh2xqaoK7u7vSiFd7tbW1SjugOtJ0385Ys2YNXn755a4dvERHR+Pq1auora1FQkICnnzySZw4cUJjANPxQ5YvbNT04a9cuRLLli0zuF9sNpu+QAghxEAdgxcWiwWZTKb4gteUN6V79+64fv06gLbCjSUlJQgLC0N9fb3OkZeSkhLk5ORoDUy2bNmC9PR0/PDDD+jdu7deP0tRURFCQkL0aqvLp59+il69eikFEhwOR+OOpvaL9n19fVFZWQk/Pz+N9+dyuRCJRCqjWuoW68q5ubmhsbFRa7BpKIZhkJGRYVCQaA5m3yrN5XIRFRWFuLg4rFy5ErGxsfjqq6/Utg0MDFTZdldeXg4OhwMfHx+11yxZsgR1dXWKl65kP4QQQozXMXjx8fFBVVWV4n3HbdJy7bdLHzlyBBMnTgQAvaaNSktLtW4VBtq+K3bv3o2XX34ZN2/e1OtnEYlE4PF4erXVZdiwYSojIJoyEANtswTy9t26dVO7aLexsVGxaysiIkKRLK+99jleOvLz89N7TZK+MjIyEB0dbfU1NRbP88IwjMZkO/Hx8Sqplg8dOoS4uDiN6114PJ5iK7b8RQgh9sxaeU/00TF46bjjSN30P9C2+7SoqAgAsGvXLjz00EMA2oIfbVul5etCdAUvQNsyg19++QWvv/46UlJS9P6ZzEUgEODOnTtqz5WXl8Pf3x+A5uAlLy9PEZhERUUhMzNTpU1OTo7GkRdDFlTr68CBA5g2bZpJ72kMs04bLV26FFOnTkVoaCgaGhrw66+/IikpSVFGfMmSJSguLlZsJ1uwYAG+/vprLFq0CM899xzOnTuHTZs24ZdffjFnNwkhxGY0NjbioYceUvw9aWvq6+tVgpfS0lLExMQAaBt56devn8p1bDYbMpkMQqEQ5eXlil0+HA5HZdFvew0NDYiNjUVubq7GNq2trXBycgLQNtqwc+dOPPLII/jf//6ncUGupoXFphQUFKRX8BIaGqo203BeXp5iC3RUVBRu376t0iY3NxcPP/yw0jGGYSCUyODs4YnbOUXwDOuDZpEELSIpWsT/vERStEpkaBVJ0SqWolUiRatYhlaxFELJv/8USmQQiqUQSWUQSWS4uCkB0Y++h89uJeHYG2M69wF1glmDl7KyMsybNw8lJSXg8/mIiYlBYmKiYriwpKREKdqMiIjAgQMH8Prrr+Obb75BUFAQ1q5da5Jt0oQQYg9OnTpls0UMgbZgov0Id8dEdZpGXoC2QGX//v06a/W0V1JSgoiICKSlpWls03HqxNvbG7t27cKjjz6KAwcOqL3GlOtdNAkKCkJ+fr7acx1HXg4ePAgAkMkYNLRKUNcixpmU2/ANDsf+6yXIbnLGHycuQ9I7HQ2tEtS3itHYKsGBMzdQ+GceGsW5aBZK0SSUoEkkgYwBmtNLIBNmwu2ak0l+HpmwGc1iBkUNEjg2aw44LcGswcumTZu0nt+6davKsXvvvRdXrlwxU48IIcS2HTlyxOrrCbRpaGhAWFiY4r1AIFCqYqxpzQvQ9iW9evVq7Nq1S+/nyRcAp6Wladw5k52drTJ1wufz4eDgoPEafXO8dEZQUBBOnzmL0rpWVDQIUdnY9qpuEuHvY9fAcvXBma2XUF5bjwuHk3H0o0OoaxFD9s9a3uqjZ+E+gA/HjCtgZFJU3MhA7jHlZHXVja24fqdR7fMdXDzA1JfA25ULZ0c2XLhtLydHNpy5bDg7tv27k6MDeJx//93JkQ0ep+0Yj+MALscBPI4DLiYlom7uLDzyZDy4Vt7wYrE8L4QQQnS7desW+vfvD6FQaLLFpKakbs1L+40W6hKmyUVFReHq1asICgpSOu7i4qJSUkCutLQUAoEAPj4+qK6uVrt5Izs7W+0OI/mCWXW7nwoLC9GtWzfNP6gODMOgplmMO7UtKKlrRUldC8rqW1FWL0RZfSvK64UoLilFXuJF7JCqVniuvZoF5wgP8NC2WLmxqRk1zf+m73B2ZIPTUoXePSLh7eEGdycOjh/h4vH4MLg7ceDGc4Qr1wFfn/PEiifi4MrjwN2JAxcuG648Dlx5HBTkdMfmTT9g1XsTjf452/v9q7NYunQpwsO8TXK/zqDghRBCbERFRQV8fX3h5uaGyspKBAcHW7tLKjoGL5rqG6kzfPhwtaMy8oWl6oKXkpIS9O/fH5GRkcjJydEYvNx///0qx3v16oW0tDSNwcuYMWO09reuWYy8qibkVzejqKYZRTUt/7yacae2Ba1i9blX5BiWM2Qt9WA7sODjyoWvGw8+bm3/PHeLwbSpg9CjeyS8XLj4+JIXNr8+Gp7OjuC7OILHYeP+lK+xb/G/gcesP7zx3rRoxQaWO3fu4EyvCEzoo36kKzDAdAt2GYZBfn6+ShkCa6HghRBCbMTx48cxbtw45ObmoqKiwiaDl44Ldl1cXPSuFzdkyBC1C2jl26XVfTGWlpZi4sSJiuBF3fUFBQVqR1HkwYu6IEU+bdQikiK3sgk5lY3ILm/7Z15lE/KqmlHXojuRqa8bD0GeTgj0cIKA7wR/DycEeDghwIMHf3cnvHbzW+z/eCocHJSnrubsluDNmUMVo1Tf8V0R7sXTuog4PDwc+fn5iIqKAqBaTbojT09P1NTU6PwZ9HHjxg3FomxbQMELIYTYiKNHj2Lp0qX4888/Tb7F1VQ6Lthtr7W11aipLm25XuRrXiIjI3Hy5Em1bSQSidp0Gr169cKRI0cAAGKpDDkVTUgrrUdGWQMOXEjFre3pKKy7Dm1F3v3deQjzcUGotwtCvFwQ4uXc9vJ0QQCfBx5H+9oPjgNLJXABVKfXQkJCUFxcrAhGamtr4enpqXSNPFeOPHjRluMFgGLNjynYyhZpOQpeCCHERuTn5yMsLMwsycVMRV15APmCWE3rS3TRFrxUV1fD29sbkZGRajd5yGQylQW5rWIpUkvqcbUE2H/mGrLXnUZ6WQNEkn+neSrqm8Gqbcs55uniiO5+boj0dUWknxsifF0R7uuCbt4ucOF27mvS0dFRbWbcjuS5XuTBiLpK0T169EBGRgamTJkCoG3kxVIBxenTp7Fo0SKLPEsfFLwQQogNyMvLU+zi8ff31ztLrKWpW0js6uqKxsZGrTuNtPH390d6errac/LdQprW1hQVFcHDxx+/Xy7C1cIaXC+qQ2pJPcTSthGH8pomiIvrAABuPA6iA90RHeiOA6fcsf65YegZ4A5fN/MtjJYvaO44rdVxRCQ0NFQpQ3z7BHVyUVFRSlu/dY28mEpNTQ1cXV3NnhfHEBS8EEKIDTh69CjGjx8PoC3Rmq6U+dbUcaRDXtpFU10jXfQpESB/plTG4NadOlzIqUZyfjWSkpJQns/gzK5rSu19XLnoH8JHisAbH8zuhcHdAxHq5QIHBxYaGxuR/asvRnT3NbivhpInqmsfvEgkEpXaet26dcPx48cV73NzcxEbG6vUJiwsTClvTFVVlcbSOXKurq4ad3Lp69ChQwbl5rEECl4IIcQGHD9+HF9++SUA89SkMRV1ayjkJQI6M/KiLniRFzVMK63H2awq3LjTiNgPDqCx3TraxpJCOHsHIS7MCwO7eWJAqBdiQvgI8XIGi8XCf7PjEMltQJhPd8U1lsjxIqcuy666IowdSwTk5eVh1qxZSm3khR7b05UTSP5nqTPBy99//40VK1YYfb05UPBCCCFWxjAMqqqqFF9o8irD9kIevLQvE2AIDw8P1NXVKd7XNYtxOqsS+87dxJliMaasaUuCV+PgAUl5CbwDQzE0whtDIrxxseEoFjw1G0PjBqm9t3zH0YABAxTHLB28yAtSyrXPrivXcVqsoKBAbR/ZbDYkEgmkUqle0zjyUTxjtzjn5OSgtLRUJTePtVHwQgghVnbr1i307dtX8d7R0dFmizOq+00/MDAQqampRi/YZbFYaBFJ8cOpHBy6VYbk/GrIGEBYmgWRowe8HR0wNMIHrU2xGB/vg2cfmQT2Pzt4kjbcQe/oHhrv3atXL+zbt0/pmKWDl467pNQFL/LaT3KadlCFhYWhoKAAUqlUKdOxJsYWZ2QYBt9++y327duHDRs2GHy9uVHwQgghVtZ+vYs9EggEOHbsmMHTRumlDdh//Q4O3S7D2exKZO3/tzhhlL8bgpx48IsYgP8tnQQnRzb++qsMd+7cUQQuQFshy467n9qLjo7GqlWrlI4VFhZaLGeJumkjdcFLe9q2N8u3SzMMo9diXWOmIAsLC/Hyyy9j7Nix2Lt3LxwcHAy63hIoeCGEECs7efIknnnmGWt3QydNlZjl00btqztrklPRiH3XS/DXtTvILP+3Jg+LxUJ8hBcm9RNgQu8AhHq7YPPmbHh794STY9vi1sjISJw+fdqgPru7u6OxUbn2jyVHXry9vVFVVaV0rLy8HP3791dpy+fzUVdXB5FIBF9f9YuJ5cGLg4ODSj0ndfz8/JCTk6N3f/fv349169bhq6++QnR0tN7XWRoFL4QQYkUSiQStra1wc3NTOs7hcCAWi9VOHViLuhwvQNsXdHV1tcbraptF2HvtDn6/XITrRf+ubeGyHTC6px+m9Q/ErsKe+ObhXvD2/rduTklJCfr06aN4HxERofRFXFNTo5LITR0HBwfIZDLFCIKukQ9TYrFYKiMpmp4vX7Tb2tqqcVQlKioKhw8fBpfLxejRo3U+35Cda83Nzfjiiy+QmJhoU3/u1KHghRBCrOjSpUtqU97LF+0KBAIr9Eq9+vp6tdl11a2DkcoYnMyswO/JRTh8uwwiadt6DrYDCyOjfHF/jACT+gaC79z2JXkuKBDl5eVKwUvHrdcuLi5oaWlRvM/JyUH37v/uItJEvk5EvmiVYRirToXoCl6ampo0LrANDw9HXl4eHB0d9VqEa8ial3Xr1uHFF1+0+cAFoOCFEEKsSl7PqCP5WgVbCl40jbwAbenuvby8UNkoxM5LhdhxoQDFtf8GGr0FHnh4cAhmDgiCj5qkcPLt0r169VIc05U3Jjs7W6/gRb7jKDw83GTp8g3B4/GUptQqKyvV5meRJ6qrr6/H8OHD1d6Ly+VCLBZDIpForN7dnpeXl171jWpqanD8+HEsXrxYZ1tbYHurcAgh5C5y+fJlxMXFqRy3xVwv2oIXERyRUeeAESuP4fOD6SiubYGniyPmjwjH/ldG4u9XR+HpkRFqAxdAfa4XdWto2n8ZGxq8AG01g7y8vHReY0pBQUEoLS1VvJdKpeBwVMcO5CMv6koDtGdIzSJ923766adYvHixzrwxtoKCF0IIsRKGYdDc3Kz2N2h7CF6kMgaJN0sxe/0ZXKlkkN/ChUgqQ2wIH6sejsX5JePx4Yy+6BvE13lvfbLsAlBUlwbaghd9Fq22D14suVhXTt2OI3XkwUtxcbHWiuKhoaF6rfWR0xW8FBcXIy0tTe0IoK2iaSNCCLGS4uJihISEqD3n7++v+MK1FfLgpVUsxe+Xi7DpdC5yK5sAAFx3HwzuHYHPX7wHsaGeBt/b398fFy5c0NlOHrwMHjxY77wyQUFBKC4uBmDbwYt8txHDMCrlA9qLiorSukDaUB9//DHee+89k93PEih4IYQQK9G0WBdoG3npmNzM2qpqanG7Glj72XFUNLRVZPZw4uDx4WEIGjQbPcNDjApcANWRF035WyIjI3H27FnFe32mOdq3sUbwIhAIkJralsOmqalJ61oVeSFKbe69917U1tbq/XwXFxeNI3wZGRloaGjA4MGD9b6fLaDghRBCrOTSpUt4+OGH1Z6zpWmjVrEU2y8UYMWfKRDxQ+EcFowgvhOeGx2J/8SFwpXHAdBL53206bilt6SkRO2oSmRkJLZt26Yx54wmnp6eqKmpQWFhIeLj4zvVV0MFBQXh6NGjAICKigqt27RlMpnWKSMAGDRIfSkETeR/ltRl5F22bBk+/PBDg+5nCyh4IYQQK7l58yaWLVum9pwt1DeSSGXYmVyINUcyUdEgRENDA7pFeOKdB/rjocEh4HJMt2zSyckJQqFQ8V7TTiOBQIA7d+4gLy/PoHo9vXr1Qnp6utWnjcrLy7VmIe7WrZvJ+6cpeCkqKgKbzUaPHprLK9gqCl4IIcQKGIbRWL8GaNteKxKJLNyrNgzD4OCtMnx2MA05FW1rWoI9ndEvzBUrXhiDmH7dzN4HTcGLfPeMvot15eSLdqurqy2+20g+6gPoTpDXrVs3o4soaqJpFO/ixYt6JbqzRbTbiBBCrCArK8smf+NNKajBQxvPYcG2y8ipaIK3KxcfTu+D42+OQYgbCz5euncOmUJJSYnGHDdsNhsZGRl6bZOWa7/jyNLbgVksluKZZWVlWoOXp556Cvfff79Jn69pJ1dycrLabfr2gIIXQgixAm2Lda2hslGIt3ZdwwPrz+Jyfg2cHB3w0tgonHhrDObfEwEux0Fjhl1TcXR0VIw2aUtQ161bNyQlJRkUvERFRSEjI8PqeUx0jbwIBAKTf8aaRl5u3rypVM3cntC0ESGEWMGlS5ewYMECrW3YbDYkEonahGamIpHK8PP5fHxxOAMNrRIAwEODQ/DW5GgEeCgniGtqaoKrq6vZ+uLn54fKykpFUjdNwUtkZCR27NihdgGqJjweD3fu3FGqlWRJzs7OaG5utmhdJTl1wYuuaUtbR8ELIYRYQUZGhs5pIx8fH1RVVWld4NkZVwpqsHT3DaSVNgAA+gV74KOZ/TCom/o1IeauCSSf3ggKCkJVVZXaFPpAW/Di4eFh0G4joO1L3NKLdeWCgoJQUlJiM8FLbm6uQWuGbA0FL4QQYmESiQQODg46AwH5l46pg5cmoQSrDqVj69k8MAzg6eKIxZN7Yc6QULAdrDet0n5tRvsq0B1FRkYaNGUk16tXL6sFL/JdUvX19RpLLJiLuqrf9rzeBaDghRBCLO7WrVt6rTUwda6XoqIifLtzH46JeyqKJs4eFIz37usDL1fDRjHMQd8SAT169MBLL71k8P1nzJgBX19fY7rWae23S1t63Q2bzYZMJlM6lpycjCeeeMKi/TAlCl4IIcTC9F2sa8rgpVEowQurt+Pw7u0IfHQlgj2dsWJ2f9zb00/ve5i7IrO/vz9SUlIglUq1psd3cnLCrFmzDL7/qFGjOtG7zgkKCsK1a9es9vyO0tLSlCp42xvabUQIIRamb/Di7+9vkuAlOa8aU786iZPnL0NaX4l5Q4Jx6PXRBgUuliAfeamsrISfn231rbPk9ZWsvdsJaJuSk8lkZl0Ibm4UvBBCiIXpm+W1Y8p8Q4kkMnx+MA3/+fYcCqtbwKm/g7mPPoL7glv/SemvP12jIaYgD140lQawZ0FBQbh16xa8vb2t8nwej4fW1lYAbdW4o6KirNIPU6HghRBCzGT16tU4ffq00rHW1lbweDy9fgPvzLRRbmUTHtxwFt8cz4aMaVvbEtfNA088dD9OnTpl8P0aGxvh5uZmVF/05e3tjaqqKq3bpO2Vu7s7srKyLL7TSK79KJ69L9YFKHghhBCzOXv2LJYtW4Zbt24pjl29ehUDBgzQ63pjg5f910swfd1p3CiuA9/ZEd88Ngjvjg+Fv68Phg0bhvPnzxt8z4aGBrPvkpEvLC0tLdWYXddeybPsWit4aT+KR8ELIYQQjcRiMbZv346XX34ZhYWFAAzLrNuxWKEuIokMH+69hRd3XEGjUIIh4V44+Npo3BcjwPXr1xEbGwtXV1e0tLSo7D7RxRLBi1xXnDYC2rZLWzN4kQfC6enp6Nmzp1X6YSpmDV5WrlyJIUOGwN3dHf7+/pg1axbS09O1XpOUlKSIUNu/5DUpCCHEnvj7++P777/H/PnzUVNTY7ayAEU1zXj423PYejYPALDg3u745bnhCOS3Zcm9du0aYmJiAAB9+vRBamqqQfc3d2mA9rpq8BIUFGT14EUqlQKAWZMNWoJZe3/ixAm8+OKLOH/+PA4fPgyJRIJJkyahqalJ57Xp6ekoKSlRvGyxgBkhhOije/fu+PzzzzF37lyUlJQYtJNGn+3JZ7MrMX3daVwrrAXf2RGbnozDO1N7gcP+9694+cgLAIwcOdLgdS+WGnlxc3NDdnZ2lwxeIiIiEBwcbJVny9e8ZGRkIDo62ip9MCWz7pNKTExUer9lyxb4+/vj8uXLOstw+/v7w9PT04y9I4QQyxk0aBAWLVqEvXv3GnSdg4ODxp0+DMNg2/l8fPjXbUhlDPoH87F+7iCEeruotC0rK1Nk6h05ciRee+01nbWV2rNU8OLv748bN27AxUX1Z7B3y5cvN/uOLU3ka166wnoXwMJrXurq6gBAr61iAwcOhEAgwPjx43H8+HGN7YRCIerr65VehBBibSKRSKXo3YQJE7B27VqD7uPj46OS2h1oW9/y7h838d6ftyCVMZg1IAi7FsSrDVwkEonSl6avry+qqqoM6oclgxexWGz251iDtQIX4N9pIwpeDMQwDBYtWoSRI0eiX79+GtsJBAJ89913SEhIwO7duxEdHY3x48fj5MmTatuvXLkSfD5f8bJW3QpCCGmvpqbGJDk91O04qm4S4fFNF7DjQgFYLOCdqb3w5ZwBcHJU/+WYnp6uMlUQGhqKgoICvfthyeClK04ZWZu8yGdWVpZRdaFsjcXS67300ku4fv26Ss6DjqKjo5X+J4uPj0dhYSFWrVqldqppyZIlWLRokeJ9fX09BTCEEKurrq42S/CSV9mE+VsuIq+qGe48DtY+OhBje2lfBNp+vYvcqFGjcOrUKcydO1evflhqwa6/v3+X2yZtC9hsNsRisV4FQe2BRX6Cl19+GXv37sXx48cREhJi8PXDhw9HZmam2nM8Hg8eHh5KL0IIsbbq6mp4eXl1+j7tg5eUghrM3nAWeVXNCPFyxp4XR+gMXIC2nUbqghddv0y2RyMv9q+oqAh9+vSxdjdMwqwjLwzD4OWXX8aePXuQlJSEiIgIo+6TkpJCkTghxK6YctooPz8fh2+X4eVfrqBVLEP/YD42zY+Dv7uTXvdIS0tTmTYKCwtDXl6e3v2wVPASFxeHbt26mf05dyNPT88usd4FMHPw8uKLL2LHjh34888/4e7ujtLSUgAAn8+Hs7MzgLZpn+LiYvz0008AgDVr1iA8PBx9+/aFSCTCtm3bkJCQgISEBHN2lRBCTMpU00b+/v74Yc9hfF6QDBkDjIn2wzePDTKoNpFYLAaXy1U6xmKxFOn4fXx8dN7DUsGLo6Oj1bYTd3WBgYFdJngx67TRhg0bUFdXhzFjxkAgECheO3fuVLQpKSlRWjQmEonw5ptvIiYmRjGsuX//fsyePducXSWEEJMyxbQRwzD4O7MJR65kQsYAjwwJxQ9PxBkUuFRWVsLX11ftuXvuuQdnz57V6z6WzLBLzGPz5s1Gz4DYGrNPG+mydetWpfeLFy/G4sWLzdQjQgixjM5OGzEMg08S0/B9ciWkzfV4ZVwUXp/YU6+Cju21z6zb0ciRI7Ft2zZMnz5d531EIpHK6A2xL66urtbugslYbLcRIYTcTTozbSSTMXjvz5vYfqEADo5OiPLmYtEk47KiqttpJNevXz/cvHlTr/sYGjQRYk72v1+KEEJskLHTRhKpDG/suobt/+RwWfFAf4T5GJ9tVt1OIzkHBwc4OzsjKyvL6PsTYg0UvBBCiBkYs0ZELJXhlV9TsCelGGwHFtbMGYDHhnXTawpek4qKCq21lFavXo0333wTn3/+OSQSidHPIcSSKHghhBAzMWSqRSKV4bWdV3HgRim4bAdsfHwwZg5o23Xj4OAAmUxm8PPFYjE4HO2rA8LDw7Fnzx74+/vjvvvuw9WrVw1+DiGWRsELIYRYmUQqw6LfrmH/9RI4slnY8PggTOwToDjv5eWFmpoag++bnp6OXr166WzHYrHw5JNP4qeffsKqVauwbt06pfOdGfkhxBwoeCGEECuSyhi8uesa9l67A44DC988NgjjewcotVFX30gf2ta7qBMQEIBt27bh4MGDqK2tVRxvaWlR5OYixBZQ8EIIISam7xSPTMbgrd+v4Y+rbYHL148NwqS+qqnxLRW8yL3wwgvYuHGj4j3leCG2hoIXQggxsYaGBp111hiGwQd7b2H3lbbFueseHYgp/dTX9ImIiEB6errB/cjIyEDPnj0Nvm7q1Kk4cuQIhEIhAApeiO2h4IUQQkxMnwR1XxzOwM/n88FiAV/OGYCp/TXXb5syZQoOHDhgcD+kUikcHR0Nvs7BwQGPP/44tm3bBoCCF2J7KHghhBAT05Wg7odTOVh3rC23yvKZ/TAjNkjr/dzd3cHhcFBdXa13H1paWsDj8fRu39Fjjz2GX375BTKZTK+RJEIsiYIXQggxMW0J6nYlF+Lj/akAgLcmR+Px4WF63XPmzJn4888/9e5DdnY2oqKi9G7fEZfLxaRJk7Bv3z7U19fTyAuxKRS8EEKIiWmaNjp4qxRvJ1wHADw3KgILx3TX+57Tp0/HX3/9pXf79PR0REcbV1JA7vnnn8e3335L00bE5lDwQgghJqZu2uhyfg1e+SUFMgb4T1wIlk7rbVASOw8PDzg4OOid78UUwQufz0efPn1w6NAhCl6ITaHghRBCTKzjtFFuZROe/fEShBIZxvfyx4oH+htV6HDmzJnYu3evXm2N3WnU0auvvoqdO3dS8EJsCgUvhBBiYu2njaoahZi/5SJqmsWICeFj3WMDwWEb91fvjBkz9A5eKisr4evra9Rz2gsJCcG8efPg7+/f6XsRYirai14QQggxmHzaqEUkxTM/JiO/qhkhXs7Y9OQQuHCN/2uXz+eDxWKhtrYWnp6eGtsxDGPSlP4bN240aqSIEHOhkRdCCDGx6upqePA98drOFFwtrAXf2RFbnxoKP3fjty7L6TP6UlVVpbWStKEocCG2hoIXQggxMYlEgi+PZuPgrTJw2Q74/ok4RPm7meTe+gQvplisS4gto+CFEEJM7E5tC749mQMA+PzhGAyN0J5t1xCenp5gGAZ1dXUa26Snp5tksS4htoqCF0IIMaHkvGqklTYAAF4eF4WZA4JN/ozp06drHX3JyMigkRfSpVHwQgghJlJU04znf74MGcNgct8AvD7BPKMfurLtZmVldSq7LiG2joIXQgjRQSqVorW1VWubJqEEz/6YjMq6Jni4OuPLOQPg4GCeha5eXl4Qi8Ua+9Ta2gonJyezPJsQW0DBCyGE6LBnzx6MGzcOLS0tas/LZAwW/XYVaaUN8GSLMC42olNbovUxatQonD59WuW4VCoFm80267MJsTYKXgghRIerV6/igQcewFNPPQWpVKpyfsOJf3cWvTepG0IDzZ/QbdKkSTh06JDK8by8PISHh5v9+YRYEwUvhBCiw+3bt/Hyyy/jvvvuwxtvvKGUAO5kRgVWHUoHACyb2RcCJ6naooym1r9/f1y/fl3luKnKAhBiyyh4IYQQHYRCIZycnDBv3jz4+fnhyy+/BAAUVjfjlV9TwDDAI0NC8ejQbip1jcyFxWIhICAApaWlSscpxwu5G1DwQgghWtTV1cHDw0PxfunSpcjIyMC2Hb/i+Z8vo7ZZjNhQTyyb2ReAcl0jc5s4cSKOHDmidIyCF3I3oOCFEEK0uH79OmJiYhTvWSwW1q1bhyX/W41bhVXwceViw9xB4HHaFsnK6xpZwoQJE3D48GGlY0VFRQgONn1uGUJsCQUvhBCixfXr1xEbG6t07NfLd1DPj4S0phDrHhuIIE9nxTlLTRsBQGBgIMrKypTW4DAMAwcH+quddG30J5wQQrToOPJys7gOy/+6Da5fBCYGCDGiu69Se0tOGwFATEwMbty4AQBobGyEq6urxZ5NiLVQ8EIIIVq0n4apbxVj4fYrEEllGBMfB29hiUp7S04bAcpbpjMzM2mnEbkrUPBCCCEayGQysFgssFgsMAyDxbuuo6C6GSFezti48D5kZmaqXNPQ0AA3N9NUkNbHyJEjcerUKQC0WJfcPSh4IYQQDXJyctC9e3cAwNazeUi8VQpHNgvfPDYIPh7OahPWAW2Lei3FyckJLBYLLS0tVE2a3DXMGrysXLkSQ4YMgbu7O/z9/TFr1iykp6frvO7EiRMYPHgwnJycEBkZiY0bN5qzm4QQopZ8vcvVwlqsOJAKAFg6rTdiQz0B/Ltg1tpGjx6NU6dOUTVpctcwa/By4sQJvPjiizh//jwOHz4MiUSCSZMmoampSeM1ubm5mDZtGkaNGoWUlBQsXboUr7zyChISEszZVUIIUXHt2jV0j+6LF7dfgVjKYGq/QMwfEa44HxMTg2vXrlmvg/+Qr3upq6sDn8+3dncIMTuzVg5LTExUer9lyxb4+/vj8uXLGD16tNprNm7ciG7dumHNmjUAgN69eyM5ORmrVq3Cgw8+aM7uEkKIklu3bqEkfBKKa1sQ6u2MTx+KUZoSiomJweXLlzFp0iQAbWtkrKFv3764efMmOBzzFoMkxFZYdM1LXV0dAGhdiX/u3DnFXwRykydPRnJyMsRisUp7oVCI+vp6pRchhJhCTmkNDqXXguPAwtpHBsLDyVHpfExMjFJ9oYaGBquMfLBYLAQGBkIgEFj82YRYg8WCF4ZhsGjRIowcORL9+vXT2K60tBQBAQFKxwICAiCRSFBZWanSfuXKleDz+YpXaGioyftOCLn7XMstQVaNBACwaFJPDOymmnjOx8cHVVVViveWTFDX0X333Yf+/ftb5dmEWJrFxhhfeuklXL9+HadPn9bZtuNKfXn2SHUr+JcsWYJFixYp3tfX11MAQwjpFKFEigVf7QHbJwz3RPlgwejuGts6OjpCJBKBy+VaPEFdew8//LDVpq0IsTSLBC8vv/wy9u7di5MnTyIkJERr28DAQJUqqeXl5eBwOPDx8VFpz+PxwOPxTNpfQsjd7bPEdGSk3oJXaBS++M8AODho3vrcu3dvpKeno3///hZPUNcRlQUgdwuz/klnGAYvvfQSdu/ejWPHjiEiIkLnNfHx8SqFxg4dOoS4uDg4OjpquIoQQkwjKb0cm07nQlSRh4+fmoYADyet7dvvOLLmtBEhdxOzBi8vvvgitm3bhh07dsDd3R2lpaUoLS1FS0uLos2SJUvwxBNPKN4vWLAA+fn5WLRoEVJTU7F582Zs2rQJb775pjm7SgghqGoU4q3f2xbgCtiNeHT8YJ3XtF+0a81pI0LuJmYNXjZs2IC6ujqMGTMGAoFA8dq5c6eiTUlJCQoKChTvIyIicODAASQlJWHAgAFYvnw51q5dS9ukCSFmxTAMluy+gYoGIaJ8XRDl76ZXptzo6GhF8k1rTxsRcrcw65qX9mXaNdm6davKsXvvvRdXrlwxQ48IIUS9XclFOHS7DI5sFt4a6Y2/8nRPcwMAh8OBRNK2K4mmjQixDFrdRQi56+VXNeHDv24BAN6YFI3WsjzExsbqfb2fnx/Ky8tp2ogQC6HghRByV5NIZXh951U0i6QYGuGN50ZFKmoa6SsmJgY3btygkRdCLISCF0LIXW1DUjauFNTCncfBF/+JBduBhRs3bmhNptlRbGwsrl27BolEQrsiCbEACl4IIXetG0V1+OpoJgDgo1l9EeLlgosXL6KlpQWurq5636djmQBCiHlRFS9CyF2pVSzFG7uuQiJjMK1/IGYNCMaBAwewYcMGbNu2zaB7+fn5oaKiAmw220y9JYS0R8ELIeSu9OWRDGSUNcLXjYuPZ/XHli1bcPz4cfz+++9GZe3mcDiQSqVm6CkhpCMKXgghd53L+dX47mQOAGDFA/2x/svPUFNTgx9//NHoFPvR0dHIyckxZTcJIRpQ8EIIuas0iyR447drYBjgwUEhyExKgIODA1avXt2p+8bExKC2ttY0nSSEaEXBCyHkrvLp32nIq2qGgO+ExRMjMffhN3Dw4MFO3zc+Ph719fUm6CEhRBcWo08aXDtSX18PPp+Puro6eHh4WLs7hBAbciarEnN/uAAA+Onpobh9LAEuLi5K9dUIIdZhyPc3jbwQQu4KjUIJFv9TdHHusG6Ij/DEsl27kJiYaOWeEUIMRXleCCF3hZUHUlFc24IQL2csndYbO3bswJw5cyipHCF2iEZeCCFd3tmsSmy/0Fa9/rMHY+DEYeGnn37CgQMHrNwzQogxaOSFENKlNQklWJzw73TRiChfJCQkYMaMGUblcyGEWB+NvBBCurRPE9NQVNOCYE9nLJnWGzKZDN9//z3+/PNPa3eNEGIkGnkhhHRZ57Kr8NO5fADApw/GwI3Hwb59+zB+/Hi4uLhYuXeEEGNR8EIIsUtvv/02rl27pvF8s0iCt/+ZLnp0aDeM7OELANi4cSMWLlxokT4SQsyDghdCiF26fv061qxZo/H8qoMZKKhuRhDfCUun9QIAyGQycDgcygFFiJ2j4IUQYpc4HA6amppw584dlXOX82uw5WwuAGDF7P5wd2rbDl1XVwdPT09LdpMQYgYUvBBC7I5MJgOLxcJLL72Eb775RumcUCLF2wnXwTDA7EHBGBPtrzhXUVEBPz8/S3eXEGJiFLwQQuxORUUF/P39MWrUKCQnJ6OpqUlxbt3RLGSVN8LXjYf37++jdF1lZSV8fX0t3V1CiIlR8EIIsTsFBQXo1q0bWCwW5s+fjx9//BEAcOtOHTacyAYALJ/ZF54uXKXraOSFkK6BghdCiFVduXLF4Jwr+fn5CAsLAwA89NBDSEhIgEjcVrtIKmMwtV8gpvYXqFxHIy+EdA0UvBBCrOrYsWNatzyrIx95AQBHR0dMmTIFr3y2Cbfu1IPv7IhlM/uqva6yspJGXgjpAih4IYRY1bVr11BXV2fQNe2DFwCY+MCj2LZ1MwDgvfv7wN/dSe11FRUVNPJCSBdA5QEIIVaVlZUFLperu2E7RUVFCAkJAQDIZAxWHC2EAz8QfRwr8eCgYI3X0bQRIV0DjbwQQqxGJBLBy8vL4JEXoVCoKKr4y6UCXMytht+wmfAsOAEWi6XxupqaGnh5eXWqz4QQ66PghRBiNampqYiNjYVQKDTq+tK6VnxyIA0AsOSRMWisqdTaXiaTwcGB/tojxN7R/8WEEKu5du0aYmNjDbqmubkZLi4uYBgG//3jJhqEEgwI9cRT90SYqZeEEFtDwQshxGqMCV7ki3UP3CjFkdQyOLJZ+PTBGLAdNE8XEUK6FgpeCCFWk5mZiR49emhdp9JRQUEB/AKD8MHemwCAF8ZEITrQHQDA4/HQ2tqq9rr262QIIfaNghdCiFUwDAOpVAoOhwOGYfS+rqCgAOfKHVDZKEKUvxteHNtdcS4wMBClpaVqr6OdRoR0HRS8EEKsoqSkBAJBWxZcJycnjSMmHZ1OScWFcgewWMCnD/YHj8NWnNMVvFCCOkK6BrMGLydPnsT06dMRFBQEFouFP/74Q2v7pKQksFgslVdaWpo5u0kIsYL26134fL5e26VbRFL8feEm2B5+eHxYGAaHeSudFwgEKCkpUXstJagjpOswa/DS1NSE2NhYfP311wZdl56ejpKSEsWrR48eZuohIcRarl+/bnDwsuZoBhrqahHk74vFU6JVzmsLXmjkhZCuw6wZdqdOnYqpU6cafJ2/vz88PT1N3yFCiM24fv06nn/+eQBtwUttba3W9jeL6/DDqVyAYfDxA/3h7uSo0iYwMBDnzp1Te31FRQWio1UDHkKI/bHJNS8DBw6EQCDA+PHjcfz4cWt3hxBiBnV1dYpfUjw9PbWOvEikMryz+zokEgkC+C6Y2CdAbTtdIy80bURI12BTtY0EAgG+++47DB48GEKhED///DPGjx+PpKQkjB49Wu01QqFQKTtnfX29pbpLCNFh7969AIAZM2YoHW9tbVXatqxr2mjzmVzcLK6Hi7QRQ2KjNLbz8/NDRUWF2nM0bURI12FTwUt0dLTSsG58fDwKCwuxatUqjcHLypUrsWzZMkt1kRBigF27dqGoqAjTp09XyuVy69Yt9OnTR/Gez+ejqqpK7T0KqprxxeEMAMBjfZzBLovU+Dw2mw2pVKr2HC3YJaTrsMlpo/aGDx+OzMxMjeeXLFmCuro6xauwsNCCvSOEaMIwDMrLyzFu3DgcOnRI6VzHzLqenp5q17wwDIOle26gVSxDfKQPergK0a1bN6P609raCmdnZ6OuJYTYFpsPXlJSUhS5INTh8Xjw8PBQehFCrC8tLQ29e/fGwoUL8c033yid6xi8aJo22n2lGKezKsHjOGDF7P4oLCzUGbywWCyNoy+EkK7BrNNGjY2NyMrKUrzPzc3F1atX4e3tjW7dumHJkiUoLi7GTz/9BABYs2YNwsPD0bdvX4hEImzbtg0JCQlISEgwZzcJIWZw7NgxjB8/Hj4+PoiIiEBycjLi4uIAANnZ2eje/d/MuOqCl8pGIZbvvw0AeG1CT0T4uqKgoAD/+c9/tD7Xz88PVVVV8Pf3N/FPRAixFWYdeUlOTsbAgQMxcOBAAMCiRYswcOBAvP/++wDaMmwWFBQo2otEIrz55puIiYnBqFGjcPr0aezfvx+zZ882ZzcJIWZw4sQJxVq1119/HV9++SWAtqkghmHg4PDvXz/qtkov33cbtc1i9BZ44NlRbRWji4qKEBQUpPW56nYcGVJ+gBBi+8w68jJmzBitf2ls3bpV6f3ixYuxePFic3aJEGIBUqkUjY2N4PP5AIDw8HCw2WxkZ2fD0dERoaGhSu35fL7STsHjaeX48+odOPxTAsCR3RboSCQSODqq5ndpLzAwECUlJUrTUu23ZRNC7J/Nr3khhNifa9euYcCAAUrHFi1ahC+++EIps64ch8NRrFNpEkrw3z/aKkY/MzICMSGeinb6VJ8WCAQq9Y0oxwshXQsFL4QQkzt69CjGjRundGzAgAEoLi7GkSNHVIIX4N+pnVWH0lFc24IQL2e8PrGn4nxDQwPc3Nx0PlvdtBFtkyaka6HghRBicmfOnME999yjcvzll1/G+vXr0b9/f7XXXSmowdazeQCAFQ/0hwv335ntgoICvbZJq6ssTQnqCOlaKHghhKjozAJXkUgEqVSqNqfKuHHj8MEHH8Dd3V31mQDeSbgOhgEeHBSC0T2Vg42CggKEhYXpfL58zUt7NPJCSNdCwQshRMXixYuNTvh48eJFDB06VO05FouFd999V+25/KpmpJc2wMeVi//e11v1fH6+XiMvLi4uaGlpUTpGIy+EdC0UvBBCVGRlZWnNbK2NPL+LQc8rb0BerQSMuBUfzOgLL1euSht9p43UoZEXQroWCl4IISoqKyuRk5Nj1LWXLl3CkCFD9G4vkzF4O+EGwHNBfKgzpseoz6it77SROjTyQkjXQsELIURFTU0NcnNzDb6uubkZjo6OOnOxtLf9Qj4u59eA5+yGF+IFGrdD19XV6V3+w83NDQ0NDYr3NTU18PLy0rtPhBDbRsELIUSFj4+PzuBl8+bNeOutt5TS+mvaZaTJndoWfPJ3GgBgbEw4eIzQuA530DHXi0wmU8roSwixb/R/MyFEiUwmA5/PVxq5UCcpKQlDhw7Fgw8+iI0bN0IikeDYsWMq+V00kVeMbhJJMTjMC2P6h6utLA20ZdblcPRPCK4u1wshpOug4IUQokTfVPrV1dV4+OGHcfDgQbDZbEyZMgXHjh1Tm4BOnT+uFiMpvQJctgM+fTAGnp6eaitLA2110LRVl+9I3XZpQkjXYdbaRoQQ+yNPpd/a2orm5ma4uLhobc9ms/Hcc89hzpw5uHjxol7TM5WNQiz7q61i9KsTeiDK3w1pfD6Ki4vVti8qKlKph6SNQCDArVu3ALTlnTFkDQ4hxPbRyAshRElVVRV8fHwQERGhcd2LvE17Hh4emDBhgl7P+HDvLdQ2i9FH4IH/Gx0JQH1labni4mKd1aTbaz/yQnWNCOl6KHghhCiRf9lrC14yMzPRo0cPo+5/6FYp9l0vAduBhc8eilFUjNY2bXTnzh0EBwfr/Yz2a15omzQhXQ8FL4QQJZWVlTpHXjIyMtCzZ0+157SpaxErKkb/3+hI9AvmK87x+XyNwYuhIy/e3t6orq4GQCMvhHRFFLwQQpRUVVXpHHkxNnj53/7bKG8QItLXFa+OVx650Ra8GDry0j5XTEVFBY28ENLFUPBCCFEiH6kICwtDfn6+2jaZmZmIiooy6L4nMirwW3IRWCzg04di4OTIVjrv7u6ucXu2IQnq5OTFJWnkhZCuh4IXQogS+WJcHo8HkUiktk1TUxPc3Nz0vmdDqxhLEq4DAJ6MD8eQcG+VNg4ODhqrWTMMozHzriaOjo4Qi8U08kJIF0TBCyFEiXzNC9AWNHQMKIwJJFb+nYY7da3o5u2CxVOiNbbTFLwYIyAgAGVlZTTyQkgXRMELIUSJWCwGl9tW1dnHx0ex8FXO0IRxZ7MqseNCAQDg0wdj4MI1LL1UfX29wVNGwL8lAqiiNCFdDwUvhBCN1C3azczM1HuxbpNQgsX/TBc9Prwb4rv76LhC1Z07dwzaaSQnz/XS0tKiM9EeIcS+UPBCCNFIXfCSkZGhd46XzxLTUFTTgmBPZ7wztbfO9g4ODpDJZErHiouLDdppJEf1jQjpuih4IYQodFzPoil40Wfk5Vx2FX4817Zb6ZMH+8ONp3u6SN2OI0O3Sct1rCxNCOk6KHghhCjU1dWBz/83cZy64CUnJweRkZFa79MolOCt368BAB4dGopRPfTb7ePp6alSIsDQBHVyVJyRkK6LghdCiELHnTlBQUG4c+eOUhuhUAgej6f1PisOpCqmi969r4/ez1eXqM7YkZeAgABkZGQoBWOEkK6BghdCiELHgotsNhtSqVTxXiqVgs1mq7tU4WRGhWJ30ecPx+g1XSSnKXgxZHeTHJfLxZ07d2inESFdEAUvhBAFdTlRWCyWYhFtQUEBunXrpvH6+lYx3lYkowvDiO6GBQ7qKkuLRCKdIz2acDgcSlBHSBdEwQshREFe16i94OBgxdSRrsW6y/+6jZK6VoT7uODtqb0Mfr62ytLGEAgENPJCSBdEwQshRhKJRCoJ3Oxd++y6cu0X7WZmZmrcJn00tQy7LrfVLlr1cKzByegA1WkjmUxmcDbf9gIDA2nkhZAuiIIXQox04MABfP3119buhkmpmzZqH7xoGnmpahTi7YQbAIBnR0YgTk3tIn10DF46W5eIRl4I6ZoM/9WIEAIAyMrK6nJbcTsu2AXagpcDBw4AAPLz8xEWFqZ0nmEYLN1zA5WNQvTwd8MbkzTXLtKl41ZpY3cayb399tu024iQLoiCF0KMlJWVhbKyMmt3w6R0TRvJZDKV3UYJV4px8FYZHNksfDlnAJwcte9G0qbjyIuxOV7kvL2NGwEihNg2Cl4IMVJxcbFKKnt7JxaLVXb2+Pj4oLKyEiKRCI6OjkrnCqub8eHeWwCA1yf2RL/gzo1yqAteOjPyQgjpmmjNCyFGYhgGDg5d/38h+YLZjpl1pTIGb/x2DY1CCeLCvPD86O6dfpaLiwuampoU740tykgI6drM+jfvyZMnMX36dAQFBYHFYuGPP/7Qec2JEycwePBgODk5ITIyEhs3bjRnFwkxilAoBJfLBdAWxHR1XC4XN2/eVFqs+8OpHFzMq4Yrl40v/jMAbAfjdwXJddxZRCMvhBB1zBq8NDU1ITY2Vu8dGbm5uZg2bRpGjRqFlJQULF26FK+88goSEhLM2U1CDJaXl4eIiAh4eXmpJFWzV9qCsLCwMBw5ckQRvNy6U4fVhzIAAO9P74NuPi4m60f7AKazu40IIV2TWde8TJ06FVOnTtW7/caNG9GtWzesWbMGANC7d28kJydj1apVePDBB83US0IMl5WVhaioKHA4HJSUlMDLy8vaXeq0+vp6eHh4qD0XERGBL7/8Eu+99x5aRFK88ksKRFIZJvQOwH/iQk3aj/ZB1N0yNUcIMYxN/a1w7tw5TJo0SenY5MmTkZycDLFYrPYaoVCI+vp6pRch5iYPXgIDA1FaWmrt7piEuuy6chERESgvL0dQUBCW77+N7Iom+Lvz8NlDMZ1KIkcIIcawqeCltLQUAQEBSscCAgIgkUhQWVmp9pqVK1eCz+crXqGhpv0tkBB15MGLQCDoMrle1G2TlouIiECPHj1w8FYZdlwoAIsFfDlnALxduSbvh6OjI8RisdK6IkIIac+mghdAdcGefAhZ0293S5YsQV1dneJVWFho9j4SUlBQgNDQ0C4VvGgbeYmMjMTEqdPxzu62oov/NzoS90SZJ3Oth4cH6uvrUVJSQjuNCCFq2VSeF3VD8OXl5eBwOBp/I+TxeEZXnCXEWFKpFBwOp0tNG6krDSDHc3JGUcR9qM2pQv9gPt6YaHwWXV3klaVLS0speCGEqGVTIy/x8fE4fPiw0rFDhw4hLi5OJTkWIdYikUgUWWZNNfJiykrKxtI2bfTtyWycy6mCC5eNrx4ZAC7HfH91yCtLd7Y0ACGk6zJr8NLY2IirV6/i6tWrANq2Ql+9ehUFBQUA2qZ8nnjiCUX7BQsWID8/H4sWLUJqaio2b96MTZs24c033zRnNwkxSEFBgaK+j7u7e6cXiZeWluLJJ580Rdc6RdO0UXJetWJb9IfT+yLSz82s/ZBn2e1saQBCSNdl1uAlOTkZAwcOxMCBAwEAixYtwsCBA/H+++8DAEpKShSBDPBvAbikpCQMGDAAy5cvx9q1a2mbNLEp8sW6gOa1WIZIS0uz2ZGXmiYRXv4lBVIZg5kDgvBwXIjZ+yEPXmjkhRCiiVnXvIwZM0Zr4qutW7eqHLv33ntx5coVM/aKkM5pH7yYQmpqqk0ELx0rSjMMgzd3XUNJXSsifV3xvwf6W2RbtLyyNI28EEI0sak1L4TYg47BC4/HQ2trq9H3S0tLg0QiMUXXOkUoFMLJyUnx/odTuTiaVg4uxwFfPzYIbjzLrO+Xj7w0NDRoTJpHCLm7UfBCiIFyc3MRHh6ueB8YGIiysrJO3S8kxPzTMYa4UlCDTxPTAADv398HfYIsF0R0rCxNCCEdUfBCiIHEYrFS8rTO7jiSb7u2FbXNIry8IwUSGYP7YgSYO6ybRZ8v3yp9NxS8JIQYh4IXQgwglUpV1n10JtdLY2MjXF1d4eTk1Kmpp86SBwoyGYPXdl5FcW0Lwnxc8Mlsy6xzac/T0xMFBQXg8/kWfS4hxH5Q8EKIAYqLi1VKUHRm5CUjIwPR0dGKrLLW0tjYCHd3d6w9lomk9ArwOA7YMHcw3J0sn1+Jz+cjNTWVdhoRQjSi4IUQA6jbaSQQCIweeUlNTUXv3r2tvs6jsrISrWwXfHU0EwCw4oH+Fl3n0h6Px0NhYSHtNCKEaETBCyEGUBe8BAYGGj3ykpaWhl69ell95OVWdhHOFArBMMDjw7vhwcHWXUDs4uJCIy+EEI0oeCHEANnZ2ejevbvSMX9/f1RUVBh1v/T0dERHR1t15KVVLMVHCRcgdnTDgFBPvHd/H6v0oz0+n08jL4QQjSh4IcQA2dnZiIyMVDrGZrMhlUqNul9TUxNcXV2tNvLCMAz++8dN5BSVwsPLG+vnDgKPw7Z4Pzri8/k08kII0YiCF0IM0NraCmdnZ5XjxmzrbV/g0VojL1vO5OH3y0VgWurx2n2DEOSp+rNZg0AggEAgsHY3CCE2ioIXQvSkLUBhsViQyWQG3S8vLw8REREAAA8PD4sHLyczKvDx/tsAgBEhPIyJ7a7jCsvZtWuXUi4dQghpj4IXQvRUWlqqcTTA19cXVVVVBt1PvtMIaBt5seS0UW5lE17acQUyBnh4cAiCncQqRRmtiQIXQog2FLwQoidtBRmNSVQn32kEWHbkpb5VjGd/vIT6VgkGdfPExw/0Q3V1tU0FL4QQog0FL4ToSVvwYkyiuvbBi6VGXqQyBq/9ehXZFU0Q8J2wcd5g8DhsjWt5CCHEFlHwQoiedAUvho68lJaWIiAgAIDlRl5WHEjFsbRy8DgO+G5eHPzdnXRfRAghNoaCF0L0lJWVpZLjRc7QRHXyxb/yukEuLi5obm7ufCe1+OlcHjadzgUArHo4Fv1DqHYQIcQ+UfBCiJ6amprg5uam9pyhIy8VFRXw9/dXvDd38cOjqWX4cO8tAMBbk6MxPfbfBHAMw1i8+CIhhHQGBS+E6MnBQfP/LoaOvLTfaWRuN4vr8PIvKZAxwJy4UCwcozx6JE+URwgh9oKCF0JMwNXV1aBpn/aLdc3pTm0Lnt56Cc0iKUb18MXHD/RTGWWprKyknUaEELtCwQshVpCWlmb2kZe6ZjGe3noJ5Q1CRAe445u5g+DIVv1fPj8/H926dTNrXwghxJQoeCFED8ak/9cmJydHkV1XzpgsvZq0iKR45sdLSCttgJ87D5ufGgIPJ0e1bdUVmySEEFtGwQshetAnD4qzs7PeU0dSqRQcDkfpmJubG5qamozuo5xYKsNLO64gOb8G7k4c/PT0UARrqVmUk5NDwQshxK5Q8ELuOsePHzc4J0tjY6PGnUZy+iaqa25uVhsImaI4o0zG4J2EGzj6Ty6XzfOHoLfAQ+s1OTk5KpWyCSHEllHwQu4qFRUVmDNnDq5fv27QdfoEL/qWCMjIyEB0dLTKcQ8Pj05n2f0kMQ0JV4rAdmDhm8cGYUi4t85rGhoa4O7u3qnnEkKIJVHwQu4qixYtwmOPPWZwkGDKkZfU1FS1O406O/LyzfEsfHcyBwDw6YMxmNAnwOh7EUKILaPghdw19u7di5CQEIwZM8bgIEHf4EWfkRdNO406UyLgh1M5+PxgOgBg6bReeGhwiF7X1dbWgs+nTLuEEPtCwQu5K9TV1WHNmjV4//33jSqCqO+0kT4jL1evXkWfPn1UjhtbnPHHs3n4eH8qAOD1CT3xf6P1X3xLO40IIfaIghdyV3j77bexbNkyODs7GzXCYaqRl6amJrBYLLULdo3p144LBfjgn7T/L42Nwivj1ReO1IR2GhFC7BEFL6TLO3bsGBwdHTFq1CgAxo1w6BO8+Pj4oLKyUmubEydOYOzYsWrPGdqvXcmFWLrnBgDg/0ZH4o1JPQ2uUUQjL4QQe0TBC+nSpFIpPv74Y6xYsUJxzFwjLw4ODjqT2SUmJmLKlClqzxnSr18vFmBxQtuOqfkjwrFkai+jiitmZ2fTNmlCiN3h6G5CiP0qLS1Fz549lbYCG7OrR5/gRR/p6eno2bOn2nP6jrxsPZOLD/+6DQCYNzwMH0zvY3RV6NLSUgQGBhp1LSGEWAsFL6RLKywsRGhoqNIxHo8HkUhk0H30DV54PB7q6urU7uCRj3JoCjT0GXnZeCIbn/ydBgB4blQElk7rrbifVCrFmDFj4OXlpWjf0NCA6dOnY9GiRRrvaWzgQwgh1kLTRqRLKygoUAlejKFv8PLII49g69atas8dPHhQ45QRoD1JHcMw+PJwhiJweWVclFLgAgCHDx/Gfffdh7179ype+/fvx+nTp9XeUygUgsvl6vyZCCHE1lgkeFm/fj0iIiLg5OSEwYMH49SpUxrbJiUlgcViqbzS0tIs0VWtKioqDE4rT6yrsLDQJBWT9Q1eZs6cib1796otsHj06FGMGzdO47VcLhdisVjluEzG4H/7U/HV0UwAwFuTo7FoUrTKiMnWrVsxf/58pWMuLi5oaWlRuxYnPz8f4eHhOn8mQgixNWYPXnbu3InXXnsN7777LlJSUjBq1ChMnToVBQUFWq9LT09HSUmJ4tWjRw9zd1Wnixcv4o8//rB2N4gB1E0bGaOhoUGv4IXD4WDixIlITExUOi4UCiESiQxOwy+SyPDazqv44XQuAOC9+/vgxbGq26HLysogk8nUrl/p3r07cnNzVY7TYl1CiL0ye/DyxRdf4JlnnsGzzz6L3r17Y82aNQgNDcWGDRu0Xufv74/AwEDFi81mm7urOgUFBaG4uNja3SAGKCwsREiIftlmtWlqaoKLi4tebZ999ll8//33SsfOnDmDkSNHGvTMhlYxntp6EXuv3QHHgYUv58TimZERatv+9NNPePLJJ9WeGzZsGM6fP69ynLZJE0LslVmDF5FIhMuXL2PSpElKxydNmoSzZ89qvXbgwIEQCAQYP348jh8/rrGdUChEfX290stcgoODcefOHbPdn5ieSCQCj8dTOc5msyGRSPS+j0wm0zuA9vX1hY+PD9LT0xXHEhMTMXXqVL2fV97QijnfnseZrCq4cNnYPH8IHhioPghjGEbreprhw4dT8EII6VLMGrxUVlZCKpUiIEC5QFxAQIDGtSMCgQDfffcdEhISsHv3bkRHR2P8+PE4efKk2vYrV64En89XvEwxRaCJr68vKioqzHZ/YjmmqOCszUsvvYRvvvlG8f769evo37+/XtdmljXgwQ1ncbukHj6uXPz6f8MxuqefxvanTp3CiBEjNAZXUVFRyMrKUjmen5+PsLAwvfpECCG2xCJbpTsuLGQYRuP2zOjoaERHRyvex8fHo7CwEKtWrcLo0aNV2i9ZskRpG2h9fb3ZAhh9kpAR2yEUCtWOugD/Bi/e3t5mefaAAQOQnZ2tGA0MCgrSa0tyTasUD6w7iUYJEObjgh+fGopwX1et12zevBkffvihxvMsFgtcLhetra1wcnJSHBeLxbTbiBBil8w68uLr6ws2m60yylJeXq4yGqPN8OHDkZmZqfYcj8eDh4eH0osQACgqKtK43sWYRHWGmj9/Pn788UedW6SBtoB+8+lcXC0Voa6hHkMjvLFn4T06A5fa2lrU1tbq3DU0aNAgpKSkKD2P8rsQQuyVWYMXLpeLwYMH4/Dhw0rHDx8+jBEjRuh9n5SUFAgEAlN3zyg8Hg9CodDa3SB60JbjxdgKzoaYNWsW/vjjDxw6dAgTJkzQ2E4slWHpnpv4aN9tsLjOmNrTA9ueGQZvV92jItu3b8fcuXN1tuu47qWkpIQy6xJC7JbZp40WLVqEefPmIS4uDvHx8fjuu+9QUFCABQsWAGib9ikuLsZPP/0EAFizZg3Cw8PRt29fiEQibNu2DQkJCUhISDB3V/USFBSEO3fuICJC/a4PYju05Xgxpr6RoRwdHTF+/HgcOHBA4/RUeX0rXtqRgot51WCxgHv7heGFEUHgcnT/XsEwDP7880/s27dPZ9uhQ4diy5Ytive0WJcQYs/MHrzMmTMHVVVV+Oijj1BSUoJ+/frhwIEDioWCJSUlSjlfRCIR3nzzTRQXF8PZ2Rl9+/bF/v37MW3aNHN3VS8UvNiPgoICjSMehoy8iEQio9eGLFiwAEOGDFF77kJOFV7ckYLKRiHceBysfXQAruxL1btfV65cQUxMjF598/T0RG1treI9BS+EEHtmkQW7CxcuxMKFC9We65hKffHixVi8eLEFemUcyvViP7QlqPPw8EB+fr5e9+lMUUZvb29MnDhR6RjDMPj+VA4+TUyHVMYgOsAdGx4fhEg/N2Sf0j+ounjxotaMvR0FBgYqCjFmZ2dj9uzZBv0shBBiK6i2kYEo14v90FYx2ZAFu6aqKA0AdS1ivLDtClYcSINUxuCBgcHY8+IIRPq13d+Q6SxDK0IPHz4cFy5cAADk5ORQdl1CiN2i4MVA8mkjYvu0JZYzJM+LqYKXi7nVmPbVKSTeKoUjm4Xls/rhi//EwoX77wCoIdNZZWVlBu3aa59pt76+Xm3la0IIsQcUvBiIpo26BkuOvIilMqw+lI5HvjuH4toWhPm44PcFIzBveJjKdmVDRl7Ky8vh7++vdz/69euHGzduGNR3QgixRRZZ89KVmDszKzGNuro6rTl/DBnhaGxsNLigolx+VRNe/fUqrhbWAgAeGhyCD2f0hRtP/f96hvRLIpHA0dFR775wOG3PrKmpMfrnIYQQW0AjLwaytcRe7bMLk3/pqibt5uaGhoYGve5lzMiLTMZgy5lcTFlzClcLa+HuxMG6Rwdi1cOxGgMXwPxbuPv164e9e/fSTiNCiF2j4MUIDMPYRJmAhoYGrF+/3qACg3cLbTlegLbCjDKZTK97GRq85FY2Yc5357Dsr9toEUsxLMIbf786CtNjg3Rea+7kecOGDcP27dspeCGE2DUKXozA5/P1/q3dnLKysiAWi5GdnW3trtgcbdl1DaVv8CKRyvDDqRxMWXMSl/Jq4MplY/msfvjlueEI8XLR61n6jgg1NTXBxUW/e7Y3bNgwHD16lHYaEULsGgUvRujsol2JRIK8vLxO9yMzMxOjR4+mRZhq6Bp5AaD36Jk+wcuVghrM+PoMPt6fCqFEhlE9fHHw9dGYNzwMDg76TzXqW/zT0J1GckFBQQgJCaGRF0KIXaPgxQidzfVy6dIlvPfee53uR1ZWFmbNmoWbN292+l5djaVGXqqbRHj79+uYvf4sbpfUg+/siE9m98dPTw/Ve7TFGMYGLwCwcuVKm6kVRgghxqDgxQidzfVy+fJlk2y3zszMxMyZM3Hr1q1O36urqampgZeXl852+oxyqAteJFIZtl/Ix7jVSdiZXAgAeHhwCI69cS8eGdrN7Au7OxO8PPbYY3BwoP/1CSH2i7ZKGyEoKAinT582+vrk5GS0trZ2uh8VFRUICwtDc3Nzp+/VFekKIJydndHa2gpnZ2et7doHLwzDICmjAisPpCKjrBEA0CvQHR/P6oe4cPXFF82hrKwMwcHBFnseIYTYEgpejNDZaaPy8nJ4e3uDYZhO/4bOYrHA4/HQ2toKJyenTt2rq5DJZHp9rvJEdfoGL7fv1GPFgVSczqoEAHi6OOLV8T0wb3gYOGzTjWTId7Np+xnKysowaNAgkz2TEELsCQUvRhAIBEYHL83NzXBxcQGPx+tUivb6+npForHevXsjLS0NAwYMMOpeXU1FRQX8/Px0tpNvS9ZVH6isqhbv7cvAvlsVYBiAy3bA/HvC8eLYKPCd9U8Spy8XFxe0tLRo3U3UmWkjQgixdxS8GMHJyQlCodCoa69evYrY2Fg0NDSguLjY6OAlKysLPXr0ANCWeOzmzZsUvPyjoKBA504jQHdCuIKqZqw9lonTGWVIv1kBALg/RoC3p/RCqLf5FuPKR4S0BS+GlgYghJCuhIIXC7t8+TLi4uKQkZGB4uJi9OnTx6j7ZGVlISoqCkBb8LJ9+3ZTdtOu6cquK6cpIVx2RSO+PZGNhCvFkMoYMAAm9PbHaxN6ol+w+YsZyvulbUeQUCikaUJCyF2LghcjsVgsyGQyg3dtXL58GXPmzEFTU1OndhxlZmZi3LhxAIDo6Gikp6cbfa+uprCwEH379tXZruPIy9XCWmxMysbB26WQb0K6t6cf8sK98cOTQ8zVXZ39IoQQooz2SxrJ19cXlZWVBl8nH+4PDg7udPAinzbicrkQiURG36ur0TfHC5/PR01NLY7cLsMj353DrG/OIPFWW+Ayobc/El6Ix49PDzXLuhZd/aLin4QQohmNvBhJvuPIkHUH7VO6dzZ4qayshI+Pj+K9vNq1tkrKdwt9po1qm0U4ldeIP5MuQ5bZtvCV48DCzAHBeP7eSPQMsF7VZV0jL62treDxeBbsESGE2BYKXowkLxHQcZFseXk58vLyMHToUJVrrl27pmgvEAhQUlJi9PNZLJbSVtq+ffvi9u3bGD58uNH37CrkO7o6YhgGKYW12HmxEH9eK0ZdQRVaauoQ5uyIR4aE4skR4Qjy1L5t2hL4fD5qa2s1nqedRoSQux0FL0bSlGV3x44d2LNnD06cOKFy7vLlyxg8eDAAwNHR0ehq0O23ScvJdxwZG7wUFRUhJCTEqGttTcf8KFWNQuxJKcbOS4XILG9UHI8K9gdf4oJdS8bDmctWey+pVGrxbLQeHh7Iz8/XeJ6CF0LI3Y7WvBhJU6K648ePIzQ0FJmZmSrnkpOTFcFLZ7TfaSQnD16MUVdXh4EDB6KmpqbTfVNH2xexqYnFYnA4HLSKpdh/vQT/91Myhq88io/3pyKzvBFOjg6YPTAYvz0fj99eHgeBC6MxcAHapvr0qShtSrrWvFDwQgi529HIi5HUVZZuaGgAm83Giy++iC1btmDFihVK5ysrK5XWyHA4HEgkEnA4hv1naL9YVy4iIgLZ2dkG/hRt/ve//yEiIgKlpaV61QMyhFgsxujRo5GXl2f2ej8iiQy7T15DbgsPg5cfRpNIqjjXP5iPOUNCMWNAEDyc2hbgtrS06FwYq09FaVPTteaFghdCyN2OghcjBQQEoKysTOnYkSNHMHHiRAwfPhzvvfeeUmDS1NSkkoY+MDAQpaWlBk/XZGZmYsKECUrH2Gw2ZDKZwT9HTk4OcnNz8cgjj6CkpAS9e/c2+B7apKeno7CwEHl5eYiIiDDpvQGgRSTFiYwKHLxViiOpZajIvAphqzP4IimCPZ0xY0AQZg4IQq9A1YXMTk5OaGlp0Xp/awQv+oy89OvXz4I9IoQQ20LBi5HUBQv79u3DBx98ABaLhcmTJ+PQoUOYNm0agLbMugMHDlRqL99xZGjwkpWVhQULFqgc9/PzMzjz6rvvvouPP/4YV65c6dQCYk2uX7+OSZMm4dy5cyYLXoprW5CUXo7jaRU4nVWBVvG//x2cRbWIH9wb/30hHoO6eWkd7WGxWDqrStPICyGE2B4KXkxEJpOhqKhIkZb+iSeewKuvvqoIXtov1pUzdrt0VVWV0jZpuX79+uHWrVt6By+nTp2Cv78/oqOjUVJSguTkZIP7osuNGzfw/PPP49ixY3jssceMukeLSIrk/GqczqxEUnoF0ssalM6HeDljSt9ATOkXiEO/XMfoUSMwOMw0FZ6tEbw4OTlprTpOwQsh5G5HwUsncDgciMViODo6IiUlRanKb0BAAGQymaJI4OXLl/Hoo48qXR8cHGxUZlxNFYfli3bHjh2r8x4ymQwff/wxfvnlFwCd37qtSVpaGpYtW4bvv/9e72uEEiluFtfhXHYVzmRV4XJ+DUTSf0dXHFjAoG5eGNvLH2Oi/dBH4AEWi4Xq6mocOXwYzzz9tMn6b43gRdfaIE1bwQkh5G5BwUsnCAQClJaWIjQ0FPv378d9992ndP7xxx/H9u3b8dprr6mtdBwUFIRjx44Z9My6ujqNiej69euHPXv26HWfbdu2Ydq0afD29lb6WUxNJBKBy+XCxcVF45duXbMYKYU1SM6rwcW8alwrrIVQojwlJ+A7YUR3X9wb7YfRPXzh6cJVOp+dnY3nnnsOq1evNumohDWCFwA6p7MIIeRuRsFLJ8h3HIWGhuL8+fNYunSp0vmpU6di2rRpeO655+Dq6qpyvTHTRu2rSau7X1FRkc57NDU14ccff0RiYqLimLu7u8lT0ldXVyt2L8XFxSE5ORlD4+9BRlkDrhXWIqWwFlcLa5FT0aRyrY8rF0PCvXFPD1/c090HEb6uGkckzpw5gw8//BA//fSTweuH2Gw2pFIp2Gz126UbGxsVAR4hhBDbQMFLJ8hzvZSVlcHb21tly7OjoyMGDhyI77//XiUTLwB4enpqzaSqjrbgRZ51V9O0ktzx48cxe/ZsODo6Kl1ramcvXgE/uDu2nMlFitAPmz/fAVn/RkhkqqMKYT4uiAvzxpBwLwyJ8EaklmClvV9//RW7du3C7t27VRL36cPd3R2NjY3g89VXi25sbFSsY7IFQqEQXC5Xd0NCCOnCKHjphKCgIGRkZODvv//G1KlT1bZ56qmnEB8fj507d6qcMyZgULdNur3Q0FAUFRVpre1z8+ZNjBw50uBna1LdJEJWeSOyKxqRVd6IjLIGpJU2IDvpdzh6BWH/X7chE3miKuMG/PpOh7crF/2C+RgQ6omBoZ6IDfWEt6vhX8ipqalISEjAzp07Dc6VIyff2aMteLHGtJGDg4PaESFDd5MRQkhXRMFLJwQFBSEpKQl5eXn4+uuv1bbp3bs34uLiNGbWZRhG50hJe5mZmXjhhRc0nu/Xrx9u3LihM3h5/vnnVY7L8550zEfDMAxqm8UorGlGXlUz8iqbkFfVhLzKJuRWNqGmWaz2OeLyXEQNG4eYHgHoG+SBX5KdsOftsQjydDbJSE9SUhLmzZtndOACtOVU0bYt2VrBi4eHBxoaGuDp6al0nHYaEUIIBS+dEhwcjPz8fNTX16vduix36NAhjfVx5AnJNP3m31F1dbXWZ8XExOD06dOKLdrq1NbWKmXSlUhlqGgUgu3mje3Hr0Lm5o87ta0oqmlBUU0zimpa0CjUXocp2NMZ3f3dEOXnhh4BbogOdMeSa98gcfnDikAl+8/eaK0uAcuru14/qy6nT5/G2rVrO3UPeTVuTawVvMj/XFDwQgghqih46QQvLy8cPXoUr732mtZ22gr7yRft6hu86BIbG4tvvvkGANAskqCqUYTKRiEqG0Uob2hFaW0TMsqb8OyPyW3v61pR2SiEjAHqclqRtOMMnELUZ9n1d+chzMcF4T6uCPd1Vfx7pJ8rXLjKf5SkUim4HLbSCEt8fDzOnTuH7t07H7wwDKMzkNOHLY+8qOsXBS+EEGKh4GX9+vX4/PPPUVJSgr59+2LNmjUYNWqUxvYnTpzAokWLcOvWLQQFBWHx4sVqM8paG4vFgqurq8oWaUPIg5c+ffpobCOUSFHfIsGd8kpI2U44nl6OumYx6lrEqG0Wo6ZZ9M9LjJomEU5fzUOv9/5WyjwrJ64qQiO8cCRVubQB24EFX/8A+LoIMSQ2CEGezgj2dEKItwtCvVwQ4uUMJ0fNBQw7ysnJUQlS4uPjsXr1ajz++ON630cTU5UbsNWRl6CgIBQUFKB///5Kx8vKyhAdHW3x/hBCiC0xe/Cyc+dOvPbaa1i/fj3uuecefPvtt5g6dSpu376tdhdHbm6uYnvxtm3bcObMGSxcuBB+fn548MEHzd1dg73zzjtKXzBSGYNWsRQtYilaRFLFvzcJpWgRS9r+KZKiSSRBk1CCS+XAsewL+LvGD01CCRqFEjS2StAglKChVYL6FrEi54mwNAstlY54asslrX0SsbhobmyEA88FPI4DfFy58PNwgp8bD1U308EbPQT3zeqHQA8nBHo4IYDPg48rD4cPOSArKwsvPTpQ6/31cf36dcTExCgdi4iIQG5ubqfvDbRlB9YWAOuLz+dr3a7e2toKHo/X6ecYavz48fj2229VAmMaeSGEEAsEL1988QWeeeYZPPvsswCANWvW4ODBg9iwYQNWrlyp0n7jxo3o1q0b1qxZA6BtwWtycjJWrVpl1eClVSzFwu1XIJRI0SqWQSiRQiiWoVUShe9XHEWruO19+0yw+hAWi9FakIlrLspfoAzDQFh0CwDAcnSCg6MTWBVZ8A8JQ48gD3i6OILv3PbycuHCy4ULTxdHeLtysUd2L8aM9sTUiWPhwlWeulm2bD+mTZuGIUPCVPoiEAhw6tQpvfsuk8lw5MgRTJo0SeXcjRs3cP/99ysdk49UNTU1qc17Y4iTJ0/iww8/7NQ9gLbgJTU1VeN5+fZzS+vVq5faflHwQgghZg5eRCIRLl++jHfeeUfp+KRJk3D27Fm115w7d07ly3Dy5MnYtGmTIhW/NTiwWDiWVm7QNU6ODnDhcuDsyIYzlw0XxYsDFy4bbjwORD15uPBHCp6Y2gtuThy48Thwd+KgsigXa/93GGPGjoNEWAlRazNaHaVYuvRFnYnYRGNHIDMzFa68iSrnbt26hTfffFPtdYaWCEhLS8NDDz2EoqIilay/N2/exOLFi1WuGTJkCC5duoQxY8bo/Rx1jCloqY6uIojWwmKxEBYWhry8PISHhyuOW2saixBCbIlZg5fKykpIpVKV3xQDAgI0pqIvLS1V214ikaCyshICgUDpnFAohFAoVLw3dZZYOUc2C589GAOeowN4HPY//2z7dydHBzg7suGkeDnAicOGg4Pu39jFYjHm7FuP5+9VXh+yat/PeO+dtzB58mSD+zpw4ED89ttvas81NzdrHPXw9fVFVVWV3s85f/48pk6dih9//BEvv/yy0rmWlha1pQDi4+Nx+vTpTgUvphx90LVg15pp+qdMmYKDBw+qbGu3xkgQIYTYEs3bYEyo41+2uvKaqGuv7jgArFy5Enw+X/HSlt+kM1gsFv4zJBQzBwRjSr9AjI32x4juvhgc5oW+QXxE+rkhyNMZ3q5cuHA5egUuQFsWXolEdRtyUlKS0V/wgYGBKCsrUzkuFAq1rt9wcHAw6Mv6/Pnz+OSTT/DHH39AJvt3uqyxsVFjgCQvE9AZp06dMlmSPV0Ldq1p/PjxOHLkiLW7QQghNseswYuvry/YbLbKKEt5ebnG35wDAwPVtudwOGq3xS5ZsgR1dXWKV2Fhoel+ACuprq6Gq6trpxaKcrlcpREpAEhPT9e5U8WQ4KWoqAjh4eGYMGECDh48qDh+8+ZNlV0ycs7OzmhpaenUiIapFusCth288Pl8NDc3QyxuSwIoFos7lZCPEEK6CrMGL1wuF4MHD8bhw4eVjh8+fBgjRoxQe018fLxK+0OHDiEuLk7tehcejwcPDw+ll71hs9lKoy+JiYmYMmVKp+7Zt29f3L59W+nYzZs30a9fP4P6oklDQwPc3NzAYrHw3HPP4fvvv1ecU7fTqL2oqChkZWXpfIYmmZmZ6Nmzp9HXt8fhcDT+vDKZzOpTNPLcOADUViYnhJC7kdmnjRYtWoQffvgBmzdvRmpqKl5//XUUFBQo8rYsWbIETzzxhKL9ggULkJ+fj0WLFiE1NRWbN2/Gpk2bNC4y7Qo6jjbt379fa4ZcfQwcOBApKSlKx/QJXvz9/VFernthcnJyMuLi4gC0jbD5+PggPT0dgO7gZeTIkUhKStL5DHXq6urg7u5ukaBC07odS5KvewFopxEhhMiZPXiZM2cO1qxZg48++ggDBgzAyZMnceDAAYSFtW3VLSkpQUFBgaJ9REQEDhw4gKSkJAwYMADLly/H2rVrbTLHi6nIE9UBbVMD1dXVnf6SUhe86DNtpO+OowsXLmD48OGK9y+99JKivlNeXp7iv68606ZNw759+3Q+Q52zZ89qHLUzNVvY2TNo0CBcvnwZAAUvhBAiZ5EJ9IULF2LhwoVqz23dulXl2L333osrV66YuVe2IygoSBG8nDlzxiSLUdUlhNO1YBdoC1407QRrLzk5WWmHUWxsLPLy8lBXVweGYbSWRHBzc4ObmxtKSkpUdo/pcurUKTz00EMGXaOLplEcWwheHBwc4Ovri/LycgpeCCHkHxbZbUS0az/ysm/fPpXkbsaQBw/yXUBNTU16TYHoM/LCMIzaRHNPPfUUPvroI7WZkzuaM2eOxu3c2ly/fh2xsbEGX6eNpsXDthC8AG15kQ4dOkTBCyGE/IOCFxvQPnjRtV7EEO0Xxt6+fRt9+/bVeY0+wUt+fr5S4jS5GTNm4I8//tCr/1OmTEFiYqLOdu21trbC0dERbLb+NZb0wePxVHZmARS8EEKIraLgxQbIg5eMjAz06NHDZItR26970WexLtC2eFhX8NJxvYsch8PBBx98gHvvvVfnc7hcLkJDQ5Gdna2zrdzFixcxZMgQvdvrS1OiOlsJXgIDA1FeXq42gSMhhNyNKHixAZ6enqitrTXZlJFc++Dl1q1beo28qMuz09H58+cxbNgwteeeeOIJrRWy23v00Ufxyy+/6NUWaFvvMnr0aL3b60tTrhdbCV6Atv+Wp0+fBp/Pt3ZXCCHE6ih4sQHykZbjx49j7NixJrtv7969FcX9MjMzERUVpfMaLperSIqmSXp6uknyrIwePRonTpzQuOakvr4eR44cwccff4zZs2fj3Llziu3ZpmTrIy9AW30vXZmpCSHkbkHBi42QSCRwcnKCk5OTye7p6OgIsVgMhmEglUr1zs6q7QtSKBSCy+Vq3U2kLzabjf79++PatWsq57788kvMnTsX165dw7hx47Bjxw7s27fPpJ+PnD2MvIwYMUKv6ThCCLkbUK5xG+Ht7Y1x48aZ/L7BwcG4ffs2PD09DbpO02/5165dw4ABA0zTOQCPPfYYfvnlF6V7bt++HVlZWdi7d69FRhrsYeSFy+Xi559/tnY3CCHEJtDIi414/PHHMWvWLJPfd+DAgdi2bZte613kvLy8UFNTo/actvUuxhg8eDAuX76s2NJ96NAh7Nu3D2vXrrXYFIk9BC+EEEL+RcGLjZg6darawpOdNWDAAGzfvl2vnUZy2nYcXbx4EUOHDjVV98BisXDPPffgzJkzuHz5MtauXYstW7aYfDu0NvYwbUQIIeRfNG3UxcXExKCoqMig4EWe60XdaE11dbXJg6xHH30Ub731FlpaWvD777+bZV2LNjTyQggh9oVGXro4Nzc33HPPPVprDXWkqURARUUFfH19Tdk9AECvXr0QHByMH3/80eC1OaZAIy+EEGJfKHi5Cxw/ftyg3UGasuxqSk5nChs3bkRwcLBZ7q2LppGX5uZmODs7W6FHhBBCtKHg5S6g7xZpOU3By5kzZ0y6WNdWaBp50VVgkhBCiHXQ38xEhbrgRSgU4uLFixg4cKCVemU+rq6uaGpqsnY3CCGE6ImCF6LCzc0NjY2NSsd27tyJOXPmdMmRCBaLBYlEYu1uEEII0VPX+yYiJscwDLZt24Z58+ZZuytm4+Pjg7y8PGt3gxBCiB4oeCFqtU8Ql5SUhOHDh3fpxasLFy7E+vXrrd0NQggheqDghajl7OyM5uZmAMA333yDhQsXWrlH5jV8+HBcu3ZN8TNrKhZJCCHE+ih4IWrJF+2mp6fDw8MDgYGB1u6SWbFYLMydOxc7duwAAIhEIvB4PCv3ihBCiDoUvBC15MHLV199hddee83a3bGIOXPmYOfOnWAYhhLUEUKIDaPghagVGBiImzdvori4GDExMdbujkXweDwMGzYMJ0+epOCFEEJsGAUvRC2BQIBPPvkECxYssHZXLOqFF17Ahg0bKHghhBAbRsELUUsgEMDZ2RmTJ0+2dlcsKjg4GGw2G7dv36bghRBCbBQFL0St3r17Y8+ePV0yKZ0uL7zwAj7//HO4u7tbuyuEEELUuPu+mYheOBwOevXqZe1uWMU999wDsVhMIy+EEGKjDKvYR8hdgMVi4YMPPoCPj4+1u0IIIUQNCl4IUWPWrFnW7gIhhBANaNqIEEIIIXaFghdCCCGE2BUKXgghhBBiVyh4IYQQQohdoeCFEEIIIXaFghdCCCGE2BUKXgghhBBiV8wavNTU1GDevHng8/ng8/mYN28eamtrtV4zf/58sFgspdfw4cPN2U1CCCGE2BGzJql77LHHUFRUhMTERADA//3f/2HevHn466+/tF43ZcoUbNmyRfGey+Was5uEEEIIsSNmC15SU1ORmJiI8+fPY9iwYQCA77//HvHx8UhPT0d0dLTGa3k8HgIDA83VNUIIIYTYMbNNG507dw58Pl8RuADA8OHDwefzcfbsWa3XJiUlwd/fHz179sRzzz2H8vJyjW2FQiHq6+uVXoQQQgjpuswWvJSWlsLf31/luL+/P0pLSzVeN3XqVGzfvh3Hjh3D6tWrcenSJYwbNw5CoVBt+5UrVyrW1PD5fISGhprsZyCEEEKI7TE4ePnwww9VFtR2fCUnJwNoq87bEcMwao/LzZkzB/fddx/69euH6dOn4++//0ZGRgb279+vtv2SJUtQV1eneBUWFhr6IxFCCCHEjhi85uWll17CI488orVNeHg4rl+/jrKyMpVzFRUVCAgI0Pt5AoEAYWFhyMzMVHuex+OBx+PpfT9CCCGE2DeDgxdfX1/4+vrqbBcfH4+6ujpcvHgRQ4cOBQBcuHABdXV1GDFihN7Pq6qqQmFhIQQCgV7tGYYBAFr7QgghhNgR+fe2/HtcK8aMpkyZwsTExDDnzp1jzp07x/Tv35+5//77ldpER0czu3fvZhiGYRoaGpg33niDOXv2LJObm8scP36ciY+PZ4KDg5n6+nq9nllYWMgAoBe96EUvetGLXnb4Kiws1Pldb9Y8L9u3b8crr7yCSZMmAQBmzJiBr7/+WqlNeno66urqAABsNhs3btzATz/9hNraWggEAowdOxY7d+6Eu7u7Xs8MCgpCYWEh3N3dta6tMUZ9fT1CQ0NRWFgIDw8Pk96b/Is+Z8ugz9ly6LO2DPqcLcNcnzPDMGhoaEBQUJDOtiyG0Wd8hgBt/8H4fD7q6urofwwzos/ZMuhzthz6rC2DPmfLsIXPmWobEUIIIcSuUPBCCCGEELtCwYsBeDwePvjgA9qabWb0OVsGfc6WQ5+1ZdDnbBm28DnTmhdCCCGE2BUaeSGEEEKIXaHghRBCCCF2hYIXQgghhNgVCl4IIYQQYlcoeOlg/fr1iIiIgJOTEwYPHoxTp05pbX/ixAkMHjwYTk5OiIyMxMaNGy3UU/tmyOe8e/duTJw4EX5+fvDw8EB8fDwOHjxowd7aL0P/PMudOXMGHA4HAwYMMG8HuwhDP2ehUIh3330XYWFh4PF46N69OzZv3myh3to3Qz/r7du3IzY2Fi4uLhAIBHjqqadQVVVlod7an5MnT2L69OkICgoCi8XCH3/8ofMaq3wPGlSsqIv79ddfGUdHR+b7779nbt++zbz66quMq6srk5+fr7Z9Tk4O4+Liwrz66qvM7du3me+//55xdHRkfv/9dwv33L4Y+jm/+uqrzKeffspcvHiRycjIYJYsWcI4OjoyV65csXDP7Yuhn7NcbW0tExkZyUyaNImJjY21TGftmDGf84wZM5hhw4Yxhw8fZnJzc5kLFy4wZ86csWCv7ZOhn/WpU6cYBwcH5quvvmJycnKYU6dOMX379mVmzZpl4Z7bjwMHDjDvvvsuk5CQwABg9uzZo7W9tb4HKXhpZ+jQocyCBQuUjvXq1Yt555131LZfvHgx06tXL6Vjzz//PDN8+HCz9bErMPRzVqdPnz7MsmXLTN21LsXYz3nOnDnMf//7X+aDDz6g4EUPhn7Of//9N8Pn85mqqipLdK9LMfSz/vzzz5nIyEilY2vXrmVCQkLM1seuRJ/gxVrfgzRt9A+RSITLly8rikjKTZo0CWfPnlV7zblz51TaT548GcnJyRCLxWbrqz0z5nPuSCaToaGhAd7e3uboYpdg7Oe8ZcsWZGdn44MPPjB3F7sEYz7nvXv3Ii4uDp999hmCg4PRs2dPvPnmm2hpabFEl+2WMZ/1iBEjUFRUhAMHDoBhGJSVleH333/HfffdZ4ku3xWs9T1o1qrS9qSyshJSqRQBAQFKxwMCAlBaWqr2mtLSUrXtJRIJKisrIRAIzNZfe2XM59zR6tWr0dTUhP/85z/m6GKXYMznnJmZiXfeeQenTp0Ch0N/NejDmM85JycHp0+fhpOTE/bs2YPKykosXLgQ1dXVtO5FC2M+6xEjRmD79u2YM2cOWltbIZFIMGPGDKxbt84SXb4rWOt7kEZeOmCxWErvGYZROaarvbrjRJmhn7PcL7/8gg8//BA7d+6Ev7+/ubrXZej7OUulUjz22GNYtmwZevbsaanudRmG/HmWyWRgsVjYvn07hg4dimnTpuGLL77A1q1bafRFD4Z81rdv38Yrr7yC999/H5cvX0ZiYiJyc3OxYMECS3T1rmGN70H69eofvr6+YLPZKhF8eXm5SlQpFxgYqLY9h8OBj4+P2fpqz4z5nOV27tyJZ555Brt27cKECRPM2U27Z+jn3NDQgOTkZKSkpOCll14C0PYlyzAMOBwODh06hHHjxlmk7/bEmD/PAoEAwcHB4PP5imO9e/cGwzAoKipCjx49zNpne2XMZ71y5Urcc889eOuttwAAMTExcHV1xahRo/Dxxx/T6LgJWOt7kEZe/sHlcjF48GAcPnxY6fjhw4cxYsQItdfEx8ertD906BDi4uLg6Ohotr7aM2M+Z6BtxGX+/PnYsWMHzVfrwdDP2cPDAzdu3MDVq1cVrwULFiA6OhpXr17FsGHDLNV1u2LMn+d77rkHd+7cQWNjo+JYRkYGHBwcEBISYtb+2jNjPuvm5mY4OCh/zbHZbAD/jg6QzrHa96BZlwPbGfk2vE2bNjG3b99mXnvtNcbV1ZXJy8tjGIZh3nnnHWbevHmK9vItYq+//jpz+/ZtZtOmTbRVWg+Gfs47duxgOBwO88033zAlJSWKV21trbV+BLtg6OfcEe020o+hn3NDQwMTEhLCPPTQQ8ytW7eYEydOMD169GCeffZZa/0IdsPQz3rLli0Mh8Nh1q9fz2RnZzOnT59m4uLimKFDh1rrR7B5DQ0NTEpKCpOSksIAYL744gsmJSVFsR3dVr4HKXjp4JtvvmHCwsIYLpfLDBo0iDlx4oTi3JNPPsnce++9Su2TkpKYgQMHMlwulwkPD2c2bNhg4R7bJ0M+53vvvZcBoPJ68sknLd9xO2Pon+f2KHjRn6Gfc2pqKjNhwgTG2dmZCQkJYRYtWsQ0NzdbuNf2ydDPeu3atUyfPn0YZ2dnRiAQMHPnzmWKioos3Gv7cfz4ca1/39rK9yCLYWjsjBBCCCH2g9a8EEIIIcSuUPBCCCGEELtCwQshhBBC7AoFL4QQQgixKxS8EEIIIcSuUPBCCCGEELtCwQshhBBC7AoFL4QQQgixKxS8EEIIIcSuUPBCCCGEELtCwQshhBBC7AoFL4QQQgixK/8PsFCmGdivWXkAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"y_hat = expected_value(res_nograd.x[:-1], t)\n",
"fig, ax = plt.subplots()\n",
"ax.plot(t, y_hat, label='Estimated Mean')\n",
"ax.plot(t, data, c='k', lw=0.5, label='Data')\n",
"ax.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "bed235f7",
"metadata": {},
"source": [
"# Fancy way -- use autograd from `pytensor`\n",
"\n",
"Instead of writing a numpy function, we write a pytensor computation graph."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "796e917a",
"metadata": {},
"outputs": [],
"source": [
"import pytensor\n",
"import pytensor.tensor as pt\n",
"import pymc as pm"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "eada50df",
"metadata": {},
"outputs": [],
"source": [
"# Create symbolic containers that represent values we want to compute on later\n",
"# The d in \"dscalar\" and \"dvector\" means \"double\", i.e. float64\n",
"params = pt.dvector('params') # represents δ0, δ1, δ2, and σ in a single vector\n",
"t_pt = pt.dvector('t')\n",
"data_pt = pt.dvector('data')\n",
"\n",
"\n",
"# Symbolically compute mu. Note that I'm using pt.exp, not np.exp\n",
"mu = params[0] / (1 + pt.exp(-(params[1] + params[2] * t_pt)))"
]
},
{
"cell_type": "markdown",
"id": "f21598bb",
"metadata": {},
"source": [
"`mu` isn't a number, it's a set of instructions on how to compute that number."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "ce645e0a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"True_div [id A]\n",
" ├─ ExpandDims{axis=0} [id B]\n",
" │ └─ Subtensor{i} [id C]\n",
" │ ├─ params [id D]\n",
" │ └─ 0 [id E]\n",
" └─ Add [id F]\n",
" ├─ ExpandDims{axis=0} [id G]\n",
" │ └─ 1 [id H]\n",
" └─ Exp [id I]\n",
" └─ Neg [id J]\n",
" └─ Add [id K]\n",
" ├─ ExpandDims{axis=0} [id L]\n",
" │ └─ Subtensor{i} [id M]\n",
" │ ├─ params [id D]\n",
" │ └─ 1 [id N]\n",
" └─ Mul [id O]\n",
" ├─ ExpandDims{axis=0} [id P]\n",
" │ └─ Subtensor{i} [id Q]\n",
" │ ├─ params [id D]\n",
" │ └─ 2 [id R]\n",
" └─ t [id S]\n"
]
},
{
"data": {
"text/plain": [
"<ipykernel.iostream.OutStream at 0x105841930>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pytensor.dprint(mu)"
]
},
{
"cell_type": "markdown",
"id": "6464488f",
"metadata": {},
"source": [
"PyMC uses pytensor on the backend, so we can use it to compute the negative logp of the data.\n",
"\n",
"If you were super hardcore, you could also write the function yourself using the definition of the normal PDF. Wikipedia is a good source for probability stuff."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "8df0ab20",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Neg [id A]\n",
" └─ Sum{axes=None} [id B]\n",
" └─ Check{sigma > 0} [id C] 'data_logprob'\n",
" ├─ Sub [id D]\n",
" │ ├─ Sub [id E]\n",
" │ │ ├─ Mul [id F]\n",
" │ │ │ ├─ ExpandDims{axis=0} [id G]\n",
" │ │ │ │ └─ -0.5 [id H]\n",
" │ │ │ └─ Pow [id I]\n",
" │ │ │ ├─ True_div [id J]\n",
" │ │ │ │ ├─ Sub [id K]\n",
" │ │ │ │ │ ├─ data [id L]\n",
" │ │ │ │ │ └─ True_div [id M]\n",
" │ │ │ │ │ ├─ ExpandDims{axis=0} [id N]\n",
" │ │ │ │ │ │ └─ Subtensor{i} [id O]\n",
" │ │ │ │ │ │ ├─ params [id P]\n",
" │ │ │ │ │ │ └─ 0 [id Q]\n",
" │ │ │ │ │ └─ Add [id R]\n",
" │ │ │ │ │ ├─ ExpandDims{axis=0} [id S]\n",
" │ │ │ │ │ │ └─ 1 [id T]\n",
" │ │ │ │ │ └─ Exp [id U]\n",
" │ │ │ │ │ └─ Neg [id V]\n",
" │ │ │ │ │ └─ Add [id W]\n",
" │ │ │ │ │ ├─ ExpandDims{axis=0} [id X]\n",
" │ │ │ │ │ │ └─ Subtensor{i} [id Y]\n",
" │ │ │ │ │ │ ├─ params [id P]\n",
" │ │ │ │ │ │ └─ 1 [id Z]\n",
" │ │ │ │ │ └─ Mul [id BA]\n",
" │ │ │ │ │ ├─ ExpandDims{axis=0} [id BB]\n",
" │ │ │ │ │ │ └─ Subtensor{i} [id BC]\n",
" │ │ │ │ │ │ ├─ params [id P]\n",
" │ │ │ │ │ │ └─ 2 [id BD]\n",
" │ │ │ │ │ └─ t [id BE]\n",
" │ │ │ │ └─ ExpandDims{axis=0} [id BF]\n",
" │ │ │ │ └─ Subtensor{i} [id BG]\n",
" │ │ │ │ ├─ params [id P]\n",
" │ │ │ │ └─ 3 [id BH]\n",
" │ │ │ └─ ExpandDims{axis=0} [id BI]\n",
" │ │ │ └─ 2 [id BJ]\n",
" │ │ └─ ExpandDims{axis=0} [id BK]\n",
" │ │ └─ Log [id BL]\n",
" │ │ └─ Sqrt [id BM]\n",
" │ │ └─ 6.283185307179586 [id BN]\n",
" │ └─ ExpandDims{axis=0} [id BO]\n",
" │ └─ Log [id BP]\n",
" │ └─ Subtensor{i} [id BG]\n",
" │ └─ ···\n",
" └─ All{axes=None} [id BQ]\n",
" └─ MakeVector{dtype='bool'} [id BR]\n",
" └─ All{axes=None} [id BS]\n",
" └─ Gt [id BT]\n",
" ├─ Subtensor{i} [id BG]\n",
" │ └─ ···\n",
" └─ 0 [id BU]\n"
]
},
{
"data": {
"text/plain": [
"<ipykernel.iostream.OutStream at 0x105841930>"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"negative_logp = -pm.logp(pm.Normal.dist(mu=mu, sigma=params[3]), data_pt).sum()\n",
"pytensor.dprint(negative_logp)"
]
},
{
"cell_type": "markdown",
"id": "d758e492",
"metadata": {},
"source": [
"The reason for doing all this is that now we can ask for the gradient of the negative logp with respect to the model parameters"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "244dc5b3",
"metadata": {},
"outputs": [],
"source": [
"# Vector of first derivatives of logp w.r.t parameters, shape is (4,)\n",
"d_logp = pytensor.grad(negative_logp, params)\n",
"\n",
"# Matrix of second derivatives w.r.t parameters, shape is (4, 4)\n",
"hess_logp = pytensor.gradient.jacobian(d_logp, params)"
]
},
{
"cell_type": "markdown",
"id": "2341a74d",
"metadata": {},
"source": [
"Finally, we can compile these graphs into functions to compute on"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "5e338d71",
"metadata": {},
"outputs": [],
"source": [
"f_logp = pytensor.function([params, t_pt, data_pt], negative_logp)"
]
},
{
"cell_type": "markdown",
"id": "ad474c90",
"metadata": {},
"source": [
"Sanity check -- the compiled pytensor function should be the same as the one we wrote by hand"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "01909675",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array(240.78339392)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_logp(initial_guess, t, data)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "0c79af29",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"240.7833939198833"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"logp_data(initial_guess, t, data)"
]
},
{
"cell_type": "markdown",
"id": "67982dce",
"metadata": {},
"source": [
"We can also compile the gradient and hessian functions"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "6fe9f1d7",
"metadata": {},
"outputs": [],
"source": [
"f_grad = pytensor.function([params, t_pt, data_pt], d_logp)\n",
"f_hess = pytensor.function([params, t_pt, data_pt], hess_logp)"
]
},
{
"cell_type": "markdown",
"id": "7cd756d9",
"metadata": {},
"source": [
"Note: What I did was computationally inefficient. We could compile a single function to return all 3 values. Since each graph does similar things (they all compute mu, for example), `pytensor` could figure out how to make them go faster. For this simple example, I don't care.\n",
"\n",
"Anyway, now we have all these functions and we can run the scipy optimizer."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "0bf57c8b",
"metadata": {},
"outputs": [],
"source": [
"res = optimize.minimize(f_logp, \n",
" x0=initial_guess,\n",
" args=(t, data),\n",
" jac=f_grad,\n",
" hess=f_hess,\n",
" # Trust region methods use gradient information to take more careful \n",
" # steps during optimization. NCG is one that can also use hessian information.\n",
" method='trust-ncg')"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "7840f396",
"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>True</th>\n",
" <th>Nograd</th>\n",
" <th>Grad</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>δ0</th>\n",
" <td>3.0</td>\n",
" <td>2.960324</td>\n",
" <td>2.960323</td>\n",
" </tr>\n",
" <tr>\n",
" <th>δ1</th>\n",
" <td>-5.0</td>\n",
" <td>-5.311618</td>\n",
" <td>-5.311626</td>\n",
" </tr>\n",
" <tr>\n",
" <th>δ2</th>\n",
" <td>10.0</td>\n",
" <td>10.759441</td>\n",
" <td>10.759451</td>\n",
" </tr>\n",
" <tr>\n",
" <th>σ</th>\n",
" <td>0.3</td>\n",
" <td>0.305551</td>\n",
" <td>0.305549</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" True Nograd Grad\n",
"δ0 3.0 2.960324 2.960323\n",
"δ1 -5.0 -5.311618 -5.311626\n",
"δ2 10.0 10.759441 10.759451\n",
"σ 0.3 0.305551 0.305549"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.DataFrame(np.c_[np.r_[true_params, true_sigma], res_nograd.x, res.x],\n",
" columns=['True', 'Nograd', 'Grad'],\n",
" index=['δ0', 'δ1', 'δ2', 'σ'])"
]
},
{
"cell_type": "markdown",
"id": "31314f3d",
"metadata": {},
"source": [
"Results are about the same. Like I said, nelder-mead does a good job on low-dimensional problems :)\n",
"\n",
"The advantage here is that we can now compute confidence intervals for the parameters. The standard errors of the parameters will be the square root of the diagonal of the hessian, see here for example:\n",
"\n",
"https://stats.stackexchange.com/questions/381984/standard-error-from-hessian-matrix-when-likelihood-is-used-rather-than-ln-l"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "ba02068a",
"metadata": {},
"outputs": [],
"source": [
"hess_val = f_hess(res.x, t, data) / data.shape[0]\n",
"inv_hess_val = np.linalg.inv(hess_val)\n",
"std = np.sqrt(np.diag(inv_hess_val))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "3abba60b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.7057784 , 4.39726544, 9.47485549, 0.21605585])"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"std"
]
},
{
"cell_type": "markdown",
"id": "1987ddce",
"metadata": {},
"source": [
"To get p-stats, we want to compute the probability of observing values as extreme or more extreme than our estimated parameters. To do this, we can look at the survival function of a normal(0, std), where std is our computed std. Alternatively, we can just scale the parameters by z, then look at the SF of a N(0, 1). As you like.\n",
"\n",
"Here I show that both approaches are the same, in case you doubt me."
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "c472bf7b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[1.3679193610885027e-05,\n",
" 0.8864645274542062,\n",
" 0.12806630935892627,\n",
" 0.0786496011516985]"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z_stat = res.x / std\n",
"[stats.norm(loc=0, scale=s).sf(x) for s, x in zip(std, res.x)]"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "f5688f97",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[1.3679193610885027e-05,\n",
" 0.8864645274542062,\n",
" 0.12806630935892627,\n",
" 0.0786496011516985]"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[stats.norm(loc=0, scale=1).sf(z) for z in z_stat]"
]
},
{
"cell_type": "markdown",
"id": "22b0739e",
"metadata": {},
"source": [
"Evidently, $\\delta_0$ and $\\sigma$ are precisely measured, but $\\delta_1$ and $\\delta_2$ are not."
]
}
],
"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.11.9"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment