kaspermunch/pool_nielsen.ipynb

## pool_nielsen.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "47f5be28-140e-4cad-aa03-b70004fe054f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib_inline.backend_inline import set_matplotlib_formats\n",
    "set_matplotlib_formats('retina', 'png')\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "sns.set_style(\"ticks\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "be04d271-e862-4f9c-8e08-fe282e08f981",
   "metadata": {},
   "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>years</th>\n",
       "      <th>Ne</th>\n",
       "      <th>population</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>204558</td>\n",
       "      <td>nonAfr</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1015</td>\n",
       "      <td>204558</td>\n",
       "      <td>nonAfr</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1090</td>\n",
       "      <td>204558</td>\n",
       "      <td>nonAfr</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1171</td>\n",
       "      <td>204558</td>\n",
       "      <td>nonAfr</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1257</td>\n",
       "      <td>204558</td>\n",
       "      <td>nonAfr</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   years      Ne population\n",
       "0      0  204558     nonAfr\n",
       "1   1015  204558     nonAfr\n",
       "2   1090  204558     nonAfr\n",
       "3   1171  204558     nonAfr\n",
       "4   1257  204558     nonAfr"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sampledemog_raw_data = pd.read_csv('smcpp_example_data.csv')\n",
    "sampledemog_data = pd.DataFrame(dict(years=sampledemog_raw_data.x,\n",
    "                               Ne=sampledemog_raw_data.y,\n",
    "                               population=sampledemog_raw_data.label\n",
    "                              ))\n",
    "sampledemog_data['years'] = sampledemog_data.years.round(0).astype(int)\n",
    "sampledemog_data['Ne'] = sampledemog_data.Ne.round(0).astype(int)\n",
    "sampledemog_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b9290902-6cef-463d-a96c-4abc572b7a91",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHsAAAN4CAYAAABTYuTOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAB7CAAAewgFu0HU+AAB6dklEQVR4nO3deZyddXk3/uvMmtkyIYAhCWtlUxEDRnYLCBatFAU3MCJE+SkgKqKCS5VHigV97KO2aGvVjgWxgLILSlkCFKzCKAIKRMIiYYtkny2Z7fz+iBxy7jOZJGTmLPf9fr9efXXmzD0zV1C+Zj5zfa8rl8/n8wEAAABAKtRVugAAAAAAJo6wBwAAACBFhD0AAAAAKSLsAQAAAEgRYQ8AAABAigh7AAAAAFJE2AMAAACQIsIeAAAAgBQR9gAAAACkiLAHAAAAIEWEPQAAAAApIuwBAAAASBFhDwAAAECKCHsAAAAAUkTYAwAAAJAiwh4AAACAFBH2AAAAAKRIQ6ULoDr19PTEwoULC+/vscce0dHRUcGKAAAAgE0h7GFMCxcujHnz5hXev/TSS2Pu3LkVrAgAAADYFK5xAQAAAKSIsAcAAAAgRYQ9AAAAACki7AEAAABIEWEPAAAAQIoIewAAAABSRNgDAAAAkCLCHgAAAIAUEfYAAAAApIiwBwAAACBFhD0AAAAAKSLsAQAAAEgRYQ8AAABAigh7AAAAAFJE2AMAAACQIsIeAAAAgBQR9gAAAACkiLAHAAAAIEWEPQAAAAApIuwBAAAASBFhDwAAAECKCHsAAAAAUkTYAwAAAJAiwh4AAACAFBH2AAAAAKSIsAcAAAAgRYQ9AAAAACki7AEAAABIkYZKFwAAAAAwlre97W2xaNGiiIg46qij4p//+Z836/N/9atfxfe///146KGHYvXq1bHVVlvF3Llz4xvf+MZklFs1hD0AAABA1bn//vsLQU9ExK233hovvPBCbLvttpv0+XfeeWeceuqpMTIyUnjtz3/+c9TX1094rdVG2AMAAABUnSuvvDIiImbPnh3PPvtsDA8Px09+8pM4/fTTN+nzv/3tbxeCnlNPPTXe+MY3Rj6fj+nTp09azdVC2AMUjIzm439+90zc9KsnY8ny/kqXU/W2m94WH/y718SuO0yrdCkAAJAqa9asiRtvvDEiIo444oi4//774/7774+f/OQnceqpp0Zd3cZHEP/xj3+MiIjXve518clPfnJS6602wh4gRkfzcfcDz8Z//fcjsXhJb6XLqRkvrBiIz37nrjj7/XNjv9dsV+lyAAAgNX7xi19ET09PRETsv//+MXPmzLj//vvj2WefjTvuuCMOP/zwjX6NgYGBiIjYcccdJ7XWamQbF2TYiyHPx/9pQXztkm5Bz8uwdnAkvtL167jxl09UuhQAAEiNq666KiIiWlpa4uCDD463ve1thW6eyy67bJO+Rj6fj4iIhobs9blk708MRD6fj1/9/vn4r/9+JJ54dnWly6l5o/mIf73ygfjz8v74wN++OurqcpUuCQCADNhjjz0iIuITn/hEnH766XHHHXfEFVdcEffff3+sXLkypk2bFvvuu28cf/zxcdBBB437tVatWhWXX3553H777fHYY49FX19fTJ06NXbfffd485vfHO9617uiubl50uuIiFi8eHHcc889ERHx13/919HS0hItLS1xwAEHxC9/+cu4884745lnnonZs2eXfO5VV10Vn/vc54peu/rqq+Pqq6+OiIj99tsvLrnkknj66afjiCOOiIiIr33ta7HnnnvGP/7jP8b9998fDQ0NseOOO8YnPvGJOPTQQzdabzUS9rBJHntmZTRNe6HSZTABVvUMxpW3PxqPPb1qg8/sMmtqvOPQXWNqW1MZK6std973dCz4zdNFr125YFG8sHIgzjx+n2hsSP+EfwCAjRkZzUdv/2Clyyi79tamqC/jLwDz+Xx86Utfissvv7zo9RdeeCFuuummuOmmm+KDH/xgnHPOOWN+/q233hqf//znY+XKlUWvL1u2LP73f/83/vd//ze+//3vx7e+9a3Ye++9J62OF1111VWFrpxjjjmm8Prb3/72+OUvfxmjo6NxxRVXTNgcnscffzzOO++86O196abDH/7wh5q+/pXLv/hPENbT3d0d8+bNK7y/w0GnRcv0XSpYEeWw03Yd8b6j9owD9pqpO2Uj8vl8/PimhXHZzQtLPrbXK7eOL5y8X7S3CssAgOy66/5n4rtXPRgre9dWupSym9beHB857rVxyOtKO08m0osdNdtuu2288MILscMOO8TJJ58cr3nNa2JgYCB+8YtfxBVXXFEITv7zP/8zDjjggKKvcfvtt8dHP/rRGB4ejrq6unj7298ef/M3fxPbbLNNPPvss3HttdfGbbfdFhHrrlRddtllseeee054HS8aHR2NI444Ip599tmYPn163HnnndHY2BgR62bwHHzwwdHX1xfbbrttLFiwoPCxF61cuTKee+65iIh4xzveERERhx9+eHziE5+IiIjW1tbYaaedijp7XlzF/pGPfCQOOeSQeOGFF+IPf/hDfOpTn9qc/ziqis4eIHaY0R4n/M2ecfDes4Q8myiXy8W8t+wZ227VEt/+6f0xOvpSbv77x5bF2RfdFf/nlAPiFdNbK1glAEDlXHTF76JvzXCly6iIlb1r46IrfjfpYc+LXnjhhZgzZ078x3/8R7S1tRVeP+igg2LWrFnxjW98IyLWdcysH7L09fXF5z//+RgeHo6Ghob49re/HYcddljh43vvvXe85S1viUsvvTTOO++8GBgYiLPOOit+9rOfjbkN6+XWsb5f/vKX8eyzz0ZExNFHH10U5rS0tMRRRx0VV111Vbzwwgtxyy23xFvf+taiz582bVpMmzat5LVXvepVG/znNzIyEp/5zGfilFNOKbz2lre8ZYPP1wIDmiHDZm/bFp+e9/r4l0+/Kd44Z7ag52X4m/13ii99aP9oaS6+trV4SU985l/ujMef2fB1OQAAmCif//zniwKWF73nPe8pvP3iKvIXXXXVVbFs2bKIiDjllFOKgp71zZs3L/72b/82IiIee+yxQqfPRNWxviuvvLLw9rHHHlvy8fVf29RBzRvT0NAQJ5xwwoR8rWqhs4dN0tbSGB3trqSkxcyt2+KtB+0Sh+4zO+rrZb5b6vV7zogLTj8kvvz9X8WKnpfalJevXhuf/uc7ozOrs49yuXjl7M447Z17x9adLZWuBgAoszPeMyfz17jKZerUqRucpTN9+vRob2+P3t7e6O/vL/rYnXfeWXj7fe9737jfY968eXHjjTcWPu/II4+csDpetGrVqrjlllsiImL33XePV7/61SXPvOENb4jtt98+nn766fj1r38dTzzxROyyy5aNHNl1113HDKhqmbCHTfLFD+4fc+fOrXQZULVeuf20+PrH/zr+z/f/t2iF/dDwaCxdtaaClVXW0pUDUVeXi8+fvF+lSwEAyuyQ182OA187y4DmMpg9e3bkchv+fq2trdHb2xvDw8XX6h599NGIiJg5c2bMmDFj3O+x9957R319fYyMjGywM+fl1vGin/3sZzE4uO6/L2N19USsG6fwjne8Iy666KLI5/Nx2WWXlWzf2lwzZ87cos+vRsKeDOrq6oqurq5xn3nxXzBg071iemt87Yw3xld+eE/8/rFllS6naty38M8xMpov6194AIDqUF+Xi872sdd1M3FaWsbvot5QALNixYqIiNh66603+j2amppi6tSpsWLFisLnTVQdL1r/CtdXv/rV+OpXv7rRuq655po466yzNrgWflOkrasnwsyeTOrt7Y0lS5aM+38b+pcXGF97a1Oc9+ED4/DXb1/pUqrGmsGRePrPPZUuAwCAhM1dzv3i82MNZ95SjzzySPzhD3/Y7M9buXJl/PznP9+i772xEKoW6ezJoPb29o226A0ODgp84GVqbKiPs973+jj2sF1j+epsXuH6lyt+F8vWu7726FMrY6ftplawIgAAkjo7O+PPf/5zYUjzeNauXRs9Pet+gZfcdjURfvrTnxbe/tCHPrTB2T8veuihh+K73/1uRKwb1PzimnXWEfZk0Pz582P+/PnjPtPd3R3z5s0rU0WQTrvM6oxdZnVWuoyK2GOnreKXDzxXeP/RxSviyP12rGBFAAAk7bHHHvHnP/85nnvuuViyZMm4TQEPPPBAjIyMRETEX/3VX01oHYODg3H99ddHRMSUKVPi9NNPj/b29nE/57DDDosf//jH0dPTE/fdd1888sgjseeee05oXbXMNS4AJtxuO2xV9P6ji1dWphAAADbojW98Y+HtH//4x+M+u/7HDz744Amt47bbbouVK1dGRMQRRxyx0aAnYl0o9OI6+IiI//qv/5rQmmqdsAeACbfbDtOK3n/i2dUxNDxamWIAABjTO9/5zsKVrB/84AdFq9jXd/nllxfWru+8887x5je/eULrWH8w89vf/vZN/rx3vetdhbevv/766Ovrm9C6apmwB4AJt+v204reHx4ZjSefW1WZYgAAGFN7e3tccMEFkcvlYmhoKE499dT43Oc+FwsWLIgHHnggbrrppvjoRz8aX/rSlyIiorm5Ob71rW9FY2PjhNWwZMmSuPvuuyMiYvr06ZvVNbT33nvHbrvtFhERfX19hatgCHsAmARtLY0xe9vi9ltXuQAAqs+b3vSmuOiii6KjoyNGRkbiqquuilNPPTXe/e53x8c//vG45ZZbImJdR89//dd/TfhcnKuvvrowC+htb3tbNDRs3mjh4447rvD2ZZddNqG11TJhDwCTInmV69GnVlakDgAAxnfkkUfGLbfcEmeeeWbss88+MW3atGhsbIxZs2bFIYccEl/96lfj2muvjde85jUT/r2vvvrqwtvHHHPMZn/+29/+9kKn0cMPPxz333//hNVWy3L5fD5f6SKoPsltXJdeemnMnTu3ghUBtea6Ox+L7137+8L7O23XERd95k0VrAgAALJBZw8AkyK5kWvxkp5Ys3a4QtUAAEB2CHsAmBS7zJ4adXW5wvuj+YjHnjGkGQAAJpuwB4BJMaWpIXbarqPotUcXr6hQNQAAkB3CHgAmTfIqlyHNAAAw+YQ9AEyako1c1q8DAMCkE/YAMGmSYc9zy/qip3+wMsUAAEBGCHsAmDQ7zZwaTQ3F/1OjuwcAACaXsAeASdNQXxe7zO4ses2QZgAAmFzCHgAmVcncHkOaAQBgUgl7AJhUJRu5XOMCAIBJJewBYFIlO3uWr14Ty1YNVKYYAADIAGEPAJNq9rbt0Tqloeg13T0AADB5hD0ATKq6ulzsuv20oteEPQAAMHmEPQBMutIhzTZyAQDAZBH2ADDpkkOaFz29MvL5fIWqAQCAdBP2ADDpkp09Pf1D8fyy/soUAwAAKSfsAWDSbbtVS3S2NxW99uhiV7kAAGAyCHsAmHS5XK7kKpchzQAAMDmEPQCURcmQZmEPAABMCmEPAGWRDHsee3pljIwa0gwAABNN2ANAWSSvca0ZHImnl/RUqBoAAEgvYQ8AZTGtozm23aql6DVDmgEAYOIJewAom+RVrj+a2wMAABNO2ANA2djIBQAAk0/YA0DZJDt7nnx2VQwNj1SmGAAASClhDwBls+v204reHx7JxxPPrq5MMQAAkFLCHgDKpq2lMWZv2170mqtcAAAwsYQ9AJTVbjtOK3rfRi4AAJhYwh4Ayio5t0dnDwAATCxhDwBltXtiI9fTS3piYO1whaoBAID0EfYAUFa7zO6M+rpc4f3RfMRjT6+sXEEAAJAywh4Ayqq5sT522m5q0WuucgEAwMRpqHQBAGTPbjtOi8efXVV4/4pb/hg33/OnClZUzXKx6/adccrbXxtT25oqXQwAADVA2ANA2e22w7S46VcvhTu9A0PROzBUwYqq2+IlPTE6GvHp97++0qUAAFADXOMCoOz23Gl6pUuoOX94fGmlSwAAoEYIewAou51mTo0DXzuz0mXUlOU9a2N0NF/pMgAAqAGucQFQEeecODce+dOKWNGzptKlVKX+NcPxL1f8rvD+6Gg+VvWuja2mTqlcUQAA1ARhDwAVUV9fF6/5q60rXUbVGhnNx7d/en9RN8+y1WuEPQAAbJRrXABQherrcrFVR3PRa8tX6YICAGDjhD0AUKWmJ7p4lq0W9gAAsHHCHgCoUlt3JsKeVQMVqgQAgFoi7AGAKpXs7HGNCwCATSHsAYAqNT3Z2eMaFwAAm0DYAwBVauupLUXv6+wBAGBTCHsAoEqVdPYIewAA2ATCHgCoUskBzT39gzE0PFKhagAAqBXCHgCoUlsnBjRH6O4BAGDjhD0AUKXaWhqjqbG+6LXlhjQDALARwh4AqFK5XK6ku0fYAwDAxgh7AKCKGdIMAMDmEvYAQBUr6ewR9gAAsBHCHgCoYjp7AADYXMIeAKhiyfXrZvYAALAxwh4AqGJbT20pen/ZqoEKVQIAQK0Q9gBAFUte41q+ek3k8/kKVQMAQC0Q9gBAFZueGNC8ZnAk+tcMV6gaAABqgbAHAKpYsrMnwtweAADGJ+wBgCrW3Fgf7S2NRa+Z2wMAwHiEPQBQ5WzkAgBgcwh7AKDKJef2LFsl7AEAYMOEPQBQ5bbuLF6/vlzYAwDAOIQ9AFDlkkOal7nGBQDAOIQ9AFDlSmb26OwBAGAcwh4AqHIlM3t09gAAMA5hDwBUuWRnz4rVa2J0NF+hagAAqHbCHgCocsnOnpHRfKzqW1uhagAAqHbCHgCoctM6pkRdrvg169cBANgQYQ8AVLn6ulxM60gMaTa3BwCADRD2AEANKFm/rrMHAIANEPYAQA3Yeqr16wAAbBphDwDUgNLOnoEKVQIAQLUT9gBADUiuXzezBwCADRH2AEANKLnGJewBAGADhD0AUAOmd7YUvW9AMwAAGyLsAYAakOzsWd03GEPDIxWqBgCAaibsAYAakJzZExGxfPXaClQCAEC1E/YAQA1oa2mMpobi/9m2fh0AgLEIewCgBuRyudg6ObdntfXrAACUEvYAQI2Ynly/rrMHAIAxCHsAoEZMTwxptpELAICxCHsAoEYkhzQvXy3sAQCglLAHAGqEzh4AADaFsAcAakRpZ48BzQAAlBL2AECNGKuzJ5/PV6gaAACqlbAHAGpEcvX6msGRGFg7XKFqAACoVsIeAKgRydXrEeb2AABQStgDADWiubE+2lsai15bLuwBACBB2AMANSTZ3bPM+nUAABKEPQBQQ7YuGdJsIxcAAMWEPQBQQ5KdPct19gAAkNBQ6QJ4ySWXXBLnn3/+Bj9+6aWXxty5c8tYEQDVJrmRy4BmAACShD1V5KGHHoqIiJNOOik6OjpKPj5r1qxylwRAlZk+VWcPAADjE/ZUkYcffjiam5vjnHPOifr6+kqXA0AVSoY9OnsAAEgys6dKDA4OxqJFi2L33XcX9ACwQVsnZvasWL0mRkfzFaoGAIBqJOypEo8++mgMDQ3Fq171qkqXAkAVS4Y9I6P5WNW3tkLVAABQjYQ9VeLFeT25XC7OOuus+Ou//uvYe++945hjjolLLrkkRkdHK1whANVgWntz1OWKX1vuKhcAAOsR9lSJhx9+OCIiLr/88njhhRfi6KOPjre85S2xZMmSOP/88+PMM88U+AAQ9fV1Ma2juei1ZYY0AwCwHgOaq0Qul4tZs2bFxz/+8Tj22GMLry9dujROPvnkuOmmm+KKK66I448/voJVAlANpne2xPLVL13d0tkDAMD6dPZUiS9+8YuxYMGCoqAnImKbbbaJz372sxERcc0111SgMgCqzdbWrwMAMA5hTw143eteFxERixcvrnAlAFSD6Z3WrwMAsGHCniowNDQUDz74YNx7771jfry/vz8iIpqbm8f8OADZorMHAIDxmNlTBUZHR+OEE06IkZGRuPvuu2P69OlFH7/nnnsiImLOnDkVqA6AapNcv75s1UCFKgEAoBrp7KkCzc3N8eY3vzlGR0fjwgsvLNq69dRTT8XXv/71qKuri5NOOqmCVQJQLaZPbSl6X2cPAADr09lTJT772c/GAw88ENdee20sXLgwDjzwwFi6dGnceuut0d/fH5/73OcKs3sAyLZkZ8+q3sEYGh6Jxob6ClUEAEA10dmzEV/84hdjjz32iG984xsbfXZ0dDSuvPLKOPHEE2O//faLvfbaKw4//PD41Kc+tcF5PC+aMWNGXHnllTF//vzo6+uLH/3oR3HnnXfG61//+vjhD38YJ5988gT9iQCodckBzRERK9ZbxQ4AQLbp7BnHzTffHFdcccUmPdvT0xOnn356Yb7Oi5599tl49tln44YbboiTTz65sEZ9LNOmTYvPfvaz4z4DAO0tjdHYUBdDwy9d+122ak28YnprBasCAKBaCHs24I477ohPfvKTm/RsPp+PM888sxD0HHLIIXHCCSfENttsEw8//HB873vfi2eeeSa6urpi+vTp8eEPf3gySwcg5XK5XGzdOSWeX9ZfeM3cHgAAXiTsGcMPf/jD+PrXvx5DQ0Ob9Pz1118fd911V0REHHfccXHBBRcUPjZnzpx461vfGvPmzYtFixbFRRddFMccc0xst912k1L7+l7sKno5Fi5cOMHVADCRpk8tDnts5AIA4EXCnvU8+eSTceGFF8aCBQsiIqK+vj5GRkY2+nldXV0REdHe3h7nnHNOycenTZsWX/7yl2PevHmxdu3auPjii+Pss8+e2OLHcOWVV8ZFF1006d8HgPLbutNGLgAAxmZA819ceumlcfTRRxeCnl133TW+/OUvb/TzFi9eHA899FBERBx++OExbdq0MZ+bO3du7LLLLhER8Ytf/GJiigYgs6ZPLR7SvEzYAwDAXwh7/uLBBx+MoaGhaGpqio985CNx1VVXxY477rjRz/vNb35TePuAAw4Y99n99tsvIiKeeeaZeOqpp7asYAAyLbl+ffkqYQ8AAOu4xvUXzc3N8e53vztOO+20mD179iZ/3qJFiwpv77zzzuM+u8MOOxTefvTRRzcpTNoS73znO+PAAw98WZ+7cOHCOO+88ya4IgAmSklnj7AHAIC/EPb8xbnnnht1dZvf6PT8888X3p41a9a4z86cOXPMz5sss2bN2mhNANSmks6e1QY0AwCwjmtcf/Fygp6IiFWrVhXebmtrG/fZ1tbWwts9PT0v6/sBQETE9ETYM7B2JPrXbNoWSQAA0k3Ys4UGBwcLb0+ZMmWcJ4s/vv7nAcDmSl7jinCVCwCAdYQ9W6i+vr7wdi6XG/fZfD5fePvldhIBQETElKaGaGtpLHrN+nUAACKEPVts/atZa9aM/5fstWvXFt5uamqatJoAyAZDmgEAGIuwZwutP6dnYGD84Zj9/f2Ftzs7OyetJgCyoXRIs7AHAABhzxZbf037c889N+6z6398xowZk1YTANlQ2tljIxcAAMKeLbbbbrsV3n7qqafGfXbx4sWFt3fddddJqwmAbNDZAwDAWBoqXUCtmzNnTuRyucjn89Hd3R3HHHPMBp+95557IiJi5syZsf3225erRABSautEZ8/9f3whPvvtuypUTeXU1+Vir7/aOt5z5O5RX+/3WAAAwp4tNHPmzJgzZ07cd999cdNNN8XZZ58d7e3tJc91d3fHE088ERERRx11VLnLBCCFpne2FL3ft2Y4/vD4sgpVU1kPLFoadfW5eO+Re1S6FACAivPrrwlw4oknRkTEypUr49xzz43R0dGij69atSrOPffciIhobGyM97///WWvEYD02W7r1o0/lCEPPLq00iUAAFQFYc8EeNvb3haHHHJIRET87Gc/ixNPPDFuuumm+N3vfheXX355HHvssbFo0aKIiPjYxz4WO+ywQyXLBSAldp45NV7zV1tXuoyq0ds/VOkSAACqgmtcE+Rb3/pWnHrqqXHvvfdGd3d3dHd3lzxz8sknx4c//OEKVAdAGuVyufjyhw+M3z7y51jZk73hzH96viduuPuJwvs9A4MVrAYAoHoIeyZIe3t7XHzxxXHNNdfEddddF4888kj09PTEVlttFfvss0/MmzcvDjjggEqXCUDKNDfWx4GvnVnpMirioSeWFYU9OnsAANYR9oxj//33j4ULF27y83V1dXHcccfFcccdN4lVbbmurq7o6uoa95nBQb8dBaC6tbc0Fr0/sHY4hkdGo8FGLgAg44Q9GdTb2xtLliypdBkAsEU6WptKXusbGIrO9uYKVAMAUD2EPRnU3t4eM2bMGPeZwcHBWLFiRZkqAoDN197aWPJaT/+gsAcAyDxhTwbNnz8/5s+fP+4z3d3dMW/evDJVBACbr7GhPpqb6mPt4Ejhtd4Bc3sAAFxqBwBqVnJujyHNAADCHgCghiXn9vT0WzAAACDsAQBqVnJuj84eAABhDwBQw0qvcensAQAQ9gAANSt5jcuAZgAAYQ8AUMPaEp09ZvYAAAh7AIAaVjqgWWcPAICwBwCoWR2JAc19rnEBAAh7AIDa1d5i9ToAQJKwBwCoWSWr13X2AAAIewCA2lUS9vQPRj6fr1A1AADVQdgDANSs5IDm4ZF8rBkcqVA1AADVQdgDANSs9kTYExHRayMXAJBxwh4AoGa1NjdELlf8Wu+AIc0AQLYJewCAmlVXl4v2luTcHp09AEC2CXsAgJpm/ToAQLGGShdA+XV1dUVXV9e4zwwO+osyALWhvbUxYtlL71u/DgBknbAng3p7e2PJkiWVLgMAJkTpNS6/sAAAsk3Yk0Ht7e0xY8aMcZ8ZHByMFStWlKkiAHj5kuvXe8zsAQAyTtiTQfPnz4/58+eP+0x3d3fMmzevTBUBwMvX3pro7HGNCwDIOAOaAYCa1l7S2eMaFwCQbcIeAKCmdSQ6e/pc4wIAMk7YAwDUtOSA5p4BnT0AQLYJewCAmlZ6jUtnDwCQbcIeAKCmJbdx9ZnZAwBknLAHAKhpyWtcfWuGY2Q0X6FqAAAqT9gDANS05Or1iIg+69cBgAwT9gAANS05syciotdVLgAgw4Q9AEBNa26sj6aG4r/S9Ah7AIAME/YAADUv2d3T6xoXAJBhwh4AoOYl5/ZYvw4AZJmwBwCoedavAwC8RNgDANS85Pr1Hte4AIAME/YAADUveY2r1zUuACDDhD0AQM1rbym+xmUbFwCQZcIeAKDmdejsAQAoEPYAADWvdPW6zh4AILuEPQBAzSsZ0KyzBwDIMGEPAFDzSlav6+wBADKsodIFUH5dXV3R1dU17jODg/6SDEDtSG7j6ukfinw+H7lcrkIVAQBUjrAng3p7e2PJkiWVLgMAJkwy7BkaHo21QyMxpclfdQCA7PE3oAxqb2+PGTNmjPvM4OBgrFixokwVAcCWSV7jiojoGxgS9gAAmeRvQBk0f/78mD9//rjPdHd3x7x588pUEQBsmdYpjSWv9fQPxdadLRWoBgCgsgxoBgBqXn1dLtoSG7l6+82fAwCySdgDAKSC9esAAOsIewCAVOho1dkDABAh7AEAUqI9MaS5d0BnDwCQTcIeACAVSq9x6ewBALJJ2AMApEJy/brOHgAgq4Q9AEAqtJfM7BH2AADZJOwBAFKhvSXR2eMaFwCQUcIeACAVkp09Pa5xAQAZJewBAFLB6nUAgHWEPQBAKpSsXjezBwDIKGEPAJAKydXrfWuGYnQ0X6FqAAAqR9gDAKRCcvV6Ph/Rv0Z3DwCQPcIeACAVkp09ERE9rnIBABkk7AEAUqG5qT4a6ov/atNjSDMAkEHCHgAgFXK5XOlGLuvXAYAMEvYAAKnRbv06AICwBwBIj/aWxPp1nT0AQAYJewCA1Eh29pjZAwBkkbAHAEiN5Pr1Xtu4AIAMEvYAAKlROrNH2AMAZI+wBwBIjeTMHte4AIAsEvYAAKlh9ToAQERDpQug/Lq6uqKrq2vcZwYH/SYUgNrT3mL1OgCAsCeDent7Y8mSJZUuAwAmXHtyQLPOHgAgg4Q9GdTe3h4zZswY95nBwcFYsWJFmSoCgIlRunpd2AMAZI+wJ4Pmz58f8+fPH/eZ7u7umDdvXpkqAoCJkVy9Pjg0EoNDI9HUWF+higAAys+AZgAgNZIzeyJc5QIAskfYAwCkxlhhj/XrAEDWCHsAgNSor6+L1inFt9R7ze0BADJG2AMApIr16wBA1gl7AIBUSa5ft5ELAMgaYQ8AkCodifXrBjQDAFkj7AEAUqW9pbizxzUuACBrhD0AQKq06+wBADJO2AMApEpyQLPV6wBA1gh7AIBU6WhNXuPS2QMAZIuwBwBIleQ2rt4BnT0AQLYIewCAVEnO7LF6HQDIGmEPAJAqJavXhT0AQMYIewCAVEmuXu8bGIzR0XyFqgEAKD9hDwCQKslrXKP5iIG1wxWqBgCg/IQ9AECqJFevR1i/DgBki7AHAEiVluaGqK/LFb1mbg8AkCXCHgAgVXK5XHRYvw4AZJiwBwBInbYW69cBgOwS9gAAqVOyfn1A2AMAZIewBwBInfbkNS4DmgGADBH2AACpk1y/7hoXAJAlwh4AIHVKBjTr7AEAMkTYAwCkTnuLmT0AQHYJewCA1Ele4+p1jQsAyBBhDwCQOu0txde4elzjAgAyRNgDAKSO1esAQJY1VLoAyq+rqyu6urrGfWZw0G9AAahdBjQDAFkm7Mmg3t7eWLJkSaXLAIBJ05YY0LxmcCSGhkejsUFTMwCQfsKeDGpvb48ZM2aM+8zg4GCsWLGiTBUBwMRKdvZERPQODMZWHVMqUA0AQHkJezJo/vz5MX/+/HGf6e7ujnnz5pWpIgCYWMltXBHrNnIJewCALNDLDACkTkN9XbQ01xe9Zv06AJAVwh4AIJXakuvXBwxpBgCyQdgDAKRSyfp1G7kAgIwQ9gAAqVS6ft01LgAgG4Q9AEAqJdev9wh7AICMEPYAAKlU0tljZg8AkBHCHgAgldpbkjN7dPYAANkg7AEAUqm9NXmNS2cPAJANwh4AIJVKr3Hp7AEAskHYAwCkUrKzx+p1ACArhD0AQCp1tOjsAQCySdgDAKRSW8nMnqHI5/MVqgYAoHyEPQBAKiVn9oyO5mNg7XCFqgEAKB9hDwCQSh2Jzp4I69cBgGwQ9gAAqdTS3BB1dbmi16xfBwCyQNgDAKRSLpeL9pbERi5DmgGADBD2AACpVRL2uMYFAGSAsAcASK3kkObeAde4AID0E/YAAKk11vp1AIC0E/YAAKnV0ZLo7DGgGQDIAGEPAJBayfXrBjQDAFnQUOkCAAAmS/Ia19N/7o3fPvLnClVTrLO9KXaZ1VmyHh4AYEsJewCA1EoOaP7D48vi3Mf/t0LVlHrDq2fEFz+4f+RyAh8AYOK4xgUApFbyGle1ufehJfH4M6sqXQYAkDLCHgAgtXbbYatKl7BRL6wcqHQJAEDKuMYFAKTWDjM64ox3z4mf3fV49FTJJq5VvWtjeCRfeL9/jaHRAMDEEvYAAKl21AE7xVEH7FTpMgo+95274vePLSu83zcwXMFqAIA0co0LAKCM2qYUzxHS2QMATDRhDwBAGbVOKW6s7lujswcAmFjCHgCAMtLZAwBMNmEPAEAZtbYUhz19A8IeAGBiCXsAAMqoLXGNq981LgBgggl7AADKqDVxjavPNS4AYIIJewAAysjMHgBgsjVs/BHSpqurK7q6usZ9ZnBwsEzVAEC2tLYktnENuMYFAEwsYU8G9fb2xpIlSypdBgBkks4eAGCyCXsyqL29PWbMmDHuM4ODg7FixYoyVQQA2dGaGNC8ZnAkRkZGo77e7XoAYGIIezJo/vz5MX/+/HGf6e7ujnnz5pWpIgDIjrbE6vWIiP61w9HR2lSBagCANPIrJACAMkpu44qI6BtwlQsAmDjCHgCAMprSVB91dbmi1/rXGNIMAEwcYQ8AQBnlcrlobU5s5DKkGQCYQMIeAIAya03M7el3jQsAmEDCHgCAMmubkuzscY0LAJg4wh4AgDJLDmnud40LAJhAwh4AgDJrS4Q9ZvYAABNJ2AMAUGatLcXXuPoHXOMCACaOsAcAoMx09gAAk0nYAwBQZq2JAc39BjQDABNI2AMAUGYlnT1WrwMAE0jYAwBQZq0trnEBAJNH2AMAUGZtJde4hD0AwMQR9gAAlFlryTUuM3sAgIkj7AEAKLPkzB6dPQDARBL2AACUWWtL8TWuNYMjMTIyWqFqAIC0EfYAAJRZsrMnIqJ/ratcAMDEEPYAAJRZa2JAc4T16wDAxBH2AACUWUtzQ9Tlil/rX6OzBwCYGMIeAIAyy+Vy0ZLcyGVIMwAwQYQ9AAAV0Ja4ytXvGhcAMEGEPQAAFdBa0tnjGhcAMDGEPQAAFdDWUhz29LvGBQBMEGEPAEAFJDdymdkDAEwUYQ8AQAW0Ja5x9Q+4xgUATAxhDwBABejsAQAmi7AHAKACSmf26OwBACaGsAcAoAJKt3Hp7AEAJoawBwCgAtoS17j6B4Q9AMDEEPYAAFRAaWePa1wAwMQQ9gAAVEDpzB6dPQDAxBD2AABUQHIbl7AHAJgowh4AgApoS1zjGlg7EiOj+QpVAwCkSVWEPaOjo7F06dJ47LHHYsWKFUWvAwCkUXJmT0TEgO4eAGACVCzsGR4ejquvvjpOPvnkeP3rXx9vfOMb4+ijj45rrrmm8Mx73/veOOuss+KPf/xjpcoEAJgUbS0NJa8Z0gwATITSv2WUwaOPPhof//jH48knn4yIiHx+XctyLpcrem7RokXx+9//Pv77v/87Pv3pT8fJJ59c5koBACbHlKaGyOUi8uvd3DK3BwCYCGUPexYtWhTve9/7ore3txDyNDY2xtBQ8V9uVq5cGQMDA5HL5WJ4eDi++tWvRmNjY8ybN6/cJQMATLi6uly0NjcUdfP0DQh7AIAtV9ZrXCMjI3HGGWdET09P5PP5OOSQQ+LHP/5x/OY3vyl5dtq0aXH11VfHwQcfHBHrun++/vWvx/PPP1/OkgEAJk1ryfp117gAgC1X1rDnmmuuiSeffDJyuVy8//3vj+9///ux7777RlNT05jPv+pVr4rvf//78f73vz8iItasWRNXXHFFOUsGAJg0yY1cfa5xAQAToKxhz8033xwREdtuu22cffbZm/x5Z599drziFa+IiIi77757UmoDACi31inFN+r7XeMCACZAWcOehx56KHK5XBx++OEb7OYZS1NTUxx++OGRz+cLQ50BAGpdcv26bVwAwEQoa9izYsWKiIiYPXv2Zn/uzJkzIyKir69vQmsCAKiU5DUu27gAgIlQ1rCnra0tIl5eYLNs2bKIiOjo6JjQmgAAKqW1pfgal84eAGAilDXs2X777SOfz8e99967WZ83MjISt99+e+Ryudh+++0nqToAgPIq6ewxswcAmABlDXsOOeSQiIi477774le/+tUmf953v/vdWLx4cUREHHTQQZNSGwBAuSUHNNvGBQBMhIaNPzJxjj/++Ojq6orBwcH45Cc/Gd/+9rdj33333eDza9eujW9/+9vxve99LyIiGhoa4t3vfne5yk2trq6u6OrqGveZwcHBMlUDANnV1pKc2eMaFwCw5coa9my33XbxiU98Ir72ta/FypUr4/3vf3/st99+8ZrXvKbwzKOPPho/+clP4oEHHohbbrklVq5cGfl8PnK5XHzoQx9yjWsC9Pb2xpIlSypdBgBkXuk2Lp09AMCWK2vYExHxwQ9+MJYtWxY/+MEPIp/Px69//ev49a9/HblcLiIirr766rj66qsjIiKfzxc+75hjjokzzzyz3OWmUnt7e8yYMWPcZwYHBwvb0wCAydGWuMZlZg8AMBHKHvZERHzmM5+JfffdN775zW/Go48+Ou6z2223XZx++unxnve8p0zVpd/8+fNj/vz54z7T3d0d8+bNK1NFAJBNpZ09rnEBAFuuImFPRMQRRxwRRxxxRNx3331x7733xqJFi2LVqlUxPDwcnZ2dscMOO8TcuXPjwAMPjIaGipUJADBpkjN7BtYOx8hoPurrchWqCABIg4qnKPvss0/ss88+lS4DAKDsktu4ItYFPu2JEAgAYHOUdfU6AAAvaZtSGuqY2wMAbKlJ6ez53//938n4sgUHHnjgpH59AIByaGluiFwuYr2dFDZyAQBbbFLCnvnz5xe2a020XC4XDz300KR8bQCAcqqry0VLc0P0rzeYud+QZgBgC03azJ7116YDADC21imNRQGPzh4AYEtNStjzhje8YUK+zvLly+Oxxx4rdAkJkACAtGmb0hBL13vfzB4on/41Q3H3/c/GstVrKl0KFdLUUB9vePWM2GFGR6VLgQk1KWHPJZdcssVf4/rrr4/zzz+/KOiZMWNG/MM//MMWf20AgGrRmhjS3OcaF5RFPp+P//O9X8XDTy6vdClU2I//+5H4p4//dew0c2qlS4EJU/HV60nLli2Lc889N2699daIeKmb57jjjovPf/7z0d7eXsnyAAAmVFtizXq/a1xQFkuW9wt6iIiItYMjcdf9zwp7SJWqCnuuv/76+MpXvhKrVq0qhDzbbbdd/MM//EO88Y1vrHB1AAATr3VK8V/H+lzjgrLo7ffvGi9Z0eMqH+lSFWHP0qVL40tf+lIsWLAgIl7q5nnXu94Vn/3sZ3XzAACp1TYl2dnjGheUw9DwaNH79XW52O8121WoGsrt6T/3xOIlvYX3V/cNVrAamHgVD3uuu+66+MpXvhKrV6+OiHVBz8yZM+P888+Pgw8+uMLVAQBMrpLOHte4oCyGRkaK3m9raYzPn7xfhaqh3K77n8fie9f8vvB+T7+wh3SpWNizdOnS+OIXvxi33357RLzUzfOe97wnzj77bN08AEAmlM7s0dkD5ZDs7GlsqKtQJVTC1Namovd7dPaQMhUJe6699tr4x3/8x6JunlmzZsX5558fBx10UCVKAgCoiJJtXGb2QFkIe7Ktoy0R9ujsIWXKGva88MIL8aUvfWnMbp5zzjkn2traylkOAEDFtSWucdnGBeUh7Mm2jkRnz+q+ocjn85HL5SpUEUyssoU911xzTVxwwQWxevXqQsgza9as+MpXvhIHHnhgucoAAKgqrYlrXH2ucUFZlIQ99fUVqoRKmJro7BkeGY01gyPR0lzxsbYwISb9v8l//vOf40tf+lLccccdhZAnl8vFe9/73jj77LOjtbV1sksAAKhapdu4dPZAOQwNFw9o1tmTLcnOnoh1c3uEPaTFpP43+eqrr44LL7ywqJtn++23j/PPPz8OOOCAyfzWAAA1IbmNa2DtcIyO5qOuzlUCmEwlnT2Nwp4saZ3SEPV1uRgZzRdeW90/GK+YrhmBdJiUsGfJkiXxpS99Ke68885CyFNXVxcnnHBCfPrTn46WlpbJ+LYAADUn2dmTz68LfJJbuoCJVXqNS9iTJblcLjpam2Jl79rCazZykSaTEvYcffTR0dvbW3Rt6/Wvf31MnTo1/v3f/32Lv/4nPvGJLf4aAADVIDmzJyKib82QsAcmWemAZjN7sqajrbE47LGRixSZlLCnp6cncrlc0STz7u7u6O7unpCvL+wBANJirPkQ/YY0w6SzjYvk3B6dPaTJpM3sebGrZ6JZhQcApEl9XS5amhtiYO1LAU/fgCHNMNkMaKZk/Xq/s5f0mJSw54wzzpiMLwsAkEptU4rDHhu5YPINjejsybrk+nXXuEgTYQ8AQIW1tjRGrFpTeL/PNS6YdMlrXA3CnsxxjYs0c6IBAFRYciOXzh6YfENDOnuyrqMteY1L2EN6ONEAACqsdUpxs7WZPTD5kp09TbZxZY7OHtJM2AMAUGGlnT2uccFkGxoxoDnrprYVn71m9pAmTjQAgAprbSn+gaPPNS6YdFavo7OHNHOiAQBUWFviGlf/gM4emGzCHpIze/rWDMdIYksb1ConGgBAhbXp7IGyKwl76v1olDVTE509ERE9/c5f0sGJBgBQYa22cUHZDZesXjegOWvaxwx7XOUiHYQ9AAAVVnKNy4BmmHSucdHYUBctzcXn72pze0gJJxoAQIUZ0AzlNzhsGxelc3t09pAWTjQAgAorWb0+IOyByZbs7GkS9mTS1NbE+nWdPaSEEw0AoMJak9e41g7H6Gi+QtVANpRe4zKzJ4tK1q/r7CElhD0AABWW7OzJ5yPWDJrbA5PJzB4iSq9xmdlDWjjRAAAqLDmzJyKib0DYA5NpeMTMHkrXr1u9Tlo0bPwRAAAmU3IbTMSL69dbyl8MZESys6dB2JNJBjSnz5Ll/XHd/zwWHa1Ncdxhu0ZTYzavaAp7AAAqrL4uFy3NDTGw9qVuHhu5YPKMjuZjeKR4LpbOnmxKzuxxjau2DY+Mxmcv+p9YumpNREQ880JvfOp9r69wVZXhRAMAqAJtySHNa1zjgskyPDJa8pqwJ5t09qTLwj+tKAQ9EREPLlpawWoqy4kGAFAFknN7+qxfh0mTvMIVEdFYn82rHllXMrNHZ09NW75e0BMRMaUpu/9eC3sAAKpAciNXv2tcMGkGh0dKXtPZk00dbcVnb0//YOTz+Q08TbVb3lMc9mw1dUqFKqk8JxoAQBVoTVzj6nONCybNWJ09TY1+NMqi5Mye4ZF80fw0akuys2d6h7AHAIAK0tkD5TM81jUunT2ZNDUxsyfC+vValuzsmd4p7AEAoILM7IHyGauzp6Hej0ZZ1NLcEPV1uaLXzO2pXStWJ65x6ewBAKCSbOOC8kmGPQ31dZHL5TbwNGmWy+VKNnKttpGrZi1frbPnRcIeAIAq0Jq4xtXnGhdMmmTY4wpXtiXn9ujsqV3LV68ten/61OYKVVJ5TjUAgCqgswfKZ2ikeBuXsCfbknN7enT21KS1QyMlV6Bd4wIAoKLM7IHy0dnD+jpaE+vXdfbUpOS8noiIrV3jAgCgkmzjgvIZHBL28JLkNS4ze2pTcl5Pc1N9tDQ3bODp9HOqAQBUgdbENa4+17hg0iRXrzc21FeoEqpByTWuPmF7LSoZzjx1SqYHrwt7AACqQFviGtfAmqHI5/MVqgbSzcwe1lfS2dO3dgNPUs3GCnuyzKkGAFAFktu4RvMRA2t198BkMLOH9SVXrxvQXJuWryoOe7bqyO4mrghhDwBAVUhu44qwkQsmi7CH9ZXO7HGNqxat6EmsXc/wcOYIYQ8AQFVoSXT2RET0GdIMk6Ik7Kn3Y1GWlc7s0dlTi0qucWV47XqEsAcAoCrU1+Wipbl4SGz/gM4emAylnT0GNGdZcvX6wNrhkv+OUP2SYc9WZvYAAFANknN7dPbA5HCNi/UlZ/ZERPSa21NzViTCnq2FPQAAVIOSsGdA2AOTYWjYNi5ekpzZExGxWthTU4aGR6InMWtpq6kGNAMAUAWSQ5r7dfbApNDZw/oa6uuiNXH+mttTW5avXlvy2vTOlgpUUj2cagAAVaK1JXmNy8wemAzCHpKS3T3Wr9eW5BWupoa6MbdcZolTDQCgSrQlrnHp7IHJYUAzScm5Pav7nL+1ZNkYw5lzuVyFqqkOwh4AgCqRvEZgZg9MDp09JE3V2VPTkp090zM+nDlC2AMAUDVKO3tc44LJMDRiQDPFSq5xmdlTU5Jr14U9wh4AgKrR2pLo7HGNCyaFzh6SOtqKw3adPbWlJOzpFPY41QAAqoTOHigPYQ9JyWtcq3X21JQViW1cW3Vke+16RES2x1NnVFdXV3R1dY37zOCgww0Ayq01EfaY2QOToyTsqRf2ZF1yQLPOntriGlcpYU8G9fb2xpIlSypdBgCQkFwTaxsXTI6h4eKZPQ22cWWe1eu1TdhTStiTQe3t7TFjxoxxnxkcHIwVK1aUqSIAICKitSXR2eMaF0yKZGdPU6POnqwr6eyxer1mDA2Plly7E/YIezJp/vz5MX/+/HGf6e7ujnnz5pWpIgAgonRmz8Caocjn85HL5SpUEaSTmT0kjbV63flbG1b0rCl5zYBmA5oBAKpGa+Ia12g+YmCt7h6YaGb2kJTs7BkZzRuSXyNWJK5wNTbURXuiUzaLnGoAAFWibYy/nPphAyZeaWePmT1Z19Faev6a21MbkvN6tupo1pEVwh4AgKrR2lx6w77PkGaYcK5xkdTS3BAN9cUBgfXrtWF5Yu26eT3rONUAAKpEfX1dTGkq7jDoH9DZAxNteKR4G5ewh1wuZyNXjUpe49pK2BMRwh4AgKrSOiW5kUtnD0y0ZGdPg7CHGGsjl7CnFiSvcW0t7IkIYQ8AQFVpaym+ytUv7IEJNTqaj+GRfNFrOnuIiJLOntU6e2pCycweYU9ECHsAAKpKaWePa1wwkYZGRkteazKgmYiYWtLZI2yvBcmwZ/rU5gpVUl2EPQAAVaQtEfb0D/hhAyZS8gpXhM4e1jGzpzatKBnQ3FKhSqqLUw0AoIq0Tim+xmVmD0ysoeGRkteEPUSUrl83s6f6DY+Mxqq+4rBnK509ESHsAQCoKm0tic4e17hgQunsYUOS17jM7Kl+K3vWRr54BJfV63/hVAMAqCK2ccHkGhb2sAGucdWe5LyehvpcyX+OWeVUAwCoIm2Ja1z9Azp7YCKN1dnTUO/HIqxer0XJsGdax5Soq8tVqJrq4lQDAKgiOntgciXDnob6usjl/HCIzp5atCIR9mztCleBsAcAoIq0tSQ6e4Q9MKEGEwOaXeHiRcmZPQNrR8bsBKN6LF9tOPOGONkAAKpIaWePa1wwkZI/vDc1+pGIdcaa9aK7p7olr3FtpbOnwMkGAFBF2hJhT/+Azh6YSMmwp9G8Hv4iuXo9wtyeapcMe1zjeomTDQCgirQmBzSvHY58cq8s8LKVhD0N9RWqhGpTX19XMiTf+vXqprNnw4Q9AABVpK2l+DfLo6P5WDM4soGngc2VXL3eYGYP67GRq7YkBzRPF/YUONkAAKpIcmZPhCHNMJGGRgxoZsNs5KodIyOjsaq3eECzsOclTjYAgCqSvMYVEdFnbg9MmNJrXH4k4iXJzp7VOnuq1sretTGauOVsG9dLnGwAAFWkob4umpuKZ4j028gFE0bYw3imlnT2CNur1YrE2vW6ulx0tgl7XlT6qyMAACqqbUpDrF1vTs95P/h1VayHzuVysfuO0+L0d74uOtv9hZraZBsX4zGzp3aUDGfuaI66ulyFqqk+wh4AgCrTOqUxlq/3G8tqmhmxdOVATGlqiE+esG+lS4GXZXA4ObPHNi5eYmZP7UiGPeb1FBNjAwBUmVds1VrpEsa18E8rKl0CvGwlnT1V0DVH9ZjaWjwk38ye6mUT1/icbAAAVebYw15Z1XNEDIymliVXr1fzv2uUX8k1Lp09VWtZ8hqXsKeIa1wAAFVmzu6viK4v/k08/syqyOc3/vxke25ZX/zbVQ8U3u/pH4x8Ph+5nNkI1B4zexiPa1y1IzmgWWdPMWEPAEAV6mxvjn32eEWly4iIiCXL+4veHxnNx8Da4Wid0riBz4DqZRsX4ynt7BkSblep5asHit6fbu16EScbAADj6mgtDXV6rSOmRpWGPQY085Lk6vXR0Xz0rRmuUDWMZ7nOnnEJewAAGFdLc0PJOltXG6hVOnsYT7KzJ8L69Wo0MpqPlb3FYY+ZPcWcbAAAjCuXy5V09wh7qFVDI8nV634k4iVTmuqjITHHyXlXfVb3ro3R0eKhdjp7ijnZAADYqPaW0jkWUIsGh3T2sGG5XC6mtlm/Xu2WJzZx1eXWzbrjJU42AAA2KtnZ0+s33dSo0tXrZvZQzEau6pcMe6Z1NEd9nSHa6xP2AACwUe0lP/zo7KE2DY3o7GF8JRu5dPZUneRwZvN6SjnZAADYKDN7SIuhYTN7GF+ys2e1867qrOgp7uwxr6eUkw0AgI1K/qbb6nVqlW1cbMxUnT1Vb/kqYc/GONkAANgoMyxIC2EPG1N63gm3q01yZo+wp5STDQCAjepoSQxoHvDDD7WpJOyp9yMRxUrCHp09VScZ9pjZU8rJBgDARpUOaPbDD7WptLPHNi6Klaxed95VnRXJzp4Oa9eThD0AAGyU33STFq5xsTGurVa30dF8rOgp3sY1vVNnT5KTDQCAjWov2cY1FPl8vkLVwMtnGxcbY/V6dVvdNxgjo8X/+2NmTyknGwAAG5X8TffwyGisHRzZwNNQvXT2sDHJbVxrBkdKQkIqJ7l2PZeLmNbuGleSkw0AgI3qSHT2RNhQQ+0ZHc2XdAQIe0hKhtsR67pJqA7LEmvXO9ubo96g9RL+iQAAsFGtUxqjLlf8Wu+AH36oLUMjoyWvGdBMUntrU+QS551wu3qUDmd2hWsswh4AADaqri4XbS2GllLbkle4InT2UKq+LhdtUxJzynT2VI3liWtchjOPzckGAMAmSV7l8ptuas1Yc1eEPYwlOaTZ+vXqsTxxjWsra9fH5GQDAGCTJOdY9Prhhxqjs4dNNTW5fl1nT9Wwdn3TONkAANgkyfXrBpZSa4aFPWyikvXrwu2qkezssXZ9bE42AAA2SWlnj2tc1JbBMcKeBlt8GEPy2qpwu3okZ/ZsZUDzmJxsAABskmRnj990U2uSM3saG+oil1y7BKGzp1rl8/mSbVxbu8Y1JmEPAACbpKSzZ0BnD7UlObPHFS42pHRmj/OuGqzuG4zhkXzRazp7xuZ0AwBgkyTDHr/pptYIe9hUOnuqU3I4c0TENNu4xuR0AwBgkyRnWJjZQ60pCXvM62EDkuG2mT3VITmcubO9SWi7Af6pAACwSdp19lDjSjt76itUCdWu5BqX864qLF9tOPOmEvYAALBJkp09PTp7qDHJ1esNOgLYgOQ1rt7+wRgdzW/gacplRWIT13TDmTfI6QYAwCZJXmsYHBqJtUMjG3gaqs/QSOk2LhhL8rwbzUf0rxFwV1ryGtd0nT0b5HQDAGCTJK9xRaz7bTfUisEhA5rZNB1tjSWvrXbeVdzyRGfPVlMNZ94QpxsAAJukraX0hx9XuaglyZk9TWb2sAFTmhqiKREG9hjSXHHJzp6tp+rs2RBhDwAAm6S+LlcS+BhaSi2xep3NUbp+XbhdacsTq9e3EvZskNMNAIBNltxQ4xoXtcTMHjaH9evVJZ/Px4rVBjRvKqcbAACbrN1GLmpYsrPHNi7GM7Wks0fYU0m9A0Ml/w4b0LxhDZUuAACA2pH8TbfOHmpJcvW6zh7GkzzvfnLrH+O/f/2nClVD8t/fCAOaxyPsAQBgk+nsoZaVzOypF/awYcmZPat6B2NVr4C7WnS0NkWjIesb5HQDAGCTJX/T7VoDtcSAZjbHK7ZqqXQJjGO7rVsrXUJVc7oBALDJSjt7hD3UjtKwR1cAG3boPtvHdNeEqlJdLuLYQ3etdBlVzTUuAAA2WenMHte4qB2Dw7ZxseleMb01vv2ZN8VDTyyPtYMjG/8EyiMX8crZnTFr2/ZKV1LVhD0AAGyyDp091LBkZ0+TsIeNaG9tiv1es12ly4DN5nQDAGCTtZfM7NHZQ+2weh3ICqcbAACbbKrV69Sw0tXrZvYA6STsAQBgkyUHNK8ZHImhYbMsqA1DI7ZxAdngdAMAYJMlBzRHGNJM7UgGk8IeIK2cbgAAbLL2lsaS11a7ykWNKF297schIJ2cbgAAbLL6+rponVK80FVnD7VC2ANkhdMNAIDNUrqRS2cPtWEwGfbU+3EISCenGwAAm6UjMaTZRi5qxXByZk+jbVxAOgl7AADYLB0tyc4e17ioDa5xAVnhdAMAYLN0tLnGRW0S9gBZ4XQDAGCztJdc49LZQ20oCXvM7AFSyukGAMBm6TCgmRo0OpqPkdF80Ws6e4C0croBALBZSgc06+yh+g2NjJa81thgQDOQTsIeAAA2S3tiQPNqnT3UgOQVrgidPUB6Od0AANgsVq9Ti4aGRkpeE/YAaeV0AwBgs7SXzOxxjYvqp7MHyBKnGwAAmyXZ2TOwdjiGx5iHAtXEzB4gS4Q9AABslo62ppLXDGmm2o3V2dNQn6tAJQCTT9gDAMBmSQ5ojrB+neo3NFw8s6exoS5yOWEPkE7CHgAANktjQ120NBdff9HZQ7VLdvaY1wOkmRMOAIDNVjKkeUBnD9VN2ANkiRMOAIDN1pG4ytXTJ+yhupWEPfV+FALSywkHAMBma09s5LJ+nWpX2tljExeQXsIeAAA2W0fiGlevAc1UueSA5gbXuIAUc8IBALDZSjt7hD1Ut2RnT1OjH4WA9HLCAQCw2aa2JTt7XOOiupnZA2SJEw4AgM3WnhzQrLOHKmdmD5Alwh4AADZbR/Ia14DOHqqb1etAljjhAADYbO0GNFNjhkaKBzQLe4A0c8IBALDZSjp7+oQ9VLdkZ49tXECaOeEAANhsydXrfWuGY2RkdANPQ+UNu8YFZIgTroo9/vjj8brXvS7e/va3V7oUAIAiydXrERG95vZQxQZt4wIyxAlXpYaHh+Mzn/lMrFmzptKlAACUSHb2RAh7qG7Ja1xNjbZxAekl7KlSF110Ufz+97+vdBkAAGNqaqyP5qbiH5atX6eaDQ0b0AxkhxOuCt13333x7//+73HkkUdWuhQAgA3qaCm+ytXbr7OH6mX1OpAlTrgq09fXF2effXbstNNOcdZZZ1W6HACADUquX9fZQzUrCXvM7AFSrKHSBVDsH//xH+PZZ5+Nyy67LJqbmytdDgDABiXn9gh7qGZWrwNZ4oSrIrfeemv89Kc/jVNPPTVe+9rXVrocAIBxJTdy9fS5xkX1Kl29bkAzkF7CniqxbNmy+Pu///vYa6+94rTTTqt0OQAAG5Xs7OnV2UMVGxoxswfIDidclfjCF74QfX198bWvfS0aGtyuAwCqX0eys8eAZqrY4JBtXEB2OOGqwGWXXRYLFiyIs846K175yldWuhwAgE1SMqB5QGcP1Ss5s6dJ2AOkmBaSKnDDDTdERMQFF1wQF1xwQcnHH3nkkdhjjz1i9uzZcdttt5W7PACAMbnGRS0pXb1uZg+QXsKeKnDsscfGfvvtV/L66tWr4+KLL45tttkmjj/++Ojo6KhAdQAAY3ONi1piZg+QJcKeKnDccceN+frTTz9dCHs+9rGPlbkqAIDx6eyhlgwPF8/ssXodSDMnHAAAL0ty9XrvwFCMjuYrVA2Mr/Qalx+FgPRywm3EF7/4xdhjjz3iG9/4xkafHR0djSuvvDJOPPHE2G+//WKvvfaKww8/PD71qU/FvffeW4ZqAQDKJ9nZk89H9K1xlYvqJOwBssQ1rnHcfPPNccUVV2zSsz09PXH66afHPffcU/T6s88+G88++2zccMMNcfLJJ8dnP/vZTf7+22+/fSxcuHCzagYAKJdkZ09ERE//YEkIBNWgJOypF/YA6SXs2YA77rgjPvnJT27Ss/l8Ps4888xC0HPIIYfECSecENtss008/PDD8b3vfS+eeeaZ6OrqiunTp8eHP/zhySwdAKAsmhvro7GhruiH6F5DmqlCI6P5GElcMdTZA6SZsGcMP/zhD+PrX/96DA1t2l9Wrr/++rjrrrsiYt2w5fXXp8+ZMyfe+ta3xrx582LRokVx0UUXxTHHHBPbbbfdpNS+vhe7il4OHUUAwMbkcrnoaG2M5avXFl7rMaSZKjSUGM4cYfU6kG7CnvU8+eSTceGFF8aCBQsiIqK+vj5GRkr/hyGpq6srIiLa29vjnHPOKfn4tGnT4stf/nLMmzcv1q5dGxdffHGcffbZE1v8GK688sq46KKLJv37AADZ1d7alAh7dPZQfYYTV7giIpoadfYA6eWE+4tLL700jj766ELQs+uuu8aXv/zljX7e4sWL46GHHoqIiMMPPzymTZs25nNz586NXXbZJSIifvGLX0xM0QAAFWb9OrUgOa8nIqLBzB4gxZxwf/Hggw/G0NBQNDU1xUc+8pG46qqrYscdd9zo5/3mN78pvH3AAQeM++x+++0XERHPPPNMPPXUU1tWMABAFehIDGnW2UM1GivsMbMHSDPXuP6iubk53v3ud8dpp50Ws2fP3uTPW7RoUeHtnXfeedxnd9hhh8Lbjz766CaFSVvine98Zxx44IEv63MXLlwY55133gRXBACkjc4easHQyFhhj5k9QHoJe/7i3HPPjbq6zU/3n3/++cLbs2bNGvfZmTNnjvl5k2XWrFkbrQkAYEu0J8Ke1cIeqtDY17hyFagEoDz0Lv7Fywl6IiJWrVpVeLutrW3cZ1tbWwtv9/T0vKzvBwBQTZLXuKxepxolt3E1NtRFLifsAdJL2LOFBgdf+u3VlClTxn12/Y+v/3kAALUq2dlj9TrVKNnZY14PkHZOuS1UX//SXd+N/XYgn88X3n65nUQAANWktLNH2EP1GRoS9gDZ4pTbQutfzVqzZs24z65du7bwdlNT0zhPAgDUho6WZGePa1xUn+SAZsOZgbQT9myh9ef0DAwMjPtsf39/4e3Ozs5JqwkAoFw62kq3cY2O5jfwNFTGWDN7ANLMKbeF1l/T/txzz4377PofnzFjxqTVBABQLu2Ja1yj+YiBtcMVqgbGZmYPkDVOuS202267Fd5+6qmnxn128eLFhbd33XXXSasJAKBcOlpLr6Yb0ky1EfYAWeOU20Jz5swpDGbu7u4e99l77rknIiJmzpwZ22+//aTXBgAw2aY01UdDffGSCuvXqTYlYU+9H4OAdHPKbaGZM2fGnDlzIiLipptuit7e3jGf6+7ujieeeCIiIo466qhylQcAMKlyuVzJ+vXVOnuoMqWdPQY0A+km7JkAJ554YkRErFy5Ms4999wYHS3+H5NVq1bFueeeGxERjY2N8f73v7/sNQIATBbr16l2rnEBWdNQ6QLS4G1ve1tcddVVcdddd8XPfvazeP755+MDH/hAzJgxIxYuXBjf/e5345lnnomIiI997GOxww47VLhiAICJ0279OlVuaMQ2LiBbhD0T5Fvf+laceuqpce+990Z3d/eY83tOPvnk+PCHP1yB6gAAJk9ySLPOHqrN0FBxZ0+DsAdIOWHPBGlvb4+LL744rrnmmrjuuuvikUceiZ6enthqq61in332iXnz5sUBBxxQ6TIBACZcR1vxNS6dPVSb5DWuJmEPkHLCnnHsv//+sXDhwk1+vq6uLo477rg47rjjJrGqLdfV1RVdXV3jPjM46DdyAMCmSXb2WL1OtRkaMaAZyBZhTwb19vbGkiVLKl0GAJAS7SUDmnX2UF0MaAayRtiTQe3t7TFjxoxxnxkcHIwVK1aUqSIAoJbp7KHaDQ0b0Axki7Ang+bPnx/z588f95nu7u6YN29emSoCAGpZR8k2LmEP1UVnD5A1TjkAALaIa1xUu5Kwp96PQUC6OeUAANgiY13jyufzFaoGSiXDHqvXgbRzygEAsEWSnT0jo/kYWDtcoWqgVOnMHtu4gHQT9gAAsEWmtjWVvOYqF9Uk2dnT1OjHICDdnHIAAGyRluaGqKvLFb1mSDPVxMweIGuccgAAbJFcLhcdhjRTxWzjArLGKQcAwBZrT65fH9DZQ/UoDXvM7AHSTdgDAMAWS3b29PQJe6geQyM6e4BsccoBALDF2kvWr7vGRfUYTmzjsnodSDunHAAAW6yks8eAZqqImT1A1jjlAADYYh2Jzh4Dmqkmg8IeIGOccgAAbLGOtuQ1Lp09VI9kZ0+TAc1AyjVUugAAAGpfR0vxNa4XVg7Eo4tXVKiaYm0tjTFz67bI5XKVLoUKGBnNx+hovug1nT1A2gl7AADYYskBzY8/syrO+uadFaqm1L57viK++MH9o6HeD/lZM5QYzhwh7AHSzykHAMAWS87sqTa/feTP8cCipZUugwoYTlzhihD2AOnnlAMAYIvtMmtqVPstqede6K10CVRAcl5PROjwAlLPNa4M6urqiq6urnGfGRw0VBEA2HRbTZ0SZ7x7Tlx+yx9jZc/aSpcTERHDwyOx/qiWVX3+fpNFY4U9OnuAtBP2ZFBvb28sWbKk0mUAACnzN/vvFH+z/06VLqPgop/8Lm761Z8K76/srY4QivIaGhkr7LGNC0g3YU8Gtbe3x4wZM8Z9ZnBwMFasqI4NGgAAL0dne3PR+6t7dfZk0djXuKr8ziHAFhL2ZND8+fNj/vz54z7T3d0d8+bNK1NFAAATr7O9eGi0zp5sGhwq3sbV2FAXuWofMAWwhVxWBQAglTrbEp09fcKeLEp29jSZ1wNkgJMOAIBUmpa4xrWyxzWuLEquXjevB8gCYQ8AAKk0NXGNq3dgMEbGGNZLuiUHNDfo7AEywEkHAEAqJTt78vmI1f26e7JmaLh0Zg9A2jnpAABIpY62ppLXbOTKnuTMHmEPkAVOOgAAUqmhvi7aWxqLXltlSHPmCHuALHLSAQCQWp2Jq1yrDGnOnJKwp96PQED6OekAAEitzsSQZp092TNYMrPHNi4g/YQ9AACkVrKzZ2WvsCdrSlavN/oRCEg/Jx0AAKmVDHsMaM4e17iALHLSAQCQWslrXDp7sseAZiCLnHQAAKRWZ1uis6dPZ0/WCHuALHLSAQCQWtOSM3t6dPZkzdBIMuwxoBlIP2EPAACpNTVxjWu1bVyZo7MHyCInHQAAqZXs7OnpH4rhRKcH6TZUsnrdj0BA+jnpAABIrWRnT0REj7k9mTI0pLMHyB4nHQAAqTW1tSlyueLXbOTKFjN7gCwS9gAAkFr19XXR3pKY29OrsydLzOwBsshJBwBAqnUmrnKtMqQ5U8zsAbKoodIFUH5dXV3R1dU17jODg37jBQCkQ2d7czz9597C+65xZYvOHiCLhD0Z1NvbG0uWLKl0GQAAZZHs7HGNK1tKwp56YQ+QfsKeDGpvb48ZM2aM+8zg4GCsWLGiTBUBAEyezsT6dZ092aKzB8giYU8GzZ8/P+bPnz/uM93d3TFv3rwyVQQAMHk624rDntVWr2dKadhjGxeQfmJtAABSbVriGtfKHp09WVK6et2PQED6OekAAEi1qe3Jzh5hT5YMDRVv42oQ9gAZ4KQDACDVppXM7HGNK0uS17iahD1ABjjpAABItamJa1x9A0MlAQDp5RoXkEVOOgAAUi3Z2RPhKleWGNAMZJGwBwCAVGtvbYpcrvg1G7myYWQ0H6Oj+aLXdPYAWeCkAwAg1errcjG1zUauLBoaHil5TdgDZIGTDgCA1JvaVnyVa5XOnkwYHmM2k7AHyAInHQAAqdeZGNK8uldnTxaMNYi7od6PQED6OekAAEi9zpL168KeLBjU2QNklJMOAIDU60zM7DGgORvGmtnT1GgbF5B+wh4AAFIvuX7dgOZsSF7jyuXWDewGSDthDwAAqTc1Efbo7MmGZNjTWF8XuZywB0g/YQ8AAKlX0tljZk8mlIQ95vUAGeG0AwAg9abaxpVJydXrjQ3m9QDZIOwBACD1kp09fWuGxxzeS7oMjRSHPQ06e4CMcNoBAJB6UxPbuCLM7cmCZKDnGheQFU47AABSr6O1KZJLmGzkSr/BITN7gGxy2gEAkHp1dbmY2lZ8lWuVzp7USw5obhL2ABnhtAMAIBM6E0OaVxnSnHrJmT0GNANZIewBACATOhNDmlf16uxJOzN7gKxy2gEAkAnJIc2r+3T2pF1y9bptXEBWNFS6AMqvq6srurq6xn1mcNBvugCAdEmuXzegOf2SM3sa64U9QDYIezKot7c3lixZUukyAADKamoi7LF6Pf1Kwh6dPUBGCHsyqL29PWbMmDHuM4ODg7FixYoyVQQAMPmmJQY0rzSgOfWEPUBWCXsyaP78+TF//vxxn+nu7o558+aVqSIAgMlX0tljQHPq2cYFZJVoGwCATCiZ2aOzJ/UGh2zjArLJaQcAQCYkt3ENrB0uCQNIF9e4gKxy2gEAkAnTOppLXlvlKleqJVevC3uArHDaAQCQCW1TGqOuLlf02qo+V7nSzMweIKuEPQAAZEJdXS46E1e5Vpnbk2qucQFZ5bQDACAzOhNDml3jSrehYQOagWxy2gEAkBmd7Tp7skRnD5BVTjsAADKjsy3Z2SPsSbOSsKfejz9ANjjtAADIjKmJzp7Vfa5xpZnOHiCrnHYAAGTGtMTMnpU6e1KtZGZPo21cQDYIewAAyIypibBntQHNqaazB8gqpx0AAJkxLXGNS2dPupnZA2SV0w4AgMyYmhjQvLpP2JNmOnuArHLaAQCQGdM6isOegbUjsXZoZANPU+uGRoQ9QDY57QAAyIzOtqaS16xfT6/Szh4DmoFsEPYAAJAZbS2NUV+XK3pN2JNOI6P5GB3NF72mswfICqcdAACZkcvlojMxpHmVjVyplFy7HiHsAbLDaQcAQKZ0Jtav6+xJp+QVrghhD5AdTjsAADKlsy0Z9ujsSaOxwx4ze4BsEPYAAJApUxPXuKxfTyedPUCWOe0AAMiUaYlrXCtd40olM3uALHPaAQCQKcnOHte40inZ2ZPLRckmNoC0EvYAAJApyc4eA5rTKRn2NNbXRS4n7AGyQdgDAECmTE0OaO7T2ZNGJWGPK1xAhjjxAADIFJ092TBcEvbYxAVkh7AHAIBM6UzM7Fk7OBJr1g5XqBomy2BiQHODzh4gQ5x4AABkSmeisyfCVa40Sl7jahL2ABnSUOkCKL+urq7o6uoa95nBQX/hAQDSqXVKQzTU52J4JF94bVXv2pgxvbWCVTHRzOwBskzYk0G9vb2xZMmSSpcBAFARuVwuOtubY9mqNYXXzO1JH2EPkGXCngxqb2+PGTNmjPvM4OBgrFixokwVAQCUV2dbMuzR1Zw2QyMGNAPZJezJoPnz58f8+fPHfaa7uzvmzZtXpooAAMorOaRZZ0/6DCUGNOvsAbLEiQcAQOYkhzQb0Jw+ydXrtnEBWeLEAwAgc6bq7Em9kpk99X70AbLDiQcAQOZMS3b2CHtSZ9CAZiDDnHgAAGTO1DbXuNIu2dnT1GhAM5Adwh4AADJnmmtcqWdAM5BlTjwAADKnZEBz72Dk8/kKVcNkMLMHyDInHgAAmZMMewaHRmLN4MgGnqYWJcMe27iALHHiAQCQOZ2Ja1wRrnKlTXL1umtcQJY48QAAyJyW5oaSH/6FPekyNJIMewxoBrJD2AMAQObkcrnobEsMabaRK1VKZvbo7AEyxIkHAEAmdXYkhjT36OxJE9u4gCxz4gEAkEmdbYmwR2dPqgwO6ewBssuJBwBAJiWHNJvZky7JmT1Nwh4gQ5x4AABkUnL9urAnXZLbuBoMaAYyRNgDAEAmTTWgOdXM7AGyzIkHAEAmTdPZk2q2cQFZ5sQDACCTSq9x6exJk5Kwp96PPkB2OPEAAMiksQY05/P5ClXDRNPZA2SZEw8AgExKdvYMDY/GwNrhClXDREtu4xL2AFnixAMAIJOSYU+Eq1xpMjSUHNBsGxeQHcIeAAAyaUpTfTQluj1W9RnSnAYjI6MxmriRp7MHyJKGShcAAACVkMvlorOjOV5YMVB47erbF8WM6W0VrGrz7DJrahy418yY0uyv9etLzuuJEPYA2eJ/FQAAyKzOtqaisOeXDzxXwWpentnbtsVn3j83Xrn9tEqXUjWS83oihD1AtjjxAADIrGkdUypdwhZ75oW++PQ//09cd+djton9xdidPWb2ANkh7AEAILMOfO3MSpcwIYZHRuN71/4+zv+Pe2JVr7lDrnEBWecaFwAAmfXm/XaMpsb6eODRF2IkOdG3yj3x7Kp44tnVRa/d89Dz8fF/uj0+Pe/18dpdt6lQZZU3NDxS8pqwB8gSYQ8AAJmVy+XisH23j8P23b7SpWy2oeHR+NHPH46rbl9U9Pry1WviC/92d7znyN3jhDfvEfX12Qs5kp09uVxEfV2uQtUAlF/2Tn4AAEiBxoa6mP93r4kv/38HxrT25qKP5fMRl9/8x/j8v95dNIA6K5JhT2N9XeRywh4gO3T2AABADdt3z1fEP3/qsPh/P/5t/O7RF4o+9tATy+O0r90aU9uaKlRdZZSEPY2GMwPZIuwBAIAat9XUKfHlDx8YVy54NH70i0didL35Q2sHR+KFwex196zPvB4ga5x6AACQAnV1uXj3EbvHV884JF4xvbXS5VSV6R1TKl0CQFkJewAAIEX23Gl6fOusw+KNc2ZXupSqkMtFHHf4rpUuA6CsXOMCAICUaW9pjLNPnBsnH/3qeOr5nkqXU1E7z5wa20xrqXQZAGUl7AEAgJR6xVat8YqtXOkCyBrXuAAAAABSRGdPBnV1dUVXV9e4zwwODpapGgAAAGAiCXsyqLe3N5YsWVLpMgAAAIBJIOzJoPb29pgxY8a4zwwODsaKFSvKVBEAAAAwUXL5fD5f6SKoPt3d3TFv3rzC+5deemnMnTu3ghUBAAAAm8KAZgAAAIAUEfYAAAAApIiwBwAAACBFhD0AAAAAKSLsAQAAAEgRYQ8AAABAigh7AAAAAFJE2AMAAACQIsIeAAAAgBQR9gAAAACkiLAHAAAAIEWEPQAAAAApIuwBAAAASBFhDwAAAECKCHsAAAAAUkTYAwAAAJAiwh4AAACAFBH2AAAAAKSIsAcAAAAgRYQ9AAAAACki7AEAAABIEWEPAAAAQIoIewAAAABSpKHSBVCd+vv7i95fuHBhhSoBAACA7Nhjjz2io6Nji76GsIcxLV68uOj98847r0KVAAAAQHZceumlMXfu3C36Gq5xAQAAAKSIsAcAAAAgRVzjYkxvetObit7fcccdo6WlZZM//4wzzogVK1bEVlttFRdddNFmfe/N/dyFCxcWXTP70pe+FHvsscdmfU82bkv+M6021fhnqURNk/k9J/prT8TXcy6lTzX+u/xyVeOfxbk0+V/PuZQ+1fjv8stVjX8W59Lkfz3nUnWYiH8Owh7GNHPmzJg3b97L/vympqbC/9/cu4Zb8rkR6/7F2NL7jZTa0v9cqkk1/lkqUdNkfs+J/toT8fWcS+lTjf8uv1zV+GdxLk3+13MupU81/rv8clXjn8W5NPlfz7mUHq5xAQAAAKSIsAcAAAAgRYQ9AAAAACki7AEAAABIEWEPAAAAQIoIewAAAABSRNgDAAAAkCINlS6AdJo/f3709vZGe3t7WT+XyZOm/1yq8c9SiZom83tO9NeeiK/nXEqfNP3nUo1/FufS5H8951L6pOk/l2r8sziXJv/rOZfSI5fP5/OVLgK2RHd3d8ybN6/w/qWXXhpz586tYEVA1jmXgGrjXAKqjXNpcrnGBQAAAJAiwh4AAACAFBH2AAAAAKSIsAcAAAAgRWzjoubNmjUrzjjjjKL3ASrJuQRUG+cSUG2cS5PLNi4AAACAFHGNCwAAACBFhD0AAAAAKSLsAQAAAEgRYQ8AAABAigh7AAAAAFLE6nUybc2aNXHxxRfH9ddfH4sXL47W1tbYf//94yMf+UjsueeelS4PIB5//PE49thjY+edd45rr7220uUAGXTJJZfE+eefv8GPX3rppTF37twyVgSwzr333hvf+9734ne/+10MDQ3FDjvsEMcdd1y8733vi6ampkqXV1HCHjJrcHAwTjnllLj33nvjNa95TZxwwgmxfPny+PnPfx7//d//Hd/5znfi0EMPrXSZQIYNDw/HZz7zmVizZk2lSwEy7KGHHoqIiJNOOik6OjpKPj5r1qxylwQQV1xxRXzpS1+K9vb2eOtb3xpTpkyJ22+/PS644IL43e9+F9/4xjcil8tVusyKyeXz+Xyli4BK+MEPfhBf+9rX4phjjomvfe1rhYPgwQcfjBNOOCG23nrruPXWW6OhQSYKVMY3v/nN+Nd//deIiNhzzz119gAV8Y53vCMef/zxuO+++6K+vr7S5QDE448/Hm9/+9tj2223jYsvvji23377iIhYu3ZtfOADH4jf/e530dXVFQcddFCFK60cM3vIrCeffDKmTZsWH/vYx4oS39e+9rWx6667xvPPPx/PPPNMBSsEsuy+++6Lf//3f48jjzyy0qUAGTY4OBiLFi2K3XffXdADVI2LL744BgcH4+///u8LQU9ERHNzc3zyk5+Md77znTE8PFzBCitPywKZ9Q//8A/xD//wDyWvDwwMxDPPPBONjY0xffr0ClQGZF1fX1+cffbZsdNOO8VZZ50Vt9xyS6VLAjLq0UcfjaGhoXjVq15V6VIACm677bbo6OgYc+zGAQccEAcccEAFqqouwh74i/7+/vjDH/4Q3/jGN2L16tXxkY98ZMx76QCT7R//8R/j2Wefjcsuuyyam5srXQ6QYS/O68nlcnHWWWdFd3d3rFy5Mnbeeed497vfHfPmzYu6OpcFgPJZsWJFLFmyJObMmROrVq2Kb3/723HLLbfE8uXLY4cddoj3vOc98YEPfCDzZ5OwByKiu7s75s2bV3j/fe97X5x11lkVrAjIqltvvTV++tOfxhlnnBGvfe1r4+mnn650SUCGPfzwwxERcfnll8d+++0XRx99dCxdujTuuOOOOP/88+Pee++Nb37zm5n/oQoonyVLlkTEumum73rXu2J4eDgOO+ywyOfzsWDBgrjgggviwQcfjH/6p3+qcKWVJeyBiGhoaIgPfvCD0dfXF7fffnv8+Mc/jpUrV8ZXv/rVzK/sA8pn2bJl8fd///ex1157xWmnnVbpcgAil8vFrFmz4uMf/3gce+yxhdeXLl0aJ598ctx0001xxRVXxPHHH1/BKoEs6evri4h1nYd77713/OAHP4ipU6dGRMTy5cvjxBNPjJ/97GdxxBFHxN/+7d9WstSKEsFDRMyZMyfOOeecOO+88+LGG2+MffbZJ2688cb48Y9/XOnSgAz5whe+EH19ffG1r33NJkCgKnzxi1+MBQsWFAU9ERHbbLNNfPazn42IiGuuuaYClQFZtf6w+L//+78vBD0REdOnT49PfvKTERFx3XXXlb22aiLsgYT29vb4zGc+ExERN998c4WrAbLisssuiwULFsRZZ50Vr3zlKytdDsBGve51r4uIiMWLF1e4EiBLXpyrWl9fH69+9atLPr7XXntFRMSf/vSnstZVbfzakEwaGRmJe+65J1avXh1HHXVUycd33HHHiFjXBghQDjfccENERFxwwQVxwQUXlHz8kUceiT322CNmz54dt912W7nLAzJoaGgoHnnkkVizZk284Q1vKPl4f39/RIRB8kBZ7bDDDtHY2BhDQ0MxMjISjY2NRR8fGhqKiIiWlpZKlFc1hD1kUl1dXXz84x+Pnp6euOOOO2LGjBlFH//DH/4QERE77bRTJcoDMujYY4+N/fbbr+T11atXx8UXXxzbbLNNHH/88bYEAmUzOjoaJ5xwQoyMjMTdd98d06dPL/r4PffcExHrrsMDlEtTU1Pss88+cc8998Rdd90VRx55ZNHH77///oiIeNWrXlWJ8qqGsIdMyuVyccwxx8SPfvSjuPDCC+Of/umfClsklixZUvitumGDQLkcd9xxY77+9NNPF8Kej33sY2WuCsiy5ubmePOb3xw33nhjXHjhhXHhhRcW/r701FNPxde//vWoq6uLk046qcKVAlnzgQ98IO655574v//3/8acOXNim222iYiIF154Ib71rW9FLpeL97znPRWusrKEPWTWmWeeGd3d3XHjjTfGY489FgcddFCsXLkybr311li9enWceuqpcdhhh1W6TACAivnsZz8bDzzwQFx77bWxcOHCOPDAA2Pp0qVx6623Rn9/f3zuc58rzO4BKJc3v/nNceKJJ8Yll1wSb3vb2+Itb3lLRETccsstsXTp0jj99NMzfzbl8vl8vtJFwPq++MUvxhVXXBGnnnpqYZL6hoyOjsbVV18d11xzTSxcuDD6+/tj2223jX333TeOP/74Me+Xr6+/vz++973vxY033hjPPPNMTJkyJV772tfGSSedJOgBCsp5LiU9/fTTccQRR8See+4Z11577Zb8MYAUKee5tHLlyvi3f/u3uOWWW+L555+P1tbW2HvvveNDH/pQHHjggRP5xwJqXLn/znTjjTfGpZdeGg899FDkcrnYY4894qSTTiqEP1km7KGq3HzzzXHGGWdERGz0gOjp6YnTTz+9cF88KZfLxcknn1xYCwrwcjiXgGrjXAKqkbOpurjGRdW44447Npr+viifz8eZZ55ZOBwOOeSQOOGEE2KbbbaJhx9+OL73ve/FM888E11dXTF9+vT48Ic/PJmlAynlXAKqjXMJqEbOpupTV+kCICLihz/8YXz0ox8trMnbmOuvvz7uuuuuiFg31PQHP/hBHHnkkTFnzpw44YQT4qqrropdd901IiIuuuiieP755yetdiCdnEtAtXEuAdXI2VSdhD1U1JNPPhmnnnpqXHDBBTE0NBT19fWb9HldXV0REdHe3h7nnHNOycenTZsWX/7ylyMiYu3atXHxxRdPXNFAqjmXgGrjXAKqkbOpugl7qJhLL700jj766FiwYEFEROy6666Ff6nHs3jx4njooYciIuLwww+PadOmjfnc3LlzY5dddomIiF/84hcTUzSQas4loNo4l4Bq5GyqfsIeKubBBx+MoaGhaGpqio985CNx1VVXxY477rjRz/vNb35TePuAAw4Y99n99tsvIiKeeeaZeOqpp7asYCD1nEtAtXEuAdXI2VT9DGimYpqbm+Pd7353nHbaaTF79uxN/rxFixYV3t55553HfXaHHXYovP3oo49u0gEEZJdzCag2ziWgGjmbqp+wh4o599xzo65u85vL1h/QNWvWrHGfnTlz5pifBzAW5xJQbZxLQDVyNlU/17iomJdzOERErFq1qvB2W1vbuM+2trYW3u7p6XlZ3w/IDucSUG2cS0A1cjZVP2EPNWdwcLDw9pQpU8Z9dv2Pr/95ABPJuQRUG+cSUI2cTeUj7KHmrL/SL5fLjftsPp8vvP1y02eAjXEuAdXGuQRUI2dT+fgnRs1Zv51vzZo14z67du3awttNTU2TVhOQbc4loNo4l4Bq5GwqH2EPNWf9u50DAwPjPtvf3194u7Ozc9JqArLNuQRUG+cSUI2cTeUj7KHmrL/a77nnnhv32fU/PmPGjEmrCcg25xJQbZxLQDVyNpWPsIeas9tuuxXefuqpp8Z9dvHixYW3d91110mrCcg25xJQbZxLQDVyNpWPsIeaM2fOnMIwr+7u7nGfveeeeyIiYubMmbH99ttPem1ANjmXgGrjXAKqkbOpfIQ91JyZM2fGnDlzIiLipptuit7e3jGf6+7ujieeeCIiIo466qhylQdkkHMJqDbOJaAaOZvKR9hDTTrxxBMjImLlypVx7rnnxujoaNHHV61aFeeee25ERDQ2Nsb73//+stcIZItzCag2ziWgGjmbyqOh0gXAy/G2t70trrrqqrjrrrviZz/7WTz//PPxgQ98IGbMmBELFy6M7373u/HMM89ERMTHPvax2GGHHSpcMZB2ziWg2jiXgGrkbCoPYQ8161vf+laceuqpce+990Z3d/eYdz5PPvnk+PCHP1yB6oAsci4B1ca5BFQjZ9PkE/ZQs9rb2+Piiy+Oa665Jq677rp45JFHoqenJ7baaqvYZ599Yt68eXHAAQdUukwgQ5xLQLVxLgHVyNk0+XL5fD5f6SIAAAAAmBgGNAMAAACkiLAHAAAAIEWEPQAAAAApIuwBAAAASBFhDwAAAECKCHsAAAAAUkTYAwAAAJAiwh4AAACAFBH2AAAAAKSIsAcAAAAgRYQ9AAAAACki7AEAAABIEWEPAAAAQIoIewAAAABSRNgDAAAAkCLCHgAAAIAUEfYAAAAApIiwBwAAACBFhD0AAAAAKSLsAQAAAEgRYQ8AAABAigh7AAAAAFJE2AMAAACQIg2VLgAAoNbccccd8eEPfzgiImbPnh233XbbRj/nAx/4QPz617+OiIhrrrkmXvWqVxV9/O67747rr78+uru7Y+nSpZHL5eIVr3hF7L///vGud70r9t57702q7cEHH4wbb7wxuru747nnnouVK1dGY2NjdHZ2xp577hmHH354HHvssdHU1DTm5++xxx4REfHe9743zjvvvPiv//qv6Orqiueffz622WabmDt3bnziE5+I2bNnR0TE4OBgXHfddXHLLbfE73//+1i5cmU0NzfHNttsE/vuu2/8zd/8TRx++OGbVDsAMDFy+Xw+X+kiAABqycjISBx66KHxwgsvRETEZZddFvvss88Gn1+yZEkcdthhMTo6Grvttlv87Gc/K3xs5cqV8alPfSruuuuucb/ncccdF1/+8pc3GNL09vbGOeecE7fccstG699hhx3ie9/7Xuyyyy4lH1s/7Nlll13iwgsvLPp4Y2Nj3H333dHZ2RmLFy+OU045JZ588slxv9++++4b3/nOd2KrrbbaaG0AwJbT2QMAsJnq6+vj6KOPjq6uroiIuOGGG8YNe2644YYYHR2NiIi3v/3thddXrVoVxx9/fDzxxBMREdHS0hJvetOb4pWvfGWMjIzEww8/HP/zP/8TQ0NDcdVVV8Xzzz8f3//+96O+vr7o64+OjsYpp5wS9913X+HrHHroofHKV74ypkyZEitWrIh77rknfv/730dExOLFi+MTn/hEXHPNNVFXN/at/j/96U9x1VVXlbx+wAEHRGdnZwwODsapp55aCHpmzpwZhx12WMycOTP6+/vjj3/8Y9x+++0xOjoav/3tb+PjH/94XHLJJZvyjxcA2ELCHgCAl+Ed73hHIez5+c9/Hp/73OdKQpgXXX/99RERUVdXF3/3d39XeP0LX/hCIeg5+OCD4//+3/8bW2+9ddHnPvHEE3HGGWfEokWL4pe//GX867/+a5xxxhlFz1x99dWFoGf77bePH/3oRzFz5sySOm644Yb49Kc/HaOjo7Fw4cK477774vWvf/2YNf/qV7+KiHUdRR/96Edj2rRp8Zvf/KbQWXTTTTfFokWLIiJiv/32i+9///vR3Nxc9DUeeOCBOOmkk6K/vz/uueee+O1vfxv77rvvmN8PAJg4BjQDALwMe+65Z+HK09KlSwvhSNJjjz0WDz30UERE7L///rHddttFxLog5Oabb46IiN133z3+9V//tSToiYjYZZdd4t/+7d8KQcp//Md/RE9PT9EzV199deHtL3zhC2MGPRERb3vb2+KNb3xj4f0XO3025JBDDokLLrggtt9++2hvb49DDz00DjzwwIiIuP/++wvPnXTSSSVBT0TE3nvvHSeffHJErOuGeuCBB8b9fgDAxBD2AAC8TO94xzsKb68/h2d91113XeHt9a9w/eQnPym8/aEPfWjMsORFO+ywQxxzzDEREdHX1xe33npr0cePP/74+NjHPhbvfe9749BDDx235hcDqogoCY2S3ve+923wY+t3Mb3YVTSWk046KW666aa4//77C8EPADC5hD0AAC/T3/3d3xVCj5tvvjkGBwdLnnkxBGppaYk3v/nNhdfvueeewtuvec1rNvq91r/+9Jvf/KboY0cffXScccYZcd55523wKlnEuhlBS5YsKbw/PDw87vccbw7RG97whsLb3//+9+PMM8+MBQsWRH9/f9Fz06ZNi5133jkaGxvH/V4AwMQxswcA4GXadttt4+CDD44777wzenp64s4774wjjzyy8PHf/va38fTTT0dExBFHHBHt7e0RsS5k+dOf/lR47uijj96s7/vss8+O+/FVq1bFE088EU899VQ89dRT8cQTT8QjjzwSjz32WKy/iHW8paytra0xffr0DX78TW96U+y3336F0OrnP/95/PznP4/GxsbYd999441vfGMceuihsfvuu2/Wnw0A2HLCHgCALfCOd7wj7rzzzohY18Wzftjz4mDmiOIrXKtXrx43aNmYVatWlbw2ODgYl112WfzkJz+JP/7xjxv83Pr6+hgZGdno9+jo6Bj343V1dfGd73wnzj///Lj22msLf56hoaH49a9/Hb/+9a/j61//euy8887xzne+M0488cRoaWnZ6PcFALacsAcAYAsceeSR0dHRET09PbFgwYLo6+uLtra2GB4ejp///OcREbHNNtvEwQcfXPic9a9PTZkyJT7+8Y9v1vdMDnJ+4YUX4pRTTolHHnmk6PW6urqYPXt27LbbbvHa17429t9//7jrrrviO9/5zka/R0PDxv+a2NHREV/96lfjox/9aNxwww1x2223xe9///vCmvmIiCeffDL+6Z/+KS6//PK45JJLYtasWZv4pwQAXi5hDwDAFmhubo6jjjoqfvrTn8aaNWtiwYIFcfTRR8cvf/nLWLFiRUSsu6a1/iydzs7Owttr166NE088sbDS/OX49Kc/XQh6ttlmm5g/f34ccMABsdtuu5UMfk4Od54IO+64Y5x22mlx2mmnxapVq+Lee++Nu+++O2677bZ4/vnnIyLi6aefjs997nPxn//5nxP+/QGAYgY0AwBsofW3ct12221F/z+i+ApXxLqAaNttt42IdXNzNmUl+cDAQAwMDJS8fv/99xfWvre2tsbll18ep5xySuy1115jbvh6MYCaLJ2dnXHkkUfGueeeG7fffnt8/vOfL3zsV7/6VdGAaABgcgh7AAC20Ny5c2P77bePiIg777wzhoaG4o477oiIiN122y1e/epXl3zO+tusNrS2fX3/7//9v5gzZ04cfPDB8e1vf7vw+u9+97vC2wcffHChjrHk8/miLWDrX7faHCMjI3HOOefEcccdFwcffPCYW8giInK5XJx00kmx8847F14T9gDA5BP2AABsoVwuV+je6enpiUsvvbSwMeuYY44Z83PW7wa68sorS+btrG/x4sVxxRVXRETE0qVLY6+99ip8bP2gZWNdO5dccklhO1jExlevb0h9fX088sgj8Yc//CGWLl067tWw0dHRWLlyZeH9GTNmvKzvCQBsOmEPAMAEWD+8+Zd/+ZeIWDcgeUNhz6GHHhr77LNPRKwLbE455ZT47W9/W/Lck08+GaeeemqsWbMmIiL23nvvOPTQQwsf33PPPQtv/+Y3v4mbb7655GsMDg7Gv/3bv8WFF15Y9PpY18I21Tvf+c7C2//n//yfePDBB0ueyefzceGFFxbCnte97nXCHgAog1x+S/Z+AgBQcPzxx8d9991XeP/AAw+MH/7whxt8/rnnnot3v/vd8cILL0TEug6hAw44IObMmRO5XC4WLVoUt912W6EDp7OzM6644oqia1EjIyPx9re/PR599NHCa4cccki8+tWvjubm5njmmWdiwYIFha6fxsbGGBoaioiIt771rfHNb36zqKY99tgjIiJmz55dNHcoaXBwMN71rnfFwoULC7UfdNBBsfvuu8f06dNj2bJlcdddd8WiRYsiIqKpqSn+8z//M/bdd9/x/hECABNA2AMAMEEuu+yyOPfccwvvX3jhhXHssceO+znPP/98nHnmmUUh0Vhe+cpXxje/+c3YfffdSz72xBNPxAc/+MHC1bENmTt3bnzyk5+MefPmRcS6LVrJTqBNDXsiIv785z/HaaedFr///e/HfW7bbbeNr3zlK0UdSQDA5BH2AABMkNWrV8dBBx0UQ0ND0dLSEnfddVe0t7dv0ucuWLAgfv7zn8d9990XS5cujaGhodhqq63i1a9+dRx11FFx9NFHj7ueffXq1fGjH/0obrvttnjiiSdiYGAgWlpaYubMmfGqV70q/vZv/zYOO+ywyOVy8da3vjUef/zxiIi44oor4nWve13h62xO2BOxrrPoF7/4Rdx0003xhz/8IZYtWxbDw8Mxffr0eOUrXxmHH354HHfccZv8zwEA2HLCHgCACbJ48eI48sgjIyLi7/7u7+LrX/96hSsCALLIgGYAgAly3XXXFd4+7rjjKlgJAJBlwh4AgAkwODgYP/3pTyMiYvvtt48DDzywwhUBAFkl7AEAeBkGBwcLW61WrlwZ55xzTmFA8vvf//7I5XKVLA8AyLCGShcAAFCLnnjiiXjnO98ZU6dOjZUrV8bIyEhErOvqOf744ytcHQCQZcIeAICXYbvttouhoaFYtmxZ4bXW1tb46le/Gi0tLRWsDADIOte4AABehs7Ozth///2jtbU1Ojs747DDDotLL7005s6dW+nSAICMs3odAAAAIEV09gAAAACkiLAHAAAAIEWEPQAAAAApIuwBAAAASBFhDwAAAECKCHsAAAAAUkTYAwAAAJAiwh4AAACAFBH2AAAAAKSIsAcAAAAgRYQ9AAAAACki7AEAAABIEWEPAAAAQIoIewAAAABSRNgDAAAAkCLCHgAAAIAUEfYAAAAApIiwBwAAACBFhD0AAAAAKfL/A+hYPkpuk6TkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 444,
       "width": 573
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with sns.axes_style('ticks') :\n",
    "    fig, ax = plt.subplots()\n",
    "    sns.lineplot(data=sampledemog_data, x='years', y='Ne', hue='population', ax=ax)\n",
    "    ax.set_xscale('log')\n",
    "    ax.set_yscale('log')\n",
    "    handles, labels = ax.get_legend_handles_labels()\n",
    "    ax.legend(handles=handles[1:], labels=labels[1:], frameon=False)\n",
    "    sns.despine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1d017b69-9da1-4b4d-a08c-4b5c7aca4511",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6821067136421018"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def exp_coal(g, N):\n",
    "    \"\"\"\n",
    "    Compute expected coalescence time in epoch\n",
    "    N is the number of diploid invididuals\n",
    "    g is the number of generations spanned by the epoch\n",
    "    \"\"\"\n",
    "    return 2*N - (g * np.exp(-g/(2*N))) / (1 - np.exp(-g/(2*N)))\n",
    "\n",
    "def epoch(demog, h, i):\n",
    "    \"Recursively compute expected coalescence time across all epoches\"\n",
    "    g, N = demog[i]\n",
    "    N *= h\n",
    "    if i == len(demog)-1:\n",
    "        return 2*N\n",
    "    return (1-np.exp(-g/(2*N))) * exp_coal(g, N) + np.exp(-g/(2*N)) * (g + epoch(demog, h, i+1))\n",
    "\n",
    "def pool_nielsen(gens, Ne, h):\n",
    "    \"\"\"\n",
    "    Compute expected coalescence time in units of 2N\n",
    "    Ne is the a list/series of Ne in the epoque.\n",
    "    gens is the a list/series of generation at which an each epoque begins (the last epoque lasts forever)\n",
    "    h is the relative population size, 0.75 for chrX.\n",
    "    \"\"\"\n",
    "    epoques = list()\n",
    "    for i in range(len(gens)):\n",
    "        if i == 0:\n",
    "            epoques.append((gens[i+1], Ne[i]))\n",
    "        elif i == len(gens)-1:\n",
    "            epoques.append((None, Ne[i]))    \n",
    "        else:\n",
    "            epoques.append((gens[i+1] - gens[i], Ne[i]))\n",
    "    return epoch(epoques, h, 0)\n",
    "\n",
    "gen_time = 30\n",
    "gens = sampledemog_data.years / gen_time\n",
    "Ne = sampledemog_data.Ne\n",
    "\n",
    "pool_nielsen(gens, Ne, 0.75) / pool_nielsen(gens, Ne, 1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}