{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting a Line\n", "\n", "---\n", "\n", "In this notebook, we reproduce the classic \"Fitting a Model to Data\" example from emcee, https://emcee.readthedocs.io/en/latest/tutorials/line/, fitting a line with uncertainties, using both emcee and ```approxposterior```. This notebook demonstrates how ```approxposterior``` can be used to perform accurate Bayesian inference of model parameters given data with uncertainties." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import emcee\n", "from scipy.optimize import minimize\n", "import corner\n", "import george\n", "\n", "from approxposterior import approx, gpUtils as gpu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Data**\n", "\n", "---\n", "\n", "First, we generate data according to $y_i = m * x_i + b + \\epsilon_i$ where $\\epsilon_i$ are indepedent, Gaussian errors, for each measurement $i$. Then we'll plot the data to see what it looks like, with the grey line being the true, underlying model." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAHgCAYAAABXZE86AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdfVBca34n9u/Tb7x2nwYJSQgJAUISCAESFy40dHeR0WRn1jOZ2SpXXN4qr9+26iZbuxtPah3Hs37Jlj1Vd7NJnE15t1K+NfZ6k1V57dgT72TLlVlb4053Q4NAQkJIQq9ICEkIkOgDNNCvT/5AHIOuQLx1n+7T30/V1Fz96O7zG+6d++3nec55HiGlBBERERmDSe8GiIiI6OAw2ImIiAyEwU5ERGQgDHYiIiIDYbATEREZCIOdiIjIQCx6N7Abhw8fljU1NXq3QURElBHXr1+fk1JW7OY9ORXsNTU1GB4e1rsNIiKijBBCPNvtezgVT0REZCAMdiIiIgNhsBMRERkIg52IiMhAGOxEREQGwmAnIiIyEAY7ERGRgTDYiYiIDITBTkREZCAMdiIiIgNhsBMRERkIg52IiMhAGOxEREQGwmAnIiIyEAY7ERGRgTDYiYiIDITBTkREZCAMdiIiIgMxVLD39vait7dX7zaIiIh0Y6hgJyIiyncMdiIiIgNhsBMRERkIg52IiMhAGOxEREQGwmAnIiIyEAY7ERGRgTDYiYiIDITBTkREZCAMdiIiIgNhsBMRERkIg52IiMhAGOxEREQGwmAnIiIykLwIdh7nSkRE+cJQwa6qKiYnJxEKhfRuhYiISBeGCfZQKITR0VFMTEzg8uXLDHciIspLhgl2n8+HVCoFAIjFYvD5fPo2REREpAPDBHtvby9MprX/OTabjWvqRESUlwwT7C6XCy0tLaitrcXVq1fhcrn0bomIiCjjLHo3cJAURYGiKAx1IiLKWzk1YpdS6t0CERFRVsupYF9cXEQ4HNa7DSIioqyVU8EupURfXx8mJyf1boWIiCgr5VSwA0AqlcKtW7cwOjqqPd62V9yRjoiIjCangn39cTYAePbsGfr7+7G6uqpjR0RERNklp4LdbrejqqpK+/P8/Dz8fj/evn2rY1dERETZI6eCHQDa2trQ1NQEIQQAIBqNor+/H0+fPtW3MSIioiyQc8EOAHV1dejq6oLNZgOwdlPd7du3sbS0xEfiiIgor+VksAPA4cOH4fV6oSiKVltdXYWqqlhZWdGxMyIiIv3kbLADQFFREXp6enDy5Emtlkgk4Pf7MTc3p2NnRERE+tA92IUQZiHEiBDiP+7l/WazGRcvXkRzc7O27h6LxTAwMIAnT54A2N057XwEjoiIcpnuwQ7glwDc2++H1NTUwOFwaI/ESSlx584d/OEf/iHPaScioryha7ALIU4A+AaA7x/E51mtVjidTpSVlWm1H//4xzynnYiI8obeI/Z/CeBXAGy5hZwQ4jMhxLAQYnh2dvajH2gymdDd3Y1Tp04BwKYpeqvVyml2IiIyNN2ObRVCfBPAjJTyuhCid6vXSSm/APAFALS3t+/oWTaTyYSWlhY4nU6YTCbU1NQgEongl3/5l3H48OED6T9T1r+IcKaBiIh2Qs/z2HsAfEsI8RMACgE4hBD/Tkr5M3v9wPfDr7q6Gna7HaWlpSguLkZDQwPGx8cRDodx6dKlfTVPRESUjXSbipdSfldKeUJKWQPgpwH8eD+hvpWysjI4nU5YrVatNj09jUAggGQyedCXIyIi0pXea+wZYTKZ4HA4UFdXp9WWlpYQDocRi8V07IyIiOhgZUWwSyl9UspvpvMaQgg0NTXh0qVLmx6JW1hYwP3797kVLRERGUJWBHsmnThxAm63G0VFRVheXsbs7Cx++MMf4tq1a4jH43q3R0REtC963jynG0VRYLPZ8PTpU0gp8eu//uv43ve+h0gkgkQiAYslL38tRERkAHk3Yl/X19enTb8nEgncvn0bkUgEqqoiGo3q3B0REdHe5G2w9/b2amvtNpsNra2tANbW3RcXF3Hv3j2uuxMRUc7J22B3uVxoaWlBbW0tfvzjH+Ozzz5DSUmJ9vNHjx5hcHCQd80TEVFOydtgB9bW2qurq+FyuWC32+HxeGCz2bSfz87Owu/3Q1VVHbskIiLaubwO9vdZrVbY7XYUFxdrtZWVFQSDQUxNTenYGRER0c4w2N8jhEBxcTE6Ojq0u+NTqRRGRkYwNjamnRRHRESUjRjs71FVFZOTk5iYmIDH40Fpaan2s4mJCQwMDPCueSIiyloM9g1CoRBGR0cxMTGBy5cv4/bt2/B4PKisrNRe8+bNG/j9fszPz+vYKRER0Ycx2Dfw+XzaVHssFoPP54PFYkF7ezsaGhq0162urqK/vx+Tk5N6tUpERPRBeRHsPp9vR+eZv/9s+/pZ6ABw5swZdHZ2aqfEpVIp3Lp1C6Ojo2ldd19fGgiFQmm7BhERGUdeBPtObXy2/erVq3C5XJt+fuTIEXi9XjgcDq327Nkz9Pf3Y3V19cD7eX9pgOFOREQfw2B/z8Zn2z+kuLgYbrcbVVVVWm1+fh5+vx9v377Var29vZtG/HvxoaUBIiKi7TDY98BsNqOtrQ1NTU0QQgAAotEo+vv7MTExcWDX2W5pYONr9vsFgoiIjIPHmG1jPTC3GinX1dXB4XDg+vXriMVikFJibGwM4XAYUkot9PdqfWlAVVVcuXJly1kEIiKidRyx79Phw4fh9XqhKIpWm5qagqqqSCaT+/78jy0NEBERbcRgPwBFRUXo6enByZMntVoikYCqqpibm9v2vdk0lZ5NvRAR0d4w2A+I2WzGxYsX0dzcrE3Bp1IpDAwM4PHjxzp3R0RE+YJr7AespqYGDocDJpMJqVQKUkrcvXsXqqqitbUVZrM5I31sdX/Ax+4bICKi3MYRexqUl5fD6XRqm9kAwIsXLxAMBhGJRHTsjIiIjC6vg32nO9LthclkgsPhwKlTp7TawsICAoEAZmZm0nJNIiKivA72dBNCoKWlBa2trdrz6PF4HIODg3j48CGklDp3SERERsM19gyorq6Gw+HA0NCQtvXs+Pj4gT3vnm5clyciyh0csWeI0+mE1+vFoUOHtNr09DTC4fCBPO9OREQEMNgzqqCgAF1dXairq9NqyWQS4XAY09PTOnZGRERGwWDPMJPJhKamJrS1tWmPvkkpMTQ0hPHxca67ExHRvnCNXSdVVVWw2+0wm83aVPzDhw+hqira2to2PSpHRES0Uwz292TyBjGHwwFFUbC0tKTVZmZmEAgE0NHRAbvdnrFeiIjIGDgVnyaqqmJychKhUGjb15lMJtjtdtTX12u1SCSCQCCAly9fprtNIiIyGAZ7GoRCIYyOjmJiYgKXL1/+aLgLIdDY2Ij29nZt3T2ZTOL69euIRCLbrrvv9AsEERHlBwZ7Gvh8PqRSKQBALBbb8fR+ZWUlPB4PSkpKtNrKygoWFhYQi8W+9PrdfoEgIiLjY7BvY6+j4d7eXm2nOZvNtqujUO12OzweD44ePQoA+Pzzz/Fbv/Vb8Pv9UFV102v3+gWCiIiMi8G+hf2Mhl0uF1paWlBbW4urV6/C5XJt+doPfXmwWq3o6OjA2bNntdrKygqCwSCmpqa02n6+QOy0FyIiyi0M9i3sdzSsKAqqq6u3DfXtvjwIIXDu3Dl0dHTAYll7eCGVSmFkZARjY2NIpVK7+gLxMZzWJyIyBgb7Fg56NPwhO/nycOzYMXg8HpSWlmq1iYkJhEIhRKPRHX2BOKheiIgo+zHYt3CQo+Gt7PTLQ2lpKTweDyorK7Xa27dv4ff7EY/HP/ierabVt6pn4osMERGlH4N9Gwc1Gt7Kbr48WCwWtLe3o7GxUautrq5iYWFBOzFu3VbT6pcuXcLNmzc/ON2eiS8yRESUfgx2ne32y0N9fT26urq0LWellFhaWsKtW7e0qfStptU33lX/oen2dH+RISKi9GOw56CKigp4vV44HA6tNjk5if7+fqyurm45ra4oivZ6TrcTERkTgz1HFRcXw+12o6CgQKvNz8/D7/fj3LlzH5xWVxQFJSUlu55u52NwRES5g8Gew8xmM+x2O0pKSiCEAABEo1H09/fDZrN9cFrdYrHsarqdj8EREeUWBrsBFBUVoaurCzabDcDaunskEsHi4qJ2JOxe8TE4IqLcwmA3iMOHD8Pr9cLpdGq1aDSKvr4+rKysaLWlpSXcvHlzx5/Lx+CIiHILg91AioqK0NPTg5MnT2o1VVXh9/sxNze3p8/kY3BERLnFoncDdLBMJhMuXryI0tJSRCIRAGtT6AMDA5tG7ruhKAoURWGoExHlAAa7QRUWFsJsNqOgoADRaFRbd5dSbnu+OxER5TZOxRuY1WqF1+tFeXn5pnoymdRG80REZCwMdoMrLCyEy+VCTU2NVpNSIhAIYGZmRr/GiIgoLRjsecBkMqG5uXnTCXHxeByDg4N4+PAhp+aJiAyEa+xpko3PexcWFn6pNj4+jnA4jEuXLunQERERHTQGe47b7RcIIQQsFgsOHTqEN2/eAACmp6cRCASQTCZhNpvT0CUREWUKg30b2TjqPghCCHR1deHevXt48uQJgLWNa8LhMOx2u87dERHRfnCNPcf19vbueDc4VVWRSqWQSCRgMpnQ1NSEtrY2bZQupcTCwgLGx8e57k5ElKM4YtdZpmYF1g9zWX+ePRQKweVyoaqqCna7HUNDQ/j8888BAA8fPoSqqmhra9POfSciotzAEXue2HiYy/qf1zkcDni9XlRUVGi1mZkZBAIBLCwsZLLNPdvNzAURkZFxxG5Q788ErB/mkkqlYDKZvhSCVqsVnZ2dGB8fx6NHjwAAkUgEwWAQ0Wh007nvRESUvThizxPrh7kUFBSgpaXlg/u+CyHQ2NiI9vZ2bd09mUxicXFR246WiIiyG4M9jyiKgsLCQiiKsu3rKisr4fF4UFJSotVWVlYwMDCAWCyW7jaJiGgfGOx5JpFIYHJyEqFQaNvX2e12eDweHD16VKvNzc3B7/dDVdV0t0lERHvEYM8jqqoiEolgYmICly9f/mi4W61WdHR04Ac/+IF2x/zKygqCwSCmpqYy0TIREe0Sgz2PbBxpx2KxHT1qJ4TA2bNn8emnn8JiWbvXMpVKYWRkBGNjY5vutCciIv0x2PPIxrV1m822q8fDjh49Cq/Xu2lnuomJCYRCIUSj0YNsc1t8rI2IaHsM9jyiKApKSkpQW1uLq1evfvDO+O2UlJTA7XajsrJSq719+xZ+vx/z8/MH3S4REe0Bgz3PWCwWVFdX7zrUN76/vb0djY2NEEIAAFZXV9Hf349nz54dZKtERLQHugW7EOKkEOKvhRD3hBB3hBC/pFcvtHv19fXo7OzUtpxNpVIYHR3FrVu3uO5ORKQjPUfsCQD/RErZCKALwD8UQpzXsR/apYqKCni9XjgcDq02OTmJ/v5+rK6u6tgZEVH+0i3YpZSvpJQ33v31IoB7AKr06of2pri4GG63G1VVf/O3bn5+Hn6/XzvvnYiIMicr1tiFEDUALgEY/MDPPhNCDAshhmdnZzPdGu2A2WxGW1sbmpqatHX3aDSKUCiEiYkJnbsjIsovuge7EKIUwJ8B+I6U8ktHiUkpv5BStksp2zeePkZrVFXd0U5ymVBXVweXywWbzQZg7Xz3sbExjIyMIJlM6twdEVF+0DXYhRBWrIX6FSnlD/TsJRetn7G+053kMuHQoUPwer1wOp1abWpqCn19fVheXtaxMyKi/KDnXfECwO8DuCel/B29+shlG89Y3+lOcplQVFSEnp4eVFdXazVVVREIBJCu5ZRsmrkgItKTniP2HgB/D8BXhBA33/3nJ3TsJ+esn7EO7H4nuXQzmUxobW1FS0uL1mMsFsPg4CA6OzsPtNdsnLkgItKLRa8LSymDAIRe1zeC9TPWVVXFlStX9rzpTDqdOnUKdrsdw8PDiEajkFIiEokgkUggkUho+8/vx4dmLrLxd0FElAm63zxH+6Moyo53kvP5fLh48WIGutqsvLwcXq8X5eXlWi0ajSIYDCISiezqs27evImbN29uqmXzzAURUaYx2CkjCgsL4XK5UFNTo9UWFxcRCAQwMzOzr89en7nY6x74RERGwmCnjDGZTGhubkZpaan2vHs8Hsfg4CAePHgAKeWeP3s3MxdEREbGYKeMKywshKIoKCoq0mr379/H8PAw4vG4jp0REeU+BjvpwmKxwOPx4NChQ1ptenoawWAQS0tLOnZGRJTbGOx5xufzZc3z7gUFBXC5XKirq9NqS0tLCAQCmJ6e/uB7EokEVldXt32krbe3lzfQEVHeYrCTroQQaGpqQltbG8xmM4C18B4aGsL4+PimdfdQKIRIJIJoNMrn1YmItsBgp6xQVVUFt9uN4uJirfbw4UNcu3ZNW3ffONOQTTvtERFlEwY7ZQ2HwwGv14uNh/3MzMzA7/djYWFh0/T6+8+rf+j5diKifMRgp6xitVrR2dmJM2fOaLXl5WUEg0FUV1ejpKQEBQUFfF6diGgLDHbKOkIINDQ0oL29XdtyNplM4saNGxBCaDfdERHRlzHYKWtVVlbC7XajpKREq6VSKSSTScRiMR07IyLKXgx2ymp2ux0ejwdHjx7ValJK+P1+qKqqY2dERNmJwU5Zz2q1oqOjA+fOndNqKysrCAaDeP78uY6dERFlH92ObSXaDSEEzp49C7PZjGQyCWBtWv7mzZsIh8OQUmr7zxMR5TOO2CnjVFXF5OTknjaYMZlMsFgssNvtWu3p06dIJpP7OkSGiMgoOGLPcbm2SUsoFMLo6ChSqRQuX76868fWEokEEokEzGYzKisr8erVKwBr6+6JRALz8/Ppap2IKCdwxE4Z5fP5kEqlAOx+97iNW8p+7WtfQzweR2Nj46Yp+P7+fqyurh5020REOYPBThnV29sLk2ntH7v3d4/7mA9tKVtfX4/Ozk6cPn0atbW1SKVSWFpawtLSkvYFgogonzDYKaNcLhdaWlpQW1u762n4rbaUraiogKIo2mY2ALC6uoq+vj6O3oko7zDYKeMURUF1dfWud49zuVxbbilrNpuhKAqqqqq0Wjgcht/vx5s3bw6sdyKibMdgp5xisVhQWFj4wS8FQgi0tbUhlUphbm4O4+PjiEajCIVCmJiY0KFbIqLMY7CToYRCIdy/fx+vX7/Gb/zGb2hnuo+NjWFkZER7Bj7X9fb27ur+BCLKHwx2MpSNd90nEgk8evRI+9nU1BSCwSCWl5f1ao+IKO0Y7GQI65veHDp0aNNd9z/7sz+L6upq7XULCwvw+/2YnZ3Vq1UiorTiBjWUUy5evPil2sZNb77zne/g9OnTSCQSuHLlirYW73Q6MTY2hlQqhXg8jsHBQTQ0NKC+vj7T/xOIiNKKI3bKee9vepNIJL501/2pU6fQ3d2NwsJCAGs71d27dw/Xr19HIpHQpW8ionRgsFPOe3/TG0VRPvi6srIyeDwelJeXa7WXL18iGAwiEolkpFcionRjsFPOe3/Tm62CHYD2qFxNTY1WW1xcRCAQwOvXrzPQLRFRejHYyRB2s+mNyWRCc3MzLl68qI304/E4rl27hgcPHmTtKXF8xI2IdoLBTnnr5MmT6OnpQVFRkVa7f/8+hoeHEY/HdeyMiGjvGOyU15xOJ7xeLw4fPqzVpqenEQgEsLS0pGNnRER7w2CnvGez2dDV1YXTp09rtUgkgkAgoJ33TkSUKxjsRFjbZ/78+fNoa2uD2WwGsLZz3fDwsLYtLRFRLmCwE21QVVUFt9uN4uJirfbw4UMMDg5y3Z2IcgKDnTLO5/PB5/Pp3caWHA4HvF4vKioqtNrs7Cz8fj8WFhZ07IyI6OMY7EQfYLVa0dnZiTNnzmi15eVlBINBvHjxQsfOiIi2x73iKadkcqQvhEBDQwMURcHNmzeRSCSQTCZx48YNhMNhnD9/HkKIjPVDRLQTHLETfURlZSXcbjdKSkq02pMnTzAwMIBYLKZLT+un2YVCIV2uT0TZi8FOtAN2ux0ejwfHjh3TanNzc/D7/VBVNaO9rJ9mNzExgcuXLzPciWgTBjvRDlmtVrS3t+PcuXNabWVlBcFgEM+fP89YH++fZpfNNyISUeYx2Il2QQiBzz77DJ9//jmsVisAIJVK4ebNm7h9+7YWuOn0/ml23D+eiDZisJPhZOJxOpvNBo/HA7vdrtWePn2KUCiE1dXVtF77/dPsdnLwDRHlDwY70R6VlJTA7XajsrJSq719+xaBQADz8/NpvfZuTrMjovzCYCfaB4vFgvb2djQ2NmqPvq2urqK/vx/Pnj3TuTsiykd8jp0MQe8byOrr66EoCq5fv454PI5UKoXR0VGEw2E0Nzdra+JEROnGf9sQ7VNvby96e3tRUVEBr9cLh8Oh/WxychJ9fX1YWVnZ93X47DoR7QSDnegAFRcXw+1248SJE1otHA7D7/fjzZs3e/7c959dz/Sz80SUOxjsRAfMbDbj0qVLuHDhgrbuHovFEAqF8OTJkz195vvPrmci2NdnIogotzDYidKktrYWLpcLBQUFAAApJe7cuYMbN24gmUzu6rPef3ZdUZQD75eIjIHBTpRGhw4dgtfrhdPp1GovXrxAMBjE8vLyjj/n/WfXGexEtBUGO1GaFRYWoqenB9XV1VptYWEBfr8fs7OzO/4cPrtORDvBYCfKAJPJhNbWVrS0tGhT6vF4HIODg3j06JHO3RGRkTDYiTLo1KlT6O7uRmFhIYC1dfd79+5heHgYiURC5+6IyAgY7EQZVlZWBq/Xi/Lycq326tUrBINBRCKRA7kG72gnyl8MdiIdFBQUwOVyoaamRqstLi4iEAjg9evX+jVGRDmPwU6kE5PJhObmZly8eHHTuvu1a9fw4MEDSCl17jB9OKNAlD4MdiKdnTx5Em63G0VFRVrt/v37GB4eRjwe17EzIspFDHaiXUrHnu2KosDr9eLw4cNabXp6GoFAAIuLiwd2HSIyPp7uRrQL63u2p1IpXL58GVevXj2wz7bZbOjq6sK9e/fw+PFjAEAkEkEwGMTFixc3vVbv0+yIKHtxxE60C+/v2X7QASuEwPnz59HW1gaz2QwASCQSGB4exvLysqHX3YnoYDDYiXbh/T3b03UDWFVVFdxuN4qLi7Xa8vIyFhcXue5ORNtisBPtwvt7trtcrrSdk+5wOOD1enHkyBGtFovF4Pf7sbCwcKDXIiLjYLAT7dLGPdvfPyf9oMPdarXi008/xZkzZ7Ta8vIygsEgXrx4caDXIiJj0DXYhRBfF0LcF0I8EkL8qp69EO1FutfcgbV194aGBjgcDu1892QyiRs3buDOnTtpW3dP10wEEaWXbsEuhDAD+NcA/jaA8wD+rhDivF79EO1Fptbc1z/f6XSitLRUqz158gShUAjRaPRAr5XumQgiSh89R+yfAngkpXwipYwB+PcAvq1jP0S79qE193Tx+XwIBAJwu904duyYVn/z5g0CgQDC4fCBXivdMxFElB56BnsVgOcb/jz1rkaUUzJ9TrrVakV7ezvOnTun1VZWVtDX14fnz59v886d2+lMBLeGJco+em5QIz5Q+9JioRDiMwCfAUB1dXW6eyLKCUIInD17Fk6nEzdu3EA8HkcqlcLNmzcRDochpdTW4/difSZCVVVcuXIlY19aiGj/9ByxTwE4ueHPJwC8fP9FUsovpJTtUsr2ioqKjDVHlAuOHDkCj8cDu92u1Z4+fYqFhQVtKn2vMj0TQUQHQ89gHwJwRghRK4SwAfhpAD/UsR+inFRSUgK3243jx49rNVVV8eDBA/zoRz/KaC+cmifSn27BLqVMAPhHAH4E4B6AP5FS3tGrH6JcZrFY8Mknn+D8+fO4f/8+nj59itevX+Nb3/oWfvCDH+jdHhFlkK6HwEgp/wLAX+jZA5GRnD59WltjB9b2mf/zP/9znD59Gs3NzdoNcURkXPx/OZHBfOMb39AC3GKxoLm5GZOTk+jr68PKyorO3a3h5jdE6cNgJzKY9Tvaa2pq8Hu/93toaGgAAITDYfj9frx580bX/tK9+Q3X+SnfMdiJDEhRFJw6dQo///M/jwsXLmiPvsViMYRCITx58kS33rj5DVF67SjYhRD/QgjhEEJYhRBXhRBzQoifSXdzRLR/tbW1cLlcKCgoAABIKXHnzh3cuHEDyWQy4/1kchteony00xH735JSLgD4JtaePz8L4L9LW1dEdKAOHToEr9cLp9Op1V68eIFgMIjl5eWM9pLJbXiJ8tFOg9367r9/AsAfSSnfpqkfoqzn8/lycvq4sLAQPT09m3ZwXFhYgN/vx+zsbEZ74eY3ROmz02D/f4QQ4wDaAVwVQlQAWE1fW0SUDiaTCa2trWhpadGmw+PxOAYHB/Ho0SOduyOig7CjYJdS/ioAF4B2KWUcQAQ8iY0oZ506dQrd3d0oLCwEsLbufu/ePQwPDyORSOjcHRHtx05vnisE8AsA/i8hxJ8B+K8AHNwZkUSUcWVlZfB6vSgvL9dqr169QjAYRCQS0bEzItqPnU7F/x8AmgD8LoB/BaARwP+ZrqaIKDMKCgrgcrlQW1ur1RYXF+H3+xGLxXTsjIj2aqdbyp6TUrZu+PNfCyFupaMholyTizfSbWQymXDhwgUoioLR0VGkUikkEgksLCyguLh42yNgVVWFqqoIhUK8EY4oS+x0xD4ihOha/4MQohNAX3paIiI9nDx5Em63G0VFRVpteXkZQ0NDiMfjX3p9uneQI6K92XbELoS4DUBi7XG3nxVCTL778ykAd9PfHhFlkqIo8Hq9uH79Oj7//HMAwOvXrxEIBNDR0bHp3HfuIEeUnT42Yv8mgP8CwNcBtAH4XwD8DoBLAL6R3taISA82mw1dXV04ffq0VotEIggEAnj16pVW4w5yRNlp22CXUj6TUj4D8HewdrPcYQAV7/76W+lvj4j0IITA+fPn0dbWBrPZDABIJpMYHh7GvXv3IKXctIPcuXPn8N3vfpenthFlgZ3ePPf3AXRJKSMAIIT4HwGEsHaXPBEZVFVVFex2O4aGhrStZx89egRVVfHJJ59AURQoigJg7Ua69ZvvLl++zO1iie5mEHkAACAASURBVHSy05vnBICNp0Uk39WIyOAcDge8Xi+OHDmi1WZnZ+H3+zdtZqOqKtfcibLAToP93wAYFEL8MyHEPwMwAOD309YVEWUVq9WKTz/9FGfOnNFqy8vLUFUV0WgUwNqNd1xzJ9LfjqbipZS/I4TwAXBjbaT+C1LKkXQ2RkR7l47RshACDQ0NcDqdGBkZQSKRgJQSi4uLSCQScDgcaGlpgaqquHLlim7T8O8/W7/+BYMzCJQvdrrGDinlDQA30tgLEeWAY8eOwePxYGhoSKutrKwgkUjAbrdDURTdQn392fqN6/xE+WbHwU5EtK60tBQejwc2m03bejYejyMcDsPhcHz0/ekaPfPZeqKdr7ETEW1isVhgt9tRXFys1VKplPbImx74bD0Rg52I9kEIgeLiYjgcDi1QpZS4desWbt++rY2eM2Xjs/V83I7yFYOdiPbNZrNBURRYLH+zuvf06VOEQiGsrq5mtBdFUVBdXc1Qp7zFYCeiA2E2m6EoCgoKCrTa27dv4ff78fbtWx07I8ovDHYiOjBCCNjtdpw/f1476jUajSIUCuHp06f6NkeUJxjsRHTgTp8+jc7OTthsNgBrN9Xdvn0bN2/ezPi6O1G+YbATUVpUVFTA4/Foe8kDwPPnz9HX14eVlRUdOyMyNgY7EaVNcXExenp6cOLECa0WDofh9/vx5s0bHTsjMi4GOxGlldlsxqVLl3DhwgVt3T0WiyEUCuHJkycHco3e3l4+s070DoOdiDKitrYW3d3d2l3zUkrcuXMHN27cQDKZ/Mi7iWinGOxEtGc+nw8+n0/bbU5V1W1fX15eDq/Xi7KyMq324sULBINB7bx3ItofBjsR7cv6wSsTExMYHR39aLgXFhaiu7sb1dXVWm1hYQF+vx+zs7PpbpfI8BjsRLQvGw9eWd8r/mNMJhNaW1vR2tqqbUUbj8cxMDCAhw8fpqVPrsNTvmCwE9G+bDx4paioCFeuXNnxe6urq9Hd3Y3CwkKtNj4+juHhYSQSiQPvlSgfMNiJaF/2e/BKWVkZvF4vysvLtdqrV68QDAYRiUT23d9O1/+JjILBTkT7tt+DVwoKCuByuVBbW6vVFhcX4ff78fr16z33tdv1fyIjYLATUVYwmUy4cOECLl26pE3tJxIJXLt2Dffv34eUctefuZf1f6JcZ/n4S4iIMufEiROw2+0YGhrStp598OABVFXFpUuXYLVat32/z+fT/np9/T+VSsFkMm3a3pbIqDhiJ6KsoygKvF4vDh8+rNVev36NQCCAxcXFHX/OxvX/lpYWBjvlBQY7EWUlm82Grq4unD59WqtFIhEEAgG8evVqx5+zvv7PUKd8wWAnoqwlhMD58+fxySefwGw2AwCSySSGh4dx7969Pa27Exkdg52Ist7x48fhdrtRUlKi1R49eoTBwUHEYjEdOyPKPgx2IsoJDocDHo8HR44c0Wqzs7MIBALczIZoAwY7EeUMq9WKTz/9FGfOnNFqy8vLUFUV0WhUx86IsgeDnYhyihACDQ0N6OjogMWy9sRuJBLBkydPcOXKFe25daJ8xWAnopx07NgxeDwePHv2DE+fPsXMzAx+8Rd/Ed///vc5eqe8xmAnopxVWloKVVW1u+MTiQSCwSD8fj/C4bDO3RHpg8FORDntK1/5irYFrcViQXNzM1ZXV9HX14fJyUmduyPKPG4pS0Q5bX13OVVV8bu/+7swm82Ix+NIpVK4desWlpaWUFJSAiGE3q0SZQRH7ESU89Z3l/vGN74Bj8cDh8Oh/Wx1dRULCwsIh8OYnJxEKBTSsVOi9GOwE9G++Xy+TYev6KmkpARutxvHjx/Xaqqqase3Xr58meFOhsZgJyLDMZvN+OSTT3D+/HkIIRCJRLQb7GKxWNZ8CSFKBwY7ERnW6dOn4XA4UFpaqq2xm81mHD9+nM+7k2Ex2InI0KxWK44fP47Tp0/jyJEj+N73vofy8nL09fVp570TGQmDnYgMzefzIRAI4Pjx4zhx4gQaGhoAAOFwGH6/H3Nzczp3SHSwGOxElBeEELDb7bhw4YI2LR+LxTAwMIAnT57o3B3RwWGwE1Feqa2tRXd3NwoKCgAAUkrcuXMHN27cQDKZ1Lk7ov1jsBNR3ikvL4fX60VZWZlWe/HiBYLBIJaXl3f9eb29vejt7T3ADon2jsFORHmpsLAQ3d3dOHXqlFZbWFiA3+/H7OzsgV+P4U+ZwmAnorxlMpnQ0tKC1tZWbb/5eDyOgYEBPHz4UOfuiPaGwU5Eea+6uhrd3d0oLCzUauPj4xgeHkYikdCxM6LdY7ATEQEoKyuD1+tFeXm5Vnv16hUCgQCWlpZ07IxodxjsRETvFBQUwOVyoba2VqstLS0hEAhgenpax86Ido7BTkS0gclkwoULF3Dp0iVt3T2RSGBoaAj379/X9pwnyla6BLsQ4n8SQowLIUaFEP+3EMKpRx9ERFs5ceIE3G43ioqKtNqDBw8wNDSEeDyuY2dE29NrxP6XAC5IKVsAPADwXZ36ICIDSNexsYqiwOv1oqKiQqu9fv0agUAAi4uLB349ooOgS7BLKf+TlHL9VtMBACf06IOI6GNsNhs6OztRX1+v1SKRCAKBAF6+fKljZ0QfZtG7AQC/COCP9W6CiGgrQgg0NjZCURTcvHkTyWQSyWQS169fh6qqkFJq+88T6S1twS6E+CsAxz7wo1+TUv6Hd6/5NQAJAFe2+ZzPAHwGrD1rSkSkl+PHj8Nut2NoaAiRSAQA8OjRIywuLqK0tFTn7ojWpC3YpZRf3e7nQoifA/BNAJflNreZSim/APAFALS3t/N2VCLak4Nag7fb7fB4PLhx4wZmZmYArJ0Sp6oqVFWFoigHch2ivdLrrvivA/jvAXxLSrn7ExeIiHRktVrx6aef4uzZs1otmUwiGAxiampKx86I9Ftj/1cACgD85bt1qQEp5X+tUy9ERLsmhMC5c+egKAqEEJBSIpVKYWRkBKqqorGxUXsOniiTdAl2KWX9x19FRJT9jh07BqfTiYWFBa325MkTqKqKTz75RDv3nShT+HWSiGifzGYznE4njh37m/uF37x5A7/fj3A4rGNnlI8Y7ERE+6SqKp4/f454PI6Ghgatvrq6ir6+PkxOTkJVVUxOTiIUCunYKeUDBjsR0T6EQiGMjo5iYmICX/3qVzE3N4fOzk5YrVYAQCqVwh//8R/j1q1bmJiYwOXLlxnulFYMdiKiffD5fEilUgDWHnvz+Xw4cuQIPB4PHA4HAOD27dva4THrrzGa3t5e9Pb26t0GgcFORLQvvb292t3vNptNC7eSkhK43W5UVVWhublZ25nOYrGgra1Nr3YpDzDYichwMjl6dLlcaGlpQW1tLa5evQqXy6X9zGw2o62tDT/5kz+J2tpaHDlyBL/927+NeDyOp0+fZqQ/yj/ZsFc8EVFOUxQFiqJsCvWN6urqcOzYMSwuLqKhoQFSSty+fRvhcBjNzc0wm80Z7piMjCN2IqIMsFqt2heAdc+fP0dfXx9WVlZ07IyMhsFORJQhZrMZPT09OHnypFZTVRV+vx9zc3N7+kzetEbvY7ATEWWQ2WzGxYsXN91QF4vFMDAwgCdPnujcHRkBg52ISAc1NTXo7u7WtpyVUuLOnTu4ceMGksmkzt1RLmOwExHppLy8HF6vF2VlZVrtxYsXCAaD2nnvRLvFYCci0lFhYSG6u7tx6tQprbawsIBAIKCd9060Gwx2IiKdmUwmtLS0oLW1VdvsJh6PY3BwEA8fPtS5O8o1DHYioixRXV2N7u5uFBYWarXx8XEMDQ0hkUjo2BnlEgY7EVEWKSsrg9frxaFDh7Ta9PQ0AoEAlpaWdOyMcgWDnYgoyxQUFKCrqwt1dXVabWlpCYFAANPT0xnrg8/I5yYGOxFRFjKZTGhqasKlS5e0dfdEIoGhoSHcv39fOy2O6H0MdiKiLHbixAm43W4UFRVptQcPHuDatWuIx+M6dkbZisFORJTlFEWB1+tFRUWFVpuZmUEgEOBNdfQlDHYiohxgs9nQ2dmJ+vp6rRaJRKCqKqLRqI6drVFVFZOTkwiFQnq3kvcY7ERE++Tz+eDz+dJ+HSEEGhsb8cknn2hHvUopsbi4iHv37um27h4KhTA6OoqJiQlcvnyZ4a4zBjsRGU42jh4PMvyPHz8Oj8eDkpISrfbo0SMMDg4iFosdyDV2w+fzIZVKAVg70CYTX3Joawx2IjKUfBk92u12eDwe2Gw2rTY7Owu/3w9VVTPaS29vr3bnvs1m4yNyOmOwE5Gh5NPo0Wq1QkqJcDiM8fFxAMDKygqCwSCmpqYy1ofL5UJLSwtqa2tx9epVuFyujF2bvozBTkSGkk+jx/XZiRcvXuA3f/M3tX3lU6kURkZGMDY2pn3JSTdFUVBdXc1QzwIMdiIylHwaPW6cnYjH41hYWEBpaan284mJCQwMDGTFXfOUOQx2IjKcfBk9vj878bWvfQ0ejweVlZXaa968eQO/34/5+Xm92qQMY7ATEeWoD81OWCwWtLe3o6GhQXvd6uoq+vv7MTk5qWO3lCkMdiKiHLbV7MSZM2fQ2dkJq9UKYG3d/datWxgdHc3Yujvpg8FORGRQR44cgdfrhcPh0GrPnj1Df38/VldXdeyM0onBTkRkYMXFxXC73aiqqtJq8/Pz8Pv9ePv2rY6dUbow2ImIssxBn4NuNpvR1taGpqYmCCEAANFoFP39/ZiYmNjyfdm4gx99HIOdiChP1NXVoaurS9utTkqJsbExjIyMIJlMbnptvuzgZ0QMdiKiPHL48GF4vV4oiqLVpqam0NfXh5WVFa2WTzv4GQ2DnYgozxQVFaGnpwcnT57Uaqqqwu/3Y25uDkB+7eBnNAx2IqI8ZDabcfHiRTQ3N2vr7rFYDAMDA3j8+HFe7eBnNBa9GyAiIv3U1NTA4XBgeHgY0WgUUkrcvXsXqqrC4XBAURSGeo7hiJ2IKM+Vl5fD6/WirKxMq7148QKqqn7ppjrKfgx2IiJCYWEhuru7cerUKa2WSCSgqipmZmZ07Ix2i1PxREQEADCZTGhpaYHT6cTt27cBrG1FOzg4iIaGBtTX12vr8e/jXfPZg8FORESbVFdXw+FwwGQyaY+8jY+PIxwO49KlS7BYGB3ZjH93iMhwOHrcP6fTCafTicXFRa02PT2NQCCAjo6OTee+U3bhGjsREX2QyWSCw+FAXV2dVltaWkIgEMD09LSOndF2OGInIsph6Z6dEEKgqakJTqcTt27dQjKZRCKRwNDQEM6cOYNz585tue5O+uCInYiIPqqqqgputxvFxcVa7eHDh7h27Rri8biOndH7GOxERLQjDocDHo8HFRUVWm1mZgaBQGDTWjzpi8FOREQ7ZrPZ0NnZifr6eq0WiUQQCATw8uVLHTujdQx2IqIsk+3noAsh0NjYiPb2dpjNZgBAMpnE9evXcffuXUgpde4wvzHYiYiySC6dg15ZWQmPx4OSkhKt9vjxYwwMDCAWi235vt7eXp4Wl0YMdiKiLJJr56Db7XZ4PB4cPXpUq83NzcHv90NVVR07y18MdiKiLJKL56BbrVZ0dHTg7NmzWm1lZQXBYBBTU1M6dpafGOxERFkkV89BF0Lg3Llz6Ojo0LacTaVSGBkZwdjYmDYLQenHDWqIiLKMoihZcQ76XpYBjh07Bo/Hg6GhISwtLQEAJiYmoKoq2tvbUVBQcMBd0vs4YiciogNVWloKj8eDyspKrfb27Vv4/X7Mz8/r2Fl+YLATEdGBs1gsaG9vR2Njo1ZbXV1Ff38/VldXdezM+BjsRESUNvX19ejq6oLVagWwtu6+tLSEpaUlrrunCYOdiIjSqqKiAl6vFw6HQ6tx9J4+DHYiIkq74uJiuN1uVFVVabX5+Xn4/X68fftWx86Mh8FOREQZYTab0dbWhpKSEu2o12g0iv7+fkxMTOjcnXEw2ImIKKOKiorgcDhgs9kAAFJKjI2NYWRkBMlkUufuch+DnYiIMs5qtcLr9cLpdGq1qakp9PX1YWVlRcfOch+DnYiIdFFUVISenh6cPHlSq6mqCr/fj7m5OR07y20MdiIi0o3JZMLFixfR3Nys7ZEfi8UwMDCAx48f69zdzmTbaXUMdiIi0l1NTQ1cLpe25ayUEnfv3sX169e57r5Luga7EOKXhRBSCHFYzz6IiChzVFXF5OTkl86aLy8vh9frRXl5uVZ7+fIlAoEAIpFIptvMWboFuxDiJID/HMCkXj0QEVFmhUIhjI6OYmJiApcvX/5SuBcWFsLlcqGmpkarLS4uIhAIYGZmJuumvbORniP2/xXArwCQOvZAREQZ5PP5tK1kY7HYB0+QM5lMaG5uRmtrq7buHo/HMTg4iOXlZUjJ2NiOLsEuhPgWgBdSylt6XJ+IiPTR29urhbXNZtt29F1dXY2enh4UFhZqteXlZSwuLiKRSKS71ZyVtvPYhRB/BeDYB370awD+KYC/tcPP+QzAZ8Da32QiIqPbyznoucLlcqGlpQWqquLKlSsfPXPe6XTC6/Xi+vXrePPmDYC1kX4gEEBHRwdKS0sz0XZOSduIXUr5VSnlhff/A+AJgFoAt4QQTwGcAHBDCPGhLwGQUn4hpWyXUrZXVFSkq10iIsoQRVFQXV390VBfV1BQgK6uLtTV1Wm1paUlBAIBTE9Pp6vNnJXxqXgp5W0p5REpZY2UsgbAFIA2KSX/7hAR0QeZTCY0NTXBbrdr+8wnEgkMDQ1hfHyc6+4b8Dl2IiLKGQUFBVAUBcXFxVrt4cOHuHbtGuLxuI6dZQ/dg/3dyJ17BxIR0Y5YLBZ4vV5sXJ6dmZlBIBDAwsKCjp1lB92DnYiIaLesVis6OztRX1+v1SKRCILBIF6+fKljZ/pL213xRERE6SSEQGNjI5xOp3bkazKZxPXr1xEOh9HY2Kitx6fTzZs3036N3eCInYiIclplZSU8Hg9KSkq02uPHjzEwMIBYLLbl+4y6ix2DnYiIcp7dbofH48HRo0e12tzcHPx+P1RV1bGzzGOwExGRIVitVnR0dODcuXNabWVlBcFgEFNTUzp2ltnZAQY7EREZhhACZ8+exaeffgqLZe02slQqhZGREYyNjWn71BsZg52IiAzn6NGj8Hq9sNvtWm1iYgKhUAjRaPRAr5VIJLC6uvqlk+r0wmAnIiJDKikpgdvtRmVlpVZ7+/Yt/H4/5ufnD+QaoVAIkUgE0Wj0g8fQ6oHBTkREOUNVVUxOTu44QC0WC9rb2zc9+ra6uor+/n6srq7uu5+NB/ZsdQxtpjHYiYgoJ4RCIYyOjmJiYmLXo+P6+np0dnbCarUCWFt3X1pawtLS0r7W3TfeEPexY2gzhcFOREQ5wefzaSG8l9FxRUUFvF4vHA6HVtvv6N3lcqGkpAQFBQW4evXqjk+sSycGOxER5YTe3l6YTGuxtdfRcXFxMdxuN6qqqrTa/Pw8/H6/dt77blksFhQWFmZFqAPcUpaIiDJsr+vQLpcLLS0tUFUVV65c2XOQms1mtLW1oaSkBMvLywCAaDSKUCiEpqYm1NbW7ulzswWDnYiIcoaiKFAU5UBGx0VFRbBYLLDZbIjFYpBSYmxsDOFwGC0tLTCbzQfQceZxKp6IiPKW1WqF1+uF0+nUalNTU+jr69NG87mGwU5ERHmtqKgIPT09qK6u1mqqqiIQCGB2dlbHzvaGwU5ERHnPZDKhtbUVLS0t2g16sVgMg4ODePz4sc7d7Q6DnYiI6J1Tp07B5XKhoKAAACClxN27d3H9+nUkEokPvodbyhIREWWBrXaxKy8vh9frRXl5uVZ7+fIlgsEgIpHIptdyS1kiIqIs8LFd7NafS6+pqdFqi4uLCAQCmJmZ0WrcUpaIiCgL7GQXO5PJhObmZly8eFFbd4/H4xgcHMSDBw8gpeSWskRERNlgN7vYnTx5Ej09PSgqKtJq9+/fx/DwMNrb23e0pexuD6/ZDwY7ERHlnfVd7Gpra3e0x7vT6YTH48GhQ4e02vT0NILBIMxm85Zbyvb29uLSpUt7PrxmLxjsRESUlxRFQXV19Y53sSsoKIDL5UJdXZ1WW1paQiKR2PaEOFVV93V4zW4x2ImIiHZICIGmpia0tbVt2nI2mUxifHwcUsovvUdRlH0fXrMbDHYiIqJdqqqqgtvtRnFxsVZ7+PAhrl27hng8vum1iqLsatp/vxjsREREe+BwOOD1enHu3DntRLiZmRn4/X4sLCxseu1up/33g8FORES0R1arFXa7fdPIfXl5GcFgEC9evNClJwY7ERHRPgghUFxcjPb2dlgsa6ehJ5NJ3LhxA5FI5IPr7unE89iJiChnZMPObluprKxEaWkphoaGtK1nV1ZWkEgkYLfbtRvo0o0jdiIiogNit9vh8Xhw9OhRrRaPx6Gq6paHyBw0BjsREdEBslqt6OjowLlz57RaMpmEqqp4/vx52q/PYCciIjpgQgicPXsWDodDm4KXUuLmzZu4ffv2thva7BeDnYiIKE1sNhsURdFuqgOAp0+fIhQKIRqNpuWaDHYiIqI0MpvNUBQFBQUFWu3t27fw+/2Yn58/8Osx2ImIiNJMCAG73Y7GxkYIIQAAq6ur6O/vx7Nnzw70Wgx2IiKiDKmvr0dnZyesVisAIJVKYXR0FLdu3TqwdXcGOxERUQZVVFTA6/XC4XBotcnJSfT19WF1dXXfn89gJyIi2gefz7frjXOKi4vhdrtRVVWl1cLhMPx+P968ebOvfrjzHBER5aVM7GKnqipUVYWiKFAUZdPPzGYz2tra4HQ6cffuXUgpEY1GEQqF0NTUpB0ss1sMdiIiojQIhUIYHR1FKpWCyWRCS0vLB19XV1cHRVEwPDyMWCwGKSXGxsYQDof3dF1OxRMREaWBz+fTbogTQuCnfuqntnztoUOH4PV64XQ6tdrU1NSerstgJyIiSoPe3l5t1zmbzYbe3t5tX19UVISenh5UV1fv67oMdiIiojRwuVxoaWlBbW0trl69CpfL9dH3mEwmtLa2oqWlZc+nwXGNnYiIKE3Wb5rbSahvdOrUqU2Pw+0GR+xERERZqKysbE/vY7ATEREZCIOdiIjIQBjsREREBsJgJyIiMhAGOxERkYEw2ImIiAyEwU5ERGQgDHYiIiIDYbATEREZCIOdiIjIQBjsREREBsJgJyIiMhAGOxERkYHw2FYiIqI08fl8Gb8mR+xEREQGwmAnIiIyEAY7ERGRgTDYiYiIDITBTkREZCC6BbsQ4h8LIe4LIe4IIf6FXn0QEREZiS6Puwkh/jMA3wbQIqWMCiGO6NEHERGR0eg1Yv8HAP65lDIKAFLKGZ36ICIiMhS9gv0sAI8QYlAI8f8JITp06oOIiMhQ0jYVL4T4KwDHPvCjX3t33TIAXQA6APyJEKJOSik/8DmfAfgMAKqrq9PVLhERkSGkLdillF/d6mdCiH8A4AfvgvyaECIF4DCA2Q98zhcAvgCA9vb2LwU/ERER/Q29puL/HMBXAEAIcRaADcCcTr0QEREZhl6HwPwBgD8QQowBiAH4uQ9NwxMREdHu6BLsUsoYgJ/R49pERERGxp3niIiIDITBTkREZCAMdiIiIgNhsBMRERkIg52IiMhARC49ZSaEWARwX+8+DO4wuKdAJvD3nH78Hacff8fpd05Kad/NG/R6jn2v7ksp2/VuwsiEEMP8Hacff8/px99x+vF3nH5CiOHdvodT8URERAbCYCciIjKQXAv2L/RuIA/wd5wZ/D2nH3/H6cffcfrt+necUzfPERER0fZybcRORERE28iJYBdCfF0IcV8I8UgI8at692NEQoiTQoi/FkLcE0LcEUL8kt49GZUQwiyEGBFC/Ee9ezEiIYRTCPGnQojxd/88u/TuyYiEEP/tu39XjAkh/kgIUah3T7lOCPEHQoiZdyefrtfKhRB/KYR4+O6/yz72OVkf7EIIM4B/DeBvAzgP4O8KIc7r25UhJQD8EyllI4AuAP+Qv+e0+SUA9/RuwsD+NwD/r5SyAUAr+Ls+cEKIKgD/DYB2KeUFAGYAP61vV4bwhwC+/l7tVwFclVKeAXD13Z+3lfXBDuBTAI+klE/eHff67wF8W+eeDEdK+UpKeePdXy9i7V+GVfp2ZTxCiBMAvgHg+3r3YkRCCAcAL4DfB9aOiJZShvXtyrAsAIqEEBYAxQBe6txPzpNS+gG8fa/8bQD/9t1f/1sAf+djn5MLwV4F4PmGP0+BgZNWQogaAJcADOrbiSH9SwC/AiCldyMGVQdgFsC/ebfc8X0hRIneTRmNlPIFgP8ZwCSAVwBUKeV/0rcrwzoqpXwFrA3AABz52BtyIdjFB2q8lT9NhBClAP4MwHeklAt692MkQohvApiRUl7XuxcDswBoA/C/SykvAYhgB1OXtDvv1nm/DaAWwHEAJUKIn9G3K1qXC8E+BeDkhj+fAKd80kIIYcVaqF+RUv5A734MqAfAt4QQT7G2pPQVIcS/07clw5kCMCWlXJ9t+lOsBT0drK8CmJBSzkop4wB+AKBb556M6rUQohIA3v33zMfekAvBPgTgjBCiVghhw9oNGj/UuSfDEUIIrK1L3pNS/o7e/RiRlPK7UsoTUsoarP1z/GMpJUc5B0hKOQ3guRDi3LvSZQB3dWzJqCYBdAkhit/9u+MyeJNiuvwQwM+9++ufA/AfPvaGrD8ERkqZEEL8IwA/wtqdl38gpbyjc1tG1APg7wG4LYS4+a72T6WUf6FjT0R78Y8BXHk3EHgC4Bd07sdwpJSDQog/BXADa0/UjIC70O2bEOKPAPQCOCyEmALwPwD45wD+RAjx97H2heq//OjncOc5IiIi48iFqXgiIiLaIQY7ERGRgTDYiYiIDITBTkREZCAMdiIiIgNhsBMRERkIg52IiMhAGOxEtCUhRIcQYlQIUSiEKHl3/vYFvfsioq1xgxoi2pYQ4nsACgEUYW0f9s91bomItsFg80d0VAAAAKZJREFUJ6JtvduadQjAKoBuKWVS55aIaBuciieijykHUPr/t3fHNAiFMRRGb8OMAEasIAQF6EIAAhhRRVKGt5K3/i/NOQq6fWk6NMk52+YOHJiNHdhVVa9sb2avSS7d/Vg8ErDj8N/dgHWq6p7k293Pqjol+VTVrbvfq2cD/rOxA8AgbuwAMIiwA8Agwg4Agwg7AAwi7AAwiLADwCDCDgCDCDsADPIDNmrkroPw7pAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Set seed for reproducibility.\n", "seed = 42\n", "np.random.seed(seed)\n", "\n", "# Choose the \"true\" parameters.\n", "mTrue = -0.9594\n", "bTrue = 4.294\n", "\n", "# Generate some synthetic data from the model.\n", "N = 50\n", "x = np.sort(10*np.random.rand(N))\n", "obserr = 0.5 # Amplitude of noise term\n", "obs = mTrue * x + bTrue # True model\n", "obs += obserr * np.random.randn(N) # Add some random noise\n", "\n", "# Now plot it to see what the data looks like\n", "fig, ax = plt.subplots(figsize=(8,8))\n", "\n", "ax.errorbar(x, obs, yerr=obserr, fmt=\".k\", capsize=0)\n", "x0 = np.linspace(0, 10, 500)\n", "ax.plot(x0, mTrue*x0+bTrue, \"k\", alpha=0.3, lw=3)\n", "ax.set_xlim(0, 10)\n", "ax.set_xlabel(\"x\")\n", "ax.set_ylabel(\"obs\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Inference**\n", "\n", "---\n", "\n", "Now we want to infer posterior probability distributions for our linear model's parameters, $\\theta$, i.e slope and intercept, given the data and uncertainties, D, via Bayes' Theorem: $p(\\theta | D) \\propto l(D | \\theta)p(\\theta)$ where $l(D|\\theta)$ is the likelihood of the data for a given set of model parameters, and $p(\\theta)$ is our assume prior probability of a given $\\theta$. We sample the posterior distribution using the emcee MCMC code. See the emcee example (https://emcee.readthedocs.io/en/latest/tutorials/line/) for more details of this procedure." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Define the loglikelihood function\n", "def logLikelihood(theta, x, obs, obserr):\n", " \n", " # Model parameters\n", " theta = np.array(theta)\n", " m, b = theta\n", " \n", " # Model predictions given parameters\n", " model = m * x + b\n", " \n", " # Likelihood of data given model parameters\n", " return -0.5*np.sum((obs-model)**2/obserr**2)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Define the logprior function\n", "def logPrior(theta):\n", " \n", " # Model parameters\n", " theta = np.array(theta)\n", " m, b = theta\n", " \n", " # Probability of model parameters: flat prior\n", " if -5.0 < m < 0.5 and 0.0 < b < 10.0:\n", " return 0.0\n", " return -np.inf" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Define logprobability function: l(D|theta) * p(theta)\n", "# Note: use this for emcee, not approxposterior!\n", "def logProbability(theta, x, obs, obserr):\n", " \n", " lp = logPrior(theta)\n", " if not np.isfinite(lp):\n", " return -np.inf\n", " return lp + logLikelihood(theta, x, obs, obserr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we've set up the required functions, we initialize our MCMC sampler in emcee, pick an initial state for the walkers, and run the MCMC!" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/anaconda3/lib/python3.7/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in double_scalars\n", " lnpdiff = f + nlp - state.log_prob[j]\n" ] } ], "source": [ "# Randomly initialize walkers\n", "p0 = np.random.randn(32, 2)\n", "nwalkers, ndim = p0.shape\n", "\n", "# Set up MCMC sample object - give it the logprobability function\n", "sampler = emcee.EnsembleSampler(nwalkers, ndim, logProbability, args=(x, obs, obserr))\n", "\n", "# Run the MCMC for 5000 iteratios\n", "sampler.run_mcmc(p0, 5000);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The MCMC is complete, so let's examine the joint and marginal posterior probability distributions it derived for the model parameters." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAF+CAYAAABNtpsqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXzU9bX4/9eZmUySyQoZEpJAEiAsYkAERCgiKKgFwV+V0taFatWL1ba33FattLUuXfT+tLZXb+3Vtiqudde6YJVaFtlBwk5YEwIh+zrZJ/P+/jGTMUCCEQYyMznPx2MeZvJ5z3tOYjh553zeixhjUEopFVosPR2AUkqpr06Tt1JKhSBN3kopFYI0eSulVAjS5K2UUiFIk7dSSoUgTd6qx4nIb0VkpYi8ISKOL7suIikislpElovIpyKS2hNxK9WTNHmrs0pElh33PAcYYoyZAiwFbu7G9XLgImPMVOB54JazELpSQUWTt+ppU4Alvo+XABd92XVjTJsxxuP7XByw44xHqVSQ0eTdi4lIXxF5W0TqRaRARK47SdtzfCWKGhHZJyJXH3d9mYg0iYjL98jrZhh9gBrfxzVA3+5cF5ExIrIO+CHweTffS6mwocm7d/sT0AKkANcDfxaRc49vJCI24F3gfbzJcwHwoogMO67pD40xsb7H8A6vz/Al92XAmPaPRSQDqAISfE0TgMrj+uz0ujEm1xhzIXAvsOjUvnylQpcm7zNMRPJF5C4R2eob4f7Nd8NtiYjUichSEenTA3HFAHOBe40xLmPMZ8A/gPmdNB8BpAF/8JUsPgVWddH2BMaYQ8aYacaYaUBu+8fGmEPAZ8AVvqZX+Prt6ITrIhLZ4XoN0NCdOJQKJ5q8z465wGXAMGAO3trtzwEn3v8H/9nZi0TkfRGp7uLx/mnGNAxoM8bs6fC5LcAJI29AuvhcznGfe0hEykVklYhM604QxphtQIGIrMSbnJ8Rkf4i8kBX14GxIrJCRP4NLAQe6c57KRVObD0dQC/xhDGmBMCXhEqNMZt9z98Gpnf2ImPM7DMYUyxf1JLb1eC9AXi83UApcJeI/AG4BJgK/LtDm58BO/GWYb4DvCciY4wx+zt25Bt9c9znji971AP3neT6GuDiTr8qpXoJHXmfHSUdPm7s5HlsIN/MV082XTw+8zVzAfHHvTQeqDu+P2NMK/AN4EqgGPgp8BpwuEObdcaYOmNMszFmMd7yx6xAfl1KqS9o8g5ivrq4q4vHkq5e56snSxeP9ql4ewCbiAzt8NLz6GLanTFmqzFmqjEmyRhzBTAYWH+S8A2dl1s6+zq/6iKdSSKyxrdI5xURiejO+ygVTjR5BzFjzMwOszeOf8w8zb7rgbeAB0UkRkQmA/8f8EJn7UVktIhE+ZLnnUAq8JzvWqKIXOG7bhOR6/GWNf75ZXGc4iKdAuBS3yKdA764lepVNHn3bncA0Xjr2a8Atxtj/CNv38j/576n84GjvrbTgcuMMc2+axHAb4AyvKsffwR8wxjTnbnep7JIp8gY0+j7nBvwoFQvI3oMmupJvl8OO40x74hIIvCaMeby7lwXkUHAq3gTektPxK9UT9HZJuqME5H+wBudXPomp7hIR0TigcXAfE3cqjfSsok644wxxcaYizp5FHNqi3RseMs893ezNKNU2NHkrXrUKS7SuRa4EPiVb1rkt3sofKV6jNa8lVIqBOnIWymlQpAmb6WUCkG9braJ0+k0WVlZPR3GWbF/v3dbkSFDhnR6/XBFPQADkmK61f6r9h+ONm3aVG6M6dfTcSjVK5K3iCzAuwc1GRkZbNy4sYcjCg53Pb8GgEe+O6mHIwkdIlLQ0zEoBb2kbGKMedoYM94YM75fPx00KaVCX69I3r3VokWLWLSo+4fMnOn2SqnA6RVlk95qzZo1QdVeKRU4OvJWSqkQpMlbKaVCkCZvpZQKQVrzDmMDBgwIqvZKqcDR5B3GXnzxxaBqr5QKHC2bKKVUCNLkHcYWLlzIwoULg6a9UipwtGwSxnJzc4OqvVIqcDR590JZWVkUFBQw/nsPASA3fs1/LTIysqfCUkp9BVo26YUKCgowxjBt2lSmTZuKMQZjDFOnTqW5ufnLO1BK9ThN3kopFYK0bBLGhg0b9pXbb9q0CRHpsk1mZib5+fmn1L9SKnB63RmW48ePN719P28RwRhzSvt5t7+2txKRTcaY8T0dh1JaNlFKqRCkyTuMLViwgAULFgRNe6VU4GjNO4zt2bMnqNorpQJHR95KKRWCNHkrpVQI0uStlFIhSGveYWzMmDFB1V4pFTg6z7sX0nnep07neatgoWUTpZQKQZq8w9gNN9zADTfcEDTtlVKBozXvMHb48OGgaq+UChwdeSulVAjS5K2UUiFIk7dSSoUgrXmHsUmTuj8F8Gy0V0oFjs7z7oV0nvep03neKlho2UQppUKQJu8wNnfuXObOnRs07ZVSgaM17zBWUVERVO2VUoHTK0beIrJARDaKyMaysrKeDkcppU5br0jexpinjTHjjTHj+/Xr19PhKKXUaesVyVsppcKN1rzD2PTp04OqvVIqcHSedy+k87xPnc7zVsFCyyZKKRWCNHmHsZkzZzJz5sygaa+UChyteYexxsbGoGqvlAocHXkrpVQI0uStlFIhSJO3UkqFIK15h7HZs2cHVXulVODoPO9eSOd5nzqd562ChZZNwlRWVhYi0ukjMzOzp8NTSp0mLZuEqYKCAqZOnQrAsmXLuvWaadOmndH2SqnA0ZG3UkqFIE3eSikVgjR5K6VUCNLkrZRSIUhvWIaxb33rW0HVXikVOJq8w9gdd9wR8PaZmZmIyDGf+8EPfuC/lp+f/5XeUyl1ajR5h7GGhgYAHA5HwNp3TM7Htz8+qSulzhxN3mFs1qxZQPfnYZ/p9kqpwNEblkopFYI0eSulVAjS5K2UUiFIk7dSSoUgvWEZxm666aagaq+UChxN3mFMk7dS4UvLJmGsvLyc8vLyoGmvlAocHXmHsW9+85tA9+dhn+n2SqnA0ZG3UkqFIE3eSikVgjR5K6VUCNLkrZRSIUhvWIax22+/PajaK6UCp1ckbxFZACwAyMjI6OFozp5vf/vbQdVeKRU4vaJsYox52hgz3hgzvl+/fj0dzllTWFhIYWFh0LRXSgVOrxh591bz588Huj8P+0y3V0oFTq8YeSulVLjR5K2UUiFIk7dSSoUgTd5KKRWC9IZlGPvpT38aVO2VUoGjyTuMzZkzJ6jaK6UCR8smYSwvL4+8vLygaa+UChwdeYex2267Dej+POwz3V4pFTg68lZKqRCkyVsppUKQJm+llApBmryVUioE6Q3LMPbLX/4yqNorpQJHk3cYmzFjRlC1V0oFjpZNwlhubi65ublB014pFTg68g5jCxcuBLo/D/tMt1dKBY6OvJVSKgRp8lZKqRCkyVsppUKQJm+llApBesMyjP3ud78LqvZKqcDR5B3Gvva1rwVVe6VU4GjZJIRlZWUhIp0+MjMzWb16NatXr+52f2e6vVIqcHTkHcIKCgowxnR5fdq0aUD352H//Oc/P6PtlVKBoyNvpZQKQZq8VcBkZmZ2WcbJysrq6fCUCitaNlEBk5+f3+U1ETl7gSjVC+jIWymlQpCOvMPYH//4x6Bqr5QKHE3eYWzMmDFB1V4pFThaNgljS5cuZenSpUHTXikVODryDmO/+c1vgO6feHOm2yulAkdH3kopFYJ6RfIWkQUislFENpaVlfV0OEopddp6RfI2xjxtjBlvjBnfr1+/ng5HKaVOW69I3kopFW70hmUYe+qpp4KqvVIqcDR5h7Hhw4cHVXulVOBo2SSMvffee7z33ntB014pFTg68g5jv//97wGYM2dOULRXSgWOjryVUioEafJWSqkQpMlbKaVCkCZvpZQKQXrDMoy98MILQdVeKRU4mrzD2MCBA4OqvVIqcLRsEsZeffVVXn311aBpr5QKHB15h7E///nPAHz7298OivZKqcDRkbdSSoUgTd5KKRWCNHkrpVQI0uStlFIhSG9YBrmsrCwKCgo6vZaZmXnS177xxhtf6b3OdHulVOBo8g5yBQUFGGNO6bVOpzOo2iulAkfLJmHsueee47nnngua9kqpwNHkHcaCKXlnZmYiIp0+srKyuv2eSikvLZuosyI/P7/LayJy9gJRYUdEfgtcDJQA3zXGNHS4lgB8AowEJhpjtotICvA20Aq0AdcbY46e/chPj468lVIhQ0SWHfc8BxhijJkCLAVuPu4lDcCVQMe76+XARcaYqcDzwC1nLOAzSJO3UmFIRPJFZEZPx3EWTAGW+D5eAlzU8aIxptUYU3bc59qMMR7f0zhgxxmP8gzQ5K2UOi0iMlREmkTkxS6uR4rI30SkQETqRGSziMw8rs0yXx8u3yOvm2/fB6jxfVwD9O1mzGNEZB3wQ+Dzbr5XUNGadxj78MMPg6q9Clt/Ajac5LoNKASmAoeAWcBrIjLKGJPfod0PjTF/Pf7FIpKBt7wBMKZD6eS7QBWQ4HueAFR2J2BjTC5woYh8C1gEfL87rwsmOvIOYw6HA4fDETTt1Vl3gYjsFJEqEXlWRKIC/QYi8h2gGvhXV22MMfXGmPuNMfnGGI8x5n3gIDCuO+9hjDlkjJlmjJkG5LZ/bIw5BHwGXOFregWwqhsxR3Z4WoO3Lh5yNHmHsSeffJInn3wyaNqrs+56vAltCDAM+GVnjUTkfRGp7uLxfledi0g88CDw068SlG+2xzBOrDU/JCLlIrJKRKZ1py9jzDagQERW4v1anxGR/iLyQIf3+xC4HPiLiNwEjBWRFSLyb2Ah8MhXiT9YaNkkjL322msA3HHHHUHRXp11/2uMKQT/dLon6CSBG2Nmn2L/vwb+Zowp7O50TxGJAF4CFhtjdne49DNgJ9ACfAd4T0TGGGP2HxfrtOP7NMYsOu5T9cB9Ha7P6iSUi7sVcBDTkbdS4auww8cFQFqgOhaRMcAM4A9f4TUW4AW8CfqHHa8ZY9YZY+qMMc3GmMV4yx+dJV3lo8m7h2VlZXW58lBEvnTzqXCgqy/PmI6HjGYARZ01EpElHWZ5HP9Y0tlrgGlAFnBIRIqBO4G5ItLpzA3xDs3/BqQAc40xrV8SuwG6NZwXkd+KyEoReUNEHMddSxCR9b6vJcf3uRQRWS0iy0XkUxFJ7c77BBstm/Sw09l4Klzo6ssz5ge+mnUD8HOg0wNHjTEzO/v8l3ga+HuH53fiTea3d9H+z8A5wAxjTGPHCyKSCFwILAfcwLfxljUWflkQHRfpiMj38S7S+d8OTdoX6XSsa7cv0vH4auC3AL/5svcKNjryVip8vQx8DBzwPQKWoIwxDcaY4vYH4AKa2hfE+EbzP/d9nAncBowBijuM6q/3dRfhi60Mb2L9EfANY0x35nr32kU60htGfSKyAFjgezoc6O4CgEBx4v2hDEbBGluwxjXcGBPX00EoL98viJ3GmHd8I/jXjDGXd9LuOeBRY8x23/MxwFNAInC5MabzTfODWK8omxhjnsb7Z16PEJGNxpjxPfX+JxOssQVzXD0dQ28jIv05dm+Sdt+kFy/S6RXJWykVunxlmYs6uyYin+FNvov5Cot0jDHNvqchu0hHk7dSKmQZY7b59kxZCZQC3/WN1G83xtwH/kU6Y4DhIvIUkCci/413O9gmTtyJMCRo8j47eqxk0w3BGpvGpbqlty7S6RU3LJVSKtzoVEGllApBmryVUioE9bqat9PpNLrk2utwRT0AA5JiejiS0LFp06ZyY0y/rq6H48/X/v3evaGGDBkS8H6+6s9goGIJVl/289VRr0veWVlZbNzYs1N1PR4PDQ0NOBwOLJae++PnrufXAPDIdyf1WAyhRkROupgjGH6+Qon+DB7ry36+OtKySQ9oaGigtraWhoaQnF6qlAoCIZu8JYR3LHI4HMTHx2Oz2di/fz8tLS09HZJSXVq0aBGLFh0/G69n+glULOEg5MomInKuMWaHCYM5joWFhezatYvdu3czffp0oqICfkqVUqdtzZo1QdNPoGIJByE18haRK4CXRGRoT8cC4HK5+PTTT3G5XF/pde1lk6SkJKxWK2VlZezZs+ekr/F4PLhcLjwez0nbKaV6h5AZeYvIVcA9wA+MMXtFRLo7+u64q2BGRkbAYlq/fr1/JHDppZcec+1kNyXbD+11OBxMnz6dPXv2MGzYsJO+V3vCB4iNjQ3Ul6CUClFBn7x9RyeBd7/fGmPMKt8BpvNFxIl3g/l9xpi6rvrouKvg+PHjA1ZumTBhwjH/befxeCgtLcXtdgMnJluLxeL/XFRUFDk5OTQ0NODxeLqcfdIx4SulVNAnbyDZGFMsIpOBNSLyBt6tH/8F9AF+DPyFbuwmFmixsbEnjLjBO0p2u93YbLZuJdvujKo7JnylzqYBAwYETT+BiiUcBHXyFpGvA++IyNXGmCUiciGQC7xsjHnY1+Yh4Lv0QPLuis1mo7KykmHDhh0zkna73ZSXl+N0OrHZvvjWOxwOPB6P/9H+muNLL01NTf4Si97cVGfLiy++GDT9BCqWcBC0Nyx9ifs3wNvAKBGJMsbU4z0H7/4OUwX3AS4RsfZQqCc4cuQIR48epaCggOLiYn/5pLy8nKKiIkpLS4+5+WixWLBYLJSWlrJs2TL//O/j54Pv2bOH7du3n3BzU29mKtX7BOXIW0SmAg/jPd2iFXgSeAbvGXnuDu1uxXvg6XeNMW09EWtnBg70HtodERFBUVERTU1N1NbWMnjwYNLS0oiKiuq0TLJkyRL/SHvixInH1LlbWlqwWq0MGzbshJubejNTnUkLF3rPAf7jH//Y4/0EKpZwEJTJG+8p1AuMMesBRGQn8JiI3GyMcftG2cOAq4CbjDFBdYCo3W5nyJAhuN1u7HY7RUVF5ObmsmfPHmbPno3dbj+hHp6bm0tFRQUpKSmMHj0a8I7I3W43b731FoMHD6asrIzs7OwTSiZ6M1OdSbm5uUHTT6BiCQdBWTYxxiw2xqwXkfZfLk/gHYGn+a63AfuB64wx23oozC9ls9no378/I0aMICIigpqaGvbs2eO/+dixHp6SkkJCQgKTJk06JjkvXbqUTz/9lN27d5OdnU16evoJJZLO+lMqGK1duxYR8T+WLVvOsmXLERHCbUOvMy3oRt4iYm0vgXQokewC0oE78M71xhjTAoTEunKPx8PkyZMpLS1l2LBhtLS0sGPHDg4ePEhOTg6DBw8mOjqajIwMiouLycrKIiHBe6bqjBkz8Hg8XHDBBbjdbhoaGmhqagK0RKKCT1ZWFgUFXe+tFBkZScflGe0bU214xhDCO170iKAYqonIRSIyH7yj6g5zuxERizGmEfgRMENEJnTVTzDpeBPR4XCQnJzM2LFjKS8vZ/369bz55pssXryYv//972zfvh2n00lmZiZlZWW89dZbVFVVUVJSQmxsLF//+tfZv38/27Zto6Kigvj4eH8dXPdGUcGkoKAAY0ynj6lTpzJx4sSeDjFs9OjI25ekHcBT3qcSY4z5P2OMR0Rsxhi372MLUAksAQ72ZMzd1X4TsX3qn8PhoLS0lG3btlFeXk7//v0pLi4mOjqaPXv24Ha7iY6OZvfu3ezatYuioiIuu+wywFt+aW1tZcCAAWRmZhIREUFbWxuff/45paWllJSUMHLkSKKiojhy5AgDBw7Ebrf38HdAhYsvW/17NvsJVCzhoEeTtzHGg3ea32K8Jzl/TUSijTF/6DirxNeuQkR+5xuFB732m4cej4fy8nJqa2vJzs4mKSmJoqIiWltbsVgsJCYm0tzcTF5eHtXV1QwcOJCDBw/S1tZGZWUlBw8eJDs7G5fLRUREBB6Ph/r6ev/skra2Nvbt20dERASRkZEcPXoUCN/N6oPRmdp+IVg8/XRgzlwORD+BiiUcBEvN2w1kAIuBW0XkMaDZGLNIRCb5Pv4caOrJIL+K9puIbrebnTt3UlFRQUxMDAkJCWzdupUjR45QW1uL1WplypQp1NfXU1BQwODBg7nyyiuJiopi7dq17N+/n4yMDCZPnkx9fT2ff/45SUlJ9OnTB4fDQWZmpn9/lWHDhhETE+OfqqjOjjO1/YJSJxMUNW/gXaDYGPMvYCPeuduJvmsTgCKAUNsGtqWlhc2bN1NcXExycjIDBw5kz549HDhwgNraWtLT09mwYQNvvvkmO3fuZPTo0Vx11VVcccUV2O127Ha7f5FORUUF0dHRREVFUVBQwFtvvUV1dTUej4f09HRiY2P9UxS1ZKICacGCBSxYsCAo+glULOEgWEbejcBwEfkPvAtzHgImiMh1wOOhlrTbFRYWsmHDBmpqasjKyiIiIoKkpCQsFou/nFJSUsLBgwfJysoiPj7evwpz8uTJpKenc+WVV1JXV0dLSwvFxcXs3r2bgwcPsnfvXhobG5k8eTIZGRns27cPl8vFeeedp8lbBdSXbVd8NvsJVCzhICiStzGmSEQKgXvxbvn6nohcgne3wJBI3J1tATtw4ECmTp1KXV0diYmJtLW1MXbsWG6//XZ27dpFXV0d+fn5/hLI2rVrOXjwIDabDZvNxpgxY9i0aRN9+vShoaGB6upqSktLaWhooLGxkdLSUnbu3InVauXQoUPs3buXxMREsrOzddqVUmEuKJK3z1+Ad40xm3zPl/tuVIaEzpao2+12zj33XIqLiykqKiIqKoqUlBQuv/xy4uPjeeihh2hsbKR///6kpKSwceNG3G432dnZDBgwgNdff53169eTmJhIdHQ05513HlFRUURGRnLw4EGqq6uJiori3HPPRURobW0lLi6uJ78NSqmzJGiStzGmEChsP2QhlBI3nHyJutPpBCAxMZGjR48SFRVFcXExVVVVtLa2kpaWRkxMDC6Xy3+yzkcffcTRo0ex2+20trZSVlaGxWJh0qRJrFixgsrKSg4dOsS4ceMA7w1SYwx1dXWkpKScvS9cKdUjgiZ5twuVMsnxTrbfttVqJSUlhZKSEg4cOIDL5eLw4cPMmDGDw4cPY7fbycjIYOfOndTW1lJTU8Nnn32GxWJhwIABREVFUVJSQlVVFUeOHKG+vp64uDicTqe/v5iYGEREE7cKuDFjxgRNP4GKJRwEXfIOZ06nE4/Hw7Jlyzhw4AAXX3wxsbGxrFixgpqaGsaOHUtzczMHDhwgLy+P5uZmUlJSiIyMpKGhgfj4eA4fPkxjYyMOh4PCwkJeffVV6uvr6devH+eccw5bt25l8ODBxMfHc+DAgS73/vZ4PLS1tWGxBM1OuipIBWoHv0D0o7sJfkGT91lktVpJTU1lxowZlJeXM2LECP+qSoBzzjmHffv2+VdexsbGUlxcTG1tLSkpKbS1tflH4A0NDf59w51OJ3FxcTQ0NFBeXk5LSwvl5eX+qYSdjVa8p/20YdOfAKVCkv7TPYuMMTQ0NLBv3z7y8/N55JFHcDqdjBw5ktLSUnJzc1m2bBmFhYUA5OTkUFlZSVVVFY2NjTQ1NeF2u+nTpw9ut5v6+nr27dtHW1sbqamp7Nixg9bWViIjIykuLqa+vr7TY9rAW5u32aw68lZf6oYbbgBO/xSbQPQTqFjCgSbvs8QY4z9Fp6Kign379lFRUUFzczNHjhyhpqaG8vJyGhsbiYyMBEBE/HueGGOIjIwkMjKSlpYW+vXrR3NzMwD19fUcPHiQ+vp6ioqK6NOnDxs3bsThcJCQkMCiRYuwWCxER0f7pzFaLBasVk3c6ssdPnw4aPoJVCzhQJP3GdS+s2C70tJSysvLyczMZObMmRw5coSysjK2bNkCQHR0NElJSbhcLlpbWzl69Kg/QR+v/SZlXV0d9fX17N27F/CWQ1JTU6mvr6e6upoNGzbwwgsvMG3aNAYOHKhTCZUKE5q8z6CGhgZKSkoA72ELgwcPJjk5GYfDwbx586ioqMDj8fD73/+ePXv2UFtbi9vtpq2teye61dXVHfPcZrMRFxdHamoqjY2NlJWVUV1dzdKlS2lqauLcc8/l4osvJjo6GrfbTUtLCxEREYD3cOTDhw9TWVnp36FQKRW8gmVvk7DkcDj8J+QcOHDAf/yZ92ah23/w8Pz588nJySE2NpaSkpJTPki4T58+1NfX4/F4GDFiBOeddx6RkZGsWbOGDz/8kE8++YSNGzdSX19PcXExDQ2NNDY2+Zfqr127lrVr156wBNnlcvHpp58e81eEUoGWmZl5zCk7HR96ys6JdOR9BlksFuLj48nPz2f79u1UVlaSnZ2N3W7H6XRis9nYtm0bGzduJDY2lurqahobG/2v/arKysoAWLduHXFxccTFxfluTNro06cP/fv3Jy4uzr8SVET8N1GdTicTJ06ksrLyhD2T169f79+5sKsboCp8TZo06az0k5+f3+W19u0eAhVLONDkfRYMGzbMnyDtdjvJyclYLBZqamrYvHkzRUVF2O12IiMj/T+kp3Mzsb1kEh8fj9PpJDU1lfT0dKqqqqipqaFfv34kJycTE5MP4N+PJSsrq9MRzoQJE475r+pdHnrooaDpJ1CxhANN3mdBVFQU48ePp7y8HKfTicViweVysXXrVpqbm4mPj6dv37706dOH9PR0ysrKaG1tPe33bZ8XXlxcTF1dHXFxcf6/BhITE7HZvL8gvmyUHxsbqyNupYKMJu+zpKnJW1uurKxk3759HD58mDfffJO2tjZiYmIoKiqiqamJxMRESktLsVqtpKWl+csoxxs5cuRJD3otKirCarXS1taG3W73b27V1tZGVVWVfy+UL1ZafpHAdUdC1dHcuXMBePPNN3u8n0DFEg5CNnm3b2DV03F0V/uGVVu3bmXNmjX+XQLbZ5fs3LmTffv2ERsbi8fj8a+SdLvdnfZXUlLin23SnoQ7llra2tr8s1ZaW1txuVw0NDSwbt06oqOjiY+Pp6mpCRELjY2N2O12CgsL9fxLdYKKiopjnp/shPjMzMxu9xOIWHqzkEveIhJhjGkNpcQNX2xclZOTQ1lZmf+mYHFxsf9U+JaWFhoaGkhOTu72YoTm5mYGDRrEsmXL6NOnD06n84SRs9vtpry8HIDa2lr//imtI66isbGB6upsmpub2bFjB7m5ucycObPT3RGVgi9OiFc9K+fEaX4AACAASURBVKSmCorIVcDjIrJYRM4Vkb49HdNX1dDQQEJCAkeOHCEmJoakpCSioqIYOnQokydPZtSoUURGRvrnX3fFGENZWRmFhYVs3ryZWbNmUVlZSXx8PPX19V3+42pqaiI/P5/33nuP4uKjVFfX8Mtf/pLm5mYaGxspLCxk69atZ+JLV0oFUMgkbxHJAf4MvA6UAD8Avisi6d147QIR2SgiG9un0/UUp9PpP1Vn8+bNuFwutm/fztGjRznvvPM4//zziY+PJz4+/qT91NbWUllZybx589iyZQuvvfYaL730Evn5+f5tZl0uV5dll7q6OqKjHVitFg4dOsTf//53Bg0aRHp6OiNGjDgTX7pSKoBCJnkDKcAKY8ynxpi7gQ+BAcBcEUk82QuNMU8bY8YbY8b369fvbMTaJZvNRk5ODm63m/z8fIqLi0lNTaVfv37YbDZEhMjISOx2e5f7g4O35h0bG8vjjz/uP+zhG9/4Bnl5eTz//PPk5OTQ2NhIRUUF2dnZJyTxxMREWltbsNki/PX1vXv3YrPZtK6ojjF9+nSmT58eFP0EKpZwEPQ1bxGJMsY04T1V/gERudoY87Yx5n0RsQCzASdQ3aOBfgU2m41LLrmEpKQkMjIySE9PZ+TIkWzZsoX+/ftTW1tLVVUVMTExnb7e4/FQX1/PvHnzTmgTFRXF3Llzue666zhy5AgvvfQSjz/+OMYYWltb/eWY2tpaGhubiIy0k5KSQmZmpn8b2pSUFFwu1zHncare69577w2afgIVSzgI6uQtIpcC54jIM8aYGhFZDHxNRCqNMcuNMf8QkRnAQuCHPRvtV+NwOJg4caJ/mt6qVav497//jcfjYejQoVRUVDBw4ED/DoMdVVRUUF9fz5w5cygtLe20f6vViohwww03MGXKFK677jqqqqoYOnQoQ4cOZdmyZRhjaGlppbq6mry8PCIjI8nJyeHAgQMkJiZy5MgRtmzZwuWXX05CQgKg0wg7IyILgAUAGRkZPRyN6i2CNnmLyNeB/x/vafLtk50/xVs+uUpEUowxrwFb8CZ4S6idewneWSjJyckMGzaM8vJy/9J1q9VKQUGB/3lHgwYNwm63M2HChC6Tacf9UTIzM3njjTe47rrr2LJlC8YY2tra8Hg8eDweaioqaGtr808ZrK2tJTo6mh07drBhwwbcbjfXXnvtmfkGhAFjzNPA0wDjx48Pu2kYM2fOBGDJkiU93k+gYgkHQZm8RWQ08ArwH8aYlSLiBASoAR4BrgfuFpFvAxOAK0Mtcbe0tPjnVbfXwe+8807WrFnjr4d3Nt/a7XazevVqFi1a9JVGwampqbz++uvccMMN5ObmHlNuaWxspKqqyn/48caNG0lLS2PEiBFkZ2frDcxerquFYj3RT6BiCQdBmbyBKOA1oL+IjAd+A5QBE4F7jTHPiMjrwEig0BhT1HOhnprCwkL27duHx+PxL1u/4ooruPLKK3n11VfZu3cvFRUV2Gy2Y242NjQ0YLVaueaaazrt98CBA6xdu5bKykqqq6uprq6mpaWFBQsWMHLkSP7+979z4403snXrVjweDxaLxb/Yp7y8nPz8fFwuF1u2bGHq1KkBTd5NTU3s2bOny3M1lVLdF5TJ2xizXkQigHl4E/fP8f5ZOgl4RUR2GWO2AOt6MMzTMnDgQMC7javL5aKkpIQtW7YwYMAAIiMjSUlJoays7IRZIueeey42m43k5OQT+jTG8L3vfY/9+/cD+JfFV1ZW8uGHH/LrX/+ab33rW7z00kt8//vfp8Htxmr1LpMXEVpaWqirq+PQoUPYbDZWrFjBpZdeypEjRxgyZMhp17v37NnD9u3bARg9evRp9aVUbxd0UwlExAZgjFkFvArcYox5EmgzxqwEPgLqezDE09a+/H3QoEEkJiZit9tZtWoV7777Ls8//7x/i9isrKwTTr7ZunUrxcXFnfa7ZcsW9u/fz/3338/u3bvZu3cvGzduZO3atYwdO5a77rqLhQsXIiJ8+OGHpKen0dbmYdy4cRhjqKmpYdu2bVitVhwOB/369ePw4cN88sknXb7nVzFs2DBycnJO2HJWKfXVBcXIW0QuAgYZY14wxrhFxGaMcRtj1oiIHcAYY3w17vOBkC58tY+0U1JSiIuLIzk5mauuuorGxkYyMjI4dOgQpaWlxMfHn7CHhM1m4+DBgzQ1NZ0wTfCtt97Cbrczb968Y5a3Jycn89JLL/HEE0/wxz/+ka1bt/Laa6+RnZ2Nw+HglWcXER0dTUJCAi6Xi7KyMubNmwfAhg0baGtrIz09nenTp5/W9MGoqCgdcYeg2bNnB00/gYolHPRo8vbN03YAT3mfSowx5v98Cbx9D5MW32j8OuBO4FpjzJGejDvQRASHw8Fll13GkSNH2L9/P3v27CE2NvaEI9FsNhvGGA4ePEhOTo7/862trbz77rvMmDGDhISEE07jsVqtLFy4kAsvvJAf/vCHTJ48mat/9SKpqWn84he/4Le//S1Wq5X4+HgaGxspLy+npKQEl8vFgAEDOOecc/jwww+56KKLiI+P73VzwDvuD9Pb3HnnnUHTT6BiCQc9mrx9M0RcvvnbbXjncEcbY/5gjGnt0M4tInXANcaYfT0Vb6DExsZisVhwOBz+OrLT6cTj8WC3exfNlJaWUl9fT0NDwzHHj7Un8/379x+TvFeuXElFRYV/y8wPP/ywy/e/8847eeedd/w3RR/45mx27drFW2+9RVZWFk1NTRw4cIDS0lISExO59tprOXToEPv27SM6Oppx48YB3vLPypUrmTJlCnFxcWE9B7y8vJyioiIIkr9WlQqWH0Q3kAEsBm4VkceAZmPMIhH5GlBljHm7RyMMoPYdBjuy2Wz+bVqzs7PJz89n8ODBtLW1UVpa6t9l0GKxICLs23fs77C33nqLxMRELrnkki99/7i4OP7v//6PX72xhaIjR7jjjjt48sknKS8v57PPPsPhcHDo0CEA/43MhQsXMmHCBEaOHElsbCwOh4MlS5awYsUKAGbNmhWIb03Qat+CAO/Paq8ybdo0AJYtW9bj/QQqlnAQLH/3vgsUG2P+hXcZ/O1A+34lFwBVPRXY2eRwODDGUFBQ4K97p6SkMGjQINr3ZBERBg0axN69e/2vc7lcfPTRR8yZM6fbe3FbLBbSUlMZPHgIubm5/OIXv+Dhhx/277vSfpKPMYbDhw/zj3/8g127drF582aqq707EUyZMoWLL76YKVOmBPg7EXxsNhv9+/fv6TCU8guWkXcjMFxE/gP4PvAQMEFErgMeD7W9u0+VxWIhMzOT+Ph42traSEpKIiYmhpiYGPLy8vztCgoKaGtr80/xW7duHU1NTf7VZ19F3759+NnPfsbDDz/Ms88+y5/+9CduueUWKnyrLi0WC8XFxSxfvpyCggL69u3L1q1bufHGG0lNTQ37EbdSwSooRt6+RTaFwL3AT4wxDwK/B1b2hsTt8XhwuVx4PB5sNhsTJ05k4sSJXHTRRcyYMYOqqqpj9ue22WwUFBSQm5sLfLGfhq8m+5Vdf/31XH311fztb3/j4MGDvP3229TX1zN58mT/+5aUlFBfX4/FYmHnzp0888wzrFu3jubmZpqbm9m2bRtNTU2n+Z1QSnVXsIy8Af4CvGuM2eR7vjzUlryfqoaGBv8eJrGxsYwdO5bY2FiGDBnCunXriImJITU1lfa9yCMiIhARXn/9dc4//3yys7Pp06cPGzZs4Dvf+c4pxXDPPfewceNG7rnnHrZv387vfvc7fvKTnxATE+PfHCsmJob09HQ2b94MQH19PQcOHCAqKsq/SlOnAip1dgRN8jbGFAKF7WdT9pbEDV+cb+lwOHC73VRXVzNy5EisVis5OTmsWrXqmCmDIkJ8fDzvvvsuv/rVr4iKimLChAmsW3fqC04dDgcPPfQQN954Iz/+8Y955plneP/991m9ejVutxubzcbGjRs5evQoANHR0SQlJbFt2zaSkpI455xzjll84/F4aGho6HVTCsPRt771raDpJ1CxhIOgSd7tekOZ5HgdZ5+Ulpb6yx9OpxObzcbs2bOprKzk8OHD1Nd7F5eOHj2apUuX8sEHHzBr1izOP/98/vnPf3Lo0CH69evnL6l0xmazsXv3booSxgPwl7+85782bdo0XnjhBSZPnsz//M//MGvWLFpaWhg6dChxcXEsW7bMP0VxxIgRVFRUkJCQwCWXXHLMfiXH/zWhQtcdd9wRNP0EKpZwoEOiION0OklLS8PpdFJbW8uKFSvYv38/xhiam5v9+48UFxczYMAAPvjgA/r378/ll18OeOd/9+/fH6fT2eVj8+bN7Nixg5rqasrLytiwYQNFRUUUFRUxcuRIxo4dyz333ENzczOLFy+mvr6empoaqqur8Xg8HDlyhNzcXD7++GMyMzOJioqisrLymNq9w+HwL+ZRoa2hoYGGhoag6CdQsYQDTd5Bpn1KWvtxZBUVFdTV1TF06FBSUlL8J+GICI2NjSxdupSjR48yZswYHA4Hq1evPqFPj8dDUVERa9eu5dVXX2Xt2rX861//wuWqp6WlhfXr1/un/1mtVh5//HH/3O5hw4bx5z//mR07dpCUlOQ/kae+vp69e/eyZMkS/7FpDQ0NVFdX+w+IaF+MpELbrFmzAjKrKBD9BCqWcKD/soJYZmYm559/PjabjVGjRnHBBRcQFxeH0+nEYrGQmJiIx+PhmWeeISIiggsuuMC/aKZdS0sLTzzxBM888wxLly6loqICp9PJhAkT6JvUl759k7Db7axevZrm5mYABg8ezAMPPMBnn33Gn/70J6ZPn869997L+++/T1tbG7W1tdTW1pKbm8snn3xCSUkJ/fr1w+Fw+Lew1dGRUmeWJu8gFhERQU1NjX9xjDEGu91OcnIyffv2xW63c8011/DYY4+xa9cuZs+ezY4dO9iyZYu/j/ZDjRMSEvjxj3/MHXfcwdChQ8nIyMBqsWKxWJg0aRJNTU2sXbvWf2P0uuuuY86cOfz3f/8369atY8GCBUydOtU/+gZobm7G5XJx+PBhHn74YYqLi0lOTiYxMVHLJUqdYZq8g1T7IcPjxo1j6tSpTJkyhe9973tMnjyZ0aNH+8+U3LZtG3Fxcdx+++3MmzePqKgonn32WX8/FouFiRMn+mvWnenbty/jxo2jrKzMP3IXER599FEyMjK44447qKqq4sc//jGlpaX+Ta+ioqKIjIxk7dq1LFmyhF/84hc0NDRouUSps0D/hQWp9tkaFouFyZMnExcXR3V1Needdx4TJkzgvPPOIyoqCpvNht1uZ9OmTWzdupWrr76a1157jZaWFn9fY8aMITo62l8Pj576Q2pH30RbbH/aYvtTO/om+sy5j0sXvcyWLVt46aWXgC/2QCkvL/fvbTJ58mTS0tL8+3/X1tYSFRXlr4W/9dZb5OXlHfP+Sp2uzMxMRITly5ezfPlyRMT/yMrK6unwekTQTRVUXu1lh+joaMCbSN1uN4MGDWLUqFFs2LCBmJgYoqKiiIuLo6amhnfffZebb76ZV155hY0bNzJy5Eh/f6NGjWL9+vXs378fOWdip+9pi44nJSWFe+65hwMHDjBo0CDAe+jrP/7xDxYuXMhPf/pTrrnmGlJTU0lMTPRPAxwxYgR5eXkMGjQIm82GzWZj8ODB/r7DecfBcHfTTTf1eD/5+fkAPPfccyf01Vt/tqS3TavuP+gcM/+BZ3o6jK/Me5OwDocjGrvdzpEjRVRUlPvrz21tbdhsEUycOJFNmzZiDIwYMdz/erfbzfbtO0hMTISE9C7fp6XsAE1NzRjjIaV/f2xW7+/3srIympoaGTNmDHv27MUYQ/s6qoiICCwWKy0tzTgcMfTt25fk5H5YrV+MDcLl39ejN35tkzFmfFfXx48fbzZu3Hg2QzrrRIRA5Y27nl8DwCPfnRQU8fQ0ETnpz1dHWjYJETExscTFxWIMVFfXEB8fh8Vi9V8XsdDS0kJtbS2pqWm++bBfHDhks9lwOp1UVVWe9H1EhH79+mGMobyszP+PIikpCavVxq5du0lPT6ehocF/c7O1tdW3E6GbxsZGampqqKxs34/F0NraEjb/uHqj8vLygBxEEYh+AhVLOOh1ZZMBSTGn9Vv+bDp+iXlxcTEff/wxe/fuJTo6GlfRerZv3+5fPFNfX8+E6Fv53wcfJC3tevpffjk/+MUv/P2VZjQzd+5tXPqrd7p8z5rPnmXiFVcQWbOXN954g1GjRvmPnkq1VfP000+TNncuVZs2ERcXR21tLSKC1WqlsLCQc845h6SsLKZeeSVTp5xLTEwMRUVFpKWlhcWWqo/eeOLnRGQBsAC+2CQsnHzzm98ETn8P7UD0E6hYwoGOvINY+03L9jnTTqfTv3/28OHDGT58OLNnzyYnJ8d/SMNTTz1FXl4eV199Ne+//z4rV67095ecnNztBQ5Dhw5l4sSJbNu2zb8hVmZmJnfddRevvfYaN910E7m5uVRVfbHVektLC3369GHo0KGkpaUhIjQ1NfmnNravvuy4EjMcGGOeNsaMN8aMb993XakzLWSTt/SCuxTtS8ztdrv/9PasrCwuvfRSLrzwQsaMGUNtbS2JiYnExMTgcDgYMGAA8+fP58Ybb2T48OH8+te/9m8mBd4bPc11nZ9tcfznJ06c6D/Zvt3ChQuJj4/n888/Z/78+VRXVx+zIKe4uJgPPviA119/nc2bN7N161Ysli9KOu3Lm2trayktLeXTTz895pg3FRyysrL8szmOn+GRmZnZ0+EpQjB5i8i50Ds2sGrfsKqyspKioiJ/rc9ms5GWlsaFF17IuHHjGD16NOnp6VgsFqqqqsjPz+dPf/oTv/3tbzHG8Mtf/tJ/Mk5qairLH5nPx7+aTeXBbVQe3MbHv5rNx7+azfJH5h/z/tHR0YwdO5Zdu3b5l8AnJibygx/8gLfffpvbbruNESNG4PF4cLu9p4Pt3r2bXbt2sWXLFoqKirDb7VRWVrJ69Wr/uZ3tv5S2b9/OmjVrWL9+/Vn8rqruKCgo8N2UNkydOpWpU6f6n7fP/FA9q1vJW0SiROQnIvKWiLwpIv8lIlFf/srAEpErgJdEZOjZfu+e1HGzqnbGGAYMGMDMmTOZOnUqgwYNIjo6mqioKKKiovjggw/YsWMHixYtYufOnbz88su43W5/kj2Zjz/+2P9oaGhARHjjjTfYtm0bO3fu5PLLL8fhcPDAAw/w4IMPUlNTg9PpJCEhgcTERBITE6muriY3N5fNmzezfPlyPvroI7Zu3eofvcXExDBx4kQmTZrEhAkTzuS3T6mw1N0bls8DdcATvufXAi8A885EUJ0RkauAe4AfGGP2tu/73c3XhvQNpc7OT7RYLLS1tdHa2orD4fDv5Hf06FGMMVxwwQX84Q9/YP369axcuZJXX32V++67j8TERO68804++eQTRCA1NY27777b3+8jjzxywrxZYwyVlZXs3eudIhgfH8+1117L3/72N6644gruvvtuHnzwQf8vDoCKigreeecdxowZw/XXX092drZ/3ndbWxvl5eU4nc5uHZisetbtt98eNP0EKpZw0N3kPdwYc16H5/8WkS1dtg4wX337QaDOGLNKRFKA+SLiBF4F9hlj6rp6vTHmaeBp8M7DPRsxnw2FhYVs3LiRVatW0djY6D/VvaGhgdLSUtxuN7fccgsPP/ywf3/u++67D4BJkyaxqcU7f5zok79P+zza9hPlAW699VZWrVrFAw88wMsvv8yVV17JkiVLsFqt/p0PPR4PLS0tJCQkMG3aNBITE6mrq8Plcvlr+CkpKWfmm6MC5tvf/nbQ9BOoWMJBd2vem0XEvyxPRC4EVp2kfUD5RtgXAX1F5A3gRby/ePoAPwZ65dlbqampVFZW0traisvlYsKECfTr1w8RwW63ExMTw8qVK/nXv/7FNddcwxNPPEFJSQkA8fHxxMTEUl/voqam5qTv0z4SLykp8d/8tNvtPPzww3g8Hn7+859z9913M3ToUOLj4/2zSCIiIvB4PKxatYq//vWvvPzyy95DIIqKiI6Opm/fvmfwu6MCpbCwkMLCwqDoJ1CxhIOTJm8R2SYiW4ELgdUiki8iB4E1wMVnOjgRyRGR4SIy2hjj8sUxBlhtjHnYGHMbcBT47pmOJRjV1tZy3nnnYbFYiIiIwOVy0adPH+Lj4wHvkvpvfvOb3H///Vx11VW43W5uvfVWf3L1thP+8Y9/+G9odqV99N3xqLWBAwdy9913s23bNjZv3syLL77I0aNH/fuatLa2cuTIET777DNeeeUV3nzzTR599FG2bdtGVVUVTU1N1NfX6wKeIDd//nzmz5//5Q3PQj+BiiUcfNnIezYwB/g6MAiYCkzzfXzlmQxMRGYBrwA/BZ4SkTm+BH4OcH+HqYL7AJeIWLvoKmw5nU6GDRvGvffey7hx4+jXrx9paWlkZ2cD3oS7efNmUlJS+Otf/8qjjz7K0qVLWbt2LQBWm5WkpCSOHj3Khx9+2K0k2n7QQrsrrriCpKQk3nnnHcaOHcs555zD5MmT/ddramo4ePAgJSUllJWVUVJSws6dO/3nXXacx66U6r6TJm9jTMHJHmcqKBEZDzwG3ArcBvwZmCEiNqDNd0CxEZFbgTuAZ4wxbV33GJ6sVit9+/bl8OHDnHvuuSQnJzN58mRGjRrF+eefT1paGhaLBbfbzWeffcbgwYO57rrrWL16NQcPHgQg2hHNtGnTyMvLO2ZBz/FEhIiIiBOSd0REBFdeeSUrV66kuLiYmTNnsmLFimN+EbjdbqxWK9nZ2YwZM4ZrrrmGqKgoYmNj9ag0pU5RsC6PdwK/NcasAxCRA8D38Ja/Pb4kPhT4BnCTMWZHz4XaswoLC6mqqiIxMZGIiAj27dtHfn4+CQkJ5OTkUFdXh91uJyIigt/85je8++67/POf/+S9995jxtBZ2KxWRo8eTXl5uX9E3pXIyEh27959wmk96enptLW18eyzz3LZZZfx2GOP4XQ6vZtg4S3fHDp0iJEjR3LzzTeTnJxMS0sLH330ETNmzAjYrnB6Yr3qTYIyeRtjPhKR1A6fygUaO4yu+xhjdonId3yllF6nPeENHDgQ8I7C165di8ViobW1lYqKCiorK9m9ezfgrT8XFRXx+eef8/bbbzNr1izqXS7OPTeHRXdcxp133skNN9zAypUrueSSSzqdBbJ8+fJOb24mJyeTnZ3N4sWL+dGPfkRcXBxlZWX+BNp+mk95eTk7duwgOTmZl19+me3btwMwd+7cgHxP9MR61ZsEXfIWEZsxxm2MOdrh0zZggK+uPR+YJyLzemvi7shutzNkyBBcLheRkZHs3buXvn37kpiYeEyibZ8rfv/99/Pmm2/y2GOP8fSaCgoK8oERRERE8PTTTzNhwgRWrlzJ5Zdf7r/x2S46Oto/xe94F110Ec899xwbNmxg+vTprFu3DmMMIkJlZSV2u53S0lJ27drFrl27OHr0KJmZmeTk5NDa2upfNp+cnIzNdmo/lu3lFy3DBNZPf/rToOknULGEg6D421JELhKR+QDGGLeIHB9XI3AA+BnwfeBuY4ze5fJpaWlh+fLlrF27Fo/HQ1paGpdddhmzZs3yJzIRoaSkhLVr1/Lvf/+bOXPmkJqaytGjxbz99tsAJCQkMHnyZCwWC8uXL/cfSNwuKiqK+vr6TmemjBs3jri4OBYvXsyVV17J0aNHqa+v9193uVy43W727t3L3r172bdvH7GxsRw6dIiDBw+Sn5/PgQMHKC8vx+VyndKeJ+3bCWjJJLDmzJnDnDlzgqKfQMUSDnr0p1xELCISCzwFLBKR7wN0qGvje94MpAA3AN/rzTXuzhQWFuJ2u7nggguYO3cuv/71r5k4cSJDhgxhxowZ/vnUVquVwYMHc9ddd3HkyBEyMzOJj4/jJz/5Cbm5uQDExMRw8cUX09DQwIoVK3C5XOTl5fHJJ5+wY8cOrFZrp0vsIyMjueaaa3j77beZO3cuOTk5/uPRwDtHfMuWLfTp04eRI0cyd+5cJkyYgMfjITo6mqysLAYPHkxsbCyvvPIKy5cv1z1PgkReXh55eXlB0U+gYgkHPVo2Md6jWFwishhoA74mItHGmD8YY47PEO8Bbxljdp/1QINce9174MCB2O12APr06UNKSgptbW1ER0fzwQcf4HK5KCoqIiIigmuvvZYZ//Vnhg8fTl5SEjfffDNLliwBvFMQJ02axKpVq3jvvfcA74ZUo0aNYs6cOf6j2Y43Y8YMFi9ezPbt27nvvvuYN28eKSkpxMfH43a7iY6OZt++faSmphIfH09paSmDBw8mLS0Nm81GYmIiW7duxRhDSkoKSUlJtLS0+L8m1TNuu+024PT30A5EP4GKJRwES83bDWQAi4FbReQxoNkYs0hEJgNlxpjf9WiEQay97n08h8PBsGHDqK2t9Z96Y7FYyMrKIi8vj927dzNy5Eiee+45rrrqKm655RYyMjKIiIggIyODvn37UlxcTHJysn/l5Mn2q54yZQoAK1as4L/+678YNWoULS0t/vduP7B4+fLlDB8+nJiYGNLT07Farf6ZIu1z1K1WK4cPHyY2NrbTr02p3i5YioPvAsXGmH8BG4HbgUTftfF4N8VSp6C2tpaIiAiSk5P9MzASExPJyMigrq6OvXv3MGLECB5//HE2bdrE559/7i91xMbGkp2dfcKNy66kpKQwYsQIVq5cicVi4e677yYvL89fuy4vL/dvFVtbW0t6ejppaWl4PB5KS0v9JwKNHj2aoUOHkp2d7f+rQil1rGAZeTcCw0XkP/DekHwImCAi1wGP94a9uwOtfSqh1WplwoQJZGVlsX79ev8J8m1tbfS1WKisrOKhhx7iiSeeoLCwkAcffJDrr7+em2+++YQ+P/nkky4TudVqpaqqigsvvJA33niDsrIypk2bRnZ22MghxgAAIABJREFUNna7nbq6OkQEi8VCSkoKaWlpLF26lFGjRhETE4Pb7cZms/lvsHb114RSyisoRt7GmCKgELgX+Ikx5kHg98BKTdynJyUlBafTyWWXXca0adOYNGmSf/Mqi8WCxSI8+eSTLF68mHvvvZfLL7+chx9++ISFON01ceJE6uvr2b59O1arlf/8z/9k586d/tJJU1PT/2vv3OPbqq58/12y7MiyLb8fcRI7zsMJdggJMU6AUJI0UNrybOkl9A7l2RdtL72dmT6403l0GDrtwLQzbXnNMCUFOiSltAXa0kJbyDQBQhLABCeYxPgRvxLb8UOWbVnWun+cI+EktuPEsiQ7+/v56GMd6Wifn7bPWdpn7bXXorm5mV27drFz5062bNmC2+0mIyODvLw8EyliMEyQeBl5A/wH8CtV3W1vv2RPaBomwRtvvEFDQwMZGRmce+654Xjs0FJ4h8PBunXruP322zn77LP57ne/y6ZNm7jjjjt46qmnKCkpOaXjnX++Vdx5x44drFixgquvvpp77rmH/Px8qqurEREaGhrw+/0kJiZSVlY2ZltmxWR88Dd/8zdx006ktMwE4uaKUNVGVd0dSjhlDPfkCQaDlJWVsXr1aoqKinA4HKSlpbFx40aKi4tJSHAAwp49e8jJyeGaa66hu7ubBx54AKfTyW233RYuvTZR+vv7cblc1NTUANbioC984Qu89tpr4WyGfr+fhoYGBgcHSU5Oxufz0dXVxeHDh8P7jPSDm8RVsWXjxo1s3LgxLtqJlJaZQNwY7xDGTRI5fD4fwWCQNWvWcNZZZzFnzhxSUlJ45513SExMJDk5GafTicPhoLu7m97eXm6++WaSk5N58MEH6ejo4KqrruLdd9+d0PHq6+u59tprcblc3HbbbeHXKyoqAE6oFu/1esMJspxOJ4FAIGyofT7fCX7w0fD7/Rw8eDCchtYQed54443wOoBYtxMpLTOBuDPehsgRKvTrdrtJSkpi2bJltLS00N/fT19fH7NmzQpXvUlISGB4eJjm5mZuueUWFi9ezObNmwkEAlx77bU899xz4x6rra2Nj3/84/h8PrZs2cKyZcvC74XypBz/u+z3+/n1r3/Nww8/HPZ7JyUlhZfgd3Z2kpGRMa7LpLGxkQMHDkwoQX8gEKC1tXVCdTwN7/PlL3+ZL3/5y3HRTqS0zASM8Z7BjLZcfNmyZXg8HnJycnAkJJDsfn/BjdPpZOnSpbzzzjt89rOfZcmSJfzyl79k0aJFfPGLX+QXv/jFCaNnsHJ8f+tb36K/v58nnniCs88++5j3MzMzSUpKOsF4OxwOOjs7eemll3jppZdoa2ujtbWV5uZm9u7dS0tLC01NTeN+x3nz5k04pLC9vZ3m5mYOHz6M1+sd9bsYDNOFeJqwNEwxXV1dVFdXU1FRQUZGBq/0pdDT0016eno4G9+CBQsYHh7mtdde44477uAHP/gBjz/+OHfddRdPPPEEjY2NfPrTnyYtLQ2wDPe9997L4OAgmzdvZs6cOXR2dp5w7NmzZ9Pe3k52dvYxrw0MDIRrXb711lssWbKEwsJCPB4PLS0tJzXKpxJSmJOTA1g5Wo7PPjiZydHpXuDaMD0xI+8ziI6ODl5//XX279/PqlWrSEx0hkP4Qi6W6upq2tracLvd/OEPf+Cb3/wmRUVFPPbYY9xzzz3U1tZy7733Mm/ePMrLy/nRj36EqrJ582bOP/98MjMzR33Mnj0bv99PX19f+NHU1BQuULxr165wgq1t27bh9/vJz88/7QyDoxHKrOjxeE4oAhFKJ3s6k6Oq+pCqVqhqxXgrUA2GSGJG3mcAgUCA9vZ25syZw6WXXkp/fz9lZWVk7t/F0FCA9PR0gsEgCQkJzJ49O7xcXVV56qmnePfdd/nZz37Gpk2bKC8v59Zbb+Wqq64iPT2doaEhtmzZQmFh4bgaZs+ezeLFi4+pxNPf309aWhqDg4O88sorHD16lEOHDjE0NEQgEGDdunVA5HNzh9xJIzHpZKcvxcXFYxb0KC4upq6uLrqCooQx3mcAIV8vQHl5OarK8PAww8MBVINkZ2fT29tLcnIyfX199PZa2QjcbjdOp5PGxkZWr17No48+yqpVq/jd737H5z//efbt28fWrVspLy8f1VUykuzsbFpbW09wSbhcLrq7uwkGg5x33nnk5eXR19dHYWEhfX19YVdH6AcoJycnoqPxEKMZdIPF3XdHJq1QJNoZrY3xjHOkqjTFI8Z4nwGEDGDoL1gGPRAIkJ6ewdqPfQyPx8P+/fupra3F6XTi8/l47733SEpKspbSZ2XxyU9+kieeeIKVK1eyZcsWhoaGJpTxb3h4mOeff54VK1ZQVVV1zHtHjhzh8OHDZGdn09fXx4033khTUxO9vb00NzeTkpKCx+Ohp6cnHIVSUFAQwd4xnIwLLrggbtqJlJaZgPF5nwGEfL0jR6w5OTnMmjWLtLRUPvWpT3Heeecxd+5cPB4PZWVlXHrppcyfPx94P29JVlYWmzZt4vXXX0dEJpyq9aWXXqK+vp7XXnvthPf6+/sZHByku7ub5uZmOjs7w5OG6enpuN1uenp6cLvdFBYWHvMDZIgOO3bsYMeOHXHRTqS0zATMyPsMJSEhgcREy/impqaSmZkZzkvi9/s5cODAMSXPEhISwpEfmzZtCo/AJ8Jjjz1GXl5eeHJ0NESEw4cP09bWFo4gcbvdBAIBduzYwcaNG82IO8LMnz+f+vr6Ud8rLi4OP7/zzjuByefQjkQ7kdIyEzDG+wylr6+Prq4uPJ40RIS8vDzq6uqYM2cOR48eDa/AHFkKLS8vL/x806ZN/PSnPw0bcL/fP+oqx9bWVl544YVw4qnRcDqd1NXV0dXVxRNPPEFmZiarV6+2CkW88w5vvvkm/f39XH/99VPi7z5Tqa+vPyH23jB9MFfCGYiIsHPnznCss4gwNDREWloabreb/v5+Ojs78fv9BINBent7GRoaYvv27fZE5zC5ublcc801fOITn+BLX/oSZ599NgkJCScc6+GHHw6PuMea1HQ6nfT391NVVYXL5SIhIYGEhASKioq47LLLUFUKCwuprq6mrKzMGHCDAePzPmOprKy0452txTb19fW8++67LFy4kJUrV/K1r32NdevWcdZZZ7FixYpjQugSEhLo6urC6XTy7LPPcsEFF7B+/Xq2bt16THHiQCDAj3/8YzZu3HhSg5ubm8ucOXPIycnB4/EgIuEkV9dffz25ubkMDAyccqIsg2GmMq2Mt4hcLiJ3i8gPRCRHRBJjrWm6Yvm5M8Kj5ezsbPLy8vD7/Tz99NMcPHiQyy+/nMWLF7Ny5UoyMjKO+XwotC4xMZGUlBSOHj3KTTfdxNKlS9myZQuqyoMPPkhTUxNvvz1+vWin08mKFStobm7myJEjpKWlsWDBAlatWhWul+nxeMJ1Lg0GwzRym4jIKuAB4HbgauCHwBMi8idV7T7JZ83y5ZOQkZFBWVkZ3/ve93j66afZtWsXZ511FikpKSxYsIDu7m62b99OX1/fMZ9zOBwkJyeTlZVFf38/c+fO5eabb+bmm28G4AMf+AAtLS1hF81YHDhwgJaWFkSE7OxslixZEl4639HREU4x29XVRVZWFo2NjccUXDZMHd///vfjpp1IaZkJTBvjDZQCv1fVp4GnReSzwEeBYRH5LTA8VjpZVX0IeAigoqLCzNCMgsPhYGBggAULFnD22WdTVFREXl4eubm5BINBGhsbw4t4RkNESE1NDYcUrl+/nm3btoUN8ngEAgG6u7vDk6PvvfceTz75JB/5yEcoKChg2bJlFBcXM2vWLFwuF/X19dTW1gKYUmlRYMWKFXHTTqS0zASmk/F+FbhRRC5Q1R2q+qCIOIAbgG0nG30bTk5OTk44P8maNWtoaGiwY8HT2LdvHzU1Nfh8Pl599dUx2xARcnJyeOutt8jMzJzwsb1eLyUlJQQCAVJSUqiuriYxMZGVK1eSlJREQUFBOAthS0sLBQUFJCYmhnN+G6aOF154AWDSRRAi0U6ktMwE4vqsF5EVwCAgqlotIruBi0SkU1X3q+r9InIh8NeAqY80SZxOJwMDA/T19dHa2kp6ejq7d+8mKyuLiy66iJSUFJxOJ01NTbS0tBzz2dTU1GMmK49nvBjv5ORkOjs7KSgoYNasWbS3tzMwMMDBgwcJBoPhAhKlpaX09fXR0NBAeno6L7/8Mhs3bmTRokUR6wPDidx1113A5A1mJNqJlJaZQNwabxH5MJar41fAOhH5O+Bh4GvAVSKSp6rbgJ2ASUoRIUpLS8N/fT4fc+fOJTMzk9/+9re8+uqr5OTkMGfOnLD7xOfz0drayp49e047P3ZoBN3a2kpSUhIZGRkUFhbidrs5evQomzdvJiMjgxUrVnDFFVeQkZFBfX09dXV11NbWGuNtOCOJO+Nt17BMAb4EfEFVnxaR84HHgDuBu4GbgX8QkTbgAizft2ESjMxnXVpaSk1NDYsWLWLFihW4XC46OjqoqqoiNzcXr9cLWAt9PB4PbW1tEdEwPDyMx+OhuLiY7Oxs/H4/Xq+XWbNmUVZWxoc//GEWL15MYmIiCxYsICUlhcrKyogcezT8fr+ZGDXELXFnvO1JR6+I7AI8IpKoqi+LyPXAVuArqvr3IjIXWAl8VVUbYql5JhDKZw1QW1vL3r17Aavyjs/no6ysjBtuuIFDhw6xfv16qqur6ejoID09nYyMjEmXFps1axZNTU0kJSWRn5+Px+OhpKSElStX0traytq1aykvLw+HNqakpLB+/XqGh4dpbW2dkmyDoRJrYCZGDfFH3BnvEbQCHwSeBoZUdaeI3AB8T0TeVNWDwKGYKpxBjMxnfbzrpKenh9TUVPLz8/F6vWzfvp0lS5YQCARYvHgxCxYsYNeuXfT395/28UORJocPH2bOnDl4vV7S0tKoqKigvLwcp9MZTikbKhwRDAaprq4OfzbSuU9CuVwmUmLNYIg2cWe8RSRBVYdV9T4R2QI8ICKfB3yq+j8iUgWMPftlOC1G5rN2uVwsX74ceL/iu9vtZunSpRw6dIiCgoJwlr/MzEy6urpYuHAh27ZtC4fwJScnn5IxT0tLw+fzUVxcTGJiImlpafT09NDU1MSGDRvCRRqamppITU1l0aJFeL1eBgYGcLlcU5Jt8FRKrM1kHnzwwbhpJ1JaZgJxYbxFZC1QoqqPquqwiCSpql9VrxORJ4DvA6+IiBO4GDDlv6PESKOemJjIhg0bqKmpoa6ujra2NgoLC7nqqqvo7u4mKSmJp556imAwSGpqKg0NE/dm9fX1oaocOXKE/Px88vPzSU5ODie2ysnJIT09nfz8fLKzs3G73bhcLoAxXSaTqUs5mc/ONJYsWRI37URKy0wgpmeliDhEJBV4EPiGiHwOQFX9IuKyn28C/gfIBdYBV6qqcZfEgGAwSCAQYNGiRQwNDZGens6yZcvCI+ENGzawdOlSCgoK8Pv9uN1uUlJSyM/PJy0tjaSkpFEfycnJiEh4tWZ3dze9vb2oarjWJVh3BL29vTgcDkSEhISEE/KUj2QydSkn89mZxjPPPMMzzzwTF+1ESstMIKYjb1UNYk1ObsZyhVwgIsmq+j1VHRix338BiMgsVR0coznDFNPf309PTw8NDQ14vV4qKiooLS3l6NGjpKamMjAwgMfjwev1kpeXRyAQwOfzkZCQgM/nG3N1Zl9fXziee3BwkMbGRvx+P6tWrQpHeXg8Ho4ePUpTUxOdnZ2sXbsWh8MxbjRIyI/vdDo5ePDgKUWNmJqW73PvvfcCcMUVV8S8nUhpmQnEhdsEyw1SBGwGbhORfwUGVfUbInIBMKCqe4ATE0YbokbIkC1fvjw8sZmUlITL5cLtdlNeXk5xcTFFRUXh8mahaJDQCFZEcLlcYX/4rFmzSE9PB0BVmT17Nnl5eXi9Xvbt28fixYu5+OKLqaiowO/3h1dY1tTUkJKSMm40SMjlc/DgwVOOGjE1LQ3xTrwY718Bn1DVP9irKu8CHrHfOw/YAuEwQkOMCBm0QCBAXl4eTqfzmNdUlcsvv5y2tjZefPFFEhISyMjIwOfz4XA4cDgc4eLHIVwuF5mZmQwNDeFwOCgpKaGoqIh9+/bh9XqZO3cu5557Lk6nE7fbzdq1a6mursbj8YSLO5wsGsREjRhmIvFivPuBJSLyaeBzwLeBShH5JPDvxmjHFyOr0efl5eHz+ejs7OS9994jJSWFrq4uMjMzycnJITExMbzIp6+vj+bmZvx+f9jHnZZm5ROfN28eiYmJHD16lIyMDNxuN2VlZVxyySU0NTWRnJwczn2SmprKc889N+Gl8SeLGvH5fFRVVbF8+XJcLtcpT1SarJWGWBAXxltVm0WkEfgm1qrKZ0RkPXDAGO74Y2Q1ep/PR1dXF+3t7YgIjY2N1NfX43K5yM3NxefzkZSUxMDAAAkJCSxcuJCOjg4cDke41Fpvby9+v5/8/HxUlYKCAkpLSyksLGTHjh1cfPHFzJs3j87OTrZv386f//xnAoEACxYsiMjS+KqqKnbv3g1Yi5JCi5Um6jYxWSvjl+Li4nGzWhYXF1NXVxc9QREkLoy3zX8Av1LV3fb2S/aEpiHOCFWjB8sP7vV6ycrKIj09ndWrV5OXl8eiRYvIysrizTffJBAIcODAAQoKCigpKcHpdPLiiy+SnJzMvHnz6OrqIikpiczMTNra2vB6vQwPD7Nz506am5sZHh5mzpw51NbW8vzzz9PR0cHixYtZs2YNgUCA9vb2Sa2wDMW0h0beoe9lsHj00Ufjpp1TbeNkhvlk6Yrjmbgx3qraCDSKiKiFMdzTAIfDEXadhAxeWVkZvb29XH311RQXF1NfX09SUhK1tbWkpaXR3t5Obm4uYP0QHD16lMOHD9PZ2UlKSgr79u0jJyeHWbNmsXDhQpYuXcqf//xn5s6dy+WXX86BAwdYv349DofjGBfO6a6wdLvdrFmzJrw9kyYqJ1ohfjwiNVcQiXbMvMX7xI3xDmHcJNOPkZEZXV1d1NbW4nA46O/vJzs7m6amJvx+P0eOHKGxsZGPfexj7N69G6fTSVdXFw6HA4/Hw7x58/B4PBw+fJi+vj78fj+FhYWUlJTQ0tJCbm4uCxYsICcnJxzh4nA4qKurw+PxxLgX4pNIVIjfsmULANddd13M24mUlplA3Blvw/TF6/Xy+9//nvT0dAYGBnj99de57LLLWLZsGR/4wAf4yU9+wqWXXorH48Hn84WXvG/cuJF33nmHXbt24Xa7WbRoES6Xi56eHrKzswkGg6xcuTLsa3c4HOTk5NDe3s4Pf/hDFi9eTG5urlnKPkXcf//9wOQNZiTaiZSWmcCZve7XEFF27tzJm2++ya5du+jq6qK5uZna2loKCgro7OwMV4Hv7OzE6XTS19dHRkYGg4OD5OfnMzAwQGtrK4cOHcLhcHDllVeyevVqli5dSlZWFvX19QSDQfLy8hgYGOCxxx5j+/bt7Nu3j8zMzNPOJ24wTEfMyHsGEC95OCorK3n77bfx+/3MmTOHq666igsvvJCGhgY6OztxuVykpKTQ1NTE7Nmzw3nD09PTaWtro7KykmAwSG5uLosXL2b58uX4/X56enp47rnnUFVWrlxJZ2cnxcXFnH/++bS0tPDRj36Ujo4OnE7nCe6TeOkbgyHSGOM9AxiZi3syk22TjdxITU3l1ltvDcdMj1xinpKSwuzZs3n22WfZv38/6enp+Hw+UlNTUVXmzp3L6tWrSU9PD7tM9u7dS1lZGVVVVbz99tukpaWRkZHBoUNWapu1a9dSUlJCUlISHR0do2qKVN8YDPGGMd4zgEjl4Yhk5Ibf7z8mn8iCBQtoaWkhIyODc845h9zcXHbu3El/fz8VFRW43W5aW1tJS0uju7ubZ555hqGhIVwuFxUVFbS3t7NixQoKCwuprq5m+fLl4cLEwWCQpKSkUb+/yVFimKkY4z0DiFQejpGLbybLaFVo3G43eXl5VFZW4nK5GB4epquri5qaGjIzM8NFiAsLCyktLUVVWbZsGZ2dnSQnJ4eLMIwM64Pxv7/JUTJ5nnzyybhpJ1JaZgLGeBvCjFx8M1lC8bihqjghw7to0SLcbjfBYJCioiLq6uo4dOgQg4ODnHPOOaxYsYKkpCQqKytJTEwErGyGPp+P/v5+48OOAZEqdBGJdqai6MZ0xZz9hjHxer388Y9/DBccPhVC+UQCgQA9PT14vd5wgqq9e/dy6NAh9u/fT0JCAg6HI2ycq6qqCAQClJSUkJWVhdvtpqSkhMrKSoqLizl06BBvvPEGXV1dU/CNDaPxyCOP8Mgjj8RFO5HSMhMwI2/DqAQCAZ577jnefvttADZs2HBa7bhcLrxeL8FgEK/XS0NDA7W1tSxdupRFixYRDAb54Ac/SFdXF8FgkF27dpGUlMS5555LamoqwWAQv99PSUkJPp+PhoYGGhsbyc/PJysrK5Jf2TAGIWN50003xbydSGmZCRjjPUOZrHuhvb097EKprKw8bR1er5fW1laKiorweDzH5AJvampCVcnMzGTVqlW0t7fT19dHVlYWAwMD1NTUUFhYiN9vpXEPZRrMz8+f8NJug2GmYoz3DGWyIXIh3+KaNWtOO+FTSEeoUHBopBxKBDUyz3YoWsTtdpOTk0N1dTV79+4lGAyG/eQOh4OsrCwz4jZEjPGyDsZ7xsFpa7xDCaxirSNemWyIXKQmL/Py8sJGNzRxGboTOD7PttPpJCMjg+rq6nBe7NLS0nCmv9MhElkH45nxEk/BxJNPnamMZ5zjPePgtDubRaRcVd82hnt84i1EzufzhUuhjaerpqaGPXv20NjYyCWXXDLhmpNjEYnY9XgmEomnDNOTaWW8ReRDwHdE5BOq+m6s9RhOTsh4FhQU4PF4TnonUFpaSmNjI8PDwzQ2Nk462VTI/TPayN8wMX7zm9/ETTuR0jITmDZnsYhcCfwdVqWdd+UU7mlE5DMisktEdh05cmTqRBpOICcnh8LCQtxuNzt37gyPvsfC5XJxySWXUF5ePmbu5lDkykQSUYXcP6EcKSc7vuFEQnMR8dBOpLTMBOLeeIuIQ0QcWEWJh1V1u4jkA38pIv8sIitFJG28NlT1IVWtUNWKUBEAQ3QIGc89e/bw8ssvs3PnzpN+JuQLH8tlEpqM7enpobW1lUAgMG57wWCQYDBIamrqMaGLholx3333cd9998VFO5HSMhOIe+MN5NlVdS4EMkXkSeAxLJdPJnAHsDyG+gwToLKykvPPP3/UsMOBgQGqqqoYGBiYUFsulwuHw0FXVxd79uwJf3YsQx6KeHE4HAwMDJgR+CmydetWtm7dGhftRErLTCCujbeIXAbUiciHVbUXWA2cA+xQ1X9W1c8CLcCnYqnTcHJSU1PZsGHDqJOVNTU17N27l5qamhPeCwQCJxjlgYEBgsEgg4OD9PT00NbWRk1NDc3NzbS3t5/QRmhpfuiWeyK+93hi/vz5iMioDxNNMnWEwghHe8yfPz/W8uJ3wtI23HcBvwDOFpE/qWqfiJwFDI8IFTwAlIlIgqoOx1KzYXzGWjhUWlp6zN+RtLe309jYyOHDhykrK8PpdIYNb05ODikpKeHnXV1d4QnK44818kcjnqJwJoKJKIkN8R5GGJcjbxG5GPhn4IvAPcDHgVQAVQ3YBYpVRG4Dbgf+yxjuyOHz+XjllVci7lro7OzklVdeobOz85jXXS7XMZXbRxIqRDwwMHDCqNrpdDJ37lzmzp2Ly+WioKAgHMsd8osb94hhKoiHUXlcGm9gPvAZVX1FVXcD1cC/iogTQEQS7BH4lcBNqvp27KTOPKqqqti9ezdVVVURbffQoUM0NDSEiykcz2hRJE6nk7KyMoqKisKj6vEMc8jNkpSUhMfjifoE5VtvvTXmRX26D+MaiT/q6upQ1VEfQFQMu8TT7djxqyZFxKmqARE5F/gC8A+q2mC/lwQkqeoppbwTkSPA2EvSpoYc4ERnbHwwmjYBkoF+IJIniACzgMEx2nUACcAwkDWKruP3Hc0iO4FEYAgIHNdmcJzPTZQlqnpMdJOIfAb4TOh94J1RPhcP54DREPvjn0xDsapOKCQurox3iFGMeDKW7/sNVf167JSdHiKyS1UrYq1jNOJV20zTFQ/fx2iI/fEjqSEu3CYiUi4iF4tIHsBxhtuhqv3Al4CNInL6Ke4MBoNhhhDzaBMR+TDwHaAWSBSRz6hqU+h9VQ3ai3Q6gd8C78VGqcFgMMQPMR15i8g64N+A21T1asAPnGW/lxDaT1WDqtoB3K2q03F9+0OxFjAO8aptpumKh+9jNMT++BAhDTH1edsRIwWq+icRKQD2ADuBNuBlVX3EnqxEVfeYNLAGg8FgEdORt6ruU9U/2Zu3AvfZI/BXgI+IyHzgYqDZ3t8YboPBYCBOo00AROS3wB2qeuKaaYPBYDjDiZdoEzlu++NAHtAbG0UGw9RyKimNZyrx0AfTWUNcGO+QO0REZonIrcC3gBtVtSW2yqaekROz8UI8apoOTOQiFJFyOLNdgPHQByKSGAcaJtUPcWG8RxDEyhL4MVXdG2sxU4mILAdQ1eF4MpYi8gHgRhGZFWstIxGR+SKSKiKeWGs5HhHJhGMGIaMacbEqQT0uIoujKO94DZeLyN0i8gMRyQkZsSgePx764Erg30Vks73GJOoVrSPRD3FlvFV1SFV/o6qjLS+eMYjIQmCPiGyF+DHgIvIR4HGs9AH+Ea/H9NZSrAyTTwHfB74iIsl27H/MEZHLsS7Ch0Xk0yJSZCdNcxy332lXgoqg1lXAA1gBASnAD4GPikh6lI4fD32wDLgf+BlWVNsXgE+JyJwoaohIP8TthOVMxg6L/CFWbvJaVf3QiPccdvGJaOoRwA08CDyuqr+1R7iKlT+mQ2KUcldEzgGeBG7Dyrdytap+LtRPseivEdpR24iiAAALH0lEQVRKgT8Cn8RKplYCLAO+rqoH7X4NXZhvAN2qepFYlaBuwMpxsQU4YOern2q91wOXqOot9vZngQrgaawFcMNT5Uaw++J1oDfGffBBrHUl19vblwPrgAbgJ6raNYXHDv2gR+RciIvRy5mEbWxagdeAlYCKyM9F5BwRWRgLQ2Rn2O0DDgO9IjIXeA74HrBNRFbZdwdRGymNOFYB8JKqvoSVXXKdiPwA+JdY9dcIhoHnVXWbqv4E+C8sA3W3iMy1DWE8VYJ6FSgUkQsAVPVBrLUVNwApU+n/tdteC2TFog9EJJRveBcwT0SusXU9C2zD+tHNmUoNRPhcMMY7SoSM0Qhjkwtco6qXAYuxLvpie9+ou1DESrfbC3wEuAZ4VFVvw7pD+LmIFER5cmfkiPUGEXkIqAJ+jDVSPAr8PxFJiaImS9j7I6jDwCoR+T8AqtqIdTG+g/UjE/NKUCKyQkTOEpEyVa0FdgMXichSW/P9WC6yv56i4y8TkSUistzOALoaWEF0+2ADcKuIJKtqN7AZuECsugGo6tPAAPDlKdQQ8XMh5rlNzgTsScALRKQaqFLVOuB3wEIRKcTyP74JfAP4Y7TcEyKyBusHo0FVXxaRfwdexPpRvw6si9v2lUbtB0VE1gNXi8hbqvqfIrIIKAMcqvode58OrGRlg9HSNYa227BG2kdV9VFVbRCRJqz+m00MK0GJlTfoIeBXWD8mfwc8DHwNuEpE8lR1G9aq5oiXF7LnUL4DvIzVB3er6jN2HwSi1AeXAd/F8i/32y//EcgHrhSRfFXdinX9nTUVbjiZqqpgOkZCcfOIzAPYgJW796+AXwP/ClwNZGCNKr3Ah+x9fw0URUnXJVgTk18FnrdPrjlANrAP64QvxLqlfgvIj5KuS+2T+Cv2ca+zX3fZOm+2t68HXgIyo/i/PF7btfbrVwAvAH9pb98NdGGtDl6F5a7IGaW927BGwuUR1ilYxvg3wJX2a+cDB7F+VIqBvwf+BDyB5e89O8IaKoD9WCNMwRpR/hvWgNEx1X1gt70c6w4t9H/KwbrjzcHKK38Llhvl50AjsHwKNFxsX+drIn0uROWkP5MfWMv+v2g/L7aN4X/aF9F1wEVR1iO2Ifwx1uQfWLdv+2xd8+0T/CHg27aBjPiFNYauZOAR4OP2azcD/4T1Q5MLrAdqgK323ynXNQFtFwOLgHL74vsJcAT45IjP/9h+3WlvJ2AlYHs60kbzON3fAv4CSLS3K4E6rFBcgLlYPzwRHzAAlwE3jNhea/9YJNjbTrsPnp2qPrC/74NY5RQrsOZxHgXeBTbZ+6Rh/cAUTpGGG4HKqTgXpvzEP9MfwE1Yk0KZ9naubcC/DcyzXxPsyJ8o6rrLfrjs7f+wT/Sv2duz7L/pUdb1TfsiWwt025p+Z/dXOdYP4CWhvosTbf+CdZeSgjXBmmfvH7pAz8VyVxSNaCsJSJ1ivbfbxsIz4rWLsEabC6PQX7NHPE8FfjNiOzf0+hRruBArxLTL7g+n3QeHgHOieO5E/FwwE5ZTjKo+AvwBuFNE0tVKabsd69e+zN5H1f4vTjUjoji2Y11Qj4jIfVhlxx7CWqCzRFVDvuSeKOv6b6zbzM8DD6k1kXM7ltFeqar1qvq8WpODUeEk2j6P5W76oKr2qRVJdASsYtn25/bZ+9wealNV/XqKJfxOQW+CfYz7sEJAHxCRdBFJVNX/wZr4nbJ5FXvyGz12hbQTmCtW/dmbsM479xT2QUjDdqwQvFvt/hi2++A5oG8qjj1Cw8i01hE/F4zxjg4/s/9+U0Sy1Jr1fw3rdnvKGSPE7/dYt28vYIXgXa9Wsec/YtV/BKZ2+fBoulT1gFql7n4OpNgX+EGsIhwLx/pcDLXVYkXprJcYVoISkbUicoN9/GGxaryiqtdhXeffB24RkS9guXoCYzY2+eMH5MRFVP1YBVe+BnwO+KqqnlhBOrIaQgb8ZeAZ+7mKyHVYYbr9YzYWGQ3DI/sh0ueCiTaJDq9hLXj5X8BLIvJLLJ/phVE6fgL2xWqfvE57JLDHfgBgR05cjOWiiIWukbPsr2O5J/5RROqATwCXh/aNI23JWK6xbcCPJcqVoOy2QwusRERSVPUBVfWLiEtVB1R1k4jcguXaOQdrEvPQFB8/OOI8Q1UHxVqQ8hdY8wb7InH8k2gI2HcbQ3Z/OLEWVP0V1mClabx2I6Qh3A+RPhfMCssIISIVwFF7lDjeftdhGfI3NQppAOxwsVuxJtMOqeqj9usJWHYpaD+vBH4KXKFRyCszji7HCE2X27rmAN9V1eqp1nWK2q4FrgQ+o6pbROQXwI9U9YXRwr3EijOO+GjPbvurWK6Qc4DXVfV7Y+w3a4RLLKrHF5E7gadUdX+kj38KGq4B3lLVA7HSYO836XPBGO8IICKXYE1c/QJrwm9KToxTxb4texwrLCwIfB34tareab9/TEyriORqFMrMTUBXoqoOjdg/SVX9o7UVS212nO5cVX1eYlwJSkS+AhRhuQduw1rwMaiq3xBrReXAVGo4yfEvBI7oFOfmn0AfHI3kiP80NJxvP4/I/8H4vCeJfdt8Ltat83vA34q1qOT4/T4oItFyR4RIAl5U1cdV9b+BDwGbROSfIHxLf7F9woEVjx4PuoZEZJ2I/F97/6GxGoqlNtsQvGB/LtaVoH4FtKrqH7CiST6PtZYA4LwoaBjv+BVEJzf/yfrgaIw1VBLJ/4NGKVRmJj+wEhKF7mLux1p+W3rcPh6gOMq6VmHFj2aNeC00QgwtfinHGj2e8boiqQ3Lp1k6VTpHOV4hVljgp7HimP8WK4b6k6FzcyYf/0zUEJUT60x7YE1abMaK+/3fWDlMYqXlPmDnca/dCnwpxn0Ul7pOR9vxFyXwcSx/+eyp0DeO7m9hrZa8wt5eTxTj4WN9/DNNg/F5T5KRE1MikqmqR+3n3wY+CqQDH9YoTbaN0BX2E4vIs1i3bteqaquIfB0rxvwm7DDzM11XJLSJVcDiL7CWz1+nUS4oIiLzsBYI7ba3o5ouN9bHP9M0GOM9CUb+U0TkH7EWcPxSrfjOW4B7gbWq+nYMdX0Dywf7MayoDT/WTPi1RldktYlVleYS4KDGsKDIVE6MTofjnykajPE+TY672L+LNRmxUa3Y0gysiYpfq2pVjHVdqKoX2ttnYcWidqiV2fCM1xXv2gyGsTDG+zQ47mK/B2sC6wrbcIdigY8Jd4u1rmjqmC66bD1xq81gGA9jvCeBiNyLlRHsSttwx6RUmNE1eeJZm8EwGibO+zQRkSJgCXF2sRtdp048azMYxsKMvCdBaEIi3i52o+vUiWdtBsNoGONtMBgM0xDjNjEYDIZpiDHeBoPBMA0xxttgMBimIcZ4GwwGwzTEGG+DwWCYhhjjbTAYDNMQY7wNBoNhGmKMt8FgMExDjPE2GAyGaYgx3gaDYUoRkfkiEtXCFGcCxngbDAbDNMQY7zMce1S0X0T+U0T2isjjIrJRRLaLyLsiUhlrjYYZgVNENotIlYg8KSLuWAua7hjjbQBYBPwbsBxYilXpei3wV8CdMdRlmDksAR5S1eVAD3B7jPVMe4zxNgC8p6pv2RVl3gb+YNfeewuYH1NlhplCo6put58/hjU4MEwCY7wNAIMjngdHbAcBZ/TlGGYgx+eeNrmoJ4kx3gaDIRoUicj59vPrgT/HUsxMwBhvg8EQDfYBN4pIFZAF3B9jPdMeU0nHYDAYpiFm5G0wGAzTEGO8DQaDYRpijLfBYDBMQ4zxNhgMhmmIMd4Gg8EwDTHG22AwGKYhxngbDAbDNOT/A22z4YdhxzIzAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = corner.corner(sampler.flatchain,\n", " quantiles=[0.16, 0.5, 0.84],\n", " truths=[mTrue, bTrue],\n", " labels=[\"m\", \"b\"], range=([-1.1,-0.8],[3.7, 4.8]),\n", " plot_contours=True, show_titles=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looks good!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Inference with ```approxposterior```**\n", "\n", "---\n", "\n", "Now let's see if we can derive similar constraints using ```approxposterior```.\n", "\n", "First, ```approxposterior``` requires a function that samples model parameters from the prior distributions." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def sampleFunction(n):\n", " \"\"\"\n", " docs\n", "\n", " Parameters\n", " ----------\n", " n : int\n", " Number of samples\n", "\n", " Returns\n", " -------\n", " sample : floats\n", " n x 3 array of floats samples from the prior\n", " \"\"\"\n", "\n", " # Sample model parameters given prior distributions\n", " m = np.random.uniform(low=-5, high=0.5, size=(n))\n", " b = np.random.uniform(low=0, high=10, size=(n))\n", " \n", " return np.array([m,b]).T\n", "# end function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the ```approxposterior``` parameters." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Define algorithm parameters\n", "m0 = 20 # Initial size of training set\n", "m = 10 # Number of new points to find each iteration\n", "nmax = 5 # Maximum number of iterations\n", "bounds = [(-5,0.5),(0.0,10.0)] # Prior bounds\n", "algorithm = \"bape\" # Use the Kandasamy et al. (2017) formalism\n", "\n", "# emcee MCMC parameters: Use the same MCMC parameters as the emcee-only analysis\n", "samplerKwargs = {\"nwalkers\" : 32} # emcee.EnsembleSampler parameters\n", "mcmcKwargs = {\"iterations\" : 5000} # emcee.EnsembleSampler.run_mcmc parameters\n", "\n", "# Data and uncertainties that we use to condition our model\n", "args = (x, obs, obserr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we create the initial training set by running the true forward model $m_0$ times. ```approxposterior``` learns on this training set and, each iteration, runs the forward model $m$ additional times in regions of parameter space that will most improve its owns predictive performance, iteratively improving the posterior distribution estimate." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Create a training set to condition the GP\n", "\n", "# Randomly sample initial conditions from the prior\n", "theta = np.array(sampleFunction(m0))\n", "\n", "# Evaluate forward model to compute log likelihood + lnprior for each theta\n", "y = list()\n", "for ii in range(len(theta)):\n", " y.append(logLikelihood(theta[ii], *args) + logPrior(theta[ii]))\n", "y = np.array(y)\n", "\n", "# We'll create the initial GP using approxposterior's built-in default\n", "# initialization. This default typically works well in many applications.\n", "gp = gpu.defaultGP(theta, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now initialize the ```ApproxPosterior``` object and we're ready to go!" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "ap = approx.ApproxPosterior(theta=theta, # Initial model parameters for inputs\n", " y=y, # Logprobability of each input\n", " gp=gp, # Initialize Gaussian Process\n", " lnprior=logPrior, # logprior function\n", " lnlike=logLikelihood, # loglikelihood function\n", " priorSample=sampleFunction, # Prior sample function\n", " algorithm=algorithm, # bape, agp, or alternate\n", " bounds=bounds) # Parameter bounds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run ```approxposterior```! Note that we set cache to False so ```approxposterior``` only saves the most recent sampler and MCMC chain instead of saving each full MCMC chain to a local HD5f file (see emcee v3 documentation for more details on how emcee caches data)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# Run!\n", "ap.run(m=m, nmax=nmax,estBurnin=True, mcmcKwargs=mcmcKwargs, cache=False,\n", " samplerKwargs=samplerKwargs, verbose=False, args=args, onlyLastMCMC=True,\n", " seed=seed)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, let's plot the joint and marginal posterior probability distributions to see how ```approxposterior``` did. In addition, as red points, we'll plot where ```approxposterior``` decided to evaluate the forward model to improve its own performance. As you'll see below, ```approxposterior``` preferentially runs the forward model in regions of high posterior probability density - it doesn't waste time on low likelihood regions of parameter space!" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAF+CAYAAABNtpsqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXyU1dnw8d81M9kXQhJICAkJBBJkVxEUUSiL1K22BcQVl1pa9bX1fapWbeujrdW21mp9Wqz6vFbFaq0Wxb2uKLKjIsgWEsgCSchK9pnMct4/ZjImyBLAZOYm1/fzuT/OzH3OfZ/EcOXkus8ixhiUUkpZiy3UDVBKKXX0NHgrpZQFafBWSikL0uCtlFIWpMFbKaUsSIO3UkpZkAZvFXIi8lsRWSEiL4lI7JHOi8gZIrJaRD4SkedFJCIU7VYqlDR4q14lIssPeD8GyDXGnAW8B1zbjfMlwAxjzDRgF3BRLzRdqbCiwVuF2lnAW4HXbwFTj3TeGFNujGkLfOYBfD3eSqXCjAbvPkxEkkXkZRFpEZESEbnsMGVPEpEPRKRBRApF5HsHnM8RkTdFpF5EKkXkLyLi6EYz+gMNgdcNQHJ3z4vIUOBc4PVu3EepE4oG777tr0A7kAZcDjwqIqMPLBQIwsvwB8lkYBHwrIjkdSq2GKgCBgETgGnADYH6Q0RkeSBlMqHjtYgMAeqBfoFr9APqDrj9Qc+LSCLwNHClMab9mL8DSlmUBu8eJiLFInKriGwK9HD/n4ikichbItIkIu+JSP8QtCsOmAv8yhjTbIz5BHgVuPIgxUcCGcBDxhivMeYDYOUBZYcC/zLGOI0xlcDbwGgAY0ypMWa6MWY6sLHjtTGmFPgEmBO4xpzAdTv72vnAL5PngbuNMTuO49uglGVp8O4dc4HZQB5wIf7c7Z1AKv7/Bz85WCUReV1E9h/iON5UQR7gNcYUdPrsCwIB98CmHOKzMZ3e/xm4JDAaZDD+dMbbR2qEMWYzUCIiK/AH5ydFJF1E7jnUeeBSYDJwV6AHv+BI91HqRNOdnKQ6fv9jjNkHEAhCVcaYzwPvXwZmHqySMeaCHmxTPF/lkjs0AAkHKbsdf0rkVhF5CPgW/rTIh53KfAT8EGgE7PhTGq8ceKFA7/vAz+444KMW4L8Pc35J4FCqz9Ked+/Y1+l120Hex3+TNwv0Rs0hjk8CxZqBxAOqJgJNB17PGOMGvgucD1QCPwP+BewJ3M8G/AdYCsTh/4uiP/D7b/LrUkp9RYN3GAvkxZsPcbx1qHqBfLIc4ugYilcAOERkRKeq44Eth7jmJmPMNGNMijFmDjAMWBc4nQxkAX8xxriMMbXA34Hzuvl16iQdpY6SBu8wZow51xgTf4jj3OO8dgv+nvKvRSRORM7EP9nloOkIERknItGB4HkL/lElTwWuVQPsBq4XEYeIJAFX4c+hH5ZO0lHq2Gjw7ttuAGLw57OfB643xgR73oGe/52Bt1cCFYGyM4HZxhhXp2t9H/g2UA0U4p8883+70QadpKPUMRDdBk2FUuCXw1ZjzCuBHvu/jDHndOd8YJLOC/gDuo71Vn2KjjZRPU5E0oGXDnJqHjpJR6ljomkT1eOMMZXGmKkHOSrRSTpKHRMN3iqkdJKOUsdGc95KKWVB2vNWSikL0uCtlFIW1OdGm6SmppqcnJxQN+OYFRUVAZCbm3vcZffUtgCQmRJ3VNftyz799NMaY8yAULdDqT4XvHNyctiwYUOomxEWbn1mNQAPLDwjxC2xDhEpCXUblAJNmyillCVp8LaYO+64gzvuOHCF1N4tq5QKvT6XNrG61atXh7ysUir0tOetlFIWpMFbKaUsSIO3UkpZkOa8LSYzMzPkZZVSoafB22KeffbZkJdVSoWepk2UUsqCNHhbzM0338zNN98c0rJKqdDTtInFbNy4sVvlcnJyKCnxz+T+85//3OVcdnY2xcXFx3RdpVR40OB9giopKWHatGkALF++vMs5EQlBi5RS36Q+kTYRkUUiskFENlRXV4e6OUopddz6RPA2xjxujJlojJk4YICu5qmUsj5Nm1hMXl5eyMsqpUJPg7fFPP744yEvq5QKvT6RNlFKqRONBm+LWbRoEYsWLQppWaVU6GnaxGIKCgpCXlYpFXra81ZKKQvS4K2UUhakwVsppSxIc94WM2HChJCXVUqFnhhjQt2GXjVx4kSzYcOGUDejx4kIh/p/23Hu1mf8mw4/sPCM3myapYnIp8aYiaFuh1KaNlFKKQvS4G0xV1xxBVdccQXgX/ZVRA56ZGdndyl7NNdVSoU/zXlbzJ49e4KvS0pKDpkaAZg+ffoxXVcpFf60562UUhakwVsppSxIg7dSSlmQ5rwt5owzuj+sr6fKKqVCT8d5W9jhxnJ3p56O8z56Os5bhQtNmyillAVp8LaYuXPnMnfu3JCWVUqFnua8Laa2tjbkZZVSoac9b6WUsiAN3kopZUEavJVSyoI0520xM2fODHlZpVTo6ThvC9Nx3r1Px3mrcKFpE6WUsiAN3hZz7rnncu6554a0rFIq9PpEzltEFgGLAIYMGRLi1hyftra24y6bnZ2NiDDxmvsBkKumBM/l5ORQXFx8XG1USvW8PtHzNsY8boyZaIyZOGDAgFA3J+SKi4sxxjB9+jSmT5+GMYZp06Yxbdo0SkpKQt08pVQ39IngrZRSJxoN3kopZUF9Iud9Irngggt6tOxHH3101G1SSvU+HedtYcc6zrvDwcZ5H+81T3Q6zluFC02bKKWUBWnaxGKmT58OwPLly3ukrFLKGrTnrZRSFqTBWymlLEiDt1JKWZAGb6WUsiB9YGkxF198cY+W1XHeSlmDjvO2MB3n3ft0nLcKF9rztpjW1lYAYmNje6SsUsoaNHhbzHnnnQd0b+z2sZRVSlmDPrBUSikL0uCtlFIWpMFbKaUsSIO3UkpZkD6wtJirr766R8vqOG+lrEHHeVuYjvPufTrOW4UL7XlbTE1NDQCpqak9UlYpZQ0avC1m3rx5QPfGbh9LWaWUNegDS6WUsiAN3kopZUGaNglzOTk5lJSUfO1zESE7OzsELVJKhQMN3mGupKSky+iPo9mXUil14tLgbTHXX399j5bVcd5KWUOfGOctIouARQBDhgw59WBpiHDVk+OudZz30dNx3ipc9ImetzHmceBx8E/SCXFzjktZWRkAWVlZPVJWKWUNfSJ4n0iuvPJKoHs572Mpq5SyBh0qqJRSFqTBWymlLEiDt1JKWZAGb6WUsiB9YGkxP/vZz3q0rI7zVsoaNHhbzIUXXhjyskqp0NPgbTE7duwAID8/v0fKKqWsQYO3xfzoRz8Cujd2+1jKKqWsQR9YKqWUBWnwVkopC9LgrZRSFqTBWymlLEgfWFrML3/5yx4tq+O8lbIGDd4WM2vWrB4tm52djYgcslx2djbFxcXdvq5Sqmdo8LaYjRs3AjBhwoQeKXukwHy4wK6U6j0avC3m5ptvBro3drunyiqlQk8fWCqllAVp8FZKKQvS4K2UUhakwVsppSxIH1hazH333Rfyskqp0NPgbTFTpkwJeVmlVOhp2sRiVq1axapVq0JaVikVetrztpg777wT6N547J4qq5QKPe15K6WUBWnwVkopC9LgrZRSFtQngreILBKRDSKyobq6OtTNUUqp49YnHlgaYx4HHgeYOHGiCXFzjsvDDz8c8rJKqdDrE8H7RNKd5V17uqxSKvT6RNrkRPLee+/x3nvvhbSsUir0tOdtMffeey/QvV1yeqqsUir0tOetlFIWpME7xHJychCRQx7Z2dmhbqJSKgxp2iTESkpKMMbSA2CUUiGgPW+llLIg7XlbzGOPPRbyskqp0NPgbTH5+fkhL6uUCj1Nm1jMa6+9xmuvvRbSskqp0NOet8U8+OCDAFx44YUhK6uUCj3teSullAVp8FZKKQvS4K2OSnZ29iEnFOXk5IS6eUr1GZrzVkeluLj4kOdEpPcaolQfp8HbYpYsWRLyskqp0NPgbTFZWVkhL6uUCj3NeVvMCy+8wAsvvBDSskqp0NOet8U8+uijACxYsCBkZZVSoac9b6WUsiAN3kopZUEavJVSyoI0ePeCw+2WozvlKKWOhT6w7AXf5G45L730UsjLKqVCT4O3xaSmpoa87KF0TJ0/1LnDzc5USh0dDd4W89RTTwFw9dVXh6zsoejUeaV6j+a8Leapp54KBtpQlVVKhV6f6HmLyCJgEUC/fv16vReoDyUPn1JR6niIyG+Bs4F9wEJjTGunc/2Ad4FRwOnGmC9F5AzgT0A7UB6o4+79lh+fPtHzNsY8boyZaIyZOHz4cIwxvXportefUunt73tPHCq0RGT5Ae/HALnGmLOA94BrD6jSCpwPdH4iXwLMMMZMA3YBF/VYg3tQnwjeSvU1IlIsIrNC3Y5ecBbwVuD1W8DUzieNMW5jTPUBn5UbY9oCbz2Ar8db2QM0eCuljouIjBARp4g8e5gyz4pIhYg0ikiBiFx3NOcPoz/QEHjdACQfRbuHAucCr3e3TjjpEznvE8mbb74Z8rJKHeCvwPojlLkf+IExxiUiI4HlIvK5MebTI50XkSHAM4FyEzqlThYC9UC/wPt+QF13GiwiicDTwJXGmPbu1Ak32vO2mNjYWGJjY0NaVlnGaSKyVUTqReTvIhL9Td9ARC4B9gPvH66cMWaLMcbV8TZw5HbnvDGm1Bgz3RgzHdjY8doYUwp8AswJ1JsDrOxGmx3A88Ddxpgd3ftKw48Gb4tZvHgxixcvDmlZZRmX4w9ouUAe8MuDFRKR10Vk/yGOQ6YUAr3XXwM/605jRGSxiLQC24EK4M2jOX8wxpjNQImIrAh8rU+KSLqI3NPpum8C5wBPiMjVwKXAZOAuEVkuIpZcB1n62hP0iRMnmg0bNoS6Gcds+vTpACxfvvy4y976zGoAHlh4xlFdty8TkU+NMRND3Y4jEZFi4HfGmL8F3p8H/I8xJvewFY/uHn8Gyo0xvxeRu4HhxpgrjlDHDpwBTAd+f+AQvSOdV1/RnrdSJ66yTq9LgIxv6sIiMgGYBTx0NPWMMV5jzCdAJnD90Z5XX9HgrdSJq/PGpEPwT0j5GhF5S0SaD3G8dbA6+HvGOUCpiFQCtwBzReSzbrbNQaec9zGc79z+34rIChF5SURiDzjXT0TWBb6WMYHPzhCR1SLykYg8LyIR3WxzWNHgrdSJ60YRyRSRZOBO4KCblBpjzjXGxB/iOPcQ134cf3CdEDj+BrzBVw8Pg0RkoIhcIiLxImIXkTn4884fdOf84egkHaXUieg54B38AWoXcO83dWFjTKsxprLjAJoBZ8eEmEBv/s6O4vhTIHvwD+37I3CzMWZZN88fTp+dpNPnHliKSDX+37y9KRWo6eV7dle4ti1c25VvjEkIdSOUX+AXxFZjzCsikgT8yxhzzkHKPQX80RjzZafPhuL/a2SqFcd697lJOsaYAb19TxHZEK4jFMK1beHcrlC3oa8RkXS6pj06zKMPT9Lpc8FbKWUtgbTM1IOdE5FPgDvwB2KdpKOUUlbQlyfpaM+7dzwe6gYcRri2TdulusUYc8cBH7UA/93p/HkHqbakRxvVC/rcA0ullDoRaNpEKaUsSIO3UkpZUJ/LeaemppqcnJwevYfP58Pr9WK327HZbF0+7/w+1PbUtgCQmRIX4pZYx6efflpzuOGmvfHzFe6KiooAyM098uz2g/0MHk39471/uDnSz1dnfS545+Tk0NOrCvp8PlpbW4mNjcVms+HxeKipqSE1NRWHI3y+5Z1XFVTdIyKHneDVGz9fJxL9GezqSD9fnYVPJLEgj8dDVVUVsbGxJCYmBnvVNpuN+Pj4YLmamhrKy8vxeDw4HI6wC+JKKesJn7/hj5KISCju6/P5aG5uxufzUVNTw65duyguLqa1tbVLOY/HQ2VlJR6Ph9TUVDIy/KtxlpeXU1MTjrO+lfpm3HHHHdxxx4Gj96xT3yos1/0TkdGBLZNCMsaxtbWVxsZGAFJTU/H5fAfdQqyjtw2Qnp5Oenp6l573oRyYclHKalavXm3p+lZhqegQWCryHyIyIlRt6EiRxMbG4nA4yMjIICkpCY/HQ1FREe3t/mUSOnrbnQO1w+EgPT39sCmTjl8OB/bklVKqM8sEbxH5Dv5ZUzcaY3YeTdpERBaJyAYR2VBdXX3kCofRkc8+sFdcVlZGQUEBW7dupbGxEZ/Pv8qk0+lk06ZNOJ1OjDFdjs7nOnT+5aCs4Zv8+VKqu8I+bSIiHVHyXqDBGLNSRNKAK0UkFf+SjoXGmKZDXcMY8ziBac0TJ07skXRLVlYW9fX1lJaWUlFRQU5ODm1tbezfv5/KykoAxo4d26VOQUEBX37pX6Fy3LhxwNcfdqrw1xs/X0odKOyDNzDQGFMpImcCq0XkJfxLP74P9Ad+CjxBN1YTOx5HGu4XGRlJv379WL16NW63myFDhpCRkUFeXh6FhYVkZGTgdrupqqoC/HnwvLw8gOB/lToRZGZmWrq+VYR18BaRbwOviMj3jDFvichkYCPwnDHmd4Ey9wML6eHgfeADyM46RpY4HA4mTZqE0+lk0KBBJCcnA5CRkcE777xDTU0NGRkZeDwepk6dyuDBgxkzZgytra1hN4FHqWP17LPPWrq+VYRt8A4E7nuBl4GxIvKhMaZFRE4CvCIigREnhcAoEbEbY7w91Z6OB48HGylSU1PDpk2baGxsZMqUKSQnJ9Pe3s5LL73ErFmzKCws5L333qOtrY3JkyeTl5dHU1MTTU3+TE9jYyONjY3Ex8cTHR1NXV2djgVXSh1WWEYHEZkG/A74MeAGFgNP4t8jz9Op3HX4975b2JOBG74aKXIwqampDBw4MPjgMT4+nn/+85+8/fbbtLa2cvLJJzN//nyam5uZOXMm+/bto7CwkN27dzNu3Di2b99OXFwcycnJOBwO6uvraWlpIS0tjS+//JJx48bpA0xlGTfffDMADz/8sCXrW0VYBm8gB1hkjFkHICJbgT+JyLXGGI+I2IE84DvA1caYLaFrqj+wDx8+nObmZt577z3OOecchg8fTkREBA6Hg+LiYgYMGMDs2bPxeDwUFhZSUFCAMYaqqira2tpISkoiOTmZwYMH43Q62blzJ9u2beOzzz5j3bp1XHvttfogU1nCxo0bj6leTk4OJSVfzQ7/85//HHydnZ1NcXFxj97fasIyeBtjngb/dkWBnvb/ADcCGUCpMcYrIkXAZcaY5hA2NSgxMZGysjJWrVqFx+MhLS2NhIQE2traiIyMxOv1UlpayieffEJ5eTkOh4OUlBRSUlJoa2tj8ODBtLS00N7ejtvtpry8nIyMDOrr62lubmbdunXMmDEj1F+mUj2mpKQEYwzTp08HYPny5cFzIZpQHdbCLnh3zl13SpFsAwYDNwC3B861A726cejhRpyICHPmzKG+vp5Zs2aRnJyMzWYjMzOT9vZ2fD4fy5Yto7q6mpaWFrKysoiOjsbj8ZCcnExCQgL79+9n//795OXlYbPZEBFOPfVU2tvbSU9PZ8eOHQwdOpSIiIjgPZVSfVNYBG8RmQoMNcYsCfSqbcYYX+CczRjTJiI3Ac+LyKSOdEpvO9SIk47RJsXFxQwfPhxjDA6Hg9GjR1NaWkpUVBSFhYXs3LmT1tZWBg8ejNvtprKykvr6eiIjI4mLiwuODR89ejQjR46ksrKSESNG8Omnn7Jx40ZEBIfDwbBhw0Lx5SulwkhIg3dgAk4s8Jj/rcQZY/5mjPF1pEwCr21AHfAWsDtU7T1wxInH46G8vJzS0lIaGxsxxpCSksLgwYNpb2+nrKyMDz/8kJSUFEQEu92Ow+Fgz5497Nu3j7S0NJxOJw0NDbS2tjJp0iQyMjKCPe3q6mpWrlzJxx9/zMyZM7nooouor6/H5XIRFRUVqm+DUod1vPMWQl3fKkIavAO962YReRrwAlNEJMYY81DnUSWBcrUicp8xpi1U7T1wxElNTQ0ffvgh69evZ9asWfTv35/U1NTgqoExMTGkpKSQlJREXV0dGRkZbNy4kR07dlBWVobD4SAmJobmZn/a3ufzcdpppwVHnOzatYs1a9ZQWVmJ1+ulra2N4uJiIiMjvzZbU6lw8fjjx7dHc6jrW0VYpE0ADzAEeBq4TkT+BLiMMXeIyBmB158BzsNdpLelpqbS0tJCbW0tu3btYtCgQVRUVJCbm0tUVBQjR47E5/Oxe/duCgsL8Xq9VFdX09jYSENDA3a7nba2NlpbW4MjSRoaGvB4POzYsYOamhrcbjeDBg2if//+ZGZmYrPZGD58eIi/cqVUqIVL8F4GzDfGvC8iE/BPznkqcG4S/vVLCNUysIficDi45JJLSEhIYPbs2dTW1tLQ0EBzczO1tbXs3r2bjIwMqqurKSwsZMqUKTQ2NrJt2zYAvF4vLpcLAGMMe/fuZe/evRQVFZGcnMywYcM47bTTcLvdAGzevJmEhAScTicxMTEh+7qVOpxFixYBx94DDnV9qwiX4N0G5IvID/FPzLkfmCQilwGPhFvQ7iwyMpKZM2cSGxtLS0tLcGx2amoq7e3t7N69my+++CI4CzMtLY24uDiqq6uDeXCv14uIUFFRgTGGmJgYEhISaG1tZePGjaSlpbF7924aGhpISEggJSWFyMhIoqOjcTqduva3CisFBQWWrm8VYRG8jTHlIlIG/Ar/kq+vici38K8WGLaBGwjOfPR4PLS3t1NcXMz69esZOnQo48ePZ8SIEfTr14/q6mpGjhzJnj17iIuLIy4ujpiYGOLj46mvrwcgLs6/CeuAAQNwuVxs3ryZ4uJiqqqqSE1NZcOGDWRmZjJ16lQqKiqIiIigvLycYcOG0dDQQExMDBkZGdjtdh1GqE4o2dnZh/yZPpoJPCeSsAjeAU8Ay4wxnwbef9QxXDCcdSzh2tjYGNycwefzUVlZGVzjZO3atVRWVvLFF1/Q1NREaWkp7e3teDz+Z7Je71cz+6Oiomhvb6eqqor+/fvjcDiIjY3F6/VSUlKCy+Vi69atnHXWWdTV1VFVVUVrayvV1dVERUUxfPhwRo0apeuiqBPK4YJzX+2ohM2/cGNMGVDWseCUFQJ3Zx0bNMyYMYPi4mJEBI/Hw9q1a2loaCAxMZGoqCiKioqCwbpjw4bOXC5XMCfe1NSEMYY9e/YQHR1NbGwsERERjBw5ks2bN5OUlER2djZ2u53k5GRKSkpobm6mpqbmkOuwKKVODGETvDuEe5rkUESEuLi44OYKKSkpxMTEsHfvXgYPHkxTUxMrVqygvr6+2z0Ft9tNQ0MDAG1tbbhcLtxuN6+99hoFBQWICOPHj2f06NG0tbVRXV0dHK6oVKhMmDDhkOcOXL9k4jX3AyBXTSE7O/uI9Y/3/ieSsAveVuN0Otm+fTuZmZkkJyczZswYGhsbSU1NZd26dURFRZGbm8uqVaswxgSD8bFobW1l586dOBwO2traqKqqIi4ujgEDBtDa2kpdXV2XLdWUCoXDrebXsX5Jh1uf8W8WvP7Jrz473tUAT/TVBDvoEIXjVFBQwGeffcbGjRuDmwaPGTOGhIQEEhMTGThwIMXFxWzZsqVLj+NY83R1dXWsWbOGXbt2ERUVRVZWFqeffjqZmZkkJCQQFRVFTU0NPp8vuE74wdIzSilr0573ccrLy6O1tRWXy4XNZiM6Ohog+JAxOTmZ6upqDtyY9lizQ263G7fbjcvlIjMzk4iICNavX4/dbicuLg632010dDTNzc3s27cP8D9UjY2NpbW1VYcV9gARWQQsAhgyZEiIWxN6V1xxBXDsO9qEur5VaPA+Th0PEnft2kX//v0ZO3YscXFxrFixghUrVuByuRg5ciRJSUnBnnlnI0aMCO6oc6Bhw4YdcsxqbGwsO3fupLa2lpEjRzJ58mRSUlIYOHAgzc3NJCYmkpaWFizb2tpKY2Mj4B+S2Fef0PcE3YC4qz179li6vlVo8P4GdN5I2OVyBR8m1tTUEBERQUlJCXa7Pbg0bGe1tbXBWZYH2r9/P+3tB1/11uFw4HK5qK2tpa6ujtbWVuLi4mhsbCQxMZHo6Gja29uJjY1FRILj0XVHHqVODBq8j5PP58Pj8TBmzBhsNhubNm3i008/pbKykjFjxhATE0NdXR0tLS1HzD0bY3C73Xg8HjweD7t376a2tpaTTz6ZsrIyUlJSgr8QOnSMCd+7dy/9+vWjpqaGuLg4KisrycjIIDU1lejo6OA65NrjVurEoMH7OHWM8ti3bx/Z2dnk5eVRUFCA0+kkOzubAQMGsHz5cgYNGkRDQwMtLS0HvY4xhoyMDFauXBn8rK6ujkGDBuF2u5k4cSLvv/8+DocjOHGnw759+9i8eTOFhYUMGzaM1NRURo4cCUB5eTmDBg0KzuLsSKUopaxNg/dxio2NZd++fdTX1xMTE0NiYiLnnXcehYWFtLe387//+79UV1cHd9Spra0NrmHSmdfrZeXKlcyfP58f/vCHZGVlkZqaGpwyD7Bjxw4uvfRStm/fzplnnsnOnTux2Wz4fD4KCgqw2+20tLQQFxcXHJbY0NDAlClTmDBhgo7/Vr3ijDPOsHR9q9DgfZxsNhvZ2dnExcURHR0dzDmPGTOG//znP0RHRzNp0iQiIiLIzc1l/fr1VFdXB1cK7NDxfuHChZx99tkAXxuznZ+fz+rVq/njH//IAw88QFxcHK2trTgcDkQkuE9mdHQ0SUlJuN1uysrKsNvtFBUVceGFF5KcnEx8fLymT1SPuf/++y1d3yp0zNg3oGOThsTERBITE4OjOzIyMpgyZQrz5s1j1qxZTJgwgc2bN38tcMNXwftIu4BERUXxi1/8gtWrVzN69GhcLheTJk0K5tPb29tpbGykoqKCxsZG6urqePfdd1m6dCkvvPAC+/btO+ioF6WUtWjP+xvUsUiVx+OhsbGRrKwssrOzaW5uxuVyUVZWRmJiYpe897Bhw2hra2Pnzp3ExMQwYMCAYI+7uLgYu91+yHu98sorPPPMM9x99934fD6ys7NJSkoC4IsvvsButzNkyBDy8vKIjY1lzJgxuN3uLg88lfqmzZ07F4B///vflqxvFZYN3h0LWIW6HQfTsVFx//79GTBgAG1tbSfvgFEAACAASURBVKxevZrm5ubg2g4di1Pt3r0bj8dDbW0tI0aMCE7yAf9DzI70hsvlCj7A7Dhns9m4+uqrmTNnDhdffDG7d+/G7XYHd+XZsmVLMP89YsQIioqKgrvV68JVqqfU1tZaur5VWC5tIiIREN4LWHU8aCwtLaWkpIRt27axefNmNm7cCEBiYuLX6gwcOJD8/PyDXq+hoYGrrrqKb3/72/zlL3/5Wtpl0KBBLF26lMGDB1NZWdmlZ79u3Tr++c9/8uKLL+L1etm2bRt79+6lvb0dY0yXpWmVUtZhqeAtIt8BHhGRp0VktIgkh7pNB9MxnK+4uJhNmzYRExNDTk4OaWlptLS0BFMbHXw+HyUlJYwYMeJr12ppaeEHP/gBO3bsYOrUqfztb3/jkksuYeXKlV3GjQ8YMICXX36Z1NRUqqqqupyrqqpi69atvPLKK3z22WcsXbqUzZs3U1NTQ1tbW3DDZKWUdVgmeIvIGOBR4EVgH3AjsFBEBnej7iIR2SAiGw5cY6SnZGdnM2XKFMaNG4fT6cRms+F0OikrK6OxsZGEhIQu5Ts2Iz7QRx99xPbt2/nDH/7AX//6Vx555BHq6uq46aabmDhxIo888kjwz8SMjAyWLVtGfHw88fHxXQJ4dXU1e/fupbq6GpvNxqeffkpDQwNOp0uHECplQZYJ3kAa8LEx5gNjzG3Am0AmMFdEkg5X0RjzuDFmojFm4oABA3qjrURGRpKfn092djZer5eWlpbg7MmYmJguDyJtNhszZ87kxRdf/NoszIKCAhwOB9OmTQNgxowZvPPOO9x///1kZGRw9913M2bMGG688cbgQ9JXXnmFiIgIkpKSgtfzer3s37+furo6CgoKaG5uDk720V131Ddp5syZzJw507L1rSLs/9WKSLQxxglsAO4Rke8ZY142xrwuIjbgAiAV2B/Shh5CfHw8p59+OiJCS0sLY8aMISMjg5qaGt5+++1guS+++ILKykpWrVrF1KlTg5/v2LGDoUOHEhkZGfwsIiKCOXPm8KMf/Yht27bx97//naeeeopdu3bx4osvkpuby7vvvsv06dNJTU2lpqYGl8tFZGQkLpeL1tZW7HY7MTGxREZGUFRURFZWVpd7KHWsfvWrX1m6vlWEdfAWkRnASSLypDGmQUSeBqaISJ0x5iNjzKsiMgu4Gfg/oW3t13k8HmpqavB4PFRUVFBfX8+pp54aXPVPREhPT8ftdpOamkpTUxPPP/88p59+OgBFRUVs3bqVUaNGBbdG61BdXc2WLVsAmDZtGtHR0SxevJjZs2dzyy23MG/ePJYuXcp3vvMdBg4cSHJyMlFRUURGRpKens7+/fuxJ9hxOl0UFRUB/mGL4F9r3Ol0UlBQQF5eXpcRMEqp8BC2aRMR+TbwMLDJGNORDP4AaAC+IyIXBz77AnAGeuFhpaamhj179tDS0sLo0aMZOnQoSUlJVFVVsWbNGtxuNzt37mT37t2Ulpbi9XpZv3490dHRREdHY4yhtraW/Px8+vXr1+Worq6mra0teIwbN45rr72WgoICHnzwQVpbWxk7dixLly6lvr6elpYWmpubaWxsDC5e5d8T08GXX37ZZRq+z+dj48aNbNq06ZBL0ip1KOeeey7nnnuuZetbRVj2vEVkHPA88ENjzAoRSQUEf+B+ALgcuE1EFgCTgPPDccPixMTE4HT5ji3SUlJSqK+vx+VyBcdjd3A4HBQWFuLz+bDZbOzatQv4qkd8JB0zLf/+979zxRVX8Nxzz3HyySfz/PPPM2/ePPLz8/F6vcFZlxGnphAZEUnR2ndISEjguuuuA/yLbQ0cOBA48oxPpQ50sAfvVqpvFWHXWw2IBv4FpIvIROBZ4E/AauB7xpgngW8BfwAmG2M2haylh1FRURGcXNOxgfDevXuDe07abDZEJDgRx+Fw4HQ62bt3L+CfwAPdD94Ap59+OgsXLuTDDz/kyiuvxOl0cuaZZ/L000+zZcsWEhISaG5uZuXKldTW1FBVtY+KigqcTicVFRVs2rQJm83GwIEDmTRpkqZMlApTYdnzNsasC0zGmQ/cC9yJf6eSM4DnRWSbMeYLYG0Im3lEWVlZAAwePBin00l6ejoNDQ1s2rSJzMxMoqKiMMYEVxjsGPWxbds2srKyqKurC065P5y2tjY+++yz4AJYZ555JieffDI333wzCxYs4Nlnn+Wcc87hb3/7G9ddd11wmVmfz4cx/hlp69evZ/v27YwYMQIRYdy4cT37zVFKHZew63mLiAPAGLMSeAH4gTFmMeA1xqwA3gYOvih2mImMjCQ3N5eoqCj69evH4MGDaWhoCO4zGR8f32WYXmRkJCLChg0bABg9ejQ+n48vv/zysPdZunQpTz/9NH/+85+Di05dddVVPProo6xcuZKLLrqI+vp65s6dy9y5c9m/f3/wF4YxPux2O5WVlRQVFVFTU3NUPX2lVGiERc9bRKYCQ40xS4wxHhFxGGM8xpjVIhIJ/unwgRz3yYBlk1oTJkzA5/NRX1/Pli1baG5uDo5Ksdls5Ofns27dOgDGjx9PREQE69at45RTTjno9SoqKlixYgW5ubkUFRXxhz/8gZ/85CcAXHLJJfTr14+rr76a6667jn/961/cddddvPbaa8G1VQAqKyuJi4vj9NNPJz8/n4qKCnJzc3Wj4m7SDYi7uuCCCyxd3ypC+q9TRGwiEg88BtwhIj8GCATwjjVM2kXEISILgV8AC40xe0PX6uMTGRnJsGHDiImJYejQoQwbNoyBAwcGUyOlpaVs2LABYwzR0dGMHz8+GMwPZunSpURGRnLDDTfw05/+lLq6On73u9+xefNmwP/k/YEHHuCDDz7g3nvvZciQISxatCiQrvnqOu3t7dTV1WGMCY4FV90Tiklg4eyWW27hlltusWx9qwhpzzswQqQ5MH7bi38Md4wx5iFjjLtTOY+INAHfN8YUhqq9x6rzxgc1NTXBUSRjxoyhqKiIurq64DrcPp8vuCdlVlYW55xzDvfddx8+n69Lr+7FF1/E6XRSWFjIoEGDWLJkCeCflm8vKuIfM2aQPn8+/SdPBvy7izz88MM0NDQwbtw4+vXrR//+/h3t09PT8Xq97Nmzh6KiIiZPnkxsbGwwH9/a2hrcyFg3cVAqPITL38UeIAt4GpgkIn8SkfsBRGSKiJwUmFVpucDdobGxkTfeeCPY887Ly2PChAmceuqpjBs3jpSUFIDgtPm1a/3PYs866ywAVqxY0eV6xhjKy8uJiIigc28vX4QvPB4e9Xj45fPPU/b++4B/jePs7Gyee+45GhoauOOOO6irq8fr9Y+w9Pl8OJ1Oqqur2b59O8XFxXi9XlpbW2lsbNSeuOq26dOnM336dMvWt4pwCd7LgEpjzPv4p8FfD3SsV3IaUB+qhn1TVqxYwccff8zq1avJyMggMzOTzMzM4J6TLpcL8K9zEhUVxfr16wF/TzonJ4ePP/64y/U6Amp6enqX3PR9u3djxz8o3g5c9dprvPvuu9jtdq655hoiIiJYvHgxl19+OdHR0URFRXYZ7VJRUcEzzzzDsmXL2L17N9HR0cHdgZQ6kpycHD766CM++uij4F9qnY/s7Oxv/J7Z2dld7tH5/jk5Od/4/cJFuATvNiBfRH4I/Bi4H8gSkcuAR4wxlSFt3TfgrLPO4uyzzw72pG02G1VVVRQXF7N8+fJgz1ZEOOWUU3j11VeDkw3OPvts1qxZExz3Df6d5SMiIkhO7roqbn7gl0CH0UDBG2+watUqkpKSuOaaa6ipqeHOO+8kJyeH5uaW4PrfTU1NFBcXs2/fPkpKSqioqKCqqiqYMlHqSEpKSpg2bRrTpk0Lpt06H8XFxd/4PYuLi7vco/P9S0pKvvH7hYuwCN7GmHKgDPgV8F/GmF8DDwIrwnnThaORmJjI+eefH9yIYf/+/Wzbti04TT0mJob+/fsD/kWqCgsLeeCBBwC45ppriI2N5fbbbw/2kmNiYoKrFHbWeZqpBI5/REbyxhtv0NzcTG5uLrNmzWLJkiXExsaQlNQPp9MZ3JShvr6ebdu2UVpayrp166ivr9eUiVJhKCyCd8AT+GdPvhZ4/5ExpiyUDepJ7733Hp9//jlOpxOv10u/fv2C6Q+Hw0FkZCQvvPAC77//Pmlpadx6661s3ryZ1atXAwQDfV1dXZfrvhHY6KHzb7wJbjdOp5PXX38dgPPOO4+UlBR27drFyJEjSUxM7LL+d319PXv37qWwsJClS5fy+uuv8/777/eZacdKWUFYjPMGCATqso69KcNxrZJvisfjYezYsXg8HjIzM+nXrx/GGNauXcvq1atxu91ERUWRm5vLL3/5S8aMGcMFF1zAQw89xDPPPMOUKVOIiooiLi6Ouro6Bg4cGExr/Donh/M3bqTzFsPGZmP6tGl88MEHnH766YwaNYq77rqLZz5voKmpmeeee47Zs2dzwQUX8PnnnwOwdetW9u7dS0xMDIMGDWL8+PHB5W2VOpyLL774yIXCuL5VhE3w7nCipEkOp6amhpaWFqZPn07//v3JyMggLS2Np556KriF2Z49exg6dCi7d+/mtttu4/HHH+fiiy/m0UcfZefOnTidTmJjY6murqa6upqYmJjg9Xc5HOR3SqdU2mw0NzcTERHBE088wbx587jpppt4cfsKiooKSZh2FjfffDMPPvggJ598Mvv3+5dG7xgi6PP5gm1JT08nMjKSgQMH6iYO6qBuuOEGS9e3inBKm/QZqampZGRkkJqaSmRkJDk5ORhjiIuLIzIyErfbTUtLC9u3byc+Pp7169fzxhtvcPvttxMZGcmbb77JggULuPbaa4mMjCQ1NZXLLrsseLSedFKX+6W73dSvWxe87sqVK7Hb7QwenInT6eLZZ5/l8ssv55xzzuGLL77A7XZjjMHhcODz+airq+O9995j+fLlvPjiixQWFlJSUvK1XX+UAv8v/eN5ThLq+lahwTsEHA4H6enpOBwOXC4X69ato6qqiuTkZOLi4oJrb4N/J565c+dy9913U1payoIFC1iyZAltbW1ERkYyevRotm3bhtPpDF6/NjBmHPwPLAHO6nT/4uJi6urqghN1/vrXv7Jx40buuususrOzSUhICI779nq9GGNIT08nJSWFc845h7S0NCIiIigvL+ell14K9tSVAv8zlfPOO8+y9a1Cg3eIFRQUUFhYSFVVFQkJCYwaNYrZs2cHZ1OKCOvXr2fQoEFcc801XH/99bS2tgbHgZ988sm43e5grhpg66hRuPA/tDRABP6hPB28Xi+LFy9GBLKzcxg0aBC33XYbNpstuJHDmWeeCfgfXg4ePJjs7GwmTpxIaWkpKSkpJCcns3z5cl5++WXefPPNXvleKaW+osE7xPLy8hg3bhwTJkxg8uTJTJo0iZNPPpl58+YFJ8bY7XZcLheFhYWUlZUxf/58Vq5cSXl5ORkZGQwfPpzly5ezb98+wN/z/mvg+h0974md7jlo0CCWLVuGy9WOw2Hn17/+NTU1Nbz00ksMGzaM22+/nXfeeQefz4fP56O1tZX4+Hjq6+upqqqitLSU2NhYRo0axdChQxk5cmSvfb+UUn4avEMsKiqKsWPHEhUVRUxMDFOmTCEhIYExY8YwadIkoqKiAIiNjSUrK4vFixfzyCOPEBcXx8svv4zb7ebCCy8kOjqaf//737S3twNQe8B9Ujq9zsrKQkSoqKwA/CsdnnbaaTzzzDO4XC4uv/xyANxu//IyJSUlrFmzhrVr11JWVobb7WbdunUMGzaMiy++mFGjRvXo90gp9XUavMNMRkYGU6ZMISMjg9GjRwdnZIoIzc3NfPTRR+zdu5eLLrqIuro63nnnHeLj4/ne975HXV0db775JsYYXsI/bbVj6M6PgRGB14mz/i/f+tXL+OIHsau2nZ+/tpf+F/43o695iFdffZXs7GzOPvtshg4dijGGpqYmtmzZQkFBAc3Nzaxdu5atW7eyfft2UlNTaW9vp7m5WR9gKtWLdKxXmLHZbLjd7mCv1+FwcMopp7BlyxZSU1Npa2tj8eLFpKenM2nSJNatW0d2djYjRoxgypQprFy5koyMDHbi34LoLvx7yhn8Dy13AhJ18J15ohL68+jDj3LLLbcwf/58brrpJvLz84mLiyMuLo68vDySk5NJT08nLi6OcePGBVMrjY2NAEfc9Ued+K6++mpL17cK6QPDqrtIH3qSufKeJ0PdjEPyer3B4N3Y2MT+hv04HA6aGhtpbW3DGB8iNk499RTsdjubN3+Jy+UKbtywbds2mpqaSEhIYEhjDYuf+Q02nw+fzcYNC3/F3uQ0SBx8yPvX7d5MXl4eKSnJrF27FpvNTmSkf8pPdHQ0UdHRDEgdEEzzeL1eamtrSUpKIioqihN9CZQ/XjXlU2PMxEOdnzhxounYCakvEhGOJqbc+ox/xvADC88Ii/aEmogc9uerM02bhBGv1xfYWceLzWYnKiqKqMgoBIiKisbhsCMiwTW/RWyMGDECn89HUVEhA6vKuKp6J0Ob62huboJj+KGNioqkrKwMu93BgAED8Ho9wWVjwf8Pwe1243S20dDQQE1NDbV1tZSX78Xn8x722qpvqKmpoaamxrL1raLPpU0yU+J67Lf88VqzZg0rV25h+PDhnHbaabS3t7Ns2TJWrVoFTifS3s7y5cuZPHkya94pZPny5cTGxvJM60aevOMGbouIIMLh4FKfj1EuF/379QOfl0ifD5fxcdl/HuOF3FxKZ912yDYMa9nIkiVLWLRkCbf95FymTp0anM2ZlJREdnY28xYupKGhHrfbTUZ+Bv/610qSk5M5L3sQM2bM6MXvWO/741WhbkH4mzdvHgDLly+3ZH2r0J53GBk3bhynnnoqCQkJOBwOdu3axfnnn8/3vvc9ZsyYQVxcHAkJCaxatYqKigpuuukmfD4fV155JTeOHYvb7cbW1obNZmN+WhqvNTTgM4Z2ESKNYcbevfz1448ZXHPoXeTGjh3LoEGD+Pe//02/fv148sknKS8vp6GhgZaWFj7++GMee+wx1q5di9PpJCYmhgULFjBp0iQmTZrUi98tpfo2ywZvOQEXmI6NjSU5OZnKykqWLVvG2rVr2blzJzNmzCAtLY2hQ4dit9ux2+1kZ2fz9ttvc9999yEiXPD73yMitIhggKr8fJrS0/muMUggfWLDP2HnF8/ec9D7R5h2bDYbc+bM4YMPPsDtdnPaaacxbdo0IiIiMMZQW1vLp59+ypYtWyguLqaxsZGUlBS++93vBlcm7GsjT0RkkYhsEJEN1dXVoW6O6iMslzYRkdHGmC0n6gJWeXl5gH/I4Pr16xk+fDgVFRW0trZSWFgYXDY2IyODlpYWHn30UXJzc7n00kv57IknePZHP8I5fjzVgdEgI/bvx+500vk33ZC6cia717LVfhKNjQ28+6frue222wJjyvsxZ84cnnrqKdasWcNZZ53FT3/6U77//e8H61dXVwd3u+9YGVFEGDBgQJddffrKyBNjzOPA4+B/YBni5qg+olvBW0SigRuAqfhHnX0CPGqMcR624jdMROYAvxeR+caYnb15794SHR3NuHHjADjnnHNoa2sLjCrZjIiQnJxMYmIi/fv3Z/To0YwcOZLbb7+drKwspsyZw5B77uGXv/wlubm55OTkUHPSSXg+/7zLErGfJyVRXFyMc/BQXK52IiIiqKjwT9hJS0tj5syZOBwOli5dyvDhwxk7dizDhw/3j3ppagouP1tQUMDgwYMpLy/H6XQyevRosrKyyMjICG5g3OEE/ENJqZDqbs/7GaAJ+J/A+0uBJcD8nmjUwYjId4DbgRuNMTs71v3uZt1FwCKgyw7s4c5utxMfH4/H4+GMM85ARCgoKKCkpIRVq1ZRW1uLz+cjLy+PH//4x2zYsIHbb7+djRs38vLLL3PXXXdxyimnULx8OUNvugk7/p2eH8zJYcJJJ7HJJNDS0kz//v059dRTAf+Mz/79+3Paaafx8ccfc+eddwKwaNEibrvtNqKjo4NLwba0tFBcXExkZCQiwrBhwxg8eDBJgQ0hWlpadAu1Puj666+3dH2r6G7OO98Y8wNjzIeBYxGQ15MN6yyQ3/414DXGrBSRNOBnIvI7ETlZRBIOV98Y87gxZqIxZmLnndatoqamhoaGBnJzcxk5ciSpqamMHj2auLg4bDYbe/fuxev1snDhQrxeL0888QQZGRnceuut1NXVEVdXhy0mBhvgs9tJ2LiR3/72tzQ0NODxeILBtrNvfetbbN++nfLycgAuuugiUlNTmTZtWrCM3W6nsbGRqKgoUlJScLvdwen5zc3NVFRUBFdHVH3HggULWLBggWXrW0V3g/fnIhLcQkVEJgMre6ZJXxfoYU8FkkXkJeBZ/H819Ad+CozrrbaEQsf63+PGjcMYQ2VlJUVFRcGNg+12O16vl5UrV3LfffeRmJjIgw8+yP79+7nzzjtpnDABAG9MDBGRkQy+5BLS09Oprq7G6XTRr1+/r91z+vTpwFfDraKiorjmmmv4z3/+g9frH8/tdruJi4tjzJgxnHXWWQwcOJCGhgaWLVvGl19+GQzkqm8pKyujrOzYdzAMdX2rOGzaREQ289WqogtFpDTwPhvY2tONE5ExgBuIMsZsCvzS2Aj8wxjzu0CZ+4GF9OIvk97Wsf43wBlnnEFlZSX19fVUVlaSnJxMTU0NUVFRuN1u7rnnHvLz8xk7diw///nP+c1vfsNv/vlP/vuFF4j//HNaTjkF+2ef8X8mT+b9qiTq6+sZmp/f5X7GGF5++eXgvTtceumlPPDAA3i9Xux2OzabjbKyMqqrq/nHP/7BtGnT+PDDD9m/fz+ZmZlccMEFfeahpfrKlVdeCRz7OOtQ17eKI+W8L+iVVhyEiJwH/B5YDYwVkfuMMa+JyEmAp1POuxAYJSJ2Y8wJP8UvKyuLSy65hPr6ekpKSkhLS6OxsRGXy0VcXByTJ0/m6quv5rHHHuPiiy+mtLSUp59+mvb2dn7xi18QEREBn32GiBAfH098fDwT47r+Afbwww/z6KOPctlllzF//lePNZKTk8nKygo+3HS73dTX1/Piiy8yZMgQXnnlFcaOHUt+fj4XXHABI0eO1Hy3Uj3ksMHbGFPSWw3pTEQmAn8CrgLWAVcCs0TkLfx5bxModx1wPbCwLwRu8PeEMzIyuOSSS2hubsZms/Hll18C/hEdmzdvJjc3l5/85CcsWbKE//qv/yIyMpInnniC3bt38+CDDx72+m+99Ravvvoq8+fP59577/1a8B09ejTJyckUFRV1+dzlctHU1ERtbS2zZ89m5MiR2Gw2KioqiI2NJTExscswQqXU8QnXf02pwG+NMWsDgXoX/ry2Mcb4RMQR6IF/F7jaGLMllI3tSYea9JKTk8P5558fHBJos9mw2Ww4HA4qKyuJiori+uuvZ9++fdx4443cf//9bN26lYsvvphdu3bR0tKC1+vF6/XS0tJCS0sL//nPf3j11Vc577zzuPXWW6mvr6e2tjZ41NXVcdJJJ7F9+3bS0tJIT08PHuPHj6e9vZ2WlhZWrFjBJ598QmlpKWvWrGHNmjXBVQeVUt+MsAzexpi3gfc6fbQRaOvUu+5vjNkGXGKM2dzrDexFra2tNDY2dtlQtbW1lddee40dO3bQ1tZGVlZWcORJXV0dDQ0NiAhNTU3ceuutDBkyhFtuuYVPPvmE+Ph4nnvuOWw2G2lpA0lLG8hll11Ge3s7b731FrNnz+Yvf/kL6enpDBgwoMuRmprKqFGjcLlcNDY2BoN+S0sLFRUVGGMYMGAAmZmZ7Nmzh7KyMlwuV5/ZEFap3hR2MyxFxGGM8RhjKjp97AAyRcSOP4UyPzBR54Qfh9axFVrHfwE2bdrE9u3baW1tZebMmezatYuxY8ficDjYvHkzu3fvJiIiAq/Xy5dffsmCBQtYtmwZ48ePZ/Xq1cydO5c77riD837+/8jJyea5557jnnvuYc6cOdxzzz1dHlIeqGPXHI/H06Vcxw474J/92dDQQHJyMkOHDsVms5GYmMimTZvIy8sjOjq6J75VKkz87Gc/s3R9qwiL4C0iU4GhxpglxhiPiNiMMZ3zBG34Uyc/B74D/MAY0ye6cjab7WsjNsaNG0dxcTEtLS3U1dWRmZlJfHw8ubm5XYZIRUdH4/P5ePfdd/n2t7/N888/z8CBA3nyySd56KGH2FBVRVNTI4sX/4IZM2bwyCOP0NTUdNj2nHTSSYD/YWXnIOzz+YiNjSUyMpKEhIRgjjspKYmEhAQ2b94czM13zCBVJ6YLL7zQ0vWtIqTBW0RsQCzwmP+txBlj/taR1zbGeACMMa7AxJwrgLmBlEmfFRsby/e//322b99OZGQkWVlZRERE8Oabb5KUlERKSkqXsgDr1q1j1qxZvPPOO9jtdn4xfz5v/vUl3m2OY9q0aSxevJjIyMgj3js6OpqYmJiDLjzVsTfmunXrOOmkk9i/fz+DBg0CvlqzpeO/6sS1Y8cOAPIPGIJqlfpWEdLgHehdN4vI04AXmCIiMcaYhzoCdyevAUuNMdt7vaFhKCIigtGjR9PW1saIESOoqKggOjqaoUOHkpaWxv79+3G5XGzbto3Y2FjsdjslJSXMnj2b5++5h6HXXMNVHi9XIVS9+TrewEbHR7Jq1Sra2toOOiuzsbGRPXv2/P/2zjw+jurK99/TarWkVmtfLMvWZmMLycI2BrwIDLaBsAQ8eEIe2UgIkEC2N0wmk2R4YTJhGLK9rJPHmoQtPAKTkEkAO5NgGJsYG2MTYhvLi2xLlm1kyZskq9VqtXTmj6puWorkVd3VLd3v59MfdbWq+v50VXXq1rnnnsNLL71EMBhk6tSp+P1+8vLyBuVsMYxt7rjjDuDM46ydPj5ZSJQJyxBQBjwBzBWR79uLbxCRi0Vkuqrebwz3YFwuV2SicvPmd5DLIQAAIABJREFUzWzatCkyiSkig1wo4bJle/fu5dnPfx4dGCAt2EtKX5Csp58+5TZffPFFPB5PpKp9NLt37+Yvf/kLGzZsYMWKFXR1dZlFOmOQyspKRGTYV0VFhdPyxg2JYrx/C7Sq6kpgA1bsdnhodyFWUizDCSgtLSUYDBIMBvF4PHR3d/9VhEfYgC/v6iIYCKBAig6Q89hjpK9addI2VJUXX3yRxYsXjxiznZaWRkpKCuXl5Rw+fJjnn3+e5uZm9uzZw9q1a1m/fj2BQFyTURpGmebmZlR12FdTU5PT8gZRUVEx4o1GRKisrHRa4hmTEBOWWBOS1SLyKeBO4JtYI/CPAD8eq7m7z4ahi2dqamq47bbbWLduHSUlJTz11FMcOXIkEgESTUpeHn/o62Npfz+Cle+geOVKBj76UfLz80lJSRm2zW3bttHY2EhbW9ug6BewFg+53e6ILzwcDXP8+HFaWlro7u4mFApRVFQUWWhUWFgYacusxDTEgpPdTJL5vEsI462qB0SkBbgHK+XrCyKyGGg0hvvU8Hg8vO9976OqqootW7ZQX19PVlYWKSkpvPPOO7S2tkb2TU9P5ymvl6VdXYQ7d+BDHzppGytWrAAY1mUSxu/309/fz8DAABMmTKCiooKDBw+yfft2Lr30Ui6//HLy8/Mj2QonTJhw5n+0wTCOSQjjbfMo8FtV3WhvrxoSLmg4CeHRb3t7Oy6XVVl+5syZzJs3j5/+9Ke0tbVF9l2blcXfX3ITVzT+mR/v3Mjx++7jI7t3c9111zFS2twVK1Ywa9YsDh48OKIGVSU1NZWcnBy8Xi/btm2jo6ODjo4OSkpKuOiii+jv7yc9PZ3CwsJR7wOD83zta19L6uOThYQx3qraArSEE04Zw31mVFRUsHjxYsrKynj00UfxeDwUFBRwwQUXsGrVqkF+8A0VNbxZfi6bW3dQ3tnJF77wBb72ta/xjW98g1tvvTXi1+7t7eXb3/42a9euxev1kpV1wvTpqCoNDQ00NDSQnp7OggULuOyyy1i2bBldXV289tprLFy4cET3jCG5ueKKK5L6+GQhYYx3GOMmOTOCwSAtLS2UlZUxbdo0fv7zn7Nz5062b99OUVFRpIDxzp2Dq8eJCFlZWRw5coSioiJqamq46667eOaZZ/jRj35Eb28vd955Jw0NDRQXF5+Slp6enkgRhsLCQqqrq1myZAmHDx/mtddei8ThXnvttaPbCQ6RrJWaYsXbb78NwGw7j3yyHZ8sJJzxNpwZLS0tNDY2ApCZmcnChQvp6OjA4/HQ0NDArl27Iu6KvLw8hlY5FxFSU1MZGBhg4sSJNDY2Mn++VX+jtLSUyZMnU1lZGWnjRBQWFjIwMEBKSgrnnHMOs2fPpqCggLfffpuVK1cya9YsLrroIl599VW2bNnCvHnzOP/8809pkVAiYgoQD+auu+4CzjzO2unjkwVjvMcIZWVlkZ8ul4sZM2Ywe/Zs3G433/3ud3njjTfo7u4mLy+PSZMm0dHRMez3iAg5OTmRggsiQmZm5mm5ODo7O5k0aRIzZ86ktraW2tpacnNzI3Ho4SX7f/jDH2hsbKSxsZHCwkKmTp06Kn1hMIwHjPEeI3g8nkHGr7Kykp6eHjIyMrjnnnt48skn2bBhA+3t7UyePJlgMEhABFUGGebs7Gx6enoA/qo8WnFxMYcOHRq2/bS0tEG/a29v591336WoqIg33niD2tpaPvaxj7Fr1y6OHz/O+vXrmT9/Pueffz4LFiyI3HwMBsOpYYz3GCVceX5gYACXy0VpaSnFxcV0d3fT09Njx38LoVAwYnT7+vrw+/2EQkMzE1i0t7ePmJfb7XbT19dHTk4Oubm5kTSwbW1tvPXWW+Tk5DB//nx+8pOf8Pzzz5OXl8f06dOpqKgwqzANhjPAGO8xjt/vp7Gxkb6+PiZNmkRhYSE+n4+CggK2pnjp7rayAY7kRjldRISjR4+SlpZGMBgkPT2d6dOnMzAwgNvtxuv1csstt0QmNIcu9jEYDKeGMd5jHK/XS0ZGBiJCeXk5O3bsoKysjNraWna8vBewQgHBGq17PJ4RR96nQmdnJy6XCxEhFAqRlZXFsWPHAGhsbIwkp3K5XHi9XlMabQxy//33J/XxyYIx3mMcl8tFTk4OhYWFTJ48GRGhuroaESEt7SAgFBUVEQqFCIVC5ObmIiJ0d3efcXs+n4/i4mLS0tLo6uqipaUFv9/PlVdeid/vZ9OmTRQXF1NcXGxcJmOQ+vr6pD4+WTDGexxQUlKC2+0mFApRWlpKc3MzgUCAQG8Al8vKBFdWVsaaNWvweDwUFhbS29t7RiPwUCjEhAkTmDJlCnV1dRw+fJi1a9dSV1fHli1byM3NpaGhgX379nHddddFjhsYGMDv95vR+Bjg9ddfB87ciDp9fLJgjPc4wO12U1JSQigUYt++fezZs4etW7cSSqsjJcXNxIkTmT59Oj09PQSDQd59912Ki4vp6Og45RF4dIKf7du309nZSW9vL+3t7Rw5coS0tDS8Xi91dXW0t7fj9/tpaWlh2rRpwHu1OgEzGk9y7r77buDM46ydPj5ZMMZ7HOF2u5kzZw5ZWVlMmzaNxzd24vf3UKAF+Hw+rrnmGgoKCnj99ddpaGhg06ZNg0bB6enpEf/4UDweD16vl/T09Ehq2kAgQEFBARkZGeTl5dHS0sKkSZOYN28ejY2Nw1b8MROYBsOpYYz3OMPj8TBjxgxqampY3vI6oVAf51eeT11dHRkZGezevZuBgQHa2tro6OgYFBoYCASGLX8GkJGRQV9fHwUFBZFSbBkZGVRXV9Pb20tubi7p6emsWrWKSy+9lJkzZ+JyuQYVJfZ4POzZs4eysrKkXW1pMMQL41wcp7hcLtzuFPr6QgQCgUj0R2FhIVVVVZSUlJCRkUFGRkbkmOzsbGD4HMjBYJBQKMThw4dpb2+nqamJo0eP4na7cblc1NTU4Pf72b9/P62trfh8PhobG9myZQs7duwA3lviH10ByGAwDE9SjbxF5DqgHsgCvgF0qOpfVxswnDKpqW68Hi9btmxhwYIFVFRUkJ2dTWdnJwUFBbzxxhs0NjbS399Pbm4uoVAo4hYZitfrJScnh1AoRHl5eSSKpb+/n0AgQGpqKtOmTWPOnDnAe8WIS0tLWb58OT6fj8mTJ5vVlgbDKZA0xltELgAeAj4L3AD8BPiliLyqqidcYWKyvo1MWloavkwf69atIxAIcMMNN1BdXc3Ro0fJy8sjLS2Nnp4eXC4XaWlppKWl0dvbG/GF+3y+iGulpKQEj8fDhAkT6O3tpbW1lbq6Omprazn33HNZu3Yty5Yti/i1w0WJX3nlFZYvX052djY33ngjR44cobCwELc7aU5PQxQ//OEPk/r4ZCGZro7pwB9U9XfA70TkDuD9QL+IrAD6R0ona7K+jYyIsGjRIg4dOsS8efPo6OggNzeXKVOmkJ6eTigUYv/+/UyaNImWlhYGBgbIzc2lv78fsHJ35+bmUlRURDAYpKqqitbW1si+PT09+Hw+Ghoa6O7uprm5mWnTpg2aCJ07dy6BQID8/PxBVXZKSkoc6RPD2XG2qVidPj5ZSCbj/QbwCRGpV9XXVfVhEXEBNwOrTzb6NozMgQMHyM7O5sCBA4RCITZv3kxdXR0ul4uqqioWLVrElClT2LlzJ+np6YC1knLPnj309fVRXFwcyVXS0dFBb28vnZ2dHDp0iO7ubjo7OznvvPOorKzE6/Wyc+dOpk6dGhlZ+3y+SG7vUChkquwkOS+//DJw5kURnD4+WUho4y0is4FeQFR1q4hsBBaKyBFV3aaqD4rIxcA/AuOj9lEMCPueKysree211zh69Ci7du2ip6eH5uZm+vv7OXDgACkpKUycOJGdO3fS2trK8ePHmTBhAuecc05kKbyq4nK56OrqorW1lTfffBO/38/ChQvxer1s3bqVrq4usrKy8Hq9+P1+iouLI4Y8HJNucJbKykqam5uH/V1FRcUJj73vvvuAMzeeTh+fLCSs8RaRa7BcHb8FFonI14GfAV8B/kZEilV1NbAeMKs6zoKw7xng8ssvZ8eOHUyZMoVt27axd+9eUlJSKCkpobKykurqaiorK3nnnXc4ePAgZWVlHDlyJFKrMjMzk5ycHM4555xIqOFLL72Ex+NhyZIlZGdnM2XKFLxeL01NTRw/fhyXy2UMdoLR3NyMKWqV2CSc8RYrDi0T+AJWJfnficgC4BfA3cD9wCeBb4jIQazok/c7pXesEW3IZ8+eTWFhIUeOHKG2tpa2tjaampqoqakhNzc3EkHi8XjIzMyMLIGfMWMGXV1dzJkzh5/97Gfs2rWL3/zmN+Tk5HDllVeSl5fH8ePHKS0tpaOjg/z8fIf/asN4paKiYtjQ1/Dvmpqa4ivoNEg4421POh4XkQ1AtoikqupaEfkw8BzwRVX9FxGZDJwPfFlV9zqpeazidruprKykvLwcv99PV1cXx44dY86cORQVFbFmzRo2b96Mz+ejoqKCgwcP0tfXx65duygoKCA1NZVly5bR1NREZ2cnq1evpr6+nnXr1tHS0kJFRQW9vb34/f5IiKLJa2KIJycyziMZ9UQhka+UVuByIANAVddjTU7eLSJTVXWfqr5gDHfsCecdycrKorq6mkmTJpGTk0NRURGVlZXU1NQgIvT19XH06FHKysqYMWMGaWlpFBQUcO+99zJ9+nSuuuoqzjvvPLq6umhubuZPf/oTa9asYeXKlezevXtQZfswoVCI1tbWs0pTG2tE5NMiskFENgytDWowxIqEG3mLSIqq9qvqAyLyLPCQiHwG8KvqayKyCeh3WOa4IjouOz8/n2AwSEdHB7NmzeLSSy8lPz+fPXv28PLLL9Pb20tOTg4HDhxg06ZN9PT0cMkll3DzzTeTmppKW1sbVVVVeL1eSktLOXLkCFlZWRQXF+P1egmFQhw6dCgS533o0KGEDx00oaiDefjhh5P6+GQhIYy3iFwCVKnqU6raLyIeVQ2q6k0i8kvgh8A6EXEDlwGJOwwbg4RzdKsqmZmZdHZ2EggEKC0tZeLEifT09LB//36mT59OTk4OVVVVbNy4kfLycjo6OkhNTSUvL4/Ozk5ycnLw+Xyce+65BAIBampqBhVmaGtrG2SswyGDJnQweaiurk7q45MFR423HaftBR62NiVTVR9S1aCIpKtqQFU/JCK3AqXALGCpqu5zUvd4pq+vj66uLiZPnozX60VVefPNN3nmmWdobW1l6dKl5ObmRhbrzJs3j8LCQvbu3Ut9fT1lZWUEAgHcbjfZ2dl4vV5EJOJfHGqsTehg8vHCCy8AcP311yfl8cmCo8ZbVQewJiefwHKF1ItIhqr+QFUDUfv9HEBE0lR1+JykhpgjIrS0tNDU1MTEiRPxeDy43W7S0tLYt28fe/fuZf369VxzzTX09/ezfft2jh07RmZmJg0NDZGIlGAwiNfrHXb5uzHWyc/3vvc94MyNp9PHJwuJMmEZAsqAJ4C5IvJ9EfkmgIjUi8gce7+/zoZkiBsDAwPk5eUxZcoUpk+fHhk5z549mzvuuIMlS5Zw9dVX43a7WbRoEUVFRQQCgchCn4GBAY4fPx75megTkQZDIpMQPm+shTgfVNWV9qrK+4DH7d9dBDwLkTBCg0P4/X4CgQATJ06MLJMHayLzhhtuoL6+PpLvu6Ojg8WLF7N//358Ph/9/f3k5+fj9/vx+XwcPHgwUrHejLQNhtMnUYx3D1AtIp8C7gS+iTUC/wjwY2O0E4Phqt1E154sLS0dFC0CVkX61tZWZs2aRXp6Orm5uZGCxG632yzQMRjOkIRwm6jqAaAFuAdrEc69wPeA14zhThzCUSfRC2nCMeDhGO1wtZ1QKERVVRWpqals3bqV3t5eSkpKOHbsGO+++y7t7e2ICG+99RaBQGCkJg0Gwwgkysgb4FHgt6q60d5eZU9oGhKYoaPx6J8ul4tDhw5x6NAhXn/9daqqqigsLGRgYID09HS2bdvG7t278Xq9kSX5huTnqaeeSurjk4WEMd6q2gK0iIiohTHcSUB4ND7Sdn19PQcPHsTr9bJhwwZqa2spKSnB5XIxZ84cfD5fJKuhYWxwtpWQnD4+WUgIt0k0xk2S/Pj9ftatW8fx48cB+OAHP8ikSZMIBAI0NTVFXCzhJFjRk5+G+FFZWRmJsR/6Olna1xPx7LPP8uyzzybt8clCwoy8DWOHTZs2sXHjRvx+P+eeey7Z2dlceOGFtLW14fV6B014GpwjVmlfH3zwQQBuuummpDw+WTDG2zDqhP3XdXV1wHv+79LSUidlGQynRaKnizXG2zDqeL1e5s+fPyrfNTRRlcEQLxI9XWzC+bwNhmjCWQUPHTrktBSDIaEwxtsQF4LBILt27SIYDJ5Wju7CwkJKS0tNVkGDYQjmOdQQF1paWmhsbAQgMzPzlHN0m0RVZ8fZFBI+U371q18l9fHJgjHehrhQVlYWSWwVjgM3o+nY40Qh4bP9vzp9/KmQCJOZxm1iiAsej4eJEycSCAQIBoP4fD5Wr14diQU3jB0ef/xxHn/88aQ9/lRoampCVYd9jfSkM9oY422IG16vN5JGdv369axdu5b169c7LSvhOdFimpO9YuUaORFOG994GO8TER6VD/eqrKwctXaM28QQN6KXzs+dO3fQz2RGRD4NfBogJydn1MPIKioq4u76MJw5J3KZhG/Eo4Ex3gZH8Pl8LFmyxGkZo8LQAsQbNmxwWJEhUTmZL/x0DHvSuk0kEaLkDQaDwSGSzniLyAwwCawMBsP4JqncJiJyFfBtEfmgqu50Wo/BYPhrli9fntTHJwtJY7xFZCnwVeBzqroznPf7FI+NTCiVl5fHUKXBYDjbrJFOH58sJLzbRERcIuLCKkrcr6prRGQC8A8i8i0ROV9Esk70Har6iKpeqKoXFhUVxUW3wTBeeeCBB3jggQeS9vhkIeGNN1BsV9W5GMgTkV8Bv8B6asgD/g4wNbQMhgThueee47nnnkva45OFhDbeInI10CQi16hqFzAPmAW8rqrfUtU7gHeBjzup02AwGOJNwvq8bcN9H/Ab4DwReVVVu0WkBuiP8nk3ArUikqKq/U5qNhgMhniRkMZbRC4DvgXcCfQBDwA/BwKqGora73bgM8DHjeE2GAzjiUR1m1QCn1bVdaq6EdgKfF9E3AAikmKPwJcCt6jqO85JNRgMhvgjibTWZWj4n4i4VTUkInOAzwHfUNW99u88gEdVTystnYi0A/FJ+/UehUCiloJJVG2JqqtaVUeMbjrB+ZUIf4/R4Hz7J9NQoaqnFBKXUMY7zDBGPAPL9/22qn7VOWVnhohsUNULndYxHImqbazpSoS/x2hwvv3R1JAQbhMRmSEil4lIMQxe+i4iLlXtAb4AXCEiyZ+GzmAwGM4SxycsReQa4NvAbiBVRD6tqvvDv1fVAXuRzhFgBbDHGaUGg8GQODg68haRRcCPgNtV9QYgCNTYv0sJ76eqA6p6GLhfVdud0HqWPOK0gBOQqNrGmq5E+HuMBufbh1HS4KjP244YKVHVV0WkBHgLWA8cBNaq6uP2ZCWq+tbp5DMxGAyGsYyjI29VbVDVV+3N24AH7BH4OuBaEakELgMO2Psbw20wGAwkaLQJgIisAP5OVXc4rcVgMBgSjUSJNpEh2x8AioEuZxQZDLHFVIJKjD5IZg0JYbzD7hARSROR24B7gU+o6rvOKos90ROziUIiakoGTuUiNJWgEqMPRCQ1ATScVT8khPGOYgArS+DfquoWp8XEEhGZCaCq/YlkLEXkUuATIpLmtJZoRKRSRHwiku20lqGISB4MGoQMa8TtSlBPi8i0OMobquE6EblfRP5dRArDRiyO7SdCHywFfiwiT9hrTPId0HDW/ZBQxltV+1R1uapud1pLLBGRqcBbIvIcJI4BF5FrgaexlncHoz539NHSzjD5PPBD4IsikmHH/juOiFyHdRH+TEQ+JSLlqqpD9dkG4+tEVYJyQOsFwENYAQGZwE+A94tITpzaT4Q+qAMeBP4DK6rtc8DHRWRSHDWMSj8k7ITlWMYOi/wJVm7y3ap6VdTvXHbxiXjqEcALPAw8raor7BGuYuWPOexUyl0RmQX8CrgdyABuUNU7w/3kRH9FaZsOvAJ8BCuZWhVQB3xVVXfZ/Rq+MN8GOlR1oViVoG7GynHxLNBo56uPtd4PA1eq6q329h3AhcDvsBbA9cfKjWD3xZ+BLof74HKsdSUftrevAxYBe4EnVfVYDNsO39BH5VxIiNHLeMI2Nq3Am8D5gIrIr0VklohMdcIQqUU30AZ0ichk4PfAD4DVInKB/XQQt5FSVFslwCpVXYWVXXKRiPw78F2n+iuKfuCPqrpaVZ/ESlv8Z+B+EZlsG8JEqgT1BlAqIvUAqvow1tqKm4HMWPp/7e++BMh3og9EJN1+uwEoE5Fltq4XgdVYN93CWGpglM8FY7zjRNgYRRmbImCZql4NTMO66CvsfePuQhEr3W4XcC2wDHhKVW/HekL4tYiUxHlyJ3rEerOIPAJsAh7DGikeBf6PiGTGUZMl7L0RVBtwgYj8bwBVbcG6GLdj3WQcrwQlIrNFpEZEalV1N7ARWCgi59qaH8Rykf1jjNqvE5FqEZlpZwCdB8wmvn2wBLhNRDJUtQN4AqgXq24Aqvo7IADcFUMNo34uOJ7bZDxgTwLWi8hWYJOqNgH/BUwVkVIs/+NfgH8CXomXe0JE5mPdMPaq6loR+THw31g39ZvAurhtX2ncbigishi4QUQ2q+pPReQcoBZwqeq37X0OYyUr642XrhG03Y410j6qqk+p6l4R2Y/VfxNxsBKUWHmDHgF+i3Uz+TrwM+ArwN+ISLGqrsZa1ewbzbbt9q/Fylu0FqsP7lfVF+w+CMWpD64GvoPlX+6xP34FmAAsFZEJqvoc1vVXEws3nMSqKpiqmlcMX8ASrNy9XwJeAr4P3ADkYo0qjwNX2fu+BJTHSdeVWBOTXwb+aJ9ck4ACoAHrhC/FeqTeDEyIk6732SfxF+12b7I/T7d1ftLe/jCwCsiL4/9yqLYb7c+vB14G/sHevh84hrU6+AIsd0XhMN93O9ZIeMYo6xQsY7wcWGp/tgDYhXVTqQD+BXgV+CWWv/e8UdZwIbANa4QpWCPKH2ENGF2x7gP7u2diPaGF/0+FWE+8hUAacCuWG+XXQAswMwYaLrOv8/mjfS7E5aQfzy+sZf+ft99X2Mbwp/ZFdBOwMM56xDaEj2FN/oH1+NZg66q0T/BHgG/aBnLUL6wRdGUAjwMfsD/7JPBvWDeaImAxsAN4zv4Zc12noO0y4Bxghn3xPQm0Ax+JOv4x+3O3vZ2ClYDtd6NtNIfovhf4GJBqb88FmrBCcQEmY914Rn3AAFwN3By1fYl9s0ixt912H7wYqz6w/96Hgc9j3Ux+DzwF7AQ+ZO+ThXWDKY2Rhk8Ac2NxLsT8xB/vL+AWrEmhPHu7yDbg3wTK7M8EO/Injrrus1/p9vaj9on+FXs7zf6ZE2dd99gX2SVAh63pv+z+moF1A7wy3HcJou27WE8pmVgTrMX2/uELdA6Wu6I86rs8gC/Gej9rG4vsqM8WYo02p8ahvyZGvfcBy6O2i8Kfx1jDxVghpsfs/nDbfbAPmBXHc2fUzwUzYRljVPVxYCVwt4jkqJXSdg3W3b7W3kfV/i/GmqgojjVYF9TjIvIAkI812v6EiFSratiX3BlnXc9gPWZ+BnhErYmcz2IZ7fNVtVlV/6jW5GBcOIm2z2C5my5X1W61IonaAfS9YtkN9j6fDX+nqgb1NEv4nYbeFLuNB7BCQB8SkRwRSVXV17AmfmM2r2JPfqODV0i7gcli1Z+9Beu888awD8Ia1mCF4N1m90e/3Qe/B7pj0XaUhui01qN+LhjjHR/+w/55j4jkqzXr/ybW43bMGSHE7w9Yj28vY4XgfVitYs+vAH3hnWJ5UxlOl6o2qlXq7tdApn2B78IqwjF1pOMc1LYbK0pnsThYCUpELhGRm+32+8Wq8Yqq3oR1nf8QuFVEPofl6gmN+GVn335I/noRVQ9WwZWvAHcCX1ZVf4w1hA34WuAF+72KyE1YYbo9I37Z6Gjoj+6H0T4XTLRJfHgTa8HL/wJWich/YvlML45T+ynYF6t98rrtkcBb9gsAO3LiMiwXhRO6omfZ/4zlnvhXEWkCPghcF943gbRlYLnGVgOPSZwrQdnfHV5gJSKSqaoPqWpQRNJVNaCqHxKRW7FcO7OwJjH3xbj9gajzDFXtFWtBysew5g0aRqP9k2gI2U8bfXZ/uLEWVH0Ja7Cy/0TfO0oaIv0w2ueCWWE5SojIhcBRe5R4ov1uwjLkf9E4pAGww8Vuw5pM26eqT9mfp2DZpQH7/Vzg/wPXaxzyypxAlytK03W2rknAd1R1a6x1naa2G4GlwKdV9VkR+Q3w/1T15eHCvcSKMx710Z793V/GcoXMAv6sqj8YYb+0KJdYXNsXkbuB51V122i3fxoalgGbVbXRKQ32fmd9LhjjPQqIyJVYE1e/wZrwi8mJcbrYj2VPY4WFDQBfBV5S1bvt3w+KaRWRIo1DmblT0JWqqn1R+3tUNTjcdzmpzY7TnayqfxSHK0GJyBeBciz3wO1YCz56VfWfxFpRGYilhpO0fzHQrjHOzX8KfXB0NEf8Z6Bhgf1+VP4Pxud9ltiPzXOwHp33AP8s1qKSoftdLiLxckeE8QD/rapPq+ozwFXAh0Tk3yDySH+ZfcKBFY+eCLr6RGSRiPy9vX/fSF/kpDbbELxsH+d0JajfAq2quhIrmuQzWGsJAC6Kg4YTtX8h8cnNf7I+OOqwhrmM5v9B4xQqM5ZfWAmJwk8xD2Itv50+ZJ9soCLOui7Aih/Nj/osPEIML36ZgTV6HPe6RlMblk9zeqx0DtNeKVZY4Kew4pj/GSuG+iPhc3OQ96niAAAC3ElEQVQstz8eNcTlxBpvL6xJiyew4n4/ipXDxCktDwDrh3x2G/AFh/soIXWdibahFyXwASx/+cRY6DuB7nuxVkteb28vJo7x8E63P940GJ/3WRI9MSUieap61H7/TeD9QA5wjcZpsi1KV8RPLCIvYj263aiqrSLyVawY81uww8zHu67R0CZWAYuPYS2fv0njXFBERMqwFghttLfjmi7X6fbHmwZjvM+C6H+KiPwr1gKO/1QrvvNW4HvAJar6joO6/gnLB/u3WFEbQayZ8BuNrtHVJlZVmiuBXepgQZFYTowmQ/vjRYMx3mfIkIv9O1iTEVeoFVuaizVR8ZKqbnJY18WqerG9XYMVi3pYrcyG415XomszGEbCGO8zYMjF/n+xJrCutw13OBZ4ULib07riqSNZdNl6ElabwXAijPE+C0Tke1gZwZbahtuRUmFG19mTyNoMhuEwcd5niIiUA9Uk2MVudJ0+iazNYBgJM/I+C8ITEol2sRtdp08iazMYhsMYb4PBYEhCjNvEYDAYkhBjvA0GgyEJMcbbYDAYkhBjvA0GgyEJMcbbYDAYkhBjvA0GgyEJMcbbYDAYkhBjvA0GgyEJMcbbYDAYkhBjvA0GQ0wRkUoRiWthivGAMd4Gg8GQhBjjPc6xR0XbROSnIrJFRJ4WkStEZI2I7BSRuU5rNIwJ3CLyhIhsEpFfiYjXaUHJjjHeBoBzgB8BM4FzsSpdXwJ8CbjbQV2GsUM18IiqzgQ6gc86rCfpMcbbALBHVTfbFWXeAVbatfc2A5WOKjOMFVpUdY39/hdYgwPDWWCMtwGgN+r9QNT2AOCOvxzDGGRo7mmTi/osMcbbYDDEg3IRWWC//zDwJyfFjAWM8TYYDPGgAfiEiGwC8oEHHdaT9JhKOgaDwZCEmJG3wWAwJCHGeBsMBkMSYoy3wWAwJCHGeBsMBkMSYoy3wWAwJCHGeBsMBkMSYoy3wWAwJCH/A+ON5XKbgA3EAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "samples = ap.sampler.get_chain(discard=ap.iburns[-1], flat=True, thin=ap.ithins[-1])\n", "\n", "fig = corner.corner(samples, quantiles=[0.16, 0.5, 0.84], truths=[mTrue, bTrue],\n", " labels=[\"m\", \"b\"], show_titles=True, scale_hist=True,\n", " plot_contours=True, range=([-1.1,-0.8],[3.7, 4.8]));\n", "\n", "# Plot where forward model was evaluated\n", "fig.axes[2].scatter(ap.theta[:,0], ap.theta[:,1], s=8, color=\"red\", zorder=20);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The predictions are close to the true MCMC!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }