Skip to content

Instantly share code, notes, and snippets.

@aflaxman
Created November 13, 2020 05:14
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save aflaxman/b98d3212f59ce24f667225de1fbfafbd to your computer and use it in GitHub Desktop.
Save aflaxman/b98d3212f59ce24f667225de1fbfafbd to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Thu Nov 12 21:09:45 PST 2020\r\n"
]
}
],
"source": [
"import numpy as np, matplotlib.pyplot as plt, pandas as pd\n",
"pd.set_option('display.max_rows', 8)\n",
"!date\n",
"\n",
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Poisson times Beta\n",
"\n",
"https://stats.stackexchange.com/questions/372973/pymc3-mixture-model-with-latent-variables"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
]
}
],
"source": [
"import pymc3 as pm, theano.tensor as tt, patsy\n",
"np.random.seed(12345) # set random seed for reproducibility"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import pymc3 as pm, numpy as np, matplotlib.pyplot as plt, pandas as pd\n",
"from scipy import stats\n",
"N = 100\n",
"\n",
"mu_1_true = 5\n",
"mu_2_true = 10\n",
"alpha_true = 5\n",
"beta_true = 2\n",
"\n",
"X1 = stats.poisson.rvs(mu_1_true, size=N)\n",
"X2 = stats.poisson.rvs(mu_2_true, size=N)\n",
"t = stats.beta.rvs(alpha_true, beta_true, size=N)\n",
"\n",
"df = pd.DataFrame()\n",
"df['Z1'] = X1*t\n",
"df['Z2'] = X2*t"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9ebgcV30m/J7e97vfq32xLFsWBhtbIQYclgQTMIQlhARPFggkDgn5Jvm+yTrzZMgymfkyEyDJQPAYQlhjDAGzBK8QE+PYxsi2vEi2JVmWdKUr3b1v713d1Wf++J1Tdaq6qrd7r+7V7fM+j55Wd1d316069dZ73t9yGOccGhoaGhobF4G13gENDQ0NjdWFJnoNDQ2NDQ5N9BoaGhobHJroNTQ0NDY4NNFraGhobHCE1noHvDA6Osp37dq11ruhoaGhcdHgsccem+Ocj3m9ty6JfteuXTh48OBa74aGhobGRQPG2Cm/97R1o6GhobHBoYleQ0NDY4NDE72GhobGBocmeg0NDY0NDk30GhoaGhscmug1NDQ0Njg00WtoaGhscLQlesbYZxhjM4yxZ5TXbmeMHRL/TjLGDvl89iRj7GmxnU6MXw/gHHjiS0CtstZ7oqGhcYHQScHUZwF8HMDn5Quc81+Q/2eMfQTAUovPv55zPtfrDmqsMGafA775W0A0Bex/+1rvjYaGxgVAW6LnnD/AGNvl9R5jjAH4eQA/ubK7pbFqMEr0WCuv7X5oaGhcMCzXo/8JANOc82M+73MA9zLGHmOM3dzqixhjNzPGDjLGDs7Ozi5ztzR8UReWTb26tvuhoaFxwbBcor8JwG0t3n815/waAG8G8CHG2Gv8NuSc38o5P8A5PzA25tmXR2MlYAqC10SvodE36JnoGWMhAD8L4Ha/bTjnU+JxBsAdAF7R6+9prBAkwZua6DU0+gXLUfRvAPAc5/yM15uMsSRjLC3/D+CNAJ7x2lbjAqKuFb2GRr+hk/TK2wA8DOByxtgZxtgHxFvvgcu2YYxtYYzdKZ5OAHiQMfYkgEcBfIdzfvfK7bpGTzANetREr6HRN+gk6+Ymn9ff5/HaFIAbxf9PALhqmfunsdKQwVht3Who9A10ZWy/wbJujLXdDw0NjQsGTfT9BovodWWshka/QBN9v0FaNqZW9Boa/QJN9P2Gug7Gamj0GzTR9xusylht3Who9As00fcbpGWjrRsNjb6BJvp+g+51o6HRd9BE32/QHr2GRt9BE32/wdS9bjQ0+g2a6PsNlnWjPXoNjX6BJvp+gyR4reg1NPoGmuj7DToYq6HRd9BE32/Q3Ss1NPoOmuj7Dbp7pYZG30ETfb9Bp1dqaPQdNNH3G9Q1Yzlf233R0NC4INBE32+wlDwHGvU13RUNDY0LA030/QbVstGNzTQ0+gKa6PsN9SoAJv6vi6Y0NPoBnSwO/hnG2Axj7BnltT9ljJ1ljB0S/270+eybGGPPM8aOM8b+aCV3XKNHmFUglrH/r6GhseHRiaL/LIA3ebz+Mc751eLfne43GWNBAJ8A8GYA+wHcxBjbv5yd1VgmOCdFHx2g59q60dDoC7Qles75AwAWevjuVwA4zjk/wTk3AHwZwNt7+B6NlYJZA8CBaJqea+tGQ6MvsByP/rcZY08Ja2fI4/2tACaV52fEa55gjN3MGDvIGDs4Ozu7jN3S8IW0arR1o6HRV+iV6D8JYA+AqwGcA/ARj22Yx2u+iduc81s55wc45wfGxsZ63C2NlpAZN1rRa2j0FXoies75NOfc5Jw3AHwKZNO4cQbAduX5NgBTvfyexgrBInqh6LVHr6HRF+iJ6Bljm5Wn7wTwjMdmPwKwlzG2mzEWAfAeAN/q5fc0VgiS2LV1o6HRVwi124AxdhuA1wEYZYydAfBhAK9jjF0NsmJOAvgNse0WAJ/mnN/IOa8zxn4bwD0AggA+wzk/vCp/hUZnkJ0rLUWvrRsNjX5AW6LnnN/k8fI/+Gw7BeBG5fmdAJpSLzXWCHVXMFZbNxoafQFdGdtPcHv0plb0Ghr9AE30/QTTHYzVHr2GRj9AE30/ocm60USvodEP0ETfT3Dn0eusGw2NvoAm+n6Ctm40NPoSmuj7CZaiTzmfa2hobGhoou8nSGIPxYBgVFs3Ghp9Ak30/QRJ9MEoEIrqgikNjT6BJvp+glTwIUH0WtFraPQFNNH3E+oK0Qej2qPX0OgTaKLvJ9SrAAsAgRAQimii19DoE2ii7yeYVVLyjFFAVls3Ghp9AU30/YR6lZQ8AAS1otfQ6Bdoou8n1Kuk5AGRdaOJXkOjH6CJvp9QF9YNQIped6/U0OgLaKLvJ5hVUvIAKXvdj15Doy+gib6fUDcUotcFUxoa/QJN9P2EesUm+mBEZ91oaPQJNNH3E0zD9uhDMR2M1dDoE7QlesbYZxhjM4yxZ5TX/hdj7DnG2FOMsTsYY4M+nz3JGHuaMXaIMXZwJXdcowfUK3Z6pS6Y0tDoG3Si6D8L4E2u1+4DcCXn/GUAjgL44xaffz3n/GrO+YHedlFjxaCmV+rulRoafYO2RM85fwDAguu1eznndfH0EQDbVmHfNFYa9Sp588DaB2M5B/LTa/f7Ghp9hJXw6N8P4C6f9ziAexljjzHGbm71JYyxmxljBxljB2dnZ1dgtzSaYLoLpipEuGuBF74HfGw/kDu3Nr+vodFHWBbRM8b+C4A6gC/5bPJqzvk1AN4M4EOMsdf4fRfn/FbO+QHO+YGxsbHl7JaGH+qG0gIhCoADjXrLj6wacufot/NTa/P7Ghp9hJ6JnjH2XgBvBfCLnHvLQs75lHicAXAHgFf0+nsaK4B6Rcm6EYS/VgFZWaxVzq7N72to9BF6InrG2JsA/CGAt3HOSz7bJBljafl/AG8E8IzXthoXCKahWDficc2IXvxuZWltfl9Do4/QSXrlbQAeBnA5Y+wMY+wDAD4OIA3gPpE6eYvYdgtj7E7x0QkADzLGngTwKIDvcM7vXpW/QqMzqOmVMii7Vpk3UtFXtKLX0FhthNptwDm/yePlf/DZdgrAjeL/JwBctay901g5NBrkiavBWEAreg2NPoCujO0XSOWuplcCa0f0cn+0R6+hserQRN8vkFaJWjAFrKF1oxW9hsaFgib6foEsjgq5Ff0aFU1ZHr0meg2N1YYm+n6BJNag0qZYff2C749U9Nq60dBYbWii7xfI1aTWjXWjFb2GxoWCJvp+geXRR5yPa2bdiN+9GIKxZh04es/atYvQ0FgmNNFvJBTngEdu8Sakuo+iXzPr5iJS9Me/C/zTzwMzR9Z6TzQ0eoIm+o2Ew3cAd/8hkD3d/J5feuVaLRCuevTrXSnLOEJxbm33Q0OjR2ii30iQ6thLJVvWjTsYu8YefaMO1Dy7aKwfGEV6vBhmHxoaHtBEv5FQzTkfVVjWjVwzdp0oemD9E2itTI9ex1VD4yKAJvqNhIok+nzze+suvbICsCD9f70HZOWMo6KJXuPihCb6jQSpOL0IyZ1eudbWjVkFUuP0/3Wv6CXRr/P91GjGmYPAM19f671Yc2ii30iotLJuBKE3da9cQ+vGIvp1rugNQfTaurn48MgngbtbLWndH9BEv5FgKfpWwVih6BkjG2ctrZvUBP1/vStlbd1cvKiVgdIcdW/tY2ii30hopeilcpdKHljbBcLrVZvoLxqPfp3vp0Yz6mXK7Orzc6eJfiOhlUfvTq8EiPTXogUC57Q/SbE28LpX9Drr5qJFTYz70vza7scaQxP9RoLMtmmVXhlUiD4UXZtgrFmjx0gSiKTXP9FbefSa6C861MVNuji7tvuxxtBEv1HQaNhE75l1UwUCYSCgnPK1Ino1XhAbWP/Taqno1/sNSaMZUtFrotfYEDDyAEQrAb+sGxmIlQhG18a6sTKAooLo1zmB1nTWzUULregBdLY4+GcYYzOMsWeU14YZY/cxxo6JxyGfz76JMfY8Y+w4Y+yPVnLHNVxQVbynR1+1UyslQpG1V/TxwYsoGJtb/315NJywFL326NvhswDe5HrtjwB8j3O+F8D3xHMHGGNBAJ8A8GYA+wHcxBjbv6y9XQkc/Efgh/9nrfdi5SHVZiTlr+hVfx4gol0Tor/IFL3Mo2/UbBunGxy+A5g71nqbJ28HHv5E99+t0Rpa0QPogOg55w8AWHC9/HYAnxP//xyAd3h89BUAjnPOT3DODQBfFp9bWzx1O/DEF9Z6L1YeUsUPbPP36EMuog9G1qZgSs0Aig1eHB69vEn2Yt/c8ZvAI3/fepsffAS45z8DR77Z/fdr+EN79AB69+gnOOfnAEA8jntssxXApPL8jHjNE4yxmxljBxljB2dnV/GkGMX1bxX0gqpC9EYeaJjO9+uVZqJfs6wbqehjF4eirxWBtCzu6pLoa2VSlfnz/tuUs8Dc80AgBHzzt4GFE73vq4aNRsMea33eYno1g7HM4zVfg5Nzfivn/ADn/MDY2Njq7VWtBJQXV+/71wqqogeaG5vVDQ+iXwfWTXyQblLuG9N6gVmjgpv0Znre7U1JjrVWRD/1OD2+5aNUsfzV99lKVKN3qFXfJU30vWCaMbYZAMTjjMc2ZwBsV55vAzDV4++tHIwSYBTWriJ0tSAVfWab87mE6eHRr1XBlDu9Eriwqv7MY+SJdwKZQy+reKtd7mdJuJ6F6db7AwD73w684xbg3JPAd/+0u9/RaIZF9ExbNz1+7lsA3iv+/14AXsbijwDsZYztZoxFALxHfG5tUZPFLxvMvqm6FL3bYqh7ePRr1QLBHYwFLizR/+hTwD0dNrqSwVdL0Xdp3UhFX5j277dy9iAwehnNbvbdCOx/BwVwNZpwNltGsVrvbGPr3G2iG67Z4efWCs98HfjX/7YqX91JeuVtAB4GcDlj7Axj7AMA/n8ANzDGjgG4QTwHY2wLY+xOAOCc1wH8NoB7ADwL4Cuc88Or8ld0A5lBUXLHly9yVHLk8cqOkG5F70X0wcjaNDVTe+PHBun/F/LGW82TL95JqqRMrUxvosderZtGHSh7jDnOqZXuth+zXxu5FCjOrF87aw3x7k8+hL///vHONpZEP7AdAPc+/iuNhz5OM7Je8MK/Ak98cWX3RyDUbgPO+U0+b/2Ux7ZTAG5Unt8J4M6e9245OP80kBy3g2iA8FtF+f1G8+mrOSCaURSyB9F7pVeuiXWjrHa1ForeKADcpMdouvW2FtELRd9t1o06zvLngeSo8/3Fk+Qfb73Wfi29CeAN6s+S8spzWGE8+ilg1/XA+BWr/1vLxEy+itl8h2NWplYObgfOPEr2zWoeT86B+/4E+PHfBDZf1f3njSK1BVkFbMzK2BPfB259PfD9/+58XfqtQG9EP/8C8NVf7S2XerVRyQGxDJE90ByM9UqvDEXWyLpxFUwBFzYTSo6DTsaAnAEmxwAW6N26AYCCR0D2rPDntx2wX5PxgFYBXD8snQU+97bOC4RKC8Cdv7dqSnIlYdQbqDc4ikaHMx0Z0JZ25mpn3hhFukH3ugZyrQSEEyu7TwIbj+jPPw3c/suk3N32jHoCeiH6Fx8ADn8dOPv48vZxNVDNkTqNSaJ3KWSvrBvZj97Lwjh8B93YVgN1V3olcGEVfbVAj52MATlmIgm6ifZq3QBA3iMge+YgEIoD4y+xX5NE3yqA64cT3wde/DfgfIf2gRzL6z3FFUBZEHy5U6K3FP0OelztgKyc7fUqBI0iFTyuAjYW0WcngS+9mw7W0C6nggdsdQb0RvTyRM4+2/MurhoqOSA6YCv6JuvGK48+BoCTf6yCc+DrN1MRz2rAXTAFXFiP3lL0HfymJPpwgm6ivVg38mbmqegPAlteDgQVFzW9DEU/d1T8bofHc+riIfqiQeO082CsVPSS6LtU9LUyee5eazB7QW5XK7bezg9GkQTFKmDjEH05C3zxXUTmv/Q1uou7iV49Ab0EZiR5zjzX+36uFqrCugnHKSjbSXql7H3jzqUvL1LF7PmnVmdfVUUfSdIi4RfUoxcXpPtm//Q/A399mdPOkuosnBDFXT0QfXoL3YTdir5epcDdtmudry9H0ctWC50KmU4V/bH7gKP3dr8/K4iSVPS1LhV9ehONsW4V/al/B+79L8DXfr2zFaosol+OotcefWtEUsAlrwXe8yVgYj/1OV8tRT+zDhW9DMYyJiwGr6wbV1MzSfxuopcXxMxzq+Ph1yvkdwdDtL8XujrW8EmxPf80kas6NuS2kQSRdS/WTXyIVHr+nOv3nqEb6tYDztfDcfqtXoh+XhB9JzMkzu0YQbuZyvf+vPOU1FWCtGy6VvSRJJAY6b5oSp7ro3cB9/9l++2Xa93USkBYE31rBEPAjf8L2P0T9DySpKwKFcv16Cvr3LqR/rzbYuDcu02xtHLcmTcFUf/WqFFp/krDdO3LhexgWTfs/j7uMSCJQCXzlbBu4kOk0t3EffYgPaqplRLpie6tG7Nmt0/o5HjmzlIaJ1j7G1j2NDB/fE27QErrpmuPPhSjYHq31o283i97E/CDvwae+Vrr7aWidwvMTmEUtHXTNVoRfSi+PEVfmgcK66jSjnMaZNKfdyt6swaAe1g3fopeKXQ+twr2jTun/0IqenVMuMeAJDGVzB1E36N1Ex8i+8BN3GcOUtrmgEcLqNSEfcPtFIsn7XhLJ4pe2jabX9b6+Fdy9ved+VF3+9Qp6lXq87N01ncTS9F3m3UTjlNaa7fWjSTud3wS2PFK4Bsfat2FdNnWTUlbN10jkvS3bga29q7oWZD+P3Nkefu3kjCKlBduKfoBJ1mZSiWqCrlQuLuDpbyJsSDZGSuNesWp6C9kB0sH0bt+U64rqu6LUQLA6Nj1mnUTH7QVvZrhNPUEsOUa78+lN3kHb1vBIiHWmaI/+xitOrbz+ta99peU3oSTP+xunzrF3FHqKnv8Pt9NSr1m3YRiPRJ9jizG+BDw7s/SDPfQl/y3lyKgl/RKs07XqbZuukQkRQdcrS6UwdjMVqDUi6JfIvUDALPrKCArSV0W/7gVvVqgpMJS9K7q2OIMkfyWl68S0a+lom9RS+Fp3ZRFwJjZ1k0ngTmAFGWtZCv6esX+7loZWHgB2HSl92dTExS87WahE5lxM35FZzfOqceBiZdQEVGrXvvZ0/QYigOTj3a+P91AnpcWdlVJWDeG2UDN7OAc1JR6jV6sm2qerinG6Pztfg21kfY7J5ai74HoJTdpRd8loiIfVT3olqLf3ruiH72MFOh6CshKUo/6ePRWywH3ClNCVbsDroUZujA2X0VEv9KrKtUrThvpQnr0VUXRu8nQUvQq0Rdp6g+INEnebAn6QX5/fAhIiRYK0qeffY6KayZe4v3Z1AQp0m5iAnPH6HODO9qP70YDmDoEbL3Gngn63Wwl0e+7kW4OcnH3lYQ8pjn/voclRcmXOlH1dbGOQCBAir6a665bayVnX1MA9SBaOOEvfpYTjFWD/quAjUv08s6oKriaYt0Y+e4HrMxsGb/iwij6o/d01mVRDjCZr93k0SvpjCos68Yj6yY1Bmx6Kc1isqe63/dWcBdvrYVHHxt0kqFZs/fBrehltWK0DSG6Ib9fKnrAVqzTou3TuA/Ry+278ennj9lCpNxmH+eP07jZem37orXsaRo7l99I19D0M97bLQcdKXqV6DvIvKlVgLAY8wnReqIbVV91Ef2+t9JM129xGFXRdzrrk5AiVBdMdQl5wFQFZxSJ3JKi3303KpJzO7NlbB959Ku9fujDnwAe+J/tt/NT9HKwWXnrbkXvY90UZqhP0CZhU620fdPk0Q/QzeZC9GCXRD+wzXn+S0o2iUp4RtEmeqvquEOVLYk+MawQt1D004fJChne7f1Z2ZOl08wbzoHZ56khWnyovXUjC6W2XGMTvd/flT1Ns4Qd19HzyVUIyFpEf853k7JC7h0r+pCYjclrvhufXlabSyRHqCfQkW94X/tqYVW3zQLluNQtELqEpegVoq+VaBoeF2uZd2Pf1EoU8IxmgPH9RAa9VC52g+JcZwPTUvRK1o1qMdR9FL1F9C7rRjZ/Gr+CglGdEP2z36bmWG7MPAs87FpGr8mjv4DVsZJQOiX6WtmeTvs1jPODqujd/WumDwPj+4BA0PuzbqunHUrzdPxku+N2i7mcfZwCf2OXK8e/haIf3EHHLLN1dQKy3Sr6agdEXyvbit4i+i4UvZqyLLH/7TQb8rJuVaLv1qe3Wm1oj747eFk3hihIsBppdVEdW1HIdHwf/b+TfPrZo72v4lSao4uvXdFS1UPRq6/LrJqmhUc88ug5tz36SAIY2due6KePAF/7NVrz1E2C3/8fVGij/g1eih64MPaNvBjlkovSvlMJQP0b1EZT0S73U/Zaig+RMgwnnIrez58H7DYInRK9zLgZ3dueuAHKuNlyNd1o2llS2dNYimzCsek85fyvRkBWXqfFWV9Ltdi1daPYbrJraDdFUzIYq+KKnwHAvO2b5RC9oYOxvSEiTpDDoxe9JHpR9CqZjol2ru0CsvnzwCdfCTx5W+e/I9Fo2OTTbnBWvBS98rq8gN2DyCuPvpoj4pfWweaXtc6lr5WBr32AlL9pUFxBwihS6Tzgmta6FP2F7GCpKnrAPjZS0YcTzQVTy7Vu4kOUuZESRVCFGTqnEz4ZNwCRdTDa+axRZtyM7lWOp8/4Nut0897ycvFb8gbmcfyrBaC8gO9MhvG7tx8Ctv84sHQayPlbLD3Buk65b1yie+tGERSS6Lu2blyKPjUO7Hw12Tde20sYmugvDLysG0NctPFhet4N0VeUgGdqjII77Yj+5INUwNLLRVHJklUEtB+c1RwAZufguglJ7ufoZc7PeXn0Moc+KYh+00uB3Bn/hVru+68Ur/j5z1Pxj3oBHP+urWzcWUCe1s0FUPRGAQCj/jOAPQYk0Q9f4vLoS66smy72s7xIfYdkvCi9iRS6DGaO7/f/rLwxdKzojxKpDWy3j6ffjbO8SDfzoV30vJUlJXLoT9ZH8Oy5HIoTIu//zAqrelWQ+dzcSoaJYIBZ/2+LWsU+d9EMxee6Ivp8s3UDAC95ByVjuHteVfM2t/Sq6LVH3yX8PPpIsjdFLy9ueYdXMm94ZQmH7vw0zJrLYjn5oPOz3UD1jNsNTpkGFhCnM+q6cGeOkOebHHF+LiGeq822ZFVsSniam15Kj172zfN3A4/eClz3IWDvDTStPf5dOwCuTm8dBVzG2lk3shVswnWzl7OnoV3eefRAb1k3Us0DQtGfszNuWlk3QHdtEOaPA8N7yIqR47viM77ddRfhGJGg198lUitfqI2gwYHHq9tpprHS9o16nfoEZEuGiZEkJRQUO7Fu6mV7nDHWXS593SBB4rUwzd4b6PH0Q87Xqzk7FtNtimVNZ930Bk+PXmRQRDNkNXRl3YiLQN7hx/bRHf3JL6P+t9fi6kf/E47cfavzM6f+nR57ITCV3Nu1W6i6gkZuRT99mBq9uRGOkwKcV5Zmk9NmS9G3yLx55O+Bod3AGz5Mz/e/nS6O4/fRQD96D5EP0JyR4E6vBFYmGDvzLLWr9kM1TzUW1s1e/GZpnl5LDPvn0UtC7Ma6kb8DiDYI0xTTSE00rzblRjdtEOaOkm0DtLfC5LlQScwvxVUQ/XGDbow/miyQ5bPSrRCMom23+hJ9HaMpGjdlwwQe+t9Oq9ANVdEDojq2Q6K3jtFA83tytTG1749sQyItz64VvbjR6Tz6LiHvjA5FLzIoAgGa3nazbqw7hXH8Cgrm3fEbKCe24BwfRvLEXfb2hRnbN+2J6JUB2amil4gqRG/WKe3OTz2O7LE7Hqq/JQdscpRmA142VWGaKjslae94JammI9+k9S+NAvDyX7L3UcLdYK1bpdwKX30fLefmB9kKNubysUtzNMNxtzlQA3pAd/1u3ESfmqAxc+bR9mpebt9JG4R6lfrcSGuuXRZTV0R/CjwUwymDrqeDpxYpA2elM86MIjC0k/LUW1g3o2kaa4GlU8C9fwI89jn/71QVPUB2a6fWjRR2Xoo+FKWbkjrrrpWoAE6m0XZN9LLVRrztpr2gZ6JnjF3OGDuk/Msxxn7Xtc3rGGNLyjb/dfm73CFCEerj4Q7GSh87PtRbMFaq5Ut/igJTb/s4/v21t+FfzOuwPftD+yKSaj6SXr6i78Sj91L0lRyV2ZtV/8Kckb20kpTMCy7M0Gwnodg8mS3ehCPz7SUCQSoqOXov8OSX6RhfLpYQbqXowzG6IJdL9JyTAm0VEzEKJAIse0OQYXGOiCA2SARRr1J6Yr3iJPpoFx0svRQ9QDOoTog+vYm+o13W1sKLRDJNit7PuvEheq+/KzuJRnorAIZoKIBDk1mY0XT3XTzbQa7f69X8TaBcMzGcCIMxYN/kVwDw1teGWjAFdGfduBMc3Ei62h7LYyoFUi/B2HDCtl9XGD1/K+f8ec751ZzzqwFcC6AE4A6PTX8gt+Oc/3mvv9cToqnm9Eo5NUoMdx+MZQF7pjC0C/jAvcA1v4xijeMe8wDCvAYcE4sznHqIbio7rlueR58cbz843dkB4QQpo2quvR88ciltJy2C4gyRvJrfnRpvDgqadUpPlfnJEvvfTjfUZ78F7HuL7YVLYmg0mj16YGWqY6s5UlKtspSkRy/tIkvRL9DsRQ1MSp9VnU7HumhsVs7awTnA9m8B/xuvCmsBkjb2jZw5jlxKj6Go6NDaTtG7ZoE+1o2RpgylV+4ZQckwMV+LtW6C1gsk0aU3+Vo3xaqJeCSEkYiJK6dF/KfY4tioBVMAzV5zZ+jG2A5eN0MViVFX7YUY36keFX1t9RYdAVbOuvkpAC9wzle4Vn6ZiKSclbFqqlwvil42OHKhWK3jcX4ZZvkAGke+TS+e/Hdg+ytEj40erZvoALVr6MS6UZWHbMBVEUTPglQY44VRQQ7Spy/MOlU6QGTujhNIMk25iH7X9baK3f8Op40EKDn9rirdlSB6qeRb3RilRx8M0b45rJthZ5Wo2qLYsZ+dKvoFb0UPdG7dAO0zb2QP+pE99mutqmPdwVigpUdfTlAb5Z/cR+PiVDFEGWG9tuP1grTU0pv9Fb1RRyISxNtCDyNu5qkXU6tz7Vb0L/9luhYe/kT7/XHXpriRGHH+trwxWNZNl8dmFZcRBDNureYAACAASURBVFaO6N8DwC9Z/JWMsScZY3cxxnxHN2PsZsbYQcbYwdnZFer1rvakbzTsrBuge6KXa7J6oFCto4EA7jOvJUWfOwfMHAZ2vbp3AivO0k0iOdaZdeMekNJimDlCSs/duVJixEX0xZlm8k5N0D6o/TusoK1r22AYeMk7SfHsfq2S0SEuHGu92FVQ9FIJVrL+fYzU5dpkMzXOSZ0lRp2BYS+i79S6qRs09hxEL4J4rW68KjpdO3bxJM0cYsr4bNUorlOP3igCpTnkY5SKunc8ja2DcbywJMTOSto3cqblo+g55yjVTCTDAdzE78JUZDfNHo2Ct03CebOiz2wGXvYLwBNfbL+ASjtFn3Qpenkseg7GllYt4wZYAaJnjEUAvA3AVz3efhzATs75VQD+NwCPKgMC5/xWzvkBzvmBsbExv826g9qTXvSm/tcTBXz4m8/0puh9/DpZpXd348cQqBXtZcd2Xm8rwG6bHJXmOiN62YPHPSBlY7PpZ1qrxwGRLicDsm7fHaDBy01nJbHcJ/e2APDT/x34zYfs3jrRtH3hWO0YXDeeFSF6hRBLPheyJBTAbmxWWaJ6h8SIM9XTajTVg3Vjda4ctF+LD9FNb3Sv/41XhdUGoR3Rv9jcMyc26L+f1Tzl9ztSXD2WnxTZS9kI7UcmHsKBXUM4LIdBt4uwtIJREIp+E40zV1yiWm+Ac2B35TD2Nl7Ed9Nvt8ee1/UhPx92CYpX/T/EBT/yaNehQq2b8UJimMaYtK/k+JbnuJesm1XKoQdWRtG/GcDjnPOm+SXnPMc5L4j/3wkgzBhrk1O2glCJXly0xxYbePSkCJJVc513sHRntigoVk3EwgE8wl+CajBFiiEUEy1gZWvbDleSt750jkheLpjg54fWK9RL3H0TimVoqbjsae/USolAkIqEZEBW9rlRIZ+rXrE7O0dFOG6rUcCpgldV0Sstbv2m9DIYC9j2hhUPGXUSvbowuLWfg50RnFoVK8EYMLjTrkhth+QYANbeo194kdJcVcQHWwdj3TZkbMAOQkuIYqm5EJ3LTCyMAzuHcKYcFt+z0oo+ac96XLMYuU7s1ee+giJL4v7I6+2x50n0ctERVxbL+D5aGvDRW1sHTL3sLRWJURrLkl/UGUA43vzdd/0R8JVf8U8UUN2GVcBKEP1N8LFtGGObGKPRxBh7hfi9C7fopLpAuGjsv1gLI1+p2UGyTsmluuSr6IvVOoYSEWwezuDJxHUAOPUECUV7LwYqinS/5Bj52n4XlTvtUyKaUSow2/jBo5dSrxSjQAPObcdI5aR6xZai7+C+HU0r1s0FUvReF3+j4WHdLNpEn3ATvUe1YjRDr5ttCna8iB4AfulrNOPpBMEQHd9W1o1ZA5bO2FWuErE21o2bwKyUTGWcifbU0wE6/5lYGNfuHEaeC/JcqQK3ukFiRVo3QJN9Iytht80+gEcTr8FiPdy6rYF1k/ZIV3z179A5b7VaVDVHM12/mZfMSpNjRw1wh5PNiv7oXZR2/PfXAU99pVm4Ges4GMsYSwC4AcDXldc+yBj7oHj6cwCeYYw9CeDvALyH89Xu7asgkrSVtLjDLtRCyFfq3VfHtlL0Rh3JaAh7x9O4q36AXtx1PT32QvSNBg2g5KgyPfVRqO5e9BKxDKXcAe0DfyOX0vRfqg0/Re8o4pqhC8EvWOXYlwH7QvBb1lBaDcsZHrkpO8jrZd3USgC4vSiNtO/ksU0MOxfh8FT0Hfa78SP6oZ12JlInSG1qHYzNniZbzW3dtAzG5pvPm1er4uxpIBjB+QbdBFKxEC7flEYj0mXPn3ZQV1eyFL2T6Ms1E3FUEDZLyEa3kl2a9JhpWt9J564RijWvRrXjlcDWA8DDH/dPXfVrfyDhbpKmzgDC8eZgbGUJ2PvTVOvw9V8HvudKQFzPRM85L3HORzjnS8prt3DObxH//zjn/CWc86s459dxzh/y/7ZVgGrdiDvsfDWEQrUO3i7X2I0WHn2xaiIZCeKyiRRuX7oC5rXvB666id7spRhI9rmR1g3g79O3UvQAzWoGd7T+vZG95FGfPUjPvTx6wGXdCGvJIwupCWretV/L5NhA6+XsOkH+PBZTIrjsdbzcjaPiQ6R65cWaHCVVyQLCo/dY9afTG7cf0XeLdJt+N4sn6dHLujEK3taku886oIxT5eaQPQ0MbEe+aiIdDSEYYAgGGHZsEbZctUs70g+GF9E3WzcjjH6vEh1BsWoqrYc9iF5YhN94ZgG/8H8edr7HGPC6P6Zj9+DfeO+TV9xLhVzIRBZdVnNkEwXDJAxURS/jaJuuBN5/N7DjVdQqRIW67sEqYONWxgIuj54eizwCs8FRCYmB3Ul1rDhRPJrBUrn5wilWhaKfSKHUCOPEK/6clBvQm6K3FOZo+wUT3IVcEvL5+BXtyVhm3sgiL3fWTTRD6r3g6onj3s4Pnh69h3UDLM8OyJ/HDwvjMBHwngFZZebiAo4N0s1FtkxIjIjU1AFnHr06/Xeni/phpYg+NQHMnwBeuN97trMocsK9grGA9/H0tG48jn92Ehjcjly5jnQsZL0cT8tisxVS9CrRy2CmW9EbJkZA+1aLjaBcMynQGs14n2uxiM2ZAsczZ3NoNFzHbu8bgJf8LPCDv6bKcTe8MtlUyFmZ/G31mEZcRG8USbjFBkRMbHcz71wEHv36hbpAuDjwJU4EUwiIk9iJoheLjpzIB/Bjf/ldzBWc071C1bZuAODotJK73xPRK/63JHq/gJxfvq983km+tqyoPCWUj1vRM0aq3m3duL18Pzg8+hbBWKB3om80gMJ5nDEHkQ8MeBdNWUSvKHqAMo5Ccft1GS+w0iuVC9AqtGrTl6e0QDODTqytVrjmV4jQvvAO4NbX2m2fJRZeBA/F8LmnK3C4oq1mrJ0S/cIJYGg3cpUaMvGw9XIoLj67UtaNdV5S9kLcLkVfMkyMMto3Mz5qBWeRHPW+NkQwNlcLwTAbmC14WDRv/itS0d/+neasOK9jpMKybuabt3cHY+Uxlcc4MezMYHOnfq8CNjbRqwuEC9VQAhF9jon3OiF6QVLnq1EY9QbOLzmXCSsZZN1cOp4CY8DRaWVK2wmBPfcd4MxjyhcqVoJl3fh49BXFG1QR64LoE8NEegsv2L/rhrs6tjjnnVrphViGLgTOWwdjgc6JfvGU0+YpzQGNOk7VBrEIH5Uni+csj16Q4fxx59/cRPSKou+0iEm2P1huSfuO64DfeQr4mb+lc33be5w3mcWTyMe24MPffhbPnlPHXYvGZp5E71rDoLxIZDSyB7lyDZmYTfSJaBQFHl8dRQ+IoilXMLZmYoTR7zUSo6jWGzAbXFSOewVj6RrN1qjC+8yiR4ZNahz46b8ETj8MPP5Z53uVnH9qJSC6xYYVj14l+qRzbLqJPj4sOEns0yqvLgVsdKKXB65asA58WSj6bCMBgNlE/+DHgC//ovf3COWyYIqbRMVp30jrJhYOYsdwAsdnFEXfzqOv5IB//gDwr3+hfKEk+jHy/OJD3Vs3Uq12QvQA+fQADcJguPn95LhdHSvTMDvIuLnt0dP45MMzopKy1ELRd7Gc4MyzwMcPAD/4qP1ajlIrZ/gQ5nnax7qRhKIEYwFKLVUDpJLo5YWoeqeZzY7f84W7z81yEI4B174PeNvfUSxFbRG88CJycWpRMJ1XBIi7l4+KThT9vKi2Hd6DfKWOTNy2bhKREHKIo7FSWTdNRO+h6Kt1jEKMddFuu2TUyT5skV6ZNYjiziz6xH6u/kVg92uA+z7sVOHtFD1jzupYt6KveSh6yQUyY0eq+lXuRQ9seKKXHSyLtnUjFH2+aor0ugXqS/PdPwOe+xfvAJNQLnM1Una5sjO1TmbdAFQ96FD0wVDrxmZHvkmD8vzTtgdrefRiQLQqmspOitbLLvWx96eBt34M2H6d9+fckD69V148IC4oMUWuZMnb9ttWwUMvzCt513l7SUGvFghAe0Vv1oFv/CalnJ5VZkGCGKb5EGYaaR/rRpwXN9HXSnZwTe6LVPShmFOVR9N0wV5IopfYeoAKnWQfdM6BxZNYjFKLgtm8Yk/4tSo2azTe3JaSGoQG7NndyB7kKjWkFUWfjAZR4HGY5ZVW9HKBluY2CCWDFD0PJxEV1lHZMP2vDaHoFy1F70P0jAHX/ioJJhnvACidup3tlhixvXZ1BuBH9FLMuFMza66/fxWwwYleWXzEsm5ISVoploungK//hq1i3avGAFavmtlas6Kvmw1Uag0kI0T0l02k8OJcEUZd8fxa5Ygf+id6LM3ZdkBpjj4j96lV172ZI9Qb320RRBLAgfd3bh3IPil+vrvVBsFsXoWqBY5N55HnQqlUcsv36B/6O2DqCWBgh10nAFhT/fN8CDNmGryTrBt54QHObp0q0XuprMwWZ3GWF1aD6CMJYPPVwOlH6HlxFqgVMRfxIHq/GZJfab8VhJaK/gUi/qFdwrpxKvo8EmiUV0rRu2In6c1EvEqfqnLNxAhbAlLjSESIvIuGSWOwtNBc1yAU/Xy1DdEDdlaaDMrL3vKt0isBZwdLRzDWlUfvToG2iN6l6C+CXjfrE+riI7USOBiqIPIkoh+mRTJyZ4F33kLbquQhIT36SsT+rECpRoUcySgNqL0TKdQbHKfmla6ZfkS/cILU2aVvoOfnxW8XZ52EK6tjvTBzpHXla6eQAVk/lZ4cp7z80kLHxVJ1s4ETc0XkITzuaq6FR++R3ufGzHO02Pj+twPXfZBujPKmkz8HDoY5DGCBZ8C8FlVv8ugVIlb/lqgk+rI30ac3r42iB4Cdr6SZTK1idWGcCZGdNJNTrRufYGyrHi5qdtTCC8DANjQCERSqdUcwNhkNIs8T4C6PnjcaaJgdLPHnhpdHDzjiICWjjlGWA0uOISFEFeXSjwLgzXUTQtHn6rStp0cvMbCdHkUlMGXJNFpbN4Czg6WastpJMBawP2toj355UBcIN0qoBeIAKNWQqmPFhfjaP6BUq0iaiNMNMfinqkT0OSXFUkb/VesGAI7NuDJvvIj+yS/T/twg/PnzT4kvnXNaCUnFNlFRmCXSbbX2aKeQ1o2fSpeplMUZZbnB1op+crEMo96wKymrLRS9bK3r25+lAHzjgzS9vfEj9sLa02Llq/w5VKPDqCOEecjUWdfF77YIIkmyQoBmRV8r0r54qazM1g6IPrs6RL/jVWRbTT1uWQ1TjKpJHZklwTAFBd3WTSuidyv64T0oGnU0OBzplclICHnEwVxZN49+5a9w7r9d0f3f5OXRA46AbLFqYjyQA5JjlqgqGaZSzOe6PoSir4Cu2bOtFH1yjNKHxWpanm2cvSA9ejkDUIOxjZpdwyDFS8zl0VuKvmB/bpWwwYleKvo8UCuiGohhMBFGMMBIle+6nvpe/MTv0dR1/Apa5s0NoVzOlkjVqNaNm+i3DxMxOAaWF9E3GsCTtwGXvI4UuWpFFOecCjM5LhZ0duXwy5vSShD98B6a7vt1VVSzTSzrpnV6pYxVFKBaNz6KHvC/IebOAf/4ZuDckxSQTI0pa9mKY5Y/j2KULvp5LoneZXcZebqgpSXGmE3GbqKXf2s4jrJh4guPnLJzsTNb6D2/Nghmney+VSF6EXM5/bBQ9AxToL97JudKIfSqju2E6DknRT+yBzkxe3Vm3ZCiZ4aT6IOzh7GVT6NSLqIrGAXwYAQv/Yt/xcGTC55FU5RHnwOSo5Z1UzLUoinXjFco+goiGE1FcSZbbs6llwgEqB24VPTt+txIJEfp+FbzlGygKnrAzrypLJGwkWM+NgiAKR69VvTLg2rdGCVUWAyD8TBS0RCR9fW/C/yH2ylgChDhzhxuLkypLIGzAM6W6XCp1k2xKqwbMfgysRDi4SDOq9NocQEZ9Yad63z6IVIQV4tMn00vtddlLbmJ3pWzKyGX91sG0R+dzlOJeDgG/O5TlLftBavcXMwiwJzk6AGZfWQr+jy1QAiEnQubSHgR/fRh4NNvIIV50+20ADlA09/0FvvmmDuHXIiO0wIXF5w7ruFVZi69bHd6pfhOhJO45/B5/Mk3nsGTZwRpZjbT1N4vxdIqluqi1UGnSAxTTObUw6ToM1uRq9O4bMoV92pV3EqtyuNfWqDH4T008wWc1k0khBwSCNYKzo9XSaEWlrpYohMAjCLMUAL5Sh2Hp3J2ZpMkXgBlo4ZBkKK3rJtq3TkuVdTLaAQiaCCAfZvSMOoNzBVbrNQ1sN326Nt1rpSQ41/0BGomekHglSXndwVDdG4s68Y1o1kFbFiiv/nzB/GFx8WFLjz6MqIYiIeRjoUcZG1h4kq6SN39sKs58EgaDU62TyvrhjGGzQMxZ659bAC8soTX/M/78cVHxKA49E9kFe17Cz3fdCXlcxtF0edG9eh9VMvMYRpsHWS/eOHMYglv+psH8LXHzlj76UnAgG3dFKa9V6HywLHpPLYOxhFKKH1U6lX/RlFuoi9nSclzE3j/XcBlb3Ruv+lKRdGfw0KQLjzLunETfbXQnNnQUtGfB8JxnF6gC3Za3rwzW63f9ITVJG0ViB6gXi2TPxT2ym6r4ddMruosmooNeij6FmpVVgSrGTciw8xh3URDKPA4gmbVEQdJ1Ijgy/neiB4AZvIV2rf4MCVKCLDyIoJoCKJXFb2sM3FZN7UKGkEaZ/s20d/aOiC7vXtFL8eMXLFKZr5JwraI3iMnPzHSnF6pib47NBoc339+Fg+eFhemyLop8Sgy8TDSsbClVByQytht31RyMCP2SXdYN4ZU9PaFMJGJNSv6ag7TuZJdNfv8ncAVb7U94IkrSSWeeoge3R494EH0z9I+d9JvxgMHTy6iwUEqqh2iGZp+FmeIQDu4uRybKeDS8RRiSUn0+eb1YlW4iX72eXr+1o/RakJuTFwJzD1P31uawxyGEQkGWlg3BTsQK2ERvYeib9SBSAKTguhnZFZLhhbiQO6s998hL+A2M56esfNVREZTjwNDOy2iL9dMFKqKgGmp6FtYN/OC6If3WKIm40qvzEs7TvHpk3WayZTzXazzAABGEbUgfZ+VOTS0y1bKAMIVu5W0IxgbG6BU3SbrpgQzSHGgKzbTeGhJ9AM7SMTUKu1Xl5KQ51emZboVveGj6AG6kbkVvc6j7w5zxSoMs4HzBZNsgiq13y3yiKXoc56KXhD9zGHn69UcaiEiiGgo4LJupKK31a2Xome8gSQqWCgZdPGVF52Wi/ScX7ifHpMeRK9OTxsNm+h7xGOn6II8NtNBcyrG7KKpwkzbjBuzwXF8poC94ykMpxOosJidXukOxEq4iV4Gx4Yv8d5+05VExi8+AAA4zwexbSiOJSTRQNDDuik0q6Z4C+sGAMKJZkWflkTvE5BddUUvfHreAIZ22+0A4JFi6Zt142PdGHlag5YFgaGdlqhRrZtEJOTZqjhjZsVPdK/ojQB9n3UzHdrpUPRRQ3ynkl5ZMkznuFRRr6AeIEFxuaXoW2TeDIrMm9zZ9qtLScgxs+AmeqnoFY/eS9Frj355mMrSBTmdq9oLhBsl5E0i+oyfdRMfEr5vs6KvBInodwwnXIqevicVVRT9QAzTuYod/BEnOYMSFouGPUVUu0oO7iQr54QH0ac8FP3SaSKucTvLodHgNPXtEJLoHZW8rZAaE9aNx7qyLpxZLKFab2DvRAojqQgFZLu1bqSik+lvbkyIm6Po/3LWHMJYOopIKIRSeKBZ5amrS0nEhyhfXM2pdxG9VILTMtiZGKagbluiXyVFP7gDyFA1LIZ3o1wzMZqiYzrjLpryDMYyb1KR5H/uEP1GMGxdJxlH1o2HojeKiKMi/ttlfr1RRFUQvXWjGtxJN/oGzVYShlT0Y4iHlTx6wDv9uFZGjVHGzUQmhqFEuI2iF2Mse1rx6NspenGNyg6irTx6941VLbYyCqIwr7UVuhxsUKKnEzpXqIKLDpa8VkTOInof6wYgVT/tVvRLKAfowtg5knRUxko1lVCIfvNADPUGt4M/YsBkWAkLRcNWKirRBwKkUGUmjWolRDPN01MZiFVaHNx9+Dyu/6v76TfaoFCt47nzOQwmwpgrGB19xiqa8lqFyoVjwqLaO5HGaCqKJR630yvbKXrpM2dP02zGr5BkZA+lZIqWr5O1AQwmwsjEwygEB5uD11UPRX/t+4C3fNRZWKZc4GYojqklSfTiJsoYBQx9iV5cwKsRjJXY+Up6HNqNYtXErhGX9QEQ0ddKznoC2Yvey+6TN7izj1sFdNK6UStjQ8EAKkFxHCUpKmPTLHVv3VTcRD+0i1IURRwkURPfmRxDIMAQDwdRFiKLGu650ysrMFhU7HsI24YSrVMsB5VceqnoI+08enF+/awbv2Cs/KxF9NTQ7Es/PIU/+OcnW/9mj9jQRF9vcJghWiCcG0UUeZtgLEDEOfe8M5WxkkNRKJhdIwnkKzUr6FUQWTeJsH03nsgQkU0vSaKXir5IhCotCfeqQDI3HHAGYxlrro6VN6OxfdZLsiJXesqt8ORkFg0O/OzLSRl2pOqTY7TvRqGtdSPrCC4dT2EkGUGuEUejkiPScbc/kIgNkBUjL5Ds6da99ANBmtGIGdKJagoD8TAG4mFk2YB31o17Oj5+BXDgV52vRdKw6i3MiHXfcZBoZmvrYGwovqqVjrjsTRT8G7kUZaOOXaNEvDPtqmNb9XBRF0YfFkRfqSEWDiASclKFGZYdLAUpKse60a6zpxtGARVRsT5XqFKzMtnmW4iilLmIBgJWTCUZDSqK3qNyvFZGlUUQDjJEQwFsG4q3tm4yW2lml50kQRJJt68qD4bpmMlsHWsNCCUYy7lYy8KD6OtlInmjCISTeOJ0Fg8e86mAXyY2JNGfzdp37logbqdXQnr0YVp8xKu/9/hLqCBFBqQAoJpDjscRDwcxnomiwe1pY6laRyISRCBgK6TNAzRozwklaBE9K2GxZIBnT5GF4M6z3qQQvdvfdU9PZ56lAJKiPqUqdwSCffDYqUUwBrz7QBdEnxq3p+ptrJtj03lsysSQiYUxkooiz+Ool5baK3rAtm/aET1gxzYCYZwsxzGYoHNMHSzd1k2+Mx80ELCOq+yVctlEylb0AAVkfYOxi6tn20hc+S7gD15AI5JGqWZi80AMkWDAad1Zq6ipRN9iQQ2VjISiz1fqjkCsRD2csr8PQD1n57zzbpudGUWrNUmDi3E8KIn+JADy/8shOyssHglSrxvA7nejXs/1CqqIIBUNgTEmiL7sfc0DRNrpzbQsY4tFhpqQGKGsMMA7GFuvEJ94efQAiYIapf3mXT2FVhIbkujPZe3BXgnEAaMAViujpKRXmg1uZSs44A7IikVHso04hpMRa9DLKa3a0Exik1T0khgURV8zOeoLp4jA3NNnSVrxoeYOku7mTTNHHP48APL/4SqF98FjpxZx+UQal0+kEQ8HOwvIyqIpuT8tcGymgL0TRAYjqQjykIq+jUcPENE3GqTUfYj+6TNLWCrVrGPG05tQrXMrBjPHPRqbeXn0fhD7Mm/Qub125zAWSzVU62LMpDdTnr0XcZTmgcQqFEupYAwIhlGpm+CcUh7H0lHvfjdqQLalolfITQTAqaFZqGlTLpcTFNaNsWTXFAS67VMvZtsSM/mK8MwZkD0Fs8ExyJdQjtjiJxkJKT3pxbrK6g2mVkGZRyzi3DaUQLXewFyhhUU5sI3GXLvVpVRIizUUA0JipqoGY632Bx4ePUBjxSgCkUTTAi8riQ1J9FNLZVwiprIlxIHyIhivo8RjlqIH4G3fjF5GGQfSGhGLjiyaMQwnI02fLVZNRyAWAEZSUYQCDOdk5o244MbC9JwvnvImsPH9NH1MeNgima3A7HPAuafI/pg72tTjZqHUmaJvNDgeP72Ia3YOIRBguHQ81bl1I9FidamGlXFDF8toKkKVlFZ6pY+ilxkwlSVReWp4HqdKzcS7bnkIn/rBCcvuqiWpbF5aN9Om6Bgq/em6Qd/XJdHPVIKIBAN42TbxXAZkM1up+MtrhbLS/OorevlTQqwkIsFmovdaF7kT6wawib7s7HMjwV1r5xo58sjP8yEEjS6IvtEAakXkuT0mZvNVIs2BbcDiKdHQLIdK1D6m8UiQVpkCvNc0rpdR5mHr2tw6SCpbne03YWA7zSLbrS6lQp5n9ZiqHr27c6X7c6V5y6PPV71vqiuB5S4OfpIx9jRj7BBj7KDH+4wx9neMseOMsacYY9cs5/c6xVS2jKu304Et8Ki1Ak0ZMo+eDqZnQDYUpQZfMvPGalEcw1AyYvXllpk3RWHdqAgGGMbTUYtwc6Kp156MCYAjuDRpT01VhOPUF95LLb/m9yi494V3AM99m7xsV2qltG6m3aXwLhyfLSBfqePaHUQEe8dTVvC0JTpU9GezZZRrpq3ok1EUECcC6FTRW5lJzcdJxiJOLZSsYHQ1Rhf7QJyCsefrgtCtXGVXQ7N2EBfm+TL5u5uEHddRLn1pYXUDsepPyRhRxEPRW/t4xn6tE6IPhKzjnqvUPK2beCyGKqIWkTXyM8jxOGb5AMI11+zw/DPAF9/l3QJcxGPyjShGkqSIZ9TMm8WT1NAMS6jFbKJ3KnqPdZVrFZQaYaTEtb5tmK7BtimWubNkdXWq6JMeRC+FjIPofayb8qLl0ecr9XVt3byec3415/yAx3tvBrBX/LsZwCdX4PdaolIzMVcwsHs0iYF4GDkzYgWjVOsGgHcuPUDkMX3YDqQAmK5FMZwIN1k3chlBNzYpufRnl2oo8Bh2JmrIoIhgLe/vPb/lI8AbPtz8+uAO4L3forqAr/0aveZL9K0V/cGTpPCu3UlEf+lECudzlaYFVZqgZtq08OilDXSZat3wBEJm2e7x7gV1VSQZsPZIrZSzj3PZMk2Jr/gZzE28CgAwmCBFP2WI6bO0b9ytcNtBXJjnisC24QQm0oLorerYFrn0F1LR10TWVySI8XTUGYxNTTibplTwlQAAIABJREFUdQGtiV6q2MGdVlsQIp/m8Z2IBFFkcYu8eXEW8zyDPBKI1F2EfvIHlBl1/HvNvymKhXJmxAooz6q59NlT1OeG5VCPORW9Zb1abRCUzJt6GcVG2EoLlYq+bYplo07xuY49enGTUY9pIGAvEO7XTkEKAYdH71zgZSWx2tbN2wF8nhMeATDIGNu8mj8o7ZItg3GMp6PI1u0MjxKPYiChWjc+xLbrespTf+IL1omarkYwlIwoswG6wOQygm5sGrCrY6eyZeSQwKaoge1MEI8f0e/+Cbsgxo2RPcB7vy3aD4Ts1sICix0S/WOnFjGSjGCnSMmTFssL7ewbqeKjGeqN4wM5O7h0jL43FQ2hHBAZKMU528t0Q1X0Mod+sJnoZUaPzK7CL3wRx7e/mzaPUzB2tiHbIAiV5+5c6YHJhRL+4cEXKWAn9uV0gWHHcBzjGZqFTLuJ3t2X3qzT/l8goi9WbetmPB3DQtGw10IIBOj4KYVHVnqlFwJByjaRaxOABI2XdZOMhOzaCACB4hzmMIBqMI2Y6WpqJjNijt3b/JviBrxkRjCaouvLkWKZP4dybh5pVgZXLM2kSvRe1k2tgrxpWzfpWBiDiXAbRS+uSaPN6lIqLOvGdUxlq2Krc6Wb6JXGZkYJPJxY18FYDuBexthjjLGbPd7fCmBSeX5GvNYExtjNjLGDjLGDs7M+vdc7gLz4twzSxblQsw9chUWRioSsu7xviuU176Wuknf+PnUJBDBj0NRSDnrVuvFU9Jk4zi/Rgs1ns2XkeBKjoQq2MfG3tcsm8cPYZcD77wF+/vMOC6RSM61MIPeatm5If56JYPCl40R+x9oRfTRNaYNtUiufPZfDpkwMAwk6Vowx+0KoFf0VvbrsYvY0qSUPBX5czBjO5yqom0RqWRGfkHUSC5CNzWS/cGUBah9844mz+It/OUJtKsS+zFaC2D6UwHAiglCAYVqSUGqCYjluRV/JAuCrVxXrQsmwey2NpWk8zKvNuwZ32Iq+0WhPYle81Wocxzn3tW4S0SByPGEJoWBlDvN8AI1oBomGaxyVFKJ3L8ItbsALtUhzQFnYR4Fzj9P+KHZhPBKyiT4+DIDZRM85UC8jb4Ys6wYgVd9R0RTQuUcvr4Umok86g7Hu9wNBiqGIYKwZTqBm8vXp0QN4Nef8GpBF8yHG2Gtc73s1YfHMb+Kc38o5P8A5PzA21jqjoxVksGXrYBzj6RhmDVs9skgSgQBrHYwF6CT87KfoLvy9P6NtkXAoejXrxh2MBSjFsmSYyFfrOJsto4AEEo0CdgXbKPpOMLLHboYmsCiIbttQHLlK3U49c2G+UMWLc0Uc2GlnhWwfiiMSCrQPyDJGQdg2qZWHJrO4artTwQTiykD3I/pQhKa8lWzL1Eo5Y2hw289dEudjQBRMNfW76cC6kX1i7n9+xlJgZUSxfTiBgIi7WMHYQJDIPufKpV/tqlgXJNnFw2TdAK52xSrRW3GKFkT/zlusDqbVesOXfJLREJYacUvRRyrzdMxjA0hxl6KXx6Q4C5x7wvmeIPrFOqnvsVTUad0AiJynzwSUBABS9OL6DYboeEvrxqwBvIF8PeRQyNuG4p0VTQHtO1dKeAVjAXs5QT+PXn5WWDdVRtfEulT0nPMp8TgD4A4Ar3BtcgaAOvfeBqDNig3Lw1S2DMaAiYEoxtNRTFftQRqM0kXeMhgrkRoHfu4z1E8EQJ6TqouGgo5+N8WqaTVZUjEhgnfnlyo4u1iGEU6DVZdwSXgelUCyo17lT5xexPee9WmF64L052UDJz/7Rqr2/Vts4g0FA7hkNIlj0x7BMjcuvYFmOz7IlgycnC/hqu3OLINgXBnofsFYwK6O9SH6mtnAi3NF7Bd/p5zBLZVrYAxIR0PIxENYQhKcBZutmxbB2Lwk+udUoo9gh1hjYDwTc+ape+XSW1Wxq5xeKX9OUfTSXnIEZAd30s2uWui8h4uA1dDMx7pZ4nFaZaphImosYg4DCCcGkWQV1GtKGmNxXsSTGHDUZd+I8zJfC4u/QTnGoqAwOSOIPm0nAyRE1o3VZkRNPxYB3mIj7BBhe8Zomc8pv8ybSNL2zrtNr/Qj+mqO4moyE8fx2WGaEfIGKozez6w3Rc8YSzLG0vL/AN4IwL0O37cA/IrIvrkOwBLn3KeccGUwlS1jLBVFNETpZjnTVvShGF3kiUjQXnykFXZdD7zhT2HERrCANIZFVkAmHkZOVMeSom/26DerRJ8tw4wQge0MzGE2NNFRx8mP3ncUH/jcQbuNcAt0SvSyalaSl8TeiXR76wYA3vpR4PV/7Pv2oUnyJK92EX0kqTz3U/QAEWx5kaoNPYj+1HwJ9QbHay8ndXdWIfqBeBiBAMNAPAyOAIzIkO0PG+2tm4IYDwdPLaI4dhUWkpdimg9h+5Ag+nTUVTTl0QZhhRW92eC+szOgOb0ScFXHymOolvZ3SvSyoZlPMJaWE6T+9Qwc8zyDcJJucMWcktJZmqN40rYfA47d4/wisTB23ow2K/rUBBCKIbNAbQEiAwrRR0NU4iLrGoZ20tKcgLWKWQURx2zkP/z4DjAGfPL7SjGkG1LVd5xe6XNjiCjWTWzA+3pPjFhVtWXY7RpWA8tR9BMAHmSMPQngUQDf4ZzfzRj7IGPsg2KbOwGcAHAcwKcA/Nay9rYDnFuqYIuIsE9kYlbFHQCE43SRM8aQioZaK3qJV/8O7n7TA6giYhN9LIRcuY5yjYpVEp4evU30U9kyWJyIfgtmcY511j9eTjN//5+fxF1Pt74/SqLfv5kGnF8u/eRCCQEG6xhJ7B1P4Wy2bE+He8ShySwYA162zUn0sbRK9G0U/fxxylH3IHrpz79mLxG9DL5nS0T0AKzHcpdEn6/UEAkGYDY4vl+6BH9z2ecQjqWsWMNEJuYkUa82CCvcovgT9x/HG//m33wrOktKMNZubKace5memj3d+RJ5AtbqUl6KPhpCASLrRvSZKYaGEEyItGZ18RG5NOZlb6SF3fPKLFUo+iKI6MczURQNk1InGQMGdyBqkHiIDNjXjaODJUDFg3PHyLYRXSMr3En024YS+Llrt+P2H03aVetuSJ++U0WfGifF7u79FBYV+V59bqw/YtgaP0W+Tq0bzvkJzvlV4t9LOOd/KV6/hXN+i/g/55x/iHO+h3P+Us55U679SuNstmylUo2noygqRB+J2yfPt1WxBxZLtN2QIPp0jBR9wbXoiAo5jZ5cLGEmX0U4OQhUljBunsdko3UwE6BA2NRSGTe9YgdevmMI//HLT5Cl4LePLkXftKycwORiGZsH4ggHnad+73gKnAMnZrtcBs6FQ5NZXDaebopbJNNKcDLYhujnjtH/PXLopT9/1fYBZGIhh3UzKAhJEtNs4jLqBpqdVIKxrT36q7cPYiAexv3Pz2ByoWQtDQkAE5kosqUaKrJQJ7OFpuZqfvgKtyg+Op3H5EIZJ+a8z4vsnpqIhBAOBjCcjLisG3GzlIVAQPfWjY9Hn0ccgVrRWvKvFBlGSMzcynlxHBomzdCSo8Den6bXjt9nf5Eg+hKPUTA25bKfhH1T4lHEk/YNyl5lShL9fmqCNv+CQ9Gnok7i/NDr96DBub+ql0TfaXplJAn8+veaV2YLJxRF7/NdiRHIkKWsDF6Pin7dgXOOqWzZsk3GMzHrTgkA0YRK9C06WLqwUDRIXChEkqvUrUHmlV4ZDQUxkozgidNZcA7EUkMAbyDGy3ix3l7tEaE0cOl4Cp9534/hktEUPvytw77bLxQNBBiplqalDBWcXig12TaAmnnTgU/vA845npzMNtk2AJAeUDzrdope9g7xUPTHZgrYOhhHIhLClsG4RfTZcg0DCboRpyIhBBhw/1aRCHb3HwlCYS0Xd8hX6hhIhPHay8bw/ednccp1rMbFLM0iIasvvaLqS/NkTa3QIhJzYnlA2VLajbJhIhoKICh6LTXn0o/T/iye7MG6aV4vViIprBsAlmVSjYwgKqybiuxJX1oAZSGNUruK9BbgqGLfiJlWETGkorb9ZC2LKG728zzjaBxoKXpRR2C1A5k5Yit6hJsEx7ahBN59YDu+/KiPqh/sUtEDtCiOW0BYefQtFL1SVJdv2CJyNbChiH5RkOMWH0WfSNrT9m4U/ULRwEA8jJBQwZlYCPlya0UPUC699KyTGZvcj1aHrbRAP0jvecsAtW143b4xK13Tcx9LBgYTEQQDDJtEP3wv+BH9zhEaqKfm23e+9MPphRIWS7WmQCwADA4Mos7FcGvn0Vsf8i6WkhW3RPSi+rhsWzcys+oMHwVe+wfAc/8CPPttuhhbdCQsVOtIR0N4/b4xzBWqODFbdCh6mdXSlEuvBmRLi3QB97jqlxvzojfL4z5E7+611FQdK+wPp3Wz/GBsQqwbC8Ai+lp8BLE0Eb3Vk15mPSVHaF/23kCL68jWFEYRHMxS33ImbM1IRebNAstY1x/9vuhJLxX9yF5qHzLzrK9HL/FbryNVf4uXqt9xHfUxGtrd7vC0RkQpmPK1bmxOkLFEreg7gJpDDwgCFg2GKjyMTMImGN/FRzywUDIwnLCDumTd1C1/0Cu9EiCfXt4MMkO2XXOGjyFbbj2bUAu/AGAsFYVhNnxvTovFGoaEl9wUNBQoGyZm81VsH27OAIiEAsjEQpYF1Av8ArEAMJqmNggA2it6gC4Cl0oyGxwvzNKqVQCwZTBm9YrPlgxrxgWQT58r14DrPkStnGefbdvnplCtIxUL4TV7xyyedlo3rjYI8kYk+5EDK14VKxX9QR+iLxmmtRAH4EH0QM9EL68P7/RKRdHPH0cdQSA2iLgQNHXZk96yssT433sD5fKffYyeW+vFMiSjQcW6cWbeLDLnmJI3NyumFI5Ra2VV0XNvot8+nMC7rtmG2340adtwEluvBf7Tc8u33hzWTXuiXzIjlL3skcG3EthQRK/m0Esk0nSQZfsDiW6sm8WiYQViASATDyFXqdmLjnhYNwCs/igAMDjsJPp2hCqnlZsH6TtGUvT78wVv732+WMVIMmr9rle/G1kVuN1D0QMUg1gsdXZMvPDE6Szi4aDV+kDFaCpqE0Mnit7DtrFWrRKVvJsH4siW6DwsKYoeoHO0VK5Rbv5bPgKAAn5+4JyjUKGaiJFUFFeJYPL2IXssTbi7kma2kZCYPWp/UWl+xfz5mtnAYqmGZCSI4zMFqyhMRalqOpaxHE/HMJt3LRI+SK0Eesm6CYlFPtxIREL2jXvhBLJsAOl4GMkM/e2NkqgIlcFwWVi0TWRgW0RfQF0sDJ6KhjAkCtPc1k0u4CR6uU+ODrTjVzQpej8R9so9IzDqDWuZyBWHFYzN+ge/FaLP1sNkOQZWZiboxoYielvR20SSTNuFLw4i6EbRFw0rEEufDcOoN6xMl1aKHiCFHRHeZT2cRg7Jtis6nc2WEQ4yjAryliQ+7/O5xWINQ0k7O+R8rtnmmWxH9ImIVXjVCw5NZvHSbQOOKbb63bai92mBALQkequ1griRyBv60ek8Gpz63EgMxMNWERV2XY9/SbwT31zag3876l11Xa03UG9wq5Ly9ZdTFoW0tOhvCCMcZPZNNBAAxi6n2YJEeaGJ6H/98wdxxxPtU2TdkGPkdftoX5443bygR9GoO+o4xtM081tSZ4yDOyggmjtLN6YOl6yT7Q+Yhw2Viirrxi6exAIGkIqGkMrQOLd60kvrRpJaeoICnlNU7UrrxdJ4TEaJ6EZT0SbrphBy1iU0KXqAArILJ6xunRVEHJWxKi4Zo/PaS/JBo8HxsfuOtm41Ek4A4HTTcXeulFCIfqEW8bTIVgobjuijoYBDfY8OpGHwIMq8WdH7Lj7iwkLRad3ILAQZ8PRKrwRsRb9lMG4RWD2z3frOVjiXrWDzQNy6w7dT9Asle9YxkYnBqLsudgCn571z6CWGk5HOlhT0QLVu4shUDi/3sG0AsoZKYjnGXhW9umoVYNtaR85RNknGbd2IG7nZ4Pj9/Hvwn2u/hv942xPWcVBh2RTiXP7q9bvwkXdfZS3RB1Ba7ng65uz3P7YPmH3efu6ybuYKVdx3ZBrfPeKfMeUHadv85OXjCAaYZ0C2bJiOGWXLXPqZI10FGf0amgHUAsFaN7ZRx0wjg1Q0jGCIbgBMZvgUPeoKtl7jsG6MIJ1HScpj6ait6GMD+PTI7+H+5Judv+9OrwREQJYD5yjv3gzGEA1539R2iwZqJ+Y6XC9ZwbGZAv72e8fwjSd8Fp4BnMH4VumVAgtGaNX8eWCjEf1SBVsH4w4FQgHZuId102LxEQWccyyWDAynVOuGvkfOIPx8NUn0W4fi9l1dTEUX2ijnc0t29hAAy7v0WjiBc47FooGhhCR62tadeXN6oYxEJGi1g3VjKBFBtkfr5rlzeRhmwzMQK2GIBdY7I/rm1MrjMwVMZKJWFog8Ps8Koh90zNhsRX9itoByzcTv/NRecM7xG198rKkIScZSJNlkYmG869ptTWp2POPKahm7nHKhy1mRSph1ZFMcmaJ980uPbAV5rneMJLB/c8aT6IuGszJ73yYi8ruetld8so7ldHdE79fnBqClMy0rDsBMI20duyJLIiiJvjRH51RdSGfrtZQFVKQ+L7IqNCn+DkerCQD3Rt6AXMzZIssiehGMzZYMFAcvozenqJI21GIpx3QsjLF0FC/2oOjPZkkotGwZEvEn+j/91mH8jzufFYVURMFzmug7x1S23FQIJDNv3NZN2343AvlqHTWTu4KxQtEvSUXvrRokEW0bjFu5tKERUlftPPqpbMXxt0jraM5D0ecqddQb3FL09gpXzm0nF0vYPpTwnIoDZE30quhbBWIlrOXnWgVjB3YAYI5FzyWOz+Qtfx6gGyljwLPnyHseVM6Rat08dYZshLe+bDP+7qaX47nzOXz4W84iblkV6867dmMi7cpokml9s8+LJfu4Q73K2caLcwW7XL9DzIkbymgqimt3DuHQZLYpW6tsONdD2DuRxg37J/DpB0/YMzq5/mqN1sytmw386j8+iu881boIj6wbb/IJBQOohuxYzBwfsGZD5UAKoZpU9HPNC+lsvZYepx6n9WJZDPFw0EoRdSh6NM9aACWP3jBxer6EN37sAfx/383TesTnngIAhGOtU1x3jyZ7ugHLxmjHZ1sQvY+in1wo4fMPn8T3n5+1G5sFwsgaq5daCWxIoneqxfFMFCVO/zIuRQ+06XcDm5DdHj1AmTGRUKCp+Ehi62ACQ4kwVYkGw8BbP4bQK34dqWgIC0X7d932kdngOJ+rOP6WcDCAwUTYSrfz2kfVugGAaVcXS3cBkBtDyQjKNbM5E6EDPHkmi7F01DELcaMRESTdStGPXkpZDztf5XiZc45jMwXLtgHomEykY3hOkKkzGEtxlErNxNNnl5CIBHHJWAqvu3wc73z5Vtx3xNlDSI4Dv3iLxETGldE0djk9zj7n2f7gsFD0lVqjo7V8Vcib+mgqgmt2DqFcM/HceWedQ9FwBmMB4HffsBf5Sh3/8OCL9v5I4ommMZWt4P7nZ/H/3n4Ij5yY9/ztSs3EyfmSVW3rhXAkDhP023M8Yyn6cjCFqOxJX5pr7na6+WpSsmcfs9aLVVNEx9NRzMtFwiHiEK7zEgzQot+nFor45c/8EDP5Kp6eKtAKcSYdt3ZEv2csiRd7IHpZsX58puBv/aq9bZSCqc8/fBINrtQJiOyy1VxGENhARN9ocIymog4iACgL4W/r78JnzBstxQHYRN8ul37BIlEniQBkjbQihngkiMf/5Abc+FJa5g4H3g+M7sVQMowFpZXsL9z6CP7s23Yx1GyeBvnmAefsZCQZcbaglftYct6Mmnqng4jSL4deQt4oegnIHpsuYN+mtO9sAQBYjIi+1Ajhqwcn8Yf//BRufeAFPHR8zhlPSG9q+uxMvoqSYWLPmDPlcvNgzGrPrAZjrXbS/7e9M4+Oq7rz/PdX+6aqklQlyVpsLTbe5BUbjM1mVkMgIYQkhAMh6SSEbjIkk5xs3ckkMyeTM52eyaFDdxbCkoROdzINPsDQJA1ZWAwJARtjbIN3Y0mWta+lWlV3/rjvvnrv1atSlaqE7NL9nONjqeqp6t1afu93v78tmsRbPWNY3ehXPcbmoBtj0aTOwxYNzWb6stX5XRiPpTIXw8Bi3rp54JCm/UEmcHjw9Jh6XsUalaFIAg6bBT6nTe02apRv+HB6/Tmvbgzgus4GPLzrBM/UEbn0AOCsUrOvHDYLPvvoblMJYueeHgxHEvjIpuxaBoHHZUNUkeOGWED9LiTsVXCKnvRTw9kevdMHhJarhj7CXLrXPVzlRJpl2i1HE9O6Yin1+R1W7NzTg4GJOHasbsDpsRiStSsAAClY4XHl2TkCaA/5MBxJmGYz5UN49BOxVHYqq8DEo4/EU/jVa12wEDJzA9w1msHg0tDPiMVC+I97L8Fdl3bobq+rcuI/0luwz7lRl7o04/ARBWH0aryZD414Q4YjiZyplQIiyjJ+NV4nhhUtfGAijr+cGMYLhzLZID0m2UMA38KbafTDym1Ce3farKj22HUe5HAkganEtGkOvUDk4Rcr34gZscaLrBGL0sFyxz+9ii8/tg/P7O/Fd595B7c9+Co2f+d3eQ2haMbWbLhQaeWtgCEYC/CL4MHT41jTlJGUAh4H0gy6iVqTeXLGtWS1ArZY+IwAE49+KpHC8cEIdqzmF65iZYLBiTjCPieICI1BNxYFXDpDzxjDVDJb1gCAz1+1DJPxFB58SfHqhU7v9KuG6id3nA+7lfDJn/1FJwlOpxl++tJxrGkKYGtH7poAr8OGKPH3YwgZjz5pr4Jb9KSPDGbG7WlpOl8x9JOIMKduV6JWxypGdMpEugFE2wfCT+44Hzdt4Br+oIfPuk3AMaMUkgnIFve+dI9G4bJz05mzEaCJoX98TzcmYil8aGMzAOVC5qsDcwXmdIwgUEGGPheibD1gSF2acfiIgpBK9Fk3mceaaatvRo3Hrsotrxzj6WfHByMYU4y/yKE3xhtCPqdp1o3q0WvOsd6vz6U/laNrpRbx9yOR4gKyp8f4jNiZDH1V22YcwWJcsmYpHrv7Iuz71jXY/Y2r8N0PrkFiOp03uKWmhlbrz1+kWDptFrg0Xp94f984NYpochprmjPbZ3FB09YMqMHYGaUbpVmdTr5ZqRh60aKYB2Pf7p0AY8D2FXVw261FB/4GJuMIaZIANi6p1hn6WDLNm+qZJAOsaPDjfWsX4ZGXT/ALt8GjtxCwubUGD965GQMTcXzmF6+ru5RnD5zBicEI7r6sI+8OjY8TVMb/aTT6tL2K96RnTMlCMunt1LSR3xcbw0TaqQZiASCsjG3sGYlif88YphIpuE3W+KVrzsPDn9iMS5aF1c/eCeLrjMGh28GbMdsUy56RKWzt4GvK+Zk1BGPTaYZHXj6JdS1BXKNc+Acm4sBV30b8ffcjlZ67oSPAAjD0fpcNTpsly9AXGox97eQwHDaL6mUAmTbH4udiqdakMb5yNKORvtnNA5q9Sll/lnTjc5h69EaNHhCGPmOMuhQvLp9GP1vpRnzYtYFSM9ZdfjOWffst/M8Pb8Km1hoQEWp9Tlyh5Inny0vuGubn31ytf01ETMD4/orfdx3hF1KtRy8uaNotuzHrJhfi9Ts5pDEO4eU8R33kJP9d8ehFILazKYDWkBcnikzlG5xM6DTy1Y1+PsQmLsZY5i/Yu+uSdkQS03jx8IDB0EfR4HfBYbNgfUsQ9310Pd44NYovP7YP6TTDj184hiW1HuzozJbQtPDGZopHr9Ho064AfGyKF02lk+aVwiIgC2Bs2qkzcmLXdNeju3HD/buQnGamDsrNG5txidLFdEmtB3YrYV+Se/b5cugFLTUe2CyE45qg6nSa92vKhZhJvaEliCqnLbehFx49WQCHDy8cHsCJwQj+alurfsdS24Hx6k4AcxuMnbtLyFkCEaHO7zQx9DMHY98diuDxPT24Y8sSuDVfJiJClcuG0alkzj43+ajVGPqXjw1ia0ct/nR8CHu7RnHpeWH0jEbhdVizugbWep0YiyaRSKXhsGWu0cORBJw2i+4L3+B3qWmHQEb6MHrEWkTWymwN/UwefS5CPl7+3Z9L7wQ//7oqp85rBzK7Hq0+D2Q0+leODcLrsKI9lNH2xbHaVNKJWAoOqyVn3rWgpdoNu5X0MlOY68I49SfemVNp3SD0+caAC+0hLw6cHsv72EYGJ+NY25TJ2BAXtf7xGHxhn64XvRnLG6pAxNNL0aQ39M2az8GOzkX4yo7l+N5vDyGWnMab3WP4zk2dqjOTC68jUzQ1DL+6GyJXABZiiAwchRcwHz1Zv5q/VtNxjKftuu9RU9CNuy5th81CWNXox+rGgK6ewQy71YLWWi92j7rB7B5MJewzesh2qwWLazy69/Kx3V346uNv4V8/c6HqtWsRsmpzjRsddb48hl5xSJx+gAi/+NNJ1PuduH7NIlWSEp/3TPM46dGXxJ0XteLmjdl5uDMNH7n/D0dhsxD+5vKOrPuEfDMb6UZktxzum0D3SBTXrm7A0rBPTVHsHeNposZtc6jK3BAPKy0atMfXB1wYnIyr6XinhqYQrnLqLlhGgjk0+u8/dzhvKt6xgUnUeB26HUUx2KwW1Hqd+kIkA10j5hlDQroJuvXPLS7sI1NJrG4K6OIzZhe0yXhyRg9QnOviGo/OC0SdYui7X+dFMMr7cPD0OFYt8oOI0BbyomskmhncPQPpNMNwJKEWygE8tRPIpM1GEvmb6rnsVjRXu7kGLVIsFemm2RCr+evLOnDL+c147mAfQj4Hbjm/ecZz9DitGEu7Ebf5kIAdVUpqqsXNd0/xvqMA+EXg/zx7SM2iAcCz0BatAwCMpRy6NVgshL+9fiW+smMFbljbiLaQN6+EJFhW78PRwSmw0ApEWXaLYjPaQl6ddPO08jl/5OWTpseL+EY1/mCdAAAcFklEQVRztQdL63y5UyyFR+8KgDGG3e+O4KqV9bBbLep7Kgy+cDaldFMin76kHTdv1H9wjcNHDvdN4Ik3etRMjBODEezc043btyxRdX4tIr/YTB+dCaH3P/0mn060bWkt1rcE8WbXKBhj6B2LYVEwO2gq2iAYc+lHphI6fR7gaYDamao8hz53IBbgHo5f2aloeWTXCfz69a4cf8UzbpaGZ+fNa883v0cfNT1/4eUay8e1cZQ1TfqCFVONPmY++9eM9rBPr+sGl/CU0em4KlOkptN458wEVisjG9tCXkynmRprmInRaBLTSiaZoM7Qa0edF5vn4t0eUs61bhXQeQsSiy/BmfGYzqMH+Pfhux9cg49sasY3b1iVtXMyw+uw4dfsSuxqvZf/rgRUbR6lCnyQG/rnTqZw/x+O4p0z4/oHUOSb4VR2O+HZsDTsw7tDEQxd+GX8IHVzQRfu9rAXJ4YiSKcZxqaS+NOxIQTcdvzu7T51F6xFpFY2Bd1YWufDwEQ8qwIdgM7Q943HMR5LYblSzOa0WRH02DWGPnc76HKxIAx9LkSr4t6xKG776av4wq/34s5H/oK+8Rju//0ROGwW3H1ZtjcPQPVezMYIzoRIg3x6Xy/qqpzoCPuwriWIoUgC3SNRXixlko8uAnNGnX7I4PkB3LgRAV/f+RbiqekZUyu156b16MdjSUzEUziWY4vKGMPRgUl0zFK2ERhjClqS02n0jkVNPfoar8M0BuOwWdTGV2ub9Ybe77LDQtkafeGG3ot3h6YyHqrFykflAeqs2OODEcRTaaxu5M/dpgT+Cg3Iqjn0mtiQqLRWDb06DyH3ebcrueJpiwO45SGcsTUjzbJjHQB/zb53yzp8YH2TySNl43Fa8UpyGV6teT/cdqva48ju5cFoGuYZP7sH+e0nBw2GUzH0oylnWQx9R50PaQa8bt2A59KbCpJC2sM+JFJp9IxG8bu3+5BKM/z9h9bCSoSfv3Iy6/jukSnYLIR6v0t1bkzlG5sDsNgAVwCHlVnM2hhWnabLaKZLqDT0c0KVy47ByTg+++huRBMpfPHq8/D6yRFce9+LeGIv1+a1QVgtqkc/S40e4MZg29IQiEitKH31xDAGJ+NZGTcAUKt4d8bMG237A8Ha5iD+181r8MLhAXzuX99A71gsbyBWYGxsJto89IxGTWeXDkUSGJ1KzlqfF/DWyuYefe9oDGlmHl8g4tv8Wy/IzvcWxr/T4NGLubLadU7EUgV5gADQEfIhMZ1W89EBZHR6xaMXerwYwi5iBIXm0meqYjPvq89pg9dhVV+nmYKxADdk0eS0miUkztnM0BeLz2FDIpXGSCShe+1ET3r7KDf0L5/msosugA0Ay65GYs1t2JNeNqtYlxFhSIUEWsjFo03zvvxm/xksCrhwzap6XLdmEX79epfaoVbQMxrFoqALVgupn/mjuYb12D06Q6/t6hqucqojH8eldDO3VLlseOnIIPZ1j+G+Wzfg3iuX4el7L0ZLtQc+pw2fzeHN878tTaMXiDzl5Q1VcNkt+O1+3qPErMI0pDY2M9fojXx082J8+8ZVeO5gH6bTrCBDb2xsJgw9YN4AKpNxU6Kh97swFImbDmTpmsE43bm1FZtbs1sD+902+Jw2tNVmjw809vWZjKcKDoaZpuUZDX3POJw2i2rggx4ewyg0Z3tQeQ+Mlan1fhf6JvTSTT5D3xHSn6vQmPMF5QtFODl9E3FdKqNLGRvpmjiJtN2DHmXJWRc5dxADV34fY/DNamdspD3sBVGmy2ch303xXr7VM4YXjwzg2tUNsFgIn9jaiolYCjsNjcu6R6JoDvLXrqXGA4fNkjsg6w4C3hAO900g5HOojhrAe1eJ6tizWqMnohYi+iMRvU1EB4jo8ybHXE5EY0S0V/n330o73fIivthfvPo8XL2KT5jvCPvwxD3b8NJXr8hb/i30NLMxgjOhzcnftpRH9u1WCzobA3jxCC+cMvPofU4bHDYLBjXVsUllGInRoxd8YlsbvrqDGyHR8CofRgPYM5qRU46ZyA6lZtwI6v1OMGbetG2mPvq5aKn2YHNrtWmP74DHnpV1U4xGD/AgtIpq6LmRO9g7jhUNVbqWzW1FpFhq+9xoqfM71dYWMwVjtecqnlfk0DfkaVVRKOKz3z8e03n0XqVVsSsxjJid71T9LhtOmlzkIjNMaSsGl92KlmoP9vVwQ1+IFBL2cdnoX/78LhKpNK5TUko3Lg5iXXMAP3v5hK7NQc9IlDcpBG/D0B7y5jb0H/0lcNnXcLhvMiv1uM6fmRswEUuBKL8EVyqlePQpAF9ijK0EsAXAPUS0yuS4lxhj65V//6OE5ys7N6xtxKcvbsPnti/V3W5Vtvb5KEW6Cbi5RtwW8uoM+vqWoJqVYebRExFCXofOo1crd325M17++vIO7Pnm1bznzgwYG5udHo3CZiFYyFyLPNo/Ca/DmrfHTSHUVYnpTdk6fddwFFYLFf0c9926Hj/42AbT+4wSlZguVQg1XgeCHrveOxfNzTwhMMZwsHccqxr1khE39IVr9FYL6TpyAjxtVnj00QKCsfV+J7wOq3qR7s4xHH42qB69oRWIL6CZnEQBOKwWXLWyPlu6QUafLodGD/CdZSzJv0OFeMhEhPawF71jMYR8DmxSdoZEhDu3tuLYQAR/Ps4L4RKpNPomYrrBRsbMmyN9E5ke+YvWglU14Gj/pBqIFYR9TsSSaUzEU6qTMVdDR4ASDD1jrJcxtkf5eQLA2wAKi+KcJdy0oQnfuGHVrF7gUqQbi4XQGvKquwiBtsWvmUcPcJ1em3Ujqlhrcnj0gkJTH42NzURH0JYaj96DVTjazwOxhaS/5aNe7c+TrdN3jUxhUcBlOtAkH1Uue06vLmjw6HnWTeHBsPaQV59iWdMO3PQjYO1H0DUcxehUEp1N+slCbSEv+sbjWbqvGYOTcdR6HVmfTVHxzBhT56Wa9YEREBHawpkujd0aj7RUhNwyYqgncThdmGJKq+yUDysb/VhWX4XByURW3UqkwIrkQtHuLAt9TKHTX72qQVc7sKOzQZFTecpl71gUzBDIXlrnQ/dIFLHkNB7900lcc9+L+N5vM/MJTo/FMBlPqXOOBdqiqXztoMtFWTR6ImoFsAHAqyZ3X0REbxLRb4gou/ds5jHuIqLXiej1gQHzKUBnE0L2mU1lLAA89bmL8eVrl+tuEwHZGq8jZ3pbyKf36EXjp2pveT4oxurY06O8L35H2GeaeXO0v/TUSiDj0Ztl3nQNT5VFU9ai9ejjqWkkptNFaaRZKZZEwPrbAE8N3ujibQqMLZuLCcgOTSZ0mq6gTjNUZiqZgtNmmfEC2BbyqRel7pGpsgRiAX1qsbHdQETpgfNuzI31zQG0hZSKYkPmTTmlGwC67K9Cd2jtIf431xkqgT0OGy47L4zfHjiDdJplUisNhp4x4KuP78M3nzwAu8WC5w72qXLP4TMiEGvw6DWGPt+Al3JRsqEnIh+AxwF8gTFmSJTFHgBLGGPrANwP4Ilcj8MYe4AxtokxtikcDpd6WnOOyNuerSfic9qyts/N1W6EfI68EkWtod+N8Ohrvfk79RWKsbHZ6VG+VRUtXbVFLxOxJM6Mx0pOrQTyV8d2jUTzNmObDdUeO6YS04inpjW96At/L9tCXvRPxE0rq/d28dm5yw1fbjXFsgBDP2jocyPQzhqYips3+zLSHvKiZzSqvl/GHPrZotWUjUZ1ysI/E/3TVVi/OIhWcZEzyDeF9hgqFJEU4LLnbh9u5MZ1i/CJra24yKSB247OBvSNx7G3e9Q0kC12EE/uPY0PbWzGN29chZ7RqCpzqhk3Ro1eZ+jntnMlUKKhJyI7uJH/JWNsp/F+xtg4Y2xS+fkZAHYiMqmHPve4eGkI916xtCDdu1CICJ/c1pZV3KWl1ufAYCShegyZFsXl8ei1jc1S02mlL74bHWEf4qm0LgtH6L6lBmIBXnHKZ4XqPfpYchoDE/Gye/QBZZ1jU8lZGZuOPEZ7b9co1jRlz85trS3G0CfUqWJatNPDjPNic9Ee9oIx4M/Hh3Pm0M8G7cAd42sXVQz9MKvCuuagunZjQLbchl44HUXJcGEfvv3+1aYXhitW1MNuJfzn/jPoHo1mBbLbQl6sWuTH3Zd14B9uWYurVvK+TX88xEdHHu7jU9EChhYd2pGPc925Eigt64YAPATgbcbY93Mc06AcByK6QHk+80kH5xhepw1fvGa5rudMObhn+1J86uK2nPeHvE4kUmn1C3KsfxJOmyVn1k2xaKWbfqUvfqNSBQjop+qUK7VSUFeVXR0724ybmdBWx6oBwSKlGyC782E8NY0DPeNYvzjbAXDZrWipceO1k8N5H5sxxjtXmtRw1GuqY80mL5nRoZzri8pg9HIZel0A1vDaielTU3Zu5F12KxoDrixDX27pxu+yo8HvKlvfmIDbjq0dIfxm/xl0j0yhwe/SXRCcNiue+fwl+Np1K2CxEBYF3FjRUIU/vsNf68N9E1myjXhch9VyTkg32wDcAeAKTfrk9UR0NxHdrRxzC4D9RPQmgB8AuJUVMo1bkhPR70akIT5/qB8XddSWJYsC0PeBybRLdqnGQqvTH+2fhEPp/VIOzKpjRdfK8ks3mXWKi+ZMbW21LKn1wCIahml4W5mdm2uk4kc3teClI4N5G5xNxFNIpNKm0o3qCY7H+LzYIoqCROpuuXZH2ouM8bVL2blxC4YWqQHl1pDXRLqZhsNqKavDtLrRn7PQcTbs6GzAqeEpvHh4sKBA9vYVdXjt5DDGokkc7c9OrQT47j2sVMeOn83SDWNsF2OMGGNrNemTzzDGfswY+7FyzD8xxlYzxtYxxrYwxl4p36kvTIQWPzQZx4nBCE4OTaltfsuBtrGZyKFvCrpRrTQt0+bSHzg9htaQp+hsmFzwMX3Zc26B8hknQaaDZUIzdKTw7bPTZkVztQfHDB7q3lPmgVjBHRe1wue04YfPH8v52EPqIJlsY+Wy86EyfeNxRBOpguo4vE4bGvwuvDtUvhx6QB+MNXr0KQfPOKpvyCTitYa8ph69cRRiqfz9LWtzptXOhqtX1cNCPG5SSHxj+/I6pNIMv37tFKLJaSxvMN/xhpTq2IlY6tzIupG8d9Rq+t384R2uA25fXj5Dr21sJvR40WCtI+xVUyxPDkaw6+hgVopoKYSrsqtju0eicBrmAZSDjEefLLgXvZH2sDerd83erlHU5ZmdG3DbccdFS/DMW71ZuwGBWZ8bLfV+F9foCwzGinMFULYceoDXm4hJS0ZNPO3khr61ZYl6W1utFyNTSXXADlBc/UKhhHxOVeIq1+OJyuumHGnPWjYuDqLKZcPDu04C4APbzQj7nDg1zHsmnbUavWR+EJWSQ5E4nj/Uj6V1vvLr10obhNOjUQTcmc6CHeFMmt5Du07AbrHgzq2tZXtes+rYruEpNFVnt2wulaCq0Wdyu4sNCLaHfLxhmCYTaW/XKNa3BPOe719ta4PDasFPXjhuer9ZnxstdX4X+sdjmCowGAtkDH25cugFIvPG+NqRrw5JZsXypZk2ImaZN5Px1JxWhJYLMYSlkPiGzWrBpeeF1f5CuWJY4Sqn2iHzrJVuJPODCJaeGp7Cq8eHsX15+VNRRY65KJYSdIR9GJxM4MRgBP++uws3bWhU89/LgVl1LG+vXN4LGQC47VY4bBaMTSULHgxupD3s1TUMG4kkcHJoyjQQqyVc5cRHN7dg5xvdahxEi/DozbJuAKBBkbhyzVI1PVclV7xcgViBCKIaX7vOG/8Ljt347wjVZlIWM7n0GUMfKaJr6Hxyw9pGXNBagy3tuWfoahG77KagO6e3XqcMQQekoZcYsFstCHrsePrNXiSm09heRn1eUOPlhr5nNIYmzYDyjjrukX37qQOIJdP49CXtZX1es+rYruHy59ADPBhW7eEdLCdjKdgsBGeRAUHhJe/v4YHVvcooyFz6vJbPXNKONENmeLcGsaOpzlHNXO93YWAyjsl48R59uXLoBeJCYzTWvqogVmy6UndbSw0PYGvTSyfjqbJl3Mwl4Son/u/dF6m7kpm47DzugBkrYo2PKZAavSSLWq8DPaNR+Jw2bFqS3bWxVIIeO0YiySyPfmmYa40vHB7A9uVh07SxUqg3DNYYjyUxFk3OiUcPiJ1LUtWJi5WH1jUH0RR04+s738LxgUnsPTUKIhRUW9FS48HVK+vx/948rZN+AO7RV3vsObX0Or8L02mGqcR0wYHMFQ1+2CxUUGO7YhBGuhCd3WmzojHo1vW8mQuN/mwgXOXEXZe24yObsttna48RSI9ekoUojb94aajsefwA75vTNx7DWDSpG1DeVO1Wn+8zZfbmAX4B01bHvnKUD/Zeucif789mDe93kyhqupQWr9OGRz91AQDgjof+gucP9eO8uqqCH+vaznr0T8TVofCCN06NqimRZtRrDES+hmZaGgIuvPiV7dixOv/A72LJ5dHnos2QeROJp+A7BzT62fC316/E9WsW5by/TmfopUcvMSCCdOVMq9RS7XUgpXiZjRrpxmohLK+vQmeT37RcvFSM1bGP7e5BuMqp9uwvN0E39+gnStCJ28M+/OyTF2AsmsSb3WMFyTaCK5bXw2ohPHuwT73tSN8EDvaO44a1jTn/TpseWUwgszHoLnuHRN7Ko3DZq7WWN1gT5TSR+PQ5Id3MBdKjl+RFZN5cPgeBWEDf6dKYTvaj2zfi4Ts3lz0LRiBmxw5N8qyiD25oKluevpFqb8ajL+WLtqY5gJ9+fBPcdisuXlZ4h4+Ax44t7TV49sAZ9bYn956GhYAb1uX2BLWpg4V69HOFx2GD11m47NUa8mIilsJwJIF0mikjHOd3DfOFdtaANPSSLG67cDG+c1On6dDyclCt6cthbJfcXO2Zs+cFeOZN33gMT+49jVSa4UN5+v6USlAZsjIRT5a8db6ooxZ7v3U1blyX2xM345pVDTg2EMGxgUkwxvDkmz3YtjSUN5up1uuAcMznOzXxtgtb8JVrVxR8vBjU/sPnj2FKaYVdiRp9IbjsVgTc9jkfOgJIQ39OsqLBj9u3LJn5wFkiiomsFtLpiO8Fojr28T3d6GzyZw1sKCfVHjtSaYYzY/GypPg5bcV7plcpBWfPHezDnlOj6BqOzjic22bNFJB55tkbPn9JDW67cHHBx1/QVoNPbG3FQ7tO4OFdPONooUo3AJdv5nroCAAs3FdYkhMh3TT4ix/2USrhKhcGJ+MYnIzjWzeaDSwrH6Kvz+BkfN68yqagG51Nfjx74AxOj/Iq4GtXz1xtLAaQ5Bs6crbyjfetxPHBCL7/3GEA5etceS4S9jnVSWFzifToJVkIA6gNxL5XiFx6m4Xw/iJlkGLRjukrpqFZublmVQPe6BrFE2/04KpV9QXJSELaORe9YZvVgvs/tkHN7V/Ihv7iZaE5SzbQIg29JAvRHiDXOMO5pF4xYNtX1JlOWCon2oKk+TQ216yuB2PAeCyFm2aQbQTigjjfwdjZEnDb8fCdm3HNqnpVt1+I3LN9Kf7hw+vm/HkW7qVUkhO71YKNi4NqI6f3kvawFxYCbrugcN13tmiDzvMZEFxeX4WWGjfGoym1onImxKSp+Q7GlkJryIsHPr5pvk9jQXDufkokc8rOv9k2L8/bHvZh9zeuzln+X06CnrPDoycifOemNYgmpgsugLt8eR3e6hnL2fhMItEiDb3krOO9MPIAlw8Ec53HPBOFevKCNc0B6Q1LCkZq9JIFi91qUYOwxcwYlUjONaShlyxogspQ9fn26CWSuUQaesmCRhSHLdTqTMnCoCRDT0Q7iOgQER0loq+Z3E9E9APl/n1EtLGU55NIyo0IyM5nHr1EMtfM2tATkRXAPwO4DsAqAB8jImMp43UAlin/7gLwo9k+n0QyF4iiKenRSyqZUjz6CwAcZYwdZ4wlAPwKwAcMx3wAwC8Y588AgkSUuy2fRPIeU+2xw0J8tKBEUqmU4sY0AejS/N4N4MICjmkC0Gt8MCK6C9zrx+LFc18sI5EAwIc3tWBxrXfO2i5LJGcDpRh6s28Gm8Ux/EbGHgDwAABs2rTJ9BiJpNx0NgXQuYBL8CULg1Kkm24A2oGIzQBOz+IYiUQikcwhpRj61wAsI6I2InIAuBXAU4ZjngLwcSX7ZguAMcZYlmwjkUgkkrlj1tINYyxFRJ8D8J8ArAAeZowdIKK7lft/DOAZANcDOApgCsAnSz9liUQikRRDSTlljLFnwI259rYfa35mAO4p5TkkEolEUhqyMlYikUgqHGnoJRKJpMKRhl4ikUgqHGnoJRKJpMIhHi89uyCiAQDvzvLPQwAGy3g65wILcc3Awlz3QlwzsDDXXeyalzDGTCfYnJWGvhSI6HXG2IIavbMQ1wwszHUvxDUDC3Pd5VyzlG4kEomkwpGGXiKRSCqcSjT0D8z3CcwDC3HNwMJc90JcM7Aw1122NVecRi+RSCQSPZXo0UskEolEgzT0EolEUuFUjKGfaVB5pUBELUT0RyJ6m4gOENHnldtriOg5Ijqi/F893+dabojISkRvENHTyu8LYc1BInqMiN5R3vOLKn3dRPRflc/2fiL6NyJyVeKaiehhIuonov2a23Kuk4i+rti3Q0R0bTHPVRGGvsBB5ZVCCsCXGGMrAWwBcI+y1q8B+D1jbBmA3yu/VxqfB/C25veFsOZ/BPBbxtgKAOvA11+x6yaiJgD3AtjEGOsEb4F+KypzzT8DsMNwm+k6le/4rQBWK3/zQ8XuFURFGHoUNqi8ImCM9TLG9ig/T4B/8ZvA1/tz5bCfA7hpfs5wbiCiZgDvA/Cg5uZKX7MfwKUAHgIAxliCMTaKCl83ePt0NxHZAHjAp9JV3JoZYy8CGDbcnGudHwDwK8ZYnDF2AnzGxwWFPlelGPpcQ8grGiJqBbABwKsA6sX0LuX/uvk7sznhPgBfAZDW3Fbpa24HMADgEUWyepCIvKjgdTPGegD8bwCnAPSCT6V7FhW8ZgO51lmSjasUQ1/wEPJKgYh8AB4H8AXG2Ph8n89cQkQ3AOhnjO2e73N5j7EB2AjgR4yxDQAiqAzJIieKJv0BAG0AGgF4iej2+T2rs4KSbFylGPoFNYSciOzgRv6XjLGdys19RLRIuX8RgP75Or85YBuA9xPRSXBZ7goi+hdU9poB/rnuZoy9qvz+GLjhr+R1XwXgBGNsgDGWBLATwFZU9pq15FpnSTauUgx9IYPKKwIiInDN9m3G2Pc1dz0F4E7l5zsBPPlen9tcwRj7OmOsmTHWCv7e/oExdjsqeM0AwBg7A6CLiJYrN10J4CAqe92nAGwhIo/yWb8SPA5VyWvWkmudTwG4lYicRNQGYBmAvxT8qIyxivgHPoT8MIBjAP5uvs9nDtd5MfiWbR+Avcq/6wHUgkfpjyj/18z3uc7R+i8H8LTyc8WvGcB6AK8r7/cTAKorfd0A/juAdwDsB/AoAGclrhnAv4HHIZLgHvun8q0TwN8p9u0QgOuKeS7ZAkEikUgqnEqRbiQSiUSSA2noJRKJpMKRhl4ikUgqHGnoJRKJpMKRhl4ikUgqHGnoJRKJpMKRhl4ikUgqnP8P3pZNt+8tEnUAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(df.Z1)\n",
"plt.plot(df.Z2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is a sort of \"approximation Bayesian computation\" (ABC) approach, that includes X1 and X2 as latent variables:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Multiprocess sampling (2 chains in 2 jobs)\n",
"CompoundStep\n",
">NUTS: [t, beta, alpha, mu2, mu1]\n",
">CompoundStep\n",
">>Metropolis: [X2]\n",
">>Metropolis: [X1]\n"
]
},
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
" </style>\n",
" <progress value='3000' class='' max='3000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 100.00% [3000/3000 00:09<00:00 Sampling 2 chains, 43 divergences]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sampling 2 chains for 1_000 tune and 500 draw iterations (2_000 + 1_000 draws total) took 11 seconds.\n",
"There were 43 divergences after tuning. Increase `target_accept` or reparameterize.\n",
"The acceptance probability does not match the target. It is 0.5502104637542914, but should be close to 0.8. Try to increase the number of tuning steps.\n",
"The rhat statistic is larger than 1.4 for some parameters. The sampler did not converge.\n",
"The estimated number of effective samples is smaller than 200 for some parameters.\n"
]
}
],
"source": [
"with pm.Model() as model:\n",
" mu1 = pm.Uniform('mu1', lower=0, upper=20)\n",
" mu2 = pm.Uniform('mu2', lower=0, upper=20)\n",
" alpha = pm.Uniform('alpha', lower=0, upper=20)\n",
" beta = pm.Uniform('beta', lower=0, upper=20)\n",
" \n",
" X1 = pm.Poisson('X1', mu=mu1, shape=N)\n",
" X2 = pm.Poisson('X2', mu=mu2, shape=N)\n",
" t = pm.Beta('t', alpha=alpha, beta=beta, shape=N)\n",
"\n",
" Z1 = pm.Normal('Z1', mu=X1*t, sd=1, observed=df.Z1)\n",
" Z2 = pm.Normal('Z2', mu=X2*t, sd=1, observed=df.Z2)\n",
" trace = pm.sample(500, cores=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But I'm pretty sure that after conditioning on $t$, $Z1$ and $Z2$ are themselves Poisson, so the X1 and X2 terms seem unnecessary:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [t, beta, alpha, mu2, mu1]\n"
]
},
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
" </style>\n",
" <progress value='3000' class='' max='3000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 100.00% [3000/3000 00:08<00:00 Sampling 2 chains, 4 divergences]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sampling 2 chains for 1_000 tune and 500 draw iterations (2_000 + 1_000 draws total) took 8 seconds.\n",
"There were 3 divergences after tuning. Increase `target_accept` or reparameterize.\n",
"There was 1 divergence after tuning. Increase `target_accept` or reparameterize.\n",
"The estimated number of effective samples is smaller than 200 for some parameters.\n"
]
}
],
"source": [
"with pm.Model() as model:\n",
" mu1 = pm.Uniform('mu1', lower=0, upper=20)\n",
" mu2 = pm.Uniform('mu2', lower=0, upper=20)\n",
" alpha = pm.Uniform('alpha', lower=0, upper=20)\n",
" beta = pm.Uniform('beta', lower=0, upper=20)\n",
" \n",
" t = pm.Beta('t', alpha=alpha, beta=beta, shape=N)\n",
" \n",
" Z1 = pm.Normal('Z1', mu=mu1*t, observed=df.Z1)\n",
" Z2 = pm.Normal('Z2', mu=mu2*t, observed=df.Z2)\n",
" trace = pm.sample(500, cores=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And as long as I am trying to make it work, let's suppose we know alpha and beta, and just go after mu1 and mu2:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [t, mu2, mu1]\n"
]
},
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
" </style>\n",
" <progress value='3000' class='' max='3000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 100.00% [3000/3000 00:05<00:00 Sampling 2 chains, 0 divergences]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sampling 2 chains for 1_000 tune and 500 draw iterations (2_000 + 1_000 draws total) took 5 seconds.\n"
]
}
],
"source": [
"with pm.Model() as model:\n",
" mu1 = pm.Uniform('mu1', lower=0, upper=20)\n",
" mu2 = pm.Uniform('mu2', lower=0, upper=20)\n",
" alpha = alpha_true\n",
" beta = beta_true\n",
" \n",
" t = pm.Beta('t', alpha=alpha, beta=beta, shape=N)\n",
" \n",
" Z1 = pm.Poisson('Z1', mu=mu1*t, observed=df.Z1)\n",
" Z2 = pm.Poisson('Z2', mu=mu2*t, observed=df.Z2)\n",
" trace = pm.sample(500, cores=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Which we can recover, although not perfectly, from these 100 samples:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAD4CAYAAAC0VQLEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd5gUxdaHfzWziSXnDEvOcZGsoKASFBBFUDGnaw5X/VAMXCN6zdcMZkVRwYAoShARRHRBcpCcwwILCyxsmvr+6K6Z6u7qMGHDLOd9nn12Old3V9epE+oU45yDIAiCIOIFX0kXgCAIgiDCgQQXQRAEEVeQ4CIIgiDiChJcBEEQRFxBgosgCIKIKxKK82I1atTgaWlpxXlJgiCIuGfp0qUHOec1S7ocpYViFVxpaWnIyMgozksSBEHEPYyx7SVdhtIEmQoJgiCIuIIEF0EQBBFXkOAiCIIg4goSXARBEERcQYKLIAiCiCtIcBEEQRBxBQkugiAIIq4gwUWUfXZlAHtXlHQpCIKIEcU6AJkgSoTJA7T/E46WbDkIgogJpHERBEEQcQUJLoIgCCKuIMFFEARBxBUkuAiCIIi4ggQXQcgc3QUczwwt550AfrgfyD1ecmUiCMJA2RVcPz8CLPhvSZeCiDdeagc83zy0vOQt4M93gMWvlVyZCIIwUHbD4X9/Vft/1v0lWw4ivinM1/5zXrLlIAgiSNnVuIjSz6HNmimuNMMD2n/GSrYcBEEEIcFFhEdhAZCXE5tz/a8r8OmlsTlXUREUXPSpxDVHd5d0CYgY4vo1MsbeY4wdYIytltZVY4zNZoxt1P9XLdpiEqWGzy8Hnq4b2bGBQmD3UuO67QujL1NRIkyE8apx/f4a8MvTJV2KkmX778BLbYGVX5R0SYgY4aUb+QGAQaZ14wDM5Zy3ADBXXyZOBzb+FPmxC18EJp0D7PwrfnxG8a5x/Twe+PXZ2Jxr99KQzy+e2L9G+7/jj5ItBxEzXIMzOOcLGGNpptXDAfTXf38IYD6A/3O92sGNwPtDjevajQC636iZnz4dZT2m8+VAlyuAE4eAL66ybj/jOqD9xVoY8/Sbrds3/Ai0Gqxde8bd1u1n3Qc0OxvYuxKY9aB1+4BHgUY9gB1LgLmPW7cPegao2xHY/Auw4Hnr9gtfBmq00MrxuyIybeTbQOUGwOppwF/vWbdf+hFQvjrw96fA8inW7Vd8CSSlAn9OAtZ8Y91+7Uzt/6JXgX9MQicxBRg7Tfv963PAll+N21OrAqM/0X7PmaAJHMH7Q4FK9YCLJ2nLP44D9q0yHl+9GTBMD5L57k5g3Qz99+1AanXjvtNuBLL3GNc1PAMYOEH7PXUskJNl3N60H9DvAe33JxcD+aeM21ueD/S501hmGae69/5Qre4JwbX0Q2DTPOM+bnWv9+0lX/fs7h0Ir+79NQnY87f2zqs21baXVN0Dwqt7q77U/q+fCWRu0H7X6QAMnqj9Luq6p3r20bZ7pzmRdiNrc873AoD+v5bdjoyxmxhjGYyxjPz8OOytlQVKYwBEade4hMAS/xGnpsJYITSt3FJYl1wR766U1znCM4x7aEB0jet7znl7ffkI57yKtD2Lc+7q5+rWrRvPyMiIvLThMKGy/v80zwiefxJ4qg7Q+w7gvCejP594ro8dCd/vM2UM8M+PwJgpQLNztHIBRf+OwqkLYt9blwC1WgM/jdfGcJ33pPYM441YfQdbfgU+GgaknQlc83305SpOMt4Dvr8H6Hp1SAuzg3NNO2xxHuArPeZhxthSznm3ki5HaSHSN7OfMVYXAPT/B2JXpDLGjj+AZ5sAJ4+UzPXFdX//X2zPGygI/xheqP1nPi1QQ2bDLODwlujLFSuO7ND+i45d5oaSm9Nr4UvAzj+t63OPA9/fC+QeK/4yxRVhdLCWfwp8Nhr4+yPn/Za8TcEeJUikgus7AFfrv68G8G1siuORQKHWk/z1uRidL6DZ6QtyozsP59YGef5E4ORhYLeuae5YApzMsh5bVOTHKHTdTESCSwp04Kbn9Nlo4NWu0ZcrVgTEwGO9zH9/DLx9VsmUZc4E4N1zgfU/GNf/+TaQ8W5sOiVznwA+v8J+u9CuS7uJ1xEPZRcdlmP7nPf78QFg+o3RF4mICC/h8J8BWAygFWNsF2PsegATAZzLGNsI4Fx9ufgo0B2hv70Ym/Ot+hL44b7oBeG3twGPVzOtFB8L08ZAvXce8Mkl0V0nHIpKcEUSXRYUXH7Jd2TYIaoiuSLnIHTD4uMqBXx+mbFBDehli6QTYea354H1TibACHx8OYcjLk5MERGhXt6leJbMX3TlIaLGVXBxzi/jnNflnCdyzhtwzt/lnB/inA/gnLfQ/xdvDRWNpi9GGavy9ASqOQejO8/yT+23MRbSMorT5JR/smjOG5XGxUKNbnEi5yB0I6g5lzYNQxIgQS0oymc5eWAYO3t8HpvnAc81AVZ9Fb4Ae7wGMOOu8I5xIvicPOwr3ruPBFdppvR4H8OhME/7749AcB3cCHx9i6b9CEQlLYretWxaMZsRBQtf0j5WmTVfa+bQrG3RXT/SiMJNc7TrH7dxX0YjuMCtpkIz2XuBWQ/ZPzPB359q2dtjTWnUuAAYWt9YCa5df7nvE24gzi7dND7tek2AhUMgH1j6QXjHOBJGVKGo1/7EGF6fiDXxKbiEL8qrxrX0Q60RzjuhfUgrpgD7Voa2C1NCIKAJigmV1eNSIkIyFdo11nMmhHwqgpX62BPz+JRwiVTj+uMt7f+ev9XbIzIV6s8iUOgukGbcBfzxOrD1V+f9vr1Vy95eGANzmYyd4IrWDxouZs3UUJ4S8Dt5vVZJDdhe+qE1U0g4QlfUy0hMhftWlR7zaBknPgWX0Li8Cq6FL2n/j+0LaSBJ5UPbmaRxCUERq4ghOWWQW2P95bUhs02snOGyj6sgD3iumaaluCIJXBVuGlfucavQFPe//XegwEWgivN7NSk+aTuU0MiJg1rHZK1LPJGd4DJH8M1+DHi7n7drCz4eqYXZe8Hc2ZHLExQOHDi2H3i+FXBgnbfzzrhL298zYY6FUpnazPkC83KAN/taBxdHQl6O1qmYcaciU0gY35J43k5ti+o8h7cCb/XVOqFEkXN6CC6xX2FeqDGVjw06bwtD66NxeCs/EOaeaXzN9JDZJhyHshOyhnAyS/Pjzfy39+Ptyiqez9Hd6nDsZ+oDr3Q2rhP3suhlLdtBcL3ieYVrBnMzPQoOrNX+L3nH5Xw2gstcLxa9DOxd7u3ags1zvc/vZb6eSnBxDqyfARzfp4Vpe2HpB9r+XjVn8T52LNYicJd/BjzTSE+6fCJUBzjXMsSYTcyrp2n5Arf+Flq3ZxmwfxUw+xFvZXDi6brAazbDnMKpS+J5O/m4VM9MZN7YMt/9GkTUxKfgCtdUKPbLPR7SuORKLCppoFDSvhwawqO7nX1Hqg/Ei8Zl3t/uXHasnq6F2xvKIl1z5r36uRWv/bs7NW0seJxL7zRQoO3zUlvgnf6h9ScOar1PQGsYDWWR7kUOUFHeY4z8N2a8TlMSKNTub8VU0/oYmyTdMNcZeVnWyoWpVNT1CZXVqYbMPFHDfR8zP9wH/Ph/QO5RIDcbeLoe8EwDbdvhLVp+xD/eMB4jtCr5vYtZBmTrh4ybVmxGhLIDmnUhiElbPLoLeLULcGSn9RxeBFeBlN5p5Rfa9DziuITk8MpMRET8Ca6CvFBIsFfBJYI4crNDGpdoAHIOhxppHjAKMTteagt8cIFxnexjESbH9TNNwRliHw82d9kM5JWvrtXC7WXkhlaEO6sE17IPTVGVLqbCwnzgP3qylEObQutfaA282ll9jCyE7H4LIrl/LwR9GC5VnweAbQuBQpNPK5zOh0xhgWaiDTea0pPGFQj5SP1Joe1y5v1Dm7W8d7FCfCfznwmtK8i1F+xi/8NbQhqLiOZNTA3tJz+fL64KbwiDjFyXzZ2UZR9p5fj7E+tx4vpObYssuKbfqJkIw/m2iaiJvxmQPxuthdoC9pUre69xOahxHQs1RG/0AC5+VwvWEAlfeaEkuAq0rBe12gApla3X2LPMuLxpdug3D2iV2YBDcIYK2QzkxN6VQIVaQMU66u2qhtZLKpugb85me6AASqFiDjIxnlT6KTVQqjIW1YBX0XHxIrhys63rRQOVexz4aLjzOU4e0RrHnrdogSazHw1fgzTvb3gekibhFmn7v65aPR63Q73dFVNFEN/Uhlmhdaey7e9P7J/xrvaNpV+rdbQAIKlCaD/zN2LuOHhl/UzNKtJXSm4sBwfJZZIJFNhvC5bJVMfzc6QOEQmu4iD+NC4htAD7j/TF1sZl0UiZzXuix5Wj90R5AMEPdNtvwHvnA1OvtC/L9/dqIdsTGxt7YXYNcTiV283HJRrgt8/UzB52OJnh7JhxN7DlF+d9ozWZedW4whH2XhC9fFfBVagWmqKs234LZUOx4/PLNbPZvlWhwIQTDtnR9q/ROiIyThqX6AjIGVt8DmHcpyLIV1iQq2V3N5dD1OEcSYsT5mMVsiDYusDo40uSNS4H02g4/HAfMOcx7ffXeuZ+8eyCwklRB+w0p+MHQgFbqvoaHLgcf01qPBLfT1l8DCcOqe3VAlHRzL03c88pELA2lE7h6Bnvaj3pU0eM57ITFnaNsDn0NnuvVXDtWhpKFbXicy1B7cGN2rJTdgylEHV47ZwDS98PLX8yUru25bxhCq7CfOMkkgbB5VDGWM//JIIINs8Ffn7Yfj+7DkOgQIv8NE8x8/09ofchEPNAgYfqnlPY/pu9tY6I4Xo2UYVHdgI/PRRaJ55TrMcf/fGmFjyRYZr2RLz/fKkzGMh317gAICEFSK4UWpZ9XOa64KWeOc3IbRi+IAS9Q8i72Ga+j08u1syCJ7PU9ZVMhcVK2RBcL7QEXm5vv5+oVOZG0CzIeMDaUHj1o8nntpuwTrbfH9wIPN1AGzdmdgS/2NoouDgHJp8DfDRCWyd8VXIU1Z+T1NdUfWROjmdVQzFPMRdUuALFXD63gdmiRx9rwSVrHU45/uw0h0CBFvm57jvj+oz3tPdhOE6YpgKhYIFwBb4lHF5fPiaZw3OPST4uF8H1w/3Oz9TsgxMa6prpxvU5Cn9ZYT5sfZIWwSWZB2VToaOGqSAvx3lGbvk5cel9AOrvQDwb+bpZ20PjPjl30bhIcBUH8SW4Dm02LqtC11UNjqio5sGj5owBvND64XgWXNK5P73Yup1L2lzBKa2RyzumRQIqNSApqk58KE5h1z/cZ1w+eQTYvUwtFE5k2g+ULDhlXbdlvjWEPtwG+JQpO77cIDuZMyP1cdhxfL9xeeUXwLNpVk0oUAi1D8/lvlXP+9Cm0H3Ix2+aqwUfrJgKbF9sPY5zYOHLpnX6s0osJ93D5yETo5OpENAGazsNrreYBMNoIgKFDhqXJCQSko3CKiFZy9SSvUdhKnQbL+iSGV8+35rpRquK+Lbl8YbCVyjfh/DFifXKNkY/55GdWsRlSaQ0O42IH8F1MktzMMs4OVcF028OTZdRmGfdX4YHrMeb/WjmmVIFTuYKu3MDWg9NJSzk4AzzxyxmEnZiymhg0tn2wRKb52lC00yBzTP6a7Ix/ZQ5o4Vbbjnzx+7Vx2WXqeLgJvWswG6Y398P92l1yxyIwQPqa7tl6DB0ovT/X98U6nQtlBJDfzISeH+Qtv39QdZzZW3TZh42lwuw1gmhEXjpaE2/wX6bXP4jOzwOVhfH5tv7pORyJZYzCkQe0Exx7w1SjJtz8XG5aTjmb2v/amN2jL0rNLP7Ot2KERRchZpJfsMsY50JFKjLFAzaOQoseSv8sX1EWMRPVGHuces6X4JxrMrxTGMqJ0DrjQrcBNemOdoEcuZr5B4Dkitq/i5LtKDOis/V6wUqMyRgn4tw+Seh42TtJON99f5mhD/JbrzZtOu1/+1HGtc7aTj714Z+i2wkgqUfABe+4lCgMAWX6KGrhDoAfHoJkLXV4Xo2ZO9Wr9+3CvjmltDyTw8Cw9+w7uemARpMe9I92zVk8lACM6p3t/h1rcE9ZoqcVZm4IkEu/3uD7J+XisJ893B4QOvEHNkeWhZa05Ht1uPNyyePaPsf2QGsnKoN43DCXH8KTkmug1xgj/5epl6hTbZZIGlcIqhDHmIQKHA2FRLFQvwILlXPypcgRb9B0zCOOgRpHN3lfp0fHzAuH96iDa689GPnqKz9LjkFzQJI4JZM1Czwvr/bfl8Zf5LWA3ZLsjuhMlC5UWj5xTb2+yZXdD6XyoSStQ0oX0sR1u0SDi9MYVnbtPRI/R80+kVUORj3r9ECHG77C6jZUltnDiow514U73Thi9ZGWiWk3HIVinvJy1GH04eDqjFcOdW6Tt437wSw4cforsm5ljopHKEF6DkobRpwuW5s/Mm4TUx2mpCiDWA3HGeqG2/0tAptJ/JNgitQEDrnzw8Dw183bhfv3NCxMo3F9CS4StusAmWL+DEVKnvlJrnrJLQA52lH3NgyP7oJIDmP0O7NI+vNCSe9l+zwRz2O7XF7vis+s657pZMWFu5kKjwmmWKy9J64MMn9NVkLnRYaXiAAbJwDZcPws546SAROBAq1aD8vqFL1zJ5gXeemtU8eoGkDGe96u64T4QwFEEEU858GPhsT+TWFn2qZywzAKnYusT/OqQ4L/yfnwFt9rOWRCUdoAVaNa98qYxlFMmmBeL9236qdH8+cKUREOR/dpc3T97vHFF+EJ+JH41JFQkUyrUmk+PxRCi4bjcvTcREIPJF6Jk9hYo0U2ZQWzvYtvwD1bPIWAsCkc0K/PxoG3LXC6psTGlbGu9ZAFMHmudp/oZ3bmRm9kqvQsN0E18F/tOzk1cOY+8uOSMcwRXvNSM1ePz3ocF7FORNStHckNC6VhvvRcKDDKO2vyZnW7W6Y64DZomK2lAhT4Y82U+XIGpvMwX+My19eDey+UzNlCq2+x7+Kt80qw8SPxhXueKRY40uIbjZhOx+X63GK4Aw3vrwmFD3ndT6uWDS0TjhpXDIiLZG5oyIaC9k34kZRTEHiZRDvySPR+ZrMGR6KE7vgg6jP62AOdjKp5udojf+HFxhzEXolnLnaOHf3YdqZClX8/qqxvkQyAJxQEj/iXxUdV5yT/P3xBtDFIYuGG3ZRhV6OC1dTW/N16LcqqEWFW1hxtHh9Vz6fJnDM46QKcoFpNwInw5jvKCaCi8FglvzuDvdDEpKja/wDBZqptyQc/mumF813pYyo9QOwiapVIaeX8ko4frpAofu4wXAEl5mTWUD56pEdSxiID8EVKNTyoFnWF/NYCbtJFb0QjakwmgbMq6nQPL4p5nidgNAPLHjeuv7AWs2H4u0k2r+VLpGeXkgsF76mnZAS2bsWLHhe8+tVdBhYW1Q4ZROJBrkzJfD5NauJ3RAMM3bmu1hRmOs+EWSgIPJ2JxpXA2EgPkyFP9ynHueyYWbxlsMudN0LEZsKA8BX10V+3f2rIz82lmya620/n18d/elZaEEzj+afjM2kfgkpkR0nT98RLotf1zochza67xsLGvUu+muo7iUhWRNcKvNcbYdMOEXFqWzrBKfJlYzWHqdB1m6EYy0gHIlKcDHG7mKMrWaMrWGMeYzTjgB52oNIOf9p933ciCbQYc/fkVX4vJzwGu3SyoG11nWpirmgTmSGgiwi5bfntUGldoRTn+QMFV5Z+Tmw8efwjxPkFbHZ1kys8xt6JSHFXuOqmlbsxVH6T5MqmKZdicJU2PTsyI4jLEQsuBhj7QHcCKA7gE4ALmCMtYhVwQxE0niYkQcRlgQLX7TOCuuFX56MfVlKC3adiaI2W8rphtwoDRMDNh9YtOe3+zaK+rpOGlesp7Pxwp+KWbGZTzNf1+2kLa/9Dvh1YmTnTyjhNqgMEY3G1QbAH5zzHM55AYBfAVwUm2KZiIXgKg0NULgDOss6lRuUzHXD8T9FaiqMJbEwmyU4fEN230Ys7r2bg5nbn6wFZqgS9kaq1fS8DWg1JLJjV0+zrsverWnA5fRJU5e8qR7zRxQr0Qiu1QDOYoxVZ4ylAhgCoKF5J8bYTYyxDMZYRmZmhLOZmk07o1zSvKgoaY0LUJvLwqUkHPZFRWq1krluOL7G0tDhsZvaPhzudah7toIrBvfulDsxIQW2QTs9/xXZ9bpdG+OoSL18XrRPc7o4osiIWHBxztcBeBbAbACzAKwAYAl/45y/wznvxjnvVrNmzcguZta47IRQ3c7qnuWAR0uH4FJFVoVLhVrRnyPWpEUwMBSIjSYdLrU7hNewRTJ2KBpqtLKui0hwSSnSRn3o3Emwew+x0LjkbPVDTNGiTqazShFq48xXNMMIKtWzrqtpSo9G+QqLjaiCMzjn73LOu3LOzwJwGEDRhEGZNS47Z3L5Gupt6dcWcc85ijl4KtS2rqvfzbpOEGvTlTyhX6REGjzjZL4qKloMDK+BUZmxVIyOIp2YjKqeehVcVRpLC5Im4zT/GmD//mKicUnXNnceneqyW5ntYL5QdpMOl0Z2DhUqv2jD7sZleXZ2okiJNqqwlv6/EYCRABTJ6mKAV8HlS1RvS62m2dNVmJNsRkLEDXeKt6lZZGKtOZarEv05Em0aoMoWy7ERr6bC8i5aZh+XKVVkmg0IT+M69wlv+3k551ki3ZBDR0eVDcZrMIldEmS3qU6qNlGvt/tmwkH+HnNMCXSdBGM0gktMMVStaWTnUJ/YuqrrVaHf18+J4bUIN6IdxzWNMbYWwAwAt3HOi2aEXUpl47LdZHk+P1Cng3GdsDvbmSW6jI2ubIB9w62ifC3g7tVaBvO7VwHtFPEsTsEDnRQJVM+OYtBoagQj+c3ObzvNqVZb5/N4DcM25zk0IyK+ZDqMAu5YZl3f5MzwfFzlbczbaWcCF0iTPHoRXOWqAAMnAHcqyiVQNdheBZddCjTVFPWyttPjZqDf/yn2ibGPK7kScIY0F5ijxhVhbgSfPzRspbw03GLAo5GdD9CeayWFb5n5gHpdreuJIidaU+GZnPO2nPNOnPMoB984UMtkS7bTOvyJwMXvAXU6htb11xN/xqL3aIf5A+x1u30jUrstUKWhNu1GhVrAuYrJEOWG9SxTUtC6nYA2w4BkSZj3u995jMhYxYSRgubn2m9TcfG71nuzE9yqqWjMY7eSK1v3sZzHpZqqGuaUykD1Zur9nToGTc4yLtsJ1/YXA52vkM7pQXAxP9D3HqsmUK+LcR8zXk2Fdo29WRgOf10ToPL2jqOtx8VCcMlpt7pdBwx9wf78D0pRt5GaxJkvlJ9T7nS0V8xK7pWkilqH+IqvjOt9/tA9uCVfJmJKfGTOKF8DOOPG0LJdhmVfgpYLTBYG4qONRWSWHeF84OaGW25UWl+g/Zf9TkkmMyQPAKM/Bh40BQ1UaQRbmg9Qr7/mB6BRD+fyCoQmWzXNOsZGZSpNO1MtcG5ZZFy+YTZwVpSpfHx+4OrvjetUnZur9ZmjnYTM1abZpZnPem5xDlmoeRFcqncNGAfgqoS9/HxTKgPtL1E/W5XWBFgFWpex1jGFquelaoyvDTNfYOZ6qRwmAWruTMrzrbnN/WYH84fKLZui3Xy5TgkKRLka9TRdyxc6LwVmFCvxIbgAYOjzQDm9IvqTgDGSO018xMKEKH/UogcrfxSxyMQhY+4dOmkITj6l0Z8A5z0FjJJmOW7Y035/mXBNfjVbA2l97FMImrW0YMOmalglU+HYadpMstd8b923dnugoimjRc1WwDkRmjqFZs38QP10zRnf/hJtnbmxbtzXqk3ZIUe/+RKAagofEA8YhYxZmFdwyNwBaJqr6MHLA3BVDaD8fCvWAy55V92wtxqk1thU5sfGpjRPKm1NldGiYh0tMlPmjmVA95ut+wLOCXSdOnyRZvNgPu259r7TGPiU4uLLlcvfzjQruHjWqu982P+A3ncAaTYzoxNFQvwILiDUUPgSjQJANBpCE5MbFCFEkqQPvU+Ms1OZP0BVr1ngZLtnDOh9uzH0tnGv0O+Oo41mUBkxQNIzehkb91abGWvYJEFhgEXayT6u+unSvqbnIO69WjMtSMIrdlkURIPMfJpmevEkTSCLdYZzyOZBlyhQ2cTkS1C/M1GmIc8DN8yzmh99CcCVpuEPstksMSUUlSYHl6iEhVzX3Rr0+zcB95jGbKk6Ui08mIhbnq81yuZz+Uznq97M3qLhNK1OjwjHapnpK00W6vMDdTsC5z0Rela+BGuZzcjC3dyREgLLXA+YH6hYGzjvyciDSYiIiC/BJfAnmiqR3ogEgzakhsmn0Lj625hUIi5POL6ACELnRS+6+032QtGr4Lp5gX5OIdBTgau+0X4bgi4cymk2i8kagWxyMjeYYtudy4ArbfxuIhBn4H+AQS6pdcRzkRsNUWzztWW/4U3znc8r1y1fglqLEXWu+41Ag3TrM/H5gAamcGlzAteUypp22lEKuCnMBVoO1nxBQWtBRWCEPlNvsJNk835SqwGV65uu62F6F5XgSSynNcoDHgutY77Q8xnyPDD4v3oZ9e/L3LiPeFPT0G6abz2/XE5VgI1X5Ocsv/egBUZ/jrcqcn42PVuLqpS/K3+iNvZNdMJE3TF/e+Y6drdpUkqiyIgzwaVXHH+isbESpr82F+q7KTQuJ7NEuAEKZnb+YVyONOWMHaJRdOpti4/MLrw/tYbuV9GfjfkjHL9fM1UK7MydHFpwiIwcnCGbYd0+dBWFuqms02XuiVZFIymPXxLKmZPGVa+zs0/QILj8Nr1p072ZfUS+BOtxdgJE3q/gFHD551r0ndwZaXOh1sCK8PwLX3G+B8N1bcx1l34Uqi8plazPTJTrzHtD6xgDLngJaNJPm5+ux03aeqF1p19rPEeNFsAtC40BKCpEh+qCl4HLvzRuG78fqFRfewajPrAea+i4SPcQ1Lj07bVaW4+94ivgjqWm8yUC7UYAPW/Vlu0CL8zPy+l9uN0/ERbxMR+XIGgqNJlv2o4Azhkfyn1n8HGZKpdqRP7oT4CpY4FNs72XpcOlwKovtN8plY2zm5qduFGjt8YqO71otGu3BR7YqjV2a7+1Zid/YLP2PzjdhqnhNUcGOgmZLlcA394qHasLq7Qz1Z0GgZMDu/NYrcxizqKkVB6I9tAAACAASURBVGsZzVz4CrB9kRahKRCNjLnnbw6Bd8rhKs8G7U+0CuDkSkBX06SibUdogQhrvgYObdI1NdP955s0LoG8n2wqvPZH7V0KH8tdy0Pb2o3Q/iZ4iMq0E1xthxuXKzUAjkpBPyoTKfNp2tHV31nXA96CVC773D4DTDdJ8N36h9YhSEwJpazaqBgvZdCyJCFm1rhUqAK9xH0LE60suO7fArzZGzi+z9klYC7fJe9525fwRHxqXJbfMCVstWk8b10i9exeAq7UTWSJKVbzihvVm4UGIP5rofO+Q563d16Hg9kcOH4/cHtGaDm1mvYx9R8X/bVsxwQp1lVuoAn/MVOcdz62z/56I14H7t8YioA0BNBw7cO/8JXQqpqttR70GdcbzyMaabMgtoTAO0iuGs21RrPZAKDBGdZnMexVRRqyBM03IjoXKhOjncYlN4ByRGzNltpQh2gxmyzt6HOn9r9xH/36Cg3FrjMh/HV2EawyrQYbfaF21GpjHdKg0n7tOktivbzuxnmhNsAOIcxEBLCw5ABa1LKI9PWSwf6GecC4nTEeDE3El+BKv1r7b46ocgqdlSt6rdahqbO7XQc0k4MS9Epe3ePMLDwADH0RuG+TZiK40SHdS/cbJW0wivRQ5vtMTFEPrHbKriE+NrdymBtrMVC6kkLA+5O1jzvFFHJsPke2YoJIM6M+AO7829pAtb8YSL8mtHztj+rjheAyD4q2m7X2qu/UEZm12mh+uMRy1ufpFJUajHD1aw3g8DeAM/9tLJsTnRTjqaLh4UzvnbLuN2o+t8unAveuU5vX7To09btqDXTroZEPHvaCUgv0q38Lv10vyTpQP93dnyY0tXqdgetna/5Ww/XC0C4bpBv960RMiC/BdfZ44OEDeoWUejvmUfGqcHiv9Lwl9Pu+TcC//1Hv1/J8zYxUQY9As+tBiqSpXvw7bngVenJDa/YTiY8tXMHV63ZtgKg5nB2wz0oiriEc+F5ILOfcOx23U6sDdumihDnO3OiaQ+GFAPclhMxzl35sX6Yb5oXCq8s5pKoSAlc0sF2uCIXou0XyFUWvPJJot+SK6qSygHM9Fh2XcL+5cFD5eYWGCBjLl5CsCWKvFogRbwI9bjEGqjTsbo1IdBJcD2zV2g2iSIk/H5dokETDU6eDtUfj5GcJhwo26X4mHFWvT64U8m9d+hEw94lQNFU0mtatS4Csbd73lxuOm38zbhNCzW3cl1ze2/7SluXn3PVqYJk+vYxdlgPx7M2amFdq6ymjOowKrXM7V9BUaNK4zjPnHBSCyx86xqkn3kDqmDhFcArBKde72m01s65TarCbF7jndoyEWHSYxHnMY9fs8CWoJ4eMBSqNS+44RROW3vly7c8N8X2pMrCIDtXNv0WWTo3wRHwJLs8owuGjYfjrmr3/1BFnM8iDO0O/21xotI1HM6NrrdbqiCg7RMM66FlrQ1+7neZzU+VIBLQwYDk3XkplY/CDoNXgkOCyNU3q74EHNB9VuFNVVGkEPJrlPgZHJmgqNAkJc089aDL1A4F87bfbwPRcPQeek+DaqwdQ7DZFqrnls4w0HHzoC8DMf9tvlwVN84HWcnklqSKQexSehnNUrAMc3hzZddxw+56j6SB6pWF3IHOdczaOujbjLYmYEL+CS6jzqpDpqDQuhYCJRSLe4qR8dc23oTLhMab5MuxoN0L7f0LP5G2X0NhsknHah/PIc8WFI7QAbQbcDT8azUcquowFFjxnFEJu84M1PwdYN8NbRv2KNqa2WHPGDc6CS2asYoZfr1SqC2Qe9ZaT76pvgU1zjBGCsUJ0HBPKWcfFFRdD/qsFBVUpAg2Z8ET8Cq4aLTRznDK5bIxMhfGM0yR9XlA5tg3IAzbtBJekcRUXjXsBj3qYQ+uc8VrDKvty3ATXyElajj8vWrxTsE6seWCr9RkPfRHIiGEI9hVfAWumq7Okm6nSsGiEFiBZPEwdzJvma8KyOEhIjm7ANBE18Su4AOs4FIEhaj4CwXXVt8DmXyIqkjseTRn3rgdys4uoDB5ILGfvywNMGpdLcEZxCq5wMAcguE64WA6o2th5H0H5Gu77xApVoMoZ11uHCkRDlYbhzXvmlUET3fMIytiZ6ut1oUG+pxHxLbjssBuQ6IScj6xpf+0vpoTp46pUF4CH3m1JwbxoXGGEDRcV3a4Hju503ueKacA/NuH1kRJpktjTDTmK1wvie47GZyzodFn05yBKhLIpuJwGKttx9kNaY9PJQ1RRrLlsavxFIBkEl53GJToQMWhkIuWCF933aTFQ+yNKP3amwnCZcDQ2wo8oEcqo4JLwGmWUUkkRMl1MtBpUMteNBjmrgl0AhdDETlc/IxF7Yjm4uTgiEIkioWwKrmgGIBPeqFTP2QcGaAEQPj+ZZIjYIQQXaUunNWWzKywc6B3HuI+fIYqOlMrAoGdiMwV8vFBcYfCnK0WZToqIG8pmLRDzHBFEcXNHhrexTkRkBIOtSOM6nYlK42KM3cMYW8MYW80Y+4wxRuqNG2RXL9sklY9gNmrCM2QqJBCF4GKM1QdwJ4BunPP2APwAxjgfdRpDHxpBRE8wkwt9T6cz0ZoKEwCUY4zlA0gFsCf6IhEEQdggNK7a7YD8U6apiYjThYgFF+d8N2PseQA7AJwE8DPn/GfzfoyxmwDcBACNGnmcapwgCEKFzwdcPQOo1S40tx5x2hGNqbAqgOEAmgCoB6A8Y8ySjZZz/g7nvBvnvFvNmjbThBAEQXilyVkktE5zognOGAhgK+c8k3OeD2A6gN6xKVZZhGzyBEEQsSAawbUDQE/GWCpjjAEYAGBdbIpFEARBEGoiFlyc8yUAvgKwDMAq/VzvxKhcBEEQBKEkqqhCzvljAB6LUVlOD2gcF0EQRFSUzZRPBEEQRJmlbKZ8Ko10vRrYsQToXQST8REEQZxGkOAqLspVAS6bUtKlIAiCiHvIVEgQBEHEFSS4CIIgiLiCBBdBEAQRV5DgIgiCIOIKElwEQRBEXEGCiyAIgogrSHARBEEQcQUJLoIgCCKuIMFFEARBxBUkuAiCIIi4ggQXQRAEEVeQ4CIIgiDiChJcBEEQRFxBgosgCIKIK0hwEQRBEHEFCS6CIAgirohYcDHGWjHGlkt/2Yyxu2NZOIIgCIIwE/EMyJzzDQA6AwBjzA9gN4CvY1QugiAIglASK1PhAACbOefbY3Q+giAIglASK8E1BsBnMToXQRAEQdgSteBijCUBGAbgS5vtNzHGMhhjGZmZmdFejiAIgjjNiYXGNRjAMs75ftVGzvk7nPNunPNuNWvWjMHlCIIgiNOZWAiuy0BmQoIgCKKYiEpwMcZSAZwLYHpsikMQBEEQzkQcDg8AnPMcANVjVBaCIAiCcIUyZxAEQRBxBQkugiAIIq4gwUUQBEHEFSS4CIIgiLiCBBdBEAQRV5DgIgiCIOIKElwEQRBEXEGCiyAIgogrSHARBEEQcQUJLoIgCCKuIMFFEARBxBUkuAiCIIi4ggQXQRAEEVeQ4CIIgiDiChJcBEEQRFxBgosgCIKIK0hwEQRBEHEFCS6CIAgiriDBRRAEQcQVJLgIgiCIuCIqwcUYq8IY+4oxtp4xto4x1itWBSMIgiAIFQlRHv8KgFmc80sYY0kAUmNQJoIgCIKwJWLBxRirBOAsANcAAOc8D0BebIpFEARBEGqiMRU2BZAJ4H3G2N+MscmMsfLmnRhjNzHGMhhjGZmZmVFcjiAIgiCiE1wJALoCeJNz3gXACQDjzDtxzt/hnHfjnHerWbNmFJcjCIIgiOgE1y4AuzjnS/Tlr6AJMoIgCIIoMiIWXJzzfQB2MsZa6asGAFgbk1IRBEEQhA3RRhXeAeBTPaJwC4Broy8SQRAEQdgTleDinC8H0C1GZSEIgiAIVyhzBkEQBBFXkOAiCIIg4goSXARBEERcQYKLIAiCiCtIcBEEQRBxBQkugiAIIq4gwUUQBEHEFSS4CIIgiLiCBBdBEAQRV5DgIgiCIOIKElwEQRBEXEGCiyAIgogrSHARBEEQcQUJLoIgCCKuIMFFEESR8cb8TXjx5w0lXQyijEGCiyCIIuO5WRvw6rxNJV0MooxBgosoMb5bsQeTf9tS0sWICbkFhcgvDJR0MQjitIAEF1Fi3PnZ33hy5rqSLkZMaPXwLPR77pdiveZTM9di2tJdxXpNgigNJJR0AYqLj//YjpoVkjGofZ2SLgpRRtlz9FSxXm/Sb1sBABenNyjW6xJESRM3GtfuIyfxz/5jnvd/afY/SBs3EwW6+eaRb1bjX58sLariEUSxwjkv6SIQOpMWbMGurJySLsZpRVSCizG2jTG2ijG2nDGWEatCqegzcR7Oe2mB5/3fXrAZAJBHfgeiDJJ9qqBErvv6L5tw6duLS+TapZED2afw1A/rcPV7f5Z0UU4rYqFxnc0578w57xaDc8UMBgYAoI4pURbJOpFX7NfcdOAY/vvTBvy59XCxX9uOtHEzMfHH9SVXAK2ZwdGT+SVXhtOQYvVxbck8gdGm3toFHeviyl5pOJlXiGvet/ZaLklvgFHdGgaXzceP7dkYF3aqhz1HTuKeqcuD63MLCgEA89bvx4Wd6tsef8c5LdC3RQ2s2XMUj89Yi9yCQmSdyEedyikAgAcGtUJ642pYuv0wnptlHY/y6IVt0a5eZSzceBD/m7fRsv3pkR3QrGYFzFm7H5MUEXQvje6MelXKYcaKPfjkj+2W7W+OTUe18kn4MmMnvlI44j+4tjvKJfnx8eJt+H7lXsv2qTf3AgC8s2Az5q47YNiWkujHh9d1BwC8OncjFm06aNheNTUJb12ZDgB4dtZ6LNueZdhet3IKXh7TBQDwnxlrsHZPtmF705rl8czIjgCAB6evxJbME5byCe7+/G/sNfmIujauiv8b1BoA8K+PlyIrx9hY92leA3cOaAEAuPq9P3Eqv9CwfUCbWrjprGYArO8d8F73Dp/Iwy0KM7Nc9wTydW48sykGtq2NzZnH8dD0VZbjbz+nOf7aehhdGlXBW79a64ZT3TuZF7rXSOvevee2xDfL96Bn02qYsmSHZbu57skCa/Tbiz3XPbG/TDh174nv12L17qOG7XLde+vXzfh7h7FuutW9tvUq4bEL2wGIru75mSa5snLyDfdY1HXvdCdajYsD+JkxtpQxdpNqB8bYTYyxDMZYRn5+8fdKAmFaCtfuOYbth3Nc/WlZOXn4c+thS2MZDpszj2PYawsjPr4kyC0IoCDgrsaezC/Ehn3GZxjgHJnHcouqaKUe81PbfigHr87bhP/+FP4AXfO5CgIca/ZkI7fAe4V/5NvV+OzPHdidddJ9Z8U1i4N56/fj3YVbcUxhGg14qIdFjShBYYCTn6sYYdE4eRlj9TjnexhjtQDMBnAH59zWEdWtWzeekRGZKyxt3EwAwBMj2uOj37dh9r39HPfv8NhPOJZbgKUPD0T1CsnB47dNHOp4XLtHZ+GE3pt12vfcF3/FxgPH8dPdZ6FVnYrh3EoQUaYbz2yCh4a0AdN7b6WZtHEzUSU1EcsfPc91P8D4DB/7djU+XLwdH17XHf1a1vT8TuIBt3tZtOkgrpi8BNNu6Y30xlUBAH9uPYxL316M9MZVMe2W3mFdb9Wuo7hQ7/RsmzgUHyzaigkz1tqWYefhHOw4nIM+zWsE1w199Tes2ZON7+/oi5REP2pXSkbFlETXe3S6T7tjIn3Hj89Yi/cWbcX4IW1w41lNDdtO5Rei9SOzojp/tBzIPoXuT88NLhdVORhjS0ubO6YkiUrj4pzv0f8fAPA1gO6xKJQTj3yzGhsPHA9GC9qRq28vDLNX5nV3IWO4qR96PDd8p/mk37Zi9xFvvV4vrNubjZW7joR1zO+bD+Ljxds87XskJ6Q5/7LhAN6Yb58ZQe4Y7cvWzDEn80omsKA4WLHziNL/tGBjJgAYzG05+nNITfKHfZ0CkynhRJ6z5t/vv7/gislLDOvkPuvAF3/FZZP+CLscRU1AL6TPZ+3Uydrl2j3ZSBs3E2//urlIy5NfGMDcdful8hXp5QgbIhZcjLHyjLGK4jeA8wCsjlXB3Lj5Y3u779o92cjTK3VhmBplwOP+IvhDZsmWQ2j/2E/49Z/MsK4JIGptK68ggIe+XoUD2acw+JXfMOy1RWEdf/mkJXjk2zVhX/fa9/9S+v4EcsdB9DXiQbOUCccqMfz1RRj9jiLqTnEK4acql+hdcC3adBDbDp6w1NMTNh2mPUdO4pU5G5UNrFj1lt7Yr96dbd0pRqzefRTTl4U/WFrUnwSF4JIzlQx59TcAwDNFHKjx2rxNuP7DDMzfoPmLze9Briu/bzqI/dnFO7bvdCEajas2gIWMsRUA/gQwk3M+KzbFcmfueq3i3PvFcoMJAwhVYsCbxrX7yMmgvVzVRq3fl40W439wtWFn6MELizcfcr0mAHRqWCX4O9qm/JcNBzBlyQ488m2x9R08IXccxEftVwiurQdPxDRa7bsVe3DpW9GHbS/bkYUmD/6Av7bZl23VLmPgwD/7j1v2EU8hwDken7EW+46ewkndPxqOxnXF5CXo//x8FBQaK2qOpHHtysrBB4u0wcm3TVmGl+b8YznPip1Hgv5GVWBFrLngfwtx7xcrwj6uUNK4OOcGv1aejT/PHKgRS0QgjhBI5vZFXrx88hKc/7L3ITyEdyIWXJzzLZzzTvpfO875U7EsmFemL9sd/H37lGWWTNRugmv7oRPoM3EeXvtFM3epNLTPluxAfiHHnLX7LdtUgo6DI7egEHuPGs1/h0/kYfzXq4IRj3LzHa0SIspRXKYLcQ9uiHLtPJwT7GwczsmzaDFnPz8fl769GFP/ska3RcKdn/2NPx2EjVcW6Nrzwo2hqLdRb/2Olg//iEPHc8E5xy2feo/6yth2GO8t2or7v1qBU/law5ucEL6p0FyvZY2r77O/YMKMtTiSk4fjNuO9hr++CAePRx4o8/mfOwwD/IuKwsKQxjX+m9Vo+tAPAIAvM3bigE2gz0Vv/I7hRRT0lJigNZl5heqOrvm9yGZ1InbETeYMFUdNleL7lXstmajdBJcIgxUNk6z6//uLFYaoQflMQR8Xt64DgLs/X45ez8xDXkEg2AN8auY6fLpkB35ctU8/VtZGHIvpSqysb30mzgMA9Hx6LiYtMIZQy73d9xdtUx5/+ESeYVyNeP7TJDPRA1+tVIZ/A8Bj34VvrnQi2sgzUX6/ZKr6a1sW8goCSH9yDj7+Yzt2maLyFFat4LtO8Guf3IncAuTpwl/lv5m9dr8WKGRjApQ7WGMnL8H2w1ZrQGGAW0xZmw4cx2uK0Hk7Hpy+Uun7HKeH9xfl+KVtB09gasZOAJqWLkL2t2Qex/1frcS/HNwFK0xacDgEAhycc+zPPmWpP0n6+8vXtT3z8xXLlNmkaIlrwTV7nVUDMmMWXAWFAcdKJW+atmwX5m84EPTJLHDwXS34JzPYu5277gB+XK0Jp3HTVuKiN37H2j3ZQS3l7qnLsSXzuKHxMYeOR4r51vILA2GZTkSQyD49I4B2To60cTMNQRgq2332qXw8NH1V0GcC2PsY59i8u1hrjF5C982s3n0Uj367GlsPngge71dJI6gFuBBOKpKCPfZAMKuLavf7vlyBE3mFtlqFfF8LNx1UmlnzCgOW+jDwxV/x/M9W06GZQIBj4o/r8dmfOx19nyIYadGmg8GAoBO5BbbfWDgN+ihp/JMs3EVQxr4o/Uf/7D9mGXsIAE0f+gGXTfoDPZ6ei1fmGoW8/P4Aa/1+c/5mFAZ42EFhRHjEteDygrliNR//I/6jhw0D7r4lvy/0iH7ZYBVcHBwHj+fiqvf+xBvztQZ704GQj2P635op8/CJPINf4pr3/zI4w6/94C/3m9F55od1uN60/ws2k/W9NPsfXPTG71izR90DzSsIuE7HIRpJucFTNdgdJ/yMWWv2GdZ50Xjk5xXJYKGjOfm2Zq/CAMeyHVmWAaxOXPC/hfho8XY8/9OGYPkT/eqasvWgdVB1okLIiWro0ztBeQWBoI/Gp1CXhSYjtkz8cT3Gfx0axFxY6P6g8gusGpcXsk/lY92+bEMHxE7gTF+2G2njZuKKyUsw7LVFOJKTh3aP/YT/2czBtTnzhG0U6tu/bkbauJlIGzcTf+/IwhFp0O99X4b8Y5Hck4rzXlpg8IcD2rgxAPhji9YRMHewRD0QGpf5ubwydyPW7smOqMNEeCcuBJehYQuTrzJ2WXo/H/y+zbKfOaxdkKBosAoDHOt1DYlzeBr06fcxQwjzDoVpR2bn4Rz86+OlhgwJgrcXbMHc9QcMGqAqIAAANurPbqfN9bo9ORu9nplrWGf+GM2BAKr9JtiY+Ox6nvJTnS35Du3eg4oZK/Zgf/YpdHr8Z9spRQoCAYx843dc8L+FnoYcyEL88Ik85Ov3rhIudgg/iIpCvQ7kF/Kgn8R8Zjm451RBIXZl5eCtXzfjUym7hZdo2bzCwrCjagGtA3Lt+8aOkZ3mZ87IIToQ3yzfrdodY95ZjOdmbcAXGTstdWPywq3B3xe98bsychfwblYX5v9T+YV48vu1wc7AlszjtvX1ug+M40zNnbpEXT1eoicfUFXvgkDAcG99n52HGz8q0lSupx1xIbg+VAgawPjBqxp4QPsYXjFFVYXjD0r0WR+RWejYRTfJXDbpj7CyGjw5cy1mrdkXDLuds3Y/9mefMpgUb5+yzPU8FZK1rF4v6tnyl+lmw1P5hSgoDCD7VAEOHjeOOzJ/jOYxQ4J8SaCpOgOAZs7ZeTjH0gjJ70DWZsS1Oed4afY/2CZpNHKmkqXbD+OOz/7GID1q60ReofI9yEK3z8R5Bp/l0Zx8fPzH9qAAPplXiBbjfwxuL+Q8KGie/mGd54kiEyXbH+ccR0/mB8WxKE9eQQCv6maoQs7xzd+7gx00eTzVoJd/Q99nrULZiylq4IsLcPh4ZDkNzYJKlbkCMEYzauXS/m/JPKF8H4f1MW4PfLUSz/3kHLpulyDbqzYz9t0lOJVfiG+X79bbAe15X/TG77b11Ux+IceLP2/AUj1iWLzbhZsOovUjs5QZdgoC3FDGXVknDZ0zInriYj4uO/+CXH2dsjObAzZU4dh2qHqs2w+FGtMA5zh2ypuD+kC2tyiu//60Hj+t0Sp6XmEAgQDHDYoemzpDeKi8E39cj691U6XQyEa+8Tu2TRyK1o/MQo8m1ZTXNwuqFTutZrYEH8PiLe5h/0Nf1aK77hnY0rD+r20hv1uypKEIM9CurJN4Ze5G/LRmH2bdfRZ2HMrBpW8vxsVdG+CFSzsFo7WypACd/8xYg6cu6oAtmSHt09zIrdx1FN31+x7/zSp8v3Iv2tWrhK+W7sKPq4xh4YEAD77/AAdenvMP6lQu53rPsqnwjfmb8d+fNuDCTvUM5TkhDcIuKOS4W8+z2bVRFXjBa+PtNjDZK146Z4Ax2vR5hflaLvbcdQfw4OA2wWWvmpTXsgCav1LIP/GdmgNK/tl/DC1qVVCOL8zNL8Sr8zbh1XmbsG3i0GBwhuD2KX9bjtmVlaP0nRGxIy4El51/QSac0OeCgBauzjkw+h2td7sl84QytFclEA9LmREW/JPpydmtXdf5g/t2+W70blYDr/8S8i3kFQQcp2Z5ec4/OKd1reDyHCmR7lsuWQSW2IybMhdz7LtLLPv4fCxmUznImijnwL1Tl+O6vk0AIDjW6Viu1tis2XMUs1bvU86ttnR7Fj77c4chYapZMzmeG2q0RAN27FSBMslsIecGjU1+L04IU2HbR2cFNRIRzLJ+r9ag5eaH7lkWQst2eMt4Ymf6jRY7n6TX6YFkDcwpqTKgBXG8v2grhnWqh+oVkj2XMRzB9fKcjVioJ/AVdaFqaqKhw3PeSwvw+PB2uKpXmuX4I5KQ+/iP7cEMKE7cMzX88WpEeMSF4PIrzHVAdIN2Wz1sHCt96EQemktmIifkEOUFGw867Glks8uHnLEty9JL1XwhToJrI16e4z28GXAfg+UmYIHwGg83npy5zrA8/e/dGN5Fy+jv9zH8te0wRumDiQOc47M/7cd6PWjKwm6+l8xjuej25Gx0aVQ12MO3S5Qc4JFFJVbS8/0ZzGj6acQsySela3o1QcpEkpjXC7/bDJ4f8bq3TCw5kibpFkSx9+gp/GfGWiz4JxPvX9vd87iycBJbL5SyzucVBvDMD+sMQkswf0OmUnDJ7/CRb0rX4P7TmbgQXKp0L4AW/VTcBAIcxyTBZTYdRMPHimlNjp7Mtx1EGintHv3JcfsLkgYZqwGm4Y4zO3xCa8S2ZJ7AXZ+FzDEBbn8uVUNpDixZsycbB4/nGXwOds83M/sU1u0N3+TTsnZFiznqpENjG80MA7Em2vGAcnCDWwCSICsn3zbqVUU4vmKZjfuP22YJMU9rQpRu4kNw2ZgKl24vutQuduQVBpB9MtTQeTFjRsOzs9a7mvzCxU2LkB3XblqiV6b+tTOs/eWe7h6pUQkEOOYrhiUA6sjK6z80Rsepett2jvM9ETZmHNwyB9Mqh3D8k/mx016j5TWbMPZI8BoNfPhEXtAX6oW8wsgE/QaHqYrW7c0uFdOkEN6ID8Flo3H9FoaZLlZM/HG9oWFPcgh9duN/l3XBHZ9ZnbtmSnJ2VfMEepESbvb7UzaNebimO7Pg3a8QRuaxZ9GSfTI/OFzCC04D24sbLwE3scarZibILSJBvykz8mE3Xli8+RB6Nase9XmWLl1aKyEhYTKA9oiTyPAwCQBYXVBQcEN6evoB1Q7xIbhszHEl0aCbw2iTIsgzJxDzMpVmSmKKeMDeD2cOvw4Xp153rJizTvmtETEiUlOhG3Yh/7EiVmmgEhISJtepU6dNzZo1s3w+X5lTEwOBAMvMzGy7b9++yQCGqfaJC2ltp3GVBgJSpFK4VCufFOvixJxbPrWOFasRRgRYpNhpXNEkhgVKVnuNNxpVSw1rC3RNaAAAGWNJREFU/4+vd56O792rYzMPotcEz6+M6RzWee3yQsaK5MSYNbfta9asmV0WhRYA+Hw+XrNmzaPQNEr1PsVYnogpzYJrpj7257azm3s+pmpqIj6+vjtSwpiHKVqu6Z2G1tJMzd1txnB5wc3cMbxzvYjPLcgtRQELRcn957cq6SLY0rJ2BU/79W5WHe9e3Q1NapS33ad6+SS0q1fZsC7Sz9qrqbBB1dCYu4rJ7salGSv2RFYgj0QyC4ANvrIqtAT6/dnKp7gQXH6FqbC4ZNldA1p42i81yZvV9eGhbTDr7rNwZouatvs8dZFtRwMAMLh9HU/XkpkwrB36tQpdM1KzRa2Kya7HpsTgA10TxgBOeV6z4qRLoyq4smfjqM5R1ME90VDBQ2MPAH2a18CANrVRMdlqdRBaW7kkv8XC0LyWN8Fo5luPAqZyudD1aldOcd3/y6XhT3QZDimx07jinlGjRqVVq1atU4sWLdpFcnxcPElV0tJwBiyGQ9/mNYK/bzyzCW46q6mn4/o0t9dCLuveKPh7WOd6qF3J/iP65b7+GNy+ruO1zHLDTSOddJVmornt7ObBcnIODO1YN2xzpd/HXDMcRBOwIvgjjCABL0L4ZtN7HNIhfOFvJhDgeGJEe0y7pXdwXbiNk90YxdJAuB2CCilWQXdhJ60ul0v0IynBh3n/7hfcJuqliik39LDd5jVasZJUnvJJfqx7fBBW/+d8wz41Knir/5/f1BPntq3taV87YqhxxT3XXXfdwe+++y68AagSpferkZBTPoleWvUwG9yRXetb1sl+KbFdTJ9wfrvaGD+0Lcp77HXWqmgvjMb2DAkup3FfqyachyY1yrvOiHvr2c0My5XL2fvX5tx7VvCDq5SSiHvP1VIvBTjH65d3xbJHzjXsf5vp3GZ8jLkmwjWn6CrvYYZfs7ALJ3rQS+fivHbGRufW/mrT7uIHz/F83SdHdABgHPtUQaF1OBGuxtWzacjEO+WGHphs0/i/c2W6p/N9dJ3ml1LVy66NquKPBwfYHntmC62TF5zZ2scwcWQHwz4iR6Wo001rhrSsOpVTMLpbw+Dy0I6hDlssou8qpoTexS39m6Fckt+iRX59ax90Tws9U9mcLtOzaXVP9diJGPq4SpwNGzYkNWnSpN3o0aMbt2jRot2wYcOafPPNNxW7du3aunHjxu1/+eWX1Hvvvbfeo48+GvzwWrRo0W7Dhg1JADB48ODjNWvWjNipGBdRhXLSUiFsvJox2tWrhMu6N1KO+fru9r44U88q/vDQtqiYnICujatGFJ5sl09xZNf6KCf5shJNDcTzozoFp2wQ25IVGsvwzvXw7XLNRNKxQagnPPmqbnj0W/WI/iS/D81rmT9ErZx2cuG+81rhm7/32IavJ/jdNS4zshDq3LAKlu80pjV69uIO+G3jwYinkL+gYz0s2nQIn/25Axd2qofz29W25JAz93brmExHY3s2wsKNB1HXlIvwyRHt8bCeMaFicgIevqAN/m/aKrSuUxEdGmg+GzlzfIVkPw4e1wJvzmldC7f2b4ZzXvhVWe7q5ZNwUZf6eNRhviuZBlXL4fObeiFt3EwAWuNunsTS7v5UTLmhR3BgdEqiz5KhpXK5RNSpnILK5RKVQS3CDCjXhzHdG2H8N6tRGOAYP6RNMD9gOUWjn+jz4cpejYOTRb5+eVfMXKndmypvoIpf7++Pfv+dH1xe8eh56PT4z7imd5rhmoNsrBgNq6ViZNf6+HPbYVzUpT5eGt05+HzNqDLkd21UxXOarqLQuO7/akXDf/YdCy+KxoWWdSrm/PeSTq4DL3fu3JkyderULenp6ds7duzY5tNPP62ekZGxfsqUKVWeeuqpuh07dgxvDEwYxEUXQBYKI7s2AAC0qVtJua85MGDmnWdibM/GhlQ0goZS1FS18kn4z/D2wZ6n3ZQKALDpqcGWdYl+ZpsgVTahmDWLelIDIwSX6qN98VJjhFTGwwPx5/gBGNi2Nvw2vXbRsMqIR2lnXmOMBT/4J0dYfW1+5iy4GlVLtWRfKAxwTLiwLZrUKI+v/tXLckyllETHhupi/Z0DwNW91D6lJjVC7/KCjtbgELkRW/7ouaiWqmnsNSok4blLOuLJER0w//6zDcdseXoIxko+rM9v7onRZzTCN7f1wWc39gyubyvVRaGhV0pJwPOjOhk0jGBZ9I7MR9d3N2gFr1/eFYPa2ZswzQNkGWPKDtMbV3Q1dG7saFG7YjDdlCpQSAi/5Y+ei4bVyqFWRaN5XpTdXB1EkQa1rxNM8ntIkaXe52OWDqhKe25YrRzutPE1166UYrA4VE5NxNZnhuCxC9sq9weA1y7vAgAYpic+FmWspDB1ymw/ZB1v5lRvt00cajAvqjqk8Uz9+vVzu3fvftLv96Nly5YnzznnnGyfz4euXbvm7Nq1q0hDj6PWuBhjfgAZAHZzzi+IvkhWhDmlf6uaGHNGQ9SplIJ+LWsqUyR1b1INWzJPWDIVeB3/k56mja26pk+a7T6qcWWMMXxyQw+0NaVTYmDB3HWA1R8lf/RyIzT91t4Y+cbvAICRXepbGig5JF3u8T8xvJ3jjLXC1Hr7OfZBJyJ7fr0q1l67z+dsKpx6c09M/m2rYV1BgOOaPk1wTR8tce60W3ohNSkB936xAuv2ZrumGXpgUCus3n0UG/Yfw5jujfDhYut7r6ILIrtTpST6USU1EUdy8oP7bps4VLnvvH/3w/ZDOYZZdwEEo+I6m3w/SQk+bHl6CHILAhijT0ki53KcdFU3pCb5ccVkLVmxuN/ypoCequUT8eplXTD23SWGGY3rVymH3UdOKmcqUAmuIR3qBu9j8ZZDGP+1WiP3+1hwTjCV4BLrGGP47QHNhDpu2kp8rmdBEaYzc2dM6/Rx+H0sOA6wW5p6zKLZL/bQkDZ4aEgbw7r0RlUxuH2d4DQwglv6N0NKoj/4Tb0wqlOwvE5c0LEezm9XJ1jPs3Vt0snkDgB1K6dYLBHmx18hOSE4KzQAjOhcP5iZpSgElxfNqKhISkoKVkifz4eUlBQOAH6/H4WFhSwhIYEHpFyhubm5MYtEioWp8C4A6wCoVaAYcDxXEzoVkhPAGMPZUjZ0QXrjqli6PcvWBLbPlDHhl/v6K/erVTHF0qCN7FI/OJOxE6rIQg5uaBS8mkC6Nqoa/AiET2vSVd2UgQ/yGd0c6hVTEm0bbIFosFWTJ7apW8kQqj6sUz18J0V5Jfh8rsmP0xtrPoUGVcvpgsuo31ZMSTAMBmUINfZyPsKXR3dGTV0LcJuqpkJSAubf199TB6ZpzQpKTckJn0/TVMVYIDldlNmpL4RaarJRWNSrXA5JCT6DZvPkiPYY0KYWej0zLzg9x5x7+wUbRzsTtbgPYaZtVrM8cgsCBtOinzG0qq2Zkge1r4N3FmiTQs68sy92KLQLALiub5Og4LrhzKY4nleAa3qnGfYRryLBx5CvX79Hk5DPauLIDpi2TIvg82Ly9/t8ys6N0MTFM+gZhl9MNtkLbbxBVbXFTQRwvH1lOtKfnGPYZv6eF/3fOcgtLESCHnQjvx6v335ZIS0tLfeHH36oAgALFy5M3b17d8y0sKi6AIyxBgCGApgcm+KoEVN42zlOAaCTZBpR1Y8Jw9oZAjTSqns3C784ujMW/t/Z7jtGgGiHeys+OpGjUXxk57atjX4trWH04oO4pX8zQ0MQach7Od2J7PcxfHDtGQYH+rMXdwjqW4zBEryS6GeW528XQCDKpwV8aLwypnOwkTuvbW30aFIN1conBe/Lxxjm39cf393eByO61EcfPQq0iu77rGvj20lN9qNKahLqVXGfT0tFU4cxSjJuWfOHdqwbFCZmjauxXifFs7hnYEuM7dk4+P6F0G5eq0JQ63MT2MIUmOj34cPrugcDKgDA72doVaciVv/nfIzoHPo22tWrjMEd1D6hlrVD32C5JD8eHNzGoq2JDo/fx4JJmmVBMaZ7I3z5Ly0S04sW8tCQ1soOhxBYQuNS1fffHjgbC+53/nav69MEz4zsgEvSGxjWJyX48MSI9nhrrBbooopkNvcbKqcmolbFlGC0rvg2z4syIjEeueqqq7KysrL8rVu3bvvaa6/VbNy4cbA3d+GFFzbp27dv661btybXrl2740svvVTD6VxmotW4XgbwAAB7iRIDLu/RCMdOFeCGM+2jx67s1Rgn8wtxbpvamKNImtqneQ30aV4D05dpmlO4vZ9IQ7ydfGVA6ONTaWviYzcHdFivoVGjQrLrvl4QZsjjpwowuENd9G9VC18s3QnOtXKK9uGtsenBGZoFZvPaUxe1x0Cbj1Zox+aPPylBWzGgTS2MPkOLyPzf5V3wxV+70LpOReW7O6d1LTw/qhMu6KhucKN5LnP/3c9zthCnbPprHz8fSX4fZupBKOVMDb64r6N6IuAaFbXGTwgn1azH5udtRryrlEQ/mtWsgI+v7xEMPhDnrZCc4Ki5mfnyX70cJ0/1BTUuX3CWbLtE2V6+w+oVklE1NQm3n90cny7ZHkyULATVdX2b4MmZ61A11Rpp3NBD9o+kBJ9hyEqwbIBlnN43t/XBpAVbMHPVXjSunur6fYtnUdby97Zq1Spv48aNQZ/EtGnTtqm2LVq0SBnyPmPGjK2q9V6JWHAxxi4AcIBzvpQx1t9hv5sA3AQAjRpZK4cXUpMScM+5LR33aVKjPJ7RQ3FfGt0ZXZ+Y7enc7197hmXqCxWVUhJRrXwSLuve0LQ+wWYmYiOz7j5TaRbp0aQa7hrQAlcqgg7E+DW3eY3ydTtyapI/JoJLmN/k9ErLHzkveB2hEzBYI6UYEJwluG/zGrhc0SAIxH2Z2y5xDwnSGKe6lcvhroH2fjnGmKXHHCuahWE2FKaxl0db0w2JzsmUG3rg+1V7g0LnrbHpBk3xyEnNwiCEZeVyiTizRQ3cfJZ1qIKbwGlbtxJuOqsprlLUL3kIWTjV5ow056wrQRMzC82H5hT236NJNZylsCQYzuljuO/8Vpizbn9QcIkoyBvObOrYqY0UlUzt3LAKXr+iK16HJjjlSVYv76EQfszbN0yERzQaVx8AwxhjQwCkAKjEGPuEcz5W3olz/g6AdwCgW7duRfL2zDnVnAbV/vbA2Yb5j85uZfWXqUhJ9CNj/EBLD3fxgwMck3OKyt+6jtoF6PMxW6F898CWeGDaStdBwieDUVGJtj3bcLjn3JbIysnDMMl8VFka8yZ8AVVSk/Dv81qiQnICLuxUD3PW7UfFlERc0aMRth08gZvOaurYoxbfsnkfIbi8TGhZ2hAal5OG3rt5DfSWBroPMmVC6dGkOlbvzg4GC/h8DB9frx6Q62Yq9PmYJdhBdazKnxkpXfQhJYl+FuwUJjgMtJ56szXS1I7q0oBht2CKaHHTphhjhud2z0DrdxzSuEhwxZKIu+ec8wc55w0452kAxgCYZxZaRY3It/fDXWd6PqZhtVS0qB2ZZVNllimfnKAcM1Nf96X0bBr5QMpLz2iIbROHuqaTEvb/SuUSDFGLXj+VVy/rYliuUSEZb1yRbtswjBvcGm+NTUf3JtVQMSUR953fCq3qVAzma0xJ9OOJEe1dzTRBjcu0XjT6eR40YTueuqi950G4sUQ01NFMMPrg4NZ4+8p09PCQT1I8K1Wj6YasrYVjKnTjjSu64tvb+iA1KcHgY4sFLfRxiU+MaG8ZcxdrvCQZfvSCtmhZuwKevqhD0FIh4wv632JevNOauBiAbMe7V3fDtoM5ngcjFyeLxp2DA8dOOWbUiBU5ksbl5vNQMaxTPdzpYV4wQUqi36IlRIO5tx/UNKJoS6/ooZnG5v67X1AjLQ5E4EU0aa8S/D6c7zCeS8bvY8Eo0Zfm/OOytxFWRBpXheSEYHRrSHCFf/5f7+9vEXjjBrdGi9oVcNkZkbkdvPDZjT3x9oLNeO7ijq77tqhdET/f0892u49MhUVCTFp8zvl8APNjca5wqJiSqBxkC2hRg8Ux/YaZ6bf2xk59YrziEFqAFm25ft8xVC6XiFQ59L5Yrh45so/rwcGtUVAYwHlt6+DctrWR6GMYld7Q5QzuhOOfigXCvBlvg01jqXHJCEFuN6eeE42rWyM5UxL9wU5JUdGrWfWYpJwCEBzonqa4FyJySp+qEiPMWRCKi66NqqJro+KdIHLKjT1DUU6MYdotvXDxm4s9mwoBLfVUm7pFGhxqQXRCfYyhXpVyeHNsyLR373mld7oPJ0QUXSwSDYfLK2M6O04tInhgUCu88ctmw7qimjooGFVYiqcmKko6NKiM9689A72icBkQVuKrW0goqVY+CVf2bBzVAMdL0htY5ksqasQYrEjHVpVmSkJwDe9c31Oqp1v7N7dkSY/ExOwF1Tiu042zW9Uq1rn3SjubNm1K7NGjR8umTZu2a968ebsnnnjCW4ScRJnVuIjSzy39muGiLvXLpOCKO1NhEWV1eP2Krpi0YEvEc28RZY/ExES88MILu/r27ZuTlZXl69KlS9shQ4Zkp6enn3I/WiO+vq5SxLRbemHKjfZzBpUkwrcWTURjceDzsTIptABrRpHSTlFpXC1rV8R/R3UqMh8aUTJEM61J48aN8/v27ZsDAFWrVg00a9bs5I4dO8Kapyq+vq5ShMi3VxppWC0V8+/r7ylrAFE01CyBwKBoKM0zMRMOfHNbQxxYG9sPvVbbHIx4vVimNdmwYUPS2rVrU/v16+dtdlAd0rjKKGk1ylMvtwSJJIquJHEbK0gQZqKd1uTo0aO+kSNHNps4ceLOatWqhZVtgGorQcSQV8Z0xraD6szqBBFzPGhGRUU005rk5uayoUOHNhs1atThq6++2ttMnBIkuAgihgyX0mTFG2+NTXedG40gvGI3rUkgEMCYMWMat2zZ8tSECROsGdE9EF/2DIIgioxB7et4zthBEG7YTWsye/bsCt988031hQsXVmzdunXb1q1bt506dWpYY3FYpHM2RUK3bt14RkZGsV2PIAiiLMAYW8o57wYAK1as2NapU6eDJV2mombFihU1OnXqlKbaRhoXQRAEEVeQ4CIIgiDiChJcBEEQRFxBgosgCCK+CAQCgTId/6nfn+3YLhJcBEEQ8cXqzMzMymVVeAUCAZaZmVkZwGq7fWgcF0EQRBxRUFBww759+ybv27evPcqm8hEAsLqgoOAGux1IcBEEQcQR6enpBwAMK+lylCTFOo6LMZYJYHuEh9cAUObHLpigez49oHs+PYjmnhtzzmvGsjDxTLEKrmhgjGWIAXinC3TPpwd0z6cHp+M9FxVl0T5KEARBlGFIcBEEQRBxRTwJrndKugAlAN3z6QHd8+nB6XjPRULc+LgIgiAIAogvjYsgCIIgSHARBEEQ8UVcCC7G2CDG2AbG2CbG2LiSLk8sYIw1ZIz9whhbxxhbwxi7S19fjTE2mzG2Uf9fVTrmQf0ZbGCMnV9ypY8OxpifMfY3Y+x7fblM3zNjrApj7CvG2Hr9ffc6De75Hr1er2aMfcYYSylr98wYe48xdoAxtlpaF/Y9MsbSGWOr9G2vMkbzULvCOS/VfwD8ADYDaAogCcAKAG1LulwxuK+6ALrqvysC+AdAWwDPARinrx8H4Fn9d1v93pMBNNGfib+k7yPCe78XwBQA3+vLZfqeAXwI4Ab9dxKAKmX5ngHUB7AVQDl9+QsA15S1ewZwFoCuAFZL68K+RwB/AugFgAH4EcDgkr630v4XDxpXdwCbOOdbOOd5AD4HMLyEyxQ1nPO9nPNl+u9jANZB++CHQ2vooP8fof8eDuBzznku53wrgE3Qnk1cwRhrAGAogMnS6jJ7z4yxStAauHcBgHOexzk/gjJ8zzoJAMoxxhIApALYgzJ2z5zzBQAOm1aHdY+MsboAKnHOF3NNin0kHUPYEA+Cqz6AndLyLn1dmYExlgagC4AlAGpzzvcCmnADUEvfraw8h5cBPADjlAVl+Z6bAsgE8L5uHp3MGCuPMnzPnPPdAJ4HsAPAXgBHOec/owzfs0S491hf/21eTzgQD4JLZe8tMzH8jLEKAKYBuJtznu20q2JdXD0HxtgFAA5wzpd6PUSxLq7uGZrm0RXAm5zzLgBOQDMh2RH396z7dYZDM4nVA1CeMTbW6RDFuri6Zw/Y3ePpcO8xJx4E1y4ADaXlBtDMDnEPYywRmtD6lHM+XV+9XzcfQP9/QF9fFp5DHwDDGGPboJl8z2GMfYKyfc+7AOzinC/Rl7+CJsjK8j0PBLCVc57JOc8HMB1Ab5TtexaEe4+79N/m9YQD8SC4/gLQgjHWhDGWBGAMgO9KuExRo0cOvQtgHef8RWnTdwCu1n9fDeBbaf0YxlgyY6wJgBbQnLpxA+f8Qc55A855GrT3OI9zPhZl+573AdjJGGulrxoAYC3K8D1DMxH2ZIyl6vV8ADQfblm+Z0FY96ibE48xxnrqz+oq6RjCjpKODvHyB2AItKi7zQDGl3R5YnRPfaGZBFYCWK7/DQFQHcBcABv1/9WkY8brz2AD4jzyCEB/hKIKy/Q9A+gMIEN/19/g/9u7dxsEghiKorcFeqEjKqEHOiCgJRIq2YRgVmSbkBnOkVyAJ3majzx1+oOer9Wz9YvtvfWa7qd6rh6tO7yttXO6fNNjdd7X6VXd2icaqeMy8gmAUSYcFQLAh+ACYBTBBcAogguAUQQXAKMILgBGEVwAjPIGCpepPT6nwpoAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(trace.mu1, label='mu1')\n",
"plt.plot(trace.mu2, label='mu2')\n",
"plt.axhline(mu_1_true, color='C0', linestyle='dashed')\n",
"plt.axhline(mu_2_true, color='C1', linestyle='dashed')\n",
"plt.legend(loc=(1.01, .01));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Maybe that was just bad luck. Try it again with a different random seed:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [t, mu2, mu1]\n"
]
},
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
" </style>\n",
" <progress value='3000' class='' max='3000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 100.00% [3000/3000 00:05<00:00 Sampling 2 chains, 0 divergences]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sampling 2 chains for 1_000 tune and 500 draw iterations (2_000 + 1_000 draws total) took 5 seconds.\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAD4CAYAAAC0VQLEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dd3wUxfvHP3OXhJBAqKGXUEIvApGiqBRRithRVBR7735V1K+KvZcvYvkhChbEhqCCIgiigAgEFOm9hZZQA4SUu5vfH7tzt2V2b68kl0ue9+t1r7vbOrs7O888ZZ5hnHMQBEEQRLzginUBCIIgCCIUSHARBEEQcQUJLoIgCCKuIMFFEARBxBUkuAiCIIi4IqEsT1a3bl2ekZFRlqckCIKIe1asWHGQc54e63KUF8pUcGVkZCA7O7ssT0kQBBH3MMZ2xroM5QkyFRIEQRBxBQkugiAIIq4gwUUQBEHEFSS4CIIgiLiCBBdBEAQRV5DgIgiCIOIKElwEQRBEXEGCiyAqE5wDq74Eik7EuiQEETYkuAiiMrFjETD9NmDuk7EuCUGEDQkugqhMnDigfJ86EttyEEQEkOAiiMpESYHynZga23IQRASQ4CKIykSxKriSUmJbDoKIABJcBFGZKFaDMhJJcDnm69HAO1mxLgWhoUyzwxMEEWOEqTCJTIWOWTcj1iUgDJDGRRCVCW+x8u1Oim05CCICSHARBEEQcQUJLqJysPdv4OiuWJeiHMBiXQCCiJjKI7i8HsoWUJmZ0A94u3OsS1GO4GV4Kg6s/AzwFJfdOYkKTeURXNNvA15qHOtSlF/+/hzYtyrWpSBKGxYDjWvtd8APdwN/vFr254413hJg+h3AkR2xLkmFIqjgYox9zBjLZYyt0SwbwRhbyxjzMcbiI050zbexLkHo5KwoOy3x+7uA/zu7bM4VDF6G2kB54tsbgS+vKZtzleU9Flk6Th4su3OWF7b/Aaz6AvjxvliXpELhROOaDGCwYdkaAJcC+CPaBSJUio4DEwcA34yOdUnKlgPrgGdqApt+iXVJyp4104ANM0v5JGWocZ06Asx/AfD51FNXRv+a6CBUxmsvPYIKLs75HwAOG5at55xvLLVSlSbx0pv3eZTvnOWxLUdZs3Ox8l0ZBVeZUgbvwc+PKubBTbPVBRWo8T55MLR8j5VSaJcepe7jYozdyhjLZoxl5+Xllfbp5Oz7N/Cb+0r3XNsWAMs+jPw4TH00Pm/kx4onPIXKd0JybMtRUSnLBrQwX/nmFVDjeq0V8EpG8O3ipJ8cb5S64OKcT+CcZ3HOs9LT00v7dHJO5moKVMqC69OLgJ/+42xbTzHwYmPg32/M66bfrnxXWsFVgQfIeooUp30sKQvLg0+9RpdI0FOBBJdjyFRYGlSOqEKfRliVtuAKhVNHlNxxvzxuXrfxJ+WblzPBlbteLmijhadI+a7IGtfz9YD3z4zRydUGtCwEl9cguCqSxiXj+H7FRyujol97GVM5BJe28d+1BJhyRaCBjCWiMtsJU+HrKi+81xv47ubSO754LgteAha8UnrniTUHY+QidlLnooWouy535Mda9iEwtkYgu3155O0uwPt99MvixaceZzgJh58KYAmAtoyxHMbYTYyxSxhjOQD6AJjFGCvfnnStue3Ti4DNvwQm1Islwo9lZwiPdgOz+Vdg3feGc2jOv2NxdM8XKtoOxYIXY1eOio7TevV2F2DmA+Gdw6hxRWIuWzxO+daa/csbXllnmEyFpYGTqMKrOOcNOeeJnPMmnPOPOOfT1d9VOOf1Oefnl0Vhw0ZmbgvmO/KVQY9UCAxjI1Ka555yGfD1dfJyAMDRnWoZvMAnFyrjUMoS4eMiSgerOmfF0Z1A9sfhncvv41I1rkjMZW5V+EXD5+v1KNfkjYI14+BmYLeDyF8yFUaVymEqlJnb7F7cXX8Bz9YCdv5ZemXSlsFoTihr86BMsJ/MA7b/Dky7GdizUrHfW/HLE85eXid4KS1QyPh8wNyngI+HKOY0O/x1rgw6ZkIw+M8VQeMttLZoBLUs/1DRIv/8H/BW58je8/FZwEfnWq8nU2GpUEkEl+Qlteu5bf1N+d72e/jnPLYn+Db+F1pTuU8dLfvGW3cvDI0L58CH/YEPzrLed8l4+5fXKf98ARzLifw48cb0O4DF/wt//5zlyv67NA3woa3Aqy2Bec8ZNg5R44oEoXGJ+hWJ1uFK1B8zEk6qw3J2/gkc2wX8OjbyYwaFNK5oUjkEl0yj4F7g6G5g5oPm5J9Ox514PYrtvfiked1bHZyXi3Pl8831wCvNlcAEGeu+BzbOlq+LBG0j5r9m4cRXy2j0LYieZLSEbHEBMOMOYMfC6Bwvnlj1haIxhYvRvLp/DfBOd6DgELDwdf26cDUub0no2o7YPhoal99UGKI1oui4uZMqjnFoa+TlCgppXKVB5RBcMu2K+4DPLgGyPwJy1xpXKl/MBeRuUMwvW341H2PVVGDuk9a95em325sKtKZC7gPWTlf+r5km3/7r64CpV1ofL1y4jcZVeMxiH7XsIpiCRViVIhGAXg9weFtk54+EZR/anz/7Y+CVFqXnuzTeO7uyhOrjEryeCbzczPn2p44EtKNoDEB2heHjKjkFvNREMWUDylAOrSA7sl1erjlPyt93p+iG34i2hDSuaFJxBdficcCKycpvq+CMQ5uV31YaFxiwZ4Xyc7UkSa8IZLBqtFdNtW+Qtb1fbU+yrAenSjUuFVEu44y54uUXZWUu4OQhc1Jgnw+Y/ZjixBb/Cw7DRDjDE357UTEvzn8WGNctNvNteYqUAecfG9N5apj5AHDqcGhmriM7FWuAz6s0fnY+RmMds9NKtIJr9bfOAy9OHQFKHIai525QskqI5xENs6QwFYbybpScUr5XTVW+3+sNfH6Z5BiGOv/nOGU7I39PkVtXjEifMwmuaFJxBdfcJwMZmaUal2ZZiaEyantJiepAWFm0m7CVp9SxLod4eWSIchkFl7Hie0sC2phTvCVA/l7gj9eCv+xSH5dBU6ySpv9fdFzRdESjyVzAay2B98/Qb3d4G/DXe8DUq5T/v70AvNoC2DxX0WTFgE2raEK/OUfC768o5kUR3u8k+3g0IskWvgl8OEDJxSeem/bcBYflwjmURve7WxVrQM5yYPlE4I22wAGjZUAcNwTB5fdxcWDaTfah7uEGFuStN5QnCj4ud4Q+LnEtu5eGf4zv71Q6YcHQPWcyFZYGCcE3KUdsnA3UbqkIjIwQMg/INC6tQDEOatSaNkQGhxJJwyo0NaM2otvGJrzb3zBwvfAwNq5/jgPmPWve/0QeMHkYcNVUoE4r/boZdwKrv1Z+p9QFsm6wKYekR2xc5jJUlddaAi3O0TSmaqMktFABM/jKxPQyq75UvtfNAOp3sNa4juwwX5tsG8CZudJzCnBXD76dHfOeUb73rACaq/VQW8debaF8jzWYWUNpMMXxPtaMNDm8Dajf0byt0WJgW+dsfFzeEuDnR4C+DwI1m0YvutV/byIxFaoh9XbCf9+/QPUGQLV6kjJoBIjxuqwEKufmdY46R5rnIc676efg+xGOiR+N60Su4t9593Rg8tDQ9pX5Fk4dDfzWmkC8JcDit5XfzBUQXB6J5uR/mWzMgUJAFh4z92C1wRk6jcvwYp2wGHS5boaSgWHJePM6IbSA4GY4J4JLxvbfgQL1RQ4mNA5vA15sEhAyYgD4kvfUMlo0trJGpaRQHvbNmKLpHN1tXY7iguiGKGsF9faF9uHoRs0/1HIwl9I4G83Wxvon62QZzyl7vlvnK6bDnx9RjxuhdiMIVePyec3n9ofD27xr/3eWYg70l0MTtau9XpMFhgGvtgKW/p9+8dZ55nMwFvy5WQl8u+dChER8CK4l7wKfX2q9vjAfWDrBukLJNC7tlATFGr/MjkWB38wVeGFkjb+Tl8lTCOxfrTi2sz8ylEt9mbxFemewcQR+QhXr44tyajGalIIJFe2LbJUS6ISNj8XuHNrlxccDv0X0oFhmJVxly09JzHDiXG91BN7uZF3ON9oAC17WLzuwVpla3ue1NiWumKwIJeNUFlrNfaOkV601GRob4zn/tS6nrC4zlzI0YdpN+mMZ60u4Gpd4D3xe4M/xFpkgNBQdt2iMjR00h1GFXo9yXZ9dDDxXV79O+LgWvQXsXGJ9jIJDgd9+Uzz01ysTyAUHAwJbIPNnudzA0g+szw8YhnRo7kV5zvoRZ8SH4Dq8XWn8tWhf7FkPAT8/rBc6WmQ+Lp3gKlBCiE1OcBYQeruWmCuyUeMyNoiAMp7rg77Kb2PDpn2ZvrtFXnYAcBsEl7cEeC4dWPmJ8n/5RCB/n/K75JTZz7TtN2DiIOtG2anGJWuYBcbe9MwHFV+QE/NdwWG5RgsAU0cG31+LLIBg63z9/0Vv6v9/0FeZWn7caQEznxER6GOM2CsOMkO1yPIPmHvixnLpkAkud+AYItgFCM1UaDeOS5jMN/8CzHlCHpCk5aUmwAv1ga+uNZzCQnBp68jPY4B3e+m3ez0TeKOdRbYWjY9qkhoIc3ibEpSTv1cu6LXXqO28Gp+D9r92MLLRPA4oz+CfKebl2g7Wh/2BvI3Ku68V7CdiNK1TBSQ+BJcsSWfh0cB8P8Jc9ckF8v1lqnuhxlS4/kfggzMVM8nxfYHlzKWv/MYcf35trBjYtVQ+/mqKJDrJXy6Hob3GKT4K8xVhqRXm639QvmW9xI0/ATnL9NemxWlKLBFhKcPYgGd/pPiCivKt9xG82iK0qEIrDVdb5tdaB4TG1Kvt9xfP+Ogu6/ImpirfxvtbGOT6tDkxt/+ub2BlGfBFnbLSuATFJ5QIzt3LzEMW7AKCZBrXwS3qMQ3XVnQcjhB1b+EbwBaJec1/DzSCa+n7QN4GxbwqOHU48C4LjuxU/LubJOMXl7ynCK9131tkx9HUB52p0LCttj5MGhL4LRVcLkV4GTFOfPpuT2UspzYhNWlcUSM+gjNkFei1TMXZPdbgOzq6W3Esa5E1zNqXdPdfyveWX4EGnQPLGdP7x5JrysvnLbY2X+kwaCW7lzrYB2aNS9bLF/fATgB8bJFSUtvg+30gkobT6Vip11oHfgttMxjBchQe2gqkNQISq1pHGmqf88k8JQz64veDjxFjbv2+nmJzZyFJFVy5hog5u8b97yl6YfP9XYHjAEBiinkff6CPTHBp6o/PC8x+FPj7c6BRd/12TgKCtA35+B5AUnWgs6GTZeWrKcwHktPMy0UA0WUGk/ihLYHfRj/SF1cAT1h0qADgnR7WQS2i81m1llxY+6wEl8WAZCMyAeVyy311X19rXmakjc2QCSIk4kPjkpmbtJVZWyl/uFvf6O79Rx6cYaWZ/K6ZSoO59A2a0dfkH8sU5hQpTiecNPbMpYJLDAi2abTyNWmoxD0pOKyYyPzH0YToG3E6KPNkGCaRPyUBJoJdS5VMEN+okZFW/s4PB5iXnTgg77gAgXtgvL9znzRvKwSO0Q9ip1F+f6e57mr9H4kSjUuEfQfTuLg3kMlfKxgAZWyblp8fDVyrlY+r+Diw8lP9st9eMJdh7Qzg5abA3r/N6/xlswle+PkR/T0U1yQbPsC5fSSmMPcnVpXXe0uNy3BMK8HlkrQ7kQy0p0HIUSM+BFfQ+Xw0Lwrn+h72hHMswuGdDKZk+t7Zl6rJ6a8PgJwVgeOW9oBht0HjNA7yBQIvptNBouJl3WVwdPtsBFdpzoUktF4ZH5+nfG/6OfTs4G+0lS8fW0NJpLzrL7NGL2uUkyTaEWDQuIIIG0B/ri2/KnVJtz5RibyUZSwxalxiG6Pw1JrBASWYwJ8dxkajdhJJunmu8j2hn/U2uRaTKcreQ3F/Xs80rwumKR9RIzo9xRaCSzPcxNZUaPH+yiw9e1bYC22iTIgPwSVT2bVofT3cZ67Esoq5Wp3Fd7DNZIUnDgB/vR/4L447+1Fg4oBAI+opkjtso4WxsS6WmKfmPAFMGuY85Fa8vMZ7tUyNzpQ1YtFIcBop0c7Yv2Q80NZgwjGacFd/q5jkZATzA+Us0/83NobzDUlwuQ/4X9dAOiLdOm0HzWudjsuOcFM+Cay0E23ZjMEv/m1skl3LjhvM7ymGIngK5fVeHLso3zBO0qhxWdTrT4abzeMHN9mXiSgT4kNwyXo+gpJCc4iyMcrKLq1Nq/7W6xa9Cey0iFQElMghQIm2W/+j9XaRYnzhZeloAKWsTjUuEZFofOH3/aNoYWWRPTwcrAJwwsVKSHOuRGGeyAWm32a9v5PgEy1GDSzZMO7LzleqNQn6vNYmUDsiFVxW59yzMvi+MpO97VCSIIJLCLt9/wB7NecfW0PxpWnL+sfr5v38ZbAZaD2um30ZiJjgZAbkjxljuYyxNZpltRljcxljm9XvWqVbShuN68939P+5z9wTtWsMtM5yJ2hfPmHeCjbGSXBgLbDhJ6XRyZ7k/JxG34MddhFlWmaPUZKPfn+neV3ehvAaxXhkw0x5b19kkXg905zuSotW43Iyq3YwwWWH1idqjGILhhBYO9RQ82hrXBMl/kVTGWTRqyXWA+yDTiqqmk6XTzR3LlZM1mtZWv+s0TxOc8DFHU40rskAjOEwYwDM45xnApin/i897ByixvE/3Accs8mcYCRUwaUNlw2V43uBL69StJ2Z9zvf78Ca4NsIQskeL8u4ASgabHnVuEoDWR7IDwcEBozLIugEWm3fKqu/FuNYOLt0YXYs/zD0fU4ejDzxbST+XKtzTrvJ4lw2AmXv3/bBDsxtOJ9NwEi0TOBXfRmd4xBBCRoOzzn/gzGWYVh8EYB+6u9PACwA8GjQsx3crPhhtHS8GOh5i+L4nzLCvM9pV9trXHkGm/OuJcrIe6fk2IxNkmEXROCU3yzm2yovZE8yj1mrbBzQ+E1F1oZosNmgKZWVZvvDPXqT+7YF4R1nw8zwy7DiE/nyAxbBHHZm788usY9ezN+rv7d2E5QaXQ3hIssnqkXb9jlp97pdo8y48PV10SlfBSJcH1d9zvk+AFC/JVktFRhjtzLGshlj2SUlYfZs7HxcQc0JwY4dLGIxTOwCSqKVvLS04D7bDmqlozSjKctqanfOoXuo0WqsQ8Fq2Ihx0LFA5LWUUXzSXugXHonOTAAhQeHuZUWpD0DmnE8AMAEAsrKyOG6YJd8wKQWW60QiVhnV60dWQLvgjEgY8or1OC1Hg5UdcNNc4KNBoe3TZrA8C4GWdsOAziNCP3ZFwJVg7liUBEnrFAllZZLdvwq4cgrw1TXqecupD9OVGDDd2fkMg/mluA/46NzA/+oN5ZGa0eS85+xzqsraN7t2DwBS6yjrbyShqCVcjesAY6whAKjfpZvLxE4rsgpTdsKD64NvEy7h+i5CoWnP0PeRpRkyUnQiNj3y8kCNpuZlwdI6RUKo49IiIW9D2Z0rXJpkBX5HM5t6aQgtMaWNoLSsN4SJcAXXDwBGq79HAyhdh0ik08Jbkdao9M5Z2pV45BfBt5EhSzNkZO10JRWPkUFBbPgVgYZdJAujYM6rb5Gx/uDGyI/tFKdDJWJJzWaB37LMHeWJdgZ/fWm1U4QJJ+HwUwEsAdCWMZbDGLsJwMsABjHGNgMYpP4vxVLGIKWi0UeVdWOoB4haUaQYXxqnBJsiBbDO1F69YXjnjAV9Hwxvv7YhzvXmlK4hZrkvDcLx1TXtFXybaJLWOPC7vAtaYxtBgqvMCHqnOedXcc4bcs4TOedNOOcfcc4Pcc4Hcs4z1e8oOW0siEWFMIbIRmL6qxahH87IoOeCb2OFE43LimgKrprNgfbDo3c8I3VaB9/GyC3zQ/c5NTsj+DYA0Ofu0MsTDnZBQSWS/JyCqrXly8vSlAkANRoH36a8YLSqBMvwQ0SN+OgiaB3JZ9xTuudKrgGMlmTBcKKplDZdr1K+IylLSgRjxVPrBt9m7DGg910OjpXuvDNwy29AF834tJvmAqODhGV3HQkMfhmoa5GrUEbjHqEJLuZybhJmDOh0ufNjh0tPm3ndcrKt11kJDKsgDmNG+miRml46xy0NjB1q0rjKjPi409qw1tbnWm/nhHZBUgY9sA5ocbbSiGkxTi0SDN3gyCiZDUW0WyTanzEdVihYBXY0N0xd4rQxd5qxwJUAXPC2Ml3G1V8rQSktzjJv94Bm5meXG+h9h5I5PBRCElwWU1xYcflHwbeJFLt7apX8FlCmNZEh07gGPKnvSDjh4iCzBgvS2zk/Zv8nAhpvyKZ8CZeHkM0GMNdzlxsY+jrQ77HIy0LYEieCS/Mycg6c9VD4x0qpY79eNM7GOYW08zPdthBI0WgfdhVV62wORloQM4nIWhCJxtWyX/j7WgmBJllAt2sD5jCZ4LrhZ7Np0GkWBpdbCRvufDnQxmJOMUDuCw01SCYUweVKMPeyo6ExNDwt+DZWhJu+KLWu0ikwIrsfCclA2xAzyDjJUHPjHCC9LXCnw0H+bYcEptBhbiC9fWhl0uJOcv5euRKAKz+X+LiYovH2K91EQkS8CC6jv2ngU0rPximP7gj8rt3SflsxhYjRL6XVcuq00h/HqqI+uAG4fbHjYqLN+cClNql8ItW4/psLZJwZfDsrrDQu7gMuGg+cr0aByUwmzAXUN0zS6XTWY6vgHONy2XbDLDKVA0CCRBAHqx+687nNjZfsmKFi54dsIIt61CA6A2lNQjsnY/JOgVRwVQFqNbefWcFIlWrBt6nVXPm2C8Yaocm+kZgCHNqs/M5dD9w8V3nnwuHJPKBxVvDtAKDjJUonjHxcMSM+BJe2Zy5C2HveAqSqCTt63mq/f1WNX0fmI5OZCIzzL4nGuNPlSu/RSQLQtIb2ee4A4IK39PtY0efuwOSX2qSvTn1+yTWD9ygfWAdknme93krjMjVukutgLuDs/wDnPR9Y1uN6+/IIrBqy2/4wbCdpOBqdBgz/n3n5NdOAhzYADxtmU25xtlzzkMHcZiEtmxyyUYgZxt2J1gE9ds8HCGhcsnJoMfkILeqezFQo6oE7hFRYtTKCl0M0/HaCq6MmnVtSasCC4CsBqlRX3rlwqV7fWXovUT4nUYUdbQYkE2ETJ4JLfRk7j1BMCQLhOA4lUk7WuHWyqFyPaActqi+2aFBaD7Q/j1PfRxdNmLTVPu2HK42+yESuFYbnPQ9c/1Pw8zhJLVSjsf7+GrHS9IzH1l6Hfx+m3HtteHXHi52Z1qyc3kYzrGWjI7mvjU4DqtaUB5xY9byNofIul7ls9Tvq/9+9QgkmEVwrSegLKH6jLDXZ7MFNwI0W2U3ODjJrtujkCeGSdRPwny1mH60xK72s7jGXtakQMAuujpdYl0sM7DY+I22WEpeF4LrOYphocs1AhK1sctWwcPCe+Msp8XEJet2ufDc5PTrFInTEieBSX8Z6Bhu2eKm0gkuYas5+RPGFGXvlADByKjDqO/2yWxeYzXQptYEbfwF63aE0ckDAnDFAMr27YzSNRFIKMOQ1+ToAGDFZsadXrRUQXFUMjnQn49yc+m7sTF1WgtVO4/KX1Wbm3WBYXZ/d7MIymvZWGpIH1uoF1qhpepOulbYiTMJiWIDWxzX8f8qzOssgWOq21jfwrQYokZdPHgxEPF7/kyKQhEZ9fF+gN++uEvDpJqaYr/HGOYr57OZ5wP1rAu9Ki3OU7zaDgWqSCE6TtmR4tn0fUI4py6spNHftMS/7yN6E7U5UhPhlEwPLRkzW+1ytBFfdNvJjJiYH6pdsctVwcPKeiOhM26hCcT8p6WdpEIORvWEgenimeZHUyqE1YZ39EDD/eUWLMgo6QTvJINNG3eQmnWa9lU/xSSWwQ5hqxEvmF5oMlpVU2+C3u0CJgNLlNBP7MXMvrlaLwO9itVcZluByOB4n1GleACDTEOmpvd4qaUDBoeDH6HO3st2qqeZ1VhpXqIKrbmvgonfNy42Rqlrh/dBGYOGbQOYgoGFX4PG9wKmjwFsd9KbClDqKZpynyYRx1VfWZXEnBszRQhAI8/O5zwSOm5SqCMOFbyjRlEbzVJXqQDONFivufeuBisBLqR04n+78BiFj7JScO1b5ls2ynFTNfEzGrDVeUea6rYHDqmm29SCzhiaen7GsduY70Yk4877Asl63A8s+BC75APj+buvkvjKcCC6ROsvk49LUR3E/yyqJciUjPgRX3/uVSm30iTCJ4Go3XMmaEO2US0mp5kiqR7YHzpNUzVmvb+QU63WMKY1f7zuBv9TEwtqXoVkfYM235sGiTq7VscblIJehEdMQBa3gcuCUB4A+dyn+S6ngstD0tPemy5XW98HfiDgrClwu4Mz7lewk1RsAQ18NrEtKBQrU8fYut7ItEPAFiUa6RjOgrXEaOwtE4yaCVepmBq7NnaQIuJvnA+ltAufzl9XwCl/wFvBnSyDj7ECgESARBsZX3+Iei4TQSdUD9bvRaYGyafcXZubLJwHf3qA5l+a5iHooe1bMwgTntmmmEpMVDVbLkFeUDwB0ugx41mJwdbiI6FiTj0v7n5LilibxYSpMrAqc87D55RONhVZLcLnLLtllSu2Ar2DYG8q3P5pLU3G1AQkytL0ydyIwWDNfl7Zxvmi8Ym4xCgM7TUOUx+n0L8Ec7he8rUTqZUjGUQm0ZRaC0NTzNPRIw4nI0u5z0bs2fsUwGpFBz1gnMRbPvNuowLWKBjmUcV2t1ez7YoYDkWorIVkTQao+jyY9Apr2hZpZv411Pa2RUn+Mjb0xWtL4nIOVu706/vHSiYFgJ63gYkzRmq+dYdaktM+pTqbyLQunF/XYqGEZ6/d5LwBn3GtfXv++EbYFTdQ6oO3Qiftv0rhk95A0rtIgPgSXFaLR0wquWI1e73ql0vNr3se8rrMkY0KnyzR/NKZCI9qXIbGqYm4xYvVyDh8HXOMwQk5gJbiEiS3rBuD0m4DrbTJXCB9Rn7vhWGiI5/bQJtlK+30AhxFuUWpEktOAJ/Yr4/eMgiuUhrLfGOD+1YEgE5ENPbGqWXBp6X5dQHg4re/GIAGjqdDKrD7sDcWEPvx/wJ1LgS6aSQ91ZWOKNtiqv7kB196Tuq2BMbuA7qNhwsrHZRRkZ9ytTCFSmty+SPGJDnlF6axpfeVV1I6LsbOlvU4yFZYq8WEqtEL4bbTh7qH0eKsxe0UAACAASURBVEsTYznaD9enH7rso8AgZ7/WYRHZFfxk8sU9RocebZV5vhKtdeXnwMwHAuNkuo0yb3vbQrkpUHQkirW58bjh24Lq9YHEVH1ePatB47L7ddG75qjA0qgTwjwtGi9u7Hw4jE7TRkZqNa60Rkq9Pi9IhnSngtLouxTCITFVCU6xSqZ7+s3KBwDqGbJaaLNcWNXTdhcAvW7TLzNGNPqPod47kz8uijNQA4oPU5tIWpY1JKkacJM6W3VjNb3VkNeUAc/CXGscXiK9ByS4SoP4FlzCVKgNYIj1IMB+jwGHt5nH21xpmDdM2piGKbjsGuZgg5WrNQBO7A/8r14fGLNT+d33AeD7O633lU4BAqUxBJTs3pb+KZsGXiu0ajTVZy2RHkODTMCWJuc9p2hbIlt/JEJS+LgSkhXBqB04b4XT+m4cgyfqBffKrQROSGukaIy/PKEEr8iw8+la4XIrnafCo8rkl9E0/Y+apqSJelEN6nh0h0XGFcmyXobxosak07FueyoRcW4qVAVXlWqBdEmxTnRZp5WSZVyEz4dL+wuVb6veqQ4xxqwBcO/f+lWy3mrmecBp1wBXf6OU1Qrx8gq/hFNkGpdRIxEvecOuyrdWOGk102u+Ce3cdpSG2SatETBiknnAejjnEpkrQplNwGmjbhzrKOqFURO76ivgpl+dn79mM+DKz0LPCRkMESUYjWlV7loeSFaQmKp/VlVrmaN0AWdannGwM0UVlhnxrXG1OR9Y/2MgPBcoP6bCUGiuJgrVpty5bCJwIleJaguGuOakVLMTXnY/nAoD4eBv0Nl+OyOiYSg+CZMW2SRLyR7f+w7l/4hPlOSvWgF9w0/Aa62U31a+l5AoyzE1EdS/gU8Dfe5Rpmt3fLowTIWdRygN8/kvmTUlp5GQpU3fB5TsODKhEirpbZRIyM1zFC3OCU6GmFSprgjW3UvVfWRRhSS4SoP41rgunQjct0rfO4rHHk6j04CnjyqObUFCFaCmZBp5KYaX5IbZSnRXKJjGyCHQKBpzRQZDTHkhsgdoy+ZyA4NfDFxbcpoyTk6Lk+lTQiEmnZlwBlq7lcHCoeDUwiB8nT2uDwwC7nOnEnpfHmEsOkJLIDpATucXczp57U1zAsEj2mchIiszbZJCE2ETkcbFGLsPwC1QWs4POedvR6VUTklMts6BFm9E0rgKs2SrAcp3qD6L/2yW5zEUL2+okwmm1A6MrVkyPrR9S5N47NQEw2m9EdvV71R6ZSnP9P8vULuVPm1Xy/7W24cy63pCFaC4RC+4Gnc3jy8jokbYgosx1gmK0OoJoBjAbMbYLM755mgVLsyCxfT0MSG1rqJ5BpsWxYpq9eTLhSYbySy4l3wA/Dk+dF9Fh4uAff+Gf14dZVgnROMVzkDu0qTH9Yqfq+vIoJuWObf9ARywmSsMUMLTI+l4JCQpUbaCR3cEgohkhBLJ6A90CXH2bCJsIjEVtgfwF+e8gHPuAfA7AJssm2WEVeVOrmE/aDbeqZUR/bBh4Tt0FCBiQY0mwJCXQ48Mu+JT4L5/wj+vlDLQuNIaKRMcjvq2dM8jkjM7DYpwJwLdrim7wfmh0LArcNpV9ts06GwdxRoOVWtZR6sCoWlcIkl3eZglvZIQialwDYAXGGN1AJwCMBSAaW5wxtitAG4FgGbNQphUMVSSawD5e6xt/mN2ld65KyrNz1DmPes8Ivi25Zmy1MIZA855pPTPc/4LQP/Hox/NRyiEIuAHv6Ik9Y6mT46wJWzBxTlfzxh7BcBcACcArAJgSiXNOZ8AYAIAZGVllV6X9+qvgfU/KFNzENFBzOhaUahIPi6XO/hcb7Gk1YBAkE5Fx50QSNtFlAkRBWdwzj8C8BEAMMZeBJATjUKFRc2mSqJWonLB3Eq2ffuN1O8KJLjKO1bzjpV3jFlbiHJJpFGF9TjnuYyxZgAuBRDmEHyCCJOnDwffpjIG7BDhccdiZVwhUa6JdADyNNXHVQLgLs75kSiUiYg2I7+gdDQE4YTaLZRPNLnqS+t8m0RYRGoqrMBhehUIkUevslORfFxE/CCbwoWIiPjOnEEQjiAfF0FUJEhwEQRBEHEFCS6i4kOZugmiQkGCiyAIgogrSHARBEEQcQUJLqLiI6ZJqVmKKccIgigz4nsiSYJwQst+wMipQOtzY10SgiCiAAkuonLQbmjwbQiCiAvIVEgQBEHEFSS4CIIgiLiCBBdBEAQRV5DgIgiCIOIKElwEQRBEXEGCiyAIgogrSHARBEEQcQUJLoIgCCKuiEhwMcYeYIytZYytYYxNZYwlR6tgBEEQBCEjbMHFGGsM4F4AWZzzTgDcAEZGq2AEQRAEISNSU2ECgKqMsQQAKQD2Rl4kgiAIgrAmbMHFOd8D4HUAuwDsA3CMcz7HuB1j7FbGWDZjLDsvLy/8khIEQRAEIjMV1gJwEYAWABoBSGWMjTJuxzmfwDnP4pxnpaenh19SgiAIgkBkpsJzAWznnOdxzksAfAfgjOgUiyAIgiDkRCK4dgHozRhLYYwxAAMBrI9OsQiCIAhCTiQ+rqUAvgWwEsBq9VgTolQugiAIgpAS0USSnPOnATwdpbIQBEEQRFAocwZBEAQRV5DgIgiCIOIKElwEQRBEXEGCiyAIgogrSHARBEEQcQUJLoIgCCKuIMFFEARBxBUkuAiCIIi4ggQXQRAEEVeQ4CIIgiDiChJcBEEQRFxBgosgCIKIK0hwEQRBEHEFCS6CIAgiriDBRRAEQcQVJLgIgiCIuCJswcUYa8sY+0fzyWeM3R/NwhEEQRCEkbBnQOacbwRwGgAwxtwA9gCYHqVyEQRBEISUaJkKBwLYyjnfGaXjEQRBEISUaAmukQCmylYwxm5ljGUzxrLz8vKidDqCIAiishKx4GKMJQG4EMA3svWc8wmc8yzOeVZ6enqkpyMIgiAqOdHQuIYAWMk5PxCFYxEEQRCELdEQXFfBwkxIEARBENEmIsHFGEsBMAjAd9EpDkEQBEHYE3Y4PABwzgsA1IlSWQiCIAgiKJQ5gyAIgogrSHARBEEQcQUJLoIgCCKuIMFFEARBxBUkuAiCIIi4ggQXQRAEEVeQ4CIIgiDiChJcBEEQRFxBgosgCIKIK0hwEQRBEHEFCS6CIAgiriDBRRAEQcQVJLgIgiCIuIIEF0EQBBFXkOAiCIIg4opIJ5KsyRj7ljG2gTG2njHWJ1oFIwiCIAgZEU0kCeB/AGZzzi9njCUBSIlCmQiCIAjCkrAFF2MsDcDZAK4HAM55MYDi6BSLIAiCIOREYipsCSAPwCTG2N+MsYmMsdQolYsgCIIgpEQiuBIAdAfwPue8G4CTAMYYN2KM3coYy2aMZefl5UVwOoIgCIKITHDlAMjhnC9V/38LRZDp4JxP4Jxncc6z0tPTIzgdQRAEQUQguDjn+wHsZoy1VRcNBLAuKqUiCIIgCAsijSq8B8AUNaJwG4AbIi8SQRAEQVgTkeDinP8DICtKZSEIgiCIoFDmDIIgCCKuIMFFEARBxBUkuAiCIIi4ggQXQRAEEVeQ4CIIgiDiChJcBEEQRFxBgosgCIKIK0hwEQRBEHEFCS6CIAgiriDBRRAEQcQVJLgIgiCIuIIEF0EQBBFXkOAiCIIg4goSXARBEERcQYKLIBzAOceUpTtxtKA41kUhiEoPCS4iJqzYeQSf/7Uz1sVwzNq9+Xhi+ho8/O2/sS4KQVR6IppIkjG2A8BxAF4AHs45TSpJOOKy9/8EAIzq3TzGJXFGQbEXAEjjIohyQESCS6U/5/xgFI5DEOUWj88HAHAxFuOSEAQRDcHlmG15J3Hl/y3RLbugS0Nc2ycDp4q9uH7SMtM+l/doghFZTXH4ZDHu+HyFaf2o3s0xvGsj7D16Cg989Y9p/S1ntcS5Hepja94JPP7datP6ewZkom9mXazdewzP/rjOtP6RwW3Ro3ltrNh5GK/O3mha/9TwDujYqAYWbT6Id+ZvNq1/8dLOaJVeDb+uO4APF24zrX/rytPQqGZV/Lhqr9R09v6oHqidmoRvsnfj2xU5pvWTb+iJqklufLZkB2b+u8+0/qvb+gAAJvyxFfPW5+rWJSe68cmNPQEA4+ZtxuIt+v5HrZQkfHBtDwDAK7M3YOXOI7r1DWsk4+2R3QAAz/y4Fuv25uvWt0xPxUuXdgEAPPbdv9iWd9JUPsH9X/6NfccKdcu6N6+FRwe3AwDc/tkKHDFoO2e2rot7B2YCAEZ/vAyFJV7d+oHt6+HWs1sBgKneif1zjxfi0fPb4eZPs03rtXVP1I31+/L9xypPde+V2RvgdgFVEtz+9VT3rOteh0ZpeHp4RwCxqXuRtnuVnUh9XBzAHMbYCsbYrbINGGO3MsayGWPZJSUlEZ6OMHKyyBPrIsQt01bm4PO/duGf3UeDbsvVb1YGGldhiRfHToX2rqzecwz/7D5WSiUinOD1cXAefDsichiP4E4zxhpxzvcyxuoBmAvgHs75H1bbZ2Vl8exsc882lhw+WYzfNuTish5N/Muem7kOHRul4dLuTWz2dMbuwwVoWjvF0bacc4yfvwWX9WiCRjWrBt3+P9+swrcrcrDj5WGRFrPMyRgzCwBiWvYRH/yJ5TuO4Ovb+qBni9q2287fcAA3Ts5G/7bpmHRDz1ItVzj3pjzcz3DYuP845m/IxR39WsW6KBGTMWYW+raui89v7hX1YzPGVlAMQYCINC7O+V71OxfAdACl+0aXAjd/shwPfbMKuccDpoKPFm3Hg1+vivjYM//di7Ne/Q1/bMpztP2W3BN4Y+4m3DFlpaPthfkmks5HpOw7dgrXfrQ0pkELnHO8NXcT9h49FeJ+zrf1Ki4uuF2xDcR9ZfYGzF6zP6ZliCYXv7sYr8zeAK8vvlUV8Q4u2kLu/rIg7LeQMZbKGKsufgM4D8CaaBUsVA6fLFZVdY4ez83Fl8t2OdpvrWobZ1BMQNF8gVbvOaY7RzDEqU8Vh2b+K/L4Qto+mnywYCsWbj6I71buiVkZ3pizCf+btzlsX4AT659XDc5wx3gAyfsLtuL2CuTzOKX6hnxxbmM7ZfBxEaVLJK9hfQCLGGOrACwDMItzPjs6xQqNE0UedH9uLp6buQ4eH8ehk8V4YoYzGSoafSGwDp8MXXMYM+1fXPreYhzIL8RZr87H9oOKI9jNhDAMTbCE+g5f+t6f0uVeH8fcdQdKVSOrlqzE90TL11ZY4sXj01fj0Ikiy22W7ziMzmN/wdECpbMy/rctAJR6UFjixQmHZQnlWZd4lXuYEGONKxgZY2b5619ZsyX3uOmePj59Ne6cElzQxrvGlX+KfM1lSdhvIed8G+e8q/rpyDl/IZoFC4XjhYoj++c1+/wvgLETfe1HSzH2h7UA5KY1r7pMa/IqdqjJfLl8N1buOoo5a/dj9+FTmPCHEsGV4FJK4dG8lJxzFFhoVOrmOFnkwZ4QzF7r9uVj37FT+Dp7t275pMXbccun2fhREvHV9Zk5+O+M1fD5ON6YsxGHThSBc47Ji7fjoEFonCr2Sk2BXh/HN9mKudKJsPD6OPq99htmScojmLYyB18s3YV35ivCaFveCcz4W6/NjZ+/BccLPVix84iup+71cQx++w90evoX03GXbT+MbXkn/P+LPT5sC6GBF3XB5TKrZ1tyT1g+01gwf0Nu8I1KgXPf/APnvx1wcXPO8cXSXfhpdXDT5s5DBVi71zq4JPd4ITLGzMIva8vWTPpvzlFkjJmF3YcL/MuKPT7c9lk2NuwPWFLyCwPBNMPfWWSKMiSiS/nuPlqw+3ABRk1c6hdYAgaGEm9gvM2KnUfw8aLtAICFmw9i8p87MGXpTrR47CdTj96r9qi1ZrdJi7eHVK706lUAAHmqv0z4Q7S9yW+yc9DhqV+w61CBaX9hstp7rBBnvjw/pHOP/ngZHvn2XxwrCNyT/WqI7/5jeiHIOcexUyX4/K9daPn4T3hn/hY8MX0NtuadwNgf1+GeL/7WbX/Je4tx2rNzASgCKueIUvYZf+9B7nHlPp4q8ZqEg5FP/tyBHYcK8Oi0QPYJYydC9FyTEpR7N/jthbjfEGouZIeP6+9tiZdjh+S+AsAV/7cEA9743f9/20HrcsooVuuV2yC3vD6Oc9/8HYPfXuhflnOkABMl4ee2x/f40GXsL/j+n8hNrrH0eeYdD7xXt33m3KR5/tt/YNi4RZbr1+87DgD4bIk5bP9kkQfXT1qGLbmhPVMniPZj2fbD/mXr9uXjl7UH8Kgmi0q+Jgp09Z5j2GrzHhCRE5eC6825m7Boy0GMm7cZW3JP6Exr/oaMKdkZnp2pHx/zyZ87AAC5x4t0vSKhcWkF1/78QvzuMLACCIyhOXhC0U4S3Ga/2d+7lfEoc9bJeo7hh1pvOqC8KD7OMWnxdmSMmeXXgrwGxfFksbk3WOjx+n1seQahvmH/cf/vkROWoO8rv8Hj9eleTreLmYRDfmGJX3Oct/6A/1lotSRjG3uiSGkAUpMUE2SxsfAIDALmnJs0LqfkHJZrtD+v3icddyM0LmNwhugo7dL0yG+YtBzPz1qP3Hz92CDB19m70eGp2dh04Li/zHknipBf6MHLP2/wb7f/mHx/AJi8eDumLJWnzIq23Fq1+6jJBHiiyIPVOfbh93PWHbBdH4pVIdBZUS5u79FTmP63ou3/uv4AFmzMw1tzNzk+nlNEx0x0SrVlWJVzzP/88g2daBqoXrrEleD64Petup7shwu349w3f/en49mfX4icI8rLYFVtThYp2x4v9GDUxKX+5V6fD9sPntTZ4yct3oHRHy/zaxhGOOdYsDFglpmsCsX8whLMW38AbpdZcNVJVV6AAonwMDqo52+wf/FleDn3DyYVAlQcd8HGXOw+XGDpPxqvmufsBMCaPYp55P6v/sF7C7b6lydoTGi/qffkgnGL/JqjtpHSXqbxmsXzSUly6xpGrRYhzHXFXh+0RfWG0GJrOyii5D4fxx1TVmLp9sPwGe5BkUcpl9sFTFuR4/cjye6VaMS05SnyeHFEbfxX7T6KgmIvho1biFaPK9q/6LFrB8L2fmme7rjaezD2x3V4Yrrcj8sRHcnl83F88PtWXPTuYlz87mLdujs+X4Hh4wMmsWCmMZnZXSaYPV6fyZICBO7zTlWjHjnhLzzw1SoUebzYcVBZ5nTYSSgIs3miJipH+xw++H0rij0+3DhZP8xH9n4T0SOuBNfLP2/A87PWm5af+2agl/+06sfS9ni0FU0M7Bw3bzOyNaPxvT7g0Wn/4kC+uVE/ZVEJ56w7gOsnLff/F9rZtryTuOmTbBxUe2taH5ddR8zj1Tc4N07O1pkU1+3Nxwuz1mFL7nFc+9FS4+4AFA1gq5olIMEgOK+ftBz9X18gdSTvP1aIH1bt9R9DhtZcYsyUoG2/b5i0HEu3HdJpITWqJvp/ayOw2j45Gz+tDhxL3B+Pj+uEnfYean2BWsHhsSi3DCGIBIdOFKHl4z9JzwdoNS6Gh75ZhfPfUnw5MmHJJN2mu6asRLfn5oJz7heaIuDjSEEJFm0OHkbtVKGcu+4A3v1tC77/Zw82arTlYCzYmOs3Na/fl4+Wj//k1wC1zxIA/t6lDNoWGvH9X5ozh2hxOqD67i/+Ruexc/z/t+adwKETRf53cM/RU/hh1V5/eQqLfXjrV0XT0taxxVsO2taHF2atwwXvLLRcLxDPStvB0j6H3YcL8Ne2Q6b9ypPPsyJSpimfygLRZGgbR23jJsxnxgAErVM5GMUeH+ZvyMWmA/aNgmjUPF4fZvy9B+d1rK9bX1jiRbsnZ+PJCzrgpr4tpCHBoryFJV5c8t5iFHl8WLj5oM58p2X0x4H0MTJTpcfH8cdms/lTezyjABVcITGhCXYe0gc6HDKYlt60MON4fRyvzt6AAe3qITnR7Y/E/GvbId399fo4EtVsRi7/NofRrkFaoNwhmAqNQwiMJmWjJmXUGESD7bW4V9pjPPjVP/hVTXn0+V87TRGYM//di7d/NacLA5ROl8jW4ZFEpz4+fTXGDGmnW7Z8xxEs3xHolF3QpSF+WbsfM+46Ex0b1ZCe51hBCa6ftBx9WtbB1Ft749MlOyyvCwi8Z6LKzg4SNHH6C79i+0tDwRjD6pxjGD5+EV6+tLNpO3EcxQwMDNSYngX/ajKd7NZYQ0SnZ8nWQ7hm4lLcOzATDw5qIy3PhwsD/uupy3ahfcM0nNa0pmk78S5sPnAcZ7auCwA6bdzlYjptTHCyyINdhwqQe7wQM//dh6eHd4ha1pUVK1bUS0hImAigE+JM+XCID8Aaj8dzc48ePaSRRnEpuOyippIT3aZlMj+JVcMvw1jf3vp1E95fsBVt6lez3U9U+uydR/DJkp0Y0aMJGtRI9q8/qvZu/+/3rbipbwup2UloPwNeX+BvbIU5VIbwdQGBBt7r4zqn/2u/mPPe6codYvg+AH/DLNBeS5HH6zfxyNhxqADtnpyNL27u5Rf2Rt+iVigJE+z0v/dguibi8HhhQCAcLShGzZQkAHoz5ZtzNuK+c9voBJGPm7XMEp8PST4Xbv00Gx0bpaFIXV9Yot/OzjwpNKrvNGWctXqfqaET5lcZ7/++FQs25mH/sUJc3K2xaf0XS3dhWOeGlvsDAe14x8ECneC69dNsHMgvxPd39/VroBvVzkLDGvLMLcUen//+AzCZVO3IP+XB3mOn8N4CxSSt1eCNeHwcpz0zR7puo6ZDI4tm3Z+vPO9x8zabBNesf/fpTJE7Dp7EY9+tRsMayfj5vrNQMyUJU5ftwoZ9+Xjmok7+d2Hsj+swtEtD1KuerKvbCS7mDyTS8v7v27BKI2DvG5iJWqlJltcbCgkJCRMbNGjQPj09/YjL5YrvcQQSfD4fy8vL67B///6JAC6UbRM30lr7gtiZHWQj10s8kT3bcfO2YLPmZdmjCg475zkQEADC3KdtvLmmscw9roSiyxrBO9UsGns153I6Tkkcben2Q7gviClHS4mXY8fBkxGNzdJqj23/62x435JthyzNlMdOlWChqilamW613Dg5YMLVRmiOm78F8zfk6kyFV/zfElPIttfLcfoLv2LehlyMm7/FH9Vm9OXIOht+c6fkWprWSjHdV6trBoBXZ2/Esu2HsetwAcbNk2tlTjXNEq8Pl73/JzLGzMLSbYcwZ90BrMo5hk0HjuMLdcC+KLNbEvYPAD2en4vrPl7qV7l+tsjiIYts3Jx7HEP+t9C/j50GUuTxSYOIAGCpRuCNnPCXab32dr4zb7O/7dh9uAB3fbESYzQJj/eqEbf7jhXitGfnYsnWQ3jsu9X4RI1e1FofhBm1SHOC3zfl+afo0bLKkP+yJIzOoA2d0tPT8yui0AIAl8vF09PTj0HRKKXEjcZ1MgKbsUzjCoUfVu3FD6v2YtodZ6BH81qOMi0Agd75cbWhWrbjMJbtUF46r8+nc+DmHDkl7b2GEnllRAiPUAdHerw+9Ht9AQDg9RFdwzp3iY0JzQoXY5adDCF8pt1xBuY5GKe0ctdRTF22C1f1bGZad4skE7yREp9PF0kn0nZpBdfsNftMY+R+Wr3PH2BR4uUmIXuy2OMPQBGEErkqw6nWU+zxYYXq171S0+Cf91bATC6EoCwq7kB+IY4XerB4S8Cn8/j01bi6l/4eFxR7pBYNo6Vg2kpzxnmBXWi71fhKH+c4UeTB4ZMBN8AbczehapIbo8/IMEX+AeZB6Fd9GLgvs/7dp+tUiPdZ64+0syRosTK/h4mrogotgXp9lopVXAiu7//Zgx/VwIFwkE03Eg4/rd6H2z7LRpNaSvRSfqG9QLDTWMbN34IeGYHErl4ft+w5Tw5xPJlADPQ1TskQjAJN4/yfb8LL2RjOfgkuFrRnKuvdWvHYd6tRTxPGHApWkZULNY3W7Z+bc0reqckz6fH5dIFDgBI16VRjdsoNGu3SDidpiQqKvf6EvUZ6vThPunz63zmoW60KDp4oQuOaVXHVhL+wShIqb/Qr22GMYnSC18vR56V5OpMxAGVowvEiDOnUwLTPg19Z19O7vtA/3yKPFyVeHz5aFPr7+PGi7bixbwtHybMrAyNGjMiYN29ejTp16ng2b968NtT948JUuPnACZMPJRQ+lQxaDIePFm3HwRPFjqbBAIKPY9EGUhQUey17zmMlczWFQu5x5w0GEP1xQE55Y+4mHLAY+xQuN30S3mwEVj3kUAJAlm0/bNKYTxZ5IrIeRIKdhhMJD3y1yi+U9hw9JRVaADBxYXgdMKeM/22LSWgJVucck1oBQrHGFJb4sH6fs7yjRiYu2k6DkjXceOONB3/44YewNYq4EFzdmpmjfcozRtOJEwqKPSGNQypPJEUx86zWDBUJvYJMUxKMUASUFbKhGyeLvTGbQ+3fIAOGS5v9DjolCRa+NSfYJZt2u5jjFG7Wx/fiwvGha4KCprWiP84sVmzcuDGpRYsWHa+88srmmZmZHS+88MIWM2bMqN69e/d2zZs37/Tbb7+lPPjgg42eeuopfyh1ZmZmx40bNyYBwJAhQ06kp6eH/SLEhamwWpXSKWb/tumonVol6j3RG8/MwBdLnWWnFyzdfjjadvAyo0qCK2I/YrS58vSmyC/0hN1DvvkTZ+Y3Le+qyX7tOHiiKCz/X2XB6Z0Z0aMJvpHMymwFY8Aoi7GPANC8TkpQf9XdhlRoRkb1bobP/7J+7xvWTLZcFy4Pf7uq6ab9x6MqEds0qF7w2uVddwfbbvfu3clfffXVth49euzs0qVL+ylTptTJzs7e8MUXX9R84YUXGnbp0iV8B30Q4kLjqppkDnGPBoM7NcDoM5r7/5/Rqk5Ujls7NXS/ymu/bPQPpJQfMzqhtJfZTI7ZsEZ4L5YsHDjWuF0M1aqEX2+2aqZ679JECSEPpg0EG2YA6HP5hYpsnFEs6NG8Vqkd22naLrFVzwxnmvVCBwO8gxHMR1itSqLtepESrqLQuHHjop49e55yu91o06bNuL3vRQAAFrtJREFUqQEDBuS7XC507969ICcnJzznskPiQuOSjc2K1nG10VM9W9TGn1sjN1VFS8hoMZo5BrSrF1YW8Kt7NbXUMKuGeZ+Ng41DwUlP18jUW3qjQY1kDHrzd0uTntvFkBolTb26OnVLx8Y1sC33hD9KFAAS3axMNCgXA4Z2buDYv2rHpd0a68aWhUpKKXUkAeU6tY/00u6NpXO9nda0Jg6eKMJTF3TQ5ceMBi9d2hmPaULmndKuQXXLdaLzE22caEalRVJSkv9JuVwuJCcncwBwu93wer0sISGB+zTBVkVFRVFL4Fj+usoSkiPoqQzp1AAN0uSaRI2qibrQdm1whGxUfyg4NW/2alEbF3SxH0AKAPcObK37b5cB3MonuOPlYejR3LqHKpuyo7QZ3qURbu7bwvH2TWtXRZ9WddCibiq2vDhUt+7M1gGN2cWY7fioUNA+yyqJ+leme7Pg2kdSggtPDG0vjWpzyh+P9A/bl9jT4O/77wUdQj5Gd02dKk3BdW3v5mhRN9X/3+o9qproxuQbeqJlun0SgFCZdMPp0iEUTmjboDqu7d0ct57d0r9M+ForY9LdjIyMon/++ScVABYtWpSyZ8+eqGlh8SG4Ep0X05hFoFZqEsZeKH9Ra6cm6UwT2p5e9wjNId/deYaj7To1roH2DdNst9nx8jBdaiM7nr+4E165rIujbY24o/By7Xh5GJrU0of8Xt2rGa7Maqpr/AQuBlPKIjsu6aY3dWqLfMc5AeHuYgx7j0YnQrF6smoC4txk7rEaqKuFc45bzm6J285pFXYZmtRKweVZTTGsS0N0ahy8vgBK/f7m9j748Nos/7pHBrdFrRSzSWvJYwN0gt/I6DMy/L9TkpxrsqEGKjWoURWvXh6ov1Zas1W37RJDdhFt/kInCHOwU+tDq/SAkE10Mzx3cSedSVdo607qSUXjuuuuO3LkyBF3u3btOowfPz69efPm/hdy+PDhLfr27dtu+/btVerXr9/lrbfeqhvKsSMWXIwxN2Psb8bYzEiPZUWVEExY46/upvt/ekYt01QUglopSTpTkzaqT9ar/Pb2Pvjy1t66Zc9frB/cPbxrIwBAm/rWZgMtqUluVLHxES18pD8AoKahsalbzdx5+b9re2BU7+ZBM3W/eYV8UHEwufXj3X3tN7A4Ttv61fHK5V1Qr7pZ83W5mOmlTq9eBf3apkuPff/ATP3+mpNptaHmdVLw0HnyPHXKfparTGh7/S3VhureAa3x1a29HTVIohHsGqG5qFqVBLx7dXdk1EkNuu36Zwfjr8cG4vSM2qih1p0EF8Od/VqDMabTCgAgJVE59gejeuiWi7RmWgU/FJ/zoPb18d413YOmpRLUrZakS4mVYvHuW1kcXr28C2bdG6inbR2+h5n1lO2EmT/7v+dabjvznr5Y9vhALHq0P36+72z/ctHOaOuk6BhHo1NYnmjbtm2xdvzVtGnTdtxwww1HtOuqVavGFy9evHnDhg3rvvzyy53btm1b27Zt22IA+PHHH7fn5eX96/F4Vh44cODfBx54ICQnZDQ0rvsAmON+o4iVxiWrC9o0MiufHISLT2ts6VSvkZKo17g0v1MlvcqsjNq6MOuFj/TH5T30GoDXZgDtNoNpCwD6tatnGdwwokcT/1QNxp5jH0kgiTAltUyvZttjvNQiQCOYOaOz2vC6XQwjepiPsXjMAOl+drZ/F2O6Z/b+Nd3xx8P9/dpv39b6jpjRnKltELSmtPYN06SCcvtLQ7H9paHY9tIwyzIBQOOaVf2NWJraa+YAxo3shlvOaoE7+7dGr5Z1ggquh89vi29uV7RvxhgahzEA1Rg008QmrFqUp2qSW1evZt3bV/d86hvM58lJLtRMScLgTg3QX9Np6NRYeeYlXp//fbMSJjIS3AxDOzfE2W3MHWpZh61VvWq652g0zQqaWUxhkuh26fIxul0M467qhofPb4uXbMz/Tw/vgE9u7Onf184/2rhmVdRLS0aTWim6eyzamUEdAsm0hfuzgsmtmBOR4GKMNQEwDMDE6BRHjrYir3rqvMD5g+xXOzUJjDFdYzdmSDss+E8/vHNVN6QlJ+rGVmiFWDW1sRrdJxB1COgFY9PaKabAETtHvcvFdIJuywtDbH0kCZrpduunJaNutSp46dLO+OHuM00mESDQaFWrkoD1zw3GwHb1LI8t09hGZFlHHAoz0vaXhmLLC0MwsmdT0zaiUdZO69GlSQ30aqnsq315L1Q1U2PDP6RzQ1RNcvt71ELD6dmiNp6S+Ga0xzQey9gwtmtQHcwgKK1oXifFX+9qqMl6PV6OWqlJeGJYB/9zl/WkF48ZgMt7NMEbI7rirv6t0dZGcL+qMesOsHhecx44W/e/ZV25xnVT3xbYKukcAUDHRjV0wkqU+tz29fDqZV10JtCPRp/u/52oahEeH/d3hsS7oUVYGoyIenZ5D319uapnU8y+X39dX9zcC92b1UJSQuCeJrhceOUys8ARdSoYyYkuXNi1Ee7q31rnu3r/mu74RXP+alUScE4buZZvxErjFPXP7WK47RxFoxWd4cpoKixNIg27ehvAIwAs30zG2K0AbgWAZs3Cc3oyxvD08A7o3bKO3+whlndulIbVe+wHVnZqFPAJ3K76GTLUl79BjWQ8OrgdXpm9QefjSnS7sOaZ85GS6PYn3BQ8fH5b7NVkRHh8aDu8+JMyb1Gwyv/6iK64rk9zVElwI0FtGK3iLLSVPTnRbWu+AMzh2hOuy4KPc2Q+8bNp21/uPws9nv9Vt+z6MzJwXZ8M5BwpwDmvLfAvn3bHGeio3kPR6GsF7nvXdNdNXa6lhaSRfe+a7v4kpELLe/aijujaJOAbuKx7EyzcfBB39GuFK09vig4N06QCRz/vmn6dtjfcsVEa/jdSb0YWMGbed1CH+nhfnSizbjVFcBnn8AICGuDpGbWwfMcRtG+YhsY1q1rmeGxdr5oum4bonAzu2ADvXtMdrTRzggn8PjaV/hIBN/Oevo7N04Ay4zWgaDhXnK4XKtqOniifxxvIrSl7pk8MbW9Ky3ZHv1Z+bVtbl+c9dA5aSYIqzlC1a62pMNHNcOXpzZBRJ1WXX9EpdSQdNEDpIGlxu50LFivTvsyy0715LRR7fXgyjIAYwpqwBRdj7AIAuZzzFYyxflbbcc4nAJgAAFlZWWHHDd9wpjnyrGqiG9/e0Sdo9vE61ar4HdYyREU0zodlFdF0V399hF+C2is9v2N9XCNxRr97dXfdrLRdmuiDFGTzcAGh2cXrp1XB6YboMbeLwQ2GJY8NME1uKHuhGWNwM6C5wYdSr3oVk2apFSJDOjXQ/dcWWxYJx3lAuxXv+nV9MnTbXNytsX8aD6spNgCjcNefSzzX1CQ3Zt17luUxkhPc/jE6y54YiKMFJcisVw3v/qYILmEylGVmEM9IyUF3xBSYYmTcVd2wcucRf37Bge3ro2uTGnh4cFu4XQwrnxyE52etk4aAC9KrV0H2f8/FX9sO+QfFCpOeU0TCWKuI3R/v7ouTxR7MVjO5l3g5BndsgNlr95vqR/uGabrpegSPDpYH3cjqhDbyUdvhEJpVr5Z1cE6b9JATEtep5mxoSqLED35WZl2kJSdiljrRaeOaVbE/v9DUgUpwMXh8XKpVVU104+vb+oRUZiI4kWhcZwK4kDE2FEAygDTG2Oec81HRKZo1ogJPu+OMqAzqE6p/uANpRX2tnVpFqhUMCxLuLuTWdX2aI7N+dTw5Q5mS3SqoxMjMe/raNlxWDf9/zmuDH1ft081vJCOY/DRec/sGaf6xWR012q5u6IF6zZGaUMTeV/dqhkyDxiGeZ7BUWlUSXX7BVa96st83JsyVwr8oFVxq+UXGd7tAG3EsrcZUo2oivtcEvdROTZL6V43UrVZF50sJlSL1eq3GSAp/5rz1Sr7NEq8Pb488DYdPFpumFaqdqtyfUb2b4buVe/DmFV2lkysKZKa2LE0Ur3ZfrRb5f9f2QLsnzZ3UQR3qY+66AxjV29xpvMNhJKesHn52Uy8Ue3x+wTXr3r7SDDFVk9w4XujRaf8d1EjhzHrRDdcnFMIWXJzzxwA8BgCqxvWfshBaADD5BsUGLxrMszLrRjQy/rLuTZBzpAB39muNbk1rml66/w5r76+IMoRpxWgq6NgoDWv3Bk85JJpVF2PopgmlTXBovgi1ty24e0Am7h6QKc0G3rlxDb8J1son9MjgtnhzjjnbxxtXdMV1Oc2RnOjWXY9WfggtM9LxLWL3izQ+FhERKjo1VvEy8x46BwPf+B3JCW58f1dP/3QkgvvOzcRT36/1axNFkmhNca7W9arh1/UHcLZDP8nMe/r6p583YqWBG4kkUq3QL7jsBa2oWy3qpiI50Y1GNavqci1ef0aG3zf6/MWd8fzFwcc/yiJ27z83EAFqJfSshOyH12VJl7eom+qfTDQYVgFcWoFmdayRpzfFhwu36zq+F3ZthPYN00Iy3xLOiYvMGUaCOdeF6u6UpAQXHj5fMWsYbd8AcPNZLU3LtHgtHLDf3N4HJ4JMfQIEevR1qyWhU+Ma/pxnsRy0OOWWXvj0zx14fc4m1LHIBHJnv9a4s19r0/LUKgk4o5X9sIxrejXDtytyMDiCQblA4J6LzsMHo7qjQ0OlsQ2mcQnzUPXkBHRtWhNdDfEm1/VRfH4+H0dW81q4s7+59/7fYR1Qt3oVPDSoDW44M8MUrWdFp8Y1LDsczes4Sz0XibbqNxUGiRC86LTGaFO/um6sofa8Yy/sGPK5tebJ5U+ca5pF2G6g9bMXdXTUUVv11HlSC8pDg9pIIxWtBt87ucWPDWmPuwdk6iIRGWMktCzYsmVL4jXXXNMiLy8v0eVyYfTo0XlPPvlkSGmAoiK4OOcLACyIxrHCwWieWfnUoDJNWFtiMWtsSlKCo8Gal3ZrjBKvzx9xWF81VUWSKTtS0pIT/RpZtOjcpAZ+XrMfjWomI7N+dax55vyIjymEu/ge3CnQ8bDyXQqa1q6Ku/u3xhVZ5ghJ3TlcDN/eIR9QXiMl0e/LcSq0gnFT35Zo1yAN12mmvZHhJDrSitvOaYlVOUdxfsfgHQfjAPkEhyZsK7RCIl0yX5oQOLIhF0ZfqBU1JIOsAeAewzjAKTf3sh2X5uQeu1ws5IHOlZnExES88cYbOX379i04cuSIq1u3bh2GDh2a36NHD8cZA+JS4zLy4iWd8ev6wCR3acllW4m2H1RMPuFOEudyMV2orn8W2goWQnv72a1wTpt03TibSBHmVNmtEoKrdwt56DRjDP85v23UyhIt3C7m2OQYLi3Tq5nC0Z0SSgReWMd3Max6+rxSmxVCy5mtgydsOK1pTVxnGBZT2dm4cWPS4MGDM3v27Hli5cqV1dq3b19w4403Hnz22WcbHzp0KGHy5Mnbfvzxx5rVqlXzPvvsswcAZVqTmTNnbm7btm1x8+bNSwCgVq1avlatWp3atWtXUqUTXPWi1NMNl7Mz62Lqsl0Y1D58Z7mWRLVhCJYTrmvTmshwaFayY8rNvVDLoS8gElwuFlWhBSiBJwfyi6RzLTHG8Mv9Z6NxkEi/8krVRLcpY0p5oCwsAeVJg5lx15mxLoI1M+5qitx10Z3oq16HAlz8bplMa7Jx48akdevWpZxzzjkhzbJZIQRXrBnSuSG2vTg0ahrSzWe1REGxF9dr8sPJ+D5KL5STXmd5pW/rurYZ0+0G/5Z3Vo89L/hGMSDcpNdnZdbFOgfBSkR8IKY1AWCa1uT5559vFExwHTt2zHXppZe2evnll3fXrl07pIzYJLiiRDTNesmJbjxiMQaG0PPAoDY4vUVtx5kU4okEB9ngh3VpiPMiCIsPhxopiXj4/LZS/5Qdn93Uq5RKVIlxoBmVFpFMa1JUVMSGDRvWasSIEYdHjx4d8lw9FUZwfXhdlu1UH0TFxO1ijlP1VETevbp7TM5rHIRPEEYyMjKKfvrpp5qAfloTn8+HkSNHNm/Tpk3h2LFjD4Rz7LiY1sQJgzrUx3kOIqQIgiCI0sdqWpO5c+dWmzFjRp1FixZVb9euXYd27dp1+Oqrr0JyfrOy1FKysrJ4dnZ2mZ2PIAiiIsAYW8E5zwKAVatW7ejatWv4GRfihFWrVtXt2rVrhmxdhdG4CIIgiMoBCS6CIAgiriDBRRAEQcQVJLgIgiDiC5/P56tYaXUMqNdnObaLBBdBEER8sSYvL69GRRVePp+P5eXl1QCwxmqbCjOOiyAIojLg8Xhu3r9//8T9+/d3QsVUPnwA1ng8nputNiDBRRAEEUf06NEjF8CFsS5HLCnTcVyMsTwAO8PcvS6ACj92wQBdc+WArrlyEMk1N+ecV94UMQbKVHBFAmMsWwzAqyzQNVcO6JorB5XxmkuLimgfJQiCICowJLgIgiCIuCKeBNeEWBcgBtA1Vw7omisHlfGaS4W48XERBEEQBBBfGhdBEARBkOAiCIIg4ou4EFyMscGMsY2MsS2MsTGxLk80YIw1ZYz9xhhbzxhbyxi7T11emzE2lzG2Wf2updnnMfUebGSMnR+70kcGY8zNGPubMTZT/V+hr5kxVpMx9i1jbIP6vPtUgmt+QK3XaxhjUxljyRXtmhljHzPGchljazTLQr5GxlgPxthqdd04xliFTOUUVTjn5foDwA1gK4CWAJIArALQIdblisJ1NQTQXf1dHcAmAB0AvApgjLp8DIBX1N8d1GuvAqCFek/csb6OMK/9QQBfAJip/q/Q1wzgEwA3q7+TANSsyNcMoDGA7QCqqv+/BnB9RbtmAGcD6A5gjWZZyNcIYBmAPgAYgJ8BDIn1tZX3TzxoXD0BbOGcb+OcFwP4EsBFMS5TxHDO93HOV6q/jwNYD+WFvwhKQwf1+2L190UAvuScF3HOtwPYAuXexBWMsSYAhgGYqFlcYa+ZMZYGpYH7CAA458Wc86OowNeskgCgKmMsAUAKgL2oYNfMOf8DwGHD4pCukTHWEEAa53wJV6TYp5p9CAviQXA1BrBb8z9HXVZhYIxlAOgGYCmA+pzzfYAi3ADUUzerKPfhbQCPQD9lQUW+5pYA8gBMUs2jExljqajA18w53wPgdQC7AOwDcIxzPgcV+Jo1hHqNjdXfxuWEDfEguGT23goTw88YqwZgGoD7Oef5dptKlsXVfWCMXQAgl3O+wukukmVxdc1QNI/uAN7nnHcDcBKKCcmKuL9m1a9zERSTWCMAqYyxUXa7SJbF1TU7wOoaK8O1R514EFw5AJpq/jeBYnaIexhjiVCE1hTO+Xfq4gOq+QDqd666vCLchzMBXMgY2wHF5DuAMfY5KvY15wDI4ZwvVf9/C0WQVeRrPhfAds55Hue8BMB3AM5Axb5mQajXmKP+Ni4nbIgHwbUcQCZjrAVjLAnASAA/xLhMEaNGDn0EYD3n/E3Nqh8AjFZ/jwbwvWb5SMZYFcZYCwCZUJy6cQPn/DHOeRPOeQaU5zifcz4KFfua9wPYzRhrqy4aCGAdKvA1QzER9maMpaj1fCAUH25FvmZBSNeomhOPM8Z6q/fqOs0+hBWxjg5x8gEwFErU3VYAT8S6PFG6pr5QTAL/AvhH/QwFUAfAPACb1e/amn2eUO/BRsR55BGAfghEFVboawZwGoBs9VnPAFCrElzzMwA2QJnF9jMo0XQV6poBTIXiwyuBojndFM41AshS79NWAOOhZjSij/WHUj4RBEEQcUU8mAoJgiAIwg8JLoIgCCKuIMFFEARBxBUkuAiCIIi4ggQXQRAEEVeQ4CIIgiDiChJcBEEQRFzx/4qKjElmoYTLAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"np.random.seed(12345+1)\n",
"N = 100\n",
"\n",
"mu_1_true = 5\n",
"mu_2_true = 10\n",
"alpha_true = 5\n",
"beta_true = 2\n",
"\n",
"X1 = stats.poisson.rvs(mu_1_true, size=N)\n",
"X2 = stats.poisson.rvs(mu_2_true, size=N)\n",
"t = stats.beta.rvs(alpha_true, beta_true, size=N)\n",
"\n",
"df = pd.DataFrame()\n",
"df['Z1'] = X1*t\n",
"df['Z2'] = X2*t\n",
"\n",
"with pm.Model() as model:\n",
" mu1 = pm.Uniform('mu1', lower=0, upper=20)\n",
" mu2 = pm.Uniform('mu2', lower=0, upper=20)\n",
" alpha = alpha_true\n",
" beta = beta_true\n",
" \n",
" t = pm.Beta('t', alpha=alpha, beta=beta, shape=N)\n",
" \n",
" Z1 = pm.Poisson('Z1', mu=mu1*t, observed=df.Z1)\n",
" Z2 = pm.Poisson('Z2', mu=mu2*t, observed=df.Z2)\n",
" trace = pm.sample(500, cores=2)\n",
" \n",
"plt.plot(trace.mu1, label='mu1')\n",
"plt.plot(trace.mu2, label='mu2')\n",
"plt.axhline(mu_1_true, color='C0', linestyle='dashed')\n",
"plt.axhline(mu_2_true, color='C1', linestyle='dashed')\n",
"plt.legend(loc=(1.01, .01));"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Thu Nov 12 21:13:10 PST 2020\r\n"
]
}
],
"source": [
"!date"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "dismod_mr",
"language": "python",
"name": "dismod_mr"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment