From b622feb6ad690d745000471fd1230a73cbc57ea8 Mon Sep 17 00:00:00 2001 From: lgrcia <20612771+lgrcia@users.noreply.github.com> Date: Tue, 30 Jul 2024 20:30:08 -0400 Subject: [PATCH] feat: fast injection-recovery --- docs/index.md | 1 + .../notebooks/tutorials/GP_optimization.ipynb | 37 ++- docs/notebooks/tutorials/analytical-ir.ipynb | 314 ++++++++++++++++++ docs/notebooks/tutorials/exocomet.ipynb | 1 - docs/notebooks/tutorials/tess_search.ipynb | 48 ++- nuance/core.py | 23 +- pyproject.toml | 2 +- 7 files changed, 385 insertions(+), 41 deletions(-) create mode 100644 docs/notebooks/tutorials/analytical-ir.ipynb diff --git a/docs/index.md b/docs/index.md index 5b9c5c6..bda974e 100644 --- a/docs/index.md +++ b/docs/index.md @@ -71,6 +71,7 @@ citation.md :caption: Tutorials notebooks/tutorials/GP_optimization.ipynb +notebooks/tutorials/analytical-ir.ipynb notebooks/tutorials/ground_based.ipynb notebooks/tutorials/tess_search.ipynb notebooks/tutorials/exocomet.ipynb diff --git a/docs/notebooks/tutorials/GP_optimization.ipynb b/docs/notebooks/tutorials/GP_optimization.ipynb index a606d59..fb1dc96 100644 --- a/docs/notebooks/tutorials/GP_optimization.ipynb +++ b/docs/notebooks/tutorials/GP_optimization.ipynb @@ -1,39 +1,40 @@ { "cells": [ { - "cell_type": "code", - "execution_count": 1, - "metadata": { - "tags": [ - "hide-input" - ] - }, - "outputs": [], + "cell_type": "markdown", + "metadata": {}, "source": [ - "import os\n", + "# GP optimization\n", "\n", - "os.environ[\"XLA_FLAGS\"] = f\"--xla_force_host_platform_device_count={os.cpu_count()}\"\n", + "nuance requires a Gaussian Process (GP) of the light curve to be built and optimized before searching for transits.\n", "\n", - "import jax\n", "\n", - "jax.config.update(\"jax_enable_x64\", True)" + "In practice, any `tinygp.GaussianProcess` object can be provided. Here is an example of how to build and optimize a custom GP on the light curve of the active star [TOI 451](https://ui.adsabs.harvard.edu/abs/2021AJ....161...65N/abstract).\n", + "\n", + "```{note}\n", + "This tutorial requires the `lightkurve` package to access the data\n", + "```\n", + "\n", + "In order to run this tutorial on all available CPUs, we set the `XLA_FLAGS` env variable to" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 1, "metadata": {}, + "outputs": [], "source": [ - "# GP optimization" + "import os\n", + "import jax\n", + "\n", + "jax.config.update(\"jax_enable_x64\", True)\n", + "os.environ[\"XLA_FLAGS\"] = f\"--xla_force_host_platform_device_count={os.cpu_count()}\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "nuance requires a Gaussian Process (GP) of the light curve to be built and optimized before searching for transits.\n", - "\n", - "In practice, any `tinygp.GaussianProcess` object can be provided. Here is an example of how to build and optimize a custom GP on the light curve of the active star [TOI 451](https://ui.adsabs.harvard.edu/abs/2021AJ....161...65N/abstract).\n", - "\n", "## Loading data\n", "\n", "As in previous tutorials, we will download light curves of TOI 451 using the `lightkurve` package." diff --git a/docs/notebooks/tutorials/analytical-ir.ipynb b/docs/notebooks/tutorials/analytical-ir.ipynb new file mode 100644 index 0000000..be22977 --- /dev/null +++ b/docs/notebooks/tutorials/analytical-ir.ipynb @@ -0,0 +1,314 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fast injection-recovery\n", + "\n", + "Because nuance provides a way to perform the full-fledged modeling of transits in the presence of stellar variability and systematic signals, it can be used to analytically estimate the detectability of planetary signals in a given light curve.\n", + "\n", + "Traditionally, this task is done empirically by injecting and blindly recovering synthetic transits. Although nothing can formerly replace this process, we can estimate the signal-to-noise ratio of an injected signal given our model of the light curve (Gaussian process + linear systematic model).\n", + "\n", + "```{note}\n", + "This tutorial requires the `lightkurve` package to access the data\n", + "```\n", + "\n", + "In order to run this tutorial on all available CPUs, we set the `XLA_FLAGS` env variable to" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import jax\n", + "\n", + "jax.config.update(\"jax_enable_x64\", True)\n", + "os.environ[\"XLA_FLAGS\"] = f\"--xla_force_host_platform_device_count={os.cpu_count()}\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The dataset\n", + "\n", + "We will run this tutorial on the light curve of [TOI 451](https://ui.adsabs.harvard.edu/abs/2021AJ....161...65N/abstract), for which a GP was optimized in [the previous tutorial](./GP_optimization.ipynb). Let's download the data" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAEiCAYAAABkykQ1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAClN0lEQVR4nO2deXwb9Zn/PzptOYcijLEjSJSAICJBQAxtCC1sISluILAUdhfacqUtW7qllLuk8ONaCi1n6ZZCdznC1UK7UJaGgH/g5QglhTYmoDiYYpooFMWJMcJOYtnWMb8//Hsm3xnNSDPS6PTzfr3ySiKPR3N8j+d+bJIkSWAYhmEYhmEYhikCe6UvgGEYhmEYhmGY2ocVC4ZhGIZhGIZhioYVC4ZhGIZhGIZhioYVC4ZhGIZhGIZhioYVC4ZhGIZhGIZhioYVC4ZhGIZhGIZhioYVC4ZhGIZhGIZhioYVC4ZhGIZhGIZhisZZ6QuoVTKZDGKxGKZNmwabzVbpy2EYhmEYhmEYy5EkCTt37oTf74fdntsnwYpFgcRiMcyaNavSl8EwDMMwDMMwJeejjz7Cfvvtl/MYViwKZNq0aQAmHvL06dMrfDUMwzAMwzAMYz3Dw8OYNWuWLPvmghWLAqHwp+nTp7NiwTAMwzAMw9Q1RkL/K5q8/dprr+Hkk0+G3++HzWbDM888k/d3XnnlFbS3t6OhoQHBYBCrVq0yfU5JknDttddi5syZ8Hg8WLp0KT744ANrbophGIZhGIZhJiEVVSx2796Nww47DPfcc4+h4zdv3oyTTjoJxx13HDZs2ICLL74Y3/72t9HZ2WnqnLfeeit+/vOf47777sObb76JKVOmoKOjA6Ojo0XfE8MwDMMwDMNMRmySJEmVvghgwr3y+9//HqeeeqruMT/84Q/x3HPPYePGjfJnZ555Jj777DO88MILhs4pSRL8fj8uu+wyXH755QCAoaEhtLa2YtWqVTjzzDMNXe/w8DC8Xi+GhoY4FIphGIZhGIapS8zIvDWVY7Fu3TosXbpU8VlHRwcuvvhiw+fYvHkz+vv7Fefxer1YtGgR1q1bZ1ixYBiGYRiGqXbS6TSSyWSlL4OpYlwuFxwOhyXnqinFor+/H62trYrPWltbMTw8jEQiAY/HY+gc9Hvq89DPtBgbG8PY2Jj8/+HhYTOXzjAMwzAMUzYkSUJ/fz8+++yzSl8KUwPMmDEDbW1tRfdmqynFopLccsstuOGGGyp9GQzDMAzDMHkhpWKfffZBU1MTN/NlNJEkCSMjI9ixYwcAYObMmUWdr6YUi7a2Nmzfvl3x2fbt2zF9+nRD3go6B/2e+PC2b9+Oww8/XPf3Vq5ciUsvvVT+P9X0ZbKJRqPo6+tDMBhEIBCo9OUwDMMwzKQinU7LSkVzc3OlL4epckiG3rFjB/bZZ5+iwqIqWhXKLIsXL0ZXV5fisxdffBGLFy82fI65c+eira1NcZ7h4WG8+eabOc/T0NAg96zg3hW56evrQyKRQF9fX6UvhWEYhmEmHZRT0dTUVOErYWoFGivF5uNUVLHYtWsXNmzYgA0bNgCYSKzesGEDtm7dCmDCS3DOOefIx19wwQX429/+hiuvvBK9vb345S9/id/+9re45JJLDJ/TZrPh4osvxk033YRnn30WkUgE55xzDvx+f86KVIxxgsEgPB4PgsFgpS+FYRiGYSYtHP7EGMWqsVLRUKi//OUvOO644+T/U6jRueeei1WrVmHbtm2yQgBMeBuee+45XHLJJbj77rux33774f7770dHR4fhcwLAlVdeid27d+Nf//Vf8dlnn+GLX/wiXnjhBTQ2NpbydktCNYYdBQKBqrkWhmEYhmEYpjxUTR+LWqNa+lh0dXXJFbGWLFlSsetgGIZhGKY6GB0dxebNmzF37tyaNJpOdiRJwne+8x3893//N+LxON5++21cfPHFOPzww/Gzn/2sJN+Za8yYkXlrKseCyYbDjhiGYRiGqRfOO+882Gw22Gw2uN1uBINB3HjjjUilUgCAV155Rf653W6H1+vFwoULceWVV2Lbtm1Z5xseHsbVV1+NUCiExsZGtLW1YenSpXj66aeRy7b+zjvv4JRTTsE+++yDxsZGzJkzB2eccYZcPUnklltugcPhwG233Zb1s1WrVimud7/99sOKFSs0z0O88MILWLVqFVavXo1t27bhkEMOMfLoqoKaqgrF7EEMgWJPBcMwDMMw9cJXvvIVPPTQQxgbG8OaNWvwve99Dy6XCytXrpSPef/99zF9+nQMDw+ju7sbt956Kx544AG88sorCIfDACCHuw8NDeGmm27C5z73OTidTrz66qu48sorcfzxx2PGjBlZ3z8wMIAlS5Zg+fLl6OzsxIwZM7BlyxY8++yz2L17d9bxDz74IK688ko8+OCDuOKKK7J+Pn36dLz//vvIZDJ45513sGLFCsRiMXR2dmre/4cffoiZM2fi6KOPLvAJVg72WNQovb29SCQS6O3trfSlMAzDMAzDWEZDQwPa2toQCATw3e9+F0uXLsWzzz6rOGafffZBW1sbDjroIJx55pn44x//iJaWFnz3u9+Vj/nRj36ELVu24M0338S5556L+fPn46CDDsL555+PDRs2YOrUqZrf/8c//hFDQ0O4//77sXDhQsydOxfHHXcc7rrrLsydO1dx7KuvvopEIoEbb7wRw8PDeOONN7LOZ7PZ0NbWBr/fj2XLluGiiy7CSy+9hEQikXXseeedh+9///vYunUrbDYb5syZo3mNNpsNzzzzjOKzGTNmyPnEjzzyCKZOnYoPPvhA/vm//du/IRQKYWRkRPOcVsCKRY2STqcVfzMMwzAMw9QjHo8H4+PjeY+54IIL8Mc//hE7duxAJpPBE088gW984xvw+/1Zx0+dOhVOp3bgTltbG1KpFH7/+9/nDJcCgAceeABf+9rX4HK58LWvfQ0PPPCAofvJZDJyeJfI3XffjRtvvBH77bcftm3bhj//+c95z6fFOeecgxNPPBHf+MY3kEql8Nxzz+H+++/H448/XtIyxKxY1CjUvCSTySAajVb4ahiGYRiGqVei0Si6urrKLm9IkoSXXnoJnZ2dOP744/MeHwqFAABbtmzBJ598gng8Ln9mhqOOOgo/+tGP8PWvfx177703li1bhttuuy2rSfPw8DD++7//G2eddRYA4KyzzsJvf/tb7Nq1S/fcH3zwAe677z4ceeSRmDZtWtbPvV4vpk2bBofDgba2NrS0tJi+fuJXv/oVtm3bhosuugjf+ta3cP311+OII44o+HxGYMWiRhEnSjka0VVqUWEYhmEYprKUu/Ht6tWrMXXqVDQ2NmLZsmU444wzcP311+f9PfIu2Gy2vJ4G4uabb8bUqVPlP9Tm4Mc//jH6+/tx3333YcGCBbjvvvsQCoUQiUTk3/3Nb36DAw44AIcddhgA4PDDD0cgEMCTTz6p+I6hoSFMnToVTU1NmDdvHlpbW/H4448bur5i8Pl8eOCBB3DvvffigAMOwFVXXVXy72TFokYJBAIIh8OWV4TSUyA4p4NhGIZhJiflrkB53HHHYcOGDfjggw+QSCTw8MMPY8qUKXl/77333gMAzJkzBy0tLZgxY0ZeueWCCy6QGytv2LBBETbV3NyMf/7nf8btt9+O9957D36/H7fffrv88wceeAA9PT1wOp3yn02bNuHBBx9UfMe0adOwYcMGbNy4Ebt378Zrr72Ggw46yMwjyUJLedLqmv3aa6/B4XBg27ZtmonnVsOKRY3S3d2NSCQCn89naTM6o1YJ9mBUH/xOGIZhmFIQCASwZMmSsjW/nTJlCoLBIGbPnq2bB6EmkUjgP//zP3HssceipaUFdrsdZ555Jh5//HHEYrGs43ft2oVUKoW99toLwWBQ/qP3fW63GwcccIAsnEciEfzlL3/BK6+8olBMXnnlFaxbt06h0NjtdgSDQey///7weDwFPJFsWlpaFOV1P/jgg6yk7DfeeAM//elP8Yc//AFTp07FhRdeaMl354IVixqFJonWZCkGPatEKBSCx+ORQ7DK7RZl8sPvhGEYhpks7NixA/39/fjggw/wxBNP4Atf+AI++eQT3HvvvfIxP/7xjzFr1iwsWrQIjzzyCDZt2oQPPvgADz74IBYuXKibC7F69WqcddZZWL16Nf7617/i/fffx+233441a9bgH//xHwFMeCs+//nP49hjj8Uhhxwi/zn22GPxuc99zlASdzEcf/zx+MUvfoG3334bf/nLX3DBBRfA5XLJP9+5cyfOPvtsXHTRRVi2bBkef/xxPPnkk/jv//7vkl4X97GoUagDotfrtfS8gUBA0yKh/jwYDMp9NJjqgN8JwzAMM1mYN28ebDYbpk6div333x8nnHACLr30UrS1tcnH7LXXXvjTn/6En/zkJ7jpppsQjUbh8/kQDodx22236cpQ8+fPR1NTEy677DJ89NFHaGhowIEHHoj7778fZ599NsbHx/HYY4/hhz/8oebvn3766bjjjjtw8803l+TeAeCOO+7AihUrcMwxx8Dv9+Puu+/G+vXr5Z//4Ac/wJQpU+RrCIfDuPnmm/Gd73wHixcvxr777luS67JJRrNbGAVm2puXgq6uLiQSCXg8nqpokCc27CuXq5RhGIZhmGxGR0exefNmzJ07F42NjZW+HKYGyDVmzMi8HApVo5QqkarQOH0Ow5l8cE4HwzAMwzAirFgwCtQKQj7hkX7u8/ng8Xjg8/lY2Kwg5RT2WZlkGIZhGEaEFYsapVTlX9WeEBIeI5GIprBKP4/H41iyZAni8TgLmxXESmE/n5JS7vKDDMMwDMNUN6xY1CjpdFrxt1UMDg4ikUhgcHAQABRCoyis6nkq2HNRWawU9vMpKeUuP8gwDMMwTHXDikWNkslk5L+7u7stO69YxpYSsv1+f5awSh6TgYEBhaeCPRflQc+bYKWwzx6J2oJzXhiGUcP1eRijWDVWWLGoA0gJsALqOOn3+7PCnLSE1WQyiWg0miWEslBaWqzOb9ASStkjUVtwzgvDMAT1M1A3TGMYPWisiL0wCoH7WNQJfX19BQmA6jKx7e3tACaUFa/Xi1QqhVQqhWg0qjh/KBRCJBKRv1sUQLn0bOmxumeFKJSK74zfZe3AfUwYNTx/Jy8OhwMzZszAjh07AABNTU2w2WwVviqmGpEkCSMjI9ixYwdmzJgBh8NR1PlYsahR/H6/LPyPj48XLExoCZQUDjU0NASPx6MpcAYCAQwODiIWiyGRSGDt2rWytpvJZJBOp9Hb28ubWYnQa2RIaAkUuYQMPaFUT+Fgqo98Y4KZfPD8ndxQozhSLhgmFzNmzFA0FywUVixqlObmZgwMDGBkZAShUEhz0zBirdISKElpcTqdSCQSsNlsip9Ho1H09vYimUzKnw0NDWWdO51Oo6uri61lFUBLoMglZOgJpWwFr2/Yol2/RKNRpFIpuFwunr+TFJvNhpkzZ2KfffZR7NcMo8blchXtqSBYsahRRMFezzNgxFqlFiij0Sji8bisXADZCT19fX1Zi1RjYyPS6bTic4fDwdayCkEKAVXnCgaDBSkJbAWvL0iR8Pl8iMfjSKVSSCaTPEfrEFqnPR4Pv9tJjsPhsExoZJh8sGJRx2gJl/k2GFJYEomE4vNIJILe3l60tLQglUpl/V46nVbkXQBAS0sLYrEYfD6fNTfEmGZgYEAWHI0mYpNHCph4h/F4nC3aNYbaE0H/J0WC5rfL5eIiC3WKniGBvVQMw5QSVixqlFAoJAt/oVBI/ly9aQQCAXR1dclN7gDk3EyoL4bNZoPT6UQqlZI9FslkUvZiuFwuhXcik8koqtF4vV7E43EAkP9mygd5qwoRHEWPFL1v0aItjjH6GQsp1YXaW6keD+Sx4PdWv+h5GznvgmGYUsKKRY1CmwYJefSZ1qYRDAYVFZy0NpPu7m7EYjG5aoQkSVnhThSvOTAwAGBCeaDcikwmo/BkUEI5x+eXHi0LpPjsaZwY9VoFg0FNjwV9lziWALCQUoWo5556PDD1TS7ln9dlhmFKCSsWNYxayAsEAvD5fEgkEorwIzF5V52ETaFSevkUIpIkob29XfaAiAnboiLicrng8/lYkCkTojJJlbr8fj98Ph8ikQi2bt0qv6tcCoAojHR0dOh+FxEMBuVO7RzuVl2I1moOfZl8qHuaiMo/500xDFNKuEFejdLd3a3IZxgfH5cTr4Hs8CN1szNSShKJhKxUABNeCfJa2Gw22O12+f/UPI+a3+nhdDq583aJERvaic0Ixc7pYtlgwufz6XZoztdgTawyEw6HEQgENMcbd4CuLqxonMfvtLYQ1wRuVsowTDlhxaJGEZUBYCI3gqySRjYRUcjw+/3weDxwuVyQJEn2WkiShEwmo/BidHV1AUBOC7XP55MFUEocZ4EkP2aENz1hUeyc7vV6AUxU7CL6+/tlhVL9u/ROSUlVI1YiE0OuxPEmKqysVFaeXCVHrRhvTHUihjwBMFy4gWEYplhYsahRSIAkzIYfkUAYDofR3NwMYCKe3m7XHxLUDK+vry9LsRGvgyoRsefCHGaEN1IC6J3T77W3t2P58uVob2/H+Pg4gAnPk8vlAjCRC0OoBU3yOpCSqoYS++lvINsTpg6VYioLJeI7nc6sstJmFEC2etcetC5EIhE27DDsdcwDPx/rYMWiRiFlgEgmk9i2bZthQUEUCGkDisViWbWuKQxKRF2KVn0dwISCkUql4PP5WCAxiFHhLRqNyoodJVaLv9fd3Y3Vq1fD7XbLFYCAiXfi9/vhcrlkRUNcTOk4h8Oh6WmisZGrHjpdi9/vR19fHy/SFUZvTJlVANUKJFNdaAlF4ntlww7DXsfc8POxDk7erlGoao+IJEmaQoRe8iZ97na7ZWVB3aMiVzK3HmI1qXg8jiVLlpg+x2TEaFKlWihUewvE3Irly5ejq6tL8U7o35FIBA6HI8tDkU6n5XP09PQgEonA7/cjFArpVpMRx9iSJUvkBH+uFlU+tOZ5IBDA4OAgIpEIBgcH0dzcLBds0OpHw9QmJBT19vZmjQG9oh2czD+54GpgueHnYx2sWNQRXq8XxxxzjPx/dVMsEvKoARoJmKIHohBFIhepVArRaJQ3sCJQCwJapUNFa0tjYyNGR0cBAGvXrlUIkFq5OaIyKhYEAPaETsViMTQ3NyvitsV3SuOJusDzIl1+9PoTiAn9FJoITBRZYOWvPqD5lkqlFD2LqGrb1q1b5fnIfSwmJ/kMV5Nd4eRqadbBikWNEgqFsHHjRlkRcLlcCqUC0G+SJioVwES+BoXUqAVLALJVmwiHw4rv1sNutyuETaYw1IKAegGMRqMYGxsDALjdbkUVKPq3zWbTfF92u13Rod3lciGdTsPhcCAUCmWVqgUgW0Z7enqQyWSy8n0AXqQrgajMiUKC3++Xe9SMjY3BbrcjlUqhqamJSwXXCTTf1q5dK89XcZ0X5zAr/ZODfIqC+udaCudkVDbEMvzcRLQwWLGoUQKBgGLjSCaT6O7uRnt7u3yMKDAODAxohk9R2VCC+iCIiEqF1+tFT09PXqVCjMNPJpPo7OxEKBTiCVoAakFAa0Mgz4KoVIjovS+HwyErCGTJpvyL3t5ehbdDFEbFjuyxWAzhcJiFlQojeq/IS9nb2wun0wmXy4VkMim/s0wmI48VdWlqpnbRm/8EzWcOT61/1DkDel5uMiZqGSbU0Q6TAXoutB9Opnu3ClYs6ggKVxEFvFQqpVAUenp6sGDBAvkYir8G9vSwsNvtiupBIvk2LkJURgBMusWpEPSsQ2rrv9aGoOVp0vJS0DsmxU/0XNntdjgcjqyO60QymZTfv3hem83GHooqQfRSkscwmUzCbrfLCmNLSwsGBgZkz5Taw8HvsXYh7xQwMV8dDgemTp2KoaEheW7zOly/iPNYVBS0vBHUTBeYWDfE4gyUI6eOdpgMqPfTyXTvVsGKRQ3T0tKiUBoaGxvlCaEOdyIymQx6e3uRTqezjhF7WFgFCTMAT9B8qBd/yoURQ5Mof4He88aNG9HY2KgQKAitd0nvWEtx1KoApkbre+j36HoBsHeqQojChLg5ZjIZOBwONDU1Kd5fJpNReClZ6Kxt2tvbFXk06XRaLjtN6wGHvtUv4h6iruJGHmjKeRQ9leq9WSuPrx4wYkBRFz2op/svFzbJaklykjA8PAyv14uhoSFMnz697N9PdeiN4HK54HA4MDo6mmXF1ou9twp1qBWjjzq2k9zQIi6XC01NTZqeI3UuTCG4XC44nU7dZO5wOJylkNrtdjQ0NCiu1+PxcLhFmenu7kYsFoPf70d7eztWr15t+hw8X2ufaDQqhzc6HA64XC65mAMwsU7Mnz+fBac6JJfgTF4IWpsno5dS/QzUTMZnYhQzMi97LGoUsTxoLuUgHA4D2CMgSpIkN8ErhYdCDZW5FHM/GG0onEh0Q6sRw5HU6IWvmaGlpQXt7e2ykKoeW+RBUX+vurcJW0XLi9jbhKo/5cPr9WJkZETRe4Y309pF7eGUJAnpdDprvqbTaXk/6Onp4XdeR+QKSeWk/fzPgCumWUNFG+S99tprOPnkk+H3+2Gz2fDMM8/k/Z1XXnkF7e3taGhoQDAYxKpVq7KOueeeezBnzhw0NjZi0aJFeOuttxQ//9KXviTHmtOfCy64wKK7Kg/BYBAulyunUuH3++UkbzUNDQ0lVyqIWCyGaDTKnS0NQk3NQqEQwuEwPB6P4uder1fz94p5n6RsDgwMYO3atbKQqj5nJpMxpMBwQnDhFDJP1E2dEolE3tA2tYIaCoWMXyRTccRx0t3djUgkgmQyiUwmo0jUz4UVxgimNlA3uezp6ZFz9To7O+X1pp736XyNPo02qWVyU1HFYvfu3TjssMNwzz33GDp+8+bNOOmkk3Dcccdhw4YNuPjii/Htb38bnZ2d8jFPPvkkLr30Ulx33XXo7u7GYYcdho6ODuzYsUNxrvPPPx/btm2T/9x6662W3lupCQQCCIVCupuHzWZDLBbDSy+9lBVOQxZmEibp71ISiUQQiUS4s6UB8i1+U6ZMQTgclhN0i4Wsm0BujwiQnZSvh9vtrusNqpQU0gFWayOUJAkulwuNjY26v0d9THJ1Suf3WJ309vbKpZ/VeU/M5MToXI1GowqlkpL6AWWzxck27/PtvYwxKhoKtWzZMixbtszw8ffddx/mzp2LO+64AwBw8MEH4/XXX8ddd92Fjo4OAMCdd96J888/HytWrJB/57nnnsODDz6Iq666Sj5XU1MT2traLLyb8pNL8CBBUYytFbHZbPLCUm6rFYfJ6KNXW1wkFoth9+7dmsn5Xq8Xu3fvNtVVudi8DC2GhobYrVwghYYsaFVzo6pQejgcjryd0vk9VjfJZFJuimk2Z04sC87UPnq9KNRFQNSyg81mk9cbdbPFep/3nFdhPRX1WJhl3bp1WLp0qeKzjo4OrFu3DgAwPj6O9evXK46x2+1YunSpfAzx+OOPY++998YhhxyClStXYmRkpPQ3YDFut1v+t8PhMGW9rmTOPofJ6COWkl29erWsVIhWZ6fTqetVGBoaMqVUlAq/389u5QIxazWjQg5mDAQulwsul0sOf8r1rvg9Vg49C3R3d7dCYdQzIIlo7Q/z588v/iKZqkFrrvb19SlC5Hp7e+WGqoQkSfJ6Q+tPKBSaFPO+EA8xk5uaSt7u7+9Ha2ur4rPW1lYMDw8jkUggHo8jnU5rHiPmGXz9619HIBCA3+/Hu+++ix/+8Id4//338fTTT+t+99jYmGIyDg8PW3RXhSMKl5lMRqEsqK2XZMkqdRWoXKgrDjHZ6PWkSKfTcjJ3Lgt0NeByueRkfbYAlZ5cGyJV8VB/dswxx8iWzN7eXoRCId0qXtyjpHLoeYvE0CeHw4HW1lZFmVk15NEQsdvt8tjh91sfqOdqNBpFKpWSlUryUJFsQHKC3+9XWO4ByBUK63mM0PNxuVy6DWgZ89SUYmEV//qv/yr/OxwOY+bMmViyZAk+/PBDHHDAAZq/c8stt+CGG24o1yUaQtws1MpCW1ubYvOhalCVTNZLpVLc3yAPgUAAW7duzRIGRWWiEsqhVv8KPZqamtDV1SUvzLxQlxYxdIGa4eXqxD5lyhS56hfR29vL76YK0QuLE+djOp3GwMAAWlpadBULUamg9YNy7eo91GUyI5YGpzAo6ltDZYfpOLGXDQDN7tPqkui1vqaTN8fj8WR1JOd5UTg1FQrV1taG7du3Kz7bvn07pk+fDo/Hg7333hsOh0PzmFz5FIsWLQKQ2/K3cuVKDA0NyX8++uijIu7EGnK5v7WamJVTqdCqSCNJkmaFKkYZ8kANrfSQJEm3MlSpaG5uNnzs0NCQwrXMrubSQGMGgCJ0IV/uWCwWy1ofksmkojKM1vdMpiTOaqe9vV2xBiSTScOKv3ptdrvd/H4nAel0Gn19fRgYGAAw4a2gqpGi4Wp8fBw+n08u6iCGQ9FaHovFamZNz7V+aYWOcehn8dSUYrF48WJ5IyVefPFFLF68GMDEAnnEEUcojslkMujq6pKP0WLDhg0AgJkzZ+oe09DQgOnTpyv+VBq/31/pS9DN68hlUWdBJRtR+DayoA0NDeWs9mM1RpsxiojJgLxQWw9VBdq4cSNWr16NwcFBLFmypOAcJoq/VsOKYeXQeva0fu7cuVNxLJVtz4fawKQ2BDC1h96eGgqFFGFQRtbgdDqN/v5+ABMGJQqR7Orq0lU4qpl861cqlUJvby+6u7vR1dWFwcHBMl9h/VHRUKhdu3YpXvbmzZuxYcMG7LXXXpg9ezZWrlyJjz/+GI888ggA4IILLsAvfvELXHnllfjmN7+J//3f/8Vvf/tbPPfcc/I5Lr30Upx77rk48sgj8fnPfx4/+9nPsHv3brlK1Icffohf//rXOPHEE9Hc3Ix3330Xl1xyCY499lgceuih5X0ARdLc3GzYSlWq0BmjXhAKz2hpaWFXowY+nw+JRAI+n083HEqNkYRNozidTkuTvr1eryIZkN+zdVA4As09mtuxWAzt7e26eTpG0BoD3Fircmg9e7WFGZgQGnfv3l3QOu/1ejE+Ps7vt4bRqwbV19eHBQsWAICcT9XS0oKBgQGk02l0dnaiqakpa68Rw+SAPYalVColV+CsFfTWL3U4KP1bHf7FmMdSxYKSg43yl7/8Bccdd5z8/0svvRQAcO6552LVqlXYtm0btm7dKv987ty5eO6553DJJZfg7rvvxn777Yf7779fMdDPOOMMDAwM4Nprr0V/fz8OP/xwvPDCC3JCt9vtxksvvSQrHLNmzcLpp5+Oa665ptjbLytUCabSUEJxKBTCxo0bdTc2EoIGBgbkcne8kU0gdk2Ox+OIRqN5lQqrsbqS1NDQEKLRKOdYlAASLO12OzweDyRJwujoKBobG9HV1WVauBSVSvXvijHV9ZzEWa1oJeNqFW9Ip9MFrxlDQ0NcdrbGEQ1ThKhsUP4VMLEHO51OuVKU1rix2WzyurBp06ay3Ucp0DNs6RllnU4nbDYbUqmUvIcx5jCtWJx33nm45557MGXKFMXnW7Zswdlnn421a9caPteXvvSlnJugVlftL33pS3j77bdznvfCCy/EhRdeqPmzWbNm4dVXXzV8jdVKNbit/X6/XP0HMB4uwxZsJeK7dLvdVaEwipArXZIkU0JrT0+PYlNjC1BpoD4khXqwcimVJJyQFY+TvCtLqdZ9ir/nd1ubUPgj/S1WO/L5fAohOplM5kz0BybWei0FlspTE7VsNKICCA6HQ9HLiZ4bhYbW2n1VA6ZzLN555x0ceuihir4QDz/8MA477DDsvffell4co081WPsHBgawZs0adHZ24qWXXjL0O7QwcZ7FHigHwe/3G7I6GvUKOp3WOCQdDgfcbrdcWcwo5E4HUDPxuNWG1jxpaWkBMDEOxOdrRUK/+hxkATXjiWasJxqNorOzE2NjY6bmoBl4ftYuwWAQLpcLqVQK3d3diEQicnNMypcQGRgYKGhOq/MPaiH/Sk/WaG9vRzgc1jSWkVKVTCYNySgszygxLXm89dZb+NGPfoQvfelLuOyyy9DX14fnn39e7njNlAetag6lRl2ulr6bGu/kgxayrq4utmIL0P0b9VQY9RoYCW/Sqm+vRuzerJdTo1UC1+v1yo0njVi0atn6VSrUsdNiXDBZ2VpaWtDe3p5V2KIQ1BXJqIKMJElwuVxoaWlRlBJmygOVxQSUBTPU1lYzqNfz3t5eDA4O1kUZ0clGIBBQVGwS0Vqz0+l0Qfk4sVhMDmcOBAI1kX9Fz4VyTKgDeUtLi6EcVXWpXa25wXmjSkybPlwuF2677TZcddVV+MlPfoJnnnkG//f//l9WKiqA2i1ZavKVsTRCJBJBIpGQhV4xJnQy09PTU5Hv1VMqKHfGKFqb1NDQkKyUGLFo1YL1q9yoK2ppbYQk/Lvdbku+LxfUhC3fOyILu14JW8YcwWBQVihEQTGfUpGrcqBa4KSStWbnIFtrqwP13M3lkXA4HKY8FuKx4npOXbqrWZimNRSAogO50cI3wWBQzmnVmxtc+VCJacUimUzisssuw09/+lOsXLkSixcvxmmnnYY1a9aU4vqYHJRzMjscDmzbtq2oc4jCJ/270NKY9UYlGxeqocQ+KxG7mpIQohZIeHHORr1x64XBFJrwny+sRjReJJNJjI2NKbrU6kEWdqNKJZObQCCAhoYGAPnD0sR3ambNdjgcBZURZYNA9eByueByueD3+3N6JEKhkCJU1uPxaBqTHA4HvF6v3DuJzl9LazStoRRCStD9iKjnFhU1EMe21r3XgoJVTkyHQh155JEYGRnBK6+8gqOOOgqSJOHWW2/Faaedhm9+85v45S9/WYrrZHTwer2WVhCiiZROpxXhLYW629Wo3fi1tEBZjehaNRKSVC6srhBls9k0u5oCULiPOak/PwsWLFCEzFF4UqEJ/2qFtqenR/EO1CGXmUzGkBIcDAbR09ODTCbDXskioNA3r9crJ5UCyKn4i+/HTLiL2+1Ge3u7vC4ZDYuqhXCYyYAYLkdeTD3U68X4+DhaW1uzrPhitbGhoSEsX75c/lktha6KlRcJrUa06vmSTqfl3j4ul0sOAWNyY9pjceSRR2LDhg046qijAEwIDT/84Q+xbt06vPbaa5ZfIKNPKcqSSpIkKxGSJMHv9+uGxTgcDtOJhJlMBpIkyZvfZJ6kopBdLUpFKZAkSQ6JEb0S7KEwTyAQkMNbGhsbkUwmi/YkiohCaXd3N1avXo2mpqas4yKRSM7QF9HCzl7JwiFhiMIKSem3Opne4XDA5/Ohq6tLbrxoNCyKrbWVhTy/agXezBhJp9N5Q4PUJYlryVMlNv4UvXJqI4nWM6PQKafTyWPcIKYViwceeCCr1CwALFy4EOvXr7fkohhj6E3ocDhc8DnVEy0ej+vmcrjdbixYsAAul8uUglGqZn21hihYW1HRp5oRE/VJCGGBxDyi5Y2UUSvnkxiTLwq14XA4a47nUy5YcSwedY4ElQE1+87zCZnpdBrbt29XVBmrpe7K9YaRvBU6hhTBeDwujxeHw6E5RswoG+qqgul0Gt3d3fL/a2V+q3u/xGIxjI6Oore3N+sZif8nRYpkIva8Gsd0KJTYsE6L2bNnF3wxjDncbnfOWtRWfYdemIXoAqUmWkYTogBkxTwSteRitYLBwcGSNMRzuVxwOp0lHyNGKXQDmmzjQQ+9pphWhtE1NzfL/6Y6736/X1F1RgyR1Or0O9nfUzGon2F7ezva29uzugSbxYgikk6nZUGR319lMVJliI6hvIdUKiWHQOmtB0YVUq0qf8CeqlAtLS2KMLlqnft6a6Zenw4RUijoOcRiMTQ3N1fV/VUrpj0Wc+bMwdy5c3X/MOVDTxi1ssGa3ndQeJRYLWFgYAAul8uwVUQvRKKWXKzFoFce0CqSySQSiYRl/SyKQSw1bLZ6zGQZD7nQ2yDtdrulYXRidbL29nZZueju7pYtlKLAQd1pgez3xO/NPOpnRlbp5uZm05XagInxYcabrM6v4opPlcGIN4COoURssSx4sYgh0WrE6mE9PT1YvXo1Nm3aVHVzXWvNLDayoprur5oxrVi8/fbb6O7ulv+8+eabuO+++3DQQQfhd7/7XSmukdEhVynBUpNMJtHT05M1cZPJJBobGw1dG7kWJ2N1ILEzqhVhUHa7XVehS6VSiqQ7NerYWWBibFnZiEuSpJzl+nIxGcZDPsRnJo4XM9XEjCj8mUxGMRdJ6Y3FYrJFUpzbYtUn9Xui+c0hBMahZ+h2uxUCWyQS0fXw5sJosj0hvk/RaNTb28sKRhnJFyYqeggGBwfL4pVWrx8ej0ceW6K3q1rQ2meKrb7Ia5kxTJsyDzvssKzPjjzySPj9ftx222047bTTLLkwJj/t7e3o7+8vWalSscFZOp3O+h7x/36/X+7mmUgkDC105LFQu30nQ3UgquDh8XjkZwxMxLUWUpVJ/W68Xi+Gh4flMoFibCxht9s1GwX5/X7E4/GSjSu9zUfPnT4ZxkMuRCU0FAopEhFF7HY7JEnSDXcwGgahp/yRgKmG3qf6PdH85uRt49AzXL16NQBlNb5SeTYp7MVms8Fms8leKPUYIAt1JBKB3+9He3t7Sa6nHFRr6I5RqFJbOZvkiuuHzWaDz+fD2NgYMplMVY4HqlZmpdLFa5kxLDNJzps3D3/+85+tOh1jkEKb1hkJjxkfH0dHRwc6OjrkJG0Ru90Oj8eDcDiMeDxuKqnQZrPJAslktEgHg0E5LlYU4K0q9ToyMoKTTjoJ4XAY4+PjCqHE5XIhHA6joaEByWQS8Xhc4Z3o7+8vWZUqsnZr9bLg0BltxDKSPT09uoIEVVwrBq/XK89FLSUilUopvj8cDusKZuyxKBwxCbcYHA5H3rWexgxV60smk4hEInLPknA4jFAopLBQl0rJKRf1stak02lLPctGkSQJsVgMmUwGLpcL8Xi8arxZtKcA1q49tda/o5KY9lgMDw8r/i9JErZt24brr78eBx54oGUXxhijUA06lUrB4/FoavNerxfj4+OKSUSWtDVr1siby4IFCxRChVHrCVnGenp60Nvbi1AohCVLlhR0H7WK2B+gFBtDKpVCV1dXliAI7FEqSXhMpVJoa2uThQW1p8Lj8WB0dDRLaM1nIQeU1lD6XgqfE3tZkBWUrovZA1neyDpYKlwuF4455hj5/1qKhfiu7Xa7ogOvGvZYFEY0GkU8Hper+xVjlc5X4IM8zVrnz2QycDgcindLvUkqGYZrBdXce8OINyUUCqGvr09zfS8V6oRukhNSqZS8htN+XkkvECmNVnpzbDabXB2zq6urZj1d5cK0RDNjxgz4fD75z1577YX58+dj3bp1uPfee0txjUwOilkYtTYch8OBY445Rje+kyxo5L0ga3MgEDCcJEyWMbKO1brVqFhKISxKkqQoHSkiWuoo4S8Wi6GxsRFAdoIbJYCrY2wzmUzeuP2ZM2fKCgh9D0G9LIh4PM7lZ3UotVIBTIwFMWROq3+FCIU99vX1aSb5TkZPpBWQYBSJRDA4OFjQOWgO+3y+nDlcsVjMkPC1du1aRCIRZDIZeDyeqgt7sZpKJq1reVOi0Sg6OzvlfkCBQECeVy6XS167bTZb1vu2ynClNiCNjIxgyZIlinL01bCfUzSAlQqXJElyGGA9eLpKjWmPxcsvv6z4v91uR0tLC4LBYFVUn5lskBBmRSUoh8OB+fPnA9C2mtAiS7He6twIsoCbmdDq7tu1HvtqBvEZWo3dbkdDQ4P8bLU6rW7atEnxGYU/iQIsKRTpdFrTM6EWdtULeq6Qia1bt2J8fFzO6WABVJu+vj5LlAqyMqvD4tTvq7m52dA8TqfT8vzVKo8plqIW/8/kxufzyWtCLBYrqBoUjRc9b4RRWlpashqx1kNoW75yrkbKvZYKLW+KGA7Z29ur8Fa4XC5FTxt1R+lSGySAiXVEXA8qCZXGttqTIz7HSt9jtWNaE/iHf/iHUlwHUwSBQACDg4OaQpzorqSJplf3XkwUFJPDRAGBEo4ByAmlYvImYE7JsdvtioW7kgt6uaH3ZrViQYofAHkDUqNXSlANhczlW6QpvM3MYi4KK2Io3GRSLo0gCprF0N/fD0AZ0qD1vszMXxpH6rUAUJZ7FNcRJjdi6JhWcYVyMjAwkBXKVkv1/PXWklyhUGKxhEoIkFrFKmgNoDU2mUzK3mL1+j42Nlbya6Rx2dnZKa8hLpcLHR0dJf9uNVrvuBAjp1EKUfQnG4YUi2effdbwCU855ZSCL4YpDLEbr5rZs2cjEAgoGizlSszNJVSIMfk0aT0eT5ZiUAzigl/vAmau91YIpFBQwyJ6ly6XSxGSQrHzWouuVmMkGi/k0VIv2B6PB0uWLEFXV5cpAVgrlweYXMqlEajpVbGYtVy6XC40NTUpFEB1XpbX61UYHPTWgnJYTesFUZFMp9MFrxGUo6GOxfd6vTkbcoqV6ZLJJJqamrLmtdirpJrXZ721hIw6kUgEmzZtwvz58zUNaOW+L709j9YAcW0Wk+5FyjHX0ul00d4wq/Z3vXfsdDrlRn7j4+OGjWn5UBtcRepdZjGKIcXi1FNPNXQyCplgyotamBc9EpFIRE64I6jMqFYJWWCP4ClavglxISGLjjiZSDHw+XwYGBjQ/Q6iqalJkQwlWmtIUK1XAdOKOE1SBKjcH8UGi1Yseof0Xun/NC5IkPB6vZg9e7YiaVes9OV2uxXeMbViYMayTuFPPp8PfX19GBwclMOhqjmxcrJAQqna0ED5NjS+xsfH5ffudrsVc1kcD8VWq5pMWJHsbrfb5WaHDodDVhAdDgdmz56d04CktoCrlRB1+BuF5lSjMJVrLSGFLZ1Oy/kstCbR75YatSCaz6hC+RJWKQ96HbaNQIVHqICH2T4rVhmQtN4xGb/GxsbQ0NBQsFyqfj52ux2ZTAbpdFoziVt9T5NV0TCU1UOJtvn+sFJRGcQJpdWJV6v/REtLCxoaGhS/J5JMJuF0OnN6I+jnWslmzc3NeZUKYGLTokRFSkwT76uekz/dbrfm57ma3YmIi14sFlOUbAUglwIWK1BRch3F7Xs8HnzlK1/B8uXLFRWBgIn4apfLBbvdrggLIOsZJe8BwJo1a0xZVqlzq/pvWpA5iXsPVOrTLGLlHpvNZtqFryd8ioJnMBiUBWGay7QOiAIyKxbGoeTTYhD35WQyKSsH6XQafX19htYXPWifp/UZ0O99UmnEtURMyI5Go1l7Hq1B5Swiod479fY8WgMWLFig2LeJQt/n9OnTixprYjnjgYGBrCTzXFi1v2u9YxqjmUymqDBScd1yOBzynk3n7enpydnc18qyxpUsKGAWQx6LvfbaCx988AGam5vxzW9+E3fffTemTZtW6mtjDCJakY1aMmKxmLywUslBqklNScU0OUhg9fl8SKVSWUlaosVAnEhmrSpUP53uqV4ao+lZLfTCEbSem5ZlSf1/8R3kspD4fD45QTeVSqG7u1u21InKQTwel2NmxYZZdH30d29vb8EWNEoy9/l8ssciGo0qvCv1MAasQC83SgubzYb29nY0NzcrxoM4l3OFMhjZCG02G3p7e+VQGSopHAwG0d3drdjQWbHID70bt9tdVIgJvQf1OUiAzNVvyKgFW9wfKNykWg1A9FwpHEwM5RSFavV9lMPaLIYXU7Unre8Sw7S03k+h84v2ILLEm0XM9chkMllJ5rmem9X7uxj+C0y825aWFsuaCKfT6SzjOSkYes19rfK+i/dWCxEchhSL8fFxDA0Nobm5GQ8//DB++tOfsmJRJdDip5Wwpd4kHA6HXEHCZrPJky0ejyuUCfXkIGUhkUhoNsTSOj4YDGaFYBmlFiaOGfRcvn6/35CV3263K/pM6EGLF+XAAHs2pFAoJC9MAwMDCqGSzisKgurqHuI92O12Sxoz0cJP19PS0pIVD55rLEwWN3MhSYg078W5uXbtWlmQyNeI0UhYG5UQpnNKkiRvvOqxWokmXrWAOIbFdbZYQqGQoqAHhTnmq0JnVEAdHx9XjMtq7kMkXqe6wARFWmiV0C1HrpeR8Cf19agx0k/I6/VieHhY95hiBG86p91uRzAYlPeZTCZTtjVarVQAExEV7e3tchf7fIghnkZRV19UY5XyJBp6qlWBFzGkWCxevBinnnoqjjjiCEiShIsuukjXNf/ggw9aeoFMbnJtFOIiQjH41JFSkiQ514Li3NWlZcW8CaPasnoi9fb2GgqJIqgbNVnR1ZbWWhQi1VaLaDSKjRs3Gt7EqctpLkjh6+rq0rUYidY5cXOlfAcxwY3ijimxUX0PZPXu6uoyZRVSJ4Zu3749Sxh1uVxy+ULRkgdoC2L1poiqKaQbu1bvAtFDlk9RUY83o5bsTCajyNEh67k6V4uZQBzDhRpj1NZmSZLQ19eHJUuWIB6PI5FIYHx8PK9SYQYyLNRaqKp6DKfTad3qT+XK9TL6PXSc2+1WzOW2tra8idR63vFicizUqOc4NdC0Yo3Ot/+Lgjd5dhOJBNauXWv4O8zMO3pubW1tZenpYiQSoZowZEZ67LHHcOKJJ2LXrl2w2WwYGhpCPB7X/MOUF4rp8/v98t9aMZMDAwNYs2ZNllW6o6ND3nzEyaleEMLhsOlNJBAIoKOjQ26qlw86jhqpUWJgV1cXent7qzaONx/qnAE9d7Ye6m6nwIQyQDH0fr9fUWZPC3JRk5BHwns4HEZ7ezuWLFmC+fPnZ42ddDqdlZxJ90PjJh6Pa8b9aqEWkrXyspqamtDR0YGGhoashktqQazWBBujUDxtd3d3QRv/rl27smJyczVKy4fZa/B4PHC5XLJXo9BGb/UKvRufzwePxwOfz1dwSKH6d2w2myJMleYI/TscDsvJ+UbQ608VCoVkT0U1x37TeqeHOpeQKEeulxmDWSAQyFIqAGNNDvU6pVulVLhcLjmXjwiFQpat0flyFSgvSZ1jmqv6mRozc488sJTbCJQ2B6LW8g4NeSxaW1vxk5/8BAAwd+5cPProo2hubi7phTHGoIFGE5osF1qdJ/UmjpbFpBKVecRGP+RiTKVSSCQSWSVTaxkzi7nH45E3E7/fL+dGUOiQx+NRzEUx30as0iFW6+rr68vKXaANjkKTUqkUJEmSlT0tq5N6jBTSpJH6X4hjc2hoCNFoVA7HERtyqS03tbLQmqXYsJh0Oo1NmzbJiiEwEb5CAiXFdZci96GlpUWuUEZjIhaL1X23ZjOI3gMq1WxVzX1JkjA4OJhVqS8YDCoUASPkyuvZunWrIsyoWj2HtE6oQ2JcLhecTqfsea2ENdhoFSF1YQ6z7N692/TviN6MfJ4NaqQojuHe3l5ZuSCFoNDnq7UXiFBIWSn6VmghGsTo3miti0Qi8r1X43woB6Yb5G3evLkU18EUgTihyBINIKujMbnZ1eEJWgKaXp5FIZuH0QZP6g0slUopEurqZZIaTcAFJp4BWfkHBgYUz5H+rX4n5DkUPYj0PqmEr9oLoRZkHQ4Hli1bBkBpVRPRGyNm0IsNFi1TWvdRz4gNuooR/sWwNno3kUhE0+hAoWdEMTHX9L5EL4WexbTcVEtIpVopF8NNgfy9JvJBa4P4vmlOkYFhbGws73vOtVapr6+ajD7ie966davms6S9hcI4K9HEUSvElOYpCadAYUYbEaNjSSzeIua65VuDYrFYVr8dccwVGw6ltacRYtEDq5vNGkFUnAi693rfq/QwrVgw1QdVlgAgLwpi2TzC6lrRRjEaIkfWMbHHhvoe6gEzQgN5n9TCILli1UnWQO53RT8jTxAtfmLOAqC0yBgV5mkcmrUaqRVdUSBIJBKQJEkOGxEV5WoQEEuBaCiguWAWm80mJ9nTGBGFE9ETJs671tbWghuyiYmMYkNOAFXj4a6WvBxxTtGzEq3CRtYHrffV2NgIm82myJciw4Tb7TbtCTMag099M6rFUiu+Z717VeclJJNJRT4XUHpFVKuKEM1TqpJodafnXBWgJElCJBKRDQFm1gLaq8jrDUzcz+DgYE5vgxFy7Wm055hRKgpJ1NZCLGbT19eXNe/U42myYJO4DmBBDA8PywLi9OnTK3otYgnJeDyuCJ2phvADrYoNWlCoE3lWaOOsN48FgCzBSw/RgiQunNQBmdz5Zp+P3oZJ16U1doxuskbftxq1h03vHGLtfOr6XU+Iz6/QMpDAnmR9sRu7+P6KDa/QYvny5QCQFXYiesAqSbV4LESMVq3RQquR4fLlyxUlm0XhORwOF239FqFQSVH5rdScFO+ZGgL6/X7s3r1bU1EjpUmtPJGwKM7Dct6T0b2hWqE9i+ZXZ2enPAa1qkpqYWSeihUxy9FtXAv1/ahLG9fT/mRG5mWPRR2gtkLR39SwJlfFJyuE0XzQsfkqnoyOjio2PbfbLVuVKuGmLgW0+eWz6qsXZwAKhYusP+Pj43JoUyHvRB372t7erquM5rP2igpuIYihXSJerxcjIyNyWBCdv1grWLUi5k2Z9f6QIkJlNSmpks6rFbpWjPKihmLV1aWUq6V5arWE0olrKXmlCqnQ09vbm5VcTe+AQlDJIyUWeShkbGlB103oVVgqB6KnjxSJWCyGcDisqUDTs5YkSeER3rhxIyKRiKLoSCmsz+oQU/o3efjFeameo2bGipljzfTK0UPsR0X/J4x6Co14FukYl8sl90IqVCErtDqW2LFdHPek2I6Pj2tWNSQjsOiFz/dcqtEoogcXF68DqPKCVq14rWoKVGGpt7fXVCWDfJUZchEIBLBgwYKc9ezVE7uaYnatwoxS0dfXJ7+XQCCAE088ER6PB+l0Wg4RK6Y/gNn3SUK8njBP5zO7uKu7xrrdbvk7vF4vdu7cqWjqJVagq+dKdIUIfg0NDQiHw7qVfERo3ViwYIFlORC0rvT39ys+r5Yci2pBnHtUjc2MYEPzXquzMJ2X5uHo6CiWL18uGwyoWp+ZylC5oOt2uVzo6OiomNAjVgYS18WNGzfmFZbFylF0P6QMk5JudUVCcQyI/6a1T1wX1WtkoVUFtcKqqFqcw+HAgQceWNC4cDqdCIfDiudOeXzidxvd041Uk6JjqEhJc3NzwZXvigncEStY0nsUO91TZUvRQ0y/Q38bGVvFyF/lxpBUMjw8bPgPU36oFNmCBQsUnxsp92ZmsBZSOo4Ul7Vr1yISiZi2jIoLfrWWMywFeuV11e8gFArB4/EU1CfA7PvMJ8yLpY/NxAWrF/WhoSFZKBoaGsoaM9QBvpLW0VJjZD5qKZU0ZsyOC6sVNPGdUUljZg/i3CtEEBct2R6PB42NjQD2lIYdHx/X/V1akwFtQbMQSPCt5BpNClNDQ0NWX49cgmNjYyMikYiuIi9JUtY6aUVpUa1SwGqPhXgNhSAqJA6HAy0tLVnrhigEU9K4WVKpFAYHBxXPPZlMyuWUw+GwKaXTaHnVVColC+eRSAQjIyOmrjufYc6M4S6ZTGJsbAwul0veA10ul6z8RyIRuN1u+XiHwyErQvk872JBj1rY8wyFQs2YMSNLY9ajWlzekxExvEXcsNSTU+yyLR5v5PxmN0GjyYJarsje3l50dHRUTcKlFRipkJVMJuUFTb3gqN9BMaEdZn83XwK/+nwUW2tlqI3L5crbDKrWiUajGBsby3tcrvLRVHaYXPVagj3NKxIkSMjUerZGk8jtdruiy7vNZqv5OVsKrAjJEsMlSVGghFRqYJfJZLK8RVrv3eFwFBUCQ/1KqmGNVjeSyxfmku++SWkjxNyLYu5Xay0nqLEsUej6SfdO8zcejxs6F/VWMep9ttlsmsdu27ZNLn8MKGUTs6E96uO1ysuaTchesGCBZs4R5S6azT/LZDJoaGhAc3OzZtNYMdcnnU7L/89l2FHn+lR6fhnBkGLx8ssvy//esmULrrrqKpx33nlYvHgxAGDdunV4+OGHccstt5TmKhnDGNmwci1oehQa30eLvCRJORdwrYWfJmUlempYifjsjFqGSZFXl/CrJGaFoVAolLVBFot6I6EFtxYWW6PkapTmdDqxzz776G74drtdfhZ0jF4PCbFCmKgAagliRt9hJpPBBx98IMes22y2ivUIqAXE/jHUcNIo4lxQC9MOhwPz58/XfObq964uNSxituxtNVTCoXVKLENaTOle2rdIiRC9iUb3JDOFL+j8emuAUUONeBx1GDdSZliszkfvP9840FPc6HMqRysqn2YNhuroilQqlZWAb8az4/f7NT3DDodDVtiff/550/sXzS/1M86lwOeaN4WMt0pjyM/zD//wD/KfRx55BHfeeSduueUWnHLKKTjllFNwyy234Pbbb8dDDz1U6utlDGJ1F8hC4/vIpSlOTpfLlWUF0iKdTqO7u1tekAcHB7F69Wp0d3ebvv5KIj47igPWg7phF5M7UW60xlp3d7ccXlBsVaN81ELMqVHUTabUpFIptLe3Y/ny5Zp5C+KzFjuz56KlpUV22QPGwi5y5XCMjo6iqalJvp5qjQsuZadco98pxlz7fD5D66IYPUDx22QRplAQt9utK6zRmkzhcoB+Po8ZgVzdL6Pc6D3bYpQKwm63y8Kf2MHcqAJldP+k43JZ3hcsWKAbQUL7htfrRUNDg+JnTqdT9lhQw1k6nkLq6J4oFJfC23ft2mXoPnORSqUUIWXqUNx881E8nrwVYgJ+PsRnRlUPg8GgYq91uVxwu90YHBxEZ2en6RA0r9cr58mIJeEpFExUKrxer/wOqECN3n2L/clqAdPSy7p163DkkUdmfX7kkUfirbfesuSimOKxOtGnkPwKPVpaWgy73cX4SdECW0vQwkCbhZ5Q5vf75TA16qpdSO5EudEaa8W+IwrzyKeE1Us3dsJIfHNnZye6u7s1vV/iJkkKCHXBVm/a9N7i8bihZG+RfCEH1dw8jahEMqT6O8XnEovFDK2LVMWI5oaYBAog75wQcyxIwSgWr9crKyput7siBiB6tpFIBNFoVBHPng+xAhQhzqVMJqOoTmgk/l9EvX92d3dj9erVWLt2rTwvxTh6tUArCsUbN27UFXjJODA8PJw1R4PBoBxa29LSIudlUgEHrXtSJ7IXgyRJ8ndoeXDyzUfxuecz0GlBz9Bms2Hbtm1YvXo1BgcHsxQwmk9kFDOaBmC32zE0NCSvqfRsW1tb5UIsouFgfHwcwWBQfmfUR0VE9GBRta1ayDU1rVjMmjUL//Vf/5X1+f33349Zs2ZZclFM8eRTBMxa6wpZTEXEzSuf0JnPWl9rFWYCgQCcTqe8MekJD7FYTLYUUdWnWggf0Rpr4jvyer0FWVsGBwd1Lankrq6VZ5QPEjSMxAgnk0mFIEnYbLasAg6E1qatTh6tJS+ZFeRaI2l97O7uttSrof7OQCBgej1zOBzo6OhAR0eHXBGHBBYSGHPNCXEsiIJLMQwNDWHTpk0IBoOKUq+lQmv/EvPRIpGIKU9Fa2tr1mdk2bfCUqzeP8XiFGJlKFrvclWBMmJFp5wXYI8XPBAIZBXg0NrX6ZkajSzIhbimiP/OtR75fD7Fu1W/a/KGm82noGcqhkyRpxCYmEPJZDJLyTTqtRC9xWLDWbHyk7j3a3Xs1vp/IpHIKtlb7ZjuY3HXXXfh9NNPx/PPP49FixYBAN566y188MEHeOqppyy/QKYw8sXDlzshWqyfTtYPvRAZ8XPqYUAdhKuhq2shGO1KTT+vpf4MWmOtvb1djhkfGRmRFSujkPCsR65Qj1qE7tWs612svy9Jku4z0cpT0utrUQqqsQ9NrjVSXXTCqnVStM5S/XsKmzASMkhNQylvhQREEliM5HDRWPD5fJY2y6OqQmLfDCsRrdy0ltK4ikajRSky8XhcM7dIkqQsAZa8hcXkDVGfF3pW4rwrpjCF1j04nU75Oo3kK5KC4/F4iu5Onclk4PV6Zes8kes6tm/fLpdpVediDA4OFrxWUkEDUbGgpqzAnjmUTqez+vCYwW63a+YWjo+PK5LpgWzPr7rymFYRj2r0/qoxbaI68cQT8de//hUnn3wyPv30U3z66ac4+eST8de//hUnnnhiKa6RKQFWhjYZhSz3mUzGcNz9zp07kU6nkclk0NLSUnXCiVHo3o1SD/0ZaIyl02nLBVaKd64XSAjzer2GPAdkSRXLR+YS5Ix4HGtJmS0GI95aCrUgK6eVz0arlj2th/nq8EuSpKibT6Ezfr/f8HpOY0FvjTEa+qGHum+GVeQKlTFrxbXZbIr363a7NQVVUfBzOBwKb2ExlmMKUyz2WavRuge1McHoOjA+Pm5JNb+hoSH5nDTvtK6D3q9aIBdlFS1hXyuMTQ+HwyF7YTwej5xnIUL5F3pzMZcXR6ysp352YjPJdDqNTZs25VQi9RLAa4GCfN+zZs3CzTffjKeffhpPP/00fvzjHxcUBvXaa6/h5JNPht/vh81mwzPPPJP3d1555RW0t7ejoaEBwWAQq1atyjrmnnvuwZw5c9DY2IhFixZl5X6Mjo7ie9/7HpqbmzF16lScfvrp2L59u+nrr2WKDW0qFLHXgREBSlRCai23QoTiZ/NtJDabrS7yBsQOo1aVmhVJJpPYuHEjOjs70dnZWfNKBgkaRjZzm82GdDotCzkLFixAOBxGPB4v6jnkqkBWrPWZQnYqDZVuFMutakGGALEpo1XkamiaL3xH3ScgmUzC6XSivb29oLh/LawIiSvFnBQFTMrnaGlpQWdnp9w/wGiDNAoVovdrJGyKjqWk20LWaLVSWyplnvYZsdu60euisZ5OpzUVlUKUITHMl0Lw1Mq9KBvo5Reqx6bL5cLUqVPzfr8o0JMXwO12y7lGYqELUoj1xkSuPChxTGkhKkFaeSukrJK8oL7fXEne1URBK8jatWtx1lln4eijj8bHH38MAHj00Ufx+uuvmzrP7t27cdhhh+Gee+4xdPzmzZtx0kkn4bjjjsOGDRtw8cUX49vf/jY6OzvlY5588klceumluO6669Dd3Y3DDjsMHR0d2LFjh3zMJZdcgj/84Q/43e9+h1dffRWxWAynnXaaqWtnCoMUmvb2dlOWBqJWBUh1FQs1VJXjkEMO0RUQKlHFplBEq2ypIOGgkpVorMRo/wpJkhQCprpzb6Hksp6pDS9mhYtqKZssPp9kMpkzwVisxmKloq/X0DQXZCXVeu6FXlsgEEA4HJa7VZPnQy9R12yzMKvnpGgMowReMck2mUyabpBmBppz4+PjijXazLqsnqdmFVbRy6KH1+uVx4vR84uJ75T07nA4FHs07VEzZ840dc0ipJCRkiEKyaJsID5fse+Kw+FQ3H8ymTSkFM6cORMejwcOh0N+j2J+i1bzV7V8YrPZDM2BXIpFrtAyMb+EQv0aGhpqLqcUKECxeOqpp9DR0QGPx4Pu7m55IxwaGsLNN99s6lzLli3DTTfdhK9+9auGjr/vvvswd+5c3HHHHTj44INx4YUX4p/+6Z9w1113ycfceeedOP/887FixQrMnz8f9913H5qamvDggw/K1/nAAw/gzjvvxPHHH48jjjgCDz30EN544w386U9/MnX9jDGKFYjVSV+1SL6u1AsWLJDjTvWeUyWq2BQKWeK8Xq/l7n4REoiMliysZnL1r9DD4XBkde4tBaIbH7DGql0J1BbiWCyma12nTs5mOgabwUziNllJGxsbZWVATMotBlEwVxsCCjH+EKX2umqtg1Y3ztTygKjHkJl1WUxQ7uzsNB0ims8i7vf7MT4+Dp/PJ1ciNLIWiu9K7MRNpaipB0MikSjYSNDU1IRUKoWenh75PaVSqbwVxOiZAfnvX4+BgQEEg0FZaXY4HIrwQa31c/78+Yq9urGxMauClFlyXTvlTKmru4lhWTabrSYqRZpO3r7ppptw33334ZxzzsETTzwhf/6FL3wBN910k6UXp2bdunVYunSp4rOOjg5cfPHFACYsCevXr8fKlSvln9vtdixduhTr1q0DAKxfvx7JZFJxnlAohNmzZ2PdunU46qijSnoPkxH1wkshMkZK2FH5QlrsyH1Zaw231Imiq1evzjomX0J9LTUKJMvP+Pg4GhsbkUgk4PF48jZnMguFRIpjq166tBtBbIhX7P1SIqmRn9ntdlMlKKslFEpLKKIyjkB5Gy2aSTi22+1oaGiAz+dDX1+fZUUscgnDJMzQ2mtm3lI9fwAKy34hTVb1IMs3gKzGY2JRg2LQsoarvQBm1uXBwUEkEgmMjo4WJCDng8ZTIpFQ9BUx0jQX2LM303kGBgbQ0dEBYM/7K7RDu9azFKsz6eXkqBsemukITohzHJiQI9Tfp35G6u8NBoPYunWrIWXQ4XBg6tSp8j0baTRJ3yGOJZpD9LvUxbza9zbTisX777+PY489Nutzr9eLzz77zIpr0qW/vz+rLFxrayuGh4flEp3pdFrzGFqA+vv74Xa7MWPGjKxj+vv7db97bGxMEaZAjWOY/IiThVx8Ri01kiQpJiT9uxqrzGgRjUbR09MjV8cgV71aUItEIrL1Um+DMtv5upJoLZDBYBCbNm2y7DvIWq9uQFgryhchdl+m7rTFnqsQ4S2XZ0kcq4V4oKqlGIH4bB0Oh0I5KrcyatTzSD1dgD2d5q261mAwKK9PaoppLEe/KypsViv9ouAnCo1k1f3ggw8KFoJzoe6SbGRdFvcBwHxFI6OI1Yyoq7p6LdRbI8T70FqH6OdaRjEjqIVrcf6JHlBRgRgYGJBzPSRJgt/vR3NzM/r7+4syUFHIkc/ny1nhS3xWgHJOaFXgItRyy86dO/Ne0/j4OAYHB5FKpWR5VSukMJcSVi2Y9me3tbVpLoivv/469t9/f0suqhq55ZZb4PV65T/cs8M4xSSK620MWs1kqhGxssPQ0JCcE0D3JS6o+XpX1FKYj947L3ZDFYVa2pREN3alihIUAwlc/f39ppUKtSegmHA5PWVMrDoUDofleu9mqMaKUxTqoA6nE6H+IlY3e4tGo4aF3nQ6LQvmhFWKcyAQwIknnojly5dbWnFGDPGi7uA+n89QuJ7ZdU6dzCpJEiKRiKVKRWNjIzwej2Z3cSPXW0iYo1kcDgeam5sVe4pYapYwskZQcrxW2E2hMf9qYyztBQ6HQ5FvRPkXYv6M6NnQqpZkBpfLJYcc5avwJeZ3qEszH3LIIbrzRn19Rq43nU7L90zrq9Z8qYWcC9OKxfnnn48f/OAHePPNN2Gz2RCLxfD444/j8ssvx3e/+91SXKNMW1tbVhLh9u3bMX36dHg8Huy9995wOByax7S1tcnnGB8fz/KuiMdosXLlSgwNDcl/PvroI2tuapJBC5ZWjLbZ7r+1kGsQDAZlYbixsVGOV6X7dzgcCIfDhjbcWsqxIMQqPFZct1oxiUQi6O3trbnQOBESuo1ulmLTLq2wjEJzLbRi/sn6SwmVALKqldhstrxejGrxWFB4T2Njoyyg0HPXGj9k/bW6CEFfX5/h6jF0ffRurcirALIFYlFJLSaHxm63IxaLydWTAJhq+qle57SaFYrXTs+nlLlc1GuiqalJMb+KWd/8fj/C4XDB16R+R+l0Gn19fYqxo7UOGFkjKDmecv6i0ahcga+5ubmg69brJk79T9TKu9YY9Pv9BTXzpGdNlcTylWmm8UXKsBaRSASrV68uifcpVwPY5uZmy7/PakyHQl111VXIZDJYsmQJRkZGcOyxx6KhoQGXX345vv/975fiGmUWL16MNWvWKD578cUXsXjxYgATcXNHHHEEurq6cOqppwKYWJS7urpw4YUXAgCOOOIIuFwudHV14fTTTwcwEd61detW+TxaNDQ0FJ24w+xxp3Z3d2dt1maa8agtjFbH71qFGAIwOjoqu0/FsCgAstCWi1oM8xE321I1XxMbZdUiZoRum82Gjo4ORbiA2qVvZCzpoQ4HlSQJvb29Cquw2lths9lyKkUUslYNTJkyRfYcimFQ6XRaM3eLQku8Xq/h3C4ja1EwGMTGjRt1hRLx2kiIsjoUUi3Ai+uxWD3H6XSaWptFD+3y5cuzwknyoRVGqW5WCECuYkSUKrxIZGhoSBY01SFYwWAw693T/5uamjA0NKQIn4nFYjnDr/OhNefGx8fR2tqaM7zH6DhSjw+a95FIBOFwuKgmckB2udVYLIZ4PC6HhALKLtkejwfNzc0FeX/i8Tja29sRCATQ1dUlNwDUCykSxxw11dOj2CaCWlDunJacVAt7nWmzhM1mw9VXX41PP/0UGzduxJ/+9CcMDAzg3//9301/+a5du7BhwwZs2LABwEQ52Q0bNmDr1q0AJrwE55xzjnz8BRdcgL/97W+48sor0dvbi1/+8pf47W9/i0suuUQ+5tJLL8V//dd/4eGHH8Z7772H7373u9i9ezdWrFgBYMJi9a1vfQuXXnopXn75Zaxfvx4rVqzA4sWLOXG7jBRqwSQLv7pSSzVb88UNVaybTt2ojV5zLYb5iBU9SkmpwwxKCZU1NcIhhxwCYM9YGBgYMOTSN4Jo/RUh1zwJFuqSl3rPno4RE8wrDW3SaqGGKt6on19zczM8Hg9GRkYMP1+ja5FRQTiXJ70YRMu1eK3hcFj2Xni9XoXgZGSckueZvENm1y318WqFhKr4VAKbzSa/W/UzU+eRAHvGgph8K2L1upVOpw17hvIhjg/1GtXX11ew5VzPs0TPVswREZ+X2+2Wy7abRd31W89jI3oqCDEUq1xQKWwtxa0UiozVmFYsHnnkEbz33ntwu92YP38+Pv/5z2Pq1KkYHR3FI488Yupcf/nLX7Bw4UIsXLgQwIRSsHDhQlx77bUAgG3btslKBgDMnTsXzz33HF588UUcdthhuOOOO3D//ffLVQsA4IwzzsDtt9+Oa6+9Focffjg2bNiAF154QZHQfdddd2H58uU4/fTTceyxx6KtrQ1PP/202UdRM1RjbL4ZYUoknU6jp6dH83zV2lguEAho3itZ4Kox/twqSEgQy+WVglotfwoY78qu1eyKBGQrGisaabzU0tKiaBwH6D/7cm/GRhDHoRgWlUwmFZ4VCv2gMBdAO95ZCyNrkRkFMFdJ3GIQBXh1mBUZftQJ3EYEYRJ8RkZGLNt3tMrrVmLOi/NM/czEPjS0ptMx6rK1pQzbsjL/hryfVFxC7OtSaKM2SZI03524XlDiuRiaOTQ0VLDcIJJL0SVFMB6Py2WdrULrnTudTt2xEIvFNPMpqnFdVWOTTF6l3W7HlClTsGrVKjmUCJjIUcjVXKfeGB4elqscTJ8+vdKXk5Ouri653GcxYRJWQi5isZSsGcRk3WpFdIsPDg7quo2r6b2UGnV1FCuguvCi+79aQ+O0iEajOcNiCJfLpTCiABMdjmkTVv/MLHQuEa/Xi/HxcaRSKTl8IBgMZiUy5sKqnAAroLVQD7/fj3g8Lt8vUew9qMejnjUyF1RyttRjOhqNytX7rKCY9U0MN6LzqEOQzGK324tef8QOzYQ4ttT3rDfu9K6FBG+z10kVxKwcH+prp3vTWi+sRl02mJ57IZWpjMxhrX3D6D5uFWLFSK1wM6/Xi2OOOaak16CFGZnXdI4FANxwww04++yzEYlEcP311xdyCqaMVGNsvtmysyIUttDb21vVwiPdY76Y0Gp6L6VCXJytVCpsNhtSqZS8+Pb29ioW/1roZ5EvkTcXoVDIsrkdCoUUSp/H45E3MK0NVxTu9HpgaHlZKonb7VasOep68zRuxAR5KwS1XPkMWogeD3rOmUymLD1aKNzEql4QVnnRqKlcMdfk8XgsyfWKxWJobm6W85zi8Th8Pp/srVHfs3rcEZlMRpF3QfH8+dZIUSER35PRvhVmIPlBNN4AE/PCSgVUi2QyaVleitF+Hnr9LOjfpVQsHA6Hoq1BLBbLKmu7a9eukn2/VRSkWJx11lk4+uij8dWvfhUbN27Eo48+avV1MRZSS/0PctWGJqZOnSpbUau5GRot/LkW3mqy5paSUuTA0FgRn6+oZADVr7R1d3cbbrikVfpRbGwl/r8QxBhxYKISDtXr11pDRMFbKyyisbER8XhcUfO/0qi9o+l0GrNnz5avv6WlRSE8WTVeKXRkdHTUkLVdbSzREuysQq00+nw+JBIJZDIZ2YNj1lpLiiYplsV6EKlEaD4hNleTR8C6MBKqKiYmlVPzT/U9+3y+nF558ZqMGl3oOFoX1AqOlU1kc8kPpfZYAIXlpWjJEcXOG/K2lwIqlmAk4qcWooIKSt4GgKOOOgpvvvkm+vr6cPTRR2PLli1WXxtTx4RCIUXZVcJI5S1apCnmkxoWVRv5NrFqs+aWEjHuXCyHWAxaz1eSJEXM6uDgIIDqzDMCjJcxdbvdAKB5D1YqbWK+D1WEUkOhKMlkUrNOPkElOqupoIJ67DkcDmzatEkWkMSEVKqnr1UK0yyUR2NUsBWfGcWEU8lfq9cM9fih/Ip0Oo1t27bJx9HnHo8nZy19u90uC/exWEwWsAsZC7RPABNjM1dugt/vz5u7YFVvi9HRUYyNjcFut8t9UEjgFb2lYjJyKSAhUxwf1KOh1PNOXBvUOSQOh6Pgfgv0PIvp13DIIYfIORKUmzM4OFhUT5pi+2fkorW1NWezPRH1s65GTCsW4k3Onj0bb7zxBubMmYMvf/nLll4YU9/QRqueqEYX/kQiIXcHpXCjWoKaGVWrwGs1YsLc/PnzS/pd4hql3uSrScgF9mwSuZK3tTqMi1hZuIDKPOZC/P5UKiVv1C6XS3PTI6WoGpg/f76ccOvxeJBOp2XhLJlMysnaPT09CmssjaNi5quR9+NyucpehEI9ftRV7GjMkdLpdrsxMDCgK8Sr13TyVKj7Pxh5joFAAOl0GslkErFYLKdiRs3FygX1QCFLMz0PmgN6idtWU8r1QKuHiBqXy4Xx8XHFZ9TsrRAymUyW59kMZLALBALo6OiQK0gW05MmGo0ilUqVLOnezDWpn3U1YlqxuO666zB16lT5/01NTfj973+PSy65BMcee6ylF8fULzRRiUImrLiBVVsnbr3SnQQ1BbKyeVytIFbJsmKhznUO9SZfbaFR1MskVwlBt9utqNyjvodSliHWCr8SvRok8JGCrxX/W0hxhlKhroSkh3ruUh8L8mIUMl/FBoR6FY0ymQzcbjcikQief/75sqxp4jMh7wLNG2pyFwwGFZWi1KVAc+Hz+bLGqBFFnypzmbESl6sUJ1nVxSRreh4k+NE9l1IQ1PLYW7keiF4X9fsSu3Nbva4WE7JGXjKRaDQqvyuay2bmFuWSOJ1OLF++XLORqJUVpOicWp9V2x6mhekci+uuu07z8xtuuKHoi2EmD2KCIGA+VtPv92NgYEARb1gtuRbqqiVGqpDUwmJRDOoYa3VccDEJlbk2IXGTr4axIRKNRg2Ne1KaS12oIBqNKuaTXqherh404u9TnHMxIQ2lhO5NXaWMwiZisRgcDgdaW1tli6Ldbi9IQRWroVFlFS3S6bT8M+qkXM5xK+bYLF++XP6cDEEul0tuYGZ0zaZnR1WUxHPlanKq7lngcDgsjy8vpslbR0eHXN2Lus87HA45x6HQiodGofyXUuYZ6iVuA9praqmTudXo7a3q50FhTB6PB+Pj40U/s+bmZvT398vfrc71s4Lp06djfHwcPp9PHqOSJGFwcLDq9jI1hhSLZ599FsuWLYPL5cKzzz6re5zNZsPJJ59s2cUx9YtYqUqr6U0+Yby5uRkDAwOKjY7CMnJ1HS0HolXHyMY1GRK4RQuluCEVUzbSCNWssBm1eqdSqaznV6rrEeecngKhrmGvt6E6nU7LS19agZhUOzAwAIfDAbfbjdHRUTQ2NsoKAJW27OzslH9XT4nNp/iJzzaXsKmuxFTu8atXQZDWaOp+bCRkTiQWi6G9vT3L6EJV3ERFhcZWKpWC3W6XE5QDgQDWrl1blLBOVmB6jzt27JA/N2Mlz2Qy6Orqkj0kkiShsbFRUYrZiupTetdGXaPNdjU3i55BRmu8aymDpUYtI9CzUveG0urmbuaZiYYweu+lbspK41ysEgXsmUvVjKE+Fna7Hf39/dhnn31yNqWx2Ww1kbFuBbXUx6LaEZP7jGCz2eSu1VRTW6/WdjkRBRayZuSzslnRf6AWUG9ExdaiN0K1P1szvQzC4XBZPBbiOzFT952E9FQqlSUEVZvibKZkqV7JVXVPCb1eQWb79SxfvlxRo1/0GpQL6mMB7Cm1K87fQq3SooU9Fy6XC06nUz6Ofk+03BaD0TK6ZpUNs8cXgtFeJqX2cKrHu9bel49SeKCIUu3/dN/5DJ+5vJJGcTqdmqF9ddPHQnyApdbSmMmL0ckouh0lSUJXV5dsoShVSUYjiMoRzZN0Op1zEdKKYa9H1JYv0VpPVklA3/pN9d3NbNzJZBKrV6+Wm7xVW7+TXCFFWpRaUQ4EAnIIkNGKZWLMPFn11IJjtYQoAhPKXC6hUr0G6R2r7imhZ+k306+H8hrIy1mpEDLR6kzzVBRSCy2UQfeUSCTksrAUakZhreSdGBwclJ8ZKRNWeAAaGxtx4IEHGjJqmFUSytERecGCBQAmwvcikYhmoz4g20NsNeJ4Fw0SZiz5pXxe5FlQN04tFrrvfPc5MjJStKKply9UC8nbBfWxYBirUC9KZqEqUqlUquLWaXGxFTsp6y1AXq+3agSuciMu0MlkUg49IQGDnhkpZcUYNEhQrCYBF0DeDtaikFuua29vbzftZichpre3V2GBFK2r1UI+i7dZKyMlzubL4TGiGA8NDSEajRb0DqwkGAzK4WA+ny9LSBX7lpBhxIig7nA45BAqWrfdbrfiXvW6frtcLjQ1NRVtBbbZbKbCMIsRDq1qMKhGvG69sJhSN8UVx3tXV5f8uXi/osVd61mUykjt9XqzlFGr1k86Rz7lOtd7p31O6xgjxtVqWk/1MKRY/PznPzd8wosuuqjgi2EmH6Ll2uhCo7XYJ5PJiudXiIttT09P3g2pmqrllBt6VmI4UDwel0PcrFAoCNrg1HG3lUIMG8jlzRLHRzVvJnrhMdOmTauIyz4XVoQo0DlsNptc6lpvvRGFcCNCZk9PT8WVX7FRoug5oTGoVqKi0aghAZy8tx6PR9ezrBejn0wmTb83rVAbsXlcOBzOUi5orRBzMQpRLhwOhyWKkPqc6rwsPa9WOYtVUENFNS6XCwcffLCpUMBioBBpqrQHQNHg0SqK6Uxvs9my8rZEpkyZkvM51UrvK0OKxV133WXoZDabjRULxhSi1TadTsvJ2Nu2bZMXc/XCrrfIiz0LKj35qJNmvmMmO/39/fK/SfBPJBKW1gsnq9nAwIClHWkLRd2xNx+VyFEwG6OttcFWo+IsCh16qAVS9fpD92VU2DQjfGQyGaxevVo3xKWUiO9cXJfj8bhuGJ7ZXCnycOhZfIvJ4RDREvrsdruiipLa0CAqGkb2m1zY7faixz8pDVSdbP78+XKIWCXGhx7qkE6yyNO6m88zaxUUIi32gil1GJi6opyRa6TwXK3xkc+jKu6X1YwhxWLz5s2lvg5mEqPeBOLxeEELeymsE2YRyw/mo7W1tQxXVN2Ii7K4qBaT1OdwODB16lQMDQ3JMdzxeBypVKrklZWMQBvT6OhozrHtcrkqVlXJTIy2ugoaJW9WY5lZI+GW4tijXAAt1KFeZgsU5LKEV6Lyi/jOSZHIF06jzpXKJWS5XC45h0f9++JzUz87sT+AUS+A1nW0tbWhublZ/i7x3dhsNkuF32QymXPs5MPr9crvXxwH6u7olYZKB1OeXEtLi3xt5M0rV5Uoo4ntxSKOUYfDId+3mcIChSqdtZLjzDkWTEVRu77T6bQcPmBmQQqHwwAqv+CK9abzUelrrTTUtEhvsaSfmQ1fEa2F6XRajusWyxJTXHwloI3p+eefz6lAtbS0VOwazcRoi8dW2lOYD7PW51zFF6ZNm4aRkRHZ+k5rVk9Pj5w/RGidQ30toqJRCaXM6DtXezbEXKlcUDNFwuVywefzyQK9KNj39fXJSqp4XqPrgNa8isViaG5ulpUmsVeJ3W63XAAeHR0tOPSO8m3U86nUuRNmEffvBQsWyMqpy+UqWY6JHtQ8rhxrkKj4UglmMcqiVJS6i7tVGCo3q+bvf/87nn32WWzdujUrQ/3OO++07OKqGS43aw20SY2NjSGTySisDmbcjFSOU6sEXjkXG7pmCoXKNb2qrQxnuVGXCFZDymJfXx/Gx8eLLk0olrGsRDliQJlfYcTCVYlyo8WgVaq0mjBS4lf0eKk/F5vlGcHlcimKERihUmNTjV4JXXHe0hpmpnQyMCEgjYyMZAmeHo8HAOTvtaoXhHh+rWdbyvLXhTbho/WqmhV2Pc9SKBQyHSZkBeUqMS7OAao6WOhYzefpE8dPJUuom5F59ZtS6NDV1YV58+bh3nvvxR133IGXX34ZDz30EB588EFs2LCh0GtmJimBQABLlizBggUL4PF44HA4ZHe8mRyESCQCn8+n6IoruvbLATXBcrlceZUKZsL6Ru8cmFg0w+EwPB6PLLD09vYikUhAkiTZCparl47R76yUxU9Mis1HLebgkAWTwiCqjebm5rzH6MXHp9NpxONxefzZbDaFMKWH2XyharBG63XHBpTXR2vr9u3bTZ1/aGhIoVT4/X7Y7XZ5rgMTyoXTaS6oQv0+Ghsb4fF44Pf74XK5ZG9lNBpFZ2cn1qxZIyfSkiHDagrtvUEljSORCKLRqPx5NBpFV1eX4rNKEQgE5DUbmPBKOZ1OBAKBkioVVubgFYI4B4aGhhAMBuV1waxXgbzyelgVIlxOTO/QK1euxOWXX45IJILGxkY89dRT+Oijj/AP//AP+Od//udSXCMzCaBa8ADkzSwUCsnCpJHQgIGBASxZskS27pRbiHS73QCg2SRMi3IpPNUKKZWUa0KhP2RR7Orqkjcnh8OBjo6Ogvp+2Gw2eDwe2YJO4QSV2JjzjUUSgjweD+bPn1+mq7IO8f6qMR4435zLFxKTSCTk+3I6nejo6EA4HFasU6Jwm0wmTQkD1eLFJAWRhEQRUZik912MwEPPi56rmJeQLydGLWCqj7fZbFiyZAna29vlinN9fX3y/WUyGfmzQCAg7zNer9eQ0mgVjY2NWZ+Je4g4bsnYUi2KO63ZoVBIMSbIMFIKJUC9v9rtdjknrRyIY4X+pvFLFePMGMCMhsvVirHJtGLx3nvv4ZxzzgEAOaxg6tSpuPHGG/HTn/7U8gtkJgfkUhU3s0AgIG8GAwMDeSdqMplUCIu04JWyU7FoOTJTLaaSVvNKomVto1wTMeeELPv0LJuamtDV1YXe3l7TAuvMmTPlcUDjLJFIoKenp6yWPwqDyrU5+Hw+tLe3Z43barJS5iIQCMjztBo9dsFgMKfXy8zYamlpkWv4d3R0oKOjQxZgC2VwcLAq3nM+o4zaAFBM7Lc698IM6jFG/xcrAxHiPYnjQPTKtLe3Y/ny5RgfHzecH6ClFJhlbGws6zNRIK+FvUK9386fPx8ej6eo+WAUaqxYTqW8vb0d4XAYAwMDWWF0kiSVTAmo9NpgBNNvfMqUKXJexcyZM/Hhhx/K3SA/+eQTa6+OmTSIFhmx1rjZ8oO5aspbiRhTTNeer965mIxcbfX9y4W62pBeyIW6wykpbSQMmBEARYVFtPKpuyeXGjGxUc/CSwmmdD2kjFBybKUrWhmB8ouq1bpmNudBD6p+RUKFejw7HA7TVYGqpWS2kT4I6ryEQnMJ1BTbsRiYUPqoAhSw537Ee8p1f7T+GMmFKrTyk4jW/drtdrjd7qwcC+pwX43KhpjXSFClKOpdUoqeFvn6yZSKXLJJMpnMO5bNJrjXyh5g2mNx1FFH4fXXXwcAnHjiibjsssvw4x//GN/85jdx1FFHWX6BzOSArEnhcFhRa5y8FtWGuNmkUin09vbm3QxJmFEXPJhMqC2heiEXontdtCy2tLSYEgptNhsSiQSef/55RKPRrM7QWjHkpYLuPRQK5YznFpVsUkaA2vFyUUhEucISzEB5UFYgnofeGQkayWTSkMCptnZTGFwtvGdxnAaDQdN5FoTdbldY560InRkYGCgqx45CJuPxeMHemELvw263yx41rcTtUnvii4GeeSQSkUO2qBJff38/xsbGqrK/TalwOp1Z3lFxvJutmlUra4Npie3OO+/Erl27AAA33HADdu3ahSeffBIHHnjgpKkIxViP2pokWj30OnuKUFMeAGUpJUrWOeq+a4Zq6QBdCei90Gafr3widQFOJpNoaGgwVaJXtBal02lEIhGFFT2TyaChoaFsG7RY/7yvry/LmkXVRbQ8N9VcGYYQrZXVUNVICyNriVHE96c1fo1Y3kXlw+PxVE3TMyOox2ahVZUaGhoU78QKxY/KluutLUYqBpptZKlG791Td29gT3WsdDqt6ANRSx5KEXV/EI/HY6gUcaFQlSIS1Mux94tjJxQK5Rz3et5RrbFBhjPyhKqppqaI+TCtWOy///7yv6dMmYL77rvP0gtiGGJwcFAuRZuPeDyOTCaDdDpdFpdoe3s72tvb0dnZKS+alOyXL3l7svevUDfiyveu1MqH0dA4rXdAC7ZaiC9naWLRCyEyNDSUlbxrJCSlWjDTVK9SWDn3qFKZGNudT9DIRbVZIvPNCfXYLLRnQyKRKDr8yWazwWazKYo96M0dMYQr11g12shSD70GeaRUeDweOZeDigZs374d6XRaM0ekFhANRzRu9EoRiwqWGqrgpffcbTYbDjnkENlARV23xbDEUqHVSFIdEigqh+I1S5KkqTTQ3NELu6uWog5GKbxuI4Bdu3ZheHhY8YdhioVcqLFYTFGFJReJRKIipdjEylWhUAhOpzPnJkRNfCYzZqt1ia7/QCAgV+LxeDymww3o3QwNDWFsbAy9vb2yAFWu0sS57ltdWrKWqHQpXyNQ5TaRfGPIyBijhOtAICAbGMwIo9UoOJidE8WEeEqSVFQIVGNjIxoaGhSf6c0jdQiXXmEECocqVOHJFwo3Pj4ul0inPS6dTsthhGJyfC2tCerEfookAJRzKZVKyeu4Oh8rmUzmfO6SJMk5MFRJjyj1Gq5e58QqaeFwWC6+IcoG+UoZqxVydfjd4OCgtTdRYkwrFps3b8ZJJ52EKVOmwOv1wufzwefzYcaMGZM6xIOxHrvdrii/aZSWlpYSXpUSEnRbWlrkakO5qMZKOeVGHSNstuKRaEktJv9GLDVZTUJxrSoX1Rz7TWhZ1CVJyhlHn2vOUtIoVRlbvXp13hKpWlTjMzM7J6jSUqGYXRspF4HypOh6KSFWT8AU8/ko1FJPgbJSSFXnklBflCVLlihKl4pzqNy9mKxEyzMrrtc2my2r5LgZxPyNgYGBsuUnaa1zep+RsREwN75nz56t+L8VRRHKield+ayzzoIkSXjwwQfR2tpa8UYlTP2hjo2lfw8ODhqaYGK4QylCXMi1S/GxgLkkrGoOFakEZkNoKBSqt7dXHiuSJCkshEYqR9ntdjgcDnlslOudGKk/z2OkNOhVLiomoZTmPo03swKyFeVKS4HZOSHmEJWqi7UI5STEYjFs2rQJ8+fPx5IlSzQrE2ldJ5Erz0vMGaD1otB8ARofdB46P7AntFbr+6u1AlQ+1Nfe19enqAglSRK6u7sxMDCg+0zzrePi75GSVk2IyewUGkdygxjqpQ4jVM8fI328qgnTisU777yD9evXY968eaW4HoZRLPxdXV3yxAyHw7qJTSLqzrBWx32TYCIuBEZLoFaLVbyaKHTzTKfTuiUh872LSiXCRaPRvIJJOStVTTba29tzCjJGKFbAVFNvxrlAIGDYCFQIYj4GfQetBVYr44FAQDZkOBwOua9SMVDPBSNeiFrKsVKjVd539erVimPUY0Sdl2I0kb8a10wqPU2Mjo5i+fLlsuIL7FGMxsfHc+Yo0XOq2+Ttz33uc/joo49YsWDKgmgxMlLSlY4jtHokGPVikGdCLYSS1VO0PGQymbwJiNVqmaw0ZjdP2pRTqVTBVVtisRhisRgcDgdaW1sRj8dLnrhtxpJbq8JELVBMgjUwIexYUblIqwpYvVDKAhW0xra0tGD37t2yMObz+XImZuut+/mMT2oPutmxIxqdbDabfL5c31nOYhLlRO0xtNvtkCRJDkc0kqdD+6zNZoPT6Sx7Yzw96J1Rvw5K3qb3T14HsZ8RhQ4a6e2h7nFUzZhWLO6//35ccMEF+Pjjj3HIIYdkxVQeeuihll0cw4gWIyA7zEBLmE8mk+jp6UFDQwOSySRcLhd6e3sVoTNGvBi0AMZiMYViQW5rWkj0rk0NWWLKUbmiHlGXMxUX8kKt0Ol0Wn7Ppa4mVotx0kzpGBkZQUdHR6UvoySoy46Wgm3btinW3IGBAYXQqlbY9NZ9IyWv6fd9Pl/eikVqRCW0sbExq3KSFurcilpXMsS1u7m5WTb+iUo+VcXr7e2VC7GolY7x8XHFzyh/gRrqVvL5aJUnVleNA7RLiKs9Obm+oxbGgGnFYmBgAB9++CFWrFghfyZqkJWozMPUL5TEShNU7WLXW9wzmYyiezMJnTShjYTekHVFHd9IiyQtckZLJYrWmlpZIKoJMbdC3Ul3zZo1RZ/fqsZpegSDQfT09OT9nmpsLFdPsIJXOmhtLEeRCvE7XC6XYl5pVdkqJl+h2J4WwEQIlM/nk4VgvXwAsXs7XXM1l3E20xNEDAHyeDwIBALymmiz2dDX1yevf6JiOjIyglAohJ6eHgATMmdjY2NVPR+tbu3U+FXtzVA/K5INHA5HThm6VrybphWLb37zm1i4cCF+85vfcPI2U3LII0CLEGCuQgJZtslCYiZRVy+hTl3tglyyuaxYtNE1NzfXbDJetRKNRgtWCsQwBXWHVCuhjUXrOo3m5zDFo457Bow1sisF9ahAGu0xYzUdHR3o7OyUe0AY6btBqAVTLUFZFBrj8bhmXlcuKLfOSP8MrT2vmvcMI4K92+1GIpGAJElIp9Ow2+1IpVLo7u6WE9kzmQwSiQQ2btyYFTZMVb4cDgcymQycTqdCMauG5yOOL/U+r1ZM1c9KjILQMj45HA7Mnz+/KhVLLUwrFtFoFM8++2zFXyIzOdCqLGEUMeTIzITU2ljEz9SbTL7Ooqx8W4M61pneSSHlPQnaqOj8pUKvKR6Q7SmptOWtnhFDFwmjSoXYn6KYMccUh7qxGr2XXJ22c6G1xxixgJsxCIyPjyty/3KV5ldfT7UncBvxBFH+gDoxW0s5kyQJPp8va72kkFeKXiCq8fnoVR4TFdOurq4sDwb9ef755xWeC7fbXXX3mAvTisXxxx+Pd955hxULpiyoJ6jWgiNCybhiyBGFTxmtBKQOuaHzaHWLJgsDoG/5pGY+6prptbRQVAPiu6C/aSzQZmM2rpuETHWDJqvJF3PucrnQ0tIibzJMadB6D1SJRqxIozWXnU6nLNxYQalzeiqBVqlwI4mpZkilUnIFHUqGBgoXMI2Un1VbnM2GRKXTaYWgGI/HdUOIqlFQzoWR6xULnlBjQBKw1R4gp9OZpXB4vV75M9GTUyuonxFVu6RxJK5JYo4JUWt7gmnF4uSTT8Yll1wil/9UJ2+fcsopll0cw6jJV22EknH9fn/WgqVOws5HOp2W42FJoZEkCatXr5bzLozkewB7LFS1XJe8GhAVM1FIpFjWrVu3YmhoyHR4UTqdRm9vb8mSJOl8WsqF1+vFMcccY+n3MdrQe9i4caM8X0mZEK2pWmVFi42xV5NMJuVu3fWCVqlwK58ZMOEpKOV80RKUteLnC4HkpVQqJRuwJoORSS+sWCQejyORSGR5o9SleWt179SLehCLh1QijLAU2CSTwaW54pAnU/L28PCwbDWZPn16pS9n0mCmZCct4jRZjXosKCeDfo+6fhezQXo8nqpr3lOLqK184v+tqp1fyneVa/xWqrfGZMRoFZZCMKPU1vO6QHPTao+FVqUd9XeWqkJQoc3/tJKNXS4XnE5nxasZVQtU3p0Q33M9lN8lRVs958WQXpI5yIMtemkqvU6YkXlNeyw4yZCpFajULCXOGV2QKGSJJjlZSIxU9NEjV0wtYxytpku5vAFGUCughfQ8MYOe4GnWo8YUDoVmqBtyFYuZ84n9depBcFJT7LzUI5eVv9Shpvly/CjpWD0GbDabZrJxvbzrQlCPeXU0gvie1et+Lc4XvWgFujete6rVYi+myqAkk0k4nU5s3LixVNfDMDkxWyqyqalJXtC7uroQjUYRjUbR2dmJzs5OuZytmmAwCI/HI1dzCgQCppQKStimv0vZMGqyE41G0dXVVXCSPDVZEt83oS6TWAzd3d2IRCJZ44iuu9R5Hswempub4fF45IZcdrvdkqpg+ZQKm80mj7OOjg7NHK56ohT3k6vLMq3bpRLE9M7r9/vh8Xiw1157aY6BSlQdq3bUY179bKlfiNYeXYvzJRAIKPIzc/2cZJTe3t6aUp4IUx4Ll8uF2bNnT5pwJ6b6MBvrOjQ0hK6uLtnNKNbRBiDH1WtVZxgcHEQkEsHWrVuxc+dOU9cpbiSl3OiY7CTulpYWDAwMGG5iRRU3tBZvyq2xwuOUb7y2trYW/R2MMdRVuhwOh2Z8M7n+rcJutyMYDCoadgYCgbrKvVLHkoueXqtK++YSzgoRwoxawLVydCifL1c+iXjPZrwqtWiZN4pW9SvKkaMk70QiodlQtp7mixZixEQt5uCYNtFcffXV+NGPfoRPP/3Ukgu45557MGfOHDQ2NmLRokV46623dI9NJpO48cYbccABB6CxsRGHHXYYXnjhBcUxO3fuxMUXX4xAIACPx4Ojjz4af/7znxXHnHfeebDZbIo/X/nKVyy5H6Z0iIusGQ8AJYSRgB8MBuFyueQQmEQigVgsJi/2ZAEnQXBoaKjgEChJknJaKZjiETcXSt5PJpOGBRgtwZHGAFUAssLjRAn/6hrtdJ3s1SofZNkmS3MoFJLfj9frlb0KxxxzjGaRkkKhCnHJZFJh6Mhnzawl1IKz6ImrVsu9GQt4IBDASSedJOfeGankJnrDzHhVatEybxStMU8exJGREUUSN90/rcsAanq+0H1QBAX9myAZhfp96EVWVCumcyx+8YtfoK+vD36/H4FAAFOmTFH8vLu72/C5nnzySVx66aW47777sGjRIvzsZz9DR0cH3n//feyzzz5Zx19zzTV47LHH8F//9V8IhULo7OzEV7/6VbzxxhtYuHAhAODb3/42Nm7ciEcffRR+vx+PPfYYli5dik2bNmHfffeVz/WVr3wFDz30kPz/hoYGs4+CKTPqikDqhKdckICvPp/a8yFW66BYeNHKVojFrZ6tTpVEfK7hcDhnr4h8dHd3o7m5WS71R4qkmKNTLO3t7WhubtaMOc8V3sFYj1bjsXg8LofCUQEH8ip0dHRgzZo1RecYOhwO2WMB1G6Fm1yIXuXOzk7F+pxr/TS6tpai30whFnCx6VtfX5/s3RK9XHa7HQ0NDYrzmvGq1LtlXo16XxcT3MXE+Vq04ouoFUa1B0vMT8pkMjV3v6arQt1www05f37dddcZPteiRYvwuc99Dr/4xS8ATCSGz5o1C9///vdx1VVXZR3v9/tx9dVX43vf+5782emnnw6Px4PHHnsMiUQC06ZNw//8z//gpJNOko854ogjsGzZMtx0000AJjwWn332GZ555hnD16qGq0KVHxIkxbAlwFhyoM1mwyGHHJJVCtHlcmUpJqRQkDuWqpqQ5c1sKKDH49GsBsEUh1aVDXVlEQCa71gLraRqdc6FGbQUyueeey5LeHK5XOjo6CjoO5jCofGjhhK7CXXX5GIoZjzVGnrPt1iWL19u+TnNoK4aqEa93kymd24VWmunOJ5q/ZmK9wdoJ/OLSnk13G9Jq0KZURxyMT4+jvXr12PlypXyZ3a7HUuXLsW6des0f2dsbCwrjMDj8eD1118HMGFtTqfTOY8hXnnlFeyzzz7w+Xw4/vjjcdNNN6G5udmKW2NKBFl6aIEhq59elZ1wOCwLA5IkKWI1KXZe/D2KlU2lUshkMrLViRazQnKLxBKDk8XqVC60nmt7e7sc70zoxc+rEceC3W4vOplaK5Zay47T0tJS1PcwhaGnLKhLXtI4KwaHw4H58+cDgNwbp9KCQqkhz4zoASwWdYJ9JbzB6u7t6qZv6q7s9dgIsdTk6iVSD3NHq7ohoVV2t9bu17RiQaxfvx7vvfceAGDBggVyKJJRPvnkE6TT6ayExdbWVllgVNPR0YE777wTxx57LA444AB0dXXh6aeflgW+adOmYfHixfj3f/93HHzwwWhtbcVvfvMbrFu3TiF8fOUrX8Fpp52GuXPn4sMPP8SPfvQjLFu2DOvWrdMVJsbGxjA2Nib/f3h42NT9MtZhNgxKhIQ8ime32+1wu92KxarQWuVaHHjggXWzGNYK6u7sZsqJig2ZRKWgEAFGrfjoxclybkVloPdI+41YNx6YUAYo7CaVShWVeEwNGIHc5VLrCVF40vIkOhwOZDIZQ89UK6QIKH15WS3EUDatfhrRaFSRsM7Fbqyh1jqSF4p6npQi9K/UmFYsduzYgTPPPBOvvPIKZsyYAQD47LPPcNxxx+GJJ54oqfXt7rvvxvnnn49QKASbzYYDDjgAK1aswIMPPigf8+ijj+Kb3/wm9t13XzgcDrS3t+NrX/sa1q9fLx9z5plnyv8Oh8M49NBDccABB+CVV17RDVW55ZZb8oaBMeVBrPucq7dET09P1me0MYkeC7WwODg4aMl1+v3+imx8kwW9Z1uMoE4dvOn86jhxM++Rjuvp6cmpqLInq7KITcp2794teyrT6XRWFTk9HA5HXgFS3VdnMqE1J80I3G1tbZrnqIQ3OF9PBfpjRU7OZIHzECeIRqNyBAZ5OWvxeZiuCvX9738fO3fuRE9PDz799FN8+umn2LhxI4aHh3HRRRcZPs/ee+8Nh8OB7du3Kz7fvn072traNH+npaUFzzzzDHbv3i3HOU6dOhX777+/fMwBBxyAV199Fbt27cJHH32Et956C8lkUnGMmv333x977713Tnf3ypUrMTQ0JP/56KOPDN8rUxoCgUDOpHv1ok4FB4A9G50oPBBWdG+m87jd7kkpSJQDveoq9LnX6zV9TrfbragfHo/HFUKlz+fLquCRi76+vpzCRWNjY01uHPWCOolSXSGMqshRBSmq1KLGiJBM3rBarmZTCNFoFKlUqqg+IQMDA5rVkSpdTYt604jXpu6rw/1p9qBVAYkiBOq1+pUZaL/weDxYtmxZza4Tpj0WL7zwAl566SUcfPDB8mfz58/HPffcgxNOOMHwedxuN4444gh0dXXh1FNPBTAhCHZ1deHCCy/M+buNjY3Yd999kUwm8dRTT+Ff/uVfso6ZMmUKpkyZgng8js7OTtx666265/v73/+OwcFBzJw5U/eYhoYGrhxVRYiJ3KlUKmfPAgphIGWiu7sbiURCLjWcSqXQ3d0tlxa1qtY6MCGoVDrZsF7Rc41rhWA4nc6s2GctRMEyEonA7/fLv6cVIpWPYDCoqHmvxmazsbWugqgt3pS4LeZEiO+G3pXZpGS/3z9pu6pTTgIZAcSQMBGHw4HW1lZZmRfnKyluVvSTsYpoNKowQo2PjyvGh8vlksNsmQm01k9RmZjsz6pe8jFNKxaZTEazprfL5TLt9rv00ktx7rnn4sgjj8TnP/95/OxnP8Pu3buxYsUKAMA555yDfffdF7fccgsA4M0338THH3+Mww8/HB9//DGuv/56ZDIZXHnllfI5Ozs7IUkS5s2bh76+PlxxxRUIhULyOXft2oUbbrgBp59+Otra2vDhhx/iyiuvRDAY5MosNYS4uXd0dOSsQEKKQiKRUMT6SpIESZKQyWRMeymMKh9UG5+pDKQsGlEq9H5fvS7oLfx6CkKuceLz+ThcroKoldP29naFAkDrCgk/heZeTeY8GnXSLY139RrqdrsBQDYUAZDLjVJBjWp6jmrrOuXRhEIh9PT0yL1KBgcHeV7/f7QEZ62k7MlqbKmXPBLTvsnjjz8eP/jBDxSC2Mcff4xLLrnEdCnNM844A7fffjuuvfZaHH744diwYQNeeOEFOaF769at2LZtm3z86OgorrnmGsyfPx9f/epXse++++L111+Xcz2ACYvj9773PYRCIZxzzjn44he/iM7OTlkZcjgcePfdd3HKKafgoIMOwre+9S0cccQRWLt2LXskagBypfp8PkUYTC4NX3RFx2IxWdgnVzUA02Ez+dz6Ho8Hy5cvn7RWymqhUIWCSCaT6OzslF33uUIvtJpZ6RWiIPr7+001zGKsRys8gyALuc/nU7zLXAYDcS1pbGyc9O9WPWdovDudSrsmNSoVPRmUs9nS0lJ1z1HvWgKBgMLIalVobT2gtX5qfVbPjQEnA6b7WHz00Uc45ZRT0NPTg1mzZsmfHXLIIXj22Wex3377leRCqw3uY1EZ1L0LRMuGViK31+vF7Nmz5XAUsn45HA7ZW2Fl6BNRDXWnJzNaISvURG98fLygXiT5DCeU9wXsqRajbhCmBYfKlR9x3aBxovWOxfVGbNrV0dGhWT2OSo+K447KWE8266sW6vr9uTxA5K2o5h5AlKBts9nk0uKBQABr165VNMk78cQTK3yltcVk9VhUMyXtYzFr1ix0d3fjpZdekjfRgw8+GEuXLi3sahnGBGpXKgkFkUgEjY2NWaVFR0ZGEAgE0NPTA0mSZAu2KFgW2wlWjd1u58WwwtC4ECtsDA4OyiEVRvD7/XIolZ51Ur0BUuOs3t5eDA4O5lUqOFSuMogWUa2ywPR/rfWGKoWl0+msHjpDQ0NZ3kyyWHM/A+V6HQ6H5ZwWsVs1MLH2hkIhDA4OIpFIVFVuhQiVzHU6nQrFZ3x8XP73ggULKnFpNU29hAQVQj0oVaY9FswE7LGoDtRWQ61meepNyyhGSkhqQRZNyueYzImblUJ89uqGeUaw2WxwOp26deq1rN0+n08WIo10+xZr89fqBlKr6G3e4nqiZyU30lFazEOkBnGTtcO6npeCQpvoZ4ODg4r5Q95lSvyuRo9FrnGUq9cFw+ihjsqoFkrqsQAmbryrqws7duzIEuLEnhIMU2oCgQDee+892ROhZY02qlRQzgXp2maVCpfLJf+OWDEkFouxYlFmKMmzv79f9liYeZ+SJCGZTCISiaCnpwcOhwMtLS1yZ11qciYKRqKCa+S7HA4HJ21XCD2LqJEKNWKlL5vNhunTp2etMU1NTRgfH1d4OqopP6CcqL1DpHT5fD55zpAHo7m5WQ5XJI9jteVWiOSqTMdFGRg9cnkl6qEylOnk7RtuuAEnnHACurq68MknnyAejyv+MEy50UvQtdlscLlcaGxszPn7LpcLLpdLrhKlPodRkskkMpmMLHRSmAuHu5QfShAVu9+aeZci9E5jsRgSiQTS6bQs7IiJh+J7zhduRYpKNQtNkxEaN3o5UiQQ0DohSRKGhobg9/vh8Xjk44aGhhRC5WTrXSEiFieg0rNOp1MOMyQoVCwYDMqKucPhqNlnR0pUKpUy3PeGmRzkSk6vh/XCtMfivvvuw6pVq3D22WeX4noYxjQUp6uGrM65hDxy7TU2NmqGrpiNFFSHt7CnojKQJVFMoqScGUqwdbvdhrxZFBaVyWSQTqdlYYeqCZFioBaUcpFOpxGPx6vK1c3sGTfiuxU3eMqhUedfxWIxRW6FuA5MdtRWfbLG6lVME4WtUChU8usrFey1YPSoB69ELkx7LMbHx3H00UeX4loYpiDa29sRDofh8Xg0u5zqhaXYbDZZsFQnfReKmCDOVB4xiZLGxvj4OJYsWYKRkRFD55AkCS0tLbDb7XL3ZEBpdSJLLGGky3C9bir1gNlyl5TES0iShLGxMfT29rK1WoA8Er29vUin03C5XHJ53kwmg2g0mtdrVEtwKWlGi3rwSuTCtGLx7W9/G7/+9a9LcS0MUzA0UY0IdEShoTG5IC8J19+uDkQhJRQKKTZ5Mw09qb6+0+lU1OOnUAefzyeH1BlpFupyuep2U6kH9ARCGkOityIcDstdugkqZU0Vwpg9kBJO1ZRI+U+n03UXOlZP98IwRjEdCjU6Oor//M//xEsvvYRDDz00qwv3nXfeadnFMYwZyDJopCIPYE6wNEu1lkecbKjDMMR/2+12wwnddrsdDodDjpem85JlOx6PIxQKyeVI8zXFosZfTHWSq9zl2NiY4v/UbVks3gCUdn2pRSg/xefzyV5dUtxI+WLLPsPUPqYVi3fffReHH344AGDjxo2Kn5XCAswwRhHDUajsLJUsLAekZCeTSS5kUOVEo1FTVaLa2tpkZSESicjCJMXKihVu1IKnFjw+ag+thnjAhKVd/JxCeNRlVic7YsNKddldtugzTP1gWrF4+eWXS3EdDFM0YslPshaWS6nwer045phjWJioEfr6+kxZlNWKAIW4kIAkCpZGzit6PpjaQB3eSEUA1E0XqbrRZG7ypYVewqroyeAO5QxT+xTUx4JhqpFAIICtW7dmVfrRappnNRQnzMJEbUAJpEZC5vx+v1xff3R0NKvPSSEx9NQjY3BwkCuHVRijnW5FwVjdDE0cG8lkkpVGDfQqbpEng7wZXEGJYWob08nbDFPNHHPMMVl5P1qVokScTm392kxoH+dU1BaBQAAdHR0Ih8N5j922bRt6e3vh8/kUSbuSJKG7u9uQciIi9lXJl4vBlB6jFaD0EnHpc3Ed6evrkwVorgqlpLe3F4lEQpFX4fF45F4gteLt5ffLMNqwYsHUHerE2HyCn164lFYPCz1lIxaL8QZTgxixjJIVWq0ESJJUkGIgljbm5omVp9CSoGrBkhK4XS6X7N0wU7J2skKKWXt7e01VUOL3yzDacCgUU3eUMjGWlA2x2RqFXrELv/YotzLocDjQ2trKseRVRKHhiyRYUn6N1nk43yobqp5W689FDI0zGk7HMJMBViyYuiIajSKVSsHlcqGpqSlvZ2Ux/8LhcMBmsxlK+CYFY2RkBOFwuC42ysmGXpUfI1D5WbNhUG63G+3t7bIgAnBFnFpFLBZB+TKiwsj5VtrUy3MR76Orq4s7bDPM/4dDoZi6gkrOOp1OzJ49O+/xYlJ3Op02pFSIORzpdFruJktx1UxtoBfCYKTJotvtlseZGdxut/zdHEZRu5BiKIayxWIxfqeTFO6wzTB7YMWCqQso3tnn88kLfCk2eJvNphAmKTGcBcXagzpnqzFSQYzyJFKplKk8CfKgsSBS24iNEcPhcE0mHzPWwR22GWYPHArF1AViyUJqUDU4OCiXMCwEyqMQkSQJbrdb9myEQiEA+jXameqFhIBCw6GI7du3Gz7W4XAoOncztYE6hl5depbfpXmi0ahcGSoUCvEzZJg6gT0WTF0gCvTkNSgmidvhcGDmzJlyWIxYInRoaAhOp1OxGbLFqjYp1sNks9lMdfCmLs3d3d1FfS9TXtQeSZ7vxUNhq8lkkj29DFNHsGLB1AWBQEAOSSAlQy/UxQjpdBqxWExRBYpCXhwOB4c91Qk0RoyME0ruF5EkyVS/E4L7V9QWHLpmPeLc4+fKMPWDTdIq1s/kZXh4WC41On369EpfDqMiGo2ip6cHmUwGXq8XIyMjpiv4iPj9frlDMoVF+Hw+LhtaR3R3dxcs8Itlh7Wg3ByqVCaOJ4ZhGIapZszIvJxjwdQlfX19chJuvpKzRti+fbssePr9fixZsoRLDNYBYuw8Cfq5lAun06lZOWxoaCinciFJEtLpNI455hhrLpxhGIZhqhAOhWLqEp/PJ/+bKjcVA4VGAXu6bGuFR6i78TLVC/WxIOXQiMciVznifApsJpPh3AqGYRimrmHFgqlLKHHb4/GgtbVV/lzsUWCz2QznYPj9fkVZUfJSqBM4uexs7SC+o2AwWJa8B86tYBiGYeoZViyYukT0JojVocQeBYcccgg6Ojry9iEIh8Nob29He3t7VoJ4ru9lqht6V1SeWKz8ZTWk0JrpecEwDMMwtQbnWDB1iVhbnvpZqPtSRCIRDA4Oorm5GQMDA0ilUll9K9RhVPlq1nNN+9pB/a6o6R2RLyHbDAsWLOBxwTAMw9Q9rFgwdYmYlEseC60CaLFYDP39/bInw+/3Ix6Pw+fzIRaLIZ1OyyEzYkMspn6gCmKEw+FAa2urZWFL5BFhGIZhmHqHQ6GYukTMdaCQF7/fD4/HA6/XqzhWDI/atm0bEomEQqhMpVLo7e3l3Ik6RawgBgDLli3DwMCAZefnhngMwzDMZIE9FkxdEgwGFR4G0WLc1dWl+3tqr4bL5ZL7X3Ajp/pC7EcyNjaGTCYDv9+PaDRaVM8TLWKxGPetYBiGYeoeViyYukQr10EUJBOJhObvNTY2yrH2fr8fzc3NiEQiACZ6GHBIS/1AXi0AOPHEEwHsKUFbCrq6ujiUjmEYhqlrOBSKqVvUPSVIkIzH4wiHw5q/k06n5RK0AwMDCAQCeStBMbWJVgUvMdRNLE1cKDabDcBE3gaH0jEMwzD1DnssmLpFzLMIBAJZ4VEAsqzTYghMMplENBrlSk91itZ7FcfIxo0bi/4OSZJkJZbOyzAMwzD1CnssmLolGAzC5XIhlUrJCoLY0C4QCMjeCT3rdG9vb9mul6kOKFlfq4pYIeg1U2QYhmGYeoMVC6ZuCQQCcDqdSCaT6O3tVYRFEaFQCC6XCw6Hw5LQF6a26evrQzKZRDKZLGo8+P1+DqFjGIZhJh0sSTF1DcXRA1DEuFP+BQBZ+XA4HHJZWpfLBZfLhZaWFk2FhKlPgsGgrFC0tbVh+fLluvk4uYjH4+ylYBiGYSYdnGPB1DUURy82zAMmQpySySQikQi8Xi8SiQRaWlqySoJ2dXUp8jSY+iYQCCiS/GncmGV8fFwOv2MYhmGYyQJ7LJhJQS7r8dDQEICJKlBq74RW5SCmvvH5fPLfpGSYDYsSO7YzDMMwzGSh4orFPffcgzlz5qCxsRGLFi3CW2+9pXtsMpnEjTfeiAMOOACNjY047LDD8MILLyiO2blzJy6++GIEAgF4PB4cffTR+POf/6w4RpIkXHvttZg5cyY8Hg+WLl2KDz74oCT3x1QX0WgUnZ2dSKfT8mc2m00zXArIrZAw9Uk8HgcwoWimUim4XC4sWLAAfr8fALI6t+vByijDMAwz2aioYvHkk0/i0ksvxXXXXYfu7m4cdthh6OjowI4dOzSPv+aaa/CrX/0K//Ef/4FNmzbhggsuwFe/+lW8/fbb8jHf/va38eKLL+LRRx9FJBLBCSecgKVLl+Ljjz+Wj7n11lvx85//HPfddx/efPNNTJkyBR0dHXJjNKZ+oeTcTCYjf+Z0OrFkyRKEQiH2TjCKvJxkMik3RmxubobH48Hs2bPhcDh0f9/hcCAcDrMyyjAMw0w6bJJVNRULYNGiRfjc5z6HX/ziFwCATCaDWbNm4fvf/z6uuuqqrOP9fj+uvvpqfO9735M/O/300+HxePDYY48hkUhg2rRp+J//+R+cdNJJ8jFHHHEEli1bhptuugmSJMHv9+Oyyy7D5ZdfDmAiFKa1tRWrVq3CmWeeaejah4eH4fV6MTQ0hOnTpxfzGJgyEo1G5RKyTU1NGBoagt/vz8qtYCY34jgJhUIIBAJyvo3H40EqlVL0PFHjcrnk32MYhmGYWsaMzFsxj8X4+DjWr1+PpUuX7rkYux1Lly7FunXrNH9nbGwMjY2Nis88Hg9ef/11ABP159PpdM5jNm/ejP7+fsX3er1eLFq0SPd7mdpHrALV0dGBjo4OjI+PA9gT+sIwwMRYiUQiCm8FsMeT4fP5FKF0WiSTSc6xYBiGYSYdFVMsPvnkE6TTabS2tio+b21tRX9/v+bvdHR04M4778QHH3yATCaDF198EU8//TS2bdsGAJg2bRoWL16Mf//3f0csFkM6ncZjjz2GdevWycfQuc18LzCh1AwPDyv+MLWD2IWb4MRsRgv1GCEo3yYejytC6URcLpec6E1J4AzDMAwzWah48rYZ7r77bhx44IEIhUJwu9248MILsWLFCkXFlkcffRSSJGHfffdFQ0MDfv7zn+NrX/ta0c3PbrnlFni9XvnPrFmzir0dpoxoKRGcmM1oQWNFK08iGo0ilUrBbrfLXdsJh8OBUCgk51/EYjHuf8IwDMNMKiqmWOy9995wOBzYvn274vPt27ejra1N83daWlrwzDPPYPfu3XIM9NSpU7H//vvLxxxwwAF49dVXsWvXLnz00Ud46623kEwm5WPo3Ga+FwBWrlyJoaEh+c9HH31U0H0zlSGfEkGhUiwIMrnGipj8HwqFFD9zu90AoMi94HAohmEYZjJRMcXC7XbjiCOOkOPegYnk7a6uLixevDjn7zY2NmLfffdFKpXCU089hX/8x3/MOmbKlCmYOXMm4vE4Ojs75WPmzp2LtrY2xfcODw/jzTffzPm9DQ0NmD59uuIPUz9ohUoxjBrR49XX1yeXoAUm1jT1+OEwO4ZhGGYyUdHO25deeinOPfdcHHnkkfj85z+Pn/3sZ9i9ezdWrFgBADjnnHOw77774pZbbgEAvPnmm/j4449x+OGH4+OPP8b111+PTCaDK6+8Uj5nZ2cnJEnCvHnz0NfXhyuuuAKhUEg+p81mw8UXX4ybbroJBx54IObOnYv/83/+D/x+P0499dSyPwOmOggGg4rO3AyjBXkxaKwEAgHEYjEAE9XlwuEwenp6kMlkcpakZRiGYZh6pKKKxRlnnIGBgQFce+216O/vx+GHH44XXnhBTqzeunWrIjdidHQU11xzDf72t79h6tSpOPHEE/Hoo49ixowZ8jFDQ0NYuXIl/v73v2OvvfbC6aefjh//+MeKeOgrr7wSu3fvxr/+67/is88+wxe/+EW88MILWdWkmMlDIBDgXAvGEOqx4vf7EYvFZO8FJXZT920eVwzDMMxkoaJ9LGoZ7mPBMJObaDSq8HL19fVhbGxMViy4lwXDMAxTD5iReSvqsWAYhql2RAVCVBIoLycSicButytK0LpcLnR0dFTichmGYRimYtRUuVmGYZhyo07spwpiYp8KUirsdjs8Hk9WxSiGYRiGmQywYsFMarjMLJMPdQ8UUjTEvArC4XDIx/G4YhiGYSYbrFgwkxouM8vkQ93XQqwcFo/HEQ6H4fF44HK5kEwm0dfXx+OKYRiGmZRwjgUzqeEys4xZxJKzPp8vK4Fb698MwzAMMxngqlAFwlWhGIbp6upCIpEAAITDYa4AxTAMw9QdZmReDoViGIYpEHUnboZhGIaZzHAoFMMwTIGoO3EzDMMwzGSGFQuGYZgi4K7tDMMwDDMBh0IxDMMUCJcrZhiGYZg9sGLBMAxTIFxWlmEYhmH2wIoFwzBMgaib5zEMwzDMZIZzLBiGYQqE8ysYhmEYZg/ssWAYhmEYhmEYpmhYsWAYhmEYhmEYpmhYsWAYhmEYhmEYpmhYsWAYhmEYhmEYpmhYsWAYhmEYhmEYpmhYsWAYhmEYhmEYpmhYsWAYhmEYhmEYpmi4j0WBSJIEABgeHq7wlTAMwzAMwzBMaSBZl2TfXLBiUSA7d+4EAMyaNavCV8IwDMMwDMMwpWXnzp3wer05j7FJRtQPJotMJoNYLIZp06bBZrNV+nJqguHhYcyaNQsfffQRpk+fXunLYcoAv/PJB7/zyQe/88kFv+/JhyRJ2LlzJ/x+P+z23FkU7LEoELvdjv3226/Sl1GTTJ8+nRejSQa/88kHv/PJB7/zyQW/78lFPk8FwcnbDMMwDMMwDMMUDSsWDMMwDMMwDMMUDSsWTNloaGjAddddh4aGhkpfClMm+J1PPvidTz74nU8u+H0zueDkbYZhGIZhGIZhioY9FgzDMAzDMAzDFA0rFgzDMAzDMAzDFA0rFgzDMAzDMAzDFA0rFowpXnvtNZx88snw+/2w2Wx45plndI+94IILYLPZ8LOf/Uzz52NjYzj88MNhs9mwYcMGxc/effddHHPMMWhsbMSsWbNw6623WncTjCnK8c63bNkCm82W9edPf/qTtTfDGMKKdz5nzpys9/mTn/xEcQzP8+qhHO+c53n1YNW6/txzz2HRokXweDzw+Xw49dRTFT/funUrTjrpJDQ1NWGfffbBFVdcgVQqZe3NMFUFKxaMKXbv3o3DDjsM99xzT87jfv/73+NPf/oT/H6/7jFXXnml5s+Hh4dxwgknIBAIYP369bjttttw/fXX4z//8z+Lvn7GPOV458RLL72Ebdu2yX+OOOKIgq+bKRyr3vmNN96oeJ/f//735Z/xPK8uyvHOCZ7nlceK9/3UU0/h7LPPxooVK/DOO+/gj3/8I77+9a/LP0+n0zjppJMwPj6ON954Aw8//DBWrVqFa6+91vL7YaoIiWEKBID0+9//Puvzv//979K+++4rbdy4UQoEAtJdd92VdcyaNWukUCgk9fT0SACkt99+W/7ZL3/5S8nn80ljY2PyZz/84Q+lefPmleAuGDOU6p1v3rw56zOmOij0neuNA4LnefVSqnfO87w6KeR9J5NJad9995Xuv/9+3fOuWbNGstvtUn9/v/zZvffeK02fPl0x75n6gj0WjKVkMhmcffbZuOKKK7BgwQLNY7Zv347zzz8fjz76KJqamrJ+vm7dOhx77LFwu93yZx0dHXj//fcRj8dLdu1MYVjxzolTTjkF++yzD774xS/i2WefLdUlM0Vi5J0DwE9+8hM0Nzdj4cKFuO222xQhEDzPawsr3jnB87z6yfe+u7u78fHHH8Nut2PhwoWYOXMmli1bho0bN8rHrFu3DuFwGK2trfJnHR0dGB4eRk9PT1nugyk/zkpfAFNf/PSnP4XT6cRFF12k+XNJknDeeefhggsuwJFHHoktW7ZkHdPf34+5c+cqPqOFqb+/Hz6fz/LrZgrHinc+depU3HHHHfjCF74Au92Op556CqeeeiqeeeYZnHLKKSW+A8Ys+d45AFx00UVob2/HXnvthTfeeAMrV67Etm3bcOeddwLgeV5rWPHOeZ7XDvne99/+9jcAwPXXX48777wTc+bMwR133IEvfelL+Otf/4q99toL/f39CqUCUM5xpj5hxYKxjPXr1+Puu+9Gd3c3bDab5jH/8R//gZ07d2LlypVlvjqmFFj1zvfee29ceuml8v8/97nPIRaL4bbbbmOBo8ow8s4BKN7noYceCrfbje985zu45ZZbuGNvjWHVO+d5XhsYed+ZTAYAcPXVV+P0008HADz00EPYb7/98Lvf/Q7f+c53yna9THXBoVCMZaxduxY7duzA7Nmz4XQ64XQ6EY1Gcdlll2HOnDkAgP/93//FunXr0NDQAKfTiWAwCAA48sgjce655wIA2trasH37dsW56f9tbW3luyEmL1a9cy0WLVqEvr6+ctwGYwIj71yLRYsWIZVKyR4rnue1g1XvXO8YnufVhZH3PXPmTADA/Pnz5d9raGjA/vvvj61btwLgOT5ZYY8FYxlnn302li5dqviso6NDrhoBAD//+c9x0003yT+PxWLo6OjAk08+iUWLFgEAFi9ejKuvvhrJZBIulwsA8OKLL2LevHkcHlFlWPXOtdiwYYO8eTHVg5F3rsWGDRtgt9uxzz77AOB5XktY9c71juF5Xl0Yed9HHHEEGhoa8P777+OLX/wiACCZTGLLli0IBAIAJub4j3/8Y+zYsUMeAy+++CKmT5+uUEiY+oIVC8YUu3btUliXNm/ejA0bNmCvvfbC7Nmz0dzcrDje5XKhra0N8+bNAwDMnj1b8fOpU6cCAA444ADst99+AICvf/3ruOGGG/Ctb30LP/zhD7Fx40bcfffduOuuu0p5a4wO5XjnDz/8MNxuNxYuXAgAePrpp/Hggw/i/vvvL9l9MfoU+87XrVuHN998E8cddxymTZuGdevW4ZJLLsFZZ50lKw08z6uLcrxznufVQ7Hve/r06bjgggtw3XXXYdasWQgEArjtttsAAP/8z/8MADjhhBMwf/58nH322bj11lvR39+Pa665Bt/73vc4HLKeqXRZKqa2ePnllyUAWX/OPfdczeMLLT/4zjvvSF/84helhoYGad9995V+8pOfWHcTjCnK8c5XrVolHXzwwVJTU5M0ffp06fOf/7z0u9/9ztobYQxT7Dtfv369tGjRIsnr9UqNjY3SwQcfLN18883S6Oio4vd4nlcP5XjnPM+rByvW9fHxcemyyy6T9tlnH2natGnS0qVLpY0bNyqO2bJli7Rs2TLJ4/FIe++9t3TZZZdJyWSyRHfFVAM2SZKkMukwDMMwDMMwDMPUKZy8zTAMwzAMwzBM0bBiwTAMwzAMwzBM0bBiwTAMwzAMwzBM0bBiwTAMwzAMwzBM0bBiwTAMwzAMwzBM0bBiwTAMwzAMwzBM0bBiwTAMwzAMwzBM0bBiwTAMwzAMwzBM0bBiwTAMw5jmlVdegc1mw2effVaR7+/q6sLBBx+MdDoNALj++utx+OGHF3Qum82GZ555RvfnW7Zsgc1mw4YNGwAAmzZtwn777Yfdu3cX9H0MwzD1CisWDMMwTE6+9KUv4eKLL1Z8dvTRR2Pbtm3wer0VuaYrr7wS11xzDRwOh+HfOe+883Dqqadmfb5t2zYsW7bM8Hnmz5+Po446Cnfeeafh32EYhpkMsGLBMAzDmMbtdqOtrQ02m63s3/3666/jww8/xOmnn27J+dra2tDQ0GDqd1asWIF7770XqVTKkmtgGIapB1ixYBiGYXQ577zz8Oqrr+Luu++GzWaDzWbDli1bskKhVq1ahRkzZmD16tWYN28empqa8E//9E8YGRnBww8/jDlz5sDn8+Giiy6Sw5cAYGxsDJdffjn23XdfTJkyBYsWLcIrr7yS85qeeOIJfPnLX0ZjY2PWz371q19h1qxZaGpqwr/8y79gaGgIwESo1MMPP4z/+Z//ke+DvkcdCvXWW29h4cKFaGxsxJFHHom3334763u+/OUv49NPP8Wrr75q7oEyDMPUMc5KXwDDMAxTvdx9993461//ikMOOQQ33ngjAKClpQVbtmzJOnZkZAQ///nP8cQTT2Dnzp047bTT8NWvfhUzZszAmjVr8Le//Q2nn346vvCFL+CMM84AAFx44YXYtGkTnnjiCfj9fvz+97/HV77yFUQiERx44IGa17R27Vp8/etfz/q8r68Pv/3tb/GHP/wBw8PD+Na3voV/+7d/w+OPP47LL78c7733HoaHh/HQQw8BAPbaa6+sc+zatQvLly/Hl7/8ZTz22GPYvHkzfvCDH2Qd53a7cfjhh2Pt2rVYsmSJ4efJMAxTz7BiwTAMw+ji9XrhdrvR1NSEtra2nMcmk0nce++9OOCAAwAA//RP/4RHH30U27dvx9SpUzF//nwcd9xxePnll3HGGWdg69ateOihh7B161b4/X4AwOWXX44XXngBDz30EG6++WbN74lGo/LxIqOjo3jkkUew7777AgD+4z/+AyeddBLuuOMOtLW1wePxYGxsLOd9/PrXv0Ymk8EDDzyAxsZGLFiwAH//+9/x3e9+N+tYv9+PaDSa85kwDMNMJlixYBiGYSyhqalJVioAoLW1FXPmzMHUqVMVn+3YsQMAEIlEkE6ncdBBBynOMzY2hubmZt3vSSQSmmFQs2fPlpUKAFi8eDEymQzef//9vEoR8d577+HQQw9VnH/x4sWax3o8HoyMjBg6L8MwzGSAFQuGYRjGElwul+L/NptN87NMJgNgIuzI4XBg/fr1WdWdRGVEzd577414PG7RVRfOp59+qlCkGIZhJjucvM0wDMPkxO12KxKurWLhwoVIp9PYsWMHgsGg4k8uD8PChQuxadOmrM+3bt2KWCwm//9Pf/oT7HY75s2bZ/g+Dj74YLz77rsYHR1VnEeLjRs3YuHChTnPxzAMM5lgxYJhGIbJyZw5c/Dmm29iy5Yt+OSTT2SPQ7EcdNBB+MY3voFzzjkHTz/9NDZv3oy33noLt9xyC5577jnd3+vo6MDrr7+e9XljYyPOPfdcvPPOO1i7di0uuugi/Mu//IuspMyZMwfvvvsu3n//fXzyySdIJpNZ5/j6178Om82G888/H5s2bcKaNWtw++23Zx23ZcsWfPzxx1i6dGkRT4BhGKa+YMWCYRiGycnll18Oh8OB+fPno6WlBVu3brXs3A899BDOOeccXHbZZZg3bx5OPfVU/PnPf8bs2bN1f+cb3/gGenp68P777ys+DwaDOO2003DiiSfihBNOwKGHHopf/vKX8s/PP/98zJs3D0ceeSRaWlrwxz/+MevcU6dOxR/+8AdEIhEsXLgQV199NX76059mHfeb3/wGJ5xwAgKBQBF3zzAMU1/YJEmSKn0RDMMwDGOGK664AsPDw/jVr35V1HnGxsbQ2NiIF1980bD3YXx8HAceeCB+/etf4wtf+EJR388wDFNPsMeCYRiGqTmuvvpqBAKBosKyhoeH8Zvf/AZ2ux2hUMjw723duhU/+tGPWKlgGIZRwR4LhmEYZlJyySWX4Ne//jWuuOIKXH755ZW+HIZhmJqHFQuGYRiGYRiGYYqGQ6EYhmEYhmEYhikaViwYhmEYhmEYhikaViwYhmEYhmEYhikaViwYhmEYhmEYhikaViwYhmEYhmEYhikaViwYhmEYhmEYhikaViwYhmEYhmEYhikaViwYhmEYhmEYhikaViwYhmEYhmEYhima/wePG5zcvm5C2gAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "import lightkurve as lk\n", + "import numpy as np\n", + "\n", + "# single sector\n", + "lc = lk.search_lightcurve(\"TOI 451\", author=\"SPOC\", exptime=120)[1].download()\n", + "\n", + "# masking nans\n", + "time = lc.time.to_value(\"btjd\")\n", + "flux = lc.pdcsap_flux.to_value().filled(np.nan)\n", + "error = lc.flux_err.to_value().filled(np.nan)\n", + "mask = np.isnan(flux) | np.isnan(error) | np.isnan(time)\n", + "time = time[~mask].astype(float)\n", + "flux = flux[~mask].astype(float)\n", + "error = error[~mask].astype(float)\n", + "\n", + "# normalize\n", + "flux_median = np.median(flux)\n", + "flux /= flux_median\n", + "error /= flux_median\n", + "\n", + "# plot\n", + "plt.figure(figsize=(8, 3))\n", + "plt.plot(time, flux, \".\", c=\"0.7\", ms=2, label=\"PDC-SAP flux\")\n", + "plt.xlabel(\"time (btjd)\")\n", + "plt.ylabel(\"normalized flux\")\n", + "plt.legend()\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and rerun the GP optimization." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import jax.numpy as jnp\n", + "from nuance.core import gp_model\n", + "from nuance.utils import minimize\n", + "from tinygp import kernels, GaussianProcess\n", + "\n", + "# rotation period\n", + "star_period = 4.9555\n", + "\n", + "# gp\n", + "initial_params = {\n", + " \"log_period\": jnp.log(star_period),\n", + " \"log_Q\": jnp.log(100),\n", + " \"log_sigma\": jnp.log(1e-1),\n", + " \"error\": np.mean(error),\n", + "}\n", + "\n", + "\n", + "def build_gp(params, time):\n", + "\n", + " kernel = kernels.quasisep.SHO(\n", + " jnp.exp(params[\"log_sigma\"]),\n", + " jnp.exp(params[\"log_period\"]),\n", + " jnp.exp(params[\"log_Q\"]),\n", + " )\n", + "\n", + " return GaussianProcess(kernel, time, diag=params[\"error\"] ** 2, mean=1.0)\n", + "\n", + "\n", + "# optimization\n", + "mu, nll = gp_model(time, flux, build_gp)\n", + "params = minimize(nll, initial_params)\n", + "gp = build_gp(params, time)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Estimating the SNR" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The analytical injection-recovery presented in this tutorial will consist in computing the recovered SNR of injected transit signals given a range of planetary radii and periods.\n", + "\n", + "We will compute the duration and depths assuming a circular orbit at 0 impact parameter. Let's instantiate a [Star](nuance.Star) object to hold the stellar parameters of TOI 451" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "from nuance import core, Star\n", + "\n", + "star = Star(radius=0.950, mass=0.879)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We then generate `n` random epochs and `jax.vmap` the [core.solve](nuance.core.solve) function" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "n_epochs = 10\n", + "\n", + "t0s = np.random.uniform(0, 12, n_epochs)\n", + "solve_function = jax.vmap(\n", + " jax.vmap(jax.jit(core.solve(time, flux, gp)), in_axes=(None, None, None, 0)),\n", + " in_axes=(0, None, None, None),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```{note}\n", + "The [core.solve](nuance.core.solve) function returns a function that computes the depth of a transt given its epoch and duration, it has the signature `function(epoch, duration, period=None, depth=None) -> (log_likelihood, weights, variance)`. When the `depth` parameter is set, the function computes the depth of a transit **injected** with this `depth`. Usually `depth=None`, as this feature is only useful for injection-recovery estimates.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now define `snr` the function that will:\n", + "1. compute the duration of the transit given its period\n", + "2. Compute the `N` depths of the transit given a set of `N` radii\n", + "3. Compute the `n_epochs` fitted depths of the injected transits\n", + "4. Compute the SNR of the injected transits" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "@jax.jit\n", + "def snr(period, radii):\n", + " duration = star.transit_duration(period) # 1\n", + " depths = star.transit_depth(radii) # 2\n", + " _depth, _var, _ = solve_function(t0s, duration, period, depths).T # 3\n", + " return _depth / jnp.sqrt(_var) # 4" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's compute the SNRs for a set of periods and radii" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 83/83 [00:11<00:00, 7.00it/s]\n" + ] + } + ], + "source": [ + "from tqdm import tqdm\n", + "\n", + "# We downsample the period grid to speed up the calculation here\n", + "periods = star.period_grid(np.ptp(time), oversampling=0.1)\n", + "radii = np.linspace(0.5, 10.0, 50)\n", + "\n", + "snrs = np.array([snr(period, radii) for period in tqdm(periods)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now plot the resulting recovery assuming a lower limit on the SNR for a signal to be detected" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAG2CAYAAABh8Lw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+XklEQVR4nO3deVhWdf7/8RegAirggoALKrmUu6ZJZH6zJCnNMr81tkyaW5vmQpnLhLaKWjqamZbl0pRpU6NTWabDNy0NN9DUXDMTMkFNA4ES4j6/P/pxT7eAwuGGw+F+Pq7rXCOf+5zPeR9svN/X+7McL8MwDAEAANiIt9UBAAAAlBYJDAAAsB0SGAAAYDskMAAAwHZIYAAAgO2QwAAAANshgQEAALZDAgMAAGyHBAYAANgOCQwAALAdEhgAADzMl19+qf79+6tRo0by8vLSmjVrLnvNxo0bdfXVV8vX11ctW7bUsmXLyj3OSyGBAQDAw2RnZ6tTp05asGBBic4/duyY+vXrpxtvvFG7d+/WuHHjNGLECH3++eflHGnxvHiZIwAAnsvLy0urV6/WgAEDij1n4sSJWrt2rfbt2+dsu+eee/TLL79o3bp1FRBlYdUsuatFHA6HfvrpJwUEBMjLy8vqcAAAlZRhGDp//rwaNWokb+/yG6z47bfflJub65a+DMMo9N3m6+srX1/fMvedmJio6Ohol7aYmBiNGzeuzH2bZlQSmzZtMm677TajYcOGhiRj9erVLp87HA4jLi7OCAsLM/z8/IzevXsbhw8fLtU9UlNTDUkcHBwcHBwlOlJTU934Tefq119/NcLCwtwWa+3atQu1TZs27bJxSIW/cy/WqlUrY/r06S5ta9euNSQZOTk5ZfgtmFdpKjAF43HDhg3TwIEDC30+a9YsvfLKK1q+fLkiIiIUFxenmJgY7d+/X35+fiW6R0BAgCQpNTVVgYGBbo0fAFB1ZGZmKjw83Pm9UR5yc3OVlpbmlu+kgngv7ssd1ZfKqtIkMLfeeqtuvfXWIj8zDENz587V008/rTvuuEOS9Pbbbys0NFRr1qzRPffcU6J7FJTWAgMDSWAAAJdVEdMN3PmdVF7fb2FhYUpPT3dpS09PV2BgoPz9/d1+v5KwxSqkY8eOKS0tzWX8LSgoSJGRkUpMTCz2ugsXLigzM9PlAACgMjEMwy1HeYqKilJCQoJL24YNGxQVFVWu970UWyQwaWlpkqTQ0FCX9tDQUOdnRYmPj1dQUJDzCA8PL9c4AQCwg6ysLO3evVu7d++W9EehYPfu3UpJSZEkTZ48WYMHD3ae/8gjj+j777/XU089pYMHD+q1117T+++/r/Hjx1sRviSbJDBmTZ48WRkZGc4jNTXV6pAAALDczp071aVLF3Xp0kWSFBsbqy5dumjq1KmSpJMnTzqTGUmKiIjQ2rVrtWHDBnXq1EmzZ8/Wm2++qZiYGEvilyrRHJhLCQsLk/THeFvDhg2d7enp6ercuXOx17lr+RgAAOXFHUNApb2+V69el7ymqF12e/XqpV27dpU2tHJjiwpMRESEwsLCXMbfMjMztW3bNkvH3wAAgDUqTQUmKytL3333nfPngvG4evXqqWnTpho3bpxeeOEFtWrVyrmMulGjRpfcORAAAFRNlSaB2blzp2688Ubnz7GxsZKkIUOGaNmyZXrqqaeUnZ2thx56SL/88ouuv/56rVu3rsR7wAAAUBlZMYRUFXjUu5AyMzMVFBSkjIwM9oEBABSrIr4vCu5x9uxZt2xkV69ePY/6frPFHBgAAIA/qzRDSAAAeCKGkMyhAgMAAGyHBAYAANgOQ0gAAFiIISRzqMAAAADboQIDAICFqMCYQwUGAADYDgkMAACwHYaQAACwEENI5lCBAQAAtkMCAwAAbIchJAAALMQQkjlUYAAAgO2QwAAAANshgQEAALZDAgMAAGyHSbwAAFiISbzmUIEBAAC2QwUGAAALUYExhwoMAACwHRIYAABgOwwhAQBgIYaQzKECAwAAbIcEBgAA2A5DSAAAWIghJHOowAAAANshgQEAALbDEBIAABZiCMkcKjAAAMB2SGAAAIDtMIQEAICFGEIyhwoMAACwHRIYAABgOwwhAQBgIYaQzKECAwAAbIcEBgAA2A5DSAAAWIghJHNIYAAAsBAJjDkMIQEAANshgQEAALbDEBIAABZiCMkcKjAAAMB2qMAAAGAxT6yglBUVGAAAYDtUYAAAsBBzYMyhAgMAAGyHCgwAABaiAmMOCQwAABYigTGHISQAAGA7JDAAAMB2GEICAMBCDCGZQwUGAADYDhUYAAAsRAXGHCowAADAdkhgAACA7TCEBACAhRhCMocEBgAAC5HAmEMCAwCAhUhgzGEODAAAHmjBggVq3ry5/Pz8FBkZqe3bt1/y/Llz5+rKK6+Uv7+/wsPDNX78eP32228VFG1hJDAAAHiYVatWKTY2VtOmTVNycrI6deqkmJgYnTp1qsjzV6xYoUmTJmnatGk6cOCA3nrrLa1atUpTpkyp4Mj/iwQGAAALFQwhlfUojTlz5mjkyJEaOnSo2rZtq0WLFqlmzZpasmRJked//fXX6tGjh+677z41b95cffr00b333nvZqk15IoEBAKCKyMzMdDkuXLhQ6Jzc3FwlJSUpOjra2ebt7a3o6GglJiYW2e91112npKQkZ8Ly/fff69NPP1Xfvn3L50FKgEm8AABYyJ2TeMPDw13ap02bpmeeecal7cyZM8rPz1doaKhLe2hoqA4ePFhk//fdd5/OnDmj66+/XoZh6Pfff9cjjzxi6RASCQwAABZyZwKTmpqqwMBAZ7uvr2+Z+i2wceNGTZ8+Xa+99poiIyP13XffaezYsXr++ecVFxfnlnuUlm2GkPLz8xUXF6eIiAj5+/urRYsWev755z1y6RgAAEUJDAx0OYpKYIKDg+Xj46P09HSX9vT0dIWFhRXZb1xcnB544AGNGDFCHTp00J133qnp06crPj5eDoejXJ7lcmyTwMycOVMLFy7Uq6++qgMHDmjmzJmaNWuW5s+fb3VoAADYRo0aNdS1a1clJCQ42xwOhxISEhQVFVXkNTk5OfL2dk0ZfHx8JFm3B41thpC+/vpr3XHHHerXr58kqXnz5nrvvfcsnQENAEBZWbGRXWxsrIYMGaJu3bqpe/fumjt3rrKzszV06FBJ0uDBg9W4cWPFx8dLkvr37685c+aoS5cuziGkuLg49e/f35nIVDTbJDDXXXed3njjDR0+fFitW7fWN998o82bN2vOnDnFXnPhwgWXGdiZmZkVESoAAJXaoEGDdPr0aU2dOlVpaWnq3Lmz1q1b55zYm5KS4lJxefrpp+Xl5aWnn35aJ06cUIMGDdS/f3+9+OKLVj2CvAybTCJxOByaMmWKZs2aJR8fH+Xn5+vFF1/U5MmTi73mmWee0bPPPluoPSMjw2WSEwAAf5aZmamgoKBy/b4ouMeePXsUEBBQpr7Onz+vjh07etT3m23mwLz//vt69913tWLFCiUnJ2v58uV6+eWXtXz58mKvmTx5sjIyMpxHampqBUYMAADKi22GkCZMmKBJkybpnnvukSR16NBBx48fV3x8vIYMGVLkNb6+vm5bQgYAQHngZY7m2KYCU9wMaKuWbwEAAOvYpgJTMFmoadOmateunXbt2qU5c+Zo2LBhVocGAIBpVGDMsU0CM3/+fMXFxemxxx7TqVOn1KhRIz388MOaOnWq1aEBAGAaCYw5tklgAgICNHfuXM2dO9fqUAAAgMVsMwcGAACggG0qMAAAVEUMIZlDBQYAANgOFRgAACzmiRWUsqICAwAAbIcEBgAA2A5DSAAAWIhJvOaQwAAAYCESGHNIYAAAsBAJjDnMgQEAALZDBQYAAAtRgTGHBAYAAAuRwJjDEBIAALAdKjAAAFiICow5VGAAAIDtkMAAAADbYQgJAAALMYRkDgkMAAAWIoExhyEkAABgO1RgAACwEBUYc0hgAACwEAmMOSQwAABYiATGHObAAAAA26ECAwCAhajAmEMCAwCAhUhgzGEICQAA2A4VGAAALEQFxhwqMAAAwHZIYAAAgO0whAQAgIUYQjKHBAYAAAuRwJjDEBIAALAdKjAAAFiICow5JDAAAFjMExOQsiKBAQDAQlRgzGEODAAAsB0qMAAAWIgKjDkkMAAAWIgExhyGkAAAgO2QwAAAANthCAkAAAsxhGQOFRgAAGA7VGAAALAQFRhzSGAAALAQCYw5DCEBAADboQIDAICFqMCYQwIDAICFSGDMIYEBAMBCJDDmMAcGAADYDhUYAAAsRAXGHCowAADAdkhgAACA7TCEBACAhRhCMocEBgAACzkcDjkcjjL34WkYQgIAALZDBQYAAAsxhGQOCQwAABYigTGHISQAACxUkMCU9SitBQsWqHnz5vLz81NkZKS2b99+yfN/+eUXjRo1Sg0bNpSvr69at26tTz/91OxjlxkVGAAAPMyqVasUGxurRYsWKTIyUnPnzlVMTIwOHTqkkJCQQufn5ubq5ptvVkhIiD744AM1btxYx48fV506dSo++P+PBAYAAAsZhlHmVUSlrcDMmTNHI0eO1NChQyVJixYt0tq1a7VkyRJNmjSp0PlLlizR2bNn9fXXX6t69eqSpObNm5cp5rJiCAkAAAu5cwgpMzPT5bhw4UKh++Xm5iopKUnR0dHONm9vb0VHRysxMbHIGD/66CNFRUVp1KhRCg0NVfv27TV9+nTl5+eXzy+lBEhgAACoIsLDwxUUFOQ84uPjC51z5swZ5efnKzQ01KU9NDRUaWlpRfb7/fff64MPPlB+fr4+/fRTxcXFafbs2XrhhRfK5TlKokxDSHl5eUpLS1NOTo4aNGigevXquSsuAAA8gjtXIaWmpiowMNDZ7uvrW6Z+CzgcDoWEhOiNN96Qj4+PunbtqhMnTuill17StGnT3HKP0ip1AnP+/Hm98847WrlypbZv367c3FwZhiEvLy81adJEffr00UMPPaRrrrmmPOIFAKBKcWcCExgY6JLAFCU4OFg+Pj5KT093aU9PT1dYWFiR1zRs2FDVq1eXj4+Ps61NmzZKS0tTbm6uatSoUab4zSjVENKcOXPUvHlzLV26VNHR0VqzZo12796tw4cPKzExUdOmTdPvv/+uPn366JZbbtGRI0fKK24AAGBCjRo11LVrVyUkJDjbHA6HEhISFBUVVeQ1PXr00Hfffecy2fjw4cNq2LChJcmLVMoEZseOHfryyy+1fft2xcXFKSYmRh06dFDLli3VvXt3DRs2TEuXLlVaWpoGDBigr776yq3BnjhxQn/9619Vv359+fv7q0OHDtq5c6db7wEAQEWyYh+Y2NhYLV68WMuXL9eBAwf06KOPKjs727kqafDgwZo8ebLz/EcffVRnz57V2LFjdfjwYa1du1bTp0/XqFGj3Pq7KI1SDSG99957JTrP19dXjzzyiKmAinPu3Dn16NFDN954oz777DM1aNBAR44cUd26dd16HwAAKpIVL3McNGiQTp8+ralTpyotLU2dO3fWunXrnBN7U1JS5O393xpHeHi4Pv/8c40fP14dO3ZU48aNNXbsWE2cOLFMcZeFl+HG/Yezs7N1+PBhhYeHKzg42F3dSpImTZqkLVu2lKmqk5mZqaCgIGVkZFx2jBAA4Lkq4vui4B4rVqxQzZo1y9RXTk6O7rvvPo/6fnPbRnbTp0/Xxo0bFRAQIB8fH9WuXVsLFiyQv7+/W/r/6KOPFBMTo7vvvlubNm1S48aN9dhjj2nkyJHFXnPhwgWXNfCZmZluiQUAAHfhXUjmuGUfmCVLlujcuXNav369jhw5ovfff1+33nqrxo0b547uJf2xBn3hwoVq1aqVPv/8cz366KMaM2aMli9fXuw18fHxLuvhw8PD3RYPAADuYNW7kOzOLQnMqlWr1K5dOz377LM6deqUnnvuOZ09e1Y7duxwR/eS/hjfu/rqqzV9+nR16dJFDz30kEaOHKlFixYVe83kyZOVkZHhPFJTU90WDwAA7kACY45bEpgaNWooLCxMXbt2lZ+fn7p166a2bdu69RfasGFDtW3b1qWtTZs2SklJKfYaX19f55r4kqyNBwAA5eOLL75wa39uSWBGjBihd955R7feeqsCAgLUt29f7dy5U7fccos7upf0xxr0Q4cOubQdPnxYzZo1c9s9AACoaAWrkMp6VHa33HKLWrRooRdeeMEtIyJuSWDuuOMO9e3bV71791aXLl3Up08fZWRkuPUdCePHj9fWrVs1ffp0fffdd1qxYoXeeOMNS9egAwBQVp4yhHTixAmNHj1aH3zwga644grFxMTo/fffV25urqn+3LqMWpLOnj2runXrysvLy53dSpI++eQTTZ48WUeOHFFERIRiY2MvuQrpYiyjBgCUREUuo162bJlbllE/+OCDtvl+S05O1tKlS537y913330aPny4OnXqVOI+3LaMukB5vtDxtttu02233VZu/QMAUNE8cRn11VdfrbCwMNWvX18zZszQkiVL9NprrykqKkqLFi1Su3btLtuHqSGk48ePa/369cW+dvunn34y0y0AAB7HU4aQJCkvL08ffPCB+vbtq2bNmunzzz/Xq6++qvT0dH333Xdq1qyZ7r777hL1VeoE5r333lPLli11yy236IorrtA//vEPSX9sOzxjxgxFRkaqadOmpe0WAABUYY8//rgaNmyohx9+WK1bt9auXbuUmJioESNGqFatWmrevLlefvllHTx4sET9lXoI6fnnn9fjjz+u4cOHa8qUKXr00Ud1+PBhzZgxQy1atFDv3r01adKkUj8YAACeyFOGkPbv36/58+dr4MCB8vX1LfKc4ODgEi+3LnUCc/ToUY0dO1bNmjXTggUL1LRpU23ZskV79uxRmzZtStsdAAAezYqXOVa0vLw8NWvWTNdee22xyYskVatWTTfccEOJ+iz1EFJeXp7z/UZNmjSRn5+fXn75ZZIXAABQpOrVq+vDDz90a5+mJvGuWLHCOUbl4+OjunXrujUoAAA8hadM4h0wYIDWrFnjtv5KPYTUs2dPTZs2TU888YTq1q2r3377TfPmzdN1112n9u3bq3Xr1qpWze2rswEAqLLskICUVatWrfTcc89py5Yt6tq1q2rVquXy+ZgxY0rVX6kzjU2bNkmSjhw5oqSkJCUnJys5OVlvv/22fvnlF9WoUUOtW7fWnj17Sts1AAAex1Mm8b711luqU6eOkpKSlJSU5PKZl5dX+ScwBVq1aqVWrVrpnnvucbYdO3ZMO3fu1K5du8x2CwAAqqBjx465tb9SJTApKSmX3OMlIiJCERERzk1oTpw4ocaNG5ctQgAAqjBPWIX0Z7m5uTp27JhatGhRpiknpZrEe8011+jhhx/Wjh07ij0nIyNDixcvVvv27d0+4xgAgKrGUybx5uTkaPjw4apZs6batWunlJQUSX9scDdjxoxS91eq1Gf//v168cUXdfPNN8vPz09du3ZVo0aN5Ofnp3Pnzmn//v369ttvdfXVV2vWrFnq27dvqQMCAABVz+TJk/XNN99o48aNuuWWW5zt0dHReuaZZ0q9CW6pEpj69etrzpw5evHFF7V27Vpt3rxZx48f16+//qrg4GDdf//9iomJUfv27UsVBAAAnspTJvGuWbNGq1at0rXXXisvLy9ne7t27XT06NFS92dq8Mnf31933XWX7rrrLjOXAwCA/89TEpjTp08rJCSkUHt2drZLQlNSpjayAwAAKI1u3bpp7dq1zp8LkpY333xTUVFRpe6PHecAALCQp1Rgpk+frltvvVX79+/X77//rnnz5mn//v36+uuvnXvMlQYVGAAALFSwjLqsR2V3/fXXa/fu3fr999/VoUMHrV+/XiEhIUpMTFTXrl1L3R8VGAAAUCFatGihxYsXu6UvKjAAAFjIU/aBiY6O1rJly5SZmemW/kwnML/++qtycnKcPx8/flxz587V+vXr3RIYAACewFMSmHbt2mny5MkKCwvT3XffrX//+9/Ky8sz3Z/pBOaOO+7Q22+/LUn65ZdfFBkZqdmzZ+uOO+7QwoULTQcEAIAn8ZQEZt68eTpx4oTWrFmjWrVqafDgwQoNDdVDDz1UsZN4k5OT1bNnT0nSBx98oNDQUB0/flxvv/22XnnlFbPdAgCAKsrb21t9+vTRsmXLlJ6ertdff13bt2/XTTfdVOq+TE/izcnJUUBAgCRp/fr1GjhwoLy9vXXttdfq+PHjZrsFAMCjeNrLHCUpLS1NK1eu1DvvvKM9e/aoe/fupe7DdAWmZcuWWrNmjVJTU/X555+rT58+kqRTp04pMDDQbLcAAHgUTxlCyszM1NKlS3XzzTcrPDxcCxcu1O23364jR45o69atpe7PdAVm6tSpuu+++zR+/Hj17t3buYve+vXr1aVLF7PdAgCAKig0NFR169bVoEGDFB8fr27dupWpP9MJzF133aXrr79eJ0+eVKdOnZztvXv31p133lmmoAAA8BSeshPvRx99pN69e8vb2z07uJRpI7uwsDCFhYW5tJkZxwIAwFN5SgJz8803S/rjpY6HDh2SJF155ZVq0KCBqf5MJzDPPffcJT+fOnWq2a4BAEAVk5OTo9GjR+vtt992Tjr28fHR4MGDNX/+fNWsWbNU/ZlOYFavXu3yc15eno4dO6Zq1aqpRYsWJDAAAJSQHSooZTV+/Hht2rRJH3/8sXr06CFJ2rx5s8aMGaMnnnii1HvImU5gdu3aVagtMzNTDz74IHNgAAAoIU8ZQvrwww/1wQcfqFevXs62vn37yt/fX3/5y19KncC49V1IgYGBevbZZxUXF+fObgEAgM3l5OQoNDS0UHtISIjLq4lKyu0vc8zIyFBGRoa7uwUAoErylH1goqKiNG3aNP3222/Otl9//VXPPvuscyuW0jA9hHTx6wIMw9DJkyf1j3/8Q7feeqvZbgEA8CieMoQ0b948xcTEqEmTJs7tV7755hv5+fnp888/L3V/phOYv//97y4/e3t7q0GDBhoyZIgmT55stlsAADyKpyQw7du315EjR/Tuu+/q4MGDkqR7771X999/v/z9/Uvdn+kE5tixY2YvBQAAHqhmzZoaOXKkW/oqVQITGxur559/XrVq1VJsbOwlz50zZ06ZAgMAwBN4SgUmPj5eoaGhGjZsmEv7kiVLdPr0aU2cOLFU/ZUqgdm1a5fy8vKcfy6Ol5dXqYIAAMBTecrbqF9//XWtWLGiUHu7du10zz33lG8C88UXXxT5ZwAAgEtJS0tTw4YNC7U3aNBAJ0+eLHV/bl9GDQAASs5TllGHh4dry5Ythdq3bNmiRo0albq/Us+BKSnmwAAAcHmeMgdm5MiRGjdunPLy8nTTTTdJkhISEvTUU0/piSeeKHV/pZ4D82fJycn6/fffdeWVV0qSDh8+LB8fH3Xt2rXUgQAAgKprwoQJ+vnnn/XYY48pNzdXkuTn56eJEyea2n7F9ByYOXPmKCAgQMuXL1fdunUlSefOndPQoUPVs2fPUgcCAIAn8pQKjJeXl2bOnKm4uDgdOHBA/v7+atWqlXx9fU31Z3oOzOzZsxUfH+9MXiSpbt26euGFFzR79myz3QIA4FE8ZQ5MgbS0NJ09e1YtWrSQr6+v6dhNJzCZmZk6ffp0ofbTp0/r/PnzZrsFAABV0M8//6zevXurdevW6tu3r3Pl0fDhw03NgTGdwNx5550aOnSo/vWvf+nHH3/Ujz/+qA8//FDDhw/XwIEDzXYLAIBH8ZQKzPjx41W9enWlpKSoZs2azvZBgwZp3bp1pe7P9KsEFi1apCeffFL33Xefc3O7atWqafjw4XrppZfMdgsAgEfxlDkw69ev1+eff64mTZq4tLdq1UrHjx8vdX+mE5iaNWvqtdde00svvaSjR49Kklq0aKFatWqZ7RIAAI/jKQlMdna2S+WlwNmzZ01N5C3zRna1atVSx44d1bFjR5IXAABQpJ49e+rtt992/uzl5SWHw6FZs2bpxhtvLHV/piswBfbv36+UlBTnmu4Ct99+e1m7BgCgyvOUCsxLL72km266STt37lRubq6eeuopffvttzp79myRO/RejukE5vvvv9edd96pvXv3ysvLy/nLK3iRY35+vtmuAQDwGJ6QwOTl5WnMmDH6+OOPtWHDBgUEBCgrK0sDBw7UqFGjinxH0uWYTmDGjh2riIgIJSQkKCIiQtu3b9fPP/+sJ554Qi+//LLZbgEAQBVTvXp17dmzR3Xr1tXf/vY3t/Rpeg5MYmKinnvuOQUHB8vb21ve3t66/vrrFR8frzFjxrglOAAAqjrDMORwOMp0VPYKjCT99a9/1VtvveW2/kxXYPLz8xUQECBJCg4O1k8//aQrr7xSzZo106FDh9wWIICKl5WVZXUIpVK7dm2rQ6hU7Pb3VxlV5O/QE4aQJOn333/XkiVL9J///Eddu3YttPCntC+BNp3AtG/fXt98840iIiIUGRmpWbNmqUaNGnrjjTd0xRVXmO0WAABUQfv27dPVV18t6Y+XP/9ZwfzZ0jCdwDz99NPKycmRJD333HO67bbb1LNnT9WvX1+rVq0y2y0AAB7FUyowf34htDuYmgOTl5enWbNmqX379pKkli1b6uDBgzpz5oxOnTqlm266ya1BAgBQVVn1KoEFCxaoefPm8vPzU2RkpLZv316i61auXCkvLy8NGDCg1Pd0J1MJTMFs4ovVq1fPVBkIAABUnFWrVik2NlbTpk1TcnKyOnXqpJiYGJ06deqS1/3www968skn1bNnzwqKtHimVyG5ezYxAACeyIoKzJw5czRy5EgNHTpUbdu21aJFi1SzZk0tWbKk2Gvy8/N1//3369lnn60Uc11Nz4Fx92xiAAA8kTvnwGRmZrq0+/r6FnrPUG5urpKSkjR58mRnm7e3t6Kjo5WYmFjsPZ577jmFhIRo+PDh+uqrr8oUrzuYTmDcPZsYAMxi2bCrI0eOWB2C7dl1GXV4eLhL+7Rp0/TMM8+4tJ05c0b5+fkKDQ11aQ8NDdXBgweL7H/z5s166623tHv37jLF6U6mExh3zyYGAABlk5qaqsDAQOfPZt7yfLHz58/rgQce0OLFixUcHFzm/tylzC9zBAAA5rmzAhMYGOiSwBQlODhYPj4+Sk9Pd2lPT09XWFhYofOPHj2qH374Qf3793e2ORwOSVK1atV06NAhtWjRokzxm2F6Ei8AACi7ip7EW6NGDXXt2lUJCQnONofDoYSEBEVFRRU6/6qrrtLevXu1e/du53H77bfrxhtv1O7duwsNW1UUKjAAAHiY2NhYDRkyRN26dVP37t01d+5cZWdna+jQoZKkwYMHq3HjxoqPj5efn59z37cCderUkaRC7RWJBAYAAAtZsRPvoEGDdPr0aU2dOlVpaWnq3Lmz1q1b55zYm5KSIm/vyj1IU7mju4QZM2bIy8tL48aNszoUAABMs2on3tGjR+v48eO6cOGCtm3bpsjISOdnGzdu1LJly4q9dtmyZVqzZo2Jp3UfWyYwO3bs0Ouvv66OHTtaHQoAALCA7YaQsrKydP/992vx4sV64YUXrA4HQCVw7tw5q0OoVE6ePGl1CLZX8LLiiuBwOJyresrSh6exXQVm1KhR6tevn6Kjo60OBQCAMrNqCMnubFWBWblypZKTk7Vjx44SnX/hwgVduHDB+fPFWywDAAB7sk0FJjU1VWPHjtW7774rPz+/El0THx+voKAg52HVWnUAAIpDBcYc2yQwSUlJOnXqlK6++mpVq1ZN1apV06ZNm/TKK6+oWrVqys/PL3TN5MmTlZGR4TxSU1MtiBwAgOKRwJhjmyGk3r17a+/evS5tQ4cO1VVXXaWJEyfKx8en0DVFvYUTAIDKxIp9YKoC2yQwAQEBhXb8q1WrlurXr2/pToAAAKDi2SaBATxNVlaW1SHYxsXVWU+XkpJidQi29+uvv1bYvajAmGPrBGbjxo1WhwAAQJl5YgJSVraZxAsAAFDA1hUYAADsjiEkc0hgAACwEAmMOQwhAQAA26ECAwCAhajAmEMCAwCAhUhgzCGBASqpc+fOWXbvtWvXWnZvM5o2bWp1CJXK9OnTrQ7B9hwOh9Uh4DJIYAAAsBAVGHNIYAAAsBAJjDkkMAAAWIgExhyWUQMAANuhAgMAgIWowJhDAgMAgIVIYMwhgQEqqb1791p271WrVll2bzMGDRpkdQiVSmpqqtUhAOWOBAYAAAtRgTGHBAYAAAuRwJjDKiQAAGA7VGAAALAQFRhzSGAAALAQCYw5DCEBAADboQIDW8rKyjJ9be3atd0SQ3kvVd2yZUu59n8pR48etezeZvD2ZdgZFRhzSGAAALAQCYw5JDAAAFiIBMYc5sAAAADboQIDAICFqMCYQwIDAIDFPDEBKSuGkAAAgO1QgQEAwEIMIZlDAoPLutSeKyXZU6Use7YU59y5c27vs7QGDx5crv1fd9115dr/pZT3HjcA/svhcMjhcJS5D0/DEBIAALAdKjAAAFiIISRzSGAAALAQCYw5DCEBAADboQIDAICFqMCYQwIDAICFSGDMIYGBU3HLnS+1ZLkky5nXrl1rOqbiNG3a1PS1W7ZscUsMR48edUs/ADwbCYw5zIEBAAC2QwUGAAALUYExhwQGAAALkcCYwxASAACwHSowAABYiAqMOSQwAABYiJc5msMQEgAAsB0qMJXEn/dgqV279iXP3bVrl/PPXbp0KbafAu+8847Lz5GRkaWKbdu2bcV+VpL9WFatWnXZczZu3FiakNSrV69SnW9Hpf2dALAnhpDMIYEBAMBCDCGZwxASAACwHSowAABYiCEkc0hgAACwkGEYZR4CIoEBAAAVigqMOcyBAQAAtkMFxmIFy57PnTtX6LMjR44Uec3Jkyedf754GXVR/Vy8jHnr1q2ljrM4JVkiffToUbfdr0BZlhiHh4e7JYbU1FS39APAs1GBMYcEBgAAC7GM2hyGkAAAgO2QwAAAYKGCIaSyHqW1YMECNW/eXH5+foqMjNT27duLPXfx4sXq2bOn6tatq7p16yo6OvqS51cEEhgAACxUMIRU1qM0Vq1apdjYWE2bNk3Jycnq1KmTYmJidOrUqSLP37hxo+6991598cUXSkxMVHh4uPr06aMTJ06441dgCgkMAAAeZs6cORo5cqSGDh2qtm3batGiRapZs6aWLFlS5PnvvvuuHnvsMXXu3FlXXXWV3nzzTTkcDiUkJFRw5P9FAgMAgIXcOYSUmZnpcly4cKHQ/XJzc5WUlKTo6Ghnm7e3t6Kjo5WYmFiimHNycpSXl6d69eq555dgAquQSikrK0u1a9dWVlaWzp07V+yS3IvfLn3xcum9e/dKkjp06CBJWrt2rfP8gjc8b9myxdn25/v8efnun8+RVKJy3vLly4tsd9fyYgBAyblzGfXF/45PmzZNzzzzjEvbmTNnlJ+fr9DQUJf20NBQHTx4sET3mzhxoho1auSSBFU0EhgAAKqI1NRUBQYGOn/29fV1+z1mzJihlStXauPGjfLz83N7/yVFAgMAgIXcuQ9MYGCgSwJTlODgYPn4+Cg9Pd2lPT09XWFhYZe89uWXX9aMGTP0n//8Rx07dixTzGXFHBgAACxU0cuoa9Sooa5du7pMwC2YkBsVFVXsdbNmzdLzzz+vdevWqVu3bmV6ZnegAgMAgIWs2Ik3NjZWQ4YMUbdu3dS9e3fNnTtX2dnZGjp0qCRp8ODBaty4seLj4yVJM2fO1NSpU7VixQo1b95caWlpkv6Y41m7du0yxW4WCQwAAB5m0KBBOn36tKZOnaq0tDR17txZ69atc07sTUlJkbf3fwdpFi5cqNzcXN11110u/RQ1SbiikMAAAGAhq17mOHr0aI0ePbrIzy5+Ye4PP/xgIqryRQIDAICFeBu1OSQwF/nzPi9F/e+XX34pSfqf//kf7d27V3Xr1tWRI0e0bds29evXT5K0aNGiQmvxC/Z2eemllyRJR48elSRNmTJFkrR161bnuatWrXI552ItWrRw/vniLLkoJd3f5c/7ywAAUJnZZhVSfHy8rrnmGgUEBCgkJEQDBgzQoUOHrA4LAIAyMQyjzO9B8sQKjG0SmE2bNmnUqFHaunWrNmzYoLy8PPXp00fZ2dlWhwYAgGlWvY3a7mwzhLRu3TqXn5ctW6aQkBAlJSXpf/7nfyyKCgAAWME2CczFMjIyJOmSL5K6cOGCy4usMjMzyz0uAABKw+FwyMvLq8x9eBrbDCH9mcPh0Lhx49SjRw+1b9++2PPi4+MVFBTkPHhZIQCgsmEIyRxbJjCjRo3Svn37tHLlykueN3nyZGVkZDgPVtkAAFA12G4IafTo0frkk0/05ZdfqkmTJpc819fXt8g3cWZlZeno0aPq0qWLUlNTVbduXZ07d05r166VJEVGRurkyZPasmWLs2rz52XO0n+XQz/yyCN64IEHdOLECec5y5cvV3h4uDNhKqryU/BZwTX/93//V+w5ZUXiBgCVF0NI5tgmgTEMQ48//rhWr16tjRs3KiIiwuqQAAAoMzayM8c2CcyoUaO0YsUK/fvf/1ZAQIDzRVJBQUHy9/e3ODoAAMwhgTHHNnNgFi5cqIyMDPXq1UsNGzZ0HgW71gIAAM9hmwqMJ2aXAICqjzkw5tgmgQEAoCpiCMkc2wwhAQAAFPDICsyqVavk7++vbdu2aevWrbr22mu1atUq55udw8PDddNNN0mSpk+fLknq1auXy5ufw8PD1aJFC6WmpjrP+fNy6T8vXf7zny9eUn38+PFC51wOy6LLht8fgMqEISRzPDKBAQCgsmAIyRyGkAAAgO1QgQEAwEJUYMwhgQEAwELMgTGHISQAAGA7VGAAALCYJw4BlRUJDAAAFmIOjDlehgc9dWZmpoKCgiQV3tflYuHh4eWyX0jBPjDsRQIAlV9GRoYCAwPLpe+C76TmzZvL27tsMzocDod++OGHco23smEODAAAsB2GkAAAsBBDSOaQwAAAYCF3LIFmGTUAAIANUIEBAMBCDCGZQwIDAICFSGDM8dgE5ujRo5f8vLyWObN8GgCAsvPYBAYAgMqACow5JDAAAFiIBMYcViEBAADboQIDAICFqMCYQwIDAICF2MjOHBIYAAAsRAXGHObAAAAA2/HYCgz7sQAAKgMqMOZ4bAIDAEBlQAJjDkNIAADAdqjAAABgISow5pDAAABgIRIYcxhCAgAAtkMFBgAACxmGUeaN6DyxAkMCAwCAhdyRfHhiAsMQEgAAsB0qMAAAWIgKjDkkMAAAWIgExhwSGAAALEQCYw5zYAAAgO1QgQEAwEJUYMwhgQEAwEIOh0NeXl5l6sMTExiGkAAAgO1QgQEAwEIMIZlDAgMAgIVIYMxhCAkAANgOFRgAACxEBcYcEhgAACxEAmMOQ0gAAMB2qMAAAGAhKjDmkMAAAGAhNrIzhwQGAAALUYExhzkwAAB4oAULFqh58+by8/NTZGSktm/ffsnz//nPf+qqq66Sn5+fOnTooE8//bSCIi0aCQwAABYyDMMtR2msWrVKsbGxmjZtmpKTk9WpUyfFxMTo1KlTRZ7/9ddf695779Xw4cO1a9cuDRgwQAMGDNC+ffvc8SswxcvwoLpTZmamgoKCrA4DAGATGRkZCgwMLJe+y+M7qaTxRkZG6pprrtGrr74q6Y95OOHh4Xr88cc1adKkQucPGjRI2dnZ+uSTT5xt1157rTp37qxFixa57wFKwaMqMB6UqwEA3MBu3xuZmZkux4ULFwqdk5ubq6SkJEVHRzvbvL29FR0drcTExCL7TUxMdDlfkmJiYoo9vyJ4VAJz/vx5q0MAANhIeX5v1KhRQ2FhYW7rr3bt2goPD1dQUJDziI+PL3TemTNnlJ+fr9DQUJf20NBQpaWlFdl3Wlpaqc6vCB61CqlRo0ZKTU1VQEBAmZeslbfMzEyFh4crNTW13MqXlQHPWbXwnFWLJz+nYRg6f/68GjVqVG739fPz07Fjx5Sbm+uW/gzDKPTd5uvr65a+KyOPSmC8vb3VpEkTq8MolcDAwCr9D0cBnrNq4TmrFk99zoqYM+nn5yc/P79yv8+fBQcHy8fHR+np6S7t6enpxVaEwsLCSnV+RfCoISQAADxdjRo11LVrVyUkJDjbHA6HEhISFBUVVeQ1UVFRLudL0oYNG4o9vyJ4VAUGAABIsbGxGjJkiLp166bu3btr7ty5ys7O1tChQyVJgwcPVuPGjZ1zaMaOHasbbrhBs2fPVr9+/bRy5Urt3LlTb7zxhmXPQAJTSfn6+mratGlVevxS4jmrGp6zauE5q65Bgwbp9OnTmjp1qtLS0tS5c2etW7fOOVE3JSVF3t7/HaS57rrrtGLFCj399NOaMmWKWrVqpTVr1qh9+/ZWPYJn7QMDAACqBubAAAAA2yGBAQAAtkMCAwAAbIcEBgAA2A4JTCUSHx+va665RgEBAQoJCdGAAQN06NAhq8MqdzNmzJCXl5fGjRtndShud+LECf31r39V/fr15e/vrw4dOmjnzp1Wh+VW+fn5iouLU0REhPz9/dWiRQs9//zztnuHTFG+/PJL9e/fX40aNZKXl5fWrFnj8rlhGJo6daoaNmwof39/RUdH68iRI9YEWwaXes68vDxNnDhRHTp0UK1atdSoUSMNHjxYP/30k3UBm3S5v88/e+SRR+Tl5aW5c+dWWHwoHRKYSmTTpk0aNWqUtm7dqg0bNigvL099+vRRdna21aGVmx07duj1119Xx44drQ7F7c6dO6cePXqoevXq+uyzz7R//37Nnj1bdevWtTo0t5o5c6YWLlyoV199VQcOHNDMmTM1a9YszZ8/3+rQyiw7O1udOnXSggULivx81qxZeuWVV7Ro0SJt27ZNtWrVUkxMjH777bcKjrRsLvWcOTk5Sk5OVlxcnJKTk/Wvf/1Lhw4d0u23325BpGVzub/PAqtXr9bWrVvL9TUCcAMDldapU6cMScamTZusDqVcnD9/3mjVqpWxYcMG44YbbjDGjh1rdUhuNXHiROP666+3Ooxy169fP2PYsGEubQMHDjTuv/9+iyIqH5KM1atXO392OBxGWFiY8dJLLznbfvnlF8PX19d47733LIjQPS5+zqJs377dkGQcP368YoIqB8U9548//mg0btzY2Ldvn9GsWTPj73//e4XHhpKhAlOJZWRkSJLq1atncSTlY9SoUerXr1+hV7RXFR999JG6deumu+++WyEhIerSpYsWL15sdVhud9111ykhIUGHDx+WJH3zzTfavHmzbr31VosjK1/Hjh1TWlqay3+/QUFBioyMVGJiooWRlb+MjAx5eXmpTp06VofiVg6HQw888IAmTJigdu3aWR0OLoOdeCsph8OhcePGqUePHpbudFheVq5cqeTkZO3YscPqUMrN999/r4ULFyo2NlZTpkzRjh07NGbMGNWoUUNDhgyxOjy3mTRpkjIzM3XVVVfJx8dH+fn5evHFF3X//fdbHVq5SktLkyTnzqUFQkNDnZ9VRb/99psmTpyoe++9t8q94HHmzJmqVq2axowZY3UoKAESmEpq1KhR2rdvnzZv3mx1KG6XmpqqsWPHasOGDRX+FtaK5HA41K1bN02fPl2S1KVLF+3bt0+LFi2qUgnM+++/r3fffVcrVqxQu3bttHv3bo0bN06NGjWqUs+JPyb0/uUvf5FhGFq4cKHV4bhVUlKS5s2bp+TkZHl5eVkdDkqAIaRKaPTo0frkk0/0xRdfqEmTJlaH43ZJSUk6deqUrr76alWrVk3VqlXTpk2b9Morr6hatWrKz8+3OkS3aNiwodq2bevS1qZNG6WkpFgUUfmYMGGCJk2apHvuuUcdOnTQAw88oPHjxztfAldVhYWFSZLS09Nd2tPT052fVSUFycvx48e1YcOGKld9+eqrr3Tq1Ck1bdrU+e/S8ePH9cQTT6h58+ZWh4ciUIGpRAzD0OOPP67Vq1dr48aNioiIsDqkctG7d2/t3bvXpW3o0KG66qqrNHHiRPn4+FgUmXv16NGj0DL4w4cPq1mzZhZFVD5ycnJcXvomST4+PnI4HBZFVDEiIiIUFhamhIQEde7cWZKUmZmpbdu26dFHH7U2ODcrSF6OHDmiL774QvXr17c6JLd74IEHCs3Hi4mJ0QMPPOB8QzMqFxKYSmTUqFFasWKF/v3vfysgIMA5jh4UFCR/f3+Lo3OfgICAQvN6atWqpfr161ep+T7jx4/Xddddp+nTp+svf/mLtm/frjfeeMPS18+Xh/79++vFF19U06ZN1a5dO+3atUtz5szRsGHDrA6tzLKysvTdd985fz527Jh2796tevXqqWnTpho3bpxeeOEFtWrVShEREYqLi1OjRo00YMAA64I24VLP2bBhQ911111KTk7WJ598ovz8fOe/TfXq1VONGjWsCrvULvf3eXFiVr16dYWFhenKK6+s6FBRElYvg8J/SSryWLp0qdWhlbuquIzaMAzj448/Ntq3b2/4+voaV111lfHGG29YHZLbZWZmGmPHjjWaNm1q+Pn5GVdccYXxt7/9zbhw4YLVoZXZF198UeT/J4cMGWIYxh9LqePi4ozQ0FDD19fX6N27t3Ho0CFrgzbhUs957NixYv9t+uKLL6wOvVQu9/d5MZZRV25ehlEFtssEAAAehUm8AADAdkhgAACA7ZDAAAAA2yGBAQAAtkMCAwAAbIcEBgAA2A4JDAAAsB0SGMBGHnzwwTLv8rpx40Z5eXnpl19+ueR5CQkJatOmzSXfTfXMM884t9GvCJMmTdLjjz9eYfcDUHnxKgHARubNm6eK2nvyqaee0tNPP12p3k315JNP6oorrtD48eN1xRVXWB0OAAtRgQFsID8/Xw6HQ0FBQapTp06532/z5s06evSo/vd//7fc71UawcHBiomJ0cKFC60OBYDFSGAAN+vVq5dGjx6t0aNHKygoSMHBwYqLi3OpnFy4cEFPPvmkGjdurFq1aikyMlIbN250fr5s2TLVqVNHH330kdq2bStfX1+lpKQUGkK6cOGCxowZo5CQEPn5+en666/Xjh07XOL59NNP1bp1a/n7++vGG2/UDz/8cNlnWLlypW6++Wb5+fm5tM+YMUOhoaEKCAjQ8OHD9dtvv7l8vmPHDt18880KDg5WUFCQbrjhBiUnJzs/HzZsmG677TaXa/Ly8hQSEqK33npLkvTBBx+oQ4cO8vf3V/369RUdHa3s7Gzn+f3799fKlSsv+wwAqjYSGKAcLF++XNWqVdP27ds1b948zZkzR2+++abz89GjRysxMVErV67Unj17dPfdd+uWW27RkSNHnOfk5ORo5syZevPNN/Xtt98qJCSk0H2eeuopffjhh1q+fLmSk5PVsmVLxcTE6OzZs5Kk1NRUDRw4UP3799fu3bs1YsQITZo06bLxf/XVV+rWrZtL2/vvv69nnnlG06dP186dO9WwYUO99tprLuecP39eQ4YM0ebNm7V161a1atVKffv21fnz5yVJI0aM0Lp163Ty5EnnNZ988olycnI0aNAgnTx5Uvfee6+GDRumAwcOaOPGjRo4cKBL8te9e3f9+OOPJUrEAFRhlr5KEqiCbrjhBqNNmzaGw+Fwtk2cONFo06aNYRiGcfz4ccPHx8c4ceKEy3W9e/c2Jk+ebBiGYSxdutSQZOzevdvlnCFDhhh33HGHYRiGkZWVZVSvXt149913nZ/n5uYajRo1MmbNmmUYhmFMnjzZaNu2rUsfEydONCQZ586dK/YZgoKCjLffftulLSoqynjsscdc2iIjI41OnToV209+fr4REBBgfPzxx862tm3bGjNnznT+3L9/f+PBBx80DMMwkpKSDEnGDz/8UGyfGRkZhiRj48aNxZ4DoOqjAgOUg2uvvVZeXl7On6OionTkyBHl5+dr7969ys/PV+vWrVW7dm3nsWnTJh09etR5TY0aNdSxY8di73H06FHl5eWpR48ezrbq1aure/fuOnDggCTpwIEDioyMdLkuKirqsvH/+uuvhYaPStJXenq6Ro4cqVatWikoKEiBgYHKyspSSkqK85wRI0Zo6dKlzvM/++wzDRs2TJLUqVMn9e7dWx06dNDdd9+txYsX69y5cy738Pf3l/RHhQqA52IVElDBsrKy5OPjo6SkpEIrfGrXru38s7+/v0sSVJGCg4MLJQ4lMWTIEP3888+aN2+emjVrJl9fX0VFRSk3N9d5zuDBgzVp0iQlJibq66+/VkREhHr27ClJ8vHx0YYNG/T1119r/fr1mj9/vv72t79p27ZtioiIkCTn8FiDBg3c8KQA7IoKDFAOtm3b5vJzwXwQHx8fdenSRfn5+Tp16pRatmzpcoSFhZX4Hi1atFCNGjW0ZcsWZ1teXp527Nihtm3bSpLatGmj7du3F4rlcrp06aL9+/e7tLVp06bI5/qzLVu2aMyYMerbt6/atWsnX19fnTlzxuWc+vXra8CAAVq6dKmWLVumoUOHunzu5eWlHj166Nlnn9WuXbtUo0YNrV692vn5vn37VL16dbVr1+6yzwGg6qICA5SDlJQUxcbG6uGHH1ZycrLmz5+v2bNnS5Jat26t+++/X4MHD9bs2bPVpUsXnT59WgkJCerYsaP69etXonvUqlVLjz76qCZMmKB69eqpadOmmjVrlnJycjR8+HBJ0iOPPKLZs2drwoQJGjFihJKSkrRs2bLL9h0TE6Ply5e7tI0dO1YPPvigunXrph49eujdd9/Vt99+67IfS6tWrfSPf/xD3bp1U2ZmpiZMmOAc8vmzESNG6LbbblN+fr6GDBnibN+2bZsSEhLUp08fhYSEaNu2bTp9+rTatGnjPOerr75Sz549i+wXgAexehIOUNXccMMNxmOPPWY88sgjRmBgoFG3bl1jypQpLpN6c3NzjalTpxrNmzc3qlevbjRs2NC48847jT179hiG8cck3qCgoEJ9/3kSr2EYxq+//mo8/vjjRnBwsOHr62v06NHD2L59u8s1H3/8sdGyZUvD19fX6Nmzp7FkyZLLTuL9+eefDT8/P+PgwYMu7S+++KIRHBxs1K5d2xgyZIjx1FNPuUziTU5ONrp162b4+fkZrVq1Mv75z38azZo1M/7+97+79ONwOIxmzZoZffv2dWnfv3+/ERMTYzRo0MDw9fU1WrdubcyfP9/lnCuvvNJ47733io0dgGfwMowK2tYT8BC9evVS586dNXfuXKtDKZMJEyYoMzNTr7/+utv7zsrKUuPGjbV06VINHDiwxNd99tlneuKJJ7Rnzx5Vq0YBGfBkzIEBUKS//e1vatasmRwOh9v6dDgcOnXqlJ5//nnVqVNHt99+e6muz87O1tKlS0leADAHBkDR6tSpoylTpri1z5SUFEVERKhJkyZatmxZqRORu+66y63xALAvhpAAAIDtMIQEAABshwQGAADYDgkMAACwHRIYAABgOyQwAADAdkhgAACA7ZDAAAAA2yGBAQAAtkMCAwAAbOf/ASqw05TAiA3WAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "snr_limit = 6\n", + "\n", + "ax = plt.subplot(111, xlabel=\"period (days)\", ylabel=\"radius ($R_{\\oplus}$)\")\n", + "plt.pcolor(periods, radii, np.mean(snrs.T > snr_limit, 0), cmap=\"Greys_r\")\n", + "_ = plt.colorbar(label=\"recovery\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is a good way to see what can be detected with this optimized GP model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "nuance", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/notebooks/tutorials/exocomet.ipynb b/docs/notebooks/tutorials/exocomet.ipynb index df5a9af..a049139 100644 --- a/docs/notebooks/tutorials/exocomet.ipynb +++ b/docs/notebooks/tutorials/exocomet.ipynb @@ -26,7 +26,6 @@ "metadata": {}, "outputs": [], "source": [ - "# in order to run on all CPUs\n", "import os\n", "import jax\n", "\n", diff --git a/docs/notebooks/tutorials/tess_search.ipynb b/docs/notebooks/tutorials/tess_search.ipynb index 883df44..b886b93 100644 --- a/docs/notebooks/tutorials/tess_search.ipynb +++ b/docs/notebooks/tutorials/tess_search.ipynb @@ -229,7 +229,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -251,15 +251,22 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 8, "metadata": {}, "outputs": [ { - "name": "stderr", - "output_type": "stream", - "text": [ - "100%|██████████| 14370/14370 [00:19<00:00, 720.80it/s]\n" - ] + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "2ec1c2b5f51f4bc2aee273b4ee308dbb", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/14370 [00:00