{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bf128f9e",
   "metadata": {},
   "source": [
    "# Tutorial\n",
    "\n",
    "In this section, the user can find a tutorial on how to use `EpiK` for inference of sequence-function relationships on a small dataset, how to interpret the inferred hyperparameters of the different kernels, and how to use `EpiK` to make calibrated phenotypic predictions in unobserved sequences."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5909eeca",
   "metadata": {},
   "source": [
    "### Import required libraries and functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "58f88b18",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[KeOps] Warning : Cuda libraries were not detected on the system or could not be loaded ; using cpu only mode\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from scipy.stats import pearsonr\n",
    "from epik.utils import encode_seqs, split_training_test\n",
    "from epik.model import EpiK\n",
    "from epik.kernel import ConnectednessKernel, ExponentialKernel, JengaKernel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd35935",
   "metadata": {},
   "source": [
    "### Load SMN1 data\n",
    "\n",
    "We will use a previously published dataset, in which the activity of nearly every possible GU/GC 5´splice site sequence was measured in parallel in a single experiment ([Wong et al. 2018](https://pubmed.ncbi.nlm.nih.gov/30174293/)). The measured phenotype `y` is the percent spliced in, this is, the percentage of times exon 7 is included in the mature mRNA. The assay was performed in different biological replicates, providing a measure of uncertainty of our estimates given by `y_var`, as previously analyzed ([Zhou et al. 2022](https://pubmed.ncbi.nlm.nih.gov/36129941/))\n",
    "\n",
    "We start by loading the data from the pre-processed CSV file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6a86484d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>y_var</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>seq</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AAACAAAA</th>\n",
       "      <td>0.223881</td>\n",
       "      <td>0.001738</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAACAAAC</th>\n",
       "      <td>0.514961</td>\n",
       "      <td>0.009195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAACAAAG</th>\n",
       "      <td>0.208009</td>\n",
       "      <td>0.001500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAACAAAU</th>\n",
       "      <td>0.235378</td>\n",
       "      <td>0.002109</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAACAACA</th>\n",
       "      <td>0.357162</td>\n",
       "      <td>0.004423</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UUUUUUGU</th>\n",
       "      <td>0.199017</td>\n",
       "      <td>0.001373</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UUUUUUUA</th>\n",
       "      <td>0.074571</td>\n",
       "      <td>0.000536</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UUUUUUUC</th>\n",
       "      <td>0.193181</td>\n",
       "      <td>0.001294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UUUUUUUG</th>\n",
       "      <td>0.121894</td>\n",
       "      <td>0.000536</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UUUUUUUU</th>\n",
       "      <td>0.152185</td>\n",
       "      <td>0.000803</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>30732 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 y     y_var\n",
       "seq                         \n",
       "AAACAAAA  0.223881  0.001738\n",
       "AAACAAAC  0.514961  0.009195\n",
       "AAACAAAG  0.208009  0.001500\n",
       "AAACAAAU  0.235378  0.002109\n",
       "AAACAACA  0.357162  0.004423\n",
       "...            ...       ...\n",
       "UUUUUUGU  0.199017  0.001373\n",
       "UUUUUUUA  0.074571  0.000536\n",
       "UUUUUUUC  0.193181  0.001294\n",
       "UUUUUUUG  0.121894  0.000536\n",
       "UUUUUUUU  0.152185  0.000803\n",
       "\n",
       "[30732 rows x 2 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('../data/smn1.csv', index_col=0)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a23b05c9",
   "metadata": {},
   "source": [
    "### Prepare data for model fitting\n",
    "\n",
    "`EpiK` models assume that sequences are provided in a one-hot encoding `X`, so we provide the utility function `encode_seqs` to generate this encoding from the raw sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0e1ae656",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = encode_seqs(data.index.values, alphabet='ACGU')\n",
    "y = data['y'].values\n",
    "y_var = data['y_var'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc8f2ae5",
   "metadata": {},
   "source": [
    "We next split the data into training and test subsets with `split_training_test`. To keep this tutorial as simple and fast to run as possible, we will use only a small subset of the dataset for both training and testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6c25c176",
   "metadata": {},
   "outputs": [],
   "source": [
    "splits = split_training_test(X, y, y_var, ptrain=0.075, dtype=torch.float32)\n",
    "train_x, train_y, test_x, test_y, train_y_var = splits\n",
    "\n",
    "# Subsample test data\n",
    "idx = np.random.uniform(size=test_y.shape[0]) < 0.05\n",
    "test_x, test_y = test_x[idx, :], test_y[idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83dbed32",
   "metadata": {},
   "source": [
    "### Fit Connectedness model regression\n",
    "\n",
    "To define an `EpiK` model, we first need to specify the kernel function that we want to use for Gaussian process regression. In this case, we will use the `ConnectednessKernel`, in which the predictability of mutational effects decreases depending on the sites at which two sequences differ. \n",
    "\n",
    "As the number of parameters is equal to sequence length, this kernel provides a good compromise in datasets with relatively few measurements, such as this one. For every kernel class, we need to specify the configuration of sequence space, consisting of the number of alleles `n_alleles` and sequence length `seq_length`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "070f6a1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = ConnectednessKernel(n_alleles=4, seq_length=8)\n",
    "model = EpiK(kernel, track_progress=True)\n",
    "model.set_data(train_x, train_y, train_y_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b9f96ad",
   "metadata": {},
   "source": [
    "After defining the `EpiK` model, we load the training data in the model and optimize the kernel hyperparameters by maximizing the marginal likelihood using the Adam optimizer for a number of predefined iterations. The argument `track_progress=True` allows us to follow the progress of the optimization and track how the loss function decreases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3c4d89c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Optimizing hyperparameters: 100%|██████████| 100/100 [00:14<00:00,  6.91it/s, MLL=-8918.879, Mem(alloc/res)=0.00/0.00MB]\n"
     ]
    }
   ],
   "source": [
    "model.fit(n_iter=100, learning_rate=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a5b5ee7",
   "metadata": {},
   "source": [
    "### Make predictions on test sequences\n",
    "\n",
    "Once we have optimized the hyperparameters of the kernel function, we can make predictions by computing the posterior distribution over the sequence-function relationship evaluated at sequences of interest. This allows us to obtain point estimates `test_y_pred` that maximize this posterior probability. However, we can additionally compute the posterior variance `test_y_var` quantifying the uncertainty of the predictions by setting `calc_variance=True`. \n",
    "\n",
    "Lets now make predictions in the held-out test sequences. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2771e770",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>stderr</th>\n",
       "      <th>lower_ci</th>\n",
       "      <th>upper_ci</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.116530</td>\n",
       "      <td>12.501842</td>\n",
       "      <td>-22.887154</td>\n",
       "      <td>27.120216</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.571670</td>\n",
       "      <td>12.174002</td>\n",
       "      <td>-22.776333</td>\n",
       "      <td>25.919674</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.843107</td>\n",
       "      <td>12.838286</td>\n",
       "      <td>-26.519680</td>\n",
       "      <td>24.833466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.053030</td>\n",
       "      <td>12.459997</td>\n",
       "      <td>-21.866964</td>\n",
       "      <td>27.973024</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-2.211174</td>\n",
       "      <td>12.297050</td>\n",
       "      <td>-26.805273</td>\n",
       "      <td>22.382925</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1492</th>\n",
       "      <td>0.437592</td>\n",
       "      <td>12.282001</td>\n",
       "      <td>-24.126410</td>\n",
       "      <td>25.001593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1493</th>\n",
       "      <td>-2.043804</td>\n",
       "      <td>12.466360</td>\n",
       "      <td>-26.976524</td>\n",
       "      <td>22.888916</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1494</th>\n",
       "      <td>-3.416660</td>\n",
       "      <td>12.331523</td>\n",
       "      <td>-28.079706</td>\n",
       "      <td>21.246386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1495</th>\n",
       "      <td>-0.705414</td>\n",
       "      <td>12.133599</td>\n",
       "      <td>-24.972612</td>\n",
       "      <td>23.561785</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1496</th>\n",
       "      <td>1.321175</td>\n",
       "      <td>12.068440</td>\n",
       "      <td>-22.815706</td>\n",
       "      <td>25.458055</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1497 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          coef     stderr   lower_ci   upper_ci\n",
       "0     2.116530  12.501842 -22.887154  27.120216\n",
       "1     1.571670  12.174002 -22.776333  25.919674\n",
       "2    -0.843107  12.838286 -26.519680  24.833466\n",
       "3     3.053030  12.459997 -21.866964  27.973024\n",
       "4    -2.211174  12.297050 -26.805273  22.382925\n",
       "...        ...        ...        ...        ...\n",
       "1492  0.437592  12.282001 -24.126410  25.001593\n",
       "1493 -2.043804  12.466360 -26.976524  22.888916\n",
       "1494 -3.416660  12.331523 -28.079706  21.246386\n",
       "1495 -0.705414  12.133599 -24.972612  23.561785\n",
       "1496  1.321175  12.068440 -22.815706  25.458055\n",
       "\n",
       "[1497 rows x 4 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred = model.predict(test_x, calc_variance=True)\n",
    "pred"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83f8a8a9",
   "metadata": {},
   "source": [
    "Once we have the predictions, we can compare them with the measurements for these held-out sequences. We can evaluate the accuracy of the point estimates by computing the $R^2$ and uncertainty calibration by calculating the percentage of times the measurement is within the 95% posterior predictive interval. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ef545c87",
   "metadata": {},
   "outputs": [],
   "source": [
    "r2 = pearsonr(pred['coef'], test_y)[0] ** 2\n",
    "sd = pred['stderr']\n",
    "obs = test_y.numpy()\n",
    "calibration = ((obs > pred['coef'] - 2 * sd) & (obs < pred['coef'] + 2 * sd)).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3189342f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0.5, 0, 'Predicted PSI (%)'),\n",
       " (-20.0, 200.0),\n",
       " Text(0, 0.5, 'Measured PSI (%)'),\n",
       " (-20.0, 200.0),\n",
       " None]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAFSCAYAAADBzXWtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABk0UlEQVR4nO2deVhT1/b3v0mAECAkMgYUcQQcGLUidVbqWIrVDlp7Heq1rdVatbbWtmq1Vlvbawfrre3v1tre62wdqu21dSiOQAuIAyoOVRElOACJDIHkZL9/8OZcwpyQkEDW53ny6DlnZ5+1kvDNztprry1gjDEQBEEQdoPQ1gYQBEEQxpAwEwRB2BkkzARBEHYGCTNBEISdQcJMEARhZ5AwEwRB2BkkzARBEHYGCTNBEISdQcJMEARhZ5AwEwRB2Bk2FeZVq1bhkUcegVQqhZ+fH8aOHYvs7GyjNhqNBrNmzYK3tzc8PDwwfvx45OfnG7XJycnBmDFj4ObmBj8/P7zxxhvQ6XTN6QpBEITFsKkwHz16FLNmzUJKSgoOHjwIrVaL4cOHo6SkhG8zb9487Nu3Dzt27MDRo0dx584djBs3jr/OcRzGjBmDiooKnDp1Ct9//z02btyIJUuW2MIlgiCIpsPsiLt37zIA7OjRo4wxxoqKipizszPbsWMH3+bixYsMAEtOTmaMMfbLL78woVDIlEol3+arr75inp6erLy8vHkdIAiCsABOtv5iqIpKpQIAeHl5AQDS09Oh1WoRHx/PtwkLC0P79u2RnJyMvn37Ijk5GeHh4fD39+fbjBgxAjNnzkRWVhaio6Nr3Ke8vBzl5eX8sV6vR0FBAby9vSEQCKzlHkEQrQTGGB4+fIjAwEAIhZYPPNiNMOv1esydOxf9+vVDz549AQBKpRIuLi6Qy+VGbf39/aFUKvk2VUXZcN1wrTZWrVqFZcuWWdgDgiAcjVu3bqFdu3YW79duhHnWrFk4f/48Tpw4YfV7LVq0CPPnz+ePVSoV2rdvj5ycHHh6elr8fowxqFQqyGQyhxmRO6LPgGP67Ug+a7QcUq8XoKevEzp37AipVGqV+9iFMM+ePRv79+/HsWPHjL59FAoFKioqUFRUZDRqzs/Ph0Kh4Nv88ccfRv0ZsjYMbaojFoshFotrnJfJZFYTZkP/rf2Da8ARfQYc029H8Vmj5bA/Ixd+Uim85JX6YS1/bZqVwRjD7NmzsXv3bhw5cgQdO3Y0ut6rVy84Ozvj8OHD/Lns7Gzk5OQgLi4OABAXF4dz587h7t27fJuDBw/C09MT3bt3bx5HCIJo1ZTrOPyYkQs/qSviu/lZ/QvIpiPmWbNmYfPmzdi7dy+kUikfE5bJZJBIJJDJZJg+fTrmz58PLy8veHp64tVXX0VcXBz69u0LABg+fDi6d++Ov/3tb1i9ejWUSiXeffddzJo1q9ZRMUEQhKk4C4WIbCdHj0DP5vlVYKt0EFb5+6fWx3fffce3KSsrY6+88gpr06YNc3NzY08++STLy8sz6ufGjRts1KhRTCKRMB8fH/b6668zrVbbaDtUKhUDwFQqlaVcM0Kv17PCwkKm1+ut0r894og+M+aYfrdmn8sqdOynzNuspNxYT4qKiqyqGQLGaDNWtVoNmUwGlUpFk38WwhF9BhzT79bqs0Zbd/hCpVJBLpdbTTOoVkYzkZubiyFDhqB79+6IiIjAjh07bG0SQRB1oOP02JVxu9liytWxi6wMR8DJyQmffvopoqOjoVQq0atXL4wePRru7u62No0giCowxuAkEuLRzt4I9nazya8AGjE3EwqFAlFRUfz/fXx8UFBQYFujCIIwQqPlsDM9FwUlFejg426z0AwJs4UYNGgQBAIBBAIBXFxc0K1bN2zevLnWtunp6eA4DkFBQWbfb926dejQoQNcXV0RGxtbI5e7Nt577z3eRsMjLCyMv85xHBYvXoyOHTtCIpGgc+fOeP/990HTEIQjYIgpy91c0MbN2aa2kDBbAMYYTp8+jU8++QR5eXnIzs7GyJEjMXnyZFy/ft2obUFBASZPnoxvvvnG7Ptt27YN8+fPx9KlS5GRkYHIyEiMGDHCKJe7Lnr06IG8vDz+UXWl5UcffYSvvvoKX375JS5evIiPPvoIq1evxtq1a822lSBaAowx7M20XUy5OiTMFuDKlSt4+PAhRo4cCYVCgY4dO2L69OngOM6ovnR5eTnGjh2Lt956C48++qjZ91uzZg1mzJiBadOmoXv37li/fj3c3NywYcOGBp/r5OQEhULBP3x8fPhrp06dQmJiIsaMGYMOHTrgqaeewvDhwxs1GieIlopezyAQCDAkzM8uRBkgYbYI6enpaNOmDb/SMDc3F++88w7EYjEiIiIAVH4jT5s2DUOHDsXf/va3Gn2sXLkSHh4e9T5ycnJQUVGB9PR0o4p7QqEQ8fHxSE5ObtDWK1euIDAwEJ06dcKkSZOQk5PDX3v00Udx+PBhXL58GQBw5swZnDhxAqNGjWrS60MQ9opGy2Hrn7eQW1gKP6mrXYgyQFkZFiEjIwMqlQpSqRQcx0Gj0UAikWD9+vUIDAwEYwwpKSnYtm0bIiIisGfPHgDAv//9b4SHhwMAXn75ZTzzzDP13icwMBB3794Fx3G1VtS7dOlSvc+PjY3Fxo0bERoairy8PCxbtgwDBgzA+fPnIZVK8dZbb0GtViMsLAwikQgcx+GDDz7ApEmTzH9xCMJOqZqn3FYusbU5RpAwW4CMjAzMmjULc+bMQVFRERYsWIB+/fph6tSpfJu4uDhwHFfnN7KXlxdfh9paVB35RkREIDY2FsHBwdi+fTumT5+O7du3Y9OmTdi8eTN69OiBzMxMzJ07F4GBgZgyZYpVbSOI5ubAeaXdxJSrQ8JsATIyMjBjxgx06dIFAPDPf/4TERERmDFjBjp06NCoPlauXImVK1fW2+bChQtQKBQQiUQ19j2sWnGvscjlcoSEhODq1asAgDfeeANvvfUWJkyYAAAIDw/HzZs3sWrVKhJmotVQodPDWSTA0G5+kIqd7E6UAYoxN5m//voLRUVFfHF/AOjevTs6d+5cZ7pcbbz88svIzMys9xEYGAgXFxf06tXLqOKeXq/H4cOH+Yp7jaW4uBjXrl1DQEAAAKC0tLTGbgwikQh6vd6kfgnCXtFoOexIv4XL+cXwdHW2S1EGaMTcZNLT0+Hs7IyQkBCj88OGDcPu3bvx9ttvN6ofU0IZ8+fPx5QpU9C7d2/06dMHn332GUpKSjBt2jS+zZdffondu3cbCfiCBQuQkJCA4OBg3LlzB0uXLoVIJMLEiRMBAAkJCfjggw/Qvn179OjRA6dPn8aaNWvwwgsvNMougrBnqsaUQ/w9bG1OvZAwN5GMjAx07doVLi4uRufj4+Oxfv165Obmom3btha957PPPot79+5hyZIlUCqViIqKwoEDB4wmBO/fv49r164ZPS83NxcTJ07EgwcP4Ovri/79+yMlJQW+vr4AgLVr12Lx4sV45ZVXcPfuXQQGBuKll16iHceJVsHRy/fsNqZcHaouB6ouZw0c0WfAMf22d581Wg4ioQB6xuAiElrERqouRxAEYSaG8MW52yqInUR2+cVRGyTMBEG0SqrGlKOD5LY2xyRImAmCaJWk3yxsMTHl6tDkH0EQrQqNloNOz9C3kzeEAuvtZG1NaMRMEESrwRC+OJtbBJFQ0CJFGSBhJgiilVA1phzXydvW5jQJmwrzsWPHkJCQgMDAQAgEAr64j4HqRd0Nj48//phv06FDhxrXP/zww2b2pOls3LgRcrmcP37vvff4HU8AYOrUqRg7dmyz2tShQwd89tlnzXpPgjCXK/nFLTamXB2bCnNJSQkiIyOxbt26Wq9XLeiel5eHDRs2QCAQYPz48Ubtli9fbtTu1VdfbQ7zeZRKJV599VV06tQJYrEYQUFBSEhIMFp1ZyoLFizAb7/9ZkEr66b6l4KBP//8Ey+++GKz2NAUzp49iwEDBsDV1RVBQUFYvXp1g885fPgwHn30UUilUigUCixcuBA6nc6oDWMMn3zyCUJCQiAWi9G2bVt88MEH1nKDMBONlsNdtQbh7WStQpQBG0/+jRo1qt5av9WL8uzduxdDhgxBp06djM4b/rhswY0bN9CvXz/I5XJ8/PHHCA8Ph1arxa+//opZs2Y1WIqzLg4dOoTMzExERUWZPVKuqKiosSLRFAwrAu0ZtVqN4cOH8ystz507hxdeeAFyubzOL5UzZ85g9OjReOedd/DDDz/g9u3bePnll8FxHD755BO+3WuvvYbffvsNn3zyCcLDw1FQUED7NNoZGi2HXRm3ESh3hZ+n/dRTbjLMTgDAdu/eXed1pVLJnJyc2KZNm4zOBwcHM39/f+bl5cWioqLY6tWrmVarrfdeGo2GqVQq/nHr1i0GgBUVFTG9Xm/SY9SoUaxt27bs4cOHNa4VFBQwvV7POI5j77//PuvZsydzc3Nj7dq1Yy+//DJTq9V82w0bNjCZTMb0ej3TarVs0KBBzN/fny1dupRptVo2ZcoUlpiYyJYuXcp8fHyYVCplL774ItNoNHwfgwYNYq+88gqbM2cO8/b2ZoMHD2Z6vZ598skndd77yJEjDIDRY8mSJUyv17Pg4GC2Zs0avv8bN26wJ554grm7uzOpVMqefvpplpeXx19fsmQJi4yMZN9//z0LDg5mUqmUPfPMM0ylUpn8ujb2sW7dOtamTRuj1+HNN99koaGhdT7nrbfeYr179zY6t3fvXubq6srbmpWVxZycnNjFixdNsofjOFZQUMA4jrOaz/b2sJXPpeVa9p+UG+zX83nNfu/CwkIGgKlUKkvIXw1aTLrc999/D6lUinHjxhmdnzNnDmJiYuDl5YVTp05h0aJFyMvLw5o1a+rsa9WqVVi2bFmN8yqVyqSNRwsLC3HgwAG8++670Ol0UKlURtcFAgHfp1arxQcffIAOHTrgxo0bWLBgAebOnYt//OMfAICysjJ+aStQ+Wvh8uXLiIqKQnFxMSoqKnD48GEIhUL89NNPyMnJwezZs+Hh4YHFixcDAHQ6HX744QdMmzYN//3vf3mfKioqsHLlSgQHB9e4d48ePbBq1SqsXLkSf/75JwDA3d0dKpUKer0eGo2G/39CQgLc3d2xf/9+6HQ6vPHGG3jqqaewf/9+AJVbZ127dg07d+7Epk2bkJeXh9mzZ2PZsmW8jdW5detWg1Xx5s2bh9dff73Wa8eOHUNcXBzKyspQVlYGAOjfvz9Wr16Nmzdv1hqiefjwIZycnIzeL4Ovx44dQ//+/bFjxw506NABO3fuxP/93/+BMYbBgwdj2bJlaNOmTZ22MsZQXFwMoGWmaZmDrXzOKSyDVMThkUAx1Gp1s90XQI2/dYtjFbk3AzQwYg4NDWWzZ89usJ9vv/2WOTk5MY1GU2cbS42YU1JSGAD2448/mjyi2L59O/P29uaPq46Y9frK0WdERAR/PGXKFObl5cWKi4v5c//85z+Zh4cH0+l0TK+vHDFHR0c3aHdD9zY8qo6Yf/31VyYSidjNmzf56+fPn2cAWGpqKm+zm5sbU6lUvM8LFixgsbGxddpSUVHBLl++XO/j/v37dT7/scceYzNmzDA6Z7ArKyur1uccOHCACYVCtmnTJqbVatmtW7fYgAEDGAC2adMmptfr2YsvvsjEYjGLjY1lR48eZUeOHGFRUVFsyJAhJr/Xrf3R3D6XlmvZ1Xx1s9yrrgeNmAEcP34c2dnZ2LZtW4NtY2NjodPpcOPGDYSGhtbaRiwWQywW1zhvyOowlcY87+jRo1i7di0uXboEtVoNnU4HjUaDsrIyuLm58c+v+m/1fiMjI+Hu7s4fP/rooyguLkZubi6Cg4MBAL169aphy6FDh7Bq1apG37s23y5duoSgoCC0b9+ev9ajRw/I5XJcunQJffr0gUAgQIcOHeDp6QnGKje4DAgIwN27d+t8fZydndG1a9d6X7uGqP46VX8NqzNixAh8/PHHmDlzJiZPngyxWIzFixfj+PHjEIkq6ykwxlBeXo4ffviBL+n67bffolevXrh8+XKdn62q93WUETPQfD5rtBx2nb4Nf6krOvtJrXqv+rC2ny0ij9nwBxEZGdlg28zMTAiFQvj5+Vndrq5du/KiVR83btzAhAkTEB4ejh9//BHp6el8JkpFRYVFbaoq3IZ7P/7444iIiLD6vYFKoa2KQCCot9B+Tk5Og5vQ1rezi0KhqHU3F8O1upg/fz6KioqQk5OD+/fvIzExEQD4ieWAgAA4OTkZ1dnu1q0bbzPR/FTNUx7Wzfp/37bEpiPm4uJiflsjALh+/ToyMzPh5eXFj8zUajV27NjBx2KrkpycjNTUVAwZMgRSqRTJycmYN28enn/++XrjgJbCy8sLI0aMwLp16zBnzpwaolhUVAS5XI709HTo9Xr84x//gEgkAgBs377d5PudOXMGZWVlkEgqN45MSUmBh4cHgoKC6nxO1Xsbdiepfm8XFxdwHFfvvbt164Zbt27h1q1b/P0uXLiAoqIifndwcwgMDERmZma9berbQCAuLg7vvPMOtFot/6Vw8OBBhIaGNvgZEAgECAwMBABs2bIFQUFBiImJAQD069cPOp0O165dQ+fOnQGA3z3c8OuEaF4eanRo18YNA7v6tP5fI1YJkDSS33//vUZGAAA2ZcoUvs3XX3/NJBIJKyoqqvH89PR0Fhsby2QyGXN1dWXdunVjK1eurDe+XBsqlcrseNG1a9eYQqFg3bt3Zzt37mSXL19mFy5cYJ9//jkLCwtjjDF2+vRpBoB9+umn7Nq1a+yHH35gbdu2ZQBYYWEhY4yx7777jslkMr7fpUuXssjISP54ypQpzMPDg02cOJFlZWWxn3/+mfn7+7O33nqLbzNo0CD22muvGdmXmZnJALDPPvusznufPHmSAWCHDh1i9+7dYyUlJYyxyoyXTz/9lDHGmF6vZ1FRUWzAgAEsPT2dpaamsl69erFBgwbVarMhDrdmzRoWHBxs8uvaWIqKipi/vz/729/+xs6fP8+2bt3K3Nzc2Ndff8232bVrFwsNDTV63urVq9nZs2fZ+fPn2fLly5mzs7PRHAfHcSwmJoYNHDiQZWRksLS0NBYbG8see+yxeu0x+K3X6y3qpz1jbZ/LKnTsdI59vaZFRUVWjTHbzeSfLWmKMDPG2J07d9isWbNYcHAwc3FxYW3btmVPPPEE+/333xljlR/cDz74gAUEBDCJRMJGjBjBfvjhB5OFOTExkS1ZsoR5e3szDw8PNmPGDKMvodqEmTHG1qxZU++9GWPs5ZdfZt7e3gwAW7p0KWPMWJgZY+zmzZs10uWUSmWtNjeXMDPG2JkzZ1j//v2ZWCxmbdu2ZR9++KHR9e+++45VH4MMGTKE/0KPjY1lv/zyS41+b9++zcaNG8c8PDyYv78/mzp1Knvw4EG9tpAwW5ayCh37T8oN9luW0q5eU2sLM+1gAtrBxBo4os+AY/ptLZ+rxpTtbUWftXcwaRFZGQRBOB6cnqGzrwdiO3rZlSg3By0iK4MgCMdBo+Vw/Mo9uDqL0LeTt8OJMkDCTBCEHWEIX2i0eggdT495SJgJgrAL7Dmm3NxQjJkgCLvASShAj0AZIts5zsRpXdCImSAIm6LRcjhwPg86PUNUkNzhRRkgYSYIwoYYwhcioRBiJ5IjA/RKEARhEyp0eoop1wHFmAmCsAnOIgH6dPBCFz8PEuVq0IiZIIhmpXI7qFyoNTp09ZeSKNcCCTNBEM2GIaYsdXWGpyv9YK8LEmaCIJoFvZ5h9+nbFFNuBPSVRRCE1WGMQSgUYFCILwJkrWg3aytBI2aCIKyKRsth25+3kK/WIFAuIVFuBCTMBEFYDUNM2dtDDD9pzX02idohYSYIwiowxrD/bB7FlM2AYswEQVgcHaeHk0iI+G5+kEmcSZRNhEbMBEFYFI2Ww7a0W7h2rxhyNxcSZTOwqTAfO3YMCQkJCAwMhEAgwJ49e4yuT506FQKBwOgxcuRIozYFBQWYNGkSPD09IZfLMX36dBQXFzejFwRBGKhaurOTj3vDTyBqxabCXFJSgsjISKxbt67ONiNHjkReXh7/2LJli9H1SZMmISsrCwcPHsT+/ftx7NgxvPjii9Y2nSCIWjh88S7FlC2ATWPMo0aNwqhRo+ptIxaLoVAoar128eJFHDhwAH/++Sd69+4NAFi7di1Gjx6NTz75BIGBgRa3mSCImpTr9OD0DEPD/ODqLCRRbiJ2P/mXlJQEPz8/tGnTBkOHDsWKFSvg7e0NAEhOToZcLudFGQDi4+MhFAqRmpqKJ598stY+y8vLUV5ezh+r1WoAlbPI1tg03NCvI21I7og+A47pd1mFDvvO5eORLgJEtJMDQKv339r+2bUwjxw5EuPGjUPHjh1x7do1vP322xg1ahSSk5MhEomgVCrh5+dn9BwnJyd4eXlBqVTW2e+qVauwbNmyGudVKpXVhNkQ93aUkYQj+gw4nt8aLYf95+/CXcQhyJ1BpVLZ2qRmwdp+2rUwT5gwgf9/eHg4IiIi0LlzZyQlJWHYsGFm97to0SLMnz+fP1ar1QgKCoJMJoOnp2eTbK4Ng9jLZI6zZY4j+gw4nt+Z2ffQ3k+O3gFiyOW0+4ilsGthrk6nTp3g4+ODq1evYtiwYVAoFLh7965RG51Oh4KCgjrj0kBl3FosrrkKyZD5YQ2qZpa0FjiOg0gkqvN6a/S5MTiC3xotB8aA/l19IBRUDm5au89VsbafLSqPOTc3Fw8ePEBAQAAAIC4uDkVFRUhPT+fbHDlyBHq9HrGxsbYy0yHYs2cP3n///RopjkTrx5ASdza3CE4imuizBjYV5uLiYmRmZiIzMxMAcP36dWRmZiInJwfFxcV44403kJKSghs3buDw4cNITExEly5dMGLECABAt27dMHLkSMyYMQN//PEHTp48idmzZ2PChAmUkWFFOI7j37PMzExwHNcs9yRsT9U85T4dvWxtTqvFpsKclpaG6OhoREdHAwDmz5+P6OhoLFmyBCKRCGfPnsUTTzyBkJAQTJ8+Hb169cLx48eNwhCbNm1CWFgYhg0bhtGjR6N///745ptvbOWSQyASiRAVFQUAiIqKqjecYQlodG4/ZN1RUZ5yMyBgrT2vpRGo1WrIZDKoVCqrTf6pVKpWNyFUX4zZUj5zHIf333+fP168eLHVvwiaQmt9rzVaDqUVHNq4OQMwjrG2Vp/rQ6VSQS6XW00zWlSMmbAvmkMgm3t0TtTEEL7IuqNyqAk+W9KisjIIx2Ts2LFISEggUbYBVWPK/bv42Noch4GEmWgRkCjbhtzCUvhLXTGMYsrNCgkzQRA10Gg53HtYji5+UnTxk9raHIeDYswEQRhhCF9cufvQ1qY4LCTMBEHwVI0pDwn1a/gJhFWgUAZBEDxFpVoEyiUYHOJLMWUbQsJMEAQ0Wg7X7hWjR6AMCpmrrc1xeMwS5uvXr+P48eO4efMmSktL4evri+joaMTFxcHVld5Uwro0VDyJMI2q4YvuAZ40UrYDTBLmTZs24fPPP0daWhr8/f0RGBgIiUSCgoICXLt2Da6urpg0aRIWLlyI4OBga9lMODB79uxBZmYmoqKiMHbsWFub0+KpKsq0zNp+aPTkX3R0NL744gtMnToVN2/eRF5eHtLT03HixAlcuHABarUae/fuhV6vR+/evbFjxw5r2k3YAFsXErJF8aTWTrlOjw7e7iTKdkajR8wffvghX9WtNsRiMQYPHozBgwfjgw8+wI0bNyxhH2En2MNI1bA822AHhTPMR6PlcDqnCLEdvdCPVvTZHY0W5vpEuTre3t78vnxEy6f6SNWWy6NpeXbTqRq+oEGyfdLkrIyff/4ZSUlJ4DgO/fr1w/jx4y1hF2FH2NtI1db3b8lQTLll0CRhXrx4MXbt2oUxY8aAMYZ58+YhKSkJa9eutZR9hJ1AI9XWgVAgQJhCipj2bUiU7RiThDktLQ29e/fmj7dt24YzZ85AIpEAAKZOnYrBgweTMLdSDKJM6WotD42Ww4kr9zEwxBe9gmnnEXvHpCXZL7/8MubOnYvS0lIAlZuj/uMf/0B2djbOnTuHr776CiEhIVYxlLAPaDeRlochfMEAOItolNwSMEmYU1NTERAQgJiYGOzbtw8bNmzA6dOn8eijj2LAgAHIzc3F5s2brWUrYWMoXa3lUa6jmHJLxKRQhkgkwsKFC/H0009j5syZcHd3x5dffkkbnzoI9jYJSDSMi0iImPZtEKaQkii3IMyqLtepUyf8+uuvePLJJzFw4ECsW7fO0nYRdsrYsWOxePFiWnVn52i0HPZm3kZJBYdutMy6xWGSMBcVFeHNN99EQkIC3n33XTz55JNITU3Fn3/+ib59++LcuXMm3fzYsWNISEhAYGAgBAKBUdxSq9Vi4cKFCA8Ph7u7OwIDAzF58mTcuXPHqI8OHTrw+5AZHh9++KFJdhCm4Ygj5ZYUtjHElN1cnODu4njvVWvAJGGeMmUKUlNTMWbMGGRnZ2PmzJnw9vbGxo0b8cEHH+DZZ5/FwoULG91fSUkJIiMjax1xl5aWIiMjA4sXL0ZGRgZ27dqF7OxsPPHEEzXaLl++HHl5efzj1VdfNcUtgqiXljThyekZdmXcpphyC8ekGPORI0dw+vRpdOnSBTNmzECXLl34a8OGDUNGRgaWL1/e6P5GjRqFUaNG1XpNJpPh4MGDRue+/PJL9OnTBzk5OWjfvj1/XiqVQqFQmOIKQTQKe1r12BCMMYiEAvTv4oMgLwmJcgvGJGHu2rUrvvnmG/z973/HwYMHa1SQc3V1xcqVKy1qYFVUqsrt0+VyudH5Dz/8EO+//z7at2+P5557DvPmzYOTU92ulZeXo7y8nD9Wq9UAKj/YjDGL223o1xp924LG5DG3Fp+FQqHRhKdQKKzXJ1v5rdFy2HfmDoaE+SHIS8Lb0hy0lvfaFKztq0nCvGHDBvz973/HunXrEBUVhX/961/WsqsGGo0GCxcuxMSJE+Hp6cmfnzNnDmJiYuDl5YVTp05h0aJFyMvLw5o1a+rsa9WqVVi2bFmN8yqVymrCXFxcDAAtfhSTlJSE7OxshIaGYvDgwXW2a00+Dx48GAMGDIBIJIJKpaq3rS381mg57M+6Bx93ZzjpyqBSaZrlvgZa03vdWBr6HDQVAbOTrzmBQIDdu3fXOtuv1Woxfvx45ObmIikpyUiYq7Nhwwa89NJLKC4uhlgsrrVNbSPmoKAgFBUV1du3uTDGoFKpIJPJWvQHl+M4rFixgj9+99136xw5txafTaW5/WaMYWd6LuRuLjaLKTvie61SqdCmTRuoVCqraEajR8yMMZu86FqtFs888wxu3ryJI0eONPgixMbGQqfT4caNGwgNDa21jVgsrlW0DVkd1qBq1khLxcnJyehnfX3hIqB1+GwOzeU3p2cQCYUY2s0f3u4uNn2dHe29trafjc7K6NGjB7Zu3YqKiop62125cgUzZ860SMqaQZSvXLmCQ4cONaqUaGZmJoRCIfz8aIdfa0B5zPaBRsth25+3kPOgFD4eYocRREeh0SPmtWvXYuHChXjllVfw2GOPoXfv3ggMDISrqysKCwtx4cIFnDhxAllZWZg9ezZmzpzZYJ/FxcW4evUqf3z9+nVkZmbCy8sLAQEBeOqpp5CRkYH9+/eD4zgolUoAgJeXF1xcXJCcnIzU1FQMGTIEUqkUycnJmDdvHp5//nm0adPGjJeDaAz2mpXgKFQt3WmY6CNaGcxEjh8/zmbPns0iIyOZXC5nYrGYtW3blj3++ONs7dq1rKCgoNF9/f777wxAjceUKVPY9evXa70GgP3++++MMcbS09NZbGwsk8lkzNXVlXXr1o2tXLmSaTQak3xSqVQMAFOpVCY9r7Ho9XpWWFjI9Hp9k/rR6XQWssj6/VrK55ZGc/i953Qu+y1LaTevrSO+10VFRVbVDLuZ/LMlarUaMpnMaoF8ZoHJEWtt7WStfi3hc0vEmn6X6zi4iIQoqeDg7iKym9fVEd9rlUoFuVxuNc0wq1YG0bxYq6obVYtrOWi0HHam5+KS8iE8xE4OI4COCglzC8BQ1Q2ARau6WatfwrJUjSmHKaS2NodoBiiUgZYRygCst3OINfp1xJ+3gHX8PnQhHwyw29oXjvheWzuU0eTNWInmw1ojWhop2ycaLQehQICBIb5wFjlOjjBBoQyCsEsM4Ytzt4vg4iQkUXYwTBoxG4r9NIQ1hvYE4ShUjSnHtKd8fEfEJGGWy+X1fnOz/79sm2b3CcJ8TucUUT1lB8ckYf7999+tZQdB2ITGTHxaa9K1Ohoth3KdHrEdvaDXcyTKDoxJwjxo0CBr2UEQzU5jFtdYawFOdQzhiw7e7rh3/kSz3JOwX0ya/NPpdEblMgEgPz8fy5Ytw5tvvokTJ05Y1DiCsBaNWVzT1AU4jW1fNaYc20FOi34I04R5xowZmDNnDn/88OFDPPLII1i3bh1+/fVXDBkyBL/88ovFjSQIS9OYxTVNWYCzd+/eRu8TeO1eMR9TNpRWNeeeROvBpAUmISEh+PLLLzF8+HAAwLp167By5UpcuHABMpkMCxcuxB9//NHiYtEtZYGJvVJbDNYUn60dw62vf0vHmBljKCgowNq1a/lzixcvrvX5Gi2HolItFDLXGvXOmyuubQla++e7NuyqVsbt27fRtWtX/vjw4cMYP348ZDIZgMpdtLOysixrIWHXNHUHaWvvQF1X/4YQQXXxqy10YKpANmakrdFy2JVxGxeVlSmo1QWtpYgyYR1MEmZXV1eUlZXxxykpKYiNjTW6btj7i7Au9hB7tEQM1prx1Lr6r0usLfklkZiYWOeGAgZR9pWKMTjEt8n3IlofJglzVFQU/v3vfwMAjh8/jvz8fAwdOpS/fu3aNQQGBlrWQqIG1h5lNpamFkGydhGl2vqvS6zrO9+U+9fGvYflUMjElKdM1IlJ6XJLlizBqFGjsH37duTl5WHq1KkICAjgr+/evRv9+vWzuJHE/6guIAkJCTb92Tt27Ngm2dDQ85saa63ev0GsDelo9Z23dKqcRssht7AUXfykCPJya3J/ROvF5DzmtLQ0HDx4EAqFAk8//bTR9aioKPTp08eiBhLG1CUstrbJGs+3lDBW77+uL4Oq5y39BWhIifOXuqKLH5XuJOrH5Opy3bt3R+fOnaHT6SAUGkdCXnzxRYsZRtRNU0epLQFr/zKoq6+GRtbmUDVPeVg32iSYaBiTYsz37t3DqFGj4OHhAU9PT/Tt29doM1Wi+bAXUbbWJKQ9FPFvyo7gVV+X0goO7b3cKKZMNBqThHnhwoXIzMzE8uXL8cknn6CoqAgzZsww++bHjh1DQkICAgMDIRAIakxmMcawZMkSBAQEQCKRID4+HleuXDFqU1BQgEmTJsHT0xNyuRzTp0+nzJBmwtqTkE0RRktgbnw7KSkJK1aswPZde5B+sxBt3JwxoKsviTLRaEwS5oMHD2Ljxo1YtGgR5s2bh3379uH48eM1lmk3lpKSEkRGRmLdunW1Xl+9ejW++OILrF+/HqmpqXB3d8eIESOg0Wj4NpMmTUJWVhYOHjyI/fv349ixYxRSaQaaa79AW/0yMPdLh+M4ZGdnQ8cE2J1xG/cfahp+EkFUwyRhvnPnDiIjI/njrl27QiwWIy8vz6ybjxo1CitWrMCTTz5Z4xpjDJ999hneffddJCYmIiIiAj/88APu3LnD/7FcvHgRBw4cwL/+9S/Exsaif//+WLt2LbZu3Yo7d+6YZRPROEwJNdhDzrUpNOVLRyQSoWPXUFzjfNCjUzsM76GgkTJhMiZP/lX/AxSJRLDGtoHXr1+HUqlEfHw8f04mkyE2NhbJycmYMGECkpOTIZfL0bt3b75NfHw8hEIhUlNTaxV8wnI0ZhJy7969Fks5a65lyk2d+BswYADahXOI60LhC8I8TBJmxhhCQkKMPmzFxcWIjo42ytAoKChosmFKpRIA4O/vb3Te39+fv6ZUKuHnZzzL7eTkBC8vL75NbZSXlxuFXww7szDGrPIlY+i3Ne57KxQKa/WLMQadTmc08nz88cfNFtaqAp+YmNgUkxtFYmIib29j3zeNlkPKXw/Q09cZfTv7AECrfM+r05o/33VhbV9NEubvvvvOWnY0K6tWrcKyZctqnFepVFYTZsOEpKOMoBhjKCsrQ1RUFLKzsxEaGmr2pKwhbiuRSJCdnY2CggK7yUoxoNFy2H/+HrzdnVHm5gyV0HE2T3XEz7dKpbJq/yYJ85QpU6xlRw0UCgWAynrPVVcX5ufn87FNhUKBu3fvGj1Pp9OhoKCAf35tLFq0CPPnz+eP1Wo1goKCIJPJrFZdDoBDVd8y+PzEE09Ar9c3WUhDQ0ORkZGBmJgYeHl5WcJEi1Gu5fBzxm2095NjWJgvX62wtve6JVWNayyO+Pm2NibHmJuLjh07QqFQ4PDhw7wQq9VqpKamYubMmQCAuLg4FBUVIT09Hb169QIAHDlyBHq93qi4UnXEYjHEYnGN8wKB9UY5hr4d6YNr8NfJybIfM3t7DZ2dRIhoJ0fPtpVf6nW91821G4otcLTPt7X9NCkrw9IUFxcjMzOTj0Nev34dmZmZyMnJgUAgwNy5c7FixQr89NNPOHfuHCZPnozAwED+Q92tWzeMHDkSM2bMwB9//IGTJ09i9uzZmDBhAhVTakUYsiSEQiEyMzNRUVFha5MAVIYvfj6bhwqdHuHt6h8tNld6IdE6sKkwp6WlITo6GtHR0QCA+fPnIzo6GkuWLAEAvPnmm3j11Vfx4osv4pFHHkFxcTEOHDgAV1dXvo9NmzYhLCwMw4YNw+jRo9G/f3988803NvGnJdGShKFqal55eTlWrlxp88p6hmXWLk5CuDo3/GdkDysZiZaDSTuYtFYcbQeT5vhJbQ2fKyoqsHLlSv64rp1BrI2W02N72i1+O6iq/jXkd2uNMdvT57s5sKsdTIiWT0v+Se3i4oKIiAgAtht1MsbgLBKibydvs2pftDZRJqxDo2dlqmYxNMSaNWvMMoawPvZYNrSx7NmzB2fPnkV4eLhNJs80Wg77z+YhvpsfOvt6NPv9Cceh0cJ8+vRpo+OMjAzodDqEhoYCAC5fvgyRSMRnRxD2S0ssG1p1pH/u3DmMHTu2We2vWrpTJnFutvsSjkmjhbnqztdr1qyBVCrF999/jzZt2gAACgsLMW3aNAwYMMDyVhIWpyWJMmDbkT5jDHtO3641pkwQ1sCsyb+2bdvit99+Q48ePYzOnz9/HsOHD29xBYQcbfLPFOraTbohrOVzc0+e6fUMQqEA+WoN/KTiBn1pye+1uTiiz3Y5+adWq3Hv3r0a5+/du4eHDx822SjCPtizZw8mTZqESZMm2Tw9zUBzhy+2pd3CnaIy+Hu6Wl10WtJELGFdzBLmJ598EtOmTcOuXbuQm5uL3Nxc/Pjjj5g+fTrGjRtnaRsJG8BxHDIyMqBUKqFUKpGRkeFQwmGIKft4iBEgc234CU3EXnY+J+wDs4R5/fr1GDVqFJ577jkEBwcjODgYzz33HEaOHIl//vOflraRsAEikQgxMTFQKBRQKBSIiYlpcXHppvDLubxmiym35BRGwjo0aYFJSUkJrl27BgDo3Lkz3N3dLWZYc0Ix5rqxtxiztanQ6eEsEkCt0cHT1clk2831uyXX0Wip73VTsHaMuUnVZfLy8pCXl4eBAwdCIpGAMeYwb4yj0NpGyfVNHhrCF306eKGrv7RZ7WqJKYyE9TArlPHgwQMMGzYMISEhGD16NL+11PTp0/H6669b1EDCfmjpP7Hri+NWzVPu4mebxSMkyoQBs4R53rx5cHZ2Rk5ODtzc3Pjzzz77LA4cOGAx4wjLY664tvTJqYbiuEnZdylPmbAbzBLm3377DR999BHatWtndL5r1664efOmRQwjLE9Tdn5u6ZNTdVV302g56Dg9Bof6kSgTdoNZwlxSUmI0UjZQUFBQawF6wvY0defn1lCycuzYsVi8eDE/uabRctiVcRvn76jh6iwiUSbsBrOEecCAAfjhhx/4Y4FAAL1ej9WrV2PIkCEWM46wHE0V1+qi1lKpOlLelXEbvlIxItvJbGwVQRhjVlbG6tWrMWzYMKSlpaGiogJvvvkmsrKyUFBQgJMnT1raRsJCNHXmv6WOlGvjzxsF8JWKKXxB2CVmCXPPnj1x+fJlfPnll5BKpSguLsa4ceMwa9Yso41TCfuiNRZpNxWNlgOnZ3i0sw+EAvvbP5AgADOEWavVYuTIkVi/fj3eeecda9hEWIGWvIDBUhhS4jr7eqBvJ29bm0MQdWJyjNnZ2Rlnz561hi2EBahtUq8lZlVY2saqecqxHb0s2jdBWBqzJv+ef/55fPvtt5a2hWgidaXDtbSsCmvkTF9SPqQ8ZaLFYFaMWafTYcOGDTh06BB69epVo0aGJbeW6tChQ6250a+88grWrVuHwYMH4+jRo0bXXnrpJaxfv95iNrQEqo+Kq0/ytZQlvw35YSoaLYfich2feUGiTLQEzBLm8+fPIyYmBkDlllJVsfQH/88//zT6WXv+/Hk89thjePrpp/lzM2bMwPLly/nj2nKsWzuN2eHD3kUZsOxOJYbwRVAbNwwM8TWrD5owJWyBWcJcdZspa+Pra/wH9eGHH6Jz584YNGgQf87NzQ0KhaLZbLJXWsqouCEs4UfVmPKArj5m9UETpoStaFJ1ueamoqIC//nPfzB//nyjkfmmTZvwn//8BwqFAgkJCVi8eHG9o+by8nKUl5fzx2q1GkBl+cImVEGtE0O/1ui7OkKhsFnu0xBN9bmpftwuLIW/VIyhYb68PaZQPaTy+OOPN+qLojnfa3vBUX22JmYJ85AhQ+oNWRw5csRsg+pjz549KCoqwtSpU/lzhmL9gYGBOHv2LBYuXIjs7Gzs2rWrzn5WrVqFZcuW1TivUqmsJszFxcUAHCfGaSufNVoO94orENRGAu9AV/5L1xyioqKQnZ2N0NBQ3peGoPfaMXxWqVRW7d+sQvnz5s0zOtZqtcjMzMT58+cxZcoUfP755xYzsCojRoyAi4sL9u3bV2ebI0eOYNiwYbh69So6d+5ca5vaRsxBQUEoKiqiQvkWwhY+G5ZZ+3uKMaybv0X6NDXGTO+1Y/isUqnQpk0b+yqU/+mnn9Z6/r333mv0yMJUbt68iUOHDtU7EgaA2NhYAKhXmMVica3FlgQCgdU+WIa+HeWDCzSvzxoth12nb8PP0xXDLJgS5+Rk+p8IvdetH2v7aVYec108//zz2LBhgyW75Pnuu+/g5+eHMWPG1NvOEBekpeGOhbpMi7ZyCeUpE60Ci07+JScnw9XV8jsK6/V6fPfdd5gyZYrRCObatWvYvHkzRo8eDW9vb5w9exbz5s3DwIEDERERYXE7CPtDo+VwJb8Y4e1k8PO0/m7WBNEcmCXM48aNMzpmjCEvLw9paWlYvHixRQyryqFDh5CTk4MXXnjB6LyLiwsOHTqEzz77DCUlJQgKCsL48ePx7rvvWtwGonlpTGy3akpcT+ZJI2Wi1WCWMMtkxvVrhUIhQkNDsXz5cgwfPtwihlVl+PDhtWZLBAUF1Vj1RzSOhoTPlgsrGpM/XFWUDeELU2ymhSOEPWOWMH/33XeWtsMhsZU4NCR8tlxY0dgl2VpOj44+7ojr5A2BQGCSzbRwhLB3zJr8u3XrFnJzc/njP/74A3PnzsU333xjMcNaO5Yu1FNXVbnaztVXac7WlegaKrik0XI4efU+xCIBHu3sw4+UG2uzrf0jiMZgljA/99xz/LJspVKJ+Ph4/PHHH3jnnXeMalYQtWNpcahN5M2tNGcPlejq2sbKEL74/dhJfLBiBe+bKTbbg38E0RBmLTBp06YNUlJSEBoaii+++ALbtm3DyZMn8dtvv+Hll1/GX3/9ZQ1brYZarYZMJrNasnhtCfiW+jnNcRzef/99/tgw+Vr9XHUBsnaM2ZKLDjiOg1YP/JiRCx93F5za8TUMXVb1zR5izI642MIRfVapVJDL5fa1wESr1fILNA4dOoQnnngCABAWFoa8vDzLWdeKaUqhnqqiUlc1tqZWmrOXkaThCyw8Igrdeg9CdJAcZVdr980Um+3FP4KoDbNGzLGxsRgyZAjGjBmD4cOHIyUlBZGRkUhJScFTTz1lFH9uCdhixGwudY20axsB2jLzwBI+cxyHpctX4LZehnZCFd5b8q5Zo+PmxBFHj47os7VHzGbFmD/66CN8/fXXGDx4MCZOnIjIyEgAwE8//YQ+ffpY1EDif9QXm26p9ZfrQ6sHygOjIQAQExVp9uiYIFoaZoUyBg8ejPv370OtVqNNmzb8+RdffNEhi9Q3F5YsIm/vVOj0+DEjF8MHxGJIyGizalYQREvF7E+7SCQyEmWgchsowrq0lmL4DeEsEqB3sBdC/D0c5ucxQRgwW5h37tyJ7du3IycnBxUVFUbXMjIymmwYUTetWZQ1Wg4HzisxtJsfQhVSW5tDEDbBrBjzF198gWnTpsHf3x+nT59Gnz594O3tjb/++gujRo2ytI2Eg2DIU3YXO0Eqrn3MQAtCCEfALGH+5z//iW+++QZr166Fi4sL3nzzTRw8eBBz5syxemV/onWi1zPsyrhtVPuiOpZeLUkQ9opZwpyTk4NHH30UACCRSPDw4UMAwN/+9jds2bLFctYRDgFjDEKhAANDfOoUZVpKTTgSZgmzQqFAQUEBAKB9+/ZISUkBAFy/ft2hNmQkmo5Gy2F72i3cfahBuzZudU700VJqwpEwa/Jv6NCh+OmnnxAdHY1p06Zh3rx52LlzJ9LS0mrUaiaIujDs0ecrFcPXo+ZWX9VxlIwUgjBLmL/55hvo9XoAwKxZs+Dt7Y1Tp07hiSeewEsvvWRRAwnTsNcVcdVhjGHfmTvwlYpN2g6qJfhGEE3FLGEWCoUQCv8XBZkwYQImTJhgMaMI82gpdYZ1nB5OIiGGdfNHGzdnylMmiGqYvRnr8ePH8fzzzyMuLg63b98GAPz73//GiRMnLGYc0XhayuRYZUw5Fzful8DL3YVEmSBqwSxh/vHHHzFixAhIJBKcPn0a5eXlACoLe6xcudKiBhKNoyVMjhnylH2lYgR709J9gqgLs4R5xYoVWL9+Pf7v//4Pzs7O/Pl+/frRqj8bUleBeXvh4IX8evOUCYKoxCxhzs7OxsCBA2ucl8lkKCoqaqpNPO+99x4EAoHRIywsjL+u0Wj4yUcPDw+MHz8e+fn5Frt/S8QeR8rlOg56PcOwbn4kygTRCMzOY7569WqN8ydOnECnTp2abFRVevTogby8PP5RNYY9b9487Nu3Dzt27MDRo0dx584dStezMzRaDjvTc3EhTw03FycSZYJoBGZlZcyYMQOvvfYaNmzYAIFAgDt37iA5ORkLFizgtzaymIFOTlAoFDXOq1QqfPvtt9i8eTOGDh0KoHL37m7duiElJQV9+/a1qB2E6Wi0HH7OuA0/T1f0CLR8MXGCaK2YJcxvvfUW9Ho9hg0bhtLSUgwcOBBisRgLFizAq6++alEDr1y5gsDAQLi6uiIuLg6rVq1C+/btkZ6eDq1Wi/j4eL5tWFgY2rdvj+Tk5HqFuby8nJ+wBCp3MAEqc2utsXLR0K8jrYpkjOGPm0XwlbphWJgvf66146jvtSP6bE3MEmaBQIB33nkHb7zxBq5evYri4mJ0794dHh4eFjUuNjYWGzduRGhoKPLy8rBs2TIMGDAA58+fh1KphIuLC+RyudFz/P39oVQq6+131apVWLZsWY3zKpXKasJcXFwMAA7xU16jrUzVC/dxhsxTzH/xOQKO9l4DjumztYu1NWlbCBcXF3Tv3t1SttSgagnRiIgIxMbGIjg4GNu3b4dEIjG730WLFmH+/Pn8sVqtRlBQEGQymdX2/APgEHuiGcIXXfzcEdJG6hA+V8WR3msDjuiztTFJmF944YVGtduwYYNZxjSEXC5HSEgIrl69isceewwVFRUoKioyGjXn5+fXGpOuilgs5nf5rooh88MaVM0saY1wHAetHth1ujKm/EgHL6jV6lbtc1209ve6NhzNZ2v7aZIwb9y4EcHBwYiOjrZJPKm4uBjXrl3D3/72N/Tq1QvOzs44fPgwxo8fD6AyjS8nJwdxcXHNbpsjY1gKLu0UjYjefRHfzc/WJhFEi8YkYZ45cya2bNmC69evY9q0aXj++efh5eVlLduwYMECJCQkIDg4GHfu3MHSpUshEokwceJEyGQyTJ8+HfPnz4eXlxc8PT3x6quvIi4ujjIymhGO45B2+gw4iMCuncaQ58ZAIBA41EQQQVgak/KY161bh7y8PLz55pvYt28fgoKC8Mwzz+DXX3+1yh9ibm4uJk6ciNDQUDzzzDPw9vZGSkoKfH0rZ/k//fRTPP744xg/fjwGDhwIhUKBXbt2WdyO1oSla2ho9UB5YDQe6N0RHR1Fu1kThAUQsCYo6s2bN7Fx40b88MMP0Ol0yMrKsnhmRnOgVqshk8mgUqmsNvmnUqlsPjli6epzhtoXflJXDAnxNhJle/G5uXFEvx3RZ5VKBblcbjXNMLu6HFBZ/tPws9Veq5kRlVij+tyNByV87QsaKTcv9PfWujFZmMvLy7FlyxY89thjCAkJwblz5/Dll18iJyenRY6WHQVLVp/TaDncLipDmMKzwdoXJCCWhzalbf2YNMx55ZVXsHXrVgQFBeGFF17Ali1b4OPjYy3biAYwdbcSS2zNZAhfBMokaCuX1CvKSUlJLaJwf0ui+i8f2mqrdWKSMK9fvx7t27dHp06dcPToURw9erTWdjQBZ33MjRdbQpT9pK4YHOpbb1uO45CdnQ2ABMSSGH75GN57ek1bJyYJ8+TJkx0muG/P1DZqAqxf8vNBSQUCZRIMDvVt8HMgEokQGhpKAmIFaFPa1k+TsjJaCy0xK6PqiBmAVUMGGi2HGw9KEKZo/Gtj8NnDw8OhJgYdMUPBEX2266wMwnYYditJSEiw6l5/hvBFbkGZWbnqNKojCNMhYW7BiEQiq+71VzWmPIx2HiGIZsNxfmO2YqwVc9RoOQR7uaNfF28SZYJoRmjE3Eqw9Ej5j+sFkEmc0b+rD4kyQTQzJMwtHGvFlFVlWov2SxBE46FQRgvGmrUvrL2btamLYwjCkaARczNT2wjXnFGvNWpfAECIv9TqokxLigmifkiYm5G9e/fWECRzRcqS2RglmgocuZQPkVCARzp4WX2kbM30PoJoDZAwNxO1CVJTRcqQy2wIY5gjctt37cHfV/wfTp1KgZPQ+pN81kzvI4jWAsWYm4m6ahw0te6B4TnmxJtLy7XYnXEbEoEW2punodc/XqcNlowJ05JigqgfEuZmJDExsYYgmSNS1UXS3Ipjrs5O6NetLe5dzkR0dN1fDJaeZARoRSBB1AcJczNTmyDVdq6uEWptImlqxTGNlsPBC/kYGuaH2ZPGguPqFnIqM0kQzQ8Jsx1hEOO6RqjVRXL06NH8suzGjryrpsS5uVS2re85VGaSIJofu578W7VqFR555BFIpVL4+flh7NixfI1fA4MHD4ZAIDB6vPzyyzay2HwM2Rm7du3ixTcjI8NoQq/qxFl5eTkmT56MSZMm8RkdDYmmjtNjV8Ztk/OUq08ytmYoS4SwB+xamI8ePYpZs2YhJSUFBw8ehFarxfDhw1FSUmLUbsaMGcjLy+Mfq1evtpHF5lF1JHz27FmEh4fj0qVLuHz5Mvbt22fUNiEhAW+//TacnZ2hVCqhVCqRlpbWoKAwxuAkEuLRzt5m5Sk7wkiZ8qsJe8GuhfnAgQOYOnUqevTogcjISGzcuBE5OTlIT083aufm5gaFQsE/rFEf1ZpUTyEbO3YsQkJCEBYWZpRGt2fPHixbtgy//PILIiMjoVAooNPpcPXq1RoCXhWNlsPO9FwUlFSgg497o0XZkUaPlF9N2BN2LczVUalUAAAvLy+j85s2bYKPjw969uyJRYsWobS01BbmNYnq4YKYmBgA4AWb4zhs3boVx44dwyeffIIzZ85g/PjxGDp0KC/gFRUVNfot0VTgx4xcyN1c0MbNudH2ONrokfKrCXuixUz+6fV6zJ07F/369UPPnj3588899xyCg4MRGBiIs2fPYuHChcjOzq5338Hy8nKUl5fzx2q1GkDlz31rbOhi6Ldq37VlXQiFQqOJv3fffZdfLRgREcGPdAsLCwEA586dQ1RUFM6ePYvy8nKsXLkSUVFRSExMBADs2bMXO9Jz0b1TECZOTeBtaYjqo8fHH687v9kUn+2dxMRE3ldz7W6JfjcVR/XZmrQYYZ41axbOnz+PEydOGJ1/8cUX+f+Hh4cjICAAw4YNw7Vr19C5c+da+1q1ahWWLVtW47xKpbKaMBcXFwMABAIBkpKSkJ2dja5du2Lo0KEA/hc2yM7OhkQiQXZ2NnQ6HbZv3w5vb29cvnwZ48aNw5YtW6DVanHq1Cl07doVjDHMmjULGzZsAGMMly5dwoABA6BnDJcvZ6OruxOc7l5CYWF/k8Q1KioK2dnZCA0N5W1vis+OgiP67Yg+G369W4sWIcyzZ8/G/v37cezYMbRr167etrGxsQCAq1ev1inMixYtwvz58/ljtVqNoKAgyGQyq+35BwAymQx6vR6ZmZnIzs7GgQMH+GuG0a9hA9Pw8HCcO3cOQqEQJ06cwLlz5xAYGIjLly8DALKyslBWVoaLFy9CIBCgqKgI+/fvh1wuR6eQMFQERMK3QyhyLmYiLCqqRvinIRITE5u02q+qz47yxwo4pt+O6LO1sWthZozh1Vdfxe7du5GUlISOHTs2+BzDT/CAgIA624jFYojF4hrnDel21sDQt5OTEyIiIpCUlAQ/Pz9s27YN+fn5vL2LFy/G6NGj4eLigrS0NOTl5UEqlUIoFOLu3buQSqW4c+cOFAoFrly5Al9fX2zduhVdunSBl5cXmNAZ//w5DR0D/sITvTphyuLF/P1NxZxNVKuKedUURkfCEf12NJ+t7addC/OsWbOwefNm7N27F1KpFEqlEkDlN7NEIsG1a9ewefNmjB49Gt7e3jh79izmzZuHgQMHIiIiwsbW1824ceOQmpqKlJQUcByHdu3aQalU4umnn8a+ffuQkZGBqKgoiMViDBgwAD/++CPu378PsViM/v37o3///pBIJDh8+DBEIhFu374NgUAAjuPg0rkvnPQatBepceZMJvR6jh+NWzsPuWp83BDnJgjCdOw6K+Orr76CSqXC4MGDERAQwD+2bdsGAHBxccGhQ4cwfPhwhIWF4fXXX8f48ePrTR2zBziOg0QiwaBBg9CuXTuEhIRgzpw5GDt2LJ95sXXrVpSVleHo0aMoKSlBt27dIBKJ0KVLF7i7u2PRokXw9fXFhQsXkJ2dDSYQYciQodjx2RLMfKI/BAKgR48eOHfuHADrp4BRuhlBWA67HjE3NBEXFBSEo0ePNpM1lqPqMucJEybwS6kNYnb//n3cu3cPt2/fhpOTE9zdK3OPe/fuDScnJ5SVlWHYsGFIT0+vDB24SHC6SIxBodGQu7viySfH8iPl8vJyiMXiJqWANSbWTEu3CcJyCJgj5bjUgVqthkwmg0qlstrkn0qlgoeHh1HctrYqcXv37sVnn30Gf39/XLx4EV5eXggICMDGjRshkUiwfft2rF27FpcuXUJhYSH0AhE8Ix6Dn6crxvbpDDeJBBERETh79izf79tvvw0XF5c671sfplaWM/Rt8NnRJoQc0W9H9FmlUkEul1tNM+w6lNGaSEpKwooVK7Bnzx5+ZFxVHA0LOv78809wHIcLFy7Aw8MD9+/fx+3bt/Hrr7+C4zjs2rWrUpD1eri5uUHe7VFAo4Z3+R2kpqRAr9fj7Nmz6N69O4CaiyVMWThiTnjC1PKlBEHUhIS5GeA4ji++tHXrVixbtsxIGA0CqNfrcerUKRQVFeHu3bu4cOEC8vPzUVBQwBc00uv18PHxwZD44fjn+q/R00uA7nIOapUKOp0Oly9fRllZGS5cuMBPgBqE2FShteZqOEdbWUgQpkDC3AyIRCKEhoZCr9cDqFzhV10Yo6KiIBQK0bdvXwCVqxNLSkpQUFCAy5cv4/z58/jll18AAD7+Aeg4dCK69x+JuXNmI+j/53YHBQWB4zikpKTg0qVLyMzMRFpaGoD/pREaxLqxQmuNynI0UUgQ9WPXk3+ticGDB+Pxxx/H/v37kZGRwdfCqBrDNcSCdTodNm7cyE9+ikQi3Lt3D5s2bUL+/UIIO/fFjH6PICJQisjERGRmZkIgECAvLw9+fn4ICAiAUqlEcHAwbt68CQCYMGEC9u3bx1evM0Vo6xNwcxah0EQhQdQPCXMzUlWAkpOTkZGRgezsbHTr1g1bt25FRkYGysvLsWPHDpSWlkIgEEAkEsHV1RV9+/bFrVu3UC4NhIzT4OGlU1i+/QwiIyMRExMDoVBoJMSzZs1CVlYWwsLCoNfrMXr0aKxcuRJA5SrDhIQEowlBc2jKllO07x9B1A2FMpoRw0/4S5cuYdu2bTh06BCysrJw/vx5vo0hxlxRUQEXFxf06dMHX3/9Nd5f+SHGPjUB8eHt8dLjj2Lnzh3YtWsXli5dCr1ej4ULF0IikSAsLAwhISEYN24cHx+OiYmBi4uLUZH9lStXNim+a4lwBIkyQdQOCXMzIhKJEBERgfz8fLi6uuLChQvgOA737t1DcHAwLl++zFfpcnJygrOzMwIDAyF0ccWPGbm4VSKASCSEQFBZba+4uBiMMWzduhUffPABysrKoNfrERMTw283VTU+PHbsWLz99tv8cvS6SoU21hcqk0kQ1oFCGc0Ix3EQCoXQarV48OAB/Pz8UFZWBoVCAWdnZ3Tt2hV6vR7Z2dlo27YtRCIRuoZ1x+6M25jcpTcu/3Wanzh89tlncefOHRQUFODcuXO4cOEC1Go1YmNj+fg1UHNUahg5Z2ZmGpUKNWdyLyEhgcIRBGEFaMTcTCQlJeG9997Dli1boFKpUF5ejsLCQj5L4vLly6ioqMDdu3chFov54kUqJkGPTu3wWHcFoqOj+C2nBAIBAgMD4ePjg+LiYjx8+BAajQb5+fnIyMiodyRc28jZ1FCEId3N3pe/E0RLhIS5GeA4Dr/++itOnjyJM2fO4O7duwAqV0ydOnUKv/32Gy5cuICbN2/ipZdegpOTE9q27wC1VgCZXg119iksW7YMWq3WaMspAPykX58+feDt7Y2CggJkZWU1GEOuGnM2NRRB6W4EYV0olNGM6PV6ODk5oUuXLjh//jyKiorAGINerwfHcQgLC8MzzzyDtMyz+P2vEgjLinDlyhX89ttvYIxBLpdj4sSJEIvFiImJQWRkJM6cOYOYmBiMHj0aK1asAGMMx48f5+s+Vw01VE9tMzczgtLdCMK6kDA3AyKRCAEBAbhw4QLc3Nxw+/ZtyOVyAEBJSQlEIhG0Wi2ysrLwxltvI2rsS4iXuuDwvz/Hju3bUVFRAcYYysrK8Ndff6Fr165ITU2FWCzmc5I5jkNMTAwyMzMRFxcHoVBoJJp1pbYZiieZKq6U7kYQ1oOEuRnYtWsXTpw4gaysLH5LGmdnZwQFBeHevXsoLi6GVqvFzZs38c/vNsP32Fn0bivBg/v30bZtW35UHRwcjLy8POTn5+PBgwfo0aMHkpKSkJaWBolEgqioKCxevLiG2FYPPVQV1KbkIpMoE4R1oBizleE4DmfOnOHT23Q6HXQ6HSoqKlBaWgo3Nzc4OztDLxBBK/GCtvAOVBdP4uBvvyEzMxNZWVl49tln8ddff+H333/n9/VTKpXIysqCv78/UlNT+dCFgaqiWVdqW12xYooZE4RtoRGzFag6WhWJRIiJicHJkyf5rXeEQiEYY/Dw8MDdu3eh0XJw7zEErLQIggoVtFotOI6Ds7MzXF1d+RV6e/fu5SvLKRQKdO/eHaGhoejUqVON0EV1ags91BYrbsoImiAIy0DCbGGqCxvHcSgvL0dRURHfRq/XQyQS4cqVKxC5SODRcyi0Dx+g/K80eHh48ItMNBoN9Ho9vv/+e/z++++4desWysrKUFpaCrVajS5dumDp0qWNjhPXdr2qYNcW8qjreY2h6si7tlKnBEHUDgmzBeE4DhkZGfwiEL1ej23btiE5ORlubm78JB4A6HQ6AAATiFBRqITmxmkAlUX7gcrNHg2x5Xv37uHBgwe8oLu7uyM4OBjOzs78vZtSaKjq6L7qCHrfvn1mj5737t2L7OxshIaGAqgsdwpUFlOikThB1A8JswXZt28fLl++DAB4+umnkZ6eDqVSCZVKVaNgkMDJBc6+HVCRd5kX5aoYBNywUEQoFEIsFsPb2xt6vR6XL1+Gv79/gzaZGpowjKCByjrOQM0Jw4YwjLwlEgkyMjKg1+v5jXQzMjIom4MgGoAm/yyEQYxCQkIQEhKCU6dOYe/evbh69WqNvQsFTi5w7zEUInd5o/oWCAR8DvTTTz+NTp06oVu3bsjOzsbSpUvx448/8jbUZhNg2kIQkUjUpFoYVZ8bExOD3r17w8/PDwqFgq/jQRBE3bSaEfO6devw8ccfQ6lUIjIyEmvXrkWfPn2a7f4ikQjl5eVITk5GaWkpsrOzodPpUFZWBqHwf99/BlHmSgpQdvWPBvsVCoV8gX3GGFJTU5GTkwPGGDp06ICrV6/ixIkTSElJgbu7u9HIuKkLQZqSq5yYmIiCggJ4eXlh7969CAsLowlFgmgkrWLEvG3bNsyfPx9Lly5FRkYGIiMjMWLECH7ps7WoPrllCDWcOXMGarUapaWl/EQej0AIbcHtRokyAD6TQyAQoLy8HPfv34enpydEIhH69u0LpVKJiooK7NixAxcuXOBHxgbbDHUxGiuI1UfVdYlybaPz6lSdUBQKhTh79iyl4hFEI2gVI+Y1a9ZgxowZmDZtGgBg/fr1+Pnnn7Fhwwa89dZbFr2XYSJtz549yMjIQHh4OIYPHw4PDw88fPgQhw8frrWAEBM6Q9whCmXXM1F+65xJ9zOg0+lw8+ZNMMbQqVMnuLm5ITo6Grt27YKbmxvy8/Px7LPPGk3aAagRYzb4UH1SsLHx6KrtEhISap0krJqFYeqo3ZyViATRmmjxwlxRUYH09HQsWrSIPycUChEfH4/k5ORan1NeXo7y8nL+2JAJUWN0W429e/ciMzMT4eHh2LlzJ06dOoU7d+6AMQZnZ2ejibqqCF1cwTr0gSjnCoQCAALzf6hotVoAwP3793H48GHk5uaCMYbS0lIoFAoA/9vfLyMjg7cnMzOT39rKUPJTLBYjKioKiYmJNeLRjz/+eK3iWLXdtm3bkJaWhmvXriE0NLTGPaKiovDEE08gMTGR76++1xf432tssKulYfgMNeRna8JRfbYmLV6Y79+/D47jamQo+Pv749KlS7U+Z9WqVVi2bFmN8yqVqs4X3LDTtUQiwZUrVyCVSiGVShEcHFyvfUzoBNYhFj7uThBy94COHRvpWd0IBAJIpVIEBATwAiuRSBAdHY2rV68iIiICV65c4VPVDGlrKpUK2dnZfJH+mJgYZGdno6CggB/ZGtoWFxfXef+oqMryo3K5HFKpFHK5HK6urggLC+PvIZFIcOPGDRQWFpqUzWF4blW7WhKMMf61EwgENrameXBEnw2lFaxFixdmc1i0aBHmz5/PH6vVagQFBUEmk8HT07PO5xlGhVFRUQgLC8OFCxeQlpYGnU4HFxcXo1E4j0AIcakAQokON2/c4Cfy6kMgEMDJyYkfHQsEAj5nmTEGPz8/REZGwtnZGRKJBG3btkXHjh2h0WiMRsC1VZUz+NC5c2e+vZeXFwDUeF5dGEbAhpHxiBEjjEbYhnuEhoaiTZs2Jv2xVn2NDXa1JAxf7DKZzGFEyhF9tjYC1sJ/f1RUVMDNzQ07d+40iotOmTIFRUVF2Lt3b4N9qNVqyGQyqFSqeoUZQA3B4zgOZWVl8PDwAMdxyMnJQXZ2Nm7dUeLcfT1yUv+LsmI1JBIJ7ty5A4VCgYiICLRr1w4+Pj6QyWTw8vKCRCKBh4cHvLy84OLiAhcXFxQXF/P3kkgkqKio4I9dXFyM4s+mVImrK8ZsDnX1odPpUFxcbNYfa0uOMTPGoFKpHEqkHNFnlUoFuVzeKM0whxY/YnZxcUGvXr1w+PBhXpj1ej0OHz6M2bNnW/x+1WtNiEQifvGISCRCly5d0C64I37MyMVYqSvi3/o7AJj1wZXJZEbHEomkTltqO27IB0uIX119NKXvlirKBGEpWrwwA8D8+fMxZcoU9O7dG3369MFnn32GkpISPkujOdFyevyYkQs/qSviu/nxS6sJgiAaS6sQ5meffRb37t3DkiVLoFQqERUVhQMHDjRqybKlcRIKENvRC519PRzmZx1BEJalVQgzAMyePdsqoYvGotFy+OVcHoaF+aOLn9RmdhAE0fJpFSv/bI1Gy+HHjFxIXZ3hKWk133UEQdgIEuYmotcz7D592yimTBAE0RRoeNcEGGMQCgUYFOKLAJkriTJBEBaBRsxmotFy2PbnLeSrNQiUS0iUCYKwGCTMZmCIKXt7iOEnFdvaHIIgWhkkzCbCGMPPZ/MopkwQhNWgGLMJaDk9nEVCDOvmB5nEmUSZIAirQCPmRqLRctiedgvX7hVD7uZCokwQhNUgYW4Ehpiyn9QVnXzcbW0OQRCtHBLmRnDk0l2KKRME0WxQjLkeNFoOziIhhoT6wdVZSKJMEESzQCPmOjCELy7cUUPiIiJRJgii2SBhroWqMeWebS1fBJsgCKI+SJhrIeWvBxRTJgjCZlCMuQoaLQfnCg79u/hAJBSQKBMEYRNoxFyFvWdu42xuEZxENNFHEITtIGGugq+HK/p0bHk7MxME0bqgUAb+t/16rwAxHj58aJX+1Wo1BALHCY84os+AY/rtiD6r1WoAsNp+niTMAC/G7du3t7ElBEG0JB48eFBjN3tLIGC0hTP0ej3u3LkDqVRqlW98tVqNoKAg3Lp1C56ejpF+54g+A47ptyP6rFKp0L59exQWFkIul1u8fxoxAxAKhWjXrp3V7+Pp6ekwH1wDjugz4Jh+O6LPQqF1pulo8o8gCMLOIGEmCIKwM0iYmwGxWIylS5dCLHacbagc0WfAMf0mny0PTf4RBEHYGTRiJgiCsDNImAmCIOwMEmaCIAg7g4SZIAjCziBhbgbWrVuHDh06wNXVFbGxsfjjjz9sbZLFeO+99/gaCYZHWFgYf12j0WDWrFnw9vaGh4cHxo8fj/z8fBtabDrHjh1DQkICAgMDIRAIsGfPHqPrjDEsWbIEAQEBkEgkiI+Px5UrV4zaFBQUYNKkSfD09IRcLsf06dNRXFzcjF6YRkM+T506tcb7PnLkSKM2Lc3nVatW4ZFHHoFUKoWfnx/Gjh2L7OxsozaN+Tzn5ORgzJgxcHNzg5+fH9544w3odDqTbCFhtjLbtm3D/PnzsXTpUmRkZCAyMhIjRozA3bt3bW2axejRowfy8vL4x4kTJ/hr8+bNw759+7Bjxw4cPXoUd+7cwbhx42xoremUlJQgMjIS69atq/X66tWr8cUXX2D9+vVITU2Fu7s7RowYAY1Gw7eZNGkSsrKycPDgQezfvx/Hjh3Diy++2FwumExDPgPAyJEjjd73LVu2GF1vaT4fPXoUs2bNQkpKCg4ePAitVovhw4ejpKSEb9PQ55njOIwZMwYVFRU4deoUvv/+e2zcuBFLliwxzRhGWJU+ffqwWbNm8cccx7HAwEC2atUqG1plOZYuXcoiIyNrvVZUVMScnZ3Zjh07+HMXL15kAFhycnIzWWhZALDdu3fzx3q9nikUCvbxxx/z54qKiphYLGZbtmxhjDF24cIFBoD9+eeffJv//ve/TCAQsNu3bzeb7eZS3WfGGJsyZQpLTEys8zkt3WfGGLt79y4DwI4ePcoYa9zn+ZdffmFCoZAplUq+zVdffcU8PT1ZeXl5o+9NI2YrUlFRgfT0dMTHx/PnhEIh4uPjkZycbEPLLMuVK1cQGBiITp06YdKkScjJyQEApKenQ6vVGvkfFhaG9u3btxr/r1+/DqVSaeSjTCZDbGws72NycjLkcjl69+7Nt4mPj4dQKERqamqz22wpkpKS4Ofnh9DQUMycORMPHjzgr7UGn1UqFQDAy6uyRntjPs/JyckIDw+Hv78/32bEiBFQq9XIyspq9L1JmK3I/fv3wXGc0ZsEAP7+/lAqlTayyrLExsZi48aNOHDgAL766itcv34dAwYMwMOHD6FUKuHi4lKj+lZr8t/gR33vsVKphJ+fn9F1JycneHl5tdjXYeTIkfjhhx9w+PBhfPTRRzh69ChGjRoFjuMAtHyf9Xo95s6di379+qFnz54A0KjPs1KprPWzYLjWWKi6HNEkRo0axf8/IiICsbGxCA4Oxvbt2yGRSGxoGWFNJkyYwP8/PDwcERER6Ny5M5KSkjBs2DAbWmYZZs2ahfPnzxvNlzQnNGK2Ij4+PhCJRDVmbfPz86FQKGxklXWRy+UICQnB1atXoVAoUFFRgaKiIqM2rcl/gx/1vccKhaLGZK9Op0NBQUGreR06deoEHx8fXL16FUDL9nn27NnYv38/fv/9d6NywI35PCsUilo/C4ZrjYWE2Yq4uLigV69eOHz4MH9Or9fj8OHDiIuLs6Fl1qO4uBjXrl1DQEAAevXqBWdnZyP/s7OzkZOT02r879ixIxQKhZGParUaqampvI9xcXEoKipCeno63+bIkSPQ6/WIjY1tdputQW5uLh48eICAgAAALdNnxhhmz56N3bt348iRI+jYsaPR9cZ8nuPi4nDu3DmjL6WDBw/C09MT3bt3N8kYwops3bqVicVitnHjRnbhwgX24osvMrlcbjRr25J5/fXXWVJSErt+/To7efIki4+PZz4+Puzu3buMMcZefvll1r59e3bkyBGWlpbG4uLiWFxcnI2tNo2HDx+y06dPs9OnTzMAbM2aNez06dPs5s2bjDHGPvzwQyaXy9nevXvZ2bNnWWJiIuvYsSMrKyvj+xg5ciSLjo5mqamp7MSJE6xr165s4sSJtnKpQerz+eHDh2zBggUsOTmZXb9+nR06dIjFxMSwrl27Mo1Gw/fR0nyeOXMmk8lkLCkpieXl5fGP0tJSvk1Dn2edTsd69uzJhg8fzjIzM9mBAweYr68vW7RokUm2kDA3A2vXrmXt27dnLi4urE+fPiwlJcXWJlmMZ599lgUEBDAXFxfWtm1b9uyzz7KrV6/y18vKytgrr7zC2rRpw9zc3NiTTz7J8vLybGix6fz+++8MQI3HlClTGGOVKXOLFy9m/v7+TCwWs2HDhrHs7GyjPh48eMAmTpzIPDw8mKenJ5s2bRp7+PChDbxpHPX5XFpayoYPH858fX2Zs7MzCw4OZjNmzKgx2GhpPtfmLwD23Xff8W0a83m+ceMGGzVqFJNIJMzHx4e9/vrrTKvVmmQLlf0kCIKwMyjGTBAEYWeQMBMEQdgZJMwEQRB2BgkzQRCEnUHCTBAEYWeQMBMEQdgZJMwEQRB2BgkzQRCEnUHCTNg9U6dOxdixY/njwYMHY+7cuc1uR1JSEgQCQY0iNrZi4MCB2Lx5c5P66Nu3L3788UcLWURYChJmwiyq7vnm4uKCLl26YPny5SbvbWYOu3btwvvvv9+ots0tph06dOBfF3d3d8TExGDHjh389dLSUixatAidO3eGq6srfH19MWjQIOzdu5dv05gvnp9++gn5+flG5Tfnz58PLy8vBAUFYdOmTUbtd+zYgYSEhBr9vPvuu3jrrbeg1+vN9JiwBiTMhNkY9ny7cuUKXn/9dbz33nv4+OOPa21bUVFhsft6eXlBKpVarD9Ls3z5cuTl5eH06dN45JFH8Oyzz+LUqVMAgJdffhm7du3C2rVrcenSJRw4cABPPfWU0e4fjeGLL77AtGnTIBRW/gnv27cPmzdvxm+//YbVq1fj73//O+7fvw+gcieOd955p9b9+0aNGoWHDx/iv//9bxO9JiwJCTNhNmKxGAqFAsHBwZg5cybi4+Px008/Afhf+OGDDz5AYGAgQkNDAQC3bt3CM888A7lcDi8vLyQmJuLGjRt8nxzHYf78+ZDL5fD29sabb76J6uVcqo8oy8vLsXDhQgQFBUEsFqNLly749ttvcePGDQwZMgQA0KZNGwgEAkydOhVAZfnVVatWoWPHjpBIJIiMjMTOnTuN7vPLL78gJCQEEokEQ4YMMbKzPqRSKRQKBUJCQrBu3TpIJBLs27cPQOVI9+2338bo0aPRoUMH9OrVC6+++ipeeOGFxr7suHfvHo4cOWI0Ar548SIGDx6M3r17Y+LEifD09MT169cBAG+++SZmzpyJ9u3b1+hLJBJh9OjR2Lp1a6PvT1gfEmbCYkgkEqOR8eHDh5Gdnc3vkqzVajFixAhIpVIcP34cJ0+ehIeHB0aOHMk/7x//+Ac2btyIDRs24MSJEygoKMDu3bvrve/kyZOxZcsWfPHFF7h48SK+/vpreHh4ICgoiI+fZmdnIy8vD59//jmAyq3qf/jhB6xfvx5ZWVmYN28enn/+eRw9ehRA5RfIuHHjkJCQgMzMTPz973/HW2+9ZfJr4uTkBGdnZ94/hUKBX375BQ8fPjS5LwMnTpyAm5sbunXrxp+LjIxEWloaCgsLkZ6ejrKyMnTp0gUnTpxARkYG5syZU2d/ffr0wfHjx822h7ACTa6VRzgkVXdJ1uv17ODBg0wsFrMFCxbw1/39/Y12Bv73v//NQkNDmV6v58+Vl5cziUTCfv31V8YYYwEBAWz16tX8da1Wy9q1a2e0I/OgQYPYa6+9xhhjLDs7mwFgBw8erNVOQ/nKwsJC/pxGo2Fubm7s1KlTRm2nT5/O1wtetGgR6969u9H1hQsX1uirOsHBwezTTz/lfVu5ciUDwPbv388YY+zo0aOsXbt2zNnZmfXu3ZvNnTuXnThxwqiPqv7Vxqeffso6depU4/zSpUtZ586dWc+ePdmuXbtYeXk569mzJ0tLS2Nr165lISEh7NFHH2Xnz583et7evXuZUChkHMfVeU+ieaE9/wiz2b9/Pzw8PKDVaqHX6/Hcc8/hvffe46+Hh4fDxcWFPz5z5gyuXr1aIz6s0Whw7do1qFQq5OXlGe1w4eTkhN69e9cIZxjIzMyESCTCoEGDGm331atXUVpaiscee8zofEVFBaKjowFUhgaq77TR2F1XFi5ciHfffRcajQYeHh748MMPMWbMGACVmRR//fUXUlJScOrUKRw+fBiff/45li1bhsWLFzeq/7KyMri6utY4/9577xm9/suWLUN8fDycnZ2xYsUKnDt3Dvv378fkyZONdhaRSCTQ6/UoLy+nfRrtBBJmwmyGDBmCr776Ci4uLggMDISTk/HHyd3d3ei4uLgYvXr1qpExAAC+vr5m2WCOkBQXFwMAfv75Z7Rt29bomlgsNsuOqrzxxhuYOnUqPDw84O/vD4FAYHTd2dkZAwYMwIABA7Bw4UKsWLECy5cvx8KFC42+yOrCx8cHhYWF9ba5dOkS/vOf/+D06dPYsGEDBg4cCF9fXzzzzDN44YUX8PDhQ/4LsqCgAO7u7iTKdgQJM2E27u7u6NKlS6Pbx8TEYNu2bfDz84Onp2etbQICApCamoqBAwcCqNzAMz09HTExMbW2Dw8Ph16vx9GjRxEfH1/jukHoOI7jz3Xv3h1isRg5OTl1jrS7devGT2QaSElJadhJVAqnKa9L9+7dodPpoNFoGiXM0dHRUCqVKCwsRJs2bWpcZ4zhpZdewpo1a+Dh4QGO46DVagGA/7fq63H+/Hn+lwJhH9DkH9FsTJo0CT4+PkhMTMTx48dx/fp1JCUlYc6cOcjNzQUAvPbaa/jwww+xZ88eXLp0Ca+88kq9OcgdOnTAlClT8MILL2DPnj18n9u3bwcABAcHQyAQYP/+/bh37x6Ki4shlUqxYMECzJs3D99//z2uXbuGjIwMrF27Ft9//z2AyrS2K1eu4I033kB2djY2b96MjRs3Nvk1GDx4ML7++mukp6fjxo0b+OWXX/D2229jyJAhdX5ZVSc6Oho+Pj44efJkrdf/9a9/wdfXl8/a6NevH44cOYKUlBR8+umn6N69O+RyOd/++PHjGD58eJN9IyyIrYPcRMuk6uSfKdfz8vLY5MmTmY+PDxOLxaxTp05sxowZTKVSMcYqJ/tee+015unpyeRyOZs/fz6bPHlynZN/jFXuwzZv3jx+78EuXbqwDRs28NeXL1/OFAoFEwgERvv0ffbZZyw0NJQ5OzszX19fNmLECHb06FH+efv27WNdunRhYrGYDRgwgG3YsMGkyb/aWLlyJYuLi2NeXl7M1dWVderUic2ZM4fdv3+/Tv9q480332QTJkyocV6pVLLg4GB2+/Zto/PLli1jXl5eLCwsjKWmpvLnc3NzmbOzM7t161a99yOaF9rzjyBaIEqlEj169EBGRgaCg4PN7mfhwoUoLCzEN998Y0HriKZCoQyCaIEoFAp8++23yMnJaVI/fn5+jV7eTjQfNGImCIKwM2jETBAEYWeQMBMEQdgZJMwEQRB2BgkzQRCEnUHCTBAEYWeQMBMEQdgZJMwEQRB2BgkzQRCEnUHCTBAEYWf8P9IdpiW9kkKAAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 400x350 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 1, figsize=(4, 3.5))\n",
    "lims = (-20, 200)\n",
    "axes.scatter(pred['coef'], obs, s=5, c='black', lw=0, alpha=0.5)\n",
    "axes.plot(lims, lims, lw=0.75, linestyle='--', alpha=0.5)\n",
    "axes.grid(alpha=0.2)\n",
    "axes.text(0.05, 0.98, '$R^2$=' + '{:.2f}\\nCalibration = {:.2f}'.format(r2, calibration),\n",
    "          transform=axes.transAxes, ha='left', va='top')\n",
    "axes.set(xlabel='Predicted PSI (%)', xlim=lims, \n",
    "         ylabel='Measured PSI (%)', ylim=lims, aspect='equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8aa988d",
   "metadata": {},
   "source": [
    "This scatterplot demonstrates the relatively good performance of the model in predicting the phenotype of held-out sequences in terms of $R^2$. Importantly, the model provides well-calibrated uncertainties, as the posterior distribution includes the actual measurements the expected number of times. This offers an estimate of how confident we can be about predictions for specific sequences with unknown phenotypes, aiding in decision-making."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78225021",
   "metadata": {},
   "source": [
    "### Compare different models\n",
    "\n",
    "When analyzing a new dataset, we often want to fit different models and compare them. This is easy to do by just providing a different kernel to the `EpiK` model. Here, we fit models with kernels of increasing complexity and compute the $R^2$ of the predictions in held out data to track model performance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "be816a11",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Optimizing hyperparameters: 100%|██████████| 100/100 [00:39<00:00,  2.53it/s, MLL=-9160.443, Mem(alloc/res)=0.00/0.00MB]\n",
      "Optimizing hyperparameters: 100%|██████████| 100/100 [00:44<00:00,  2.26it/s, MLL=-8903.045, Mem(alloc/res)=0.00/0.00MB]\n",
      "Optimizing hyperparameters: 100%|██████████| 100/100 [00:46<00:00,  2.16it/s, MLL=-8276.521, Mem(alloc/res)=0.00/0.00MB]\n"
     ]
    },
    {
     "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>kernel</th>\n",
       "      <th>r2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Exponential</td>\n",
       "      <td>0.559659</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Connectedness</td>\n",
       "      <td>0.577626</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Jenga</td>\n",
       "      <td>0.653812</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          kernel        r2\n",
       "0    Exponential  0.559659\n",
       "1  Connectedness  0.577626\n",
       "2          Jenga  0.653812"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kernels = [ExponentialKernel, ConnectednessKernel, JengaKernel]\n",
    "labels = ['Exponential', 'Connectedness', 'Jenga']\n",
    "r2s = []\n",
    "for label, kernel_type in zip(labels, kernels):\n",
    "    kernel = kernel_type(n_alleles=4, seq_length=8)\n",
    "    model = EpiK(kernel, track_progress=True, max_n_lanczos_iterations=200)\n",
    "    model.set_data(train_x, train_y, train_y_var)\n",
    "    model.fit(n_iter=100, learning_rate=0.1)\n",
    "    test_y_pred = model.predict(test_x, calc_variance=False)['coef']\n",
    "    r2s.append(pearsonr(test_y_pred, test_y)[0] ** 2)\n",
    "results = pd.DataFrame({'kernel': labels, 'r2': r2s})\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b65fc14",
   "metadata": {},
   "source": [
    "Thus, this analysis shows that mutations at different sites affect the predictability of other mutations in a different manner, given the improved performance of the `ConnectednessKernel` over the `ExponentialKernel`. However, the `JengaKernel` provides even better predictive power, suggesting that the decay in predictability of mutational effects also depends on the specific alleles that are mutated at each site.\n",
    "\n",
    "### Interpreting Jenga model parameters\n",
    "\n",
    "Once we have settled with the Jenga model, we aim use the inferred hyperparameters of the kernel function to gain some insigths about the sequence-function relationship we are studying. We can extract the allele and site specific decay rates with the `get_decay_rates` method as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "45721d81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <th>G</th>\n",
       "      <th>U</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>-3</th>\n",
       "      <td>0.020853</td>\n",
       "      <td>0.025342</td>\n",
       "      <td>0.438904</td>\n",
       "      <td>0.170491</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-2</th>\n",
       "      <td>0.288113</td>\n",
       "      <td>0.039503</td>\n",
       "      <td>0.612494</td>\n",
       "      <td>0.129223</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-1</th>\n",
       "      <td>0.076758</td>\n",
       "      <td>0.081544</td>\n",
       "      <td>0.820145</td>\n",
       "      <td>0.430954</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>+2</th>\n",
       "      <td>NaN</td>\n",
       "      <td>0.141537</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.713322</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>+3</th>\n",
       "      <td>0.774562</td>\n",
       "      <td>0.136644</td>\n",
       "      <td>0.406120</td>\n",
       "      <td>0.128953</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>+4</th>\n",
       "      <td>0.284240</td>\n",
       "      <td>0.052869</td>\n",
       "      <td>0.256599</td>\n",
       "      <td>0.123743</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>+5</th>\n",
       "      <td>0.016301</td>\n",
       "      <td>0.034593</td>\n",
       "      <td>0.713966</td>\n",
       "      <td>0.014079</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>+6</th>\n",
       "      <td>0.013091</td>\n",
       "      <td>0.012482</td>\n",
       "      <td>0.041129</td>\n",
       "      <td>0.073866</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           A         C         G         U\n",
       "-3  0.020853  0.025342  0.438904  0.170491\n",
       "-2  0.288113  0.039503  0.612494  0.129223\n",
       "-1  0.076758  0.081544  0.820145  0.430954\n",
       "+2       NaN  0.141537       NaN  0.713322\n",
       "+3  0.774562  0.136644  0.406120  0.128953\n",
       "+4  0.284240  0.052869  0.256599  0.123743\n",
       "+5  0.016301  0.034593  0.713966  0.014079\n",
       "+6  0.013091  0.012482  0.041129  0.073866"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alleles = ['A', 'C', 'G', 'U']\n",
    "positions = ['-3', '-2', '-1', '+2', '+3', '+4', '+5', '+6']\n",
    "delta = kernel.get_delta().detach().numpy()\n",
    "delta = pd.DataFrame(delta, index=positions, columns=alleles)\n",
    "delta.loc['+2', ['A', 'G']] = np.nan\n",
    "delta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6941a926",
   "metadata": {},
   "source": [
    "We can visualize these values in a heatmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "422ee2f4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAADfCAYAAABbArz6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2cklEQVR4nO3dd1gUV/s38O8uZUE6+NAMzQYWrCgiPooGxfKzYowJKpbYsRMVE4UYu4kxlggqgj6vPbao0YgIohGxYEuiCBZAFNQgVVnKnvcPwsSFXWCXZQfw/uSa6wpnZs7cs4S9c86cOUfAGGMghBBCeCTkOwBCCCGEkhEhhBDeUTIihBDCO0pGhBBCeEfJiBBCCO8oGRFCCOEdJSNCCCG8o2RECCGEd5SMCCGE8I6SESGEEN5RMiKEkHokJiYGgwcPhrW1NQQCAY4fPy61nzGGZcuWwcrKCrq6uvD09ERiYqLUMZmZmfDx8YGhoSGMjY0xadIk5OXlqfEuKqJkRAgh9Uh+fj7at2+PrVu3yty/bt06bNq0CcHBwYiLi4Oenh68vLxQUFDAHePj44M///wTEREROHXqFGJiYjBlyhR13YJMApoolRBC6ieBQIBjx45h2LBhAEpbRdbW1liwYAH8/f0BANnZ2bCwsEB4eDhGjx6N+/fvo3Xr1rh+/TpcXFwAAGfPnsXAgQPx7NkzWFtb83IvmrxcVcUkEgmeP38OAwMDCAQCvsMhhDQQjDHk5ubC2toaQqFyHUlVfT8xxvD333/D1NRU6hoikQgikUihaz158gTp6enw9PTkyoyMjODq6orY2FiMHj0asbGxMDY25hIRAHh6ekIoFCIuLg7Dhw9X4i5rrkEko+fPn8PGxobvMAghDVRqaio++ugjpc5V9vspMDAQQUFBCp2Tnp4OALCwsJAqt7Cw4Palp6fD3Nxcar+mpiZMTU25Y/jQIJKRgYEBAODBoxQYGBjyHA1w7Wkm3yFw2lrx/3mUef0sseqD1MQ79AnfIQAAtk3txncInC72pnyHwKkr/Ru5uTlo2dSW+45RRtm52m0nQKChXWE/KylE4R9hSElJgZGREVeuaKuovmsQyais6WtgYAhDQ/6/fPX0i/gOgWNQBz6PMgX6+nyHwBFqN+I7BACAnr7yX3KqVhf+dsrUlWRUpibd/2XnCjR1INCQkWAEpV1zRkZGNf4dWFpaAgAyMjJgZWXFlWdkZKBDhw7cMS9fvpQ6r7i4GJmZmdz5fKDRdIQQog5CDfmbijg4OMDS0hKRkZFcWU5ODuLi4uDm5gYAcHNzQ1ZWFm7evMkdc+HCBUgkEri6uqosFkU1iJYRIYTUeQI5iUeiWDLKy8tDUlIS9/OTJ09w+/ZtmJqawtbWFnPnzsWKFSvQokULODg4YOnSpbC2tuZG3LVq1Qr9+/fH5MmTERwcjKKiIvj5+WH06NG8jaQDKBkRQoh6CASlm6xyBdy4cQO9e/fmfp4/fz4AwNfXF+Hh4Vi4cCHy8/MxZcoUZGVloUePHjh79ix0dHS4c/bu3Qs/Pz98/PHHEAqF8Pb2xqZNm5S7LxWhZEQIIeogr0tOwW46Dw8PVPZ6qEAgwPLly7F8+XK5x5iammLfvn0KXbe2UTIihBB1UFEyaqgoGRFCiDoIhNzIuQrlhJIRIYSohYYQ0JDRCmKUjABKRoQQoh7UTVcpSkaEEKIO1E1XKUpGhBCiDtQyqhQlI0IIUQdKRpWqM+3D2NhYaGhoYNCgQXyHQgghqlf20qusjdSdZBQaGopZs2YhJiYGz58/5zscQghRLYGceekE1DIC6kgyysvLw8GDBzF9+nQMGjQI4eHhfIdECCGqJdQAhJoyNkpGQB1JRocOHYKTkxMcHR0xZswY7Nq1q9LpLgghpN6hbrpK1YlkFBoaijFjxgAA+vfvj+zsbFy8eFHu8WKxGDk5OVIbIYTUaWpYQqI+4z0ZJSQk4Nq1a/jss88AlC5/++mnnyI0NFTuOatXr4aRkRG30ZLjhJA6j5JRpXgf2h0aGori4mKpdTQYYxCJRNiyZYvUMrxlAgICuGnTgdLFoyghEULqMoFAIHvFWOqmA8BzMiouLsaePXvw/fffo1+/flL7hg0bhv3792PatGkVzhOJRB/c+vCEkPpNIBRAIJSReGSVfYB4TUanTp3CmzdvMGnSpAotIG9vb4SGhspMRoQQUt8IhUIIhBWfjDAZZR8iXj+F0NBQeHp6yuyK8/b2xo0bN3D37l0eIiOEENUq66aTtRGeW0YnT56Uu69r1640vJsQ0mBQN13leB/AQAghHwKhUCCnm46SEUDJiBBC1EIAeV1ylIyAOvCeESGEfAjKuulkbYooKSnB0qVL4eDgAF1dXTRr1gzffvut1GMNxhiWLVsGKysr6OrqwtPTE4mJiSq5D4lEgqioKCxfvhyTJk3CZ599htmzZyMsLAypqalK10vJiBBC1EAgFEIoY5PVdVeZtWvXYtu2bdiyZQvu37+PtWvXYt26ddi8eTN3zLp167Bp0yYEBwcjLi4Oenp68PLyQkFBgdLxv3v3DitWrICNjQ0GDhyIM2fOICsrCxoaGkhKSkJgYCAcHBwwcOBAXL16VeH6qZuOEELUQN7IOUVH0125cgVDhw7lltuxt7fH/v37ce3aNQClraKNGzfi66+/xtChQwEAe/bsgYWFBY4fP47Ro0crFX/Lli3h5uaGHTt2oG/fvtDS0qpwTHJyMvbt24fRo0fjq6++wuTJk6tdP7WMCCFEDarqpis/36ZYLJZZT/fu3REZGYmHDx8CAO7cuYPLly9jwIABAIAnT54gPT0dnp6e3DlGRkZwdXVFbGys0vGfO3cOhw4dwsCBA2UmIgCws7NDQEAAEhMT0adPH4Xqp5YRIYSoQVUto/JTmgUGBiIoKKjC8YsXL0ZOTg6cnJygoaGBkpISrFy5Ej4+PgCA9PR0AICFhYXUeRYWFtw+ZbRq1arax2ppaaFZs2YK1U/JiBBC1KDsGZGMHQCA1NRUGBoacsXypjw7dOgQ9u7di3379qFNmza4ffs25s6dC2tra/j6+tZK7PIUFxcjJCQE0dHRKCkpgbu7O2bOnAkdHR2F66JkRAghaiBv5FxZmaGhoVQykufLL7/E4sWLuWc/zs7OSE5OxurVq+Hr6wtLS0sAQEZGBqysrLjzMjIy0KFDBxXcyb9mz56Nhw8fYsSIESgqKsKePXtw48YN7N+/X+G6KBkRQogaqGoAw9u3byu0sDQ0NCCRSAAADg4OsLS0RGRkJJd8cnJyEBcXh+nTpysX/D+OHTuG4cOHcz+fO3cOCQkJ0NAoXQbDy8sL3bp1U6puSkaEEKIGVbWMqmvw4MFYuXIlbG1t0aZNG9y6dQsbNmzAxIkTS+sTCDB37lysWLECLVq0gIODA5YuXQpra2sMGzasRvewa9cu7N69Gz/99BOsra3RqVMnTJs2Dd7e3igqKsKOHTvQpUsXpeqmZEQIIWpQ1TOj6tq8eTOWLl2KGTNm4OXLl7C2tsbUqVOxbNky7piFCxciPz8fU6ZMQVZWFnr06IGzZ88q9SznfSdPnsTBgwfh4eGBWbNmYfv27fj222/x1Vdfcc+MZA26qA4BawCzkebk5MDIyAhpL7Oq1eda22If/813CJx2TSrOiM6XVykJfIfAGbDtEd8hAADCZ7nzHQKnW1MzvkPg1JUJcnJycmD1H2NkZ2cr/d1S9v1kM/UghKJGFfZLxG+RGvJpja6hbllZWVi4cCHu3LmD4OBgdOzYscZ10ntGhBCiBqqaDqguMDY2xvbt27F+/XqMGzcOX375ZY1mdwAaWDcd++cfvnWxN+E7BE5BkYTvEDg65oq9d1CbooPqRiwmjbT5DoFTIuH/b6eMZj38gq6Kqrrp+JSSkgJ/f3/cv38f7dq1w3fffYebN29i5cqVaN++PTZu3Mi9fKuo+vMpEEJIPSYQyN/qi3HjxkEoFGL9+vUwNzfH1KlToa2tjW+++QbHjx/H6tWrMWrUKKXqblAtI0IIqasEwtI1jcpj9ahJcOPGDdy5cwfNmjWDl5cXHBwcuH2tWrVCTEwMtm/frlTdlIwIIUQNhEKBnGRUf5pGnTt3xrJly+Dr64vz58/D2dm5wjFTpkxRqu56lJMJIaT+KktGsrb6Ys+ePRCLxZg3bx7S0tIQEhKisrqpZUQIIWrQEFpGdnZ2+Pnnn2ulbmoZEUKIGtT3llF+fn6tHk/JiBBC1KBsbjpZW33QvHlzrFmzBi9evJB7DGMMERERGDBgADZt2qRQ/dRNRwghaiAUyOmmqyfJKDo6GkuWLEFQUBDat28PFxcXWFtbQ0dHB2/evMFff/2F2NhYaGpqIiAgAFOnTlWofkpGhBCiBvX9mZGjoyOOHDmClJQUHD58GJcuXcKVK1fw7t07NG7cGB07dsSOHTswYMAAbhZvRVAyIoQQNZD3gms9aRhxbG1tsWDBAixYsECl9VIyIoQQNajvLaPaRsmIEELUQCAnGUkoGQGoI6Pp0tPTMWvWLDRt2hQikQg2NjYYPHgwIiMj+Q6NEEJUor6PpqttvLeMnj59Cnd3dxgbG2P9+vVwdnZGUVERfvvtN8ycORMPHjzgO0RCCKkxed109eU9o9rGezKaMWMGBAIBrl27Bj09Pa68TZs23DK6hBBS31Eyqhyv3XSZmZk4e/YsZs6cKZWIyhgbG6s/KEIIqQUNrZvu0qVLGDNmDNzc3JCWlgYA+N///ofLly8rVR+vySgpKQmMMTg5OfEZBiGE1Lr6Ph3Q+44cOQIvLy/o6uri1q1bEIvFAIDs7GysWrVKqTp5TUaMKbeypFgsRk5OjtRGCCF1WUNKRitWrEBwcDB27NgBLS0trtzd3R3x8fFK1clrMmrRogUEAoHCgxRWr14NIyMjbrOxsamlCAkhRDUEkLPSqxJ1paWlYcyYMTAzM4Ouri6cnZ1x48YNbj9jDMuWLYOVlRV0dXXh6emJxMREld1LQkICevbsWaHcyMgIWVlZStXJazIyNTWFl5cXtm7dKnOGV3k3FRAQgOzsbG5LTU2t5UgJIaRmNIQCuZsi3rx5A3d3d2hpaeHMmTP466+/8P3338PExIQ7Zt26ddi0aROCg4MRFxcHPT09eHl5oaCgQCX3YmlpiaSkpArlly9fRtOmTZWqk/fRdFu3boW7uzu6du2K5cuXo127diguLkZERAS2bduG+/fvVzhHJBJBJBLxEC0hhChHXuJRdAaGtWvXwsbGBmFhYVzZ+8t/M8awceNGfP311xg6dCiA0kXxLCwscPz4cYwePVrJO/jX5MmTMWfOHOzatQsCgQDPnz9HbGws/P39sXTpUqXqVKplVFxcjPPnzyMkJAS5ubkAgOfPnyMvL0/hupo2bYr4+Hj07t0bCxYsQNu2bdG3b19ERkZi27ZtyoRHCCF1j0D2iLqyfrryz8HLBgWU98svv8DFxQWffPIJzM3NuQlKyzx58gTp6enw9PTkyoyMjODq6orY2FiV3MrixYvx+eef4+OPP0ZeXh569uyJL774AlOnTsWsWbOUqlPhllFycjL69++PlJQUiMVi9O3bFwYGBli7di3EYjGCg4MVDsLKygpbtmzBli1bFD6XEELqA6FAAKGMYdxlZeWffQcGBiIoKKjC8Y8fP8a2bdswf/58LFmyBNevX8fs2bOhra0NX19fpKenAwAsLCykzrOwsOD21ZRAIMBXX32FL7/8EklJScjLy0Pr1q2hr6+vdJ0Kt4zmzJkDFxcXvHnzBrq6ulz58OHDafoeQgiRo6pnRqmpqVLPwgMCAmTWI5FI0KlTJ6xatQodO3bElClTMHnyZKUaAsqaOHEicnNzoa2tjdatW6Nr167Q19dHfn6+0pMVKJyMLl26hK+//hra2tpS5fb29tyLT4QQQqTJHEn33rIShoaGUpu85+JWVlZo3bq1VFmrVq2QkpICoHRwAQBkZGRIHZORkcHtq6ndu3fj3bt3FcrfvXuHPXv2KFWnwslIIpGgpKSkQvmzZ89gYGCgVBCEENLQqWo0nbu7OxISEqTKHj58CDs7OwClgxksLS2leqpycnIQFxcHNze3Gt1DTk4OsrOzwRhDbm6u1DOuN2/e4Ndff4W5ublSdSv8zKhfv37YuHEjtm/fDqC07zAvLw+BgYEYOHCgUkEQQkhDp6q56ebNm4fu3btj1apVGDVqFK5du4bt27dLfSfPnTsXK1asQIsWLeDg4IClS5fC2toaw4YNq9E9GBsbcwMvWrZsWWG/QCDAN998o1TdCiej77//Hl5eXmjdujUKCgrw+eefIzExEY0bN8b+/fuVCoIQQhq6qgYwVFeXLl1w7NgxBAQEYPny5XBwcMDGjRvh4+PDHbNw4ULk5+djypQpyMrKQo8ePXD27Fno6OjU6B6ioqLAGEOfPn1w5MgRmJqacvu0tbVhZ2cHa2trpeoWMCXm5CkuLsaBAwdw9+5d5OXloVOnTvDx8ZEa0KBOOTk5MDIywrOXb2BoaMhLDO8rkSg3zVFtKCiS8B0CJ6+gmO8QOHVlbkqTRtpVH6Qmmhp15EMBoFlHpsjJycmB1X+MkZ2drfR3S9n3k3fwJWjpVhxtVvQuD0em/bdG11C35ORk2NjYQChU3bwJSr30qqmpiTFjxqgsCEIIaejkPR+qjyu9lj2fevv2LVJSUlBYWCi1v127dgrXWa1k9Msvv1S7wiFDhigcBCGEfAjqSou8pl69eoUJEybgzJkzMvfLGuRWlWolo+o+9BIIBEoFQQghDZ28lpGio+nqgrlz5yIrKwtxcXHw8PDAsWPHkJGRgRUrVuD7779Xqs5qJSOJpO48dyCEkPpIVQMY6oILFy7gxIkTcHFxgVAohJ2dHfr27QtDQ0OsXr0agwYNUrjOGj19UtUMsIQQ0tCVJSNZW32Tn5/PvU9kYmKCV69eAQCcnZ3Vt55RSUkJvv32WzRp0gT6+vp4/PgxAGDp0qUIDQ1VKghCCGnoGtLieo6OjtyLt+3bt0dISAjS0tIQHBwMKysrpepUOBmtXLkS4eHhWLdundSUQG3btsXOnTuVCoIQQho6Vc3AUBfMmTMHL168AFA6oeuZM2dga2uLTZs2Kb3suMJDu/fs2YPt27fj448/xrRp07jy9u3bK7xiKyGEfCjen4eufHl98/6rPZ07d0ZycjIePHgAW1tbNG7cWKk6FU5GaWlpaN68eYVyiUSCoqIipYJQleDYJ9DR439+vJiEv/kOgXPz2hO+Q+BE+7flOwROVLYe3yEAADqZG/MdAmfq7htVH6Qml7/qw3cIAIDCEtUN3tIQCKAhI/PIKqvLioqK4OTkhFOnTqFVq1YAgEaNGqFTp041qlfhbrrWrVvj0qVLFcp//vlndOzYsUbBEEJIQyWU00VX354ZaWlp1crgNYVbRsuWLYOvry/S0tIgkUhw9OhRJCQkYM+ePTh16pTKAySEkIZAKCjdZJXXNzNnzsTatWuxc+dOaGoqNZFPBQrXMnToUJw8eRLLly+Hnp4eli1bhk6dOuHkyZPo27evSoIihJCGpiG99Hr9+nVERkbi3LlzcHZ2hp6edLf30aNHFa5TqZT23//+FxEREcqcSgghH6SGlIyMjY3h7e2t0jpV074ihBBSKSFkP6RX3bzX6hMWFqbyOquVjExMTCCo5oiPzMzMGgVECCENkYZATsuono2mqy3VSkYbN26s5TAIIaRh0xCWbrLKSTWTka+vb23HQQghDZpAIHtSVGoYlapWMsrJyal2hfVlpUJCCFEnahlVrlrJyNjYuMpnRowxWs+IEELkaCgzMADA48eP0bRpU5XWWa1kFBUVVa3K7t27V6NgCCGkoaqtl17XrFmDgIAAzJkzh3u+X1BQgAULFuDAgQMQi8Xw8vLCTz/9BAsLi5pd7B/NmzdHr169MGnSJIwcORI6Ojo1rrNayahXr15y9+Xm5mL//v3YuXMnbt68CT8/vxoHRQghDU1tvGd0/fp1hISEoF27dlLl8+bNw+nTp3H48GEYGRnBz88PI0aMwO+//670td4XHx+PsLAwzJ8/H35+fvj0008xadIkdO3aVek6le6tjImJga+vL6ysrPDdd9+hT58+uHr1qtKBEEJIQ6bqJSTy8vLg4+ODHTt2wMTEhCvPzs5GaGgoNmzYgD59+qBz584ICwvDlStXVPYd3aFDB/z44494/vw5du3ahRcvXqBHjx5o27YtNmzYwC22pwiFklF6ejrWrFmDFi1a4JNPPoGhoSHEYjGOHz+ONWvWoEuXLgoHkJ6ejjlz5qB58+bQ0dGBhYUF3N3dsW3bNrx9+1bh+gghpC4q66aTtSlj5syZGDRoEDw9PaXKb968iaKiIqlyJycn2NraIjY2tia3UIGmpiZGjBiBw4cPY+3atUhKSoK/vz9sbGwwbtw4bs2j6qh2Mho8eDAcHR1x9+5dbNy4Ec+fP8fmzZuVuoEyjx8/RseOHXHu3DmsWrUKt27dQmxsLBYuXIhTp07h/PnzNaqfEELqCuE/AxjKb2XDvXNycqQ2sVgst64DBw4gPj4eq1evrrAvPT0d2traMDY2liq3sLBAenq6Su/pxo0bmDFjBqysrLBhwwb4+/vj0aNHiIiIwPPnzzF06NBq11Xt6YDOnDmD2bNnY/r06WjRooVSgZc3Y8YMaGpq4saNG1IT7TVt2hRDhw4FY0wl1yGEEL5VNbTbxsZGqjwwMBBBQUEVjk9NTcWcOXMQERGhkoEDytiwYQPCwsKQkJCAgQMHYs+ePRg4cCCEwtKbcXBwQHh4OOzt7atdZ7WT0eXLlxEaGorOnTujVatWGDt2LEaPHq3wTZT5+++/uRZR+Rlfy1R3CiJCCKnrhO+1gsqXA6VJ5v33NEUikcx6bt68iZcvX0otZldSUoKYmBhs2bIFv/32GwoLC5GVlSXVOsrIyIClpaVK7mXbtm2YOHEixo8fDysrK5nHmJubIzQ0tNp1Vrubrlu3btixYwdevHiBqVOn4sCBA7C2toZEIkFERARyc3OrfVEASEpKAmMMjo6OUuWNGzeGvr4+9PX1sWjRIpnnisXiCk1aQgipy8rmpquw/ZOMDA0NpTZ5yejjjz/GvXv3cPv2bW5zcXGBj48P9+9aWlqIjIzkzklISEBKSgrc3NxUci+JiYkICAiQm4gAQFtbW6HZexSetVtPTw8TJ07ExIkTkZCQgNDQUKxZswaLFy9G37598csvvyhapZRr165BIpHAx8dHbp/p6tWr8c0339ToOoQQok4agtJNVrkiDAwM0LZtW6kyPT09mJmZceWTJk3C/PnzYWpqCkNDQ8yaNQtubm7o1q2bsuHL9PbtW6SkpKCwsFCqvPxQ8+qo0UQUjo6OWLduHZ49e4b9+/crdG7z5s0hEAiQkJAgVd60aVM0b94curq6cs8NCAhAdnY2t6WmpioVPyGEqItAIJC7qdoPP/yA//u//4O3tzd69uwJS0tLpRa8k+fVq1cYNGgQDAwM0KZNG3Ts2FFqU4ZKZkXS0NDAsGHDFGoVmZmZoW/fvtiyZQvy8/MVup5IJKrQpCWEkLpM1kg6eVMEKSo6OlpqdQUdHR1s3boVmZmZyM/Px9GjR1X2vAgA5s6di+zsbMTFxUFXVxdnz57F7t270aJFC6V7x3hdXO+nn36Cu7s7XFxcEBQUhHbt2kEoFOL69et48OABOnfuzGd4hBCiMgI57xTVx3FaFy5cwIkTJ+Di4gKhUAg7Ozv07dsXhoaGWL16NQYNGqRwnbwmo2bNmuHWrVtYtWoVAgIC8OzZM4hEIrRu3Rr+/v6YMWMGn+ERQojKyOuSq4+jhvPz82Fubg6gdPHVV69eoWXLlnB2dkZ8fLxSdfK+7LiVlRU2b95c4xdoCSGkLmtIs3Y7OjoiISEB9vb2aN++PUJCQmBvb4/g4OBKR9hVhvdkRAghH4Kq3jOqT+bMmcNN9RMYGIj+/ftj79690NbWRnh4uFJ1UjIihBA1EEL2iLH6uLbemDFjuH/v3LkzkpOT8eDBA9ja2qJx48ZK1UnJiBBC1KAhtYzKa9SokdSMEMqoj0mZEELqndoc2q1u3t7eWLt2bYXydevW4ZNPPlGqTkpGhBCiBgKB/K2+iYmJwcCBAyuUDxgwADExMUrVSd10hBCiBg2pmy4vLw/a2toVyrW0tJSeK5RaRoQQogYCCCCUsQlQ/5KRs7MzDh48WKH8wIEDaN26tVJ1UsuIEELUQCgs3WSV1zdLly7FiBEj8OjRI/Tp0wcAEBkZif379+Pw4cNK1UnJiBBC1KAhvfQ6ePBgHD9+HKtWrcLPP/8MXV1dtGvXDufPn0evXr2UqpOSESGEqIFATpdcfeymA4BBgwYpNQedPPWwgUgIIfWPUPDvIAbpje/IlJOVlYWdO3diyZIlyMzMBADEx8cjLS1NqfqoZUQIIWrQkLrp7t69C09PTxgZGeHp06f44osvYGpqiqNHjyIlJQV79uxRuE5qGRFCiBo0pPeM5s+fj/HjxyMxMRE6Ojpc+cCBA+k9I0IIqcs0IKdlVA+fGV2/fh0hISEVyps0aYL09HSl6mxQyWiIoyX0Dfhf9fWzdh/xHQLnkJNykxbWhtEHHvMdAueXOf/lOwQAshdb48u1ZZ58h8DJExfzHQIAoKhYorK65LWC6mPLSCQSyXy59eHDh/jPf/6jVJ3UTUcIIWogqOSf+mbIkCFYvnw5ioqKAJQuEJiSkoJFixbB29tbqTopGRFCiBqULTtefquPLaPvv/8eeXl5MDc3x7t379CrVy80b94cBgYGWLlypVJ1UjIihBA1KFt2XNamiNWrV6NLly4wMDCAubk5hg0bhoSEBKljCgoKMHPmTJiZmUFfXx/e3t7IyMhQ2b0YGRkhIiICJ0+exKZNm+Dn54dff/0VFy9ehJ6enlJ1NqhnRoQQUlep6pnRxYsXMXPmTHTp0gXFxcVYsmQJ+vXrh7/++otLBPPmzcPp06dx+PBhGBkZwc/PDyNGjMDvv/+ugjv5V48ePdCjRw+V1EXJiBBC1EBVyejs2bNSP4eHh8Pc3Bw3b95Ez549kZ2djdDQUOzbt4+bNy4sLAytWrXC1atX0a1bN2VvAQAgkUgQHh6Oo0eP4unTpxAIBHBwcMDIkSMxduxYhVt6ZaibjhBC1ED27Av/LiuRk5MjtYnF4mrVm52dDQAwNTUFANy8eRNFRUXw9Px3dKSTkxNsbW0RGxtbo3tgjGHIkCH44osvkJaWBmdnZ7Rp0wbJyckYP348hg8frnTd1DIihBA1EPyzySoHABsbG6nywMBABAUFVVqnRCLB3Llz4e7ujrZt2wIA0tPToa2tDWNjY6ljLSwslH4HqEx4eDhiYmIQGRmJ3r17S+27cOEChg0bhj179mDcuHEK103JiBBC1EDeYIWystTUVBga/vuepEgkqrLOmTNn4o8//sDly5dVF2gl9u/fjyVLllRIRADQp08fLF68GHv37lUqGVE3HSGEqIGsYd1lGwAYGhpKbVUlIz8/P5w6dQpRUVH46KN/X7S3tLREYWEhsrKypI7PyMiApaVlje7h7t276N+/v9z9AwYMwJ07d5Sqm5IRIYSog6CSTQGMMfj5+eHYsWO4cOECHBwcpPZ37twZWlpaiIyM5MoSEhKQkpICNze3Gt1CZmYmLCws5O63sLDAmzdvlKqbuukIIUQN3h+sUL5cETNnzsS+fftw4sQJGBgYcM+BjIyMoKurCyMjI0yaNAnz58+HqakpDA0NMWvWLLi5udV4JF1JSQk0NeWnDQ0NDRQXKzeVEyUjQghRA1UN7d62bRsAwMPDQ6o8LCwM48ePBwD88MMPEAqF8Pb2hlgshpeXF3766SclopbGGMP48ePldiFWdwSgLLwmIw8PD3To0AEbN26UKg8PD8fcuXMr9HkSQkh9paqVXhljVR6jo6ODrVu3YuvWrQrVXRVfX98qj1Fm8AJALSNCCFGL9wcrlC+vL8LCwmqtbkpGhBCiBlUN7f7QUTIihBB1kDdDN+UiAJSMCCFELRrS4nq1oV4mI7FYLDVqQ9aKg4QQUpeoamh3Q8XrS6+GhobcJH/vy8rKgpGRkdzzVq9eDSMjI24rP6cTIYTUNSp657XB4jUZOTo6Ij4+vkJ5fHw8WrZsKfe8gIAAZGdnc1tqampthkkIITWmqsX1Gipeu+mmT5+OLVu2YPbs2fjiiy8gEolw+vRp7N+/HydPnpR7nkgkqtYkgoQQUlcIIWdot9ojqZt4TUZNmzZFTEwMvvrqK3h6eqKwsBBOTk44fPhwpZPxEUJIfUMDGCrH+wCGLl264Ny5c3yHQQghtYreM6oc78mIEEI+BA1hBobaRMmIEELUQFVz0zVUlIwIIUQN6JlR5SgZEUKIGtBLr5WjZEQIIeog7w1XykUAKBkRQoha0ACGylEyIoQQNRDI6aajod2l6OVfQgghvKOWESGEqAENYKgcJSNCCFEDemZUOeqmI4QQdVDxGhJbt26Fvb09dHR04OrqimvXrqkqUl5QMiKEEDUo66aTtSnq4MGDmD9/PgIDAxEfH4/27dvDy8sLL1++rIXI1YOSESGEqEHZDAyyNkVt2LABkydPxoQJE9C6dWsEBwejUaNG2LVrl+oDV5MG8cyIMQYAyMvL5TmSUoVaGnyHwCnIrxufCQCUiPP5DoGTm1s3lqqvS88LdFB31gjLFxfzHQIAIDe39O+n7DtGGdz3U26uzHno8v65RvlVr+Wt21ZYWIibN28iICCAKxMKhfD09ERsbKzScfKtQSSjsv9g+rg48hwJqS86fMd3BKQ+yc3NhZGRkdLnAkALB5tKj7O1tZX6OTAwEEFBQRWOe/36NUpKSmBhYSFVbmFhgQcPHigVY13QIJKRtbU1UlNTYWBgoPQLZDk5ObCxsUFqaioMDQ1VHCHFQrFQLPUxFsYYcnNzYW1trXQcVX0/Mcbw999/w9TUFELhv09OPrTVrBtEMhIKhfjoo49UUpehoSHvf0RlKBbZKBbZKBbZahqLsi2iMtX5flLkGo0bN4aGhgYyMjKkyjMyMmBpaalUjHUBDWAghJB6RFtbG507d0ZkZCRXJpFIEBkZCTc3Nx4jq5kG0TIihJAPyfz58+Hr6wsXFxd07doVGzduRH5+PiZMmMB3aEqjZPQPkUiEwMDAOtFPS7FQLBRLw4pF1T799FO8evUKy5YtQ3p6Ojp06ICzZ89WGNRQnwhYTcYsEkIIISpAz4wIIYTwjpIRIYQQ3lEyIoQQwjtKRoQQQnhHyegfQUFBcHJygp6eHkxMTODp6Ym4uDi1x1FUVIRFixbB2dkZenp6sLa2xrhx4/D8+XO1xwIAR48eRb9+/WBmZgaBQIDbt2/zEsf7nj59ikmTJsHBwQG6urpo1qwZAgMDUVhYyEs8Q4YMga2tLXR0dGBlZYWxY8fy9vsqIxaL0aFDB15/Z/b29hAIBFLbmjVreIkFAE6fPg1XV1fo6urCxMQEw4YN4y0WUhElo3+0bNkSW7Zswb1793D58mXY29ujX79+ePXqlVrjePv2LeLj47F06VLEx8fj6NGjSEhIwJAhQ9QaR5n8/Hz06NEDa9euVfu1PTw8EB4eXqH8wYMHkEgkCAkJwZ9//okffvgBwcHBWLJkidpjAYDevXvj0KFDSEhIwJEjR/Do0SOMHDmSl1jKLFy4sEZT2KgqluXLl+PFixfcNmvWLF5iOXLkCMaOHYsJEybgzp07+P333/H555/XWixECYzIlJ2dzQCw8+fP8x0Ku3btGgPAkpOTeYvhyZMnDAC7deuW2q7Zq1cvFhYWVq1j161bxxwcHOpELCdOnGACgYAVFhbyEsuvv/7KnJyc2J9//lnrv7PKYrGzs2M//PBDrV27urEUFRWxJk2asJ07d6otFqI4ahnJUFhYiO3bt8PIyAjt27fnOxxkZ2dDIBDA2NiY71DqrOzsbJiamvIdBjIzM7F37150794dWlpaar9+RkYGJk+ejP/9739o1KiR2q9f3po1a2BmZoaOHTti/fr1KC5W/9IQ8fHxSEtLg1AoRMeOHWFlZYUBAwbgjz/+UHssRD5KRu85deoU9PX1oaOjgx9++AERERFo3LgxrzEVFBRg0aJF+Oyzz+rMxJN1TVJSEjZv3oypU6fyFsOiRYugp6cHMzMzpKSk4MSJE2qPgTGG8ePHY9q0aXBxcVH79cubPXs2Dhw4gKioKEydOhWrVq3CwoUL1R7H48ePAZQ+F/76669x6tQpmJiYwMPDA5mZmWqPh8jBd9OMD//v//0/pqenx20xMTGMMcby8vJYYmIii42NZRMnTmT29vYsIyODl1gYY6ywsJANHjyYdezYkWVnZ9dqHFXFoo5uupUrV0pdXygUMpFIJFVWvqvy2bNnrFmzZmzSpEm8xvLq1SuWkJDAzp07x9zd3dnAgQOZRCJRayw//vgjc3d3Z8XFxYyx2vmdKfM7KhMaGso0NTVZQUGBWmPZu3cvA8BCQkK4cwsKCljjxo1ZcHCwSmIhNfdBTgeUm5srNf16kyZNoKurW+G4Fi1aYOLEiVIrKqorlqKiIowaNQqPHz/GhQsXYGZmVmsxVBULUDqCzcHBAbdu3UKHDh1q5fqZmZlS/6fq4+MDb29vjBgxgiuzt7eHpmbplIrPnz+Hh4cHunXrhvDwcKm1YNQdy/uePXsGGxsbXLlyRSWzKFc3lpEjR+LkyZNSa+aUlJRAQ0MDPj4+2L17t9pikfW5/Pnnn2jbti0ePHgAR8eaL4RZ3VguXbqEPn364NKlS+jRowe3z9XVFZ6enli5cmWNYyE190FOlGpgYAADA4Mqj5NIJBCLxWqPpSwRJSYmIioqSi2JSF4s6mRqair13EdXVxfm5uZo3rx5hWPT0tLQu3dvdO7cGWFhYSpNRIrGUp5EIgEAlf23U91YNm3ahBUrVnA/P3/+HF5eXjh48CBcXV3VGosst2/fhlAohLm5uVpj6dy5M0QiERISErhkVFRUhKdPn8LOzk4lsZCa+yCTUXn5+flYuXIlhgwZAisrK7x+/Rpbt25FWloaPvnkE7XGUlRUhJEjRyI+Ph6nTp1CSUkJ0tPTAZT+8Wlra6s1nszMTKSkpHDvzSQkJAAALC0teVvIKy0tDR4eHrCzs8N3330nNfxe3THFxcXh+vXr6NGjB0xMTPDo0SMsXboUzZo1U/vaMuWXrdbX1wcANGvWTGWLT1ZXbGws4uLi0Lt3bxgYGCA2Nhbz5s3DmDFjYGJiotZYDA0NMW3aNAQGBsLGxgZ2dnZYv349AKj975vIR8kIgIaGBh48eIDdu3fj9evXMDMzQ5cuXXDp0iW0adNGrbGkpaXhl19+AYAK3WFRUVHw8PBQazy//PKL1Bopo0ePBgAEBgYiKChIrbGUiYiIQFJSEpKSkip8yaq717lRo0Y4evQoAgMDkZ+fDysrK/Tv3x9ff/11g1y6oLpEIhEOHDiAoKAgiMViODg4YN68eZg/fz4v8axfvx6ampoYO3Ys3r17B1dXV1y4cEHtiZHI90E+MyKEEFK30NBuQgghvKNkRAghhHeUjAghhPCOkhEhhBDeUTIihBDCO0pGhBBCeEfJiBBCCO8oGX3ACgoKsHLlSiQmJvIdCiHkA0fJ6AM2e/ZsJCYmokWLFrV2jfHjx0st7+zh4YG5c+fW2vUqEx0dDYFAgKysLF6uXxNBQUFVTlD79OnTOrM0PCGKomTUwAQFBUEgEEhtTk5OFY7bu3cvUlJSsHPnTrXGd/ToUXz77bdqvWaZ7t2748WLFzAyMgIAhIeH15sFC/39/REZGcn9XD7JA4CNjQ1evHiBtm3bqjk6QmqO5qZrgNq0aYPz589zP8uazt/Hxwc+Pj7qDAsAeF2NVVtbm7fJXWtKX1+fm/hUHg0NjXp7f4RQy6gB0tTU5GbVtrS0rNFqtdHR0ejatSv09PRgbGwMd3d3JCcnA/i36ygkJAQ2NjZo1KgRRo0ahezsbLn1le+mE4vFWLRoEWxsbCASidC8eXOEhoZy+//44w8MGDAA+vr6sLCwwNixY/H69Wu59ScnJ2Pw4MEwMTGBnp4e2rRpg19//ZW7l7JuuujoaEyYMIFb0l0gEHATv4rFYvj7+6NJkybQ09ODq6sroqOjK/2cBAIBtm3bhgEDBkBXVxdNmzbFzz//LHXMvXv30KdPH+jq6sLMzAxTpkxBXl6eQp912b/v3r0bJ06c4GKPjo6W2U138eJFdO3aFSKRCFZWVli8eLHU0t8eHh6YPXs2Fi5cCFNTU1haWvI2AS75sFEyaoASExNhbW2Npk2bwsfHBykpKUrVU1xcjGHDhqFXr164e/cuYmNjMWXKFKnF25KSknDo0CGcPHkSZ8+exa1btzBjxoxqX2PcuHHYv38/Nm3ahPv37yMkJIRrAWRlZaFPnz7o2LEjbty4gbNnzyIjIwOjRo2SW9/MmTMhFosRExODe/fuYe3atTJbFN27d8fGjRthaGiIFy9e4MWLF/D39wcA+Pn5ITY2FgcOHMDdu3fxySefoH///lUO9Fi6dCm8vb1x584d+Pj4YPTo0bh//z6A0mVKvLy8YGJiguvXr+Pw4cM4f/48/Pz8qv1Zl/H398eoUaPQv39/Lvbu3btXOC4tLQ0DBw5Ely5dcOfOHWzbtg2hoaFSax4BwO7du6Gnp4e4uDisW7cOy5cvR0RERKX3SojK8brOLFG5X3/9lR06dIjduXOHnT17lrm5uTFbW1uWk5OjcF1///03A8Cio6Nl7g8MDGQaGhrs2bNnXNmZM2eYUChkL168YIwx5uvry4YOHcrt79WrF5szZw5jjLGEhAQGgEVERMis/9tvv2X9+vWTKktNTWUAWEJCgsxznJ2dWVBQkMx9UVFRDAB78+YNY4yxsLAwZmRkJHVMcnIy09DQYGlpaVLlH3/8MQsICJBZL2OMAWDTpk2TKnN1dWXTp09njDG2fft2ZmJiwvLy8rj9p0+fZkKhkKWnp1frs27fvj33c/nPlbGKy4wvWbKEOTo6Si1/vnXrVqavr89KSkoYY6W/jx49ekjV06VLF7Zo0SK590pIbaBnRg3MgAEDuH9v164dXF1dYWdnh0OHDmHSpEkK1WVqaorx48fDy8sLffv2haenJ0aNGgUrKyvuGFtbWzRp0oT72c3NDRKJBAkJCVU+v7h9+zY0NDTQq1cvmfvv3LmDqKgomS2bR48eoWXLlhXKZ8+ejenTp+PcuXPw9PSEt7c32rVrV91bxr1791BSUlKhbrFYXOWKu+UX03Nzc+O6zO7fv4/27dtDT0+P2+/u7s59Vj179qzys1bU/fv34ebmJtW6cnd3R15eHp49e8Ytxlf+87GyssLLly+Vvi4hyqBuugbO2NgYLVu2RFJSklLnh4WFITY2Ft27d8fBgwfRsmVLXL16VSWx6erqVro/Ly8PgwcPxu3bt6W2xMRE9OzZU+Y5X3zxBR4/foyxY8fi3r17cHFxwebNm6sdU15eHjQ0NHDz5k2pa96/fx8//vijQvenqNr8rCujpaUl9bNAIOCWTidEXSgZNXB5eXl49OhRjf4Pu2PHjggICMCVK1fQtm1b7Nu3j9v3/pLkAHD16lUIhUI4OjpWWa+zszMkEgkuXrwoc3+nTp3w559/wt7eHs2bN5fa3m9hlGdjY4Np06bh6NGjWLBgAXbs2CHzOG1tbZSUlFS415KSErx8+bLCNatq6ZVPHFevXkWrVq0AAK1atcKdO3eQn5/P7f/9998rfFaVfdZVxV5eq1atEBsbK7X67e+//w4DAwO1L0NOSFUoGTUw/v7+uHjxIp4+fYorV65g+PDh0NDQwGeffaZwXU+ePEFAQABiY2ORnJyMc+fOITExkfuCBQAdHR34+vrizp07uHTpEmbPno1Ro0ZVa4ixvb09fH19MXHiRBw/fhxPnjxBdHQ0Dh06BKB0MEJmZiY+++wzXL9+HY8ePcJvv/2GCRMmyP0injt3Ln777Tc8efIE8fHxiIqKkoq3/PXz8vIQGRmJ169f4+3bt2jZsiV8fHwwbtw4HD16FE+ePMG1a9ewevVqnD59utL7OXz4MHbt2oWHDx8iMDAQ165d4wYo+Pj4cJ/VH3/8gaioKMyaNQtjx46FhYVFtT7r8rHfvXsXCQkJeP36NYqKiiocM2PGDKSmpmLWrFl48OABTpw4gcDAQMyfPx9CIf3pkzqG74dWRLU+/fRTZmVlxbS1tVmTJk3Yp59+ypKSkpSqKz09nQ0bNoyrz87Oji1btox7+F32UP2nn35i1tbWTEdHh40cOZJlZmZydVQ2gIExxt69e8fmzZvHXaN58+Zs165d3P6HDx+y4cOHM2NjY6arq8ucnJzY3LlzpR7Kv8/Pz481a9aMiUQi9p///IeNHTuWvX79mjFWcQADY4xNmzaNmZmZMQAsMDCQMcZYYWEhW7ZsGbO3t2daWlrMysqKDR8+nN29e1fuZwWAbd26lfXt25eJRCJmb2/PDh48KHXM3bt3We/evZmOjg4zNTVlkydPZrm5uQp91mVevnzJ+vbty/T19RkAFhUVVWEAA2OMRUdHsy5dujBtbW1maWnJFi1axIqKiuT+PhhjbOjQoczX11fuvRJSGwSMvdeGJ0QBQUFBOH78OE0/g9LnLMeOHaswKwIhpHqorU4IIYR3lIwIIYTwjrrpCCGE8I5aRoQQQnhHyYgQQgjvKBkRQgjhHSUjQgghvKNkRAghhHeUjAghhPCOkhEhhBDeUTIihBDCO0pGhBBCePf/ASSqcLDQBWEpAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 450x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 1, figsize=(4.5, 2.))\n",
    "axes.set_facecolor('lightgrey')\n",
    "sns.heatmap(delta.T * 100, cmap='Blues', ax=axes,\n",
    "            cbar_kws={'label': 'Decay rate (%)'},\n",
    "            vmin=0, vmax=100)\n",
    "axes.set(xlabel='5´ splice site position', \n",
    "         ylabel='Allele')\n",
    "axes.set_yticklabels(alleles, rotation=0)\n",
    "sns.despine(left=False, bottom=False, top=False, right=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7a19307",
   "metadata": {},
   "source": [
    "The decay rate of a mutation is the percent decrease in the predictability of other mutatiob. These effects are constant and thus accumulate as discount coupons when two sequences differ by sevelal mutations. In the Jenga model, a mutation decreases the predictability of other mutations by accumulating the decay rates associated to the alleles that are exchanged. Next, we show examples of mutations, even at the same position, that have very different effects on the predictability of other mutations under this prior. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5537c593",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For instance, a G-1U mutation is expected to decrease the predictability of other mutations by 89.77%\n"
     ]
    }
   ],
   "source": [
    "d = (1 - (1 - delta.loc['-1', 'G']) * (1 - delta.loc['-1', 'U'])) * 100\n",
    "print('For instance, a G-1U mutation is expected to decrease the predictability of other mutations by {:.2f}%'.format(d))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9253948b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In contrast, a mutation at the same position, A-1C, is expected to decrease the predictability of other mutations by only 15.20%\n"
     ]
    }
   ],
   "source": [
    "d = (1 - (1 - delta.loc['-1', 'A']) * (1 - delta.loc['-1', 'C'])) * 100\n",
    "print('In contrast, a mutation at the same position, A-1C, is expected to decrease the predictability of other mutations by only {:.2f}%'.format(d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "176ae219",
   "metadata": {},
   "source": [
    "In general, we see that the positions close to the exon-intron boundary, which are known to be more relevant for the recognition of the 5´splice site during the splicing reaction, have stronger decay factors. But we also see that, within each position, alleles with weak decay factors are more likely to be generally exchangeable without affecting the effects of other mutations. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}