{ "cells": [ { "cell_type": "markdown", "id": "21394ae3", "metadata": {}, "source": [ "# Nested sampling and Bayesian evidence\n", "\n", "Nested sampling computes the **Bayesian evidence** (marginal likelihood)\n", "\n", "$$\\mathcal{Z} = \\int \\mathcal{L}(\\theta)\\, \\pi(\\theta)\\, d\\theta$$\n", "\n", "and produces posterior samples as a by-product. The evidence is what MCMC does not give\n", "you: it is the quantity needed to compare models.\n", "\n", "This notebook uses `psiphy.mcmc.NestedSampler` with the\n", "[nautilus](https://nautilus-sampler.readthedocs.io) backend, which needs few likelihood\n", "calls (useful for expensive likelihoods) and can evaluate likelihoods in batches.\n", "The [dynesty](https://dynesty.readthedocs.io) backend is a drop-in alternative if it is\n", "installed; the interface and the result format are identical:\n", "\n", "```python\n", "NestedSampler(log_likelihood, prior, backend=\"dynesty\") # instead of \"nautilus\"\n", "```\n", "\n", "Contents:\n", "1. Checking the evidence against an analytic value\n", "2. Posterior samples and getdist export\n", "3. Model comparison with Bayes factors" ] }, { "cell_type": "code", "execution_count": 1, "id": "6a3ae90c", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T11:44:47.395787Z", "iopub.status.busy": "2026-10-14T11:44:47.395686Z", "iopub.status.idle": "2026-10-14T11:44:49.035442Z", "shell.execute_reply": "2026-10-14T11:44:49.034866Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from getdist import plots\n", "\n", "from psiphy.mcmc import NestedSampler\n", "from psiphy.utils.prior import UniformPrior\n", "from psiphy.toy_models import NoisyLine\n", "\n", "np.random.seed(1)" ] }, { "cell_type": "markdown", "id": "c67cbd09", "metadata": {}, "source": [ "---\n", "## 1 Evidence of a Gaussian likelihood in a box\n", "\n", "Take a normalised 2-D Gaussian likelihood and a uniform prior on a box much wider than the\n", "Gaussian. The integral of the likelihood is 1, and the prior density is $1/V$, so\n", "$\\mathcal{Z} = 1/V$ exactly." ] }, { "cell_type": "code", "execution_count": 2, "id": "bf14c550", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T11:44:49.037020Z", "iopub.status.busy": "2026-10-14T11:44:49.036848Z", "iopub.status.idle": "2026-10-14T11:44:53.343033Z", "shell.execute_reply": "2026-10-14T11:44:53.342659Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "log Z (nautilus) = -5.988\n", "log Z (analytic) = -5.991\n", "effective sample size = 10066\n" ] } ], "source": [ "mean = np.array([0.5, -1.0])\n", "cov = np.array([[1.0, 0.3], [0.3, 0.5]])\n", "icov = np.linalg.inv(cov)\n", "norm = -0.5 * np.log(np.linalg.det(2 * np.pi * cov))\n", "\n", "def log_like_gauss(theta):\n", " d = theta - mean\n", " return norm - 0.5 * d @ icov @ d\n", "\n", "box = UniformPrior({'x': (-10, 10), 'y': (-10, 10)})\n", "ns = NestedSampler(log_like_gauss, box, backend='nautilus', n_live=500, seed=0)\n", "res = ns.run(dlogz=0.01).get_result()\n", "\n", "print(f\"log Z (nautilus) = {res['log_evidence']:.3f}\")\n", "print(f\"log Z (analytic) = {np.log(1 / 400):.3f}\")\n", "print(f\"effective sample size = {res['ess']:.0f}\")" ] }, { "cell_type": "markdown", "id": "b3d3dc69", "metadata": {}, "source": [ "`get_result()` returns the same keys for every backend:\n", "\n", "| Key | Content |\n", "|---|---|\n", "| `log_evidence`, `log_evidence_err` | $\\log\\mathcal{Z}$ and its error (`None` for nautilus) |\n", "| `samples`, `weights` | weighted posterior samples (weights aligned with samples) |\n", "| `equal_weight_samples` | resampled, unweighted draws |\n", "| `log_likelihood`, `ess`, `backend` | per-sample log-likelihood, effective sample size, backend name |" ] }, { "cell_type": "markdown", "id": "a09bb1de", "metadata": {}, "source": [ "---\n", "## 2 Posterior: fitting a noisy line\n", "\n", "The likelihood function takes only `theta`; the data are captured in a closure. The prior\n", "object supplies the unit-cube transform and the parameter names." ] }, { "cell_type": "code", "execution_count": 3, "id": "b8a9d727", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T11:44:53.346276Z", "iopub.status.busy": "2026-10-14T11:44:53.346137Z", "iopub.status.idle": "2026-10-14T11:44:58.995996Z", "shell.execute_reply": "2026-10-14T11:44:58.995631Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "posterior mean: [-0.793 4.893] truth: [-0.8 5. ]\n", "log Z = -53.90\n" ] } ], "source": [ "line = NoisyLine(Nx=40, noise_sigma=0.8)\n", "theta_true = np.array([-0.8, 5.0])\n", "y_obs = line.simulate(theta_true)\n", "\n", "prior_line = UniformPrior({'m': (-3, 3), 'b': (0, 10)})\n", "ns_line = NestedSampler(lambda t: line.log_prob(t, y_obs), prior_line, backend='nautilus', seed=0)\n", "res_line = ns_line.run().get_result()\n", "\n", "post_mean = np.average(res_line['samples'], weights=res_line['weights'], axis=0)\n", "print('posterior mean:', post_mean.round(3), ' truth:', theta_true)\n", "print(f\"log Z = {res_line['log_evidence']:.2f}\")" ] }, { "cell_type": "code", "execution_count": 4, "id": "6b5f7a1e", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T11:44:58.997108Z", "iopub.status.busy": "2026-10-14T11:44:58.997033Z", "iopub.status.idle": "2026-10-14T11:44:59.211162Z", "shell.execute_reply": "2026-10-14T11:44:59.210773Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Removed no burn in\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAHqCAYAAAAgWrY5AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWF9JREFUeJzt3QeYVNXdP/Dv1N2d7cuy1GUpywJLL1IUkCKIUkRArG80xWgs+FcTTTQxmhhLfI0xQaOJJa+JQUEQFAsoRYr0XqT3pbftu9Pu/zlnHUJZYGd3Zu65934/zzMPsMzu3rkzd75z2u/YNE3TQERERMqx630AREREVD2GNBERkaIY0kRERIpiSBMRESmKIU1ERKQohjQREZGiGNJERESKYkgTEREpygmLCAaDOHjwIJKTk2Gz2fQ+HCLDEvWPiouL0bhxY9jt/JxPFE2WCWkR0NnZ2XofBpFp7N+/H02bNtX7MIhMzTIhLVrQoTeWlJQUvQ+HAHzxxRe47rrreC4MpqioSH7gDV1TRBQ9lgnpUBe3CGiGtDpDEHwujIvDRkTRxwEl0o0Y0yQiootjSJNuOnfuzLNPRHQJDGnSdUyaiIgujiFNRESkKIY06aZnz548+0REl8CQJt2UlZXx7BMRXQJDmnSzceNGnn0ioktgSBMRESmKIU26GTNmDM8+EdElMKRJN3PmzOHZJyK6BIY06aawsJBnn4joEhjSpJusrCyefSKiS2BIk264TpqI6NIY0qSbmTNn8uwTEV0CQ5qIiEhRDGnSTbdu3RAIath0yI/ZW3w4WRrks0FEdBbn2f8giqXV+7wYPbUQQQ1ommbHmoIAbuvmxuOD45Gd7uCTQUSWx5CmmNM0DS/OqcD2b9fgTze1wbB2LthsNhwpDuLVbyrQ589F+Pq+FLRtwKAmImtjdzfF3JvfVuL9VV70benEdfluGdBCg2Q7nhvhwQsjPbjm9SLsPxXgs0NElsaWNMXU5sMBPPVFORZNSEHTpNHV3ueOHnE4VBTEmHdKsHBCCuJdVSFORGQ1bElTzFT4NNz6XgleHOlBXpYDixYtuuh9fz4wHi3rOfD4p9zOkoisiyFNMfPaogrkpNtxV0+3/PeJEycuel/RBf7GeA8+Xu/DnG0+PktEZEkMaYqJogoNf5xTIVvRoTHojIyMS35PuseOv9/swc+mlMIf0PhMEZHlMKQpJl6ZX4Hr811o1/C/M7b79+9/2e8b1s6NZul2/HN5ZZSPkIhIPQxpiroyrya7up8YknDO16dPn16j739uuAe/m1Uhx7SJiKyEIU1R9+6ySvRr6UTr+rVb99wzx4kezRz422K2ponIWhjSFFXBoIZXvqmQs7XP16lTpxr/nN9fl4AX55SjuIKtaSKyDoY0RdVXW31IS7ChTwvXBf/ncl34tYtp38iJoW1cmLiwIsJHSESkLoY0RdXfl1Tinivjqv2/VatWhfWzHhscj4mLKuD1szVNRNbAkKaoOVQYxLwdftzStfqQDleHRk50bOTEB6u9Efl5RESqY0hT1Ly7vBLju7iRHF99Wc/hw4eH/TMfvjoOf2GXNxFZBEOaorbT1VtLK/HTPhdvRa9cuTLsnzukjQtHizVsPOSv4xESEamPIU1RsWJfAAkuG7plX3wPlyNHjoT9c+12G/7nCjf+bzm7vInI/BjSFBVT13kxptOlZ2+npKTU6mf/oEcc3l9VyVKhRGR6DGmKSlf3tPVejO1UtZHGxQwdOrRWP79NA4fcqEMs7yIiMjOGNEXchkMBiP0wOje5dIWxjz76qNa/466ecfgnu7yJyOQY0hRx09ZVtaJDu11Fw01d3Ji1xYfSSq6ZJiLzYkhTxE1b78PYzpfu6hbat29f69+RkWhH7+ZOfPEdu7yJyLwY0hRR248FcLIsiJ7NLr+ZRlJSUp1b01PWcpY3EZkXQ5oi3tV9Yye3XCp1OcuWLavT77qhg0t2eYutMImIzIghTRE1Y6MPoztevqs7EjKT7OiZ45BBTURkRgxpipjTZUFsOhxA35YXL2BytmuvvbbOv3NcZ3Z5E5F5MaQpYuZs98uAjnPWbFb3xo0b6/w7Rav9yy0++MSaLyIik2FIU8SIbudr29R8j+iCgoI6/86sZDvy6tuxaBdreROR+TCkKWJVxmRIt615SCcmJkbkdw/Pd+OzzRyXJiLzYUhTRGw9GoTo5M7LqvlLasSIERH53SPau/DZJi7FIiLzYUhTRMze4sPQtq6wqox9+OGHEfndXZo4UFSpYefxQER+HhGRKhjSFBGzt4bX1R1J4oPB9e3c+GwTu7yJyFwY0lRnYmb1wl1+DMyt2dKrkDZt2kTs7A/Pd2HmZnZ5E5G5MKSpzlbu86NVpl3W0w5HZmZmxM7+NW1cWLrHjxJuuEFEJsKQpjqbu92PQa3D7+pevHhxxM5+UpwNfZo78TX3mCYiE2FIU53N2+HDoNbhdXVHg1iK9Tl3xSIiE2FIU51U+DQs2+tHv5bht6SvueaaiJ59MXFNrNUWa7aJiMyAIU11IsaBOzZyIjm+5kuvQrZt2xbRsx9ao73taDCiP5eISC8MaaqTudt9GFjLru59+/ZFfClWqDVNRGQGDGnSZdKYEBcXF/GzP7SNC7M4eYyITIIhTbUmljutO+jHlc1r15IeM2ZMxM/+4DwnFu7krlhEZA4Maaq1Rbt86JHtRII7/PHoSJYFPVu6x47W9R1y7TYRkdExpKnW5u/wY0Bu7UuBBoPRmeAljkkcGxGR0TGkqdZEKdB+LWu/PrpVq1ZROfsDcp2Yv4OTx4jI+BjSVCvlXg3rCvzoXcvxaKFJkyZROfvig8OSPX5ZU5yIyMgY0lQry/eJ9dEOeGo5Hi0sWLAgKmc/jePSRGQSDGmqfVd3K322pqwJjksTkRkwpKlWxDKnuoxHCwMGDIja2ee4NBGZAUOawuYPaFi6N4C+dQzpSFccOxvHpYnIDBjSFLa1BQHkpNvlmuS62LVrV9TOPselicgMGNJUy/Houm9N6XQ6oz4uPY/rpYnIwBjSFLaFu+o+Hi3cdNNNUT37V+eKEqEsakJExsWQprCIvZoXySImdZ/ZPXXq1KiefVFTXKyXDgS5XpqIjIkhTWERezUnxdnQJK3uLx2v1xvVs5+ZZEfjVBs2HgpE9fcQEUULQ5rC8u2e2u96db6cnJyon/2+LV2y5U9EZEQMaQqL6D7uE6GQjlbt7rP1beGUE92IiIyIIU26hfTcuXOjfvbFWm4x0U2MpRMRGQ1DmmqssDyIPScD6NTYYZiz1qKeHSKf956MzraYRETRxJCmGlu+L4DuTZ1wOmq/qcbZ+vbtG/Wzb7PZqsald7PLm4iMhyFNNbZkd+S6uoWjR4/G5OyLcWlOHiMiI2JIU3jj0S0iF9Lbtm2LydkX49IMaSIyIoY01UgwqGHZPj9650S3lGc0iDH0facCOFnKcWkiMhaGNNXI1qNBZHhsyEqO3Evm1ltvjcnZF2PovZs75RpvIiIjYUhTzJdehcyYMSNmZ79vCxY1ISLjYUiTbiFdVlYWs7PPcWkiMiKGNOkW0k2bNo3Z2e+V48TaAj8qfCxqQkTGwZCmGhUxEROvOjaKbBGT/Pz8mJ39xDgb8hs6sHI/x6WJyDgY0lSzIibZkStiEjJ79uyYnn2ulyYio2FIU8yLmOiFO2IRkdEwpEmX8WihT58+MT37V7WsWoYl1nwTERkBQ5pqVsQkCiFdWFgY07PfINmOzEQbNh8JxPT3EhHVFkOaLlvEpJ7HhvpJkX+pbN68OeZnX+4vvZOTx4jIGBjSpEtXt57j0ou5IxYRGQRDmnQL6ZtuuinmZ18UNWFIE5FRMKRJt5D+8ssvY372W9e3o9Sr4cBpbrZBROpjSNNli5h0iHARk5Di4uKYn32bzYarWjixeJcv5r+biChcDGm6qGV7A+jRLPJFTEIaNmyoy9mXdbw5Lk1EBsCQJt0mjXXr1k2Xsy92xOK4NBEZAUOadAvpzz//XJez37WpA9uPBVBUwaImRKQ2hjRdvIjJXj9655hn+VWI22lD96ZOLN3D9dJEpDaGNFVry9Eg6ifZkBmFIiYhV1xxhW5nv2opFiePEZHaGNKkWxGTiooK3c6+nDy2iy1pIlIbQ5qq9e1uH66Mckhv2LBBt7MvPoCs2OeHL8BxaSJSF0OaqvXtbj+ubGG+8eiQ1AQ7WtRzYF0BN9sgInUxpOkCJ0uDOFSkIb9hdIqYhNx44426nn12eROR6hjSdIGle/3omeOAwx6dIiYh8+fP1/Xsix2xFnHyGBEpjCFN1Xd1x2Dnq1OnTul69kObbWgax6WJSE0Maap2ZveVLVxRPzOZmZm6nv3sdAdcdht2neBmG0SkJoY0ncMf0OSs51450R2PFvr06aP72ee4NBGpjCFN59hwKICcDIec/Rxtn376qe5nnyFNRCpjSJMu49GqYOUxIlIZQ5piXmkspGvXrrqf/fYNHThYqOF4CceliUg9DGk6x7dy0ph1WtJimZn4UCIeNxGRahjSdMbBwiCKKzS0rh+bl8WaNWvU6fJmHW8iUhBDms4Qa4ZFYNls0S1ioho5eWw3W9JEpB6GNJ2xaJdPBlasjBw5Uomzf0UzJ9YV+FHhY1ETIlILQ5rOWLjLj74xKGISsmTJEiXOvsdtQ4dGDqzcz9Y0EamFIU1SUYWG7ccC6No0+kVMQo4fP67M2e/b0oWFOxnSRKQWhjRJS/f40SPbCbczduPR6enpypz9q1s5MX+HT+/DICI6B0OadBmPFgYMGKDM2e/X0inXiPsCHJcmInUwpEkSs5v7tYzdeLTw8ccfK3P20zx25GU5sHIfu7yJSB0MaZKtRxFOvS1UDrQ6A3JdmL+DIU1E6mBIE1bvDyC3vgMp8bFdH92xY0elzv6AXCfmcVyaiBTCkCYs2u1DXx1KgcbHxyt19sU5EBPovH6OSxORGhjShEVifXSMJ40JK1asUOrsnxmX5nppIlIEQ9riNE2TIX1VDIuYqGwgx6WJSCEMaYvbdjSI5HgbmqTF/qVw/fXXQzVyXHo710sTkRoY0hYnll7pMR4trF69GqoR3f7L9nJcmojUwJC2uIU7Y1/EJOTw4cNQTWqCHW2yHFjB9dJEpACGtMXpUcQkJDk5Gaqul57H9dJEpACGtIUdLAzidLmGNln6vAyGDRsGFQ1szXFpIlIDQ9rCxIYSYmMJuz22RUxCpkyZAhWJHbFEdzf3lyYivTGkLUyUwBzYmkuvzicqr3Vq7JRL04iI9MSQtjCx1EgsOdJLfn4+VDWkjQtfbeVSLCLSF0Paog6cDqKwQkP7hg7djiE1NRWquibPia+3MaSJSF8MaYuPR9ts+oxHC0uWLIGqeuY4sfN4EMdLgnofChFZGEPawl3dHI++OJfDhqtznZjD1jQR6YghbVFiHbCoU62noUOHQmVyXHobJ48RkX4Y0ha092QAZV4NbRvo+/Rv3rwZKhuSVzV5TGxCQkSkB4a0RZdeiapaeo5HCwcOHIDK8rLsEPm8/RjHpYlIHwxpi04a03PpVYjH44HKxIeYIW04y5uI9MOQtqB52/UfjxZuuOEGqI7rpYlITwxpi9l9IgBvQJNduXqbNGkSVDc4zyWHB/wBjksTUezp/05NMW9FqzAebRT1k+xoUc+O5dy6koh0wJC2GDFb+Zo8/bu6hby8PBjB9e1c+Gwzq48RUewxpC0kENRkqcuhbdUI6aysLBjBiPZufLqRIU1EsceQtpA1BwJokGxH0zQ1nvZFixbBCHo2c+BIcVCuLyciiiU13q0pJmZt8eFaRVrRRiL22x6e78LMTWxNE1FsMaQtFtJD2+i/Pjpk0KBBMArZ5c2QJqIYY0hbRFGFhrUFfvRvpU5LeufOnTAKMY6/ZI8fxRVcikVEscOQttCuV32aO5HgVmfp1d69e2EUSXE2ef7E7HgiolhhSFuEiuPRbrcbRjKivRiX9up9GERkIQxpi1AxpMeOHQsjGdm+ar10MMgubyKKDYa0Bew4FkCFX0N+QwdUMmXKFBhJToZDLmFbvo9LsYgoNhjSFiBmJQ/PdytXCtTv98NoRJf3p+zyJqIYYUhbgBhHFV21qmnZsiWMZnRHNz5ez3FpIooNhrTJFZYHsXJ/QO7mpJpmzZrBaK5o5kBJJbD5MLu8iSj6GNIm9+V3PvRv5YRHoaVXIfPnz4fRiCGDsZ1dmLqOrWkiij6GtMnN3OxTsqvbyMZ1duMjhjQRxQBD2sT8AQ1ffFc1aUxF/fv3hxGJoibHS4PYfoxd3kQUXQxpExNlLHPS7WiiyK5X5ysoKIBRN9wY08nNLm8iijo1370pYkuvVO7qNlLt7vPd1MWND9dwXJqIooshbVKapuHTjV6M7KBmV7dgtxv35de3hejy1vAdZ3kTURQZ912SLmnT4QDKfUC3pmpVGTvbzTffDKMSXd63dnPjP6sr9T4UIjIxhrRJTVnrlV2yqlUZO9u0adNgZLeJkF7llb0WRETRwJA2IREak9d4Mb6rul3dQmWlsVuhnZs4EOcElu/lLG8iig6GtEm7uiv8QI9sdbu6jVpx7Gyil+K27nHs8iaiqGFIm5BoRave1S3k5eXB6MS4tJjlLdakExFFGkPajF3d349Hq+7rr7+G0bXKdKB5hh1ztxtvRy8iUh9D2mRWHwggEFS/q9tMfnBFHN5dZuzxdSJSE0PaZP69shJ39IhTvqtbuOqqq2AGYpb3l1t8OF4S1PtQiMhkGNImIsZFP1jjxe3d1e/qFo4fPw4zSPPYMaqDC/9eyQpkRBRZDGkTmbPNh+bpduTWN0ZX99atW2EWd/eJwz+WVnLNNBFFFEPaRP69yiu7uin2rmrhRDCoYekeTiAjoshhSJtEUYWGmZt8yhcwMUtZ0POJOQA/Ea3pJZxARkSRw5A2iUmrKnFNnhP1k4zzlM6cORNmm+X98QYfCss5gYyIIsM47+h0SWI89O4+8YY6S6WlpTAT8QHp+nwX3l3GCWREFBkMaRNYc8CP4yWabEkbSZMmTWA2D/WPx18WViAQZAUyIqo7hrQJiHHQH/eOk9snGkmHDh1gNj1znGiUYpfzA4iI6oohbXDFFZqsHf3DXsab1T1r1iyY0f+7Oh5//qZC78MgIhNgSBvcO8sqMaSNC03T+FSq4saOLuw6EcTaA1yORUR1w3d2g1cYEy22Rwcaa8JYSK9evWBGTocND/SLY2uaiOqMIW1gYrlPdpodVzQz1oSxkJKSEpjVT3rHYeZmH/adCuh9KERkYAxpA29J+fK8csO2ooVNmzbBrNI9dvy0TxxenMOxaSKqPYa0Qc3b7kdhuYaR7V16HwpdxMMD4vHBai8OFrK4CRHVDkPaoH43qxy/HppguGVXZxs3bhzMTBQ3uatnHF6aW673oRCRQTGkDeibHT4cLAriZgPV6a7O7NmzYXY/HxiP91Z4UXCarWkiCh9D2oB+P7scTw5JkLOIjayoqAhm1yjVLrexfGYWW9NEFD6GtMHM3ebD3pNB3Nbd2K1ooUGDBrCCxwfHY/oGL7Ye4UxvIgoPQ9pAxH7FP/+kDC+M9MBl8Fa00KNHD1iBmOktur2f/LxM70MhIoNhSBvIf1Z7keCyYUwnc8zo/uyzz2AVD/aLx6r9AczbzpreRFRzDGmDKPdqePKzcrw0KgE2m/Fb0VaT4LbhL2M8uG9KKbx+7pBFRDXDkDbQkqurWjhxZQtztKKF7t27w0pGdnAjL8uBP81ngRMiqhmGtAFsOOiXG2m8cqMHZuLzWa/r99UbPfjfeRXYe5KTyIjo8hjSigsENfz0w1L8YXgCGiSb6+lav349rKZ5PQceGRCPh6ZxEhkRXZ653vVN6I9zKuR45o8MuF80VU/UW99yNIBp67w8RUR0SQxphS3b48erCyrw3u1Jhi7/eTGjR4+GFcU5bfI5ve+jUhxgJTIiugSGtKJOlwVx279K8Ob4RDRNM+fTtGDBAlhVzxwn/t/V8bjjXyVySIOIqDrmfPc3OH9Aw83/V4JRHVy4oaPxK4tdzMmTJ2Fljw2Kh9MBPP0lS4YSUfUY0gp6ZHoZHHYbXhplrtnc56tXrx6sTAxhTPpBEv61wsvxaSKqFkNaMX+cU4652/2Y9INEw2+gcTl9+/aF1YntLKf9KAn3TinFugK/3odDRIphSCvk1W8q8I8llfjqZ8lITTD/UzNjxgy9D0EJ3bKd+NtNiRj2ZrFcE09EFOI88zfS1RuLK+RM7m8eSJbbG5K1jO3slhPIhr5RjNn3JqNjY16aRMSQ1p2maXjmy3K8t9KLOfclIzvdAavo0qWL3oeglPFdq9bCXyuC+mfJ6NCIQU1kdXwX0FFppYZ7Jpdi5/EAlv6/FGSZrKLY5djt1nq8NQ1qTQOG/q0Yn/80GV2a8hIlsjK+S+pk7QE/ur9cCLcDmHu/9QJaWL16td6HoKSbu8Vh4thE2fX97rJKvQ+HiHTEj+kxVubV8IevyuUEMbF14S3dWO6TLjSmsxv5DR0Y924JFuz04S9jEpEcb+7Z/kR0Ies133Qcexa1mtu/UIi9J4NY/1iq5QN6xIgRej8tSmvbwIFlD6fA5bDJ182MDaz1TWQ1bEnHIJw/2ejDM7Oqqkq9fUsiBuWZZ0/ouli+fDkGDx6s92EoLTHOhr/fnChb0/dOLpUrAH5/XQKuasnXEJEVsCUdJYeLgnjx63K0fa4Qv5tVjqeHJWDVoynnBHRlZSWefvpp+afViMe8cOFCyz72cJ/3/q1csvflzivi8IP3SzHk9SJ8tNaLSj/rfhOZmU0TTT0LKCoqQmpqKgoLC5GSkhKV31FwOojPNnvx6SYfFu/2Y2wnN37SOw49cxyw2Wy6HJOqxGO/5ZZb8MEHH1jysdflefcFNExZ68U7yyqxriCAcZ3dGNnBhYG5LrmtabRZ+XVLFGvs7q4lr1/DtmMBLN3jx5Lvb0dLNAxr68Jt3dyyJnNSHCf6XMqcOXNqe/otTYxR39Y9Tt72nAjgo3Veue/4re+Voke2A71ynPLWI9uJxqm2aj8gEpExMKQv0lI5Vfb9rVzDoaIg9p8K4kBhELtOBLD5cAC7TwbRPN0utxzs09yJCf3j0aGRQ26MQTUzduxYnqo6al7PgZ8PSpA3sb3p0r1+LNsbwJvfVuLeglKUVGpoXd+BvPoOtMq0o3GKHQ3FLdkml/2lJdiQEm+Te1wTkXosE9KBQED+efPfDyLoKpZjeZUBoNKnodIPVPg1VPir9nH2BoDUBBvS4m1I89jQMNmOpql2NE61Y0SuHT+/SrzhOb5/YwuKdrW8lZaE32149p9WIh6zz+ez7GM/+89ITjC5sknVDVeKrzhQVKHJYjk7jnnlB8t1u4OYVRzEkRINx4qDKKrUZJAHg5BLvFLibIh32eB2Qq7hdztEgFe13sWf4t+BikL5+6y+1ShRXYiR5uLiYjRu3PiShZ0sMya9YsUK9OzZU+/DICIiOmP//v1o2rQpYPWWdG5urvxz06ZNykx2EZ+i8vPzsXnzZiQnJ8NKxGN/8MEH8de//tWSj93Iz3tBQQF69+4t31xUuZao7r744gtcd911PJUxInrSsrOzL/seYJmQdjiqNq4Qn1hUeWMJdXc2adJEmWOK5WPPysqy7GMXjP7YxbEb+fjpXMFgkM+nDi43sZPrpEk3Bw8e5NknUoQYGyX1MKRJN+vWrePZJ1JE586d9T4EqgZDWkdxcXH47W9/K/+0GvGYH3roIcs+dqs+76T2mDSpxzJj0ioSb9KiPKRVH3ubNm0sGVRWft6JKDxsSZNuuCSOSB28HtXEkCbdlJWV8ewTKYLXo5oY0qSbjRs38uwTKYLXo5oY0kRERIpiSJNuxowZw7NPpAhej2piSJNuuFUlkTp4PaqJIU26KSys2k2JiPTH61FNDGnSjajdTURq4PWoJoY06YbrMonUwetRTQxp0s3MmTN59okUwetRTQxpIiIiRbF2N+mmW7duPPtEUXCsJIg1BwKo8GtonmFHuwYOuByX3reY16OaGNKk6ybzRBQ5J0uDeO7rCry7rBJdmzrgcdmw80QQXr+GV8d4cH2+m9ejwTCkSTdr165Fu3bt+AwQRcCWIwEM/Vsxhue7sOmXqWiY8t/RzK+3+vDDSaV4caSG27pXv/Mcr0c1MaSJiAxuz4kABr9ehBdGePA/V1wYwte0cWH2vcnyPnFOG8Z2vniLmtTCiWOkmxtuuIFnn6iOKnwaxv2zBI8NSqg2oEPaNXTgi3uS8bMppdh8OMDr0SAY0qSbRYsW8ewT1dGvPy9HbqYDE/pfPKBDOjdxytb2be+VoNKv8Xo0AIY06ebEiRM8+0R1sOaAH/9eWYmJYz2w2S49ezvkh73cyMmw4+V5FbweDYAhTbrJyMjg2SeqpUBQwz2TxWQwDzKTav5WLsL8L2M8eGV+hRzL5vWoNoY06aZ///48+0S1JJZZiSVWP7gi/ElgORkOPNgvHr/9spzXo+IY0qSb6dOn8+wT1YIvoMn10M+PSKhxN/f5Hro6Dp9v9mHX8arWNK9HNTGkiYgM5j+rvGhZz44+LVy1/hmpCXb87Ko4vDjn3LFpUgtDmnTTqVMnnn2iWoxFP/dVOX4zNKHO525C/3hMWeeVlcp4PaqJIU26cblq3wogsqopa73ISrajf6u616ISE85GtXfhnWWVvB4VxZAm3axatYpnnygMmqbhD19VyFZ0bceiz3d/v3j8bXElVq7k9agihjQRkUHM3+GHiOYhbSJX0fmKZk5kJtpwtIQb3qiIIU26GT58OM8+URjeWFwpJ3tFqhUdcn/feHwTGMjnQkEMadLNypUrefaJauhIcRBfbfPh9h6XL/8ZrnFd3Kg8uFb+DlKLciF91113we12Iykp6cxtyZIldb4vqefIkSN6HwKRYfzf8kqM7eRGSnxkW9GCx21Dx9QTmLTKG/GfTSYLaeG+++5DSUnJmVufPn0icl9SS0pKit6HQGSYCWP/XF4p625HS6vGaXhvZWXUfj6ZKKTJGoYOHar3IRAZwop9AQSCQJ/mkZswdr47x16Lo8VBbD1y4TaWpB8lQ/q9996Tmy+0b98eL7/8MoLBYETuKxQVFZ1zq6zkJ0e9fPTRR7r9bqo5cY2cfc0UFxfz9MWYaEXf2TPyE8bONm3aVNzYyY2PN7DLWyXR+1hWSxMmTMBLL70kg3fFihUYP3487HY7Hn744TrdNyQ7O/ucf//2t7/F008/DVWIiRvPf1WBU+VB+clZ3Ow2wOkA3A4bEt023NXTLYsZNEi2wSH+kyiKnn/+eTzzzDM8xzrW6f5onRcrHon+8NCYTm48/mkZfnlN3auZUWTYNDHYobDXX39dtpaXLl1ap/uKFkBqair2799/zlhoXFycvOmtpFLDvZNLEee0oUmqXa5bdDmqAlo8Q74g4PUDpV4NxZUaiiqqbuI+qfE2pCXY8NDV8aiXaIvqp+1IWr9+PUsRGqQlfXaPU0FBAfLz81FYWMh5BTHwxWYvnv2qAosfSon69ZjfviMaPXUaq3+egux0R1R/n9UVfZ9Jl7uOlGtJn0+0jCN5X3EyVJuwFAxqeHBqKdpkOdAg+cLHIDI3zg7EOYHkeBsanvV/Xr+GwgoNp8o0PP5puQxxEdoi5H81JAFJceoGtpiNT+o7/4OseHOh2Jm02otbu0VvwtjZ16PTYcOoDi5M3+DDg/0Z0ipQbkx68uTJ8k1ANPDFOtoXXngBY8eOrfN9VfaLT8qRGi+6r8N/OtxOG+on2ZGX5UCvHCcG5DrRKtMObwD42eRS3PV+CZ6cWSbDXDXLli3T+xCIlFbu1fDpJh9u6uKO2fUourynree4tCqUC+mJEyeiWbNmSE5Oxu233y6XWD366KPy/+699155q8l9jaLMq+FwURBtsyLzVNhtouvbjtb1HejXyoXOTZyo8AN3vl8qu9P3nAjIDzVEpL7PNvtwRbPqe9iiZXCeC2sOBHCilIVNVKD8mHSs+/9jbdvRACYurECHRtEdeRBP84kyDXtPBuUHg5b1HPjtsHjYdZx4dvLkSTnpj4zlwIEDcgKmateSGY3/ZwmGtXXhR72jP2/m7OvxxreLZev9tu76z9exeiYp15K2mje/rUR6QvSfBjGZLDPRju7ZTnTLduJkWRA/eL8Uv/uyXLeW9caNG3X5vURGUOHTMHuLT44Rx/p6vD7fhc83+2Lye+nSGNIKXIgJMd5WWSzj6tjYie5NnThcHMT//LsUv58V+7AWs4SJqHpzt/vQuYlD7vkc6+vxunZuzNriQyBoiY5WpTGkdSZmbut1GSTG2dClSVVY7z8dxI8nlcqKQzH7/YmJMftdREYzY4MPoztGf8JYdddj0zQ7GqfaZaUz0hdDWoGQ1vvDqgjrHtlO5NV34JHpZXjk49jMBh8xYkTUfweREYllmZ9s8uKGGHV1V3c9Dpdd3pzlrTeGtM5EFTGxXEoFGYl29G3hlAVSxGzwp78oj+rv+/DDD6P684mMavm+AOon2tEyM3Zrlc+/HsW4tJhdTvpiSOvs9u5uOS6tCjHbu1WmA1e2cKKgMCiXbYm1mkQUO9M3eDG6Y4wnq5ynd44Tu04EYzoERhdiSOssNcGGMoVCOkSUJ+2R7UD9RBt+NCk6reo2bdpE/GcSmSekYzceXd31KKqPXd3KiXnb2ZrWE0NaZ6JIQVG5eiEdWrYl6veK7fEOFgVx9welKK6I3LFmZmZG7GcRmYXYKrLMC3RtGtuynNVdj4NauzBnuz+mx0HnYkjrTNTWFmPSKteUiXdVTSwTMz5/+mHV2upIWLx4cUR+DpGZfLbZKydtxXqjnOqux8F5TszZxpa0nhjSigS12NFKdY1S7LJVvfdU1Vi1qFxGRJH15RYfrsvXdzw6JL+hQw7HiXLCpA+GtALqJ9lwrNQYgSda1T2bOeQx/+SDUuw/VfuL95prronosREZXWmlhqV7/LKbOdaqux5Fa14cy1x2eeuGIa2AXwyKN9QMSjlWneaQXeC/mlmOx2aUyXWd4dq2bVtUjo/IqObtEBtqOHXZYvZi1+NgOS7NLm+9MKQVkO6xIxAEyhWc5X0p4o3kqhZily1NtqpFKyAc+/bti9qxERnRF9/5cF07fbq6L3Y9hsalVZ43Y2YMaUWISVliXbLROOw2uYNXs3QH7v4wvO7vuDjusEMUIkJQz5C+2PWYk+GAx23D1qPGe38yA4a0In49NAEFp4OG/bQq6vyGur8f/6SsRo9jzJgxMTk2IiPYdjQIf6BqspYeLnU9ivXS3+xgl7ceGNKKEPWzxe2Uomumw+n+LvWK7u+yy87+ZllQonNndQ9rF/ulVzW5Hq9u5cI3O7leWg8MaYX8fGC83I3KyET3d6fGTmSn2eU49eGiiz+eYNDYj5UokvTs6r7c9Xh1rhPf7OS4tB4Y0grJzbTjVJmGyhjsQBWL7u9uTZ14dHoZfv1Z9d3frVq10uXYiFQjep2+3e3D4Dz9QvpS12PzDDvsNht2HucH61hjSCtEbG4hLgZR1N4MUuJtcqOOYyUaHphadsEG8k2aNNHt2IhUMn+HD92znfKa0culrkfRBS/HpdnlHXMMacX87roE2UUc7f2c75lcesEtGlwOG7o1dSDBabtgmdaCBQui8juJjNjVPaytvlXGLnc9yi5vTh6LOWfsfyVdith5pmU9O3YcDyC/YeSfnkuFcXX/9+b4xDr/TvEpPC/LgUNFNrlM60+jPWiYws+HRGeH9LQfJSl9QsTksWdnV+h9GJbDd0oF/e76BBwp1iLemq5NazmSrWxR+1uMUz8yvQy/n1WOAQMGROTnEhnZjmMBWcioYyN9ll6FXO56bF3fLufL7KtDKWAKH0Na0RnSORl27D6pzth0pMJajLmJzeS3Hwtg4uc7DbsunChSRDWvIXn6Lb2qaQVAcXx9WzqxeBeXYsUSQ1pRz16fgENFwYiVCo1UazgSYS026RC7aQVP78GEqbWr+01kFmLzCj021Djfrl27Lnufvi1cWLSbIR1LDGmFx6bz6juw5YiaXUt1DWvRW5DqccrH+dMPy6I+UY5IRaInSczsHtha/+lBTuflj+EqtqRjjiGtsN8Oi5ebV5wsVafbO5JhndH5RrRr4ECD5AtnfhNZwabDAaQm2JCdru94tHDTTTdd9j5dmjiw60QAheXqvieZDUNaYWIM6PkRHnx3JKD82G1twvrUhk/OFPDPzXTgpx+WKv2BhCjS5m5To6tbmDp1ao2WVIqtNJfsYZd3rDCkFSeWKiW4bThcrHZIh4QT1Jrfe87j7NTYgfs/KsNBA+4GRlQbc7f7MDBX/65uwev97/V4KWLy2CJOHosZhrQB/HGUR45N+wLGCeqahLU7PfuCfbXFp/SfzyjDs7PKo3iERPoTFfhEBa8BirSkc3JyanS/vi1dDOkYYkgbgFi2JAqcbFZ0EtnFXC6o4+q1qHYnLbFEa8vRAH7zeVkUj45IX2sOBOQ+8g2S1XgbrmktfXF9rjrg52TPGFHj1UGX9ezwBJRUajhZVruu4EhUDot0UBfvWHDJJVoHTgfxy08Z1GTeru5BCszqDpk7d26N7pccX7XyRHzIoOhjSBtoEpno9t50KICg4pPIIjGpTExQ6dnMiROlGh75uPpdtIiMbJ4i66NrQ+wbv2i3T+/DsASGtIGIbrF0jw17FapEFo7zgzqpRe/LrqXunu2Qy9AmTGNQk3mIugDf7vHLnaVU0bdv35rfl5PHYoYhbTCiNb33VLBWa4r16vK+WFD7S45f9v5iD9vOjR0QBRPFzG9WJyMzWLHPj7z6dqR51HkLPnr0aI3ve1ULFxbv9rOHKwbUeYVQjYjxWlGIf+1B414goe7vimM7atzV376hA3FOG+6dcuG+1ERGLAU6ULGu7m3bttX4vk3S7HKS57ajxuzVMxKGtAH9akiCvED2nw4asjUdsv5gzSeeiKAW1clS422458My+A2yHI2oOgt2+jBAkfXRdRmXFl32FF0MaYP631Ee7DoRRKnXuGE1qXRE2BPKWtd3oF6iDfdMZlCTMYl6B8v2+mXIqeTWW28N6/5iBca33Gwj6hjSBiWqkHUQ3d4H/Iab7R1yg+dr+We4Qd0q04H6STb8dHKZYQq8EIWs3h9Abn0HUhPUevudMWNGWPe/sjnLg8aCWq8SCssTQxKQkWjHjmNBQ3Z5e2wVZ/4eblC3rOdAgyTR9V3Krm8ylG92+pSa1R1SVhZeTQLRSBBDbqdrWbuBaoYhbXAvjarad7q4QjNcUO8PNDzn3+EGdYt6DmQl29n1TYayYKcf/VupNWlMaNq0aVj3F9vMXtHMgaV7OS4dTQxpgxMXygsjE7D+oPG6vTd7cy/4Wm1a1JmJNtw7mbO+SX1iZYJYuiTWGasmPz8/7O/huHT0MaRNQOxFK/ak3XHcWN3e1yYsqvbrtRmjzvAwqEl9YkVD41Q76iep99Y7e/bssL+H49LRp94rhWrllRs9OFocrHVtb6OXEhUTccRGJPdNYcETUrurW8Xx6Nrq3dyJ5fsCrF0QRQxpkxAlNF8e7cGGg+FdMHq2ppdUdrnsfcIJalHByeO2ycpkRi30QuZfH91f0ZDu06dP2N8jtpcVO3ltPMTNNqKFIW0iogstO904F0yqvbhG96tpUIuCJ22z7HA5gAlTGdSkFvHBUbSk+7VUb9KYUFhYWKvvE+PSS1jUxDohfdddd8HtdiMpKenMbcmSJdXe1+fz4YEHHkB6ejoyMjLw4IMPwu+39kzD50ckwBuA3OZR9dZ0vmtnje8bTlCLEqJi+fTDH5ezRU3K+O5IUM4dESU1VbR58+ZafZ8Yl2ZRk+hR8tVy3333oaSk5MztYt0wzz77LBYtWiRfXJs2bcLChQvx3HPPwcpESP1ljAc7jgfCqkam9ySySAd1p8YOlPs0/OKT8qgfF1GNu7oVnNVdV1eyPKj1Qrqm3nnnHfz6179Go0aN5O3JJ5/E22+/Das7U42sQO1lWVNKh0XtZ4vds7o2ceBUmYZffRpekQaiqE0ay1Wzq1u46aabavV9Yi6IuM7ExFWySEi/9957svu6ffv2ePnllxEMXvjknzp1CgcOHECXLv+dfCT+vm/fvkuOrRQVFZ1zq6yshFmrkaUl2LArjGVZsW5NX+tZGPb3hDORzG63oUe2QxZ7eepztqhrS1wjZ18zxcU1m0tA549HqztpTPjyyy9r9X3iOhOzvDkuHR3KvWImTJiAl156SYb0ihUrMH78eNjtdjz88MPn3E90gwtpaWlnvhb6u3gTSU1Nrfbni5/jcrkwdepUDB48GD179kS/fv3knzNnzpT36datm/xgsHbtWvnvG264QXarnzhxQh5X//79MX36dPl/nTp1kj9v1apV8t/Dhw/HypUrceTIEaSkpGDo0KH46KOP5P+JDx1ijH3ZsmXy39deey02btyIgoICJCYmYsSIEfjwww/l/7Vp0waZmZlYvHix/Pc111wjt5ITH0Li4uIwZswYeV9xnK1atUKTJk2wYMECed8BAwbI+/Wq2IXVB4FTiaOBnTOh+b1wp2cjrl4LFO+oum9Si95yX+fQtpFvjh+Hz2fOkCU7RUUwUXAktJ5ZzMYWk71CY8miJSyCNsVWikOBTKypbI/rPd/I/1tR2RFxtkp0cldtfzetdAgGJixDur0Ix4LpWFLRFe1cu5CSWIo13nbyPl3d38k/PykbhD7xa1DffgqngimYV94LYxK/kv+33puHl6bF4cctqsbPUtsOQdnB9fAVHYEjPhkpeYNxan3VcxPfoA0c8Sloc3oFTh0M4nfeQRhYb5f8cOfxeOTzOmnSJHnfvLw8ZGVlyedZGDRoEHbu3Im9e/fKORJjx47FlClT5JyHli1bolmzZpg/f768r3g9iOdQ3F+8Vm+++WZMmzZNhpu4n/jZX39dVaf8qquuwvHjx7F161b5b3Ff8borLS2Vz2GHDh0wa9Ys+X+9evWSr3MxlCOMGzdOrmUVQdmgQQP06NEDn332mfy/7t27yzka69evl/8ePXq0fD2cPHkS9erVQ9++fc/UZhYfZsVxrl69Wv5bvO6WL18u9xMW1424LsTxC+J4xLUS2sbwiy++kK9NCo/YDEcMwzTPULJdJNXlw1doXPqGju6IHhMBNk3xtSqvv/66bFkvXbr0gpa0CMwdO3bIkBLE31u3bo3Tp09fENLijU18bf/+/TI8Q0TgiZtZnSgN4sGpZXLHHZfDVqPvCbeYSG0NiF+K+RW9a/394bT8K/0alu7xo30jh+xloJoTHzbO7nESH0hEdSrRY3X2tUQX9+6ySny9zYf3/ydJ2dM0b948DBw4sHbfu92Hp74ox8IJfD3UVCiTLncdqfux7nviE391xIxuUWs21NoVxN+zs7Mv2ooWxMk4+2bmgBbqJdrlp/fNhwPKdXuLlnddhPNhIs5pQ89mTrk87VAhx87CIa6Rs6+Z5OTkWjxb1ia6ulUvYiJ6EGvrimZOrDngh9evdJvPkJQL6cmTJ8tPGKKBL7qNX3jhBdnVWJ0f/vCH+MMf/oDDhw/Lm5jZ/ZOf/CTmx6y6PwxPQGUA2H9araAOdY3HKqjFhLruTZ14ZHoZTpmkMhsZg6qbapzt888/r/X3JsXZkJclJqsao0aDkSgX0hMnTpTjeOLT+u233y6XYz366KPy/+699155C/nNb34jl2e1a9dO3sR43xNPPKHj0atJjIVNHOuRk8iKKsy1LCvcoE6Ot6FzYwce+KgMpZX81E/RJ2oWFFdqaJOl3NttRHG9tEXHpGPd/29mh4uCeHR61fi02D1L7/HpXOde7PDnROznhfOhQpwLsZb87VsTazxWT1XExDsxrGTlaykc/1lViWnrvfjoh2oPE4g5Pbm5F+5MF87jnL7Bh8l3qTvurhLTjElT5DRMsSMnw471hwI1rsQVzda0mP0dSeF8oBDnQuxGJFrUFvmcSrpuqqF2V7dQUVFRp+9nedDoYEhbcHxa2HNS//XToeVZegV1iww74pyifCiLnVD0qL4+OmTDhg11+n4xQdUf1LD/FMellQhpsYb5Zz/7GYYNGybX7D799NNn1miS+mVD958OhjV5yijj0+GWD81v6ECpF/j1ZwxqijxRhetwsSYrAJqduJ44Lq1ISH/wwQfo3bs33n33XRw6dAi7d+/Giy++iK5du8rAPnbsWOSPlCLG7bThzzd6sO5gAD6xE4VORIETvYnyoV2aOLD7RDCsTUmIamLhLr+cAyK2klXdjTfeWOefwTreioT0M888g7Zt28qJBuvWrcOaNWtkZaO//OUvctmUqJQkKm6RujKT7GiVaZdBrdf4tKhAFi3hdHuLiWM9sp34xYwynObSLLJgV7cQqqBXFxyXViSkRdUusTRKFBMJSUhIwP333y/LY4o3/V/96leRPE6Kgt9f75F7L2/Xqb63KBEaTeEEdWKcDR0bOTBhWhn8OvYukLkYZdJYqIpjXXVr6pSFk8rC2IGPohDSV1xxhSy9WZ369evLsA7VCSa1ifHpY8VBHCuJfVCLGt7RFk5Qi96FBsl2PDyd49NUd2LOx87jAXRtaozxaLFXQF3Fu8Q2sU6s3O+PyDFRDUNaVP0SmwCIsWfh+eefxz/+8Q/Zoq6O2CyCy1qMQYyV/XVcoiyXGev9p8UmG7EQTlDnZtpR4eNEMqq7xbv9cncoo6zDF4WhIvJzvt9sg2IY0mLPZrFjkFjoLiqBiZ2qRPe2qPX66quvyl19QkSQi7HpO+64I0KHSNGWEm/Di6MSsGq/P6yJZHUN6lGeuVBxhmpoIpkoeEJk5lKgZ/v0008j8nPE5DFuWxnjkBbb5YklV++88w7uueceubmFqJIitm4UWz+KbfNEBaLmzZvLbflEkLOGtrE0S3egRT07Vh+o+UQyIy3NCnciWdemTlmdrcLHsTUy76Ya0RBqSbM3NTJq9AqKj4+Xrebzd0kRIS32QxaL4EM3MctbzPAWy7ECAS5qN9pEsoemleK7I0G5fjjaQvtIxzKoa/qhQvQutMy0y20+/36zR7awiWqqpFLDpsMBuTuUUYj37EgQlfzEhhs7jgfRur4xxuNNW3FMbCZ/9dVX44EHHsCbb76Jb7/9VtYj3bVrF6ZPnx65o6SYeWW0R45N7zmp1o5ZerSos9McsiLZY5+UR/WYyHxEd6+Y6SwmUlmRXIrFcWl1y4KKbu+RI0dG40dTlNntNrw2ziMLexwpju6M767u76B6UHdo6JBVo56dxaAmc66PDhH1LiKFRU0ih7W7qdqKZK+NS8R3RwKm3Noy3A8t3bKdsuuSW1uSGddHRwPLg0YOQ5ouuu/yKzd6sPqAP6zJU+EE9Sdlg3Q7++G0phNcNrRt4MCDU0sRDHIiGV2auF7EdSOWXxlJJHs/OzV2yCGzwnKukKgrhjRdlCjskd/AgRX7/GFV4appUPeJj1z3WrSDulGKXU4m+wXHp+kylu/zo31Dh5w8ZSRLliyJ2M8S+9X3aMalWJHAkKZL+vW1CWiWbsfqggCCEV6aVd9e9zKEsdQ2q2p8uoAbcZCJ1keHnF3vIhL6tXTKDUaobhjSdFl/GOFBktsma/KG43JBfSqYovvZD6c1LcanOzdx4rFPyuD1s9ubqvfNDmOujxb1LyKpX0sXFjGk64whTTXyyo0JqPBD1iKOVFDPK++lxNkPJ6hFl3dOhh0PTWN9b7qQqNi3bG/V9pRGM2DAgIj+PDEmL8bmK/mBtk4Y0lQjopjH6+M8OFwcxMHCYESCekziV8qc/XCCOifdDjFE/5vPGdR0rlX7A8it70Cax3hvrR9//HFEf54Yk2/XwIGV+9jlXRfGeyWRbsRkkDduSsT2YwGcKI1MUBv1A0vnxg7sPB4M+zyQBdZHtzReKzpaRJc3x6XrhiFNYUlw2zBxnAfrD4a/vOL8oF7vzTN0fe9OjRx4+OMyBLgsi773jVgfnWu8SWNCx44dI/4z+3LyWJ0xpCls6R47/jJWrKEOv8DH2UFdqcUpd/bDCeqMRDvqJ9nx6HRWIyPID2uLd/nlrGYjEns0RCOkxWYb/CBbewxpqhURTn8clSA3dy/31S6or4jbYPiz37q+HYUVLBtKkL1LTdPsyEwy5tuq2OkwGu8TDZNtcr96qh1jvppICdnpDrlb1vK94c/gVHmMOqxlWWL/6cZObD4S4CxWi/tmh1gfbcxWdDT1a8WlWHXBkKY6eXJogiyZKZadhLt2uN+gYcqe/bDKhrptaJ5hxyMfc7a3lS3YZbxNNc52/fXXR+XnsqhJ3TCkqc5+c20CcjMdshyiWCdaU2UH1yvdog6HWJZV6Qee+ZLj01YkarobtdJYyOrVq6MY0j5oYVQspP9iSFNEPH1dAlrWq2pR1zSofUVH5J8iqFUM63Ba02JZlthUYMuRAMq8fDOyGrFjXHqCDY1TjfuWevjw4aj8XFH8RwwL7T7B5Yq1YdxXFCkZ1C0yah7Ujvjkc/6tYlCHI95lQ179qmVZZMGtKQ269CokOfnc6zFSxAfYvi1Yx7u2GNIUUc9c/9+gvtwYdUre4Au+plpQh9OaFhqn2iCWTbPb23rro41exGTYsOjNEenXiiFdWwxpikpQh7q+LzXr+9T66dV+XbWgDsfZ3d7cS9caxFirqDR2da6xQ3rKlClR+9liXHrRLl/Ufr6ZMaQpal3feVkOLN0T/jpq1YI63NZ0nNOGDo0c+H8fl3GyjAXsOB6UFehyMhx6H4qyxP7ax0o0HCnmuHS4GNIU1VnfYh21aFFXV5ksvkGbS36/ShPKwg1qUcQhwWXDU19wtrfZVc3qNnYrWsjPz4/azxbbvIrW9PztbE2HiyFNUV9H/dIoj1yedX73ryO+ZvtJqxLU4cpv4MCOY0GUhFk6lYzFqPtHny81NTWqP39Qaxfm7eCOWOFiSFPUNUmz469jPVhTEMDxkv8GdenempchVCGow21Nu502tMly4JHpnO1t5vFoETwDWxt7ZrewZMmSqP78ga2dmMeWdNgY0hQTop7x327yyPKZBWHuR61SUIdLrJsVNRx+xyInprT9WBB2G9CyHt9KazIufbJMQ8FpjkuHg68sipnUBDv+cXMi9p0KyD2pk/MGhv0zVBqnrikxiUwUuwi3vjmpb+52Hwa1dspZ/UY3dOjQqP58MS4tehzm7eC4dDgY0hRTos71W7ckynHaPTu2yHKKtaFXUIfb5S2ICWSi6tLjn3ASmdnM3e7H4Dzjd3ULmzdvjvrvGJjLLu9wMaQp5hx2G16/yYP4ioNYWosdtIzYqm6ebseJ0iAO1rKrn9QjPmCKMdaBBq80FnLgwIGo/w7RkhYfbKjmGNKkC9E92KNlsix6smTPhTO/wxHroK5Na1p09YkiJ499Ulbr3gNSy8bDAdRLtMmJkWbg8Xii/jvaZImNaDTsOcH9pWvKHK8uMqQbbrgBTw1LwJ9v9GBtQaBOE0qM0KoWY/IZHht+8zm7vc1g7ja/XFZkpusxFh/Oq8al2ZquKYY06WbSpEnyz6xkuxynPlAYxKbDfgTrsKWd6kEtNuDYdSLISWSmmTTmMt31GItxaXHuqGYY0qTMhLJ/3OyBw2ar0eYcRm1Vi7XTzTPsstubjMsf0LBwlx8DDF6vWw/X5Lnw9TbuL11TDGnSTV5e3gXjtn+60YPmGQ58W8dxapXDWsz0PlWm4RAnkRnW6gMB5KTb5fp/s16P0dK8ngPJcTZsPMRx6ZowzyuMDCcrK6var/92WAJe+X6cWqyprivVgtpus6FjIwce/5STyIxqzraq9dFWuB6jYWgbF2ZvZZd3TTCkSTeLFi266P81SLbjndsScbREw7oCPwJ1nBEdalWrEthiEllKvA3PfFmh96FQLYgxVbOsj67J9RhpQ9u6MHsLQ9rQIV1eXo7c3FykpaVd9D533XUX3G43kpKSztyiXX+WYkds+fjmeA+S4mz4drc/YhtV1DWsIxX0res7sON4gEuyDKbCp8l5E/1amSukY2lArksuvSz3cjmiYUP6qaeeQk5OzmXvd99996GkpOTMrU+fPjE5Pqq7QYMG1WjJxoujPHh+pAer9vsj0v1dl9Z1JFviohJZ/SQbnvyMS7KMtjVllyZO2RNitesxUsS5E+dQTL6jS1NyUGXVqlX48ssv8fLLL2P8+PF6Hw5Fyc6dO9GgQYMa3bdpWlX394NTy3C81C/HdF2OyL1J6tUNLlrTi3ZVVV0TPQekvi+/82JYO5elr8eIdXlv9ck/yUAtab/fj7vvvhuvvfaa7Mq+nPfeew8ZGRlo3769DPVg8NIzgouKis65VVZWRvDoKRx79+6tVfd3ZqJNBtuxs7a9NCrxQaNZuh1PzFS3NS2ukbOvmeLiYljZrK1+DDNhsIR7PdbVtRyXNmZIv/TSS+jatSv69+9/2ftOmDABW7duxbFjx/D222/j1VdflbdLyc7Olpubh27PP/98BI+ewlGTD2HVdX8/O9wja39vOxaQu0vVpfiJCsS6aVHTW4x1qkhcI2dfM/n5+bAqMdwi9kTv0sQBs6nN9VgX3Zo6cLAoyKWIl2HTxK7litixYwcGDx6MNWvWyNbx/PnzMXr0aJw+fbpG3//666/LlvXSpUsv+D/RAhBvMPv370dKSsqZr8fFxckbGY+ogf3ojHKcKguia1OnHOM1KrF1p8MOvDAy+vWTa9OSPrvHqaCgQAZ1YWHhOdeSFfxjSYUck/7XHUl6H4op3PzPEgzPd+EHPa33Hlz0fSZd7jqyq7YE4MiRI3JRfWZmpqwlKx6I+PuyZcsu+/12++UfjjgZZ98Y0PqZMmVKnb5fFD8R66lzMx1Yusdv6B2mRGGMfaeCdV5qFg3iGjn7mklOToZVzdriM2VXdySux7qMSxOMEdJikphoTa9du1be3nrrLfmGIP4uusDPN3nyZBniojNg5cqVeOGFFzB27Fhdjp1qN/8gEn59bUJV7e/TQaw/WPc11XqVCxVj7b+bxXXTqvIFNMzZ5jftRKdIXY/hGNLGia+2+rgM0SghLbZKa9q06Zlb/fr15Rik+LsYL7n33nvlLWTixIlo1qyZDPLbb79dLsd69NFHdX0MVHMtW7aM2OlKjLPhH7d4kOi2YfFuP07XsaSoHkQ5VG7hpy6xNrpVph31TVQKNFrXY001S3fI87liH0uEGmoJVsiAAQPOGY9+4403zvn/BQsW6HBUFCniA1YkhdZUHykO4ufTy5CVrKF1fbssw2kEomiLGJcWxy8qrpFavvzOJ2ckm1Wkr8eauqGjCzM2etGrudJxpBu+E5BuxMTAaBAB98/bE2W3txirLjNQVaPGqXa8NIdd3iqaucmHEfnmDeloXY+Xc0MHN2Zs5Lj0xTCkyZQcdhv+PCZRFgsR3ZT7Twfk3AXVZSXZZUua1Ft6dagoiJ45bO1FWo9shxyeEisc6EIMadJNTdbC19WTQxPwzq2JOFqsye0FxeQflcW7bBBHyJrGavl0o08uFRIf/swqFtfjxVZpjGrvxowNXl1+v+oY0qQbsd42FhLcNrwx3iPrZH9rgElloq7xsVK1j9FqPtnoxcgOsS32Ydbr8eLj0uzyrg5DmnStFRwroUplYl31hoPfVypTdKlWcpwNbyxmuVpVFFVoWLo3gCFtzDseHevr8XwDW7uw8VDAFKV+I40hTbqpSfGZSMv6flKZ2Jvj2z1+lEZo+8tI8rhtKGejQhli3+OrWjjl7Hsz0+N6PLsuv/gQJCbn0bkY0qSbm2++WZffK8YV/3e0B+0aOLBiv79qT2eFJpWJDxBGLMhi5q7uUR3M3YrW83oMuaGDi+PS1WBIk26mTZum69n/1ZAE/N/tifD6gcW7/CiuUCMYxVppxee3WYY/oOGL73wY0d7c49EqXI/X57swb4exlkzGAkOadKPCNqFiq8g/j/HIjS1WH/Bj1wn9l2oFglWtadKfGBIRW4mK/czNTu/rMd1jxxXNHPh6G7u8z2b+Vx4pS68KR9VpklY1Vl3qFZOE9G1V+4KAmymthClrvRjX2fytaFWuR1HYZPp6LsU6G0OadCN2O1OJaFX/dWwinh6WgDUFfmw9qs8McLGvtJs1M3Qnnvup67y4qYs1QlqF61EsxZq52SeHGagKQ5p08/XXXyt59nMyHHKsWli024+jMa4AJpb8/Owq6+2vqxqxUYsoMZtb3wErUOF6FBtutKxnxxx2eZ/BkCa6yAzwl0d78OoYD/aeCsoucBGesSB+T8MUXpoqdHWPt0grWiW3d4/D+6vY5R3CdwLSzVVXXaX82a+XaMc/bknEb6+Nl3tVryvwy+7oaBFj4S5HVdc76dvV/ZGFurpVuh5v7urGp5t8nOX9PYY06eb48eOGOfuhLvDMRLtsVW885JeTzCJNzC5vWc8a3asqW7jLL3szrNLVrdL1KAoO9c5x4tONbE0LDGnSzdatWw119kVp0WeuT8C/7khEeoIdaw74sWKfX+5aFYliKGJ96OlyDb8eGh+R46Xa+9dKL+7obp1WtGrX4+093Ozy/h7nkBLVYrxahLVYT32wUMNzX5dj82ENWck21E+0o16iLezdkgrLg1h3MID2jRxyVyDSj9iBbNp6L579ZSqfBh2XYj3wURlOlAblkJOVMaTJsmUII9GybpJmw2vjRNUyDbtOBPHawgpsPqLJYiRiNytxS3DZEOesqk8sqonJsp8aUOmvaj3vORWUgf+H4QmyW530NWOjV3a3Wm3ynkrXY3K8Dde1c8nJe/deZe2eJYY06WbmzJkYNWqUKZ4Bt9OGtg0c+Ou4xDOtMdEN/ua3lThVrsnJZpWBqglJoqKY2Msg3lkV3r+/LkEWUyE1vLfCix9cYa2ubhWvR/EcPDOrnCGt9xNB1lVaWgqzEntYN6/nwPMjPXofCoXhUGHVcrupP0yy3HlT7Xoc2taFeyaXyVUVnRpbN6r48Z1006RJE559UspbSytxS1e3/JBlNapdj2Jex497x+EfS/Sv8a8nhjTppkOHDjz7pAxfQMOb31bg/r7WrPam4vX4o15u/GeV19JrphnSpJtZs2bx7JMyZmzwITfTgfaNrNm1quL1mJ3uwJUtnHICmVUxpImIALy2qAIP9LP2TGIV3XNlHF5fVAGrYkiTbnr16sWzT0rYdMiP7ccCchcmq1L1ery+nQsnyjQs3eOHFTGkSTclJSU8+6SE1xZV4qd94i1dM13V69Fut+HBfvF49RtrtqYZ0qSbTZs28eyT7sRWpB+u8eKnV1pzwpgRrse7eroxa6sPBadju22sChjSRGRpf/6mArd2c1uuwpiRpCbYZS11MW/AaviqJN2MGzeOZ590daosiL8vqcRjgzlhTPXr8dGB8fK5EvW8rYQhTbqZPXs2zz7pauLCSoxq70KzdNZMV/16zMlwYFxnN16eZ63WNEOadFNUVMSzT7o5WRrEXxdW4FdDEvgsGOR6fGJIvKyHf7zEOq1phjTppkGDBjz7pJtnZ5fjpi5utK7PVrRRrsdm6Q7c3NWN/7VQa5ohTbrp0aMHzz7pYtfxAN5b6cVvr2Ur2mjX4xNDEmQ972MWaU0zpEk3n332Gc8+6eKJz8rx8NXxyErmW6DRrsemaXY5G/+ludZoTfMVSkSWsmyPH4t2+fDwAM7oNqpfXZOAd5dVYs+JAMyOIU266d69O88+xXynq3sml+LFkR54LLgdpVmuxyZpdvkh65EZZTA7hjTpxufz8exTTIkuUvEGf1t3N8+8wa/HRwfGY8PBAGZvMdZxh4shTbpZv349zz7FzHeHA3hlfgXeuMkDm42taKNfj3FOG14d48EDU0tRbuL9phnSRGR6FT4Nt75XgudGJMg9iskcrs9344pmTvz683KYFUOadDN69GiefYqJn88oQ5ssB37S29qbaJjxevzLGA8+WFMpJwOaEUOadLNgwQKefYq6KWu9+OI7H/5+M7u5zXg91ku0442bEnHn+6WyipzZMKRJNydPnuTZp6hauc+P+z8qxbQfJcmdlMic1+PIDm7c2MmN2/9VikDQXOPTfNWSburVq8ezT1Gz/1QAN75TgrdvSUTnJk6eaZNfjy+MSEBlQMPTX5prfJohTbrp27cvzz5FxaHCIAa/XozHBsXLVhaZ/3p0Omz48M4kvLfCi4/WemEWDGnSzYwZM3j2KeKOFIuALsKPe8Xhwf6sKmal67F+kh2f/CRJLsv6YrM5gpohTUSmseNYAFe9WoTbu8fh8Wu4eYYVdW7ixIwfJ+Ou/5Tiq63Gn/HNkCbddOnShWefIubb3T70/2sRfjk4Hk8OZUBb+Xrs1dyJj3+chP/5dwkmraqEkXE2BenGbudnRKo7TdPwp/kVsuTnP29LxLB2HIPm9Qhc2cKFOfcl47o3S7DzRBBPXBMPu914leb4Lkm6Wb16Nc8+1XmC2I1vl8i10MsfSWFA83o8R/tGTix9OAWztvgw6q0SnDDgOmqGNBEZTjCo4Y3FFej8UiE6NHJgwYMpaMZyn1SNxql2zL0/Wb5OOr5YKD/Qid4Xo7BpRjraOigqKkJqaioKCwuRkpKi9+EQgOLiYiQnJ/NcGMyBAweQnZ2t27U0f7sPj31aBrfDhjfHe2RrierOCtfjsj1+/OiDUuTVt8vtSvOyHMpnElvSpJvly5fz7FONW86fb/bimteL8JMPS/HzgfFYOCGZAc3rMewJZat/noLeOU70/UsRfvifErkiQGUMaR1VVlbi6aefln9ajXjMCxcutOxjt+rzHq6jxUFMXFiB/BcK8dQX5fhhzzhs/mUqxneNi/l2k2Z+3sx6PVZW85yJLS7F8rztT6aieYZDLtkb8nqR7AZXcctLdnfryMpd8OKx33LLLfjggw8s+diN/LxHu7t736mAnOgzeY0XK/cHMDzfhXuujEPflk5d94E2+vNmxeuxqAbPWaVfw/T1Xry1tBLL9wUwqLUTN3Rw49q2LjRKtev+euJgDulmzpw5PPsWV1qpYeOhANYf8mPpHj/m7/Cj1KthUGsX7usbj2FtXUhwG2/ZjBFZ9XqMc9pwc7c4eRO7aIkd02Zs9OEXn5QhKc6GPs2d6JXjlBPP2jd0oEGyLaYfFhnSpJuxY8fy7JuUmI9a7gNOlWk4XR7EiTINB04Hse/Uf2/bjgXk19o1cKBzYwe6Zzvx6MAEtGtg17XFbFW8HoGMRDtu7xEnb+I1vON4EEt2+7F8nx/TN3ix6XAAgSDQop4dzdLtaJZW9aeYQZ7hsSPdY0N6gk3+mZZggyMC67ItE9KhSeyii0EVoWNR6ZhiRTxmn89n2cd+9p9GEzru6/9agICrCF6/hkp/VbdhpQ8o92soLNfgdEC+UaXGV715NUm1o2maHa3S7BjQzC7HA3MzHXJjBECsXxW1lr0oLoaSjP68WfF6LKrjc9YgDhjdtupWxYETJUHsOy0+YPrkh8xdh4JYvDWI02VBnK7QcLpMw6lyDYUVGhJcQLzLhninuAEJ4u+uqj/tvhL5Ey+3wMoyY9KhcTQiIiJV7N+/H02bNr3o/1smpIPBIA4ePCjXAbIrjaj2AoEAduzYgdzcXDgc+q0zJTIyEb1ibXrjxo0vWSLZMiFNRERkNFwnTUREpCiGNBERkaIY0kRERIpiSOtALHV44IEHkJ6ejoyMDDz44IPw+/0Xvf/OnTtx3XXXyfs3adIEf/zjH2GVx15QUIDRo0ejXr16yMzMxPjx43Hs2DGY/bEnJSWdc3O5XOjUqVPMj9mqysvL5cS4tLS0i97nrrvugtvtPud5WrJkCVQUzrGGe40a6bHdZaDnLIQhrYNnn30WixYtwubNm7Fp0yZZM/e555676EzaUaNGoVu3bjh69Cjmzp2LiRMn4j//+Q/M/tiF+++/X/65d+9e7N69GxUVFZgwYQLM/thLSkrOubVr106WbaTYeOqpp5CTk3PZ+913333nPE99+vSBqmp6rOFeoyq4L4znwUjPmSRmd1NsNW3aVJsyZcqZf0+ePFlr1qxZtffdtGmT5nA4tMrKyjNfe/rpp7Wrr75aM/tjFzp27Ki9//77Z/7973//W2vfvr1mhccesmzZMvkaKCgoiPIRkrBy5UqtQ4cO2qxZs7TU1NSLnpQ777xTe+ihhwxx0sI51tq+To3w2O400HMWwpZ0jJ06dUoWVunSpcuZr4m/79u3TxZar259t3D2SjnxtfXr18Psj1145JFHMGXKFPn/p0+fxqRJkzBy5EhY4bGHvP3223K4Q6ynpOgS3bp33303XnvtNdktejnvvfee7BJu3749Xn755TPXq4pqcqx1eZ3q6b0wngcjPWeS3p8SzGT48OEiSS962717t7Zv3z7592PHjp35vqNHj8qv7d+//4Kf6fV6tVatWmmPPfaYVlFRoW3cuFF+0hUtK7M/dmHbtm3alVdeqdlsNnkTfy8sLNSs8NiFkpISLSUlRZs+fXqMHo21Pffcc9qPfvQj+fd58+ZdsiW9atUq+Rz6/X5tyZIlWnZ2tvanP/1JU1FNj7W2r1M9rQrjeTDScxbCkI4gER7ixX2xWyAQ0E6ePClf8Dt27Djzfdu3b5dfO336dLU/VwTzkCFDtHr16mmdO3fWfvOb32hZWVma2R+7+J6cnBztF7/4hVZaWipv4u+9evXSrPC8C++++67WsGFDzefzxejRWJd4PkS37okTJ2oU0ud77bXXlHtthnustX2dquS1MJ4HIzxnltlgIxZqsg+rmDEp6rSuXbsWrVq1kl8Tfxd1xcXeotUR3TKzZ88+8+/HH38cV199Ncz+2E+ePCknjImJYh6PR35NzDR96aWXcPz4cTnb28zPu/DWW2/hzjvvhNPJSzXaxGSpI0eOIC8v78wsZ1G2UbzOPvvsM/Tq1euS33+p0o6qudix1vZ1qhJ7GM+DIZ4zvT8lWJFoCXft2lU7dOiQvIm/P/PMMxe9/7p162S3p5g8NnXqVC0zM1N+zQqPPTc3V/vlL3+plZeXy9vjjz8uu/ut8Ni3bNkiu/hFlz9Fn+ipEV26oZuYMCWGGsTfz564GfLhhx/KXpRgMKitWLFC9vr88Y9/VPKpCudYw32dGumxfWig5yyEIa0DMc583333aWlpafL2wAMPnNOdec8998hbyJNPPqllZGRoHo9H69Onj7Zo0SLNKo9dzG4fOnSofPzi/gMHDtRWr16tWeGxi679/v3763S0dH539/nPT79+/eT/JyYmanl5edqLL74ohzZUdKljPf9xXe51auTH1s9Az1kIN9ggIiJSlAE65ImIiKyJIU1ERKQohjQREZGiGNJERESKYkgTEREpiiFNRESkKIY0ERGRohjSREREimJIExERKYohTUREpCiGNBGRSZSXl8PhcOC5557T+1AoQhjSREQmsXnzZgSDQXTo0EHvQ6EIYUgTEZnEhg0bzuxBT+bAkCYiMlFIezwetGzZUu9DoQhhSJNuhgwZgj59+mDJkiUYMGAAEhMTkZubi88//1z+v/izd+/e8utdunTBqlWr+GwRXSak27Vrh6VLl565ptq2bYtJkybxvBkU95Mm3TRo0AD16tVDSUkJfvzjHyMtLU1OeKmsrMSLL76Il156CXfffbccYxNfF62DdevW8RkjuohGjRohNTUVZWVl8toR19Trr7+OLVu2yODu1asXz53BOPU+ALKmo0ePypvNZsOaNWvkm4tgt9sxYcIE/PnPf8bq1auRkpIiv378+HG88sorMsDj4uJ0Pnoi9Yhr5PDhw/B6vbLXqXnz5vLro0ePlj1U77zzDkPagNjdTbpYv369/PPpp58+E9BCUlKS/FO0okMBLYjWgQhwcSOiC23cuFH++eSTT54JaCE7Oxtt2rTB7t27edoMiO94pOss1FGjRp3z9a1btyIhIUGOV59t27ZtaNWqFVwuV0yPk8ho19Qdd9xxwf+JD7diQhkZD0OadGtJixZ048aNz/m6GHMWazzP79IWX+/UqVOMj5LIWCEtrqesrKxzvi7Gp8WYNK8fY2JIk24h3blz5wu+LsL4/K/7fD7ZwuabDNGlQ1pUGzvfu+++K+dyjBkzhqfPgBjSFHOBQEBWRjo/jMXEl0OHDl3w9e+++04GNUOaqHqapskx6YKCAuzcufPM148dO4YXXnhBDiuJZYxkPJzdTTG3fft2VFRUXBDGoeVV5389NMmMIU1UvT179siljN27d8fw4cPlCgnRzT1x4kQ5Hv33v/+dp86gGNKk2wSX80P3YmEs7i9mep89Y5WILrym3nrrLXl74okn5PJG0YIWLWlRk4CMicVMiIiIFMUxaSIiIkUxpImIiBTFkCYiIlIUQ5qIiEhRDGkiIiJFMaSJiIgUxZAmIiJSFEOaiIhIUQxpIiIiRTGkiYiIFMWQJiIiUhRDmoiISFEMaSIiIqjp/wPwCFd/0hcd7QAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mc = ns_line.to_getdist(labels=['m', 'b'], label='nautilus')\n", "g = plots.get_subplot_plotter(width_inch=5)\n", "g.triangle_plot(mc, filled=True, markers={'m': theta_true[0], 'b': theta_true[1]})\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b7d879f8", "metadata": {}, "source": [ "---\n", "## 3 Model comparison: is there curvature?\n", "\n", "Compare a straight line ($M_1$: $y = m x + b$) with a quadratic\n", "($M_2$: $y = a x^2 + m x + b$). The **Bayes factor**\n", "\n", "$$\\ln B_{21} = \\ln \\mathcal{Z}_2 - \\ln \\mathcal{Z}_1$$\n", "\n", "weighs goodness of fit against model complexity: the extra parameter $a$ spreads prior\n", "mass over values the data don't support, which is penalised (Occam's razor).\n", "A rough guide (Jeffreys/Kass–Raftery): $|\\ln B| < 1$ inconclusive, $1$–$3$ positive,\n", "$3$–$5$ strong, $>5$ very strong.\n", "\n", "We test two data sets: one generated from a straight line and one with mild curvature." ] }, { "cell_type": "code", "execution_count": 5, "id": "edb1e08e", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T11:44:59.212227Z", "iopub.status.busy": "2026-10-14T11:44:59.212148Z", "iopub.status.idle": "2026-10-14T11:45:24.027878Z", "shell.execute_reply": "2026-10-14T11:45:24.027455Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "data ln Z (line) ln Z (quad) ln B_21\n", "straight line (a = 0) -52.46 -54.83 -2.37\n", "curved (a = 0.08) -55.96 -47.58 8.37\n" ] } ], "source": [ "x = line.x\n", "sigma = line.noise_sigma\n", "\n", "def make_loglike(y, model):\n", " def loglike(t):\n", " mu = model(t)\n", " return -0.5 * np.sum((y - mu) ** 2) / sigma**2 - 0.5 * len(y) * np.log(2 * np.pi * sigma**2)\n", " return loglike\n", "\n", "linear = lambda t: t[0] * x + t[1]\n", "quadratic = lambda t: t[0] * x**2 + t[1] * x + t[2]\n", "\n", "prior_M1 = UniformPrior({'m': (-3, 3), 'b': (0, 10)})\n", "prior_M2 = UniformPrior({'a': (-0.3, 0.3), 'm': (-3, 3), 'b': (0, 10)})\n", "\n", "datasets = {\n", " 'straight line (a = 0)': -0.8 * x + 5.0 + np.random.normal(0, sigma, x.size),\n", " 'curved (a = 0.08)': 0.08 * x**2 - 0.8 * x + 5.0 + np.random.normal(0, sigma, x.size),\n", "}\n", "\n", "rows = []\n", "for name, y in datasets.items():\n", " r1 = NestedSampler(make_loglike(y, linear), prior_M1, backend='nautilus', seed=0).run().get_result()\n", " r2 = NestedSampler(make_loglike(y, quadratic), prior_M2, backend='nautilus', seed=0).run().get_result()\n", " rows.append((name, r1['log_evidence'], r2['log_evidence'], r2))\n", "\n", "print(f\"{'data':24s} {'ln Z (line)':>12s} {'ln Z (quad)':>12s} {'ln B_21':>9s}\")\n", "for name, z1, z2, _ in rows:\n", " print(f\"{name:24s} {z1:12.2f} {z2:12.2f} {z2 - z1:9.2f}\")" ] }, { "cell_type": "markdown", "id": "98d2a371", "metadata": {}, "source": [ "For straight-line data the Bayes factor favours the simpler model even though the\n", "quadratic fits at least as well; for curved data the quadratic wins clearly. The posterior\n", "on the curvature parameter shows why:" ] }, { "cell_type": "code", "execution_count": 6, "id": "1facf860", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T11:45:24.028973Z", "iopub.status.busy": "2026-10-14T11:45:24.028893Z", "iopub.status.idle": "2026-10-14T11:45:24.104831Z", "shell.execute_reply": "2026-10-14T11:45:24.104517Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAFUCAYAAAA5ywiWAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAARgFJREFUeJzt3Qd4U2X7+PG7dAFllL1X2RvZKMgGcTBEkSUgvCACIiAq+CpLX0H5KS4UF+AGFEEFRYYVAQE3yB6yQUD2aoGS/3U/+SckbVqSNmlGv5/rOpzk5OTkec45TW6eGWaxWCwCAACANGVL+2UAAAAQNAEAALiJkiYAAAA3EDQBAAC4gaAJAADADQRNAAAAbiBoAgAAcANBEwAAgBsImgAAANxA0AQAAOAGgiYAAAA3EDQBAWDChAkSFhaWaZ83e/Zs83l79+6VrH7+Qu1cZBTnA0gdQROynM2bN0vv3r2lRIkSEh0dLcWLFzfPt2zZIqHmp59+MgHF6dOn/Z0UeAnXFPAfgiZkKV988YXUrVtXVqxYIQ888IC88cYbMmDAAPn+++/N9i+//FJC7Qd24sSJKYKm+++/Xy5duiRlypTxW9rg3WvqLdwbQOoi0ngNCCm7d+82PwhxcXHy448/SqFCheyvPfLII9KsWTNT4rRx40YpV66cBKILFy5ITExMho8THh5uFiD5veXte8Nb9ywQCChpQpYxdepUuXjxorz99ttOAZMqWLCgvPXWW3L+/Hmzn02/fv2kbNmyN2xDs2/fPhkyZIhUrlxZcuTIIQUKFJB7773XZTuZ1atXS4MGDSR79uxSvnx587mu2D5Dqw179uwp+fLlk6ZNm7r9efr+xx57zDzWIFCPZWu7k1q7lUOHDpmSN62y1KpLfd9DDz0kly9fTvW82tK5Y8cOE3TmzZvXnN+nn35aLBaLHDhwQDp16iR58uSRokWLyosvvujyOH/88Yd06NDB7JcrVy5p3bq1rFu3Lt3nT/PSv39/KVKkiMlL9erVZebMmeKObdu2yf79+2+4ny3vun+3bt1M2vVaaBCekJCQrjyeO3dORowYYe47TXfhwoWlbdu28vvvv6d5TT3Jd1r3Vmr3hjtpT+u4qfHkuLt27TJ/k7GxseY+09Ji/Zt2x9atW2Xw4MHmbyZnzpzmXuzVq5ccOXLErfcDipImZBlff/21+SHSEiVXbr31VvO67qfVdp745ZdfTLVJ9+7dpWTJkuYH580335QWLVqYHxD9klZ//fWXtGvXzgQV+kNw9epVGT9+vPmBS40GQxUrVpTnnnvOBCHuft7dd99tAplPP/1Upk2bZgJDlTxgtDl8+LA0bNjQVPsMGjRIqlSpYn6AP//8c/PDFBUVleY5uO+++6Rq1aoyZcoUWbx4sTz77LOSP39+E9S0atVKnn/+efn4449l9OjRJujR8+3Yzkyvi/5wPv744xIZGWnep/lZuXKlNGrUyKPzd/ToUWncuLH5oR02bJjZ/9tvvzUB4dmzZ01QkhbNR/PmzeWHH34Qd2jApPfO5MmTzQ/+q6++KqdOnZIPPvjA4zzqD7uec013tWrV5MSJEyZQ1B/9G11TT/Pt6t5yxd20+/q4ep41WNTzrEHku+++a4JKvbduRAPHPXv2mNJmDZi0RFk/S8/Z8uXLb/h+wLAAWcDp06f1m9vSqVOnNPfr2LGj2e/s2bPmed++fS1lypRJsd/48ePNfjYXL15Msc/atWvNPh988IF9W+fOnS3Zs2e37Nu3z75ty5YtlvDwcKfjOX5Gjx49Uhzb3c+bOnWq2bZnzx6nfWfNmpVie58+fSzZsmWz/PLLLymOfe3atRTbkqdz0KBB9m1Xr161lCxZ0hIWFmaZMmWKffupU6csOXLkMOfVkZ6XqKgoy+7du+3bDh8+bMmdO7fl1ltv9fj8DRgwwFKsWDHLv//+6/Q53bt3t+TNm9fp/Lk6F/q8efPmqeY5ed71vnE0ZMgQs33Dhg0e51HTN3To0FQ/M7Vr6km+07q3XJ0Pd9Oe1nFd8fS4/fv3d3p/ly5dLAUKFHDrsy5cuJBi21NPPWXu+UuXLrl1DIDqOWQJWuWhcufOneZ+ttdt+7tLq8hsrly5YkoHKlSoYKoR9H/EKikpSb777jvp3LmzlC5d2qlUo3379qkeW0se0vN5nrh27ZosXLhQ7rrrLqlfv36K193pzv+f//zH/ljbxOhxNP7QUg4bTZ9Wj/z999/2bXpeli5das6LtjezKVasmKni0VIWLSVx9/zpZ86fP9/kRR//+++/9kX3O3PmzA3Pkb7P3VImNXToUKfnDz/8sFl/8803HuXRdo7Wr19vSv48kZ58u7q3kvMk7Zl9XC2l0ns/+X6u2Ep7lZam6nnR6kO997XEEnAHQROyBHeDIX1dAwRbtYe7tCfauHHjpFSpUqYdib5fq0b0y1l/rNTx48fNflptkZwGEqlx1Sjdnc/zhKZNf3hq1Kgh6eUYyChtc6LtjpKfS92uVVeOn63Vf67OgQZE+qOm7aLcPX+6n54HW9s1x0XbwKhjx46JNyVPk7a1ypYtm71dkLt5VC+88IJs2rTJXFutLtVqSMcgMzXpybc7HR48Sbuvj5v8HtOgRzneT6kFaFo9pz1kNXjS9+l5efTRR839qW2pAHfQpglZgv5Qa+NmbceQFn1d2wjZ2u+kVsKiX8LJSxZmzZpl2ow0adLEfJ6+V9sc6Zd/RjiWKmXG56WXqx5XqfXCSqudS0bZ8q+N0vv27etyn1q1aokvZWSgUm23oyUoCxYsMCUx2jFB2+zocBnaYNqb+XZ1b3mDr46b3vupT58+Mm/ePHNe9G9GAyX9Gx8+fLgp2QLcRdCELEOrLbThpxb7u+rRs2rVKlMyMGrUKPs2/R+pq/FwtPeaI224q1/Ijj3DtPeU43v1f7b6Y7Jz584Ux9u+fbtHeXHn8zz58da0aWNcLeHIbPrZ+r9/V+dAe6VpiY2Wumi3dXfOnx5PSxY1sG3Tpo1kBk2TY+mK9vLSIMbW89LdPNroD7n2jtRFS4e0hOR///ufCZpSu6a+yrenaff3cZPTjhGffPKJaTjv2BBee81po/rbbrstw5+BrIPqOWQZ2mtLv6QffPBB0w7C0cmTJ017CQ0ctNeRYzWLVnc5llDpl62WAiT/H3Dy/+2+9tprTiVSuo+2LdG2Q47d2bVXlLbV8YQ7n6ds4+PcaCBE/YHStiXac/DXX3/N1JIhzYv2iNOBRR27uWuvJv2x0wBXr4u750/369q1q2nf4yoI1Gohbw05YDN9+vQU10LZSobczaNev+TVq9o7TEtJExMT07ym3si3K+6mPVCOm9zBgwfNWtv8ObYD1DZ4er59XeqI0EJJE7IM/dLULuA9evSQmjVrmgbKWjqgX9jvvfeeaRcxZ84cpxIDre564oknpEuXLqYoX9tgaNf+SpUqOTWqvfPOO+XDDz801WTaTXzt2rWmG7OO2eNIR3JesmSJqX7RUgRtgKo/sDqWzo2qDh25+3n16tUz6//+978mL9qlW0vcXNHu4VodpF3tdcgBbVeiAeJnn31mSue0gbKv6PAEy5YtMz+Uel4iIiJMqaAGCtrGx9Pzp8MexMfHmy7rAwcONOdIA2O9Znqe9LE3hxzQruwdO3Y0pRZ6LT766CPTmLl27doe5VHb1Gn18D333GPeq21tNL06xIStVDG1a6rBVEbzndHrEyjHdVSnTh3znyUdO0uvk7aL078dW1Wf4zUCbogOhMhq/vrrL0vPnj0tRYsWNd2N9c9Au7Fv3rzZ5f5Lly611KhRw3SNrly5suWjjz5KMeSAdqV/4IEHLAULFrTkypXL0r59e8u2bdvMcAXJu9evXLnSUq9ePXO8uLg4y4wZM1IcT9m2HT9+PEWaPPm8Z555xlKiRAl7XrUruatu5Uq78uvQA4UKFbJER0eb9Gn398TExFTPZ2rp1HTExMSk2F+78levXj3F9t9//93kQ/OTM2dOS8uWLS0//fRTiv3cPX9Hjx41aS9VqpQlMjLSXO/WrVtb3n77baf9vDHkgA57cM8995iu8vny5bMMGzbMZTf2G+VRz/Njjz1mqV27tjmWnj99/MYbb9zwmnqS77TurdTuDXeuT1rHTU1GjptaWpNbtGiRpUqVKubvvFatWuZ9EyZMsERERFgSEhLcTisQpqfgxqEVELq09ElHGdYGtI6DEQI3oj3btPRLq7487XEJIPhQPYcsT3vWaDXUmDFjTNWIVlMBAJAcQRMgYtot6QIAQGroPQcAAOAG2jQBAAC4gZImAAAANxA0AQAAuCHkG4LrVAY6W7hOL5CR+aAAAEBo0tGXdHBZHX1fZ0jIskGTBkzemL8IAACEtgMHDpihZ7Js0KQlTLYT4Y15jJB1nD171gTc3DvIUo5sFJl9u0i/b0SKMS8bstb3fe7/HzNk2aDJViWnARNBE9KDewdZyvlcItFhIrlz6c3v79QAmepGzXhoCA4AuK5ABZEBy6xrAFmrpAkA4IHoXCKlGnLKABcImgAA1505JLJ2ukiToSJ5SwTEmUlKSpIrV674OxkIAZGRkRIeHp7u9xM0AQCuu3BcZN10kVrdAiJoOn/+vBw8eNB0CQe80WZJe8flypUrXe8naAIABCQtYdKAKWfOnFKoUCHG2kOGaOB9/Phxc09VrFgxXSVOBE0AgICkVXL6Q6cBU44cOfydHIQAvZf27t1r7q30BE30ngMABDRmc0Cg3EsETQCA63IWEGnwH+sabtGSixkzZqT7bH311VcycuRItz4nNjY21dcnTJggCQkJbr23Tp06ZtoQbzpy5Ig0btzYTF/maxcvXpQePXpIhQoVpFKlSvL555/bXxs9erR88sknPvncMEuIt67TUT7z5s0rZ86cYXBLcO8AQUQDgD179ki5cuUke/bsZtu0ZTt88lkj21ZK93t/+OEHGTFihPz5558uX7969apERGS8NYwGPhrsnD59OtVSlFOnTrkMrG70Xm8YMmSINGnSRO6//37xtUmTJsnff/8ts2fPNvdIo0aNZOvWrVKgQAE5duyYNG3a1DxPXgXn6p7yJFagpAkBSb8Y3VkAeNnliyKH/7Su4eTSpUty3333SbVq1aR27drSrl07s33w4MGyfft2E5R07NjRbCtbtqw88cQT0rBhQ+nbt6/8888/0rJlS6lXr55Ur15dhg0bZi+R0R/+zp072z9n/PjxpgSlQYMG8tRTT5ljOdLX9Ti6zzfffGNPg2rWrJlJhwYOadEAyxZA6fHHjRtnAh4NJp599ln7fprubt26mXzUrFnTpMcVDUbmzp0rXbt2tW/r1auX1K9fX2rVqiV33HGHOZa36GfZ8qxpbtGihSxYsMA8L1y4sJQvX16WLl0q3ubXoGny5MnmptC5XjSTetPojedIT4ReXMfFdqIAAF727w6Rt5tb13CyZMkSE2hs2bJFNmzYIHPmzDHbtWqucuXKpqRJq9psTpw4IevXr5ePP/7YlP58/fXX8ttvv8nGjRtNyc+8efNSnOHFixfL/Pnz5Y8//pCff/5ZDh065PS6loRoEKLHef311+3VerbqwVWrVpl06G+qJzRfa9eulV9++UWmTp1q/1wN+IYOHWrS8scff8ivv/4qn332WYr36/s0eNGejjYvv/yy2V/zq8GcVh+6oudHAz1Xy/Tp012+Z//+/VKmTBn7cw38dJuNBoArVqwQb/Nr77mVK1eai6GBkxZfPvnkkyZy1xsyJibGvt/AgQNNUZyN40UBACAzaOmSVvloNVTz5s3l9ttvT3P/fv362Rsea6mSljytXr3a9AjUkqAaNWpI9+7dnd6jP/T33nuvfeLYAQMGSHx8vP11rVK6++677YHB7t27vZK3nj17mnXBggUlLi7OVGFpoKfpOXr0qNO4WduTFW4o7cZfpEgRp23arujDDz80pVC66LFd0RIpXbypaNGiJpYIqaBJo3ZHWkSp0bFG0LfeeqtTkKQnAAAAf9FgQn+Iv//+e1m+fLk8/vjjqbZjUo4DKL700ksmUNKSJw18Ro0alWqj7bR6e0VHR9u3aXsdHcvKGxzb9+hxtSDD1uR53bp1Tq+7or/TjvnR4PDVV181pVf6u64lcFoFmFpJk5ZuuaKFJlq4klzp0qVl3759UqxYMfNcS+5s1aVK0+KLYSoCqk2TFjuq/PnzpzihGqFqVD527FjTah4AgMykpSkasGi7pf/7v/8zQcWBAwdMw2Hb71dqtIG2/udfgw9t2+Oqiku1atXKVM9piY4ef+bMmW6nT0unbpQOT2jQp+2wpkyZYt92+PBhcx6S0ypDxxIoza+mRxtmX758Wd56661UP0dLmTT4dLW4CpiUlsbZqiS1VEwb4zu2C9MSQS0Z9LaAGdxSiy6198Ett9xigiPHIkOttyxevLipF9XiTb0wX3zxhcvjJCYmmsWxRTwCCw24gUwSP/nG+7Qc6/w8LJtIVG7rGk7++usv8x93DWa0JEZ7iWmwoI+1cbf+dmlplGO7JptHHnlE7rnnHrOf/p61adPG5dm98847TWmUtufR6jGtBkxrmAFHjz76qLRt29aU+mgjaE/bNbmihRZaKqZ5CwsLM01nNADSqUgcaXsmrZ7bvHmzyeNtt90mH330kWnrpYGT5jd5+6yMeOyxx6R///6mwbeWjGn7Llv1n14frVYcM2aMhOyQAw899JB8++23pkgv+cVwpMWirVu3ll27dpmTlZw2NJs4cWKK7Qw5EJpBU0a6Cd8Iw1UgSwZNASS17uGhTsdP0lIa/XnWQEh77b355psS6D777DNT4pNa4+3MbPqjAZsuyYXEkAPa9XLRokWmsVtaAZPSsRiUBk2u6P8CNNO2RYtOAQAIFn369JGbbrrJDG2gPcKeeeYZCQb33nuvVK1aNVMGt0yL/va/8MILPjm2X6vnNIp++OGHzdgKGp1q5HcjtkZ3tsZfyWkjOV0AAOlwbJvIZ31F7n1fpHAVTqEf2MYbCkbDhg3zdxLMWFq+4tegSRt4aZfEL7/80hRF2ga+0iIybfWuXSn1de3WqXWi2qZJx6TQnnVajwwA8LKrCSLHt1nXAAInaLLV0eoAlo5mzZplxreIiooy3Tp1gKwLFy5IqVKlzGijqY1ICgAA4Ct+r55LiwZJOgAmAACAvwXMkAMAAC/2igPgdQHRew4AECDylRXp/ql1DcAJJU0AgOtyxIpUSXtOtZAsZQvg8apSM3r0aDNqd2oT4eoQADowpc5R52vvvfeeGTlchxvQUc3feOMNiYyMTLGfzuih8+npBL/ZsmWT5557zgz6qXbu3CmDBg0yo4nreEp33HGHmV5F99PBK3UsJZ2j1p8oaQIAXHfuqMiqF61reJ2OHp4Zfv75Zzl58mSmBEx79uyRp59+WlatWmXGUNQJft9++22X++r0MzoskO733XffmcmPT5w4YR/lu0uXLvYpVHRUc9sctRpMaWDmzWli0oOgCQCyMi25cVqeE1kxybp23A5DJ6Bt2rSpmddMh77RIXNU2bJlnSbvrV+/vhl/0NZDfPjw4SaA0UlldRJaDR4cgw6dl+7KlStm0ek/GjZsaKZS6datmyl5UUeOHJH27dubQS91WhJXc8DZ6FQnOg2ZjQ7fo4ND66CZmvavv/7aa1f0888/N/PxaR50qpXBgwfLp59+6nLfuXPnmteVjs2o58Y2LpW+1xYU6Sjoei5sYzJqb3o9d5oPfyJoAgDADVpyo5PCTp48WTZs2GCCpGbNmrl17nbs2CE//vijmQrsgQcekNmzZ9tf08c6aa1WZ2l1lM7vpiVFevyaNWvah9nRwEuDqS1btsj7779v5ldLjQZsthk0lAZb69atkz/++MMEehq4Oc7T6jiFiwZrrpb27du7/CwdtVzniLXRAFK3ebqvDi+kU7Ho3Hy62EZGt9GgM608ZwbaNAEA4GYpk05AawuUtK1N/vz53Tp3vXv3trfxufnmm001nbbr0RKpDz74wF7ys3DhQlPaMn/+fPP88uXLJrBQGjDYSqhKlChhSndSo6VQOoGuY2mWBma6PSIiwgSAuq1KFedR33WgaccSs8z0xhtvSI8ePcx0aMeOHZOWLVtKgwYNzCTESkuy0ipdywwETQAAZPTHNCJCkpKS7M+1IbMjbbDtSEubdCDn8+fPS8GCBaVGjRr28Qtfe+01UxV1I1qdlZqcOXM6paF79+6mobat0bUGe8nTaCtpSq30rEiRIqYdUnKlS5c2M3jY7N2712xzRbfv27fPXu2m+9ryqhP9aomcKly4sJkNREvMbEGTpldnC/EnqucAANdFZBcpWNm6hhMtIdIeXtrgWWlPMS2xURUqVJD169ebx1q1tn379jTP3v3332+qombMmCH9+/e3b9fqv2nTppleZkrXmzdvNo+1HdPMmTPt7Zu++uqrVI+v7a0c06Dtomzzu3700Uf2dlLJ2UqaXC3fuQiYlM7UoWnRqdA06NM8aZCWWo8+fV1pSZcGRZpnFRcXZ2/4rbOAxMfH24NJtXXrVtMey58ImgAAzkMOVO9iXcNJvnz5TKNlbaitQUndunVlzZo15rVnn33WlJToj7oGNtWrV0/z7GmbHW2fpMGGVknZPPHEE6ZKStsj6Wc0btzYXl32yiuvmHZJ2hBc2/to1/7UaImSY5Cj79Vt2kZI2zWlVhKUHnFxcTJx4kS55ZZbTPBYqFAhefDBB81rhw8fNu2hbLSHnDbyLl++vGkjpUMJaEmb0nZa2kNOz6FWW7Zu3dop+NKAylZS5i9hlhvNZRLkdFwHnQBY64jz5Mnj7+RARKYtsxa/esPItpV8dk65dxCwfNmb7VqSyJULIpExItnC/TqOkVbHaGmElpBkz07Jlye02k9LxrQdljYsD3ZbtmwxgZitlM/b95S73/eUNAEArrtwXGTdG9Y1gpa2odJqPg0QQsGBAwfMMAr+RkNwBGQJEgAgY7R6K1S0T2W4g8xGSRMAAIAbKGkCgEDByNsuhXjTWwTRvUTQBAAISDoYpI5FdPz4cdMjK61xiQB3Aia9l/Q+cjWZsDsImhDy7ah82cMOCDm5iog0Gy0S5tBzzk/Cw8OlZMmSZhRoHQQRyCgNmPSe0nsrPQiaAACOvyoiYREB1QusYsWKZvJWIKO0hCm9AZMKnL8MAID/XTwpsuNbkUodRHK6N6+ar+mPXEZ+6ABvofccAOC6pMsiZw5Y1wCcEDQBAAC4gaAJAADADQRNAAAAbiBoAgBclz2PtRG4rgE4ofccAOC6yJwixWpzRgAXKGkCAFx35aLIkQ3WNQAnBE0AgOsSzlrHadI1ACcETQAAAG4gaAIAAHADQRMAAIAb6D2HkDdt2Q639hvZtpLP0wIEvPAokbylrGsATgiaAADX6SS9dXqlPCPxk298llqO5UwipFE9BwC4zmIRuXbVugbghKAJAHDd+aMiq/7PugbghKAJAADADQRNAAAAgR40TZ48WRo0aCC5c+eWwoULS+fOnWX79u1O+yQkJMjQoUOlQIECkitXLunatascPUqxMQAAyEJB08qVK01AtG7dOlm2bJlcuXJF2rVrJxcuXLDvM3LkSPn666/ls88+M/sfPnxY7r77bn8mGwAAZEF+HXJgyZIlTs9nz55tSpx+++03ufXWW+XMmTPy3nvvySeffCKtWrUy+8yaNUuqVq1qAq3GjRv7KeUAEKJiCok0HiISGePvlAABJ6DaNGmQpPLnz2/WGjxp6VObNm3s+1SpUkVKly4ta9eu9Vs6ASBkZQsXic5jXQMIzMEtr127JiNGjJBbbrlFatSoYbb9888/EhUVJbGxsU77FilSxLzmSmJiollszp5lpm4AcNul0yJ/x4vEtRTJ4fzdC2R1AVPSpG2bNm3aJHPmzMlw4/K8efPal1KlSnktjQAQ8q4miPy73boGEHhB07Bhw2TRokUSHx8vJUuWtG8vWrSoXL58WU6fPu20v/ae09dcGTt2rKnmsy0HDhzwefoBAEDo82v1nMVikYcfflgWLFggP/zwg5QrV87p9Xr16klkZKSsWLHCDDWgdEiC/fv3S5MmTVweMzo62izwz6S3AACEKo+DJh0OICYmxmtVctoz7ssvvzRjNdnaKWm1Wo4cOcx6wIABMmrUKNM4PE+ePCbI0oCJnnMAACCgq+e0EXb//v1l9erVGf7wN99801ShtWjRQooVK2Zf5s6da99n2rRpcuedd5qSJh2GQKvlvvjiiwx/NgDAhehcIuWaW9cAMhY0ffTRR3Ly5EkzblKlSpVkypQpZsDJ9FbPuVr69etn3yd79uwyffp085layqUBU2rtmQAAGRSVS6R0E+saQMaCJp3qZOHChXLo0CEZPHiwqV4rU6aMKQ3SgObq1aueHhIAEFC953bSew7wZu+5QoUKmbZGGzdulJdeekmWL18u99xzjxQvXlzGjRsnFy9eTO+hAQD+HKdp83zrGoB3es9pt//333/fTH2yb98+EzBpo+2DBw/K888/b6Y5Wbp0aXoPDwAAENxBk1bB6fxv3333nVSrVk2GDBkivXv3dhq1++abbzbzwwEAAGTZoOmBBx6Q7t27y5o1a6RBgwYu99Equv/+97/eSB8AIFjET3Zvv5ZjfZ0SIDCCpiNHjkjOnDnT3EfHWBo/fnxG0gUA8IdsESI5C1rXADLWEFwHoTx27FiK7SdOnJDwcGbFBoCgFlNQpMF/rGsAGQuadBwlVxITEyUqKsrTwwEAAAQFt8tfX331VbMOCwuTd999V3Lluj7wWVJSkvz4449SpUoV36QSAJA5zh8V+fNjkTq9RHIV4awD6QmadDoTW0nTjBkznKritISpbNmyZjsAIIhpbULSZesaQPqCpj179ph1y5YtzbAD+fLlc/etAAAAQc/j7hHx8fG+SQkAAECwB006XcozzzwjMTEx5nFadEoVIBhNW7bD6XnChfNmPf37XZI9xtqGb2TbSn5JG7LQGEYAgjto+uOPP+TKlSv2x6nRRuIAgCCWs4BI3X7WNQDPgybHKjmq5wAghIVHiuQu6u9UAKExTlNyZ8+elYULF8q2bdu8kyIAgP8knBHZudS6BpCxoKlbt27y+uuvm8eXLl2S+vXrm201a9aU+fPne3o4AEAguXJJ5PDv1jWAjAVNOohls2bNzOMFCxaYcZtOnz5tBr989tlnPT0cAABAaAZNZ86ckfz585vHS5Yska5du5oJfO+44w7ZuXOnL9IIAAAQfEFTqVKlZO3atXLhwgUTNLVr185sP3XqlGTPnt0XaQQAAAi+wS1HjBghvXr1MnPPlSlTRlq0aGGvttN2TQCAIBaVU6REA+saQMaCpiFDhkjDhg3lwIED0rZtW8mWzVpYFRcXR5smAAh20XlEKrT2dyqA0AialPaY08WRtmkCAAQ5naz3/HGRXIVEwqP8nRoguIOmpKQkmT17tqxYsUKOHTsm165dc3r9+++/92b6AACZ6eJJkT8/tI4KziCXQMaCpkceecQETVqyVKNGDaZOAQAAWYLHQdOcOXNk3rx5cvvtt/smRQAAAKEw5EBUVJRUqFDBN6kBAAAIlaDp0UcflVdeecWMBA4ACDFh2UQic1jXADJWPbd69WqJj4+Xb7/9VqpXry6RkZFOr3/xxReeHhJ+NG3ZDs4/gOtyFRa5+RHOCOCNoCk2Nla6dOni6dsAAACyVtA0a9Ys36QEAOB/F46LbJovUqOrSEwhf6cGCCjpqrS+evWqLF++XN566y05d+6c2Xb48GE5f/68t9MHAMhM15JEEk5b1wAyVtK0b98+ue2222T//v2SmJhoplLJnTu3PP/88+b5jBkzPD0kAABA6JU06eCWOoXKqVOnJEeOHPbt2s5JRwkHAAAIRR6XNK1atUp++uknM16To7Jly8qhQ4e8mTYAAIDgLWnSueZ0/rnkDh48aKrpAABBLEc+kZrdrGsAGQua2rVrJy+//LL9eVhYmGkAPn78eI+nVvnxxx/lrrvukuLFi5vjLFy40On1fv36me2Oi7anAgD4SES0SP446xpAxoKmF198UdasWSPVqlWThIQE6dmzp71qThuDe+LChQtSu3ZtmT59eqr7aJB05MgR+/Lpp596mmQAgLsSz4vsXWVdA8hYm6aSJUvKhg0bzMS9GzduNKVMAwYMkF69ejk1DHdHhw4dzJKW6OhoKVq0qKfJBACkx+XzIvvWiBSoKBKdi3MIZCRoMm+KiJDevXtLZvjhhx+kcOHCki9fPmnVqpU8++yzUqBAgUz5bAAAAI+Cpq+++krc1bFjR/EWrZq7++67pVy5crJ792558sknTcnU2rVrJTw83OV7dKwoXWzOnj3rtfQAAICsy62gqXPnzk7PtUG2xWJJsU256lmXXt27d7c/rlmzptSqVUvKly9vSp9at27t8j2TJ0+WiRMnei0NAAAvi598431ajuW0IzgbguswA7Zl6dKlUqdOHfn222/l9OnTZtHHdevWlSVLlvg0sXFxcVKwYEHZtWtXqvuMHTtWzpw5Y18OHDjg0zQBQEiJyC5SuLp1DSBjbZpGjBhhpkpp2rSpfVv79u0lZ86cMmjQINm6dav4io4FdeLECSlWrFiaDcd1AQCkQ45Ykap3ceoAbwRN2rYoNjY2xfa8efPK3r17PTqW9rxzLDXas2eP/Pnnn5I/f36zaDVb165dTe85/dzHH39cKlSoYII0AIAPXLsqknhOJDq3SLZ09RUCQpbH4zQ1aNBARo0aJUePHrVv08ePPfaYNGzY0KNj/frrr3LTTTeZRelx9fG4ceNMQ28d0kAblleqVMkMa1CvXj0zjQslSQDgIxf+Ffn5LesagBOP/xsxc+ZMMzlv6dKlpVSpUmabthuqWLFiihG9b6RFixYpGpQ7+u677zxNHgAEXqNmAFkzaNLqMS0BWrZsmWzbts1sq1q1qrRp08begw4IVdOW7XBrv5FtK/k8LQCAzJWuCmsNjnQOOl0AAACyAo/bNAEAAGRFdI0AAFyXu6hI8zGcEcAFSpoAAAC8HTRdvXpVPvjgA6fhBgAAIeTiCZHfP7CuAaQ/aIqIiJDBgwdLQkKCJ28DAASLpCsi5w5b1wAyVj2nA1jqqN0AAABZiccNwYcMGWJG7tYBLXWE7piYGKfXa9Wq5c30AQAABGfQ1L17d7MePny407hNOrK3rpOSkrybQgAAgGAMmnRSXQBAiMqeV6TKndY1gIwFTWXKlPH0LQCAYBGZQ6RIDX+nAgidwS13794tL7/8smzdutU8r1atmjzyyCNSvnx5b6cPAJCZLl8UOb5VpFBVkaicnHsgI73nvvvuOxMk/fzzz6bRty7r16+X6tWrm0l8AQBBLPGsyK5l1jWAjJU0jRkzRkaOHClTpkxJsf2JJ56Qtm3benpIAACA0Ctp0iq5AQMGpNjev39/2bJli7fSBQAAENwlTYUKFTKDW1asWNFpu24rXLiwN9OGDJq2bAfnEAAAfwVNAwcOlEGDBsnff/8tN998s9m2Zs0aef75582glwCAIBYeJZKvnHXtT/GT3duv5VhfpwRIf9D09NNPS+7cueXFF1+UsWOtN2vx4sVlwoQJTgNeAgCCUM78IrXu83cqgNAImnTUb20Irsu5c+fMNg2iAAAhwHLNOllveKRImMfNXoGQlqG/CA2WCJgAIIScPyayZpp1DcDzkqa6devKihUrJF++fHLTTTeZ0qbU/P777+4cEgAAIPSCpk6dOkl0dLR53LlzZ1+nCQAAIDiDpvHjx5t1UlKStGzZ0owCHhsb6+u0AQAABGebpvDwcGnXrp2cOnXKdykCAAAIhYbgNWrUMGM0AQBCUEwhkSbDrWsAGQuann32WRk9erQsWrRIjhw5ImfPnnVaAABBLFu4SFRO6xpAxsZpuv322826Y8eOTr3oLBaLea7tngAAQerSKZHdK0TKtxbJkc/fqQGCO2iKj4/3TUoAAP53NVHkxC6RMk39nRIg+IOm5s2b+yYlQBabLHlk20qZkhYAgB9HBF+1apX07t3bTNh76NAhs+3DDz+U1atXeylZAAAAQR40zZ8/X9q3by85cuQwo38nJiaa7WfOnJHnnnvOF2kEAAAIzt5zM2bMkHfeeUciIyPt22+55RamUAGAYBedWySulXUNIGNtmrZv3y633npriu158+aV06dPe3o4AEAgiYoRKdXQ36kAQqOkqWjRorJr164U27U9U1xcnLfSBQDwhysJIse3WdcAMhY0DRw4UB555BFZv369GZfp8OHD8vHHH5sBLx966CFPDwcACCQJp0W2LLSuAWSsem7MmDFy7do1ad26tVy8eNFU1UVHR5ug6eGHH/boWD/++KNMnTpVfvvtNzO6+IIFC6Rz585OA2bqZMHafkqr/rTd1JtvvikVK1b0NNkA4Ln4yZw1AOkvadLSpf/+979y8uRJ2bRpk6xbt06OHz8uzzzzjKeHkgsXLkjt2rVl+vTpLl9/4YUX5NVXXzUNz7VkKyYmxvTcS0ig2BgAAAR40NS/f385d+6cREVFSbVq1aRhw4aSK1cuEwDpa57o0KGD6Y3XpUuXFK9pKdPLL78sTz31lHTq1Elq1aolH3zwgakOXLhwoafJBgAAyNyg6f3335dLly6l2K7bNKjxlj179sg///wjbdq0ceqh16hRI1m7dq3XPgcA4CBbhEiuItY1ACdu/1WcPXvWlP7ooiVN2bNnt7+mk/R+8803UrhwYfEWDZhUkSJFnLbrc9trruhgm7YBN23pBgC4KaagSL0HOF1ARoKm2NhY055Jl0qVUs6ZpdsnTpwo/jZ58uSASAcAAMiiQVN8fLwpZWrVqpWZSiV//vz217R9U5kyZaR48eJeS5iOB6WOHj0qxYoVs2/X53Xq1En1fWPHjpVRo0Y5lTSVKlXKa+kCgJB27h+RPz4Uuel+kdzW72EAHgZNzZs3t7c1Kl26tClZ8qVy5cqZwGnFihX2IEkDIO1Fl9Z4UDr8gS4AgHSyJHHqAG80BN+6dausWbPG/lyHC9CgpmfPnnLq1CmPjnX+/Hn5888/zWILyPTx/v37TVA2YsQI07vuq6++kr/++kv69OljSrMcx3ICAAAIyKDpscceszeu1kBGq8Juv/12E/A4Vou549dff5WbbrrJLErfr4/HjRtnnj/++ONmwMxBgwZJgwYNTJC1ZMkSp0boAAAAmcHjPqUaHOn4TErbNt11113y3HPPye+//26CJ0+0aNHCtJNKjZY2TZo0ySwAAABBFTRpo2+dPkUtX77cVJkpbRhO9/7MMW3Zjkz6JABZTs4CIvUHiGSP9XdKgOAPmpo2bWqq0XQeuJ9//lnmzp1rtu/YsUNKlizpizQCADJLeKRITCHON+CNNk2vv/66REREyOeff24mzy1RooTZ/u2338ptt93m6eEAAIEk4YzI9m+sawAZK2nS4QYWLVqUYvu0adM8PRQAINBcuSTyz0aR4nVFsuf1d2qAgJKuyYV02hSdNFeHH1DVq1eXjh07Snh4uLfTBwAAEJxB065du0wvuUOHDknlypXtU5foqNuLFy+W8uXL+yKdAAAAwdWmafjw4SYwOnDggBlmQBcdjFJH8NbXAAAAQpHHJU0rV66UdevWOc09V6BAAZkyZYrpUQcACGJRMSKlGlvXADIWNOm8bufOnUuxXUfr1jGcAABBLDq3SFwLf6cCCI3quTvvvNNMa6IT5+po3rpoydPgwYNNY3AAQBC7mihyep91DSBjQdOrr75q2jQ1adLEzAGni1bLVahQQV555RVPDwcACCSXTols+NS6BpCx6rnY2Fj58ssvZefOnWbIAZ0frmrVqiZoAoCgED/Z3ykAkFXGaVIVK1a0B0oaOAEAEJABcMuxmZESZAEeV8+p9957T2rUqGGvntPH7777rvdTBwAAEKwlTePGjZOXXnpJHn74YdOuSa1du1ZGjhxpxmuaNGmSL9IJAMgMYdlEonJb1wAyFjTpJL3vvPOO9OjRw75Ne83VqlXLBFIETQAQxHIVFmky1N+pAAKSx/+VuHLlitSvXz/F9nr16snVq1e9lS4AAIDgDpruv/9+U9qU3Ntvvy29evXyVroAAP5w/pjI2unWNYCM957ThuBLly6Vxo0bm+c60KW2Z+rTp4+MGjXKvp+2fQIABBHLNZHL56xrABkLmjZt2iR169Y1j3fv3m3WBQsWNIu+ZsMwBAAAIEsHTfHx8b5JCQAAQACjTykAAIAbCJoAANflyCdSu4d1DcA706jAN6Yt28GpzSLcvdYj21byeVoAu4hokdgynBDABUqaAADXJZ4T+fsH6xqAE4ImAMB1ly+IHFhnXQNwQtAEAADgBoImAAAANxA0AQAAuIGgCQBwXWQOkaK1rGsAThhyAABwXfa8IpVv54wALhA0AQCuS7oiknBaJHusSHhkaJyZ+Mnu7ddyrK9TgiBH9RwA4LqLJ0R+fc+6BuCEoAkAAMANBE0AAABuIGgCAAAI9qBpwoQJEhYW5rRUqVLF38kCgNAWFu7vFAABKeB7z1WvXl2WL19ufx4REfBJBoDglbuoyK2P+TsVQEAK+AhEg6SiRYv6OxkAACCLC+jqObVz504pXry4xMXFSa9evWT//v3+ThIAhK4L/4r8Nsu6BhA8JU2NGjWS2bNnS+XKleXIkSMyceJEadasmWzatEly587t8j2JiYlmsTl79mwmphhA0AxkCNeuXRU5f9S6BhA8QVOHDh3sj2vVqmWCqDJlysi8efNkwIABLt8zefJkE1wBAABkqeo5R7GxsVKpUiXZtWtXqvuMHTtWzpw5Y18OHDiQqWkEAAChKaiCpvPnz8vu3bulWLFiqe4THR0tefLkcVoAAABCOmgaPXq0rFy5Uvbu3Ss//fSTdOnSRcLDw6VHjx7+ThoAhCadqLdaZ+saQPC0aTp48KAJkE6cOCGFChWSpk2byrp168xjAIAPRGYXKcQgwkDQBU1z5szxdxIAIGu5fEHk6GaRItVFomL8nRogoAR00ARAZNqyHTc8DSPbVuJUwTsSz4n8/b1IbGmCJiAZgqYA+uEDAACBK6AbggMAAAQKgiYAAAA3EDQBAK6LiBYpUMG6BuCENk0AgOty5BOpcQ9nBHCBoAlAcGAi3sxxLUnkaqK1pClbeCZ9KBAcqJ4DAFx34bjI2letawBOCJoAAADcQNAEAADgBto0AQDgbru5lmM5V1kYJU0AAABuoKQJAHBdrsIit4wUCY/krADJEDQBAK4Ly8bAlkAqqJ4DAFx38aTIxrnWNQAnBE0AgOuSLouc2mNdA3BC0AQAAOAG2jR5wbRlO7xxGAAAEMAoaQIAAHADJU1AFirtHNm2ks/TgiAXnUekQlvrGoATgiYAwHVROUVK1OOMAC4QNAEIjukrkDmuXBI5uVskf3mRyByc9fTeq0y3EpJo0wQAuC7hjMi2RdY1ACcETQAAAG6geg7IQtxpME5jcQBwjZImAAAANxA0AQCuC48UyV3cugbghOo5AMB1OQuI1O3DGQFcoKQJAADADQRNAIDrzv0jsnKKdQ3ACUETAACAG2jTBMB3GOkbyPjfB6OLBwxKmgAAANxASRMAjwfAVAyCCaSBUtaQRNAEALgupqBIwwdFonNzVoBkCJoApA//kw5N2SJEcuTzdyqAgBQUQdP06dNl6tSp8s8//0jt2rXltddek4YNG/o7WUCWtvbvEzfcp0lcgUxJC7zo0mmRvatEyjYTyRHLqQWCKWiaO3eujBo1SmbMmCGNGjWSl19+Wdq3by/bt2+XwoUL+zt5QMhpvP9tfycB/nQ1QeTYZpGSDbgOQLAFTS+99JIMHDhQHnjgAfNcg6fFixfLzJkzZcyYMf5OHgAA/q8KZ1iCTBHQQdPly5flt99+k7Fjx9q3ZcuWTdq0aSNr1671a9oAeKcKT1GNB2QQYz5lioAOmv79919JSkqSIkWKOG3X59u2bXP5nsTERLPYnDlzxqzPnj3rs3QmXDjvs2PDfxIunnda+0uDg7My9fMuSOZbvvlwpn5ew7L5M/XzgsrFRJFEi3WdLcHfqYG3LRrvvWPd+qiECluMYLFYgjdoSo/JkyfLxIkTU2wvVaqUX9KD4DepV3N/JwHwg7c467iBSSF3hs6dOyd58+YNzqCpYMGCEh4eLkePHnXars+LFi3q8j1alacNx22uXbsmJ0+elAIFCkhYWJhPolMNyA4cOCB58uSRUJaV8prV8kteQxfXNnRxbb1HS5g0YCpevHia+wV00BQVFSX16tWTFStWSOfOne1BkD4fNmyYy/dER0ebxVFsrO+7zeqPaqj/sGbFvGa1/JLX0MW1DV1cW+9Iq4QpKIImpaVGffv2lfr165uxmXTIgQsXLth70wEAAGSGgA+a7rvvPjl+/LiMGzfODG5Zp04dWbJkSYrG4QAAAFk6aFJaFZdadZy/aVXg+PHjU1QJhqKslNesll/yGrq4tqGLa5v5wiw36l8HAAAAycY5AAAAuDGCJgAAADcQNAEAALiBoCkZHQizV69eZtwLHd9pwIABcv582tNoJCQkyNChQ80Amrly5ZKuXbumGJBTB9ZMvsyZM8dpnx9++EHq1q1rGvdVqFBBZs+eLcGW1w0bNkiPHj3MoJA5cuSQqlWryiuvvJIin67Oh/aO9Kbp06dL2bJlJXv27NKoUSP5+eef09z/s88+kypVqpj9a9asKd98843T69r8T3txFitWzORN50DcuXNnhs9poOX1ypUr8sQTT5jtMTExZrC3Pn36yOHDzlOd6Oclv4ZTpkyRzODta9uvX78UebnttttC7toqV397ukydOtXv19aTvG7evNl8/9jSqsPRpOeY7nx/B0t+dUaMBg0aSO7cuaVw4cJmfMPt27c77dOiRYsU13bw4MESbHmdMGFCinzofe/za6sNwXHdbbfdZqldu7Zl3bp1llWrVlkqVKhg6dGjR5qnaPDgwZZSpUpZVqxYYfn1118tjRs3ttx8881O++ipnjVrluXIkSP25dKlS/bX//77b0vOnDkto0aNsmzZssXy2muvWcLDwy1LliwJqry+9957luHDh1t++OEHy+7duy0ffvihJUeOHCY/NvHx8eZ8bN++3el8JCUleS1vc+bMsURFRVlmzpxp2bx5s2XgwIGW2NhYy9GjR13uv2bNGnO+X3jhBXP+n3rqKUtkZKTlr7/+su8zZcoUS968eS0LFy60bNiwwdKxY0dLuXLlnK5jes5poOX19OnTljZt2ljmzp1r2bZtm2Xt2rWWhg0bWurVq+d0nDJlylgmTZrkdA3Pnz/v07z6Ir+qb9++5to55uXkyZNOxwmFa6sc86iLHjssLMz8vfrz2nqa159//tkyevRoy6effmopWrSoZdq0aek6pjvf38GS3/bt25vfmU2bNln+/PNPy+23324pXbq007Vr3ry5+SzHa3vmzJmgy+v48eMt1atXd8rH8ePHnfbxxbUlaHKgXzL6Y/7LL7/Yt3377bfmC+XQoUMuT6D+wOiX0meffWbftnXrVnMc/bGxn2gRy4IFC1K9EI8//ri5ARzdd9995o8g2PKa3JAhQywtW7ZMETSdOnXK4iv6Iz906FD7cw3Iihcvbpk8ebLL/bt162a54447nLY1atTI8uCDD5rH165dM3+8U6dOdTof0dHR5g87vec0EPOa2peY5m3fvn1OP6yuvsx8zRf51aCpU6dOqX5mKF9bzXerVq2ctvnj2nqaV3fSe6Njpvc7LVDzm9yxY8dMXlauXOkUND3yyCOWzNTQB3nVoEn/E5MaX11bquccrF271hS76+jjNloFky1bNlm/fr3LkrrffvvNVGfofjZaRFi6dGlzPEdaTKjz6enI5jNnznSaTVn3dTyGat++fYpjBEteHZ05c0by5085q7wOVKpVXW3btpU1a9aIt1y+fNmk1TGdmi99nlo6b3T+9+zZY6oPHffRIfe1mNm2T3rOaSDmNbVrqMXfyack0iobLfq+6aabTPXO1atXxZd8mV+tNtYqjcqVK8tDDz0kJ06ccDpGKF5brapYvHixqWpMLjOvbXry6o1jpvc7LRDzm9rfrUr+/fvxxx+b36IaNWqY+VovXrwowZjXnTt3muYDcXFxpup8//799td8dW2DYnDLzKI/ivql6SgiIsLccKm1t9HtOkde8h8THbHc8T2TJk2SVq1aSc6cOWXp0qUyZMgQ0x5i+PDh9uMkH+Vcn+uEjJcuXTJtaIIlr45++uknmTt3rvlittFAacaMGeYHKDExUd59911Tz64/PtqmK6P+/fdfSUpKcnk+t23blmreXO1vy5dtfaN9PD2ngZjX5LRdgLZx0rZqjnPw6b2r10vzp9dZv3yPHDkiL730kviKr/Kr7ZfuvvtuKVeunOzevVuefPJJ6dChg/ly1UnDQ/Xavv/++6b9i+bdUWZf2/Tk1RvHTM93WqDmNzmdp3XEiBFyyy23mODIpmfPnlKmTBkTbGzcuNH8bWu7py+++EKCKa+NGjUy7X71Pzl6b06cOFGaNWsmmzZtMve0r65tlgiaxowZI88//3ya+2zdutWnaXj66aftj/V/bjp/nv7vzRY0hVJebfTm7dSpkxlVu127dvbtepPrYnPzzTebH6pp06bJhx9+mClpg3v0f2rdunUzpaJvvvlminkhbWrVqmW+oB588EHTGDXYRlHv3r27/bE2ntb8lC9f3pQ+tW7dWkKVlnjr/9C1cW6oXtusSms29Dt49erVTtsHDRrkdK/rf2L1HtfvYL3ng0WHDh2c7lENojQYnDdvnsuSU2/JEkHTo48+anrHpEWL94oWLSrHjh1z2q5F0tprRl9zRbdr8ePp06edIlot9k7tPUov8DPPPGNKWvRLSPdN3qpfn+v/7D0pZQqUvG7ZssX8Ieof6FNPPXXDdGuVZfI/7vTSYmctHXB1PtPKW1r729a6Tb9kHPfRakbbPp6e00DMa/KAad++ffL99987lTKldk9rfvfu3esUFAdLfpP/jehn7dq1y9zHoXZt1apVq0wJg5YE34ivr2168uqNY6b3+zsQ8+tIpx1btGiR/Pjjj1KyZMkbXlul97ovgqaCPs6rjV6/SpUqmXz48tpmiTZNhQoVMnWZaS36P6kmTZqYE6x1oTb6Y6HFnLYbK7l69epJZGSkrFixwr5Nv4i0blWPl5o///xT8uXLZ/9fm+7reAy1bNmyNI8RqHnV7qItW7aUvn37yv/+9z+30q3nwzEYyQjNn6bVMZ2aL32e2vm80fnXahv9Q3PcR6tOtUrRtk96zmkg5tUxYNI2A8uXLzdtW9y5htpWIXk1VjDkN7mDBw+aNk22ezKUrq3Ne++9Z45fu3Ztv1/b9OTVG8dM7/d3IOZXaYmwBkwLFiww96d+b7lzbZW3vn8zK6/JaXMXLS2z5cNn1zbdTchDlHYrvummmyzr16+3rF692lKxYkWnbsUHDx60VK5c2bzu2K1Ru3V+//33pltjkyZNzGLz1VdfWd555x3T7Xfnzp2WN954wwwvMG7cuBRDDjz22GOmhf/06dMzZcgBb+dV81ioUCFL7969nbqCai8OG+0Jod329Vzo/tqTI1u2bJbly5d7tYur9mybPXu26fk0aNAg08X1n3/+Ma/ff//9ljFjxjh11Y6IiLD83//9nzn/2jPD1ZADeowvv/zSsnHjRtPryNWQA2mdU1/wdl4vX75shlMoWbKk6bbseB0TExPNPj/99JO5jvq6dlX/6KOPzHXv06ePT/Pqi/yeO3fOdG/WHjV79uwx92HdunXNtUtISAipa2ujXcz1++bNN99M8Zn+urae5lXvxT/++MMsxYoVM9dQH+v3irvHdOc7LZjy+9BDD5lhUXTIF8e/24sXL5rXd+3aZYaS0Hzqva7fZXFxcZZbb7016PL66KOPmnxqPvS+12FSChYs6PRb44trS9CUzIkTJ8wXYa5cuSx58uSxPPDAA+ZL1UYvkMaa2m3eRn80tVt9vnz5zBdRly5dzI3q2DW5Tp065pgxMTGmm+SMGTNSjEukx9T9dDwLvZF1vI1gy6t+Set7ki/abdTm+eeft5QvX96SPXt2S/78+S0tWrQwN7W36dhQ+gej51O7vOr4Oo7dbrWbuaN58+ZZKlWqZPbX4R8WL17s9LoOO/D0009bihQpYr4AWrdubcaa8uSc+oo382q77q4W273w22+/ma7s+gWt17Fq1aqW5557zinICJb86g9Ku3btTGCgAYbeqzqOjOMPa6hcW5u33nrLjJ+m3bKT8+e19SSvqd2nup+7x3TnOy2Y8pva363tt2T//v0mQNLvXf0O07HG9D/qvh6nyRd51SF5NKDS45UoUcI816DQ19c2TP/xRtEYAABAKMsSbZoAAAAyiqAJAADADQRNAAAAbiBoAgAAcANBEwAAgBsImgAAANxA0AQAAOAGgiYAAAA3EDQBAAC4gaAJQJbTokULGTFihL+TASDIEDQBCBqXL1+WQBOIaQLgGwRNAHzm2rVr8sILL0iFChUkOjpaSpcuLf/73//Ma2XLlpWXX37Zaf86derIhAkTnEqEhg0bZkqFChYsKO3bt5e3335bihcvbo7tqFOnTtK/f3/zeMmSJdK0aVOJjY2VAgUKyJ133im7d+82r/Xr109Wrlwpr7zyioSFhZll7969GUqTpmXy5MlSrlw5yZEjh9SuXVs+//zzG56f8ePHS82aNSUmJkaKFCkiDz30kFy5ciUdZxpAZiBoAuAzY8eOlSlTpsjTTz8tW7ZskU8++cQEB554//33JSoqStasWSMzZsyQe++9V06cOCHx8fH2fU6ePGkCpV69epnnFy5ckFGjRsmvv/4qK1askGzZskmXLl1McKPBUpMmTWTgwIFy5MgRs5QqVSpDadKA6YMPPjCPN2/eLCNHjpTevXub4Cw1Ole6Lm+99ZY5N7Nnz5b58+fLu+++61FaAGSeiEz8LABZyLlz50yA8vrrr0vfvn3NtvLly5sSIE9UrFjRlFY56tChgwnAWrdubZ5rqY6W+rRs2dI879q1q9P+M2fOlEKFCpngpEaNGibgyZkzpxQtWjRdeXNMU2Jiojz33HOyfPlyE4ypuLg4Wb16tQmImjdv7vIYWsI1adIk+/MyZcpImzZtZPv27elKEwDfo6QJgE9s3brVBBS2wCa96tWrl2KblihpqYweX3388cfSvXt3U6Kkdu7cKT169DDBS548eUy1m9q/f3+G0uIqTbt27ZKLFy9K27ZtJVeuXPZFS55sVYKu7Nu3T4YOHWqCuHz58pn3zJs3T0qWLOmVNALwPkqaAPiEtu1JiwY4Wj3lyFV7Hm3vk9xdd91l3rt48WJp0KCBrFq1SqZNm+b0upbcvPPOO/b2Txqc3KjRdnrSdP78ebPWtJQoUcJpP23H5crx48dNulu1aiUvvfSSeV9SUpLUr1/ftIcCEJgImgD4hFZhaeCkbYr+85//pHhdq8u0PZHN2bNnZc+ePW4dO3v27HL33XebEiYt6alcubLUrVvXvKbtnbSKSwOmZs2amW1aVeZIq+c0SPFGmqpVq2aCIy3FSq0qLrmvv/7afP6nn35qqumUVmNqgKYNzwEEJoImAD6hgc0TTzwhjz/+uAlSbrnlFlPCog2lBwwYYEpZtPGzlgppL7dx48ZJeHi428fXKjrtFafH00bXNlrVpT3mtJddsWLFTDAzZswYp/dqdd369etNrzmtFsufP78pZUpPmnLnzi2jR482jb+1REvbbJ05c8Y0EteqQVt7LkeaPg3IvvrqKxN0aRCljcm1xEkDNwCBiaAJgM9or7mIiAgTfBw+fNgEMYMHD7b3rNNSHA188ubNK88884zbJU1KAxwNdrRUqWfPnvbtGvzMmTNHhg8fbqrktBTq1VdfNUMF2GiQo8GMBiyXLl0yn6uBVHrTpPtpsKOBz99//20CLi35evLJJ13ur0GZBo7333+/KY3ToK9bt26mnROAwBVmSV6BDwAAgBToPQcAAOAGgiYAAAA3EDQBAAC4gaAJAADADQRNAAAAbiBoAgAAcANBEwAAgBsImgAAANxA0AQAAOAGgiYAAAA3EDQBAAC4gaAJAABAbuz/AUntd+bRONmAAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(6, 3.5))\n", "for (name, _, _, r2), c in zip(rows, ['C0', 'C1']):\n", " ax.hist(r2['samples'][:, 0], weights=r2['weights'], bins=50, range=(-0.05, 0.15),\n", " density=True, alpha=0.5, color=c, label=name)\n", "ax.axvline(0, c='k', lw=0.8)\n", "ax.axvline(0.08, c='C1', ls='--', lw=1)\n", "ax.set(xlabel='curvature $a$', ylabel='posterior density', title='Quadratic model: posterior on $a$')\n", "ax.legend(fontsize=8)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "de632c8c", "metadata": {}, "source": [ "---\n", "## Summary\n", "\n", "| Task | Code |\n", "|---|---|\n", "| Run | `NestedSampler(log_likelihood, prior, backend='nautilus').run(dlogz=...)` |\n", "| Results | `get_result()` → `log_evidence`, `samples`, `weights`, `equal_weight_samples`, ... |\n", "| Plot | `to_getdist()` → getdist `MCSamples` |\n", "| Compare models | $\\ln B_{21} = \\ln\\mathcal{Z}_2 - \\ln\\mathcal{Z}_1$ |\n", "| Other backend | `backend='dynesty'` (also gives `log_evidence_err`) |\n", "\n", "The likelihood must be **normalised** (include constants such as\n", "$-\\tfrac12 N \\ln 2\\pi\\sigma^2$) when evidences of models with different data\n", "normalisations are compared; constants common to all models cancel in the Bayes factor.\n", "\n", "Parallel evaluation: pass `pool=` (any object with a `map` method, e.g. a\n", "`multiprocessing.Pool`)." ] } ], "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.12.11" } }, "nbformat": 4, "nbformat_minor": 5 }