{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyPPbUR8ohtkP9FJR6sP2FYt"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["# Sampling standard Gaussian random variables"],"metadata":{"id":"856tFswmENXK"}},{"cell_type":"code","execution_count":1,"metadata":{"id":"V0_GOOUdY9X_","executionInfo":{"status":"ok","timestamp":1697483838271,"user_tz":-120,"elapsed":399,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"outputs":[],"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","import math\n","from scipy.stats import norm # To access the CDF of the standard Gaussian distribution"]},{"cell_type":"markdown","source":["Throughout this tutorial, we are only allowed to use Numpy's random number generator to sample uniform random variables."],"metadata":{"id":"OfdXU_IGvjye"}},{"cell_type":"code","source":["def Unif():\n","  return np.random.uniform()"],"metadata":{"id":"65qU4uPjvxW0","executionInfo":{"status":"ok","timestamp":1697483844843,"user_tz":-120,"elapsed":9,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":2,"outputs":[]},{"cell_type":"code","source":["Unif()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"-LmH8Ve0v0y9","executionInfo":{"status":"ok","timestamp":1697483847416,"user_tz":-120,"elapsed":580,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"90b0b048-a8d1-4877-d291-61733357bbc7"},"execution_count":3,"outputs":[{"output_type":"execute_result","data":{"text/plain":["0.27943576228382394"]},"metadata":{},"execution_count":3}]},{"cell_type":"markdown","source":["The purpose of the tutorial is merely to implement the Box--Muller and Ziggurat algorithms, we will not try to optimise their efficiency."],"metadata":{"id":"x9GuSRDixM9M"}},{"cell_type":"markdown","source":["## 1. The Box--Muller method\n","\n","We recall that the Box--Muller method is based on the remark that if $(U_1, U_2)$ are independent $\\mathcal{U}[0,1]$ variables, then\n","$$X_1 = \\sqrt{-2\\log(U_1)} \\cos(2\\pi U_2), \\qquad X_2 = \\sqrt{-2\\log(U_1)} \\sin(2\\pi U_2),$$\n","are independent $\\mathcal{N}(0,1)$ variables.\n","\n","Complete the function `BoxMuller()` so that it returns the variable $X_1$ as above."],"metadata":{"id":"e9XL05DXL_CB"}},{"cell_type":"code","source":["def BoxMuller():\n","  U1 = Unif()\n","  U2 = Unif()\n","  return (-2*np.log(U1))**.5*np.sin(2*math.pi*U2)"],"metadata":{"id":"Mklav__ovfcr","executionInfo":{"status":"ok","timestamp":1697483869727,"user_tz":-120,"elapsed":827,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":4,"outputs":[]},{"cell_type":"code","source":["BoxMuller()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"hfv2ACrNwcgN","executionInfo":{"status":"ok","timestamp":1697483872554,"user_tz":-120,"elapsed":259,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"ec79fc0e-20f8-4a66-d7f9-6da70cb86d34"},"execution_count":5,"outputs":[{"output_type":"execute_result","data":{"text/plain":["0.25125458138828327"]},"metadata":{},"execution_count":5}]},{"cell_type":"markdown","source":["We check the validity of the function `BoxMuller()` by plotting the histogram of 1000 realisations together with the density of the $\\mathcal{N}(0,1)$ law."],"metadata":{"id":"vpE-9zbFxjH8"}},{"cell_type":"code","source":["# Sample\n","N = 5000\n","Sample = np.zeros(N)\n","for k in range(N):\n","  Sample[k] = BoxMuller()\n","\n","plt.hist(Sample, bins=50, density=True)\n","\n","# Theoretical density\n","x_max = max(abs(Sample))\n","xx = np.linspace(-x_max, x_max, 500)\n","den = np.exp(-xx**2/2)/(2*math.pi)**.5\n","\n","plt.plot(xx,den,color='red')\n","\n","plt.show()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":430},"id":"C9FKpNdbx4To","executionInfo":{"status":"ok","timestamp":1697487029771,"user_tz":-120,"elapsed":3163,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"5e8cb37e-e9d7-4af2-bed3-28e51b0d6af5"},"execution_count":40,"outputs":[{"output_type":"display_data","data":{"text/plain":["<Figure size 640x480 with 1 Axes>"],"image/png":"iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOnklEQVR4nO3de1hVZcL+8e8GBEQFDyiIknjKM1KgiIfUotBqpnyrsX41Fm+vzVjNWy9TjTap01SDljVO5Wg542Rnp6ax41BGopl4CCXNU2qSooJgAYoJCvv3x2JvJVHZCDz7cH+ua1/7tPbi3oLsm/U8ay2b3W63IyIiIuLG/EwHEBERETkfFRYRERFxeyosIiIi4vZUWERERMTtqbCIiIiI21NhEREREbenwiIiIiJuT4VFRERE3F6A6QCNobq6mgMHDtCmTRtsNpvpOCIiIlIPdrudI0eOEBUVhZ/fubeheEVhOXDgANHR0aZjiIiISAPs27ePrl27nnMZrygsbdq0Aaw3HBoaajiNiIiI1EdZWRnR0dHOz/Fz8YrC4hgGCg0NVWERERHxMPWZzqFJtyIiIuL2VFhERETE7amwiIiIiNtTYRERERG3p8IiIiIibk+FRURERNyeCouIiIi4PRUWERERcXsqLCIiIuL2GlRY5s2bR0xMDMHBwSQmJrJu3bp6ve7NN9/EZrNx/fXX13rcbrczY8YMOnfuTMuWLUlOTmbnzp0NiSYiIiJeyOXCsmTJEtLS0pg5cyYbNmxg8ODBpKSkcOjQoXO+Li8vjwceeIBRo0ad8dyTTz7Js88+y4IFC1i7di2tWrUiJSWF48ePuxpPREREvJDLheWZZ55h8uTJpKam0r9/fxYsWEBISAiLFi0662uqqqq49dZbefTRR+nRo0et5+x2O3PnzuWRRx7huuuuIzY2lpdffpkDBw6wdOlSl9+QiIiIeB+XCktlZSU5OTkkJyefWoGfH8nJyWRnZ5/1dX/84x/p1KkTd9555xnP7dmzh4KCglrrDAsLIzEx8azrrKiooKysrNZFREREvJdLZ2suLi6mqqqKiIiIWo9HRESwffv2Ol+zatUq/v73v5Obm1vn8wUFBc51/HSdjud+Kj09nUcffdSV6CLiaTZuhH/9C3JzoaQE2rWD+Hi48UYYONB0OhFpZi4VFlcdOXKEX/7ylyxcuJDw8PBGW++0adNIS0tz3i8rKyM6OrrR1i8iBu3cCb/5DXz88ZnPffABPPooXHcd/OUv0K1b8+cTESNcKizh4eH4+/tTWFhY6/HCwkIiIyPPWH737t3k5eXxs5/9zPlYdXW19YUDAtixY4fzdYWFhXTu3LnWOuPi4urMERQURFBQkCvRRcQTLFgA990HlZXQooVVTK64Ajp2hIICWLbMKi3vvguffAILF8Ktt5pOLSLNwKU5LIGBgcTHx5OZmel8rLq6mszMTJKSks5Yvm/fvmzevJnc3Fzn5ec//zljx44lNzeX6OhounfvTmRkZK11lpWVsXbt2jrXKSLeJWbqh3T/3fssGHYjTJkClZWs6H4pY+54npiedxCTF03M+mBi9sUQ03cyV97+LGujB8KPP8JttzHnsl+C3W76bYhIE3N5SCgtLY3bb7+dhIQEhg4dyty5cykvLyc1NRWASZMm0aVLF9LT0wkODmbgT8aa27ZtC1Dr8fvvv5/HH3+c3r170717d6ZPn05UVNQZx2sREe/0uxWL+fXafwEwZ9RtPJ80EWy2Opfd2bEbN9/yJ36X9RK/XvcOD3z+KqT3g4cfbs7IItLMXC4sEydOpKioiBkzZlBQUEBcXBwZGRnOSbN79+7Fz8+1vaUfeughysvLueuuuygpKWHkyJFkZGQQHBzsajwR8TCTct53lpWHxv0v/xx81XlfY7f5MWvsf1PUqh3Tl/8dfv976NoVJk1q6rgiYojNbvf8ballZWWEhYVRWlpKaGio6TgiUl8rV1I1Ziz+9mqevGwSf036hcurmLp8Eb9e94415yU729qTSEQ8giuf3zqXkIiY8cMPcOut+Nur+dfAy/nrsJsatJrZY+6ACRPgxAm45RY4erRxc4qIW1BhEREzfvUryM/n23ZRTL9yylnnrJyP3eYHf/ubNSS0cyfcf3/j5hQRt6DCIiLN7/334a23ICCA+372IMcCW17Y+tq3h1dftUrP3/8OK1Y0Tk4RcRsqLCLSvI4dg//9X+t2WhqbO/dunPWOHg133WXdvvtua4hIRLyGCouINK/0dMjLg+homD69cdf9pz9BeDhs3WodCVdEvIYKi4g0nwMH4Omnrdt//jO0bt2462/fHp580rr9xBPWxF4R8QpNei4hEfFtMVM/rHX/8Y/ncduPP/Jll37cuC4I1n94lldegEmTrDK0eTPMng2zZjX+1xCRZqctLCLSLC764SATN30CwJOjb2/wXkHn5e9vbV0BePZZOHiwab6OiDQrFRYRaRb/u/pNWlRXsbxHPOuiB57/BRfi2mth+HDrfEPawiLiFTQkJCJNLqrsENdtzQJg7oj/1+jr/+nQE8CI6Kt5jdX8+NcXGGFPZMOzjf91RaT5aAuLiDS5yev+TYvqKr7oFstXUX2a5Wt+0W0wX0X2puXJCm7Peb9ZvqaINB0VFhFpUu2OlXLzV9bclb8Oc/1cQQ1mszF/2I0A3JHzPhw50nxfW0QanQqLiDSp/5ebQcuTFWyK7MUX3QY369f+pPcwdrfvQlhFOfzjH836tUWkcamwiEiT8a+u4tbc/wCwKOG6ptsz6Cyq/fz5R8J11p3nn4fq6mb9+iLSeFRYRKTJXLFrHVFHijncMpT/9BlhJMM7A8ZSFhhinRhx2TIjGUTkwqmwiEiTuW3jRwD8M/YqKgICjWQ4FtiStwclW3eee85IBhG5cCosItI0vvmGy/I2Uo2N1y4ZbzTKK5deY9346CPYs8doFhFpGBUWEWka8+cD8FnPBPLDIoxG2dO+CyQng90OixcbzSIiDaPCIiKN78cf4aWXAHj1kmvMZnFITbWuX3pJk29FPJAKi4g0vnffhZIS8kM7saLHpabTWCZMgLAw+O47WL7cdBoRcZEKi4g0vpdfBuBfAy/HbnOTXzMtW8Itt1i3dUwWEY/jJr9JRMRrFBTAJ9aRbf89YKzhMD/hGBb617+gpMRoFBFxjQqLiDSuN96AqioYNoy89l1Mp6ltyBAYMACOH4clS0ynEREXqLCISOOqGQ7il780m6MuNtuprSwaFhLxKCosItJ4Nm+G3Fxo0QImTjSdpm633gp+frB2rY7JIuJBVFhEpPG88op1fe210KGD2SxnExkJY8ZYt//5T6NRRKT+VFhEpHFUV1vzV8A9h4NO59j6o3ksIh4jwHQAEfES69ZBfj60bg3jzR6Kvy4xUz903m53rDXrbX4EbNzI2LtetI6EC+TNcpOD3InIGbSFRUQax1tvWdc/+xkEB5vNch4/hISxKuYSAK7dttJwGhGpDxUWEblwdju8/bZ1+6abzGappw/6jQLg2u2fG04iIvWhwiIiF279eti7F1q1gnHjTKepl096D6PCP4A+xXvpXfSd6Tgich6awyIiDXL6nJBpyxfxK+D96Ev5zaOfmQvlgrLg1qzsfilX7lrHtds/588du5mOJCLnoC0sInJh7Hau3vEFAB/2HWk4jGs+6mPlvWrnGsNJROR8GlRY5s2bR0xMDMHBwSQmJrJu3bqzLvvOO++QkJBA27ZtadWqFXFxcbziOFZDjTvuuAObzVbrMs5DNiuL+LpBBbuILi3kWIsgsnrEm47jks96DuGkzY9+RXlElxSYjiMi5+ByYVmyZAlpaWnMnDmTDRs2MHjwYFJSUjh06FCdy7dv357f//73ZGdns2nTJlJTU0lNTeXjjz+utdy4ceM4ePCg8/KG43gOIuLWxn9jbV1Z3iOB4y3ce++gnypt2YZ10QMBbWURcXcuF5ZnnnmGyZMnk5qaSv/+/VmwYAEhISEsWrSozuXHjBnDhAkT6NevHz179uS+++4jNjaWVatW1VouKCiIyMhI56Vdu3YNe0ci0qyu2GVtYf344iTDSRrmk97DABUWEXfnUmGprKwkJyeH5OTkUyvw8yM5OZns7Ozzvt5ut5OZmcmOHTu47LLLaj2XlZVFp06d6NOnD1OmTOHw4cNnXU9FRQVlZWW1LiLS/KJLCuhTvJeTNj+yeiSYjtMgy2oKS0L+ViguNpxGRM7GpcJSXFxMVVUVERERtR6PiIigoODs47+lpaW0bt2awMBArrnmGp577jmuvPJK5/Pjxo3j5ZdfJjMzk9mzZ7NixQrGjx9PVVVVnetLT08nLCzMeYmOjnblbYhII0netRaA9dEDKAtubThNw+wP68SWTj3wt1fDBx+YjiMiZ9EsuzW3adOG3Nxcjh49SmZmJmlpafTo0YMxNScgu/nmm53LDho0iNjYWHr27ElWVhZXXHHFGeubNm0aaWlpzvtlZWUqLSIGOArLp70SDSe5MJ/0HsaAQ9/C0qVwxx2m44hIHVzawhIeHo6/vz+FhYW1Hi8sLCQyMvLsX8TPj169ehEXF8dvf/tbbrzxRtLT08+6fI8ePQgPD2fXrl11Ph8UFERoaGiti4g0r9DjRxm6bwsAn/YaajjNhXEMC/HJJ3DsmNkwIlInlwpLYGAg8fHxZGZmOh+rrq4mMzOTpKT6T7irrq6moqLirM/n5+dz+PBhOnfu7Eo8EWlGo7/NoUV1FTs7RPNduyjTcS7I1k7dyQ/tBD/+CJ9+ajqOiNTB5b2E0tLSWLhwIYsXL2bbtm1MmTKF8vJyUlNTAZg0aRLTpk1zLp+ens6yZcv49ttv2bZtG08//TSvvPIKt912GwBHjx7lwQcfZM2aNeTl5ZGZmcl1111Hr169SElJaaS3KSKN7Yrd1t5BmR6+dQUAm43MXkOs2//5j9ksIlInl+ewTJw4kaKiImbMmEFBQQFxcXFkZGQ4J+Lu3bsXP79TPai8vJy7776b/Px8WrZsSd++fXn11VeZOHEiAP7+/mzatInFixdTUlJCVFQUV111FY899hhBQUGN9DZFpFGdOMHY3V8CsMzD5684LO+RwO0bPoSPPrJO5mizmY4kIqex2e12u+kQF6qsrIywsDBKS0s1n0WkOaxcCaNHc7hlKEPufYVqP3/TiS5Y8InjbP/rbXD8OHz9NQwYYDqSiNdz5fNb5xISEddlZACwsvulXlFWAOsovWPHWnc++shsGBE5gwqLiLiuprCs8LBzB53X1Vdb1yosIm5HhUVEXFNYCBs3AvB5zCWGwzSy8eOt61WroLTUbBYRqUWFRURc88knAGyO6MnhVm3NZmlsPXvCxRfDyZPavVnEzaiwiIhras60vrL7pYaDNBENC4m4JRUWEam/6mpnYfG6+SsOjsLyn/9YuzeLiFtQYRGR+tu40TqjcZs2bIjqazpN07jsMmjZEg4ehC1bTKcRkRoqLCJSfzV7B3HFFZz0b5Zzpza/oCCrtAAsW2Y2i4g4qbCISP3VDAfh7afNuPJK61qFRcRtqLCISP2UlsLq1dZtby8sycnW9YoVUFlpNouIACosIlJfn30GVVXWbr/du5tO07QGDYJOneDYMcjONp1GRFBhEZH6csxfGTfObI7m4Od3aiuLhoVE3IKXzpoTkUZXc8A4bx4Oipn6ofP2TT905Clg46K3mXAyCYC8WdcYSiYi2sIiIuf37beQlwcBAaf2oPFyjtMOxBbsJPT4UcNpRESFRUTO77PPrOthw6B1a7NZmklBaDi72nfF315N0nebTMcR8XkqLCJyfpmZ1vXll5vN0cw+725tZRmVt9FwEhHRHBYRqZNzPofdzvr3MugI/GJXCOtOm+fh7b7oFkdqzvuM+C7XdBQRn6ctLCJyThcXf0fHYyX8GBBEbuc+puM0qzUXDeKkzY/uPxyka2mh6TgiPk2FRUTOaXjN/I31XftTGdDCcJrmdTQohI0150wauUfDQiImqbCIyDmN+O4rAFZ3G2w4iRlfxFjv2/HvICJmqLCIyFn5V1eRuHczAKu7xRpOY4ajqA3btxnsdsNpRHyXCouInNXAgl2EVh6jLKgVX0f0NB3HiNzOfTgeEEjH8hLYvt10HBGfpcIiImc1fK81f2XNRYOo9vM3nMaMyoAW5HSx5rGQlWU0i4gvU2ERkbManmfN2/jCR+evOKyJHmTdWL7cbBARH6bCIiJ1Cjx5giH7twKw+iLfnL/ikO2Yv5OVpXksIoaosIhInS49sI3gk5UUtWrLzvCLTMcxalPkxfwYEARFRbB1q+k4Ij5JhUVE6pT0nbV3UPZFsWCzGU5jVmVAC77s0s+6o3ksIkaosIhInYbmfw1YE27ltH8HzWMRMUKFRUTOVFHBJQd2ALA2eqDhMO4h2zGPZ8UKqK42G0bEB6mwiMiZ1q+35q+EtGV3+66m07iFzZ17QUgIFBfDli2m44j4HBUWETnTihUArIse4PPzVxxO+LeAkSOtO5rHItLsVFhE5EwrVwKwTsNBtY0ZY11rHotIs1NhEZHaTpyAL74ANH/lDGPHWteaxyLS7BpUWObNm0dMTAzBwcEkJiaybt26sy77zjvvkJCQQNu2bWnVqhVxcXG88sortZax2+3MmDGDzp0707JlS5KTk9m5c2dDoonIhdqwAcrLKQluzY6O3UyncS/x8dCqFXz/PWzebDqNiE9xubAsWbKEtLQ0Zs6cyYYNGxg8eDApKSkcOnSozuXbt2/P73//e7Kzs9m0aROpqamkpqby8ccfO5d58sknefbZZ1mwYAFr166lVatWpKSkcPz48Ya/MxFpmJrhoPVdB2C3aSNsLS1awIgR1u3PPzebRcTHuPzb6JlnnmHy5MmkpqbSv39/FixYQEhICIsWLapz+TFjxjBhwgT69etHz549ue+++4iNjWXVqlWAtXVl7ty5PPLII1x33XXExsby8ssvc+DAAZYuXXpBb05EGqBmwu0aDQfVbdQo61qFRaRZuVRYKisrycnJITk5+dQK/PxITk4mOzv7vK+32+1kZmayY8cOLrvsMgD27NlDQUFBrXWGhYWRmJhYr3WKSCOqqoKaPyY04fYsHIVl5UqdV0ikGQW4snBxcTFVVVVERETUejwiIoLt27ef9XWlpaV06dKFiooK/P39+etf/8qVV14JQEFBgXMdP12n47mfqqiooKKiwnm/rKzMlbchImezaROUlkKbNmyN6GE6jXsaOhQCA6GgAHbvhl69TCcS8QnNMkDdpk0bcnNzWb9+PU888QRpaWlkXcBxDNLT0wkLC3NeoqOjGy+siC+rmb/CiBFU+fmbzeKuWraEIUOs2xoWEmk2LhWW8PBw/P39KSwsrPV4YWEhkZGRZ/8ifn706tWLuLg4fvvb33LjjTeSnp4O4HydK+ucNm0apaWlzsu+fftceRsicjY181cYPdpsDnd3+rCQiDQLlwpLYGAg8fHxZGZmOh+rrq4mMzOTpKSkeq+nurraOaTTvXt3IiMja62zrKyMtWvXnnWdQUFBhIaG1rqIyAWy2099ANfMMZOzcPz7aAuLSLNxaQ4LQFpaGrfffjsJCQkMHTqUuXPnUl5eTmpqKgCTJk2iS5cuzi0o6enpJCQk0LNnTyoqKvjoo4945ZVXmD9/PgA2m43777+fxx9/nN69e9O9e3emT59OVFQU119/feO9UxE5t61b4fBha8gjIQHeW2Y6kfsaPtw6ZcHu3XDgAERFmU4k4vVcLiwTJ06kqKiIGTNmUFBQQFxcHBkZGc5Js3v37sXP79SGm/Lycu6++27y8/Np2bIlffv25dVXX2XixInOZR566CHKy8u56667KCkpYeTIkWRkZBAcHNwIb1FE6sWxdWX4cGtSqZxdWBgMHgy5udZWltN+n4lI07DZ7Z6/X15ZWRlhYWGUlpZqeEikoW6+GZYsgUcfhRkziJn6oelEbidv1jWn7tx3Hzz7LNxzDzz/vLlQIh7Mlc9vHcZSRGrPX9GE2/rRAeREmpUKi4jAnj1w8KB16PmhQ02n8QyOwrJ5M/zwg9ksIj7A5TksIuKFas7OTHy8NelW6vTTYbLP2kXR44cD/PfkuXzWyyp6tYaNRKTRaAuLiJwqLI4T+0m9OE5fMDR/i+EkIt5PhUVEThWW4cPN5vAw67sOAGDovq8NJxHxfiosIr6upAS21Gwh0BYWl6yLtgrLoIJdBJ84bjiNiHfTHBYRH3T6XIzR3+aw2G4nr21nxvz5S4OpPM++sAgOtu5A56OHueTAN2R3izUdScRraQuLiI+Lz98KQE7XfoaTeCCbjfXRGhYSaQ4qLCI+LmH/NgC+7NLfcBLPtK5mHssQTbwVaVIqLCI+LKDqJHEHdwDwZRdtYWmIL7taRe+SAzvwr64ynEbEe6mwiPiwfof2EHKigtKgVuwKjzYdxyN9E34RZYEhtDpxnL6H9piOI+K1VFhEfFjCfmv+yoYufbHb9OugIar9/NnYpS9wanhNRBqffkOJ+LD4/dsBzV+5UI7htISaCcwi0vhUWER8ld2uPYQaiWMeS7y2sIg0GRUWER/VpayIzkcPc8LPn9zOF5uO49G+6nwxJ21+RB0phn37TMcR8UoqLCI+Kr5m/sqWiB4cbxFsOI1nOxbYkm2dult3HKc5EJFGpcIi4qMS8q3hixzNX2kUjmEhFRaRpqHCIuKjHHsI5ej4K43C+e+owiLSJFRYRHxQ64pj9Cn6DtAB4xqLc0+rr76CI0fMhhHxQiosIj4o7sAO/O3V7AuL4FCbDqbjeIWC0HDyQztCdTWsXWs6jojXUWER8UGO4SBtXWlczvlAGhYSaXQqLCI+KN4x4barJtw2pi+7ah6LSFNRYRHxNSdPcolOeNgknFtY1qyBKp0IUaQxqbCI+JrNm2ld+SNlgSF8E36R6TReZXvHbtCmjTXpdvNm03FEvIoKi4ivqRmuyI3qQ7Wfv+Ew3qXazx+GDbPuaFhIpFGpsIj4mtWrgdMOdCaNa8QI61qFRaRRqbCI+JqaD1LNX2kiKiwiTUKFRcSX5OfD3r2ctPmRG9XHdBrvlJgIfn6wd6/17y0ijUKFRcSX1PzVv61Td44FtjQcxku1aQODB1u3a4bfROTCqbCI+BLHcJDmrzQtDQuJNDoVFhFfUvMBuiGqr+EgXk6FRaTRqbCI+IqjR60T86EtLE3OUVhyc61/dxG5YCosIr5i3Trr6KvR0RwM7Wg6jXeLjoauXa1/73XrTKcR8QoqLCK+wjE84fjrX5qWhoVEGlWDCsu8efOIiYkhODiYxMRE1p3jL4iFCxcyatQo2rVrR7t27UhOTj5j+TvuuAObzVbrMm7cuIZEE5GzUWFpXiosIo3K5cKyZMkS0tLSmDlzJhs2bGDw4MGkpKRw6NChOpfPysrilltuYfny5WRnZxMdHc1VV13F/v37ay03btw4Dh486Ly88cYbDXtHInKmqirIzrZuq7A0D8e/c3a2ToQo0ghcLizPPPMMkydPJjU1lf79+7NgwQJCQkJYtGhRncu/9tpr3H333cTFxdG3b1/+9re/UV1dTWZmZq3lgoKCiIyMdF7atWvXsHckImfasgXKyqBVKxg0yHQa3xAbC61bW//uW7aYTiPi8QJcWbiyspKcnBymTZvmfMzPz4/k5GSyHX+9ncexY8c4ceIE7du3r/V4VlYWnTp1ol27dlx++eU8/vjjdOjQoc51VFRUUFFR4bxfVlbmytsQ8T2OYYlhwyDApf/24qKYqR86b7/aoScjj37F7x98gdcuudr5eN6sa0xEE/FoLm1hKS4upqqqioiIiFqPR0REUFBQUK91/O53vyMqKork5GTnY+PGjePll18mMzOT2bNns2LFCsaPH0/VWTajpqenExYW5rxER0e78jZEfI/jiKsaDmpWOV2s3cfj928znETE8zXrn1qzZs3izTffJCsri+DgYOfjN998s/P2oEGDiI2NpWfPnmRlZXHFFVecsZ5p06aRlpbmvF9WVqbSInIumnBrxJddrRNMJuRvNZxExPO5tIUlPDwcf39/CgsLaz1eWFhIZGTkOV87Z84cZs2axSeffEJsbOw5l+3Rowfh4eHs2rWrzueDgoIIDQ2tdRGRszh4EPbssU7IN2yY6TQ+ZWNUX6qxcVFpIR2Pfm86johHc6mwBAYGEh8fX2vCrGMCbVJS0llf9+STT/LYY4+RkZFBQkLCeb9Ofn4+hw8fpnPnzq7EE5G6OLauDBoEKvfN6mhQCDs6dgO0lUXkQrm8l1BaWhoLFy5k8eLFbNu2jSlTplBeXk5qaioAkyZNqjUpd/bs2UyfPp1FixYRExNDQUEBBQUFHK05XPXRo0d58MEHWbNmDXl5eWRmZnLdddfRq1cvUlJSGultivgwDQcZ5TgNguaxiFwYl+ewTJw4kaKiImbMmEFBQQFxcXFkZGQ4J+Lu3bsXP79TPWj+/PlUVlZy44031lrPzJkz+cMf/oC/vz+bNm1i8eLFlJSUEBUVxVVXXcVjjz1GUFDQBb49Ed9z+l4qAEuXfEgccN/elrz7k+ek6X3ZpR+/3PgRCSosIhekQZNu7733Xu699946n8vKyqp1Py8v75zratmyJR9//HFDYojIeQSfOM6Awm8ByNEJD41w/LsPKNxN8InjHG8RfJ5XiEhddC4hES82+OBOWlRXUdC6PfmhnUzH8Un5oZ0oaN2eFtVVDD6403QcEY+lwiLixRzzJr7s0h9sNsNpfJTNRk4Xa/dmzWMRaTgVFhEv5tgzJafmeCBihg4gJ3LhVFhEvJTNXl17C4sY4ziAXPz+bdjs1YbTiHgmFRYRL9WreB9hFeUcaxHEtk7dTcfxaVs79eBYiyDaHj9Kz8P5puOIeCQVFhEv5diN9qvOF3PSXyc8NOmkfwBfdb4Y0AHkRBpKhUXESyXstz4YNRzkHhwTb3U8FpGGUWER8VKX1nwwOj4oxawva74Pl6qwiDSICouIFwov/4HuPxykGhsbuvQ1HUeADTWFpccPB+DQIcNpRDyPCouIF3LsHfRN+EWUBbc2nEYAyoJbsyP8IuvO6tVmw4h4IBUWES8Un18zHKTjr7gVx/FYnCekFJF6U2ER8UKOCbeav+JeHMdjUWERcZ0Ki4iXCTpRwcCC3YD2EHI3zgKZkwPHj5sNI+JhVFhEvMygwl0EVp+kqFVb9raNNB1HTvNd284UhbSFykqrtIhIvamwiHiZhHyd8NBt2Wyn5hVpWEjEJSosIl4m3nnAOM1fcUfO74sKi4hLVFhEvIndTvz+7QDkdNX8FXfk3FNo9Wqw282GEfEgKiwi3mTHDtr/WMbxgEC2RPQwnUbqsCWiJwQFQXExfPON6TgiHkOFRcSb1AwzfBXZmxP+LQyHkbpUBrSAIUOsOzqAnEi9qbCIeJOawqIDxrm5ESOsa81jEak3FRYRb7JqFQDruw4wHETOSYVFxGUqLCLe4tAh2LkT0BFu3V5SknW9fTscPmw2i4iHUGER8RY1f63v0AkP3V94OPTpY93WPBaRelFhEfEWNYXlS+3O7Bkcw0IqLCL1osIi4i00f8WzaB6LiEtUWES8wbFjznPT6Ai3HsJRWNavt84tJCLnpMIi4g3Wr4eTJyEqivywCNNppD4uvtiay3L8OGzYYDqNiNtTYRHxBjXDQYwcqRMeegqbDYYPt25rWEjkvFRYRLyB4wPPMcwgnsFRWDTxVuS8VFhEPF1V1akPvJEjzWYR15w+8VYnQhQ5JxUWEU+3ZQuUlkKrVhAbazqNuCIhAQIDobAQvv3WdBoRt6bCIuLpHMNBSUkQEGA2i7gmOBji463bmscick4qLCKe7vQJt+J5dDwWkXppUGGZN28eMTExBAcHk5iYyLp168667MKFCxk1ahTt2rWjXbt2JCcnn7G83W5nxowZdO7cmZYtW5KcnMzOmnOiiMh5aMKtx4mZ+qHz8qvdgQBsf+fjWo+LSG0uF5YlS5aQlpbGzJkz2bBhA4MHDyYlJYVDhw7VuXxWVha33HILy5cvJzs7m+joaK666ir279/vXObJJ5/k2WefZcGCBaxdu5ZWrVqRkpLC8ePHG/7ORHxBfj589x34+0Niouk00gCOE1VeXLyX0ONHDacRcV8uF5ZnnnmGyZMnk5qaSv/+/VmwYAEhISEsWrSozuVfe+017r77buLi4ujbty9/+9vfqK6uJjMzE7C2rsydO5dHHnmE6667jtjYWF5++WUOHDjA0qVLL+jNiXg9x9aVwYOhTRuzWaRBilu1Y0+7zvhh59L9203HEXFbLhWWyspKcnJySE5OPrUCPz+Sk5PJzs6u1zqOHTvGiRMnaN++PQB79uyhoKCg1jrDwsJITEw86zorKiooKyurdRHxSZq/4hVyulgnrIzfv81wEhH35VJhKS4upqqqioiI2of+joiIoKCgoF7r+N3vfkdUVJSzoDhe58o609PTCQsLc16io6NdeRsi3kOFxSs4zv+UsH+r4SQi7qtZ9xKaNWsWb775Jv/+978JDg5u8HqmTZtGaWmp87Jv375GTCniIcrKYNMm67Ym3Hq0L7taW1gGH/yGgKqThtOIuCeXCkt4eDj+/v4UFhbWerywsJDIyMhzvnbOnDnMmjWLTz75hNjTDm7leJ0r6wwKCiI0NLTWRcTnrFkD1dXQvTtERZlOIxdgd4eulAa1IuREBf0O7TEdR8QtuVRYAgMDiY+Pd06YBZwTaJOSks76uieffJLHHnuMjIwMEhISaj3XvXt3IiMja62zrKyMtWvXnnOdIr7KsdvrXx59CYB3QmK0O6yHs9v8nHsLaVhIpG4uDwmlpaWxcOFCFi9ezLZt25gyZQrl5eWkpqYCMGnSJKZNm+Zcfvbs2UyfPp1FixYRExNDQUEBBQUFHD1q7b5ns9m4//77efzxx3nvvffYvHkzkyZNIioqiuuvv75x3qWIF0rYvwU4NZwgns3xfYzP18Rbkbq4fBzviRMnUlRUxIwZMygoKCAuLo6MjAznpNm9e/fi53eqB82fP5/KykpuvPHGWuuZOXMmf/jDHwB46KGHKC8v56677qKkpISRI0eSkZFxQfNcRLxZQNVJLjmwA4D1XVRYvIFjC8uQ/Vt1IkSROtjsds//n1FWVkZYWBilpaWazyJeL2bqhww6uJP3X/4/SoNaEXffG9htOsuGpws6UcGmv0wkqOoko+96kRUvTDYdSaTJufL5rd9yIh5oSL41z+HLrv1VVrxERYsgvup8MQBD931tOI2I+9FvOhEPFF8zMdMxjCDeYV30QAAS920xnETE/aiwiHgau52hNR9o6zXh1qus6zoAgCH5KiwiP6XCIuJheny/n47HSqjwb8FXnfuYjiONKKdLP6psfnQrKbBObCkiTiosIh7GMb9hY1QfKgNaGE4jjak8KISvI3padz7/3GwYETejwiLiYYbWDBesrZnvIN5lXbQ1LMTKlWaDiLgZFRYRD+PYwrJOhcUrOb+vKiwitaiwiHiSvDy6lhVxws+fDVF9TaeRJuCcSL11KxQVmQ0j4kZUWEQ8Sc1f3Zsje/FjoI4E7Y1KWoayI/wi686qVWbDiLgRFRYRT1JTWDQc5N00LCRyJhUWEU9S8wGmCbfezXE8FhUWkVNUWEQ8xcGDsHMn1dh0hFsv59xTKDcXysqMZhFxFyosIp6i5rgc2zp1pyy4teEw0pQK24RDz55QXQ2rV5uOI+IWVFhEPIVz/soAw0GkWVx2mXWtYSERQIVFxHOsWAFo/orPUGERqUWFRcQTHD4MX1sHjFvfVVtYfIKjsKxbBz/+aDaLiBtQYRHxBI7jcfTrx+FWbY1GkWbSvTt06QInTsDatabTiBinwiLiCRzDAo6/usX72WwaFhI5jQqLiCdQYfFNo0ZZ1yosIiosIm7vyBHYsMG67fgAE9/gKKirV0NlpdksIoapsIi4uy++sI7HERMD0dGm00hz6tcPOnSwJt3m5JhOI2KUCouIu/vsM+t67FizOaT5+fnB6NHW7eXLzWYRMUyFRcTdOT6oVFh8k+P7rsIiPi7AdAAROYeSklPzV1RYfErM1A8B6FXsz6fA8ayVxD6wlMqAFgDkzbrGYDqR5qctLCLubOVKa/5K797QtavpNGLArg7RFLVqS/DJSuIO7jAdR8QYFRYRd+YYBrj8crM5xBybjeyLYgFI+m6T4TAi5qiwiLgzTbgVcBaW4XtVWMR3qbCIuKviYthU8wE1ZozRKGJW9kWDAIg7sJ3gE8cNpxExQ4VFxF3VnJ2ZAQMgIsJsFjEqr10UB9qEE1R1kvj9203HETFCewmJuBHHniEAf/zkH0wCXgrqzh9Oe1x8kM1G9kWDuGHLcpL2buKLmDjTiUSanbawiLippJr5CtndYg0nEXeQ3W0wAMO/+8pwEhEzVFhE3FDHoz/Q+/A+qrGxJnqQ6TjiBhwTb2MP7qRVxTHDaUSanwqLiBtybF3Z1qk7pS3bGE4j7mB/WCf2hkUQYK9mSP4W03FEml2DCsu8efOIiYkhODiYxMRE1q1bd9Zlt2zZwg033EBMTAw2m425c+eescwf/vAHbDZbrUvfvn0bEk3EKwyrKSyrNRwkp1ldMyyUtHez4SQizc/lwrJkyRLS0tKYOXMmGzZsYPDgwaSkpHDo0KE6lz927Bg9evRg1qxZREZGnnW9AwYM4ODBg87LqlWrXI0m4jWG1xwgzDEMIAKndm9O0vFYxAe5XFieeeYZJk+eTGpqKv3792fBggWEhISwaNGiOpcfMmQITz31FDfffDNBQUFnXW9AQACRkZHOS3h4uKvRRLxC57IiYkoOctLmx7rogabjiBtxFNiBBbvhhx8MpxFpXi4VlsrKSnJyckhOTj61Aj8/kpOTyc7OvqAgO3fuJCoqih49enDrrbeyd+/esy5bUVFBWVlZrYuIt3D89fx1ZC+OBoUYTiPu5FCbDuxu3xU/7NZ5pkR8iEuFpbi4mKqqKiJ+chCriIgICgoKGhwiMTGRl156iYyMDObPn8+ePXsYNWoUR44cqXP59PR0wsLCnJfo6OgGf20Rd6PhIDkXx7CQ87QNIj7CLfYSGj9+PDfddBOxsbGkpKTw0UcfUVJSwj//+c86l582bRqlpaXOy759+5o5sUgTsdudx9nQhFupi2PiLZmZZoOINDOXjnQbHh6Ov78/hYWFtR4vLCw854RaV7Vt25aLL76YXbt21fl8UFDQOefDiHiqnt/nE3WkmAr/FqzrOsB0HHFDq7vFUo0Nvy1b4MABiIoyHUmkWbi0hSUwMJD4+HgyT2v21dXVZGZmkpSU1Gihjh49yu7du+ncuXOjrVPEE4zasxGA9V37U9FCpVzOVNIylM2Rvaw7y5aZDSPSjFweEkpLS2PhwoUsXryYbdu2MWXKFMrLy0lNTQVg0qRJTJs2zbl8ZWUlubm55ObmUllZyf79+8nNza219eSBBx5gxYoV5OXlsXr1aiZMmIC/vz+33HJLI7xFEc8xMs8qLJ93v8RwEnFnzp8PFRbxIS6f/HDixIkUFRUxY8YMCgoKiIuLIyMjwzkRd+/evfj5nepBBw4c4JJLTv3ynTNnDnPmzGH06NFkZWUBkJ+fzy233MLhw4fp2LEjI0eOZM2aNXTs2PEC356IB6msZNi+rwFYFaPCIme3KiaOe7P/CZ9+CnY72GymI4k0OZvdbrebDnGhysrKCAsLo7S0lNDQUNNxRBpm5UoYPZrikDCG3PsKdptbzIkXNxR48gTfzL8Vjh2Dr76CWE3QFs/kyue3fiOKuIuazftfdItTWZFzqgxoAWPGWHc++cRoFpHmot+KIu6iprCsiokzm0M8w5VXWteaxyI+QoVFxB388AOsXw/A55q/IvXhKCwrV8Lx42aziDQDFRYRd/DZZ1Bdzc4O0RSE6jxaUg/9+1vHYDl+HL74wnQakSanwiLiDjQcJK6y2U5tZdE8FvEBKiwi7qDmA0fDQeISzWMRH6LCImLa7t2wZw+0aMEax4ntROojOdm63rgRiorMZhFpYiosIqY5NucnJXEssKXZLOJZIiJgcM3JELWVRbycCouIaf/5j3U9bpzZHOKZHD83GRlmc4g0MRUWEZOOHwfHyUTHjzebRTyT4+fmP/+B6mqzWUSakAqLiEkrV1qHV+/c+dSmfRFXDB8OoaFQXAxffmk6jUiTUWERMckxHDR+vE5gJw3TosWpvYUcP08iXsjlszWLSCP66CPr+uqrzeYQjxMz9UPn7ZuOdeEpIPeF17n+xwTn43mzrjGQTKRpaAuLiCm7d8M330BAwKndU0UaYEX3eABiD+6k/bFSw2lEmoYKi4gpjs33I0ZAWJjZLOLRDrXpwJZOPfDDzmV7NpiOI9IkVFhETNFwkDSi5T2toaCxuzXxVryTCouICT/+CMuXW7e1O7M0gqwe1rDQZXs24FddZTiNSONTYRExISvLOgZL164wcKDpNOIFNkb1pTSoFe2OHyHu4Dem44g0Ou0lJNJMTt+rY+anL5AKvN5hAA9P+8hcKPEaVX7+fN79Uq7d/jmjv81hQ5d+piOJNCptYRFpbna7c55BVo+E8ywsUn/La36exn6reSzifVRYRJpZz+/ziSk5SKVfAF9009FtpfGs6HEp1diILdhFpyOHTccRaVQqLCLN7MqdawHI7hZLeVCI4TTiTYpbtSM36mIAknevM5xGpHGpsIg0syt3rgFgWa9Ew0nEG31a83Pl+DkT8RYqLCLNKLz8By45sAM49cEi0pg+6T0MgOHffQVHjhhOI9J4VFhEmtHlu9bjh51Nkb0oCA03HUe80K4O0exp15mgqpPw8cem44g0GhUWkWZ05S4NB0kTs9lY1svaysJ775nNItKIVFhEmknLyuOMyssF4NPeKizSdJw/Xx9+CCdPmg0j0khUWESaycjvcgk+WUl+aCe2dexuOo54sZwu/fi+ZSh8/z2sWmU6jkijUGERaSbOvYN6J4LNZjiNeLMqP38+6znEuvPuu2bDiDQSFRaR5lBVxeW71wOn9uIQaUrLHMNC774LdrvZMCKNQIVFpDmsWUP4sVJKg1qxvusA02nEB3wecwkEBcGePbBli+k4IhdMhUWkOdRsll/eM4GT/jrnqDS9Y4EtITnZurN0qdEsIo1BhUWkqdnt8PbbAHzSO8lwGPEpEyZY1//6l9kcIo2gQYVl3rx5xMTEEBwcTGJiIuvWnf2cFVu2bOGGG24gJiYGm83G3LlzL3idIh5lwwbYs4cfA4KcZ9MVaRbXXQf+/pCbC7t2mU4jckFcLixLliwhLS2NmTNnsmHDBgYPHkxKSgqHDh2qc/ljx47Ro0cPZs2aRWRkZKOsU8SjvPUWAJ/1TODHwGDDYcSnhIfD5Zdbt2t+DkU8lcuF5ZlnnmHy5MmkpqbSv39/FixYQEhICIsWLapz+SFDhvDUU09x8803ExQU1CjrFPEYdrvzg+KjPiMNhxGfdNNN1rUKi3g4lwpLZWUlOTk5JDsmcgF+fn4kJyeTnZ3doAANWWdFRQVlZWW1LiJuaeNG+PZbaNny1HExRJrThAnWsNDGjbB7t+k0Ig3mUmEpLi6mqqqKiIiIWo9HRERQUFDQoAANWWd6ejphYWHOS3R0dIO+tkiTc/xVe/XVGg4SM8LDYcwY63bN5G8RT+SRewlNmzaN0tJS52Xfvn2mI4mc6bThIOdmeRETNCwkXsClwhIeHo6/vz+FhYW1Hi8sLDzrhNqmWGdQUBChoaG1LiJuJzfX2gQfHAzXXGM6jfiyCRPAzw9ycqwhShEP5FJhCQwMJD4+nszMTOdj1dXVZGZmkpTUsONLNMU6RdyC46/Z8eOhdWuzWcQnxUz90Lo8s54vogcCkH7nY6cen/qh4YQi9efykFBaWhoLFy5k8eLFbNu2jSlTplBeXk5qaioAkyZNYtq0ac7lKysryc3NJTc3l8rKSvbv309ubi67TjsmwPnWKeJxNBwkbsaxl9rV278wnESkYVw+RvjEiRMpKipixowZFBQUEBcXR0ZGhnPS7N69e/HzO9WDDhw4wCWXXOK8P2fOHObMmcPo0aPJysqq1zpFPM5XX1kH6goKgmuvNZ1GhIyLh/PHZQsYXLCT6JIC9rVt2DC+iCk2u93zT+NZVlZGWFgYpaWlms8ixpy+eX3a8kX8at07/Ofi4UyZ8LDBVCKnvPLmI4z6Lpc5o27j+eE3A5A3S/OrxBxXPr89ci8hEXfmV13Fz7euAGBp/zFmw4icZumAsQBM2JJlDVuKeBAVFpFGlrhvC52PHqY0qBVZPXXuIHEfH1+cxI8BQfT8Pp9BBTq3kHgWFRaRRnb9luUAfNRnBBUBgYbTiJxyNCiEZb0TAfivLZ8ZTiPiGhUWkUYUdLKS8TusvTDeHTDGbBiROrxTMyz0s20rCag6aTiNSP2psIg0orG71xNaeYwDbcJZW3PcCxF3sirmEopDwgg/VsrIvFzTcUTqTYVFpBFdvzULgPf6j8Zu038vcT8n/QN4v99lAEyoGb4U8QT6jSrSSNr+WMbY3esB7R0k7u3fNcNCV+1cA0eOGE4jUj8qLCKN5PotWQRVneTriJ5s79TddByRs9oU2Zvd7bvQ8mQF/PvfpuOI1IsKi0gjuWnzpwD8c1Cy4SQi52GzndoK+PLLRqOI1JcKi0hj2LiRAYe+pcI/gHc1HCQe4N8DL7duZGbCnj1mw4jUgwqLSGP4xz8AWNY7idKWbQyHETm//LAIPu8WZ92p+fkVcWcqLCIX6vhxePVVQMNB4lmWDL7KurFoEVRVmQ0jch4qLCIX6r334IcfONAmnFUxcabTiNTbJ72ToH172L8fPv7YdByRc1JhEblQf/sbAO8MvJxqP3/DYUTqrzKgBUyaZN2p+TkWcVcqLCIXYtcuWLYMbDbejL3KdBoR1915p3X9/vtQWGg2i8g5qLCIXIgXXrCux40jv22k2SwiDTFwICQmwsmT8NJLptOInJUKi0hDHT9+au+KKVPMZhG5EHfdZV0vWKDJt+K2VFhEGurtt+HwYYiOhquvNp1GpOFuvhnatYO8PPjoI9NpROqkwiLSUPPnW9d33QX+mmwrHiwk5NRclnnzzGYROYsA0wFEPNKmTbB6NQQEnPpFL+KBYqZ+CED0kb6swIbfxx8z5q4XyWvfBYC8WdeYjCfipC0sIg3x7LPW9fXXQ+fORqOINIZ9bSPJ6hEPwG0bNSwk7keFRcRVhw45j2zL/fcbjSLSmF6+9FoAfrH5U1pWHjecRqQ2FRYRVy1YABUVMGQIDB9uOo1Io1nR41Ly2nYmtKKcG77ONB1HpBbNYRGpB8c4f+DJE3yx4M90BP43cgzvTdOmc/EedpsfixJ+zh8/fYH/Wb+U1+PGmY4k4qQtLCIu+Pm2FXQsL+Fg6w581GeE6Tgije6tQVdSEtyamJKDXLlzrek4Ik4qLCL1Zbfz31++C8DL8ddy0l8bKMX7/BgYzCuXWHsGTV7/b8NpRE5RYRGpp1F5G+l/aA/HWgTx+mBtKhfv9fKl11LhH0DC/m3W7vsibkCFRaSe7sn+JwBvDB5Hacs2htOINJ2i1u1Y2n+sdWfOHLNhRGqosIjUQ3z+Vobt+5pKvwAWDplgOo5Ik1s4tObnfOlS2LrVaBYRUGERqZe717wFwL8GXk5BaLjhNCJNb1f4RWRcnAR2OzzxhOk4IiosIueVm8sVu9dTZfNjwbAbTacRaTbPJU20brz5Jnzzjdkw4vNUWETOJz0dgA/6juK7dlGGw4g0ny2RveDaa6G6Gv70J9NxxMepsIicy6ZN8E9rsu1fk24yHEbEgOnTretXX4VvvzWbRXxagwrLvHnziImJITg4mMTERNatW3fO5d966y369u1LcHAwgwYN4qOPah8d9I477sBms9W6jBun3UbFDdT8sn6/7yh2dIwxm0XEhKFDISUFqqo0l0WMcrmwLFmyhLS0NGbOnMmGDRsYPHgwKSkpHDp0qM7lV69ezS233MKdd97Jxo0buf7667n++uv5+uuvay03btw4Dh486Ly88cYbDXtHIo1l7Vp47z3w8+PPI281nUbEnJkzreuXXoLt241GEd/lcmF55plnmDx5MqmpqfTv358FCxYQEhLCokWL6lz+L3/5C+PGjePBBx+kX79+PPbYY1x66aU8//zztZYLCgoiMjLSeWnXrl3D3pFIY3FsCp80iW87dDWbRcSkpCS47jprLsvvf286jfgolwpLZWUlOTk5JCcnn1qBnx/JyclkZ2fX+Zrs7OxaywOkpKScsXxWVhadOnWiT58+TJkyhcOHD581R0VFBWVlZbUuIo1qxQpYtgxatDj116WIL3viCfDzg3fegfNMAxBpCi4VluLiYqqqqoiIiKj1eEREBAUFBXW+pqCg4LzLjxs3jpdffpnMzExmz57NihUrGD9+PFVVVXWuMz09nbCwMOclOjralbchcm7V1fDAA9bt//kfiIkxGkfELQwYAJMmWbenTrWOzyLSjNzi7G0333yz8/agQYOIjY2lZ8+eZGVlccUVV5yx/LRp00hLS3PeLysrU2mRxvPaa/Dll9CmjbauiM+Lmfqh83aX4Mv4zP9VgpYv545f/JGsngkA5M26xlQ88SEubWEJDw/H39+fwsLCWo8XFhYSGRlZ52siIyNdWh6gR48ehIeHs2vXrjqfDwoKIjQ0tNZFpFEcOwbTplm3H34YfrJ1UMSX7Q/rxOJLfwbA9M/+RkDVScOJxJe4VFgCAwOJj48nMzPT+Vh1dTWZmZkkJSXV+ZqkpKRaywMsW7bsrMsD5Ofnc/jwYTp37uxKPJEGiZn6ofPyTMpdsH8/+aGd6FPcx/m4iFieHz6R4pAwen6fz+0bPjAdR3yIy3sJpaWlsXDhQhYvXsy2bduYMmUK5eXlpKamAjBp0iSmOf5CBe677z4yMjJ4+umn2b59O3/4wx/48ssvuffeewE4evQoDz74IGvWrCEvL4/MzEyuu+46evXqRUpKSiO9TZHziywr5tdr3wZg1pg7qAgINJxIxP2UBbfmqcusuSz3rXqdDuUlZgOJz3C5sEycOJE5c+YwY8YM4uLiyM3NJSMjwzmxdu/evRw8eNC5/PDhw3n99dd58cUXGTx4MG+//TZLly5l4MCBAPj7+7Np0yZ+/vOfc/HFF3PnnXcSHx/P559/TlBQUCO9TZHzm5H5IiEnKviySz8+6DvKdBwRt/XWoGQ2R/QktPIYD6x82XQc8RE2u93zp3qXlZURFhZGaWmp5rOIy2KmfsjY3ev5x9uPctLmx7V3/IXtnbqbjiXi1uLzt/Kv1x6iGht+q7+wjtUi4iJXPr91LiHxecEnjvPHZQsA+PuQ61VWROohp2t/3hqYjB92mDwZKitNRxIvp8IiPu9/V79JdGkh+9t05C8jbjEdR8RjPHH5f1McEgZbtsCTT5qOI15OhUV8W04Od619B4CZV/6aY4EtDQcS8RwlLUP54xV3WXcee0znGZImpcIivuv4cZg0iQB7NR/0GcmnvRNNJxLxOO/1uwzGj7eGhO680zqrs0gTUGER3zVjBmzdSlGrtky/aorpNCKeyWaD+fOtI0OvXg2zZ5tOJF5KhUV806pVMGcOANNSfsMPIWGGA4l4sG7dYN486/bMmdapLUQamQqL+J7vv4dbb7VO3nbHHRoKEmkMt90GN90EJ09at48dM51IvIwKi/gWux1SU2HvXujZE+bONZ1IxDvYbLBgAURFwY4dUHM0c5HG4hZnaxZpNnPnwnvvQWAg/POfEKahIJELdfr5tpIuu5dXlzyC/z/+wUMH2/DPwVcBOqOzXDhtYRHfsXo1PPSQdfvPf4ZLLzWbR8QLZXeL5elRtwHw2LL5DCzYZTiReAttYRGv5vjLr0vpIZa+nEbHkyf5oO8o7v3uItBZmEWaxPxhN3LJgR1cuWst85em8/NJz5iOJF5AW1jE64VU/sjCdx6j47EStnbqzoPj77PG20WkSdhtfvz2mv8jr21noksLefGdJ6zjHolcABUW8Wp+1VX8+YOn6X9oD0Wt2vI/N0znx8Bg07FEvF5ZcGvuvGEGZUGtGLJ/K/z3f1uT3kUaSIVFvJfdzmPL5pOycw0V/gH8asLvORDayXQqEZ+xOzyaX014mBN+/vDGGzB9uulI4sFUWMR7TZ/OrbkZVGPj/659gA1d+plOJOJzsrsN5uGUml2cn3jCecBGEVepsIh3evpp65cj8EjK3XzUd6ThQCK+663YK+FPf7LuPPigdSh/ERepsIj3efJJeOAB6+Zlk3g9brzhQCLCtGnw8MPW7bvvhr//3Wwe8TjarVm8y2OPWSc1BJgxg78eH2I2j4ic8vjjUF4Of/kL/M//wJEjcP/9plOJh9AWFvEO1dXw29+eKitPPAGPPqrdl0XcRMzUD4mZ9hExQcm8MPS/rAf/7//4y4j/R8zvPjAbTjyCtrCIR4uZ+iFBJyr48wdPc/U3qwF4bOyd/L1ssA4MJ+KObDbSx6RSGtyah1a+zH2r36DzkSJ47CrrlBkiZ6EtLOLROh05zBtvPszV36ymwj+A//3Zg/x96ATTsUTkXGw2/pr0Cx656m6qbH78YvOnkJwMxcWmk4kbU2ERz5WVxYcv3celB3ZQEtyaX058nPf6jzadSkTq6dVLrua/b5xJWWAIfP45DBkC69aZjiVuSoVFPE9VFaSnwxVX0PFYCds6xnDdpGdYFz3QdDIRcdGKHvH81y/nQI8ekJcHI0daJyfVUXHlJ1RYxLPs3g2jR1u7R1ZX868BY5nwyzl81y7KdDIRaaBd4RdBTg7ccAOcOAFpaTB+POzbZzqauBEVFvEMJ09au0IOHgxffAFt2sDf/85vr0njeAudG0jE47VtC2+9BfPmQVAQfPwxDBgAL75o7QUoPk+FRdxfdjYkJFjHaygvhzFjYNMm62Rq2m1ZxCs4d3ve240rbpvLhqg+1nFafvUrNkT35+e3/9l0RDHMZrd7/kBhWVkZYWFhlJaWEhoaajqONJKxd73Ib1e+wrU7VgFQEtyaJ0ffzhuDU7Db1LVFvJlfdRWpOe+T9vmrtDpx3HrwjjusYy117240mzQeVz6/VVjE/XzzDTz1FCf/vogAezXV2HhrUDKzx9zB9yFhptOJSDPqdOQwv1u5mBu+/sx6ICDAKi4PP6zi4gVUWMQzrV4NTz0F777r3EPg055DmHPZJLZ30i8mEV8Wd2AHSws/hk8+sR4ICIBbboF774WhQ82GkwZTYRHP8cMP8PrrsGgRbNhw6vGf/YwbQkeR07W/uWwi4lbyZl1j/WHz6KOnigtYc9zuuQduvBFatzYXUFymwiLu7ehRptz2BFfv+IKrdq4hqOoEABX+ASztP5YXh/4Xu8OjDYcUEXc2+MAOJm38kGu3rSSo6iQAPwYE0fKG660tL+PGWXsbiVtTYRH3YrfDrl2wfDl88IH1l1FFhfPpbR1jWBJ7Fe/2H80PmqMiIi5of6yUiZs+4RebPqH7DwdPPdG6tXW4/6uvto7p0rWruZByViosYtbJk7B9O6xfb5WU5cshP7/WInltO5NxcRIf9LuMryN6avdkEbkwdjuDCnbxfoe9sGQJHDhQ+/levWDECOtIuiNGQN+++r3jBpq8sMybN4+nnnqKgoICBg8ezHPPPcfQc0x6euutt5g+fTp5eXn07t2b2bNnc/XVVzuft9vtzJw5k4ULF1JSUsKIESOYP38+vXv3rlceFRZD7HYoKIAdO6w9ezZtsuah5ObCjz/WWrTCP4DcqL580W0wH1+cxI7wbvplISJNwmavZkDht4zdvZ6x335J3IFv8OMnH3VhYRAbe+oyaBD06QPt2ul3UzNq0sKyZMkSJk2axIIFC0hMTGTu3Lm89dZb7Nixg06dOp2x/OrVq7nssstIT0/n2muv5fXXX2f27Nls2LCBgQOtc7/Mnj2b9PR0Fi9eTPfu3Zk+fTqbN29m69atBAef/yimKixN5PhxOHQI9u+3tpA4rvPzYedOq6QcPVrnS48GtmRLRE/Wd+1P9kWx5HTpqyPSiogRocePcun+bSTs30ZC/laGFe8+448qp7Aw67xGPXta1926QWRk7UtISPO+AS/WpIUlMTGRIUOG8PzzzwNQXV1NdHQ0v/nNb5g6deoZy0+cOJHy8nI++OAD52PDhg0jLi6OBQsWYLfbiYqK4re//S0PPPAAAKWlpURERPDSSy9x8803N+ob9ip2u3XeDcelsrL2/RMnrNJRXn7uS1kZHD5sXb7//tTts/2HPk2VzY99YRHsaR/Frg7RbI7sxdcRvdjTPkoHdxMRtxRQdZJeh/fRtyiPvof20K8ojz5FeUQe/b5+K2jTBiIirNMJtG1rlZzTb4eFQatW0LLluS+Bgdbu2QEB0KKFde3v71NbeFz5/A5wZcWVlZXk5OQwbdo052N+fn4kJyeTnZ1d52uys7NJS0ur9VhKSgpLly4FYM+ePRQUFJCcnOx8PiwsjMTERLKzs+ssLBUVFVScNmmztLQUsN54ozpxwpqsZbefOnOo4/aFXBq6npMna5eRZji/xgmbH4dat6ewdXsOte5AQZv2FLZqT37bSPLadia/bQQn/Fuc+cLK402eTUSkISqBrWERbA2LgF6JzseDThynS1kR0aWFdC0tJLq0kIgjhwkvLyH8WAkdy0sIrjphnTLgyJGmC+jvf6rInH7x9wc/P6vQnH7xq/nj8KeP1bXc2V57up8WJsf9wED46KNGfauOz+36bDtxqbAUFxdTVVVFRERErccjIiLYvn17na8pKCioc/mCggLn847HzrbMT6Wnp/Poo4+e8Xh0tHaFbXT2ajhSbF1ERLzcLtMBAKqqrMtpf5i7jbCm2ZPzyJEjhJ1n3S4VFncxbdq0Wlttqqur+f777+nQoQM2L9qUVlZWRnR0NPv27fOtoS43ou+BefoeuAd9H8zzxu+B3W7nyJEjREVFnXdZlwpLeHg4/v7+FBYW1nq8sLCQyMjIOl8TGRl5zuUd14WFhXTu3LnWMnFxcXWuMygoiKCfHBCobdu2rrwVjxIaGuo1P5yeSt8D8/Q9cA/6Ppjnbd+D821ZcXBpVmRgYCDx8fFkZmY6H6uuriYzM5OkpKQ6X5OUlFRreYBly5Y5l+/evTuRkZG1likrK2Pt2rVnXaeIiIj4FpeHhNLS0rj99ttJSEhg6NChzJ07l/LyclJTUwGYNGkSXbp0IT09HYD77ruP0aNH8/TTT3PNNdfw5ptv8uWXX/Liiy8CYLPZuP/++3n88cfp3bu3c7fmqKgorr/++sZ7pyIiIuKxXC4sEydOpKioiBkzZlBQUEBcXBwZGRnOSbN79+7F77RZx8OHD+f111/nkUce4eGHH6Z3794sXbrUeQwWgIceeojy8nLuuusuSkpKGDlyJBkZGfU6Bos3CwoKYubMmWcMf0nz0ffAPH0P3IO+D+b5+vfAKw7NLyIiIt5NR/YSERERt6fCIiIiIm5PhUVERETcngqLiIiIuD0VFg9UUVFBXFwcNpuN3Nxc03F8Rl5eHnfeeSfdu3enZcuW9OzZk5kzZ1JZWWk6mlebN28eMTExBAcHk5iYyLp160xH8hnp6ekMGTKENm3a0KlTJ66//np27NhhOpZPmzVrlvNwIL5GhcUDPfTQQ/U6jLE0ru3bt1NdXc0LL7zAli1b+POf/8yCBQt4+OGHTUfzWkuWLCEtLY2ZM2eyYcMGBg8eTEpKCocOHTIdzSesWLGCe+65hzVr1rBs2TJOnDjBVVddRXl5ueloPmn9+vW88MILxMbGmo5ihHZr9jD/+c9/SEtL41//+hcDBgxg48aNZz2FgTS9p556ivnz5/Ptt9+ajuKVEhMTGTJkCM8//zxgHVk7Ojqa3/zmN0ydOtVwOt9TVFREp06dWLFiBZdddpnpOD7l6NGjXHrppfz1r3/l8ccfJy4ujrlz55qO1ay0hcWDFBYWMnnyZF555RVCQkJMxxGgtLSU9u3bm47hlSorK8nJySE5Odn5mJ+fH8nJyWRnZxtM5rtKS0sB9DNvwD333MM111xT6/+Dr/HIszX7Irvdzh133MGvf/1rEhISyMvLMx3J5+3atYvnnnuOOXPmmI7ilYqLi6mqqnIeRdshIiKC7du3G0rlu6qrq7n//vsZMWJErSOVS9N788032bBhA+vXrzcdxShtYTFs6tSp2Gy2c162b9/Oc889x5EjR5g2bZrpyF6nvt+D0+3fv59x48Zx0003MXnyZEPJRZrPPffcw9dff82bb75pOopP2bdvH/fddx+vvfaaz5+uRnNYDCsqKuLw4cPnXKZHjx784he/4P3338dmszkfr6qqwt/fn1tvvZXFixc3dVSvVd/vQWBgIAAHDhxgzJgxDBs2jJdeeqnWubOk8VRWVhISEsLbb79d60Sot99+OyUlJbz77rvmwvmYe++9l3fffZeVK1fSvXt303F8ytKlS5kwYQL+/v7Ox6qqqrDZbPj5+VFRUVHrOW+mwuIh9u7dS1lZmfP+gQMHSElJ4e233yYxMZGuXbsaTOc79u/fz9ixY4mPj+fVV1/1mV8UpiQmJjJ06FCee+45wBqWuOiii7j33ns16bYZ2O12fvOb3/Dvf/+brKwsevfubTqSzzly5AjfffddrcdSU1Pp27cvv/vd73xqeE5zWDzERRddVOt+69atAejZs6fKSjPZv38/Y8aMoVu3bsyZM4eioiLnc5GRkQaTea+0tDRuv/12EhISGDp0KHPnzqW8vJzU1FTT0XzCPffcw+uvv867775LmzZtKCgoACAsLIyWLVsaTucb2rRpc0YpadWqFR06dPCpsgIqLCL1tmzZMnbt2sWuXbvOKInaUNk0Jk6cSFFRETNmzKCgoIC4uDgyMjLOmIgrTWP+/PkAjBkzptbj//jHP7jjjjuaP5D4NA0JiYiIiNvTbEERERFxeyosIiIi4vZUWERERMTtqbCIiIiI21NhEREREbenwiIiIiJuT4VFRERE3J4Ki4iIiLg9FRYRERFxeyosIiIi4vZUWERERMTtqbCIiIiI2/v/XD9BqNtCLuMAAAAASUVORK5CYII=\n"},"metadata":{}}]},{"cell_type":"markdown","source":["## 2. The Ziggurat algorithm\n","\n","This is adapted from Wikipedia's page on the [Ziggurat algorithm](https://en.wikipedia.org/wiki/Ziggurat_algorithm), which was published in 2000.\n","\n","The algorithm primarily generates a positive random variable $X'$ distributed as $|X|$ with $X \\sim \\mathcal{N}(0,1)$, so that $X'$ has density\n","$$\\sqrt{\\frac{2}{\\pi}} \\exp\\left(-\\frac{x^2}{2}\\right)$$\n","on $[0,+\\infty)$.\n","\n","To generate a realisation of $X$ from a realisation of $X'$ you may use the following function."],"metadata":{"id":"zpKIMSDyPWoK"}},{"cell_type":"code","source":["def RandomSign(x):\n","  return (1-2*(Unif()<.5))*x"],"metadata":{"id":"EoF3c8Jo6-O7","executionInfo":{"status":"ok","timestamp":1697484078940,"user_tz":-120,"elapsed":235,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":7,"outputs":[]},{"cell_type":"markdown","source":["The Ziggurat algorithm does not depend on the density of $X'$ being normalised so we set\n","$$f(x) = \\exp\\left(-\\frac{x^2}{2}\\right).$$"],"metadata":{"id":"R_RHpB327eOr"}},{"cell_type":"code","source":["def f(x):\n","  return np.exp(-x**2/2)\n","\n","y_max = f(0)\n","\n","def f_inv(y): # inverse function\n","  return (-2*np.log(y))**.5"],"metadata":{"id":"fizYUEu1_OUq","executionInfo":{"status":"ok","timestamp":1697484080715,"user_tz":-120,"elapsed":15,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":8,"outputs":[]},{"cell_type":"markdown","source":["### 1. Precomputing the ziggurat\n","\n","It may help to have in mind what a [ziggurat](https://en.wikipedia.org/wiki/Ziggurat) actually is.\n","\n","The basis of the Ziggurat algorithm consists in covering the graph of $f$ with a predefined number $n$ of horizontal layers defined as follows.\n","\n","The base is defined as all points inside the distribution and below $y_1 = f(x_1)$. This consists of a rectangular region from $(0, 0)$ to $(x_1, y_1)$, and the (infinite) tail of the distribution, where $x > x_1$ (and $y < y_1$).\n","\n","This layer (call it layer $0$) has area $A$. On top of this, add a rectangular layer of width $x_1$ and height $A/x_1$, so it also has area $A$. The top of this layer is at height $y_2 = y_1 + A/x_1$, and intersects the density function at a point $(x_2, y_2)$, where $y_2 = f(x_2)$. This layer includes every point in the density function between $y_1$ and $y_2$, but (unlike the base layer) also includes points such as $(x_1, y_2)$ which are not in the distribution. Further layers are then stacked on top, until $y_n > y_\\mathrm{max} = f(0)$ in which case we set $x_n=0$.\n","\n","The following function takes $x_1$ as an argument and returns a pair $(x,y)$ with $x$ and $y$ are two lists which respectively contain $(x_1, \\ldots, x_n)$ and $(y_1, \\ldots, y_n)$."],"metadata":{"id":"1rW988tY_O2X"}},{"cell_type":"code","source":["def layers(x1):\n","  y1 = f(x1)\n","  A = x1*y1 + (2*math.pi)**.5*(1-norm.cdf(x1))\n","\n","  x_c = x1 #current x_i\n","  y_c = y1 #current y_i\n","  y_n = y_c + A/x_c #next y_i\n","\n","  x = [x_c]\n","  y = [y_c]\n","  while (y_n < y_max):\n","    y_c = y_n\n","    x_c = f_inv(y_c)\n","    x.append(x_c)\n","    y.append(y_c)\n","\n","    y_n = y_c + A/x_c\n","\n","  x.append(0)\n","  y.append(y_n)\n","\n","  return [x,y]"],"metadata":{"id":"dDPIJXz17dl6","executionInfo":{"status":"ok","timestamp":1697484458856,"user_tz":-120,"elapsed":232,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":9,"outputs":[]},{"cell_type":"markdown","source":["The covering of $f$ by the ziggurat is visualised as follows."],"metadata":{"id":"RvgsjTDPFW3M"}},{"cell_type":"code","source":["x1 = 2.5\n","plt.plot(np.linspace(0,1.2*x1,100),f(np.linspace(0,1.2*x1,100)),color='red')\n","\n","[x,y] = layers(x1)\n","n = len(x)\n","\n","\n","plt.plot([0,1.2*x1],[0,0],color='blue')\n","plt.plot([0,x[0]],[y[0],y[0]],color='blue')\n","plt.plot(np.linspace(x1,1.2*x1,20),f(np.linspace(x1,1.2*x1,20)),color='blue')\n","\n","for i in range(n-1):\n","  plt.plot([x[i],x[i]],[y[i],y[i+1]],color='blue')\n","  plt.plot([0,x[i]],[y[i+1],y[i+1]],color='blue')\n","\n","print(\"Number of layers:\", n)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":448},"id":"1e8n6102FuWC","executionInfo":{"status":"ok","timestamp":1697486732112,"user_tz":-120,"elapsed":331,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"b0bf5014-1864-4c14-d186-2efe20e08d0f"},"execution_count":22,"outputs":[{"output_type":"stream","name":"stdout","text":["Number of layers: 11\n"]},{"output_type":"display_data","data":{"text/plain":["<Figure size 640x480 with 1 Axes>"],"image/png":"iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA77UlEQVR4nO3deVxVdf7H8dcFBEQFNWRRKTOzzVIjJWzToijbtKm0JjVbZmy0n8pUapmOOROZtmc52WY15lZaWVlmallWE8qk5pZLWgpqJchFIeH8/vgmRIFyiXu/d3k/H4/z6HA4554Pp+vlzTnfxeU4joOIiIiIJWG2CxAREZHQpjAiIiIiVimMiIiIiFUKIyIiImKVwoiIiIhYpTAiIiIiVimMiIiIiFUKIyIiImJVhO0CaqO8vJwdO3bQpEkTXC6X7XJERESkFhzHYd++fbRs2ZKwsJrvfwREGNmxYwcpKSm2yxAREZE62L59O61bt67x+wERRpo0aQKYHyY2NtZyNSIiIlIbhYWFpKSkVPwer0lAhJFDj2ZiY2MVRkRERALMkZpYqAGriIiIWKUwIiIiIlYpjIiIiIhVCiMiIiJilcKIiIiIWKUwIiIiIlYpjIiIiIhVCiMiIiJilcKIiIiIWKUwIiIiIlYpjIiIiIhVATE3jTc4DhQX264itP36/0F8PBxmdmkREQliHoeRjz76iIkTJ5KTk8POnTuZO3cuvXr1OuwxS5YsISsrizVr1pCSksLo0aO58cYb61hy/SguhsaNrZYgv5KfDwkJtqsQEREbPP5b1O1207FjRyZPnlyr/bds2cKll15Kjx49yM3NZdiwYdxyyy289957HhcrIiIiwcfjOyOXXHIJl1xySa33nzJlCsceeywPPfQQACeddBLLli3jkUceITMz09PT15uYGCjalG++aNgQoqOhQQM4wjTHUn/cbkhMNOsxMXZrERERe7zeZmT58uVkZGRU2ZaZmcmwYcNqPKakpISSkpKKrwsLC+u9LpcLGg24GpYtq9wYFmaCyVFHQcuWZklOhtat4cQT4ZRToG1bCA+v93pCnTKgiEjo8noYycvLI/HQn7+/SExMpLCwkP3799OwYcPfHZOdnc24ceO8XZppQflr5eXmz3W3G7Ztq/6Y6GgTTDp3hrPPNsvxx+u3qYiISB35ZW+aUaNGkZWVVfF1YWEhKSkp9X+iZctMICkpgQMHYP9+07J1zx7YscMsO3fC1q3w9dewdq3ZLzfXLC+8YF6nRQsTSi6+GC67zNxRERERkVrxehhJSkoiPz+/yrb8/HxiY2OrvSsCEBUVRVRUlLdLM1wuc7cjOhqaNjXbjjuu+n3LymDLFli9Gj7/3ISZL76A3bth7lyzAKSmmlDSqxd07Ki7JiIiIofh9TCSnp7OO++8U2XbwoULSU9P9/ap6194OLRrZ5ZD3ZkPHICcHFi8GObPN+EkJ8cs48ZBhw7Qrx/8+c/QqpXV8kVERPyRx117i4qKyM3NJTc3FzBdd3Nzc9n2SxuLUaNG0b9//4r9Bw0axObNm7nrrrtYt24dTz31FLNmzWL48OH18xPYFh0NZ50Fo0fDZ5+ZxzrPPQdXXgmRkeYuyogRkJICF14Is2fDwYO2qxYREfEbLsf5bSvOw1uyZAk9evT43fYBAwbw4osvcuONN7J161aWLFlS5Zjhw4fz9ddf07p1a+69916PBj0rLCwkLi6OgoICYmNjPSnXrp9+MuHj5Zer9tpp3Rr+9je49VYz9GiIcrsrB54rKoJGjezWIyIi9au2v789DiM2BGwY+bXNm02D13//27QxAYiKMo9vRo40PXJCjMKIiEhwq+3vb80G4itt28L48abL8LRpppFrSQk8/7zpKtyvH6xbZ7tKERERn1MY8bXoaOjfH/77X/jkE9PrprwcXnkFTj4Z+vY1XYhFRERChMKILS4XdOsGb70FX35peuc4DsycCaeeCoMGmdnjREREgpzCiD9ITTVjlOTmml44ZWWmbUm7dubRjtttu0IRERGvURjxJx07wrx5sHQpdOliWnWOGQPt25s7Jv7f1lhERMRjCiP+6NxzzZglr74KbdqYYen79jXDzW/caLs6ERGReqUw4q/Cwiobs/7jH6Yb8Pvvm/Yk//iHGflVREQkCCiM+LvoaBg7FlatgosuMt2Bx42DTp3M3RMREZEApzASKI4/HhYsgFmzIDkZ1q83w9CPGKG7JCIiEtAURgKJywXXXANr1phB0srL4cEH4fTTzbglIiIiAUhhJBA1awYvvQRvvAFJSaZdSXo63Hef6RYsIiISQBRGAtkVV5hZga+7zoSQsWPh/PPhu+9sVyYiIlJrCiOB7qijYPp0MzNw48bw0UeV45WIiIgEAIWRYHHDDbByJZxxBvz4I/TuDUOGQGmp7cpEREQOS2EkmLRrZybfu/NO8/XkyXDeeXpsIyIifk1hJNhERpoeNm+/DU2bmrFITj8dFi+2XZmIiEi1FEaCVc+ekJNjBkfbvRsyMmDiRM1vIyIifkdhJJi1bWse2wwYYMYkuesuuP562L/fdmUiIiIVFEaCXUwMvPACPP00RETAjBmmHcnOnbYrExERARRGQoPLBYMGwcKF0Ly5Ga21a1fT+0ZERMQyhZFQ0r07fP45nHii6WFz9tnw+uu2qxIRkRCnMBJq2rUzPWwyM6G4GK6+Gh591HZVIiISwhRGQlFcHMyfD4MHm941w4dDVpZp5CoiIuJjCiOhKiICnngCJkwwXz/yCPTtCwcO2K1LRERCjsJIKHO5THff//wHGjSA2bPhoovgp59sVyYiIiFEYUTM2CMLFkBsLHz8MZx7rrr+ioiIzyiMiHH++SaIJCfD6tVwzjmwZYvtqkREJAQojEil006DZcvg2GNh0yY46yxYs8Z2VSIiEuQURqSqtm1NIDnlFPOo5txzzSBpIiIiXqIwIr/XsiUsXWpGaf3xx8pHOCIiIl6gMCLVO+ooWLTIBJGiIrj4Yli82HZVIiIShBRGpGaNG5vB0Q6N1tqzp5nfRkREpB4pjMjhNWwI8+bBpZeaAdEuvxzefdd2VSIiEkQURuTIoqPNhHpXXgklJdCrF7z1lu2qREQkSCiMSO1ERpoRWq++GkpL4U9/gnfesV2ViIgEAYURqb0GDeDVV00g+flnuOoqtSEREZE/TGFEPBMRAdOnVz6yueIK9bIREZE/RGFEPNegAcycWdmo9bLLNA6JiIjUmcKI1E1UFMyZY2b5PdTt9/PPbVclIiIBSGFE6i462nT7PTQw2iWXwFdf2a5KREQCjMKI/DENG8Ibb0B6Ovz0k7lTsnGj7apERCSAKIzIH9e4Mbz9NnTsCPn5kJEB27fbrkpERAKEwojUj2bN4L33oH172LYNLrwQdu2yXZWIiAQAhRGpP4mJZtyRlBRYv95MrldYaLsqERHxcwojUr+OPho++ABatICVK6F3bzMeiYiISA0URqT+tW9vJtNr3Bg+/BBuuAHKymxXJSIifkphRLwjNRXmzjUDpM2ZA0OHguPYrkpERPyQwoh4T0YGvPIKuFwweTL885+2KxIRET+kMCLede218PjjZn3MGHj2Wbv1iIiI31EYEe8bMgTuucesDxpk2pOIiIj8QmFEfGP8eOjf3zRkveYayMmxXZGIiPgJhRHxDZcLpk417UjcbjPj77ff2q5KRET8gMKI+E5kJLz2Gpx2mhk2vndv2xWJiIgfiLBdgC2OY2a+Fx8Lj4U570L37rg3VM5f43ZbrOlXfv2+iI+HMMV1ERGvC9kwUlxsxuQSG1oCG6psSUy0U8nh5OdDQoLtKkREgl+d/u6bPHkybdq0ITo6mrS0NL744ovD7v/oo49ywgkn0LBhQ1JSUhg+fDgHDhyoU8EiIiISXDy+MzJz5kyysrKYMmUKaWlpPProo2RmZrJ+/XoSqvkzcvr06YwcOZLnn3+ebt26sWHDBm688UZcLhcPP/xwvfwQdRETA0VF1k4v/PJI5JlX4O/DiaEY17PPQd++Vmtyuyvv0sTEWC1FRCRkuBzHszG609LS6NKlC08++SQA5eXlpKSkcPvttzNy5Mjf7T9kyBDWrl3LokWLKrb9/e9/5/PPP2fZsmW1OmdhYSFxcXEUFBQQGxvrSbkSCEaMgAcfNA1cP/gAzjnHWilud+Xju6IiaNTIWikiIgGvtr+/PXpMU1paSk5ODhkZGZUvEBZGRkYGy5cvr/aYbt26kZOTU/EoZ/Pmzbzzzjv07NnTk1NLMMvOhj/9CUpLTQ+bb76xXZGIiPiQR49p9uzZQ1lZGYm/aW2YmJjIunXrqj3m+uuvZ8+ePZx99tk4jsPBgwcZNGgQd999d43nKSkpoeRX084XFhZ6UqYEmrAweOkl2L4dvvgCLr8cli+Hpk1tVyYiIj7g9Y6LS5Ys4f777+epp55ixYoVvP7667z99tuMHz++xmOys7OJi4urWFJSUrxdptgWEwNvvAGtW8O6ddCnDxw8aLsqERHxAY/CSHx8POHh4eTn51fZnp+fT1JSUrXH3HvvvfTr149bbrmFU089ld69e3P//feTnZ1NeXl5tceMGjWKgoKCimX79u3V7idBJikJ3nzTBJP334esLNsViYiID3gURiIjI0lNTa3SGLW8vJxFixaRnp5e7THFxcWE/WbkqPDwcABqajsbFRVFbGxslUVCROfO8PLLZv2JJ2DKFLv1iIiI13n8mCYrK4upU6cybdo01q5dy2233Ybb7WbgwIEA9O/fn1GjRlXsf/nll/P0008zY8YMtmzZwsKFC7n33nu5/PLLK0KJSBVXXQX/+pdZHzIEPvzQbj0iIuJVHo8z0qdPH3bv3s2YMWPIy8ujU6dOLFiwoKJR67Zt26rcCRk9ejQul4vRo0fz/fff06JFCy6//HL+deiXjUh1Ro2Cr7+G//wHrr7aNGxt1852VSIi4gUejzNig8YZCVEHDkD37vD553DyyaaHjZf//2ucERGR+uOVcUZEfCo6Gl5/HZKTzV2SG26AGho9i4hI4FIYEf/WsiXMmwdRUfDWWzB2rO2KRESknimMiP/r2hWmTjXr//wnzJ5ttx4REalXCiMSGPr1g7//3awPGAC5uVbLERGR+qMwIoFjwgTIzIT9+6FXL9izx3ZFIiJSDxRGJHCEh8Orr5ouvt9+qyHjRUSChMKIBJZmzWDuXNPn9sMPYeRI2xWJiMgfpDAigadDB3jxRbP+0EMwfbrVckRE5I9RGJHAdPXVZpRWgFtuUYNWEZEApjAigWv8eLj4YtOgtXdv+OEH2xWJiEgdKIxI4AoPN49ojjsOtm6F66+HsjLbVYmIiIcURiSwNWtmhoyPiYH339cIrSIiAUhhRALfaafBs8+a9X/9C954w249IiLiEYURCQ7XXQdDh5r1/v1hwwa79YiISK0pjEjwmDgRzjkHCgtNg9aiItsViYhILSiMSPBo0ABmzYLkZPj6a7j5ZnAc21WJiMgRKIxIcElKgjlzICLCBJPHH7ddkYiIHIHCiASfbt3g4YfN+h13wCef2K1HREQOS2FEgtOQIaZR68GDcO21kJ9vuyIREamBwogEJ5cLnnkGTj4ZduyAvn01w6+IiJ9SGJHg1bgxvPaa+e+SJTB6tO2KRESkGgojEtxOPBFeeMGsT5gA8+ZZLUdERH5PYUSC39VXw/DhZv3GG2HTJqvliIhIVQojEhomTDC9bAoK4Jpr4MAB2xWJiMgvFEYkNDRoADNnQnw8rFxZOXS8iIhYpzAioaN1a5g+vbKnzUsv2a5IRERQGJFQc+GF8I9/mPVBg2DVKqvliIiIwoiEotGjITMT9u83jVv37bNdkYhISFMYkdATFgavvGIe22zYAH/5iybUExGxSGFEQlN8vGnQGhEBM2bAlCm2KxIRCVkKIxK6unWDBx4w68OGQU6O1XJEREKVwoiEtqwsuPJKKC0144/s3Wu7IhGRkKMwIqHN5TLDxbdpA1u2mB42IiLiUwojIs2awezZEBkJ89+yXY2ISMiJsF2ALY4DxcW2qxC/cdIZkP047r9Xzuzrdvvm1L9+L8bHm84+IiKhJGTDSHGxmVlepNJff1mMxETfV5CfDwkJvj+viIhN+htMRERErArZOyMxMVBUZLsK8TeOA8WffQWXXELMwb24/nW/1yfVc7sr78LExHj1VCIifilkw4jLBY0a2a5C/FHjjNPg8dHwt7/B2DugR1dIT/fJuV0un5xGRMSv6DGNSHUGDYI+feDgQfPfH36wXZGISNBSGBGpjssFzzwD7drB9u0wcKDmrxER8RKFEZGaxMbCrFlm/JG33oJHHrFdkYhIUFIYETmczp3h0UfN+ogR8NlnVssREQlGCiMiRzJokJm35uBB6NsXfvrJdkUiIkFFYUTkSFwumDoV2raFb79V+xERkXqmMCJSG3Fxle1H3ngDHnvMdkUiIkFDYUSktlJT4aGHzPpdd8F//2u3HhGRIKEwIuKJwYPhqqvg55/N+CMFBbYrEhEJeAojIp5wueC556BNG9iyBW65Re1HRET+IIUREU81bQozZ0JEBMyZA1Om2K5IRCSgKYyI1EXXrjBhglkfPhxyc62WIyISyBRGROpq+HC4/HIoKYFrr4V9+2xXJCISkBRGROrK5YIXX4SUFNi40QyOpvYjIiIeUxgR+SOaN4dXX4XwcJg+HZ5/3nZFIiIBR2FE5I866yz45z/N+u23w+rVdusREQkwCiMi9eGuuyAzE/bvN+OPuN22KxIRCRh1CiOTJ0+mTZs2REdHk5aWxhdffHHY/ffu3cvgwYNJTk4mKiqK9u3b884779SpYBG/FBYGL70Eycnw9dfwf/9nuyIRkYDhcRiZOXMmWVlZjB07lhUrVtCxY0cyMzPZtWtXtfuXlpZy4YUXsnXrVubMmcP69euZOnUqrVq1+sPFi/iVhATTbiQszLQdeeUV2xWJiAQEl+N41vw/LS2NLl268OSTTwJQXl5OSkoKt99+OyNHjvzd/lOmTGHixImsW7eOBg0a1KnIwsJC4uLiKCgoIDY2tk6vIeIz48bBP/4BjRpBTg6ccMJhd3e7oXFjs15UZA4TEQkGtf397dGdkdLSUnJycsjIyKh8gbAwMjIyWL58ebXHvPnmm6SnpzN48GASExPp0KED999/P2VlZTWep6SkhMLCwiqLSMAYPRp69DAp49prTTsSERGpkUdhZM+ePZSVlZGYmFhle2JiInl5edUes3nzZubMmUNZWRnvvPMO9957Lw899BD/PNT7oBrZ2dnExcVVLCkpKZ6UKWJXeDj85z/QogV89RVkZdmuSETEr3m9N015eTkJCQk888wzpKam0qdPH+655x6mHGY+j1GjRlFQUFCxbN++3dtlitSv5OTKNiNTpsCsWXbrERHxYx6Fkfj4eMLDw8nPz6+yPT8/n6SkpGqPSU5Opn379oSHh1dsO+mkk8jLy6O0tLTaY6KiooiNja2yiASciy6CUaPM+i23wKZNdusREfFTHoWRyMhIUlNTWbRoUcW28vJyFi1aRHp6erXHnHXWWXzzzTeUl5dXbNuwYQPJyclERkbWsWyRAHHffWZQtH37zPgjJSW2KxIR8TseP6bJyspi6tSpTJs2jbVr13LbbbfhdrsZOHAgAP3792fUob8Ggdtuu40ff/yRoUOHsmHDBt5++23uv/9+Bg8eXH8/hYi/iogww8U3b2561tx1l+2KRET8ToSnB/Tp04fdu3czZswY8vLy6NSpEwsWLKho1Lpt2zbCwiozTkpKCu+99x7Dhw/ntNNOo1WrVgwdOpQRI0bU308h4s9SUmDaNDPD7+OPm542vXrZrkpExG94PM6IDRpnRILCHXfAQw9B06awciW0aQNonBERCV5eGWdERP6A+++Hrl1h717o2xdqaMAtIhJqFEZEfCUyEmbONHdGPv8c7rnHdkUiIn5BYUTEl9q0MfPWAEyaBPPnWy1HRMQfKIyI+Frv3nD77WZ9wAD47ju79YiIWKYwImLDxImQmgo//gg33mi7GhERqxRGRGyIijLtR2Jj4bPqJ5kUEQkVHo8zEiwcB4qLbVchIS3pOJj8Iu5+f6nY5HZ773S/fs/Hx0OY/hQRET8RsmGkuLhybAcRe3r/shi/mRDba/LzISHBN+cSETkS/W0kIiIiVoXsnZGYGDPapYhtjgPFqzdDRgYx7nxcd94FY8fW+3nc7so7LzEx9f7yIiJ1FrJhxOXSsNviPxqf2Raef8DM7DtpHFxwJmRmeu18LpfXXlpExGN6TCPiL669Fm67zdwqueEG+P572xWJiPiEwoiIP3n4YejcGfbsgeuug4MHbVckIuJ1CiMi/iQ6GmbNgiZN4OOPYcwY2xWJiHidwoiIv2nXDp591qxnZ8O779qtR0TEyxRGRPzRtdfC3/5m1vv1g+3b7dYjIuJFCiMi/uqhh+D00+GHH6BvX/j5Z9sViYh4hcKIiL861H4kNhY+/RTuucd2RSIiXqEwIuLPjjsOXnjBrE+cCG++abceEREvUBgR8XdXXQXDhpn1AQNg61ab1YiI1DuFEZFAMGECpKXB3r2mcWtJie2KRETqjcKISCCIjISZM6FZM/jvf+HOO21XJCJSbxRGRALFMcfAyy+b9SeeMI1bRUSCgMKISCC59FIYOdKs33wzrF9vtx4RkXqgMCISaMaPh/POg6IiuPpqKC62XZGIyB+iMCISaCIi4NVXITERVq82I7U6ju2qRETqTGFEJBAlJ8OMGRAWBtOmwfPP265IRKTOFEZEAlX37uaRDcCQIbBypdVyRETqSmFEJJCNHGkatR44YNqP7N1ruyIREY8pjIgEsrAweOklaNMGNm82I7SWl9uuSkTEIwojIoGueXOYM8cMjPbmm2YOGxGRAKIwIhIMUlPNQGgAd98NS5ZYLUdExBMKIyLB4tZbKx/T9OkDO3bYrkhEpFYURkSChcsFTz0Fp50Gu3aZCfV+/tl2VSIiR6QwIhJMYmJM+5HYWPjkE02oJyIBQWFEJNgcf7zpYQPw2GNmtFYRET+mMCISjK68EkaNMuu33AJr1titR0TkMBRGRILV+PGQkWEm0rvqKigosF2RiEi1FEZEglV4OEyfDikpsGED/PWvtisSEalWhO0CbHEczbwuISCmBbz8Olx4Ie75iyo2u931f6pf/5uKjzeDw4qI1EbIhpHiYmjc2HYVIr5wBvBTlS2Jid49Y34+JCR49xwiEjz0t4uIiIhYFbJ3RmJioKjIdhUivuPsP0BxxhXwv5XEdGiH68NF5h9CPXG7K++41OPLikgICNkw4nJBo0a2qxDxoUbRNH7rOTOPzerP4O+DYNo084+hnnnhJUUkiOkxjUgoSUmBmTNN69KXX4Ynn7RdkYiIwohIyOnRAx580KxnZcHSpXbrEZGQpzAiEoqysuC66+DgQbjmGti2zXZFIhLCFEZEQpHLBc8+C506we7dZoTW/fttVyUiIUphRCRUxcTA3Llw1FGQkwODBpmRy0REfExhRCSUtWkDs2aZoeNfegmeeMJ2RSISghRGRELd+efDpElmPSsLPvzQbj0iEnIURkQEhg6Ffv2grMw0aN2yxXZFIhJCFEZExDRo/fe/oUsX+PFHuPJKDVEsIj6jMCIiRsOGpkFrUhKsWgX9+0N5ue2qRCQEKIyISKVWrUwgiYw0/73vPtsViUgIUBgRkarOPBOmTDHr48bB66/brUdEgl6dwsjkyZNp06YN0dHRpKWl8cUXX9TquBkzZuByuejVq1ddTisivjJwoGnUCqZha26u1XJEJLh5HEZmzpxJVlYWY8eOZcWKFXTs2JHMzEx27dp12OO2bt3KHXfcwTnnnFPnYkXEhyZNggsvhOJiuOIKyM+3XZGIBCmPw8jDDz/MrbfeysCBAzn55JOZMmUKMTExPP/88zUeU1ZWxp///GfGjRtH27Zt/1DBIuIjERFmht/27WH7dujdGw4csF2ViAQhj8JIaWkpOTk5ZGRkVL5AWBgZGRksX768xuPuu+8+EhISuPnmm2t1npKSEgoLC6ssImJBs2bw1lvQtCksXw5//auGjBeReudRGNmzZw9lZWUkJiZW2Z6YmEheXl61xyxbtoznnnuOqVOn1vo82dnZxMXFVSwpKSmelCki9al9e5g9u3LI+EOjtYqI1BOv9qbZt28f/fr1Y+rUqcTHx9f6uFGjRlFQUFCxbN++3YtVisgRZWTAY4+Z9REj4M037dYjIkElwpOd4+PjCQ8PJ/83Ddny8/NJSkr63f6bNm1i69atXH755RXbyn8ZRCkiIoL169dz3HHH/e64qKgooqKiPClNRLztb3+D1atNt9/rr4dly6BTJ9tViUgQ8OjOSGRkJKmpqSxatKhiW3l5OYsWLSI9Pf13+5944omsWrWK3NzciuWKK66gR48e5Obm6vGLSCBxueDxx00PG7cbLrsMduywXZWIBAGP7owAZGVlMWDAAM444wy6du3Ko48+itvtZuDAgQD079+fVq1akZ2dTXR0NB06dKhyfNOmTQF+t11EAkCDBjBrFqSnw7p1Zg6bpUshJsZ2ZSISwDwOI3369GH37t2MGTOGvLw8OnXqxIIFCyoatW7bto2wMA3sKhK0mjaF+fMhLQ2+/NIMijZ7NhrQWUTqyuU4/t9Pr7CwkLi4OAoKCoiNjbVdjogAfPwxXHAB/PwzjByJe3Q2jRubbxUVQaNGdssTEftq+/tbf8qISN2ccw48+6xZf+ABePFFq+WISOBSGBGRuuvfH0aPNuv/9392axGRgOVxm5Fg4Thmyg0R+YNG3Acbvsc9662KTW53/Z/m1/9m4+NBTdNEgkfIhpHiYiqeb4vIH+ECqs5N9ZtBmutdfj4kJHj3HCLiO/rbQkRERKwK2TsjMTGmxb+I1A/HgeLVm+HSnsT8uB3XxZfAjBlm9t964HZX3nHRsCYiwSVkw4jLpa6HIvWt8ZltYf4LcP75sOA1uGuwGT7e5arX89Tzy4mIZXpMIyL1Kz0dpk83ieGZZyA723ZFIuLnFEZEpP717m3msQG45x54+WW79YiIX1MYERHvGDIE7rzTrN90E3zwgd16RMRvKYyIiPc88AD07QsHD8JVV8HKlbYrEhE/pDAiIt4TFmaGie/eHfbtg0sugc2bbVclIn5GYUREvCsqCubNg44dzWhlmZmwa5ftqkTEjyiMiIj3xcXBu+9CmzbwzTfQs6e5UyIigsKIiPhKcjK8956ZWCYnB/70JygttV2ViPgBhRER8Z327eGdd8yIgwsXwoABUFZmuyoRsUxhRER8q0sXeO01aNDADBd/++1mLHkRCVkKIyLie5mZZiA0lwuefhrGjrVdkYhYpDAiInb06QNPPWXWx4+Hxx6zW4+IWKMwIiL2DBoE//ynWR82DF56yWo5ImKHwoiI2HX33TB8uFm/6SaYO9duPSLicwojImKXywWTJsGNN5qeNX36mC7AIhIyFEZExL6wMJg6Fa65Bn7+2cz6+9FHtqsSER9RGBER/xARAa+8ApdeCvv3w2WXwX//a7sqEfEBhRER8R+RkTB7NvToYYaLz8yEr76yXZWIeJnCiIj4l4YN4Y034Mwz4aefICMD1q61XZWIeJHCiIj4nyZNzLDxnTvD7t1w/vmwcaPtqkTESxRGRMQ/NWtm5q859VTIyzMz/YpIUIqwXYAtjgPFxbarEJHDij4K3vgAevbEve7bis1ut/dO+evPhvh409FHRLwrZMNIcTE0bmy7ChE5sgTgyypbEhN9c+b8fEhI8M25REKZMr+IiIhYFbJ3RmJioKjIdhUiUluOA8Ubv4fevYn5dg2uY9rAu+/C0UfX63nc7so7LzEx9frSIlKDkA0jLhc0amS7ChHxROPOrWDZ69C9O2z6GnqeB4sXQ5s2Xjmfy+WVlxWR39BjGhEJLK1bw5Il0K4dbN1qgsnWrXZrEpE/RGFERALPoUBy/PHw7bdw3nnwzTe2qxKROlIYEZHA1KqVCSQnnADbtsG558K6dbarEpE6UBgRkcDVsiUsXQodOsDOneYOyapVtqsSEQ8pjIhIYEtMNI1YO3eGXbtMG5KcHNtViYgHFEZEJPDFx8OHH0JaGvz4I1xwAXzyie2qRKSWFEZEJDg0bQrvvw/nnAMFBXDRReZrEfF7CiMiEjxiY2HBArj4YjPnw2WXwWuv2a5KRI5AYUREgktMDLzxBlxzDfz8M1x7Lbzwgu2qROQwFEZEJPhERsKrr8LNN0N5Odx0EzzyiO2qRKQGCiMiEpzCw2HqVMjKMl9nZcGoUWaSGxHxKwojIhK8XC6YNAnuv998/cAD5m7JwYN26xKRKhRGRCS4uVzmjshzz0FYmGk/0ru3aeAqIn5BYUREQsNNN8HcuRAdDfPnw4UXmjFJRMQ6hRERCR1XXAELF5oxST79FLp1gy1bbFclEvIURkQktJx9NixbBikpsH49pKfDl1/arkokpCmMiEjoOeUUWL4cOnaE/Hwzwd78+barEglZCiMiEppatYKPPjLDxhcXw5VXwpQptqsSCUkKIyISumJjzR2RgQPN4Gi33QYjR9quSiTkRNguwBbHUc8+EQFoAI8/B61PgPH34X7yOeABANxu31Tw68+j+HjTA1kklIRsGCkuhsaNbVchIv7BBYz4ZamUmOj7SvLzISHB9+cVsUn5W0RERKwK2TsjMTFQVGS7ChHxN44Dxeu3ww1/JmZdDq7IKJg8Ga67zmvndLsr78LExHjtNCJ+q05hZPLkyUycOJG8vDw6duzIE088QdeuXavdd+rUqbz00kusXr0agNTUVO6///4a9/cVlwsaNbJagoj4qcapKfD5fLjhBnjrLbj1eli/wsxtEx7u1XO7XF59eRG/5PFjmpkzZ5KVlcXYsWNZsWIFHTt2JDMzk127dlW7/5IlS7juuutYvHgxy5cvJyUlhYsuuojvv//+DxcvIuI1sbEwbx7cc4/5etIkuOwy2LvXZlUiQcnlOJ7Np52WlkaXLl148sknASgvLyclJYXbb7+dkbXoEldWVkazZs148skn6d+/f63OWVhYSFxcHAUFBcTGxnpSrojIHzdzpun+u38/HH+8mePmlFPq7eXd7soG9UVFumsrwaO2v789ujNSWlpKTk4OGRkZlS8QFkZGRgbLly+v1WsUFxfz888/07x58xr3KSkpobCwsMoiImJNnz7wySdmCPmNGyEtDWbPtl2VSNDwKIzs2bOHsrIyEn/T3y0xMZG8vLxavcaIESNo2bJllUDzW9nZ2cTFxVUsKSkpnpQpIlL/OneGnBw4/3xzK+Paa+HOO+HgQduViQQ8n3btfeCBB5gxYwZz584lOjq6xv1GjRpFQUFBxbJ9+3YfVikiUoMWLeC990wIAdOO5KKLoIY2cyJSOx6Fkfj4eMLDw8nPz6+yPT8/n6SkpMMeO2nSJB544AHef/99TjvttMPuGxUVRWxsbJVFRMQvRETAgw/CrFmmccfixeauyccf265MJGB5FEYiIyNJTU1l0aJFFdvKy8tZtGgR6enpNR734IMPMn78eBYsWMAZZ5xR92pFRPzFNdfAF1/ASSfBjh3Qo4cJKeXltisTCTgeP6bJyspi6tSpTJs2jbVr13LbbbfhdrsZOHAgAP3792fUqFEV+0+YMIF7772X559/njZt2pCXl0deXh5FGnFMRALdySebQPLnP0NZGYwYYWb//fFH25WJBBSPw0ifPn2YNGkSY8aMoVOnTuTm5rJgwYKKRq3btm1j586dFfs//fTTlJaWcvXVV5OcnFyxTJo0qf5+ChERWxo3hpdfhmeegagoMwtw586m942I1IrH44zYoHFGRCQg5OaaxzfffGNGav3HP2DUqCOO2qpxRiRYeWWcEREROYxOnUz330OPbe69FzIyQCNOixyWwoiISH2KjYVXXoFp08wtjiVLoGNHePNN25WJ+C2FERERb+jfH1asMO1HfvjBNGy99VZNFy5SDYURERFvad8eli+HO+4w0/E++6x5lFPL6TNEQoXCiIiIN0VFwcSJ8OGHZm6bTZvg7LNhzBgoLbVdnYhfUBgREfGF7t3hq69M49bychg/3ky497//2a5MxDqFERERX2na1DRunTkTmjc3XYG7dIEHHrBdmYhVITvOiONAcXG9vJSIiOfy82HoUJj/Fm5iSGR3xWab44z8+rMxPh7C9Cer/AG1/f0d4cOa/EpxceUgQyIivpcIzPj91kTfV1KT/HxISLBdhYQCZV4RERGxKmTvjMTEqLu/iPgHx4FitwNvvUXM3cNw7c4337j1L2ZI+bg4n9XidlfenYmJ8dlpJcSFbJsRERG/9NNPcOed8Nxz5uukJHjkEejTx4xV4mWaJ0fqk+amEREJRM2amcHRFi0yg6bl5cF118FFF8HGjbarE/EKhREREX90/vlmXJL77jMDp33wAXToYAZLU1dACTIKIyIi/ioqysz8u3q1uTNSWmoGSzvpJJgzxzQ2EQkCCiMiIv6uXTtYsABmz4ajj4Zt2+CaayAjwwQVkQCnMCIiEghcLrj6ali71jyqiYoy89106gSDB8Pu3bYrFKkzhRERkUASEwPjxplQ0rs3lJXBU0+ZuycPPggHDtiuUMRjCiMiIoHo2GPh9ddh8WLo3BkKC2HECNOeZMYMMxmfSIBQGBERCWTdu8OXX8K0adCqFWzdaroCd+kC77+vRq4SEBRGREQCXVgY9O8PGzaYrsBNmsCKFZCZCRdcAF98YbtCkcNSGBERCRYxMaYr8KZNMGwYREaaxzhpaaZ9yVdf2a5QpFoKIyIiwaZFCzOE/IYNMGCA6Ykzbx507AjXXgtr1tiuUKQKhRERkWB1zDHw4osmfPTpY7bNng2nngrXX69QIn5DYUREJNgd6mHz1Vdw1VWmUeurr5rh5f/0J9O+RMQihRERkVBx6qnw2muwcqUZQM3lMt2DU1OhZ09Ytky9b8QKl+P4/zuvtlMQe8JxNNeUiIS4tWvhoYdg5kxwzLgk7tPPJXHFuwDk50OjRjYLrJtff77Hx5vORmJHbX9/R/iwJr9SXAyNG9uuQkTEppOAZ39ZfvGrJzaJib6up/7l50NCgu0q5EiUF0VERMSqkL0zEhMDRUW2qxAR8S+OA8W79sF//kPM1Mdxbf/WfCM8wjR+ve02M7qry2W30MNwuyvv6sTE2K1Faidk24yIiMgRlJXBm2/Co4/CRx9Vbu/cGf72NzPsvB82KnG7Kx/DFxX5ZYkho7a/v/WYRkREqhcebkZuXboUcnLMAGpRUaY3zq23mrlwhg6F1attVyoBTmFERESO7PTTzQBq338PEyfCccdBQQE8/rjpMpyeDs89p+ffUicKIyIiUntHHQV33GGGmn/3XXPnJCICPvsMbrkFkpPNfz/+WGOWSK0pjIiIiOfCwuDii82gadu3w4QJcPzx5s7Ic8/BuedCu3Ywbhxs2WK7WvFzasAqIiL1w3FMQ9dp08wcOL9+ZHPWWabB6zXXeH3gDzVg9R9qwCoiIr7lcsF558Hzz0NeHrz8MmRkmO2ffAJDhkDLlnDJJSaw7N1ru2LxE7ozIiIi3vX99zBrFkyfDl9+Wbm9QQO48EIzT86VV0Lz5vVyOt0Z8R+1/f2tMCIiIr6zcaOZMXjWLFizpnJ7RAT06GFCyZVXQuvWdT6Fwoj/UBgRERH/tnYtzJlj2pesWlX1e2ecYULJpZdCp04ejfiqMOI/FEZERCRwbNwIb7wB8+bBp59W7RbcsiX07GmCyQUXQJMmh30phRH/oTAiIiKBKT8f3noL5s+HDz4w6eKQiAjTM+eiiyAz0wxNH1a1L4bCiP9QGBERkcB34IAZjv7tt82yeXPV78fHm7YmF1wA558P7drhLnYpjPgJhZEjcBwoLq6XlxIREV/ZvNncLfngAxNS3L8Zfr5Va9zpF5A45ykACgscmsT67wzDwU5h5Ah+fRtPRESCU4Irn/NabeactBLO7p3AqX9qT0R0hO2yQobCyBEojIiIhJ5GFNGl6Ua6nVzAmT0a0rVvWxI7tLBdVtBSGDkCPaYREQlOhz7f9+89wOa31/Hp2z/x8f+asHzP8RQS97v9jwn/jq7J2+lyWgmp58dx+tVtaXrM7/cTzymMiIiI/Er5wXLWvr2Z5a99z6efhfP5tiTWlrTFqWZmlOMivuX0pB10OrmEzmc1otMVR5N0WgKuMLU/8YTCiIiIyBEUfldIzsxv+HxhIV+ujiYnvxVbD6ZUu28L125Obbqd044p5NROYXQ49yhOvuQYGifpmX9NFEZERETq4IeNP7LitS2s/GgfuV83IHdnEutL21BOeLX7HxP+Hac038FJRxdzUodwTkyL48QLWnFU+6N8XLn/URgRERGpJ8V7ivn67S2sWvojX+WWs+rbJqzZ24q88sQajznK9QPtG+2gfWIB7Y/9mXanRNGuS1OOO6cVcUeHRpsUhREREREv+2HDD3z9/nesWV7Iuq/LWbu9Eev2JrGt7PAT/R3l+oG2MXm0bV7Asa1KadsujGNObkSbzs04Oi2Z6GYNffQTeJfCiIiIiCXu3cV8s3g7G5b/wIZVJWzY2oBNu2L5xp1EfnnCEY9PDNvF0Q13c3RcISkJJaQcDa3bRtH6xMa0OqUpLTsn0qBRpA9+kj9GYURERMQP7duxj03LdrJlxU9sWbufzVvC2JLXkG8Lm7K1JBk3R24Q66KcFq4fSI76kZaNC0ludoDkFgdJahlG0tGRJB4bQ+LxsSSc2Jy4lFhrvYAURkRERAKMU+7w4+a9fPt5HtvXFLJtwwG2b3PYlh/J93sb8X1xM74/mEgpUbV+zUhKSAj/gRaRBbSIcdOiSQktmh0k/iiHo1qEEZ/cgPiUhnS+6th6b8tS29/fdRoTd/LkyUycOJG8vDw6duzIE088QdeuXWvcf/bs2dx7771s3bqV448/ngkTJtCzZ8+6nFpERCRoucJcHNWuGUe1a8bpNexTXuawZ/0edq7aw84N+9ix+QA7tpeRt8tF3g+R5O1rxM79Tdl1sBlFNKGUKL4ra8l3+1vCfuAHYOvvX3dp+P849/aO3vvhDsPjMDJz5kyysrKYMmUKaWlpPProo2RmZrJ+/XoSEn7/HOzTTz/luuuuIzs7m8suu4zp06fTq1cvVqxYQYcOHerlhxAREQkVYeEuEk6OJ+HkeI4UHYr3FLN7/Y/kf7OP3VuK2P1dCbt3HmTPHvhhbzh7CiPZ427IDyWNSGzXxCf1V8fjxzRpaWl06dKFJ598EoDy8nJSUlK4/fbbGTly5O/279OnD263m/nz51dsO/PMM+nUqRNTpkyp1Tk1HLyIiIh3xcSAq56blnjlMU1paSk5OTmMGjWqYltYWBgZGRksX7682mOWL19OVlZWlW2ZmZnMmzevxvOUlJRQUlJS8XVhYaEnZdZKcbEmyhMRETmkqAgaNbJz7t8PyH8Ye/bsoaysjMTEqoO8JCYmkpeXV+0xeXl5Hu0PkJ2dTVxcXMWSklL90LwiIiIS+OrUgNXbRo0aVeVuSmFhYb0HkpgYkwJFRETE/F60xaMwEh8fT3h4OPn5+VW25+fnk5SUVO0xSUlJHu0PEBUVRVRU7bst1YXLZe92lIiIiFTy6DFNZGQkqampLFq0qGJbeXk5ixYtIj09vdpj0tPTq+wPsHDhwhr3FxERkdDi8WOarKwsBgwYwBlnnEHXrl159NFHcbvdDBw4EID+/fvTqlUrsrOzARg6dCjnnXceDz30EJdeeikzZszgyy+/5Jlnnqnfn0REREQCksdhpE+fPuzevZsxY8aQl5dHp06dWLBgQUUj1W3bthEWVnnDpVu3bkyfPp3Ro0dz9913c/zxxzNv3jyNMSIiIiKAhoMXERERL6nt72+P2oyIiIiI1DeFEREREbFKYURERESsUhgRERERqxRGRERExCqFEREREbFKYURERESsUhgRERERqxRGRERExCqPh4O34dAgsYWFhZYrERERkdo69Hv7SIO9B0QY2bdvHwApKSmWKxERERFP7du3j7i4uBq/HxBz05SXl7Njxw6aNGmCy+Wqt9ctLCwkJSWF7du3a86bI9C18oyuV+3pWtWerlXt6VrVnjevleM47Nu3j5YtW1aZRPe3AuLOSFhYGK1bt/ba68fGxurNWku6Vp7R9ao9Xava07WqPV2r2vPWtTrcHZFD1IBVRERErFIYEREREatCOoxERUUxduxYoqKibJfi93StPKPrVXu6VrWna1V7ula15w/XKiAasIqIiEjwCuk7IyIiImKfwoiIiIhYpTAiIiIiVimMiIiIiFVBH0YmT55MmzZtiI6OJi0tjS+++OKw+8+ePZsTTzyR6OhoTj31VN555x0fVWqfJ9fqxRdfxOVyVVmio6N9WK09H330EZdffjktW7bE5XIxb968Ix6zZMkSTj/9dKKiomjXrh0vvvii1+v0B55eqyVLlvzufeVyucjLy/NNwRZlZ2fTpUsXmjRpQkJCAr169WL9+vVHPC4UP7Pqcq1C9TPr6aef5rTTTqsY0Cw9PZ133333sMfYeE8FdRiZOXMmWVlZjB07lhUrVtCxY0cyMzPZtWtXtft/+umnXHfdddx8882sXLmSXr160atXL1avXu3jyn3P02sFZrS+nTt3VizffvutDyu2x+1207FjRyZPnlyr/bds2cKll15Kjx49yM3NZdiwYdxyyy289957Xq7UPk+v1SHr16+v8t5KSEjwUoX+Y+nSpQwePJjPPvuMhQsX8vPPP3PRRRfhdrtrPCZUP7Pqcq0gND+zWrduzQMPPEBOTg5ffvkl559/PldeeSVr1qypdn9r7ykniHXt2tUZPHhwxddlZWVOy5Ytnezs7Gr3v/baa51LL720yra0tDTnr3/9q1fr9AeeXqsXXnjBiYuL81F1/gtw5s6de9h97rrrLueUU06psq1Pnz5OZmamFyvzP7W5VosXL3YA56effvJJTf5s165dDuAsXbq0xn1C+TPr12pzrfSZValZs2bOs88+W+33bL2ngvbOSGlpKTk5OWRkZFRsCwsLIyMjg+XLl1d7zPLly6vsD5CZmVnj/sGiLtcKoKioiGOOOYaUlJTDJu1QF6rvqz+iU6dOJCcnc+GFF/LJJ5/YLseKgoICAJo3b17jPnpvGbW5VqDPrLKyMmbMmIHb7SY9Pb3afWy9p4I2jOzZs4eysjISExOrbE9MTKzx+XNeXp5H+weLulyrE044geeff5433niDV155hfLycrp168Z3333ni5IDSk3vq8LCQvbv32+pKv+UnJzMlClTeO2113jttddISUmhe/furFixwnZpPlVeXs6wYcM466yz6NChQ437hepn1q/V9lqF8mfWqlWraNy4MVFRUQwaNIi5c+dy8sknV7uvrfdUQMzaK/4nPT29SrLu1q0bJ510Ev/+978ZP368xcokkJ1wwgmccMIJFV9369aNTZs28cgjj/Dyyy9brMy3Bg8ezOrVq1m2bJntUvxeba9VKH9mnXDCCeTm5lJQUMCcOXMYMGAAS5curTGQ2BC0d0bi4+MJDw8nPz+/yvb8/HySkpKqPSYpKcmj/YNFXa7VbzVo0IDOnTvzzTffeKPEgFbT+yo2NpaGDRtaqipwdO3aNaTeV0OGDGH+/PksXryY1q1bH3bfUP3MOsSTa/VbofSZFRkZSbt27UhNTSU7O5uOHTvy2GOPVbuvrfdU0IaRyMhIUlNTWbRoUcW28vJyFi1aVOOzsvT09Cr7AyxcuLDG/YNFXa7Vb5WVlbFq1SqSk5O9VWbACtX3VX3Jzc0NifeV4zgMGTKEuXPn8uGHH3Lsscce8ZhQfW/V5Vr9Vih/ZpWXl1NSUlLt96y9p7zaPNayGTNmOFFRUc6LL77ofP31185f/vIXp2nTpk5eXp7jOI7Tr18/Z+TIkRX7f/LJJ05ERIQzadIkZ+3atc7YsWOdBg0aOKtWrbL1I/iMp9dq3Lhxznvvveds2rTJycnJcfr27etER0c7a9assfUj+My+ffuclStXOitXrnQA5+GHH3ZWrlzpfPvtt47jOM7IkSOdfv36Vey/efNmJyYmxrnzzjudtWvXOpMnT3bCw8OdBQsW2PoRfMbTa/XII4848+bNczZu3OisWrXKGTp0qBMWFuZ88MEHtn4En7ntttucuLg4Z8mSJc7OnTsrluLi4op99Jll1OVahepn1siRI52lS5c6W7Zscb766itn5MiRjsvlct5//33HcfznPRXUYcRxHOeJJ55wjj76aCcyMtLp2rWr89lnn1V877zzznMGDBhQZf9Zs2Y57du3dyIjI51TTjnFefvtt31csT2eXKthw4ZV7JuYmOj07NnTWbFihYWqfe9Q99PfLoeuz4ABA5zzzjvvd8d06tTJiYyMdNq2beu88MILPq/bBk+v1YQJE5zjjjvOiY6Odpo3b+50797d+fDDD+0U72PVXSegyntFn1lGXa5VqH5m3XTTTc4xxxzjREZGOi1atHAuuOCCiiDiOP7znnI5juN4996LiIiISM2Cts2IiIiIBAaFEREREbFKYURERESsUhgRERERqxRGRERExCqFEREREbFKYURERESsUhgRERERqxRGRERExCqFEREREbFKYURERESsUhgRERERq/4fP12aBVRPoCIAAAAASUVORK5CYII=\n"},"metadata":{}}]},{"cell_type":"markdown","source":["### 2. Rejection sampling\n","\n","To sample from the (normalised) ziggurat, we first draw a layer $i \\in \\{0, \\ldots, n-1\\}$ uniformly. Then we draw a point $(x,y)$ uniformly in the layer with the following algorithm.\n","\n","* If $i \\geq 1$, the layer is a rectangle. We first draw $x$ uniformly in $[0,x_i]$. If $x \\leq x_{i+1}$ then whatever the choice of $y$ the point $(x,y)$ will be below the graph of $f$ so there is no actual need to generate $y$ and we may return $X=x$ directly. If $x \\in (x_{i+1}, x_i]$, then we draw $y$ uniformly in $[y_i,y_{i+1}]$. If $y \\leq f(x)$ then we return $X=x$ ortherwise we reject the draw and restart from the beginning.\n","\n","* If $i=0$, we first define $x_0 = A/y_1 > x_1$ and draw $x$ uniformly on $[0,x_0]$. If $x \\leq x_1$ then we return $X=x$. Otherwise we need to draw $X$ in the tail of the distribution, namely according to the distribution $$\\frac{f(x)1_{x>x_1}}{\\int f(x)1_{x>x_1} dx}.$$ To do so we use the following *fallback algorithm*: let $x'=-\\log(U_1)/x_1$, $y'=-\\log(U_2)$; if $2y'>x'^2$ return $X=x'+x_1$ otherwise reject.\n","\n","Implement this algorithm in the next function."],"metadata":{"id":"z8nH7VaHLsXS"}},{"cell_type":"code","source":["x1 = 3.5\n","[x,y] = layers(x1)\n","n = len(x)\n","\n","def Ziggurat():\n","  while True:\n","    i = math.floor(n*Unif())\n","    if (i>=1):\n","      X = x[i-1]*Unif()\n","      if (X<=x[i]):\n","        return X\n","      else:\n","        Y = y[i-1] + (y[i]-y[i-1])*Unif()\n","        if (Y < f(X)):\n","          return X\n","    else:\n","      A = x[0]*(y[1]-y[0])\n","      x0 = A/y[0]\n","      X = x0*Unif()\n","      if (X<=x[0]):\n","        return X\n","      else:\n","        X = -np.log(Unif())/x[0]\n","        Y = -np.log(Unif())\n","        if (2*Y>X**2):\n","          return X+x[0]"],"metadata":{"id":"A842qW3WEjq-","executionInfo":{"status":"ok","timestamp":1697486979932,"user_tz":-120,"elapsed":211,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":34,"outputs":[]},{"cell_type":"code","source":["Ziggurat()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"dZcvfzwwF0py","executionInfo":{"status":"ok","timestamp":1697486853462,"user_tz":-120,"elapsed":217,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"23e3d6ee-feaa-49bb-cba0-2cc533c41383"},"execution_count":31,"outputs":[{"output_type":"execute_result","data":{"text/plain":["0.1964378888490481"]},"metadata":{},"execution_count":31}]},{"cell_type":"markdown","source":["We draw the associated histogram."],"metadata":{"id":"4cPnyQeuGkDj"}},{"cell_type":"code","source":["# Sample\n","N = 5000\n","Sample = np.zeros(N)\n","for k in range(N):\n","  Sample[k] = RandomSign(Ziggurat())\n","\n","plt.hist(Sample, bins=50, density=True)\n","\n","# Theoretical density\n","x_max = max(abs(Sample))\n","xx = np.linspace(-x_max, x_max, 500)\n","den = np.exp(-xx**2/2)/(2*math.pi)**.5\n","\n","plt.plot(xx,den,color='red')\n","\n","plt.show()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":430},"id":"aFFN8pN6A7iU","executionInfo":{"status":"ok","timestamp":1697487040492,"user_tz":-120,"elapsed":3782,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"606202ba-b4a0-4eb2-a4c1-67cafa53dbdc"},"execution_count":41,"outputs":[{"output_type":"display_data","data":{"text/plain":["<Figure size 640x480 with 1 Axes>"],"image/png":"iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNfUlEQVR4nO3deVxVdeL/8dcFBEQFF5RFUcTdFMkNLbcmEstKm5qxphmLbz/7ZsvUl6aSqbSmZtAycypHG2ectKZ0aspmWigjtUVckVxSc0NUvIgaXEVlvb8/DvcquXEROHd5Px+P87j3Hs49vK8ZvD3ncz7HYrfb7YiIiIi4MT+zA4iIiIhcigqLiIiIuD0VFhEREXF7KiwiIiLi9lRYRERExO2psIiIiIjbU2ERERERt6fCIiIiIm4vwOwA9aGqqor8/HxatGiBxWIxO46IiIjUgt1u5/jx40RHR+Pnd/FjKF5RWPLz84mJiTE7hoiIiNTB/v376dChw0W38YrC0qJFC8D4wKGhoSanERERkdqw2WzExMQ4f49fjFcUFsdpoNDQUBUWERERD1Ob4RwadCsiIiJuT4VFRERE3J4Ki4iIiLg9FRYRERFxeyosIiIi4vZUWERERMTtqbCIiIiI21NhEREREbenwiIiIiJur06FZc6cOcTGxhIcHExiYiJr166t1fsWL16MxWJh/PjxNdbb7XamTp1KVFQUTZs2JSkpiZ07d9YlmoiIiHghlwvLkiVLSE1NZdq0aWRnZ9OvXz+Sk5M5fPjwRd+Xm5vL7373O4YPH37O11544QVeeeUV5s2bx5o1a2jWrBnJycmcPn3a1XgiIiLihVwuLLNmzWLSpEmkpKTQu3dv5s2bR0hICAsWLLjgeyorK7nzzjt59tlniYuLq/E1u93O7Nmzeeqppxg3bhzx8fEsWrSI/Px8li5d6vIHEhEREe/jUmEpKytjw4YNJCUlndmBnx9JSUlkZWVd8H1/+MMfaNeuHffcc885X9u7dy9Wq7XGPsPCwkhMTLzgPktLS7HZbDUWERER8V4u3a35yJEjVFZWEhERUWN9REQE27dvP+97vvnmG/7+97+Tk5Nz3q9brVbnPn66T8fXfio9PZ1nn33Wlegi4ml274Z33oE1a+DHHyE8HIYMgV/9Cjp2NDudiDSyBr1K6Pjx4/zmN79h/vz5hIeH19t+09LSKC4udi779++vt32LiMmOH4f//V/o3h2efho++gi+/RY+/BDS0iAuDh54AHRkVcSnuHSEJTw8HH9/fwoKCmqsLygoIDIy8pztd+/eTW5uLjfddJNzXVVVlfGNAwLYsWOH830FBQVERUXV2GdCQsJ5cwQFBREUFORKdBHxBDk58POfw969xuvRo+GmmyAqCg4cgKVLYcUK+MtfYPlyo8z8ZFyciHgnl46wBAYGMmDAADIzM53rqqqqyMzMZOjQoeds37NnTzZv3kxOTo5zufnmm7nmmmvIyckhJiaGzp07ExkZWWOfNpuNNWvWnHefIuKdbrz7zxQPGQZ793IgtC233/EnYq/8LbEHOhO7LpjYQ12NkvLll9C+PWzbBoMHGyVHRLyeS0dYAFJTU7nrrrsYOHAggwcPZvbs2ZSUlJCSkgLAxIkTad++Penp6QQHB9OnT58a72/ZsiVAjfWPPPIIzz//PN26daNz5848/fTTREdHnzNfi4h4qS1b+OfiJwkrLWF9+178z23TsAU3P2ez2CkfA9Bu3J+Y//7z9LPupGD4tdzym5nkh7YDIHf62EaNLiKNw+XCMmHCBAoLC5k6dSpWq5WEhAQyMjKcg2bz8vLw83NtaMzjjz9OSUkJ9957L0VFRQwbNoyMjAyCg4NdjScinqa4GH7+c8JKS1jXvjd3/+IZSoJCLvqWwy3a8Ovbn+fdtx6n55F9vPGvZ7jlNzMv+T4R8VwWu91uNzvE5bLZbISFhVFcXExoaKjZcUSktux2uPVW+OADDoS25aa7ZvNjSFit3x5lK+SDNx8l8sQx3okfTdr1v9URFhEP4srvb91LSETMs2gRfPABBAZy//g0l8oKwKHQtjx802NUYeGOTZ+TtHNNAwUVEbOpsIiIOQoLITXVeP7ss2yK6l6n3azp2Je/Dr4FgOkZr0BRUT0FFBF3osIiIuZITYVjxyA+Hh599LJ2NWv4b9jZJobwk8WgSSVFvJIKi4g0vrVr4a23wGKB+fOhSZPL2l1ZQBP+cO0k48VrrxmXPIuIV1FhEZHGZbfDlCnG84kTjblU6sHXnfuzrGsiVFTA735XL/sUEffh8mXNIiJ14ZhDZcSeDSxavpxS/wB+1nQkB6vX14fnf3YP1+1dD598AqtXG/ceEhGvoCMsItJ47HYe/fotABb1v5GDYe3qdff7WkXDXXcZLzSWRcSrqLCISKO5at939LPu5FRAEHOH/KJhvsmTT4K/P2RkGEdZRMQrqLCISKOZvPo9ABb3G80xF+dcqbW4OGNsDEB6esN8DxFpdCosItIo+h7ayfB9OVRY/PjboFsa9ps98YTx+N//wq5dDfu9RKRRqLCISKOYtO4DAP7Te2S9j105R48ecMMNxhVJf/5zw34vEWkUKiwi0vCsVsbsWAXA3weOa5zv+X//Zzz+4x+a/VbEC6iwiEjD+9vfCKyqIDu6B1sjuzbO97z2WujbF0pKYMGCxvmeItJgNA+LiDSsigp4/XXAuJS5ocWeNa/LnVHD+ePmzez648skFXQ3ZtYF3dFZxAPpCIuINKyPPoIDBzjaNJRPe1zdqN/6w96jONkkiK7HDjDw4PeN+r1FpH6psIhIw1q4EID3+iZRGhDYqN/6RFAI/+05AoA7vvusUb+3iNQvFRYRaThHjsDHximaf/f5mSkRFvdLBmDs9m8IPX3ClAwicvlUWESk4SxZAuXl0L8/P7SNNSXCxugebA/vRHBFGTds/8aUDCJy+VRYRKThVJ8Ocs48awaLhQ/6XAPALd+vMC+HiFwWFRYRaRjbtsG6dcZ9fe64w9QoH/YaRRUWEvdvoX3xYVOziEjdqLCISMN4803j8frroV0Dz2x7CdbQcLI69QVgnI6yiHgkFRYRqX9VVfDWW8ZzM08HnWVp7+rTQluXG1P2i4hHUWERkfq3ahXs3w+hoXDTTWanAeDTHldT6t+Ebkf3w5YtZscRERepsIhI/fv3v43HceMgONjcLNVOBIXwVef+xgtHPhHxGCosIlK/7PYzheDWW83N8hOf9rjKePLee+YGERGXqbCISP1at844HdS8OYwebXaaGr7omkiZXwBs3Qrbt5sdR0RcoMIiIvXLcfTixhuhaVNzs/yELbg5qzr1M17otJCIR1FhEZH648angxycp4VUWEQ8igqLiNSfnBzYs8c4snL99WanOa/Puw0xJrPbuNHIKiIeQYVFROqP46jF9ddDs2bmZrmAH0PCYORI44WOsoh4jACzA4iI54udYtyR+fO/vUl34Lflcfynep1buu02+PJLeP99eOwxs9OISC3oCIuI1IuYIivdj+ZRYfFjRdxAs+NcnGMyuzVroLDQ3CwiUisqLCJSL362ex0A6zv0xhbc3OQ0l9ChA/TrZwwS/vRTs9OISC3UqbDMmTOH2NhYgoODSUxMZO3atRfc9v3332fgwIG0bNmSZs2akZCQwJuOm6JVu/vuu7FYLDWWMWPG1CWaiJjEUVgyuww2OUkt3Xij8fixG5+6EhEnlwvLkiVLSE1NZdq0aWRnZ9OvXz+Sk5M5fPj8t2xv3bo1Tz75JFlZWWzatImUlBRSUlL47LPPamw3ZswYDh065Fzeeeedun0iEWl0IWWnGJK3CYAvuwwyOU0tjR1rPH72GZSXm5tFRC7J5cIya9YsJk2aREpKCr1792bevHmEhISwYMGC824/atQobrnlFnr16kWXLl14+OGHiY+P55tvvqmxXVBQEJGRkc6lVatWdftEItLort73HUGVFexrGcnuNh3MjlM7gwdDeDgUF8O335qdRkQuwaXCUlZWxoYNG0hKSjqzAz8/kpKSyMrKuuT77XY7mZmZ7NixgxEjRtT42ooVK2jXrh09evRg8uTJHD169IL7KS0txWaz1VhExDzXVJ8O+rLLILBYTE5zabFTPib2yQz+HdEXgNefeNVYV72IiPtxqbAcOXKEyspKIiIiaqyPiIjAarVe8H3FxcU0b96cwMBAxo4dy6uvvsp1113n/PqYMWNYtGgRmZmZzJgxg5UrV3L99ddTWVl53v2lp6cTFhbmXGJiYlz5GCJSn+x25/gVjzkdVO3L6vE21+668Dg8EXEPjTIPS4sWLcjJyeHEiRNkZmaSmppKXFwco0aNAuD22293btu3b1/i4+Pp0qULK1as4Nprrz1nf2lpaaSmpjpf22w2lRYRs+TkEHniGCVNglkT09fsNC75uvOVVFj86HrsAB1/PEReqyizI4nIBbh0hCU8PBx/f38KCgpqrC8oKCAyMvLC38TPj65du5KQkMCjjz7KbbfdRnp6+gW3j4uLIzw8nF27dp3360FBQYSGhtZYRMQk1VfZfBObQFlAE5PDuMYW3Jx1MVcAcM2e9SanEZGLcamwBAYGMmDAADIzM53rqqqqyMzMZOjQobXeT1VVFaWlpRf8+oEDBzh69ChRUfrXjojb+/xzAFbGDTA5SN2s7GzkHr432+QkInIxLl8llJqayvz581m4cCHbtm1j8uTJlJSUkJKSAsDEiRNJS0tzbp+ens6yZcvYs2cP27Zt46WXXuLNN9/k17/+NQAnTpzgscceY/Xq1eTm5pKZmcm4cePo2rUrycnJ9fQxRaRBHD8O1QPuv4q90uQwdfN1ZyP30LzNNKnU5c0i7srlMSwTJkygsLCQqVOnYrVaSUhIICMjwzkQNy8vDz+/Mz2opKSE+++/nwMHDtC0aVN69uzJW2+9xYQJEwDw9/dn06ZNLFy4kKKiIqKjoxk9ejTPPfccQUFB9fQxRaRBrFgBFRXktoziQMsLnxZ2Z9+368yRkDDCTxZzZf4O1sb0MTuSiJyHxW63280OcblsNhthYWEUFxdrPItIY3roIXjtNd688gaeHn2/2Wnq7M//eZFx21by6tAJvDTiN+ROH2t2JBGf4Mrvb91LSETqbtkywBhw68kcp4WG5240OYmIXIgKi4jUTV4e7NgB/v5kdYw3O81l+bq6cMUf2knYqePmhhGR81JhEZG6qT66wuDB7n935ksoaBHOD2064oedq/Z9Z3YcETkPFRYRqRtHYRk92twc9cRxWmt4bo6pOUTk/FRYRMR1VVXwxRfG87Nus+HJvjp7HIvnX4sg4nVUWETEdRs3wtGj0KKFcddjL7Ampi9lfgHEFBfA7t1mxxGRn1BhERHXLV9uPI4cCU08azr+CzkVGMzG9j2NF47PJyJuQ4VFRFy3YoXxeM01psaob6sdN290fD4RcRsqLCLimooK+Oor43n1Hde9xeqOZxUWjWMRcSsqLCLimo0bjXsItWwJ/fqZnaZeZbfvSal/E8jPh507zY4jImdRYRER1zhOl4wYAf7+pkapb6UBgWyM7mG80GkhEbeiwiIirnEMSPWy00EONU4LiYjbUGERkdqrqICvvzaee9mAWwfnbQaWL9c4FhE3EmB2ABFxb7FTPnY+T8jfwdITJygKbs6V7+zHvvigickaRk50DwgKAqsVfvgBevQwO5KIoCMsIuKCIXmbAVgT0we7xTt/fJQGBMJVVxkvNB+LiNvwzp84ItIghuZtAs4a5+GtHONzNI5FxG2osIhIrQRUVjDwwPeAjxUWjWMRcQsqLCJSK/HWnTQrP82PwS3Y3jbW7DgNKzERgoOhoAC2bzc7jYigwiIiteQcv9LRe8evOAUFnRnHotNCIm7By3/qiEh9cRQW5/12vJ3jtJAG3oq4BRUWEbkk/6pKBhzcBvjA+BUHR2H5+muNYxFxAyosInJJvQ7vpVn5aWxBzdjRtpPZcRrHoEEQGGjMx7J7t9lpRHyeCouIXNLg/VsBWN++l/ePX3EIDjZKC8A335ibRURUWETk0gYeqC4sHXqbnKSRDRtmPKqwiJhOhUVELs5uZ9BBY/6Vdb5aWBz3TxIR06iwiMhFdSo6RNuSIkr9A9gU1d3sOI3r6quNxx9+gMOHzc0i4uNUWETkohzjVzZFdjfus+NLWrWCPn2M599+a24WER+nuzWLyEU5puP3pfErZ9+h+vmAGH7NFuanL+KPa84UttzpY82IJuKzdIRFRC5qYPX4lbUxV5icxByOzz2oeuCxiJhDhUVELuzwYbocOwjAhva9TA5jDseRpT7W3TQtO21yGhHfpcIiIhdWPW5je3gnbMHNTQ5jjvzQdhxs0ZYAexUJh3aYHUfEZ6mwiMiFVV/O60vjV87H8fkdA5BFpPGpsIjIhVVPmOar41cc1lV/fscAZBFpfHUqLHPmzCE2Npbg4GASExNZu3btBbd9//33GThwIC1btqRZs2YkJCTw5ptv1tjGbrczdepUoqKiaNq0KUlJSezcubMu0USkvpSUQHY2oCMsjgnz+udvx7+q0uQ0Ir7J5cKyZMkSUlNTmTZtGtnZ2fTr14/k5GQOX2BSpdatW/Pkk0+SlZXFpk2bSElJISUlhc8++8y5zQsvvMArr7zCvHnzWLNmDc2aNSM5OZnTpzXATcQ0a9ZAZSUHW7QlP7Sd2WlM9UN4R4qDmtGs/DS9C/aYHUfEJ7lcWGbNmsWkSZNISUmhd+/ezJs3j5CQEBYsWHDe7UeNGsUtt9xCr1696NKlCw8//DDx8fF8U32o2W63M3v2bJ566inGjRtHfHw8ixYtIj8/n6VLl17WhxORy1A94NbXj64A2C1+zqukBum0kIgpXCosZWVlbNiwgaSkpDM78PMjKSmJrKysS77fbreTmZnJjh07GDFiBAB79+7FarXW2GdYWBiJiYm12qeINJDVqwHIbt/T5CDuwVFY+udvNzmJiG9yaabbI0eOUFlZSURERI31ERERbN9+4f+Ji4uLad++PaWlpfj7+/OXv/yF6667DgCr1ercx0/36fjaT5WWllJaWup8bbPZXPkYInIpdrtxSgjI8bX7B12Ao7hdeVCFRcQMjTI1f4sWLcjJyeHEiRNkZmaSmppKXFwco0aNqtP+0tPTefbZZ+s3pIicsXs3HD0KQUF8HxFndhq38F1UdyotfrQ/XkjE8SNmxxHxOS6dEgoPD8ff35+CgoIa6wsKCoiMjLzwN/Hzo2vXriQkJPDoo49y2223kZ6eDuB8nyv7TEtLo7i42Lns37/flY8hIpdSfTqI/v0p929ibhY3cTKwKdvbxgLQX0dZRBqdS4UlMDCQAQMGkJmZ6VxXVVVFZmYmQ4cOrfV+qqqqnKd0OnfuTGRkZI192mw21qxZc8F9BgUFERoaWmMRkXpUfTqIxERzc7gZx2khjWMRaXwunxJKTU3lrrvuYuDAgQwePJjZs2dTUlJCSkoKABMnTqR9+/bOIyjp6ekMHDiQLl26UFpayieffMKbb77J3LlzAbBYLDzyyCM8//zzdOvWjc6dO/P0008THR3N+PHj6++TikjtOY6wDBkCG82N4k6yo3vym42f6AiLiAlcLiwTJkygsLCQqVOnYrVaSUhIICMjwzloNi8vDz+/MwduSkpKuP/++zlw4ABNmzalZ8+evPXWW0yYMMG5zeOPP05JSQn33nsvRUVFDBs2jIyMDIKDg+vhI4qIS06dgpwc43liImzUdPQOjiMsfQp2QWkpBAWZnEjEd1jsdrvd7BCXy2azERYWRnFxsU4PiVyuVavg6qshIgIOHSI27ROzE7kPu50Nr95Jm1M248/JhVPhInIuV35/615CIlLT2eNXLBZzs7gbi4Xs6vlY0DxRIo1KhUVEajp7/IqcwzmRngqLSKNSYRGRmnSF0EVlR6uwiJhBhUVEzrBaYd8+41TQoEFmp3FLmyK7UWHxg4MHQXNAiTQaFRYROcNxdOWKK6BFC3OzuKlTgcFsa9fZeLFqlblhRHyICouInKHxK7WyQQNvRRqdCouInKHxK7WigbcijU+FRUQMlZWwbp3xXEdYLso58HbjRjh92twwIj6iUe7WLCLuKXbKx87nPQpz+ezECU4ENiV+0R6q/PaZmMy9HQiLMCbWKyiADRuMifZEpEHpCIuIAHBl9f1xvovqRpWfv8lp3JzFcmaWWw28FWkUKiwiAkDCoR8AyInqYXISD+EoLBrHItIoVFhEBIAr840jLBsd4zPk4q66ynjMygLPvyWbiNtTYRERmpeepNsRYxK0nOjuJqfxEP37g7+/MdmeJpATaXAqLCJC/KEf8MPO/rAIjjRrZXYczxASAvHxxnPH5eAi0mBUWESEK/N3ALAxWuNXXOKYr0aFRaTBqbCICAmHjMKiAbcuUmERaTQqLCK+zm7XEZa6ckywt2EDlJebm0XEy6mwiPi4DsUFhJ8spswvgO8j4syO41m6d4ewMDh1CrZsMTuNiFdTYRHxcf2rj658HxFHaUCgyWk8jJ8fDB5sPNdpIZEGpcIi4uMSdDro8mgci0ijUGER8XEav3KZVFhEGoUKi4gPC6wop/fh3YCuEKozR2HZvh2Ki83NIuLFVFhEfFjvw3sIqqzgaNNQ8lpGmh3HM7VtC507G9Pzr1tndhoRrxVgdgARMY9j/EpOdA/jDsRSa7FTPnY+fyU4hpvZy4t/WMScL0qd63OnjzUjmohX0hEWER+m8Sv1w3E6zTEBn4jUPxUWER+mGW7rh+OGkQn5P+jOzSINRIVFxFcVFtKpyEoVFr7THZovy9aILpT5BdD2ZBEdbIfNjiPilVRYRHxV9WW4u9t04HhQM5PDeLbSgEC2tesMnBkXJCL1S4VFxFetXg1o/Ep9cZwWulKFRaRBqLCI+KrqIyw5Kiz1wjnwVoVFpEGosIj4oqoqWLsW0BGW+uIofn0KdtOkUnduFqlvKiwivmj7drDZONkkiB/CO5mdxivsbRVNUXBzgirL6Xk41+w4Il5HhUXEF1WPX9kU2Y1KP3+Tw3gJi4Xvoqovb9Z8LCL1ToVFxBc5B9z2NDmId9E4FpGGU6fCMmfOHGJjYwkODiYxMZG11efCz2f+/PkMHz6cVq1a0apVK5KSks7Z/u6778ZisdRYxowZU5doIlIbzgG3mn+lPm10TCB36AeTk4h4H5cLy5IlS0hNTWXatGlkZ2fTr18/kpOTOXz4/JMlrVixgjvuuIPly5eTlZVFTEwMo0eP5uDBgzW2GzNmDIcOHXIu77zzTt0+kYhc3IkTsGULABs1w229cpwS6nLsIKGnT5icRsS7uFxYZs2axaRJk0hJSaF3797MmzePkJAQFixYcN7t//nPf3L//feTkJBAz549+dvf/kZVVRWZmZk1tgsKCiIyMtK5tGrVqm6fSEQubv164yqhmBgOt2hjdhqv8mNIGLktowCdFhKpby4VlrKyMjZs2EBSUtKZHfj5kZSURFZWVq32cfLkScrLy2ndunWN9StWrKBdu3b06NGDyZMnc/To0Qvuo7S0FJvNVmMRkVqqHr/CkCHm5vBSOTotJNIgXCosR44cobKykoiIiBrrIyIisFqttdrHE088QXR0dI3SM2bMGBYtWkRmZiYzZsxg5cqVXH/99VRWVp53H+np6YSFhTmXmJgYVz6GiG+rHr9CYqK5ObyUBt6KNIyAxvxm06dPZ/HixaxYsYLg4GDn+ttvv935vG/fvsTHx9OlSxdWrFjBtddee85+0tLSSE1Ndb622WwqLSK1YbfXPMLy3yJT43gjxwRyCYeq79xssZicSMQ7uHSEJTw8HH9/fwoKCmqsLygoIDIy8qLvnTlzJtOnT+fzzz8nPj7+otvGxcURHh7Orl27zvv1oKAgQkNDaywiUgv794PVCgEB0L+/2Wm80vft4ij1D6D1KRvs2WN2HBGv4VJhCQwMZMCAATUGzDoG0A4dOvSC73vhhRd47rnnyMjIYODAgZf8PgcOHODo0aNERUW5Ek9ELsVxdKVfP2ja1NwsXqosoAnb2sUZLxyn30Tksrl8lVBqairz589n4cKFbNu2jcmTJ1NSUkJKSgoAEydOJC0tzbn9jBkzePrpp1mwYAGxsbFYrVasVisnThiX/J04cYLHHnuM1atXk5ubS2ZmJuPGjaNr164kJyfX08cUEUDjVxqJ8/5MKiwi9cblMSwTJkygsLCQqVOnYrVaSUhIICMjwzkQNy8vDz+/Mz1o7ty5lJWVcdttt9XYz7Rp03jmmWfw9/dn06ZNLFy4kKKiIqKjoxk9ejTPPfccQUFBl/nxRKQGXSHUKHKq52NRYRGpPxa73W43O8TlstlshIWFUVxcrPEsIhdSVgZhYXD6NOzYAd27EzvlY7NTeaVOP+az8q/3QmAg2Gygf3yJnJcrv791LyERX7Fpk1FWWrWCbt3MTuPV9rWM4ljTUKMkfved2XFEvIIKi4ivOHv8ii61bVgWC99FVZdCnRYSqRcqLCK+QuNXGpVjAjkVFpH6ocIi4it0hVCjckwg5yyKInJZGnWmWxFpXI5BtS1P2cjZuROAfp8WUbxCg20bmvNKod274cgRCA83N5CIh9MRFhEfkJBv3Ihvd+v2FDdtYXIa31DctAV0ry4ta9eaG0bEC6iwiPiAK6tvxOc8TSGNw3H6TeNYRC6bCouID0g4ZBSWjVEqLI1KhUWk3qiwiHg5i72KhOojLBt1hKVxOQrL2rXGnZtFpM5UWES8XNyxg4SVlnAqIIgdbWPNjuNb4uONWW5//BGqBz2LSN2osIh4OceA282RXajw14WBjSowEPr3N57rtJDIZVFhEfFyV+ZvB2BjdE+Tk/gox0R9mo9F5LKosIh4uYRDxhEW57wg0rg08FakXqiwiHix4PLT9Dy8F9ARFtM4Cst338GpU+ZmEfFgKiwiXqyvdRcB9iqszVtjDdVMq6bo1AnatYOKCti40ew0Ih5LhUXEi13pvJxZR1dMY7HotJBIPVBhEfFiCc4ZbjV+xVQqLCKXTYVFxIvpCIubcFwppMIiUmcqLCLe6sABok4cpcLix+aIrman8W2DBhmnhnJzoaDA7DQiHkmFRcRbVf9rfkfbWE4FBpscxseFhkKvXsZzHWURqRNNeynirap/MWr8inlip3zsfP6CfzS/5Htem/5PZq7yd67PnT7WjGgiHkdHWES8VfXMqhq/4h5yqm886RgILSKuUWER8UYVFbB+PQAbo3SHZneQU/3fIf7QTiz2KpPTiHgeFRYRb7R5M5w6hS2oGXvatDc7jQA72nbiZJMgQstO0uXoAbPjiHgcFRYRb+QYvxLVHbtF/5u7g0o/fzZHdgPOXG4uIrWnn2Qi3qh6/IpueOheNlb/91BhEXGdCouIN6ouLNntNeDWnTgH3h5SYRFxlQqLiLc5dgx2OKbk14Bbd+IYeNujcB9Ny06bnEbEs6iwiHgbx8Rk3btT1DTU3CxSgzU0HGvz1vjbq+hbsMvsOCIeRYVFxNtUnw5y3r9G3IrmYxGpGxUWEW+TlWU8qrC4JcdpIRUWEdeosIh4k6qqM6eEhg41N4uc10bnwNsfTE4i4llUWES8yfbtYLNBSAj06WN2GjmPzZFdqbT4EX38CBHHj5gdR8Rj1KmwzJkzh9jYWIKDg0lMTGTt2rUX3Hb+/PkMHz6cVq1a0apVK5KSks7Z3m63M3XqVKKiomjatClJSUns3LmzLtFEfJvjdNCgQRCge5u6o5OBTfkhvCMACfk6yiJSWy4XliVLlpCamsq0adPIzs6mX79+JCcnc/jw4fNuv2LFCu644w6WL19OVlYWMTExjB49moMHDzq3eeGFF3jllVeYN28ea9asoVmzZiQnJ3P6tC77E3GJY8CtTge5NcdpoSs1H4tIrblcWGbNmsWkSZNISUmhd+/ezJs3j5CQEBYsWHDe7f/5z39y//33k5CQQM+ePfnb3/5GVVUVmZmZgHF0Zfbs2Tz11FOMGzeO+Ph4Fi1aRH5+PkuXLr2sDyfic3SFkEfQwFsR17lUWMrKytiwYQNJSUlnduDnR1JSElmOQ9GXcPLkScrLy2ndujUAe/fuxWq11thnWFgYiYmJF9xnaWkpNputxiLi84qLYetW47kKi1vLiTam6O9r3QWVlSanEfEMLhWWI0eOUFlZSURERI31ERERWK3WWu3jiSeeIDo62llQHO9zZZ/p6emEhYU5l5iYGFc+hoh3WrcO7Hbo3Bl+8v+TuJddbWI4HtiUZuWnz5RMEbmoRr1KaPr06SxevJgPPviA4ODgOu8nLS2N4uJi57J///56TCnioXQ6yGNU+fmzKcq4czO1PDot4utcKizh4eH4+/tTUFBQY31BQQGRkZEXfe/MmTOZPn06n3/+OfHx8c71jve5ss+goCBCQ0NrLCI+TxPGeZQN0b2MJyosIrXiUmEJDAxkwIABzgGzgHMA7dCLXJXwwgsv8Nxzz5GRkcHAgQNrfK1z585ERkbW2KfNZmPNmjUX3aeInMVu1xVCHia7fXVhWbXK3CAiHsLliRpSU1O56667GDhwIIMHD2b27NmUlJSQkpICwMSJE2nfvj3p6ekAzJgxg6lTp/L2228TGxvrHJfSvHlzmjdvjsVi4ZFHHuH555+nW7dudO7cmaeffpro6GjGjx9ff59UxJvt2mXcpTk4GPr1MzuN1EJ2+57Gk507obAQ2rY1N5CIm3O5sEyYMIHCwkKmTp2K1WolISGBjIwM56DZvLw8/PzOHLiZO3cuZWVl3HbbbTX2M23aNJ555hkAHn/8cUpKSrj33nspKipi2LBhZGRkXNY4FxGf4jitMGAABAaam0VqxRbcnJ1tYuh2dL9xdOymm8yOJOLW6jQV5oMPPsiDDz543q+tWLGixuvc3NxL7s9isfCHP/yBP/zhD3WJIyIacOuRNrTvZRSWVatUWEQuQXN3i3io2CkfO59/9P7n9AEm72rCp2etF/e2oX1Pbt/0ucaxiNSCbn4o4uGalp2m5+G9wFnjIsQjZDuuFFq3DsrLzQ0j4uZUWEQ8XLx1JwH2KvJbhFPQItzsOOKCPW3aQ6tWcOoUfPed2XFE3JoKi4iHu7L6fjSOG+qJ57Bb/M5chq7TQiIXpcIi4uH6528HIDtap4M8kgqLSK2osIh4MrudK6sLy0YVFs901VXGo2a8FbkoFRYRDxZTXEDbkiLK/ALYGhFndhypi8GDwc8P8vLgwAGz04i4LRUWEQ828MD3AGyJ7EJpkyCT00idNG8Ojvur6SiLyAWpsIh4sIEHjcKyvn1vk5PIZdFpIZFLUmER8WADDmwDYEOHXiYnkcviKCwaeCtyQSosIh4q9PQJuh/JA4wp3sWDOa4Uys6G06fNzSLiplRYRDxU/4Pb8cPO3lZRHGnWyuw4cjk6d4aICGO22w0bzE4j4pZUWEQ81ICD1aeDNH7F81ksOi0kcgkqLCIe6syAW50O8gqaQE7kolRYRDxReTkJ+T8AsL6DjrB4hbOvFLLbzc0i4oYCzA4gInWQk0PTilKKgpuzu00Hs9PIZYid8jEAQRVlbPYLILCggOGTF7C/ZSQAudPHmhlPxG3oCIuIJ/r2W8C4Oshu0f/G3qA0IJCtEV2AM+OTROQM/aQT8URnFRbxHhvaG/eDUmEROZcKi4insdudhUXjV7yLo4A6brkgImeosIh4mtxcOHSIMr8AvovsZnYaqUfrYq4AoEfhPkJPnzA5jYh7UWER8TTVR1e2RuiGh97mSLNW7GkVjR92nRYS+QkVFhFP4zwdpPEr3mhdB+Moy+D9W01OIuJeVFhEPI2jsGiGW6/kOC006IAKi8jZVFhEPElREWzZAugOzd5qbfURlvhDOwkqLzU5jYj7UGER8SSrVxtXCXXpohseeqm8lpEUNG9NYFUFCYd+MDuOiNtQYRHxJNWng7j6anNzSMOxWJzjWHRaSOQMFRYRT/LVV8bjsGHm5pAGtbZ6fh0NvBU5Q4VFxFOcPg1r1hjPR440N4s0KMfA2/7526GiwuQ0Iu5BhUXEU6xbB6WlEBEB3TRhnDfbEd4JW1Azmpedgu++MzuOiFtQYRHxFI7TQSNGgMVibhZpUFV+/qx33Cfq66/NDSPiJlRYRDzF2YVFvJ7jtJAKi4hBhUXEE5SXn7lCSIXFJzjmY+Hrr41L2UV8nAqLiCfYuBFKSqBVK+jTx+w00gg2R3aj1L8JFBbCD5qPRaROhWXOnDnExsYSHBxMYmIia9euveC2W7du5dZbbyU2NhaLxcLs2bPP2eaZZ57BYrHUWHr27FmXaCLeyXE6aPhw8NO/M3xBWUATcqJ7GC90WkjE9cKyZMkSUlNTmTZtGtnZ2fTr14/k5GQOHz583u1PnjxJXFwc06dPJzIy8oL7veKKKzh06JBz+eabb1yNJuK9NH7FJ9U4LSTi41wuLLNmzWLSpEmkpKTQu3dv5s2bR0hICAsWLDjv9oMGDeLFF1/k9ttvJygo6IL7DQgIIDIy0rmEh4e7Gk3EO1VVnfmFpcLiU9ZVTyCnwiLiYmEpKytjw4YNJCUlndmBnx9JSUlkZWVdVpCdO3cSHR1NXFwcd955J3l5eRfctrS0FJvNVmMR8VqbNxs3PWzWDK680uw00og2tO9lnALcuxcOHDA7joipXCosR44cobKykoiIiBrrIyIisFqtdQ6RmJjIG2+8QUZGBnPnzmXv3r0MHz6c48ePn3f79PR0wsLCnEtMTEydv7eI23OcDrr6aggIMDeLNKqSoBDo39944fh7IOKj3OKn3/XXX+98Hh8fT2JiIp06deJf//oX99xzzznbp6WlkZqa6nxts9lUWsSrxE752Pl8ztIljAVePBXBnLPWi48YNQrWr4fly+FXvzI7jYhpXDrCEh4ejr+/PwUFBTXWFxQUXHRAratatmxJ9+7d2bVr13m/HhQURGhoaI1FxCvZ7c4b4K11TCQmvuWaa4zH5cvNzSFiMpcKS2BgIAMGDCAzM9O5rqqqiszMTIYOHVpvoU6cOMHu3buJioqqt32KeKIuxw7Q9mQRpf5N+C6qh9lxxAzDhoG/P+zeDfv3m51GxDQuXyWUmprK/PnzWbhwIdu2bWPy5MmUlJSQkpICwMSJE0lLS3NuX1ZWRk5ODjk5OZSVlXHw4EFycnJqHD353e9+x8qVK8nNzWXVqlXccsst+Pv7c8cdd9TDRxTxXEPzNgPG4MuygCYmpxFThIbCgAHGcx1lER/m8hiWCRMmUFhYyNSpU7FarSQkJJCRkeEciJuXl4ffWRNb5efnc+VZVzbMnDmTmTNnMnLkSFasWAHAgQMHuOOOOzh69Cht27Zl2LBhrF69mrZt217mxxPxbEP3GXfqzerY1+QkYqprroG1a43CMnGi2WlETGGx2z3/JhU2m42wsDCKi4s1nkW8QuyUj7HYq9jw6q9pfcrGrXe+wAbHnBziU3Knj4XPPoMxY6BTJ8jNNTuSSL1x5fe35vgWcVM9CvfR+pSNkibBbIrqZnYcMZPjkvZ9+1RYxGepsIi4qav2bQJgXYcrKPfX+BWf1rw5DB5sPNc4FvFRKiwibmponjF+ZVWneJOTiFsYNcp4VGERH6XCIuKG/KsqSczbAsCqTv1MTiNu4ez5WDx/6KGIy1RYRNzQFQW7CS07SXFQM75v19nsOOIOrroKmjQx7im0e7fZaUQanQqLiBtyjF9Z07EvVX7+JqcRtxASAkOGGM91Wkh8kAqLiBu6qnr+lVUdNX5FzuI4LVQ9h5WIL3GLmx+KyFnKyhh04HtAA26l5o0wh+QFsxgo/PATBj3xEVgsQPVcLSJeTkdYRNzNmjU0rSjlSEgYP4R3MjuNuJHs6F6cbBJE25IiehzZZ3YckUalwiLibqrHJ2R1jHf+C1oEoCygCWs79AFg2N6NJqcRaVwqLCLupvpu6LqcWc7n69gEAIbn5piaQ6SxqbCIuJPjxyErC4BvVVjkPL7pbNxMNnH/FgIryk1OI9J4VFhE3MnKlVBeTm7LKPJaRZmdRtzQjvBOHG7WiqYVpQw4uM3sOCKNRoVFxJ189hkAX1f/K1rkHBYL31SfFhqWq3Es4jtUWETcyeefA/B1rAqLXNiZwpJjag6RxqTCIuIucnPhhx/A358szb8iF/FNpwQA+lp30fKUzdwwIo1EhUXEXVQfXWHIEI4HNTM3i7i1wy3asCO8I37YuTr3O7PjiDQKFRYRd+EoLMnJ5uYQj/BN9WlDjWMRX6HCIuIOKiqc868werS5WcQj1JiPxW43NYtIY1BhEXEH69dDURG0bAkDB5qdRjzAmpi+lPkF0MF2GHbtMjuOSINTYRFxB9WXM5OUBP7+5mYRj3AqMJjs9j2NF47TiSJeTIVFxB1o/IrUwcq4AcaTTz81N4hII1BhETFbURGsWWM8v+46U6OIZ1keV3368Msv4dQpc8OINDAVFhGzffklVFZCjx7QqZPZacSDbG8bS36LcKOsrFxpdhyRBqXCImK2Tz4xHnU6SFxlsbDCcZTl44/NzSLSwFRYRMxUVXXmF81NN5mbRTzS8i7VheWTT3R5s3g1FRYRM2Vng9UKzZvDiBFmpxEP9G2nfhAYCHv2GLd2EPFSAWYHEPFFsVOMoyoPf/M2/wd8Gh3P5KnLzA0lHulkYFMYORKWLTOO1vXoYXYkkQahIywiJvrZ7nUAfNllkMlJxKPdcIPx6BgPJeKFVFhETNL2xDH6WXcCZ41DEKkLR2H56is4ftzcLCINRIVFxCTX7F4PQE5UN440a2VyGvFo3btD165QXg5ffGF2GpEGocIiYpJrd68F4Msug01OIl5Bp4XEy6mwiJggsKKcYbk5AGRq/IrUh7MLiy5vFi9Up8IyZ84cYmNjCQ4OJjExkbVr115w261bt3LrrbcSGxuLxWJh9uzZl71PEU+XuH8zzcpPY23emq0RXcyOI95g5Eho1gzy843L5UW8jMuFZcmSJaSmpjJt2jSys7Pp168fycnJHD58+Lzbnzx5kri4OKZPn05kZGS97FPE0zmuDloeNxAsFpPTiFcIDoYxY4znH3xgbhaRBuByYZk1axaTJk0iJSWF3r17M2/ePEJCQliwYMF5tx80aBAvvvgit99+O0FBQfWyTxGPZrdz7a7q8StdNX5F6tH48cbj0qVmphBpEC4VlrKyMjZs2EBSUtKZHfj5kZSURFZWVp0C1GWfpaWl2Gy2GouIx9iyhY7FBZT6NzFmKRWpL2PHQkAAbN0KO3eanUakXrlUWI4cOUJlZSURERE11kdERGC1WusUoC77TE9PJywszLnExMTU6XuLmKL6cP1Xna80ZikVqS+tWsGoUcZzHWURL+ORU/OnpaWRmprqfG2z2VRaxHO8/z4An3W/yuQg4i0ct3oA+LWlG8/zBRtmL+DWo72d63OnjzUjmki9cekIS3h4OP7+/hQUFNRYX1BQcMEBtQ2xz6CgIEJDQ2ssIh5hzx747jsqLH58ofEr0gCWdUsEYED+dtqeOGZyGpH641JhCQwMZMCAAWRmZjrXVVVVkZmZydChQ+sUoCH2KeK2qk8HrenYh6KmKtpS/wpahJMT1R2A63atMTmNSP1x+Sqh1NRU5s+fz8KFC9m2bRuTJ0+mpKSElJQUACZOnEhaWppz+7KyMnJycsjJyaGsrIyDBw+Sk5PDrl27ar1PEa9RXVg+66YyLg3ns+7G36/kH+p2MYSIO3J5DMuECRMoLCxk6tSpWK1WEhISyMjIcA6azcvLw8/vTA/Kz8/nyiuvdL6eOXMmM2fOZOTIkaxYsaJW+xTxClYrrFoFwOcqLNKAPu82hCdWLmTovk20KC3heFAzsyOJXDaL3e75czjbbDbCwsIoLi7WeBZxX6+/DvfdB4MHE3vNVLPTiJf7Yv59dD12gN/e9Bj/6T1Sg27FLbny+1v3EhJpLP/+t/H485+bm0N8wufdhwCQ/MMqk5OI1A8VFpHGUFgIX35pPL/1VnOziE/4pMcwAH62ez0hZadMTiNy+VRYRBrDv/8NlZUwYAB07Wp2GvEBWyK6kNsyiqYVpc5bQYh4MhUWkcawZInxOGGCuTnEd1gsfNRrOAA3bf/a5DAil0+FRaShHToEK1caz3/5S3OziE/5qKdRWEbuWQ/FxSanEbk8KiwiDe2998Buh6FDoVMns9OID9neNpadbWIIqqyADz80O47IZVFhEWloOh0kZrFYnEdZnH8PRTyUCotIQ9q/H779FiwW+MUvzE4jPsgxjoXPP4djureQeC4VFpGG5PhX7fDhEB1tbhbxSbvbxLCtbSxUVJyZC0jEA6mwiDSkN980Hn/1K3NziE9besUo44nj76OIB3L5XkIicnGxUz4GoNfhPXy6aROl/gEM2tYSW/V6kcb2Ya9RpK1cCF9/DXv2QFyc2ZFEXKYjLCIN5OdbjJltv+iaiC24uclpxJdZQ8Ph2muNF2+9ZW4YkTpSYRFpAP5VlYz73ph75f0+PzM5jQgwcaLx+OabxmX2Ih5GhUWkAQzLzaFdyY8cbRrKys4DzI4jArfcAiEhsGsXrF5tdhoRl6mwiDQAx+mg//QeSYW/hoqJG2je/MyNNxctMjeLSB2osIjUs+alJ0nemQXA+1fodJC4EcdpocWL4fRpc7OIuEiFRaSejft+BcEVZexsE8PmSN2ZWdzINddAhw5QVAQffGB2GhGXqLCI1LPbv/sMgMX9ko0ZbkXchb8//M//GM/nzzc3i4iLVFhE6lN2Nn0LdlPqH8D7V1xjdhqRc/3P/xhFevlyYwCuiIfQaECR+lT9r9bPul/FjyFhJocROSP2rIkL34jtz6i9G5j76zRmjLrbuT53+lgTkonUjo6wiNSXkhL45z8BeKdfsslhRC7M8ffzti1fEFBZYXIakdpRYRGpL0uWwPHj5LaMYnXHvmanEbmgzK6DKWzWkrYlRVy7a63ZcURqRYVFpL789a8ALOk3GrtF/2uJ+6rwD+DdvkkA/Oq7DJPTiNSOfqqK1Id162DNGmjShPf6JJmdRuSS3uk3hiosjNybTedjB82OI3JJKiwi9eHVV43HCRMobN7K3CwitbC/ZSSZXQcBMDH7I5PTiFyaCovI5SooMMavAPz2t+ZmEXHBwv43AfCLzV/QvPSkyWlELk6FReRy/fWvUFYGQ4bAoEFmpxGptW9iE9jZJobmZae4bfMXZscRuSgVFpHLUVYGc+cazx96yNwsIq6yWFg4wDjKMjH7I6iqMjmQyIWpsIhcjvffh0OHIDISbrvN7DQiLnv/imuwBTUj7sd8yNAVQ+K+VFhE6spuhxdfNJ5PngyBgebmEamDk4FN+Vf1Jc7MnGluGJGLUGERqavMTMjOhpAQuP9+s9OI1NmCQeMo9/M37i+0VhPJiXtSYRGpq+nTjcf/9/8gPNzcLCKXIT+0HR/2HmW8mDHD1CwiF6KbH4q4wHEDub6HdvLfzEzK/fwZVZbAwbNuLCfiieYl3sptWzLhgw9gxw7o0cPsSCI11OkIy5w5c4iNjSU4OJjExETWXuIQ4rvvvkvPnj0JDg6mb9++fPLJJzW+fvfdd2OxWGosY8aMqUs0kUZx35r3APhPrxEcDGtnchqRy7crvCPcfHPNsVkibsTlwrJkyRJSU1OZNm0a2dnZ9OvXj+TkZA4fPnze7VetWsUdd9zBPffcw8aNGxk/fjzjx49ny5YtNbYbM2YMhw4dci7vvPNO3T6RSAPrfOwg1+9YBcDribeanEakHk2ZYjwuWgQHDpibReQnXC4ss2bNYtKkSaSkpNC7d2/mzZtHSEgICxYsOO/2f/7znxkzZgyPPfYYvXr14rnnnqN///689tprNbYLCgoiMjLSubRqpenNxT09/O3b+GFnWdfB/NA21uw4IvVn6FAYORLKy+FPfzI7jUgNLhWWsrIyNmzYQFLSmZu7+fn5kZSURFZW1nnfk5WVVWN7gOTk5HO2X7FiBe3ataNHjx5MnjyZo0ePXjBHaWkpNputxiLSGLoV7uPm778CYPawO01OI9IAnn3WePzb32DfPnOziJzFpcJy5MgRKisriYiIqLE+IiICq9V63vdYrdZLbj9mzBgWLVpEZmYmM2bMYOXKlVx//fVUVlaed5/p6emEhYU5l5iYGFc+hkidPVJ9dOXT7lexNaKL2XFE6t/IkXDttcZRlueeMzuNiJNbXNZ8++23c/PNN9O3b1/Gjx/PRx99xLp161ixYsV5t09LS6O4uNi57N+/v3EDi2/KyWHsjm+pwsLLw35ldhqRhuMoKm+8Abt2mRpFxMGly5rDw8Px9/enoKCgxvqCggIiIyPP+57IyEiXtgeIi4sjPDycXbt2ce21157z9aCgIIKCglyJLnL5pk0D4KNewzV2RbxS7FmX5/8jbgDX7NnA++PvJfXGRwHInT7WrGgirh1hCQwMZMCAAWRmZjrXVVVVkZmZydChQ8/7nqFDh9bYHmDZsmUX3B7gwIEDHD16lKioKFfiiTScr7+G//yHSosff776DrPTiDS4WcN+DcD4rSvoXbDH5DQidTgllJqayvz581m4cCHbtm1j8uTJlJSUkJKSAsDEiRNJS0tzbv/www+TkZHBSy+9xPbt23nmmWdYv349Dz74IAAnTpzgscceY/Xq1eTm5pKZmcm4cePo2rUrycnJ9fQxRS5DVRWkpgKwuN9odrfRmCnxfpujuvGfXiPww86Ty/9mzM8iYiKXZ7qdMGEChYWFTJ06FavVSkJCAhkZGc6BtXl5efj5nelBV111FW+//TZPPfUUv//97+nWrRtLly6lT58+APj7+7Np0yYWLlxIUVER0dHRjB49mueee06nfcQ9vP02rF8PLVrwsq4MEh/ywsi7SP4hi6v3beLa3WuBG82OJD7MYrd7fm222WyEhYVRXFxMaGio2XHEm5w8aUxRfuAApKcTW9TX7EQijerxlW9w/+r32N26A12se6BJE7MjiRdx5fe3W1wlJOK2Zs0yykqnTvDII2anEWl0fxnyS46EhNHl2AGYO9fsOOLDVFhELmTPHvjjH43n06dDcLC5eURMcCIohFnDjQG4PP00HDpkbiDxWSosIudjt8P998Pp08YkWhMmmJ1IxDSL40eTE9UNbDYdaRTTqLCInM+//gWffQaBgfCXv4DFYnYiEdNU+fnzZPKD4Odn/L+RkWF2JPFBLl8lJOKtHJNmtSgtIXP+ZNoBLw+6jT8v2AnsNDWbiNm2RnSBhx+Gl182jj5u2QIhIWbHEh+iIywiP/F05nzalfzI7tbtmTvkF2bHEXEff/gDdOgAe/fCU0+ZnUZ8jAqLyFmu27maX27+giosTBnzEGUBuoRTxKl5c3j9deP5yy/D8uXm5hGfosIiUq1NSRHpGa8C8NfEn7Mupo/JiUTc0A03wL33Gs/vvhuKi02NI75DhUUEwG7nT5+9RvjJYraHd3LeR0VEzuOllyAuDvLyjHEtIo1Ag25FAObMIXnnasr8Aki98VGdChI5j7Pv5jxg6H28u+cJ/BYu5P+OtOGDPj8DdEdnaTg6wiKyZo3z5oYzRt3N9xFxJgcScX8bOvTmlatvB+BPn82he2GuuYHE66mwiG87cgR+8QsoL+eT7lfx94HjzE4k4jFeuep2voq9kqYVpcxdmk7z0pNmRxIvpsIivqu8HH71K9i/H7p14/EbHtEEcSIuqPLz55Gbfkd+i3C6HDvIzE9ehqoqs2OJl1JhEd9kt8NDD8GyZcbkV++9x4kgTYIl4qpjIWE8OO4JyvwCGPNDFqSlmR1JvJQKi/imWbOM+SQsFnj7bYiPNzuRiMfKbt+Lx2+ovlrohRdg/nxzA4lXUmER3/Puu/DYY8bzl16CcRq3InK5ll5xDbOvvsN4MXkyfPqpuYHE61jsdrvd7BCXy2azERYWRnFxMaGhoWbHETcVO+VjfrZrLa9/8EeaVFWy6MqxTL3uPo1bEakvdju5+UvgrbcgOBg++QSuucbsVOLGXPn9rSMs4jOG7d3I3KV/oklVJR/2GskzSfeqrIjUJ4sF/v53uOkmOH3aePz2W7NTiZdQYRHfkJHB/PefJ6iygozuQ3l07P9R5edvdioR7xMYCP/6F1x3HZSUGFP5f/WV2anEC2imW/F+77wDEyfStKKCzC6DeOjmx6nw1199kYbgmA03uO//8o8dVobmbeb0tdfxwLgnyOyaCGg2XKkbHWER7/bqq3DnnVBRwX96jeC+W35Pub+m3RdpaKebBHP3bc+wrOtggivKeP39P3Lb5i/MjiUeTIVFvFNZGdx3H/z2t8acKw88wMM3/U5lRaQRlTYJ4r5bnuS9PtcSYK9i5iezSVu+ACorzY4mHkiFRbxPQQEkJZ2ZZyU9HV59FbtFf91FGlulnz+P3fAwrw6dAMD/rn0fxo6FH380OZl4Gl3WLF7Bcd58xJ4NvPTJy7QtKcIWGMLDNz/G8i6DTE4nIgBjt33Ni5/OJqS8FGJijMufR4wwO5aYSJc1i88JKi/l6cz5LHp3Gm1LitgR3pFbJr6ksiLiRj7uNZzb7nwRunY17uE1ahQ8+aRxClfkElRYxPMtW8bnCx7gnvUfAvBG/xu5eeLL7G4TY3IwEfmp7yPiYONGSEkxxpf96U+QkABff212NHFzKiziufbtM+62PHo0nYqsHGrehpTbpvHMdfdR2iTI7HQicgGxz68ktt2t3D9uCoUhLWHbNhgxgn/FX8fgBxaZHU/clMawiOc5dswYSPvKK8ahZD8/Flx5Iy8N/zUluuOyiEcJO3WcJ1a+wa+++wyAUwFBNJ3ymHG/L/0893qu/P5WYRHPsX8/vPwy/PWvxgyaAD/7Gbz4IrH/OmRuNhG5LP0PbuP3yxcw8OA2Y0VYmHETxd/+FqKizA0nDUaFRbyH3c7PJ77EnTmfcvP3K2lSZczfsLVdHC+MvIuVnfvrfkAi3sJuJ3lnFq9v/8A4TQTGVP8TJxrzKvXX/+/eRoVFPF9+vnHJ44IFsGOHc/WqjvHMS7yVr1RURLxW7p+uh48+ghkzYNWqM1/o29cYrHvnndCunXkBpd6osIjnsduNYrJ0qbGsWeP80skmQXzSYxiL+o9lU1R30yKKSOMbeGArv974CdfvWEVQZTkAVVjwG3Y1jB8P48YZl0mLR1JhEfdnt8O+ffzugT8zJG8zQ/I20cFWWGOT9e178a++1/Fxz2EaTCvi40JPn+DmbV9x2+ZlJBzaWfOLcXHGnC7XXAMjR0KHDjoC6yEavLDMmTOHF198EavVSr9+/Xj11VcZPHjwBbd/9913efrpp8nNzaVbt27MmDGDG264wfl1u93OtGnTmD9/PkVFRVx99dXMnTuXbt261SqPCoubKyuDPXtgyxbYsMFYsrPh6NGam/kFkNUpns+7DWFZ10QOt2hjUmARcWfRtsNct3MN1+1czZC8zQTYq2p8vTCkJVsjujDqV2NgwAC44gro3Bma6F5i7qZBC8uSJUuYOHEi8+bNIzExkdmzZ/Puu++yY8cO2p3nnOKqVasYMWIE6enp3Hjjjbz99tvMmDGD7Oxs+vTpA8CMGTNIT09n4cKFdO7cmaeffprNmzfz/fffExwcXK8fWBpAWZkx5uTAgTPL/v2waxd7szYSU1Rwzg8UgHI/fzZHdmV1x75kdYxnffvenAq89H9vERGH5qUnGXhgK0PyNjM0bzNXFOw+788bAgKMIzHduxunkGJioH1742hMhw4QHa1CY4IGLSyJiYkMGjSI1157DYCqqipiYmJ46KGHmDJlyjnbT5gwgZKSEj766CPnuiFDhpCQkMC8efOw2+1ER0fz6KOP8rvf/Q6A4uJiIiIieOONN7j99tvr9QP7LLsdKiqMpbzcWM5+fvq0canwxZaiIuOoyNnLsWPG+ks4EdiU3a07sDWiC1siu7A5ois72sZSFqAfECJSf4LKS+lVmEufgt30se6iT8Fu4o4dMO5fdDEWC7RsCW3aQOvWxqNjCQuD5s2hWbMzy9mvmzY1rmYKDDRKz9nPmzQBP83ReiGu/P4OcGXHZWVlbNiwgbS0NOc6Pz8/kpKSyMrKOu97srKySE1NrbEuOTmZpUuXArB3716sVitJSUnOr4eFhZGYmEhWVtZ5C0tpaSmlpWf+8hUXFwPGB69X5eXGXUXtdmOBms9/+tqV5668x9X9nV1GHI8NfDv3Ur8ArC3aUNi8NdbmrSlo3ob9LSPZ1yqSva3acySk5bnnlCvLjUVEpJ6cArLbxJDdJgZ6jzJW2u20O3GM2KJ8Yn88RIfiAiKOHyPixFFjOX6MQHulcQfphriLtL//mSLj72/8LPTzu/jy023Ofm2x1Px5+tOfrZd6/VO1fX9gIHz8ce0/dy04fm/X5tiJS4XlyJEjVFZWEhERUWN9REQE27dvP+97rFbrebe3Wq3OrzvWXWibn0pPT+fZZ589Z31MjO4dY5qqCiguMBYRETezH9hg1jevrDSW06fNSlB/wsIaZLfHjx8n7BL7dqmwuIu0tLQaR22qqqo4duwYbdq0weLDI8NtNhsxMTHs379fp8YamP6sG4/+rBuX/rwbj/6sjSMrx48fJzo6+pLbulRYwsPD8ff3p6Cg5r+iCwoKiIyMPO97IiMjL7q947GgoICos6ZfLigoICEh4bz7DAoKIiio5s3tWrZs6cpH8WqhoaE++5e/senPuvHoz7px6c+78fj6n/Wljqw4uDQSKDAwkAEDBpCZmelcV1VVRWZmJkOHDj3ve4YOHVpje4Bly5Y5t+/cuTORkZE1trHZbKxZs+aC+xQRERHf4vIpodTUVO666y4GDhzI4MGDmT17NiUlJaSkpAAwceJE2rdvT3p6OgAPP/wwI0eO5KWXXmLs2LEsXryY9evX89e//hUAi8XCI488wvPPP0+3bt2clzVHR0czfvz4+vukIiIi4rFcLiwTJkygsLCQqVOnYrVaSUhIICMjwzloNi8vD7+zLuG66qqrePvtt3nqqaf4/e9/T7du3Vi6dKlzDhaAxx9/nJKSEu69916KiooYNmwYGRkZtZqDRc4ICgpi2rRp55wuk/qnP+vGoz/rxqU/78ajP2vXeMXU/CIiIuLdNJuNiIiIuD0VFhEREXF7KiwiIiLi9lRYRERExO2psPiA0tJSEhISsFgs5OTkmB3H6+Tm5nLPPffQuXNnmjZtSpcuXZg2bRplZWVmR/MKc+bMITY2luDgYBITE1m7dq3ZkbxOeno6gwYNokWLFrRr147x48ezY8cOs2P5hOnTpzun95CLU2HxAY8//nitpj2Wutm+fTtVVVW8/vrrbN26lZdffpl58+bx+9//3uxoHm/JkiWkpqYybdo0srOz6devH8nJyRw+fNjsaF5l5cqVPPDAA6xevZply5ZRXl7O6NGjKSkpMTuaV1u3bh2vv/468fHxZkfxCLqs2ct9+umnpKam8u9//5srrriCjRs3XvCWB1J/XnzxRebOncuePXvMjuLREhMTGTRoEK+99hpgzKwdExPDQw89xJQpU0xO570KCwtp164dK1euZMSIEWbH8UonTpygf//+/OUvf+H5558nISGB2bNnmx3LrekIixcrKChg0qRJvPnmm4SEhJgdx6cUFxfTunVrs2N4tLKyMjZs2EBSUpJznZ+fH0lJSWRlZZmYzPsVFxcD6O9wA3rggQcYO3Zsjb/fcnEeebdmuTS73c7dd9/Nfffdx8CBA8nNzTU7ks/YtWsXr776KjNnzjQ7ikc7cuQIlZWVzlm0HSIiIti+fbtJqbxfVVUVjzzyCFdffXWNGcml/ixevJjs7GzWrVtndhSPoiMsHmbKlClYLJaLLtu3b+fVV1/l+PHjpKWlmR3ZY9X2z/psBw8eZMyYMfziF79g0qRJJiUXqbsHHniALVu2sHjxYrOjeKX9+/fz8MMP889//lO3n3GRxrB4mMLCQo4ePXrRbeLi4vjlL3/Jf//7XywWi3N9ZWUl/v7+3HnnnSxcuLCho3q82v5ZBwYGApCfn8+oUaMYMmQIb7zxRo17aonrysrKCAkJ4b333qtxI9S77rqLoqIiPvzwQ/PCeakHH3yQDz/8kK+++orOnTubHccrLV26lFtuuQV/f3/nusrKSiwWC35+fpSWltb4mpyhwuKl8vLysNlsztf5+fkkJyfz3nvvkZiYSIcOHUxM530OHjzINddcw4ABA3jrrbf0A6eeJCYmMnjwYF599VXAOF3RsWNHHnzwQQ26rUd2u52HHnqIDz74gBUrVtCtWzezI3mt48ePs2/fvhrrUlJS6NmzJ0888YROw12ExrB4qY4dO9Z43bx5cwC6dOmislLPDh48yKhRo+jUqRMzZ86ksLDQ+bXIyEgTk3m+1NRU7rrrLgYOHMjgwYOZPXs2JSUlpKSkmB3NqzzwwAO8/fbbfPjhh7Ro0QKr1QpAWFgYTZs2NTmdd2nRosU5paRZs2a0adNGZeUSVFhELtOyZcvYtWsXu3btOqcM6gDm5ZkwYQKFhYVMnToVq9VKQkICGRkZ5wzElcszd+5cAEaNGlVj/T/+8Q/uvvvuxg8kch46JSQiIiJuT6MCRURExO2psIiIiIjbU2ERERERt6fCIiIiIm5PhUVERETcngqLiIiIuD0VFhEREXF7KiwiIiLi9lRYRERExO2psIiIiIjbU2ERERERt6fCIiIiIm7v/wMNIR5AFp0xwQAAAABJRU5ErkJggg==\n"},"metadata":{}}]}]}