{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Dependent density regression\n", "Author: [Austin Rochford](https://github.com/AustinRochford/)\n", "\n", "In another [example](dp_mix.ipynb), 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": true }, "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": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/fonnescj/anaconda3/envs/dev/lib/python3.6/site-packages/pandas/io/parsers.py:2108: FutureWarning: split() requires a non-empty pattern match.\n", " yield pat.split(line.strip())\n", "/Users/fonnescj/anaconda3/envs/dev/lib/python3.6/site-packages/pandas/io/parsers.py:2110: 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": {}, "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", " \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": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFpCAYAAACrqZC7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcU+W9P/BPkpksszEr60CFYbEFB0coq5QyTFFRvFQQ\nfhfBy7W1tnXrtfVX9WWV1qoItHVrX7Wtrbf6w9KXWG29KoqgRZHiAox4gYFBZXe2zEYmyyTn98eY\nkMmcLcnJOSfJ5/2XTDLJ8yTj+Z7neb7P97EIgiCAiIiIdGc1ugFERETZikGYiIjIIAzCREREBmEQ\nJiIiMgiDMBERkUEYhImIiAySo/cbNjd36f2WmikpyYPb7TG6GbrIlr5mSz8B9jVTZUtf072fFRWF\noj/nSDgOOTk2o5ugm2zpa7b0E2BfM1W29DVT+8kgTEREZBAGYSIiIoMwCBMRERmEQZiIiMggDMJE\nREQGYRAmIiIyCIMwERGRQRiEiYiIDMIgTEREZBAGYSIJvkAQTW4PfIGg0U0hogyle+1oIrMLhkLY\ntO0I9jQ0o63Th9IiB2rGV2B57VjYrLxvJSLtMAgTxdi07Qi2vn8i8u/WTl/k3yvqxhvVLCLKQLyt\np7Shx/SwLxDEnoZm0cf2NLRwapqINMWRMJmentPDHd0+tHX6RB9zd3nR0e3D4JI8Td+TiLIXR8Jk\neuHp4dZOHwScmx7etO2I5u81qMCB0iKH6GMlhU4MKhB/jIgoEQzCZGpaTQ+rncp25NpQM75C9LGa\n8eVw5GbmmaZEZAxOR5OpJTs9nMhU9vLasQD6gry7y4uSQidqxpdHfk5EpBUGYTK18PRwq0ggVjM9\nnEims81qxYq68Vgytwod3T4MKnBwBExEKcHpaDK1ZKaHk53KduTaMLgkjwGYiFKGI+Es4QsE03ZU\nl+j0sJqp7Mo422Kmz9FMbSGixDAIZzgttvcYfbFPdHo42ansaFpuk0r282RFL6LMwSCcAkYHrWjJ\nVH9K9GKfiv4n8prhqezo/ofFm+msRRUtqc9z8ZzR6PYEVPeNFb2IMgeDsIb0GKHEE4yU1kSXzK2S\nfY14L/ap6H+yr6lFpnOyn2OY1Of5dv1p+PxBVX3Tqi1EZA4MwhpK5QglkWCUzPYer7837ot9Kvqf\n7Gtqkeks9zm2daqroiUXPL3+vgQxNX1jRS+izKJ7EC4pyUNOTvreqVdUFIr+3OvvRX1jq+hj9Y2t\nuGGJC0574h/371/4SDQY5bnsuH7xBaK/UzjIhYoSF5rcPQMeKy92oeq8Msk2nW45i7Yu6Yu9zZ6L\nivL8yM9S0X+tX1MsCcvr70WvxYqSIofka8l9jhYr8M+PzuA7iy+AzSY9Mpf7PGPJ9S2Z7xSQ/vvN\nROxr5snEfuoehN1uj95vqZmKikI0N3eJPnaiuRvNIhdGAGhp70Hjp60Jj1A8vgBe+9dnoo/9c88J\nzK8ZjsI8+4DHfIEgxo0YJHrBrq4qQ1dHD8R7A5QMcqG0UDqpKegP9PssmtwezfufitcMC88s1De2\notndozizUF1VJrq2HAoBL+/8FD09fqy65PzIz2OXDYKBoOTnGW/fpNqi9J3K/f1mGvY186R7P6Vu\nILJiOlpqHVWLBKLwxfzDQ00QJJ4jlokbz3tvfP1wZMoyVnu3H2v++B6mnH8ugERPXbd2+uC0WwFY\n4A8EVa+JOu05cSU1aZmJnMrXDIt3mnt57VgEgyG8tfcUQiJf9Ft7TwEWC5bXVuG5N4+KLhtIfZ7x\n9o0VvYgyR0YHYal11KVfHyN5oQyPgtQGydiLuZjqqtLIa+XYLIpru9HvDQAHP2uTfX13d/8AEtsm\nrz8EAJg1aShWXTJB9Q3H8tqxCAkCdn50JnIT4LTbIAgCgqFQvxGjlpnIqXxNAOjy+PHBwfjWu21W\nKy6ZNgpv7jkl+nshAdj+4UkcOdGB403dkZ9HB/fY4GnPtYneXCn1jRW9iDJHRgdhqdHOoWPtshdK\nqSAZSy7ZJsxmBd79+Aze3HMKpUUO5Dlz43rvCaNK4O7yq+rvnoYWLJp1nmSbDh1rV/U6QN9aaWuH\nF4KAfoHC6w/ijQ9OwmKxDBgxpmKEpuVrhm/K3j/YhPZu8c9ULrlJbmQedrK5W/Tn4eAeHTwL8ux4\nYcfRhPsWruhFROkrY4OwXICUu1AGQwK2f3gy8rNwkAwGQ7ht5Vf7PV8uUzUsGAKCX4xEWzt9khfw\nDw81i773zv1n4LSLj5hiubu8ONHUrcmBB/WNrWhy98BqEX+e2Igx3hGamtkGLUd9amYt5KaC5Ubm\nYWJT1UD/zz46eKrpm5n2nRORtjIuCIcvWP7ekGQwkrpQtnZ6JQP3W3tPweXah29efF5kGlbNyEit\nti4fdtafTuo1SgqdqBxcoOmBB2qCSiylEVoi263EXlOrPdPRlKaCldaGrRbxz6yk0AmXIwdNbs+A\n9kp9XqyMRZT5MiYIx16wSgrtcKgcQUaTmqYMCX1ZsH5/L1bUjY8EgOqx5f1Gr8nw9YbEf+4PYvak\noTh4rF1xLbEwz57wOqraQAUARfl2uByJ/flILRN4vL2q1qzjCU5qbsrCpNa7o9ms1r4saItF9Hsf\nUVHQb7khLM+Zg5899V5cwZSVsYgyX8YE4dgLVpvKddRYFgsgSKU549y0cf2RlkiwHzm4AB5vAG2d\nPlgkRkLJKC1yYuUlEwBA1Vri8tqxCIYE7G1oQftZH0q/eHzxnNGiI7F4AlVYe7cfP3vqvYTqUEsF\n+p37z+DQMbfia6oJTonclMmtd8daUTcONqtlwHdwLunv3M/znDmSeQBS78PKWETZIe2DsC8QRLPb\nI3nBctptyHfmoK3TJ7mFKJpcAAb6po2jR0BtXX60dfkxr2Y4Lpk2ClveO644Mh5WmofTber3S0eP\nYJXWEiNrukda4O72objAjkljSiAIAu59crdslrhcoBKbZk1kZKa0jq70mmqDUzI3ZWqCnNhaNQC0\nd/mwZG5V5OcuR98ION73YWUsouyQtkE4di+sFK8/iB9fUwObxYJHnqvXZP1WTH1jG5bVjus3Qmrr\n9MJh77vARu/RXTxnNO59crdoW8I3De4uH0oKnaiuKsW8mhHwBYIDLtZia4mxwae924+39vZfa5bK\nEpcLVBdXD8VHjW64uwe2OZ6Rmdp19HCmd4+vt99NhprgNKjAkdRNWTxBzpFrQ9kgp+RWuD+/2iDZ\n10QzsaXW9pnARZR+0jYIq8l0DfvnvtNYtWCCpuu3saIvqGIjpNiLo9S67cXVw7BkbhXaOr3Y+sEJ\n1B9pwfY9p1BcYEfNuHKs+Mb4hIr7i5HKEo/mtNsw+4KhqL2oEjv2nRF9TrxBS03RitZOL9b88T20\nd/dfQ1UTnOQCtT8QxF0rLwIAyZuyeAuBqN0KF8/7yH1OE0YV9/t3IglcDNhE5pCWQTjeYPPhoSaE\nQgI+amwBcG5q1Wm3IiQA/oB4QhQAFOfb0X5WeSoz9oIaO0qNDVBy+19tViu27znZ74ahvduP7XtO\n4cjJTtyzeqroxVXNlqloatauvf4gLBYLBhVIT1XHG7TCff/wULNsPeXwqDt2iloqOI2rLAKgPIqs\n+GKbkBaFQDy+AN6WyGpXusmJfp9wUHQ5ciKj/+i/kehZlXdj1s7jSeBixjWRudjWrFmzRs839HgS\nS5iK1tbpxUs7xWspi/EFQvjsTBd6vggg4dgzp3oYbllSjX99/LlocHHabbjl6mrJi2y02RcMRc24\nCtVtslosuGBMGeZeOBwXXzAMC2d+CTXjKmC1WOALBLHx9Qb0+Aa2qfOsH909AUyuKh/wWE6OFe9+\nfEb098TbAFXr5B3dfrR1+XDkRIfo4wn3vWYE2rt8sqPF2HbMvXA4qqvK0OPrRXuXL/KdAsCJ5rN4\n44MT6O7xY2hZHj45PbDObHRbv3JeCSxWK1rae+Dz96K0yInZFwzF8tqxsFokNkjH+O9XD4m+DyD/\n2c6eNBQrF4yHIAj4yxuH8f9eO4R/7PwMr+0+htffP4FdH59Ba6cPy2vH4us1I+Du8uGT013oDfa9\nao8viKOnOtHt8WPfkRbR77yj24+ZE4ego9uHnBwrigqdePIfH2Pr+ycizw+/To+vF+NHFqOt04uc\nHCtyZA6jSAf5+Q5NrjXpIFv6mu79zM8XH6ik5UhYq/254XXcKeeLj4i8/iDe3X9G9r3KZCpqqRE7\nYvYFgjh6skO2b3sbWrBs3th+o7VgKITNbzXirDeg+r0tFqiKwm2dXuxtaBF9zGm3YfGcMarfM5oj\n14bVC8+Hy5kTmREoyrerqma1om78gOImwLkM59opI1A3tXLATENshvj1iy/AZdNGJjQ16wsEZUuK\nSu0ZLi10YOUlE2CzWrFxa4Povuzo0eySuVU4eMwt+h57DregQ+Lzip3Snz5pmOQM0tv1pzk6JjJA\nWgZhteuKSsIX9cVzRuPt+tOio+H6xjbJteTYWsxq19nEnhebaCYXH9vP+gaswUqtkTvtNsy6YCgs\nAPYebu23zzgoPQvfz6ACO9pFErKAvjXWto4edOfaVBfNiO57bJZxOJtYaa3WFwhi32HpJYl9h1vw\n8+tnDCgRGZshftOymoTLP3Z0+2RLig4pEc+Cv2hCBRy5NlXLKnsaWvC1ycMllxk6uv0oLnCIJswB\n/af0X975qeT7eP3BuM41JiJtpGUQBvrWFT3eXuzcL54spEZ0Io9PYv+ou8uLuil9J9HWN7airdMb\n2XcrdmqR3EhC7nmxQVRugFoaswYrdzHPc+Tg6q/3jZqXfj2IM21nsfaZD0Wf68i1wieyPl4zrhz1\nja2igdGea8Mjz9XL9tsXCPZLNBN7bnQgVLNW29Htk83mbu304ejJDowZMQiDS/IGjDijz2NePPs8\nydeRIzcj47TbcMeqKfjHO59I7udWs4bv7vICgiD5PqVFTlSPLVOdcCg1OhfD/chEqZe2QbivctEE\nHDrmFr04KRXdAM5d1OUupsUFjkjwaO3s23dbXVWK5bVj0RsU0NrhGbA3WGokIZVAEwyGJA+vl2t3\nmNzFvL373KjZkWvD6++dEA20QN/a+YyvDMbhE50Dk8Vs4iNtuRGU3DYyudGWmkMbBhU4UFpolwzE\nVguw4S97UVrkQHVVmeTnu2v/aVw2bWRkZBqbHCUXgORmZC6uHoZCV65sbWg1yyrhRDKp98lz5mB5\nbVW/wiFyU/rxFJLhfmSi1EvbIAzIXwQrJcoHAkBZUf+Lutzr5LtyJbOUPd4AWjt9qg45kC0yIbOu\nBwBF+bno8gT6jcCjqd1TqrSGCQCHT3SguqoMyxacD/Se25scGxiLCxzw+HpFp/DD/d78VqPikoHU\nQRBL5lbha9XDAIsFFcUu0T3SF00YLPn60Wur2yWOHwSAlvYetHV6sX3PSXx4qAltXf7IaLFMxdqo\nmhsGqeluNcsq4Ruu5bVjRbc8HW/qxnNvHlU9pV9a6MDkceWoP9I3q1OUnwuvPyh6Y5bsmc1EpCyt\ngzAgfREUKx9YXVWKuqkjUVrkHHBRF3ud6qpSyRFU9MVQzSEHcqNVuXW9siIn7lk9VXZklmOzIM+Z\nK3rRjZ3CVToWMRy0Cguc/aZpY9du/b0h3Pvkbsl+y1Uxi9bW6Y1MGztybaqm9sMj1sVzxgw471iK\n1DRsebELWz840e9GSyw5SmptNNFTnqL7AJzbriV2AwAAvUEBHomku+gbGaUp/YsmnHvNvQ0taO/2\nwWEXv8GIPgebU9JEqZH2QVjuIhjPxVHsdTq6fZKHuKsRPZJwOXIkA63cul74UIbCPLvk+2zadkR0\n1D9ycMHAKVyVWeXR07TRwhd6XyAoO/qGxaJqz7Ilatq4ZnwFBEHAGx+IT+1LnfW84cZZaOvw4vl/\nHsXeI+I3TVI3SlO/PAT/2i+/BU3N2qja5C6pm4yffXsauj0ByanweMtYit1Uzp48HItmjsKmbUf6\n/a15vzhq02m3wR8IorjAgXxXLuobWyPnYDNbmig10j4Ih8lN+cWzphX9/GS3QtWML0eOzYKNWxuw\np6FZMoP1XJLXwAMBlLY+yU1ze7y96A0KCG/5jCervKW9R3Y9UKnYRUWxS9VnFzvqdNrFA53cWc9A\n3zaeY5+L79cFgNJCOyaPq0D9kdZ+n++iOWPwikzWMNA3Wm92e1DxxYyGlmcai422xW64lPIW/L2h\nfqVNxW4qK4cX4+hnrfjgoHQS312rpgwoFMNsaaLUyZggnArxboUKF7+IXruVK68ZvTad6LRmc3tP\nwiOktk6v5KlP5cUuxfVA+aWA+PYsh0lNK8vtVQ5v45Gbav/yl0qxasEE+Ob13yJVOEj5ZkEA8MAz\nHwCwwOcPJjwyTOZkJLm/RY+vF/c+uVu0XeGbymAohN+/8BH+ueeEZNJWe7cPEPpOCEukjUQUPwZh\nBWKBJvZourC5F/adpBS+wMtddIsL7Lhn9dQBo554pzU/PNQkuZ1JLLEmNthLnfo0Y9IwxYut1I1D\n7HYgLcjtVVbaxuO02/Dv3+gbwcV+vk57jqobrfCULZD4yDDZk5Fi/xbD+73V7O9VU2tdaRmB2dJE\n2mMQVhAdaGz2XAT9AeTYLF+s64nXfQ6Tu+h2nvWjx9eLwjx7QsX01VxU5Wogh4OR1Lm41y2aiLa2\ns6raEh3Y4q3rHcspUZ9abq+y0jaei6uHIc8h/aceW8taZSGxuEeGiZyMFC36b7HZ7cEjz9WLflbv\nH2zColnnRW7w1H4nSssIzJYm0h6DsEqOXBsqyvPR3Ny37qhm6ljpoluQZ4+sF8dTLlDpohpPKU2p\n0awtwdrB8R4iEWvWBUNhtYivjUvtVY7exgPIbxcSE/sZnO0J4Od//kAxEMc7Mozn0Ai5GzNHrg32\nXJvM3nA/1vzxPUw5v+9vQOk7KSlw4MJxZZhXM+KLtiR/sAURqcMgnASlqWOli+4LO46qPv0mmtxF\n1QLg1qXVqBxcqLIX59qayDRjbLCQywKPNXJwATzeXtHZBLEbHKUgm+i6euxn4CuQzvyOlsjIUKkP\naquvKSUNurv7J61JPXdQfi6qx5b1y4SePK4c86eMiJQ5VXszQ0TxYxBOMamL7uI5oyX32SpNc8pd\ngEuL+qZmUy02WJQU2pHvssPjDUgG4PAWmOiLem9QEA2YYjcFaoNsojcU0b+vZp04kZGhUh/UHkuo\nto3hvyWp5xblO/DW3nPb8Fo7fdj2wUnUTa3Ez6+fzn3CRCnGIJxiUhfdJrcn4QQYrc7CTUZssGjr\n8kuWkAxngS+eMxrdnkC/i7rNOvCsZSXJBlk1pM7yjb2JSJRYH+LNng6///sHmxRPnlpeOxZ5Ljve\n2XdKVTGa8PsxCYsotRiEdRJ70U02SSfR9U8txJN8FZsFnufITWXTNCN28wQgoapYap8fb/Z0uI2L\nZp2HNX98T3QGIvy3ZLNaBxzbKFeMRulGMJFkQiIaiEHYIMmOZpNd/0xGPMlX0Vng6Sj25imZqlhK\nCXeJ3pgV5tklz8SO/VtSW4wmtuZ4+G/s3M4Anj1MpAUGYQNpMZrVY2o2VjyVxLJxW4vadd1YydyY\nJfK3pPR+0dXewgE3z5nbb488q2kRJYdB2EBGjmaTEU8lsWzb1pJMVSxA/LSq879UgsVzRsu+b6J/\nS3LBW+xmQurGi9W0iBLDIGwCRoxmkyUWLPJduX3Z0V2+rN3WkmxVrHAwXTxnDJ59vQEHj7nx7v4z\nOHTMrWraN96/JangHW/RFVbTIkoMgzAlRO7inU6jeq0lm3AX9sKOo3hn/5nIv1M97RsbvOMtupKN\nyw5EWmAmBSUlfPEOB9zYf2eb8FS9GLVT80pT2r6A/NnJWgjfTKiVbcsORFphECbS2PLasaibWomy\nIieslr590nVTK1VPzauZ0k41uZsJoK/wSiJ9I6L+OB1NpLFkE+60mtJO1uI5o/F2/WnRQyLCZw9X\nFLs4AiZKAkfCZCq+QBBNbo8uU66plujUvBZT2lro9gTgkzjfub3bB3uOlQGYKEkcCZMpJFrgIlMZ\nWREtzCwjcqJMxiBMppBogYtMZYY95GaoUU6U6bJviEGmY4ZsYLMyOts82SQzIpLHkTAZLtkCF5Q6\nZhiRE2UyiyAIgp5v2NsbRE4O/yemc7z+Xty4bhua3D0DHhtc4sKv/28tnHbeLxJR5tH9yuZ2e/R+\nS81UVBSiubnL6GboQu++VleVia49VleVoaujB6lqCb/TzMS+Zp5072dFRaHozzm8IFMwQzZwpsn2\nEqJE6YBBmEyBa4/akdruddOyGqObRkQxGITJVNLxRCmzkdruleeyY/Hs84xrGBENwC1KRBlEbrvX\nrv2nVW33yqSqZURmx5EwUQaR2+7V0t4ju92LVcuI9Mf/s4gyiNwRhMWFDrgc0vfd4Wns1k4fBJyb\nxt607Yjq9+comig+HAkTZRC5UpNtnT787Kn3REe3SlXLlsytkk2Uix5Ft3b6UFxgR824cqz4xniO\noolk8P8OogwTXWoyltToNp4zjMVGu9GjaABo7/Zj+55T+NlT7yMYCmnRLaKMxJEwUYYJb/daNOs8\nrPnje3B3DwyusaNbNScmSa0ZL54zWnIUfbypGxu3HsaqBRO07SRRhuBImChD9fh60S4SgIGBo1s1\nZxhLrRlvfP2waPAO25vlh3AQyWEQJspQcklaYucBy52YJLdmfPAzN4oL7JLtaD/r6xfwiegcTkcT\nZah4zwOWq1rW2uGRXDNu7/Zh2peHYNf/fi76eKlIwCeiPgzCRBkskZrcYlXLlNaMV14yASdbzuJ4\nU/eAx8UCPhH1YRAmymDRo1ubPRdBfyChgKg0qs5z5OCe1VOxceth7G1oQftZH0p5CAeRIgZhoizg\nyLWhojw/qaPglEbVNqsVqxZMwLJ5Y3kIB5FKDMJEpIrak654CAeRegzCRBQXBlki7XCLEhERkUEY\nhImIiAzCIExEuuAJS0QDcU2YiFKK5xQTSWMQJqKUCtecDgvXnAaAFXXjjWoWkSnwNpSIUkbpnGJO\nTVO2YxAmopSJ55xiomzEIExEKRPvSU5E2YZBmIhSRs05xUTZjIlZRJRSiZzkRJQtGISJKKXU1pwm\nykYMwkSkC9acJhqIa8JEREQGYRAmIiIyCIMwESUtnrrQ4ed2efysJU1Zj2vCRJSweOpCh5/74aEm\ntHX5YbUAIQEoYy1pymIMwkSUsHjqQsc+NyQo/w5RpuNtJxElJJ660HLPlfodomzAIExECYmnLrTc\nc6V+hygbMAgTUdyCoRC27D4Gi0X88di60HI1pKV+hygbMAgTUdw2bTuC7XtORdZ1Y4XrQoczoft+\nJl5DOvZ3iLIJE7OIKC5y67tWCzC3ZgSWfn0MNm5t6Jc1PXlcOeZPGYE9DS1o6/KJZkcTZRvVQbi7\nuxuhUAhFRUWpbA8RmZzc+q4gAJd8dSSee/PogKzpbR+cRN3UStz/nRno6PbB5chBj6+XtaQpqykG\n4WPHjuG2227DsWPHIAgCRowYgV/96lcYPXp0Qm9YUpKHnJz0/R+uoqLQ6CboJlv6mi39BLTpa+Eg\nFypKXGhy9wx8/RIXKocXo/6v+0R/t76xFTcsmYzK4cVJt0MJv9fMk4n9VAzC9957L7797W/j0ksv\nBQC8/PLLuOeee/D0008n9IbuL9aH0lFFRSGam7uMboYusqWv2dJPQNu+VleV9RvpRv/8xKl2NIsE\naABoae9B46etKT/Igd9r5kn3fkrdQCgmZrnd7kgABoCFCxeivb1du5YRUdpZXjsWdVMrUVbkhNUC\nlBU5UTe1Estrx8pmQsdmQMdT7pIoEymOhO12Oz7++GNMnDgRALB//364XK6UN4yIzEvujGCbtS8T\nWmykHM6AjqfcJVEmUwzCd911F26++WYUFxdDEAR0dHTgl7/8pR5tIyKTkzojOJzpvKehBe4uL0oK\nnagZXx75eTzlLokymWIQvvDCC7FlyxZ8+umnCIVCGD16NOx2ux5tI6I0JTdSVip3uWRuFbOlKWtI\nBuHHHnsMN998M+68807Rxx988MGUNYqIMoPYSFlNuctUJ24RmYVkEA6vAU+bNm3AYxapWnVERArC\niVutIoGYpSsp20gG4draWgBAU1MTbrjhhn6PcU2YiBLlyLUpJm4RZQvJILxhwwa0trZi27Zt+PTT\nTyM/DwaD2LdvH2677TY92kdEGUgpcYsoW0gG4QULFqCxsRG7du3qNyVts9nw/e9/X5fGEVFmkkvc\nIsomkkG4uroa1dXVqKurQ2HhuUofgiDgxImB00hERPGS2uJElC0Utyi9/PLLeOihh9DTc64M3YgR\nI7B169aUNoyIiCjTKZameeKJJ/Diiy9i4cKFeP3113H33Xdj8uTJerSNiCiCJS4pEymOhMvKyjBy\n5EhMmDABDQ0NuOaaa/Dss8/q0TYiIpa4pIym+Bfscrmwa9cuTJgwAdu3b0dzczO8Xq8ebSMiipS4\nbO30QcC5Epebth0xumlESVMMwj/5yU+wfft2zJkzB+3t7bjsssuwcuVKPdpGRFlOqcQlp6Yp3SlO\nR7/00kuR0pWPPfZYyhtERBTGEpeU6RRHwtu3b4cgCHq0hYion3jOJiZKR4oj4eLiYlx66aWYOHEi\nHI5zf/A8wIGIUk1tiUtfIMiiH5SWFIPwN7/5TT3aQUQZKtkAKVfiUipz+qZlNVp3gyglGISJKCW0\n2lokV+Jy49aGfqPkcOZ0nsuOxbPP07pLRJrjJjsiSgmttxaFS1xGT0FLZU7v2n+amdOUFhiEiUhz\nemwtksucbmnvQUe3+GNEZqI4HX3q1Kl+/7ZYLHA4HCgtLU1Zo4govemxtSicOd0q8j7lxS5mTlNa\nUAzCN954Iw4fPozx48dDEAQcPnwYFRUVsNlsuO+++zBz5kw92klEaUQuQGq1tUguc3rGpGHMkqa0\noDgdPWTIEPzlL3/B888/j7/97W/YvHkzJk2ahKeffhobNmzQo41ElGbCAVJM9NaiZC2vHYu6qZUo\nK3LCagHKipyom1qJ6xZN1OT1iVJNcSR88uRJTJo0KfLvCRMm4NixYxg2bBhCoVBKG0dE6Utua5FW\npDKnbTbdW0VGAAAgAElEQVSmu1B6UAzCI0eOxIYNG/Bv//ZvCIVCeOmll/ClL30Je/bsgZUnmBCR\nBLmtRVoLZ04TpRvFKLpu3ToEg0H88Ic/xJ133olQKIQHHngAx48fx09/+lM92khEaSx2a5EUufOC\neZYwZSrFkXBBQQGuv/56fPWrX0UoFMKFF16IgoICXHnllXq0j4gynFxRDwC6niXM8pekN8UgvGPH\nDtx111248MILEQqFcM899+D+++/HvHnz9GgfEWW4cFGPsHBRjzCpx1bUjdesDVpV9yKKl2IQ/tWv\nfoWNGzdi5MiRAIDjx4/jpptuYhAmoqTJFfX48FAzLBbx39vT0IIlc6s0G63K3QhoGeyJYine4vX2\n9kYCMNCXqMWsaCLSgnxRD+WCH9ESXTfWo7oXkRTFkfDw4cPx1FNPYenSpQCA5557DiNGjEh5w4go\n88kX9XDAYoFkwQ+XIwdNbg8K8ux4YcfRflPJsyePwKKZo1RNJetR3YtIimIQvv/++3Hffffht7/9\nLQRBwIwZM/Czn/1Mj7YRUYaTq3p10YS+Yh9ij+U5c/Czp95DW6cPDrsVXv+52bnWTh/+vuMoPD1+\nVVPJelT3ShYTxjKXYhAuKyvDww8/rEdbiCgLqSnqEf1YnjMHx5u6I49FB+BoateN5W4Eoqt7GREI\nk0kYY+BOD5JBuLa2FhaprAgAb7zxRkoaRETZJbqoR7PbA1gsqCh2wWa1whcIom5KJRbNOg89vl64\nHH0jYDXimUqWuxEwMnM6kYQxZnqnF8kg/PTTT+vZDiLKYsFQCJvfaowEjpJCO/Jddni8gX6BZF7N\nCMn121jxTCXLVffauLXBkMxppYQxqVE+M73Ti2QQZvIVEeklNnC0dfnR1uWP/DscSILBkOT6baxE\nDoqILX+ZaCDUQiIJY0a2N1HZPm2uuCZMRJRKcoEjVn1jG6rHlmP7hycHPOa02+APBFFS6MTsycOx\naOaopNtmZOZ0Iglj6ZTpzWnzPgzCRGQoucARy93lRd2USgDA3oYWtJ/1ofSL9dvFc0aj2xPAoAIH\nKocXo7m5K+m2GZk5rTZhLFo6ZHqHcdq8j2QQfuGFF2R/cfHixZo3hoiyj1zgiFVc4MDWD06g/kgL\n3N0+FBfYUV1VGhk95TlyNW1bIoFQS/EeB2l0e9VKx2nzVJEMwv/6178AAMeOHcNnn32GuXPnwmaz\n4e2338bYsWMTDsIlJXnIyUnfD7eiotDoJugmW/qaLf0EzNvX2ZNH4O87jio+r7jQ0W8qur3bj+17\nTqGwwInrF1/Q77la9fWmZTXIc9mxa/9ptLT3oLzYhRmThuG6RRN1Obf41n+fAq+/F+5OH0qKHHDa\nB162o/tqdHvVON1yFm1d0tPmNnsuKsrzBzxm1r/fZFgEQRDknrBq1So88sgjKC0tBQB0dHTgxhtv\nxDPPPJPQG2oxRWSUiorCtG5/PLKlr9nST8DcfT23Ptg34isucCDflQuPNwB3lw8lhU5UV5WivrFV\ndMRcVuTEz6+fHhk9paKvZk0gkuqrWdsL9LXt7t/vUvVdhpn571cNqRsIxTXhpqYmFBcXR/7tcrnQ\n3KwuiYKISA2pLULRgaSj24c395wS/X09ko5iM6fNzsztTZdpcz0oBuGvf/3r+M///E8sWLAAgiDg\nlVdewWWXXaZH24goy8QGjuh/p1PSkZa0GtGabWQc73p3plIMwnfeeSe2bNmC3bt3w2Kx4LrrrsP8\n+fP1aBsRUUQ8oyevvxdNbk8kMJsp+Kil1RYes24FkiuQkk1UbVEqLy/H2LFjsWTJEuzbty/VbSIi\nEqU0egoHnPrGVjS5e+C0WwFY4PMHTRN81NJqC4/ZtwKZedpcD4pB+L//+7+xdetWNDU14bLLLsM9\n99yDpUuX4lvf+pYe7SMiilAaPcUGnNjTlcwUfORotYWHW4HMT/F28G9/+xuefPJJuFwuFBcX47nn\nnsPmzZv1aBsRkajw6Ck6gKitvLWnoQW+QDCVzUuamspXer4OpY5iELZarbDb7ZF/OxwO2Gy8cyIi\nc1FbeSsdgk84CU1MPEloWr2O2fgCQTS5Paa/mVJDcTp62rRpeOihh9DT04OtW7di06ZNmD59uh5t\nIyJSTW3lrXQIPlpt4cm0rUBmTTJLhm3NmjVr5J4wa9YsnDx5Ej09PThw4AAuvvhifO9734M1wQ57\nPH7lJ5lUfr4jrdsfj2zpa7b0E8j8vubYrGjp8OLoqU7Z582+YChqxlXo1CplvkAQbZ1e5ORYkRNV\n0eor55Wgx9eLjm4/fP5elBY5MfuCoVheOxbWqLPelb5Xta9jdvn5Djz5j4+x9f0T6PH1jYB7fEEc\nPdWJHl8vLhhTZnAL5eXni9/4KVbM+tOf/oRrrrmm35T0+vXrcfvttyfUkHSveJLO7Y9HtvQ1W/oJ\nZEdfo7Ojm909cNj7Rnrh05XCmdSpGjXFsxdX7ahO6TXVfq9m2yccr8JBLnz3wa1xVdkyk4QrZj38\n8MN4+eWX8fjjj2PIkCEAgJ07d2rbOiIiDYSzp29Y4kLjp6267RNWE1Bjg6DarUNabeFJ961A7s70\nOaYxHopBePTo0bjhhhtwzTXXYO3atZg6daoe7SIiSpjTntPvgqzVxVlqNCkXUJfXjh0QoKurylDf\n2Cr6Htw6JK6kKDMrpikGYYvFgrq6OlRWVuLWW2/Ftddei9xcbY8LIyIyM7mRbm9QkN2LGwwJ/U5+\nau30YbtEDWwgvUd1qeS052RUklmYYhAOLxmff/75ePbZZ3HrrbfiwIEDKW8YEZFZyI1066ZUSk6T\ntnV6sbehRfQxqwUIiWTkpPOoLtUysd60YhD+9a9/Hfnv0tJS/OlPf8Krr76a0kYREZmFUtWpRbPO\nk5wmHVRgR7vEnmSxAAyk96gu1TKx3rRkEH7sscdw88034/HHHxd9/IorrkhZo4iIzEKp6lSPr1d6\nmnRcueQZyKWFDkweV476I60ZM6pLteg1+UyZrpcMwhMnTgTQV6yDiChbqTlCUW6a1GY7IhqgL5pQ\ngRV14+Gbl95bh/QQDIXw+xc+wjv7TmZMkY4wySB8/vnn49SpU6yORURZTW3VqSVzq/C16mGAxYKK\nYlfk50rrmOm+dUgPZj8JKhmSQXjlypWwWCzw+XxobW3FyJEjYbVacezYMYwaNYrrwkSUNeQCqdIe\n4Uxcx9RTpp8EJRmEt23bBgD4r//6L1xzzTWR/cH19fX4wx/+oE/riIhMQC6QbtzaoDhKS/dqVUZS\ncxJUOs8kKGZHNzY29ivQUV1djU8++SSljSIiMmPgip06VhqlLZ4zBi/sOJpRBw7oTc2afDpTDMJD\nhw7FI488goULF0IQBLz44os477zzdGgaEWWjVJyUk6qArjRKe/b1Bryz/0zkZ5m0lqmXTDsJKpZi\nEF6/fj0effRR3HbbbQCA2bNn48EHH0x5w4goOyWbhOMLBHG65SyCgSBybJa4A3o8AVtulFZc4MDB\nY27R38uEtUw9La8dizyXHe/sO5Vx27kUg/DatWsZdIlIF8kk4fQbQXf5UFroQJ4zF8ebuiPPkQvo\niYzA5UZp53+pBO9GjYKjZcJapp5sViuuX3wBLps20nRLFMlSnNtpaGjA2bNn9WgLEWU5NUk4UsIj\n6NZOHwShL+BGB+Boexpa4AsEpX8f5wL2pm1HZNu8vHYs6qZWoqzICaul71i9uqmVWPGNcSgtEl+v\nzIS1TCOE1+QzJQADKkbCVqsV8+bNw+jRo+FwnPuj+fOf/5zShhFR9kk0CUduBC0mdiSazAhcLnM6\nk9cySRuKQfj222/Xox1ERAkn4ciNoMXEBnQttsGIFd3Q+sABM2aMA+ZtVzpQDMLTpk3D//7v/8Lj\n8UAQBASDQZw4cYLlLIkoJRIJXHIjaDGxAT2REbiawJNIoQ6x15VbrzZSKjLZs41iEL777ruxe/du\ndHR0YMyYMTh48CAuuugiLF26VI/2EVGWSSRwyY2gRw4ugMfbKxvQ4xmBJ5rApTSSlntduYzxW/99\nSr/X0XNUmsnlJPWiGIR37tyJLVu24L777sO1116Lnp4erF27Vo+2EVEWi7emstwIujcoKAYmtSPw\nVAUeqdcNBkOob2wV/Z09DS3w+nsB6D8qzfRyknpRDMKDBw9Gbm4uqqqqcOjQIVx++eXo6urSo21E\nRKpFj6Bt9lwE/YFIELBZoRjQ1YzAUxV4ZF/3cAs6uv2ij7m7vHB3+pAD/UelmV5OUi+Kt0dDhgzB\nE088gZqaGvzlL3/B//zP/8DvF/+DICIymiPXhmHl+QmPwuS2wSSzhUqO3Ot2dPtRLJEVXlLoREmR\nQ/HmIHY7lhbC6+hS7eIWLHUUg/D999+PyspKVFdXY8GCBXjppZewZs0aHZpGRGQuqQo8cq9bWuTE\nhePLRR+rGV8Opz0nZTcHcsLr6FLtSqepaF8giCa3JyU3K0okp6NPnToV+e+amhqcOnUK8+fPx/z5\n83VpGBGR2aSqjrHS6/at61ok16uNOuRA6y1YejNDdnfc5wkfP34clZWV2LJliy4NJCLSUjzZw2LP\nTVXgiX7dtk4vBhXYUTOuPBIQ5NarjTrkIN3PSjZDdjfPEyairBDPqEfpuakIPDarFctrxyIYErC3\noQXt3T7UN7bCZjsSeV+5jPFkbg6S3dYUbya7GZglu5vnCRNRVpDbAnTJtFH9ApCaEVIqAs+mbUew\n/cOTsu8rJZGbAzNMxxrFLNndPE+YiDKe3Kjnrb2n8OaeU5EAtHjOaENGSFqNzOK5OTDDdKxRjFpH\nj6V4q7N+/Xp0dnbitttuww9/+EP09vbyaEMiSityo56QgH6nJm18/bDumcZKbUzF+xqxrSmV4s1w\nNkt2t+JIeNCgQfjJT36iR1uIiFIintrSBz9zGzJC0ntkZpbp2GQlM6VuhuxuxSD8/PPP46GHHkJn\nZycAQBAEWCwWHDhwIOWNIyLSglz2cKz2bh9mThyKd/afGfBYKkdIemc4m2U6NlnJTKmbIbtbMQj/\n5je/wdNPP43x47VZHygpyUNOTvqksMeqqCg0ugm6yZa+Zks/gezu603LapDnsmPX/tNodvfAYgVC\noYG/V17sws3/pwZlWw5h1/7TaGnvQXmxCzMmDcN1iybCZktdwlJ0G+N530S/19mTR+DvO46K/Hw4\nKocXJ/SaSrz+Xrg7fSgpcsBpVwxB/cT20+vvlayrXd/YihuWuFS/R2VcLdGORRAEQe4JK1aswMaN\nGzV7w+bm9K07XVFRmNbtj0e29DVb+gmwr2Hh7Thb3jveLxM5rG5qZWQE5QsE0ez2ABYLKopduo2S\n4tkylMz3em4qd+B0rNbZ0clmYov1s8ntwZ1P7IJYELNagAe+M8M0U+pSN0qKtwgTJ07ELbfcgtmz\nZ8PhODc9sXjxYu1aR0Skk3D28Iq6cbJVqIKhEDa/1WjI9h299t3qOR2bikzsTJhSVwzC3d3dyM/P\nx969e/v9nEGYiNKZUgBKx+07iRbdSHXQT1VhDKMqhWlJMQiLbUfyer0paQwRkd7EAlCyQSPZClTx\nMnvRjVRmYpshwzkZikF427ZtePjhh+HxeCAIAkKhELxeL95991092kdEpLtEg4ZRwdDso/ZUThub\nIcM5GYp/FQ8++CDuuusuVFVVYcOGDVi4cCEuu+wyPdpGRGSIRI8sDAfD1k5fvwIgm7YdSVlbvf5e\n0xfd0KMwhtw50GamGIQLCwsxY8YMTJ48GV1dXbj99tuxa9cuPdpGRGSIRIKGURWo3J2pqbSl9Rm7\ny2vHom5qJcqKnLBagLIiJ+qmVqbNtHGqKE5HO51OfPLJJ6iqqsLu3bsxY8YMBAIBPdpGRGSYeNca\njapAVVKk7VRvqqbU9Z421ntdPlGKQfgHP/gBHn74Yaxfvx6/+93vsGnTJixdulSPthERGSbeoGHU\ndhmnPUfTDGGp9WWPtxerLpmQdEBLdSa22ZPUYikG4ZKSEjzyyCMAgM2bN6Ojo4NHGRJR1lAbNIzc\nLqNVhrDclPrO/Wdw6JhbVUAzchSqJknNTKNkySD8wQcfIBQK4e6778b999+PcGGt3t5erFmzBlu2\nbNGtkURE6cCo7TJaTfXKTakDylnXRo9CldblF88Zgxd2HDXVKFkyCO/cuRO7d+9GU1NTZCQMADk5\nOVi+fLkujSMiSidGb5dJdqpX7WlTUnuljd4qpbQu/+zrDf0O5jDDVi7JIHzzzTcDAF544QVWxyIi\nioNeZSe1pva0KbFEM7lR6IeHmhOuihUPuZuI4gIHDh5zi/5eMlW7kiU7/t6+fTumTJkCANi6dSu+\n+93v4tFHH2V2NBFRhgpvJSotlE4kE0s0kxuFtnX58MyWQwiKHVulIbmtZed/qSQlW7mSJRmEn3zy\nSTz++OPw+Xw4ePAgfvSjH2H+/Plob2/HunXr9GwjERHpwBcIorXDiyVzq3D/d2Zg9qShos8TSzST\nK3ACAO/sP5PSoiVhUvuRV3xjXEIFWFJNcjr6xRdfxKZNm+ByubBhwwbU1tbi6quvhiAIWLhwoZ5t\nJCIiEb5AEKdbziIYCCY1lSqVULXq0vFwOXNUJZqpmcrWY9pXbl3ejIc9SAZhi8UCl8sFAPjXv/6F\nFStWRH5ORETG6Rc0u3woLUwuy1cpoUptotny2rHweHuxMyr5KVoqi5bEEluXN+NhD5JB2GazobOz\nEx6PBwcOHMDs2bMBACdPnkROjuL2YiIiShEts5DVnhilJnDarFasumQCDh1zm/KMX6Oz18VI3jJ9\n5zvfweLFi7Fs2TIsXboUgwcPxssvv4zVq1fjW9/6lp5tJCJKCa3rI+tB6xrVasptxkOPwxqSZabD\nHiSHtJdeeilqamrgdrtx/vnnAwDy8/Px85//HNOnT9etgUREWtOqqIQRlZe0rlGdinKbZpz2NSvZ\neeUhQ4ZgyJAhkX/PnTs35Q0iIkq1ZKdzjawMpVXQjL6B0DphyYzTvmbFxV0iyipq10DlGFkZKtka\n1WI3EJPHlWP+lBHYe7hV05FruhYt0RODMBFllWSnc7UI4slOYycz3St2A7Htg5Oom1qJn18/nSNX\nnTEIE1FWSXY6N5kgrtU0dvR0r82ei6A/oCpoqrmBSMeRq5lORYoXgzARZZVkp3OTCeJaT2M7cm2o\nKM9Hc3OXqudrndRltERuaswWsBmEiSjrJDOdm2gQ12IaO1kuRw6KCxxwi2w70mMPr9YBMJ6bGqOP\nWZTCIExEWSfZ7N1EgriRo9DoACQWgIHU7uFNRQCM96bG6GMWpTAIE1HWSjR7N5Egnor9uGrFBqBo\nZUWp38ObigAYz02NGWYhpBg3BiciSnNylZdiq3EZVUlKLgAVF9hxz+qpWFE3PmVTslpX+AqTO7Up\n9qZG66pgWuJImIhIQ3JTr0ZUkpILQJ1n/ejx9aIwz27I+yczDR/P2ryRsxBKGISJiDSkNPWqdyUp\nowNQKt9f7U1NshnxqcQgTESkEbVrj3pWkjI6ACWTTa50oxLP2rxZ61kzCBMRacSs+3CNDkDxvL/U\ndP5Ny2okX1/NTY1Z61kzCBMRacToqV8pRgegeN5fajo/z2XH4tnnJd0Ws9WzZnY0EZFGzH6WrtHn\n6Cq9v9x0/q79p9Pq3Ge1GISJKOvEbh/S0vLasaibWomyIieslr59uHVTKw1be0xlX7V+fbnp/Jb2\nnqS2EqX6c0gUp6OJKGvoUbrQ6KnfsFT31ePrxbOvN+DgMbdmry83nV9e7EpoOt+s5SrDGISJKGvo\nWbrQ6LXHVPU1HNTerj8Frz+k6evLZVLPmDQsoZsZs5arDDP+NoCISAepqtxkRqnsazioRQdgLV9f\najr/ukUT436tdPjOORImoqxg1u1DqaC2r/GeaiQX1MRePxFS0/k2W/xjxnT4zhmEiSgrmHX7UCoo\n9bUgz46NWxviXieVC2phxQUO+HtD8AWCSa2FazGdnw7fue5BuKQkDzk5xm+QTlRFRaHRTdBNtvQ1\nW/oJsK+zJ4/A33ccFfn5cFQOL9ajWSkRb19fe/+E5F7c6xdfIPk+hYNcqChxocndI/mcHn8Q9/5x\nNyqKXZgxaRiuWzQxoVGsmET+fs3+nesehN1uj95vqZmKikI0N3cZ3QxdZEtfs6WfAPsKAItmjoKn\nxz+gctOimaPS9rOJt68Lpo7AvU/uFn2td/adwmXTRgKA5DR1dVWZaOKUzQoEQ0CPrxcA0OTuwd93\nHIWnx69JAlSif79m+c6lbiAsgiAIurUCSNs/dIAXsUyULf0E2Ndo8a6Fmlm8fW1ye3DnE7sgduG3\nAJg1aajstqNzW376glpxgQPjRhbj8Il20anqsiInfn799KQ/52T/fo3+zqWCMNeEiSjrGL19SE+x\nfZVbJ3XYbXhn/5nIv8W284glTnV0+3DnE7tE398sCVBm/c65RYmIKIvIldaUIradJ7oEZTiwizFL\nApRZMQgTEWUZsb24syYNhc8vvm82PJqVYvaa2WbG6WgioiwjNqUMAIeOuRPezmP0cYnpikGYiChL\nxa6TSpWMVDOaNUvN7HTDIExERAC0Gc2aNQHKrBiEiYgIgPxo1ugtPpmKQZiIiPqJHs2a/SjAdMcg\nTEREksx+FGC6420MERGJkjs16cNDzaY4CjDdMQgTERnAFwiiye2JBLLYf5uB3KlJbV0+PLPlEIIh\n8XOFSR1ORxMR6Sh2jbWk0I58lx0eb8B0a65yJS4B4J39Z+By5nBaOgkcCRMR6Si8xtra6YMAoK3L\nj+NN3ZF/h9dcN207YnRTVZW4FCtpSeoxCBMR6URujTWWWYLb8tqxmDVpqOTjSiUtU8mMU/jx4nQ0\nEZFO5NZYY5nl9CGb1YpVl0xIqqSl1jJp21R6tZaIKI3JnTYUy0ynD5ntgIbYKf3YKfx0GiFzJExE\npJNwMBOrzxzLbKcPmeWABqVtU8GQgPojLWkzQmYQJiLSUWwwKy5wIN+VC483AHeXz7SnD5nlgAal\nbVPbPzwZ+Xc6FBZhECYi0pFUMEuX2sxGHdDg9feiye2By5EjuW3KagFCwsDf3dPQgiVzq0z5uTII\nExEZIDaY8fShPrE3I+EkrPrGVjS7e1Ba5ECeM1c0CIsFYMA8SW5iGISJiMhwUhnPgiDgjQ/6TzG3\ndvowcnABPN7eyPp0dVUp6htbk8rgNmI2gkGYiIgMJ3VQhNMuHgw93l7cs3oqeny9kaC5cWuDaNKb\nUpKbkVueGISJiMhQchnPXr/4NiN3lxc9vt5+U8yJZnAbeVIUgzARERkqniImYWJTzIlkcMvdAOiR\n0GXOjVNERJQ15IqYSE1Hy00xh5Pc1ARPuRsAPUpyMggTEZGh5CpyzbpgKOqmVmJwiQtWC1BW5ETd\n1ErN9lHL3QDoUbWM09FERGQ4ufVcm9WKG5a40Phpq+aZy3JVzPSoWsYgTESUxsJFLMxe5EOJ0nqu\n056Tsn2+RpbkZBAmIkpDYkUszF4nWQ0jipYYWZKTQZiIKA0Zua0mUxlxA5C+t0tERFlKaVtNOhzh\nR30YhImI0ozR22pIOwzCRERpxuhtNaQdBmEiojQjt69Wj201pB0mZhERpaHw9pn6xla0tPfouq2G\ntMMgTESUhsLbalJVxIL0wSBMRJTGUlnEglKPa8JERGQ4XyCIJrcn67ZXcSRMRESGCVf+2tPQjLZO\nX8ZU/lKLQZiIiAxjdOUvXyCoe6nKaBZBEAQ937C3N4icHCYPEBFlO6+/Fzeu24Ymd8+AxwaXuPDr\n/1sLpz01Y8VgMIQ//uNj7Np/Gs3tPagodmHGpGG4btFE2Gz6jcB1Hwm73R6931IzFRWFaG7uMroZ\nusiWvmZLPwH2NVOlc1+b3B40iwRgAGhp70Hjp62RpDOt+7lxa0O/EXiTuwd/33EUnh5/SkbgFRWF\noj/P/Al3IiIyJaMqf5mp9jaDMBERGcKoyl9mqr3NxCwiIjJMuMLXnoYWuLu8ulT+Co/AW0UCsd61\ntxmEiYjIMOHKX0vmVg3IUo7OXNZSeAQevSYcpnftbQZhIiIynCPXFknCEts7PHvyCCyaOUqzvcNG\njMDFMAgTEZGpiO0d1jpzWW4EricmZhERkWlonbmsVA4zPAI36vALjoSJiMg01GQuqzmwIl3KYZqn\nJURElPW02jscntJu7fRBwLlymJu2HdGwtcljECYiItPQYu+wmYpxKOF0NBERmYpY5vLsycOxaOYo\nVb+v1ZS2HhiEiYjIEFInGIllLlcOL1ZdO9pMxTiUMAgTEZGu1CZNRe8djoeZinEoYRAmIiJd6XGG\nsFmKcShhECYiSjOpKueoB6WkqSVzqzQZqZqlGIcSBmEiojShRznHVNM7aSrRKW29pMe3RkREontf\n/77jqOn2vsox6gxhs2IQJiJKA+m091WOUWcImxWno4mI0kA67X1Vki5JU3pgECYiSgPptPdVSbok\nTemB09FERGkgE6dxjT7ByAw4EiYiShPJlnMk82EQJiJKE8mWczSKVHlKYhAmIko7Zt/7GpYuZ/oa\niUGYiIhSQo/ylOmOtyJERKS5TNnXnGoMwkREpDk1+5qJQZiIiFKA5SnVYRAmIiLNZeK+5lRgYhYR\nEaUEy1MqYxAmIqKUYHlKZQzCRESUUumyr9kIXBMmIiIyCIMwERGRQRiEiYiIDMIgTEREZBAGYSIi\nIoMwCBMRERmEQZiIiMggDMJEREQGsQiCIOj5hr29QeTksGIKERGR7hWz3G6P3m+pmYqKQjQ3dxnd\nDF1kS1+zpZ8A+5qpsqWv6d7PiopC0Z/rPhImIiKiPlwTJiIiMgiDMBERkUEYhImIiAzCIExERGQQ\nBmEiIiKDMAgTEREZhEGYiIjIIAzCREREBmEQJiIiMgiDMBERkUEYhIni8Oqrr+Kqq67ClVdeiUWL\nFuEPf/hD5LFHH30U77//vibvU1tbixMnTiT9+2+88QYeeeSRpNtzxx134Pnnn0/6dYioP90PcCBK\nV+u5V6AAAAZxSURBVJ9//jkeeughPP/88ygpKcHZs2exatUqjB49GvPnz8d7772H6dOnG93MfubP\nn4/58+cb3QwiksAgTKSS2+1GIBCA1+sFAOTn52Pt2rVwOBx44YUXsH//ftx99914/PHH0dHRgV/9\n6lfwer3o7OzEnXfeibq6Otxxxx0oKCjAxx9/jM8//xw33ngjlixZgvb2dtx+++04c+YMqqqq4PP5\nAADd3d2466678Pnnn6OpqQkzZ87E/fffj927d2P9+vUIhUIYN24c7rzzTtHff/7557F7927cdNNN\nuPHGGyN9+eSTT3Drrbdi9erVWLduHXbv3o1gMIirrroKq1evhiAIWLt2Ld58800MHjwYwWAQ06ZN\n6/d5nDhxAt/+9rdRUlICp9OJxx57TLKtTzzxBJxOJxobGzFhwgRs2LABdrsdf/7zn/HMM8+gsLAQ\nY8aMwahRo3DzzTfjn//8Jx599FH09vaisrIS9913H0pKSnT6pol0JBCRavfcc4/wla98RViyZImw\nbt064cCBA5HHVq5cKezatUsQBEG4+eabhSNHjgiCIAg7d+4UrrjiCkEQBOHHP/6xcOONNwqhUEg4\nePCgMG3aNEEQBOGnP/2p8Mtf/lIQBEHYvXu3MH78eOH48ePCP/7xD+E3v/mNIAiC4PP5hLq6OuGj\njz4Sdu3aJUyZMkXo7OyU/f3NmzcLP/7xj/v14bXXXhOuuuoqwev1Chs3bhQeeOCByOuvXLlSeO+9\n94RXXnlFWLlypeD3+4XW1lZh9uzZwubNm/u9zvHjxyPvIwiCbFsvvPBC4fTp00IwGBSWLFkivPHG\nG8KBAweEBQsWCF1dXYLX6xWuvvpq4dFHHxVaW1uFK6+8UmhvbxcEQRCeffZZ4a677kr6uyMyI46E\nieLw05/+FN///vfx9ttv4+2338ayZcuwYcMGLFiwoN/z1q9fj+3bt+PVV1/Fvn37cPbs2chjs2fP\nhsViwfjx49He3g4A2L17N37xi18AAL761a9i5MiRAIArrrgC9fX1eOqpp3D06FG0t7fD4+k7k3v0\n6NEoLCyU/f1YBw8exNq1a/H000/D4XDg3XffxYEDB7Br1y4AgMfjwaFDh9DY2IgFCxYgNzcXpaWl\n+NrXvib6emVlZaisrFRs67hx4zB06FAAQFVVFTo6OvDZZ59h3rx5KCgoAABcfvnl6OzsxL59+3D6\n9Glce+21AIBQKIRBgwap+4KI0gyDMJFKb775JjweDxYuXIglS5ZgyZIl+Otf/4rnnntuQBBesWIF\npk+fjunTp2PmzJn40Y9+FHnM4XAAACwWS+RnFosFQtTR3jabDQDw9NNPY8uWLVi2bBlmzZqFhoaG\nyPOcTqfi70dra2vDLbfcggceeADDhw8HAASDQdx+++2R9re1tSE/Px/r1q3r93o5OeKXiug2yLU1\n3OfotlqtVoRCoQGvGQwGcdFFF+G3v/0tAMDn8/W7iSHKJMyOJlLJ6XTiF7/4RSRrWRAEHDhwAF/+\n8pcB9AW+YDCI9vZ2fPrpp7j11lvxta99DW+88QaCwaDsa8+cORMvvvgiAKC+vh7Hjh0DALzzzjtY\nvnw5rrzySvh8Phw8eFA0cEn9flggEMCtt96KVatW9UsemzFjBv76178iEAjg7NmzWLFiBfbu3YuZ\nM2filVdegd/vR0dHB3bs2KH4+ahta3Sb33rrLXR3d8Pv9+O1116DxWLB5MmTsXfvXnzyyScAgN/8\n5jdYt26d4vsTpSOOhIlUmjFjBm666SZ897vfRSAQAADMmTMnkvA0Z84c3HvvvXjooYewdOlSXH75\n5cjJycGMGTPg9XojU7NibrnlFtxxxx24/PLLMWbMmMh08n/8x39gzZo1+N3vfoeCggLU1NTgxIkT\nGDVqlKrfD3v11VexZ88e9PT0YPPmzRAEAbNmzcJtt92Gzz77DN/85jfR29uLq666KhKkP/roI1xx\nxRUoLy9HVVWV4uejtq1h48ePx7XXXovly5cjLy8PJSUlcDgcqKiowAMPPIAf/OAHCIVCGDJkCNav\nX6/4/kTpyCJEzzkREenkk08+wVtvvYXVq1cDAL73ve/h6quvRm1trbENI9IRR8JEZIgRI0ZERtsW\niwUXX3wx5s2bZ3SziHTFkTAREZFBmJhFRERkEAZhIiIigzAIExERGYRBmIiIyCAMwkRERAZhECYi\nIjLI/wcrR/3zLrMINQAAAABJRU5ErkJggg==\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": true }, "outputs": [], "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": {}, "outputs": [ { "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": true, "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": true }, "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": { "scrolled": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 30000/30000 [01:49<00:00, 274.28it/s]\n" ] } ], "source": [ "SAMPLES = 20000\n", "BURN = 10000\n", "\n", "with model:\n", " step = pm.Metropolis()\n", " trace = pm.sample(SAMPLES, step, tune=BURN, random_seed=SEED)\n", " " ] }, { "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": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAFyCAYAAADyNHPCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVXXi//H3ZZPVVETL9aGWWi4pY5qVmZCVk5i4m2KL\nheZXLbMU+LoFimZqiTW2jo1mRimaTdY0oKVjSuY3LNe+VpZbSgkIyFeEe35/+PPOoFwPyoXLPb2e\nf917zofPeUvp+57lnmMzDMMQAACwLC93BwAAAFWLsgcAwOIoewAALI6yBwDA4ih7AAAsjrIHAMDi\nfNwdoKpkZ+dXeGzduoHKyTlThWlcz9Mye1peyfMye1peiczVwdPySmS+WmFhIU7XsWcvycfH290R\nrpinZfa0vJLnZfa0vBKZq4On5ZXIXBUoewAALI6yBwDA4ih7AAAsjrIHAMDiKHsAACyOsgcAwOIo\newAALI6yBwDA4ih7AAAsjrIHAMDialTZ79q1SzExMZcs37hxowYOHKihQ4fq/fffd0MyAAA8V415\nEM4bb7yh9evXKyAgoMzyc+fOae7cuVq9erUCAgI0fPhw9erVS2FhYW5KCgCAZ6kxe/bNmjXTkiVL\nLln+ww8/qFmzZrrmmmvk5+enP/3pT/r666/dkBAAAM9UY/bs7733Xh05cuSS5QUFBQoJ+fdj+4KC\nglRQUGA6X926gfLx8VbU5A9dku+jhQ+4ZB5XutzjDMvj7t/FleatCTwts6fllchcHTwtr0RmV6sx\nZe9McHCwCgsLHe8LCwvLlL8zrn6ucHZ2vkvnq6ywsBC3Zbqa7boz79XytMyellcic3XwtLwSmSuT\nwZkacxjfmVatWunnn39Wbm6uiouL9fXXX6tz587ujgUAgMeosXv2H330kc6cOaOhQ4cqLi5Oo0eP\nlmEYGjhwoBo2bOjueAAAeIwaVfZNmjRxfLUuKirKsTwiIkIRERHuigUAgEer8YfxAQBA5VD2AABY\nHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxl\nDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8A\ngMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDF\nUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2\nAABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAA\nWFyNKHu73a4ZM2Zo6NChiomJ0c8//1xm/VtvvaUBAwZo4MCB+uc//+mmlAAAeCYfdweQpPT0dBUX\nFys1NVVZWVmaN2+eli5dKkk6ffq0VqxYoc8++0xFRUXq37+/evfu7ebEAAB4jhqxZ79z50716NFD\nktSpUyft3r3bsS4gIECNGjVSUVGRioqKZLPZ3BUTAACPVCP27AsKChQcHOx47+3trZKSEvn4nI93\n3XXX6f7771dpaanGjBlToTnr1g2Uj4+3yzKGhYW4bC5XcVemq91uTfwdmvG0zJ6WVyJzdfC0vBKZ\nXa1GlH1wcLAKCwsd7+12u6PoN2/erJMnTyojI0OSNHr0aIWHh6tjx46XnTMn54xLM2Zn57t0vsoK\nCwtxW6ar2a47814tT8vsaXklMlcHT8srkbkyGZypEYfxw8PDtXnzZklSVlaWWrdu7Vh3zTXXyN/f\nX35+fqpVq5ZCQkJ0+vRpd0UFAMDj1Ig9+969e2vr1q0aNmyYDMNQcnKyli1bpmbNmikyMlJffvml\nhgwZIi8vL4WHh+v22293d2QAADxGjSh7Ly8vJSYmllnWqlUrx+uJEydq4sSJ1R0LAABLqBGH8QEA\nQNWh7AEAsDjKHgAAi6PsAQCwOMoeAACLo+wBALA4yh4AAIuj7AEAsDjKHgAAi6PsAQCwOMoeAACL\no+wBALA4yh4AAIuj7AEAsDjKHgAAi6PsAQCwOMoeAACL83G2Ij4+/rI/OHfuXJeHAQAArud0z75r\n167q2rWrCgsLdfLkSd1666264447dPr0aRmGUZ0ZAQBAJTjds4+OjpYkvfvuu0pNTZWX1/nPBX36\n9NGQIUOqJx0AAKg003P2+fn5ys3Ndbz/7bffdObMmSoNBQAAXMfpnv0FY8eOVb9+/RQeHi7DMJSV\nlaXp06dXRzYAAOACpmXfv39/3Xbbbfrmm29ks9k0a9YshYaGVkc2AADgAqaH8YuLi5WWlqaMjAx1\n795dq1atUnFxcXVkAwAALmBa9omJiTpz5oz27t0rHx8f/fLLL0pISKiObAAAwAVMy37Pnj16+umn\n5ePjo4CAAD3//PPav39/dWQDAAAuYFr2NptNxcXFstlskqScnBzHawAAUPOZXqA3atQoPfLII8rO\nztacOXOUnp6ucePGVUc2AADgAhW6Gr99+/bKzMxUaWmpli5dqrZt21ZHNgAA4AKmZT9hwgQtWbJE\n119/vWPZQw89pL/97W9VGgwAALiG07IfP3689u3bpxMnTigyMtKxvKSkRNddd121hAMAAJXntOzn\nzZun3NxczZkzR9OnT3c8/MbHx4eb6gAA4EGcXo0fHBysJk2aaNasWVq5cqUaN24su92uhQsXKi8v\nrzozAgCASjD96t0zzzyjpk2bSpIaNmyoLl266Nlnn63yYAAAwDVMyz43N1fDhg2TJPn5+WnIkCHK\nycmp8mAAAMA1TMs+ICBAX3zxheP9l19+qYCAgCoNBQAAXMf0q3fPPfecnn32WU2ZMkU2m03XXnut\n5s+fXx3ZAACAC5iW/Y033qi///3vysnJka+vr4KDg6sjFwAAcBHTw/hHjx7VI488oqFDh6qoqEij\nRo3SkSNHqiMbAABwAdOynzFjhkaPHq3AwEDVr19fffv21dSpU6sjGwAAcAHTss/JydEdd9wh6fwT\n8IYMGaKCgoIqDwYAAFzDtOz9/f3166+/Oh5r+/XXX8vPz6/KgwEAANcwvUAvPj5eY8aM0S+//KJ+\n/frp9OnTWrx4cXVkAwAALmBa9h06dNDq1at16NAh2e12tWjRgj17AAA8iGnZHzt2TElJSdq+fbt8\nfX115513KiEhQfXq1auOfAAAoJIqdG/822+/XVu2bFFGRobat2/P1fgAAHgQ07IvKCjQyJEjFRwc\nrJCQED388MM6ceKE0/HlPRHv6NGjlUsJAACummnZd+7cWR9++KHj/eeff66bbrrpknHHjx/XsWPH\nNGLECMfrY8eO6fDhwxo9erRrUwMAgAozPWf/2WefKTU1VTNnzpTNZlNRUZEkad26dbLZbNq3b58k\nKSUlRZmZmTp58qRGjBjx7w34+Oiuu+6qmvQAAMCUadlv27btkmUlJSXy8Sn7o3PnzpUkvf7664qN\njXVRPAAAUFmmZT9p0iQlJiYqJCREknTgwAHFxcVp7dq15Y4fMmSI3n33XeXk5MgwDMfy8ePHuygy\nAAC4EqZl36ZNG0VHR2vGjBnavXu3Vq9efdmr8cePH6969erphhtucNx1DwAAuI9p2Y8dO1YtW7ZU\nbGys6tevrzVr1qhhw4ZOx+fl5emdd95xaUgAAHD1TK/GX7x4sZKTk7Vo0SINGzZMI0eO1MaNG52O\nb926tXbv3u3SkAAA4OqZ7tn/8MMPSktLc9wxLzIyUvHx8YqIiCgzLiIiQjabTf/3f/+nDRs2qGHD\nhvL29pZhGLLZbMrIyHC6DbvdrlmzZunAgQPy8/PT7Nmz1bx5c8f6L774Qq+88ook6aabbnJ8MwAA\nAJgzLfuUlBSdOXNG+/fvV5s2bdS8eXOtXr36knErVqy46hDp6ekqLi5WamqqsrKyNG/ePC1dulTS\n+Zv6vPDCC1q+fLnq1aunN954Qzk5OdyuFwCACqrQV+9mzJih0tJSpaamqm/fvlq4cKHjGfcXNG7c\nWNL5799fzN/fX4WFhWrdunW529i5c6d69OghSerUqVOZ0wDffPONWrdureeff16HDx/W4MGDKXoA\nAK6AadkvWrRI7777rh5//HGFhYVp5cqVevrppy8p+wsyMjK0d+9e3X333ZLO33GvQYMGOnPmjKKi\novTwww9f8jMFBQUKDg52vPf29nZ8lz8nJ0eZmZlat26dAgMDNWLECHXq1EktWrS4bO66dQPl4+Nt\n9sersLCwEJfN5SruynS1262Jv0MznpbZ0/JKZK4OnpZXIrOrmZa93W5XWFiY4/31119/2fHZ2dla\nu3atateuLUmaMGGCxo4dq9TUVA0YMKDcsg8ODlZhYWGZbV64aU+dOnXUoUMHR4YuXbpo3759pmWf\nk3PG7I92RbKz8106X2WFhYW4LdPVbNedea+Wp2X2tLwSmauDp+WVyFyZDM6YXo1/7bXXatOmTbLZ\nbDp9+rSWLl2qRo0aOR2fk5OjoKAgx/tatWopLy9PPj4+Ti+qCw8P1+bNmyVJWVlZZQ73t2/fXt9/\n/71OnTqlkpIS7dq1y/QDBwAA+DfTPfvExETNmTNHx48fV+/evdWtWzclJiY6HX/PPffooYceUp8+\nfWS32/XZZ58pMjJS69atK3OE4D/17t1bW7du1bBhw2QYhpKTk7Vs2TI1a9ZMkZGRmjx5sh577DFJ\n0n333ef03D8AALiUadmHhoZq0aJFFZ5w8uTJ2rRpk7Zu3Spvb2899thj6tmzp7KysrRw4cJyf8bL\ny+uSDxCtWrVyvL7//vt1//33VzgDAAD4N9Oyr6g9e/aoXbt22rFjh4KDg3Xvvfc61u3YsUO33HKL\nqzYFAACugMvK/r333lNSUpJSUlIuWWez2bR8+XJXbQoAAFwB07J/8cUXNWnSJNOJkpKSJFXu5joA\nAMD1TK/G37RpU5lH1Zo5evSoHnnkEd1zzz3Kzs7WqFGjdOTIkUqFBAAAV890z75OnTq677771K5d\nO9WqVcuxfO7cueWOnzFjhkaPHq0FCxaofv366tu3r6ZOnaqVK1e6LjUAAKgw07KPjo6+oglzcnJ0\nxx13aMGCBbLZbBoyZAhFDwCAG5kexo+Ojla7du1UWFiovLw8tW3b9rIfAPz9/fXrr786bqDz9ddf\ny8/Pz3WJAQDAFTEt+3Xr1mncuHE6cuSIjh07pvHjx5f71LsL4uLiNGbMGB06dEgPPPCAnnnmGU2b\nNs2loQEAQMWZHsZftmyZPvjgA9WtW1eSNHbsWI0aNUqDBg0qd3yTJk20evVqHTp0SKWlpWrZsiV7\n9gAAuFGFHoRzoeglqV69ek7vcS9JI0eOVO3atdWzZ0/16tWLogcAwM1My75NmzaaM2eOY09+9erV\natu2rdPxGzZs0JEjR7R582YtXrxYhw4dUrdu3TRr1iyXhQYAABVnes5+9uzZ8vX1VUJCguLj4+Xj\n46OZM2c6HW+325WTk6OioiIZhqGSkhKdOnXKpaEBAEDFme7Z+/v7a8qUKRWe8JZbblFAQIAefPBB\nPfXUU5c9CgAAAKqe07KPjo7W2rVr1bZt2zLn6A3DkM1m0759+8r9ucWLF2v79u3asmWLtm7dqi5d\nuqhr1666/fbbXZ8eAACYclr2a9eulXT+q3dXsnd+xx136I477tDp06f1z3/+U6+99pqWL1+ub775\npvJpAQDAFTM9jD9p0iR98sknFZ5wwYIF2r59u/Lz89WjRw9Nnz5d3bp1q1RIAABw9UzL/vrrr9fL\nL7+sm2++Wf7+/o7lzp5PHxoaqvnz56tly5auSwkAAK6aadnn5uYqMzNTmZmZjmWXez79I4884rp0\nAACg0kzLnufTAwDg2Uy/Z8/z6QEA8GymZX/h+fSBgYFlnk/vTF5enqZNm6ZRo0YpNzdX8fHxysvL\nc2loAABQcaZlf+H59JIcz6cvKChwOn769Onq0KGDcnNzFRgYqAYNGujZZ591XWIAAHBFTMv+Sp9P\nf+TIEQ0dOlReXl7y8/PTpEmT9Ouvv7ouMQAAuCKmF+hdeD79L7/8ogceeEB5eXlavHix0/He3t7K\nz893fDg4dOiQvLxMP1MAAIAqYlr2zZs3v+T59NnZ2U7HT5w4UTExMTp+/LjGjRunrKwsJScnuzQ0\nAACoOKdlf/z4cRmGodjYWL3xxhsKCgqSJJ04cUKPP/64Pv3003J/LiwsTH/961/17bffqrS0VImJ\niapfv37VpAcAAKacln1KSooyMzN18uRJjRgx4t8/4OOju+66y+mEF26ve7kxAACg+jgt+7lz50qS\nXn/9dcXGxlZ4wiu9vS4AAKhapufs7733Xq1fv15RUVGaOXOm9uzZo+eee07t27cvd/yV3l4XAABU\nLdOyT0hI0ODBg5WRkaGffvpJ8fHxmj17tt57771yx3N7XQAAahbTsj979qz69++v//7v/1ZUVJS6\ndOmi4uJip+NjYmIcX7v7T+zZAwDgHqZl7+3trX/84x/6/PPP9eSTTyo9Pf2y35ufMGGC43VJSYky\nMjJUu3Zt16QFAABXzLTsExMT9fbbb2vmzJlq0KCBPv74Y82ePdvp+K5du5Z5f9ttt2nw4MF68skn\nK58WAABcMdOyb9OmjR5++GHt2LFDb7/9tmJjY9W2bVun448dO+Z4bRiGDh48qNzcXNekBQAAV8y0\n7NetW6eXX35Zd999t+x2u8aPH68nnnhCgwYNKnf8yJEjZbPZZBiGbDab6tWrp2nTprk8OAAAqBjT\nsl+2bJk++OAD1a1bV5I0duxYjRo1ymnZp6WlqU6dOmWWHT161AVRAQDA1TB9Qo3dbncUvSTVq1ev\n3Kvtjx8/rmPHjmnkyJGO18eOHdPhw4c1evRo16YGAAAVVqFz9nPmzHHsya9evbrcc/ZXe3tdAABQ\ntUzLfvbs2VqyZIkSEhJkGIa6deummTNnXjLucrfXLS0tdVFcAABwpUzL3t/fX48++qg6duwoHx8f\ndenSRcHBwU7H79u3T/n5+QoJCZEkHThwQHFxcVq7dq3rUgMAgAozPWf/4Ycfql+/fvr444+Vlpam\nvn376osvvnA6vk2bNoqOjtbmzZv1l7/8RU888YTGjh3r0tAAAKDiTPfsly5dqrS0NDVs2FDS+Svr\nx44dq549e5Y7fuzYsWrZsqViY2NVv359rVmzxvGzAACg+pnu2QcFBSksLMzxvnHjxvL19XU6fvHi\nxUpOTtaiRYs0bNgwjRw5Uhs3bnRNWgAAcMVM9+w7dOigxx9/XAMHDpS3t7c++eQTNWjQQOvWrZMk\n9e/fv8z4H374QWlpaapXr54kKTIyUvHx8YqIiKiC+AAAwEyFnnrXoEEDbdmyRZIUEBCggIAAx/Pq\nLy77lJSUMu9vvPFGrVmzxlV5AQDAFTIt+wtfqTMTHR2ttWvXqm3bto7b5V5gs9m0b9++q08JAACu\nmmnZV9SFr9bt37/fVVMCAAAXcFnZX3Dq1Cl9/PHHysvLK7N8/Pjxrt4UAACoANOr8bdu3XrJss8+\n+8zp+Mcff1x79+6tXCoAAOAyTvfsN2zYoOLiYqWkpGjixImO5efOndPrr7+ue+65x+mkFT3PDwAA\nqp7Tsi8sLNT//M//qLCw0HHlvSR5e3tr0qRJTie8++679cEHH+jWW2+Vt7e3Y3mjRo1cFBkAAFwJ\np2U/ePBgDR48WNu2bVP37t0dywsKCi57b/wzZ84oOTm5zGNxbTabMjIyXBQZAABcCdML9IqKivTC\nCy9o3LhxGjRokE6dOqWpU6dqwIAB5Y7ftGmTtm3bJn9/f5eHBQAAV870Ar1XXnlFUVFR2rBhgzp2\n7KiNGzfqnXfecTq+cePGl1yJDwAA3KdCX71r27atlixZon79+ikoKEjnzp1zOvbcuXO6//77dcMN\nN5S5h/7y5csrnxYAAFwx07KvX7++kpKStHv3br3wwguaN2/eZS+2u5rH2drtds2aNUsHDhyQn5+f\nZs+erebNm18yJjY2VpGRkRo+fPgVbwMAgD8q07JfuHCh0tPT9dBDDykwMFBNmza97A1yunbtesUh\n0tPTVVxcrNTUVGVlZWnevHlaunRpmTEvvfQSpwcAALgKpufsg4OD5eXlpTVr1qioqEhBQUGXvRr/\nauzcuVM9evSQJHXq1Em7d+8us/7TTz+VzWbTnXfe6dLtAgDwR2C6Z79gwQL9+uuv2rNnjx5//HGt\nWbNG+/fvV1xcnMtCXPx1Pm9vb5WUlMjHx0fff/+9/v73vyslJUWvvPJKheesWzdQPj7e5gMrKCws\nxGVzuYq7Ml3tdmvi79CMp2X2tLwSmauDp+WVyOxqpmX/r3/9S2vXrlV0dLSCg4O1bNky9evXz6Vl\nHxwcrMLCQsd7u90uH5/z0datW6cTJ07ooYce0tGjR+Xr66vGjRub7uXn5JxxWT5Jys7Od+l8lRUW\nFuK2TFezXXfmvVqeltnT8kpkrg6ellcic2UyOGNa9l5e54/022w2SVJxcbFjmauEh4dr06ZN+vOf\n/6ysrCy1bt3asW7KlCmO10uWLFH9+vU5nA8AwBUwLfv77rtPTz31lPLy8vT2229r/fr16tu3r0tD\n9O7dW1u3btWwYcNkGIaSk5O1bNkyNWvWTJGRkS7dFgAAfzSmZR8bG6stW7aoUaNGOn78uCZMmKBe\nvXq5NISXl5cSExPLLGvVqtUl4yZMmODS7QIA8EdgWvY7duyQv7+/IiIiJJ0/nP/dd9+pefPmql27\ndpUHBAAAlWNa9q+88op2796t7t27yzAMffXVV2rcuLEKCgr05JNPuvyQPgAAcC3TsjcMQ+vXr3fc\nNe/EiRNKSEjQihUrFBMTQ9kDAFDDmV5Wf/LkyTK3x23YsKFOnjyp4OBgGYZRpeEAAEDlme7Zd+7c\nWZMnT1ZUVJTsdrs+/vhjde7cWZ9//rkCAwOrIyMAAKgE07JPTEzUe++9p9TUVHl7e6t79+4aOnSo\ntm7dqvnz51dHRgAAUAmmZT9mzBi99dZbGjlyZJnlPXv2rLJQAADAdUzP2RcVFen48ePVkQUAAFQB\n0z37nJwcRUREKDQ0VLVq1ZJhGLLZbMrIyKiOfAAAoJJMy/7NN9+sjhwAAKCKmJZ9WFiYvvjiC8dT\n6UpLS3XkyBE9+eSTVR4OAABUnmnZP/3008rLy9Mvv/yiLl26KDMzU+Hh4dWRDQAAuIDpBXoHDhzQ\n8uXL1bt3bz322GNatWqVjh49Wh3ZAACAC5iWfWhoqGw2m1q0aKEDBw6oadOmOnfuXHVkAwAALmB6\nGP+GG25QUlKShg8frmeeeUYnT57kNrkAAHgQ0z37WbNmqU+fPrr++us1YcIEnTx5UosWLaqObAAA\nwAVMy97b21tdunSRJEVGRmratGlauXJllQcDAACuYVr25Vm/fr2rcwAAgCpyVWXPOXsAADzHVZW9\nzWZzdQ4AAFBFnF6NHxMTU26pG4ahs2fPVmkoAADgOk7LfsKECdWZAwAAVBGnZd+1a9fqzAEAAKrI\nVZ2zBwAAnoOyBwDA4ih7AAAsjrIHAMDiKHsAACyOsgcAwOIoewAALI6yBwDA4ih7AAAsjrIHAMDi\nKHsAACyOsgcAwOIoewAALI6yBwDA4ih7AAAsjrIHAMDiKHsAACyOsgcAwOIoewAALI6yBwDA4ih7\nAAAsjrIHAMDiKHsAACyOsgcAwOIoewAALI6yBwDA4ih7AAAsjrIHAMDiKHsAACyOsgcAwOIoewAA\nLM7H3QEkyW63a9asWTpw4ID8/Pw0e/ZsNW/e3LH+7bff1scffyxJ6tmzp8aPH++uqAAAeJwasWef\nnp6u4uJipaamavLkyZo3b55j3eHDh7V+/Xq99957Sk1N1b/+9S/t37/fjWkBAPAsNWLPfufOnerR\no4ckqVOnTtq9e7dj3bXXXqs333xT3t7ekqSSkhLVqlXLLTkBAPBENaLsCwoKFBwc7Hjv7e2tkpIS\n+fj4yNfXV/Xq1ZNhGJo/f75uuukmtWjRwnTOunUD5ePj7bKMYWEhLpvLVdyV6Wq3WxN/h2Y8LbOn\n5ZXIXB08La9EZlerEWUfHByswsJCx3u73S4fn39HO3v2rBISEhQUFKSZM2dWaM6cnDMuzZidne/S\n+SorLCzEbZmuZrvuzHu1PC2zp+WVyFwdPC2vRObKZHCmRpyzDw8P1+bNmyVJWVlZat26tWOdYRga\nN26c2rRpo8TERMfhfAAAUDE1Ys++d+/e2rp1q4YNGybDMJScnKxly5apWbNmstvt+uqrr1RcXKwt\nW7ZIkp5++ml17tzZzakBAPAMNaLsvby8lJiYWGZZq1atHK+/++676o4EAIBl1IjD+AAAoOpQ9gAA\nWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgc\nZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUP\nAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCA\nxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR\n9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYAAFgcZQ8AgMVR9gAAWBxlDwCAxVH2AABYHGUPAIDFUfYA\nAFhcjSh7u92uGTNmaOjQoYqJidHPP/9cZv3777+vAQMGaMiQIdq0aZObUgIA4Jl83B1AktLT01Vc\nXKzU1FRlZWVp3rx5Wrp0qSQpOztbK1as0Jo1a3T27Fk9+OCDuv322+Xn5+fWzI/O2+iSef4aF+GS\neQAAcKZGlP3OnTvVo0cPSVKnTp20e/dux7pvv/1WnTt3lp+fn/z8/NSsWTPt379fHTt2dFfcKudp\nHySqMq8r5uYDFYA/uhpR9gUFBQoODna89/b2VklJiXx8fFRQUKCQkBDHuqCgIBUUFJjOGRZ2/mc+\nWviA6wMJrXmUAAAMH0lEQVRX4bxVObenzVvVc1fEhf+PPIWn5ZXIXB08La9EZlerEefsg4ODVVhY\n6Hhvt9vl4+NT7rrCwsIy5Q8AAC6vRpR9eHi4Nm/eLEnKyspS69atHes6duyonTt36uzZs8rPz9cP\nP/xQZj0AALg8m2EYhrtD2O12zZo1S99//70Mw1BycrI2b96sZs2aKTIyUu+//75SU1NlGIbGjBmj\ne++9192RAQDwGDWi7AEAQNWpEYfxAQBA1aHsAQCwuD982e/atUsxMTHujlEhpaWlio+P17BhwzRi\nxAj98ssv7o5UIf3791dMTIxiYmIUHx/v7jim0tLSHHmHDBmiDh066PTp0+6O5VRxcbEmT56sIUOG\n6NFHH9WhQ4fcHemyyvs7l5ycrFWrVrkp0eX9Z96DBw9q+PDhGjZsmGbNmqXS0lI3pyvff2bes2eP\nevTo4fh/esOGDW5OV77/zDxp0iRH3oiICE2aNMnN6S518e940KBBevDBB5WUlCS73e7mdJeqEd+z\nd5c33nhD69evV0BAgLujVMiFWwW/9957yszM1Ny5cx13Gqypzp49K0lasWKFm5NU3IABAzRgwABJ\n0nPPPaeBAweqdu3abk7l3Pvvv6/AwEC9//77+vHHH5WUlKS33nrL3bHKdfHfuVOnTmnKlCk6dOiQ\nRo8e7eZ0l7o476JFi/T000/rlltuUVxcnDZu3KjevXu7OWVZF2feu3evHnnkET366KNuTubcxZlf\nfPFFSVJeXp5GjRpV43YSLs47ffp0TZs2TeHh4XrxxRf10Ucf6YEH3HuPkIv9offsmzVrpiVLlrg7\nRoXdfffdSkpKkiQdO3ZM9evXd3Mic/v371dRUZEeffRRjRo1SllZWe6OVGHfffedDh48qKFDh7o7\nymUdPHhQd955pySpZcuW+uGHH9ycyLmL/84VFhZqwoQJNe4fxgsuzrtkyRLdcsstKi4uVnZ2tkJD\nQ92YrnwXZ969e7c+//xzjRgxQgkJCRW6KVl1c/Zv8ZIlSzRy5Eg1aNDADamcuzjviRMnFB4eLun8\nV8l37tzprmhO/aHL/t5773XcvMdT+Pj4aOrUqUpKSvKIryD6+/tr9OjReuutt/Tcc8/pmWeeUUlJ\nibtjVchrr72m//qv/3J3DFM33nijNm3aJMMwlJWVpRMnTtTYw8sX/51r2rSpbr75ZjcmuryL83p7\ne+vo0aPq27evcnJy1KJFCzemK9/FmTt27KgpU6Zo5cqVatq0qV555RU3pitfef8W//7779q2bZvj\nKFtNUt7/x1999ZWk80dgi4qK3BXNqT902Xuq559/Xv/4xz80ffp0nTlzxt1xLqtFixbq16+fbDab\nWrRooTp16ig7O9vdsUydPn1aP/74o2699VZ3RzE1cOBABQcHa9SoUdq0aZPatWsnb29vd8eyrMaN\nG+uzzz7T8OHDNW/ePHfHMdW7d2+1b9/e8Xrv3r1uTlQxn376qfr27esR/y8nJyfrtddeU2xsrEJD\nQ1W3bl13R7oEZe9B1q1bp9dee02SFBAQIJvNVuP/IqxevdrxD+KJEydUUFCgsLAwN6cyt2PHDt12\n223ujlEh3333nf70pz9pxYoVuvvuu9W0aVN3R7KssWPHOi6ADAoKkpdXzf8ndPTo0fr2228lSdu2\nbVO7du3cnKhitm3b5jg9VdN98cUXSk5O1uuvv67c3Fzdfvvt7o50Cc86hv0Hd8899yg+Pl4jRoxQ\nSUmJEhISVKtWLXfHuqxBgwYpPj5ew4cPl81mU3JyskecOvnpp5/UpEkTd8eokObNm2vx4sX661//\nqpCQEM2ZM8fdkSwrNjZWcXFx8vX1VUBAgGbPnu3uSKZmzZqlpKQk+fr6qn79+o7rfmq6n376yWM+\nuDZv3lyxsbEKCAhQt27d1LNnT3dHugR30AMAwOJq/jEoAABQKZQ9AAAWR9kDAGBxlD0AABZH2QMA\nYHGUPVADHTlyRG3atNGMGTPKLN+3b5/atGmjtLQ0STK9zWx8fLyOHj1aZTmt4PDhw0pISHB3DKBK\nUfZADVWnTh1t2bKlzK1vN2zYoHr16jnef/jhh5edIzMzU3y79vKOHTumw4cPuzsGUKVq/t1NgD+o\noKAgtW3bVjt27HDctnfr1q1l7uzXpk0bHThwQBMmTFCrVq301FNP6dVXX9W+ffvUrl07nTx5UrGx\nsVq5cqUGDhyo5cuXq0mTJsrMzNTLL7+sFStWKCYmRtdcc43+93//Vy+99JKys7OVkpKikpISNWnS\nRElJSZfc/vPLL7/UvHnzZBiGGjVqpIULFyowMFDJycnatm2bbDab+vXrp9jYWGVmZurVV1+Vr6+v\njhw5ooiICAUGBio9PV2S9Prrr6t+/frq3r27evfurW+++UZBQUFasGCBmjRpoqysLM2ZM0dnz55V\n3bp1lZiYqObNmysmJkYdOnTQzp07derUKU2bNk09e/bUb7/9phkzZujXX3+VzWbT5MmTddttt2nJ\nkiU6ceKEfv75Zx09elSDBw/WE088odmzZ+vIkSN67rnnNHPmzOr7DwxUJwNAjXP48GGjV69exvr1\n641Zs2YZhmEYu3btMuLi4oypU6caa9asMQzDMFq3bm0YhmH89ttvRo8ePYxPP/3U6NWrl5GTk2MY\nhmH06tXLOHz48CWvt2/fbowcOdIwDMMYOXKkkZKSYhiGYfz+++9Gv379jNzcXMMwDGPVqlVGQkJC\nmWxnz541unfvbuzdu9cwDMNYsGCBsXz5cuOdd94xxo0bZ5SUlBhnzpwxBg4caGzatMnYvn270blz\nZ+PYsWPGmTNnjE6dOhmrVq0yDMMw4uLijLffftvxZ0lLSzMMwzCWL19ujBkzxjh79qzRq1cvY9eu\nXYZhGMaGDRuMAQMGOHLPnj3bMAzDyMjIMKKjow3DMIynnnrKSE9PNwzDME6cOGFERkYa+fn5RkpK\nijFo0CDj7Nmzxm+//WZ06tTJyMvLK/O7AKyKPXugBouIiNBLL70ku92uTz75RH369NGGDRsuGRca\nGqq4uDhNnDhRr732murUqXNF2+nYsaMkadeuXTp+/LhGjRolSbLb7brmmmvKjD1w4IAaNmyoG2+8\nUZI0efJkSdLEiRMVHR0tb29vBQQEKCoqStu2bVNERIRat26t6667TpJUt25dde/eXZLUqFEjnT59\nWpJUq1Yt9e/fX5IUHR2tRYsW6dChQ6pdu7YjX58+fTRjxgzl5+dLknr06CFJuuGGG5Sbmyvp/FGH\nH3/8USkpKZKkkpISx2H6bt26yc/PT6GhoapTp45jHsDqKHugBrtwKH/nzp3avn27Jk+eXG7ZS9KP\nP/6o0NBQ7d69W3fddVe5Y4z/f/7+4scM+/v7S5JKS0sVHh6uV199VZJ09uxZFRYWlhnr6+srm83m\neJ+fn6/CwkLZ7fZLtnXhegNfX98y68p7gJOXl5djXrvdLm9v70vmvHjeC8+G+M88drtdf/vb3xwf\neE6ePKnQ0FClp6eXeZaEzWbjegb8YXCBHlDD9enTRwsXLlT79u2dPkRo3759Wrt2rdLS0pSWlqb9\n+/dLOl+qF4qxbt26OnjwoCQpIyOj3HluvvlmZWVl6aeffpIk/eUvf9H8+fPLjGnRooV+//13x1xv\nvvmmVq1apVtvvVXr1q1TaWmpioqK9NFHH6lbt24V/nMWFRVp48aNkqS0tDTdeeedatmypXJzcx1P\nbduwYYMaNWp02SMXt956q959911J0sGDBxUVFXXZ54t7e3tf8uEHsBrKHqjhevXqpX379unPf/5z\nuevPnTunuLg4xcfH69prr9WUKVM0depUnTt3TnfddZdiY2N1+PBhTZw4UXPmzNHAgQMVEhJS7lxh\nYWFKTk7WU089paioKO3Zs0dTp04tM6ZWrVp64YUXNGXKFEVFRengwYOKjY3V0KFDde211+qBBx5Q\n//791atXL/Xu3fuK/qyffvqpoqKitGXLFiUkJMjPz08vvviikpKS1LdvX61cuVIvvvjiZeeYNm2a\ndu3apaioKE2aNEnz589XcHCw0/GtWrVSfn6+nn322SvKCngSnnoHoEa48M0CAK7Hnj0AABbHnj0A\nABbHnj0AABZH2QMAYHGUPQAAFkfZAwBgcZQ9AAAWR9kDAGBx/w+YPGkx4XpWFQAAAABJRU5ErkJg\ngg==\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": { "scrolled": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 5000/5000 [00:28<00:00, 178.45it/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": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFcCAYAAAD/HgalAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4HNW5/z8zs7vaXWmlVXdRdcUF2+BCN8aUmGKMcylp\n9As8gCmBawIkFIcECJB24RLKj9ASCAESTCcGF3ox4N6LJNuSbHVptX1mfn+sdy1pZ4tkdZ3P8+jB\n7J45c87s7nznvOctkq7rOgKBQCAQCHodua8HIBAIBALBUEWIsEAgEAgEfYQQYYFAIBAI+gghwgKB\nQCAQ9BFChAUCgUAg6COECAsEAoFA0EcIERYIepG9e/dy1FFHGb731Vdfcc4550T+PWXKFBYsWMCC\nBQuYP38+l1xyCZ9//nnUcS+++CLjx49nzZo17V5/9NFHOfbYYyN9nHvuucydO5cHHngAo8jEju3P\nPvtsbrnlFsrKyhLOq6WlhUsuuSSJKyAQCNpi6usBCAQCY4qKili6dGnk/7ds2cKVV17J448/ztSp\nUyOv/+Mf/2D+/Pk8//zzTJs2rV0fZ511FnfffXfk/5uamjj33HM58cQTOemkk6LO2bH9G2+8waWX\nXso777xDWlpazLE2NTWxfv36Ls1TIBjKiJWwQDBAOOKII7j44ot57rnnIq999dVXNDU1sXjxYj76\n6COqqqri9lFbW4vX6yUjIyOpc5533nmMHj2at956C4DXXnuNCy64gPPOO49TTjmFl156CYA77rgD\nr9fLggULUFU1ZjuBQNAesRIWCAYQRxxxREQQAV566SXmz59Pfn4+xx57LH/7299YvHhx5P13332X\nb7/9Fo/HQ1NTExMnTmTJkiVMmTIl6XOOHz+ebdu20drayquvvspTTz1FZmYma9as4fLLL+cnP/kJ\nDzzwAPPnz2fp0qVx2wkEgvYIERYIBhCSJGG1WgGoqanho48+4vXXXwdCq9Z7772X66+/HrvdDhwy\nL/v9fu677z527NjB3Llzu3TO1NRUnnjiCVatWkVZWRlbtmzB7XZHtU+2nUAgEOZogWBAsX79esaN\nGwfAP//5TwCuvfZa5s6dy0MPPYTL5eLf//531HEWi4W77roLl8vFQw891Olzjh8/nurqas477zz2\n7dvH9OnTufnmmw3bJ9tOIBAIERYIBgzr1q3j5Zdf5tJLL0VVVV599VWWLFnC8uXLWb58OStXruSa\na67hhRdeMPR+tlgs3HPPPbz00kts2rQpqXO++uqr7N27lzPPPJMNGzaQlZXFddddx4knnsiKFSsA\nUFUVk8mEqqrouh63nUAgaI8wRwsEvYzb7Y4KU/rHP/4R1a6iooIFCxYAIMsyaWlpPPLIIxxxxBF8\n+OGHaJrG/Pnz2x1z2WWX8cILL7Bq1SrDc8+YMYP58+fz61//mpdffhlJktq9H95DliQJTdMoLS3l\nhRdeICUlhRNOOIHXXnuNefPmIUkSs2bNIisri/LycoqLi5kyZQpnn302zz77LPn5+YbtRo0adTiX\nTiAYdEiilKFAIBAIBH2DMEcLBAKBQNBHCBEWCAQCgaCPECIsEAgEAkEfIURYIBAIBII+QoiwQCAQ\nCAR9RK+HKNXUtPT2KXuMzEw7DQ2DPxPQUJinmOPgYCjMEYbGPAfbHHNzHYavi5XwYWAyKX09hF5h\nKMxTzHFwMBTmCENjnkNhjiBEWCAQCASCPkOIsEAgEAgEfYQQYYFAIBAI+gghwgKBQCAQ9BFChAUC\ngUAg6COECAsEAoFA0EcIERYIBAKBoI8QIiwQCARJ8N13qznnnNNZtOhqbrjhGq6++jJeey26DnQi\nVq1aQW1tTVJtt2/fyrPPPt3pc/QmO3fuYM2a75Jq6/P5OP/8+YkbxuDdd9/iL395tMvH90d6PWOW\nQCAQDFSmT5/BkiUPAOD3+/nJT/6LH/zgbBwO42xIRrz66suUlNxJTk5uwrZjx45n7NjxXR5vb7By\n5UdkZ2czbdrRfT2UAYkQYYFAMOC4995f8dZbb0S9LssSmqZ3qc/588/j3nt/k3R7t9uNLMsoisK2\nbVv44x8fRlEULBYLt932KzIzM7n77ttpbW3F5/Ny7bU34vV62bFjG7/5zd08/vgzLF36OsuWfYAk\nSZx66hlccMGP+O1v76WpqYnm5iZ+/OOLWb78PyxZ8gD/+c97/POfL2M2mxk7djQ33ngb//nPe7zz\nzptomsaVV17DjBmzIuN74onHWLv2OzRN56KLfsrs2XNYtOhqLr/8KsaOHceNN17L73//v9x3390U\nF5dQXl4GwJIl95OdnRN1/Ny5p7Fx4wb+/OdH0HWd3Nw8fv7zxbz33tuYTGbGjTsCn8/HU089jqIo\njBgxkttu+yV+v59f//pXtLS0MHJkQdR1/PTTVXz88UruvPMeAC6//Cf84Q+P8cEHS3nnnfcIBoOk\npaXx298+HDmmqqqSe+65k6eeeg6Aq6++jCVL7sfhSOfBB39NU1MTADffvJjRo8d09qvQqwgRFggE\ngiT59tvVLFp0NbIsYzKZ+PnPF2O32/nd737L7bf/irFjx/PJJyt57LE/cMUV11BfX8ef/vQ4DQ0N\n7NlTzvHHn8iYMeNYvPhO9u7dw0cfLePxx/8fkiRx883XccwxxwKhFfdFF/2U775bDUBTUyPPPPMk\nzz77d+z2VJ5++lGWLn0dm82Ow+HgwQf/0G6cX3zxGVVV+/jLX/6Kz+fjmmsuZ+bMY7jnnt9w2203\nk52dw/XX30R+/jAAJk+ewuLFd/Kvf73Kiy8+yzHHHG94/EMP/ZYlS+6npKSUf/3rVerr6znzzHPI\nzs5mwoRJ/PjH/8Vf/vL/yMzM4umn/8K7775FIOCntHQ011xzPRs3bojMKcxxx53I44//Lx6Ph7Ky\nXYwcWUBGhpPGxkb+9KfHkWWZW25ZxObNGxN+Pi+88FemT5/FwoXns2dPBfffv4S//OWZ7vjoe4wB\nL8KaptHQ0EBWVhaSJEW9r+s6mqahKEMjD6lAMBS4997fGK5ac3MdPVokpq05ui21tTURs/HUqUfz\nxBOPMWrUaH74wwu5995fEgwGOf/8H7U7ZteunezfX81NN10LQEtLC3v37gWgqKi4XdvKyn2Ulo7C\nbk8FYObMmXz44QomTpwc1TbU9w62bt3CokVXAxAMBqmurmLs2HFMmTKNDRvWc+yxx7eZ10wAjjxy\nCp9+uorc3DzD4xsa6ikpKQXghz+8AAitZAEaGxuoq6vlrrtuB0L7v7NmHUtTUyPHHHMcAJMmTcZk\nai87iqIwZ86prFq1nA0b1jN//kJkWcZsNnPvvb/EZrNx4MABgsFgrI8FXdcj8/7uu9V89NF/Ite0\nvzPgRdjr9fLOO0ux2exkZGTgdGZiNltwuVpoaXHR2toMSBQWFjFx4iTS0zOi+ggGgzQ1NVFTcwC3\nO2Q6GjmykMLCIkNhFwgEgrbk5OSyY8d2xowZy5o131FYWMTOnTtwu1t5+OE/U1tby7XXXsEJJ5yE\nLMtomkZRUTElJaP4/e//F0mSeOWVvzNq1BhWrPgQSWrvMzt8+EjKynbj8Xiw2Wx8/fXXFBYWAUS1\nBSguLuGoo2bwi1/8Ek3TeO65/8fIkSPZsGE9u3btZNq0o3j55b/xk59cDMDWrZvJy8tn3bq1lJaO\ninl8Tk4Oe/ZUUFhYxN/+9hyFhcUH56OTkeEkLy+PBx/8A2lpaXz66SpsNjs7d+5gw4b1nHTSHLZt\n22Iopuecs4CHH76fpqZGbrnlNnbs2M6HH37I44//Fa/Xy5VX/qxde4vFQkNDA6qq4na7qaqqjMz7\njDMmcsYZ82hoqDfcsuhvDHgRBjCZzCiKgsvlwuVytXtPUUJTrKzcR1lZGXl5eQwfPpzW1laam1tw\nuZpxu92AjtlsiYjunj17WLPmewoKCpgwYRI2m61bx6zruhB4gWCQ8Itf/JI//vEhdF1HURRuv/0u\ncnJyefbZp3j//XcwmcxceeU1QMj0+5vf3MMf//gYM2bM5LrrrsTvDzBhwiRyc42dtZxOJ1dccQ03\n3ngNkiQzenQpl156TWTF15ETTpjN999/y3XX/Tcej5vZs09B03QefPA+7r//YfLzh3H11Zdx9NHT\nAXj33bd55ZWXsFqt3HXXr0lPz4g63m5PZfHiO3nggV8jyzLZ2dlceOFPMJvNPP74nykpKeWmm/6H\nxYtvQtd17PZU7rprCVOnHsUDDyzh2muvpLi4BLPZHDXeESNGAnDSSXOQZZmCgkJsNhtXXnkxFouZ\n7Oycdh7l2dk5zJw5i6uuuoSRIwspKCgE4JJLruDBB+/jzTf/hdvdyhVXXN31D7WXkPTwOr6X6G5T\nkdvt5o03XsdqtXZrv23x+wPY7TbS0hykpqbicKQjyzKZmXbq613oesghRFEUTCYzKSkpmM0WLBYz\nimJCURRaW1vZt28vDQ311NfXIcsKpaWjDM0z/Y2eNvH1B8QcBwdDYY7QvfNctOhqFi++k+Likm7p\nr7sYbJ9lrHrC/fvu30+wWMwEg0EaGxtobGyIvJ6ebqO52dOuraZp7f4AJElCkiQsFkubdirbt29l\nx45tFBYWM3HiRFJT0xKujlVVxe/3EwwGSEmxtutTIBAIBAMLIcLdjCzLyHJyOVDC7fburWDHjm0o\nioLVasVut2M2p6CqAQKBIMFg+L9+gkENXdcOmrNlzGaFlBQrKSkpAKiqfvB9DbPZjNVqw2azYbXa\n0DQNr9eLzxf6CwQCBINBVFVFVYOYzSlkZmaRmZlJQUEhmZmZwmQuEAxyHnvsqb4ewpBGiHA/IWxO\nV1X1oEdftBkmZNqOPjYYDBo6O/j9flpbW5M6v6KY0DSVuroa6upq2LRpAyaTQnp6BkVFw1FVExkZ\nGVgsKZjNpsi+TnNzMy5XCz6fH03TyM7OZvjwEYPGG90XUGly+chISyHFPDjmJBAI+g9ChAWGhM3c\nLpeLvXv30tzsIRgMRkK+NE1DksL74KbIinnr1s3oOqSnO8jIcKLrOoFAEE0LPShYLFasVuvB1bmV\nnJxcnE5nlGhrmobL1YLH40FV1TYmfj1iCdB1DVlWsNlsB60HFlJSUpK2RMRD1TReWb6D77fVUN/s\nIys9haPG5XLR3DEo3dC/QCAQgBBhQSdIxoEs3Mbj8eDxeAxaHPJe13Udv9+PJMkHHd/SIiEHXq8X\nVVVRFOXgnjpAaG9d1cCnKlhNGrKkR8zpkhRKoJCVlUV2dg4jRxaSk5MDHNqrDwu6qgYPmvhD2Xg0\nLTUyLr/fz/97cw2rdxwaa12zjw9X7yUYDHLJvIldun4CgUDQESHCgj5DkqTIXnbI8a0x8p7ZbI4K\nZdB02FRjZ7/LgicoYzNp5Kf5mZjrbte2qamJpqYmtm7dDEiATshHTicUC6AfDCWRAQld18jISMXr\n1bBYzDQ1u9hYmQdEm5+/XL8Xm2sj6Y7Ug97yaVitKVitIe95u92OyWSK+AZ0dlWu63okFlTsx8dG\n13VaW11Rr1ut4HJ1zaM2GcdIgaC7ESIsGDBsqrFT1ngoXtsTVCL/PznPHdXebE7ec9xiseD1evD7\n/WiyFU/QWDy9qoI3KKO0trbbb9c0jWAwgKbpkeQJsixFxF6WlYN7+jImkxmbzY7dbsNms2MyKZEH\nh+bmJgKBABaLmbS0dNLT03E40hk1ajSpqamGY1JVFV3XuyT6A5XWVhevvvqPyENcGIfDRkuLkQUm\nPj6fjwsu+BFpabELMfj9fu6/fwmVlftITU3lllt+QWFhEatWLef//u/P5OXlA3DlldcwfvwEbr/9\nFnw+H4sX38mYMWNZu3YN69ev4Wc/u6zT40tEc3MTX375BWecMS+p9n/+8++56KKfMmzYsG4fSyL+\n8pdHKS4uYezYcXz66cdcfvlVnHvuD3jzzQ/atXv00UexWtM477zzE/Z5552Luf/+h2O+v3Tpvzj7\n7HO7PRy0O8K7hAgLBgSqBvtdxqK632VhQo4bJYH+qBp4gzJWkxa3rdWkYTNpeILRK2GbScNq0qJe\nl2UZiyUl6vW2aJqKpqkEAgE8Hjf19dFtFEWJ7I+73a243a1UVVWyYcN6nM4M8vOHUVBQRGNjKN68\noaGB5ubmg+If7kNGUcykpFgOWhRSMJkUsrMdNDd7kWWZ3Nw8SkpKB/TKLyUlZIFoi81mIxDomfO9\n9da/sdnsPPXUc1RUlPHHPz7EH/7wGFu3buG6625kzpxTI21XrVrOiSfOZtq0o3n77aXcdNOtvPrq\ny9x11697ZGw7dmzns89WJS3CN910a4+MozN0Z4WoeAIM8OKLzzJv3tn9MidD/xuRQGCANyjHXJ16\ngjLeoEyqJVocIb4ZWzbQIEWG/DR/u1V3mPw0f0Kx725CZnvLwQT3u9m2bQsmkzmy6jXKQAQ6Pp8P\nn89HeB/e52uJxLWXle1mw4b1jB07lnHjjhgyK+jDYffu3ZF8y0VFJZSV7QZg69YtbN++lX/+82Um\nTJjEtdfegM1mP+gX4cVms7Fs2fvMnn1K1Mo9zKJFVxtWMnr00T+ybt0aAE4/fR4XXvhjVq1aziuv\n/A1dlxg+fAS/+tUSXnjhr+zYsZ2lS//Fsccez0MP3Y/f78NiSeG22+5E0zR+8Yufk56ewXHHncAX\nX3zG4sV3kpWVzX333UVrayuqqnLVVdcyffpMLr74QgoLizGbzSxZcn9knB2rKN1zz33ceuuNOJ2Z\ntLS08PDDf+L3v3+QvXv3oGkaV111LUcfPYOVKz/i+eefwenMJBAIUFxcwnffrWbp0tdZsuQB/H4/\n99xzBwcO7Gf06LHceuvt7a6PUUWntoRX0osWXc3YsePZtWsnbreL++77HatXf0V9fR333nsnDzzw\ne8O+Fi26OjIHu93OhRf+mKOOms7mzRt5/vlnuOuuX/Pgg7/B5WqhqamR+fMXsnBh4hV6MggRFgwI\nurI6DdNZMzbAxNzQ60bC3dckWnEng8lkIhDws379OrZs2Ux+/jAsFgsWSwoWi4X8/HyczsxuGO3g\nYezYcXz++SfMnj2HjRs3UFtbg6qqzJw5i5NOmsOIESN5+OH7Wbr0dRYuvIDPP/+EN954jauvvo7/\n+78/c8UVV/PQQ79l5MgCfvrTS6P671jJaObMY6mqquSpp55DVVWuvfZKpk+fybJlH3DZZZcxa9Zs\n3nvvbVpbW7nkkitYuvR1Fiz4IXfffQfnn38Rxx13AqtXf80TTzzG1VdfR319Hc888zfMZjNffPEZ\nAM8//wwzZhzDhRf+mJqaA1x33X/zyitv4PF4uOyyKxk37oh2Y+xYRamsrAwIPSCcfPIp/Pvfr5GR\n4eSOO+6mqamR66+/mr/97Z88/vj/8vTTz5OensHixTdFzd3v93HttTcybNhw7rrrdj777OPIe7Eq\nQsWq4TxhwiRuuulWnnzy/1i27AMuvvgynnvuGe699/6YfbWdwxdffMZ7773NUUdN591332b+/IXs\n3buX0047g5NPnkttbQ2LFl0tRFjQdZI1y/Ynuro67aoZW5ZCAj0hxx11rfrT9TvcsSiKgqZpkQT4\nYb7//lscjnRGjBjBmDFjcTjSu2nEA5ezzz6X8vLd3HDDNRx55FTGjz8CRVE4++wFEUE46aSTWbly\nObIsc/PNi4GQKfSCC37E888/w89/fht//euTVFSUR1U/6ljJKC8vn6lTpyFJEiaTiUmTjqSsbBc3\n3PBzXnvt77zwwt8oLi5h9uw57frZtWsHL774LH//+/PAoYiF4cNHRFlNyst3R0zYubl52O2pkayA\nRUUlUdfAqIpSqG1oLjt37mDduu/ZtGkDAKoapL6+jtTUVDIynEDoYaMjeXnDGDZseGT+FRXlhI0z\nsSpCxRLhceNCJu78/Hzq6uqiro1RX23ncMwxx/H443+mubmJdeu+5+ab/4eGhnr++c+XWLVqBXZ7\natyKTp1FiPAQorNm2Y70tfh0ZXV6OGZsCIl/+P3DvX5tOdxr2Z1jMcJiseDzedm9exdbt24hPT2D\nzEwnTmcmI0cWGFYjG+xs2bKJKVOmceONt7JlyyYqK/ei6zqXXvojnnjir+Tl5bN69TeMHz8hckxD\nQz179lRw8cWX8/LLLyLLMpIk4fVGO49FVzIq5d133+Sii35KMBhkw4Z1nHnmObz55r+54YYb0DQL\nDz30Wz7+eCXDh49A00JlAIqKSvjxj3/GkUdOpby8jO+//xaIVW2plLVr1zBu3BHU1BygpaU58tka\n+QsYVVGCQ9n/iotLyMvL45JLrsDn8/L883/F4UjH5WqloaGBzMxMtmzZFHFiC1NTs5/a2lpycnJY\nt24NZ5+9gPLy7ZE+jSo6xcJo3JIko+t63L7Cc5BlmVNOOY1HHnmQk06ag6IovPzyi0yePIWFC8/n\nu+9W88UXn8Y8f2cRItwL9KR4dabvrphlIXTD33Cg8zf87p63LMGEHDdF6V6QINWcuN/DMWN3pKvX\nry1G4pmb6qfU6cWexHy6cyzJEvYcr6ryUFlZyZo135OWlsbMmcf2iXdtmNB+d3vMZgwFrit9daSg\noIinn36Cl1/+G2lpDu644y4kSeL22+/il79cTEqKlZKSUs49d2HkmOeff4ZLLrkCgIULL+DWWxeR\nnz+MMWPGRfXfsZJRRoaT77//lmuuuZxAIMDcuacxfnxILC+//HJSU0MhcccffyJ+v59du3bwz3++\nxPXX38Tvf/8gfr8fn8/LTTf9T8w5XXLJ5TzwwK9ZufIjfD4ft932y7jOS0ZVlF599eXI+wsW/JDf\n/e43LFp0Na2tLhYuvACz2cydd97NrbcuwuHIMOw/I8PJn/70MDU1B5g8eQrHHXdCRISNKkKF6yon\ny9Sp0/if/7mRRx99Mqm+zj77XC68cAH/+Me/I2N45JEH+M9/3iMjIwNFUfD7/Z0aQyxEFaXDwKiA\nQ1t6crXS2b5VDVaWOWOIkcqcksaYArCtMZ1tB6Kdf0qcHsMbfk/M+3D63HDAbmjGbjv+RJ9lvOtn\nVVROKY19/ZIZC+iderjpymeZaI6dJRgMMmbMWI46anqve1nHihPOyXFQWzvw4oQ7G+oy2CoMGTHY\n5iiqKPUBPbla6WzfXTXLqhrsazTOmRxrT7Un5n04fSYyY6satHhDmbhiCWm86+dVZdbvtzNlWGLx\njLU/DVLSczpcE3t3YTKZ2LlzB/v37+fEE0/q1X1jSZIMY3odDgdeb68NQyA4bHpdhDMzQ8kJugu3\nW8HhCOUi7gvS043P6w3AfpexF2uNOwV7qk5XL4Nfhb3Nxiv/A63GfdtVsO/TcQeiVSLVopObmWI4\nnhavRKvfWFk8QRmT1Y7DesiYElRD8zOiq/Pujj6Pz9AJqj48AQmbOdRe0218t8fMvkaFVr9EqiWF\nkU6VowsDUWJqVyG1Uo9xLST2ttiwWk3MKj4UpBpUaXe+Fq8UUzw7M6d4Y4n3WULs7+vhofPFFysZ\nPXo0M2bM6PNYzFgrjv7MK6+8nLhRBwbiPDvLUJhjr/9aGhq6d7/K7XbT0uLtdIB+rP3KePuYHd8z\nrCd80Gxa1WLBpxqLV6tfoqbBF1mtdHbv9PuqVIKacd/ugMTnO6R2q7LwmHwxHPpy7T7crcZmSlWD\nVEuK4Q3fqmgEvW6a22yNtPplWv3GN/qO806W7u7TfXC8IdPwITN7q19i2wEZvz9ouBLNtUsxxwGw\no8ZEwB9kQq6bLbXRpvPx2W5sJouhGbmzc4o1lnifZXebozuyZs1G1q7dxBFHTGTChIl9YtodbCbM\nWAyFeQ62OQpz9EFi7S0ekWN84wybLI2OOdYRvZ3e0WxqhIxOUIOAClvrEu9zthVpgDq3UXKGMKFV\nmUk5ZNKMNSZF0ijM8MX1LlZkGOlU2XYg+ukgoElsrrW3G293OkKF6Yk+uxK6NDHXTUCV2NeSQign\ndUckypts1HtMtPgPfUZtzcyxwqzaksyc+mMcc1h0N25cx+7dO5k79/Q+s1AJBAOFISfCsfYWa1pN\ntAaMb5yA4TGWPQHGOQ+tLOLv+R1CQ+aTCieKpKPqh+704X4DqsSU/FYkKVr8s2wBvGri5XJYSML/\nNsIs60zISezYNGVEALcnSJUrpd14VV2O2sPsSjxvIktAd2ew0nRYtz+10/uqsgRT8lup95jjrmZb\n/MY/q/0uC7OLGyP/Dp0/+uIbzanjNYoXx9zXKIoJn8/HsmUfcMYZ8/rEaVIgGCgMCREO38DMshZT\nkFoDxpeiojEFi8nYgXxXrUKpA8K13uM5zEQjoerG6revJYV6jxmTrEWtqPa1KCiSFvPYQ21DQhL+\ntxFeNbl0jzVuY3N0mKoWC+Oy3FgOXsJkV2md8XiO1ef4bDet/s6J0KYaO/taYgtDvJVovAeCRHiC\nMk1eExNyQuLpDsiUNVo50Br7OiW6Rm3jmPsbwWAgIsSx0jUKBEOdQS3CHW9gKYqGL+Yq0lhkNGS8\nQWMRDmoyGw/YmTbcjaqFxD6W2bRzSAf76Ho/bYUk1phSFA2znHy6x1j4VJlV5U6GOw6JQzKrtFhW\nibAloO0xHfu0KBpb6+x8XO5MOkQr0YNYGJOsEW87c2KuG02HiiYrsb43sfhyX3q7sR6Z70bVYl+n\n3owH7gkCAT/Llr3P6acLIRYIjBiUIhy+4e5qsFLedOgG5lO7zys7TE2rhe+qoMETEnpF6vmwa1WX\nKHB4qPMkZ9LMS/W3uw5hfGrILN5WvDojVoeQ8KnR4hBvlRbPdB+2BBiJarjPjvG2scSpcw9iIVr8\nZjbV2GOKXMgs7UYCw+vqsATbWTAOIRmONdZ16o7KUf0Bv9/Phx9+wLx5Z0cqRAkEghCDRoR9Ph8f\nf/Ix5e5cnCOnkJLqRNdB7pbffOzVjk+TqWw5dCMOm4nDJuOO+77dgVXRODLfDcQ3abYVIGj7cCBF\n/hsWBF0HSaJTYmVEsuIQ33TfPma244q6M+LUcSWZ7INYMvOYlOdud83C2a9KMryUNx36TMJz6uw5\nOhsP3NdpRePh9Xr57LNPovIcCwRDnQEvws3NzSxb9j6rVq2gaPqFjJp+qKZnPENhwO/BZLZ2QxiF\n8fE+dxPfdY8DAAAgAElEQVSenW+RZtExDTuGoHUkmFMN87fqqg9JthDXBtoGvyqx8eBKzZGixTRp\nxs7OFM3e5vZOV/HFKizo0eNNNllEPI/ntuxpSonaDy3O8MYVpwaPiUxbKB4r+dV85+fR1kTuDsjs\nbrCy32WhoskaEeR8u59vqoyTWMQ6R1trRLxrtKvByqS82N778TJv9bZgS5JEdXUVGzasZ/LkI3v+\nhALBAGFAi/CuXTv5wQ9OoampkbR0J6VTTkn6WLMlBVQ/KD2zTyVbHHzz1ee4m6qBfyObLNjScig5\n6mzyR03HmpaD11VL9Y6v2fr5S0w65SqKJp8a1U/j/p2kOodhTjmU31RDpqLJxvaKeoK7/k1OTg4O\nhwNN01BVFVVV0VFozT0j6W3lRI5ebRnp8FLtSjE8JtmQoWQdnFRdxnMwvjm8Otb0+Hvv4X3XLFsg\nzmpbxyJr+DVjc35nQp8kCb6rSotyoqtoslHvNhn2b3SOgAobD9gj2ww2k4ZJ1jD+EEPhUOHntlim\n+bZWBOj5wg/xMJlMbN68kezsbIYPH9GzJxMIBggDWoTNZjMTJkwkOzuH4+fM44uqtE4cLRsKsElS\nCerGN2agg2dy7LuWzazz8xuuw93agtvtxmaz4XA4SEtzYLNLtPrq0INuAkWj8Rz3PxyoqaG6bj2q\nbSSKNQN/awOu/Zvx7l9D2vRLDM+RkjGSSkpZ8erT7V6XJJkjT7uOwnxb0m5DYXN0ItSAh70Vu5Ey\nJxq+b/FXs3nTdoLBAHZ7KgUFhdjtdsO2YS/g6hbLwbCr5EZb02qJuc/ddt81nie5VdGYPqKZ3Y3W\ndtsJYToT+rThgD3GHjC4YnjdQ2ivXpEPCeOeppSokDVQSDMHDvYTPY/qFkvMz62jFaEwS8Pvl5IS\n7J5aHSuKwpdffs4ZZ5xJamrnkvALBIORQVPAwWyxsmK3E+9hOl/ZTCrZNj97DW7MEErWP31EM99V\npcc1pZY4PUnd1DqaBVWNqH3eEMZ3WhM+clpW4W5tQVEUFEXBnToZd0ppZ6adNJoaxOduxObIiXov\n4Gtl+V+vxWyx422tRwuGUlPl5ORQWFhETk4usqygKDKSbMZid5KblUZebh4HpNFUuozFOhqdk4sb\nKW+yxo23hfBDU/TFD4uzVdEwKxoaodSVXSnvuHy3M475Xo8xttAcHCnRTmYdSVHUg/vzxv2ESO4B\nJt71sCh6r62ObTY7P/jBmZHycd3FYMuyFIuhMM/BNsdBnzFLkSHbHmBfy+GJsCcoM8LhZ2+LcfiJ\nV5WxKLE9jhVJoyA9VBYtVOnG+KYWzyxY3mSNscqLJoiFCUfObJcCc2WZE+LUnHZYAgQ1OaHAG2E1\n6chp2YbvmSw2zvjvx5FMdgi2ojbtomz1q+ypKOP777/rYJKfge7IYUdLLZ9+/BVbP/01k+deQXbR\n0VjTsvG6ajGnpGO2Rguzv7We//3Db/F5WrHYHUw68y7MNuP6tkENzN696NZsVMmKTEhww0LkVRW8\nqsKY3ACFqa5OrwK9QblLDmw2k4bdrCWV4MWnhh3ljCo4hcKpkg2Li7XtYGT2h54Lg2ptdfHVV19w\n3HEn9Ej/AsFAYdCIMMDkvFaqXZYY3sixViTtsZk0nNZgzD1Hq6Kx66ADTluPY6uikWkLMCbLQ0WT\nNWH4TKz4T00PmVuTpeO+YnyvY50Ch5cpw9zoeigHdXVrrKQVxtdrWHqQmlbZ8NpIkgzmg1sCZgdK\nzlRO/dE4JuS0smavQo3XRpD2KR/tGfmMmn4uTmcmO7/8O007l5PmzMMiBcgYfRrmgllR59m75XNq\n9ldhsVhQZRtKSuxtCF3T8acMx9tcS92eL8gunIQ9Iz+q3da9bsqq3sNqMdPicuPyS7S6XATcDQwf\nlkdBQQEjRxaQl5d3sEC4hqbpyJqOVVHxqsY/pVjhSmFzd6s/cYIXm0mL+dBnVjSy7cEuJQ9JRE+G\nQcmyzJ49FeTk5DB27PjuP4FAMEAYVCJsVqAww2d4Q0qzBHHF2LdrS36aH4spttOQWdEMb4YmRaPR\na+aTitiOXsmkktzvsnRqZdVx79KiaAfDooz3QY/MPxgPrEOjN/71UKSDq2u9vZl2Uw1J3/T3uyzo\nOlR547cfMe54zps7mdSDhe01PeSktLdZaxf2NSLNw5kXHYfy4+NRtZCIfVNJZBXXEflgXKo9Ix97\nRj6xdl8kUyrLln9K8dQzKZx8FmkWO2lA0OehYuNHfPDM0+i6saPWxJOvYNT0c6Ne9zTuo3rdC2SN\nn4clawyYUpE1L3atFv3AbjbUKqiKHYWJqFLsDF7h694xJzWEYpqzbEFKnJ6IVcWqaAQ0KYbZ2fi7\nYURPl0U0mUysXbuGrKwcsrONrSsCwWBnUIkwxE5v2LFAQ3gvMGyW7Zgy0Kif3FR/zFVqMgKfTCrJ\neKZHBQ2TouNTY6eC3FpnjxmXPMxxSLC9QTlBDupDaTVLsgJMyGqOHNvx2hyKKTYOWapOIkzIq8p8\nXO6MzAuiE2GoemjPW5IkNhw49Fl2JkFKLCcmm1nnnIvvoNVc2O51s9XOqOnzGTt2LAc2vkVdXS2S\nJCFJMrIsARJS3WqayrKx503AZM3A72miZvdq1i57IiTcq79CNlmwpmbhba1HV4NMmH0Zw8Ycg9WR\ngxrwYjZ4dtOCXhzafoqsOqqaSlAz/rwOtFqYU9LYzgdhc63xPnNBus9QzI1IlFGtO5BlmU8/XcVZ\nZ83HbE48JoFgsDFoHLM6JolPtlRhonjJtu97gzIrypx0NlVhGJtJZXZxI5tr7TFTHtpMakzTYyJn\nr/B+sJGpWJE0ThvVEMlzHa9tR1ItGrOLGmIWFTDLGp9UGPcV36koNrEciOJdn0Oe0PHOZWxmH5Mb\nYG+9HNOxL0VRObm4EUVOvtSlqgZpbm6htdVFS0sLLlcLLpcLT9oUtIwjok+i+kE2ofqaOVC+ljX/\n+QtqIFShPi1zBCdf+iiSYfYZnWNHNpNpC0bGZORzUJilUZTazCcV8RzJ2vfbG05auq6TkeHktNPO\nOOy+BpszTyyGwjwH2xwHvWNWR2KlAuz4eqIE+G3fTzbBRCzy0/xsrbNTEcfpKnzD65iJqe2NMNZ4\nWwOx9xdVXcKvypgVLTKvZAsRtPolQ7Nk22sTq6+8VD+VLcYxxfGI1T7eytok6ZhkPa7AWJXQteyY\nYWx8PuyoiW0S9qkyH+7KRJKizfNhcYr6bikmMjMzyczMPDSvyMOPwdgsCrMKmkg1a0iTJ1I++VbW\nr19LZWUlLa1ufK0NWA280jVV5Yu9Dky6l5FOncn5XsNc27tb0g8KcGyfgRDRGdWg55y0JEmivr6O\n8vJyiouLe+QcAkF/ZdCKcE/Q+Qo6oZta24o/H5c7Y7YtzvAyMTfkNFXq9DIuy01ASxy3GV71VLfE\nNvsaJZ/oaFYOES1+qRa9y/VtgW5N2xkvnaYvhrm2LcMcfibnRWcYs5ttWBUtToibhIYU0amuilM8\nxzmvKqMcrIwEMqWloygtHRV5PxTKFH2crIR+xqpkp6IZtu3YxQljTQzPz+uQazu+uTdF1kAyzpbW\n07mqTSYT27ZtFiIsGHIIEe4kRmLTseRgmKIML6MzvZEbfSJP2BKnN2bYUjySqXZklHyi42qpY8GL\nMCOdasKbr1HlJDgYLtUFYjkQ5aeF9uVjea7HDtc59JAD0atWkxIS6M56GXdWnOJZU9o+KBltkyT7\n0JSSO5lPKg4gbynjlMnpODOzkkrf6bQF2R/D56GnnbQA6uvrqampITc3t8fOIRD0N4QId5K2YmOy\n2gl6Q6bjZFIBJroBlzW2F8FkVlvx40z1pIQ8LEhGBQny0/wcXajjSrA101Y0Ui2h+NcGj6kT9ZXb\nU5Dui2mSj+WdPcwRWnkbvVeU4T1Y9CI2YStEKI924qxo0HlximdNyU/zI0m0czrr+F0Kf/caPCa+\n3Geck1qSpFAYVkY+L334NsNS6rFNvDjGiHRk9NC1jhMa15k0nl3FbDazefNGcnPn9Oh5BIL+hBDh\nLqLI4LDqNIfu+0nVz413A85NDe1TGhFvtRU/Lhhmjmgm3ZrczTNWHWBZir067OgA1NHrPBnaJg/p\nKDpG1zSW6bvtg0ZnrQnh+R+Z72Zirptmr8zqqvSEDkxdEad440+mfrAiQ6Ytdix7W4aPOYaVL/6c\nkwvmYUvPi54zOhpy+yJbBnQmjefhUFVVicvlIi2tMyloBYKBixDhbiSRkxfEvgEXZ3gPekxHE2+1\nlWh13RXzYTLzCNNRNMIZqBLTfpWu68Zex0ZjifWwECaZB6J4KDJk2jWGJ2GeTiRORmblWOPvTInG\nZP0TrI4czj3vInbsXk3R1LOi3pckYghwe3+Gjg8xPVWFyWw2s3Hjeo455rju61Qg6McIEe5l4t2A\nk9kr7Egi82ZPrl6SSbl4iEM39dxUP6VOL3Zzmxt4HK/vWMR7WOjMg0QsOj4wheORO3pHG5FMtaKO\nY+xs/eBkCmDYTBpzZh/Psa4WPlq/hoBlGNa0HCTVTUG2mX0tsT3Cw2FPAJ5A6Lua7NbL4bBnTwVH\nHz1DxA0LhgRChPsIo1CproppMubZniCRKbwjHWNZ+zuxnM2SWQEmY1buSDyrhlHijLbjW7/fblh0\nJPzdSU93cOmZeWzctIPn//4olXu2c+a8s7BPvDjmg1+GNcjmWntcJ8SeCGGSJIlNmzYydeq0bulP\nIOjPDJDb4dBgYq6bEqcHm0klZK5VKXF6Eopp+GY8p6SRU0oamVPSyOS8nq8PGxaNZLCZtAElwG0J\nPzApcvt/xyKRWVmNccnCD2JG+FSZTyqcbDhgR9Ojj5syLLnvTmHBSK676jKynBm88/YbeGo2G54v\nHNNe1mg7KNKhmOFYmbbizauzSJLE7t270LSedQQTCPoDA/CWOHg5XDFNRiC6k3ii0ZHecuzpDyRj\nVo5FxwexQxxKnLGpJrqyVGe+O1lZWdxyy2JycnJ457m7oWFjlHiPz3Z3Yqsh8bw6SyDgY+fO7d3W\nn0DQXxkit8WBRW+L6eHQUTSsiorDEuj0an4wEc9CkMibOiymJxU1kqIYt0u0mk7mu5OVlc3Pf76Y\nrKws3n72l6TWfNBOvP1q57YaujuESVFM7Ny5o9v6Ewj6K2JPWHBYxHM06wnv2YFAdzjLBbTYdYq7\nK3FGdnYON998K/fe+ys+eP9tjpk1C+lghYvOpmjtCUtHU1MT1dXVDBs2rHs7Fgj6EUPs9ijoKTqu\nwAbSar4n6Or+fpjDWU13htzcPGbMmElVVSWbNm2MvK7Iodj12Oj0tKXDbDazdavxnrVAMFgQK2FB\nv2Sgr6QTxTInojdDz+bOPZ2vvvqSjz5axqRJkyOvlzq9Mat9Qe94u1dXV4nkHYJBjRBhQb8imfja\ngcThxCv3VuhZcXEJY8aMZdOmDVRWVjJixAgA7Ob4iWB6w9vdbDazadNGZs06pmdPJBD0EQNwjREb\nv9+P3x9A01R8Pj9+vw9d19F1Ha/Xi6ZppKU5GD58RLu/7OxcVFUVIRH9gHB8bduwmFgewYOd3gw9\nO/XUUC3fFSs+jLwWz/u9N73dKyrKUFW1d04mEPQyA34lbDabKS0dhdPpJDc3n8zMTBRFwefz4XK5\naGioQ1U1hg8fjsORHnE86UgwGGTHjm1UVFRQV1eLLEvYbHbsdjt2eyomkwmz2YLJZEJRZGpqanG7\nG9B13bBPXdfRNA1N0w7eQNqHm4T/K0lSKHUgEpqmYjabUJT4H0u4z8GWUagzaRuHEt2R/SsRU6dO\nIycnly+//JwFCxaSlhYqQN5XiWA6smXL5namcoFgsCDpup4gdXv3EgyqmEzJeVz2FX6/H7PZHFOw\n27bbuHEj5eXl6LpOamoqaWlppKWlYbFYsFqtkT9FUdr1J8syiqJE/jRNw+12U1dXh8vlwuPxoKoq\ngUAgskq3Wq3Y7XZSU1NJTU1l//791NbWUldXh8fjifQdXv2HHhzazyMYDBIMBrFarQQCASRJwmRK\n/CwWfqhQlJ777Fq8Em9tiLUHqTN/sheHtVe/rkOK9957jxdffJELL7yQ8847r917QRU8AQmbWacv\nfr4Wi4Uf/vCHCX+TAsFAo9dFuKYmQU28AURurqPfzKelpZlgMIgsy8iyjCRJuN0eXK5mPB5v5MEi\nMzOT7OwcrFYrqqpSWbmP/furqa+vx+VyEQj40TQdXVeRJJm0tDRGjSpEUWykpFiprq5k//79+Hxe\nzGYLuq7j94dMlmaziWAw9MAgSRJmsxlZTn7pqmqh+sPGe5Aqc0oae2wlnJ5uo7nZk7jhAKKjc1ui\nOXq9Hm6/fTEpKRZ++9uHkno46y2CwSAFBYVMnXoUqampMdv1p99kTzIU5jnY5pib6zB8vf/8ygSH\nhcMRXVs2Lc1BXl50+bowiqJQWFhEYWFR5DVd1wkEAvj9PhTFhM1ma/djKC0dha7rHDiwn+rqalJT\n7WRl5eBwODCbzWiahs/nw+Px0NhYT1VVFdXVVaiqmvCm3pfFKAYTsZzbjnXEf962Wm2ccMKJfPTR\nMlav/ppjjz2+l0acGJPJRHV1FXv27KGgoICJEyeRmZnV18MSCA4bIcKCdkiShMViwWKJnbJQkiTy\n84eRnx+dREGWZWw2GzabjaysLEaNGoOmaezZs4fKyj3oOsiyEtn71jQVvz9AMBggEAhwpNyMpgap\n8dgN9yA1TSMYDOJ0ZpKTk4OiKFRWVtLS0hw15rDp3WQy9atVXU8Tq3iEZU+Acc74q/1TTjmNlSuX\n8/bbS5k+fWa/8zswm03s31/N3r17GD9+AtOmHdXXQxIIDouhc2cS9BmyLFNcXExxcXFS7U9raWbn\nrnL2VNeSIuvYrGmYzVmkpFhIT0+nsLC4neAeddR06uvr2bVrB83NzQf3zNNwOp1kZGTQ3NyMy+XC\n5/Pi8Xhwuz243a243S4CgWDMBw5N0/D7/SiK0u/EKBbxnNv2NSqMTieuRSEnJ4eT55zOl6vXsPLj\nVZx+6mk9NNLDw2w2s23bVrKysigqSu57JRD0R4QIC/odDkc606YeybSpyR+TlZVFVtasmP0Zoes6\nHo+H6upKZDlIeXkoMYTNZicjIx2n00l+/jCam5upra2hoaGRlpYmNE3DZLJgsZgwmSwEg35crla8\nXg+gYzKZ0XUN7WC5o5DbRfjfoW2AnlqZxyse0eqXaPXLpFuNPa3DZuzMo65kzpEyja46vt9nYuqI\nYNJhUaoGrX45VB/a3LOJVkwmhW+++Yrs7Jy4+8QCQX9GiLBgyCJJEna7nVGjxpCb66CkxNgJJC3N\nwYgRIxP2FwwGaWxsxOVyYTIpBz3gTe2c5WRZoqWlhcrKSurqamlqasJsNkW8z3VdQ5LkLq+8E+V8\n/qYyPWbyk7ZmbFkGW3oe+1rBXONJWCtY02HjATt7m1NQ9VDHiqRTmOHr0UQrkiTxyScr+cEPzhKe\n04IBiRBhgaCbMJlM5OTkkJOTE7ed05kZcYbzeDxUVVWiKCYsFjMWSwo+n5eKinKqqioJBAKdWjXH\nc25rm/wEaCes8czYlU0KE3KIW5hjU42d8qb251R1yfBc3Y3L5eKbb74WWbUEAxIhwgJBH2Kz2Rg1\nanTU6yNGjETXdfbsqaCiopzq6ip0XUuYyAUOJdiobrHgVWWM4q47Jj+JZ8b2qQrugEx5k9UwaYeu\nh84Vi55OtCLLMmVlu8jIyCA3d2bPnEQg6CGECAsE/RRJkigqKqaoqBhN0ygvL2PPngoOHKhGlmNn\nzAinuyzK8PJxudOwTcdyiPHM2J6WGtaW6zSS2eb4QyvqUqf3oNgb012lF+NhMplYu/Z7GhqqKS2d\nQH5+fo+dSyDoTkTkpUAwAJBlmdLSUcyePYf58xeSm5tHMBiMe0yqOflyiPHyRO/f+S017hTj91wW\nzLKGVYktsN1ZejEeJpMJl8vFqlXLWbVqBS0tgyfRg2DwIkRYIBhgpKSkcNJJJ3Pssccjywqxkt51\ntgBDrBrIUssWFGuGYT+eoExAkxnmiF17uLcTrZhMJurqalm27H0hxIJ+jzBHCwQDlKKiYkaMGMkn\nn3xMQ0OdYZvOFGCIVQM5eMzRbGupxZ4RbeINr3LDe8OxvKP7ipUrP+KMM84kJcV4JS8Q9DVChAWC\nAYzJZOLYY49j6dJ/Y7FEhzW1FVaT1U7Qm9hBqmPVpimTJ/PNW2sNRbjtKvfIfDcTc929FiecDIFA\ngBUrPuT00+f1aPERgaCrCHO0QDDAsdlsCR2RFBkcVr1LoijLMqMc9ez69k0ItNDWVN1xlavIkG7V\nSE/pewEO43K5+PjjlTHN9gJBX9JPfiYCgeBwGDVqTEJHrcPh+OOOZ9dXL/HZ329hdmEdc0oamZzX\nc0k4uhNZlqmtPcCqVSsoLy8nEAj09ZAEggjCHC0QDAKKiopYs+Y7VLVnhNhms3H88SexfPkytqz/\nilmzju10H7ESffQGimKivr6Ompr9qKqGw5FOVlYWRx45lbS0tN4djEDQBrESFggGAZIkUVBQ0KPn\nOOWUU5EkiY8++rBTpl1Nhw0H7Kwsc7KizMnKMicbDtjR+sA6HMpMZsHn81JVVcm7777Fl19+jtfr\n7f3BCAQIERYIBg0TJkzC7+85U2tubi5HHjmV8vLd7N69M+njwjmpQ4lADqXO3FRj77GxJovJZKKy\nch9vv/0G3377DXV1dWLvWNCrCHO0QDBIsNvt5OXl0djY0GPnOPXU01m3bg3Ll3/EqFFjEraPl5O6\np9NZdgZZVigvL2Pbtq2YzSYyMpxkZmZx1FHThVe1oEfpB19/gUDQXZSWju5RB61x48aTnz+MtWu/\nP1i6MT7xclKH01n2J1JSUpBlhZaWFsrKdvPFF5/19ZAEg5z+9QsQCASHRUlJCRZLzyWmkCSJmTNn\nEQgEWLt2TcL24ZzURvRWOsuuIssyVVWVbN68qa+HIhjECBEWCAYRveGgNWNGqGTg6tVfJ2zb2dSZ\n/Q1FUdiwYR3V1dV9PRTBIKWf/wQEAkFnKSgo7NFY2GHDhlFYWMTGjRtpbXUlbB8rJ3VfprPsDIqi\n8Pnnn+J2D4zxCgYWQoQFgkFGTk4uWg/H/8yYMQtNU/nuu28Ttg2nzpxT0sgpJY0DKtHHIXRWrVoh\nEn0Iuh0hwgLBIMNkMuFw9GwCihkzZgLJmaTDKHJoj9gblFH771ZwTNzuVpYufZ2vv/4SjyexU5pA\nkAwiREkgGIQ4HBnU1dX0WP/Z2TmMGjWGbdu2Ut/YRIo9M24mLE0PxQuHqzlZFY1se4DJea2YB0gE\nkCRJSJLC3r17KCsro6BgJEcfPROr1drXQxMMYIQICwSDkIyMnhVhgBkzj8FaOJvP9+WBObVdmcSO\npuZwwo4wXlVhX4tCtcsSKXc4kMzTJpNCdXU1//nPe5xyymk4HI6+HpJggCLM0QLBICQ7O7vH8kiH\ncY49nVHTzwVzGvEyYcVL2KHqcr/JntUVgsEgy5Z9QH29cT1ngSARQoQFgkFIfv4wgkG1x/pXNWjw\nG6/+9rss7fZ84yXsiHXMwEJnxYqPRBiToEskLcIul4vm5uaeHItAIOgmLBYLdntqj/XvDiSfCSte\nwo5Yxww0JEni009XsXHjBjRtwD5NCPqAhN/6iooKzj//fObOncupp57Keeedx+7du3tjbAKB4DDI\nyMjosb53N1oB403ccCYsVYNWf+gWEythR8djBjKyLLN580aWLv0Xa9eu6dH0oYLBQ0LHrHvuuYf/\n/u//Zt68eQC8++673H333bz44os9PjiBQNB1MjIyemSvUtWgptV4jxcg1+5nc+0hT2ibSSMv1U+J\n08OephRUPfrZfyBkz0oGWZbRdZ3t27eyY8c2ioqKGTduPBkZzr4emqCfkvBr39DQEBFggLPOOovG\nxsYeHZRAIDh8nM7MHjGNxt/j1VF1okoXljeFPKNPG9VAgWPgZs9KFlmWkSSJPXsqeP/9d3j//XdY\nt24tfn98i4Bg6JFwJWyxWNi4cSOTJk0CYMOGDdhstgRHxSYz047JNEACA5MgN3dohCYMhXkOtjk6\nHOPZuPG7dnGs6eld/+2GsauQWqnT6o82R5vw0eA1LiBR405hVqnO7EydoOrDE5CwmXVCt4PDH1eY\n7phj92IDdPbvr6CqqozRo0czffp0LJbY1oRkGGzfVyOGwhwlPUEF6zVr1nDLLbfgdDrRdZ2mpib+\n8Ic/MG3atC6dsKampUvH9Udycx2Daj6xGArzHKxzfOONf6FpIS/p9HQbzc3dk+lpw4H2cb9hmvd8\nS3rh0RjvF+ucUtJIqiW0Ole10Ko6XpKPztKdc+wpQrdcndLS0UyZMg2TqfPpGgbr97Utg22OsR4o\nEn7606ZN44MPPqCsrAxN0ygtLT3sJziBQNA7pKen09jY0O39hs3Hbfd9K7d9zpfv/B/n3vAcAaKz\nSIWdrzpmz4qX5GMwIkkSILF79y72769m3ryzD74mGIrEFOFHH32UG264gTvuuMPw/QceeKDHBiUQ\nCLqHjIyMHhHhcFGGCTnuyGp2ux8+C3ip37MGR+GxUceEna86rqLDST4g1OdQQZIkWlpa2LJlMxMm\nTOzr4Qj6iJgiHN4DnjVrVtR74qlNIBgYOJ1Odu/WkOWecT1WZCLm5fHjJ1BSUsrHrz3EZYufxEVW\n1Eo3Xvas/S4LE3Lcg8JLOllMJhObNm2gtHSUyEE9RIn5dZ87dy4ABw4cYOHChe3+du3a1WsDFAgE\nXWfYsBG9Vn5PkiTOPPMcdF1j06pnDEsXxvOsHugJO7qKJEl8881XfT0MQR8RcyX8yCOPUFdXx/Ll\nyykrK4u8rqoqa9eu5ZZbbumN8QkEgsMgNTWVlBRjb+We4MgjpzByZAFff/0VZ5xxJiNHFrR7P5w9\nKxS+1J7BkLCjq1RV7WPfvr1R10sw+In52HnGGWcwa9Ys7HY7s2bNivydeOKJPPnkk705RoFAcBj0\nZJjiZJYAACAASURBVOasjsiyzMKF56PrOq+++godgy8UOXb2rMGSsKMrmExmvvvu20hcd1NTI19/\n/RVvvfUm69atFakwBzExV8JTpkxhypQpnHZa+zJduq6zd+/eXhmcQCA4fByOdJqamnrtfJMnH8mk\nSZPZuHED69evY8qUqe3eN/KsDu8ZD2V8Pi+rVq3A7/fT2FiP2RzaO9+2bQu7d+9k4sRJjBkzTvjk\nDDIShii9++67/O53v8PjORR7N3LkSD788MMeHZhAIOgeMjKc7NlT0avn/K//uojNmzfx+uuvMHHi\npHaxsEae1YlWwD0RU9zfkGWZhoZ6gIgAh19XVZU1a75n+/btnHzyKUMiicVQIeHX+cknn2Tp0qWc\nddZZLFu2jF/96ldMnTo10WECgaCfMHLkyF5PlzhixAhmz57D/v37WbVqhWGbsGd1PFHV9FBI08oy\nJyvKnKwsc7LhgB0tboqhwYmiKHi9HpYt+0CkDh5EJBTh7OxsCgsLGT9+PNu2beOnP/0pW7du7Y2x\nCQSCbiAtzUFKSu+Hv5xzzgLsdjvvvPMmLlfXMh9tqrFH5aEua7SxqcbevYMdQGiaynvvvUdtbW1f\nD0XQDSQUYZvNxpdffsn48eNZsWIFNTU1eL3e3hibQCDoJpzO3nPOCpOWlsbZZ5+L2+3mrbfe7PTx\niWKK1SHuq7Rq1XKqq6v6ehiCwyShCN91112sWLGCk046icbGRs4880x+9rOf9cbYBAJBN+F0ZvbJ\neU8++RTy8/P5+OMV7NixvVPHipji+EiSxKeffszKlcvZsGE9ra2tfT0kQRdI+C1+++23ueOOO5Bl\nmUcffZTVq1dz2WWX9cLQBAJBd5GZmdUnYS4mk4mf/ewyAJ555klcLlfSx4Zjio0IxxSrGrT65SG7\nKpZlmfr6OrZu3cybb/6bt95aKlbHA4yEIrxixYqoWD+BQDCwGDGi952zIGRSHlF8BOec+0MaGhp4\n/vm/Jn0/iRdTnJfqZ3OtcNgKI0kSKSkpBAJ+PvvsE6qrq/t6SIIkSRii5HQ6mTdvHpMmTWqXeUcU\ncBAIBg4pKSmkpaVRX987peGiKiWN/wknLCzl8zf+wPLlH3Lqqacn1U+smGLAsAiEZU+Acc7+Xcqw\np5Ekic8++5gTTpjNsGHD+no4ggQkFOGFCxf2xjgEAkEPk5WV1WsiHPZqDuMJKmSWnsiU01z8619P\nM2bMWIqLSxL2YxRTDLCyzGnYfl+jwuh0Bm0scbIIIR44CBEWCIYImZmZQHmPnyeeV3Pp5JNZv+Kv\nPP30E/ziF79sl40vHm2rNbX6YztstfolvEE50nYoExbio4+eTmnp6L4ejiAGQ/x5USAYOgwbNqxX\nKirF82oOSlbmnf1f1NbW8Kc/PdKl+OF4DlupFn3IFoEwIlSh6WuWL/9QeE/3U4QICwRDhPz8/F45\nTyKv5rPOOJWTT57Lvn17+eMfH+mUxzTEd9ga6VSHvCm6IyaTicbGBt5//202bdooHG37GQm/rpWV\nle3+qqqqqK+v742xCQSCbkRRFNLTez5pR6JKSSZF4kc/+gmzZ89h3769/PnPj9Da2jkhnpjrpsTp\nwWZSAR2bSaXE6eHowt6pnTwQkSSZjRvXs3z5h71WY1qQmIR7wtdffz3bt29n3Lhx6LrO9u3byc3N\nRVEU7rvvPo477rjeGKdAIOgGnE4nHk/PVytKVClJkiR+9KOfous6n3yyij8/+r9cfd2tZKalRK1k\njYo3xCoCIUs2BLFRFOXgqvhdZs8+mYwMYwc3Qe+RUITz8/O57777mDx5MgBbt27lscce484772TR\nokW8/vrrPT5IgUDQPTidTior9/V4ObxkKiXJssxFP/oZKQWz0ewFfFGZg82sM8wRiIh1uzCnNkIu\nHxx+W4etw2EoVGkKI0kSgYCfDz/8D8cccxwFBYV9PaQhTUIR3rdvX0SAAcaPH09FRQXDhw8XhaYF\nggHGiBEjWbPm+3Yx/z1JIpHcUpeGdfj0yP97VShrPHRbMooFhpDAdwdR8cwGQj9YkSSJzz//lPz8\n4ZjNJmRZwWRSyMvLp6iouK+HN2RIKMKFhYU88sgjLFiwAE3TePvttykuLub7779Hlgf5I6NAMMhI\nT8/AYjEOH+pt4oUyVbdYiLVY3++yMCHHHVmxtl3FdhajeObuFvr+jMlkoq6upt1ru3btwuv1MG7c\nEX00qqFFQhV96KGHUFWVW2+9lTvuuANN07j//vvZs2cPS5Ys6Y0xCgSCbkKSJJzO/rEPGL9AgxS3\neEOrX6bFJ7N+f/vUlasrzEmnrhRVmowxmRTWrPmeTZs29vVQhgQJV8JpaWlcddVVzJw5E03TmDZt\nGmlpaZx77rm9MT6BQNDNZGQ4aWpq6uthREKZQrWC2+Np/v/t3Xd8VGXaN/DfmTN9MkkmyYSQRkIg\ndAKEktCTIF1UiqyC6PP6PO6uqKzu+q7u6664rhXXgq7bLKuuu49KVVcB6SAlEEPvJZBQkhBSZyZT\nzpz3jzgh5UzNlDOZ6/v5+PlIJjNznwnkOvd9X/d1VUGl0oCRdy7mwTI8DlyN/jFI35oum2wszlSx\nsFjUHs1iPenSFMqiH6Hcp5ZKpTh27AisVitycoYF980jjNsf7a5du3DHHXdg7dq1WLt2LebMmYNt\n27YFY2yEkACIiwtNR6WOXB1lqr96DGXHdwo+xvGSHwO38Hq1p7NYT7o0AQh6pyY7Dxyr8q05hT/H\nKpVKcfr0SezevYuOpQaQ25nwG2+8gX/9619IS2vJoCsvL8cjjzyCgoKCgA+OEOJ/qanpOHCgWBQ5\nHR2PMilYOxI1FqSMHIWyaq7lZoHnIWElULJ2WO0MON71uD2dxTpuAtruCTv0iLKAYVqCYbCTtnzZ\npw5UgplUKkVV1XVUVFyGVhuNlJQUZGf3h1qt9v1FSTtug7DNZmsNwEBLopYY7qIJIb6Ry+VISEgQ\nxZK0hGkJxDwPXG+Sw8xJcLVRAY5XQdbh97ytsRyc2n3WbttZrDuuzjOHImnL3T5124S0tgI9Vrlc\nDrO5GRcunMfp06eRnp6O4cNzoVQqvX4ts5VDfZMZMVEKKGSdtyIijdsgnJycjH/84x+YP38+AGDV\nqlVISUkJ+MAIIYHTs2cKamtrRTEbPlGtxqX6WwGE44Wnbg1WJSSNNVBqE1y+Xo8oi8d7qM7OM/sa\nDLvKl33qYI9VJpPi2rWruHKlHH36ZGPIkBywrPtgytnt+GzrOZSeqcbNBjPiohUYnq3HwsI+YEXw\n9zBU3F75Cy+8gEOHDmHKlCkoKipCaWkpfv/73wdjbISQAOnTp68oVrRcBZCOVNF61FQcdfJoS+nK\n7MRbhT684TjP7AhWngTDQPB0n7qtUI1VImFx7txZfP31etTX17n9/s+2nsPmgxWoaTCDB1DTYMbm\ngxX4bOu5gIwvXLidCcfHx+PNN98MxlgIIUHiWJKuq3P/yzOQXAWQjlRSYOSIXFxrMoOXtARuu82M\nlCgT+iUBapkdulgVGhq6Pi5XmdveLHd7y90+tdCMNlRjBVqqnnEchy1bvsOkSYWIj48X/D6zlUPp\nmWrBx0rP3MC8SVkRuzTtNAgXFha6LG23ZcsWn95Qp1NDKu0+H7Ze71k/1HAXCdcZadc4aFA2Dh8+\nHNIlaTUHaK7yMFjcZw8p5cA1UzTA3sqLZmVKHDiwDaoB8tY69tHR/qkfnRZnx5mqzr+r0uJagn2g\n5Gl5yMutuFLHwmBhoJHzSInlMCKNb1cbu+11hmqsbf3wwx4UFhYiKSmp02PXbhhws9Es+Lzaxmaw\nchn0CZpOj0XCv0mnQfiTTz4JyBvW1nafKjR6vRbV1d73Qw03kXCdkXiNOl1P1NZ+D5lMFsJRAXo1\nA4Olc6BgGTs4noFKaodeY0G1QXjZOr7XCPzpz49i+/aduOeehYiL80/Lxj4xJlgsnTOO+8QY/TLb\ndiU71oSs6PbnhNu2Xo6OVqGhwSSKsba1du3XyM8fhx49kmCz2WCz2WC3c1CoohCnVaCmoXMg1mmV\n4CzWTv/+utu/SWc3FE6DMCVfEdK9yeVy6PWJqKurDek4nGUo94s3wsK1BKFmmwSX64UzcdUxieg/\naASOHt6Ho0cPI3fkGEyZfhfSk/VdSkjypAlFoHhbqMNVgpnBGryxSyQS7N69EzzPw27nAbT8l5mZ\nheF99dhcUtHpOcOzEyJ2KRrwYE+YENJ9JScn4+bNmpAuSbsKdjK2ZT/T3b7nww/9N06dykfJJUCl\n74+jTQk4cqwJKXES5PS0oCu/4/3VqckTXT3v6xiro+BHKBpTCK2sVFSUo0d8M6bkpqL0bA1qG5uh\n0yoxPDsBCwv7BHZAIkdBmJAIlpXVF0eOhHZf2MFVsHOXsCRlGUA/GvGyNo8rY3DdCFSes6KXzhYW\nnZH8dd5XbI0pWJZF7c0biIsyY/kDhTCY7XRO+EdOg/C6detcPvHOO+/0+2AIIcElk8mg1yeitlb8\nZQldFdZwddSJZ2Qoq5PBZrViWIo1mEP2ir/O+4bqjLM7EokEBkMTtm7ZhD59siDt0RNynS7gva3F\nzmkQ3r9/PwDg8uXLuHTpEiZNmgSWZbF792706dOHgjAh3URycgpqam6IYjbsiqtla4PV/VGns9ea\nYa8+hhHDcoIwWu/5q6GEmBtTMAwDq9WMkydP4NChUkilLKKiopGQkIDMzCwkJAgXYunOVbacBuGX\nXnoJAHDffffhyy+/RFxcHACgvr4eS5cuDc7oCCEB17t3Fg4dCp/+4ELL1q72jB0Umjj88x+f44eD\ne7Fw4SJoteI6/uKv876hPDfsDYVCAQAwmYwoL7+MCxfOQaPRIjk5GQMHDgagjYgqW26voqqqql3/\nUZVKhepq4UPXhJDwI5PJoNPpQj2MLnHVkclBKeXQMzEWBw8ewHPPPYP9+/f6rWqYL92LOj7H1TV4\nU4rTX68TbDKZHBaLGWVlF7Fp0wY0NzdHRJUtdvny5ctdfUNFRQXef/99GI1GHDlyBK+99homTZqE\nvLw8n97QaHT9DyWcaDSKbnU9zkTCdUb6NdbX14W8elZXJaitYKQy1BkZ8AJtDtNirbijaCTUag2O\nHz+GkpIDOHr0CPR6PRIS9D69p50HjlercbxKg7M3VbjSoIDRKmkZy49D4OyAySoBy/CQMK6fo9dY\nYbUzMNsksNlbzkinRpsxUG9E261ThUIGs9nm8rPw5HXEym6343J5OXadssJk5jo9Xt9kwaRhyZCK\n9Y5CgEajEPw6w/O82y6VGzduRHFxMRiGQX5+PoqKinweSHc7fN2drseZSLjOSL/GqqoqbNv2HWQy\nz+o4i1V0tAo1tSYcr1KjxuT8eM6NG9VYv34tDhxoyX0ZOHAQ7rxzHtLTnXdpEjq7e6xKLZixnRFr\nau3E1DGRDIDT5zgyl92dE+5YrMObMYcLRq7GV8eUEOobLWGAFx/KQ6IufFoqel2so62EhAT06dMH\n8+bNw+HDh/06MEJI6On1+rAPwA4yFhjW0wjO7rzIRkKCHg8++BBuu20a1q5dhRMnjuPEieMYPjwX\ns2fPQUpKauv3Oju72y/e6DILmefRrjuU45gQywivWbfNXPbX2eRgnnH2N5WMd7q3rdMqERMlPLMM\nN26D8EcffYTNmzejqqoKM2bMwO9+9zvMnz8fDz74YDDGRwgJAoZhEBeXgJs3b4R6KH7jSQBKT++F\nZct+iRMnjuOrr9ahtLQEpaUlyM0dhVmz5iA5OdnpmVsrx7jMQr7uJEA7a9UY6sxlsZGyzs+Gd6cq\nW24XKNauXYv3338fKpUKsbGxWLVqFVavXh2MsRFCgkiv18OD3aluaeDAQfi///c3WLp0GdLTe6Gk\n5AD+8IdnsfG775zOdmuMMqdtBxWsHWbOu/VfMWUui0W/eCNStSaopBwc7Spz0uWYOyEj1EPzG7d/\nSyQSCeTyW38JFQqFRw2cCSHhJSMjExZL905Oc4VhGAwZMhRPP/1b/OxnjyAqSouNm7fD6KTDUzMn\nQbzKeRayswDNMsI3OmLOXA42Ow8cvCzDzkuxqGhUgueBFK0ZE3vVIVVxFV9/tQ5HjhwWRU/srnK7\nHD169Gi88sorMJlM2Lx5Mz777DOMGTMmGGMjhASRRqOBVquN6EAMtATjYcOGIzOzNz748AOYGquh\njuncmUkltaN/QksilVAS2Ilq4QSs1GgzGEa48hdp0bIFcKsGdTPH4kojCxnLY3BiS4b3mTOncPHi\neQwaNBh9+mSHcLRd4zY72m634/PPP8eePXtgt9uRl5eHe+65x+fZcHfKQI2EjFogMq6TrrHFvn17\ncPXqlSCNyP88zRr2FMdx+Lq4Emz8kE6PaeVW2Owt1amUrB3xaisGJxpam0XYeeB4lRrXm+Qwc52z\ntLuSuezv6xQTzg5sL4t1UmyEw+SMunafF8dxiI6OweTJha0FQMTIWXa02x/9Rx99hLlz52LlypV4\n5513sHjxYrz++ut+HyAhJPR69EgCxzk/fxoMvhS+CBSWZXF7XjKirJdhaqiCnbOhufEGFGhCo0X2\nY6BgfpypKXG6puXIjCOjusrQEoAVbEtP5LbHpByJY7QE3Z4nZTfbYlkWBkMTNm3aAIPBEIwh+pXb\nH/+bb76JRYsWobKysvVre/bsCeigCCGhkZaWjlBtszna720vi8W2slhsL4vFsSo17F3MFetqUJcw\nwORBGkzrZ4bh+MfY/skvUFsv/Mu+skkOzn6ri5EjSJs5FpfrVThRHT7nWkPFUXZTiKvkNavVgu++\n2xAWzUjachuEMzMz8dOf/hSLFi3CwYMHgzEmQkiISKVSxMWFpoRlx8DlOArkTeCycbcCrrdB3V2w\njtaqcc/8O/Dw0l9CpRVuNGCySWCwSFyeHxbDDF/MulJ20263Y9u2zbh+/XqARud/bhOzGIbBlClT\nkJqaimXLlmHJkiWCTZsJId1DQoIejY3B3R/vavs9x/JvtVEBg0UFldQOqcSORsut31XOeuo6K8bh\nrP9wv6x0XC3j0Sywaq+S2gEGou1iFC4G6o2Qy6UovynxIXmNwc6d25GWloaRI0eLPl65nQk78rb6\n9++Pf//739iwYQNOnjwZ8IERQkIjNTU96BnS3u4DduSYRRssEjhm0W0DcFsdZ6PezsBZCZAUJdyX\nuOz4Thwv3fvjudbO6CywZyQMMDLdiskZdSjIqMPkjDoMThS+KRIilbK4evUKvvpqHc6cORXYwXaR\n25nwn/70p9b/j4uLw4cffogNGzYEdFCEkNBJSEiAWq2BzSYcaAKhK+33XM2ihbSdjfo6A3fMyByz\nZzljRW3FIfyw4U84wNkwfBqLlEHTOj3P2XKqu0xpMdaADsaYulJ2k/mxU8WhQ6UoK7uI/Pzxomtf\nCbgIwm+//TYeffRRvPPOO4KPz549O2CDIoSEDsMwyMsbix07tgatMI9jH1DoXK27fUBXs2ghbYO6\nJzPwjkHAEXwGJBgxIKFNfeq+WSga8gK2bv0O32//EOZmM5Kz86HUxkMl5QWXU90thTt7PE8buspm\n3i7fh5pUKkVTUxM2bfoWI0aMRGZm71APqR2nQXjQoEEAWop1EEIiS2JiIrKz++PMmZNgWY/6vHRZ\nx9mlp/uArmbRQtoGdW9m4J4En4SEBNx99z2YNet2bNjwDbb+4xEo1DqMyc3B+DvugIRpf47VWV1q\noGXf2tnj8nIrsmNvnRMO5kzZ3ZjFimEYHDiwH9euXUVe3lhIJOJYUnBarOPq1asun5icnOzTG3an\nggiRUOABiIzrpGvsjOd5bN68MSRJWt4GFGctBdsW1HA2Y3PVjrBtUPH0+9q6cOE8Pv74Q1y/fg16\nfSIWL74f/fr1b71OV0UpJqTXYddl4cc1cjsmpteCYYI7K/W2kEZXBKogid1uh0qlwtix46HTxfn9\n9Z1xVqzDaRAuLCwEwzAwm82oqalBWloaJBIJLl++jPT0dJ/3hbvTL7pI+MUNRMZ10jUKMxgM2Ljx\nPxDq6Som7bOjmXbBiOddB3VPZrhdCT5WqxVffrkOmzdvBM/zyMsbi3nz7oZEEYNtZbEQ/mx55KU0\nYN+VaKePF2TU4WKd0usbg64wWCQux1yQUee3zO9ABWHHTZ6MsWH4sJzWm6JA87qf8NatWwEAjz/+\nOBYtWoSRI0cCAI4cOYL33nsvAEMkhIiNRqPB8OEjceDAfkilwVmW9oWEaVkKVWt4VNea2wdcxnVy\nj+O57fZ3OwRUX/aOHWQyGebNW4Dc3JH49NOPsW/fHhw5cgh3zf0JVMmznC6FRytsTpfKNXIeMom9\nS8e6fNGVBLpQE7rZOl1zFuOvXMHECRNDdpTJ7b+q8+fPtwZgABg6dCguXrwY0EERQsQjM7M3zp49\ng6Ym8a8USFnfs2nbZuJ2XBL3R/DJyMjEU089gx07tuHLL9fi039+gPF3xSM2c1yn79VrLJBLnSer\npcRysNp9vzHwVVcS6EJNaC/7UoMGzBkjGuq/RExMNHieB8/z4Dg7pk6d3pphHUhug3BSUhLeeust\nzJw5EzzPY/369cjIyAj4wAgh4hEXFxcWQdgVT/aaXS1N+yP4sCyLwsIpGD48F59++hG+X/dH5M5o\nRvqACbDwt2Zi1QY5jlWhtUtTx/GMSONRX9+1GwNfk7l8TaALJddH0RQYkGBCXV1d69es1uCdk3cb\nhFesWIGVK1fiiSeeAACMGzcOL730UsAHRggRj549k3H+/DnRVx8S4iqwdtwvdpX568/go9Pp8PDD\nj2H16s+x5Zu3wdnsSBl8m+D7Ci2VSxiVz7PSrh4x8mT5Xmy6sp0QaG6D8Msvv0xBl5AI17NnsmiO\ndHjLWWCtMUrbZU7rNRZUG1zvsfoz+EgkEixY8BMkJPZEtXqoy/d1VrTClxsDfx0x6kohDX/wZiYv\n5r1st0H4zJkzMBgM0Gg0wRgPIUSEJBIJdDodGhoaQj0Ur7hahuxYV/pyvQqAcBGMtrMlfwef0fkF\n2HYxVvAxowUov3YDGSnCDSO8nZV2tUa3GPgykxfzXrbbICyRSFBQUIDMzMx2DZM//vjjgA6MECIu\n8fEJYReEva2m5UwgZ0tKqR0qmfAszdRYjdf/8hymTpmCadNmOt0O8PTGQMzLsp7ydSYv1r1st0H4\nySef9Osb6nRqSKXBKYUXDM7OfnU3kXCddI2uDRnSD9evXxb9vnB09K1f0GoO0FzlYbB0Lcs1Lc4O\nXWznWZS/pMXZcaaq8+/F5BgbotRKfP31lygtLcH//M//IDs7G0D76/SUq89DI+eh1yng71/PNg4w\nWRmoZLzXr93xGm0cUG1UCH5vtVEBtcb1e4yN4WHjzB3G0/lztFhY6PVacWRHjx49GidOnIDRaPwx\ndZtDRUWFz+Usa2vFm0HnrUgo8ABExnXSNbrHshoYDBZIJAI9/ERCqMCDXs3AYPEsYClZOxI1Flxp\nVIDjW34BswwPi8WGunpTwGoj94kxwWIRWmJVIfe3v8e6dauxc+d2PPfcc7jjjrmYP/8uNDWZfXov\nZ5+HXm2G0eC/4hhdTQAT+lkaLBKnP0uDhUF1rdnjmbzRRQK01WpBdXWjX4Ow18U6HJ555hkUFxej\nvr4evXv3xqlTpzBixAjMnz/fb4MjhIgfwzCIi4trd5QjHAgtQ3bsNeyQpG35zczxt5ZsOZ4JeG1k\nV3u7KpUK99yzGKNGjcF77/0V69atRnl5GRYtegBqtdrro0bBWpYNRI1pMSdY+cptEN6zZw82btyI\n559/HkuWLIHJZMLLL78cjLERQkQmLi4h7IKwUIBzVnO5X7wROy8JJ0l5krjU1UYKrvZ2+/Tpi//3\n/57F++//FSUlJbhcfgWzHvg9DIj3aqYZjCNGgUoAE3OCla/cBuHExETIZDJkZWXh9OnTmDVrVtAL\nuhNCxCElJRWnTp2AXO55/16x6BjghAKRweJb4lKw2vtptVo89tgT2Ljxa5yti8cNW2Kb8Xk30wzk\nEaNAJoCJNcHKV26DcI8ePfDXv/4V+fn5WLFiBQDAYgleNRFCiHgkJCSEZQB2pmMg8nW5M5jt/SQS\nCebNX4jVJQw4gce7MtP0V0vEQC4bh2OxEFfcDv2FF15Aamoqhg4diqlTp+Lrr7/G8uXLgzA0QojY\ntOwLx4d6GAHjWO4U4my5093SKxeAyabJyoCTKJ0+1uzlsSw739KqcXtZLLaVxWJ7WSyOValhFz42\n7ZYvn6Mv7+E4t83ZW1YxAvFZB5rTmXDbfsLDhw/H1atXUVRUhKKioqAMjBAiTjpdHG7erAn1MALG\n2+XOUJy9Vcl4pzNNY0MVPv3oQyyYP9/jfrnOZvJ2Hhjaw7eZfDCWjYO1DRBIToPw4sWLBfsJl5eX\nIzU1FRs3bgzmOAkhIpGeno6TJ49BLhc+rxnuvF3uDEXGrpR1nqBkunEGP5Tsx/Fjh3D77XegoKAI\nLOt859HVTP5yvRIMgEGJroOa0DJ2MJaNPd0G8NcyeyBQP2FCiFd0ujgoFErwvI9rlWHC08SlUGXs\nOptpTp86AH2iH8CaNV9g1arPsXv3LixYsBCDBg0RfB3XVcUYXKpXgWGE97Y9mYkGKgHMkwxsZ1nw\nYpopUz9hQojXRowYidLSEnAcF5SqQmIXioxd5zNNCcaNm4CcnGFYv34tdu/eibfffhODBg3GvHkL\nkZyc3O51XM3kHZwlezmbiVo5BkN7GAI66/RkG+BinTJoCXO+on7ChBCv9eqVgZSUVPzwQwnKyi52\nq1K0vghlxq6zmWZUlBaLFi3B5MmF+OKL/8Xx48dw8uQJTJs2A7Nnz2ldonY1k3cQ2tt2NRO90qjA\nTZMsoLNOd9sAMok9LJpVuB3CihUr0NDQgCeeeAK//OUvYbPZqLUhIQRSqRSjR4/BlCm3QaVSwW4P\nw9RUP2ubsRsonB1obGY8zgROSUnFsmW/xMMPPwqdLg7ffvsfvPrqS7h+/XprVnG/eCPSY0xwErWK\n+wAAIABJREFU1kVKaG/b3TK2Y9Z5olrt+cV5wV0GttXufqYsBm5nwjExMfjtb38bjLEQQsJQXFw8\npk6dgR07tuHmzZqw7Tssdp33X+UezzQZhsHQocPQt28/fPbZv7B//z6s2lGGzCG9wEs1rcvn6dHN\nuNzg2d62J8vYQGBnna62AXgeYVHi0m0QXrNmDV555ZXWFmY8z4NhGJw8eTLggyOEhAeWZVFQUITv\nv9+Fa9eugmUje3k6EPxREESlUuGBBx5Er5F3wyDPaJ33Ol6rV4wJGbEmj/a2PVnGbnntwLVIdLkN\nwIRHiUu3Qfjdd9/FJ5980to+ixBChDAMg/HjJ+LAgf24ePECpFK3v15CSszHVjryZy1mzg7Y1WmA\nQDOsKoMMkzPqPd7bdgTn641yNHMSAJ2n5MGYdTrbFw+HEpce1Y6mAEwI8dSoUWMQGxuLGzduwGQy\nwWQywWg0QCKRiCKT2h8FHoIdwP1ZEMTVaxksDE6cuYgh/TM9er22M9GjlWpUNIpr1hkOJS7dBuFB\ngwbhsccew7hx46BQ3Dqcf+eddwZ0YISQ8NW3bz/07duv9c9WqxU7dmxDbe3NkC9Vd2VZN1QVmvxV\nEISzt/zn7LWaG6vx549exYB+2bjrrvlITU3z6HVZCTA0yQgpK85ZZyCbVXSV2yDc1NQEjUaDQ4cO\ntfs6BWFCiKdkMhmKim7D/v17cfnypZAtVXd1WbcrAbwrs+euFgTpePPAMsJZ0Gk6Bn2zeuP48WM4\nceI4Ro/Ow5w5dyI+PsHtGMNh1ilGbv8lCB1Ham5uDshgCCHdF8MwyMsbi5iYWBw7diQkM+KuLOv6\nGsD9NXvuyv5mx5sHjm95Y5axg+OZNq8lRd4vfoUTJ45j7dpV2L9/L0pKDmDSpALMmDELUVFat+8V\njFlnOO3nu+M2CG/duhVvvvkmjEYjeJ6H3W5Hc3Mz9u7dG4zxEUK6mQEDBiImJgYXLpxHbe1NNDU1\nQS6XB2W/uCvLur4GcH+1OWw705Qq1bA1e5aM5ermQSbhMS61HhpZ22DGYNCgwRgwYCCKi/fhyy/X\nYcuW7/D997sxbdoMFBZOabc1GUze3tCEQ7D2aCb8/PPP48MPP8TPfvYzbN68GSaTKRhjI4R0U8nJ\nKUhOTgEAGI1GXL58CadPn4LNZg1oMO7Ksq4vAdyfWc0OrATQKnk0eNjW3WB1fvPQzEnAMhAcg0Qi\nQV7eWOTmjsLOndvxzTdfY/36Nfj++5148MGfIjOzt9v39ncQ9PSGJpy6K7n9WLRaLfLy8pCTk4PG\nxkY8+eST2LdvXzDGRgiJAGq1Gv37D8CsWbejR48k2GwCZ2f8aKDeiIxYE1RSDgAPlZRDRqzJ7bKu\nLz1yPZk9B4qjR3BxRbTT7/Ekqcuxn/+HP7yE226bjpqaGqxY8TI2bdrgtEqav/sTA971bXYE65Yb\npsBX7+oKt38DlEolLl68iKysLBQXF8NiscBqtQZjbISQCCKVSjF+/ESMGjUmoB2aHMu6kzPqUJBR\nh8kZdRjsplWfgycBvG2DecfsWUigz886AlEz1xKIhHhzfEilUmPevAVYtuwJREVFYc2aL/CnP73V\nWshJ6L39GQQ9vaHxJliLgduP/xe/+AXefPNNFBQUYO/evRg3bhymTJkSjLERQiJQZmZvzJgxG3Fx\n8QGdFftS59lVABea/Z28ofZ69uwPrgKRN7N/If37D8QzzyzHwIGDcfz4MfzhD8/i0KFSj967K0HQ\n0xuaUK4++MLtaHQ6Hd566y3I5XKsXr0amzdvxrRp04IxNkJIhFKr1Zg8uRBjxowFy7Kiaw4hFMCd\nzf54Hj4tf3eF6+YKwKjkBo9n/0Kio6PxyCPLMHfuAhiNRvzlL+/ggw/+DoOhyecg2HYFQYin2wGh\nXH3whdPErJKSEtjtdjzzzDN44YUXWpeHbDYbli9fjo0bNwZtkISQyNSrVy+kpqaipOQgysouQCKR\niLJBhKvZX5VBjskZdUE9P+suicwfR4gkEgmmTp2OwYOH4uOPP0Bx8T6cOnUSdy9cDJWu0OMENmdJ\nVHnazlsSnhzT6uqZ6mBzGoT37NmD4uJiVFVV4a233rr1BKkUCxcuDMrgCCGEZVmMHj0GI0bk4uLF\n87hy5Sqqqq6DYZiQV99y8PT4UrCqNgUzECUnJ+PJJ5/Gd99txNdfr8d7f/8Thk+zI2VQ5xVTofd2\nlvEsL7ciO7b9SRxPC4KEQ81oB6dB+NFHHwUArFu3jqpjEUJCTiqVtpbDtNlsuHr1ChoaGmAyGWE0\nGlFTcyNkY/NXWUl/CmYgYlkW06fPxLBhw7Fly3c4uPMjmJvNSOozGiqtHkoph57RXKf3drWCcKWO\nRVa08PEpdwVBwql6l8tzwtu2bUNubi4AYPPmzVi1ahUGDhyIn//855DJZEEZICGEdCSVSpGe3qvd\n16qrq1Fa6rqIUKCKN4htCdRxnQMS/B+IXH2GSUk9sWjREsyfvxA//HAQ3+96C1crayHhm7H43kWQ\nJI5o9/3umkn42gKx7RjFWjPagV2+fPlyoQfef/99fPrppygoKEBVVRUeeugh3H///aioqMDevXsx\nceJEn97QaPTwhHkY0GgU3ep6nImE66RrDH8ajQZJSQk4ffpcp31jOw8cr1bjeJUGZ2+qcKVBAaNV\nggS1Ff6qDZKgtsJqZ2C2SWCzt5SCTI02Y6De6Lf3cFAoZDCbO2eOC12nySZBsrbrNwJWDjhaqcbJ\nG+4/Q6lUirS0dOTn5UOrluPwoR9QXLwPRqMR/foNaP35sAyPKw0K2OydB6eR88jSmbxKHvPXz9lu\n5zBo0BC/Fo7RaISrjDmdCa9fvx6fffYZVCoVXnvtNRQWFmLBggXgeR4zZ87028AIIcRfMjIyMGRI\nDo4ePdSuSYS/Ske6IoYl0EBcpyNxqrxeAY6/dUGevDbDMBg3bgIyMjLx97//BVu3bsa5c2dx553z\nMGDAQLASxukKQkos5/XnF4yfs785vUSGYaBStQx+//79mDBhQuvXCSFErAYMGICsrD7gOA5A8Is3\n+HL+2B88vU53R4E6cgS2tgHY2Ws7k5KSiqef/i3y88fh8uVLWLnydbz44nMoLt6HfnENgke4RqR5\nVxQq3Ip0ODidCbMsi4aGBhiNRpw8eRLjxo0DAFy5ciVkbcgIIcQTubmjIJPJUVFRjus1TT53Tgon\n7jK0jVYJLtUrvaqn7Lrox63XrjVJoVPZXN54KBQK3H///8HkyYX47rsNKCk5iA8++Dvi1q3B/Pl3\nY1JOLswc27qCIGE6z45d6UqHrFByGk0feugh3HnnnbDZbJg/fz4SExPxzTff4I033sDSpUt9fkOd\nTg2pVBzHCvxBr3ff2qs7iITrpGvsHhzXWFQ0HgBQXVOLo2/uRr2x8y9gjZyHXqdAOP5Kio5uH6TU\nHKC5ysNg6RxRNXIeV41RKKu7lVDbehRILsXIdCtsHGCyMlDJ+NbPo7GZcVn0w2HflWho5DxSYjmM\nSLO63McdMmQAhgwZgMrKSnz77bfYtm0b/va3PyM3Nxf/9V//BV1snNNrdMXd9Xvzc7ZYWOj12qCs\n/DK8iyKtlZWVqK2tRf/+/QEAO3bsgFKpxJgxY3x+w+rqRp+fKzZ6vbZbXY8zkXCddI3dg7Nr/Nfm\nM9h8sKLT1zNiTaLdK3QlOlqFhobO3eyOVakF91fTY0yoNsgFj1ApWQ49oiyoMnSeIfM8sL0sVvB5\nznj7mVZWXsc///kxzp49DaVShbvumocJEyYhNlYjeI2uOLt+b8dktVpw9933+jUIO7sBdpodDQBR\nUVFISEho/XNGRgZSU1O7NJDulJ3Z3bNNHSLhOukauwdn1zgwQweT2Yb6JguazTYopbaAZS4Hg7Ps\naGcZ2hmxzThXq4JQIwcbz6DeLPsxQ5mBzS5BXbMMVjuDpCgrjNaWP3fGC76e2SZBr5hmj7Oao6Ki\nkJeXj9hYHU6dOoHS0h9QUVGO3NwR8LZaqb8y1IOZHe1yJhwI3elOPBJmFkBkXCddY/fg7hrNVg71\nTWY01Fbi0A/7YbXaIJPJRFkK0xVnM2GHjmd5OburGa1wMFVJOUzOqAPDtC8rqWTtiFHaUGmQCz4P\n4FGQUefT/mt9fR0++ODvOH36FNLS0vDTny5FQoLe69fp6nnwYM6EKcOKEBIxFDIWiTo1EnWZ6JWW\niqamRtTX18NoNKKpqQmXLpWBFWtpJS90rCjlqpiIM22TmToevQKcB/WuVAiLiYnFY489ji+++Azb\nt2/FSy/9AQ899HP069ffq9dxV1FLTML/bxshhPhAJpNBp4tDRkYmBg4chNGjx2DSpAKwLBvQfsah\nItQLOT3G5HHHobZHrzztaOQLlpXiJz9ZhAcffBAmkwlvvfU61q5dJdi3uDugmTAhhPwoMTERM2bM\nxu7dO1FTc0M0DSL8wVkxkWNVEJwh6zWWkDZJKCoqQmxsPN5772/YuPFbbN26BRMnTsLUqdMRExPr\nl/cQA5eJWYHQnRJDIiHRBYiM66Rr7B78cY0syyIzszdsNhuqq6tEuV/sLDHLExIGkLN8a+JUx2Qm\nJWuHSsahySJ1WfqRYYBEjRW9YpqRFm1GnzgTkqL8VwZUoZAhKioGkyYVIDo6BuXll3DixHFs374V\ndXV16NEjCRpNlH/erANRlK0khJBIlpMzDGazGZcuXQDLdt9flR1nyBdqlbhU73npx0Dvv8rlchQU\nFGHChEnYu/d7bNz4DXbu3I5du3Zg2LARmDp1OjIzewfs/QOt+/7NIoSQLho1ajQMhibcvFnjl1lR\noLo4+QMraWnJWGUQrpB1rVGO7Dgj5CGKGlKpFBMmTMLYseNRWlqCTZs2oLS0BKWlJRg6dBgWL16C\n6OiY0AyuCygIE0KIEwzDYOLEydi48Rs0Nzf7/DqOJggd90/7xRth4cQTlF2VfjRzEuy4FIueWtel\nLgONZVmMHDkaubmjcPr0KfznP1/iyJFDeP7587jvvgcwdOiw0AzMRyL4sRNCiHixLIuCgimQSHxP\n0nI0QWg50sO0LvFuvqDDtrJYbC+LxbEqNewhTspWSu1Os6UBBmauZdwnqtVBHZfgaBgG/fsPwOOP\nP4kFC34Ck8mEd999G59++jHMZnOoh+cxmgkTQogbKpUKkycX4PLly7DZrLBareA4DlVVlbDZbC6T\nt1w1QXB0JhJLyz1PzxNXNskxIMEY1Nm7s6V8iUSCoqLb0L//AHzwwd+xa9cOHD16BFOnTsf48RMg\nlwsnRIkFBWFCCPFAbKwOsbG6dl+z2+04efIEzp49A5vNKhiMXS3xdhSK4NaR44jRtUY5zFxLKcuO\ngtmVyNlSfscl8ZSUVDz11DP4z3++wtat3+Hzz/+Nb7/9D267bRomTpwEpdK7rkzBQsvRhBDiI4lE\ngkGDBmPOnDsxcOBgtJSAbM/1Em97juAWSo5s6Um96qBgPSvkEUjOlvIdS+Jt+yPLZDLceedcvPDC\nq5g+fRasVgvWrPkCzz77DM6ePROU8XqLgjAhhHSRRCLBwIGDMG3aLKhUatjtnUtGeiKYwc0duRTo\nqQ1MVSxPuVrKv94ox9FKNbaXxXbaV9dqta3BeMaMWWhsbMAbb6zAxo3ftPvZiAEFYUII8RONRoNp\n02YgPb0XbLZbxTQ6loxkGeFAEKzg5imhUpcZsSa/VcVyxcYBtSap06X8Zk6CS/XOZ8hAy8/jjjvm\n4vHHfwWtVou1a1fjz39+G01NTQEfv6eoi1IXREJXGiAyrpOusXsQ0zVevHgBJSUH2u0TO5KL5Kwd\np2vc73M6466Lkr8F8nxzx9d27AFXGxUwWBwfhnC3JlfdnzqOs6GhAR988HecOnUCMTExmDp1BiZM\nmCiYuEVdlAghJMxlZvZGdHQMdu7cDp5vmfm2rS4lVMdZrAJRFctZwhUgXMvaU86SxqKjo/HYY49j\nw4ZvsHHjN/jii//Fpu++Q8FtszEhfzQ0amWXrsdXFIQJISRA4uPjMXXqdGzbtgVmc3OnmVU4tdzz\nN0fClYNjOdnZUr0j6U0ltUOvsaDaIPe6laJEIsHMmbMxfsIkbD/WiGZZD1ii4vHNiRqouOuYPDga\nUZrgnoEW8b0XIYSEP8c+sVYbLbqkoFBxfXba+RJwXkoDJmfUYWgPY5daKV429YAkfgjU0YmQSFio\nohMB3SB8/M1RfPXVehiNBo+vpasoCBNCSIDJZDJMmTIV2dn9odFEwWy2dMuexZ7y5uy0g0pqh05l\naw2wviaNuboB0Gfm4tuNG/D73z+L6upqr8bnK1qOJoSQIJBIJBg6NAdDh+bAaDTi3LmzuHSpDGZz\nsyjbJQaS4+y00HIyy/CCs+GOM1xn/ZHdcXUDoIpOxB1zF+HSuSNQKoNTaYuCMCGEBJlarcbQoTkY\nNGgw9u79HlevXoFU6vmvYxvXUqBC7Aldzrgqj5kabQbDoDU7um3muLPX8mZf3dUNgEpqx+TJ42Gf\nMBpabbTnF9QFFIQJISREWJbF+PETcfr0KRw5cggs67pJRPvjOyqvjzaJiSOoOjumpdbwqK41+/1G\nw9UNgGO2bef8937uUBAmhJAQ69evP+LjE3Dw4H5wHAeAAcMw4DgbjEYjZDIZAOcZxUBoGz/4wt1y\nspQNXOa4qxuAYKMgTAghIpCQkIDp02d1+vr169dw+vQpVFy55jShSAyNH3wVimNavu4nBwIFYUII\nEbGkpJ5ISuqJi1duYMsnRwS/J5hdjboTMZzTDsP7JkIIiTzJiTrERwtn7Iqp8QPxDgVhQggJAwoZ\ni+HZesHHxNb4gXiOlqMJISRMLCzsAwA4cr4G1bUmKKU29IiyhiShiPgHBWFCCAkTrESCe6dk46fz\nVDhfVoP6mms4cvgHOOsoJGaB7MwUTigIE0JImFHKpUjUqZGoy0J6WioOHNjvdcGPUHHWPSkczzr7\nQwTffxBCSPhTKBQYP34iJkyYBJlMLvomEY6zzi0Vq5jWs84nqgPfvYizt1Qa40T0ETF8kKuI22wc\npFLXVWEIIYR4z263o7i4GGfPnhXlrNjGAf85roTB0nn+p5HbMWtQMwIRHuw88EO5DFfqWBgsDDRy\nHimxHEakWQVn3xaLBUuWLOnUejIQgv5Tqq3tPgkEer0W1dWNoR5GwEXCddI1dg+RcI2A6+vMyhoE\nrVaP/fv3wWw2iao5hMEigcHSuVxky2MMqmvNred2o6NVaGgw+eV9j1WpUVYna/deZ6oksFhsgpXG\nrFYLqqsb/RqE9Xqt4NfF89MhhBDiF4mJiZg1azYyMjJhswWxELIbjuYJQgJ11tlV68LKJnnIl6Yp\nCBNCSDckkUgwcuRojBo1WjT7xI7mCUICddbZVetCR6WxUKIgTAgh3VhmZm9MmlQIsRxhGqg3IiPW\nBJWUA8BDJeWQEWvCQL0xIIlToZh9e0N8O/eEEEL8Sq/XY9q0Gdi2bTPMZnNQEo6cEWqewDCdjy2l\nxdnRJ8bU5WNLnrQuDCWaCRNCSATQaDSYPn0WevfOAsuysNlsIR2Po3kCKxE+tnSmSubTsSWh2bSr\n2Xeo0UyYEEIihFQqxfDhuRg2bATKyspw4cI51NTcCOlxJneJU562aHRXBEQsrQs7oiBMCCERhmEY\nZGZmIjMzExUV5ThwYD+CXDKilSeJU560G3TMpm89l239s+MYkhhaF3YkknsBQgghoZCamoaZM2+H\nThcHjgv+cSZ/JE6J/RiSKxSECSEkwikUChQUFGHo0GEBPc4ktF/rj2NLYj+G5AotRxNCCAEA9OvX\nH4mJidi9excsFv9lUbvbr3UkSHXOjvYsccoxm25J7GpPDMeQXKEgTAghpJVOF4eZM2dj797v/daZ\nyd1+rVDilC5WhYYGz15f7MeQXBHx0AghhDiYrRyqao0wWwO/b8uyLMaPn4icnOGw2WxdStryZr+2\n7bElb4n5GJIrNBMmhBAR4+x2fLb1HErPVONmgxlx0QqMy0nB7fnpYAPcnKFfv/7IyMjE6dOncPVq\nBerr6yCTCQdUZ/yV/eyOmI8huUJBmBBCROyzreew+WBF659rGsz4ctcFGE0W3DslO+Dvr1AoMHRo\nDoYOzUF9fR1KSg7i5s0aj7szBXu/VozHkFwJg/sEQgiJTGYrh9Iz1YKPlZ65EZSl6bZiYmIxeXIh\nNBqNx88JRdOGcBLhl08IIeJV32TGzQaz4GO1jc2obxJ+LJAkEgkmTiyANw0hwnW/NhhoOZoQQkQq\nJkqBuGgFagQCsU6rREyUIgSjaqlDPXbsOOzatQMs23mZuaNw3a8NBvoYCCFEpBQyFsOz9YKPDc9O\ngELmPgAGSlJSTwwePFSwEYSzloRdyX7urmgmTAghIrawsA+Alj3g2sZm6LRKjMtJxu356SEeGTBg\nwEA0N5tw7tzZlsIejMRlUQ7SGQVhQggRMVYiwb1TsjFvUhbqm8yIiVIgNTkW1dWNoR4aAGD48FwM\nGZKDs2dPY/2eCpTVyVofE2qiQNqjRQFCCAkDChmLRJ06pEvQzkilUvTu0x815ijBxyub5LAEOZM7\nXFAQJoQQ0mWuMrnNHIu0zP7geT5kLRPFioIwIYSQLnNkcgvRaZUYOXwIZs++Az17Jgsmc0UqCsKE\nEEK6zJNMboVCgbFjx2PMmHzwfPhUtQokSswihBDiF0KZ3MOzE1q/7tCrVwbi4uKxa9cOGI0Gj0tg\ndkcUhAkhhPiFUCa3s0QyrVaL6dNnYv/+vSgvv+yXlonhKHJvPwghhASEp5ncEokE+fnjMHLkKNjt\nkbk8TUGYEEJISPXu3QdTpkyFXK6IuOxpCsKEEEJCLjZWhxkzZqFHjx4RlT1NQZgQQogosCyL8eMn\nYciQnIgJxBSECSGEiMqAAQMxZcoUeNMuMVxRECaEECI6PXv2xPTpM6FSqbt10hbDB3kX3GbjIJWK\nr/YpIYQQ8bHb7SgpKcGpU6eCdozJYrFgyZIlLZ2hAizoB7Nqa7tPJw29XiuaTiaBFAnXSdfYPUTC\nNQKRcZ1trzEjoz8SElJx8GAxqqoqAx6MrVYLqqsb/RqE9Xqt4NdpOZoQQojoRUVFYfLkQowfPxFS\nqazbHGWiIEwIISRsJCenYMaMWYiLi+8WGdQUhAkhhIQVmUyGgoIi9O8/EBwX3n2KI7NYJyGEkLA3\ndGgO9Ho9Dh36AXY7/+MeLgObzQqTyQiZTBbqIbpFQZgQQkjY6tkzGT17Jnf6+rVrV3HmzBlUVl4D\ny7JByXT2BQVhQggh3Y4jODc3N6Ok5ACuXq0Ay4ov5NGeMCGEkG5LqVRi7NjxiI/XizKjmoIwIYSQ\nbo1hGEycOBlyuSLUQ+mEgjAhhJBuTyqVYuLESaKbDVMQJoQQEhFiYmIxZkw+bDbxHGuiIEwIISRi\npKamITc3F2q1GmazJdTDoexoQgghkSUrqy+ysvqiqakJ58+fw/XrV1FXVxeSc8UUhAkhhESkqKgo\n5OQMQ07OMFy/fh3Hjx9FTc2NoI6BgjAhhJCIl5SUhKSkJNy4cQOHD5cG7X0pCBNCCCE/SkhIQFHR\nbUF7P0rMIoQQQkKEgjAhhBASIhSECSGEkBChIEwIIYSECAVhQgghJEQoCBNCCCEhQkGYEEIICREK\nwoQQQkiIUBAmhBBCQoSCMCGEEBIiFIQJIYSQEKEgTAghhIQIw/M8H+pBEEIIIZGIZsKEEEJIiFAQ\nJoQQQkKEgjAhhBASIhSECSGEkBChIEwIIYSECAVhQgghJEQoCBNCCCEhQkGYEEIICREKwoQQQkiI\nUBAmhBBCQoSCMCFe2LBhA+bOnYs5c+bg9ttvx3vvvdf62MqVK3Hw4EG/vE9hYSEqKiq6/PwtW7bg\nrbfe6vJ4nnrqKaxZs6bLr0MIaU8a6gEQEi4qKyvxyiuvYM2aNdDpdDAYDLjvvvuQmZmJoqIiHDhw\nAGPGjAn1MNspKipCUVFRqIdBCHGCgjAhHqqtrYXVakVzczMAQKPR4OWXX4ZCocC6detw7NgxPPPM\nM3jnnXdQX1+PN954A83NzWhoaMDTTz+NKVOm4KmnnkJUVBSOHz+OyspKLF26FPPmzUNdXR2efPJJ\nXL9+HVlZWTCbzQCApqYm/OY3v0FlZSWqqqqQn5+PF154AcXFxVixYgXsdjv69u2Lp59+WvD5a9as\nQXFxMR555BEsXbq09VouXryIZcuW4YEHHsCrr76K4uJicByHuXPn4oEHHgDP83j55Zexfft2JCYm\nguM4jB49ut3nUVFRgf/+7/+GTqeDUqnE22+/7XSsf/3rX6FUKnH+/Hn069cPr732GuRyOT7++GP8\n85//hFarRe/evZGeno5HH30UO3fuxMqVK2Gz2ZCamornn38eOp0uSD9pQoKIJ4R47He/+x0/cOBA\nft68efyrr77Knzx5svWxxYsX8/v27eN5nucfffRR/ty5czzP8/yePXv42bNn8zzP87/+9a/5pUuX\n8na7nT916hQ/evRonud5/rnnnuNff/11nud5vri4mM/OzubLy8v5r776in/33Xd5nud5s9nMT5ky\nhT969Ci/b98+Pjc3l29oaHD5/NWrV/O//vWv213Dpk2b+Llz5/LNzc38v/71L/7FF19sff3Fixfz\nBw4c4L/99lt+8eLFvMVi4Wtqavhx48bxq1evbvc65eXlre/D87zLsQ4bNoy/du0az3EcP2/ePH7L\nli38yZMn+alTp/KNjY18c3Mzv2DBAn7lypV8TU0NP2fOHL6uro7neZ7/97//zf/mN7/p8s+OEDGi\nmTAhXnjuuefw8MMPY/fu3di9ezfuvvtuvPbaa5g6dWq771uxYgW2bduGDRs24PDhwzAYDK2PjRs3\nDgzDIDs7G3V1dQCA4uJi/PGPfwQAjBo1CmlpaQCA2bNn48iRI/jHP/6BCxcuoK6uDkajEQCQmZkJ\nrVbr8vkdnTp1Ci+//DI++eQTKBQK7N27FydPnsS+ffsAAEajEadPn8b58+cxdepUyGRURRqGAAAD\nI0lEQVQyxMXFYeLEiYKvFx8fj9TUVLdj7du3L5KSkgAAWVlZqK+vx6VLl1BQUICoqCgAwKxZs9DQ\n0IDDhw/j2rVrWLJkCQDAbrcjJibGsx8QIWGGgjAhHtq+fTuMRiNmzpyJefPmYd68efj888+xatWq\nTkH43nvvxZgxYzBmzBjk5+fjV7/6VetjCoUCAMAwTOvXGIYB36a1N8uyAIBPPvkEGzduxN13342x\nY8fizJkzrd+nVCrdPr+tmzdv4rHHHsOLL76I5ORkAADHcXjyySdbx3/z5k1oNBq8+uqr7V5PKhX+\nVdF2DK7G6rjmtmOVSCSw2+2dXpPjOIwYMQJ/+ctfAABms7ndTQwh3QllRxPiIaVSiT/+8Y+tWcs8\nz+PkyZMYMGAAgJbAx3Ec6urqUFZWhmXLlmHixInYsmULOI5z+dr5+flYv349AODIkSO4fPkyAOD7\n77/HwoULMWfOHJjNZpw6dUowcDl7voPVasWyZctw3333tUsey8vLw+effw6r1QqDwYB7770Xhw4d\nQn5+Pr799ltYLBbU19dj165dbj8fT8fadsw7duxAU1MTLBYLNm3aBIZhkJOTg0OHDuHixYsAgHff\nfRevvvqq2/cnJBzRTJgQD+Xl5eGRRx7Bz372M1itVgDAhAkTWhOeJkyYgGeffRavvPIK5s+fj1mz\nZkEqlSIvLw/Nzc2tS7NCHnvsMTz11FOYNWsWevfu3bqcfP/992P58uX429/+hqioKAwfPhwVFRVI\nT0/36PkOGzZsQGlpKUwmE1avXg2e5zF27Fg88cQTuHTpEu666y7YbDbMnTu3NUgfPXoUs2fPRkJC\nArKystx+Pp6O1SE7OxtLlizBwoULoVarodPpoFAooNfr8eKLL+IXv/gF7HY7evTogRUrVrh9f0LC\nEcO3XXMihJAguXjxInbs2IEHHngAAPDzn/8cCxYsQGFhYWgHRkgQ0UyYEBISKSkprbNthmEwfvx4\nFBQUhHpYhAQVzYQJIYSQEKHELEIIISREKAgTQgghIUJBmBBCCAkRCsKEEEJIiFAQJoQQQkKEgjAh\nhBASIv8fNGpI0tK9t2oAAAAASUVORK5CYII=\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.6.1" }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 1 }