notebooks/main_analyses/joint_LDstats.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:38.908240Z",
     "iopub.status.busy": "2022-01-02T14:38:38.905937Z",
     "iopub.status.idle": "2022-01-02T14:38:41.569460Z",
     "shell.execute_reply": "2022-01-02T14:38:41.568646Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "from sympy import *\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "from scipy.stats import binned_statistic, describe\n",
    "\n",
    "import sys\n",
    "sys.path.append('../../src/')\n",
    "from aDNA_coal_sim import *\n",
    "from plot_utils import *\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:41.581922Z",
     "iopub.status.busy": "2022-01-02T14:38:41.581295Z",
     "iopub.status.idle": "2022-01-02T14:38:41.595129Z",
     "shell.execute_reply": "2022-01-02T14:38:41.595743Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.rcParams['font.sans-serif'] = \"Arial\"\n",
    "plt.rcParams['figure.facecolor'] = \"w\"\n",
    "plt.rcParams['figure.autolayout'] = True\n",
    "plt.rcParams['pdf.fonttype'] = 3\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "import os\n",
    "main_figdir = '../../plots/jointLD_stats/'\n",
    "supp_figdir = '../../plots/supp_figs/jointLD_stats/'\n",
    "os.makedirs(main_figdir, exist_ok=True)\n",
    "os.makedirs(supp_figdir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deriving the joint LD statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:41.603064Z",
     "iopub.status.busy": "2022-01-02T14:38:41.602366Z",
     "iopub.status.idle": "2022-01-02T14:38:41.754312Z",
     "shell.execute_reply": "2022-01-02T14:38:41.754928Z"
    }
   },
   "outputs": [],
   "source": [
    "# Define variables in sympy... \n",
    "rho,t,eta,gamma = symbols('rho t \\eta \\gamma')\n",
    "\n",
    "gamma = rho*(1  - exp(-t *(rho/2 + 1)))/(rho + 2)\n",
    "eta = 2*(1  - exp(-t *(rho/2 + 1)))/(rho + 2)\n",
    "\n",
    "eT_AT_B_200 = (rho**2 + 14*rho + 36)/(rho**2 + 13*rho + 18)\n",
    "eT_AT_B_111 = (rho**2 + 13*rho + 24)/(rho**2 + 13*rho + 18)\n",
    "eT_AT_B_022 = (rho**2 + 13*rho + 22)/(rho**2 + 13*rho + 18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:41.856471Z",
     "iopub.status.busy": "2022-01-02T14:38:41.813892Z",
     "iopub.status.idle": "2022-01-02T14:38:43.272220Z",
     "shell.execute_reply": "2022-01-02T14:38:43.272778Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\left(\\rho^{3} e^{\\frac{t \\left(\\rho + 2\\right)}{2}} + 15 \\rho^{2} e^{\\frac{t \\left(\\rho + 2\\right)}{2}} + 48 \\rho e^{\\frac{t \\left(\\rho + 2\\right)}{2}} + 48 e^{\\frac{t \\left(\\rho + 2\\right)}{2}} - 4\\right) e^{- \\frac{t \\left(\\rho + 2\\right)}{2}}}{\\left(\\rho + 2\\right) \\left(\\rho^{2} + 13 \\rho + 18\\right)}$"
      ],
      "text/plain": [
       "(rho**3*exp(t*(rho + 2)/2) + 15*rho**2*exp(t*(rho + 2)/2) + 48*rho*exp(t*(rho + 2)/2) + 48*exp(t*(rho + 2)/2) - 4)*exp(-t*(rho + 2)/2)/((rho + 2)*(rho**2 + 13*rho + 18))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Defining the staggered statistics (from the note)\n",
    "e_TATB_200_anc = (1 - gamma)*eT_AT_B_200 + gamma*eT_AT_B_111\n",
    "e_TATB_022_anc = (1 - eta).simplify() * eT_AT_B_022 + eta*eT_AT_B_111\n",
    "e_TATB_111_mod_coupled = (1-gamma)*eT_AT_B_111 + gamma*eT_AT_B_022\n",
    "e_TATB_111_mod_uncoupled = (1-eta)*eT_AT_B_111 + eta*eT_AT_B_200\n",
    "\n",
    "# Simplifying the statistics\n",
    "e_TATB_200_anc = e_TATB_200_anc.factor().simplify()\n",
    "e_TATB_022_anc = e_TATB_022_anc.factor().simplify()\n",
    "e_TATB_111_mod_coupled = e_TATB_111_mod_coupled.factor().simplify()\n",
    "e_TATB_111_mod_uncoupled = e_TATB_111_mod_uncoupled.factor().simplify()\n",
    "\n",
    "\n",
    "e_TATB_022_anc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:43.334669Z",
     "iopub.status.busy": "2022-01-02T14:38:43.331831Z",
     "iopub.status.idle": "2022-01-02T14:38:43.371015Z",
     "shell.execute_reply": "2022-01-02T14:38:43.371831Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\left(\\rho + 10\\right) e^{- \\frac{t \\left(\\rho + 2\\right)}{2}}}{\\rho^{2} + 13 \\rho + 18}$"
      ],
      "text/plain": [
       "(rho + 10)*exp(-t*(rho + 2)/2)/(rho**2 + 13*rho + 18)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eD0DT_numerator = (e_TATB_200_anc - e_TATB_111_mod_coupled - e_TATB_111_mod_uncoupled + e_TATB_022_anc)\n",
    "\n",
    "e_D0DT_numerator = eD0DT_numerator.factor().simplify()\n",
    "e_D0DT_numerator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:43.392178Z",
     "iopub.status.busy": "2022-01-02T14:38:43.391468Z",
     "iopub.status.idle": "2022-01-02T14:38:43.395098Z",
     "shell.execute_reply": "2022-01-02T14:38:43.394429Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\left(\\rho + 2\\right) \\left(\\rho + 10\\right)}{\\rho^{3} e^{\\frac{t \\left(\\rho + 2\\right)}{2}} + 15 \\rho^{2} e^{\\frac{t \\left(\\rho + 2\\right)}{2}} + 48 \\rho e^{\\frac{t \\left(\\rho + 2\\right)}{2}} + 48 e^{\\frac{t \\left(\\rho + 2\\right)}{2}} - 4}$"
      ],
      "text/plain": [
       "(rho + 2)*(rho + 10)/(rho**3*exp(t*(rho + 2)/2) + 15*rho**2*exp(t*(rho + 2)/2) + 48*rho*exp(t*(rho + 2)/2) + 48*exp(t*(rho + 2)/2) - 4)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_D0Dt_norm = e_D0DT_numerator / e_TATB_022_anc\n",
    "e_D0Dt_norm # Note this does not actually go \"negative\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:43.421431Z",
     "iopub.status.busy": "2022-01-02T14:38:43.420480Z",
     "iopub.status.idle": "2022-01-02T14:38:44.981449Z",
     "shell.execute_reply": "2022-01-02T14:38:44.982387Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00,  6.40it/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAADQCAYAAADcQn7hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA200lEQVR4nO3dd1yVZf/A8c99NiAyxIUoKoii4ESFcOVKQ0kzZ/ZolrvSsmFP/Up7srInW2YmWtpTmjscuS1NE9yCirlwgYMhCDLPuH9/HEUJhANyOCDX+/W6Xwfucd3fc4Av172+lyTLsowgCIIVKGwdgCAIjy6RYARBsBqRYARBsBqRYARBsBqRYARBsJoqn2D69Olj6xAE4ZFV5RNMUlKSrUMQhEdWlU8wgiBYj0gwgiBYjUgwgiBYjcrWAQhCedLr9cTFxZGdnW3rUColnU6Hh4cHarXaovVFgilGQno2/xd+AoUkIUkgYX5VSBIKCSTp3vdKSUKhMM9XKiTzJEkolRIqhYRKoUClkFCrFKiVCtRKCbVSgUapQKtWoFUp0aoU6NRK7NRK7DQK7DQqHDRKHLQq1ErR4XxYcXFxODo60rBhQyRJsnU4lYosyyQnJxMXF0ejRo0s2kYkmGJkZmZyKeYgJiRkJGRARsKEAhPSnfkKTLKEEUXefCMKjCjvvCow3PkaSv9LrVUpcNSpqG6nxtlOjbO9Bmc7NW6OWmpW01LT0Ty5O9vh7qxDq1KW2efwqMjOzhbJpZQkSaJGjRokJiZavI1IMMWoJSezRTu9zNozSkpMkgqjpMaICoOkwiCpyZU05KImFw3ZqMmUNWTIWjJNGm6ZtKQYtKTLdmRk6biV6UAq1UiUq3EWB5JkJ7LQ5duPJEEtRy0eLvZ41XTAu1Y1mtRypEntatRztqvSf2BV+b0/rJJ+diLBFMNep4NazUGWQTYBd17vfn//ZDLe+dpo/tpkvPO1wTzJJpSyEaVsRE1OyQIp5ieVo7DnltKFZJy5YqzBOb0rlzPciLtdkz2X3FmJK3d7TzUcNLSq70zr+s6083ShnacLOrXo7QhlTySY4jjXh0kRZdOW6U7yMeaCUX9nygVjDhhywZANhhzQZ4I+685rJuRmQE465KSZX7NvQVYKZN40v95OQGvMpJYpk1rE4wv0VgL35YwchR1XlPU5qXfnYHYjok57MffvBuhRoVMr6NioBl19atK9WS0aujmUzfsVqjyRYMqTQgEoQGnZGXiLybI5+dxOgPRrcCsOUi9D6hVIuQBJZ9BmJOJtOoM3Z3hKvQsAg6ThjNKLrVnN+fNsSz4848UHGxW0qu/M023q0a9lXWpU05ZtrEKVIlX1inYBAQEcOnTI1mFYX0YyJJ2G6ycg/rB5Sj6bfxWFI7uMLVmtf4w9Jn9QqOnpW5vxXRvTpoGLjQIvW6dOncLX19fWYRAXF8dff/3F0KFDy6S9LVu2MGXKFIxGIy+++CLTpxd+3vBB61m6PZTsMxQJpqokmMJkpcClCDi/E87tNPd27khTOPOrviMrDF2JkRsS1LgGE7t50bmJW6U+SVpREsyPP/5ITEwMs2fPfui2jEYjPj4+bN++HQ8PD9q3b88vv/xC8+bNLVqvadOmFm1/V0k+wzK/sUKpFCcLKw07F2j2JITMgSnH4OUj0P3/wM2H6qZURim3skn7b5ZqZ8PFP/nXD/t5ev4+jsfdsnXkldrevXt57bXXWL16Na1btyY2Nvah2jtw4ADe3t40btwYjUbDsGHDWLduncXrWbp9aZT5OZgq3iGq3Gp4QZfXofM0uHYMjv0CR38mWB9FsCaKk3jz37iBhM5LZWRHT17v3RQn+zI+n1SOGk7/zSrtXvwkpMjlnTp1on379nz22Wf4+fkVuW7nzp1JT08vMP+zzz6jZ8+eAMTHx1O/fv28ZR4eHuzfv7/ANg9az9LtS6PME8w/u8+ZmZnFbqNQKNDpdMWuJ5QTSQL3Nuap23Q4uAj2f0eLzHMs0fyXHaa2zNz/HL8dv8aM0BaEtnK3dcSVzunTp2nWrFmx6+3Zs6ccorEeq19FcnV1xd3dvdCejSRJyLKM0Wjk8uXL1g5FKA17V+j6JgS9BIe+h12z6Zl7hC7a48zLDuWNX/pz4II3/9eveaW7c7i4noa1JCUl4eTkhEpl/vObMWMGKSkp1KhRg/feey/fupb0YOrVq8eVK1fylsXFxVGvXr0C2zxoPUu3Lw2rJxhfX1+OHj1a5Dpt2rSxdhjCw9LYw2Mvg/9g2PZ/aI6v5FX1GkJU+5mwfyqD424xb0Rb6rva2zrSCu/ixYu4u5t7ffHx8RgMBpydnfnrr78KrGtJD6Z9+/acPXuWCxcuUK9ePZYvX86yZcssXq9p06YWbV8aVk8whw8fLpN1hArCsQ4MWgjtRsGGKfgkn2Oj9l3evDqWfnMz+WZEGzo3qWnrKCu0Zs2akZSUhJ+fH7Vq1WLdunUkJibm60WUhEql4ptvvuGJJ57AaDQyZswYWrRokbf8ySefZNGiRbi7uz9wvaK2fxhlfplaqVRiNBrLskmrKu4ydXpuOhvOb0AhKe48SW2eFChQSOaLcApJkTcpJSWSJKGUlOZJkf9VrVCjUqjyvWqUGvOk0KBValEpVJXjUnBOOqx/BU6uBWCJoTf/lZ/jy2c70qt5bRsHV7iKcpn6rjlz5pCWlkZycjJNmzbl5ZdftnVIxSrJZ2jVHszmzZtZs2YNKSkpaDQaAJo0acKUKVOoUaNGqdo0mUxMmjSJqKgotFotixYtwtvbu8A6ISEhPPXUU0yYMOGh3kNqdiofH/j4odooKYWkQKfUoVPpsFPZYaeyw0HtQDV1NezV9lTXVKe6tjpOGiectE646Fxws3PDzc6NGroa6FTldMJc6wjP/ACejyFveZvRqm14G+OZ+PPrfDQ0kP7i5G+xpk2bZusQrMqqCUav1xMWFoZCce92m7i4OC5dulTqBBMeHk52djYRERFERkYybdq0Atfs3333XVJSUh4q9rscNA4MbWq+29Ikm5CRkWUZGRmTbDLPk2VMmDCZTOZX2YTRZMQo35nufG0wGTDKRvQmPQaTAYPJgN6kJ9eYi96kJ8eYQ44xB4PJQKYhk0xD8VfgCuOidaFutbq4O7jjXs2dhk4NaezUmMZOjXHRlfEduZIEHcYiubdF/mUYnTJO8r30CS8uf4McQxDPtPMo2/0JlYpVE0xoaCjffPMNL730EgAXLlygUaNGeHiU/pdu7969eUONBAYGFji8Wb16NQqFosyGI3HVufJu4Ltl0palDCYDOcYcsgxZZBmyyNRnkqHPyJvSctNIy03jVs4tUnNSuZl9k+SsZJKzk0nKSiIlJ4WUnBRikmMKfT/NazTHz80Pfzd//Nz8cNW5PnzQHu2Qnt+M/L9QOqSd5if1R/xr1XQkghkkkkyVZbUEk5qayrRp0zhz5gw6nQ5/f3/CwsL4/vvvAfPZ8xs3btC2bdsStZuWloaTk1Pe90qlEoPBgEql4sSJEyxbtozVq1fzwQcfPLCNsLAwwsLCAEpUPKe8qBQqVAoVDuqSP9Vskk0kZyVzNeMqV29fJf52PBduXSA2NZYLaRe4mX2TvfF72Ru/N2+bZq7NCHIP4jH3x2hTqw1aZSkfcHTzRnp+M/zYn1apsSzXfMjotf+mvusTdGhUBklMqHSsepI3NzeX7du34+bmRnR0NHXr1qVfv37k5uYSEhJCcHAw58+fJzU1leDg4CIfsLrrtddeIzAwkCFDhgDmuw7j4uIAePPNN9m9ezd2dnZcvHgRjUbD119/XWRvpio9iyTLMlczrnIi6QQnk04SnRTNiaQT5Bjv1aaxV9nT07MnIY1C6FC3AypFKf4HpV2FH0Mh+SxHTd5MVM5gxUvd8axh+zIQFe0kb2VUYU7yvv322yiVSpo3b46bmxtnz55l7NixjBw5El9fX2bMmMH48eP56aefmDx5skVtBgcHs2HDBoYMGUJkZCT+/v55yz799NO8r2fMmEGdOnXEyI33kSSJetXqUa9aPZ5o+AQAOcYcjtw4QsS1CPbF7+N0ymnWn1/P+vPrqaGrQX+v/jzr+yx1HOpYvqPq7jD6N+RF3Wlz6xzvGL7hhSWOrJkUjJNd5X20QCg5q1+mvnnzJhEREaSlpdGkSRPatWuHJEkEBATQtWtXTp48yRdffMG8efP45ptvim3/7lWk6OhoZFlm8eLFbNq0CW9vb0JDQ/PWu5tgiruKVJV6MJa4lHaJTbGb+O3Cb1xKuwSASlLRt1FfRvuNxsfFx/LGbpxE/r43Uu5tvjQ8zeFGE1g8uj0qGxYvFz2Yh2fTcg0luQ8mNjaWo0ePcuTIEaKioti4cWNZhmIRkWAKJ8syx5OO83PMz2y9tBWTbAKge/3uvN7+deo71i+mhTvObEP+ZSiSbOKV3Ml4dhvFtN5NrRh50USCeXgVMsEcOXKEOXPmkJKSgr+/PxMnTqRhw4ZluetSEQmmeHHpcfwU8xO/nvuVLEMWGoWG0X6jedH/RexUdsU3EPkdbHmLHFnNEP37vD9hJG1tVMBKJJiHZ9N6MA8ydOhQQkJCmDVrFj4+PgwYMIBt27aV1+6Fh+Dh6MHbHd9m48CN9Gvcj1xTLmHRYYSGh7L7yu7iG+g4Hto9j1bS87lqHm8v309GjsH6gQs2V24JxsXFhREjRtCmTRteeOEFdu7cyRtvvFFeuxfKQC37Wnzc+WP+1/d/+Lr6cj3jOi/9/hKfHfwMvUn/4A0lCfp8gsmtKV6Ka4xIW8SsTafKL/AKKC4ujhUrVpRZe1u2bKFp06Z4e3vzySeflHi9MWPGUKtWrWLr05RUuSUYLy8v5syZk1e2wdnZubx2LZSxNrXa8EvIL7zW7jWUkpIfY35k9JbRXL199cEbqXUoBi1EVqgZpdpO/MEN/P73jfILuoLZuXMnR44cKZO2jEYjkydPZvPmzcTExPDLL78QE1PwJsui1hs9ejRbtmwpk3juV24JJicnh/nz59OgQQP69OmDn58fPXv2JD4+vrxCEMqQUqHkeb/nWdJnCXUc6hCdGM3gDYPZF7/vwRvVbYX0+L8B+FS9gFmr9pGSkVtOEVccFa1kJkCXLl1wdS37myHLbdiStWvNT9xmZGRw/PhxoqOjiY6OZvjw4Vy9epVz586VVyhCGWpdqzWr+q3inb/e4c+4P5n8+2Rmd55N74a9C98geAryma3UvhLJa7nf8tlWL2Y93bJ8g75rhlPx65Sq3aJrFle0kpnWVO7jIjk4OBAYGEhgYGB571qwEmedM3O7z2XOoTn8L+Z/vPHnG2ToMxjYZGDBlRVKpIHfYZofTAgH2Hh4DSc6euJXz0p/7BWUKJkpCCWgkBS8HvA61TTV+PbYt7y37z0y9BmMbD6y4MqujVD0mgmbXucd1c+8tb4bP0/oVv41cIrpaVhLRSuZaU3lnmCuX79OnToluO1cqDQkSWJiq4k4qh2ZfXA2sw+ax/wpNMkEjMF4aAkeCScIiPuJ9VE+PNXaur/sFUVFK5lpTeV+z3bv3g84NhceGSObj2TmYzMB+PTgp/x++feCKymUKJ80Pzs2UbWexb/9SWZu1bg35v6Smc899xxvvfUWo0aNKnVv4v6Smb6+vgwZMqRAycyrV68Wud7w4cMJCgri9OnTeHh45FU9eFjl/qiAv78/x48fL8tdPhRxJ6/1hEWHMffoXOxUdix+YjEt3ArWeZVXjUE6uYbfjB041ekbXn/Cuo8RVLQ7eUXJzDJWKWrNCmVirP9YLqddZt35dbz0+0sse3IZdavVzbeO1PsDjKc3EcIBVuxZT1yHyXi4VJ2RCR71kpm2e6xVeORJksT7Qe/ToU4HkrKSmLRzErdzb+dfyckDZRfzH9k7iiXM//1vG0QqWItIMIJVqZVqPu/2OY2cGnEu9Ryz9s8quFLQy+irN6CpIg790RVcTi5dLWKh4in3BKNUVq7R/4SH56R14qvHv8JOZcfG2I1sit2UfwW1DnWPdwB4SbGGeTur9nNKj5JyTzDFjfIoPJoaOTXi9YDXAfgw8sOCzy35D0bv7EUDRSKK6F+4kJRhgyiFsiYOkYRyM9hnMN3qdyNdn86/9/4bo+m+q40KJeoe5ueUJit/5dsdBR/WEyofkWCEciNJEjMfm4mbnRuHbxxm8cnF+VdoMRC9a1M8pCS0x5dxLuF24Q0JlYZNE4y1H7QSKh5XnSv/Cf4PAPOOzuP0zdP3Ft7Xi5mkCufb7SdsEaJQhmyaYAYPHmzL3Qs20qleJ4Y2HYpBNvDR/o/Id6+nbyh6t+a4SzdxjFkmzsVUcla/0e7u+EX/JMsyN2/etPbuhQrqlbavsP3Sdo4kHGFj7Eb6e/U3L1AoUPd8F5aPYJJqHfP//Bczng6wbbBCqVm9B7Njxw5GjRrF5MmTC0wODrYfiEuwjeqa6rza7lUAPj/8ef4b8Jo+SY5bC2pLqeiPriD5ds4DWqm8bFUy80GlMS3dvqSsnmC6deuGo6MjXbt2zTd169aNli1tVGhIqBBCvUJpVbMVSVlJfBv17b0FkoS2y1QAnpc28lPEBdsEaEW2KJkJhZfGLMn2JWX1BLN27Vq6dOlS6LLt27eXuD2TycSECRMICgqiW7duBSrhzZs3j/bt29OhQwdWrlxZqpiF8qGQFLzT8R0UkoJlp5ZxLuW+n2WLgeTY18VbcZULEb+SrbdsrK3KwFYlM6Hw0pgl2b6kKl3BqfDwcLKzs4mIiCAyMpJp06blfRhJSUnMnz+fo0ePkp2dTfPmzRk8eLB4wLIC863hy2Cfwaw4vYKPDnzE972/N/+8lGo0nSbDtncZYQhn9eERjAz0LNN9+//oX/xKpXB8VNHVAmxVMvNBrFlK0+o9mF69erFo0SISExPzzTeZTERERDBp0iSWLFlicXt79+7NG286MDAwX6kFNzc3jh07hlqt5vr16+h0ukKTS1hYGAEBAQQEBBSISyh/L7d5GRetCwevH2R33L1xlqR2o9GrHemo+Ju/dm/FaCrTyiI2VZKSmceOHSsw3U0uFV2xPZjMzOIfPFMoFOh0ukKXhYeHs2jRIvr160dCQgIuLi5kZWWRnZ1N165dmThxIh07drQ44LS0NJyc7tVvVSqVGAyGvPKDd4vqvP/++7zyyiuFtjFu3DjGjRsHmOvBCLblpHViXMtxzD44m7lH59LFowsKSQFaR5QBz0PE14TcXs32mFD6+JVdNcTiehrWYquSmQ9izVKaxSYYV1dX3N3dKawulSRJyLKM0Wjk8uXLhW7v4ODAlClTmDJlCrm5uSQnJ6PT6XBxKd3QodWrV8/3gZtMprwf1F0vvfQS48aNo2/fvvzxxx88/vjjpdqXUH4GNx3MkpNLOJNyhm2XttGnobmXqgiaiDHyW/oqDvDyH3/Rx2+QjSN9eLYqmWmt7YtS7CGSr68vsbGxXLhwocB0d36NGjUs2plGo6Fu3bqlTi4AwcHBbNpkfho3MjISf/97x9GnT5/m6aefRpZl1Go1Wq0WhUI8DVEZaJVaxrcaD5jv8DWY7pTPrO6Oye8ZlJJMhxvLOR5nm0LdZclWJTOh8NKYxW3/MIotmWkymYr9I71/nftLZm7evJk1a9Zw8+ZNNBoNAD4+PkyZMsXipFTYviZNmkR0dDSyLLN48WI2bdqEt7c3oaGhzJw5k82bNyNJEn379i3Q5fwnUTKz4tCb9IT+Gkrc7Tg+DP6Qp7yfMi+4cRLmP0amrOVj37X8Z1inUu9DlMx8eCX5DK1ak3f9+vX0798/34nWuLg4EhISaNu2bVnuttREgqlYNpzfwL/3/pt61eqxYcAG1Eo1AFmLQrCL28ss47+YNH0OLg6aUrVf0RJMZVRmNXk3b97M2rVrSU5ORqPRIEkSTZo0sbgHkpGRwfDhw1GpVJhMJkJDQxk2bBgeHh6WvROhynmy0ZMsPL6QC7cu8Ou5XxnS1PyoiV3wBFixl+HSNlYevMT4bk1sHKlgiSITjF6vp2XLlnndttjYWDQaDZcuXbIowezatYvly5fnfT958mSGDRv2kCELjzKlQsnk1pN5fffrLIhawADvAWiUGvDpS7Z9XRpnXuP0vvUYu7yGUiHub6roHnhyJTU1lXXr1rFq1SoWLVrEgQMHmDVrFh4eHhYf3uTm5rJ161ZiYmLYtm0b2dnZZRa48Ojq5dkLHxcfErIS2Bi70TxTqULT8UUAnszawK7TCTaMULDUAxOMs7Mz33//PW+99Rb+/v5ERUUxYMCAEjU+b948EhMTWb9+PYmJicydOzdvWXx8fJk9iyE8WhSSguf9ngdg8YnFmGSTeX67URglNd0Vx/htT2Sp2y/j045VSkk/uyJP8ur1etRqtcXzIf9J3getl5ubS0hICMHBwZw/f57U1FSCg4OZPn16iYIvC+Ikb8WkN+kJWRvCtYxrfPn4l/Ro0AOA3FVj0ZxcyXeGfvSesoDGNauVqN0LFy7g6OhIjRo1xCMkJSTLMsnJyaSnp9OoUSOLtinyHMz06dNRKpU0b94cZ2dnLly4QExMDM8999wDH2C0ZPuRI0fi6+vLjBkzGD9+PD/99BOTJ0+27F0KVYJaoWZUi1F8cuATfjjxA93rd0eSJDRBE+DkSoYp/2DevtO881S7ErXr4eFBXFyceESklHQ6XYku0hR7mfrmzZtERESQlpZGkyZNaNeuXZGZ/59Dxz5o+4CAALp27crJkyf54osvmDdvHt98843FgZcV0YOpuDL1mfRe05tbObdY0mcJ7Wqbk0nmN52xT4rm/5jAv//9EXYaMRRORVWi+2AsuemuuLGp7xcbG8vRo0c5cuQIUVFRbNy40dJQyoxIMBXbvGPz+C7qO7p6dOWbHnf+AR1bBuETOWFqyKn+GxjcvoFtgxQeqNj76FNSUhgxYgROTk5otVq8vb2ZPn06qampJdrRkSNHePbZZ3nyySd56623uHjxIo0bN2bQoEHMmjXLJslFqPiGNxuOTqljd9xuzqacNc9s8TQ5aif8FBfZ/9dO2wYoFKnYBDN9+nQaNGhAXFwcGRkZec8BBQUFcf36dYt3NHToUEJCQpg1axY+Pj4MGDCAbdu2lT5yoUpw1bkywHsAAEtOLjHPVOtQtHkWgICkcE5erfzPJz2qij1Eat26NceOHSsw/5dffmH37t189913+eY/6BCpQ4cOHDhwIO/75ORkunfvTlRUVClDLxviEKniu5J+hX6/9kMhKdg2aBs17WtC0ln4JoBMWctn/ut475kgW4cpFKLYHsyDzrkMHz68RH+YXl5ezJkzJ+86urOzs8XbClVbfcf6PF7/cQwmA6vOrDLPdGtCpvtj2Es5KKJXcjvHYNsghUIVm2ASEhJYvXo1p06dKtAzKcl9BDk5OcyfP58GDRrQp08f/Pz86NmzJ/Hx8SWPWqhynvU1HxKtPL0SvVEPgP1jYwF4hu2sPyp+jyqiYgtOTZs2jc2bN/Pf//6Xs2fP4u7uTosWLWjevDkJCZbfrr127VrA/ADk8ePHiY6OJjo6muHDh3P16tUCxbsF4X4BtQNo4tKEsyln2XppK/0a94Nm/cjW1qBZzhUW/bWF4R3HiZvnKpgSl2u4cOECJ06cyJuWLl2ab3lJLlNXBOIcTOWx5swaZkTMwN/Nn2Uh5oprhm3vo9r3JWuMnfAet5RW9Z1tG6SQj1Xrwdy1f/9+li1bxr59+7h27Rr29vb4+vrSt29fhg8fnq/GbnkTCabyyDJk0XNVT9Jy01j65FJa1mwJNy/A163JkdV87LuWGcOKv8NcKD/F1oNZs2YNKSkpeRXpSlIPBqBfv354eHjQv39/pk+fTs2aNcnOzubcuXPs3r2bZ555hpdffpnQ0NCHfzfCI81OZccgn0EsPrGYpafuJBjXRmQ26Ib95V3oYlZyKysIJ7vCn5MTyl+RPZj169fTr1+/fFeSiqtI988eTGpqarFXjCxZx1pED6ZyuXr7Kn3X9s1/yfrURljxLOdNddn7xGZGBVv2IJ5gfUVeRQoNDS1wmbok9WDAssvR4pK1YCn3au50r989/yVrnz5k62rhpbjGiX2bRDmGCqTIQ6Tly5cTHh6OUqlEluW8kpeWetiSm4JQmBG+I9hxeQerzqxibMuxqJVq1O1Hw55P6ZK2niOXh9PO07XYdgTrKzLB/PHHHw9V8vJhS24KQmECagfg7ezNudRz/H75d55o+ATKgFGY9nzGE4qDfLg3inaeYiysiqDIQ6T7S15u3bq1RCUvS1py82EHABeqDkmS8oqBrzi9wjzTyYPsRj3RSEYc/15BamauDSMU7ioywXz11VckJiayYcMGkpKS8uq16PX6YhsuruTmDz/8wJQpUxg0aBAtWrTg66+/tihgk8nEhAkTCAoKolu3bgVu0Pviiy/o2LEjHTt2ZObMmRa1KVQ+/Rv3x05lx8HrBzmfeh4A+yDznb1DpJ2sPlT4SKNC+SoywcycOZOoqChq166Ng4MD3333HWPHjiUiIsKixvV6PSEhIXTs2JGxY8fSv39/wNwzCg8PJyAggDVr1tClSxe+/PJLi9oMDw8nOzubiIgIPvnkE6ZNm5a3LDY2lqVLl7Jv3z4iIyPZtm0b0dHRFrUrVC7VNNXMd/NifnwAAO8eZNm746lI4EzEBnGytwIo8hzMnDlz8lWk69y5M1OnTrX4duyiSm6uX7+eVatWMXr06EIH936QvXv30qePedziwMDAfJeY69evz5YtW1AqzRXO9Ho9Op3O4raFymVo06GsOrOK9efXM6XtFOzV9mg6PA+7ZvH47d+IjP0XQV7iXJ8tWb2iXXElN7Ozs/nss8+Ij49n/vz5xcbw4osvMmjQIPr27QtAgwYNiI2NRaW6lytlWeaNN94gPT2dBQsWFGgjLCyMsLAwABITE7l06VKx+xUqpuc2PcexxGO8H/Q+z/g8A+nXMX7eHNkk857XCj76V29bh1ilWb2inaurKyEhIQwfPpyAgIACvR+dTse7777Lu+++a1F71atXz9fjMZlM+ZJLdnY2zz77LOnp6Xz77beFtjFu3DgOHTrEoUOHqFmzpkX7FSqmoc2GAuaTvbIsg2Mdcr36oJJM1DyzkhtpYiwuWyq3inb/LJl54cKFfMvr1atnUTvBwcF5MURGRuLv75+3TJZlnnrqKVq1asWCBQvyDpWER1dvz964aF34++bfRCWai5fZ3TnZO0y5k18ixdVJm5KL0apVq0LnL1u2TB4/fnyB+QqFotD1vb295aVLl8pHjhyRFy1aJLdu3VreunVrcbsvwGg0yuPHj5eDgoLkwMBA+dSpU/KcOXPkdevWyWvXrpW1Wq3ctWvXvGnfvn1FtteuXbsSxyBULJ8f+lz2W+InT/9zunmGySRnzmkty+9Xl9/44AM5R2+0bYBVWLHnYNq2bfvAERgLe45HlMwUyltcehxPrn0SlULFzsE7cdG5IEd+h7TlLSKMzUl6Zg39W7nbOswqqdwq2omSmYK1eDh60NmjM3qTnl/P/QqA1HoEeqU9QcoY/vhzt40jrLqKTTB3K9qNHj2amjVr4ufnx9ChQ5k5c2aJKtqJkpmCNQ1taj7Zu/L0SowmI+iqQ0vzvDY3VnPqWpotw6uyyr2i3T9LZkZHR9u0ZKY4RHo0GE1GQn4NIf52PPN6zKOLRxdIOAXfBpIha/mv3zpmDBYjD5S3cqloV5GJBPPo+OHED3xx+It8o0BmhvXB/moEH5pG8/Lbn4liVOWs2EOkrl27FtuIJesIgrUN9B6IRqHhz7g/ib9tPvS27zQRgOFsZdVB8XxSeSt2VIFDhw7RoUOHBy6XZZkbN25YvMPr169Tp04di9cXBEu56Fx4ouETbIjdwKrTq5jabio0DSHbrjZeWdeYu2cd+uCpqJXF/l8VykixCSYmJoadO3fSo0ePB65TkhvaevfuLR5AFKxmaLOhbIjdwNqza5nUehIapQZN4Fj440Oeyg5n0/GhPNXasps6hYdXbCr39PQkMjIST09PPD09SUtLy/takiQ8PT3x8PCweIdlfMpHEPJp6daSZq7NSMlJYevFrQAo2r+AQaHjcWUUW//4Q/wOliOL+or3/0Dmzp2b9/WHH35Y4h2KgbEEa5IkKe+S9fK/71RjtHeFNuaRIbsmryIiNtlW4VU5FiUYo9HI5cvmE2T3Jxvxn0CoiEIah+CkdSI6KZroRPPhuOqxychIDFDuZdUf4qphebEowUiSxLlz51i+fDk5OTmEh4eTlpYmeiNChWSnsmNQk0EALD115z6tGl7omzyJVjLgdXEZZ29YXoNIKD2LEsw777zDxYsXSU9Pp1evXtSqVYs//vijVDfHiSechfIwrOkwlJKSbRe3kZBpvuNc02UqACOVO/hxd4wNo6s6LEowjRs3ZsyYMYwdO5bnnnsOBwcHDAYDLVq0KPEOjx49WuJtBKGk6larS/cG3THIhnslNet3ILtOO5ylDDTHl5GQLmrFWFupbgho1aoVgwYNynfCVxAqmmd9zSd2V51ZRY4xBwDdnV7MaOk3Fu0+a6vQqgxxx5HwyGpbqy3NXJtxM/smWy5sMc9sFkJOdU8aKBJJOrBS9GKszKYJZv/+/bbcvfCIkyQprxez9NRS81VPhRJtl1cBmMgawnaJXow12TTBDB482Ja7F6qAvo364qpz5dTNUxxJuFM4rfWz5FbzoIkinpsHVpAg6vZaTbGPCjysIUOGFDpflmVu3rxp7d0LVZxWqeUZn2cIiw5jyckltKvdDlQaNI+/CRteYZK0hu92DeG90Ja2DvWRZPUEs2PHDn766SeqVauWb74sy/z555/W3r0gMKLZCH48+SO7ruziXMo5vF28ofUIcnd9hnf6ZW4dXMGNbj7Uri7G0CprVj9E6tatG46OjnTt2jXf1K1bN1q2FP81BOurYVeDgd4DAVh8crF5plJt7sUAk6Q1LPjjtK3Ce6RZPcGsXbuWLl26FLps+/bt1t69IAAwqsUolJKSTbGbuHr7qnlmq2HkVm+Al+Ia6YeWE5+aZdsgH0HiMrVQJXg4etCnUR8MsoH/xfzPPFOpRtP9bcDci/l00wkbRvhosnqC6dWrF4sWLSIxMTHffJPJREREBJMmTWLJkiXWDkMQGOM3BoA1Z9aQkp1inuk/BL1zYxopbuB4cikHLogLD2XJ6gkmPDycjIwM+vXrR6NGjWjbti2+vr54eXmxYMECRo0axejRoy1uz2QyMWHCBIKCgujWrVuhz0MlJibi4+NDdra4/Cjc4+PiQxePLmQbs1n29zLzTKUK9RMfADBNtYo56yIxmkSVgDJT1iO5PWhkR1mW5ZycHPnq1avyzZs3S93+mjVr5FGjRsmyLMsRERFyaGhovuVbtmyRW7duLTs6OspZWVnFtidGdqxaDl8/LPst8ZMfW/aYnJGbYZ5pMsmGH/rJ8vvV5R/eGSz/sv+SbYN8hJTrORiNRkPdunVxcXEpdRt79+6lT58+AAQGBhYYEUChULBjxw5cXV0f2EZYWBgBAQEEBAQUOHQTHm1ta7elTa02pOWm3SvlIEkon/wEGQXPKbezdst20rL1tg30EWHVBLN582ZefPFFnn76aYYNG8awYcN47733SE4ufUWxtLQ0nJyc8r5XKpUYDIa873v16kWNGjWKbGPcuHEcOnSIQ4cOUbNmzVLHIlROk1pPAmDxicXcyrllnlm7BQSMQSWZeEX/PXN3nLFhhI8Oq95op9frWbhwYb7CVHFxcVy6dKnYJPAg1atXJz39XrEgk8mESmX1+wWFR0hg3UAC6wYSeS2SH078wKvtzM8mSd3fwXB8FZ1yTvJzxBpOtvWghbtTMa0JRbFqDyYjI4Phw4czcuRIRowYwfLly/Hw8KBt27albjM4OJhNmzYBEBkZib+/f1mFK1QhU9pOAcwPQd7IuDPsjr0rqu7vAPC28if+vWI/OYbKM4hgRWTVBLNr1y6WL1/Ozz//zLJly9izZ89Dtzlw4EB0Oh2PPfYYr776Kl988QWff/4569evL4OIharCz82PXp69yDHmsCB6wb0FAWMw1fTFU5HAU8nf8+UO8bT1w7DqsUVubi5bt26lfv36xMXFFXnZODY2lsaNGxfbpkKh4Lvvvss3r1mzZgXWu3jxYonjFaqWl9q8xM7LO1l7di2jWozCs7onKNUoBs5HXtiDMaotjPizHYd9a9POs/QXJqoyq/Zg5s2bR2JiIuvXrycxMTFfBbwffviBKVOmMGjQIFq0aMHXX39tzVAEoYDGTo0Z4D0Ao2xk3tF59xa4t0Hq8gYAn6oX8N7KCDJzDQ9oRSiKJMtlO/aIUqnEaDQft+r1etTqgoON5+bm8swzzzB48GCee+45Jk6cyPz588syDIsFBAQUuNQtVB3XM64TsjaEXFMuS59cSsuadx7ANeoxLeyJ4voxlhu6EdN+Fh885WfbYCshq/Zgpk+fzptvvsmSJUsIDw/niy++YOzYsURGRrJ+/Xp0Oh2jR48mISHBmmEIwgPVcajDc82fA+CDiA8wmO70VJRqFE9/h0mpZZhqF1f2hxN+NN6GkVZOVu3BANy8eZOIiAjS0tLw8fGhXbt2+dbPzs7ms88+Iz4+3ia9GNGDEbIMWQxcN5D42/G8HvA6o1qMurfwr69h+/+RKFdnkPEj5o7vT6v6zjaLtbKxeoK56/Dhw8yZM4fU1FT8/f2ZMGECjRo1ylseHx9PvXrlPyi5SDACwJ64PUzaOQk7lR3hT4XjXs3dvMBkRP5pINKF3Rw3NWSS5iNWv9JDFKeyULk9KjBs2DD69evHrFmz8PHx4emnn2bbtm15y22RXAThrs4enXmi4RNkGbL4aP9H94ZFViiRBi9BdmmEv+Iib+Z8zbgfD5KtF/fHWKLcejAdOnTgwIEDed8nJyfTvXt3oqKiynL3JSZ6MMJdiZmJhIaHclt/my+6fUFPz573Fib8jWlRDxS5t/mvfgiXWkzkq2FtUCrE8MlFKbcejJeXF3PmzMn7z+Ds7FxeuxYEi9S0r8nUtlMB+Hj/x6Rmp95bWKsZikHfIyPxhnol+hPreW3lMQxGk01irSzKLcHk5OQwf/58GjRoQJ8+ffDz86Nnz57Ex4sz80LFMbjpYFrXbE1CVgJv730bk3xfAmnaB6nn+wB8rf6G1OjNvLoySiSZIpTbIdJdGRkZHD9+nOjo6Lzp6tWrhRaOKg/iEEn4p2u3rzF442Bu5dzi5TYvM67luHsLZRk2vQ4HF5ErqxinfxWHFk/y5bDWqJWiAu0/lXuCqWhEghEKc/eqkkJSENYrjI51O95bKMuw+U04EEaurGKCfir49OHLYa2prit4Y2lVZvWU27Vr1zJZRxDKU2ePzoz1H4tJNvHmn2+SkHnfzaCSBH0/hY4T0UgGFmi+QHNmIwO++YuzN9If3GgVZPUejIODAy1atHjg+rIsc+PGDS5fvlyWYVhM9GCEBzGajIzfPp791/fTplYbFvRagJ3K7t4Ksgxb34FI83NMXxsGsFAxlP8OaUMfv7o2irpisXqCuXTpEjt37qRHjx5FbuPh4VGWYVhMJBihKElZSQzdOJSEzASC6gYxt8dctErtvRVkGf76CnnnTCTZxB/GVkzRT+bpx/x4/YmmVNNW7WJoVj9E8vT0JDIyEk9PTzw9PUlLS8v7WpIkPD09bZZcBKE4bnZuLOy1EFedKxHXInj1j1fJNebeW0GSoNNUpJFrke1ceVwZxQbtuxyO+J1en+9me8wN2wVfAZTLae/7O0n3l2z48MMPy2P3gvBQGjs3ZlHvRThrndkTv4dpu6ehN/6jKLjX40jjdkGdlnhKCYRr3+PFjIVM+d9eJv58mCs3M20Su62VS4IxGo1551juTzZlfHQmCFbTxKUJC3svpLqmOruu7GLqrqmk5ablX8nFE17YBoGTUUjwgmozO7VvYIjZyOOf7eLN1VFcSs6wSfy2Ui4JRpIkzp07x/Lly8nJySE8PJy0tLR8xcAFoaJr5tqMsN5hOGoc+TPuT4ZuGEpMckz+ldR20OcjpLF/gHsb6krJLNR8zlLVTOKPbKb7nF28tvIYx66kVol/sOVyH0xsbCy7du3CaDSi0+nw8vIiMTGRr776it9//70sd19i4iSvUFJX0q4wbfc0Tt08hUah4a0ObzHYZ3DBf5gmIxxYCLs+gmzz8ChHTE2YaxjALlMrfGo7MaR9fQa2qYerg8YG78T6yjzBKBQKTKaib52Oiori3Llz7Nq1K985GVsQCUYojRxjDrMPzGbVmVUAdPXoypS2U2ji0qTgytlpcHAhRMyDTPOYYNdwY6WhM2uMXYiX6tC+oQs9mtWmh28tGtesVp5vxarKPMFUNiLBCA9jY+xGPoj4gCxDFhIS/b36M7n15Hv1ZO6XmwGHfoD9YXDr3n1fh00+7DC2ZZepFafkBtR3tae9pyvtGrrQztMFn1qOKCrpU9siwYgEIzykpKwkwqLDWHVmFQaTAbVCTY8GPejTqA+d6nXKf98MgMkEl/bCsWUQsw70964wJeDCHqMfR03eHDN5cVpugFqjpUltR5rVdqRpHUe8alXD09Ued2c7NKqK/fxTpUswJpOJSZMmERUVhVarZdGiRXh7e+ctX7hwIQsWLEClUvHuu+/Sr1+/ItsTCUYoK1fSr/DtsW/5LfY3ZMx/VtXU1Xi8/uO0q90O3xq+NHFuglp53/NKOelwbiec225+Tb+Wr80c1Jw2eXBedifWVJfzsjuX5Vpcl2uQIjlSx8mBes521KqupXZ1HbWra6nhoMXFQY2LvQYXew3V7dQ46lQ2eRiz0iWYtWvXsn79epYsWUJkZCQff/wx69atA+D69ev06tWLQ4cOkZ2dTadOnTh06BBarfaB7YkEI5S1+NvxbL24lS0XtnDq5ql8y1QKFV5OXtRxqIObnRu17GvhqnPFTmWHnUqHLj0BbeIZlElnUSadQXHrCgoZJMwTyNw9WNLLSlJwJFWuxm3sSMOeNNmeLLRkyxrzK2oMsgo9KlCqUKo0qFQqlCoVapUKtUqJUqFCpZRQq5SolAqUkgKVUkKpML8+/8RInKs5luqzqHT3Me/du5c+ffoAEBgYmC85HDhwgODgYLRaLVqtFm9vb6Kjo2nfvr2twhWqoHrV6jHGbwxj/MZw8dZFdsftJiY5hpjkGC6lXeJ0ymlOp5wuviEHwKGOBXs0ArfvTA/p7gXg+4aB6nezE87VSjdEc6VLMGlpaTg53RuQXKlUYjAYUKlUBZY5Ojpy69atAm2EhYURFhYGQGJiovWDFqqshk4NaejUMO/7DH0G51PPk5iZSGKWeUrJTiHbkE2WIYssYxa5xlyMJiNG2YhJNmGSTcjIBe+bkU1g1IPJALLR/GoymuffnUwmwGTeVr77KpsP4WTg7sHcP9q+f1/V7OxL/f4rXYKpXr066en3Hok3mUyoVKpCl6WnpxdamnPcuHGMG2cuIhQQEGDdgAXhPg5qh3uDu1UBFfsUdCGCg4PZtGkTAJGRkfj73+u6dejQgT179pCdnc2tW7c4deoUfn5iND5BsJVK14MZOHAg27dv57HHHkOWZRYvXsznn3+Ot7c3oaGhvPLKK3Tu3BmTycSsWbPQ6cT4NYJgK5XuKlJZE1eRBMF6Kt0hkiAIlYdIMIIgWE2VP0Ryc3OjYcOGtg6jWImJidSsWdPWYdiEeO8V/727ubmxZcuWAvOrfIKpLKryuSLx3ivvexeHSIIgWI1IMIIgWI1IMJXE3TuPqyLx3isvcQ5GEASrET0YQRCsRiQYQRCsRiQYQRCsRiSYSmznzp2MHTuWZ599lqioKFuHU2727dvHqFGjGDVqFKmpqbYOp1xVtp+5SDCVWGZmJmFhYbz++uts27bN1uGUm7CwMBYsWMALL7zAihUrbB1OuapsP/NKV66hKvvyyy/ZsWMHAEFBQbzzzjtkZGTw9ddfM3v2bBtHV37uDuBXt25dmw/cV9769+9fqX7mIsFUIlOnTmXq1Kl53yclJfHmm2/ywQcfUKtWLdsFVs7s7e3Jycnh2rVr1KljSc3aR0dl+5mLQ6QKYv/+/XTr1g0wlwGdMGECQUFBdOvWjXPnzhW6zWuvvca1a9d4++23Wb16dTlGaz2WfA7jxo1j/PjxLFiwgJEjR9ow2rJlyXuvdD9zWbC52bNny35+fnLHjh1lWZblNWvWyKNGjZJlWZYjIiLk0NBQG0ZXfqry5/CovnfRg6kAvLy8WLt2bd73RQ3N8iiryp/Do/reRYKpAAYNGoRafW+0vwcNzfKoq8qfw6P63kWCqYCKGpqlKqnKn8Oj8t5FgqmAihqapSqpyp/Do/LeK19KrAIKG5qlKqrKn8Oj8t5FuQZBEKxGHCIJgmA1IsEIgmA1IsEIgmA1IsEIgmA1IsEIgmA1IsEIgmA1IsEIgmA1IsEIgmA1IsEIgmA1IsEIFdr+/fsJCgrCzs4OFxcX/vOf/9g6JKEERIIRKqwdO3YQEhLCCy+8QFRUFG+++SbvvfceR44csXVogoXEs0hChZSbm0vTpk157733eP755/Pm161blzfeeIPXXnvNhtEJlhI9GKFC2r17N6mpqQVq7qrVarRarY2iEkpKJBihQvr999/x9/fPV+UtISGB+Ph42rZta8PIhJIQCUaokI4ePUpOTk6+ed9++y2NGzcmMDDQRlEJJSUKTgkV0t0E8/3339O5c2fWrVvH7Nmz2b59O5Ik2To8wUKiByNUOFevXiUhIYFly5Yxd+5c/P39WbZsGevWraNTp062Dk8oAXEVSahwNm3axLPPPktKSoqtQxEekujBCBXO0aNHadmypa3DEMqASDBChSMSzKNDHCIJgmA1ogcjCILViAQjCILViAQjCILViAQjCILViAQjCILViAQjCILViAQjCILV/D9Q1bGk0X1BbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def eD0Dt_norm_eval(r,ta):\n",
    "    \"\"\"Function just to evaluate the expression...\"\"\"\n",
    "    return(e_D0Dt_norm.evalf(subs={rho: r, t:ta}))\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(4,3))\n",
    "rhos = np.logspace(-3, 3, 50)\n",
    "tas = [0, 1e-2, 1e-1]\n",
    "for ta in tqdm(tas):\n",
    "    cur_d0dt_norm = np.array([eD0Dt_norm_eval(r,ta) for r in rhos])\n",
    "    ax.plot(rhos, cur_d0dt_norm, lw=2, label=r'$t_a = %0.2f$' % ta)\n",
    "    \n",
    "ax.legend(fontsize=10)\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel(r'$\\rho$', fontsize=14)\n",
    "ax.set_ylabel(r'$\\frac{\\mathbb{E}[D^{(0)}D^{(t)}]}{\\mathbb{E}[p^{(0)}_A(1 - p^{(t)}_A)p_B^{(0)}(1 - p_B^{(t)}]}$', fontsize=14)\n",
    "debox(ax);\n",
    "plt.savefig(main_figdir + 'eD0Dt_theory.pdf', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Corroborating theory w/ simulation results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:44.988857Z",
     "iopub.status.busy": "2022-01-02T14:38:44.988155Z",
     "iopub.status.idle": "2022-01-02T14:38:45.046374Z",
     "shell.execute_reply": "2022-01-02T14:38:45.047009Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>d0dt</th>\n",
       "      <th>std_d0dt</th>\n",
       "      <th>pApB</th>\n",
       "      <th>std_pApB</th>\n",
       "      <th>rec_dist</th>\n",
       "      <th>ta</th>\n",
       "      <th>scenario</th>\n",
       "      <th>maf</th>\n",
       "      <th>seed</th>\n",
       "      <th>ed0dt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.001312</td>\n",
       "      <td>0.005836</td>\n",
       "      <td>0.005246</td>\n",
       "      <td>0.011681</td>\n",
       "      <td>0.000050</td>\n",
       "      <td>0</td>\n",
       "      <td>SerialConstant</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11</td>\n",
       "      <td>0.250146</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.000558</td>\n",
       "      <td>0.002567</td>\n",
       "      <td>0.004923</td>\n",
       "      <td>0.010981</td>\n",
       "      <td>0.000149</td>\n",
       "      <td>0</td>\n",
       "      <td>SerialConstant</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11</td>\n",
       "      <td>0.113276</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.000378</td>\n",
       "      <td>0.001777</td>\n",
       "      <td>0.004866</td>\n",
       "      <td>0.010911</td>\n",
       "      <td>0.000249</td>\n",
       "      <td>0</td>\n",
       "      <td>SerialConstant</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11</td>\n",
       "      <td>0.077621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.000297</td>\n",
       "      <td>0.001414</td>\n",
       "      <td>0.004862</td>\n",
       "      <td>0.010842</td>\n",
       "      <td>0.000348</td>\n",
       "      <td>0</td>\n",
       "      <td>SerialConstant</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11</td>\n",
       "      <td>0.061098</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.000229</td>\n",
       "      <td>0.001069</td>\n",
       "      <td>0.004790</td>\n",
       "      <td>0.010660</td>\n",
       "      <td>0.000448</td>\n",
       "      <td>0</td>\n",
       "      <td>SerialConstant</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11</td>\n",
       "      <td>0.047914</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       d0dt  std_d0dt      pApB  std_pApB  rec_dist  ta        scenario  maf  \\\n",
       "0  0.001312  0.005836  0.005246  0.011681  0.000050   0  SerialConstant  0.0   \n",
       "1  0.000558  0.002567  0.004923  0.010981  0.000149   0  SerialConstant  0.0   \n",
       "2  0.000378  0.001777  0.004866  0.010911  0.000249   0  SerialConstant  0.0   \n",
       "3  0.000297  0.001414  0.004862  0.010842  0.000348   0  SerialConstant  0.0   \n",
       "4  0.000229  0.001069  0.004790  0.010660  0.000448   0  SerialConstant  0.0   \n",
       "\n",
       "   seed     ed0dt  \n",
       "0    11  0.250146  \n",
       "1    11  0.113276  \n",
       "2    11  0.077621  \n",
       "3    11  0.061098  \n",
       "4    11  0.047914  "
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ld_raw_test = pd.read_csv('../../results/ld_stats_raw/ld_stats_time_sep_raw.csv.gz')\n",
    "ld_raw_test['ed0dt'] = ld_raw_test['d0dt']/ld_raw_test['pApB']\n",
    "ld_raw_test.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T14:38:45.078424Z",
     "iopub.status.busy": "2022-01-02T14:38:45.077223Z",
     "iopub.status.idle": "2022-01-02T14:38:46.349474Z",
     "shell.execute_reply": "2022-01-02T14:38:46.350136Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 70.96it/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABIEklEQVR4nO3deVxVdf748dc5914u+w7uO7lCoqCCuKDmPlnmVFqapelYM01Z1qTNTNOqNWXfan6V6KRNpm1amuKupZSguGGK+4ob+84F7j2f3x/ARQSUnSt+nj1m8pzzOed+7knenM/5fD7vjyKEEEiSJNUDtbErIElS0yUDjCRJ9UYGGEmS6o0MMJIk1RsZYCRJqjcywEiSVG9sJsBomsasWbMIDQ0lPDycU6dOlSuTlJRE586dMZlMACxYsIDw8HDCw8MJDAykefPmAHzwwQf06NHDeuz48eMN+l0kSSombMSqVavE1KlThRBC7N69W4wbN67M8Y0bN4rAwEDh4uIi8vLyyp0/duxYsWnTJiGEEI8++qiIjY2t9zpLknRzNvMEExUVxahRowAICQkhNja2zHFVVdm6dSuenp7lzl29ejUeHh6MGDECgH379jF//nwGDBjA/Pnz67/ykiRVyGYCTGZmJm5ubtZtnU6H2Wy2bg8fPhwvL68Kz50/fz6vvvqqdXvixIl89tlnbN++naioKNatW1fheREREQQHBxMcHEyPHj3q6JtIklTCZgKMq6srWVlZ1m1N09Dr9bc87+jRo7i7u+Pn5weAEILnnnsOb29v7OzsGDt2LAcOHKjw3JkzZxIbG0tsbCwODg5180UkSbKymQATFhZGZGQkANHR0QQEBFTpvK1btzJ69GjrdmZmJv7+/mRnZyOEYPv27QQFBdVLnSVJurlbPyI0kPHjx7Nlyxb69++PEIKlS5eycOFC/Pz8GDduXKXnHT9+nOHDh1u33dzcePvttxkyZAhGo5Fhw4YxZsyYW35+gUWrk+8hSVIpRQg5mxqgfdcAzh073NjVkKQmxWaeYBpbRl5hY1dBqieFhYUkJCRYx09JNWdvb0/r1q0xGAxVKi8DTDFzgYnL6Xm0dJcve5uahIQEXFxcaN++PYqiNHZ1bltCCFJSUkhISKBDhw5VOsdmXvI2Njdy2PD71cauhlQPTCYTXl5eMrjUkqIoeHl5VetJUAaYYm5KDpGHrzR2NaR6IoNL3ajufZQBppgD+SReOMaVjLzGrookNRkywFxnjLqHjbKZJEl1RgaY64zWxchmkiTVIRlgSig6AtUzXDl/gmuZsjtTqnsJCQl88803dXa9jRs30qVLF/z8/FiwYEGdXbcuyQBTwt4VgFGymSTVk23btrF///46uZbFYuHPf/4zGzZs4OjRo6xcuZKjR4/WybXrkhwHU8LeHUhijC6GBYevMLV/+0aukFQf2r+8vl6ue27B2Jsej4qK4vnnn8fd3Z1NmzaxevVqOnbsWOPP27NnD35+ftZrTJw4kTVr1tC9e/caX7M+yABTwt4Voc+mt/kUCedOkpjZC19X+8auldREDBgwgD59+vDee+/h7+9/07IDBw4sk1mgxHvvvcc999wDwKVLl2jTpo31WOvWrYmJianbStcBGWBKKCrKXcMhfm1RM+nIIB4Lbd/YtZLq2K2eNOrT8ePH6dq16y3L7dq1qwFq0zBkgLlej/shfi2jdTG8f/iKDDBSnUlOTsbNzc2a4+hf//oXaWlpeHl58c9//rNM2ao8wbRq1YqLFy9ajyUkJNCqVat6/AY1IwPM9e4agdDb08d8gvNnT5GU1RsfF2Nj10pqAs6dO0fLli2BouaN2WzG3d2dX3/9tVzZqjzB9OnTh5MnT3L27FlatWrF119/zYoVK+q83rUlA8z1jC4ofvfAsXWMUPey8cggpoS0a+xaSU1A165dSU5Oxt/fH19fX9asWUNSUlKZp5Dq0Ov1/Oc//2HkyJFYLBamTZtmk2lfZYC5Uff74dg6xuhi+PDwFRlgpDrh7OzMnj17AHj//fd57733SElJoVevXjW+5pgxY6qUTK0xyQBzo84jETo7+ojjnDpzmuTsXng7y2aSVHdeeOGFxq5Cg5ED7W5k74rSaRiqIhih7mXTETnoTpJqSgaYivS4H4AxagwbDssAI0k1JQNMRTqPQqgG+qnxnDhzhtScgsaukSTdlmSAqYiDO0qnIegUwT2KbCZJUk3JAFOZ7vcDMFqVKRwkqaZsJsBomsasWbMIDQ0lPDycU6dOlSuTlJRE586drTlBhRC0atWK8PBwwsPDmTt3LgA//fQTffr0ITQ0lMWLF9esQl1GI1Q9oepR4k+fJU02kySp2mymm/rHH3/EZDKxe/duoqOjeeGFF1izZo31+KZNm3j55Ze5erW0uXL69Gl69+7NTz/9ZN1XWFjI7Nmz2bt3L05OToSFhTFu3DiaNWtWvQo5eqJ0GIz+9DaGKbFsPtqfh/u0rfX3lKQ7ic08wURFRTFq1CgAQkJCiI2NLXNcVVW2bt2Kp6endd++ffu4dOkSQ4YMYcyYMRw/fpz4+Hj8/Pzw8PDAzs6OAQMGsHPnzlt+vqCC9eeu602KlL1JklRtNhNgMjMzcXNzs27rdDrMZrN1e/jw4Xh5eZU5p0WLFsydO5cdO3Ywb948Jk+eXO46Li4uZGRkVPiZERERBAcHExwczOXky+ULdBmLUHT0V4/w+6lzpOfKZpJUc3Wd0W7atGn4+vpWmP6hsmx3DZ0Fz2YCjKura5kZpJqmWWeeViY4OJj77rsPKMq3cfnyZVxcXMpcJysrC3d39wrPnzlzJrGxscTGxqLYV7Acg5MXSoeBGBQLQ5RYNh+9Vv0vJknF6jKjHcDjjz/Oxo0by+2vLNtdY2TBs5kAExYWRmRkJADR0dEEBATc8pzXXnuN//u//wPg0KFDtGnThu7du3Py5ElSU1MpKChg586dhIaG3vJamQWZFGoVLB9b3Js0RvYmSbVQktHu+++/JzAwkDNnztT6moMGDSrzyqDE9dnu7OzsrNnuKttfn2zmJe/48ePZsmUL/fv3RwjB0qVLWbhwIX5+fowbN67Cc15++WUmT57M+vXr0ev1LFu2DIPBwMKFCxk5ciSapjFt2rQq5cmwCAuxV2MJbXlDMOr6B8T65xmgHmbOqfNk5PbCzbFq6/JKNuhfbrcuU6PrVtwML1HXGe1uprJsd42RBc9mAoyqqnz22Wdl9lWU/evcuXPWP3t4eLB+ffkcq/feey/33ntvteuw5fyW8gHG2QelXRh253YxWOxjS3wofwxqXe1rS5LMaHeH23ZhG6/0ewWdqit7oMf9cG4XY3QxrDw8QQaY29ktnjTqS11ntLuZyrLdNUYWPBlgitnp7Eg1pbI/cT99mvcpe7DrvYj1cxikxvHSyfNkmgJxtZfNJKnq6jqj3c1Ulu2uS5cuDZ4Fz2Ze8jY2V7uidZE2n9tc/qBLM5R2/TEqZgaK/WyVvUlSNV2f0W7KlCn87W9/Y+rUqbV6gpg0aRKhoaEcP36c1q1b89///hcom+2uW7duPPTQQ/To0aPS/fVJEUJUMMLsztMjsAfqbBVvB2+2PbgNVbkh9sZEwIYX2WQJ5ju/BSyZ2qfiC0k2Jz4+nm7dujV2Nazef/99MjMzSUlJoUuXLjzzzDONXaVqqc79lE2kYg56B7ydvbmUfYmDiQfp3ax32QLd7oUNLxKuHuLlExfJMgXiIptJUg3IjHZ3qBHtRgCw+XwFzSTXFtAmBKNSSJjYz7b4xAaunSTdfmSAuc7wdsOBou5qTWjlC5TMTdLFsF4OupOkW5IB5jr+3v40d2pOYm4icUlx5Qt0KxpbM0Q9yJ4TF8nON5cvI0mSlQww11EUpcxTTDluraF1HxyUAvprB9gWL3uTJOlmZIC5Qcl7mC3nt1BhB1vx3KSxOjk3SZJuRQaYG9ztcze+Dr5cybnCkZQj5Qt0L5oXNUQ9QPTxBHJkM0mSKiUDzA1UReWedkXDsSvsTXJvCy1746TkE6IdZPsx2ZskSZWRAaYC1vcw5yppJl3XmySbSZJUORlgKtDLtxde9l4kZCdwLPVY+QLdippJw9T9/Ho8gdwC2UySpIrIAFMBnaqzNpMq7E3y7AAteuKsmOhrOcSOY0kNXEPpdmTLKTPrK5WmDDCVKGkmbT6/+aa9SaNlM0mqIltNmVmfqTRlgKlEULMgPIwenM88z8n0k+ULdC/KBTxc3ceuY5fIK7CUOaxpgqSs/IqDk3THseWUmfWZSlNOdqyEXtUzrN0wvj/xPZvPbaazR+eyBbw6QbMAXK8dJshyiJ+PBzM6oAUAZrPGgxG7iUvIILidBytnhKCqFSQVlxpcwBe3zvVcE4enHr7pcVtOmVmfqTRlgLmJ4e2G8/2J79lyfgt/6fWX8gV63AfXDjNGjWH94TGMDmiBpgkeitjNgQvpAMSeTyMlpwAfF2PDVl6yOTJlplRGn+Z9cDO6cSbjDKfTT9PJvVPZAt3vh+1vMkIXy5vHLmMq7EmWycyhi+nWIne3dsPb2a5B6y1V7lZPGvXFllNm1mcqTRlgbsKgGhjaZig/nPqBzec385T7U2ULeN8Fvt1xSzxKYGEcPx8PZmR3H4a1Udh2UaNnGw++nxWKosjm0Z3OllNm1mcqTRlgbmF4u+H8cOoHtpzfwlM9nypfoPv9kHiU0WoMkXGjGLV3GouSYyjs1AfDtPUoqnyPLpVNmenr68uaNWtISkoq8+RQXZMmTeLnn38mOTmZ1q1b89prrzF9+vQyqTEtFgvTpk2zpsas7v7aspmUmZqm8fTTT3Po0CGMRiNLlizBz8+vTJmkpCTCwsKIi4vD3t6ejIwM63KxBQUFLFy4kNDQUH744QfmzJljfXH12muvMXjw4Jt+fnBwcLn1sAEKLYUM/mYwWYVZrL1/LR3cOpQtkHgMPulHmnAmXHxGjG4G9iIPoeph9lEUl2a1uzFSrcmUmXXrtkyZ+eOPP2Iymdi9ezfR0dG88MILZbrKNm3axMsvv8zVq6WL0C9cuJBhw4bx3HPPcfz4cSZNmsT+/fvZt28f7777LhMmTKh1vQw6A0PaDmHt6bVsPb+VGXfPKFvAtyt4d8Ej+TgBhUf4xXsCw1JXss/sx8KvTrNypq/sQZLKkCkzG0FUVBSjRo0CICQkpNzThKqqbN26tUy//+zZs/nTn/4EgNlsxt7eHoB9+/bx+eefM3DgQF544QXM5toN5b9pjhgonZukxvCDx+PcW/A2Dxf8nX0X0knJKajVZ0vS7cxmAkxmZiZubqXLeup0ujKBYfjw4Xh5eZU5x93dHQcHB65evcrkyZOZP3++tezHH3/Mzp07yc7OLrdiZImIiAiCg4MJDg4mKany4f6hLUNxMjgRnxrPxcwK2szFg+5G6vayO/483xte5Ru7NwluK3uQpDubzQQYV1fXMl1zmqZZu/Ru5vDhwwwbNoy3337b+p5l2rRpdOzYEUVRuO+++zhw4ECF586cOZPY2FhiY2Px8fGp9DOMOiODWxdde8uFCp5ifLuDlx9eShbdOcNurTt9dCf5alJH2YMk3dFsJsCEhYURGRkJQHR0NAEBtx5xefToUR588EFWrFjB6NGjARBCcPfdd5OQkAAUzf8ICgqqdf2sme7OVRBgFMU6N2mMGsM6Sz8Om9sQ+MEBzOYKkodLDc5G+jJue9W9jzYTYMaPH4+9vT39+/dn9uzZfPDBByxcuJC1a9dWes7cuXMxmUw8++yzhIeHc99996EoCkuWLOGBBx5g8ODB5ObmMmPGjEqvUVVhrcJw0Dvwe8rvXM6+XL5AcTNplG4v27VedFES+FVM49S19Fp/tlQ79vb2pKSkyCBTS0IIUlJSrO86q8JmuqkbW2Xd1Neb88scNp3bxJzgOUztMbXsQSHgo16QdpaH8//Bn/Q/MUR3EPGnX1Fb3HzuiVS/CgsLSUhIwGQyNXZVbnv29va0bt0ag6Fqiw7aTDf17WBEuxFsOreJzec3lw8wilLUmxT1AaN1May39KOf/gQOPrYz/uJOZTAY6NChw60LSnXOZppIt4MBrQZgr7MnLimOqzlXyxcobiaN1u3hRy2MgLwIUvJktjvpziUDTDU4GhwZ2HogAFvPby1foEUgwq0tzZR0eimnMBr0eDrK9aulO5cMMNV000F3QqBZCgG4V7cbU6GF1NzChqyeJNkUGWCqaVDrQdipdhxIPEBi7g1LluQkomYXpc8cr4vCgXxSc/IboZaSZBtkgKkmJ4MTYa3CEAi2Xdh2w1Gl+H/gquTxkO5nFu8629BVlCSbIQNMDVTaTHL0Rtg5U9LxP023gbUHLpKYKbtHpTuTDDA1EN4mHL2qZ9+1fSTnJZceyEtBKcxFUYqGxbRVkxjGHpb9dq7R6ipJjUkGmBpwsXOhf8v+aEJj+4XtpQecfBBt+mEWOs6JojwwM/TrWR59nmy5hrV0B5IBpoYqbCYpCidGrSQk/z+MKXibdOFEb/UUfvlH+XZvzTOXSdLtSgaYGhrSZgh6Rc/eq3tJM6VZ93du5oyzvY487PlWKZogOUO/nv9GncVskRMfpTuLDDA15GZ0o1+LfliEhR0XdxTt1DTUL+5lh/IUR9p/yJPPvYnQ2TFSF4su4xwbfq9g9K8kNWEywNTCiPZFTyibz28u2pGTCBd+QxEWnK7uQdXpUQIeREUwTbeBiJ1n5Ixe6Y4iA0wtDGkzBJ2iI+ZyDBn5GaDdEDw0AaF/BuBh/S9cuHSJ6DOpjVBTSWocMsDUgoe9B32a98EszOy4uANNgWRVxRpmVBWa9YBOQ3Egn0d121i8q/ZrEkvS7UIGmFoq6U3afG4z036byz1tW/NE82ZobUPB2beoUGjRsrOP6zez69hlTl4rv2qfJDVFMsDU0tC2Q1EVld1XdnPg2gEsChx0dCR14vKiHDEAnYaCb3d8lTTGqb+xRE4fkO4QMsDUkreDN0HNgjBrZtq6tkWn6Aj0DcTL4boVEBTF+i5mhn49PxxIkNMHpDuCDDB1oKSZ1MGtA1sf3MrSkUvLryYQ8CA4N6OrepG+Io4vdp9r+IpKUgOTAaYODGs7DAWFXy/9ioPeoeKlSvRG6FuUfHyGbj3Loy+QI6cPSE2cDDB1wNfRl16+vSjQClh/Zn3lY12Cp4PegcG6OJqZzvBdrJw+IDVtMsDUkXva3QPAG9Fv8MSmJ9BEBdMCHD2h16MAPKmLZImcPiA1cTYTYDRNY9asWYSGhhIeHs6pU6fKlUlKSqJz587W5Sfy8vKYMGECAwcOZMyYMdblX3/66Sf69OlDaGgoixcvbpD6h7QIsf75wLUDpJoqGVAX8jQChfv1v5KfdoWNR+T0AanpspkA8+OPP2Iymdi9ezcLFizghRdeKHN806ZNjBgxgqtXS38gP/30UwICAti1axePPfYYb775JoWFhcyePZvNmzfzyy+/EBERwbVr1+q9/n7ufngYPQDwcfTBy96r4oJenVC6jsUOM4/pN7NYTh+QmjCbCTBRUVGMGjUKgJCQkHKLoKmqytatW/H09KzwnNGjR7N161bi4+Px8/PDw8MDOzs7BgwYwM6dOyv8zIiICIKDgwkODrY+/dSUoih8NPQjAHLNueSZ8yovXDzwbop+KycSrrHnrJw+IDVNNhNgMjMzcXNzs27rdDrM5tJeluHDh+Pl5VXpOS4uLmRkZJS7Tsn+isycOZPY2FhiY2Px8fGp9XcI9A0k0CeQrIIsfjj1Q+UF24ZAqyDcyWaCbhcRO+X0AalpspkA4+rqSlZW6RB6TdPQ62++8OT152RlZeHu7l7uOiX7G0rJio/Ljy7HolkqLqQo1qeY6foNbD92lVOJcvqA1PTYTIAJCwsjMjISgOjoaAICAqp1zoYNGxg4cCDdunXj5MmTpKamUlBQwM6dOwkNDa3XupfQhEaAdwCtnVuTkJ1QmiemIt3GgVtbOihXuUfdL6cPSE2SzQSY8ePHY29vT//+/Zk9ezYffPABCxcuZO3atZWe89RTT3HkyBEGDBhAREQEr776KgaDgYULFzJy5EhCQ0OZNm0arVq1qvf6mzUzUyKnMHLVSOtAuy+OfFH5CTo9hDwFwAx9JKv3XyIxS04fkJoWRcguDACCg4PLvViuKk1oTF4/mcMph637nA3OZBdms3zMcnr69Kz4xPwsWNgD8jMYl/8Gg8JHMmdklxrVQZJskc08wdzOkvOSywQXgD90/AMA/zvyv8pPNLpAUNE7mxn69XwZfZ7cAjl9QGo6ZICpA+oNt9Hfy59p/tPQK3q2XtjKpexLlZ/cbxaoesbo9uCcd5nvYhPqubaS1HBkgKkDXg5eBDULQqfo6O7ZneWjl9PCuQWjOoxCExpfxX9V+cluraDHA+jQeEK/kSVRZ7DcmHpTkm5TMsDUAUVRWDJiCT28enAs9RjTt0xHExqPdX8MgNUnV5NVcJNu6P5FXdaT9D+TnprCRrn6gNREyABTR9Lz0zmScgQNjYOJB0k1pdLNqxt9m/clpzCH1SdXV35yi57QfiBO5DFRt52InafLTR/QNEFSVr6cViDdVmSAqSNe9l4E+gaWZrQrnotU8hSzPH45hVph5Rfo/wwA0w2bOJKQwt5zpYu5aZpg0uJoQudvY2JENJpsQkm3iZsPlZWqTFEUPh/5OammVLzsvaxjYQa2Hkh71/acyzzH1vNbGd1hdMUX8BsO3p1pnnyCMWoMETtb0bdD0byrlJwCYs+nYdEEsefTSMkpwMfF2FBfTZJqTD7B1CFVUfF28C6T0U5VVKZ0nwIUDbyrtImjqhDyNAAzDZFsjb/KqcRsADwdDTja6QBwtNPh6Wiox28hSXVHBpgGcG+ne3E3unMk5Qj7E/dXXrDnRHD0xl85S4gaz3+jiiZBpuYWklucXjM330xq7k2aWpJkQ6oUYPLy8rh0qfxYjiNHjtR5hW5nmtBIzksu95TioHfg4S4PA7cYeGdwgD5PAvCkbj2r9l8iKSsfb2c7gtt7olcVgtt74u1sV2/fQZLq0i0DzPfff89dd93F2LFjufvuu4mJibEemzJlSr1W7naiCY3HNzzO0G+HMnXD1HIpMyd2nYhBNbDj4g7OZ56v/EJ9ngSdkXt0B2hjuciXu8+hKAorZ4Swe+4wvp4ZUnFScUmyQbcMMG+++Sb79u3j4MGDLF26lOnTp7NixQqAcr+p8/JukmSpiUvMTeRA0gEEggNJB0jMTSxz3NvBmz90/AMCwfKjyyu/kLNPUVMJmK7bwP+iz5NXYEFVFXxcjDK4SLeVWwaYwsJCmjVrBkBQUBA7d+5k0aJFvP766+X+sg8YMKDc+ceOHaujqtq2zPzMm24D1pe9a06vISO/4iRYgDVXzAR9FGpuCt/tk6sPSLenWwYYX19f4uLirNuenp5s2bKF+Ph46/6ffvqJd955h+zsbC5eLPvD8PDDD9dxlW2Tn7sfzgZnoGgmtZ+7X7kyd3ncRVjLMPLMeXx34rvKL+bTGe4aiZECpui2sGTXWTl9QLot3TLAfPnll/j6+pbZl5CQwMqVK/nll18A8Pf3x8XFheTkZKZOnUqnTp0YNGgQDz/8MAbDndGlqqoquybuYtW4Vfw68VdUteJbWzLwbkX8CgotNxt4V/QU87hhK9dS09lw+ArHrmSiaXKZE+n2UaN8MDqdjq+++oqJEyeW2b9z504GDRoEwKVLlzh//jz+/v64urrWTW3rUW3ywVSHEIIH1j7AqfRTvDXgLcZ1GldZQVg0CK7G8XLhk3yrDUUT4GKv58Dfh6PXyxEGku2r0d9SIQQffvghXbp0oWvXrkyZMoUtW7ZYgwtAq1at6N+//20RXBqSoijWp5j/Hflf5QPvFMU6fWCmYQNCaIAgy2TmVFJ2A9VWkmqnxr8GL1y4wIQJE5gyZQrZ2dncd999PPnkk/IRvgrGdhyLl70Xx9OOE3M1pvKCPcaDa0s6colw9RB9lGPoFQt+Pk4NV1lJqoUaB5gVK1bw9ttv88orr/DDDz8QFxdHVFQU77zzTl3Wr0my09kxqesk4BYD73QG6PkIADN164gVXXhf9wlpuQUNUU1JqrUaBRhvb+9yL379/Pz48MMPWbJkSZ1UrKl7qMtD2Ovs2XVpF2fSb7Iukt9wBBCqi6eHco7vtcF4cZMubkmyITUKMIGBgURERJTb365duwqnFEjledh7WF/w/u/oTZ5i2vQFpagn7hX9V+zSAth8XnZZS7eHGgWYN998k4iICB566CF+/vlnUlNTuXTpEm+88QYdO3as6zo2WZO7Twbgp9M/kZKXUmm5Qp8eCFH0FDNc3cdr645iKqxkUTdJsiE1CjB9+/YlJiaGlJQUhg8fjo+PD23btmXNmjUsXLiwRhXRNI1Zs2YRGhpKeHg4p06dKnN88eLFBAcHExISwrp16wB47rnnCA8PJzw8nK5duxISEgLAs88+S1BQkPVYZUvHNrYObh0Ibx1OgVbAt8e/rbhQ5mUMiQcpGTT9D/2XpGZk8tkvpxuuopJUU6KWkpOTxaZNm0RkZKRISUmp8XVWrVolpk6dKoQQYvfu3WLcuHHWY1euXBH+/v7CZDKJ9PR0659LFBQUiL59+4q4uDghhBBhYWEiKSmpWp8fFBRU47rXxp4re4T/Mn8x6OtBwmQ2lS9weqfQXnUV4lVXof2z6N/vzntSdH4lUlxIyWn4CktSNdR6tJaXlxcjRoxg9OjReHp61vg6UVFRjBo1CoCQkJAyg9727NlDWFgYRqMRNzc3/Pz8ykxf+PjjjxkxYgQBAQFomsbJkyeZOXMmYWFhfP7555V+ZkREBMHBwQQHB5OUlFTjutdGcLNgunl2I9WUyrrT68oXaBOCAgiA4qeYv9qtxdOcxFvr4xuwppJUfTYzHDQzMxM3Nzfrtk6nw2w2V3jMxcXF2uwpKChg0aJFzJkzB4CcnByeeeYZli9fzsaNG/nkk0/KBKPrzZw5k9jYWGJjY/Hx8amvr3ZTiqLwWI/igXdHKxh4l5+GQCmOLQqi8yiMwsTfjSvZeOQqu042TmCUpKqwmQDj6upKVlbp0h6apqHX6ys8lpWVhbu7OwBbt25l0KBB1gDk6OjIs88+i6OjIy4uLgwdOpRDhw413BepgZHtR+Lr6MuZjDNEXYoqc0xz8OaooQeFQiXG0oU/p05E6O0Zq/xGXyWef609QoFZDm6UbJPNBJiwsDAiIyMBiI6OJiAgwHqsb9++7Nq1C5PJREZGBvHx8fj7+wNFAWb06NJE2idOnCAsLAyLxUJhYSFRUVH07t27Yb9MNRlUA492exQo32WdklvIuKyXCc3/f0ws/AebLhvJ7VM0EfJth+WcTcrii9/ONXSVJalKbCbAjB8/Hnt7e/r378/s2bP54IMPWLhwIWvXrqV58+b89a9/ZeDAgQwdOpS33noLe3t7AI4fP16ma7xbt25MmTKFkJAQBg8ezGOPPUaPHj0a62tV2YS7JuCgdyD6SjTHU49b93s6GnAwGkjGDVAIauuO45Dnwa0NftpZJum28+G2kyRmmhqv8pJUiRrNpm6KGmo29c0s2LOAr+K/Ylyncbw14C0AkrLyCZ2/BXctk1TFlei59+Drag9HfoTvppKtuhCW+z7DenVh4cOBjVp/SbqRzTzBSPBot0dRFZXIs5HWlJveTnrWOi8g2vhn1rkswMe5OL9O9/ug/UCctSxetPue1Qcuse98aiPWXpLKkwHGhrRxacOwtsMwa2a+PvY1AEpuCt3M8egVjW7meJTc4hG/igKj3wFFxyPqVroqF/jnmiNlMt9pmuBaponETJNcclZqFDLA2JiSXDHfnviW3MJccPJBadMPVH3Rv52u605v1gP6TEdF4y375Ry5nMHKPReAouAyMWI3/d7eRt+3t/HQZ7vlkrNSg5MBxsYE+gZyt8/dZORnsPb02qInlanr4Pl4eHw93LiqQPhccPAkSPzOGDWG9zYfJy2noGi52evWt957Po2k7PwG/jbSnU4GGBtU8hTz5dEvi9ZXUlVw9i0fXAAcPWHYPwB43WElptxs3t9yHG9nO3q2cStTVC54IjU0GWBs0LC2w2jl3IoLWRf4+eLPtz6h91RoHoC3JYmnDT/xVcwFjlzO5PtZ/enVxh2dqtC3gyc+Lsb6rroklSEDjA3Sq/pKB95VSNXB6H8D8JR+Ha1I4l9rj6CqCque6k/03GF8I1eElBqBDDA2arzfeJwNzuy7to8jyVVYA7xdKAQ8iEEU8Jr9SmLPp/HjQZn8S2pcMsDYKGc7Z/7Y+Y8AfHzg46p1M9/zGhgcGSai6a/+ztvrj/HkJ5GEzt/KxIho2YskNTgZYGzYI10fQUHh18u/MmHthKIXvjfj1goGvgDAAoevSMvOocOV9SzXvc6+8ymk5Mhk4VLDkgHGRmlCY87OOYiiTDCcTD9pHd17U6F/AY/2tLWcZ4puC19YRuJOFq0M2Xg63hmrbEq2QwYYG5VqSi337mXTuU23PtFgDyPnA/CS4VtcyWWO+SnO5zuTLJ9gpAYmA4yN8rL3ortX9zL7Pj30Kcl5ybc+ucto6DQMB/J5Sf81v4uOYE1aJUkNRwYYGyUQ2Kl21m1XO1dyCnP4cP+HVThZIHJSEMDDup/pqZxCAa7JlA5SA5MBxkalmlI5mHQQAJ2i4/8N/X8YVAM/nvqRuKSKU4Ba5SbD1UMoFA3+/a/dv3Ekj79+fZDsfHO9112SSsgAY6O87L3o6dMTFZXOHp0J8A6wTiGYHzP/5j1KTj7QIpCSTmlvJYsPnJdzNjmHV9dUYUyNJNURGWBslCj+R0MjPjWesK/DmOY/DV9HX35P+Z0fT/1Y+cmKQuqw9ymJMELACPMOJtjtZtX+BH48IAfgSQ1DBhgblZKXwoHEA9btHHMOV3Ku8EJQ0TiXD/d/SGZBZqXne7bzB6UouJS83Z1vt5TWSiKv/HCYc8k59Vl9SQJkgLmteBg9GN1hNEHNgkg1pfLJwU8qLSuSTgClE7BF+4HYmbP5n9tiTAUF/PXrA2VWI9A0QVJWvkxMJdUpGWBslLeDN718elm3e/n0wsfRB0VRmNt3Lqqi8vWxrzmZdrLC81NxtTaREJA66C1waUlH0xH+7rSWuIQM3ttclFxc0wSTFkcTOn+bnFIg1SkZYGyUoigsG72MbQ9uY8eDO/hi9BfW2dBdPLvwYOcHsQgL8/fMr/Cpw6t5G+INPTALhXhDDzzbB8ADEYDC49oqQnTHiNh5hp+PJxYlpzqfhlkTxJ5Pk1MKpDpjMwFG0zRmzZpFaGgo4eHhnDp1qszxxYsXExwcTEhICOvWFS2xmpqaire3t3WR+w8//LDSsrcjVVHxdvCuMFPUM72ewd3ozt6re9l0vvwIX0VV6Tp3F2mz4vB+ZltRW6nDQBj4PIrQWOy8CFeymfPdISwWDUc7HQCOdjo5pUCqO424LnYZq1atElOnThVCCLF7924xbtw467ErV65YF7xPT0+3/nnLli3iL3/5S5nrVFb2VoKCgur0+9QFi2YRUzdMFT2/6CmmbpgqLJqlzPFvjn0j/Jf5i2HfDhM5BTnlz7do4qHPfhOd5q4XD332m7BYNCHMBUJEDBXiVVcRPX+MaPe3n8SDn/4qgv/2lWj3t59Eh7+tE4mZt75fklQVNvMEExUVxahRowAICQkps0bRnj17CAsLw2g04ubmhp+fH3Fxcezbt499+/YxePBgHnzwQa5cuVJp2dtRqimVg4kHsQgLBxMPkmoquyzJhLsm0M2zG9dyr7Hk8JJy56fkFLCvuOmzr6TpozPAhCVg50I/UxTTHXay51waT+g28LXhDZyMinyCkeqMzQSYGxe41+l0mM3mCo+5uLiQkZFB165def311/nll1+4//77eeaZZyotW5GIiAiCg4MJDg4mKcn2FpH3svci0DcQnaIj0DcQL3uvMsd1qo55/eYBsOzIMi5mXixz3NvZjqB2HuhVhaB2Hng7F0898OwAf1gIwDx1GZ2USyy0/BEdFoz56XJSpFRnbCbA3LjAvaZp6PX6Co9lZWXh7u7O0KFDGTJkCFC09OyBAwcqLVuRmTNnEhsbS2xsLD4+PhWWaUyKovD5yM/Z+uBWlo5cWmHKy0DfQO7teC+FWiHv7n233PkrZ4Swe+4wvr4xZebdD8HdE9FZ8vnC7l1UBH8ufI5kXGVXtVRnbCbAhIWFERkZCUB0dDQBAQHWY3379mXXrl2YTCYyMjKIj4/H39+fJ598klWrVgGwbds2goKCKi17uyp50XuzfLqzg2bjZHDi54Sf2ZWwq+z5qoKPi7Hi88f8G4tbO1orSbyrX0QiHgAUmm+R2EqSqshm1qbWNI2nn36auLg4hBAsXbqUyMhI/Pz8GDduHIsXLyYiIgJN05g3bx4TJkzg7NmzTJs2DSEETk5OLFmyhBYtWlRY9lZsYW3q2lj2+zLe3/c+7VzbsXrcaux0drc+CRAnNsOKB1GA5wtmsVobxIjuzfhschCqKhM8SLVjMwGmsd3uAabQUsgDax/gXOY5nuv9HNMDplftREsh4g0fFAQWoTC+4DXihB9/GeLHnJFdrMU0TZCSU4C3s51cnUCqMptpIkm1Y9AZmNt3LgCL4hZxLeda1U5MPolSPORXpwi+tXuD5koa/9lxim/2li5DK0f6SjUhA0wT0r9Vf4a2GUqeOa/SEb7leLRHUDSrQAiwVwpZ7/P/cMDEvB9+Z+eJJDnSV6oxGWCamBeCX0BBYduFbVVaiUA7vhkoHSycJFzxyjzKmhbLEJqFp7/aT1KmCQdD0V8VB4Mqx8lIVSYDTBNj1BnLrERwOv30TcunevQsnRQJLHSZi7B3o3PaTha3WEN2vpnHl+0lO98CQHa+RY6TkapMBpgmJj0/vcz22zFv37Sp5NWyHUfVrhQKiNH8OOXUE/HQclANDEv7jnk+v5KYlY8eM56kA0ImD5eqTAaYJuYu97tw0jtZt2OvxfLl0S8rLa+oKt7P7iCs8FMmFr7GgYsZpPj0g3EfATAj+zMeMkRhRs/d6lm+Mb6Fl6O+3r+H1DTIANPEqKrKqnGryixS8sH+D/g9+fdKz/F1c6BDuw7oVbV0SkHgIzDoRRRhYYH6KX2UY/ys9eKLwntISrzcEF9FagLkOJhit/s4mBJmzczArweSXZgNgI+DD0l5SbRybsW3936Lq51rhedVOM5FCPh+GhxZTYHQMzJ/AWdpyegezfjokd4YdPL3k3Rz8m9IE3Mm44w1uAB8NOQjunt151L2Jf71278qfR9T4ZQCRYFh/0IAdoqZSOM8vEhnw5FrPPv1AQotRT1UMt2mVBkZYJqYjq4dcdQ7AuBscKa7V3feG/QeTgYntpzfwjfHv6neBQ321j86KAWstHubNva5RB6+yl9XHiC/0MLEiN2EvL2VhxftloPwpDJkgGlCNKHx5JYnyTXnAkUvfFGgjWsb/tX/XwC8u/dd4lPiq35RJ2+EzhEhwCKgs5pApPt7tLbPY8PvV/nT8n3sOZeGRcCec2kkZefXwzeTblcywDQhJQmqShxMKk1SNar9KB7q/BCFWiFzfplDTmEVly3JvIJiyUVRiv6ynNN8cUk/xhav92ljn8fPx8vm0dEscia2VEoGmCbEy96LQJ9A67aDzgEPo4d1+8U+L9LZozMXsi7w2u7XqvbOJLnsQL35BQ+T79YRh5SjbPb6gJbGvLLFc+QTjFRKBpgmRFEUFgxaYN3OteSSYkqxbtvr7Xlv8Hs46B3YcHYDq06uuvVFvfzKbB43dsUwbR14dsIh5XeWqW8WD8ADEOQUj/iVJJABpsnRK2UHwak3/Cfu4NaBf4T8A4AFexZwIu3ELS6oLx1Ro8D254ehurWCx9dhdu9AZ3GWX43P0oULgMLcHw5zPqWo+WU2axy7kommyWbTnUoGmCbGy8GLoGZB6BQdQc2C8HLwKlfm3k73Mt5vPPmWfOb8MofcwtzKL+jsC237g6JDadsf1bV50X7XluiG/RMBOCiFrDPOY5ByiLPJuTzwyW/sO5dK4BubGfXhLu5+bTNmmSXvjiQH2hVrKgPtoKg3KdWUipe9V6XJoXILc5m0fhJnMs4wtM1QFoYvRKfqKrmgBrnJ4ORTuhYtQNpFxIf+1iecfKHnH74f8e1Fd4x6lfzrgkrkXwfQvaUb0p1FPsE0QVXJ42uvt8deXzTGZfvF7Yz4fkTlqR1UtehJ5sbr5aVb/ygEGBUz72S8xLwuV4uDS0mmGawpHuSgvDuLDDB3qFRTKsdTj1u3E/MS+eTgJ9W6hubgWbqhgNZ5NEpBNjPOv0CE/r2inSiAwGTWZGa8O5AMMHcoD6MHDjqHMvsWxS1izak1Vb5Gap65NJeMgNTB8yF4OgqCEfr9bLD7Gw7kAQpP/m8fhxLSZWa8O4wMMHeoVFMq2ebSOUsls69f/e1VdibsrNI1PBztECXzIhXwcLKHse8jnFsgBHRTLxJj9xdcyeJUYjaPfb4HffFKBY4GuQb2nUAGmDvUje9bHHQOTPOfhkVYeOHnF8qMCK5MqurBXq0LhUJhr9aFVNUDFIX0Qa9Zy7iqeWx2+AeT/CxkmczWdzNZ+WaZGe8OYDMBRtM0Zs2aRWhoKOHh4Zw6darM8cWLFxMcHExISAjr1q0D4MKFC9xzzz2Eh4czePBgjh8veqfwwQcf0KNHD8LDwwkPD7ful0rplLI9RnmWPCZ3m8x4v/GYLCb+vO3Pt0y36e1i5IOWHzCg8BMWtvw/vF2MALgH3g9K0YtfATQXibyV9AwL9f9BRaPkvUxCaq586dvE2Uw39erVq1m7di3Lli0jOjqa+fPns2ZN0fuAq1evMnz4cGJjYzGZTAwYMIDY2FhmzpzJ+PHjuf/++9m0aROLFi1i9erVTJ48mdmzZxMUFFTlz29K3dRVYdEsDPh6gDW1g6POkd8e+Q2BYPbPs/n54s80c2zG8jHLae7UvNLrVJhHJjsR8V5nlOLswEqnoXB6O0LAReHDxIK/cxkfHA062jsXkpSWSevWbfn+qTB0MsdMk2Iz/zWjoqIYNWoUACEhIWV+2Pfs2UNYWBhGoxE3Nzf8/PyIi4vj/fffZ+zYsQCYzWbs7Yu6Xfft28f8+fMZMGAA8+fPb/gvcxtIy08rM+Exz5JHWn4aelXPvwf9m96+vbmWe42ZW2ZyOv10tfLIaPae5CoOCAG5iiPawysxO/oC0FZN4hfj83QggdxCC0fTVHqpJ3jq6j/546dRsmepibGZAJOZmYmbW+lALJ1Oh9lsrvCYi4sLGRkZeHt7YzAYOH78OHPmzOHVV18FYOLEiXz22Wds376dqKgoa5PqRhEREQQHBxMcHExSUlKFZZoqD6MHTobS3L2BPoF42ReN+rXX2/PR0I/wc/fjbMZZ7l9zP5MjJ99yCZQSqclXMGomFAWMmom0hHh0uckoxc0mg2LhB+NrrDS8jgu5bNb68jfznzidcI3L6bnsPpWMxSLnNDUFNhNgXF1dycrKsm5rmoZer6/wWFZWFu7u7gDs2LGD+++/ny+//JIuXboghOC5557D29sbOzs7xo4dy4EDByr8zJkzZxIbG0tsbCw+Pj719+VsUFp+GnnmopnQOkXHwvCFZZ5C3IxuvNH/Det2XHIcR5OPVunaXr6tOGHsgVmonDD2wLPd3dD87qLmUnGQcVdyCNUdY7fxLwxRDpCGC5k4MvTdLUxaEs1df99Ifr5ZvqO5zdlMgAkLCyMyMhKA6OhoAgICrMf69u3Lrl27MJlMZGRkEB8fj7+/Pzt27ODZZ59l48aNBAcHA0VPO/7+/mRnZyOEYPv27dV6F3On8LL3ItA3EJ2iI9A3sNycJU1ovL779TL7ntnxDGczzt7y2oqq0vVvv5D+VBzdXt6JotOhPFR2ZYPvHCYiAGfFxH/t/k1f5XdcFBMFGLCnAE1o3PufXUyMkAPzbmc285JX0zSefvpp4uLiEEKwdOlSIiMj8fPzY9y4cSxevJiIiAg0TWPevHlMmDCBnj17kp+fT/PmRS8hu3TpwqJFi/jyyy/56KOPMBqNDBs2jNdee+0Wn37nveSFm89ZSsxNZNh3w8qd42H04NN7PqWHd4/qfdjVI4jP+hf3HwGPrIIVE8qssZQr7Ph74TRWa4OK9wjcyCEDJ1QUds8bhqooZV8oSzbNZgJMY7sTA8zNJOcmM+S7IWX29Wvej5irMTjqHflo6Ef0a9GvytfTzBby3mqNo5ZLruqIw8SvUFaOtwYYIUqbT3FaR/5q/gvnRdEvjr5KPLGiC0HtvThwIZ2gdh6snBGCqpYPMhX2akmNxmaaSJJt8bD3wFHnaN121Dnyn2H/YUyHMeSac3lq61N8f+J7knOTq/R+JCXPTK/8RYzMf4ee+YtJdepgPSYEnNe8re9oeurOsNXuRV7XL8UBE3tENzQUDp+7ilnTiD2fRnIFuX/lXCfbIwOMVKEUUwq5ltI8MbmWXDILMpk/cD6Tuk6iUCvktd2vMeS7ITy+8fFb9jB5O9vRq50PZ9S2BLXzwsu1dB6UAGYV/qVMeYNiYYp+CzF2f+YRdQsqAhNG7CjEolmY8cVeLqflEn85A7PZQlJWPsnZ+ewrnuu0T851sglyDVCpQjdmwivZpyoqT/R4gpXHVlr370/cz+Xsy7R2aV3p9RRFYeWMkNLmCyDa9sd8IZp94i66t/SFovzkCAEZwhF3NRdXNY+37ZbyrPYD/zI/xgYtBIBDlzK4550N5GLAV8kmRXGjd1sPerf1YP+FNILaeeDpaCApK182lxqRfIKRKuTl4EWQb2nvW2/f3taepqyCrHLl/7r9r1zMvHjTa5YZlKcoKI+vR33+GB3n7OT9px8CRVe0PAoKoYXvFk01KG7lNFPT+cTwEdvtnidE+R1QyMUeV/J4WN3GIvVdYs+lsPDBu/lyWl+WT+vDI0tiCJ2/jYcX7eZahkl2dTcC+ZK3mHzJW54mNFLyipKGX5/AStM0+q3oh8liKlPexeDC62Gv42rnSpBvEDpdJRnyKpKdiHivCwoahQI+yH+AF42rrS9+UbC+ELYIhcNae94yT2Gv6AqAE3lYUBF6I/lmcLRTyTcLLMXvYXSqQvBNXg5L9UM+wUiVUhUVH0cffBx9blhSFrp5ditTdkibIWQVZjH759lM3zydXst7UWCpxjsQJx9oF4oZHftFV/Y1n2ydMFkSWUp+FeoUQaDuLN8ZX2e94WXClDhycMBEUXAByC3Q0CulvzstMgdNo5ABRqq2lLwUDiSVjo7e8sAWPhzyIRPummDdJxBsOLuh6hdVFJSp61Cfj6fjnF/4erp/ye7ri5RpNgH00F1gud0CdhieY5CyH4XSl81mi4XSjFhgr1fQNE02lRqQDDBStd34A6qqKoqi8ErfV8rs/8ev/+CLI19g0ao4r0hVUV2b4eNqj+Lsi9IuDFQ9SvNeZZZOuX50Xsn4mQ66RL6we49Ddk/yvO5r3MnCgg5Kh/aRU6ARMn87D3zyK/n5Zo5dybT2QMmgUz/kO5hi8h1M1SXlJjH0u6HW7e0PbsfH0YerOVcZ/v3wcuV7+fZiuv90BrYaiKpW43dayWoGjt7w31FwKQZa9kVo+XD1kDV03PhGpeRF8RYtiGXmUcSI7tZjKhY0VOtZqlKUzS+onTsrngwhLa9Q9jrVIRlgiskAU3VCCB7f+DgHEw8S6BvIslHLUBSFI0lHmBg50VqulVMr8i35JJuSATCoBnY9vAsnO6fKLl2565dO0SyIBe2gOJeNdTQwFQebq8KNZZbhrLEM5irXz7kSOGEiB3tAoYuvIyeTcvFv4cL8B+7G28WIr6u9DDa1IANMMRlgqufGeUya0Ji0bhJHU0tnXH879luyCrKYvmW6dZ9BNfDe4PcY0mZI7X5wLWY4vxv+94dyh0Tx/yk3NKUEEKt1ZpllFNu1XpgoysBXkhirfHiCu1u5sPixPng62nEyKRsvJzsZdKpBBphiMsDUTnJeMvd8dw8WUfS+xVHnyO5HdhN9OZo/bftTufL9W/bnpT4v0cm9U80/VAhYOgYuxkCrPnD1MJhzEOgB802bUAXo+MXSk1XaIHZogRRgd9OPKn2TA33aefDRI73QKeWTbUllyQBTTAaY2hFC8MSmJziYeJAunl1YMXoFOp2OhIwERv842lquhWMLss3ZZBVkoaAwuv1onu71NO1c29Xsg29oOpF8Ary7ID4fBZf2AMXB4YYnmtJ6Q74w8JVlEJ9YHiAFd8qGpIpCVKlOnvb836ReeLvYo1NVGXBuIANMMRlgaq+qzaZXo14lPiPeuk+n6BjTYQyjOoxiQMsB1XsRXGllNMhJgpWPIi7vBW4eaKD05fDvWlt+1PqzUwvijGhZrY/183Zk4cM98XI2kp5biKeT3R0deGSAKSYDTN2rqNn00/ifGPZ92TwzynVvQXSKjiUjllifaG61BO4taRpkXYMPuiOKx8hU5Wolgeia5sY27W6itEB+1QLIwLlG1ejewoW/jerMXb4u6HW6OybgyABTTAaYuldRs0lVVR784UGOZ5UuJfP3oL/z5r43K7xGoHcgi0cuJi4prvrTD65nLoALMbDjLcTF3WUOVeXJRlFAE5CgefOzFsgvoid7tK5kcWOPmECPBbN1DE7FujZz5OUx3fB0tEOv09G1uUvdPLnZGBlgiskAUz8qajZN2ziNfYn7APB39+eR7o8w77d5VbpezCMxJGQn4OfmV7MfSE2DrKuQmwJGN/i0H6KwKC3F9YHmZgGnhBBwTbizX7uLjVovIrVQLBgQ5cavFv2IGSkkHwOVBZ63xnUmI89Cj1ZudPBx5kJKHh6OBro0dyUtz3xbjs+RAaaYDDAN4/pmk07RsfGBjYz7cRx5lrwqna9DhwULTnonoiZFoVdrmXHEYobEY7B0FKKCWeKV9UTdqOSnKF/oOSNasl/zY6/owhHRnjOiVfHgvnJnVeHKpe7yceDfD/bETqejawtXa4C15Sx+MsAUkwGmYVzfbAr0DWRu37n88ac/Wo/bYUcBVZuQOLbDWDq6dWRy18nkWIrWePK09yQtP63CPMM3ZTEX9UB5dIRLsdC8J7zvB2aTtXu6qsHmekUvjlUuCS+OaO3Zq3UlXrTllGhFMu7VuFJ5ozo7o6Hj0JVcrmUV4uWgY9qAdkzt155DV7Lx83XC29me1NzGG50sA0wxGWAazvXNJiEEYV+HkV2YjaPekV0P7mLQd4PIMefc+kIVMGCgkEJ6+fTi81Gfk56fjofRo3ZBx70DRIRDyrEywQYot10d+ULPVeHBGdGCw1pHDolOnBItSRC+xfOoKlPdMFc06bB3Szs8nO1p4+GESRP4+TrTxt2JxJwChnT2xcfFgeOJWQgh0CkKnZu51LppJgNMMRlgGo9ZM3Mm44z1vYpZM3Mi9QQJWQm8tOslay9UdTnpnMix5FiDjr+7P3/s9kfGtBvDxdyL1XuPUzLext4dIoZC4mHw6QkBf4TtrwJlU4ZWPwQUn1c84jgPIynClcvCizOiBce0NsSLdpwTzUnGrYL3PPWriwcUWsCiwV3NXejZ2oN2vm4EtHHj0PkMMkz5TB1wV7nzZIApJgOM7bm+OdXCsQVXcq/UONhUxIiR54Ofp4dnDzLNmegVPW1c2pCQk0Brp9Zcyr1Eb5/eZBRm4GH0INVUlNPT2+iJkpdSNLhPUYqedK4dBXM+XNwL6Zdg73/Kfx9qFnRupAkowEC2cCANZ5KEO5eFFwnCm/OiGWe0FpylOZm41ODq19ey7FwtEDhiIrd4W8FCaxJRgAS8OLNgfLmr2UyAKVkX6dChQxiNRpYsWYKfn5/1+OLFi1m0aBF6vZ6///3v/OEPfyA5OZlHHnmEvLw8WrZsydKlS3F0dKyw7K3IAGObSppTnkZPpm2exr5r+xq7SnjrvHki4Anc7d2xCAsmi4lmDs1IMiXR3bM7qk5F1QTe2Ul4F5goTD5NTOoJsk5vpWXuNa7o9bQ1m3HXNM7p9ThrGuf1enIUcBYCH4uGkxA0t1iwA+yEwA6umwNeNRpQiIpZGMgRdpzXOaFpDhRqjqQLZ1Jw5SqupAlXkoU7ScKDDJzIFo5kY48Juyp/4hTfmbzx/KVy+20mwKxevZq1a9eybNkyoqOjmT9/PmvWrAHg6tWrDB8+nNjYWEwmEwMGDCA2NpYXX3yR3r178/jjj7NgwQKMRiOTJk2qsKzRaLzp58sAY/tKUngKIXhp50vWrm5bVp2X1tZZmpXQCYGrpuFmNuNhtuChgFuhmWZmCy6ahhPgJAQumsBJaNbQUJOnppJkgAUo5Cs6CtGRj44C9MV/NlAgDBRgoAA7Pm2ZweHHfy93HZtZVSAqKopRo0YBEBISUuaHfc+ePYSFhWE0GjEajfj5+REXF0dUVBTz5hWNnxg9ejTz5s2jU6dOFZbt06dPo3wvqe6UpPAE+HzU59Z8wa52ruy9upeF+xZyMv0kAHr0mDE3Wl1LVDm4QHEkqDwcWBSFNFUlTa/nXG0rVicsxf8zUVm9bSbAZGZm4ubmZt3W6XSYzWb0en25Yy4uLmRkZJTZX9G+6/dXJCIigoiICACSkpLq42tJ9eT6YAMwoPUA+rfqT1JuEun56fi5+TF9y3QOJh6km0c3nun1DB8f/JjfU36nm0c3evv05qsTXzXiN7gNCIFCadNMLbddlOpCBTIqGWFtMwHG1dWVrKzSgU6apqHX6ys8lpWVhbu7u3W/g4NDuX03lq3IzJkzmTlzJlDURJJub6qi0sypGc2cmgHw+cjPy4wiDmkVYt0GiE+P52DSQQI8Arjvrvv4Jv4bjmcep4WxBY4GR05nn27Mr1M/RLk/3LK49bX6jV3V1z9xVfKmxWYCTFhYGD/99BMPPfQQ0dHRBAQEWI/17duXV155BZPJRH5+PvHx8fj7+xMWFkZkZCSPP/44GzZsYODAgZWWle48qqLi7eBd6fbSUUvLBKAJXSaUjs9BkJSbRJopDReDCxeyL9DBpQMWLJzNOIumaegUHU56J1759RUu5FywTtpUUXnirie4knuFvII8VFWlhXMLHPWOfH38azLJxAEHurp15S63u0jNT8XZzrmoy1wU/XJNzEuksLAQM2ayTFkYdAYKzAWgwun80sDngw8CQT75aGjo0WMp/qeon8eChoaKigEDZsWM6SZNmhrJASis8JDNvOQt6UWKi4tDCMHSpUuJjIzEz8+PcePGsXjxYiIiItA0jXnz5jFhwgSuXbvG1KlTycrKwtvbmxUrVuDk5FRh2VuRL3mlmirp6XKzc+Ns5tmbjq+5cW5WbT6vNtdoKDYTYBqbDDCSVPea3vxwSZJshgwwkiTVGxlgJEmqNzLASJJUb2SAkSSp3sgAI0lSvZEBRpKkeiPHwRRzdXWlc+fO1u2kpCR8fHxuckbdysjIKDOHqr6vUZWytypT2fGq7q+o3O1236t7fm3ve3WPVWVfXdxzb29vNm7cWP6AkIQQQsyYMaPMdlBQUKN+fn1foyplb1WmsuNV3V9Rudvtvlf3/Nre9+oeq8q++rznsolU7N57773tP78616hK2VuVqex4Vfc39j2vizpU9/za3vfqHqvqvvoim0iVkFMHGoe87w2vPu+5fIKpREkaB6lhyfve8OrznssnGEmS6o18gpEkqd7IACNJUr2RAUaSpHojA0wV/fbbb0ydOpWpU6eSnp7e2NW5o2zfvp0nn3yysatxx9i2bRszZszg0Ucf5dChQ7W6lgwwVRQREcGiRYuYPn0633zzTWNX545x6tQpDhw4gMlkauyq3DFyc3OJiIhgzpw5bN68uVbXkgGmiiwWC/b29rRo0YIrV640dnXuGH5+frzwwguNXY07yr333ktubi4fffQRU6dOrdW1ZICpIkdHR/Lz87ly5QrNmzdv7OpIUr1JTk7mmWee4fXXX8fX17dW15IBBoiJiSE8PBwoWt1g1qxZhIaGEh4ezqlTp4CiwUh/+tOfWLRoEZMnT27E2jYdVbnvUt2qyj1//vnnuXLlCnPnzuX777+v3QfW2yyn28Q777wj/P39Rb9+/YQQQqxatUpMnTpVCCHE7t27xbhx4xqxdk2XvO8NrzHu+R3/BNOpUydWr15t3b7ZGtlS3ZH3veE1xj2/4wPMhAkTMBgM1u3K1siW6pa87w2vMe75HR9gbnSzNbKl+iPve8NriHsuA8wNSta7BsqtkS3VH3nfG15D3HP5K+IG48ePZ8uWLfTv39+6RrZU/+R9b3gNcc9lugZJkuqNbCJJklRvZICRJKneyAAjSVK9kQFGkqR6IwOMJEn1RgYYSZLqjQwwkiTVGxlgJEmqNzLASJJUb2SAkWxaTEwMoaGhODg44OHhwRtvvNHYVZKqQQYYyWZt3bqVsWPHMn36dA4dOsRLL73EP//5T/bv39/YVZOqSM5FkmxSQUEBXbp04Z///CdPPPGEdX+LFi148cUXef755xuxdlJVyScYySb98ssvpKenl8t/bDAYMBqNjVQrqbpkgJFs0vbt2wkICCiTgS0xMZFLly7Ru3fvRqyZVB0ywEg26cCBA+Tn55fZ98knn9CxY0dCQkIaqVZSdcmEU5JNKgkw//3vfxk4cCBr1qzhnXfeYcuWLSiK0tjVk6pIPsFINufy5cskJiayYsUKPv74YwICAlixYgVr1qxhwIABjV09qRpkL5JkcyIjI3n00UdJS0tr7KpItSSfYCSbc+DAAe6+++7GroZUB2SAkWyODDBNh2wiSZJUb+QTjCRJ9UYGGEmS6o0MMJIk1RsZYCRJqjcywEiSVG9kgJEkqd7IACNJUr35/9ALvdvQa3DzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def eD0Dt_norm_test(r,t):\n",
    "    \"\"\" Revamped version of that function ... \"\"\"\n",
    "    return(((r+2)*(r+10))/((r**3 + 15*r**2 + 48*r + 48)*np.exp(t*(r+2)/2) - 4.0))\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(4,4))\n",
    "for x in tqdm([0,100, 1000]):\n",
    "    ld_raw_df_filt = ld_raw_test[(ld_raw_test['ta'] == x) & (ld_raw_test.maf == 0.00)]\n",
    "    ax.scatter(ld_raw_df_filt.rec_dist*1e4, ld_raw_df_filt.ed0dt, s=5)\n",
    "    z = x/(1e4)\n",
    "    # Dealing with some rescaling between sims vs. theory (factors of two here ... )\n",
    "    cur_d0dt_norm = np.array([eD0Dt_norm_test(y,z/2) for y in np.unique(ld_raw_df_filt.rec_dist*1e4*4)])\n",
    "    ax.plot(np.unique(ld_raw_df_filt.rec_dist*1e4), cur_d0dt_norm, lw=2, label=r'$t_a = %d$' % int(z*1e4))\n",
    "    \n",
    "ax.legend()\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel(r'$\\rho$', fontsize=14)\n",
    "ax.set_ylabel(r'$\\sigma^2_t$', fontsize=14)\n",
    "debox(ax);\n",
    "ax.set_xlim(0.9, 1.5e2)\n",
    "ax.set_ylim(-0.01,0.175)\n",
    "plt.savefig(main_figdir + 'eD0Dt_theory_sims.pdf', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}