{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Dependent density regression\n", "Author: [Austin Rochford](https://github.com/AustinRochford/)\n", "\n", "In another [example](http://pymc-devs.github.io/pymc3/notebooks/dp_mix.html), we showed how to use Dirichlet processes to perform Bayesian nonparametric density estimation. This example expands on the previous one, illustrating dependent density regression.\n", "\n", "Just as Dirichlet process mixtures can be thought of as infinite mixture models that select the number of active components as part of inference, dependent density regression can be thought of as infinite [mixtures of experts](https://en.wikipedia.org/wiki/Committee_machine) that select the active experts as part of inference. Their flexibility and modularity make them powerful tools for performing nonparametric Bayesian Data analysis." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "from IPython.display import HTML" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from matplotlib import animation as ani, pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import pymc3 as pm\n", "import seaborn as sns\n", "from theano import shared, tensor as tt" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "plt.rc('animation', writer='avconv')\n", "blue, *_ = sns.color_palette()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "SEED = 972915 # from random.org; for reproducibility\n", "np.random.seed(SEED)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use the LIDAR data set from Larry Wasserman's excellent book, [_All of Nonparametric Statistics_](http://www.stat.cmu.edu/~larry/all-of-nonpar/). We standardize the data set to improve the rate of convergence of our samples." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.5/site-packages/pandas/io/parsers.py:1961: FutureWarning: split() requires a non-empty pattern match.\n", " yield pat.split(line.strip())\n", "/opt/conda/lib/python3.5/site-packages/pandas/io/parsers.py:1963: FutureWarning: split() requires a non-empty pattern match.\n", " yield pat.split(line.strip())\n" ] } ], "source": [ "DATA_URI = 'http://www.stat.cmu.edu/~larry/all-of-nonpar/=data/lidar.dat'\n", "\n", "def standardize(x):\n", " return (x - x.mean()) / x.std()\n", "\n", "df = (pd.read_csv(DATA_URI, sep=' *', engine='python')\n", " .assign(std_range=lambda df: standardize(df.range),\n", " std_logratio=lambda df: standardize(df.logratio)))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rangelogratiostd_logratiostd_range
0390-0.0503560.852467-1.717725
1391-0.0600970.817981-1.707299
2393-0.0419010.882398-1.686447
3394-0.0509850.850240-1.676020
4396-0.0599130.818631-1.655168
\n", "
" ], "text/plain": [ " range logratio std_logratio std_range\n", "0 390 -0.050356 0.852467 -1.717725\n", "1 391 -0.060097 0.817981 -1.707299\n", "2 393 -0.041901 0.882398 -1.686447\n", "3 394 -0.050985 0.850240 -1.676020\n", "4 396 -0.059913 0.818631 -1.655168" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We plot the LIDAR data below." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n", " (prop.get_family(), self.defaultFamily[fontext]))\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFpCAYAAACrqZC7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0W+WZP/CvJFuLYzteopDFoRA7TmiDQ4ibDdJAGtaZ\n/OISSIAS6DCF/soAbWmnJCkTyvIjZetpC+10OWylEwgHAqWFQgmGlCUhe0ymTYyd0ux4kxdFtiRL\n9/eHexVZvpukq3uvpO/nnDnTWLb0Xtnoue/zPu/z2gRBEEBERESGs5s9ACIionzFIExERGQSBmEi\nIiKTMAgTERGZhEGYiIjIJAzCREREJikw+gXb2/uMfkndlJcXwecLmD0MQ+TLtebLdQK81lyVL9ea\n7dfp9ZZIfp0z4SQUFDjMHoJh8uVa8+U6AV5rrsqXa83V62QQJiIiMgmDMBERkUkYhImIiEzCIExE\nRGQSBmEiIiKTMAgTERGZhEGYiIjIJAzCREREJmEQJiIiMgmDMJGMYDiCNl8AwXDE7KEQUY4yvHc0\nkdVFolFsaGzB7uZ2dPUGUVHqwsxaL1YsqoHDzvtWItIPgzBRgg2NLdi040js3529wdi/r11ca9aw\niCgH8baesoYR6eFgOILdze2Sj+1u7mBqmoh0xZkwWZ6R6eEefxBdvUHJx3x9A+jxBzG2vEjX1ySi\n/MWZMFmemB7u7A1CwKn08IbGFt1fa3SxCxWlLsnHykvcGF0s/RgRUSoYhMnS9EoPa01luwodmFnr\nlXxsZu0YuApz80xTIjIH09Fkaemmh1NJZa9YVANgKMj7+gZQXuLGzNoxsa8TEemFQZgsTUwPd0oE\nYi3p4VQqnR12O65dXItlC6vR4w9idLGLM2Aiygimo8nS0kkPp5vKdhU6MLa8iAGYiDKGM+E8EQxH\nsnZWl2p6WEsquyrJsVjpfbTSWIgoNQzCOU6P7T1mf9inmh5ON5UdT89tUum+n+zoRZQ7GIQzwOyg\nFS+d7k+pfthn4vpTeU4xlR1//aJkK5316KIl9342LDgT/kBY87WxoxdR7mAQ1pERM5RkgpHamuiy\nhdWKz5Hsh30mrj/d59Sj0jnd91Ek936+33QcwVBE07XpNRYisgYGYR1lcoaSSjBKZ3vPQGgw6Q/7\nTFx/us+pR6Wz0vvY1auti5ZS8BwIDRWIabk2dvQiyi2GB+Hy8iIUFGTvnbrXWyL59YHQIJpaOyUf\na2rtxDeWeeB2pv52/+aVjyWDUZHHiZsazpb8mZLRHnjLPWjz9Y94bEyZB9VnVMqO6XjHSXT1yX/Y\nO5yF8I4ZFftaJq5f7+eUKsIaCA1i0GZHealL9rmU3kebHfjLxydwc8PZcDjkZ+ZK72cipWtL53cK\nyP/95iJea+7Jxes0PAj7fAGjX1I3Xm8J2tv7Rnw9GI7g4NEetEt8MAJAR3c/Wj/tTHmGEgiG8eeP\n/iH52F92H8GXZ05ASZFTclxTJo6W/MCuq65EX08/Rl7NkPLRHlSUyBc1RULhYe9Fmy+g+/Vn4jlF\nYmahqbUT7b5+1cxCXXWl5NpyNAq8/uGn6O8PYeUl02JfT1w2iIQjsu9nstcmNxa136nc328u4rXm\nnmy/TrkbiLxIR8uto+pZpdrZG4TdBgjCyO+TqsRN5rXXv/VJLGWZqNsfwg+f3I5Z004FkMRxuZ12\nADaEwhHNa6JuZ0FSRU16ViJn8jlFyaa5VyyqQSQSxeY9xxCV+B1v3nMMsNmwYlE1Xnz3oOSygdz7\nmey1saMXUe7I6SAst4565QWTZT8oxVmQliCZ+EEu9eEMAHXVFbHnKnDYVNd2418bAPb/o0vxOn3+\n4QEkcVwDoSgAYP70cVh5yVTNNxzxH/ZdvQMYXezEzCnSH/Z6ViJn8jkBoC8Qws79ya13O+x2XDL7\ndLy7+5jkz0UF4J1dR9FypAeH2/yxr8cH98Tg6Sx0SN5cqV0bO3oR5Y6cDsJys50Dh7oVPyjlgmQ8\npUKbeA47sOV/T+Dd3cdQUepCkbswqdeeeno5fH0hTde7u7kDS+afITuuA4e6NT0PMLRW2tkzgIYF\nkxGJCtjT3IFufxBNrZ1wOFok07aZmKHp+ZziTdmO/W3o9ku/p0rFTUozc9HRdr/k18XgHh88i4uc\neOW9gylfm9jRi4iyV84GYaUgqfRBGYkKeGfX0djXxCAZiURxx3VfjH1dqUo1XiQKRP45E+3sDcp+\ngO860C752h/uOwG3U3rGlMjXN4AjbX5dDjwQ10pdTntsJi2OSS5tm+wMTUu2Qc9ZX+JNmRSlVLDS\nzFwklw2Jf+/jg6eWa7PSvnMi0lfOBWHxAys0GJUNRnIflJ29A7KBe/OeY/B49uIr558Bh92uaVaU\njK6+ID5sOp7Wc5SXuFE1tljXAw/iA3A8pT2pajO0VLZbST2nXnum46mlgtXWhu026b+v8hI3PK4C\ntPkCI8Yr936xMxZR7suZIJz4gVVe4oRL4wwynlyaMioMVcGGQoOxGeC008vxwb4TaY9dFByUDnjB\nUATnTR+H/Ye6VdcSS4qcKa+jag1UwND+2HZfAFVjk98yILdMEBgY1LRmnUxw0nJTJnI7HRAEAZFo\nVDbIOez2oSpom21Y1kI00Vs8bLlBVOQuwL1Pb08qmLIzFlHuy5kgnPiB1aVxHTWRTabCWSSmjZta\nOkZUHpcVuxAIDiYd+NVUlLpx3SVTAUDTWuKKRTWn1nFPBlHxz8cbFpwpORNLJlCJBAA/fbEppT7U\ncoH+w30ncOCQT/U5tQSnVG7KBkIRvL3zKGw2m2qQu3bxFDjsthG/g1NFf6e+XuQukK0DkHsddsYi\nyg9ZH4SD4QjafQHZDyy304FR7gJ09QZhk0kVxlMKwMBQ2jh+BpRYefzS5lbVdcfxFUU43qV9v3T8\nDFZtLTG2ptvSAZ8/iLJiJ6ZPLocgCLj7iW2KVeLJZg9SmZmpraWrPafW4JTOTZmWICe1Vg0A3X1B\nLFtYHfu6xzU0A072ddgZiyg/ZG0QTtwLK2cgFMGdX52JQP8gHnl+T8bGI1YeJ27rcTmHPmDj9+g2\nLDgTdz+xTXLc4k2Dry+IsmIXpn2uHA0LzpR8Tam1xMTg0+0PYfOe4WvNclXiSoHKVWhHMJz8+nAi\nrWvpYqV3f3Bw2E2GluA0util6aZM7n4rmSDnKnSgcrRbdivcb99olr3WVCux5db2WcBFlH2yNghr\nqXQV/WXvcSy/sEbXQqpE8R+oUjOkxA9HuXXb8+vGo2HBZDz3VjP2H/Jhy74T2P+PLkz7XAWuvWgK\nilyFsmNIZk0XkK8Sj+d2OnDe2eNwft143PvUDsnAlWzQ0tK0orN3AD98cju6/cPXULUEJ6VAHQpH\nsOa6cwEMpdP1aASidStcMq+j9D5NPb1s2L9TKeBiwCayhqwMwskGm10H2hCNCjg5EE76tcpGOdF9\nUj2VmfiBmjhLTQxQSvtfNzS2DCv46uoL4cN9J7CruR3n142X/XDVum1KpJaaB4YyCTabDd4yj2yq\nOtmgJV77rgPtiv2Uff6hxxJT1HLBaUpVKQD1WaT3n9uE9GgEEgiG8b5MVbvaTU7864hB0eMqiM3+\n5bIqWxLWzpMp4GLFNZG1ZGUQTjbY9JwMD7UVTOB2OjD7LC+aWrskq6LdTgduuWI61j27SzZ1KUq2\ng5Pc/le103aU1kuT3TYlt50mkbh/Wm6tOJ1r/92bBzRXmItpb7kgvvWvbdjT0onzzh6Hc6aMwds7\nR1Yvx491xaIaFHmc+GDvsZQbgSi1FFV6b8+bPu6fBXRDQXHXgTZ09YViv5PKuOAo9T7F71+XO+Qi\nMaUPKBe1sQMXkfGyMgjrtUe3yFWAaxZPhbNQuphqIBTBln0nFF+rUqajllaJM+b27n7VGwypNdhI\nNIqXNrcmNdu32QDVuwsMbUfa09wh+Zjb6UDDgsmaXzOeq9CBr10+DR53QSwjUDrKqamb1bWLa0c0\nNwFOVTgvmjURi+urRmQaEivEb2o4G5fNnpRS8AmGI4otReVucipKXLjukqlw2O1Yv6lZsvVpYnDc\nf8gn+Rq7P+lAj8z7lZjSnzN9vOwN3vtNxzk7JjJBVgZhreuKarr9QfT4g2hYcCbebzouOaNpau1C\nXc0YyT2hib2YtayzyX1P/IxILS5KrcHKrZG7nQ7MP3scbAD2fNI5bJ9xRLrOaoTRxU50++XXWLt6\n+uEvdGhumhF//YkZAbGaWG2tNhiOYO8n8ksSez/pwP03zR3RIjKxQvzW5TNTbv/Y4w8qthQ9rVy6\nCv7cqV7VrIdod3MHvjRjguyNWY8/hLJiVyx1nyg+pf/6h5/Kvs5AKJLUucZEpI+sDMLAUCoxMDCI\nD9NolhFfyBOUSSn6+gaweNbQSbRNrZ3o6h2I7buVOrVIbiah9j3JFJolrsEqfZgXuQpw1QU1cBU6\ncOUFEZzoOokf/W6X5PfKVUDPnDIGTa2dkoHRWejAT19sUj2M4tTeZunrjw+EWtZqe/xBxWruzt4g\nDh7tweSJozG2vGjEjDP+POaG886QfR4lShkZt9OBVStn4Q8f/F12P7eWZRVf3wAgCLKvU1HqRl1N\npeRNohStSxAA9yMTGSFrg/BQ56KpOHDIJ/nhpNZ0Azj1oa70YVpW7MKmnUdizTnKip2oq67AikU1\nGIwI6OwJ4M3thyX7TQOnZhJqa3HJFJolrsEqfZiLs32xZ/Fb24/IbjUKhqOY+/mx+ORI74ig4XBI\n3yTIzaASD6NIpge1lkMbRhe7UFHilA3EdhvwyPN7UFHqQl11pey66dZ9x3HZ7EmxmWlicZRSAFLK\nyJxfNx4lnkLF3tBallXEQjK51ylyF2DFouphjUOUUvpaAzDA/chERsjaIAwofwhWybQPBIDK0uEf\n6krPM8pTOCzAdvtDeGf3MbQc7UVgIBw7R1iKOJMY+t/yDSaU0o0AUFbsRO/JkGzhkNY9pWprmADw\nyZEe1FVXYvnF04DByLAiJnG8vr4Bxe5gUgdhJNOD2mG3Y9nCanypbjwgVmYnBENXoQPnTh0rmz2I\nX1t9R+b4QQDo6O5HV+8A3tl9VLE4Sm5tVMsNg1y6W8uyinjDtWJRjeSWp8Ntfrz47kHNKf2KEhdm\nTBmDppahrE7pqEIMhCKSN2bpntlMROqyOggD8h+CUu0D66orsLh+EipK3SM+1KWep666QnYGFf9h\nqHZyDgDFBhNK6cbKUjfWfq1ecWZW4LChyF0o+fOJKVy1YxHFoFVS7B6Wpk1cuw0NRnH3E9skn0Op\nkEvqe8W0savQoSm1L85YGxZMRlQQ8OHHJ1S7fMmlYceUebBp55FhNwxSxVFya6OpnvIUfw3AqUpv\nqRsAABiMCAjIFN3F38iopfTPnXrqOcXjKV1O6RuM+HOwmZImyoysD8JKH4LJfDhKPU+PPyh7iLsW\n8TMJtX2rSuugJUVOlBQ5ZV9nQ2OL5Kx/0tjikSlcjVXl8WnaeOIHfTAckX0upUKuRLa4tPHMWi8E\nQRi2tUjLWc+P/Md8dPUMYONfDmJPi/RNk9yNUv1Zp+GjfcqnV2lZG9Va3CV3k3Hv12fDHwjLpsKT\nbWMpdVN53owJWDLvdGxobJHMUridjlgP9FGeQjS1dsbOwWa1NFFmZH0QFiml/JJZ04r//nS3Qs2s\nHYMChw0bGltktw7FpxuB5A+vVyrKCgwMYjAiwGE/dW1aq8o7uvsV1wMVm10oFHIlSpx1up3SgU7p\nrGdgaF390Gd9sq9TUeLEjCleNLV0Dnt/lyyYjD8pVA0Dp06M8pYX6X6msdRsW+qGS61uITQYRTB8\navlA6qayakIZDv6jEzv3yxfxrVk5C+/sPqpa40BE+siZIJwJyW6FstuGtt1WJHS/kts6JHa/AlJP\nayrtK1abIXX1DsgeajGmzKO6Hqi0FNBytDelmxe5tLJSiltcV1dKtZ/1uQqsvHgqghcO3yJVMtqj\neqMlAHjgdzsB2BAMRVKeGaZzMpLS32IgOIi7n9gmOS7xpjISjeI3r3yMv+w+Ilu01e0PAsLQCWGp\njJGIkscgrEIq0CQeTSdaeM4EXDL7dE3dr4pcBVi2sHrEh3iyaU2lfcVShTWJwT6xsls0d/p41Q9b\nuRuH9ZuaFXsmp0Ipxa22ru52OnDNRUMzuMT31+0s0HSjpbWyW0m6JyMl/i2K+7217O/VsgWuvMQN\n2Gw8vYnIQAzCKuIDjcNZiEgoHEsxS6WO44Oq1q1DqTTT1/KhqtROUgxGcufi3rjkC+jqOqlpLPGB\nLdm+3oncMv2plVLcauvq59eNR5FL/k89sQ2mxkZiSc8MUzkZKV7832K7L4Cfvtgk+V7t2N+GJfPP\niKW1tf5OZtaOgbdMPjPAamki/TEIa+QqdMA7ZhTa24fWHbWkjtU+dIuLnFi/qTnpdoFqH6rJtNKU\nm806HKkV4CTb1zvR/LPHwW4beVOgtFc53XX1xPfgZH8Y9/92Z0qdy5Qkc2iE0o2Zq9ABZ6FD4QYv\nhB8+uR2zpg39Daj9TsqLXThnSiUunDnxn2NJ/2ALItKGQTgNaqljtQ/dV947qPn0m3hKH6o2AN+6\nsg5VY0s0XsWpsaaSZkwMFh5XgWIbxXiTxhYjMDAomU2QusFRC7KprqsnvgfBYvnK73ipzAzVrkHr\nKUdqRYM+//CiNdlK9lGFqKupHFYJPWPKGHx51sRYm9NUDrYgIm0YhDNM7kO3YcGZsvts1dKcSh/A\nFaVDqdlMSwwW5SVOjPI4ERgIywZgcQtM/If6YESQDJhSNwVag2yqNxTxP69lnTiVmaHaNWg9llDr\nGMW/JbnvLR3lGnbCWGdvEI07j2JxfRXuv2kO9wkTZRiDcIbJfei2+QIpF8DodRZuOhKDRVdfSLaF\npNihrGHBmfAHwsM+1B32kWctq0k3yGohd5Zv4k1EqqSuIdnqafH1d+xvUz15SurYRqVmNOLrsQiL\nKLMYhA2S+KGbbpFOquufekim+Kqs2Im1X6uPFQkVuQozOTTdSN08AUipK5bW70+2eloc45L5Z+CH\nT26XzECIf0sOu33EsY1KzWjUbgRTKSYkopEYhE2S7mw23fXPdCRTfNV7MoT+4KBixy8rS7x5Sqcr\nllrBXao3ZiVFTsyapu1vSWszmsSe4+Lf2KmdATx7mEgPDMIm0mM2a0RqNlEyncTycVuL1nXdROnc\nmKXyt6T2egUO24jq/SJ34bA94OymRZQeBmETmTmbTUcyncTybVtLOl2xAOnTqqZ9rhwNC85UfN1U\n/5aUgrfUzYTcjRe7aRGlhkHYAsyYzaZLKliM8hQOVUf3BfN2W0u6XbHEYNqwYDKee6sZ+w/5sGXf\nCRw45NOU9k32b0kueCfbdIXdtIhSwyBMKVH68M6mWb3e0i24E73y3kF8sO9E7N+ZTvsmBu9km67k\n47IDkR5YSUFpET+8xYCb+O98I6bqpWhNzaultINh5bOT9SDeTGiVb8sORHphECbS2YpFNVhcX4XK\nUjfstqF90ovrqzSn5rWktDNN6WYCGGq8ksq1EdFwTEcT6Szdgju9UtrpalhwJt5vOi55SIR49rC3\nzMMZMFEaOBMmywiGI2jzBQxJtxoh1dS8HiltPfgDYQRlznfu9gfhLLAzABOliTNhMl2qzS1ymZkd\n0URWmZET5TIGYTJdqs0tcpkV9pBboUc5Ua7Lz2kGWYYVKoGtzOxq83SLzIhIGWfCZKp0m1tQZllh\nRk6Uy2yCIAhGvuDgYAQFBfyPmIYMhAbxHw81os3XP+KxseUe/Pz7i+B28l6RiHKT4Z9uPl/A6JfU\njddbgvb2PrOHYQgjr7WuulJy3bGuuhJ9Pf3I5Cj4O81NvNbck+3X6fWWSH6dUwwynRUqgXNRvrcQ\nJcoGDMJkOq476ktuy9ety2eaPTQiSsAgTJaRjadJWZHclq8ijxMN551h3sCIaARuUSLKIUpbvrbu\nO65py1eudS4jsjLOhIlyiNKWr47ufsUtX+xcRmQ8/pdFlEOUjiAsK3HB45K/7xbT2J29QQg4lcbe\n0Nii+fU5iyZKDmfCRDlEqdVkV28Q9z69XXJ2q9a5bNnCasViufhZdGdvEGXFTsycMgbXXlTLWTSR\nAv7XQZRj4ltNJpKb3SZzhrHUbDd+Fg0A3f4Q3tl9DPc+vQORaFSPyyLKSZwJE+UYccvXkvln4IdP\nbofPPzK4Js5utZyYJLdm3LDgTNlZ9OE2P9Zv+gQrL56q70US5QjOhIlyVH9wEN0SARgYObvVcoax\n3Jrx+rc+kQzeoj08iINIFoMwUY5SKtKSOg9Y6cQkpTXj/f/woazYKTuO7pPBYQGfiE5hOpooRyV7\nHrBS57LOnoDsmnG3P4jZZ52GrX/9TPLxComAT0RDGISJclgqfbmlOpeprRlfd8lUHO04icNt/hGP\nSwV8IhrCIEyUw+Jntw5nISKhcEoBUW1WXeQqwNqv1WP9pk+wp7kD3SeDqOBBHESqGISJ8oCr0AHv\nmFFpHQWnNqt22O1YefFULL+whgdxEGnEIExEmmg97YoHcRBpxyBMRElhkCXSD7coERERmYRBmIiI\nyCQMwkRkCJ6wRDQS14SJKKN4TjGRPAZhIsoosee0SOw5DQDXLq41a1hElsDbUCLKGLVzipmapnzH\nIExEGZPMOcVE+YhBmIgyJtmTnIjyDYMwEWWMlnOKifIZC7OIKKNSOcmJKF8wCBNRRmntOU2UjxiE\nicgQ7DlNNBLXhImIiEzCIExERGQSBmEiSlsyfaHF7+0LhNhLmvKe5jXhQCAAACgq4poOEQ1Jpi+0\n+L27DrShqy8Euw2ICkAle0lTHlP9iz906BCWL1+OOXPmYO7cubj66qtx+PBhI8ZGRBYn9oXu7A1C\nwKm+0BsaW2S/t6svBGAoAEPlZ4hynWoQvvvuu7F8+XI0NTVh7969uOqqq7B27VojxkZEFpZMX2il\n75X7GaJ8oBqEu7q6cOWVV8Jms8Fms2HZsmXo6uoyYmxEZGHJ9IVW+l65nyHKB6pB2G634+DBg7F/\n//3vf4fDwY32RPksEo3izW2HYLNJP57YF1qph7TczxDlA9XCrO985zv46le/irPOOgs2mw1/+9vf\n8NBDDxkxNiKyqA2NLXhn9zHZx8W+0MFwJNYla2atd9i5wnI/Q5RPVIPwl770Jfzxj39EU1MTBEHA\nOeecg4qKCiPGRkQWpLS+a7cBC2dOxJUXTMb6Tc3DqqZnTBmDL8+aiN3NHejqC0pWRxPlG01blCor\nK3HhhRdmeixElAWU1ncFAbjki5Pw4rsHh816O3uDaNx5FIvrq/D/bp6LHn8QHlcB+oOD7CVNeU02\nCN9www145plnMHfuXNjiFn4EQYDNZsOWLVtSesHy8iIUFGTvf3Beb4nZQzBMvlxrvlwnoM+1loz2\nwFvuQZuvf+Tzl3tQNaEMTS/slfzZptZOfGPZDFRNKEt7HGr4e809uXidskH44YcfBgC89NJLur6g\nzxfQ9fmM5PWWoL29z+xhGCJfrjVfrhPQ91rrqisl13frqitx5Fg32iUCNAB0dPej9dPOjB/kwN9r\n7sn265S7gZCtjh47diwA4PXXX8fEiROH/d/rr7+emVESUVZYsagGi+urUFnqht0GVJa6sbi+CisW\n1ShWQidWQCfT7pIoF6muCb/++uu46aabVL9GRPlD6Yxghx2yldBiBXQy7S6JcplsEP7ggw/w/vvv\no62tbdiWJL/fb8jAiMj65M4IFiuddzd3wNc3gPISN2bWjol9XWxhKRJbVwLAtYtrDRg5kTXIBuHC\nwkKMGjUKNptt2KENY8eOxc0332zI4IgoOynNlNXaXS5bWM1qacobskF49uzZmD17Ni6++GLU1vLO\nlIiSJzVT1tLuMtOFW0RWobomXFtbi/fffx9/+9vfEAye+g/n1ltvzejAiCg3iYVbnRKBmK0rKd+o\nBuFHHnkEH3/8MVpaWvDlL38Zb7/9NubNm2fE2IgoB7kKHaqFW0T5QrUMcfPmzXjiiSdQWVmJe++9\nFxs3bkRPT48RYyOiHKW0xYkon6jOhJ1OJwoKCmCz2RAOh3HaaafhxIkTRoyNiHKUUuEWUT5RDcKj\nRo1Cf38/Zs6ciVWrVsHr9cLtdhsxNiLKcXJbnIjyhWo6+sc//jEcDgfuvPNOVFdXw2az4ac//akR\nYyMiIsppijPhSCSCn/zkJ7j//vsBALfccoshgyIiShR/NjFT15QrFIOww+HAgQMHjBoLEdEIbHFJ\nuUx1TXju3Lm499570dDQMKxzVk0NqxiJKPPY4pJymWoQfu211wAA7777buxrNpsNb7/9dsYGRUQE\nsMUl5T7VINzY2GjEOIiIRmCLS8p1XFAhIstK5mxiomzEIExEliW2uJQitrgMhiNo8wUQDEcMHh1R\n+lTT0UREZpI7m/jKCyZj/aZmyappomzBIExEGZXu/l65FpfrNzXLVk1/65pZuo2fKJNUg/Dtt98O\nm8027GslJSU455xzcMUVV8DOfXpEJEHv/b3xLS7VqqYHQoNpjZ3IKKr/JXi9Xpw4cQKzZs3CrFmz\n8NlnnwEA/vSnP+GBBx7I+ACJKDuJ+3s7e4MQcGqmuqGxJe3nVqua9sk8RmQ1qjPh/fv349lnn4XT\n6QQArFixAl/72tfwzDPPoKGhIeMDJKLsk+n9vWLVdKdEsC0vcaO81IW+nv6Un5/IKKoz4Y6ODhQW\nFsb+XVBQAJ/PB6fTGQvMRETxtOzvTYda1bTbyXIXyg6qf6mzZ8/GzTffjKVLl8Jms+HVV19FfX09\nTp48ySBMRJLUZqp67O+Vq5pmdTRlE9UgvHbtWjz//PN48803IQgCzj//fFx99dUoLCzECy+8YMQY\niSjLiDPV+Oplkbi/N11yVdNE2UQ1CBcWFmLlypVYuXKlEeMhohxh1Ew1vmqaKNuoBuGuri7cd999\n2LJlC2w2G+bPn48f/OAHqKioMGJ8RJSlkp2pKu0n5lnClKtUg/Ddd9+NmpoarFq1CoIg4IUXXsDa\ntWvx+OOPGzE+IspyajNVpf3EAAw9S5jBnoymGoQPHTqExx57LPbv22+/HUuXLs3ooIgofyidFwzA\nkLOE9W7rk9vHAAAgAElEQVQsQqSV6l9XNBpFZ2dn7N+dnZ2IRqMZHRQR5Qel/cS7DrQr7jXW88CG\nTDYWIVKiOhP+93//dzQ0NOCCCy6AzWbD5s2bcccddxgxNiLKccr7ieX3EkudJRyfSk5GphuLEClR\nDcINDQ34/Oc/j23btkEQBFx//fWoqeE+PCJKn/J+YhdsNsjuNfa4CtDmC6C4yIlX3js4LJV83oyJ\nWDLvdE2pZC2NRVh9TZmiqa1MbW0tamv1W38hIgKU9xOfO3WoI5bUY0XuAtz79HZ09QbhctoxEDq1\nRNbZG8Sr7x1EoD+kad3YiMYi6WLBWO6SDcLLli0bcXpSvBdffDEjAyKi/KJlP3H8Y0XuAhxu88ce\niw/A8bSmkrU2FjEjEKZTMMbAnR1kg/Cdd95p5DiIKE/F7ydu9wUAmw3eMg8cdjuC4QgWz6rCkvln\noD84CI9raAasRTKpZKUbATMrp5Uqx+Vm+az0zi6yQXj27NlGjoOI8lgkGsVLm1tjgaO8xIlRHicC\nA+FhgeTCmRNl128TJZNKVmossn5TsyHbpBKlWjCWSuAm8/C2iIhMl7hFqKsvhMNt/hFbhjbtOIyK\nUm2BNZUe1WJjkfgUtFHbpBKlchKVmeNNVTAcQZsvYMmxGYHnfRGRqZQCR6Km1i7U1YzBO7uOjnjM\n7XQgFI6gvMSN82ZMwJJ5p6c9NjMrp1MpGMumSm+mzYcwCBORqZQCRyJf3wAWz6oCAOxp7kD3ySAq\n/rl+27DgTPgDYYwudqFqQhna2/vSHpuZldOpnESVDZXeIqbNh8gG4c2bNyv+4MKFC3UfDBHlH6XA\nkais2IVNO4+gqaUDPn8QZcVO1FVXxGZPRa5CXcdmxJGMSpI9icrs8WrFBimn2ARBEKQeEI8uDIVC\n+Pjjj2P7hJubm1FXV4fnn38+pRccHIygoCA/3lwi0uY3r3yMV987qPp9kyeU4uCx3hFf/z8LJuOm\nhrMzMTREIlE8+Yf/xdZ9x9HR3Y8xZR7MnT4eNy75AhwOY9KmA6FB+HqDKC91we1UTmBaYbxqjnec\nxDd+tAlS0cduA365ajHGjxll/MBMIBuERXfccQduuOEGzJgxAwDQ1NSEZ555Bo8++mhKL6hHisgs\nXm9JVo8/GflyrflynYC1r/XU+uDQjK+s2IVRnkIEBsLw9QVRXuJGXXUFmlo7JWfMlaVu3H/TnNjs\nKRPXatV9t3LXatXxAkNju+s3WzX9LkVW/vvVwustkfy66prwJ598EgvAAFBXV4fm5mb9RkZEeU9u\ni1B8IOnxB/Hu7mOSP29E0ZHakYxWY+XxZkva3AiquQmPx4Pf//73sX+/+uqr8Hg8GR0UEeWnxC1C\n8f8W146lWK3oSE96beGx2lagFYtqsLi+CpWlbthtQzPgxfVVsuvduUp1Jrxu3Tr853/+J+666y7Y\n7XZMmTIFDz74oBFjIyKKSWb2NBAaRJsvEAvMVk3LKtFrC49VtwIpNUjJJ6pBuLq6Ghs3boTfP9Sr\ntbi4OOODIiKSolYtLAacptZOtPn64XbaAdgQDEUsE3y00msLj9W3Alk5bW4E1SAsCAJefPFF/OMf\n/8D3vvc9HDlyBG1tbTj33HONGB8RUYza7Ckx4CSermSl4KNEry083Apkfaq3g+vWrcPWrVuxadMm\nAMCoUaPwwAMPZHxgRERyEteOAe2dt6zavjFeKi0rM/k8lDmqQfijjz7CI488ArfbDQAoLy9HMMhf\nHBFZi9bOW9kQfPQqQsvVYjarFZmlQzUd7XK5hp0rHI1Kn91JRGQmrZ23siH46LWFJ9e2Alm1yCwd\nqkG4trYWr776KgRBwJEjR/DrX/8as2bNMmJsRESaKQWceFYLPnJNNZJtWSlHr+exAqsXmaVCtWOW\n3+/Hj370IzQ2NgIAFi1ahNWrV2PUqNRaimV7x5NsHn8y8uVa8+U6gfy41vjq6HZfP1zOoaAmnq4k\nBp9MzZqS6VKldVan9pxaf69W7qClRcloD/7vuk1JddmykpQ7ZhUXF+P+++/XfUBERHoTq6e/scyD\n1k87DdsnrCWgJgZBrbM6vbbwZPtWIF9v9hzTmAzVIDxnzhzcdtttuO6662Jfu+666/C73/0uowMj\nIkqV21kw7ANZrw9nudmkUkBdsahmRICuq65EU2un5Gtw65C08tLsOaYxGapBePTo0Xjrrbdw6NAh\nrFmzBgBijTuIiPKB0kx3MCIo7sWNRAW8s+to7GudvUG8I9MDG8juWV0muZ0FOVVkJlJdGBk1ahSe\neOIJ9PT04LbbbkMwGBxWLU1ElOvEmW5nbxACTs10NzS2KG6N6uodwJ7mDsnH7DIfo9k8q8u0XOw3\nraljVkFBAR588EH87Gc/w8qVK9HXl9vFHUREIrWuU0vmnyGbJh1d7ES3zJ7kqExJbDbP6jItF/tN\nq86E448xvP3223H11VdDpaCaiChnqHWd6g8OYmatV/LxmVPGyDbLqChx4cJzJ+bUrC7TxCYdAEZ0\nTMtWqjPhe+65Z9i/r7jiClxxxRUZGxARkZUoNQERU8dKe3EdjhbJdcxzp3px7eJaBC/M7q1DRohE\no/jNKx/jg71Hc6ZJh0g2CD/zzDO44YYb8NBDD0k+/v3vfz9jgyIisgqtXaeWLazGl+rGAzYbvGWe\n2NfVmmVk+9YhI+Rikw6RbBB2uYZSKEVF/OMgovymFEjV9gjn4jqmkXL9JCjZIHz11VcDAG699VbD\nBkNEZEVKgXT9pmbVWVq2d6syk5aToLI5kyAbhOXS0CKmo4koU6watBJTx2qztIYFk/HKewdz6sAB\no2lZk89mskFYTEMfOnQI27dvx0UXXQQA2LRpE774xS8aMzoiyiuZOiUnU0FdbZb23FvN+GDfidjX\ncmkt0yi5dhJUItkgLKahr7/+emzcuBHl5eUAgG9+85v41re+ZczoiCiv6FGAEwxHcLzjJCLhCAoc\ntqSDejIBW2mWVlbswv5DPsmfy4W1TCOtWFSDIo8TH+w9lvUnQSVS3aLU0dERC8AAUF5ejo4O6Q4w\nRESpSrcAZ9gsui+IihIXityFONx2qs2uUlBPZRauNEub9rlybImbBcfLhbVMIznsdtzUcDYumz3J\nkssU6VANwjU1NfjBD36AK6+8EgCwceNG1NRk/90HEVlLugU4UrNoqRkqIB3UU52Fy1VONyw4EwcO\n+XJ2LdMMubidSzUIP/DAA/j5z3+O++67D4IgYO7cubjzzjuNGBsR5ZF0CnCUZtFSEoN6OrNwpcrp\nXF7LJH0oBuFIJILXXnuNQZeIMi6dAhylWbSUxKCuxzYYqVmaWqOOZFm1atyq48oGikHY4XBgw4YN\nWLFihVHjIaI8lmrQUppFS0kM6qnMwrUEnlQadUg9r9J6tZkyVc2eT1TT0XPmzMEbb7yBSy+91Ijx\nEFEeS7W7lNIsetLYYgQGBhWDejKz8FQLuNRm0krPq7Re/a1rZg17HiNnpbncTtIoqkH45ZdfxlNP\nPQW32w2PxwNBEGCz2bBlyxYjxkdEeSiVAhylWfRgRFANTFpn4ZkKPHLPG4lE0dTaKfkzu5s7MBAa\nBGD8rDTX20kaRTUIv/TSS0aMg4goLfGzaIezEJFQOBYEHHaoBnUts/BMBR7F5/2kAz3+kORjvr4B\n+HqDKIDxs9JcbydpFNUgPHHiRCPGQUSkC1ehA94xo9De3pfyz8sFj0wFHqXn7fGHUFbsgs8vvV5d\nXupCR4ff8FlprreTNIpqED5+/Dgefvhh7N+/H8HgqTf77bffzujAiIisJlOBR+l5K0rdqKupxDu7\njo54bGbtGLidBabMSnOpnaSZ1d2qQXjNmjW4/PLL8be//Q2PPPIInnvuOZx++ulGjI2IyFIyFXjU\nnndoXdcmu15t1qxU7y1YRrNCdbdqEPb5fLjqqqvw29/+FjNnzsSMGTOwYsUKHnFIRFkpmVmP1Pdm\nKvDEP29X7wBGFzsxc8qYWEBQWq82a1aa7WclW6G6WzUIFxYWAhg6VenYsWMYM2YMurq6Mj4wIiI9\nJTPrUfveTAQeh92OFYtqEIkK2NPcgW5/EE2tnXA4WmKvq7Renc7NQbrp2GxsJ2mV6m7VIFxfX4/u\n7m5cc801uOKKK+B0OnHJJZdkfGBERHpS2gJ0yezThwUgLTOkTASeDY0tw9Z+k5mZpXJzYIV0rFms\nUt2tGoTFlpUNDQ2YPXs2/H4/amu5CZuIsofSrGfznmN4d/exWABqWHCmKTMkvWZmydwcWCEdaxar\nVHfLBuGWlhbJr9vtdrS0tPAkJSLKGkqznqgw9P/FABQYGDRlhmT0zMwq6Vi9JJtSt0p1t2wQvvnm\nm2Gz2SAIAo4fP47i4mIAgN/vx/jx49HY2GjIAImI0pVMb+n9//CZMkMyemZmlXRsutJJqVuhuls2\nCItB9r777kN9fT0uu+wyAMAbb7yBHTt2GDM6IiIdKM16EnX7g5j3hXH4YN+JEY9lcoZk9MzMKunY\ndKWTUrdCdbfqmvD27dvxX//1X7F/X3rppfjv//7vlF+wvLwIBQXZk+JI5PWWmD0Ew+TLtebLdQL5\nfa23Lp+JIo8TW/cdR7uvHzY7EI2O/LkxZR7cdvVMVL55AFv3HUdHdz/GlHkwd/p43LjkC3A4Mlew\nFD/GZF431d/reTMm4tX3Dkp8fQKqJpSl9JxqBkKD8PUGUV7qgtupGoKGSbzOgdCgbF/tptZOfGOZ\nR/NrVCU1Ev2ojk4QBOzYsQP19fUAgB07diAq9Zerkc8XSPlnzeb1lqTcCi/b5Mu15st1ArxWAGg4\n7wxcNnsSevxBvLn9sGQXqrrqSgT8wdj3tvsCgM0Gb5kHXV0nMz72+DGKMzOl103n97pk3ukI9IdG\npGOXzDtd97+VdCuxpa6zzRdAu69f8vs7uvvR+mmnZVLqcjdKqkH47rvvxh133AGPxwMACAaDePTR\nR/UdHRGRQcTq4WsXT1HsQhWJRvHS5lZTtu8Yte/WyHRsJiqxcyGlrmmf8KZNm/D3v/8dgiBg8uTJ\ncDqdRoyNiChj1AJQNm7fSbXpRqaDfqYqsa1S4ZwOTcnySCQCp9OJSCSCQ4cOAQC3KBFRTpAKQOkG\nDaMPBLB6041MVmJbocI5HapB+H/+53/wyCOPoKysDDabDQBgs9l4ihIR5axUg4ZZwdDqs/ZMpo2t\nUOGcDtUg/OSTT+KPf/wjzxUmoryRatAwIxgOhAYt33TDiLRxNvavBgDVWzOv18sATER5RQwaUuSC\nhloKOxiO6DpGka9XfdaeimA4gjZfQLdxr1hUg8X1VagsdcNuAypL3VhcX5U1aeNMUZ0Jz58/Hw89\n9BD+5V/+BS7Xqbs/rgkTUS5Ldq3RrA5U5aX6pnozlVI3Om1s9Lp8qlSD8CuvvAJgqFOWiGvCRJTr\nkg0aZm2XcTsLdE31yqXUAwODWHnJ1LQDWqbTxlYvUkukGoTZI5qI8pnWoGHmdhm9KoSVUuof7juB\nA4d8mgKambNQLevyVpola+4Z1tnZiWDw1B3ehAkTMjIgIqJsZdZ2Gb1SvUopdUC90MzsWajaunzD\ngsl45b2DlpolqwbhLVu2YNWqVejs7ITdbkc4HEZZWRm2bNlixPiIiLKG2dtl0k31aj1tSq7q2uyt\nUmrr8s+91TzsYA4rbOVSDf0PP/wwnn76adTU1GDv3r245557sHz5ciPGRkSUlcRgaHaqM1lKVeHx\npKqulWahuw60Z6w6PJ54EyGlrNiF/Yd8ko9lsnpdjab595lnnonBwUHYbDYsX74c7733XqbHRURE\nJhC3ElWUyBeSSRWaKc1Cu/qC+N2bBxBJ4/AfLZRuIqZ9rjwjW7nSpRqECwqGMtannXYaGhsbceDA\nAfT09GR8YEREZKxgOILOngEsW1iN/3fzXJw3fZzk90kVminNQgHgg30nsKGxRdfxSpHbj3ztRVNk\nx2fmYQ+qa8LXX389enp68K1vfQvf/e530dfXh9WrVxsxNiIiUhAMR3C84yQi4UhaqW+5gqqVl9bC\n4y7QVGimVB0uMqKDl9K6vBUPe1ANwhdccAGKi4tRV1eHt956CwDg9/szPjAiIpI2LGj2BVFRkl6V\nr1pBldZCsxWLahAYGMSHccVP8TLZtCSRVJGaFQ97UA3CK1euxMsvv6z6NSIiMoaeVchaT4zSEjgd\ndjtWXjIVBw75LHnGr9nV61Jkb5kGBwfR39+PaDSKgYEB9Pf3o7+/H21tbejv7zdyjEREGaF3f2Qj\n6N2jWku7zWSk0nfbaFaqXpedCf/yl7/E448/DpvNhnPOOSf29eLiYvzbv/2bIYMjIsoEvZpKmNF5\nSe8e1Zlot2nFtK9VyQbhW2+9FbfeeivuvfderF271sgxERFlVLrpXDM7Q+kVNONvIPQuWLJi2teq\nVNeExQDc29uLbdu2oaqqCtOmTcv4wIiIMkHrGqgSMztDpdujWuoGYsaUMfjyrInY80mnrjPXbD3j\n10iyQfh73/sevv71r2PatGno7u7G0qVLUVxcDJ/Ph+985zu46qqrjBwnEZEu0k3n6hHE001jp5Pu\nlbqBaNx5FIvrq3D/TXM4czWYbBD+61//Gpvx/v73v0d1dTWefPJJnDhxAt/4xjcYhIkoK6Wbzk0n\niOuVxo5P9zqchYiEwpqCppYbiGycuVrpVKRkyQZhl+vUH+LOnTuxePFiAMC4ceNgs9kyPzIiogxI\nN52bThDXO43tKnTAO2YU2tv7NH2/3kVdZkvlpsZqAVtxTfizzz7D6NGjsW3bNtx+++2xr8cfaUhE\nlG3SSeemGsT1SGOny+MqQFmxCz6JbUdG7OHVOwAmc1Nj9jGLcmSD8M0334yGhgYUFhZi1qxZqKkZ\n+uPcs2cPzxImoqyWbvVuKkHczFlofACSCsBAZvfwZiIAJntTY/Yxi3Jkg/Bll12G+vp6dHR0DKuG\nHj9+PO677z5DBkdElEmpVu+mEsQzsR9Xq8QAFK+yNPN7eDMRAJO5qbFCFkKO4i2I1+vFWWedNWwN\n+LTTTuNMmIgIyp2XErtxmdVJSikAlRU7sfZr9bh2cW3GUrJ6d/gSKZ3alHhTo3dXMD2p7hMmIiLt\nlFKvZnSSUgpAvSdD6A8OoqTIacrrp5OGT2Zt3swshBoGYSIiHamlXo3uJGV2AMrk62u9qUm3Ij6T\nGISJiHSide3RyE5SZgegdKrJ1W5Uklmbt2o/awZhIiKdWHUfrtkBKJnXl0vn37p8puzza7mpsWo/\nawZhIiKdmJ36lWN2AErm9eXS+UUeJxrOOyPtsVitn7V5O5SJiHKM1c/SNfscXbXXV0rnb913PKvO\nfdaKQZiI8k7i9iE9rVhUg8X1VagsdcNuG9qHu7i+ypS1x0xeZyZeQymd39Hdn9ZWIiPei1QwHU1E\necOI1oVmp34BY64zEBzEc281Y/8hn26voZTOH1PmSSmdb9V2lSIGYSLKG0a2LjRz7TGT1ykGtfeb\njmEgFNX1NZQqqedOH5/SzYxV21WKzL8NICIyQKY6N1lNpq9TDGrxAVjP15BL59+45AtJP1c2/M45\nEyaivGDV7UN6S7ancjIpc6WgJvcayZJL5zscyc8Zs+F3ziBMRHnBqtuH9KblOlNdJ1UKaqKyYhdC\ng1EEw5G01sL1SOdnw+/c8CBcXl6EggLzN0inyustMXsIhsmXa82X6wR4refNmIhX3zso8fUJqJpQ\nZsSwMiLxWtWu8zevfCy7F/emhrNlX6dktAfecg/afP2y39MfiuDuJ7fBW+bB3OnjceOSL6Q0i5WS\nyt+v1X/nhgdhny9g9EvqxustQXt7n9nDMES+XGu+XCfAawWAJfNOR6A/NKJz05J5p2fteyN1rUrX\neeRYNz7Ye1TyuT7YewyXzZ4EALJp6rrqSsnCKYcdiESB/uAgAKDN149X3zuIQH9IlwKoVP9+rfI7\nl7uBsAmCIBg2CiBr/9ABfojlony5ToDXGi/ZtVArU7pWqets8wWw+ldbIfXBbwMwf/o4xW1Hp1LZ\nQ0GtrNiFKZPK8MmRbslUdWWpG/ffNCft9zndv1+zf+dyQZhrwkSUd6zWujBTpK5TaZ3U5XTgg30n\nYv+W2s4jVTjV4w9i9a+2So7BKgVQVv2dc4sSEVEeUWqtKUdqO098C0oxsEuxSgGUVTEIExHlGam9\nuPOnj0MwJL1vVpzNyrF6z2wrYzqaiCjPSKWUAeDAIV/K23nMPi4xWzEIExHlqcR1UrmWkVpms1bo\nmZ2NGISJiAiAPrNZqxZAWRWDMBERAVCezZq9xSdXMQgTEdEw8bNZqx8FmO0YhImISJbVjwLMdryN\nISIiSUqnJu060G6JowCzHYMwEZEJguEI2nwBBMORYf/bSpROTerqC+J3bx5AJCp9rjBpw3Q0EZGB\n4tdYO3uDcDvtAGwIhiKWW29VanEJAB/sOwGPu4Bp6TSY/1smIsoj4hqrGNgGQlEMhCIQcGq9dUNj\ni7mD/CctLS6lWlqSdgzCREQGUVpjjWelwLZiUQ3mTx8n+7haS8tMsmoaPxlMRxMRGURpjTWeVU4e\nAob2Dq+8ZGpaLS31lkvbprJrtEREWUzptKF4Vjt5yGoHNMSn9KXS+Nk0Q+ZMmIjIIGIwk+rPHM+K\nJw9Z5YAGtW1TkaiAppaOrJkhMwgTERkoPph19Q7A5RwKtqFwxNInD1nlgAa1bVPv7Doa+3c2NBZh\nECYiMpDcMYLZ0pfZrAMaBkKDaPMF4HEVyG6bstuAqDDyZ3c3d2DZwmpLvrcMwkREJkgMZlYowrKC\nxIMixCKsptZOtPv6UVHqQpG7UDIISwVgwFqFbokYhImIyHRyFc+CIODtncNTzJ29QUwaW4zAwGBs\nfbquugJNrZ1pVXCbcVIUgzAREZlO7qAIt1M6GAYGBrH2a/XoDw7Ggub6Tc2SRW9qhW5mbnliECYi\nIlMpVTwPhKS3Gfn6BtAfHByWYk61gtvMk6IYhImIyFRam5jEk0oxp1LBrXQDYERBlzU3ThERUd5Q\namIil45WSjGLRW9agqfSDYARLTkZhImIyFRKHbnmnz0Oi+urMLbcA7sNqCx1Y3F9lW57qZVuAIzo\nXMZ0NBERmU5pPddht+Mbyzxo/bRT98plpS5mRnQuYxAmIspiYhOLbGj0oURtPdftLMjYPl8zW3Iy\nCBMRZSGpJhZW75OshRkducxsyckgTESUhczcVpOrzLgByN7bJSKiPKW2rSYbjvCjIQzCRERZxuxt\nNaQfBmEioixj9rYa0g+DMBFRllHaV2vEthrSDwuziIiykLh9pqm1Ex3d/YZuqyH9MAgTEWUhcVtN\npppYkDEYhImIslgmm1hQ5nFNmIiITBcMR9DmC+Td9irOhImIyDRi56/dze3o6g3mTOcvrRiEiYjI\nNGZ3/gqGI4a3qoxnEwRBMPIFBwcjKChg8QARUb4bCA3iPx5qRJuvf8RjY8s9+Pn3F8HtzMxcMRKJ\n4sk//C+27juO9u5+eMs8mDt9PG5c8gU4HMbNwA2fCft8AaNfUjdebwna2/vMHoYh8uVa8+U6AV5r\nrsrma23zBdAuEYABoKO7H62fdsaKzvS+zvWbmofNwNt8/Xj1vYMI9IcyMgP3ekskv577CXciIrIk\nszp/Wan3NoMwERGZwqzOX1bqvc3CLCIiMo3Y4Wt3cwd8fQOGdP4SZ+CdEoHY6N7bDMJERGQasfPX\nsoXVI6qU4yuX9STOwOPXhEVG995mECYiItO5Ch2xIiypvcPnzZiIJfNO123vsBkzcCkMwkREZClS\ne4f1rlxWmoEbiYVZRERkGXpXLqu1wxRn4GYdfsGZMBERWYaWymUtB1ZkSztM64yEiIjynl57h8WU\ndmdvEAJOtcPc0Nii42jTxyBMRESWocfeYSs141DDdDQREVmKVOXyeTMmYMm80zX9vF4pbSMwCBMR\nkSnkTjCSqlyumlCmuXe0lZpxqGEQJiIiQ2ktmorfO5wMKzXjUMMgTEREhjLiDGGrNONQwyBMRJRF\nzD6EPl1qRVPLFlbrcl1WacahhkGYiCgLyKVwb10+0+yhJcXooqlUU9pGYRAmIsoCcincIo8TDeed\nYd7AkpRNRVNG4D5hIiKLU0rhbt133FL7XtWYdYawVXEmTERkcUop3I7ufkvte9UiW4qmjMAgTERk\ncUop3DFlnqxL4WZL0ZQRmI4mIrI4pRTu3OnjszaAmX2CkRVwJkxElAXkUrg3LvkCurpOmjw6ShWD\nMBFRFpBL4Toc1k9oZvve5kxiECYiyiJW3/caL1vO9DUTgzAREWWEEe0psx1vRYiISHfZdKavmRiE\niYhId1raUxKDMBERZYC4t1lKPranlMMgTEREumN7Sm1YmEVERBnB9pTqGISJiCgj2J5SHYMwERFl\nVDbtbTYa14SJiIhMwiBMRERkEgZhIiIikzAIExERmYRBmIiIyCQMwkRERCZhECYiIjIJgzAREZFJ\nbIIgCEa+4OBgBAUF7JhCRERkeMcsny9g9EvqxustQXt7n9nDMES+XGu+XCfAa81V+XKt2X6dXm+J\n5NcNnwkTERHREK4JExERmYRBmIiIyCQMwkRERCZhECYiIjIJgzAREZFJGISJiIhMwiBMRERkEgZh\nIiIikzAIExERmYRBmIiIyCQMwkQS/vSnP6GhoQFLly7FpZdeiu9+97uxxx577DGEQiHdX3PRokVo\nbm5O+3kee+wxPPjggwCA5557Dk8//XTazylauXIl3nnnHd2ejyjfGX6AA5HVtbW14Z577sHLL7+M\n8ePHQxAE7N+/P/b4448/jhtvvBFOp9PEUQ4ZHBxEQYH8f8bXXHONgaM5JRKJwOHgaWlEahiEiRJ0\ndHSgoKAAZWVlAACbzYazzjoLAHDPPfcAAK6++mrY7XY8++yz2Lx5M377298iHA4DAO68807MmzcP\nwNDsdunSpfjwww/R3t6OG2+8Eddddx0AYMeOHbHn++IXv4j4s1QefPBBbNu2DeFwGOXl5XjggQcw\ncf1bct4AAAWXSURBVOJEHDlyBMuWLcMVV1yBrVu3Yvny5fjXf/1X/OAHP0BzczO8Xi/GjRuHMWPG\nABiaFQcCAdx5551Yu3Yt9u7dCwAIBALo7u7G9u3b4ff7sW7dOhw4cADBYBBz5szB6tWr4XA40NLS\ngtWrVyMQCKC2thbBYFDyPfvoo49w//33Y/r06fjrX/+Kb3/72/D7/Wm/L3PmzMHbb7+NX/3qV6it\nrcXBgwfxwAMPwOfzIRwO44YbbsCyZcvS/p0TmUYgomEikYjwzW9+U5g9e7Zw2223CU899ZTQ1dUV\ne7y2tlbw+/2xf3d1dQnRaFQQBEFobW0VFixYEHvswgsvFH70ox8JgiAIhw8fFs455xzB7/cLwWBQ\nOP/884WtW7cKgiAIr732mlBbWyscOHBAEARB6OzsjD3HCy+8IHz729+OPUdtba3w2muvxR5ft26d\nsGrVqtjPLVy4MPaaP/vZz2L/WxQKhYSVK1cKTz/9tCAIgrBmzRrh5Zdfjl37d77zHWHDhg2CIAjC\nV77yFWHjxo2CIAjC7t27hWnTpgmNjY0j3rOtW7cK06ZNE3bt2pX2+7JgwQJh+/btgiAIwp///OfY\n+xIOh4WvfOUrQktLiyAIgtDX1ydcfPHFsX8TZSPOhIkS2O12/OIXv0BzczO2b9+OTZs24YknnsAf\n/vCH2Ow43uHDh/Hd734Xn332GQoKCtDR0YH29nZ4vV4AwOWXXw4AqKqqQmlpKU6cOIFwOAyPx4M5\nc+bEvmft2rWx5/zLX/6C9evXIxAIYHBwcNjruVwuXHbZZbF/f/TRR7jrrrsAABUVFbjooosUr++u\nu+5CbW0tbrjhBgBAY2Mjmpqa8NRTTwEABgYGcNppp8Hv96O5uRlLly4FAJxzzjmora2Vfd7Pfe5z\nmDlzZtrvi9vtRn19PQDgoosuQmlpKQDg008/RWtrK+64447Ya4TDYRw8eBDV1dWK10xkVQzCRDJq\na2tRW1uLr371q7j88suxbds2XHzxxSO+74477sCqVauwePFiRKNRzJgxY1ja1uVyxf63w+FAJBJR\nfN2jR49i3bp1ePHFFzFp0iTs2rUL3/ve92KPezwe2Gy22L+FJI4E/+lPf4q+vj6sW7du2M//4he/\nwKRJk4Z9r9/vH/Y6aoqKiob9O9X3Re41BUFAeXk5fv/732seE5HVsTqaKMFnn32G3bt3x/594sQJ\ndHV1oaqqCgAwatQo+P3+2ON9fX2xx1588UVNldOTJ0/GwMAAtm/fDgB444030NfXB2Ao+BUWFsLr\n9SIajeL5559XfK558+Zh48aNAACfz4dNmzZJft/GjRvx/vvv49FHH4Xdfuo//UWLFuHXv/51LAh2\ndXXh8OHDKC4uxpQpU/CHP/wBANDU1JRU9Xaq70sgEMDOnTsBAJs2bUJvby8A4Mwzz4Tb7cYrr7wS\n+/7W1tZhvwuibMOZMFGCwcFBPPbYYzh69Cjcbjei0Si+/e1v4/Of/zwA4MYbb8T1118Pt9uNZ599\nFqtXr8Ytt9yC0aNHY8GCBZIp60ROpxM//vGPhxVmTZgwAQAwdepUXHrppbj88stRXl6OhQsXYseO\nHbLPdcstt2DNmjW49NJL4fV6Y6ncRI8//jiAoaIyYOhmYv369VizZg0efvhhLF26FDabDYWFhViz\nZg0mTZqEhx56CKtXr8ZvfvMb1NbW4uyzz9b8Pqb6vjz66KP44Q9/CACYPXs2KisrUVJSgoKCAvzy\nl7/EAw88gCeeeALRaBSVlZX4yU9+onlMRFZjE5LJZRERZZjf70dxcTEAYOvWrVi1ahUaGxuHzd6J\ncgVnwkRkKX/+85/x9NNPQxCE2MyYAZhyFWfCREREJuHtJRERkUkYhImIiEzCIExERGQSBmEiIiKT\nMAgTERGZhEGYiIjIJP8fpX1q9LE28DkAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(8, 6))\n", "\n", "ax.scatter(df.std_range, df.std_logratio,\n", " c=blue);\n", "\n", "ax.set_xticklabels([]);\n", "ax.set_xlabel(\"Standardized range\");\n", "\n", "ax.set_yticklabels([]);\n", "ax.set_ylabel(\"Standardized log ratio\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This data set has a two interesting properties that make it useful for illustrating dependent density regression.\n", "\n", "1. The relationship between range and log ratio is nonlinear, but has locally linear components.\n", "2. The observation noise is [heteroskedastic](https://en.wikipedia.org/wiki/Heteroscedasticity); that is, the magnitude of the variance varies with the range.\n", "\n", "The intuitive idea behind dependent density regression is to reduce the problem to many (related) density estimates, conditioned on fixed values of the predictors. The following animation illustrates this intuition." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n", " (prop.get_family(), self.defaultFamily[fontext]))\n" ] } ], "source": [ "fig, (scatter_ax, hist_ax) = plt.subplots(ncols=2, figsize=(16, 6))\n", "\n", "scatter_ax.scatter(df.std_range, df.std_logratio,\n", " c=blue, zorder=2);\n", "\n", "scatter_ax.set_xticklabels([]);\n", "scatter_ax.set_xlabel(\"Standardized range\");\n", "\n", "scatter_ax.set_yticklabels([]);\n", "scatter_ax.set_ylabel(\"Standardized log ratio\");\n", "\n", "bins = np.linspace(df.std_range.min(), df.std_range.max(), 25)\n", "\n", "hist_ax.hist(df.std_logratio, bins=bins,\n", " color='k', lw=0, alpha=0.25,\n", " label=\"All data\");\n", "\n", "hist_ax.set_xticklabels([]);\n", "hist_ax.set_xlabel(\"Standardized log ratio\");\n", "\n", "hist_ax.set_yticklabels([]);\n", "hist_ax.set_ylabel(\"Frequency\");\n", "\n", "hist_ax.legend(loc=2);\n", "\n", "endpoints = np.linspace(1.05 * df.std_range.min(), 1.05 * df.std_range.max(), 15)\n", "\n", "frame_artists = []\n", "\n", "for low, high in zip(endpoints[:-1], endpoints[2:]):\n", " interval = scatter_ax.axvspan(low, high,\n", " color='k', alpha=0.5, lw=0, zorder=1);\n", " *_, bars = hist_ax.hist(df[df.std_range.between(low, high)].std_logratio,\n", " bins=bins,\n", " color='k', lw=0, alpha=0.5);\n", " \n", " frame_artists.append((interval,) + tuple(bars))\n", " \n", "animation = ani.ArtistAnimation(fig, frame_artists,\n", " interval=500, repeat_delay=3000, blit=True)\n", "plt.close(); # prevent the intermediate figure from showing" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n", " (prop.get_family(), self.defaultFamily[fontext]))\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "HTML(animation.to_html5_video())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we slice the data with a window sliding along the x-axis in the left plot, the empirical distribution of the y-values of the points in the window varies in the right plot. An important aspect of this approach is that the density estimates that correspond to close values of the predictor are similar.\n", "\n", "In the previous example, we saw that a Dirichlet process estimates a probability density as a mixture model with infinitely many components. In the case of normal component distributions,\n", "\n", "$$y \\sim \\sum_{i = 1}^{\\infty} w_i \\cdot N(\\mu_i, \\tau_i^{-1}),$$\n", "\n", "where the mixture weights, $w_1, w_2, \\ldots$, are generated by a [stick-breaking process](https://en.wikipedia.org/wiki/Dirichlet_process#The_stick-breaking_process).\n", "\n", "Dependent density regression generalizes this representation of the Dirichlet process mixture model by allowing the mixture weights and component means to vary conditioned on the value of the predictor, $x$. That is,\n", "\n", "$$y\\ |\\ x \\sim \\sum_{i = 1}^{\\infty} w_i\\ |\\ x \\cdot N(\\mu_i\\ |\\ x, \\tau_i^{-1}).$$\n", "\n", "In this example, we will follow Chapter 23 of [_Bayesian Data Analysis_](http://www.stat.columbia.edu/~gelman/book/) and use a probit stick-breaking process to determine the conditional mixture weights, $w_i\\ |\\ x$. The probit stick-breaking process starts by defining\n", "\n", "$$v_i\\ |\\ x = \\Phi(\\alpha_i + \\beta_i x),$$\n", "\n", "where $\\Phi$ is the cumulative distribution function of the standard normal distribution. We then obtain $w_i\\ |\\ x$ by applying the stick breaking process to $v_i\\ |\\ x$. That is,\n", "\n", "$$w_i\\ |\\ x = v_i\\ |\\ x \\cdot \\prod_{j = 1}^{i - 1} (1 - v_j\\ |\\ x).$$\n", "\n", "For the LIDAR data set, we use independent normal priors $\\alpha_i \\sim N(0, 5^2)$ and $\\beta_i \\sim N(0, 5^2)$. We now express this this model for the conditional mixture weights using `pymc3`." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def norm_cdf(z):\n", " return 0.5 * (1 + tt.erf(z / np.sqrt(2)))\n", "\n", "def stick_breaking(v):\n", " return v * tt.concatenate([tt.ones_like(v[:, :1]),\n", " tt.extra_ops.cumprod(1 - v, axis=1)[:, :-1]],\n", " axis=1)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [], "source": [ "N, _ = df.shape\n", "K = 20\n", "\n", "std_range = df.std_range.values[:, np.newaxis]\n", "std_logratio = df.std_logratio.values[:, np.newaxis]\n", "\n", "x_lidar = shared(std_range, broadcastable=(False, True))\n", "\n", "with pm.Model() as model:\n", " alpha = pm.Normal('alpha', 0., 5., shape=K)\n", " beta = pm.Normal('beta', 0., 5., shape=K)\n", " v = norm_cdf(alpha + beta * x_lidar)\n", " w = pm.Deterministic('w', stick_breaking(v))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have defined `x_lidar` as a `theano` [`shared`](http://deeplearning.net/software/theano/library/compile/shared.html) variable in order to use `pymc3`'s posterior prediction capabilities later.\n", "\n", "While the dependent density regression model theoretically has infinitely many components, we must truncate the model to finitely many components (in this case, twenty) in order to express it using `pymc3`. After sampling from the model, we will verify that truncation did not unduly influence our results.\n", "\n", "Since the LIDAR data seems to have several linear components, we use the linear models\n", "\n", "$$\n", "\\begin{align*}\n", "\\mu_i\\ |\\ x\n", " & \\sim \\gamma_i + \\delta_i x \\\\\n", "\\gamma_i\n", " & \\sim N(0, 10^2) \\\\\n", "\\delta_i\n", " & \\sim N(0, 10^2)\n", "\\end{align*}\n", "$$\n", "\n", "for the conditional component means." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "with model:\n", " gamma = pm.Normal('gamma', 0., 10., shape=K)\n", " delta = pm.Normal('delta', 0., 10., shape=K)\n", " mu = pm.Deterministic('mu', gamma + delta * x_lidar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we place the prior $\\tau_i \\sim \\textrm{Gamma}(1, 1)$ on the component precisions." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "with model:\n", " tau = pm.Gamma('tau', 1., 1., shape=K)\n", " obs = pm.NormalMixture('obs', w, mu, tau=tau, observed=std_logratio)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now sample from the dependent density regression model." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 20000/20000 [01:42<00:00, 195.86it/s]\n" ] } ], "source": [ "SAMPLES = 20000\n", "BURN = 10000\n", "THIN = 10\n", "\n", "with model:\n", " step = pm.Metropolis()\n", " trace_ = pm.sample(SAMPLES, step, random_seed=SEED)\n", " \n", "trace = trace_[BURN::THIN]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To verify that truncation did not unduly influence our results, we plot the largest posterior expected mixture weight for each component. (In this model, each point has a mixture weight for each component, so we plot the maximum mixture weight for each component across all data points in order to judge if the component exerts any influence on the posterior.)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n", " (prop.get_family(), self.defaultFamily[fontext]))\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAFzCAYAAAA0dtAgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8zvXj//HntZPDzPHDJh+ffZxPE3LIUsYQGnMYMZJu\nSJ+VHAqZT1TIfFJpUXzlE3IspxyW4kPRQeSXjGyfIgvVNrHYCW27fn+4uT7WdnnPdl27dvV+3G83\nt9ve7+t1vTx31fb0PlusVqtVAADgT8/D1QEAAEDpoPQBADAJSh8AAJOg9AEAMAlKHwAAk6D0AQAw\nCS9XB3C28+fTizy2WrWKSkvLcmIax3K3vBKZS4O75ZXIXBrcLa9E5uKqWdPP7mts6d/Ey8vT1RFu\ni7vllchcGtwtr0Tm0uBueSUyOwOlDwCASVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCA\nSVD6AACYBKUPAIBJUPoAAJhEmSn96OhoBQcHq0+fPoW+brVaNWfOHPXo0UN9+/bVt99+W8oJAQBw\nb2Wm9AcOHKhly5bZfX3//v1KSkrSrl27NHv2bD3//POlFw4AgD+BMlP67du3V5UqVey+vmfPHvXv\n318Wi0WtW7fW5cuXlZqaWooJAQBwb27zaN2UlBQFBATYlgMCApSSkqJatWrd8n3VqlWUl5en+j69\n1SE5tr/SzyHzOMqtHqFoj6s/i+JkdjV3y+xueSUylwZ3yyuR2dHcpvStVmuBdRaLxfB9jn6u8fnz\n6Q6dryRq1vRzaZ7i/N2uzlwc7pbZ3fJKZC4N7pZXInNJMthTZnbvGwkICFBycrJtOTk52XArHwAA\n/I/blH5oaKjef/99Wa1WffPNN/Lz86P0AQC4DWVm9/5TTz2lQ4cOKS0tTZ07d9aTTz6pnJwcSVJk\nZKRCQkK0b98+9ejRQxUqVNDcuXNdnBgAAPdSZkr/1VdfveXrFotFzz33XCmlAQDgz6fMlD4KGjVv\nb4nneHtaqAOSAAD+DNzmmD4AACgZSh8AAJOg9AEAMAlKHwAAk6D0AQAwCUofAACToPQBADAJSh8A\nAJOg9AEAMAlKHwAAk6D0AQAwCUofAACToPQBADAJSh8AAJOg9AEAMAlKHwAAk6D0AQAwCUofAACT\noPQBADAJSh8AAJOg9AEAMAlKHwAAk6D0AQAwCUofAACToPQBADAJSh8AAJOg9AEAMAlKHwAAk6D0\nAQAwCUofAACToPQBADAJSh8AAJOg9AEAMAlKHwAAk6D0AQAwCUofAACToPQBADAJSh8AAJOg9AEA\nMAlKHwAAk6D0AQAwCUofAACToPQBADAJSh8AAJOg9AEAMAlKHwAAk6D0AQAwCUofAACToPQBADAJ\nSh8AAJOg9AEAMAlKHwAAkyhTpb9//3717NlTPXr00NKlSwu8/vPPP2vEiBHq37+/+vbtq3379rkg\nJQAA7snL1QFuyM3N1axZs7R8+XL5+/tr0KBBCg0NVcOGDW1jFi9erN69e2vYsGE6efKkxo4dq717\n97owNQAA7qPMbOnHx8crMDBQdevWlY+Pj8LCwrRnz558YywWizIyMiRJ6enpqlWrliuiAgDglsrM\nln5KSooCAgJsy/7+/oqPj883Zty4cRo9erRWr16t7OxsLV++3HDeatUqysvL02E5a9b0c9hcpcGZ\neYs7t7t9hpL7ZXa3vBKZS4O75ZXI7GhlpvStVmuBdRaLJd9yXFycBgwYoFGjRunIkSOaOnWqduzY\nIQ8P+zss0tKyHJrz/Pl0h87nbM7MW5y5a9b0c7vP0N0yu1teicylwd3ySmQuSQZ7yszu/YCAACUn\nJ9uWU1JSCuy+37hxo3r37i1JatOmja5evaq0tLRSzQkAgLsqM6XfsmVLJSUl6ezZs7p27Zri4uIU\nGhqab0zt2rV14MABSdKpU6d09epVVa9e3RVxAQBwO2Vm976Xl5dmzpypMWPGKDc3VxEREWrUqJFi\nY2MVFBSkbt26adq0aXr22We1YsUKWSwWzZs3r8AhAAAAULgyU/qSFBISopCQkHzrJkyYYPu6YcOG\nWr9+fWnHAgDgT6HM7N4HAADORekDAGASlD4AACZB6QMAYBKUPgAAJkHpAwBgEnYv2cvOzr7lGytU\nqODwMAAAwHnsln6bNm1ueeObhIQEpwQCAADOYbf0ExMTJUlvvvmmfHx8NGTIEFmtVm3YsEG///57\nqQUEAACOYXhMf/fu3RozZoz8/PxUuXJljR49Wrt27SqNbAAAwIEMS//KlSv68ccfbctnzpwxPN4P\nAADKHsN770+aNEkPPviggoKCJEknTpzQ7NmznR4MAAA4lmHp33///WrXrp2++eYbWa1WtWnThsfZ\nAgDghop0nf6lS5eUl5enbt26qVy5cvrtt9+cnQsAADiYYelv2bJFUVFRiomJkSSlpKRo4sSJTg8G\nAAAcy7D0V65cqU2bNsnPz0+SVL9+ff36669ODwYAABzLsPS9vb3l6+ubb52np6fTAgEAAOcwLP2q\nVavq9OnTtrvzbd26VQEBAU4PBgAAHMvw7P3o6Gg9/fTTOn36tEJDQ1W+fHktWbKkNLIBAAAHMiz9\n+vXra8OGDUpKSpLValW9evXYvQ8AgBsy3L0/YcIEeXp6qkGDBmrYsKE8PT01YcKE0sgGAAAcyLD0\nz5w5U2DdDz/84JQwAADAeezu3n/vvff07rvvKikpSYMGDbKtT09PV7169UolHAAAcBy7pd+pUycF\nBgZq9uzZmjp1qm19pUqV1KRJk1IJBwAAHMdu6depU0d16tTR9u3bbZfrAQAA92V4TH/YsGG6dOmS\nbfm3337T8OHDnRoKAAA4nmHpZ2VlqUqVKrblqlWrKiMjw6mhAACA4xmWfl5enrKysmzLmZmZys3N\ndWooAADgeIY35+nTp49GjRqlyMhISdK6desUHh7u9GAAAMCxDEv/scceU61atbR3715ZrVYNHTpU\n/fv3L41sAADAgQxLX5IGDBigAQMGODsLAABwIsNj+qdPn1ZkZKRCQ0MlSd9++60WLlzo9GAAAMCx\nDEv/hRdeUFRUlPz8/CRJzZo104cffuj0YAAAwLEMSz89PV2dO3e23aDHw8ND3t7eTg8GAAAcy7D0\nPT099fvvv9tKPyUlRR4ehm8DAABlTJHuyDdu3DilpaVp4cKFGjZsmEaNGlUa2QAAgAMZnr3fv39/\n/fWvf9XHH3+s7Oxs/etf/1K7du1KIxsAAHCgIl2y165dOzVr1kyS5Ovr69RAAADAOQx37586dUoR\nEREKDg5WcHCwBg0apFOnTpVGNgAA4ECGpR8dHa0RI0bo6NGjOnr0qEaMGKHo6Gi741988cUirQMA\nAKWrSE/Z69+/vywWiywWi/r166fs7Gy74w8fPlxg3VdffVWylAAAoMQMj+m3aNFChw8ftp289//+\n3/9TUFBQgXE7d+7Uzp079dNPP2nChAm29RkZGSpfvrwDIwMAgOIwLP3ExESNGDFCf/vb3yRJZ86c\nUePGjTVo0CBJ0saNGyVJ9erVU5cuXXTs2DF16dLF9v5KlSopODjYCdEBAMDtMCz9f/7zn0WaqGnT\npmratKlCQ0NVtWrVEgcDAACOZVj6QUFBqlixYr51KSkp8vf3L3R8bm6uXnvtNZ09e1Y5OTm29bGx\nsSWMCgAASsLwRL5BgwYpMTHRtrxv3z5FRkbaHf/kk0/qwoULCg4OVpcuXWx/AACAaxlu6c+YMUNP\nPPGERo0apeTkZH3yySf6v//7P7vjL1++rNmzZzs0JAAAKDnD0g8ODtaCBQsUGRmp6tWrKy4uTpUr\nV7Y7vlGjRrfc/Q8AAFzDsPT37NmjF198UXPmzNG3336rkSNHKjY21nY2/w3jx4+XxWJRRkaGwsPD\n1aZNG5UrV872Osf0AQBwLcPSj42N1bJly1S/fn0NGDBAu3fv1siRI/Xxxx/nG9e1a1fb13369HF8\nUgAAUCKGpb9hw4Z8W+w9evRQ8+bNC4wbMGCAY5MBAACHMix9Hx8fbdiwQUlJSZoyZYrOnTun1NRU\n1alTp9DxN3bz38zPz0+tW7fWwIED5eFheMEAAABwAsMGjomJ0Zdffqk9e/ZIuv5o3blz59odX7Nm\nTSUnJ6tt27Zq27atUlJSJF2/Te+t3gcAAJzLcEv/4MGDev/9922776tVq6arV6/aHZ+YmKhVq1bJ\nx8dHkjRkyBA98sgjWrlypfr372/3ffv379eLL76ovLw8DR48WGPHji0w5oMPPtCiRYtksVjUtGlT\nvfLKK4bfIAAAuM6w9MuVK5dvd31eXt4tx//666/y9vb+31/g5aW0tDT5+PjY/iHwR7m5uZo1a5aW\nL18uf39/DRo0SKGhoWrYsKFtTFJSkpYuXap169apSpUqunDhguE3BwAA/sew9Bs3bqxt27bJarXq\n3LlzWrp0qdq2bWt3fIcOHTR27Fj169dPFotF27ZtU7t27ZSZmWm39OPj4xUYGKi6detKksLCwrRn\nz558pf/ee+9p+PDhqlKliiSpRo0at/WNAgBgdoalP23aNM2bN0/nz5/Xgw8+qNDQUE2bNs3u+Jkz\nZ2r9+vX66KOPZLVade+992ro0KHy9vbWe++9V+h7UlJSFBAQYFv29/dXfHx8vjFJSUmSpKFDhyov\nL0/jxo1T586dDb/BatUqysvL03BcUdWs6eewuUqDM/MWd253+wwl98vsbnklMpcGd8srkdnRDEu/\nUqVKmjNnTpEn9Pb21ogRIzRixIgiv8dqtRZY98crAHJzc/Xjjz9q1apVSk5O1vDhw7Vjx45b3h1Q\nktLSsoqcoyjOn0936HzO5sy8xZm7Zk0/t/sM3S2zu+WVyFwa3C2vROaSZLDHsPSLauXKlRo5cqT+\n9a9/FShsSZo6dard9wYEBCg5Odm2nJKSolq1auUb4+/vr9atW8vb21t169ZVvXr1lJSUpDvvvNNR\n3wIAAH9qDrto/sYNfHx9fVWxYsUCf26lZcuWSkpK0tmzZ3Xt2jXFxcUpNDQ035ju3bvr4MGDkqSL\nFy8qKSnJdg4AAAAw5rAt/aFDh0qSxo0bd/shvLw0c+ZMjRkzRrm5uYqIiFCjRo0UGxuroKAgdevW\nTffdd58+//xzPfDAA/L09NTUqVNVrVo1R8UHAOBP75aln5ubq8cff/yWj9L9owsXLigmJka//PKL\n1qxZo8TERB05ckSRkZG3fF9ISIhCQkLyrZswYYLta4vFoujoaEVHRxc5CwAA+J9b7t739PTUb7/9\nZnht/s2effZZtW3bVpcvX5Yk1a9fX2vXri1ZSgAAUGKGu/dbtWqlcePGqU+fPvL19bWt/+NW+Q0p\nKSmKjIzUu+++K+n6vfu53z4AAK5nWPoJCQmSpHXr1tnWWSwWu6Xv5ZV/ysuXLxd6SR4AAChdhqW/\natWq25rw/vvv18yZM5WZmanNmzdr7dq1ioiIKHZAAADgGEU6e//TTz/VF198IYvFok6dOqlTp052\nx44ZM0bbtm3T5cuXtW/fPo0YMUL9+vVzWGAAAFA8hqX/1ltvaevWrQoLC5MkzZs3T/3799fo0aML\nHX/lyhWFh4crPDzcsUkBAECJGJb+tm3btH79elWqVEmSNGLECEVGRtot/S5duqhhw4bq2LGjOnbs\nqNatWxc4zg8AAEpfkU6rv1H4f/y6MJ9//rkmT54sDw8PxcbGqlOnTnr00UdLlhIAAJSY4SZ4UFCQ\noqOjNXjwYFksFm3YsEFBQUF2x3t6euqvf/2r7U9qaiqX7AEAUAYYlv6MGTP0xhtv2J60d8899+jx\nxx+3O75Pnz4qX768unTpokGDBmn27Nns3gcAoAwwbOOKFStqypQpRZ6wXbt2+vrrr3Xo0CF5eXnJ\n29tbQUFBbO0DAOBidkt/586d6t27t9asWVPo68OHDy90/fPPPy9JSk5O1ieffKKJEycqPT1dX331\nVcnTAgCAYrNb+t9//7169+6t48eP39aEx48f14EDB/TFF1/oxIkTat68uYKDg0scFAAAlIzd0h8/\nfrzy8vLUq1cvu7fcLcycOXPUsWNHPfbYY7rrrrvk4+PjkKAAAKBkbnlM38PDQ6+99tptlf769etL\nHAoAADie4dl1TZs2VXx8fGlkAQAATmR49v63336ryMhIBQYGqmLFirb1GzdudGowAADgWIal/+yz\nz5ZGDgAA4GSGpd+hQwdJ0sWLF1W9evUiTXrgwAGdOnVKDz30kH799Velp6erXr16JUsKAABKxPCY\n/tGjR9W1a1cNGDBAknTs2DHNmDHD7vilS5dq0aJFeueddyRJOTk5mj59uoPiAgCA4jIs/ZiYGL31\n1luqVq2aJKlly5b6+uuv7Y7fsWOHVqxYYTv+HxAQoIyMDAfFBQAAxWVY+r///rsaNmyYb523t7fd\n8eXLly/wusViKWY8AADgKIbH9H18fJSZmWkr7pMnT6pcuXJ2xwcEBOjw4cOyWCzKy8vTkiVL1KhR\nI8clBgAAxWJY+v/4xz80evRopaamatq0afr00081f/58u+NnzJihZ555Rt9//71atWqldu3a6eWX\nX3ZoaAAAcPsMSz8kJET169fXp59+KqvVqqioKAUGBhY6Ni8vTxcuXNDbb7+t7Oxs5eXlydfX1+Gh\nAQDA7TM8pv/WW2+pbt26GjZsmIYPH67AwEC99dZbhU/m4WF7DG+FChUofAAAyhDD0v/ggw+KtO6G\nwMBAnTt3rmSpAACAw9ndvf/555/rs88+U2pqql566SXbeqPL7zIzMxUeHq62bdvmu21vbGysA+IC\nAIDislv63t7e8vX1lcViyVfetWrV0tixY+1OGB4ervDwcMemBAAAJWa39Dt06KAOHTro/vvvV+PG\njYs84Y079wEAgLLF8Oz9zz//XLVr15afn5+mTJmiY8eO6dlnn9W9995b6Pjx48cXejMedu8DAOBa\nhifybd68WX5+fvryyy918eJFzZ07V6+++qrd8V27dlWXLl3UpUsXBQcH6/Lly6pZs6ZDQwMAgNtn\nuKXv6ekpSTp48KD69u2ru+66S1ar1e74P+7eHzhwoEaPHl3CmAAAoKQMt/TLly+vpUuXKi4uTp06\ndZLVatXvv/9e5L/AYrEoJSWlRCEBAEDJGW7px8TEaO3atZo8ebJq1qypM2fOqG/fvnbH33xM32q1\n6r///a+Cg4MdlxgAABSLYenXq1dP//znP5WVlaWsrCz97W9/02OPPWZ3fNeuXW1fe3p6atSoUWrd\nurVj0gIAgGIzLP0zZ85o8uTJSkhIkMViUfPmzTV//nzVrVu30PEeHh7q169fvnVbt24tsA4AAJQu\nw2P6zz33nB588EHFx8fr6NGjGjx4sGbOnGl3/IoVK4q0DgAAlC7DLf2LFy9q0KBBtuWIiAi98847\nBcYdO3ZM8fHxSktL05o1a2zrMzIybuvEPwAA4ByGpe/h4aEffvhB9evXlySdPn3adhnfzVJSUnT8\n+HFlZ2fr+PHjtvW+vr6KiYlxYGQAAFAchqU/adIkDR8+XM2aNZMkJSYm5nsAzw3du3dX9+7d9dln\nn9m9Wx8AAHAdw9Lv3Lmz4uLidPToUVmtVrVu3VrVq1e3O/7QoUPq1KmT7bK9rKwszZgxQ6+88orj\nUgMAgNtmeCJfvsEexsPPnz+vYcOGKSUlRQkJCYqIiFBgYGCxAwIAAMcw3NLftWuXZsyYoaCgIOXl\n5SkxMVGzZ89W9+7dCx0fExOj999/XwMGDJCPj4/mzZunjh07Ojw4AAC4PYalv2DBAq1fv1716tWT\nJCUlJSkqKspu6WdkZGjv3r1q1KiRfvnlFx0/fpzSBwCgDDDcX1+uXDlb4UvS3//+d5UvX97u+IiI\nCDVr1kwrVqzQxo0bdeTIEY0ZM8YxaQEAQLEZln63bt20ePFinT9/XqmpqVqyZIm6deumK1euKDs7\nu8D42bNnKyoqShaLRZUrV9Ybb7yh++67zynhAQBA0Rnu3n/jjTckSbGxsfnWL1q0SBaLRQkJCfnW\nd+jQocAcI0eOLElGAADgAIaln5iYWKSJpkyZovnz5ysiIsJ2ud7NNm7cePvpAACAwxiWflHd2Jp/\n5plnHDUlAABwIIeVflBQkKTCd+8DAADXc1jp3/DDDz9oyZIlOnPmjHJycmzr2b0PAIBrGZb+qVOn\n1KBBA8N1N0yYMEH9+vXTgAEDCn0wDwAAcA3DS/YmT55cpHU3eHl5acyYMQoODlaHDh1sf4pi//79\n6tmzp3r06KGlS5faHffhhx+qSZMmOnbsWJHmBQAAt9jSv3jxoi5evKirV6/q1KlTslqtkqT09HRl\nZWXZnfC+++7T/v371blz59sKkpubq1mzZmn58uXy9/fXoEGDFBoaqoYNG+Ybl5GRoVWrVqlVq1a3\nNT8AAGZnt/S3b9+ulStXKjU1VY8++qhtvZ+f3y3vsBccHKzHH39cHh4e8vHxkdVqlcVi0YEDB24Z\nJD4+XoGBgapbt64kKSwsTHv27ClQ+rGxsRozZozefvvtIn2DAADgOrulP3LkSI0cOVJLlizRP/7x\njyJPOHPmTMXExKhFixZFeirfDSkpKQoICLAt+/v7Kz4+Pt+YEydOKDk5WV27dqX0AQC4TYYn8vXs\n2VNXr15VuXLl9OmnnyohIUFDhgxRlSpVCh1fpUoV9erV67aD3Dh8cLObb/KTl5enmJgYxcTE3Na8\n1apVlJeX404orFnTz2FzlQZn5i3u3O72GUrul9nd8kpkLg3ullcis6MZlv7EiRO1ceNGnT17Vs89\n95w6deqkZ555RkuWLCl0fPfu3bVu3Tr17t1b5cqVs62vUKHCLf+egIAAJScn25ZTUlJUq1Yt23Jm\nZqa+++47Pfzww5Kk8+fPKyoqSosXL1bLli3tzpuWZv/8g+I4fz7dofM5mzPzFmfumjX93O4zdLfM\n7pZXInNpcLe8EplLksEew9L38PCQt7e39u3bp8jISD366KPq16+f3fGvvfaaJOmFF16QxWKxHdP/\n4z36/6hly5ZKSkrS2bNn5e/vr7i4OL3yyiu21/38/HTw4EHb8ogRIzR16tRbFj4AAPgfw9K/evWq\nfv31V3388ceaOHGipMJ3xd9Q1Hv1Fwji5aWZM2dqzJgxys3NVUREhBo1aqTY2FgFBQWpW7duxZoX\nAABcZ1j6I0eOVK9evRQcHKyWLVvq7Nmz8vNzzvGKkJAQhYSE5Fs3YcKEQseuWrXKKRkAAPizMiz9\nIUOGaMiQIbblO+64Q8uXL3dqKAAA4HiG19RlZ2frtdde09NPPy1JSkpK0v79+50eDAAAOJZh6T//\n/PPKycmxHasPCAjQokWLnB4MAAA4lmHpf/fdd5o8ebK8vb0lSb6+vsrLy3N6MAAA4FiGpX+j7G+4\nevXqLc/eBwAAZZPhiXzt2rXTkiVLdO3aNR08eFDLly9XaGhoaWQDAAAOZLilP2nSJFmtVvn6+mr+\n/Pm688479eSTT5ZGNgAA4ECGW/re3t6KiopSVFRUaeQBAABOYlj6L730UoF1fn5+at26tYKDg50S\nCgAAOJ7h7v0LFy7oo48+Um5urnJzc7Vr1y599913iomJ0eLFi0sjIwAAcADD0k9NTdXmzZsVHR2t\n6Ohobdq0SRcvXtTatWu1ffv20sgIAAAcwLD0U1JSVKVKFdtylSpVdP78eVWqVEk+Pj5ODQcAABzH\n8Jh+w4YNNWPGDA0cOFAWi0WbN29WgwYNdO3aNXl4GP6bAQAAlBGGrT137lxVqlRJs2fP1gsvvKCK\nFStq7ty58vDw0LJly0ojIwAAcIBbbunn5uYqLi5OzzzzTKGvV69e3SmhAACA491yS9/T01Pvvvtu\naWUBAABOZLh7/+6779aHH35YGlkAAIATGZ7It2XLFi1fvlzly5dXhQoVZLVaZbFYdODAgdLIBwAA\nHMSw9Ddt2lQaOQAAgJMZln6dOnVKIwcAAHAyw9L/5ZdfNH/+fCUmJurq1au29Xv27HFqMAAA4FiG\nJ/JNnz5dwcHBslqtevnll9W2bVsNGDCgNLIBAAAHMiz9tLQ0DR48WF5eXmrTpo3mzZunffv2lUY2\nAADgQIal7+3tLUmqWLGifv75Z+Xk5OjixYtODwYAABzL8Jh+u3bt9NtvvykyMlIDBw6Uj4+Pevbs\nWRrZAACAAxmW/o1b8Pbv318dOnRQRkaGateu7fRgAADAsW7rMXl33HGHGjdurL59+zorDwAAcJJi\nPRvXarU6OgcAAHCyYpW+xWJxdA4AAOBkdo/pnzx50u6bcnJynBIGAAA4j93SHzt2rN03lStXzilh\nAACA89gt/b1795ZmDgAA4GTFOqYPAADcD6UPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCASVD6AACY\nBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCASVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASl\nDwCASVD6AACYBKUPAIBJUPoAAJhEmSn9/fv3q2fPnurRo4eWLl1a4PXly5frgQceUN++fTVy5Ej9\n9NNPLkgJAID7KhOln5ubq1mzZmnZsmWKi4vTjh07dPLkyXxjmjVrpk2bNmn79u3q2bOn5s+f76K0\nAAC4pzJR+vHx8QoMDFTdunXl4+OjsLAw7dmzJ9+Yjh07qkKFCpKk1q1bKzk52RVRAQBwW2Wi9FNS\nUhQQEGBb9vf3V0pKit3xGzduVOfOnUsjGgAAfxperg4gSVartcA6i8VS6NitW7fq+PHjWr16dZHm\nrlatory8PEuU72Y1a/o5bK7S4My8xZ3b3T5Dyf0yu1teicylwd3ySmR2tDJR+gEBAfl216ekpKhW\nrVoFxn3xxRdasmSJVq9eLR8fnyLNnZaW5bCcknT+fLpD53M2Z+Ytztw1a/q53WfobpndLa9E5tLg\nbnklMpckgz1lYvd+y5YtlZSUpLNnz+ratWuKi4tTaGhovjEnTpzQzJkztXjxYtWoUcNFSQEAcF9l\nYkvfy8tLM2fO1JgxY5Sbm6uIiAg1atRIsbGxCgoKUrdu3fTSSy8pKytLEyZMkCTVrl1bS5YscXFy\nAADcR5ns7wkNAAAO6UlEQVQofUkKCQlRSEhIvnU3Cl6SVqxYUcqJAAD4cykTu/cBAIDzUfoAAJgE\npQ8AgElQ+gAAmASlDwCASVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCASVD6AACYBKUP\nAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCASVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCA\nSVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCASVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ\n+gAAmASlDwCASVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCASVD6AACYBKUPAIBJUPoA\nAJgEpQ8AgElQ+gAAmASlDwCASVD6AACYBKUPAIBJUPoAAJgEpQ8AgElQ+gAAmASlDwCASZSp0t+/\nf7969uypHj16aOnSpQVev3btmiZOnKgePXpo8ODBOnfunAtSAgDgnspM6efm5mrWrFlatmyZ4uLi\ntGPHDp08eTLfmA0bNqhy5cravXu3HnnkEb388ssuSgsAgPspM6UfHx+vwMBA1a1bVz4+PgoLC9Oe\nPXvyjdm7d68GDBggSerZs6cOHDggq9XqirgAALidMlP6KSkpCggIsC37+/srJSWlwJjatWtLkry8\nvOTn56e0tLRSzQkAgLvycnWAGwrbYrdYLLc95o9q1vSTJG1/pV8J0rmGszK7+rO48d/EnbhbZnfL\nK5G5NLhbXonMjlZmtvQDAgKUnJxsW05JSVGtWrUKjPnll18kSTk5OUpPT1fVqlVLNScAAO6qzJR+\ny5YtlZSUpLNnz+ratWuKi4tTaGhovjGhoaHasmWLJOmjjz5Sx44dDbf0AQDAdRZrGToTbt++fZo7\nd65yc3MVERGhqKgoxcbGKigoSN26ddPVq1c1ZcoUJSQkqEqVKlqwYIHq1q3r6tgAALiFMlX6AADA\necrM7n0AAOBclD4AACZRZi7Zc6Xo6Gh98sknqlGjhnbs2OHqOIauXr2q4cOH69q1a8rNzVXPnj01\nfvx4V8cyFBoaKl9fX3l4eMjT01ObN292dSS7fvjhB02aNMm2fPbsWY0fP16PPPKI60IVwcqVK7Vh\nwwZZrVYNHjy4TOYt7Odt586dWrRokU6dOqUNGzaoZcuWLk75P4Xlfe2117Rnzx55eHioRo0aiomJ\nkb+/v4uT/k9hmRcuXKj33ntP1atXlyQ99dRTCgkJcWXMfArLPHHiRJ0+fVqSlJ6eLj8/P23dutWV\nMW0Ky5uYmKjnnntOWVlZqlOnjl5++WVVqlTJxUn/wArroUOHrMePH7eGhYW5OkqR5OXlWTMyMqxW\nq9V67do166BBg6xHjhxxcSpjXbt2tV64cMHVMW5bTk6O9Z577rGeO3fO1VFu6b///a81LCzMmpWV\nZf3999+tI0eOtJ4+fdrVsQoo7Oft5MmT1lOnTlkfeugha3x8vAvTFVRY3vT0dNvXK1eutM6YMcMV\n0ewqLPPrr79uXbZsmQtT3ZrR7+GYmBjrwoULSzmVfYXlHThwoPXgwYNWq9Vq3bBhg3XBggWuimcX\nu/cltW/fXlWqVHF1jCKzWCzy9fWVdP1+BTk5OVy66EQHDhxQ3bp1VadOHVdHuaVTp06pVatWqlCh\ngry8vNS+fXvt3r3b1bEKKOznrUGDBqpfv76LEt1aYXlv3nrLzs4ucz9/7vY7Tbp1ZqvVqp07d6pP\nnz6lnMq+wvKePn1a7du3lyR16tRJu3btckW0W6L03VRubq769eune+65R/fcc49atWrl6khFMnr0\naA0cOFDvvvuuq6MUWVxcXJn6ZWNP48aNdfjwYaWlpSk7O1v79+/Pd8MrONaCBQsUEhKi7du3a8KE\nCa6OUyRr1qxR3759FR0drUuXLrk6TpEdPnxYNWrU0N///ndXR7mlxo0b254Z8+GHH9puJleWUPpu\nytPTU1u3btW+ffsUHx+v7777ztWRDK1bt05btmzRW2+9pTVr1uirr75ydSRD165d0969e9WrVy9X\nRzHUoEEDjRkzRqNGjdKYMWPUpEkTeXp6ujrWn9akSZO0b98+9e3bV6tXr3Z1HEORkZHavXu3tm7d\nqlq1amnevHmujlRkO3bscIt/eL/44otau3atBg4cqMzMTPn4+Lg6UgGUvpurXLmy7r77bn366aeu\njmLoxolONWrUUI8ePRQfH+/iRMb279+vFi1a6C9/+YuroxTJ4MGDtWXLFq1Zs0ZVq1ZVYGCgqyP9\n6fXp06dM7sb9o7/85S/y9PSUh4eHBg8erGPHjrk6UpHk5ORo9+7deuCBB1wdxVCDBg309ttva/Pm\nzQoLCyuTN4+j9N3QxYsXdfnyZUnSlStX9MUXX5TZ46E3ZGVlKSMjw/b1559/rkaNGrk4lbG4uDiF\nhYW5OkaRXbhwQZL0888/a9euXW6xdeSOkpKSbF/v3bu3zP/8SVJqaqrt6//85z9u8fMnyfb77ean\nsJZVN37+8vLytHjxYg0dOtTFiQrijny6funKoUOHlJaWpho1aujJJ5/U4MGDXR3LrsTERE2bNk25\nubmyWq3q1auXxo0b5+pYt3T27Fk98cQTkq6fj9CnTx9FRUW5ONWtZWdnq0uXLvrPf/4jP7+y+9Ss\nmw0bNky//fabvLy8FB0dreDgYFdHKqCwn7eqVatq9uzZunjxoipXrqxmzZrp3//+t6ujSio87/79\n+3X69GlZLBbVqVNHL7zwQpm6ZK+wzIcOHVJiYqIkqU6dOpo1a1aBh5q5kr3fw9OmTVOrVq0UGRnp\n6oj5FJY3KytLa9eulST16NFDTz/9dJk7yZPSBwDAJNi9DwCASVD6AACYBKUPAIBJUPoAAJgEpQ8A\ngElQ+kAZEBoaqnvvvVe5ubm2dZs2bVKTJk1sd3tbt26dVqxYYTjXihUrbNcLo+QSEhL0wQcfuDoG\n4BCUPlBG1KxZU5999plt+f3331eLFi1sy5GRkUV6VO4777xTrNLPycm57feYQUJCgj788ENXxwAc\nwsvVAQBcN2DAAG3evFkhISE6e/assrOz1bhxY9vrCxcuVFZWlp555hm98cYbSkhI0KJFi5Sdna3B\ngwdrypQpOnHihFJTUzV+/HiVK1dOr7zyipYtW6agoCA99NBDkqRp06bZlqdNmyZPT0+dPn1amZmZ\n2rp1q44ePaqXX35ZmZmZkqTx48erS5cuBfKmp6dr7ty5On78uCwWi9q1a6eZM2cqMzNTc+bMsd3m\nNTw8XGPHjpUkjRgxQi1atFB8fLx++uknPfzww/L399fq1auVmpqqKVOmqHfv3pKkJk2a6IknntCe\nPXt05coVPfXUU+rZs6ek67dHfvXVV5Wbm6vq1atr1qxZCgwM1MGDBzV37ly1atVKR44ckcVi0YIF\nC9SgQQNJ0pYtW7R27Vrl5uaqUqVKev7551W/fn1t3rxZO3bsUOXKlfX999/Lz89PCxculJeXl15/\n/XVlZGSoX79+at++vZ599lnn/A8AlAJKHygj7r77bq1du1aXLl3Sli1b1L9/fx0/frzQsVFRURoz\nZoxWrVqlEydOqHPnzgoJCVFISIg2bNig119/Pd8/GG4lISFBq1evVsWKFXX58mU999xzWrp0qWrV\nqqXU1FQNGjTIVog3mzt3ripWrKitW7fKw8NDFy9elCS9+eabysvL0/bt25WZmakhQ4aoSZMmCgkJ\nkSQlJydr9erVOn/+vO6//3498sgjWr9+veLj4zVu3Dhb6UuSh4eHtm7dqh9++EGRkZFq166dJGnq\n1KlavXq1GjZsqA0bNmjy5MnasGGDJOnkyZOKiYnRrFmztHjxYr355pt65ZVXdPjwYe3cuVNr1qyR\nj4+P9u3bp+nTp2v9+vWSpGPHjmnbtm2qXbu2nn32Wa1evVqTJk3S+PHj9cknn+j111+/jf+aQNlE\n6QNlhMViUe/evRUXF6cPPvhA69ats1v6Hh4emj9/vvr166c77rjDduvP4ujVq5cqVqwoSTpy5IjO\nnTunRx99NF+uH3/8US1btsz3vo8//libN2+Wh8f1o4TVq1eXJB04cEDTp0+XxWJRpUqVFBYWpgMH\nDthKv1evXvLw8JC/v7+qVq2q7t27S5JatGihlJQUXb16VeXKlZMk2+2w69evr+bNm+ubb76RxWJR\n06ZN1bBhQ0lSRESEXnjhBduzHerVq6fmzZtLklq3bq2PP/5Y0vV75CcmJtrmtFqttmdYSNJdd92l\n2rVrS5JatWqlL774otifKVBWUfpAGTJw4EANHjxYHTp0ULVq1W459ty5c/Lw8NClS5d05coVVapU\nqdBxnp6eysvLsy1fvXo13+s3Cl+6XoRNmjTRmjVriv09WK3WAvcbv3n5RqHfyHZj+cZjgHNycvKN\n+eO8hc1/s5sfZ+rh4WE7V8FqtSoiIkITJkwo9H1/zHXzSZXAnwUn8gFlSN26dTVp0iQ9/vjjtxx3\n6dIlTZ48Wa+++qrCwsI0Y8YM22u+vr5KT0+3Lf/tb3+zHV9PTU3VwYMH7c7bpk0b/fjjj/ryyy9t\n6+Lj41XYIzq6du2qf//737bXbuzev+eee7Rx40ZZrVZlZGTogw8+KPaDfzZt2iTp+lPtEhIS1KpV\nK7Vp00YJCQk6deqUpOvH6Zs3b273Hz03hIaGauvWrUpOTpZ0/cFP9vak3KxSpUr5Pk/AnVH6QBkz\nZMgQNW3a9JZjpk+froiICLVr105PPPGEfv31V61bt06S9PDDD2v69Onq16+fTp48qQcffFDJycl6\n4IEH9Pzzz+vOO++0O2+VKlX05ptv6o033lB4eLh69+6tRYsWFVr60dHRyszMVJ8+fRQeHq4333xT\nkvT444/LarWqb9++Gjp0qMLDw9W5c+difRa5ubnq37+/HnvsMc2aNUs1atRQ9erV9dJLL2ny5Mnq\n27evtm3bpvnz5xvO1b59e02cOFFRUVEKDw9Xnz59tGfPHsP3BQcHKzs7W+Hh4ZozZ06xvg+grOAp\newDKpCZNmujrr7+Wr6+vq6MAfxps6QMAYBJs6QMAYBJs6QMAYBKUPgAAJkHpAwBgEpQ+AAAmQekD\nAGASlD4AACbx/wHh3zqCcEurbgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(8, 6))\n", "\n", "ax.bar(np.arange(K) + 1,\n", " trace['w'].mean(axis=0).max(axis=0));\n", "\n", "ax.set_xlim(1 - 0.5, K + 0.5);\n", "ax.set_xticks(np.arange(0, K, 2) + 1);\n", "ax.set_xlabel('Mixture component');\n", "\n", "ax.set_ylabel('Largest posterior expected\\nmixture weight');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since only three mixture components have appreciable posterior expected weight for any data point, we can be fairly certain that truncation did not unduly influence our results. (If most components had appreciable posterior expected weight, truncation may have influenced the results, and we would have increased the number of components and sampled again.)\n", "\n", "Visually, it is reasonable that the LIDAR data has three linear components, so these posterior expected weights seem to have identified the structure of the data well. We now sample from the posterior predictive distribution to get a better understand the model's performance." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 5000/5000 [01:21<00:00, 61.64it/s]\n" ] } ], "source": [ "PP_SAMPLES = 5000\n", "\n", "lidar_pp_x = np.linspace(std_range.min() - 0.05, std_range.max() + 0.05, 100)\n", "x_lidar.set_value(lidar_pp_x[:, np.newaxis])\n", "\n", "with model:\n", " pp_trace = pm.sample_ppc(trace, PP_SAMPLES, random_seed=SEED)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below we plot the posterior expected value and the 95% posterior credible interval." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n", " (prop.get_family(), self.defaultFamily[fontext]))\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFeCAYAAACy1qeuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4XMXVuN9btqrbKu69V9wbYFxCJzFfiEMgDoQS+BKH\nXyAQTDck2EmAUAzhA0ICCSUQMCY2EMA2JhjjItu49ya5qFl1tf3e+/tjtWtJe3e1KxfZ0rzPo8fW\n3pm5M/eu5sw5c+YcyTAMA4FAIBAIBGccuaU7IBAIBAJBW0UIYYFAIBAIWgghhAUCgUAgaCGEEBYI\nBAKBoIUQQlggEAgEghZCCGGBQCAQCFoIIYQFAoFAIGghhBAWCM4AU6dOZdWqVVGfr1mzhgsvvDDy\n+6xZsxg6dCgjRoxg5MiR/M///A8vv/wyfr8/qu7ChQvp378/H3/8cVSbAwYMYMSIEYwYMYJLLrmE\n999/P2bfDh8+TP/+/SPlJ06cyG233cbXX3+d8PgWLlzIj370o4TLCwSCEEIICwRnGQ8//DAbN25k\n5cqV3HvvvXz00UfceuutNI6r88EHH5CZmcmiRYui2sjNzWXjxo1s2LCB+++/n4ceeoj9+/fHve+6\ndevYuHEjH374IRMnTmT27NksXLjwlI5NIBA0RAhhgeAsxel0Mm7cOF588UW+/fZbVqxYEbl25MgR\n1q1bx2OPPcbKlSspKyszbUOSJCZPnkxGRga7du1K6L45OTnccMMNzJ49myeffBJd1wF4+eWXmT59\nOiNGjODyyy/n888/B2Dfvn088sgjfPvtt4wYMYLRo0cDsGLFCmbMmMHIkSOZPHkyCxYsOImnIRC0\nToQQFgjOcjp16sSQIUPIz8+PfLZo0SKGDBnCJZdcQu/evVm8eLFpXV3XWbZsGRUVFXTv3j2p+158\n8cUcP36cAwcOANC1a1fefPNN1q9fz+zZs7nnnnsoKSmhd+/ePProo5x33nls3Lgx0k+Hw8Ef/vAH\n8vPzeemll3j77bdZunRpM5+CQNA6EUJYIDgHyM3NpaqqKvL7hx9+yJVXXgnAlVdeyQcffNCgfElJ\nCaNHj2bYsGHMnj2bOXPmMGjQoKTvCVBZWQnAZZddRl5eHrIsc/nll9O9e3c2b94cs/64cePo378/\nsiwzYMAArrjiCtauXZtUHwSC1o4QwgLBOUBxcTEZGRkArF+/nsOHD3PFFVcAISG8e/duduzYESmf\nm5tLfn4+GzZsYNasWaxevbpZ9wTIzMwEQtr39773PUaPHs3o0aPZs2cPFRUVMetv2rSJWbNmMX78\neEaNGsU///nPuOUFgraIEMICwVnOsWPH2LZtW2SvddGiRRiGwYwZM5g0aRIzZ86MfN4Yq9XK3Xff\nze7du5M2BX/++ee0b9+enj17cuTIER588EEeeugh1qxZQ35+Pn379o2UlSQpqv6vf/1rpk2bxpdf\nfsn69eu59tpro5zLBIK2jhDCAsEZIhAI4PP5Ij/BYDBueY/Hw9q1a/n5z3/OsGHDmDx5Mj6fj08+\n+YTHHnuMRYsWRX4eeughFi9ebNqm1Wrlpptu4oUXXkion2VlZbzxxhs8//zz3HXXXciyjMfjQZIk\n2rVrB8D777/Pnj17InXat29PcXFxg6NUtbW1ZGRkYLPZ2Lx5M0uWLEno/gJBW0Jt6Q4IBG2Fn/3s\nZw1+v/3225k4cWJUuccee4x58+YB0L17dy655BJuuukmZFlm6dKl2O12ZsyYgcViidS55ppreO65\n5/jqq69wOp1RbX7/+99nwYIFLF++nKlTp5r2b8yYMRiGgcPhYMiQITz77LORM8x9+vThpptu4tpr\nr0WSpIjXc5jx48fTp08fzj//fCRJYs2aNTzyyCP84Q9/4LHHHmPs2LFcdtllVFdXJ//gBIJWjGQI\n+5BAIBAIBC2CMEcLBAKBQNBCCCEsEAgEAkELIYSwQCAQCAQthBDCAoFAIBC0EEIICwQCgUDQQpzx\nI0qlpTVn+panjawsJxUV7pbuxmmlLYwR2sY4xRhbD21hnK1tjDk5aaafC034JFBVpaW7cNppC2OE\ntjFOMcbWQ1sYZ1sYIwghLBAIBAJBiyGEsEAgEAgELYQQwgKBQCAQtBBCCAsEAoFA0EIIISwQCAQC\nQQshhLBAIBAIBC2EEMICgUAgELQQQggLBIJWy4UXjuXGG69j1qyZPPjgvXi93qTbePfdt5pV7y9/\n+T/WrVuTdL2W5uOPF1NWVppUnWPHjjJr1syTvvepaudcQghhgUDQarHZbLz22lv84x/vYrFYWLTo\nvaTbePfdt5MWwpqmccsttzNmzLik6pwNNEcIC5rPGQ9bKRAIBC3B8OHnsXfvXgD++c83+OijfwNw\n1VUzmDnzOjweDw8/PIeSkhJ0XePGG2+hvLycsrJS7rjjNjIyMlmw4CXWrl3Nq6++RCDgp1OnLtx/\n/yM4nU6uueYqpk79Dvn5a7juup+wZs03TJx4PlOmTCc/fy0vvPAMmqYxYMAg7r77PqxWa1Sd6dMv\nifS3oqKCJ5+cR3FxMQB33HEXw4adxzPPPEHHjrn88Ic3sGbNN/z9739lwYKXmD//MaxWKzt37qC2\ntpZf/vJOJk26AE3T+L//e56NG9cTCPi5+uofMGPG9wF4883X+fTTj5EkmfHjJzJgwEB27drBo48+\niM1m56WX/sqBAwd4/vmncbvdZGZmcv/9c8nOzmbnzh3Mn/8YAGPHjjd95g8/fB+XXXYFEyacD8Dj\nj89l0qQL6N9/IL/97cN4vR4A7rzzNwwdOrxB3YULF7Ju3QbuuuteAH7zm19x7bU/ZuTI0THfwbmI\nEMICgeC0M3fugyxevOiUtnnVVTOYO/d3CZUNBoOsXr2KceMmsnPnDj7+eDEvv/w6hmHws5/dyHnn\njeTo0SNkZ+fwxBPPAuByuUhNTeWdd97kuedeIjMzk8rKSl5//VWeeebPOBwO3njjNd55501++tNb\nAcjIyOCvf30TgDVrvgHA5/Mxb96jPPPMn+nWrTu//e3DLFr0HjNnXhdVpz7PPvskM2dez/Dh51FU\nVMSvfz2bN998j9tv/yW3334jffsO5plnnuDJJ59DlkNGzWPHjvHKK69z5Mhh7rjjdkaPHst//vMR\nKSkp/OUvf8fv9/O//3szY8eO59Chg3z11Ze8/PLr2O12qqurSE/P4P3332X27F8xYMAggsEgzzzz\nBPPnP0VWVhbLln3Gyy+/wP33P8L8+Y/yq1/dw4gRo3jhhWdNn/u0aRezbNnnTJhwPoFAgPXr13H3\n3XMwDHj66Rew2WwUFhYwd+4DvPrqPxJ6l029g3MNIYQFgjj4AhpVLh8ZqTZslrMzlm1tbS0OhyMy\nEQtO4PP5uPHGkLAbPvw8rrzye3zwwXtceOEUHA4HAJMnT2HTpm8ZN24CL7zwLH/+83NMmnQBw4eP\niGpv27YtHDy4n//935sBCAYDDB48NHJ92rSLo+oUFByiY8dOdOvWHYDLLruShQv/FRHCZnUA8vPX\ncvDggcjvtbW1uN21OJ0p/Pa3v+XHP/4xv/zlnXTu3CVSZurU6ciyTNeu3ejUqTMFBQdZt241e/fu\nZcWK5XXtuDh8uJD8/LVcfvlV2O12ANLTM0z6fpD9+/dx552/AEDXNdq3z8blclFTU8OIEaMAuOSS\ny1m9+uuo+uPHT+TZZ5/E7/ezZs0qhg8fgc1mx+Vy8fTTf2DPnt3IskJh4SHTZ2BGU+/gXEMIYYHA\nBE3XeWf5XjbuLqW82ke7dBsj+uXww6l9UM4iYbd16xa2bt2M1Wqjc+fO9OvXn6ysdi3drSjmzv1d\nwlrrqSS8J9wQw7Rst27defXVf/DNN1/zyisvMmrUmCjtyjAMRo8ex6OPzjNtw253mHxqfr/4dcAw\ndF566a/YbPaoa7t37yY9PYOysrIGn0uS1KikhGEY3HnnPYwbN6HBlTVrVpmUb9wH6NmzFy+99LcG\nn9fU1DRZF0LPf8SIkaxd+w3Lln0eMbe/886bZGW157XX3kbXdaZNmxRVV1EUdP3Es/P5/HV9iv8O\nzjXOntlEIEgQX0CjpMKNL3D6HFneWb6XpfmHOV7twwCOV/tYmn+Yd5bvPW33TAafz8fy5cvYvn0r\nVqsVMDhy5DCfffYfPv30Ew4ePNjSXTxrGT58JF99tQKv14vH4+G///2C4cPPo6ysFJvNziWXXM6P\nfjSL3bt3AuB0OnG7awEYPHgoW7Zs4vDhQgC8Xi8FBfG1uG7denDs2NFInU8//ZjzzhvZZD/HjBnP\n+++/G/l9z55dABQVHeNvf/sbf/vbm6xevYpt27ZGynzxxVJ0XefIkcMcPXqEbt26M3bsBBYteo9g\nMAiENHOPx8OYMeP56KN/R5zOqqur6sabgtvtrut7dyorK9i6dTMQMuvv37+PtLQ0UlNT2bTpWwA+\n++yTmOOYNu0SPvpoMZs3fxtZCNTWumjfPhtZlvn0049NndI6d+7M3r270XWd4uIiduzYBjTvHZzN\nCE1YcM5wprRTX0Bj425z79CNu8v4/uTeLWqaLi4u5ptvvkbXNVS14Z+wxWKhttbF2rXfsGvXdoYP\nH0GHDh1bqKdnJ/37D+Cyy67k1lt/AoT2lvv1G8CaNd/w5z8/iyTJqKrK3XfPAeC7372au+++g/bt\ns1mw4CUeeGAuc+c+QCAQ0sxuvfV/I6ZmM2w2G/ff/wgPPXRvxDEr7BgVj1/96h7+9Kc/cMMN16Jp\nGsOHj+Duu+9j/vzf8pvf/Ibs7BzmzHmIefPm8sorfwcgL68Dt956A7W1tdx9933YbDauumoGRUXH\nuOmm6zEMg8zMLObPf4rx4yeyZ89ubrllFqpqYcKESdx22y+4/PIreeKJeRHHrN/97g8888yTuFwu\nNE1j5swf0atXb+677xHmz38MSZIYOza2F/jYseP53e8e4fzzL8RisQBw9dU/4MEHf8N//vMR48ZN\niGwN1GfUqFF07NiJH//4B3Tv3pN+/foDkJWVlfQ7OJuRDMOIbys5xZSW1pzJ251WcnLSWtV4zDib\nxvjW0t0szT8c9fn00V24bnq/k2q7/jhLKtzc99JqUyOiLMG8n40nN6tlPDEPHjzAunVrUZTEFh3B\nYIDs7FzGj59I9+55Z827PF2cTd/X04nZOB9/fG7EG7s10NreZU5OmunnwhwtOCdoSjtNxDSdqBk7\nI9VGu3Sb6bWsNDsZqebXTjc7d+5g7do1CQtgAFW1UFlZwaeffkJJSclp7J1AIGgOwhwtOCeocvko\nr/aZXquo8VLl8sXUTpM1Y9ssCiP65Zhq3SP6ZbeIKXrTpo3s3r0LVW3evQ1D57PPPqNfv6H07Nnr\nFPdOcDbwwANzW7oLgmZwzgvhYDDIzp07CAaD+P0+gsEgkiRhtztwOBxkZGSQlpaGJMlIkhTx6NM0\nDU0L4vcHkGWZzMxMFOXsPIIiOKGdHjcRxE1pp2EnqzBhJysgphn7h1P7ACEtu6LGS1aanRH9siOf\nnyl0XWf16lUcPlwYtf+bLIqisG7dWqqrqxk+/LxT1EOBQHAynPNC2O2uZfPmjTHd/AOBAJoWBE64\n00uShGHogFR3ttJAkhTS0tJIT08nJSUVXdcIBAJ1P0EsFhWLxYrNZsVisaKqKuXlGbhcAWw2O4ah\nEwwGCQQCBIMBUlJSadeufeQMHoQ8Wo8dO0p5+XGCwSDp6Rl07NiR9PSMBu7+uq6jaVrEieFUcy6c\nfW1Mc7XT5jpZKbLMddP78f3JvaOe1Zl6fjU1NXz11Ze43bUxBbCmgzcoY1d1ErFSq6rCrl07CAYD\njBo15hT3WCAQJMs5L4SbwmKxJCzMvF5PXRi14oTKFxbaqaysRdd1AGRZQpJkZFkmGAyiaUEsFgtO\nZyrBYJDa2lpUVYlo3IZRwLffrkdRFGw2O8GghqYF0DSD0MJAivQ/tACwY7fbsNvtWCxWNE0jEPDX\nLTRCnrJWqxWr1YbVasHvD+B21+LxePB4PAQ1g61FKkeqVGr9kGqT6N5eZmJfK1arBcMwIguPcP8z\nMpxUVbkxDAMkK4othY45GXTr2oWUlFT8fj8+n5eamhp8Ph9Op5OUlBSsVhsWS6hNXdfRdR3DMLDZ\nbAmdLzSjOdrpyZixIST8w9dPpXd2U4L80KGD5OevRZIk0yAcugHbS50Uu6x4gjIOVScv1c+gHDdy\nE49XVVX27duLqqqmASkEAsGZo9UL4dOJJEkoimJqxlZVNaK9hOOj2mzWqPpWa8iMGj7DpygqjZsL\nC8bwWcXmsrXEycHKEwsSlw+2HdWpdVcwJNdtWkeSglRVe+omfB1PsBaHWkOOYzf92rlCp/nRUVUL\nsizXmfk1QAJJxhOQsSkasmQgSQayLEe2CsLWi2BQwzBC9UILkWDkX1lW6qwPFiwWG11sFrL6y6iW\nVLIz7CiSh/x1a3C7PXg8bnRdJz09nbS0dDIzM0lNyyQzVaXCFYwaW1aajfQUa9TnAC5XDUePHqGq\nqgqbzUZ6egb/3elhxbdFkTIRs7YB132noVk7GAzi9XqpqanGYrHSrl270PMxEeTDerdn+uiuZKZa\nqa4sZ//+vRQWFsQ1P28vdXKw8oT1xxNUIr/Hepf1UVWVXbt2oigqQ4acu9GGBIJzHSGEzwDJmgxP\nR9uaDsUuc4FT7LIyMNsds77ZhF9Qk4asqFETvizLKKolrpamaRoulwuXyxV3bOGtAp/Ph8frY3up\nTrFLq2vTTV5qqanmV1ZWSllZKbqu4/f7SZMzqSAlqv0Uo4x/L3oXVbWSnZ1ObW0AVVXqtHovVusJ\nrT0Q1Fl9KAuzP5mV3x5EOp6PRZWRZRVdDxIIaBiGXrdAC1k10tLS2VGWytYjeqTu8WofX2w8yhcb\nj2BXNHKdXoZ08MUVwCfzLuujqio7dmxDVVUGDBjYdIWTwDAMamvjv+9kSUlJbdKq8u67b7N48QcY\nBnz3uzMioSJfffUlFi9eRGZmFgC33fZzJkw4n82bv+Wpp36PxWJl7tzH6dKlKzU1NTzyyH089dSC\nZltxmmLPnl2UlZVGEh0kSllZKc888wS/+90fT0u/kmX27J9F4k7fffcdPPLI47hcNfzmN7/iH/94\nN275RPjLX/6P4cNHxM1OtWFDPhaLJSohxKlmw4Z8/vnPN/jjH585qXaEED6NnIzJ8FS37Q3KeILm\nM7MnKOMNyqRY9ahrQS35Cf9ktTQzmtNmSOu2MyTPiyzLJs/KjyxZAAO/34/H44nUbRwq0K+reDXz\n/V+vphDEhk3WAaNOe48u66r1sK/YBpi1I+HVVApqUpEVJe5zau67NENRFDZv3oTFotK7d9+E6jSH\n2loX//rXP7HZTs3xLp/Pxw9+cC2pqeZnLwH279/L4sUf8Morf0dVVX796zuYMOF8unbtBsDMmddx\n3XWzGtT55z/f5IknnuXYsWN88MF7/PKXd/Laa39h1qyfnjYBDLBnz2527tyelBAOBoNkZ+ckJYCD\nwWDSDn7NqQPw5JPPASGr0qnilltub7LMxo3rcTicSQnh5o7xVHDG75qV5Wz2MQszrFadtDSHacSV\nM0F6euz75hdYGph/w4LDalUZ3S1wUvdNtm2nBs4jBu5A9ESSYjXIybJh9lpqvFLcCV+1O0mznwhr\nEdSg1G0+0Za6bThTDNP7xONUtDkxwyCo+fAEJByWcHlHpP0ar4QzxRGzHacGKUcNav3mE/Hh2lRG\ntQ/EXVzFe5bJjCleX+K9S4j9fd27dzt5eVn07Nmzyf41B7sdsrMzT9nfqcfjITs7jbS0aCEcDoqQ\nn1/EqFEj6do1B4BJk8azfv0qRo4cTEqKDafTFhVAISXFjtOpYLdLpKc78XgqqKmp4OKLL4rZl6lT\np3LppZfy1VdfYbPZeOqpp+jevTtHjhzh/vvvp7y8nHbt2jF//nw6derEJ598wgsvvIAsy6SlpdWF\noHwZr9fL9u1buO2227jooov47W9/y+7du9E0jdmzZzN9+nQWLlzIZ599htsd2nr5/e9/z+23386S\nJUvw+XzMnTuXrVu3oigKc+bMYfz48VF13njjjQb9X7RoEa+++iqSJNG/f3+eeOIJ5syZg9VqZceO\nHYwcOZL/9//+n2l/vF4v9913Hzt37qRXr17oepDMTCc5OWlMnTqV9957j3btUgCD3/9+Ltu3b6dv\n37784Q9/wOEIzVfh8itXrmTBggX4/X66du3K/Pnzo97PnDlzuOiii7j00kuZOnUqM2bM4IsvvqjL\n9vQMNpuNxYs/QJZlli//lIceeohevXrxyCOPcPToUQDuv/9+Ro0axYIFCygoKKCwsJBOnTpRWFjI\nvHnz6Ns3tBidNWsW9957L7quM2/ePLxeL3a7nXnz5tGrVy8yM51YrWrMIByJcsaFcEVF8zShWFRX\n11BT4yFwcjKtWaSnO6iu9kR9rulQ65cpLDcXHIXlMr3TPc02TQc02FcaHdQdoOB4dNthrdkXvS0K\nQI7Th7s2ehwAzhQHDlXHE4ye1R2qTtDrptp/4rNav0yt33yirfVLlFb4EtbSTleb7rr+JmtNyHFK\nMfohsafUgs8fZFjeie93460CTQeHajV9lsmOKVZf4r3LWN/XMJ98soxJk86nU6fOcfvXHFyuU/t3\n6vV6KCuroS7scYT6UZayszuzevVa9u4txGazs2zZFwwYMJDS0hpqa3289977vP/+Qvr3H8js2XeS\nnp7OzJk/5q677sZms/HQQ4/x+9//kVtu+d+4kZs0TUeWrfz1r2/xySdLmDv3Mf74x2d48MFHmDbt\nUi677EqWLPmQhx+ey/z5T/Hccwt46qkF5OTkUlNTQ1WVj5/+9Gfs3Lk9kjv3T396jsGDz+Ouu+6n\npqaGW2+9gX79hlFT42Xr1m0sWbKYQEDh2LGjBIMapaU1vP32G3i9Af7617c4dOggd975C95+e2Gk\nzuuvv016ekaDsezfv4/nn3+BF1/8K5mZmVRXV1FaWoPXG6C4uJQFC15BUZSY/fnww/cBhb///V32\n7t3DzTf/mMpKN6WlNWiazvHjtXg8bg4cOMA99zzAnDlzmTfvUV5++W9cd90s/P4glZVu9uwp5Lnn\nnufJJxdE0hP+7W9/44c/vKHRew9QXe2JtG+xOHj55b+zcOG/+POfX2LOnIe46qqrcTicESvHQw89\nwNVXz4xKCVlb62Pnzt28+OJfsNnsvPPOm7z//ofcfPNtlJWVcfRoEXl53amtdfH00y+iqirr1q3h\n97//I48//gSVlW78/mDCUb1iCes2a46OtZcab4+1qf3XxpN6LBqbDJtqt/H1rSUpaIZ5+15NZkux\nk2EdTgiSxqbcMIqk0zXDx6Cc2AsjVYG8VL9p/bxUf1R/7aoeV2jb1eQE8OlqE5I3cQ/KcaMbUFBl\np/6RtzAFVXYkYGCOm51l5sI91rNMdkzhd2Z2j+aiKDKrVq3kggsuIi8vr9ntnC306NGTH//4J9x5\n5y9wOJz06dM34kR59dXXcOONtyBJEq+88iLPP/8099//CH379ufll18D4NtvN5CdnYNhGDz88H2o\nqsrs2b+iXbv2UfcKZwf6zncuZcGCpwHYtm0z8+Y9AcCll17Biy+GzLNDhw7n8cfnMnXqd5g8eYpp\n39euXc3KlV/y9tshrdXv91FcHHIKHDNmHJmZmVGT/+bN33LNNT8EoHv3HnTo0JHCwoJIHbNUhRs2\nrOOii6aRmZkJNExnOGXK9MjzitWfTZs2cs011wLQp09fevc2P6mQm5vHsGGhs+mXXHI57733T+DE\nVoBZesJRo5pOcjF58lQA+vcfyJdffmFaJlZKSIDzz78wsvU0dep3uPPOX3DzzbexfPnnTJkyDQjl\nlf7d7+Zy+HABkiRFnGhPFW1OCMfSfgZkx544wbzO+LSG0YVjCbvGyBgE9ZBGu+t4bE3MrK85KX6O\nu+MduZI4XBMyqw7Jdcd14rHIBgOzm96fHpTjxjDgcLUNzQgVViQj8jzr11fk5IR2IpyONpvj3CRL\n0DvLWyeEzZA4VOWg3KNS44/eKoBo4WkmzBMZkyyF3u/AbPcpdfqTZZmVK79k6tTpZ2VKxGS58soZ\nXHnlDABeeukFcnJyARoI0u9+92p+85tfNahnGAavv/4qjz46n6ef/iM///kdHDt2lH/965/cdtsv\nou5Tf7841tZxuMw999zPtm1b+eabldx88yzTZPaGYfD443+kW7ceDT7fvn1rg9gDjWrF+JyYdQzD\nLP1hdJ1Y/YHY9eOXafi7WXrCRGJHWyyhv2FFkeviQUQTLyVk/fgSOTm5pKdnsHfvHpYv/5x77rkf\nCDmDjRw5mvnzn+TYsaP88pe3xe1TsrS52NFhQRnSqqTIBPl1Qbrp59tLnTHrbCg8MdHGm9QboyPz\nVUEmS/dnmba7qSgFTTfva0GVA6/W9GsrdlkjGnQsrdyrhTTypjCM0PhC2rcESGiGHHk+jRmU46ZH\npgeHqgEGDlWjR6YnppYWNt9rcZS/ZNuMh27A5uKUJp2bzAhr5fGo8ZuvbYtdVgwjJDwv6lHJ5O6V\ndM9IbEyxnpEiQ4r11HrdS5LEN9+s4gzndjktVFSUA1BUVMSXXy6PaKz18/D+979f0KtX7wb1Pvlk\nCRMmnE96ejperzdyXtvna2T/rmPZss/r/v2MwYOHATBkyDCWLv0UCKX6Gzo0pAkeOXKYwYOHcMst\nt5OZmUVJSXFdusQT733cuAm89947kXcQTqsYj+HDR0RSChYUHKK4uKjJzEKjRo1h+fKlVFVVAifS\nGTYmVn/q33P//r3s22ee6rO4uCiSDnHp0k8jWnEYs/SEBw4ciGonEZzOFDyeE88yVkpIM6ZNu5i3\n3vo7LpcrotW7XC5yckJ+BR9/vLhZfYpHm9CEw8LIIusxBWWsibOg0oZVNZ+M9pcp9EwDixJf2Jkj\nRbTKxhypsXHcbSGgN98bMyxITsaUqxshB7CC47aYgr+oxkq3DC8plhOCIFEtLZk92VhtajrUBpLT\nBLeXOjlSE0ubiP9c4mnlTeEJylR4VLIcQRQZ0mw6Q/PcaHrs53Q6PezjEYpE9+0pDebh85kHTTmd\nbT3wwG+iJZ0/AAAgAElEQVSorq5CUVTuuute0tPTAXjxxWfZs2c3kiTRoUNH7rnngUgdr9fLJ58s\n4emnXwDg2muv5557/h+qamHu3MdN71NTU80NN1wbOdoEoVSE8+c/xttv/4PMzCzuu+8RAF544VkO\nHy7AMAxGjRpLnz79yMvrwBtvvM6NN17HrFk3cuONN/Pss09xww3XYhgGHTt2avIozNVX/4Ann5zP\nT37yQxRF4YEH5tblmo5Nr169ueGGm5g9+2fIskK/fv1NY1DH6s/VV1/DvHmPcv3119SlGxxgep9u\n3bqzcOG/mD//MXr06MnVV1/T4LpZesJf//ouhg3Ljtt/MyZNuoCHHrqXr776kjvvvMc0JWRYy23M\nlCnTeO65p7jhhpsjn11//U/43e/m8vrrryZ9hCwRzvlUhtXVVSxZ8qFp2MrGE5hN0fFp5ibAkCkn\n1qwW+1qXNA/ndXTjD8KXhzLxxTjGcmqJ19cQDlXjoh6VdXvI5mbyzmlehuXVxhResepF94VmCYfm\n9guSF071F2JfFWTGdY5Kswa4oHt1zHHoRqjvsfaGY7+f5J9VrGfUI9OT9HGvphyzGhMMakyf/h3T\nPdBkOVPnhFsi/d0111zFX/7yj8i+6pmgtaX5M6O1jbFNOWaFJ9z9FXYOVZ2YwOILyOapFaW1VjYc\ngwqPtU7An37sSmgSL6lNbF+xf3s3QQ2Oe0Llw/u5R2pslHssDQNpJGA1aEjo3smeBY5nvjfrV2MS\ndaqKvRCLTY0/FGwk5vljCYbluZGgwfcrTJo12GBP+ATJPatTFZSjuaiqwurV33DZZVec9BlZSZLi\nnukVCNoqrUoIJ+qdfCrx6TJHa8y0RQOHqqPKeowJufnkpvojZkx3QOZgpT0ikOtrWY2fh13RSbEE\nqQ1EOw2FHDRISliZkahwiG++lxoIKjMTdKLCqbGwTtRSkcg4Bue6GzyzsONcjwwvh6rs9RZJoTEl\ne49TGZSjudTWuk65Wbq18d57p36fUNB2OKeF8IED+5k9+zZKSoqxWq3kDvkeWT0nJVTX0PxIsiW2\nK2MD4pUxv2ZTdC7oVomqNF4YRJeX0dHrHJ4ao0g6FtlosCdb7LKymZAWFW9fsbEp06spEEMIhTyf\nT9wjvrAK72BE9zdR4RBvr7o+hVW2KJNz9wxvQsIpGWe55oyj/j61OyBzoMJOsctKQZU9IpDznH7W\nHUtv1j3iPSO7Er1vfTrCoyqKwq5du+jatdspMUsLBIKGnNNC2OVycfDgAcrLyzEkma4XJp7r1eOq\nQJIVHGnJb/wngjcAL7/6Kpq3EkmSMVBQ7KlkdJ9ISu4ALI5MdH8NevVBLFWbIW8CpEeHDVRrD4Ak\ng+NEJCOfFvKSPlzmw1m2jGDAhyTJpKamkZqaSmpqKhabk6M1vaPai0UsJzEzOqb6qPRaTurcbqIO\nTpoh46k7eRDWjnWDmMLJKusENJr0DAcDq6zj180XRsmcP5Yk2HAsNepYUkGVg3K3atp+vHvUF6ax\nnlFAl9hR5ox7hC7ennMyAltVFdatW8Mll1wev6BAIEiac1oIDx06jG++Wc+SJR8SlJysOJSVcF1n\nRi6SHjA/WRe2zTb+WNcJeqtQ7ekgyXH3yTw1pWzfko8e9De8sOFrZNWKPaUd3tryyHVJWsTAC2+k\nQ5+x2FOz8brKKNq7ll2r3mLyT57FaSKrdEsmW4rsbPsi+pyhM6MDF/30BRLNsBdjyFG4Kwr5Zulf\n6TzuNixpHaKuHz+0gTdXrkCSZBRFQVXDqRhVLBZL5HdVVTGQUG39Cdo7YahNB+MPU1Qtk6FW4CF6\nAeXXZVYWZqJIBrkOF3ZFw6tFf80dashSsb3UyWGT7YRkzh9vLXHG3HJwBWL/ieWmnLiHpmO6tZCb\n4qdHpofCqoaWivARsZLiInRdx+s4oaXGMuVDfIc2w4gtmKuqqti9e2dM71eBQNA8Wo13tMXqYMXB\n+F6viZBqCRA0ZLwx2uma7qFnppd1R9Pj3qtHppuB7V0YhoFhGJFzhmFB4/V6I5Fb3G43mhZE1+u0\nOCwQdKNpGl4pHVfGxJA2bIKseejg+wb0ALW1tdTU1OBy1SLljceeNxxJTux5GIaOFOMe9XFXFVO8\nP5+eI66IulZZvI+v37oHw0hun1JWrQyZehvdhkxLqLyuBfny9TvoPvzS0KIlLQc5xjgri/eRmRdt\nEag4tJaybf9GliF74OWk5A5EsaWj+arxlu2mat9n6JqGJMnYbKEczTabre7fcM5mK4pq5Xj6RRhK\njJzEMVc3BqPaFWCVfOyraUeFP5UANsy05toj67Bk9cbqjA6e4a4qAiScGdERrmQ0bCoRYdu1nY7f\nHzTVrNOsAYK63KQmfdVVMxLOz90StDaP2li0hXG2tjHG8o5uNULYbnewudhJgYm3ajI4VI1sp5/C\navN2wkd/dpQ1HQoynmYRxswsaO5gFvv41JQelQ32FZs6WlR/wlUkI2YIzFj3s8kaPj1aw7PJfgan\nH0LTZRT8SEaQQCBAMBggEAji9/sJBoMEgwH8AR1dtuGwQorDjs1mp8ToRmUwnSA2CNZiyDYkJXpP\nV/NW4dr2d/SgH29AJ3vUrVgc5sdDAu5yqo9tIyV3ANaULLSAD5BQLDa8NWUU7V3Djv++hqSoUdaJ\nRDhhcTBfBOiahmySb9pdVcSK1+9gwKQf02vUd+Pew1NTii2lnfk9wgueBBZQEPp+Jvq+Gx+DMgyD\nzp27MH78xITqtwStbeKORVsYZ2sbY5s4otQz09uMs5sN8QRlOqb6Kaw2byfsTNO/vRtJVimqliOe\nx1mOAH3aeXBa9Cb36OKZBRMNfwnR+4rxnZEMumd4GZwbWiDU+GS+ORwdTzZc1mz8Ia9pc4Hj0y1s\nq+mFT6s3nnoxrOubXKvrTK5+Vcee6qdHjpteEnUOZl7sqs6OMpWDldH36d3BypBhobywtX6ZLw7G\nGgNYnFlcPWUIKVadLcUBDtecyCvszMij16jvMnnyFM7r7KayshZd11EUuc5qIWMYBn6/D5/vxE8g\nEMDv9+H3+wkEDWoML7pJvmIARavGUKK3SbTqg4wcOYqcwRfG7HsYR1p23XOPvmZXDSSJyL55UySz\n99/Ye1uSJAoKCujXr79w0hIIThGtSgg7LYl53MbDoepk2oNxvVL3V5zYt7MrOp1SfaiyQanbylcF\ntijt0uxcaKxzrroROnucKI33LpuK3NUryxsSihIcqHQkqQWH7ldaGysTkBQR0I3jJcfyEG/8bMJh\nGMP1DAOKXNaGgr1eWEe7qmNX9Ji5fgHWHU0Pxdz2mD/XMq8Th1PBMEJ/DuEQkUiQYtFJTU2N+0y2\nlsimi4U0a4BJfXR2lnmiF1t9++IJ9I+7gAgT3hs2O5NsUXTaO81NzCeLmfe2qirk56/l4osvO+X3\nEwjaIq1KCMfzuI0dQKEheal+rGrsdiyK3mAy9GoKR10NBUAsbSOsWYT/H6tM7PO5BjZZx6fHzppj\nVfS6RUB8j19Np4lEECHTpSRJBHUaaeokPOmH4yWbCZD6NA5/GbYUlNSGnodNCR35qW9NCJvy81LN\nBVSIEzG3YwW49wRlPAEJ3YBtJc6oRBXh7YVYnsaNkzLY6oKpDMl1xw3haZH1uJaFMOHn3jgxBIQC\ni7RzBOmR6WlwHjygS6YLLJkgOont6cby3q6srGTv3t306dMvoXYEAkFsWpUQhtgp3hpnSbIrOhZF\nN3VGidVOToo/KS21MfUTA8TSVsMCx2xiDnv0BvTY+8y7jjtjarf1tWZvUG4iEURdbGsjOoykmdCJ\nFQ7UE5QpSuCsrleT+e+hzMh7AKKCbBRUOZClaM3aruikWgJ4gkq9xUfiZleHquOwGKw97IwS5poh\nNRndKpFY2fU1/PpbEfEWXI09l4O6edmSWisX9ahscP9YPguHNi+j7/AL8RoxHMnq0c5hnvxXVVX2\n7t0jhLBAcApodUI43oQYKwGA2cRp1o43KMdJY9c0DlXHqoQmyHhlYpkew1q6FXPv43j7wYqk0799\nQzNuoqb7ck9Dzanxs4kXjznxyFsNQzoqkvkYzTTrcBCSbhkevAGJEndy7ygs9ItqYi8WjtVY6dcu\nZC5PRNDGI/aevxFZ7PXM9OKslxSjNhA/elY4MUR9Uz40XETK7gK2f/0WOb3GYE81+w6GLAXxwpqG\nqayspKKivFWkOxQIWpJWJ4TDxJoQG3/e1MRZ/3oygsuMvFQ/u47H9+AOT3iNwyEmkrA93kStGRJ+\nTcai6JFxJZoNKFZkp/rPJlZbuSl+jtbYknIICvc3Vl9iadYlTWjcjWNu13+unoAzrmXAp8ks3Z+F\nJIX61txsRvEWSuEoa1aTv8qmvnurj6RH9an+Qsmq6Byo6cIlt7yIocZeBHZI8VNUa4v8HivOtdVq\nZffuXYwbNyGRYQsEghi0WiF8Okg2jV3oOMiJCbt/ezf/PRQr00rIc9lsAm0qqlHYvBlPkzPb34tO\nMA/NjR4VaxsASNr5Kx7xNOum8ixnO/0xQ3w6LEYTDl5SKLRo3bZysgkrIn2M4zjn02QCumxq6Yj/\n3WtoRQhoUmT7ILxQCh1bs4DFEtNQb1d0qnyxcyE3jnN95MhhdF1HTjQijEAgiEII4SQxEza5KX5T\nzbV/ezd+7cRkX+tPzHO5von85MybJzCLANVY2DfOOhWvbmPMFg4AKw42L71bLOeyeN7ZdkWvO64T\nfU2RdAbXCUsz64eqQIe05PMEJ5vNKNH8zmbbJIkumhqbkQ0jsRjaGfYgxTF8HsysIZoW5NChA/Ts\nmXh4VIFA0BAhhJOkvrBR7U6C3hMTsKkHrHJi0mpqArYqIY0l2RjA8c4FJ2LKDgsls6xAXdvp9MlI\nXNMLt6XpUOFRm53Nqku6L6ZJPpZ3doe0aIeuMF0zfFia2EUIC6z63tFNOXglm80onkablxpazMX7\nDoS/exUeldVHzBNDNM5C1TMzdsILDANZCp01jiWAwdwaoigqBw4IISwQnAxCCDcTRYY0u0G1v+Fn\nTe0vx5uAdx1PLEdufZo6FzymUzXp9sQEhJk2m5XpoLq66bphzc2q6Ow6nlw6yXghE80WNrFM3/UX\nGsnup4fHPzTPzaAcN9Vemfxj6U0eH0om0UOYeP1PJE+yIkOWI/ZZ9voUu0IOZbHKBvweLDZnrNNb\nEWJZQ0pKSvB4PDgcp/6cskDQFhBC+AwTawKOt18cz+TZlHbdnHyziXr5QnTkr8TDYEYfwTHb/zbr\nS1N75snsp5uhyJDl1OmYgHk6mUQPTfU/mTzJifoneIKhfeaYZWMK39CFphYxFovKrl07OO+8kXH7\nIRAIzBFC+AwTawKOt18cz+TZlHZ9qvLKxqKx5hbbC/rEpG52BAcpccEfJt5iIZmFRCwaL5jCR3ca\ne0fHI17KwMZ9jGfVMPsOhO9dVGOtc0qL7VQ3KMeN1apSWH7C4tDO7uewbouqE2Z852qyHEHT71D9\ncRUWFgohLBA0k3NeCNvtDvr3H4iiqMiyFIn3Gwj4Q7F9AwGCQQ1ZlpHlE5mMNE1D1w10XcMw9Ei2\no/BPKD6wn0DAj6aFsiCpqoyqWrFY1EhCAq/XiyzLWCyWhFPxQfQEnKjDjhmJmGdPB/H3o6OJN6mf\njcRyNkvGYz2Z/f1kvwP1+7eluOmUjKO7Beid7mkwlmKXRBDz+4XfVX2BK0nR48pxehl59CidO3WK\n/UAEAoEp57wQtlqtjBo15rS1Hwrg70dVVZRG2XCys1MpLCzF5XJRUXEct9uDx+PG4wkJP1W1YLVa\nsVqtEcEfEv46Pp+/Lo2hi0AgiCRBtt1NoSs600aOw0Mw4MficJCS4iQlJRWHw4Gu6wSDoRSIeblB\nAkEddwCcFhlVsVBba+By1RAMBrFabZHFhdVqJS0tHYvF2mDhEAz68XhCYwgGNcDAYgGv14ckSRiG\ngdVqiRxJaWo/uj71J/VzjcYLpuZ4rCeyvx/PqpGTEtuqocgwrIMbVWl6IdZ4LF0yDdO417GcxFRZ\nbxA60xNUKKhO4e2lu7n7J0IICwTJcs4L4dONJEnYbOYmO0mScDgcOBwOcnJymn0Pj8eD2+3mO5LE\nktVH2XKgkgqXn3ZpNob1asc1F/XCaoleBCSCYRhUVVVSVFSE3W4nL69Dk040hmHg9XoJBPx06JBF\nRYUHRVEIBoMcO3aU8vLjVFVVIrvcOCwGnkDTFoCeuSrdunYhEAjULVQ8eL0edD1k4g0vBsKZNQ3j\nxP/D1ySJOmuGgiwrKIqMYUAg4EPT9LozqxK6fiKHsySB1WpLykpxsiSzt9uYaBNziNJaK1tLiKlJ\nJ3u2vP79qiorOVIl4UjLwWklrpMYJlozwIESP65aD6kpwkFLIEgGIYTPAsKCHOCGy9vhC2hUuXxk\npNqwNXWupgkkSSIzM4vMzOh0evHqhPuUkpKC2x3SnFRVpWfPXvTs2StS1r10N0vzD0e1Ybcq+AMa\nWWl2RvTL5odT+6AkENQhvB2g6yGhCkS2EML/mhEIBPD5QukFFUWuW7BIaFqQkpJiampqqKmpwev1\nYLHYcDjs2Gw2VNVCIODH6VSR5XICAR+yHKovyyqKIqPrOpqmYxghS0Zo0eSpE/BWdF0nEAigqgop\nKel4dSXu3q5kSQWtFp/PR8jaYImMV9d1+qR7CAYzOexKrVcvseAgye6FyxJM6mdlwfMLOFBYxO23\n3ES/3N5JbzV4ggrrN21n8sRRCdcRCARCCJ+V2CwKuVlNB9g/G/jh1D4AbNxdRkWNNyJ0Z1zQE5c7\nkPRCIqS9SklHYbJYLFgs5tmBMjKaDhiSbALxQCBARUUFx4+XYrVayc3NIzU1DUmS8AU0vilYzfFq\nX1S99ukOrrr8ImwWhUAgQE1NDdXV1VgsKhaLFavVgqbD+re2AtEJFMp9KXTp1gFVNuoJbo2ioqLI\ntklzuPyyy3jiifl8/NEi+v3q10ltNUBou6G89DCGMfKMWh0EgnMdIYQFJ4Uiy1w3vR/fn9w7Snt3\n2hJLmXcuYrFYyM3NJTc3N+qazaIwol+OqYVgRL/syPOxWCy0a9eOdu0aJkEoqXBT6TLPYOTy6vTo\nPTBqkWYYBgcPHmD//n2UlZXGXJDEonfvPgwaNJjt27exd+8eevbqm1Sc9LxUPwGfm8LCArp1657U\nvQWCtsw56CYjOBsJa+8naz5vLfxwah+mj+5C+3Q7sgTt0+1MH90lYjmIR0aqjXbp5n4IWWl2MlKj\nr0mSRM+evZg27TtcdNE0gsFg0n2+4orvArBkyb9R5JAzWGwMQme9NXpkehiU40ZVVfbt25v0fQWC\ntozQhAVnJadyX7wliGchaIpENelY5OXlkZ2dQ2VlRVJ9jtKGu/avS91pbl42O3JWUlKMy+UiNTXV\ntI5AIGiI0IQFZxWarvPW0t08+Mpq7ntpNQ++spq3lu5G008u8EZL0VwLwclo0gBDhw4/KW34s88+\nwWkJnVs2I9aRM4vFwo4d25K+r0DQVhGasOCs4p3lextogMerfZHfr5ver6W6dcY5GU0aICcnp9na\ncM+evdiyZTPlx0vISzXP0BUvGlthYQGjRo0RKQ4FggQQfyWCswZfQGPj7lLTaxt3l+ELaGe4Ry3P\nyey1N1cbnjJlGoZh8MUXyxmU46ZHpgeHGgreUn8POBaaprFv356k7ysQtEUkIxwR4QwRDGqo6rm3\nxyc4/Rwrq+W23y/F7BspS/B/c6bTMTvlzHfsHObTTz+lvLw8qTrBYJA77rgDn8/H888/j8PhIKiB\nJyDhsBgk8udrt9uZMWOGOK4kEDTBGTdHV1Sc3njGZ5Jkz5aei5zJMWoBjXZpNtPztVlpdjR/4LT1\npbW+y65d+7J371JUVSU93UF1tSehehdccBGLFy/is8+WMWXKtMjn7ngO0/UoL69h1ap8+vUb0Jxu\nN5vW+h4b0xbG2drGmJMTHZIYhDlacBYR9go2IxGvYEFDfAENQ00hMyv5kKoXXDAZVVX54otlkchl\nyaCqKtu2bW2WOVwgaEsIxyzBWUWsCFyJegULQh7m7yzfy8bdpZRX+8hMtZIm25lovhA3JT09ndGj\nx7J69Sq2b9/GkCFDk+6Hruts2rTxtCZYEQjOdYQQFpxVnKxXsCDaw7zCFaSCFOyFAfplJmaOBpg6\ndTqrV69i+fKlzRLCkiSxb98+BgwYREqK2MsXCMwQ5mjBWYmIwNU84nmYH66Q0RKwLGs61PplOnfp\nTp8+fdm+fStFRcea1R9VVdiwIb9ZdQWCtoAQwgJBK6LK5aPcxLENwB2QqfXH/pPXjVD+4BUHM/ni\nYCYrDmYy8pJfIEkyy5Z93uS9w8K7saA/evQIJSUlSY1DIGgrCCEsELQi4sWdBom1R9PZWuJENzkG\nFs4fHEraIOEJKrhtPRh56c9ZtWolZWVlpq02Ft5fHMhk47EUwse6LRYLGzcKbVggMEMIYYGgFRHP\nwxzAW5eXeHtpwyxM8fIHd+43EUNS+PjjxababmPh7dUUjtTYWbo/KyLwq6qqOHBg36kYokDQqhCO\nWQJBKyPsSb5hVynlNeam6WKXlYHZ7kjoyXj5g3XFQdeeA6my9mfZvjT8hgWHqpOX6qd/e3dM4a0Z\nciTk5ZBcN1u2bKF7954inKVAUA/x1yAQtDLCHua/mjk8Rv4j8ARlvPWErl2Nn6xh3MU30XPkVfgN\nK2FT9cFKB1tLUmIK7zDFLiuaDn6/j23btjZzVAJB60QIYYGglZKT6Yi5P+xQdez1hK4ih5IymLaT\n4sdI6WJ67bjbElN4hwkLfFmW2b17J4FAIMERCAStHyGEBYJWSrz94VynNyoLUqxkDT0zvTG1Xa8m\n094RP5ZlY4G/ceP6pMYhELRmxJ6wQNCKiRWBrL1+gFqXt0FZWQrt3Q7MduMNythVHUUOOW05VL3O\n8aohDlVncK4bVYHCKhuaES2s66c9lCSJgwcPMHDgINLS0k/9gAWCcwwhhAWCVkz9CGSK1YLmD2Cz\nKBw8aGXt2m9Q1egpQJEhxRptqo6VV9iihIR3//ZutpU4Oe6x4gnKEeetxmkPVVVlw4b1TJ485dQP\nWCA4xxBCWCBoA9gsCjnZKZGsND169GDLlk0EAomlRQoL0mKXFbcfPDWldGsnMSjnxBRiUeC8jm40\nvaEmbUZx8TGOHDlM587me80CQVtB7AkLBG2UPn36JpwhKWyqvqhHJd309ax4/Q42fPoison7dViT\njiWAAVTVwoYN+c3K0CQQtCaEEBYI2ij9+w9AlpOLza3IMGxgL3r16M7mzd9y6NDBZt/f5/OxadPG\nZtcXCFoDQggLBG0UWZbp2bMnhmESwzIOkiRx1VUzAFi8eFHC9RpH25Jlmb1791BdXZXU/QWC1oQQ\nwgJBG2bw4KHNMgn37z+Afv36s3XrFvbv3xu3rFliiHA4S0VRWLt2TXO7LxCc8wghLBC0YaxWK126\ndG1W3RPa8Idxy5klhqgfv7qiopx9+/Y0qw8CwbmOEMICQRtnyJChzYpi1bdvPwYOHMSOHdvZs2e3\naZl4iSHC4SwVRWHTpk0ikpagTSKEsEDQxklLS6djx87NqnvVVd8DQnvDZnvL8RJD1I9fbRi6iKQl\naJMIISwQCBgyZEizNNFevfowZMhQdu/exc6dO6KuN5UYIhzOMhxJq7KyIuk+CATnMkIICwQC2rVr\nT25uXrPqfve7VwOwaNHCKG04XmKI+uEsIRRJKz8/v1l9EAjOVYQQFggEAAwZMqxZ2nC3bt0ZNWoM\nhw4dMD33GysxRONwlgDHj5dy6NCh5nRfIDgnSVgIu91u3O7oPxqBQNA6yMnJoX377GbVveqqGUiS\nxL///UHUkaf60bam9Kjkoh6VDMl1m0bbslgsbNq0UUTSErQZmhTCBQUFzJw5k3HjxjF+/HiuvfZa\nCgsLz0TfBALBGWbw4CEEg8Gk63Xo0IEJEyZx9OhR1q0zP/ebSDhLAL/fx5Ytm5Pug0BwLtKkEH7k\nkUeYOXMmmzdvZtOmTfzgBz/g4YcfPhN9EwgEZ5iOHTuRmZnZrLpXXHEVqqqyePGHzRLkYWRZZs+e\nncLyJmgTNCmEy8vLueaaa5AkCUmS+P73v095efmZ6JtAIGgB+vcflNTecDgcZWZWNhdcMJmyslJW\nrVoZt6w/2DCEZWNkWSE/f21zui8QnFM0mcpQlmX2799Pr169ADhw4ACKklzQd4FAcO7QvXt3tmzZ\nhN/vi1tON0LRsIpdJ/IH9500i1WrvmbJkn8zdux47HZ7g7JFNVa82om1f/2cw433iIuKjnL06BE6\ndWreGWaB4FygSU34zjvv5Prrr+emm27i5ptv5vrrr+euu+46E30TCAQtRLdu3ZtM7GAWjvKoO5Np\n1z5IdXUVn3/+aVRZrxYqG/5pHMKyPqpqYf369UknmBAIziWa1IQvvPBClixZwubNmzEMg/POO492\n7dqdib4JBIIWYsCAgezatR1FMZ8i4oWjdOYOIDMrh88++w/nn38h6RlZMcuGKXZZGZjtjnLa8nrd\nbNmymWHDhjdrHALB2U5CR5Tat2/PlClTmDp1qhDAAkEbwGq10qFDp5jX44Wj9GoKl3/3BwQCfv79\n7w/ilg1TP4RlfRRFEU5aglZNzL+MG264AYDx48czYcKEyE/4d4FA0Lrp27d/TActq6KjSOZmYoeq\nM2HMCDp37sLq1asoPXYwZujK+nXsMcpIkiyctAStlpjm6CeeeAKA999//4x1RiAQnD106NCBlJRU\nUwetXcedaIb5Gj4v1Y9Flfmfa37Eq6/9g0WLFjL9R/05WBm971u/Trzzw0VFRykoOES3bt2THodA\ncDYT82ufm5sLwMcff0znzp0b/Hz88cdnrIMCgaDlMHPQircfrEg6fdu52VripMg2nik/fYFOE2dT\nWlpGj0wPdiUUujL8Ey+EZX1U1UJ+/hpqampOzcAEgrOEJveEzQSuEMICQdtg4MBBaJrW4LN4e7ya\nIYnjHPwAACAASURBVLGjnte0JCs4M/KotXTD0HWm9AyFrry4V3mTISyjkfjqqxUipKWgVRHTHP31\n11+zcuVKSkpK+OMf/xj53OVyndQNs7KcqGrrOWeck5PW0l047bSFMULbGGdzxjhoUF+Kiooivzs1\nSDlqUOuPlpxOi0G512baTmGFwbjeDrJO4s9f13V27vyWyZMnxyzTFt4jtI1xtoUxxhTCFouFlJQU\nJEnC6Tyxl5Obm8vPfvazZt+woqL1eDnm5KRRWtq6zWNtYYzQNsbZ3DFmZ3dh27bdWCyWE205JWr9\njqiyWXY/R2rMhXAQO4eOFJOTGQrgoekhrdquNh1Puj5bt+5CVZ306zcg6lpbeI/QNsbZ2sYYa0ER\nUwiPHTuWsWPHcvHFF9OvX7/T1jGBQHB206FDBzIzM6mtrY18Ft7DrR8tKy/VT//2bso9lroAHg3x\nuspY9umHzJz5w6hIW7GiZpmhqiqbNm0iK6s9OTk5p2ycAkFL0GSwjn79+rFy5Up27NiBz3fCS3L2\n7NmntWMCgeDsYcyYcSxd+nlkKymcnnBgtjtKm81L9XOwMlpLrjq6lY1fLqXnmB9QGjxxPRw1C0Jt\nJoKiyHz11ZdccsllpKSknOToBIKWo0kj0JNPPskrr7zCa6+9RklJCW+//TYHDx48A10TCARnC+3a\ntadXr15RTlFm6QkH5bjpkenBoYY8ocMe0GN7yiCrHK023xQudlljJnQwx+CLL5adVMYmgaClaVII\nf/nll7z66qu0b9+exx57jIULF1JVVXUm+iYQCM4iRowYhdVqvt9bn7CWfFGPygYe0MOGDmPQsLEo\n9gzTerGiZsXD5/OyYsVyEV9acM7S5DfearWiqiqSJBEIBMjLy2vgKSkQCNoGiqIwcuQogsHE0hw2\n1pIlSeJ7l1+Mt6bMtHy8qFmxkCSJysoK1q5dnVQ9geBsock94ZSUFDweDyNGjGDOnDnk5ORE0pMJ\nBIK2RZcuXenYsTOlpSXNqt+1S2fW7j8IGXlR15qKmhULWZYpKDiEw+Fk2rTzm9UvgaClaPIr/6c/\n/QlFUbj33nvp3bs3kiTx7LPPnom+CQSCs5AxY8adlPn3osHpHNi4GG9NKclEzYqHoijs2rWDHTt2\nNLsNgaAliCuENU3jmWeewWq14nA4+PnPf869995Lp06xs6sIBILWjd1up3fvPs0WxO3bt6O9to/l\nf/sFWZXLk4yaFRtFUcjPz+fAgf0n1xBgGAbr1q1h9epV7Nu3p8HJEIHgVBJXCIdWl7vOVF8EAsE5\nwuDBQ0+q/qWXXo5kaHy65F0wtCbLazrU+uUmvadVVWXdurUcPlzY7L4ZhsHXX/+Xgwf3c/ToETZu\n3MDChe/x0UeLWbt2NcXFRcIRTHDKUObOnTs3XoHCwkKWLl1KVlYWNTU1lJeXU15e3uy8wm63v1n1\nzkZSUmytajxmtIUxQtsY56kco6IoeDxuKisrkKTkVViHw0lFRTk7dmwnL68DXbp0MS2nG7Ct1Mm2\nkhT2lDs4Um3DHZDJdgYwu63NZiEQ0CgsLMBqtVJdXU1paQmlpaVUVVWSmZkVt7+6rvPf/66gpKQE\nRQkdpZIkCVVV0LQg1dXV7N+/lz17dlNZWU4wGCQ9PR1ZbsZm9kkgvq/nHikp5icLmnTM+uijjwBY\nsWJF5DNJkli2bNmp6ZlAIDgnGTp0OAcO7EOSmhcM+tJLr2DVqq/5+OPFjBkz1lSQba9LBhEm0cAe\nsiyzceMGJEmKtKvrOps2fUu3bt0ZPHgoDkfDgCK6rvPFF8uoqCiPK1QtFiuGYXDs2DEOHy5k3bo1\nZGdnk52dS2pqKPWj3x8gGAzg9/vx+Xx1n/nJyMhg0qQLm7VwEbROmhTCy5cvPxP9EAgE5xhWq5Vu\n3Xo02/SbnZ3NhAkT+frrr1i/fh1jxoxrcD1eysRil5WB2e643tRhTTZMWLAWFBxi37595ObmIssy\nmhYkGAzi8XgJBPxJabWKEppCKysrqaiowDCMuPWPHTvKunVrGTt2XMwygrbFmbWhCASCVsXQocMJ\nBpve043FpZdegSzLfPTR4qhoXPFSJjYnsEd9LBaViopyjh8vo7KyEpfLhaYFT8qsXF/rjoWiqBw8\nuJ+tW7c0+z6C1oUQwgKBoNk4HA66du3a7Po5OTmMHz+BoqJjUQE37KqOI0bwjnBgj0Qdts4mVFVl\n+/ZtHDiwr8HnwWCQ6moRjbCt0aQ5WiAQCOLx/9u78/ioynt/4J8zZ/bMJJNlspKFJAwhQAhbAsi+\niVALaq9bBa/aa2+91q29FbG1tVqoWvvqrf1V21ur1uuCWnbUagibsiWA7BAChJAA2fdllnPO749x\nhixn9j35vl8vXy+TmTnneSYh33me5/t8n/HjJ6C6uhoymXd/TpYt+y4OHTqIzZs3YNKkG6UxWYnj\nwyASo0w40zj4JKZp2sjIWpZKrdupOjo60d3djdbWZrS1tYPnOaSkpGLy5CnQaqND3UwSBDQSJoT4\nRKPRYMSINNHH3BmpxscnYMGCxWhpacGOHSX9HnN0GATDAFWtqm+PTGTsCVtHrshE7xGOWFaC8+fP\n4erVGnR3d0Mmk0KhUKC5uQmfffYpysoO0uEUw4DLj66PPfbYoEw+rVaLwsJC3H777UFPzSeEhJ/x\n4yegpuZT+2iYF+DRmcFLlizFvn178fnn23HTTTMRHW095EHsyEQA2FWlE21HbSuLnGh4Vf4yFBxl\nSUulLKqrL6O6+jJSU1MxcmQOkpNTgtw6Egwuf1X1ej2uX7+OyZMnY/LkyairqwMAfPbZZ1i7dm3A\nG0gICX/R0TFISbkRJGxbiwaOVE83qEVfr1KpsGzZd2E0GrFt25ZBj/c9DMJZwlaXifEpYSucMAwD\nhmFw7do17N69E1u3bsLhw2XgOO8T4Uj4cfnbevbsWbz77rtYtWoVVq1ahXfeeQcXLlzA66+/jn37\n9gWjjYSQCDBu3HiYzWaXW4scTU3PmjUbyckp2Lt3N65everwPs4StqLkgscnMUUCmUwGs9mMqqpL\n2LZti9P3h0QWl0G4sbERMtmNdRapVIqWlhbI5XLI5eL/0Aghw09cXDwSE5O83lrEslLcfvv3IAgC\nNm782OF9bAlbYtJ0XMRMRXuDYRhwnAUlJSU4dOjAoG1dJPK4XBMuKirCww8/jOXLl4NhGGzZsgVT\npkxBV1cXBWFCSD/5+eNwvW4HVFL+26no/lydGTx+/ASMHp2HEyeO4+zZ08jLyxe/z7cnLg1cc56U\nLqCzwz99CWcymQzV1ZdRV1eHxMREaDRaxMbGIj4+AQqFeHlEEp4YwUUlcrPZjA8//BCHDh2CIAgo\nLi7G3Xff3W907ImGhqHzL0Sv1w6p/ogZDn0Ehkc/g9XHkpIvsK/SLLq1KEvX47TcJGCtaLV27a+R\nkZGJ1at/7jT5k+NhT9hiJUB0tArt7T0+9yHcifXTbDaD43hotRokJCQgMTEJ6emZETtYGmr/JvV6\nrej3XY6EZTIZVq5ciZUrV/q9UYSQoSc/fywaGncDGDxSdXRmcN9gmpGRialTi1FWdhCHD5dh0uTi\nfoG2L1vCFrH+rZbJrMH42rVrqK2tRVnZIej1eowcmY2srGzazRKGXAbh5uZmvPDCC9i/fz8YhsGM\nGTPw7LPPen2KEiFkaEtNTUNcrA7j5F39thaJrdU62sr03e/ehqNHj+CbWhna4mLQy7Eutzn5YuCI\neiiQSCSQy+Voa2vD4cNlOHbsG4wYkY4JEyZG7Oh4KHL56/bLX/4SWVlZ2Lx5MzZu3IjMzEw899xz\nwWgbISRCjRplAMdZ+m0tEuNoK1OdkIkFdz+LlPxF6OWkcGebkzd4AThZr8auKh12Vumwq0qHk/Vq\n8JFReMttLCuFIAiorr6Mw4fLQt0c0ofLIFxdXY3HH38cSUlJSE5OxmOPPYYrV7w/MJsQMvSNHJkD\nmcz5aMvZVqbrHXJokseKPjZwm5Mv9aM93c8c6RiGwZUrV9Dc3BTqppBvuQzCPM+jqenGD6ypqYnS\n4gkhTjEMg4yMTKfPcbaVqZeTwMiLr5bZtjmJjWLLq2Vuj2K93c8c6aRSFkeOlIe6GeRbLteEH3ro\nIaxYsQJz584FwzDYvXs3nnrqqWC0jRASwcaMGYvz5887PNjBVnRDbCuTkuXBMHC4zUkm4XH8uho1\nHTcysHssLCrqWZhMapcZ2IB7RyWGMukrkOvUzc3NqKq6hKyskf69MPGYyyC8YsUK5Ofn27corVq1\nCrm5ucFoGyEkgqlUKiQnp6CpqUH0cWenJCVrrcU4xB6TSnjsuaxDLycemeo65RiT0O0ycDn7EOBq\nP3MgeVp3uy93A7dUKsXx48eQkZHZL2NaEASYTCbaaxxEbp09ZjAYYDAYAt0WQsgQYzAYsGfPVYd1\nBQYW3VCwNwIOxwMWDmjskaPHBPR0NCAuRosOU5TTe7o7inX2ISBJY7IHsWBnTtvWqW1s69QAHI7w\nvQncJpMRJ04cx4QJhRAEAefOnUVFxTkYjb2YMKEQBkOe3/tGBnMYhO+44w6HJ3wAwCeffBKQBhFC\nho6UlFRotVr09vaKPi5hrIFYEIDrnXIYOQnqO+XY2yOFhbdOFytZHlo04Iv3f4b59/8BMrXzIOzJ\nKNZR5a18fbdPI1JvuVqndjTC9yZwSyQSnD9/FgzD4NKlizCZjGBZFizL4tixb1BTU4MZM2ZCqVT6\n3jHikMMg/PTTTwezHYSQISorKxunT590WCjidIMal9tuBJBejkUvx/b7GpJkzPju42CVMS7v13cU\n64rYUYm2156s9zyw+cqbdWpvAzcASCQsKirOQiKRgGVvvOcsy6K1tQXbt2+BwZCHtLQR0Ol0fin2\nYTRzaOs0IkajgEI2eClguHEYhIuKioLZDkLIEDV6dB5Onz4p+pizADJQXNo4tLU1Qh2TJPKoAJWU\nR3ocj9wYzwPkwMpbvgQ2X3izTu1rgpmzwMowDM6fP4eTJ49DKmWh0UQjNjYWBQWFUKs928bF8TzW\nl1biaEUDmtuNiItWYKJBj7vm54IdxpW8hm/PCSFBwbIs0tMzRB9zFkAGMkMOqalR9LE0rRFzs1ox\nJcPsl6lib0+C8pWzE6IcjfCdHe3orwQzhUIBlpWip6cbV6/WYvv2LTh8uAwWi8Xta6wvrURJeQ2a\n2o0QADS1G1FSXoP1pZU+ty+SURAmhATc2LHWs4YHchZABj2X5ZGZMQJmYxes584IYBkeWboeTEju\n8uvINBiBzZF8fTeydD1QSTlYR/gcsnQ9DutuexO4fcWyLKqqLmHLlk04duwbtLe3OX2+0czhaIV4\nlvzRikYYzZz/Gxkh3MqO9qfYWDWk0qGzDuDoZIyhZDj0ERge/QxVH/V6LcaMyUV9ff2gx9LjeFTU\nu/6boJQD13piIOuze4YTGMjlUuhibqzdRkcPznb2hqN2pcfxiNX55x6OzIgRYOGM6DEzUMkEWP9k\n9r9n335O0wqQXzGjtpVFl4lBlFxAmo7DpHQBEiawba2vv4KamgtQKBSIj49HUlISCgoK+iX2Xmvs\nQnOHUfT1LR29YOUy6BMGJ9wNh3+TDoPw7t27nb5wzpw5Xt2wpSUwCQ2hMNSO2hIzHPoIDI9+hrqP\naWk5OHv2AqTS/n92cmN6YDLdyEJWsjxkLG/PjlZJeeijTGjoEl+jvdIsQU50j9+PMhzYLlt2dG5M\nN9rb/XILt3SLDHLF+mnQ9SAnuv92qmCerWw09qC9vQaVlZdx6tR5zJ073741jTNziNMq0NQ+OBDH\napXgTOZBv5uh/n31N4+PMvzb3/4GADCZTDhx4oR9n3BFRQUKCgq8DsKEkOEpPj4e8fHxaGvrP3Xp\nKEO57/7cXosE1W3iW2UGJh/5a1+vs8zpQPO2D2IJZsFuO8uyaG9vw2efbce8efOh1UZDIWMx0aBH\nSXnNoOdPNCQM6yxph0H43XffBQA89dRTWLNmDSZMmAAAOH78ON55553gtI4QMqQYDHnYv//rQaNh\nYHAA6fu1s6xhmDuhlPLgBaC8WoYrzQr7iDpebca4xC748jc+mGcW+2tvcij2OPfFMAwsFjO+/PJz\nTJ8+EykpqbhrvrXS4tGKRrR09CJWq8REQ4L9+8OVyzXh8+fP2wMwABQUFKCioiKgjSKEDE3p6Rk4\nceKYw+IdjjirbnX51B60JufiGpeBqtYblbl6ORa1HSyud8qRHmMMWgDyhTdFNwJ5Hd8x+OqrPcjP\nH4exY8fh3oUG3DEnh/YJ9+FygkKlUmHz5s32r7ds2QKVKrAL/YSQoct21rCnxLKGlT0XcHLnm1j/\n0ccO9/VygiQijif016lO4XY6FMuyOHXqBHbv3gmLxQKFjEVirJoC8LdcBuF169bhnXfewfjx4zFh\nwgT84x//wLp164LRNkLIEJSba4Bc7nkpRNsa7dysVszLasXcrFYsKIiBwWBAZVUNeszOh7nhfjyh\nv/Ymh2qPszNSqRSNjQ347LNt9rOMOY5De3sbqqqq0NjY+O22M3FGM4f6lu4huZXJ5XR0Tk4ONmzY\ngM7OTgCARqMJeKMIIUMXwzDIyhqJCxfOe/X6/mu0DO69dyV+s+436O1sglKrd/i6cDie0Bl/neoU\nrqdDMQwDs9mMHTtKwLIsjEYjBIGHXC6HxWKBVCpFbGws4uLiYTDkQa/XDosqWy57IQgCPv74Y7zx\nxhvQaDSoqanBkSNHgtE2QsgQlZ2dA5NJvMCEp5KSknHzokW4WrHf6fNCGYDEcDzQZZLYR+f+KroR\niuIdnmBZCQABCoUcSqUSEokEcrkcEokEbW1tuHTpIr744jPU1dUNiypbbk1HHzhwACUlJQCAqKgo\nrF27NuANI4QMXVqtFlpttN+ut3jxLWg69zkuHt4KBoMrcwGBDUADA6ozvGA9HGJXlQ47q3TYVaXD\nyXo1eMHzalmO+Os6oSIIAj7715coO31N9PGhVGXL5XT0wYMHsWnTJtx2220AgNjYWBiN4pVPCCHE\nXXq9HjU1V/xyLZlMhu/fex9+//uX0XF5N25Z9QJaegdvz/E3d7YCDdyr6ypz2R97k0O5x9lfTJwU\nbd0WAIPX+ls6etHWaURibHgn27nDZRBWKBT9yo/xfPhM5xBCIld6egaqqi6J7hn2hsEwGnPnzsWu\nXbvQdGoj5i+42W8ByFHRC2cBNV/fPShAO6v81fd0Jn/tTQ7mHmd/U8kEh2vbsVolYjQKkVdFHpe/\n/QaDAVu2bIEgCKipqcFf//pXTJ48ORhtI4QMYcnJKX4LwDb33HMPysvLsXXrJkyaNBnx8Qk+Xc/Z\nSFcQnG8FEgT0Oye5x8Kiuk0FQDwLONwTx4JNyjreGz5x1NCpsuXy8+Hq1atx6NAhNDQ04M477wTP\n8/jv//7vYLSNEDKEMQyDhATH2cze0Gq1+N737oLJZMJ77/3D55k720jXOhpj7CPd0w1ql1uBrrt5\nTrJNuCWOhQOxte3M6C7EchfAccNkTVij0eDFF18MRlsIIcNMUlIy6uvrnB4s76ni4ukoKzuEU6dO\noLT0SyxceLNX13FV9MIQ1+1wulTB8jBynvUpHDKXw40gACN1vTDEdcPM31gO6O4SUFLyBRYuXAyW\njewRscsfeXFxMf7v//6v3/fuu+++gDWIEDJ8jByZ7fcRDcMwuP/+BxAdHY2NG/+J6urLXl3H1UjX\nzEucbgVydB6xkuWRGRO5mcvBYKsDbssg31utw6VWJWzpSQzDoLOzAyUlX4j+/kTSKNllEI6JicGX\nX37Zb1uSrXAHIYT4Qi6XIzY2zu/XjY6Owf33PwSO4/Dmm3/xakeHreiFGNvU8ej4bozQDg6o4xK7\nHQboZK0J45P6V/4alxj+da2D6XSDGhX1MtFlAJu+gdhsNuPKlWocOLAfW7duwubNGyMmELsMwlFR\nUXjzzTfR1taGH//4xzAajf2ypQkhxBeJiYkBue7YseOwcOFi1NXV4aOPPvD49c6KXiRGmXCmUY09\nl3Wo6VBCEIA0rRGzM28E1Hx9NzJjeqBgxUe8tsxlmoLuz5Pa17ZAvGHDR9i37ytcvVoDs9kMQeBx\n6NCBILXYN25VzJJKpXjppZcwatQorFy5Eh0dQ+egZUJIaGVmZvWrnuVJ4QtXli+/HenpGfj66704\nfLjM49c7KnrBMOiXsGU9sUmJc03WkZotq7q+Sw4jJ4GCtW5PioSTnELN09rXDMNAJpNDJpP1+151\ndTWuXbsa0Lb6g8sg3PcYw8ceewx3332300LbhBDiidjYOKhUaqeVpLwlk8nw7w/+J2IS0vH+hx+g\nsbHBo9eLHRoxJqHb5UhtYFa1kbNuTwr3k5zCgTvLAO6QyaQoKzsY9tPSLoPw888/3+/r22+/HTt2\n7AhYgwghw09SUqLT7UDusnA3RtG2oH62ZzRmrfwfFN/5MrYdrIPJLH6MorMReN+pY1cjtS6TJKyO\nEow0/qx9bTKZcPhwuZ9aFhgOtyi98847uP/++/Hyyy+LPv6zn/0sYI0ihAwvickjcL2zS/SxvpWk\nHLFN/zZ0K9BlUkEl5SGV8OgwfTtFyQDqmCSoY5Lw6aFjWHHTiEGvdVZ6si9XpxSBgcvpVCrI4Vy+\nvhtyuRRXmiU+lR6VSCSoqrqIjIxMJCcnB6i1vnEYhBUKa0kwtZqmTwghgRUVrUevRfxkHHcCl1j5\nSEB8/6hJnoRjJ05gwvjxDl/bt5bzQLaRmlglpySNCVGy8DxKMJJIGGBKhhk50T0+lx6VSqX4+us9\nUKmioFQqoFKpoFAoIJVKIZXKwLIslEoV0tPT/bpf3e32OXrg7rvvBgA8+uijQWsMIWR40mkV0Col\n6OgdvADsKnA5y6YVo9Qk4KOPtiBjRBqiY+KcTh07GoHbRmSORs/OgvRQyYZ2VE/bn/xV+5phGPT2\ndqO3txutrS2DHuc4C44elWPkyByMHTvO7+VUnXF4J0fT0DY0HU0I8ReFjMWUvGTs/Gbw0XWuApez\nNVoxrNCLloYa/O1vf8HD//W0V1PHrk4pchakxbgKaMEIeO7ydPo+ErCsFDzP4/z5c6isrEB2dg4m\nTgzOGQkOg7BtGrq6uhplZWVYtGgRAKCkpARTp04NSuMIIcPHvYtH4/r1a6hq5DxaB3S2RismM0GC\nSYUTUF5ehvXv/Q05C37m0dTxwIDoTZC2cRXQHD0+TRu6HSqeTt9HEtt09IUL50MfhG3T0KtWrcKG\nDRsQGxsLAPjRj36Exx9/PCiNI4QMH6xEgnsXGbB7zx5wjMLtUZ+zNVqt3AwLPzC5pwej738IbW3t\nOHL4IOJHH4MiedKg1w4cgXszAnQ1neoqoDl6XH7FDIOux/79YI2UXRXScJVARwZzOfHd2NhoD8AA\nEBsbi8bGxoA2ihAyPKWljYBWEwWzWXyLiiO20bI1O5oZdOTgwAAlkcnwyCOP4tVXX0bJBy9i2b+/\nAJU+z+kI3N8jQHcOiHD0eG0ri5xogGGCOzXsTiENyvz2jMvPLLm5uXj22Wdx9OhRHD16FL/4xS+Q\nm5sbjLYRQoahjIxMjwsC2aZ/l43tHVSP2VF5SJVKjR//+AnExcVh21vPQla7zWEtZ09KKbrLVUBr\nN0odPt5lYtBrkfhlb7Un/FVIg9zgMgivXbsW0dHReOGFF/DrX/8aGo2m32EOhBDiT2PG5Htd5UjK\nelaPOSZGh8ceexIajQbvv/cWdn25CRJm8AcAT0spusNVQItWWBw+HiUXIJPwQS8K4s9CGqHiz7Ko\n/uB0OprjOGzfvh1PP/10sNpDCBnm5HI5kpNT0dhYH5T7JSUl48kn/xv/7//9EVu3bkZtbS3uW/kQ\nBFZpn8J2VaDDmxGgq/3Gcqnjx9N0HMx8aKaGPc38DhfhmtXt9HMLy7JYv359sNpCCCEAgFGjRnm8\nLuyKsxFQWtoIrF79c+SOGo1ebQE+OyPrV7+aYQIzAhx4QISS5ZCm7cXo+G7Rx20HSExKN/s8Nezt\niFCsnnYkHMUY7Kl7d7lMzCouLsbnn3+OJUuWBKM9hBCClJRUKJVqcJx4nWdPOBsB9U3aio6Oxs3f\nfw6X26Psr+2bfBWIEaAtoI2O78apejWaeuSo7VCguUdmv7bYVicJo3I5knb0wcBfI0J/FdIIhnDO\n6nYZhDdu3Ii33noLSqUSKpUKgiCAYRjs378/GO0jhAxTycnJqK2t8fk6jrKam7ql/bYv6aNMaOh2\n/ofanb2/3jjXpEZNh+PMa0cBz5sPBkN5n68j4ZzV7TII//Of/wxGOwghpJ/s7FxcunQRcrn7JSkH\ncjYCsh/uAGsgqm5TARDPyu42AS2dRiREK/w+AvRllOZuURB/3CvceLI3OhBr+v7iMginpaUFox2E\nENKPXq9HVFQUzGaz19fwtKSlIz0dDfjjq7/BA//+7xg5Mtvn6/Xlj1Gaux8MwnlE6C5vi6aEaz1v\nl0H42rVreOWVV3D27FkYjUb79+lMYUJIoKWkpKK6+rLXr/e0pKUjCksDGuqv4pVX1mHZsu9iyZKl\nYFnfrulOG/09SgvnEaG7vJ1OD9esbpdBeM2aNVi6dCnOnDmD3/3ud/jggw+QkZHh9Q1jY9WQSv3z\nyxsO9HptqJsQcMOhj8Dw6Gek9bGoqBB1dVfsR6u6Izq6/2gnPY5HRb17f3PUMgFpOguutrHWylsy\nAWk6DlMm52DSiDV4/fXXsXXrJpw9ewqPPPIIkpKSPOqPI47amB7HI1Y3ePQGDO5nIO/lCwsH9Jit\n76Wnf/oH9tHCWauiiWnoVkAd5fweM2IEWDjjgPYM7jPHcUH7t8IILkrTrFixAps2bcKtt96KrVu3\ngud53HXXXfj444+9umFDQ4dXrwtHer12SPVHzHDoIzA8+hmpfdy6dbPb25Wio1Vob+/p9z2x6Uup\nhO+3JmyTpetBvt6aqXy9Uw4j13/E1NPdhQ8+eBfl5WWQy+UoKJiAceMKMG7ceGg03v/R9nSKwT5q\nSQAAIABJREFUVayfgbqXt3y9j1gfu0wS7KzSARC7gIB5Wa1+mU7neQ7f+97dPl+nL0dB3eVIWCaz\n/qKq1WpcvXoVCQkJaG5u9mvjCCHEkdTUVFy+XOX168WSl5zVXD7doMblNkfTncBDD/0Q48dPwJYt\nm1BeXoby8jIwDIORI7MxatRo5OTkIicnB1FRGp/a6O46paeHN/hyL08EIgt7KEynD+QyCE+ZMgWt\nra245557cPvtt0Mul+Pmm28ORtsIIQS5uaNw/vw5yOXuT0mLGZi8JBaI3MseZlBcPB1FRdNw7dpV\nnDhxHCdOHMPFixdQVX0Few8cQW9XM5L0CZgxYybmzJnrdts9ybz2daQZyH2+gcrCDucEK2+5DMK2\nkpUrVqxAUVEROjs7YTAYAt4wQggBrPWdtdoYGI29fr/2wEDkSfYwwzBITU1DamoaFi2+BcevyXG9\nQwozFLD0tOJqxX5s2PAmvvjic9x88xLMnu1+MHZHOO/3DWQWdrgmWHnLYRCurKwU/b5EIkFlZSWd\npEQICZrU1FRcunQx4PfxdrrzdIMaNZ0qgLGuVsrUccgsXIb09AyUrv8NPvnkI3zxxeeYM2c+pk+f\ngbi4eJ/aaeHCe79vIKeNgzWdHiwOg/DDDz8MhmEgCAKuXbsGjca6vtHZ2YmUlBSUlpYGrZGEkOFt\n1KjRqKysAMu6nLzziTfTnc6mXqOS8vHrF17GztJ/obS0BFu3bsK2bZuRlzcGM2bMRFbWSEilMshk\nUkilMigUCjCM67nkHjMTkJGmp+vLjgRj2jiSymY64/A32hZkX3jhBUyZMgW33HILAODzzz9HeXl5\ncFpHCCEANBoNbrppNvbt+woSSWCHPZ5Od7qaemUVWixffjsWL74Fhw+XYd++r3HmzGmcOXN60POj\noqKQlZWN7OxsjByZjZEjc6BSDQ5kKpng15HmwPVlJcsjXm3GuMQuyLzcURrMaWN/fXgIBZcfK8vK\nyvCLX/zC/vWSJUvw+uuvB7RRhBAyUGpqGqZNm4H9+7/2W6EMMZ5Od7o79apSqTBz5mzMnDkb169f\nw6FDB9Ha2gyz2QKLxQyTyYS6uus4deoETp06AcC6O6WwcBKmT5+BvLx8+wcQKevfkebA9eVejkVt\nB4vrnXKkxxhdJnuJBcFgTBu7m5wWzkHaZRAWBAHl5eWYMmUKAKC8vBw8H/lTAISQyDNiRDqKi6fj\nwIH9AS/64+50pzdTr8nJKfjud1eIXq+9vR1VVRdx8eIFHD16GGVlB1FWdhA6XSzGjLEGYpmMhcnE\nIWbUYqj1eTDyUq9Hms6m0zlB4jTZy50gGMhpY1fJaeF6hnBfLoPwL3/5Szz11FP2KRGj0YhXX301\n4A0jhBAxGRmZEAQBhw4dCOiI2BP+nHqNjo5GQUEhCgoKsXz57bh06SL27/8a5eWHsH//1/2fvG8v\n1JoY3Lz0dtx0UzGUisEFSFxxp762o2SvUGZou7MN6kxj+GaQ27i1T7ikpASXLl2CIAjIzs726VQT\nQgjxVWZmFjjOgvLy8rAogxuoqVeGYZCdnYPs7Bzceec9aG5uAsBAq1Wiq8uIkydPYNu2zdj40Vv4\naud2LFt2K2JidJBIJGAYBjKZDCkpKVAqHZejdKe+tliyl7MgeK1DDkNcN+QBzKNztRbfZZKEdQa5\njVtvEcdxkMvl4DgO1dXVAEBblAghIZWdnQuTyYzjx7+BVBrYrGl3BXLqVSaTISkpGcCNko7z5i1A\nUVExtm3bgt27d+Ltt98c9DqGYZCYmIj09AxkZo7EzJmz+yV7OZtOtxFL9nIWBI2cBLsv65CiDdzU\nr6u1eDCIiBOjXP7mvvfee/jd734HnU5nT51nGIZOUSKEhFxe3hiYzWacOXMqbAJxoHE80NHLgOO/\nDfpRGtx1172YM2cejh07CouFA8/zEAQBvb09qK2tQXX1ZXuJza+/3ov/+q/HoNcn2hOWRsdbp2av\ntCnACYMDl9jatvMRNAMjF9ipX1dr8VGyyChx6fK39u9//zu2bdtG5woTQsLS+PEFMJtNqKw8H+qm\nBNTgJCN5vySj5OQUJCeniL5WEAQ0NTWitLQEpaUleOmltVjxH7+FUZbcbw17/sgWnGlQo6nH9dq2\nOyNoILBTv87W4iVMZJS4dBmE9Xo9BWBCSFibNGkKOjo6YDZ3hbopAeNLEhTDMEhI0OPOO+9BWtoI\nHKlh0cqkAZbB1ypM6QbHu7e2bQuC1zqsJ06JnW4UyKlfV2vxkVDi0mUQnjFjBl5++WUsW7as35me\ntCZMCAkn06bNwFdflbj9/HDeOzqQPw9EmDZ9Fjoqo2zxt59r7Szy4gVIWcatoGkLgoa4buy+rIOR\nC83Ur6O1+EgocekyCG/atAmAtVKWDa0JE0LCjUKhQHFxMT799AtIpY636vhj72iwA7g/D0TotUhg\ngfhBEj1mCf745zcw96apKCyc5HZ1MrkUSNGG79RvOJe4dBmEqUY0ISRSjBw5Eqmp6airu+awBrMv\n07qhKv7gjwMRbB8cZBLH1+KM7bhQcQIVp8qh1ydi0aKbUVw8vd8sqCORMPUbjtxOJ2xqaoLRaLR/\nnZqaGpAGEUKIL6ZNm47t27eC4wZPuPo6rRuq4hS+HIgg9sFBKuEBDA7Co1JVmP7scygp+RcOHNiH\n999/F+vXv4/s7Bzk5Y1BXl4+srJGihZJiYSp33DkMgjv378fq1evRlNTEyQSCcxmM3Q6Hfbv3x+M\n9hFCiEekUimKi6dhz57dgwp5+DKt62sA93UK29uRptgHB4CFVm6GhZeIjOiTcd999+PWW5djz55d\nOHnyBCorz+P8+Qps3boZKSmp+MEPfoi0tBGi9wvG1G8kree74jIIv/LKK3j77bfx5JNPYuPGjfj4\n449RW1sbjLYRQohXkpNTkJGRgdramn7T0r5M63obwP01hd13pClVqmHpdZ2M5eyDg4WXYFZGK8y8\neDCLidHh1ltX4NZbV6CrqxMVFefwzTdHcPDgAfz2ty/izjvvwcyZs906etFfIqEWtKfc+gwxcuRI\nWCwWMAyDO++8E3v37g10uwghxCeTJ0+FRNI/2NqmdcW4mta1BXAxzgK4bSRqDfyMfQr7dIParX4M\nxEoArVJwawTo6oODmbd+cHB1ragoDSZOnIwHHvgP/Od/PgqZTIb33vsH/va3v6C7O3hrvp6+lxwP\ndJkk4MIzJwuAGyNhWxWapKQklJaWIi0tDW1tbQFvGCGE+EImk6GgoABHjhzut4bp7bSuN+uy/txa\n5A1/JHQNVFg4ERkZv8Kbb/4Vhw+X4ciRcuh0OsTHJyAuLh7Z2TmYPXsuJBKJX6eNPXkvI2nE7DII\nr1q1Cm1tbXj88cfxk5/8BB0dHXjmmWeC0TZCCPFJTs4oXLx4ER0d7fbv+ZJA5GkA9+fWIk/xAnCm\nUQ0TJx51fNk6FBcXj6ee+hm+/PJznDp1Ck1Njbh48QIqK8/j0KEDKD9cjrn/9jTaLDF+C4KevJeh\nPN3JUy6D8Ny5c6HRaFBQUIAvv/wSANDZ2RnwhhFCiD8UFRXjX//6bFBtaW8SiNwJ4H1Hf4EYibpr\nYCCyYRke6TFGn7cOsSyLJUuWYcmSZQCsB/00Nzdh48Z/oldbgOu9cfbn+iMIuvtehnr2wVMum7Jy\n5Uq3vkcIIeEoJkaHnJxc8Lz/Ap4tgPf9Y84LwMl6NXZV6bCzSoddVTqcaVR7vQbtC2eBSCYRMCbB\n/9OyLMtCr0/Egw/9J3ILF4g+p65T7vX6rLvr+e6MmMOJw9ZYLBb09PSA53n09vaip6cHPT09qK+v\nR09PTzDbSAghPiksnASZLLDnoDtKGhIEIEvXA5WUAyBAJeWQpesJaBELZ4GolwtsIDJyLHhWPFHK\nWRB0J4kqX9/t8r30NoEuVBxOR7/xxhv405/+BIZhUFhYaP++RqPBAw88EJTGEUKIP7AsiwkTJuLg\nwX2QyRyXtPSWs5FnfZccc7Nag1rEIpTT4M7u3dPegIP79mLGtGIoldbpaUdJVNO0wqDXu7Mc4Eth\nk1BwGIQfffRRPProo/j1r3+N5557LphtIoQQv8vKysKFC+fR1tbq92u7mzQUrPrFoQxEzu5df7EM\nO0rfxeaNH2Hy5CmYPn0mjDETRJOo5FfMMOjEZ11dredHUglN9le/+tWvnD1hzpw5AID29nbs3bsX\nJpMJCQkJXt+wu1t8Tj8SRUUphlR/xAyHPgLDo5/URyA+Ph6VlefdPpjAXSwjoLZdAQs/+LoqKY/c\nuB6/rsEqFDIYjWLnIN2QoDbDzDMwWiSw8AxUUh4joq0JWYGur+Ho3rNGK6FUKlBfX4+KinM4WFYG\nbfZiSBWDp697zQwyonu9et8YBkiMMiMzphfp0UbkxvUgWWN2u9+CICA/f5znN3YiKkq8/rbDkfBP\nf/pT/OAHP0BeXh5aW1uxfPlyaDQatLS04Mknn8S//du/+bWBhBASaEq1FrGJmWiuvwIp679IFE5T\noH2zswNRy9mdvb+Op41j8J3vLMfSpbdag/DRk5BHxYpeo8vEeL2Fq28bw/X0JBuHQfj06dPIy8sD\nAGzevBk5OTn4+9//juvXr+OHP/whBWFCSMTgeB7rSytxtKIBze1GKKWxfi/eEOopUGcFKvwRiMwc\ncKpejaYe9wtgODznVyJBXt4YjDKMQelFHkZ+8PoxzJ1oa76OqOREt9sYSUU6bBwG4b5HVx0+fBgL\nFy4EACQnJwe1VighhPhqfWklSspr7F8HonhDqE8RClSBCltgu9KmACfc6JA/rs9KgJRoC6paByfL\nXTy+E9t2/x3p6elYsGAxioqmuVxGiKQiHTZOe1RXV4fe3l4cOnQIRUVF9u/3PdKQEELCmdHM4WhF\ng+hjvuxbdURsD3GguSpQ0bePntZTtgW2vgHY2fU9JbbtaISmHbPGKDB+fAFqa6/i7bffxIsv/grH\njn0DQRicNQ149h6EE4cj4YcffhgrVqyATCbD5MmTkZubCwD45ptv6CxhQkjEaOs0orldfODQY2HR\nZRQQrYrs2T13srNVMt7jqVpnga3v9Vt6pIhVWbz64OFoBiE6bxYmFk5BU1Mjtm/fgv379+H1119D\ndnYO7rzzHmRljfT4PQjH9WGHQfiWW27BlClT0NjYaF8bBoCUlBS88MILQWkcIYT4KkajQFy0Ak0i\ngTg+WomZ06bg+LFyAJEbiN3ZF+zNVK2zwNbXgdpon9dfHa0fx8cnYNWqB7Fo0RJs2bIRR48ewcsv\nr8WyZbdiyZJl9sM5Qrk32heM4GhsHyAWCzfooG1CCAmk/910Alv2Xhz0/e/OysZ/rBiPrq4u7Nix\nA11dXSFonX+UV8tQUT94bdWQaEZhmhnbTynRZRocUKPkPJaN7QUA9JgZqGQCbH+iLRwcvs4RQ6IZ\nUzLM3nXCDadOncIbb7yBpqYm5Obm4kc/+hFSUlIAOH8PPGkTx3G47777/NZmZ4IehBsaOoJ5u4DS\n67VDqj9ihkMfgeHRz+HcxxvZ0Y1o6ehFrFaJiYYE3DU/F+y3yT5GoxGffroNghCeI6a+oqNVaG/v\nX8jCWWZwj1mCnVU6iI/2BYzQ9jrMej5ZL34QBCCIXk8l5TA3q9XnNXGxPtp0d3fjww/fw6FDByCX\nyzFnzjyMGTMW2TmjUNmm8zk7muc5fO97d/vWgQH0eq3o9ykI+2A4/1EbaoZDP6mP1iSttk4jYjQK\nKGSDZ+Rqa2vw9dd7B524FG6cBSixfbwcD+yq0olO1bIML5p0laXrwbjE7kHBXcnyiFFaUNclh6Og\nPi+r1ef1V2d9tCkvP4T3338X3d3W6XSpVIrc3FGYOKkYhVNugkYp8erDQDCDcHj/phFCiB8pZCwS\nY8UPFwCAtLQRyMzMwpUr1X6vqhUsYmurzoqJONL32L+BiVOA46AezPXXKVOKUFAwAefPn8eZM6dw\n9uwZ+3+fbt+IRYtuxqxZc/ttuQ03FIQJIaSPqVOL0dDQCJOpN9RN8SuxYiJxKjNqO8QD1MCM4oHB\nPVwqhMnlCowdOw5jx1rLTLa0tKC0tAR79uzEJ598hM8++xQLFizC/PkL7IdGhBOXtaP9bSjVrqVa\nvEPHcOgn9dE9DMMgISEBFy5Uhu1o2J3a0QOJ1VNO0pi9rnkd6NrU3vQRAFQqFfLzx2L27DmQyeS4\ndOkCTp48gb179wAA0tMzXC43BLN2NAVhH9AftaFjOPST+ug+lco6ZV1Xdz0sA7G3AQqw7suVswIk\njPX/u80StPYOzigeEW1EssZxRrGvhyS44ksfAUAul2P06DzMnj0PCoUcFy9W4sSJY/j6673geR46\nXSzUavGlibA4wIEQQoazsWPHoa6uDi0tTUO6VK+jmtej47vRZXJdftPVsYKhplKpsHTprZg7dwF2\n7PgCO3Z8iU2b/olNm/6JzMwsTJo0GZMnT0VCgj4k7aPsaB9QtunQMRz6SX30nNlsxvbtW8Dz4RVk\n3Mkc9pQtq1rO8jjXFPpDEALRRwDo6urE0aNHcORIOc6ePQOe5yGRsN8W/1gKlmUpO5oQQsKBTCbD\njBmzsGtX6ZAvMmQb0Q7cF2yrrGXmGBQkdQU16SoQoqI0mDlzNmbOnI3Ozk4cO3YUW7duxtatm3D8\n+DE88MBDSEx0/+QmX0X420kIIYGVmJiIsWPHwWLxfn3SxtPDE4LNWa3o2g4FdlXpcLJeDT6o86eB\no9FocNNNs/Dcc79GcfF0XL58Cb/5zfPYtas0aG2gIEwIIS6MHTsOen2iwxN8XOEFa+WpXVU67KzS\n2YOZmQuvoOy8VjRjHxWfbnC81zoSqdVqPPDAD/DDHz4CpVKJjz76EPX19UG5N01HE0KIG2bNmoMv\nv/wcPT2er1M6OjzBekYvEzaHzzs7BKGvvoU8hpKJEydj1CgDrl6tDdqU9BB7CwkhJDCkUikWLrwZ\ncrnzo/0GcjbFay0XGT4jTFtlLVdshTyCKVhT+RqNFrm5owJ7kz4oCBNCiJtkMhkWLrwZUungfbWO\nuHscIBAeh8/n67uRpeuBkuVgPaRhsGCWpnQ0lT9U1qUpCBNCiAeUSiUWLlwMlnVvNc82xeuOUIww\nB5Iw1lrR80a2YoRWvHRnMEtT2qbyrVPkg2cNwj3ZzRVaEyaEEA+p1WosXLgYpaUlMJmMTot5eHJ4\nQjgdPs9KgILkbkjZwYU8bAU+As3ZVP71DjkEAajvCu1+Zl9RECaEEC9ERUVhyZKl2L17J1pamsGy\njpOZBlalYhlB9PjAYB9+4IptVNz3BKVgtI/jgY5eBl1mx1P5vZwEl9sGJ7sB1jZHCgrChBDiJZlM\nhgULFqGs7CCqqi45PBhgYDBzVpUqHAWyNGXf848ZZuDZxfJvP7C4P7SNtMxtCsKEEOIDhmFQVDQN\nWm0MTpw45rSyVt9gFooRZjjhhf4BVyXlIZXw6DDdSHrr5TyvUjbwCEZn+n4ACNX7T0GYEEL8YMyY\nMUhISMD+/V/BbDa7dehDuB9+EEhie6cB8aDLMjzkrGAP1vooExq65KL7md1ZVxf7ABCq9WQKwoQQ\n4id6vR5Ll96KAwf24erVWpfn1g5XzhKuRJ8vMJia2gZWAvuo9WQ9RJPd3FlXd1Q8BQj+evIwmwAh\nhJDAkkqlmDlzNqZMmQp+qGxm9TNP9k4D1tFtlNz6ny3A2vYzq6TW/cwqKYcsXY/LdXVnHwBCsU+b\nPqYRQkgAZGfnQiqV4eDB/U4zp4cjd8tj2oiNbr3N3Hb2AcC2nmwN7MFBI2FCCAmQjIxMFBZO9MsJ\nTEOJs/KYWrnZo9GtbV3d3cQqZ8VTQrFPm0bChBASQKNGjUZPTy/Onj3ttzViy7enL0VyVvXAvdN9\nk6MEAZAq1bD0+n+rkbPiKbYRNx/EOExBmBBCAqygYAJ6e3tx+fIln6ambVm9Dd0KdJlUEVslCnAx\nncwAWqWAdtdnSXjF2QeAYKMgTAghQVBUVIyOjna0tbV6fY1wyur1l1Bs0wpVJTDRtoTmtoQQMvxM\nn34TBMG7jOlwy+odCjxdTw4ECsKEEBIkarUaEydOgtls9vi17mT1kshDPzVCCAmi7OxcJCenePy6\ncMvqJf5BQZgQQoJs2rQZADzLpHK2rSfcTl8i7qMfGyGEBJlCocDkyVM93j9sqxJlTWRyv0oUCV+U\nHU0IISGQmZmJrq5OnDp1Eqybw1hbVq86SkBDizGi9wkTK/rxEUJIiOTnj8XixTcjKioKHOf+qFjK\nhj6r11ccby04MtyzumkkTAghIRQTo8Pixbfg5MkTOHv29JCvMx1OxwiGgwj+HEUIIUMDwzAYP74A\nCxYsgkwm83ovcSSwFRyxHt7A2AuOnG5QB/ze4Tj6ZoQg/7QtFg5S6dD+pEcIId6yWCzYvXs3rl69\nOuTOI7ZwwPZTSnSZBo//ouQ8lo3tRSDCAy8AR67IUNvKosvEIEouIE3HYVK6WXT0zXEc7rvvPv83\nRETQf8ItLUMni0+v16KhoSPUzQio4dBHYHj0k/oYOQoKiiCTncbJkydEk7aio1Vob+8JQct802WS\noMs0+OAE62MMGlqM9hKW/uzjyXo1qlpl/e5VUS+ByWQRLffJ85zff4/0eq3o92k6mhBCwtCYMfmY\nNWtOqJvhV6EoOBLu5T4pCBNCSJhKTk7GvHkLwDBD4091KAqOhHu5z6HxkyWEkCEqNjYOixcvgUwm\nc/3kCGArOKKSchhYcCQQiVPhXu5zaK36E0LIEBQVFYXFi2/Bzp070N7eFrD7CIIAnufAsoELDWLH\nCDLM4G1L6XE8cmN6fN62ZBt99z0C0iYcyn1SECaEkAigUCiwZMlSXL1ai8bGWjQ3XwDLsmAY3zfX\ncpwFMpkcGRmZUKujcPr0KQgC75drO9L3HGFr4lT/c5Ir6lmYTGq/nJNsK+sptjc51CgIE0JIBElN\nTcOECXnIyyvE6dMncf78ebfLXvZlsVjA8zxiY+OQm5uL7Oxce9DNycnF0aOHcenSJchkgQ0TrhKn\nxiR0ezRa5XjYR9i214mNvkM9ArahIEwIIRFILpejsHASUlLSsG/fXtECH2azCTKZHFKpFCwrhVTK\nQqVSIyZGh4SEBCQmJomuNctkMhQVTUNu7iicOnUSHR2d6Opqh9lsgVwuh0TivwjmTuKUbcTsjDuV\nuPqOvsMFBWFCCIlgSUlJWLJkGXbv3onOzg4wDAOO45CWlob8/HGIjY3z+tpxcfH2bVKCIKCrqxOX\nL19GRcVZcBznl2BsS5yyVtDqz5PEKVslLhtbJS4AfpnSDhQKwoQQEuFUKhUWL16C8vJDYFkWY8eO\nh1Kp9Os9GIaBRqPF2LHjMHp0Hr755gguXrwIqZQFz/PgeR56fSLUajUuX65yu9qXPxKn/D2lHUwU\nhAkhZAiQSCQoKpoWlHtJpVJMmVKEvLx8nD59EvHx8cjMHGkPvAzDoLr68qCRsth6LSCeOGXNjnZv\nBOuvKe1QoCBMCCHEKxqNRjTwT51ajObmZnR3dwFwvV4rljgVq1Ohvd29dvhrSjsUwnSATgghJFIx\nDIPZs+cCsGZEuXtyki1xytOp41BU4vKXMG4aIYQQG6OZQ31LN4xmLtRNcYtarcaMGTfBaOKCUrvZ\nWSWucEbT0YQQEsY4nsf60kocrWhAc7sRcdEK3DQhDbdOzwDrx61CgZCcnIJR+RPwZdU10cf9uV4b\nznuBnYmAJhJCyPC1vrQSJeU1aGo3QgDQ1G7Elr0Xsb60MtRNc8v4MQbERytEHwvEeq23U9qhEiHN\nJISQ4cdo5nC0okH0saMVjRExNa2QsZho0Is+lhptiZhgGSjDvPuEEBK+2jqNaG43ij7W0tGLtk7x\nx8LNXfNzsXDKCMRHKyFhgPhoJRZOGYEn7pvr9n7ioWp4954QQsJYjEaBuGgFmkQCcaxWiRiN+DRv\nuGElEty70IA75uSgrdOIGI0CCpl1O9HChTejuroKJpMJZrMZZrMZra2tMJl6Q9zq4KAgTAghYco2\nlVtSXjPosYmGBHsgixQKGYvE2P7bktRqNfLy8gc9V62WYP/+I6iru462tlawbGT11V0UhAkhJIzd\nNT8XgHUNuKWjF7FaJW6akIpbp2eEuGWOGc3coBGvp6KiojBhQiEA4PTpUzh58viQnLoeej0ihJAh\nRGwqd0SqDg0NHaFu2iBi26kmGvS4a36uT9up8vPHQhAEnDp1YsgF4qHVG0IIGaLEpnLDjW07lU1T\nu9H+9b0LDT5de+zYcRAEAadPnxxSgZiyowkhhPgsGNupxo0bj/z8cbBYzD5fK1xQECaEEOKzYG2n\nGjduPGbOnAONRgOLxeKXa4bS0BnTE0IICZlgbqdKTU1Damoaqqsv4+TJE+js7IjYKWoaCRNCCPGZ\ns8pYgdpOlZGRiaVLv4OioumIiYmByRR509SR+dGBEEJI2BHbTjXRkGD/fqBkZWUhKysLbW2tOHPm\nNGpra2A0msCyEshksoDe21cUhAkhhPiFs8pYwRATo8O0aTMAAEajEa2trWhpaUZbWxtqai4DYMAw\nTNDa4w4KwoQQQvwqHLZTKRQKJCUlISkpCQAwceIkHDt2FFVVl8Kq+hatCRNCCBny5HI5pk4txq23\nrkB6egakUhnM5tCvIdNImBBCyLChVCoxdWoxAKCurg6XL1/C1atXwXGWkExVUxAmhBAyLNmmqzmO\nw44dX6K9vQ0SH8preoOmowkhhAxrLMti4cLFiI2NA8/zQb03BWFCCCHDnkQiwbx5C5CQoA9qIKYg\nTAghhMAaiOfMmYcRI4J3TCQjCIIQtLsBsFg4SKXhkx5OCCGEhErQE7NaWrqDfcuA0eu1YXmmpz8N\nhz4Cw6Of1MehYzj0c6j1Ua/Xin6fpqMJIYSQEKEgTAghhIQIBWFCCCEkRCgIE0IIISFCQZgQQggJ\nEQrChBBCSIhQECaEEEJChIIwIYQQEiIUhAkhhJAQoSBMCCGEhAgFYUIIISREgn6AAyFt4YRkAAAG\n3UlEQVSEEEKsaCRMCCGEhAgFYUIIISREKAgTQgghIUJBmBBCCAkRCsKEEEJIiFAQJoQQQkKEgjAh\nhBASIhSECSGEkBChIEwIIYSECAVhQgghJEQoCBMi4rPPPsOKFSuwfPlyLFmyBD/5yU/sj7322msw\nmUx+v+f8+fNRUVHh83Vee+01vPTSSwCADz74AG+//bbP17RZuXIldu7c6bfrETLcSUPdAELCTX19\nPZ5//nls3LgRKSkpEAQBZ8+etT/+pz/9CQ8++CDkcnkIW2llsVgglTr+Z3zPPfcEsTU3cBwHlmVD\ncm9CIgkFYUIGaGxshFQqhU6nAwAwDIMxY8YAAJ5//nkAwN133w2JRIJ3330Xu3fvxj/+8Q+YzWYA\nwNNPP43p06cDsI5uly9fjn379qGhoQEPPvgg7rvvPgBAeXm5/XpTp05F37NUXnrpJRw6dAhmsxmx\nsbFYu3Yt0tLSUFNTgzvuuAO33347Dhw4gDvvvBPf+c538Oyzz6KiogJ6vR7JyclISEgAYB0Vd3d3\n4+mnn8Zzzz2HY8eOAQC6u7vR2tqKsrIydHZ2Yt26dTh37hyMRiOKi4vxzDPPgGVZVFZW4plnnkF3\ndzcMBgOMRqPoe3bw4EG8+OKLGDduHE6fPo0nnngCnZ2dPr8vxcXF2LFjB/7yl7/AYDDg4sWLWLt2\nLVpaWmA2m3H//ffjjjvu8PlnTkjICISQfjiOE370ox8JRUVFwo9//GPhrbfeEpqbm+2PGwwGobOz\n0/51c3OzwPO8IAiCcOHCBWHWrFn2x+bNmyf89re/FQRBEK5cuSIUFhYKnZ2dgtFoFGbOnCkcOHBA\nEARB2L59u2AwGIRz584JgiAITU1N9mt89NFHwhNPPGG/hsFgELZv325/fN26dcLq1avtr5szZ479\nnn/84x/t/29jMpmElStXCm+//bYgCIKwZs0aYePGjfa+P/nkk8L69esFQRCE2267TdiwYYMgCIJw\n9OhRIS8vTygtLR30nh04cEDIy8sTjhw54vP7MmvWLKGsrEwQBEH44osv7O+L2WwWbrvtNqGyslIQ\nBEHo6OgQFi9ebP+akEhEI2FCBpBIJPjzn/+MiooKlJWVoaSkBG+++Sa2bt1qHx33deXKFfzkJz9B\nXV0dpFIpGhsb0dDQAL1eDwBYunQpAGDEiBGIjo7G9evXYTaboVKpUFxcbH/Oc889Z7/mnj178P77\n76O7uxsWi6Xf/RQKBW655Rb71wcPHsTPf/5zAEBcXBwWLVrktH8///nPYTAYcP/99wMASktLcfz4\ncbz11lsAgN7eXiQlJaGzsxMVFRVYvnw5AKCwsBAGg8HhdTMzMzFx4kSf3xelUokpU6YAABYtWoTo\n6GgAQFVVFS5cuICnnnrKfg+z2YyLFy8iJyfHaZ8JCVcUhAlxwGAwwGAw4Pvf/z6WLl2KQ4cOYfHi\nxYOe99RTT2H16tVYuHAheJ7HhAkT+k3bKhQK+/+zLAuO45zet7a2FuvWrcMnn3yC9PR0HDlyBD/9\n6U/tj6tUKjAMY/9a8OBI8P/5n/9BR0cH1q1b1+/1f/7zn5Gent7vuZ2dnf3u44pare73tbfvi6N7\nCoKA2NhYbN682e02ERLuKDuakAHq6upw9OhR+9fXr19Hc3MzRowYAQCIiopCZ2en/fGOjg77Y598\n8olbmdPZ2dno7e1FWVkZAODzzz9HR0cHAGvwk8lk0Ov14HkeH374odNrTZ8+HRs2bAAAtLS0oKSk\nRPR5GzZswFdffYVXX30VEsmNf/rz58/HX//6V3sQbG5uxpUrV6DRaDBq1Chs3boVAHD8+HGPsre9\nfV+6u7tx+PBhAEBJSQna29sBACNHjoRSqcSmTZvsz79w4UK/nwUhkYZGwoQMYLFY8Nprr6G2thZK\npRI8z+OJJ55Afn4+AODBBx/EqlWroFQq8e677+KZZ57BI488gpiYGMyaNUt0ynoguVyO3//+9/0S\ns1JTUwEAo0ePxpIlS7B06VLExsZizpw5KC8vd3itRx55BGvWrMGSJUug1+vtU7kD/elPfwJgTSoD\nrB8m3n//faxZswavvPIKli9fDoZhIJPJsGbNGqSnp+Pll1/GM888g//93/+FwWDA+PHj3X4fvX1f\nXn31VfzqV78CABQVFSE+Ph5arRZSqRRvvPEG1q5dizfffBM8zyM+Ph5/+MMf3G4TIeGGETyZyyKE\nkADr7OyERqMBABw4cACrV69GaWlpv9E7IUMFjYQJIWHliy++wNtvvw1BEOwjYwrAZKiikTAhhBAS\nIvTxkhBCCAkRCsKEEEJIiFAQJoQQQkKEgjAhhBASIhSECSGEkBChIEwIIYSEyP8Hxa9vs0xocdYA\nAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "\n", "ax.scatter(df.std_range, df.std_logratio,\n", " c=blue, zorder=10,\n", " label=None);\n", "\n", "low, high = np.percentile(pp_trace['obs'], [2.5, 97.5], axis=0)\n", "ax.fill_between(lidar_pp_x, low, high,\n", " color='k', alpha=0.35, zorder=5,\n", " label='95% posterior credible interval');\n", "\n", "ax.plot(lidar_pp_x, pp_trace['obs'].mean(axis=0),\n", " c='k', zorder=6,\n", " label='Posterior expected value');\n", "\n", "ax.set_xticklabels([]);\n", "ax.set_xlabel('Standardized range');\n", "\n", "ax.set_yticklabels([]);\n", "ax.set_ylabel('Standardized log ratio');\n", "\n", "ax.legend(loc=1);\n", "ax.set_title('LIDAR Data');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model has fit the linear components of the data well, and also accomodated its heteroskedasticity. This flexibility, along with the ability to modularly specify the conditional mixture weights and conditional component densities, makes dependent density regression an extremely useful nonparametric Bayesian model.\n", "\n", "To learn more about depdendent density regression and related models, consult [_Bayesian Data Analysis_](http://www.stat.columbia.edu/~gelman/book/), [_Bayesian Nonparametric Data Analysis_](http://www.springer.com/us/book/9783319189673), or [_Bayesian Nonparametrics_](https://www.google.com/webhp?sourceid=chrome-instant&ion=1&espv=2&ie=UTF-8#q=bayesian+nonparametrics+book).\n", "\n", "This example first appeared [here](http://austinrochford.com/posts/2017-01-18-ddp-pymc3.html)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 0 }