diff --git a/tests/TFR_tests.ipynb b/tests/TFR_tests.ipynb index 6664dcb..8476acd 100644 --- a/tests/TFR_tests.ipynb +++ b/tests/TFR_tests.ipynb @@ -2,10 +2,29 @@ "cells": [ { "cell_type": "code", - "execution_count": 39, + "execution_count": 40, "id": "f235d0c0-d7fb-46fd-89df-3d66851f457f", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/bartlett/fsigma8/borg_velocity/tests/tfr_inference.py:310: SyntaxWarning: invalid escape sequence '\\s'\n", + " sigma_m = np.median(df['e_Kmag'])\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 40, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", @@ -15,7 +34,12 @@ "import astropy.cosmology\n", "from astropy.coordinates import SkyCoord\n", "import borg_velocity.poisson_process as poisson_process\n", - "from astropy.cosmology import LambdaCDM, z_at_value" + "from astropy.cosmology import LambdaCDM, z_at_value\n", + "import jax.numpy as jnp\n", + "\n", + "import importlib\n", + "import tfr_inference\n", + "importlib.reload(tfr_inference)" ] }, { @@ -327,10 +351,175 @@ "print(L / N)" ] }, + { + "cell_type": "markdown", + "id": "d0e2f630-355c-477c-8b35-ec570eb97382", + "metadata": {}, + "source": [ + "# Check MB" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "id": "f6103c09-b314-4898-b615-33676f7fc66d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[STD ] Building gravity model lpt\n", + "[STD ] | | ___________ \n", + "[STD ] | | /-/_\"/-/_/-/| __________________________ \n", + "[STD ] | | /\"-/-_\"/-_//|| \u001b[34;1mBORG3\u001b[39;0m model\n", + "[STD ] | | /__________/|/| (c) Jens Jasche 2012 - 2019\n", + "[STD ] | | |\"|_'='-]:+|/|| Guilhem Lavaux 2014 - 2019\n", + "[STD ] | | |-+-|.|_'-\"||// __________________________ \n", + "[STD ] | | |[\".[:!+-'=|// \n", + "[STD ] | | |='!+|-:]|-|/ \n", + "[STD ] | | ---------- \n", + "[STD ] | | \n", + "[STD ] | | Please acknowledge the following papers:\n", + "[STD ] | | - Jasche & Lavaux (A&A, 2019, arXiv 1806.11117)\n", + "[STD ] | | - Jasche & Wandelt (MNRAS, 2012, arXiv 1203.3639)\n", + "[STD ] | | - Jasche & Kitaura (MNRAS, 2010, arXiv 0911.2496)\n", + "[STD ] | | - Lavaux, Jasche & Leclercq (arXiV 1909.06396)\n", + "[STD ] | | - And relevant papers depending on the used sub-module/contribution\n", + "[STD ] | | \n", + "\n", + "[STD ] | | This is BORG version 6b1404bfd8011f03ee1b8635275f68d14d16bef3\n", + "[\u001b[35;1mWARNING\u001b[39;0m] | | | Entering [/home/bartlett/borg/libLSS/physics/forwards/primordial_as.cpp]virtual void LibLSS::ForwardPrimordial_As::updateCosmo()\n", + "[\u001b[35;1mWARNING\u001b[39;0m] | | | | Sigma8 is set, but sigma8 normalization is not supported for PRIMORDIAL_AS. Ignoring the supplied sigma8.\n", + "[\u001b[35;1mWARNING\u001b[39;0m] | | | Done (in 2.1512e-05 seconds) (ctx='[/home/bartlett/borg/libLSS/physics/forwards/primordial_as.cpp]virtual void LibLSS::ForwardPrimordial_As::updateCosmo()')\n", + "[STD ] | Smoothing field with rsmooth=4\n", + "[STD ] | | Smoothing field with rsmooth=4\n", + "[STD ] | | Smoothing field with rsmooth=4\n", + "[STD ] | | Smoothing field with rsmooth=4\n", + "[STD ] \tMade 668 of 2000\n", + "[STD ] \tMade 1345 of 2000\n", + "[STD ] \tMade 2000 of 2000\n", + "[STD ] Obtaining peculiar velocities\n", + "[STD ] Radial projection\n", + "[STD ] Making MB data\n", + "26.779654 162.9076 97.677765\n", + "27.107544 164.64821 92.31098\n", + "16.885567 111.88842 61.495857\n", + "0\n", + "0\n", + "Low: 14.289446\n", + "High 16.209558\n", + "0.078125 250.0\n" + ] + }, + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 44, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMgklEQVR4nO3de3hU1b0//veeBEKSJmMukGQkgRDBqISQQopAysULiiJFPBXEC7Xql5aLxgsSsLbSnzIJPUVb8ALUox49XL7PQ0K1goJFAjk5X4iBAEEBjQMEmJyYGieQjAlm1u+POJu57JlkkrnsmXm/nifPafbsmax9zqnzdq3P+ixJCCFAREREpCKaQA+AiIiIyBEDChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREalOZKAH0BsWiwUXLlxAXFwcJEkK9HCIiIioB4QQuHjxInQ6HTQa93MkQRlQLly4gPT09EAPg4iIiHqhvr4egwcPdntPUAaUuLg4AF0PGB8fH+DREBERUU+0tLQgPT1d/h53JygDinVZJz4+ngGFiIgoyPSkPINFskRERKQ6DChERESkOgwoREREpDoeB5R9+/bhrrvugk6ngyRJ2L59u9M9X3zxBWbOnAmtVou4uDjceOONOHv2rPx6e3s7lixZguTkZMTGxmLmzJk4d+5cnx6EiIiIQofHAaW1tRW5ublYt26d4ut1dXUoKChAdnY29u7diyNHjuD555/HgAED5HsKCwtRVlaGLVu2oKKiApcuXcKMGTPQ2dnZ+ychIiKikCEJIUSv3yxJKCsrw6xZs+Rrc+fORb9+/fDuu+8qvsdkMmHgwIF49913MWfOHABX+prs2LEDt912W7d/t6WlBVqtFiaTibt4iIiIgoQn399erUGxWCz48MMPMWLECNx2220YNGgQxo0bZ7cMVF1djcuXL2PatGnyNZ1Oh5EjR6KyslLxc9vb29HS0mL3Q0RERKHLqwGlsbERly5dQnFxMW6//Xbs2rULd999N2bPno3y8nIAQENDA/r374+EhAS796akpKChoUHxc/V6PbRarfzDLrJEREShzeszKADwi1/8Ak8++SRGjx6NoqIizJgxA2+88Ybb9wohXDZuWb58OUwmk/xTX1/vzWETERGRyng1oCQnJyMyMhLXX3+93fXrrrtO3sWTmpqKjo4ONDc3293T2NiIlJQUxc+NioqSu8ayeywREVHo82pA6d+/P/Lz83Hy5Em766dOncKQIUMAAGPGjEG/fv2we/du+XWj0Yja2lpMmDDBm8MhIiKiXjCazKisa4LRZA7YGDw+i+fSpUv46quv5N8NBgNqamqQmJiIjIwMLF26FHPmzMGkSZMwdepUfPTRR/jggw+wd+9eAIBWq8UjjzyCp59+GklJSUhMTMQzzzyDnJwc3HLLLV57MCIiIvLc1qqzWF56DBYBaCRAPzsHc/Iz/D4Oj7cZ7927F1OnTnW6Pn/+fLz99tsAgP/4j/+AXq/HuXPncO2112LlypX4xS9+Id/7/fffY+nSpdi0aRPMZjNuvvlmvPbaaz0ufuU2YyIiIu8zmsyYWLwHFptkECFJqCiaijRtdJ8/35Pv7z71QQkUBhQiIiLvq6xrwryNB5yub37sRozPSurz5wesDwoREREFr8zkWGgcNtRGSBKGJsf4fSwMKERERAQASNNGQz87BxE/tv2IkCSsmj3SK8s7nvK4SJaIiIhC15z8DEwaMRCnm9owNDkmIOEEYEAhIiIiB2na6IAFEysu8RAREZHqMKAQERGR6jCgEBERkeowoBAREZHqMKAQERGR6jCgEBER+ZgaDt9Twxg8wW3GREREPqSGw/fUMAZPcQaFiIjIR4wmsxwMAMAigBWltX6dxVDDGHqDAYWIiMhHDE2tdicDA0CnEDjd1Nbrz/R0qcYXY/AHLvEQERH5iPXwPduA0JfD93qzVOPtMfgLZ1CIiIh8xJuH7/V2qUZNBwB6gjMoREREPuStw/fcLdV095lqOQDQEwwoREREPuaNw/f6ulSjhgMAPcElHiIiIh/yVv+RYF2q6S3OoBAREfmIbVGrBOCxn2fi4YLMXoeKYFyq6S1JCCG6v01dWlpaoNVqYTKZEB8fH+jhEBEROTGazJhYvMepbiRYGqX5giff31ziISIi8gGlolbAefdNsLWg9xcu8RAREfmAUlGrlXX3zb5T39j1NVl2ezZyBmuRmRwb0ss3PcEZFCIiol7obubDWtSq9EWrkYCvvrmIom32fU30O09g3sYDmFi8B1urzvpu8EGANShEREQe8qSjq9FkxlsVp7Fx/9ewfuFKALr78o2QJFQUTQ2pmRTWoBAREXnAkzoQTzu6pmmjkTUoFj/uDgbQfTgBguO8HF9iDQoREYU1T8+38bSjq2Og6algOC/HlziDQkREYas359tYi19tuQsTrnbzWGkkYN19eVh+R3bYNGHrCc6gEBFR2Ort+TaPFmTib/sNsKD7MOFuN4/1vTNydQCAmbk6l03YjCYzDE2tYbPDhwGFiIjClqfn2zh2hr1/XAYW33SN28Bg3c2zfNsxWNC1dLHsjmyMuvoqpyDi6rwcT5ehQgGXeIiIKGx1d76NbfGs43KQAPBfB87ivf8507M/Jl35n1dF98P4rKQezYT0ZhkqFHAGhYiIwpqr820cZy0eLchUXKZ5dW8d4mP6YcGkLMXPdxUwJo0Y2KOA0ttlqGDHGRQiIgp7adpouxkNpVDxt/0GSC7eX7LzhMsZDXcBoyc8LcoNFQwoREREDpRChQXAvHHKdR8WAZeBo68Bo7tlqFDFJR4iIiIHSsWzGgm4cVgiIjUS3nGoO3EXOKwBY0VpLTqF6FXAcLUMFcrY6p6IiEjB1qqzcqiQJMD6bSkBmD4yFR8db4BFXJnR6G5XjdFkDquAocST728GFCIiIheMJjMOnWnGok2H7a5LALYvmoBzzWZYhMDYoYlhGzo84dOzePbt24e77roLOp0OkiRh+/btLu9dsGABJEnCK6+8Yne9vb0dS5YsQXJyMmJjYzFz5kycO3fO06EQERH5jLUx2r9a251eEwC2Vp3Dks2HsWRzDU8f9gGPA0praytyc3Oxbt06t/dt374dBw4cgE6nc3qtsLAQZWVl2LJlCyoqKnDp0iXMmDEDnZ2dng6HiIjI67ZWncXE4j2Yt/EA/vD3zxXv2XzwbNj1JvEnj4tkp0+fjunTp7u95/z581i8eDE+/vhj3HnnnXavmUwmvPnmm3j33Xdxyy23AADee+89pKen45NPPsFtt93m6ZCIiIgUedoe3mgyo/pMM4q2HZNPHFaqg5AUrodDbxJ/8vouHovFggcffBBLly7FDTfc4PR6dXU1Ll++jGnTpsnXdDodRo4cicrKSgYUIiLyCk/bw9ve745GApZNz0bJzhM9bpFPnvN6QCkpKUFkZCQef/xxxdcbGhrQv39/JCQk2F1PSUlBQ0OD4nva29vR3n5lDbClpcV7AyYiopDjafdWx/uVaCTgr3PzMGZoAtK00bgqul+ftg6Te14NKNXV1fjLX/6CQ4cOQZJc9dtTJoRw+R69Xo+VK1d6Y4hERBQGPGkPbzSZ8Y+jF7qdObEIIOknUfL7w7E3iT95tZPs/v370djYiIyMDERGRiIyMhJnzpzB008/jaFDhwIAUlNT0dHRgebmZrv3NjY2IiUlRfFzly9fDpPJJP/U19d7c9hERBRietq91VoM+9KHJ7r9TKX3O7bIJ+/xakB58MEHcfToUdTU1Mg/Op0OS5cuxccffwwAGDNmDPr164fdu3fL7zMajaitrcWECRMUPzcqKgrx8fF2P0RERK6kaaOxbHq2HFKUlmB6sqxjJUngEo6febzEc+nSJXz11Vfy7waDATU1NUhMTERGRgaSkpLs7u/Xrx9SU1Nx7bXXAgC0Wi0eeeQRPP3000hKSkJiYiKeeeYZ5OTkyLt6iIiI+mJr1Vm5iFUC8Ozt1zoVyL5VYehROAEASQCTRgz0/kDJJY9nUD777DPk5eUhLy8PAPDUU08hLy8Pv//973v8GS+//DJmzZqFe++9FxMnTkRMTAw++OADREREeDocIiIiO44zIwLA6o9O2vUoOVLfjA37DT3+TAtcHwZIvuHxDMqUKVPgSXf806dPO10bMGAA1q5di7Vr13r654mIiNzqrkB2a9VZLNt2zKPP5BZi//NqDQoREZG/GE1mVNY1OXVvdVcgazSZUeRhOAGAWXk61p/4GQMKEREFHdtW9I7n4KRpo3F33tV291sDxmenv1XsDOvIsenF9sMX5CDkKhiRd3m9URsREZEvddeEzWgyo+zwebv3lB0+jynXDsTZf3VfR+Kujf2+U9941J2Weo8BhYiIgkp3NSafnf7W6XWLAJZsrunR5yvNsERIEmL6azzqTkt9wyUeIiIKKko1JhoAQ5NjsL68rsdBpDvWL0hrD5XWjk6XwYi8jzMoREQUVNK00dDPznE6cXhF2TF8euKbHn2G0jKOrQhJQunC8WjrsMht7I0mMzQSeECgnzCgEBGRKhlNZhiaWpGZHOu0hOLYNE0APQ4ntuHjv7/6Bq/urYNwCB2rZo9Ebrr9obbWYMQDAv2DAYWIiFRn/b46FO88AeGiGLX6THOPduM40gBy+DCazLj/b/bhBAB+M3mYy8JXHhDoPwwoRESkKuvL66DfeeXwPosAirYdQ8cPFowarEVrRye+bW3v1Wf/dkqWHD6Uim0B4PW9dXhg/BCX4SNNG81g4gcMKERE5Bfulmxs7yne6XyysADw/N+Py787Fsn21BvlX8vhIzM5FpIEpxkUa1t7hpDA4i4eIiLyOXeN1ayMJjP+cfRCj5ZuenrInyPbXTdp2mgUTc92uoeFr+rAgEJERD7lqrGabSdWa4B56UPn2RN3ejORcvT8d/J/XjApC8unZzttKebsSeBxiYeIiHyqu8ZqjgGmpzQAfjs1C6/vrfPovat3nsTM3Ctn6yyYnIWZo3UsfFUZzqAQEZFPuTu8D3BdrNodC7oKWqdcO7Dbe20pNVdL00ZjfFYSw4mKMKAQEZFPpWmjsex218somcmxvVqqAbqWi3ra/8RKAlhjEgQYUIiIyKe2Vp1FyUcnYEFXOJiTP9ip0VpfeDz50ts0RH7FgEJERD7jWF8iAGw6WI8J+is7eQxNrb1qutZbQoDn5wQBBhQiIvIKo8mMyromu905b1UYFOtLBK7s5FGqUfGEBGDR1Kwe389txMGBAYWIiPpMqc+J0WTGhv0Gl++x3cmjn53Tq5UXDYDie3LwwI1DFF+fNy4dy6dnI0Lq+nRuIw4e3GZMRER94qrPyR9mXtftez84eh5Dk2MwJz8D7T904vd//7zHf1cDoGzRBOSmJ6CyrknxnrtGXY3xWUncRhyEGFCIiKhPXPU5+deljm7fu+lAPTYfrMejBZnQRvfz6O9aALR1WABc2cpscTiV2LqUw/Nzgg+XeIiIqE9c9TkZNVjbo/cLAWzcb8C/7zrl0d91DCD62TlcygkhnEEhIqI+sYaDFaW16BQCGgCPFAxFW0enz/6mBnAKIHPyMzBpxEAu5YQISQjHcxzVr6WlBVqtFiaTCfHx8YEeDhERAThS34wN+wzYUWuEEF27a3zxBSMB2P5j7QkFF0++vzmDQkREfba16iyKth2zCyS++rff4ntyGE7CAAMKERH1iXUXjz+m4//OmZOwwSJZIiLqk94c9neDLh5SLxqfWHftUOjjDAoREXnMaDLD0NSKzORY+bA/TzLK8QstADyrU9EAiOnPf68OFyySJSIij2ytOis3ZpMALJyShW/bOrD5YH2vPk+SurYa94RGAvSzczAnP6NXf4sCy5Pvb0ZRIiLqEaPJjH8cvYCibfaH/726t67X4QRwHU40EvDm/DF2S0HWLrW25/1QaOISDxERdWt9eR2Kd57w+6nDn3zR6BRgbM/wodDFgEJERG6t31cH/c4Tfv+7AlCcmeFpxOGBSzxEROTSkfpm6Hd4N5z05YtHqYMshSYGFCIiUrR+Xx1mvVrp1c9cd18eyhZNcDq7p6fWzstjgWyYYEAhIiIn68vroN/h/ZqTT082Ijc9we5gP43UtRuoO5IE/HQIm7SFC9agEBGRnSP1zT6rOdl26DweGj/E6WC/fae+sTtscO64dGw+UG8fkIKuKQb1BQMKERHJ1pf7viD2s9PNyE1PQJo2Wq4lcQwshqZWbDpgXyArAO7eCSMeL/Hs27cPd911F3Q6HSRJwvbt2+XXLl++jGXLliEnJwexsbHQ6XR46KGHcOHCBbvPaG9vx5IlS5CcnIzY2FjMnDkT586d6/PDEBFR7/lrt87YocrLNGnaaIzPSkKaNhqZybFOdSrcvRNePA4ora2tyM3Nxbp165xea2trw6FDh/D888/j0KFDKC0txalTpzBz5ky7+woLC1FWVoYtW7agoqICly5dwowZM9DZ2dn7JyEioh4xmsyorGuya3ZmNJlR7KetxIPiB3R7T5o22q5OJUKSuHsnzPSp1b0kSSgrK8OsWbNc3lNVVYWf/exnOHPmDDIyMmAymTBw4EC8++67mDNnDgDgwoULSE9Px44dO3Dbbbd1+3fZ6p6IqHds29Tbto2vrGvCvI0H/DKGzY/diPFZST2612gyy8s+DCfBT1Wt7k0mEyRJwlVXXQUAqK6uxuXLlzFt2jT5Hp1Oh5EjR6KyUnk7W3t7O1paWux+iIjIM0aTWQ4ngH3beKUlFV/wdJnGdtmHwotPi2S///57FBUVYd68eXJSamhoQP/+/ZGQYL8GmZKSgoaGBsXP0ev1WLlypS+HSkQU8gxNrXI4seoUAh8eNaJTiB4f2Ofozflj8P1lC4QAzpvMWL3zJDqFQIQkYVaeDtsPX5B/5zIN9ZTPAsrly5cxd+5cWCwWvPbaa93eL4SAJCnH9+XLl+Opp56Sf29paUF6errXxkpEFA6OnTcpXn/xwy/69LkHDN9ixR3XA+iapdFpB0AjSfjpkK6dOs/cdi1ON7Uhpr8GrR2dMJrMDCnULZ8ElMuXL+Pee++FwWDAnj177NaZUlNT0dHRgebmZrtZlMbGRkyYMEHx86KiohAVFeWLoRIRhQWjyYxiL7est9q4z4CHJ2Zi36lvFOtb0rTRLl8jcsXrNSjWcPLll1/ik08+QVKSfSHUmDFj0K9fP+zevVu+ZjQaUVtb6zKgEBGRZxx36lSfafZZnzMBoPp0s8v6Fne1L0SueDyDcunSJXz11Vfy7waDATU1NUhMTIROp8O//du/4dChQ/jHP/6Bzs5Oua4kMTER/fv3h1arxSOPPIKnn34aSUlJSExMxDPPPIOcnBzccsst3nsyIqIw5bhTZ9n0bDSYvvfp39z9+f8q1recbmrDpydcv8alHnLF44Dy2WefYerUqfLv1tqQ+fPn44UXXsD7778PABg9erTd+z799FNMmTIFAPDyyy8jMjIS9957L8xmM26++Wa8/fbbiIiI6OVjEBERoLxTx9unESv5+5ELkGDfjT5CkhDTX4ON+w1O92sksOkaueVxQJkyZQrctU7pSVuVAQMGYO3atVi7dq2nf56IiBwYTWYYmlqRmRyruFPHnzRSVyiy7thp7ehUXFp6tGAYZ0/ILZ7FQ0QUxGyXcyQAD40f4jST4U2SBJfbkQWAtXPzkPSTKLmxmtFklkOLlQbAwwVDfTRCChU+b9RGRES+YTSZUbTtynKOAPDO/5zxWTiJkCQUTc+W288rSU+0b6ym1LJef08OZ0+oW5xBISIKUp98/r8+CyNWGgAWXFmymZOfgZm5Onx41KjYP6Wtw+J0zfGkYoYT6gkGFCKiILR+X53Pi19/NX4IFkzJcgoWadpo3DkqDat2fGG3dOOujX2aNprBhDzCJR4ioiCzvtz34QQA6pvNLs/C4WnD5Gt9Os04UHiaMRGFK6PJjAn6PT5f2rFaNCULE4cnIzM5VjF88LRh8oQn399c4iEiUinb7cPWL39DU6vfwgkAvLq3Dq/urXPZnp5LN+QrDChERCrk2A3WGg427vs6IOOxtqefNGIgAwn5BQMKEZHKKHWDXV56DN9casenJ78J2LjYnp78iQGFiEglrEs6/7rU7tQN1iKAf//4VGAG9iN3u3SIvI0BhYhIBbZWnUXRtmMQ6OoI68tusL3BXTrkbwwoREQBZu0Iaw0k1v+phpCiAbB2Xh5+OiSB4YT8igGFiChArEs6dd9cUgwiw1Nicep/W/0+LisJgP6eHNw5ShewMVD4YkAhIgoA2106rgQynADAL0brnLYVE/kLO8kSEfmJ0WTGP45ewH/+j6HbcKIG79dcgNFkDvQwKExxBoWIyA9si2CDhQXgtmIKGM6gEBH5mLWvSaDCidSDezQKN3FbMQUSAwoRkY8ZmloDtpwza7QOUg8SyqMFw1ByDw//I/XgEg8RkY9lJsdCIyEgIeX9Ixe6/bsaAA8XDEWaNhqTRgzk4X+kCpxBISLysTRtNPSzcwLyty3CeYlHwpUlnQhJgv6eHDmMpGmjMT4rieGEAo4zKEREXuR4ArHRZMZnp79FbFQk3pw/Bo+8U+3X8UgAFk7JwhvlX6NTCGgk4JGCTMwYlYa2DgtnSki1JCFEMBWVAwBaWlqg1WphMpkQHx8f6OEQEQFwPoH47ryrse3Q+UAPCxoJWDY9G/+61I6N+wwQsD8hmchfPPn+ZkAhIuol29kSAJhYvEe1vU00EiCEfev8CElCRdFUzqCQ33jy/c0lHiKiHrIu13zXdhlfGFuwpapervG4MydVteEEUC7Q7RSCfU5ItRhQiIh6wF2jNQHgH8ca/D0kj7iaQWGfE1Ir7uIhIuqG42nDambdnSNJV3bvREgS9LNzUMw+JxREOINCRNQNQ1OrR+FEAvweZjToOnnYto8JAKeeJuxzQsGCAYWIqBuZybEehQ5/h5N5P8vAkpuvsetlYuUYQtK00QwmFBS4xENE1I00bTSKpmcHehgu3ZWrY+igkMMZFCKibvzp4xN49dM6AF3LN78YrcMt16XAaDLjpR0n/DYOpVkcFrpSqGJAISJyY+F71dhRe2WHjkDX+TZGkxkHDM1+G8e8celYctNwvH/kAkp2noBFsNCVQhsDChGRC0fqm+3CiZVFwOvhpLsaly0H65E7+CosmJSFmbk6FrpSyGMNChGFNaPJjMq6JhhNZrv/DAAHT3/rt3H8fHiy29ctAlhRWgujycwD/SgscAaFiMLW+vI6FO88AYErPUOs59Qsm56NQ6f9t4Sz/8umbu9h51cKJwwoRBR2jCYz1u75EpsO1MvXbJdXLALQ+7H41fHvu8KCWAonDChEFFZsTxwOJiyIpXDDgEJEYcNoMqs6nGgkyLtznp1+LUZdfRVi+mvQ1mFhQSyFHY+LZPft24e77roLOp0OkiRh+/btdq8LIfDCCy9Ap9MhOjoaU6ZMwfHjx+3uaW9vx5IlS5CcnIzY2FjMnDkT586d69ODEBF1x9DUqtpwIgEoWzgBmx+7ERVFU7FgUhbGZyUhNz2BBbEUljwOKK2trcjNzcW6desUX1+9ejXWrFmDdevWoaqqCqmpqbj11ltx8eJF+Z7CwkKUlZVhy5YtqKiowKVLlzBjxgx0dnb2/kmIiLqRmRwrH6anNrN/ejXDCJENSQjR63+fkCQJZWVlmDVrFoCu2ROdTofCwkIsW7YMQNdsSUpKCkpKSrBgwQKYTCYMHDgQ7777LubMmQMAuHDhAtLT07Fjxw7cdttt3f7dlpYWaLVamEwmxMfH93b4RBRijCYzDE2tyEyOBdA1YxLbPwL1zWYIIZCRGIO3/vs0ttdcCPBInUVIEiqKpjKcUEjz5PvbqzUoBoMBDQ0NmDZtmnwtKioKkydPRmVlJRYsWIDq6mpcvnzZ7h6dToeRI0eisrJSMaC0t7ejvb1d/r2lpcWbwyaiEGBb/Gq7ZThYcAsxkT2vNmpraOjquJiSkmJ3PSUlRX6toaEB/fv3R0JCgst7HOn1emi1WvknPT3dm8MmoiDnWPwqoI5wsnx6Nh77eabia44rTdxCTGTPJ51kJcn+v3pCCKdrjtzds3z5cphMJvmnvr5e8T4iCk9qLH59cdYNWDA5C78uyFQMI0XTsxHx4z/zuIWYyJlXl3hSU1MBdM2SpKWlydcbGxvlWZXU1FR0dHSgubnZbhalsbEREyZMUPzcqKgoREVFeXOoRBRCrMWvagopv9t+HK0dnVgwKQvF9+RgRWktOoWQw8ic/AzMHM0zdYhc8eoMSmZmJlJTU7F79275WkdHB8rLy+XwMWbMGPTr18/uHqPRiNraWpcBhYioO48WZKrucDH9jhNYX16HOfkZqCiaKm8hnjRiICrrulrbc9cOkTKPZ1AuXbqEr776Sv7dYDCgpqYGiYmJyMjIQGFhIVatWoXhw4dj+PDhWLVqFWJiYjBv3jwAgFarxSOPPIKnn34aSUlJSExMxDPPPIOcnBzccsst3nsyIgoLtsWxGgmYl5+BzQfPqqIGBQBKdp7AzNE6pGmjkaaNdhqvfnYO5uRnBHqYRKrjcUD57LPPMHXqVPn3p556CgAwf/58vP3223j22WdhNpuxcOFCNDc3Y9y4cdi1axfi4uLk97z88suIjIzEvffeC7PZjJtvvhlvv/02IiIivPBIRBQuHItjLQKqCicAYAHk3TlK411RWotJIwZyFoXIQZ/6oAQK+6AQEQBU1jVh3sYDgR5Gt/6+aAJy0xNcjnfzYzdifFZSAEZG5F8B64NCROQP1oZssf0jIMF+S7Hj72rQ1mEBoFzMy+3FRMoYUIgoqDg2ZHMMI2oLJ7YBJE0bDf1s5x09XN4hcsaAQkRBw2gyo2jbMTmEqC2M/H3RBJxouOg2gMzJz8CkEQO5vZioGwwoRKR6RpMZn53+FofOfKe6UGIlScCg+AHITU/oNoBYd/QQkWsMKESkalurztrNmqiVEFd26zCAEPUdAwoRqYrjicTLth0L8IiUOda/sNiVyLsYUIhINWxnSyQABcOTAz0klx77+TC8WWFgsSuRjzCgEFFA2W4ZdiyA3f9lUyCH5lKEJOHhgqF4uGAoi12JfIQBhYgCxnbLcLBwnC1hMCHyDQYUIgoIx7bvweD5O6/DHaPSGEqI/EBth38SUZgwNLUGVTiJkCSGEyI/YkAhooDITI6FFOhB9BCLYIn8j0s8ROQztluGrV/utkWxamQ9KydCkvDs9Gsx6uqrWARLFAAMKETkE7YFsBoJ0M/OAQC35+gEWoQkoXTheLR1WBhKiAKMAYWIvM6xANYigOU/biFW2zk61qCkAfDs7dciNz0hwCMiIoA1KETUR0aTGZV1TTCazPI1pQJYC9QTSqwiJAmLpmZBkrrGV/LRCWytOhvoYREROINCRH2gtIwzJz9Dsb5ELUs6GnSFkQhJwrO3X4uSj05A2Mz0rCitxaQRA7m8QxRgDChE1CtKyzgrSmvxXdtl6HeecLr/9pGp+Ki2IeAhpWzRBLnGRGmmp1MI+dA/IgocLvEQUa+4+nJXCicA8HFtA4ruyPbDyNxr67BgfFYS0rTRyEyOhcZhrzMP/SNSBwYUIuoVpS93dywAvrhg8tl4ekIjwS58pGmjoZ+dgwip60HY74RIPbjEQ0S9Yv1yX1Fai04hoAGQm67F4XrXIWR7jdF/A1TwaMEwp/AxJz8Dk0YM5KF/RCrDgEJEvWb9cn+r4jT+VvG123ASaBoADxcMVXwtTRvNYEKkMlziIaI++1vF1wE/V+dXE4a4fX3Z9GyGEKIgwoBCRB5x7HvyVoUh4OFEAvB25Rm394wafJVfxkJE3sElHiLqMdu+JxKA+36Wjk0H6wM9rG63LnNnDlHw4QwKEfWIY98TAaginHSHO3OIghNnUIjILevpw/+61B7wpRxPaACsnZeHnw5JYDghCkIMKETkkuOSTrCwzprcOUoX6KEQUS8xoBCRIqUlnWDw/J3X4Y5RaZw1IQpyrEEhIkVKrezVRILzrE6EJDGcEIUIBhQiUuRpK3t/0kjA9kUTUHwP29QThSou8RCRojRtNG6/IRU7ahsCPRQ71iCSm56A3PQEtqknClEMKERhzrpLJzM51u4L3mgyY+dxdYQTCUDR9GyMGnyVUxBhm3qi0MSAQhTGbHfpaCRAPzsHc/IzAHTVoAgV1KA8cdM1mDsugyGEKMywBoUoTDnu0rEIYEVprdzCPjM5NuBbiyMkieGEKEwxoBCFKaVdOp1C4HRTm/x7oCdQnr39WoYTojDl9YDyww8/4He/+x0yMzMRHR2NYcOG4Y9//CMsFot8jxACL7zwAnQ6HaKjozFlyhQcP37c20MhIjeUdunYnlnzyRf/G4BR2eMBf0Thy+sBpaSkBG+88QbWrVuHL774AqtXr8af/vQnrF27Vr5n9erVWLNmDdatW4eqqiqkpqbi1ltvxcWLF709HCJyYD2NGOiqObFu09VIwG+nDIOhqRVGkxmNLd8HcpiQAB7wRxTGJCG8WwY3Y8YMpKSk4M0335Sv3XPPPYiJicG7774LIQR0Oh0KCwuxbNkyAEB7eztSUlJQUlKCBQsWdPs3WlpaoNVqYTKZEB8f783hE4U0paLY79ouo3jnCaflnGnXp2DX54GbRVl+RzYWTMoK2N8nIu/z5Pvb6zMoBQUF+Oc//4lTp04BAI4cOYKKigrccccdAACDwYCGhgZMmzZNfk9UVBQmT56MyspKxc9sb29HS0uL3Q8R9Yx1xuRIfbNTUezy0mPQK4QTAD4PJ5KbCtzl0xlOiMKd17cZL1u2DCaTCdnZ2YiIiEBnZydeeukl3HfffQCAhoauvgopKSl270tJScGZM2cUP1Ov12PlypXeHipRyHM87M8xiASylb2rudt19+VhRi4P+SMKd16fQdm6dSvee+89bNq0CYcOHcI777yDf//3f8c777xjd5/k8K9PQgina1bLly+HyWSSf+rr6709bKKQc6S+GUXbguuwvwhJwpihCYEeBhGpgNdnUJYuXYqioiLMnTsXAJCTk4MzZ85Ar9dj/vz5SE1NBdA1k5KWlia/r7Gx0WlWxSoqKgpRUVHeHipRyNpadRbLth1TfE0jBXbmxJEGgAU8S4eI7Hk9oLS1tUGjsZ+YiYiIkLcZZ2ZmIjU1Fbt370ZeXh4AoKOjA+Xl5SgpKfH2cIhCkmN7etvfAaDIVTgBcF1qPI4b/V/HpbTEFCFJKF04Hm0dFp6lQ0R2vB5Q7rrrLrz00kvIyMjADTfcgMOHD2PNmjX49a9/DaBraaewsBCrVq3C8OHDMXz4cKxatQoxMTGYN2+et4dDFHIcd+LcnXc1yg6fl3+fm5/ucjnHAgQknADA2vvy0NrxA1aU1qJTCLtD/4iIHHl9m/HFixfx/PPPo6ysDI2NjdDpdLjvvvvw+9//Hv379wfQVW+ycuVKrF+/Hs3NzRg3bhxeffVVjBw5skd/g9uMKVwZTWZMLN6jqiWanoiQJFQUTZVne3j6MFF48uT72+sBxR8YUCjcWJdwvm3twOJNhwM9nB6x1rpYZ0qshxASUfjy5PubpxkTqZzjVmGlWo5Am3LtQOw7+Y1dseukEQM5U0JEvcaAQqRijicOCyDgJww7ktDVkRaAUyBhMCGi3mJAIVIxpROH1TZ7ItAVTMZnJTGQEJHXMKAQqZC15sTc8UOgh9It2xOQiYi8hQGFSGXsak58uJ7jjVoWNlcjIl9hQCFSEaeaEx+u51yvi8Pgq2LwcS8PBdRIQOnC8exjQkQ+wYBCpCJKNSe+cvzCRRy/cLHX77cIoK3D4sURERFdwYBCpALBVHNixdoTIvIlBhSiAHPsc6Jm1roV1p4Qka8xoBD5kdIhf459TtRs7X15SPpJFJuvEZHPMaAQ+YhjGHE85E8/OwfpiTGqPFfH1cnDY4YmMJgQkV8woBD5gGMYWXZ7Nko+OiGHEYsAVpTWonTheEiSb3fr9IZtx1ou6RBRIDCgEHmZ47KNRQAlO0/Acb9LpxCoOt2MRVOysO7TOr+PszsCwKvz8pAYyyUdIvI/BhQiL1PaKmwBFGdKXvzwC2gkYHT6Vaip/85fQ+yRCEnCT4dwSYeIAkMT6AEQhZrM5Fin3TiSBEwfmap4v0VAleGESzpEFEicQSHyAyGAj2obAj0Mt16cdQNyrtaircPCJR0iCjgGFCIvMzS1Km4XVuNuHat5P0vHzdelMJQQkWpwiYfIy5SWeNRu08F6TCzeg61VZwM9FCIiAAwoRH1mNJlRWdcEo8kc6KF4xDFEWbc+B9tzEFFo4hIPUR+4ar6m4tUczPtZOu7KvRr/am3H4k2H7V7rFAKnm9q41ENEAccZFCIXupsZUep3sqK0FuaOHyCpdI1HIwFLbh6O8VlJGDMkARqHcfIAQCJSCwYUIgVbq85iYvEezNt4wGVthlK/k04h8Og71arrDAt0bXXWz86RZ0fStNHQz85BxI9piluLiUhNuMRD5MDVzMikEQPtvrxj+0covl+F2QQAsH3hBOSmJ9hdm5OfgUkjBuJ0Uxu3FhORqjCgEDlwNTPiWJvR2tHp55H1TVuHY7P9LmnaaAYTIlIdBhQiB5nJsdBI9n1LrLUZticUK92nVqwtIaJgw4BC5MBam7GitBadQsi1GftOfSMv/UgScMfIVMzM1WF7zYVAD9mOBGD2T6/G9sMX7MbPWRIiCiaSEGos53OvpaUFWq0WJpMJ8fHxgR4OhRjrLEls/wi57TsATCzeo/rZkv/z82F4uGAo0rTRMJrMrC0hIlXx5PubMyhENpT6mozPSkJlXZPqw4kEyOEEYG0JEQU3bjMm+pHS7p3lpcdwpL4ZmcmxgR3cjyQAb84fo9hnRQA43dTm7yEREfkEAwrRj96qMDjNklgEMOvVSvx/H3wemEE5KJqejZuvS0XR9Gyn1yQJLIQlopDBJR4idM2ebNhvUHxNANhR2+DfATmQABTdkY0Fk7Jc36TyJSgiIk8woFBYst0unKaNxmenvw30kNzavuhKkzWjyYzinSec7rEu8bDuhIhCAQMKhR2lQtjYKPX+V0ECcKLhohxQDE2tiq30NeASDxGFDtagUFhRKoQt2nZM1TMoAl2t9q2HFrpqsf/bKVmcPSGikMGAQiHJ1UnESm3sBYC3K8/4b3C9YG21D7husV8wfKA/h0RE5FPqndcm6iWlJZw5+RkAXM8+qJ3GZoeOu1b8REShgjMoFFJcnURsnUk5+21w9gnJH5po14BNPzsHET82Q2EreyIKRT4JKOfPn8cDDzyApKQkxMTEYPTo0aiurpZfF0LghRdegE6nQ3R0NKZMmYLjx4/7YigUZtydRAwA/1P3rwCMqnsSgFmjdS5frzr9rd1y1Zz8DFQUTcXmx25ERdFUeYaIiChUeD2gNDc3Y+LEiejXrx927tyJzz//HH/+859x1VVXyfesXr0aa9aswbp161BVVYXU1FTceuutuHjxoreHQ2HGuvxhy/Yk4k0H6wMzsG4snJqFe/PTXb5uEc5dYtO00RiflcSZEyIKSV6vQSkpKUF6ejreeust+drQoUPl/yyEwCuvvILnnnsOs2fPBgC88847SElJwaZNm7BgwQJvD4nCiHX5w7YGxbr8UbjlUKCH59Lre+sQH93PqbbEijUmRBRuvD6D8v7772Ps2LH45S9/iUGDBiEvLw8bN26UXzcYDGhoaMC0adPka1FRUZg8eTIqKysVP7O9vR0tLS12P0TuWPuEWP/n+n112F5jDNyAHEg//lhZBLB650ksuz1bri2xYo0JEYUjr8+gfP3113j99dfx1FNPYcWKFTh48CAef/xxREVF4aGHHkJDQ1fL8JSUFLv3paSk4MwZ5a2eer0eK1eu9PZQKQQZTWYUbTsmd30XAJZtOxbIISkqGJ6M/V822V3rFAKjBl+FiqKpON3Uhpj+GrR1WDA0OYbhhIjCjtcDisViwdixY7Fq1SoAQF5eHo4fP47XX38dDz30kHyf5PBviUIIp2tWy5cvx1NPPSX/3tLSgvR01+v1FL7+o8IQFEfSVHzZ5HKrcJo2moGEiMKe15d40tLScP3119tdu+6663D27FkAQGpqKgDIMylWjY2NTrMqVlFRUYiPj7f7IXJkNJnxNxcH/qmNAPBowTBuFSYicsHrMygTJ07EyZMn7a6dOnUKQ4YMAQBkZmYiNTUVu3fvRl5eHgCgo6MD5eXlKCkp8fZwKAxYD/6ra7wUFLMnQNe/GTxcMBQPFwzF6aY2LuMQETnwekB58sknMWHCBKxatQr33nsvDh48iA0bNmDDhg0AupZ2CgsLsWrVKgwfPhzDhw/HqlWrEBMTg3nz5nl7OBTibLvGBgsJgP6eHLvGa0REZM/rASU/Px9lZWVYvnw5/vjHPyIzMxOvvPIK7r//fvmeZ599FmazGQsXLkRzczPGjRuHXbt2IS4uztvDoRDmWBAbDDQSULZwgnwyMRERKZOEUDq4Xd1aWlqg1WphMplYjxLGVpQeVW3jNXc2P3YjxmclBXoYRER+58n3N8/ioaC0vrxOleFEkq70N9HAvtcJwIZrREQ9xdOMKegcqW+GfueJQA/DiXX5ZlD8ALnwdd+pb7CitBadQnCnDhGRBxhQSHWsu3Iyk2Odvsy3Vp1FkQobrwFdPU3aOix2fUzm5Gdg0oiB3KlDROQhBhRSFWsAEehaHim6Ixs5V2uRmRyLxpbvUVSq3qJYV8s3bLxGROQ5BhRSDaU29fodXUs50o+/q5XtoYRERNR3DCikGp+d/tZlCFFzOAGAv87Nw4xcXaCHQUQUMriLh1TD1VlMaqeRgDFD2deEiMibGFDI74wmMyrrmmA0me2ujxmS4LQtV+0kCdDPzuHSDhGRl3GJh/zKtjW9BKBoejYWTM4C0FVMWnxPjt3rgP+Wd+4YmYqdxxvgrnVhhCShdOF41H9rhiQBPx2SwHBCROQDDCjkN0aT2e7cHAFAv/MEznzbiiU3DUeaNhrftV22O1en6I5s/OtiBzbs/9rn49tZ2+A2DFn7mOSmJ7BVPRGRjzGgkFe562FiaGpVPNRv04F6bDlYj9tHpmLHsQb5ugBQsuMEhJ/WfZTCyfN3XoexQxPQ1mFhHxMiIj9iQCGvsVu+kX5cvpmUJYeW881tLt9rEbALJ/J1wKtrPOMyE3DTdSko3nnC7VIO0DVjcseoNIYSIqIAYEAhr3BavhFdPUxqzn6Hj483KM6cBMJnp7/DK3PzcGNmIma9WmmXfSR0BSuLANvSExEFGAMKeYWr5Zudtc6zIoHUKQRON7VhfFYSiu/JcTonh23piYjUgQGFvCIzORaShG6XTQLNth29q3NyGEyIiAKPfVDIK9K00Sianh3oYbil1I4+TRuN8VlJDCVERCrDgEJeMzNXh5+ptKOqBKBs4QTMyc8I9FCIiKgHuMRDXrF+Xx2Kd5xQ5Zk5kgQUz85h7xIioiDCgEJ9tr68DvqdJwI9DEVP3HwN5v4sg0s4RERBhgGFes1oMuOz09+iWKXhJEKSGE6IiIIUa1CoRxwP+NtadRYTi/dgyeaabpd1NABenHWDXw8C1MC5IJaIiIIHZ1CoW7YdYjUSsOz2bJR8dKLHzdcsAL6/bPFbfYokdRXEsuaEiCh4cQaF3HLsEGsRQMlO1+FEaZYkQpKQPzQBGj9NoRRNz2Y4ISIKcgwo5JZSh1hLN++556dXI0LqSiO2JwDrZ+f4NKRoACz/8fwfIiIKblziIbcyk2Oh+fF8GisNgFuvT8HHn/+v0/0CwPbDF1C6cLziCcDePpNn0vBkPD1tBE8bJiIKMZxBIZljISzQ1WlVPztHnhEBumZQdimEE6tOIdDWYbHr0PrPLxqwbNsxr4/5v7/6FwbFD2A3WCKiEMMZFALgXAirn50jd12dk5+B7NQ4zHqtUj5rx91EiEaCfN6N0WRG4ZbDOGBo9sm4rYf/MZwQEYUWzqCQYiHsitJau5mU1o7OHh8EODe/q/fI1qqzmKDf47VwMi4z0akI1/bwPyIiCh0MKKRYCGudmbCy1qLYclXwuuTma+TQ05eSEwnAm/PHYMaoVADAAcO38nXgSgEuZ0+IiEIPl3hIuRBWAmL6X8mv1lqUFaW16BRCDgcAsHzbMXlnz6KpWWhs+R7/97P6PhfECgDJP4nCjmMNdtc0ErB2bh7GDE1gOCEiClEMKOQUPoCusHL3a5WKtSj//KIRA+OiMGnEQKRpo1F+8hvsqO0KEa9+WodXP63z2tj2nGh03uYsgKSfRDGcEBGFMAYUAtAVPuq/bcM6m3BhrUXJTo1Da0cnjp0z2R0K+Pzfj2Px1Cw5nPhC8k+inGZ3WHdCRBT6GFAIQFeh7Kt7nWc+OoXArFcrXdaSrPPibIkjCcAt16egf6TGaWmJsydERKGNASWEGU1mGJpakZkc2+0XuqGp1eUuHV+eoSNJUPy7GgD6e3KQpo3GnPwMTBoxEKeb2tiMjYgoTDCghCh3fU0A5/CiVCjrSxKAx36eiQ37DU6vPX/ndbhjVJpdEEnTRjOYEBGFEQaUINOTWRFXfU2sRa224UVC1+F6CyZnORXKepOErh0+EZKEgfFRuPm6FADA3yoMTvUljuGEiIjCj8/7oOj1ekiShMLCQvmaEAIvvPACdDodoqOjMWXKFBw/ftzXQwl6W6vOYmLxHszbeAATi/dga9VZxfvc9TVxDC8CgH7nCazfV4dJIwbilbm5eHHWDZC8fKifJAH33zgET067Fg/cOFSeEbFto8/6EiIisvLpDEpVVRU2bNiAUaNG2V1fvXo11qxZg7fffhsjRozAiy++iFtvvRUnT55EXFycL4cUtLqbFbF17LzJ6f0SutrPK4UXACjecQLFO07IfUZ+fk0y9n3Z5LXxWwQUW9KzvoSIiJT4bAbl0qVLuP/++7Fx40YkJCTI14UQeOWVV/Dcc89h9uzZGDlyJN555x20tbVh06ZNvhpO0OtJt1egK8gU7zgBVzKTY53axQNdMynWj7cIeDWcAO63Bqdpo3nYHxER2fFZQFm0aBHuvPNO3HLLLXbXDQYDGhoaMG3aNPlaVFQUJk+ejMrKSl8NJ+gptZpX+tKvPtOsuOtGANh8oGtJqGh6tm8G6YJGApduiIjIIz4JKFu2bMGhQ4eg1+udXmto6GrqlZKSYnc9JSVFfs1Re3s7Wlpa7H7CTU/qNbZWncXiTYddfsZf93yF8fo9+Neldvxq/BCfj1n+u3Pz7HYQERERdcfrNSj19fV44oknsGvXLgwYMMDlfZJDFaYQwumalV6vx8qVK706zmDkrl7DaDKjaNuxHn3Ohv0GxWUeX9AAGDM0odv7iIiIbHl9BqW6uhqNjY0YM2YMIiMjERkZifLycvz1r39FZGSkPHPiOFvS2NjoNKtitXz5cphMJvmnvr7e28MOGq7qNVwt7bjip3YnWDY9m0s7RETkMa8HlJtvvhnHjh1DTU2N/DN27Fjcf//9qKmpwbBhw5Camordu3fL7+no6EB5eTkmTJig+JlRUVGIj4+3+6ErjCYzTjYEftnrjpxUeQlKIwHL7+jqr0JEROQpry/xxMXFYeTIkXbXYmNjkZSUJF8vLCzEqlWrMHz4cAwfPhyrVq1CTEwM5s2b5+3hhLytVWdRtO2Yz2dEJLieddGga6ZkweQsGE1mbhkmIqI+C0gn2WeffRZmsxkLFy5Ec3Mzxo0bh127drEHioesdSeehBNXZ9+4owEwd1w6Nh1wXlp7/KZrcN+4DDmMsCU9ERF5gySED/qa+1hLSwu0Wi1MJlNYL/d8cOQ8lmyucbo+O0+H+ROG4kTDRadTgLNT45xOJ9YAmJ6Tig+POe+iun9cBhbfdA0AYIJ+j937JAmoLLqJgYSIiHrEk+9vn7e6J99xteup9PAF3P1aV0+ZiqKp2PzYjagomoo5+RnITU/AwqlZ8i6eCEnCsunZ2FnrHE4WTcnCS3fnyLMixffkyL1YNBJQPDuH4YSIiHyChwWqnLvDAccMcb1919oKv6JoKsZnJcnXn/6/Ndh26Lz8+203pCBnsFax/f3lTvuLbEtPRET+whkUFXN3OKDRZEb1mWa373dshX+kvtkunADAjtoGmDt+UOyL8reKr2E0me2usS09ERH5AwOKSrk6HNBoMsvBxV3XWKBr+SamvwaVdU0wmsw4ePpbxftON7XhsZ9nOl23HvBHRETkb1ziUSlXhwNWn262Cy6uSAAmXpOEu1+rhEV01YwsnKLck2Ts0AQMih+Av1UY7D7X3QF/REREvsQZFJVSOnU4QpIACW7DiQRg1NVaCHSdSGw7A/P63q8xfWSq3f33/PRq5KYn9OisHyIiIn/hDEqAuCt+BYB9p76x+10C8Oz0a5GeEO22aZoAcPS8SfG1TiHw0Pih+M3kYfjsdDPGDk1AbvqVQlsWwRIRkVowoATA1qqz8jKNRgL0s3PsTvs9Ut+MolL7BmwCgH7HiT79XeuSTZo22i6Y2GKjNSIiUgMu8fiZu+JXoCu8zHq10uNur93RSOCSDRERBQ3OoPiZq+JX626Z5aXePVdHAvDYpEw8PDGT4YSIiIIGA4qfZSbHQuNQ6GpdelEKL56YNy4dWw7Wy0tHjxYMw8MFQxlMiIgo6HCJx8/c7ZaxhhdHyg3tHe6RgM0HusKJBGDZ7dlYced1DCdERBSUeFhggBhNZsXdMlurzsoH/GkAPPrzTIwblohH3ql2+VmOMzJAV/CpKJrKgEJERKrhyfc3A4oKGU1mvFVxGn+r+LprRkRCr4pmNz92o905PERERIHE04xDgDWcAL0LJ+wCS0REwYwBxY+MJrN8Lo47fS2WZRdYIiIKdtzF4yfdNWezpbTTpyc0EvDXuXkYMzSB4YSIiIIaZ1D8oLvmbI7StNGY+zPl8OKKBsCjBZkMJ0REFBIYUPygu+ZsSsYPS3T5mgbApOHJ8pZk69k8G/YbMLF4D7ZWne3zmImIiAKJAcUPlPqbdFfEesH0vdM1CcCr8/JQtmgCKr66clKxwJXDA7ubnSEiIgoGDCheplQI6645m6vPKNnpfDBg0fRs3DlKh9aOTrf1Kd3NzhAREakdi2QdGE1mGJpakZkc63Etx/ryOhTvPAEB50LYOfkZmDRioGJzNkeudvGMGnwVgO6LaLnFmIiIgh0Dig1Pdto4Wr+vDnqbWQ/rUkt2ahxaOzrlwGMbTFyFIaUAogEQ079rwss6I2PtOCv9WIQiwC3GREQUGthJ9kdGkxkTi/c4HeLXk3bxRpMZE4r3KDZUsxawaiRg2fRs5FytRWZyLPad+gZF27pOLpYAFN9jH4ZsW95bOYYm23b5AHo0O0NERBQonnx/cwblR+522nT3hW9oanXZ7dW2eFW/o2uGRbK5br2naNsxTBoxUP5bc/IzkJ0ah1mvVcqfbZ2Vsd7nOCPDYEJERKGCRbI/6s1OG3fvdUcpywgA1aeb7a61dnQ6BR8WwBIRUThgQPmR404bjQQ8O/3aHs1KKL138dQsj0IL0HUooK2+hCYiIqJgxiUeG3PyM/Bd22UU7zwBiwBKdp7AVdH9elQoq7RLJz0xxqmOxBVJAn46JMHummMxLAtgiYgoXLBI1kZfCmXdfebppjYcPfcdVn90Ug4as/J0KDt8vkc7hmyLYRlOiIgoWLFItpf6UijrirWQdXxWEmaO1tkFjWduuxanm9oQ01+D1o5OGE1mxb/jWAxLREQU6hhQbCj1H/FmzYfSrpt9p77pde8VIiKiUMUiWRuetqTvKaX299brnpxyTEREFC44g+LAk5b03TGazHirwoCN+w2K7e99saREREQUChhQFHij5sO2bb6VY6M1Xy8pERERBSsu8fiA49KNLdtGa75aUiIiIgp2nEHxAVenEQPOMyTeXFIiIiIKFQwoPqC0dAN0TVcpzZBwGzEREZE9ry/x6PV65OfnIy4uDoMGDcKsWbNw8uRJu3uEEHjhhReg0+kQHR2NKVOm4Pjx494eSsAotb7/Pz8fhv9efhO3EBMREfWA12dQysvLsWjRIuTn5+OHH37Ac889h2nTpuHzzz9HbGwsAGD16tVYs2YN3n77bYwYMQIvvvgibr31Vpw8eRJxcXHeHlJAcOmGiIio93ze6v6bb77BoEGDUF5ejkmTJkEIAZ1Oh8LCQixbtgwA0N7ejpSUFJSUlGDBggXdfqavWt0TERGR73jy/e3zXTwmkwkAkJiYCAAwGAxoaGjAtGnT5HuioqIwefJkVFZWKn5Ge3s7Wlpa7H58zVVzNSIiIvI9nxbJCiHw1FNPoaCgACNHjgQANDQ0AABSUlLs7k1JScGZM2cUP0ev12PlypW+HKod2x4mbD9PRETkfz6dQVm8eDGOHj2KzZs3O70m/VhAaiWEcLpmtXz5cphMJvmnvr7eJ+MF2H6eiIhIDXw2g7JkyRK8//772LdvHwYPHixfT01NBdA1k5KWliZfb2xsdJpVsYqKikJUVJSvhmqH7eeJiIgCz+szKEIILF68GKWlpdizZw8yMzPtXs/MzERqaip2794tX+vo6EB5eTkmTJjg7eF4zNrDxBbbzxMREfmX1wPKokWL8N5772HTpk2Ii4tDQ0MDGhoaYDZ3LZFIkoTCwkKsWrUKZWVlqK2txa9+9SvExMRg3rx53h6Ox9h+noiIKPC8vs3YVR3JW2+9hV/96lcAumZZVq5cifXr16O5uRnjxo3Dq6++KhfSdscf24yNJjN7mBAREXmRJ9/fPu+D4gv+7INiNJlhaGpFZnIsgwoREVEfePL9zbN43OB2YyIiosDweaO2YMXtxkRERIHDgOKCu+3GRERE5FsMKC5wuzEREVHgMKC4wO3GREREgcMiWTfm5Gdg0oiB3G5MRETkZwwo3UjTRjOYEBER+RmXeIiIiEh1GFCIiIhIdRhQiIiISHUYUIiIiEh1GFCIiIhIdRhQiIiISHUYUIiIiEh1GFCIiIhIdRhQiIiISHUYUIiIiEh1GFCIiIhIdYLyLB4hBACgpaUlwCMhIiKinrJ+b1u/x90JyoBy8eJFAEB6enqAR0JERESeunjxIrRardt7JNGTGKMyFosFFy5cQFxcHCRJ8uvfbmlpQXp6Ourr6xEfH+/Xv60GfH4+P5+fzx+Ozx/Ozw547/mFELh48SJ0Oh00GvdVJkE5g6LRaDB48OCAjiE+Pj4s/5/Uis/P5+fz8/nDUTg/O+Cd5+9u5sSKRbJERESkOgwoREREpDoMKB6KiorCH/7wB0RFRQV6KAHB5+fz8/n5/OH4/OH87EBgnj8oi2SJiIgotHEGhYiIiFSHAYWIiIhUhwGFiIiIVIcBhYiIiFSHAcUFvV6P/Px8xMXFYdCgQZg1axZOnjxpd48QAi+88AJ0Oh2io6MxZcoUHD9+PEAj9h29Xg9JklBYWChfC/VnP3/+PB544AEkJSUhJiYGo0ePRnV1tfx6KD//Dz/8gN/97nfIzMxEdHQ0hg0bhj/+8Y+wWCzyPaH0/Pv27cNdd90FnU4HSZKwfft2u9d78qzt7e1YsmQJkpOTERsbi5kzZ+LcuXN+fIrec/f8ly9fxrJly5CTk4PY2FjodDo89NBDuHDhgt1nhOrzO1qwYAEkScIrr7xidz3Un/+LL77AzJkzodVqERcXhxtvvBFnz56VX/fV8zOguFBeXo5Fixbh//2//4fdu3fjhx9+wLRp09Da2irfs3r1aqxZswbr1q1DVVUVUlNTceutt8pnBYWCqqoqbNiwAaNGjbK7HsrP3tzcjIkTJ6Jfv37YuXMnPv/8c/z5z3/GVVddJd8Tys9fUlKCN954A+vWrcMXX3yB1atX409/+hPWrl0r3xNKz9/a2orc3FysW7dO8fWePGthYSHKysqwZcsWVFRU4NKlS5gxYwY6Ozv99Ri95u7529racOjQITz//PM4dOgQSktLcerUKcycOdPuvlB9flvbt2/HgQMHoNPpnF4L5eevq6tDQUEBsrOzsXfvXhw5cgTPP/88BgwYIN/js+cX1CONjY0CgCgvLxdCCGGxWERqaqooLi6W7/n++++FVqsVb7zxRqCG6VUXL14Uw4cPF7t37xaTJ08WTzzxhBAi9J992bJloqCgwOXrof78d955p/j1r39td2327NnigQceEEKE9vMDEGVlZfLvPXnW7777TvTr109s2bJFvuf8+fNCo9GIjz76yG9j9wbH51dy8OBBAUCcOXNGCBEez3/u3Dlx9dVXi9raWjFkyBDx8ssvy6+F+vPPmTNH/u++El8+P2dQeshkMgEAEhMTAQAGgwENDQ2YNm2afE9UVBQmT56MysrKgIzR2xYtWoQ777wTt9xyi931UH/2999/H2PHjsUvf/lLDBo0CHl5edi4caP8eqg/f0FBAf75z3/i1KlTAIAjR46goqICd9xxB4DQf35bPXnW6upqXL582e4enU6HkSNHhtz/PoCufxZKkiTPKIb681ssFjz44INYunQpbrjhBqfXQ/n5LRYLPvzwQ4wYMQK33XYbBg0ahHHjxtktA/ny+RlQekAIgaeeegoFBQUYOXIkAKChoQEAkJKSYndvSkqK/Fow27JlCw4dOgS9Xu/0Wqg/+9dff43XX38dw4cPx8cff4zf/OY3ePzxx/Gf//mfAEL/+ZctW4b77rsP2dnZ6NevH/Ly8lBYWIj77rsPQOg/v62ePGtDQwP69++PhIQEl/eEiu+//x5FRUWYN2+efGBcqD9/SUkJIiMj8fjjjyu+HsrP39jYiEuXLqG4uBi33347du3ahbvvvhuzZ89GeXk5AN8+f1CeZuxvixcvxtGjR1FRUeH0miRJdr8LIZyuBZv6+no88cQT2LVrl906o6NQfHag698axo4di1WrVgEA8vLycPz4cbz++ut46KGH5PtC9fm3bt2K9957D5s2bcINN9yAmpoaFBYWQqfTYf78+fJ9ofr8SnrzrKH2v4/Lly9j7ty5sFgseO2117q9PxSev7q6Gn/5y19w6NAhj58lFJ7fWhj/i1/8Ak8++SQAYPTo0aisrMQbb7yByZMnu3yvN56fMyjdWLJkCd5//318+umnGDx4sHw9NTUVAJwSYmNjo9O/bQWb6upqNDY2YsyYMYiMjERkZCTKy8vx17/+FZGRkfLzheKzA0BaWhquv/56u2vXXXedXLUeyv+3B4ClS5eiqKgIc+fORU5ODh588EE8+eST8mxaqD+/rZ48a2pqKjo6OtDc3OzynmB3+fJl3HvvvTAYDNi9e7c8ewKE9vPv378fjY2NyMjIkP9ZeObMGTz99NMYOnQogNB+/uTkZERGRnb7z0NfPT8DigtCCCxevBilpaXYs2cPMjMz7V7PzMxEamoqdu/eLV/r6OhAeXk5JkyY4O/hetXNN9+MY8eOoaamRv4ZO3Ys7r//ftTU1GDYsGEh++wAMHHiRKct5adOncKQIUMAhPb/7YGunRsajf0/GiIiIuR/mwr157fVk2cdM2YM+vXrZ3eP0WhEbW1tSPzvwxpOvvzyS3zyySdISkqyez2Un//BBx/E0aNH7f5ZqNPpsHTpUnz88ccAQvv5+/fvj/z8fLf/PPTp8/epxDaE/fa3vxVarVbs3btXGI1G+aetrU2+p7i4WGi1WlFaWiqOHTsm7rvvPpGWliZaWloCOHLfsN3FI0RoP/vBgwdFZGSkeOmll8SXX34p/uu//kvExMSI9957T74nlJ9//vz54uqrrxb/+Mc/hMFgEKWlpSI5OVk8++yz8j2h9PwXL14Uhw8fFocPHxYAxJo1a8Thw4flXSo9edbf/OY3YvDgweKTTz4Rhw4dEjfddJPIzc0VP/zwQ6Aeq8fcPf/ly5fFzJkzxeDBg0VNTY3dPwvb29vlzwjV51fiuItHiNB+/tLSUtGvXz+xYcMG8eWXX4q1a9eKiIgIsX//fvkzfPX8DCguAFD8eeutt+R7LBaL+MMf/iBSU1NFVFSUmDRpkjh27FjgBu1DjgEl1J/9gw8+ECNHjhRRUVEiOztbbNiwwe71UH7+lpYW8cQTT4iMjAwxYMAAMWzYMPHcc8/ZfSGF0vN/+umniv9dnz9/vhCiZ89qNpvF4sWLRWJiooiOjhYzZswQZ8+eDcDTeM7d8xsMBpf/LPz000/lzwjV51eiFFBC/fnffPNNcc0114gBAwaI3NxcsX37drvP8NXzS0II0bc5GCIiIiLvYg0KERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpzv8PmbdL1KcRBpsAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGiCAYAAAD5t/y6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAq1UlEQVR4nO3df3DU1b3/8deaH0uIyV6SmGxWAsQxttUN1BscINcrv0Iwlx8qzgXFUpgyjlaIpoAIcmdK79gE6Qh6ZUqvDiMK0jB3JNZbKBIuEmUCNQRzTWhr8TZo0KxpNWwSjJsYzvcPv37GhYBdCISzeT5mPjPuOe/97HmjJC/Pfj67LmOMEQAAgEWu6u8FAAAARIoAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsc1EBpqysTC6XSyUlJc6YMUarV6+Wz+dTQkKCJkyYoKNHj4Y9LxQKqbi4WGlpaUpMTNTMmTN14sSJi1kKAAAYQC44wNTU1Oi5557TyJEjw8bXrl2rdevWacOGDaqpqZHX69WUKVPU3t7u1JSUlKiiokLl5eU6cOCAOjo6NH36dPX09Fx4JwAAYMC4oADT0dGh++67T88//7yGDBnijBtj9PTTT2vVqlWaNWuW/H6/XnzxRX3++efatm2bJCkYDGrTpk166qmnVFBQoJtvvllbt25VfX299u7d2zddAQCAqBZ7IU9atGiRpk2bpoKCAj3xxBPOeGNjowKBgAoLC50xt9ut8ePHq7q6Wg888IBqa2vV3d0dVuPz+eT3+1VdXa2pU6ee9XqhUEihUMh5fPr0aX322WdKTU2Vy+W6kBYAAMBlZoxRe3u7fD6frrrq4i7DjTjAlJeX68iRI6qpqTlrLhAISJIyMjLCxjMyMvTBBx84NfHx8WE7N1/XfP38M5WVlelnP/tZpEsFAABXoKamJg0dOvSizhFRgGlqatIjjzyiPXv2aNCgQeesO3NXxBjzrTsl56tZuXKllixZ4jwOBoMaNmyYmpqalJycHEEHAACgv7S1tSkrK0tJSUkXfa6IAkxtba1aWlqUl5fnjPX09OjNN9/Uhg0b9N5770n6apclMzPTqWlpaXF2Zbxer7q6utTa2hq2C9PS0qL8/PxeX9ftdsvtdp81npycTIABAMAyfXH5R0RvQE2ePFn19fWqq6tzjtGjR+u+++5TXV2drrvuOnm9XlVWVjrP6erqUlVVlRNO8vLyFBcXF1bT3NyshoaGcwYYAACAb4poByYpKUl+vz9sLDExUampqc54SUmJSktLlZOTo5ycHJWWlmrw4MGaO3euJMnj8WjhwoVaunSpUlNTlZKSomXLlik3N1cFBQV91BYAAIhmF3QX0vksX75cnZ2deuihh9Ta2qoxY8Zoz549Ye93rV+/XrGxsZo9e7Y6Ozs1efJkbd68WTExMX29HAAAEIVcxhjT34uIVFtbmzwej4LBINfAAABgib78/c13IQEAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdfr8u5AA4GKNWLGzv5cQseNrpvX3EoABhR0YAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1+BwYIAJ8PgkAXBnYgQEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsE1GA2bhxo0aOHKnk5GQlJydr3Lhx+t3vfufML1iwQC6XK+wYO3Zs2DlCoZCKi4uVlpamxMREzZw5UydOnOibbgAAwIAQUYAZOnSo1qxZo8OHD+vw4cOaNGmS7rjjDh09etSpuf3229Xc3Owcu3btCjtHSUmJKioqVF5ergMHDqijo0PTp09XT09P33QEAACiXmwkxTNmzAh7/POf/1wbN27UoUOHdNNNN0mS3G63vF5vr88PBoPatGmTtmzZooKCAknS1q1blZWVpb1792rq1KkX0gMAABhgLvgamJ6eHpWXl+vUqVMaN26cM75//36lp6frhhtu0P3336+WlhZnrra2Vt3d3SosLHTGfD6f/H6/qqurz/laoVBIbW1tYQcAABi4Ig4w9fX1uvrqq+V2u/Xggw+qoqJCN954oySpqKhIL7/8svbt26ennnpKNTU1mjRpkkKhkCQpEAgoPj5eQ4YMCTtnRkaGAoHAOV+zrKxMHo/HObKysiJdNgAAiCIRvYUkSd/5zndUV1enkydP6pVXXtH8+fNVVVWlG2+8UXPmzHHq/H6/Ro8ereHDh2vnzp2aNWvWOc9pjJHL5Trn/MqVK7VkyRLncVtbGyEGAIABLOIAEx8fr+uvv16SNHr0aNXU1OiZZ57Rf/7nf55Vm5mZqeHDh+vYsWOSJK/Xq66uLrW2tobtwrS0tCg/P/+cr+l2u+V2uyNdKgAAiFIX/TkwxhjnLaIzffrpp2pqalJmZqYkKS8vT3FxcaqsrHRqmpub1dDQcN4AAwAA8E0R7cA8/vjjKioqUlZWltrb21VeXq79+/dr9+7d6ujo0OrVq3X33XcrMzNTx48f1+OPP660tDTdddddkiSPx6OFCxdq6dKlSk1NVUpKipYtW6bc3FznriQAAIBvE1GA+eSTTzRv3jw1NzfL4/Fo5MiR2r17t6ZMmaLOzk7V19frpZde0smTJ5WZmamJEydq+/btSkpKcs6xfv16xcbGavbs2ers7NTkyZO1efNmxcTE9HlzAAAgOkUUYDZt2nTOuYSEBL3++uvfeo5Bgwbp2Wef1bPPPhvJSwMAADj4LiQAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOvE9vcCACAajFixs7+XcEGOr5nW30sALgg7MAAAwDoEGAAAYB3eQgKinK1vbQDA+bADAwAArBNRgNm4caNGjhyp5ORkJScna9y4cfrd737nzBtjtHr1avl8PiUkJGjChAk6evRo2DlCoZCKi4uVlpamxMREzZw5UydOnOibbgAAwIAQUYAZOnSo1qxZo8OHD+vw4cOaNGmS7rjjDiekrF27VuvWrdOGDRtUU1Mjr9erKVOmqL293TlHSUmJKioqVF5ergMHDqijo0PTp09XT09P33YGAACilssYYy7mBCkpKfrFL36hH/3oR/L5fCopKdFjjz0m6avdloyMDD355JN64IEHFAwGdc0112jLli2aM2eOJOnjjz9WVlaWdu3apalTp/5dr9nW1iaPx6NgMKjk5OSLWT4QEa4nQbThNmpcTn35+/uCr4Hp6elReXm5Tp06pXHjxqmxsVGBQECFhYVOjdvt1vjx41VdXS1Jqq2tVXd3d1iNz+eT3+93anoTCoXU1tYWdgAAgIEr4gBTX1+vq6++Wm63Ww8++KAqKip04403KhAISJIyMjLC6jMyMpy5QCCg+Ph4DRky5Jw1vSkrK5PH43GOrKysSJcNAACiSMQB5jvf+Y7q6up06NAh/fjHP9b8+fP1hz/8wZl3uVxh9caYs8bO9G01K1euVDAYdI6mpqZIlw0AAKJIxAEmPj5e119/vUaPHq2ysjKNGjVKzzzzjLxerySdtZPS0tLi7Mp4vV51dXWptbX1nDW9cbvdzp1PXx8AAGDguujPgTHGKBQKKTs7W16vV5WVlc5cV1eXqqqqlJ+fL0nKy8tTXFxcWE1zc7MaGhqcGgAAgG8T0SfxPv744yoqKlJWVpba29tVXl6u/fv3a/fu3XK5XCopKVFpaalycnKUk5Oj0tJSDR48WHPnzpUkeTweLVy4UEuXLlVqaqpSUlK0bNky5ebmqqCg4JI0CAAAok9EAeaTTz7RvHnz1NzcLI/Ho5EjR2r37t2aMmWKJGn58uXq7OzUQw89pNbWVo0ZM0Z79uxRUlKSc47169crNjZWs2fPVmdnpyZPnqzNmzcrJiambzsDAABR66I/B6Y/8Dkw6C98DgyiDZ8Dg8vpivgcGAAAgP5CgAEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsE9vfC8DANWLFzv5eAgDAUuzAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHUiCjBlZWW65ZZblJSUpPT0dN1555167733wmoWLFggl8sVdowdOzasJhQKqbi4WGlpaUpMTNTMmTN14sSJi+8GAAAMCBEFmKqqKi1atEiHDh1SZWWlvvzySxUWFurUqVNhdbfffruam5udY9euXWHzJSUlqqioUHl5uQ4cOKCOjg5Nnz5dPT09F98RAACIerGRFO/evTvs8QsvvKD09HTV1tbqtttuc8bdbre8Xm+v5wgGg9q0aZO2bNmigoICSdLWrVuVlZWlvXv3aurUqWc9JxQKKRQKOY/b2toiWTYAAIgyF3UNTDAYlCSlpKSEje/fv1/p6em64YYbdP/996ulpcWZq62tVXd3twoLC50xn88nv9+v6urqXl+nrKxMHo/HObKysi5m2QAAwHIXHGCMMVqyZIluvfVW+f1+Z7yoqEgvv/yy9u3bp6eeeko1NTWaNGmSs4MSCAQUHx+vIUOGhJ0vIyNDgUCg19dauXKlgsGgczQ1NV3osgEAQBSI6C2kb1q8eLHeffddHThwIGx8zpw5zj/7/X6NHj1aw4cP186dOzVr1qxzns8YI5fL1euc2+2W2+2+0KUCAIAoc0E7MMXFxXrttdf0xhtvaOjQoeetzczM1PDhw3Xs2DFJktfrVVdXl1pbW8PqWlpalJGRcSHLAQAAA0xEAcYYo8WLF2vHjh3at2+fsrOzv/U5n376qZqampSZmSlJysvLU1xcnCorK52a5uZmNTQ0KD8/P8LlAwCAgSiit5AWLVqkbdu26Te/+Y2SkpKca1Y8Ho8SEhLU0dGh1atX6+6771ZmZqaOHz+uxx9/XGlpabrrrruc2oULF2rp0qVKTU1VSkqKli1bptzcXOeuJAAAgPOJKMBs3LhRkjRhwoSw8RdeeEELFixQTEyM6uvr9dJLL+nkyZPKzMzUxIkTtX37diUlJTn169evV2xsrGbPnq3Ozk5NnjxZmzdvVkxMzMV3BAAAop7LGGP6exGRamtrk8fjUTAYVHJycn8vBxdoxIqd/b0EYMA7vmZafy8BA0hf/v7mu5AAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYJ2IAkxZWZluueUWJSUlKT09XXfeeafee++9sBpjjFavXi2fz6eEhARNmDBBR48eDasJhUIqLi5WWlqaEhMTNXPmTJ04ceLiuwEAAANCRAGmqqpKixYt0qFDh1RZWakvv/xShYWFOnXqlFOzdu1arVu3Ths2bFBNTY28Xq+mTJmi9vZ2p6akpEQVFRUqLy/XgQMH1NHRoenTp6unp6fvOgMAAFHLZYwxF/rkv/71r0pPT1dVVZVuu+02GWPk8/lUUlKixx57TNJXuy0ZGRl68skn9cADDygYDOqaa67Rli1bNGfOHEnSxx9/rKysLO3atUtTp0791tdta2uTx+NRMBhUcnLyhS4f/WzEip39vQRgwDu+Zlp/LwEDSF/+/r6oa2CCwaAkKSUlRZLU2NioQCCgwsJCp8btdmv8+PGqrq6WJNXW1qq7uzusxufzye/3OzVnCoVCamtrCzsAAMDAdcEBxhijJUuW6NZbb5Xf75ckBQIBSVJGRkZYbUZGhjMXCAQUHx+vIUOGnLPmTGVlZfJ4PM6RlZV1ocsGAABR4IIDzOLFi/Xuu+/q17/+9VlzLpcr7LEx5qyxM52vZuXKlQoGg87R1NR0ocsGAABR4IICTHFxsV577TW98cYbGjp0qDPu9Xol6aydlJaWFmdXxuv1qqurS62treesOZPb7VZycnLYAQAABq6IAowxRosXL9aOHTu0b98+ZWdnh81nZ2fL6/WqsrLSGevq6lJVVZXy8/MlSXl5eYqLiwuraW5uVkNDg1MDAABwPrGRFC9atEjbtm3Tb37zGyUlJTk7LR6PRwkJCXK5XCopKVFpaalycnKUk5Oj0tJSDR48WHPnznVqFy5cqKVLlyo1NVUpKSlatmyZcnNzVVBQ0PcdAgCAqBNRgNm4caMkacKECWHjL7zwghYsWCBJWr58uTo7O/XQQw+ptbVVY8aM0Z49e5SUlOTUr1+/XrGxsZo9e7Y6Ozs1efJkbd68WTExMRfXDQAAGBAu6nNg+gufAxMd+BwYoP/xOTC4nK6Yz4EBAADoDwQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALBOxAHmzTff1IwZM+Tz+eRyufTqq6+GzS9YsEAulyvsGDt2bFhNKBRScXGx0tLSlJiYqJkzZ+rEiRMX1QgAABg4Ig4wp06d0qhRo7Rhw4Zz1tx+++1qbm52jl27doXNl5SUqKKiQuXl5Tpw4IA6Ojo0ffp09fT0RN4BAAAYcGIjfUJRUZGKiorOW+N2u+X1enudCwaD2rRpk7Zs2aKCggJJ0tatW5WVlaW9e/dq6tSpkS4JAAAMMJfkGpj9+/crPT1dN9xwg+6//361tLQ4c7W1teru7lZhYaEz5vP55Pf7VV1d3ev5QqGQ2trawg4AADBw9XmAKSoq0ssvv6x9+/bpqaeeUk1NjSZNmqRQKCRJCgQCio+P15AhQ8Kel5GRoUAg0Os5y8rK5PF4nCMrK6uvlw0AACwS8VtI32bOnDnOP/v9fo0ePVrDhw/Xzp07NWvWrHM+zxgjl8vV69zKlSu1ZMkS53FbWxshBgCAAazPA8yZMjMzNXz4cB07dkyS5PV61dXVpdbW1rBdmJaWFuXn5/d6DrfbLbfbfamXCgADzogVO/t7CRE7vmZafy8BV4BL/jkwn376qZqampSZmSlJysvLU1xcnCorK52a5uZmNTQ0nDPAAAAAfFPEOzAdHR16//33nceNjY2qq6tTSkqKUlJStHr1at19993KzMzU8ePH9fjjjystLU133XWXJMnj8WjhwoVaunSpUlNTlZKSomXLlik3N9e5KwkAAOB8Ig4whw8f1sSJE53HX1+bMn/+fG3cuFH19fV66aWXdPLkSWVmZmrixInavn27kpKSnOesX79esbGxmj17tjo7OzV58mRt3rxZMTExfdASAACIdi5jjOnvRUSqra1NHo9HwWBQycnJ/b0cXCAb33sH0P+4BsZeffn7m+9CAgAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1Ig4wb775pmbMmCGfzyeXy6VXX301bN4Yo9WrV8vn8ykhIUETJkzQ0aNHw2pCoZCKi4uVlpamxMREzZw5UydOnLioRgAAwMARcYA5deqURo0apQ0bNvQ6v3btWq1bt04bNmxQTU2NvF6vpkyZovb2dqempKREFRUVKi8v14EDB9TR0aHp06erp6fnwjsBAAADRmykTygqKlJRUVGvc8YYPf3001q1apVmzZolSXrxxReVkZGhbdu26YEHHlAwGNSmTZu0ZcsWFRQUSJK2bt2qrKws7d27V1OnTr2IdgAAwEDQp9fANDY2KhAIqLCw0Blzu90aP368qqurJUm1tbXq7u4Oq/H5fPL7/U7NmUKhkNra2sIOAAAwcPVpgAkEApKkjIyMsPGMjAxnLhAIKD4+XkOGDDlnzZnKysrk8XicIysrqy+XDQAALHNJ7kJyuVxhj40xZ42d6Xw1K1euVDAYdI6mpqY+WysAALBPxNfAnI/X65X01S5LZmamM97S0uLsyni9XnV1dam1tTVsF6alpUX5+fm9ntftdsvtdvflUqPOiBU7+3sJAABcNn26A5OdnS2v16vKykpnrKurS1VVVU44ycvLU1xcXFhNc3OzGhoazhlgAAAAviniHZiOjg69//77zuPGxkbV1dUpJSVFw4YNU0lJiUpLS5WTk6OcnByVlpZq8ODBmjt3riTJ4/Fo4cKFWrp0qVJTU5WSkqJly5YpNzfXuSsJAADgfCIOMIcPH9bEiROdx0uWLJEkzZ8/X5s3b9by5cvV2dmphx56SK2trRozZoz27NmjpKQk5znr169XbGysZs+erc7OTk2ePFmbN29WTExMH7QEAACincsYY/p7EZFqa2uTx+NRMBhUcnJyfy/nisA1MAAGiuNrpvX3EnCB+vL3N9+FBAAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHVi+3sBAABEYsSKnf29hIgdXzOtv5cQddiBAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALBOnweY1atXy+VyhR1er9eZN8Zo9erV8vl8SkhI0IQJE3T06NG+XgYAAIhil2QH5qabblJzc7Nz1NfXO3Nr167VunXrtGHDBtXU1Mjr9WrKlClqb2+/FEsBAABR6JIEmNjYWHm9Xue45pprJH21+/L0009r1apVmjVrlvx+v1588UV9/vnn2rZt26VYCgAAiEKXJMAcO3ZMPp9P2dnZuueee/SXv/xFktTY2KhAIKDCwkKn1u12a/z48aqurj7n+UKhkNra2sIOAAAwcPV5gBkzZoxeeuklvf7663r++ecVCASUn5+vTz/9VIFAQJKUkZER9pyMjAxnrjdlZWXyeDzOkZWV1dfLBgAAFunzAFNUVKS7775bubm5Kigo0M6dOyVJL774olPjcrnCnmOMOWvsm1auXKlgMOgcTU1Nfb1sAABgkUt+G3ViYqJyc3N17Ngx526kM3dbWlpaztqV+Sa3263k5OSwAwAADFyXPMCEQiH98Y9/VGZmprKzs+X1elVZWenMd3V1qaqqSvn5+Zd6KQAAIErE9vUJly1bphkzZmjYsGFqaWnRE088oba2Ns2fP18ul0slJSUqLS1VTk6OcnJyVFpaqsGDB2vu3Ll9vRQAABCl+jzAnDhxQvfee6/+9re/6ZprrtHYsWN16NAhDR8+XJK0fPlydXZ26qGHHlJra6vGjBmjPXv2KCkpqa+XAgAAopTLGGP6exGRamtrk8fjUTAY5HqY/2/Eip39vQQAwDkcXzOtv5dwRejL3998FxIAALAOAQYAAFiHAAMAAKxDgAEAANYhwAAAAOsQYAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BBgAAWIcAAwAArEOAAQAA1iHAAAAA6xBgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrxPb3Aq5EI1bs7O8lAACA82AHBgAAWIcdGAAALjEbd/aPr5nW30s4L3ZgAACAdQgwAADAOgQYAABgHQIMAACwDgEGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALBOvwaYX/7yl8rOztagQYOUl5ent956qz+XAwAALNFvAWb79u0qKSnRqlWr9M477+if//mfVVRUpA8//LC/lgQAACzhMsaY/njhMWPG6B//8R+1ceNGZ+x73/ue7rzzTpWVlYXVhkIhhUIh53EwGNSwYcPU1NSk5OTkPl+b/6ev9/k5AQCwScPPpvb5Odva2pSVlaWTJ0/K4/Fc1Lli+2hNEenq6lJtba1WrFgRNl5YWKjq6uqz6svKyvSzn/3srPGsrKxLtkYAAAYyz9OX7tzt7e12Bpi//e1v6unpUUZGRth4RkaGAoHAWfUrV67UkiVLnMenT5/WZ599ptTUVLlcrku+3vP5Ok1eqt2gKx390z/90z/90//f278xRu3t7fL5fBf9+v0SYL52ZvgwxvQaSNxut9xud9jYP/zDP1zKpUUsOTl5QP4H/DX6p3/6p/+Biv4j6/9id16+1i8X8aalpSkmJuas3ZaWlpazdmUAAADO1C8BJj4+Xnl5eaqsrAwbr6ysVH5+fn8sCQAAWKTf3kJasmSJ5s2bp9GjR2vcuHF67rnn9OGHH+rBBx/sryVdELfbrZ/+9KdnvcU1UNA//dM//dM//feHfruNWvrqg+zWrl2r5uZm+f1+rV+/Xrfddlt/LQcAAFiiXwMMAADAheC7kAAAgHUIMAAAwDoEGAAAYB0CDAAAsA4BRl9919Itt9yipKQkpaen684779R7770XVmOM0erVq+Xz+ZSQkKAJEybo6NGjYTWhUEjFxcVKS0tTYmKiZs6cqRMnToTVtLa2at68efJ4PPJ4PJo3b55Onjx5qVv8u5WVlcnlcqmkpMQZi/beP/roI/3gBz9QamqqBg8erO9///uqra115qO5/y+//FL/9m//puzsbCUkJOi6667Tv//7v+v06dNOTTT1/+abb2rGjBny+XxyuVx69dVXw+YvZ68ffvihZsyYocTERKWlpenhhx9WV1fXpWjbcb7+u7u79dhjjyk3N1eJiYny+Xz64Q9/qI8//jjsHNHa/5keeOABuVwuPf3002Hj0d7/H//4R82cOVMej0dJSUkaO3asPvzwQ2f+iurfwEydOtW88MILpqGhwdTV1Zlp06aZYcOGmY6ODqdmzZo1Jikpybzyyiumvr7ezJkzx2RmZpq2tjan5sEHHzTXXnutqaysNEeOHDETJ040o0aNMl9++aVTc/vttxu/32+qq6tNdXW18fv9Zvr06Ze133N5++23zYgRI8zIkSPNI4884oxHc++fffaZGT58uFmwYIH5/e9/bxobG83evXvN+++/79REc/9PPPGESU1NNb/97W9NY2Oj+a//+i9z9dVXm6efftqpiab+d+3aZVatWmVeeeUVI8lUVFSEzV+uXr/88kvj9/vNxIkTzZEjR0xlZaXx+Xxm8eLF/db/yZMnTUFBgdm+fbv505/+ZA4ePGjGjBlj8vLyws4Rrf1/U0VFhRk1apTx+Xxm/fr1YXPR3P/7779vUlJSzKOPPmqOHDli/u///s/89re/NZ988skV2T8BphctLS1GkqmqqjLGGHP69Gnj9XrNmjVrnJovvvjCeDwe86tf/coY89Vf/ri4OFNeXu7UfPTRR+aqq64yu3fvNsYY84c//MFIMocOHXJqDh48aCSZP/3pT5ejtXNqb283OTk5prKy0owfP94JMNHe+2OPPWZuvfXWc85He//Tpk0zP/rRj8LGZs2aZX7wgx8YY6K7/zN/gF/OXnft2mWuuuoq89FHHzk1v/71r43b7TbBYPCS9Hum8/0C/9rbb79tJJkPPvjAGDMw+j9x4oS59tprTUNDgxk+fHhYgIn2/ufMmeP83e/NldY/byH1IhgMSpJSUlIkSY2NjQoEAiosLHRq3G63xo8fr+rqaklSbW2turu7w2p8Pp/8fr9Tc/DgQXk8Ho0ZM8apGTt2rDwej1PTXxYtWqRp06apoKAgbDzae3/ttdc0evRo/eu//qvS09N188036/nnn3fmo73/W2+9Vf/zP/+jP//5z5Kk//3f/9WBAwf0L//yL5Kiv/9vupy9Hjx4UH6/P+wbeadOnapQKBT29mV/CwaDcrlczpfnRnv/p0+f1rx58/Too4/qpptuOms+mvs/ffq0du7cqRtuuEFTp05Venq6xowZE/Y205XWPwHmDMYYLVmyRLfeeqv8fr8kOV86eeYXTWZkZDhzgUBA8fHxGjJkyHlr0tPTz3rN9PT0s77Y8nIqLy/XkSNHVFZWdtZctPf+l7/8RRs3blROTo5ef/11Pfjgg3r44Yf10ksvSYr+/h977DHde++9+u53v6u4uDjdfPPNKikp0b333isp+vv/psvZayAQOOt1hgwZovj4+Cvmz+OLL77QihUrNHfuXOebhqO9/yeffFKxsbF6+OGHe52P5v5bWlrU0dGhNWvW6Pbbb9eePXt01113adasWaqqqpJ05fXfb9+FdKVavHix3n33XR04cOCsOZfLFfbYGHPW2JnOrOmt/u85z6XS1NSkRx55RHv27NGgQYPOWReNvUtf/V/H6NGjVVpaKkm6+eabdfToUW3cuFE//OEPnbpo7X/79u3aunWrtm3bpptuukl1dXUqKSmRz+fT/Pnznbpo7b83l6vXK/nPo7u7W/fcc49Onz6tX/7yl99aHw3919bW6plnntGRI0ciXkM09P/1hft33HGHfvKTn0iSvv/976u6ulq/+tWvNH78+HM+t7/6ZwfmG4qLi/Xaa6/pjTfe0NChQ51xr9crSWclw5aWFidFer1edXV1qbW19bw1n3zyyVmv+9e//vWsNHq51NbWqqWlRXl5eYqNjVVsbKyqqqr0H//xH4qNjXXWFY29S1JmZqZuvPHGsLHvfe97zlX30fzvXpIeffRRrVixQvfcc49yc3M1b948/eQnP3F246K9/2+6nL16vd6zXqe1tVXd3d39/ufR3d2t2bNnq7GxUZWVlc7uixTd/b/11ltqaWnRsGHDnJ+FH3zwgZYuXaoRI0ZIiu7+09LSFBsb+60/D6+k/gkw+ir1LV68WDt27NC+ffuUnZ0dNp+dnS2v16vKykpnrKurS1VVVcrPz5ck5eXlKS4uLqymublZDQ0NTs24ceMUDAb19ttvOzW///3vFQwGnZrLbfLkyaqvr1ddXZ1zjB49Wvfdd5/q6up03XXXRW3vkvRP//RPZ90y/+c//1nDhw+XFN3/7iXp888/11VXhf8YiImJcf5vLNr7/6bL2eu4cePU0NCg5uZmp2bPnj1yu93Ky8u7pH2ez9fh5dixY9q7d69SU1PD5qO5/3nz5undd98N+1no8/n06KOP6vXXX5cU3f3Hx8frlltuOe/Pwyuu/7/7ct8o9uMf/9h4PB6zf/9+09zc7Byff/65U7NmzRrj8XjMjh07TH19vbn33nt7vb1y6NChZu/evebIkSNm0qRJvd5eNnLkSHPw4EFz8OBBk5ub2++30p7pm3chGRPdvb/99tsmNjbW/PznPzfHjh0zL7/8shk8eLDZunWrUxPN/c+fP99ce+21zm3UO3bsMGlpaWb58uVOTTT1397ebt555x3zzjvvGElm3bp15p133nHusrlcvX59G+nkyZPNkSNHzN69e83QoUMv+W205+u/u7vbzJw50wwdOtTU1dWF/SwMhUJR339vzrwLyZjo7n/Hjh0mLi7OPPfcc+bYsWPm2WefNTExMeatt966IvsnwJivbifr7XjhhRecmtOnT5uf/vSnxuv1GrfbbW677TZTX18fdp7Ozk6zePFik5KSYhISEsz06dPNhx9+GFbz6aefmvvuu88kJSWZpKQkc99995nW1tbL0OXf78wAE+29//d//7fx+/3G7Xab7373u+a5554Lm4/m/tva2swjjzxihg0bZgYNGmSuu+46s2rVqrBfWNHU/xtvvNHr3/X58+cbYy5vrx988IGZNm2aSUhIMCkpKWbx4sXmiy++uJTtn7f/xsbGc/4sfOONN6K+/970FmCivf9NmzaZ66+/3gwaNMiMGjXKvPrqq2HnuJL6dxljzN+/XwMAAND/uAYGAABYhwADAACsQ4ABAADWIcAAAADrEGAAAIB1CDAAAMA6BBgAAGAdAgwAALAOAQYAAFiHAAMAAKxDgAEAANb5f62eRuRSmYpfAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMgklEQVR4nO3de3hU1b0//veeBEKSJmMukGQkgRDBqISQQopAysULiiJFPBXEC7Xql5aLxgsSsLbSnzIJPUVb8ALUox49XL7PQ0K1goJFAjk5X4iBAEEBjQMEmJyYGieQjAlm1u+POJu57JlkkrnsmXm/nifPafbsmax9zqnzdq3P+ixJCCFAREREpCKaQA+AiIiIyBEDChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREalOZKAH0BsWiwUXLlxAXFwcJEkK9HCIiIioB4QQuHjxInQ6HTQa93MkQRlQLly4gPT09EAPg4iIiHqhvr4egwcPdntPUAaUuLg4AF0PGB8fH+DREBERUU+0tLQgPT1d/h53JygDinVZJz4+ngGFiIgoyPSkPINFskRERKQ6DChERESkOgwoREREpDoeB5R9+/bhrrvugk6ngyRJ2L59u9M9X3zxBWbOnAmtVou4uDjceOONOHv2rPx6e3s7lixZguTkZMTGxmLmzJk4d+5cnx6EiIiIQofHAaW1tRW5ublYt26d4ut1dXUoKChAdnY29u7diyNHjuD555/HgAED5HsKCwtRVlaGLVu2oKKiApcuXcKMGTPQ2dnZ+ychIiKikCEJIUSv3yxJKCsrw6xZs+Rrc+fORb9+/fDuu+8qvsdkMmHgwIF49913MWfOHABX+prs2LEDt912W7d/t6WlBVqtFiaTibt4iIiIgoQn399erUGxWCz48MMPMWLECNx2220YNGgQxo0bZ7cMVF1djcuXL2PatGnyNZ1Oh5EjR6KyslLxc9vb29HS0mL3Q0RERKHLqwGlsbERly5dQnFxMW6//Xbs2rULd999N2bPno3y8nIAQENDA/r374+EhAS796akpKChoUHxc/V6PbRarfzDLrJEREShzeszKADwi1/8Ak8++SRGjx6NoqIizJgxA2+88Ybb9wohXDZuWb58OUwmk/xTX1/vzWETERGRyng1oCQnJyMyMhLXX3+93fXrrrtO3sWTmpqKjo4ONDc3293T2NiIlJQUxc+NioqSu8ayeywREVHo82pA6d+/P/Lz83Hy5Em766dOncKQIUMAAGPGjEG/fv2we/du+XWj0Yja2lpMmDDBm8MhIiKiXjCazKisa4LRZA7YGDw+i+fSpUv46quv5N8NBgNqamqQmJiIjIwMLF26FHPmzMGkSZMwdepUfPTRR/jggw+wd+9eAIBWq8UjjzyCp59+GklJSUhMTMQzzzyDnJwc3HLLLV57MCIiIvLc1qqzWF56DBYBaCRAPzsHc/Iz/D4Oj7cZ7927F1OnTnW6Pn/+fLz99tsAgP/4j/+AXq/HuXPncO2112LlypX4xS9+Id/7/fffY+nSpdi0aRPMZjNuvvlmvPbaaz0ufuU2YyIiIu8zmsyYWLwHFptkECFJqCiaijRtdJ8/35Pv7z71QQkUBhQiIiLvq6xrwryNB5yub37sRozPSurz5wesDwoREREFr8zkWGgcNtRGSBKGJsf4fSwMKERERAQASNNGQz87BxE/tv2IkCSsmj3SK8s7nvK4SJaIiIhC15z8DEwaMRCnm9owNDkmIOEEYEAhIiIiB2na6IAFEysu8RAREZHqMKAQERGR6jCgEBERkeowoBAREZHqMKAQERGR6jCgEBER+ZgaDt9Twxg8wW3GREREPqSGw/fUMAZPcQaFiIjIR4wmsxwMAMAigBWltX6dxVDDGHqDAYWIiMhHDE2tdicDA0CnEDjd1Nbrz/R0qcYXY/AHLvEQERH5iPXwPduA0JfD93qzVOPtMfgLZ1CIiIh8xJuH7/V2qUZNBwB6gjMoREREPuStw/fcLdV095lqOQDQEwwoREREPuaNw/f6ulSjhgMAPcElHiIiIh/yVv+RYF2q6S3OoBAREfmIbVGrBOCxn2fi4YLMXoeKYFyq6S1JCCG6v01dWlpaoNVqYTKZEB8fH+jhEBEROTGazJhYvMepbiRYGqX5giff31ziISIi8gGlolbAefdNsLWg9xcu8RAREfmAUlGrlXX3zb5T39j1NVl2ezZyBmuRmRwb0ss3PcEZFCIiol7obubDWtSq9EWrkYCvvrmIom32fU30O09g3sYDmFi8B1urzvpu8EGANShEREQe8qSjq9FkxlsVp7Fx/9ewfuFKALr78o2QJFQUTQ2pmRTWoBAREXnAkzoQTzu6pmmjkTUoFj/uDgbQfTgBguO8HF9iDQoREYU1T8+38bSjq2Og6algOC/HlziDQkREYas359tYi19tuQsTrnbzWGkkYN19eVh+R3bYNGHrCc6gEBFR2Ort+TaPFmTib/sNsKD7MOFuN4/1vTNydQCAmbk6l03YjCYzDE2tYbPDhwGFiIjClqfn2zh2hr1/XAYW33SN28Bg3c2zfNsxWNC1dLHsjmyMuvoqpyDi6rwcT5ehQgGXeIiIKGx1d76NbfGs43KQAPBfB87ivf8507M/Jl35n1dF98P4rKQezYT0ZhkqFHAGhYiIwpqr820cZy0eLchUXKZ5dW8d4mP6YcGkLMXPdxUwJo0Y2KOA0ttlqGDHGRQiIgp7adpouxkNpVDxt/0GSC7eX7LzhMsZDXcBoyc8LcoNFQwoREREDpRChQXAvHHKdR8WAZeBo68Bo7tlqFDFJR4iIiIHSsWzGgm4cVgiIjUS3nGoO3EXOKwBY0VpLTqF6FXAcLUMFcrY6p6IiEjB1qqzcqiQJMD6bSkBmD4yFR8db4BFXJnR6G5XjdFkDquAocST728GFCIiIheMJjMOnWnGok2H7a5LALYvmoBzzWZYhMDYoYlhGzo84dOzePbt24e77roLOp0OkiRh+/btLu9dsGABJEnCK6+8Yne9vb0dS5YsQXJyMmJjYzFz5kycO3fO06EQERH5jLUx2r9a251eEwC2Vp3Dks2HsWRzDU8f9gGPA0praytyc3Oxbt06t/dt374dBw4cgE6nc3qtsLAQZWVl2LJlCyoqKnDp0iXMmDEDnZ2dng6HiIjI67ZWncXE4j2Yt/EA/vD3zxXv2XzwbNj1JvEnj4tkp0+fjunTp7u95/z581i8eDE+/vhj3HnnnXavmUwmvPnmm3j33Xdxyy23AADee+89pKen45NPPsFtt93m6ZCIiIgUedoe3mgyo/pMM4q2HZNPHFaqg5AUrodDbxJ/8vouHovFggcffBBLly7FDTfc4PR6dXU1Ll++jGnTpsnXdDodRo4cicrKSgYUIiLyCk/bw9ve745GApZNz0bJzhM9bpFPnvN6QCkpKUFkZCQef/xxxdcbGhrQv39/JCQk2F1PSUlBQ0OD4nva29vR3n5lDbClpcV7AyYiopDjafdWx/uVaCTgr3PzMGZoAtK00bgqul+ftg6Te14NKNXV1fjLX/6CQ4cOQZJc9dtTJoRw+R69Xo+VK1d6Y4hERBQGPGkPbzSZ8Y+jF7qdObEIIOknUfL7w7E3iT95tZPs/v370djYiIyMDERGRiIyMhJnzpzB008/jaFDhwIAUlNT0dHRgebmZrv3NjY2IiUlRfFzly9fDpPJJP/U19d7c9hERBRietq91VoM+9KHJ7r9TKX3O7bIJ+/xakB58MEHcfToUdTU1Mg/Op0OS5cuxccffwwAGDNmDPr164fdu3fL7zMajaitrcWECRMUPzcqKgrx8fF2P0RERK6kaaOxbHq2HFKUlmB6sqxjJUngEo6febzEc+nSJXz11Vfy7waDATU1NUhMTERGRgaSkpLs7u/Xrx9SU1Nx7bXXAgC0Wi0eeeQRPP3000hKSkJiYiKeeeYZ5OTkyLt6iIiI+mJr1Vm5iFUC8Ozt1zoVyL5VYehROAEASQCTRgz0/kDJJY9nUD777DPk5eUhLy8PAPDUU08hLy8Pv//973v8GS+//DJmzZqFe++9FxMnTkRMTAw++OADREREeDocIiIiO44zIwLA6o9O2vUoOVLfjA37DT3+TAtcHwZIvuHxDMqUKVPgSXf806dPO10bMGAA1q5di7Vr13r654mIiNzqrkB2a9VZLNt2zKPP5BZi//NqDQoREZG/GE1mVNY1OXVvdVcgazSZUeRhOAGAWXk61p/4GQMKEREFHdtW9I7n4KRpo3F33tV291sDxmenv1XsDOvIsenF9sMX5CDkKhiRd3m9URsREZEvddeEzWgyo+zwebv3lB0+jynXDsTZf3VfR+Kujf2+U9941J2Weo8BhYiIgkp3NSafnf7W6XWLAJZsrunR5yvNsERIEmL6azzqTkt9wyUeIiIKKko1JhoAQ5NjsL68rsdBpDvWL0hrD5XWjk6XwYi8jzMoREQUVNK00dDPznE6cXhF2TF8euKbHn2G0jKOrQhJQunC8WjrsMht7I0mMzQSeECgnzCgEBGRKhlNZhiaWpGZHOu0hOLYNE0APQ4ntuHjv7/6Bq/urYNwCB2rZo9Ebrr9obbWYMQDAv2DAYWIiFRn/b46FO88AeGiGLX6THOPduM40gBy+DCazLj/b/bhBAB+M3mYy8JXHhDoPwwoRESkKuvL66DfeeXwPosAirYdQ8cPFowarEVrRye+bW3v1Wf/dkqWHD6Uim0B4PW9dXhg/BCX4SNNG81g4gcMKERE5Bfulmxs7yne6XyysADw/N+Py787Fsn21BvlX8vhIzM5FpIEpxkUa1t7hpDA4i4eIiLyOXeN1ayMJjP+cfRCj5ZuenrInyPbXTdp2mgUTc92uoeFr+rAgEJERD7lqrGabSdWa4B56UPn2RN3ejORcvT8d/J/XjApC8unZzttKebsSeBxiYeIiHyqu8ZqjgGmpzQAfjs1C6/vrfPovat3nsTM3Ctn6yyYnIWZo3UsfFUZzqAQEZFPuTu8D3BdrNodC7oKWqdcO7Dbe20pNVdL00ZjfFYSw4mKMKAQEZFPpWmjsex218somcmxvVqqAbqWi3ra/8RKAlhjEgQYUIiIyKe2Vp1FyUcnYEFXOJiTP9ip0VpfeDz50ts0RH7FgEJERD7jWF8iAGw6WI8J+is7eQxNrb1qutZbQoDn5wQBBhQiIvIKo8mMyromu905b1UYFOtLBK7s5FGqUfGEBGDR1Kwe389txMGBAYWIiPpMqc+J0WTGhv0Gl++x3cmjn53Tq5UXDYDie3LwwI1DFF+fNy4dy6dnI0Lq+nRuIw4e3GZMRER94qrPyR9mXtftez84eh5Dk2MwJz8D7T904vd//7zHf1cDoGzRBOSmJ6CyrknxnrtGXY3xWUncRhyEGFCIiKhPXPU5+deljm7fu+lAPTYfrMejBZnQRvfz6O9aALR1WABc2cpscTiV2LqUw/Nzgg+XeIiIqE9c9TkZNVjbo/cLAWzcb8C/7zrl0d91DCD62TlcygkhnEEhIqI+sYaDFaW16BQCGgCPFAxFW0enz/6mBnAKIHPyMzBpxEAu5YQISQjHcxzVr6WlBVqtFiaTCfHx8YEeDhERAThS34wN+wzYUWuEEF27a3zxBSMB2P5j7QkFF0++vzmDQkREfba16iyKth2zCyS++rff4ntyGE7CAAMKERH1iXUXjz+m4//OmZOwwSJZIiLqk94c9neDLh5SLxqfWHftUOjjDAoREXnMaDLD0NSKzORY+bA/TzLK8QstADyrU9EAiOnPf68OFyySJSIij2ytOis3ZpMALJyShW/bOrD5YH2vPk+SurYa94RGAvSzczAnP6NXf4sCy5Pvb0ZRIiLqEaPJjH8cvYCibfaH/726t67X4QRwHU40EvDm/DF2S0HWLrW25/1QaOISDxERdWt9eR2Kd57w+6nDn3zR6BRgbM/wodDFgEJERG6t31cH/c4Tfv+7AlCcmeFpxOGBSzxEROTSkfpm6Hd4N5z05YtHqYMshSYGFCIiUrR+Xx1mvVrp1c9cd18eyhZNcDq7p6fWzstjgWyYYEAhIiIn68vroN/h/ZqTT082Ijc9we5gP43UtRuoO5IE/HQIm7SFC9agEBGRnSP1zT6rOdl26DweGj/E6WC/fae+sTtscO64dGw+UG8fkIKuKQb1BQMKERHJ1pf7viD2s9PNyE1PQJo2Wq4lcQwshqZWbDpgXyArAO7eCSMeL/Hs27cPd911F3Q6HSRJwvbt2+XXLl++jGXLliEnJwexsbHQ6XR46KGHcOHCBbvPaG9vx5IlS5CcnIzY2FjMnDkT586d6/PDEBFR7/lrt87YocrLNGnaaIzPSkKaNhqZybFOdSrcvRNePA4ora2tyM3Nxbp165xea2trw6FDh/D888/j0KFDKC0txalTpzBz5ky7+woLC1FWVoYtW7agoqICly5dwowZM9DZ2dn7JyEioh4xmsyorGuya3ZmNJlR7KetxIPiB3R7T5o22q5OJUKSuHsnzPSp1b0kSSgrK8OsWbNc3lNVVYWf/exnOHPmDDIyMmAymTBw4EC8++67mDNnDgDgwoULSE9Px44dO3Dbbbd1+3fZ6p6IqHds29Tbto2vrGvCvI0H/DKGzY/diPFZST2612gyy8s+DCfBT1Wt7k0mEyRJwlVXXQUAqK6uxuXLlzFt2jT5Hp1Oh5EjR6KyUnk7W3t7O1paWux+iIjIM0aTWQ4ngH3beKUlFV/wdJnGdtmHwotPi2S///57FBUVYd68eXJSamhoQP/+/ZGQYL8GmZKSgoaGBsXP0ev1WLlypS+HSkQU8gxNrXI4seoUAh8eNaJTiB4f2Ofozflj8P1lC4QAzpvMWL3zJDqFQIQkYVaeDtsPX5B/5zIN9ZTPAsrly5cxd+5cWCwWvPbaa93eL4SAJCnH9+XLl+Opp56Sf29paUF6errXxkpEFA6OnTcpXn/xwy/69LkHDN9ixR3XA+iapdFpB0AjSfjpkK6dOs/cdi1ON7Uhpr8GrR2dMJrMDCnULZ8ElMuXL+Pee++FwWDAnj177NaZUlNT0dHRgebmZrtZlMbGRkyYMEHx86KiohAVFeWLoRIRhQWjyYxiL7est9q4z4CHJ2Zi36lvFOtb0rTRLl8jcsXrNSjWcPLll1/ik08+QVKSfSHUmDFj0K9fP+zevVu+ZjQaUVtb6zKgEBGRZxx36lSfafZZnzMBoPp0s8v6Fne1L0SueDyDcunSJXz11Vfy7waDATU1NUhMTIROp8O//du/4dChQ/jHP/6Bzs5Oua4kMTER/fv3h1arxSOPPIKnn34aSUlJSExMxDPPPIOcnBzccsst3nsyIqIw5bhTZ9n0bDSYvvfp39z9+f8q1recbmrDpydcv8alHnLF44Dy2WefYerUqfLv1tqQ+fPn44UXXsD7778PABg9erTd+z799FNMmTIFAPDyyy8jMjIS9957L8xmM26++Wa8/fbbiIiI6OVjEBERoLxTx9unESv5+5ELkGDfjT5CkhDTX4ON+w1O92sksOkaueVxQJkyZQrctU7pSVuVAQMGYO3atVi7dq2nf56IiBwYTWYYmlqRmRyruFPHnzRSVyiy7thp7ehUXFp6tGAYZ0/ILZ7FQ0QUxGyXcyQAD40f4jST4U2SBJfbkQWAtXPzkPSTKLmxmtFklkOLlQbAwwVDfTRCChU+b9RGRES+YTSZUbTtynKOAPDO/5zxWTiJkCQUTc+W288rSU+0b6ym1LJef08OZ0+oW5xBISIKUp98/r8+CyNWGgAWXFmymZOfgZm5Onx41KjYP6Wtw+J0zfGkYoYT6gkGFCKiILR+X53Pi19/NX4IFkzJcgoWadpo3DkqDat2fGG3dOOujX2aNprBhDzCJR4ioiCzvtz34QQA6pvNLs/C4WnD5Gt9Os04UHiaMRGFK6PJjAn6PT5f2rFaNCULE4cnIzM5VjF88LRh8oQn399c4iEiUinb7cPWL39DU6vfwgkAvLq3Dq/urXPZnp5LN+QrDChERCrk2A3WGg427vs6IOOxtqefNGIgAwn5BQMKEZHKKHWDXV56DN9casenJ78J2LjYnp78iQGFiEglrEs6/7rU7tQN1iKAf//4VGAG9iN3u3SIvI0BhYhIBbZWnUXRtmMQ6OoI68tusL3BXTrkbwwoREQBZu0Iaw0k1v+phpCiAbB2Xh5+OiSB4YT8igGFiChArEs6dd9cUgwiw1Nicep/W/0+LisJgP6eHNw5ShewMVD4YkAhIgoA2106rgQynADAL0brnLYVE/kLO8kSEfmJ0WTGP45ewH/+j6HbcKIG79dcgNFkDvQwKExxBoWIyA9si2CDhQXgtmIKGM6gEBH5mLWvSaDCidSDezQKN3FbMQUSAwoRkY8ZmloDtpwza7QOUg8SyqMFw1ByDw//I/XgEg8RkY9lJsdCIyEgIeX9Ixe6/bsaAA8XDEWaNhqTRgzk4X+kCpxBISLysTRtNPSzcwLyty3CeYlHwpUlnQhJgv6eHDmMpGmjMT4rieGEAo4zKEREXuR4ArHRZMZnp79FbFQk3pw/Bo+8U+3X8UgAFk7JwhvlX6NTCGgk4JGCTMwYlYa2DgtnSki1JCFEMBWVAwBaWlqg1WphMpkQHx8f6OEQEQFwPoH47ryrse3Q+UAPCxoJWDY9G/+61I6N+wwQsD8hmchfPPn+ZkAhIuol29kSAJhYvEe1vU00EiCEfev8CElCRdFUzqCQ33jy/c0lHiKiHrIu13zXdhlfGFuwpapervG4MydVteEEUC7Q7RSCfU5ItRhQiIh6wF2jNQHgH8ca/D0kj7iaQWGfE1Ir7uIhIuqG42nDambdnSNJV3bvREgS9LNzUMw+JxREOINCRNQNQ1OrR+FEAvweZjToOnnYto8JAKeeJuxzQsGCAYWIqBuZybEehQ5/h5N5P8vAkpuvsetlYuUYQtK00QwmFBS4xENE1I00bTSKpmcHehgu3ZWrY+igkMMZFCKibvzp4xN49dM6AF3LN78YrcMt16XAaDLjpR0n/DYOpVkcFrpSqGJAISJyY+F71dhRe2WHjkDX+TZGkxkHDM1+G8e8celYctNwvH/kAkp2noBFsNCVQhsDChGRC0fqm+3CiZVFwOvhpLsaly0H65E7+CosmJSFmbk6FrpSyGMNChGFNaPJjMq6JhhNZrv/DAAHT3/rt3H8fHiy29ctAlhRWgujycwD/SgscAaFiMLW+vI6FO88AYErPUOs59Qsm56NQ6f9t4Sz/8umbu9h51cKJwwoRBR2jCYz1u75EpsO1MvXbJdXLALQ+7H41fHvu8KCWAonDChEFFZsTxwOJiyIpXDDgEJEYcNoMqs6nGgkyLtznp1+LUZdfRVi+mvQ1mFhQSyFHY+LZPft24e77roLOp0OkiRh+/btdq8LIfDCCy9Ap9MhOjoaU6ZMwfHjx+3uaW9vx5IlS5CcnIzY2FjMnDkT586d69ODEBF1x9DUqtpwIgEoWzgBmx+7ERVFU7FgUhbGZyUhNz2BBbEUljwOKK2trcjNzcW6desUX1+9ejXWrFmDdevWoaqqCqmpqbj11ltx8eJF+Z7CwkKUlZVhy5YtqKiowKVLlzBjxgx0dnb2/kmIiLqRmRwrH6anNrN/ejXDCJENSQjR63+fkCQJZWVlmDVrFoCu2ROdTofCwkIsW7YMQNdsSUpKCkpKSrBgwQKYTCYMHDgQ7777LubMmQMAuHDhAtLT07Fjxw7cdttt3f7dlpYWaLVamEwmxMfH93b4RBRijCYzDE2tyEyOBdA1YxLbPwL1zWYIIZCRGIO3/vs0ttdcCPBInUVIEiqKpjKcUEjz5PvbqzUoBoMBDQ0NmDZtmnwtKioKkydPRmVlJRYsWIDq6mpcvnzZ7h6dToeRI0eisrJSMaC0t7ejvb1d/r2lpcWbwyaiEGBb/Gq7ZThYcAsxkT2vNmpraOjquJiSkmJ3PSUlRX6toaEB/fv3R0JCgst7HOn1emi1WvknPT3dm8MmoiDnWPwqoI5wsnx6Nh77eabia44rTdxCTGTPJ51kJcn+v3pCCKdrjtzds3z5cphMJvmnvr5e8T4iCk9qLH59cdYNWDA5C78uyFQMI0XTsxHx4z/zuIWYyJlXl3hSU1MBdM2SpKWlydcbGxvlWZXU1FR0dHSgubnZbhalsbEREyZMUPzcqKgoREVFeXOoRBRCrMWvagopv9t+HK0dnVgwKQvF9+RgRWktOoWQw8ic/AzMHM0zdYhc8eoMSmZmJlJTU7F79275WkdHB8rLy+XwMWbMGPTr18/uHqPRiNraWpcBhYioO48WZKrucDH9jhNYX16HOfkZqCiaKm8hnjRiICrrulrbc9cOkTKPZ1AuXbqEr776Sv7dYDCgpqYGiYmJyMjIQGFhIVatWoXhw4dj+PDhWLVqFWJiYjBv3jwAgFarxSOPPIKnn34aSUlJSExMxDPPPIOcnBzccsst3nsyIgoLtsWxGgmYl5+BzQfPqqIGBQBKdp7AzNE6pGmjkaaNdhqvfnYO5uRnBHqYRKrjcUD57LPPMHXqVPn3p556CgAwf/58vP3223j22WdhNpuxcOFCNDc3Y9y4cdi1axfi4uLk97z88suIjIzEvffeC7PZjJtvvhlvv/02IiIivPBIRBQuHItjLQKqCicAYAHk3TlK411RWotJIwZyFoXIQZ/6oAQK+6AQEQBU1jVh3sYDgR5Gt/6+aAJy0xNcjnfzYzdifFZSAEZG5F8B64NCROQP1oZssf0jIMF+S7Hj72rQ1mEBoFzMy+3FRMoYUIgoqDg2ZHMMI2oLJ7YBJE0bDf1s5x09XN4hcsaAQkRBw2gyo2jbMTmEqC2M/H3RBJxouOg2gMzJz8CkEQO5vZioGwwoRKR6RpMZn53+FofOfKe6UGIlScCg+AHITU/oNoBYd/QQkWsMKESkalurztrNmqiVEFd26zCAEPUdAwoRqYrjicTLth0L8IiUOda/sNiVyLsYUIhINWxnSyQABcOTAz0klx77+TC8WWFgsSuRjzCgEFFA2W4ZdiyA3f9lUyCH5lKEJOHhgqF4uGAoi12JfIQBhYgCxnbLcLBwnC1hMCHyDQYUIgoIx7bvweD5O6/DHaPSGEqI/EBth38SUZgwNLUGVTiJkCSGEyI/YkAhooDITI6FFOhB9BCLYIn8j0s8ROQztluGrV/utkWxamQ9KydCkvDs9Gsx6uqrWARLFAAMKETkE7YFsBoJ0M/OAQC35+gEWoQkoXTheLR1WBhKiAKMAYWIvM6xANYigOU/biFW2zk61qCkAfDs7dciNz0hwCMiIoA1KETUR0aTGZV1TTCazPI1pQJYC9QTSqwiJAmLpmZBkrrGV/LRCWytOhvoYREROINCRH2gtIwzJz9Dsb5ELUs6GnSFkQhJwrO3X4uSj05A2Mz0rCitxaQRA7m8QxRgDChE1CtKyzgrSmvxXdtl6HeecLr/9pGp+Ki2IeAhpWzRBLnGRGmmp1MI+dA/IgocLvEQUa+4+nJXCicA8HFtA4ruyPbDyNxr67BgfFYS0rTRyEyOhcZhrzMP/SNSBwYUIuoVpS93dywAvrhg8tl4ekIjwS58pGmjoZ+dgwip60HY74RIPbjEQ0S9Yv1yX1Fai04hoAGQm67F4XrXIWR7jdF/A1TwaMEwp/AxJz8Dk0YM5KF/RCrDgEJEvWb9cn+r4jT+VvG123ASaBoADxcMVXwtTRvNYEKkMlziIaI++1vF1wE/V+dXE4a4fX3Z9GyGEKIgwoBCRB5x7HvyVoUh4OFEAvB25Rm394wafJVfxkJE3sElHiLqMdu+JxKA+36Wjk0H6wM9rG63LnNnDlHw4QwKEfWIY98TAaginHSHO3OIghNnUIjILevpw/+61B7wpRxPaACsnZeHnw5JYDghCkIMKETkkuOSTrCwzprcOUoX6KEQUS8xoBCRIqUlnWDw/J3X4Y5RaZw1IQpyrEEhIkVKrezVRILzrE6EJDGcEIUIBhQiUuRpK3t/0kjA9kUTUHwP29QThSou8RCRojRtNG6/IRU7ahsCPRQ71iCSm56A3PQEtqknClEMKERhzrpLJzM51u4L3mgyY+dxdYQTCUDR9GyMGnyVUxBhm3qi0MSAQhTGbHfpaCRAPzsHc/IzAHTVoAgV1KA8cdM1mDsugyGEKMywBoUoTDnu0rEIYEVprdzCPjM5NuBbiyMkieGEKEwxoBCFKaVdOp1C4HRTm/x7oCdQnr39WoYTojDl9YDyww8/4He/+x0yMzMRHR2NYcOG4Y9//CMsFot8jxACL7zwAnQ6HaKjozFlyhQcP37c20MhIjeUdunYnlnzyRf/G4BR2eMBf0Thy+sBpaSkBG+88QbWrVuHL774AqtXr8af/vQnrF27Vr5n9erVWLNmDdatW4eqqiqkpqbi1ltvxcWLF709HCJyYD2NGOiqObFu09VIwG+nDIOhqRVGkxmNLd8HcpiQAB7wRxTGJCG8WwY3Y8YMpKSk4M0335Sv3XPPPYiJicG7774LIQR0Oh0KCwuxbNkyAEB7eztSUlJQUlKCBQsWdPs3WlpaoNVqYTKZEB8f783hE4U0paLY79ouo3jnCaflnGnXp2DX54GbRVl+RzYWTMoK2N8nIu/z5Pvb6zMoBQUF+Oc//4lTp04BAI4cOYKKigrccccdAACDwYCGhgZMmzZNfk9UVBQmT56MyspKxc9sb29HS0uL3Q8R9Yx1xuRIfbNTUezy0mPQK4QTAD4PJ5KbCtzl0xlOiMKd17cZL1u2DCaTCdnZ2YiIiEBnZydeeukl3HfffQCAhoauvgopKSl270tJScGZM2cUP1Ov12PlypXeHipRyHM87M8xiASylb2rudt19+VhRi4P+SMKd16fQdm6dSvee+89bNq0CYcOHcI777yDf//3f8c777xjd5/k8K9PQgina1bLly+HyWSSf+rr6709bKKQc6S+GUXbguuwvwhJwpihCYEeBhGpgNdnUJYuXYqioiLMnTsXAJCTk4MzZ85Ar9dj/vz5SE1NBdA1k5KWlia/r7Gx0WlWxSoqKgpRUVHeHipRyNpadRbLth1TfE0jBXbmxJEGgAU8S4eI7Hk9oLS1tUGjsZ+YiYiIkLcZZ2ZmIjU1Fbt370ZeXh4AoKOjA+Xl5SgpKfH2cIhCkmN7etvfAaDIVTgBcF1qPI4b/V/HpbTEFCFJKF04Hm0dFp6lQ0R2vB5Q7rrrLrz00kvIyMjADTfcgMOHD2PNmjX49a9/DaBraaewsBCrVq3C8OHDMXz4cKxatQoxMTGYN2+et4dDFHIcd+LcnXc1yg6fl3+fm5/ucjnHAgQknADA2vvy0NrxA1aU1qJTCLtD/4iIHHl9m/HFixfx/PPPo6ysDI2NjdDpdLjvvvvw+9//Hv379wfQVW+ycuVKrF+/Hs3NzRg3bhxeffVVjBw5skd/g9uMKVwZTWZMLN6jqiWanoiQJFQUTZVne3j6MFF48uT72+sBxR8YUCjcWJdwvm3twOJNhwM9nB6x1rpYZ0qshxASUfjy5PubpxkTqZzjVmGlWo5Am3LtQOw7+Y1dseukEQM5U0JEvcaAQqRijicOCyDgJww7ktDVkRaAUyBhMCGi3mJAIVIxpROH1TZ7ItAVTMZnJTGQEJHXMKAQqZC15sTc8UOgh9It2xOQiYi8hQGFSGXsak58uJ7jjVoWNlcjIl9hQCFSEaeaEx+u51yvi8Pgq2LwcS8PBdRIQOnC8exjQkQ+wYBCpCJKNSe+cvzCRRy/cLHX77cIoK3D4sURERFdwYBCpALBVHNixdoTIvIlBhSiAHPsc6Jm1roV1p4Qka8xoBD5kdIhf459TtRs7X15SPpJFJuvEZHPMaAQ+YhjGHE85E8/OwfpiTGqPFfH1cnDY4YmMJgQkV8woBD5gGMYWXZ7Nko+OiGHEYsAVpTWonTheEiSb3fr9IZtx1ou6RBRIDCgEHmZ47KNRQAlO0/Acb9LpxCoOt2MRVOysO7TOr+PszsCwKvz8pAYyyUdIvI/BhQiL1PaKmwBFGdKXvzwC2gkYHT6Vaip/85fQ+yRCEnCT4dwSYeIAkMT6AEQhZrM5Fin3TiSBEwfmap4v0VAleGESzpEFEicQSHyAyGAj2obAj0Mt16cdQNyrtaircPCJR0iCjgGFCIvMzS1Km4XVuNuHat5P0vHzdelMJQQkWpwiYfIy5SWeNRu08F6TCzeg61VZwM9FCIiAAwoRH1mNJlRWdcEo8kc6KF4xDFEWbc+B9tzEFFo4hIPUR+4ar6m4tUczPtZOu7KvRr/am3H4k2H7V7rFAKnm9q41ENEAccZFCIXupsZUep3sqK0FuaOHyCpdI1HIwFLbh6O8VlJGDMkARqHcfIAQCJSCwYUIgVbq85iYvEezNt4wGVthlK/k04h8Og71arrDAt0bXXWz86RZ0fStNHQz85BxI9piluLiUhNuMRD5MDVzMikEQPtvrxj+0covl+F2QQAsH3hBOSmJ9hdm5OfgUkjBuJ0Uxu3FhORqjCgEDlwNTPiWJvR2tHp55H1TVuHY7P9LmnaaAYTIlIdBhQiB5nJsdBI9n1LrLUZticUK92nVqwtIaJgw4BC5MBam7GitBadQsi1GftOfSMv/UgScMfIVMzM1WF7zYVAD9mOBGD2T6/G9sMX7MbPWRIiCiaSEGos53OvpaUFWq0WJpMJ8fHxgR4OhRjrLEls/wi57TsATCzeo/rZkv/z82F4uGAo0rTRMJrMrC0hIlXx5PubMyhENpT6mozPSkJlXZPqw4kEyOEEYG0JEQU3bjMm+pHS7p3lpcdwpL4ZmcmxgR3cjyQAb84fo9hnRQA43dTm7yEREfkEAwrRj96qMDjNklgEMOvVSvx/H3wemEE5KJqejZuvS0XR9Gyn1yQJLIQlopDBJR4idM2ebNhvUHxNANhR2+DfATmQABTdkY0Fk7Jc36TyJSgiIk8woFBYst0unKaNxmenvw30kNzavuhKkzWjyYzinSec7rEu8bDuhIhCAQMKhR2lQtjYKPX+V0ECcKLhohxQDE2tiq30NeASDxGFDtagUFhRKoQt2nZM1TMoAl2t9q2HFrpqsf/bKVmcPSGikMGAQiHJ1UnESm3sBYC3K8/4b3C9YG21D7husV8wfKA/h0RE5FPqndcm6iWlJZw5+RkAXM8+qJ3GZoeOu1b8REShgjMoFFJcnURsnUk5+21w9gnJH5po14BNPzsHET82Q2EreyIKRT4JKOfPn8cDDzyApKQkxMTEYPTo0aiurpZfF0LghRdegE6nQ3R0NKZMmYLjx4/7YigUZtydRAwA/1P3rwCMqnsSgFmjdS5frzr9rd1y1Zz8DFQUTcXmx25ERdFUeYaIiChUeD2gNDc3Y+LEiejXrx927tyJzz//HH/+859x1VVXyfesXr0aa9aswbp161BVVYXU1FTceuutuHjxoreHQ2HGuvxhy/Yk4k0H6wMzsG4snJqFe/PTXb5uEc5dYtO00RiflcSZEyIKSV6vQSkpKUF6ejreeust+drQoUPl/yyEwCuvvILnnnsOs2fPBgC88847SElJwaZNm7BgwQJvD4nCiHX5w7YGxbr8UbjlUKCH59Lre+sQH93PqbbEijUmRBRuvD6D8v7772Ps2LH45S9/iUGDBiEvLw8bN26UXzcYDGhoaMC0adPka1FRUZg8eTIqKysVP7O9vR0tLS12P0TuWPuEWP/n+n112F5jDNyAHEg//lhZBLB650ksuz1bri2xYo0JEYUjr8+gfP3113j99dfx1FNPYcWKFTh48CAef/xxREVF4aGHHkJDQ1fL8JSUFLv3paSk4MwZ5a2eer0eK1eu9PZQKQQZTWYUbTsmd30XAJZtOxbIISkqGJ6M/V822V3rFAKjBl+FiqKpON3Uhpj+GrR1WDA0OYbhhIjCjtcDisViwdixY7Fq1SoAQF5eHo4fP47XX38dDz30kHyf5PBviUIIp2tWy5cvx1NPPSX/3tLSgvR01+v1FL7+o8IQFEfSVHzZ5HKrcJo2moGEiMKe15d40tLScP3119tdu+6663D27FkAQGpqKgDIMylWjY2NTrMqVlFRUYiPj7f7IXJkNJnxNxcH/qmNAPBowTBuFSYicsHrMygTJ07EyZMn7a6dOnUKQ4YMAQBkZmYiNTUVu3fvRl5eHgCgo6MD5eXlKCkp8fZwKAxYD/6ra7wUFLMnQNe/GTxcMBQPFwzF6aY2LuMQETnwekB58sknMWHCBKxatQr33nsvDh48iA0bNmDDhg0AupZ2CgsLsWrVKgwfPhzDhw/HqlWrEBMTg3nz5nl7OBTibLvGBgsJgP6eHLvGa0REZM/rASU/Px9lZWVYvnw5/vjHPyIzMxOvvPIK7r//fvmeZ599FmazGQsXLkRzczPGjRuHXbt2IS4uztvDoRDmWBAbDDQSULZwgnwyMRERKZOEUDq4Xd1aWlqg1WphMplYjxLGVpQeVW3jNXc2P3YjxmclBXoYRER+58n3N8/ioaC0vrxOleFEkq70N9HAvtcJwIZrREQ9xdOMKegcqW+GfueJQA/DiXX5ZlD8ALnwdd+pb7CitBadQnCnDhGRBxhQSHWsu3Iyk2Odvsy3Vp1FkQobrwFdPU3aOix2fUzm5Gdg0oiB3KlDROQhBhRSFWsAEehaHim6Ixs5V2uRmRyLxpbvUVSq3qJYV8s3bLxGROQ5BhRSDaU29fodXUs50o+/q5XtoYRERNR3DCikGp+d/tZlCFFzOAGAv87Nw4xcXaCHQUQUMriLh1TD1VlMaqeRgDFD2deEiMibGFDI74wmMyrrmmA0me2ujxmS4LQtV+0kCdDPzuHSDhGRl3GJh/zKtjW9BKBoejYWTM4C0FVMWnxPjt3rgP+Wd+4YmYqdxxvgrnVhhCShdOF41H9rhiQBPx2SwHBCROQDDCjkN0aT2e7cHAFAv/MEznzbiiU3DUeaNhrftV22O1en6I5s/OtiBzbs/9rn49tZ2+A2DFn7mOSmJ7BVPRGRjzGgkFe562FiaGpVPNRv04F6bDlYj9tHpmLHsQb5ugBQsuMEhJ/WfZTCyfN3XoexQxPQ1mFhHxMiIj9iQCGvsVu+kX5cvpmUJYeW881tLt9rEbALJ/J1wKtrPOMyE3DTdSko3nnC7VIO0DVjcseoNIYSIqIAYEAhr3BavhFdPUxqzn6Hj483KM6cBMJnp7/DK3PzcGNmIma9WmmXfSR0BSuLANvSExEFGAMKeYWr5Zudtc6zIoHUKQRON7VhfFYSiu/JcTonh23piYjUgQGFvCIzORaShG6XTQLNth29q3NyGEyIiAKPfVDIK9K00Sianh3oYbil1I4+TRuN8VlJDCVERCrDgEJeMzNXh5+ptKOqBKBs4QTMyc8I9FCIiKgHuMRDXrF+Xx2Kd5xQ5Zk5kgQUz85h7xIioiDCgEJ9tr68DvqdJwI9DEVP3HwN5v4sg0s4RERBhgGFes1oMuOz09+iWKXhJEKSGE6IiIIUa1CoRxwP+NtadRYTi/dgyeaabpd1NABenHWDXw8C1MC5IJaIiIIHZ1CoW7YdYjUSsOz2bJR8dKLHzdcsAL6/bPFbfYokdRXEsuaEiCh4cQaF3HLsEGsRQMlO1+FEaZYkQpKQPzQBGj9NoRRNz2Y4ISIKcgwo5JZSh1hLN++556dXI0LqSiO2JwDrZ+f4NKRoACz/8fwfIiIKblziIbcyk2Oh+fF8GisNgFuvT8HHn/+v0/0CwPbDF1C6cLziCcDePpNn0vBkPD1tBE8bJiIKMZxBIZljISzQ1WlVPztHnhEBumZQdimEE6tOIdDWYbHr0PrPLxqwbNsxr4/5v7/6FwbFD2A3WCKiEMMZFALgXAirn50jd12dk5+B7NQ4zHqtUj5rx91EiEaCfN6N0WRG4ZbDOGBo9sm4rYf/MZwQEYUWzqCQYiHsitJau5mU1o7OHh8EODe/q/fI1qqzmKDf47VwMi4z0akI1/bwPyIiCh0MKKRYCGudmbCy1qLYclXwuuTma+TQ05eSEwnAm/PHYMaoVADAAcO38nXgSgEuZ0+IiEIPl3hIuRBWAmL6X8mv1lqUFaW16BRCDgcAsHzbMXlnz6KpWWhs+R7/97P6PhfECgDJP4nCjmMNdtc0ErB2bh7GDE1gOCEiClEMKOQUPoCusHL3a5WKtSj//KIRA+OiMGnEQKRpo1F+8hvsqO0KEa9+WodXP63z2tj2nGh03uYsgKSfRDGcEBGFMAYUAtAVPuq/bcM6m3BhrUXJTo1Da0cnjp0z2R0K+Pzfj2Px1Cw5nPhC8k+inGZ3WHdCRBT6GFAIQFeh7Kt7nWc+OoXArFcrXdaSrPPibIkjCcAt16egf6TGaWmJsydERKGNASWEGU1mGJpakZkc2+0XuqGp1eUuHV+eoSNJUPy7GgD6e3KQpo3GnPwMTBoxEKeb2tiMjYgoTDCghCh3fU0A5/CiVCjrSxKAx36eiQ37DU6vPX/ndbhjVJpdEEnTRjOYEBGFEQaUINOTWRFXfU2sRa224UVC1+F6CyZnORXKepOErh0+EZKEgfFRuPm6FADA3yoMTvUljuGEiIjCj8/7oOj1ekiShMLCQvmaEAIvvPACdDodoqOjMWXKFBw/ftzXQwl6W6vOYmLxHszbeAATi/dga9VZxfvc9TVxDC8CgH7nCazfV4dJIwbilbm5eHHWDZC8fKifJAH33zgET067Fg/cOFSeEbFto8/6EiIisvLpDEpVVRU2bNiAUaNG2V1fvXo11qxZg7fffhsjRozAiy++iFtvvRUnT55EXFycL4cUtLqbFbF17LzJ6f0SutrPK4UXACjecQLFO07IfUZ+fk0y9n3Z5LXxWwQUW9KzvoSIiJT4bAbl0qVLuP/++7Fx40YkJCTI14UQeOWVV/Dcc89h9uzZGDlyJN555x20tbVh06ZNvhpO0OtJt1egK8gU7zgBVzKTY53axQNdMynWj7cIeDWcAO63Bqdpo3nYHxER2fFZQFm0aBHuvPNO3HLLLXbXDQYDGhoaMG3aNPlaVFQUJk+ejMrKSl8NJ+gptZpX+tKvPtOsuOtGANh8oGtJqGh6tm8G6YJGApduiIjIIz4JKFu2bMGhQ4eg1+udXmto6GrqlZKSYnc9JSVFfs1Re3s7Wlpa7H7CTU/qNbZWncXiTYddfsZf93yF8fo9+Neldvxq/BCfj1n+u3Pz7HYQERERdcfrNSj19fV44oknsGvXLgwYMMDlfZJDFaYQwumalV6vx8qVK706zmDkrl7DaDKjaNuxHn3Ohv0GxWUeX9AAGDM0odv7iIiIbHl9BqW6uhqNjY0YM2YMIiMjERkZifLycvz1r39FZGSkPHPiOFvS2NjoNKtitXz5cphMJvmnvr7e28MOGq7qNVwt7bjip3YnWDY9m0s7RETkMa8HlJtvvhnHjh1DTU2N/DN27Fjcf//9qKmpwbBhw5Camordu3fL7+no6EB5eTkmTJig+JlRUVGIj4+3+6ErjCYzTjYEftnrjpxUeQlKIwHL7+jqr0JEROQpry/xxMXFYeTIkXbXYmNjkZSUJF8vLCzEqlWrMHz4cAwfPhyrVq1CTEwM5s2b5+3hhLytVWdRtO2Yz2dEJLieddGga6ZkweQsGE1mbhkmIqI+C0gn2WeffRZmsxkLFy5Ec3Mzxo0bh127drEHioesdSeehBNXZ9+4owEwd1w6Nh1wXlp7/KZrcN+4DDmMsCU9ERF5gySED/qa+1hLSwu0Wi1MJlNYL/d8cOQ8lmyucbo+O0+H+ROG4kTDRadTgLNT45xOJ9YAmJ6Tig+POe+iun9cBhbfdA0AYIJ+j937JAmoLLqJgYSIiHrEk+9vn7e6J99xteup9PAF3P1aV0+ZiqKp2PzYjagomoo5+RnITU/AwqlZ8i6eCEnCsunZ2FnrHE4WTcnCS3fnyLMixffkyL1YNBJQPDuH4YSIiHyChwWqnLvDAccMcb1919oKv6JoKsZnJcnXn/6/Ndh26Lz8+203pCBnsFax/f3lTvuLbEtPRET+whkUFXN3OKDRZEb1mWa373dshX+kvtkunADAjtoGmDt+UOyL8reKr2E0me2usS09ERH5AwOKSrk6HNBoMsvBxV3XWKBr+SamvwaVdU0wmsw4ePpbxftON7XhsZ9nOl23HvBHRETkb1ziUSlXhwNWn262Cy6uSAAmXpOEu1+rhEV01YwsnKLck2Ts0AQMih+Av1UY7D7X3QF/REREvsQZFJVSOnU4QpIACW7DiQRg1NVaCHSdSGw7A/P63q8xfWSq3f33/PRq5KYn9OisHyIiIn/hDEqAuCt+BYB9p76x+10C8Oz0a5GeEO22aZoAcPS8SfG1TiHw0Pih+M3kYfjsdDPGDk1AbvqVQlsWwRIRkVowoATA1qqz8jKNRgL0s3PsTvs9Ut+MolL7BmwCgH7HiT79XeuSTZo22i6Y2GKjNSIiUgMu8fiZu+JXoCu8zHq10uNur93RSOCSDRERBQ3OoPiZq+JX626Z5aXePVdHAvDYpEw8PDGT4YSIiIIGA4qfZSbHQuNQ6GpdelEKL56YNy4dWw7Wy0tHjxYMw8MFQxlMiIgo6HCJx8/c7ZaxhhdHyg3tHe6RgM0HusKJBGDZ7dlYced1DCdERBSUeFhggBhNZsXdMlurzsoH/GkAPPrzTIwblohH3ql2+VmOMzJAV/CpKJrKgEJERKrhyfc3A4oKGU1mvFVxGn+r+LprRkRCr4pmNz92o905PERERIHE04xDgDWcAL0LJ+wCS0REwYwBxY+MJrN8Lo47fS2WZRdYIiIKdtzF4yfdNWezpbTTpyc0EvDXuXkYMzSB4YSIiIIaZ1D8oLvmbI7StNGY+zPl8OKKBsCjBZkMJ0REFBIYUPygu+ZsSsYPS3T5mgbApOHJ8pZk69k8G/YbMLF4D7ZWne3zmImIiAKJAcUPlPqbdFfEesH0vdM1CcCr8/JQtmgCKr66clKxwJXDA7ubnSEiIgoGDCheplQI6645m6vPKNnpfDBg0fRs3DlKh9aOTrf1Kd3NzhAREakdi2QdGE1mGJpakZkc63Etx/ryOhTvPAEB50LYOfkZmDRioGJzNkeudvGMGnwVgO6LaLnFmIiIgh0Dig1Pdto4Wr+vDnqbWQ/rUkt2ahxaOzrlwGMbTFyFIaUAogEQ079rwss6I2PtOCv9WIQiwC3GREQUGthJ9kdGkxkTi/c4HeLXk3bxRpMZE4r3KDZUsxawaiRg2fRs5FytRWZyLPad+gZF27pOLpYAFN9jH4ZsW95bOYYm23b5AHo0O0NERBQonnx/cwblR+522nT3hW9oanXZ7dW2eFW/o2uGRbK5br2naNsxTBoxUP5bc/IzkJ0ah1mvVcqfbZ2Vsd7nOCPDYEJERKGCRbI/6s1OG3fvdUcpywgA1aeb7a61dnQ6BR8WwBIRUThgQPmR404bjQQ8O/3aHs1KKL138dQsj0IL0HUooK2+hCYiIqJgxiUeG3PyM/Bd22UU7zwBiwBKdp7AVdH9elQoq7RLJz0xxqmOxBVJAn46JMHummMxLAtgiYgoXLBI1kZfCmXdfebppjYcPfcdVn90Ug4as/J0KDt8vkc7hmyLYRlOiIgoWLFItpf6UijrirWQdXxWEmaO1tkFjWduuxanm9oQ01+D1o5OGE1mxb/jWAxLREQU6hhQbCj1H/FmzYfSrpt9p77pde8VIiKiUMUiWRuetqTvKaX299brnpxyTEREFC44g+LAk5b03TGazHirwoCN+w2K7e99saREREQUChhQFHij5sO2bb6VY6M1Xy8pERERBSsu8fiA49KNLdtGa75aUiIiIgp2nEHxAVenEQPOMyTeXFIiIiIKFQwoPqC0dAN0TVcpzZBwGzEREZE9ry/x6PV65OfnIy4uDoMGDcKsWbNw8uRJu3uEEHjhhReg0+kQHR2NKVOm4Pjx494eSsAotb7/Pz8fhv9efhO3EBMREfWA12dQysvLsWjRIuTn5+OHH37Ac889h2nTpuHzzz9HbGwsAGD16tVYs2YN3n77bYwYMQIvvvgibr31Vpw8eRJxcXHeHlJAcOmGiIio93ze6v6bb77BoEGDUF5ejkmTJkEIAZ1Oh8LCQixbtgwA0N7ejpSUFJSUlGDBggXdfqavWt0TERGR73jy/e3zXTwmkwkAkJiYCAAwGAxoaGjAtGnT5HuioqIwefJkVFZWKn5Ge3s7Wlpa7H58zVVzNSIiIvI9nxbJCiHw1FNPoaCgACNHjgQANDQ0AABSUlLs7k1JScGZM2cUP0ev12PlypW+HKod2x4mbD9PRETkfz6dQVm8eDGOHj2KzZs3O70m/VhAaiWEcLpmtXz5cphMJvmnvr7eJ+MF2H6eiIhIDXw2g7JkyRK8//772LdvHwYPHixfT01NBdA1k5KWliZfb2xsdJpVsYqKikJUVJSvhmqH7eeJiIgCz+szKEIILF68GKWlpdizZw8yMzPtXs/MzERqaip2794tX+vo6EB5eTkmTJjg7eF4zNrDxBbbzxMREfmX1wPKokWL8N5772HTpk2Ii4tDQ0MDGhoaYDZ3LZFIkoTCwkKsWrUKZWVlqK2txa9+9SvExMRg3rx53h6Ox9h+noiIKPC8vs3YVR3JW2+9hV/96lcAumZZVq5cifXr16O5uRnjxo3Dq6++KhfSdscf24yNJjN7mBAREXmRJ9/fPu+D4gv+7INiNJlhaGpFZnIsgwoREVEfePL9zbN43OB2YyIiosDweaO2YMXtxkRERIHDgOKCu+3GRERE5FsMKC5wuzEREVHgMKC4wO3GREREgcMiWTfm5Gdg0oiB3G5MRETkZwwo3UjTRjOYEBER+RmXeIiIiEh1GFCIiIhIdRhQiIiISHUYUIiIiEh1GFCIiIhIdRhQiIiISHUYUIiIiEh1GFCIiIhIdRhQiIiISHUYUIiIiEh1GFCIiIhIdYLyLB4hBACgpaUlwCMhIiKinrJ+b1u/x90JyoBy8eJFAEB6enqAR0JERESeunjxIrRardt7JNGTGKMyFosFFy5cQFxcHCRJ8uvfbmlpQXp6Ourr6xEfH+/Xv60GfH4+P5+fzx+Ozx/Ozw547/mFELh48SJ0Oh00GvdVJkE5g6LRaDB48OCAjiE+Pj4s/5/Uis/P5+fz8/nDUTg/O+Cd5+9u5sSKRbJERESkOgwoREREpDoMKB6KiorCH/7wB0RFRQV6KAHB5+fz8/n5/OH4/OH87EBgnj8oi2SJiIgotHEGhYiIiFSHAYWIiIhUhwGFiIiIVIcBhYiIiFSHAcUFvV6P/Px8xMXFYdCgQZg1axZOnjxpd48QAi+88AJ0Oh2io6MxZcoUHD9+PEAj9h29Xg9JklBYWChfC/VnP3/+PB544AEkJSUhJiYGo0ePRnV1tfx6KD//Dz/8gN/97nfIzMxEdHQ0hg0bhj/+8Y+wWCzyPaH0/Pv27cNdd90FnU4HSZKwfft2u9d78qzt7e1YsmQJkpOTERsbi5kzZ+LcuXN+fIrec/f8ly9fxrJly5CTk4PY2FjodDo89NBDuHDhgt1nhOrzO1qwYAEkScIrr7xidz3Un/+LL77AzJkzodVqERcXhxtvvBFnz56VX/fV8zOguFBeXo5Fixbh//2//4fdu3fjhx9+wLRp09Da2irfs3r1aqxZswbr1q1DVVUVUlNTceutt8pnBYWCqqoqbNiwAaNGjbK7HsrP3tzcjIkTJ6Jfv37YuXMnPv/8c/z5z3/GVVddJd8Tys9fUlKCN954A+vWrcMXX3yB1atX409/+hPWrl0r3xNKz9/a2orc3FysW7dO8fWePGthYSHKysqwZcsWVFRU4NKlS5gxYwY6Ozv99Ri95u7529racOjQITz//PM4dOgQSktLcerUKcycOdPuvlB9flvbt2/HgQMHoNPpnF4L5eevq6tDQUEBsrOzsXfvXhw5cgTPP/88BgwYIN/js+cX1CONjY0CgCgvLxdCCGGxWERqaqooLi6W7/n++++FVqsVb7zxRqCG6VUXL14Uw4cPF7t37xaTJ08WTzzxhBAi9J992bJloqCgwOXrof78d955p/j1r39td2327NnigQceEEKE9vMDEGVlZfLvPXnW7777TvTr109s2bJFvuf8+fNCo9GIjz76yG9j9wbH51dy8OBBAUCcOXNGCBEez3/u3Dlx9dVXi9raWjFkyBDx8ssvy6+F+vPPmTNH/u++El8+P2dQeshkMgEAEhMTAQAGgwENDQ2YNm2afE9UVBQmT56MysrKgIzR2xYtWoQ777wTt9xyi931UH/2999/H2PHjsUvf/lLDBo0CHl5edi4caP8eqg/f0FBAf75z3/i1KlTAIAjR46goqICd9xxB4DQf35bPXnW6upqXL582e4enU6HkSNHhtz/PoCufxZKkiTPKIb681ssFjz44INYunQpbrjhBqfXQ/n5LRYLPvzwQ4wYMQK33XYbBg0ahHHjxtktA/ny+RlQekAIgaeeegoFBQUYOXIkAKChoQEAkJKSYndvSkqK/Fow27JlCw4dOgS9Xu/0Wqg/+9dff43XX38dw4cPx8cff4zf/OY3ePzxx/Gf//mfAEL/+ZctW4b77rsP2dnZ6NevH/Ly8lBYWIj77rsPQOg/v62ePGtDQwP69++PhIQEl/eEiu+//x5FRUWYN2+efGBcqD9/SUkJIiMj8fjjjyu+HsrP39jYiEuXLqG4uBi33347du3ahbvvvhuzZ89GeXk5AN8+f1CeZuxvixcvxtGjR1FRUeH0miRJdr8LIZyuBZv6+no88cQT2LVrl906o6NQfHag698axo4di1WrVgEA8vLycPz4cbz++ut46KGH5PtC9fm3bt2K9957D5s2bcINN9yAmpoaFBYWQqfTYf78+fJ9ofr8SnrzrKH2v4/Lly9j7ty5sFgseO2117q9PxSev7q6Gn/5y19w6NAhj58lFJ7fWhj/i1/8Ak8++SQAYPTo0aisrMQbb7yByZMnu3yvN56fMyjdWLJkCd5//318+umnGDx4sHw9NTUVAJwSYmNjo9O/bQWb6upqNDY2YsyYMYiMjERkZCTKy8vx17/+FZGRkfLzheKzA0BaWhquv/56u2vXXXedXLUeyv+3B4ClS5eiqKgIc+fORU5ODh588EE8+eST8mxaqD+/rZ48a2pqKjo6OtDc3OzynmB3+fJl3HvvvTAYDNi9e7c8ewKE9vPv378fjY2NyMjIkP9ZeObMGTz99NMYOnQogNB+/uTkZERGRnb7z0NfPT8DigtCCCxevBilpaXYs2cPMjMz7V7PzMxEamoqdu/eLV/r6OhAeXk5JkyY4O/hetXNN9+MY8eOoaamRv4ZO3Ys7r//ftTU1GDYsGEh++wAMHHiRKct5adOncKQIUMAhPb/7YGunRsajf0/GiIiIuR/mwr157fVk2cdM2YM+vXrZ3eP0WhEbW1tSPzvwxpOvvzyS3zyySdISkqyez2Un//BBx/E0aNH7f5ZqNPpsHTpUnz88ccAQvv5+/fvj/z8fLf/PPTp8/epxDaE/fa3vxVarVbs3btXGI1G+aetrU2+p7i4WGi1WlFaWiqOHTsm7rvvPpGWliZaWloCOHLfsN3FI0RoP/vBgwdFZGSkeOmll8SXX34p/uu//kvExMSI9957T74nlJ9//vz54uqrrxb/+Mc/hMFgEKWlpSI5OVk8++yz8j2h9PwXL14Uhw8fFocPHxYAxJo1a8Thw4flXSo9edbf/OY3YvDgweKTTz4Rhw4dEjfddJPIzc0VP/zwQ6Aeq8fcPf/ly5fFzJkzxeDBg0VNTY3dPwvb29vlzwjV51fiuItHiNB+/tLSUtGvXz+xYcMG8eWXX4q1a9eKiIgIsX//fvkzfPX8DCguAFD8eeutt+R7LBaL+MMf/iBSU1NFVFSUmDRpkjh27FjgBu1DjgEl1J/9gw8+ECNHjhRRUVEiOztbbNiwwe71UH7+lpYW8cQTT4iMjAwxYMAAMWzYMPHcc8/ZfSGF0vN/+umniv9dnz9/vhCiZ89qNpvF4sWLRWJiooiOjhYzZswQZ8+eDcDTeM7d8xsMBpf/LPz000/lzwjV51eiFFBC/fnffPNNcc0114gBAwaI3NxcsX37drvP8NXzS0II0bc5GCIiIiLvYg0KERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpDgMKERERqQ4DChEREakOAwoRERGpzv8PmbdL1KcRBpsAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "importlib.reload(tfr_inference)\n", + "\n", + "# Get some parameters from the data\n", + "sigma_m, sigma_eta, hyper_eta_mu, hyper_eta_sigma = tfr_inference.estimate_data_parameters()\n", + "\n", + "# Other parameters to use\n", + "L = 500.0\n", + "N = 64\n", + "xmin = -L/2\n", + "R_lim = L / 2\n", + "Rmax = 100\n", + "Nt = 2000\n", + "alpha = 1.4\n", + "mthresh = 11.25\n", + "a_TFR = -23\n", + "b_TFR = -8.2\n", + "sigma_TFR = 0.3\n", + "sigma_v = 150\n", + "Nint_points = 201 \n", + "Nsig = 10\n", + "frac_sigma_r = 0.07 # WANT A BETTER WAY OF DOING THIS - ESTIMATE THROUGH SIGMAS FROM TFR\n", + "interp_order = 1\n", + "bias_epsilon = 1.e-7\n", + "\n", + "cpar, dens, vel = tfr_inference.get_fields(L, N, xmin)\n", + "\n", + "RA, Dec, czCMB, m_true, eta_true, m_obs, eta_obs, xtrue = tfr_inference.create_mock(\n", + " Nt, L, xmin, cpar, dens, vel, Rmax, alpha, mthresh,\n", + " a_TFR, b_TFR, sigma_TFR, sigma_m, sigma_eta, \n", + " hyper_eta_mu, hyper_eta_sigma, sigma_v, \n", + " interp_order=interp_order, bias_epsilon=bias_epsilon)\n", + "\n", + "\n", + "plt.figure()\n", + "rtrue = jnp.sqrt(jnp.sum(xtrue ** 2, axis=0))\n", + "robs = czCMB / 100 \n", + "plt.plot(rtrue, robs, '.')\n", + " \n", + "MB_pos = tfr_inference.generateMBData(RA, Dec, czCMB, L, N, R_lim, Nsig, Nint_points, sigma_v, frac_sigma_r)\n", + "r = jnp.sqrt(jnp.sum(MB_pos ** 2, axis=0))\n", + "\n", + "rtrue = jnp.sqrt(jnp.sum(xtrue ** 2, axis=0))\n", + "\n", + "for i in range(3):\n", + " print(r[i,0], r[i,-1], rtrue[i])\n", + " \n", + "# Check truth always inside the range of the MB pos\n", + "dr_low = rtrue - r[:,0]\n", + "dr_high = r[:,-1] - rtrue\n", + "m = dr_low < 0\n", + "print(m.sum())\n", + "m = dr_high < 0\n", + "print(m.sum())\n", + "print('Low:', dr_low.min())\n", + "print('High', dr_high.min())\n", + "print(r.min(), r.max())\n", + "\n", + "czCMB_new = ((1 + zcosmo) * (1 + vr_noised / tfr_inference.utils.speed_of_light) - 1) * tfr_inference. utils.speed_of_light\n", + "\n", + "plt.figure()" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "3f0bd23a-9040-4d90-9abf-bdc9fa69b4b0", + "id": "67759f3d-0cb3-4771-bb07-f3d5bfc6848c", "metadata": {}, "outputs": [], "source": [] diff --git a/tests/corner.png b/tests/corner.png new file mode 100644 index 0000000..0bb7ee4 Binary files /dev/null and b/tests/corner.png differ diff --git a/tests/likelihood_scan_a_TFR.png b/tests/likelihood_scan_a_TFR.png new file mode 100644 index 0000000..47ef585 Binary files /dev/null and b/tests/likelihood_scan_a_TFR.png differ diff --git a/tests/likelihood_scan_alpha.png b/tests/likelihood_scan_alpha.png new file mode 100644 index 0000000..cf9f26d Binary files /dev/null and b/tests/likelihood_scan_alpha.png differ diff --git a/tests/likelihood_scan_b_TFR.png b/tests/likelihood_scan_b_TFR.png new file mode 100644 index 0000000..75cdd4b Binary files /dev/null and b/tests/likelihood_scan_b_TFR.png differ diff --git a/tests/likelihood_scan_sigma_TFR.png b/tests/likelihood_scan_sigma_TFR.png new file mode 100644 index 0000000..8eda8fd Binary files /dev/null and b/tests/likelihood_scan_sigma_TFR.png differ diff --git a/tests/likelihood_scan_sigma_v.png b/tests/likelihood_scan_sigma_v.png new file mode 100644 index 0000000..d73ddba Binary files /dev/null and b/tests/likelihood_scan_sigma_v.png differ diff --git a/tests/test_ll.png b/tests/test_ll.png new file mode 100644 index 0000000..560c2ad Binary files /dev/null and b/tests/test_ll.png differ diff --git a/tests/tfr_inference.py b/tests/tfr_inference.py index 3258587..aaa16ff 100644 --- a/tests/tfr_inference.py +++ b/tests/tfr_inference.py @@ -3,18 +3,20 @@ import numpy as np from astropy.coordinates import SkyCoord import astropy.units as apu import astropy.constants -from astropy.cosmology import LambdaCDM import pandas as pd -from scipy.interpolate import interp1d import jax.numpy as jnp import jax.scipy.special +import matplotlib.pyplot as plt +import corner + import numpyro import numpyro.distributions as dist -from jax import lax +from jax import lax, random import borg_velocity.poisson_process as poisson_process import borg_velocity.projection as projection +import borg_velocity.utils as utils # Output stream management cons = borg.console() @@ -166,22 +168,16 @@ def create_mock(Nt, L, xmin, cpar, dens, vel, Rmax, alpha, mthresh, # Initialize lists to store valid positions and corresponding sig_mu values all_xtrue = np.empty((3, Nt)) + all_mtrue = np.empty(Nt) + all_etatrue = np.empty(Nt) all_mobs = np.empty(Nt) all_etaobs = np.empty(Nt) + all_RA = np.empty(Nt) + all_Dec = np.empty(Nt) # Counter for accepted positions accepted_count = 0 - - # Cosmology object needed for z <-> r - cosmo = LambdaCDM(H0 = cpar.h * 100, Om0 = cpar.omega_m, Ode0 = cpar.omega_q) - - # Precompute redshift-distance mapping - z_grid = np.logspace(-4, -1, 10000) # Covers z ~ 0 to 0.1 - dL_grid = cosmo.luminosity_distance(z_grid).value # Luminosity distances in Mpc - - # Create an interpolation function: distance -> redshift - dL_to_z = interp1d(dL_grid, z_grid, kind="cubic", fill_value="extrapolate") - + # Bias model phi = (1. + dens + bias_epsilon) ** alpha @@ -207,12 +203,11 @@ def create_mock(Nt, L, xmin, cpar, dens, vel, Rmax, alpha, mthresh, r_hat = np.array(SkyCoord(ra=RA*apu.deg, dec=Dec*apu.deg).cartesian.xyz) # Compute cosmological redshift - # zcosmo = z_at_value(cosmo.comoving_distance, rtrue * apu.Mpc / cpar.h).value - zcosmo = dL_to_z(rtrue / cpar.h) + zcosmo = utils.z_cos(rtrue, cpar.omega_m) # Compute luminosity distance # DO I NEED TO DO /h??? - dL = (1 + zcosmo) * rtrue / cpar.h # Mpc/h + dL = (1 + zcosmo) * rtrue / cpar.h # Mpc # Compute true distance modulus mutrue = 5 * np.log10(dL) + 25 @@ -232,18 +227,28 @@ def create_mock(Nt, L, xmin, cpar, dens, vel, Rmax, alpha, mthresh, # Apply apparement magnitude cut m = mobs <= mthresh + mtrue = mtrue[m] + etatrue = etatrue[m] mobs = mobs[m] etaobs = etaobs[m] xtrue = xtrue[:,m] + RA = RA[m] + Dec = Dec[m] # Calculate how many valid positions we need to reach Nt remaining_needed = Nt - accepted_count selected_count = min(xtrue.shape[1], remaining_needed) # Append only the needed number of valid positions - all_xtrue[:,accepted_count:accepted_count+selected_count] = xtrue[:,:selected_count] - all_mobs = mobs[:selected_count] - all_etaobs = etaobs[:selected_count] + imin = accepted_count + imax = accepted_count + selected_count + all_xtrue[:,imin:imax] = xtrue[:,:selected_count] + all_mtrue[imin:imax] = mtrue[:selected_count] + all_etatrue[imin:imax] = etatrue[:selected_count] + all_mobs[imin:imax] = mobs[:selected_count] + all_etaobs[imin:imax] = etaobs[:selected_count] + all_RA[imin:imax] = RA[:selected_count] + all_Dec[imin:imax] = Dec[:selected_count] # Update the accepted count accepted_count += selected_count @@ -266,11 +271,15 @@ def create_mock(Nt, L, xmin, cpar, dens, vel, Rmax, alpha, mthresh, np.zeros(3,) )) # km/s + # Recompute cosmological redshift + rtrue = jnp.sqrt(jnp.sum(all_xtrue ** 2, axis=0)) + zcosmo = utils.z_cos(rtrue, cpar.omega_m) + # Obtain total redshift vr_noised = vr_true + sigma_v * np.random.randn(Nt) - zCMB = (1 + zcosmo) * (1 + vr_noised / astropy.constants.c.to('km/s').value) - 1 + czCMB = ((1 + zcosmo) * (1 + vr_noised / utils.speed_of_light) - 1) * utils.speed_of_light - return zCMB, all_mobs, all_etaobs, all_xtrue + return all_RA, all_Dec, czCMB, all_mtrue, all_etatrue, all_mobs, all_etaobs, all_xtrue def estimate_data_parameters(): @@ -305,66 +314,414 @@ def estimate_data_parameters(): return sigma_m, sigma_eta, hyper_eta_mu, hyper_eta_sigma - -def likelihood(a_TFR, b_TFR, sigma_TFR, eta_true, m_true): +def generateMBData(RA, Dec, cz_obs, L, N, R_lim, Nsig, Nint_points, sigma_v, frac_sigma_r): + """ + Generate points along the line of sight of each tracer to enable marginalisation + over distance uncertainty. The central distance is given such that the observed + redshift equals the cosmological redshift at this distance. The range is then + +/- Nsig * sig, where + sig^2 = (sig_v/100)^2 + sig_r^2 + and sig_v is the velocity uncertainty in km/s - loglike = 0 + Args: + - RA (np.ndarray): Right Ascension (degrees) of the tracers (shape = (Nt,)) + - Dec (np.ndarray): Delination (degrees) of the tracers (shape = (Nt,)) + - cz_obs (np.ndarray): Observed redshifts (km/s) of the tracers (shape = (Nt,)) + - L (float): Box length (Mpc/h) + - N (int): Number of grid cells per side + - R_lim (float): Maximum allowed (true) comoving distance of a tracer (Mpc/h) + - Nsig (float): ??? + - Nint_points (int): Number of radii over which to integrate the likelihood + - sigma_v (float): Uncertainty on the velocity field (km/s) + - frac_sigma_r (float): An estimate of the fractional uncertainty on the positions of tracers + + Returns: + - MB_pos (np.ndarray): Comoving coordinates of integration points to use in likelihood (Mpc/h). + The shape is (3, Nt, Nsig) + + """ + + myprint(f"Making MB data") + + # Convert RA, DEC to radial vector + r_hat = np.array(SkyCoord(ra=RA*apu.deg, dec=Dec*apu.deg).cartesian.xyz).T + + # Get min and max distance to integrate over + # cz = 100 h r, so sigma_v corresponds to a sigma_r of ~ sigma_v / 100 + robs = cz_obs / 100 + sigr = np.sqrt((sigma_v / 100) ** 2 + (frac_sigma_r * robs)**2) + rmin = robs - Nsig * sigr + rmin = rmin.at[rmin <= 0].set(L / N / 100.) + rmax = robs + Nsig * sigr + rmax = rmax.at[rmax > R_lim].set(R_lim) + + # Compute coordinates of integration points + r_integration = np.linspace(rmin, rmax, Nint_points) + MB_pos = np.expand_dims(r_integration, axis=2) * r_hat[None,...] + MB_pos = jnp.transpose(MB_pos, (2, 1, 0)) + + return MB_pos + + +def likelihood(alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, m_true, eta_true, + dens, vel, omega_m, h, L, xmin, interp_order, bias_epsilon, + cz_obs, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh): + """ + Evaluate the likelihood for TFR sample + + Args: + - alpha (float): Exponent for bias model + - a_TFR (float): TFR relation intercept + - b_TFR (float): TFR relation slope + - sigma_TFR (float): Intrinsic scatter in the TFR + - sigma_v (float): Uncertainty on the velocity field (km/s) + - m_true (np.ndarray): True apparent magnitudes of the tracers (shape = (Nt,)) + - eta_true (np.ndarray): True linewidths of the tracers (shape = (Nt,)) + - dens (np.ndarray): Over-density field (shape = (N, N, N)) + - vel (np.ndarray): Velocity field (km/s) (shape = (3, N, N, N)) + - omega_m (float): Matter density parameter Om + - h (float): Hubble constant H0 = 100 h km/s/Mpc + - L (float): Comoving box size (Mpc/h) + - xmin (float): Coordinate of corner of the box (Mpc/h) + - interp_order (int): Order of interpolation from grid points to the line of sight + - bias_epsilon (float): Small number to add to 1 + delta to prevent 0^# + - cz_obs (np.ndarray): Observed redshifts (km/s) of the tracers (shape = (Nt,)) + - m_obs (np.ndarray): Observed apparent magnitudes of the tracers (shape = (Nt,)) + - eta_obs (np.ndarray): Observed linewidths of the tracers (shape = (Nt,)) + - sigma_m (float): Uncertainty on the apparent magnitude measurements + - sigma_eta (float): Uncertainty on the apparent linewidth measurements + - MB_pos (np.ndarray): Comoving coordinates of integration points to use in likelihood (Mpc/h). + The shape is (3, Nt, Nsig) + - mthresh (float): Threshold absolute magnitude in selection + + Returns: + - loglike (float): The log-likelihood of the data + """ + + # Comoving radii of integration points (Mpc/h) + r = jnp.sqrt(jnp.sum(MB_pos ** 2, axis=0)) + + # p_r = r^2 n(r) N(mutrue; muTFR, sigmaTFR) + # Multiply by arbitrary number for numerical stability (cancels in p_r / p_r_norm) + number_density = projection.interp_field( + dens, + MB_pos, + L, + jnp.array([xmin, xmin, xmin]), + interp_order, + use_jitted=True, + ) + number_density = jax.nn.relu(1. + number_density) + number_density = jnp.power(number_density + bias_epsilon, alpha) + zcosmo = utils.z_cos(r, omega_m) + mutrue = 5 * jnp.log10((1 + zcosmo) * r / h) + 25 + muTFR = m_true - (a_TFR + b_TFR * eta_true) + d2 = ((mutrue - muTFR[:,None]) / sigma_TFR) ** 2 + best = jnp.amin(jnp.abs(d2), axis=1) + d2 = d2 - jnp.expand_dims(jnp.nanmin(d2, axis=1), axis=1) + p_r = r ** 2 * jnp.exp(-0.5 * d2) * number_density + p_r_norm = jnp.expand_dims(jnp.trapezoid(p_r, r, axis=1), axis=1) + + # Peculiar velocity term + tracer_vel = projection.interp_field( + vel, + MB_pos, + L, + jnp.array([xmin, xmin, xmin]), + interp_order, + use_jitted=True, + ) + tracer_vr = projection.project_radial( + tracer_vel, + MB_pos, + jnp.zeros(3,) + ) + cz_pred = ((1 + zcosmo) * (1 + tracer_vr / utils.speed_of_light) - 1) * utils.speed_of_light + d2 = ((cz_pred - jnp.expand_dims(cz_obs, axis=1)) / sigma_v)**2 + scale = jnp.nanmin(d2, axis=1) + d2 = d2 - jnp.expand_dims(scale, axis=1) + + # Integrate to get likelihood + p_cz = jnp.trapezoid(jnp.exp(-0.5 * d2) * p_r / p_r_norm, r, axis=1) + lkl_ind = jnp.log(p_cz) - scale / 2 - 0.5 * jnp.log(2 * np.pi * sigma_v**2) + loglike_vel = - lkl_ind.sum() + + Nt = m_obs.shape[0] # Apparent magnitude terms norm = 0.5 * (1 + jax.scipy.special.erf((mthresh - m_true) / (jnp.sqrt(2) * sigma_m))) - loglike += - 0.5 * jnp.sum((mobs - m_true) ** 2 / sigma_m ** 2) + loglike_m = ( + 0.5 * jnp.sum((m_obs - m_true) ** 2 / sigma_m ** 2) + jnp.sum(jnp.log(norm)) + Nt * 0.5 * jnp.log(2 * jnp.pi * sigma_m ** 2) + ) # Linewidth terms - loglike += - 0.5 * jnp.sum((etaobs - eta_true) ** 2 / sigma_eta ** 2) + loglike_eta = ( + 0.5 * jnp.sum((eta_obs - eta_true) ** 2 / sigma_eta ** 2) + Nt * 0.5 * jnp.log(2 * jnp.pi * sigma_eta ** 2) + ) + + # loglike = - (loglike_vel + loglike_m + loglike_eta) + loglike = - (loglike_eta + loglike_m) + + return loglike + + +def test_likelihood_scan(prior, alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, m_true, eta_true, + dens, vel, omega_m, h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh): + """ + Plot likelihood as we scan through the paramaters [alpha, a_TFR, b_TFR, sigma_TFR, sigma_v] + to verify that the likelihood shape looks reasonable + + Args: + - prior (dict): Upper and lower bounds for a uniform prior for the parameters + - alpha (float): Exponent for bias model + - a_TFR (float): TFR relation intercept + - b_TFR (float): TFR relation slope + - sigma_TFR (float): Intrinsic scatter in the TFR + - sigma_v (float): Uncertainty on the velocity field (km/s) + - m_true (np.ndarray): True apparent magnitudes of the tracers (shape = (Nt,)) + - eta_true (np.ndarray): True linewidths of the tracers (shape = (Nt,)) + - dens (np.ndarray): Over-density field (shape = (N, N, N)) + - vel (np.ndarray): Velocity field (km/s) (shape = (3, N, N, N)) + - omega_m (float): Matter density parameter Om + - h (float): Hubble constant H0 = 100 h km/s/Mpc + - L (float): Comoving box size (Mpc/h) + - xmin (float): Coordinate of corner of the box (Mpc/h) + - interp_order (int): Order of interpolation from grid points to the line of sight + - bias_epsilon (float): Small number to add to 1 + delta to prevent 0^# + - cz_obs (np.ndarray): Observed redshifts (km/s) of the tracers (shape = (Nt,)) + - m_obs (np.ndarray): Observed apparent magnitudes of the tracers (shape = (Nt,)) + - eta_obs (np.ndarray): Observed linewidths of the tracers (shape = (Nt,)) + - sigma_m (float): Uncertainty on the apparent magnitude measurements + - sigma_eta (float): Uncertainty on the apparent linewidth measurements + - MB_pos (np.ndarray): Comoving coordinates of integration points to use in likelihood (Mpc/h). + The shape is (3, Nt, Nsig) + - mthresh (float): Threshold absolute magnitude in selection - # Peculiar velocity terms + """ + + + pars = [alpha, a_TFR, b_TFR, sigma_TFR, sigma_v] + par_names = ['alpha', 'a_TFR', 'b_TFR', 'sigma_TFR', 'sigma_v'] + + orig_ll = likelihood(*pars, m_true, eta_true, + dens, vel, omega_m, h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh) + + for i, name in enumerate(par_names): + + myprint(f'Scanning {name}') + + if name in prior: + x = np.linspace(*prior[name], 20) + else: + pmin = pars[i] * 0.2 + pmax = pars[i] * 2.0 + x = np.linspace(pmin, pmax, 20) + + all_ll = np.empty(x.shape) + orig_x = pars[i] + for j, xx in enumerate(x): + pars[i] = xx + all_ll[j] = likelihood(*pars, m_true, eta_true, + dens, vel, omega_m, h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh) + pars[i] = orig_x + + plt.figure() + plt.plot(x, all_ll, '.') + plt.axvline(orig_x, ls='--', color='k') + plt.axhline(orig_ll, ls='--', color='k') + plt.xlabel(name) + plt.ylabel('Negative log-likelihood') + plt.savefig(f'likelihood_scan_{name}.png') + fig = plt.gcf() + plt.clf() + plt.close(fig) return -def likelihood_model(): +def run_mcmc(num_warmup, num_samples, prior, initial, dens, vel, omega_m, h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh, + m_true): + """ + Run MCMC over the model parameters - # TO DO: Sort out these priors - a_TFR = numpyro.sample("a_TFR", dist.Uniform(*alpha_prior)) - b_TFR = numpyro.sample("b_TFR", dist.Uniform(*alpha_prior)) - sigma_TFR = numpyro.sample("sigma_TFR", dist.Uniform(*alpha_prior)) - - eta_true = numpyro.sample("eta_true", dist.Normal(mean, sigma), sample_shape=(Nt,)) - m_true = numpyro.sample("m_true", dist.Normal(mean, sigma), sample_shape=(Nt,)) - - numpyro.sample("obs", TFRLikelihood(a_TFR, b_TFR, sigma_TFR, eta_true, m_true)) - - return - - -class TFRLikelihood(dist.Distribution): - support = dist.constraints.real - - def __init__(self, a_TFR, b_TFR, sigma_TFR, eta_true, m_true): - self.a_TFR, self.b_TFR, self.sigma_TFR, self.eta_true, self.m_true = dist.util.promote_shapes(a_TFR, b_TFR, sigma_TFR, eta_true, m_true) - batch_shape = lax.broadcast_shapes( - jnp.shape(a_TFR), - jnp.shape(b_TFR), - jnp.shape(sigma_TFR), - jnp.shape(eta_true), - jnp.shape(m_true), - ) - super(TFRLikelihood, self).__init__(batch_shape = batch_shape) + Args: + - num_warmup (int): Number of warmup steps to take in the MCMC + - num_samples (int): Number of samples to take in the MCMC + - prior + - initial + - dens (np.ndarray): Over-density field (shape = (N, N, N)) + - vel (np.ndarray): Velocity field (km/s) (shape = (3, N, N, N)) + - omega_m (float): Matter density parameter Om + - h (float): Hubble constant H0 = 100 h km/s/Mpc + - L (float): Comoving box size (Mpc/h) + - xmin (float): Coordinate of corner of the box (Mpc/h) + - interp_order (int): Order of interpolation from grid points to the line of sight + - bias_epsilon (float): Small number to add to 1 + delta to prevent 0^# + - cz_obs (np.ndarray): Observed redshifts (km/s) of the tracers (shape = (Nt,)) + - m_obs (np.ndarray): Observed apparent magnitudes of the tracers (shape = (Nt,)) + - eta_obs (np.ndarray): Observed linewidths of the tracers (shape = (Nt,)) + - sigma_m (float): Uncertainty on the apparent magnitude measurements + - sigma_eta (float): Uncertainty on the apparent linewidth measurements + - MB_pos (np.ndarray): Comoving coordinates of integration points to use in likelihood (Mpc/h). + The shape is (3, Nt, Nsig) + - mthresh (float): Threshold absolute magnitude in selection - def sample(self, key, sample_shape=()): - raise NotImplementedError + """ - def log_prov(self, value) - return likelihood(self.a_TFR, self.b_TFR, self.sigma_TFR, self.eta_true, self.m_true) + Nt = eta_obs.shape[0] + + def tfr_model(): + + alpha = numpyro.sample("alpha", dist.Uniform(*prior['alpha'])) + a_TFR = numpyro.sample("a_TFR", dist.Uniform(*prior['a_TFR'])) + b_TFR = numpyro.sample("b_TFR", dist.Uniform(*prior['b_TFR'])) + sigma_TFR = numpyro.sample("sigma_TFR", dist.HalfCauchy(1.0)) + sigma_v = numpyro.sample("sigma_v", dist.HalfCauchy(1.0)) + +# # Sample the means with a uniform prior +# hyper_mean_m = numpyro.sample("hyper_mean_m", dist.Uniform(*prior['hyper_mean_m'])) +# hyper_mean_eta = numpyro.sample("hyper_mean_eta", dist.Uniform(*prior['hyper_mean_eta'])) +# hyper_mean = jnp.array([hyper_mean_m, hyper_mean_eta]) + +# # Sample standard deviations with a 1/sigma prior (Jeffreys prior approximation) +# hyper_sigma_m = numpyro.sample("hyper_sigma_m", dist.HalfCauchy(1.0)) # Equivalent to 1/sigma prior +# hyper_sigma_eta = numpyro.sample("hyper_sigma_eta", dist.HalfCauchy(1.0)) +# hyper_sigma = jnp.array([hyper_sigma_m, hyper_sigma_eta]) + +# # Sample correlation matrix using LKJ prior +# L_corr = numpyro.sample("L_corr", dist.LKJCholesky(2, concentration=1.0)) # Cholesky factor of correlation matrix +# corr_matrix = L_corr @ L_corr.T # Convert to full correlation matrix + +# # Construct full covariance matrix: Σ = D * Corr * D +# hyper_cov = jnp.diag(hyper_sigma) @ corr_matrix @ jnp.diag(hyper_sigma) + +# # Sample the true eta and m +# x = numpyro.sample("x", dist.MultivariateNormal(hyper_mean, hyper_cov), sample_shape=(Nt,)) +# m_true = numpyro.deterministic("m_true", x[:, 0]) +# eta_true = numpyro.deterministic("eta_true", x[:, 1]) + + hyper_mean_eta = numpyro.sample("hyper_mean_eta", dist.Uniform(*prior['hyper_mean_eta'])) + hyper_sigma_eta = numpyro.sample("hyper_sigma_eta", dist.HalfCauchy(1.0)) # Equivalent to 1/sigma prior + eta_true = numpyro.sample("eta_true", dist.Normal(hyper_mean_eta, hyper_sigma_eta), sample_shape=(Nt,)) + + # Evaluate the likelihood + numpyro.sample("obs", TFRLikelihood(alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, eta_true), obs=jnp.array([m_obs, eta_obs])) + + + class TFRLikelihood(dist.Distribution): + support = dist.constraints.real + + def __init__(self, alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, eta_true): + self.alpha, self.a_TFR, self.b_TFR, self.sigma_TFR, self.sigma_v, self.eta_true = dist.util.promote_shapes(alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, eta_true) + batch_shape = lax.broadcast_shapes( + jnp.shape(alpha), + jnp.shape(a_TFR), + jnp.shape(b_TFR), + jnp.shape(sigma_TFR), + jnp.shape(sigma_v), + # jnp.shape(m_true), + jnp.shape(eta_true), + ) + super(TFRLikelihood, self).__init__(batch_shape = batch_shape) + + def sample(self, key, sample_shape=()): + raise NotImplementedError + + def log_prob(self, value): + loglike = likelihood(self.alpha, self.a_TFR, self.b_TFR, self.sigma_TFR, self.sigma_v, + m_true, self.eta_true, + dens, vel, omega_m, h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh) + return loglike + + rng_key = random.PRNGKey(6) + rng_key, rng_key_ = random.split(rng_key) + values = initial + myprint('Preparing MCMC kernel') + kernel = numpyro.infer.NUTS(tfr_model, + init_strategy=numpyro.infer.initialization.init_to_value(values=initial) + ) + mcmc = numpyro.infer.MCMC(kernel, num_warmup=num_warmup, num_samples=num_samples) + myprint('Running MCMC') + mcmc.run(rng_key_) + mcmc.print_summary() + + return mcmc + + +def process_mcmc_run(mcmc, param_labels, truths, obs): + + # Convert samples into a single array + samples = mcmc.get_samples() + + samps = jnp.empty((len(samples[param_labels[0]]), len(param_labels))) + for i, p in enumerate(param_labels): + samps = samps.at[:,i].set(samples[p]) + + # Trace plot of non-redshift quantities + fig1, axs1 = plt.subplots(samps.shape[1], 1, figsize=(6,3*samps.shape[1]), sharex=True) + axs1 = np.atleast_1d(axs1) + for i in range(samps.shape[1]): + axs1[i].plot(samps[:,i]) + axs1[i].set_ylabel(param_labels[i]) + axs1[i].axhline(truths[i], color='k') + axs1[-1].set_xlabel('Step Number') + fig1.tight_layout() + fig1.savefig('trace.png') + + # Corner plot + fig2, axs2 = plt.subplots(samps.shape[1], samps.shape[1], figsize=(10,10)) + corner.corner( + np.array(samps), + labels=param_labels, + fig=fig2, + truths=truths + ) + fig2.savefig('corner.png') + + # True vs predicted + for var in ['eta', 'm']: + vname = var + '_true' + if vname in samples.keys(): + xtrue = obs[var] + xpred_median = np.median(samples[vname], axis=0) + xpred_plus = np.percentile(samples[vname], 84, axis=0) - xpred_median + xpred_minus = xpred_median - np.percentile(samples[vname], 16, axis=0) + + fig3, axs3 = plt.subplots(2, 1, figsize=(10,8), sharex=True) + plot_kwargs = {'fmt':'.', 'markersize':3, 'zorder':10, + 'capsize':1, 'elinewidth':1, 'alpha':1} + axs3[0].errorbar(xtrue, xpred_median, yerr=[xpred_minus, xpred_plus], **plot_kwargs) + axs3[1].errorbar(xtrue, xpred_median - xtrue, yerr=[xpred_minus, xpred_plus], **plot_kwargs) + axs3[1].set_xlabel('True') + axs3[0].set_ylabel('True') + axs3[1].set_ylabel('True - Predicted') + xlim = axs3[0].get_xlim() + ylim = axs3[0].get_ylim() + axs3[0].plot(xlim, xlim, color='k', zorder=0) + axs3[0].set_xlim(xlim) + axs3[0].set_ylim(ylim) + axs3[1].axhline(0, color='k', zorder=0) + fig3.suptitle(var) + fig3.align_labels() + fig3.tight_layout() + fig3.savefig(f'true_predicted_{var}.png') + + + return def main(): + myprint('Beginning') + # Get some parameters from the data sigma_m, sigma_eta, hyper_eta_mu, hyper_eta_sigma = estimate_data_parameters() @@ -372,22 +729,83 @@ def main(): L = 500.0 N = 64 xmin = -L/2 + R_lim = L / 2 Rmax = 100 - Nt = 2000 + Nt = 100 alpha = 1.4 mthresh = 11.25 a_TFR = -23 b_TFR = -8.2 sigma_TFR = 0.3 sigma_v = 150 + Nint_points = 201 + Nsig = 10 + frac_sigma_r = 0.07 # WANT A BETTER WAY OF DOING THIS - ESTIMATE THROUGH SIGMAS FROM TFR + interp_order = 1 + bias_epsilon = 1.e-7 + num_warmup = 1000 + num_samples = 1000 + prior = { + 'alpha': [0.5, 2.5], + 'a_TFR': [-25, -20], + 'b_TFR': [-10, -5], + 'hyper_mean_eta': [hyper_eta_mu - 0.5, hyper_eta_mu + 0.5], + # 'hyper_mean_m':[mthresh - 5, mthresh + 5] + } + initial = { + 'alpha': alpha, + 'a_TFR': a_TFR, + 'b_TFR': b_TFR, + 'hyper_mean_eta': hyper_eta_mu, + 'hyper_sigma_eta': hyper_eta_sigma, + # 'hyper_mean_m': mthresh, + 'sigma_TFR': sigma_TFR, + 'sigma_v': sigma_v, + } + + # Make mock + np.random.seed(123) cpar, dens, vel = get_fields(L, N, xmin) - - zCMB, mobs, etaobs, xtrue = create_mock( + RA, Dec, czCMB, m_true, eta_true, m_obs, eta_obs, xtrue = create_mock( Nt, L, xmin, cpar, dens, vel, Rmax, alpha, mthresh, a_TFR, b_TFR, sigma_TFR, sigma_m, sigma_eta, - hyper_eta_mu, hyper_eta_sigma, sigma_v) + hyper_eta_mu, hyper_eta_sigma, sigma_v, + interp_order=interp_order, bias_epsilon=bias_epsilon) + MB_pos = generateMBData(RA, Dec, czCMB, L, N, R_lim, Nsig, Nint_points, sigma_v, frac_sigma_r) + # Test likelihood + loglike = likelihood(alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, m_true, eta_true, + dens, vel, cpar.omega_m, cpar.h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh) + myprint(f'loglike {loglike}') + + # Scan over parameters to make plots verifying behaviour + test_likelihood_scan(prior, alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, m_true, eta_true, + dens, vel, cpar.omega_m, cpar.h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh) + + # Run a MCMC + mcmc = run_mcmc(num_warmup, num_samples, prior, initial, dens, vel, cpar.omega_m, cpar.h, L, xmin, interp_order, bias_epsilon, + czCMB, m_obs, eta_obs, sigma_m, sigma_eta, MB_pos, mthresh, + m_true) + + param_labels = ['alpha', 'a_TFR', 'b_TFR', 'sigma_TFR', 'sigma_v', 'hyper_mean_eta', 'hyper_sigma_eta'] + truths = [alpha, a_TFR, b_TFR, sigma_TFR, sigma_v, hyper_eta_mu, hyper_eta_sigma] + param_labels = ['hyper_mean_eta', 'hyper_sigma_eta'] + truths = [hyper_eta_mu, hyper_eta_sigma] + obs = {'m':m_obs, 'eta':eta_obs} + process_mcmc_run(mcmc, param_labels, truths, obs) if __name__ == "__main__": - main() \ No newline at end of file + main() + +""" +TO DO + +- Fails to initialise currently when loglike includes the BORG term +- Runs MCMC with this likelihood +- Add bulk velocity +- Deal with case where sigma_eta and sigma_m could be floats vs arrays + +""" \ No newline at end of file diff --git a/tests/trace.png b/tests/trace.png new file mode 100644 index 0000000..4a29ba6 Binary files /dev/null and b/tests/trace.png differ diff --git a/tests/true_predicted_eta.png b/tests/true_predicted_eta.png new file mode 100644 index 0000000..2bcd5da Binary files /dev/null and b/tests/true_predicted_eta.png differ