Estimation Example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "56da4b6f-d30b-439b-bed2-72310265891c",
   "metadata": {
    "id": "56da4b6f-d30b-439b-bed2-72310265891c"
   },
   "source": [
    "# Estimate recombination and mutation rates from a simulated dataset\n",
    "This notebook shows how to use recombulator-x to estimate recombination and mutation rates from a simulated dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "95cfc9f8-f5d2-4fae-80cb-f06d6bc66b7e",
   "metadata": {
    "id": "95cfc9f8-f5d2-4fae-80cb-f06d6bc66b7e"
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "import pandas\n",
    "import recombulatorx\n",
    "import recombulatorx.testing"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f0f47596-4e80-47b8-86aa-37a31545cc19",
   "metadata": {
    "id": "f0f47596-4e80-47b8-86aa-37a31545cc19"
   },
   "source": [
    "## Simulation of STR pedigree dataset\n",
    "Recombulator-x includes functions to simulate pedigrees and their STR genotypes, a useful feature when testing new approaches.\n",
    "In this section we generate a pedigree file simulating 100 type I and 100 type II families with given recombination and mutation rates for six markers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f3c705b-fd86-4831-b407-4b900074c635",
   "metadata": {},
   "source": [
    "We define the simulation parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3861fb77-6661-4f89-81bf-e9e39980af35",
   "metadata": {
    "id": "3861fb77-6661-4f89-81bf-e9e39980af35"
   },
   "outputs": [],
   "source": [
    "my_recomb_rates = numpy.array([0.05, 0.3, 0.01, 0.2, 0.01])\n",
    "my_mut_rates = numpy.array([1e-4, 1e-4, 1e-2, 1e-3, 1e-4, 1e-4])\n",
    "n_fam_I = 100\n",
    "n_fam_II = 100\n",
    "ped_path = 'testsim.tsv'\n",
    "\n",
    "numpy.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81a1834c-baf9-4519-92f7-dc1e2b2e68ad",
   "metadata": {},
   "source": [
    "The dataset is generated through the functions *generate_complex_families* and *individuals2ped*. We use the **generate_complex_families** function to simulate random individuals from a given number of informative families, following the given rates and the **individuals2ped** function to write the generated individuals to a pedigree file in TSV format, which will be saved to *ped_path*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "017fcd1e-5043-4c7e-8221-52ffa212370a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('FAM_0_I', 'GRANDFATHER', None, None, array([[12., 29., 39., 43., 59., 68.]]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "individuals = list(recombulatorx.testing.generate_complex_families(n_fam_I, n_fam_II, my_recomb_rates, my_mut_rates))\n",
    "individuals[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b388aac9-ca79-405e-bf47-ec31f557cf74",
   "metadata": {
    "id": "b388aac9-ca79-405e-bf47-ec31f557cf74"
   },
   "outputs": [],
   "source": [
    "recombulatorx.testing.individuals2ped(\n",
    "    path=ped_path, \n",
    "    marker_names=[f'M{i+1}' for i, _ in enumerate(my_mut_rates)], \n",
    "    individuals=individuals,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "oO5dHzIFyMNO",
   "metadata": {
    "id": "oO5dHzIFyMNO"
   },
   "source": [
    "The created pedigree file is a table in textual format that can be read with the pandas module. See the online documentation (Usage section) for a detailed description of the PED format. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5d7fedbd-0190-4302-8e5f-15128af9af00",
   "metadata": {
    "id": "5d7fedbd-0190-4302-8e5f-15128af9af00",
    "outputId": "9246d532-81fb-4682-fb66-868dad6b416d"
   },
   "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>FID</th>\n",
       "      <th>IID</th>\n",
       "      <th>PAT</th>\n",
       "      <th>MAT</th>\n",
       "      <th>SEX</th>\n",
       "      <th>PHENO</th>\n",
       "      <th>M1-A1</th>\n",
       "      <th>M1-A2</th>\n",
       "      <th>M2-A1</th>\n",
       "      <th>M2-A2</th>\n",
       "      <th>M3-A1</th>\n",
       "      <th>M3-A2</th>\n",
       "      <th>M4-A1</th>\n",
       "      <th>M4-A2</th>\n",
       "      <th>M5-A1</th>\n",
       "      <th>M5-A2</th>\n",
       "      <th>M6-A1</th>\n",
       "      <th>M6-A2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>FAM_0_I</td>\n",
       "      <td>GRANDFATHER</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-9</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>68.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>FAM_0_I</td>\n",
       "      <td>MOTHER</td>\n",
       "      <td>GRANDFATHER</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>-9</td>\n",
       "      <td>12.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>66.0</td>\n",
       "      <td>68.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>FAM_0_I</td>\n",
       "      <td>SON_1</td>\n",
       "      <td>0</td>\n",
       "      <td>MOTHER</td>\n",
       "      <td>1</td>\n",
       "      <td>-9</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>66.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>FAM_0_I</td>\n",
       "      <td>FATHER_1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-9</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>46.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>58.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>66.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>FAM_0_I</td>\n",
       "      <td>DAUGHTER_1</td>\n",
       "      <td>FATHER_1</td>\n",
       "      <td>MOTHER</td>\n",
       "      <td>2</td>\n",
       "      <td>-9</td>\n",
       "      <td>14.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>46.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>58.0</td>\n",
       "      <td>66.0</td>\n",
       "      <td>66.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>FAM_0_I</td>\n",
       "      <td>FATHER_2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-9</td>\n",
       "      <td>18.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>67.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>FAM_0_I</td>\n",
       "      <td>DAUGHTER_2</td>\n",
       "      <td>FATHER_2</td>\n",
       "      <td>MOTHER</td>\n",
       "      <td>2</td>\n",
       "      <td>-9</td>\n",
       "      <td>12.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>67.0</td>\n",
       "      <td>68.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>FAM_1_I</td>\n",
       "      <td>GRANDFATHER</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-9</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>53.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>61.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>FAM_1_I</td>\n",
       "      <td>MOTHER</td>\n",
       "      <td>GRANDFATHER</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>-9</td>\n",
       "      <td>14.0</td>\n",
       "      <td>15.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>53.0</td>\n",
       "      <td>53.0</td>\n",
       "      <td>61.0</td>\n",
       "      <td>61.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>FAM_1_I</td>\n",
       "      <td>SON_1</td>\n",
       "      <td>0</td>\n",
       "      <td>MOTHER</td>\n",
       "      <td>1</td>\n",
       "      <td>-9</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>53.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>61.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       FID          IID          PAT     MAT  SEX  PHENO  M1-A1  M1-A2  M2-A1  \\\n",
       "0  FAM_0_I  GRANDFATHER            0       0    1     -9   12.0    0.0   29.0   \n",
       "1  FAM_0_I       MOTHER  GRANDFATHER       0    2     -9   12.0   16.0   27.0   \n",
       "2  FAM_0_I        SON_1            0  MOTHER    1     -9   12.0    0.0   29.0   \n",
       "3  FAM_0_I     FATHER_1            0       0    1     -9   14.0    0.0   21.0   \n",
       "4  FAM_0_I   DAUGHTER_1     FATHER_1  MOTHER    2     -9   14.0   16.0   21.0   \n",
       "5  FAM_0_I     FATHER_2            0       0    1     -9   18.0    0.0   25.0   \n",
       "6  FAM_0_I   DAUGHTER_2     FATHER_2  MOTHER    2     -9   12.0   18.0   25.0   \n",
       "7  FAM_1_I  GRANDFATHER            0       0    1     -9   15.0    0.0   26.0   \n",
       "8  FAM_1_I       MOTHER  GRANDFATHER       0    2     -9   14.0   15.0   26.0   \n",
       "9  FAM_1_I        SON_1            0  MOTHER    1     -9   14.0    0.0   29.0   \n",
       "\n",
       "   M2-A2  M3-A1  M3-A2  M4-A1  M4-A2  M5-A1  M5-A2  M6-A1  M6-A2  \n",
       "0    0.0   39.0    0.0   43.0    0.0   59.0    0.0   68.0    0.0  \n",
       "1   29.0   34.0   39.0   43.0   49.0   54.0   59.0   66.0   68.0  \n",
       "2    0.0   34.0    0.0   49.0    0.0   54.0    0.0   66.0    0.0  \n",
       "3    0.0   37.0    0.0   46.0    0.0   58.0    0.0   66.0    0.0  \n",
       "4   27.0   34.0   37.0   46.0   49.0   54.0   58.0   66.0   66.0  \n",
       "5    0.0   36.0    0.0   48.0    0.0   54.0    0.0   67.0    0.0  \n",
       "6   29.0   36.0   39.0   43.0   48.0   54.0   59.0   67.0   68.0  \n",
       "7    0.0   36.0    0.0   49.0    0.0   53.0    0.0   61.0    0.0  \n",
       "8   29.0   36.0   39.0   44.0   49.0   53.0   53.0   61.0   61.0  \n",
       "9    0.0   39.0    0.0   44.0    0.0   53.0    0.0   61.0    0.0  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pandas.read_csv(ped_path, sep='\\t').head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "-Zwk97emzgUR",
   "metadata": {
    "id": "-Zwk97emzgUR"
   },
   "source": [
    "## Dataset preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7768e352-dfde-4136-b54e-cb3c1b31661d",
   "metadata": {
    "id": "7768e352-dfde-4136-b54e-cb3c1b31661d"
   },
   "source": [
    "The initial steps of the Python module recombulator-x consist in reading the PED file and identifying the informative families for the estimation of recombination rates using the function **ped2graph**. This function takes a ped file as input and build a graph with the relationships. It returns a list of tuples, each composed by the graph, a dictionary (with iid as key and their tab row as value) and the family identifier. \n",
    "\n",
    "The pedigree file can be a *.tsv* (tab as separator value), a *.xlsx* or whatever format with a space as separator value.\n",
    "\n",
    "For recombination, informative subfamilies are either those with:\n",
    "\n",
    "* a phased mother and at least one son or phased daughter, called type I families,\n",
    "* an unphased mother and at least two between sons and phased daughters, called type II families.\n",
    "\n",
    "Notably, females can be phased when their father is available: in this way, they will be virtually transformed into males, thus being allowed to take part to informative families."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "Km7-Rkcp54vD",
   "metadata": {
    "id": "Km7-Rkcp54vD"
   },
   "source": [
    "**Important**: Genetic markers (from the 7th column on) must be provided according to their physical genomic position. Indeed, the algorithm will infer the recombination rate between M1 and M2, M2 and M3 and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cd09f894-fdf8-47f0-bb16-4bb2cd67b0ba",
   "metadata": {
    "id": "cd09f894-fdf8-47f0-bb16-4bb2cd67b0ba",
    "outputId": "4c60a42f-d8a3-4ed4-bd76-694973b6a67d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "200"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "family_graphs, marker_names = recombulatorx.ped2graph(ped_path)\n",
    "len(family_graphs)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2o8aOirR3Y68",
   "metadata": {
    "id": "2o8aOirR3Y68"
   },
   "source": [
    "The function **plot_family_graph** can be used to graphically represent the reported relationships between individuals within the same family. Note that plotting requires matplotlib to be installed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cfdc1f5b-8f85-49ab-b1af-0bfe20e3ca30",
   "metadata": {
    "id": "cfdc1f5b-8f85-49ab-b1af-0bfe20e3ca30",
    "outputId": "482a0c28-c51b-437c-97e3-4dbb7c2d0c34"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG+CAYAAAB/H2v/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABa+0lEQVR4nO3dd3gU1dvG8e9sKmmUEHpvoVfp0lUU7EhRRF+UJjYEu4Kg4E9UEARBQBALIKFJEUGKVEEsNCmht1BCJz3Z3Xn/iKwiLUCS2c3en+vikuzOzjwJMXvnOWfOMUzTNBERERHxQjarCxARERGxioKQiIiIeC0FIREREfFaCkIiIiLitRSERERExGspCImIiIjXUhASERERr6UgJCIiIl5LQUhERES8loKQiIiIeC0FIREREfFaCkIiIiLitRSERERExGspCImIiIjXUhASERERr6UgJCIiIl5LQUhERES8loKQiIiIeC0FIREREfFaCkIiIiLitRSERERExGspCImIiIjXUhASERERr+VrdQEiIiKZwuGA+CSIT4C4RIhLgDQ7OJ3pz9ts4OcLYcEQEgyhQRAcBD7qCXgzwzRN0+oiREREboppwoUEOBoLJ8+mfwxgABl5dzMMKJAPihRID0jidRSERETE85gmHD8NR45DYnLGg8+1BOeCYgWhYHh6QBKvoCAkIiKeJTEZdu5PH/rKCmEhULE05ArImvOLW1EQEhERz2CaEHMC9sUA5q13gK7FZkCZ4lAkQt2hHE5BSERE3J/DAdv2wtkL2Xvd8NxQuWz6RGvJkRSERETEvdntsGV31g2FXU/uEKhWHnx8rLm+ZClFXBERcV8OB2y1MAQBnI9Pr+HibfiSoygIiYiIezLN9EnRFywMQRedj4foA1ZXIVlAQUhERNxT7Bk4dc7qKv4RewZOnbW6CslkCkIiIuJ+UtNg9yGrq7hc9AFIS7O6CslECkIiIuJeTBN2HUifH+Ru7A73DGhy0xSERETEvVyIh9Pnra7i6k6etXbytmQqBSEREXEvMbHpW2a4K4P0vc0kR1AQEhER95Galj4h2Z1XuDOBE2fSd7YXj6cgJCIi7uPYKfcOQReZJpw4bXUVkgkUhERExH2cPGN1BRkX60G1ylUpCImIiHtwOiEhyeoqMi4+Mb0zJB7N1+oCREQkZylVqhQBAQEEBgaSkJBAlSpVeO2112jUqJHrmGbNmhETE8Pu3bsx/t7d/cC2HdRs24JzP/x8yfmM5nU5O385eUJDSbPbef/bL5m2bDE+Pjb8ff0oWbAwA/+vOzXLR7Ji4x/0GT2MTROnul5/4NhRanbrzLkffqbm048BkGq3E334INVKlwUgskRJhvZ4nrKdH3I9BhDg78+vYydz4NhR13MXo8/73XrTtk5lCM511a9F165dWbt2Lbly5SIkJIQRI0ZQt27dW/nySiZTEBIRkUw3ffp0atasCcDs2bNp06YNixcvpn79+uzevZvdu3eTN29eVq5cSfPmzdNflHj9blDXDwYRn5TEujGTyBsaBsDS338l+vBBapaPvO7rLwaki+Hov4EpNFfQJY/927+fW/DLah59723OdnkUn2sEoYceeogJEybg6+vLggULaN++PQcOHLhunZJ9FIRERCRLPfzww2zYsIGPP/6YGTNmMGnSJB5//HEKFy7MxIkT/xWEkq952/zuI4eYs2YFh6MWuEIQwB231c/aT+AKWtWpS1xiAmdijhFRtOBVj7v//vtdf2/QoAExMTHY7XZ8ffX26y70LyEiIlmufv36zJs3D4fDwVdffcXy5cvJnz8/gwYN4vz58+TOnRuc5jXvGNu4O5pyRYuTLyz3Na8VffiQawgM0ofBMiouKfGS11YpXYYpbw++7LiZK5fTsnZdIvLkyfC5R44cSZs2bRSC3Iz+NUREJMuZf08qXrhwIaVKlaJixYoA3HHHHUydOpVnnnkG4xoTjy/OI/q3vTFHaDfgVZJSU2hUpTpfvv4OAJHFS1xxjlBGXGto7GJIOhN3gVPnz7F8+Nj08JYB3377LVFRUaxatSpDx0v2URASEZEs99tvv1G1alUmTpzIrl27KFWqFABJSUkcOHCAZ555hoh84cQnJZFmt+P3d9fk1LlzBPoHEBYcTK3ykeyJOczZuAvkDQ2jbNFibJo4lck/zuf7NSuz/HO4GJJM0+S9r7+g03tvsfP2nwm8zuumT5/OoEGDWLZsGQULXn0YTayh2+dFRCRLzZ07l7Fjx9KvXz+WLVvGnj17OHDgAAcOHODYsWMcPXqUzZs3ExQcTMMq1ZiwYI7rtWPnzeSOOnUxDIPyxUrwQONmPP3he5yLi3Mdk5CcvbfcG4ZB/ye6kT93HsZ+9+01j42KiuLtt99m6dKllChRIpsqlBuhjpCIiGS6jh07um6fr1y5MgsXLmTlypXcdddd5PnXvBqbzUanTp2YOHEin/Z9nW/eGsQLn37M+PlzMIEKxUow/uW3XMdPfv0dhnw7ifq9/w9fHx/yhoQRkScvrz32RKbU/d85QgCrR0247DjDMBj2zIt0HNKfnm+/QVBQ0BXP17lzZwoVKsQDDzzgemzZsmWEh4dnSr1y6wzT1GpQIiLiBk6fg7/2WF3FjaleAfKGXf84cVsaGhMREfcQGmx1BTcu9MqdIPEcGhoTERH34O8Hfr6es6t7gD/4+vLuu+8ye/bsy56eNWsWZcuWvcILxZ1oaExERNzH9n1w6oz770BvAAXzQ2QpqyuRW6ShMRERcR9FItw/BEF6jUUirK5CMoGCkIiIuI/cIZDreivzuIGQIM+c0ySXURASERH3YRhQrIDVVVxfUQ+oUTJEQUhERNxLgXBw5/24/P0gIp/VVUgmURASERH34usDkSWtruLqIkuBj94+cwr9S4qIiPvJnxcKuGHXpVB+yJfb6iokEykIiYiIeypXIn1dIXfh7wdli1tdhWQyBSEREXFPfr5QpWz6BGqr2QyoUi592E5yFAUhERFxX7lDoWq59AUMrWIYULU8hOl2+ZxIK0uLiIj7O3sB/toNzmx+y7LZoFp5yBOavdeVbKMgJCIiniEuIX0LjuSU7LleUCBUKpO+eKLkWApCIiLiOZxOOHAUDh/P2uuUKAwlC6d3hCRHUxASERGPs3fjFkqmmPimpKXPH7rVd7KL5wgJggoltX2GF1HUFRERj+F0OnnnnXcoV7sGAY1rQ81IyJ/v1iZTG0b6StG1KkLtSgpBXkYdIRER8QgxMTHcf//9/PnnnwAYhoHT6Ux/MjUNYs/Ahfj0uUTJqdc+WWBA+l1gYSFQIC/4+WVx9eKuFIRERMTtzZgxgy5dupCS8s9EaT8/P1JTrxJ47A6ITwS7HRzO9I6RzZa+h1lIkNYDEhc3WrJTRETkcn369GHkyJGXPe7jc40w4+ujW94lQzRHSERE3NqRI0eu+LjhDitOi8dTEBIREbc2c+ZMJk2ahN9/5vHYdGu7ZAJ9F4mIiNurW7cuaWlp+Pj4kCtXLgB8fTW7Q26dgpCIiLi9vn37YhgGr7zyCvHx8XzyySdMmDDB6rIkB9BdYyIi4tbWrFlDkyZNCAoK4siRI+TNm9fqkiQHUUdIRETclmma9OnTB8MwePXVVxWCJNOpIyQiIm5r3rx5PPDAA4SGhnLkyBHCwsKsLklyGHWERETELTkcDl5++WUMw+DNN99UCJIsoY6QiIi4pcmTJ9O1a1fy5cvHoUOHCA7WHmCS+dQREhERt5OcnMwbb7yBYRi8/fbbCkGSZRSERETE7YwZM4YTJ04QERFBr169rC5HcjAFIRERcSvnz5/n3XffxTRNBgwY4FpAUSQraFlOERHJEgnJTg6edHD8nIM0u4ndkb4Xqp+vQaE8PpSM8CE48PLfxz/88EPi4+MpWrQo3bp1s6By8SaaLC0iIpkiJc3ktz2pbD2Yyv4Tds4mpL+9GIBh/P0XE0wTLr7x5A02KF3Ql2ol/albzp8zp45TunRpUlJSmDBhgoKQZDkFIRERuSXHzzlYuS2ZNdtTSE5LDz038s5y8fgAP0g9tobZ414l1D+FXbt2XbbRqkhmUxASEZGbEp/sZNrqRDbsTsVmgDMT3k2cDjs2H18ifGJ488kqhFxh6EwkMykIiYjIDdu4L5WvViSQlGJmSgD6L5sBufwNnmgRTO0y/pl/AZG/KQiJiEiGpaSZfP1zAhv2pF6c8pNlLp6/bjk/nmwRQoCfkYVXE2+lICQiIhmSmOJkxPw4Dpx03NAcoFtlGFAqwocX7w294l1mIrdCQUhERK4rKdXkozkXiDnjyJKhsOuxGVAknw+vPBhKUIDCkGQefTeJiMg1pdlNPl0QZ1kIgvSJ2EfPOPj0h3jS7Pr9XTKPgpCIiFzT3N+S2HvcblkIushpwr7jdub+lmRtIZKjKAiJiMhV7T1u56eNyVk6KfpGmMBPG5PZe9xudSmSQygIiYjIFaXZTSYti09fFdqNGAZMWqYhMskcCkIiInJFP/6ZxMnzTsuHxP7LacLJ804W/qkhMrl1CkIiInKZNIfJsi0pbjMk9l8msHxLirpCcssUhERE5DJ/7E0lMdW9Q0Ziqskfe1OtLkM8nIKQiIhcZvmWZLebG/RfhgHLtiZbXYZ4OAUhERG5xNEzDvbHZu/q0TfDNOFArIOYM7qDTG6egpCIiFxiz7E0q0u4IXuPKQjJzVMQEhGRSxw86cDm5sNiF9lscOCkw+oyxINprzERES9XqlQpAgICyJUrFwC2fNWp88gIVn7dmwObf+Cx97fjFxDMos86knj+OABnYv4ib5FKGIYPfoEh3Nf3B754NpwuH+0jICi369zf9a/JnT2+Ibx4NVZ+/SwxO1cQGJLf9XzVFj2o0LDzJc857MlElKjF7Y8Nx9c/6Kp1H41exW9z38W0J1Igtw9t27blgw8+wGbT7/iScb5WFyAiItabPn06NWvWxO4weXb8WZITL3Bo62LyFa3K/j/nUqHhY9z97HTX8V88G869L/1wSejJiOp3PE/Vlr2u+ZwjLYUfPn2QbSu/oMadL1z1XP5BeWjx1BfkiSjFx0/k4p7Wd/L111/zf//3fzdUk3g3xWYREXGJTzZxmrD399kUiWxGtZbPEL1uSrbW4OMXQKGy9Yk/c/iax+UvXp2w/KVwmuAggJo1a3LgwIHsKVJyDHWERESEjh07kitXLuwOKNywL7vWfUvttm9QtGJT1n73CudO7CZPwfLXPc+CT9piGD6ujy8OpV20Zekodq2f5vq4UYcPKFSu4SXHpCZd4NiutdR9oH+G6z8Sc4yZM2eyYMGCDL9GBBSERESEf4bGYs87eObDdSSeP0GxSi0wbDbK1WvPrnVTqPfgwOue57/DZd/1r3nJ89caGtuydBTR66ZwIXYvxaveSeEKTTJUe2rSBR7v1J5XX32V2267LUOvEblIQUhERFx8fQyi131LWko809+pDYDTkYZpmtx239vYfLLubeNiSIo/c4QFw9uyY/WXVG761DVfk5ocx6LPOvB05/vo27dvltUmOZeCkIiIuPiSxp4NM7j/5cXkKVTB9fjcD+/k8F8/UbJGmyyvISRfMRp2+IA10/pRocGj+PrnuuJxacnxLPqsA8Uqt6J//4wPo4n8myZLi4iIy6KFc8mTv/glIQigbN1HMmXS9Jalo5j9fjPXny1LRl3xuJLV7yFPwfJsXzXpquf6a8U4Th74k8NbFtCwXi1q1qzJkCFDbrlG8S5aR0hERC7x+eI4/tyX5vZbbED6fmO1y/jRq3Wo1aWIh1JHSERELlEywhcPWVgag/R6RW6WvntEROQSpQr44nSTblBS3El+HPXIZY8Xrdic+g8Pwmmm1ytys/TdIyIil6hQxJfQXAZxSdanoVyhETz85sqrPh+ay6BCEb2Vyc3T0JiIiFzCx2bQomoghpuPjxkGNK8aiI+n7BArbklBSERELtOkcoBHzBNqWjnA6hLEwykIiYjIZfIE26hVxg93bbbYDKhV2o88wXobk1uj7yAREbmidg2CsLnpu4TNBo80DLK6DMkB3PRbXERErBaR24f2jdwzbLRvFEREbp/rHyhyHQpCIiJyVc2rBlC+sK/bDJHZDChX2JfmVTU3SDKHgpCIiFyVzTDo2jKYAD/D8jBkM8DfF7q2DMbm7re0icdQEBIRkWuKyO3DS/eF4uODZbfUGwb42OCl+8MooCExyUQKQiIicl2lC/rS595Q/H3I9s6QzQA/H+hzXyhlCmrxRMlc2nRVREQy7OBJOyPmx5GYYmbLNhw2A4ICDPrcG0pJbaUhWUBBSEREbkhCspOpqxPZsDsV0+nAsGX+UJUBmEC9cv482jSIkEANYEjW0HeWiIjckOBAGw/WTGL1190xzJRMX4HaMCAo0KD33SF0vytEIUiylPqMIiJyw95//32O7VjKgHaB/HU8F8u3JnMuwcRmcFNDZhdflyfYoGW1QJpWDiBYAUiygYbGRETkhhw4cIDIyEgGDBjAW2+9BYDTabL1UBo/b01mxxG7KwxdLRj9+3GbARWL+dKyWiDVSvhhs/o+ffEqCkIiInJDnnjiCZYsWcKePXsIDg6+7Hm7wyTmjIODJ+0cjHVw9IyDVLtJmsPEz8fA39egSD4fSkb4UDLClyLhPvj5KPyINRSEREQkw7Zs2ULNmjX57LPPeOaZZ6wuR+SWKQiJiEiG3XvvvURHR7N9+3b8/PysLkfklmmytIiIZMjq1av54Ycf+O677xSCJMdQR0hERK7LNE0aN25MamoqGzZswGbTHV2SM6gjJCIi1zV37lzWrVvHkiVLFIIkR1FHSERErslut1O9enWKFi3KkiVLrC5HJFOpIyQiItf09ddfs2PHDr755hurSxHJdOoIiYjIVSUlJVGhQgUaNWrE9OnTrS5HJNNpoFdERK7qs88+4/jx4wwePNjqUkSyhDpCIiJyRefOnaNMmTJ06tSJMWPGWF2OSJZQR0hERK5o6NChpKSkMGDAAKtLEckyCkIiInKZo0ePMnLkSF566SUKFSpkdTkiWUZDYyIicpmePXsya9Ys9u7dS+7cua0uRyTL6PZ5ERG5RHR0NBMnTuTDDz9UCJIcTx0hERG5RPv27dmwYQPR0dEEBgZaXY5IllJHSEREXDZs2MDMmTOZPHmyQpB4BXWEREQESN9YtVWrVpw8eZJNmzbh4+NjdUkiWU4dIRERAeCnn37i559/Zt68eQpB4jXUERIREZxOJ3Xq1CE4OJjVq1djGIbVJYlkC3WERESE6dOns2nTJtasWaMQJF5FHSERES+XmppKpUqVqFq1KnPnzrW6HJFspY6QiIiXGz9+PAcOHGDevHlWlyKS7dQREhHxYvHx8ZQtW5Y2bdrw5ZdfWl2OSLbTXmMiIl5s+PDhnD9/nkGDBlldiogl1BESEfFSJ0+epEyZMvTo0YNhw4ZZXY6IJdQREhHxUkOGDMFms/Hmm29aXYqIZRSERES80P79+xkzZgyvvfYa4eHhVpcjYhkNjYmIeKEuXbqwdOlS9uzZQ3BwsNXliFhGt8+LiHiZLVu2MGXKFMaMGaMQJF5PHSERES/Ttm1bdu/ezbZt2/Dz87O6HBFLqSMkIuJFVq1axcKFC5k+fbpCkAjqCImIeA3TNGnUqBFpaWls2LABm033y4ioIyQi4iW+//571q9fz9KlSxWCRP6mjpCIiBew2+1Uq1aN4sWL89NPP1ldjojbUEdIRMQLfPXVV+zcuZMpU6ZYXYqIW1FHSEQkh0tKSqJ8+fLcfvvtfPfdd1aXI+JWNEgsIpLDjR49mhMnTjB48GCrSxFxO+oIiYjkYGfPnqVs2bI8+uijfPbZZ1aXI+J21BESEcnBhg4dSkpKCv3797e6FBG3pCAkIpJDxcTEMHLkSPr27UuhQoWsLkfELWloTEQkh+rRowezZ89m3759hIWFWV2OiFvS7fMiIjlQdHQ0kyZN4qOPPlIIErkGdYRERHKgRx55hN9++43o6GgCAwOtLkfEbakjJCKSw/z666/MmjWLr776SiFI5DrUERIRyUFM06Rly5acOnWKTZs24ePjY3VJIm5NHSERkRxk8eLFrFixgvnz5ysEiWSAOkIiIjmE0+mkdu3ahIaGsmrVKgzDsLokEbenjpCISA7x3XffsXnzZtauXasQJJJB6giJiOQAqampVKpUiapVqzJ37lyryxHxGOoIiYjkAOPHj+fAgQPMmzfP6lJEPIo6QiIiHi4uLo6yZcty7733MmnSJKvLEfEo2mtMRMTDDR8+nAsXLjBo0CCrSxHxOApCIiIeLDY2lo8//pjnnnuO4sWLW12OiMdREBIR8WBDhgzBZrPxxhtvWF2KiEdSEBIR8VD79+9n7NixvPbaa4SHh1tdjohH0mRpEREP1aVLF5YtW8bu3bsJDg62uhwRj6Tb50VEPNDmzZuZMmUKY8eOVQgSuQXqCImIeKA2bdqwZ88etm3bhp+fn9XliHgsdYRERDzMypUr+fHHH4mKilIIErlF6giJiHgQ0zRp2LAhDoeDDRs2aE8xkVukjpCIiAf5/vvv+fXXX1m6dKlCkEgmUEdIRMRD2O12qlWrRokSJVi8eLHV5YjkCOoIiYh4iMmTJ7Nz506mTp1qdSkiOYY6QiIiHiApKYny5cvTpEkTpk2bZnU5IjmGVpYWEfEAo0aN4sSJE7z33ntWlyKSo6gjJCLi5s6ePUuZMmV47LHH+Oyzz6wuRyRHUUdIRMTNDR06lNTUVPr37291KSI5joKQiIgbi4mJYeTIkfTr149ChQpZXY5IjqOhMRERN9a9e3e+//579u7dS1hYmNXliOQ4un1eRMRN7dy5k0mTJjFs2DCFIJEsoo6QiIibateuHX/88QfR0dEEBARYXY5IjqSOkIiIG/r111+ZPXs2X331lUKQSBZSR0hExM2YpknLli05deoUmzZtwsfHx+qSRHIsdYRERNzM4sWLWbFiBQsWLFAIEsli6giJiLgRp9NJ7dq1CQsLY+XKldphXiSLqSMkIuJGpk2bxubNm/nll18UgkSygTpCIiJuIjU1lYoVK1K9enW+//57q8sR8QrqCImIuIlx48Zx8OBBFixYYHUpIl5DHSERETcQFxdH2bJluffee5k0aZLV5Yh4De01JiLiBoYPH86FCxcYNGiQ1aWIeBUFIRERi8XGxvLxxx/z/PPPU7x4cavLEfEqCkIiIhYbPHgwPj4+vPHGG1aXIuJ1FIRERCy0b98+Pv/8c1577TXy5ctndTkiXkeTpUVELPT444+zfPly9uzZQ1BQkNXliHgd3T4vImKRzZs3M3XqVMaOHasQJGIRdYRERCzSpk0b9u7dy19//YWfn5/V5Yh4JXWEREQssGLFCn788UdmzJihECRiIXWERESymWmaNGjQANM0+fXXX7WnmIiF1BESEclmc+bMYcOGDSxbtkwhSMRi6giJiGQju91O1apVKVmyJIsXL7a6HBGvp46QiEg2mjx5MtHR0UybNs3qUkQEdYRERLJNYmIiFSpUoGnTpkydOtXqckQErSwtIpJtRo0axYkTJ3jvvfesLkVE/qaOkIhINjh79ixlypShc+fOjB492upyRORv6giJiGSDDz74gLS0NPr37291KSLyLwpCIiJZ7MiRI3z66af07duXggULWl2OiPyLhsZERLJY9+7d+f7779m7dy9hYWFWlyMi/6Lb50VEstDOnTuZNGkSw4cPVwgScUPqCImIZKGHH36YjRs3snPnTgICAqwuR0T+Qx0hEZEssn79eubMmcPXX3+tECTiptQREhHJAqZp0qJFC86cOcPGjRvx8fGxuiQRuQJ1hEREssCiRYtYuXIlCxYsUAgScWPqCImIZDKn00nt2rUJCwtj5cqV2mFexI2pIyQiksmmTZvG5s2b+eWXXxSCRNycOkIiIpkoJSWFihUrUrNmTebMmWN1OSJyHeoIiYhkonHjxnHo0CEWLlxodSkikgHqCImIZJK4uDjKli3Lfffdx8SJE60uR0QyQHuNiYhkkmHDhnHhwgUGDhxodSkikkEKQiIimSA2NpZhw4bx/PPPU7x4cavLEZEMUhASEckEgwcPxsfHhzfeeMPqUkTkBigIiYjcon379vH555/z+uuvky9fPqvLEZEboMnSIiK3qHPnzqxYsYLdu3cTFBRkdTkicgN0+7yIyC3YtGkTU6dOZdy4cQpBIh5IHSERkVtwzz33sG/fPrZt24avr363FPE0+r9WROQmrVixgkWLFjFjxgyFIBEPpY6QiMhNME2TBg0aALB+/XrtKSbiofQrjIjITZg9ezYbNmxg+fLlCkEiHkwdIRGRG2S326latSqlSpVi0aJFVpcjIrdAHSERkRv05ZdfEh0dzbRp06wuRURukTpCIiI3IDExkfLly9OsWTOmTp1qdTkicou0srSIyA0YNWoUsbGxvPfee1aXIiKZQB0hEfEOpgkpqRCXCPEJEJ8EaXZwOsEwwGaDQH8IDYKQYAgJAl+fS05x5swZypYty+OPP86oUaMs+kREJDNpjpCI5GyJSXD0JJw4DXZH+mMGcKVfAS8AsWf++Tg4FxSJgALh4OvDBx98QFpaGm+//XY2FC4i2UFBSERyHtOEU+cg5gScj7/C8xk8T0IS7D4Ee48QHxLI4jlz6devHwULFszMakXEQhoaE5GcJSkFovdfOQDdAqdp4nQ6sRctQGD5kulDaSLi8RSERCRnMM30IbB9h9P/npU/2YICoVKZ9HlEIuLRFIRExPPZHbBtD5yLy97rli8BRQpk7zVFJFMpCImIZ0uzw5ZdEJ9ozfVLFYWSha25tojcMg1yi4jnsjusDUEAB2Lg8HHrri8it0RBSEQ8k2mmD4dZGYIu2nck/fZ8EfE4CkIi4pliYrN/TtC17DoIyalWVyEiN0hBSEQ8T2JyehfGnTidsOtAeqdKRDyGgpCIeBbTTF8nyB0Dx9kLcPyU1VWIyA1QEBIRz3LyLFxIsLqKq9t7GBwOq6sQkQxSEBIRzxJzwuoKrs3hTA9rIuIRFIRExHPEJ7p3N+iiI24e1kTERUFIRDzH0ZPpO8e7u4QkuJC5e52JSNZQEBIRz2CaEHs6a/cQyywGEHvG6ipEJAMUhETEMySlpM+/8QQmnjGEJyL4Wl2AiEhGlIqsQGJ8PDEzF+Lnm/6j6+eNv9PypWd4sV0nRjzfjx0H9/Py2JHsPHQAgEolSvPxMy9SsWQpBkz6nHlrVwGwJ+YIEXnykjs4GIDp77xPz2H/o88jj/Jgk+aua/7f/wZSs1wF+rR/jMk/zufF0cMoXaiI6/nG1WrwWZ/XLnnO7nAQkScv415+kwq1KoLhCWN5It5LQUhEPINpUqJgIeatXUW7Zi0BmPjDXG6LrATA0VMnafZiT0Y815fH7rgbgGnLFtO8T082fTGFd5/qxbtP9QKg+Ys9Lws9GdGi5m18P+Tj6z738pgR9Bk1jIUP3gPBuW7msxWRbKKhMRHxDE4nXe++l0k/zgPgfHw867f/xd31GgIw5vuZNK9ZxxWCAB5t1ZqmNWozZu7MbC21VZ26HDxxzD32QRORa1JHSEQ8RuOqNRgzdyZHT51k3tpVtG/eCh+bDwB/7t7JnbfVv+w1DatUY+kfGzJ0/pc++4SBk8e7Pj4Ue4Ka5Sq4Pv550+/UfPox18dP3n0vL7V/7JJzOJ1O5qxeQaeWd4FdCyuKuDsFIRHxDH/fLdblzjZMXrSA79esYMrb7zFlyaJMu8Qnz7502Ryhf7vW0NjFkHQo9gT5QsP4deyX6fuPiYhb09CYiHiGv+ccP9G6LZ/O+o5A/wDKFyvherp2+Yqs27b1spet27aV2uUjs7y8FjVvY9PEqRyZ8QMVipeg94ihnrHmkYiXUxASEY9SJH8E/+vxLEN7PnfJ48880I6fN/7O1KX/dIimLVvMik1/0PvB9tlWX1BgIF+88jYL1//Cxu3bsu26InJzNDQmIp7hX7ehd73n/sueLhpRgBUjPuflsSN5e+JYDMMgsnhJVo4cR+Hw/JlSwn/nCEWWKMn0d/532XFF8kfwcsfHGTDsI+bf1SpTri0iWcMwTdMT1mkVEW934CgcOuoZK0tfVL86BPpbXYWIXIOGxkTEM4QGeVYI8vWBAD+rqxCR61AQEhHPEBpsdQU3JjRYq0qLeAAFIRFxe6dPn2biN18Te/4cHjOanzvE6gpEJAMUhETELZ09e5Yvv/ySe+65h0KFCtGjRw/m/vaL59ySXihzJmiLSNbSZGkRcRvnzp1j7ty5REVFsWTJEux2O02bNqVDhw48/PDDFMoXDus3u/dcIQPIlweqlrO6EhHJAN0+LyKWOn/+PPPmzSMqKorFixdjt9u5/fbbGT58OO3ataNw4cKXviAiH8SesabYjDCBogWsrkJEMkhBSESy3YULF5g/fz5RUVEsWrSI1NRUGjduzMcff0y7du0oWrTo1V9ctKB7B6GgQMgTanUVIpJBGhoTkWwRFxfHggULiIqK4scffyQlJYWGDRvSoUMHHnnkEYoVK5bxk+05BDGxWVfsTXI4nczev4N2/9cFm01TMEU8gYKQ3BLTNDl1wcnBkw4OnbKTmGKS5jAxAD8fg9zBNkrk96FkAV9yB+mNwdvEx8dfEn6Sk5OpX7++K/yUKFHi+ie5EocDftsGKamZW/AtMIFF2zbS5tkeNGrUiAkTJlC5cmWryxKR61AQkhuWZjf5Y28qv0SnsP+EneS09Md9/s45F7+jDCP9786/Pw7NZVC+sC9NKgdQubgfNq2xkiMlJCSwcOFCpk+fzsKFC0lKSqJu3bqu8FOqVKnMudC5ONgcnTnnygy5AuG2yqxas4bu3buzf/9+3nzzTd544w0CAgKsrk5ErkJBSDLs5AUHq7alsGpbComppivo3AibkR6MwkNttKwWQOOKAQQHqlPk6RITE/nxxx+JiopiwYIFJCYmUqdOHTp06ED79u0pXbp01lz44NH0rTes5mODmhUhJAiA5ORkhgwZwgcffED58uWZMGECjRs3trhIEbkSBSG5rjSHyfzfkli0MRmDfzo8t8oA/Hyh4+1BNKkUgKEOkUdJSkpi0aJFREVFMX/+fBISEqhVq5Yr/JQtWzbrizBN2HvY2vlCNgOqR15xAcWtW7fSvXt3fv31V5555hn+97//kTt3bguKFJGrURCSa9p/ws7EZfHEnnNm6dItlYr58mSLYMJDfbLwKnKrkpOTLwk/8fHx1KhRwxV+ypcvn/1FmSbsj4HDx7P/2j42qFbhmqtIOxwOxowZw5tvvklYWBhjxozhgQceyMYiReRaFITkikzT5Ic/kpm3IQnDyLwu0NXYjPQ9Kru2DOG2ctqt252kpKSwePFioqKimDdvHnFxcVSrVs0VfiIjI60uMd2RE7DvyI2P196swACoUtY1HHY9hw4donfv3vzwww+0a9eOUaNGXb5GkohkOwUhuYzTNJm2KpEV21IsuX6XZkE0rRJoybUlXUpKCkuWLCEqKoq5c+dy4cIFqlSp4go/lSpVsrrEK0tMgp37IS4xa69TvBCUKgI3eIu8aZpERUXxwgsvkJKSwocffki3bt10q72IhRSE5BKmaTJlVSIrLQpBF3VpHkTTygpD2Sk1NZWlS5cSFRXF999/z/nz56lUqRIdOnSgQ4cOnnMruGmmd4f2x2Rad8jk7y3OcgVAxdIQdmsbqp45c4ZXXnmFSZMm0bRpU8aPH+8+nTURL6MgJJeYuyGRBb8nW10GAM+0DqF2WQ2TZaW0tDSWLVtGVFQUc+bM4dy5c0RGRrrCT5UqVTx3EnuaHU6cgiOxN7XekGmCk/SdqbfGB7EhMTe28DCaVg6kbCHfTPm6LF++nB49enD48GH69+/Pq6++ir+/vudFspOCkLjsOprGR9/HWV2GS6AfvPdYHvIEa9ggM6WlpfHzzz+7ws+ZM2coX768K/xUq1bNc8PPlZgmnL1A4pHT2M/GE2KmYjP+CTqubg/g8/dfUpwGB5MC2JUYyJqzoZxO8wP+Wf6hSD4braoFUq9CAIF+t/a1SkpK4t133+Wjjz6iUqVKTJgwgQYNGtzSOUUk4xSEBICUNJN3vjvP2Xhnlk+MziibAVVK+PF8m5Cc9cZsAbvdzooVK4iKimL27NmcPn2asmXLusJPjRo1cuzXODnVZOYvCazcnh6AfHFSLDCV4oGp5PJx4meYOE1INQ3Op/lyMNmf2FQ/TK7+9TBID1AhgQZPNA+mVplb7+Js2rSJ7t2788cff/Dcc88xZMgQQkO1Z5lIVlMQEgC+W5PA8q0p2XbDzY14qlUwDSO1Mu+NstvtrFq1iqioKGbNmsWpU6coU6aMa8JzrVq1cmz4uWjHkTQmLYvnfKKZJd/bFwNRvXL+PNo0iJBbXBzUbrczatQo3n77bcLDwxkzZgz33ntvptQqIlemICQcPmXn3agLVpdxVYF+8NGTeQn0z9lv2pnB4XCwevVqV/iJjY2lVKlSrs5P7dq1c3z4AXA6TaavTWT51pSbWgH9RtkMyBVg8EzrECKL+t3y+Q4cOECvXr1YvHgxHTt2ZOTIkRQsWDATKhWR/1IQEr7+OZ61O1PdZkjsSh5vFkQz3VJ/RQ6Hg7Vr1xIVFcXMmTM5ceIEJUqUcIWf2267zSvCz0V2h8mEJfH8uS8tW69rGOmBqFfrEGqWvvWhMtM0mTp1Kn369MHhcPDxxx/TtWtXr/q3FMkOCkJeLjHFSb/J57A7rK7k6gygUF4bgzrl1pvA35xOJ7/88osr/Bw7doxixYq5wk+9evW88mvlcJp8vjiezfvTsnQl9GuxGfDsPSFUL5U5d3+dOnWKfv368fXXX9OyZUvGjRtHuXLlMuXcIpJ+Z6h4sXXRqW4dgiB9Dsaxs072HrdbXYqlLoafPn36ULx4cZo0acLs2bPp2LEjv/zyCwcPHmTYsGHUr1/fK0MQwDcrEthkYQiC9GG4MYvi2XMsczpS+fPn56uvvmLx4sXs37+fatWq8cEHH5CWlr0dL5GcSh0hL/fOd+c5esbNkxDpv2U3qhjAky2CrS4lW5mmya+//kpUVBQzZszgyJEjFC5cmPbt29OhQwcaNmyoVYn/9tueFMb/lGB1GUD6MFmeIIN3H8tzy7fX/1tCQgIDBw5k+PDhVKtWjQkTJlC3bt1MO7+IN9JPUC+WkmZyzANCEKSv3bLrqHf8BmyaJhs2bODll1+mVKlSNGzYkGnTpvHggw+yatUqjhw5wsiRI2ncuLFC0N8uJDr5dmUWb6txA0wTziWazPolc2sKDg7mo48+YsOGDdhsNho0aEDfvn1JSHCPACjiifRT1IsdPmW3dAjhRp087yQlzZMqzjjTNPn999959dVXKV26NPXr1+ebb77h3nvvZcWKFRw5coRRo0bRpEkThZ8rmLIqgeRU9/reME1YsS2FHUcyP8DXqVOHDRs28MEHHzB27FiqVKnCokWLMv06It5AQ2NeplSpUgQEBJArVy7ik018w6vTpPNIVn7dmwObf+Cx97fjFxDMos86knj+OABnYv4ib5FKGIYPfoEh3Nf3B754NpwuH+0jICi369zf9a/JnT2+Ibx4NVZ+/SwxO1cQGJLf9XzVFj2o0LDzJc857MlElKjF7Y8Nx9f/6rt4x50+xMpvniP++FbKlinNpk2bsuxrlF1M02Tjxo1ERUURFRXF/v37iYiIoF27dnTo0IGmTZvi4+NjdZlu769DqYxcEG91GVdkGJAvxMb7nXNjs2XNvK29e/fSs2dPli1bRufOnfnkk0+IiIjIkmuJ5ES+Vhcg2W/69OnUrFmTScvi+XVXKsmJFzi0dTH5ilZl/59zqdDwMe5+drrr+C+eDefel364JPRkRPU7nqdqy17XfM6RlsIPnz7ItpVfUOPOF656Lr/AUOre+yb1y6Qw44tBN1SHOzFNk82bN7vCz969ewkPD3eFn2bNmuHrq/8tb8TSzcmurS/cjWnC6Tgnfx1Ky7S7yP6rbNmyLFmyhK+//pq+ffuyaNEihg8fTpcuXbx20rzIjVCP3YvFJaVvp7H399kUiWxGtZbPEL1uSrbW4OMXQKGy9Yk/c/iaxwUG56VIhQaYPrmyqbLMczH8vP3220RGRlKrVi3GjRtHixYt+Omnnzh27Bjjxo2jVatWCkE3KPa8g22H7W4Zgi6yGbB8a9ZuZGwYBk8++SQ7duzgrrvu4sknn6R169bs27cvS68rkhPop64X6tixI7ly5SL2vJMqd77MrnXfUrvtGxSt2JS1373CuRO7yVOw/HXPs+CTthjGP0M3F4fSLtqydBS71k9zfdyowwcUKtfwkmNSky5wbNda6j7QP0O1OzxjbjemafLXX3+5Oj+7du0ib968PPTQQ4waNYqWLVvi53frKxB7u1XbUty2G3SR04Rth+2cPO8gInfWDnUWKFCAqVOn0qVLF3r16kXVqlV599136dOnj0K2yFXo/wwvdHFo7KM5F1j/+xYSz5+gWKUWGDYb5eq1Z9e6KdR7cOB1z/Pf4bLv+te85PlrDY1tWTqK6HVTuBC7l+JV76RwhSYZqt2N3+8A2LZtmyv87Ny5k9y5c/PQQw8xYsQIWrVqhb9/1gyPeCPTNFm9I8WtQ9BFNgPWRadwf72rz4PLTPfccw/btm2jf//+vPrqq0ybNo0JEyZQu3btbLm+iCfR0JgX8/eD6HXfkpYSz/R3avNd/5rs/X0Wu3+NwunI2sULq9/xPO3eWk37dzZw6uAmdqz+MkOv83XDucM7duzg3XffpWrVqlStWpURI0ZQr149FixYQGxsLF9++SX33HOPQlAmO3nBSWKKB6Qg0rtC2b0gaEhICJ988gnr168nLS2NevXq8eqrr5KY6D7LDIi4A3WEvFguXzt7Nszg/pcXk6dQBdfjcz+8k8N//UTJGm2yvIaQfMVo2OED1kzrR4UGj+Lrf/U5QE4n5PJ3j+weHR3NjBkziIqKYuvWrYSGhvLAAw/wv//9j7vuuouAgACrS3RraWlpvP/++0ybNg0fHx/8/f0pWbIkAwcO5Ny5c9xzzz1ERkbidDrJlSsXI0eOpEGDBq7X79u3j3LlylG77evUuudl1+O71k1l1bfP0+yJMZSv3xGAQ1sXs2XZZ9zbZx6QPvk/b5FKgIHDnkL+4jWodU8/8hauCMAfPwxlx6qJBOUu7Dpv697fEZynMEd3rWbhyAdd59+x+kt2rJ4MQMLZI/j4BbrulGzwyGB2r/+O8GLVqNqyF/tjHZimyaBBgzh37hwjRoxgxYoVrs/1otKlSzNnzpxLnrva1yEj6tWrxx9//MHHH3/MoEGDmDVrFuPGjeOOO+64ofOI5FQKQl7s8JaFhOYrfkkIAihb9xGi10255SD03zlC5eo+QvU7n7/suJLV7+Gv5Z+zfdUkqt/x7BXPZU9NJGpQPQJsacTHnadYsWJ06dKF//3vf7dU443YvXu3a9hry5YthISEcP/99/Pee+/RunVrAgO1KWxGde3alfj4eNatW0fevHkBWLp0KdHR0RQsWJDIyEjXEgmjR4/mqaeeYvv27a7XT5o0iap1mrFr3VRq3t3vkrujQvIV548fPqBMnYfw8b1yF+7isK7pdLJz7VfMH9aGh17/mdD8JYH0/wcaPvL+Za+L/uVbikQ2JXrdFMrX70ilJl2p1KQrACu/ftYVei7avf4719+TUk3OxDsvO+e/P9drPXelr0NG+fn58cYbb9CuXTt69uzJnXfeyZNPPsmwYcMIDw+/4fOJ5CQKQl7mwIEDrr93ebwjx0MuDztVW/Skaouero+7fXb6smOu9Fin9za5/t7sic+uWsOVnmvbZ+5Vjwfw9Q/isSF/Map73kzdsuB69uzZ4+r8bNq0ieDgYO677z4GDhzI3XffTa5cnncXm9V2797NnDlzOHz4sCsEAa4OxYoVKy45vlWrVrz22muujx0OB5MnT+bR12cz6YOnOBq9iqIVm7meDy9WFcPmw/ZVE6nW8plr1mLYbFRq0pVju9eyfdUk6j989aUZUhLPc/ivJTzSfx2z32/K+dh95C5Q5kY+dQ6evPnZ/v/9OtyMChUqsHz5ciZNmsTLL7/MwoULGTFiBI8++qhutRev5R7jDGKJEhGelYMjwmzZEoL27dvH0KFDqVOnDuXLl2fw4MFUqFCBmTNnEhsby7Rp03jooYcUgm7Sxo0bKVeuHPny5cvQ8TNnzqRTp06ujxcvXkyxYsUIDK9AhYad2XWFJR9uu/9ttvz0KalJFzJ0jYhSdTh7bKfr472/zWT2+82Y/X4zVn7zXPpjv8+kWOWWBOUu6LqpICO2LB3lOlfHNnX5/PPPL3k+OjqamjVruv688sorVzzPf78ON8swDJ5++ml27NhB8+bN6dy5M23btuXgwYO3fG4RT+RZ74SSqQL9DArlsXH83OXteiskxZ3kx1GPXPZ40YrNadhuEOUKZ92364EDB1ydn99//51cuXLRtm1bXn/9ddq0aUNwsHdt9pqd9u7dS7t27UhKSqJRo0Y8+eSTrnBw/Phx7HY7v/76q+v4iRMn8tRTT7HPDuXqtufPHz4gJfEcAUF5XMfkKVieEtXvZvOSTylYOgObkv5ngf0rDY1F/zLFtcxDhYadWTT6Eerc9yY227Vn8F+8e9IwoGPjIFbP/oBz5865nr/W0Ni1vg63qlChQkRFRTF//nx69+5NlSpVGDx4MM8//7xWNBevoiDk5ZpUDmDmL0lucVt6rtAIHn5z5RWfc5rQuGLmTkA+dOiQK/xs2LCBwMBA2rRpw8svv0zbtm0JCQnJ1OtJulq1arFnzx7Onj1L3rx5KVu2LJs2bWLy5Ml8//33wD/hIC0tjd69e9O5c2fWrVvHqVOn+OGHH9iwYQMXkgfjdILTkT7pv0rz7pdcp3abV5nzv2aEhpe4bk0nD278ewL1lZ0+vJUzMdtYM6VP+r4ZQHLCGY5sW0qJaq0z/Lnf6IZGV/s6ZOYw1n333UezZs1466236Nu3L1OnTuWLL76gevXqmXYNEXemoTEv16hiAJ6wh2fB3DYqFLn13H748GE++eQTGjZsSMmSJXnrrbcoUqQIU6dOJTY2llmzZtGxY0eFoCxUvnx5HnjgAZ5++ulLOiNX2kHdz8+PkSNHcuTIEb7//nu+/vprHnzwQQ4fPswLn26h03ubaNXtyysOUwXnKUxko8fZtPiTq9aSPln6a45sX0blvyc9X0n0um+p1qo3nQZvptN7m+j03iYatBtC9LpvM/x5myb43+S38H+/DpktLCyMUaNGsXbtWhISEqhTpw5vvvkmSUlJmX4tEXejjpCXCwm0Ub+8P+t3pbr1wnQtqwfe9G/BMTExzJw5k6ioKH755Rf8/f1p3bo13377Lffddx9hYWGZXK1cz+TJkxkyZAj169fH19eXvHnzEhERwWuvvUZKSsolxwYFBTFkyBAGDhxIamoqH374IQBhgTaO46Roxeas/OY5Th3afNl1atzVh51rv77s8QWftOWf2+erc1+/ha47xv7LnpbM3t9m0rbP/EseL1P7QTbMeYfEC7EEhRXI0OcdHHj5bx0Xh78uCg0NZfXq1Zcd9++vw4MPPpglk5sbNmzIxo0bGTp0KIMHD2bGjBmMHz+eFi1aZPq1RNyFdp8Xdh1J5KN5WbsX0q0I8IUPn8xDUEDGW1dHjx5l1qxZREVFsWbNGvz8/GjdujUdOnTg/vvvJ3fuG9tAVtxP1NpElm9NxuEeU9wy5P3OubN8m43MsnPnTrp3786aNWt4+umn+eijjy65y08kp/CAQRHJKqZpMmfOHO5uWpWdaybjrhtYdLw9KEMh6Pjx43z22Wc0a9aMYsWK0bdvX8LCwpg8eTKxsbHMnz+fLl26KATlECUL+HhUCAr0g/xhnvMjt2LFiqxcuZLPP/+cGTNmUKlSJaKiotDvzpLTeM7/lZKptm3bxp133snDDz9MZGQkn73ZljzBNtxpKRGbAZWK+XJ7patPkj5x4gRjx46lRYsWFClShD59+hAcHMzEiROJjY3lhx9+4MknnyRPnjzZV7hki5IetvxDqQK+HrdWj81mo2fPnuzYsYPGjRvTsWNHHnjgAQ4fPmx1aSKZRkHIy5w9e5YXX3yRGjVqcPDgQebPn8/ChQupXjWSp1qF3PBdLVnJ1weebBF82ZvHyZMnGTduHK1ataJIkSI8//zzBAQE8MUXX3DixAkWLlxI165d1cbP4QrkthHoZ3UVGWMzoExBzwpu/1akSBFmzZrF7Nmz+f3336lcuTKjR4/G4bj5BSJF3IWCkJdwOByMGzeO8uXLM2nSJN5//33++usv7r33XlfQqFTMj9Y13WebiCdaBBMemj6f4tSpU0yYMIE777yTwoUL8+yzz+Lj48O4ceM4fvw4ixYt4qmnnsrwIn3i+WyGQeNKAdg8oMniNKFhpOfvP/fQQw+xfft2OnfuzPPPP8/tt9/Otm3brC5L5JYoCHmB1atXc9ttt9GrVy/atm3Lrl27ePXVV6+4MejDDXPRMNL6XdI73R5EuXzxTJw4kdatW1OoUCF69eqFaZqMGTOGY8eO8dNPP9GtWzfy589vdblikeZVAt36bkdI7wZFFvWlUF7PmCR9PXny5OHzzz9n1apVnDt3jlq1ajFgwIDL7vYT8RQKQjnY4cOH6dSpE02bNsXPz4/169fz1VdfUbhw4au+xmYYPNkimEYVrQtDpfy28vFrD1GoUCG6d+9OWloao0eP5tixYyxdupQePXoQERFhWX3iPgrl9SGyqK9bd4WcJrSq5j6d1szSpEkTNm3axBtvvMEHH3xAjRo1rnjbv4i7UxDKgZKSknjvvfeIjIxkxYoVfPnll6xfv5769etn6PU+tvQwdE+t9B/e2fMm48R0pLF6yvP079mCpKQkRo4cydGjR1m+fDm9evWiQIGMrdUi3uWO6u7bFTKAPMEG1Ut5yGSmGxQQEMCgQYPYuHEj+fLlo2nTpvTq1Yvz589bXZpIhmkdoRzENE1mz57Nyy+/TExMDH369OHtt9++pQUDdx1NY9KyBM7EO7N0IvWJvb+SsONLHmrTlHbt2l2zayXyb6Zp8ukP8Ww/nOaWgeiFtiFUK2n9cHNWczqdjB07ltdff53Q0FBGjx7Nww8/bHVZItelIJRDbN26lRdffJGff/6ZNm3a8Mknn1ChQoVMOXdKmsmc9Yks25qCzSDT3mycTgeYDoqyme4PVaV4saKZc2LxOmfjnfSfdo6UNKsr+YfNgPoV/HmqlXdt13L48GGeffZZ5s+fz0MPPcTo0aMpUqSI1WWJXJWCkIc7c+YM77zzDmPGjKFcuXJ88skntGnTJkuudeyMgxXbklmzI4VUe3rb/0a/eZwOOzYfX3zMJJpU9KFtgwLkCdYIrdy6X3am8OXyy/crs4IBhOYyeO+x3De0InpOYZomM2fO5PnnnycpKYmhQ4fSo0cPbJ6wsaF4HQUhD+VwOBg/fjz9+/cnNTWVAQMG8MILL+Dvn/Ut+OQ0kw27Uli7M4VDpxzY/15KxMeWvrGkaYKJiel0Yjqd2HzT50fYzDRKFzC4q04eapTyw8edZ7iKxzFNk/E/xfPH3jTL10i3GdDnvlAqFcuZc4My6uzZs7zyyitMnDiR22+/nQkTJlCxYkWryxK5hIKQB1q5ciUvvPACW7ZsoWvXrrz//vsUKlTIklocTpPjZx3sOpLImj/3sWdfDMdiT2FPSyNv7hDKlSpE83rluK1yYfKH2TxuZV3xLGkOk9E/xLEjxm7p4qA97gqmbjnPXzcos6xYsYIePXpw8OBB3nrrLV5//fVs+aVNJCMUhDzIoUOHeOWVV4iKiqJBgwZ8+umn1K1b17J6EhMT+fHHH4mKimLBggUkJiZSp04dOnTowCOPPEKZMmUsq028V0qaydhFcWw/bM/WzpABGAY8dUcw9csrBP1XUlISgwcP5sMPPyQyMpIJEybQsGFDq8sSURDyBElJSXz44YcMHTqU3LlzM3ToUB5//HFLxtuTkpJYtGgRUVFRzJ8/n4SEBGrVqkWHDh1o3749ZcuWzfaaRP7L7jD5ekUC66JTb2ou242yGeDvC93vDKF6KXU6rmXLli1069aN33//nd69e/P+++/f0p2tIrdKQciNmabJrFmz6NevH8ePH+ell17irbfeIjQ0NFvrSE5OZvHixURFRTFv3jzi4+OpUaOGK/yUL18+W+sRyaiN+1L5akUCSSlmlt5aX7WEH0+2CNbE/wxyOByMHj2at956izx58jBmzBjuv/9+q8sSL6Ug5Ka2bNnCiy++yIoVK7jvvvsYPnw45cqVy7brp6Sk8NNPPxEVFcXcuXOJi4ujWrVqrvATGRmZbbWI3Ir4ZCfTViWyYU9qpi7/YAABfvBok2AaRvpr/ttNOHjwIM888ww//vgj7du359NPP7VsvqN4LwUhN3P69GkGDBjA559/Tvny5RkxYgR33313tlw7NTWVJUuWEBUVxffff8+FCxeoUqWKK/xUqlQpW+oQyQr7T9j5+a9kNuxOxem8ueGyi0EqPNRGy2oBNKoYQEigukC3wjRNvvvuO1588UXS0tL46KOPePrppxUsJdsoCLkJu93OuHHjGDBgAHa7nXfeeYfnnnsuy++sSE1NZdmyZa7wc+7cOSpWrEjHjh1p3749VapUydLri2S3+GQna3eksH5XKkfPOFwdoovLPzjN9EnPF9+GLz4f5G9QoYgvzaoGULm4Hza9UWeq06dP069fP7766iuaN2/O+PHjNewu2UJByA38/PPPvPjii/z111889dRTDBkyhIIFC2bZ9dLS0li+fDlRUVHMmTOHs2fPUqFCBTp27EiHDh2oUqWKfhsTr5DmMIk57eDgSTtHTjtISjFJc5jYDAM/X8gbbKNkhC8lC/iQL0TLP2SHpUuX0rNnT2JiYhgwYACvvPIKfn7evR6TZK2cG4RME5JTITkFHE4wnWDY0n/tCwqEAOvv7Dh48CAvv/wyM2fOpGHDhnz66afcdtttWXItu93Ozz//TFRUFLNnz+bMmTOUK1fOFX6qVaumH/Ii4hYSExMZOHAgw4cPp3LlynzxxRfUq1fP6rIkh8o5Qcg04XwcnLkAcQnpfxzOqx/v5wuhwel/wvNAaFC2lZqYmMjQoUP58MMPyZcvHx9++CGPPfZYpgcRu93OypUrXeHn1KlTlClTxhV+atSoofAjIm5r48aNdOvWjY0bN/LCCy8wePBgQkK8a+82yXqeH4Tsdjh+GmJi07s/N7toSGgQFC0IEXkhi9bnMU2TqKgoXnnlFU6cOEG/fv148803M/V/bIfDwapVq4iKimLWrFmcPHmS0qVL06FDBzp06ECtWrUUfkTEY9jtdkaOHEn//v2JiIhg7NixWbafongnzw1CTiccOApHTpCpa+n7+ECpIlC0QPqMyUyyadMmXnzxRVatWsUDDzzAsGHDMm3xQYfDwZo1a1zh58SJE5QsWdIVfurUqaPwIyIebd++ffTq1YslS5bw6KOPMmLECAoUKGB1WZIDeGYQuhAPO/and4CySlgwVCwNuQJv6TSnTp2if//+jB8/nsjISEaMGMFdd911y+U5nU7Wrl1LVFQUM2fO5Pjx4xQvXtwVfurWravwIyI5immafPvtt7z00kuYpsmwYcN48skn9bNObolnBSGnE/bHpHeBsoNhQJli1+0OHT9+nKNHj1K7dm3XY3a7nbFjxzJgwABM02TgwIE8++yzt3T3g9PpZN26da7wc/ToUYoVK0b79u3p0KED9erVs2TbDRGR7HTy5EleeuklpkyZQqtWrRg3bpy295Gb5jlByOGE7XvSJ0NntyIRUK7EFcNQWloatWrVYu/evezevZtixYqxbNkyXnzxRbZv3063bt0YMmQIERERN3Vpp9PJr7/+SlRUFDNmzCAmJoYiRYq4wk+DBg0UfkTEKy1evJiePXsSGxvLwIED6du3L76+vlaXJR7GM4KQwwlbd8H5eOtqKBgOkaUuC0ODBw9mwIAB2Gw27rnnHvz9/Zk9ezaNGzfm008/vaRLlFGmabJhwwZX+Dl8+DCFChVyhZ9GjRop/IiIAPHx8QwYMICRI0dSo0YNJkyYQJ06dawuSzyI+wchpxO27YUz562uJH2IrGxxVxjatm0bNWvWxG63uw7Jnz8/I0eO5NFHH72hcWvTNPn9999d4efgwYMULFiQRx55hA4dOtC4cWN8fHwy/VMSEckJfvvtN7p3787WrVt56aWXGDRoEMHBwVaXJR7A/YPQ/hg4dMzqKv5RsTQUDMfhcFCvXj02bdqE05m+XpFhGFSqVIktW7ZkKLSYpsmff/5JVFQUUVFRHDhwgIiICFf4adKkicKPiEgGpaWlMXz4cAYOHEjhwoX5/PPPM+XmFMnZ3Ht8JS7BvUIQwO5DkJrGq6++yp9//ukKQZAebLZv387o0aOv+nLTNNm4cSNvvPEG5cqV47bbbmPSpEm0bt2aZcuWcfToUcaMGUPz5s0VgkREboCfnx+vvfYaW7dupVSpUrRu3ZonnniCU6dOWV2auDH37Qg5nfDHdkhMtrqSy4XnJvC2aqSk/HP7fnBwMMHBwZw9e5bw8HCOHj3qGhozTZMtW7a4Oj979uwhPDychx9+mA4dOtC8eXNN8BMRyUSmaTJ58mT69euHj48PI0aMyJIV/MXzuW8QOng0fcFENxWbP5Rzfgb58+cnd+7cnDhxgiZNmrB//37XkJevr68r/OzatYu8efO6wk+LFi20kaCISBY7ceIEL774ItOnT6d169Z8/vnnlCpVyuqyxI24ZxByOGHdZnA4rK7k6oJzQZ3KYBgcO3aM22+/nYMHD+JwODAMg3z58nH69Gny5MnDQw89RIcOHWjVqpXCj4iIBRYsWMAzzzzDmTNneO+993jhhRfUiRfAXYPQ8VMQfcDqKq6vVkWOJcRRv359jhw5wr+/lMHBwUyfPp0777wTf3/rd7oXEfF2cXFxvPXWW4wePZo6derwxRdfUKNGDavLEou552TpmFirK7g+Azh6ksqVK3P48GH+mycTEhIoUqSIQpCIiJsIDQ3l008/5ZdffiE5OZk6derw+uuvk5SUZHVpYiH3C0Lxiel/3J0JxJ7hqS5P0LRpU2rWrEloaOglhyxZssSa2kRE5KoaNGjAH3/8wcCBA/nkk0+oVq0ay5cvt7ossYj7DY0dPg77jlhdRcZVLQfheVwfnjlzhj179rB//36aNWtGoUKFrKtNRESuKTo6mh49erBq1Sq6du3Kxx9/TL58+awuS7KR+wWhHfsg9ozVVWRcqSJQsojVVYiIyE1yOp1MnDiRV155hYCAAEaOHEnHjh11q72XyLIgVKpUKQICAggMDCQhIYEqVarw2muv0ahRI9cxzZo1IyYmht27d7u+4Q7M+ZGaXTpw7oefLy20eV3Ozl9OntBQ0ux23v/2S6YtW4yPjw1/Xz9KFizMwP/rTs3ykazY+Ad9Rg9j08SprtcfOHaUmt06c+6Hn6n59GMApNrtRB8+SLXS6bsWR5YoydAez1O280OuxwAC/P35dexkDhw76nru4hft/b4v07Z3j2t+Ld5//32++uordu/ezezZs3nwwQdv8qsqIiJZ5dixYzz//PPMmjWLtm3bMmbMGEqUKGF1WZLFsvTewenTp1OzZk0AZs+eTZs2bVi8eDH169dn9+7d7N69m7x587Jy5UqaN28OdgekpF33vF0/GER8UhLrxkwib2gYAEt//5XowwepWT7yuq+/GJAuhqP/BqbQXEGXPPZv/35uwS+refSVvpzt+fQ1V4G+44476NSpE0899dR1axMREWsULlyYmTNnMnfuXHr37k2VKlUYMmQIzz77rFb6z8GybbL0ww8/TK9evfj4448BmDRpEo8//jjdunVj4sSJ6QelpF73PLuPHGLOmhVMeq2/KwQB3HFbfTq2zN49ZVrVqUtcYgJnrrN8e7169ShTpkw2VSUiIrfigQceYPv27XTp0oUXX3yRxo0bs3XrVqvLkiySratJ1a9fn3nz5uFwOPjqq69Yvnw5+fPnZ9CgQZw/f57cxvXL2bg7mnJFi5MvLPc1j4s+fMg1BAbpw2AZFZeUeMlrq5Quw5S3B1923MyVy2lZ+zYiwvNn+NwiIuL+cufOzZgxY+jcuTPdu3endu3avPbaa7z99tsEBgZaXZ5komwNQhenIy1cuJBSpUpRsWJFIH3oaOrUqTzzWJdrTk670nN7Y47QbsCrJKWm0KhKdb58/R0AIouXuOIcoYy41tDYxZB0Ju4Cp86fY/nwseBm881FRCRzNG7cmI0bN/LBBx8wZMgQZsyYwfjx42nWrJnVpUkmydZ1hH777TeqVq3KxIkT2bVrF6VKlaJUqVKsXr06fXjMZhCRJy/xSUmk/auDc+rcOQL9AwgLDqZW+Uj2xBzmbNwFAMoWLcamiVN547H/42xcXJZ/DhdD0sHp83n9sSfp9O5bJKemXP+FIiLikQICAnjnnXfYtGkT+fPnp3nz5vTo0YNz585ZXZpkgmwLQnPnzmXs2LH069ePZcuWsWfPHg4cOMCBAwc4duwYR48eZfO2bQQFBtKwSjUmLJjjeu3YeTO5o05dDMOgfLESPNC4GU9/+B7n/hV8EpKzd2VQwzDo/0Q38ufOw9hx47L12iIikv0qV67M6tWr+eyzz/juu++oVKkSM2fOvGxnAfEsWRqEOnbsSI0aNShXrhwTJ05k4cKFrFy5krvuuos8efL8U4TNRqdOnZg45RswDL55cxCLNqyj5tOPUePpx9iydw/jX37Ldfzk19+hWply1O/9f1T5vw7c/lw3lv7xG6899kSm1H1x+Ovff+ISEy47zjAMhvV9haEffkhi4tVXwx48eDDFihVj3bp1dOvWjWLFinHy5MlMqVVERLKPzWajd+/ebN++nfr169O+fXseeughYmJirC5NbpL7Laj4x3bP2GID0vcbKxAOFUtbXYmIiGQz0zSZPXs2zz33HAkJCQwdOpSePXtis7nf7lVyde73rxUWkh4wPIEJhAZbXYWIiFjAMAzatWvH9u3b6dSpE71796Zp06Zs377d6tLkBrhfEAoNAvfqUV1baBAA7777LjVr1rzsz969ey0uUEREslLevHkZP348K1as4OTJk9SsWZOBAweSkqIbaTyB+w2NpabB+s2eEYb8/aBBddB+NCIiAiQnJzN48GCGDh1K+fLlmTBhAo0bN7a6LLkG9+sI+ftBhIfs/FskQiFIRERcAgMDGTx4MH/++SdhYWHcfvvt9O7dm/Pnz1tdmlyF+wUhgCIFrK7g+gygcITVVYiIiBuqVq0aa9euZeTIkXz99ddUrlyZuXPnWl2WXIF7BqGwYAhy4yXMDSB/3vTulYiIyBX4+PjwwgsvsH37dmrVqsWDDz7II488wrFjx6wuTf7FPYOQYUDZ4lZXcQ0GlCpqdREiIuIBSpQowfz58/nuu+9YvXo1lSpVYvz48TidTqtLE9w1CAHkyw2F3HQz09JF3btjJSIibsUwDDp27MiOHTt4+OGH6dmzJy1atCA6Otrq0rye+wYhgLLFwD9b94W9NoP02+WLFbS6EhER8UD58uVj0qRJLF26lJiYGGrUqMHgwYNJTU21ujSv5d5ByNcXIt1o1WbDgIpldKeYiIjcklatWrF161b69OnDwIEDqVOnDuvXr7e6LK/kfusIXcmxk7DroLU1GEC1CpA3zNo6REQkR9m0aRPdunXjzz//5LnnnmPIkCGEhoZm2/XTHCYxpx0cPGnn6BkHKWkmdgf4+oCfj0FEbhslI3wpEeFLLv+c1wjwjCAEEHMC9hy25toGULkc5M9jzfVFRCRHs9vtfPrpp/Tv35/w8HDGjBnDvffem2XXO5fgZNX2FP7cm8qxsw6cfycBHxv8OxUYBjjNfx7LH2ajWgk/mlUNoGg+N5q6cgs8JwgBHD8F0Qey95o2G1Qtp06QiIhkuf379/PMM8+wePFiOnbsyMiRIylYMHPmpZqmSfRROz9vTWbj/rS/H7vx89j+DkflCvnSsnoAtUr74+vjuZ0izwpCAHEJsGM/JCVn/bVyh6TPUcoVkPXXEhERIT2wTJkyhT59+uB0Ovn444/p2rUrxi3MTz0d5+DL5QlEx9hdQeZWGUZ6kIoIs/H0HSGULeSZHSL3nix9JaHBcFtlKFEoa85vkB53y5WAGpEKQSIikq0Mw+Dxxx9n586d3HvvvTz99NPccccd7Nmzx3VMYmIib731FkePHr3muUzTZOW2ZPpPO8/uo3Ygc0JQ+rnT/3s6zsnQ2ReY+UsiaXbP6q2AJ3aE/i0uAQ4fh5NnM+d8NgMK5ofihRSARETELfz000/07NmT48eP884779CvXz/eeustPvroIx588EHmzJlzxdclpjgZuyienTH2bKnTACJy23i+TSiF8vpkyzUzg2cHoYtSUtPnD8XEQtrf/+AXe3ZXY/DPDve5AqBoQSiYL/2WfRERETeSkJDAO++8wyeffELZsmXZu3eva2XqRYsW0bp160uOv5DoZPi8uEsmQmcHmwGB/gb97g+lRIRnvJ/mjCB0kWlCYnJ6pyguEeLiISkFnM700GMY6f9KwbnSh9hCgyEkKD0IaW0gERFxc+vXr6dFixYkJ6fPk7XZbJQsWZIdO3YQEJA+khGfnD5UFXvema0h6CKbkb4W8msPh1Es3P3DkPtXeCOMv0NOcC7IoilEIiIiVlm7dq0rBAE4nU7279/P8OHDeeONN0izm4yYH2dZCIL0OUipdhg2N46324cRHurew2Q5qyMkIiKSQzmdTkJCQkhOTsbHxwe7/Z+5PzabjS1btrDjQml+2piMO7yx2wwoV9iXfg+EYnPjUZec1RESERHJoQzD4JtvvmHTpk3ExMRw5MgRDh48yKFDh0hOTmbyjOWcL9DZLUIQpHeGdh21s2pbCs2ruu9G5eoIiYiIeLjEpFSGzE7k1AXrhsSuxs8X3u2Um/xh7jlE5nnrCImIiMglft7m4KSF84KuxeGA79YkWl3GVSkIiYiIeDC7w2TZVveYF3QlThM2H0jj1AWH1aVckYKQiIiIB9tyII24JHeNQelsBqzanmJ1GVekICQiIuLBlm9Nxua+N2UB6V2hldtSSHO4X2BTEBIREfFQZ+IcRB+1u+XcoP9KTDHZejDN6jIuoyAkIiLiofbHuue8myvxscH+E9mz79mNUBASERHxUAdj7W4/LHaRw+meQUgLKoqIiHio/bFXHhbbv2k+mxZ9gul04LAnE5S7EG2enwOGwZalo9i1bgoXdx+PbNSFaq2exfh79ecvng2nSGRT2rzwz67237xangdfW0ZoeImr1rJ0wv8Ru38DiedP0OWjfQQE5b7smIMn7Zim6bqWO1AQEhER8VAHT14+NJZ4/jhrpvblwdeWExpeHIBThzaDYfD7vMEc37OO+/ouJDAknOT40ywZ14XUpAvcdt+brnNcOHWAI9uXU6xyywzXUqnJ/9G400dMeb3iVY9JToNTcU4i3GhxRQ2NiYiIeCCnaZKUenk7KOnCSQybDwHBeVyP5S9RA3tKAn8tH8vtjw0nMCQcgMCQcG5/bDhbl31GWkqC6/g6bV/nt7nvciObTxSt2JxcoRHXPS7ezW71VxASERHxQParzJPOV7QKhcrWZ3r/miwZ34UtS0aRcO4oZ49HY/P1J2/hSzs2eQtXxObjx7nj0a7HSlS7G9+AYPb+NjPT63a3W+gVhERERDzQ1Zo1hs3GHd2/4r5+P1KsUitO7PuVWYMbpz93A3Nz6j34Dn8s+B8Oe2pmlOu2FIREREQ8kN91ZvnmKVSBSk3+jzt7fkuBUrdxZPtyHGkpnD2285Ljzh7bidORRp5CkZc8XrBMPfIVrcyO1V9mbt0+7jNRGhSEREREPJLNMAj0vzxUJJw7yvG9v7o+Tkk8R9zpg4QXrUKVFj1YM60fyfGnAUiOP8Oaaf2o2vIZ/AKCLzvXbfe/zebFI3DaM297jJBA9wpCumtMRETEQ5XM70P00UvX5nE6HGz88SPiTh/C1y8XptNO+fqdKFmjDSWq3c2WpZ8yb9g9GIYNTCcVGj1O9Tuev+L58xauSPGqd/59u/21LR7TidMxfwEwa3AjwgqU5d4+8y45JsAP8oe5Vw/GMG9kSriIiIi4jVnrElmyORmH0+pKMiayiC8vPxhmdRmXcK9YJiIiIhlWMsLXY0KQjw1KFXS/gSj3q0hEREQypHTB7F2Y8KfPOxN/5sgljwUE5aFtn7nXfa3DCWUKuF/scL+KREREJEPCQ30oX9iXvcezZwf6u3pdf67Q1QT5G1Qr6ZeJ1WQODY2JiIh4sFbVA7MlBN0KmwFNqwTg5+ted4yBgpCIiIhHq1HKj1A3uyX9v5xmehByRwpCIiIiHszXx6Bl9UDcNQrZDKhe0s+tNlr9NwUhERERD3dXzUDCw2zcwA4a2cbHBp1uD7K6jKtSEBIREfFw/r4GT7cKvur+Y1Zq3yiIiNzu2Q0CBSEREZEcoVxhP+6q4T5DZDYDyhf2pVlV95wbdJGCkIiISA7xQP1cFMvvg83iNGQzIJe/QddWwdjccbzuXxSEREREcgh/X4OX7gslf5jNsjBkGODnC/0eCHXbCdL/piAkIiKSg4TmsvHqQ2EUzJP9YchmQKCfwcsPhFE8v2es2axNV0VERHKghGQnYxbFs+s/u9NnFcOA8BAbL7QNpXA+9+8EXaQgJCIikkM5TZOVf6Uw45dEHE6yZAVqm5F+3jtqBPBgvSAC/Nx7TtB/KQiJiIjkcCfPO/hyeQK7j9ldweVWGQaYJoSH2nj6jmDKF3a/fcQyQkFIRETECzhNk51H7CzfmsyWA2nwd5C5UReDVOmCPrSqFkjtsv74+XhWF+jfFIRERES8zJk4Byu3p/Dn3lROnHdimmAANlt6ODLBtR6RyT+BKW+IjWol/GheNcBjJkNfj4KQiIiIF0tJMzly2sGhk3aOnnGQYjexO9K3xvD3NYjIbaNkhC8l8vsQHJjzbjZXEBIRERGvlfOinYiIiEgGKQiJiIiI11IQEhEREa+lICQiIiJeS0FIREREvJaCkIiIiHgtBSERERHxWgpCIiIi4rUUhERERMRrKQiJiIiI11IQEhEREa+lICQiIiJeS0FIREREvJaCkIiIiHgtBSERERHxWgpCIiIi4rUUhERERMRrKQiJiIiI11IQEhEREa+lICQiIiJeS0FIREREvJaCkIiIiHgtBSERERHxWv8P5Bk16obUdasAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "recombulatorx.families.plot_family_graph(family_graphs[0][1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "PMTL-bbf3U17",
   "metadata": {
    "id": "PMTL-bbf3U17"
   },
   "source": [
    "The fuction **preprocess_families** will then check the consistency of each family graph and raise errors whenever necessary. For instance, an error is raised when more than two parents or same-sex parents are present in the same family. Unconnected individuals are also flagged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c1ac7846-b84e-4f97-b5ab-fadd1491b0eb",
   "metadata": {
    "id": "c1ac7846-b84e-4f97-b5ab-fadd1491b0eb",
    "outputId": "23abcbfe-9992-499d-a644-21861bd01345"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "200"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "processed_families = recombulatorx.preprocess_families(family_graphs)\n",
    "len(processed_families)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0986dd9d-8c60-4ae7-8bf1-5c7cc45b0c21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ProcessedFamily(fid=('FAM_0_I', 'MOTHER'), is_mother_phased=True, mother=array([[12., 16.],\n",
       "       [29., 27.],\n",
       "       [39., 34.],\n",
       "       [43., 49.],\n",
       "       [59., 54.],\n",
       "       [68., 66.]]), maternal_haplotypes=array([[12., 29., 34., 49., 54., 66.],\n",
       "       [16., 27., 34., 49., 54., 66.],\n",
       "       [12., 29., 39., 43., 59., 68.]]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "processed_families[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tbunWmze3tMa",
   "metadata": {
    "id": "tbunWmze3tMa"
   },
   "source": [
    "## Estimation\n",
    "\n",
    "The estimation of recombination and mutation rates can be launched with the function **estimate_rates**, which estimates recombination and optionally mutation rates from a set of families. The function takes the following parameters: \n",
    "\n",
    "*   the families, \n",
    "*   the initial recombination rate,\n",
    "*   which mutation rate needs to be estimated (no: no mutation rate estimation, one: just one mutation rate for all markers, all: a mutation rate for each marker),\n",
    "*   the type of implementation (the default implementation is the one using dynamic programming)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "58ed4780",
   "metadata": {},
   "source": [
    "While the numba package is not required for recombulatorx, when it is available the estimation will be much faster!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d88ff711-f4d9-4ba4-a201-6a4c86f7967f",
   "metadata": {
    "id": "d88ff711-f4d9-4ba4-a201-6a4c86f7967f"
   },
   "outputs": [],
   "source": [
    "est_recomb_rates, est_mut_rates = recombulatorx.estimate_rates(processed_families, 0.1, 0.1, estimate_mutation_rates='all')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "P8nOUjcR6Tba",
   "metadata": {
    "id": "P8nOUjcR6Tba"
   },
   "source": [
    "The output will be returned as a tuple of arrays, where the first stores the recombination rates and the second the mutation rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a0da2bb3-e67f-4f7a-b821-b7b944168b71",
   "metadata": {
    "id": "a0da2bb3-e67f-4f7a-b821-b7b944168b71",
    "outputId": "9085f5ef-c02d-4684-ab45-190262f999ab"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.03694946, 0.32475622, 0.012345  , 0.1856328 , 0.00980138]),\n",
       " array([1.00000000e-08, 1.00000000e-08, 8.32904925e-03, 1.00000000e-08,\n",
       "        2.15161190e-03, 1.00000000e-08]))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est_recomb_rates, est_mut_rates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "B8VfzmiU8FuI",
   "metadata": {
    "id": "B8VfzmiU8FuI"
   },
   "source": [
    "In case the dataset has been simulated, a direct comparison of initial and estimated recombination rates is possible as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "677dd957-ec90-4875-86e9-9eca01c13010",
   "metadata": {
    "id": "677dd957-ec90-4875-86e9-9eca01c13010",
    "outputId": "04788c6c-355e-437c-c931-3ea2984e4dff"
   },
   "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>simulated</th>\n",
       "      <th>estimated</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.05</td>\n",
       "      <td>0.036949</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.30</td>\n",
       "      <td>0.324756</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.01</td>\n",
       "      <td>0.012345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.185633</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.01</td>\n",
       "      <td>0.009801</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   simulated  estimated\n",
       "0       0.05   0.036949\n",
       "1       0.30   0.324756\n",
       "2       0.01   0.012345\n",
       "3       0.20   0.185633\n",
       "4       0.01   0.009801"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGYCAYAAABhxLkXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtPElEQVR4nO3dfVRVdb7H8c8B5Ul5UpQHhzyYlnpHQCUJe7BVZwSnmetDD+pUENerN701GqkTpoDXCjJz0DKZ7JJUt9FqpbWmIos7OFaEBam31FKvDpoC6lw5CgkG5/7R8jhnRPMgDz/h/Vprrzh7//Zvfzdb1vm092/vbXE4HA4BAAAYzKOjCwAAAPgpBBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPG6dXQBraGpqUmHDx+Wv7+/LBZLR5cDAAAugcPh0MmTJxURESEPj4ufQ+kUgeXw4cOKjIzs6DIAAEALHDx4UD/72c8u2qZTBBZ/f39JP+5wQEBAB1cDAAAuhd1uV2RkpPN7/GI6RWA5exkoICCAwAIAwBXmUoZzMOgWAAAYj8ACAACMR2ABAADG6xRjWAAAnUNTU5MaGho6ugy0ou7du8vT0/Oy+yGwAACM0NDQoP3796upqamjS0ErCwoKUlhY2GU9K43AAgDocA6HQ0eOHJGnp6ciIyN/8iFiuDI4HA7V1dWpurpakhQeHt7ivggsAIAO98MPP6iurk4RERHy8/Pr6HLQinx9fSVJ1dXV6tu3b4svDxFhAQAdrrGxUZLk5eXVwZWgLZwNoWfOnGlxHwQWAIAxeB9c59Qax5XAAgAAjEdgAQCgDdx///2aMGFCm28nKytLsbGxbb6dC7FYLNq4cWObb4dBtwAAY1kffbddt3cg5/ZW62vFihVyOByt1l9rslgs2rBhQ7sEqtZCYAEAoA0EBgZ2dAmdCpeEAAC4DG+++aaGDRsmX19f9e7dWzabTbW1teddErrlllv00EMPac6cOQoODlZoaKjWrFmj2tpapaamyt/fXwMHDtT777/vXGft2rUKCgpy2d7GjRsvOoj1888/1y9+8QuFhIQoMDBQY8aMUXl5uXO51WqVJE2cOFEWi8X5WZLefvttjRgxQj4+PhowYIAWL16sH374wbl8z549uvnmm+Xj46OhQ4fqww8/bNkvrQU4wwJcSFY7/N9RVk3bbwNAmzly5IimTp2qpUuXauLEiTp58qS2bNlywUtBBQUFmj9/vrZu3ar169dr5syZ2rBhgyZOnKgFCxbo97//ve677z5VVFS0+Hk0J0+eVEpKip599lk5HA4988wz+uUvf6k9e/bI399fn3/+ufr27auXXnpJSUlJzueibNmyRcnJyVq5cqVuuukm7du3TzNmzJAkZWZmqqmpSZMmTVJoaKhKS0tVU1OjOXPmtKjGliCwAADQQkeOHNEPP/ygSZMmqX///pKkYcOGXbB9TEyMFi5cKElKT09XTk6OQkJCNH36dElSRkaGVq9erR07duj6669vUU233nqry+cXXnhBQUFB2rx5s371q1+pT58+ks49Lv+sxYsX69FHH1VKSookacCAAVqyZInmz5+vzMxMffTRR9q9e7c++OADRURESJKefPJJjRs3rkV1uovAAgBAC8XExOi2227TsGHDlJiYqLFjx+rOO+9UcHBws+2jo6OdP3t6eqp3794uASc0NFSSnI+yb4mqqiotXLhQxcXFqq6uVmNjo+rq6lRRUXHR9bZv365PPvlETzzxhHNeY2OjTp8+rbq6Ou3atUuRkZHOsCJJCQkJLa7TXQQWAABayNPTUx9++KE+/fRTbdq0Sc8++6wee+wxlZaWNtu+e/fuLp8tFovLvLNjU86+ANLDw+O8y0s/9bTYlJQUHT9+XCtWrFD//v3l7e2thISEn3wL9qlTp7R48WJNmjTpvGU+Pj4XXbc9EFgAALgMFotFN9xwg2644QZlZGSof//+2rBhQ6v03adPH508eVK1tbXq0aOHJGnbtm0XXeeTTz7R888/r1/+8peSpIMHD+rYsWMubbp37+58HcJZI0aM0DfffKOBAwc22++QIUN08OBBHTlyxPkSw88++6wlu9UiBBYAAFqotLRURUVFGjt2rPr27avS0lIdPXpUQ4YM0Y4dOy67//j4ePn5+WnBggX67W9/q9LSUq1du/ai6wwaNEivvPKK4uLiZLfbNW/ePOcLCM+yWq0qKirSDTfcIG9vbwUHBysjI0O/+tWvdNVVV+nOO++Uh4eHtm/frq+++kqPP/64bDabrrnmGqWkpOjpp5+W3W7XY489dtn7eKkILAAAY7Xmg9zaQkBAgP7yl78oNzdXdrtd/fv31zPPPKNx48Zp/fr1l91/r1699Oqrr2revHlas2aNbrvtNmVlZTnv3mnOf/7nf2rGjBkaMWKEIiMj9eSTT2ru3LkubZ555hmlpaVpzZo16tevnw4cOKDExET96U9/0n/8x3/oqaeeUvfu3TV48GD967/+q6QfL09t2LBB06ZN06hRo2S1WrVy5UolJSVd9n5eCovD1MfwucFutyswMFA1NTUKCAjo6HLQWXBbM9BuTp8+rf379ysqKsqI8RJoXRc6vu58f/PgOAAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAAAMcv/992vChAkdXYZb2qNmAgsAAB3gwIEDslgs573McMWKFT/5vqDWcKUFI94lBAAwV3u8IsNlex3/uozAwHbe5ysEZ1gAALgMTU1Nys7OVlRUlHx9fRUTE6M333xTkvR///d/uueee9SnTx/5+vpq0KBBeumllyRJUVFRkqThw4fLYrHolltukXT+mY9bbrlFDz30kObMmaPg4GCFhoZqzZo1qq2tVWpqqvz9/TVw4EC9//77znUaGxs1bdo0Z03XXnutVqxY4VyelZWlgoICvf3227JYLLJYLCouLpYkHTx4UHfffbeCgoLUq1cvjR8/XgcOHHDpOy0tTUFBQerdu7fmz5+v9ngtIYEFAIDLkJ2drZdffll5eXn6+uuv9fDDD+vee+/V5s2btWjRIu3cuVPvv/++du3apdWrVyskJESStHXrVknSRx99pCNHjuitt9664DYKCgoUEhKirVu36qGHHtLMmTN11113afTo0SovL9fYsWN13333qa6uTtKPIepnP/uZ3njjDe3cuVMZGRlasGCBXn/9dUnS3LlzdffddyspKUlHjhzRkSNHNHr0aJ05c0aJiYny9/fXli1b9Mknn6hnz55KSkpSQ0ODpB/f9Lx27Vrl5+fr448/1t/+9jdt2LChLX/FkrgkBABAi9XX1+vJJ5/URx99pISEBEnSgAED9PHHH+sPf/iDTp06peHDhysuLk6SZLVanev26dNHktS7d2+FhYVddDsxMTFauHChJCk9PV05OTkKCQnR9OnTJUkZGRlavXq1duzYoeuvv17du3fX4sWLnetHRUWppKREr7/+uu6++2717NlTvr6+qq+vd9n2q6++qqamJr344ouyWCySpJdeeklBQUEqLi7W2LFjlZubq/T0dE2aNEmSlJeXpw8++OByfo2XhMACAEAL7d27V3V1dfrFL37hMr+hoUHDhw9XVlaW7rjjDudZkAkTJmj06NFubyc6Otr5s6enp3r37q1hw4Y554WGhkqSqqurnfNWrVql/Px8VVRU6Pvvv1dDQ4NiY2Mvup3t27dr79698vf3d5l/+vRp7du3TzU1NTpy5Iji4+Ody7p166a4uLg2vyzUoktCq1atktVqlY+Pj+Lj452ntZrz1ltvKS4uTkFBQerRo4diY2P1yiuvuLRxOBzKyMhQeHi4fH19ZbPZtGfPnpaUBgBAuzl16pQk6d1339W2bduc086dO/Xmm29q3Lhx+utf/6qHH35Yhw8f1m233aa5c+e6vZ3u3bu7fLZYLC7zzp4NaWpqkiStW7dOc+fO1bRp07Rp0yZt27ZNqampzss6F9ufkSNHuuzLtm3b9O233+o3v/mN23W3JrcDy/r165WWlqbMzEyVl5crJiZGiYmJLqnu7/Xq1UuPPfaYSkpKtGPHDqWmpio1NdXl9NHSpUu1cuVK5eXlqbS0VD169FBiYqJOnz7d8j0DAKCNDR06VN7e3qqoqNDAgQNdpsjISEk/XvpJSUnRq6++qtzcXL3wwguSJC8vL0k/DmJtbZ988olGjx6tWbNmafjw4Ro4cKD27dvn0sbLy+u8bY8YMUJ79uxR3759z9ufwMBABQYGKjw8XKWlpc51fvjhB5WVlbX6PvwjtwPL8uXLNX36dKWmpmro0KHKy8uTn5+f8vPzm21/yy23aOLEiRoyZIiuvvpqzZ49W9HR0fr4448l/Xh2JTc3VwsXLtT48eMVHR2tl19+WYcPH9bGjRsva+cAAGhL/v7+mjt3rh5++GEVFBRo3759Ki8v17PPPquCggJlZGTo7bff1t69e/X111/rT3/6k4YMGSJJ6tu3r3x9fVVYWKiqqirV1LTeLdWDBg3SF198oQ8++EDffvutFi1apM8//9yljdVq1Y4dO/TNN9/o2LFjOnPmjO655x6FhIRo/Pjx2rJli/bv36/i4mL99re/1aFDhyRJs2fPVk5OjjZu3Kjdu3dr1qxZOnHiRKvVfiFuBZaGhgaVlZXJZrOd68DDQzabTSUlJT+5vsPhUFFRkb755hvdfPPNkqT9+/ersrLSpc/AwEDFx8dfsM/6+nrZ7XaXCQCAjrBkyRItWrRI2dnZGjJkiJKSkvTuu+8qKipKXl5eSk9PV3R0tG6++WZ5enpq3bp1kn4c+7Fy5Ur94Q9/UEREhMaPH99qNf3bv/2bJk2apMmTJys+Pl7Hjx/XrFmzXNpMnz5d1157reLi4tSnTx998skn8vPz01/+8hddddVVmjRpkoYMGaJp06bp9OnTCggIkCQ98sgjuu+++5SSkqKEhAT5+/tr4sSJrVb7hVgcboySOXz4sPr166dPP/3UORpakubPn6/Nmze7nCL6ezU1NerXr5/q6+vl6emp559/Xv/yL/8iSfr00091ww036PDhwwoPD3euc/fdd8tisWj9+vXn9ZeVleUy+vnvt3P2F4rOzfrou22+jQM+7XC91oCHVAEmOH36tPbv36+oqCj5+Ph0dDloZRc6vna7XYGBgZf0/d0uz2Hx9/fXtm3b9Pnnn+uJJ55QWlqa8wE1LZGenq6amhrndPDgwdYrFgAAGMet25pDQkLk6empqqoql/lVVVUXvYfcw8NDAwcOlCTFxsZq165dys7O1i233OJcr6qqyuUMS1VV1QVvv/L29pa3t7c7pQMAgCuYW2dYvLy8NHLkSBUVFTnnNTU1qaioyOUS0U9pampSfX29pB8fZhMWFubSp91uV2lpqVt9AgCAzsvtB8elpaUpJSVFcXFxGjVqlHJzc53vM5Ck5ORk9evXT9nZ2ZJ+fGRxXFycrr76atXX1+u9997TK6+8otWrV0v68d7xOXPm6PHHH9egQYMUFRWlRYsWKSIi4op6iyQAAGg7bgeWyZMn6+jRo8rIyFBlZaViY2NVWFjofMpeRUWFPDzOnbipra3VrFmzdOjQIfn6+mrw4MF69dVXNXnyZGeb+fPnq7a2VjNmzNCJEyd04403qrCwkIFXAABAkpt3CZnKnVHG6By4SwjoXM7eRWK1WuXr69vR5aCVff/99zpw4ID5dwkBAHAxnp6ekvSTj47HlensW6T/8RUD7uDlhwCADtetWzf5+fnp6NGj6t69u8vQAly5HA6H6urqVF1draCgIGcwbQkCCwCgw1ksFoWHh2v//v3661//2tHloJUFBQVd9PEnl4LAAgAwgpeXlwYNGsRloU6me/ful3Vm5SwCCwDAGB4eHtwhimZxkRAAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeC0KLKtWrZLVapWPj4/i4+O1devWC7Zds2aNbrrpJgUHBys4OFg2m+289vfff78sFovLlJSU1JLSAABAJ+R2YFm/fr3S0tKUmZmp8vJyxcTEKDExUdXV1c22Ly4u1tSpU/XnP/9ZJSUlioyM1NixY/Xdd9+5tEtKStKRI0ec0x//+MeW7REAAOh03A4sy5cv1/Tp05WamqqhQ4cqLy9Pfn5+ys/Pb7b9f/3Xf2nWrFmKjY3V4MGD9eKLL6qpqUlFRUUu7by9vRUWFuacgoODW7ZHAACg03ErsDQ0NKisrEw2m+1cBx4estlsKikpuaQ+6urqdObMGfXq1ctlfnFxsfr27atrr71WM2fO1PHjxy/YR319vex2u8sEAAA6L7cCy7Fjx9TY2KjQ0FCX+aGhoaqsrLykPn73u98pIiLCJfQkJSXp5ZdfVlFRkZ566ilt3rxZ48aNU2NjY7N9ZGdnKzAw0DlFRka6sxsAAOAK0609N5aTk6N169apuLhYPj4+zvlTpkxx/jxs2DBFR0fr6quvVnFxsW677bbz+klPT1daWprzs91uJ7QAANCJuXWGJSQkRJ6enqqqqnKZX1VVpbCwsIuuu2zZMuXk5GjTpk2Kjo6+aNsBAwYoJCREe/fubXa5t7e3AgICXCYAANB5uRVYvLy8NHLkSJcBs2cH0CYkJFxwvaVLl2rJkiUqLCxUXFzcT27n0KFDOn78uMLDw90pDwAAdFJu3yWUlpamNWvWqKCgQLt27dLMmTNVW1ur1NRUSVJycrLS09Od7Z966iktWrRI+fn5slqtqqysVGVlpU6dOiVJOnXqlObNm6fPPvtMBw4cUFFRkcaPH6+BAwcqMTGxlXYTAABcydwewzJ58mQdPXpUGRkZqqysVGxsrAoLC50DcSsqKuThcS4HrV69Wg0NDbrzzjtd+snMzFRWVpY8PT21Y8cOFRQU6MSJE4qIiNDYsWO1ZMkSeXt7X+buAQCAzsDicDgcHV3E5bLb7QoMDFRNTQ3jWboI66Pvtvk2Dvj8ps23oayatt8GABjKne9v3iUEAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA47UosKxatUpWq1U+Pj6Kj4/X1q1bL9h2zZo1uummmxQcHKzg4GDZbLbz2jscDmVkZCg8PFy+vr6y2Wzas2dPS0oDAACdkNuBZf369UpLS1NmZqbKy8sVExOjxMREVVdXN9u+uLhYU6dO1Z///GeVlJQoMjJSY8eO1Xfffedss3TpUq1cuVJ5eXkqLS1Vjx49lJiYqNOnT7d8zwAAQKdhcTgcDndWiI+P13XXXafnnntOktTU1KTIyEg99NBDevTRR39y/cbGRgUHB+u5555TcnKyHA6HIiIi9Mgjj2ju3LmSpJqaGoWGhmrt2rWaMmXKT/Zpt9sVGBiompoaBQQEuLM7uEJZH323zbdxwOc3bb4NZdW0/TYAwFDufH+7dYaloaFBZWVlstls5zrw8JDNZlNJSckl9VFXV6czZ86oV69ekqT9+/ersrLSpc/AwEDFx8dfcp8AAKBz6+ZO42PHjqmxsVGhoaEu80NDQ7V79+5L6uN3v/udIiIinAGlsrLS2cc/9nl22T+qr69XfX2987Pdbr/kfQAAAFcetwLL5crJydG6detUXFwsHx+fFveTnZ2txYsXt2JlANA5tMvl0pzb23wbwD9y65JQSEiIPD09VVVV5TK/qqpKYWFhF1132bJlysnJ0aZNmxQdHe2cf3Y9d/pMT09XTU2Nczp48KA7uwEAAK4wbgUWLy8vjRw5UkVFRc55TU1NKioqUkJCwgXXW7p0qZYsWaLCwkLFxcW5LIuKilJYWJhLn3a7XaWlpRfs09vbWwEBAS4TAADovNy+JJSWlqaUlBTFxcVp1KhRys3NVW1trVJTUyVJycnJ6tevn7KzsyVJTz31lDIyMvTaa6/JarU6x6X07NlTPXv2lMVi0Zw5c/T4449r0KBBioqK0qJFixQREaEJEya03p4CAIArltuBZfLkyTp69KgyMjJUWVmp2NhYFRYWOgfNVlRUyMPj3Imb1atXq6GhQXfeeadLP5mZmcrKypIkzZ8/X7W1tZoxY4ZOnDihG2+8UYWFhZc1zgUAAHQebj+HxUQ8h6Xr4TksQPMYdIsrSZs9hwUAAKAjEFgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOO5/Wh+AEAXlxXYxv3zBGicjzMsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGC8FgWWVatWyWq1ysfHR/Hx8dq6desF23799de64447ZLVaZbFYlJube16brKwsWSwWl2nw4MEtKQ0AAHRCbgeW9evXKy0tTZmZmSovL1dMTIwSExNVXV3dbPu6ujoNGDBAOTk5CgsLu2C///RP/6QjR444p48//tjd0gAAQCfldmBZvny5pk+frtTUVA0dOlR5eXny8/NTfn5+s+2vu+46Pf3005oyZYq8vb0v2G+3bt0UFhbmnEJCQtwtDQAAdFJuBZaGhgaVlZXJZrOd68DDQzabTSUlJZdVyJ49exQREaEBAwbonnvuUUVFxQXb1tfXy263u0wAAKDzciuwHDt2TI2NjQoNDXWZHxoaqsrKyhYXER8fr7Vr16qwsFCrV6/W/v37ddNNN+nkyZPNts/OzlZgYKBzioyMbPG2AQCA+Yy4S2jcuHG66667FB0drcTERL333ns6ceKEXn/99Wbbp6enq6amxjkdPHiwnSsGAADtqZs7jUNCQuTp6amqqiqX+VVVVRcdUOuuoKAgXXPNNdq7d2+zy729vS86HgYAAHQubp1h8fLy0siRI1VUVOSc19TUpKKiIiUkJLRaUadOndK+ffsUHh7ean0CAIArl1tnWCQpLS1NKSkpiouL06hRo5Sbm6va2lqlpqZKkpKTk9WvXz9lZ2dL+nGg7s6dO50/f/fdd9q2bZt69uypgQMHSpLmzp2rX//61+rfv78OHz6szMxMeXp6aurUqa21nwAA4ArmdmCZPHmyjh49qoyMDFVWVio2NlaFhYXOgbgVFRXy8Dh34ubw4cMaPny48/OyZcu0bNkyjRkzRsXFxZKkQ4cOaerUqTp+/Lj69OmjG2+8UZ999pn69OlzmbsHAAA6A7cDiyQ9+OCDevDBB5tddjaEnGW1WuVwOC7a37p161pSBgAA6CKMuEsIAADgYggsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYLwWBZZVq1bJarXKx8dH8fHx2rp16wXbfv3117rjjjtktVplsViUm5t72X0CAICuxe3Asn79eqWlpSkzM1Pl5eWKiYlRYmKiqqurm21fV1enAQMGKCcnR2FhYa3SJwAA6FrcDizLly/X9OnTlZqaqqFDhyovL09+fn7Kz89vtv11112np59+WlOmTJG3t3er9AkAALoWtwJLQ0ODysrKZLPZznXg4SGbzaaSkpIWFdCSPuvr62W3210mAADQebkVWI4dO6bGxkaFhoa6zA8NDVVlZWWLCmhJn9nZ2QoMDHROkZGRLdo2AAC4MlyRdwmlp6erpqbGOR08eLCjSwIAAG2omzuNQ0JC5OnpqaqqKpf5VVVVFxxQ2xZ9ent7X3A8DAAA6HzcOsPi5eWlkSNHqqioyDmvqalJRUVFSkhIaFEBbdEnAADoXNw6wyJJaWlpSklJUVxcnEaNGqXc3FzV1tYqNTVVkpScnKx+/fopOztb0o+Danfu3On8+bvvvtO2bdvUs2dPDRw48JL6BAAAXZvbgWXy5Mk6evSoMjIyVFlZqdjYWBUWFjoHzVZUVMjD49yJm8OHD2v48OHOz8uWLdOyZcs0ZswYFRcXX1KfAACga7M4HA5HRxdxuex2uwIDA1VTU6OAgICOLgftwProu22+jQM+v2nzbSirpu23gS6lU/xt8HfRZbjz/X1F3iUEAAC6FgILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGK9FgWXVqlWyWq3y8fFRfHy8tm7detH2b7zxhgYPHiwfHx8NGzZM7733nsvy+++/XxaLxWVKSkpqSWkAAKATcjuwrF+/XmlpacrMzFR5ebliYmKUmJio6urqZtt/+umnmjp1qqZNm6Yvv/xSEyZM0IQJE/TVV1+5tEtKStKRI0ec0x//+MeW7REAAOh03A4sy5cv1/Tp05WamqqhQ4cqLy9Pfn5+ys/Pb7b9ihUrlJSUpHnz5mnIkCFasmSJRowYoeeee86lnbe3t8LCwpxTcHBwy/YIAAB0Om4FloaGBpWVlclms53rwMNDNptNJSUlza5TUlLi0l6SEhMTz2tfXFysvn376tprr9XMmTN1/PjxC9ZRX18vu93uMgEAgM6rmzuNjx07psbGRoWGhrrMDw0N1e7du5tdp7Kystn2lZWVzs9JSUmaNGmSoqKitG/fPi1YsEDjxo1TSUmJPD09z+szOztbixcvdqf0VmF99N0238aBnNvbfBsAAFxp3AosbWXKlCnOn4cNG6bo6GhdffXVKi4u1m233XZe+/T0dKWlpTk/2+12RUZGtkutAACg/bl1SSgkJESenp6qqqpymV9VVaWwsLBm1wkLC3OrvSQNGDBAISEh2rt3b7PLvb29FRAQ4DIBAIDOy63A4uXlpZEjR6qoqMg5r6mpSUVFRUpISGh2nYSEBJf2kvThhx9esL0kHTp0SMePH1d4eLg75QEAgE7K7buE0tLStGbNGhUUFGjXrl2aOXOmamtrlZqaKklKTk5Wenq6s/3s2bNVWFioZ555Rrt371ZWVpa++OILPfjgg5KkU6dOad68efrss8904MABFRUVafz48Ro4cKASExNbaTcBAMCVzO0xLJMnT9bRo0eVkZGhyspKxcbGqrCw0DmwtqKiQh4e53LQ6NGj9dprr2nhwoVasGCBBg0apI0bN+rnP/+5JMnT01M7duxQQUGBTpw4oYiICI0dO1ZLliyRt7d3K+0mAAC4krVo0O2DDz7oPEPyj4qLi8+bd9ddd+muu+5qtr2vr68++OCDlpQBAAC6CN4lBAAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwnhEvP8TfyQpsh23UtP02AABoRZxhAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA43Xr6AIA4CdlBbbDNmrafhsAWowzLAAAwHgEFgAAYDwuCQG4bNZH323T/g/4tGn3QJto678LSTrg85s234Ypl0s5wwIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgtCiyrVq2S1WqVj4+P4uPjtXXr1ou2f+ONNzR48GD5+Pho2LBheu+991yWOxwOZWRkKDw8XL6+vrLZbNqzZ09LSgMAAJ2Q24Fl/fr1SktLU2ZmpsrLyxUTE6PExERVV1c32/7TTz/V1KlTNW3aNH355ZeaMGGCJkyYoK+++srZZunSpVq5cqXy8vJUWlqqHj16KDExUadPn275ngEAgE7D7cCyfPlyTZ8+XampqRo6dKjy8vLk5+en/Pz8ZtuvWLFCSUlJmjdvnoYMGaIlS5ZoxIgReu655yT9eHYlNzdXCxcu1Pjx4xUdHa2XX35Zhw8f1saNGy9r5wAAQOfg1oPjGhoaVFZWpvT0dOc8Dw8P2Ww2lZSUNLtOSUmJ0tLSXOYlJiY6w8j+/ftVWVkpm83mXB4YGKj4+HiVlJRoypQp5/VZX1+v+vp65+eamh8famO3293ZHbc11de1af+SZLc42nwbauPfU3vgWJilrY8Hx+LSdYq/DY7FJbvS/zbOfm87HD+9H24FlmPHjqmxsVGhoaEu80NDQ7V79+5m16msrGy2fWVlpXP52XkXavOPsrOztXjx4vPmR0ZGXtqOGKwdXvEm5bTLVq54HAtzcCzM0ua/KY7FJessfxsnT55UYODFt3NFPpo/PT3d5axNU1OT/va3v6l3796yWCwdWNnlsdvtioyM1MGDBxUQENDR5XRpHAtzcCzMwvEwR2c4Fg6HQydPnlRERMRPtnUrsISEhMjT01NVVVUu86uqqhQWFtbsOmFhYRdtf/a/VVVVCg8Pd2kTGxvbbJ/e3t7y9vZ2mRcUFOTOrhgtICDgiv3H19lwLMzBsTALx8McV/qx+KkzK2e5NejWy8tLI0eOVFFRkXNeU1OTioqKlJCQ0Ow6CQkJLu0l6cMPP3S2j4qKUlhYmEsbu92u0tLSC/YJAAC6FrcvCaWlpSklJUVxcXEaNWqUcnNzVVtbq9TUVElScnKy+vXrp+zsbEnS7NmzNWbMGD3zzDO6/fbbtW7dOn3xxRd64YUXJEkWi0Vz5szR448/rkGDBikqKkqLFi1SRESEJkyY0Hp7CgAArlhuB5bJkyfr6NGjysjIUGVlpWJjY1VYWOgcNFtRUSEPj3MnbkaPHq3XXntNCxcu1IIFCzRo0CBt3LhRP//5z51t5s+fr9raWs2YMUMnTpzQjTfeqMLCQvn4dK13ynt7eyszM/O8y11ofxwLc3AszMLxMEdXOxYWx6XcSwQAANCBeJcQAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAHCF6kr3zVyRj+bvLI4dO6b8/HyVlJQ435sUFham0aNH6/7771efPn06uEIAgMm8vb21fft2DRkypKNLaXPc1txBPv/8cyUmJsrPz082m835HJuqqioVFRWprq5OH3zwgeLi4jq4UkjSwYMHlZmZqfz8/I4updP7/vvvVVZWpl69emno0KEuy06fPq3XX39dycnJHVRd17Nr1y599tlnSkhI0ODBg7V7926tWLFC9fX1uvfee3Xrrbd2dIldwt+/P+/vrVixQvfee6969+4tSVq+fHl7ltWuCCwd5Prrr1dMTIzy8vLOe2Gjw+HQAw88oB07dqikpKSDKsTf2759u0aMGKHGxsaOLqVT+/bbbzV27FhVVFTIYrHoxhtv1Lp165zvGauqqlJERATHoZ0UFhZq/Pjx6tmzp+rq6rRhwwYlJycrJiZGTU1N2rx5szZt2kRoaQceHh6KiYk57715mzdvVlxcnHr06CGLxaL//u//7pgC2wGBpYP4+vrqyy+/1ODBg5tdvnv3bg0fPlzff/99O1fWNb3zzjsXXf6///u/euSRR/iibGMTJ07UmTNntHbtWp04cUJz5szRzp07VVxcrKuuuorA0s5Gjx6tW2+9VY8//rjWrVunWbNmaebMmXriiSckSenp6SorK9OmTZs6uNLOLycnRy+88IJefPFFl4DYvXt3bd++/byzkZ2SAx3CarU6CgoKLri8oKDA0b9///YrqIuzWCwODw8Ph8ViueDk4eHR0WV2en379nXs2LHD+bmpqcnxwAMPOK666irHvn37HJWVlRyHdhQQEODYs2ePw+FwOBobGx3dunVzlJeXO5f/z//8jyM0NLSjyutytm7d6rjmmmscjzzyiKOhocHhcDgc3bp1c3z99dcdXFn74C6hDjJ37lzNmDFDs2fP1jvvvKPS0lKVlpbqnXfe0ezZs/XAAw9o/vz5HV1mlxEeHq633npLTU1NzU7l5eUdXWKX8P3336tbt3P3AlgsFq1evVq//vWvNWbMGH377bcdWF3XdPaStYeHh3x8fBQYGOhc5u/vr5qamo4qrcu57rrrVFZWpqNHjyouLk5fffXVeUMKOjPuEuog//7v/66QkBD9/ve/1/PPP+88xe3p6amRI0dq7dq1uvvuuzu4yq5j5MiRKisr0/jx45tdbrFYutTtgx1l8ODB+uKLL8674+G5556TJP3zP/9zR5TVZVmtVu3Zs0dXX321JKmkpERXXXWVc3lFRYVzfBHaR8+ePVVQUKB169bJZrN1qcujjGExwJkzZ3Ts2DFJUkhIiLp3797BFXU9W7ZsUW1trZKSkppdXltbqy+++EJjxoxp58q6luzsbG3ZskXvvfdes8tnzZqlvLw8NTU1tXNlXVNeXp4iIyN1++23N7t8wYIFqq6u1osvvtjOlUGSDh06pLKyMtlsNvXo0aOjy2lzBBYAAGA8xrAAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMb7f/wg5tAAS/TjAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "recomb_cmp = pandas.DataFrame({'simulated': my_recomb_rates, 'estimated': est_recomb_rates})\n",
    "recomb_cmp.plot.bar()\n",
    "recomb_cmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "724aa103-82fa-4447-b469-db4d2cb28890",
   "metadata": {
    "id": "724aa103-82fa-4447-b469-db4d2cb28890",
    "outputId": "0fa7c7a4-5e32-4758-ec77-ab0c822202fd"
   },
   "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>simulated</th>\n",
       "      <th>estimated</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0001</td>\n",
       "      <td>1.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0001</td>\n",
       "      <td>1.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.0100</td>\n",
       "      <td>8.329049e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0010</td>\n",
       "      <td>1.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0001</td>\n",
       "      <td>2.151612e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.0001</td>\n",
       "      <td>1.000000e-08</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   simulated     estimated\n",
       "0     0.0001  1.000000e-08\n",
       "1     0.0001  1.000000e-08\n",
       "2     0.0100  8.329049e-03\n",
       "3     0.0010  1.000000e-08\n",
       "4     0.0001  2.151612e-03\n",
       "5     0.0001  1.000000e-08"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGYCAYAAAC6fxM1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwgElEQVR4nO3de1hVdb7H8c8GlY0XtgpyUxQqCh0NEhIhy06ROGNzRJtSpxR9PDrdPBpeHjFFnSw6moamRdZ4a8aDOSet44UynG5KmHhprDQ1DZ90o44jGCYoe58/fNydPSK5cSvy8/16nvUQa33Xb/3W2vjsT791szidTqcAAAAaOJ/67gAAAIA3EGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEZoVN8duFYcDocOHz6sFi1ayGKx1Hd3AADAZXA6nTp16pTCw8Pl41P7WMwNE2oOHz6siIiI+u4GAACog0OHDqldu3a11twwoaZFixaSzh+UgICAeu4NAAC4HOXl5YqIiHB9j9fmhgk1F045BQQEEGoAAGhgLufSES4UBgAARiDUAAAAIxBqAACAEW6Ya2oAAGaorq7W2bNn67sb8JLGjRvL19fXK20RagAADYLT6ZTdbtfJkyfruyvwspYtWyo0NPSKnyNHqAEANAgXAk1wcLCaNm3Kg1QN4HQ6dfr0aR09elSSFBYWdkXtEWoAANe96upqV6AJDAys7+7Ai/z9/SVJR48eVXBw8BWdiuJCYQDAde/CNTRNmzat557garjwuV7ptVKEGgBAg8EpJzN563Ml1AAAACPUKdQsWLBAkZGRslqtSkxM1JYtW2qtX7lypWJiYmS1WtWlSxetW7fObfk777yjXr16KTAwUBaLRTt27LiojTNnzuipp55SYGCgmjdvroceekilpaV16T4AAPVu6NChSktLu+rbmTZtmuLi4q76di7FYrFo9erV12RbHl8ovGLFCmVkZCg3N1eJiYnKyclRamqq9uzZo+Dg4IvqN2/erEGDBik7O1sPPvigli9frrS0NG3btk2dO3eWJFVUVKhHjx565JFHNGLEiBq3+8wzz2jt2rVauXKlbDabnn76afXv31+bNm3ydBcAAAaJnLj2mm7v4It9vNLO3Llz5XQ6vdKWt1ksFq1ateqahC5v8nikZs6cORoxYoSGDRumTp06KTc3V02bNtWiRYtqrJ87d6569+6t8ePHq2PHjnruuefUtWtXzZ8/31UzePBgZWVlKSUlpcY2ysrK9Kc//Ulz5szRfffdp/j4eC1evFibN2/W559/7ukuAABQ72w2m1q2bFnf3TCKR6GmqqpKxcXFbuHDx8dHKSkpKiwsrHGdwsLCi8JKamrqJetrUlxcrLNnz7q1ExMTo/bt21+yncrKSpWXl7tNAABca3/961/VpUsX+fv7KzAwUCkpKaqoqLjo9NO9996rUaNGacyYMWrVqpVCQkL0xhtvqKKiQsOGDVOLFi10yy23aP369a51lixZclEwWr16da0X3n7xxRd64IEHFBQUJJvNpp49e2rbtm2u5ZGRkZKkfv36yWKxuH6XpHfffVddu3aV1WrVTTfdpOnTp+vcuXOu5Xv37tU999wjq9WqTp06acOGDXU7aHXk0emn48ePq7q6WiEhIW7zQ0JCtHv37hrXsdvtNdbb7fbL3q7dbleTJk0u+uBqayc7O1vTp0+/7G0AprgaQ/HeGm4HbjRHjhzRoEGDNHPmTPXr10+nTp3Sp59+esnTTkuXLtWECRO0ZcsWrVixQk888YRWrVqlfv36adKkSXr55Zc1ePBglZSU1Pn29lOnTik9PV2vvPKKnE6nZs+erd/85jfau3evWrRooS+++ELBwcFavHixevfu7XpuzKeffqohQ4Zo3rx5uvvuu7V//36NHDlSkjR16lQ5HA71799fISEhKioqUllZmcaMGVOnPtaVsXc/ZWZmqqyszDUdOnSovrsEALjBHDlyROfOnVP//v0VGRmpLl266Mknn1Tz5s1rrI+NjdXkyZMVHR2tzMxMWa1WBQUFacSIEYqOjlZWVpb+8Y9/6Msvv6xzn+677z499thjiomJUceOHbVw4UKdPn1aH3/8sSSpTZs2kn5+dcGF36dPn66JEycqPT1dN910kx544AE999xzev311yVJH374oXbv3q1ly5YpNjZW99xzj1544YU697MuPBqpCQoKkq+v70V3HZWWlio0NLTGdUJDQz2qv1QbVVVVOnnypNtoTW3t+Pn5yc/P77K3AQCAt8XGxur+++9Xly5dlJqaql69eul3v/udWrVqVWP97bff7vpvX19fBQYGqkuXLq55F858XHitQF2UlpZq8uTJ+uijj3T06FFVV1fr9OnTKikpqXW9nTt3atOmTXr++edd86qrq3XmzBmdPn1a33zzjSIiIhQeHu5anpSUVOd+1oVHIzVNmjRRfHy8CgoKXPMcDocKCgou2fGkpCS3eknasGGDRzsaHx+vxo0bu7WzZ88elZSUXPMDBgDA5fL19dWGDRu0fv16derUSa+88opuu+02HThwoMb6xo0bu/1usVjc5l24VsbhcEg6f13rv57K+qWn8qanp2vHjh2aO3euNm/erB07digwMFBVVVW1rvfjjz9q+vTp2rFjh2v6+9//rr1798pqtda67rXi8S3dGRkZSk9PV0JCgrp166acnBzXRUySNGTIELVt21bZ2dmSpNGjR6tnz56aPXu2+vTpo7y8PG3dulULFy50tXnixAmVlJTo8OHDks4HFun8CE1oaKhsNpuGDx+ujIwMtW7dWgEBARo1apSSkpLUvXv3Kz4IAABcLRaLRXfddZfuuusuZWVlqUOHDlq1apVX2m7Tpo1OnTqliooKNWvWTJJqfNbb/7dp0ya9+uqr+s1vfiNJOnTokI4fP+5W07hxY1VXV7vN69q1q/bs2aNbbrmlxnY7duyoQ4cO6ciRI64XU17rO5Q9DjUDBgzQsWPHlJWVJbvdrri4OOXn57uGxEpKSuTj8/MAUHJyspYvX67Jkydr0qRJio6O1urVq13PqJGk9957zxWKJGngwIGSzl94NG3aNEnSyy+/LB8fHz300EOqrKxUamqqXn311TrtNAAA10JRUZEKCgrUq1cvBQcHq6ioSMeOHVPHjh2v6LqYCxITE9W0aVNNmjRJ//mf/6mioiItWbKk1nWio6P11ltvKSEhQeXl5Ro/frzrpZIXREZGqqCgQHfddZf8/PzUqlUrZWVl6cEHH1T79u31u9/9Tj4+Ptq5c6d27dqlGTNmKCUlRbfeeqvS09M1a9YslZeX69lnn73iffREnd7S/fTTT+vpp5+ucdlHH3100byHH35YDz/88CXbGzp0qIYOHVrrNq1WqxYsWKAFCxZ40lUAgOGu57vzAgIC9MknnygnJ0fl5eXq0KGDZs+erV//+tdasWLFFbffunVr/fnPf9b48eP1xhtv6P7779e0adNcdyXV5E9/+pNGjhyprl27KiIiQi+88ILGjRvnVjN79mxlZGTojTfeUNu2bXXw4EGlpqZqzZo1+uMf/6j/+q//UuPGjRUTE6P/+I//kHT+VNiqVas0fPhwdevWTZGRkZo3b5569+59xft5uSzO6/Vxhl5WXl4um82msrIyBQQE1Hd3gKuGW7phojNnzujAgQOKioq6bq7fgPfU9vl68v1t7C3dAADgxkKoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAEADMnToUKWlpdV3NzxyrfpMqAEA4Dp08OBBWSyWi15QOXfu3F98v5M3NMTwVKd3PwEAcN2YZrvG2yu7ttv7FzbbNd7fBoSRGgAAriKHw6Hs7GxFRUXJ399fsbGx+utf/ypJ+uc//6lHH31Ubdq0kb+/v6Kjo7V48WJJUlRUlCTpjjvukMVi0b333ivp4hGUe++9V6NGjdKYMWPUqlUrhYSE6I033lBFRYWGDRumFi1a6JZbbtH69etd61RXV2v48OGuPt12222aO3eua/m0adO0dOlSvfvuu7JYLLJYLK4XVh86dEiPPPKIWrZsqdatW6tv3746ePCgW9sZGRlq2bKlAgMDNWHCBF2r10wSagAAuIqys7O1bNky5ebm6quvvtIzzzyjxx57TB9//LGmTJmir7/+WuvXr9c333yj1157TUFBQZKkLVu2SJI+/PBDHTlyRO+8884lt7F06VIFBQVpy5YtGjVqlJ544gk9/PDDSk5O1rZt29SrVy8NHjxYp0+flnQ+aLVr104rV67U119/raysLE2aNElvv/22JGncuHF65JFH1Lt3bx05ckRHjhxRcnKyzp49q9TUVLVo0UKffvqpNm3apObNm6t3796qqqqSdP4N30uWLNGiRYv02Wef6cSJE1q1atXVPMQunH4CAOAqqays1AsvvKAPP/xQSUlJkqSbbrpJn332mV5//XX9+OOPuuOOO5SQkCBJioyMdK3bpk0bSVJgYKBCQ0Nr3U5sbKwmT54sScrMzNSLL76ooKAgjRgxQpKUlZWl1157TV9++aW6d++uxo0ba/r06a71o6KiVFhYqLfffluPPPKImjdvLn9/f1VWVrpt+89//rMcDofefPNNWSwWSdLixYvVsmVLffTRR+rVq5dycnKUmZmp/v37S5Jyc3P1/vvvX8lhvGyEGgAArpJ9+/bp9OnTeuCBB9zmV1VV6Y477tC0adP00EMPuUZT0tLSlJyc7PF2br/9dtd/+/r6KjAwUF26dHHNCwkJkSQdPXrUNW/BggVatGiRSkpK9NNPP6mqqkpxcXG1bmfnzp3at2+fWrRo4Tb/zJkz2r9/v8rKynTkyBElJia6ljVq1EgJCQnX5BQUoQYAgKvkxx9/lCStXbtWbdu2dVvm5+eniIgIff/991q3bp02bNig+++/X0899ZReeuklj7bTuHFjt98tFovbvAujKg6HQ5KUl5encePGafbs2UpKSlKLFi00a9YsFRUV/eL+xMfH6y9/+ctFyy6MLNUnQg0AAFdJp06d5Ofnp5KSEvXs2bPGmjZt2ig9PV3p6em6++67NX78eL300ktq0qSJpPMX3nrbpk2blJycrCeffNI1b//+/W41TZo0uWjbXbt21YoVKxQcHKyAgIAa2w4LC1NRUZHuueceSdK5c+dUXFysrl27enkvLkaoAQDgKmnRooXGjRunZ555Rg6HQz169FBZWZk2bdqkgIAA7d+/X/Hx8frVr36lyspKrVmzRh07dpQkBQcHy9/fX/n5+WrXrp2sVqvXbueOjo7WsmXL9P777ysqKkpvvfWWvvjiC9cdV9L563vef/997dmzR4GBgbLZbHr00Uc1a9Ys9e3bV3/84x/Vrl07ff/993rnnXc0YcIEtWvXTqNHj9aLL76o6OhoxcTEaM6cOTp58qRX+v1LuPsJAICr6LnnntOUKVOUnZ2tjh07qnfv3lq7dq2ioqLUpEkTZWZm6vbbb9c999wjX19f5eXlSTp/Lcq8efP0+uuvKzw8XH379vVan/7whz+of//+GjBggBITE/WPf/zDbdRGkkaMGKHbbrtNCQkJatOmjTZt2qSmTZvqk08+Ufv27dW/f3917NhRw4cP15kzZ1wjN2PHjtXgwYOVnp7uOrXVr18/r/W9Nhbntbp5vJ6Vl5fLZrOprKzskkNmgAkiJ671epsHX+zj9TYBT5w5c0YHDhxQVFSUrFZrfXcHXlbb5+vJ9zcjNQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQA0GDfIDbs3HG99roQaAMB178Ij/y+8ZRpmufC5/uvrHjzFE4UBANc9X19ftWzZ0vVCxqZNm7reZ4SGy+l06vTp0zp69KhatmwpX1/fK2qPUAMAaBBCQ0Mlub9pGmZo2bKl6/O9EoQaAECDYLFYFBYWpuDgYJ09e7a+uwMvady48RWP0FxAqAEANCi+vr5e+xKEWbhQGAAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACM0qu8OAGgAptmuQptl3m8TwA2NkRoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBHqFGoWLFigyMhIWa1WJSYmasuWLbXWr1y5UjExMbJarerSpYvWrVvnttzpdCorK0thYWHy9/dXSkqK9u7d61bz7bffqm/fvgoKClJAQIB69Oihv/3tb3XpPgAAMJDHoWbFihXKyMjQ1KlTtW3bNsXGxio1NVVHjx6tsX7z5s0aNGiQhg8fru3btystLU1paWnatWuXq2bmzJmaN2+ecnNzVVRUpGbNmik1NVVnzpxx1Tz44IM6d+6cNm7cqOLiYsXGxurBBx+U3W6vw24DAADTWJxOp9OTFRITE3XnnXdq/vz5kiSHw6GIiAiNGjVKEydOvKh+wIABqqio0Jo1a1zzunfvrri4OOXm5srpdCo8PFxjx47VuHHjJEllZWUKCQnRkiVLNHDgQB0/flxt2rTRJ598orvvvluSdOrUKQUEBGjDhg1KSUn5xX6Xl5fLZrOprKxMAQEBnuwy0KBETlzr9TYPWn/v9TZ5ojCAy+HJ97dHIzVVVVUqLi52CxE+Pj5KSUlRYWFhjesUFhZeFDpSU1Nd9QcOHJDdbnersdlsSkxMdNUEBgbqtttu07Jly1RRUaFz587p9ddfV3BwsOLj42vcbmVlpcrLy90mAABgLo9CzfHjx1VdXa2QkBC3+SEhIZc8DWS322utv/CzthqLxaIPP/xQ27dvV4sWLWS1WjVnzhzl5+erVatWNW43OztbNpvNNUVERHiyqwAAoIFpEHc/OZ1OPfXUUwoODtann36qLVu2KC0tTb/97W915MiRGtfJzMxUWVmZazp06NA17jUAALiWPAo1QUFB8vX1VWlpqdv80tJShYaG1rhOaGhorfUXftZWs3HjRq1Zs0Z5eXm666671LVrV7366qvy9/fX0qVLa9yun5+fAgIC3CYAAGAuj0JNkyZNFB8fr4KCAtc8h8OhgoICJSUl1bhOUlKSW70kbdiwwVUfFRWl0NBQt5ry8nIVFRW5ak6fPn2+sz7u3fXx8ZHD4fBkFwAAgKEaebpCRkaG0tPTlZCQoG7duiknJ0cVFRUaNmyYJGnIkCFq27atsrOzJUmjR49Wz549NXv2bPXp00d5eXnaunWrFi5cKOn89TJjxozRjBkzFB0draioKE2ZMkXh4eFKS0uTdD4YtWrVSunp6crKypK/v7/eeOMNHThwQH369PHSoQAAAA2Zx6FmwIABOnbsmLKysmS32xUXF6f8/HzXhb4lJSVuIyrJyclavny5Jk+erEmTJik6OlqrV69W586dXTUTJkxQRUWFRo4cqZMnT6pHjx7Kz8+X1WqVdP60V35+vp599lndd999Onv2rH71q1/p3XffVWxs7JUeAwAAYACPn1PTUPGcGtwoeE4NAJNctefUAAAAXK8INQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABghDqFmgULFigyMlJWq1WJiYnasmVLrfUrV65UTEyMrFarunTponXr1rktdzqdysrKUlhYmPz9/ZWSkqK9e/de1M7atWuVmJgof39/tWrVSmlpaXXpPgAAMJDHoWbFihXKyMjQ1KlTtW3bNsXGxio1NVVHjx6tsX7z5s0aNGiQhg8fru3btystLU1paWnatWuXq2bmzJmaN2+ecnNzVVRUpGbNmik1NVVnzpxx1fzP//yPBg8erGHDhmnnzp3atGmTfv/739dhlwEAgIksTqfT6ckKiYmJuvPOOzV//nxJksPhUEREhEaNGqWJEydeVD9gwABVVFRozZo1rnndu3dXXFyccnNz5XQ6FR4errFjx2rcuHGSpLKyMoWEhGjJkiUaOHCgzp07p8jISE2fPl3Dhw+v046Wl5fLZrOprKxMAQEBdWoDaAgiJ671epsHrVfhfyCmlXm/TQDG8eT726ORmqqqKhUXFyslJeXnBnx8lJKSosLCwhrXKSwsdKuXpNTUVFf9gQMHZLfb3WpsNpsSExNdNdu2bdMPP/wgHx8f3XHHHQoLC9Ovf/1rt9Gef1VZWany8nK3CQAAmMujUHP8+HFVV1crJCTEbX5ISIjsdnuN69jt9lrrL/ysrea7776TJE2bNk2TJ0/WmjVr1KpVK9177706ceJEjdvNzs6WzWZzTREREZ7sKgAAaGAaxN1PDodDkvTss8/qoYceUnx8vBYvXiyLxaKVK1fWuE5mZqbKyspc06FDh65llwEAwDXmUagJCgqSr6+vSktL3eaXlpYqNDS0xnVCQ0Nrrb/ws7aasLAwSVKnTp1cy/38/HTTTTeppKSkxu36+fkpICDAbQIAAObyKNQ0adJE8fHxKigocM1zOBwqKChQUlJSjeskJSW51UvShg0bXPVRUVEKDQ11qykvL1dRUZGrJj4+Xn5+ftqzZ4+r5uzZszp48KA6dOjgyS4AAABDNfJ0hYyMDKWnpyshIUHdunVTTk6OKioqNGzYMEnSkCFD1LZtW2VnZ0uSRo8erZ49e2r27Nnq06eP8vLytHXrVi1cuFCSZLFYNGbMGM2YMUPR0dGKiorSlClTFB4e7noOTUBAgB5//HFNnTpVERER6tChg2bNmiVJevjhh71xHAAAQAPncagZMGCAjh07pqysLNntdsXFxSk/P991oW9JSYl8fH4eAEpOTtby5cs1efJkTZo0SdHR0Vq9erU6d+7sqpkwYYIqKio0cuRInTx5Uj169FB+fr6sVqurZtasWWrUqJEGDx6sn376SYmJidq4caNatWp1JfsPAAAM4fFzahoqnlODGwXPqQFgkqv2nBoAAIDrFaEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxQp1CzYMECRUZGymq1KjExUVu2bKm1fuXKlYqJiZHValWXLl20bt06t+VOp1NZWVkKCwuTv7+/UlJStHfv3hrbqqysVFxcnCwWi3bs2FGX7gMAAAN5HGpWrFihjIwMTZ06Vdu2bVNsbKxSU1N19OjRGus3b96sQYMGafjw4dq+fbvS0tKUlpamXbt2uWpmzpypefPmKTc3V0VFRWrWrJlSU1N15syZi9qbMGGCwsPDPe02AAAwnMehZs6cORoxYoSGDRumTp06KTc3V02bNtWiRYtqrJ87d6569+6t8ePHq2PHjnruuefUtWtXzZ8/X9L5UZqcnBxNnjxZffv21e23365ly5bp8OHDWr16tVtb69ev1wcffKCXXnrJ8z0FAABG8yjUVFVVqbi4WCkpKT834OOjlJQUFRYW1rhOYWGhW70kpaamuuoPHDggu93uVmOz2ZSYmOjWZmlpqUaMGKG33npLTZs29aTbAADgBuBRqDl+/Liqq6sVEhLiNj8kJER2u73Gdex2e631F37WVuN0OjV06FA9/vjjSkhIuKy+VlZWqry83G0CAADmahB3P73yyis6deqUMjMzL3ud7Oxs2Ww21xQREXEVewgAAOqbR6EmKChIvr6+Ki0tdZtfWlqq0NDQGtcJDQ2ttf7Cz9pqNm7cqMLCQvn5+alRo0a65ZZbJEkJCQlKT0+vcbuZmZkqKytzTYcOHfJkVwEAQAPjUahp0qSJ4uPjVVBQ4JrncDhUUFCgpKSkGtdJSkpyq5ekDRs2uOqjoqIUGhrqVlNeXq6ioiJXzbx587Rz507t2LFDO3bscN0SvmLFCj3//PM1btfPz08BAQFuEwAAMFcjT1fIyMhQenq6EhIS1K1bN+Xk5KiiokLDhg2TJA0ZMkRt27ZVdna2JGn06NHq2bOnZs+erT59+igvL09bt27VwoULJUkWi0VjxozRjBkzFB0draioKE2ZMkXh4eFKS0uTJLVv396tD82bN5ck3XzzzWrXrl2ddx4AAJjD41AzYMAAHTt2TFlZWbLb7YqLi1N+fr7rQt+SkhL5+Pw8AJScnKzly5dr8uTJmjRpkqKjo7V69Wp17tzZVTNhwgRVVFRo5MiROnnypHr06KH8/HxZrVYv7CIAALgRWJxOp7O+O3EtlJeXy2azqaysjFNRMFrkxLVeb/Og9fdeb1PTyrzfJgDjePL93SDufgIAAPglhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMEKj+u4AAAAem2a7Cm2Web9NXFN1GqlZsGCBIiMjZbValZiYqC1bttRav3LlSsXExMhqtapLly5at26d23Kn06msrCyFhYXJ399fKSkp2rt3r2v5wYMHNXz4cEVFRcnf318333yzpk6dqqqqqrp0HwAAGMjjULNixQplZGRo6tSp2rZtm2JjY5WamqqjR4/WWL9582YNGjRIw4cP1/bt25WWlqa0tDTt2rXLVTNz5kzNmzdPubm5KioqUrNmzZSamqozZ85Iknbv3i2Hw6HXX39dX331lV5++WXl5uZq0qRJddxtAABgGovT6XR6skJiYqLuvPNOzZ8/X5LkcDgUERGhUaNGaeLEiRfVDxgwQBUVFVqzZo1rXvfu3RUXF6fc3Fw5nU6Fh4dr7NixGjdunCSprKxMISEhWrJkiQYOHFhjP2bNmqXXXntN33333WX1u7y8XDabTWVlZQoICPBkl4EGJXLiWq+3edD6e6+3yVA/rginn24Ynnx/ezRSU1VVpeLiYqWkpPzcgI+PUlJSVFhYWOM6hYWFbvWSlJqa6qo/cOCA7Ha7W43NZlNiYuIl25TOB5/WrVtfcnllZaXKy8vdJgAAYC6PQs3x48dVXV2tkJAQt/khISGy2+01rmO322utv/DTkzb37dunV155RX/4wx8u2dfs7GzZbDbXFBERUfvOAQCABq3B3dL9ww8/qHfv3nr44Yc1YsSIS9ZlZmaqrKzMNR06dOga9hIAAFxrHoWaoKAg+fr6qrS01G1+aWmpQkNDa1wnNDS01voLPy+nzcOHD+vf/u3flJycrIULF9baVz8/PwUEBLhNAADAXB6FmiZNmig+Pl4FBQWueQ6HQwUFBUpKSqpxnaSkJLd6SdqwYYOrPioqSqGhoW415eXlKioqcmvzhx9+0L333qv4+HgtXrxYPj4NbpAJAABcRR4/fC8jI0Pp6elKSEhQt27dlJOTo4qKCg0bNkySNGTIELVt21bZ2dmSpNGjR6tnz56aPXu2+vTpo7y8PG3dutU10mKxWDRmzBjNmDFD0dHRioqK0pQpUxQeHq60tDRJPweaDh066KWXXtKxY8dc/bnUCBEAALixeBxqBgwYoGPHjikrK0t2u11xcXHKz893XehbUlLiNoqSnJys5cuXa/LkyZo0aZKio6O1evVqde7c2VUzYcIEVVRUaOTIkTp58qR69Oih/Px8Wa1WSedHdvbt26d9+/apXbt2bv3x8I50AABgKI+fU9NQ8Zwa3Ch4Tg1uCDyn5oZx1Z5TAwAAcL0i1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiN6rsDAHC9ipy41uttHnyxj9fbBHAeIzUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEaoU6hZsGCBIiMjZbValZiYqC1bttRav3LlSsXExMhqtapLly5at26d23Kn06msrCyFhYXJ399fKSkp2rt3r1vNiRMn9OijjyogIEAtW7bU8OHD9eOPP9al+wAAwECNPF1hxYoVysjIUG5urhITE5WTk6PU1FTt2bNHwcHBF9Vv3rxZgwYNUnZ2th588EEtX75caWlp2rZtmzp37ixJmjlzpubNm6elS5cqKipKU6ZMUWpqqr7++mtZrVZJ0qOPPqojR45ow4YNOnv2rIYNG6aRI0dq+fLlV3gIrkzkxLVeb/Pgi3283mZD4O1jeaMeR+B65PV/31avNtdg8J1TO49HaubMmaMRI0Zo2LBh6tSpk3Jzc9W0aVMtWrSoxvq5c+eqd+/eGj9+vDp27KjnnntOXbt21fz58yWdH6XJycnR5MmT1bdvX91+++1atmyZDh8+rNWrV0uSvvnmG+Xn5+vNN99UYmKievTooVdeeUV5eXk6fPhw3fceAAAYw6ORmqqqKhUXFyszM9M1z8fHRykpKSosLKxxncLCQmVkZLjNS01NdQWWAwcOyG63KyUlxbXcZrMpMTFRhYWFGjhwoAoLC9WyZUslJCS4alJSUuTj46OioiL169fvou1WVlaqsrLS9XtZWZkkqby83JNd/kWOytNebU/yfh8bCm8fS46j95RbnF5vUw3g8+Hft/d4/d83f5Nec73/TV7on9P5y5+5R6Hm+PHjqq6uVkhIiNv8kJAQ7d69u8Z17HZ7jfV2u921/MK82mr+9dRWo0aN1Lp1a1fNv8rOztb06dMvmh8REXGp3btu2HLquwdm4Dh6j+1qNPriVWn1usffpXfwN+k9DeVv8tSpU7LZav+MPL6mpqHIzMx0GyFyOBw6ceKEAgMDZbFY6rFntSsvL1dERIQOHTqkgICA+u5Og8Vx9B6OpfdwLL2D4+g9DeFYOp1OnTp1SuHh4b9Y61GoCQoKkq+vr0pLS93ml5aWKjQ0tMZ1QkNDa62/8LO0tFRhYWFuNXFxca6ao0ePurVx7tw5nThx4pLb9fPzk5+fn9u8li1b1r6D15GAgIDr9g+sIeE4eg/H0ns4lt7BcfSe6/1Y/tIIzQUeXSjcpEkTxcfHq6CgwDXP4XCooKBASUlJNa6TlJTkVi9JGzZscNVHRUUpNDTUraa8vFxFRUWumqSkJJ08eVLFxcWumo0bN8rhcCgxMdGTXQAAAIby+PRTRkaG0tPTlZCQoG7duiknJ0cVFRUaNmyYJGnIkCFq27atsrOzJUmjR49Wz549NXv2bPXp00d5eXnaunWrFi5cKEmyWCwaM2aMZsyYoejoaNct3eHh4UpLS5MkdezYUb1799aIESOUm5urs2fP6umnn9bAgQMvazgKAACYz+NQM2DAAB07dkxZWVmy2+2Ki4tTfn6+60LfkpIS+fj8PACUnJys5cuXa/LkyZo0aZKio6O1evVq1zNqJGnChAmqqKjQyJEjdfLkSfXo0UP5+fmuZ9RI0l/+8hc9/fTTuv/+++Xj46OHHnpI8+bNu5J9vy75+flp6tSpF506g2c4jt7DsfQejqV3cBy9x7RjaXFezj1SAAAA1zne/QQAAIxAqAEAAEYg1AAAACMQagAAgBEINQCABod7XFATY1+T0BAcP35cixYtUmFhoesdVqGhoUpOTtbQoUPVpk2beu4hAFyf/Pz8tHPnTnXs2LG+u4LrCLd015MvvvhCqampatq0qVJSUlzP+SktLVVBQYFOnz6t999/3+3N5Ki7Q4cOaerUqVq0aFF9d+W699NPP6m4uFitW7dWp06d3JadOXNGb7/9toYMGVJPvWtYvvnmG33++edKSkpSTEyMdu/erblz56qyslKPPfaY7rvvvvru4nXv/7/D7/+bO3euHnvsMQUGBkqS5syZcy27ZYSKigq9/fbb2rdvn8LCwjRo0CDX8WyoCDX1pHv37oqNjVVubu5FL9h0Op16/PHH9eWXX6qwsLCeemiWnTt3qmvXrqqurq7vrlzXvv32W/Xq1UslJSWyWCzq0aOH8vLyXO9lKy0tVXh4OMfxMuTn56tv375q3ry5Tp8+rVWrVmnIkCGKjY2Vw+HQxx9/rA8++IBg8wt8fHwUGxt70bv7Pv74YyUkJKhZs2ayWCzauHFj/XSwAenUqZM+++wztW7dWocOHdI999yjf/7zn7r11lu1f/9+NWrUSJ9//rmioqLqu6t1RqipJ/7+/tq+fbtiYmJqXL57927dcccd+umnn65xzxqm9957r9bl3333ncaOHcuX8S/o16+fzp49qyVLlujkyZMaM2aMvv76a3300Udq3749ocYDycnJuu+++zRjxgzl5eXpySef1BNPPKHnn39ekpSZmani4mJ98MEH9dzT69uLL76ohQsX6s0333QLgI0bN9bOnTsvGk3Epfn4+Mhutys4OFiPPfaYDhw4oHXr1slms+nHH39Uv3791KZNGy1fvry+u1p3TtSLyMhI59KlSy+5fOnSpc4OHTpcuw41cBaLxenj4+O0WCyXnHx8fOq7m9e94OBg55dffun63eFwOB9//HFn+/btnfv373fa7XaO42UKCAhw7t271+l0Op3V1dXORo0aObdt2+Za/ve//90ZEhJSX91rULZs2eK89dZbnWPHjnVWVVU5nU6ns1GjRs6vvvqqnnvWsFgsFmdpaanT6XQ6b7rpJucHH3zgtnzTpk3OiIiI+uia13D3Uz0ZN26cRo4cqdGjR+u9995TUVGRioqK9N5772n06NF6/PHHNWHChPruZoMRFhamd955Rw6Ho8Zp27Zt9d3FBuGnn35So0Y/3z9gsVj02muv6be//a169uypb7/9th571/BcOLXs4+Mjq9Uqm83mWtaiRQuVlZXVV9calDvvvFPFxcU6duyYEhIStGvXrotO2+PyXDhuZ86ccZ1WvqBt27Y6duxYfXTLa7j7qZ489dRTCgoK0ssvv6xXX33VNZzv6+ur+Ph4LVmyRI888kg997LhiI+PV3Fxsfr27VvjcovFwi2glyEmJkZbt2696I6S+fPnS5L+/d//vT661SBFRkZq7969uvnmmyVJhYWFat++vWt5SUnJRV8quLTmzZtr6dKlysvLU0pKCqdA6+j+++9Xo0aNVF5erj179ri9XPr7779v8BcKE2rq0YABAzRgwACdPXtWx48flyQFBQWpcePG9dyzhmf8+PGqqKi45PJbbrlFf/vb365hjxqmfv366b//+781ePDgi5bNnz9fDodDubm59dCzhueJJ55w++L9/18ekrR+/XouEq6DgQMHqkePHiouLlaHDh3quzsNytSpU91+b968udvv//u//6u77777WnbJ67hQGAAAGIFragAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAI/wf3OolmYtsCdcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mut_cmp = pandas.DataFrame({'simulated': my_mut_rates, 'estimated': est_mut_rates})\n",
    "mut_cmp.plot.bar();\n",
    "mut_cmp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abc577f5-663d-4aa2-8a29-e1d40ef7a5a1",
   "metadata": {},
   "source": [
    "## Documentation\n",
    "All functions used in this notebook have a docstring documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5781dac8-bceb-4556-bd40-69283ecfe01b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function ped2graph in module recombulatorx.io:\n",
      "\n",
      "ped2graph(path)\n",
      "    This function takes a pedigree file as input and encodes it as graphs.\n",
      "    \n",
      "    The pedigree file can be a table in .ped format (headerless and space-separated),\n",
      "    in .tsv format (with header and tab-separated) or in .xlsx format (experimental).\n",
      "    The format is inferred by the extension.\n",
      "    \n",
      "    The graphs are DiGraphs objects from the networkx module, where each node is an\n",
      "    individual and edges go from parent to child. The nodes include individual information\n",
      "    (sex and genotypes) as attributes.\n",
      "    \n",
      "    IMPORTANT: markers must be sorted by genomic position, since we assume that recombination\n",
      "    happens between adjacent markers.\n",
      "    \n",
      "    path: str\n",
      "        path to the pedigree file\n",
      "    \n",
      "    Returns: tuple\n",
      "        a list of tuples (family ID, DiGraphs) for each family ID in the pedigree file and \n",
      "        a list of marker names.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(recombulatorx.ped2graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d6198a59-1fb0-41ea-95ff-741bc1fb37fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function plot_family_graph in module recombulatorx.families:\n",
      "\n",
      "plot_family_graph(G, title=None)\n",
      "    Plot a family graph.\n",
      "    \n",
      "    Requires the matplotlib module.\n",
      "    \n",
      "    G: networkx.DiGraph\n",
      "        family graph\n",
      "    title: str\n",
      "        plot title\n",
      "    \n",
      "    Returns: matplotlib.Figure\n",
      "        the family plot\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(recombulatorx.families.plot_family_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "39cd8626-ec03-4b99-a75b-4caf14ccc871",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function preprocess_families in module recombulatorx.families:\n",
      "\n",
      "preprocess_families(family_graphs)\n",
      "    Preprocess family graphs for rates estimation.\n",
      "    \n",
      "    Does the following preprocessing steps:\n",
      "    - checks the consistency of the parsed family graphs,\n",
      "    - finds all the informative subfamilies from each family,\n",
      "    - phases all females using their father when available,\n",
      "    - packs the inherited maternal haplotypes into a matrix\n",
      "    \n",
      "    family_graphs: sequence of tuples (str, networkx.DiGraph)\n",
      "        sequence of tuples of family IDs and family graphs\n",
      "    \n",
      "    Returns: list of ProcessedFamily\n",
      "        the preprocessed families\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(recombulatorx.preprocess_families)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "359c6ee0-ac75-42f5-86ed-a45c3529df70",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function estimate_rates in module recombulatorx.estimate:\n",
      "\n",
      "estimate_rates(families, starting_recombination_rates=0.05, starting_mutation_rates=0.001, estimate_mutation_rates='no', implementation=None, optimization_method='L-BFGS-B', maxiter=1000)\n",
      "    Estimate recombination and optionally mutation rates from a set of families.\n",
      "    \n",
      "    Estimation is done by expectation maximization, that is by finding recombination \n",
      "    and mutation rates that maximize the likelihood of observing the given families.\n",
      "    Implemented using scipy.optimize.minimize to find recombination and mutation rates \n",
      "    that minimize the negative log likelihood.\n",
      "    \n",
      "    families: sequence of families\n",
      "        the preprocessed genetic data of all the families\n",
      "    starting_recombination_rates: array of shape (n_markers - 1, )\n",
      "        starting probability of recombination between adjacent markers as initial \n",
      "    starting_mutation_rates: array of shape (n_markers, )\n",
      "        probability of unit mutation (+1 or -1) for each marker\n",
      "    estimate_mutation_rates: string\n",
      "        one of 'no', 'one', 'all'. If 'no' \n",
      "    implementation: string\n",
      "        the name of the implementation to use for the likelihood computation\n",
      "    optimization_method: string\n",
      "        the name of the method to be passed to scipy.optimize.minimize\n",
      "    maxiter: int\n",
      "        the maximum number of optimizing iterations (passed to scipy.optimize.minimize)\n",
      "      \n",
      "    Returns: arrays\n",
      "        an array of estimated recombination rates if estimate_mutation_rates is 'no',\n",
      "        a tuple of an array of estimated recombination rates and a scalar mutation rate if estimate_mutation_rates is 'one',\n",
      "        a tuple of an array of estimated recombination rates and an array of mutation rates if estimate_mutation_rates is 'all'\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(recombulatorx.estimate_rates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6089e59b-3626-45a1-b66e-8cbfc57cd24b",
   "metadata": {
    "id": "6089e59b-3626-45a1-b66e-8cbfc57cd24b",
    "outputId": "57e144ab-f130-4942-d54d-d86a8462142f"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function generate_complex_families in module recombulatorx.testing:\n",
      "\n",
      "generate_complex_families(n_fam_I: int, n_fam_II: int, recombination_rates, mutation_rates)\n",
      "    Generate randomly simulated individuals from a given number of type I and II families.\n",
      "    \n",
      "    For each family, generate a random number of children (from a Poisson distribution with a\n",
      "    mean of 2.0), with a minimum of 1 or 2 for type I or II families, respectively. A random \n",
      "    number of children (from a binomial distribution with a mean of 0.5) are set to be male and \n",
      "    the rest female. For each daughter, a father and his paternal haplotype are generated randomly.\n",
      "    \n",
      "    n_fam_I: int\n",
      "        number of type I families\n",
      "    n_fam_II: int\n",
      "        number of type II families\n",
      "    recombination_rates: array of shape (n_markers - 1, )\n",
      "        recombination rates between adjacent markers used to generate children\n",
      "    mutation_rates: array of shape (n_markers, )\n",
      "        mutation rates of markers used to generate children\n",
      "    \n",
      "    Returns:\n",
      "        an iterator to the generated individuals\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(recombulatorx.testing.generate_complex_families)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8cf50660-a283-4a93-b6d3-b5fa05754bbb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function individuals2ped in module recombulatorx.testing:\n",
      "\n",
      "individuals2ped(path, marker_names, individuals)\n",
      "    Write individuals to a pedigree file.\n",
      "    \n",
      "    If path ends with `.tsv' write a tab-separated table with an header row,\n",
      "    otherwise a headerless space-separated table in `.ped' style.\n",
      "    \n",
      "    path: str\n",
      "        path to the pedigree file\n",
      "    marker_names: sequence of str\n",
      "        names for the markers\n",
      "    individuals: sequence\n",
      "        individuals as generated by generate_complex_families\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(recombulatorx.testing.individuals2ped)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "venv",
   "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.9.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "ed4a9856b4a60b97a7fb96513744f09baa8490c3ea00c622e8f393bc9b623388"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}