342 lines (341 with data), 79.0 kB

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting an Epigenetic Clock from Targeted Bisulfite Sequencing Methylation Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import stdlibs and third party libraries\n",
    "import os\n",
    "import subprocess\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "from scipy import optimize\n",
    "import scipy.stats as stats\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.linear_model import ElasticNet\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use latex formatting for figures\n",
    "rc('text', usetex=False)\n",
    "\n",
    "params = {'legend.fontsize': 'x-large',\n",
    "         'axes.labelsize': 'x-large',\n",
    "         'axes.titlesize':'x-large',\n",
    "         'xtick.labelsize':'x-large',\n",
    "         'ytick.labelsize':'x-large'}\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import Methylation Matrix\n",
    "- read methylation matrix into *pandas* dataframe\n",
    "- drop sites with missing data\n",
    "- drop sites that show no variation (sites not useful for model fitting)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in methylation data\n",
    "meth_matrix = pd.read_csv(f'TBS_methylation_matrix.tsv.gz', sep='\\t', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# age data for model fitting \n",
    "sample_ages =  {'S10': 37,'S11': 42,'S12': 45,'S13': 53,'S14': 38, 'S15': 46,'S16': 49,\n",
    "                'S17': 76, 'S18': 62, 'S19': 67, 'S1': 37, 'S20': 39,'S21': 45,'S22': 76,\n",
    "                'S23': 32, 'S24': 52, 'S25': 30, 'S26': 63, 'S27': 47, 'S28': 45,'S29': 32,\n",
    "                'S2': 64, 'S30': 33, 'S31': 52, 'S32': 54, 'S33': 57, 'S34': 47, 'S35': 34,\n",
    "                'S36': 46, 'S37': 51, 'S38': 34, 'S39': 34, 'S3': 44, 'S40': 38, 'S41': 43, \n",
    "                'S42': 32, 'S43': 44,'S44': 51, 'S45': 41, 'S46': 44, 'S47': 51, 'S48': 59,\n",
    "                'S4': 63, 'S5': 45, 'S6': 65, 'S7': 33, 'S8': 41, 'S9': 64}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop rows (CpG sites) with null data points \n",
    "meth_matrix = meth_matrix.dropna(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop row that have little variation \n",
    "meth_matrix = meth_matrix.loc[meth_matrix.var(axis=1) >= 0.005]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(16549, 48)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check size of final methylation matrix \n",
    "meth_matrix.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA\n",
    "- perform principal component analysis as a QC step\n",
    "- check for sample outliers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA(whiten=False)\n",
    "\n",
    "pcs = pca.fit_transform(meth_matrix.T.values)\n",
    "\n",
    "variance_explained = pca.explained_variance_ratio_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Close sample spread with little variation captured by 1st and 2nd PCs, keep all samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAucAAALBCAYAAAD78g/jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzde3icZ3nn8e/t0SSaABaHFCuOiJMgCIGWDbFIsUPBLtBAdhca11lgAVMaCFSm29ICbQNtgUJKCC0BaheynBoXStbGtFsOKZvEBuLgEokzmIMTYuw4MgFSmcQzyWj87B+SHFmRZkbWHF5rvp/rmssz7/uOdHOw/JtH93s/kVJCkiRJUvstancBkiRJksYZziVJkqSMMJxLkiRJGWE4lyRJkjLCcC5JkiRlRFe7C8iSk08+OZ1++untLkOSJEkL3PDw8M9SSr8y/bjhfIrTTz+doaGhdpchSZKkBS4i9sx03LYWSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScqIrnYXIEkAxXKFg8UyY5XDdOUWsbiQp5DPtbssSZJaynAuqe1GRkts2LabLcP7KJYrFPI51i7vY/3qfnp7uttdniRJLWM4l9RWI6Ml1mzcwf7R0pFjxXKFTTv3cMOuA2wdXElvT6GNFUqS1Dr2nEtqm2K5woZtu48K5lPtHy2xcfutlMqVFlcmSVJ7GM4ltc3BYpktw/uqXrN5aB+jxXKLKpIkqb0M55LaZqxymGKNVfFiucJY5XCLKpIkqb0M55Lapiu3qOZElkI+R1fOH1WSpM7gv3iS2mZxIc/a5X1Vr7l4oI+eQr5FFUmS1F6Gc0ltU8jnWL+6n6WzjEtc2tPN4Kp+up13LknqEIZzSW3V29PN1sGVrFux7EiLSyGfY92KZWwdPN8555KkjhIppXbXkBkDAwNpaGio3WVIHalUrjA6ZYfQnkLeFXNJ0oIVEcMppYHpx92ESFImdOdzhnFJUsezrUWSJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIzIbDiPiNsjIs3w+O4c33NTK+uWJEmSjlVXuwuo4qlAbsrrhwLfAj5Z431XAFdNeX1/g+uSJEmSmiKz4TyldNfU1xHxKiAPfKjGW+9JKY00rTBJkiSpSTLb1jKDVwP/llK6s8Z1r42In0fEdyPifRHxqFYUJ0mSJM1XZlfOp4qIAWA58KYal74f+CZwAHgC8Hbggog4J6VUnOVrXwpcCnDaaac1rGZJkiRpriKl1O4aaoqI/w08C3hsmkPBEXEmsBt4aUrpE7WuHxgYSENDQ8deqCRJklSHiBhOKQ1MP575tpaIWAy8GLh6LsEcIKV0G+Or6Kc3oTRJkiSpoTIfzoGXAicAH53rGyPiVODRwN5GFyVJkiQ12vEQzl8N/EtK6cDUgxFxUUR8fyKAExErIuL1EXFuRCyLiAuAzwA/AT7d+rIlSZKkucl0OI+IpwFPBj44w+ke4CzGxysC3AesAa4HfghsBHYCK1JK9zS/WkmSJGl+Mj2tJaW0E4hZzn0M+NiU118DVrakMEmSJKkJMr1yLkmSJHUSw7kkSZKUEYZzSZIkKSMM55IkSVJGZPqGUEmSpONRsVzhYLHMWOUwXblFLC7kKeRz7S5LxwHDuSRJUgONjJbYsG03W4b3USxXKORzrF3ex/rV/fT2dLe7PGWc4VySJKlBRkZLrNm4g/2jpSPHiuUKm3bu4YZdB9g6uJLenkIbK1TW2XMuSZLUAMVyhQ3bdh8VzKfaP1pi4/ZbKZUrLa5MxxPDuSRJUgMcLJbZMryv6jWbh/YxWiy3qCIdjwznkiRJDTBWOUyxxqp4sVxhrHK4RRXpeGQ4lyRJaoCu3KKaE1kK+RxdOeOXZuf/OyRJkhpgcSHP2uV9Va+5eKCPnkK+RRXpeGQ4lyRJaoBCPsf61f0snWVc4tKebgZX9dPtvHNVYTiXJElqkN6ebrYOrmTdimVHWlwK+RzrVixj6+D5zjlXTZFSancNmTEwMJCGhobaXYYkSTrOlcoVRqfsENpTyLtirqNExHBKaWD6cTchkiRJarDufM4wrmNiW4skSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRmR2XAeEW+JiDTDo7/Ke/IR8a6IuDMiihFxU0Qsb2XdkiRJ0rHKbDifcDtwyrTHj6tcfyVwCfBq4KnAbcD1EdHb3DIlSZKk+ct6OK+klEamPSozXRgRi4HXAH+eUvq/KaXvAK8A7ps4LkmSJGVa1sN5X0Tsm3h8PiJWVrl2OXAicN3kgYkg//+Apze5TkmSJGneshzOv8r4yvd/BV4M/Bz4ckQ8Z5brT5n4c2Ta8ZEp5x4kIi6NiKGIGLrrrrvmWbIkSZJ07LraXcBsUkqfm3boyxHRB7yB8dXwRn2fq4GrAQYGBlKjvq4kSZI0V1leOZ/JV4DTZzl358Sf02/+XDLlnCRJkpRZx1s4PxfYO8u5YcZv/rxg8kBELAKeDdzU/NIkSZKk+clsOI+Iv4uI34yIMyPinIjYADwHuGri/EUR8f2IOBUgpXQQ+ABweUT8t4h4EvARoAB8sE3/MSRJkqS6ZbbnnPGbOK8BfgUYBb4FPDuldOPE+R7gLCA/5T1vAO4HPgQ8nPHV9OeklGxrkSTpGBTLFQ4Wy4xVDtOVW8TiQp5CPtfusqQFK1LyHshJAwMDaWhoqN1lSJKUCSOjJTZs282W4X0UyxUK+Rxrl/exfnU/vT3d7S5POq5FxHBKaWD68SyvnEuSpDYZGS2xZuMO9o+Wjhwrlits2rmHG3YdYOvgSnp7Cm2sUFqYMttzLkmS2qNYrrBh2+6jgvlU+0dLbNx+K6XyjJt2S5oHw7kkSTrKwWKZLcP7ql6zeWgfo8VyiyqSOofhXJIkHWWscphijVXxYrnCWOVwiyqSOofhXJIkHaUrt6jmRJZCPkdXzhghNZp/qyRJ0lEWF/KsXd5X9ZqLB/roKeSrXiNp7gznkiTpKIV8jvWr+1k6y7jEpT3dDK7qp9t551LDGc4lSdKD9PZ0s3VwJetWLDvS4lLI51i3YhlbB893zrnUJG5CNIWbEEmSdLRSucLolB1Cewp5V8ylBnATIkmSNGfd+ZxhXGoh21okSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScqIrrm+ISLOAZYBt6eUvtn4kiRJkqTOVHXlPCLeHBHPmXj+yIj4IvA14NPA1yLihoh4eAvqlCRJkha8Wm0trwZ+NvH8ncDDgacCJwHLgYdNHJckSZI0T7XaWh4N/Hzi+bOAV6aUhidefyMiXgt8qlnFSZIkSZ2k1sr5fuDxE8/zwD3Tzt8LPLLRRUmSJEmdqFY4vxb4m4g4GbgG+OuI6AGIiIcBbwd2NrdESZIkqTPUamt5K/DrwG7gq8BvACMRcSdwCnAI+M2mVihJkiR1iKor5ymlIuO95q8D7gN+DOyZeFwJPMlxipIkSVJj1JxznlI6DHx04iFJkiSpSdwhVJIkScqImuE8Iv4gIj4VEesnXr8iIvZHxC8i4u8iIppfpiRJkrTwVW1riYjXAX8NfAH4i4ndQP8EeBeQA94I/BD4QJPrlCRJkha8Wj3nrwJ+P6W0KSLOA24CXptSuhogIn4KXIrhXJIkSZq3Wm0tpwNfBkgpfRVIwM1Tzn8ROLMplUmSJEkdplY4LwEnTnl9P+OzzSeVge5GFyVJkiR1olrh/DbgrCmvT2V81vmkM4E7Gl2UJEmS1Ilq9ZxfBRQnX6SUDk47/xvAZxpdlCRJktSJqobzlNI/1Tj/loZWI0mSJHUwNyGSJEmSMqKeTYh+NSI+EhE/jIhfTjx+OHHs11pRpCRJktQJam1C9Bzg34BvAZ8ERiZOLQGeC3w1Iv57Sun6plYpSZIkdYBaN4ReAbw3pfSnM5z7q4h4J+O7hZ7b8MokSZKkDlOrreWJwEeqnP8ocHbjypEkSZI6V61wPgIMVDk/ABxoXDmSJElS56rV1rIRuDoingD8Ow/0nPcCFwCvA97avPIkSZKkzlFrzvm7IuJ+4PXAm4A0cSqA/cCbUkrvbW6JkiRJUmeotXJOSukq4KqIOB04ZeLwnSml25tXliRJktR5aobzSRNh/PamVSJJkiR1uHntEBoRAxHxjEYVI0mSJHWyulfOZ7EJeDyQa0AtkiRJUkebbzhfB5zUiEIkSZKkTjevcJ5SuqVRhUiSJEmdru5wHhFLOXpay/7mlCRJkiR1ppo3hEbEpRFxG7AX+CpwC7A3Im6LiFc1u0BJkiSpU1RdOY+I1wGXA+8HruOBHUKXAM8D3hsRJ7kRkSRJkjR/tdpa/gB4dUrpmmnHvwdsi4jvAm8BDOeSJEnSPNVqa1kK/EeV8zuB3saVI0mSJHWuWuF8F/CyKudfBny/ceVIkiRJnatWW8ubgH+JiFUc3XPeC1wAnAf8dtOqkyRJkjpI1XCeUvpcRKwEXge8kimjFIEdwOtSSkPNLVGSJEnqDDXnnE+E75e0oBZJkiSpo9Wccy5JkiSpNQznkiRJUkYYziVJkqSMMJxLkiRJGWE4lyRJkjIis+E8It4QEV+JiLsj4j8j4qaIeG4d77s9ItK0x02tqFmSmqVYrnDgYIk77j7EgYMliuVKu0uSJDVBzVGKkyLi28CFKaW9U583rzR+E/gIcAtwiPE565+JiGemlHbUeO8VwFVTXt/fnBIlqflGRkts2LabLcP7KJYrFPI51i7vY/3qfnp7uttdniSpgeoO58DpQH6G502RUnretENvnFg5X8P4BkjV3JNSGqlxjSRl3shoiTUbd7B/tHTkWLFcYdPOPdyw6wBbB1fS21NoY4WSpEbKbFvLdBGxCFgM3FvH5a+NiJ9HxHcj4n0R8agmlydJDVcsV9iwbfdRwXyq/aMlNm6/lZItLpK0YBw34Ry4DHg4cHWN694PvBRYBbwFuAC4OSJmXFqKiEsjYigihu66667GVStJ83SwWGbL8L6q12we2sdosdyiiiRJzTaXtpa2iYhBxsP581NKVf+lSin97ZSX346IYWA3cBHwiRmuv5qJwD8wMJAaVrQkzdNY5XDNGz+L5QpjlcMtqkiS1GyZXzmPiNcDVzIezK+f6/tTSrcBBxjvk5ek40ZXbhGFfK7qNYV8jq5c5n+US5LqlOmf6BHxNuCvGJ8MM+dgPvE1TgUeDTRzsowkNdziQp61y/uqXnPxQB89habeny9JaqFjDedNb/+IiKuANwAvA34QEb0Tj54p11wUEd+fCOBExIqIeH1EnBsRyyLiAuAzwE+ATze7ZklqpEI+x/rV/SydZVzi0p5uBlf1011jdV2SdPw41nAeDa1iZn8IdDMequ+c8njvlGt6gLN4YKzjfYyPWrwe+CGwEdgJrEgp3dOCmiWpoXp7utk6uJJ1K5YdaXEp5HOsW7GMrYPnO+dckhaYSKm+RfCIeAxwR0rp8MTz/SmlBTW/a2BgIA0NDbW7DEl6kFK5wmixzFjlMF25RfQU8q6YS9JxLCKGU0oD04/XPa1l6m6gTd4ZVJI0TXc+ZxiXpA6Q6RtCJUmSpE5iOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEbNOa4mIy+r9IimlyxtTjiRJktS5qo1SfNW0178CnASMTrzuAQ4BPwUM55IkSdI8zdrWklI6Y/IBvBH4HvBrKaVHpJQeAfwa8B3gz1pTqiRJkrSw1dtz/nbgD1JK3508MPH8j4B3NKMwSZIkqdPUG86XAffOcPwQ0Ne4ciRJkqTOVW84/xrwNxGxePLAxPN3TJyTJEmSNE/Vbgid6veBzwL7ImLXxLGzgXuA5zWjMEmSJKnT1BXOU0rfjIh+4KWMh3KAq4FPpJSKzSpOkiRJ6iT1rpyTUioBH2piLZKkCcVyhYPFMmOVw3TlFrG4kKeQz7W7LElSk9UdziNiNfAHQD9wYUppX0RcAtyaUtrepPokqeOMjJbYsG03W4b3USxXKORzrF3ex/rV/fT2dLe7PElSE9V1Q2hEXAR8HrgbeDxwwsSpAuMz0CVJDTAyWmLNxh1s2rmHYrkCjK+ib9q5hzUbdzAyaiehJC1k9U5reTPw2pTSJUB5yvGbgXMaXpUkdaBiucKGbbvZP1qa8fz+0RIbt99KaSK0S5IWnnrD+ROA62c4fjfwyMaVI0md62CxzJbhfVWv2Ty0j9Fiueo1kqTjV73h/G7glBmOPxm4o3HlSFLnGqscPtLKMptiucJY5XCLKpIktVq94fxTwDsi4mETr1NEPBG4Ari2KZVJUofpyi2qOZGlkM/Rlav3R7ck6XhT70/4y4AADgAnAUPAt4E9wFubU5okdZbFhTxrl/dVvebigT56CvkWVSRJarV6NyG6F1gdEauAAcZD/VBK6cYm1iZJHaWQz7F+dT837Dow402hS3u6GVzVT7fzziVpwap7zjnAxDzz7U2pRJJEb083WwdXsnH7rWweemDO+cUDfQyucs65JC10kVKq78KIfmA1sIRp7TAppbc1vrTWGxgYSENDQ+0uQ5IolSuMTtkhtKeQd8VckhaQiBhOKQ1MP17XynlEvBz4EFAEfgpMTfQJWBDhXJKyojufM4xLUgeqt63lL4F3AX+RUnKGlyRJktQE9U5r6QU+bDCXJEmSmqfecH4j8JRmFiJJkiR1unrbWq4BroiIPuCbwP1TT6aUbm50YZIkSVKnqTecT+4C+p4ZziXAu5YkSZKkeao3nJ/R1CokSS1TLFc4OGVM4+JCnoKTYSQpE+rdIXRPswuRJDXfyGiJDdt2s2X4gQ2O1i7vY/1qNziSpCyYNZxHxErgP1JKlYnns7LnXJKyb2S0xJqNO9g/WjpyrFiusGnnHm7YdYCtgyvp7Sm0sUJJUrWV85sYH6H404nnCYgZrrPnXJIyrliusGHb7qOC+VT7R0ts3H4rl114tpsfSVIbVQvnZwB3TXkuSTpOHSyW2TK8r+o1m4f2sX51v+Fcktpo1nA+tc/cnnNJOr6NVQ5TLFeqXlMsVxiruNecJLVTvdNaAIiIU4BlwAlTj6eUvtTIoiRJjdWVW0Qhn6sa0Av5HF25evemkyQ1Q13hPCJ6gX8GnjF5iPFe80n+DlSSMmxxIc/a5X1s2jn7L0IvHuijp5BvYVWSpOnqXSJ5D5AHBoAi8GzgZcAPgec1pzRJUrFc4cDBEnfcfYgDB0s1W1NmU8jnWL+6n6WzjEtc2tPN4Cr7zSWp3epta1kF/HZK6esRcRjYm1K6MSIOAW8GvtCsAiWpUzV6JnlvTzdbB1eycfutbB564GtePNDH4CrnnEtSFkRKqfZFEb8EfjWltCcifgK8MKX0lYg4HfhuSukhzS2zNQYGBtLQ0FC7y5CkGWeST1o6EbKPdSZ5qVxhdMoOoT2FvCvmktRiETGcUhqYfrzetpbdwGMnnn8PeFlEnAi8BPhZY0qUJEH9M8lLx9ji0p3PsWRxN6c+4iSWLO42mEtShtQbzj8KPGni+TuBlwOHgLdOvJYkNUi9M8lHi+UWVSRJapW6es5TSu+b8nx7RDwBeCrwo5TSt5tVnCR1ImeSS1LnmtOc80kppb3A3gbXIkmi8TPJi+UKB6f0mC8u5CnYyiJJmTRrOI+I/1nvF0kpfaIx5UiSGjmTvNETXyRJzTXrtJaJkYn1SCmlBbEE47QWSVlRe1rL+TXDdTMnvkiS5mfO01pSSovqfCyIYC5JWTI5k3zdimVHWlAK+RzrViyrK5g3e+KLJKk5jqnnXJLUfL09BS678GzWr+6f80zyeie+rF/trqCSlCV1h/OI+C/AH/PASMXvAX+bUvpmMwqTJI3PJD+W8OzEF0k6PtV1q39EvBAYBvqBGyceZwLDE+d0jIrlCgcOlrjj7kMcOFiq+Y+pJNVjcuJLNXOZ+CJJao16V84vB/4mpfQXUw9GxNsmzl3b6MI6gVMUJDVLIye+SJJap94lk1OAa2Y4vmninOZocorCpp17jqyWF8sVNu3cw5qNOxgZLba5QknHs0I+x/rV/Syd5YP+0p5uBlfZby5JWVNvOL8ZWD7D8QHgPxpXTmdwioKk6ZrR4jbfiS+SpNart63lo8C7I+IJwM6JY08DLgH+LCJWTl6YUrq5sSUuPE5RkDRVM1vc5jPx5Vi4G6kkzU+94XzTxJ9/WeUcQAL8KVyDUxQkTZppo6DJFrcbdh1oyEZBxzrxZa68j0aS5q/ecH5GU6voMJNTFKoFdKcoSAtfvS1ul114duZ/i9aKDxmS1AnqSn8ppT2zPYC9016rhskpCtU4RUFa+OptcRstlltU0bHxPhpJapx655z/U0ScNMPx04AvNbyqBc4pCpJg4bS4LZQPGZKUBfX2TZwDfCMizp08EBFrgW8CB5tR2ELnFAVJrdwoqJkbni2UDxmSlAX19pwPAO8Hbo6INwNnAeuAN6WU3t2s4ha6Vk9RkJQtrdooqNk3anofjSQ1Tl3hPKVUAl4VEfuBdwFjwLNTSra0zFOrpihIyp7JFrcbdh2YsV+7ES1u87lRs96xiO5GKkmNU/cyRkQMAm8A/g9wO7AxIp7UpLokqSM0s8VtPjdqjoyWuPyzu1h15XbOv2Ibq67czuWf3cXIDF/L+2gkqXEipVT7oohPA88CXptSuiYiHgJ8AFgDvD6l9A/NLbM1BgYG0tDQULvLkNSBSuUKo1NWqRvR4nbgYIlVV26v2W6y/Q2rWLL4gWA902r7pKUTHyZmWm0fGS2ycfutbB56oH3m4oE+Blc551ySpouI4ZTSwPTj9facLwOWp5R+BJBSuhd4WUR8gfFe9AURziWpXZrR4nYsN2rOZ/a699FI0vzVG86fllK6f/rBlNKmiLi5wTVJs3JrcKl+x3KjZr1jEdevnrlNxftoJGl+qobziHhkSukXMwXzifM5oKcplUnTuDW4OkEjP4AuLuR5ya+fRvlw4rlP6mVRwOEE1313hK3D+/jlfWMPulHTsYiS1F61Vs7viohTUko/BYiIrwP/PaU0uaxyMnAL4DKJmsqtwdUJ5voBtFaQL+Rz/N7Tz+B9N/yI3/vYLUe+5vPPWcqmS87j3f/+gwfdqOlYRElqr1rhPKa97gdOqHGN1FDz6YGVjhdz/QA6Nch3LQpe+NTHcNG5p/KwE7s4MZ9jcSHP6KEya//h5gd9zWtv2cuXf3gXm1+z8kGh37GIktRejVj6qD3uRZoHtwbXQjfXkYeTQX7Tzj085pEFNl1yHr+8b4y1//AVnnHldlZduZ0vfGek5tf84JcePEbRsYiS1F7+XlKZZw+sFrq5fACdGuQfdmIX7/qdJzP48a9x7S17j/w9KZYrPHpx9zF/qJ3L7PViucKBgyXuuPsQBw6Wav5dlSRVV6utJXH0yvj011LT2QOrhW4uH0CnBvk1y/v451v2zrg6viiY14faesYiepO2JDVePT3nOyNi8if8Q4AbI2JyqcXfa6rp7IHVQjeXD6BTg/xzn9TL733slhmvP5yY94faamMRvUlbkpqjVjh/a0uqkKqY7IG9YdeBWXcstAdWx7O5fAAdLZaPhO5qq+PXfXeE55+zlGtv2Vvza86VN2lLUvNUDecpJcO5MmGyB9atwbUQzeUDaIIjQb7a6vjW4X1suuQ8vvzDuxr+oXa+GxVJkmYXKdlCPmlgYCANDQ21uwxVUSpXGJ0y19mtwbWQjIwW6/oAOtlS8pwn9VKaGI04k8cveSjvf/G5bBneS0rwrLOXkFsEDzmxi0c/tJuHFbqO+vtU74ZHd9x9iPOv2Fbzuh1/uppTH3FS/f8FSFIHiYjhlNLA9OO12lraKiIuBC4HzgbuBN6XUvq7Ot73RmA90At8D/jTlNIXmlmrWsOtwRqw3Y0AACAASURBVBe+Ru6Qebyp5ybM8evGf5P04Zt+zEVPOW3W1fF7SmM86iEncMnTz2DDtluP2oho7fI+Lnn6GfzZp77Fzh//Yk43c3qTtiQ1T2ZXziNiAPgK8G7gGuDXgQ8Af5RS+kCV9/0R8DfAqxnfvfQVwB8CT00pfava93TlXGovp3/MTalc4Z77xrivXOGDX7yNzcNHr7ivX9VPggfduDlpaU83n3z1Cv7XP3+Nb+wdPXKs1s2cxXKFyz+7q2qP/LoVy+w5l6QqZls5z3I4/wRwekpp5ZRjVwIXp5ROn+U9AewD/jGldNmU47cA300p/W6172k4l9pnpukfk+oJjJ1uppavBDVD9Iuf+hhe9Ywzef7f7+Ce+8aA+oJ17f+9zvcDlSRVMVs4z/LvHM8Hrpt27DpgWUT0zfKe04Gls7zv6Q2tTlLDzHWHzIVovpv5dOdzLFnczamPOIkli7vpzufqunHzX76xn+L9FS59xplHjm0e2sfP77m/ah1z2ahIklS/unvOI+IsxltFHge8KqU0EhHPB/aklL7ZhNpOAUamHRuZcm6mf3FOmXbd1Pedwgwi4lLgUoDTTjvtmAqVND+dPv2jWe089W5udO/9Y7zgnKVc/aXbCMY3Nxo7fJg7/rPEQ07s4uSHnkjXouCkE7uO6v+vt0deklS/usJ5RPwG8AXgZsZXoCdvv38i8HLgd5pSXQuklK4GrobxtpY2lyN1pLnskLnQHMtmPvXeNFvvjZtjhxO/LI1x6TPO5BmPO5l/vmUvz73qy0c+KLzgnKWsX93PluF9vOCcU4/6wOBN2pLUWPW2tVwOvC2l9Czg/inHbwTOa3hV4+5kfNrKVEumnJvtPczyvtneI6nNJkNkNQtx+sextPOMjJa4/LO7WHXlds6/YhurrtzO5Z/dxcgMX2Nyc6NqXnDOUm743gF+WSrz2+csZfDjX+PaW/YeCfTFcoVP3rKXF37wK6w669G89EM7GRktzuM/tSSpmnr/pfsvwLUzHD8A/ErjyjnKDuCCaceey3gbzWy//74d2D/L+25qaHWSGqaeEHmsu1lmWb3tPKPFMvDAKvumnXuOCs+bdu5hzcYdDwrNhXyOwdWPZeksrTFLe7p50VMfw9jhxOe/M8LB0ljVDwqbdt7Oyv6TF3z/vyS1U73hvAT0zHD88cBdjSvnKO8BzouId0TEEyLi5cAfAO+cvCAizouI70fEeQBpfPTMlcDrIuKlE+97J+MfLt7TpDolzdPkDpnVQuSx7maZZXNp5znWm2ZP6Smw+TUrefFTH3PUjZsveupj2PiSc3n3v/+Ai55yKp/79p38slSuWsunhu/gub/ae9QHBklSY9V7Q+jngD+PiBdNvE4RcTLwduD/NqOwlNItEfHbjLfUvJ7xmzrfNG3G+UnAWTzQA09K6aqIOHHifUuAXcDzm3TTqqQGmZz+Uc8OmQtFvT3hEcEv7rnvmG+afdRDT2BwdT8vedoyDt0/xtjhxA3fO8Cnv34Hr7/gLN74qW/xvF/t5fPfmX4v/dGK5Qq5iAXb/y9JWVBvOH8jsI3xtpFu4F+AM4EfA29uSmVASumzwGernN8OxAzHrwCuaFZdkpqj06Z/TLbzVJtDvuYpp/LhL9/Gbz2p95hvmu3O5zght4jrvnMnK/tPJhfBs564hOu+M8JLP/xVFnd3cekzH8tFG3ZU/fqFfI5KSguy/1+SsqKucJ5S+mlELAdeBAww3g7zXuDjKaX7mlifFrBO3qZds+uk6R+T7Tw37Dow62Y+Fw/08dIPf5VnP7G3rlX22ULzkp5uXvq0ZbP+ZuKE3CKeffYSPnnL3lm//gvOWcp13xlZkP3/kpQVmd0htB3cIbR13KZdesDIaPFBoXnNU07l4oE+3vipb/HDA/fw8pWnUypXuLZKeK5nZ8+ZdhKdvH7/fxZZ+w83z/pBYeNLzuX1m7/JP73yafSclPfDtSTNw2w7hNYVziPiz4C7Ukofnnb8EuBRKaV3NazSNjKct4bbtEsPNjU0RwQf/vJtXDu0j3vuGwPgYSd2semS8xj8+Neq/N2Z/86cd/5nkX/Yfiubp3xwfsE5S3nReadx04/uYu3yx9C1KHjvDT/yw7UkzcN8w/ltwLqU0k3Tjq8ENqWUHtuwStvIcN58xXKFyz+7q2qPbT2rf9JCdsfdhzj/im0POv74JQ/lit95Mtfespd//cb+pt00O/lBoVw5zKIIuhYBBIsLeUYPlbnID9eSNG+zhfN6bwhdCsw0JmA/cOp8ClNn6fRt2qV6zDbF5YcH7uFlH/4qa849lY/87gCnPqLAiV25Y7pptto9H7P1/RfLFf6+jnGOfriWpGNXbzj/KfBrjE9rmerJwM8bWZAWtk7epl2qV7UpLvfcN8Y1Xxk/fqwh+Fjv+fDDtSQ1X72zsLYC74mIp0weiIhzgb8FtjSjMC1MnbpNuzQXzdyUaa67jE7lh2tJar56E9CbGG9rGYqIn0XEz4BbGG9ruaxZxWnh6dRt2qW5mtyUad2KZUft7LluxbJjvvHzWHcZneSHa0lqvnrnnN8LrIqIZwHnThweTind2LTKtCDVM9d5IW7TLh2LRm/KNN+2lHo2TfLDtSTNT7095wCklG4AbmhSLeoQnbhNu3SsGrkp03zbUvxwLUnNV3c4j4h+YDWwhGntMCmltzW4Li1wnbZNu9QuU6eyRMS8dhkFP1xLUrPVO+f85cCHgCLjk1umvimllB7fnPJayznnkhaS6VNZXvn0Mzh0/xif+Gr1XUbfcMFZHLq/UnX3z2o7jUqSapvvnPO/BN4F/EVKydvwJSnjZtqJ99pb9rLpkvPY/oO7Zm1Lec0zH8t7r/8RH/+Pn1Qds9jIdhtJ0gPqvaW+F/iwwVySsm+2qSy/vG+MN37qW2x4ybm85LzTHjQFZsvvr+SNm7/Jh2768ZzHLEqSGqPelfMbgacAtzWxFklSA1SbyjK5y+gLB/q44U+eSUqJrtwiTjohx3uv/xE33TrzvnLu/ilJrVHvyvk1wBUR8YcRsSoiVk59NLNASdLc1JrKcs99Y3x4x+2klDj1ESexZHE3h+6v8PH/+EnVr7t5aB+jxXKjy5UkTVHvyvm1E3++Z4ZzCXAZRZIyYnKzoLlMZXH3T0nKhnrD+RlNrUKS1DDHslnQsQR6SVLj1fVTNqW0p9qj2UVKkuo3uVnQ0llmjs+0WdBkoK/G3T8lqfnmsglRF/BUYBlwwtRzKaVrGlyXJGke5rpZkLt/SlI21LsJ0eOAzwL9PLABUQCHgcMppRNme+/xxE2IJC00c90saGS06O6fktQCs21CVG84/wwwBvwesAcYAB4FXAX8SUrpy40ttz0M51K2Td2KfradKzV/7v4pSc033x1Cfx14dkrpFxGRAFJKN0fEnzMe0Jc3rlRJerDpW9HPtnOl5s/dPyWpfeq97T4PjE48/xmwZOL5j4GzG12UJE01uRX9pp173LlSkrSg1RvOvw88aeL5N4D/FRFnAX8C7G1GYZIEs29FP2ly58pSjRndkiQdD+oN5+8FTp54/jbgN4DvAa8A/rwJdUkSUH0r+knuXHn8KZYrHDhY4o67D3HgYKnmBkiS1Cnq6jlPKf3zlOffiIjTGW9n2ZNS+nlzSpMkd65ciLx/QJJmV/ec86lSSkXgaw2uRZIexJ0rF5bJ+wemtilN3j9ww64DbB1cSW9PoY0VSlJ7zRrOI+Iy4D0ppeLE81mllC5veGWSxLFtRa9sqvf+gcsuPNtpMZI6VrWV81cBHwSKwKU8sPnQdAkwnEtqCneuXDjqvX9g/Wr/95TUuWYN5ymlM6Y8P70l1UjSDOa6Fb2yyfsHJKm2mj3nEZEHvgS8IqX0/eaXJEkP1ttT4LILz2b96n53rjxOef+AJNVW8ydgSqkMPA4Ya345kjS77nyOJYu7OfURJ7FkcbfB/Dgzef9ANd4/IKnT1bs8cS3wP5tZiCRpYZu8f2DpLG1I3j8gSfWPUvwF8McR8RvAV4F7p550WoskqR7ePyBJ1UVKsw1hmXJRxI+rnE4ppTMbV1L7DAwMpKGhoXaXIUkLXqlcYbRY9v4BSR0rIoZTSgPTj9e7Q+gZta+SJKk+3fmcYVySZuAt8ZIkSVJG1NtzTkT0AxcDy4ATpp5LKf1eg+uSJEmSOk5d4TwiLgD+Ffg+8ETgm8CZjK+839K06iRJkqQOUm9by18D70opnQPcB7wQOI3xzYm2Nqk2SZIkqaPUG87PBq6ZeD4GFFJK9wJ/BbyxGYVJkiRJnabecH6IB1pgRoDTJ56PAUsaXJMkSZLUkeq9IXQYOI/xnvNtwOUR0Qe8BPh6k2qTJEmSOkq9K+dvAvZMPP9LYB9wJVAAXt2EuiRJkqSOU3XlfGJn0A8BH0kpfR0gpfQz4L+2oDZJkiSpo9RaOf888HpgT0RsjYjntqAmSVKGFMsVDhwsccfdhzhwsESxXGl3SZK0YFVdOU8pDUbEnzA+OvGVwOciYg8PrKbf2YIaJUltMjJaYsO23WwZ3kexXKGQz7F2eR/rV/fT29Pd7vIkacGp2XOeUiqmlD6WUno68CTg08AfMb6a/mlX0yVpYRoZLbFm4w427dxzZLW8WK6waece1mzcwchosc0VStLCU+8NoQCklHallP4YWAq8AngG8NlmFCZJap9iucKGbbvZP1qa8fz+0RIbt99KyRYXSWqoOYVzgIg4E3gb8G6gB/j3RhclSWqvg8UyW4b3Vb1m89A+RovlFlUkSZ2hrjnnEZEH1gCvAlYBBxjvO/9QSmlPlbdKko5DY5XDNW/8LJYrjFUOt6giSeoMtUYpPoHxQP4y4JHAF4C1wL+llPxdpiQtUF25RRTyuaoBvZDP0ZWb8y9gJUlV1Pqp+j3gRcAHgcemlC5MKf2LwVySFrbFhTxrl/dVvebigT56CvkWVSRJnaFWOF8DnJZS+gvbVySpcxTyOdav7mfpLOMSl/Z0M7iqn+58rsWVSdLCVjWcu0ouSZ2rt6ebrYMrWbdiGYWJEF7I51i3YhlbB893zrkkNUGklNpdQ2YMDAykoaGhdpchSZlSKlcYLZYZqxymK7eInkLeFXNJmqeIGE4pDUw/Xte0FklS5+rO5wzjktQi3mYvSZIkZYThXJIkScqImuE8xj0pIs6Y4Vw+Ip7RnNIkSZKkzlI1nEfEUuAbwLeB3RFxbUQ8bMoljwS2NbE+SZIkqaGK5QoHDpa44+5DHDhYqrkjcivVuiH07UAJWA48Avhb4PqIeE5K6eDENdHE+iRJkqSGGRktsWHbbrYM76NYrlDI51i7vI/1q/szMSK2VlvLc4DXpZS+nlK6EXg6cC9wXUQ8ZOIaZzFKkiQp80ZGS6zZuINNO/ccWS0vlits2rmHNRt3MDJabHOFtcP5o4A7J1+klO4F/htwGPgccFLzSpMkSZIao1iusGHbbvaPlmY8v3+0xMbtt1Jqc4tLrXC+D3jC1AMppUPAhcBDgf/TpLokSVLGZLlPV6rlYLHMluF9Va/ZPLSP0WK5RRXNrFbP+Y3Ai4HPTz2YUjoYEb8FfLlZhUmSpOzIep+uVMtY5XDND5TFcoWxyuEWVTSzWuH8HcBZM51IKf08In4TuKDhVUmSpMyY7NOd2g4w2ad7w64DbB1cSW9PoY0VSrV15RZRyOeqBvRCPkdXrr3bAFX97imlvSml66ucH0kp/WPjy5IkSVlwvPTpSrUsLuRZu7yv6jUXD/TRU8i3qKKZ1Zpz/pSI+GJE9Mxw7uERsT0inty88iRJUjsdL326Ui2FfI71q/tZOksb1tKebgZX9dOdz7W4sqPVWrf/Y+CLKaXR6SdSSv8JbJ+4RpIkLUDHS5+uVI/enm62Dq5k3YplFCZCeCGfY92KZWwdPD8T90/U6jl/GvCeKuf/FdjcuHIkSVKWHC99ulK9ensKXHbh2axf3c9Y5TBduUX0FPJtXzGfVOtvUh/w8yrnfwGc2rhyJElSlhwvfbrSXHTncyxZ3M2pjziJJYu7MxPMoXY4HwUeW+X8YyeukSRJC9Dx0qcrLRS1wvnNwCurnH/VxDWSJGmBOh76dKWFolbP+ZXAlyLibuAdKaX9ABGxFHgzsBZ4ZnNLlCRJ7Zb1Pl1poagazlNKX4mIS4GNwGsi4uDEqcXA/cBrUkqunEuS1AG68znDuNRktVbOSSl9NCL+HfgfwOOAAH4AbJ5cSW+GiHgDsAZ4wsT3/A7w9pTSdTXedzuwbNrhHSmlpzejTkmSJKlRaobzCT8H/ndK6d5mFjPNbwIfAW4BDjHe+/6ZiHhmSmlHjfdeAVw15fX9zSlRkiRJapyq4TwiHgn8I/BcYFFE7AReklK6vdmFpZSeN+3QGyPiuYyvptcK5/eklEaaU5kkSZLUHLWmtbwD+HXgr4A3AEuADzS7qJlExCLGe93rWb1/bUT8PCK+GxHvi4hHVfm6l0bEUEQM3XXXXQ2rV5IkSZqrWm0tzwMuSSn9G0BEXAd8KyK6UkpjTa/uaJcBDweurnHd+4FvAgcY71d/O3BBRJyTUipOvzildPXk1xwYGEgNrViSJEmag1or56cCw5MvUkrfY7x/e+mxfLOIeEtEpBqPt8zwvkHGw/nalNK+at8jpfS3KaXrU0rfTiltZvwDxuOAi46lZkmSJKlVaq2c54DytGOViePH4u+BT9a45mdTX0TE64G3As9PKV0/12+YUrotIg4Ap8/1vZIkSVIr1TOtZXNETJ120g1cExFHWkRSSr9VzzdLKf2MaeG7moh4G/A64MKU0hfrfd+0r3Eq8Ghg77G8X5IkSWqVWuH8H2c49k/NKGS6iLgKeDXwYuAHEdE7caqYUhqduOYi4G+AZ6WU7oiIFcD5wI2Mj398AvBO4CfAp1tRtyRJknSsau0Q+opWFTKDP5z4c3qo/kfgdyee9wBnAfmJ1/cxPmrxMuAhwD7gC8BbU0r3NLNYSZIkab7q3YSo5VJKUcc1HwM+NuX114CVzatKkiRJap5a01okSZIktYjhXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjOhqdwGSpGwqliscLJYZqxymK7eIxYU8hXyu3WVJ0oJmOJckPcjIaIkN23azZXgfxXKFQj7H2uV9rF/dT29Pd7vLk6QFy3AuSTrKyGiJNRt3sH+0dORYsVxh08493LDrAFsHV9LbU2hjhZK0cNlzLkk6oliusGHb7qOC+VT7R0ts3H4rpXKlxZVJUmcwnEuSjjhYLLNleF/VazYP7WO0WG5RRZLUWQznkqQjxiqHKdZYFS+WK4xVDreoIknqLIZzSdIRXblFNSeyFPI5unL+8yFJzeBPV0nSEYsLedYu76t6zcUDffQU8i2qSJI6i+FcknREIZ9j/ep+ls4yLnFpTzeDq/rpdt65JDWF4VySdJTenm62Dq5k3YplR1pcCvkc61YsY+vg+c45l6QmipRSu2vIjIGBgTQ0NNTuMiQpE0rlCqNTdgjtKeRdMZekBomI4ZTSwPTjbkIkSZpRdz5nGJekFrOtRZIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZYThXJIkScoIw7kkSZKUEYZzSZIkKSMM55IkSVJGGM4lSZKkjDCcS5IkSRlhOJckSZIywnAuSZIkZURXuwuQWqFYrnCwWGascpiu3CIWF/IU8rl2lyVJknQUw7kWvJHREhu27WbL8D6K5QqFfI61y/tYv7qf3p7udpcnSZJ0hOFcC9rIaIk1G3ewf7R05FixXGHTzj3csOsAWwdX0ttTaGOFkiRJD7DnXAtWsVxhw7bdRwXzqfaPlti4/VZK5UqLK5MkSZqZ4VwL1sFimS3D+6pes3loH6PFcosqkiRJqs5wrgVrrHKYYo1V8WK5wljlcIsqkiRJqs5wrgWrK7eo5kSWQj5HV86/BpIkKRtMJVqwFhfyrF3eV/Waiwf66CnkW1SRJElSdYZzLViFfI71q/tZOsu4xKU93Qyu6qfbeeeSJCkjDOda0Hp7utk6uJJ1K5YdaXEp5HOsW7GMrYPnO+dckiRlSqSU2l1DZgwMDKShoaF2l6EmKJUrjE7ZIbSnkHfFXJIktU1EDKeUBqYfdxMidYTufM4wLkmSMs+2FkmSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKCMO5JEmSlBGGc0mSJCkjDOeSJElSRhjOJUmSpIwwnEuSJEkZYTiXJEmSMsJwLkmSJGWE4VySJEnKiMyG84h4S0SkGR79Nd6Xj4h3RcSdEVGMiJsiYnmr6pYkSZKOVWbD+YTbgVOmPX5c4z1XApcArwb+f3v3HyZZVR54/Psy00AbdkZUZIAJw/pMYhBXkWlUwHUZlYDosoLGrBEJbBZwGU2UFdYoJmgM8UeMRIUoohJRJNElKioEQd3wQ4SZICiCCjojA8zIIM6IdENP8+4f5/ZQlF1VXd3VXXe6v5/nuQ9T955z76nTh1tvnXrvvQcCPwGujIglM9dMSZIkafrqHpyPZeaGpmWsVeGIWAS8HvjzzPxyZn4fOAF4uFovSZIk1Vbdg/OlEbG+Wi6LiIM7lF8B7ARcPr6iCua/DrxgogoRcVJErI6I1ffdd1/PGi5JkiR1q87B+Q2UWe+XAa8B7geujojD2tTZo/rvhqb1Gxq2PU5mnpeZQ5k5tNtuu02zyZIkSdLUzWpw3uYiz8blTIDM/Fpmfi4zb87MqzPzWOBq4LTZbLMkSZI0WxbO8vE+AlzcocymNtu+DbyyzfZ7q/8uAX7WsH73hm2SJElSLc1qcJ6Zm2gffHdyAHBXm+1rKBd/Hg58HCAidgBeApw3jeNKkiRJM662OecR8XcR8aKIeFpE7B8R5wCHAWc3lDk6Im6PiL0AMnML8FHgrIh4eUTsB3wSGAQ+1oe3IUmSJE3abKe1dGMP4NPAbsBm4BbgJZn5jYYyi4GnAwMN604DHgHOB55ImU0/LDNNa5EkSVKtRWb2uw21MTQ0lKtXr+53MyRJkjTHRcSazBxqXl/btBZJkiRpvjE4lyRJkmrC4FySJEmqCYNzSZIkqSYMziVJkqSaMDiXJEmSasLgXJIkSaoJg3NJkiSpJgzOJUmSpJowOJckSZJqwuBckiRJqgmDc0mSJKkmDM4lSZKkmjA4lyRJkmrC4FySJEmqCYNzSZIkqSYMziVJkqSaMDiXJEmSasLgXJIkSaoJg3NJkiSpJgzOJUmSpJowOJckSZJqwuBckiRJqgmDc0mSJKkmDM4lSZKkmjA4lyRJkmrC4FySJEmqCYNzSZIkqSYMziVJkqSaMDiXJEmSasLgXJIkSaoJg3NJkiSpJgzOJUmSpJowOJckSZJqwuBckiRJqgmDc0mSJKkmDM4lSZKkmjA4lyRJkmrC4FySJEmqCYNzSZIkqSYMziVJkqSaMDiXJEmSasLgXJIkSaoJg3NJkiSpJgzOJUmSpJpY2O8GaPs2PDrGluFRto49ysIFO7BocIDBgQX9bpYkSdJ2yeBcU7Zh8wjnfPMOvrBmPcOjYwwOLOBVK5ayauVylizeud/NkyRJ2u4YnGtKNmwe4Zhzr+WezSPb1g2PjnHh9eu46raNXHLKwSxZPNjHFkqSJG1/zDlX14ZHxzjnm3c8LjBvdM/mEc791p2MjI7Ncsvmp+HRMTZuGeHuBx5i45YRhu13SZK2W86cq2tbhkf5wpr1bct8fvV6Vq1czs7mn88oU4skSZpbDM7Vta1jj3acnR0eHWPr2KOz1KL5ydQiSZLmHtNa1LWFC3boeEeWwYEFLFzg8JopphZJkjQ3GT2pa4sGB3jViqVty/zB0FIWDw7MUovmn8mmFm0eHp2lFkmSpF4wOFfXBgcWsGrlcvZskdO85+KdOeVQ881nkqlFkiTNTQbnmpIli3fmklMO5riDlm1LcRkcWMBxBy3jklMO8WLEGWZqkSRJc5MXhGrKliwe5G1H7suqlcu3PSF08eCAM+azYDy16MLr17UsY2qRJEnbH4NzTcvOAwsMxvtgPLXoqts2TnhRqKlFkiRtn/zNW9pOmVokSdLcE5nZ7zbUxtDQUK5evbrfzZC6MjI6xubh0VqmFg2PjrGloW2LBgc65spLkjQfRMSazBxqXm9ai7Sdq2tqkU8vlSSpewbnknrOp5dKkjQ15pxL6imfXipJ0tQZnEvqKZ9eKknS1BmcS+opn14qSdLUGZxL6imfXipJ0tT56Sipp8afXtqOTy+VJGliBueSemr86aV7trhdok8vlSSpNYNzST3n00slzUXDo2Ns3DLC3Q88xMYtIx2vr5GmwieENvAJoVJv1fnppZLUDR+spl7zCaGSZl1dn14qSd3wwWqaTaa1SJIkteCD1TTbDM4lSZJa8MFqmm0G55IkSS34YDXNNoNzSZKkFnywmmabI0mSJKkFH6ym2WZwLkmS1IIPVtNsq21wHhFrIyInWG6dQr1rZqvdkiRpbvHBappNdb7P+YFA49fQXYBbgIsnUfe9wNkNrx/pYbskSdI8s2TxIG87cl9WrVzug9U0o2obnGfmfY2vI+JEYAA4fxLVH8zMDTPSMEmSNC/5YDXNhtqmtUzgZODSzLx3EmXfEBH3R8StEfGhiHjyTDdOkiRJmq7azpw3ioghYAXw9kkU/zBwM7AR+D3g3cDhEbF/Zg5PsO+TgJMA9t577561WZIkSerWrM6cR8SZLS7ybFzOnKDqycBPgSs6HSMzP5CZV2bm9zLz88BLgd8Bjm5Rv1wmwwAAE4JJREFU/rzMHMrMod12220a706SJEmantmeOf8InS/o3NT4IiIWAa8B3p2Z2e0BM/MnEbER2KfbupIkaXYNj46xZXh020WXiwYHOj4ESJpLZjU4z8xNNAXfk3AssCPwqakcMyL2Ap4K3DWV+pIkaXZs2DzCOd+8gy+sWc/w6BiDAwt41YqlrFq53NsVat7YHi4IPRn4YmZubN4QEUdHxO1VAE5EHBQRb4mIAyJiWUQcDnwF+BnwL7PbbEmSNFkbNo9wzLnXcuH16xgeHQPKLPqF16/jmHOvZcPm37hsTJqTah2cR8TzgWcBH2tRZDHwdMotFgEeBo4BrgR+BJwLXA8clJkPzmxrJUnqreHRMTZuGeHuBx5i45aRbUHrXDM8OsY537yDezaPTLj9ns0jnPutOxmZo+9falTru7Vk5vVAtNl+AXBBw+t/Bw6e8YZJkjTD5lOKx5bhUb6wZn3bMp9fvZ5VK5d7n3HNebUOziVJmo/GUzwaZ5LHUzyuum0jl5xyMEsWD/axhb21dezRjr8KDI+OsXXs0VlqkdQ/tU5rkSRpvpmPKR4LF+zQ8Y4sgwMLWLjAsEVzn6NckqQamWyKx+bh0Vlq0cxbNDjAq1YsbVvmD4aWsnhwoG0ZaS4wOJckqUbmY4rH4MACVq1czp4tcun3XLwzpxxqvrnmB4NzSZJqZL6meCxZvDOXnHIwxx20bNv7HxxYwHEHLeOSUw6ZcxfBSq14QagkSTUynuJx4fXrWpaZqykeSxYP8rYj92XVyuXbnhC6eHDAGXPNKwbnkiTVyHiKx1W3bZzwotC5nuKx88CCOfvepMmYW7+JSZI0B5jiIc1fkZn9bkNtDA0N5erVq/vdDEmSABgZHWPz8KgpHtIcFBFrMnOoeb1pLZIk1ZQpHtL8Y1qLJEmSVBMG55IkSVJNGJxLkiRJNWFwLkmSJNWEwbkkSZJUEwbnkiRJUk0YnEuSJEk1YXAuSZIk1YTBuSRJklQTBueSJElSTRicS5IkSTVhcC5JkiTVhMG5JEmSVBMG55IkSVJNGJxLkiRJNWFwLkmSJNWEwbkkSZJUEwbnkiRJUk0YnEuSJEk1YXAuSZIk1YTBuSRJklQTBueSJElSTRicS5IkSTURmdnvNtRGRNwHrOt3O/rsKcCmfjdijrAve8e+7B37snfsy96wH3vHvuyd2ejLZZm5W/NKg3M9TkSszsyhfrdjLrAve8e+7B37snfsy96wH3vHvuydfvalaS2SJElSTRicS5IkSTVhcK5m5/W7AXOIfdk79mXv2Je9Y1/2hv3YO/Zl7/StL805lyRJkmrCmXNJkiSpJgzOJUmSpJowOJckSZJqwuB8HouItRGREyy3TqHeNbPV7jqKiDNb9OXyDvUGIuJ9EXFvRAxHxDURsWK22l1HEXFaRHw7Ih6IiF9WfXLEJOrN63EZEUdGxHcj4uGqL06dZL3TI2JdVe+miPj9mW5rnTn+emcq50XPiRObyue1Y7KIiBdGxJeq81xGxBkTlHleRFwXESPV2PubiFgwiX0fHxE/rM6ft0fEa3vRZoPz+e1AYI+G5XeAYeDiSdR9b1Pdo2aojduTtTy+T/YAftqhzvuBPwFOpvw9fgJcGRFLZq6Ztfci4JPASuC5wHXAVyLikEnUnZfjMiKGgC8BlwH7A2cCZ0XE6zvUexPwTuAdVb2vA5dGxLNmtMH15vjrrbV0d170nDixqX5eOyZhF+AHwOnAhuaNEfHblHPfD4EVwP+ijL+/brfTiHgF8Ango8CzgfOBT0fES6fd4sx0cSEzAU4ERoE9OpRbC5zR7/bWaaEEQ3d0WWcRMAKc1LBuAeXkcWa/31OdFuAW4AMdyszbcQlcBFzXtO79wNo2dQK4Gziraf2NwAX9fk91Whx/U+63rs6LnhO76tuOn9eOycn1CXAWsB7YoWHdKuDXwG+12dd1wEVN6z4PfGu67XTmXI1OBi7NzHsnUfYNEXF/RNwaER+KiCfPdOO2A0sjYn21XBYRB3covwLYCbh8fEVmjlG+wb9gBtu5XYmIHSgf2r+eRPH5Oi4PoWEcVS4HlkXE0hZ19gH2bFHP8Vdx/E1bN+dFz4mTN9nPa8dkZ4cAV2Tmow3rLgeeADxnogoRsSPl14yJzp/Pn0xKTDsG5wK2/Sy+AvjYJIp/GDgWOJQyM3I4cF1EDM5U+7YDNwAnAC8DXgPcD1wdEYe1qbNH9d/mn9k2NGwTvA14Ip0fCDGfx+UeTDyOxre1qtNYrrGe4+8xjr+p6/a86JichC4+rx2TkzOV8+dTgIUt6u0EPGk6DVo4ncqqn4g4E/jLDsXemZlnNq07mZIHeEWnY2TmBxpefi8i1gB3AEdTfl6fE7rpy8z8WtP6q6sZy9Mosz7z2lTHZUScQgmOjsrM9e0qz5dxqdnj+Jsez4szZlKf147J7ZfB+dzzETpfILKp8UVELKLMarw7q6SpbmTmTyJiI+Vn8rmk675s8m3glW22j/8cuQT4WcP63Ru2zRVTGZdvoVyseFRmXtntAefwuJzIvZRx1Gj3hm2t6lDV+1FTvbk2/rrm+Jsx7c6L8+mcOCXT+bx2TLY0lfPnJmBri3oPA7+YToMMzueYzNxE+4BxIscCOwKfmsoxI2Iv4KnAXVOpX1dT7MtGB9C+T9ZQ/ic+HPg4bMtvfQmdf0LfrnTblxHxLuDNwJGZ+f+mcsy5Oi5buJYyjt7VsO4IYF2bGd+1wD1VvX9rqjfvbrfWyPE3o9qdF+fNOXEapvx57Zhs6VrgdRGxQ0Pe+RHAQ8BNE1XIzEci4kbKWP10w6YjgOurayWmrt9Xzrr0fwFuBv65xbajgduBvarXBwFvoZxgl1UD8ybKT2y79Pu99LEP/45yC7anUW5Jdw7wKPBfW/Vlte5s4D7g5cB+wAXAA3S4Y85cXqo+GQZeQZmVGF8Wt+rL+T4uKRcmjVJu/fV7wB9Xffj6hjLPrfrsuQ3r3kT5ADq2qvceSnD07H6/J8ff9r90Oi96TpxSn074ee2YbNtnu1Tjb3/KhMRHqn8vr7b/NrCFclvE/Si3m7wfeE/DPvaq+vfohnWvoMye/xnwdODU6vVLp93mfneaS38X4PlAAi9usf34avs+1esDKLcP+kX1IX4n8A/Akn6/lz734+cot2J6GPg5cCXwonZ9Wa0bAN5HuYhkhPINfqjf76fPfZktlgta9aXjMqFcdHdz9f7XAac2bT+06rNDm9b/H0oKwcPAd4HD+/1e+tyPjr/e9WXb86LnxK77s+XntWOybb+Nn/ual2819e111ZjbAPwNsKBh+z5VneMn6PcfAY9Q7pN+bC/aHNXOJUmSJPWZt1KUJEmSasLgXJIkSaoJg3NJkiSpJgzOJUmSpJowOJckSZJqwuBckiRJqgmDc0maARGxNiLO6Hc76iYiMiKO7aL8PlWdF8xkuxqO599NUl8ZnEvqi4i4oAq6MiK2RsS6iPhoRDy5qdxBEfEvEbExIkYi4s6I+ExEHNBQ5u0RcXVEbKn2t7TDsY+uyj2jxfZzq/ZM5xx5IPDBadSfNRFxZsPfonl5Sp+bdxewB/CdPrejpyLiyIj4bkQ8XH0hOHWS9U6vxubDEXFTRPx+m7IvioixiLijaf23Wvytfz3d9yVp+gzOJfXT1ZTAax/gT4FXAp8e3xgRJ1RlHgFeC+wL/CGwFvj7hv3sBHyZ8vj6ybgUuBc4sXlDRDwB+CPgE5n5aDdvpqq/I0Bm3peZ21Ows5byt2he7u9jm8jMsczckJmj/WxHL0XEEPAl4DLKY8TPBM6KiNd3qPcm4J3AO6p6XwcujYhnTVB2CfCPwBUT7OoYHv833hO4G7h4au9IUi8ZnEvqp0eqwGt9Zn4JOBs4IiIGI2JPyuOmz8/MP8zMKzPzp5m5OjPPAI4a30lm/kVmvp9Jzq5m5lbgk8DrImKnps2vBnYBPhERB0TEZRHx84h4MCJujIgjGgtXs57vrmbb76d8mfiN9IiI+KOI+E5EbI6ITRHx1Yj43Ybt4+kbr46Ir0TEQxHxk4g4vul4u0TE2RFxV8Os69satu9e/SpxX0T8KiKujYgXTqJbxoPg5iUjYqdqlvaLDccZjIjvR8RFTe0/NiKuiojhqv3/vd1BI+LPqhnkByNiQ0RcHBF7TNAvL5hCP/19RNxdlbkpIo5pKvPsiLiu6scfR8SrJ9FPvXAqcGNm/nlm3paZFwAfBt7aqkJEBHAa8MHM/HRV73Tglmp/jWV3AD4DnMME/09k5i8a/8bAM4G9gI/25u1Jmg6Dc0l1Mkw5Ly2kBMk7Ae+eqGBmPjDNY50P7EqZRWx0IvDVzLwbWAT8E7ASOAD4V+DLjUF15U+BnwMHASe0ON74ezkAOAwYA746PtPe4D2UXw+eRZnJPH/8eFWA9hXKF5M3Un5JOA64r9o+CHwT+A/AS4HnAF8Dvh4R+3bskRYy82HKLxYviYg3VKs/BOwMnNxU/H2ULz77AxcBn42I53Q4xFuA/wQcDezN5GZwO/XTpcCzq3Y/k/JF7+KIeHFVZpDSN78Enkvpx9OAp7Y7aETsXX2RaLfc2qHthwCXN627HFgWrVOy9qHMcE9Urzkf/x1AAu/t0I5xrwduyswbJ1le0kzKTBcXF5dZX4ALgCsbXj8DuBO4vnp9LrC5y30eSglKlk6y/GXANxpe71vVf1mbOjcDb294vRa4aoJya4Ez2uznSdWxDqle71O9PrWhzALgV8DJ1esXV2WGWuzzeGA9sLBp/TeAs9u05UzgUeDBpuXmpnJ/DIwA76KkGh3YsG28/X/VVOc64MKG1wkc26Ytz6nK7NW03xd00U+HVu1c3LTvTwJfrP79P6v3uGvD9mdW+273d1sILO+wLOsw7h4BTmpat1917ANb1Dm42v67TetXAb9ueL2SkrK1pOFve0ebtuwBjI73nYuLS/+XhUhS/xwaEQ9SgqudgKt4bCY2ZuH45wH/NyKWZ+YdlFnzn1GCdiJiN0qO74uAJZTAbGdgWdN+buh0oIjYH/hLyozyU3js/S0Drm0o+t3xf2TmWET8HNi9WrUCeCAzV7c4zIFVO39ZJo+32Ynyq0Q7d1GC/0aPNL7IzH+MiCMpM7NvzYlnWr/d9PraCfa7TUQcCvw55cvZE3nsF91llDzoVtr104HAjsDdTf2wI/Dj6t/PAG7Lhl9gMvP7EbG5zTHJkhJ1R7sy/RLl4t3PACdkSVeZjP9B+SJz0Yw1TFJXDM4l9dN3KLOxW4F7MrMxGPwhsCgilmbm+hk6/qXABuDEiHgHJbXhQ/nYhaAXUNIsTgd+SglwL6YEeY3aXvgZ5SLTK4BrKGkvG6tNt06wr0eaXieTT0HcAbiNkh7S7KEOdUerLygtRcQulLScMaA5tadrEbE3JbXkQsps/CZgKXAlv9kvzdr10w7AZkqQ3qleV6o2/6BDsXWZuV+b7fdSvkQ12r1hW6s6VPV+1FRvfNszKakvX2n4UrJDaXZsBY7LzG1BeJWbfiLw2cz8VZv2SppFBueS+mm4TUD4eUpe8RmUnNjHiYhdc5p555m5NSI+SQlQvkeZuf1EQ5EXAqdn5perY/4W8DTg+10eal9gN0o6zG3Vvg6m+18H1gC7RsRQi9nz1ZQvGFsy8+dd7nsy/oGSAvES4IqI+NfM/OemMs+nBNzjDqZ1MHsgMAi8KTOHASJiRQ/auZryt9w5M1v9rX4AnBQRT8zMX1bH3g9Y3GHf91B+/Win051lrgUOp3whGXcEJahv9UV0bXXsw4F/a6p3TfXvGym5+41OAV4OHEn5daTREZRfKD7Wob2SZpHBuaRaysy7q4sPPxYRTwQ+TslJfxLw3yi5tS+EbbOZT6Lk+wI8o/qJ/2eZ+YsOhzqfklbxIR67EHTcD4HXRsQ1lNSbd1X/7dY64GHgjRHxAUre9Hsos73d+AblbjD/FOW+2LdQZkr3zczzgc8Cb6ZcaPp2ygzr7pS0nNsy84sT7xaABVFuv9dsU/Ul5nXAq4DnZeYt1f7Pi4gbMnNtQ/k/iYjbKQHysZSLZN/Y4pg/pvTB/46Iz1Iu4PyLzt3Q0Tcos++XRMT4HU12pXxRGMnMj1PSOP4K+Ez1XgYpt+dsm/7To7SWDwLXRcRfU341eB6lj948XiAinku54PW4zLwhMzMi3k+55eJtlP49ntJnJ1Zt+zVNXxyrdJ9HWnxJOZly15ibpvl+JPWQd2uRVFtVwPlfKHnen6MEy18A/iPlDinj3gXcRAngodxV5SYabrfY5hhrKSknu1Jy0BudQDlP3gB8kXJnjK7vaJGZmyiB6mGUVJa/pdyhpKv7qGdmAi+jzEx/lNIfn6HksJOZI5T+Wg18ihKcX0K5G8m6Drvfh5Ie0bzsHxHLKbflOy0zb6nK/y1wPXBRRDRO9LwVOIkSEL+OcvHnv7d4P7dQgtKTKTPZbwHe1Lkn2qv66SjKe/8gcDvwVUrf3VmVeYgym/xkyt/3s1XZmfjFobl9NwKvoMxo30wZv2/PzMZbGT4BeHr13/F6Z1OugTirqncEcFRm3txtGyJiL0p/OGsu1UyUc5gkSVMXEftQ8vL/c2Ze0760JKkVZ84lSZKkmjA4lyRJkmrCtBZJkiSpJpw5lyRJkmrC4FySJEmqCYNzSZIkqSYMziVJkqSaMDiXJEmSauL/A2AuBNPpEIZ6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "sns.scatterplot(x=pcs[0], y=pcs[1], ax=ax, s=100)\n",
    "ax.set_xlabel(f'PC1 Variance Explained = {variance_explained[0]:0.3f}')\n",
    "ax.set_ylabel(f'PC2 Variance Explained = {variance_explained[1]:0.3f}')\n",
    "#plt.savefig('pca.png', dpi=200, bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit epigenetic clock using penalized regression\n",
    "- under for loop selected sample is held out as the testing sample\n",
    "    - model is trained using all other samples\n",
    "    - age prediction is made and saved for the held out sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "84e6391765264ad5b183798cf300eeab",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=48.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "lm = ElasticNet(l1_ratio=.7, normalize=False, fit_intercept=True)\n",
    "\n",
    "predicted_ages = {}\n",
    "\n",
    "# tqdm for interactive progress bar \n",
    "for index in tqdm(range(len(sample_ages))):\n",
    "    # leave on out cross validation \n",
    "    training_samples = [list(meth_matrix)[count] for count in range(len(sample_ages)) if count != index]\n",
    "    \n",
    "    # retrieve training values and ages for fold \n",
    "    training_ages = np.array([sample_ages[sample] for sample in training_samples])\n",
    "    training_values = meth_matrix[training_samples].values\n",
    "    \n",
    "    # fit linear model\n",
    "    lm.fit(training_values.T, training_ages)\n",
    "    \n",
    "    test_sample = list(meth_matrix)[index]\n",
    "    \n",
    "    # predict test age\n",
    "    predicted_ages[test_sample] = [sample_ages[test_sample],  lm.predict(meth_matrix[[list(meth_matrix)[index]]].values.T)[0]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Evaluation \n",
    "- fit a trendline for the of known vs predicted ages\n",
    "- evaluate the fit of the expected vs predicted ages based on the trendline using Pearson correlation ($R^{2}$)\n",
    "- Plot data and trendline w/ equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc('text', usetex=False)\n",
    "\n",
    "def r2(x,y):\n",
    "    # return r squared\n",
    "    return stats.pearsonr(x,y)[0] ** 2\n",
    "\n",
    "def plot_known_predicted_ages(known_ages, predicted_ages, label=None):\n",
    "    # define optimization function\n",
    "    def func(x, a, b, c):\n",
    "        return a * np.asarray(x) + c\n",
    "    # fit trend line\n",
    "    popt, pcov = optimize.curve_fit(func, known_ages, predicted_ages)\n",
    "    # get r squared\n",
    "    rsquared = r2(predicted_ages, func(known_ages, *popt))\n",
    "    # format plot label\n",
    "    plot_label = f'$f(x)={popt[0]:.2f}x {popt[2]:.2f}, R^{{2}}={rsquared:.2f}$'\n",
    "    # initialize plt plot\n",
    "    fig, ax = plt.subplots(figsize=(12,12))\n",
    "    # plot trend line\n",
    "    ax.plot(sorted(known_ages), func(sorted([1 + x for x in known_ages]), *popt), 'r--', label=plot_label)\n",
    "    # scatter plot\n",
    "    ax.scatter(known_ages, predicted_ages, marker='o', alpha=0.8, color='k')\n",
    "    ax.set_title(label, fontsize=18)\n",
    "    ax.set_xlabel('Chronological Age', fontsize=16)\n",
    "    ax.set_ylabel('Epigenetic State', fontsize=16)\n",
    "    ax.tick_params(axis='both', which='major', labelsize=16)\n",
    "    ax.legend(fontsize=16)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/u/home/c/colinpat/.local/lib/python3.7/site-packages/scipy/optimize/minpack.py:808: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtoAAALXCAYAAACzcOftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxU5fXH8c9hCZCwSERtXGqiAgoiqNGKS5HULVbrRtXihnW3qSvRn2LRuiAarVZTbdWquIsrboNbBAUUjILKjpKwjgoSUBIgLM/vj2dCkzCBBJLcWb7v12tek3nunXtPJvNKTp459zzmnENERERERBpXi6ADEBERERFJREq0RURERESagBJtEREREZEmoERbRERERKQJKNEWEREREWkCSrRFRERERJqAEm0RkRhjZoPMzJnZkUHH0lTM7MjI9zgo4DicmT3ZxOfIjJznlqY8j4jEHiXaIhJzIklJfW+Z1RKZ6rdKM5tnZi+b2W+inGMPM3vEzGaaWYWZlZnZDDMbYWb96xFjtHPWvu3aNK9QfDCzPmZ2i5llNvN5dzGzu83sazP7xczWmFmpmT1jZr9rzlhEJLm1CjoAEZEozqn1+AjgYuAR4JNa25YAO0S+fh94KvJ1G6Bb5HknmdlhzrlJAGaWDYwF1kb2nwa0A7oCxwC/AB/VM9bq56xtWT2PUdvTwAtA5VY+P1b0AW4GxgCltbZ9jH/N1zbmCc3s98Dz+J//S/j3zCogEzgZ+MDMfu+ce6cxzysiEo0SbRGJOc65Z6o/NrNW+IT509rbIturEu3ZUZ47HhgFnAVMigzfDKQCfZxzX0U53q8aEO4m59xWzrn1wPrGPGascc5tAFY35jHNrCc+uV4GHO2cm1Fr+1D8+6BRzysiUheVjohIolscua8+O9wV+Clakg3gnPu+sYOoqgU2s6PM7LNIucr3ZvZPM2tfa9+oNdqRcpVXzOznyG2UmWVFyiLGRDnnUWb2npktN7PVkVKKS6PsV2pmY8xsbzN7O1JusSJSdrPJPx1m1snM7jKzbyNlGUvM7Hkz26PaPrcAT0QeflStnObJyPaoNdrmXWRmE81sZeT2jZndWo+X+Vb8LPmFtZNsAOc945wr2tKBzOxCM/vSzFZFXov3zOzwOvbtH3ndfoq8znPN7L9m1mUL5zg28lp/Ymad6/H9iUic0Yy2iCSSttWSmxR86cidQAW+HKPKd0B3MzvVOfdqI56zunXOueW1xg4ABgCP4stN+gNXAPua2dGRWd6ozGx7fNnMTsC/gRn4kpqPgLQo+18c2e8z4A6gHDgaeNjM9nTO5dd6yi74Eo/XgHygN3AJ0BFfTlN13E7ABODXwOP4spsM4HJgopllO+fmAa9Gxi8GhkXiBf/ab87T+FnniZG4lwN741+3oXU9yczaAr8HFjjnRm/hHJtlZncB1+E/AbkR6ID/Pj4ys5Oql52Y2SXAw8CiyP08/GtzIrArsLSOc5wHPAa8CQx0zmmWXSQROed000033WL6BgwCHDCoju2Zke3RbguAQ2vt3xc/w+2A2fiE8TJgnwbEtLlzOmBqrf2rxk+uNf7PyPiZUb7fI6uN3R0ZO6vW86vGx1Qby8CXRzwXJe5/4stS9qg2Vho5xum19v1XZLx7reevAnrX2nd34Gfgyc19H9W2HVn7ZwqcHhl7GmhRa/8WtY9Ra3uvyHPfaOB7y9WKuTuwARgHpFQb3xmf9JcCLSNjuwJrgOnAdlGO3aLWe+WWyOMbIo8f2tL3pZtuusX3TaUjIpJIRuFnbY8GjsfPFlcCb5jZ/lU7Oec+BQ4ERgCdgPPxSc90M/u4eglEA89Z/XZhlH1nOederzU2PHJ/yhbOcyIQxl/oV909UfYdgL8Y8L9m1qX6DT+D2gI4qtZzFjvnRtYaqyqx6Aq+rAM/2/wxsKjWccvxs+fHsPXOitwPdrVm92s/jqJj5P7nbTg/wEmAAXc75zaWGznnFuNLYXYHqt5Lf8R/cvJ3t+mnF9FibmFmhfgZ/r855y6vx/clInFMpSMikkgWOuc+qD5gZm8As/Af6x9SNe6c+wY/44qZ7Q70wyfHRwCjzOzA6olWQ865GdHqhsNmthzYUnKfBUyKkoD+GHl+dftE7jcX1061Hs+Nss9PkfvtI/c7RL4+Bt/tJZptSRy7AmHn3A9b8dyqBLvDNpwf/OsMviSmtqqxPYBiIv+AAJPreeyr8PENcc4N2+oIRSRuKNEWkYTmnJtnZjOB35hZmnOuPNo+wFNm9jS+Dvow4GB8+UA8ssj9ufhZ8GhqJ9ab63Jite4/AO7autCazBx8GUefoAPZjPeB3wIXm9kLzrlo/9yISAJRoi0iyaB15L49vsQhKuecM7OJ+ER7lyaIY5/aA2aWAWxH9Bnl6kqBvcysRfVZbTPbMfL86uZE7pc2YLa9Ppbg65Q71vO4roHHn43veb5TQ2e1nXOrzewd4BQzO8Y5914Dz12l6ufQk00v3OxRa5/Zkfs+1b7enG/wF3QWAWPNLMc5N2cLzxGROKYabRFJaGbWA999ZFFV8mZmR5vvzV1733b8r8Z4ehOE093MTq41dn3kvnbtdm1v4i9y/FOt8cFR9h2Jn939e+R7qiHSnq9NPeKtIZLgPwscbGYDou0TSfyrrIzcp9fzFM9G7u82sxp/nyL14VsyFH+h5mNm1r2O+AaaWc5mjvEG/h+EfDOr+get6h+i8/FdRapKRV7GXwNws5l1rH2gaDE756bhy5Ra4pPtvevxfYlInNKMtogkkm5mdnbk69bAnvgWda34X0ILcB+wfaR++xt8+7/dgIH4pPypSA13Q89Z2weuZk/ub4BnzOxR/Kxzf/yFi2OBF7dwnrsi8T1hZgcDM/H15IfiW8htnD12zi00s8vw7eNmREpi5uFrrHvhV0jswaarNdbHEPyM/0gzG4m/ALISf5Hg8cAXRGrfgc/xNdtDIn2iy4ES59zEaAd2zr1kZi/iS166Rn4+ZfifybHAvpsLzDk31cz+iL9g9KtIfBPxyffu+AsdewO5mznGLDMrwLf3+zgST1V7v/b4ri/rI/suNLOr8N1ZvjGzp/Cv8y6Rc/0ZmBLlHDPNrB9+ZnuMmf0ukoCLSIJRoi0iiaSq4wf4xHM5Ptm7xzn3frX9rsEnQocDp+FLL1YAX+MT2ie38pzRtlVPtL+MnPsO4FL8BXyFwI1b6j7hnFsaWTDlXnwC5/A9tPvjv8dVtfZ/wsxm42e8L8F/j0vxF4b+rVZc9eacW2FmhwHX4tvxnQSsAxbia9ofq7bvfDP7M/6fnIfx//yMwCe/dRmIr5O/AD9DvR4owa/4WJ/43jazffAXHh4HnBo57+JIfFc558Zs4RjXm9m3+N7gw/H/SEzE97v+pNa+D5vZd/je41fgu70sBj7Et5as6xxzqiXbH5nZUc65r+vzPYpI/DDnGlpCJyIiDWVmDhjhnBvUyMfdHp9A/8c5t8mqjyIiEhzVaIuIxIlo9dbA/0Xu34+yTUREAqTSERGR+PGOmc3Dl6C0AH4HnIBfEn1LF1OKiEgzU6ItIhI/3sJfKHgK0A5fF30vfmXCzfXBFhGRAKhGW0RERESkCSTsjHaXLl1cZmZm0GGIiIiISAL74osvljrndoi2LWET7czMTIqLi4MOQ0REREQSWOTamajUdUREREREpAko0RYRERERaQJKtEVEREREmoASbRERERGRJqBEW0RERESkCSjRFhERERFpAgnb3q8+VqxYwdKlS6msrAw6FBGJcykpKXTp0oVOnToFHYqIiMSIpE20V69ezQ8//MCuu+5Ku3btMLOgQxKROOWcY9WqVSxcuJA2bdrQtm3boEMSEZEYkLSlI0uWLGGHHXYgNTVVSbaIbBMzIzU1lS5durBkyZKgwxERkRiRtIn26tWrad++fdBhiEgC6dChA6tXrw46DBERiRFJm2ivW7eOVq2StnJGRJpAq1atWLduXdBhiIhIjEjaRBtQyYiINCr9ThERkeqSOtEWEREREWkqSrRFRERERJqAEm0RERERkSagRFtEtlpZWRknnHAC3bp1o3fv3hxzzDF8++23QYclIiISE5Roi8hWMzOuuuoqZs+ezVdffcUJJ5zAhRdeGHRYIiIiMUGJdgJatGgRp556KjvssAMtWrTg6quv5oorruCEE05o8LHuv/9+evXqxYYNG5og0poWLFjAgAED6NSpEx07duTUU09l/vz5W3zewoUL+etf/0rfvn03LkBUWlraqOeoj+OOOw4z46abbqox/vLLL3Paaaex++67065dO7p3784NN9zAL7/8ssVjjhkzBjPb5LbddtvV2K8hr0F9PPbYYzXOl5qaSu/evXn++edr7Lfddttx1FFHbXx86KGHbtN566Mxf4ZN8TMTERGpokQ7AQ0aNIhZs2YxYsQIJkyYwIABA/j3v//NLbfc0uBjXXLJJSxZsoQRI0Y0fqDVVFRUkJOTw8yZMxkxYgRPP/00c+bMoX///pSXl2/2ud9++y0jR46kc+fOHHHEEU1yji15/vnn+eqrr6Juu+eee2jZsiXDhg1j9OjRXHbZZTz88MMcffTR9f4H5oEHHuDTTz/dePvggw9qbK/va1BfkydPpk2bNhvP9+KLL9KiRQvOOussPv744zqfd//993PSSSdt8/nr0pg/w6b+mYmIiOCcS8jbgQce6DZn+vTpm90er8LhsDMzN2LEiI1jeXl5Ljs7e6uPmZ+f73r06NEY4dXp/vvvdy1atHBz5szZODZ37lzXsmVLd++99272uevXr9/49aOPPuoAV1JS0qjn2Jxly5a5nXbayT333HMOcEOGDKmx/ccff9zkOSNGjHCA+/DDDzd77I8++sgB7v3339/sfvV9Derr0EMPdb17964x9vnnnzvAXX/99VGfc8stt7i+ffu68vLyrT7vljTWz7Apf2aJ+rtFRESiA4pdHfmoZrQTyGmnnUZGRgbOOc477zzMjBtuuIFnnnmGgQMH1tj322+/pXXr1gwdOrTG+GWXXUaHDh0oLi7eOHbmmWcyffp0JkyY0GSxv/HGGxxyyCHstddeG8eysrI47LDDGDVq1Gaf26JF/d7G9T1HQ14bgOuvv559992XP/3pT1HPu8MOO2wydtBBBwG+zKcx1Pc1qA/nHF9//TU9e/asMb7TTjsBRF1R9fbbb+edd94hFAqRmpraaLHUti3vk+pi4WcmIiKJT4l2Arn11lu58MILad++/caP/I877jiWL1++STnBXnvtxYUXXsj999/PTz/9tPH5jz/+OK+99hrZ2dkb9+3Tpw8dOnRg9OjRUc/rnGPdunVbvK1fv77O2KdNm8a+++67yXjPnj2ZPn361rwcW32Ohrw248aN46mnnuJf//pXg2IZO3YsAPvss0+99j/rrLNo2bIl22+/PQMHDmy0uvJo5syZw8qVK+nRo0eN8ap68ZNPPrnG+N///nfefPNN3nvvPTp16rTZY2/re6Ux3ifN9TMTERHZdGoq2R155KZjp58Ol18OFRVw/PGbbh80yN+WLoUBAzbdftllcMYZsGABnHPOptuvvRZOPBFmzYLu3bc69J49e7J8+XL2228/DjnkEADuuusuzIz99ttvk/2HDh3KU089xfDhw+nevTt///vfef7552tc3AZ+trR379589tlnUc87duxY+vfvv8X4+vXrx5gxY6JuW7ZsGZ07d95kPD09nbKysi0euz4aco76vDaVlZVccsklDB48mO4N+LktWrSIoUOHctRRR9VI2qPp1KkT1157Lf369aNjx45MnjyZYcOG0bdvXyZPnsyOO+5Y7/PW15QpUwDYe++9WbduHeXl5bz//vvceOONPPjggzVinjZtGrfccgt77rkn/fr1A/yMd+1Z/yrb+l7Z1vdJc/zMREREqijRTjBTpkzh2GOP3fh48eLFdOzYkZSUlE32zcjI4KqrruLee+9l3bp1PPDAA5x++ulRj7vDDjswe/bsqNsOPPBAPv/88y3G1qFDh3p+F8Grz2tz9913s2rVKoYMGVLv465cuZKTTjqJVq1a8cQTT2xx//3335/9999/4+N+/frx29/+loMPPpgHHniA22+/vf7fVD1VJdoDav3TWFBQwF/+8pcaYz179sSXp9VP0O+V5viZiYiIVFGiXVsdM64ApKZufnuXLpvfvttum9++DbPZAL/88gvfffddjcRs9erVtGnTps7ndO3alTVr1nD44YdvkkRV165dO1atWhV1W/v27enTp88W4zOzOrd17tw56oxkXTOYW6Oh59jcazN//nzuuOMOHnvsMdasWcOaNWs2bluzZg3Lly+nQ4cOtGzZcuP4qlWrOPHEE5k7dy5jx45l11133arv44ADDqBbt271Sli3xuTJk9l+++0ZPXo0zjlKS0sZPHgwQ4YMYeDAgey8885bfextfa9sy/skyJ+ZiIgkJ9VoJ5ApU6bgnKuRaG+//fYsX7486v4ffvghl1xyCX379mX8+PF8/fXXdR572bJldOnSJeq2sWPH0rp16y3efve739V5/J49ezJt2rRNxqdPn75JrfDWasg5tvTazJ07l9WrV3P22WfTuXPnjTfwreE6d+7MN998s3H/tWvXMmDAAIqLi3nnnXfo1avXNn8/m/vHZVtMmTKF7OxssrOzOeigg/jjH//IQw89RGVl5SZ9tBtqW98r2/I+iYWfmYiIJBfNaCeQyZMn07p16xoXi+29995UVlaycOHCGrNxX375JaeccgoXXngh9913H926deOGG27g7bffjnrskpISDj744KjbGqMc4A9/+AODBw9m7ty57LHHHgCUlpYyfvx4hg8fvsVj10d9z1Gf16ZPnz589NFHm5yjf//+nH322VxwwQUbO2Ns2LCBs846i6KiIt56662N9fNbq7i4mFmzZm1S2tEYfvjhB77//nv+/Oc/1xjPzc1lxx135LXXXuPaa6/d6uNv63tlW94nQf7MREQkSdXV9y/eb8nYR3vQoEGb9D4uKSlxgHvllVc2js2ZM8ftuOOO7rTTTtvYf/nxxx93gBs7duwmxy0rK3Nm5h599NEmi33lypVuzz33dPvuu697/fXX3ahRo9x+++3nsrKy3C+//LJxvzFjxriWLVvW6BPunHMvvfSSe+mll9yll17qAPfQQw+5l156yY0ZM6ZB52joa1MbUXoyV8U0ZMgQ9+mnn9a4LViwYLPf28CBA92QIUPcK6+84j788EN3zz33uO23397ttttubsmSJQ1+DZz733vi5ptv3iT+UCjkAPfqq69usu2CCy5wLVq02OS8zWlb3yfRbMvPLJpE/N0iIhLTfvzRufHjAzs9m+mjHXhC3FS3ZEy0e/fu7c4777xNxg8++GA3aNAg55xf0CYrK8v169fPrV69euM+69atc3vvvbfr27fvJs9/5plnXJs2bdzSpUubLHbnnJs3b5479dRTXYcOHVz79u3dSSedtMmiK1ULuDzxxBM1xoGot379+tX7HFvz2tQWLWnbfffd64yverIb7XsbNmyY69Wrl+vYsaNr1aqV23XXXd1FF13kFi9eHPXc9XkNpk6d6gD38MMPb3KMO++80wFu/vz5m2x78803o772zW1b3ifRbMvPLJpE/N0iIhKTKiudu+8+5zp1cm633fzjAGwu0Ta/PfFkZ2e7ulqMAcyYMSNp+uE++eSTXHnllYTD4a1aTCQ3N5cuXbrw9NNPN0F00tweeeQRhgwZwrx585p0cZlklUy/W0REAlNUBHl5MGMGHHss3HcfBPS718y+cM5F7f2qiyGTwNlnn83OO+/MQw891ODnTpkyhaKiIm6++eYmiEyCMHbsWK6++mol2SIiEr++/BIqK+GNNyAUCizJ3hJdDJkEqvr/fvnllw1+7vfff8+TTz5ZY8lriW/PPvts0CGIiIg0zMqVMGwY9O7tFwG84gr4619hMy2MY4ES7SRxyCGHbFXnhOOOO64JohERERGpB+fg2Wfh+uth8WLIz/eJdpSF+GKRSkdEREREJPZMngyHHQbnnAO77AKffgp33x10VA2iGW0RERERiT3ffQdz58Ljj8N550GL+JsfVqItIiIiIsGrrIQHH4SWLeGqq+C00+C446B9+6Aj22rx96+BiIiIiCSWUAj22w8GD4YJE3xttllcJ9mQ5Il2ovYQF5Fg6HeKiEgDffcdnHACHH88bNgAb78NI0f6JDsBJG2i3bp1a1atWhV0GCKSQFatWkXr1q2DDkNEJH6UlcG4cVBQAFOn+oQ7gSRtjfaOO+7IokWL2GWXXWjXrh2WIP85iUjzc86xatUqFi1axE477RR0OCIisWvDBnjqKZg5E4YPh+xsWLAAOnQIOrImkbSJdseOHQFYvHgxa9euDTgaEYl3rVu3Zqeddtr4u0VERGqZONEvNDNpEvTtC2vW+AVnEjTJhiROtMEn2/qjKCIiItKEfvwRrrsORoyAX/3K3599dly262uoxP8ORURERCQ4lZXwxht+dcfZs+Hcc5MiyYYkn9EWERERkUbmnO8e8vrr8OijsOuuMG9eQpeI1CU5/p0QERERkaY3c6bvHHLiiTB+PCxZ4seTMMkGJdoiIiIisq1WroRrr4VevfyCM//4B3z9Ney4Y9CRBUqlIyIiIiKy7V56CQYNgjvuSPoEu4pmtEVERESk4SZMgDPO8Bc7tm8P06b5mmwl2Rsp0RYRERGR+lu0yLfnO+wwX4f97bd+PEnrsDdHibaIiIiIbNnatTBsGHTvDi+/DEOG+Isfe/QIOrKYpRptEREREdmyFi3glVfgmGPgnntgjz2CjijmaUZbRERERKKbPh0GDICffoKWLWHMGHj1VSXZ9aREW0RERERqKiuDK6+E/faDDz+Eb77x46rDbhAl2iIiIiLiOQePPALdusGDD8JFF8GcOXDkkUFHFpdUoy0iIiIinhm8+y7ssw888AD06RN0RHFNM9oiIiIiyWzBAt+ub+ZM/3jECBg7Vkl2I1CiLSIiIpKMVq2C226Dvff23US+/NKPt2/vZ7ZlmynRFhGJUaFQiJycHLKyssjJySEUCgUdkogkitdf9/2vhw6F44+HGTNg4MCgo0o4SrRFRGJQKBQiLy+PcDhMeno64XCYvLw8Jdsi0jjGjfMdRIqK4KWXIDMz6IgSkjnngo6hSWRnZ7vi4uKgwxAR2So5OTmEw2HS0tI2jpWXl5ORkUFRUVGAkYlIXFq2DG6+GU4+GX73O1820ro1tFJfjG1lZl8457KjbdOMtohIDCopKSE1NbXGWGpqKqWlpcEEJCLxaf16ePhh367voYegahKyXTsl2c1AibaISAzKysqioqKixlhFRQWZ+nhXROrrk0/ggAPg8suhVy+YPBmuvz7oqBpdLF/P0uyJtpkdaWYuym15rf06m9ljZrbUzMrN7AMz69Xc8YqIBCE/P5/KykrKy8txzlFeXk5lZSX5+flBhyYi8eLrr2HFCl+DXVTkV3lMMLF+PUuz12ib2ZHAR8AVwOfVNq1zzhVH9jHgEyATyAfKgBuAnkAf59zCLZ1HNdoiEu9CoRAFBQWUlpaSmZlJfn4+ubm5QYclIrGqogLuvhv22APOPRfWrYO1a32ZSIKKhetZNlejHWRxzgzn3Gd1bPsDcBiQ45z7CMDMPgVKgOvwSbqISELLzc1VYi0iW+YcvPwyDB4M8+f7UpFzz/U12Aleh11SUkJ6enqNsVi6niVWa7T/ACyuSrIBnHMrgDeBkwKLSkRERCSWTJ0K/fvD6adD584wZgz8619BR9VsYv16liAT7WfNbL2Z/WRmz5nZr6tt6wlMjfKcacCvzax984QoIiIiEsPmz/fJ9sMPwxdfQL9+QUfUrGL9epYgEu0VwL3AhUAOcBtwFPCpme0Y2ScdX5dd27LIfedoBzazi82s2MyKlyxZ0rhRi4iIiARt3TooLIS77vKPjz8eSkrg0kuhZctgYwtAbm4uhYWFZGRkUFZWRkZGBoWFhTFTdhcTC9aY2QHAJGC4c+4mM5sNfOmcO7PWfhcCjwK/ds4t2NwxdTGkiIhI86m6eLekpISsrCxdvNsUiorgyiv9DPbvfw9vvglmQUeV9GJ+wRrn3JfAbOCgyFAZ0Wet06ttFxERkRgQ6y3W4t78+XDaaX5Fx5Ur4dVXlWTHiZhItKupml6fhq/Trq0HMN85t7L5QhIREZHNKSgoICUlhbS0NMyMtLQ0UlJSKCgoCDq0xLBiBXzwAdx+O8yYAaecoiQ7TsREom1m2UB3fPkIwBvALmbWr9o+HYETI9tEREQkRpSUlJCamlpjLJZarMUd5+D55+Gaa/zjXr1g4UIYMgTatg02NmmQZm+uaGbP4vthfwksB/bHL0azCHggstsbwKfAM2ZWfcEaA+5u7phFRESkbllZWZssGhJLLdbiyuTJcMUVMG6cXz69ogJSU6FDh6Ajk60QxIz2VHyf7CeAd4GrgFeB3zjnlgI45zYAJwDvAw8BrwHrgf5bughSREREmlest1iLC8uWwSWXwIEHwsyZ8MgjMGmST7IlbsVE15GmoK4jIiIizaeq60hpaSmZmZnqOtJQP/4IPXrAOefAzTfDdtsFHZHU0+a6jijRFhEREQnC++/DM8/AE09Aixbwyy8qEYlDMd/eT0RERCRpfPcdnHwyHHMMjB8Pixf7cSXZCUeJtoiIiEhzqKiAG2/0JSIffAB33gnTpsGuuwYdmTSRZu86IiIiIpKUWrSAF1+EM86A4cNh552DjkiamGa0RURERJpKcTGcfjqsWuV7YE+ZAk89pSQ7SSjRFhEREWlsP/wAF1wABx8MH3/sW/aB6rCTjBJtERERkcayfj3cey906wZPPw3XXguzZ8P++wcdmQRAibaIiCScUChETk4OWVlZ5OTkEAqFgg5JkkWLFvDaa3D44TB1KhQUQMeOQUclAVGiLSIiCSUUCpGXl0c4HCY9PZ1wOExeXp6SbWk6c+b4OuzvvwczCIXg7bf9rLYkNSXaIiKSUAoKCkhJSSEtLQ0zIy0tjZSUFAoKCoIOTRLNzz/DdddBz54werS/0BFUhy0bKdEWEZGEUlJSQmpqao2x1NRUSktLgwlIEtOIEcKGCJQAACAASURBVNC9uy8NOftsX4d93HFBRyUxRom2xC3VYIpINFlZWVRUVNQYq6ioIDMzM5iAJDF98AHsvjtMnAiPPw6/+lXQEUkMUqItcUk1mCJSl/z8fCorKykvL8c5R3l5OZWVleTn5wcdmsSzcBjOPx++/to/fvhhmDDBt+8TqYMSbYlLqsEUkbrk5uZSWFhIRkYGZWVlZGRkUFhYSG5ubtChSTyqrPTlId26wXPPweef+/H27X2HEZHN0BLsEpdKSkpIT0+vMaYaTBGpkpubq8Ratt3o0XDFFb6ryIknwj/+AXvtFXRUEkeUaEtcysrKIhwOk5aWtnFMNZgiItKoJkzws9ajR8OxxwYdjcQhfeYhcUk1mCIi0uhWrIDBg30PbIAbb/Q12UqyZSsp0Za4pBpMEZEtU3emetqwwXcO6dbNl4dU1WG3bQspKcHGJnHNnHNBx9AksrOzXXFxcdBhiIiIBKKqO1NKSgqpqalUVFRQWVmpSYnaJk6EvDwoLoa+feGBByA7O+ioJI6Y2RfOuahvGs1oS500EyIiEr/ioTtTTPydmTYNFi+GZ56B8eOVZEuj0oy2RKWZEBGR+JaVlUV6ejpmtnHMOUdZWRlz584NMDIvsL8zq1fDffdBly5w0UW+bKSiwrfrE9kKmtGWBouHmRAREalbrK+Q2ex/Z5yDUaOgZ09/keOECX68RQsl2dJklGhLVCUlJaSmptYYU59qEZH4EevdmZr178ysWb5zyMknQ5s28N578MQTjX8ekVqUaEtUsT4TIiIimxfr3Zma9e/M4sW+k8j998NXX8HRRzf+OUSiUKItUcX6TIiIiGxZbm4uRUVFzJ07l6KiophJsqGJ/86sXw+PPgq33eYf9+8P8+bBlVdC69bbfnyJKTFxUW0dlGhLVLE+EyIijSeW/0hJ4mqyvzPjxsFBB8HFF8OYMT7pBujYcZtjlthTdVFtOBwmPT2dcDhMXl5ezPweU9cREZEkpg5DkjAWL/arOj7/POy6KxQUwBlnQLWuK5J4cnJyCIfDpKWlbRwrLy8nIyODoqKiZolBXUdERCQqdRiShLFypV86/W9/g5kz4cwzlWQngVhv3tAq6ABERCQ4JSUlpKen1xiLpT9SInVyDl57DT76CB580C+fvmCBSkSSTFZW1iYz2rHUvEEz2iIiSUwdhiQuTZ0KRx0Fp53m67BXrPDjSrKTTqw3b1CiLSKSxGL9j5RIDStWwF//Cn36wOTJfiZ78mTo1CnoyCQgsd68QRdDiogkuVAoREFBAaWlpWRmZpKfnx8zf6REali2DPbZx89k33YbbL990BGJbPZiSCXaIiIiErvGjoXHHoMnn4SWLeHnn1UiIjFFXUdEREQkvsybB6efDkceCR9/7B+DkmyJK0q0RUREJHasWQN//zvsvTe89RbccgvMmAF77BF0ZCINpvZ+IiIiEjvM4IUX4A9/8IvO/PrXQUckstU0oy0iIiLB+vprXybyyy+QkgITJ8KLLyrJlrinRFtERESC8dNPcPnlsP/+UFQE06b5cdVhS4JQoi0iIiLNa8MGKCyErl3hkUd8sj17NhxySNCRiTQq1WiLiIhI8zKD11+HAw6A+++HffcNOiKRJqEZbREREWl6paUwcCDMn+8T7ddeg/ffV5ItCU2JtoiIiDSd8nIYOtSv6DhqFHz5pR/v0MEn3CIJTIm2iIiINI2RI30/7Ntug1NPhVmz4OSTg45KpNmoRltERESaxgcfwI47+r7Yhx0WdDQizU4z2iIiItI4liyBSy/1fbAB7rsPJk1Ski1JS4m2iIiIbJu1a+GBB6BbN/jvf+Hzz/14Whq0bBlsbCIBUqItIiIiW6+oCPr0gSuvhIMP9qs85uUFHZVITFCNtoiIiGy9iRNh9WrfF/sPf1AnEZFqNKMtIiIi9bdyJQwZAq+84h9fc41fOv2kk5Rki9SiRFtERES2zDl49lno3h2GDfMXOQK0aQNt2wYbm0iMUqItIiIimzd5Mhx+OJx9NmRkwIQJcNddQUclEvNUoy0iIiKbN2sWfPut7ygyaBC00DydSH0o0RYREZGaKiuhsBBSUnwHkTPOgN//3i+bLiL1pn9JRURE5H9Gj4b99oNrr4WPP/ZjZkqyRbaCEm0RERGBuXPhxBMhNxc2bIC33oKRI4OOSiSuqXRERERE4Mcf/Qz23Xf7xWdSUoKOSCTuKdEWERFJRhs2wNNP+4scb7sNDjkEFiyAjh2DjkwkYah0REREJNlMmgSHHuo7iHz4ob/4EZRkizQyJdoiIiLJ4scf4fzz4Te/gXnz4MknYdw4lYmINBEl2iIiIsli1Sp47TW47jrfG/u889QTW6QJqUZbREQkkb39Nrz5Jjz8MOy+O8yfrxIRkWYS+L+xZjbazJyZ3V5rvE9k20oz+9nM3jCzvYKKU0REJK7MmgXHHw8nnABjxsBPP/lxJdkizSbQRNvM/gT0jjLeFfgE6AScBZwPZAIfm9mOzRmjiIhIXFm5EgYPhn33hfHj4d574euvoUuXoCMTSTqBJdpm1hm4D7gmyubrgfVArnNulHPuFeB4IB0Y3HxRioiIxBnn4Pnnff317NlwzTW62FEkIEHOaN8FTHXOPR9l2yHAp8655VUDzrmFwFTglGaKT0REJD58+ikMHAhr1/ql0mfMgMceg512CjoykaQWSKJtZocD5wJ/qWOX9UBllPE1wJ5m1rapYhMREYkbixfDOef4nthjx/rFZ0B12CIxotkTbTNLAf4D3OOcm1XHbrOAA82sdbXndQB6AgZ0ruPYF5tZsZkVL1mypJEjFxGR6kKhEDk5OWRlZZGTk0MoFAo6pOSxdi0MHw7dusHIkXDjjf7ix332CToyEakmiBnt64B2wB2b2ecBYBfg32a2i5ntDjwBtI9s3xDtSc65R5xz2c657B122KExYxYRkWpCoRB5eXmEw2HS09MJh8Pk5eUp2W4uLVrACy/AUUfB9Olwxx3Qvv2WnycizapZE20z+zUwBPgb0MbMtjOz7SKbqx63dM6Nw5eVDAAWAqX4DiQj8CUly5ozbhERqamgoICUlBTS0tIwM9LS0khJSaGgoCDo0BLXjBlw5plQVgYtW8LHH8Prr8OeewYdmYjUoblntPcA2gLPAGXVbuC7iZQBvQCccw8BOwL7Ar92zh0N7AxMdM6tbea4RUSkmpKSElJTU2uMpaamUlpaGkxAiWz5crj6athvPxg92rfqA9Vhi8SB5l4ZcgrQP8r4R/jk+7/At1WDzrk1wDQAM+sFHIW/iFJERAKUlZVFOBwmLS1t41hFRQWZmZnBBZVonPOdQ4YMgaVL4aKL4PbbQaWRInGjWRPtSLu+MbXHzQxgnnNuTOTxrsBlwAR8p5Fs4Abg1TraAYqISDPKz88nLy8P8DPZFRUVVFZWkp+fH3BkCcTML53evTu8+y7sv3/QEYlIAwW+BHsd1gK/wddkv42fxb4Vv0qkiIgELDc3l8LCQjIyMigrKyMjI4PCwkJyc3ODDi2+LVwI5577vzZ9zz7ra7GVZIvEpeYuHYnKOWe1Hv+ALxMREZEYlZubq8S6saxe7ZdKHzYM1q+H446Dvfbyi8+ISNyK1RltERGR5DBqFPToATfd5BPsGTP8Ko8iEvdiYkZbREQkaRUVQVoafPgh5OQEHY2INCIl2iIiIs2prAxuuQVOPRX69YM774SUFGilP8kiiUalIyIiIs1h/Xr4z3+ga1coLIRJk/x4aqqSbJEEpURbRESkqY0bB9nZcOml0LMnfPklqBWiSMLTv9AiIiJNrbgYfvoJXnwR/vhH3yNbRBKeZrRFREQa26pVcOut8Nxz/vFf/gIzZ8LppyvJFkkiSrRFREQai3Pw8suwzz5w880wYYIfb93a12KLSFJRoi0iIrINQqEQOTk5HLvLLkxOT/elIZ06wUcf+YseRSRpKdEWERHZSqFQiLy8PMLhMD3btiXr55/5W5cujL7jDjjyyKDDE5GAmXMu6BiaRHZ2tisuLg46DBERSVTr1vHPHj1YtWIFL++2GzhH2oYN/Lh6NRkZGRQVFQUdoYg0AzP7wjmXHW2bZrRFpF6qPh7PysoiJyeHUCgUdEgiwfnoIzjgAK6cM4fDVq/2tdlmlLdsSWpqKqWlpUFHKCIxQIm2iGxR9Y/H09PTCYfD5OXlKdmW5DN/PgwY4JdK/+UXbt53Xy7OyKjRSaSiooLMzMzgYhSRmKFEW0S2qKCggJSUFNLS0jAz0tLSSElJoaCgIOjQRJrXTz/Bu+/CbbfB9OkccvfdVK5dS3l5Oc45ysvLqaysJF+L0YgISrRFpB5KSkpIrdWaTB+PS1JwDl54Aa6/3j/ef39YuBBuugnatSM3N5fCwkIyMjIoKysjIyODwsJCcnNzg417G6lUTKRx6GJIEdminJwcwuEwaWlpG8fKy8t1wZcktsmT4cor4ZNP4IAD/DLq7doFHVWTqyoVS0lJITU1lYqKCiorKxPiHwiRpqCLIUVkm+Tn51NZWamPxyU5LFsGl14KBx4IM2bAf/4DkyYlRZINKhUTaUxKtEVkixL143GRqNasgZEj4YorYPZsuPhiaNky6KiajUrFRBpPq6ADEJH4kJubq8RaEtcHH/ha7EcfhYwMKCnxqzsmoaysrE1KxdRJRWTraEZbRESS19y5cMopcPTRvjd2OOzHkzTJBpWKiTQmJdoiIpJ8KipgyBDo0QPefx+GDYNp02DnnYOOLHAqFRNpPOo6IiIiyaeiwifZRxwBw4fDLrsEHZGIxCl1HREREfniC/jTn2D1akhNha++gqefVpItIk1GibaIiCS2H3+ECy+Egw6CoiKYNcuPJ3Edtog0DyXaIiKSmNatg/vug65dYcQIuOYa366vd++gIxORJKH2fiIikphatPAt+w491Cfce+8ddEQikmQ0oy0iCSEUCpGTk0NWVhY5OTmEQqGgQ5IgfPstDBwIS5b4RPu99+Cdd5Rki0gglGiLSNwLhULk5eURDodJT08nHA6Tl5enZDuZ/PIL/N//Qc+e8Oab8OWXfrxTJzALNjYRSVpKtEUk7hUUFJCSkkJaWhpmRlpaGikpKRQUFAQdWrNLypn9Z56B7t3hrrt8V5HZs+HYY4OOSkREibaIxL+SkhJSU1NrjKWmplJaWhpMQAFJ2pn9t96C3XaDzz6DJ5/0S6iLiMQAJdoiEveysrKoqKioMVZRUUFmZmYwAQUkaWb2v/8eLrjAr+QI8Oij8Omn8JvfBBuXiEgtSrRFJO7l5+dTWVlJeXk5zjnKy8uprKwkPz8/6NCaVcLP7FdWwj33QLdufqGZiRP9eIcO/sJHEZEYo99MIhL3cnNzKSwsJCMjg7KyMjIyMigsLCQ3Nzfo0JpVQs/sjx4NvXpBfj789rcwdSr8+c9BRyUislnqoy0iCSE3NzfpEuva8vPzycvLA/xMdkVFReLM7I8d6+/feQeS/OcsIvFDM9oiIjGqoR1EEmpm/+ef/ez1u+/6x0OHwjffKMkWkbhizrmgY2gS2dnZrri4OOgwRES2SlUHkZSUlBqz03GbONfXhg3w1FO+J/YPP8Att8DNNwcdlYhInczsC+dcdrRtmtEWEYlBSdNBpLpJk6BvXzj/fMjK8o+VZItIHFONtohIDCopKSE9Pb3GWEJ1EIlmyhRYsMB3FBk4UJ1ERCTu6beYiEgMSugOIlXWrPGrOT7xhH98wQUwaxacfbaSbBFJCPpNJiISgxK6N7hz8Oab0LOnr8X+5BM/3rKl74ktIpIglGiLiMSghOogUt2sWb5zyB/+ACkpvqvI448HHZWISJNQop0kGtomTCRIer96ubm5FBUVMXfuXIqKiuI/yQaYPx8++wzuuw+++gqOOSboiEREmoza+yWBpG0TJnFJ79cEs369r8FesgRuuMGPrVgBnToFG5eISCNRe78kl5RtwiRu6f2aQMaPh4MPhosugg8+8D2yQUl2HNCnSiKNQ4l2EigpKSE1NbXGWMK3CZO4pfdrAli8GM46Cw4/3C8689xzPtFWJ5G4UPWpUjgcJj09nXA4TF5enpJtka2g33pJICnahEnC0Ps1AaxYAW+8ATfd5C9+/NOfwCzoqKSe9KmSSONRop0EErpNmCQcvV/jkHPw2mtwzTX+8T77wMKFcNttkJYWbGzSYPpUSaTxKNFOAgnbJkwSkt6vcWb6dN855NRT4b334Jdf/LjqsOOWPlUSaTzqOiIiIg23YgUMHQr/+pdfZObWW+Gyy6BVq6Ajk22kzj8iDaOuIyIi0rjWroVnn/UdRebMgb/+VUl2gtCnSiKNRzPaIiJSPx9/7Hti//e/voOI+mGLiGhGW0REtsGCBXDmmdCvH3z4Icyb58eVZEs16r0tsikl2iIiEt2aNb72unt3GDUKbr4ZZs6ErKygI5MYo97bItEp0RYRkbo9/TSccIJPsG+5BWq1fRMB9d4WqYsSbRER+Z9vvoGBA6G8HNq0geJiGDkSdt896Mgkhqn3tkh0SrRFRASWLYO8POjTB959F6ZN8+Oqw5Z6UO9tkeiUaIuIJLMNG+Chh6BrV3j4Ybj8ct+u7+CDg45M4ohWdBWJTom2iEgyM/OlIb17w5Qp8OCDkJ4edFQSZ9R7WyQ69dEWEUk28+bBTTfB8OGwyy6+H3bHjj7pFhGRBlEfbRERgYoK36Jv773hlVf8hY7g67CVZIuINDol2iIiyeCll3yCfeutcNJJvl3fSScFHZWISEILPNE2s9Fm5szs9lrjPc3sVTNbbGblZjbNzAabWaugYhURiVvvvONrr8eOhRdegF//OuiIREQSXqCJtpn9CegdZXxnYAywB3AVcCLwOnA3cEczhigiEp+WLoXLLoMvvvCPH3jAf/3b3wYbl0hAtES8BCGwRNvMOgP3AddE2XwC0AU43Tk30jlX5JwbAowEzm3GMEVE4su6db5zSNeu8Oij8NlnfrxDB2jZMtjYRAKiJeIlKEHOaN8FTHXOPR9lW0rk/uda48uJgXIXEWlamnnaSh995BecueIKyM6Gr76Cv/wl6KhEAqcl4iUogSStZnY4fma6rr8ALwFLgUIzyzKzjmZ2CnAOcG8zhSkiAdjamScl58Ann/jOIq+9Bu+9Bz17Bh2RSEzQEvESlGZPtM0sBfgPcI9zbla0fZxzPwB9gX2AucAK4BXgLufc3Zs59sVmVmxmxUuWLGn84EWkyW3NzFPSfixcXu77Yb/xhn983XUwfTqcfLLa9YlUoyXiJShBzGhfB7RjMxc1mtkOwKtAOTAA6A/cDtxkZtfX9Tzn3CPOuWznXPYOO+zQuFGLSLPYmpmnpPtY2Dl47jno3h3uuAMmTPDjbdv6m4jUoCXiJSjNmmib2a+BIcDfgDZmtp2ZbRfZXPW4JT4ZzwSOdc694pwb45wbChQAt5lZl+aMW0Saz9bMPCXVx8KTJ8MRR8BZZ8GvfgXjx/sVHkWkTloiXoLS3D2p9wDaAs9E2TY4ctsf6AV865wrq7XPJKA1sBe+hltEEkx+fj55eXmAT5YrKiq2OPOUlZVFOBwmLS1t41jCfiw8dSrMng2PPQaDBqmTiEg95ebmKrGWZtfcpSNT8GUgtW/gk+/+wLfA98BekRaA1f0mcr+o6UMVkSBszcxTQn8svHYt3Hcf/Pvf/vFZZ8G338IFFyjJFhGJceacCzoGzMwBdzjnboo8PgT4BJ+YFwA/AUcC1wNvOedO3dIxs7OzXXFxcZPFLCKxJRQKUVBQQGlpKZmZmeTn58f/7NW778JVV/nl0k8/HV58MeiIRESkFjP7wjmXHW1bTC5n7pz7zMyOAIYC/wQ6AqXArai9n4hEkVAfC8+dC1df7buJ7Lmnvz/hhKCjEhGRBoqJRNs5t0kfKufcZ8DxAYQjIhKsRYv84jPDh/sZ7TZtgo5IRES2Qkwk2iIiSW3DBnj2WSgpgaFDfVeR+fNhu+22/FwREYlZWs5cRCRIn38Ohx0G554Lo0f7ix9BSbaISAJQoi0iEoQff4Q//xkOPtjPZD/xBIwbB61bBx2ZiIg0EiXaIiJBWLkSXn4Z8vN9X+xBg6CFfiWLiCQS1WiLiDSXd96BUAgefBD22EN12CIiCU7TJyIiTW32bPj97/3tvffgp5/8uJJsEZGEpkRbRKSprFwJ110H++4Ln3wC99wD33wD228fdGQiItIMVDoiItJU1q2Dp56Cs8+GYcPgV78KOiIREWlGSrRFRBrTxInw8MPw2GO+NGTmTJWIiIgkKZWOiIg0hnDYdw455BBfh/3dd35cSbaISNJSoi0isi3WroW774Zu3eD55+H//g9mzYLu3YOOTEREAqbSERGRbfXUU5CTA/feC3vtFXQ0IiISIzSjLSLSUDNnwllnwc8/+5Ucx42DUaOUZIuISA1KtEVE6mvFCrj2WujVC956C6ZM8eOqwxYRkSiUaIuIbIlz8N//+jrs++6D88+HOXPgt78NOjIREYlhSrRFROrj5Zeha1coLoZHHoEdd2zyU4ZCIXJycsjKyiInJ4dQKNTk5xQRkcajRFtEJJpFi3y7vpISMIMXXvCrOx5wQLOcPhQKkZeXRzgcJj09nXA4TF5enpJtEZE4okRbRKS61av9Ko7du/vketIkP96pk0+4m0lBQQEpKSmkpaVhZqSlpZGSkkJBQUGzxSAiIttGibaISJU33oCePWHIEDjmGJg+Hc44I5BQSkpKSE1NrTGWmppKaWlpIPGIiEjDqY+2iEiV0aOhbVt4/3046qhAQ8nKyiIcDpOWlrZxrKKigszMzOCCEhGRBtGMtkic0wVz26CsDK68EiZM8I/vvtu37As4yQbIz8+nsrKS8vJynHOUl5dTWVlJfn5+0KGJiEg9KdEWiWO6YG4rrV/vO4d06wYPPvi/RLt9e78ATQzIzc2lsLCQjIwMysrKyMjIoLCwkNzc3KBDExGRejLnXNAxNIns7GxXXFwcdBgiTSonJ2eT8oLy8nIyMjIoKioKMLIYNn48/PWvMHkyHHEEPPAA9OkTdFQiIhKnzOwL51x2tG2q0RaJYyUlJaSnp9cY0wVzWzBhAixd6juKnH56s3YSERGR5KLSEUkqiVbPnJWVRUVFRY2xprpgLm5fu1Wr4Lbb4KWX/OMrr4SZM303ESXZIiLShJRoS9JIxHrm5rpgLi5fO+fglVegRw8YOtQvNgOQkgK12uaJiIg0BSXakjQScQGQ5rpgLu5eu2nT4He/gwEDoEMHKCrytdgiIiLNSDXakjQStZ45Nze3yTtRxN1rN2MGfPUV/OtfcPHF0Eq/6kREpPlpRluSRnPWMyeamH/t1q2Dhx7636z1aafBd9/B5ZcryRYRkcAo0ZakoQVAtl5Mv3ZjxsCBB8Jf/uJXdHTOX+S43XZBRyYiIklOibYkDS0AsvVi8rWbPx/++Efo3x9WrICXX4Y33lAnERERiRlasEZE4tPnn/sk+7rrID8f2rULOiIREUlCWrBGROKfczByJHzzDdx+Oxx0ECxcqBIRERGJWSodEZHYN2UKHHkknHkmhEKwerUfV5ItIiIxTIm2iMSun36Cyy7zFztOmwb/+Q9MmgRt2wYdmYiIyBYp0RaR2FVRAc89B3l5MGeO74ndsmXQUYmIiNSLarRFJLZ8+KHvIPLQQ7DbbjBvnkpEREQkLjVoRtvMWpjZvmbWz8zSmiooEUlCJSVw6qlw1FHw7rvw/fd+XEm2iIjEqXon2mb2F+B74CugCOgeGX/dzK5omvBEJOFVVMBNN8E++/gE+447YPp0yMgIOjIREZFtUq9E28wuAv4JvA6cAVRfEeIT4LTGD02k8YVCIXJycsjKyiInJ4dQKBR0SLJhAzz5pF82fdYsuPFGXewoIiIJob4z2tcA9zrnLgZeq7VtJpHZbZFYFgqFyMvLIxwOk56eTjgcJi8vT8l2ECZPhnPPhcpKaN8epk6FZ5+FXXcNOjIREZFGU99EOwt4t45t5YCKKCXmFRQUkJKSQlpaGmZGWloaKSkpFBQUBB1a8liyxHcOOfBAGD0aZs7046rDFhGRBFTfRHspkFnHtu7AokaJRqQJlZSUkJqaWmMsNTWV0tLSYAJKJuvWwT//CV27whNPwFVXwezZsN9+QUcmIiLSZOqbaL8FDDWzPaqNOTPrAlyNr90WiWlZWVlUVFTUGKuoqCAzMzOYgJKJGTz1FBxyCHz9NfzjH5rFFhGRhFffRPsmYA0wFfgAcMADwAxgPXBrk0Qn0ojy8/OprKykvLwc5xzl5eVUVlaSn58fdGiJ6bvv4Oyz/eqOLVv6/tihkO8uIiIikgTqlWg755YC2cCdQGvgO/xiN4VAX+fciiaLUKSR5ObmUlhYSEZGBmVlZWRkZFBYWEhubm7QoSWWlSt955AePeD11+GLL/z4dtv5mW0REZEkYc65oGNoEtnZ2a64uDjoMESSh3O+c8j118PixXDOOTB8OOy8c9CRiYiINBkz+8I5lx1tW337aM81s951bNvXzOZuS4AisUy9t+vJDF591SfWn37qa7KVZIuISBKrb412JtCmjm1tgd0bJRqRGNPcvbfjLqn/4Qe46CK/0Az4jiITJ/qLHkVERJJcvZdgx18AGU02sLwRYhGJOc3ZezuuFtSprIR774Vu3WDECD+DDdCpE7RoyK8VERGRxFXnX0Qzu9rM5pvZfHyS/WbV42q3JcC/gNHNFbBIc2rO3ttxs6DO6NG+//XgwXD44X5Vx0GDgo5KREQk5rTazLa5wIeRr88DioEltfZZA0wHHmv80ESCl5WVRTgcJi0tbeNYU/XeLikpIT09vcZYTC6o8/77sGEDvP02HH980NGIiIjErDoTbefcKGAUgPmWXLc650qaKS6RmJCfn09eXh7gk96Kioom673dnEl9g/z8M9x+Oxx3HOTkwK23wp13QkpKsHGJiIjEuPr20T5fSbYko+bsvR1zMB4oogAAIABJREFUC+ps2ABPPgndu0NBAYwb58fT0pRki4iI1EO9+2ibWQqQC3THdxqpzjnnbmvk2LaJ+mhLPAqFQhQUFFBaWkpmZib5+fnBLKjz+ef/396dx1dW1/cff30YyGAClImixgUSVHBDKkaoFq2kKsat+nOpK7hQcQkqaLAKKlSoxSAqRi1UFCuKSlVkuygKuOBCB6U6iIxgAlJC2YYtAQLM5/fHuaOZa5K5yeTkZHk9H4/zSO73fO+5nwyHmfd853u+X+jrg4svhr32guOPhz33nP86JEla4KZbR3u6OdoTL/Aw4CcUy/wlsGF7t4kpfUEFbWkx6u3tXRg7VV58MVxzTbGiyOte50oikiTNQrN/eg5QPAi5I0XI3gvYGTgauLL+vaTF6p574GMfgy9/uXh94IGwdi3st58hW5KkWWr2T9BnAB8Hrqu/Xp+Zw5n5IeC/gOPLKE5SyTLhrLPgiU8stk6/8MKifcstYdttKy1NkqTFrtmg/UDgusxcD4wCqyacOx941hzXJalsa9cWy/O96EVFsD73XDjppKqrkiRpyWg2aF8LPKj+/VXAcyec2xO4ey6LkjQPrrwSfvpTOO44+PWvYd99q65IkqQlpamHIYELgL8DTgdOAD4TEX8N3AvsW2+TtJCtXw9f/CKsW1fs6vj858PwMKxatcm3SpKkmWt2RPtw4HMAmfk54F1AK9ABfAx4z2wLiIhzIyIj4qgJbSfX2yY7fjfbz5KWrZ/+tFie74ADoFYrQjcYsiVJKlFTI9qZeRNw04TXnwY+vbkfHhGvBnaf5NRHgH9vaOsETgXO2NzPlZaN666DQw+Fr3wFHv7w4uurXw0Rm36vJEnaLE2NaEfE+RHx2CnO7RIR58/0gyNiFfAJ4JDGc5l5VWb+fOIBPKZ++ksz/Sxp2br5Zvj2t+Gww+B3v4PXvMaQLUnSPGl2jvazgO2mOLctxfztmToGWJOZp0bEV5vovx9wSWZeNovPkpaHTPjOd+BnP4NjjoHddoNrr3WKiCRJFZjJThRT7dX+KODOmXxoROxNEZzf0WT/vwUejaPZ0tR++9ti5ZCXvhTOPhvurP9vaciWJKkSU45oR8QbgTfWXyZwYkTc0dDtAcATgR80+4ER0UKxSsmxmXlFk2/bj2KFk1M3ce23AG8B2HHHHZstSVrcbrsNPvxhGBwsNpk5/nh429uKtbElSVJlphvRXg/cXz+i4fWG42aK1UjePIPPPJQioB/dTOeI2Bp4JXBW/aHMKWXmiZnZnZndO+ywwwxKkhaxu++G//zPYkWRtWvhoIMM2ZIkLQBT/mmcmV+iPlUjIi4A3paZm7W0XkTsCBwGHACsjIiVE06vjIjtgTsy8/4J7S8GtsdpI9Kf/eQnRbg+4QR4yEPgqqucIiJJ0gLT1BztzNxnc0N23c7A1sApwLoJB8B769/v1vCe/SmWFjxnDj5fWtyuvbZYOeQZzyi2TL/mmqLdkC1J0oIzZdCOiEdFxIsmae+JiIsj4s6I+H19XnSzLgX2meSAInzvA1w54bMeQrHz5Fcz894ZfI60tNxzDxx1FOy6a7Fc34c+VCzXt9NOVVe2SbVajZ6eHrq6uujp6aFWq1VdkiRJ82K6iZwfpFi7+swNDRGxK3AWxfzs7wK7AJ+LiBsz89ub+rDMvBW4sLE9inV9r87MxnOvBVbgtBEtd5nwhS8U26YPDEBnZ9UVNaVWq9HX10dLSwvt7e2MjIzQ19fH4OAgvb29VZcnSVKppps6shdwWkNbH9AC/H1mvoxiV8cf1NvLsD/FWtu/LOn60sK1Zg3stx/cdRdsvTVccgmcdtqiCdkAAwMDtLS00NbWRkTQ1tZGS0sLAwMDVZcmSVLpphvRfhhweUNbL/CrzLwYIDPXR8TnKVYembXMnHSrusycbHt2aWm75ZZiub7PfQ62264I3E996qKchz00NMSKFStYu3Yt99xzDytXruTBD34ww8PDVZcmSVLpphvRDoopIsWLiAdTPMx4UUO/64Bt5r40aZlZv74I17vsAp/9LBx4IPz+90XIXqS22247rrnmGsbHx1mxYgXj4+Ncc801bLvttlWXJklS6aYL2n+gmD6ywXMoNq65oKHfgylWBZG0OSLglFOKbdN/9Sv4zGfggQ+suqrNlplExJ+OzKk2mZUkaWmZLmh/CXhfRPRFxCuAj1AE6u819HsW8PtyypOWuKuvhv33h+uvL4L2WWfB+efDk55UdWVz4vbbb2ennXZiq6224r777mOrrbZip5124o47GjeZlSRp6ZlujvZngGcAx9df3w68NjPv2tAhIlqBVwOfKq1CaSkaG4OPfQyOOaYI2C97Gbz4xYtyHvZ0urq6GBkZYZdddvlT2+joKB0dHRVWJUnS/JhyRDszxzPz/wGPAp4KPDwzz57k/c8DPl1eidISc9pp8LjHwZFHFuH6d78rvi5B/f39jI+PMzo6SmYyOjrK+Pg4/f39VZcmSVLpNrkzZGYOZeYlmTk6ybk76+duK6c8aQk6/fRi5PrCC+HrX4cdd6y6otL09vYyODhIR0cH69ato6OjwzW0JUnLRizVB5O6u7tz9erVVZchwc03Fzs5HnhgMff69tuhrQ1WrKi6MkmStJki4pLM7J7s3CZHtCXN0n33wac/DY95DJxwAlxUXxlzu+0M2ZIkLQMGbakMF1wAT34yvPOdsMcecOml8La3VV2VJEmaR9OtOiJpts4/H+68E771LXjJS4qVRSRJ0rLiiLY0F0ZH4YMfhHPOKV5/4APw29/CS19qyJYkaZlqKmhHxBsj4ogpzh0REfvPaVXSYpEJp54Kj30sHHUU/OhHRfsDHlAckiRp2Wp2RPtdwM1TnLsBePfclCMtIpdeCs98JrzmNbDDDvDjH8O//VvVVUmSpAWi2TnajwYum+Lc5RSb2kjLyy9/WWw2c+KJ8KY3uZKIJEnaSLNB+z7gQVOc22GOapEWtnvvhc9+FrbZBt78ZnjDG4qt0//qr6quTJIkLUDNTh25GHjrFOfeCvz33JQjLVDnnQe77w7vfnfxPcAWWxiyJUnSlJod0T4a+H5E/AL4PPC/wMOBA4A9gOeUU55UsT/8AQ45BL7zHdh55+Lri15UdVWSJGkRaCpoZ+YPI+LlwCeBEyacGgZelpkXzn1p0gJw9dXwgx/ARz8KBx8MK1dWXZEkSVokmt6wJjO/A3wnInYFHgjclJlrS6tMqkImfOUrcM01xVrY++xTfL9qVdWVSZKkRWbGG9Zk5hWZ+VNDtpac1avhb/8WXv96OPtsuO++ot2QLUmSZmHKEe2I2A84OzNvrn8/rcz8zzmtTJovN9wA738/fPGL8OAHwxe+APvvXzzsKEmSNEvTTR05Gfgbio1qTt7EdRIwaGtxuvVW+PrX4T3vKbZR3267qiuSJElLwHRDdl3ApRO+n+7YucQapbl37rnFaiIAu+wCf/wjDAwYsqdRq9Xo6emhq6uLnp4earVa1SVJkrSgTRm0M/PqzByf8P20x/yVLG2G3/++WJ6vtxfOOgvWrSvanYc9rVqtRl9fHyMjI7S3tzMyMkJfX59hW5KkaTQ1CTUi7o+IPac495SIuH9uy5Lm2J13wvveB094Alx4IRxzDKxZY8Bu0sDAAC0tLbS1tRERtLW10dLSwsDAQNWlSZK0YDW7vF9Mc24FxRxtaeEaHy8ecnzta4s1sR/60KorWlSGhoZob2/fqK21tZXh4eFqCpIkaRGYdkQ7IraIiBUb+tZfTzzagF7gptIrlWbq4ovhzW+G+++H9nZYu7ZYWcSQPWNdXV2MjY1t1DY2NkZnZ2c1BUmStAhMGbQj4sPAvcA4xYj1RfXXE4/bgQ8Bp5VeqdSs66+HN74R9toLzjkHrrqqaHeayKz19/czPj7O6Ogomcno6Cjj4+P09/dXXZokSQvWdFNHLqx/DYowfRJwbUOfe4DfAmfNeWXSTN17L3zyk/CRj8Ddd8Ohh8Lhh8O221Zd2aLX29vL4OAgAwMDDA8P09nZSX9/P729vVWXJknSghWZm55eXR/d/o/MvK78kuZGd3d3rl69uuoyNJ/Gx+FJTyqW6/v4x+Exj6m6IkmStMRFxCWZ2T3ZuaZWHcnMIzPzuvq87CdGxN/V52dL1briCthvP7jjDmhpgZ/9DM44w5C9gLj+tiRpuWp6j+mIeAdwPfBr4Hxg13r76RHxznLKk6Zw223w3vfCbrvBd74Dl9b3VnIe9oLi+tuSpOWs2XW0/wn4FHA68Eo2Xu7vx8DL5r40aRKZxTJ9u+wCxx1XjGavXQvPeEbVlWkSrr8tSVrOml1H+xDg45n5vgnL/W3wO8ClBzR/vvpVeNSjihVFnvKUqqvRNFx/W5K0nDU7daQL+O4U50aB7eemHGkS110Hb3oTXHMNRMBpp8FFFxmyFwHX35YkLWfNBu2bgM4pzu0K/O+cVCNNdPfdxS6Ou+wCX/kK/PznRfuqVUXg1oLn+tuSpOWs2aB9FvChiNh5QltGxIOAgynmbktz54wz4AlPgA98AJ7zHLj8cnjlK6uuSjO0Yf3tjo4O1q1bR0dHB4ODg66/LUlaFpqdo304sA+wBvgFxU6RxwOPBW4A/qWU6rR8nXkmrFwJ3/teEbS1aPX29hqsJUnLUrPraN8EdAMfBbYCrqII6YPA0zLzttIq1PJw661w8MFw8cXF6+OOg//5H0O2JElatJod0SYz7wA+Uj+kuXH//XDSSXDYYXDzzfCIR8Cee7ptuiRJWvSa3rBGmnMXXQRPfSoceCA89rFwySXwnvdUXZUkSdKcaHpEOyL2B14N7Ahs3XA6M/NRc1mYloEf/hBuvBFOPRX+8R9dSUSSJC0pTQXtiPggcCTFw5CXAveUWZSWqLvvhmOPLVYTeelLi9Hrd70L2tqqrkySJGnONTui/WbgU5l5cJnFaInKhG9/uwjWw8Nw0EFF0F65sjgkSZKWoGbnaD8QOLPMQrREXXYZPPvZ8LKXwTbbwA9+AMcfX3VVkiRJpWs2aP8Q2L3MQrRE/eY38Ktfwac/XXzt6am6IkmSpHnR7NSRdwPfioibgXOAWxo7ZOb6uSxMi9T998OJJxbTRd7+9uIhx333LbZNlyRJWkaaHdFeCzwR+CLwf8C9Dcd4KdVpcfnhD2GPPYqAfe65RdiOMGRLkqRlqdkR7X+h2HZd+kt//CO8973wjW/AjjvCaacVc7Jdrk+SJC1jTQXtzDyi5Dq0mF17LZx1FhxxBPT3Q2tr1RVJkiRVbsY7Q0bENhGxU0RsVUZBWgQyi1HrI48sXj/tacWo9oc/bMjeTLVajZ6eHrq6uujp6aFWq1VdkiRJmqWmg3ZEvDAifgncBvwB2K3e/vmIeE1J9Wmh+fWvYZ994JWvhDPOgHvqexe1t1db1xJQq9Xo6+tjZGSE9vZ2RkZG6OvrM2xLkrRINRW0I+IlwHeAm4D3ARMn3w4B+899aVpQbrmleMjxyU+GNWvgc5+Diy92w5k5NDAwQEtLC21tbUQEbW1ttLS0MDAwUHVpkiRpFpod0f4w8MXMfC7wyYZzayhWJNFSdvvt8OUvwzveAWvXwlvfCitWVF3VkjI0NERrw9Sb1tZWhoeHqylIkiRtlmaD9uOAr9e/b1x9ZB3FzpFaai64AN75zmJOdmcnXH11sauj00RK0dXVxdjY2EZtY2NjdHZ2VlOQJEnaLM0G7duBB01xrhO4cU6q0cIwPAwvf3mxi+OZZ8INNxTtBuxS9ff3Mz4+zujoKJnJ6Ogo4+Pj9Pf3V12aFhEfqJWkhaPZoH0e8P6I2H5CW0bESqAP8HfypeCuu4qVQx73OKjV4Kij4PLL4SEPqbqyZaG3t5fBwUE6OjpYt24dHR0dDA4O0tvbW3VpWiR8oFaSFpbI3PQ+NBHRCVxMMW3kHGA/4L+AJwF/BXRn5nWlVTkL3d3duXr16qrLWFxuvx123bVYVeRjH4NHPKLqiiTNQE9PDyMjI7S1tf2pbXR0lI6ODs4///wKK5OkpSsiLsnM7snONTWinZnDwB7AWcBzgPuBZwI/B/ZaaCFbM3DppfCmN8G998J228Fll8FXv2rIlhYhH6iVpIWl6XW0M/PazHxzZj4iM1sysyMz35iZfyyzQJXkppuKlUOe8pRiHvYVVxTtzsOWFi0fqJWkhWXGO0NqkbvvvmLlkMc8Bj7/eTjooGK5vie6QqPK4cN588cHaiVpYWl2jvYXpjm9nmK3yEuAb2Xm3TMqIOJcYF/g6Mw8vOHc3wBHAH8DbEWxI+XRmfm1TV3XOdpTuP9+2GOP4gHHT34SHv/4qivSErbh4byWlhZaW1sZGxtjfHzchzxLVKvVGBgYYHh4mM7OTvr7+/21lqQSTTdHu9mgPUTx0OP2wH0UO0Q+CNgSuLXebXvgKmCfzLy2ycJeDRwHPJSGoB0RLwC+DXwV+AYwDjweuD0zT97UtQ3aE/zhD3DkkUWwXrWq2OVx1SqI2PR7pc3gw3mSpKVusx+GBF5DMWr9MmDrzHwYsDXwCoo1tl8E7Flv+2iTRa0CPgEcMsm5bYEvAp/NzDdk5jmZ+f3MPL6ZkK26O++Eww4rRq2/+U245JKivb3dkK154cN5kqTlrNmg/QngmMz8dmauB8jM9Zn5TeAY4BOZuZoiZD+nyWseA6zJzFMnOfcKYAfg401eSxNlwle+UizV96//Cq94RfGw47OfXXVlWmZ8OE+StJw1G7R3p5gWMpmrgA1P0v0WWLWpi0XE3hRrcb9jii57A7cAu0XEbyLivoj4Y0R8OCJWNFnz8hUBX/sadHTAT38KX/4yPPzhVVelZciH8yRJy1mzQft64OVTnHsF8H/177cD1k13oYhoAU4Ajs3MK6bo9jCglWJ+9snAs4EvAR8Ejp3m2m+JiNURsfrGG5fZrvA33AAHHghXXlm8/vKX4eKL4WlPq7YuLWvudilJWs62bLLfp4DjIuJhFDtC3gA8mCJkPx94d73fM4BfbeJahwIPAI6eps8WFPO9D8vM4+ptF0bEA4F3RMQRmXlb45sy80TgRCgehmzmB1v0xsdhcLB42HFsDJ7+dHj0o2H77auuTAKKsG2wliQtR00F7cz8ZETcCXwIeMGEU9cC/5SZJ9Vffwa4a6rrRMSOwGHAAcDKiFg54fTKiNgeuAO4ud52XsMlvge8FXgC8NNmal/SvvtdeNe7ivnXz3tesarIrrtWXZUkSZJofkSbzPx8RJwEPALoAEaAa3PC+oD1rdqnszPFSPUpk5x7b/14MnDZJq6zvsmyl7ZzzinWxT7zTHjBC1xJRJIkaQGZ0c6QWfhjZl5c/zrT6RmXAvtMckARvvcBrgROr7ft2/D+5wF3A2tm+LlLwx13wPveBz/6UfH66KNhzRp44QsN2ZIkSQvMlCPaEbEfcHZm3lz/flqZ+Z9N9LkVuHCSzwK4OjM3nFsTEScD/xIRWwC/pHgg8gDgI5l556Y+a0lZvx5OOaUI2ddfD9tsA898ZvFVkiRJC9J0U0dOptj6/Ob699NJYJNBe4YOBP4XOAh4CDAMHJKZn5rjz1nY/vu/4aCD4Be/gD33hNNPh732qroqSZIkbcJ0QbuLYh72hu9Lk5l/Me8hM8eBw+vH8nXRRXD11XDyyfD618MWM5rtI0mSpIrEzKdZLw7d3d25evXqqsuYufFx+NSn4JGPhFe9Cu69F+66C7bbrurKJEmS1CAiLsnM7snOzWh4NCIeHRGviYj++tdHzU2JAuDss+GJT4RDD4Xvf79o22qrjUJ2rVajp6eHrq4uenp6qNVqFRUrSZKk6TQVtCNi64j4AnA5xeogx9S//i4iPt+wHrZmau1aeP7zi9VDttgCajX4/Of/olutVqOvr4+RkRHa29sZGRmhr6/PsC1JkrQANTuifSzwWuDDwKOBbetfjwBeDwyUUdyycfnlxVzsj38cfv3rYvOZSQwMDNDS0kJbWxsRQVtbGy0tLQwM+MsvSZK00DS7Yc2rgCMz818ntP0BOLq+NN/BwDvnuLala/364uHG22+Hd78bXvxiGBqC9vZp3zY0NER7Q5/W1laGh4fLq1WSJEmz0uyI9krg4inO/QJomZtyloGf/axYnu/Nb4azzoLMYrOZTYRsgK6uLsbGxjZqGxsbo7Ozs6RiJUmSNFvNBu3vA8+d4txzgfPnppwlbGSkWJ7v6U+H664rNqA577wZ7ejY39/P+Pg4o6OjZCajo6OMj4/T399fYuGSJEmajWaD9nHAKyPiMxHxrIh4XP3rZ4FXAsdGxM4bjvLKXcSuvx6++U34wAfgiivgta+d8bbpvb29DA4O0tHRwbp16+jo6GBwcJDe3t6SipYkSdJsNbWOdkSsn/By4htikjYyc8Xml7Z5Kl9HOxPOPBMuvhiOOqpou+WWpqaISJIkaXGYbh3tZh+GfOMc1rP0XX558ZDj974HT3hCMYrd2mrIliRJWkaaCtqZ+aWyC1kSbrsNjjwSPv1paGuDT34S3v72YtMZSZIkLSsz2hlyMhGxRUQ4VAswOgonnQRvehP8/vfwrncZsiVJkpapKYN2RNwSEXtMeB0RccYkDzs+FbixrAIXlYc9rFgP+4QTYIcdqq5GkiRJFZpuRHt7Np5asgXwwnq7puI8bEmSJDEHU0ckSZIk/SWDtiRJklQCg7YkSZJUgk0t7/fwCQ8/rpjQduuEPo+Y+7IkSZKkxW1TQfu/Jmk7veF10LAzpCRJkrTcTRe03Q1SkiRJmqUpg7a7QUqSJEmz58OQkiRJUgkM2pIkSVIJDNqSJElSCQzakiRJUgkM2tIyVKvV6Onpoauri56eHmq1WtUlSZK05Bi0pWWmVqvR19fHyMgI7e3tjIyM0NfXZ9iWJGmOGbSlZWZgYICWlhba2tqICNra2mhpaWFgYKDq0iRJWlIM2tIyMzQ0RGtr60Ztra2tDA8PV1OQJElLlEFbWma6uroYGxvbqG1sbIzOzs5qClpknN8uSWqWQVtaZvr7+xkfH2d0dJTMZHR0lPHxcfr7+6subcFzfrskaSYM2tIy09vby+DgIB0dHaxbt46Ojg4GBwfp7e2turQFz/ntkqSZ2LLqAiTNv97eXoP1LAwNDdHe3r5Rm/PbJUlTcURbkprk/HZJ0kwYtCWpSc5vlyTNhEFbkprk/HZJ0kxEZlZdQym6u7tz9erVVZchSZKkJSwiLsnM7snOOaItSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgrSnVajV6enro6uqip6eHWq1WdUmSJEmLhkFbk6rVavT19TEyMkJ7ezsjIyP09fUZtiVJkppk0NakBgYGaGlpoa2tjYigra2NlpYWBgYGqi5NkiRpUTBoa1JDQ0O0trZu1Nba2srw8HA1BUmSJC0yBm1Nqquri7GxsY3axsbG6OzsrKYgSZKkRcagrUn19/czPj7O6Ogomcno6Cjj4+P09/dXXZokSdKiYNDWpHp7exkcHKSjo4N169bR0dHB4OAgvb29VZcmSZK0KBi0tUmZWXUJkiRJi45BW5NyeT9JkqTNY9DWpFzeT5IkafMYtDUpl/eTJEnaPAZtTcrl/SRJkjaPQXuZqNVq9PT00NXVRU9PzybnWru8nyRJ0uYxaC8Ds3mw0eX9JEmSNk9UvXRbRJwL7AscnZmH19s6gaEp3rIqM2/d1HW7u7tz9erVc1XmotbT08PIyAhtbW1/ahsdHaWjo4Pzzz+/wsokSZIWt4i4JDO7Jzu35XwXM1FEvBrYfZouHwXOaGi7o7yKlqahoSHa29s3avPBRkmSpHJVFrQjYhXwCeBg4KtTdPtDZv58/qpamrq6uv5iRNsHGyVJkspV5RztY4A1mXlqhTUsCz7YKEmSNP8qCdoRsTewH/COTXT9aETcFxG3RcQZEbHbPJS35PhgoyRJ0vyb96kjEdECnAAcm5lXTNHtnnqf7wE3Ao8FPgD8NCL2zMzLp7j2W4C3AOy4445zXfqi1tvba7CWJEmaR1WMaB8KPAA4eqoOmTmSmW/NzG9l5o8z8z+AZwIJHDbN+07MzO7M7N5hhx3mvHBJkiSpWfM6oh0RO1IE5QOAlRGxcsLplRGxPXBHZt7f+N7M/GNE/AR46vxUK0mSJM3efI9o7wxsDZwCrJtwALy3/v2m5mFXu/C3JEmS1IT5nqN9KbDPJO0XUITvk4ArJ3tjfTR8b+D00qqTJEmS5si8Bu36jo4XNrZHBMDVmXlh/fXHKUbbf0bxMOSuwPuB9Uwzt1uSJElaKCrdGXIalwFvA94AbAPcDJwPHDnNSiWSJEnSgrEggnZmRsPrLwBfqKgcSZIkabNVuTOkJEmStGQZtCVJkqQSGLQlSZKkEhi0JUmSpBIYtCVJkqQSGLQlSZKkEhi0JUmSpBIYtCVJkqQSGLQlSZKkEhi0JUmSpBIYtCVJkqQSGLQlSZKkEhi0JUmSpBIYtCVJkqQSGLQlSZKkEhi0JUmSpBIYtCVJkqQSGLQlSZKkEhi0JUmSpBIYtOdIrVajp6eHrq4uenp6qNVqVZckSZKkChm050CtVqOvr4+RkRHa29sZGRmhr6/PsC1JkrSMGbTnwMDAAC0tLbS1tRERtLW10dLSwsDAQNWlSZIkqSIG7TkwNDREa2vrRm2tra0MDw9XU5AkSZIqZ9CeA11dXYyNjW3UNjY2RmdnZzUFSZIkqXIG7TnQ39/P+Pg4o6OjZCajo6OMj4/T399fdWmSJEmqiEF7DvT29jI4OEhHRwfr1q2jo6ODwcFBent7qy5NkiRJFYnMrLqGUnR3d+fq1aurLkOSJElLWERckpndk51zRFuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqUEU9WAAAWB0lEQVSBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqgUFbkiRJKoFBW5IkSSqBQVuSJEkqQeVBOyLOjYiMiKOm6fPv9T6nzGdtkiRJ0mxVGrQj4tXA7pvo87fA64Db56UoSZIkaQ5UFrQjYhXwCeCQafpsBZwAHA2sm6fSJEmSpM1W5Yj2McCazDx1mj79wArg2PkpSZIkSZobW1bxoRGxN7Af00wbiYhHA4cDL8jMeyNivsqTJEmSNtu8j2hHRAvFdJBjM/OKabp+DvhWZl4wg2u/JSJWR8TqG2+8cXNLlSRJkmatiqkjhwIPoJh3PamIeB3wVOA9M7lwZp6Ymd2Z2b3DDjtsXpWSJEnSZpjXqSMRsSNwGHAAsDIiVk44vTIitgcSOI5iDvc99TYo/lKwVf31aGbeO4+lS5IkSTMy3yPaOwNbA6dQrCKy4QB4b/37LmAH4F8b+jwSeGX9+xfMa9WSJEnSDM33w5CXAvtM0n4BRfg+Cbhyij5fA35DMeVkTVkFSpIkSXNhXoN2Zt4KXNjYXl9R5OrM3HBusj53A/83oY8kSZK0YFW+BbskSZK0FFWyjnajzNzkItmZ2TkPpUiSJElzwhFtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtLVq1Wo2enh66urro6emhVqtVXZIkSdKfGLS1KNVqNfr6+hgZGaG9vZ2RkRH6+voM25IkacEwaGtRGhgYoKWlhba2NiKCtrY2WlpaGBgYqLo0SZIkwKCtRWpoaIjW1taN2lpbWxkeHq6mIEmSpAYGbS1KXV1djI2NbdQ2NjZGZ2dnNQVJkiQ1MGhrUerv72d8fJzR0VEyk9HRUcbHx+nv76+6NEmSJMCgrUWqt7eXwcFBOjo6WLduHR0dHQwODtLb21t1aZIkSQBEZlZdQym6u7tz9erVVZchSZKkJSwiLsnM7snOOaItSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgLUmSJJXAoC1JkiSVwKAtSZIklcCgLUmSJJWg8qAdEedGREbEURPanlJv/9+IuDsiro+IcyLiaVXWKkmSJDWr0qAdEa8Gdp/k1PbAlcB7gH2Bg+ptP4yIPeevQkmSJGl2tqzqgyNiFfAJ4GDgqxPPZeYPgB809D8XuAl4PXDxPJUpSZIkzUqVI9rHAGsy89Qm+48C9wD3lVeSJEmSNDcqGdGOiL2B/Zh82sjEflsAK4AO4J/rzf9RbnWSJEnS5pv3Ee2IaAFOAI7NzCs20f0bwDhwNfAy4PmZ+dtprv2WiFgdEatvvPHGOatZkiRJmqkqpo4cCjwAOLrJvntShOw1wFkR0T1V58w8MTO7M7N7hx12mJNiJUmSpNmY16AdETsChwEfBFZGxPYRsX399IbXKzb0z8w/ZOZ/Z+a3gF7gBuCov7iwJEmStMDM94j2zsDWwCnAugkHwHvr3+822Rszcxz4NfDo8suUJEmSNs98Pwx5KbDPJO0XUITvkyjWz/4LEdEKdAObmtctSZIkVW5eg3Zm3gpc2NgeEQBXZ+aF9dcnALcAqynWzt4J6KNYfeT181OtJEmSNHuVbVizCb8ADgDeArQB/1tve3Nm/qbKwiRJkqRmVLoF+waZGZl5+ITXX8jMp2fmAzNz68x8VGa+xpANtVqNnp4eurq66OnpoVarVV2SJEmSJrEggraaU6vV6OvrY2RkhPb2dkZGRujr6zNsS5IkLUAG7UVkYGCAlpYW2traiAja2tpoaWlhYGCg6tIkSZLUwKC9iAwNDdHa2rpRW2trK8PDw9UUJEmSpCkZtBeRrq4uxsbGNmobGxujs7OzmoIkSZI0JYP2ItLf38/4+Dijo6NkJqOjo4yPj9Pf3191aZIkSWpg0F5Eent7GRwcpKOjg3Xr1tHR0cHg4CC9vb1VlyZJkqQGkZlV11CK7u7uXL16ddVlSJIkaQmLiEsys3uyc45oS5IkSSUwaEuSJEklMGhLkiRJJTBoS5IkSSUwaEuSJEklMGhLkiRJJTBoS5IkSSUwaEuSJEklMGhLkiRJJTBoS5IkSSUwaEuSJEklMGhLkiRJJTBoS5IkSSUwaEuSJEklMGhLkiRJJTBoS5IkSSUwaEuSJEklMGhLkiRJJTBoLzK1Wo2enh66urro6emhVqtVXZIkSZImYdBeRGq1Gn19fYyMjNDe3s7IyAh9fX2GbUmSpAXIoL2IDAwM0NLSQltbGxFBW1sbLS0tDAwMVF2aJEmSGhi0F5GhoSFaW1s3amttbWV4eLiagiRJkjQlg/Yi0tXVxdjY2EZtY2NjdHZ2VlOQJEmSpmTQXkT6+/sZHx9ndHSUzGR0dJTx8XH6+/urLk2SJEkNDNqLSG9vL4ODg3R0dLBu3To6OjoYHBykt7e36tIkSZLUIDKz6hpK0d3dnatXr666DEmSJC1hEXFJZnZPds4RbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEkZlV11CKiLgRuLqCj34QcFMFn6uFyftBjbwnNJH3gybyflicdsrMHSY7sWSDdlUiYnVmdlddhxYG7wc18p7QRN4Pmsj7Yelx6ogkSZJUAoO2JEmSVAKD9tw7seoCtKB4P6iR94Qm8n7QRN4PS4xztCVJkqQSOKItSZIklcCgLUmSJJXAoD0DEbFvRJwfEddHxD0RcW1EfCMiHt/Q75ER8V8RcVtE3B4R34qIHauqW/MnIs6NiIyIoxraV0XE5yPipogYjYjvR8RuVdWpckTEs+r//RuPWxv6eT8sIxHx/Ij4UUTcWf8zYXVE9Ew47/2wDETEhVP8/pARce6Eft4PS8iWVRewyLQDlwCfBW4EdgT+Gfh5ROyWmVdHRCtwPnAPsD+QwFHABRHxpMwcraZ0lS0iXg3sPkl7AGcCncBBwDrg/RT3xF9n5rXzWafmxTuB/57w+r4N33g/LC8RcSAwWD8+QjHA9ddAa/2898Py8XZgu4a2pwHHAWeA98OSlJkem3EAu1KE6ffUX78LuB949IQ+XRR/0B5Sdb0epd0Hq4DrgVfX74ejJpz7h3rbPhPa/gq4BTi+6to95vQ+eFb9v/Wzp+nj/bBMDoqwdBfwbu8Hjyn++59EMTDX7v2wNA+njmy+m+tfN4xYvRj4eWZeuaFDZg4BF1H8D6Sl6RhgTWaeOsm5FwPXZeYFGxoy8zaKUQvvieXH+2H5eBOwHvj3afp4PyxT9X8BfwVwZmbeUm/2flhiDNqzEBErIqIlIh4DnEAxkrkhYD0BWDPJ2y4DHj9Juxa5iNgb2A94xxRdprsndoyIbcqqTZX5SkTcHxE3R8RXG57R8H5YPvYGfge8KiKuioj7IuLKiJj4e4X3w/L1UmBb4EsT2rwflhiD9uz8guKfetYCTwJ6MvOG+rl2ijlVjW6hmF6gJSQiWij+snVsZl4xRbfp7gnwvlhKbgM+DhwA9FDMyX028LOIeHC9j/fD8vEw4DHAAPBvwHOB84DBiHhXvY/3w/K1H3ADUJvQ5v2wxPgw5Oy8nuKBhp2B9wLnRcTemTlcaVWqwqHAA4Cjqy5E1cvMXwG/mtD0w4j4EXAxxQOSh1dSmKqyBcWI5Rsy81v1tvMjohN4f0QcX1VhqlZEPIziL+Gfysz7NtVfi5cj2rOQmZdn5i/q83H/HtiGYvURKP4mOtnfOKf6W6oWqfp0gMOADwIrI2L7iNi+fnrD6xVMf0+A98WSlpm/pPjXr6fWm7wflo8Nz/Cc19D+PeAhQAfeD8vV6ygy2Jca2r0flhiD9mbKzFuBK4FH15suo5hj1ejxwG/nqy7Ni52BrYFTKH7z23BA8S8d64DdmP6euCYz7yy/VC0AWf/q/bB8XLaJ8+vxfliu9gf+JzP/p6Hd+2GJMWhvpoh4CPBY4Kp60xnA30TEzhP6dAJ/Wz+npeNSYJ9JDijC9z4Ufwk7A3h4RPzdhjdGxHbAi/CeWPIioptiGdCL603eD8vHt+tf921ofx5wbWZej/fDslP/PeHx/OVoNng/LDmRmZvuJQAi4tvAL4FfA7cDuwAHAw8F9szMtRHRBvwPxdqph1OMYn2EYp7ek/zb6NIXEQkcnZmH119vAfwEeCTQz583IHgSsHtm/rGqWjW3IuIrwBDF7xO3Ak+m+G89BuyRmTd5Pywf9c1HfkCxkdVhwB8olnM7AHhjZp7s/bD81Ofmvw14+ISFFDac835YYhzRnpmfAy+h+Fvo2cAhwA+Bv87MtQBZ7PzYQzEn88vAhj94ewzZy1NmrgdeSDFP87MUo1z3U2xI4G+aS8sainVwvwh8F3g38C1gr8y8CbwflpMsRrJeAnwNOBI4C9gLeG1mnlzv4/2wjETEVhQbm53bGLLB+2EpckRbkiRJKoEj2pIkSVIJDNqSJElSCQzakiRJUgkM2pIkSVIJDNqSJElSCQzakiRJUgkM2pKWnYh4WkR8IyKui4jxiLg5Is6LiP0jYkW9zxsiIiPi0VXXO1P1uo8o8fpH1DdmKuPaG37dO2fwnrX19/xDGTVJ0mwZtCUtKxHxbuAioB14H/Bs4E0Um0x9jmKzCE3v88DTqi4CICKeDjym/nK/KmuRpEZbVl2AJM2XiHgmcBwwmJnvbDj9nYg4Dmibg8/ZCrgvl+iOYJl5LXBt1XXU7Q/cB5wPvDAi2jPzloprkiTAEW1Jy8v7gFuAQyc7mZlXZeavG5ofFBFfiYjb61NNjo+IrTecjIjO+rSFt0fExyLiOuAeYPsoHBwRV9SnqIxExGBEbDfxA+rvPyoi3hkRQxFxR0T8MCKe0NCvqetNJiKeFxE/i4i7IuK2iDg9InZt6LOiXsdIRIxFxPkR8djGqSiTTR2JiC0j4n0R8duIuDsiboyIcyPisfXzW0fEJyJiTUTcGRHXR8SZG87PRv2/wyuB7wEDQAvF9taN/Zr6uep9d4+IMyJiXf3X6qKIeMZsa5S0vBm0JS0L9bnX+wDfy8y7Z/DWLwNXAf+PYmrJO4D3T9LvMGAX4C3AS4G7gaMpRtDPA14EfAx4A3B2RDT+/vs64AXAu4A3AjtSjLJP/JfHmVzvTyLiecDZwJ3APwJvA54I/CQiHj6h65HAB4D/BP6BIsCeMdV1G3ytXt85wEuAfwJ+C3TUz68EtgWOqv+cbwO2Bn4WEQ9t8jMa/QOwfb3e8ylG2SebPtLUzxURewA/pZhW9E/Ay4Cbge9HxFNmWaOk5SwzPTw8PJb8ATwESOCjTfZ/Q73/kQ3tZwFrJ7zurPf7JRAT2tspRrZPbnj/6+r9XzyhLYHfA1tNaHt5vf3ps7zeERNer65ff8sJbV3AvcBx9derKIL4Zxuuf8gk1zui+OPjT6976n3eOYP/HiuAVuAO4OBJft07m7jGOcCtwNb11x+tv/exE/rM5Of6AXA50NJQ5+XA6VXfwx4eHovvcERbkqZ3dsPr31CMNjc6PTMnTqf4G4qpDKc09PsaxZziv2toPy8z7234HCZ81kyvB0BEtAF7AF/PzPs2tGfmEMVDoRvetxvF/PTTGi7xX5Ndt8FzKULrf0zXKSJeGRG/iIhb6zWPAtsAu073vimu9dD6556Wf/4Xii/Vv+4/oWtTP1dEPIDi1+I0YH19KsyWQADfB5450xolyaAtabm4GbgL2GmG72t8sO4eimkQjUYaXrdP1l4PuzdPOD/d50AxvWI219tgFUVYbKwP4PoJ79swxeOGhj7/N8V1J3ogcEtm3jVVh4h4EfB1itHh1wB7AU8FbuTPP+NMvI5itPk7EbF9RGxP8fNcCrxuwlSaZn+u9vr1Pkgx0j/x6ANWTTc9R5Im46ojkpaFzLwvIi4EnhMRKzPznk29Z6Yf0fB6Q3B+KHDZhsb6KOkD+ctgvSmzvd66em2TzYN+6IT3bQjiD554fYopN5tyE9AeEQ+YJmy/CrgyM98wofatmPovCJuyYdT6zCnO91CMRDf7c90KrAc+QzGX+y9k5vpZVSpp2fJv55KWk3+jCKUfm+xkRHRFxJPm6LN+DoxTBMyJ/pFikOPC+bheZo4ClwCvqD8QCkBE7AQ8fcL7fkMxleMVDZdofD2Z71GMmh8wTZ9WiukiE72eYhR5RuoPLT4ROIHiAdeJx74U/xqwIYg39XPVf51+DOwO/DIzVzceM61TkhzRlrRsZOaPIuIQ4LiIeDxwMnANxfSKv6cIiq8BGpf4m81n3RIRHwfeHxGjFA/uPY5i1Y2f8Jdzv8u83gfr58+KiM9SzIs+ErgN+Hj9+usi4pPAByLiDorR4D2AN9evMeVobmZeEBHfpPh1fSTFCiBbUcxrPjszLwTOBV4SEZ+geKC0GziIYiR5pvanGKU/pj7XfCMRcTrw0ojYZoY/1yHAj4DvRsRJFKPhD6r3X5GZ/zyLWiUtY45oS1pWMvOTwN4UAe9YilB4MkVoPZCppyLMxmEU4a2XIlz+M8W0hBfMchrCrK6XmedSLKm3PfAN4N8p5krvnZnXTej6YYqVO/anWP6ul2IVEChC+XReRbEayUvq7/0C8AT+PHXjPyiW//tHil/j51MsUbip626kPt3kNcAFk4XsupMoHoB8ef11Uz9XZv6SYt74zcDxFCP1n6J4oPJHM6lTkqC+FJUkSZOJiJdTrMTxzMz8cdX1zJWl+nNJWlgM2pIkACJiL4qR719QbLjzFIpR8yso1vNelH9gLNWfS9LC5xxtSdIGd1LMq34HsB3FknjfAN6/yMPoUv25JC1wjmhLkiRJJfBhSEmSJKkEBm1JkiSpBAZtSZIkqQQGbUmSJKkEBm1JkiSpBP8f5TTISyoZvBYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_known_predicted_ages([age[0] for age in predicted_ages.values()], [age[1] for age in predicted_ages.values()], 'TBS Epigenetic Clock')"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}