{ "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": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ " \n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVyVdb7A8R/LgXNAZMcFFBVcQ1ChMU2n3DG3UlzAUMfGmhatyW7NlNWo49xpuXOz9GXlzbJRK3PJzEQdlzKNJlEccaFUUhGU/YB4OPv9A8YFHMEDnB/w+7z/6CW/nvOcL4eH8+E85xFd7Ha7AABAVa6yBwAAQCZCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQWmOF0GAwZGZm1n17q9Vqs9kaaZhmwW63WywW2VNIxmHAYSCEsFgsiv/qRw4D4dzDoLFCmJmZOXXq1Lpvf+XKlYqKikYaplmwWq3FxcWyp5CstLTUZDLJnkImi8VSUlIiewrJOAzMZrNer5c9hWR6vd5sNjvnvjg1CgBQGiEEACiNEAIAlEYIAQBKI4QAAKURQqCZycnJeeKJJ/r16xcXF/fMM8/k5+fLngho3ggh0JwcO3asV3SfFef9jox9Ny1+2dITLnf1jj5z5ozsuYBmzF32AADuwAsvvKAfMl+Merbq406x+R66l19+ed26dVLnApoxXhECzYbdbt+7d68YMP2m1Xum7927V9JEQEtACIFmw263W6xW4eF106qnt+K/hwWoJ0IINBuurq59YmLE8V03rWbs6Nevn6SJgJaA9wiB5mThwoXjkh6xabSid7yw20X6l5ovFrzy1SbZcwHNGK8IgebkgQceSPn8477fLXSdF+T2dFD/w3/bs23z4MGDZc8FNGO8IgSamREjRowYMeLq1auurq5arVb2OECzRwiBZsnLy6v2jQDUAadGAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlOZ4CAsKCp577rlp06a9//77DTgQAADO5HgIH3vssRkzZqxdu/bixYupqakNOBMAAE7j7tjNsrKyOnbsGB0dLYRYuHDhLbcpLy//7LPPaq536dKld+/e1RYrKiqsVqurq7qnai0Wi9ForKiokD2ITEaj0cXFRfYUMpnNZg4Do9Ho5uZmt9tlDyKNyWTiMDAajRqNxmaz1XF7jUbj5ubm2H05GMJTp06dPn167NixWq3W399/2bJlnp6e1ba5cuXKmjVrat521KhRkZGR1RYrQ6jyk6DFYjEYDAaDQfYgMlV++io/A5rNZg4Dg8Hg4uJS92fAlsdkMnEYGAwGNzc3q9Vax+3d3NycHcKrV6+2atVq3bp1Li4uW7duXb58+bPPPlttmzZt2mzdurWOO3R1ddVoNF5eXo7N0wJYLBYhhL+/v+xBJNPpdFqtVvYU0pjNZldXV8UPA7vd7u3tXfNna3WYTCZ3d3fFDwObzebj4+Ph4eGE+3LwVKSPj8/gwYMrX8CNHDkyLS2tQacCAMBJHAzhPffcs3PnzspTWDt37uzTp0+DTgUAgJM4eGq0devWDz/88Lhx4zw9PUNCQt56662GHQsAAOdwMIRCiISEhISEhAYcBQAA51P3rysAACAIIQBAcYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClucseAGjhcnNzV65cmZmZ2bZt24kTJ957772yJwIa0tGjR9etW3fhwoUuXbrMnj27S5cu9dxhcXHx+++/f+jQoTZt2owfP37kyJENMudt1OsVYVFR0cqVKxtqFKDlSUlJ6REV8+q3+eu84v92LmjwhMT58+fLHgpoMK+//nrsfaNez9R94j16yWHzXf36r1mzpj47TE1N7daj5x+2/bTBd+zygs6jZjw1c+ZMm83WUAPfkovdbnf4xsnJySaT6bPPPqv5v9LT02fNmpWenl7HXen1eo1G4+Xl5fAwzZ3FYikuLg4ODpY9iEzFxcU6nU6r1coepGEYjcbwTp0uT10pokZVLZUXiUVxuz9fPXTo0FvexGw26/X6oKAg503Z9BQVFXl7e3t6esoeRBqTyVRWVhYYGCh7kFocP348+p5f2xb8IALDq5bOHfZ+Z3TWz5mOPZXZ7faePXtm3vuCGPBw1ZLpqlgycN2bCxITExto6ltw/NToBx98cPfddx84cOA/bWA0Gn/88cea6yEhIe3bt6+2aDabr/1XTRaLxWw2q/wICCHMZrO7u7ubm5vsQRrGd999d1kTcr2CQgjvADHoN5s3bx48ePAtb2L+NyeN2CRVPgKurupewdBcDoMtW7bYYiddr6AQIrxfeYdf7dixY+rUqQ7sMDMzMzO35HoFhRAeXmLI45s3b05ISLj9bd3c3Bw+ZhwM4ZkzZ1JSUtavX3+bEObn58+ePbvmekJCwpNPPlltsbS0VKPRNP0vfOOxWCxlZWUeHh6yB5GptLTUbDYbjUbZgzSMixcvCt921Vf92ufn/6TX6295E7PZXFZWptFoGn24Jqy0tNRqtar8vWAyma5cueLu3tSv4cjLyxO+bauv+rXLycn5T0f47V24cEH41fiW8Q8tOF5Q6w5bt27t8DHjyANtsVieeeaZlStXuri43GazsLCwup8a1Wg0nBp1c3NT/JyYm5tbSzo1evfdd4sLLwiLUbjfcJbv7D97x/f+T19os9ms0WgUPwxcXV05Nerp6dn0T43GxMSI3etvWrLbRdah2NiHHTuG4+LiXC//ZDOUCl3r66tn/xkVFdWo3xSOvJA8efKkTqebN2/elClTvv/++ylTpmRkZDT4ZEBz17Nnz5H9o8XaucJsqFr68XPf41t+85vfSJ0LaBgJCQmhxcfEnuWi8loTq1lsfrlPsPv999/v2A6Dg4OnT5kkPnxEGEqrlo7v8vx2xe9+97uGmfg/cOQVYe/evdevr/opYOrUqbe8WAaAEOKTTz6ZO3fuJ3/oam9/lyi5GB2ifX/71zXfIweaIx8fn127ds2ZM+fAH98SIZEi98SYe/u9u3Vrfd7mf/fdd73nz1/5Undr+yhxpTBSV7Fi84ZevXo14Ng11fcctLe3d4PMAbRIAQEBa9eufSMnJzMzs127dl27dm0xlwIBQoiePXvu37//zJkz58+fj4iICA8Pr/02t+Xl5bVixYpFixalpqaGhoZGR0c74b3S+t7BqlWrGmQOoAVr3749rwLRUrm4uERGRkZGRjbgPoODgwcOHOjj4+OcK4bUvUAZAABBCAEAiiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaYQQAKA0QggAUJq77AGExWL55JNPvv/+e29v71GjRg0fPrz++9y/f/+OHTvKyspiYmKmT5/u6elZ/31CERkZGVu2bMnNze3WrVtycrK/v7/siZwhLS3tq6++Kigo6NWrV3JycqtWrWRPhMaye/fu3bt3GwyGvn37JiYmajQa2RPJJ/kVYXZ2dt++fWf89eMVhV3ePOM7IvmphIQEk8nk8A6tVuvs2bN/PXHGkuOat/M6PfL2pqioqJ9//rkBZ0YLtmTJkr6Dhy84ULq8uNvTnx3q3rPXvn37ZA/V6P7rv/7rVyMm/OmQeVlh5BMf7u3Rs2daWprsodDwzGbzlClThj/8xH+f0r2VGzbz9TV9+vS5cOGC7Lnkc7Hb7Y2x3/T09FmzZqWnp99+swcffHCLoYuY/FrVxxaT+J+Rb/5u4vz58x2731WrVj3y6lvij/uFh1fV0rb//nXR7m+++caxHTqNxWIpLi4ODg6WPYhMxcXFOp1Oq9VKuffU1NQBI8eLV34Ufu2rlg5/0X7LM2fPnnXaSQWz2azX64OCgpxzd0KIbdu2jZ35hHj5n8I7oGpp/6oe/3zrxIkTLi4uThvjRkVFRd7e3iqfyDGZTGVlZYGBgQ2727feeuv3yz4Tz/1DuHtULW18caz7qa1btzbsHTWIwsJCHx8fDw+P2jetN5mvCA0Gw7avt4uxL15fcvcQo5/fuHGjw/vcuHGjGDX/egWFEPHP7U89lJeXV49JoYSNGzeKwbOvV1AI0e/BHNfA1NRUeUM1uo0bN4qhT16voBBi8OxTl0qPHz8ubyg0io0bN4rRz1+voBBi7IvbU3aUl5fLG6pJkBlCvV5v8WgldL43rQaGFxQUOLzPgoICEdjxpiU3jd23TWFhocP7hCIKCgpEQMfqqwEd63NANn0FBQUiMLz6amAL/6zVdIsj3LOVVedXUlIiaaKmQmYIg4KCWrtZROG5m1bPH4mIiHB4nxEREeL8kZuWyos8yi6FhYU5vE8o4hYHj80qsv/VpUsXSRM5wy0+a3OFyDnZsj9rNd3ia12U3cpuCAkJkTRRUyEzhO7u7nPmzBGrHxNX//3zyOWfxaYFjz/+uMP7fOyxx1y//qvIPlb1sfGKWP3ow0mJPj4+9Z4XLdyMGTO80zeIY9urPrZZxIY/Duge1qdPH6lzNa5HHnlE88274vSBqo8tRrFu3pihgzp2rPHiGM3c448/Lja/Ii79VPWxQS9WP/rb3/6WC0eFvXEcOXIkJiam1s2MRuOTTz7p5hsiYsaInsN8A4NXrFhRz7v++9//HhTSVnS/T/QZ59K6zaxZs8rLy+u5Tycwm815eXmyp5CsqKjIYDBIHGDPnj0RERGi892i34MiqNOoUaMuXrzozAFMJlN+fr4z79Fut3/55ZdhHTqIyIGi7wThHzpp0qTCwkInz3CjwsLCiooKiQNIZzQaCwoKGmPP7733nl9QiOg5VMSMcfMNefzxx41GY2PcUf0VFBQ4bTbJV41Wys7OPnDggI+Pz8CBA/38/Op/72VlZYcOHSotLY2JienUqVP9d+gEXDUqZF81WsloNB4+fDg3N7d79+533XWXk+/d+VeNVjIYDGlpafn5+VFRUV27dnXyvVfDVaONdNVoJb1ef+jQoatXr/bp06dDhw6NcRcNwplXjcr/C/VCiLCwsPj4eI1G4+XlVfvWdeDj4zNkyJAG2RVU4+npOWDAANlTOJtOpxs0aJDsKeAMvr6+w4YNkz1F08KvWAMAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNEIIAFCau8O3zMrKWrp06eXLl2NjY+fNm+fh4dGAYwEA4ByOvyKcO3fu3Llz16xZ06FDh9dee60BZwIAwGkcf0U4YcKEiIgIIcSIESN27NhRcwObzVZUVFRz3dPTU6fT1dy4ksPzNHc8AoIHgUdACMGDwCMghLjzB8HV1fHXdY6HcM6cOZV/ePXVV6/9+UYXLlzo3LlzzfWZM2e+9NJL1RZLS0s1Gk3NQKrDYrHo9XrZU0im1+u1Wq2np6fsQaQxm81lZWV2u132IDKVlJQYDAaV320xm81XrlxRPITFxcVGo1Gj0dRxez8/P4efOhwPYaUPPvigTZs2AwYMqPm/wsPD09PT67gfrVar0Wi8vLzqOU/zZbFYPDw8goODZQ8ik4eHh06n02q1sgeRxmw2a7XaoKAg2YPIpNFovL29Vf55yGQy6XS6wMBA2YPI5O7u7uPj45yfh+oVwo0bN6alpS1fvryhpgEAwMkcD+GmTZtSUlLee+89FxeXBhwIAABncjCER48eTU5OHj169LRp04QQrVq1WrVqVYMOBgCAMzgYwpiYmPLy8oYdBQAA5+M3ywAAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQmrvsAZqHK1eubN++/dy5c+Hh4aNHj27VqpXsiYAmzWg0fv3112fOnAkLCxs5cmRAQIDsiZyhvLx8+/btv/zyS3h4eHx8vI+Pj+yJUCeEsHa7d++eMXNmTmAf0babuLQv9NlnP169eujQobLnApqoH3/8MTEx8YxHJ9EhWhQcCJz3zMr3Vjz00EOy52pc+/btS54xI9u/t2jbXVze1/73z67+6MPhw4fLngu1I4S1KCwsTJg8pST5Q9F7dOXKxX9tS5g85fTPPynyQy5wRyoqKiZOnJg9YqEYmFy5Unjm+6QZ409l9AsPD5c7W+MpLi6elDC5KPF9ETO2ciUnI6XyiSIoKEjubKgV7xHW4quvvirpMOBaBYUQInpMceivvvrqK3lDAU3X3r17s93bXqugEEJEDKiIefDzzz+XN1Sj27ZtW1G72GsVFEKIqHh9+L1bt26VNxTqihDW4uLFi6Jd9+qrbbvn5OTIGAdo6nJycm75LXPx4kUZ4zhJTk6OaFvjs27Xo2V/1i0GIaxFaGioyM2svnopMzQ0VMY4QFOn5rdMaGiouFTjs8491bI/6xaDENZi3Lhx/tmp4tj260v/2haQ8+PYsWP/840Add1///0dbXni4N+vL535Xnv0i8mTJ8sbqtGNGTMm8NJhcfSGd0wyUvzOHxw/fry8oVBXXCxTi4CAgI0bPk+eMePivhjRtpu49FNY8b/+vuFzf39/2aMBTZFWq920adO0adNOp64THaJFQVZQdur/rVndgq+UEUL4+flVPlFc+Pb/RNvu4vJPoYX/Wv35+sDAQNmjoXaEsHZDhgzJPHWq8u8Rdup0b3x8vLe3t+yhgKYrNjY2IyMjJSXl9OnTYWH9R45cpcIPjvfdd9/JEydSUlJ++eWX8PCB8fHx/IXj5oIQ1om3t3dCQoLsKYBmw9PTc8KECbKncDZvb+9JkybJngJ3jPcIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACjNXfYA6jpy5MjRo0dbt249aNCgkJCQ+u8wOzv7+++/NxgMsbGxd911V/13ePXq1f3791+4cCEiImLw4MHu7g1wtPzwww8nTpwICAgYPHhwQEBA/XeYlZX1ww8/WK3WuLi47t2713+HaLKOHz+elpam1WoHDhwYFhZW/x3m5eV99913paWl0dHR/fr1q/8OG0N6enp6erqPj8+gQYPatGkje5wWyu6o8+fPz507d+rUqW+++abZbK72f48cORITE1P3vZWUlJSXlzs8TPNSVFQ0btw4EdRJDEwWfce38g9atmyZ2WzOy8tzeJ+LFy/W+gaJ2IniniSXgLCkpKR6Pp7/+Mc/OnTsKHrcLwbNEp1io6Ojjx07Vp8d5ubmDh06VLTtJgbOENFj/INCVq9eXW2boqIig8FQxx1ardbnn39e4xsi7p4s+k9z9Wv36KOPmkym+gwpnclkys/Plz2FZIWFhRUVFTeuXL16dfr06S7+oeKeJBE70dM3aNGiRfW8l+XLl7fyDxJ9x4uBySK485gxY4qKiuq5z4ZiNBoLCgpKSkomTJggAsPFgIdF3wnefkFLly6VPZrzFBQUGI1G59yXgyG02WwPPPDA6dOnrVbr2rVrFy9eXG0DQngbkydPFgOmixVXxEqTWGkSfzklWrfZtm2bwyFcs2aNCO4i3jhXtcPlJSJ6zBNPPOHwhDk5Oa1a+4qnv6za4UqTmPa3Ll26VHt6uiPDhg0TQ58U7xmqdvinw67e/j/88MON29xRCN9++20R1lv8b07VDt8uEN0GL1iwwOEJmwJCaL9VCJ966ikR/YBYXlL1tX7znAiJ+Pjjjx2+i507dwqfELHkRNUOV1wRA5MnTZpU79kbRmUIExMTRf9pYkXZv58oMoVvu23btsmezkmcGUIH3yPMyMiIi4uLiIhwdXVNSkpKS0tr2NepLVhJScmmLVtF0tvC3aNqKbiLGPPH1atXO7zPVatWiYcWC792VR97eInp73z88ccWi8WxHa5fv/5Kz9EiKv760rCnzlp89+3b59gOz507t/vgITHldeHqVrUUGmUbNvejjz5ybIei8rOe/JpoFVT1sa61mPa/q1atcniHaJqsVuvq1atF0tvCw6tqybedeOjP9flaf/jhh2LMH0RIZNXH7h4i6e0vtm4rKiqq97wNo6ysbMOmzWL6O8Lds2opuLMY9xJHeGNw8F2frKysrl27XvtQo9HU3ObcuXN+fn4112fOnPniiy9WWywtLdVoNDqdzrF5mpFTp05Z/cKE1uem1bDep3d/nJ+fb7PZHNjn2bNnxX1RNy0FhF2xuZ86dSo4ONixIUVoVPXVsKiMjIw+ffo4sMP09HTRtrtwu/k4CYv+KXPF5cuXry3o9XqtVuvp6Vn99reSlZUlHr55yNBeObmXsrOzb3lANgtms7msrMxqtcoeRKbK80MeHlU/KRYWFpZZXERgx5s2Cos6u/3sjQfPHTl9+rS4b/ZNS57e1oCOR48e7dWrl2P7bEBms/nEiRNmn7ZC53vT/wjrfXr7Soc/6+aluLjYYDDU/XvZz8+vjk8dNTkYQovFUuulEx06dNi7d2/Nda1WWzN4Hh4eGo3Gy8ur5vYtj0vJRbvFdP0VoRAi70xYWFhgYKBj3QoNDf0l/6xo3/P60pVCT1tFZGTktWeTO9K5c2dxPLP6at7ZyMhJjk3YvXt3UZAl7Hbh4nJ9Nf9Mx44db9yhu7u7TqfTarV12Wf79u31+WeFb9vrSwW/BAb4t2/f3oEJmwiz2ezh4REUFFT7pi2Xm5ubt7f3tSc1Pz8/rd1UUZYvfG449vLOhIaGOnY0CiHCwsJ+zDsjut93fclqdinK7tGjh8P7bEAmkykyMtK19JLNXCE0N3w75J0JCwtrChM6gaurq4+PT92fwVxdHf9LEA6GMDg4+PDhw9c+vOXrGFdX18DAwDru0PXfHJunGWnTps3w+wbt+nKRmPjnqqWrxWL769PeWuTwI5CYmHjgjcWi55Cqc0d2u9j00qSJE+tYlJoSEhJeXtTPNOQJ0SG6aunoV22unBk+fLhjE/bo0SO2W8e0nX8To+ZXLekviV1Lkz5ffeMO7+gwSExMfGXdQvH01qofKWxWsWlBYmJisz6K1PlGuI1qD4Knp2dCQsKaTQvEjBXCxVUIIUxXxdY/Jz47w+EHKikpafOTfxCxDwkv/6qlLxcPHTygXbt2t72dk1Q+eY4aNmT7lj+JhL9WrRr04uvXkv7nFUUOD6d+Lzj21mJ5efn48eNtNpvdbj937tzs2bOrbcDFMrdx4cKFvn37iq6DxPhXxIhn3PzbP//88/W5atRiscyZM8clsKOIf06MWyA6xQ0aNKigoKA+Q65Zs8YnIFgM+o14cKGIS2jXPnTfvn312WFmZmaPHj1Er+Fiwqti2FMa35AlS5ZU2+aOLpYxGo1Tp04VIZFi9PNizB9FWO/4+PiysrL6DCkdF8vYb3WxTGFh4a9//WvRKVaMWyDin3MJCn/kkUcsFkt97uWFF15w82snRjwjxr8iug3u06fP+fPn6zd4g6m8WCY7Ozs2NlZE3ivGvyxG/t7Vv/2zzz5b+ayrAmdeLONit9sdK+iGDRs++ugjnU5nNpvfeeedDh063Ph/09PTZ82alZ6eXse96fV6dU6NCiGsVuvWrVvT09N9fX1HjBgRFRVlsViKi4vrc9IjLS1tz549FRUVcXFx8fHxLjeehHRIbm7ul19+mZ2dHRERMWnSJB8fn9pvc1tms/mLL744duxYUFBQfHx8t27dqm1QXFxc91OjlQ4ePLh//36r1dq/f7eyqKwAAARWSURBVP9hw4bVc0LpzGazXq9X/NRoUVHRjadGK9nt9pSUlLS0NE9PzyFDhsTFxdX/jo4fP75z5069Xh8TEzN+/Hg3N7fab+MUJpOprKwsMDDQZrNVPlH4+PiMGDGid+/eskdznsLCwjs6NVovjRTYO31F+O677+7Zs6eRhmkWcnJy5s+fL3sKyZYuXXrw4EHZU8iUlZX14osvyp5Cstdffz0tLU32FDKdPHly4cKFsqeQbNGiRRkZGc65r6Zyrvnbb7/96aefZE8hk16v37Rpk+wpJNuzZ8/Zs2dlTyFTQUHBli1bZE8h2a5du86fPy97CpkuXbq0bds22VNIlpKSkpub65z7aiohBABACkIIAFAaIQQAKM3xq0Zv7+TJkwkJCbGxsXXc/tChQ4GBgZ07d26MYZqFsrKygwcPjho1SvYgMqWmpoaGhla7AlkpJSUlhw4dGj58uOxBZDpw4ECXLl2ayF/pk6KgoCAjI+P++++XPYhM33zzTc+ePev+L/M8/fTTdS9ONY0VQiHEvn37FH/HGwDgHEOGDHH4Z+hGDCEAAE0f7xECAJRGCAEASiOEAAClEUIAgNIc/GeYGs/KlSt37dpV+efhw4c/+uijcudxvqysrKVLl16+fDk2NnbevHlO+p2zTdLq1atHjx5d9+unm7sLFy688cYbeXl5d99999NPP13rP/nZgqn2pb8RzwDC+SFwzq80vVM2my0pKen48eOyB5FgzJgxp0+ftlgsn3766aJFi2SPI8fVq1cXLFgQGhqqzjFgs9keeOCB06dPW63WtWvXLl68WPZEcij4pa+GZ4BrnBaCJvoj5xtvvDFy5MhevXrJHkSCCRMmRERECCFGjBixY8cO2eNIM3bs2Kbzz+I4QUZGRlxcXOWXPikp6aGHHpI9kTSqfemr4RngGqeFoIm+R3jy5MkNGzZ8+umnsgeRYM6cOZV/ePXVV6/9WTU6na5///6yp3CqrKysrl27XvtQo9FIHEYiBb/01fAMcI3TQiDnFWFGRsaiRYuqLb7yyitRUVGVf/7www/tdvv48eOnTZvm9OmcodZH4IMPPmjTps2AAQOcPprz1PogKMVisaj8piCqUeEZoFZOC4Gcb7yoqKj169fffhsXF5c7+pfKm5fbPwIbN25MS0tbvny5M0dyvrocBuoIDg4+fPjwtQ9tNpvEYSCXIs8AdeGcEDStU6Opqal/+ctfKv+s1+vNZrPceaTYtGlTSkrKsmXLXFxcZM8C54mNjd2zZ4/dbhdCnD9/3tfXV/ZEkINnAOeHoGmdirnnnnt27tw5ceJEd3d3k8m0ZMkS2RM529GjR5OTk0ePHl15KqBVq1arVq2SPZQ0np6e6lw77uXllZycPG7cOJ1OZzab33nnHdkTyaTUl/5GPAMIGSHgl24DAJTWtE6NAgDgZIQQAKA0QggAUBohBAAojRACAJRGCAEASiOEAAClEUIAgNIIIQBAaf8Pcjy3Ijb6pFcAAAAASUVORK5CYII=" }, "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": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "GP Monte Carlo object:\n", " Dim = 1\n", " Number of observations = 20\n", " Mean function:\n", " Type: MeanZero, Params: Float64[] Kernel:\n", " Type: Mat32Iso, Params: [0.0, 0.0] Likelihood:\n", " Type: PoisLik, Params: Any[] 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": { "collapsed": false }, "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", " 5.084223 seconds (26.67 M allocations: 2.109 GiB, 6.69% 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": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeWATZd448O/kPtsm6QUt930XyumFiLqvouKBCrry7rpeu+u9/t7V99VVVvdW99Zdz9X1dnXXC1BQRFFRaKHQAi2F3nfuZJJM5nh+fwRrbUubTCaZSfL9/AXpZOYpzMx35nm+z/OlCCGAEEII5SqV3A1ACCGE5ISBECGEUE7DQIgQQiinYSBECCGU0zAQIoQQymkYCBFCCOU0DIQIIYRyGgZChBBCOQ0DIUIIoZyGgRAhhFBOS2sg7Onp6e7ujnNjQgjHcSltj8LxPC8IgtytkA2eABzH5fIKiIIg8DwvdyvkxLKs3E2QUzpvgGkNhH/7298ef/zxODcWBMHtdqe0PQrn9/uj0ajcrZBNNBr1+/1yt0JOHo8nlyNBJBIJBoNyt0JOTqdT7ibIiabpcDicnmNh1yhCCKGchoEQIYRQTsNAiBBCKKdhIEQIIZTTMBAihBDKaRgIEUII5TQMhAghhHIaBkKEEEI5DQMhQgihnIaBEKHM0BWClxuF/S4S4nJ33TWEUkEjdwMQQnE55hfqfaTeR1QUVW4m0/Nhtk1l18vdLIQyHwZChDJD89frbgqEtAahNQjbO/giA8yxq+bYqCKDrI1DKJNhIEQoA/AEWodbgLovAh93Ch93QpGRWjOOmmil0t40hDIejhEilAE6aBLlRxoa7AuTel/amoNQVsFAiFAGaAqMvk1XCJNoEBIDAyFCGaApMHqQ6w5BLhfyRUg0DIQIKR0nQDs9eoSL8MQbxTFChBKGgRAhpWsNEk6Ia0vsHUVIBAyECCld83D5osPqDqeyHQhlKZw+gZDSNfnjfc/rChEA7B09qZdeeul3v/ud3K2IF8/zarVa7lbIYMaMGa+88ko6j4iBECFFiwrQGYp34+64t8xNvb29s2fPvuuuu+RuCDqphoaGX/7yl2k+KAZChBStJUD4uHNBAywJsmDRprRFma2oqGjhwoVytwKdlEolw4AdjhEipGjxDxDGdIcxXwahxGAgREjR4plBOBAmjiKUKAyECClXhE942A+HCRFKFAZChJSrJQhCgovFdGEgRChBGAgRUq5E+0UBwBOFyIjLcyOEBsFAiJBytQQTnhRICOnBafUIJQIDIUIKFeKp3oiYL2K+DEIJwUCIkEK1hVTiqklgvgxCCUk2ED733HO9vb39f21qarr99ts3bNjw8MMPR6PRJHeOUC5roUUuloZvhAglRHwgDIfD99133//93/85nc7+D2+55ZZbbrnlhRdeGDdu3G9+8xspWohQjmoNibw8+yIQZ7UKlE2qq6urq6tjf66tra2trR34CRpBUkusXXDBBYPWhF27du2UKVMA4Jxzznn//feTahpCOSzAgieqMhjEfFcg0BshY024+nZuOXbsGAAsWrQIAObOnQsAr7/+ev8naATiA6HRaFy2bNmWLVsGfnj99dfH/nD//ff3/3mgp556atBXYp544ony8vKBn/A87/F4ZFl3TiG8Xq/RaNTr9XI3RB7RaJSmaUHI0VebOp86HGYFQRB3CRzu4HQ2XvJWpVM4HGZZlmVZCfdJ07SEe5OKIAiffPJJIBCYPXu21WqNfVhcXAwABw8enDdvXv82Xq93+vTphYWFsW0aGxuPHDmyZs2akpKS3t7e3bt3A8DkyZMrKysHjlj1f93v959xxhkFBQUA0NLSUl1dbTAYzjjjDLPZnM7fd1Q8zzudzkAgoFarw+F4c6Dz8vJ0Op24I6Zk0e2nn366pKRkxYoVQ3907rnnfv/73x/6+ZQpUwzffvrleZ7juPz8/FS0MCMIgmA0Gg3iXgoyH8MwarU6Z08Alx90OtpgMIgLhEE1ZPq/nFarZVlW2hNg2Kvp6h08zUl4kJG8uEpt/vZNlxCyYcOGZcuW2e32iy+++Gc/+1ns88svvxwAHnrooVdffbV/m5KSko0bN/70pz8FgPvvv/+qq66qqKj4wQ9+8O677w460M6dO4ceorCw8Morr3z11VebmpoefPDBdevWud3udevWDftyIiOVSpWfn09RlFqtjj9IJ1OySvpA+MYbb1RVVf31r38d9qfjx48/44wz4tmPSqXSarVabe4upK/9mtwNkYcgCLn867fSvFqtVqlU4gJhX5TSajO7NyX2Xy/tCTDsvfKqqSpWSFN6kW7I/8nu3bsrKiruvPNOADCZTMN+a+A2/S89d95553XXXQcAzz//PAAUFxcvX74cACorK4d+nWXZcePGAcDkyZN37dpVWlqq1WqLiormzJmzdOlS6X4/aVAUFbv21Wp1eu4AEgfCN998c+vWrX//+98pCscnEBLJw4A3uZzrnjARCFHhZRiHNeMoGasZ9/T0TJgwIfbnk902B27Tr/9dedSbbU9PT//A01lnnbV06dLi4uIHHnhg586db775plar/dOf/iT+F8gKyT4z6vX6/ieUmpqaa665xuPxrF+//oorrrj22muTbh5CuUjEymqDsAK4RE3GR2lWUVGxffv22J+7u7sBQKvVMgwDAE6ns6amZtA2brf7ZLvSaDTDDqlWVFQEg8HLL7/88ssvN5lMPT09+/btCwQCN9100+OPP97R0SH5L5Vxkn0jvOeee/r/vGDBAmWORSOUWZqDEvTUdYehyJj8blBqTZw4ccWKFRdeeKFOp2tvb7/rrrtWrly5YcOGDz74gKbpsrKygdvo9fr29vaBd10A6B9FO/3006+99tpt27Y9++yzA9NGBh6ivLz8kUce8Xg8N910k1qtFgRh/Pjx6fx9lYkSt3SFOJs2bRIEYdOmTfFszPO8y+WKpU7lJo/Hk+PJMjRN2+12uRsig0cPCv4oCYVCopNlAOCUEurc8gweJgyFQpIny/zhD39obm7+wx/+IOE+JRSb7XD55ZcTQkKhkNKSOdOjpqZm48aNNTU1fr8/oWSZZGTwdYJQVnJGwB+V4PEU6zFlLoqicjMKyiUl0ycQQqIlP0AY0x0GQgimrWWQlStXyt2EHIWBECFlkSoQhjnii1IFOboeQ0bK5ZEgeWHXKEIKQghpCUq2t24sTIhQHDAQIqQgPWGgWcny17AMBULxwECIkII0S/c6CPhGiFB8MBAipCBNfinf4fCNMIslVGIpVpUppe3JaJgsg5BSCIS0SDGVvp8/CkGWWLSYOJqFBhZdGlWsKhM6GQyECClFVwgiUpdO6gmDJUfXLc8Mvb29gUDg0KFDeXl5p59+emz9hFGrJg0qujT0K729vRzHNTY2ejyetWvXxqoyFRcXj7zZwFYBQGNjY2Nj45o1a4xG49atWzUazQUXXKBSqYY2b9QqUen+Z00QBkKElELaAcKYrhCZkodvhCPxvfMMYZn0HCv/wh9Q2m/VzNu5c+eDDz545513NjQ0PP744y+//DIAxKomORyO9evXP/bYYz6fb+SqSUMLLe3cufPuu++++eabJ0+eDF9XZVq3bt2gPVdVVQ3cbGCrYmWe5s6de8UVV+h0ug0bNvT19f3ud7/7n//5n0HHys/PT7RKlNJgIERIKdok7ReNwXyZUWkcpYSTsvzvSIZbMO/ee++94oorAOChhx7as2cPz/P9RZeWL1/++OOPX3nllYOqJg0qujS00BIAXHfddXfcccfAAw0s5xTb89KlS4duFnPHHXfEiqu/+OKLTzzxhM1mI4SsX79+6LEcDseoVaIUDgMhQkrRQUsfCDFfZlTmU86XtwH9q/9MnDixq6uLENJfdCn2yeLFi0eumjS00NLOnTunTp06dLNBewaAoZvFxPo8AUClUtlstv52Dj3W559/nkyVKCXArFGEFCHIQiAFryXuCDA8xkJF6+npif1h27ZtixYtGlh0aevWrUuWLBm2atLAoktDCy0Ne6ChexbR2qHHirNKlJLhGyFCitCRmlc3AtAThvGWVOwbSeP555/fsWNHNBpdu3ZtrMuxv2qS3W7/y1/+4vf7h1ZNGlh0aWihpePHjw88RKzHcuBmsT1v3bp12CYNrOI0cPlvs9k89FgajSaeKlFKhmWYlAvLMOVUGaYdnWRnlzDwkyTLMPU7f7xqaVEGdE8NkiNlmPpLL8ndEKXAMkwI5a7OlA3mdaVg6BGhbIJdowgpQmfKwlVXGAOhcmHpJSXAQIiQ/LwMoblU7bwvDJxANKrM6x3NBbk8+qMc2DWKkPw6U1lNnifQF0nh/hHKdBgIEZJf6gYIY3BaPUIjwECIkPw66dTuH/NlEBoBBkKEZEYISXU+C74RIjQCDIQIyczNQDhlmTIx3aE0zhdGCUqosuCg7RP9LhoWBkKEZJbSTJmYqACuCIZChTp27FisuKCI7RP9LhoWTp9ASGapzpSJ6Q5DYY4uUjSKfze8F+XTVH3i0ulrtOpv1YccVFnwZHUEAeDgwYPz5s0bWomQELJ58+bW1tZ169b11wJECcFAiJDMUp0pE9NFk7k2nEo4DFfYw/DR9ByLwCgPPcPWEYwtwPbQQw+9+uqrQ79y3333/ehHP5o1a9YPfvCDt956KwWtzn4YCBGSE0lXJgvmy5zMdQu+K+PRB1UWhOHqCI68/Z133nnjjTcCwFNPPZXy5mYpHCNESE7OdJVJwjn1meJkBQJPpn9h+oyo/KdMGAgRklMqivEOK8BCFAsTKtLAyoKDaLVahmEAwOl01tTUjLo9EgcDIUJy6kp9ymgMIcQTxTcGJTr99NNfeeWV73//+wCg0+kG1gJcuXLlCy+8sHHjxhtvvLGsrGzk7TOi8p8y4RghQnJKT8pojCtCSowYCxWnsLDw7bffjv157dq1A39ks9m2bNkSCoUGBrmTbf/MM8+kvrHZCd8IEZKNQKA7XW+EAOBm0ncsJBWKovBVL9UwECIkm94wsEJa3wjTdiyEMggGQoRkk85+UQBwp2myHEIZBscIEZJNGhZXG8iFMyiks2vXrueff76urk4QhGnTpq1fv/7888+Xu1FIJHwjREg2aX4jDLIkgjMoksYwzPe+973vX3pd3rGS/57w42sn3zauZ8b/3nDf+eef7/V65W4dEgPfCBGSB0+gN+2rvXgYaowp3QfNMt/73vfYI6qnr35NRZ14kZhVMve8OWtf3vvs2rVrP/zwQ40G76sZBt8IEZJHd4hwacyUicF8mSRt2bLl8KcNt626uz8K9tuw+PuGXmsy65yNUFMJyy2lFAZChOSR5gHCGJxBkaSnnnrq8kXXUDD8dMzLF13z9NNPi975CDWVsNxSSuErPELy6EzX4moDuRgCJ7mJo3hUVVVtOPvGk/10evGsA68f4HlerVbHv09BED755JNAIBAIBKxWa/8nXq93+vTpscpKg0ovDdogVqcJiYZvhAjJI82ZMjGYOJqkaDSqU+tO9lMVpaJAxXFc/DskhGzYsKG6utrlcv3+978f+AlN0xs3boxVYhr2KyfbACUK3wgRkgEryFMOArtGkzR58uR2b6vdPHz9295Ad3FpkV6vj3+Hu3fvrqiouPPOOwHAZDIN+iS2juig0ktDN0BJwjdChGTQFSJpT5QBAAhxJJzA6woa7LLLLtt25L2T/XTbkc2XXXZZQjvs6emZMGFC7M+xOkoDPxn1K0gSyQbC5557rre3t/+vbW1tt9566/r16x955JGE+gcQyimy9IvGuBlMHBXvpptuOhTev6fli6E/auyrf+/YGz/96U8T2mFFRcX27dtjf+7u7h70idvtjv1hYOmlYTdAyRAfCMPh8H333fd///d/Tqcz9gkh5KabbrrtttteeumlMWPG/PrXv5aokQhlm05atkPjMGEyjEbjW2+/9eeqXz3zxWOe0IkgRDPB1/e9cN+2O1567cXS0tKEdjhx4sQVK1ZceOGFl1122WuvvabT6fo/Wbdu3bPPPhvr/BxYemnYDVAykhojvOCCCwYmR9XW1i5evHjKlCkAcNVVV11yySXJtg6hLCXnGyGuOJqcmTNnfrX3qwcffPDaFy9TRTUqShUG+pJLL/nsuV3jx48XscPrr7/++uuvH/mTgaWXht0AJUN8IDQajcuWLduyZUv/J01NTdOmTev/q1arHfqtN998s6GhYejnmzZtKikpGfgJz/N+vz+hYecs4/f7o9ForD51DmIYJhQKJZSGnikYnnR4tUCNMo2BYRiKolQqiQfy29yCz8RLu89UCIfDktdhj0SkeR0uLCz84x//+Oijj3Z1dXEcV1ZWNuztDokjCILP5wsEAmq1Ov4hNpPJJPp/QcqsUY7jRl1baNKkSatXrx76eX5+/qDfQaVSabXaXD69tF+TuyHyEAQhW3/9TgZU6tHDm+pr0h7dx6m02gzIkotFQWlPAGn/MdVqdXl5uYQ7RDEURcWu/VgUiPNbyfznShkIi4qKBi4CJAjC0G0WLlx43XXXxbM3nufD4XAsnzg3MQxjNBoNBoPcDZGHWq0mhGTlCeD2C1rt6F2j/fcCaY8eJJTJlAGBEABYlpX2BMDhtIxAUZTJZOI4Tq1Wp+cOIOX1UFlZ+dFHHxFCAKC1tTU/P1/CnSOUNbrSvtb2QBGe0CwmjiL0jWTfCPV6ff9Dlslkuuaaay688EKj0ciy7J///Oekm4dQFpIxZTTGHQVzFnY5IyRSsoHwnnvuGfjXdevWrVu3Lsl9IpTFwhx45M7bdEVgnFnmNiCkHJkxVIBQ1ugMkdjwgYxwoTWEBsJAiFBayVJ9aRCsSojQQLjoNkJpJeNU+n6uHF5lraenZ8+ePXK3Ap3U0aNH039QDIQIpZXsmTIA4GYoQgg12oz+7HPxxRe/8cYbP/rRj+RuSFzimZmdlebOnZvmI+bivzJCcqFZ4ovK/zYW5UmQo6y5lzg6ceLETz/9VO5WxKurq2vMmDFytyIn4BghQunToYABwhjMl0GoHwZChNJHCZkyMViDAqF+GAgRSp9OWv5+0Rh8I0SoHwZChNKnSzlvhDmcOIrQIBgIEUoTL0MCilnkE6cSItQPAyFCadISlLsFA3gYCiMhQjEYCBFKk1bFDBACACsQvwImciCkBBgIEUqTloCyAg/myyAUg4EQoXSgWeJilLWSCw4TIhSDgRChdGilQfaiE4O48I0QyapH1grVA+ESawilg9L6RQG7RpGsdnUL2zvItHzqgvGqfJ3MjcE3QoTSoVUBa20Pgl2jSBYCgXdahO0dBACO+shf6vhd3YIga38JBkKEUi4qQLcCqi8N4mFA3rsPykFRAV4+JlQ5vznxWAG2d5Bn60mffMv+YSBEKOVag0RQXsThCfiicjcC5ZIAS56t54/6hrkY2mjyt0PC9g7Cy3GlYCBEKOUUOEAYg8OEKG36IvB0PRlhlUGekF3dwhOHhY60z7jFQIhQyilwgDAGE0dRerQFybP1vDeOFW57wuSpI8I7LUJUSEO7TsBAiFBqcQJJ/xNunDBfBqVBnYc8d5SEuHi3JwBVTvL8cXUqG/UtOH0CodTqoIFL47NtQrBrFKXa7l7h/TYxSVkBlgJI04MaBkKEUkux/aKAb4QolQghW9vJl70ZcI5h1yhCqaXYTBkA8EaBxxkUKDUOuCEjoiBgIEQopQRC2hX8RigQ8GLvKEoBhofYlPmMgIEQoRTqCVMRWSZGxQ2HCVEqfNQpKKcM9agwECKUQi1Bpd8LcAYFklxvmOztk7sRicBAiFAKKXmAMAbzZZC0CCGb20hmjT1jIEQoVQghrYp/I8SuUSStA25oVvzz3yAYCBFKFRdD0XFPIpYLvhEiCTE82d6h1GmzJ4eBEKFUUX6/KAD4osqd748yzo5OEmDlbkTiMBAilCqtSl1ZbSAC4IljBUiERtUbJnsyKkemHwZChFKlJSh3C+KDiaNIElsyLUemHwZChFLCH4V41tpXAsyXQck74CJNmTAWMCwMhAilRAYlzmG+DEoSw5NtGZgj0w8DIUIpoeS1tgfBN0KUpI+7MjJHph8GQoRSoiWQMQ/IrgzpwkXK1BeBr3rlbkRyMBAiJL0QB06GkrsV8QqwVFTZC6IiJdvcKmRojkw/DIQISa81SEjm3BoIIZ6o3I1AmemgO4NzZPphIERIespfa3sQV0TuFqAMxPBkW+bUWhoBBkKEpNeaITMI+2G+DBLhix7ij2IgRAgNERWgK5RhdwecQYESJRCyz5Ulpw0GQoQk1hYkQqbdH9w4RogS1OgHX7acNhgIEZJYxg0QAr4RosRVO7PnnNFIu7tPP/30ueee8/v9EydOvOuuu4qLi6XdP0LKl3EDhAAQZIHhiV6dMVM+kLyCLBz1yd0I6UgcCH/xi1+8/PLLNputqqrqgQceeOyxx6TdP0IKxxPoyISiE0O5GWqMSe5GoAyxz5XxcwcHkjgQUhRlNpsBwGg0qtXqoRvs2LGD5/mhn99www12u33gJzzP0zQdDGbg07VEaJrmeZ7jFF/aNTUYhgmFQjqdTu6GJKaNBjoyzJkvQjQaValUKlWaxi9a3IJVSWOboVCI47hhbyM5QrE3QELI7g51NJra/gOKi9J0NP75uAaDQaMRGdEkDoTnnXdeaWnpjBkz6uvr9+zZM3QDQRCGvbMTMngCMvmatC3MIDn+L5Chv34rnam9iz1hmJOvrH/tTDwBJKTYX7+FpjwpjoKQ3juAlIGwt7f3gw8+qK2tHTNmzM6dOx944IF//vOfg7ZZvXr1pk2b4tkbz/PRaNRqtUrYwszCcZzRaDQYDHI3RB46nU6lUmXcCeDsEXQ6aS5djuNi/wiS7G1UfqCsVgVlz6nVapZlM+4EkFAwGFTmr1/vlOwkHwGlIhaLLtbFmGpSnve7d+++8sorx44dS1HUmWeeGQ6HJdw5QspHANoyMGU0pjvT5j4iWYQ4ctibbaeKlIHQbre3tLTE/szzfCSCqzah3NIdIpFhRsAzA81BkM22GxyS3AE3cBlTWCVeUnaNnnrqqe++++6ll16qVqsjkcjNN98s4c4RUr7WjH0djOkJg0UrdyOQslU7sy4MShsIKYr69a9/LeEOEcosbUpM8UtAd4hMycvUZB+UBm1B0puNQ14KGhtHKNO1Z+YMwn492XiPQxLKptVkBsJAiJA0aA68Gb70IgZCNIKoAIe8cjciNTAQIiSNzM0X7dcXIZyS5tQjRTnoJgyfnacHBkKEpNFBy92CpAkEnBEcI0TDy9Z+UcBAiJBUMn2AMKY7nA2/BZJcTzhTF9GNBwZChCQgENKZFRPSezAQouFU9WXhrIl+GAgRkkBfhGIydir9QJgvg4ZiBXLQI3cjUgkDIUISyI5+UQDoDsndAqQ8hzwQ5rLkDB8WBkKEJNCe+SmjMSGOBFi5G4EUJovTZGIwECIkAcnfCBnO6WcO0UwrQLrHZnD1bTSQKwKtmZ8RPTKJ6xEilIMYHpzSrTDvpvcfd74QiboM2rG8x8cJoXG2C8fZL6KoNJWo7QmTafk4iQKdUO1SaFlECWEgRChZ7bQ09wlCyLG+5/oCu6eVXOswLwmHwwaDIcJ2HO39R1/wi3ll9+g0NimOMwrMl0H9BAI1rmzOF43BrlGEkiVJvygBcqT7z4HIsSUTHym0LKWoE+9kJv24+ePudZgXV7feG+XSkbqHXaOo3xEvCebAmDEGQoSSJcmaMs19r4SjXfPL79WoB5fkpoCaWHhFaf7KA+2/EEjK1zN1MVlYcA6Jk/VpMjEYCBFKCiEk+TdCV7Cqy//xvPJ71Cr9ybaZ6LjCpCtr6H4yyWONSiDQh0W1EQDDk6aA3I1ICwyECCXFE6VCXFJ7YHl/ffdjc8bcplXnjbzljNIfesN1ruCepI4XB+wdRQDQHAQ+29NkYjAQIpSU5F8HG3v/UWw9Ld80e9Qt1SrDrNKb67uf4IUTCS20u+fwR/8+sPnFvuOHkmzGQJgvgwDguD8noiBgIEQoSUlOpfeH6z30wUlFG+LcPt8022ZZ0Ox8HQC+ePH3f7t0btffH6dfeu1fPzj3X/97FRuWZsIXLr2NACBH+kUBp08glKT2JEIPAXK099kpxVerVYb4vzWl8Oqvmm6nD5B9Tz787JWvlOaNBQCGY36+5adbHr7jovueEN+gr+HS2yjIktwZKsY3QoTE44SkYoYruEcQosXWlQl9S6exldnO3/3GI9edcnMsCgKAXqO/a/XPare8zEUluHuFOfCnPDsVKdrxAGT9PPp+GAgREq8zRESX7CZAmvpemVy0oX/KYPzG2S8MdHdOKZw+8EObyZ6ntQSd3SIb9G3YO5rjcmeAEDAQIpSMZDJlXMEqALBbFov4rkZlsjrKOn3tAz8MsyEf4zMVFIpu0kCYL5PjmoNytyCNMBAiJF4yU+lbXW9OcFxGgchVPRedf+PzXz0RYb+JV09+9uepp52nM1nEt2kAnEGRy1wR8DI5dAJgsgxC4ol+I/SHj0Y5T3HeCtGHrrz4poaaN65+fu05M84zac27m3e5CtQb7vmP6B0Ogvkyuex4ILfWFsJAiJBIAZb4xGaUtLvfLrOdn0yXDKVSXXLvix/vuL67azIbjsy57L7pp19AqSTr43EzFCuAFvuMclLuTJyIUe5pLhDyi31/oFksmI0USvTEiSjncdM1YwtWJ9kAo7a0fPaKKRcvXnXTAzNWXiRhFAQAgZBefCnMSQSgOSD/f/2R3ieP+1rScyzlBkIVRQU5y8N7P8yduSwos4ieSt/p21aUd4paZUq+DWW28zo8W5Lfz7AwXyY3ddEkyVUDk8ewvX30V2PNpek5nHIDIQBMt6/6suP9x+r45xuEOg8R5H9GQegbHaJ6Kwghnd4Py2zfkaQNNlMFx4cC4WOS7G0QDIS56bgCXgc7fR+WWE4zaE66Br20FB0IS80zCOF94aPHA+T148LvD/IfdxFa7kcVhABAIKRTVKaMJ1Sj01gt+kmSNIOiqDH5q7t82yXZ2yCYOJqbFDBAKHR5Pxqbn+zYQfwUHQgpoMYUnNXp2xb7a4CFjzuFRw/wrx8XlPDMgnJZTxiiohLrunwfjs0/W8KWjFzcDd4AACAASURBVMlf1ev/jBcYCfcZ0xMmubO2CIrhCbRJs2CteG56v15rs+jK03ZERQdCACjJO9Pp3z2wGClPoM5Dnm8QnjgsRHgZm4ZyWoeo10GOD7rpfcXW0yRsiU7rsBqnOYNfSrjPmAgPotNiUYZqC5Ko6NWSJNLl21Gaf1Y6j6j0QKjX2K2mqX3Br4b+qDNEtrbl1mQXpBziUkZ7A5/ZzRUatTRz3vuNyV/V7ftY2n3G9GCqWo45Lne/KMfTbnpfiaQPi6NSeiAEgBLryp6TXOT7XeSIF7tukAzETaXv9n1cknem1G2BQuvSQLgxynkk3zMOE+Ya2ZcY7Q18bjctkPxhcWQZEAiLrct84SMs7xv2p++1kjCmz6D0ivDgSvxViWF7Q9Euh3mh5O1RUTqHdXGPf5fke8bE0ZzC8KRT7kefHv8nJXmnp/mgGRAIVSqDw1zZF/hi2J8GWLK1HTtIUVq1B8XkkPQEdhVbV1CUWvoGAZTkndHr/1Ty3WIgzCnNQZB3lhrDuelIi8NcmebjZkAgBIDivNN7fCd92q1xkcMe7MBB6SOuX7TH/2lxyh517eb5YbYvHO2SdrfuCBGXHIsykQL6RXcVWpdSKm2aj5sZgdBurghGWxnOfbINNrdhBylKHxFT6UPRjigfyDfOSkFzYlTFeSv6gsN3nIhGAHChtdwhfyD0f16cd2r6j5sZgVBFaYosS07WOwoAAZZswQxSlBaEEBHVl3oDnxVbVoiowRu/Yuspvf7PJN8t9o7miCBLnEwKz89RRdi+cLTbZp6f/kNnRiAEgCLrKb2Bz0fY4ICbHMYMUpR6boYKcQmfaX3+z4vzTklFe/rlG+cwnEfy3lFMHM0RxwMg7/oJfYHPi6xLKUjJIPrIMiYQ2k3zaaZ15ATxd1sJzeJFi1JLxABhONoV5f15hpmpaE8/iqKKLMv6grul3S2+EeYI2ftF+wK7i6zLZTl0xgRCSqV1mCv7Rlw+g2bJ1nYMhCi1RATCvuAXRdblKe0XjSnKWzHCCII4PWGZXxRQejQH5Tx6lPOEmI4Ckwz9opBBgRAAiqzLnYFR1pE66CaHMIMUpZKYQJiuR90C45xQtDvKuiTcJ8MTb1TOoSOUBq4IeBk575zO4JcO6yIVJU+t+EwKhA5zhS/cwPGjJCq814YVKlCqsALpDScWFRjOHYp2FxjnpKhJA1GUutCyeOSOExF6MHE02x0PyJxs2Bf4stAiT78oZFYgVKkMNtMcF1018mY0ZpCilOmkgU+wn9AZ/MphrkzRPPqhiixLh12bNxmYL5P15C29xAshX7jBYa6QqwESB0Kn03nXXXetX7/+iSeekHbPMUWWZc44LvJaN6nDDlKUAiL6RZ3Br4qsS1PRmGHZzRX+8FFeEFU1+CRw6e3sRgCaZS1s56SrbabZKpVBrgZIHAhvvPHGjRs3vvjiix0dHbt3S5y9BgAO62I3vZ8I7Khbbm4VMIMUSS7RqfScEPKH6u0pWF/0ZFQqQ4FpjpOulnCfPVJGVaQ4XTQJyTqc5PJ/6bCk72FxKClHJpuamsaPHz9//nwA2LRp07Db1NTU/OMf/xj6+Zo1a6xW68BPeJ5nWZbjvvX/Q4HZpC1zBg/YTAtGboyPgx1t5OyxCf0GysIwTBryDBWLYRiGYSIRZb2MNPsoLpFbhjO4N88wkwgaTkj4TsNxHMdxKlXCT6t2U2Wf/yuHUbIRlz6e8tG8Xp3WszESibAsq9fr03lQRUnb+V/vTuyslhYhvJuumej470F3e2DZSIRTq+MdU9BqtfFvPIiUgfDIkSONjY0XXHCBwWCw2Wx/+ctfhp7ER48efe+994Z+d/ny5RrNtxrD83w0Gh16+ygwLnIGv7LqRk89+KqHLM5j9TLMzpRGOBwGyN3MdYZhwuGwwSBbb8lQvii4Q7qEvtIX2FtgXDj4Co+P6ECYp1/Q5HyBZRkJBybbfGyZMa2nYiwQ6nSJ/YNnk0gkErsJpFq9S8NxsuWL+MJ1em2JCiyDLxOOi0Qi8Z//arVaEYEwFApZLJaXXnqJoqh33nnnr3/965133jlom3Xr1p3sZXEQnufNZr+BHXwfLFWdcqDtQYPhh/HspEUwLyvM4Jcqo9GoqEiQTgzD6HQ6m80md0O+0eYiBkNCeViCL1IzvfS7eq2Y/0RBEAwGg4hAaIBSo66UIU0FxrkijjussNZos6X1XqnX61mWzc/PT+dBFSUSiaTh/OcJuFsEg0G2B+52f01R3pKhNzpKQxUUGMxmcxraIOWZbbVaTz/99Fhv3rnnnltVNUp6pzhm3TgK1DTTGs/Gu3uJkKtvVEhybQkuMeoLN+i1dr22ODXNGUmhpdIZ3CvhDnsSnDSCMkVrkER5eWcQVhValsjYAJA2EC5fvvyDDz6IdeV98MEHFRWpyoW1WyqdwT3xbOlhSL03Ra1AOac1mNj9whXc6zDLc4U7LEtcQSmfRDtFVZ5CyifvxIlQtEMQomb9RDkbIW0gzMvL++53v3vhhRdedtllmzdvvvXWWyXc+UCF5kpX3Elxu/vwAkYSCHMk0YJEruDeQku6S4zGWAyTOYGOsN1S7bA7DBFeqp0hBZF3iVEXXeWwVFIgc3+DxOvZrFu3bt26ddLuc6gC85xg58OsENCqrKNu3BIg7TQpN2PHDkpKGw0J3TAYzslwnjzjtFQ1aEQUUA5zpTNYVW5bI8kOBUJag2R6Pl5HWYXhSaesqyW4AlXldmlO0WRk0soy/VSUrsA0xxOsiXP7L3vxpRAlK/F+0WqHZaGMl5jDssgVlHI2YbOsfWgoFZqDIMh3dxSESCDSaJNpoe2BMjIQAkChZVH8vaN1HuKLprQ5KPu1Jrg2v5uuTuc8+qFs5gX+8BFeYKTaobyLj6BUaPTJeXR36IDVOE0t34Iy/TI1ENrNi9x0dZxz7AQCX/Xi6qNIPE5IrAdJIJybPmgzybZ2IgBoVCaLfpIvXCfVDnGYMPs0+uW8MbqC1Q7LIhkb0C9TA6FBW6xR5wWZ43FuX+UkUQyFSKzOEHCJnD++0GGzvlynkXkOnMO6yCld76hASAu+FGYRZwQ8kvUXiOGmqx1mDITJcZgXjVqJol+Eh31OjIRIpJYEBwjdoX0OWftFY+zmSreki462yFq7FUnrqE/OWyLNtAFQJl25jG3ol9GBcKE7uC/+7Xf3Ak6uR+K0JRgAXMEqu1meiRMDmfXjeRINR7uk2mFzgg8ESMka/XIeXTmvg5DRgTDfNItmWket09vPw5AGWUeGUYYihLQlMp08yrqinDfPOCV1TYoTBZTDVOGmE3heHFl3CMIcxsJsEOWJvO/3bnq/3SJ/r0lMBgdCFaXLN87yhA7E/5UvevAaRgnrjUA4kUWzXaH9NvMChVxcdvNCF71fqr0JhLQmuM4cUqamAHDyzZzgBcYfbrCZ5snVgEEUca2KZrdUuBJ52m0Jkg5cKQolKNEZhG5aEQOEMTbLfG+oTiCSVdnBSRTZQd5+UV+4zmKYrFYZ5WzEABkeCBMcJgSAL3HFNZSgBGcQCm76gM08Sr3MtNGqrGZ9uT98RKod4rT67CDzxAl6n90s5+SiQTI7EJp0ZRSlCjFt8X+lzg2+KMZClICE3gj94WMGjV2vsaeuPYlymBe6EnxeHEF3iIRwmDDDKWDiRI28y00MktmBEADs5go3He9aawDAE7IHXwpR3PxRSGhZIjddbVfGHOF+NvNCCfNlSOKL7CClkXfiBMM5Wc5v0U+SsQ2DZEcgTCwXoMoJOLkexSnhGYT0fkX1+QBAvnEawzpZXrKCZDibMNPJPXFiv91cEatcqxAZHwhtpvm+8GEisPF/JcyR/Ti5HsUnobcfTgjRTGu+cWbKmiOOqsA8J6GOk5E1BfDyyWBRniT6eCctV3CfXTGD6DEZHwg1arNJP86bYC7Al32JldRBOSuhAUIPfTDPNENF6VLXHnESHUEYWU8IcJgwczUFElsvUGqClz6otF6TjA+EAGA3L/SEErvIXRFS78UrGY0iwkNCxXjd9H67KVWPumqeVRGRNzC7eaGb3k8kevzDYcKMJm+/qD/cqNPadBqbnI0YQuLCvLKwmxc0dD8xuei7CX3ry14ys0BBndRIgdqCiS3K5wnVlNvOS0lTCDl3x/1jug9GjDbaaA8b7UFzUchoC5kKQ0Y7bS4MGW2M7qR1qg3aYrXKSEdaLIaJkjSnOQAzCyTZE0q3RllL0rtD+xWVLxqTDYEwzzgtwvVGOV9Ci/03B4g/CnmK68RCCpLQe0+Y7RYExqQfl4qWTDv+oT4aem7Dm3rGbw67jGGPhe4zhT2lPQfNYbeJdpoibkZn/dfavwuUetg92M0L3KH9kgXCIAHA58jM44yAh5E1ENI1Ex2Xy9iAYWVDIKRAbTPN9YRqSvLOiP9bBKDOQ1aU4MWMTqo1kXWI3PR+m3kBlYLwYGD8S/f94/2zNvFqbcjkCJkcw2625oO7J7XsOjZx5bA/tZkWdHq3jrdfLEmTesIQ4sCUDfeP3CLvxAlOCAUjTQWm2TK2YVjZMEYIADbTAk/iuQAH3ThMiE6KExJbkM9D19hSM0C4Ys/fj05e7bSPsor3gTnr5te9DifpzbWb5/nC9QJJZFLkyRGsTZiZ5B0g9Ibq8gzTFJhNliWB0G5e4EowXwYAOkPEFUlFc1A2SLAYr+AJ1drN8yVvRnlndUnf4X3zN4y6ZdvYSpUgjO0Zfhl6tcpk0Y/3SbjWGpZkyjSyT5zw0PvtFmXli8ZkSSA06saqQE0nstZaTK0HZ0Sh4SU0QOgPNxo0Dslz4TQcc8pXj3267GZWYxh9a4o6MPvS+YfeONnPbeaFia4+MQJcfTvjNAflnTgBbvqA0iZOxGRJIAQAmxl7R5GUEppB6Kb321JwhS8+8GJP0ayOsfGu2XZs0pkFnlaH5/iwP7WbK9xByQJhb4QKSVbTAqXDUZ+ct7sI28fyAbNuooxtOJnsCYQOc0WiswkBwBmBnnAqmoMyW6LFeN30fskXyyj0Nk9r+viryh/EszEjgIshgkpzeOaauYf+M+w2ecapsfxqSZpHiMz9bChRx2QdIPSEamzm+YpaWa1f9gRCm2meN3RIRN21g27sHUWD9UWo+Ivx8kI4yDTnG6XMhVMR4ezqx79c9L2wYfj5el6G1PvIzm7yWhP5Ux15tFZ4up5wAhyavmZcxx4L3Tv0KxSoC4yzPeGDUjUSSzJlEGcE3HJPnHCkbLmJJGVPINSorUb9GH+4IdEvHnQTkti0aZT9WoIJPB55Q7V5xmlqlV7CBsw59O+oxnh04qr+T/ws2e8iWzvIc0eF3xwQnmmEKifhBJhTAFdPhbvnqcrNUOclrNbUOHn1nPp3h92tiEXqR4DDhBlE3okThBAPfdBmwUCYenaTmN5RXxTa6VQ0B2WwtkROCRddI22/qDXYO//QGx9W3kS+7kdiBXj+KDQGIE8Lp5eqbpmtunMOddUU1eqx1Bwb5dBTFAWVhVSVkwDAwdkXT2vcposOk+1jMy/wBCVbdLQ3AjQOE2YIeSdOBJkmrdqq1xTK2YiTy65AaK4QV3et1oMPtuhbEkoZ9YQOSDuDcMWex2tnX+Izl/Z/sr2DlJlg3UTqlGJqsnX4mezT8yk/S3WHgTYVtpdVzjz6/tBtTLoySLCW9QhwNmGmkH/iRGi/TWEVJwbKqkCYZ5wZYjpYIeGBi1o3EfByRl/zR4k37tEUhnOynM+inyzV0ac27bDQvQdnXdL/SVMAjvjgvPJRsgwogAo77HMRADgw5/I5R95SCcO8r9lN892h4ecaioC9oxlBARMn9jswEKaHitLkGWd66dpEv0hzeD2jbyRUeNYTrLGbJMuF00cDy6qe/mzZLYLqxEsfI8BbreTC8ZQhjvXMFhXCQTdhBXAXTPAWTJja9PHQbezmBRKWZMJp9RlB3okTvMD4w40FprkytmFkWRUIAcBunu+hxTztYu8o6pfwDELpFstYvvepYxNX9hR9U9p3azuZlgdT8+L6ep6WGmeBOk/spfCyeYfeHLrimt28wBuqI4SXpMF9EQqHCZVP3okTvvAhi2GSWmWUsxEjyr5AuNAlKinukIdw2D2KACCRQEiAeMIHbSZpVlYb01Nb2nOguuKbgmINfmgOwDllCbxuVjqoKhcBgI7SCl6lKe+sHrSBRm016cb4w/WStBmHCZXPGSHyTpzw0AdSV6dTEtkWCE36cQJhwmx3ol+M8HAMJ0UhgAhP+iLxBh6aadaoLAZtUfLHVfPsabv/9MXSH0U1Jx6cwzy1uR0unkDpErlMp+VTAZbqDgMAHDzJimt2c4U78fzqk2nGIr3KJm++KAC46X3KXFmtX7YFQgoo0TOlanG5NQTQFgQh7nmlEi4oM6PxfV9eWWvZkv5PtvVo59pggiWx/VAAFQ6odhEAaJpwmjXYVehqHLSNzbQAZxPmDnkHCFneG+FcVuMotVPklW2BEJIYJjziJVFcZCbnJTRxQqolRtU8u6D29X3zr+r/pNZDnAx1ZomYvS1yQK2bsAIIKs2hmRcNfSnMN80KMe0cL82rXF8EgizGQoWSfeKEK1hjM82jYPh60QqRhYHQZlrgCR0ESDimsQI0yProhJQg/mK8vMD4w0dtUuTCzTy6xemY2ueYFvtrgCVbO8j5Y1iNqAs0T0uNs1CxlJkj084b030gL9A1cINYfrUnlHB+9bAIIQnl2aJ0kn3iRGyJUTlbEIcsDIQ6jc2gcfjDg7uD4nHQhYEwpyVUjFeqXDg1H51f98a+eev7P3mnlSxxUKUG8TewSge110UAgNUYGqaeO+fI24M2wLXWcoS8/aIEiJve71D2ACFkZSAEALvYumuNfhLm8JLOXV3hBB6f3XSNJLlwsxq29BVOc379OljtJAEWTi9NamLitDwIfp0yUzvroqnHdxiYb6VM2C2SBkJ8I1QkQoi8gTDEtKkpvUFbOvqmssrWQCjyIucJHPZK3hyUMZr8Cdw1PPT+5HPh1Hx03qE39827OvZXL0M+7CKXTFSpkpugT1Gw8OuUmbDB1jJuxcyjmwduYNKVE8KGo10n2UFinBEIspLsCUmpIwTeqJwNUH6+aEx2BsJ80yyaaeGEkIjvYu5oLmuIu1RflPNEOHeecWqSR5zd8F5v4QyXfRIAEAJvtZLTSqjiOMrRj6rSQR30nMj/OjD3sjlH3tXw39wUY/nVIhapHxYhBMfXFajWLXMD3CGJ16NPkewMhCpKl2ec7hWVC9AUIAF8ts1JNAfxDxC66f0207wkryANx8yre3PfghOvg186CU9gebE0q7VZtDDBfCJlxptX7rRPmdSya+AGdnOFS7q11upwbSaFIQCHvHL+pwgk6g/VKz9TBrI1EMKJ3lExFzkBOOTBWRS5qMGbQF1KT6gm+T6fWfXvdhfPdhdMAABvlHzaTS6ZIGUB7/7CTABwdMo5U5t2DPypzTzfG6olIM1aa00BgmutKUpzgPijcgZCX/iIWT9erTLJ2IY4ZXEgFJkvAwC1HmnbgjJD/J17BIgrWJPkavpaLjL/8L/3z9sQ++u2DlheRNn0EsZBmGqFIEd1hQAAWsqXFjqPmkOubxqgzjdqS0TUsh6WQOCIrO8faJA6ue9jbikG0dMjjgXtM5NJP04QImG225h4wlI7DR6GSHtLQgrHCeR43Gvs0UyzVm3Ra4sTOgQh5PBHbx7/cjsT9JfOrPjv2aVdJfPctokA0E5DOw0XT5D4lKMoWOSAahdZY6J4ta55/IpJLZ/UDijwFOs4yTfOkuRwdW5SWYhXjSIIBA6nq2eLjzLV/3m6vfZLlVozvuLUBRdsVKk1AOCm980o+ZGIHXLRSPW/n+o48MWxYt2qVas2btyoVqd2Pr70b4Rut/vJJ5+UfLeJOpELIOqlkBAi+8MUSrPmADB8/AOECefC8VHmpdsuOPDoprP8Y9dpKrUf7r7rrh9vNlcAAAHY0i6cUwbaFHTQLHJQtV8vmXRs0qqpx7/VO2o3V7jpwatyi9YcxIL1SnHMn6ae6kBvx+PrFwbe+M9FwpzzmKldzz/71H+vCPvdUc4TYV15ia+s5utufeyKBaH/vLdWtbCsZ9pf7v/7ihUrvN7UZvNLf+Xddttt27dvl3y3IiQzZRirMuWahJIeXcF99gRLL+3519/MLX1PbHj5onnrVs8473/PfWjDktte/PPPAeCAm6gomGNLybtULGUmdj53lswzRIM2b0v/T/OMM2mmg+OlWW9eIOQwXjjKkLbcpW1/unt14bJfXvTH78y68Pw5Fz966RPzqLKP//5zN73fbp4vIsR88Pv/d17pyocuePTcmWvWzLnkj+uesvrsDz30UCoa30/irtGnn356yZIln3322ck2aG9vH/an8+fPNxi+lTPO8zzHcYIg/u0+3zivvufvPM9SVMKv1Z1B6AzwRVJksYvGsqxGo0l1n4BisV9Lz+EOu4kgxBWKeCESjBzP189K6OSs/+Td6xZdo6K+uS9cPP/Kx//6qLeve0dv8SXjgQiDy4DxPJ/M+d9vkR129kCFjQBA4/jTJx/fsadi49c/VBUYZzuD+4utpyZ/IAA44BQWFEiyp3SfAAok+tfnBFLngjjP5yTVf/Lupqv/PfCT9Ys23vnR/xt/1fgC84JET2AiCA2fvver720Z+OGVlRsffvuBX/3qVyN/V61Wq1QiX+2kDITHjh3bunXra6+9NkIg3LJly/79w7ylPf300+Xl5QM/4Xk+HGYiQjKxSKdXlzj9B636maNvO8TnbfyZxXJ29Pj9fpZlGYaRsQ0yikajNE2n5zmgj6F6Aro4N/aGq026ydEoAETiP0Q44LFPcQz8RK1S5xttnx53lxU6itRsZMjOotEoRVGir+1+ZVqgWX2Lly0xCLVlp6797Be7ZlxOvs5OternOANVedrKJI8SUx+BDhtjkeK+Eg6HczkKAkAgEDCZxKRcHg2ofCGt5O0ZSuA5LhIqMNoGfmg3F0aCPk+opizv8sjQ03pEHBMGjssftENTodfr9flGmeSbl5en08V7FQ8iWSDkOO72229/8sknqRHTv6+//vpNmzbFs0Oe5609fhObVOptUd5imjtUYlsk4rvHWbjYrtIkucJHEtRqtdFoHPSinDsYhjEYDHa7PQ3HOtxNTKZ4H13bA4eLrIsSvUMVT55d21Uzb+zC/k/6gj1eLhjUTrxhnMakG/62ZTAYkg+EALCokNQF9ZPsVMQ0ParPmxxs6io5sVb4GO3y6pZ3jCYjBdKc6n1qy0QpUmZCoRDLsvn5+cnvKkOxLFtYWCjiix/7BZMpTV2jBeWTD3UfXFD2zYNUbVeNbcIknTq/wDou4d2ZTOaSsiM9dXPGfDP7sK6rZvbs2eL+KeIk2Rjh4cOHjUbjrbfeesUVV3zxxRdXXHFFba00a9snw25e6AqKzAUIcQmsM4IyWkIDhG56n8OS8PvTsitv/ueep2s6qk7sJOT6+Za788+9YckYY74u5Q9biwqpOt+JlJnGSWdOGTCh0KAtVan0NNMq1bFwZr280lxFZ8VVt/1u+8/bvSfOn+POxj9//Nvpa8+wW8S8fsR2+NvtD3T5OmJ/beyr/8snD99+++3SNPckJHsjnDdv3muvvRb785VXXvnqq69Ktedk5BmnMayT5b1atZiBi30uMjs1KQxIOWiWdNDxbhyOdvEkatIn/KhbNnfphb96/r5f32xjdSaducnfMu2SmyJn3X+KROvIjMyigUkW6qCbVBZSxyadedk7P9695EZOfaIfyWFe6KarLfoJkhyrJUACLFjT0TOHhlHvS2td1cXrbuKizPVPbSw3lgpE6OLcq+54iCxocJgXjv7l4SxdfwvPsd9/5urxpjEc4QOU55G//u7iiy+WttmDpGQeodlsTsVuRaBAbTPPcwX3leavEvH1Rh/xRSFfZLczygxH/QmUpHfR1Q7TQnG9iNNOO/+yBzf5jtd9Nu3iCybPeqnLutpB6dKVC7WkkHq/Q6gspEImh9M+ZVzHnqbxJxJk7OZFre7/jLdfMvIe4kQADnvJ0iJ8gpRH+ldLXn7VbZWXXt93/LBKrS6cNBPU3BfHbsw3iZycSlHUKdf8ZMm6m/qO1nxvhrqyslKv10vb4KFSsrLMM888k4rdimO3LHKL7R0lAPuxQmG2S6xfNLjPbhH5qKtj6XmNW7zn3lI+f3kjmxflYX4a+xsmWUEg0BoEAGicvGpg76jNNCcYOc4LYamOVYcr18uE4aExkQoqUtEaTGNnV5bOqNDoDJ7QgQLTbBWV1AuE1mgum125YMGCNERByOIl1vrZzQtdoRoRBetjqp0k/tcFlHF4Asf8o28WI5CoL3xYdA3C+XVvtpUt81nLOALbu+C/ylVSrisah8pCao+TAEDLhFPHdh/UR09MH1SpDHnG6W76gFQHag3KvMplzjrsJfLWowcAJ11tF9svKpfsD4R6jd2gsYsrWA8Avihpkma2MVKi5gCJf0EZL11nMUzUqMX0/Fvo3pkNm6sXbACA3b1krBEmWETsJikLHVRswZGoxtg2tnJy86f9P7JbFrrpfVIdiADW9ZSH7FXkCBBPcJ9DbKaMXLI/EAKAw7LYRVeJ/vo+7B3NXvWJrBPtoqvtZrG5cHv+fnD2JUFzMc3CF71k9Vhxu0mKTgWzbVSsWu+xyasGFqNwmCslXGsNMHdUDiEOZH9qp5lmtcqo/JL0g+REILSbFyUTCA97CM3iVZ2djsbdLwoALrrKYVks4ijlHVUFvrbYatcfdZGFdsou05Luy4qovX0gEGgfuyjP32UNdsc+N+nKANQSTqJoo8Ena230HHTIQ3i5x3GcwT0Oc1yTi/TRwJie2rxAp4pIUwgsGVlbfWKgfNOMCNMb5Tw6jW30rYfgCRxwkxUlmAWXbXrD4GHivXGEo508iZr14xM9ikrglu99YveSG3i1ticMDT5yuSM+KAAAIABJREFU82zZHkCLDFCgJw1+mJmvPj7x9ClNH++ftz72o0LLIhe9V8QvOCxCyGEvkarIMIqHElZIdgerJxVuONlPjWHPmJ7a0t660p4DFtrpsU0whdymsDtgLfXmlfvyyr15Zb68cb68Mkaf1pGDnAiEFKhtlvnuYHVpwWpxe6hykhUl0jYKya/el0BegStYXWheJGLixNzD//HllbWNXQwAH3QIZ46l9LIuH7ukkNrTR2bmU42TV63c9Uh/IHRYKltcb4y3XyrVgercGAjTJ8ieSAmWEccHgkxbgXH2wA8tdO+YntrSntqS3joD4+8pnt1dMrdhyp1u+2SBUgGAmo/m+zsL/O15gfay7v2z698t8Hdwaq07fwLMvT89Lc+JQAgADnOlK7hXdCB0RqCdJuVmvKqzSkMiCR1Oem9ZwXmJHsIUcs2ve+Pt8x4FgHofCbKwyCHzWTSrgHq/nbgYAo7pFCFFrqN9jmkAUGCaW9fxCMfT4rKBhmoPgZchBVjXMy1qPYLs+e0uep/NNJdSaQFAJXALal+bfmybhme7SuZ2F8+tnbXWkz8ehqRK82qd2zYxVpiznynstvg75qar5bkSCO3mRUd7nxYIp6JE/srVTgyEWYVmSUco3o15IRwIH7WX3Z3oUZZVP3Nk+vl+6xgCsL0TvlMm39q1X1NTsNBB7XHCf5VB4+QzpzTtiAVCFaUrMM120/uK806T5ECEkMNeCrtS0kMJJVRdwarY6oN2b8sZnz0cMjreP2uTN19MZ3vIaA9rzQBpCu05kSwDADpNvklX5gsdEr2HWncCefZI+RJaUMZN788zzlCrElsAvbTvUEnvoZq5lwPAfhcxa8jUvITbmQqLi+CAi7ACHJu0akrTzv5sBbul0hkUn1Y2FOaOpocvCu1xLxOYMoKL3ldkXrig7vXztt1zZNr5H5z1gLgomH65EggBwGGpdNF7RX89KijimQtJJcFKvFWJTo1SEWHFV499ufgHrMbAEdjZDWePVcrllqelJliog27it47xW8eUdZ2YQVhoWeym94lefWKojhB4405HQqIddAtE7n5Rb+iIRVWwbtsvx3bVvLXmj0emJzyOICOlXJlpUGhZ4gzsSWYP1U68pLNEQgvKEEJcdHVhghMnZja8x+jymsafBgBf9ZGxJihXyhK8AABLiqivnAQAGietmnL849iHek2hXmvzheulOgoh5FAiMzWROPI/oxNC2l49o73r2KQzt6x+MGgqkrtBicmhQGjWTxQIF2LaRO+hnSY9ki3HiOSU0IIyAaZRq7IkNEdYHw0sPPDq7iU3AgDDw+c9ZNUYMe1MnUlW4ASqjYbjk84Y17lHy544swvNS50B8R0nQx2S/R6d7ZwR0hWS82nDGuxas+3urugR39yf1M5aOzQdZhAC0BuGaifZ2kF29ZBaD+mggZa1BnOuJMsAAAVUoXWJk947PvEaOv32u8h3ymVPd0DJSqjSpCuwt9C6NKH9L6n+x7FJZ7oLJgDA571kWj5VZFDWaUMBLCmEPX1k3ERrT9GciW1fHJ18FgA4rIuPdP15Clwj1YHaaeJhwJaOlZNzlJzLqhEy8+jWypp/7pp9tlfvohzLTrZhiIN2Gtpp0h4inTRYtFS5GUqMQLNwOATeqOCNAitAgQ5seijQUQU6sGuotCXL5FAgBIBCy5Jm17+SKTdT4yKrx4KMZeuRJI4mMkDoDH45o/RH8W9f6Do6vmPPvy56HABoDvb0kRtnKvGEWeCgPu4RaJZqnHzWjMb3Y4HQapjK8XSE7ZZwlaxDHuHU0hzqfEozuTKSdCx99se/UPHsO//1cB27xxFdPGiWrZ8lR3zQQUM7DSGOlJmh3EydUqwqM4Nx8FRaCgCiAngZ8ESJlwEPCz20SsLh6pHlViC0GefWMY+wvF+rFpm9F+LIES811y5tu1Ba9YaJO+4MjgjbG+X9VsO0ePdOyCl7Ht9bsTGqNQPAJ12kwkGloQa9CAY1zC6gql0kr3zpqbv/Ygq7Q0Y7BZTDsrgv+NU420VSHajOA6dm2NqTGaM7RPoiMhxXy4a+8+HPXPZpXyy5kVCUs2XP+MJvXjAIQJWT7OgiM/KpSRY4rQQKDaO/PehUUGyEYuOJDSkufUuv5dZjGqXS2kwLnMGkhkBwDe5Ml1C/aF/gy0LLYirugknTj20nQDVMXg0AXoYc9JBTFby6ytJCaq8TOLWuefyKyc07Yx8WWpY6/V9JeJTOUAJPHighH3fJUX2Qi3xnxyZPwcTPl9xIKIrjA8Fos800L/ZTb5T8s1GocZPvTaMuGk9VOKgiQ5oLjiVMuYGQZVlvX7eEOcECxwb6OousS52BL5PZz/EA4FWd0RKaONHj/tzMTo1zYx1LL97//BdLfhjLF/ioC5YXU2atmEamR4kR8nWkwUcaJ6+e0bgNCAEAMzXJ46pn+UTWIx/NIZxQmAJNAXJEbFJuJOBlgok8En5Nx4W/s/0+T8GEXctujp3nzuBem2meitLFXgSfrCdTrNS101VKGxcfgRK7Rjs6On7yk5+8/Z+3zVprQIgsXX/zGT+4R6M3it6hv6ft/Ufvavx0c57OGiDhMf81fsYdP9YZRPaOEkJqXNQqOcrooOTFBu3j0dNwYPPDt/QerDJqn2aNupU33Fd56Q0jvxouqnmpddwyp2MqAPSG4XiAXDBeuc+aMUuKqD19ZMbUubxKw378/HPP/d3VUKtX6/dZZ55z82/nn/9dSY6yp48sLSI6dcbcGZVPILClTUwUPLpr8wd/+J9IVwcBYi6fcO7tv5uy4tw4v6vlIud8tMmXP+7zr5/2AKAv+FWxdZmHIW+1EoHAtdMpR6atq6e4QEjT9MqVK5fZVr59wycGjaEv2Pvbbff/p+3oul+9LG6H0VDgmevOXDv+O3+68VO9Rt8b6P71tp+9sWnDhl+9J7qR+1zCyjGqdNcXR1I46iPxLCjjbmt87oazbl5+2wU/fkxFqY45G37+5N1hn/v0a+852VcK/O1Tmna8cdHfYn/d3imcMYbSKT0OwuwC6oN24mLIe+YFf7731jvP+tnqc56igDrcffDnj94dDdOLL7sx+aP4orCzi5yDGdfS2dMn9IYTDoT1n7zzwX3X3/udXy6ZsAIAdjfv+sVPv3vBb16IJxZq2fB3Prrfkz/+s2U/7o+CvMB46AO05keftpBTiqlTiuMfRlAQxV2mL774ooMr/eHpdxg0BgAoshRvWvNw52cfOpuPiNvh/nefn2ucdP0pt+g1egAotpY+uObR7k93e9qPiW6kP5rAdGykKPXx9YvufumPF8+85KJ561SUCgCmFE5/6ILf7/rHb/koc7KvLK165sCcdRF9HgC0BqGPoSrtGXBLUFOwsJDa64RXP/jwqsXXnz3j/Fju36zSeff+168+eeoXUg1PfNGL03AlE+ZgZ5eYL37y1EN3rLonFgUBYPnE02454/998vQvRv2ihmPO+fjnvryyzwdEQQBo8u4PClMbg+YbZlCnlmRkFAQFBsKampr+/6QYk848Z8yCnoYD4nbY03Bg8fjlAz8x6y2zSud1H60R30qAPX045pF5GJ7E+QTT3VCzZPy3zsNxtgl2jdXTcXzY7Ut7au3e5kMzLoz9dXunsKoU1Iq7vIZX6YADbtLVcHDQrzx7zDzB5w95eiU5ikBgc5v8K4Flh486hRCX8L8kIaTn6MEl408Z+OHiCSu6R7u7ajjm3B0PBM2lu5bfSgbEuk+6yf6e3Q7Lsu9OUSkzNTpOirtSDQYDHR1cVivIBDT6xNY77qfRGULRwYNCNBNgSLe4HcY0+MQPUyO51LghzgVlNHpDiP3WaUOAhKL08OchIcuqnvpq4X/zai0A1PtIVIB5toy5L+TrqHFmitXoB/3KHM9FeUajE3npDdUSIDVuqXaWu3rDpErUco8URam1g/+XQ9HgyHdXDcecu2NT0Fz66YrBUbDex5VoqhaWnHQefaZQXCA8++yzdzS8H+G+mRrT4j5+1NM4bsEpI3xrBJOXrd5evyXKfdOjddx59Jj/uHrCSfu44rSlDetRZJi9ffHOz528dPWWurfIgIUtPjv2saaoKH/MxKEbT2neCUAdn3AGABCAjzrJ6rEZNoK8tIjiZqzeXPefgR9uO/Je8cwFeku+hAf6oJ2EOAn3l4vebyeC2BvP5GWr3/v2//LmurcmLz1podYTUdBSPCgKHnCTfS44f0ytSTdWp3WIbI1iqB944IG0HWznzp2EkFWrVo2wzbRp076o/vy5LU/bTHaGZ75o2vnrbfev+OF9ExefKe6gjgnT6w/s/OCT52wmO8OGPzu+8zcfblrx4/8Jlxwc57iQSuJRgOGBE2BqfqpueJFIRKvVajSKS2hKD57nWZY1GsVnCw/SFCCf98R7/yiZPufdf997qKEm32gLRPxbD7/9+Bd/uvDnz9jLJw/aUiVwZ3/yy8+X/ThoKQaAA27Sx8DqsRKcFSzLajSa9Iy62PSw37qwefvDx1v3WfV53ojn7YOvP1v1t1X33ldcViHhgVgBIjzMKBj9l2JZVhAEg0Gy99GMEwwGrVbroA8PeciubvHP32NnLnzpH3e7vV1mnbkn0PXy3mc3t26/9KF/GqwFQzfWcMw5H28KmYs+XXH7wCjYEoS3W8l3p6j8wX8XmGbnG2eKbs8IKIE/pRh0Ol0qdj6IEm+yr7zyyksvvfTKK6/UVrVZxs286E9vls8T/+pNUdQVv339wOYXnt72erC+q3DizEv++tbY2Yv3Nv/EFzpcYEqqBvKXvWSuHQv2ZoaEhnUDwtGVv7k++kXRnz97JRL0jZm58HsvfG4rnzJ0yzmH3/IUTOwungMAPIGdXXDJRMV1tMTjlEmOg7/4ktv9p2c//E0oyhYsOf/CW/9IrGKmmo2s2kkqHGScBa+ahHEC2d6R1Kpj9vHTbnp1/67nfvfLmscplWp8xak/fHS/Ic827ManfPVYxGD75NtR0BkhrzeTyyaqioz80Y6vFk/4bTLtUQglBkKKoq6++ur169e/WOc/zkpQyZSiqAVrrlmw5lvrCBdZT+kNfJ5kICQA77QIN85S4+KjChdgSb03ge37Ap+X2k4vu/q8FVffMcJmeiY4/9Cb7537q9hf9/SRYiOMU1K5pfjNt1M7usznbPjfaVf8f/bOOzyO4v7/n9nrvemaeu9yl5tccQEXMMZgYyCEGlKAEBKS/IBvCCQkISGhmhrANIdibGNww8a9yt2yZKv3673f7e3O748zQrZl+U46Vev1+NGj21vtzp135j0zn3bfzZse+/yWv7lRsLzxsVwcJlA8BwoM8G0LfqgAjfSaWDloBFtvTTrAkyjmPfqPq56W1npYa6xct/i1ziroC8P/GmBeIsoUgd1byWWqOCxVbxs0CBiSU9e4oBaXmd2He5/U1eiHI6Z+ygw7Qo85bgEqan9FGocsnuMq8ZSrnjm2Yk1j2vRIGe4gBQeM+DrtUB3dCQSlSjhswm6hpj1xbG7ddxymXMBJsft65V/dJUY/Lo/aXjtCBFcI7zf005fGDbrKjqzaU/YbkvWjbYKk4dN6epQcRssRAJjcB5TiHrpuDDauXSHksjQcptzhq+r9pXbqsH0k6doghsaxFVW2ek6KuBksRhdWk86IPPrshl0nR62IvNxlwPlSpIqbTXMAmKBA1S7sIvHpouUlVRsIOqwSlZlcB/riXrt04AqN9JoY2KHDof6aPJQdXlWTNTey4R8BY1jXjOUcNEODAAADZXYfUYtGhHDooxKVGV37e38dkobNPcp1NEL/cM6BYxpzza79KtG0q5428fjqisKlfq4MAMwBOGvHs4fscjAClwmjZeioGazyDIckJbtxt0o0xeI5SuNQ3O8VpPC2tpFeEy2tHlzRX5EnOfU7pM6Wk6NWdj64rR37w/jm1Av7pHbvGR5bMzz2ReEaF8LI7ijGcSj2UevEI2mFBy0xucnQdMDqO6kUTe7+NJWlWmU5X5m/JPJycys9S4P4g9HmHhuTlXDCikMUnClaVlL5FYchFXEzbJ4TfXGvSjuOqTDkNQsG2NqG+ycXgcBnmXhi9a7pf6AYP7prHjbhBjfcnkF05Igwufaro5gsDhUGtRByPRaJq51F9lVeJg5LxWdrbb5Tcbna5lY60H/1s0aIFnMAmi/N0NDt+d5jYl7uVSpWYjzp2LvHxt4dZnIA4KwdBygYnzC0l4MRpByULkQnrbhdO5ZisFJ0x9SiaUZ3HDZOumRLKx0esRVejVMW3O7tlxkDxjMO/udswU02WXrHsVoXHDLBnVnA/WGeR+OQxX102BgIYXB6jXaQWLO7pO4gN2BHmA5wpT6eLMCVBLgSH08e4Ej8XIlOO8bPvYohp3tU4ukm1z6FYHzvW+sh4ft2vCh1OIyGw4lyU2xTaaNzr0Y8o/tzMloPMqhgXcZ1ABCiYHs7LEsfPi6QZWr0ZRMuVcKZomWjKtfWz326zvwRRfsZRPztn7Yg7DfQsxIH9Yx8YAlSeKeun9bNRdXfMMOhiqJbO47ofLChmb4jk5B0Cuezek4IeRkc5vApUD6ohbBh3K07Su4DAGY4yAvY+X4HJ+jkBRw8v13kMahNlVOOvl2TPf9M0bJIpuMeoBKVNVrW0HSAIOIQt3vMTI+SEyMBUoOHIIXP2GIYRMKU2+GrLErsLmSCoMOlJ1YfmPRwxK18rwFniiFV2NumDh4S+SBm4XMOYKZNKz35UZJVJ+UVmt2HNZLuUmH0mP1GXCLHiqFTu66f2aPHbrI/hFDqahtT8dk3N7xIowvzEkcIf9YAS1KJpIsjgoyuPWrx9H5oUr8xqIWwgzCT4xZq3ELNJcd5AUdx1YZbv36oPmPWqZLlEbeFmGAzJWJunsVzVBWP/1cM8E0L/VABwRha+bWGL6etsaXBM7kPJgjGdb/0Kaz+1iFOateOBgBrEJ+04V8WDLcFzVQV2mvARTJGReHNo6rWVo6brXNs7yMhDNOwuRX/JGeky1yKm0SHm+jT1v5QQQJTMw7858Tou1yiC6VWSRo+a8BlapR7cYq9MOW1+yryNQ/3Q6v6jaHdgf1c6dFx92xY/CoALNv4y9ITqznBWMxBAACglczSO3fHq0kmP0Sfx2uEvuaYJbbzDc7dasmsbk5ghzyjz649Ou7eyMutbXiGGgmGxnwyBnIlKESjZg/UZM1Xm6qyKLXbXx8K2/vodvUuHFN8y7AnTMMePX6/idc/KggAoys+J1n8czkLOo5sasVqHpqkvHSCYnIflPFHMxlDM23EFRjaQhjBw1ceKn1ow6JXmFTwtq8fLD2xmn1Z/YpuUAgnugI1cezke/VgGwkrHAQ0unFMlUsDpMEXMsiF3aXWHHP2i+bUyXZpGgCcc2BnCEovGymGAQhgkhIOmXCYyTmXt2jU+c0JoklG196+u+PmVlrvG+k1AAANbvzWOWqXjo69zlIPSbA1FFZ/u2/qrzsKDR42Y6MfL07p4tk2uHZru50sDkWGgxBG8AhUh0of+nrhKxzSs2LD/dHLIYPgxLeTkzTe1DLSpQeeYzHWjNQ7d2sk0xEwrnSC0GPMrdt+YtQdABCm4bt2WJhC9JuTjJAFXEb/ie4YBWr3YUsAV+bfmN5yMJ03Xu/c2Xe3C9OwtpG+xsu52IJ4TR39UQ1tCVz95HjBoMgZB148XPqQh6+MHGnxwAEDXpGJWJfpQ4A0+EPtcsG4/mtfvzCohfC8i9HqhZgKjriFqv2THt54w3/4fust3zwsczRH81da8XXx7eT1Lnx0JIPUgOIm4XwsyUUxxgbnru7NYBNPrq7Mv9HHkwPAPiNOFkB6v/jIcBkwJwn9upjxeAlanIrU/ZK8holgfAI6bIYgW9SQPrOs8TxNh9yB+r67ozUAG5uvUSEkabxLh9+oomv6PbBywqmPHJKU+vSZkZcuEn/ZRN+STki7KrSrd+5WiacjdMXJ4hBlUAuhIUBsaaP/dZb+ohEft2Bn1MlBnOKkPWW/LR9/78LtT2lMlVc9X8IvoOiQ2x/PTr61FVo812ivHgwct+Dok4sCgNNfySQEQk7GlU7QGs+qLNURz3J7CB+z4PmJcWhn9zAQjE9AjxQzpmsIFgFsBpqgJH5RyPhZPjE+ATH7uPtOTECVDuwloaLw5oLabUmiMkNfLgoBoNKOy6+9zL3VDvxGFd6jH4CQSo3xbGbTnoOTfhV5GcbweQOeokIZl1Z/AvhhsqiVXNevTewXBrWVf5aKTJXxfGFo8uAGN+w1AAPhTDFkiiBLjDhXGwUa0mcGuZI5e/62f/IjzSndJQpBgLTS2TrnjjxeF3V2egaF8VeN8FDBcMg2MuSgMD4e44pc79yhlV6xPCmB6SlH3zwy/r4wgw0AW9tgqhqJupoyxwsCoRJJeFE2kvO7eNATBShRgGYnwnELfdzSV0k7+UwokqCjVjxLo2lJnjy/zfGu6Fi28h5EsPridhG+a4dkIU7kD0PL6+W0e/GOdtzoHpgZM4v0zTz40oHJj3SEn21uxRI2mqLq+st3+CtYDGE3k8Why6BeEUbgM6FQihanoMeK0fJMJGPDUTP+TwX9UR190oa7tye3a8Zsmfe3KeVvFdRu6f4uWsl1JveB+OZUdIbwlw003S+JkUbozHkHuMkYzg9TXovnuFpyxTj6/JrNAY6kMXUaANS6wBqAyX3pI5MpQg8VEDdow+Jui5IKWTBTS/ymhFiZRWSKUF+U8J2ihmNmHKahfPy9YxuOyBkas+dw3O/SmTCNP6+n/cO9ir0lgL9soP9bPWAqCABTjr7dljiuJak08vKYBbf78NK0Kz5Gevt2rWRuvzWvPxkCQtgBAtDwoEyNfppD/K6EmKRE5+z45bP0Lj32XHnUs0nTNl3/j+Kq9aUnVndzcQ4zQcLNNbsPxrfNjW68Wx/fS45wdWJKLgoABtceuWAsi+hqPwiAE3KPPfPZ4dKHAICiYWsbXpCC+shtJV2EHsgn7s4lojcEIoA8Kbo7l/hZPpES7xrRCg5KFqBTNhzgiI+N/enCNpvOuT2+t7gcZwi+bqb7J7Vm/+MKwcZm+o0qutLeT+lDuySt9ZDGVFk+7r7Iy1Yv7NbjFRldOMhECFNuq++U+mpJl4YoQ0kIO8MiIE+C7sgi7s9FFIY3z9NfNuI2b9cnu4Wab69/IUl/oqz8DXTlJ08rnddu/y7uTd1nwNWO4dmlBycmP26KcZatd3yXJJ1/pXcnnPy4PmOGTZoGAAfNWM2DrK4Vs1cggGUZxD25RHJPxUzLh/vy0NJ0QsSKpxxOVaFDZsAYarLm5oQSAt46f0gXx+t3yXkHLrcM1dHpSvjC+Ls2+rVK+oQFx+QDGHd4AUdZ+Zu7y56IlBv0hGFtI16SRsg5V3xy9K5dCcLSYRY+2MGQf9RkHDQ3ET1aRKQK4Ktm/E41fdrWxUPm58o2zX9B7NbP3fM8g+p6/zNBNMEfNnlDrfFtIcZ4fRM9ElnYb8S6HHT6q2kclvCLunxX5mhObzl4quR2AHCR+LAJz0+KQyMvZ04SKpH3VsAQQqMV6JEiNCuRYMYpsCNVCHwGrnZhADgy6eGZ5pDF/HVcrtw9Ow2ozTdMLIUkDfsN9Ktn6YNGTA6sBgIAxtMPvnw++3qTMg8AaAxfNtLjEyDnynkqMWCd/bvEK08WhzpxFsLGxsbHHnts5cqVL774YigU/xpmV4JDwCQVerQQTVejE1b8ahU+aMKX1IIgWbxts5+hGKwbvv8Tm+xi8YiAkSiZq7NvjXvzAhR80UAPfAe4BghS+EyMZdt0jm2J0nkIuh5zp5a/cXzMXRFvgq1tMEmJunQr7yXFMlSmjttl2Qw0S4t+XkBkieNzzclKdMiEAcAlSkyU32Bw7uyLCoWXQGPY2Mb0DXFjIUlDuYl+5Sy1o/3SEWmgKKzZxA/YT5VcKCi9tR3zGGi6prtHxeE9SyCmhJffLw0cAOIshI888sgjjzzyySefpKSkvPDCC/G9+FVBAAVSdG8OsSIDmfywqoq+pP4OTTB3Tvu9RZFz49YnBD7r5VdIkswzuPZRdPzDWQ0+GImy72swxt+3x5ZclKRcFs9RjbRrj/Cspj0s0leTPR8Aal1g8kPZFRzqeoOWj5akE3F3dUngwk9yiBVZhLRbj5toKJQhVwhFelNz4b2pASY0fdj7Fl4VFwnrm4bq/NEZgh3t+D8V1ObW7pwY+hmxWzf2zP92Tf89TTAB4JQVN7hgafpVHj6dY2uS7Ib+aeGAEGfX/iVLlmRlZQHAvHnztm3bdvkJPp/PbDZfflwmkxHERapM0zRN99BgruHBklRocKMvG+nrtGis4qJ3D4+7bwxbsHjb7zfNfd4tVHd+i8WUy/jFBufuROn1Pbhv95y04CQ+jr5qHf0DcW/JkCDWj09jvKkVYk1ZqbNvTxBOZCLh5U8aMxwsPbF6Z9nvKEBhCm9pxYtSEIHi7N/AZ6Jb0zEDuvigcXkA8sSQWQBHTHifAUK9uNK8RLypBf8sHzEQoVDf1mJakxNaGWL1ocUIY4wxrnHQe3X0dM1QMuLovLjcAhVd2WhiIvINxKlRAAAEHZ69758nRt3pECYCxpYAbNfhe7IRu9unOhi22L0VeZpf9bdrD8YxPf+XKEhMxFkIH3zwwcgvzzzzTMfvnVm1atV///vfy49v2bIlI+Oi8BSKotzukA/3PGJJw4CVKWh9O6fZRc1VkZ0tJgczb/QA54btT62d+az34poVCbzrmmyrpexpcIW9st7wVTWwfAEtL6rnyel0crlcDocT92YMCUKhkM/no6iotpNIGn+r49R7Y014QbfZt+Yqf+Pz+S5/b8q5z9vkeY3CdPD59lmYag6hYYS6OrHnEAA3pQTDLrqLuSGA1WoNh8NMZhw6aR4D1Cq0ppnto3o4WKSxQczi7GmnJ8tJLJlvd6zVnHzpfPHjvW/blSBJkqZpiqI21WFBMJSPc3WwAAAgAElEQVTCH+wzwjCNq93M43aGKRifxCs+n683g/vlTDn3uYclOp40A/t8NIb1rZxpCopPh7t/qlud38r5U4MBDBDXp/+qkAGLxdtl3+wSqVTa49GyT4K933vvPbVaPWXKlMvfeuKJJ5599tloLkJRlMTiEpC9mnIKAB4Uw1eN6GsD89YMgtvp+awvWSpE5LKDf9s0/x+dyxkKBBNaXZ+GoF4mGN2bW1+JXW7hQ8lENFH2bDabx+NxuXEolDgUCQaDXq9XLr968c8QDZ/V0QbAghgfFrP7EI+jUkoLL39L5NGPati+fvHrAr7AGsQVLvzzfCRgxTmQfHEqGq+8oosCQRAymSwuQggAaoAHZHh1Nd3jdJ43puG3zxNjVCwFB6nVy857Ps/0tRiVBXFp3uVEhDAytG1zoDIumqbu82Q6PcND4lNWXB7JbMCEeFUjwRgLYn2mr4zafL6kaee6Ra/xeUIA2KXHfDaenMhE0J140Dhk0e0dn/Y3Hru//UVRmKlSCeL4DXRD/B+rr7766vjx40899VTcr9wzOASszCI0fPRuDbYELhoDThcvb0maeP33z7DD/s7HU6SLW2wb+6g9zhBeN2TNHoMQDwnvV1MNPYpKbrVtTJEt7vKtycfeqyhc6uUrAGBTK56pRfGNSQCACUo0Qdmv43qKAM1L7vmnkLDRNA36tgVjALV8wWkJs+TYywTdH94sJI136+g3qqjafk/F2T22IP66iX6pgt7Rjvsov09cYIUDMw7+5+DEn/t5MgBo88IJC7457eqOxUbnbjEvh8fW9kMjB5A498N169Zt3br19ddf74skFz0GAcxNRGVKWF2LG90XvXV07E8tiuy5u57rHFOhFs/wBBu9wTjHUXRQ58Q72wdvnxlC2IPwfjVt6NGGjctfGwo7lKJJl7+VpD8tczSdLVgKAGfsOEBBadSW3ShJE6IFXdW46Wsmq4jeBGlMVqIABWftmEHwVfK53ynRqKp1cWxe99iC8Gkd/Vk97RgEwUjWAN7QRL9eSZ+04sFfM2Py0XcMquJIaiSShg3NeEEKEl5t5YoBt9i/TZXf1B9NHFDiKYSnT5/+yU9+Yrfbb7/99uXLl993331xvHjvGZeAlmcQ65rpI53jzBA6OPGXfq7sur3/6JjbIoKVJLuh1daHwVL7DPRu3WC3eQxyTH78QU3PAzRbbRuS5Ysv7wIEpiYfe/vwhJ9RDFaQgh06WJQc51JLEjYszyIYAzRZvCkNafk9/FsCwY2pxPZ2CIQhWbb4kNibfX691NUW1wZehfMOvKoK79bj8ABtrNiD8E0zvaqKPmUd4Lj4KEltO5JoOHOk9ILTxtZ2nCqEQunVHz+b5wQDsaX84j5u4MATTyEcPXq01+tdu3btF1988cUXX7z//vtxvHhcSBXCfbnouAW+bf3xCcYI7Sl7HAGecfDljrwzSbIFFnd5iOwixCJe7NbjHSPrwp7S5sWra3q+GeUP6Rz+ykRJF1m2C89v9PHkLckTAWCHDudLICmuRgoWgVZkEQNY1J5FoBVZDD6zhzKcyId8KezQYy5LJROO+19xUdnh16F//Ql/2Cml+3mn1ByA9U30a5XU8YFODRM9vIBz2pFVe6Y+FmLyAKDOBY0uuD66HfJm27oUxc193MBBwaA0PfclMja6Pw+5QviTejr4g0MiTTB3zvh/Ap9lavmbkSMsQqSRzG6xf9OnjdlvoHeNrAtjp9qBP6zBvl4U8G6xbUiULiCISx2ReAHH6LNfHpr4cwBo90G1E67TxnPdhhBako4GvLSClH310LFumJOIalzQ4oFU+dKjzEYU9mU37o5j86IkslP6RQPt6vvUHSY/rG2k36ikTg+RVWAH0w69UpM1z6AuBgBfGDa24JvTrl66BwCc/nNB0qYWTY17k2QcmKBEKzKJh4sYYxMGhRXtmhNC+MF9RsFBaxp+LAAWZrC/m/0npbVmwqmPIkdSFDfqnTvDlPuKF4oHe/R4p25IdayB5pQVf967ND3BsMXkPpwiW3j5WxNOfliXNdchTsYAm1rpeUnAjWsJ0jI1FMsGQ8eHHAma2VON5xBwQxJ820rzOOlCTsaaMZMmnXiPG3TFt4VRUmXHr1fR37TQ5x2Y7INZZbMbf15Pv1lFnbUNvSzgxec2CPzWE6NWRl5ubMGj5ZAaXTXpZstXaYql8dIIDgPypWhRKnq0mPh1MWNxKlEgQwlcWJJGPFTAyBQNcKe4FoUQABDAwhQk56Avm+iOEZVk8bfOeS699dCoyq8AgMNMUImmttj6dlEIAHv19Pcje6RR4AvDukZ6Q6f/sp7RbN2QKJ3LZFyaOTtFdyzRcOpEye0AUG7CbCKeosVAaJYWzUkaRD1uZiKRK+nhByyUIhkHHTTi9ITbKv17a9JnTil/K77Ni54QhY+b8Wf19AunqU/r6KNmOvoi3leCpPEJC37rHP1BDX3OMfQkEABSdMdGVa37fsZTkSQyJ6zYGcKzopv9uAP1nmCzVjy7Nw1AAMkCNFNL3JdH/GE04/YsolTZRV5vDR/uziXuzCaUAxcpNoi6ZT+DAG5KRQyENrT8mDAhwBFvnvPXgpot+bVbASBNcUu7YytJ9+2iEAD2Gejv2oZiX+s/Ku14VRV9xtbbbylEWo2uvSmyJZcc5/tt0w++srvsCZLF94RhnwkvTonbpo2ah+7PQ7MS4+x000sQwLIMlNDT0WdhMhw2Y4rI5bG132YmS1ytxef6Ixl3N4RpqHXiTS34pQp6VSW1ox03uGPeyXSTsFuPX6qgNzbTBt9Q7ZUyZ8uMAy99P+NJt1AFAI4Q/l6Hl6YRUZYPa7R8nqZY2uMKzEwCRivQL4sYD+QTsxNRqhBd9dHPkaBfFjGWphPCPqz6fEWu6erpCOCWNPRJPb2lHRb+YD328RVb5v5l0Xd/DLH4DekzVKIprdavM5V39XVjDhppADQ/+dqdmlwJZwhvasE1cXKLaLKtS5TOZTMlnQ8ijGfu//e53IVGVSEAbG3D4xUogRsH2WIgNF2LpmtgoHxEu4fDQLdlEu9V41DsEQASNipTwaZWfGPSyirdf7bP+PPN256wyTJ0mlF90dRYMQfAbKD3G0DIQpkiSOCCnItkbJBzEK+rYQ9j3OiBchOuccJQL6bNCzjn73zu8IQHjcp8AMAYNjTj6Wqkiq7IpTtQ5wk0liQ+0YNbS9gwUYnGK4ke2BQQwGgFypcSB4xwyIj7M936NS2EAMAkYGUW8WEtvdcAM37Iv+4SabfOeW7h9qdCbGFQdVt54+Mp8sUshrSvG3PQiDHQ149o4Q9ggHIT/b0OYsqj3Q0B0mR07Zuc+dolx0dVfkng8OmSFQDQ4IZ2HyxJjYNuafjo5jRC09NYhf5BzUM3pcHahp58w5NVqMKOWwO5AnZKTfjkrmlPzNr/r29u+HdkFTJI8JAdBUkufEYeE8nYIOeCjH1BHc0BKDdjs39o618EBkXO3fPXusxZ9RmzIkf2GzGBYJIy2ke6wfxpesKtsS4HkwVoihoVSK+++OseDgNdlwgTEtB3DXRf5LnskpExFzgE3JVFnLHBEdOP3cAuTfvuuj/NPPBiqs2klcxqtKztn8YcMuJvW4Ztbe6YsAXh82bW5tbYqkl0T6Pls2TZIhbjouWg0lpbfG7jnrLf0oigaNjcihcmX7FOd5QQCKZpiAfzB7sKRiiWoSk9KgIVCSv8rh2SFHc2W79qU+WdKVo2d89zzHAw7o2MI/4w1vnwWRveZ8BfN9Gra+hNLfTwUEEAmFr+po8rOzHqzshLgx8Om/GS1GhznDh8Z/2kUSuZG+XtEECuBN2XRzyQTxTJequCHYjZMD+x/zzqR4QQAIDPhLtz4JAZztp/7AxmRe7u6b+fu+f5MYxSk2tfgDT0T2OOmfHWtjhnnR9aUBjv1eO3q1GzN57zQU+w2eY9dUmaDBbpn7XvX/snP+wRqABgnwmred1VKI0GDR89mE/MTUJR2mMGA/OSeli8MBJWeMiSLuEXtdg2ni242SrLmX741bi3cIRoGHP2C4Wtbm/Z4xghAKBoWN9ML0hGkuiKaGLAdaYPM5R3IHT1nU0xG03XEI+VEHdkE6nCofOsd8WIEF5AzEJ3ZKKtbbi+kxN4u2bM4dKf3bj7xSzB7HrTx/3WmCMm/Folvd/MsA3qiXX8wRg3uPE75+idup7nhr4S9aYP0xW3MoiL7CTTjrzeljS+OWUyALR54ZgF39CLAvQMhGZq0YP5hHagIwVjhUBwexaR0SMv9jlaVOMGnuDONvumUNh+YNIvxK724nPr497IEbonveVQQfW322f/iWRyAQADrGvGGh6K3vnZ7NoPACpRWTfnsAhUIkd3ZROPFaM5SdFK7CDnWrcRdkbFgxWZxGcN9MpMIvmHZCL16TPZIc8vjm34YyE4/eckvL7KtX8JtiDonIwjdpQuoUvkUCwfyFwk/YCbxKes+KQFepwyrXus3uMB0pwou6jMZG79d3J749cLXwYAXxjWNuGbUpCopx1bwEIrs1CyYKiOCywC7sgm1tTRjTFmMOcwYEESfNOmmief3WBek6/91feznlqy+TcOcWpb0vg+au0Il5Bga5h2+NVtc57z8hMAAAN824IDFL4jLdrVDo1DdeZPihIfQ11Z5hBCKQIYo0CFMuAOob2O6BhZEV5EigCWpROfNdLGTuUozuUuak2btaKdqjO8A9DfiWDavHhLK/73GfqjGvq0FfemtuogBAM0uPGXDfTLFfj7dtxHKkjjcJ3xgxz1vQh+3PCRuNpLT3z0/Yw/hhlsjGFdEz1KBj2Oq5Nx4L5cYuiqYISIFmbGvkdaIEWTVbDbfpvFc9wdqPfwld/PfHLmoZfEbn1ftHOES+D7bXN3/+XgxF+YFTmRIzt02BjAKzIIRtRjfIt1nYSXe/lcX8yGaRri4UJ0Xx4xLgENPxWEESG8nEwR3JCE1tRjR6eY3OOj79QIJimdJr1ty4C0isa4wY3XN9EvnqbWN9F1ziGW5+lyHCHYpcMvVdAf1dCVdkz1pU201fY1j50oF4zrOMKgyNn7/nFs7N0OSSoA7NRjGmB2Yg97eLIAPZBHKIZF4UgWASt7tEc6SYnGKfnN5J3nDO9gjA3KwlPFy+fu+SsrHOiLdo7QAZMKzd391/O5CxrSZ0SO7DPgOie+M4tgRx3DECCNbfYt2aqfdj6YwEX35BK/KSbmJiFFPKKJBi0jQtgFxTI0TQ0f1kLn+oXl43+2wJ/TZlgdDtsGsG0hGk5b8Sd19L/PUJtb6BbPEHOr8YXhrA1/XEu/UkHt0dP9UMItSJpabRtz1fd3PjjxxPtuYWJ19nwAqHbiCjvcmt5Df7d8KfppLiGId7XCASSyLkyPXQsnK1GKfLYlQDTbtwNAZf5NZkXuzAP/7ueU3NcWGE8/+LJTnHSqeHnkwHErPmmDn2QTvFgi+WqM76YqlnCYCR1H0oTo/jwiXYQGVU29PmJECLumVIlmaWF1LW72XDiCETo36ekJXqm96snB0LG9YSg34/er6ZfP4i2tdItn4JvUJTTGJj8+bcXfNNOrKql/naHXNtL1rv77BquN76TIbuSy1B1H0tqOpLaV75vyCADYgnhjC16Wjvg9MsGOS0DLM4lexloMQlgE3Jndk3XhDA3BE/68zrTGFbQDwMGJvxB4zf1Zs/DaAuOJJz8QeM37Jz8aOXDOgXfr8Z1ZEFN+FrP7YCBkSpb96FBdKEN35XSdeWBYMiR7MItAPCYSs/s22HK0HC1LJ75sojtiKmiCyRj9ogHZZBX/7Ms7x4YzhI+Y8PvV9BtV9F49bR8EjqYeEs478I52/EE1/fdT9BtV9Pom+rgFmwPQz+tXo2tvkLSkKpZ2HBH4rNMOv7Z72hMhtjBMw5eNeLYWpcReaAkhNDcJ3ZRGDK60afGDRcDKLNSDdeHs5DSCM3dP47uBMFAM9o5ZTxef/zpJd6IvGnktw/PbF+z8k8p07vtZT1MMFgA0uGFTG74zi1Bcls+zG0jaXWt8L1/7KwJd0L3pGnRbBmIN1ye7Kwa14pfKqdlSgoGAy0QIMJcBDAKxO2m3m4QqO11lhxZPn4yvGSK4O5tYU48tQTxLgwCAZkuyUn+3q+WF2859Xluwog/u2XNMfrzTDzt1VLIAJQlAwEQiNvAZIGCBkAkCFvTdkx2kQOfDOh+0e7HOB4OhgDgAhML2WuMHo1Oe7giKQhjPPPBiZf5NkdRTm9qwkosmxF59noHg5nTUm1LvQwI2A92RBWvqoSlGP9LZaSt21z3+dcOBpdllwE/YOf0Pc/b8ff+UR5qTJ/VRU681UtqPTT/0SnXODSdH3U4jBgC0eeGrJnpFBqGJLo9aB3XG91TiMjEvDwAYCC1KReNi7xFDnUEthAo2rfoxTrOL/xsRCyapiEkqcJP4vAMq7bjFE+c8gSou3J8La+qxPww3JCMEwJdMlMpnn2pbd/0R87Exdwc4vYu+7gPavLjNCx0JpTpgEyBiIwET+Ezo3vOLxwQBE/GZwGeCgAkCFvAv/NWFPwvTYPDjiOy1e7E1cNnNBgHVhjeSZPNF3KwLrzGecvRtjNDpotsA4KgZ63z4gdyYN0W4DFjR05C7IUfPtBARrHGpjx5v+fv/6vJXZikMqqLd056YcuytwupvD41/0CFN7bsGD3sYFFl6cnV6y8Gd0/8QqTIIACY/fNZAL00joiyx1IHZfcjlr5uQ/m8A4DBgeSbqWV6Foc6gFsLoEbFQqRJKlcgbhnN2XOXATe64KaKIhe7NRV820l80wC3piEWASnz72q0ffHD4e6/9t4nJaeMeeLbw+tsHv0k5RIM1gK3Rnt7FtxcRSAJhawD1qZ9n79E5tgbD9mLFBQ8CwHjKsXdUlnNb5jyPEWr3wR4Dvi835lRqEja6M5uIMnnx8CCihZ/WQXMsdmgJNydDfr2v+q+vvG5E5/eEQ4GkgnH33Ji3YsdTLSmTjo/+iZ8rufpVhiCWpvPfv/5Uy6kDmKaTSybN+dXz6ty4ZSGXOltn7f+nR6hev/jVIPtCHTFbEH/SAAtTUHaMc/IQaa02vDMq5UkGwZGw0R3ZhPpaerA7MyRthN0gYMIEJbo7h/htCVqYSsSlgAAAsAm4PZNgM+DDOtrpDXz0s3m5damvLfrkm58fenT8k8dfevbofx6Oy40GOf4wWALY5IdBroLeYEuD+bPCxN9c2BTFeMqxdxKsNZvm/T3IEfopWNuIF6Wgy0ujdU+KAN2fF20K/+EEm4HuykETlLENF2LfmMqn37sLaT699X8b7t5yh2Lef/7x978l3EwRrFu++UXxua8JTPVRgwcKQ83p1fdMn4fzPr993Zd3fr2EW/rpg3NaTh2Iy8VzGr5fvO33dZlzd8x8ukMF3SH8cR3M0kChNNaxjq7UvZQiWyTm5mj48ED+tauCMGxWhJcjYKGJSihNgEYPOmzEta7eumkwENychnbp4d0PP8kL8J669flI/oXxqZP+efOqez6++cFidc3MXw/Xee4QgqL9Z9v/ma26l89OAgDAeOrRtxW22q1z/kKy+BhgbSNdLEMFsQwcCGCiCs1PRoOzmlI/wCLQ4lTIEhEbW7A/HFVX2vf+C8sL7/rJxAcjL68vuJHPFrz0xrPpX5w5l7tw8vF382u3HB7/4HDKPrPrrWfuHffA8nE/ibxcMuo2JoO5ZtVT9767uzeXZZPessOr5PbGzfP/YZOmdRxvcMM3LTBRCeMUMT+WDZb/EQQrVbEsS4yWZxKc2KsmDSeG24rwEhBCmSJ0Rzbxq0KiVHmRo01PrgZwnRbJ2w7NyJ7TOQuRRpyYpSyo1Ft+mOcOr+wvQwoM+Jz+VRm/RCOZCdChgnURFQSAXTqMAa6LJXZewkb35BELUohrVgU7KJChn+VHm0au7czhmdkXFTGYnD7d0VwX9LockpSt1z13dNy9ZeVvzt/1rMgzTBLQtJ05PDPnoo88I2uO7uzR3szCVZbqmzf9mmTzv174cocK+sKwvgl/04IXpaApqpgfS4un3ODcW5L8m+la4s7sa10FYRivCC8hgQuLUom5SXDWjg8ZcedI+VhRcBG+LAwcY3w+d9GWgryy8lXZjbsOTPpVR66jEfqTZssXobC9SPs4AADGZeVvyuyN2+b+hWTyAOC0DZ+xw8/yYnCfLZChm1KJayeg6qrIOOi+PLTXgPforp7dqEsBwD9MIpuTJ7Vpx5ZUrV+y+fGa7HkVhcuG/IYKQr0PEGKH/XJ7o9zWIDbVaNwtfJ91/5RHW5JKI+9igDM2vKMdl8jRLwt6Ui/MG2ypNrxxc8HTd+TJZJxeNnaYcG31bw4DxiegsQqodqAjZhyrU3iE1NFT937w7rKxd3QsCg0uXb29dkHRBJtE8c38f+Y07pq367m2pPEnS+4YVBVKhz1G936dc2dp+guIYAHGZeVvyOxN3819LsTkYYDdenzGBndmQZSx8ywCbkghxl97ruRXhUAwS4uS+cT6Jtp75TriKWOm7qnbkacu7DhysHEPSsr9pF2wNA1HYt0oBvtUyYqarLljz35+69cPVedcX1F4y9CVw5RRU3bXbV85/p6OI3vqtieVTOrek07gsyrsDXJ7vdzWmGBr4PltDkmqVZGlk6Y35i+wyTIi1SQAwBbE37biAAV3ZsccJhEhFLaf0/3t52MeWJaX15O/H6Yw/vznP/fbzfbs2YMxnj17djQnY4z9fr9AEHuo89VACCl5aIwC5UsJEoM5gGLSQ1Vm4b4t/62uK09XZLEZ7FNtx/+67UnRfE3OtEUcpgIQsskyqnOul9uapx15TeCzWuVZYVZPntlwOEwQBEEM8+3rK4ExpiiKxYo2Q4bDd7ba8PqYlGe4LDVgXHZkldzZHFFBioYNLVjvw3fnENLoikuoeOiuHCKnpzm444LP5+PxeIP2AZBz0ZgEwuSHKxULU2cXr3nv91QokCRNpejwzuptL+/9x9Jn3pIlZ61rxgCQIrggECSL35pU2pA5M9FwuuzwKm7IZZNlBAk2xpjJHEqTdVVW0cfvPs7CRKIkKUgFt1ZtfPPwazf9+T2xOrnzaQSmFLaGzOZ9JVXrph57O79mi8hrCrFFes2os/lLjk544Hzugpbkia2CFFKWRBNMAKAxlJvxhhY8Wo5uTiNEsVWPv4CM7a/WPbc4a8bKwkVx+bx9SjAYJAiCzWb3w73isJCPnmeffZam6WeffTaakymKslqtKlWfr6hsQXzEhI9bcDhq017I59m/+oWqHV95rAZ5RgFx/W84Y9PSGW+OS31WwEnpOI0bdJVUrsur21afMetUyQo/VxpTwwKBAJPJHFoDQRyhKCoUCvF4Uc0h3IH6061/LU76rZRfDBiXHX1TZmuMqKAvDJ830kImWpqGmFFoCkJonAJuSBn4zBpms1kmkw3yBwBjfMSMt7d17UjsMrZ+v+r/mo7vIQO+xOJx6sWSvPErUuVLHSG8vhkzESxJQ+KL07TyAvbiqq/z6rbVps04lr+Ukmj666PEB1tr3c43/q/55H5MUymjplz3y78oMwsBgBkOKmz1GnOl2lSlNleFWCKDqsCoLDSqi+ySrgMrPR6PUCgEAL0PvmmlhUxYlNLD+n8SNkxVhz+v+GuyKPHxib/ozQfsN1wuF4PB6Iu10OWMCOEFXCF8yISPm3tY56jCjo/p9qWxPhqf+qyEe1FpV4HPMubMZxktB6ryFp8tuiXEjHZ1OCKEUQqhJ9B0uu25fM0vFMJShHHZkdelztZt1z1LsngmP/yvEY+RwwxNVI4uAibcnD7AC8EOhoQQRmjx4A1NV6+iFQzbTrY8nSxbmCxbTGM4YMJHTPiG5C4qxwp8lpIzX2Q376nNuf5M4S2xTiIHHISxwG8RuQ0ij17maFWbK2X2Zpssw6gqMqoKDarCjviHbvB4PCRLcMgM5x1wfRIURV1ftzN8JszQojFy6pl9fxeyBU9NfZwYIj5fI0II0O9CGMEfhnIzXW4GLxnz1+ILw66WnazQmnT1M7mylEveFXn0406vSdadqCi8pSpvcZh5dSP1iBBGI4TuQP2Ztudz1Q8ms0s05srs+p28oGPb7D+TLF6tCzY00wu6Gmcvh0mgiUqYrkE85mAZJoaQEAIAhaHcRO81YP+VrYYAECDNp1r+lCRfkCK7CQD0PljfjNU8WJiCLqmWQJIkz2uZWPN1VuPumqy5zWlTTYrcyD7hoIJBkWKPXuTWi916kccodutEHqPIYwhyxC6hxiXSusRJRmWhOSE3zIhql4/C0OyBWieudtIkRkUyNEuDuLE7dgqYMFmNJioJgND/7f07j8n907TfMdCQ8RAdEUKAARLCCCQNJyz0IVNPcmaeMe412FeHuH+cm5zLvazPSh0t409/oracP1W8vDZzDtmt7XBECK8qhD7nkdOGV2/05U9vN4k9BqMy36AqrixYQjK5R8z4gBGWZ6Dkq/UjhFCRFOYkIVmM8fV9zdASwgj+MN5vhCMmuhtDQzBsOdXyZ5W4LCNhJQCEMezU4VNWPEaBJimhY+uPJEmapjkcjsBnLaj+Nll/UuJqNynz9eoSvXqUOSGHHqAxXeg1ye1NckeT3N4gszeLPAaPUO0SatwirVuodou0LpHWJdRQ0cleB+4QrnVBrRsa3VjJhRwxSmL6MxP4PXgoJWw0VY3GJSAWAZ6Q96k9z6sEyj9OeXQIqSCMCGGEARTCCBigyo5362hzjIVFje7jVbrX6sK/yJRNnKxCgsvGsQRb/dgz/9OYKlqSJ1dnzTOoiqCrzYoRIexSCPk+q8Z0VmusbPeXf5XgvN2ZzZNONqiLLYrsyMhIY9jShls8+I6sqxtUkgVofjJKFQ4uCYwwFIUwgisEe/T0SesV0xySlPN0618F3LR89S8i2X9cJD5igpNWnCVGU1VIy/9RCDv+ik16tcazGkOF1nha4jYYlIV6TYleMzmtdE8AAB36SURBVMoqz6ZRX7kUEXRYYW9Q2Brk9ia5o1FubwozOTZpmk2WaZOl26QZDklyj9epfgoMPmhw4zoXdoYgW4xyxJAlvlAUrMNGGD0yDpqkggkJRMQcbvSa/rD7LxM0o381/n7Ut9V64s+IEAIMAiGMgDGuccI+A27zxvBFuf31p9v/4ScWnPQsKZIRU1Vw+WqDF7BnN+zKrd/BoEK1mXNqs+Z4BBd92BEh7BBCBhXSGM+m6E4k649z/Q6dumijKlzFbC1J+T8eN73zXwUp+LKRRghuzSA43Y6NCi66LhEVSmHQJokdukIYweTH29txrbPrjkPTgSr9y2TYU5z8BItxIV4iSMMJCz5iBikbT1RQmXyay+3aiMAJejSmCq2xItFwWugxOyTJJIsfZAtIloBk8UMsHsnmh1gCksUPcoRBltAtUkdjlovADAdVlvMaU6XGWJFgqXWJtTZ5ll2aZpNmWuUZPc6zT9FgDoLRj00BMPqwKYBIGqt4kCZAORKULLhUqWISQjUPTVWjEjnq8PGqMFf9ed8/VxYuuzX/xp41eGAZEUKAQSOEHbR48H5DDKnaQqS1QvcCk1DYGA8fs/IyRKhMhbT8Ls5MsNbm1u/IbNprlWfUZs5rSp0asSCOCCHf1pxtPZvcflxjrrLJ0luTxrdpJxgkykr9qzQdLEp6gs38MeDMTeLTNjhuhVzxhTohV4LHhDI1mqxCzIH2C+2eoS6EERrc+Ls2bPB1FVmPcZPlM4NrV1HibyNlgCLQGM458AEjHaRgqoYYLbuKuy836BK79SzSxwl5WaSXHfKyST+L9LFDvshLbtAl9hgogukUJzrFyS5Rokuc5BQnOkVJHdZ6dsijMZ/TGM+qTWfljmabLEOvLjGqCo2qougd3C7BEcTGAJj8EPlpD2E5B9RcpOKBmodUXNz9jkWUQqjlo2mai6Z0GPBX57/9pPLLJ6c8NjFxXM8aP+CMCCHA4BPCCDof3m/A5+xRfWuYJmvNH1o9x3LUj9X58w4ZIYGLp6mJjK4mpgyKTGs7klO/XWU535RS1p44zknwSL40zJMGOeJI4c1rARbp0xorUtqPJ+uOAxXWJY1vSxrfrhkdYgsBwO47c073mkYyM0O5EgEDACgM1U58yopbfVAoRWPl3RkF2QRMVBHTNMDtvgzV4GB4CCEA0BifscEZG252dxFlYfEcO69flSJfmKZY1jnpI0mSzR58zM5q8+LxCShbjBJ5wOjFDigvYJc42yVundilk7jbJa52sVvv50hc4kROyCV2G00JeQZVkUFdYkrIjdXCBwBBGkx+MAWwwQemADb6gUMgFQ+reUjNAxUXKbkQ09SreyGUsFGBFIrll2a8swec/zz8qj3g/PP032sEg2v8jIkRIQQYrEIYwRKAg0b6vAN8UWQftnqOVhveUotnpCpur3KyDxiBReAyNcoToy57Nd9nzW7cpbJUs/wOHunlBt3coItisAIcUZAtCnIlAbYowJXYJalWRZZdkhaNA+rghEGRIrde4m6XuHRit07iape4deyQx5hQ0JY0vlkz1sjTdNgIKdpfb/7Y4j5akPiIjD8KAIx+OGnFFXas4sLYBFQg6S7dlJiNJiphfMIgcgq9KsNGCDsIUBFnSFznhAD1Y98JkdYqw2sU7c/XPNwRidthI7QG8TEzNHuxJQBaPqTwUYoQUgQoygxB3YAwFnhNErcuxBZY5Vkxed+QNFgCYAlisx9MATD6wRvGSi5oeBcWfGoe8Hrnm9KlEIrZUCgliuSQzO9iV39X8/7Xjr+7IHPOvaPuZBJDyTXmckaEEGBwC2EEDKD34gY3bnBBswdTV/4iScpVa3zP5a/JVt+vEE6oceIjZqz3QYYI5UkgR9x1l+68NcoO+9kBFy/o4gRc3JCLF3DIHM1yW4PU1ebhK23yDJss0yrLtMkyvHxFn33iq8MKB6TOVjbpZZIBBh0GTLFDPgBgUQGCDiOaZpF+VtgfmZXz/TYPX+mUJLlESU5xokuU5BRrPXxlxHWos43Q7D5Ya/xAIRybpfwpCYIKOz5lxb4wGiOH0QqQdbu/pOGjUiUaLY8qmn5QMfyEsAMa4zYv1DjhnIO2BgAAMGC947tGy/+00nnp8mUEwb3cWSZEQbsPWry41YvbPCBio2QBpAohhQ8Kbt9aen1hMAXAGsSWAJj92BJEvjBWcCCBixK4oOSChgsyTpzNzZ2FUMxGBVIokqEUYdc30XkMrxx9x+gz/37SI4UJufFsxwAxIoQAQ0EIOxOgcKMb6l1Q56Qdoa7PsXtP1xj/y2Ups5R3C7npfgpqXbjaAQ1urOJBrgTliaFzAcVobIQEpiTONoW9QW5vVNga5PZGhGmrLDPIEQJAmMmnGT2ZFYYZHI9A6RUkePlKj0Dp48q69GslMC1262SOJrm9WeZoktub+H6rQ5wcZIsoJifMYAEiQmxB5IIUwcIIkWw+RbCdkiSnKNEjVHczB48IIYna6o0fkrRXJHrQSBY0e7DeDzliNFaBMoRdNuoCCKFsMUxRocwhW3F7GAthZ/Q+XOOEOhfWecFP2upMHzp8lZkJtysE02gaOgthZzCAyQ8tHtzqg1YveEJYxEIiNhazkJAJEg4ImSBmIREbi1goyl0ADOAhsTOEXCR2hsAZAkcInCHsDAEAKHmQwEUJHFByUQIXS1h97mVF+10ZKomWhwplcCX9AwBXyP3p2bVbGr5fWXjLbflLhvpCsIMRIQQYakLYGWsA6lx0gwtavZfunWJMtTu2NVvXSvlF6YrbBJxUAKBoaPJCtRPXOIFBQJ4Y8iQoRQChYE+cZfh+m9zexCZ9AMAM+wiqJ7VPWVRA4DULvWaB1yz0WThBt5ev9AgUHoHKy1eRLI7U1S5zNEmdrT6e3CZLt0vT7dJ0mzTNJU6MV3SX3l3bYP4ySDaaYXljcLaaR6QKUKoQpQqhe3dQJgGjFWiyEil5Q1UCI1wjQtgBSUOrBzd78DFDzYHmD4NhZ7JkaZJiVjTV4sIY3CHsJpGLxG4SXCS4SXCT2BVC7jBmIWARF+SQxcARAzGHuGCx4zKQL4ydJHKRmMcACQskbCRmgZQNEjZI2EjCjjZRe48RMEHOQQouyDmg4CI5B+QcsJkMWq22m79yhdxfnf92fc2mWall95TcLufJ+raV/cuIEAIMZSHsAGNsDqBmD93igWYPuH4o3kTTgTb7llb7NyJuZor8Jim/5MdCFn6oceJqJzYFQMbCSi5S81ECF1RckHPQQDk5MihS4LMIfBah1yTwmtmkzylOsskyHJLUjrz4vcdDgjUIlgBtdJ2kQt+wcZufuDFBckOaiHNVLwmEkJwDWh4kCWCUgrg8dnMocq0JYWfCNP62pnxtzVeOkCNZdpNSPItB9LyAeoACksKRjDchCkXKRwVpuPALhbkMJGFjCatrs31cIBASMLGYjYQsELFAxEIiFghZIGaBlNO195Zer7+SEOo8hnXV325r3DU9efJdxbclCodYRtZo6E8hvBb7WL+BEFLxQMUjSpUAAPYgNHtwsxu3eLkEsTRZvsjo3F1nfJ/GVKJ0rlo8g82UaXig4aEZGhTG0O4MOsIMK8mosIMpAK4QLWODkotUXEjggpiNJGwsZPaHOlIMlkukdYm6m5zGRCAM1iBYg9gaBGsQbEFsDQCfsCWzd0thh4AhTExYnCaZQoVpHq/rnTECIQUXtHzQ8kDLR1o+GikuOpxgEmh+aslsbX4bqf/i/MZjjf8bqy0rVl0n5uc5gtgRws4QCnZjlr8YLqN7P2HU6WdPQAB8JvCZwGchPhP4DBBEfmECjwFCFhKysJAVh51UkiIPtJdvqt9ea6tfkDn3g4WvJgyoT8CwYUQI+w8ZB2QcNEaBAMBDYqOfa/Jfbw5cX2E+V2nccaTxUTE3RyWemiCcyGKImQjUXJzEBOYP9g2KRpYgmAPY5IcqB7hI2kkiH0nzmCBmgYiFxGwQMUHMBjEL8ZjAYWAOgdgE9LOTCEmDn8K+MPKFwRvG/jD4wuCjwBcGL4l9FHhIoGiQc0HBQQoO5AjtLP5Rv3+/P9isFE9NlP5OxM0GAIqiKLhgbmUSIOMgKRtkHKTgYC0fafiIPdScX0boAUUJ+c9Oy7f67Vvqd2yqeQUDXJc2bX7KlGxZpj8MjhB2hsARwu4Q+CnwhyM/sZ+CAIVCUStlNxAIeAzgM0HAQgImCFiIz8B8JghZSMACPvPCy6sZDHslgiRFnjCe2d1yYH/bkRxZ5oKsOc/PeJIde4DHCFdiZGt0UIAxNvlD3zcfPdh+sNpyQsxLlfLH8Rj5EkEuq9voCIzBE8YuEnlI7CTBHQJXCNxh7A9DkEJBGodooDBwGBARRRaBuQxgMy4kpCIQdCykEMDliyoWAZc7GmCAAAWBMIRoCNI4SEGQgiCFAhQO0kAg4P0wOgiYiM8AHgsuTJOZiM8EIQv4DMoVqLF5T1s9xwMho1w4TimakiAYx2KyhEyQcpCUDUIizA770lUSGQf1rPraUOda3hoFAJ/PR5KkRHJRkd5qa92ulv17Wg5SmJqcOKFUO2aMukTE7jrYLkzjC+oYBgqApCEc+YeBpPEPv6MwjRkE4hA4snDkMiMryAu/D+B861TDmVZaX647cdxwOlOaPiN1ynWp066dJeCIjRDgGhPCzpAUedpUWa4/caTthMlvyZLnpUnyNMI8qSCTBomXBDeJPSQEKYiyYhQGCIYhIoohGgUpCP4Q1UxhCP3gTEPjLi5I0nB5qCQC4DIi4gocBnAYiMMADgNzGajDAeFyaNruCzZ4AjV2/zmbt17GS8yWjSpWjytSFEq5DAETCVm48/5VMBj0er1yuTyqDzkcGRHCy4Wwg2Zn6xHd8aOGU5Xm81qhukRZWJiQl6/IThYlDZUyQ5cTpEL19sZz1tpKy/kzpqowFS5NHDsxcWypdqyU0/X3MIwZsRFe07AYrAnaMRO0Y1ZmLA0zqVpXQ6Wl+qR+Q7WtjsfkZUrTMqVpheLkZHFSkiiRxxQHKBSkIEDhIAVBOrI4gwAFABCmMUkDAIRooDECgAAFGENkoyZEQ49qLwIAMBEwEXCZwETAIoDDACZBsAnMZiAmAg4DAiG7LaCz+PQGT1urq7nJ2UhS4XxFTqkyp0i5tERZKGBdnm5uqI5fI/Q/aZKUNEnK8oKbKUxVW+vPWs4daj/6wZk19oAjU5qeKU1Lk6SkiZOTxYlqgXJwVl0IhANtbn2bW9fiamt0tDQ6mnUeQ7okNV+RM1E77v5RdxIe6N5rdIR4MSKEgxoRSzg1aeLUpImRl3qPscHR3ORsOWk6u7Fua5tbT1KkRqjWClQJfEUCT57AV8i4EhlHLBGIJRzxlbaM4oIr6HYGXa6Q2xFw2QN2nd9m9lktPpvRa9J7jVwGN0mkTZUkp4mTJycuzpSmqQXKvmvMCNcsDMQoTMjtCCH3kr56e2Ojs6XJ2XJYd6zdrbf4bAl8uZqvVAmUSr5CwZPJuTIZTyrhiCVskYgjYvVZjcNAOOAKeZxBl93vcASdtoDD7LNYfDaTz6z3mHykL1GkTRElpoqTpiaX3lm0LE2S2rkxeo++jxo2wiWMCOFgweFw7Nu3z2w25+fnT5kypUvbu1ao1grVZckTO454SZ/BazJ6zWafxeq3VZrP2wMOR9DlCDitdquuopXhJ5TpqpSiNC6Ty2NyeUwuk2DyWDwmYgAAn8XrcrJMYcpH+gEgTIf94UDkp8PuaD7d6Ha4WFouM40rZAukHImYI5JxJTKuVMGT58mzpyUnaARKjVDNiyKswmKx7N+/3263FxcXl5aW9vy7+wGz2bxv3z6n0zlq1Kjx48f3/oIjDEJaWloOHz4cDAYnTJhQUFBwybsCFn+UqmiUqqjjSJimTD6z0Ws2ec1mn7XNrT9jqrIFHK6gyxl0u4LusC0UbPSygJGUn5qUmcwkmDwml0WwuEwOALAZbM4V3FK8pI/GNAC4Q57IT5IOB8NBd8jjI/0e0sMkmCK2SMIRybhSGVcq40o0AnWJslDJV2gEakXUYX/19fXl5eUY44kTJ2ZnZ/fgSxuhe+IshK2trf/6179MJlNpaemvf/3ra9a8EStffPHFIw8/kinIlQsSzhn+nlSgWbNmTTSJ5wUsfpY0PUuafsnxTz/99DdP/CZbVCDjyys/P+MuMP3jrX+KEyQ+0k9hyk/6w5iCTj35EhiIkSLmAQATMXgsHpNg7ly//cM/r8mXliTw5JW609ljM9/79B2NpufRS++9994ff//HfFmJiCv+U/uzBaW5n376qVLZ81Xj22+//eQfnypUjBJxxP+ne2bUlOKPP/5YobhWPAuuBTDGzzzzzOsvrxqjncAiWI+3/W7xrQvfeuutK2WficAkGIlCTZeRdhjjJ5988u3X3xmbWMogGDvf3TztpmmPP/87CtEhigxSQQAIUaEg1XWyKI1QRSACAIQsAUJIxBYyCSaXyRGxhTwmT8jis3qdK5+iqN/+9rcfv/fp2KRSAHi0/dcr77n95ZdfZvQoY9QIVyKezjIY48WLF7/66qsZGRmfffZZQ0PD008/3fmEEWeZLqmoqJg6qeylpe/mq4sAgMb02/tfaWBWbd68mcfjcbkxR6yfOHFi1rTZryx7L0eZH7ngqr0vGoTNe/bs6VkLy8vL586a9/qtqzMTciIXfHX3Cw6FcceOHT274IEDBxbOW7Tqtg/TFVkAQNHUS7ueDyV5Nm/e3HFOTM4yu3btumXxrauWf5gqSweAMB1+ccdzzGx6w4YNPWvhYGDEWeYSZ5n//ve/f/3D31+97X05XwEAftL31DePzbh16osvvtizW7zxxhv//tNLr9z2vpQnAwBfyPv/Nj56/V1z/va3v8XlI/QSvV7/8ccfr/7Pxy/e8paYKwEAT9D9u/U/v+PhFU8++eRAt67P6U9nmXi6Bp/9/+3de1QT2RkA8MskEIPyEMMqbyQLsrzOrqAguwVRpIIPErpdSbb0sLXaPVbB7fqoZ30dXGwVaUX0dFvXBETPciwi/mEJ0XoWughq4RAMKBiMEGE3BoWIRZtn/0jNRkTBSZxJmO/3V3LJ+fIxM3e+zJ2ZO1JpXFwcm83GMIzP57e0tNgw+BRWXl6+OvLnpiqIEMKcsN98kN99Q9bV1YU7YFYMz1QFTQE3/ORzybV2mUyGL2BZWdlH7+aYqqAp4Mbkrc3fXe3t7cUXUCgU8mI/MVVBhBANo+Uv/v23/6wfGBjAHfAXcWtNVRAhRMfony35QvSPOpVKhS8gsEMCgeDTDzZ7Pbt/gOns+lnKF0KhEPeveYFAsCHpc89nQ5SuLtPzF+8QCAS2SdcWBALBxuStpiqIEJrBcNuUvN2uMpwabPljUy6Xh4aGmt86O48zLHDw4MGSkpIX20Ui0dy5cy1b9Hr90NAQkXd3kKW7u/td7/ctWzAnLNiL3dHR4evr++phn5cFTGQts2yhYbSgWSESicTNbbIP6R4TcKn3assWOkYPmhnS1taG44AVIXT79u2VrGzLFmeai79nUFtbm3nMR6PRmI4JJhlwjV+SZQuDzvB195dIJNHR0TgytAcPHjzQaDSUPSJ88uSJVqt9+vSpuUUul7Pfee65Cv4zAx8/+k9PTw++DVsul4fMfy5g8KyQQdVgX18fjn5ncyqV6u7du+ylz2XIZoX29vYqlUqysiLMyMgIjUZzdR3vaebj8fT0xL3WbNnHdDrdhJ1206ZN27dvf7Hdw8MDw547PNXr9RiGWXPSyFEEBQUNtCjGNA6oFWw2m8Vi4ag0gYGB/Z3PBTQi44D6XlhYGL7lGRgY2N/zYkDFvHnz8AUMCAjo7++zbDEYDd+r74WHh5sDvtbQaEBAwL2HfXGBCeYWvUH/w6MBy4COCIZGLYdGfX19+9UKHw8/c8vg4/su05yDg4PH7D0mydfXd0Ct8J7x4/mXHx4NuLm7+fv7W5O5reh0Oh8fn361IuytH68J6lcr5syZ49Bb9SQxGIzXGhrFtw2Y2LKPeXt7t7a2mt8aDONchcFkMid5/YLRaMQwzJr/zVHweLyVglUZkZy33P5/Pv+cpNLL3zM6OhrfEuDz+VkrPlz+zirWsx7+99ZTfmyfqKgofE+O4fP5fO7HaeErvKazTC2V/y5jR4SEh4fjiGYK+KvsX6eGZ5hHpU5d+zp6flRISIj5M9gzkwy4IXdjSliaxzRPU0vZ1a8WJi6wkz0aPq+1BKaeF/99Ho934s/HYnzfc6EzEEIGo+Gr7w6vWbMG928FHo/39V+O/inrr840F1PAvzWW8Hg8O1nmGIbxeLzj3xz5Q2YpHaMjhPQG/fHGUvvJ8I0icvu3ZSGMjY09dOhQXl6ek5NTX1/fy6aEAGMkJibu+nLnuj3Zi0PSvKbP6vi+XUUbqK6uxr0FJCcnb9+9de2XH6Wwf+rp6iUdaBtyuV9TU4P7+Wmpqambd+R/8scPl7y93IM5s72/ZYQ5dP78eXzREEIZGRmf/m59bnHW0rB0t2kebfeuP3V7bE1ADodz/bfXc0uzloQud2O4tyqu6b3+e/4b/AGBHdq8efONGzdyTnJSwtKcMecr8gb/KJ/i4mLcAbds2dLR0fHLk9zFoctoGL3xzrch7wUdOHDAhjlbadeuXbdufZxb8bOkt5cihP7Vcznm/ci9e/eSnddUY+Mp1qqqqsrKyphMplarLS0tDQgIsPwrXDX6Cnfu3KmtrVUqlREREVwul8FgDA0N4btq1EQmk9XW1qpUqsjISC6X6+Ji7RS93d3dIpFocHAwOjqaw+GMew74tdy6dauuru7hw4cxMTGZmZljftfjmGKts7NTLBYPDw+bAjr6JeZw1ei4U6xdvXq1vr5eo9EsXLgwLS3N+i9qampqaGjQ6XTx8fGpqanWB7QV82OYGhoaGhsbEUKJiYnJyclk50UQmGsUIYRUKlVhYeHhw4ffdFZ2q7i4OCkpySZ3mjui9vb2Cxcu7Nixg+xESLNt27a8vDyHHt21hlgsVigUa9euJTsR0uTm5gqFQtwDOY7u5MmTLBYrIyODgO+y34HmkZGRqqoqsrMg06VLl3DfnzAF9Pf3i0QisrMgU01NzfDwMNlZkKazs7OpqYnsLMhUUVEx7pUWFNHc3CyVSon5LvsthAAAAAABoBACAACgNCiEAAAAKI3Qi2WOHDly6tSpSd58Njo6Wl9fn56e/qazsltXrlwJCgry8/Ob+KNTkVKp7OrqSkpKmvijU5RYLE5ISHB3dyc7EXL09PSo1er58+eTnQhpzp07x+FwKHuxjEQiYTKZYWFhE38UIYRQfn4+7mfOEFoIEUIVFRVUmDUNAAAAkVJSUsbcsDd5RBdCAAAAwK7AOUIAAACUBoUQAAAApUEhBAAAQGlQCAEAAFCaA8zne/z48YsXL5pep6amrl+/ntx8iCeXy0tKSpRKZWxsbF5envXTZzuo8vLy9PR0iszDrlAoioqK7t+/v2DBgvz8fMpOvU2plW6J4r2e6N2+0UEYDAY+n9/R0UF2IiRYsWKFTCbT6XSVlZUFBQVkp0OC0dHRnTt3+vn5UWQDMBgMGRkZMplMr9efPn163759ZGdEAqqt9DGg1xsJ3O07zM/MoqKitLS0iIgIshMhQWZmJpvNRggtW7asrq6O7HTIsXLlSkd/rNLkSaXSuLg400rn8/lcLpfsjMhBqZU+BvR6ROBu32HOEd68ebOqqqqyspLsREiwbt0604s9e/aYX1MKk8mMj48nOwviyOXy0NBQ81vrH/3oiKi20seAXo8I3O3byxGhVCotKCgY07h79+6oqCjTa6FQaDQaV69enZ2dTXh2RJhwCZw4cWL27NmLFi0iPDWCTLgEqEOn01H2pCCwNOV7/asRttu3l84WFRV15syZV3/GyckJ9+Pa7d+rl8DZs2dbWlqOHTtGZEoEm8w2QBHe3t6tra3mt1R+KB2VUaHXT4iY3b69D402Nzfv37/f9FqtVmu1WnLzIUV1dbVIJDp69Chlp9+lmtjY2MuXLxuNRoRQX1+fh4cH2RkBolG51xO/27eXI8KXSUhIEIvFWVlZdDpdo9EUFhaSnRHRJBJJTk5Oenq6aXBgxowZAoGA7KTIwWAwKHIRuaura05OzqpVq5hMplarLS0tJTsj0lBnpVuieK8nfrcPk24DAACgNHsfGgUAAADeKCiEAAAAKA0KIQAAAEqDQggAAIDSoBACAACgNCiEAAAAKA0KIQAAAEqDQggAAIDSoBACAACgtP8BrWrXNldFJiQAAAAASUVORK5CYII=" }, "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 }