From 8289f499bfa09d444c454a297f2f835ca3de0f23 Mon Sep 17 00:00:00 2001 From: EnceladeCandy Date: Thu, 19 Oct 2023 22:47:03 -0400 Subject: [PATCH 1/9] New drp.py script (added choice for the seed, the normalization and bootstrapping) --- notebooks/test.ipynb | 196 +++++++++++++++++++++++++++++++++++++++++++ src/tarp/drp.py | 143 ++++++++++++++++++++++++++++++- 2 files changed, 337 insertions(+), 2 deletions(-) create mode 100644 notebooks/test.ipynb diff --git a/notebooks/test.ipynb b/notebooks/test.ipynb new file mode 100644 index 0000000..0ac69da --- /dev/null +++ b/notebooks/test.ipynb @@ -0,0 +1,196 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import sys \n", + "sys.path.append(\"../src/tarp/\")\n", + "from drp import old_get_drp_coverage, get_drp_coverage, bootstrapping" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Changes:\n", + "\n", + "The `get_drp_coverage` function now has a `seed` argument and a `norm` argument. \n", + "You can still use the old function but its normalization can cause the tarp test to fail in specific cases, so we'd recommend to use the new one. \n", + "\n", + "Just to make sure there's no mistake, let's compare the old function and the new one..." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Checking old - new\n", + "Expected coverage: 0.0\n", + "Credibility: 0.0\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAGHCAYAAACj5No9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACHLUlEQVR4nOzdd1QU59vG8e/u0jsKggoqWMBewN57iSj23qJGjbHGWBKjxmg0scbeNcVeY4yxJDH2GAsqscQuFhQFBaSzO+8f/OQNsYRVYBa4P+fsOe7M7O5Fu7139pnn0SiKoiCEEEIIIUQ2pFU7gBBCCCGEEG9KmlkhhBBCCJFtSTMrhBBCCCGyLWlmhRBCCCFEtiXNrBBCCCGEyLakmRVCCCGEENmWNLNCCCGEECLbkmZWCCGEEEJkW9LMCiGEEEKIbEuaWZGj/PHHH3To0IH8+fNjYWGBu7s77du35/jx42mOW7NmDRqNhlu3bv3nc9arV4969eplTmAhhMihntdZKysrbt++/cL+evXqUaZMGRWSiZxGmlmRY8yfP5+aNWty9+5dvvrqK3755RdmzpzJvXv3qFWrFgsWLFA7ohBC5DoJCQmMHz9e7RgiBzNTO4AQGeHo0aMMHz6cFi1asH37dszM/v9Xu3PnzrRp04Zhw4ZRsWJFatasqWJSIYTIXZo1a8a6desYNWoU5cuXVzuOyIHkzKzIEaZNm4ZGo2Hx4sVpGlkAMzMzFi1ahEajYfr06a98DkVR+OqrryhcuDBWVlZUqlSJn3/+ObOjCyFEjjZ69Gjy5s3LmDFjXnucoigsWrSIChUqYG1tjbOzM+3bt+fGjRupxyxcuBCtVktYWFjqtlmzZqHRaBg8eHDqNoPBgLOzMx9++GHGf0HC5EgzK7I9vV7PgQMH8Pf3x8PD46XHeHp64ufnx2+//YZer3/pMZ999hljxoyhcePG7Nixg0GDBtG/f3/+/vvvzIwvhBA5mr29PePHj2fv3r389ttvrzxuwIABDB8+nEaNGrFjxw4WLVrEhQsXqFGjBg8fPgSgUaNGKIrCr7/+mvq4X375BWtra/bv35+67dSpUzx9+pRGjRpl3hcmTIY0syLbe/z4MbGxsXh5eb32OC8vL2JjYwkPD39h39OnT/nyyy9p06YNK1asoFmzZvTv359Nmzbx4MGDzIouhBC5wsCBA/H29mbMmDEoivLC/j/++IPly5czffp0Zs6cSdOmTenSpQv79+8nKiqK2bNnA+Dj44OHhwe//PILAImJiRw+fJihQ4dy9epVQkJCgJQG19zcnDp16mTdFylUI82syDWeF1CNRvPCvuPHjxMfH0+3bt3SbK9RowaFCxfOknxCCJFTWVhYMGXKFE6dOsWmTZte2L9r1y40Gg3du3cnOTk59ebu7k758uX5/fffU49t2LBhajN77NgxYmNjGTlyJC4uLqlnZ3/55ReqV6+Ora1tlnx9Ql3SzIpsz8XFBRsbG27evPna427duoWNjQ158uR5Yd/zs7Xu7u4v7HvZNiGEEMbp3LkzlSpV4pNPPiEpKSnNvocPH6IoCm5ubpibm6e5/fHHHzx+/Dj12EaNGhESEsLVq1f55ZdfqFixIvny5aNBgwb88ssvxMXFcezYMRlikIvIbAYi29PpdNSvX589e/Zw9+7dl46bvXv3LqdPn6Z58+bodLoX9ufNmxfgpUMKHjx4QJEiRTI8txBC5CYajYYvv/ySxo0bs2zZsjT7XFxc0Gg0HD58GEtLyxce+89tDRs2BFLOvu7fv5/GjRunbh8/fjyHDh0iISFBmtlcRM7Mihxh3LhxKIrC+++//8IFXnq9nkGDBqEoCuPGjXvp46tVq4aVlRVr165Ns/3YsWMvnexbCCGE8Ro1akTjxo2ZPHkyz549S93esmVLFEXh3r17+Pv7v3ArW7Zs6rH58+enVKlSbN26ldOnT6c2s40bN+bRo0fMnj0bBwcHKleunOVfn1CHnJkVOULNmjWZO3cuw4cPp1atWnzwwQcUKlSIkJAQFi5cyIkTJ5g7dy41atR46eOdnZ0ZNWoUU6ZMoV+/fnTo0IE7d+4wadIkGWYghBAZ6Msvv8TPz4+wsDBKly4NpNTw9957jz59+nDq1Cnq1KmDra0toaGhHDlyhLJlyzJo0KDU52jYsCHz58/H2to6de5wLy8vvLy82LdvH61atXphmkaRc8mZWZFjDBkyhKNHj+Lh4cGHH35IgwYNGDlyJPnz5+fIkSMMGTLktY+fPHky06ZNSy2E8+fPZ8mSJfj4+GTRVyCEEDlfxYoV6dKlywvbly5dyoIFCzh06BCdO3fmnXfeYcKECcTExFClSpU0xz4fQlCrVi2srKxe2C5DDHIXjfKyOTKEEEIIIYTIBuTMrBBCCCGEyLakmRVCCCGEENmWNLNCCCGEECLbkmZWCCGEEEJkW9LMCiGEEEKIbEuaWSGEEEIIkW3luhmFDQYD9+/fx97eHo1Go3YcIYRKFEUhOjqaAgUKoNXK+3pTInVaCGFMjc51zez9+/fx9PRUO4YQwkTcuXMHDw8PtWOIf5A6LYR4Lj01Otc1s/b29kDKN8fBwUHlNEIItURFReHp6ZlaE4TpkDothDCmRue6Zvb5R1YODg5SJIUQ8jG2CZI6LYR4Lj01WgaKCSGEEEKIbEuaWSGEEEIIkW1JMyuEEEIIIbKtXDdmNj0URSE5ORm9Xq92FJGBzM3N0el0ascQQmQAqdM5j06nw8zMTMaxC6NJM/sviYmJhIaGEhsbq3YUkcE0Gg0eHh7Y2dmpHUUI8RakTudcNjY25M+fHwsLC7WjiGxEmtl/MBgM3Lx5E51OR4ECBbCwsJB3iDmEoig8evSIu3fvUrx4cTlDK0Q2JXU6Z1IUhcTERB49esTNmzcpXry4LGYi0k2a2X9ITEzEYDDg6emJjY2N2nFEBnN1deXWrVskJSVJMytENiV1OueytrbG3Nyc27dvk5iYiJWVldqRRDah6tueQ4cOERAQQIECBdBoNOzYseM/H3Pw4EH8/PywsrLC29ubJUuWZHgueTeYM8nZGyGMJ3VaZCX5uYo3oepvTUxMDOXLl2fBggXpOv7mzZu0aNGC2rVrExQUxMcff8zQoUPZunVrJicVQojcSeq0EMLUqTrMoHnz5jRv3jzdxy9ZsoRChQoxd+5cAEqWLMmpU6eYOXMm7dq1y6SUQojsYHdwKHltLajqnVftKDmK1GkhREYwGBQWHLhGj2qFcbbN2Av8stX5/OPHj9OkSZM025o2bcqpU6dISkp66WMSEhKIiopKc8ut6tWrx/Dhw03+OYUw1vp9xxn8/Um6rTjB2TtP1Y6Tq0mdfnNSo0VO1nbiKmbvv0KHpcdJ0hsy9LmzVTP74MED3Nzc0mxzc3MjOTmZx48fv/Qx06ZNw9HRMfXm6emZFVGzXO/evQkMDFQ7hhBZ7v7TOGaeiEbR6HBNDqNsQUe1I+VqUqdfTmq0yM12BN3jbJI7AE08wFyXse1ntmpm4cWLeBRFeen258aNG0dkZGTq7c6dO5meUQiRNWITk+n3zSnCY5Lwdbdn57i26LRyoZ/apE4LIZ47E/KE0VvPA/BudU9Gd6yb4a+RrZpZd3d3Hjx4kGZbWFgYZmZm5M378nFylpaWODg4pLm9iZiYmFfe4uPj031sXFxcuo59WzExMfTs2RM7Ozvy58/PrFmzXjgmMTGR0aNHU7BgQWxtbalatSq///576v7w8HC6dOmCh4cHNjY2lC1blvXr1xudZefOnfj7+2NlZYWLiwtt27ZN3ff999/j7++Pvb097u7udO3albCwsNT9T548oVu3bri6umJtbU3x4sVZvXp16v579+7RqVMnnJ2dyZs3L61bt+bWrVtGZxTZT1JSMp1m7+ZiaBQudhas6OWPq/Ob/X2LjJNb6vTbkhotcoMTwVfos/I4ickGGpdyY3xA2Ux5nWzVzFavXp39+/en2bZv3z78/f0xNzfP1Ne2s7N75e3fFzXky5fvlcf++0KKIkWKvPS4t/XRRx9x4MABtm/fzr59+/j99985ffp0mmP69OnD0aNH2bBhA+fPn6dDhw40a9aMq1evAhAfH4+fnx+7du3ir7/+4r333qNHjx6cOHEi3Tl++ukn2rZtyzvvvENQUBC//vor/v7+qfsTExP5/PPPOXfuHDt27ODmzZv07t07df+nn37KxYsX+fnnn7l06RKLFy/GxcUFgNjYWOrXr4+dnR2HDh3iyJEj2NnZ0axZMxITE9/iuyeyg1bjlxH8VIdW0bO0hz8ezjLnqCnILXX6bUmNlhqd0z14/IQuiw4SmaBQ0BbmdqqANrM+OVNUFB0drQQFBSlBQUEKoMyePVsJCgpSbt++rSiKoowdO1bp0aNH6vE3btxQbGxslBEjRigXL15UVq5cqZibmytbtmxJ92tGRkYqgBIZGfnCvri4OOXixYtKXFzcC/uAV95atGiR5lgbG5tXHlu3bt00x7q4uLz0OGP16tVLad26taIoKd9XCwsLZcOGDan7w8PDFWtra2XYsGGKoijKtWvXFI1Go9y7dy/N8zRs2FAZN27cK1+nRYsWyocffph6v27duqnP+TLVq1dXunXrlu6v488//1QAJTo6WlEURQkICFD69Onz0mNXrlyp+Pj4KAaDIXVbQkKCYm1trezdu/eF41/38xXZy/A5a5XCY3YphcfsUkYt2PxGz/G6WiD+n9Tpl9dpY0mNTvG6Gq0oUqdzisTEJKXMoPlK4TG7lEJD1yknL143+jmMqdGqTs116tQp6tevn3p/5MiRAPTq1Ys1a9YQGhpKSEhI6n4vLy92797NiBEjWLhwIQUKFGDevHlZMt3Ls2fPXrnv36tJ/fMjmH/794TQmfFxy/Xr10lMTKR69eqp2/LkyYOPj0/q/TNnzqAoCiVKlEjz2ISEhNSPAvV6PdOnT2fjxo3cu3ePhIQEEhISsLW1TXeWs2fP0r9//1fuDwoKYtKkSZw9e5aIiAgMhpQrHENCQihVqhSDBg2iXbt2nDlzhiZNmhAYGEiNGjUAOH36NNeuXcPe3j7Nc8bHx3P9+vV0ZxTZy8ofDrDtrg0aMyiru8+Mwa/+/RJvT+p0xtdpqdFSo3O6lh8vJdrBCyU5ianNCuFf0jtTX0/VZrZevXqpFwa8zJo1a17YVrduXc6cOZOJqV7OmOKQWcem1+u+p88ZDAZ0Oh2nT59+ocg//wht1qxZzJkzh7lz51K2bFlsbW0ZPny4UR8PWVtbv3JfTEwMTZo0oUmTJnz//fe4uroSEhJC06ZNU1+jefPm3L59m59++olffvmFhg0bMnjwYGbOnInBYMDPz4+1a9e+8Nyurq7pziiyjwt/X+Wz/XfQ2uXFOSaE7XPeUztSjid1OuPrtNRoqdE52SfzvuVvXREAOhfV0715zUx/zWw1ZlakT7FixTA3N+ePP/5I3fbkyROuXLmSer9ixYro9XrCwsIoVqxYmpu7e8r0GYcPH6Z169Z0796d8uXL4+3tnTpWK73KlSvHr7/++tJ9ly9f5vHjx0yfPp3atWvj6+v70rMlrq6u9O7dm++//565c+eybNkyACpVqsTVq1fJly/fC1+Do6NMz5TTREZGEjh6Hlq7vGjjnrLv826Ymen++4FCmBip0VKjc6rfDhxg9dlIAHx0YXw5KGsWSpFmNgeys7Ojb9++fPTRR/z666/89ddf9O7dO81HZyVKlKBbt2707NmTbdu2cfPmTU6ePMmXX37J7t27gZSCu3//fo4dO8alS5cYMGDAC1cp/5eJEyeyfv16Jk6cyKVLlwgODuarr74CoFChQlhYWDB//nxu3LjBzp07+fzzz9M8fsKECfzwww9cu3aNCxcusGvXLkqWLAlAt27dcHFxoXXr1hw+fJibN29y8OBBhg0bxt27d9/mWyhM0K5fDpFQNGVKlwmB5XF1kpkLRPYkNVpqdE41Y/NBLPJ5ozMksm5c1yx7XWlmc6gZM2ZQp04dWrVqRaNGjahVqxZ+fn5pjlm9ejU9e/bkww8/xMfHh1atWnHixInUCcs//fRTKlWqRNOmTalXrx7u7u5GT/pdr149Nm/ezM6dO6lQoQINGjRIvdLW1dWVNWvWsHnzZkqVKsX06dOZOXNmmsdbWFgwbtw4ypUrR506ddDpdGzYsAEAGxsbDh06RKFChWjbti0lS5bk3XffJS4u7o2n9hGm6y9NEbQW1hTLY07PuqXVjiPEW5EaLTU6p4lL1PPEoxYAI5qUJK+dZZa9tkZJz+CdHCQqKgpHR0ciIyNf+GOKj4/n5s2beHl5YWVlpVJCkVnk55s9JScncyM8jmZzD2FQYNOA6lTxyvPWz/u6WiDUJXU695Kfb/aj1+vRarUs+O0as/ZfoaCTNb9+WBcr87cbBmZMjVb1AjAhhHid/fv3M2z4CEq8Nw+DAk1Lu2VIIyuEECJjDB06lLBYAxc8WwMwprnvWzeyxpJhBkIIk3T58mU6dOjAfaeynA2Nw8pcy7jmJdWOJYQQ4n8WLlzIoqXLOKL4EJuop3IRZwLK5c/yHNLMCiFMTnh4OAEBASQ4FcGpVjcApgSWpYhLxk9lJ4QQwnj79u1j2LBhONfvh6V7cZxtzPm6c0U0mkxa5es1pJkVQpiUxMRE2rdvz80HEbgFjgGNhk7+nrT381A7mhBCCODSpUt06NABy+I1cPBrCcDsThUo4PTqeYszkzSzQgiToSgKH3zwAQeP/oFb+wlg5YCvuz2ftZbZC4QQwhSEh4fTsmVL4m3dyddyBAAf1C9GfZ98qmWSZlYIYTKWL1/O8lWrcW37CeZuxXCyMWdRt0pZfjGBEEKIFymKQufOnbkTmUT+Tp+j6CyoXdyF4Y2Kq5pLZjMQQpiMgFatmfNnNHEuvthY6FjTpwrernZqxxJCCAFoNBr6Dx/D1b0PwdKWCp5OLOnuh5lO3XOjcmZWCGES4pP0TPntHnEuvljotCzv6U8FTye1YwkhhPifqw+jmfeXFmycKOFmx+relbG1VP+8qPoJhBC5WlhYGNv3/MbuKA/O3Y3ETKthXpeK1CzmonY0IYQQpMxcEKZz4ctDj4hOSKZwXhu+61sVZ1sLtaMB0swKIVSUkJBAy27vcb94IGYOkTjZmLOkux/VvPOqHU0IIQRw/vx5un66CNu6fdFotVQu4syS7n5Zulztf5FhBjlE79690Wg0aDQazM3NcXNzo3HjxqxatQqDwZB6XJEiRVKPs7a2xtfXlxkzZvDPVY1v3bqVeoxGo8HZ2Zk6depw8OBBNb40kUMpikL3gSMI9WmPmYMrno7m7Hi/pjSyIseSOi2ym4cPHxIw7Avs6vdHo9XStmIBvu9X1aQaWZBmNkdp1qwZoaGh3Lp1i59//pn69eszbNgwWrZsSXJycupxkydPJjQ0lEuXLjFq1Cg+/vhjli1b9sLz/fLLL4SGhnLw4EEcHBxo0aIFN2/ezMovSeRgn34xk+OWfujsnPG007BrWH1ZFEHkeFKnRXYRHx9P035joXJXAHpWKcCsjhWwNDO92WVkmMF/UBSFuCS9Kq9tba4zaiUNS0tL3N3dAShYsCCVKlWiWrVqNGzYkDVr1tCvXz8A7O3tU4/r168fixcvZt++fQwYMCDN8+XNmxd3d3fc3d1ZunQpHh4eLz1OCGN9t3Ebq27aYuHijrNZEtuGNcfRxlztWCKbUqtOG1ujQeq0yB4URaHd+5/wxKc1Gq2OZiUc+KxNBVVW90oPaWb/Q1ySnlIT9qry2hcnN8XG4u1+RA0aNKB8+fJs27YttUg+pygKBw8e5NKlSxQv/vo54mxsbABISkp6qzxCHDh+inH7H2DhUhgrJZ6dI5vjam9aH1mJ7EWtOp0RNRqkTgvTM/CzeQQ710RrZo6/m46FvWuZbCMLMswgV/D19eXWrVup98eMGYOdnR2WlpbUr18fRVEYOnToKx8fExPDuHHj0Ol01K1bNwsSi5wqNDKOkbvuYOZSGF1SDDuGNcQzj43asYRQndRpYSp+Pn+PvbFeaM0tKWYTz9ohjdFpTbeRBTkz+5+szXVcnNxUtdfOCIqipHlH9dFHH9G7d28ePXrEJ598QoMGDahRo8YLj6tRowZarZbY2Fjy58/PmjVrKFu2bIZkErnPg8h4Oi49zhO9Bc4WCt8OqINvASe1Y4kcQK06nVE1GqROC9Pwc3AoH2w4B1od5ZyS2TKqDRZmpn/eU5rZ/6DRaDLkYyQ1Xbp0CS8vr9T7Li4uFCtWjGLFirF161aKFStGtWrVaNSoUZrHbdy4kVKlSuHk5ETevHKFuXhzBoOBST8EcycijsJ5bVjbryoeznJGVmQMqdNSp8Xbux8eydhtwegNCm0rFeSrduVUX9krvbJHSvHGfvvtN4KDg2nXrt1L9zs7OzNkyBBGjRqVZtoXAE9PT4oWLSoFUry19yfMZM/FMDTAku5+0sgK8Q9Sp4XaYmJiaPDBV0TGJeHjZseM9uWzTSML0szmKAkJCTx48IB79+5x5swZvvjiC1q3bk3Lli3p2bPnKx83ePBg/v77b7Zu3ZqFaUVusXbtOn64nfJxbCXnBErmd1A5kRDqkTotTI3BYKBT38HEeVQGYEDVfCY/RvbfpJnNQfbs2UP+/PkpUqQIzZo148CBA8ybN48ffvgBne7VY7tcXV3p0aMHkyZNSjNxtxBv6/jx4wyevgLLgr7oFD2LBrVQO5IQqpI6LUzNp59+yon4/Gh05pR1NaNtjZJqRzKaRvn3ZxY5XFRUFI6OjkRGRuLgkPYMUXx8PDdv3sTLywsrKyuVEorMIj/frHX79m2qVKuBWevJmDu5M6xBMUY08VE7VqrX1QKhLqnTuZf8fLPW999/T7+xU3HvPhMNCnuG18XH3V7tWIBxNVrOzAohMlx0dDStWrUi0bcp5k7uuNpZMKBeUbVjCSGE+J9jx47Rb8D75Gk2BIDOVQqZTCNrLGlmhRAZbtiwYVyNtcKxasoFLZ8Hlsn2V5sLIURO8ezZM9q2bYtd/f5YuBQmn70lH5rQJ2fGkv9dhBAZ7r0R4zjgfA4FeLemF83K5Fc7khBCiP+xs7Oj56QlbLpljlYD87tUxMUu+67EKGdmhRAZKiFZz7SDD1HMrang6cTY5r5qRxJCCPEPF+5HsuNuSvM6qqkPVb2z99Rucmb2JXLZNXG5hvxcM9fhw4d58uQJ16x9OXc3EicbcxZ2q5QtVo8R2Y/8PedM8nPNXPPmzaNlq9Z8uOM2ickG6vu4MrBO9r+eQZrZfzA3NwcgNjYWa2trldOIjJaYmAjw2ulvxJu5fv06bdq0IVprh0ffBYCGKYFlKOgkf0ciY0mdztliY2OB//85i4yzcuVKhg0bxlc/ncesYhucbcyZ2aE82mw2p+zLSDP7DzqdDicnJ8LCwgCwsbFJs1a2yL4MBgOPHj3CxsYGMzP5tc9IkZGRBAQEEB7xhKIDppKsaGhU0o13yso4WZHxpE7nTIqiEBsbS1hYGE5OTnLSIYP9/vvvDBw4ELM8HlhUbI0BmBBQirzZeJzsP8n/6v/i7u4OkFooRc6h1WopVKiQ/MeXgZKTk+nUqROXLl3Co0EPkh09sLc0Y0pgGfk+i0wjdTrncnJySv35ioxx7do12rVrR3KynjI9JhONlno+rgRWKKh2tAwjzey/aDQa8ufPT758+UhKSlI7jshAFhYWaLUyfjMjjRw5kr179+JY3B+bap1J0CuMbeGLu6NMdi4yj9TpnMnc3FzOyGawJ0+e0LJlSyIinlCy81iirfJhY6HLcSccpJl9BZ1OJ39UQrzG4sWLmT9/PnblmpCnxRAS9Aq1irnQpXIhtaOJXELqtBCvlpSURMeOHbly/RaenScRW9gPgI9blMTD2UbldBlLmlkhxBv5+8oVnOr1wbFqOwwKtCyXP8dcTCCEENnd06dPeRgZS/7uX6F1K4q5TsPUwLJ0rOypdrQMJ82sEOKNVOwwjB2WFwAY3qg4wxoWz1EfWwkhRHaWJ68LPn2+5OTtSJxtzFnS3S/bzyf7KtLMCiHSLTIyEhsbGy48iGHKTxcBGP9OSfrV9lY5mRBCCIDw8HDy5s3LnP1XOHk7ElsLHZsHVqdYPnu1o2UaaWaFEOmSmJhIYGAgWNgQX3cESXqF5mXc6VvLS+1oQgghgEuXLlGjRg06DZ/MnriUkwzT25XL0Y0syHK2Qoh0UBSFwYMH8/vvv/O3UxVCoxIokteGL9uXk6EFQghhAh4/fkzLli2J1uvYF5UyvVmPaoUJKF9A5WSZT5pZIcR/mjt3LitWrMCubCMsvPyxMNOysFslHKxklR4hhFBbYmIi7dq148aNG3gEfoTB3IayBR0Z37Kk2tGyhDSzQojX+umnn/jwww/R2jhR4J0hAIxsXILSBRxVTiaEEEJRFAYNGsShQ4fIW7EpFCiDhU7LnE4VsDTLHVPXSTMrhHil4OBgOnfujKIo+L83nQRFR+kCDvSTcbJCCGESZs2axapVqzCzdcT9fycchjQoRrF8dionyzrSzAohXspgMNCtWzeePXtGlcB3eWjpgU6r4ct25TDTSekQQgi1XbhwgdGjRwPQ6KPFPEsCHzd7BtQtqnKyrCWzGQghXkqr1bJu3ToGjf+KZxU7Qmwy/Wp7UaagDC8QQghTULp0aRYvXsKOv59xKc4BjQa+bF8OC7PcdcJB9a920aJFeHl5YWVlhZ+fH4cPH37t8WvXrqV8+fLY2NiQP39++vTpQ3h4eBalFSJ3uZ7szMMyXXgSm0zpAg4Mb1hC7UhCBVKnhTBNCcl6rjhX5ZKFDwBD6hejgqeTuqFUoGozu3HjRoYPH84nn3xCUFAQtWvXpnnz5oSEhLz0+CNHjtCzZ0/69u3LhQsX2Lx5MydPnqRfv35ZnFyInGvevHkcPHSI2fv+ZtiGsyQmG2hUMh8bB1TH2iJ3XEwg/p/UaSFMS3x8PMOHD+dqSCjdlp9g65m76LQaJrcuzcgmPmrHU4VGURRFrRevWrUqlSpVYvHixanbSpYsSWBgINOmTXvh+JkzZ7J48WKuX7+eum3+/Pl89dVX3LlzJ12vGRUVhaOjI5GRkTg4OLz9FyFEDrJt2zbad+qKa8CHWJeoAcCAOt6MbuaLTpuz5pOVWpA+UqeFMB2KotC9e3c27zuKR9cpGKydsbc0Y2G3StQp4ap2vAxlTB1Q7cxsYmIip0+fpkmTJmm2N2nShGPHjr30MTVq1ODu3bvs3r0bRVF4+PAhW7Zs4Z133nnl6yQkJBAVFZXmJoR40ZkzZ+g1cChuXadjXaIG5joNX7Uvx7gWJXNcIyvSR+q0EKZl6tSpbPvjb9y7z8Bg7UyhPDZsH1wjxzWyxlKtmX38+DF6vR43N7c0293c3Hjw4MFLH1OjRg3Wrl1Lp06dsLCwwN3dHScnJ+bPn//K15k2bRqOjo6pN09Pzwz9OoTICe7fv0+rwDbYB4zDMn9xnG3MWduvGh395e8lN5M6LYTp2Lx5M58vWEO+tp+itbShqlcefhhcM8cvVZseql8A9u+lMBVFeeXymBcvXmTo0KFMmDCB06dPs2fPHm7evMnAgQNf+fzjxo0jMjIy9Zbej7mEyC1iY2Np3bo1zwrVxMK1CM425vwwuBZVvPKoHU2YCKnTQqjr1KlT9OrzLnmbD0OjM6NxKTe+61sVZ1sLtaOZBNWm5nJxcUGn073w7j4sLOyFswDPTZs2jZo1a/LRRx8BUK5cOWxtbalduzZTpkwhf/78LzzG0tISS0vLjP8ChMgBDAYDvXv35tytMPL3Hg/A5NZlKJTXRuVkwhRInRZCfXfv3qVVq1aYl22OhZs3zjbmTG9bNtdNv/U6qn0nLCws8PPzY//+/Wm279+/nxo1arz0MbGxsWi1aSPrdClXV6t4HZsQ2VZycjJanRku/3u336hkPlqWe7HZELmT1Gkh1BcTE4OVayGca3UF4NOWpchrJ2/+/knVtn7kyJGsWLGCVatWcenSJUaMGEFISEjqx1Hjxo2jZ8+eqccHBASwbds2Fi9ezI0bNzh69ChDhw6lSpUqFChQQK0vQ4hsy8LCgqYfTMWigA92lmZ8HljmlR8fi9xJ6rQQ6ipevAQV+n8FOnPqlHClTcWCakcyOaquANapUyfCw8OZPHkyoaGhlClTht27d1O4cGEAQkND08xl2Lt3b6Kjo1mwYAEffvghTk5ONGjQgC+//FKtL0GIbOn27dt4enpyMzyWL/dcBmBMc1/yO1qrnEyYGqnTQqjj5s2beHl5seroTc7cjcbaXMdUOeHwUqrOM6sGmb9Q5Ha3b9+mSpUq1KpTl8TaHxB8P5paxVz49t0qaHPRFFxSC0yX/GxEbvfdd9/x7rvvMnHOUr59kJ/EZANT25ShW9XCakfLMsbUAVXPzAohslZ0dDQBAQGEhYXxl+JBwv1oHKzMmNGhXK5qZIUQwlQdO3aMfv36kazAlru2JGKggW8+ulYppHY0kyWXwgmRS+j1erp160ZwcDDupauRXLwhAFPalJXhBUIIYQJu3bpFYGAgiYmJVO7zGU+xI4+tBdPblZXhBa8hzawQucTYsWP58ccfsbS0xLfbRPQKtCpfgFbl5aIcIYRQW1RUFAEBATx69IiyNRoR7lIegC/alCWfvZXK6UybNLNC5AKrVq1i5syZAHw45ztuRuqxtdAxMaCUysmEEELo9Xq6du3KX3/9hbu7O5XfnYRegUYl89GsjLva8UyeNLNC5HAREREMHz4cgPGfTuBEXD4A+tb2lrkKhRDCBGzatImffvoJKysrvv5mK79efYpGAx828VE7WrYgF4AJkcPlyZOHPXv28P3331O+VT++23IeR2tz+tX2UjuaEEIIoHPnzty6dQtvb2/23jcHIKBcAUrml9k80kOaWSFygRo1auBfpRoNZv0OwKB6RXGwMlc3lBBCCAA0Gg3jxo3j9O0Ixiw+jk6rYUTjEmrHyjZkmIEQOVBycjJ9+/bl3LlzqdvWnrjN3SdxuNpb0qt6EfXCCSGE4Nq1a/Tu3Ztnz54BYDAofPnz3wB08PPAy8VWzXjZipyZFSIHGjlyJKtWrWLXrl3cvHmTsFiFr/akFMmhDYtjbaFTOaEQQuReT548oWXLlvz999/odDpWrlzJt8dv8eetCKzMtQxpWFztiNmKnJkVIodZvHgx8+fPB2DJkiVYWFoxYtNZ4pL0VPfOSzeZeFsIIVSTlJREx44d+fvvv/H09GTq1KlcfRjNtJ9Tlhb/pEVJCjrJ3N/GeONm9tq1a+zdu5e4uDgActmquEKYpP379zNkyBAAvvjiC9q0acPi368TFPIUe0szZnYsLyt95SJSp4UwPcOHD+eXX37B1taWnTt3ksclHyM2nSUh2UCdEq50r5Z7lqzNKEY3s+Hh4TRq1IgSJUrQokULQkNDAejXrx8ffvhhhgcUQqTP5cuX6dChA3q9np49ezJ27Fj+vBnB179eBWByYGl5t59LSJ0WwjQtWLCARYsWodFoWLt2LeXLl+fLPZf5614UTjbmzGhfTlb6egNGN7MjRozAzMyMkJAQbGxsUrd36tSJPXv2ZGg4IUT6hIeHExAQQGRkJDVr1mTZsmX8cPY+3VeeINmg8E65/ARWKKh2TJFFpE4LYXr27t3LsGHDAJg+fTrN3mnJh5vPsfLITQCmBpbFzUFW+noTRl8Atm/fPvbu3YuHh0ea7cWLF+f27dsZFkwIkX5mZmZ4e3uTnJzMlq1bmf/7LRYcuAZA41Ju8m4/l5E6LYTpcXZ2xtXVlWbNmtFn0FC6LT/BqdtP0Gk1TAooxTvl8qsdMdsyupmNiYlJ807/ucePH2NpKasJCaEGR0dHfvrpJ+7fv8+2S9GpjezAukUZ3dRHxsnmMlKnhTA9VapU4fTp0zg556HrqlOcuxuJvZUZi7pVonZxV7XjZWtGDzOoU6cO3377bep9jUaDwWBgxowZ1K9fP0PDCSFeLygoKPWiHjMzM5Ks8zD3l5QxspMCSjG2ua80srmQ1GkhTENiYiIXLlxIvV+wYEHWnQpNbWS3v19DGtkMYPSZ2RkzZlCvXj1OnTpFYmIio0eP5sKFC0RERHD06NHMyCiEeIldu3bRqlUr3nvvPRYuXIhGo2XctmAS/3dFbK8aRdSOKFQidVoI9SmKwsCBA9mwYQPr1q0jMDCQ2+ExzNqfMuf3+HdKUiyfvcopcwajz8yWKlWK8+fPU6VKFRo3bkxMTAxt27YlKCiIokWLZkZGIcS/BAcH06VLFxRFwWAwoNVqWX8yhD9vRmBtrmNqYBkZI5uLSZ0WQn2zZs1i9erVJCQkYGVlhaIojNsWTHySgereeeno76l2xBzjjVYAc3d357PPPsvoLEKIdAgLCyMgIIBnz55Rr149Fi5cyL2ncUzfnTLh9qimPnjmeXG8pMhdpE4LoZ6dO3cyevRoAObMmUOzZs1Y/2cIx66HY2mmZVrbsnLCIQMZ3cyeP3/+pds1Gg1WVlYUKlRILjAQIpPEx8fTpk0bbt++TbFixdi6dStRCQZ6rvqT6IRkyns60VuGF+R6UqeFUM/58+fp2rUriqIwYMAAhgwZwu9/h/Hpjr8AGNm4BEVcbFVOmbMY3cxWqFAh9d3E8wtP/vnuwtzcnE6dOrF06VKsrGS+NCEyiqIo9O/fn2PHjuHk5MSuXbswt7Gny/I/uPEohgKOVizqVgmdXPCV60mdFkIdDx8+JCAggJiYGBo2bMj8+fM5ffsJA78/TbJBIaB8AfrV9lY7Zo5j9JjZ7du3U7x4cZYtW8a5c+c4e/Ysy5Ytw8fHh3Xr1rFy5Up+++03xo8fnxl5hci1Tp06xbp169DpdGzevJkiRYvR75tT/HUvijy2FnzXr6qs8CUAqdNCqGXu3LmEhIRQokQJNm/ezPXwON5dc5L4JAN1S7gyq0N5OeGQCYw+Mzt16lS+/vprmjZtmrqtXLlyeHh48Omnn/Lnn39ia2vLhx9+yMyZMzM0rBC5WeXKlfnxxx8JDQ2lUaNGzNz7NyduRmBnaca371ahqKud2hGFiZA6LYQ6pkyZgk6no1evXtjaOzJszRGi4pPxL+zMku5+WJgZfQ5RpIPRzWxwcDCFCxd+YXvhwoUJDg4GUj7ier4WuBAi47Ro0QKAS6FRLDl4HYAZ7ctRpqCjmrGEiZE6LYQ6dDodU6ZMAWDer1f5+2E0eW0tWNbTH2sLncrpci6j3yL4+voyffp0EhMTU7clJSUxffp0fH19Abh37x5ubm4Zl1KIXOr+/fs0adKEGzdupG7TGxTGbj1PskGhaWk3mpeVJRBFWlKnhcg6W7ZsoV+/fmn+3q6FRbPgt5SVGCcElCKPrYVa8XIFo8/MLly4kFatWuHh4UG5cinrvZ8/fx69Xs+uXbsAuHHjBu+//36GhxUiN4mNjaV169acOnWKPn36cPDgQQBWH72ZunrM5NZlVE4pTJHUaSGyxqlTp+jZsydxcXFUqFCBDz74AINBYczWYBL1Bhr45qNV+QJqx8zxNMrzS12N8OzZM77//nuuXLmCoij4+vrStWtX7O1NfyWLqKgoHB0diYyMxMHBQe04QryUwWCgc+fObN68GRcXF06cOIG3tzc3Hj3jnXlHiEvSM61tWbpUKaR21Gwrp9cCqdNCZK67d+9SpUoVQkNDadGiBTt37kSn07HyyE0+33URWwsd+0bWlQtz35AxdeCNFk2ws7Nj4MCBbxROCPHfJk2axObNmzE3N2fbtm14e3sTl6jn/bVniEvSU907L51k9RjxGlKnhcg8MTExtGrVitDQUMqUKcP69evR6XQEhTxh+s+XABjb3Fca2SzyRs0swMWLFwkJCUkzRgSgVatWbx1KiNxs3bp1fP755wAsW7aM2rVrAzBx519cfhCNi50FX3eugFamdxH/Qeq0EBnPYDDQs2dPgoKCcHV15ccff8TBwYGnsYl8sC6IJL1C8zLudK/24kWYInMY3czeuHGDNm3aEBwcjEajeWFCbr1en7EJhchFTpw4wbvvvgvA6NGj6d27NwCbT91h06m7aDUwr3NF8jnIRPfi1aROC5F5JkyYwLZt27CwsGD79u0UKVIEg0Hhw03nuPc0jiJ5bfiyfTlZrjYLGT2bwbBhw/Dy8uLhw4fY2Nhw4cIFDh06hL+/P7///nsmRBQi9yhcuDDly5endevWTJs2DYBrYc/49IeUZRBHNCpBjWIuakYU2YDUaSEyT506dXB0dGTFihXUrFkTgFVHb/Lr5TAszLQs7FYJBytzlVPmLkafmT1+/Di//fYbrq6uaLVatFottWrVYtq0aQwdOpSgoKDMyClEruDu7s7vv/+OwWBAq9WSpDcwYuNZ4pMM1CrmwuD6xdSOKLIBqdNCZJ4mTZpw7do1XFxSTixcfhDFV3v+BmBCy1KULiDzfmc1o8/M6vV67OxSVhpycXHh/v37QMoZpb///jtj0wmRC+j1en777bfU+9bW1tja2gIw/9erBN+LxNHanJkdyss4WZEuUqeFyFi3bt3i6tWrqfefN7IJyXqGbzhLot5AQ998dKsqM8yowehmtkyZMpw/fx6AqlWr8tVXX3H06FEmT56Mt7d3hgcUIqcbO3YsDRs2ZPLkyWm2nwl5woIDKZNuT21TBndHGScr0kfqtBAZJyoqioCAAKpWrcqxY8fS7Ju9/wqXH6Ss8jW9nYyTVYvRwwzGjx9PTEwMkLIGccuWLalduzZ58+Zl48aNGR5QiJxs5cqVzJw5E4DixYunbo+OT2LExrMYFGhTsSAty8mk2yL9pE4LkTH0ej1du3blr7/+wt3dnUKF/v/M6/Hr4Sw7lLI647S2ZXG1t1QrZq5ndDPbtGnT1H97e3tz8eJFIiIicHZ2lnckQhjh4MGDDBo0CEi5OrZLly4AKIrC2K3B3A6PpYCjFZNalVYzpsiGpE4LkTFGjx7NTz/9hJWVFTt37sTDwwOAsOh4hqwPQlGgk78nTUq7q5w0dzNqmEFycjJmZmb89ddfabbnyZNHCqQQRrh+/Tpt27YlKSmJjh07MnHixNR93x6/zU/BoZhpNSzoVglHa7kqVqSf1GkhMsaKFSuYPXs2AN988w2VK1cGQG9QGLo+iMfPEvBxs5cTDibAqGbWzMyMwoULyxyFQryFp0+f0rJlSyIiIvD392f16tVotSl/iufuPGXKTxcBGNeiJJUKOasZVWRDUqeFeHsHDhxI/eRs0qRJdOzYMXXfnP1X+ONGBLYWOhZ2q4S1hU6tmOJ/jL4AbPz48YwbN46IiIjMyCNEjvfjjz9y+fJlChYsyM6dO7GxsQHg2LXH9Pv2FEl6hWal3Xm3ZhF1g4psS+q0EG9nzpw5JCcn07lzZyZMmACknJGd9+vV1Atzv2hblmL57NSMKf7H6DGz8+bN49q1axQoUIDChQunTiH03JkzZzIsnBA5UY8ePdBoNJQuXZr8+fOTmGxg9v4rLD10HUUBHzd7vuogV8WKNyd1Woi3s2nTJmbMmMGoUaPQaDTcexrH8A1BnLz1BIB+tbxoXaGgyinFc0Y3s4GBgZkQQ4icT1GU1Aa1e/fuAIRFxfPed6c5e+cpAF2qFOLTliWxsTD6T1OIVFKnhTDeP2u0lZUVn376KQCHrjzig3VniIpPxs7SjM8DS9OmooeaUcW/aJTni3bnElFRUTg6OhIZGYmDg4PacUQusX//fiZPnszmzZtxd0+56vXC/Uj6fXOK0Mh4HK3N+bJdOZqVkStis4rUAtMlPxuhhsGDB2Nvb88XX3yReh3DN8duMXnXRfQGhQqeTszrXJFCeW1UTpo7GFMHjB4zCykXsKxYsSLNmKwzZ85w7969N3k6IXK0y5cv06FDB44cOcKMGTMA+Dk4lA5LjhMaGU9RV1t+GFxTGlmRoaROC5F+CxYsYNGiRXz11VecOHGChGQ943cEM3HnBfQGhXaVPNg4oJo0sibK6M8yz58/T6NGjXB0dOTWrVv079+fPHnysH37dm7fvs23336bGTmFyJbCw8Np2bIlkZGR1KhRg08nfc7H24NZdyIEgNrFXVjQVabfEhlL6rQQ6bd3716GDRsGwPTp03EtWpY2C49xMTQKjQbGNPNlQB1vuY7BhBl9ZnbkyJH07t2bq1evYmX1/8trNm/enEOHDmVoOCGys8TERNq3b8/169cpUqQIs1aso/3yk6mN7IA63qzuXVkaWZHhpE4LkT4XL16kY8eOGAwG+vTpQ5H6nWk5/zAXQ6NwtjFnRU9/BtYtKo2siTO6mT158iQDBgx4YXvBggV58OCB0QEWLVqEl5cXVlZW+Pn5cfjw4dcen5CQwCeffELhwoWxtLSkaNGirFq1yujXFSIzKYrC4MGD+f3337G3t2frjp2M/vEG18Kekc/eku/7VmVci5KY6d5opI8QryV1Woj/9vjxYwICAoiKiqJ27dr0GzeN0VuDiU8yUKuYC3uG16FhSTe1Y4p0MHqYgZWVFVFRUS9s//vvv3F1dTXquTZu3Mjw4cNZtGgRNWvWZOnSpTRv3pyLFy+mWf/4nzp27MjDhw9ZuXIlxYoVIywsjOTkZGO/DCEy1cKFC1mxYgVarZYNGzZw+JEVIRGxuDlYsntobfLayRreIvNInRbi9QwGA+3atePGjRt4eXmxactW+qy/DEDbigWZ2aE8Wq2cjc02FCP1799fCQwMVBITExU7Ozvlxo0byu3bt5WKFSsqw4YNM+q5qlSpogwcODDNNl9fX2Xs2LEvPf7nn39WHB0dlfDwcGNjp4qMjFQAJTIy8o2fQ4j/cu3aNcXX11eZM2eOEhIeo5T4ZLdSeMwu5Yez99SOJv4nJ9cCqdNC/Lf169crbm5uyoULF5Q1R28qhcfsUspN2quEP0tQO5pQjKsDRn/GOXPmTB49ekS+fPmIi4ujbt26FCtWDHt7e6ZOnZru50lMTOT06dM0adIkzfYmTZpw7Nixlz5m586d+Pv789VXX1GwYEFKlCjBqFGjiIuLe+XrJCQkEBUVleYmRGYrWrQop0+fZtiwYXy+6yIJyQaqe+cloFx+taOJXEDqtBD/rXPnzty4cYN8hYoya9/fAIxq6kMeWwuVkwljGT3MwMHBgSNHjvDbb79x5swZDAYDlSpVolGjRkY9z+PHj9Hr9bi5pR2P4ubm9soxXTdu3ODIkSNYWVmxfft2Hj9+zPvvv09ERMQrx2NNmzaNzz77zKhsQryJsLAwzp07R+PGjQGwsbFh/8WH7Lv4EDOths9al5aLCESWkDotxMsdOHAAHx8fChQoAICVlTXjN58jKj6ZMgUd6Frl5UNnhGkzupm9desWRYoUoUGDBjRo0OCtA/z7P3flHytw/JvBYECj0bB27VocHR0BmD17Nu3bt2fhwoVYW1u/8Jhx48YxcuTI1PtRUVF4enq+dW4h/ik+Pp7AwEBOnDjBqlWr6N6jJ4t/v8acX64C0LtGEUq42aucUuQWUqeFeNG5c+cICAjA0dGRw4cPY+9akFGbz3HwyiMAPmtVBp2Mk82WjB5m4O3tTa1atVi6dGnqRNxvwsXFBZ1O98K7+7CwsBfOAjyXP39+ChYsmFogAUqWLImiKNy9e/elj7G0tMTBwSHNTYiMpCgK/fr14/jx4zg4OFCifGW6LP+DmfuuoDcotCpfgA+b+KgdU+QiUqeFSOvBgwcEBAQQExNDyZIluRVvTbO5hzh45RGWZlq+bFcWv8LOascUb8joZvbUqVNUr16dKVOmUKBAAVq3bs3mzZtJSEgw6nksLCzw8/Nj//79abbv37+fGjVqvPQxNWvW5P79+zx79ix125UrV9BqtXh4yDrJQh3Tpk1j7dq16HQ6Vq3bzCe/hPHnzQhsLXTM7lierztXwNpCp3ZMkYtInRbi/z3/5OzOnTuUKFGCodOX8d73QYTHJOLrbs+uIbXoVFmGF2Rrb3qVmcFgUH777TelX79+irOzs+Lg4KD06dPHqOfYsGGDYm5urqxcuVK5ePGiMnz4cMXW1la5deuWoiiKMnbsWKVHjx6px0dHRyseHh5K+/btlQsXLigHDx5UihcvrvTr1y/drylXyYqMtGXLFgVQAGXOwiVKi68PKYXH7FJqTPtVufnomdrxxGvkhlogdVrkdgaDQenSpYsCKM7Ozsq2Q+cUn/Eps8sMXntaiUtMVjuieAVj6sAbN7P/dPr0aaVChQqKVqs1+rELFy5UChcurFhYWCiVKlVSDh48mLqvV69eSt26ddMcf+nSJaVRo0aKtbW14uHhoYwcOVKJjY1N9+tJkRQZ5dSpU4q1tbUCKIOHDlc6LDmmFB6zS6k0eZ9yPSxa7XjiP+S2WiB1WuRGkydPVgDFzMxMWbPjF6XMxD1K4TG7lN6rTiiJyXq144nXMKYOaBRFUd7kjO6dO3dYv34969atIzg4mOrVq9OtWzcGDRqUAeeLM09UVBSOjo5ERkbKuCzxVsaPH8/UqVNp2qwZPj2n8sO5UOwtzVj/XjXKFHT87ycQqsoNtUDqtMjN4uPjqVKlCsHBwcxZvIINEUV4EBVP5SLOfPtuVRn+ZeKMqQNGz2awbNky1q5dy9GjR/Hx8aFbt27s2LGDIkWKvGleIbKlzz//nKJFi2LwqsHne66h02pY2tNPGlmhOqnTQqSshHfkyBE2bd7CQUrzIOox3q62rOhVWRrZHMboZvbzzz+nc+fOfP3111SoUCETIglhugwGAwaDATMzMzQaDZWbtKXt4pTJ40c18aFGUReVEwohdVrkbklJSZibmwMpcy7HFqnN4f1XsDLXsribH47W5ionFBnN6GY2JCREJn4XudbEiRM5efIkGzduBAsbBq87Q2KygYa++RhQx1vteEIAUqdF7hUTE0O9evXo1KkTH374IceuhzPnlysATAksi4+7zPedExndzGo0Gp4+fcrKlSu5dOkSGo2GkiVL0rdv3zTzCgqR06xbt44pU6YAsGP3PraF5+d2eCwFnayZ1bE8WplsW5gIqdMiNzIYDPTs2ZNTp05x+/ZtyjZsw8jtV1EU6OjvQXs/mRoup3qjeWaLFi3KnDlziIiI4PHjx8yZM4eiRYty5syZzMgohOqOHz/Ou+++C8DQ0Z+wIcydc3cjcbIxZ1lPP5xsZC1vYTqkTovcaPz48Wzbtg0LCws+W7qJYduuEJOop7p3Xia3LqN2PJGJjJ7NoHbt2hQrVozly5djZpZyYjc5OZl+/fpx48YNDh06lClBM4pcJSuMdfv2bapUqUJYWBjN23UlsXp/rj2KwcXOkrX9qsrHVtlUTq4FUqdFbvPdd9/Rs2dPAMYuWM/W+44k6g3U93FlcXc/rMzlgq/sxpg6YHQza21tTVBQEL6+vmm2X7x4EX9/f2JjY41PnIWkSApjREdHU7NmTYKDgyld5x0sGw0lPCaJ/I5WrO1XFW9XO7UjijeUk2uB1GmRmxw9epQGDRqQmJhI67ELOEcRFAWal3Hn684VsTAz+kNoYQKMqQNG/4QdHBwICQl5YfudO3ewt5czVCJnee+99wgODiZ/tQASa71PeEwSvu72bBlUQxpZYbKkTovcIjw8nDZt2pCoN1Bx4GzOKimNbPdqhZjfRRrZ3MLoC8A6depE3759mTlzJjVq1ECj0XDkyBE++ugjunTpkhkZhVDNmDFjOPlIQ7J/NxL1Co1K5mNu54rYWRr9pyNElpE6LXKLvHnzMmHCROafiSPCsQRaDUwMKE2vGkXUjiaykNH/I8+cORONRkPPnj1JTk4GwNzcnEGDBjF9+vQMDyiEmvTOhdFU7Q56hZ7VCzMxoDQ6mbVAmDip0yI3iS3WkIS717HQaVnW0496PvnUjiSy2BsvZxsbG8v169dRFIVixYphY2OT0dkyhYzFEv/l8OHD6HQ6CpWsQODCozx+lkiTUm4s6e4n02/lILmhFkidFjnV8uXLad++PQduxjBy0zkA5nQqT5uKMv1WTpEpY2b1ej3nz58nLi4OABsbG8qWLUu5cuXQaDScP38eg8HwdsmFUNn169cJDAykfqMmdFl0kMfPEimZ34E5nSpIIytMntRpkRssX76c9957j6rNOzJm63kABtcvKo1sLpbuZva7777j3XffxcLixfk0LSwsePfdd1m3bl2GhhMiKz19+pSWLVsSERFBicZduRNtwMXOghW9/LGVMbIiG5A6LXK6AwcO8P777wPg2qAPSXqFJqXc+LCxj8rJhJrS3cyuXLmSUaNGodO9OFebTqdj9OjRLFu2LEPDCZFVkpOT6dSpE5cvX6ZgwYIUqp9ykUzfWt4UdLJWOZ0Q6SN1WuRkV69epV27diQnJxPYrS/3lZTV7D55p6R8cpbLpbuZ/fvvv6lWrdor91euXJlLly5lSCghstqIESPYt28fNjY2zP9uGxcexGKm1cjyhyJbkTotcqonT57QsmVLnjx5QtWqVanaZQSKArWKuVA4r63a8YTK0t3MxsTEEBUV9cr90dHRJj8RtxAvs2jRIhYsWADA999/z9molItkmpR2w9XeUs1oQhhF6rTIiZKSkujYsSNXrlzB09OTLdu2s+1sKABdqhRSOZ0wBeluZosXL86xY8deuf/IkSMUL148Q0IJkVUURUld2nPatGk0e6cV24LuAVIkRfYjdVrkRKGhoVy9ehVbW1t+/PFHLj7REBadQF5bCxqXclM7njAB6W5mu3btyvjx4zl//vwL+86dO8eECRPo2rVrhoYTIrNpNBrWrVvH5s2bGTNmDD8FhxIdn4xnHmtqFnVRO54QRpE6LXKiQoUK8eeff/Ljjz9Svnx51v+Zsrpde38PWeFLAEYsmjBixAh+/vln/Pz8aNSoEb6+vmg0Gi5dusQvv/xCzZo1GTFiRGZmFSLDxMTEYGNjg0ajQavV0r59e4DUItm5ciG5oEBkO1KnRU7y7Nkz7OxSlg3Ply8f+fLl497TOA5eeQSk1GkhwIgzs+bm5uzbt4+pU6cSGhrKsmXLWLJkCaGhoUydOpV9+/Zhbm6emVmFyBCJiYm88847dOvWLXU+ToBd5+9z+vYTzLQaOvjLhV8i+5E6LXKKixcv4u3tzZo1a1K3KYrCF7svYVCgRtG8eLnIhV8ixRuvAJZdycoyuZuiKLz33nusWLECe3t7/vzzT3x9fTl35ykdlx4nIdnAgDrejGtRUu2oIpNJLTBd8rPJ3R4/fkzVqlW5ceMGtWvX5sCBA+h0Oub+coW5v1zFXKdhff9q+BfJo3ZUkYkyZQUwIXKCOXPmsGLFCrRaLRs2bMDX15cHkfH0//YUCckG6vu4MrqZr9oxhRAiV0pMTKRt27bcuHEDLy8vtm3bhk6n48dz95n7y1UApgaWlUZWpCHNrMg1du3axahRowCYNWsWLVq0ICFZT/9vTxEWnUAJNzvmdamITsbKCiFEllMUhYEDB3L48GEcHBzYtWsXLi4u/HUvklGbzwHQv7YXHSt7qpxUmBppZkWuEBwcTJcuXVAUhf79+zNs2DAAZu+7QvC9SJxtzFnRszL2VjKeUAgh1DBz5kxWr16NVqtl48aNlCpVivgkPcM2BKV+cja2uQwBEy+SZlbkeMnJybRr145nz55Rv359Fi5ciEaj4Y8b4Sw7fAOAL9uVo1BeG5WTCiFE7vTnn38yZswYAObOnUuzZs0AmP7zZa4/iiGfvSWzO1aQT87ES6V7ai4hsiszMzOWL1/OuHHj2LJlC+bm5kTFJ/HhpnMoCnTy96RJaXe1YwohRK5VuXJlJkyYQFhYGB988AEAh68+Ys2xWwB81b4czrYWKiYUpixdzezIkSPT/YSzZ89+4zBCZJa6dety9OhRNJqUd/Wf7bzIvadxeOax5tOAUiqnE+LtSZ0W2ZlGo2HSpEkoioJGoyEyNil1nGyPaoWp55NP5YTClKWrmQ0KCkpz//Tp0+j1enx8fAC4cuUKOp0OPz+/jE8oxBtavHgxderUoXTp0gCpjeyx64/ZeuYuWg3M6VgBO0v5gEJkf1KnRXYTFxfHtGnTGDt2LDY2KcO8ntfpGfsu8zAqAW8XWz6WqRLFf0jX/+IHDhxI/ffs2bOxt7fnm2++wdnZGYAnT57Qp08fateunTkphTDS1q1bef/997G3t+fChQt4eqZc/ZqkNzDxhwsAdK9WWKZ3ETmG1GmRnSiKQt++fVm/fj0nTpxg7969qfv+uhfJ2hMpqzF+0bYs1hY6tWKKbMLoC8BmzZrFtGnTUgskgLOzM1OmTGHWrFkZGk6IN3H69Gl69OgBQJ8+fVIbWYBvjt3iatgz8tha8GFjH7UiCpGppE4LUzdlyhTWr1+PmZkZY8eOTd1uMCh8+sNfKAq0Kl+Aat55VUwpsgujm9moqCgePnz4wvawsDCio6MzJJQQb+revXu0atWKuLg4mjVrluY/7rCo+NRJt8c088HRRqbhEjmT1GlhyjZt2sSECRMAWLRoEfXr10/dt/XMXYJCnmJroeOTd2R4gUgfo5vZNm3a0KdPH7Zs2cLdu3e5e/cuW7ZsoW/fvrRt2zYzMgqRLrGxsbRu3Zr79+9TqlQpNmzYgJlZykgag0Fh0o8XeJaQTAVPJzr4yaTbIueSOi1M1cmTJ+nVqxcAI0aMoH///qn7wqLj+XLPZQCGNSqOm4OVKhlF9mP0lS9Llixh1KhRdO/enaSkpJQnMTOjb9++zJgxI8MDCpEeBoOBXr16cfr0aVxcXPjxxx9xdHQEIFlvYMzWYHYHP0CrgcmtS6OVuQpFDiZ1Wpiiu3fv0rp1a+Lj42nRokWa38V7T+PotvwPHj9LpFg+O3rX8FIxqchujG5mbWxsWLRoETNmzOD69esoikKxYsWwtbXNjHxCpEtsbCwRERGYm5uzbds2vL29gZQLvoZvPMtP50PRaTXM6lCech5O6oYVIpNJnRam6MGDByQnJ1OmTBnWr1+PTpdyYdft8Bi6Lj/BvadxFHSyZmUvfyzMZE0nkX5vPCdRaGgooaGh1KlTB2tr69S54YRQg52dHXv27OHPP/+kZs2aKIrCoauPmfvLFYJCnmKu0zC/SyWalZHFEUTuIXVamBJ/f3/+/PNPABwcHIhNTGb9n3dYdOAa4TGJeLvY8n2/qhRwslY5qchujG5mw8PD6dixIwcOHECj0XD16lW8vb3p168fTk5OcqWsyFIPHz7Ezc0NAHNzc2rWrMnJWxFM/vEiwfciAbAy17K4ux/1ZdJtkUtInRam5J91ukiRIiiKworDN1j0+3UiYhIB8HW359u+VchnL+NkhfGMPo8/YsQIzM3NCQkJSZ3kGKBTp07s2bMnQ8MJ8Tq3b9+mXLlyDBkyhOTkZABO346gx8oTBN+LxNpcR99aXvw+qr40siJXkTotTMV3331HsWLF+PHHH1O3zdp3hSk/XSIiJpFCeWyY1rYsP3xQUxpZ8caMPjO7b98+9u7di4eHR5rtxYsX5/bt2xkWTIjXiY6OJiAggLCwMA4fPkx8fDx3og30WX2S+CQDdUu4MqdTBfLIWt4iF5I6LUzB0aNH6devH4mJifzxxx8EBASw4vANFhy4BsAnLUrSp2YRzHQyPla8HaOb2ZiYmDTv9J97/PgxlpaWGRJKiNfR6/V07dqV4OBg3Nzc2LlzJ+EJGnqu+pOo+GT8CzuzpLufrBojci2p00Jtt27dok2bNiQmJtK2bVs+//xzNp26w5SfLgHwUVMf+tfxVjmlyCmMfjtUp04dvv3229T7Go0Gg8HAjBkz0kx8LERmGTt2LLt27cLS0pIffvgBT09PPtp8nkfRCZTM78DK3pWlkRW5mtRpoaaoqChatmzJo0ePqFixIt9++y2hUQl8sj0YgPfqePN+vaIqpxQ5idFnZmfMmEG9evU4deoUiYmJjB49mgsXLhAREcHRo0czI6MQqVauXMnMmTMBWLNmDVWrVuX49XD+vBWBhU7Lyl7+OFrLyl4id5M6LdSi1+vp0qULFy5cIH/+/OzcuRNbW1um7QgmSa9QzTsP45r7yqwaIkMZfWa2VKlSnD9/nipVqtC4cWNiYmJo27YtQUFBFC0q77RE5rl//z6DBw8GYOLEiXTu3BmAeb+mLFHbuYqnTOkiBFKnhXqWL1/O7t27sbKy4ocffsDDw4MHkfFsOnkXgOGNSkgjKzKc0WdmQ0JC8PT05LPPPnvpvkKFCmVIMCH+rUCBAmzcuJEff/yRiRMnAnDyVgTHb4RjrtMwsK78Jy0ESJ0W6unXrx8XLlygdu3aVK5cGYAlB6+TqDdQpUgeqnnnVTmhyImMbma9vLwIDQ0lX760Ux2Fh4fj5eWFXq/PsHBC/Fvr1q1p3bp16v3nZ2Xb+8lZWSGekzot1GJmZsb8+fNT74dFx7P+zxAAhjYsrlYskcMZPczgVSvIPHv2DCsrmSNOZKzk5GSGDh360umEzoQ84fDVx5hpNXIxgRD/IHVaZKWrV68yatQokpKSXti3/NANEpINVCrkRM1iclZWZI50n5kdOXIkkHJV7Keffppm2he9Xs+JEyeoUKGC0QGerx8eGhpK6dKlmTt3LrVr1/7Pxx09epS6detSpkwZzp49a/Triuxh+PDhLFy4kF27dnH58mUsLFLmjVUUhS9/vgxA20oF8czz4jREQuQ2UqdFVnvy5AktW7bkypUrGAwGZs+enbrvTkQs3x5PORExtGFxGSsrMk26m9mgoCAgpYkIDg5ObSoALCwsKF++PKNGjTLqxTdu3Mjw4cNZtGgRNWvWZOnSpTRv3pyLFy++dkxXZGQkPXv2pGHDhjx8+NCo1xTZx8KFC1m4cCEAs2bNSvM7t/PcfU7cjMDKXCsfXQnxP1KnRVZKSkqiQ4cOXLlyBU9PT0aPHp1m/+e7LpKQbKCadx7qlnBVKaXIDTSKoijGPKBPnz58/fXXODg4vPWLV61alUqVKrF48eLUbSVLliQwMJBp06a98nGdO3emePHi6HQ6duzYYdQ7/qioKBwdHYmMjMyQr0Fkjn379tGiRQv0ej3Tpk1j7Nixqfui45NoOOsgYdEJfNi4BEOkmRVvICfXAqnTIrMpisL777/PkiVLsLW15ejRo5QvXz51/+9/h9F79Ul0Wg27h9bGx91exbQiOzKmDhg9Znbu3LkkJye/sD0iIoKoqKh0P09iYiKnT5+mSZMmabY3adKEY8eOvfJxq1ev5vr166lXs/+XhIQEoqKi0tyEabt8+TIdO3ZEr9fTs2dPxowZk2b/vF+vEhadQOG8NrKCjBAvIXVaZLYFCxawZMkSNBoN69atS9PIJiTrmbTzAgC9axSRRlZkOqOb2c6dO7Nhw4YXtm/atCl13s/0ePz4MXq9Hjc3tzTb3dzcePDgwUsfc/XqVcaOHcvatWsxM0vfCIlp06bh6OiYevP09Ex3RpH1wsPDadmyJZGRkdSsWZNly5alGWd15WE0q4/eAmBSQGmszGWlLyH+Teq0yEx79uxh+PDhAEyfPp1WrVql2b/i8E1uhcfiam/J8EbyyZnIfEY3sydOnHjpcoj16tXjxIkTRgf494DwV12Fq9fr6dq1K5999hklSpRI9/OPGzeOyMjI1NudO3eMziiyTkJCAo6OjhQpUoTt27enWUc+LCqeft+cItmg0KikG/V9873mmYTIvaROi8ykKAq2trb06dOHjz76KM2+/RcfMnv/FQA+buGLvZWsyCgyn9HzzCYkJLz046ukpCTi4uLS/TwuLi7odLoX3t2HhYW9cBYAIDo6mlOnThEUFMQHH3wAgMFgQFEUzMzM2LdvHw0aNHjhcZaWlmkaImHaChQowKFDh3jw4AGurv9/wUBkbBI9V/1JSEQshfLY8EWbMiqmFMK0SZ0Wmal58+acPHmSIkWKpHlTc/x6OIPXnUFvUGhXyYPACgVVTClyE6PPzFauXJlly5a9sH3JkiX4+fml+3ksLCzw8/Nj//79abbv37+fGjVqvHC8g4MDwcHBnD17NvU2cOBAfHx8OHv2LFWrVjX2SxEm5Nq1a6n/trW1TbPkZmxiMn3W/MnlB9Hks7fk+75Vyecgc2UK8SpSp0VGS0xMTDPft4+PT5o3IOfvPqX/t6dITDbQuJQbX7YrK1NxiSxj9JnZqVOn0qhRI86dO0fDhg0B+PXXXzl58iT79u0z6rlGjhxJjx498Pf3p3r16ixbtoyQkBAGDhwIpHz0dO/ePb799lu0Wi1lyqQ9G5cvXz6srKxe2C6yl127dhEYGMiECRP49NNP0xTA6Pgk3l1zkjMhT3GwMuPbvlUolFfmlBXidaROi4ykKAoDBw5k586dbN++/YU5hoNCntBr1Z88S0imunde5nepiJnO6HNlQrwxo5vZmjVrcvz4cb766is2bdqEtbU15cqVY+XKlRQvbtxA706dOhEeHs7kyZMJDQ2lTJky7N69m8KFCwMQGhpKSEiIsRFFNnL+/Hm6dOmCXq/n7t27afY9jU2k16o/OXc3EnsrM9a8WwVfd5mmR4j/InVaZKSZM2eyevVqtFotMTExafaduBHOu2tOEpOox6+wM8t7+cuFuSLLGT3PbHYn8xeajocPH1KlShVCQkKoX78+e/fuxdw85WKBiJhEui7/g8sPonG2Mee7vlUpU9BR5cQiJ5FaYLrkZ2M6du7cSWBgIIqiMG/ePIYMGZK678jVx/T79iTxSQZqFM3Lil7+2FgYfY5MiJfK1HlmAa5fv8748ePp2rUrYWFhQMpUHRcuXHiTpxO5UHx8PG3atCEkJITixYuzZcuW1EY2Oj6JXqtSxsi62luy4b3q0sgKYSSp0+JtnT17lq5du6IoCoMGDUq9qA/g9O0I+n97ivgkA/V9XFnVu7I0skI1RjezBw8epGzZspw4cYKtW7fy7NkzIOXj4vROkC1yN0VR6NevH8ePH8fJyYldu3aRJ08eAOKT9PT/9hTB9yLJY2vB+v7VZMJtIYwkdVq8rQcPHtCqVStiYmJo1KgRX3/9der1DJdCo+iz+iRxSXrqlnBlaQ8ZWiDUZXQzO3bsWKZMmcL+/fvTrPtdv359jh8/nqHhRM504MAB1q5di06nY8uWLanzUSYmGxiyPog/bkRgZ2nGt+9WoVg+O5XTCpH9SJ0Wb2vSpEncuXMHHx8fNm3alPrJ2a3HMfRc9SdR8cn4FXZmSXc/LMzkYi+hLqM/EwgODmbdunUvbHd1dSU8PDxDQomcrUGDBqxatYrk5OTUK62fxCQyaO1p/rgRgYWZlhW9/GVogRBvSOq0eFtz5sxBURRGjRqFs7MzACdvRTDgu9NExCTi627Pql6VsbaQM7JCfUY3s05OToSGhuLl5ZVme1BQEAULygTJIn369OmT+u9rYc/o981JboXHYmuhY2G3SlTzzqtiOiGyN6nT4m1ZW1uzdOnS1PtbTt9l3LbzJOkVyhR0YFXvyjjayOpewjQY/dlA165dGTNmDA8ePECj0WAwGDh69CijRo2iZ8+emZFR5AD37t2jY8eOPHr0KM3264+e0XbRUW6Fx+LhbM2292tSz0eWqRXibUidFm9i06ZNfPrppxgMhjTbvz1+i1Gbz5GkV2hexp3NA2qQz14WrhGm440WTejduzcFCxZEURRKlSqVuh73+PHjMyOjyOZiY2Np3bo1p0+fJjY2ll27dgGQpDcwcuNZouKTKe/pxMpe/rjYyZKWQrwtqdPCWCdPnqRXr17Ex8dTrFgxevXqBcDfD6KZ8tMlAAbWLcropj5otbKylzAtbzzP7PXr1wkKCsJgMFCxYkWjJ+JWi8xfmLUMBgMdO3Zk69atuLi4cOLECby9vQGYs/8KX/96FUdrc/YOr4O7o7zTF1knN9QCqdMiPe7evUuVKlUIDQ2lRYsW7Ny5E51OR0KynsCFx7gUGkUD33ys7OUvS9SKLGNMHXjjSeGKFi2a2pTIL7d4lQkTJrB161bMzc3Ztm1b6u9MUMgTFhy4BsCUwDLSyAqRCaROi/8SExNDq1atUld3W79+PTpdykVdc/Zf5VJoFHlsLZjerqz8DgmT9UbzaaxcuZIyZcpgZWWVuub2ihUrMjqbyObWrl3L1KlTAVi+fDm1a9dGURQOXnnE0A1B6A0KrSsUIKB8AZWTCpHzSJ0W/8VgMNCjRw+CgoJwdXXlxx9/xMHBgdjEZJYfusHSQ9cB+KJNWRkjK0ya0WdmP/30U+bMmcOQIUOoXr06AMePH2fEiBHcunWLKVOmZHhIkf0cP36cvn37AjBmzBg6de3OrvP3WXrwBsH3IgEo6GTN5FZl1IwpRI4kdVqkx/jx49m+fTsWFhZs374dmzzuLPjtKquO3iIiJhGAjv4eNCvjrnJSIV7P6DGzLi4uzJ8/ny5duqTZvn79eoYMGcLjx48zNGBGk7FYWePKlSu0bNkSr0p1qNB+CLv/ekB0fDIA1uY6ulYtxIC63vJuX6gmJ9cCqdMiPb7//nv69evPBzPW8NCmCEevPcbwv46gUB4bBtYtSkd/D8x0siiCyHqZOmZWr9fj7+//wnY/Pz+Sk5ONfTqRQ5UoUYIp3/zExz9e5e9TdwHI72hFBz8PetUoQl6ZtUCITCN1WqRH5y5dOZLszZbLT4CUNziVCjnRs3oRWpbLL02syDaMbma7d+/O4sWLmT17dprty5Yto1u3bhkWTGQ/er2eoKAg/P392XbmLmN3XsGgQH0fV96rU5SqXnlkShchsoDUafEqt2/fxtLSkryu+Ri+8Sx7Lj9Bp9Xwfr2itPfzoHBeW7UjCmG0N5rNYOXKlezbt49q1aoB8Mcff3Dnzh169uzJyJEjU4/7dyEVOdvYsWOZM2cO/aevZk94HhQFOvh5ML1dOXTSxAqRpaROi3+LiorinXfeIepZDNVGreTPu3GY6zTM71KRZmXyqx1PiDdmdDP7119/UalSJSBlDkNIWe/b1dWVv/76K/U4mcIjd1m5ciUzZ87E0qM0ex87owA9qxdmUkBpORsrRBaTOi3+Ta/X06VLFy5cuIBHwDD+vBuHpZmWJT38qC+rLopszuhm9sCBA5mRQ2RjBw8eZODAgWhtnPDq9jlxaGhTsSCftSot/1kKoQKp0+LfPvroI3bv3o1jmbroSjUGYH6XitLIihzB6NHdDx8+fOW+8+fPv1UYkf1cu3aNtm3bkqw3ULrfDOKwoHg+O6a2KSONrBAqkTot/mn58uXMmTMHM6f8uAeMAmBAHW+alJYpt0TOYHQzW7ZsWXbu3PnC9pkzZ1K1atUMCSWyh6dPnxIQEEBERAQlO4wkyjo/NhY6FnevhI3FGy8uJ4R4S1KnxXMHDhzg/fffB50ZZQd9TbxBQ+Uizoxq6qN2NCEyjNHN7JgxY+jUqRMDBw4kLi6Oe/fu0aBBA2bMmMHGjRszI6MwUd988w2XL1/Go2x1ErzrASkrxRTLZ69uMCFyOanTAkBRFMaPH09ycjI1+39OhMGGPLYWzO9SCXOZdkvkIEYvmgBw7tw5unfvTnx8PBEREVSrVo1Vq1bh5uaWGRkzlEzGnXEURWHmrNn8rC/LjSdJvFMuPwu7VlI7lhDpktNrgdRpASmfoH00ZQ4HLKqSbFBY1K0SLcrKzAXC9BlTB97orZm3tzelS5fm1q1bREVF0bFjx2xRIEXG0mg02FcO5MaTJBytzZkUUFrtSEKI/5E6LQDsHRx55NWEZINCk1JuNJelaUUOZHQze/ToUcqVK8e1a9c4f/48ixcvZsiQIXTs2JEnT55kRkZhQvbt20ebNm2Ijo7m1uMY5uy/AsD4d0riai+reglhCqRO516KojBkyBDmzZuHoiisOXaLs3eeYm9pxuTWcmGuyJmMvkqnQYMGjBgxgs8//xxzc3NKlixJ/fr16dGjB2XLluXu3buZkVOYgEuXLtGxY0ciIyP58qsZ/F2wGQnJBmoWy0t7Pw+14wkh/kfqdO61YMECFixYgEajoWjFGszcl7JM7bgWJXF3tFI5nRCZw+hmdt++fdStWzfNtqJFi3LkyBGmTp2aYcGEaQkPDycgIIDIyEhq1qyJddWOnDxyGztLM6a1KSfv9oUwIVKnc6c9e/YwfPhwAKZM/4qvT8UQl6SnRtG8dK7sqW44ITLRG10Alp3JhQXGS0xMpHHjxhw6dIgiRYowa/1eRu64CsDCrpV4p5xcTCCyH6kFpkt+Nsa7ePEi1atXJyoqij59+pCn+VC2nL6Hq70lPw2tRT57OSsrspdMuQCsRYsWREZGpt6fOnUqT58+Tb0fHh5OqVKljE8rTJqiKAwaNIhDhw5hb2/Pqg07+GzvLQB61ygijawQJkTqdO70+PFjWrZsSVRUFLVr16bxe5+y5fQ9tBqY17miNLIix0t3M7t3714SEhJS73/55ZdERESk3k9OTubvv//O2HRCdXPnzmXVqlVotVo2bNjAkrOxRMYlUd7TiY9blFQ7nhDiH6RO5z5JSUm0bduWmzdv4uXlxYLV6/nsp8sAfNjEh+pF86qcUIjMl+5m9t+jEXLZ6IRcq3bt2hQoUIBZs2aRVKA8J25GYGWuZUGXiliYyaTbQpgSqdO5j5mZGW3btsXZ2Zldu3ax8Fgo8UkGqnrlYVDdomrHEyJLyJqj4rX8/f3566+/MLO2o9HsgwAMrlcMzzw2KicTQgih0WgYPnw4vXr14vyjZPZe+BOdVsPngWXQauXCXJE7pPvUmkajeeGKdbmCPWd6+PAhp06dSr3v7OzM/N+u8TAqgcJ5behfx1vFdEKIV5E6nXscOXIkzfhoG3sHJu28AKRcz1DCTZYVF7lHus/MKopC7969sbRMmRg/Pj6egQMHYmtrC5BmnJbIvuLj42nTpg1BQUFs2rSJgIAAroVFs+rITQAmBpTCylynckohxMtInc4dzp07R7NmzShUqBC//fYb7u7urDxyk5uPY3Cxs2R4o+JqRxQiS6W7me3Vq1ea+927d3/hmJ49e759IqEaRVHo168fx48fx8nJCR8fH+49jaPfN6dINig0KpmPBr6yHKYQpkrqdM734MEDAgICiImJoWDBguTNm5fdwaHM3peyGuPHLXyxtzJXOaUQWSvdzezq1aszM4cwAV988QVr165Fp9OxZcsWLPIUpOOS49x7GodnHmsmty6jdkQhxGtInc7Z4uLiCAwM5M6dO/j4+LB582Z2nn/IR1vOYVCgdYUCtKlYUO2YQmQ5uQBMALB161bGjx8PwMKFCylctgodlx4nLDoBbxdb1vavSn5Ha5VTCiFE7qQoCn379uXEiROpMxfs/juKj7cHA9DJ35Mv2paVMdIiV5JmVnD69Gl69OgBwLBhw2jRoQcdlqQ0sr7u9nzXtyqu9pYqpxRCiNzr888/Z/369ZiZmbFt2zb+irbm4+1ngZQLvia0LCWzF4hcSyYKFSxbtoy4uDiaN2/OuM++oMfKEzyIiqd4PjvW968mjawQQqgoOjqaFStWALB48WIM7iX5cPM5AHpVL8zEAGlkRe4mZ2YFixYtonjx4nTp+S591pzmVngsHs7WfNe3Ks62FmrHE0KIXM3e3p4TJ06wdetWyjVsS4+VJ9AbFNpULMjEgNIytEDkenJmNpcyGAypqwPpdDpGjRrFvEN3uRgahYudJd/3rYq7o6znLYQQajEYDKn/zp8/P737D2TI+jMkJBtoVDIfX7UvJ2dkhUCa2VxrwoQJ9O7dO3XeyaCQJ2w8dQeAJd0rUcTFVs14QgiRq8XExFCrVi3Wrl2bum3+r1dTF69Z0LUS5jr5L1wIkGEGudL333/P1KlTAWjfvj0t3mnJhB9SVo5p7+eBf5E8asYTQohczWAw0L17d44fP87169dp2bIljxK0rPzf4jWTAkrL4jVC/IO8rctljh8/Tt++fQEYM2YMAQEBbDgZQvC9SOytzBjTzFflhEIIkbuNHz+eHTt2YGFhwbZt23BwcGDizgupi9fU982ndkQhTIo0s7nI7du3CQwMJDExkcDAQL744guexCQyY+/fAIxsXEJmLhBCCBV9++23TJs2DYCVK1dSs2ZNdgc/4Oi1cCzMtExoWVrlhEKYHmlmc4no6GhatmxJWFgYFSpU4LvvvkOj0fDRlnM8jU3C192eHtUKqx1TCCFyrSNHjtC/f38APv74Y7p37879p3GM35GyMMKgukUplNdGzYhCmCRpZnOJXr168ddff+Hu7s7OnTuxs7Nj+eEb/HIpDAudlpkdymMmFxMIIYQqQkNDadOmDYmJibRt25bPP/+cJL2BIeuDeBKbRJmCDgyqV1TtmEKYJNW7l0WLFuHl5YWVlRV+fn4cPnz4lcdu27aNxo0b4+rqioODA9WrV2fv3r1ZmDb7+uCDDyhQoAA//PADnp6enLwVwZd7UoYXTAgoRZmCjionFEKYKqnTmc/NzY333nuPSpUq8e2336LVavlqz2VO336CvZUZi7r6yUVfQryCqs3sxo0bGT58OJ988glBQUHUrl2b5s2bExIS8tLjDx06ROPGjdm9ezenT5+mfv36BAQEEBQUlMXJs58GDRpw/fp1qlSpwoPIeD5Ydwa9QaF1hQJ0q1pI7XhCCBMldTpraLVapk6dytGjR7G1teWn86EsP5wye8GM9uVleIEQr6FRns+cr4KqVatSqVIlFi9enLqtZMmSBAYGpg6A/y+lS5emU6dOTJgwIV3HR0VF4ejoSGRkJA4ODm+UO7s4evQoefPmxdf3/2couBMRS7cVJwiJiKWoqy07P6iFraXM0CZyn9xUC96G1OnMtW7dOtq2bYuV1f8vUvPjufuM2HiWZINC31pefNqylIoJhVCHMXVAtTOziYmJnD59miZNmqTZ3qRJE44dO5au5zAYDERHR5Mnz6vnRU1ISCAqKirNLTe4du0arVq1olq1apw5cwaAm49j6LT0OCERsRTKY8OaPlWkkRVCvJLU6cy1fPlyunXrRoMGDUhKSgJgy+m7DNsQRLJBIbBCAcY1l+kShfgvqjWzjx8/Rq/X4+bmlma7m5sbDx48SNdzzJo1i5iYGDp27PjKY6ZNm4ajo2PqzdPT861yZwdPnz4lICCAiIgISpQoQcmSJbnyMJqOS49zPzKeoq62bBpQHc888rGVEOLVpE5nngMHDvD+++8D0Lx5c8zNzfn+j9uM2nwOgwJdqngyq2MFuTBXiHRQ/a9Eo0m7rrSiKC9se5n169czadIkNm7cSL58r55Aety4cURGRqbe7ty589aZTVlycjIdO3bk8uXLeHh48MMPP3A9IpFOS4/zKDoBX3d7Ng6ojruj1X8/mRBCIHU6o125coV27dqRnJxM165dGT9+PMsP3WD8jr8A6FOzCF+0KYtO+9/fYyGEisvZuri4oNPpXnh3HxYW9sJZgH/buHEjffv2ZfPmzTRq1Oi1x1paWmJpmXsWAhg+fDj79+/HxsaGnTt3cj/Rit6r/yA6PpnyHo58824VnGws1I4phMgGpE5nvCdPnhAQEMCTJ0+oWrUqK1asYN6v15jzyxUABtcvyqgmPul6syCESKHamVkLCwv8/PzYv39/mu379++nRo0ar3zc+vXr6d27N+vWreOdd97J7JjZysKFC1m4cCEajYZ5K9ey7bY5XZanNLKVizjzfb+q0sgKIdJN6nTGSkpKon379ly5cgVPT08+mfct7353LrWRHdWkBB819ZVGVggjqXr1z8iRI+nRowf+/v5Ur16dZcuWERISwsCBA4GUj57u3bvHt99+C6QUyJ49e/L1119TrVq11LMF1tbWODrm7nlSDQYDW7duRWNhQ9MxS5h23oJkw20AGvjmY0HXithYyMVeQgjjSJ3OONeuXePs2bM4FC6N7+A5DNl2FQAzrYaPW5Tk3VpeKicUIntStbvp1KkT4eHhTJ48mdDQUMqUKcPu3bspXDhlWdXQ0NA0cxkuXbqU5ORkBg8ezODBg1O39+rVizVr1mR1fJOi1WqZunQ9g745waVEHaBQq5gL79cvSnXvvPJOXwjxRqROZxwfH19GLf+ZFScfcSU8EUszLV2qFKJ/HW8KOlmrHU+IbEvVeWbVkNPmL0xISMDS0pJFv19j1r4r6A0Knnmsmdm+PFW986odTwiTldNqQU6S0342CQkJxBu0fLDuDIevPgagaWk3Pg8sQz57uRhXiJcxpg7I587ZWGJiIk2bNsW5ciuCdD4ABFYowOeBZbC3Mlc5nRBCiIsXL9KseQtKv7+QS0/AylzLhJal6VLFUz4xEyKDSDObTSmKwqBBgzhx6yn5qhRDA4xsXIKhDYurHU0IIQTw6NEjWrZsybMSzbj0BKzNdWwaUJ2yHrl77LAQGU31eWbFm5k9ezbfbt+Da+A4NFodgRUKMKRBMbVjCSGEIGVoQdu2bXnk6IuDXwAAcztXkEZWiEwgzWw29OOPPzL286/I1+EztFZ2VCzkxPR25eQjKyGEMAGKojBw4EDOPFLI0zhl1oePmvrQtLS7ysmEyJlkmEE2c/78eboPGoFb1+mYOeTD09maZT38sTLXqR1NCCEEMGPGDDafDMGl1Wg0Wh3t/Tx4v15RtWMJkWNJM5uNxMfH06rHAOwDJ2JmlwevvDase68arva5Y+UcIYQwdQcPHmTK+t9xaTkSgE7+nnzRtqx8ciZEJpJmNhtJUHQ4tf6Yp4lairnasP69GtLICiGECUl2LUGepinz6/aqXpiJAaXRaqWRFSIzSTObjUz84S+eJmopnNeGLYNqytK0QghhQp7EJDJ2+wUAulT2ZFKr0nJGVogsIM1sNrB69Wq0RfzZcfY+Oq2GuZ0qSCMrhBAmIi4ujkWLFnMlXx0eRSdQLJ8dE6WRFSLLSDNr4rZs2UL/oR9RsN8iNJa2DK5fjIqFnNWOJYQQgpSZC/r27cvO8w9xaemD2f9OOMhFuUJkHZmay4SdPn2aXv0G4hI4Fo2lLeU9HGUuWSGEMCFTpkxh68Gg1Cm4RjQuQZmCMpesEFlJzsyaqHv37tGqfRccWn+CZf4S2FuZMbtTBcx18v5DCCFMwaZNm5i6fDNunT5Ha2FNde+8DKwrU3AJkdWkmTVBsbGxBHToitJgGJauRXC2NuO7ftUo6mqndjQhhBDAyZMn6T9hDvk6TEJrbkmtYi4s6+mHTmYuECLLSTNrYgwGAx36DyOsXE8snNzJa2PGhgE1KO5mr3Y0IYQQwN27d2kzaiZOAWPR6Mxo6OvKwm5+Mk5WCJVIM2tidp2+yYX8zTAzsyKfjZZN79eiiIut2rGEEEIAeoPCZzv/wqx6TwBalHZlbhd/LMxkCJgQapFm1oRsO3OXUdsug5kVPs4aNnzQEGdbmYJLCCFMgcGgMHjtGfaH6AF4t7Irn7atLFNwCaEyaWZNxLkbD/hk+18YFOjo78GUwLLyTl8IIUzIkt8usefCAyzMtMzuWJ6W5QqoHUkIgTSzJuHajZsETN2F1tWbGkXzMr1tOVn+UAghTMiMZd+z4IodGjNzPm1ZShpZIUyInPpTWXR0NC0+/BqtqzeapDimBPhIIyuEECbk90OHmXM8Ao2ZOQU1T+letZDakYQQ/yDNrIr0ej2BfYaSULwBAJ+2KIG3u6zuJYQQpuLmzZt0/WItFm5F0enj2TqmjYyRFcLESDOroo9Gj+aiTVk0Wh01PCzp06Cs2pGEEEL8T1RUFO907oN52eYAzOzkh7uTjcqphBD/Js2sSlasWMHSn09hWcAHC63C171qy7t9IYQwEXq9ns6dO/O4QE00OjOqF3agjX8RtWMJIV5CmlkVXLt2jUEfDMW5bi8AhjbywdXeUuVUQgghnps7dy6/BYdg41MDLTC5bQW1IwkhXkFmM1BBsWLF6P7Zcg48ccLdwYq+tbzVjiSEEOIf3hswgJV3XIgFOlctJKswCmHC5MysCsKfJXAq1gWAj5r6YG0hSyAKIYQp+e1aJLFWLtha6BjRqITacYQQryHNbBZJTk7m448/JuzRY0ZvOU90QjKlCzjQpmJBtaMJIYQArl69ypdffsmdiBgm/3gRgPfrF5NhYEKYOBlmkEWGDx/OwoUL2XrdQIJXbSzNtExrW1bmlBVCCBPw5MkTWrZsydVbd9gQWZQnBmvKFHSgby0vtaMJIf6DnJnNAgsXLmThwoXYlW1EgldtAGZ2KE85Dyd1gwkhhCApKYkOHTpw5eo1PDtO4InBmnz2lizv6Y+VuQwDE8LUSTObyfbt28ewYcOwKlIR1xZDARjWsDgB5WUpRCGEUJuiKAwdOpRff/2VfM0GoylYFkszLct7+pPf0VrteEKIdJBhBpno0qVLdOzYEQsvP9zafoIBLe+Uy8+whsXVjiaEEAJYsGABS5YuI2+zIViXawqkfHJW3tNJ3WBCiHSTZjaThIeHExAQQFKB8uQLGIWi0dK8jDtzOlaQcbJCCGEC9uzZw/ARI3F5ZyS2peuh1cCX7crJJ2dCZDPSzGaSiIgIDK7Fca39Pmi0tKlYkBnty2Gmk5EdQghhCh49ekSeRgOwLV0PM62GOZ0qSCMrRDYkzWwmKVjYi/ytPuReZAJtKxZkZofyckZWCCFMSPGaLbC9kAeAhd0q0bS0u8qJhBBvQk4TZrD79+8DMHvfFe5FJlDQyZrJgWWkkRVCCBOQkJDA48ePiU/SM25bMABdqxaSRlaIbEya2Qz0448/4u3tzecLv2HV0ZsATGlTBjtLOQEuhBBqUxSFgQMHUqVKFT7dcJybj2Nwc7BkbHNftaMJId6CNLMZ5Pz583Tt2pWEJD0bb1pgUCCwQgHq++RTO5oQQghgxowZrFmzhvtxOrZejATg89ZlcLAyVzmZEOJtSDObAR4+fEhAQADPYmIp1fdLnpk54GxjzqctS6kdTQghBLBjxw7Gjh2Lzt6VYu/OxKBAi7LuNJHhBUJke9LMvqX4+HgCAwMJCQnBq8NYYvL4YK7TMK9LRfLayXreQgihtrNnz9K9e3c01g74DJjPM70ZRV1tmRJYVu1oQogMIM3sW1AUhb59+/LHH3/g1uQ9DF410Grg684VqV3cVe14QgiR6z148IBWrVoRm6RQrO8cYnR2FHSy5vt+Vclja6F2PCFEBpBm9i3s2LGDdevWYV+2IVYVWwEwtU1ZWpTNr3IyIYQQACNHjuTOnTsU6TyBBJt8uNhZ8F3fKrJUrRA5iFxm/xYCAwMZ+dlMdsaVIEmBoQ2K0aVKIbVjCSGE+J8FCxZw27o491zLYGGmZVXvyni72qkdSwiRgeTM7FuITdRzzs6fJEVLrWIuDGtUQu1IQggh/uF6FDxwqwrAxIBSlPNwUjeQECLDSTNrpHv37jFw4EBCw5/y0ZZzXAt7hpuDJXM7V0AnCyMIIYTqNm/ezOLFi7nyMJoP1p1Bb1BoXaEAXeWTMyFyJBlmYITY2FgC2nflpkUR9n95AL3WAp1Ww4KulXCRmQuEEEJ1J0+e5N1Rn2HtF8iXtw8BUNTVli/alEWjkRMOQuRE0symk8FgoH2/YTyuMRRHcyv0QAk3O8Y296VykTxqxxNCiFzv7t27tPlkKXm7fAmARgPNSrvzcYuS2MpKjELkWKoPM1i0aBFeXl5YWVnh5+fH4cOHX3v8wYMH8fPzw8rKCm9vb5YsWZIlOSdOnMjPG1ZhiH6Mt6OOZT382DOsDg183bLk9YUQQi3ZoU7HxMTQqlUrws79BgY9rcrmY/+IOizu7odnHptMf30hhHpUbWY3btzI8OHD+eSTTwgKCqJ27do0b96ckJCQlx5/8+ZNWrRoQe3atQkKCuLjjz9m6NChbN26NVNzrl27lilTpoBi4ONqNvw6tilNSrujlTGyQogcLjvUaYPBQM+ePQkKCsKZGLb3Kc28bpUpls8+015TCGFCFBVVqVJFGThwYJptvr6+ytixY196/OjRoxVfX9802wYMGKBUq1Yt3a8ZGRmpAEpkZGS6jj927JhiaWmpAK/MJYTIfoytBblVdqjT48aNUwDFwsJCOXr0aLpfRwhhuoypA6qdmU1MTOT06dM0adIkzfYmTZpw7Nixlz7m+PHjLxzftGlTTp06RVJS0ksfk5CQQFRUVJqbMbRaLU5OTrRp04apU6ca9VghhMjOskuddnJyQqPRsGrVKmrUqGHUY4UQ2Z9qzezjx4/R6/W4uaUdc+rm5saDBw9e+pgHDx689Pjk5GQeP3780sdMmzYNR0fH1Junp6dROatWrcrJkyf57rvv0GpVH2IshBBZJrvU6dGjR3Px4kW6detm1OOEEDmD6t3Zv6dKURTltdOnvOz4l21/bty4cURGRqbe7ty5Y3RGT09PbG1tjX6cEELkBNmhTvv6+hr9GCFEzqDaXCUuLi7odLoX3t2HhYW98K7+OXd395ceb2ZmRt68eV/6GEtLSywtZQ5YIYQwltRpIUR2oNqZWQsLC/z8/Ni/f3+a7fv373/lmKfq1au/cPy+ffvw9/fH3Nw807IKIURuJHVaCJEdqDrMYOTIkaxYsYJVq1Zx6dIlRowYQUhICAMHDgRSPnrq2bNn6vEDBw7k9u3bjBw5kkuXLrFq1SpWrlzJqFGj1PoShBAiR5M6LYQwdaouidKpUyfCw8OZPHkyoaGhlClTht27d1O4cGEAQkND08xl6OXlxe7duxkxYgQLFy6kQIECzJs3j3bt2qn1JQghRI4mdVoIYeo0yvOR+blEVFQUjo6OREZG4uDgoHYcIYRKpBaYLvnZCCGMqQOqz2YghBBCCCHEm5JmVgghhBBCZFvSzAohhBBCiGxLmlkhhBBCCJFtSTMrhBBCCCGyLWlmhRBCCCFEtqXqPLNqeD4TWVRUlMpJhBBqel4DctnshNmC1GkhhDE1Otc1s9HR0QB4enqqnEQIYQqio6NxdHRUO4b4B6nTQojn0lOjc92iCQaDgfv372Nvb49Go/nP46OiovD09OTOnTvZbvJuya4Oya4OY7MrikJ0dDQFChRAq5URV6ZE6nT2INnVkVuyG1Ojc92ZWa1Wi4eHh9GPc3BwyHa/NM9JdnVIdnUYk13OyJomqdPZi2RXR27Int4aLacjhBBCCCFEtiXNrBBCCCGEyLakmf0PlpaWTJw4EUtLS7WjGE2yq0OyqyM7ZxdvJzv/7CW7OiS7OjIre667AEwIIYQQQuQccmZWCCGEEEJkW9LMCiGEEEKIbEuaWSGEEEIIkW1JMyuEEEIIIbItaWaBRYsW4eXlhZWVFX5+fhw+fPi1xx88eBA/Pz+srKzw9vZmyZIlWZT0RcZk37ZtG40bN8bV1RUHBweqV6/O3r17szBtWsZ+3587evQoZmZmVKhQIXMDvoax2RMSEvjkk08oXLgwlpaWFC1alFWrVmVR2rSMzb527VrKly+PjY0N+fPnp0+fPoSHh2dR2hSHDh0iICCAAgUKoNFo2LFjx38+xpT+TsXbkzqtDqnTUqfTS9U6reRyGzZsUMzNzZXly5crFy9eVIYNG6bY2toqt2/ffunxN27cUGxsbJRhw4YpFy9eVJYvX66Ym5srW7ZsyeLkxmcfNmyY8uWXXyp//vmncuXKFWXcuHGKubm5cubMmSxObnz2554+fap4e3srTZo0UcqXL581Yf/lTbK3atVKqVq1qrJ//37l5s2byokTJ5SjR49mYeoUxmY/fPiwotVqla+//lq5ceOGcvjwYaV06dJKYGBglubevXu38sknnyhbt25VAGX79u2vPd6U/k7F25M6LXXaWFKnc1edzvXNbJUqVZSBAwem2ebr66uMHTv2pcePHj1a8fX1TbNtwIABSrVq1TIt46sYm/1lSpUqpXz22WcZHe0/vWn2Tp06KePHj1cmTpyoWpE0NvvPP/+sODo6KuHh4VkR77WMzT5jxgzF29s7zbZ58+YpHh4emZbxv6SnSJrS36l4e1KnpU4bS+p07qrTuXqYQWJiIqdPn6ZJkyZptjdp0oRjx4699DHHjx9/4fimTZty6tQpkpKSMi3rv71J9n8zGAxER0eTJ0+ezIj4Sm+affXq1Vy/fp2JEydmdsRXepPsO3fuxN/fn6+++oqCBQtSokQJRo0aRVxcXFZETvUm2WvUqMHdu3fZvXs3iqLw8OFDtmzZwjvvvJMVkd+YqfydircndVrqtLGkTue+Om2WkcGym8ePH6PX63Fzc0uz3c3NjQcPHrz0MQ8ePHjp8cnJyTx+/Jj8+fNnWt5/epPs/zZr1ixiYmLo2LFjZkR8pTfJfvXqVcaOHcvhw4cxM1Pv1/ZNst+4cYMjR45gZWXF9u3befz4Me+//z4RERFZOh7rTbLXqFGDtWvX0qlTJ+Lj40lOTqZVq1bMnz8/KyK/MVP5OxVvT+q01GljSZ3OfXU6V5+ZfU6j0aS5ryjKC9v+6/iXbc8KxmZ/bv369UyaNImNGzeSL1++zIr3WunNrtfr6dq1K5999hklSpTIqnivZcz33WAwoNFoWLt2LVWqVKFFixbMnj2bNWvWZPm7fjAu+8WLFxk6dCgTJkzg9OnT7Nmzh5s3bzJw4MCsiPpWTOnvVLw9qdNSp40ldTr31OlcfWbWxcUFnU73wrudsLCwF94tPOfu7v7S483MzMibN2+mZf23N8n+3MaNG+nbty+bN2+mUaNGmRnzpYzNHh0dzalTpwgKCuKDDz4AUgqPoiiYmZmxb98+GjRoYJLZAfLnz0/BggVxdHRM3VayZEkUReHu3bsUL148UzM/9ybZp02bRs2aNfnoo48AKFeuHLa2ttSuXZspU6aY7BlOU/k7FW9P6rTU6czODlKn1ZCRf6e5+syshYUFfn5+7N+/P832/fv3U6NGjZc+pnr16i8cv2/fPvz9/TE3N8+0rP/2Jtkh5Z1+7969WbdunWrjaYzN7uDgQHBwMGfPnk29DRw4EB8fH86ePUvVqlWzKvobfd9r1qzJ/fv3efbsWeq2K1euoNVq8fDwyNS8//Qm2WNjY9Fq05YJnU4H/P87aFNkKn+n4u1JnZY6bSyp07mwTht9yVgO83wKjJUrVyoXL15Uhg8frtja2iq3bt1SFEVRxo4dq/To0SP1+OdTSYwYMUK5ePGisnLlStWnfElv9nXr1ilmZmbKwoULldDQ0NTb06dPTT77v6l5layx2aOjoxUPDw+lffv2yoULF5SDBw8qxYsXV/r162fy2VevXq2YmZkpixYtUq5fv64cOXJE8ff3V6pUqZKluaOjo5WgoCAlKChIAZTZs2crQUFBqVPVmPLfqXh7UqelThtL6nTuqtO5vplVFEVZuHChUrhwYcXCwkKpVKmScvDgwdR9vXr1UurWrZvm+N9//12pWLGiYmFhoRQpUkRZvHhxFif+f8Zkr1u3rgK8cOvVq1fWB1eM/77/k5pFUlGMz37p0iWlUaNGirW1teLh4aGMHDlSiY2NzeLUKYzNPm/ePKVUqVKKtbW1kj9/fqVbt27K3bt3szTzgQMHXvu7a+p/p+LtSZ3ulfXBFanTUqfTT806rVEUEz4HLYQQQgghxGvk6jGzQgghhBAie5NmVgghhBBCZFvSzAohhBBCiGxLmlkhhBBCCJFtSTMrhBBCCCGyLWlmhRBCCCFEtiXNrBBCCCGEyLakmRVCCCGEENmWNLPC5Gg0Gnbs2AHArVu30Gg0nD17FoDff/8djUbD06dPX/n4NWvW4OTklHp/0qRJVKhQIfV+7969CQwMzPDcpqBIkSLMnTtX7RhCiBxMavSbkxqdOaSZFf/pwYMHDBkyBG9vbywtLfH09CQgIIBff/0101/b09OT0NBQypQpk+7HdOrUiStXrrxy/9dff82aNWtS79erV4/hw4e/RcoU/y7IQgiRFaRGp4/U6JzLTO0AwrTdunWLmjVr4uTkxFdffUW5cuVISkpi7969DB48mMuXL7/0cUlJSZibm7/16+t0Otzd3Y16jLW1NdbW1q/c7+jo+LaxhBDCJEiNFkLOzIr/8P7776PRaPjzzz9p3749JUqUoHTp0owcOZI//vgj9TiNRsOSJUto3bo1tra2TJkyBYAff/wRPz8/rKys8Pb25rPPPiM5OTn1cVevXqVOnTpYWVlRqlQp9u/fn+b1//0R1nNHjx6lfPnyWFlZUbVqVYKDg1P3/fsjrH/750dYvXv35uDBg3z99ddoNBo0Gg03b96kWLFizJw5M83j/vrrL7RaLdevXzfmW5jq3r17dOrUCWdnZ/LmzUvr1q25desWAHv37sXKyuqFj+aGDh1K3bp1U+8fO3aMOnXqYG1tjaenJ0OHDiUmJuaN8gghsj+p0f9PanTuJc2seKWIiAj27NnD4MGDsbW1fWH/v4vRxIkTad26NcHBwbz77rvs3buX7t27M3ToUC5evMjSpUtZs2YNU6dOBcBgMNC2bVt0Oh1//PEHS5YsYcyYMenK9tFHHzFz5kxOnjxJvnz5aNWqFUlJSUZ/jV9//TXVq1enf//+hIaGEhoaSqFChXj33XdZvXp1mmNXrVpF7dq1KVq0qNGvExsbS/369bGzs+PQoUMcOXIEOzs7mjVrRmJiIo0aNcLJyYmtW7emPkav17Np0ya6desGQHBwME2bNqVt27acP3+ejRs3cuTIET744AOj8wghsj+p0VKjxf8oQrzCiRMnFEDZtm3bfx4LKMOHD0+zrXbt2soXX3yRZtt3332n5M+fX1EURdm7d6+i0+mUO3fupO7/+eefFUDZvn27oiiKcvPmTQVQgoKCFEVRlAMHDiiAsmHDhtTHhIeHK9bW1srGjRsVRVGU1atXK46Ojqn7J06cqJQvXz71fq9evZTWrVun3q9bt64ybNiwNDnv37+v6HQ65cSJE4qiKEpiYqLi6uqqrFmz5pXfg3+/zj+tXLlS8fHxUQwGQ+q2/2vnfkKi6uIwjn+dYDCjINFkQsYiZVIiUHPIRVBRCLoIpAgcDLHBEHFcFET0hwqCaJEQlLuCYCgXIYgQEi50FiqhaCSXHMFQqU30Z+PMBN7zLoYGJ0cb5l3Efd/nA7Nwzj3Hc2fx8Lv3nHsTiYTZvn27GR4eNsYYEwqFzMmTJ1Ptw8PDxu12m69fvxpjjGltbTUdHR1p40YiEeNyuUwsFjPGGFNWVmZ6e3s3naOI/Hcoo5XRkqQ9s7IpYwyQXJ7KxpEjR9L+npqa4u3bt6mrfEheycbjcVZXV7EsC6/XS2lpaaq9vr4+q/+1/rjCwkJ8Ph+WZWXVNxsej4empiaePn2K3+9naGiIeDzOuXPnchpvamqKhYUFdu7cmfZ9PB5PLYkFAgHq6+v59OkTe/fuJRwO09jYyO7du9PGCIfDqf7GGGzbZnFxkcrKyhzPVkScSBmtjJYkFbOyqYqKCvLy8rAsK6vXpPy+zGXbNnfu3KG5uXnDsfn5+akgXi/bUM7k3/TNJBgM0traSm9vL8+ePeP8+fMUFBTkNJZt29TW1qaF3C/FxcUA+P1+Dhw4wMuXL+ns7GRgYCBtGc22bS5dukQoFNowhtfrzWleIuJcymhltCSpmJVNFRYW0tDQwOPHjwmFQhuC8Pv371tu4q+pqeHDhw+Ul5dnbK+qqmJpaSl1lQswPj6e1dwmJiZS4fDt2zfm5+c5ePBgVn1/53a7WVtb2/B9Y2MjO3bsoK+vj9evXzM2NpbT+JD8Lfr7+9mzZw+7du3a9LiWlhbC4TClpaW4XC6amprSxpibm9v09xSR/xdltDJakvQAmGzpyZMnrK2t4ff7efXqFdFoFMuyePTo0R+Xm27dusXz58+5ffs2c3NzWJZFf38/N27cAODUqVP4fD4uXLjA7OwskUiE69evZzWvu3fvMjIywvv372lra6OoqCjnl2zv27ePyclJPn78yJcvX7BtG0i+cqatrY1r165RXl6e1fJaLBZjZmYm7bOwsEAgEKCoqIgzZ84QiURYXFxkdHSUnp4eVlZWUv0DgQDT09Pcu3ePs2fPkp+fn2q7evUq4+PjdHV1MTMzQzQaZXBwkO7u7pzOW0ScTxmtjBYVs/IH+/fvZ3p6mhMnTnD58mUOHTrE6dOnGRkZoa+vb8u+DQ0NDA0N8ebNG+rq6jh69CgPHz6krKwMAJfLxcDAAIlEAr/fTzAYTNu7tZX79+/T09NDbW0tnz9/ZnBwELfbndM5XrlyhW3btlFVVUVxcTFLS0uptosXL/Lz50/a29uzGmt+fp7q6uq0TzAYpKCggLGxMbxeL83NzVRWVtLe3k4sFku7C1BRUUFdXR3v3r1LPSH7y+HDhxkdHSUajXLs2DGqq6u5efMmHo8np/MWEedTRiujBfJMpk0xIgIk35V4/PhxVlZWKCkp+dvTERGRdZTRAipmRTJKJBIsLy/T0dGBx+PJ+FCAiIj8HcpoWU/bDEQyePHiBT6fjx8/fvDgwYO/PR0REVlHGS3r6c6siIiIiDiW7syKiIiIiGOpmBURERERx1IxKyIiIiKOpWJWRERERBxLxayIiIiIOJaKWRERERFxLBWzIiIiIuJYKmZFRERExLH+AWTor5+pDvpTAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def generate_psd_matrix(n):\n", + " # generate random array of appropriate size\n", + " arr_size = int(n * (n - 1) / 2)\n", + " arr = np.random.rand(arr_size)\n", + "\n", + " # convert array to symmetric matrix\n", + " mat = np.zeros((n, n))\n", + " triu_indices = np.triu_indices(n, k=1)\n", + " mat[triu_indices] = arr\n", + " mat += mat.T\n", + "\n", + " # check if matrix is positive semidefinite\n", + " eigenvals = np.linalg.eigvalsh(mat)\n", + " if np.all(eigenvals >= 0):\n", + " return mat\n", + " else:\n", + " # if not, add identity matrix to make it PSD\n", + " mat = mat + np.eye(n) * abs(eigenvals.min()) * 2\n", + " return mat\n", + " \n", + "\n", + "def generate_correlated_samples(num_samples, num_sims, num_dims):\n", + " \"\"\" Generate samples and true parameter values \"\"\"\n", + " theta = np.random.uniform(low=-5, high=5, size=(num_sims, num_dims))\n", + " cov = [generate_psd_matrix(num_dims) for _ in range(num_sims)]\n", + " cov = np.concatenate(cov).reshape(num_sims, num_dims, num_dims)\n", + " samples = [np.random.multivariate_normal(mean=theta[i], cov=cov[i], size=num_samples) for i in range(num_sims)]\n", + " samples = np.stack(samples)\n", + " samples = samples.transpose(1, 0, 2)\n", + " theta = [np.random.multivariate_normal(mean=theta[i], cov=cov[i], size=1) for i in range(num_sims)]\n", + " theta = np.stack(theta)[:,0]\n", + " return samples, theta\n", + "\n", + "\"\"\" Main function \"\"\"\n", + "samples, theta = generate_correlated_samples(num_samples=1000, num_sims=800, num_dims=10) # You can decrease the number of simulations for faster computation\n", + "old_ecp, old_alpha = old_get_drp_coverage(samples, theta, references='random', metric='euclidean', seed = 5)\n", + "ecp, alpha = get_drp_coverage(samples, theta, references='random', metric='euclidean', norm = True, seed = 5)\n", + "\n", + "fig, axs = plt.subplots(1, 2, figsize=(8, 4))\n", + "\n", + "def plot_coverage(ecp, alpha, ax_idx, title = \"\"):\n", + " \"\"\"\n", + " Just a small utility function to plot the coverage\n", + " \"\"\"\n", + " axs[ax_idx].plot([0, 1], [0, 1], ls='--', color='k', label = \"Ideal case\")\n", + " axs[ax_idx].plot(alpha, ecp, label='DRP')\n", + " axs[ax_idx].legend()\n", + " axs[ax_idx].set_ylabel(\"Expected Coverage\")\n", + " axs[ax_idx].set_xlabel(\"Credibility Level\")\n", + " if title: \n", + " axs[ax_idx].set_title(title)\n", + " \n", + "\n", + "\n", + "plot_coverage(old_ecp, old_alpha, 0, title = \"Old\")\n", + "plot_coverage(ecp, alpha, 1, title = \"New\")\n", + "plt.subplots_adjust(wspace=0.4)\n", + "\n", + "print(\"Checking old - new\")\n", + "print(f\"Expected coverage: {np.abs((old_ecp - ecp).max())}\")\n", + "print(f\"Credibility: {np.abs((old_alpha - alpha).max())}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Test bootstrapping: \n", + "\n", + "The new `bootstrapping` function is only implemented for `get_drp_coverage`. This provides a method to offset the unpredictability of each tarp test. Bootstrapping seemed to be better suited than jackniffe for this situation because a large number of simulations is needed to have a reliable tarp test." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 800/800 [01:19<00:00, 10.09it/s]\n" + ] + } + ], + "source": [ + "# This might take some time...\n", + "ecp_mean, ecp_std, alpha_mean = bootstrapping(samples, theta, references = \"random\", metric = \"euclidean\", norm = True)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFzCAYAAADSc9khAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5JUlEQVR4nO3dd3gU1RrH8e9m0yukEEIaJXQEJQgC0iEUpYgKCtIEFUWqgCKK/aKgIIogKEWUphRFpUVUeidIi/QQSiAkIb3vzv1jkpWQTbIb0vN+nmefezM7M3smkl9Ozpx5j0ZRFAUhhBCVgkVpN0AIIUTJkdAXQohKREJfCCEqEQl9IYSoRCT0hRCiEpHQF0KISkRCXwghKhEJfSGEqEQsS7sBJU2v13Pjxg2cnJzQaDSl3RwhhLhviqKQkJBAjRo1sLDIvy9f6UL/xo0b+Pr6lnYzhBCiyF29ehUfH59896l0oe/k5ASo3xxnZ+dSbo0QQty/+Ph4fH19DfmWn0oX+tlDOs7OzhL6QogKxZQha7mRK4QQlYiEvhBCVCIS+kIIUYlUujF9UyiKQmZmJjqdrrSbIkSZo9VqsbS0lCnP5ZSE/j3S09OJiIggOTm5tJsiRJllb2+Pl5cX1tbWpd0UYSYJ/bvo9XouX76MVqulRo0aWFtbS29GiLsoikJ6ejq3b9/m8uXL1K1bt8CHgUTZIqF/l/T0dPR6Pb6+vtjb25d2c4Qok+zs7LCysuLKlSukp6dja2tb2k0SZijVX9G7du2id+/e1KhRA41Gw88//1zgMTt37iQwMBBbW1tq167N119/XeTtkp6LEPmTn5Hyq1T/yyUlJdGsWTPmz59v0v6XL1+mV69etGvXjpCQEN58803GjRvH+vXri7mlQghRMZTq8E7Pnj3p2bOnyft//fXX+Pn58fnnnwPQsGFDjhw5wqeffsqTTz5ZTK0UQoiSse30TZxsLGkT4F5sn1Gu/kbbv38/QUFBObZ1796dI0eOkJGRYfSYtLQ04uPjc7wqq44dOzJhwoQyf04hKqPgk1cYs/IoI5Yf4syNuGL7nHIV+jdv3sTT0zPHNk9PTzIzM4mKijJ6zMyZM3FxcTG8KmqFzeHDh9OvX7/SboYQwhy6NLi1kz//WMIrq0+QqYcgn1jquRVfNJer0IfcBYUURTG6Pdu0adOIi4szvK5evVrsbRRCiAIlhsGl5WwI3sqoYHcy9BY8VuM6c3u7Y2lTcLXMwipXoV+9enVu3ryZY1tkZCSWlpa4ubkZPcbGxsZQUfN+KmsmJSXl+UpNTTV535SUFJP2vV9JSUkMHToUR0dHvLy8+Oyzz3Ltk56eztSpU/H29sbBwYFWrVrx999/G96Pjo7m2WefxcfHB3t7ex544AFWr15tdls2bdpEixYtsLW1xd3dnf79+xve++GHH2jRogVOTk5Ur16dQYMGERkZaXj/zp07DB48GA8PD+zs7Khbty7Lli0zvH/9+nUGDhxI1apVcXNzo2/fvoSFhZndRiFKjC4NIrbDxW/ZcCCMSf+0Qq+xpHrSGT7vcAtLj4eL9ePLVei3bt2a4ODgHNu2b99OixYtsLKyKtbPdnR0zPN1703katWq5bnvvTeua9asaXS/+zVlyhT++usvNm7cyPbt2/n77785evRojn1GjBjB3r17WbNmDSdOnODpp5+mR48enD9/HoDU1FQCAwP57bffOHXqFC+++CJDhgzh4MGDJrfj999/p3///jz22GOEhISwY8cOWrRoYXg/PT2dDz74gH/++Yeff/6Zy5cvM3z4cMP7b7/9NmfOnGHLli2EhoaycOFC3N3Vm1zJycl06tQJR0dHdu3axZ49e3B0dKRHjx6kp6ffx3dPiGKSGglhP8D13/gx1IFJJ1qB1gquHuWHLmewqvYIWFct3jYopSghIUEJCQlRQkJCFECZM2eOEhISoly5ckVRFEV54403lCFDhhj2v3TpkmJvb69MnDhROXPmjLJkyRLFyspKWbduncmfGRcXpwBKXFxcrvdSUlKUM2fOKCkpKbneA/J89erVK8e+9vb2ee7boUOHHPu6u7sb3c9cw4YNU/r27asoivp9tba2VtasWWN4Pzo6WrGzs1PGjx+vKIqiXLhwQdFoNMr169dznKdLly7KtGnT8vycXr16Ka+99prh6w4dOhjOaUzr1q2VwYMHm3wdhw4dUgAlISFBURRF6d27tzJixAij+y5ZskSpX7++otfrDdvS0tIUOzs7Zdu2bSZ/pjBffj8rIg93TinKqZmKcnis8sOauYr/5J8V/9d/U/wHvadc2DJJUU68pyhpMYU6dX65dq9SnbJ55MgROnXqZPh60qRJAAwbNozly5cTERFBeHi44f1atWqxefNmJk6cyFdffUWNGjX44osvSmS6ZmJiYp7vabXaHF/fPTxxr3sfaimOoYiLFy+Snp5O69atDdtcXV2pX7++4etjx46hKAr16tXLcWxaWpphqEyn0/Hxxx+zdu1arl+/TlpaGmlpaTg4OJjcluPHj/PCCy/k+X5ISAjvvvsux48fJyYmBr1eD0B4eDiNGjXi5Zdf5sknn+TYsWMEBQXRr18/2rRpA8DRo0e5cOFCrtWCUlNTuXjxosltFKJY6TPg1t9wawcoCisu1GPGMV/QWqK5cog/X7KhlkM6uJdAL59SnqffsWNHw41YY5YvX55rW4cOHTh27Fgxtso4c4KuuPY1VX7f02x6vR6tVsvRo0dz/dLKHl767LPPmDt3Lp9//jkPPPAADg4OTJgwwayhEzs7uzzfS0pKIigoiKCgIH744Qc8PDwIDw+ne/fuhs/o2bMnV65c4ffff+ePP/6gS5cujBkzhk8//RS9Xk9gYCArV67MdW4PDw+T2yhEsUmLgRu/Q/RhsPVk8w0fZhyrChZaNGH7+etlJ2q6pgEe4NqiwNMVBam9UwEFBARgZWXFgQMH8PPzA9QboufOnaNDhw4APPTQQ+h0OiIjI2nXrp3R8+zevZu+ffvy3HPPAeovivPnz9OwYUOT29K0aVN27NjBiBEjcr3377//EhUVxccff2yYSnvkyJFc+3l4eDB8+HCGDx9Ou3btmDJlCp9++inNmzdn7dq1VKtWTZa+FGVP/Hm4/iskXQbHAP6OqML4v93AQoP11QP8McYJv2p2kHANvPuAjWuJNKtc3cgVpnF0dGTkyJFMmTKFHTt2cOrUKYYPH55jaKlevXoMHjyYoUOHsmHDBi5fvszhw4f55JNP2Lx5M6D+8ggODmbfvn2Ehoby0ksv5Zo9VZB33nmH1atX88477xAaGsrJkyeZNWsWAH5+flhbW/Pll19y6dIlNm3axAcffJDj+BkzZvDLL79w4cIFTp8+zW+//Wb4pTN48GDc3d3p27cvu3fv5vLly+zcuZPx48dz7dq1+/kWClF4igJRB+DyCki5AS6N2B/pwkt/uJGh1/BYrWROve2NX3UXSLgAVR+Eau1LrHkS+hXU7Nmzad++PX369KFr1648+uijBAYG5thn2bJlDB06lNdee4369evTp08fDh48aOh1v/322zRv3pzu3bvTsWNHqlevbvYDYB07duSnn35i06ZNPPjgg3Tu3Nkw+8fDw4Ply5fz008/0ahRIz7++GM+/fTTHMdbW1szbdo0mjZtSvv27dFqtaxZswZQa7rv2rULPz8/+vfvT8OGDXn++edJSUmRnr8oHfpMdTpm+E+gsQDnevx6MpVnNzmRptPQ2TeVuR3vYG1pAYmXwN4bvPuC1qbEmqhRTBkArkDi4+NxcXEhLi4uVzCkpqZy+fJlatWqJeVihciH/KwYkZkC13+D23vArjrYuPHLP8mM2+uDxtoOm5iz/DPRAVsrDaTcAl0K1BoKLg3u+6Pzy7V7yZi+EELcr/Q4uLoBYo6BYy2wcmR9SDKTDviisbbFIjKULcNRAz8zCdKiwKdfkQS+uST0hRDifqRGQvg6iAsF57qgtWXtsWSmHvJDY2WDxa3T/Pm8BTW9nEHRQcIlcG8N1YxPoChuEvpCCFFYSeFq4CeFqb12Cys2hCQy9VAtNJbWaG+e5K8XrPDzzHqWJOEiONYE755gUTrxK6EvhBCFEX9WDfy0KHBpCBotodGWTDnkpwZ+xAn+ftEa3+zAT40EC2vwfrxEHsLKi4S+EEKYK+YYXP0Z9CngXB80FoTHaxm61R2dVotLylV+vTvwM1PU0Pfuq+5fiiT0hRDCVIoCt/fCjd9AYwlOdQE4F5HCC7trcjtFS4OqGawdosUluzyyolPn47u1AM+Sm4+fFwl9IYQwhT4Tbv4BEcFg7aJOywT2n7nNwN9csXC1xNcpkxU9o3CxuWsmfPx5dUaPd2+wKN5qwKaQ0BdCiILoUuHar1lz8L0MJRP2n45k4O9uWLj6oUmN5Zu+CVSzv+uZ16RwsHYG3/4lVmahIBL6QgiRn4x4uLoRoo+oM2+s1GGbfWcieeauwF/fJ4YGHnc9WZsWpS6Y4ttfPa6MkDIMFcTw4cPRaDRoNBqsrKzw9PSkW7duLF261FCuGNRFW7L3s7Ozo0GDBsyePTtHZc6wsDDDPhqNhqpVq9K+fXt27txZGpcmROlJvQ2Xf1AD3znAEPj7T9/kmV9dswI/jnV9Ymjud1fgZyZByk2o3hWqPlRKjTdOQr8C6dGjBxEREYSFhbFlyxY6derE+PHjefzxx8nMzDTs9/777xMREUFoaCiTJ0/mzTffZPHixbnO98cffxAREcHOnTtxdnamV69eXL58uSQvSYjSk3wNLn8P8f+qc/C1apnwA6dvMuCXKli4+RsCP/DuwNdnPYDl1gqqd4Y81u8uLRL6FYiNjQ3Vq1fH29ub5s2b8+abb/LLL7+wZcuWHGsTZK9HW7NmTUaNGkXTpk3Zvn17rvO5ublRvXp1mjZtyqJFi0hOTja6nxAVTsJFtYefHK7Owc+6AZucoWHyQT+01eqgSU/kpz53CPSzvufYC+BUB7wfKxM3bu8lY/oFUBSFlAxdqXy2nZUWzX32Ejp37kyzZs3YsGEDo0aNyvGeoijs3LmT0NBQ6tatm+957O3tAcjIyLiv9ghR5sWehqvrISMOnBuo1TKB1Ex4IdiVa3pbbDXpLOsVTQu/e0I9+QZY2oNPH7CuUvJtN4GEfgFSMnQ0mrGtVD77zPvdsbe+//9EDRo04MSJE4avX3/9dd566y3S09PJyMjA1taWcePG5Xl8UlIS06ZNQ6vVGhZhEaJCijmm3rTVp4NTPcPQTMi/N5hzsSF7I2xxsNLzfc9Ymle752czI179ReH3NDjWLoXGm0ZCvxJQFCXHXwxTpkxh+PDh3L59m+nTp9O5c2fDurN3a9OmDRYWFiQnJ+Pl5cXy5ct54IEHSrLpQpScqANw7RfQaNXhmSxHToXTd5UFVrUdsdXqWRIUTfNq9/zFq0uFxCvqYijuj5Rww80joV8AOystZ97vXmqfXRRCQ0OpVauW4Wt3d3cCAgIICAhg/fr1BAQE8Mgjj9C1a9ccx61du5ZGjRpRpUoVw2LpQlQ4igK3d8O138DSTl3YJMvhU+H0+z4Tq7rtQJ/B3A6RPOJ1z/H6TEg4D64PquP4mrJ9q1RCvwAajaZIhlhKy59//snJkyeZOHGi0ferVq3K2LFjmTx5MiEhITn+IvD19aVOnTpGjxOiQlD0cPMviNgMVlXAztPw1uGTV+i7PAXrhl1Ar+OL9rfpee+tL0WvFl5zqge+T6rj+WVc2f6VJMySlpbGzZs3uX79OseOHeN///sfffv25fHHH2fo0KF5HjdmzBjOnj3L+vXrS7C1QpQyRa+WVLjxO1i75gz8E2H0+SZWDXxFz6ePRtKnwT2LDCqK2sO391bH8cvojdt7ld8urMhl69ateHl5YWlpSdWqVWnWrBlffPEFw4YNy7Eo+r08PDwYMmQI7777Lv379y/BFgtRSvQ6iNiqhr6dJ9j8N3x55GQYvRdEYNviKQA+euQ2TzXS5z5H0hX1YS2/pwx1eMoDWSP3LrLupxCmKdc/K/oMuLEFbu4A+xq5att/+LfCtxd8AJgeeJsXHkrPfY7U2+psnZrPQtUHS6DR+ZM1coUQwpjMZLi+CW7vA3ufXEMyP4Ta8+0F9ZfAhGbRxgM/IyGrNn6fMhH45pLQF0JUDulx6kNXMSHgVBssHQxv7Q+5yOZrbnx/S52583KzBCY8nJr7HLp0SAxTp2Z6diyZdhcxuZErhKj4Um9D2Eo18J3r5gj8fccu0vuTf1hxsyEAIxonMrVFfO5zKDpIOAdVm2aVWCiaKdUlTXr6QoiKLSUCrqyFxEuGxcuz7T16gT6fHMOx23g0GgueqRvHjEcSc9dIUxSIOwsONcHniXIxNTMvEvpCiIor6aoa+MlXDYuXZ9tz5AJ9Zx7CsfskNBoLng6I53/tjQQ+QOJFsPVQp2baupdc+4uBDO8IISqmxMvqkE7K9VyBv/vweXq/twvHoAloLLT0q53AJx0SsDAW+MnXwMJGXQzFwbfk2l9MJPSFEBVP3L8Qtkody3eun6M0wunzN+g9YwdOvaag0VrR3S+RzzrGGw/81NuQmaKub+vSsOTaX4xkeEcIUbHcOa5WytSlgHO9XIuYpDv54tbvLXRaGzp4J/Nllzi0xrq/GfGQdhtqPA5uD5dI00uChL4QomJQFIjaD9d/VXv2TgG5djkVZcXwYHd0WgtaeqawqNsdrI1NwtGlqFUzq3cCz05lbvWr+yGhL4Qo/xQFInepgW/lCHY1crz914GzfPvHVU75PUdcugWB1dJY2v0OtsYSUJ8B8RfALRBq9Cq3UzPzImP6Qojy7/ZeNfCtXXIF/p/7/6Xvm7+xy6U/celamldLY3mPaBytjVSgya6a6VxfnZqpLWclJkwgPX1TZSSqCyWUFK2t2mMRdOzYkQcffJDPP/+8tJtSpkRHR9OwYUMOHTpEzZo1i+y8Tz31FG3atGHSpElFds5idXsfXNsEVs5g65njrT/2hvLEm79S5akP0No50cw9le96xOBkNPAViD8H9n5qETVrlxK6gJIloW+KjES4sBjSo0vuM63dIOBFs4I/r3D8+eefeeKJJyjt2nqFDe8NGzZgZWXeAtOV4RfFzJkz6d27tyHwFy5cyMKFCwkLCwOgcePGzJgxg549e5p13hkzZtCpUydGjRpVYPGuUhd1UF3tyso+V6XL4D1neOL1DVR5+iO0ds40c0/j+155BD5A4gW12qb/0+Wqaqa5JPRNoUtVA19rB9oSeBJPl6x+ni61wvT209ONFK4ykauraxG2pGJISUlhyZIlbN682bDNx8eHjz/+mIAA9Qbmd999R9++fQkJCaFx48Ymn7tp06bUrFmTlStX8vLLLxd524tM1CG4tlH9q/ieIZ1tu0/zxOS1uA74H5aOrtSvms6KntE45xn4YerPtt9T4OBf/G0vRTKmbw6tvRrCxf0q5l8sHTt2ZNy4cUydOhVXV1eqV6/Ou+++a3hfr9fzySefEBAQgI2NDX5+fnz00UeG9xVFYdasWdSuXRs7OzuaNWvGunXrcn3Gq6++yqRJk3B3d6du3brs3LmTefPmodFo0Gg0hh7p1q1befTRRw3LMj7++ONcvHgxx7kmTJhgcvuHDx+e67Pef/993NzcSEtLy9HOJ598Ms8FZsLCwgzH3/3q2LEjoC5aM27cOKpVq4atrS2PPvoohw8fNut7ber3815btmzB0tKS1q1bG7b17t2bXr16Ua9ePerVq8dHH32Eo6MjBw4cyHHsgQMH6NKlC+7u7rmuLTY2FoA+ffqwevXqfNtQqrID38Imx/KGAHfikhg4dTVVnngHSxdP/J0y+KFnNC42eQR+yg1AAd8n1LH8Ck5Cv5L67rvvcHBw4ODBg8yaNYv333+f4OBgAKZNm8Ynn3zC22+/zZkzZ1i1ahWenv+Nlb711lssW7aMhQsXcvr0aSZOnMhzzz3Hzp07c32GpaUle/fuZdOmTbRu3ZoXXniBiIgIIiIi8PVVn25MSkpi0qRJHD58mB07dmBhYcETTzyBXm9k4QoT2j9v3rxcn/Xaa6+h0+nYtGmT4RxRUVH89ttvjBgxwuhn+Pr6Go6PiIggJCQENzc32rdvD8DUqVNZv3493333HceOHSMgIIDu3bsTExNjclvN+X7ebdeuXbRo0SLP93U6HWvWrCEpKSnHL4Z//vmHjh070qxZM3bt2sXWrVtxdXWlU6dOrF27lipVqgDQsmVLDh06lOuXZJkQfTgr8K1yBT6Aha0jD7z8OdbuflS3z2Rlr2g87PP4t5QWrZZK9u5dLsskF4YM71RSTZs25Z133gGgbt26zJ8/nx07dvDII48wb9485s+fz7BhwwCoU6cOjz76KKAG9Jw5c/jzzz8NYVK7dm327NnDokWL6NChg+EzAgICmDVrluFra2tr7O3tqV4953jpk08+mePrJUuWUK1aNc6cOUOTJk3Man+3bt1wcXEx+lmDBg1i2bJlPP300wCsXLkSHx8fQ8/9Xlqt1nB8amoq/fr1o3Xr1rz77rskJSWxcOFCli9fbhgz/+abbwgODmbJkiVMmTLFpLaa8/28W1hYGDVq1Mi1/eTJk7Ru3ZrU1FQcHR3ZuHEjjRo1Mrw/btw4+vbty5w5cwBo1KgRzz77LAcPHmTAgAGG/by9vQ3Lb/r7l6HhDkMP30qth3+XjAwdyXpLhmxx52q6NW62On7oFY2Pk874uTISIOUm1OgJ7q2N71MBSehXUk2bNs3xtZeXF5GRkYSGhpKWlkaXLl2MHnfmzBlSU1Pp1q1bju3p6ek89NBDObbl1xO928WLF3n77bc5cOAAUVFRhh5+eHh4vqFvrP35eeGFF3j44Ye5fv063t7eLFu2jOHDh+dYDD4vI0eOJCEhgeDgYCwsLLh48SIZGRm0bdvWsI+VlRUtW7YkNDTU5Laa8/28W0pKitEVq+rXr8/x48eJjY1l/fr1DBs2jJ07d9KoUSNu3brFnj17+PPPP3Mc4+DgkOt7YGdnB0BycnKebShRivLftEytTa7A/+3PE0yY9St1R84mNM4aV1sdq3pFEVAl0/j5dCmQFAbVOkD1rhXq4auCSOhXIM7OzsTFxeXaHhsbm2sWxr2zYTQaDXq93vDDnpfsQP7999/x9s75p7WNjU2Orx0cHDBF79698fX15ZtvvqFGjRro9XqaNGmS783fvNqfn4ceeohmzZqxYsUKunfvzsmTJ/n1118LbN+HH37I1q1bOXToEE5OTgCGmVD3hqWiKLm25ddWc76fd3N3d+fOnTu5tltbWxtu5LZo0YLDhw8zb948Fi1axNGjR9Hr9TRr1izHMUePHs31Czp7iMrDwyPPNpQYRQ+3/lYXMLdyyTWz5tcd//DUpO+o+sQMQuMcqGKj54eeUdR3zSPw7374yvtxsKhcMShj+hVIgwYNOHLkSK7thw8fpn59025Q1a1bFzs7O3bs2GH0/UaNGmFjY0N4eDgBAQE5Xtlj9HmxtrZGp8v5p3Z0dDShoaG89dZbdOnShYYNGxoNM3MZ+yyAUaNGsWzZMpYuXUrXrl0LbPP69et5//33+fHHH6lTp45he0BAANbW1uzZs8ewLSMjgyNHjtCwoemFuQr7/XzooYc4c+ZMgedXFMUwLp/9CyYlJcXw/smTJ9m1axfPPfdcjuNOnTqFj48P7u6lXEZYr4OI7XDjN7B2zRX4m3b8w5MTllCl9zRsfRrjbK0GfiO3PAJf0akPX7k0VKtmVsCHrwpSuX7FVXCvvPIK8+fPZ8yYMbz44ovY2dkZxpi///57k85ha2vL66+/ztSpU7G2tqZt27bcvn2b06dPM3LkSJycnJg8eTITJ05Er9fz6KOPEh8fz759+3B0dDTcBzCmZs2aHDx4kLCwMBwdHXF1daVq1aq4ubmxePFivLy8CA8P54033rjv74Wxz7KwsGDw4MFMnjyZb775hhUrVuR7jlOnTjF06FBef/11GjduzM2bNwH1F4qrqysvv/wyU6ZMwdXVFT8/P2bNmkVycjIjR440uZ2F/X52796dadOmcefOHapWVdd0ffPNN+nZsye+vr4kJCSwZs0a/v77b7Zu3QpAq1atsLOzY+rUqUyfPp2LFy8yduxYRo8eTZs2bXKcf/fu3QQFBZl8HcVCnwk3tsKtP8C2OtjknLr7yx/HeXr8t1TpNQW7mg9ib6nnux5RNHHPMH6+7IevHGqqdfGtyvgzCMWk1Hv6CxYsoFatWtja2hIYGMju3bvz3X/lypU0a9YMe3t7vLy8GDFiBNHRJfTQlC5ZfVCruF+6wo2j1qxZk927d3Px4kWCgoJ4+OGHWb58OcuXLzfcvDTF22+/zWuvvcaMGTNo2LAhAwcOzDFe/sEHHzBjxgxmzpxJw4YN6d69O7/++iu1atXK97yTJ09Gq9XSqFEjPDw8CA8Px8LCgjVr1nD06FGaNGnCxIkTmT17dqGuv6DPAnUI7Mknn8TR0ZF+/frle44jR46QnJzMhx9+iJeXl+HVv39/AD7++GOefPJJhgwZQvPmzblw4QLbtm0zhLCpCvP9fOCBB2jRogU//vijYdutW7cYMmQI9evXp0uXLhw8eJCtW7ca7hd4eHjw448/cujQIZo2bcq4ceMYPXp0rgfYUlNT2bhxIy+88IJZ11Gk9Jlw/Xe4GQx2XrkC/+fg4zw1djHOXcdgX78N1hYK3wZF81C1PAIf1IVQbNzVufi2ZWDYqpRolFJ8THPt2rUMGTKEBQsW0LZtWxYtWsS3337LmTNn8PPzy7X/nj176NChA3PnzqV3795cv36d0aNHU7duXTZu3GjSZ8bHx+Pi4kJcXFyuce7U1FQuX75s+CVkUE6eyBWm6datGw0bNuSLL74o7abcl82bNzN58mROnTqFhUXR9d+++uorfvnlF7Zv357nPnn+rBQFXbo6nHPrb7D3zVUOQa/X0+rpT7joEYRziz5oNQoLu8YQ5J9PmZTkq6AANQerSyZWMPnl2r1KdXhnzpw5jBw5klGjRgHw+eefs23bNhYuXMjMmTNz7X/gwAFq1qzJuHHjAKhVqxYvvfRSjmmBxcLKUQ1gqb1TrsXExLB9+3b+/PNP5s+fX9rNuW+9evXi/PnzXL9+vcB7E+awsrLiyy+/LLLzmUWXptbRub0bHPyMDsFYWFjQe8J7LD3rBsCn7e/kH/hpUepCKP4DK2Tgm6vUQj89PZ2jR4/mGr8NCgpi3759Ro9p06YN06dPZ/PmzfTs2ZPIyEjWrVvHY489lufnpKWl5XjAJD7eyCr3psh+WlaUW82bN+fOnTt88sknJt/YLuvGjx9f5Od88cUXi/ycJtGlqoufRO1Xx93v+Xm7fDWKWr7uLPzHkaVn1d7/B21ieaJuipGTZclIhJRb4P0YuJo2hbiiK7Ux/aioKHQ6XY4nPQE8PT0NN8zu1aZNG1auXMnAgQOxtramevXqVKlSJd9eycyZM3FxcTG8irJHJMqXsLAw4uLimDx5cmk3RdwrMxnC10HUPnCslSvwf9x8hHpBb/Piigg+OawG/hsPxzGkUVLe59SlQeIl8GgLnp0r1Vz8/JT6jVxT5jlnO3PmDOPGjWPGjBkcPXqUrVu3cvnyZUaPHp3n+adNm0ZcXJzhdfXq1SJtvxDiPmUkwpUf1YqZjgFgmfP5jrW/H2bQpCXYNOjI9nS1tz72wXhGN0vM+5x6HSScg6rNKuVc/PyU2nfC3d0drVabq1cfGRmZq/efbebMmbRt29bwiHvTpk1xcHCgXbt2hhkW97Kxscn3IRchRCnKiIfwnyAmRF3P9p5586t/PcRzk5diW7cN7r0mADCicSKTAhPyPqeih4R/wbGOOhffsgQq45YjpdbTt7a2JjAwMEfhKYDg4OBcc4azJScn55qloNWqS5mVdq14IYSZ0mMhbHVW4NfPFfirNqmBb10zkGp9p4LGgoH1kpjxSFz+IzUJF8DWS52Lb+NWrJdQHpXq8M6kSZP49ttvWbp0KaGhoUycOJHw8HDDcM20adNylL3t3bs3GzZsYOHChVy6dIm9e/cybtw4WrZsabT4lBCijEqLhrBVEHdKnVGjzfnX+A+/HGDIlKVYeTfG68npKBotj9dO5n+PxuYf+IlhYOmozsU3UoFTlPKUzYEDBxIdHc37779PREQETZo0YfPmzYaqfhEREYaHakCtk56QkMD8+fN57bXXqFKlCp07d+aTTz4p0nbJXw1C5O++fkaSb0D4j+pNVucGasXMe1y4chvL6vXwfuY9dBZWdPFNYW7HO2jz66amRKDWxe8HTgGFb18FV6oPZ5WG/B5i0Ol0nDt3jmrVquHmJn8WCpGX6OhoIiMjqVevnmGI1SQJF9Ux/NSb4FQfLIwfezrKkqc2uZKit6JtjVSWBEVjm18XNS0K0mLUhVA82uazY8VUbh7OKmu0Wi1VqlQxlBywt7c3qeyuEJWFoigkJycTGRlJlSpVzAv82FNwdb1689a5AWhydts3/32SDi3rEZHmwJCt7qTotbTwTOObbjH5B356LKREqnPx3Y3fDxT/kdC/R/aiGQXVZheiMqtSpUquxXDyFXtSnZapZBpdknDZur2MfPN72nZoTWrHqcSkannAPZ2l3aOxt8pnMCIjEZKvqTXxZS6+SST076HRaPDy8qJatWpkZORTvEmISsrKysq8Hn7cGfXBKyVTffDqHkt+2sML03/AwtGVqMBXSEmxpH7VDFb0yGchc8haCOUyeDyqrn6Vx1CRyElCPw9arda8f9hCiNziz6k3bXWp4FQn19vf/riHF6Z/j4W9C/VfmEOSZRVqOmfyfc8oqtrmsyiOLg3iz4N7K/DuY/RmsDCu1J/IFUJUUAkX1CGdzCSjgb94zS418G0caDx6LklWbtRwyOSHnlFUy2shc8ha+eocuD6U9fBV/qu9iZwk9IUQRS/+PIStgYxYtbTCPZb+tJeX3l6JxtqOZq/MJd6qGu52Olb2isp7IXNQ6+zHn4UqjcH3qVwlG0TBZHhHCFG04s/BlTWQkQBOdY3eXH2wkS9VXV2oPexjoqxr4GKj5/seUdRyySfwFb0a+E4B4DcgV519YRrp6Qshik7cv2rgZyao4ZzHbJrGDfzo9ubXRFn74milZ0WPKBrmta5ttoQLYF8jq7yCa/77ijxJ6AshikZcKISvgcxEdUjnnsBftHoX+45dJEMPY/905eBtJ+ws9SzrHk0zjwJmyiWFq0M5vk+qyyeKQpPhHSHE/YsLhfC1al18IyUQ5n//F2PfX4OToz1PfrCIv27aYa1V+LZbDA9XT8//3Km31Ju3/k9KeYUiUOiefvYi0Ckp6qo1layagxAiW9wZuLJWXZLQSCh/8d2fjH1/DWgsaD7qff66WRUrC4Wvu0TT1jvNyAnvkn4H0u6o8/BdHyqmC6hczA796OhounbtSr169ejVqxcREREAjBo1itdee63IGyiEKMNiT6vTMnUpRqdlfr7sD8Z/uBbQ0GnCx4TZNMBSozC/cwyd/QoI/Ix4tThb9S5QrX3xtL8SMjv0J06ciKWlJeHh4djb/7c4wcCBA9m6dWuRNk4IUYbdOa7etNUbf/Bq7rI/mPi/nwDoOvF/XLJuhIVGYV6nGLrXzGchc1CHiRKvQLV24NU9V50eUXhmj+lv376dbdu24ePjk2N73bp1uXLlSpE1TAhRhkUdhGs/g0YLjrVzvb1h2zEmZQV+0MT/cdb6ATQofNb+Do/VLiDwdWlqNU73R8C7tyx1WMTM/m4mJSXl6OFni4qKkmUJhajoFAVu74Hrv4LWLs+FSnp1fIAeHZqgCXyWM1ZNAfikfSxP1E3J//z6DHVtW9cH1br4WsmUomb230zt27dnxYoVhq81Gg16vZ7Zs2fTqVOnIm2cEKIMUfRw80+49ou6OpWRwM+e0GFjbUW7UW9yxqoZADMfvcOAesn5n1+vUx++cm4oT9sWI7N7+rNnz6Zjx44cOXKE9PR0pk6dyunTp4mJiWHv3r3F0UYhRGnT6+BmMERsV9edtfXItcvHi7YSGR3PZ9OeZu4xZxafUhfz+KBNLM82KCDwFT0knFWHiuRp22Jldug3atSIEydOsHDhQrRaLUlJSfTv358xY8bg5SUPTQhR4egz4cYWuPmH+mCUkadhZ369hTc/+1n9ovHjbLit3vOb8UgsQxol5X9+RVFr9dh5qYFv617EFyDuJsslCiHyps9Qx+9v7VSHc6yr5NrlowWbeWvuLwAMfG0aByzV5QqnPhzHK80SC/6MhItg6QS1Bhutty8KVqzLJZ44ccLodo1Gg62tLX5+fnJDV4iKQJcO1zdB5C5w8AOr3GHy4Ve/8/bnmwAY8tpUdmUF/riH4k0L/OSrai18v/4S+CXE7NB/8MEHDevGZv+RcPc6slZWVgwcOJBFixZha2tbRM0UQpQoXbp6w/b2bnDwByunXLu8/+VvvPPFrwAMmzKVnRbtAHjhgQQmNk8o+DNSb6nTM/0GgEujIm2+yJvZs3c2btxI3bp1Wbx4Mf/88w/Hjx9n8eLF1K9fn1WrVrFkyRL+/PNP3nrrreJorxCiuOnS4NpGuL0LHGoaDfyTZ6/z3vzfAHh+ymR2WrRDQcPQRom82TK+4KVqDeUVeoFbi6K/BpEns8f0W7ZsyQcffED37t1zbN+2bRtvv/02hw4d4ueff+a1117j4sWLRdrYoiBj+kLkQ5cKVzfC7X3qcIuVY567/vDLAXZGurHDogN6RcPgBol82Dau4MDPiIekq+AVpNbUkadt75s5uWb2d/vkyZP4+/vn2u7v78/JkycBdQgouyaPEKKcuDvwnWrnCnxFUYhP+O/hKtemHfkzK/CfqZ/EB6YEfmYSJF1Ra+lIeYVSYfZ3vEGDBnz88cekp/9XDjUjI4OPP/6YBg0aAHD9+nU8PT2LrpVCiOKlS4XwDRCVFfj3PBilKApvzf2Flk/N5ObtOHaE2zLuT1d0ioan6ibxv0djsSgo8HWpkHAJ3FuD9+NSXqGUmP1d/+qrr+jTpw8+Pj40bdoUjUbDiRMn0Ol0/PabOsZ36dIlXnnllSJvrBCiGGSmwNUNEH1AfTjKSOBPn/MzM79WCyp+teM2a5Mbkqlo6FMnmU/amRL46epcfNfm4NNPyiuUokLN009MTOSHH37g3LlzKIpCgwYNGDRoEE5OuW/4lDUypi/EXTKT4ep6iDqkVsq0zFlXS1EUpn26kU8WbwNgyptj2UQ3UnUWdPdPYX6XGKwKGi/QZ0D8v+DSBGo+a3Tqp7g/xTpPH8DR0ZHRo0cXqnFCiDIiMwnC10P0YXXxE0u7HG8risIbszcw65vtALz51hh+oRupGRZ09Enli84mBn7cv1ClMfgPlMAvAwo9qHbmzBnCw8NzjO0D9OnT574bJYQoZhmJcHUdRB8F5wC1YuZdFEVh6ifr+XRJMABvv/0KG+lOYroFrb3S+LprNDbaAj5Dn5lVQK1BVj2dKsVzLcIsZof+pUuXeOKJJzh58iQajSbXA1o6na5oWyiEKFrpcWrgxxw3GvgAcQkp/PzHPwC8O+NlNtKDhDQLHvZMY0lQNLYFJYeiU4d0nALUHr6Rej2idJg9e2f8+PHUqlWLW7duYW9vz+nTp9m1axctWrTg77//LoYmCiGKTFoMXFmdFfj1jAY+QBVne/76fhKffTyOTZoexKZZ8KBHOku7R2NvVcBtQEUP8efUB7v8n5ECamWM2aG/f/9+3n//fTw8PLCwsMDCwoJHH32UmTNnMm7cuOJooxCiKKRGQthKdSFzlwa5ZtAoisKx0+H/7W7rwarULkSnamnils53PaJwsi4o8BVIOA921cF/ANhWK44rEffB7NDX6XQ4OqoPbbi7u3Pjxg1AfTjr7NmzRds6IUTRSI1S17NNuKCOsVtY5XhbURQmfPgjLZ+cyY+bj3Ah1pJnNrsTlaKloWs63/eMxsXGhIl+iRfBqir4PQ32PgXvL0qc2WP6TZo04cSJE9SuXZtWrVoxa9YsrK2tWbx4MbVr514rUwhRytJiIHxtVuA3BIucd2AVRWHcB2uY//3fAFyOs2L272rgN3DNYGWvaKra6gv+nMTL6nCR31NG180VZYPZof/WW2+RlKQuivDhhx/y+OOP065dO9zc3Fi7dm2RN1AIcR/S4yD8R3XapEsDo4H/6nurWbByJxqNhk/ef5mf9N25naKlQdUMVvWKwtWUwE8KB42lGvguDYrpYkRRKJJFVGJiYqhatWqOEstllTycJSqNjAQ18GOOg0t9sLDO8bZer+fV99awcJUa+LM+eJmfdN25laylflbgu9mZEPjJ19X5+H5Pg+tDxXMtIl/FVnAtMzMTS0tLTp06lWO7q6truQh8ISqN7MC/c1ydpWMk8F95Z7Uh8Gd/+ErhAj/lplpTx6efBH45YVboW1pa4u/vL3PxhSjL7g58p3pG69xoNBq0Wg0ajYbP/vcKP2UGcStZS90qGaw0NfBTb6tlkr17g3vLor8OUSzMnr3z1ltvMW3aNGJiYoqjPUKI+2FC4IMa+vPfeZafV7zNTxnduJkV+Kt6ReFuSuCn31FD36sHeLQt2msQxcrsG7lffPEFFy5coEaNGvj7++PgkLMi37Fjx4qscUIIMxQQ+Hq9nkWrdzPy6bZYW1tyPdGST8KacyPJktouag/fw96EwM+IV8fxqwdB9c4UXERflCVmh36/fv2KoRlCiPuSkQjhP+Ub+C9M/4Gl6/ayY/+/zJv5MoM2u3M90ZJazpms7hVFNVMCPzMJEq+AZweo0UMWQSmHzA79d955pzjaIYQorOzAjzmWVVohZ+DrdHpGvbmC5Rv2Y2GhoVvX1gze4kF4giV+Tpmseuw2ng4mBL4uRV0ExaO1Oo4vi6CUS4X6NR0bG8u3336bY2z/2LFjXL9+vUgbJ4QoQGYShK+7K/Btc7yt0+kZOU0NfK3WgkWzX2FtRlfC4i3xccxk9WNReJkU+KkQf0G9YevzhCyCUo6Z/av6xIkTdO3aFRcXF8LCwnjhhRdwdXVl48aNXLlyhRUrVhRHO4UQ98pIzAr8I3kG/og3vuP7nw+g1VqwePbLrEnvyqU4K2o4qIHv7WjCTDxdulpAzbU5+PbPVXdflC9m9/QnTZrE8OHDOX/+PLa2//0j69mzJ7t27SrSxgkh8pCRoJZWyCPwAca8u9oQ+Es+e5mfMrpw7o4VnvY6VvWKwtfJhMDXp0PCWajygPrw1T1LKYryx+zQP3z4MC+99FKu7d7e3ty8edPsBixYsIBatWpha2tLYGAgu3fvznf/tLQ0pk+fjr+/PzY2NtSpU4elS5ea/blClFsZ8TkfvDIS+ACDej+Mi5Mdy+a8wk/pXQiNscbdTg38mi6mBH4GxJ1Vlzn0HwhWZX85VFEws4d3bG1tiY+Pz7X97NmzeHh4mHWutWvXMmHCBBYsWEDbtm1ZtGgRPXv25MyZM/j5+Rk9ZsCAAdy6dYslS5YQEBBAZGQkmZmZ5l6GEOVTRjxcWQt3Thi9aXu39i3rcWL7x7y6x5dT0da42epY3SuKOlVM+HnRZ6irXrk0VEsky6pXFYbZtXdefPFFbt++zY8//oirqysnTpxAq9XSr18/2rdvz+eff27yuVq1akXz5s1ZuHChYVvDhg3p168fM2fOzLX/1q1beeaZZ7h06RKuroVbiUdq74hyKyM+a5bOcXCuD9qcpRUyM3WM+2Ato59tT9MGPiSkaxiyxZ3jt61xzQr8+q6mBH5m1qpX9dWFzGXVqzKv2GrvAHz66afcvn2batWqkZKSQocOHQgICMDJyYmPPvrI5POkp6dz9OhRgoKCcmwPCgpi3759Ro/ZtGkTLVq0YNasWXh7e1OvXj0mT55MSkpKnp+TlpZGfHx8jpcQ5U5GQoGBP/i1JSxctZOeo74kKiGD4VvdOH7bmio2en7oaWrgZy9zWFeWOaygzB7ecXZ2Zs+ePfz5558cO3YMvV5P8+bN6dq1q1nniYqKQqfT4enpmWO7p6dnnvcGLl26xJ49e7C1tWXjxo1ERUXxyiuvEBMTk+e4/syZM3nvvffMapsQZUpG4j1j+DkDPyNDDfyfthzFykrL3HeG8NJfXhyNtMHZWg38Rm4mBL6ig/jQrHVtZZnDisrs0A8LC6NmzZp07tyZzp0733cD7q3OqShKnhU79Xo9Go2GlStX4uLiAsCcOXN46qmn+Oqrr7Czyz2VbNq0aUyaNMnwdXx8PL6+vvfdbiFKRGayOi0zjydtMzJ0DJr0Leu2HsPKSsvKea+wNrWjIfBX9oqiiXtGwZ+j6NSa+461swLfvPtzovwwe3indu3aPProoyxatOi+iq65u7uj1Wpz9eojIyNz9f6zeXl54e3tbQh8UO8BKIrCtWvXjB5jY2ODs7NzjpcQ5YIuFa5ugJij6nCLkcB/ZsI3rNt6DGsrS1Z/MYYfUzty5NZ/PfwHTAp8vXrT1sE/K/BlXduKzOzQP3LkCK1bt+bDDz+kRo0a9O3bl59++om0tDSzzmNtbU1gYCDBwcE5tgcHB9OmTRujx7Rt25YbN26QmJho2Hbu3DksLCzw8ZH1OEUFokuHaz9D9EFwqmN0WuZHCzezYXsI1laWrJk/hrWpHTh8ywYnaz3f94yiqYcZgW/nrQa+XfWivxZRppgd+s2bN2f27NmEh4ezZcsWqlWrxksvvUS1atV4/vnnzTrXpEmT+Pbbb1m6dCmhoaFMnDiR8PBwRo8eDahDM0OHDjXsP2jQINzc3BgxYgRnzpxh165dTJkyheeff97o0I4Q5ZIuHa5vgtv71OEWS3uju732fDe6tGnAj1+NYV1qOw7dtMHJSs/3PaJoZlLgK+q6ubbV1Ju29jWK+EJEWVQkyyUeO3aMkSNHcuLECbMXWFmwYAGzZs0iIiKCJk2aMHfuXNq3bw/A8OHDCQsL4++//zbs/++//zJ27Fj27t2Lm5sbAwYM4MMPPzQ59GXKpijT9Blw7ReI3KUOt9zzQJROp0er/a+vlpqp8NIf7uy8Zou9pZ7ve0YT6Jlu2mclXABLR6g5WP1rQpRb5uRaoUP/6tWrrF69mlWrVnHy5Elat27N4MGDefnllwvV6JIioS/KLH0GXNsEkTuNBn5aWgZPj1tMq2a1mP5KL9J18MoOV/4It8NWq+e7HtG08jIx8BPDwMJKHdJxaVj01yJKlDm5ZvbsncWLF7Ny5Ur27t1L/fr1GTx4MD///DM1a9YsbHuFEPoMuP5rVuD7GQ38p8Yu4re/ThK8N5SnH2/FrH/r8ke4HTZahSVBMaYHfnLWpAffJyXwKyGze/q+vr4888wzDB48mAcffLCYmlV8pKcvypzsMfzIXVmBn/PfZWpaBk+O+ZrNO09ha2PFhoWvsi65DTuuqoH/Tbdo2vuYOJEi5aZajtnvKXB7uBguRpSGYu3ph4eH5zmPXghhJl2aOqRzexc41MzVw09Ny6D/mK/ZsvMUdrZWbFg4lrVJbfjrqm1WDz+aR71NDPzU25ARp9bDl8CvtMwOfY1GQ2xsLEuWLCE0NBSNRkPDhg0ZOXJkjvnzQogC6FLVm7a394BDLbByzPF2aloGT7yykK27TmNna8XPX49lTaIa+LZaPUuCYmhrcuDfgrRY8H5MFjKv5Ao1T79OnTrMnTuXmJgYoqKimDt3LnXq1JFF0YUwlS4Vrm6EyD3qtMx7Ah9g667ThsD/ZfE41iW35s+swF/aPdr0wE+5Celx4NMbPDvJQuaVnNlj+u3atSMgIIBvvvkGS0v1D4XMzExGjRrFpUuXyvxCKjKmL0pd9pO2Ufuz5uHnvTDJVz/8RYMAbzYkteL3y/ZYaxWWdIumnalj+MnX1c/z6QPubSTwK6hinbJpZ2dHSEgIDRo0yLH9zJkztGjRguTkZPNbXIIk9EWpKiDwk1PSSUvPoKqLul2nhym7qrLhgj1WFgqLu0XTydfEwE+6Ckom+PQF90eK+kpEGVKspZWdnZ0JDw/Ptf3q1as4OcnKOkLkyYTA7zP6K7oN/5w7cUno9DB1dxU2XLBHq1GY3znG9MBPDAMUdZaOBL64i9k3cgcOHMjIkSP59NNPadOmDRqNhj179jBlyhSeffbZ4mijEOWfLk0dw88j8JOS0+j90lf8deAsjg42nAu7zbo4b9afd0CrUfiiUwzda6YW/DmKAomX1Fo9vk9C1abFdEGivDI79D/99FM0Gg1Dhw41LFNoZWXFyy+/zMcff1zkDRSi3MsunnZ7HzgZD/zHX5zP3wfP4eRgy+Yl49gQ9yA/nVMDf16nOzxW28TAT7igTvv0e1oevBJGFboMQ3JyMhcvXkRRFAICArC3N14UqqyRMX1RonTpWdMydxmdlpmUnMZjL8xn5yE18LcsHcdvSYGsPuuAhUbh84536FMn75XhcsiupeP/jLrYiqg0imVMX6fTceLECcPShPb29jzwwAM0bdoUjUbDiRMn0Ov199dyISoSXbpaWuH2bqOBn5iUSq9RX7Lz0DmcHW3ZtnQ8m5ObGwJ/bgdzAv+i+heE/0AJfJEvk0P/+++/5/nnn8fa2jrXe9bW1jz//POsWrWqSBsnRLmVmQxX10Hk31mlFXLPw4+MTuBC+O2swJ/A1tSHWPmvIxoUPmt/h74BJgZ+9hi+3wB1/Vwh8mFy6C9ZsoTJkyej1WpzvafVapk6dSqLFy8u0sYJUS5lJMCVH+F21k1bK+N/btf28+Cv7yexfdkEtqc9yHdn1MCf3T6WJ+qaGviXQWOtBr6M4QsTmBz6Z8+e5ZFH8p769fDDDxMaGlokjRKi3EqPhbBVcOcoONfNddM2ITGVPUcuGL6uW9OTrSkPsvSU+pfAx+1ieaqeic+6JF4CC2vwHwBVGhfVFYgKzuTQT0pKIj4+Ps/3ExISyvyDWUIUq/RYuLIG4k6pwyz3LHEYn5BCj5Ff0HXYXLbtPo2iwHsHXFh2Wg38mY/eYWB9cwLfRp2lU6VJEV+IqMhMDv26deuyb9++PN/fs2cPdevWLZJGCVHupMfBlbUQdxqcG6g98LvEJaTQ/fl57Dt2EXs7a9yqOPLufheWn1aHdD5pd4dnG5gY+AkXswJ/gAS+MJvJoT9o0CDeeustTpw4keu9f/75hxkzZjBo0KAibZwQ5UJGPIT/CLFZPXwLqxxvxyWk0H3EPA4cv0xVF3uCl09kU0JTwxj+J+1iTevhZ8/Dt7RTZ+nIkI4oBJPn6WdkZBAUFMSePXvo2rUrDRo0QKPREBoayh9//EHbtm0JDg7Gysqq4JOVIpmnL4pUepw6SyfmeNaQTs4efmx8Mt1HzOPQiTBcqzgQvGwCm+IfYHl24LePZYApY/iKAokXwNIp66Ztg4KPEZVGsRVcy8jIYO7cuaxatYrz58+jKAr16tVj0KBBTJgwweh0zrJGQl8UmfRYdUgnu4d/T+DHJ6TQdfjnHM4K/D+WT2RTfBOWnS5E4CecV2cByYNXwogSWRi9vJLQF0UiLUa9aRsfmjWkk7vDk5mpY+iUZWzfe4Y/lk/k94TGfHNSLUr4STsTb9oq+qzAr5r14JXcNxO5FetyiUJUeqm31R5+wtmsm7bGhzQtLbWsmD2CK9djWHuzjiHw/2fqLB1FD/FnwdZTnZbpWLsor0JUUmaXVhaiUku+AWErIeEcODfMFfgxsUl8MP93dDq1JIlWq+XHm3VYnBX4H7a9wyBTZukoOoj/F+xqQM1BEviiyEhPXwhTJV1Re/gpN9SnXzU5n06PiU2i67C5hJy5yp34JD6bNoBPDjuzKCvwP2gTy3MNTQh8fVbgO/irY/j2NYrjakQlJaEvhCniz6vTMtOi1SEdTc4/kqPvJNJ12OccD71KNTcnnn/qUT4+7MyiE/8F/pBGSQV/jl6n3idwrA3+z4KdZ3FcjajEJPSFKEj8WbWHn5Gg3rS9Z53ZqJhEug6byz//XsPT3Zk/vpvIT1EN+e6M+qTte60l8EXZYVLoT5o0yeQTzpkzp9CNEaLMiftX7eFnJoJTgNHA7zJsLieyA3/FJL6LaMhP5xzQoPDRo7GmjeFnD+k41oGaz4JttWK6IFHZmRT6ISEhOb4+evQoOp2O+vXVMq7nzp1Dq9USGBhY9C0UorTE/QvhayEzSQ38e+j1enqO+oIT/16juoczf3z3Gl9fbcivl+yx0KjlkU2qlpl909axtgS+KHYmhf5ff/1l+P9z5szBycmJ7777jqpVqwJw584dRowYQbt27YqnlUKUtLjQrB6+8cAHsLCw4J1XH+fV99ew+dtxfHWlIb9ftsfKQl3TtmctU5Y41N9103agBL4odmY/nOXt7c327dtp3Dhn3Y9Tp04RFBTEjRs3irSBRU0ezhIFunMCrq4HXSo41Slw9+SUDN44UI1NF9XAX9glhq7+pgb+2f+mZdp7F0HjRWVULMsl3n3yW7du5doeGRlJQkKCuacTomyJOab28PXpRgP/VlQ8PZ6fx8UrtwHQ6WH6QTXwLTUKX5kT+Annsx68GiiBL0qM2aH/xBNPMGLECNatW8e1a9e4du0a69atY+TIkfTv37842ihEyYg6BOE/qf/fsVaut2/ejqPTc3PYtvsMQ6YsRadXmLq7ChsvqIE/v0sMQeYEvrWr+qStg18RX4gQeTN7yubXX3/N5MmTee6558jIyFBPYmnJyJEjmT17dpE3UIgSEXMUrm1Ua+gY6XVHRMbRecgc/r10E5/qVVk+awRv7qnK+vMOaDUKX3SOoUdNUwJf+S/waz4rT9qKElfogmtJSUlcvHgRRVEICAjAwcGh4IPKABnTF7nEnlTn4YPRXndEZBydhnzG2Uu38PWqyp8rXmPptbqs+tcBC43CvI536F3HlFk6ilq+QQJfFLFiHdPPFhERQUREBPXq1cPBwYFKVqxTVBTxZyF8vTpt0kjg37gVS8fn/gv8v75/je9vBLDqX3Ue/mftTQ38rJu2EviilJkd+tHR0XTp0oV69erRq1cvIiIiABg1ahSvvfZakTdQiGKTeBmu/AS6JKNj+ADjP1zLucu38Kvhyl/fv8YPNwIMK17Nah9r3jx8O0+oOVgCX5Qqs0N/4sSJWFlZER4ejr29vWH7wIED2bp1a5E2Tohik3xdnaWTfgccjc/DB1j43mD6dGnGX9+/xuqIOnctYh7L06YsgKLPVOf82/tBzefAsWYRXYAQhWP2jdzt27ezbds2fHx8cmyvW7cuV65cKbKGCVFsUiOzqmVGZBVPy1laISU1HTtbdVEUd1dHfl74Ch8fcubbU//Vw3/GpNIKGeqQjlMA+A2UWjqiTDC7p5+UlJSjh58tKioKGxubImmUEMUm/Q5c+RGSwrKKp+X8EbgaEUPTxz9g0epdgHrv9d7yyKbV0slQyzg4N1CHdCTwRRlhdui3b9+eFStWGL7WaDTo9Xpmz55Np06dirRxQhSp7MCPP5sV+Dnr4YffiKHj4M+4cCWS2d9uJzklnf8dcubrrPLIplfLTFfH8Ks0Vm/a2rgVx9UIUShmD+/Mnj2bjh07cuTIEdLT05k6dSqnT58mJiaGvXv3Fkcbhbh/qZE5A/+eFa+uXI+m03NzuHwtijp+Hvz5/SRmH/cwjOGbXA9fl65+RpUm6pO21lWK4WKEKDyze/qNGjXixIkTtGzZkm7dupGUlET//v0JCQmhTp2C65QIUeKSr8PlH9Q58i4NQJtzEfOwa1F0fO6zuwL/Nb65XMsQ+P979I6JgZ+iBn7VZuqKVxL4ogwy++Gs8PBwfH190dxz8yv7PT+/sv1IuTycVckkXYEra7Ju2uYe0lEDfw5XrkcT4F+NHSsm8fWlWoZ5+J+0i2WAKYuYZyZB4iVwfRh8+4OVYzFdkBC5FevDWbVq1eL27du5tkdHR1OrlvG5zkKUiuRrWYF/K2uWjjbXLj8HH+fK9Wjq1qzGn99PYuHF/wJ/dnsTAz8jXp3z79EW/J6WwBdlmtlj+oqiGO3lJyYmYmtrWySNEuK+pURkBf5No7N0so0f3gULCwueCHqI+edrsfacg3kLoGTEQ1I4eHaEGo/nGjoSoqwxOfSzl0zUaDS8/fbbOaZt6nQ6Dh48yIMPPmh2AxYsWMDs2bOJiIigcePGfP755yYtxrJ37146dOhAkyZNOH78uNmfKyqw7Hn4SVfBpWGuwL9yPRr3qo442Nug0Wh4dWhn3thdhR+zAn9uhzv0DTAn8DuBd2+wkCWnRdln8r/S7CUTFUXh5MmTWFv/16OxtramWbNmTJ482awPX7t2LRMmTGDBggW0bduWRYsW0bNnT86cOZPvvYG4uDiGDh1Kly5djNb2F5VYaqTaw0+8ZDTwL1yJpNNzcwjw9+C3xa9ia2vD67ursO58VuB3vENfU2rpZCT818P3flwCX5QbZt/IHTFiBPPmzSuSm6CtWrWiefPmLFy40LCtYcOG9OvXj5kzZ+Z53DPPPEPdunXRarX8/PPPZvX05UZuBZZy87/Ad26QK4jPh92i03NzuH4rloZ1vPjj+0l8esafDeft0WoUPje1WmZGIiRdBo/24NNXhnREqSvWG7mff/45mZmZubbHxMQQHx9v8nnS09M5evQoQUFBObYHBQWxb9++PI9btmwZFy9e5J133jHpc9LS0oiPj8/xEhVQ8g0IW6XeUHVuaDTwO2YFfqOArMA//V/gf9EpxrTAz0zOumnbTgJflEtmh/4zzzzDmjVrcm3/8ccfeeaZZ0w+T1RUFDqdDk/PnI+ne3p6cvPmTaPHnD9/njfeeIOVK1diaWnan9MzZ87ExcXF8PL19TW5jaKcSL4OV1ar0zNdGoBFzlk6Zy/dpMPgz7hxK5bGdWvwx/eTmHXKnw0X1MD/snMMj9U2YQEUXSokXAT3R8CnjwS+KJfMDv2DBw8aLbfQsWNHDh48aHYD7p0JlNfsIJ1Ox6BBg3jvvfeoV6+eyeefNm0acXFxhtfVq1fNbqMow5KvQ9hqSM6+aZs78DsNmUNEZBxN6tUgeMVEPjnpz89Za9rO7xxDr1qmBH66uuKVayD4PgFamakmyiez7z6lpaUZHd7JyMggJcWEP4+zuLu7o9Vqc/XqIyMjc/X+ARISEjhy5AghISG8+uqrAOj1ehRFwdLSku3bt9O5c+dcx9nY2EghuIoqO/BTrmXNw8/dh0lOSSc1LYMH6nuzbdlEPj7pxy93BX4PkwP/LFRpCn5PgmXugoNClBdm9/QffvhhFi9enGv7119/TWBgoMnnsba2JjAwkODg4Bzbg4ODadOmTa79nZ2dOXnyJMePHze8Ro8eTf369Tl+/DitWrUy91JEeWYI/Ot5Bj7AQ439+Ov7SWxfPomZp+4K/C6mBn6qWjzNpUnWg1dORXwhQpQss3v6H330EV27duWff/6hS5cuAOzYsYPDhw+zfft2s841adIkhgwZQosWLWjdujWLFy8mPDyc0aNHA+rQzPXr11mxYgUWFhY0adIkx/HVqlXD1tY213ZRwSVfg7A1WYGf+8GrM+dvEJuQQpvmai2oRvV9mfBXVX6/rAb+V11i6G7KIua6FIi/AK7Nwe8psJLZXqL8Mzv027Zty/79+5k1axY//vgjdnZ2NG3alCVLllC3bl2zzjVw4ECio6N5//33iYiIoEmTJmzevBl/f39AXYc3PDzc3CaKiiz5WlYP/4bRwD99/gadnptDaloGf34/kWaNazLuT1e2hNlhZaGwsEsMXf1NCPzMJEi4BO4t1Vo6lg7FdEFClCyz5+mXdzJPvxxLvqZOyzQUT8sZ+KfOXafzkLncjkngoUa+bF46gRnH/Nh+xQ5rC4Wvu0bT2S+t4M/JDnyP1uDzBFjaFdMFCVE0inWePsDFixd56623GDRoEJGRkQBs3bqV06dPF+Z0QhSsgMA/efY6nZ6bw+2YBJo39mPLsom8nR34WoVF3UwM/IxEdR5+tbZZPXwJfFGxmB36O3fu5IEHHuDgwYOsX7+exMREAE6cOGHyA1NCmCUlImsM33jgn/j3Gp2HzCHqTiKBTfzYvHQCbx/1JTgr8L/pFk0nXxMDPykMPB5Ve/gyLVNUQGaH/htvvMGHH35IcHBwjvo7nTp1Yv/+/UXaOCEMxdOSrxoN/HOXb9F5qBr4LR7wV4d0jvoaevjfdIumg48JgZ8emxX47cCnH2hlmq+omMy+kXvy5ElWrVqVa7uHhwfR0dFF0ighAEiLVgM/8bL6pK2RaZn+NVxp1bQWt+8k8vu343nriA/bssbwF3c1MfBTb6uv6t2gRo9cSykKUZGYHfpVqlQhIiIi14IpISEheHt7F1nDRCWXFqMGfsI5tZaOkQVQAGxsrFj/1WjiknVMOejN39dssbZQx/A7mjKkk3wNMlPUsgqeHfOc7y9ERWH2v/BBgwbx+uuvc/PmTTQaDXq9nr179zJ58mSGDh1aHG0UlU3qbfWmbfy/RqtlHjsdzttzfyF74pnOwpoJ+3z4+5ottlo9S4JMHMNPvAyKXn3oqnpnCXxRKRTq4azhw4fj7e2Noig0atTIUBfnrbfeKo42isok5WbWkM5Fo4F/9NQVug3/nDtxyXi6OzN0YGeGb3XjWKQNjlZ6lnaPpmX19Pw/Q1Eg8QJoHdSyClUeKMYLEqJsKfQ8/YsXLxISEoJer+ehhx4y+8Gs0iLz9Muw5BtqPfykK1mBn3NI58jJMLoNn0dsfDKtH6rNjwsnMGaPD//ctsbFRs+KHlE088jI/zMUvVo4zaoq+D+t3hwWopwzJ9cKvdxPnTp1qF27NpC7UqYQZste8So5POumbc7AP3wijG7DPycuq7zC2oUTeGWXNyeirKlqo2NlrygaueUuBJiDoof4s2DrCf4DwLF2MV6QEGVToQYxlyxZQpMmTbC1tTXUvvn222+Lum2issi+aZsUljUtM2fgH/rnMl2HzSUuIYW2gWrgv5wV+K62OlaZFPg69R6BXQ2oOUgCX1RaZvf03377bebOncvYsWNp3bo1APv372fixImEhYXx4YcfFnkjRQWWHgfhP+Y5Syc2PpkeI78gPjGVdi0CWDV/PKN3enMq+r/Ab+BqQuDH/QsOfuD/DNjLLDNReZk9pu/u7s6XX37Js88+m2P76tWrGTt2LFFRUUXawKImY/plSEY8hP8EMcezVrwyPj/+h18OsHTdPlbMe5UXd3pzJtoat6zAr19Q4OuzevgO/lDzWbDzKvrrEKKUFeuYvk6no0WLFrm2BwYGGl1cRQij0mLUwI89Bc71cgX+3SuoPdf3Ebp3fYQh26rxb4wV7rY6Vj0WRb2qBQV+phr4jnXUHr5d7sV5hKhszB7Tf+6551i4cGGu7YsXL2bw4MFF0ihRwaVGqvPwY0+BS/1cJQ/2HbtIyydnEhEZB0BUigXPbVUD38NOxxqTAj9DDXyneuoYvgS+EEAhZ+8sWbKE7du388gjjwBw4MABrl69ytChQ5k0aZJhvzlz5hRNK0XFkXwDwrNLKzTMNQ9/79EL9Bj5BYlJacyYt4mPpg9j8GZ3zsdaUc1ex+peUdSpUkDgZy9v6NwI/AeCjWsxXpAQ5YvZoX/q1CmaN28OqHP1Qa274+HhwalTpwz7yTROkUvyDbiyGpLCjS5ivvvweXqO+pKk5DQ6t67P9EmDeOZ3dy7FWVHdXsfqx25Ty0WX/2foUiH+vPrAlf8AsK5SfNcjRDlkduj/9ddfxdEOUdElX1cDP/ma0Xn4uw6do9cL80lKTqNrm4YsnjOWETtqcDnekhoOmax+LAp/54ICP3t5wwez1rOVG/VC3MvsMf1bt27l+d6JEyfuqzGigspexDz5WtYi5jkDf+fB/wK/W9uGLJozluFZge/tmMnax00I/MxkNfDdAsFvoAS+EHkwO/QfeOABNm3alGv7p59+SqtWrYqkUaICMaxpmx34Of/J6fV6Jnz0I0nJaQQ92ogFn45j2B81uBJvia9TJmsfi8LXqaDAT1Jr9bi3Ar8BYOVYjBckRPlmdui//vrrDBw4kNGjR5OSksL169fp3Lkzs2fPZu3atcXRRlFeJYVD2EpIuW408AEsLCz4ddEYXnqmPfNnj2PYH15cTbDE31kNfB9TAj/hEri3Ab+nwNK+mC5GiIqhUAXX/vnnH5577jlSU1OJiYnhkUceYenSpXh6lv1pcfJwVglJDFNr6aTeMrriVVRMIu6u//XIw+K0DNrszo0kS2o5Z7L6sdtUd9Dn/xl3r2fr84SsdiUqrWJfGL127do0btyYsLAw4uPjGTBgQLkIfFFCEi+rN21TI40G/h97Q6nV+U3W/n4YgPN3LBnwuwc3kiyp45LB2sdNDPykMKj2qAS+EGYwO/T37t1L06ZNuXDhAidOnGDhwoWMHTuWAQMGcOfOneJooyhPkq6oPfy0KPVJ23sCP3jPGXq/9BWJSWms+f0Ip6MseeZ3dyKTtdSvmsGax6KoZm9K4F8Gj7aynq0QZjI79Dt37szAgQPZv38/DRs2ZNSoUYSEhHDt2jUeeEAWo6jUkq9lDencBqe6uQJ/+2418FPTMujduSlvTn+FQVs8iE7V0sQtnTWPReFRYODHZy1g3l4CX4hCMHue/vbt2+nQoUOObXXq1GHPnj189NFHRdYwUc4k34CwNerKV0aGdLbtPk3f0QtIS8+kb9dmvD59LCP+qEZihgXNq6WxrHs0LjYF3F5Kj1V/sXh2Au/HZQFzIQqh0CtnlVdyI7cYpNxUa+kkhxudpbN11yn6vbyQtPRM+nV7kIlvvMroP6uRnGlBq+ppLAmKxtG6gH+GaTGQEgHVu0KNnrnKNwhRmRXLjdxevXoRFxdn+Pqjjz4iNjbW8HV0dDSNGjUyv7WifEuJUAM/6YrRHj5A8N5Q0tIzeaLbg7w6dRwv7PAkOdOCdt6pLO9hQuCn3lZnAXl1hxq9JPCFuA8m9/S1Wi0RERFUq1YNAGdnZ44fP25YMvHWrVvUqFEDna6AedWlTHr6RejuNW2NlFbIpigKy9bto0ZgJ8b97U66XkNXvxTmd47BtqD8TrkJGXHg1QuqdzL6S0WIyq5Yevr3/m6oZKNC4l7ZxdPyWNN2f8hF0tPVapgajQafll0YmxX4PWumsKCLCYGffE19+Mq3P1TvLIEvRBGQnyJhPkPgXzW6pu2mHf/QYfBnDBi/mPT0TLaG2fLKH65k6DU8XjuZLzvHYG38j4L/JIapyxz6PaVOzZSqrUIUCZMHRzUaTa5yyVI+uRIqIPB/Dj7OgPGLycjQYWNtSfBVe8b/7UqmoqF37WTmdryDZX5dDUVR6+ho7cD3SajatHivR4hKxuTQVxSF4cOHY2OjzotOTU1l9OjRODg4AJCWllY8LRRlR0pE1hi+8SGdjdtDGDB+MZmZep557GGeGTOGcX+7oVM09K2TzGcdCgp8PSScB6uq4P+0+ktFCFGkTA79YcOG5fj6ueeey7XP0KFD779FomxKualWy8zjpu2GbccYOOEbMjP1PPv4wzw5egwTd7mhVzQ8EZDMp+3voM0v8PU6dbUr2+rqaleOtYr3eoSopEwO/WXLlhVnO0RZlnIr31k667eqga/T6RnUuyV9XxzD5N2uKGh4ul4SHz8aW0Dgp0P8OXCsrS5+Yu9dvNcjRCUmE55F/lJvZ61pe8noEocAHq6O2Fhb0j/oIXo8P4bX91RFQcOzDZL4qG0sFvnd+tGlQMIFcGmi3rS1cSu+axFCSOiLfKRFw5W1aig7Gw98gPYt63Fo/TQOpwUwfV9VAIY2SuS91nH5T7ox1MJvqdbRsXIq+msQQuQgUzaFcel3sgL/rFpawSL3kM6Jf68Zvt6fWo93D6iBP6pJQsGBf3ctfL+nJfCFKCHS0xe5pcepgR8Xqo7h31P2YM1vhxn82hJcXRw4svFNfrvtz+wjLgCMaZbA5BbxBQd+0mXweFQqZQpRwiT0RU4Z8RD+I8Sdzgr8nJUsV206xJApS9HrFXp3aca6CH++PK4G/sTm8Yx7KMGEwA8Dj3bg01cCX4gSJqEv/pORCOE/wZ0T6hz5ewJ/5S8HGTp1GXq9wvNPt6VO3zF8eVyt8/H6w3G83CyxgPPHq3P8DYFvXVxXIoTIg4S+UGWmwNX1EHMcnOvmCuQffjnAsKnL0esVRg54FO/HxrD4pDoO/84jsYxokpT/+Q218DtCjccl8IUoJRL6AnRpcG0jRB8B5wDQ2uZ4e/PfJxk6ZTmKovDCwHa49xjDd2cc0aDw0aOxDGqQnP/5DbXwu0ktfCFKmfz0VXb6DLi2CW7vB6faas2be7RrUZc2zWvTpJ43Vbu9wqp/1cCf3T6Wp+oVFPhR6lz/Gj2gelCuWUBCiJIloV+Z6TPh+u9we7da9sDSwehuTo62bF06gY+OerD6rBr4czrc4Ym6KfmfP/U2pEervfvqXaU0shBlgPwUVlb6TLixGW79Cfa+YOWY4+2lP+3lowWb1V0V+OhoNTMD/xakx4DXYxL4QpQhpf6TuGDBAmrVqoWtrS2BgYHs3r07z303bNhAt27d8PDwwNnZmdatW7Nt27YSbG0FodfBja1wcwfY+4C1S463v/1xDyPfXMFbc39h864zTN5ZldVnHUwP/JSb6lx/796y+IkQZUyp/jSuXbuWCRMmMH36dEJCQmjXrh09e/YkPDzc6P67du2iW7dubN68maNHj9KpUyd69+5NSEhICbe8HFP0ELENbv0B9jXAukqOt79Zu5sXpn8PwKtDu/Jrels2XLBHq1H4vJMJgZ90VS2v4NMPqnWQxU+EKGNMXiO3OLRq1YrmzZuzcOFCw7aGDRvSr18/Zs6cadI5GjduzMCBA5kxY4ZJ+1fqNXIVPURsV3v5dtXBxjXH24vX7OKlt1cC8OqwbiQ9/CJ/XrXD2kJhfpcYgvxT8zm3ohZls7AG337gGliMFyKEuFuxrJFb1NLT0zl69ChBQUE5tgcFBbFv3z6TzqHX60lISMDV1TXPfdLS0oiPj8/xqpQUPdz8Q+3l23nmCvyvV+00BP7YEd250/wl/rxqh41W4Zug6AICXw8J59Qbwf7PSOALUYaVWuhHRUWh0+nw9PTMsd3T05ObN2+adI7PPvuMpKQkBgwYkOc+M2fOxMXFxfDy9fW9r3aXS4oCN/9Se/i2HrnKF586d51X3l0NwLjnexDV7AX23LDF3lLPdz2i6OCTz6po+kyI/xdsPKDmYKjSuDivRAhxn0p9yua96+wqimLS2rurV6/m3Xff5ZdffqFatWp57jdt2jQmTZpk+Do+Pr7yBX/kTnWmjo0b2LjnertJPW/mvTWASzfiiXhgJAcibHGw0rO8ezQPV0/P+7y6dIg/qz7B6/uUeo9ACFGmlVrou7u7o9Vqc/XqIyMjc/X+77V27VpGjhzJTz/9RNeuXfPd18bGxrCub6UUdVCdi29TRe3l3yU9PRNra/WfwPBnujByuxuHImxwtNLzXY9oAj3zCfzsWvhVm6oLmNvkPcQmhCg7Sm14x9ramsDAQIKDg3NsDw4Opk2bNnket3r1aoYPH86qVat47LHHiruZ5dudE3DtF7C0B9ucv0i/+O5P2gz8hDtxSdxJtWDwZncO3bTByUrP9z2j8g/8jET1pq17a6g5SAJfiHKkVId3Jk2axJAhQ2jRogWtW7dm8eLFhIeHM3r0aEAdmrl+/TorVqwA1MAfOnQo8+bN45FHHjH8lWBnZ4eLi0uen1MpxZ+HqxsAJdeas58v+4OJ//sJgMW/nuIP296cj7XC1VbHih7RNHHPyPu86XGQfBU82oNPHymNLEQ5U6qhP3DgQKKjo3n//feJiIigSZMmbN68GX9/fwAiIiJyzNlftGgRmZmZjBkzhjFjxhi2Dxs2jOXLl5d088uupHC4+hNkJoJzvRxvzV32B5OyAn/8KwP43aoP4bGWeNrrWNkzioCqmXmfN/0OJN8Az07g/Xiu0stCiLKvVOfpl4YKP08/+TqErVSrWjrXz/E07GdLgpn88ToAxo15lkOeA4lIssTPKZOVvaLwddLlfd7swmmeXdTiaVIpU4gyw5xck5/ciiTllrrMYcoNdV3buwL/02+3M+WT9QCMHTuEvW5PEZWkJaBKBit7RuHpoM/7vMnX1Ru3NXqBZ2eplClEOSahX1GkRkH4Wki6oi5zeFfgx8Yn8/nyHQC8Om44u6r0JzbFgkZu6XzfIxo3uzwC3/CUraW6eLlbSymrIEQ5J6FfEaTHqevaJlzI6uHn7IlXcbbnrx8msXhHBFutepKQZsGDHul81yMKF5s8RvcUvXoz2NpFnZIpD10JUSFI+cPyLjMJrq6DuNCsdW3/+z1+Kfy24f/ftvbld21PEtItaFk9jR96FhT4Z9V5/bWGSOALUYFI6JdnulS4ujFrXdt6OWbTfLRgM416vsu23af5M9yG4dvcSc60oJ13Kt/1iMbRuoDAt/NS5+A71iqZaxFClAgZ3imv9Blw7Vf1iVungBzz5T+Y/zsz5m0CYEOoJTvOuZGpaOjmn8L8zjHY5HUfVq+DhLNqjX3/Z9T/FUJUKBL65ZE+E65vvmuZw//WtX3vy19594vfABg+ZQrB2vboFQ196iTzWYc7WOX1t50uVR3Dd6oDfgOkjo4QFZSEfnmTverVPcscKorCu1/8yvvzfwdgyJQ3+cuiDSgwsF4S/3s0Fm1egZ8Rr876cX1IXfxEyioIUWFJ6Jcnih5uBsOtHVmrXqmlJxRF4Z15v/LBV2rgPzt1Brs0LQEY3jiRGY/EYZHXTMvUSHW6p2cndR6+1rYkrkQIUUok9MuL7Jr4EdvV4mnWVe96S+Fi1kydp17/iH00A+CVZglMaRGf99T65KugSwOfvlCtvTx0JUQlIKFfXkQfhIgtWTXxcw6/WFhYsOyT4Vi2HsbOOLVu0ZQWcYx5MNH4uQwPXVmp4/duLYq79UKIMkKmbJYHd07AtU3qcoRZNfEVRWHdlqPodHoy9TBtn5sh8D9oE5tP4N+9tOGzEvhCVDLS0y/rEi7AtY3cXSJZURTe/OxnPl60laH922LfbTxbwuzRahQ+bX+HJ+qmGD+XPlOdg29fQy2r4Fi75K5DCFEmSOiXZcnXIHwdZCQYSiQrisIbszcw65vtYGHJ9XqDuRBmj7WFwvwuMXkvYK5LhYTz4FQP/J5SH74SQlQ6EvplVUqEWjEz9ZZaXgE18Kd+sp5PlwSD1pLOU+dzQe+DtVZhUddoOvnmsYB5RiIkXgbXB9U6OtZVSuwyhBBli4R+WZRdIjkpHFwagsYCRVGY8sl6PlsSjMbShg5TvuSivgY2WoVvu0XTziePwE+LUUstV3sUvHurSycKISotCf2yJrtEcuIlQ+ADvDF7A58tCcbC1olWE7/kst4dO0s9S4KiaVMjj/Vsk6+rq2fV6AnVu8pKV0IImb1TpuRTIrldi7rYuVbnwQlfc0PvjouNnu975BH4il49h6JTb9h69ZDAF0IA0tMvOzKT1YXM40LVRVDuWY6w/kPNaTK2C5Gp1ng5ZPJdj2jqGVvPVp8B8efAzhN8nlD/WhBCiCzS0y8LdOlwfRPEHAPnumBhZailc/HKbU5HWzHwd3ciU60JqJLB+t5RxgM/M1n9peFUF2oNl8AXQuQiPf3SpughYivc3gdOtUFri6IojH1/DV/98DfLd97Eue/bxKdraeqeznc9oqlqa2R5w+wbtu6PgE8fsKqAi74LIe6bhH5pyq6nc/NPtXa9pQN6vZ5X31vDwlU7sfFphHXP6cSna2leLY3lPaJxNrb4ScoNSE8Ar+5QvRtorUv+WoQQ5YKEfmmKPqTW07H1AOsq6PV6xry7mq9X78KuVnN8BswgHUtaVU9jSZCR1a4UBZIuAxbg9yS4t5aFy4UQ+ZLQLy2xp9Rx/Kx6Onq9npdnrGLx2t04NO5Etccnko66vOHibjHYWd4b+Dr1hq21G/j2gypNSuUyhBDli4R+aUi8pM7UUfTgoNbT+XjRVhav3Y1Lyyeo0mkkeqBPnWQ+bX8H63srHutSIP6CumqW35Pg4F/ilyCEKJ8k9Etayk0I/wky4tQ6OFleHNieFVf8SQ0IAuD5Jom81crI4idpMZByHVwD1Ru2ssqVEMIMEvolKT1WDfzkG+DSEAXQAHoF5p7xJjVArbHz+sNxjG6amHt4PilcLZzm1UNu2AohCkVCv6RkpqhDOgnnwLkhOj28MH0FgU1qcd67Lz9ftEeDwv8ejeXZBsk5j1X0kHBRrZtT81mo2lxu2AohCkVCvyToM+H6rxATAs510SkWjJy2gu82HeF3pRN2qfZYahQ+63iHvnXuqYWv10HCWXWJRN8nDSWWhRCiMCT0i5uiqOva3t4LjrXQYc2IN77jh80nqD7wfWx8H8Baq7Cgcwxd762Fr0tXA9+xjloD396ndK5BCFFhSOgXt6h9cHMH2Hmhs7Bn+OvLWf3HOaoP/hjrarVxtNKzuJuRwmnZNfCrNFED38atdNovhKhQJPSLU+xJuP47WDuTqXVh2JRl/LjnGtWfm41lleq42+lY3j2aJu4ZOY9Li1Jr6ld7FLwfV+fyCyFEEZDQLy6JYXA1a21bOy82bTvGT/tvUX3wJ2gdquLvnMn3PaLwc9blPC4pHHRp4P0YeHbOVW1TCCHuhyRKcUiLUWfqpMcaljoMaN6KWs93I01jS0PXdFb0iMbD/q7Cadk18C0doOYz6jx8IYQoYhL6RU2XBtc2QeJlMh3rk56awek4R0ZscyNNY8GDHul81yMKF5u7yiroMyD+rHqj1rc/OAWUXvuFEBWahH5RUhSICIY7R8mwq83g15YTpvEj4cERpOgsaFk9jaX3Fk7LTFLLMjg3VKdk2nmWXvuFEBWehH5RijkKkX+RYeXFs1N+YHuUL1U7j0Cjs6CDTypfd72ncFpKhDoU5NYafHqDlVPptV0IUSlI6BeVpHC48RsZehueeXMjf+tb4drlMQAGNUjivTaxWGWvU6bP/G/83u9pdeETuWErhCgBkjRFISMerv1MelIMA97Zzz6nXjgFtESDwpst4xn1wF11dLKHc5zqgndvtVKmEEKUEAn9+6XXwY0tpMf8y1PvHuWQ+1PY+zbBSqNjfpdYute86ylbw5KGrdXAl+EcIUQJk9C/X1H7IGo/xyPdCfEZgq1HLewsMljRK5aHq9/1lG3yNbWXX6MnVO8KFlal12YhRKUloX8/4s9DxDbC06rzxpkGaD2scNKmsaZPHI3dsp6y1WdkVch0BL+B4PawVMgUQpQaCf37kHZ1Kz8fSWTW1ZZEp2rxdszkh553qOWS9ZRt+h1IuqZWxvR+DBxrl26DhRCVnoR+IaWmptL17X1crT0AjaWWJm7pLAmKxtNBrz5dm3hZnbfvFQTVO0v9HCFEmSChXwgpKSm0f/l/3K73HBrgIafbrHw8A3srRV2/NuGC+nRtjcfApbEM5wghygyLgncpXgsWLKBWrVrY2toSGBjI7t27891/586dBAYGYmtrS+3atfn6669LqKWq1NRU+j01kOuWNQDoXuUs655OVwM/9bY6fl81EGqPVMsiS+ALIcqQUg39tWvXMmHCBKZPn05ISAjt2rWjZ8+ehIeHG93/8uXL9OrVi3bt2hESEsKbb77JuHHjWL9+fYm0NyUlhb59+7J9868k/P4JI6qfYFGvOLT6FIj7FzIT1amYNQeBrXuJtEkIIcyhURRFKXi34tGqVSuaN2/OwoULDdsaNmxIv379mDlzZq79X3/9dTZt2kRoaKhh2+jRo/nnn3/Yv3+/SZ8ZHx+Pi4sLcXFxODs7m9zW7MAPDg7G3t6ezZs306HaUXVlK0sHtXZOtY7gXNfkcwohRFEwJ9dKraefnp7O0aNHCQoKyrE9KCiIffv2GT1m//79ufbv3r07R44cISMjw+gxaWlpxMfH53gVhoWFBVZWVjg4OLBlyxY6dOighr1LY6g1DGqPkMAXQpR5pXYjNyoqCp1Oh6dnzqqSnp6e3Lx50+gxN2/eNLp/ZmYmUVFReHl55Tpm5syZvPfee/fdXhsbG9avX8+///7Lgw8+qG70HwhaO9Da3Pf5hRCiJJT6jVzNPTc6FUXJta2g/Y1tzzZt2jTi4uIMr6tXrxa6rba2tv8FPoB1FQl8IUS5Umo9fXd3d7Raba5efWRkZK7efLbq1asb3d/S0hI3N+MLh9vY2GBjI8EshBBQij19a2trAgMDCQ4OzrE9ODiYNm3aGD2mdevWufbfvn07LVq0wMpKatkIIURBSnV4Z9KkSXz77bcsXbqU0NBQJk6cSHh4OKNHjwbUoZmhQ4ca9h89ejRXrlxh0qRJhIaGsnTpUpYsWcLkyZNL6xKEEKJcKdUncgcOHEh0dDTvv/8+ERERNGnShM2bN+Pv7w9AREREjjn7tWrVYvPmzUycOJGvvvqKGjVq8MUXX/Dkk0+W1iUIIUS5Uqrz9EtDYefpCyFEWVUu5ukLIYQoeRL6QghRiUjoCyFEJSKhL4QQlYiEvhBCVCIS+kIIUYlUupWzsmeoFrbaphBClDXZeWbKDPxKF/oJCQkA+Pr6lnJLhBCiaCUkJODi4pLvPpXu4Sy9Xs+NGzdwcnLKt5rnveLj4/H19eXq1asV9qGuin6NFf36oOJfo1yfcYqikJCQQI0aNbCwyH/UvtL19C0sLPDx8Sn08c7OzhXyH9vdKvo1VvTrg4p/jXJ9uRXUw88mN3KFEKISkdAXQohKRELfRDY2NrzzzjsVekGWin6NFf36oOJfo1zf/at0N3KFEKIyk56+EEJUIhL6QghRiUjoCyFEJSKhL4QQlYiE/l0WLFhArVq1sLW1JTAwkN27d+e7/86dOwkMDMTW1pbatWvz9ddfl1BLC8+ca9ywYQPdunXDw8MDZ2dnWrduzbZt20qwteYz979htr1792JpacmDDz5YvA28T+ZeX1paGtOnT8ff3x8bGxvq1KnD0qVLS6i1hWPuNa5cuZJmzZphb2+Pl5cXI0aMIDo6uoRaa55du3bRu3dvatSogUaj4eeffy7wmCLPGUUoiqIoa9asUaysrJRvvvlGOXPmjDJ+/HjFwcFBuXLlitH9L126pNjb2yvjx49Xzpw5o3zzzTeKlZWVsm7duhJuuenMvcbx48crn3zyiXLo0CHl3LlzyrRp0xQrKyvl2LFjJdxy05h7fdliY2OV2rVrK0FBQUqzZs1KprGFUJjr69Onj9KqVSslODhYuXz5snLw4EFl7969Jdhq85h7jbt371YsLCyUefPmKZcuXVJ2796tNG7cWOnXr18Jt9w0mzdvVqZPn66sX79eAZSNGzfmu39x5IyEfpaWLVsqo0ePzrGtQYMGyhtvvGF0/6lTpyoNGjTIse2ll15SHnnkkWJr4/0y9xqNadSokfLee+8VddOKRGGvb+DAgcpbb72lvPPOO2U69M29vi1btiguLi5KdHR0STSvSJh7jbNnz1Zq166dY9sXX3yh+Pj4FFsbi4opoV8cOSPDO0B6ejpHjx4lKCgox/agoCD27dtn9Jj9+/fn2r979+4cOXKEjIyMYmtrYRXmGu+l1+tJSEjA1dW1OJp4Xwp7fcuWLePixYu88847xd3E+1KY69u0aRMtWrRg1qxZeHt7U69ePSZPnkxKSkpJNNlshbnGNm3acO3aNTZv3oyiKNy6dYt169bx2GOPlUSTi11x5EylK7hmTFRUFDqdDk9PzxzbPT09uXnzptFjbt68aXT/zMxMoqKi8PLyKrb2FkZhrvFen332GUlJSQwYMKA4mnhfCnN958+f54033mD37t1YWpbtH4XCXN+lS5fYs2cPtra2bNy4kaioKF555RViYmLK5Lh+Ya6xTZs2rFy5koEDB5KamkpmZiZ9+vThyy+/LIkmF7viyBnp6d/l3lLLiqLkW37Z2P7Gtpcl5l5jttWrV/Puu++ydu1aqlWrVlzNu2+mXp9Op2PQoEG899571KtXr6Sad9/M+e+n1+vRaDSsXLmSli1b0qtXL+bMmcPy5cvLbG8fzLvGM2fOMG7cOGbMmMHRo0fZunUrly9fZvTo0SXR1BJR1DlTtrs3JcTd3R2tVpurNxEZGZnrt2y26tWrG93f0tISNze3YmtrYRXmGrOtXbuWkSNH8tNPP9G1a9fibGahmXt9CQkJHDlyhJCQEF599VVADUlFUbC0tGT79u107ty5RNpuisL89/Py8sLb2ztHyd2GDRuiKArXrl2jbt26xdpmcxXmGmfOnEnbtm2ZMmUKAE2bNsXBwYF27drx4Ycflrm/uM1VHDkjPX3A2tqawMBAgoODc2wPDg6mTZs2Ro9p3bp1rv23b99OixYtsLKyKra2FlZhrhHUHv7w4cNZtWpVmR4nNff6nJ2dOXnyJMePHze8Ro8eTf369Tl+/DitWrUqqaabpDD//dq2bcuNGzdITEw0bDt37tx9rylRXApzjcnJybkWDdFqtYBpSweWdcWSM4W+BVzBZE8VW7JkiXLmzBllwoQJioODgxIWFqYoiqK88cYbypAhQwz7Z0+lmjhxonLmzBllyZIl5WbKpqnXuGrVKsXS0lL56quvlIiICMMrNja2tC4hX+Ze373K+uwdc68vISFB8fHxUZ566inl9OnTys6dO5W6desqo0aNKq1LKJC517hs2TLF0tJSWbBggXLx4kVlz549SosWLZSWLVuW1iXkKyEhQQkJCVFCQkIUQJkzZ44SEhJimJJaEjkjoX+Xr776SvH391esra2V5s2bKzt37jS8N2zYMKVDhw459v/777+Vhx56SLG2tlZq1qypLFy4sIRbbD5zrrFDhw4KkOs1bNiwkm+4icz9b3i3sh76imL+9YWGhipdu3ZV7OzsFB8fH2XSpElKcnJyCbfaPOZe4xdffKE0atRIsbOzU7y8vJTBgwcr165dK+FWm+avv/7K92eqJHJGSisLIUQlImP6QghRiUjoCyFEJSKhL4QQlYiEvhBCVCIS+kIIUYlI6AshRCUioS+EEJWIhL6otO5euSgsLAyNRsPx48cB+Pvvv9FoNMTGxuZ5/PLly6lSpYrh63fffTfHylvDhw+nX79+Rd7usqBmzZp8/vnnpd0MUQgS+qLMuHnzJmPHjqV27drY2Njg6+tL79692bFjR7F/tq+vLxERETRp0sTkYwYOHMi5c+fyfH/evHksX77c8HXHjh2ZMGHCfbRSde8vFyHMIVU2RZkQFhZG27ZtqVKlCrNmzaJp06ZkZGSwbds2xowZw7///mv0uIyMjCIpcKfVaqlevbpZx9jZ2WFnZ5fn+3dXtxSirJCevigTXnnlFTQaDYcOHeKpp56iXr16NG7cmEmTJnHgwAHDfhqNhq+//pq+ffvi4ODAhx9+CMCvv/6aY/Ho9957j8zMTMNx58+fp3379tja2tKoUaNclQvvHd7JtnfvXpo1a4atrS2tWrXi5MmThvfuHd65193DO8OHD2fnzp3MmzcPjUaDRqPh8uXLBAQE8Omnn+Y47tSpU1hYWHDx4kVzvoUG169fZ+DAgVStWhU3Nzf69u1LWFgYANu2bcPW1jbXsNW4cePo0KGD4et9+/bRvn177Ozs8PX1Zdy4cSQlJRWqPaJskdAXpS4mJoatW7cyZswYHBwccr1/b7C+88479O3bl5MnT/L888+zbds2nnvuOcaNG8eZM2dYtGgRy5cv56OPPgLUOvn9+/dHq9Vy4MABvv76a15//XWT2jZlyhQ+/fRTDh8+TLVq1ejTp0+hlqmbN28erVu35oUXXiAiIoKIiAj8/Px4/vnnWbZsWY59ly5dSrt27ahTp47Zn5OcnEynTp1wdHRk165d7NmzB0dHR3r06EF6ejpdu3alSpUqrF+/3nCMTqfjxx9/ZPDgwQCcPHmS7t27079/f06cOMHatWvZs2ePYd0BUc7dV7k2IYrAwYMHFUDZsGFDgfsCyoQJE3Jsa9eunfK///0vx7bvv/9e8fLyUhRFUbZt26ZotVrl6tWrhve3bNmSY2Hqy5cvK4ASEhKiKMp/1RDXrFljOCY6Olqxs7NT1q5dqyiKWtbXxcXF8P69VTqHDRum9O3b1/B1hw4dlPHjx+do540bNxStVqscPHhQURRFSU9PVzw8PJTly5fn+T3IrxrokiVLlPr16yt6vd6wLS0tTbGzs1O2bdumKIqijBs3TuncubPh/W3btinW1tZKTEyMoiiKMmTIEOXFF1/Mcd7du3crFhYWSkpKiqIoiuLv76/MnTs3zzaKskvG9EWpU8xc/q1FixY5vj569CiHDx829OxB7b2mpqaSnJxMaGgofn5+ORYOad26tUmfdfd+rq6u1K9fn9DQUJOONYWXlxePPfYYS5cupWXLlvz222+kpqby9NNPF+p8R48e5cKFCzg5OeXYnpqaahguGjx4MK1bt+bGjRvUqFGDlStX0qtXL6pWrZrjHCtXrjQcrygKer2ey5cv07Bhw0JerSgLJPRFqatbty4ajYbQ0FCTpjjeOwSk1+t577336N+/f659bW1tja6gdD/rGBf1GsijRo1iyJAhzJ07l2XLljFw4EDs7e0LdS69Xk9gYGCOwM7m4eEBQMuWLalTpw5r1qzh5ZdfZuPGjTmGmPR6PS+99BLjxo3LdQ4/P79CtUuUHRL6otS5urrSvXt3vvrqK8aNG5cr1GNjY/O9Ydq8eXPOnj1LQECA0fcbNWpEeHi4oWcLsH//fpPaduDAAUPQ3blzh3PnztGgQQOTjr2XtbU1Op0u1/ZevXrh4ODAwoUL2bJlC7t27SrU+UH9XmQvXu/s7JznfoMGDWLlypX4+PhgYWGRYynM5s2bc/r06Ty/n6J8kxu5okxYsGABOp2Oli1bsn79es6fP09oaChffPFFgUMxM2bMYMWKFbz77rucPn2a0NBQ1q5dy1tvvQVA165dqV+/PkOHDuWff/5h9+7dTJ8+3aR2vf/+++zYsYNTp04xfPhw3N3dC/3AVc2aNTl48CBhYWFERUWh1+sBdbro8OHDmTZtGgEBASYNPaWkpORY3/f48eNcuHCBwYMH4+7uTt++fdm9ezeXL19m586djB8/nmvXrhmOHzx4MMeOHeOjjz7iqaeewtbW1vDe66+/zv79+xkzZgzHjx/n/PnzbNq0ibFjxxbqukXZIqEvyoRatWpx7NgxOnXqxGuvvUaTJk3o1q0bO3bsYOHChfke2717d3777TeCg4N5+OGHeeSRR5gzZw7+/v4AWFhYsHHjRtLS0mjZsiWjRo3KMf6fn48//pjx48cTGBhIREQEmzZtwtraulDXOHnyZLRaLY0aNcLDw4Pw8HDDeyNHjiQ9PZ3nn3/epHOdO3eOhx56KMdr1KhR2Nvbs2vXLvz8/Ojfvz8NGzbk+eefJyUlJUfPv27dujz88MOcOHHCMGsnW9OmTdm5cyfnz5+nXbt2PPTQQ7z99tt4eXkV6rpF2SLLJQpRBuzdu5eOHTty7do1PD09S7s5ogKT0BeiFKWlpXH16lVefPFFvLy8jN6AFaIoyfCOEKVo9erV1K9fn7i4OGbNmlXazRGVgPT0hRCiEpGevhBCVCIS+kIIUYlI6AshRCUioS+EEJWIhL4QQlQiEvpCCFGJSOgLIUQlIqEvhBCViIS+EEJUIv8HMkSw6UYra9AAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "k_sigma = 3\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", + "ax.plot([0, 1], [0, 1], ls='--', color='k', label = \"Ideal case\")\n", + "ax.plot(alpha, ecp, label='DRP')\n", + "ax.fill_between(alpha_mean, ecp_mean - k_sigma * ecp_std, ecp_mean + k_sigma * ecp_std, color = \"orange\", alpha = 0.4, label = f\"Uncertainty zone ({k_sigma}\" + r\"$\\sigma$)\")\n", + "ax.legend()\n", + "ax.set_ylabel(\"Expected Coverage\")\n", + "ax.set_xlabel(\"Credibility Level\")\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "base", + "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.9.17" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/src/tarp/drp.py b/src/tarp/drp.py index 71ceabf..db45662 100644 --- a/src/tarp/drp.py +++ b/src/tarp/drp.py @@ -1,16 +1,20 @@ from typing import Tuple, Union import numpy as np +from tqdm import tqdm +from fast_histogram import histogram1d __all__ = ("get_drp_coverage",) -def get_drp_coverage( + +def old_get_drp_coverage( samples: np.ndarray, theta: np.ndarray, references: Union[str, np.ndarray] = "random", metric: str = "euclidean", -) -> Tuple[np.ndarray, np.ndarray]: + seed: Union[int, None] = None, + ) -> Tuple[np.ndarray, np.ndarray]: """ Estimates coverage with the distance to random point method. @@ -30,6 +34,8 @@ def get_drp_coverage( Credibility values (``alpha``) and expected coverage probability (``ecp``). """ + np.random.seed(seed) + # Check that shapes are correct if samples.ndim != 3: raise ValueError("samples must be a 3D array") @@ -92,3 +98,136 @@ def get_drp_coverage( dx = alpha[1] - alpha[0] ecp = np.cumsum(h) * dx return ecp, alpha[1:] + + +def get_drp_coverage( + samples: np.ndarray, + theta: np.ndarray, + references: Union[str, np.ndarray] = "random", + metric: str = "euclidean", + norm : bool = True, + seed: Union[int, None] = None + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Estimates coverage with the distance to random point method. + + Reference: `Lemos, Coogan et al 2023 `_ + + Args: + samples: the samples to compute the coverage of, with shape ``(n_samples, n_sims, n_dims)``. + theta: the true parameter values for each samples, with shape ``(n_sims, n_dims)``. + references: the reference points to use for the DRP regions, with shape + ``(n_references, n_sims)``, or the string ``"random"``. If the later, then + the reference points are chosen randomly from the unit hypercube over + the parameter space. + metric: the metric to use when computing the distance. Can be ``"euclidean"`` or + ``"manhattan"``. + norm : whether to apply or not the normalization (Default = True) + + Returns: + Expected coverage probability (``ecp``) and credibility values (``alpha``) + """ + np.random.seed(seed) + + # Check that shapes are correct + if samples.ndim != 3: + raise ValueError("samples must be a 3D array") + + if theta.ndim != 2: + raise ValueError("theta must be a 2D array") + + num_samples = samples.shape[0] + num_sims = samples.shape[1] + num_dims = samples.shape[2] + + if theta.shape[0] != num_sims: + raise ValueError("theta must have the same number of rows as samples") + + if theta.shape[1] != num_dims: + raise ValueError("theta must have the same number of columns as samples") + + # Reshape theta + theta = theta[np.newaxis, :, :] + + # Normalize + if norm : + low = np.min(theta, axis=1, keepdims=True) + high = np.max(theta, axis=1, keepdims=True) + samples = (samples - low) / (high - low + 1e-10) + theta = (theta - low) / (high - low + 1e-10) + + # Generate reference points + if isinstance(references, str) and references == "random": + references = np.random.uniform(low=0, high=1, size=(num_sims, num_dims)) + else: + assert isinstance(references, np.ndarray) # to quiet pyright + if references.ndim != 2: + raise ValueError("references must be a 2D array") + + if references.shape[0] != num_sims: + raise ValueError("references must have the same number of rows as samples") + + if references.shape[1] != num_dims: + raise ValueError( + "references must have the same number of columns as samples" + ) + + # Compute distances + if metric == "euclidean": + samples_distances = np.sqrt( + np.sum((references[np.newaxis] - samples) ** 2, axis=-1) + ) + theta_distances = np.sqrt(np.sum((references - theta) ** 2, axis=-1)) + elif metric == "manhattan": + samples_distances = np.sum(np.abs(references[np.newaxis] - samples), axis=-1) + theta_distances = np.sum(np.abs(references - theta), axis=-1) + else: + raise ValueError("metric must be either 'euclidean' or 'manhattan'") + + # Compute coverage + f = np.sum((samples_distances < theta_distances), axis=0) / num_samples + + # Compute expected coverage + h, alpha = np.histogram(f, density=True, bins=num_sims // 10) + dx = alpha[1] - alpha[0] + ecp = np.cumsum(h) * dx + return ecp, alpha[1:] + +def bootstrapping(samples, theta, references = "random", metric = "euclidean", norm = True): + """Estimates uncertainties on the expected probability and credibility values calculated with the get_drp_coverage function + using the bootstrapping method + + Args: + samples: the samples to compute the coverage of, with shape ``(n_samples, n_sims, n_dims)``. + theta: the true parameter values for each samples, with shape ``(n_sims, n_dims)``. + references: the reference points to use for the DRP regions, with shape + ``(n_references, n_sims)``, or the string ``"random"``. If the later, then + the reference points are chosen randomly from the unit hypercube over + the parameter space. + metric: the metric to use when computing the distance. Can be ``"euclidean"`` or + ``"manhattan"``. + norm : whether to apply or not the normalization (Default = True) + + Returns: + Mean and std of the expected coverage probability, and mean of the credibility + """ + num_samples = samples.shape[0] + num_sims = samples.shape[1] + num_dims = samples.shape[2] + + boot_ecp = np.empty(shape = (num_sims, num_sims//10)) + boot_alpha = np.empty(shape = (num_sims, num_sims//10)) + for i in tqdm(range(num_sims)): + idx_remove = np.random.randint(num_sims) + idx_add = np.random.randint(num_sims) + + # Replacing one simulation and its samples by another and its associated samples + samples[:, idx_remove, :] = samples[:, idx_add, :] + theta[idx_remove, :] = theta[idx_add, :] + + boot_ecp[i, :], boot_alpha[i, :] = get_drp_coverage(samples, theta, references = references, metric = metric, norm = norm) + + ecp_mean = boot_ecp.mean(axis = 0) + ecp_std = boot_ecp.std(axis = 0) + alpha_mean = boot_alpha.mean(axis = 0) + return ecp_mean, ecp_std, alpha_mean From 2806b4b545a5d5d0820d3cfc0da42fe15182d90f Mon Sep 17 00:00:00 2001 From: EnceladeCandy Date: Fri, 20 Oct 2023 11:14:06 -0400 Subject: [PATCH 2/9] Added description for bootstrapping function --- src/tarp/drp.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tarp/drp.py b/src/tarp/drp.py index db45662..a943a4f 100644 --- a/src/tarp/drp.py +++ b/src/tarp/drp.py @@ -194,7 +194,8 @@ def get_drp_coverage( return ecp, alpha[1:] def bootstrapping(samples, theta, references = "random", metric = "euclidean", norm = True): - """Estimates uncertainties on the expected probability and credibility values calculated with the get_drp_coverage function + """ + Estimates uncertainties on the expected probability and credibility values calculated with the get_drp_coverage function using the bootstrapping method Args: From e089486feb416fe1882ef945c096f7ad01a7d19b Mon Sep 17 00:00:00 2001 From: EnceladeCandy Date: Fri, 20 Oct 2023 11:18:50 -0400 Subject: [PATCH 3/9] Fixing minor mistake --- notebooks/test.ipynb | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/notebooks/test.ipynb b/notebooks/test.ipynb index 0ac69da..b5e0d92 100644 --- a/notebooks/test.ipynb +++ b/notebooks/test.ipynb @@ -27,7 +27,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -110,8 +110,8 @@ "plt.subplots_adjust(wspace=0.4)\n", "\n", "print(\"Checking old - new\")\n", - "print(f\"Expected coverage: {np.abs((old_ecp - ecp).max())}\")\n", - "print(f\"Credibility: {np.abs((old_alpha - alpha).max())}\")" + "print(f\"Expected coverage: {np.abs((old_ecp - ecp)).max()}\")\n", + "print(f\"Credibility: {np.abs((old_alpha - alpha)).max()}\")" ] }, { @@ -125,14 +125,21 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 800/800 [01:19<00:00, 10.09it/s]\n" + " 0%| | 0/800 [00:00" ] From bc95cb8a2f761be65ed05b31d7995fdf0514b0f7 Mon Sep 17 00:00:00 2001 From: EnceladeCandy Date: Fri, 20 Oct 2023 11:19:16 -0400 Subject: [PATCH 4/9] Added tqdm to the requirements --- requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements.txt b/requirements.txt index 38f9a5f..f5a029b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,3 +7,4 @@ Sphinx==5.3.0 sphinx_autodoc_typehints==1.22 sphinx-pyproject==0.1.0 torch==1.13.1 +tqdm==4.65.0 From 1abfa927e8a11dda9968baa6dad920a703616d94 Mon Sep 17 00:00:00 2001 From: EnceladeCandy Date: Fri, 20 Oct 2023 18:01:46 -0400 Subject: [PATCH 5/9] Removed fast_histogram import --- src/tarp/drp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tarp/drp.py b/src/tarp/drp.py index a943a4f..1172cbb 100644 --- a/src/tarp/drp.py +++ b/src/tarp/drp.py @@ -2,7 +2,7 @@ import numpy as np from tqdm import tqdm -from fast_histogram import histogram1d + __all__ = ("get_drp_coverage",) From 6647105de451bdc78650534e35097dbdf72ffcb8 Mon Sep 17 00:00:00 2001 From: Pablo-Lemos Date: Tue, 21 Nov 2023 09:50:41 -0500 Subject: [PATCH 6/9] Reformatting --- src/tarp/drp.py | 181 +++++++++++++++++++++--------------------------- 1 file changed, 78 insertions(+), 103 deletions(-) diff --git a/src/tarp/drp.py b/src/tarp/drp.py index 1172cbb..5bba20e 100644 --- a/src/tarp/drp.py +++ b/src/tarp/drp.py @@ -1,22 +1,38 @@ from typing import Tuple, Union - import numpy as np from tqdm import tqdm +from warnings import deprecated __all__ = ("get_drp_coverage",) +@deprecated("get_drp_coverage", "get_tarp_coverage") +def get_drp_coverage( + samples: np.ndarray, + theta: np.ndarray, + references: Union[str, np.ndarray] = "random", + metric: str = "euclidean", +) -> Tuple[np.ndarray, np.ndarray]: + return get_tarp_coverage(samples=samples, + theta=theta, + references=references, + metric=metric, + norm=True, + bootstrap=False, + seed=None) + -def old_get_drp_coverage( +def _get_tarp_coverage_single( samples: np.ndarray, theta: np.ndarray, references: Union[str, np.ndarray] = "random", metric: str = "euclidean", - seed: Union[int, None] = None, - ) -> Tuple[np.ndarray, np.ndarray]: + norm: bool = True, + seed: Union[int, None] = None +) -> Tuple[np.ndarray, np.ndarray]: """ - Estimates coverage with the distance to random point method. + Estimates coverage with the TARP method a single time. Reference: `Lemos, Coogan et al 2023 `_ @@ -29,11 +45,12 @@ def old_get_drp_coverage( the parameter space. metric: the metric to use when computing the distance. Can be ``"euclidean"`` or ``"manhattan"``. + norm : whether to apply or not the normalization (Default = True) + seed: the seed to use for the random number generator. If ``None``, then no seed Returns: - Credibility values (``alpha``) and expected coverage probability (``ecp``). + Expected coverage probability (``ecp``) and credibility values (``alpha``) """ - np.random.seed(seed) # Check that shapes are correct @@ -57,10 +74,11 @@ def old_get_drp_coverage( theta = theta[np.newaxis, :, :] # Normalize - low = np.min(theta, axis=1, keepdims=True) - high = np.max(theta, axis=1, keepdims=True) - samples = (samples - low) / (high - low + 1e-10) - theta = (theta - low) / (high - low + 1e-10) + if norm: + low = np.min(theta, axis=1, keepdims=True) + high = np.max(theta, axis=1, keepdims=True) + samples = (samples - low) / (high - low + 1e-10) + theta = (theta - low) / (high - low + 1e-10) # Generate reference points if isinstance(references, str) and references == "random": @@ -100,18 +118,16 @@ def old_get_drp_coverage( return ecp, alpha[1:] -def get_drp_coverage( - samples: np.ndarray, +def _get_tarp_coverage_bootstrap(samples: np.ndarray, theta: np.ndarray, references: Union[str, np.ndarray] = "random", metric: str = "euclidean", - norm : bool = True, + norm: bool = True, seed: Union[int, None] = None - ) -> Tuple[np.ndarray, np.ndarray]: +) -> Tuple[np.ndarray, np.ndarray]: """ - Estimates coverage with the distance to random point method. - - Reference: `Lemos, Coogan et al 2023 `_ + Estimates uncertainties on the expected probability and credibility values calculated with the + _get_tarp_coverage_single function using the bootstrapping method Args: samples: the samples to compute the coverage of, with shape ``(n_samples, n_sims, n_dims)``. @@ -123,80 +139,50 @@ def get_drp_coverage( metric: the metric to use when computing the distance. Can be ``"euclidean"`` or ``"manhattan"``. norm : whether to apply or not the normalization (Default = True) + seed: the seed to use for the random number generator. If ``None``, then no seed Returns: - Expected coverage probability (``ecp``) and credibility values (``alpha``) + Mean and std of the expected coverage probability, and mean of the credibility """ - np.random.seed(seed) - - # Check that shapes are correct - if samples.ndim != 3: - raise ValueError("samples must be a 3D array") - - if theta.ndim != 2: - raise ValueError("theta must be a 2D array") - - num_samples = samples.shape[0] num_sims = samples.shape[1] - num_dims = samples.shape[2] - - if theta.shape[0] != num_sims: - raise ValueError("theta must have the same number of rows as samples") - - if theta.shape[1] != num_dims: - raise ValueError("theta must have the same number of columns as samples") - - # Reshape theta - theta = theta[np.newaxis, :, :] - # Normalize - if norm : - low = np.min(theta, axis=1, keepdims=True) - high = np.max(theta, axis=1, keepdims=True) - samples = (samples - low) / (high - low + 1e-10) - theta = (theta - low) / (high - low + 1e-10) - - # Generate reference points - if isinstance(references, str) and references == "random": - references = np.random.uniform(low=0, high=1, size=(num_sims, num_dims)) - else: - assert isinstance(references, np.ndarray) # to quiet pyright - if references.ndim != 2: - raise ValueError("references must be a 2D array") - - if references.shape[0] != num_sims: - raise ValueError("references must have the same number of rows as samples") + boot_ecp = np.empty(shape=(num_sims, num_sims//10)) + boot_alpha = np.empty(shape=(num_sims, num_sims//10)) + for i in tqdm(range(num_sims)): + idx_remove = np.random.randint(num_sims) + idx_add = np.random.randint(num_sims) - if references.shape[1] != num_dims: - raise ValueError( - "references must have the same number of columns as samples" - ) + # Replacing one simulation and its samples by another and its associated samples + samples[:, idx_remove, :] = samples[:, idx_add, :] + theta[idx_remove, :] = theta[idx_add, :] - # Compute distances - if metric == "euclidean": - samples_distances = np.sqrt( - np.sum((references[np.newaxis] - samples) ** 2, axis=-1) - ) - theta_distances = np.sqrt(np.sum((references - theta) ** 2, axis=-1)) - elif metric == "manhattan": - samples_distances = np.sum(np.abs(references[np.newaxis] - samples), axis=-1) - theta_distances = np.sum(np.abs(references - theta), axis=-1) - else: - raise ValueError("metric must be either 'euclidean' or 'manhattan'") + boot_ecp[i, :], boot_alpha[i, :] = _get_tarp_coverage_single(samples, + theta, + references=references, + metric=metric, + norm=norm, + seed=seed) + + # ecp_mean = boot_ecp.mean(axis=0) + # ecp_std = boot_ecp.std(axis=0) + # alpha_mean = boot_alpha.mean(axis=0) + # return ecp_mean, ecp_std, alpha_mean + return boot_ecp, boot_alpha - # Compute coverage - f = np.sum((samples_distances < theta_distances), axis=0) / num_samples - # Compute expected coverage - h, alpha = np.histogram(f, density=True, bins=num_sims // 10) - dx = alpha[1] - alpha[0] - ecp = np.cumsum(h) * dx - return ecp, alpha[1:] - -def bootstrapping(samples, theta, references = "random", metric = "euclidean", norm = True): +def get_tarp_coverage( + samples: np.ndarray, + theta: np.ndarray, + references: Union[str, np.ndarray] = "random", + metric: str = "euclidean", + norm: bool = False, + bootstrap: bool = False, + seed: Union[int, None] = None +) -> Tuple[np.ndarray, np.ndarray]: """ - Estimates uncertainties on the expected probability and credibility values calculated with the get_drp_coverage function - using the bootstrapping method + Estimates coverage with the TARP method. + + Reference: `Lemos, Coogan et al 2023 `_ Args: samples: the samples to compute the coverage of, with shape ``(n_samples, n_sims, n_dims)``. @@ -207,28 +193,17 @@ def bootstrapping(samples, theta, references = "random", metric = "euclidean", n the parameter space. metric: the metric to use when computing the distance. Can be ``"euclidean"`` or ``"manhattan"``. - norm : whether to apply or not the normalization (Default = True) + norm : whether to apply or not the normalization (Default = False) + bootstrap : whether to use bootstrap to estimate uncertainties (Default = False) + seed: the seed to use for the random number generator. If ``None``, then no seed Returns: - Mean and std of the expected coverage probability, and mean of the credibility + Expected coverage probability (``ecp``) and credibility values (``alpha``). + If bootstrap is True, both arrays have an extra dimension corresponding to the number of bootstrap iterations """ - num_samples = samples.shape[0] - num_sims = samples.shape[1] - num_dims = samples.shape[2] - - boot_ecp = np.empty(shape = (num_sims, num_sims//10)) - boot_alpha = np.empty(shape = (num_sims, num_sims//10)) - for i in tqdm(range(num_sims)): - idx_remove = np.random.randint(num_sims) - idx_add = np.random.randint(num_sims) - - # Replacing one simulation and its samples by another and its associated samples - samples[:, idx_remove, :] = samples[:, idx_add, :] - theta[idx_remove, :] = theta[idx_add, :] + if bootstrap: + ecp, alpha = _get_tarp_coverage_bootstrap(samples, theta, references, metric, norm, seed) + else: + ecp, alpha = _get_tarp_coverage_single(samples, theta, references, metric, norm, seed) + return ecp, alpha - boot_ecp[i, :], boot_alpha[i, :] = get_drp_coverage(samples, theta, references = references, metric = metric, norm = norm) - - ecp_mean = boot_ecp.mean(axis = 0) - ecp_std = boot_ecp.std(axis = 0) - alpha_mean = boot_alpha.mean(axis = 0) - return ecp_mean, ecp_std, alpha_mean From 0f5624d1faf2b5d35bcaa4d0db660c0a9922a48b Mon Sep 17 00:00:00 2001 From: Pablo-Lemos Date: Tue, 21 Nov 2023 10:07:53 -0500 Subject: [PATCH 7/9] Fix to deprecation decorator --- requirements.txt | 1 + src/tarp/__init__.py | 2 +- src/tarp/drp.py | 11 ++++++++--- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/requirements.txt b/requirements.txt index f5a029b..4f144f3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -8,3 +8,4 @@ sphinx_autodoc_typehints==1.22 sphinx-pyproject==0.1.0 torch==1.13.1 tqdm==4.65.0 +deprecation==2.1.0 \ No newline at end of file diff --git a/src/tarp/__init__.py b/src/tarp/__init__.py index f68186c..8f97047 100644 --- a/src/tarp/__init__.py +++ b/src/tarp/__init__.py @@ -1 +1 @@ -from .drp import * +from .drp import get_drp_coverage, get_tarp_coverage diff --git a/src/tarp/drp.py b/src/tarp/drp.py index 5bba20e..b61945c 100644 --- a/src/tarp/drp.py +++ b/src/tarp/drp.py @@ -1,13 +1,18 @@ from typing import Tuple, Union import numpy as np from tqdm import tqdm -from warnings import deprecated +import deprecation -__all__ = ("get_drp_coverage",) +__all__ = ("get_tarp_coverage", "get_drp_coverage") -@deprecated("get_drp_coverage", "get_tarp_coverage") +@deprecation.deprecated( + deprecated_in="0.1.0", + removed_in="0.2.0", + current_version="0.1.0", + details="Use get_tarp_coverage instead", +) def get_drp_coverage( samples: np.ndarray, theta: np.ndarray, From eabc573ec9570105372d60b746d998a1bb3bc227 Mon Sep 17 00:00:00 2001 From: Pablo-Lemos Date: Tue, 21 Nov 2023 10:20:38 -0500 Subject: [PATCH 8/9] Changes to unit tests --- tests/test_drp.py | 49 ++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 42 insertions(+), 7 deletions(-) diff --git a/tests/test_drp.py b/tests/test_drp.py index fc9e39f..06b5f3f 100644 --- a/tests/test_drp.py +++ b/tests/test_drp.py @@ -1,11 +1,11 @@ +import unittest import numpy as np +from tarp import get_tarp_coverage -from tarp import get_drp_coverage - -def test(): - num_samples = 10 - num_sims = 10 +def get_test_data(): + num_samples = 100 + num_sims = 100 num_dims = 5 theta = np.random.uniform(low=-5, high=5, size=(num_sims, num_dims)) log_sigma = np.random.uniform(low=-5, high=-1, size=(num_sims, num_dims)) @@ -14,8 +14,43 @@ def test(): loc=theta, scale=sigma, size=(num_samples, num_sims, num_dims) ) theta = np.random.normal(loc=theta, scale=sigma, size=(num_sims, num_dims)) - get_drp_coverage(samples, theta, references="random", metric="euclidean") + return samples, theta + + +class TarpTest(unittest.TestCase): + def test_single(self): + samples, theta = get_test_data() + ecp, alpha = get_tarp_coverage(samples, + theta, + references="random", + metric="euclidean", + norm=False, + bootstrap=False) + self.assertAlmostEqual(np.max(np.abs(ecp - alpha)), 0., delta=0.25) + + def test_norm(self): + samples, theta = get_test_data() + ecp, alpha = get_tarp_coverage(samples, + theta, + references="random", + metric="euclidean", + norm=True, + bootstrap=False) + self.assertAlmostEqual(np.max(np.abs(ecp - alpha)), 0., delta=0.25) + + def test_bootstrap(self): + samples, theta = get_test_data() + ecp, alpha = get_tarp_coverage(samples, + theta, + references="random", + metric="euclidean", + norm=False, + bootstrap=True) + ecp_mean = np.mean(ecp, axis=0) + ecp_std = np.std(ecp, axis=0) + alpha = np.mean(alpha, axis=0) + self.assertAlmostEqual(np.max(np.abs(ecp_mean - alpha)/ecp_std), 0., delta=10.) if __name__ == "__main__": - test() + unittest.main() From ad55f6d20a07dbaf6194bfb77ec08bb84945c7f1 Mon Sep 17 00:00:00 2001 From: Pablo-Lemos Date: Tue, 21 Nov 2023 10:20:52 -0500 Subject: [PATCH 9/9] Updated version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 518791f..6c1dd10 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "tarp" -version = "0.0.2" +version = "0.1.0" authors = [ {name = "Adam Coogan", email = "dr.adam.coogan@gmail.com"}, {name = "Pablo Lemos", email = "plemos91@gmail.com"},