{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Importance Sampling and SMC\n", "\n", "This notebook demonstrates the importance sampling (IS) and Sequential Monte Carlo (SMC)\n", "tools in `psiphy.mcmc.importance_sampling`, applied to the built-in toy models.\n", "\n", "## What is importance sampling?\n", "\n", "We want to draw samples from a target $p(\\theta | x)$ but can only evaluate it up to a constant.\n", "IS draws samples from a tractable proposal $q(\\theta)$ and corrects with weights:\n", "\n", "$$w_i = \\frac{p(\\theta_i | x)}{q(\\theta_i)}, \\qquad \\theta_i \\sim q$$\n", "\n", "The **effective sample size** (ESS) measures how much information the weighted set carries:\n", "\n", "$$\\text{ESS} = \\frac{1}{\\sum_i w_i^2} \\in [1, N]$$\n", "\n", "Low ESS means the proposal is a poor match for the posterior — most weight falls on very few particles.\n", "\n", "**SMC** fixes this by annealing from the prior to the posterior through a sequence of tempered targets,\n", "resampling and applying MCMC moves when ESS drops too low." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import norm\n", "\n", "from psiphy.mcmc import (\n", " UniformPrior,\n", " eff_sample_size,\n", " importance_sampling,\n", " sequential_importance_sampling,\n", " SMC,\n", ")\n", "from psiphy.toy_models import GaussianSignal, NoisyLine\n", "\n", "np.random.seed(42)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 1 UniformPrior helper\n", "\n", "`UniformPrior` wraps a box prior and provides `sample(n)` and `log_prob(theta)`\n", "— the two methods expected by `importance_sampling` and `SMC`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ndim : 2\n", "sample: [[-1.25459881 4.7585001 ]\n", " [ 2.31993942 3.03342657]\n", " [-3.4398136 0.86437315]]\n", "log_prob in support : -3.891820298110627\n", "log_prob out of support: -inf\n" ] } ], "source": [ "prior = UniformPrior({'mean': (-5, 5), 'sigma': (0.1, 5)})\n", "print('ndim :', prior.ndim)\n", "print('sample:', prior.sample(3))\n", "print('log_prob in support :', prior.log_prob([0.0, 1.0]))\n", "print('log_prob out of support:', prior.log_prob([10.0, 1.0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 2 Importance Sampling — GaussianSignal\n", "\n", "Target: posterior $p(\\mu, \\sigma \\mid x)$ for a Gaussian signal with $N=40$ observations." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x_obs mean=1.250 std=0.734 (true: [1.5 0.8])\n" ] } ], "source": [ "gs = GaussianSignal(N=40)\n", "theta_true = np.array([1.5, 0.8])\n", "x_obs = gs.simulate(theta_true)\n", "\n", "print(f'x_obs mean={x_obs.mean():.3f} std={x_obs.std():.3f} (true: {theta_true})')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Broad proposal — ESS = 11.4 / 5000 (0.2%)\n" ] } ], "source": [ "# Broad uniform prior — deliberate mismatch to show ESS degradation\n", "prior_broad = UniformPrior([(-5, 5), (0.1, 5)])\n", "log_posterior = lambda t: gs.log_prob(t, x_obs) + prior_broad.log_prob(t)\n", "\n", "proposal_broad = prior_broad.sample(5000)\n", "s_broad, w_broad = importance_sampling(log_posterior, proposal_broad)\n", "ess_broad = eff_sample_size(w_broad)\n", "print(f'Broad proposal — ESS = {ess_broad:.1f} / 5000 ({100*ess_broad/5000:.1f}%)')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tight proposal — ESS = 114.1 / 5000 (2.3%)\n" ] } ], "source": [ "# Tighter prior centred near the posterior — higher ESS\n", "prior_tight = UniformPrior([(0.0, 3.0), (0.3, 2.0)])\n", "log_posterior_tight = lambda t: gs.log_prob(t, x_obs) + prior_tight.log_prob(t)\n", "\n", "proposal_tight = prior_tight.sample(5000)\n", "s_tight, w_tight = importance_sampling(log_posterior_tight, proposal_tight)\n", "ess_tight = eff_sample_size(w_tight)\n", "print(f'Tight proposal — ESS = {ess_tight:.1f} / 5000 ({100*ess_tight/5000:.1f}%)')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABEEAAAGYCAYAAACtaPyUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjRpJREFUeJzt3Qd4U9X7wPG3e1FG2XtvkCUKyBARFBD3HqDiQMWFOMA98e8EFXFP3IAKggo/ZYMKshQEkb336IDO/J/3lBvSkpamTZpxv5/nSZMmNzc3Jzc597z3PeeEORwOhwAAAAAAAIS4cH9vAAAAAAAAQGkgCAIAAAAAAGyBIAgAAAAAALAFgiAAAAAAAMAWCIIAAAAAAABbIAgCAAAAAABsgSAIAAAAAACwBYIgAAAAAADAFgiCAAAAAAAAWyAIAgAh5qOPPpKwsDBZvHixBKvPP/9cRo8eLaHC4XDIl19+Kd26dZMqVapIbGys1KpVS8455xx57733JNDVq1dPrr/+euf/GzduNPuY7mt299xzz8l3333nk3VTzgAAeB9BEABAwAm1IMiIESPkqquukubNm5ugx48//ijPPPOMVK1aVb7//nsJNtWrV5eFCxdK//79xe58GQShnAEA8L5IH6wTAIBiSUtLk/j4+JAqvSNHjpiAzsCBA+Wdd97J85hmV+Tk5EiwiYmJkU6dOvl7M0JWdna2ZGVleb2cdV/ULCTN4gEAwK7IBAEAG9DGdpkyZWT16tWmC0ZCQoI5y/z888+bx3/77Tfp2rWrub9Jkyby8ccfu+1iM2PGDLnhhhskKSnJLDtgwABZv379Ca/3wQcfSJs2bUyDS5e96KKL5J9//nG7TX/99Zf06dNHEhMTpVevXnLmmWfK1KlTZdOmTeY1rYvlySeflNNPP92st2zZstK+fXt5//33TZeT/F04zjvvPPnpp5/MMnFxcdKsWTOzbflt27ZNbrnlFqldu7ZER0dLjRo15NJLL5Vdu3Y5lzl8+LAMHz5c6tevb5apWbOm3HPPPZKamlpo2evj6enpprzdCQ/PWxV7+v5++OEHadeunXl/mmmi/1ufmf6vn9Npp512Qvcoq/xXrlxpyl2Xq1y5sgwdOtQEozztpvHEE0+Y+3R9mvVSrlw5k+ly4403yqFDh/I8/+DBgzJ48GDzHnUbNKNE9yN9vq6nOKzXX7p0qVx88cWm7HQbrr32WtmzZ0+eZTXw9MILL5j9QQMN2kVJg1Rbt27Ns5yuS8tYH9fldL/QbbWW09fTz1e/L9Z+qvuvZefOnXLrrbeark+6z+i+o5+vBjjyl6Vuj2YH6TL6WjNnziywO8y8efPMZ6bfGQ0adunSxXxn3H1np0+fbj4D/Wx1Wd0XC7J582ZTXtb71f3n5ZdfzhOos7bppZdekldeecVsr36GnTt3Nr8jJ2Nt16+//io333yzVKxY0XxWWv5allpml19+uZQvX958Z/Q7l5mZmWcdGRkZpqysz0/fm/4u5f+cv/rqK/Pbouuxvh8PPfTQCd9Z67vw33//Sb9+/cxt/S247777Ci0vAEBwIhMEAGxCGxLaOBwyZIjcf//9psuJdtPQxv3EiRPlwQcfNI21119/3TQKWrVqJR06dMizDm249u7d2zx3y5Yt8sgjj5hG34oVK0yjRY0aNUpGjhxpGsJ6e9++faaBqo2kRYsWSePGjfM0Zs4//3zTUNTGiTYOdRs0ILFu3Tr59ttvT3gf2gjT5evUqWP+14bXnXfeaQIZjz32WJ5lly9fbhoyum5tkGtXFH0PjRo1ku7du5tl9HkdO3Y05aPbfcopp5ht/vnnn+XAgQPmeRoU6NGjh2n8WstoY19fT4M4//vf/wo8u16pUiXzem+++aZpXGojq2nTpgUu7+n708/w4YcfNg1+bWDrZ6z3/fLLL6arhr6OfrbamN+wYYNpDLruE7o9VvkvWLDANC41ADVlyhQpjksuuUSuuOIKU85aNrotygo+aYNag2calNH9QoM82rXm3HPPdbs+3X4t+1mzZhXp9TXgpo1o3c/1M3r00Udl1apV8vvvv0tUVJRZ5rbbbjNZORrw0XLRMtfl9DWWLFliPjNtKOu+ro38sWPHmv1AG+ganEhOTjbr0e0+66yzpGfPnub5Shv0SpfV4JMGufRza9iwoVley1df78MPP8yz3a+99poJQGpwQdfh+j1xNXv2bLNdug9qcEyDALpvaZl+8cUXpuxdaQBEAzeffvqpeU9WGeSnAQQNpuh38umnnzZBNg2oaRBCv4v6Gq60TDQIYXVb0/ev+5LuY7ovnsxNN91k9lUdK0eDTfq90u//mjVrzP36G6Dfq//7v/8zwadhw4Y5958LLrhA5s6dKw888IDZZt1fH3/8cfNbpPuVtY+vXbvWbJMGKzXIp0FgXd8ff/xhgjCu9Lugv0W63+pvxpw5c0w56HvJ/70DAAQ5BwAgpHz44YeaMuBYtGiR875BgwaZ+yZOnOi8LzMz01G5cmVz/5IlS5z379u3zxEREeEYNmzYCeu86KKL8rzW/Pnzzf3PPPOM+f/AgQOOuLg4R79+/fIst3nzZkdMTIzj6quvPmGbPvjggxPeQ//+/R1169Y96XvNzs427+Opp55yVKxY0ZGTk+N8TJ8fGxvr2LRpk/O+I0eOOJKSkhy33nqr874bb7zRERUV5Vi1alWBrzNq1ChHeHh4njJVEyZMMO9h2rRphW7nH3/84ahTp45ZVi+JiYmO8847z/HJJ5/k2WZP35+W9datW533LVu2zKy/evXqjtTUVOf93333nbl/8uTJJ5T/mDFj8rzms88+a+6fN29entfS5S0bNmwwy+h+YXn88cfNfS+88EKe9d1+++3mc7C2ferUqWa5cePGnVDGer+ux5Xui2eddVaBZZT/9e+9994893/22Wfm/vHjx5v///nnH/O/bper33//3dw/cuRI8//ixYvN/1p2hUlISMhTNhbdx8qUKZNn/1MvvfSSWe/KlSvzlGXDhg0dGRkZeZZ1V86dOnVyVKlSxZGcnOy8Lysry9GqVStHrVq1nOVsfWcHDhzoKIqHHnrILK/l4Oq2225zhIWFOdasWZNnm1q3bm1e13Uf1/u/+OKLQl/H2q4777wzz/0XXnihuf+VV17Jc3/btm0d7du3d/6v68//W6b0u6n3v/nmm25fV8tFv0uzZ882yy1fvvyE78LXX3+d5zn6O9a0adNC3w8AIPjQHQYAbELPqOtZUUtkZKTJUNBUce1OYdEuCpqxoGdX87vmmmvy/K9nYevWrWvOjis9063jDrjOJKI0tVzPmGt2grvMAU/oGdyzzz7bnKGNiIgwZ7b1TK1mb+zevTvPsm3btnVmVCjtnqNn213fmw5SqmfyNVW+IHpGXDNjdH16ttq6aNciLdeTZSlopomm2mvXHD3jrVkxWhbaBUDPPrt2dfH0/Wm3HIv1HvSMuOvYKtb9RflMr776anNtfaae0vfjSjMWjh496tx2zWRQmq3hSjOH3NFydrffFCT/+9HX0X3dej/Wdf59VLM2tJys19LvRoUKFUwWzVtvvWWySTyh+4zuV5rF4LrP9O3bN085uJZbQVkaFs3k0IwW7aqlXTYsup9cd911JlNJMymK8/3S/a5FixamHFxpOen+mT9zQrNL9HVdP+eC9jF3NAPHlbWP5h9sV+93XaeWq2adaeaLa7nqd6FatWp5vovaxUr3Z73f+i5pVpHK3z1Pv8e6Tlf6nor6fgAAwYMgCADYhDaKNQjgSscp0KBHfnq/Nlzz08aEu/u0ga6sa3fjX2hj0HrcdZus7gNFoWns2sdfvfvuuzJ//nzTxUa7gygNwLjS8Qby0+4DrstpNwDtglMYHRtEu/xoI8r1omMyaANx7969J912XV6DJs8++6zpaqPdiTRYoY06DcQU5/3l/+z0cyvs/vyfqQYH8peR9Rnn/6yKKv/6tLxdt13Xq6+bfxu1u4k35N9Hrffo6T6qQSgNVGjjWgNXLVu2NI9rt4v8Y1QUtM9ol6L8+4yuR+XfZwoaM8aVds/S/a2gbXd9f56s13qeJ+s92ed8Mp7su677rZarjimj9+cvW+2CZJVrSkqKmZJag0baBUmDI/pdmjRpktvtdPf7qO/J3e8gACC4MSYIAKDItJHh7j49a+7aMNqxY8cJy23fvt2MteDK01kqdPwAbexo4MC1wVKSKUp1UMX8A2Lmp9ut4wy4G1TVetxTWlY6VoE2zv7++2+TpeOL91cYPYOujVvXBq31GbsLIHmDrldfd//+/XkavO72reLQ9bhmx+R/j677aP7gV/59tHXr1uYz0cCDBsF0UM+nnnrK7As6hkphdD2aSaBBL3es4IIn3wXNTNExRgr6flmv6+l6rXLxZL3+otuh26pZVe5oYFJp5opuu36/rOwPpQEUAIC9kQkCACiyzz77LM//OpCmpotbM2JoNw9tII4fPz7Pchpk0EaJzmhRFPmzNVwbdHpm3zUNX5fTQR+LS7snaBeJ/N0I8qfu6+CQ2vg69dRTT7joIJIF0ayBgrIqrJR8q0Hsi/fn6Weqg94q11lOvMlqkOrMHa402OCL9/P111+bQIj1frRblsq/j2qWgH4e7vZR/Vx0tqNXX33VdMXQwVNPtq/qPqPBLR0Q1d0+kz8IUhQ6uKfOHKTZDK6vqYOF6vvRoI529yoOfd/a5cf1valPPvnEvH/t2hMItFz1+6TTCLsrVx102DX4Y2WoWN5++22/bDcAIHCQCQIAKDKdeUFndbjssstMdw7tpqFn3W+//XbzuDYQdZYI7T6g413oOA/aYNFZSzSzQbsSFIWegdeG3rhx48wMNXr2Wxs4Ol6ATsup/fx19ghdt86mkb+h4wk9s6/dUXS2GN1ufW09W6xnmnVGCp0BQzM2dAYdXebee+81Z/i14alTiuoUpDqbhDZO3dHpYTVIomWmY33o+Ciaqq9nqMeMGWPGPNDZMJQv3l9htEuBToGq26Pjllizw2hgSKdM9gWdBeaMM84wZaYzE+nnq2PJaGPb3ZTBGhTSwElRxwXR/UafozOoWLPDaADDGoNEG8latjoLkr6Wvldrdhj9bPTzVZqNozOiXHjhhdKgQQOTDaLr1n1D123R/UU/S+36ot1JNBNBX0P3K51SWsfNueuuu8x92rVCX2vatGlmnJGTdcNyR2dc0tfXoITO3KKfoW6nBlx0dhhPs6ss+r71M9B9ULddx/rRaXd13TqbTnGDK9525ZVXmkCXZk7dfffdZgwTzZ7SQKsGM3XmGJ0hSMtdM2d0liD93dFl9Hk6oxIAwN4IggAAikyn5NSsBG2IpKenm4aYNuRduzXolKg6sKpO+aln+zUzRM/C63StBU37mZ82brQBq0EJDSJoA1QvehZfu6ToNJc6iKEGYG6++Wbzejq1ZXHoOnQsDm0oPf/88ybwoF1kNAhgvS89A69TcurjOrWqNdWsDrqqgY3CMkF0zBMNAmkjXt+PjmmgDVWdelWDKzrwpjWIqS/eX2GsrjfaSNfgh74nfb0XX3xRfEUDDxow0CCIlqdOyapBEc1k6NSpk3OqZYue8ddLUWmgQqfe1QCaNdilTuNqjTmh9DHN0ND9Wad61fE/NDijAQaru4zuq7otL7zwgulWoc/XQIZ2iRk0aJBzXbr/33HHHeY7YU2lrEERDYho0FCnWdXy1Ea6Bkj0c9fX0gZ6cej6NatK91cdtFSDcRrkmTx58gmDjXpC93kNgun315o6W4M/+v6t6WkDgWZJ6XvVctffIv3MNOilASUtGw1KKf0cNYij+9m1115rvsMaINHfJJ2WGQBgX2E6RYy/NwIAENi04XfDDTeYLgOakYHgpw3oCRMmmCyQQKDdcHRmFx0MVs/ie0oDHxps0oFuA2X8CgAAEHjIBAEAAKVKu21s27bNnLXXzJDffvvNZEtod6PiBEAAAACKiiAIAAAoVdotRAdC1S44qamppuuIZqbo/wAAAL5EdxgAAAAAAGALTJELAAAAAABsgSAIAAAAAACwBYIgAAAAAADAFgiCAAAAAAAAWyAIAgAAAAAAbIEgCAAAAAAAsAWCIAAAAAAAwBYIggAAAAAAAFsgCAIAAAAAAGyBIAgAAAAAALAFgiAAAAAAAMAWCIIAAAAAAABbIAgCAAAAAABsgSAIAAAAAACwBYIgAAAAAADAFgiCAAAAAAAAWyAIAgAAAAAAbIEgCAAAAAAAsAWCIAAAAAAAwBYIggAAAAAAAFsgCAIAAAAAAGyBIAgAAAAAALAFgiAAAAAAAMAWCIIAAAAAAABbIAgCAAAAAABsgSAIAAAAAACwBYIgAAAAAADAFgiCAAAAAAAAWyAIAgAAAAAAbIEgCAAAAAAAsAWCIAAAAAAAwBYIggAAAAAAAFsgCAIAAAAAAGyBIAjgJx999JGEhYXJ4sWL89z/888/S58+faRGjRoSExNjrs8880x5/vnn/fZZ6eu3atVKgsGsWbNMuep1aahXr55cf/31pfJaAADY9Vjk888/l9GjR7t9TN/DE088Uex6/LzzzjvpcqtWrTKvsXHjRilNHGcA3kcQBAggb731lpx77rlStmxZeeONN8xByP/93/9J8+bNZcKECf7ePLjx7bffyqOPPkrZAABCQqAeixQWBFm4cKHcdNNNPn19DYI8+eSTpR4E4TgD8L5IH6wTQDGNGjVKunfvfsJBxnXXXSc5OTlBWa7Z2dmSlZVlziSFkiNHjkhcXJy0a9fOa+sM1bICAASPYDwW6dSpk4QajjMA3yETBAgg+/btk+rVq7t9LDzc/1/XuXPnmgMNbfzXrFnTZEBow92iZ0c0JfWFF16QZ555RurXr28a9DNnzjSPT548WTp37izx8fGSmJgovXv3NmdvXP33339yww03SOPGjc1y+joDBgyQv/7664TtWb16tTlbpctVqlRJhgwZIsnJyUV6L5rSqtu6dOlSufjii80Zr3Llysm1114re/bscZsqO2nSJBP0iI2NNWeDCuoOs3nzZrOeKlWqmPevZ89efvnlPAePJysrAAD8IRCPRbQrztSpU2XTpk2m7rQuhXWHmTdvnjnm0DrbOmZ57733zLLusjl++uknad++vTnGadasmXzwwQd5ug1ddtll5nbPnj2dr6/3F4TjDCBwkQkCBBCtrCdOnGgqzosuusiMwxERESGBYOfOnXLllVfKQw89JE899ZQ5GNHG+4EDB0y6rKvXXntNmjRpIi+99JIJLmhAQ9NYr7nmGtPH+IsvvpD09HQTANADm19++UW6du1qnrt9+3apWLGi6XdcuXJl2b9/v3z88cdy+umnm4BF06ZNzXK7du2SHj16SFRUlLz55ptStWpV+eyzz2To0KEevS8t58svv9wEUFauXGkOkjTl9ffffzfrtixZskT++ecfeeSRR0zAIiEhwe36NIDSpUsXycjIkKefftoESX744QcZPny4rFu3zmzrycoKAAB/CcRjEa07b7nlFlOPaveQk1mxYoU50aL1qx5D6MkS7eYzfvx4t8svX75c7rvvPnOMo8cTGiwZPHiwNGrUyGTF9O/fX5577jkZOXKkjB071gRLVMOGDU+6LRxnAAHIAcAvPvzwQ4d+BRctWuS877///nO0atXK3K+XuLg4R69evRxvvPGGIyMjw2+fVI8ePcz2fP/993nuv/nmmx3h4eGOTZs2mf83bNhglmvYsGGe7c3OznbUqFHD0bp1a3Pbkpyc7KhSpYqjS5cuBb52VlaWWVfjxo0d9957r/P+Bx980BEWFuZYtmxZnuV79+5ttmHmzJmFvqfHH3/cLOe6TvXZZ5+Z+8ePH++8r27duo6IiAjHmjVrTliPPjZo0CDn/w899JB5/u+//55nudtuu81sr7WOgsoKAIDSEkzHIv379zd1rju6nVqvWy677DJHQkKCY8+ePc779PijRYsWZlmtgy26ztjYWOexjDpy5IgjKSnJceuttzrv++abb4p0fGHhOAMIXP7PrwfgpGcU9GzE7NmzTXeLs88+WxYtWmSyG/TMzNGjRwssLT0G0PEkinMpSh9f7b5y/vnn57nv6quvNs+dM2dOnvt1OdcsijVr1pgMD+1P7JpKW6ZMGbnkkkvkt99+k7S0NHOfbo+ebWnRooVER0dLZGSkuV67dq3JxLBot5GWLVtKmzZtTtgmT2h2iivNCtHXzN8t5ZRTTjFnlE7m119/Ndt+2mmn5blfu8zoZ6SPF1ZWAAD4UyAfixSVbvtZZ51luspa9PhD63h32rZtK3Xq1HH+r11otM7X7jclxXEGEHgIggABRitpTb187LHHzBgaGjy44oor5M8//8zTPzU/TffUxnRxLjfeeONJt0vTQ/OrVq2as/+wq/x9ia3H3fUx1mn39MBHu9WoYcOGmS4pF154oUyZMsV0S9GDLw126CBhruu0Xt/dNhVV/uU1AKLdcU72njztS63v03q8OOsFAMDuxyJFpXWtu+MWd/cprffz03G6XI87iovjDCDwMCYIEOB07IkRI0bIV199JX///XeBy+ngoRosKA7XMyUF0TE43I0T4u7gwXWwMtfHd+zYccI69MBKD7YqVKhg/tf+ugMHDjTZIK727t0r5cuXz7NO6/XdbVNR6fI6YJpFz0bpwdPJ3lNB9HkFvU93ZV3U9QIAYPdjkaLSuriw45bSxHEGEHgIggABRBvP7jIDrG4gVjZBQRW+uzMZ3qKzrujZINcuMTrYqXW2qDA6mKkGGnR5HSDUavinpqaawdesGWOUPpZ/ilgdhHXbtm1mgDKLjs6uA6tqyq5rlxh9DU/oYKodOnRw/v/111+bQIgO2FocvXr1MtML6kCq1sBp6pNPPjHvTbcbAIBAFajHIp5kZujA6dOmTTMnUKzgimadfvPNNyV6feVpdgjHGUDgIQgCBBAd40Ib0X379jV9crXfrXYH0elVNYVTRyr3Fz2oue2228z0r9pPVg8u3n33XXOfaz9adzRQogEL7RerU83eeuutZnaYF198UQ4ePGhmgrHo4zrlnE5Pp+NwaOqtLlerVq0867znnntMSq6O2K6z1Fizw+i0uZ7QaW+1C4yOIm/NDqNBlYL6DZ/MvffeawIeul06i07dunVNEEdHtteyKsq4IgAA+EugHou0bt3a1Nnjxo0zJy/02OLUU091u+zDDz9sutTq+9DbOu2tzg6jJ1+KO9WvzpKj3nnnHTNOmo4borPFnSzow3EGEHgYEwQIIBoM0DMVzz77rPTr189kXWiDWgf7XLx4sV/Hj9A+rZplof19dbs0Y0KnitMpXotC38N3331nuppov+IbbrjBTAmrA5Ba0+OqMWPGyLXXXmuyKTStVrNP9AAi/zR0uj068JkOQqrBBX2OHpDkn673ZHTdGji5+OKLTd9nfc3p06ebwViLQ6f1XbBggRmQTVOHNajz888/myDQ66+/Xqx1AgBg92ORu+++Wy699FJz7NGpUyfp2LFjgcvqyYwZM2aY4Id2sdXpdTW4c/vtt5vHy5Ur5/Hra8Bj9OjRJgNVs0X19TXQcjIcZwCBJ0yniPH3RgBAaXviiSfMqPd79uzxaj9kAAAQmPr06SMbN26Uf//91+evxXEGELjoDgMAAAAgpOhsc+3atZPatWvL/v37TZdZzQ55//33/b1pAPyMIAgAAACAkJKdnW26uersLDowuXaf/fTTT033WQD2RncYAAAAAABgCwyMCgAAAAAAbIEgCAAAAAAAsIWgGBNEp+navn27mZNb+/QBAIDApxPQJScnS40aNSQ8PHjOu3DcAQBA6B53BEUQRAMgOrIzAAAIPlu2bJFatWpJsOC4AwCA0D3uCIogiGaAWG+mbNmy/t4cAABQBIcPHzYnMax6PFhw3AEAQOged3gUBBk3bpy5bNy40fzfsmVLM/VU37593S4/a9Ys6dmz5wn3//PPP9KsWbMiv67VBUYDIARBAAAILiXpyuqPYw+OOwAACN3jDo+CIJpS8vzzz0ujRo3M/x9//LFccMEFsnTpUnNQUpA1a9bkCV5UrlzZk5cFgNCXni4ydGju7TfeEImJ8fcWAQGBYw8AXkNdC0CDJA4dPaQEkpKS5MUXX5TBgwcXeDbmwIEDUr58+SKvMz093Vzyp7UcOnSITBAAoSk1VaRMmdzbKSkiCQn+3iKgxLT+LleunNfrb28fe3DcAdgEdS0Q0op63FHsMUGys7Plm2++kdTUVOncuXOhy7Zr106OHj0qLVq0kEceecRtmqqrUaNGyZNPPlmk7dAYTlZWltkeIFhERUVJRESEvzcDgSQqSuSZZ47fBlBqxx6eHHcACGLUtXnQjoJd21IeZ4L89ddf5sBDDyzKlCkjn3/+ufTr16/AbjBz5syRDh06mLMsn376qbz11lvmLE337t1LfEYmIyNDduzYIWlpaZ68BSAg+qlpird+hwAgVHkrE8TXxx5kggCwG9pRCMW2VFGPOzwOgugXZvPmzXLw4EGZOHGivPfeezJ79mxzpqUoBgwYYDZ68uTJRX5Nd28mJydH1q5dayJAOsZIdHR0iQZeA0qLfuX27NljgneNGzcmIwRAyPJWEKS0jz181Y0HAAIB7SiEalvKZ91hNNhgDYx66qmnyqJFi2TMmDHy9ttvF+n5nTp1kvHjx0tJ6QGRfoE1QyQ+Pr7E6wNKkwbudKaDzMxMgiDIpfHovXtzb1eqpCFuSgYIsGMPAEGOutagHQW7t6WKPSaIayTGtevKyehMMtWrVxdvCQ8P99q6gNJC1hJOoN36qlTJvc3AqEBAH3sACFLUtXnQjoJd21IeBUFGjhwpffv2NdkXycnJ8uWXX5o+tj/99JN5fMSIEbJt2zb55JNPzP+jR4+WevXqmelzNeKoZ2E0jVUvAAAAHHsAAIDS5FEQZNeuXXLdddeZwUi1r80pp5xiAiC9e/c2j+v92mfXooGP4cOHm8BIXFycCYZMnTq1wMHMvGHEpL98st5RF7cu0nIa9Pnhhx+kVatWZiC3u+++W/bt22dGtNcy+PDDD81j7s5qde3a1Qz2VrduXTnzzDNN2Z133nk+eDcn3/aSuOmmm2TQoEHSrVs3r22bpjtpKrTrtmk3qAULFpjbTz31lHz99dcmHUrPDmr/b50+UemAeG+++aaJGOpjOljeZ599VuTX/uCDD+TVV1+Vf/75xwT2hg4d6nxM9+fHH3/cfNZ33nmnvPTSS87HpkyZYsqzqOnasDmdErdkM5YDISkYjj0ABAnq2sJNuds35T5gTKm2oz766CPp0qWLNGnSxDyu40HNnTvX2TYoiJ7c1/pj8eLFJzymY1K988478sADD4g3PPbYY6Z+uuKKK8SbtL3TunXrPFk++r50yvjC2kSTJk2SZ5991pS1PlajRg2ZMWNGkbOFCmsTaRdW/SyXLVtm6uIJEyac8Hwd+Lx9+/amfWeV/86dO+WCCy6Q+fPnS2RkiTuw5OHR2t5///1CH9cdzpXuJN7aUYLR1VdfbXam888/3/y/ZcsWiYmJcbusTvnXtGlT88X1hE4P7O2doiT0i6MD1nmqKO9Dv7z65clPM4t+/vln8wXTH0hd18qVK81j+iXSL+Eff/whSUlJ5kdS06I9oT8QGmDRKRTz08F49Huhn59+eV1pIEZ/DNatWycNGzb06DUROjQYpvvAXXfdRdopUAwcewDBR4+3dCprZsGDP9pR2iatVKmSMwii67DWU1waBHnhhRe80rbVtoqewC3O84rS7tOTxPm/e4W1iXbu3ClDhgwxbSmrDJcsWeJRl5PC2kTaHVVPJOvraWDFnYcfftjMArd8+XLnfdWqVZPTTz/d9Ca5/vrrxZsYUMOH9MyUTt1j0W5EVaw+//lotsA111yT577//e9/JiNEd6r777/f7KxK79MdpVevXnLOOeeYwINGLTUyqheNvumZMKURUd152rVrJ23btpVp06Y5168RUY0UnnbaaSa7oaCJgvSHRM+4XXLJJWYdPXr0cJ5108fOPfdcGThwoBmsTr9Yun0axbXO4F100UXmdXTbNILqGu3VH7eePXuazJGSlLP+0MXGxpr/9cehTZs2zh9MPXNojQ6sX2aNMnpC19W8eXO3DVj9cdXHC/pBuvzyy00mCexJ90393tx7773OboIAAIQ6PU7VY0OgtNtRejJWG/x68slq+2h75dJLL82zf2qGubaRtI2lbRjXQMPtt99uju81U8PKStAggQZCdJ2uy7tmW+hzbrjhBnMCVZexGvT6mD5Pt0kb+t9++61p1L/xxhvm8ZSUFLnxxhudbbknn3zSud787b7iKqxNtGPHDtOWqVixonN5fcyTIEhhbSL9HLW9WVAQS9ukOuurZn26C4a9++674m0EQXxI05y6d+9udlrdeQvKQNBRbTVip18KV6tWrTLRMv0CzZw500TWLJoRoenAv/zyiwks/Pnnn+ai9+tZZx01X+mX5bfffjOv/d1335muKvp6muZ05ZVXyuuvv24CF7qdrunE+c2bN0+ee+45s/7+/fubHwLXxx599FHzI5H/PeiXvVmzZiY16tdff5Wnn37avJ5FX1PvL0r3FOuHx7pYletVV11lvjgNGjQw92nQ4ciRI873r9kh+sOp71d/bA4cOOBcp/4guq7T9aI/FiWlqXj6GcGe9LPX75uVJlgoHeTxnntyLx4M+AgAQCDR+k6zZ/W4c/v27c5uzQFxPERdG/LtKG3raADitddec3a/cGV1V9f21cKFC027yZVmk2tAQh/XE8v62kq7klhZ6e66y6gVK1aYE7vaJtOMEW3Auz6mJ0f1NS+77LI8z9P2kZ7A1mV+//13890pqN1X1PaH1Z45++yzT9omatOmjSm/OnXqmJPX2m1Iu5RarICSu4tub0loxtg999wj48aNc/u4BpT0s9flvIkgiA/dd9995oulX8b9+/ebMTK++uqrE5bbu3evmf4v/1S/+iWKiooy91977bUmM8SikTJ9TOn9gwcPNtE1jb7dfPPNzmU3bNhgBrPVqOKFF15oXmvTpk2yZs0as16NLir9Ump0sCDaz07TzNQtt9xigjJW5og+ptkq7uh23HHHHea2Rm8vvvjiPF9gjZYWNcpo/fBYF+vMuqZKaZBFAymacaJ93fTLrz8m+h41uqhRYL1P+7tpf3L9PJT2SXNdp+tFfyRKSrdt69atJV4PgpNrha0VW6GyskQ0eKkXvQ0AQJDJycmRBx980NzWMQB0XIHVq1dLx44dTYNSs5f9irrWNu2ogmgbRts9CQkJJss7fza6tnesTA8NDOQPkhRGs0tc21YaSLACgZopoW2mgtpLeoJZt0e3S0/qFtTuKwoNClntGWs9hbWJwsPDzfAC+jzN8NcxODQL5r///jPPtQJK7i6aTVMSmomjbcWaNWu6fVzft7YBNVvFmwJnMIkQVbVqVZOpoBftY6UN9fwD4OhOmb/vlDuuwQLXfl4ajMgfSLD+10if9v/SAIjSPmD6WgV1fSmOk/X3LGjbivLcotIBUfULrReNVmq5//33385ULu0OpBetgFu0aGHS0jQgo5kg1hc8P40UlzQQomWtUVfYk44LoxlMn376qaxfv96crSiwEtP7R448fhsAgCDMAtEz6Zpy/8QTT5j7NBCiXQx0gHlteOU/C16qqGtt244qrN3kyupeb7UvdN8tCeu1CmvzFNaWO9lzPd2WgtpESrP39XLrrbeaYIgOKDts2DDTtpozZ44U1BWpJIEQ7VGggRkdI0U/R81O0QCMNb6jr9pTZIL4kPb3slLh9QukZ4LdDZCpGRg6YIw2klxpw0mfp107dGwPK50pPx2vQ/u6aeaDLq+D0ljL6o6kY28oHVTGSnvSHVzXa+3QmhFx6NChAt+LRgT//fdfZ1+7s846q0gZHLod1jgge/bsMWWizy2Ibpdr+lVRaEqaa5RWzzhouWsAQ2+7noHXLi66Hdp1pjQyQbTCt8Yngf1Y303dr7WCOHz4cMELR0eLPPts7kVvAwAQZKwZ8bQRZWUYa0DEmllPZ/fwK+paW7SjdJ8rqF2jYxFqV5O0tDSTuaTtraLQdepzCguK6IlV17aVZjfothXleFHHvbAGFNY2W0HtPqVdhFyHFyiKwtpE27ZtM209i7YXtTeBVd6+zATRbdLucnr58ssvTVa/awBEx5fUng4aTPWmkMsEKepUtqVBo90PPfSQ6aai6X86IIzrQDeudPDEH3/80dl1RGkWg34BdMfUTA7XAX1cafcUDQJYg9toGpZG7JSODaJ9u/RLaPX1UrpNX3zxhRn4RyNr+hzrMXd0MFSN6Os4JfpjU9RBHvVLo+ldmm6lPzTar07LwR39Iuo0WJqtUtiYIK60X50+RytXfVzfi0ZtNWhUuXJl8wXXQSl11GN9TH9cnn/++RPWUxj9IdLPUX8Qvv/+e/N8zRLRKKpGT7WrkjZudd365dXuONYI1Np/Tz9b2JemGGowDACAUKbHSdOnT3d2d3al6fzPPPOMeVyP9/QYDQGoiFPZBno7SttG2p1Gx7bQMQ1d6TG6dvvQk5TasO7UqVOe8QILou0THXxVG+naZcXduCDavtC2gGZPaLtA2yNFoWMramaGrltptlRB7T4tCx2vxHXQ2Pw0M951Qgft6qIBnILaRJs2bTKZGBr40MwaDfRoNyGdnraoCmsTaTtV25K6DZrVods+cuRI0w49GW1LaVvWk0FaiyLM4c1+ET6ihakNb43oWSPaagHqB1W/fv08aUvBSnc+3dl1cJmizsdcWvRMtg4g5G5OZ2//2GmQ5ZFHHpFQoH0UNVKr001pX0VXobb/4kRagWiQTseF0T6iJ/3x1p/itLTc29qv1cs/9kCg1N/BIFi3G/AnHZhex6jThpy7cbB0XBBtOGrD6LbbbvPLNlLXhuZxqKftqOTkZElMTDQnaHXMEQ2GaJCuJDQIoLN1FjRoqrfo1LX6HdLMfDvo1q2byZLRrOqi7MNFrb8Dq7VtY9rPTQeSsgbPsSPtjxYqARClUU8d6Th/AAT2oJF2HdVc+1sWKdasARDt86kXKxgCAECQ0AxZPZbVmR7csc5sa/9/v6GuDUmetqN04FHdX/UYTRvTOpNLsNDMf7sEQHbt2mUCpvkDIN5AJgjgB6EWgceJdHo0HV1c+4JqlzAdG0dTHa1R80+gU39ZA1+lpIgkJFCsCHrBmlERrNsNBDIdP0AbntqVQLsfeDLbhddQ1xochyLYkQkCAAEavbZGNk9PTzdj+xQ6XbJ2gdHgh16KOM0bAADBQseH0y7WOk6WXwIgiroWQCgOjAoAgRYE0XFBXO9zS8cAIfsDABCEdMBTHQtLBz/UTA93dKwGHWzRr6hrARAEAQDf0FmLVKVKlZwzHhVl9HEAAILKlLvlhUe/lV+Wb5Fxt/WUIX1bB8XsIwDsi4FRAcAHrPnpy5cvby5Kp3EuUEaGyMMP5170NgAAQSA7O0d+/3enud2lWfVCl9XuoS+99JJcfvnlkpmZKaWOuhYAQRAA8G0QRAdXLFIQRA8GdS57vfjjwBAAgGL4e/M+STmSKYlxUdKyTm7mY0F0LJDnnntOvvnmG1m+fHnplzd1LQCCIIFPB5CyphQrzHfffSd//PGH83+do/qaa64Rf3riiScko5hntHWKq549e3p9m4DS0q9fP3nsscekV69eRQuCREaK3H137kVvAwAQBBauzs0COa1JNYmIKDzJXMcF6dy5s7k9f/58KXXUtbZCOwoFCbkj7VSd+qoAEREReaYjLWxZ/ZGOi4uTYKFBEJ2O87TTTjP/6+3PPvvMr9v05JNPyvDhwyU6Otqj5+nAWjVq1JCZM2d6/Jr63EgakAgA5557rrmonTt3SvPmzc34IAWKiREZPbr0NhAAAC/4bU1uEKRz09xBwE+mU6dOMm3aNHPCrtRR1xaKdhTtqEibtKNCbkyQMmXKFHi55JJL8ixbpUqVApft27dvsV7/2muvNQEInQbsvPPOk927d5v7Z82aJW3btpXbb79d2rRpIy1btnT++GvD/ZxzzjHP0/s1gyMtLe2Edffv31+++OIL5/8///yznH766aYimTx5sjz//PPmNd577z3zero+y9SpU6Vjx47mtXWZ33//3W3mhvbR1DPYrVq1kvPPP985kGNKSorceOON5n69aIDD8swzz5gGnq5XL5s2bZIhQ4aYx7p06WLu03JITk6Wm2++2QRqtHx0Gas/6JlnnikPP/ywOWuuZbFx48Y8DcaffvpJ2rdvb56nI4+vWrUqT7nedddd5szCt99+W6zPDfAlnR1G99k5c+ZQ0ACAkLJkXe6xbsfGVYu0vB7PmectWeLT7YLnaEfRjrINRxA4dOiQQzdVry1HjhxxrFq1yly70uUKuvTr1y/PsvHx8QUu26NHj2Jt6549e5y3R40a5bjjjjvM7ZkzZzoiIyMdixYtMv+PGzfO0adPH3M7JyfHsXfvXuftIUOGOF588UXz/4cffui45JJLzO3p06c7zjjjDOf6zzvvPMcnn3xibg8aNMjx+uuvOx/T1+vQoYO5vWbNGkfVqlXNtcrIyHAcPHjwhG1//PHHHdWqVXPs3LnT/H/bbbeZi3rggQcc11xzjSM7O9uRkpLiaNu2rePrr7927N+/31GuXDlHWlqaWS41NdX5mWg5JicnO9d/8803O7dX3+fgwYMdr7zyivlfy1s/H902tWHDBkfFihXN7V27dpnbK1asMP+PHz/e0bJlS+f7DAsLc8ydO9cRTArafxE6dH/Vz9j6bgB25K7+DgbBut1AadPjmIjwMPN92fLBDQ7H5LsKvzgcjh07dpjlw8PDzTElSh/tKNpRwdyOKmwfLmr9HXL5LpqxUFh3GFdWlkZB3WGKQ7ugfPrpp2b06yNHjpgzwJamTZs6szM0a0FHx1YaL3j11VdNtoZmheiAit27dz9h3b1795Z77rnHDCRVtmxZk0kyYcKEk27TjBkzTHZHkyZNnINS6WCN7mj2StWquZH8W265xWSGqP/9738yZswYUy46//vAgQPNfRdffLE0btzYZMD06dPHZKvUqlWrwC47v/32m7z88svmfy0f164y1113ndm2/DRrRbM9WrfOnXJNM2XuuOMO2bFjh/lf31fXrl1PWg5Aabrwwgtl/fr1Mm/ePDnjjDNO/gTtnlemTO5t/R1LSPD5NgIAUBJ6HLd09FWyYuNeqVnxWB12EnpsXL16dXMcp8e0mjVcaqhrC0U76kS0o0JTyAVBtIHui2WLQhs7b7zxhixYsEAqV65suqg89dRTzsddxyPRgIwGPNTnn38us2fPNqnyiYmJ8tprrxWYNq/dPsaOHWuCGNo9JUb7NvpQWFiYM1Bj3XZ9TN+HBjb0PWvXFO3nqV12unXrdsK6dB0aCGnQoEGBKXjuuHtt120r6HlAoMwOowYMGCDLli2Tr7/+2jkoHAAAwUxPjrWuV8lcPNGhQwfTrVu7P5dqEASFoh3lXbSjAlfIjQniTzp+hmZoJCUlmVlR3n777SI/r2LFiiYAouNm6EjGBdFsiR9//FE+/vhj57gbSl/XanTlp2Ns6HP+/fdf87+Ow1HQspqNYmXIvP/++3L22Wc7s1DeffddE5DQQZPGjx9vHtPt3bVrlwl6PProoyYjY+nSpeY5+n5cX0fHGNFxS6zgj77v//7776Tlow1GbTz+888/5v8vv/zSZJu4ZtkAgUS/J4cPH84TBNHBUbdu3Sr79u1z/6T4eE1Py73obQAAQpSOX6dZB1dffXXpvjB1bcCiHUU7qjQRBPEiHUy1UaNG0qxZMxN40C4cRaFdS7QiaNGihele4i6LwhIfH2/S7HWZ2rVr5wmOaEaJNTCqK90mDWhcddVVZmBRHZh0zZo1btevA5MOHjzYDH6qA5zqoKdKAxwazdQuKToYqwY0dOpeDXLoNuv9um4NsAwaNMg857777pOzzjrLOTDq6NGjzYjD+r8uq0EUPQNwMppVo12MtBuMDuw6btw4czYdCFRHjx51DvprBUE0KFjoyOua2VS5cu7FTeYTAACB5sUXX5Qxk5fJjv0Fz7jojna99nT2QK+grg1YtKNoR5WmMB0YRAKcnlHVhoQ2uDXjwWpkbNiwQerXr5+nm0moy87ONqNqa7ebwoIlxaGzw2gwxhqrBL5j1/3XLjTrQ/s7a+BQM580XVgDh1OmTDFBSg00Anbgrv4OBsG63UBp0iaEzuS3f/9+My5I2waVT/6kAWNKY9NwEnY9DqUdFfr78OEi1t9kggQRHWNEx9PQvpPeDoAA8B6rG5j++FqDLFtj12gXMrcyMkSefTb3orcBAAhge/bsMQEQTa5oWrOCxwEUPSHQsWNH50D3pYK61rZoRyGkB0YNZXomWS++opkgALw/KKprEKTAkde1+8wjj+TevuceHXKfjwIAELBWrVplrutXLStxMZ41KTRTcu7cubJ27VqzHs2eLBXUtbZFOwquCIIAgJfpoL2PPfaYGcOnyEGQyEiRm246fhsAgABmDVjfonZSsZ7fsmVLZxBEx6QrFdS1AAiCAID31alTR5588sk89+mMRs2bNzezR7ml012/+y4fBwAgqDJBWtSuWKzn64QA3333naxcuVJKDXUtAE/HBNFZOXRWD+3nrhedulSnXi3M7NmzzVzgOmCJjmfx1ltvUfAAbGfYsGHmgPGBBx7w96YAQYVjDyDQgyBJxQ6CuK4HAAIyCKJnMp9//nlZvHixuej0pxdccEGBEVwdsbVfv35mEM+lS5fKyJEj5a677pKJEyd6a/sBICAHi9M0Yb0GUDIcewAB3h2mTvG7wyhtRwTBZJUA7BoEGTBggAlqNGnSxFyeffZZ08/9t99+c7u8Zn1oWvjo0aNNGvhNN90kN954Y0hPwaqDi2YwswNgax9++KE5wzV8+PCiPyk1VSQhIfeitwEYHHsAgWn16tXy+++/S+u6xesO07RpUzNAqs4ws3v3bikV1LUBjXYUSkt4SeZZ/vLLLyU1NdV0i3Fn4cKF0qdPnzz3nXPOOSaLJFNHZy5Aenq6mePX9RIsdByAgoIgWVlZpb49AEqf/i6qBA1oHDNr1ixp1aqVXHrppQU/MS0t9wKgVI89gvm4A/AX7Rp/2mmnSWx08QbzjouLk8aNG0ujRo1KLwiiqGsDFu0oBGwQ5K+//jLZHzExMTJkyBD59ttvnX368tu5c6dUrVo1z336vwYD9u7dW+BrjBo1ykwtaV1q165d9A3UxodeXNPqNCih96Wnu182J+f4fXqApPcdPSqe0vJQXbp0kbZt25of9Ouvv950ATr33HOlTZs25nGNervOEFGpUiXZuHGjua2jZPfv39/Mm67Lv/nmmwVGSq+66io577zzTOVx+eWXmy5H2kVJx17R8QdcPwd9XCsqHdNFZ62w3H///ea1dHt79OhhXl/p9uh26bI6pou+xrRp0zwuE8CO3AVBtPGlKb/Wd+wEcXHahzD3orcBlNqxR4mOOwAUm1Uvtm7dunRKkbq2cLSjaEfZRHhxUteWLVtmusDcdtttMmjQoEIHNNIGvyurz1/++12NGDFCDh065Lxs2bKl6Buo01DqxfVA58UXc+8bOjTvslWq5N6/efPx+8aOzb1v8GDxlDXo64IFC0wZVdH1i8i8efNkwoQJJx39Ws9wXX311fLyyy/LokWLzNksXeeSJUvcLq9ntT777DNZs2aNuTz00ENmoFo9WBw/frz8+++/Zjn9jIYOHSp//PGHWZde6wGkevDBB81r6fbq53nvvfc6179v3z4TAPnzzz/ljTfeyPMYAM+CICedIjc8XKRevdyL3gZQasceJTruAGzo7bffljvvvFPmz59fovVElvaU8NS1haMdRTvKJjz+5YmOjjZZAerUU081DegxY8aYH8P8qlWrZs7IuNLsCP3Bq1ix4P6DeqZHL6FCszCsBlBhNJChgZIrr7zSeV9ycrI50Gvfvv0Jy2t6r56xUprhoZkjVtnpAeP69eulZs2a8uuvv8quXbucz9NGmPbjVNOnT5fXX3/dvE5OTk6eFGBtwOnAt0rTjtetW1fCkgDsoVhBEAB+O/YIteMOwNcmT55sMoS1m+cZNShv+A7tKPhCicOvenZF+9K6ow3nKVOm5LlPG916ABMVFSU+YTUw4uOP33f//SL33KPh5rzLWv0PXVPP77hD5OabRSIivLZJ+QMgERERJuvDcvRY1xstS+2Come7ikKnHXZdZ/7/NfVXAxt65ksPGPOX+ebNm01XHc0M0S40K1asMN1pClq/6zYD8HIQRLviaSaa9Tvkq99IIAQE3LEHYDPWiTEd00NSiz/FrZ6U00kTrExqn6OuLRztKIN2VOjzKOdap7idO3euGS9Cu1w8/PDDZrC/a665xplOOnDgQOfy2m9306ZNZnwKnUbrgw8+kPfff9+zGRM8Zc2u4JryGh2de1/+szzWsq6p53qApPe5BAA8kZiYaFJpC9OwYUMzmraaNGmSs8Gk2Rvx8fHyySefOJf977//zKjZxaXbo1MU69TGlu3bt8vWrVvNdurZNT1rpgeU2uUFgG+DIGlpae4Dijp2kXY50wszTAHBdewB2IjWYZptbB3TloTWjdr9W0/WlcoEAtS1haMdlQftqNDlUSaIdqm47rrrZMeOHaYbhnbB+Omnn6R3797mcb1fswss9evXN6lyOpbE2LFjpUaNGvLaa6/JJZdcIqHqvvvuM9kUOuK1nnlyR6cMvuOOO8yYIT179nSm52qqrp690vLSaYS1kqlcubIZ96Mk9Pl6MGgNOqUVjo41ot1nLrvsMjNPu05lbH2OAEpGf+OaNWtm0oTdZYRpkERH1c9Ds8+uvvr4bQAGxx5AYNETaTrYt55Iq1WrlsiK4q9L2waaeaxZ0dqG0Mxkn6KuDWi0o1BawhzWaGEBTMep0KCLZi5YDQf9sdywYYMJtLh22wCCAfuv/ehPbb169Uy215w5c0yAEwh17urvYBCs2w2Uhl9++UXOPvtsk8FsxpibcnfRnzxgzAl36ck4Hf/u559/PmF6a/gGx6EI1X24qPU3UxAAQCnQsXk0RV/T8wmAAACCfTwQa7DikrLWo13AAaA0EAQBAB8cIG7btq10+jcDAFCKrK7vJR0PxEIQBEBpIwgCAF6ms1BoP2mPzmrpYKraRUYvxwZWBQAg0Dz99NOyd+9eM0ixN4MgVoaJT1HXAgiFIIhOAQsEmyAYigdenh1G6cDSOiBxgdMA7t2bewEAIIC7d+qg/jrAvzc0adLEjJnlrfWdFHWtE+0o2LUt5dHsMIFER6QODw83071q/3r9X3+UgWD40u7Zs8fsr1E6JTNCio6Yrxd3QZC1a9fKihUrzBm0E8TFifz99/HbAADYQK9evcwAh6WCutagHQW7t6WCNgiiARAdDVan5dVACBBM9Eur3SUimAo1ZLNA3AVBdOpsdeTIkROfGB6uQ+T7fgMBACgmbXjccMMNJnvj5Zdf9vwEpCczyRQwm0yJUNceKwbaUbB3WypogyBWFLNOnTpm8MHs7Gx/bw5QZBq1JAAS2kEQ/Xz1N8qVTo+r0tLS/LJtAACUxL///itTp06Vv//+W1555RWfnOEls7t00I6CndtSQR0EUVYaDN0KAARSEKRMmTInHMgVGgTRLjQffZR7+/rr9de9FLYWAICiswb89tb0uJb7P5wn42etlqeu7iQ3n9PKdx8JdW0etKNgV0EfBAGAQJKSkuK2K8xJgyAZGSK33JJ7++qrCYIAAAKONYOLt6bHtaRnZsvOA2ny345D4lPUtQAIggCAdyUlJcmdd97pNghS6JggmtJ3wQXHbwMAYJMgSP2qZc31xt2HxaeoawEQBAEA79Jp/l577TW3j1WqVElq1KghsbGxJz6o9333HR8HACBgbdy40Vzr5ATeVK9KKQVBqGsBEAQBgNLzzDPPmAsAAMEcBNGAvzfVq5KYu/5dPg6CAIDOkEQpAIB3B0bduXOnJCcnU6wAgJCRmZnpHPfK60GQY91hdh86ImnpmV5dNwDkRxAEALzom2++kerVq8vll1/u2RN1sFQ9qNQLU+gCAAKMzsR48OBB2bdvn+ne6U3lE2KkbHzutPKbdvvwJAJ1LQCCIADgXUePHs0zCKqrH374Qc444wy57777TnyiwyGyaVPuRW8DABCAU6rqAOD5p4D3xnq7tqghPVvXkqzsHPEZ6loAjAkCAN5lzfziLgiyf/9+WbBggSQm5vZ9PmGwtj/+OH4bAAAbmfrY+b5/EepaAGSCAEDpBUHi4+PNdZq77i46bV/HjrkXpsgFAASYsWPHygUXXGC6fQYt6loABEEAoPSCINZ9boMgAAAEsIULF8rkyZNlw4YNPn2dbF92hwEAgiAAUPqZINYyeWRliXz2We5FbwMAYIPpcS2z/toqNa5/X7o+NEF8hroWAGOCAECAdIdJTxe59trc2xdeKBIZyUcDALBNEKRcfLTs2J8qOTk+HBycuhYAQRAA8K5OnTqZQEj79u1PeMwKjLjNBAkPFzn77OO3AQAIEBkZGbJ9+3afBkHqVS1rrncdTJMj6VkSF+ODkwHUtQAIggCAd1133XXm4k5CQoKZGaZMmTInPqgBkhkz+DgAAAFny5Yt4nA4TDC/cuXKPnmN8gkxUjY+Wg6nZcjG3Yelee0k778IdS0AxgQBgNLTsGFDOXz4sPz3338UOwAgKLvChIWF+eQ1dL31quRmg2gQBAB8hZxrAPAiDXKkpKRIdnY25QoACAkHDhww2Yy+6gpjqVcl0Vxv3EUQBIDvEAQBAC8677zzTJeXSZMmefZEHSy1ZcvcC1PoAgACyKWXXirJyckyceJEn76ONS7Ixt3JvnkB6loAjAkCAKU3O4zq16+fyRTRIEmlSpWOP+BwiKxadfw2AAABRLurFFS3eUubepXkzNY1pe6xjBCvo64FQBAEAEo3CDJnzhxJTU013WbyBEFiY0Vmzjx+GwAAm7mxd0tz8RnqWgAEQQCgdIMger8GQdLyd3mJiBA580w+DgBAwBkwYIC5HjNmjDRo0ECCFnUtAE/HBBk1apR07NjR9HevUqWKXHjhhbJmzZpCnzNr1iyTPpf/snr1aj4AACHn6NGj5jq2gGyO+Pj4PMESAIXj2APwr5ycHJk+fbr88MMPEqFBhFKQkZltpuQFAL8HQWbPni133HGH/PbbbzJjxgzJysqSPn36mLOaJ6PBkh07djgvjRs3Lsl2A0BQZoJYQZATMkGyskS++y73orcBGBx7AP61Z88eycjIMCcxa9So4dPX0sBH3cEfSuylY2Xr3hTvvwB1LQBPu8P89NNPef7/8MMPTUbIn3/+Kd27dy/0ubpc+fLlKXQAYvfuMK7LOaWni1x0Ue7tlBSRSI9+noGQxbEH4F9btmwx19WrV5eoqCifvpYGWsLDwsz4pVv2pkhtb78AdS2Ako4JcujQIXOdlJR00mXbtWtn0sRbtGghjzzyiPTs2bPAZdPT083FogMIAkAwuPzyy02WR7ly5TzLBAkPF+nS5fhtAKV27MFxB3DyIEjt2l4PSbhVu1IZ2bj7sGzZ64NpcqlrAZQkCKLpasOGDZOuXbtKq1atClxOo8bvvPOOdOjQwRxkfPrpp9KrVy8zVkhB2SPa//fJJ5/kAwIQdD7++ONCH09ISDCBEO1OmIdmiMyf79uNA4Kcr449OO4ACrZ58+ZSDYLUqZw7Pe7mPT4IglDXAihJEGTo0KGyYsUKmTdvXqHLNW3a1FwsnTt3NhHll156qcAgyIgRI8xBjmsmSGn98AKAr1P7Nd0XQOAce3DcAQRWJoh5XV+MCQIAng6Marnzzjtl8uTJMnPmTKlVq5bHz+/UqZOsXbu2wMdjYmKkbNmyeS4AEAwj6GvqfWEj2hMAAYrHl8ceHHcABcvOzjZZjCGRCQIAngZB9MBez8JMmjRJfv31V6lfv36xCnHp0qUmVRUAQsn69evNwKcFjQdSKB0otWPH3AvT5wJOHHsA/vXqq69KcnKy3HXXXaXyerWPBUF8MiYIdS0AT7vD6PS4n3/+uXz//feSmJgoO3fuNPfrAb8144GmlG7btk0++eQT8//o0aOlXr160rJlSzO91vjx42XixInmAgChxJrxJTY2tsBl3nvvPRNI1gFUr7/++uMP5OSILF58/DYAg2MPwP80izEiIqJUXqtR9XJyZuua0rpuJe+vnLoWgKdBkHHjxpnrM88884Spcq2D+R07djgHUFIa+Bg+fLgJjGigRIMhU6dOlX79+vEBAAgp2hWmsOlx1erVq+XHH388cVDHmBiRH344fhuAwbEHYC/NaiXJzGcv8c3KqWsBeBoEKayfu+Wjjz7K8/8DDzxgLgBgl0yQwoIg1mPWsk6RkSL9+/t2A4EgxLEH4D96ErN///7SoEEDk8Ud9ONaUdcCKMnsMAAALwZBAAAIMJs2bZLly5fLwYMHSz0AkpGZLTlHjxbaxRQASm12GABA8cYEsYIgaWlpeR/IzhaZMSP3orcBALDZ9LiWK1/8UWIvHStffPGFd1dMXQuATBAACJBMEB1PpE+f3NspKSIJCXw0AAD/mHK3udry8xJzXTt8t/O+0pAYFy3aC98KwngNdS0AgiAA4D01atSQiy666MRBT4sSBAkPF2nT5vhtAAD8zJqmtnal3GlrS0vtSmXMtetkC15BXQuAIAgAeE+PHj3MpTBWECQrKyv/AyLLlvFxAAACxpY9KXmCEqWlTuXcoIvXM0GoawEQBAGA0nXJJZeYAEhERARFDwAIkkyQ0g2CWJknXs8EAQCCIADg3ak8TzZ6PsEPAECwiI+JlITYKKl9LDPDH5kgRalbAcATdDwHAC955JFHJDIyUu6//37Pn6xjhJx5Zu6F6XMBAAFg9qhLJfmrIdKuQeVSfd1aFXMzT1JTU+XAgQPeWzF1LQAyQQDAe3Sw0+zs7EKzPTZt2iQPPPCAGRvko48+Ov5ATo7I7NnHbwMAEAD8kYURFxMpF5zeQMo17SqZmZneWzF1LQCCIABQulPkpqWlyddffy0VKlTI+0BMjMjXXx+/DQCAjX338HkiA8Z4d6XUtQAIggCA94MgsbGxnk+RGxkpctllfBwAgIAwacF/8tRXf8h5HevLM9d2lpBAXQuAMUEAoHQzQeLj48310aNHzWBvAAAEojXbDsjyDXtly57cGWL8IT09XQ4ePOi31wcQmhgYFQBKMQji+pgGQpyys0Xmz8+96G0AAPxoy94Uc13aM8NY3vhhucmsHDp0qPdWSl0LgO4wAOA9VlCjqEEQDZo4/9fndu2aezslRSQhgY8GAOA3W/bmZoDUrpQ7U0tpq1Q2zjlNrtdQ1wIgCAIA3tOhQwfTxaVWrVoFLqNT6OolKysr77ggOvp+o0bHbwMAEAiZIJX8kwlS61jwxatBEOpaAARBAMB7Ro0aVaTlNPsjOTk5bxBExwpZu5aPAwAQYEEQ/2SCWK+7detWycnJkfBwL/Tip64FQBAEAErftm3bTD/nqKgoih8AEHDS0jNlf3JuF886fhoTpEZSgoSFhUlmZqbs3r1bqlWr5pftABB6GBgVAEpZYmIiARAAQMA6kJIuzWtXMIGIcgkxftmGqMgIqV69uve7xACwPYIgAOAl9evXl/Lly8vKlSuLN1hb//65F9dZYwAAKGU1K5aRVWOvk20fDfZr2deuXdvZJcYrqGsB0B0GALznwIEDcujQoZNmeTz99NPyzz//yPDhw6V9+/bHp+2bNu34bQAAbK5fv37StGlT73WFoa4FQBAEALzHGuhUx/sozE8//SQLFiyQSy+99HgQJDpa5MMPj98GAMDmHnvsMe+ukLoWAEEQAPCO7OxsycjIcM7+Uhjr8Tyzw2j2yPXX83EAAPzuvvfnyi8rtsiDF3eQq3o0lZBBXQuAIAgAeEd6errzdv4gyIhJf+X5f+PBTHP92fy18ndc3sdGXdyajwQA4FcrN++T5Rv2ytHM7ICoX/fv3+8cJBUASoqBUQHAC1yzOk6WCRIVndtdJjPjeOAkLDtbqm9YLbJsGWOCAAD8asveFHNdu1IZv27H4sWLTRfT008/3Xtjgmg9S10L2FqkvzcAAEIpCKKDokZERBS6bGRM7nSDWRnHZ4GJzEyXu4ZfnvtPSopIQoIvNxcAgAJt2ZtsrmtXSvRrKVnZH9u3bzfdTk9WvxZpdph27XJvU9cCtkUQBAC8IDw8XHr27ClhYWEn/+E9lgmS5ZIJImFhciipipSLjTS3AQDwB53lLPlIbrfN2pX9mwmis8JERkZKVlaW7NixQ2rVqlWyFWr9WqPG8dsAbInuMADgBTVq1JBff/1Vfvnll5Mua3WHyXIZRyQzJk6ef/d/Itu2icTH85kAAPxiy5Yt5jopMVbiYwqf8t3XNPND61fX7SoRrV+1nqWuBWyNTBAAKGVdr7hNulx6i0THEuwAAAQWK9jg7/FALLVr15bNmzfL1q1b/b0pAOyYCTJq1Cjp2LGjJCYmSpUqVeTCCy+UNWvWnPR5s2fPlg4dOpiBjRo0aCBvvfVWSbYZAIJaTHwZiUssJxE6VR+AQnHsAZQuh8MhzWtXkGa1KgRMEMRrmSAA4GkQRIMZd9xxh/z2228yY8YM0z+vT58+kpqaWuBzNmzYIP369ZNu3brJ0qVLZeTIkXLXXXfJxIkT+QAAhAz9TaxUqZL07du3WM+PzEiXq1+6T+Syy3IHbgNgcOwBlC49bl819jr58v7i1WcBHQTR+lXrWepawNY86g7z008/5fn/ww8/NBkhf/75p3Tv3t3tczTro06dOjJ69Gjzf/Pmzc10Vy+99JJccsklBc4HrhfL4cOHPdlMACh1ycnJsm/fviL9Xm35Z4msnP2DVKrTSE7td7W5LywnW1ovnJG7wEcf+XpzgaBRGsceHHcAgatz584ycOBA70yTq1PkTpiQe5u6FrCt8JKOHq2SkpIKXGbhwoUmW8TVOeecYw5GMjNzR552l/parlw558WKAANAoDp6LHsjLi7upMse2LFZls2YIBuWznfelx0ZJd/fNFLkjTdEoqN9uq1AMPPFsQfHHUDguuiii+Tjjz+WK6+8suQr0/pV61nqWsDWwkvSX3DYsGHStWtXadWqVYHL7dy5U6pWrZrnPv1fu9Ls3bvX7XNGjBhhDnKsC30AAQS6I0eOFDkIEhkdY64zM453e8mJjJLf+l4pcscdIowVApTqsQfHHcBxnTp1krZ3fy5/b9oXesWi9avWs9S1gK0Ve3aYoUOHyooVK2TevHknXTYs3zzcehDj7n5LTEyMuQBAaAZBjk2Rm3G82x8A/x17cNwBHP+eLF++3GQ3xscEziSS2mVt27ZtUrduXTNtLgCUeibInXfeKZMnT5aZM2dKrVq1Cl22WrVq5oyMq927d0tkZKRUrFixOC8PAEEdBImyMkHSj2eChOXkSMXtm0TWrhXJyfHhlgLBiWMPwPd0bCure2fNigkBE5jR7vENGzYs+TS5Wr9qPUtdC9iaR0EQ/RHSszCTJk2SX3/9VerXr1+kwYx01gRX06dPl1NPPVWiSPkGEGJBEJ0KvKjdYVwzQSIzjsrwOweINGmiK/PhlgLBhWMPoPRs3rzZXFctHy8xUYGRCaLZWzVq1DC3S9xFXutXrWepawFb8ygIotPjjh8/Xj7//HNJTEw0GR56sQ7+rX61OoKzZciQIbJp0ybTh/eff/6RDz74QN5//30ZPny4d98JAPiRZr1pcLdBgwYedIfJOxXukfhEkXLlfLaNQDDi2AMoPVaQoXalMgFV7F6dJlfrWepawNY8CoKMGzfODFR65plnSvXq1Z2Xr776yrnMjh07nFFkpdki06ZNk1mzZknbtm3l6aefltdee63A6XEBIBjdcsstsmjRInnwwQeLlQmSGRsvT306X+TgQZGEwEhBBgIBxx6AP4IgiQEZBClxdxitX7Wepa4FbM2jPDdrULHCfORmzu0ePXrIkiVLPNsyAAhRSdXrym1v/SxRMSfvOgPYHccegB+CIJUDJBNkyt3mqvbR1eZ6y5zPRJoVEAgZMKY0twxAEAuMzn4AYCMRUVFSrnJ1f28GAAB5VKhQQZo3by5NapQPqJKpXTk3M2XL3mR/bwoAu84OAwA4sTtMvXr15LPPPitW0URkZsilrz8icv31OhcgxQsAKHUPPfSQrFq1Su7o3yagSr9WxdzMlC17U0q2Iq1ftZ6lrgVsjUwQAPCC7du3m0GgrakFC+PIyZFfP3nFDIza87phEh0XL+HZWdJh1uTcBcaOFYnJHTcEAAC7a147SQb2bCYt61Qs2YqyskQ+/jj3NnUtYFsEQQDAC6xZsuLi4k6+cFiYLJ463gRDulx6iwmC5EREybTr7pV+rauLMH04AABOjWuUl4/v7VPyEtH69YUXjt8GYEsEQQDAC6wMkKIEQcLCwswMMZlHjzhniMmOipK5F94g/S5uzecBAPBLRmOrVq3MzI6LHz/D1FUhJzpa5P77/b0VAPyMMUEAoLQzQfQEVHTuzDDaJQYAgECYGebAgQOye/fugAyApGdmyfqdh+RgCuNmASgZgiAA4IcgSKQzCJJ7MBeWkyNl9+0S2bZNJCeHzwQA4JfpcevUqROQJd//qSnS8JaPZcqi9cVfidavWs9S1wK2RncYAPBLECR34NPM9NxMkMiMozLilt65D6akiCQk8LkAAEo9CFK7du2ALHXnDDF7SjBDjNbVtWrl3qauBWyLIAgAeEHjxo1NAKRs2bJF+/GNyZsJorIjIiUi8DKQAQA2sHnzZpcgSIYEmtqVcoMgW/eVcJrcSJo/gN3xKwAAXjBjxgyPlo86lgliBUEyY+Plka+XyCgGRgUA+D0TZF3AfQa1KyWa6y17k4u/Es2yzMz03kYBCEoEQQDADwbc87xIjkMSKlSk/AEAfhfoQZBaxzJBtuwtYSYIANsjCAIAflC+Sk3KHQAQMJo0aSIpKSlmilzZMksCtTsMQRAAJUUQBABKSA8a27RpY8YE+fPPPyUmJreriyciMjOk/0cvivxSUeSVV0SKsQ4AAIrr008/Pf5PblJIQHaH2Z98VNLSMyU+JsrzlaSniwwblnubuhawLYIgAFAMIyb95bydemifrF+fO2Xf4z+skbCwk49uuub3X2T7vyukfpvOUu+UThKenSWdf/oq98EXXiAIAgCAi3IJ0XJTn5ZStXy8ZGUXcyr5rCyRN9+krgVsjiAIAJRQlp5ZOjbtbVECIGrD0vmybMYEiYqJM0GQnIgo+d/lQ+Ts5lVFoopxdgsAgGJyOBxFrr/8Rbfv3aG9SrYSrV8ff/z4bQC2FO7vDQCAYGfN8KJBkKKylrWemx0VJb9ccbvIE0+IREf7aEsBADjRF198IUlJSTJ48ODQLh6tX7Wepa4FbI1MEAAoocyMo7k/qNGxRf/xPbZs1rHnAgDgz5lhDhw4IBkZGQH9IaRnZsnWvSkSFRkhdSrnjhECAJ4iEwQAAiATRBwOiU09LHLwoLkNAIB/pscNXKO+WSyNbv1Env16UfFWoPWr1rPUtYCtkQkCACVkZXNERhW9G0vUsSCIlUUSlX5EHh/YNffBlBSRhAQ+FwBAqQZB6tSpE9Albs0Qs2VvcvFWkJYmUqFC7m3qWsC2yAQBgJL+kEZESoVqdaRclZrF6A5zLBMEAAA/CZZMkNqVy5jrLXtT/L0pAIIYmSAAUEJ1Wp4qt479wbMf35i8QZDMmDh5+Ks/5dmLWotE8tMMACg9mzdvDoogSK2KZUqWCRIfL2KNe0JdC9gWR9oA4AdNTusp1Ru2lNjEsrl3hIVJTmQUU/YBAEpVWlqa7Nu3LyiCIFZ3mEOpGZKcliGJ8R7OpqbTADM1LmB7BEEAwA/iEsubCwAA/pSamirnnHOO7N27V8qXD+x6SYMe5RKiTRBEs0Fa1Kno700CEIQYEwQASuivWZPl/WGXypwvxhZ7HRGZmdL345dF7r//eKouAAA+VrlyZfnpp59k8eLFEqaZEgHu+OCoxRgXROtXrWepawFbIxMEAEooZf9u2bPpX6nWsEWRn3N4305Z8ct3Eh0bL6edP1DCszOl++SPcx984gmRaA9TfAEAsIEbejWXAynpUrdKbjDEI5mZIi+9lHubuhawLYIgAFBCmenpeaa9LYrUA3tl3ldvStlK1U0QJCciSuacP0i6N6lMf2UAQKnJzs6WiIiIoCnxYRe2L/6TdTyQ4cOP3wZgS3SHAYASyso4aq4jPQiCWMtaz82OipIfB90n8uKLZIEAAErNbbfdJklJSfLmm2+GfqlrlqXWs9S1gK15HASZM2eODBgwQGrUqGH6DX733XeFLj9r1iyzXP7L6tWrS7LdABAwrGluI6Nzp70tCmvZzGNBEADucdwB+H563AMHDkhMTNED+f6UlZ0j63cekmXr9/h7UwDYJQiiI0i3adNG3njjDY+et2bNGtmxY4fz0rhxY09fGgACPAhSnEyQdHE4HCIOh4RnZeb2V9b/ARgcdwC+tWnTJnNdr169oCjqhat3SMNbPpZLn5/m+ZO1ftV6lroWsDWPxwTp27evuXiqSpUqAT/tFgCUrDtM0TNBoo4t68jJkZysLInNzpSnrumU+2BKikhCAh8GwHEH4FMahLeCIHXr1g2K0q5bpay53rw3WbKzcyQiwoNzumlpImXK5N6mrgVsq9TGBGnXrp1Ur15devXqJTNnzix02fT0dDl8+HCeCwAEquj4MpJQvpLEJhR9pHrXrBG6xADex3EHcHJ79uyRI0eOmK7qtWvXDooiq5GUIJER4ZKZlSM7DqT6e3MABCGfB0E08PHOO+/IxIkTZdKkSdK0aVMTCNE+vgUZNWqUlCtXznkJlh9lAPZ07q2Pyp3v/yqnnHVhkZ8TERUtEhbm7BKTGRMnT34yT+TAAZH4eB9uLRDaOO4Aim7jxo3mWsf6C5YxQTQAUqtibjbHpt3Jnj1Z61etZ6lrAVvz+RS5GvTQi6Vz586yZcsWeemll6R79+5unzNixAgZNmyY83/NBCEQAiCU6Fm36577RCIioyWuTDkTEDmaUFaEboNAiXDcAYTueCCWelUTZePuw+ZyRosaRX+innygngVsz+dBEHc6deok48ePL/BxjUQHSzQaAIqrZpM2FB5QCjjuANzT8frOPfdcOeWUU4KqiOpW1nFBtnmeCQIApZEJ4s7SpUtNuioAhILvX3lAkvftkrMHPyjVGrQo1joiMjPlzEnviqyoKjJypEh0tNe3E7ArjjsA93r37m0uwaZeldwxuDQTxCMZGSLPPZd7m7oWsC2PgyApKSny33//Of/fsGGDLFu2TJKSkqROnTqmK8u2bdvkk08+MY+PHj3apNi1bNlSMjIyTAaIjg+iFwAIBTvWrZSDO7c4p8otqhW/fi8pB3ZLy279pVLZ8nL212/lPnD//QRBgGM47gCQ31mn1JbsHIec0dzDk6o6Ne6TT1LXAjbncRBk8eLF0rNnT+f/1tgdgwYNko8++kh27Nghmzdvdj6ugY/hw4ebwEhcXJwJhkydOlX69evnrfcAAH5lBT9cZ3wpikVTPpE9m9dK9UatJKlCJVl47hXSuUFFkUi/JOkBAYnjDsB30tLSJD4IB+Pu3qqmuXhM69fbbz9+G4AthTl0gvAApwOj6iwxhw4dkrJlc+cGBwB/GjHpL+ft0YO6ytGUw3LTmO+kUq0GRV7HJw9dI9vX/iWXPDRGGnfMDS6Puri1T7YX8Idgrb+DdbsBT2gTQPfvqKgo02Wsbt26xx+ccnfwFeaAMf7eAgBBUn8TAgUAL2WCREXHevS8yGPLe9qNBgCAktq/f7/pbqaqVq0adAW6YechMybIaU2qSUJslL83B0AQCff3BgBAsJ9Jc3aH8XBWK6v7TGY6QRAAQOnauHGjua5WrZrExnoWxA8E3UdOlLMe+Vb+3rTP35sCIMgQBAGAEnDN4rAyOzwNgmRlHJWoo2nyzOXtRaKiRFJT+UwAAD61adMmc60TGASjelVyU9037fFghhitX7Wepa4FbI0gCACUQHZWhsSXrSDRsfES5eHAqFExebvDRGRniWRl8XkAAEotEyRYgyB1Kx+bJndXsmdP1HqWuhawNcYEAYASiE0oK3d9OLt4P8BWd5iMo5IVHSuj3pkhI/o1F4mL4zMBAJRKECTPgKhBpF7VYmSCaP26devx2wBsiSAIAPjJaQMGSaseA6Rc1ZriCA+XwxWritQsxpR/AADYrDuMMxNktweZIOHh1LMACIIAgL9UrFVfKkp9PgAAQKk77bTTJD09XVq2bBncY4Ls9iATBADIBAGAktmzea1Mf2+UVKhaW/rd8WSx1xORmSldpo4XWVdd5O67RaKj+WgAAD7z8MMPB3Xp1q1yPBNEZ2oLK8qTMjJExozJvU1dC9gW3WEAoARSD+2XLSsXy9Hkgx4/d/fGf2XzykVSvlptad7yVOn36au5D9x+O0EQAAAKUadyojx8eUepVyVRcnIcElGU0srMFHngAepawOYIggBACWQePZL7Y3psphdPbF61WP73wf9Js859pOkpneTPM8+XDnUriETy0wwA8J2jR49KZmamJCbmZlMEo9joSHnm2s6ePUnr10GDjt8GYEtMkQsAJZCVcdRcR0V7Psp8lMvsMNlR0TLhzmdEPvpIJMazqXYBAPDEzz//LGXLlpWzzz7bXgWn9avWs9S1gK0RAgWAEsg4lgkSFet5ECQyOjd7JCsjnc8AAFDq0+NWqFAhqEt954FU+WfLfqlSPl6Cc3hXAP5AJggAlEBWupUJ4nl3mMhjmSAEQQAApSnYp8e1jJ68TM565Ft5+6e//b0pAIIIQRAAKIHM9JJkglhBkKMSdTRNHrvuDJHy5UVSU/lMAAA+zwSpW7duUJdy/aq50+Su33moaE/Q+lXrWepawNboDgMAJaDT8mkwI6oYA6NGxcTlyQSJS0vmswAA+Nz69evNdf369YO6tBtWK2eu1xU1CKIOebAsgJBEEAQASqDzxYPNRYMhHv8AOwdGTZes6Fh56fUpMvycpiJxnmeVAABQFFpfrVu3ztxu2LBhUBdag2NBkA27DktOTo6Eh58kyV3r13//PX4bgC0RBAEALwgLC/P4OUk16sllI8dKTEIZcYSHy74adUUaN+bzAAD4zJ49eyQlJcXUW8GeCVKncqJERoRLema2bN++XWrVqlX4EzRIQj0L2B5BEADwk9iERGnYoRvlDwAoNdnZ2XLrrbfK4cOHJSbIp2TXAEjdyommO4xmt5w0CAIADIwKACUz7+tx8s1zQ2X90vklWk94VqZ0+vFLkbFjRTIz+VgAAD5RvXp1eeutt+Tzzz8PiRJuUC13cFSri0+htH7Vepa6FrA1MkEAoAS2//uXrF86T5p2Otvj52ZnZsrKuVPN7DAdu50nF7z3XO4D118vEhXF5wIAwEnc3u8UubxrY+nRo8fJyyojQ2ToUOpawOYIggBACWRmHM0z04sncnKyZdrYx8ztNp37yF+de0vrmuVEIiL4TAAAPrFlyxapUKGClClTJiRK+MJOxwZ3Lcogr1q/Xnrp8dsAbOkkQygDAAqTefSIuS7OFLnW7DDqqDjk8+Evi3zzjUis5+sCAKAorrzySklMTJRvv/3WfgWm9avWs9S1gK0RBAGAEshKL34miI7M75wm99h6AADwJWvsjDp16oREQWdkZsvsv7fKp59+6u9NARAkCIIAQAlkpB/LBIn1PAiirCBIVkY6nwMAwKd0atxdu3aZ2w2L0n0kCOj0uGeOnCQDBw6UQ4cO+XtzAAQBgiAAUAJZ6cXvDqMio3OfF55ySB66+WyRmjVF0tL4TAAAXrd+/XpznZSUJOXLlw+JEk6Mj5Yq5eLyvL8Caf2q9Sx1LWBrBEEAoAR0cNPidodxzQTJzkiXcvt3i2zfLuJw8JkAALzO6goTKlkglgbVyhVtmlytX7Wepa4FbI3ZYQCgBO75eJ7kZGdJWHjxRpmPOhYEScvOkdde+lru6tWYgVEBAD4RqkGQhtXKyW9rdp48CKIDoy5devw2AFsiCAIAJRQeUfyf0l7XP2Cm2a3coJnsKJck0rY1nwcAwPum3C3rf51pbjbMWWf+DxUNqxcxE0SnxW3btnQ2CkDodIeZM2eODBgwQGrUqGFmNvjuu+9O+pzZs2dLhw4dJDY2Vho0aCBvvfVWcbcXAEJKvTadpHHHMyVeAyAATsBxB+A9PVvXklvOaSXdWtQIqWJtXD13fJO1a9f6e1MAhGIQJDU1Vdq0aSNvvPFGkZbfsGGD9OvXT7p16yZLly6VkSNHyl133SUTJ04szvYCQMBIPbhPJjx/l/zw+iMlXld4Vqa0//V7kY8+EsnM9Mr2AaGA4w7Aey7r2ljevuMsOad93ZAq1iY1c4Mga9asKXxBrV+1nqWuBWzN4xzuvn37mktRadaHzkM+evRo83/z5s1l8eLF8tJLL8kll1zi9jnp6enmYjl8+LCnmwkAPnck+aD8t2iWxJbJTcMtjh3//S17t66XWtXrymVjH82987LLRKKivLehQBDjuAPAybSonSTjxo2Tpk2bFr5gRobIDTfk3qauBWzL57PDLFy4UPr06ZPnvnPOOccEQjILONs5atQoKVeunPNSu3ZtX28mAHgsM/2ouY6OLd7MMGrZjIky9fVH5L+l82R1+24i/frl9lkGUCwcdwDuHUxJlxUb9kpaeuhlG5aJi5YhQ4ZIz549C19Q61etZ6lrAVvzeRBk586dUrVq1Tz36f9ZWVmyd+9et88ZMWKEHDp0yHnZsmWLrzcTADyWmX7EXEdGF3+E+ejYeHOdkpkhHz88VmTqVEasB0qA4w7Avf8t3yxt7v5cej3yrX2LSGeE0XqWuhawtVKZHUYHUHXl0Dm63dxviYmJMRcACIYgSFQJMkGi43KDIJlH07y2XYDdcdwBnGjNtoPmukmN3PEzQo0Oijp//nyTQd6rVy9/bw4AO2eCVKtWzZyVcbV7926JjIyUihUr+vrlAcD3QZCYkmeCZBAEAbyC4w7AvX+3HTDXTWtWCMki0hkrb7jhBnn//ff9vSkA7B4E6dy5s8yYMSPPfdOnT5dTTz1Vohj4D0AQyzxqBUGKnwkSdSwIEp6SLPfdcZ5I48YiaWSFAMXFcQfg3ppjQRBrJpVQ06RJE3P977//FryQ1q9az1LXArbmcRAkJSVFli1bZi7WFLh6e/Pmzc7xPAYOHOhcXgcp2rRpkwwbNkz++ecf+eCDD0yEdvjw4d58HwDgx4FRcwMZxWE9V7vDVNq5WeS//7TPoNe2EQh2HHcAJadd0a3uMKGaCWIFQXSaXKvr/Qn0fq1nqWsBW/N4TBCd1cV15GUNbqhBgwbJRx99JDt27HAGRFT9+vVl2rRpcu+998rYsWOlRo0a8tprrxU4PS4ABIv2514hbXpdLNnZxR9pPzouwVynZByVt579WIb0aMjAqIALjjuAktPJCA6mposOx9eoemhmgjRs2FDCw8NN4FS74levXt39wKjz5h2/DcCWwhwFhkoDx+HDh81UuTpTTNmyZf29OQAgIyb95ZVSOLxvp2xZ9aeUKV9Z6rY+TUZd3JrSRcgI1vo7WLcbKIgOGNq1a1epWyVRNr53Q2gW1IAx0qhRI1m3bp3MnDlTzjzzTH9vEYAArb9LZXYYAIB7ZStWk5bd+lM8AACfqVmzpowa2EUiI3w+HKDfu8RoEETHBSEIAqAgBEEAoJj+nPaFbF+7Qlp27y8N2nUtUTmGZ2dJi99/FcleLXLRRSKR/DwDALyjXr168tClp4Z8cTZt2lR+/PFHMy6IW1lZIt9+m3ubuhawLY6yAaCYtBvL6oXTpUbjU0TaFW8dWRnpsn7pfAlPPSTPjn08986UFIIgAAB4SKfI7du3r7Rt29b9AunpIpdfTl0L2BxBEAAopoyjqeY6Oq74s8Nkph+RSS/cIzo827oWHaRh5TIi4aGdrgwAKF06SUH9LfulSY3yEhHCXWJOOeUUcymQ1q89ehy/DcCWCIIAQDFlHE0z11ElmCLXeq5OtvvGQ6/Jq9d14fMAAHjNkSNHZMCAAZKTkyM7Ph4s1SrkzkpmS3FxIrNm+XsrAPgZIVAAKKaMI7lBkJhj09wWR0RklIRH5MajM48FVQAA8JbVq1ebAEhSYqxULV/8oH2wmDx5sjz88MMFjwsCwPbIBAGAYso4klriTJCwsDDTneZoymFnZgkAAN6ycuVKc92qTpKpc0Lda6+9Jr/88os0bNjQDJQKAPmRCQIAJQyClGRMEPP82HgzJsgjo+4S0cHcjhzhMwEAeMXff/9trlvWqWiLEm3ZsmWe4E8eWr9qPUtdC9gamSAAUEwZR3ODFdEl6A5jZZJkiUidHZtE9JKTw2cCAPBqEKRVXXsEQVq1alVwEETr1+XLj98GYEsEQQCgmO75eK7pwhKbULbEmSAHROTpK++QR2+8QCRW80IAACg5KxjQsk5SaBfnlLvNVcsDO8z134vnOe9zys4RefJCkc63U9cCNkYQBACK+wMaHWMuJdXl0pslPS1VdrY6VaR3bz4PAIBXpKSkyMaNG23VHaZF7dxgz7Z9qXIwJV3Kl3Gpp3V64HZ1qGsBmyMIAgB+1rhjT39vAgAgBEVERMiECRNk3bp1UqnsNrEDDXrUqlRGtu5NkVVb9kuX5tX9vUkAAgwDowJAMRzeu1OmjBkpMz95xTs/xtlZ0vTPOSJTp4pk6QghAACUTFxcnFxyySXywAMP2KooWx7LBlm5ed+J3WEWbaCuBWyOTBAAKIaUA3tk5ZwfpGyl6tJz4LASleG+rRskddMaefaVYwepKSkikfw8AwBQHK8M7iax0ZFSr0q+Mbsys0WenqIDiFDXAjbGUTYAFIMOiOqN6XHV8l8myorJn8iq8hWlRaP6IuEk6QEASu7LL7+UMmXKSNeuXaW8jQq0RUHjn4SFiTSqIlK+DnUtYGMcaQNAMaSnpZjrmPjEEpdfTEJZOSoiV596psiiRZq/zGcCACgRh8Mhd911lwwYMEDWrl1LaZoKN1LklSupawGbIwgCAMWQnppsrmPiy5S4/Kx1WIEVAABKatu2bbJnzx4zOGrr1q1tV6Avf7tELnt+mmzafdjfmwIgwBAEAYBiSE/zQRDkWGAFAICSWrJkiblu2bKlxMbG2q5Av5j7r0xY8J8sWrvL35sCIMAQBAGAYjh6LGARm1Dy7jC6Dj08/WrNMpEzzhA5coTPBADglSBI+/btbVmS7RpUNtdL1+85fmd6lsgD31DXAjbHwKgAUKLuMF4YEyQ+0USkO+pgqwsWiOTk8JkAAErE7kGQ9seCIEvWuQRBHA6R1TtEZAd1LWBjBEEAoBh6XjdMOl9yk0RERHmlO0y6iFwZnyhffvaJSEwMnwkAoETsHgSxMkGWrNttBokN05lhoiJERvYX6XgTdS1gYwRBAKAYIqKiJKFcAVPweahc5Rpy1k0j5WjZ8iIXXsjnAQAokV27dpmBUbXh36ZNG1uW5in1K0l4eJjsPnREduxPlRoVy4hEhIt0aigygLoWsDOCIADgZ7FlykqHvlf6ezMAACGiUqVK8tdff8maNWukTJmSD+AdjOJjoqRZzQqyast+My6ICYIAAEEQACieOV+8IUdTDsup/a+WpBr1SlyMYdnZUu+fJSJJ+0S6dROJiOCjAQAUi06L26pVK3OxM+0S89+Og7JtX2ruHdk5Iqu2iyTOoq4FbIxMEAAohlVzp8nBXVulRbd+klSj5EW4e+Viee7Jm3P/SUkRSUjgcwEAoATG3NxdPrjrbInWsUBUZrbIw5NEZBJ1LWBjBEEAoBjS01K8NkWumvL6w7JSRBrUry9xOngbAACFVhx3F/jQ7eNmSovaSTLorOaSGB9t23KsWDYu7x1avdZOEkmsJkJdC9iWzsoIAPCAjjJ/1DlFrnf6GOckJIomLS94912R+Hg+DwBAsWzZkyzjfvxL7nlvjhkYFC5iokTGXiuyciV1LWBjBEEAwENpaWniyMk2t2PivZMJEptQ1lwfPHiQzwMAUGwLVu8w120bVJaE2JJP4x7sXv52ibS563P5fPYaf28KgGAOgrz55ptSv359iY2NlQ4dOsjcuXMLXHbWrFlmeq78l9WrV5dkuwHAb6xARVh4hETF5ku1LSYro+TQoUNeWR8QSjjuADwPgnRpVp1iE5Ht+1Nlxca9Mk8HRAWA4gRBvvrqK7nnnnvk4YcflqVLl0q3bt2kb9++snnz5kKfp1N07dixw3lp3LgxHwCAoGQFKjRwoUFdbygXlyDTReSs558XOXLEK+sEQgHHHYBn5v9DEMRVl+bVjweH0rNEHv1WpHdv6lrAxjwOgrzyyisyePBguemmm6R58+YyevRoqV27towbN67Q51WpUkWqVavmvOjUXQVJT0+Xw4cP57kAQKDYt2+fuY4vW95r64xLKCu9RaTe2rUiOTleWy8Q7DjuAIou9WimLFu/J0/j3+46N61mrv/atE+SU9NFlm8R+d//qGsBG/MoCJKRkSF//vmn9OnTJ8/9+v+CBQsKfW67du2kevXq0qtXL5k5c2ahy44aNUrKlSvnvGiQBQACRZcuXeSuD2bJ5Y++5bV1RpWvKNeIyHtnnikSE+O19QLBjOMOwDOL1u6S7ByH1KyYILUreWfg7mBXo2IZqVslUXJyHLJo416RYX1Exo+nrgVszKMgyN69eyU7O1uqVq2a5379f+fOnW6fo4GPd955RyZOnCiTJk2Spk2bmkDInDlzCnydESNGmHRz67JlyxZPNhMAfEoz2eLLJUn5KjW9ts66HbrJnsEjpNnTT4tEMns5oDjuADyzdvtBiQgPkzOa1/Bad81QYI2PMv/fXSJnNhO55hrqWsDGinWknf9HVaeLLOiHVoMeerF07tzZBDVeeukl6d69u9vnxMTEmAsA2EX1Rq3MpWvX1v7eFCDgcNwBFM3N57SSq7o3kUNpGRRZviDIF3P+lTkrt1EuADzLBKlUqZI5A5o/62P37t0nZIcUplOnTrJW+70DQBAaP368TH/3OdmwfKHX1hmWnS21/vtbZNEikezc6XcBu+O4A/BcmbhoqVmRrjCuzjqlltSpnCjNa5YXWbuLuhawOY+CINHR0WZK3BkzZuS5X//XPvJFpbPKaDcZAAhG+pu35KcvZdeGf7y30tRkuePBq0VOO03k6FHvrRcIYhx3AEWnmdlwr3ntJNn43vXy2vVdRe77iroWsDmPu8MMGzZMrrvuOjn11FNN1xYd70Onxx0yZIhzPI9t27bJJ598Yv7X2WPq1asnLVu2NAOc6RlUHR9ELwAQzLPDxCVW8No6j6Qcko3HbtfVzBCvrRkIbhx3AEXz1Jd/yNTFG2X4Re3k8q5NKDZ3Xer0qkqiSFyS3kkZATblcRDkiiuuMA2Ap556Snbs2CGtWrWSadOmSd26etgu5j4Nilg08DF8+HATGImLizPBkKlTp0q/fv28+04AoBQHa1Rxid6bIjeyUjWpf+z2gcxM8d6ageDGcQdQND8t2WRmh0lOy6TICpAdGSH/PHq+tBo6njICbCzMEQS5c4cPHzZT5epMMWXLlvX35gCwuUaNGsm6devk2mc/llrN2nltvS9fc7pkHj1ixkzS1wCCXbDW38G63bCZKXc7bx5MSZeK175jpoHd9P4NZvwL5HUoNV3q3/yRHExNl71790lSUhJFBISYotbfHo0JAgDwTSaIWV+Z8nm62wAAUBQz/9pqAiBNapYnAFKAcgkxUr1Cgujp35kzZ7JjATZGEAQAPJCZmWmiy17vDpORLl+mJcu3GgTZxhR+AICim7Estyt67zZ1KLZC9G1dy9SzDe+7j0HIARsjCAIAxcgCCQsPl9gE76XJh+VkS9+0FLlQX2PXLj4TAECRaM/2aX/mDq3du21tSq0Q57SpZerZtps2SU4mY6cAdkUQBAA8UK1aNdm/f7/c9Oq3Eh4R4bWyy46MkiebtZWbRWTr7t18JgCAIlmybo9s2p0s8TGR0rsdmSCF6dG2jtwVHWHq2t+XLmUPA2zK49lhAMDu0+xVqFBBKtay5nLxjpzIKFl9+e3StUqGdOnSxavrBgCErujIcLmyWxOJi46U+Jgof29OQIuOi5Z9nRvJ57PXSPkpU6Rz9+7+3iQAfkAmCAAEiHptOskdd9wh7dp5b8YZAEBoa12vknxx/7nywd1n+3tTgsLFnRua60mTJpmuRADsh0wQAPDAl19+KXPnzpW9Sa2kQbuuXiu7sJwcqbx1vcjKcJHmzUXCiVEDAOBVOQ7pWzlRXr3pJukyeDCFC9gUR9kA4IFffvlF3nzzTdmxbpVXyy0y46jce+/FIq1aydQJE/hMAAAnNfvvrbJy8z4yGooqI0vi7/tK7nnvPTmtdWvTxRWA/RAEAQAPbN++3VyXqVDZ6+WWXKac7BGRq666SnJycvhcAAAF0q4ct46dKa2GfiYT5v9HSRVV2ViRSpUoL8DGCIIAQDGCIIlJ3g2CZMbGyzMfzJSqYWGSnJMje/ZoOAQAAPfmrtwua7YdkITYKDm3fV2KqShio0TG3yKyZ4/MW7pUrr32Wpk5cyZlB9gMQRAA8MC2bdvMdZkKVbxebuERkVK1atU8rwMAgDvvTv/bXF/VvYkkxkdTSB764osv5LPPPpO3336bsgNshoFRAaCIjhw54szQSKzo/SCIqlOnjuzcuVM2btwo7du357MBAJxg//798s2xLjC3nNOKEvLUlLvlpsb75U0R+XbiN7L3s7JSqWyc+2UHjKF8gRBDJggAFNH69evNdfny5SW2TDmvlltkRrpcMfoheWXXLokRkXXr1vG5AADcGj9+vKRnZkub+pXk1Ea+CcqHpIwskZd/Npd2tZOkfcPKkpGVI5/OXO3vLQNQigiCAEARbdiwwVw3bNjQ6yPKh+VkS9u50+SMTZskQkT++49B7gAAJ8rOzpaxY8c6s0CY4cQDOQ6R2WtyLzkOublPbhbNa1OWS2ZWNrsbYBMEQQCgiM477zw5cOCAfPPNN14vs+zIKPnhhvvlj6uukgwyQQAAhWQlpqSkSFJirFx3ZjPKyROR4SKDu+VeIsNl4FnNpEq5ONm4+7B8PvtfyhKwCcYEAQAPaFcYvcjSv7xabjmRUTL/vOvkpjbx8mqXLtK6dWs+FwDACRo3bmy6TK4cdwMDonoqMkLkgnbOf+MjI+S+C9vLgx/Pl+cmLJJrz2wqERGcIwZCHUEQAAgg2tVm6NCh/t4MAEAAi42NlQ6MBeIVt/VtLZ/PWSODe7eU7ByHRGifVAAhjVAnABTRwIEDTYBi+/btXi+zsJwcKb97m8jGjSI5OXwmAIA80tPT5YMPPjBjgqAEY4LsOpx70ds621t8tCwdfZXceV4biY4iAgLYAZkgAFAEhw4dkk8//dTcfuqpp7z/Y5xxVB68ra+5vX7FClmwfLk0b95cOnTowOcDAHYx5e4CH/q/L3+Xxz//Xb5962mZ8uj5pbpZITU7zM0f5d7++jaR2Chz03Vw2ezsHLrEACGOTBAAKILly5eb6zp16khSUpJPyiwjJlYkPl7eeOMNue666+SLL77gswEAyNrtB+W5bxabkrimB4OhlkhMZO4lH4fDIZMW/CdNb/tU1mw9wF4HhDCCIABQBMuWLTPXbdu29Ul5ZcbGy+Of/yGSmiotTzvN3LdkyRI+GwCwOZ26dfBr/5P0zGzp066OXNGtsb83KXhp5sc3t+dejmWBuHpvxkpZt/OQ3PDaDMnIpNsREKoIggBAESxYsMBcn3rqqT4vr/bt2zuDIHpmCgBgXw98NF/mrtouZeKi5M0hZ+bpugHv0XIde+uZUi4hWhau3inD3p9L8QIhiiAIAJyEBiLmzs09GOrWrZvPy6tly5YSHx9vxiGxuuEAAOzn05n/yOjJuZmIH9/dWxpWL+/vTQpp9auVk/HDzjG3x05bIR/MWOnvTQLgAwRBAOAk1qxZY2aEiYqKktNPP90n5RWRmSEXjXtC5OabJdrhkLPOOsvc/+OPP/L5AIANpaVnykMf52YhPnx5R7m4SyN/b1Lwy8wSeeOX3IveduO8jvXl0Styu6XePPZX56DoAEIHQRAAOInNmzdLzZo1TWAiLi7ONz/G2Vly2v8mibz3nkhWlvTtmztTDEEQALCn+Jgo+fHxC8zUrU9e5ZsAvO1kO0Smr8y96O0CPHHV6XJTn5aSk+OQQYMGyYoVK0p1MwH4FlPkAsBJ9OnTR7Zs2SIHDvhutPiciCj5+aqh5vbcH9bIvvAG5vYffy6Vez+dL7EJZd0+b9TFrX22TQCA0pWcliEL1+yQPu3qmv9PqV9JXrulBx+Dt0SEi1zb+fjtAoSHh8nbt58lURHhsiehqbRuTV0LhJIwRxCMunf48GEpV66c6R9ftqz7hgAAlMSISX8FXAGuXjhd6p3SqcAAiCIIgkAWrPV3sG43gld2drZ89NFH8tgDd8vew0fktxcvl3YNq/h7s2xPm0lHM7Il7tiUulv2JMtPSzbJjWe3kIiCgigDxti+3IBAr7+L1R3mzTfflPr160tsbKx06NDBOWBgQWbPnm2W0+UbNGggb731VnFeFgBKVXZmpqz49Xtz7Q/NOvcpNAAC2AXHHQhV27Ztk+eff14aN24sN910k2zfnyo1K5aRjKwcf28ajs0YYwVA1MPjF8otY3+Vhrd+LM99vcgERQDYoDvMV199Jffcc485IDnjjDPk7bffNn3XV61aJXXq1Dlh+Q0bNki/fv3k5ptvlvHjx8v8+fPl9ttvl8qVK8sll1zirfcBAF6VnpYiP457UlYv+FmW/vyVDBw1XsLCfTiMksMhCYdzu9uklq2gR17HH8rJkTlfvC6tepwvFWvV9902AAGI4w4EvSl3mysdX0K7Wag9h9Kk92PfyfINe52LJSXGysOXnSp39D9FYqLose4TmgB/+Eju7bJxeeraoji1URWZtnijbNqdbAIiemlTv5L0aFVTOjWpJlf1aOqb7Qbg3+4wOjNC+/btZdy4cc77mjdvLhdeeKGMGjXqhOUffPBBmTx5svzzzz/O+4YMGWKmfVy4cGGRXpO0VADF8fvvv5tBTXNyckxKq167Xq655hoz44u6+qkPZN/W9ZJ2+IDs37ZB1i2dJ+mpyRIeESmXPDhGGnbw7dS4UUfT5KlrOpnbj332m2TGxjsfWzDhHZnzxRtmW+q2Pk2qNWghiRWrSERUjFx7RiPp3bu3CSxbM9n89VfBXXu6d+8uVarkplivX79eli3LnXrRnc6dO0v16tXN7U2bNsmff/7pdrmmTZuaaX0BX9TfHHfYV1pamhkc2jpUzX+t2cWnnnqquZ2eni6TJk0qcFnNYNaTd1bXEz0xV9CydevWlV69ejm347333jN1hrWM6/K1a9eW8847z7nsY489JikpKZKammq2f//+/bLtnz9k274U6daipkwa2f/YNuRI4hVvydHMLOnSrLoZhPPyro3NYKjwoaOZIpcfa8N8fZtIrOflfSQ9S76Zv1bem75S5q/eYYJbql2DyrJk9FXO7jAXXXSR2S/Lly9vLvpbqIOr63FHvXr15KqrrnKuc8KECWZZ1wwU6zopKUnOOSd32l71ww8/mH3Letx1Wf2d1WMCy4wZMyQ52X22Snx8vJx77rnO/3/99Vc5ePCg22VjYmKkf//cfVfNmTNH9u49HsBzFRERIRdccIHzfz0BvmvXrgLL8+KLL85z3KbZUQU5//zzJTIyN0C4ePFic4xXEN1e3W61dOlSc2K+IFq+CQkJ5rYeQ61du7bAZc8++2xnfabtW9c2rraR9bNF4B93eBRmzsjIMAfBDz300AmDBi5YkDuFV34a6NDH8+9o77//vmRmZjobIK70R8D1h0DfhPWmAKCoXnjhBXNQXJSKbPkvk2TVnKl5Hi9ftbb0HvyQ1GrezmSG+FLO0SNi/cKlp6VK5rEDbtW0c2/ZsupP2bB8oWxYtsBcLD++mXuQc9ppudP5ffPNN/Loo48W+DoalO7RI3eQve+++07uu+++Qs/AWwdIP/30kwlgu6Pr0AN/ID+r3i7u8GMcd9ibDkh96aWXFvj4DTfcIKNHjza3Ndhw9dVXF7js5Zdf7hzc8ujRo3L99dcXuKwGNTp27Oj8X7OZC9LrlFrSPfUH5/8vv/C+pKW7n3p10+7Dcjjt+PHtp8N6m4ZzJc1IEJGs7Jw8j8MHjrp8NmkZmp5TrNVc2KmBuew7fFR+WbFVlqzbLTWSEnI/v6+GmNX+MOV7ySpgBpozmleX/uGzRfq+YP6/9dZbzT7sjjas9aSE5bbbbpOtW7cWeFLijz/+cP4/dOhQ+ffff90uqwG8v//+2/n/8OHDTbDAnYoVK5oTJ5aRI0ea4IY7GujZuXOn8/8nnnhC/ve//0lBrHae0hPq33//fYHL7tixwwRv1CuvvCJffPFFgcuuW7dOKlWqZG6PHTvWtD0LorP/aPBTvfvuu/L6668XuKwGapo1a2Zuf/LJJ6Y7m0Vf59prry3wuQig4w6HB7Zt26Zrc8yfPz/P/c8++6yjSZMmbp/TuHFj87grfb6uZ/v27W6f8/jjj5vHuVAG7APsA+wD7APsA8G/D2zZssWTww2OOwLgM+NCGbAPsA+wD7APSIgedxSrw6GVduUSSDnhvpMt7+5+y4gRI2TYsGHO/zUFUaOjGoUs7HXsEt3SyK2eHWHEeso6FLBPU86hhP1ZTqjvNRW7Ro0aJSrXUDruYB+hrHyB/YqyYr/yL76DgVFeRT3u8CgIoilF2s/LNcVJ7d69W6pWrer2OdWqVXO7vPbn0oMLd7T/ltWHy6J96XCc7iwEQUoHZU05hxL2Z8q5tGnf3OIK5eMOvouUFfuVf/EdpKzYr0Lze1iU4w6PpjqIjo42U91q/3NX+n+XLl3cPkf7sOVffvr06WYgK3fjgQAAAHDcAQAAfMHj+R41XVRHyf7ggw/MaLj33nuvGZnXGjBPU0oHDhzoXF7v11kF9Hm6vD5PB6bRwXcAAAA47gAAAKXF4zFBrrjiCtm3b5889dRTZoTeVq1aybRp05wj6up9rtMV6ZRk+rgGS3TEXO2f89prr8kll1zi3XdiE5qu+/jjj5+QtgvKOlixT1POoYT92ftC7biDfYSyYr/iOxgs+L2irEJ13wrT0VH98soAAAAAAACB3B0GAAAAAAAgGBEEAQAAAAAAtkAQBAAAAAAA2AJBEAAAAAAAYAsEQUJEenq6tG3bVsLCwmTZsmX+3pyQsnHjRhk8eLCZcSAuLk4aNmxoRjPOyMjw96YFvTfffNOUa2xsrHTo0EHmzp3r700KOaNGjZKOHTtKYmKiVKlSRS688EJZs2aNvzcr5Mtcf4vvuecef28KguS3bfbs2WY5Xb5Bgwby1ltviV14UlazZs0y3638l9WrV0uomzNnjgwYMMDMdqTv+bvvvjvpc+y6X3laVnbdr4p7fGDH/ao4ZWXX/UqNGzdOTjnlFClbtqy5dO7cWX788UcJpP2KIEiIeOCBB8yPPbxPf6xycnLk7bfflpUrV8qrr75qvpgjR46kuEvgq6++Mo3Ehx9+WJYuXSrdunWTvn375pnqEiWnlcodd9whv/32m8yYMUOysrKkT58+kpqaSvH6wKJFi+Sdd94xlT/sydPftg0bNki/fv3Mcrq81i133XWXTJw4UUJdcesBbXzo1MjWpXHjxhLq9De7TZs28sYbbxRpeTvvV56WlV33q+IcH9h1vyrJsZTd9itVq1Ytef7552Xx4sXmctZZZ8kFF1xg2lEBs1/pFLkIbtOmTXM0a9bMsXLlSp3u2LF06VJ/b1LIe+GFFxz169f392YEtdNOO80xZMiQPPfpfvzQQw/5bZvsYPfu3eZ3Yvbs2f7elJCTnJzsaNy4sWPGjBmOHj16OO6++25/bxKC4LftgQceMI+7uvXWWx2dOnVyhDpPy2rmzJnm9+vAgQMOO9My+Pbbbwtdxs77ladlxX5V9OMD9quilxX7VV4VKlRwvPfee45A2a/IBAlyu3btkptvvlk+/fRTiY+P9/fm2MahQ4ckKSnJ35sRtLQr0Z9//mmi6K70/wULFvhtu+yy7yr2X+/Ts0T9+/eXs88+2wdrR6j+ti1cuPCE5c855xxz9iwzM1NCVUnqgXbt2kn16tWlV69eMnPmTB9vaXCy635VEnbfr4pyfMB+VfSysth9v8rOzpYvv/zSZM1ot5hA2a8IggQxDXBff/31MmTIEDn11FP9vTm2sW7dOnn99ddNuaN49u7da34Uq1atmud+/X/nzp0Uqw9/M4YNGyZdu3aVVq1aUc5epBX8kiVLTL9h2Fdxftv0fnfLa7q1ri9UFaestCGh3c00RXrSpEnStGlT07DQMSCQl133q+Jgvyr68QH7VdHLyu771V9//SVlypSRmJgY02b69ttvpUWLFgGzX0X6ZK0okSeeeEKefPLJk/Y71zMlhw8flhEjRlDiPixn1wDT9u3b5dxzz5XLLrtMbrrpJsq9hHSAqPwVS/774D1Dhw6VFStWyLx58yhWL9qyZYvcfffdMn36dDOgF+Dpb5u75d3db/ey0kaEXix6VlG/fy+99JJ0797d59sabOy8X3mC/cqz4wO771dFLSu771dNmzY1k3UcPHjQBIIGDRpkxlYpKBBS2vsVQZAA/XJdeeWVhS5Tr149eeaZZ8wAPRphc6WN9muuuUY+/vhjH2+pPcrZNQDSs2dP8yOmkV0UX6VKlSQiIuKEs327d+8+IRIM77jzzjtl8uTJ5gyEDlgF79GUft13dVRzi57h1rLWgfl09i7d3xH6ivPbVq1aNbfLR0ZGSsWKFSVUease6NSpk4wfP94HWxjc7LpfeYud9itPjg/svl+V9FjKTvtVdHS0NGrUyNk21RPLY8aMMRNNBMJ+RRAkQA8M9HIyr732mgmEuDbStf+UjrZ++umn+3gr7VPOatu2bSYAoo2cDz/8UMLD6UlW0h9GLUsdYfuiiy5y3q//6+jR8B6NpGulrWmIOl2bTkUJ79L0Vk37dHXDDTdIs2bN5MEHHyQAYiPF+W3TwPqUKVPy3KdZRXrQGBUVJaHKW/WAziSgaefIy677lbfYYb8qzvGBXfcrbx1L2WG/KqwM9aRQwOxXPhtyFaVuw4YNzA7jA9u2bXM0atTIcdZZZzm2bt3q2LFjh/OC4vvyyy8dUVFRjvfff9+xatUqxz333ONISEhwbNy4kWL1ottuu81Rrlw5x6xZs/Lsu2lpaZSzDzE7jH2d7LdNZz657rrrnMuvX7/eER8f77j33nvN8vo8ff6ECRMcoc7Tsnr11VfNTB///vuv4++//zaP66HsxIkTHXaYfUpn/9OLvudXXnnF3N60aZN5nP2q+GVl1/2qKMcH7FfFLyu77ldqxIgRjjlz5pi26YoVKxwjR450hIeHO6ZPnx4w+xVBkBBCEMQ3PvzwQ/Oj5e6Ckhk7dqyjbt26jujoaEf79u2ZttUHCtp3db+G7xAEsbfCftsGDRpk9g9XemDdrl07s3y9evUc48aNc9iFJ2X1f//3f46GDRs6YmNjzXSLXbt2dUydOtVhB9Z0m/kvWkaK/ar4ZWXX/aooxwfsV8UvK7vuV+rGG290/q5XrlzZ0atXL2cAJFD2qzD945scEwAAAAAAgMDBwAYAAAAAAMAWCIIAAAAAAABbIAgCAAAAAABsgSAIAAAAAACwBYIgAAAAAADAFgiCAAAAAAAAWyAIAgAAAAAAbIEgCAAAAAAAsAWCIAAAAAAAwBYIggAAAAAAAFsgCAKg1M2bN0+ioqIkPT3ded+GDRskLCxMNm3axCcCAAA49gDgEwRBAJS6ZcuWSfPmzSUmJibPfeXLl5e6devyiQAAAI49APgEQRAApW758uXSrl27PPdpEKRNmzZ8GgAAgGMPAD5DEARAqdOAR9u2bfPct3TpUoIgAACAYw8APkUQBECpys7OlpUrV56QCbJkyZITAiMAAAAcewDwJoIgAErVmjVr5MiRI1KjRg3nfQsXLpRt27aRCQIAADj2AOBTBEEAlHpXGPX666/L2rVr5ccff5SBAwea+1xniwEAAODYA4C3EQQBUOpBkN69e5spcVu1aiUjR46U559/XsqWLStjx47l0wAAABx7APCZMIfD4fDd6gEgr3POOUfat28vo0aNomgAAIDPcewBwBWZIABKfXrcU045hVIHAAAcewAodQRBAJSanTt3yq5duwiCAAAAjj0A+AXdYQAAAAAAgC2QCQIAAAAAAGyBIAgAAAAAALAFgiAAAAAAAMAWCIIAAAAAAABbIAgCAAAAAABsgSAIAAAAAACwBYIgAAAAAADAFgiCAAAAAAAAWyAIAgAAAAAAbIEgCAAAAAAAEDv4fzoD1yr7/dBjAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Analytic posterior on mean for comparison (sigma plugged in as sample std)\n", "post_mean, post_sigma = gs.analytic_posterior(x_obs)\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "\n", "mu_grid = np.linspace(0.5, 2.5, 300)\n", "\n", "for ax, samples, w, label, color in [\n", " (axes[0], proposal_broad, w_broad, 'broad prior', 'C0'),\n", " (axes[1], proposal_tight, w_tight, 'tight prior', 'C1'),\n", "]:\n", " ess = eff_sample_size(w)\n", " ax.hist(samples[:, 0], weights=w, bins=40, density=True,\n", " alpha=0.6, color=color, label=f'IS ({label}, ESS={ess:.0f})')\n", " ax.plot(mu_grid, norm.pdf(mu_grid, post_mean, post_sigma),\n", " 'k--', lw=1.5, label='analytic posterior')\n", " ax.axvline(theta_true[0], c='red', lw=1.5, ls=':', label='true mean')\n", " ax.set(xlabel=r'$\\mu$', title=f'IS — {label}')\n", " ax.legend(fontsize=8)\n", "\n", "plt.suptitle('Importance Sampling: posterior on mean', y=1.01)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Key takeaway:** the broad uniform proposal samples mostly from regions of low likelihood,\n", "yielding very few effective samples. A tighter proposal that overlaps the posterior better\n", "gives a far higher ESS and more accurate estimate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 3 Effective Sample Size as a diagnostic" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArIAAAFUCAYAAADYjN+CAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATQRJREFUeJzt3QmcVfP/x/HPLDXt+yYqS2lHKmnRRvn98ksJIVL2pZBQCGWpFBKyhuRHhJ/9b8lPRCpaJC0qKhUlbdM+NTP3/3h/69zfvXf2Zqa5Z+b19Dju3HPPPds9t/s5n/P5fk9MIBAIGAAAAOAzsQW9AgAAAMDhIJAFAACALxHIAgAAwJcIZAEAAOBLBLIAAADwJQJZAAAA+BKBLAAAAHyJQBYAAAC+RCALAAAAXyKQBTLxyiuvWExMTIbD119/HZx29+7dNmbMGDv55JOtXLlyVrZsWTvhhBOsd+/eNmPGjLD5Llu2zPr27WvHH3+8lShRwqpUqWKnnnqqDRw40Hbs2MFnEiVGjBjhPuf8cuyxx1r//v2znE7HWeTxNmXKFBs/fnyaadesWeOmffTRR3O1bq+++qpVrVrVdu7cGba+mvf111+f4Tq+8847llvvvvuuXXLJJVa3bl0rWbKkW+6ll15qK1euTHf6//73v9a6dWsrVaqU+y5pn27atCnNdAcOHLD777/fzS8hIcEaNGhgTz31VLrzXLVqlfXq1csqVKhgZcqUsS5dutiCBQss2uzZs8cdp6HHRuS/Xzomoum4b9++vQ0aNCjf5o+iJb6gVwDwg0mTJrkfvUiNGjVyjykpKda1a1f7+eef7Y477rDTTjvNjdcP70cffWTffvutdejQwY378ccfrW3bttawYUO777773I/q5s2b7aeffrI333zTbr/9dhcIo/B77733DvuzViC7ePHifAkIFBzdfffdNnToUHdCFumll16yW2+91erXr2/5QSeENWrUsGHDhrmTvXXr1tmoUaPcyd6cOXOscePGwWl1kvjPf/7TzjnnHPvggw9cAKv1PvPMM23evHkuYPXceOON9u9//9sefPBBa9mypX3++ed2yy23uGBd2+v5+++/7YwzzrCKFSvayy+/7E42R48ebR07drS5c+fm23Yf7mel4Fy0fqG0T2bPnm1HHXWURRPtf50Y3HDDDVG1L+FTAQAZmjRpUkBfk7lz52a6l6ZPn+6me/nll9N9PSUlJfj35ZdfHihdunRgx44d6U6bmprq209E675nz55AYTF8+HD3uRa0r776yq2HHj3nnHNOoE6dOmmmXb16tZv2kUceOezlPfPMM4ESJUoEtm3bFjZey2vdunWgfPnygV69eqW7jm+//XYgt/7666804/74449AsWLFAldddVXY+JYtWwYaNWoUOHDgQHDcd99959ZF2+FZvHhxICYmJjBq1Kiw919zzTWBkiVLBrZs2RIcd8cdd7hlrVmzJjguMTExUKVKlUDv3r0D0fRd+/vvv9226lj103HfpEkTt++B3KK0AMgDW7ZscY8ZZT5iY2PDplUWTpcr05PZJb3333/fvf7ll1+mee3ZZ591ry1atCh4afTiiy+2mjVruqxU9erVXZZq4cKFmW6LLstq3ZYsWeKmL126tLvErLIHZX8i11Xjn3vuOZdh1nImT57sXps5c6Z7vzJ6uuTbpk0b+7//+790L31+8cUXdsUVV1ilSpXc8rp37+7WP5KyYyrdUIZM05533nmuTCNUdrZ76tSpLoOuz0uXrrXud955pysPySltk7ZBmTrPf/7zHzdOGbFQJ510kp1//vmZlhb88ssv9o9//CN4mVyX8UMv73uZNy33999/Dyt1iTRu3Dg77rjj3OepS+/KZmaHjiV9BrqsHkn7XftKl/+zO7+cqlatWppx+jyPOeYYl531/PHHH26/q0wnPv5/Fxh1rJ144oku4x363QkEAu44C6Xne/futc8++yw4Tu/r3Lmz1alTJzhO31mVGugKS3Jycqbrr8/1X//6l5uPPnMdr8osP/nkk2HT7du3z2677TY75ZRTrHz58m7f6nNSZjlSRt81fTdFWVnvOPCOqYxKC7St+k5omTrOND9lnLOi743WT99RHVNnn322u8IUKrv/7ugz01WFyGMbyCkCWSAbVDqgH6/QQeM8LVq0sGLFirnLlK+//rpt2LAhw3nph0Cvq+ZPl0X1I5pd+nHUj7xKHSLpR0uXXvXDKd26dbP58+fb2LFjXaCo4KRZs2a2ffv2LJejWkK9Xz9ACgD0A/r888/bRRddlGZava55q0xCl2p1SVbbpUAgMTHRXYZ+4403XECr4Eg/hpGuuuoqF+x7dZ8//PCDC9ZC11U/tJpOl5UVRD3xxBMuaNf+DK2dzM52a3pNp3XTj7ouz7/11ltu/XJKJSP67FWn6dHfCpC1H7QvRZe8VQpw1llnZTivv/76y81P0z3zzDPuMviuXbvc/g+l11SeosvvunTsDaGefvppt/3anzomFaRrm/WZZGb9+vWuRKZTp04ZTqPj/Oijj7YhQ4ZkOi8FjpHfm4yGrChAUuAeWlag/STeMR9K47zXvWkV9GmfRU4XOi99H3/77bcM56nX0zvJiqTATceVSjAU0Cq41n4LrV1OSkqyrVu3unIifY/0PWnXrp0LmFWjnNV3Tce+F4Dru+EdB/fee2+G66VjXsdBamqqC4oVmN98883uc8+MSjtUt6xyKn1XdGwqCNX3fenSpcHpsvvvjr7fOibTq+0FciTXOV2gCJQWpDfExcWFTfvSSy8FypQpE3z9qKOOcmUE33zzTdh0+/btC/Ts2TNsPs2aNQsMGzYssGnTpizXafDgwe5S6Pbt24Pjli5d6ub11FNPueebN292z8ePH5/jbe7Xr5977xNPPBE2fuTIkW78zJkzg+P0XJeZt27dGjbt6aefHqhWrVpg586dwXHJycnucuIxxxwTLJ/w9u95550X9n7v0vBDDz3knusSt7a5W7duYdOtXbs2kJCQEOjTp89hb7fWRZelZ8yY4d77008/5fgSa7t27QKdO3cOPq9bt667PB0bG+vmK6+//rqb14oVK8Iu1Wt/e4YOHeoufy9cuDBs/l26dMlxaUHTpk3dPvf88MMPbvwbb7yR6bZMnTrVTTdnzpw0r2l5Wq5MnDjRTffRRx9lWFrgjcvOoPXOiD6fjh07BsqVK+c+c4+3T2fPnp3mPddee22gePHiYfuwfv366c5f02l6r4RB8xw9enSa6aZMmeJemzVrVobr6u2njD5HbcPu3bvTfZ8+L22ryif0b0KojL5rmZUWeN8vb9/q+6jl63jNrIQp8rjXPo+Pjw/cdNNNYdNpfjVq1AiWW+Tk+7d//363j3TMA7lBRhbIBmVHdAkzdPj+++/DprnyyitdVkNZRWU4atWqZa+99prLsD3yyCPB6XS5TRkaZTEef/xxdxlOjUtGjhzpLvEtX74803XRcpQVCs1sKkOr+fbp08c91yVK9Zig5erysi7/KQOTE8oYh/Lm/dVXX4WNV+ZVjWI8yrJo31xwwQVh5RNxcXHucqL2UeQ2Ri5L2Std1vWWpSyTtjnyMrz2sZbvlVpkd7uVUdP2KDun9VJG1WuMF1mqkB3KXH/33XduHZU1/PXXX93nqkvGykp5WdratWtbvXr1MpyPtlcZR5VPpLfvc0JlDdo2j5dh1Ppl5s8//8zw8n7kJXll51RmkNGx1bx58zTfm4wGXYpOj2I4ZRvVYFLfQ33m2S3HiRyfWdlObqZNT0afo3olCe394O2333bZdX1XVB6hY1FZ0/SOw8jvWk7NmjXLLV+N3nLSK4Gyv8qaX3755WFZdJVM6HvjZVVz8u+OtlOlKyoPAXKDQBbIBgWYKh8IHfQjHUk1Z7r8psveCuZ06Vs1Ymp9HXlpTfPUpUcFu2vXrnX/8Kt+NrPLgt4PpFpce+UFKnHQPHr06OF+SMSro1UNmy7xqeRAl1UVYGenJk0/qJUrVw4b512S9eqBPZF1wdu2bXPBR3r1wl6wEjmPyMu93jhvusxqkDVP7/XsbLcu1etyqD6fhx56yP0IK5BSuYLkpNTDo3IBXSZWXbACV9W26nKqxnslB1qvzMoKvO3MaF/kVOTn57Xez2r7vNcVpGRGQbIuN6uW2quLjqTgTMF8dobixYuneb+Oo6uvvtod3yqd0TGe3jZGHk+iS/be98GbNr3pdOK1f//+4LQKFHUcZTRPCZ1vRjL7HL1565hT93wq09A26oRNx6JOVlU/Gym3vQ/ohFlUa5wTKnkR/bujADR00Am1el05nH93dIwdzvcNCEX3W0A+UtCpzJzqFFesWBHsliuSfgBUS/fAAw+E1fVllg1TVkVZG2UXVXMb2YhFGU1ldkTLVl2b+ofUj7Zq4zKjbIt+bEODoY0bN6YbIEVmdhQIqN41vTphL9unQC+UN+/IcepHNHSZGc0zdH5Zbff06dPdexTAellYyU7tcEZatWrlgjYFrWpYowyt9oseH3vsMRec6GQlq0BW25nRvjhSvH2poC2rwEmBpbKJw4cPtxdeeCHN66oRzqzWNtTq1atdI6nIIFYnbPo8L7vssjTvadKkiXtUTa9qM0NpnPe6NG3a1HVvp30ZGmRqutB5qbZZx503PnKeel0Nt7KS2efoHc8KXtUYT8Fg6PdIJ0XpyW3frl7DsKzqYTM6JtRHcGgDuPTk5N8dnfRG/lsA5BQZWSAPKOjTP9TpUSv00GxkRg3BFFzpsl9Gl1hDKeurbIayVBqU0VEr/IyoBfc999zjfsyz26m7GgiFUslEen1VRlKLZgV2yjaFZlt0iVE/3MoGaX0yW5YugeoSuLcsNWpRAKH3h9IPsgJTBYzZ3W4vGAjtX1TUmO1wKTOlTt6VjdX6qI9MUeZX2W2tgxfYZkZBnzKc6lM4vX0fSuufH9ksr79kNXjKbp+v6kkgskV+bkoLFMRec801LojV5xJ5kubRca+TQx0XoY0v1ZuCylfUaCo06NZnEJk91vdHx5Z6ivCoNwx9jqE9JCijqGP63HPPDeshISMZfY5q9KhMpWh9lIkODVAV7KbXa0FGsptp90p2dNVIAeXBstvsUYZV26xjIvLKlDfk9N8d/XunrLPXFzdwuMjIAtmgLGl6LatVD6Ysh2ob1SJZtZ76sVDGRa3U1QpZrYpVW+Zdzrv22mtd9k/dMCkLpEu0CnZVL6tMpjpzz4pqy/Rjqx9hzUutnkO7+FJJg1q6X3jhha4mUz+W+mHWeNU0ZkXTK5Ooy/C6nKjAUpfh1fG8WlVnRT0MKJhTYKZ10/zU0l77UfskMrOkjuuVfdP6KnhQKYaCFGWdve1VyYU6rde+VCCvkwd1OaSAXhnB7G63Ph9ljdWtld6nIFSBdGTQkVMKUtWVkniZVwVIWt60adNcjWpWdacqNVEXY6pv1f5WWYrWzTsZCqXgQIGVWoUrYNTnn1FAkRM6CdF6KxhU0JYVZWQVJKYXfCloO5x10qVoZfV0iV3bGdrNlwI3lW2EBtI61vSZ63jR906ftb5boQGwro6o1lafub5zOq71uSiTrH0dWi6gY1at8vU56CqJlvnwww+7wEvZxexQYK79p+mV2VawrRMdra+6vPJ6IdFnqPVWTbmOfd0sQNNndBez9PaxsqDa/zoGtR3KcoZmtz26aqDvtb5rOkZ1sqBjTDXdOv4nTJiQ7jI0L+0HfS91BUhBv75DKjlQDyM6edV3MSf/7nifaXYz9kCGctVUDCjCvRZoUMttWbduXeCee+4JtG3b1rXiVQvfsmXLBlq1auV6EghtPf75558HrrzySteJu1oha1r1cKAO5tNrfZ2RadOmBdcjtCW816F8//79Aw0aNHA3X1BvCieddFLg8ccfD1uX9KgVvd6zaNEi11JcvQVUqlQpcMMNNwR27doVNq2WPWDAgHTn8+2337qW/JqX5qGeDLwW7pH7V9vSt2/fQIUKFYK9E6xcuTLNPF988UW3HWplrn3Xo0ePwJIlS3K83Wp1ro79S5UqFahatWrg6quvDixYsMCti9bpcDqGV28HmrZevXrp9vag3iYiRfZa4PVAodbtuiGB9rtasH/wwQdpei1Q6/ULLrjA7TO1/vbWM7MbImS343x9Fjo+M+u1IHKd1ftGXt0QQcvJ6DuXXk8NOn50fHn7TL2FpHdTBbWU1/bXrl3bHUMnnnhi4Mknn0x3HX799VfXu4ha+es4OfPMMwPz58/P9vprP73zzjuBxo0bu2Ude+yxgXHjxqWZ9uGHH3avqfeNhg0bun9T0jvuMvuu/fe//3W9HGgems47piJ7LfB88skngQ4dOrjviLZNn/WYMWOyPO7ff//9QKdOndw+0bK0nToGtfyc/rujY0w9awC5FaP/ZRzmAihq1DOAauGUjc1vyigra6ZLy3mRTUTeUIZcGUtlzZShRc4og6mM8Mcff8yuS4dXQqWrUMoKA7lBjSwAIIxOKtSaXpe5gbymAFZd0WVU+wzkBIEsACAN1VIqK8stRJHXdLtfXY3JTqM5ICuUFgAAAMCXyMgCAADAlwhkAQAA4EsEsgAAAPAlKq0P3XFIdxlRx9K5vQUgAAAADp96hlVDU3XTFnqzn/QQyB66VV6tWrVyscsBAACQl3S3O++umBkhkD10iz9vh6lbEAAAABTcTTOUYPTis8wQyKoPskPlBApiCWQBAAAKXnbKPWnsBQAAAF8ikAUAAIAvEcgCAADAlwhkAQAA4EsFGsh+88031r17d9dPmAp633///TT9iI0YMcK9XrJkSevYsaMtWbIkbJqkpCS76aabrEqVKla6dGk799xzbf369Ud4SwAAAFCkAtndu3fbySefbBMmTEj39bFjx9q4cePc63PnzrUaNWpYly5dXCe5nkGDBtl7771nb775ps2cOdN27dpl//rXvywlJeUIbgkAAACOtJiA0p5RQBlZBaQ9e/Z0z7VaysQqUB06dGgw+1q9enUbM2aMXXfddZaYmGhVq1a1f//733bRRReF3dzgk08+sbPPPjvb/ZWVL1/ezY/utwAAAApOTuKyqK2RXb16tW3cuNG6du0aHJeQkGAdOnSwWbNmuefz58+3AwcOhE2j4LdJkybBaQAAAFA4Re0NERTEijKwofT8999/D05TvHhxq1ixYpppvPenR5ldDaGRf+ijFCtWzNXl7t271wXLocG0BpVFhJYvlChRwq2LShtSU1OD40uVKmXx8fFh8xbV8+r+waFlEqK7WOj9mn8onZEkJyfbnj17guP0/jJlytj+/ftt3759wfFxcXFu/pHbyTbxOXHs8X3i3wj+Lef3id/caI8jIsdnKhAltCrvvfde8Pl3333nxv35559h01199dWBs88+2/39+uuvB4oXL55mXmeddVbguuuuy3BZw4cPd/PObLjqqqvctHoMHa/3SteuXcPGT5w40Y1v1KhR2PjPPvvMjS9btmzY+MWLFwcSExPTLFfj9FroOL1XNK/Q8VqWaNmh47Vu6W0n28TnxLHH94l/I/i3nN8nfnPNJ3GE/s5K1NbIrlq1yk444QRbsGCBNWvWLDhdjx49rEKFCjZ58mSbPn26nXnmmbZ169awrKwakGk+999/f7YzsqqrXbduXbAWg+wl2Ushc87VgILOTHDVhitRHHt8n4ravxF//PGHi8uyUyMb9Y29br31VhsyZIgbpx1drVq1NI29XnvtNevdu7ebZsOGDXbMMcfQ2AsAAKCQN/Yq0BpZRei//vprWAOvhQsXWqVKlax27dqux4JRo0ZZvXr13KC/FcX36dPHTa+NvOqqq+y2226zypUru/fdfvvt1rRpUzvrrLMKcMsAAACQ3wo0kJ03b5516tQp+Hzw4MHusV+/fvbKK6+4TKzS1zfeeKNt27bNWrVqZdOmTXNpc8/jjz/u0tPKyGpalRrovUqLAwAAoPCKmtKCgkQ/sgAAANGhUPQjCwAAAGSGQBYAAAC+RCALAAAAXyKQBQAAgC8RyAIAAMCXCGQBAADgSwSyAAAA8CUCWQAAAPgSgSwAAAB8iUAWAAAAvkQgCwAAAF8ikAUAAIAvEcgCAADAlwhkAQAA4EsEsgAAAPAlAlkAAAD4EoEsAAAAfIlAFgAAAL5EIAsAAABfIpAFAACALxHIAgAAwJcIZAEAAOBLBLIAAADwJQJZAAAA+BKBLAAAAHyJQBYAAAC+RCALAAAAXyKQBQAAgC8RyAIAAMCXCGQBAADgSwSyAAAA8CUCWQAAAPgSgSwAAAB8iUAWAAAAvkQgCwAAAF8ikAUAAIAvRXUgm5ycbPfcc48dd9xxVrJkSTv++OPtgQcesNTU1OA0gUDARowYYTVr1nTTdOzY0ZYsWVKg6w0AAIAiHsiOGTPGnnvuOZswYYItW7bMxo4da4888og99dRTwWk0bty4cW6auXPnWo0aNaxLly62c+fOAl13AAAAFOFAdvbs2dajRw8755xz7Nhjj7ULLrjAunbtavPmzQtmY8ePH2/Dhg2zXr16WZMmTWzy5Mm2Z88emzJlSkGvPgAAAIpqINuuXTv78ssvbcWKFe75Tz/9ZDNnzrRu3bq556tXr7aNGze64NaTkJBgHTp0sFmzZhXYegMAACD/xVsUGzp0qCUmJlqDBg0sLi7OUlJSbOTIkXbJJZe41xXESvXq1cPep+e///57hvNNSkpyg2fHjh35tg0AAAAoghnZqVOn2muvvebKBBYsWODKBh599FH3GComJibsuUoOIseFGj16tJUvXz441KpVK9+2AQAAAEUwkL3jjjvszjvvtIsvvtiaNm1qffv2tVtvvdUFoqKGXaGZWc+mTZvSZGlD3XXXXS7T6w3r1q3L5y0BAABAkQpk1WgrNjZ8FVVi4HW/pW65FMx+8cUXwdf3799vM2bMsDZt2mQ4X9XRlitXLmwAAACAv0R1jWz37t1dTWzt2rWtcePG9uOPP7qutq688kr3usoHBg0aZKNGjbJ69eq5QX+XKlXK+vTpU9CrDwAAgKIayKq/2HvvvdduvPFGVy6gmx5cd911dt999wWnGTJkiO3du9dNs23bNmvVqpVNmzbNypYtW6DrDgAAgPwVE1DLqCJOvRao0ZfqZSkzAAAA8EdcFtU1sgAAAEBGCGQBAADgSwSyAAAA8CUCWQAAAPgSgSwAAAB8iUAWAAAAvkQgCwAAAF8ikAUAAIAvEcgCAADAlwhkAQAA4EsEsgAAAPAlAlkAAAD4EoEsAAAAfIlAFgAAAL5EIAsAAABfIpAFAACALxHIAgAAwJcIZAEAAOBLBLIAAADwJQJZAAAA+BKBLAAAAHyJQBYAAAC+RCALAAAAXyKQBQAAgC/FH86b1qxZY99++6173LNnj1WtWtWaNWtmrVu3thIlSuT9WgIAAAC5CWSnTJliTz75pP3www9WrVo1O/roo61kyZK2detW++2331wQe+mll9rQoUOtTp06OZk1AAAAkD+B7KmnnmqxsbHWv39/e+utt6x27dphryclJdns2bPtzTfftBYtWtgzzzxjF154Yc7WBgAAAMimmEAgEMjOhP/3f/9n55xzTrZmunnzZlu9erW1bNnS/GDHjh1Wvnx5S0xMtHLlyhX06gAAABRZO3IQl2U7I5vdIFaqVKniBgAAACCqGntFZmq//vprS0lJsbZt29r555+fN2sGAAAA5Ff3W/fee68NGTLEYmJiTBUKt956qw0cODA3swQAAADytkZW5s+fb82bNw8+P/HEE+2nn35yPReI/u7YsaNt27bN/IQaWQAAAP/FZTnKyF577bU2aNAg13esHH/88TZu3Dhbvny5/fzzz/bss8+64BYAAADIbzkKZNV/bI0aNVxXXB999JG9/PLLtmDBAmvTpo2dccYZtn79etfXLAAAABBVpQWeVatW2Q033GClS5e2CRMmWM2aNc3PKC0AAAAo5KUFHpUUfP7559azZ09r3769Pf3004e7rgAAAMBhyVEgq8hYt5/t3r273XPPPdarVy/7/vvvXcnB6aef7upk89off/xhl112mVWuXNlKlSplp5xyimt05lFCecSIES4rrEZnamy2ZMmSPF8PAAAA+DiQ7devn82ZM8fdHEENvFReoABz8uTJNnLkSOvdu7cLdPOKej9Q37TFihWzTz/91JYuXWqPPfaYVahQITjN2LFjXYMzlTjMnTvX1fB26dLFdu7cmWfrAQAAAJ/XyJYtW9Z+/PFHq1u3rrsBgh51K1rP3r177cEHH7RRo0blycrdeeed9t1339m3336b7utadWVi1ZOCF0AnJSVZ9erVbcyYMXbddddlaznUyAIAABTyGtl69erZCy+8YCtWrLDnnnvO6tSpE/a6Lu3nVRArH374obVo0cIuvPBCq1atmjVr1swmTpwYfF1B9MaNG61r167BcQkJCdahQwebNWtWhvNVsKudFDoAAADAX3IUyKq7renTp7uAUt1sqd/Y/KTeEbQMBdBqXHb99dfbzTffbK+++qp7XUGsKAMbSs+919IzevRoF+l7Q61atfJ1OwAAAJD34nMysRpazZs3z46U1NRUl5H1srwKoNWQS8Ht5ZdfHpxOt8iNLDmIHBfqrrvussGDBwefKyNLMAsAAOAvh9X9VnYcRve0aRx11FHWqFGjsHENGza0tWvXur/VsEsis6+bNm1Kk6UNpfID1VyEDgAAACikgawCSJUT7N+/P9PpVq5c6XozUGOr3FKPBeodIZTqc73a3OOOO84Fs1988UXwda3fjBkz3N3GAAAAUHhlu7RANz1QzwADBgxwjat0yV89BpQoUcJ1k6WusWbOnOkeBw4caDfeeGOuV+7WW291AalKC9S1l/qrVWMzDaLyAfVYoNdVR6tBf6u/2T59+uR6+QAAAChEt6hVbwBTp061b775xtasWeO63KpSpYqrXz377LPdzQtC+3nNrY8//tjVtCrTqwysaluvueaa4Ota/fvvv9+ef/55F1C3atXKBd1NmjTJ9jLofgsAACA65CQuy3EgWxgRyAIAABTyfmQBAACAaEEgCwAAAF8ikAUAAEDhvyECAADwP91wKKvuNIH8UqxYMYuLi8uTeRHIAgBQhCiAXb16tQtmgYKiHq50L4DM7sSaHQSyAAAUEeqoaMOGDS4bpluzx8ZSYYgjfwzu2bPH3YXVu4trbhDIAgBQRCQnJ7sgQjc00s2DgIJQsmRJ96hgtlq1arkqM+BUDACAIiIlJcU9Fi9evKBXBUVcqUMnUgcOHMjVfAhkAQAoYnJblwhEyzFIIAsAAABfylGN7JVXXpmtCPull17KzToBAAAcMSNGjLBnn33W1Wy+99571rNnzyPai0SjRo1s8uTJ1rZtWysMLrjgAmvTpo0NHjw435eVo4zstm3bMhw2b95sb775pr3yyiv5t7YAAAB5aNmyZXb//ffb888/73p0+Oc//5kngfEpp5ySrWlfeOEFq1OnTlgQq6RgeoPiLI/W9+STT7bSpUu7rqyaNWtmY8aMCb6+e/duGzp0qB1//PFWokQJq1q1qnXs2NE+/vjjw94urYPm9fvvv4eNV+Dfv3//4PP77rvPRo4caTt27LCoysjqLCU9H3zwgd19992WkJDgVh4AACDaG74pMPvtt9/c8x49ehRI7fBTTz3lAt9IkyZNsn/84x9h4xSwiq58K9v55JNPWocOHSwpKckWLVpkS5cuDU57/fXX2w8//GATJkxwGd8tW7bYrFmz3GNuaB8p1lMGOSMnnXSSHXvssfb666/bDTfcYPkqkAszZ84MtG3bNlCqVKnAkCFDAlu3bg34UWJiYkC7Qo8AABRWe/fuDSxdutQ9+kmHDh0CAwYMcEP58uUDlSpVCgwbNiyQmpoanCYpKSlwxx13BGrWrOniktNOOy3w1VdfBV+fNGmSe+9HH30UaNiwYSAuLi5w+eWXu9//0MHz8ssvBxo0aBBISEgI1K9fP/D000+HrdO6desCF110UaBixYpuec2bNw/MmTPHLSdynhqXnvnz5wdiY2PTxB96z3vvvZfh/ujRo0egf//+me4zbesrr7wSyEtaL+1jrfOiRYvC1qdfv35h044YMSJwxhlnHNaxmJO47LAaey1ZssS6d+/uUtT169e35cuXu3R2xYoV8z7SBgAARZ4ygPHx8fb999+7TOTjjz9uL774YnC/XHHFFfbdd9+5y+/KTl544YUuo7ly5crgNOpDd/To0e59imU0H2U+RWUFGmTixIk2bNgwd3lcpQejRo2ye++9N5iF3LVrl8uE/vnnn/bhhx/aTz/9ZEOGDHF3S7vooovstttus8aNGwfnqXHp+eabb+zEE0+0cuXK5ejz1R2x5syZk+YSf+Q0n3zyie3cuTPDaZS1LVOmTKbD2rVrw96j2td//etfdtddd2W6jqeddprLCCtbHDWlBevWrXPp5Ndee81thA6Uhg0b5t/aAQCAfKdgIzTgKFasmOu0fu/evWH9fKqEUIPqL70+aUV1k+qbVgFe6K1v1Veogs/IWknVdea0E3zdiUzBqy5tK4n2888/u+fXXHONKw944403bP369e5mD3L77bfbZ5995gJVBaKibXnmmWdcbWnk5XoFfp4HH3zQHnvsMevVq5d7ftxxx7nL9qpL7devn02ZMsX+/vtvmzt3rlWqVMlNU7du3eD7FQBqu0PnmZ41a9YE1zfSJZdckmYfKe5Szevw4cPduh177LEuEG7durV169bNNbLy7tam2ttLL73UKleu7La3Xbt27vXQWtwHHnjA7afMpLd+OhlQ+cC3335rZ5xxRrrvO/roo90xtXHjRlcDnF9ylJHVgfPWW2+5Mw2d+egsR2cikQMAAPAPBSbly5cPDjfddJMbr8fQ8ZpOFESFjn/11Vfd+FatWoWN//LLL934Y445Jmz8L7/8kuN1PP3008NqWBW8KQ5RQL1gwQJ361MFdaHZxBkzZgRrYEXBtgKwzChAVeLuqquuCpvXQw89FJzXwoULXeMqL4g9XDpR0ElAehSkazmhg4J577aus2fPdsH8zTff7AJ0BdjKQHsnEu3bt7dVq1a5z+D88893GWgFnQrSPbqrlgLwzAYF5JFUc3v55Ze7xmRZ3b1LWfCoycju27fPPY4dOzbDaXSQhZ6lAQCA6KbLxKFdJSkj6zVEGjduXHC8srHy7rvvpsnIii77R2ZkRZnSyIxsXtIylb2cP39+miymgtDQ4CqrBl3e+qu8QIF5KG/eXpCWW1WqVHHBaHqUzQ3N8qanSZMmbhgwYIDNnDnTBaoK3jt16hT8HDVOw5133umCcWVhFYAqqFdpga6yZ0aZ6Nq1a6cZr54edOLw/vvvp/u+rVu3ukf1lhA1gWzowQkAAAoHr2QgkgK29IK2jALR0KAxVE5rQNOjmtDI5/Xq1XPBpbKjCqzVD2xGl7qzq3r16u6yuLKZujSfHmV1VWerYC29rKyCxOwk9bTe6r9W2eTc9pjQqFEj96iyj8ymSU5OdolJrePhlhaIssMDBw50vVadcMIJaV5fvHixy8QrWI+aQBYAAKAg6HK/ssbXXXedKyVQtlh1rKLMoIJOXe7WOAWI6t9++vTp1rRpU1c/mhPqDkuX7BWAq19Z1XrOmzfP9ZuvdVD9qupu1X+qyi10qf/HH390QZ9KHlS7unr1alcOoGCubNmy6Z4oKHOqwFOX/ZVZDbV9+3ZXXxpK89FJhLq00rI6d+7s5q8GZcq2Kvup5Ysa5Gs9W7Ro4epklVlV0KlleicWKi3QkJtMvjLX2tbIBm2qn+3atavlN25RCwAAop6CVNWUqjW8LqWrfvfaa68Nvq5GXZpG7XjUpufcc891pQ5eXWlOXH311S7jqps8KRBWDwX6W42+RNnMadOmuSBQQbKmefjhh4OlB6pJVb2qgkYFl2qIlh4FmKo3Vn+rkdQWSQFy6KDgXc466yyXkb7wwgtdEK/lqbxD9bCap5x99tmulwUFk2qYr/2lcWrrlFeUjVaZgld66tFz3XtADfHyW8yhfsGKNLWmVPF5YmJinlz+AAAgGinAUPZMAVlGjYyikbKLulPW+PHjrbBRjawC019//dVlXAuDp59+2t0sS8H+4RyLOYnLyMgCAAAUEGVz1YheXXEVFsWKFQtmj/NbjmpkV6xY4VLYAAAAyBvqOqswuTak5COqAlkVT6sLBtWd6J7EursDAABAfvr666/Zwch9acGWLVtc+luPKk5WFxXqMFg3QYgs9AUAAACiJpBVMW737t1dSz519aAWaWqNp0521UpOWdqXX37Z9eMGAACiE+28UViOwcNu7KWOe1VaoO4m1DeZ+krT7dDUPYW6ulCLNQAAED287qH2799f0KuCIm7PoVvXeneRi6rut1R6oLtd6I4bfkD3WwCAokA/+WvXrrUDBw64DvVjY+m8CEf+GFQQq6v3FSpUcP3j5iYuy5c7e6nMwOuQFwAARAddTVXgoP47f//994JeHRRhFSpUsBo1auR6PtyiFgCAIkR3pdIVU8oLUFBUTuCVueQWgSwAAEWMSgr8dGcvICMUxwAAAMCXcp2RVf+xU6dOtd27d1uXLl1808ALAAAARSiQveOOO1xNzRNPPOGe6+/WrVvbkiVLrFSpUjZkyBD74osv3DgAAAAgakoLPv30UzvzzDODz19//XXX6nHlypW2bds2u/DCC+2hhx7Kj/UEAAAADj+QVd9zjRo1Cj6fNm2aXXDBBVanTh3Xpcctt9xiP/74Y77t4tGjR7vlDBo0KKw/shEjRrj+8EqWLGkdO3Z0GWIAAAAUbrE5beUYev+EOXPm2Omnnx7WJ5gys/lh7ty59sILL9hJJ50UNn7s2LE2btw4mzBhgptGfZKpVnfnzp35sh4AAADwYSDboEED++ijj9zfynoqQ9upU6fg6yozqF69ep6v5K5du+zSSy+1iRMnWsWKFYPjFVSPHz/ehg0bZr169bImTZrY5MmT3R0jpkyZkufrAQAAAJ8Gsmrsdeedd7o6WQ3dunWz4447Lvj6J598Yqeddlqer+SAAQPsnHPOsbPOOitsvO5MsnHjRuvatWtwXEJCgnXo0MFmzZqV4fySkpLc7c9CBwAAABTiQPb88893waou7996662u261Q6rngxhtvzNMVfPPNN23BggWuPjaSgliJzALrufdaejQv3cPXG2rVqpWn6wwAAIAo7EdWWdHIzKhn+PDhlpfWrVvnGpCpUVlmdyBRA7BQKjmIHBfqrrvussGDBwefKyNLMAsAAFCIM7Jbt2619evXh41TrewVV1xhvXv3zvO61Pnz59umTZusefPmFh8f74YZM2bYk08+6f72MrGR2Ve9J7NaXZUflCtXLmwAAABAIQ5kVauqHgJCA8YzzjjD9RagutP+/fvbv//97zxbOdXh/vzzz7Zw4cLg0KJFC9fwS38ff/zxrpcC3YTBo5s0KNht06ZNnq0HAAAAfF5aoO62Jk2aFHz+6quvWqVKlVxQqQzpo48+ak8//bT17ds3T1aubNmyrieCUKVLl7bKlSsHx6tP2VGjRrlb42rQ36rV7dOnT56sAwAAAApBIKtL+KG9FEyfPt3OO+88F8TKueeem26jrPyk2+Lu3bvXNTJTH7atWrVyNbUKggEAAFB45SiQVS3p9u3b3Z285IcffrCrrroq+LoaWKnEID99/fXXYc+1TN3ZSwMAAACKjhzVyKqPWDW0Sk1NtXfeecfdPatz587B11esWEHrfwAAAERfRvbBBx90XW+99tprlpycbHfffXfYnbbU56tuRgAAAABEVSB7yimn2LJly9xds9RbgOpRQ1188cXWqFGjvF5HAAAAII2YgO4eUMTphgi6w1diYiJ9ygIAAPgkLstRjWy3bt3cTD0jR450jb88W7ZsISMLAACAIyJHgeznn38e1ivBmDFj3N2+PKqbXb58ed6uIQAAAJDbQDayCoGqBAAAAPgikAUAAAB8Gcjq5gMaIscBAAAAUd39lkoJ+vfvbwkJCe75vn377Prrr7fSpUu75/l9Vy8AAADgsALZfv36hT2/7LLL0kxz+eWX52SWAAAAQP4HspMmTTq8pQAAAAB5jMZeAAAA8CUCWQAAAPgSgSwAAAB8iUAWAAAAvkQgCwAAAF8ikAUAAIAvEcgCAADAlwhkAQAA4EsEsgAAAPAlAlkAAAD4EoEsAAAAfIlAFgAAAL5EIAsAAABfIpAFAACALxHIAgAAwJcIZAEAAOBLBLIAAADwJQJZAAAA+BKBLAAAAHwpvqBXoKhJSQ3Y4rVbbeuufVapTAlrUruSxcXGFPRqAQAA+A6B7BE0c9kGe/bzpbZ5577guCplS9gNZzeydg2POpKrAgAA4HuUFhzBIPbBdxaEBbGi5xqv1wEAAJB9BLJHqJxAmdjMPDdtqZsOAAAAhSCQHT16tLVs2dLKli1r1apVs549e9ry5cvDpgkEAjZixAirWbOmlSxZ0jp27GhLliyxaKKa2MhMbKS/d+xz0wEAAKAQBLIzZsywAQMG2Jw5c+yLL76w5ORk69q1q+3evTs4zdixY23cuHE2YcIEmzt3rtWoUcO6dOliO3futGihhl15OR0AAACivLHXZ599FvZ80qRJLjM7f/58a9++vcvGjh8/3oYNG2a9evVy00yePNmqV69uU6ZMseuuu86igXonyMvpAAAAEOUZ2UiJiYnusVKlSu5x9erVtnHjRpel9SQkJFiHDh1s1qxZFi3UxZZ6J8hMyeJx1qhWxSO2TgAAAH7nm0BW2dfBgwdbu3btrEmTJm6cglhRBjaUnnuvpScpKcl27NgRNuQn9ROrLrYys3d/ij3+0SJLTknN13UBAAAoLHwTyA4cONAWLVpkb7zxRprXYmJi0gS9keMiG5GVL18+ONSqVcvym/qJvfeCU9NkZquWK2HntqhjsTEx9uXPf9j9b82zffuT8319AAAA/C4moKgvyt100032/vvv2zfffGPHHXdccPyqVavshBNOsAULFlizZs2C43v06GEVKlRw9bIZZWQ1eJSRVTCr0oVy5coVyJ29vl/5l418Z4ElJadaw6Mr2AOXtLRyJYvn67oAAABEG8VlSjRmJy6L6oysYmxlYt99912bPn16WBAreq5eCtSjgWf//v2ut4M2bdpkOF/V0WrHhA5HioLWk4+tbJ2aHO0evdvTtqpX3UZf1srKlChmy/7Ybre9Mtv+3rH3iK0XAACA30R1IKuut1577TXXA4H6klXdq4a9ew8GeCofGDRokI0aNcree+89W7x4sfXv399KlSplffr0Mb9pXKuSPdavtSs/WLt5l906aZat/Tt6uhEDAACIJlFdWpBRnau64VLAKlr9+++/355//nnbtm2btWrVyp5++ulgg7C8TmEfCZsS99pdr39v67fstrIli9lDl7S0BkfTowEAACj8duQgLovqQPZIibZAVhL37Ld735hry//cbgnF4lxDsZZ1qxX0agEAAOSrQlMjW5SVL1XcxvRtZc1PqGpJB1Js+NR5Nv3nPwp6tQAAAKIGgWwUK1k83u6/qIV1alLT9XYw5v2F9u73qwt6tQAAAKICgWyUKxYXa0N6nmLntTrYY8Pz05bay1/+4mqDAQAAijICWR/QzRKu69LQruxc3z2fOus3e/zjRZaSyl3AAABA0UUg6xPqweGitnXt1n81NXU9+/nC9fbA2wtc/SwAAEBRRCDrM/9oVtvuvbC5FY+PtTkr/nLddO3ce6CgVwsAAOCII5D1oTb1a9ioS1tZ6YR4W7Jum90+ebZt2bmvoFcLAADgiCKQ9ammtSvZo/1aW6UyCbbm753uLmDrt+wq6NUCAAA4Yghkfez46uXs8f5t7OhKpe2vxL02+JXZtuLP7QW9WgAAAEcEgazP1ahYysb1b231jirv7gY25N9zbMGqzQW9WgAAAPmOQLYQqFA6wcb2Pd1OOa6y7d2fYve+8YN9veTPgl4tAACAfEUgW0iUSoi3By9uae0bHWXJqQF7+N0f7YO5a9xruivYT2u22FeL/3CPeg4AAOB38QW9Asg7xePj7M7zmln5UsXto3m/2zOfLbGff99iy9Zvt80hvRpUKVvCbji7kbVreBS7HwAA+BYZ2UImLjbGBvyjsfXtcKJ7/u2yjWFBrOj5g+8ssJnLNhTQWgIAAOQegWwhpLuAXdKurpUpkXnC/blpSykzAAAAvkUgW0gtXrvVdu1LznSav3fsc9MBAAD4EYFsIbV1V/bu9MUdwQAAgF/R2KuQqlSmRLammzT9F3czhU5NalqNCqXyfb0AAADyCoFsIdWkdiXXO0FkQ69Im3bss1e+Wu4GvefMpkfbGQ2PsrIlix2xdQUAADgcMYFAoMh3Krpjxw4rX768JSYmWrly5aywUK8E6p0gI0N6nOz6nJ3+88H+Zb0DoVhcrJ1Wr5qd1fRoa1G3quvWCwAAINriMgLZHO4wPwazz36+NCwzW7VcCbu+a3g/spsS97q7gX256A9b8/fO4PgyJYpZh8ZHuUxto2Mquh4RAAAA8guBbD7uMD/SnbzUO4EagKl2ViUE6m82PUrQr/prp01f/Ie7E9iWnUnB12pUKGmdmx5tnZscbbWqlDmCWwAAAIqKHWRk82+HFSXerW1VejDzlw22d39K8LUTa5Z3WdqOjWtahdIJuQ6gAQAAhEA2hwhks7bvQIrNXr7RBbXzfttsqYdKq2NjYqzFCVXszKbH2On1q1uJYnEZljRwa1wAAJAVAtkcIpDNme27kw7W0/78h634MzE4vmTxOGvX4CirXqGkvfbNygzff+8Fp4bV5wIAAHgIZHOIQPbwrdu8y2Vpv1z8h/21fW+23qPGZpNv6kyZAQAAyFVcxp29kCtq9NWvU32bPLCTPdavtbWqVy3L93BrXAAAkBe4IQLyhLrlUmOuv3fste9Xbspy+gffmW8n1als9WuWtxNrVrATjypvpUtwEwYAAJB9BLIokFvj7tx7wL77ZaMbPMdUKu16Q3CBbc3yVrdGeUs41Hgsu+gpAQCAooNAFkf81rh6/Y4eJ9vKjYmusdiKP7fbxu17bf3W3W6YvvjPYI8Ix1Yr64La+oeytnoeH5d+RQw9JQAAULRwZy8aex3xW+Om12tB4p79LqD1AtsVGxJt667/3YzBUzw+1k6oXi6YtdXjMZVL26xfNuZ4mXmFLDAAAHmHXgvycYchb2+NmxHdYUzvVWC7/FCAu3LDdtu1LznNtOr2KzklYAdSUo94TwkFmQUmgAYAFEYEsvm4w1BwgZZuwrBh6x4X2HrB7W8bEy0pOeMANpQyuNXLl7JSCXFWsni8lSoebyUT4q1UwqG/3fM497fGuWkS4t1NHtSYLS8yz34PoAmeAQD5jUA2H3cYoktKaqq9PWuVTfpqeb4tQyFsyWCwG+f+Llkszpb9sd32ZxJEVyhV3Eb2Oc3KlizmemRQUKy639wqqAC6qAXPBO0AEP1xGY294GtxsbHW8JiK2Zq2d5sTrEq5ErY3Kdn27E+2PUnJtnd/8qHnKe7vPYde0zg9Tw2Y6Wa8bnxS2rKGzGzfs98GvDgzTUBc2g0KbuNdgHvwuff3ofEJ8VbmUPDrjdPzYnGxLpjMzHPTllrr+jXyNNjLKHhWUKvxhS14JmjP/5OFgjhR4KQIKHxo7EVG1vf043T5k9Mz7SnhcGpkVaersoU9SQdsb1LKwQD3ULA777e/7aN5v2c5D5UxHEjOvH43J7T2Cqyz0qlJTatdpYzrvqx4fJwrj1BDOT13fxeLs4T4OEsodnDcwb8PThNZRpFf+zeaM89FKeNdEMstKsssyOUWpaC9qCyzoJabUgDLLJKlBc8884w98sgjtmHDBmvcuLGNHz/ezjjjjGy9l9IC/zvSwcdPa7bYkH/PyXK6sX1Pt5OPrWz7k1NcALxr3wHbnZRsu/cl2+6kA7Y7zfNDj27c/x7VyE01wkdCQnzswUD3UIDrapO37clW8Ky+gNU9WlxcjMse62/3GBtjcd7fcTFufHzswb+96TSNe9T7Y2Lsppdm2pZ0eq4IDQaev769C+5TAgFLTQ24f3C1vu4xo79Dxh1838ESFY1PTg7Y4/+3yPVznJEKpYvbmMtOd2Umbh8dOlnIzT/sRSloLyrLLOjlFpWgvagss6ht646iFshOnTrV+vbt64LZtm3b2vPPP28vvviiLV261GrXrp3l+wlkC4fc9pSQE0c6S+mywwdSbO6vf9tD/8n4h9HTtkF1K1uyuO0/kOLep8zywceDz/d7zw8NyaqhQK4oIPey3l52OzTbHXmCoNeVHY+Pj7W3vvvNnbRkVm89vHcLKxYf6+qsdUzF6uQg9G89ZvI8sj47J8fwwelT3XvUQ4j+1mPyoXEpKanuGDr4+qHp9Jo3jXtPwPanpNgzny11J3QZUQnNNV0auBMdrbMbYmNMqx8X8rc3PnSa2DTjD17CGPbGD7Zt9/4Ml1mpTIK7xbbel1d0ojR48mzblsnJGFcycoeTov8pTCeARTKQbdWqlZ166qn27LPPBsc1bNjQevbsaaNHj87y/QSyhceRvARSEF/y/AqgFZgkHUgb7O47kGK/rN9mL375S7aC5wqlE1zQolKK5EPDgUOBzsFxBwMd/a1t+d90BwMe7+/DLcXQFocGbun97QVDwUcv6ImNcdnvP7ORfVZJhjK5eVUyciRo34RuryK8fQf8s/6FlT4PndzoREgnKu7x0ElRcFx8XPrj4/73t8ZrXm/M/DXdbgo95UoWs4H/bOKOAf34KwQIBgGuTcDBZ15k4IUIB6f933w0nfdcVzNe/O8vmZ6gaLm3nNM05GQsNt2TrjR/x0aOP/g+LfraZ2dkfvOdciXs5Rs7uhMfrau2JTXN48HtCH9Mf5z+bbpryve2PZOTIl2xeeCilpmeFKX3Stp2wDFhJ0X3vPlDpsutWDrBRvU5Lc1yI+cbk8GLkYvX9uqqY2YngPnVrWWRC2T3799vpUqVsrffftvOO++84PhbbrnFFi5caDNmzEjznqSkJDeE7rBatWrRawGiOgtcUAF0QdTILly92Ya+9n2W0z10SUtXuuH9yKXXTVp+lozoB85lvUMy3Coj2ReR9XZ/HzpBCMuIJ6fY75t22eJ1W7NcZrlSxVwJQ3olFKHPjxQvuFDmVKUk7tGVkIT87ZWLHCot2blnv/2+eVeW8z6+eln3w5ziBRKHts8Nqf8LLLzt9qY5OP3BIOXgtArWVdeekuUyVeaiICmveFlroDAbe+jfwrxWpHot2Lx5s6WkpFj16tXDxuv5xo0b032PsrT333//EVpDFGYKGNVDwJEshNcyFaweqQBa26J6qMyCZy03L7e5aZ3KWd7qWNt76vFV83S52bnFspar6UTBc4ni8VaieP4Hz/ec3zxbPxgZBbiRz5eu22pjP/gpy/npWDupTmUXoCrQc0Gpu7wfk2/ben3Xxnn245jdZY7q0ypPf5Czu9y7ezWzujXKu+y+Bp0Iub+Tvef/+zt0fOg4957kVFu/ZZctXb89y2Xqbog6URDvY/Q+T/d/b5z+iwnJ1ulk8X9/Hpo2xrbs3Ge/bdyR5XJrVirlSke8YzG8Vj3kMfQ1V7/u1bFbvtE/I9oHB0+Ivb/DH3XFaO/+rE+K1OWisuxh0ln3kFz4wecZbJ9OejMrPQptXFwsLmK5IULzloFMl3vwqpiu0mVFv3sFzfeBrCfyH1V9YBn9Q3vXXXfZ4MGD02RkgcOhH/X8OCONpgC6KATPBbXcnAbPWXFlEnEqKM18umrlS9rL05dnudy87Motr7c1WpeZk+Xqu5NX+ze7wfPN3ZoWSNA+6JyTcrVcrxRAQe7CNVvs3jfmZvmeEb1b2El1KoUFpF6Ndei4vNzOey/I3klndmV3uSN6H7xCdSSXqd+egpZ311EKSJUqVSwuLi5N9nXTpk1psrSehIQEl6oOHQC/BtCdmhwdvLyen/SD++rNnd2lpDvPO8U9qpwgv0oovOBZwUDkj39+NjI40sv1gufM5GfQfiSXW1SWWVDL9YLnzORn0J7fy1XAqasCKrFpfnzVbC3ztHrVgjek0dUTr5cR1ztKDq8sFPb9Gw3bejh8XyPrNfZq3ry567XA06hRI+vRoweNvQCfKyr9NRZEvXVBLbeoLLMglluUuhorKsssqOXOpNeCI9/91nPPPWetW7e2F154wSZOnGhLliyxOnXqZPl+ei0AEA2KStBelJZZEMstKkF7UVpmUdvWHUWp1wKPsrFjx451N0Ro0qSJPf7449a+fftsvZdAFgBQmBSVoL0oLbOglpvCnb2iH4EsAACA/+Iy3zf2AgAAQNFEIAsAAABfIpAFAACALxHIAgAAwJcKzZ29csPruEHFxQAAACg4XjyWnY61CGTNbOfOnW5ncJtaAACA6InP1HtBkehHNjdSU1Ptzz//tLJly+bodnWIjrM2nYCsW7eOWw0XMny2hRefbeHE51p47TjCv7UKTRXE1qxZ02JjM6+CJSOrQuHYWDvmmGPy/YNB/tEX60h8uXDk8dkWXny2hROfa+FV7gj+1maVifXQ2AsAAAC+RCALAAAAXyKQha8lJCTY8OHD3SMKFz7bwovPtnDicy28EqL4t5bGXgAAAPAlMrIAAADwJQJZAAAA+BKBLAAAAHyJQBa+9M0331j37t1dZ8m6icX7779f0KuEPDB69Ghr2bKluzlJtWrVrGfPnrZ8+XL2bSHw7LPP2kknnRTsh7J169b26aefFvRqIZ++x/p3edCgQexfnxsxYoT7LEOHGjVqWDQhkIUv7d69204++WSbMGFCQa8K8tCMGTNswIABNmfOHPviiy8sOTnZunbt6j5v+JtuOvPwww/bvHnz3NC5c2fr0aOHLVmypKBXDXlo7ty59sILL7iTFhQOjRs3tg0bNgSHn3/+2aIJd/aCL/3zn/90AwqXzz77LOz5pEmTXGZ2/vz51r59+wJbL+SerqCEGjlypMvS6qRFP5Twv127dtmll15qEydOtIceeqigVwd5JD4+PuqysKHIyAKIWomJie6xUqVKBb0qyEMpKSn25ptvuky7SgxQOOhqyjnnnGNnnXVWQa8K8tDKlStdGd9xxx1nF198sa1atcqiCRlZAFEpEAjY4MGDrV27dtakSZOCXh3kAV2SVOC6b98+K1OmjL333nvWqFEj9m0hoBOTBQsWuNICFB6tWrWyV1991U488UT766+/XKa9TZs2riSocuXKFg0IZAFEpYEDB9qiRYts5syZBb0qyCP169e3hQsX2vbt2+0///mP9evXz9VFE8z627p16+yWW26xadOmWYkSJQp6dZCHQkv4mjZt6k5ETzjhBJs8ebJLNEQDAlkAUeemm26yDz/80PVOoUZCKByKFy9udevWdX+3aNHCZe+eeOIJe/755wt61ZALqmHftGmTNW/ePKx8RN9fNchNSkqyuLg49nEhULp0aRfQqtwgWhDIAoiqcgIFsbrk/PXXX7uaLBTuz1tBDvztzDPPTNOS/YorrrAGDRrY0KFDCWILkaSkJFu2bJmdccYZFi0IZOHb1rG//vpr8Pnq1avdJUs1Cqpdu3aBrhty11hkypQp9sEHH7i+ZDdu3OjGly9f3kqWLMmu9bG7777bXaasVauW7dy509VU6mQlsqcK+I++q5F17MrcqYaS+nZ/u/32212PI/pdVdZdNbI7duxwZUHRgkAWvqR+KDt16hR87tXq6Mv1yiuvFOCaITfUHZN07NgxTTdc/fv3Z+f6mBqK9O3b1/VDqRMT9TOqILZLly4FvWoAMrB+/Xq75JJLbPPmzVa1alU7/fTTXZd5derUsWgRE9C1HQAAAMBn6EcWAAAAvkQgCwAAAF8ikAUAAIAvEcgCAADAlwhkAQAA4EsEsgAAAPAlAlkAAAD4EoEsAAAAfIlAFgB8QHesq1ChQoavr1mzxmJiYtytmnNDd1UbNGhQ8PmePXvs/PPPt3Llyrn5b9++PV+3Q0aMGGGnnHJKptPk1fYC8DcCWQCFhm5jq+BGQ7Fixez444939wrfvXt3Qa+ab02ePNm+/fZbmzVrVvD2srlx0UUX2YoVK3L8ufbs2TNXywVQOMUX9AoAQF76xz/+YZMmTbIDBw64AOzqq692geyzzz6bZlpNo4AXGfvtt9+sYcOG1qRJkzzZTSVLlnQDAOQFMrIACpWEhASrUaOG1apVy/r06WOXXnqpvf/++2GXrF9++WWXrdW0gUDA1q5daz169LAyZcq4S+i9e/e2v/76KzhP733PP/+8m2+pUqXswgsvDLvMnpqaag888IAdc8wxbr6a/rPPPgu+vn//fhs4cKAdddRRVqJECTv22GNt9OjRwdfHjRtnTZs2tdKlS7tl3HjjjbZr164cb/+qVausU6dObh1PPvlkmz17dvC1LVu22CWXXOLWUa9reW+88UamZQaPPfaYffPNNy7LreeRPvroI1cqoO0XXerXtHfccUdwmuuuu84tN6PSgocfftiqV69uZcuWtauuusr27dsXtu+VFf7ggw+C2favv/46W9sLoPAjkAVQqCn7p8yr59dff7W33nrL/vOf/wTrK3XZeuvWrTZjxgz74osvXBZSl8BDee9T4KYAVe8dMGBA8PUnnnjCBX2PPvqoLVq0yM4++2w799xzbeXKle71J5980j788EM3j+XLl9trr73mgllPbGysm2bx4sUucJs+fboNGTIkx9s7bNgwV06h9TvxxBNdAJmcnOxeU4DYvHlz+/jjj91yrr32Wuvbt699//336c7r3XfftWuuucZat27tygr0PFL79u1t586d9uOPP7rn2odVqlRxjx4Fnh06dEh3Gdofw4cPt5EjR9q8efNcoP/MM88EX9e26MRCmXatg4Y2bdpka3sBFAEBACgk+vXrF+jRo0fw+ffffx+oXLlyoHfv3u758OHDA8WKFQts2rQpOM20adMCcXFxgbVr1wbHLVmyJKB/Hn/44Yfg+zTNunXrgtN8+umngdjY2MCGDRvc85o1awZGjhwZtj4tW7YM3Hjjje7vm266KdC5c+dAampqtrblrbfecuvumTRpUqB8+fIZTr969Wq3zi+++GKa7Vi2bFmG7+vWrVvgtttuCz7v0KFD4JZbbgk+198al5lTTz018Oijj7q/e/bs6fZD8eLFAzt27HD7J3QdIrejdevWgeuvvz5sfq1atQqcfPLJGX6uudleAIULGVkAhYqyjSoR0OV7ZRKVMXzqqaeCr9epU8eqVq0afL5s2TJ3KV+Dp1GjRu7yt17z1K5d212S92jeupyu7OqOHTvszz//tLZt24ati55781CDJWUN69evbzfffLNNmzYtbNqvvvrKunTpYkcffbS7xH755Ze7UoCcNlQ76aSTgn8ruymbNm1yjykpKS7zqWkqV67s9pPWQ6UV2aGaY73HG15//XU3XiUHyrqqTEPTqExDNbUzZ85026WygQYNGqQ7T+0f7ctQkc8Pd3sBFH409gJQqKheUg271IirZs2aaRpzqQY1lIIv1V1Gymi8x3stdJrI6UPnceqpp9rq1avt008/tf/+97/ucvlZZ51l77zzjv3+++/WrVs3u/766+3BBx+0SpUquSBQ9aKhZRHZEbq93rK9+lWVPjz++OM2fvz4YD2uutpS/W52tGjRIqy7KwWoXiD70ksv2U8//eRKJHQioFIClRds27Ytw7KCvJDZ9gIo/MjIAihUFJzVrVvXZV6z0yOBgi5lJNetWxcct3TpUktMTHSt9T2aRllXjxoVKWhTXaYaiCloVvAZSl1Whc5D06n2duLEiTZ16lRXp6vaXNWGqq5Tgebpp5/u5hm6rLziZUsvu+wy1zBKDd68Gt7s1htr33qDMsehdbIKkBW0KqDUo7K0mdXHivbPnDlzwsZFPi9evLjLJgNAJDKyAIo0ZUV1eVq9GygQU0CpHgMUfCkD6VGpQr9+/VxjLpUSqDxAWVX1kCBqpa9GSyeccILrsUBdgCl76V1+VyZUl771mgLgt99+271XJQx6j5arEoju3bvbd999Z88991yeb6uCTwXPCrArVqzoekrYuHFjWLB9ONS3rLZLDdjU6M0LbtWzgzLK6fV24LnlllvcftW+bteundtfS5YscUG2R43iPv/8c1fGoZKI3PZlC6DwICMLoEhT9lDdcymwU/ClwFZBlDKmkUFgr169XAlA165dXQ1oaOt6Bba33XabG3TZXj0bqJeCevXquddVUzpmzBgXsLVs2dLdmeqTTz5xQa2CQAWVel3zVTAX2jVXXrn33ntdiYN6VFBwqUA6r240oJIOZU29oFX7U9lu1SNnFigrQ33ffffZ0KFDXY8KKrO44YYbwqZRzwmqLda+0/wU6AOAxKjFF7sCADKmvkwV7HI7VACILmRkAQAA4EsEsgAAAPAlSgsAAADgS2RkAQAA4EsEsgAAAPAlAlkAAAD4EoEsAAAAfIlAFgAAAL5EIAsAAABfIpAFAACALxHIAgAAwJcIZAEAAGB+9P+EOa5Q3q9GrgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Show how ESS changes with proposal width\n", "half_widths = np.linspace(0.3, 5.0, 20)\n", "ess_list = []\n", "\n", "for hw in half_widths:\n", " p = UniformPrior([\n", " (theta_true[0] - hw, theta_true[0] + hw),\n", " (max(0.1, theta_true[1] - hw), theta_true[1] + hw),\n", " ])\n", " prop = p.sample(2000)\n", " _, w = importance_sampling(lambda t: gs.log_prob(t, x_obs) + p.log_prob(t), prop)\n", " ess_list.append(eff_sample_size(w) / 2000 * 100)\n", "\n", "fig, ax = plt.subplots(figsize=(7, 3.5))\n", "ax.plot(half_widths, ess_list, 'o-', c='steelblue')\n", "ax.set(xlabel='Proposal half-width', ylabel='ESS / N (%)',\n", " title='ESS vs proposal width (N=2000 particles)')\n", "ax.axhline(100, ls='--', c='k', lw=0.8, label='perfect (ESS=N)')\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 4 Importance Sampling — NoisyLine\n", "\n", "Two parameters: slope and intercept." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NoisyLine IS ESS=4.2/8000 slope=-0.832 (true -0.8) intercept=5.280 (true 5.0)\n" ] } ], "source": [ "nl = NoisyLine(Nx=40, noise_sigma=0.5)\n", "theta_nl = np.array([-0.8, 5.0])\n", "y_obs = nl.simulate(theta_nl)\n", "\n", "prior_nl = UniformPrior([(-3, 3), (0, 10)])\n", "log_post_nl = lambda t: nl.log_prob(t, y_obs) + prior_nl.log_prob(t)\n", "\n", "prop_nl = prior_nl.sample(8000)\n", "s_nl, w_nl = importance_sampling(log_post_nl, prop_nl)\n", "ess_nl = eff_sample_size(w_nl)\n", "print(f'NoisyLine IS ESS={ess_nl:.1f}/8000 '\n", " f'slope={s_nl[:,0].mean():.3f} (true {theta_nl[0]}) '\n", " f'intercept={s_nl[:,1].mean():.3f} (true {theta_nl[1]})')" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAGZCAYAAAB/r+nWAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZedJREFUeJzt3QmcTfX/x/HP2JeQncleKiEppbSgUNKqUlFREtmSiqSyVIRIpbSjRbQgpQUV5afFVkmlFKJIi+y783+8v/7ndu+dO2NmmmNm7n09H4/DnXPPPfds93zP57smeZ7nGQAAAAAACESeYFYLAAAAAAAIvAEAAAAACBgl3gAAAAAABIjAGwAAAACAABF4AwAAAAAQIAJvAAAAAAACROANAAAAAECACLwBAAAAAAgQgTcAAAAAAAEi8AaQ44wfP96SkpKsUKFCtnr16hTvN2nSxOrUqZOpdeuzmoLc7oULF6a6zKpVq9wyWjan6NChgx122GER8/bs2WNPPfWUnXzyyVaqVCkrUqSIVa1a1S6++GKbOnWq5UTz58+3gQMH2j///BPYcapWrVog687tdE3r2AdJ69f3BO2GG26w8847L8VvNrUpfL89z7NJkybZmWeeaeXKlXP3sEqVKtm5555rzz77bMT3/PXXX9avXz877rjjrGjRolaiRAk79thj7dprr7Wvv/46y/fr7rvvdtub1r3z0UcftdKlS9vevXvd3ytWrHDbU6VKFStcuLAdeeSR1rt3b7ft0X7++Wdr3bq1HX744e5+0rx5c1u8eHHM79ExOuGEE9zxSU5Otl69etnWrVtTLKd5ek/LaFl9Rp+NdtZZZ7nlACAny5fdGwAAqdm1a5d7WHzxxRez7CA98cQT2XrAK1asaJ9++ql7gM3J9LA9ZcoU9zA7aNAgK1iwoHuwfu+99+z999+3Sy+91HJi4K1tVYCsh/+sds8999gtt9yS5euNB7qmFWDmdkuWLLEJEybY559/nuK9Hj16WNu2bVPMD99vBdLDhg2zTp062R133GHFihVzmYcffvihvfnmm3bjjTeGAspTTz3V/a/l6tWrZzt27LAffvjB/e6+/PJLO/7447Nsv7S+hx56yMqXL5/mcm+88YbLXMuXL5/98ccfbhuLFy9u9913nwu+dXwGDBhgH330kS1atMjy5DlQfqNlldlQsmRJe/75512QPHToUJfJuWDBAjvmmGNC3/Hyyy/bNddc447Fww8/7Pa5b9++9u2339rMmTMjtkeBvD7/4IMP2tFHH20TJ060q6++2vbv3x9xLrR9CvRvvvnmiO8CgBzFA4AcZty4cZ5uT+edd56XJ08e78svv4x4v3Hjxl7t2rW9nLrdCxYs8HKT9u3be0WLFg39/fPPP7v9uPfee2Muv2/fPi8nGjFihNvulStXZul6t23b5gVh+/btXk60e/dub8+ePV5OM2DAAHd+g9SmTRvv1FNPjZin60nfq+vrYOezYMGC3nXXXXfQ383zzz/v1vnhhx8edNn/SufyhBNO8Hr27JnmvXP9+vXufvv222+7v5955hm3jbNnz45YbsiQIW7+4sWLQ/PuuOMOL3/+/N6qVatC8zZt2uSVKVPGHVPf3r17vYoVK3otWrSIWOfLL7/s1vnOO++E5s2YMcPNmzhxYsSyzZs395KTk926wtWpU8fr1KlTBo8OABw6VDUHkGP16dPHVXtUacjB7Ny505U2Va9e3QoUKGBHHHGEdevWLUW141hVzceOHetKnFQ9UiVUqu551113haqZqvRHpTfRPv74Y1d187XXXkv3PsWqau5XoV22bJkrzVGVU5VMqcrrpk2bIj6vqqwqtVeVS1X9VAnT5Zdf7kqjs4pfjVSl87H4pVxp0f50797dVVdXSZVKzFWlNlY10W+++caVsmlf/OqkKnUMpxKu+++/35Vmab9Voq0SwUceeSR0DFVyKLoG/GrAc+bMCa1j8uTJdtppp7lqvTrXqv6rErxY1e6XLl1qLVq0cNfDOeeck2pV8/Red/rcBRdc4Eoz69ev7/ZTpfOp8ZtTqCS5UaNGbp+1jnHjxrn3Z8yYYSeeeKJrAlC3bl1XEyGcqghff/31VrNmTbeMtuvCCy90+xVOx0fHSbVKbrvtNreczpU+L88880zE+VOJY6zjEF3l2m92oZJRlUKWKVPG/ZZVgvnbb79FfFbnRcda15v2s1atWnbnnXfatm3b7GBUkqxjpXXrsyqVveyyy2z79u2WUb///rtrRqHaHpmh7VUtnfT8brLiN5ZeKi3++++/7YEHHkhzOe27rv1mzZq5v/Pnz+/+1/0onF+bRNdw+GfPPvts1xzFp5Jyne+33norVHX9s88+s3Xr1rlrM9wVV1zhvju8GYu/PXovnD6rayi6VoLOm67PLVu2pPPIAMChReANIMdS0KOq5qrarAfs1CgYveSSS1xVSj18KShRO0QFb3oY1MNwahQIdu3a1Ro3buwe9KZNm2a33npr6KFfAcZFF11kTz75pO3bty/is2PGjHFtD7Oq2rUCBgU5qu6pwEMPkdqWcJ07d3bVv/VwrG1VEK6AXcGZAofogCoz7W4V+OjhWoHh008/7TILMmP69OmuzejgwYPt9ddfdw/lyljQa9/y5cvdtmsftKwCUwV4Cu6GDx8eWk6vtS/6vM6vgrWOHTuGAlxVW1VVYNE6FLBqUnAqQ4YMcZ/Vul999VUXaOoBXdVjVcU13O7du90517Wj6sGpBcgZve7U3lWZAz179nSBss53WtavX++CDO2btkMBtjJjdDwV7CtjSteKghNtR3hAq9cKRhV06bsef/xxl4HUsGFDd8yjaX2//PKLu84VKKl9ss79TTfd5DI4dEz1W9SxCM/MOBhtuwI4Xcs6h/qsqhmH+/HHH+3888+35557zm2rrm+dI2UUpEXXZatWrVyGh6o367PaX2Ws6BxmlKo5q2+Dpk2bxnxfmT8KIKMnnzIXjjrqKPebHDVqlH3//ffuGolFGUBy3XXXud9xrDbTPq0j1vemtS0+XdvKsFLmYnQ/DtF0LSlzSJksomtKGRnKkNHvU9XildmoY6xzo/uEqIr8Tz/9FLNqvObpfT9jUJls/vxwukaU4em/7y+r79B1G73O8HX5lAGj+3ZGrk8AOKQOYek6AGS4yvauXbu8GjVqeA0aNPD279/v3o+uLvnee++55YcPHx6xnsmTJ7v5Tz/9dGiePqvJ1717d+/www9Pc3s++ugjt56pU6eG5v36669evnz5vEGDBsXc7tT41Va1bHQV2ujt79q1q1eoUKHQfn/66aduuZEjR0Yst2bNGq9w4cJenz59QvPmzJnj5c2bN2L70lvV3K/mqWqi+j5NpUuX9q644gpv+vTpXnroM9omVV/1qWroscce6x111FGheVdddZWrnvvLL79EfL5ly5ZekSJFvH/++cf9fcEFF7jqspmpaq5161z16NEjYv6WLVu8ChUqRFSF1bHQOlQVOJreq1q1aqauO31O52P58uVeeuga1ToWLlwYmvfXX3+5dei46vrzqSmGln300UdTXZ+OvaqQ16xZ07v11ltTXNtnnXVWiqrOOjYNGzaMmL969WpXpTj8OIjWoes4+regazicjpXmr1u3LuZ26lpX1ei5c+e65b766qtUq5q//vrr7u/opiiZdfPNN7tj6//eon+zqU2ffPJJaNkvvvjCq1KlSui9YsWKuWv3hRdeSLHewYMHewUKFAgtW716da9Lly4R+xx+jtIzhV/7Ooc6f1dffXVoXmpVzf/880/3G3njjTci5v/222/eaaedFvEdug/s3LkztIyuRc0fOnRoivWqmrjemz9/vvv7gQceSPX8q/r50UcfHfpb1+q5556bYjltk9ahKu/hdH0nJSV5ffv2TfEZAMgJKPEGkKOpNEslNuopXKVgsfil4SolDacqiir9+uCDD1Jd/ymnnOJKTVUaqlLFP//8M8UyKklRVXSVGvpUMqgSZZUIZhWVskaX7Kgq84YNG9zfb7/9tvtOlRiGl3JVqFDBbV94SY9K8PXevffem6ltUQmkSkBVC+D222+32rVru5I5baOqkKeHqmiHd+aUN29eu/LKK1015rVr14bOnZarXLlyxGd1LlVdWKXW/nn66quvXO0E1YDYvHlzuvdFy+tYqHQx/LipqqyOU6wSsoOVRvvb7m9req47nU/VaEgvVUM+6aSTQn+rd3mVRKsqvmpa+PySx/ARALR/KuVXCb9+Qyo11P8qXf7uu+8Our8qFVeJe5s2bSLmqwT09NNP/0/XdPS2qjRUHWXpOtY1otJPnReJta0+HQftk36DqmXwX5tbqJZA2bJlU+05XR3rqaOv6Enb4dMoALq+Vfqu5ioq2dZ1oGtPxyK8BFyd9ek3ptJ61WRRibTuKzrnr7zySmg5/R3re2NN4deFSt11vkePHn3Qfde9T8cyvDf3jRs3uiYg+q2pQzSVdqs0f968eW5fokvY0+pxPvq91JZN73Kx3tN1o5o6v/76a6qfAYDsRK/mAHK8q666ylXn7d+/v2szGE3VNBVY6KE5+sFMD/NpVeNUFWE9QKotq4IPVSfVw7OCffWS61P1YFWbVUBSo0YNt7zaVmv9WUVVg8P5VT5VVVNUlVwP7qn1TKztykpqM6vqpppEQULLli1dBoTa7SoYT0usY+PP0zlRb9D6P1Y7Vz+A8M+dqkIrmH3ppZdccKIATUMIqQfpBg0apLkdfhV8ndf0tKdVm2i1Tz2YjF53qbXnTY0C7WgKjqLna54ok8anKu86T+ofQUGs2s9rP3UN+9dTWtvmb3usa03zVq5cmSXXtKovq7q/MkH0m1PGhI7/mjVr3G891rb6NDLA7NmzXRV2tatXNWP9BvRbzUzv8/qu8HbL0XS9Huxa8wNA9R+gyT+Wulco4+zdd991mVrhx1LNCfw2zwpu9RvT9iszUBSQhwf3afGrZeu3qkw3VQvX9eE3ydC9Tvc4/a1zod+4qPmHvlfH3qfflnpDVyaJf33oXKlKuJpSKBhv3769u7Z0zce6z6ptufjXrH89aNnoa0vLhl/bWjY96wyn85fWNQMA2YkSbwA5nh7q9BCodoRqdxrNH3dWQ9qEU5CqUju1vUyLHno1FJU6MlM7XX1ObR3DS+VUIqfvUTCjztS0Xj3sH0raDx0LlTjFKu1SiXSQVNrpl/CrzefB6BilNs9/ANf/6mwpmt9e2T93CigUTKqdtB68VSKo4EzBzcE60vLXoeAi1nGL7qQpvWNFZ/S6OxRjUPuUQaFSVpV66xipxoCCxlg1OmJtm39+wvsNSOu8ZpZqDehcq9RXmQLKTNF2qn+H9FAgqDbp+u2q4y6VMKuNeKxO/A5G58sP6rKSjqU/xnR0u+Ro2n91NKdryq/pMnfuXBfMp2fy+2NQ6b8CUAXwCoz96X//+5+rRaDXyswSHTuVykfXelDQrc72ojNl/Awsf18UvKtte3THfaJ5et/PFFQ/Bf78cPodqU18+BjjWlbbGl2y7n821njkKqU/2P0eALILgTeAXEGdiakEWh1LqZQsnN/rtIKN6M6CVArmv38wKlFVqY9K1tU5U3hwqZIUv0qrqnCqBCojVW6zgjIDFNSpKqWCk+jJf6j9r9TpWPQx9vlVf8OrtKZGD/PhgZs6p1OnaCqp9Mc+1rnxg69wL7zwgit90zjC0VSdVCWIyvhQoOQHG9GlqT4FngrclXET67ilpxQzlqy67oKgQNo/Hj5lKqW3Gq56j1epfXTzDpWkKpMqK7dTordVveFnhGpAqOM4vzmIMmgySiW5KmGNHkkgvdQxW2q1a6J/N/pdqOQ5mn4jqh6ua9/vPTwzVc11f1KP8tGTmqSow0i99puMKONC50H3l3Bal5qERF8zfvOP8PHL1cGkfsfKDAu/j6hTPlVL90vidY4UyIeP6uBniumeE16jSevUPP2ewukerG3TusLpHqJaH2peAQA5EVXNAeQaKvXWQ6hKgsKrOSsgV3ClarVqj6iA+Ouvv7YBAwa4oZvSGh6oU6dOrkRGn9EDoUrzNHSYhtCJrpqs9sWq1rpo0SJ79tlnU12nHkBj9QQeXsU0M7SNCv5VQq827yodU2aBSoxVCq7AW1XA/VIyBX6qbprRdt6qTq/jqSr+qqas46KSJAVuqnGgNu/qifxgVPKkKqlqy6rtVPtQlWqFl0bqHKkKrnqS1naq+qiqsOq7dKz9oYzUi7JKuBQkq2q3aiOo7ap6SteQWeJnPGiIMVWBVQmgAkgFGsqwUYaKSgLVjlUlfgp+vvjiC7dtaQ3tlZr/ct0FTUGUghsFk2pXrWt2xIgREcFSWlQtXcdEbY+VyaHe1FU9WfN0PWTVcFe6jnQuunTp4o6bzpnOv9rzH4yaHOi3pp7NVRtDQZdKzsUfEisjdF0rY0s1IFTqHE2ZDipVj6brUZlJCth1ramNv75f/RYocFQfArom1RbfDyzVq74yF1STRvcZXecKcnVfUYaffgt+EwKV/mc0c0hBe/Swif58lSCHv6egV9dydC0DZWzpXOg9jbKg/VEpt5oEqJp4u3btQsuqHwjtk86FfmvKSFE1d52T8JEVlEGi37V+G7q2VJ1eGQ3qoV/fE97GXJmgmqd7mn5fKlVXTRe1n1dml9YVzj83qfVKDwDZLrt7dwOAaGn1Dt62bVv3XnTPvDt27HC92aq3ZfW6XLFiRddL8caNGyOWi+7VfMKECV7Tpk298uXLux6Gk5OTXS/XX3/9dcwT06RJE69UqVLe9u3bU93utHocTqtX8z/++CPm+qJ76VaP2+qtWD2RqxfmI4880rvuuusiesD2e0IO72k6vb2a65jdf//93tlnn+0dccQR7rjoffUqrvmx9j2avrtbt27eE0884bZP50Q9mr/88sspll26dKl34YUXeiVKlHDfVa9evYjjI+rJvVGjRq6ndS2jnqM7duzorVq1KmK5fv36uXOYJ08etw06Dr5p06a5c128eHHXk7qulcsvv9ybPXt2qsci+jhF9+ad3utO77dq1cpLr9R6n05tPf7x9un7dXzKlSvneoc/44wzXO/b0de/f5289tprMbdDPbOrF3odc/U4rWvv4osv9urXr5+uXs2jf8P+94WfF/V4rZ6ztZ1ly5b1brzxRm/x4sWp/k586uX/0ksvdcdE51M972vf0tvzfjT1Al6tWrUUPbEfrFfzdu3aueU0AsNDDz3keuTX9alt0qgEtWrVciMOqFd637fffuvddtttbrQG7bN6FC9ZsqTb/hdffNELSvR1tXXrVreN0b83n86DjnGlSpXc/miECZ2f6FEIZMWKFd4ll1zifl86l+ecc463aNGimOtVb+fHH3+8u67Ue37Pnj3dKAPRNE/vaRktq8+88sorMdd57bXXenXr1s3A0QCAQytJ/2R38A8AuYFK2lXCqvGiw8eYRkqquqoSM411jvihUm91gKYO92L1t5DbjRw50h544AFXvdrveCyeqSmBSq5V+yNWZ2W5hUrEVf384YcfdrWYACAnoo03AByEqoCqt+GOHTu6KraZ6TEZyG3U7EKZTGqnq6YLanevarxquxuvvwFlFqnad/jQgfFMw8WpbXpuDrpFAbeaG/i9wwNATkTgDQAHoXaXahOptpdq86iefoF4p3a66qtAfRuora2G6VLbXrVZPthQcrmVOlFUW+Xozt6Qs2n4P/Vp4HfiBgA5EVXNAQAAAAAIECXeAAAAAAAEiMAbAAAAAIAAEXgDAAAAABAgAm8AAAAAAAJE4A0AAAAAQIAIvAEAAAAACBCBNwAAAAAAASLwBgAAAAAgQATeAAAAAAAEiMAbAAAAAIAAEXgDAAAAABAgAm8AAAAAAAJE4A0AAAAAQIAIvAEAAAAACBCBNwAAAAAAASLwBgAAAAAgQATeAAAAAAAEiMAbAAAAAIAAEXgDAAAAABAgAm8AAAAAAAJE4A0AAAAAQIAIvAEAAAAACBCBNwAAAAAAASLwBgAAAAAgQATeAAAAAAAEiMAbAAAAAIAAEXgDAAAAABAgAm8AAAAAAAJE4A0AAAAAQIAIvAEAAAAACBCBNwAAAAAAASLwBgAAAAAgQATeAAAAAAAEiMAbcWf8+PGWlJRkhQoVstWrV6d4v0mTJlanTp1MrVuf1RTkdi9cuDDVZVatWuWW0bKHWmrb9/7771uLFi0sOTnZChYs6P7XMXrwwQctJ6tWrZp16NAhuzcDAHI90t1gj2tuSHeVnipdzYz58+fbwIED7Z9//rF4sX37drdPc+bMye5NQQ5C4I24tWvXLrv77ruzdJ1PPPGEm7JLxYoV7dNPP7VWrVpZTvDkk0/aeeedZ8WLF7cxY8a4h4Fhw4ZZrVq17PXXX8/uzQMAHEKku4mb7t5zzz02derUTAfegwYNirvAW/tE4I1w+SL+AuKIEqaJEyfa7bffbvXq1cuSdR533HGWnZSzfeqpp1pOMXToUDvrrLNSJPbXXnut7d+/P9u2CwBw6JHuJm66e+SRR1pODH6LFCmS3ZsBhFDijbjVp08fK126tPXt2/egy+7cudP69etn1atXtwIFCtgRRxxh3bp1S5H7Gquq+dixY11gf9hhh1mxYsXs2GOPtbvuuitUNTxfvnwuoYz28ccfuypkr732Wrr3KVZVc1Vl0rxly5bZ1VdfbSVKlLDy5cvbDTfcYJs2bYr4vOd5rsT+hBNOsMKFC1vJkiXt8ssvt59//tky46+//nKl8LHkyZN9t5clS5bYBRdcYOXKlQtVw1MtgbVr16b5uV9++cWuueaa0OdUgjBy5MiIhxn/HAwfPtweeOABq1KlimvW0KBBA/vggw9SrPPHH3+0tm3bRqzz8ccfD2S/ASA7ke4mbrobq6q50sru3bvbiy++6NI+BcF6Xnr77bcjnmHuuOMO91rPYPqMpvCS4smTJ9tpp51mRYsWdc9a5557rkvno79f7y1dutRVw9fz2DnnnBOqiTF48GC3DUqv9WzYtGlTV9Ke0fPkN1f85JNPXEGIltUzo0r89+3bF3pOKFu2rHutUm9/n2jeBgJvxC3ddFXVXNWwPvzww1SX0832kksusYceesjlGM+YMcN69+5tEyZMsLPPPtvdsFMzadIk69q1qzVu3NhVsZo2bZrdeuuttm3bNve+EqGLLrrIVQ3zb8g+VRFTQHjppZdmyf5edtlldvTRR9sbb7xhd955pyvt17aE69y5s/Xq1cuaNWvmtlWJjAL2Ro0a2e+//57h71RCqO9TwvnVV1+l2MfsoGPfvHlztz8KcGfNmmWjR492AfKWLVtS/dwff/zhjsPMmTPtvvvus+nTp7vjpBoTenCIpvP33nvvuXW/9NJL7oGnZcuWrimA79tvv7WTTz7ZvvnmGxfA62FDGQA9e/Z0iTEAxBPS3cRMd9OiZyqllwp8td2lSpVyzz1+QHvjjTdajx493OspU6a4NFTTiSee6OYNGTLEFSqoxuGrr77qgnil5WeeeaZLY8Pt3r3bPXPp2e3NN9906ezevXtd2qx0XRnyelZT4YWOvzLbM3Oe1q9fb1dddZW1a9fOfY8C9Pvvv99uueUW974yRvR8IB07dgztk4JzJDgPiDPjxo3zdGkvWLDA27Vrl1ejRg2vQYMG3v79+937jRs39mrXrh1a/r333nPLDx8+PGI9kydPdvOffvrp0Dx9VpOve/fu3uGHH57m9nz00UduPVOnTg3N+/XXX718+fJ5gwYNirndqVm5cqVbRsv6BgwYEHP7u3bt6hUqVCi0359++qlbbuTIkRHLrVmzxitcuLDXp0+fNPcj1vatWLHCq1OnjpuvSes555xzvDFjxni7d+/2ssPChQvdtkybNi3N5apWreq1b98+9Pedd97pPvf5559HLHfzzTd7SUlJ3vLlyyPOQXJysrdjx47Qcps3b/ZKlSrlNWvWLDTv3HPP9SpVquRt2rQpYp26bnRu/v777/+8vwCQ3Uh3EzvdFaWnSlfDafvKly/v0kff+vXrvTx58nhDhw4NzRsxYoRbVulruF9++cU9K/Xo0SNi/pYtW7wKFSp4bdq0ifh+reP555+PWPaFF15w85955plUtz0j50nPgFr2zTffjFi2U6dObr9Wr17t/v7jjz/ccnpGA3yUeCOuqdq4ciHVI6hySmPxS8OjqwBdccUVrlpTrOrDvlNOOcVVR1durHI9//zzzxTLqFqSqlaFVy9WCbiqHd10002WVZTLG+744493Veg3bNjg/lZpq75TVamVA+xPFSpUcNuXmQ5A1KZLOe5z5851OcvKKV6wYIErIVauvL4/NUqTw7cjI1Na7diOOuooV0VMTQx0nKNzxFOj60A56jqn4XRdaFuja020bt3aVVkLL+m58MILXRMClUBo33XtKGdf1evCt//8889373/22Wfp2jYAyC1IdxMv3U2LqnQrffSpKZyaXsUadSaaaizqu6+77rqIbVHaq5qGsY6fav+Fe/fdd93yan6XmoyeJ+1P9DOXmpTpGOkZAEgNgTfinqoDqcpS//79bc+ePTHbS6kdtt8ex6ebsG66ej81qpr+/PPPuwREN3slJg0bNnTVm8OparGCsOXLl7tteOaZZ1zVJK0/q6jNUji1J5YdO3a4/1VVSomuEr38+fNHTAoAY2UapIeqWKujl3vvvddVz/7tt9/syiuvtEWLFrljkxpV5Y/ejvROaSWgauOuBxK101Jb+9q1a7sq/QMGDIh5/g/Wbk6f9d8PF+vcaZ6qum3dutUtr4T7scceS7H9Crwls8ccAHIy0t3ESncz8mziP5/4zyZp8at4q8lW9Pao3Xf08VMmt3p7j25GpnQ8rfbvGT1PWi6a/0yQ1jMjQK/miHsKoDXUhtr9Pv300zETBQVIujmHB9+6Casdj274abn++uvdpLbFyulUgKd2RD/88INVrVo1lBOqEliVeqszDq1XnbcdSmXKlHHHQh2C+EF5uFjzMkO1BNRRnRJFtW1OjUqHlUuf2X1JS926dV37e53Dr7/+2rXnUvsydYKi9u+x6DpYt25divl6oIn1nTqH0TRPpT3q4EUJdt68eV3mTGrnWh3JAEC8Id1NvHQ3CP53qgd3/3kqLTrW0fRcN2/ePFcanVrwndHzFKttvv9MECujAfAReCMhqCqWAm8FX5UrV454T71eqodqdZAV3imKOgFRMO33ipmehE8deKjEU521qVMOP6FQNSdVK1cHI+pFU6Wxp59+uh1Kygx48MEH7ddff7U2bdpkyToVqMYqJf7uu+8iSotjUeIUdAKlhFTVxB5++GEXfC9evDjVZXWe1fu8lvE7dZEXXnjBrUfV5cKpE5gRI0aEqpurs5e33nrLdfiigFs57/qMel5VtX8F5ACQKEh3EzPdzYzoGno+9V6uGok//fRTiirk6aXnsldeecU9A6RWap/R86T0XjUNwqubq0NbvyZCWvuExEbgjYShUu+TTjrJtb1S9WOfAnLd3FUivXnzZhcQq5RUJdf169d3JZap6dSpkytF1WeUECrHU8GbqjtHl5Sr93MF+KoK9uyzz6a6TrUl1lAU0fzqyZmlbVTwr9J5tXlX4qDMAiXiyg1WKfHNN9+coXXqOCpgVcKmdmdqW/b555+7HrxVFUu9eR5qaqul3kiV+VGjRg1X6q0gWW3xda5To0wXBdnqdVwZNMo0UW+sWpeOi3qMD6fgWutTD/jKSdf1pesnvLfyRx55xM444wwXjGsd6uVeCfaKFStckJ5Wb/sAkNuR7iZGuvtf6Tj4aWb79u1djbFjjjnGpZlKj9VUUL2ga5x49eGiEucvvvjCHcuDjRCiPnjGjRtnXbp0cc39lCGuNFvHTMOLqVlERp+PlHmhv9Urup4N3nnnHdeEUPM0gorfDlzPEer/R+dLvbmrZD16yDUkmFA3a0CcSKt38LZt27r3wns1F/VO3bdvX9cjZ/78+b2KFSu63qw3btwYsVx0r+YTJkzwmjZt6nrtLFCggOvpWr1sfv311zG3rUmTJq7n6+3bt6e63alN6u0zrV7N1YNmrPVF9xKqHj8bNmzoFS1a1PXWeeSRR3rXXXed6w08o8f1qaee8lq3bu16ji9SpIg7Blpfly5dXG+g2eH777/3rr76arcd2r8SJUp4p5xyijd+/Pg0ezUX9Uaqa6R06dLuOjjmmGNcb6v79u0LLeOfg2HDhrle6dVrufa7fv363vvvv59ie7T8DTfc4B1xxBFunWXLlvUaNWrk3X///QEeBQA4dEh3EzvdTatX827duqVYNlb6269fP/cMpZ7B9TmNCOPTKCV61ipevLhXsGBB9/nLL7/cmz17dsT36/jGome8e++916tZs6Y7Xkrjzz77bG/+/PkZPk/+yDhz5sxxI+Zoe/TMeNddd3l79uyJWJ+2T88GWkb7FL3PSDxJ+ie7g38gEaikXbmfGq9SJd/InVQbQW2zVc1cY3wDAHIm0l1kNY1Uo87W0mpLD6SGquZAwNauXeuqSClQU/ufW265hWMOAADpLoAEwnBiQMDUnls5pOps7eWXX7YjjjiCYw4AAOkugARCVXMAAAAAAAJEiTcAAAAAAAEi8AYAAAAAIEAE3gAAAAAABChX9mquge9/++03Nzh9UlJSdm8OAABZRqN8btmyxZKTk91ICPGAdBsAkOhpdq4MvBV0V65cObs3AwCAwKxZs8YqVaoUF0eYdBsAkOhpdq4MvFXS7e9g8eLFs3tzAADIMps3b3aZy35aFw9ItwEAiZ5m58rA269erqCbwBsAEI/iqSkV6TYAINHT7PhoPAYAAAAAQA5F4A0AQdm506xVqwOTXgMAuF8CSEi5sqo5AOQK+/aZvfPOv68BANwvASQkAm8AyCb79u2zPXv2cPwTTP78+S1v3rzZvRkAgAwOi7h7926OWYLJn4VpNoE3AGTDmI/r16+3f/75h2OfoA4//HCrUKFCXHWgBgDxSgH3ypUrXfCNxHN4FqXZBN4AcIj5QXe5cuWsSJEiBF8Jlumyfft227Bhg/u7YsWK2b1JAICD3LfXrVvnSj01bFSePHSRlSi8LE6zCbwB4BBXL/eD7tKlS3PsE1DhwoXd/0rIdR1Q7RwAcq69e/e64Cs5OdllliOxFM7CNJssGwA4hPw23STeic0//7TxB4Ccn2EuBQoUyO5NQS5Pswm8ASAb0LY3sXH+ASB34b6duLLq3FPVHACCUrSoGghxfAGA+yWABEfgDSBd+k1Zmq7lhrauyxFF3Bg/frz16tWLHugBxCXSdsST8Tk8zSbwBoBc9gCUVTKaSdKhQweXmE2bNi3U0cg999xj7777rv3+++9WsmRJq1evng0cONBOO+00yw6rVq2y6tWr25IlS+yEE074z+u78sor7fzzz8+SbQMAxA/S7P9uVYKl2QTeABCUnTvNrr32wOsXXzQrVCiujvVll13mOhqZMGGC1ahRwwXfH3zwgf39998WD7Rv6s3U79H0v6wnf/78WbZdQFyK8/slkN1Is7M/zaZzNQAIinpCff31A9P/94oaL1TyPW/ePBs2bJg1bdrUqlataqeccor169fPWrVqlWap+SWXXGKDBg1yw3IUL17cOnfubLt37w4ts2vXLuvZs6d7v1ChQnbGGWfYggULQu9v3LjR2rVrZ2XLlnVBcc2aNW3cuHHuPeWcS/369V1nKE2aNAl9TsvUqlXLrfPYY4+1J554IiLXXcu/+uqr7jNa5qWXXnLV1g4//PCIfRg7dqwdeeSRrofbY445xl5UkBBG63nyySft4osvtqJFi9r999//n441kBDi+H4JZDfS7CNzRJpN4A0AyLDDDjvMTap2rkA5I1Qq/t1339lHH31kr7zyik2dOtUF4r4+ffrYG2+84UrSFy9ebEcddZSde+65oZJ0VW//9ttvXRV3rUeBcJkyZdx7X3zxhft/9uzZtm7dOpsyZYr7+5lnnrH+/fvbAw884D4zZMgQtx59R7i+ffu6oF/L6DujaVtvueUWu+222+ybb75xmQbXX3+925dwAwYMcIn40qVL7YYbbsjQ8QEAICuRZt+WI9JsqpoDADKeeOTL50qDO3Xq5HKKTzzxRGvcuLFdddVVdvzxx6f5WeU6P//8825czNq1a9vgwYPtjjvusPvuu8927NjhAmmtu2XLlqGgedasWfbcc8+55X755RdXot2gQQP3frVq1ULrVim4lC5d2ipUqBCar3WPHDnSWrduHSoZV/D+1FNPWfv27UPLqVMWf5lYHnroIVdq37VrV/d379697bPPPnPzVfLva9u2LQE3ACBHIM3umiPSbEq8AQCZbi/222+/2fTp013p8Jw5c1wArqA5LeqATUG3Tx2xbd261dasWWM//fSTa191+umnh95XWytVY1cptNx88802adIk1xGLSsfnz5+f5vf98ccfbt0dO3YM5fprUnUyfV84P5hPjbYhfNtEf/vblt71AABwKJFmZ3+aTeANAMg0tYVu3ry53XvvvS4AVmmwqmxlhtpZef8/7rleh9N8f55KwlevXu1KpxX4n3POOXb77benut79+/eHSs6//PLL0KRqZ8r5Dqf2XenZztS2LSPrAQDgUCLNtmxNswm8AQBZ5rjjjrNt27alucxXX33lqpT7FPyqBLpSpUquPbeqoqvjNp9KwBcuXOg6RguvUq4gXx2gjR492p5++mk3X5+VfWGdM5UvX96OOOII+/nnn936wye/M7b00jaEb5sowyF82wAAyA1Isw8t2ngDADLsr7/+siuuuMK1iVKb7mLFirngePjw4a6DkrSoB3NV+7777rtdybVKyLt372558uRxuc6qSq623KVKlbIqVaq4dW7fvt19RlS6ftJJJ7n24erY7e233w4FvuoJXT2dv/feey6QV+5+iRIl3Nji6jRNvairxFyf0/aqh3S1+UovbVebNm1clXqVtL/11luuAzd15gYAQE5Emn1ijkizCbwBIChqx7x167+v44hKqBs2bGgPP/xwqF125cqVXWdrd911V5qfVeKnIcDOOussFwCrQzYFxr4HH3zQVQ+/9tprbcuWLa7t1fvvv28lS5YMlWpr2DINAaYg+8wzz3Rtvv0OZB599FHXYZsCdL2ntuc33nija1c+YsQI1y5cAX7dunVddfWM0FBojzzyiFuPAnmVmGuYsvBhywBkQhzfL4HsRpo9Ikek2Ume36AuF9m8ebMrwdi0aZMrvQAQvH5TlqZruaGt6wa+LbnZzp07beXKle7mr9LYRKPq4RpPVMOQJbK0roN4TOPicZ+AeEDafnCJnG6TZmdtmk0bbwAAAAAAAkTgDQBB2bVL2cUHJr0GAHC/BJCQaOMNAEHZu9dswoQDrx9/3KxgwYQ/1gcb4xtAguJ+CeQ4pNlZixJvAAAS1NChQ+3kk092vdKrR3h1Hrd8+fKIZdQVjDq/S05Odp3ZqVOaZcuWZds2AwCQGxF4AwCQoObOnWvdunVzY6nPmjXL9u7day1atIgYi13DuY0aNcrGjBljCxYssAoVKljz5s1dj/MAACB9qGoOAECC0njn4TTMikq+Fy1a5IZ7U2n36NGjrX///ta6dWu3zIQJE6x8+fI2ceJE69y5czZtOQAAcV7i/fHHH9uFF17oqpwlJSVFDAmjcVz79u3rxkbVGKla5rrrrrPffvstYh0at7VHjx5WpkwZt9xFF11ka9euzZo9AgAAmaLhUKRUqVLufw2fsn79elcK7itYsKA1btzY5s+fn+p6lM5riJXwCQCARJbhwFvVz+rVq+eqnEXbvn27LV682O655x73/5QpU+yHH35wgXW4Xr162dSpU23SpEk2b94827p1q11wwQW2b9++/7Y3AAAgU1S63bt3bzvjjDOsTp06bp6CblEJdzj97b+XWttxjWvqT5UrV+asAAASWoarmrds2dJNsShxVRuxcI899pidcsop9ssvv1iVKlVcbvpzzz1nL774ojVr1swt89JLL7lEefbs2Xbuuedmdl8AAEAmde/e3b7++muXIR5NNdyig/ToeeH69evngnifSrwJvgEAiSzwztUUaCtxPvzww93fajemKunh1dZUJV2562lVWwOAXKdIEbMNGw5Meo0cIbqZFMw1/5o+fbp99NFHVqlSpdAhUUdqEl26vWHDhhSl4OFUHb148eIRE5Am7pcA4jzNDjTw3rlzp915553Wtm3bUKKrxLtAgQJWsmTJdFdbo60YgFxJJYJlyx6Y0igdzC00jJSaCiF+qORaJd1qGvbhhx9a9erVI97X3wq+w2uz7d692/WG3qhRo2zYYsStOLtfAtmNNDuBAm+Val911VW2f/9+e+KJJw66fFrV1mgrBgC5g+7lGpIKuYOGElNzL/VQrrG8lQGuaceOHe59pcvKbBkyZIjrm+Wbb76xDh06WJEiRVymOgAg9yLNjoPAW0F3mzZtXG+oyiUPr2KmnHPllm/cuDHd1dbUVkxV1v1pzZo1QWw2AGStXbsU2RyY9DoXU7ClUs5HHnnEBWOaVq1aZXPmzHGv33//fWvQoIGrYvzJJ5+45S+55JKIdSiAUw58eIKvMaJr1KhhhQsXdh13vv7666lug9KCU089NcX8448/3gYMGOBea5xpjTGtUTPU74h631Znn6nxt/+ff/4Jzfvyyy9D++dTUygNr6XtVFvlnj17Rox1nVuNHTvWpas6LxUrVgxNkydPDi3Tp08fd+66du3qzvGvv/5qM2fOdIE6kGXi6H4JZDfS7Pk5Ms3OE1TQ/eOPP7rO0kqXLh3x/kknnWT58+ePqLa2bt06l4ueWrU12ooByJVU8qsaP5rSUwqsRCG1aefO9C/7/6WVB102AxRwn3baadapUyd3z9YU3lmWgjPVTvruu+9cIJwed999txs3WsHfsmXL7NZbb7VrrrnGBfixtGvXzj7//HP76aefQvP0uaVLl7r3ZMuWLda+fXsX/H/22WdWs2ZNO//88938zNL61fGnxrFW52MKStUBmapo53bK/Ig16aHNp0yIgQMHunOuJmQ6P36v50C23S+B7EaaTZoddK/mGvprxYoVob9Vqq3SAY35qU7SLr/8cle68Pbbb7vhwfx223pfbbtVAtGxY0e77bbbXFCu+bfffrsb+9vv5RwAEtJhh6X+3vnnm82Y8e/f5cppDMfYyzZurKLcf/+uVs3szz9TLud56d403bt1D1cVY7/DrXCDBw92Jc3ppZznUaNGuXbFCuhFJd8KaJ966ilXUh1NwZ6CelWL1rCV8vLLL9vJJ59sRx99tPv77LPPjviM1qU+RRQsatjKzBgxYoSrVu23b1cw/+ijj7ptVKZBoUKFMrVeAEAuRppNmh10iffChQutfv36bhINF6LX9957r61du9b1iqr/TzjhhIhqa+E9lj/88MOuCqJKxk8//XT3IPfWW29Z3rx5M7o5AIAcQFWQM+Lbb791pacK1g877LDQ9MILL0SUaEdTybaCbVHJ7CuvvBIq7fabLXXp0sUF4v4Y0sow1pCWmaXROMaPHx+xnSoBVx8mynwGACA3Ic3OJSXeagemh53UpPWeT6UDGt9bEwDg/23dmvqhiM6Y1BBlqckTlaca1lY5KEWLFo3ahDwp0gM1RfIpaJUZM2bYEUcckaJ5UWpU8qzRMlSzSh2Aqc8PdeTpUxXpP/74w0aPHm1Vq1Z161KJuvoWiUXbKeHbGr6d/rZ27tzZtRGLVqVKlVS3FQAQx0izHdLsAANvAEBAooLXbFk2DapqriZE6VG2bFnXd0c4NUtSHx9y3HHHucRaJdGxqpWnRmNMq8MUlXor8FYTpfCOOdW2WyNpqF23KDD/M1Y1+7DtFLVf9oe51HaGO/HEE11b8qOOOird2wkAiHOk2QdFmn0Ix/EGAMSPatWquc7N1Nu3glm/1DoWtbVW0yRVHVdnm+p1PDwQV4/Y6t9DHapNmDDBVS9fsmSJPf744+7vtKhq+aRJk+y1115zHbuEU3D84osvuk7etK1aVr2apkbLq5M4dR72ww8/uBL4kSNHRizTt29f+/TTT93QWwrKtT9qVtWjR490HDUAAA490uxuOS7NJvAGAKSLAmX1xaHSapUUp9VuWm2g1QGaejtX52fqVfy6666LWOa+++5z/YOoN/RatWq5z6i/j+rVq6e5HVdccYX99ddftn379hRDlj3//PNuuEr1PXLttde66uHl1BFdKlQCr3bi33//vRvObNiwYXb//fdHLKMO3dQ5mxLvM888061b+6b+SwAAyIlIs3/McWl2kpeeRtk5zObNm12HORp7NHyMcADB6TdlabqWG9q6LqfBpxJhPzhVW+A8eVyHYuqQS8ElvWEnrrSug3hM4+JxnxD8/RLBI20/ONJt7MyiNJs23gAQFD04aigvAAD3SwAJjexEAAAAAAACROANAEHREFZ33HFgSmU4KwAA90sA8Y/AGwCCovGgH3rowBQ1NjQAgPslgMRB4A0AAAAAQIAIvAEgG6Q1BjbiH+cfAHKXXDgQFHJYmk2v5gBwCBUoUMDy5Mljv/32mxsLW38nJSVxDhLowW337t32xx9/uOtA5x8AkHPlz5/fpdO6byvdJs1OHF4Wp9kE3gBwCOnGrXEg161b54JvJKYiRYpYlSpV3PUAAMi58ubNa5UqVbK1a9faqlWrsntzkIvTbAJvADjElGOqG/jevXtt3759HP8EfIjLly8fpSYAkEscdthhVrNmTdtDR6kJJ28WptkE3gCQDXQDV/U1TQAAIOcHYJqAzCLwBoCgFC5s9s03/74GAHC/BJCQCLwBIChqC1S7NscXALhfAkhw9OoCAAAAAECAKPEGgKDs3m02ZMiB13fdpV7VONYAwP0SQAIi8AaAoKj300GDDry+4w4CbwDgfgkgQVHVHAAAAACAABF4AwAAAAAQIAJvAAAAAAACROANAAAAAECACLwBAAAAAAgQgTcAAAAAAAFiODEACEqhQmZffPHvawAA90sACYnAGwCCkjev2cknc3wBgPslgARHVXMAAAAAAAJEiTcABGX3brNHHjnw+pZbzAoU4FgDAPdLAAkowyXeH3/8sV144YWWnJxsSUlJNm3atIj3Pc+zgQMHuvcLFy5sTZo0sWXLlkUss2vXLuvRo4eVKVPGihYtahdddJGtXbv2v+8NAOQke/aY9elzYNJrAAD3SwAJKcOB97Zt26xevXo2ZsyYmO8PHz7cRo0a5d5fsGCBVahQwZo3b25btmwJLdOrVy+bOnWqTZo0yebNm2dbt261Cy64wPbt2/ff9gYAAAAAgNxe1bxly5ZuikWl3aNHj7b+/ftb69at3bwJEyZY+fLlbeLEida5c2fbtGmTPffcc/biiy9as2bN3DIvvfSSVa5c2WbPnm3nnnvuf90nAAAAAADis3O1lStX2vr1661FixaheQULFrTGjRvb/Pnz3d+LFi2yPXv2RCyjaul16tQJLQMAAAAAQLzI0s7VFHSLSrjD6e/Vq1eHlilQoICVLFkyxTL+56OpTbgm3+bNm7NyswEAAAAAyF3DianTtegq6NHzoqW1zNChQ61EiRKhSdXSAQAAAABIuMBbHalJdMn1hg0bQqXgWmb37t22cePGVJeJ1q9fP9c23J/WrFmTlZsNAAAAAEDuCLyrV6/uAutZs2aF5inInjt3rjVq1Mj9fdJJJ1n+/Pkjllm3bp198803oWWiqZ148eLFIyYAyPEKFTL76KMDk14DALhfAkhIGW7jraG/VqxYEdGh2pdffmmlSpWyKlWquKHChgwZYjVr1nSTXhcpUsTatm3rlldV8Y4dO9ptt91mpUuXdp+7/fbbrW7duqFezgEgLuTNa9akSXZvBQDkfNwvAcS5DAfeCxcutKZNm4b+7t27t/u/ffv2Nn78eOvTp4/t2LHDunbt6qqTN2zY0GbOnGnFihULfebhhx+2fPnyWZs2bdyy55xzjvtsXt10AQAAAABI5MC7SZMmriO01KiDtIEDB7opNYUKFbLHHnvMTQAQt/bsMXv66QOvb7rJLH/+7N4iAMiZuF8CiHNZOpwYACDM7t1m3bsfeN2hA4E3AKSG+yWAOBfIcGIAAAAAAOAAAm8AAAAAAAJE4A0AAAAAQIAIvAEAAAAACBCBNwAAAAAAASLwBgAAAAAgQATeABCUggXN3n77wKTXQA708ccf24UXXmjJycmWlJRk06ZNi3i/Q4cObn74dOqpp2bb9iJOcb8EEOcYxxsAArvD5jNr1Yrjixxt27ZtVq9ePbv++uvtsssui7nMeeedZ+PGjQv9XaBAgUO4hUgI3C8BxDkCbwAAEljLli3dlJaCBQtahQoVDtk2AQAQb6hqDgBB2bPHbPz4A5NeA7nUnDlzrFy5cnb00Udbp06dbMOGDdm9SYg33C8BxDlKvAEgKLt3m11//YHXV1xhlj8/xxq5jkrDr7jiCqtataqtXLnS7rnnHjv77LNt0aJFriQ8ll27drnJt3nz5kO4xciVuF8CiHME3gAAIFVXXnll6HWdOnWsQYMGLgifMWOGtW7dOuZnhg4daoMGDeKoAgDw/6hqDgAA0q1ixYou8P7xxx9TXaZfv362adOm0LRmzRqOMAAgoVHiDQAA0u2vv/5ygbQC8NSoCnpq1dABAEhEBN4AACSwrVu32ooVK0J/qx33l19+aaVKlXLTwIED3TBjCrRXrVpld911l5UpU8YuvfTSbN1uAAByEwJvAAAS2MKFC61p06ahv3v37u3+b9++vY0dO9aWLl1qL7zwgv3zzz8u+NaykydPtmLFimXjVgMAkLsQeAMAkMCaNGlinuel+v77779/SLcHAIB4ROANAEFRG9dXX/33NQCA+yWAhETgDQCB3WHzHRi/GwDA/RJAQmM4MQAAAAAAAkSJNwAEZe9es6lTD7xWD9AqAQcAcL8EkHB4CgSAoOzaZdamzYHXW7cSeAMA90sACYqq5gAAAAAABIjAGwAAAACAABF4AwAAAAAQIAJvAAAAAAACROANAAAAAECACLwBAAAAAAgQw4kBQFAKFDAbN+7f1wAA7pcAElKWl3jv3bvX7r77bqtevboVLlzYatSoYYMHD7b9+/eHlvE8zwYOHGjJyclumSZNmtiyZcuyelMAIHvlz2/WocOBSa8BANwvASSkLA+8hw0bZk8++aSNGTPGvvvuOxs+fLiNGDHCHnvssdAymjdq1Ci3zIIFC6xChQrWvHlz27JlS1ZvDgAAAAAA8VXV/NNPP7WLL77YWrVq5f6uVq2avfLKK7Zw4cJQaffo0aOtf//+1rp1azdvwoQJVr58eZs4caJ17tw5qzcJALLH3r1m779/4PW555rlo3UPAHC/BJCIsvwp8IwzznAl3j/88IMdffTR9tVXX9m8efNcsC0rV6609evXW4sWLUKfKViwoDVu3Njmz58fM/DetWuXm3ybN2/O6s0G4k6/KUvTtdzQ1nUD35aEpfvWBRcceL11K4E3AHC/BM8oSFBZHnj37dvXNm3aZMcee6zlzZvX9u3bZw888IBdffXV7n0F3aIS7nD6e/Xq1THXOXToUBs0aFBWbyoAAAAAALmvjffkyZPtpZdectXGFy9e7KqRP/TQQ+7/cElJSRF/qwp69Dxfv379XDDvT2vWrMnqzQYAAAAAIHeUeN9xxx1255132lVXXeX+rlu3rivJVql1+/btXUdqfsl3xYoVQ5/bsGFDilLw8KromgAAAAAAsEQv8d6+fbvlyRO5WlU594cT0zBjCr5nzZoVen/37t02d+5ca9SoUVZvDgAAAAAA8VXifeGFF7o23VWqVLHatWvbkiVL3NBhN9xwg3tf1cl79eplQ4YMsZo1a7pJr4sUKWJt27bN6s0BAAAAACC+Am+N133PPfdY165dXfXx5ORk11P5vffeG1qmT58+tmPHDrfMxo0brWHDhjZz5kwrVqxYVm8OAAAAAADxFXgreNbQYf7wYbGo1HvgwIFuAoC4VaCA2Zgx/74GAHC/BJCQsjzwBgD8v/z5zbp143AAwMFwvwQQ57K8czUAAAAAAPAvSrwBICj79pl98smB12eeqSEeONYAwP0SQAIi8AaAoOzcada06YHXW7eaFS3KsQYA7pcAEhBVzQEAAAAACBCBNwAAAAAAASLwBgAAAAAgQATeAAAAAAAEiMAbAAAAAIAAEXgDAAAAABAghhMDgKDkz282fPi/rwEA3C8BJCQCbwAISoECZnfcwfEFAO6XABIcVc0BAAAAAAgQJd4AEJR9+8wWLz7w+sQTzfLm5VgDAPdLAAmIwBsAgrJzp9kppxx4vXWrWdGiHGsA4H4JIAFR1RwAAAAAgAAReAMAAAAAECACbwAAAAAAAkTgDQAAAABAgAi8AQAAAAAIEIE3AAAAAAABYjgxAAhK/vxmAwb8+xoAwP0SQEIi8AaAoBQoYDZwIMcXALhfAkhwVDUHAAAAACBAlHgDQFD27zf77rsDr2vVMstDXicAcL8EkIgIvAEgKDt2mNWpc+D11q1mRYtyrAGA+yWABETxCwAAAAAAASLwBgAggX388cd24YUXWnJysiUlJdm0adMi3vc8zwYOHOjeL1y4sDVp0sSWLVuWbdsLAEBuROANAEAC27Ztm9WrV8/GjBkT8/3hw4fbqFGj3PsLFiywChUqWPPmzW3Lli2HfFsBAMitAgm8f/31V7vmmmusdOnSVqRIETvhhBNs0aJFoffJPQcAIGdo2bKl3X///da6desU7ym9Hj16tPXv39+9X6dOHZswYYJt377dJk6cmC3bCwBAbpTlgffGjRvt9NNPt/z589u7775r3377rY0cOdIOP/zw0DLkngMAkPOtXLnS1q9fby1atAjNK1iwoDVu3Njmz5+frdsGAEBC92o+bNgwq1y5so0bNy40r1q1aqnmnotyz8uXL+9yzzt37pzVmwQAADJBQbcojQ6nv1evXp3q53bt2uUm3+bNmzn+AICEluUl3tOnT7cGDRrYFVdcYeXKlbP69evbM888859yz5V4K9EOnwAgx8uf3+z22w9Meg3kUup0LZwy0aPnhRs6dKiVKFEiNClDHkgT90sAcS7LA++ff/7Zxo4dazVr1rT333/funTpYj179rQXXnjhoLnn/nvRSMAB5EoFCpiNGHFg0msgl1FHahKdPm/YsCFFOh6uX79+tmnTptC0Zs2awLcVuRz3SwBxLssD7/3799uJJ55oQ4YMcaXdqjreqVMnF4xnNvecBBwAgEOvevXqLvieNWtWaN7u3btt7ty51qhRo1Q/p5psxYsXj5gAAEhkWd7Gu2LFinbcccdFzKtVq5a98cYbKXLPtWx6cs+VgGsCgFxl/36zX3458LpKFbM8jOCInGfr1q22YsWKiCZhX375pZUqVcqqVKlivXr1cpnpqsmmSa81Yknbtm2zdbsRZ7hfAohzWR54q0fz5cuXR8z74YcfrGrVqilyz1UiHp57ro7ZACBu7Nihm96B11u3mhUtmt1bBKSwcOFCa9q0aejv3r17u//bt29v48ePtz59+tiOHTusa9eubuSShg0b2syZM61YsWIcTWQd7pcA4lyWB9633nqrq36mHPE2bdrYF198YU8//bSbRNXJyT0HACBnaNKkiWvulRql2wMHDnQTAADIIYH3ySefbFOnTnXtsgcPHuxKuDV8WLt27ULLkHsOAAAAAEgUWR54ywUXXOCm1JB7DgAAAABIFPT0AwAAAABAgAi8AQAAAAAIEIE3AAAAAAC5rY03AEB32HxmXbv+/92W2y0ApP5Eyv0SQHzjSRAAglKwoNnjj3N8AYD7JYAER1VzAAAAAAACRIk3AATF88z+/PPA6zJlNJYixxoAuF8CSEAE3gAQlO3bzcqVO/B661azokU51gDA/RJAAqKqOQAAAAAAASLwBgAAAAAgQATeAAAAAAAEiMAbAAAAAIAAEXgDAAAAABAgAm8AAAAAAALEcGIAENgdNp9Z+/b/vgYAcL8EkJB4EgSAoBQsaDZ+PMcXALhfAkhwVDUHAAAAACBAlHgDQFA8z2z79gOvixQxS0riWAMA90sACYgSbwAIioLuww47MPkBOACA+yWAhEPgDQAAAABAgAi8AQAAAAAIEIE3AAAAAAABIvAGAAAAACBABN4AAAAAAASIwBsAAAAAgAAxjjcABCVvXrPLL//3NQCA+yWAhETgDQBBKVTI7LXXOL4AwP0SQIKjqjkAAAAAAAEi8AYAAAAAIDcH3kOHDrWkpCTr1atXaJ7neTZw4EBLTk62woULW5MmTWzZsmVBbwoAHFrbtpklJR2Y9BoAwP0SQEIKNPBesGCBPf3003b88cdHzB8+fLiNGjXKxowZ45apUKGCNW/e3LZs2RLk5gAAAAAAED+B99atW61du3b2zDPPWMmSJSNKu0ePHm39+/e31q1bW506dWzChAm2fft2mzhxYlCbAwAAAABAfAXe3bp1s1atWlmzZs0i5q9cudLWr19vLVq0CM0rWLCgNW7c2ObPnx9zXbt27bLNmzdHTAAAAAAAJOxwYpMmTbLFixe7auTRFHRL+fLlI+br79WrV6faTnzQoEFBbCoAAAAAALmrxHvNmjV2yy232EsvvWSFNIZtKtThWjhVQY+e5+vXr59t2rQpNOk7AAAAAABIyBLvRYsW2YYNG+ykk04Kzdu3b599/PHHrjO15cuXh0q+K1asGFpGn4kuBQ+viq4JAAAAAABL9MD7nHPOsaVLl0bMu/766+3YY4+1vn37Wo0aNVwv5rNmzbL69eu793fv3m1z5861YcOGZfXmAED2yZvX7Pzz/30NAOB+CSAhZXngXaxYMddTebiiRYta6dKlQ/M1pveQIUOsZs2abtLrIkWKWNu2bbN6cwAg+6i5zYwZnAEA4H4JIMEF0rnawfTp08d27NhhXbt2tY0bN1rDhg1t5syZLmgHAAAAACCeHJLAe86cORF/qxO1gQMHugkAAAAAgHgW2DjeAJDwtm1TW5sDk14DAGLjfgkgzmVLVXMASBjbt2f3FgBA7sD9EkAco8QbAAAAAIAAEXgDAAAAABAgAm8AAAAAAAJE4A0AAAAAQIAIvAEAAAAACBC9mgNAUPLkMWvc+N/XAADulwASEk+CABCUwoXN5sw5MOk1kAsNHDjQkpKSIqYKFSpk92Yh3nC/BBDnKPEGAABpql27ts2ePTv0d968eTliAABkAIE3AABI+2EhXz5KuQEA+A+oag4AQdm2zaxs2QOTXgO51I8//mjJyclWvXp1u+qqq+znn3/O7k1CvOF+CSDOUeINAEH680+OL3K1hg0b2gsvvGBHH320/f7773b//fdbo0aNbNmyZVa6dOmYn9m1a5ebfJs3bz6EW4xci/slgDhGiTcAAEhVy5Yt7bLLLrO6detas2bNbMaMGW7+hAkTUv3M0KFDrUSJEqGpcuXKHGEAQEIj8AYAAOlWtGhRF4Sr+nlq+vXrZ5s2bQpNa9as4QgDABIaVc0BAEC6qQr5d999Z2eeeWaqyxQsWNBNAADgAEq8AQBAqm6//XabO3eurVy50j7//HO7/PLLXZvt9u3bc9QAAEgnSrwBAECq1q5da1dffbX9+eefVrZsWTv11FPts88+s6pVq3LUAABIJwJvAAhKnjxmDRr8+xrIhSZNmpTdm4BEwP0SQJwj8AaAoBQubLZgAccXALhfAkhwFMEAAAAAABAgAm8AAAAAAAJE4A0AQdm+3axatQOTXgMAuF8CSEi08QaAoHie2erV/74GAHC/BJCQKPEGAAAAACBABN4AAAAAAASIwBsAAAAAgADRxhvIhH5TlqZruaGt63J8AQAAEhjPjRBKvAEAAAAAyE2B99ChQ+3kk0+2YsWKWbly5eySSy6x5cuXRyzjeZ4NHDjQkpOTrXDhwtakSRNbtmxZVm8KAGSvpCSz4447MOk1AID7JYCElOWB99y5c61bt2722Wef2axZs2zv3r3WokUL27ZtW2iZ4cOH26hRo2zMmDG2YMECq1ChgjVv3ty2bNmS1ZsDANmnSBEzZSpq0msAAPdLAAkpy9t4v/feexF/jxs3zpV8L1q0yM466yxX2j169Gjr37+/tW7d2i0zYcIEK1++vE2cONE6d+6c1ZsEAAAAAED8tvHetGmT+79UqVLu/5UrV9r69etdKbivYMGC1rhxY5s/f37QmwMAAAAAQPz0aq7S7d69e9sZZ5xhderUcfMUdItKuMPp79WrV8dcz65du9zk27x5c5CbDQBZY/t2s5NPPvB6wQKqmwMA90sACSrQwLt79+729ddf27x581K8lxTV0ZCC9Oh54R22DRo0KLDtBIBAeJ7Zt9/++xoAwP0SQEIKrKp5jx49bPr06fbRRx9ZpUqVQvPVkVp4ybdvw4YNKUrBff369XNV1v1pzZo1QW02AAAAAAA5O/BWybVKuqdMmWIffvihVa9ePeJ9/a3gWz2e+3bv3u16Q2/UqFHMdaoNePHixSMmAAAAAAASsqq5hhJT7+RvvvmmG8vbL9kuUaKEG7Nb1cl79eplQ4YMsZo1a7pJr4sUKWJt27bN6s0BAAAAACC+Au+xY8e6/5s0aZJiWLEOHTq413369LEdO3ZY165dbePGjdawYUObOXOmC9QBAAAAAIgn+YKoan4wKvUeOHCgmwAAAAAAiGeB9moOAAlNIzVUrfrvawAA90sACYnAGwCCUqSI2apVHF8A4H4JIMEFNpwYAAAAAAAg8AYAAAAAIFCUeANAUHbsMDv55AOTXgMAuF8CSEi08QaAoOzfb7Zw4b+vAQDcLwEkJEq8AQAAAAAIEIE3AAAAAAABIvAGAAAAACBABN4AAAAAAASIwBsAAAAAgADRqzkABKlMGY4vAHC/BJDgCLwBIChFi5r98QfHFwC4XwJIcFQ1BwAAAAAgQATeAAAAAAAEiMAbAIKyY4dZkyYHJr0GAHC/BJCQaOMNAEHZv99s7tx/X2eDflOWpmu5oa3rBr4tAOLXf77X5ID7JQAEiRJvAAAAAAACROANAAAAAECACLwBAAAAAAgQgTcAAAAAAAEi8AYAAAAAIED0ag4AQSpShOMLANwvASQ4Am8ACErRombbtnF8AYD7JYAER1VzAAAAAAACROANAAAAAECACLwBICg7d5q1anVg0msAAPdLAAmJNt4AEJR9+8zeeeff1wAA7pcAEhIl3gAAAAAAxGvg/cQTT1j16tWtUKFCdtJJJ9knn3ySnZsDAABSQZoNAEAuDLwnT55svXr1sv79+9uSJUvszDPPtJYtW9ovv/ySXZsEAABiIM0GACCXBt6jRo2yjh072o033mi1atWy0aNHW+XKlW3s2LHZtUkAACAG0mwAAHJh4L17925btGiRtWjRImK+/p4/f352bBIAAIiBNBsAgFzaq/mff/5p+/bts/Lly0fM19/r169PsfyuXbvc5Nu0aZP7f/PmzYdga4GUdm3fmq7Dkp3XaFZvY27Y5xxn27Z/X+u4ZEPP5py33Mf/DXmeZzlBRtNsId1OPP/5XpMD7peJKLvSiERLmxJtfxPJ5gyk2dk6nFhSUlLE39rg6HkydOhQGzRoUIr5qpoO5GQPW+JtY27Y52yRnGw5Gect59myZYuVKFHCcor0ptlCuo3/dK/J4ffLRJRdaUSipU2Jtr+JlmZnS+BdpkwZy5s3b4qc8g0bNqTIUZd+/fpZ7969Q3/v37/f/v77bytdunSqiX5uzzlRpsKaNWusePHilgjYZ85zvOLa5trOKAW0SsCTc0jwkdE0OxHT7UT+vWcljh/HMLtxDXL8gkyzsyXwLlCggBs+bNasWXbppZeG5uvviy++OMXyBQsWdFO4ww8/3OKdEu1ES7jZ58TAeU4MnOfMy0kl3RlNsxM53U7kaz8rcfw4htmNa5DjF0SanW1VzZUTfu2111qDBg3stNNOs6efftoNJdalS5fs2iQAABADaTYAAP9NtgXeV155pf311182ePBgW7dundWpU8feeecdq1q1anZtEgAAiIE0GwCA/yZbO1fr2rWrmxBJ1fMGDBiQoppePGOfEwPnOTFwnuMTafbBJeK1n5U4fhzD7MY1yPELUpKXU8YrAQAAAAAgDuXJ7g0AAAAAACCeEXgDAAAAABAgAm8AAAAAAAJE4H0IbNy40Q2dpjHeNOn1P//8k+Znfv/9d+vQoYMbjL1IkSJ23nnn2Y8//hixzK5du6xHjx5WpkwZK1q0qF100UW2du3a//zdWSEz37t161br3r27VapUyQoXLmy1atWysWPHht5ftWqVJSUlxZxee+210HLVqlVL8f6dd96Z6/ZXmjRpkmJfrrrqqv/83Tl1n//++293TR9zzDHuuq9SpYr17NnTNm3aFLGe7DjHQZ7nePstp/Y7HTFiRI7/LQe1zzn994yMGzp0qDuHvXr1SnO5l19+2erVq+fuaRUrVrTrr7/ejeqSaAYOHJji+q9QoUKan5k7d64bQ75QoUJWo0YNe/LJJy2RZfQYTpkyxZo3b25ly5Z141Jr+N7333/fElVmrkHf//73P8uXL5+dcMIJlqgyc/z0fNO/f383alXBggXtyCOPtOeff94SljpXQ7DOO+88r06dOt78+fPdpNcXXHBBqsvv37/fO/XUU70zzzzT++KLL7zvv//eu+mmm7wqVap4W7duDS3XpUsX74gjjvBmzZrlLV682GvatKlXr149b+/evZn+7uzaZ7nxxhu9I4880vvoo4+8lStXek899ZSXN29eb9q0ae597de6desipkGDBnlFixb1tmzZElpP1apVvcGDB0csF/5+btlfady4sdepU6eIffnnn3/+83fn1H1eunSp17p1a2/69OneihUrvA8++MCrWbOmd9lll0WsJzvOcZDnOd5+y9G/0+eff95LSkryfvrppxz/Ww5qn3P67xkZo7S5WrVq3vHHH+/dcsstqS73ySefeHny5PEeeeQR7+eff3Z/165d27vkkksS7pAPGDDA7Xv49b9hw4ZUl9fxKlKkiDu+3377rffMM894+fPn915//XUvUWX0GOrYDRs2zF2vP/zwg9evXz93DJXOJKKMHj+f7tM1atTwWrRo4dLmRJWZ43fRRRd5DRs2dM83K1eu9D7//HPvf//7n5eoCLwDpsRC+RufffZZaN6nn37q5imgjmX58uXu/W+++SY0Tw+qpUqVcgmPfxPQzXPSpEmhZX799VeXwL/33nuZ/u6skNnv1Y9ZD9nhTjzxRO/uu+9O9TMnnHCCd8MNN0TM08P6ww8/7B0qQe6vHtTTeqhLhHP86quvegUKFPD27NmTbec4yH2Ox99ytIsvvtg7++yz01wmJ/yWg97nnPp7RsYo80cZgnqQPNg5HTFihHtgD/foo496lSpVSsiH9owELX369PGOPfbYiHmdO3d2BROJKqPHMJbjjjvOZXQmoswevyuvvNKl2Vlx/HOzjO7/u+++65UoUcL766+/At2u3ISq5gH79NNPXXXBhg0bhuadeuqpbt78+fNTrZYhqlrly5s3rxUoUMDmzZvn/l60aJHt2bPHWrRoEVpG1dLr1KkTWm9mvjsrZPZ7zzjjDJs+fbr9+uuvyhCyjz76yH744Qc799xzYy6vY/Dll19ax44dU7w3bNgwK126tKsS9MADD9ju3bstt+6vqimqCnLt2rXt9ttvty1btvzn784t51hUzVxV5FTFK7vOcZD7HI+/5ehmMzNmzIj5O81pv+VDsc858feMjOnWrZu1atXKmjVrdtBlGzVq5JqNvPPOO+73r2vj9ddfd59PRGoyp/tb9erVXTOLn3/+OdVl9XsIvy+K7psLFy5098xElZFjGG3//v3unlOqVClLVBk9fuPGjbOffvrJBgwYcMi2MV6On559GjRoYMOHD7cjjjjCjj76aJfu7dixwxJV5JMsstz69eutXLlyKeZrnt6L5dhjj3VtIfr162dPPfWUa/M5atQot/y6detC61UgXrJkyYjPli9fPrTezHx3Vsjs9z766KPWqVMn1xZWQVaePHns2WefdYFLLM8995xrL6sHm3C33HKLnXjiie7YfPHFF+44rly50q0rt+1vu3bt3M1NbWi++eYbty9fffWVzZo16z99d245x2oHed9991nnzp2z9RwHuc/x+FsON2HCBCtWrJi1bt061WVyym856H3Oqb9npN+kSZNs8eLFtmDBgnQtr2tamS1XXnml7dy50/bu3ev6cHjssccS7rArQ+mFF15wD9/KgLj//vvd8Vm2bJnLXIuma173wXD6W8fwzz//dO3lE01Gj2G0kSNH2rZt26xNmzaWiDJ6/BRkql+RTz75JEXmfyLK6PFTUK4CQxUkTp061f1uu3bt6vrzSdR23pR4Z2EHA9GTcmVFr6Mp5zvWfMmfP7+98cYbrlRMuZLqkGXOnDnWsmVLV/Kdluj1ZvS7s2uf/QDls88+czlkKgFTAqEf6OzZs1Msq9yyiRMnxixRuvXWW61x48Z2/PHH24033ug6Y9GDfUY7s8kJ+6uATaUqKv1UzqJKSvS+Hvzi/Rxv3rzZlQodd9xxKXKas+oc57R9Tmu9uek8h1PiqoAzvAbPof4t55R9PtS/Z2StNWvWuMygl156KdXrOdq3337rOoi899573e//vffec5lHXbp0SbjTo2eYyy67zOrWret+B6oV4mdUpSb6utdvIdb8RJGZY+h75ZVX3H1w8uTJMTP4EkFGjt++ffusbdu2NmjQIBdoIuPXn2pY6LeqzMdTTjnFzj//fFeQOH78+IQt9Sb7JpPUS3F0b7TR1CPv119/7XKFov3xxx8pcnLDqRdPVb1UNVtVrVSPlMppUpUNUYmJ5qsH3PCSsg0bNoRKjbRMZr47O/ZZP8C77rrL5Yj5VfD0sK1j8NBDD6Wo0qcH1u3bt9t111130O1WVU1ZsWJFunKEc+L++lT6p4wZ5cLqdbyeY1WFU0/+hx12mFte+xzEOc4J+xxvv+VwKiVYvny5e9BLzaH4Lee0fT5Uv2dkLQXO+l0qfQ5/OP/4449tzJgxrplYdOa4ej4//fTT7Y477gj9/lWL7cwzz3SlRYlYauvTcdADfPSILT79HqJreuj4q+Qxo7//RD2GPt2PlLmpUSPS00QiUaR1/PQcoszYJUuWuPTDDySV+aNrcObMmXb22WdbIjvY9af7m6qYq7mUr1atWu4YqglOzZo1LdEQeGeS2uhpOhgN3aDgWdUkldsjn3/+uZsXXa0yFv9i1UWtG4Cq3YoSfj2wqYqiX2VI1dBVfVFtKbLiuw/lPqu9liZVww2nhxjd6KKp1EvV9ZQhcTC6aUpGH3By0v76VJ1Hn/P3JR7PsUq61Y5Pw06olDg9JUuZPcc5YZ/j7bcc/TvV/mkopbSWCfq3nNP2+VD9npG1zjnnHFu6dGnEPA0NpuZhffv2jVkjTZlK0VVU/eX80ttEpYyK7777zmVCxKLfw1tvvRUxT8GOCiAOlhmbKA52DP2S7htuuMH9n6h9C2Tm+Klvmejf+xNPPGEffvihyzBWs6FEd7DrT5mOyuzRsKoqSBHV5tUzkZrfJaTs7t0tEWh4GA05ot5pNdWtWzfF8DDHHHOMN2XKlIienDX8kIai0bBD6t1XwyyF0xBE6hl19uzZbmgI9aAbawiig313Ttln9Q6rHqC13xpGZNy4cV6hQoW8J554IuJzP/74oxumR70lRtPwO6NGjfKWLFni1jF58mQvOTnZDWeQ2/ZXw2mp59EFCxa4IRhmzJjhenitX79+3J7jzZs3u2EntC7tf/iQFf4+Z9c5Dmqf4/G3LJs2bXJDAY0dOzbVdefE33JQ+5zTf8/InOheze+8807v2muvDf2t33u+fPnc713p+bx587wGDRp4p5xySsId8ttuu82bM2eO+z2r535d18WKFfNWrVoV89j5w4ndeuutrsf/5557LuGHE8voMZw4caK7/h5//PE0hzFMFBk9ftESvVfzjB4/jQChZ5vLL7/cW7ZsmTd37lw3IoSGWU1UBN6HgLrRb9eunbs4Nen1xo0bI0+EmUugfRrzUxerhhnS+N0axmDXrl0Rn9mxY4fXvXt3N8xY4cKF3Q/gl19+yfB355R9VmLQoUMH93CtwEQPtiNHjnTjmofTOJQ6Nvv27UvxvYsWLXKBm4Yv8NehG+W2bdty3f7qXJ511lnu/Go4LY0F3bNnzxTDMsTTOVZwqs/EmhSsZOc5Dmqf4/G3LBqvXPuS1gNeTvwtB7XPOf33jKwJvNu3b+/mRQ8fpiGcdG1UrFjRndO1a9cm3CHXkEzafz3X6H6owgQ9jKd17PSQr8wp/WY0bnpaGXmJIKPHUK9jpadaLhFl5hoMl+iBd2aO33fffec1a9bM3f8qVark9e7d29u+fbuXqJL0T3aXugMAAAAAEK/o1RwAAAAAgAAReAMAAAAAECACbwAAAAAAAkTgDQAAAABAgAi8AQAAAAAIEIE3AAAAAAABIvAGAAAAACBABN4AAAAAAASIwBtIYB06dLBLLrkkuzcDAIC41KRJE+vVq1d2bwaAHCBfdm8AAAAAEI+mTJli+fPnT9eyq1atsurVq9uSJUvshBNOsNySgf/PP//YtGnTsntTgByPwBsAAAAIQKlSpbLluO7ZsyfdAT+AQ4Oq5kACeP31161u3bpWuHBhK126tDVr1sy2bduWYrldu3ZZz549rVy5claoUCE744wzbMGCBaH358yZY0lJSTZjxgyrV6+eW6Zhw4a2dOnSiPXMnz/fzjrrLPd9lStXduuM9X0AACRKVfNq1arZkCFD7IYbbrBixYpZlSpV7Omnnw4tq9JuqV+/vktr9VnfuHHjrFatWi7dPfbYY+2JJ56IKCnX8q+++qr7jJZ56aWX3HvPP/+81a5d2woWLGgVK1a07t27hz63adMmu+mmm1yaX7x4cTv77LPtq6++Cr0/cOBAV/L+1FNPubS8SJEidsUVV7gSbv/9CRMm2Jtvvum+X5OeEwDERuANxLl169bZ1Vdf7RL67777ziWKrVu3Ns/zUizbp08fe+ONN1xCunjxYjvqqKPs3HPPtb///jtiuTvuuMMeeughF5Qrwb7oootc7rooCNdn9B1ff/21TZ482ebNmxeR2AMAkIhGjhxpDRo0cNXJu3btajfffLN9//337r0vvvjC/T979myXdquaujzzzDPWv39/e+CBB1w6ruD9nnvucWl1uL59+7qMbi2jdHjs2LHWrVs3F1wrbZ4+fbpL10XPAK1atbL169fbO++8Y4sWLbITTzzRzjnnnIg0f8WKFS6gf+utt+y9996zL7/80q1Tbr/9dmvTpo2dd955bns1NWrU6JAdSyDX8QDEtUWLFinC9latWpXivfbt23sXX3yxe71161Yvf/783ssvvxx6f/fu3V5ycrI3fPhw9/dHH33k1jVp0qTQMn/99ZdXuHBhb/Lkye7va6+91rvpppsivueTTz7x8uTJ4+3YsSOw/QQAIKdp3Lixd8stt7jXVatW9a655prQe/v37/fKlSvnjR071v29cuVKl8YuWbIkYh2VK1f2Jk6cGDHvvvvu80477bSIz40ePTpiGaXf/fv3j7ldH3zwgVe8eHFv586dEfOPPPJI76mnnnKvBwwY4OXNm9dbs2ZN6P13333Xpefr1q1L8RwBIG208QbinKqEKwdbVc2VA96iRQu7/PLLrWTJkhHL/fTTT67U+vTTTw/NU/uwU045xeWehzvttNMi2q8dc8wxoWWUa64c8pdffjm0jHLW9+/fbytXrnRV5QAASETHH3986LWqZleoUME2bNiQ6vJ//PGHrVmzxjp27GidOnUKzd+7d6+VKFEiYlmVpPu0zt9++82l/7Eord66datrfhZux44d7nnAp+rwlSpVikj/lZ4vX77cbTuA9CPwBuJc3rx5bdasWa7d9cyZM+2xxx5zVdY+//zziOX8qud6EIieHz0vFn8ZJcidO3d21d2iKQEHACBRRXd4prRT6WZq/PdU3Vx9qkSn7+GKFi0aeq0+VtKi9arNd6w22YcffvhB0/r0PBcAiETgDSQAJZAqydZ07733WtWqVW3q1KkRy6jdV4ECBVx77LZt27p5KgFfuHBhijFIP/vss1AQvXHjRvvhhx9cZy+iNmLLli0LtSMDAAAHpzRY9u3bF5pXvnx5O+KII+znn3+2du3apfswqvM2deb2wQcfWNOmTVO8r7Ra7bvz5cvnlkvNL7/84krOk5OT3d+ffvqp5cmTx44++ujQNodvL4DUEXgDcU4l20p4VcVcHaHpb1VdU5VvdX4WnlOuTl7UcZqqjyuwHj58uG3fvt1VcQs3ePBgVz1NDwQqPS9Tpoxdcskloc5dTj31VNf5iqrFab2qhq5Sd5W2AwCAlJRGq6RanZiperd6J1d1cvUerlpk6nm8ZcuWbgQSZYor47t3796pHkp9rkuXLm69+tyWLVvsf//7n/Xo0cONbqJq40q7hw0b5pqMKcBWR2ua51db1za0b9/edai6efNmtx3qUM2vZq6g/f3333dVz/VcoO1lGDMgNno1B+KcEuqPP/7Yzj//fJdDfffdd7teVZUIR3vwwQftsssus2uvvdblhqutthLU6PbgWu6WW26xk046yfViqp5S/Zx6tV+bO3eu/fjjj3bmmWe6YVHU+6qqtAEAgNhU+vzoo4+64btUwnzxxRe7+TfeeKM9++yzNn78eNdfS+PGjd1rf/ix1ChgHj16tBt6TEOKXXDBBS5t9mvCKcjW0J8a9UTPB1dddZUbmkyZ6j7VXtMoJXqGUAZ+nTp1IoYyUwa7gnYF6mXLlnWBPYDYktTDWirvAUAEtQVTlTXlsqfVBgwAAORuKjGfNm2aG0IMwH9HiTcAAAAAAAEi8AYAAAAAIEBUNQcAAAAAIECUeAMAAAAAECACbwAAAAAAAkTgDQAAAABAgAi8AQAAAAAIEIE3AAAAAAABIvAGAAAAACBABN4AAAAAAASIwBsAAAAAgAAReAMAAAAAYMH5P1XTBm5BQhK/AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", "\n", "for ax, idx, name, true_val in [\n", " (axes[0], 0, 'slope', theta_nl[0]),\n", " (axes[1], 1, 'intercept', theta_nl[1]),\n", "]:\n", " ax.hist(s_nl[:, idx], weights=w_nl, bins=40, density=True, alpha=0.6, label='IS posterior')\n", " ax.axvline(true_val, c='red', lw=1.5, ls='--', label='true value')\n", " ax.set(xlabel=name, title=f'NoisyLine IS — {name}')\n", " ax.legend()\n", "\n", "plt.suptitle(f'NoisyLine: IS posterior marginals (ESS={ess_nl:.0f}/8000)', y=1.01)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 5 SMC — annealing from prior to posterior\n", "\n", "SMC overcomes the low-ESS problem by bridging from the prior to the\n", "posterior through tempered targets $p_t \\propto \\text{prior} \\times \\mathcal{L}^{\\beta_t}$.\n", "When ESS drops below a threshold, particles are resampled and a Gaussian\n", "random-walk MCMC move is applied to maintain diversity." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "IS ESS=1.5/1000\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "d6e52c707f6049dca20bb490a2fc7c13", "version_major": 2, "version_minor": 0 }, "text/plain": [ "SMC: 0%| | 0/10 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "\n", "for ax, name, idx, true_val in [\n", " (axes[0], r'$\\mu$', 0, theta_true[0]),\n", " (axes[1], r'$\\sigma$', 1, theta_true[1]),\n", "]:\n", " ax.hist(s_is[:, idx], weights=w_is, bins=30, density=True,\n", " alpha=0.5, label=f'IS (ESS={eff_sample_size(w_is):.0f})', color='C0')\n", " ax.hist(s_smc[:, idx], weights=w_smc, bins=30, density=True,\n", " alpha=0.5, label=f'SMC (ESS={eff_sample_size(w_smc):.0f})', color='C2')\n", " ax.axvline(true_val, c='red', lw=1.5, ls='--', label='true value')\n", " ax.set(xlabel=name, title=f'IS vs SMC — {name}')\n", " ax.legend(fontsize=8)\n", "\n", "plt.suptitle('Broad prior: IS (poor ESS) vs SMC (annealed, high ESS)', y=1.01)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 6 2-D posterior scatter: SMC on NoisyLine" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "6304d2f106e44252b840afffc3dd1a99", "version_major": 2, "version_minor": 0 }, "text/plain": [ "SMC: 0%| | 0/10 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "\n", "sc = axes[0].scatter(smc_nl[:, 0], smc_nl[:, 1],\n", " c=w_nl2, cmap='viridis', s=8, alpha=0.7)\n", "axes[0].scatter(*theta_nl, c='red', s=100, zorder=5, label='true theta')\n", "plt.colorbar(sc, ax=axes[0], label='weight')\n", "axes[0].set(xlabel='slope', ylabel='intercept',\n", " title='SMC posterior — NoisyLine')\n", "axes[0].legend()\n", "\n", "axes[1].plot(nl.x, y_obs, '.', alpha=0.6, label='data')\n", "for s, i in smc_nl[np.random.choice(len(smc_nl), 80)]:\n", " axes[1].plot(nl.x, s * nl.x + i, alpha=0.05, c='C1', lw=0.8)\n", "axes[1].plot(nl.x, theta_nl[0] * nl.x + theta_nl[1], 'r-', lw=2, label='true line')\n", "axes[1].set(xlabel='x', ylabel='y', title='Posterior predictive draws')\n", "axes[1].legend()\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Summary\n", "\n", "| Method | When to use | ESS behaviour |\n", "|---|---|---|\n", "| `importance_sampling` | Proposal already close to posterior | Degrades quickly if proposal is too broad |\n", "| `sequential_importance_sampling` | Moderate mismatch, iterative refinement | Better than plain IS; no annealing schedule |\n", "| `SMC` | Broad prior, multi-modal, or expensive likelihood | Maintains high ESS via annealing + MCMC move |\n", "\n", "All three accept any `log_target_fn` / `log_likelihood_fn` and a `UniformPrior` (or any object with `.sample(n)` and `.log_prob(theta)`)." ] } ], "metadata": { "kernelspec": { "display_name": "beorn", "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.11.14" } }, "nbformat": 4, "nbformat_minor": 5 }