{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Approximate Bayesian Computation: Likelihood-free rejection sampling"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Florent Leclercq
\n",
"Imperial Centre for Inference and Cosmology, Imperial College London
\n",
"florent.leclercq@polytechnique.org"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"init_cell": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.stats import norm, invgamma, gaussian_kde\n",
"from matplotlib import pyplot as plt\n",
"from cycler import cycler\n",
"np.random.seed(123457)\n",
"%matplotlib inline\n",
"plt.rcParams.update({'lines.linewidth': 2})\n",
"colors=[plt.cm.Set2(i) for i in np.linspace(0, 1, 8)]\n",
"plt.rcParams.update({'axes.prop_cycle': cycler('color', colors)})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook we try to infer the unknown variance $\\sigma^2$ of a zero-mean Gaussian distribution.\n",
"\n",
"The data are $N_\\mathrm{samp}$ samples from that distribution."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"Nsamp=100"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"groundtruth=2.\n",
"likelihood=norm(loc=0.,scale=np.sqrt(groundtruth))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"(dmin,dmax)=(-5,5)\n",
"data=likelihood.rvs(size=Nsamp)\n",
"lh_data=likelihood.pdf(data)\n",
"x_arr=np.arange(dmin,dmax,(dmax-dmin)/100.)\n",
"f_arr=likelihood.pdf(x_arr)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAF1CAYAAADSlV/tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl8lOW9///XNTOZTDaSsO+LIogia0QRFRS3LmpPPbZ20y6/UrdzvudrT0/9nva4n7p00dbSU2l73KptVWqrFTcCYTGjEgQEwQURJARkDVkmyWzX748ZxgQCTGCSe5b38/HIg3ufd4BMPnPd131dxlqLiIiIiHQvl9MBRERERHKBii4RERGRHqCiS0RERKQHqOgSERER6QEqukRERER6gIouERERkR6goktEusQYc44x5r1265uNMRccw3VuM8b8Mb483BjTZIxxx9erjDH/X+pSHzbDN40xy7v7dUREQEWXiBzG4Yopa+0ya+3YVL6WtfZja22xtTaSyuseL2NMvjHmMWNMvTFmpzHml0mcs9kY02KMaYyfV22MudYYk9T7rTFmpDHGGmM8x/8diEg60Q+1iMjhfROYAowCgsC0JM+71Fq70BhTCswEfgmcAXyrO0KKSGZQS5eIdIkxZpYxpvYw+042xnxkjLkqvj7YGDPfGLMrvv1fD3NeZ607I4wxr8VbjF4xxvRtd/xlxph34i1JVcaYce32jYtvq48fc1m7fX2MMc8ZYxqMMW8CJx7l2w0D+621+6y1zdbaxUf/G/qUtXa/tfY54MvANcaY8fEcnzPGrIrn2GqMua3daUvjf9bHb7lON8acaIxZZIzZY4zZbYx5whhT1pUsIuI8FV0ikhLGmCnAK8C/WGv/HL+d9jywBhgCzAb+zRhzcZKX/CqxlqH+gBf49/jrjAH+BPwb0A9YADxvjPEaY/Lir/lK/Lx/AZ4wxhy4HToXaAUGAd+Ofx1JDXCmMeaOJDN3ylr7JlALnBPf1AxcDZQBnwOuM8Z8Ib7v3PifZfFbrn7AAHcDg4FxwDDgtuPJJCI9T0WXiKTCOcBzwDXW2n/Et50O9LPW3mGtDVprNwG/A65K8poPW2vft9a2AE8Bk+Lbvwy8YK191VobAn4GFABnAWcCxcA98ddcBPwD+Eq8k/4VwC3xVqt1wKOHe3FjTG9iBdzngIvbt0YZY7YZY05L8vs4oA7oDWCtrbLWrrXWRq21bxMrImce7kRr7cb499tmrd0F/OJIx4tIelKfLhFJhWuBJQfdfhsBDDbG1Lfb5gaWJXnNHe2WA8SKKYi19mw5sMNaGzXGbCXWmhYGtlpro+3O3RLf14/Ye97Wg/YdzpXAZmvtS/FbkUuNMQCPEPvAui7J7+OAIcBeAGPMGcA9wHhirXj5wNOHO9EY0x/4FbHitiT++vu6+Poi4jC1dIlIKlwLDDfG3N9u21bgI2ttWbuvEmvtZ4/zteqIFXQAmFglNAzYFt837KAnBYfH9+0iVpQNO2jf4Xjix2Ot3QtcBFwDvAzcZa21yQY2xpxOrOg6MDzFk8RaBodZa0uB3xK7hQjQ2XXvjm+fYK3tBXy93fEikiFUdInIkeQZY3ztvg7XOt4IXAKca4y5J77tTaDBGPNDY0yBMcZtjBkfL0COx1PA54wxs+N9uL4PtAHVwBvE+kv9hzEmzxgzC7gU+HN8OIq/ArcZYwqNMacQK6IOZwFwujHme/HXCcVfYwwQPcJ5CcaYXsaYzwN/Bv5orV0b31UC7LXWthpjphHrv3bArvj1T2i3rQRoIta5fgjwg2ReX0TSi4ouETmSBUBLu6/bDnegtbYeuBD4jDHmzniRcymxvlgfAbuB3wOlxxPIWvsesZaeB+PXvJTYEA1Ba20QuAz4THzfb4CrrbXvxk+/kdhtyh3EbhM+fITX+Sh+navj1/LHv49zgfuMMZccIebzxphGYq19PyLWB6v9cBHXA3fEj7mFWCF54HUDwH8Dr8WfwDwTuJ3Y0BX7gReIFY8ikmFMF1rIRUREROQYqaVLREREpAeo6BIRERHpASq6RERERHqAii4RERGRHqCiS0RERKQHpN2I9H379rUjR450OoaIiIjIUa1cuXK3tbZfMsemXdE1cuRIampqnI4hIiIiclTGmCNNJ9aBbi+KiIiI9AAVXSIiIiI9QEWXiIiISA9Iuz5dIiIicnxCoRC1tbW0trY6HSVr+Hw+hg4dSl5e3jFfQ0WXiIhIlqmtraWkpISRI0dijHE6Tsaz1rJnzx5qa2sZNWrUMV9HtxdFRESyTGtrK3369FHBlSLGGPr06XPcLYcqukRERLKQCq7USsXfp4ouERERSbni4uJDtv32t7/lscceA2DWrFldGpez/fGf/exnqa+vZ/PmzYwfPz41gdvpruuqT5eIiIj0iGuvvTYl11mwYAEA9fX1KbleT1FLl4iIiPSI2267jZ/97GcdtkWjUa655hp+/OMfA/DKK68wffp0pkyZwpVXXklTU9Mh1xk5ciS7d+8GIBKJ8N3vfpdTTz2Viy66iJaWFgBWr17NmWeeyYQJE/inf/on9u3bd8TtK1euZOLEiUyfPp25c+d2y/evli4REZEstuOny7vlugN/cPZxXyMcDvO1r32N8ePH86Mf/Yjdu3dz1113sXDhQoqKirj33nv5xS9+wS233HLYa3zwwQf86U9/4ne/+x1f+tKXmD9/Pl//+te5+uqrefDBB5k5cya33HILt99+Ow888MBht3/rW99KbP/BD35w3N9bZ9TSJSJZ6Yc//CHDhg1j5syZ+P1+p+OISCe+973vJQougNdff53169czY8YMJk2axKOPPsqWLUee2nDUqFFMmjQJgKlTp7J582b2799PfX09M2fOBOCaa65h6dKlSW//xje+0S3fr1q6RCSrVC9/jeuuu463160FYuMVnXP2OTz/n7/nzMnTcJf5cJf5qPloLcvXvEG/IQPZs3cPs2bNYvr06Q6nF0m9VLRIdZezzjqLxYsX8/3vfx+fz4e1lgsvvJA//elPSV8jPz8/sex2uxO3F7vCWtsjT3uq6BKRrBDeE6DqyRe48F+/dMi+SDTCsho/k3uNBqBmyzounXtdh2PcbjeDBg3iq1/9Kvfee2+PZBbJdd/5zndYunQpV155Jc8++yxnnnkmN9xwAxs3bmT06NEEAgFqa2sZM2ZMl65bWlpKeXk5y5Yt45xzzuHxxx9n5syZh91eVlZGaWkpy5cv5+yzz+aJJ57olu9XRZeIZKzq5a9R+dcXmVZ+MhN9I1m0aGGnx7ndbi657krKxp5MuL6VFb999pBjIpEItbW13HfffWzbto1TTz1VrV8ixyEQCDB06NDE+k033dTpcTfddBP79+/nG9/4Bk888QSPPPIIX/nKV2hrawPgrrvu6nLRBfDoo49y7bXXEggEOOGEE3j44YePuP3hhx/m29/+NoWFhVx88cVdfr1kGGttt1z4WFVUVNiujNshIrmnurqah+//Lb9/5nEAfHn5PHP9r/COLOei/3toS1d1dXWH4snv93PWWWcd9XUKCgqorKxU4SUZZ8OGDYwbN87pGFmns79XY8xKa21FMuerI72IZJTf/uo3zJgxI1FwAYSiYdb02cWF/3Zlp+ccXDQlW0QF24IsfrXy2MOKiLSjoktEMoK1lqonXuC6/3PDIfu8Xi/nXXB+Sl7na1/7WmI5z+1hYsNA2jbtTcm1RSS3qU+XiKS9aCDE/lc28sofn+t0fypvAf7xj39MdKJ99ke/Y1LhKPbNX89aTx0rmj/gvNnn63ajiBwTFV0ikrb8fj+Vz7/M5LahTB14MjPGToWXDj2uu4qgi//r6wRqtrH4jy9w6YOx6Us8t3uYO3cuc+bM6ZbXFJHspaJLRNJS+87uvrx8nv3R77j49qvhl9/rsQzGZSiaNpTVz+1IbAuHw9xwww2cdtppavESkS5Rny4RSTt+v59/ve7GxHooEuYtz1Y8pT5H8sz+3EUd1iORCLf8vx9rpHsR6RIVXSKSVg60cNWseSuxze1xc9755zmW6eAWLWstC5csYvb5s1V4iRyG2+1m0qRJnHrqqUycOJFf/OIXRKPRI56zefNmnnzyyR5K2PNUdIlIWvnf+397yLZvf/vbaXkrL9jWxsL5C5yOIZISfr+fu+++O2UfJAoKCli9ejXvvPMOr776KgsWLOD2228/4jkqukREesiiR57j908/dsj2q6++2oE0R5fnyWNy21Ba3t3ldBSR43Kghfk///M/mT079S24/fv3Z968efz617/GWsvmzZs555xzmDJlClOmTKG6uhqAm2++mWXLljFp0iTuv//+wx6XqdSRXkTSQuv7u1n4lxc63ZeOrVwAzz/wJOOb+1P5m7+yurCO2ZdfkrZZRY6kqqoqsRwMBqmqqkr5/+UTTjiBaDTKzp076d+/P6+++io+n48PPviAr3zlK9TU1HDPPffws5/9jH/84x9AbCqhzo7LVCq6RMRxoU+a2L/gfc46cbLTUbrkguu/yMJfPZ0YTuInv/qppg2SjDRr1qzEstfr7bCeSgemHgyFQtx4442sXr0at9vN+++/3+nxyR6XKXR7UUQcFWlqY99f12NDUc7+bGpGle9JbzZ9kFgOtgVZvHCRg2lEjk37Dwrd9cFh06ZNuN1u+vfvz/3338+AAQNYs2YNNTU1BIPBTs9J9rhMoaJLRBxjQxH2PbuBaFOQvCG9KL1otNORuuy88z59qjLP7WGqayQ2ah1MJHJ8uqPg2rVrF9deey033ngjxhj279/PoEGDcLlcPP7440QiEQBKSkpobGxMnHe44zKVii4RcUR1dTW3fucHvP7G67hL8yn/wjiMJ/Pektr/gnrmXx4ktL2R2+b8h4aSkJzX0tKSGDLiggsu4KKLLuLWW28F4Prrr+fRRx/lzDPP5P3336eoqAiACRMm4PF4mDhxIvfff/9hj8tU5sD91XRRUVFhM7mTnIgcXfvR5r3uPBY+/xLnfCZ2a9EYc8jx7d+njnf/kY47cOyBfcle98C+JX97hZlfiA2kWpDvo3LxIvXvEkds2LCBcePGdemczv7fS0ed/b0aY1ZaayuSOT/zPlaKSMa7+7/uTCwHIyGefO5pB9OkzmvrP/3AGAwGqfx7JxNFiqQpa60Krm6moktEetTyhUt4vvJFp2N0i/ZPfOV58phshxNtCzsXSETSioouEekx1lpe/sNfO2xzu91pO/hpV7W/lfjXmx/CNofUv0tEEjROl4j0mEBNHWf06dgf4je/+U1W9nsqnjGCcy+JPdn4s6f+h8pF6t8lPat9/0Q5fqm49aqWLhHpEcteqOS/b7vzkO1z5sxxIE33W/7Wp61bwWCQRS8tdDCN5Bqfz8eePXvURytFrLXs2bMHn893XNdRS5eIdLt5//MQ37s+Nmp7Qf7xvWllig79u9x5VLhHYaMW41LLg3S/oUOHUltby65dmhc0VXw+H0OHDj2ua6joEpFu5ff7EwUXQFsos0eUTlbH8bt+zcSCkQRqtlE07fjetEWSkZeXx6hRo5yOIQfR7UUR6VYP/3peh3WXK/fedi644YsALH78Be760W3qWC+So3Lv3U9EesxDv/kffvfkIx22zZ0715kwDso/oTfrvNu59NfX8l8/uZ2ZM2cyb968o58oIllFRZeIdAu/38+1N1x/yPZs7Th/NG82vp9YDoVC3HDDDWrxEskxKrpEpFtUPv+y0xHSynmzz++wHo1GqaqqciaMiDhCRZeIpJyNWgp3RJyOkVYOHqMr3+Nl5syZDqURESeo6BKRlFv82PN8/+G7nI6R1p767v1M6TfW6Rgi0oNUdIlISkWaglQ+k51zK6ZSxYjxNFRuIhpUi6BIrkiq6DLGXGKMec8Ys9EYc3Mn+28yxqw3xrxtjKk0xoxoty9ijFkd/3ouleFFJL34/X5uv+Fmyn0lTkdJe54BxUQb22iq/tjpKCLSQ446OKoxxg3MBS4EaoEVxpjnrLXr2x22Cqiw1gaMMdcB9wFfju9rsdZOSnFuEUkzfr+fs846CwBfXr7DadJf6UUnsufxNQRqtlFwan/y+hU5HUlEulkyLV3TgI3W2k3W2iDwZ+Dy9gdYaxdbawPx1dcBDbkskmMefeTRxHIoEnYwSWbIG1hC4aRBYKHhlY2aI08kByRTdA0BtrZbr41vO5zvAO07dPiMMTXGmNeNMV84howikub8fj8PzXsose7J0wxjySg+ZwSuwjxCdY20rNvpdBwR6WbJFF2dzc7a6UcyY8zXgQrgp+02D7fWVgBfBR4wxpzYyXlz4oVZjSbnFMk8i16p7LD+rW99y6EkmcXl81AyKzY/XtOyLepUL5Llkim6aoFh7daHAnUHH2SMuQD4EXCZtbbtwHZrbV38z01AFTD54HOttfOstRXW2op+/fp16RsQEedNK+049MHVV1/tUJLM4zulX6xTfXOQQM02p+OISDdKpuhaAZxkjBlljPECVwEdnkI0xkwGHiJWcO1st73cGJMfX+4LzADad8AXkQwX3hPgtPCgDtsOHghUDs8YQ8mskQA0v7mNSHPQ2UAi0m2OWnRZa8PAjcDLwAbgKWvtO8aYO4wxl8UP+ylQDDx90NAQ44AaY8waYDFwz0FPPYpIhmtcsvkwHQ4kWfnDy8g/sTc2FKHpNQ0hIZKtkurtaq1dACw4aNst7ZYvOMx51cBpxxNQRNLX0r+9yku/e4YZY6c6HSXjlcwcSdumvbS8vYOiqYPx9Cl0OpKIpJgeMRKRYzJv3jy+973vAeB5VW8lx8vTp5CCCQNpWbODxiWbKf/iKU5HEpEU0zRAItJlfr8/UXABhDUuV0oUnzUck+fmtUVLueuHt+L3+52OJCIppKJLRLps8aLFTkfISu5iL+sKd3Dp3Ov4r/vu4Oyzz+aHP/yh07FEJEVUdIlIl5058FSnI2StN3Z++qxRNBrlvvvuY968eQ4mEpFUUdElIl0SbQlxanP/DtseeuihwxwtXXXe7PMP2faHP/zBgSQikmoqukSkS5pr6rBtHUdOnzNnjkNpsk9nY5z5fD4HkohIqqnoEpGkRQMhAisPmZBCutm4k8c5HUFEUkBFl4gkrXnFNmwogndUudNRcobXnceXpn/W6RgikgIqukQkKZHmIIFVsVau4hnDHU6TO/77C//Go394lGvnfE9DSIhkOI1oKCJJibVyRck/oRzvoBKn4+SMH8z/aWxhOTz86CNUVVVpbkuRDKWWLhE5qkhTkMCq7QAUzxjhcJrcFQqFqKqqcjqGiBwjFV0iclTNb9ZCOEr+6N7kDSx2Ok7OstbSp08fp2OIyDFS0SUiR7S8cgn3/Pw+araso/gs9eVyksEw/+ln1LdLJEOpT5eIHJbf7+ecC2YBUOD1Ufm1iUwfoP5ETrFYXln4KsteW05lZaX6dolkGLV0ichhLXp5YWI5GFF/onQRDAb1byGSgVR0ichhnVF+cmLZ6/Uya9Ys58JIgteTp38LkQykoktEOrW8cimLXqlMrOt2Vvp4+oZfccbUaU7HEJEuUp8uETlErC/XzA7bVHClD9sW4c6bfsxFX7tM/y4iGURFl4gcYvHCRU5HkCO4dO51ANz7v79UC6RIBtHtRRE5xBkDTnE6giRBHepFMouKLhHpIBqMMD44yOkYkoQ8l4eZ5848+oEikhZUdIlIBy1rd2Bbwk7HkCQ8PecBpvQd43QMEUmSii4RSbCRKM0rtjkdQ5JUMWI8Ta9vxVrrdBQRSYKKLhFJaFm/k2hjEE+fQqejSBJchXmEdzYT3FzvdBQRSYKKLhEBwEYtzW/EWrmKzhjqcBpJRmHFEACa3tjqcBIRSYaKLhEBoO2DPUT2teAuzcc3rp/TcSQJhZMGYvLd+JdXc9f/u1UTYYukORVdIoK1lqbXY60lRacPxbiMw4kkGa58D2u927l07nX81z13MHPmTObNm+d0LBE5DBVdIkJwcz3hnc24CvMoOG2A03GkC978ZH1iORQKcf3116vFSyRNqegSkcQTi4VTB2M8elvIJOddNLvDeiQS4bHHHnMojYgcid5dRXJcaGcTwS31mDwXhZM0KGqm0RRAIplDRZdIjjvQylVw2gBcPk3Hmuny87xcffXVTscQkU6o6BLJYZHGNlrf3Q0GCqcOcTqOpMD8Gx7kjKnTnI4hIp1Q0SWSwwJvbYeoxTemL54yn9NxJAWmDj6FlrWfOB1DRDqhokskR0WDYQJrtgOfDrIp2SGwchs2qqmBRNKNii6RHFX1xAJ++eIjrA58hHdwidNxJEXcZT4i+9to+2CP01FE5CAqukRyUPVr1cz+9uXc/dI8vnjPHI3rlEWK4q2WzStqNRG2SJpR0SWSgxbOX5BYDoZCVFVVORdGUqpgfH+Mz0NoexOhbQ1OxxGRdlR0ieQYay3Tik9KrHu9XmbNmuVcIEkpk+emcPIgaras466bb1Mrpkga0aA8Ijlm2XMLWbqiOrFeWVmpATazzJq2LVw69zoAfvbXh/RvLJImVHSJ5BC/38/ML1zUYZt+GWefZW++llgOtgWpqqrSv7NIGtDtRZEcsujFV52OID2g/e3iPLeHc6ef7VwYEUlQ0SWSQ87od4rTEaQHtG/VuvOyf+XVp15Q3y6RNKDbiyI54rUly1j8cqXTMaSH/WD+TwG475FfqW+XiMNUdInkAL/fz9mzznU6hjgoGFTfLhGn6faiSA5YvHix0xHEYXluj4YGEXGYii6RHFBqC5yOIA57+rsPUDFqvNMxRHJaUkWXMeYSY8x7xpiNxpibO9l/kzFmvTHmbWNMpTFmRLt91xhjPoh/XZPK8CKSnB3rNjsdQRxWMWI8zSvrnI4hktOOWnQZY9zAXOAzwCnAV4wxBz8CtQqosNZOAJ4B7ouf2xu4FTgDmAbcaowpT118ETma0M5mzugzzukY4jQDre/uJtIUdDqJSM5KpqVrGrDRWrvJWhsE/gxc3v4Aa+1ia20gvvo6MDS+fDHwqrV2r7V2H/AqcElqootIMgJv1VExQreVcl3+6D4QtQRWb3c6ikjOSqboGgJsbbdeG992ON8BXjzGc0UkhaKBEC3rdzodQ9JA0dTBALSs2YENRx1OI5Kbkim6TCfbbKcHGvN1oAL4aVfONcbMMcbUGGNqdu3alUQkEUlGYM0OiFjyT9Bd/VyXN7QXnv5FsUL8Xb3PijghmaKrFhjWbn0ocEhvTGPMBcCPgMustW1dOddaO89aW2GtrejXr1+y2UXkCGwkSmBV7FZSYbyVQ3KXMSbR2hWoqcPaTj87i0g3SqboWgGcZIwZZYzxAlcBz7U/wBgzGXiIWMHV/l7Gy8BFxpjyeAf6i+LbRKSbtb63m2hzEE+fQrwjypyOI2nAd3I/XIV5hHc1E6ptcDqOSM45atFlrQ0DNxIrljYAT1lr3zHG3GGMuSx+2E+BYuBpY8xqY8xz8XP3AncSK9xWAHfEt4lINwu8FW/lmjIIYzq70y+5xnhcFEwcSM2Wddz1n7dpPkaRHmbSrYm5oqLC1tTUOB1DJKOFdjSy5/E1mHw3/a6dhsvrPmzhdbj3gIOPP/i49vuttUdcP/g6Bx97pNc9+Jjj3X+k4w7Ol+x12x9/uOWDr3Gk995ki+QjXfdw/37LK5dwzgWzACgoKNB8jCLHyRiz0lpbkcyxGpFeJAs1x1u5Ck4bgMvrdjiNpJNlb1YnloNtsfkYRaRnqOgSyTKR5iCt8afTCierA7101H7+xTy3h5lnayJ0kZ6ioksky7Ss/SQ2TMSJ5XjKfE7HkTTT/lbi03MeYHL5aAfTiOQWFV0iWcS2G3FcrVxyNBUjxtP81nYNHyHSQ1R0iWSR39x1P1/6xY08sfZFvCM1TIQcmSnwEP6kiVBdo9NRRHKCx+kAIpIa8+bN48Zbvw/Akg9WUHL2CObMmeNwKklnhRMG0vxGLYFV2/EO6eV0HJGsp5YukSzx9J+e6rA+f/58h5JIpiicNBBMbCDdSFPQ6TgiWU9Fl0iW+PyEWR3Wr7jiCmeCSMZw9/KRP7oPRC0ta3Y4HUck66noEskC0ZYQXx55fmL9oYce0q1FSUrhlEEABNZsx0aiDqcRyW4qukSyQMvaTyD86S9MFVySLO+wUjx9C4k2h2h9f4/TcUSymooukQzXfpgIka4yxlA4Od7a9Vadw2lEspuKLpEM1/bRPiL723CX5jsdRTKU75T+mHw3obpGQp80OR1HJGup6BLJcAdaJwonDXI4iWQql9dNwfgBAARWqdVUpLuo6BLJYOF9LQQ314PHRcFpA5yOIxnswC3Glg27iLaEHE4jkp1UdIlksAOtEgXj+uEqyHM4jWQyT3kB3lHlEI7GHswQkZRT0SWSoaLBCC3rYr8cD7RSiByPRIf61duxUc3HKJJqKrpEMlTrhl3Ytgh5g0vIG1DsdBzJAvmjynGX5vPG2yu56z9uwe/3Ox1JJKto7kWRDGSt/bQDvVq5JEWMy7DWXcelc68D4O7f/JzKykqmT5/ucDKR7KCWLpEMFNrWQHh3AFdhHr4xfZ2OI1nk9W1rE8vBYJCqqirnwohkGRVdIhko8Fa8A/2EgRiPfowldc67cHZi2evJY9asWc6FEckyercWyTCRpjZaP9gDBgonDnQ6jmSZ9rcSn77ul5wxdZqDaUSyi4oukQwTWLMDopb80X1w99Io9NJ9pg4+hdYNu5yOIZI1VHSJZBAbidKyZgcAhVPUgV66X2DVdqzV8BEiqaCiSySDtH6wh2hzCE+fQrzDSp2OI1nOVZhHeFczoW0NTkcRyQoqukQySGIE+kkDMcY4nEayXcGEWJ9BzccokhoqukQyRGhnM6HaBozXTcH4/k7HkRxQOHEgGGh9fw+RpqDTcUQynooukQyRaOU6tT8ur8Y1lu7n7pVP/ug+ELW0vL3D6TgiGU9Fl0gGiLaGad2wE9AI9NKzCicPombLOu6+9x6ql7/mdByRjKaPyyIZoGXdJ9hQFO/wUjx9Cp2OIzlk5bb1iWmB7n/lYSoXLdK0QCLHSC1dImnOWktgdezWolq5pKctWbIksRwMhjQtkMhxUNElkuaCm+uJ7GvFVRLvXyPSg9pPA+Q2Lja/9yF+v9/hMWhpAAAgAElEQVS5QCIZTEWXSJpLtHJNHIhxaZgI6VntbyUGIyHmPfoHZs+ercJL5Bio6BJJY+H9rbRt3AtuQ8GEAU7HEQEgGAzqNqPIMVDRJZLGWlbHHtP3jemLu8jrcBqRGK8nr8NtRxFJjooukTRlw1ECa+PzLKoDvaSRZ276DWeeeabTMUQyjooukTS15M8v8ssXHmZV40byBpc4HUckYUr5SQS37Hc6hkjG0ThdImnI7/dz3jc+D8ADVY9R+ZmxGhtJ0kpgVR35I8ucjiGSUdTSJZKGHvmf3yeWgyGNjSRpxmVo+3AvkYZWp5OIZBQVXSJpxu/3M+/x/02sezwedVqWtOIb2xcsBFZrPkaRrlDRJZJmFr1S2WH9W9/6lm4tSlo58GBH4O0d2HDU4TQimUNFl0iaOXPw+A7rV199tUNJRDqXN7gEz4AibEuY1nd3OR1HJGOo6BJJIzZqOS0yuMM2tXJJujHGUDg59v+0+a3tWGsdTiSSGVR0iaSRJU+/xAN//4PTMUSOquDkvpgCD+FPmghtb3I6jkhG0JARImnC7/dz3lWfdTqGSFJMnpvC0wbS/GYtgVV1eAePdTqSSNpTS5dImqh84RWnI4h0ScGkgWCg9d3dRJqDTscRSXsqukTSxBl9T3E6gkiXeEp95J/YG6KWlrc1fITI0SRVdBljLjHGvGeM2WiMubmT/ecaY94yxoSNMf980L6IMWZ1/Ou5VAUXySbR1jATGOJ0DJEuO9ChPrB6Bzai4SNEjuSofbqMMW5gLnAhUAusMMY8Z61d3+6wj4FvAv/eySVarLWTUpBVJGu1vPMJNqRfWJJ5vCNKcfcuILK3hbaNe2MDp4pIp5Jp6ZoGbLTWbrLWBoE/A5e3P8Bau9la+zag3xoiXWStJbBqu9MxRI5JbPiI2GCpzavqHE4jkt6SKbqGAFvbrdfGtyXLZ4ypMca8boz5QpfSieSA4OZ6IvtacfXKdzqKyDEpGN8f43UT2tpAaGez03FE0lYyRZfpZFtXRsIbbq2tAL4KPGCMOfGQFzBmTrwwq9m1S6MbS24JvBVrHSicNNDhJCLHxuX1UHBqfwACau0SOaxkiq5aYFi79aFA0j9V1tq6+J+bgCpgcifHzLPWVlhrK/r165fspUUyXnhfC22b9oHbUHiaii7JXAduMbas30W0JeRwGpH0lEzRtQI4yRgzyhjjBa4CknoK0RhTbozJjy/3BWYA6498lkjuONCXq+CU/rgK8xxOI3LsPH0K8Y4sg3CUlrWfOB1HJC0dteiy1oaBG4GXgQ3AU9bad4wxdxhjLgMwxpxujKkFrgQeMsa8Ez99HFBjjFkDLAbuOeipR5GcFQ2GE7+cDrQSiGSyoqnx+RhXbcdGNR+jyMGSmgbIWrsAWHDQtlvaLa8gdtvx4POqgdOOM6NIVqp6cgELX3qBc888m0sGnO10HJHj5h1VjrvMR6S+lbaNe/CN0fARIu1pRHoRB1RXVzP7W5dz90vz+OI9c/D7/U5HEjluxhgKp8Rbu95Sh3qRg6noEnFA5bMvJpaDoRBVVVXOhRFJoY7DRzQ5HUckrajoEnFAcYM7sez1epk1a5ZzYURSyJXvoWB8fPiItzTor0h7KrpEetiylxdz07zbE+sPPPAA06dPdzCRSGoduMXYsmEX0YCGjxA5QEWXSA9b+NQLHdb37NnjUBKR7uEpLyD/hHIIRwms3eF0HJG0oaJLpAdF28JMKxnTYZtuLUo2OtDaFdDwESIJKrpEelDL2k+YOuSUDtt0a1GykXdkGe7eBUQbg7R9oNZcEVDRJdJjbNSqY7HkDGMMhZMHUbNlHf99y50aFkWEJAdHFZHj1/bRPiL7W3H1ync6ikiPWN20iUvnXgfAz5/7HZWLKtWyKzlNLV0iPWTJkwv41aLHWefRoJGSG5a+tiyxHAwGNR6d5Dy1dIn0gGUvVHLJLd8A4IEljzucRqRntH9IJM/t4ZxpM5wLI5IG1NIl0gMWPv3p1KXBYNDBJCI9p/2txKfnPMBE7wgH04g4T0WXSDeLNAU5vdenw0R4vV4H04g4o2LEeAJrtmNDEaejiDhGRZdINwus3k7FsFMT65WVlQ6mEXFG3sBibEuYlvW7nI4i4hgVXSLdyIajtKzpOCK3nt6SXFRYMQSA5pXbsFaDpUpuUtEl0o1a1u8kGgjh6V/kdBQRR/nG9MFV4uWNt1Zw13/cqnG7JCep6BLpJtXV1fzkjruo2bKOovinfJFcZdwu1nq2c+nc67jlZ3cye/ZsFV6SczRkhEg38Pv9zJgRezzel5dP5VWnOZxIxHlvbF+XWD4wbpdut0suUUuXSDdoPwhkKBJmybKlzoURSRPnXXB+YtnrydNk75JzVHSJdIOzJ5+ZWPbme/XLRYSOD5E8NecBzphY4WAakZ6nokukG0wwwxLLlZWab07kYBXDTiWwWhPAS25R0SWSYtFAiJb1OxPrKrhEOhdYtR0bjjodQ6THqOgSSbHmVXWgXyQiR+TpXxT7gPLOzqMfLJIlVHSJpJANRQis0i0TkaMpOj0+WGqNBkuV3KGiSySFWtbtxLaEyRtY7HQUkbTmG9sXV698IntbaNu41+k4Ij1CRZdIitiopblmGwBF04Y6nEYkvRm3i6KpgwFoXrHN4TQiPUNFl0iKtH2wh0h9K+5SH/kn9XE6jkjaK5gwEJPvJrStgeC2BqfjiHQ7FV0iKWCtperJF/jVosdZ663DuIzTkUTSnsvrpnDSIECtXZIbVHSJpMCyvy/kM3d8k7tfmsfnr79Kc8qJJKlwymBwG/7wx4e56PwLmDdvntORRLqN5l4USYGFz7yQWD4wp5yIHJ272MtftizmB/N/CsCriysBmDNnjpOxRLqFWrpEjlN4d4BpZScn1r1eTfsj0hUvrK3qsD5//nxngoh0MxVdIsepuWYbFSPGJ9Y17Y9I1/zzVV/qsD5p0iSHkoh0LxVdIsch0hTsMOUPaNofka46+Fbigw8+qH6RkpVUdIkch8DKOohYDREhkkKtra089thjTscQSTkVXSLHKNoaJrA6NuWPBkMVSR1rLQ8//LBauyTrqOgSOUaBVduxwQje4aV4B5c4HUckq4TDYT0FLFlHRZfIMYgGIzSvjE/5c+Ywh9OIZJ88l4eZ5850OoZISqnoEjkGLWt3xCa2HlSCd3ip03FEss7Tcx5gSt8xTscQSSkVXSJdZCNRmt880Mo1FGM05Y9IqlWMGE/T61ux1jodRSRlVHSJdFHLOzuJNgXx9C0k/8TeTscRyUquwjzCO5sJfrTP6SgiKaOiS6QLbNTS/EYtAEVnqJVLpLsUnT4EgKbXax1OIpI6KrpEuqD1/d1E6ltxl/rwndzP6TgiWatg0kBMvpvQtgaCtfudjiOSEiq6RJJkraX59a1AvJXLpVYuke7i8noonDIYUGuXZA8VXSJJatu0j/CuAK4iLwWn9nc6jkjWK5o6GJPnIvjRPkKfNDkdR+S4qegSSUKHVq7Th2A8+tER6W6ugjwKJg4CoCn+8yeSyfSbQyQJwY/3E6prxPg8FEwc6HQckZxRdPpgcBva3t9DaHez03FEjktSRZcx5hJjzHvGmI3GmJs72X+uMeYtY0zYGPPPB+27xhjzQfzrmlQFF+kp1lqaqj8GoKhiCC6v2+FEIrnDXZxP4WmxDzrNfrV2SWY7atFljHEDc4HPAKcAXzHGnHLQYR8D3wSePOjc3sCtwBnANOBWY0z58ccW6TnBj/cTqm3A+DwUThnkdByRnFN05lBwG1rf3a3WLsloybR0TQM2Wms3WWuDwJ+By9sfYK3dbK19G4gedO7FwKvW2r3W2n3Aq8AlKcgt0iOqq6u56z9upWbLOopOH4Ir3+N0JJGc4y7Jp3BCvLWrWq1dkrmS+Q0yBGj/v7yWWMtVMjo7d0iS54o4yu/3M2PGDAB8efksvPJUZqDJrUWcUHTGUAJv76D1vd2EdjWT16/I6UgiXZZMS1dngxElOxlWUucaY+YYY2qMMTW7du1K8tIi3Wvx4sWJ5VAkzNLq5Q6mEcltHVq71LdLMlQyRVctdPh4PxSoS/L6SZ1rrZ1nra2w1lb066dRviU9zDi5IrHszfcya9Ys58KICEVnxPt2xVu7RDJNMkXXCuAkY8woY4wXuAp4LsnrvwxcZIwpj3egvyi+TSStWWsZ3/rp0BCVlZVMnz7dwUQiotYuyXRHLbqstWHgRmLF0gbgKWvtO8aYO4wxlwEYY043xtQCVwIPGWPeiZ+7F7iTWOG2Argjvk0krQW37Ce0rSGxroJLJD0caO1a/koVd/34dvx+v9ORRJKW1KNY1toFwIKDtt3SbnkFsVuHnZ37v8D/HkdGkR5lraXptS1OxxCRTrhL8lmXV8elc68D4Ce/uFct0ZIxNCK9yEGCW+pjo88XaHgIkXT0xp53E8vBYJCqqirnwoh0gYoukXastTQtj48+f7pGNxFJR+dffEFi2evO00MukjFUdIm007ZxL6HtjbgK8yicPNjpOCLSifa3Ep/67v1MHTnewTQiyVPRJRJno5amZbG+XEXTh2mORZEMcedNP6a6utrpGCJHpU4rInGt63cS3hPAXfrpY+kikt4OdKj/+d9+R+UidaiX9KaWLhHAhqM0vhbry1U8YwTGox8NkUwSDAY7zCIhko70m0UECKzZTrShDU/fQnzjNCuCSKbJc3s464TJTscQOSIVXZLzosEwTfHRrYvPGYFxdTZlqIiks6fnPMCpTf2x0WSnBhbpeSq6JOcFVtRhW8LkDS4h/8TeTscRkWNwxsQKIvtaaFn3idNRRA5LRZfktGggRPOKbQCUzByJMWrlEslExWePAKDptY+xoYjDaUQ6p6JLclrT61uxoQj5J5TjHVrqdBwROUa+k/vi6V9EtClIYNV2p+OIdEpFl+SsyP5WAqtjb87F54xwOI2IHA9jDCXnjgSg6s8v8pM7/1uTYUva0ThdkrMal22BiMU3rh95/YudjiMix8k7sozVLR/x+fvnAFBwd4Emw5a0opYuyUnBugZaN+wCj4sStXKJZAVjDDWtHybWNRm2pBsVXZJzrLU0LvoIgKKKwbhLfQ4nEpFUmX3pxYllTYYt6UZFl+Sc1nd3xya1Lsqj6IyhTscRkRQ6ZDLsoac4mEakIxVdklNsKELjks1A7BFzl1fdGkWyVcWI8TQs3oS1GjBV0oOKLskpix5+ngf+/ntW7d9IwfgBTscRkW7kKvYS/qSZlnd2Oh1FBFDRJTlkeeVSLvjeP3H3S/O44ufX8/obrzsdSUS60YEhJJqWbiEa1ICp4jwVXZIzXv3j3xPLwZCeahLJdr5T+uEZWEy0OUjzm7VOxxFR0SW5IfRJE6cXnZRY93q9eqpJJMsZY+h13igAmldsI9LQ5nAiyXUquiTrWWtpWLyJihHjE9s0YKJIbvAOLcU3ti+EozQu3ex0HMlxKrok67Vu2EVoawOm4NMnFVVwieSO4nNHgsdF64ZdtH1c73QcyWEquiSrvbZkGXfdfBs1W9ZRMnOk03FExAGeMh/F8TH5GhZ+iI1EHU4kuUqDFEnW8vv9nD3rXAB83nwqvzrB4UQi4pSiaUNpeWcnkT0tNNfUJYowkZ6kli7JWpX/eDmxHIqEWbJkiYNpRMRJxuOi1wUnAtDs/5hIQ6vDiSQXqeiSrGStpcJzQmJdTyuKSP6ocnxj+2JDURri86+K9CQVXZKVWt7+hEnFnxZdelpRRABKzhuFyXPR9sEe2jbtdTqO5BgVXZJ1ooHQIY+Gq+ASEQB3ST7FM4YD0LBwEzakkeql56jokqzTuHQztjWMd3ip01FEJA0VThmMp28hkf2tNL2hkeql56jokqwS3NZAy9pPwGXodeGJTscRkTRk3K7E+0Pzm7WE97Y4nEhyhYouyRo2HOXB/3cvV/3uJp7asRRP70KnI4lImvIOLaVgfH+IWPa//AHWWqcjSQ7QOF2SNR68+V7+/fGfALDkZyvwndTH4UQiks5KZo2ibdM+QrUNBFZtp2jKYKcjSZZTS5dkhdCOJv76/LMdts2fP9+hNCKSCVwFeYnbjE1LNxPer7G7pHup6JKMZyNR9r/4Pp8/bVaH7VdccYUzgUQkY/jG9P107K6XN+o2o3QrFV2S8ZpfryW8O8A1F38pse2hhx5izpw5DqYSkUxRMvsETIGH4Jb62IM4It1ERZdktNDOZppe3wpAr4tHJ7ar4BKRZLmLvPSaHbvN2Lj4IyKNbQ4nkmyloksyVvVr1dx+w83UfLSWwkmDyB9e5nQkEclQvpP7kj+6NzYYoeEV3WaU7qGnFyUj+f1+Zpw9AwBfXj4Lv/QqM9C4XCJybIyJje23e+t+2jbto3X9LgpO7e90LMkyaumSjFT5j1cSy6FomKWvLXcwjYhkA3dxPiXnx+Zsbaj8kEiDbjNKaqnokoxjQxGm2hGJda/Xy6xZs5wLJCJZo+DU/uSfWI5ti1D/wnvYqG4zSuqo6JKM01i1mcnln3aar6ys1ITWIpISxhhKLxmDqyiPUG0DzW9qbkZJHRVdklFaN+4hsHo7uEximwouEUklV2EepZ8ZA0DT8i0E6xodTiTZQkWXZIxIU5D9L30AQMm5I45ytIjIscsfVU7h1MFgYf8/3iMaDDsdSbKAii7JCNZa9r/4PrYljHdEGYUVQ5yOJCJZruTckXj6FRHZ30rDwk1Ox5EsoKJLMkKgpo7g5npMgYfSz56EMeboJ4mIHAfjcVF26VjwuGh9ZyctG3Y5HUkyXFJFlzHmEmPMe8aYjcaYmzvZn2+M+Ut8/xvGmJHx7SONMS3GmNXxr9+mNr7kgtAnTTQu2wxA6SUn4S7OdzaQiOQMT59Cep03CoCGVzdqUmw5LkcdHNUY4wbmAhcCtcAKY8xz1tr17Q77DrDPWjvaGHMVcC/w5fi+D621k1KcW3JEtDVM/XPvQsRSMHEgvtF9nI4kIjmmYOJA2j7aR9vGvSz8+V94K38b551/nh7ikS5LpqVrGrDRWrvJWhsE/gxcftAxlwOPxpefAWYb3f+R42StZf+C94jUt+LpX5T4tCki0pNiw0icxFu73+Mzd36TH/34R5x77rnMmzfP6WiSYZIpuoYAW9ut18a3dXqMtTYM7AcONEmMMsasMsYsMcacc5x5JYc0v76Vtg/3YXweyi4fh8lzOx1JRHKUqyCPtzyf/ioMh8PceOON+P1+B1NJpkmm6OqsxergIXoPd8x2YLi1djJwE/CkMabXIS9gzBxjTI0xpmbXLnVUFGj7aB9Nyz8GoOzzY/GU+RxOJCK5bvbnL+6wHolEqKqqciaMZKRkiq5aYFi79aFA3eGOMcZ4gFJgr7W2zVq7B8BauxL4EBhz8AtYa+dZayustRX9+vXr+nchWSVc30r9P94DoHjGcPJHlTucSETk0IGYXbjoXVLmUBrJRMkUXSuAk4wxo4wxXuAq4LmDjnkOuCa+/M/AImutNcb0i3fExxhzAnASoMFO5LBsKEL93zdgW8Pkn1hO0fRhRz9JRMQB4WiY//N//41rr71WtxklKUctuuJ9tG4EXgY2AE9Za98xxtxhjLksftgfgD7GmI3EbiMeGFbiXOBtY8waYh3sr7XW7k31NyHZwVpLw8IPCe9sxl3mo/SzYzUel4iktbZwkIceeojZs2er8JKjOuqQEQDW2gXAgoO23dJuuRW4spPz5gPzjzOj5IhAzTZa1u0Ej4uyy8fh8iX131NExHHBYJCqqioNIyFHpBHpJS20vrebxqrNAJR9dgx5/YucDSQi0gV5Lg8zxlU4HUPSnIoucVywroH6Be8DUHzuSHxj+zqcSESka56e8wBjt5YQ3htwOoqkMRVd4qhwfSv1z26AcJSCCQMpmqaJrEUk88yYfS62Ncy+Z9YTDYScjiNpSkWXOCbaEmLf/HeIBkJ4R5bR68IT1XFeRDJS6efG4hlQTGR/K/v+uh4bijgdSdKQii5xhA1H2fe3DUT2tuDpV0jZZSdjXCq4RCQzubxuyr94Cq5e+YS2N1K/4H2sPXgcccl1Krqkx1W/Vs2t3/53/K9V4yr2Un7Fqbjy9aSiiGQ2d7GX8itOwXjdtL2/h8bFH6nwkg70m056VHV1NTPOngGALy+fV559gXNK8h1OJSKSGnl9iyi7fBz75r9DYGUdxuum5OwRTseSNKGWLukx1lpe/sNfE+uhaJjlb7/pYCIRkdTLH1lG2aVjwUCzfytNr289+kmSE1R0SY+w1tK0dDOnF4xObPN6vcyaNcu5UCIi3cQ3pi+ln41NNdy0bAvNNdscTiTpQEWX9Ihm/1aa39xGxajTEtsqKys1erOIZK2CU/rT6+LYB83FjzzP7f/6n5oqKMep6JJu1/xmLU2vfQwGyj4/NrFdBZeIZLvCCQNZ328vl869jtsevJvZ55+vwiuHqeiSbtX8Zi2NSzYDUHrJSRptXkRyzhvb30ksB9uCvPqXfziYRpykpxelW1RXV/PKI3+jwnsCFSPG0+ui0RSMH+B0LBGRHte+72qeJ48p4WE0r9xG0VTNwJFrVHRJylW/1nFYiBcfmc+siQMdTiUi4oz2XSle+J+nOGVvbxoXfUS0NUzxWcM1E0cO0e1FSSkbifLSQ88k1kPRMP4tbzuYSEQkfZz/ncvodcno2HAS1VtpXLRJA6jmEBVdkjI2FGHfsxuY1mtMYpuGhRAR6ajwtIGUXXoyuA2Bt7az/8UPsFEVXrlAtxclJSLNQer/toFQXSOnnzwpsV3DQoiIHMo3ti/l+W7q/7aB1nd2YltClF46FpdXv5azmVq65LiFdjWz549rCNU14irx0ucrExL7VHCJiHQuf2Q55VeOx/g8tG3ax94n3ia8v9XpWNKNVHTJcWnduIe9T7xNtKGNvEHF9Pn6JDx9Cp2OJSKSEbxDetHn6xNx9y4gvDvAnsdXE6zd73Qs6SYquuSYWGtpemMr9c9uwIYi+Mb1o/eXT8Nd7HU6mohIRvGUF9DnaxPxjizDtoTZ+5d1BNbucDqWdAMVXdJlNhxl/4L3aVq6BYDic0ZQ+rkxmDy3w8lERDKTy+eh/IpTKZw6GKKWhpc20rB4kzrYZxn12JOk+f1+Fi14hSl2BJNLR2PyXJR+biy+k/o4HU1EJOMZl6HX+Sfg6VNIw8IPCdTUEdrRRNnnx+IuyXc6nqSAii5JSnV1NTNmfDrg6fybfsOFN32JvP7FDicTEckuhRMH4uldQP3z7xKqbWD3I6so/cxJ+EbrA26m0+1FOapoMNxxwNNImNVF21VwiYh0E++wUvpcMxnvqHJsa5j6ZzfQsGgTNhx1OpocBxVdckShT5rY8+hqppW0G/A038t5F5zvYCoRkeznLvJSfsUplMwcCS5DYGUde55YQ3hfi9PR5Bip6JJO2XCUxuVb2PPHNUTqWzmz4ozEPg14KiLSM4wxFE0bSu+vTsBd6iO8s5k9j66iuWabOtlnIBVdcojgtgZ2P7aKZv9WiFoKpwyiz9cnJvar4BIR6VneQSX0uWYSvnH9sKEojYs/Yu+TawjtanY6mnSBOtJLQjQYpmnpFgKrtgPg7l1A6cWj8Q4tdTiZiIi48j2UfX4sreP60fDqRkLbm9jz2GqKzhhK8ZnDMB61o6Q7FV2CtZa2jXtpqNxEtLENXLHm7OLp+iEWEUk3vhN74x06hcalm2lZvYNm/1Za399NrwtOJH94mdPx5AhUdOW40CdNNFZ9RPDj2LQTngHFlF4yWk8mioikMVe+h9ILR1Mwrh/7X9pIZE8L+/6yjvzRvSmZNQpPeYHTEaUTKrpyVKQpSNPyLbSs/QQA4/NQfNZwCicPwriMw+lERCQZ3qGl9P3mZJrfrKX5zVraNu6lbdM+CicPovis4bh8+jWfTvSvkWOiwTCBlXU0v1GLDUXBZWI/nNOH4SrIczqeiIh0kfG4KD5rOAUTBtC0bAst63YSWFlHyzs7KT5rGIUTB6mrSJpQ0ZUjom1hAm/V0VxTh20NA6gZWkQki7iL8yn9zBgKpwymcfFHBLfup3HRRzS/sY2iaUMonDhQc+Q6TEVXlou2hGheWUfgrTpsWwSAvCG9KD57uDpciohkobwBxZR/eTxtH+6lafnHhHc107j4I5rfqKXo9CEUTBqEy6viywkqurJUuL6VwOrttKzegQ3Fii3v8FKKpg/DO6wUY9RvS0QkWxlj8I3uQ/6JvWPFl38r4R1NNC7ZTNObtRROHkThxEG4i71OR80pKrqyiLWW4OZ6Aqu20/bh3sR278gyiqcP03hbIiI5pn3xFdxcT1P1x4TqGmmu3krz67X4xvShcMpg8gaX6MN4D1DRlQWiLSFaNuwisGo7kb3xObncBt/YvhRNGUzeoBJnA4qIiKOMMeSPKsc7sozg1v0E3tpO28Y9tL67m9Z3d+PpX0Th5EH4Tu6Ly6vSoLvobzZD2XCUtk17aVm/K9aqFZ+Dy1XspXDSIAomDMBdpGZjERH5lDGG/OFl5A8vI9LQSmD1DgJv7yC8s5mGlzfSULkJ30l9KDi1P94RZRpCKMVUdGUQG7WE6hpoWb+L1vd2J55CxIB3RBmFEweSf1If/ZCIiMhRuXv5KDl3JMVnDaf13V0E3v6E0LYGWjfsonXDLlxFefhO7kfBuH54Bhbr9mMKqOhKczYUoW1LfWzAuw/3Eg2EEvs8/QopOKU/vlP64S7OdzCliIhkKuNxUTB+AAXjBxCub6V1/U5a1u8ksq+VwMo6AivrcBV7Y33DRvfGO7wU49a4X8dCRVeasdYSqW8luKWeto/20ba5HsLRxH53qY/8MX0oOKU/ef2LHEwqIiLZxlPmo/is4RRNH0ZoeyOt63fR+sEeok1BAqu3E1i9HeN1k39COd5R5eQPL8PdSx/6k6WiKw1EmoIEt+6PFVpb6ok2tHXY7xlQjO+k3uSP7oOnb6GaeEVEpFsZY/AO7sO1k8wAAAiSSURBVIV3cC9KZp9AeEcTrRv30rZxD+HdgUQHfAB3eQHeEaXkjyjDO6yUN1bXUFVVxaxZs5g+fbrD30l6UdHVw2wkSnhXgGBdA6G6RkJ1DUT2dyyyjM+Dd3jsP3D+Cb31KUJERBxjjCFvUAl5g0ooOWcE4X0ttG3aS3DLfoJb9xPZ9/+3d/cxctR1HMff35nZh7u96xPtIS21YoREVBJTUoVifEBN1Qr/mIBGQmKTJgYjJBIi8o/xHzUmPkRIlCCJDyRIFGODkgoifykKYjHWCjaKtmBT2mto727vdmfm6x8zd13u9u527647t+XzSjY7s/Ow373fbvZz85udX536qTr1A8d45j9/4+P3fBaAgUqV/Q//imt2vV/nGucUus4hbybEJyZoHh8nPj5O8/gY8Svj2ZiHLawUUNq8hvK2dVS2rSMaqekNKiIiq1K0foBo+xZq27fgSUrz2NhMT83vnzwws16j0eDR7/6Uy/45QHRhjdLIENFIjdKmWtZr8zo8L0yhawWkUzHxaJ3k5ATxaD27nZwgOVUHn7t+uK5KafMw5c1rKG0eJtqkkCUiIv3HwoDyljWUt6xh6Oo3snvrTXz1198HoFwqs/OKHXgjoXnkNM0jp89uGBjRhgGiCwYJNwycnV4/cF4PUaTQ1YG0kZCOTZGcbpC8OtlymyJ5dfI1vyh8DYNo4yDRphqlkRrRyBClkRrBYKm3L0BERKQHrr5m58z0b598gquuuopkrEF8fKyl12ec5FSd+MQE8YmJOfsIamXCtRXCtVXCddXsfm2FcLhCOFTu60G7OwpdZrYL+A4QAve5+9dmLa8APwK2AyeBG9z9xXzZncAeIAE+7+77V6z6JXJ3iFPSekxab5JOTN8apONNkokm6ViDdKxBcmYKbyQL7zAKiNZXiTbkif2CwSy1bxjo6zeHiIjIUk2fRB8OlQmHNlB584aZZWkjyXqERuvEoxPEJ+vZ9Kk66XiDdLxB8+Uzbfdr1YhwqEwwXCGslQgGywS1EsFgPj1YIhiICAZKWLS6ujAXDV1mFgL3AB8CjgJPm9k+d/97y2p7gFPu/hYzuxH4OnCDmV0O3Ai8DdgMPG5ml7n7IilmLneHxPE4xZtJfp9PN5KW+5S0keBTMT6VkE7FeCMmnUrwyTgLWpNNSNr0+80nCvIGLhOuqRJNp+48gQdDZf2iUEREpENBOaR80TDMGqbOUyc9M0Xc0ps03bOUjk2RjDXwyZh4MoY2R8nmiAKCakQwEGHViKAcYZWQoJLdWzkiKIdYOcRKQX6fz0dB9liULSOwZX/Xd3Kkawdw2N3/BWBmDwLXA62h63rgy/n0z4C7LavseuBBd58C/m1mh/P9/WG+J4tPTnDi/mezUJWkeJxCks45+XzZQiOo5ml4sERQy9NxfguHygRDFcLhMlaNFKpERETOMQss706stl3u7lnPVN4TlY639FRNNM/OTzZJ63HWq5X3XC2/uOzyGJv2bF/yLjoJXVuAIy3zR4F3zbeOu8dm9ipwQf74U7O23bLQk3kzJT45T3oN7GzqnE6gpfYJNaiEWCVPsHmizdJuiaAaQRQoSImIiPQRMyOslQlrZUoXDi24rrtnPWJ5AEsnY3wqbtMblveavabnLMXjll61OM3GOE676CVro5PQ1S6ZzH7W+dbpZFvMbC+wN58du+iO9zzfQV39aiNwougilqObsLrQuisVetvtZ/Zji8138RyLtl+3z71YLUv5O61UG3W6XjevaSltsdy/Ub68o89e6766fd75tu3kPboU3ex3KTV0s6/FrMD7fsH2O1d/406txOfofDPrtfb9d9+MvXMe2dbppp2ErqPA1pb5i4GX51nnqJlFwFpgtMNtcfd7gXs7Lbqfmdkz7n5l0XXI0qj9+pfarr+p/fqb2i/TyWn9TwOXmtklZlYmOzF+36x19gE359OfAJ5wd88fv9HMKmZ2CXAp8KeVKV1ERESkfyx6pCs/R+tzwH6yS0bc7+4HzewrwDPuvg/4AfDj/ET5UbJgRr7eQ2Qn3cfALUv55aKIiIhIv7PsgJT0ipntzbtTpQ+p/fqX2q6/qf36m9ovo9AlIiIi0gOr61KtIiIiIucpha4CmdntZuZmtrHoWqRzZvYNM/uHmf3VzH5hZuuKrkkWZma7zOx5MztsZl8suh7pnJltNbPfmdkhMztoZrcWXZN0x8xCM/uLmT1SdC1FU+gqiJltJRta6b9F1yJdewx4u7tfAbwA3FlwPbKAlqHMPgJcDnwyH6JM+kMMfMHd3wq8G7hF7dd3bgUOFV3EaqDQVZxvAXfQ5mKxsrq5+2/cPc5nnyK7/pysXjNDmbl7A5geykz6gLv/z92fzafPkH15LziyiaweZnYx8DHgvqJrWQ0UugpgZtcBL7n7c0XXIsv2GeDRoouQBbUbykxf2n3IzN4EvBP4Y7GVSBe+TXaAYYUHUO5PnVyRXpbAzB4H3tBm0V3Al4AP97Yi6cZC7efuv8zXuYus6+OBXtYmXetoODJZ3cxsCPg5cJu7ny66Hlmcme0Gjrv7n83sfUXXsxoodJ0j7v7Bdo+b2TuAS4Dn8nGpLgaeNbMd7n6shyXKAuZrv2lmdjOwG7jWdd2V1a6j4chk9TKzElngesDdHy66HunYTuA6M/soUAXWmNlP3P3TBddVGF2nq2Bm9iJwpbufHwOBvg6Y2S7gm8B73f2VouuRheXjwb4AXAu8RDa02afc/WChhUlHLPvv9IfAqLvfVnQ9sjT5ka7b3X130bUUSed0iXTvbmAYeMzMDpjZ94ouSOaX/+hheiizQ8BDClx9ZSdwE/CB/PN2ID9yItJ3dKRLREREpAd0pEtERESkBxS6RERERHpAoUtERESkBxS6RERERHpAoUtERESkBxS6RERERHpAoUtERESkBxS6RERERHrg/8Uxc5bLj3vQAAAAAElFTkSuQmCC\n",
"text/plain": [
"