{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Poisson Regression\n", "\n", "Gaussian process models can be incredibly flexbile for modelling non-Gaussian data. One such example is in the case of count data $\\mathbf{y}$, which can be modelled with a __Poisson model__ with a latent Gaussian process.\n", "$$\n", "\\mathbf{y} \\ | \\ \\mathbf{f} \\sim \\prod_{i=1}^{n} \\frac{\\lambda_i^{y_i}\\exp\\{-\\lambda_i\\}}{y_i!},\n", "$$\n", "where $\\lambda_i=\\exp(f_i)$ and $f_i$ is the latent Gaussian process.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ " \n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1zUVf748TPAMNxG7l4AxQAVDUHFaiW1vELrrRRLNC/ZtlvKd7vYV83KtrJfl93tq33zYVu7munPtjUlb0l5w7xEmyauaOKmpnIxuQ6gw9yY7x+wqWCCA8wBzuv5R484zpx5M36YF/OZATV2u10AAKAqF9kDAAAgEyEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACitpUJoNBpzcnIaf3mbzVZdXd1Cw7QJdrvdarXKnkIyDgMOAyGE1WpV/Fc/chgI5x4GLRXCnJychx56qPGXr6ysrKqqaqFh2gSbzVZaWip7CsnKy8vNZrPsKWSyWq1lZWWyp5CMw8BisRgMBtlTSGYwGCwWi3Nui1OjAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBNqY/Pz8OXPmDBgwYODAgU899VRhYaHsiYC2jRACbcmxY8f6xPZbcd7vyNj3Die9u+yE5va+sadPn5Y9F9CGuckeAMAtWLBggWHYPJH4TO3H3eML3T1ffPHFdevWSZ0LaMN4Rgi0GXa7fc+ePWLQtOtWfzVtz549kiYC2gNCCLQZdrvdarMJd6/rVnXeiv8eFqCJCCHQZri4uPSLixPHd1y3mv3FgAEDJE0EtAe8Rgi0JS+//PK4qY9Waz1E3yRht4uszdrPXli8daPsuYA2jGeEQFvy61//On39R/33v+zy+yDXJ4Pu+u7t3dvShgwZInsuoA3jGSHQxowaNWrUqFFXrlxxcXHx8PCQPQ7Q5hFCoE3y8vJq+EIAGoFTowAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEpzPIRWq3XOnDnBwcF33313Xl5eM84EAIDTOB7CpUuXlpeXnzt3LiEh4aWXXmrGmQAAcBo3h6+5bt26VatWeXl5LV68+NSpU/UvcPny5U8++aT+ekRERN++fessVlVV2Ww2Fxd1T9VarVaTyVRVVSV7EJlMJpNGo5E9hUwWi4XDwGQyubq62u122YNIYzabOQxMJpNWq62urm7k5bVaraurq2O35XgIz5079/HHHw8bNiwiImLVqlX1L1BZWbl27dr664mJiVFRUXUWa0Ko8oOg1Wo1Go1Go1H2IDLVfPoqPwJaLBYOA6PRqNFoGv8I2P6YzWYOA6PR6OrqarPZGnl5V1dXCSEsLy+32+3Hjx9fvnz5Y489lpmZWecCnTp12rJlSyN3c3Fx0Wq1Xl5eDs/T1lmtViGEv7+/7EEk8/T09PDwkD2FNBaLxcXFRfHDwG63e3t763Q62YNIYzab3dzcFD8Mqqur9Xq9u7u7E27L8VORwcHBTz31VJcuXVJTU7Ozs5txJgAAnMbxECYmJn744Ycmk+n9998fOHBgM84EAIDTOB7C119/fffu3Z06ddq1a9df//rXZpwJAACncfw1ws6dO+/YsaMZRwEAwPnU/XEFAAAEIQQAKI4QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNDfZAwDtXEFBwQcffJCTk9O5c+eJEyfefffdsicCmtPRo0fXrVt34cKFiIiI2bNnR0RENHHD0tLS999//9ChQ506dRo/fvzo0aObZc6baNIzwuzsbG9v7+YaBWh/0tPTo2PiXvqqcJ1X0tvngoZMSJk3b57soYBm89Zbb8Xfk/hWjufH3ve99p3l9gF3rV27tikbZmZm9ozuvXDbqU99xy4vui1xRurMmTOrq6uba+AbcvwZocFgmDVr1pUrV5pxGqA9MZlMsx55pHz6ShGTWLNiHzL77VcGjhkzZvjw4XJnA5ru+PHjz736ZvUL34jA8JqVqviJj6fel5iYGBwc7MCGdrt91qxZRWNfF4Merl0a+puPXktI+uSTlJSU5hq7PgdDWDPuwoULJ0+e/EuXMZlM3377bf31jh07hoSE1Fm0WCw//1dNVqvVYrGofA8IISwWi5ubm6urq+xBmsf+/ft/0nb8uYJCCOEdIAY/kpaWNmTIkBtexfIfThqxVaq5B1xc1H0HQ1s5DDZt2lQdP+nnCgohRPiAy13v/OKLLx566CEHNszJyckpKLtaQSGEu5cY9kRaWlpycvLNr+vq6urwMeNgCN98883IyMibT1ZYWDh79uz668nJyXPnzq2zWF5ertVqW/9ffMuxWq0VFRXu7u6yB5GpvLzcYrGYTCbZgzSPvLw84dul7qpfSGHhKYPBcMOrWCyWiooKrVbb4sO1YuXl5TabTeWvBbPZXFlZ6ebW2t/DcenSJeHbue6qX5f8/PxfOsJv7sKFC8Kv3peMf2jR8aIGN+zQoYPDx4wjd/SePXvS09N37Nhx84uFhYVlZWU1ck+tVqvVar28vByYp32wWq2urq5BQUGyB5HJ1dXV09PTw8ND9iDN44477hAXFgirSbjprq6e+WffpL6/9BdtsVi0Wq3ih4GLi4u3t7dOp2v4ou2U2WzW6XSBgYGyB2lAXFyc2PWP65bsdnH2UHz8w44dwwMHDnT56VS1sVx4dri6euafMTExLfpF4cgTyV27du3du9fd3V2j0QghNBrN/v37m3swoM3r3bv36Ltixf//L2Ex1i59u973+KZHHnlE6lxA80hOTg4tPSZ2Lxd2uxBC2Cwi7cV+wW733nuvYxsGBwdPe3CSWPWoMJbXLh3foftqxeOPP948E/8CR0K4ZMkS+38IIex2++DBg5t7MKA9+Pjjj6d2M2sW9hB/ThQvxsRmvvnF9s/rv0YOtEV6vX7Hjh13n/9UPNdTvH2fWBg1xuX4li1bmvIy/3vvvff4XV1cn+8l/jRK/GFA1NYnt6Z92qdPn2Ycuz5NTcwcv77mxjtkZWXNmjWr8adGDQYDp0ZLS0sde6tVu1FaWtqeTo3+LD8/Pycnp0uXLj169Lj5Y4TFYjEYDIqfGi0pKeHUaEVFRes/NVrDbrefPn36/PnzkZGR4eHhDV+hEQoLCzMzM0NDQ2NjY53wWmlTb6CJHQVUEBISwrNAtFcajSYqKioqKqoZ9wwODk5ISNDr9c55x5C6b1AGAEAQQgCA4gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlOYmewBhtVo//vjjr7/+2tvbOzExceTIkU3fc9++fV988UVFRUVcXNy0adN0Ol3T94QisrOzN23aVFBQ0LNnz+nTp/v7+8ueyBkOHz68devWoqKiPn36TJ8+3cfHR/ZEaCm7du3atWuX0Wjs379/SkqKVquVPZF8kp8R5ubm9u/ff8YbH60ojvjTad9R01OTk5PNZrPDG9psttmzZw+dOOO149p3LnV/9J2NMTEx//73v5txZrRjr732Wv8hI184UL68tOeTnxzq1btPRkaG7KFa3H//93/fOWrCHw5Z3i2OmrNqT3Tv3ocPH5Y9FJqfxWJ58MEHRz485/WTnksLwma+tbZfv34XLlyQPZd8Grvd3hL7ZmVlzZo1Kysr6+YXu//++zcZI8TkN2s/tprFn0f/6fGJ8+bNc+x2V65c+ehLS8Vz+4S7V+3StteHluzau3evYxs6jdVqLS0tDQ4Olj2ITKWlpZ6enh4eHlJuPTMzc9Do8WLxt8IvpHbpu89CNj115swZp51UsFgsBoMhKCjIOTcnhNi2bdvYmXPEi/8U3gG1S/tWRv9z6YkTJzQajdPGuFZJSYm3t7fKJ3LMZnNFRUVgYGDzbrt06dKn3/1EPLtTuLnXLm1YNNbt5JYtW5r3hppFcXGxXq93d3dv+KJNJvMZodFo3Pb5djF20dUlN3dx3/wNGzY4vOeGDRtE4ryrFRRCJD27L/PQpUuXmjAplLBhwwYxZPbVCgohBtyf7xKYmZkpb6gWt2HDBjF87tUKCiGGzD55sfz48ePyhkKL2LBhg7hv/tUKCiHGLtqe/sXly5flDdUqyAyhwWCwuvsIT9/rVgPDi4qKHN6zqKhIBHa7bslVa/ftVFxc7PCeUERRUZEI6FZ3NaBbUw7I1q+oqEgEhtddDWznn7WabnCE63xsnn5lZWWSJmotZIYwKCiog6tVFJ+7bvX8kcjISIf3jIyMFOePXLd0ucS94mJYWJjDe0IRNzh4qm0i918RERGSJnKGG3zWliqR/337/qzVdIO/65JcH7uxY8eOkiZqLWSG0M3N7bHHHhOrfyeu/Of7kZ/+LTa+8MQTTzi85+9+9zuXz98QucdqPzZVitW/fXhqil6vb/K8aOdmzJjhnfWpOLa99uNqq/j0uUG9wvr16yd1rpb16KOPave+J344UPux1STW/X7M8MHdutV7cow27oknnhBpi8XFU7UfGw1i9W9/85vf8MZRYW8ZR44ciYuLa/BiJpNp7ty5rr4dRdwY0XuEb2DwihUrmnjTa9asCerYWfS6R/Qbp+nQadasWZcvX27ink5gsVguXbokewrJSkpKjEajxAF2794dGRkpbrtDDLhfBHVPTEzMy8tz5gBms7mwsNCZt2i32zdv3hzWtauIShD9Jwj/0EmTJhUXFzt5hmsVFxdXVVVJHEA6k8lUVFTUEjv/5S9/8QvqKHoPF3FjXH07PvHEEyaTqSVuqOmKioqcNpvkd43WyM3NPXDggF6vT0hI8PPza/qtV1RUHDp0qLy8PC4urnv37k3f0Al416iQ/a7RGiaT6bvvvisoKOjVq9ftt9/u5Ft3/rtGaxiNxsOHDxcWFsbExPTo0cPJt14H7xptoXeN1jAYDIcOHbpy5Uq/fv26du3aEjfRLJz5rlH5P1AvhAgLC0tKStJqtV5eXg1fuhH0ev2wYcOaZSuoRqfTDRo0SPYUzubp6Tl48GDZU8AZfH19R4wYIXuK1oVfsQYAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0x0O4adOmmJgYPz+/oUOHnjp1qhlnAgDAaRwM4fnz5x9++OEPPvigoKBg/PjxjzzySPOOBQCAc7g5drUzZ85MmTJl0KBBQoiZM2e+8cYb9S9TXV1dUlJSf12n03l6eta/cA3H5mkHuAcEdwL3gBCCO4F7QAhx63eCi4vjJzg1drvd4SsLIWw2W2pqqouLy/Lly69dz8rKGjZs2A0/h5kzZz7//PN1FsvLy7Vabf1AqsNqtRoMhsDAQNmDyGQwGDw8PHQ6nexBpLFYLBUVFQEBAbIHkamsrMzLy8vd3V32INJYLJbKykp/f3/Zg8hUWlrq4+Oj1WobeXk/Pz+HHzocfEZYY+fOnfPnzx89evSSJUvq/2l4eHhWVlYjt/Lw8NBqtV5eXk2Zp02zWq3u7u7BwcGyB5HJ3d3d09PTw8ND9iDSWCwWDw+PoKAg2YPIpNVqvb29Vf5+yGw2e3p6Kv5tsZubm16vd873Qw6G0G63L1q06MCBA3//+9979uzZvDMBAOA0Dobw4MGDaWlpmZmZbm5ulZWVQggfH59mHQwAAGdwMIQZGRk5OTnXnsJu4muNAABI4eDbbJ5//nn79Zp3LAAAnIPfLAMAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaW6yB2gbKisrt2/ffu7cufDw8Pvuu8/Hx0f2RECrZjKZPv/889OnT4eFhY0ePTogIED2RM5w+fLl7du3//jjj+Hh4UlJSXq9XvZEaBRC2LBdu3bNmDkzP7Cf6NxTXMwIfeaZj1avHj58uOy5gFbq22+/TUlJOe3eXXSNFUUHAn//1Ad/WfHAAw/InqtlZWRkTJ8xI9e/r+jcS/yUEfL0M6s/XDVy5EjZc6FhhLABxcXFyZMfLJu+SvS9r2Yl71/bkic/+MO/TynyTS5wS6qqqiZOnJg76mWRML1mpfj011NnjD+ZPSA8PFzubC2ntLR0UvLkkpT3RdzYmpX87PSaB4qgoCC5s6FBvEbYgK1bt5Z1HfRzBYUQInZMaeidW7dulTcU0Hrt2bMn163zzxUUQojIQVVx969fv17eUC1u27ZtJV3if66gEELEJBnC796yZYu8odBYhLABeXl5okuvuqude+Xn58sYB2jt8vPzb/glk5eXJ2McJ8nPzxed633WXaLb92fdbhDCBoSGhoqCnLqrF3NCQ0NljAO0dmp+yYSGhoqL9T7rgpPt+7NuNwhhA8aNG+efmymObb+69K9tAfnfjh079pevBKjr3nvv7VZ9SRxcc3Xp9NceRz+bPHmyvKFa3JgxYwIvfieOXvOKSXa63/mD48ePlzcUGos3yzQgICBgw6frp8+YkZcRJzr3FBdPhZX+a82n6/39/WWPBrRGHh4eGzdunDJlyg+Z60TXWFF0Nig3869rV7fjd8oIIfz8/GoeKC589VfRuZf46VRo8b9Wr/9HYGCg7NHQMELYsGHDhuWcPFnzc4Tdu9+dlJTk7e0teyig9YqPj8/Ozk5PT//hhx/Cwu4aPXqlCt843nPPPd+fOJGenv7jjz+GhyckJSXxA8dtBSFsFG9v7+TkZNlTAG2GTqebMGGC7Cmczdvbe9KkSbKnwC3jNUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASnOTPYC6jhw5cvTo0Q4dOgwePLhjx45N3zA3N/frr782Go3x8fG333570ze8cuXKvn37Lly4EBkZOWTIEDe3ZjhavvnmmxMnTgQEBAwZMiQgIKDpG549e/abb76x2WwDBw7s1atX0zdEq3X8+PHDhw97eHgkJCSEhYU1fcNLly7t37+/vLw8NjZ2wIABTd+wJWRlZWVlZen1+sGDB3fq1En2OO2U3VElJSVjx4719/cfN25cSUlJnT89cuRIXFxc43crKyu7fPmyw8O0LSUlJePGjRNB3UXCdNF/vI9/0LvvvmuxWC5duuTwnq+++qqHb5CInyh+NVUTEDZ16tQm3p87d+7s2q2biL5XDJ4lusfHxsYeO3asKRsWFBQMHz5cdO4pEmaI2DH+QR1Xr15d5zIlJSVGo7GRG9pstvnz52t9O4o7Jou7prj4dfntb39rNpubMqR0ZrO5sLBQ9hSSFRcXV1VVXbty5cqVadOmafxDxa+miviJOt+gV155pYm3snz5ch//INF/vEiYLoJvGzNmTP3HMVlMJlNRUVFZWdmECRNEYLgY9LDoP8HbL2jZsmWyR3OeoqIik8nknNtyPIQLFiyYO3duVVXV3LlzFy5cWOdPCeFNTJ48WQyaJlZUig/M4gOz+H8nRYdO27ZtcziEa9euFcER4o/najdcXiZix8yZM8fhCfPz8306+IonN9du+IFZTHk7IiKizsPTLRkxYoQYPlf8xVi74R++c/H2/+abb669zC2F8J133hFhfcX/5Ndu+E6R6DnkhRdecHjC1oAQ2m8UwtTUVBH7a7G8rPbv+k/nRMfIjz76yOGb+PLLL4W+o3jtRO2GKypFwvRJkyY1efbmURPClJQUcdcUsaLiPw8UOcK3y7Zt22RP5yTODKHjrxGmpaWlpqbqdLrU1NSNGzc245PU9q2srGzjpi1i6jvCzb12KThCjHlu9erVDu+5cuVK8cCrwq9L7cfuXmLa/3700UdWq9WxDf/xj39U9r5PxCRdXRqResbqm5GR4diG586d23XwkHjwLeHiWrsUGlM94r8+/PBDxzYUNZ/15DeFT1Dtx54dxJT/WblypcMbonWy2WyrV68WU98R7l61S75dxANLmvJ3vWrVKjFmoegYVfuxm7uY+s5nW7aVlJQ0ed7mUVFR8enGNDHtf4WbrnYp+DYx7nmO8Jbg+Ks+eXl54eHhQojw8PCCgoL6Fzh37pyfn1/99ZkzZy5atKjOYnl5uVar9fT0dHietuLkyZM2vzDhob9uNazvD7s+KiwsrK6udmDPM2fOiHtirlsKCKusdjt58mRwcLBjQ4rQmLqrYTHZ2dn9+vVzYMOsrCzRuZdw1V6/YeypnBU//fTTzwsGg8HDw0On09W9/o2cPXtWPHz9kKF98gsu5ubmarXaX7hSa2exWCoqKmw2m+xBZKo5P+TuXvudYnFxcYVVIwK7XXehsJgz289ce/Dckh9++EHcM/u6JZ23LaDb0aNH+/Tp49iezchisZw4ccKi7yw8fa/7g7C+P2z/wOHPum0pLS01Go2N/1r28/Nr5ENHfY6H0G63azSamv+54ddt165d9+zZU3/dw8OjfvDc3d21Wq2Xl1f9y7c/mrI8u9V89RmhEOLS6bCwsMDAQMe6FRoa+mPhGRHS++pSZbGuuioqKurnR5Nbctttt4njOXVXL52Jiprk2IS9evUSRWeF3S40mqurhae7det27YZubm6enp4eHh6N2TMkJMRQeEb4dr66VPRjYIB/SEiIAxO2EhaLxd3dPSgoqOGLtl+urq7e3t4/P6j5+fl52M1VFYVCf82xd+l0aGioY0ejECIsLOzbS6dFr3uuLtksmpLc6Ohoh/dsRmazOSoqyqX8YrWlSmiv+XK4dDosLKw1TOgELi4uer2+8Y9gLi6On+B0PIQhISEXLlzo0aNHXl5eaGjoDccKDAxs5G4u/+HwPG1Fp06dRt4zeMfmV8TEJbVLV0rF9remLH3F4XsgJSXlwB9fFb2H1Z47stvFxucnTZzYyKLUl5yc/OIrA8zD5oiusbVLR7d2qjw9cuRIxyaMjo6O79nt8Jdvi8R5tUuGi2LHsqnrV1+74S0dBikpKYvXvSye3FL7LUW1TWx8ISUlpU0fRep8IdxEnTtBp9MlJyev3fiCmLFCaFyEEMJ8RWxZkvLMDIfvqKlTp6bNXSjiHxBe/rVLm18dPmRQly5dbno9J6l58EwcMWz7pj+I5DdqV40G8fmbU/+8WJHDw6lfCw6/uvj0008vXLiwurp64cKFzz77bJ0/5c0yN3HhwoX+/fuLHoPF+MVi1FOu/iHz589vyrtGrVbrY489pgnsJpKeFeNeEN0HDh48uKioqClDrl27Vh8QLAY/Iu5/WQxM7hISmpGR0ZQNc3JyoqOjRZ+RYsJLYkSq1rfja6+9Vucyt/RmGZPJ9NBDD4mOUeK++WLMcyKsb1JSUkVFRVOGlI43y9hv9GaZ4uLioUOHiu7xYtwLIulZTVD4o48+arVam3IrCxYscPXrIkY9JcYvFj2H9OvX7/z5800bvNnUvFkmNzc3Pj5eRN0txr8oRj/t4h/yzDPPVFdXy57OSZz5ZhmN3W53rKBlZWXTpk07evTogAED1qxZ4+t73bnsrKysWbNmZWVlNXI3g8GgzqlRIYTNZtuyZUtWVpavr++oUaNiYmKsVmtpaWlTTnocPnx49+7dVVVVAwcOTEpK0lx7EtIhBQUFmzdvzs3NjYyMnDRpkl6vb/g6N2WxWD777LNjx44FBQUlJSX17NmzzgVKS0sbf2q0xsGDB/ft22ez2e66664RI0Y0cULpLOZEsHgAAAReSURBVBaLwWBQ/NRoSUnJtadGa9jt9vT09MOHD+t0umHDhg0cOLDpN3T8+PEvv/zSYDDExcWNHz/e1dW14es4hdlsrqioCAwMrK6urnmg0Ov1o0aN6tu3r+zRnKe4uPiWTo02SQsF9lafEb733nu7d+9uoWHahPz8/Hnz5smeQrJly5YdPHhQ9hQynT17dtGiRbKnkOytt946fPiw7Clk+v77719++WXZU0j2yiuvZGdnO+e2Wsu55q+++urUqVOyp5DJYDDwUyi7d+8+c+aM7ClkKioq2rRpk+wpJNuxY8f58+dlTyHTxYsXt23bJnsKydLT02/48wgtobWEEAAAKQghAEBphBAAoDTH3zV6c99//31ycnJ8fHwjL3/o0KHAwMDbbrutJYZpEyoqKg4ePJiYmCh7EJkyMzNDQ0O7du0qexBpysrKDh06NHLkSNmDyHTgwIGIiIhW8iN9UhQVFWVnZ997772yB5Fp7969vXv3bvy/zPPkk082vjh1tFQIhRAZGRmKv+INAHCOYcOGOfw9dAuGEACA1o/XCAEASiOEAAClEUIAgNIIIQBAaa0ohOnp6X369PHz8+vTp8+XX34pexwJNm3aFBMT4+fnN3ToUMV/4ZzNZouOjpY9hbOVlpaOGzcuICBg/PjxpaWlsseRSc0DoAaPA8L5OXDOrzRtkM1mCwgI2Llzp81mW79+fUhIiOyJnO3cuXM+Pj4HDx68cuXKH//4x4SEBNkTSbN06dI777yz9RycTrNgwYK5c+dWVVXNnTt34cKFsseRRtkDwM7jgN1ul5GD1vKM0Gq1rlmzZvjw4ZcvX9bpdH5+frIncrYzZ85MmTJl0KBBnp6eM2fOzMmp9w/EKyM2NvbFF1+UPYUEaWlpqampOp0uNTVV5d/AruwBIHgcEELIyEHr+jnCyspKvV6v0Wj279+fkJAgexw5bDZbamqqi4vL8uXLZc8ik0bTug5OJ/Dx8SksLPT09DQajZ06dSovL5c9kUwKHgDX4nHAmTmQ+YwwOjpao9Fc++/H+vj4VFZWLlmy5Mknn5Q4mNPUvwd27tx5xx13+Pr6Llu2TOJgTlb/flCT3W6vuRNqzg7JHgfSqPk4UIdTc9DS514b6ezZs88++2zN/1+8eNHb21vuPM5XXV29cOHCIUOG5OTkyJ6lVWg9B6fTREVFnTp1ym63nzp1qkePHrLHkUzBA8DO44DdbpeRg9byGmFISMjf/va3vXv32u32Tz75pH///rIncraDBw+mpaVt3rw5JCSksrKysrJS9kRwtnHjxq1cudJut69cuXLChAmyx4EEPA4IGTlwa+kbaCR3d/e0tLRnnnnm7Nmz0dHRK1eulD2Rs2VkZOTk5Pj7+/+8Ylf4BRI1LV68eNq0aV27dh0wYMCaNWtkjwMJeBwQMnKg9MvRAAC0llOjAABIQQgBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCU9n+NZbV3xsTQdAAAAABJRU5ErkJggg==" }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Load the package\n", "using GaussianProcesses, Random, Distributions\n", "\n", "#Simulate the data\n", "Random.seed!(203617)\n", "n = 20\n", "X = collect(range(-3,stop=3,length=n));\n", "f = 2*cos.(2*X);\n", "Y = [rand(Poisson(exp.(f[i]))) for i in 1:n];\n", "\n", "#Plot the data using the Plots.jl package with the GR backend\n", "using Plots\n", "gr()\n", "scatter(X,Y,leg=false, fmt=:png)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "GP Monte Carlo object:\n", " Dim = 1\n", " Number of observations = 20\n", " Mean function:\n", " Type: MeanZero, Params: Float64[]\n", " Kernel:\n", " Type: Mat32Iso, Params: [0.0, 0.0]\n", " Likelihood:\n", " Type: PoisLik, Params: Any[]\n", " Input observations = \n", "[-3.0 -2.68421 … 2.68421 3.0]\n", " Output observations = [3, 3, 1, 0, 0, 0, 0, 0, 3, 4, 7, 3, 1, 0, 0, 1, 0, 3, 4, 4]\n", " Log-posterior = -65.397" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#GP set-up\n", "k = Matern(3/2,0.0,0.0) # Matern 3/2 kernel\n", "l = PoisLik() # Poisson likelihood\n", "gp = GP(X, vec(Y), MeanZero(), k, l)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of iterations = 10000, Thinning = 1, Burn-in = 1 \n", "Step size = 0.100000, Average number of leapfrog steps = 10.029000 \n", "Number of function calls: 100291\n", "Acceptance rate: 0.801400 \n", " 21.801519 seconds (26.67 M allocations: 2.109 GiB, 3.62% gc time)\n" ] } ], "source": [ "set_priors!(gp.kernel,[Normal(-2.0,4.0),Normal(-2.0,4.0)])\n", "@time samples = mcmc(gp; nIter=10000);" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd4AU5fk48Ge2tyu7e5WjHlXK0VGwgVgiFqyJRoPmm2gSNVG/mmhEMRgi0RhbioligcTy+ybGEgEVEFBUkHZ33AF3tOt1++7s7uyU9/fH4nle3dmd3ZndfT5/wezsO+/dzcwz877P+74UIQQQQgihbKWSuwIIIYSQnDAQIoQQymoYCBFCCGU1DIQIIYSyGgZChBBCWQ0DIUIIoayGgRAhhFBWw0CIEEIoq2EgRAghlNUwECKEEMpqKQ2EnZ2dHR0dMe5MCOE4Lqn1UTie5wVBkLsWssETgOO4bJ4BURAEnuflroWcWJaVuwpySuUNMKWB8G9/+9sLL7wQ486CILhcrqTWR+F8Pl8kEpG7FrKJRCI+n0/uWsjJ7XZncyQIh8OBQEDuWsjJ4XDIXQU50TQdCoVScyxsGkUIIZTVMBAihBDKahgIEUIIZTUMhAghhLIaBkKEEEJZDQMhQgihrIaBECGEUFbDQIgQQiirYSBECCGU1TAQIpQe2oPw5nGh0kmCXPbOu4ZQMmjkrgBCKCYnfEKdl9R5iYqiRprJpDyYalXZ9HJXC6H0h4EQofTQ8PW8mwIhTQFoCsDWVr7QANNsqmlWqtAga+UQSmcYCBFKAzyBpoEmoO4Ow442YUcbFBqpy0ZRY3OolFcNobSHfYQIpYFWmkT4oboGu0Okzpuy6iCUUTAQIpQGTvmH36c9iEk0CMUDAyFCaeCUf/gg1xGEbF7IF6G4YSBESOk4AVro4SNcmCeeCPYRIiQaBkKElK4pQDghpj2xdRShOGAgREjpGgbKFx1QRyiZ9UAoQ+HwCYSU7pQv1ve89iABwNbRQb3xxht/+MMf5K5FrHieV6vVctdCBpMnT37rrbdSeUQMhAgpWkSAtmCsO3fEvGd26urqmjp16v333y93RdCg6uvrH3/88RQfFAMhQorW6Cd8zLmgfpYEWLBok1qj9FZYWDh79my5a4EGpVLJ0GGHfYQIKVrsHYRRHSHMl0FIHAyECClaLCMIe8PEUYTEwkCIkHKFedHdfthNiJBYGAgRUq7GAAgiJ4tpx0CIkEgYCBFSLrHtogDgjkB4yOm5EUJ9YCBESLkaA6IHBRJCOnFYPUJiYCBESKGCPNUVjueLmC+DkCgYCBFSqOagKr7VJDBfBiFREgqEPM9PmTKl95b33ntv+vTp+fn55513Xn19fWJ1QyirNdJxTpaGb4SoB0VRff6B+os/ED733HOLFi2qq6vr2dLU1HTzzTe/9NJL7e3tV1555Q9/+EMpaohQlmoKxnl5dochxtUqUKa68MILo/944IEH5K1JWog/EFZUVDzyyCO9t5w8efKGG25YuHCh0Wi85ZZbesdIhJAofhbckTgvT4FAVxhfCrPatm3bov/4/e9/L29N0kL8c40uWbKkz5bFixcvXrwYAHieX7Vq1fe+973+31q3bt3mzZv7b3/xxRdHjhzZewvP8263W5Z55xTC4/EYjUa9Xi93ReQRiURomhaELH21qfWqQyFWEIT4LoEjrZzOykteq1QKhUIsy7IsK2GZNE1LWJokvF7v3XffvXHjRqvV+stf/vL222+PdgxTFNXTQ9zz7/fff3/VqlUNDQ1Go/G+++6Lzh5OUdSbb765du3alpaWRx555J577rnqqqsAYNasWZWVlb3LAYCOjo6f//znO3futFgs55xzzhNPPFFaWgoAr7zyysqVK1mWXbVq1S9+8YvU/x5643ne4XD4/X61Wh0KxZoDnZubq9Pp4jui9JNub9269Ve/+tXFF1+8Zs2a/p9efPHFAzaZjh8/3mAw9N7C8zzHcXl5eZLXMF0IgmA0Gvv8WrIHwzBqtTprTwCnD3Q62mAwxBcIA2pI99+cVqtlWVbaE2DAq+mm7TzNSXiQoby+RG3+9k33vvvuY1n2xIkTarX6zjvvHPrrjzzyyA9+8IN77723urr6rLPO6llGo6mpqbKycvv27Zdffvk999zz7rvvUhRVWVnZv4Tbbrvt5ptv3rBhA8Mwzz333G233fbBBx9Eq/Hpp5/q9fo777xT9kCoUqny8vIoilKr1WazOcZvJbJklZSBkBDy0EMPff7552+99dakSZMG3Gf06NHnnXdeLKWpVCqtVqvVZu9E+tqvyV0ReQiCkM0/fhPNq9VqlUoVXyDsjlBabXq3pkT/9NKeAAPeK78/QcUKKWpJ1vX7m7z//vvV1dW5ubkAsHbt2vXr1w/x9YMHD+7du3f9+vU7d+6MRCI923/2s59RFHXBBRcM+/60ffv2aOSLKiwsjP7jvPPOe+ihh1asWDFgi12KURQVvfbVanVq7gBSBsIvvvjinXfe2b17t0ajCQQCAGCxWCQsH6Es4WbAExl+tyF0hohAiAoTBWNw2ShKxtWMe7f9Dhin/X5/z7+/+93v6nS6G2+8ce3atRs2bOjZnpOTE+PhrFZrdXV1eXk5ANA07XQ6o9vffffdjz/++LXXXnvxxRe3bNkSxw+S1qR8ZtyxY0ddXZ3Vas35moSFI5Q94phZrQ9WAGdcg/FRil100UX333+/3++nabp3hqder9++fTsh5IUXXujZuGXLlpUrV15++eUffvghAHDcUE26A3avXnPNNWvXrg0Gg11dXcuXL1+7dm10+7hx48aNG7dq1ar9+/dL84OllUQDYe9u2JUrV5JvS7BwhLJTQ0CCa6cDJ1pLB3/+85/D4fDo0aMrKirmzJnTs33NmjXXXnttRUVFcXFxz8bHH3988eLFM2bMcDqdl1xyyY9+9KPBil22bNn48eP7b1+zZg3P8+PGjZs2bdrYsWOfeuqp6PZf/epXCxcuXLx48R/+8Afpfri0QaUyXK1evVoQhNWrV8eyM8/zTqezqKgo2bVSLLfbneXJMjRN22w2uSsig6cPCb4ICQaDcSfLAMCiYurikWncTRgMBiVPlnn22WcbGhqeffZZCcuUVp8kzyxUVVW1YsWKqqoqn88nKlkmEWl8nSCUkRxh8EUkuBXiekwIxQgDIULKkngHYVRHCLL83SId4Z9MFhgIEVIWqQJhiCPexFJPEcoSGAgRUhBCSGNAstIwXwahWGAgREhBOkNAs5I1juEyFAjFAgMhQgrSIN3rIOAbIeoFl2EaAgZChBTklE/Kdzh8I8x4PcstDQvXYxoCBkKElEIgpFGKofQ9fBEISNfQihSoZ7mlYeF6TEPAQIiQUrQHISz10kmd2DqqeBRFPfnkk0VFReeee25DQwMAuFyuFStWlJaWjhgx4pZbbnG5XNE9X3nlldLS0oKCgueffx4AepZb6ujouP7664uKisrLy1esWNHe3t5T8j/+8Y+SkhL4uml0sJJ779m7Yg899FBpaeljjz32m9/8ZtKkSXl5eY8//jgADHbE999/f9asWfn5+aWlpT3T1lAU9dZbb82cOdNutytzNgPpl2FCCMVH2g7CqPYgGZ+LnUND8f73FcIyqTlW3hU/orQDrJkXCATa2tpWrVp17733vvPOO/fcc49Opzt58iQA/OIXv7jvvvteffVV6LdYUs9yS1dcccWAiysBwJ49e3q/NQ5Wcv89o2bMmPHxxx9XVFQ888wzNTU1O3fuvPLKKx966KHBlnOKcaEoaX+ricMp1pQLp1jLtinW3johHPV8cz0mOMVa1HQbdd24tGz4SdkUa/QXmwgn5fK/QzCffRml7vv6QVHUiRMnysvLHQ7H5MmTnU5nQUHB4cOHo3e/zs7OioqKzs5OAFi+fLkgCCtWrLj66qs1Gg18PSWbxWLpveZwYWFhV1dX9NOurq7oWkvRPQcrufeevSvGMIxWq1WpVCzLajQaQohKpRriiIIg7N27t7a2dufOnRs2bOhZZNjn80WXYRh2DjlZpljDN0KElKKVlv6pFPNlhmVetEzuKpymUql4/nTjeE+SJ0VRPRsHWyxpsMWVoNeKgz0GLHnAPQGgZ833nrg79BETXyhKFmn5qIhQ5gmw4E/Ca4krDAyPsVDpXn31VY7jnnnmmei65ZdeeunKlSvD4XAoFFq5cuWyZadD9YCLJbEsO9jiSv0NVrJYgx1R1EJRyoGBECFFaE3OqxvBfJl0wLJsSUnJJ5988txzzwHAs88+GwqFxo4dW15eHolEetpy+y+WFF1uabDFlfobrGSxBjti7AtFKQr2ESoX9hFmVR/h9jays13ovUWSPkIAWDZataAw/fJlsmcZJlx6qTdchgmh7NWWtM689iR0PSKUSTAQIqQIbUkLV+0hDISK9uabb8pdhWyHgRAh+XkYQictq6A7BJyAsVC5brjhBrmrkO0wECIkv7ZkribPE+gOJ7F8hNIdBkKE5Je8DsIoXIYCoSFgIERIfm308PskAvNlEBoCBkKEZEYISXY+C74RZhhcXFBaGAgRkpmLgVCS59/oCOI4NeWKfU3B2PdEomAgREhmSc2UiYoI4AxjKFSo2NcUjH1PJApOuo2QzJKdKRPVEYKCLJ2kaBjv1G+M8ClafeKaSZdp1dreW3rWFKysrKQoasOGDb/85S87Ojp6TzcT/XfvPQHgT3/600svvdTa2vrQQw/dd999qal/psJAiJDMkp0pE9VOk+lW7FgagDPkZvhIao5FoO9DT8+agtH/Drgo4IB7BoPB6urqAwcOnH322RgIE4SBECE5kVRlsmC+zGB+PPNmuavwjUcffXTA5ZD6u+OOOwBgzpw54TCOEk0U9hEiJCdHqpZJwjH1aaF/FPT7/QPuqfAV/tILBkKE5JSMxXgH5GchggsTKhXL9u2k1Ov127dvJ4S88MILQ++JEoeBECE5tSc/ZTSKEOKOYB+hEkXXFOyzcc2aNddee21FRUVxcfHQe6LEYSBESE6pSRmNwhEUyrRx48ampiYA6L0q4f333+9yuQ4dOnTLLbf0bB9wT1zLMHEYCBGSjUCgI1VvhADgYlJ3LITSCAZChGTTFQI2hQsk4RshQgPCQIiQbFLZLgoArhQNlkMozeA4QoRkk4LJ1Xpz4ggK6ezatWvDhg21tbWCIEycOPGGG25YtmyZ3JVCccI3QoRkk+I3wgBLwjiCImEMw9x6660/vObHuSeKbxlz5/+U3z2qc/JDtz+ybNkyj8cjd+1QPPCNECF58AS6Uj7bi5uhSk2pPmiGufXWW9mjqpdv+j8VdfpF4ozi6ZdOW/7mvleXL1++bds2jSb++2rvKUZRyuAbIULy6AgSLoWZMlGYL5OgzZs3H/ms/u4lD/ZEwR43zvuhoStn3bp1yTs6LsOUJBgIEZJHijsIo3AERYLWrVt3/ZwfUDDw1ATXz/nByy+/nLyj4zJMSYKBECF5tKVqcrXenAy+ESZk//79k4umDvbppKIzqquqeZ6PvUCXy3XTTTfZ7fYJEyY8//zzPdvff//9WbNm5efnl5aWPvXUU9BrwaYBP0WJwECIkDxSnCkThYmjCYpEIjq1brBPVZSKAhXHcbEXePfddwPAqVOnqqqqDhw40LP9kUceufnmm51O56ZNm1auXAkA7777LgBEl2Hq/ylKBCbLICQDVpBnOQhsGk1QeXl5i6fJZi4Y8NMuf0dRSaFer4+9wE2bNtXW1ubm5gLA73//+/Xr10e3Hzx4cO/evevXr9+5c2ck0ncE6NCfIrHwjRAhGbQHScoTZQAAghwJiXhdQX1de+21W45uHOzTLUc3XXvttaIKVKm+uQlT1Dddj9/97nefe+65wsLCtWvX9v/W0J8isRIKhDzPT5kypfcWt9t9xRVX2Gy2K6+80u12J1Y3hDKWLO2iUS7sJkzAT3/608Ohyr2NX/b/6Hh33cYTbz/wwAOiCly2bNn999/v9/tpmv71r3/ds33Lli0rV668/PLLP/zwQwDoaW6NLsM02KcoPvEHwueee27RokV1dXW9Nz7xxBNjxoxpb28fPXr0k08+mXD1EMpMbbRsh8ZuwkQYjcb33n/vT/vXvvLlX91BV3QjzQT+dfCfj2y5943/e72kpERUgc8884wgCGPHjq2oqDj//PN7tj/++OOLFy+eMWOG0+m85JJLfvSjH0GvZZgG/BTFLf7Bm9u3b6dp+oorruhdwuTJk997770pU6YcPXp0+fLlfcLk6tWrBUFYvXp1LOXzPO90OouKiuKrXgZwu91Go9FgMMhdEXkwDEPTtM1mk7siSfHnWt4xXEAKBoMGg6F305kkFo9QLS5Ng4UJg8Egy7J5eXkSlvnss882NDQ8++yzCZbjcDh++9vfvvH6G6qIRkWpQkBffc3Vq1evHj16tCT1zGZVVVUrVqyoqqry+XxqtdpsNqfgoPEnyyxZsqT/xtbW1jFjxgBA9L2w/w7/+c9/6uvr+29fvXp17/UnAYDneZ/PJ6rbOcP4fL5IJMIwWZrewDBMMBhUq9VyV0R6DE9aPVqgholGDMNQFCV5IGx2CV6TiPx+uYRCIclXYw+HpXkdLigoeO65555++un29naO48rKyrRarSQlIwAQBMHr9fr9frVaHXurr8lkivuvIHHWKCEk2t9LCBlwMM24ceOWLl3af3teXl6fn0GlUmm12mw+vbRfk7si8hAEIVN//DYGVOrhw5vqa9Ie3cuptNo0yJKLRkFpTwBpf5lqtXrkyJESFoiiKIqKXvvRKBDjtxL540ocCEeMGNHc3Dxx4sTW1taysrL+O8yePfvHP/5xLEXxPB8KhUym7J0YkWGYbG4aVavVhJCMPAFcPkGrHb5LoudeIO3RA4QymdIgEAIAy7LSngA63aBDAJFyUBRlMpk4jlOr1am5A0h8PVxxxRWvvPIKIeSVV15Zvny5tIUjlBnaUz7Xdm9hntAsJo4i9A2JA+GqVauqq6tHjRpVW1v78MMPS1s4QplBxpTRKFyhF6HeEm0a7ZN0mp+fv3HjoKNNEUIhDtxyxyFnGEalIhcPofSQHl0FCGWMtiCRfcE5nGgNod4wECKUUrKsvtQHrkqIUG846TZCKSXj5Go9snkxps7Ozr1798pdCzSoY8eOpf6gGAgRSinZM2UAwMVQPUN+s8pVV1319ttv33HHHXJXJCYcx2k02XiLnj59eoqPmI2/ZYTkQrPEG5H/bSzCkwBH5WTgXAXDGDt27GeffSZ3LWLV3t5eWloqdy2yAvYRIpQ6rQroIIzCfBmEemAgRCh1lJApE4VrUCDUAwMhQqnTRsvfLhqFb4QI9cBAiFDqtCvnjTCLE0cR6gMDIUIp4mGIXzGTfOJQQoR6YCBEKEUaA3LXoBc3Q2EkRCgKAyFCKdKkmA5CAGAF4lPAQA6ElAADIUIp0uhXVuDBfBmEojAQIpQKNEucjLJmcsFuQoSiMBAilApNdN81y2TnxDdCJKtOWVeo7g2nWEMoFZTWLgrYNIpktatD2NpKJuZRl49W5elkrgy+ESKUCk0KmGu7D2waRbIQCPy3UdjaSgDgmJf8uZbf1SEIsraXYCBEKOkiAnQoYPWlPtwMyHv3QVkoIsCbJ4T9jm9OPFaAra3k1TrSLd+0fxgIEUq6pgARlBdxeALeiNyVQNnEz5JX6/hj3gEuhmaa/O2wsLWV8HJcKRgIEUo6BXYQRmE3IUqZ7jC8XEeGmGWQJ2RXh/DiEaE15SNuMRAilHQK7CCMwsRRlBrNAfJqHe+JYYbbzhBZd1T4b6MQEVJQr9MwECKUXJxAUv+EGyPMl0EpUOsm64+RIBfr/gRgv4NsOKlOZqW+BYdPIJRcrTRwKXy2FQWbRlGy7e4SPmqOJynLz1IAKXpQw0CIUHIptl0U8I0QJRMh5MMWsqcrDc4xbBpFKLkUmykDAJ4I8DiCAiVHtQvSIgoCBkKEkkogpEXBb4QCAQ+2jqIkYHiIDplPCxgIEUqizhAVlmVgVMywmxAlwydtgnKWoR4WBkKEkqgxoPR7AY6gQJLrCpF93XJXQgwMhAglkZI7CKMwXwZJixCyqZmkV98zBkKEkoUQ0qT4N0JsGkXSqnZBg+Kf//rAQIhQsjgZio55ELFc8I0QSYjhydZWpQ6bHRwGQoSSRfntogDgjSh3vD9KO9vbiJ+VuxLiYSBEKFmalDqzWm8EwB3DDJAIDasrRPamVY5MDwyECCVLY0DuGsQGE0eRJDanW45MDwyECCWFLwKxzLWvBJgvgxJX7SSn0qEvYEAYCBFKijRKnMN8GZQghidb0jBHpgcGQoSSQslzbfeBb4QoQTva0zJHpgcGQoSSotGfNg/IzjRpwkXK1B2Gr7rkrkRiMBAiJL0gBw6GkrsWsfKzVETZE6IiJdvUJKRpjkwPDIQISa8pQEj63BoIIe6I3JVA6emQK41zZHpgIERIesqfa7sPZ1juGqA0xPBkS/qstTQEDIQISa8pTUYQ9sB8GRSHLzuJL4KBECHUT0SA9mCa3R1wBAUSSyDkoDNDThsMhAhJrDlAhHS7P7iwjxCJdNwH3kw5bTAQIiSxtOsgBHwjROIdcGTOOSNxINy5c+esWbNycnJmzZr16aefSls4Qmkh7ToIASDAAoMjKFDMAiwc88pdCelIHAhvvvnmlStXulyuhx566Oabb5a2cISUjyfQmg6LTvTnSp+Bj0h2B51pP3awN420xeXm5nq93kAg4Pf7LRZL/x22b9/O83z/7bfffrvNZuu9hed5mqYDgTR8upYITdM8z3Oc4pd2TQ6GYYLBoE6nk7si4jTTQIfVkhQViURUKpVKlaL+i0aXkKOkvs1gMMhxnFotzS8zHSn2BkgI2d2qjkSS++REcRGajsQ+HtdgMGg0cUY0iQPh+vXr58+ff9tttwHA3r17++8gCMKAd3ZC+g5AJl+TtoZpJMt/A2n64zfR6fpe1RmCaXnK+m2n4wkgIcX++I005U5yFITU3gEkDoQPPPDAr371q3vuueeZZ5558MEHt27d2meHpUuXrl69OpaieJ6PRCI5OTnS1jCNcBxnNBoNBoPcFZGHTqdTqVRpdwI4OgWdTppLl+O46C9BktKG5QMqJ0dB2XNqtZpl2bQ7ASQUCASU+ePXOSQ7yYdAqYjFojObzck+EEgeCPfs2fP666+XlJQ88MADY8eOlbZwhBSOADSnYcpoVEe6jX1Esghy5Ign004ViR8AKyoqXn755UAgsGHDhpkzZ0pbOEIK1xEk4QF6wNMDzUGAzbQbHJJctQu4tFlYJVYSB8JXXnll06ZNpaWl//73v9etWydt4QgpXFPavg5GdYbkrgFSvAOOjAuDkjeNTpky5fPPP5e2TITSRbMSU/xE6AiS8bnpmuyDUqA5QLoy8WlJQX3jCKW7lvQcQdgD3wjR0DJpNpneMBAiJA2aA0+aT72IgRANISLAYY/clUgODIQISSN980V7dIcJp6Qx9UhRDrlIps7Dh4EQIWm00nLXIGECAUcY+wjRwDK1XRQwECIklXTvIIzqCGXCT4Ek1xlK10l0Y4GBECEJCIS0ZcSA9E4MhGgg+7szcNREDwyECEmgO0wxaTuUvjfMl0H9sQI55Ja7EsmEgRAhCWRGuygAdATlrgFSnsNuCHEZcoYPCAMhQhJoSf+U0aggR/ys3JVACpPBaTJRGAgRkoDkb4QM5/Axh2mmCSDVfTM4+zbqzRmGpvTPiB6axFOsIZSFGB4cYclKc9GVJx3/DEecBu0I3u3lhOAo6xWjbFdSVIqWqO0MkYl5OIgCnXbAqdBlESWEgRChRLXQ0twnCCEnutd3+3dPLP4fu3l+KBQyGAxhtvVY12vdgS9nlP1ap7FKcZxhYL4M6iEQqHJmcr5oFDaNIpQoSdpFCZCjHX/yh0/MH/vHAssCijr9TmbSj6oY9bDdPO9A08MRLhWpe9g0inoc9ZBAFvQZYyBEKFGSzCnT0P1WKNJeMfJhjbrvktwUUGMLvluSd351y+8EkvT5TJ1MBi44h+KT8WkyURgIEUoIISTxN0JnYH+7b8eMkb9Wq/SD7TPW/l2Trqy+46UEjzUsgUC3dF2eKH0xPDnll7sSKYGBEKGEuCNUkEuoBJb31XX8dVrp3Vp17tB7Ti75mSdU6wzsTeh4McDWUQQADQHgMz1NJgoDIUIJSfx18HjXa0U55+SZpg67p1plOKPkrrqOF3nhdEIL7eo88sk71Zte7z55OMFq9Ib5MggATvqyIgoCBkKEEpTgUHpfqM5NHxpXeGOM++eZplotMxsc/wKAL19/5m/XTG//+wv0G//37x9d/O+Hvs+GpBnwhVNvIwDIknZRwOETCCWoJYHQQ4Ac63p1fNFNapUh9m+NL7jpq1P30NXk4EtPvfq9t0pyRwAAwzGPbX5g81P3XvnIi/FX6Gs49TYKsCR7uorxjRCh+HFCQjHDGdgrCJGinPNFfUunsZZZl+1++48/XnRXNAoCgF6jv3/pqprNb3IRCe5eIQ58Sc9ORYp20g8ZP46+BwZChOLXFiRxL9lNgJzqfqu88MaeIYOxG2W7wt/RNr5gUu+NVpMtV2sJODrirNC3YetolsueDkLAQIhQIhLJlHEG9gOAzTIvju9qVKYce1mbt6X3xhAb9DJeU35B3FXqDfNlslxDQO4apBAGQoTil8hQ+ibnf8bYr6Ugzlk95yz7yYavXgyz38Srlz7/04RzLtWZLPHXqRccQZHNnGHwMFl0AmCyDELxi/uN0Bc6FuHcRbkL4z703Kt+Wl/19k0bll80+VKT1ry7YZczX33jr9+Nu8A+MF8mm530Z9fcQhgIEYqTnyXeeDNKWlzvl1mXJdIkQ6lUVz/8+o7tt3W0l7Oh8LRrH5l07uWUSrI2HhdDsQJosc0oK2XPwIko5Z7mAiG/O/gszeKC2Uih4h44EeHcLrpqRP7SBCtg1JaMnLpw/FXzlvz0N5PPv1LCKAgAAiFd+FKYlQhAg1/+P/3RrpdOehtTcyzlBkIVRQU4y1P7tmXPWBaUXuIeSt/m3VKYu0itMiVehzLrpa3uzYmXMyDMl8lO7TRJcNbAxDFsVzf91QhzSWoOp9xACACTbIbL3D4AACAASURBVEv2tH7011p+Q71Q6yaC/M8oCH2jNa7WCkJIm2dbmfUSSepgNc3i+KA/dEKS0vrAQJidTirgdbDNu63Yco5BM+gc9NJSdCAsMU8mhPeGjp30k3+dFJ45xO9oJ7TcjyoIAYBASFtcmTLuYJVOk2PRj5OkGhRFleYtbfdulaS0PjBxNDspoINQaPd8MiIv0b6D2Ck6EFJAleZf0ObdEv2vn4UdbcLT1fy/TgpKeGZB2awzBJG4EuvavdtG5F0oYU1K85Z0+T7nBUbCMqM6QyR75hZBUTyBZmkmrI2fi67Ua60W3ciUHVHRgRAAinMXO3y7ey9GyhOodZMN9cKLR4QwL2PVUFZrjet1kOMDLvpgUc45EtZEp7XnGCc6AnskLDMqzEPcabEoTTUHSCTu2ZIk0u7dXpJ3QSqPqPRAqNfYckwTugNf9f+oLUg+bM6uwS5IOeJLGe3yf24zz9KopRnz3qM0b0mHd4e0ZUZ1Yqpaljkpd7sox9Mu+mCxpA+Lw1J6IASA4pzzOwe5yCud5KgHm26QDOIbSt/h3VGcu1jqukBBzgJ/6HiEc0teMnYTZhvZpxjt8n9hM82U/GFxaGkQCItyzvSGjrK8d8BPNzaREKbPoNQK8+AU/6rEsF3BSLvdPFvy+qgonT1nXqdvl+QlY+JoVmF40ib3o0+n79Pi3HNTfNA0CIQqlcFuntvt/3LAT/0s+bAFG0hRSrUE4skh6fTvKspZSFFq6SsEUJx7XpfvM8mLxUCYVRoCIO8oNYZz0eFGu3luio+bBoEQAIpyz+30Dvq0W+UkR9zYgINSJ7520U7fZ0VJe9S1mStCbHco0i5tsa4wiS85FqUjBbSL7irIWUCptCk+bnoEQpt5ViDSxHCuwXbY1IwNpCh14hhKH4y0Rnh/nvGMJFQnSlWUu7A7MHDDSdwIAE60lj3kD4S+L4pyz079cdMjEKooTaFl/mCtowDgZ8lmzCBFKUEIiWP1pS7/50WWhXGswRu7opxFXb7PJS8WW0ezRIAlDiaJ5+ewwmx3KNJhNVek/tDpEQgBoDBnUZf/iyF2qHaRI5hBipLPxVBBTvSZ1u37oih3UTLq0yPPOI3h3JK3jmLiaJY46Qd550/o9n9RmLOAgqR0og8tbQKhzVRBM01DJ4h/0ERoFi9alFxxdBCGIu0R3pdrmJKM+vSgKKrQcmZ3YLe0xeIbYZaQvV2027+7MOcsWQ6dNoGQUmnt5rndQ06fQbPkwxYMhCi54giE3YEvC3POSmq7aFRh7sIhehDi0xmS+UUBpUZDQM6jRzh3kGnNN8nQLgppFAgBoDDnLId/mHmkDrnIYcwgRckUTyBM1aNuvnFaMNIRYZ0SlsnwxBORs+sIpYAzDB5GzjunI7DHnjNHRcmzVnw6BUK7eZY3VM/xwyQqbGzGFSpQsrAC6QqJiwoM5wpGOvKN05JUpd4oSl1gmTd0w0kcOjFxNNOd9MucbNjt31NgkaddFNIrEKpUBqtpmpPeP/RuNGaQoqRpo4EX2U7oCHxlN89N0jj6/gotCwacmzcRmC+T8eRdeokXgt5Qvd08S64KSBwIOY674447CgsLzz777NbWVmkLB4BCy5mOGC7yGhepxQZSlARxtIs6Al8V5ixIRmUGZDPP8oWO8UJcqwYPAqfezmwEoEHWhe0c9AGraapKZZCrAhIHwmeffdbn8zU2Ni5atOjRRx+VtnAAsOfMc9GVRGCH3XNTk4AZpEhyYofSc0LQF6yzJWF+0cGoVIZ80zQHfUDCMjuljKpIcdppEpS1O8np22O3pO5hsT+JeybfeOONV1991WQyrVq1qr6+vv8OVVVVr732Wv/tl112WU5OTu8tPM+zLMtx3/r7UGA2acscgWqraebQNfFysL2ZXDhC9I+gHAzDpCDPULEYhmEYJhxW1stIg5fixNwyHIF9uYYpRNBwgug7DcdxHMepVKKfVm2mud2+r+xGyXpcunnKS/N6dUrPxnA4zLKsXq9P5UEVJWXnf51L3FktLUJ4F1011n5Ln7s9sGw4zKnVsfYpaLXa2HfuQ+JA2NjY+Oabby5ZsqS8vPzVV1/tv8OxY8c2btzYf/tZZ52l0XyrMjzPRyKR/rePfOMcR+CrHN3wqQdfdZJ5uaxehtGZ0giFQgDZm7nOMEwoFDIYZGst6c8bAVdQJ+or3f59+cbZfa/w2MQdCHP1M085/smyjIQdk81etsyY0lMxGgh1OnG/8EwSDoejN4Fkq3NqOE62fBFvqFavLVaBpe9lwnHhcDj281+tVislEPp8PkJIbW3tX/7yl9tuu2337r5je6+77rrVq1fHUhTP82azz8D2vQ+WqBZVN//WYPhZLIU0CuYzC9L4pcpoNCoqEqQSwzA6nc5qtcpdkW80O4nBICoPS/CGqyaV3KzXxvNHFATBYDDEEQgNUGLUlTDkVL5xehzHHVBIa7RaU3qv1Ov1LMvm5eWl8qCKEg6HU3D+8wRcjYLBINsDd4uvqjB3fv8bHaWh8vMNZrM5BXWQ+MwuLCy85557SktL77rrrpqaGmkLjzLrRlGgppmmWHbe3UWEbH2jQpJrFjnFqDdUr9fa9Nqi5FRnKAWWuY7APgkL7BQ5aASli6YAifDyjiDcX2CZL2MFQPJAeMkll7z22msMw7z44ovz5s2TtvAeNstcR2BvLHu6GVLnSVItUNZpCoi7XzgD++xmea5wu2W+MzDMQCNR2uJaeQopn7wDJ4KRVkGImPVj5ayE5IFw7dq1n3zySXFx8bZt29atWydt4T0KzHOdMSfF7e7GCxhJIMQRsQsSOQP7CiypXmI0ymIo5wQ6zHZIVWBHCMK8VIUhBZF3ilEnvd9umUuBzO0NEgfCkpKSLVu2eDyenTt3TpgwQdrCe+SbpwXCDawQ05NMo5/Et4wqQr010yDqNGI4B8O5c40Tk1WhIVFA2c1zHdK9FAqEiH0hRsrH8KRN1tkSnP79dsscGSsQlU4zy/RQUbp80zR3oCrG/fd04QWMEiW+XfSA3TJbxkvMbpnjDEg5mrBB1jY0lAwNARDkuzsKQtgfPm6VaaLt3tIyEAJAgWVO7K2jtW7ijSS1OijzNYmcm99FH0jlOPr+rOaZvtBRXmCkKlDeyUdQMhz3ynl0V7A6xzhRLd+EMj3SNRDazHNc9IEYx9gJBL7qwtlHUfw4QVwLkkA4F33IapJt7kQA0KhMFv04b6hWqgKxmzDzHPfJeWN0Bg4ooV0U0jcQGrRFGnVugDkZ4/77HSSCoRDFqy0InJjzxxs8YtaP1GlkHgNnz5njkK51VCCkEV8KM4gjDG7J2gvi4aIP2M0YCBNjN88ZdiWKHmEeDjowEqI4NYrsIHQFD9plbReNspnnuiSddLRR1rVbkbSOeeW8JdJMMwBl0o2UsQ490joQznYFDsa+/+4uwMH1KD7NIgOAM7DfZpZn4ERvZv1onkRCkXapCmzAxNEMctwn59GV8zoIaR0I80xn0EzTsOv09nAzpF7WnmGUpgghzWJG4ERYZ4Tz5BrHJ69KMaKAsptmuWgRz4tD6whCiMNYmAkiPJH3/d5FV9os8reaRKVxIFRRujzjGe5gdexf+bITr2EkWlcYQmImzXYGK63mmQq5uGzm2U66UqrSBEKaRM4zh5TplB84+UZO8ALjC9VbTTPkqkAfirhW42azzHKKedptDJBWHFyPRBI7gtBFK6KDMMpqqfAEawUi2So7OIgiM8jbLuoN1VoM5WqVUc5K9JLmgVBkNyEA7MEZ15BIIkcQCi662moeZr3MlNGqcsz6kb7QUakKxGH1mUHmgRP0QZtZzsFFfaR3IDTpyihKFWSaY/9KrQu8EYyFSARRb4S+0AmDxqbX2JJXH7Hs5tlOkc+LQ+gIkiB2E6Y5BQycqJJ3uok+0jsQAoDNPMtFxzrXGgDwhOzFl0IUM18ERE1L5KIP2JQxRriH1TxbwnwZIn6SHaQ08g6cYDgHy/ks+nEy1qGPzAiE4nIB9jsAB9ejGIkeQUhXKqrNBwDyjBMZ1sHyki1IhqMJ053cAycqbeZZFKWgFS7TPhBaTRXe0BEisLF/JcSRShxcj2Ij6u2HE4I005RnnJK06sRHlW+eJqrhZGin/Hj5pLEIT8Q+3knLGThoU0wnelTaB0KN2mzSj/KIzAXY0y1uSR2UtUR1ELrpQ7mmySpKl7z6xEdsD8LQOoOA3YTp65Rf3HyBUhM89CGltZqkfSAEAJt5tjso7iJ3hkmdB69kNIwwD6IW43XRlTZTsh511TyrInHewGzm2S66kkj0+IfdhGlN3nZRX+i4TmvVaaxyVqIfjdwVkIDNPLO+48XywptFfWtPF5mSr6BGaqRAzQFxk/K5g1UjrZcmpSqEXLz90dKOQ2GjlTbaQkZbwFwYNFqDpoKg0UabC4JGK6PLGezbBm2RWmWkw40Ww1hJqtPghyn5kpSEUu24rEvSu4KVisoXjcqEQJhrnBjmuiKcV9Rk/w1+4otAruIasZCCiHrvCbEdgsCY9KOSUZOJJ7fpI8H1N/5Hz/jMIacx5LbQ3aaQu6TzkDnkMtEOU9jF6HL+vfzvAqUesASbeaYrWClZIAwQAHyOTD+OMLgZWQMhXTXWfr2MFRhQJgRCCtRW03R3sKo497zYv0UAat1kYTFezGhQTWLmIXLRlVbzTCoJ4cHA+BYcfO2jC1bzam3QZA+a7APudtnHD45r3HVi7PkDfmo1zWzzfDjadpUkVeoMQZADUybcP7KLvAMnOCEYCJ/KN02VsQ4DyoQ+QgCwmma6xecCHHJhNyEaFCeIm5DPTVdZk9NBuHDv34+VL3XYhpnFu3radRW1/4JBWnNt5hneUJ1AxAyKHBzBtQnTk7wdhJ5gba5hogKzyTIkENrMM50i82UAoC1InOFkVAdlApGL8QruYI3NXCF5NUa2HSjuPnKw4sZh92weMVclCCM6B56GXq0yWfSjvRLOtYZLMqUb2QdOuOlKm0VZ+aJRGRIIjboRKlDTYuZai6px44goNDBRHYS+0HGDxi55LpyGYxZ99dfPzryL1RiG35uiqqdeU3H47cE+t5pni519Ygg4+3baaQjIO3ACXHS10gZORGVIIAQAqxlbR5GURI0gdNGV1iRc4fOqX+8sPKN1RKxztp0Ytzjf3WR3nxzwU5t5lisgWSDsClNByda0QKlwzCvn7S7MdrO836wbK2MdBpM5gdBuniV2NCEAOMLQGUpGdVB6E7sYr4uulHyyjAJPw8RTO76a+6NYdmYEcDJEUGmOTLls+uF3B9wn1zghml8tSfUIkbmdDYl1QtYOQnewymquUNTMaj0yJxBaTTM8wcNxrLt2yIWto6iv7jAV+2K8vBAKMA15Rilz4VREuPDAC3vm3BoyDDxez8OQOi/Z2UH+7xR5vpY8XSO8XEc4AQ5PumxU614L3dX/KxSo841T3aFDUlUSl2RKI44wuOQeOGFP2nQTCcqcQKhR5xj1pb5QvdgvHnIRIm7YNMp8jQERj0eeYE2ucaJapZewAtMOvxPRGI+NXdKzxceSSif5sJWsPyY8US28chz2OwgnwLR8uGkCPDhDNdIMtR7Cak3Hy5dOq/tgwGLjmKR+CNhNmEbkHThBCHHTh6wWDITJZzPF0zrqjUALnYzqoDTWLOaUcNJV0raL5gS6Kg6/vW3uT8nX7UisABuOwXE/5Grh3BLVz6eq/nca9f3xqqUjqGlWyq6nKArmFlD7HQQADk29auLxLbrIANk+VvNMd0CySUe7wkBjN2GakHfgRIA5pVXn6DUFclZicJkVCM2z4lt3rcaND7boW0SljLqD1dKOIFy494WaqVd7zSU9W7a2kjITXDeWWlRElecMPJJ9Uh7lY6mOENCmgpayuVOOfdR/H5OuDESuZT0EHE2YLuQfOBGstCpsxYneMioQ5hqnBJlWVhDdcVHjIgJezuhrvgjxxNybwnAOlvNa9OVSHX3Cqe0WuuvQGVf3bDnlh6NeuHTkMFkGFMAsGxx0EgConnb9tKPvqYQB3tdspgpXcOCxhnHA1tG0oICBE5V2DISpoaI0ucYpHrpG7BdpDq9n9A1RC8+6A1U2k2S5cPqI/8z9L39+5s8F1emXPkaA95rIFaMpQwzzmc0pgEMuwgrgyh/jyR8z4dSO/vvYzDMlXJIJh9WnBXkHTvAC4wsdzzdNl7EOQ8uoQAgANnOFm47naRdbR1EP0SMIpZss46x9606MPb+z8JulfT9sIRNzYUJuTF/P1VKjLFDrjr4UXjvj8H/6z7hmM8/0BGsJ4SWpcHeYwm5C5ZN34IQ3dNhiGKdWGeWsxJAyLxDOdsaVFHfYTThsHkUAICYQEiDu0CGrSZqZ1Uo7a0o6qw/M+mZBsXofNPjhojIRr5tz7dR+JwGA1pJZvEozsu1Anx006hyTrtQXqpOkzthNqHyOMJF34ISbrk7eOp2SyLRAaNKPEggTYjvEfjHMwwkcFIUAwjzpDscaeGimQaOyGLSFiR9XzbPn7H7+ywV3RDSnH5xDPLWpBa4aQ+nEXKYT8yg/S3WEAAAODTLjms08yyU+v3owDbhIr7LJmy8KAC76oDJnVuuRaYGQAirukVI1ON0aAmgOgBDzuFIJJ5SZfPwjb25ZU9n8ni1bOrXTrTDGIq4cCmCWHQ44CQCcGnNOTqC9wHm8zz5W00wcTZg95O0gZHlPmHPmGIdZO0VemRYIIYFuwqMeEsFJZrKeqIETUk0xqubZmTX/Oljx/Z4tNW7iYKjFxfGUNscONS7CCiCoNIenXNn/pTDPdEaQaeF4aV7lusMQYDEWKpTsAyecgSqraQYFA68XrRAZGAitppnu4CEA0TGNFaBe1kcnpASxL8bLC4wvdMwqRS7clGObHfYJ3faJ0f/6WfJhK1lWymriukBztdQoCxVNmTk68dLSjupcf3vvHaL51e6g6PzqARFCROXZolSSfeBEdIpROWsQgwwMhDqN1aCx+0J9m4NicciJgTCriVqMV6pcODUfqah9++CMG3q2/LeJzLdTJYb4b2Bz7dQ+JwEAVmOon3DxtKPv99kB51rLEvK2ixIgLrrSruwOQsjIQAgAtnjXXTvuIyEOL+ns1R4S8fjsoqskyYU7o35zd8FEx9evgwccxM/CuSUJDUycmAuBr1Nmas64csLJ7QbmWykTNoukgRDfCBWJECJvIAwyzWpKb9CWDL+rrDI1EMZ5kfMEjngkrw5KG6d8Iu4abroy8Vw4NR+Zcfg/B2fcFP2vhyHb2snVY1WqxAboUxTM/jplJmSwNo5aOOXYpt47mHQjCWFDkfZBChDHEYYAK0lJSEqtQfBE5KyA8vNFozIzEOaZzqCZRk4IxvFdzB3NZvUxL9UX4dxhzpVrnJDgEafWb+wqmOy0jQMAQuC9JnJOMVUUw3L0w5prpw65T+d/VU+/dtrRDzT8NzfFaH51HJPUD4gQgv3rClTjkrkCrqDE89EnSWYGQhWlyzVO8sSVC3DKT/z4bJuVaA5i7yB00ZVW04wEryANx8yo/c/BmadfB/c4CE/grCJpZmuzaGGM+XTKjCd3pMM2flzjrt472MyznNLNtVaLczMpDAE47JHzjyKQiC9Yp/xMGcjUQAinW0fjucgJwGE3jqLIRvUeEetSuoNVibf5nFH3QUfRVFf+GADwRMhnHeTqMVIu4N2zMBMAHBt/0YRT23t/ajVXeII1BKSZa+2Un+Bca4rS4Ce+iJyB0Bs6ataPVqtMMtYhRhkcCOPMlwGAGre0dUHpIfbGPQLEGahKcDZ9LReuOPJO5Ywbo//d0gpnFVJWvYRxECbkQICj2oMAAI0jFxQ4jpmDzm8qoM4zaovjWMt6QAKBo7K+f6A+auW+j7mk6ERPjRgmtE9PJv0oQQiH2A6j+ISlFhrcDJH2loQUjhPIyZjn2KOZBq3aotcWiToEIeTIJ/85uWcrE/CVTJl1y9SS9uIZLutYAGihoYWGq8ZIfMpRFMyxwwEnucxE8Wpdw+iF4xo/rem1wFO04STPeIYkh6t1kbkFeNUogkDgSKpatvgIc+Ddl1tq9qjUmtGzzp55+QqVWgMALvrg5OI74iiQi4QPvLOutfrLE0W6JUuWrFixQq1O7nh86d8Ia2pqzGaz5MWKdToXIK6XQkKI7A9TKMUa/MDwsXcQis6F4yPMG3dfXv306gt8I67TzNVu233//XduMs8CAAKwuUW4qAy0SWigmWOnar6eMunEuCUTTn6rddRmnuWi+87KHbeGAC5YrxQnfClqqfZ3tb5ww2z/2+9eKUy7lJnQvuHVdbcsDPlcEc4dZp254mdW83Y0/fW7M4Pvblyuml3WOfHPj/594cKFHk9ys/klvvK8Xu+tt94aDMaTrim5RIYM46pM2UZU0qMzcNAmcumlvf/+m7mx+8Ub37xyxnVLJ1/60MVrbpx/9+t/egwAql1ERcE0a1LepaIpM9Hzua14hiESsHoaez7NNU6hmVaOl2a+eYGQI3jhKEPKcpe2PP/g0oIzH7/yuUvOuGLZtKuevubFGVTZjr8/5qIrbeaKOELMx8/88tKS89dc/vTFUy67bNrVz123LsdrW7NmTTIq30PKplFCyK233vrggw9ef/31g+3T0tLy+eef999eUVFhMHwrZ5zneY7jBCH+t/s844y6zr/zPEtRol+r2wLQ5ucLpchijxvLshqNJtltAorFfi01hzviIoIQUyjihXAgfDJPf4aok7Pu0w9+POcHKuqb+8JVFd974S9Pe7o7tncVXT0aiNB3GTCe5xM5/3vMscHOTphlJQBwfPS55Se375214usPVfnGqY5AZVHO2YkfCACqHcLMfElKSvUJoEBx//icQGqdEOP5nKC6Tz9YfdM7vbfcMGfF/37yy9HfH51vnin2BCaCUP/ZxrW3bu698XtzVzz1/m/Wrl079HfVarVKFeernZSB8Iknnhg/fvx11103xD6bN2+urBzgLe3ll18eOXJk7y08z4dCTFhIJBbp9Opih+9Qjn7K8Pv280Uzv7hIzoYen8/HsizDMDLWQUaRSISm6dQ8B3QzVKdfF+POntABk648EgGAcOyHCPndtvH23lvUKnWe0frZSVdZgb1QzYb7FRaJRCiKivva7lGmBZrVN3rYYoNQU3b28s9/t2vy9eTr7NQc/TSHf3+udm6CR4mqC0OrlbFIcV8JhULZHAUBwO/3m0zxpFwe86u8Qa3k9elP4DkuHMw3WntvtJkLwgGvO1hVlnt9uP9pPSSOCQHH5fUp0FTg8Xi83mEG+ebm5up0sV7FfUgWCLdv3/7hhx9u2bJl6N1uu+221atXx1Igz/M5nT4Tm1DqbWHuPJo7XGydE8d3T7JwlU2lSXCGjwSo1Wqj0djnRTl7MAxjMBhsNlsKjnWkg5hMsT66tviPFObMEXuHKiqfWtNeNWPE7J4t3YFODxcIaMfePkpj0g182zIYDIkHQgCYU0BqA/pxNipsmhTR55YHTrUXn54rvFR71oHG/xpNRgqkOdW71ZaxUqTMBINBlmXz8vISLypNsSxbUFAQxxd3+ASTKUVNo/kjyw93HJpZ9s2DVE17lXXMOJ06Lz9nlOjiTCZzcdnRztpppd+MPqxtr5o6dWp8v4oYSdZHuG3btp07d+p0OoqiAICiqF27dg37rWSzmWc7A3HmAgQ5EfOMoLQmqoPQRR+0W0S/P535vbv+sfflqtb9pwsJOh/b/GDexbfPLzXm6ZL+sDWngKr1nk6ZOT5u8fheAwoN2hKVSk8zTVIdC0fWyyvFq+gs/P7df9j6WIvn9Plz0nH8TzuenLT8PJslntePaIFPbv1Nu7c1+t/j3XV//vSpe+65R5rqDkKyN8I1a9b09GdSFEVEDE1OolzjRIZ1sLxHq46n4+Kgk0xNTgoDUg6aJa10rDuHIu08iZj0oh91y6YvuGLthkd+f5eV1Zl05lO+xolX/zR8waOLJJpHZmgWDYyzUIdcZG4BdWLc4mv/e+fu+T/h1Kfbkezm2S76gEU/RpJjNfqJn4WcVLTMoQHUeVO6ruq8637KRZjb1q0YaSwRiNDOuZbcu4bMrLebZw//5YEsuOHnPMf+8JWbRptKOcL7Kfcf//KHq666Stpq95Gx4wijKFBbzTOcgYMleUvi+PpxL/FGIC/OZmeUHo75RCxJ76QP2E2z42tFnHjOsmt/u9p7svbziVddXn7GG+05S+2ULlW5UPMLqI9ahbkFVNBkd9jGj2rde2r06QQZm3lOk+vd0barhy4hRgTgiIcsKMQnSHmkfrbks75/99xrbus+eUSlVheMmwJq7ssTP8kzxTk4laKoRT+4b/51P+0+VnXrZPXcuXP1er20Fe4vKTPLKOR1MMpmmeOKt3WUAFTiCoWZTly7aOCgzRLno66OpWcc3+y5+OcjK846zuZGeKhIYXvDuBwQCDQFAACOly/p3TpqNU0LhE/yQkiqY9XizPUyYXg4LmYFFaloDaYRU+eWTJ6l0Rncwep801QVldALhNZoLps6d+bMmSmIgpDBU6z1sJlnO4NVcSxYH3XAQWJ/XUBphydwwjf8blECiXhDR+Jeg7Ci9j/NZWd6c8o4Alvb4TsjVVLOKxqDuQXUXgcBgMYxZ4/oOKSPnB4+qFIZco2TXHS1VAdqCsg8y2XWOuIh8q5HDwAO+oAt3nZRuWR+INRrbAaNLb4F6wHAGyGnpBltjJSowU9in1DGQ9daDGM16ngmTrLQXVPqNx2YeSMA7O4iI4wwxhJHMQmZbaeiE45ENMbmEXPLGz7r+chmme2iD0p1IAK4rqc8ZF9FjgBxBw7a482UkUvmB0IAsFvmOen9cX/9ILaOZq46MfNEO+kDNnO8uXB7/35o6tUBcxHNwpddZOmI+IpJiE4FU61UdLXeE+VLei9GYTfPlXCuNcDcUTkEOZD9qZ1mGtQqo/KXpO8jKwKhzTwnkUB44I933AAAIABJREFUxE1oFq/qzHQs5nZRAHDS++2WeXEcZWTr/nxvc3S260/ayWwbZZNpSvczC6l93SAQaBkxJ9fXnhPoiG436coA1BIOomimwSvr2uhZ6LCb8HL34zgCe+3mmAYX6SP+0s6aXH+bikizEFgiMjxrNCrPNDnMdEU4t05jHX7vfngC1S6ysBiz4DJNVwjcTKw3jlCkjScRs3602KOoBO6sfS/unn87r9Z2hqDeS+6aKtsDaKEB8vWk3gdT8tQnx547/tSOyhk3RD8qsMxx0vvi+AEHRAg54iFSLTKMYqGEGZJdgQPjCm4c7FNjyF3aWVPSVVvSWW2hHW7rGFPQZQq5/DklntyR3tyRntwyb+4ob24Zo09pz0FWBEIK1FZLhStwoCR/aXwl7HeQhcXSVgrJr84rIq/AGThQYJ4Tx8CJ6Ufe9eaWNY+YBwAftwqLR1B6WaePnV9A7e0mU/Ko4+VLzt/1x55AaLfMbXS+Pdp2jVQHqnVhIEydAHs6JVhGHO8PMM35xqm9N1rortLOmpLOmuKuWgPj6yya2lE8vX78/7ps5QKlAgA1H8nzteX7WnL9LWUdlVPrPsj3tXJqrStvDEx/NDU1z4pACAB281xnYF/cgdARhhaajDTjVZ1R6sUkdDjofWX5l4o9hCnorKh9+/1LnwaAOi8JsDDHLvNZdEY+9VELcTIE7JMoQgqdx7rtEwEg3zS9tvWPHE/Hlw3UX0sQPAzJx3U9U6LGLcie3+6kD1pN0ymVFgBUAjez5v8mndii4dn24ukdRdNrzljuzhsN/VKlebXOZR0bXZizhynksvhap6eq5tkSCG3mOce6XhYIp6Li/JEPODAQZhSaJa0xLxfGCyF/6Jit7EGxRznzwCtHJy3z5ZQSgK1tcEmZfHPXfk1NwWw7tdcB3ymD4+WLx5/aHg2EKkqXb5rqog8W5Z4jyYEIIUc8FDalpIYSllB1BvZHZx+0eRrP+/ypoNH+0QWrPXnxNLYHjbaQ1gyQotCeFckyAKDT5Jl0Zd7g4bhLqHGJyLNHyidqQhkXXZlrnKxWiZsAvaT7cHHX4arp1wNApZOYNWRCruh6JsO8Qqh2ElaAE+OWjD+1sydbwWaZ6wjEn1bWH+aOpoY3Ai0xTxOYNIKTPlhonj2z9l+Xbvn10YnLPr7gN/FFwdTLlkAIAHbLXCe9L+6vRwRFPHMhqYhciXe/2KFRKiIs/Oqve+b9iNUYOAI7O+DCEUq53HK11BgLdchFfDmlvpzSsvbTIwgLLPNc9MG4Z5/orzUInpjTkVDcDrkE2efz8gSPWlT51215fER71XuXPXd0kuh+BBkp5cpMgQLLfId/byIlHHDgJZ0hRE0oQwhx0gcKRA6cmFK/kdHlnhp9DgB81U1GmGCkNF1v0phfSH3lIABwfNyS8Sd3RDfqNQV6rdUbqpPqKISQw2JGaqL4yP+MTghp/n/ntbSfGLd489LfBkyFcldInCwKhGb9WIFwQaY57hJaaNIp2XSMSE6iJpTxM8e1KouoMcL6iH929f/bPf8nAMDw8EUnWVIaTz2TZ1wOcALVTMPJceeNaturZU+f2QXmBQ5//A0n/R2W/R6d6Rxh0h6U82kjJ9B+2ZYH2yNHvdPvqzljef90mD4IQFcIDjjIh61kVyepcZNWGmhZ12DOlmQZAKCAKsiZ76D3jRa/hk6PSie5ZKTs6Q4oUaJWmnT69xXkLBBV/vwDr50Yt9iVPwYAvugiE/OoQoOyThsKYH4B7O0mo8bmdBZOG9v85bHyCwDAnjPvaPufxsMPpDpQC03cDFhTMXNylpJzWjVCphz7cG7VP3ZNvdCjd1L2MwfbMchBCw0tNGkJkjYaLFpqpBmKjUCzcCQInojgiQArQL4OrHrI11H5OrBpqJQly2RRIASAAsv8Bue/E1lupspJlo4AGZetR5I4JqaD0BHYM7nkjtj3L3AeG926999XvgAANAd7u8lPpijxhJlpp3Z0CjRLHS+/YPLxj6KBMMcwgePpMNsh4SxZh93C2SVZ1PiUYnJlJOlY+sIdv1Px7H+/81Qtu9cemddnlK2PJUe90EpDCw1BjpSZYaSZWlSkKjODse9QWgoAIgJ4GHBHiIcBNwudtErC7uqhZVcgtBqn1zJ/ZHmfVh1n9l6QI0c91HSbtPVCKdUVIq6YMzjCbFeE9+UYJsZaOiGL9r6wb9aKiNYMAJ+2k1l2KgVr0MfBoIap+dQBJ8kdueDs3X82hVxBo40Cym6Z1x34apT1SqkOVOuGs9Ns7sm00REk3WEZjqtlg5dsW+W0Tfxy/k8IRTka944u+OYFgwDsd5Dt7WRyHjXOAucUQ4Fh+LcHnQqKjFBkPL0jxaVu6rXsekyjVFqraaYjkFAXCM7Bne5EtYt2+/cUWOZRMS+YNOnEVgJUfflSAPAw5JCbnK3g2VUWFFD7HMCpdQ2jF5Y37IxuLLAscPi+kvAobUERTx5IlB3tcqw+yIUv2b7anT/2i/k/IRTF8f5ApMFqmhH91BMh/zguVLnIrROpK0dTs+xUoSHFC46JptxAyLKsp7tDwpxggWP93W2FOQsc/j2JlHPSD3hVpzVRAyc6XV+Y2Qkx7qxj6XmVG76c/7NovsAn7XBWEWXWxlPJ1Cg2Qp6O1HvJ8fKlk49vAUIAwEyNczvrWF7MfOTDOYwDCpPglJ8cjTcpN+z3MAExj4Rf03GhS7Y+4s4fs+vMu6LnuSOwz2qaoaJ00RfBl+rI+BzqfyaplNYvPgQlNo22trbed99977/7vlmb4xfCC26467wf/VqjN8ZdoK+z+aOn7z/+2aZcXY6fhEq/M3ryvXfqDHG2jhJCqpzUEjmW0UGJi3bax6KzvnrTUz/vOrTfqH2ZNerOv/2RudfcPvSr4ZyqN5pGnemwTwCArhCc9JPLRyv3WTNqfiG1t5tMnjCdV2nYHRvWr/+7s75Gr9YfzJly0V1PViy7WZKj7O0mCwqJTp02d0blEwhsbo4nCh7btenjZ38Vbm8lQMwjx1x8zx/GL7w4xu9qufBFn6z25o364uunPQDoDnxVlHOmmyHvNRGBwP9MouzpNq+e4gIhTdPnn3/+mdbz37/9U4PG0B3oenLLo+82H7tu7ZvxFRgJ+l/58eLloy95/ief6TX6Ln/H77esenv1jTeu3Rh3JQ86hfNLValeXxxJ4ZiXxDKhjKv5+PrbL7jrrLsvv/OvKkp1wlH/2EsPhryuc//n14N9Jd/XMv7U9rev/Fv0v1vbhPNKKZ3S4yBMzac+biFOhmw0z/zTw7/43wtWLb1oHQXUkY5Djz39YCREz7v2J4kfxRuBne3kIsy4ls7ebqErJDoQ1n36348fue3hSx6fP2YhAOxu2PW7B26+/Il/xhILtWzokk8edeeN/vzMO3uiIC8wbrqa1tzxWSNZVEQtKoq9G0FBFHeZvv7663au5Gfn3mvQGACg0FK0+rKn2j7f5mg4Gl+BlR9smG4cd9uin+s1egAoyin57WVPd3y2291yIu5K+iIihmMjRamLrV109xvPXTXl6itnXKeiVAAwvmDSmsuf2fXak3yEGewrC/a/Uj3turA+FwCaAtDNUHNtaXBLUFMwu4Da54D/9/G278+77cLJy6K5f2eUzHj4O2s/Xfc7qbonvuzCYbiSCXGwsz2eL366bs29S34djYIAcNbYc35+3i8/ffl3w35RwzEX7XjMm1v2Ra8oCACnPJUBYcLxgPn2ydTZxWkZBUGBgbCqqqrnjxRl0pmnlc7srK+Or8DO+up5o8/qvcWst5xRMqPjWFX8tQTY2419HumH4UmMTzAd9VXzR3/rPBxlHWPT5LhbTw64f0lnjc3TcHjyFdH/bm0TlpSAWnGX18Dm2qHaRdrrD/X5kaeWzhC8vqC7S5KjCAQ2Ncs/E1hm+KRNCHKif5OEkM5jh+aPXtR747wxCzuGu7tqOObi7b8JmEt2nfUL0ivWfdpBKjt32y1n3jxepczU6Bgp7ko1GAx0pO+yWgHGr9GLm++4h0ZnCEb6dgrRjJ8hHfEVGFXvjb+bGsmlygUxTiij0RuC7LdOGwIkGKEHPg8JOXP/uq9m38KrtQBQ5yURAWZY0+a+kKejRpkpVqPv8yNzPBfhGY0uzkuvv0Y/qXJJVVj26gqR/XFN90hRlFrb968cjASGvrtqOObi7asD5pLPFvaNgnVerlizf3bxoOPo04XiAuGFF164vf6jMPfN0JhG18lj7uOjZi4a4ltDKD9z6da6zRHumxatk45jJ3wn1WMGbeOK0eZmXI8izezrjnV8bvmCpZtr3yO9Jrb4/MQOTWFhXunY/juPb9gJQJ0ccx4AEIBP2sjSEWnWg7ygkOImL91U+27vjVuObiyaMlNvyZPwQB+3kCAnYXnZ6KMWIsR74yk/c+nGb/+VN9W+V75g0IVaT0dBS1GfKFjtIgedsKy0xqQbodPa46yNYqh/85vfpOxgO3fuJIQsWbJkiH0mTpz45YEv1m9+2WqyMTzz5amdv9/y6MKfPTJ23uL4DmofM6mueufHn663mmwMG/r85M4ntq1eeOevQsWHRtmvoBJ4FGB44ASYkJesG144HNZqtRqN4hKaUoPneZZljcb4s4X7OOUnX3TGev8onjTtg3cePlxflWe0+sO+D4+8/8KXz1/x2Cu2keV99lQJ3IWfPv7FmXcGLEUAUO0i3QwsHSHBWcGyrEajSU2vi1UPlTmzG7Y+dbLpYI4+1xN2v3/oX6/u/9uShx8pKpsl4YFYAcI8TM4f/odiWVYQBINBsvfRtBMIBHJycvpsPOwmuzrif/4eMWX2G6896PK0m3XmTn/7m/te3dS09Zo1/zDk5PffWcMxF+1YHTQXfrbwnt5RsDEA7zeRm8erfIF38k1T84xT4q7PECiBX1QEOp0uGYX3ocSb7FtvvfXGG2+89dZbNfubLaOmXPn8f0bOiP/Vm6Ko7z75r+pN/3x5y78Cde0FY6dc/Zf3Rkydt6/hPm/wSL4poTWQ93SR6TZcsDc9iOrW9QvHzn/itsiXhX/6/K1wwFs6Zfat//zCOnJ8/z2nHXnPnT+2o2gaAPAEdrbD1WMV19ASi0Xj7If+P3vnHR9Hcff/7+z13nRNvXe5y02uuIALxRgwxoRQQ0gICRB4+EGeEAgEkicJvYfQgkN3AWODjXuVu2XJVpdO0vXe7/Z25/fHGSFbsnwnnar1fvmll25vtTt33pnPzHzbs4ciB19+74e/+sOktHzp1Q+8hEV9CTXrnWNWPEGB04RjvSZuIjTe1tGvrGPy9Lxffnpi7wf/95eTbyCCSJ9Qcd8/T3DFsh5Pnln5epAr232+ClqD+PMWvDKTUPKo+o7KKRl/6097hgnDUQgRQmvWrLn55ps/rnY3kQmoZIoQGr/sZ+OXnZdHWCmaafbs76cQYoCvW+l7ixhjyUeHOR4S1zrjON/i2a+RzU5Zs2TGmgd7OY0T8o6r+WrT4ueiLw9bsIoHacOp3FLsjJOjHQbBotWP591053Wbfvfp9X/xoFBl8+/ycYRAiRwoMMA3OnxvERrrNfGy3wT2/pp0gCdRLHrg+UueltF2UGuq/mr5K11V0B+B/zbBomSULQKHr5rLVHFYqv42aBgwIqeuCUEtrrB4DvY/qaspAIfMg5QZdow+c9QKVMz+ijQOW71HVeIZlzxzYtXa5ozZ0TLcIQr2mfAV2pE6uhMIypVw0Iw9Qk1H8sT8hu85TLmAk+bw98u/ukdMAVwZs712jCjuMN5rHKQvjRtyVxx6bVfFgyTrJ9sEScPHjfQ4OYyXIwAwe/YpxX103RhuXL5CyGVpOEy501/T/0tt12PHWNK1YQyN4yuqbPMeF3GzWIwerCZdEXkNuU07jo9bFX25w4gLpUiVMJvmEDBFgWrd2E3ikyU3ldWsJ+iISlRhdu8biHvt0IM7PNZr4mCbHocHa/JQcfC1upyF0Q3/KBjDV61YzkFzNAgAMFAWzyG1aEwIRz4qUYXJvbf/1yFp+LZPuY7GGBzOOHFcY67FvVclmnXJ06Yefb+qeEWAKwMASxBOO/D8EbscjMJlwngZOmwBmzzLKUnLbd6pEs2weg/TOJzwe4Uo/F37WK+JlTYvrhqsyJO8xm1Sl+74uNVdD37XgQMRfF36uX1Sh+8Uj60ZHfuicJkLYXR3FOMEFPuod+GxtMLDlrjcZGg6aPMfV4qm936aylqrsp6tLrw2+vLbNnqeBvGHo809PqYr4ZgNhyk4VbKyrPpLDkMq4mbZvccG4l7VDhxXYcjLFgywpR0PTi4Cgd869dj7O2b/D8X4yV3zoBk3eeDmLKIzR4TZvVcdw2RxpDCshZDrtUrcHSxyoPIycVgqPltr959IyNW+baODg1c/a4xYsQSh9cIMDb2e7zsi5uVfomIlxtOOvHNk4m0RJgcATjtwkILJSSN7ORhFykGZQnTchju0EykGK01/RC2aZfIkYOOkRza30ZExW+GlOGHFHb5BmTFgPGf/P08XXWOXZXYeq3fDATOsyQHuj/M8GoetnsOjxkAIw9NrtJPkup1lDfu5QQfCdJAr9fNkQa4kyJX4efIgRxLgSvTaCQHuJQw5vaMSzza79ygEk/vfWi8JP3TgZemjYTQcTVSa45tKm1y7NeI5vZ+T1bafQYUasq4AgDAFWztgZebocYGsUKPPW3C5Ek6VrBxX/UXjwj80WD6k6ACDSLz90x6CvUZ6XvKwnpEPLSEKb9cP0rq5pPZrZiRcVXJD5xG9H9a30rdkE5Iu4Xw27zEhL4vDHD0Fyoe1EDZNumFb2Z0AwIyEeEEHP+DkhFy8oJMXcIi8RrW5esbht+pyF58qWRnNdNwHVKKKZutamg4SRALido9Y6HFyYixAavgQovApexyDSITyOP3VJcm9hUwQdKT82Pv7pt0fdSvfbcTZYkgX9repw4dkPohZ+IwTmBmzyo9/mGLTS3nFFs9BjaS3VBh9Zq8Jl8mxYuTUrhtkdhmwhxwMIZS62ydUffL1VX+n0bl5iTOMP2mCa9OJlPMjgkzuXWrx7EFo0qAxrIWwkwiT4xFqPELNBcd5QWdpzfobNtzbmDXvRNlNUbeFuGAzJWJugdV7WJWI/1cM8LWOvreIYIys/Fqjl5O2+NLgmT37kwSTel/6FNd+4xSndGjHA4AthI/b8a+KRtuCZqYK7TbiEhmjqvi6cTVfVE+ar3duHSAhjNDwbRv+Wd5Yl7kQD4kOttAnbYOhggSm5uz757Hxt7pF50qtkjR80oQr1Cj//BR7Ecrn8FcVau4fhFYNGiO7Awe40sOTbl+//GUAWLnxV+XH3ueE4jEHAQCAVjLP4NqZqCaZAxB7Hq8xBpoj1vjON7p2qiXzejmBHfaOP/3F4Ul3RF9uacdz1EgwMuaTcZAvQWEatXqhLmex2lyTQ6k9gcZwxDFAt2t047jiW0Y9ERp2GfC/W3iDo4IAML7qU5LFP5O3pPPIpjas5qFpygsnKGbPfhl/PJMxMtNGXISRLYRRvHzlgfJ71y97iUmFbtxwT/mx99nd6lf0gkI41R2sS2An320A+1hY4TCg2YPjqlwaJI3+sFEu7C215oTTn7WmT3dIMwDgjBO7wlDebaQYBSCAaUo4YMYRJudMwbJxZ79NEk0zuXcP3B2/baMN/rFeAwDQ5MFvnqF26On46yz1kSR7U3HtN3tm/raz0OBBCzYF8PK0Hp5to3unttfJ4khkNAhhFK9AdaD83g1LX+KQ3lXr74pdDhkEJ7GdnKTxJt1Ylx56jsRZM9Lg2qmRzEbAuNgJQq8pv2HrsXG3AECEhu87YGkaMWhOMkIWcBmDJ7oTFKjDj61BXF14daZufyZvssG1feBuF6Hhi2b6Mi/nYg/htQ30h3W0NXjpkxMFgyLn7Pv7wfJ7vXxl9IjOC/uMeFU2YnXThyBpDIQ75IJJg9e+QWFYC+FZN6PNB3EVHPEIVXun3b/xqn/yA7brv75f5myN5a+04isS28kb3fjwWAapIcVDwtl4kotijI2uHb2bwaYef7+68Go/Tw4Ae0w4VQCZg+Ijw2XAghT021LGQ2VoeTpSD0ryGiaCyUnooAVCbFFT5tyK5rM0HfYEGwfujrYgbGy9TIWQpPEOPX69hq4b9MDKKSc+dErSGjPnRl+6Sfx5C319JiHtqdCuwbVTJZ6N0EUniyOUYS2ExiCxuZ3+v9P0Z834qBW7Yk4O4hKn7Kp4uHLyHUu3PqExV1/yfAm/iKLDnkAiO/mWNtB5L9NePRw4asWxJxcFAFegmkkIhJysi52gNZ1WWWujnuWOMD5ixYuTE9DO3mEgmJyEflPKmK0hWASwGWiKkrivmPGLQmJyEmIOcPedmoSqndhHQlXxdUX136WIKowDuSgEgGoHrrz8MvfWOvHrNXiXYQhCKjWm09ktu/ZP+3X0ZQTDp014hgplXVj9CeDHyaJWcsWgNnFQGNZW/nkqMl3G80egxYubPLDbCAyEs8WQLYIcMeJcahRoypwb4koW7PrL3um/aU3rLVEIAqSVzte7thXweqiz0zcojL9shnuLRkO2kREHhfHROFfkBtc2rfSi5UkJTM84/MahyXdGGGwA2NIOM9VI1NOUOVEQCJVJIstykZzfw4OeLEDJAjQ/GY5a6aPWgUrayWdCiQQdtuF5Go0udfriduc7oiO5ytsRwRqI20X5vgNShTiZPwotr93p8OFtHbjZMzQzZhbpn7v/hX3Tf9MZfvZtG5aw0QxVz1++M1DFYgh7mSyOXIb1ijAKnwnFUrQ8Df2uFN2UjWRsOGzB/6yiP2ygj9tx7/bkDs2EzYv+MqPyzaL6zb3fRSu5wuzZl9iciq4w/ryJpgclMdIYXTnrBA8Zx/kRymf1HlVLLhpHX1j3bZAjaU6fBQD1brAFYfpA+shki9C9RcRV2oi416KkQhbM1RIPlhGrc4hsERqIEr4z1HDEgiM0VE6+Y2LTITlDY/EeTPhduhKh8aeNdGC0V7G3BvHnTfS/aodMBQFgxuG32pMn6VLKoy+PWHGHH6/IuOhjZHBs1UoWDlrzBpMRIISdIAANDyrU6Od5xO/LiGlKdMaBXzxN7zBg78VHPbs0Y9OVz5fWrCs/9n4vF+cwkyTcfItnf2Lb3OzBOw2JveQYlyau5KIAYHTvkgsmsoie9oMAOGHPxFOfHCy/FwAoGra04yVpaIDcVjJF6O5C4rZ8InZDIAIokKLb8olfFBJpia4RreCgVAE6YcdBjvjIxJ8vbbfrXVsTe4vuuMKwoZUenNSag487DBtb6ddr6GrHIKUP7ZGMtgMac3XlpDujL9t8sNOAV2X14CATJUJ5bP4T6kslXRqhjCQh7AqLgAIJuiWHuCsfURjeOEt/3ozbfT2f7BFqvrnyrymGYxWVr6OLP3la6aIOx/cJb+oeI651js4uPTwxB3BLnLNsg/P7FOnii7075fhHjVlz7NIMANhvwWoe5PSsmP0CAazMIm7PJ1L7KmZaPtxZgFZkEiJWIuVwpgodsADGUJezMC+cFPQ1BML6BF6/R846caV1pI5OF8Mfwd+3069U08esOC4fwITDCzorKt/YWfFItNygNwJfNONrMwg556JPjsG9I0lYPsrCBzsZ8Y+ajIMWJqMHSoh0AXzZit+upU/ae3jIAlzZpsV/FXsMC3c9y6B63v9MEk0JRMy+cFtiW4gxXtdCj0UWDhrxLgddgVoaRyT8kh7flTlbM3X7T5TdDABuEh8048UpCWhkdxakoDJ5fwUMITRegX5TguYlE8wEBXakC4HPwLVuDACHpt0/1xK2WjYk5Mq9s92I2v2jxFJI0rDXSL98mt5vwuTQaiAAYDx7/4tnc680KwsAgMbweTM9OQnyLp6nEgPWO75PvvhkcaSTYCHcsGFDaWmpVCqdM2dOXV1dYi/eCxwCpqnQA8Vothods+GXa/B+M76gFgTJ4n03/0mKwbrqhz+yyR4WjwgYyZKFeseWhDcvSMFnTfTQd4DLgBCFT8VZtk3v/C5ZughBz2PuzMrXj064NepNsKUdpilRj27l/aRUhirUCbssm4HmadEvi4gccWKuOV2JDpgxALhFycnyq4yu7QNRofACaAwb25n+EW4sJGmoNNMvnaa2dVw4Ig0VxXWb+EHHibJzBaW3dGAeA83W9PaoOH2nCcSU8AoHpYFDQCKFUKfT3Xrrre+8847BYLjmmmvuuOOOBF48FhBAkRTdkUesykLmALxWQ19Qf4cmmNtnPWpV5F295RGB39b9CimSRUb3HopOfDir0Q9jUfYDDcb4h474kouSlNvqPayR9uwRntOyi0X663IXA0C9G8wBqLiIQ11/0PLRtZlEwl1dkrjwszxiVQ4h7dXjJhaKZcgdRtHe1Fp8R3qQCS0f9L+Fl8RNwrqWkTp/dIVhWwf+ZxX1bVtvTgyDjNijn3jqvztmP0oTTAA4YcNNbliReYmHT+/ckiK7anBaOCQk0rW/qanp5ptvnjFjBgD8/Oc/f/7557uf4/f7LRZL9+MymYwgzlNlmqZpuo8Gcw0Prk2HJg/6vJm+QosmKs579+CkOyewBcu/e3TTwmc9QnXXt1hMuYxfanTtTJZe2Yf79s5xK07h49ir1tE/kvCWjAji/fg0xpvaIN6UlXrH1iThVCYSdn/SmJFQ+bH3t1f8ngIUofDmNrwsDREowf4NfCa6IRMzoIcPmpAHoEAM2UVwyIz3GCHcjystSsabdPgXhYiBCIX6Rp15bV54dZg1gBYjjDHGuM5J79bTszUjyYij9+FKK1T1ZKOJi+g3kKBGAQAQdGT+nr8dG7fGKUwGjK1B2KrHt+cidq9PdShidfiqCjS/HmzXHozjev4vUJC4SKQQzps3b968eQBAUdQf//jHVatWdT/ntdde+9e//tX9+ObNm7OyzgtPoSjK4wn7cd8jljQMWJ2G1nVwWt3UQhXZ1WKyP/v7/McQAAAgAElEQVRqL3Cu2vrEF3Of8p1fsyKJd0WL/X0pexZcZK+sP3xZCyx/UMuL6XlyuVxcLpfD4SS8GSOCcDjs9/spKqbtJJLG3+g5jb54E17Q7Y4t+coH/X5/9/dmnPm0XV7QLMwEv3+PlanmEBpGuKcT+w4BcE1aKOKme5gbAthstkgkwmQmoJMWMECtQmtb2X6qj4NFBhvELM6uDnq6nMSSxQ7nF5rjL5wtfaj/bbsYJEnSNE1R1KYGLAiF0/jDfUYYoXGth3nUwTCHEpN4xe/392dw786MM596WaKjKXOw309jWNfGmaWg+HSk96e6zfWNnD8zFMQACX36LwkZtFp9PfbNHpFKpX0eLRMf7L1t27ZHH3108eLFzzzzTPd3H3nkkaeeeiqW61AUJbG6BWS/ppwCgHvE8GUz2mBk3pBFcLs8n41lK4SIXLn/L5sWP9+1nKFAMKXN/XEYGmWC8f259cXY4RHem0rEEmXPZrN5PB6Xm4BCiSORUCjk8/nk8ksX/wzT8EkDbQQsiPNhsXgO8DgqpbS4+1sir2Fc09Z1y18V8AW2EK5y418WIgErwYHky9PRZOVFXRQIgpDJZAkRQgBQA9wtw+/X0n1O53l1Bn7rLDFBxVJwkFq98qz302y/zqQsSkjzuhMVwujQ9p0TVXDRLPWAJ9PpG14Sn7DhymhmAyYkqhoJxlgQ7zN9cdSWs2Ut279a9gqfJwSAHQbMZ+PpyUwEvYkHjcNW/e7JGX/hsQfbXxRFmCqVIIHfQC8k8rHCGP+///f/nn766U8++eT5559PVAfuJxwCVucQGj56pw5bg+eNASdLb9KlTL3yhyfZkUDX42nS5Tr7xgFqjyuMvxqxZo9hiJeEf9dSTX2KSm6zb0yTLe/xrelH3q0qXuHjKwBgUxueq0WJjUkAgClKNEU5qON6mgAtSu37p5Cw0SwN+kaHMYBavuSkhFl25EWCHgxvFpLGO/X06zVU/aCn4uwdewhvaKFfqKK3deAByu+TEFiR4Jz9/9w/9ZcBngwA2n1wzIqvy7i0Y7HJtVPMy+OxtYPQyCEkkf1w//7969at27hxY3Jystfr9XrjLg04QCCAhcmoQgnv1+Nmz3lvHZ74c6sid+GOp7vGVKjFc7yhZl8owXEUnTS48PaO4dtnRhCOEPy7ljb2acPGHagPR5xK0bTub6UYTsqcLaeLVgDAKQcOUlAes2U3RjKEaElPNW4Gmukqoj9BGtOVKEjBaQdmEHyVfOH3SjSu5qsENq937CH4uIH+pJF2DoNgJFsQr2+hX62mj9vw8K+ZMf3w20ZVaTQ1EknD+la8JA0JL7VUwYB1jm/S5dcMRhOHlEQK4c6dO2tra2UymehHEnjx/jMpCd2URXzVSh/qGmeG0P6pvwpwZVfsfr5zbosIVorsqjb7AAZL7THSO/XD3eYxzDEH8Ht1fQ/QbLOvT5Uv794FCExNP/LWwSm/oBisEAXb9LAsNcGlliRsuCmHYAxARrRYuCYDafl9/FsCwdXpxNYOCEYgVbb8gNiXe3ad1N2e0AZegrNO/FoN3mnAkSHaWHGE4OtW+rUa+oRtiOPiYyS9/VCy8dSh8nuiL7d04HQhFEsv/fjZvccYiC3llw5wA4eeRArhE088gc8ngRdPCOlCuDMfHbXCN20/PcEYoV0VDyHAc/a/2Jl3JkW2xOqpDJM9hFgkip0GvG1sXdhX2n34/bq+b0YFwnpnoDpZ0kOW7eKzG/08uS51KgBs0+NCCaQk1EjBItCqHGIIi9qzCLQqh8Fn9lGGk/lQKIVtBsxlqWTCSf8tLak4+CoMbmf/caeUHuSdUksQ1rXQr1RTR4c6NUzs8IKuWYde2zXzd2EmDwAa3NDshitj2yFvtX+VprhugBs4LBiWpueBRMZGdxUgdxj/p5EO/eiQSBPM7XP+n8BvnVn5RvQIixBpJPN1jq8HtDF7jfSOsXVh/NQ68Qd12N+PAt46+/pk6RKCuNARiRd0jj/9+YGpvwSADj/UuuAKbSLXbQihazPRkJdWkLIvHTrWCwuSUZ0bdF5Il684zGxGEX9u884ENi9GojulnzXR7gEP7gdzAL5opl+vpk6OkFVgJ7MOvFSXs8ioLgUAfwQ26vB1GZcu3QMArsCZEGlXi2YmvEkyDkxRolXZxP0ljIlJQ7Qxcj6XnRDCj+4zCg5a2/RTAbAIg/39/D8qbXVTTnwYPZKmuNrg2h6hPBe9UCLYZcDb9SOqYw01J2z40/6l6QlFrGbPwTTZ0u5vTTn+QUPOQqc4FQNsaqMXpQA3oSVIK9RQKhsOHR/yJGhuXzWeQ8BVKfBNG83jZAo5WWsnTJt27F1uyJ3YFsZIjQO/WkN/raPPOjE5ALPKVg/+tJF+o4Y6bR9+e1yXovTMekHAdmzc6ujLjTo8Xg7psVWTbrV+maFYkSiN4DCgUIqWpaMHSonfljKWpxNFMpTEhWsziHuLGNmiIe4Ul6MQAgACWJqG5Bz0eQvdOaKSLP6WBU9nth0YV/0lAHCYSSrRTJ19YBeFALDbQP8wtkcaA/4IfNVMr+/yX9Y3Wm3rk6ULmYwLbdhp+iPJxhPHym4GgEozZhOJFC0GQvO0aEHKMOpxc5OJfEkfP2CxFMk4aL8JZybdWB3YXZc5d0blm4ltXuyEKXzUgj9ppP96kvq4gT5soWMv4n0xSBofs+I3z9Dv1dFnnCNPAgEgTX9kXM1XP8x5IppE5pgNu8J4XmyzH0+w0Rtq1Yrn96cBCCBVgOZqiTsLiP8Zz7g5hyhX9pDXW8OH2/KJNbmEcugixYZRtxxkEMA16YiB0HrdT9bMIEf87YJniuo2F9ZvAYAMxfUdzi0kPbCLQgDYY6S/bx+JfW3wqHbg12roU/b+fkth0mZy706TXXvBcX7APnv/SzsrHiFZfG8E9pjx8rSEbdqoeeiuAjQvOcFON/0EAazMQkl9HX2WpsJBC6aIfB5b+012qsTdVnpmMJJx90KEhnoX3qTDL1TRr1VT2zpwkyfunUwPCTsN+IUqemMrbfSP1F4pc+nm7HvhhzmPe4QqAHCG8Q96vCKDiLF8WLP10wzFij5XYGYSMF6BflXCuLuQmJ+M0oXoko9+ngT9qoSxIpMQDmDV54syLEL9hgoEcH0G+k8jvbkDlv5oPfbzFZsX/nnZ94+FWfymzDkq0Yw224Zs5a0D3Zj9JhoALU69fKcmF8MVxpt0uC5BbhEt9q+SpQvZTEnXgwjjuXv/cSZ/qUlVDABb2vFkBUriJkC2GAjN1qLZGhgqH9He4TDQjdnEu7U4HH8EgISNKlSwqQ1fnbK6Rv/PrXP+dN13j9hlWXrNuIFoarxYgmAx0nuNIGShbBEkcUHORTI2yDmI19OwhzFu9kKlGde5YKQX0+YFXYu3P31wyj0mZSEAYAzrW/FsNVLFVuTSE2zwBpvLkh/pw60lbJiqRJOVRB9sCghgvAIVSol9JjhgwoOZbv2yFkIAYBKwOof4oJ7ebYQ5P+Zfd4u0WxY8vXTrE2G2MKS6sbL5oTT5chZDOtCN2W/CGOgrx7TwRzBApZn+QQ9x5dHuhSBpNrn3TM9+5YLj46o/J3DkZNkqAGjyQIcfrk1PgG5p+Oi6DELT11iFwUHNQ9dkwBdNffmGp6tQlQO3BfMF7LS6yPEdsx6Zt/f/vr7qH9FVyDDBS3YWJDn3GXlMJGODnAsy9jl1tASh0oItgZGtf1EYFLlw1zMN2fMas+ZFj+w1YQLBNGWsj3ST5ePMpBviXQ6mCtAMNSqSXnrx1zscBroiGaYkoe+b6IHIc9kjY2MucAi4NYc4ZYdD5p+6gUOa8f0Vf5y77+/pdrNWMq/Z+sXgNOaACX+jG7W1uePCHoJPW1nftsVXTaJ3mq2fpMqWsRjnLQeVtvrSMxt3VTxMI4Ki4ds2vDT1onW6Y4RAMEtD3FM43FUwSqkMzehTEahoWOH3HZCiWNNq+7JdVXCqZOXCXU8zI6GENzKBBCJY78en7XiPEW9ood+vozfp6NGhggAws/INP1d2bNya6EtjAA5a8LXpsVY3cfpPB0iTVrIwxtshgHwJurOAuLuQKJH1VwU7EbNhcfLgedSPCSEAAJ8Jt+XBAQucdvzUGSyK/J2zH12469kJjHKze0+QNA5OY45Y8Jb24RiFOWhQGO824LdqUasvkfNBb6jV7jtxQZoMFhmYt+f/9k6/3ytQAcAeM1bzeqtQGgsaPrqnkFiYgmK0xwwHFqX0sXhhNKzwgDVTwi/R2TeeLrrOJsubffDlhLdwjFiYcPozhb1hd8VDGCEAoGhY10ovSUWS2IpoYsAN5g+ylLcgdOmdTTEbzdYQvysjbskl0oUj51nviTEhPIeYhW7JRlvacWMXJ/AOzYSD5b+4euffcwTzG80fDVpjDpnxK9X0XgvDPqwn1okHY9zkwW+fobfr+54b+mI0mj/IVNzAIM6zk8w69Gp7yuTWtOkA0O6DI1Z8VT8K0DMQmqtF9xQS2qGOFIwXAsHNOURWn7zYF2hRnQd4gjXtjk3hiGPftPvE7o7SM+sS3sgxeidTd6Co9put8/9IMrkAgAG+asUaHord+dni3gsAKlFFL+ewCFQmR7fmEr8rRQtSYpXYYc7lbiPsiooHq7KJT5ro1dlE6o/JRBoz57LD3vuOrH+sGFyBMxLeQOXavwB7CPQuxiEHypTQZXIolQ9lLpJBwEPiEzZ83Ap9TpnWOzbf0SBpSZadV2Yyv/F7uaN5w9IXAcAfgS9a8DVpSNTXji1godU5KFUwUscFFgG35BJrG+jmODOYcxiwJAW+blctks9vsqwt1P76h3lPXPvtg05xenvK5AFq7RgXkGRvmnXw5e8WPO3jJwEABvhGh4MUviUj1tUOjcMNlv+UJP8O9WSZQwilCWCCAhXLgDuC9jpiY2xFeB5pAliZSXzSTJu6lKM4k7+sLWPeqg6qwfg2wGAngmn34c1t+B+n6A/r6JM23J/aqsMQDNDkwZ830S9W4R868ACpII0jDab38tR3IPhpw0fi7ig/9uEPcx6LMNgYw1ct9DgZ9DmuTsaBO/OJkauCUaJamB3/HmmRFE1XwU7HjVbvUU+w0ctX/jD38bkHXhB7DAPRzjEugB+wL9z55/1T77Mo8qJHtumxKYhXZRGMmMd4ne0rCS+/+1xfzIZZGuL+YnRnATEpCY0+FYQxIexOtgiuSkFrG7GzS0zu0fFrNIJpSpfZYN88JK2iMW7y4HUt9N9PUuta6AbXCMvz1B1nGHbo8QtV9Id1dLUDUwNpE22zb+Cxk+WCSZ1HGBQ5f8/zRybe5pSkA8B2A6YB5if3sYenCtDdBYRiVBSOZBGwuk97pNOUaJKS30quOWN8G2NsVBafKL1p4a5nWJHgQLRzjE6YVHjhzmfO5i9pypwTPbLHiBtceE0OwY45hiFImtodm3NVP+96MImLbs8nHiwlFqYgRSKiiYYtY0LYA6UyNEsNH9RD1/qFlZN/sSSQ1258PxKxD2HbwjSctOH/NND/OEV9q6N13hHmVuOPwGk7/qiefqmK2mWgB6GEW4g0t9k35qvv6npw6rF/e4TJtbmLAaDWhasccENmH/3dCqXo5/mEINHVCoeQ6LowM34tnK5EafL51iDR6tgKANWF11gU+XP3/WOQU3JfXmA8e/+LLnHKidKbogeO2vBxO/wsl+DFE8lXZ3onXXEth5nUeSRDiO4qIDJFKFZ/05HMmBD2TLkSzdPC+/W49ceiihihM9P+MMUnddQ8Phw6ti8ClRb871r6xdN4cxut8w59k3qExtgcwCdt+OtW+rVq6v9O0V80043uwfsGa01vp8mu5rLUnUcy2g+lt1fumfEbALCH8EYdXpmJ+H0ywU5KQjdlE/2MtRiGsAhYk9uXdeEcDcET/rLBvNYdcgDA/qn3CXyWwaxZeHmB8dTj7wl8lr3TH4geOOPEOw14TQ7ElZ/F4tkfDJtTZT85VBfL0K15PWceGJWMyB7MIhCPicTsgQ22HC9HKzOJz1vozpgKmmAyxv/diOyyqr8N5J3jwxXGh8z437X06zX0bgPtGAaOpl4Szjrxtg78Xi393An69Rp6XQt91IotQRjk9avJvTtEWtMVKzqPCPy2WQdf2TnrkTBbGKHh82Y8X4vS4i+0hBBamIKuySCGV9q0xMEiYHUO6sO6cH5qBsFZuKv5nWAEKAZ727w/lJ7dkKI/NhCNvJzhBRxLtv9RZT7zw7w/UAwWADR5YFM7XpNDKLrl8+wFkvbUm94t1P6aQOd0b7YG3ZiFWKP1ye6JYa345XJqvpRgIOAyEQLMZQCDQOwu2u0hocZB1zhA5x2Q8TVLBLflEmsbsTWE52kQANBsSU7673fo/nrjmU/ri1YNwD37jjmAtwdgu55KFaAUAQiYSMQGPgMELBAyQcCCgXuyQxTo/Vjvhw4f1vthOBQQB4BwxFFvem982h86g6IQxnP3/b268Jpo6qlN7VjJRVPirz7PQHBdJupPqfcRAZuBbsmBtY3QEqcf6fyMVTsbHtrQtG9FbgXwk7bP/p8Fu57bO+M3ranTBqiplxtpHUdmH3ipNu+q4+NuphEDANp98GULvSqL0MSWR62TBtO7KnGFmFcAAAyElqWjSfH3iJHOsBZCBZtW/RSn2cP/jYgF01TENBV4SHzWCdUOrPMmOE+gigt35cPaRhyIwFWpCAHwJVOl8vkn2r+68pDlyITbgpz+RV8PAO0+3O6DzoRSnbAJELGRgAl8JvTu+cVjgoCJ+EzgM0HABAEL+Of+6tyfRWgwBnBU9jp82BbsdrNhQK3x9RTZYhE359xrjGccfgsjdLLkRgA4bMF6P747P+5NES4DVvU15G7E0TctRARrUvoDR3XP/behcHWOwqgq2TnrkRlH3iyu/ebA5Huc0vSBa/Coh0GR5cffz9Tt3z77f6JVBgHAHIBPmugVGUSMJZY6sXgOuAMNUzL/AQAcBtyUjfqWV2GkM6yFMHZELFSuhHIl8kXgjAPXOHGLJ2GKKGKhO/LR5830Z01wfSZiEaAS3/zFlvfeO/iDz/FwcmrGpLufKr7y5uFvUg7TYAtiW6yn9/DtRQWSQNgWRAPq59l/9M4toYijVHHOgwAwnnHkbZX1zOYFz2KEOvywy4jvzI87lZqEjdbkEjEmLx4dRLXw4wZojccOLeHmZcmv9Nc+89KrJnR2VyQcTCmadPvVBau2PaFLm3Z0/M8CXMmlrzICsbac/eHVJ3Qn9mGaTi2btuDXz6rzE5aFXOpqm7f3b16het3yl0Psc3XE7CH8nyZYmoZy45yTh0lbrfHtcWmPMwiOhI1uySXUl9OD3ZURaSPsBQETpijRbXnEw2VoaTqRkAICAMAm4OZsgs2ADxpoly/44S8W5Tekv7LsP1//8sADkx8/+sJTh/95f0JuNMwJRMAaxOYADHMV9IV0TZZPipMfPLcpivGMI28n2eo2LXouxBEGKPiiGS9LQ91Lo/VOmgDdVRBrCv/RBJuBbs1DU5TxDRdi/4TqP7x7K9J8fMN/19+2+RbFon8+/9xfkq6jCNb1X99XemYDgakBavBQYaw7+f7tsxfhgk9v/urzNRuu5ZZ/fM8C3Yl9Cbl4XtMPy797tCF74ba5f+hUQU8Yf9QA8zRQLI13rKOr9S+kyZaJuXkaPtxdePmqIIyaFWF3BCw0VQnlSdDsRQdNuN7dXzcNBoLrMtAOA7zzwX8Kgrwnbng2mn9hcvq0v1332u0fXXdPqbpu7m9H6zx3BEHRgdMdf8tV3cFnpwAAYDzz8FsKe/2WBX8mWXwM8EUzXSpDRfEMHAhgqgotTkXDs5rSIMAi0PJ0yBERG3U4EImpK+35919vKr71Z1Pvib68suhqPlvwwutPZX526kz+0ulH3yms33xw8j2jKfvMjjefvGPS3TdN+ln05bXjbmQymGtfe+KOd3b257Js0ldx8DW5o/nbxc/bpRmdx5s88LUOpiphkiLux7LJ+l+CYKUrVuaI0U3ZBCf+qkmjidG2IrwAhFC2CN2SS/y6mChXnudo05erAVyhRfL2A3NyF3TNQqQRJ+coi6oN1h/nuaMr+8uIAgM+Y3hZxi/TSOYCdKpgQ1QFAWCHHmOAK+KJnZew0e0FxJI04rJVwU6KZOgXhbGmkWs/dXBu7nlFDKZnzna2NoR8bqckbcsVTx+edEdF5RuLdzwl8o6SBDTtpw7OzTvvI8/JWaA/fbg/s3CVtfa6Tb8l2fwNS1/sVEF/BNa14K91eFkamqGK+7G0eiuNrt1lqQ/O1hJrci93FYRRvCK8gCQuLEsnFqbAaQc+YMJdI+XjRcFFuFsYOMb4bP6yzUUFFZWv5Tbv2Dft1525jsYYTFqtn4UjjhLtQwAAGFdUviFzNH+38M8kkwcAJ+34lAN+URCH+2yRDF2TTlw+AVWXRMZBdxag3Ua8S3/p7EY9CgD+cRLZmjqtXTuxrGbdtd8+VJe7qKp45YjfUEGo/wFC7EhA7miW25vE5jqNR8f32/bOeECXUh59FwOcsuNtHbhMjn5V1Jd6Yb6Qrtb4+nVFf7ilQCbj9LOxo4TLq39zGDA5CU1UQK0THbLgeJ3Co6SPn7n7vXdWTrylc1FodOsbHfVLSqbYJYqvF/8tr3nHoh1Pt6dMPl52y7CqUDrqMXn26l3byzP/iggWYFxR+brM0fL9wqfDTB4G2GnAp+ywJgdijJ1nEXBVGjH58nMlvyQEgnlalMon1rXQvovXEU+bMHNXw7YCdXHnkf3Nu1BK/n86BCsycDTWjWKwT5StqstZOPH0pzdsuLc278qq4utHrhymjZuxs2Hr6sm3dx7Z1bA1pWxa7550Ar9N4WiSOxrl9uYkexMvYHdK0m2KHL00s7lwiV2WFa0mAQD2EP6mDQcpWJMbd5hElHDEcUb/l19OuHtlQUFf/n6UwvjTn/40aDfbtWsXxnj+/PmxnIwxDgQCAkH8oc6XAiGk5KEJClQoJUgMliCKSw9V2cV7Nv+rtqEyU5HDZrBPtB995rvHRYs1ebOWcZgKQMguy6rNu1Jub5116BWB32aT50RYfXlmI5EIQRAEMcq3ry8GxpiiKBYr1gwZTv/pWuOrE9Ke5LLUgHHFodfkrtaoClI0rNdhgx/flkdIYysuoeKhW/OIvL7m4E4Ifr+fx+MN2wdAzkUTkghzAC5WLEydW7r23UepcDBFmk7Rke213724+/kVT74pS835qhUDQJrgnECQLH5bSnlT9txk48mKg69xw267LCtEsDHGTOZImqyrcko+euchFiaSJSkhKrSlZuMbB1+55k/vitWpXU8jMKWwN2W37imr+WrmkbcK6zaLfOYwW2TQjDtdeO3hKXefzV+iS53aJkgjZSk0wQQAGkOlBa/X4fFydF0GIYqvevw5ZOxArf7p5TlzVhcvS8jnHVBCoRBBEGw2exDulYCFfOw89dRTNE0/9dRTsZxMUZTNZlOpBnxFZQ/hQ2Z81IojMZv2wn7v3vf/WrPtS6/NKM8qIq58kDMxI5PxxqT0pwSctM7TuCF3WfVXBQ3fNWbNO1G2KsCVxtWwYDDIZDJH1kCQQCiKCofDPF5McwhPsPFk2zOlKQ9L+aWAccXhN2T25qgK+iPwaTMtZKIVGYgZg6YghCYp4Kq0oc+sYbFYZDLZMH8AMMaHLHhre8+OxG5T2w+v/W/L0V1k0J9cOkm9XFIweVW6fIUzjNe1YiaCazOQ+Pw0rbygo7RmQ0HDd/UZc44UrqAkmsH6KInB3taw/fX/bT2+F9NU2rgZV/zqz8rsYgBgRkIKe6PGUq0216gtNWGWyKgqMimLTeoSh6TnwEqv1ysUCgHA4Iev22ghE5al9bH+n4QNM9WRT6ueSRUlPzT1vv58wEHD7XYzGIyBWAt1Z0wIz+EO4wNmfNTSxzpHVQ58RL8ng/Xh5PSnJNzzSrsK/NYJpz7J0u2rKVh+uuT6MDPW1eGYEMYohN5gy8n2pws19ymE5QjjikOvSl1t313xFMnimQPw32Y8QQ5zNDE5ugiYcF3mEC8EOxkRQhhF58XrWy5dRSsUsR/X/SFVtjRVtpzGsM+MD5nxVak9VI4V+K1lpz7Lbd1Vn3flqeLr451EDjkIY0HAKvIYRV6DzNmmtlTLHK12WZZJVWJSFRtVxZ3xD73g9XpJluCABc464coUKIm5vm5X+EyYo0UT5NSTe54TsgVPzHyIGCE+X2NCCDDoQhglEIFKC11pAR8Z99fij8AO3XZWeG2m+sl8WdoF74q8hkkn16bqj1UVX19TsDzCvLSRekwIYxFCT7DxVPuz+ep7UtllGkt1buN2Xsj53fw/kSxevRvWt9JLehpnu8Mk0FQlzNYgHnO4DBMjSAgBgMJQaaZ3G3Hg4lZDAAiSlhO6P6bIl6TJrgEAgx/WtWI1D5amoQuqJZAkyfNZp9ZtyGneWZezsDVjplmRH90nHFYwKFLsNYg8BrHHIPKaxB69yGsSeY0hjtgt1LhFWrc4xaQstiTlRxgx7fJRGFq9UO/CtS6axKhEhuZpEDd+x04BE6ar0VQlARD+393P8ZjcP876PQONGA/RMSEEGCIhjELScMxKHzD3JWfmKdNuo+P9MPexhan53G59VurUTT75H7X17InSm+qzF5C92g7HhPCSQuh3HTppfPlqf+HsDrPYazQpC42q0uqia0km95AF7zPBTVko9VL9CCFUIoUFKUgWZ3z9QDOyhDBKIIL3muCQme7F0BCKWE/o/qQSV2QlrQaACIbtenzChico0DQldG79kSRJ0zSHwxH4bUW136QajkvcHWZloUFdZlCPsyTl0UM0pgt9ZrmjRe5skTuaZI5WkdfoFardQo1HpPUI1R6R1i3SuoUaKjbZ68QTxvVuqPdAswcruZAnRinMQHYSvw8PpYSNZqrRpCTEIsAb9j2x61mVQPnYjAdGkArCmBBGGUIhjIBy/yQAAB8TSURBVIIBahx4p562xFlY1OQ5WqN/pSFyX7Zs6nQVEnQbx5LsjRNP/VdjrtKlTq/NWWRUlUBPmxVjQtijEPL9No35tNZU3RGo/DLJdbMrlyedblSXWhW50ZGRxrC5Heu8+JacSxtUUgVocSpKFw4vCYwyEoUwijsMuwz0cdtF0xySlOtk2zMCbkah+r5o9h83iQ+Z4bgN54jRTBXS8n8Sws6/YpM+rem0xlilNZ2UeIxGZbFBU2bQjLPJc2k0UC5FBB1ROJoU9ia5o0XubJY7WiJMjl2aYZdl22WZdmmWU5La53VqgAKjH5o8uMGNXWHIFaM8MeSIzxUF67QRxo6Mg6apYEoSETWHm3zm/9n55yma8b+efBca2Go9iWdMCAGGgRBGwRjXuWCPEbf74viiPIHGkx3PB4glx73XlsiImSrovtrgBR25TTvyG7cxqHB99oL6nAVewXkfdkwIO4WQQYU1ptNp+mOphqPcgFOvLtmoitQw28rS/pfHzez6VyEKPm+mEYIbsghOr2OjgouuSEbFUhi2SWJHrhBGMQfw1g5c7+q549B0sMbwIhnxlqY+wmKci5cI0XDMig9ZQMrGUxVUNp/mcns2InBCXo25SmuqSjaeFHotTkkqyeKH2AKSJSBZ/DCLR7L5YZaAZPFDHGGIJfSI1LGY5aIwIyGV9azGXK0xVSVZ691irV2e45Bm2KXZNnlWn/PsUzRYQmAKYHMQTH5sDiKSxioeZAhQngSlCi5UqriEUM1DM9WoTI46fbyqLDV/2vO31cUrbyi8um8NHlrGhBBg2AhhJzov3muMI1VbmLRV6f/KJBR2xv1HbLwsEapQIS2/hzOTbPX5jduyW3bb5Fn12Yta0mdGLYhjQsi3t+baTqd2HNVYauyyzLaUye3aKUaJstrwMk2HSlIeYTN/CjjzkPikHY7aIF98rk7IxeAxoUKNpqsQc6j9QntnpAthlCYP/r4dG/09RdZj3GL9xOjeUZL8cLQMUBQawxkn3meiQxTM1BDjZZdw9+WG3GKPgUX6OWEfi/Sxwz42GWCRfnbYH33JDbnFXiNFMF3iZJc41S1KdotTXOJklyil01rPDns1ljMa02m1+bTc2WqXZRnUZSZVsUlVEruD2wU4Q9gUBHMAoj8dYSzngJqLVDxQ85CKi3vfsYhRCLV8NEtz3pQOA/7y7Df/qf788Rm/m5o8qW+NH3LGhBBg+AlhFL0f7zXiM46YvjVMk/WWD2zeI3nq3zUECg6YIImLZ6mJrJ4mpgyKzGg/lNe4VWU925JW0ZE8yUXwSL40wpOGOOJo4c3LARbp15qq0jqOpuqPAhXRp0xuT5ncoRkfZgsBwOE/dUb/ikYyN0u5GgEDACgMtS58wobb/FAsRRPlvRkF2QRMVRGzNMDtvQzV8GB0CCEA0BifssMpO2719BBlYfUeOWt4LU2+NEOxsmvSR5IkW734iIPV7sOTk1CuGCXzgNGPHVBe0CFxdUg8erFbL/F0SNwdYo8hwJG4xcmcsFvsMZmTCoyqEqO6zJyUH6+FDwBCNJgDYA5iox/MQWwKAIdAKh5W85CaByouUnIhrqlX70IoYaMiKZTKL8x45wi6/nbwZUfQ9afZj2oEw2v8jIsxIQQYrkIYxRqE/Sb6rBP8MWQftnkP1xrfVIvnpCturnGx95mAReAKNSoQox57Nd9vy23eobLWsgJOHunjhjzckJtisIIcUYgtCnElQbYoyJU4JOk2RY5DkhGLA+rwhEGRIo9B4umQuPVij17i7pB49Oyw15RU1J4yuVUz0cTTdNoIKTrQaPnI6jlclPwbGX8cAJgCcNyGqxxYxYWJSahI0lu6KTEbTVXC5KRh5BR6SUaNEHYSpKLOkLjBBUHqp74TJm01xlcoOlCoub8zErfTRmgL4SMWaPVhaxC0fEjjozQhpAlQjBmCegFhLPCZJR59mC2wyXPi8r4habAGwRrClgCYg2AKgC+ClVzQ8M4t+NQ84PXPN6VHIRSzoVhKlMghld/Drv6O1r2vHH1nSfaCO8atYRIjyTWmO2NCCDC8hTAKBjD4cJMHN7mh1Yupi3+RJOWuN73rDtTlqu9SCKfUufAhCzb4IUuECiSQJ+65S3fdGmVHAuygmxdyc4JubtjNCzplzla5vUnqbvfylXZ5ll2WbZNl22VZPr5iwD7xpWFFglJXG5v0Mckgg44ApthhPwCwqCBBRxBNs8gAKxKIzsr5AbuXr3RJUtyiFJc42S1KcYm1Xr4y6jrU1UZo8eyvN72nEE7MUf6cBEGVA5+wYX8ETZDDeAXIet1f0vBRuRKNl8cUTT+sGH1C2AmNcbsP6lxwxknbggAAGLDB+X2z9b9a6aJM+UqC4HZ3lglT0OEHnQ+3+XC7F0RslCqAdCGk8UHBHVhLrz8C5iDYQtgaBEsAW0PIH8EKDiRxURIXlFzQcEHGSbC5uasQitmoSAolMpQm7Pkmeq/xpcNvm/yWR6f9pjgpP5HtGCLGhBBgJAhhV4IUbvZAoxsaXLQz3PM5Dt/JOtO/uCxljvI2ITczQEG9G9c6ocmDVTzIl6ACMXQtoBiLjZDAlMTVrnA0yR3NCnuT3NGMMG2TZYc4QgCIMPk0oy+zwgiD4xUofYIkH1/pFSj9XFmPfq0EpsUevczZIne0ypwtckcLP2BzilNDbBHF5EQYLEBEmC2IXpAiWBghks2nCLZLkuISJXuF6l7m4FEhJFF7o+kDkvaJRPeYyKJWLzYEIE+MJipQlrDHRp0DIZQrhhkqlD1iK26PYiHsisGP61zQ4MZ6HwRIe4P5A6e/OjvpZoVgFk1DVyHsCgYwB0DnxW1+aPOBN4xFLCRiYzELCZkg4YCQCWIWErGxiIVi3AXAAF4Su8LITWJXGFxhcIbBFcauMACAkgdJXJTEASUXJXGxhDXgXlZ0wJ2lkmh5qFgGF9M/AHCHPR+f/mJz0w+ri6+/sfDakb4Q7GRMCAFGmhB2xRaEBjfd5IY234V7pxhTHc7vWm1fSPklmYobBZx0AKBoaPFBrQvXuYBBQIEYCiQoTQDhUF+cZfgBu9zRwib9AMCM+AmqL7VPWVRQ4LMIfRaBzyL0Wzkhj4+v9AoUXoHKx1eRLI7U3SFztkhdbX6e3C7LdEgzHdJMuzTDLU5OVHSXwVPfZPk8RDZb4Kbm0Hw1j0gXoHQhShdC7+6gTALGK9B0JVLyRqoERrlMhLATkoY2L2714iPGun2tH4QirlTJihTFvFiqxUUweMLYQyI3iT0kuEnwkOAhsTuMPBHMQsAizskhi4GjBmIOcc5ix2UgfwS7SOQmMY8BEhZI2EjMAikbJGyQsJGEHWui9j4jYIKcgxRckHNAwUVyDsg5YDcbtVptL3/lDnu+PPvNurpN89Irbi+7Wc6TDWwrB5cxIQQYyULYCcbYEkStXlrnhVYvuH8s3kTTwXbH5jbH1yJudpr8Gim/7KdCFgGoc+FaFzYHQcbCSi5S81ESF1RckHPQUDk5MihS4LcK/FahzyzwWdik3yVOscuynJL0zrz4/cdLgi0E1iBtch+nwl+zcXuAuDpJclWGiHNJLwmEkJwDWh6kCGCcgugeuzkSudyEsCsRGn9TV/lF3ZfOsDNVdo1SPI9B9L2AepACksLRjDdhCkXLR4VoOPcLhbkMJGFjCatns31CIBASMLGYjYQsELFAxEIiFghZIGaBlNOz95bBYLiYEOq9xq9qv/muecfs1Om3lt6YLBxhGVljYTCF8HLsY4MGQkjFAxWPKFcCADhC0OrFrR6s83EJYkWqfJnJtbPB9G8aU8nShWrxHDZTpuGBhofmaFAEQ4cr5IwwbCSjygHmILjDtIwNSi5ScSGJC2I2krCxkDkY6kgxWG6R1i3qbXIaF8EI2EJgC2FbCGwhsIewLQh8wp7K3imFbQKGMDlpeYZkBhWhebyed8YIhBRc0PJBywMtH2n5aKy46GiCSaDF6WXztYXtpOGzsxuPNP93oraiVHWFmF/gDGFnGLvCKNSLWf58uIze/YRRl599AQHwmcBnAp+F+EzgM0AQ/YUJPAYIWUjIwkJWAnZSSYrc11G5qXFrvb1xSfbC95a+nDSkPgGjhjEhHDxkHJBx0AQFAgAviU0BrjlwpSV4ZZXlTLVp26HmB8TcPJV4ZpJwKoshZiJQc3EKE5g/2jcoGllDYAlicwBqnOAmaReJ/CTNY4KYBSIWErNBxAQxG8QsxGMCh4E5BGITMMhOIiQNAQr7I8gfAV8EByLgj4CfAn8EfCT2U+AlgaJBzgUFByk4kCd0sPiHA4G9gVCrUjwzWfp7ETcXACiKouCcuZVJgIyDpGyQcZCCg7V8pOEj9khzfhmjD5QkFT41q9AWcGxu3Lap7iUMcEXGrMVpM3Jl2YEIOMPYFQZnGHvCEKAgEIn+xAEKghQKx6yUvUAg4DGAzwQBCwmYIGAhPgPzmSBkIQEL+MxzLy9lMOyXCJIUecx0aqdu3972Q3my7CU5C56d8zg7/gCPMS7G2NbosABjbA6Ef2g9vL9jf631mJiXLuVP4jEKJYJ8Vq/RERiDN4LdJPKS2EWCJwzuMHgiOBCBEIVCNA7TQGHgMCAqiiwCcxnAZpxLSEUg6FxIIYDuiyoWAd0dDTBAkIJgBMI0hGgcoiBEQYhCQQqHaCAQ8H4cHQRMxGcAjwXnpslMxGeCkAV8BuUO1tl9J23eo8GwSS6cpBTNSBJMYjFZQiZIOUjKBiERYUf8mSqJjIP6Vn1tpHM5b40CgN/vJ0lSIjmvSG+trWGHbu8u3X4KU9OTp5RrJ0xQl4nYPQfbRWh8Th0jQAGQNESi/zCQNP7xdxShMYNAHAJHF45cZnQFee73IZxvnWg61UYbKvXHjhpPZksz56TPuCJ91uWzBByzEQJcZkLYFZIiT5qrKw3HDrUfMwesOfKCDEmBRlggFWTTIPGR4CGxl4QQBTFWjMIAoQhERTFMoxAFoR+jmikM4R+daWjcwwVJGrqHSiIALiMqrsBhAIeBOAzgMDCXgTodELpD0w5/qMkbrHMEzth9jTJecq5sXKl6UomiWMplCJhIyMJd969CoZDP55PL5TF9yNHImBB2F8JOWl1th/RHDxtPVFvOaoXqMmVxcVJBoSI3VZQyUsoMdSdEhRsdzWds9dXWs6fMNREqUp48cWryxHLtRCmn5+9hFDNmI7ysYTFYU7QTpmgnrM5aEWFS9e6mamvtccP6WnsDj8nLlmZkSzOKxamp4pQUUTKPKQ5SKERBkMIhCkJ0dHEGQQoAIEJjkgYACNNAYwQAQQowhuhGTZiGPtVeBABgImAi4DKBiYBFAIcBTIJgE5jNQEwEHAYEww57UG/1G4ze9jZ3a4urmaQihYq8cmVeiXJFmbJYwOqebm6kjl9jDD4ZkrQMSdpNRddRmKq1NZ62njnQcfi9U2sdQWe2NDNbmpEhScsQp6aKk9UC5fCsuhCMBNs9hnaPXudub3bqmp2teq8xU5JeqMibqp1017g1hBd69xodI1GMCeGwRsQSzkyZOjNlavSlwWtqcra2uHTHzac3Nmxp9xhIitQI1VqBKomvSOLJk/gKGVci44glArGEI77YllFCcIc8rpDbHfY4g25H0KEP2C1+m9VvN/nMBp+Jy+CmiLTpktQMcer05OXZ0gy1QDlwjRnjsoWBGMVJ+Z0h5D7S3+hobnbpWly6g/ojHR6D1W9P4svVfKVKoFTyFQqeTM6VyXhSCUcsYYtEHBFrwGocBiNBd9jrCrkdAacz5LIHnRa/1eq3m/0Wg9fsJ/3JIm2aKDldnDIztXxNycoMSXrXxhi8hgFq2BgXMCaEwwWn07lnzx6LxVJYWDhjxowebe9aoVorVFekTu084iP9Rp/Z5LNY/FZbwF5tOesIOp0htzPosjls+qo2RoBQZqrSSjK4TC6PyeUxuUyCyWPxmIgBAHwWr8fJMoUpPxkAgAgdCUSC0Z9Oh7P1ZLPH6WZpucwMrpAtkHIkYo5IxpXIuFIFT14gz52VmqQRKDVCNS+GsAqr1bp3716Hw1FaWlpeXt737+5HLBbLnj17XC7XuHHjJk+e3P8LjjEM0el0Bw8eDIVCU6ZMKSoquuBdAYs/TlUyTlXSeSRCU2a/xeSzmH0Wi9/W7jGcMtfYg053yO0KedwhT8QeDjX7WMBIKUxPyU5lEkwek8siWFwmBwDYDDbnIm4pPtJPYxoAPGFv9CdJR0KRkCfs9ZMBL+llEkwRWyThiGRcqYwrlXElGoG6TFms5Cs0ArUi5rC/xsbGyspKjPHUqVNzc3P78KWN0TsJFkKHw3Hbbbft27dv1qxZH3zwgUw2qgI8B47PPvvsN/f/JluQLxcknTE+l1KkWbt2bSyJ5wUsfo40M0eaecHxjz/++MFHHswVFcn48upPT3mKzM+/+TdxksRPBihMBchABFPQpSdfAAMx0sQ8AGAiBo/FYxLM7eu2fvCntYXSsiSevFp/Mndi9rsfv63R9D166d13333s0ccKZWUirviPHU8Vled//PHHSmXfV41vvfXW4489UawYJ+KI/1f/5LgZpR999JFCcbl4FlwOYIyffPLJV198bYJ2CotgPdT+++U3LH3zzTcvln0mCpNgJAs1PUbaYYwff/zxt159e2JyOYNgbH/n21nXzHro2d9TiA5TZIgKAUCYCoeonpNFaYQqAhEAIGQJEEIitpBJMLlMjogt5DF5Qhaf1e9c+RRFPfzwwx+9+/HElHIAeKDjt6tvv/nFF19k9Clj1BgXI8HOMo899pjX6/3HP/7x8MMPi0Si5557ruu7Y84yPVJVVTVzWsULK94pVJcAAI3pt/a+1MSs+fbbb3k8Hvf/t3f3QU2caQDAX5aPGA0mILY1EIIgNEVgKlfRchVFrV5tEfC0GuhNvfNmrlOj9MMK2tHp3Gg7rdgCBWurLmicuyqMsd5UYyY3Gk5AbPWIH9SkIGAMEIPyFUCBzd4fOWkErodL2Dewz+8v5yU+ebK77z677767O+mJ71i/cuXKopcSc39/KHy6zBGwoDS7WdCg0+mYZXjp0qWli17OX10UGhDuCJh3/tO2aRatVsssYFlZ2YqXXy1YczhkWhhCiLJTX5zb3RtoO3369MBnnmiyzLlz51a9trrg9cPBfiEIoX57f7b2r16z7CdPnmSWoTuAyTKDJsscPHhwV+YneWtI/8nTEEI9fd0f/uOdhNXx2dnZzL5i3759e3d+kbuGFPH9EELdvV3bTm1e/saSjz/+2CU/YZSampqUSmXR58rsVfunThIihGwPO7eo3kpTrN2+fTvu7MYcm5NlXDw1WKVSKRQKHo+nUChOnDjh2uAT1eHDh1fOXuOoggghwoP4y0sZxms1BoOBccBVMXJHFXQEfHvB+/pLV2tqapgFLCoqev35PziqoCOgYuEHFy9UNjQ0MAtYWFgo/80fHVUQIeRJeGYsyjr/T11jYyPjgG+8sMFRBRFCXoTXu4s/VJ8+a7VamQUEbogkybdeesf/0f0DfO/J7yZ+WFhYyPhoniTJtxPeFz0aopzsMyVj0TaSJF2TriuQJKlY+IGjCiKEBDzfTQsz3SrDicHFB5tms1kqlSKEpFJpU9MwV3o/++yz3Nzcoe1qtXrmzJnOLRRFtba2snl3By5Go/H56b91biE8iBD/sBs3bojF4l8f9vlfAeMDXnZu8SQ8pdNC9Xq9r+9IX9I9KOCS6SudW7wIL6lfaFVVFYMTVoTQzz///FrAOucWb0+fIJG0qqpqYMynt7fXcU4wwoBrAxOcW3hePPHUIL1eHx0dzSBDd3Dv3r3e3l7OnhH29PT09fU9ePBgoKWuri7sucfeqxDkF2zr6KqtrWW2YdfV1YXGPhYwZFpoi7Xl9u3bDPqdy1mt1vr6+rAlj2UYFhDe0NBgsVhwZcWazs5OT0/PyZOHe5v5cEQiEeO15uI+RtO0Y5YHTdPUcI973rRpU2Zm5tB2oVBIEI+dnlIURRDEaC4ajRdSqbTxsmlQY2O7KSwsLCAggEGlCQ4ONlc/FpBGdGP7nYiICGbLMzg42Fw7NKDp2WefZRZQIpGYzbedW+y0van9jkwmGwj4REOjEonkzv3bLwTPH2ih7FRzR6NzwPEIhkadh0bFYrG53TRDGDjQ0mK76zPJOyQkZNDeY4TEYnFju2m64JfrL80djb5TfYOCgkaTuav09/fPmDHD3G6KeOqXOUHmdtMzzzwzrrfqEeLxeE80NMpsG3BwcR8Ti8Umkyk8PNxsNgcGBg79AJ/PH+H8BZqmCYIYzW8bL+Ry+Wtk0orZKU/5/vd6vkr/rX+QKDo6mtkSSEtLW/Xq6t89lxTwqIcXXzkaGDYjKiqK2Ztj0tLS0lLTl8le9Z8S4Gj59seisMhQmUzGIJoj4J/W/XmpbMXAqNTRSwejY6NCQ0MHPkM8MsKAb69XJEYsE04SOVqKKvfHxc91kz0aM0+0BCaeoT9fLpcf+qIgRjzHx4uHELLT9v0XctauXcv4WEEulx/8Kv/zVV97e/o4An5TliuXy91kmRMEIZfLD/w975PkL70IL4QQZacOlH3pPhmOKTa3f8+PPvrIheHq6+uvXbu2ZMmS7OzsiIiIZcuWOf9Vp9PRNJ2YmDiSUDRN9/T0sHOlFC+JROIz2eeD/Izmtkbj3eq//Vj4786LxcXFQqHQ29ubQScPCQnx8EJb971jaW+6ebf66A+HrvdcLikpYXwUGRoa2of6Mve/e7ej+ablxpFL3xh6r5WUlDCekxkeHm572Ln96y1W292fLNeLKvfX0Ybjx487n/9RFNXX1zfwhvpfJ5PJ7nfe23Ew09pp+an5Glmxr9Gr/tixYyKRiFmG7qC7u5vP53Nhlzcsx4t5nUdE4uLifrheuff4J/e6rPo7l/f9a68gbBJJkszG5xFC8+bNK6+6kFOy516XterOj/mle6bJhAcOHHCHcVGEkM1mW758uaZc/dWpvJYu62VTZZ7u05kvSPLz8729J/5TBx8+fEgQhI8PG49UdfGs0ba2tvT0dL1eHxsbq1QqBz0eCWaN/opbt26dOXPGYrFERkampqbyeLzW1lZms0Ydampqzpw5Y7VaZ8+enZqaOvrtyWg0qtXqlpaW6OjolJSU0XfFmzdvnj179v79+zExMcnJyYNKPoNHrFVXV2s0mra2NkfA8T7FHGaNDvuItcrKSp1O19vbGxcXN+hQm5mKiorS0tL+/v558+YtXbp09AFdZeA1TKWlpWVlZQih+Pj4hQsX4s6LJfCsUYQQslqtu3fvzsnJGeus3NbevXsTEhJccqf5eHT16tXvv/9+27ZtuBPBZuvWrZs3bx7Xo7ujodFoTCbThg0bcCeCzfr16wsLC5ldzpgAjhw5EhAQsGLFCha+y31HXTo7O0tKSnBngZNWq2V8f8IEYDab1Wo17ixwOnnyZFtbG+4ssKmurq6oqMCdBU5KpdJuZ/w84HHv4sWL169fZ+e73LcQAgAAACyAQggAAIDToBACAADgNFYny+Tl5R09enSEN591d3frdLpXXnllrLNyW+Xl5VKpdNjbMbnAYrEYDIaEhIT//9EJSqPRzJ8/f+rUqbgTwaO2tra9vT02NhZ3ItioVKqUlBTOTpbR6/V8Pj8iIuL/fxQhhFBGRgbjd86wWggRQkqlkgtPTQMAAMCmxMREiUTC7P+yXQgBAAAAtwLXCAEAAHAaFEIAAACcBoUQAAAAp0EhBAAAwGluXQjVanVkZKRIJIqMjNRoNLjTweC7776LiooSiUQJCQlGoxF3OnhQFMX4fU/jV2tra1JSkr+//8qVK1tbW3Gngwc3V70Dx/s+2zt/2l1RFOXv76/VaimKKi4uFovFuDNiW0NDg0AgKC8v7+7u3rNnT3x8PO6MMMjJyYmLi3PnDXWMZGZmbty48cGDBxs3bszKysKdDgacXfU05/s++zt/9z0j7O/vVyqVixcv7urq4vF44/rFcszcunVr3bp1L774Ip/Pf/PNNw0GA+6MMIiJidmxYwfuLDBQqVQKhYLH4ykUihMnTuBOBwPOrnrE+b7P/s7f3e8jtNlsvr6+Hh4eFy5ciI+Px50OHhRFKRQKgiAKCgpw54KHh4e7b6guJxAIrFYrn8/v6el5+umnOzo6cGeEBwdXvTMu9302d/7udUYok8k8PDycHykkEAhsNtuuXbsyMjIwJsaaoUtAq9XOnTtXKBTm5uZiTIw1Q5cAN9E07VgIjpEi3OkADLjW9wdhdec/1mOvjNXV1W3ZssXx7+bm5ilTpuDNh312uz0rK2vBggUGgwF3Lpi584Y6RmbNmmU0GmmaNhqN4eHhuNPBhoOrnuZ832d/5+9eZ4TOxGLxoUOHdDodTdPHjh2bM2cO7ozYVl5erlKpTp06JRaLbTabzWbDnRFgT1JSEkmSNE2TJJmcnIw7HcAqjvd99nf+XmP9BYz5+PioVKr33nuvrq5OJpORJIk7I7adP3/eYDD4+fkNtNAcvljCNTt37kxPT5dIJLGxsUqlEnc6gFUc7/vs7/w5fSEaAAAAcN+hUQAAAIAFUAgBAABwGhRCAAAAnAaFEAAAAKdBIQQAAMBpUAgBAABwGhRCAAAAnAaFEAAAAKdBIQQAAMBp/wEHeZsHiwVI5gAAAABJRU5ErkJggg==" }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Sample predicted values\n", "xtest = range(minimum(gp.x),stop=maximum(gp.x),length=50);\n", "ymean = [];\n", "fsamples = Array{Float64}(undef,size(samples,2), length(xtest));\n", "for i in 1:size(samples,2)\n", " set_params!(gp,samples[:,i])\n", " update_target!(gp)\n", " push!(ymean, predict_y(gp,xtest)[1])\n", " fsamples[i,:] = rand(gp, xtest)\n", "end\n", "\n", "#Predictive plots\n", "\n", "q10 = [quantile(fsamples[:,i], 0.1) for i in 1:length(xtest)]\n", "q50 = [quantile(fsamples[:,i], 0.5) for i in 1:length(xtest)]\n", "q90 = [quantile(fsamples[:,i], 0.9) for i in 1:length(xtest)]\n", "plot(xtest,exp.(q50),ribbon=(exp.(q10), exp.(q90)),leg=true, fmt=:png, label=\"quantiles\")\n", "plot!(xtest,mean(ymean), label=\"posterior mean\")\n", "xx = range(-3,stop=3,length=1000);\n", "f_xx = 2*cos.(2*xx);\n", "plot!(xx, exp.(f_xx), label=\"truth\")\n", "scatter!(X,Y, label=\"data\")" ] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Julia 0.7.0", "language": "julia", "name": "julia-0.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.7.0" } }, "nbformat": 4, "nbformat_minor": 2 }