{
"cells": [
{
"cell_type": "markdown",
"id": "ac19c62d-0304-40ef-9ce2-32f4f340488a",
"metadata": {},
"source": [
"# Sensitivity Analysis\n",
"\n",
"The following tutorial demonstrates how one may perform a sensitivity analysis of a simulation model via `calisim`. We will first import our required dependencies."
]
},
{
"cell_type": "code",
"execution_count": 99,
"id": "3a462176-a0ab-4cd7-93c4-7653db5b3292",
"metadata": {},
"outputs": [],
"source": [
"from calisim.data_model import (\n",
"\tDistributionModel,\n",
"\tParameterDataType,\n",
"\tParameterSpecification,\n",
")\n",
"from calisim.example_models import SirOdesModel\n",
"from calisim.sensitivity import (\n",
"\tSensitivityAnalysisMethod,\n",
"\tSensitivityAnalysisMethodModel,\n",
")\n",
"import numpy as np\n",
"import pandas as pd\n",
"from scipy.integrate import solve_ivp"
]
},
{
"cell_type": "markdown",
"id": "1e8112ea-b09f-4cc9-8ad7-94f8cc927101",
"metadata": {},
"source": [
"## SIR Model Parameters and Initial Conditions\n",
"\n",
"We next define our forward model. We will use an SIR (susceptible, infected, and recovered) compartmental model, combined with SciPy's solver for ordinary differential equations. The SIR model is expressed as a system of ordinary differential equations where:\n",
"\n",
"| Parameter | Value | Description |\n",
"|-----------|-------|-------------|\n",
"| β (beta) | 0.4 | Infection rate: probability of transmission per contact per time unit |\n",
"| γ (gamma) | 0.1 | Recovery rate: fraction of infected recovering per time unit |\n",
"| | | **Average infectious period = 1 / γ = 10 time units** |\n",
"\n",
"With the following compartments:\n",
"\n",
"| Compartment | Symbol | Initial Value | Description |\n",
"|------------|--------|---------------|-------------|\n",
"| Susceptible | S0 | 999 | Individuals who can catch the disease (N - I0 - R0) |\n",
"| Infected | I0 | 1.0 | Individuals currently infected and can spread the disease |\n",
"| Recovered | R0 | 0 | Individuals recovered or removed; no longer infectious |"
]
},
{
"cell_type": "code",
"execution_count": 100,
"id": "adb3cf44-f6dd-4ddb-8152-f710902a07a7",
"metadata": {},
"outputs": [],
"source": [
"def sir_simulate(parameters: dict) -> np.ndarray | pd.DataFrame:\n",
" def dX_dt(_: np.ndarray, X: np.ndarray) -> np.ndarray:\n",
" S, I, _ = X\n",
" dotS = -parameters[\"beta\"] * S * I / parameters[\"N\"]\n",
" dotI = (\n",
" parameters[\"beta\"] * S * I / parameters[\"N\"] - parameters[\"gamma\"] * I\n",
" )\n",
" dotR = parameters[\"gamma\"] * I\n",
" return np.array([dotS, dotI, dotR])\n",
"\n",
" X0 = [parameters[\"S0\"], parameters[\"I0\"], parameters[\"R0\"]]\n",
" t = (parameters[\"t\"].min(), parameters[\"t\"].max())\n",
" x_y = solve_ivp(\n",
" fun=dX_dt, y0=X0, t_span=t, t_eval=parameters[\"t\"].values.flatten()\n",
" ).y\n",
"\n",
" df = pd.DataFrame(dict(dotS=x_y[0, :], dotI=x_y[1, :], dotR=x_y[2, :]))\n",
" return df"
]
},
{
"cell_type": "markdown",
"id": "93a95b93-b79e-45b1-aa0e-760c380b686e",
"metadata": {},
"source": [
"We will perform a simulation study with the following ground-truth parameters:"
]
},
{
"cell_type": "code",
"execution_count": 101,
"id": "95c54892-4ce1-41bc-85c0-62f910ccabdf",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" beta | \n",
" gamma | \n",
" N | \n",
" I0 | \n",
" R0 | \n",
" S0 | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 0.4 | \n",
" 0.1 | \n",
" 1000 | \n",
" 1.0 | \n",
" 0 | \n",
" 999.0 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" beta gamma N I0 R0 S0\n",
"0 0.4 0.1 1000 1.0 0 999.0"
]
},
"execution_count": 101,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model = SirOdesModel()\n",
"pd.DataFrame(model.GROUND_TRUTH, index=[0])"
]
},
{
"cell_type": "markdown",
"id": "8b19c4fd-760a-4350-998b-1a1df38364f0",
"metadata": {},
"source": [
"When supplied to our forward model, these ground-truth parameters will generate the observed data below:"
]
},
{
"cell_type": "code",
"execution_count": 102,
"id": "97550d3b-4961-466d-aff2-f59ebd81c8af",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" dotS | \n",
" dotI | \n",
" dotR | \n",
" day | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 999.000000 | \n",
" 1.000000 | \n",
" 0.000000 | \n",
" 0 | \n",
"
\n",
" \n",
" | 1 | \n",
" 998.534208 | \n",
" 1.349201 | \n",
" 0.116592 | \n",
" 1 | \n",
"
\n",
" \n",
" | 2 | \n",
" 997.906105 | \n",
" 1.819995 | \n",
" 0.273899 | \n",
" 2 | \n",
"
\n",
" \n",
" | 3 | \n",
" 997.059813 | \n",
" 2.454180 | \n",
" 0.486007 | \n",
" 3 | \n",
"
\n",
" \n",
" | 4 | \n",
" 995.919926 | \n",
" 3.308098 | \n",
" 0.771976 | \n",
" 4 | \n",
"
\n",
" \n",
" | 5 | \n",
" 994.385263 | \n",
" 4.457212 | \n",
" 1.157524 | \n",
" 5 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" dotS dotI dotR day\n",
"0 999.000000 1.000000 0.000000 0\n",
"1 998.534208 1.349201 0.116592 1\n",
"2 997.906105 1.819995 0.273899 2\n",
"3 997.059813 2.454180 0.486007 3\n",
"4 995.919926 3.308098 0.771976 4\n",
"5 994.385263 4.457212 1.157524 5"
]
},
"execution_count": 102,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"observed_data = model.get_observed_data()\n",
"observed_data.head(6)"
]
},
{
"cell_type": "markdown",
"id": "c639dcdd-3dcd-47df-88ea-52033c479677",
"metadata": {},
"source": [
"Let's view the trajectory of infected individuals over time in days."
]
},
{
"cell_type": "code",
"execution_count": 103,
"id": "61b7585f-01c1-4303-9a9c-ad22e3fca433",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 103,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAGwCAYAAABCV9SaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5OUlEQVR4nO3df3RU9Z3/8dcQkmASZjA/h5QAgaRgFJCihCk1RaEERKuF7lF0JWIWDmyASiql6ddqxa6haNWaVdiza0F2JbR0RY+4/kCEZIsRFEhFVJpAILiQZJI0GZJIAuF+//AwdRBCft+Z3OfjnDmH+dw7M+97zz0nLz4/7rUZhmEIAACgj+tndgEAAAC9gdADAAAsgdADAAAsgdADAAAsgdADAAAsgdADAAAsgdADAAAsob/ZBfiD8+fP6+TJkxo4cKBsNpvZ5QAAgHYwDEOnT59WfHy8+vW7cj8OoUfSyZMnlZCQYHYZAACgE06cOKEhQ4ZccT9Cj6SBAwdK+uqk2e12k6sBAADt4fF4lJCQ4P07fiWEHsk7pGW32wk9AAAEmPZOTWEiMwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsAQeQwEAQB931N2g47VNGh4VrsTocLPLMQ2hBwCAPqquqUXL8otVWOL2tqUlxyhv7ng5woJNrMwcDG8BANBHLcsv1u7Sap+23aXVWpp/oEPfc9TdoJ2Hq1RW3did5fU6vwk9q1evls1m04MPPuhtO3PmjLKyshQVFaWIiAjNmTNHlZWVPp8rLy/XrFmzFBYWptjYWK1YsULnzp3r5eoBAPAvR90NKixxq9UwfNpbDUOFJe52BZi6phbNe3Gvbvltgeav/1A3P7VL817cq/qmsz1Vdo/yi9Dz4Ycf6t/+7d80duxYn/bly5fr9ddf15YtW1RQUKCTJ09q9uzZ3u2tra2aNWuWWlpa9P777+ull17Shg0b9Mgjj/T2IQAA0GM609NyvLapze3Haq78Xd3VU+QvTJ/T09DQoHvvvVf//u//rl//+tfe9vr6er344ovatGmTbrnlFknS+vXrdc011+iDDz7QpEmT9M477+jTTz/Vu+++q7i4OF1//fV6/PHHtXLlSv3qV79SSEiIWYcFAECXdWVOzrDIsDa3D49qe0LzhZ6ii329pyjQJkWb3tOTlZWlWbNmadq0aT7t+/bt09mzZ33aR48eraFDh6qoqEiSVFRUpDFjxiguLs67T3p6ujwejw4dOnTZ32xubpbH4/F5AQDgb7rS0zIiJkJpyTEKstl82oNsNqUlx1wxsHRHT5G/MTX0bN68Wfv371dubu43tlVUVCgkJESDBg3yaY+Li1NFRYV3n68HngvbL2y7nNzcXDkcDu8rISGhi0cCAED36o45OXlzx2tyUrRP2+SkaOXNHX/Fz3a1p8gfmTa8deLECf3kJz/R9u3bNWDAgF797ZycHGVnZ3vfezwegg8AwK+0p6flSr01jrBgbcycqLLqRh2raezQfXou9BTtLq32CV5BNpsmJ0W3+3v86R5BpoWeffv2qaqqSt/5zne8ba2trSosLNS//uu/6u2331ZLS4vq6up8ensqKyvldDolSU6nU3v37vX53guruy7scymhoaEKDQ3txqMBAKB7dWdPS2J05wJH3tzxWpp/wGduT3t7ivzxHkGmhZ6pU6fq4MGDPm3z58/X6NGjtXLlSiUkJCg4OFg7duzQnDlzJEmHDx9WeXm5XC6XJMnlculf/uVfVFVVpdjYWEnS9u3bZbfblZKS0rsHBABAN+qunpau6EpPUVvzkTZmTuyJcq/ItNAzcOBAXXfddT5t4eHhioqK8rZnZmYqOztbkZGRstvtWrp0qVwulyZNmiRJmj59ulJSUnTfffdpzZo1qqio0MMPP6ysrCx6cgAAAa8rPS3dqaM9Rf668sv0JetteeaZZ9SvXz/NmTNHzc3NSk9P1wsvvODdHhQUpG3btmnx4sVyuVwKDw9XRkaGVq1aZWLVAAB0j670tJipO+Yj9QSbYVw0LdyCPB6PHA6H6uvrZbfbzS4HAICAdtTdoFt+W3DZ7TsfmtItoaejf79Nv08PAABW0FeeX9UeXb1HUE/x6+EtAAACnT+uYuoN/jIf6esY3hLDWwCAnjPvxb2XXYFl1iqm3tST85E6+vebnh4AAHqIv65i6k2dvUdQT2BODwAAPaQvPr8qkBF6AADoIX3x+VWBjNADAEAP8ddVTFZF6AEAoAd15Unn6F5MZAYAoAcF6l2V+yJCDwAAvcCfVjFZFcNbAADAEgg9AADAEgg9AADAEgg9AADAEgg9AADAEgg9AADAEliyDgBAOxx1N+h4bRP32QlghB4AANpQ19SiZfnFPk9LT0uOUd7c8XKEBZtYGTqK4S0AANqwLL9Yu0urfdp2l1Zraf4BkypCZxF6AAC4jKPuBhWWuNVqGD7trYahwhK3yqobTaoMnUHoAQDgMo7XNrW5/VgNoSeQEHoAALiMYZFhbW4fHsWE5kBC6AEA4DJGxEQoLTlGQTabT3uQzaa05BhWcQUYQg8AAG3Imztek5OifdomJ0Urb+54kypCZ7FkHQCANjjCgrUxc6LKqht1rKaR+/QEMEIPAADtkBhN2Al0DG8BAABLIPQAAABLIPQAAABLIPQAAABLIPQAAABLMDX0rF27VmPHjpXdbpfdbpfL5dKbb77p3T5lyhTZbDaf16JFi3y+o7y8XLNmzVJYWJhiY2O1YsUKnTt3rrcPBQAA+DlTl6wPGTJEq1evVnJysgzD0EsvvaQ77rhDBw4c0LXXXitJWrBggVatWuX9TFjY328J3traqlmzZsnpdOr999/XqVOnNG/ePAUHB+uJJ57o9eMBAAD+y2YYFz061mSRkZF68sknlZmZqSlTpuj666/Xs88+e8l933zzTd122206efKk4uLiJEnr1q3TypUr5Xa7FRIS0q7f9Hg8cjgcqq+vl91u765DAQAAPaijf7/9Zk5Pa2urNm/erMbGRrlcLm/7yy+/rOjoaF133XXKyclRU9Pfn3hbVFSkMWPGeAOPJKWnp8vj8ejQoUOX/a3m5mZ5PB6fFwAA6NtMvyPzwYMH5XK5dObMGUVERGjr1q1KSUmRJN1zzz0aNmyY4uPj9fHHH2vlypU6fPiwXnnlFUlSRUWFT+CR5H1fUVFx2d/Mzc3VY4891kNHBAAA/JHpoWfUqFEqLi5WfX29/vSnPykjI0MFBQVKSUnRwoULvfuNGTNGgwcP1tSpU3XkyBGNHDmy07+Zk5Oj7Oxs73uPx6OEhIQuHQcAAPBvpg9vhYSEKCkpSRMmTFBubq7GjRun3/3ud5fcNzU1VZJUWloqSXI6naqsrPTZ58J7p9N52d8MDQ31rhi78AIAAH2b6aHnYufPn1dzc/MltxUXF0uSBg8eLElyuVw6ePCgqqqqvPts375ddrvdO0QGAAAgmTy8lZOTo5kzZ2ro0KE6ffq0Nm3apF27duntt9/WkSNHtGnTJt16662KiorSxx9/rOXLlystLU1jx46VJE2fPl0pKSm67777tGbNGlVUVOjhhx9WVlaWQkNDzTw0AIAfOupu0PHaJg2P4onpVmRq6KmqqtK8efN06tQpORwOjR07Vm+//bZ+8IMf6MSJE3r33Xf17LPPqrGxUQkJCZozZ44efvhh7+eDgoK0bds2LV68WC6XS+Hh4crIyPC5rw8AAHVNLVqWX6zCEre3LS05Rnlzx8sRFmxiZehNfnefHjNwnx4A6NvmvbhXu0ur1fq1P3lBNpsmJ0VrY+ZEEytDVwTsfXoAAOgJR90NKixx+wQeSWo1DBWWuFVW3WhSZehthB4AQJ92vLapze3Hagg9VkHoAQD0acMiw9rcPjyKCc1WQegBAPRpI2IilJYcoyCbzac9yGZTWnIMq7gshNADAOjz8uaO1+SkaJ+2yUnRyps73qSKYAbTH0MBAEBPc4QFa2PmRJVVN+pYTSP36bEoQg8AwDISowk7VsbwFgAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsIT+ZhcAAEB7HXU36Hhtk4ZHhSsxOtzschBgCD0AAL9X19SiZfnFKixxe9vSkmOUN3e8HGHBJlaGQMLwFgDA7y3LL9bu0mqftt2l1Vqaf8CkihCICD0AAL921N2gwhK3Wg3Dp73VMFRY4lZZdaNJlSHQEHoAAH7teG1Tm9uP1RB60D6EHgCAXxsWGdbm9uFRTGhG+xB6AAB+bURMhNKSYxRks/m0B9lsSkuOYRUX2s3U0LN27VqNHTtWdrtddrtdLpdLb775pnf7mTNnlJWVpaioKEVERGjOnDmqrKz0+Y7y8nLNmjVLYWFhio2N1YoVK3Tu3LnePhQAQA/Kmztek5OifdomJ0Urb+54kypCIDJ1yfqQIUO0evVqJScnyzAMvfTSS7rjjjt04MABXXvttVq+fLneeOMNbdmyRQ6HQ0uWLNHs2bO1e/duSVJra6tmzZolp9Op999/X6dOndK8efMUHBysJ554wsxDAwB0I0dYsDZmTlRZdaOO1TRynx50is0wLpoOb7LIyEg9+eST+vGPf6yYmBht2rRJP/7xjyVJn3/+ua655hoVFRVp0qRJevPNN3Xbbbfp5MmTiouLkyStW7dOK1eulNvtVkhIyCV/o7m5Wc3Nzd73Ho9HCQkJqq+vl91u7/mDBAAAXebxeORwONr999tv5vS0trZq8+bNamxslMvl0r59+3T27FlNmzbNu8/o0aM1dOhQFRUVSZKKioo0ZswYb+CRpPT0dHk8Hh06dOiyv5WbmyuHw+F9JSQk9NyBAQAAv2B66Dl48KAiIiIUGhqqRYsWaevWrUpJSVFFRYVCQkI0aNAgn/3j4uJUUVEhSaqoqPAJPBe2X9h2OTk5Oaqvr/e+Tpw40b0HBQAA/I7pj6EYNWqUiouLVV9frz/96U/KyMhQQUFBj/5maGioQkNDe/Q3AACAfzE99ISEhCgpKUmSNGHCBH344Yf63e9+p7vuukstLS2qq6vz6e2prKyU0+mUJDmdTu3du9fn+y6s7rqwDwAAgOQHw1sXO3/+vJqbmzVhwgQFBwdrx44d3m2HDx9WeXm5XC6XJMnlcungwYOqqqry7rN9+3bZ7XalpKT0eu0AAMB/mdrTk5OTo5kzZ2ro0KE6ffq0Nm3apF27duntt9+Ww+FQZmamsrOzFRkZKbvdrqVLl8rlcmnSpEmSpOnTpyslJUX33Xef1qxZo4qKCj388MPKyspi+AoAAPgwNfRUVVVp3rx5OnXqlBwOh8aOHau3335bP/jBDyRJzzzzjPr166c5c+aoublZ6enpeuGFF7yfDwoK0rZt27R48WK5XC6Fh4crIyNDq1atMuuQAACAn/K7+/SYoaPr/AEAgPkC9j49AAAAPYnQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALIHQAwAALKG/2QUAAKzlqLtBx2ubNDwqXInR4WaXAwsh9AAAekVdU4uW5RersMTtbUtLjlHe3PFyhAWbWBmsguEtAECvWJZfrN2l1T5tu0urtTT/gEkVwWoIPQCAHnfU3aDCErdaDcOnvdUwVFjiVll1o0mVwUoIPQCAHne8tqnN7cdqCD3oeYQeAECPGxYZ1ub24VFMaEbPI/QAAHrciJgIpSXHKMhm82kPstmUlhzDKi70CkIPAKBX5M0dr8lJ0T5tk5OilTd3vEkVwWpYsg4A6BWOsGBtzJyosupGHatp5D496HWEHgBAr0qMJuzAHAxvAQAASyD0AAAASyD0AAAASyD0AAAASyD0AAAASyD0AAAASyD0AAAASyD0AAAASzA19OTm5urGG2/UwIEDFRsbqzvvvFOHDx/22WfKlCmy2Ww+r0WLFvnsU15erlmzZiksLEyxsbFasWKFzp0715uHAgAA/Jypd2QuKChQVlaWbrzxRp07d06/+MUvNH36dH366acKD//73ToXLFigVatWed+Hhf39ab2tra2aNWuWnE6n3n//fZ06dUrz5s1TcHCwnnjiiV49HgAA4L9shmEYZhdxgdvtVmxsrAoKCpSWlibpq56e66+/Xs8+++wlP/Pmm2/qtttu08mTJxUXFydJWrdunVauXCm3262QkJBvfKa5uVnNzc3e9x6PRwkJCaqvr5fdbu/+AwMAAN3O4/HI4XC0+++3X83pqa+vlyRFRkb6tL/88suKjo7Wddddp5ycHDU1NXm3FRUVacyYMd7AI0np6enyeDw6dOjQJX8nNzdXDofD+0pISOiBowEAAP7Ebx44ev78eT344IOaPHmyrrvuOm/7Pffco2HDhik+Pl4ff/yxVq5cqcOHD+uVV16RJFVUVPgEHkne9xUVFZf8rZycHGVnZ3vfX+jpAQAAfZffhJ6srCx98skn+vOf/+zTvnDhQu+/x4wZo8GDB2vq1Kk6cuSIRo4c2anfCg0NVWhoaJfqBQAAgcUvhreWLFmibdu2aefOnRoyZEib+6ampkqSSktLJUlOp1OVlZU++1x473Q6e6BaAAAQiEwNPYZhaMmSJdq6davee+89JSYmXvEzxcXFkqTBgwdLklwulw4ePKiqqirvPtu3b5fdbldKSkqP1A0AAAKPqcNbWVlZ2rRpk1577TUNHDjQOwfH4XDoqquu0pEjR7Rp0ybdeuutioqK0scff6zly5crLS1NY8eOlSRNnz5dKSkpuu+++7RmzRpVVFTo4YcfVlZWFkNYAADAy9Ql6zab7ZLt69ev1/33368TJ07oH//xH/XJJ5+osbFRCQkJ+tGPfqSHH37YZ2na8ePHtXjxYu3atUvh4eHKyMjQ6tWr1b9/+zJdR5e8AQAA83X077df3afHLIQeAAACT0DfpwcAAKCnEHoAAIAlEHoAAIAl+M3NCQEAgeGou0HHa5s0PCpcidHhV/4A4CcIPQCAdqlratGy/GIVlri9bWnJMcqbO16OsGATKwPah+EtAEC7LMsv1u7Sap+23aXVWpp/wKSKgI4h9AAAruiou0GFJW61XnSXk1bDUGGJW2XVjSZVBrQfoQcAcEXHa5va3H6shtAD/0foAQBc0bDIsDa3D49iQjP8H6EHAHBFI2IilJYco6CLHh8UZLMpLTmGVVwICIQeAEC75M0dr8lJ0T5tk5OilTd3vEkVAR3DknUAQLs4woK1MXOiyqobdaymkfv0IOAQegAAHZIYTdhBYGJ4CwAAWAKhBwAAWAKhBwAAWAKhBwAAWAKhBwAAWAKhBwAAWAKhBwAAWAKhBwAAWAKhBwAAWEKH7sicnZ3drv2efvrpThUDAADQUzoUeg4cONBTdQAAAPSoDoWenTt39lQdAAAAPapTc3pWrVqlpqamb7R/+eWXWrVqVZeLAgAA6G42wzCMjn4oKChIp06dUmxsrE97TU2NYmNj1dra2m0F9gaPxyOHw6H6+nrZ7XazywEAAO3Q0b/fnerpMQxDNpvtG+1/+ctfFBkZ2ZmvBAAA6FEdmtNz9dVXy2azyWaz6dvf/rZP8GltbVVDQ4MWLVrU7UUCAAB0VYdCz7PPPivDMPTAAw/osccek8Ph8G4LCQnR8OHD5XK5ur1IAACArupQ6MnIyJAkJSYm6rvf/a6Cg4O79OO5ubl65ZVX9Pnnn+uqq67Sd7/7Xf3mN7/RqFGjvPucOXNGP/3pT7V582Y1NzcrPT1dL7zwguLi4rz7lJeXa/Hixdq5c6ciIiKUkZGh3Nxc9e/focMDAAB9WKdSwfe//321trbqv//7v/XZZ59Jkq699lr98Ic/VFBQULu/p6CgQFlZWbrxxht17tw5/eIXv9D06dP16aefKjw8XJK0fPlyvfHGG9qyZYscDoeWLFmi2bNna/fu3ZK+GlabNWuWnE6n3n//fZ06dUrz5s1TcHCwnnjiic4cHgAA6IM6tXqrtLRUt956q/7v//7P2ytz+PBhJSQk6I033tDIkSM7VYzb7VZsbKwKCgqUlpam+vp6xcTEaNOmTfrxj38sSfr88891zTXXqKioSJMmTdKbb76p2267TSdPnvT2/qxbt04rV66U2+1WSEjIFX+X1VsAAASeXlm9tWzZMo0cOVInTpzQ/v37tX//fpWXlysxMVHLli3rzFdKkurr6yXJuwJs3759Onv2rKZNm+bdZ/To0Ro6dKiKiookSUVFRRozZozPcFd6ero8Ho8OHTp0yd9pbm6Wx+PxeQEAgL6tU8NbBQUF+uCDD3yWp0dFRWn16tWaPHlypwo5f/68HnzwQU2ePFnXXXedJKmiokIhISEaNGiQz75xcXGqqKjw7vP1wHNh+4Vtl5Kbm6vHHnusU3UCAIDA1KmentDQUJ0+ffob7Q0NDe0aTrqUrKwsffLJJ9q8eXOnPt8ROTk5qq+v975OnDjR478JAP7kqLtBOw9Xqay60exSgF7TqZ6e2267TQsXLtSLL76oiRMnSpL27NmjRYsW6Yc//GGHv2/JkiXatm2bCgsLNWTIEG+70+lUS0uL6urqfHp7Kisr5XQ6vfvs3bvX5/sqKyu92y4lNDRUoaGhHa4TAAJdXVOLluUXq7DE7W1LS45R3tzxcoR1bUUu4O861dPz3HPPaeTIkXK5XBowYIAGDBig7373u0pKStKzzz7b7u8xDENLlizR1q1b9d577ykxMdFn+4QJExQcHKwdO3Z42w4fPqzy8nLv/YBcLpcOHjyoqqoq7z7bt2+X3W5XSkpKZw4PAPqsZfnF2l1a7dO2u7RaS/MPmFQR0Hs61dMzaNAgvfbaayotLfUuWb/mmmuUlJTUoe/JysrSpk2b9Nprr2ngwIHeOTgOh0NXXXWVHA6HMjMzlZ2drcjISNntdi1dulQul0uTJk2SJE2fPl0pKSm67777tGbNGlVUVOjhhx9WVlYWvTkA8DVH3Q0+PTwXtBqGCkvcKqtuVGJ0uAmVAb2j3aEnOzu7ze07d+70/vvpp59u13euXbtWkjRlyhSf9vXr1+v++++XJD3zzDPq16+f5syZ43NzwguCgoK0bds2LV68WC6XS+Hh4crIyOBp7wBwkeO1TW1uP1ZD6EHf1u7Qc+CAb9fn/v37de7cOe99ev76178qKChIEyZMaPePt+cWQQMGDNDzzz+v559//rL7DBs2TP/zP//T7t8FACsaFhnW5vbhUQQe9G3tDj0X9+QMHDhQL730kq6++mpJ0t/+9jfNnz9fN910U/dXCQDoshExEUpLjtHu0mq1fu0/nUE2myYnRdPLgz6vU3dk/ta3vqV33nlH1157rU/7J598ounTp+vkyZPdVmBv4I7MAKyivumsluYfYPUW+oSO/v3u1ERmj8cjt/ubk+Hcbvcl798DAPAPjrBgbcycqLLqRh2radTwqHB6eGAZnQo9P/rRjzR//nz99re/9blPz4oVKzR79uxuLRAA0P0Sowk7sJ5OhZ5169bpoYce0j333KOzZ89+9UX9+yszM1NPPvlktxYIAADQHTo1p+eCxsZGHTlyRJI0cuRIhYcH5v8amNMDAEDg6ZU5PReEh4dr7NixXfkKAACAXtGpx1AAAAAEGkIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwhP5mFwAA6Lij7gYdr23S8KhwJUaHm10OEBAIPQAQQOqaWrQsv1iFJW5vW1pyjPLmjpcjLNjEygD/x/AWAASQZfnF2l1a7dO2u7RaS/MPmFQREDgIPQAQII66G1RY4larYfi0txqGCkvcKqtuNKkyIDAQegAgQByvbWpz+7EaQg/QFkIPAASIYZFhbW4fHsWEZqAthB4ACBAjYiKUlhyjIJvNpz3IZlNacgyruIArIPQAQADJmztek5OifdomJ0Urb+54kyoCAoepoaewsFC333674uPjZbPZ9Oqrr/psv//++2Wz2XxeM2bM8NmntrZW9957r+x2uwYNGqTMzEw1NDT04lEAQO9xhAVrY+ZE7XxoitbPv1E7H5qijZkTWa4OtIOpoaexsVHjxo3T888/f9l9ZsyYoVOnTnlf+fn5PtvvvfdeHTp0SNu3b9e2bdtUWFiohQsX9nTpAGCqxOhw3TwqliEtoANMvTnhzJkzNXPmzDb3CQ0NldPpvOS2zz77TG+99ZY+/PBD3XDDDZKkvLw83XrrrXrqqacUHx/f7TUDAIDA5Pdzenbt2qXY2FiNGjVKixcvVk1NjXdbUVGRBg0a5A08kjRt2jT169dPe/bsuex3Njc3y+Px+LwAAEDf5tehZ8aMGdq4caN27Nih3/zmNyooKNDMmTPV2toqSaqoqFBsbKzPZ/r376/IyEhVVFRc9ntzc3PlcDi8r4SEhB49DgAAYD6/fvbW3Xff7f33mDFjNHbsWI0cOVK7du3S1KlTO/29OTk5ys7O9r73eDwEHwAA+ji/7um52IgRIxQdHa3S0lJJktPpVFVVlc8+586dU21t7WXnAUlfzROy2+0+LwAA0LcFVOj54osvVFNTo8GDB0uSXC6X6urqtG/fPu8+7733ns6fP6/U1FSzygQAAH7I1OGthoYGb6+NJJWVlam4uFiRkZGKjIzUY489pjlz5sjpdOrIkSP62c9+pqSkJKWnp0uSrrnmGs2YMUMLFizQunXrdPbsWS1ZskR33303K7cAAIAPm2Fc9LjeXrRr1y7dfPPN32jPyMjQ2rVrdeedd+rAgQOqq6tTfHy8pk+frscff1xxcXHefWtra7VkyRK9/vrr6tevn+bMmaPnnntOERER7a7D4/HI4XCovr6eoS4AAAJER/9+mxp6/AWhBwCAwNPRv98BNacHAACgswg9AADAEgg9AADAEgg9AADAEgg9AADAEgg9AADAEgg9AADAEgg9AADAEvz6KesA0JcddTfoeG2ThkeFKzE63OxygD6P0AMAvayuqUXL8otVWOL2tqUlxyhv7ng5woJNrAzo2xjeAoBetiy/WLtLq33adpdWa2n+AZMqAqyB0AMAveiou0GFJW61XvTYw1bDUGGJW2XVjSZVBvR9hB4A6EXHa5va3H6shtAD9BRCDwD0omGRYW1uHx7FhGagpxB6AKAXjYiJUFpyjIJsNp/2IJtNackxrOICehChBwB6Wd7c8ZqcFO3TNjkpWnlzx5tUEWANLFkHgF7mCAvWxsyJKqtu1LGaRu7TA/QSQg8AmCQxmrAD9CaGtwAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCUQegAAgCXwwFEA6KSj7gYdr23iKelAgDC1p6ewsFC333674uPjZbPZ9Oqrr/psNwxDjzzyiAYPHqyrrrpK06ZNU0lJic8+tbW1uvfee2W32zVo0CBlZmaqoaGhF48CgNXUNbVo3ot7dctvCzR//Ye6+aldmvfiXtU3nTW7NABtMDX0NDY2aty4cXr++ecvuX3NmjV67rnntG7dOu3Zs0fh4eFKT0/XmTNnvPvce++9OnTokLZv365t27apsLBQCxcu7K1DAGBBy/KLtbu02qdtd2m1luYfMKkiAO1hMwzDMLsISbLZbNq6davuvPNOSV/18sTHx+unP/2pHnroIUlSfX294uLitGHDBt1999367LPPlJKSog8//FA33HCDJOmtt97Srbfeqi+++ELx8fGX/K3m5mY1Nzd733s8HiUkJKi+vl52u71nDxRAQDvqbtAtvy247PadD01hqAvoJR6PRw6Ho91/v/12InNZWZkqKio0bdo0b5vD4VBqaqqKiookSUVFRRo0aJA38EjStGnT1K9fP+3Zs+ey352bmyuHw+F9JSQk9NyBAOhTjtc2tbn9WE1jL1UCoKP8NvRUVFRIkuLi4nza4+LivNsqKioUGxvrs71///6KjIz07nMpOTk5qq+v975OnDjRzdUD6KuGRYa1uX14FL08gL+y5Oqt0NBQhYaGml0GgAA0IiZCackx2l1ardavzQ4Istk0OSmaoS3Aj/ltT4/T6ZQkVVZW+rRXVlZ6tzmdTlVVVflsP3funGpra737AEB3y5s7XpOTon3aJidFK2/ueJMqAtAeftvTk5iYKKfTqR07duj666+X9NWEpT179mjx4sWSJJfLpbq6Ou3bt08TJkyQJL333ns6f/68UlNTzSodQB/nCAvWxsyJKqtu1LGaRu7TAwQIU0NPQ0ODSktLve/LyspUXFysyMhIDR06VA8++KB+/etfKzk5WYmJifrlL3+p+Ph47wqva665RjNmzNCCBQu0bt06nT17VkuWLNHdd9992ZVbANBdEqMJO0AgMTX0fPTRR7r55pu977OzsyVJGRkZ2rBhg372s5+psbFRCxcuVF1dnb73ve/prbfe0oABA7yfefnll7VkyRJNnTpV/fr105w5c/Tcc8/1+rEAAAD/5jf36TFTR9f5AwAA8/WZ+/QAAAB0J0IPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBEIPAACwBL99DAUA9LSj7gYdr23iMRKARRB6AFhOXVOLluUXq7DE7W1LS45R3tzxcoQFm1gZgJ7E8BYAy1mWX6zdpdU+bbtLq7U0/4BJFQHoDYQeAJZy1N2gwhK3Wi96Ak+rYaiwxK2y6kaTKgPQ0wg9ACzleG1Tm9uP1RB6gL6K0APAUoZFhrW5fXgUE5qBvorQA8BSRsREKC05RkE2m097kM2mtOQYVnEBfRihB4Dl5M0dr8lJ0T5tk5OilTd3vEkVAegNLFkHYDmOsGBtzJyosupGHatp5D49gEUQegBYVmI0YQewEoa3AACAJRB6AACAJRB6AACAJRB6AACAJRB6AACAJbB6C0BAO+pu0PHaJpadA7giQg+AgFTX1KJl+cUqLHF729KSY5Q3d7wcYcEmVgbAXzG8BSAgLcsv1u7Sap+23aXVWpp/wKSKAPg7Qg+AgHPU3aDCErdaDcOnvdUwVFjiVlk1T0oH8E2EHgAB53htU5vbj9UQegB8E6EHQMAZFhnW5vbhUUxoBvBNhB4AAWdETITSkmMUZLP5tAfZbEpLjmEVF4BLIvQACEh5c8drclK0T9vkpGjlzR1vUkUA/J1fh55f/epXstlsPq/Ro0d7t585c0ZZWVmKiopSRESE5syZo8rKShMrBtBbHGHB2pg5UTsfmqL182/UzoemaGPmRJarA7gsv79Pz7XXXqt3333X+75//7+XvHz5cr3xxhvasmWLHA6HlixZotmzZ2v37t1mlArABInR3JQQQPv4fejp37+/nE7nN9rr6+v14osvatOmTbrlllskSevXr9c111yjDz74QJMmTertUgEAgB/z6+EtSSopKVF8fLxGjBihe++9V+Xl5ZKkffv26ezZs5o2bZp339GjR2vo0KEqKipq8zubm5vl8Xh8XgAAoG/z69CTmpqqDRs26K233tLatWtVVlamm266SadPn1ZFRYVCQkI0aNAgn8/ExcWpoqKize/Nzc2Vw+HwvhISEnrwKAC05ai7QTsPV3FDQQA9zq+Ht2bOnOn999ixY5Wamqphw4bpj3/8o6666qpOf29OTo6ys7O97z0eD8EH6GU8OwtAb/Prnp6LDRo0SN/+9rdVWloqp9OplpYW1dXV+exTWVl5yTlAXxcaGiq73e7zAtC7eHYWgN4WUKGnoaFBR44c0eDBgzVhwgQFBwdrx44d3u2HDx9WeXm5XC6XiVUCuBKenQXADH49vPXQQw/p9ttv17Bhw3Ty5Ek9+uijCgoK0ty5c+VwOJSZmans7GxFRkbKbrdr6dKlcrlcrNwC/Fx7np3FMnQA3c2vQ88XX3yhuXPnqqamRjExMfre976nDz74QDExMZKkZ555Rv369dOcOXPU3Nys9PR0vfDCCyZXDeBKeHYWADPYDOOi/mUL8ng8cjgcqq+vZ34P0EvmvbhXu0urfYa4gmw2TU6K1sbMiSZWBiBQdPTvd0DN6QHQd/DsLAC9za+HtwD0XReenVVW3ahjNY0aHsXjJAD0LEIPgC476m7Q8dqmTgUXnp0FoLcQegB0GjcYBBBImNMDoNO4wSCAQELoAdAp3GAQQKAh9ADolPbcYBAA/AmhB0CncINBAIGG0AOgU0bERCgtOUZBNptPe5DNprTkGFZkAfA7hB4AOupu0M7DVR2eh8MNBgEEEpasAxbW1SXn3GAQQCChpwewsO5acp4YHa6bR8USeAD4NUIPYFEsOQdgNYQewKJYcg7Aagg9gEWx5ByA1RB6gD6ioyuwWHIOwGpYvQUEuK6swMqbO15L8w/4fJYl5wD6KpthXDSL0YI8Ho8cDofq6+tlt9vNLgfokHkv7tXu0mqfCclBNpsmJ0VrY+bEdn0HS84BBKKO/v2mpwcIYBdWYF3s6yuw2hNiEqMJOwD6Pub0AH6iM3dFZgUWALQfPT2AyboyJ4cVWADQfvT0ACbryl2RWYEFAO1H6AG6UUeHqLrjrsg89BMA2ofhLaAbdHaIqj1zcq7UW8NDPwGgfejpAb6mM5OJpc4PUXXnnBwe+gkAbaOnB1DXJhN3Zdn4hTk5l7vPDgEGALoPPT3oU3q7p0bq+rJx5uQAQO+gpwd+56i7Qcdrmzo0N8Wsnhqp60NUzMkBgN5B6EG360xokboWXNrqqbnSoxi6Opm4u4aouCsyAPQshrf6qM4O83Tl83VNLZr34l7d8tsCzV//oW5+apfmvbhX9U1n2/X5zg4xdXXZd3dMJmaICgD8Hz09PaizPR5d+XxXeku6+vmu9LZ0ZYjJH3pqGKICAP/XZ3p6nn/+eQ0fPlwDBgxQamqq9u7da1otXe3x6MrnuzIhtyuf72pvS1cmA/tTTw3LxgHAf/WJ0POHP/xB2dnZevTRR7V//36NGzdO6enpqqqqMqWeQA0eXfl8V1cwdSW4dMejGC701Ox8aIrWz79ROx+aoo2ZE9vVOwYACAx9IvQ8/fTTWrBggebPn6+UlBStW7dOYWFh+v3vf9/rtQRy8DCzt6WrwYWeGgDAlQT8nJ6Wlhbt27dPOTk53rZ+/fpp2rRpKioquuRnmpub1dzc7H3v8Xi6rZ6uzi/pyue7Gjy6o7elK/Ni8uaO19L8Az5ze9obXJhTAwC4koAPPdXV1WptbVVcXJxPe1xcnD7//PNLfiY3N1ePPfZYj9QTyMGjq5/vSmiRuie4sOwbAHA5fWJ4q6NycnJUX1/vfZ04caLbvrurwzRmD/N05fPdNS+GISYAQE8I+J6e6OhoBQUFqbKy0qe9srJSTqfzkp8JDQ1VaGhoj9XU1R4PM4d56G0BAPRVNsO4aMZsAEpNTdXEiROVl5cnSTp//ryGDh2qJUuW6Oc///kVP+/xeORwOFRfXy+73d5tdXV1fgnzUwAAuLyO/v0O+J4eScrOzlZGRoZuuOEGTZw4Uc8++6waGxs1f/58U+vqao8HPSYAAHSfPhF67rrrLrndbj3yyCOqqKjQ9ddfr7feeusbk5sBAIB19Ynhra7qqeEtAADQczr699uSq7cAAID1EHoAAIAlEHoAAIAlEHoAAIAlEHoAAIAlEHoAAIAlEHoAAIAlEHoAAIAl9Ik7MnfVhfszejwekysBAADtdeHvdnvvs0zokXT69GlJUkJCgsmVAACAjjp9+rQcDscV9+MxFPrqqewnT57UwIEDZbPZuu17PR6PEhISdOLECR5v0QGct87hvHUc56xzOG+dw3nrnLbOm2EYOn36tOLj49Wv35Vn7NDTI6lfv34aMmRIj32/3W7nAu8EzlvncN46jnPWOZy3zuG8dc7lzlt7enguYCIzAACwBEIPAACwBEJPDwoNDdWjjz6q0NBQs0sJKJy3zuG8dRznrHM4b53Deeuc7jxvTGQGAACWQE8PAACwBEIPAACwBEIPAACwBEIPAACwBEJPD3r++ec1fPhwDRgwQKmpqdq7d6/ZJfm1X/3qV7LZbD6v0aNHm12WXyksLNTtt9+u+Ph42Ww2vfrqqz7bDcPQI488osGDB+uqq67StGnTVFJSYk6xfuRK5+3+++//xrU3Y8YMc4r1E7m5ubrxxhs1cOBAxcbG6s4779Thw4d99jlz5oyysrIUFRWliIgIzZkzR5WVlSZV7B/ac96mTJnyjett0aJFJlXsH9auXauxY8d6b0Docrn05ptverd317VG6Okhf/jDH5Sdna1HH31U+/fv17hx45Senq6qqiqzS/Nr1157rU6dOuV9/fnPfza7JL/S2NiocePG6fnnn7/k9jVr1ui5557TunXrtGfPHoWHhys9PV1nzpzp5Ur9y5XOmyTNmDHD59rLz8/vxQr9T0FBgbKysvTBBx9o+/btOnv2rKZPn67GxkbvPsuXL9frr7+uLVu2qKCgQCdPntTs2bNNrNp87TlvkrRgwQKf623NmjUmVewfhgwZotWrV2vfvn366KOPdMstt+iOO+7QoUOHJHXjtWagR0ycONHIysryvm9tbTXi4+ON3NxcE6vyb48++qgxbtw4s8sIGJKMrVu3et+fP3/ecDqdxpNPPultq6urM0JDQ438/HwTKvRPF583wzCMjIwM44477jClnkBRVVVlSDIKCgoMw/jq2goODja2bNni3eezzz4zJBlFRUVmlel3Lj5vhmEY3//+942f/OQn5hUVIK6++mrjP/7jP7r1WqOnpwe0tLRo3759mjZtmretX79+mjZtmoqKikyszP+VlJQoPj5eI0aM0L333qvy8nKzSwoYZWVlqqio8LnuHA6HUlNTue7aYdeuXYqNjdWoUaO0ePFi1dTUmF2SX6mvr5ckRUZGSpL27duns2fP+lxvo0eP1tChQ7nevubi83bByy+/rOjoaF133XXKyclRU1OTGeX5pdbWVm3evFmNjY1yuVzdeq3xwNEeUF1drdbWVsXFxfm0x8XF6fPPPzepKv+XmpqqDRs2aNSoUTp16pQee+wx3XTTTfrkk080cOBAs8vzexUVFZJ0yevuwjZc2owZMzR79mwlJibqyJEj+sUvfqGZM2eqqKhIQUFBZpdnuvPnz+vBBx/U5MmTdd1110n66noLCQnRoEGDfPblevu7S503Sbrnnns0bNgwxcfH6+OPP9bKlSt1+PBhvfLKKyZWa76DBw/K5XLpzJkzioiI0NatW5WSkqLi4uJuu9YIPfAbM2fO9P577NixSk1N1bBhw/THP/5RmZmZJlaGvu7uu+/2/nvMmDEaO3asRo4cqV27dmnq1KkmVuYfsrKy9MknnzDHroMud94WLlzo/feYMWM0ePBgTZ06VUeOHNHIkSN7u0y/MWrUKBUXF6u+vl5/+tOflJGRoYKCgm79DYa3ekB0dLSCgoK+MbO8srJSTqfTpKoCz6BBg/Ttb39bpaWlZpcSEC5cW1x3XTdixAhFR0dz7UlasmSJtm3bpp07d2rIkCHedqfTqZaWFtXV1fnsz/X2lcudt0tJTU2VJMtfbyEhIUpKStKECROUm5urcePG6Xe/+123XmuEnh4QEhKiCRMmaMeOHd628+fPa8eOHXK5XCZWFlgaGhp05MgRDR482OxSAkJiYqKcTqfPdefxeLRnzx6uuw764osvVFNTY+lrzzAMLVmyRFu3btV7772nxMREn+0TJkxQcHCwz/V2+PBhlZeXW/p6u9J5u5Ti4mJJsvT1dinnz59Xc3Nz915r3TvXGhds3rzZCA0NNTZs2GB8+umnxsKFC41BgwYZFRUVZpfmt376058au3btMsrKyozdu3cb06ZNM6Kjo42qqiqzS/Mbp0+fNg4cOGAcOHDAkGQ8/fTTxoEDB4zjx48bhmEYq1evNgYNGmS89tprxscff2zccccdRmJiovHll1+aXLm52jpvp0+fNh566CGjqKjIKCsrM959913jO9/5jpGcnGycOXPG7NJNs3jxYsPhcBi7du0yTp065X01NTV591m0aJExdOhQ47333jM++ugjw+VyGS6Xy8SqzXel81ZaWmqsWrXK+Oijj4yysjLjtddeM0aMGGGkpaWZXLm5fv7znxsFBQVGWVmZ8fHHHxs///nPDZvNZrzzzjuGYXTftUbo6UF5eXnG0KFDjZCQEGPixInGBx98YHZJfu2uu+4yBg8ebISEhBjf+ta3jLvuussoLS01uyy/snPnTkPSN14ZGRmGYXy1bP2Xv/ylERcXZ4SGhhpTp041Dh8+bG7RfqCt89bU1GRMnz7diImJMYKDg41hw4YZCxYssPx/UC51viQZ69ev9+7z5ZdfGv/8z/9sXH311UZYWJjxox/9yDh16pR5RfuBK5238vJyIy0tzYiMjDRCQ0ONpKQkY8WKFUZ9fb25hZvsgQceMIYNG2aEhIQYMTExxtSpU72BxzC671qzGYZhdLLnCQAAIGAwpwcAAFgCoQcAAFgCoQcAAFgCoQcAAFgCoQcAAFgCoQcAAFgCoQcAAFgCoQcAAFgCoQdAnzFlyhQ9+OCDZpcBwE8RegAAgCUQegAAgCUQegAEpMbGRs2bN08REREaPHiwfvvb3/ps/8///E/dcMMNGjhwoJxOp+655x5VVVVJkgzDUFJSkp566imfzxQXF8tms6m0tLTXjgNA7yH0AAhIK1asUEFBgV577TW988472rVrl/bv3+/dfvbsWT3++OP6y1/+oldffVXHjh3T/fffL0my2Wx64IEHtH79ep/vXL9+vdLS0pSUlNSbhwKgl/CUdQABp6GhQVFRUfqv//ov/cM//IMkqba2VkOGDNHChQv17LPPfuMzH330kW688UadPn1aEREROnnypIYOHar3339fEydO1NmzZxUfH6+nnnpKGRkZvXxEAHoDPT0AAs6RI0fU0tKi1NRUb1tkZKRGjRrlfb9v3z7dfvvtGjp0qAYOHKjvf//7kqTy8nJJUnx8vGbNmqXf//73kqTXX39dzc3N3hAFoO8h9ADocxobG5Weni673a6XX35ZH374obZu3SpJamlp8e73T//0T9q8ebO+/PJLrV+/XnfddZfCwsLMKhtADyP0AAg4I0eOVHBwsPbs2eNt+9vf/qa//vWvkqTPP/9cNTU1Wr16tW666SaNHj3aO4n562699VaFh4dr7dq1euutt/TAAw/02jEA6H39zS4AADoqIiJCmZmZWrFihaKiohQbG6v/9//+n/r1++r/cUOHDlVISIjy8vK0aNEiffLJJ3r88ce/8T1BQUG6//77lZOTo+TkZLlcrt4+FAC9iJ4eAAHpySef1E033aTbb79d06ZN0/e+9z1NmDBBkhQTE6MNGzZoy5YtSklJ0erVq7+xPP2CzMxMtbS0aP78+b1ZPgATsHoLgKX97//+r6ZOnaoTJ04oLi7O7HIA9CBCDwBLam5ultvtVkZGhpxOp15++WWzSwLQwxjeAmBJ+fn5GjZsmOrq6rRmzRqzywHQC+jpAQAAlkBPDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsARCDwAAsIT/Dza4v9vyNVYOAAAAAElFTkSuQmCC",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"observed_data.plot.scatter(\"day\", \"dotI\")"
]
},
{
"cell_type": "markdown",
"id": "d6d122c0-2139-4b89-8293-33076561bd78",
"metadata": {},
"source": [
"## Sobol Sensitivity Analysis\n",
"\n",
"Next, let's use `calisim` to perform a sensitivity analysis of the end-of-simulation number of infected. To start with, we'll need to define our `ParameterSpecification` parameter specification:"
]
},
{
"cell_type": "code",
"execution_count": 104,
"id": "849fa816-b458-4db3-9bb7-cd527cf9db87",
"metadata": {},
"outputs": [],
"source": [
"parameter_spec = ParameterSpecification(\n",
"\tparameters=[\n",
"\t\tDistributionModel(\n",
"\t\t\tname=\"beta\",\n",
"\t\t\tdistribution_name=\"uniform\",\n",
"\t\t\tdistribution_args=[0.3, 0.5],\n",
"\t\t\tdata_type=ParameterDataType.CONTINUOUS,\n",
"\t\t),\n",
"\t\tDistributionModel(\n",
"\t\t\tname=\"gamma\",\n",
"\t\t\tdistribution_name=\"uniform\",\n",
"\t\t\tdistribution_args=[0.05, 0.15],\n",
"\t\t\tdata_type=ParameterDataType.CONTINUOUS,\n",
"\t\t),\n",
"\t]\n",
")"
]
},
{
"cell_type": "markdown",
"id": "ed4e07cc-5eee-4cd7-8a83-e990744e52cb",
"metadata": {},
"source": [
"This contains information concerning the various parameter names, probability distributions, ranges, distribution parameters, and data types.\n",
"\n",
"We next need to create a wrapper function around our forward model to ensure there's compatibility with the `calisim` API."
]
},
{
"cell_type": "code",
"execution_count": 105,
"id": "efbf14e8-3c14-449a-bbee-6d55e7c2390d",
"metadata": {},
"outputs": [],
"source": [
"def sensitivity_func(\n",
"\tparameters: dict, simulation_id: str, observed_data: np.ndarray | None, t: pd.Series\n",
") -> float | list[float]:\n",
" simulation_parameters = model.GROUND_TRUTH.copy()\n",
" simulation_parameters[\"t\"] = t\n",
"\n",
" for k in [\"beta\", \"gamma\"]:\n",
" simulation_parameters[k] = parameters[k]\n",
"\n",
" simulated_data = sir_simulate(simulation_parameters).tail(1).dotI.item()\n",
" return simulated_data"
]
},
{
"cell_type": "markdown",
"id": "1e922ada-5b5c-4a41-a2cc-424d6f729671",
"metadata": {},
"source": [
"The last step is to create a `SensitivityAnalysisMethodModel` specification for the calibration procedure itself, which we then supply to a `SensitivityAnalysisMethod` calibrator. We'll use the `Sobol` method via the [SALib engine](https://salib.readthedocs.io/en/latest/), and calculate 95% confidence intervals for the indices."
]
},
{
"cell_type": "code",
"execution_count": 106,
"id": "0b1cf8b6-48df-409a-91fa-d156cd939068",
"metadata": {},
"outputs": [],
"source": [
"specification = SensitivityAnalysisMethodModel(\n",
"\texperiment_name=\"salib_sensitivity_analysis\",\n",
"\tparameter_spec=parameter_spec,\n",
" observed_data=observed_data.dotI.values,\n",
"\tmethod=\"sobol\",\n",
"\tn_samples=256,\n",
"\toutput_labels=[\"Number of Infected\"],\n",
" calibration_func_kwargs=dict(t=observed_data.day),\n",
"\tmethod_kwargs=dict(calc_second_order=False, scramble=True),\n",
"\tanalyze_kwargs=dict(\n",
"\t\tcalc_second_order=False,\n",
"\t\tnum_resamples=300,\n",
"\t\tconf_level=0.95,\n",
"\t),\n",
")\n",
"\n",
"calibrator = SensitivityAnalysisMethod(\n",
"\tcalibration_func=sensitivity_func, specification=specification, engine=\"salib\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "d37e2a18-3ae0-4e12-b572-047b644b12c4",
"metadata": {},
"source": [
"Finally, we'll run the calibration procedure. This is composed of 3 steps:\n",
"\n",
"1. **Specify**: Define your calibration problem: Parameter distributions, observed data, objective/discrepancy function, and calibration settings (like algorithm, directions, iterations)\n",
"2. **Execute**: Run the actual calibration process (simulation + optimization/inference)\n",
"3. **Analyze**: Process, summarize, and optionally save plots/metrics of the calibration results\n",
"\n",
"Or **SEA**."
]
},
{
"cell_type": "code",
"execution_count": 107,
"id": "a2aad087-c580-488e-8424-39bc276be210",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/mnt/c/Users/james/projects/model-calibration-toolbox/.venv/lib/python3.10/site-packages/SALib/util/__init__.py:274: FutureWarning:\n",
"\n",
"unique with argument that is not not a Series, Index, ExtensionArray, or np.ndarray is deprecated and will raise in a future version.\n",
"\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHVCAYAAACXAw0nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAmqElEQVR4nO3dfXBV9Z0/8E8SSALlScUkyAaprQ/FCriwUMq2rh0Uu64PO9Wh6ggidVdXfCB1FrBKfthqXK3I1qVFqQiO60p11tWpFnWpdLcVZQTd6o7iqkUQDQ8yEE0KSHJ/f1giKcFyk5B78+X1mjkz5Nxzcj65c/nM+5zz/Z5bkMlkMgEAQJdWmOsCAABoP6EOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJKBbrgs4EE1NTfHee+9F7969o6CgINflAHkmk8nEhx9+GEcddVQUFnatc1X9DfhTDrTHdYlQ995770VlZWWuywDy3Pr16+PP/uzPcl1GVvQ34ED9qR7XJUJd7969I+KTP6ZPnz45rgbIN3V1dVFZWdncK7oS/Q34Uw60x3WJULfnlkSfPn00PWC/uuLtS/0NOFB/qsd1rcEnAAC0SqgDAEiAUAcAkIAuMaYOiGhsbIyPP/4412XkTPfu3aOoqCjXZQAHgf7WMf1NqIM8l8lkora2NrZt25brUnKuX79+UVFR0SUnRAD70t8+1RH9TaiDPLen4ZWVlUXPnj0PyUCTyWSioaEhNm3aFBERAwYMyHFFQEfQ3zq2vwl1kMcaGxubG94RRxyR63JyqkePHhERsWnTpigrK3MrFro4/e1THdXfTJSAPLZnjEnPnj1zXEl+2PM+HMpjbyAV+ltLHdHfhDroAg7FWxKt8T5Aevy//kRHvA9CHQBAAoQ6AIAEmCgBXdTgGU906vHW3npmpx4POHTpb23jSh1w0GzevDmuuOKKGDRoUJSUlERFRUWMHz8+br755igoKPjMZfny5bkuH2C/9tfffvOb30RExD333BN/9Vd/FX369ImCgoJOeRafK3XAQfOtb30rdu3aFYsXL45jjjkmNm7cGMuWLYsTTzwx3n///ebtrrnmmqirq4v77ruved3hhx+ei5IBDsj++tsHH3wQERENDQ1xxhlnxBlnnBEzZ87slJqEOuCg2LZtW/z3f/93LF++PE455ZSIiDj66KNj1KhR+2zbo0eP2LlzZ1RUVHR2mQBZO5D+du2110ZEdOpdB7dfgYOiV69e0atXr/iP//iP2LlzZ67LAegw+drfhDrgoOjWrVssWrQoFi9eHP369YuxY8fG9ddfH7/97W9zXRpAu+RrfxPqgIPmW9/6Vrz33nvx+OOPxxlnnBHLly+PP//zP49FixblujSAdsnH/ibUAQdVaWlpnHbaaXHjjTfGc889F5dccklUV1fnuiyAdsu3/ibUAZ1qyJAhUV9fn+syADpcrvub2a/AQfHBBx/E+eefH5deemkMHTo0evfuHS+++GLcdtttcc455+S6PIA2O5D+VltbG7W1tfHmm29GRMQrr7wSvXv3jkGDBh20RzYJddBFvXXzGfHSSy9FRMTJJ58cRUVFOa6opV69esXo0aPjzjvvjLfeeis+/vjjqKysjMsuuyyuv/76XJcH5LG1t54ZjY2NedvjDqS/zZ8/P2bPnt28z9e//vWIiLjvvvvikksuOSh1FWQymcxB+c0dqK6uLvr27Rvbt2+PPn365Loc6DQ7duyI3/3ud/H5z38+SktLW7yWzw3vYNnf+9GVe0RXrh3a47P6W8Sh1+M+6/040D5hTB0AQAKEOgCABAh1AAAJEOqgC+gCQ187hfcB0uP/9Sc64n0Q6iCPde/ePSIiGhoaclxJftjzPux5X4CuS39rqSP6m0eaQB4rKiqKfv36xaZNmyIiomfPnlFQUBARn8wM22PHjh1JzwzLZDLR0NAQmzZtin79+iX9t8Kh4rP6W8Sh0+M6sr8JdZDnKioqIiKaG98eTU1NsWXLloiIWLt2bRQWpn/hvV+/fs3vB9D17a+/RRx6Pa4j+ptQB3muoKAgBgwYEGVlZfHxxx83r29oaIgzzzwzIiJWr14dPXv2zFWJnaJ79+7JnqnDoWp//S3i0OpxHdXfhDroIoqKilr8p29sbIx33nknIiJKSkpafXgnQFfwx/0tQo9ri7SvZQIAHCKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOYC/z5s2LwYMHR2lpaYwePTpWrlz5mdvPnTs3jj/++OjRo0dUVlbGtGnTYseOHZ1ULcCnhDqAP1iyZElUVVVFdXV1rF69OoYNGxbjx4+PTZs2tbr9gw8+GDNmzIjq6up47bXX4t57740lS5bE9ddf38mVAwh1AM3mzJkTl112WUyePDmGDBkS8+fPj549e8bChQtb3f65556LsWPHxoUXXhiDBw+O008/PS644ILPvLq3c+fOqKura7EAdAShDiAidu3aFatWrYpx48Y1ryssLIxx48bFihUrWt3nq1/9aqxatao5xL399tvx5JNPxl//9V/v9zg1NTXRt2/f5qWysrJj/xDgkNUt1wUA5IMtW7ZEY2NjlJeXt1hfXl4er7/+eqv7XHjhhbFly5b4y7/8y8hkMrF79+64/PLLP/P268yZM6Oqqqr557q6OsEO6BCu1AG00fLly+OWW26JH//4x7F69er493//93jiiSfi+9///n73KSkpiT59+rRYADqCK3UAEdG/f/8oKiqKjRs3tli/cePGqKioaHWfG2+8MS6++OL4zne+ExERJ510UtTX18ff/d3fxfe+970oLHTeDHQeHQcgIoqLi2PEiBGxbNmy5nVNTU2xbNmyGDNmTKv7NDQ07BPcioqKIiIik8kcvGIBWuFKHcAfVFVVxaRJk2LkyJExatSomDt3btTX18fkyZMjImLixIkxcODAqKmpiYiIs846K+bMmRMnn3xyjB49Ot5888248cYb46yzzmoOdwCdpU1X6jycE0jRhAkT4oc//GHMmjUrhg8fHi+//HIsXbq0efLEunXr4v3332/e/oYbbojvfve7ccMNN8SQIUNiypQpMX78+Lj77rtz9ScAh7CCTJb3CJYsWRITJ06M+fPnx+jRo2Pu3Lnx8MMPx5o1a6KsrGyf7R988MG49NJLY+HChfHVr3413njjjbjkkkvi29/+dsyZM+eAjllXVxd9+/aN7du3G1QMf1BfXx+9evWKiIiPPvooPve5z+W4otzpyj2iK9cOB5Me96kD7RNZX6nzcE4AgPyT1Zi6PQ/nnDlzZvO6A3k45wMPPBArV66MUaNGNT+c8+KLL97vcWpqamL27NnZlAY5MXjGEzk7dtOuT4cwfOnGpVFYXJqzWtbeembOjg3AJ7IKdR7OCQCQnw76I008nBMA4ODL6kqdh3MCwKEhl8NLIgwxaYusEpWHcwIA5KesHz7s4ZwAAPkn61A3YcKE2Lx5c8yaNStqa2tj+PDh+zycc+8rczfccEMUFBTEDTfcEBs2bIgjjzwyzjrrrLj55ps77q8AADjEtelrwqZOnRpTp05t9bXly5e3PEC3blFdXR3V1dVtORQAAAfALAUAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkoFuuCwDaprC4NI6e/vNclwFAnnClDgAgAUIdAEAChDoAgAQIdQAACRDqAAASYPYrAJB3zPDPnit1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAexl3rx5MXjw4CgtLY3Ro0fHypUrP3P7bdu2xZVXXhkDBgyIkpKSOO644+LJJ5/spGoBPuW7XwH+YMmSJVFVVRXz58+P0aNHx9y5c2P8+PGxZs2aKCsr22f7Xbt2xWmnnRZlZWXxyCOPxMCBA+Odd96Jfv36dX7xwCFPqAP4gzlz5sRll10WkydPjoiI+fPnxxNPPBELFy6MGTNm7LP9woULY+vWrfHcc89F9+7dIyJi8ODBnVkyQDO3XwHik6tuq1atinHjxjWvKywsjHHjxsWKFSta3efxxx+PMWPGxJVXXhnl5eXx5S9/OW655ZZobGzc73F27twZdXV1LRaAjiDUAUTEli1borGxMcrLy1usLy8vj9ra2lb3efvtt+ORRx6JxsbGePLJJ+PGG2+MO+64I37wgx/s9zg1NTXRt2/f5qWysrJD/w7g0CXUAbRRU1NTlJWVxT333BMjRoyICRMmxPe+972YP3/+fveZOXNmbN++vXlZv359J1YMpMyYOoCI6N+/fxQVFcXGjRtbrN+4cWNUVFS0us+AAQOie/fuUVRU1LzuS1/6UtTW1sauXbuiuLh4n31KSkqipKSkY4sHiDZeqTPlH0hNcXFxjBgxIpYtW9a8rqmpKZYtWxZjxoxpdZ+xY8fGm2++GU1NTc3r3njjjRgwYECrgQ7gYMo61O2Z8l9dXR2rV6+OYcOGxfjx42PTpk2tbr9nyv/atWvjkUceiTVr1sSCBQti4MCB7S4eoCNVVVXFggULYvHixfHaa6/FFVdcEfX19c2zYSdOnBgzZ85s3v6KK66IrVu3xjXXXBNvvPFGPPHEE3HLLbfElVdemas/ATiEZX37tTOm/O/cuTN27tzZ/LPZYUBnmDBhQmzevDlmzZoVtbW1MXz48Fi6dGnz5Il169ZFYeGn58KVlZXx1FNPxbRp02Lo0KExcODAuOaaa2L69Om5+hOAQ1hWoW7PlP+9z1SzmfL/2GOPxZFHHhkXXnhhTJ8+vcU4lL3V1NTE7NmzsykNoENMnTo1pk6d2upry5cv32fdmDFj4vnnnz/IVQH8aVndfu2sKf9mhwEAZOegz37de8p/UVFRjBgxIjZs2BC33357VFdXt7qP2WEAANnJKtR11pR/AACyk9XtV1P+AQDyU9aPNDHlHwAg/2Q9ps6UfwCA/NOmiRKm/AMA5Jc2fU0YAAD5RagDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AHsZd68eTF48OAoLS2N0aNHx8qVKw9ov4ceeigKCgri3HPPPbgFAuyHUAfwB0uWLImqqqqorq6O1atXx7Bhw2L8+PGxadOmz9xv7dq1cd1118XXvva1TqoUYF9CHcAfzJkzJy677LKYPHlyDBkyJObPnx89e/aMhQsX7nefxsbGuOiii2L27NlxzDHH/Mlj7Ny5M+rq6losAB1BqAOIiF27dsWqVati3LhxzesKCwtj3LhxsWLFiv3ud9NNN0VZWVlMmTLlgI5TU1MTffv2bV4qKyvbXTtARBtDnTEnQGq2bNkSjY2NUV5e3mJ9eXl51NbWtrrPr3/967j33ntjwYIFB3ycmTNnxvbt25uX9evXt6tugD2yDnXGnABEfPjhh3HxxRfHggULon///ge8X0lJSfTp06fFAtARsg51nTHmBKCz9e/fP4qKimLjxo0t1m/cuDEqKir22f6tt96KtWvXxllnnRXdunWLbt26xf333x+PP/54dOvWLd56663OKh0gIrIMdZ015sRAYqCzFRcXx4gRI2LZsmXN65qammLZsmUxZsyYfbY/4YQT4pVXXomXX365eTn77LPj1FNPjZdfftlYOaDTdctm488ac/L666+3us+eMScvv/zyAR+npqYmZs+enU1pAO1WVVUVkyZNipEjR8aoUaNi7ty5UV9fH5MnT46IiIkTJ8bAgQOjpqYmSktL48tf/nKL/fv16xcRsc96gM6QVajLVlvHnMycOTOqqqqaf66rq3PWCxx0EyZMiM2bN8esWbOitrY2hg8fHkuXLm0+kV23bl0UFnpoAJCfsgp17RlzskdTU9MnB+7WLdasWRNf+MIX9tmvpKQkSkpKsikNoENMnTo1pk6d2upry5cv/8x9Fy1a1PEFARygrE45jTkBAMhPWd9+NeYEACD/ZB3qjDkBAMg/bZooYcwJAEB+cUkNACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AGQc/X19VFQUBAFBQVRX1+f63KgSxLqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoa6L8GXXAMBnEeoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAR0y3UBXcngGU/k7NhNu3Y0//tLNy6NwuLSnNWy9tYzc3ZsAKB1rtQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJ8Jw6gL3Mmzcvbr/99qitrY1hw4bFXXfdFaNGjWp12wULFsT9998fr776akREjBgxIm655Zb9bp/PcvkczgjP4oSO4EodwB8sWbIkqqqqorq6OlavXh3Dhg2L8ePHx6ZNm1rdfvny5XHBBRfEs88+GytWrIjKyso4/fTTY8OGDZ1cOUAbQ928efNi8ODBUVpaGqNHj46VK1fud9sFCxbE1772tTjssMPisMMOi3Hjxn3m9gC5MmfOnLjsssti8uTJMWTIkJg/f3707NkzFi5c2Or2//qv/xr/8A//EMOHD48TTjghfvrTn0ZTU1MsW7askysHaEOocyYLpGjXrl2xatWqGDduXPO6wsLCGDduXKxYseKAfkdDQ0N8/PHHcfjhh+93m507d0ZdXV2LBaAjZB3qnMkCKdqyZUs0NjZGeXl5i/Xl5eVRW1t7QL9j+vTpcdRRR7UIhn+spqYm+vbt27xUVla2q26APbIKdc5kAVp36623xkMPPRSPPvpolJbuf5D/zJkzY/v27c3L+vXrO7FKIGVZhTpnskCq+vfvH0VFRbFx48YW6zdu3BgVFRWfue8Pf/jDuPXWW+Ppp5+OoUOHfua2JSUl0adPnxYLQEfo1NmvzmSBfFVcXBwjRoxoMTRkz1CRMWPG7He/2267Lb7//e/H0qVLY+TIkZ1RKkCrsnpOXUecyf7nf/7nAZ3JlpSUZFMaQLtVVVXFpEmTYuTIkTFq1KiYO3du1NfXx+TJkyMiYuLEiTFw4MCoqamJiIh/+qd/ilmzZsWDDz4YgwcPbr5j0atXr+jVq1fO/g7g0JRVqNv7TPbcc8+NiE/PZKdOnbrf/W677ba4+eab46mnnnIm20aFxaVx9PSf57oMSNqECRNi8+bNMWvWrKitrY3hw4fH0qVLm4ecrFu3LgoLP73B8ZOf/CR27doV5513XovfU11dHf/v//2/ziwdIPtvlHAmC6Rs6tSp+z1JXb58eYuf165de/ALAjhAWYc6Z7IAAPmnTd/96kwWACC/+O5XAIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAArrlugAAKCwujaOn/zzXZUCX5kodAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEiDUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASIBQBwCQAKEOACABQh0AQAKEOgCABAh1AAAJEOoAABIg1AEAJECoAwBIgFAHAJAAoQ4AIAFCHQBAAoQ6AIAECHUAAAkQ6gAAEtCmUDdv3rwYPHhwlJaWxujRo2PlypWfuf3DDz8cJ5xwQpSWlsZJJ50UTz75ZJuKBTjY9Degq8o61C1ZsiSqqqqiuro6Vq9eHcOGDYvx48fHpk2bWt3+ueeeiwsuuCCmTJkSL730Upx77rlx7rnnxquvvtru4gE6kv4GdGUFmUwmk80Oo0ePjr/4i7+If/mXf4mIiKampqisrIyrrroqZsyYsc/2EyZMiPr6+vj5z3/evO4rX/lKDB8+PObPn9/qMXbu3Bk7d+5s/nn79u0xaNCgWL9+ffTp0yebcjvUl6ufytmx88mrs8fnuoS84TPxiVx/Jurq6qKysjK2bdsWffv2bfPv0d+IyP3nOV/4THwq15+JA+5xmSzs3LkzU1RUlHn00UdbrJ84cWLm7LPPbnWfysrKzJ133tli3axZszJDhw7d73Gqq6szEWGxWCxZLevXr8+mpelvFoulSy1/qsd1iyxs2bIlGhsbo7y8vMX68vLyeP3111vdp7a2ttXta2tr93ucmTNnRlVVVfPPTU1NsXXr1jjiiCOioKAgm5KTsiep5/qMnvzhM/GJTCYTH374YRx11FFt/h36W+75PPPHfCY+caA9LqtQ11lKSkqipKSkxbp+/frlppg81KdPn0P6w82+fCaiXbddO5P+9qf5PPPHfCYOrMdlNVGif//+UVRUFBs3bmyxfuPGjVFRUdHqPhUVFVltD5AL+hvQ1WUV6oqLi2PEiBGxbNmy5nVNTU2xbNmyGDNmTKv7jBkzpsX2ERHPPPPMfrcHyAX9Dejqsr79WlVVFZMmTYqRI0fGqFGjYu7cuVFfXx+TJ0+OiIiJEyfGwIEDo6amJiIirrnmmjjllFPijjvuiDPPPDMeeuihePHFF+Oee+7p2L/kEFBSUhLV1dX73Lrh0OUz0bH0t9zyeeaP+Uxk6QAmhe3jrrvuygwaNChTXFycGTVqVOb5559vfu2UU07JTJo0qcX2P/vZzzLHHXdcpri4OHPiiSdmnnjiibYcFuCg09+Arirr59QBAJB/fPcrAEAChDoAgAQIdQAACRDqAAASINQBACRAqAMASEBefvcr+2poaIh169bFrl27WqwfOnRojioilx555JH42c9+1upnYvXq1TmqCtpOj2MP/a3tXKnLc5s3b46/+Zu/id69e8eJJ54YJ598couFQ8+PfvSjmDx5cpSXl8dLL70Uo0aNiiOOOCLefvvt+OY3v5nr8iArehx709/aR6jLc9dee21s27YtXnjhhejRo0csXbo0Fi9eHMcee2w8/vjjuS6PHPjxj38c99xzT9x1111RXFwc//iP/xjPPPNMXH311bF9+/ZclwdZ0ePYm/7WTrn+Sgs+W0VFReaFF17IZDKZTO/evTNr1qzJZDKZzGOPPZYZO3ZsLksjR3r06JFZu3ZtJpPJZI488sjMyy+/nMlkMpk33ngjc/jhh+eyNMiaHsfe9Lf2caUuz9XX10dZWVlERBx22GGxefPmiIg46aSTjC04RFVUVMTWrVsjImLQoEHx/PPPR0TE7373u8j41j+6GD2Ovelv7SPU5bnjjz8+1qxZExERw4YNi7vvvjs2bNgQ8+fPjwEDBuS4OnLhG9/4RvNtqcmTJ8e0adPitNNOiwkTJsTf/u3f5rg6yI4ex970t/YpyIi+ee2BBx6I3bt3xyWXXBKrVq2KM844I7Zu3RrFxcWxaNGimDBhQq5LpJM1NTVFU1NTdOv2yeT1hx56KJ577rk49thj4+///u+juLg4xxXCgdPj2Jv+1j5CXRfT0NAQr7/+egwaNCj69++f63IAOpQeB20n1OW5m266Ka677rro2bNni/W///3v4/bbb49Zs2blqDJyaceOHfHb3/42Nm3aFE1NTS1eO/vss3NUFWRPj+OP6W9tJ9TluaKionj//febBxLv8cEHH0RZWVk0NjbmqDJyZenSpTFx4sTYsmXLPq8VFBT4TNCl6HHsTX9rHxMl8lwmk4mCgoJ91v/P//xPHH744TmoiFy76qqr4vzzz4/333+/efzJnkXDo6vR49ib/tY+viYsTx122GFRUFAQBQUFcdxxx7Voeo2NjfHRRx/F5ZdfnsMKyZWNGzdGVVVVlJeX57oUaDM9jtbob+0j1OWpuXPnRiaTiUsvvTRmz54dffv2bX6tuLg4Bg8eHGPGjMlhheTKeeedF8uXL48vfOELuS4F2kyPozX6W/sYU5fnfvWrX8XYsWObp3dDQ0NDnH/++XHkkUfGSSedFN27d2/x+tVXX52jyiB7ehx709/aR6jrAt56662477774q233op//ud/jrKysvjFL34RgwYNihNPPDHX5dHJ7r333rj88sujtLQ0jjjiiBa3rQoKCuLtt9/OYXWQPT2OPfS39hHq8tyvfvWr+OY3vxljx46N//qv/4rXXnstjjnmmLj11lvjxRdfjEceeSTXJdLJKioq4uqrr44ZM2ZEYaG5TnRtehx709/axzuW52bMmBE/+MEP4plnnmnxJO1vfOMbzd+Jx6Fl165dMWHCBA2PJOhx7E1/ax/vWp575ZVXWv2+u7Kyslaf40P6Jk2aFEuWLMl1GdAh9Dj2pr+1j5Gpea5fv37x/vvvx+c///kW61966aUYOHBgjqoilxobG+O2226Lp556KoYOHbrPQOI5c+bkqDLInh7H3vS39hHq8ty3v/3tmD59ejz88MNRUFAQTU1N8Zvf/Cauu+66mDhxYq7LIwdeeeWVOPnkkyMi4tVXX23xWmsPcYV8psexN/2tfUyUyHO7du2KK6+8MhYtWhSNjY3RrVu32L17d1x00UWxaNGiKCoqynWJAG2mx0HHEeq6iPXr18crr7wS9fX1cfLJJ8cXv/jFXJcE0GH0OGg/t1+7gHvvvTfuvPPO+L//+7+IiDj22GPj2muvje985zs5roxc2LFjR9x1113x7LPPxqZNm6KpqanF66tXr85RZdA2ehx76G/tI9TluVmzZsWcOXPiqquuav7KnBUrVsS0adNi3bp1cdNNN+W4QjrblClT4umnn47zzjsvRo0aZZwJXZoex970t/Zx+zXPHXnkkfGjH/0oLrjgghbr/+3f/i2uuuoqU/4PQX379o0nn3wyxo4dm+tSoN30OPamv7WP59TluY8//jhGjhy5z/oRI0bE7t27c1ARuTZw4MDo3bt3rsuADqHHsTf9rX2Eujx38cUXx09+8pN91t9zzz1x0UUX5aAicu2OO+6I6dOnxzvvvJPrUqDd9Dj2pr+1jzF1eaiqqqr53wUFBfHTn/40nn766fjKV74SEREvvPBCrFu3zjOcDlEjR46MHTt2xDHHHBM9e/bc5+GcW7duzVFlcGD0OPZHf2sfY+ry0KmnnnpA2xUUFMQvf/nLg1wN+WbcuHGxbt26mDJlSpSXl+8zkHjSpEk5qgwOjB7H/uhv7SPUQRfTs2fPWLFiRQwbNizXpQB0KP2tfYypgy7mhBNOiN///ve5LgOgw+lv7SPUQRdz6623xne/+91Yvnx5fPDBB1FXV9diAeiq9Lf2cfsVupjCwk/Oxf54rEkmk4mCgoJobGzMRVkA7aa/tY/Zr9DFPPvss7kuAeCg0N/ax5U6AIAEuFIHXVRDQ0OsW7cudu3a1WL90KFDc1QRQMfQ39pGqIMuZvPmzTF58uT4xS9+0errxpwAXZX+1j5mv0IXc+2118a2bdvihRdeiB49esTSpUtj8eLFceyxx8bjjz+e6/IA2kx/ax9X6qCL+eUvfxmPPfZYjBw5MgoLC+Poo4+O0047Lfr06RM1NTVx5pln5rpEgDbR39rHlTroYurr66OsrCwiIg477LDYvHlzREScdNJJsXr16lyWBtAu+lv7CHXQxRx//PGxZs2aiIgYNmxY3H333bFhw4aYP39+DBgwIMfVAbSd/tY+HmkCXcwDDzwQu3fvjksuuSRWrVoVZ5xxRnzwwQdRXFwcixcvjgkTJuS6RIA20d/aR6iDLq6hoSFef/31GDRoUPTv3z/X5QB0GP0tOyZKQBdTVVXV6vqCgoIoLS2NL37xi3HOOefE4Ycf3smVAbSP/tY+rtRBF3PqqafG6tWro7GxMY4//viIiHjjjTeiqKgoTjjhhFizZk0UFBTEr3/96xgyZEiOqwU4cPpb+5goAV3MOeecE+PGjYv33nsvVq1aFatWrYp33303TjvttLjgggtiw4YN8fWvfz2mTZuW61IBsqK/tY8rddDFDBw4MJ555pl9zlL/93//N04//fTYsGFDrF69Ok4//fTYsmVLjqoEyJ7+1j6u1EEXs3379ti0adM+6zdv3hx1dXUREdGvX799vjMRIN/pb+0j1EEXc84558Sll14ajz76aLz77rvx7rvvxqOPPhpTpkyJc889NyIiVq5cGccdd1xuCwXIkv7WPm6/Qhfz0UcfxbRp0+L++++P3bt3R0REt27dYtKkSXHnnXfG5z73uXj55ZcjImL48OG5KxQgS/pb+wh10EV99NFH8fbbb0dExDHHHBO9evXKcUUAHUN/axuhDgAgAcbUAQAkQKgDAEiAUAcAkAChDgAgAUIdAEAChDoAgAQIdQAACfj/wNRNvFLR92AAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/mnt/c/Users/james/projects/model-calibration-toolbox/.venv/lib/python3.10/site-packages/SALib/util/__init__.py:274: FutureWarning:\n",
"\n",
"unique with argument that is not not a Series, Index, ExtensionArray, or np.ndarray is deprecated and will raise in a future version.\n",
"\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnoAAAJOCAYAAAAgQaq3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1qklEQVR4nO3deXQUZb7/8U8lIQmBdAAhCcawLwoKUXITGRRQg8D1jrhhXIEouIyIkMv8BBxAHDEoVxYdNOIAYeCcEYerI3NAlBPNuIAywnDBBUQWE4GEgEIgDIl01e8PpMeYoOlQSaeffr/OqXPs6uqnvt0M+p3PU0+V5TiOIwAAABgnLNAFAAAAoH7Q6AEAABiKRg8AAMBQNHoAAACGotEDAAAwFI0eAACAoWj0AAAADEWjBwAAYCgaPQAAAEPR6AFBqqCgQJZlaeXKlYEupVZKSkp0yy236LzzzpNlWZo3b169nWvZsmW68MIL1aRJE7Vo0aLezlMfBg4cqIEDBwa6DACGoNEDfkZeXp4sy1J0dLT27dtX7f2BAwfq4osvDkBlwWfChAl66623NHnyZC1btkxDhgw567GWZWns2LF1Os/27ds1atQode7cWS+//LIWLlxY15LPav369Xr88cd15MgR18cGADdFBLoAIBhUVFRo1qxZev755wNdStB65513NGzYME2cOLFez1NQUCDbtjV//nx16dKlXs6xfv16zZgxQ6NGjQq6xBBAaCHRA2ohJSVFL7/8svbv3x/oUhpceXm5K+McPHiwQZqigwcPShINGACIRg+olSlTpsjr9WrWrFk/e9zevXtlWZby8vKqvWdZlh5//HHf68cff1yWZenLL7/UXXfdpbi4OLVp00ZTp06V4zgqKirSsGHD5PF4lJiYqGeffbbGc3q9Xk2ZMkWJiYlq1qyZrr/+ehUVFVU77uOPP9aQIUMUFxenmJgYDRgwQB9++GGVY87U9Pnnn+uOO+5Qy5YtdcUVV/zsd969e7eGDx+uVq1aKSYmRpdffrlWr17te//M9LfjOFqwYIEsy5JlWT875k+duR7x1Vdf1cyZM3XBBRcoOjpa11xzjb766ivfcR06dND06dMlSW3atKn2m7/55pu68sor1axZM8XGxuq6667TZ599Vu1827dv16233qo2bdqoadOm6t69ux577DHfb/Tb3/5WktSxY0ff99m7d6/v88uXL1efPn3UtGlTtWrVSrfddluNfyYLFy5U586d1bRpU6Wlpen999/363cBgF9CowfUQseOHTVixIh6SfUyMzNl27ZmzZql9PR0Pfnkk5o3b54GDRqkpKQkPf300+rSpYsmTpyo9957r9rnZ86cqdWrV+vRRx/VuHHjtG7dOmVkZOhf//qX75h33nlH/fv3V1lZmaZPn66nnnpKR44c0dVXX62NGzdWG3P48OE6ceKEnnrqKY0ZM+astZeUlOhXv/qV3nrrLf3mN7/RzJkzdfLkSV1//fV6/fXXJUn9+/fXsmXLJEmDBg3SsmXLfK/9NWvWLL3++uuaOHGiJk+erI8++kh33nmn7/158+bpxhtvlCS9+OKLWrZsmW666SZJpxdoXHfddWrevLmefvppTZ06VZ9//rmuuOKKKk3a1q1blZ6ernfeeUdjxozR/PnzdcMNN+hvf/ubJOmmm27S7bffLkmaO3eu7/u0adNG0uk/jxEjRqhr166aM2eOxo8fr/z8fPXv37/KNX2LFi3S/fffr8TERD3zzDPq16/fWZt0AKgzB8BZLVmyxJHk/OMf/3B27drlREREOOPGjfO9P2DAAKdnz56+13v27HEkOUuWLKk2liRn+vTpvtfTp093JDn33Xefb9+pU6ecCy64wLEsy5k1a5Zv/3fffec0bdrUGTlypG/fu+++60hykpKSnLKyMt/+V1991ZHkzJ8/33Ecx7Ft2+nataszePBgx7Zt33EnTpxwOnbs6AwaNKhaTbfffnutfp/x48c7kpz333/ft+/YsWNOx44dnQ4dOjher7fK93/ooYdqNe5Pjz3zXS+66CKnoqLCt3/+/PmOJGfbtm3VvkNpaWmVmlq0aOGMGTOmynmKi4uduLi4Kvv79+/vxMbGOl9//XWVY3/8282ePduR5OzZs6fKMXv37nXCw8OdmTNnVtm/bds2JyIiwre/srLSiY+Pd1JSUqp8n4ULFzqSnAEDBvzSTwQAtUKiB9RSp06ddPfdd2vhwoU6cOCAa+OOHj3a98/h4eFKTU2V4zi69957fftbtGih7t27a/fu3dU+P2LECMXGxvpe33LLLWrbtq3WrFkjSdqyZYt27typO+64Q4cPH9ahQ4d06NAhlZeX65prrtF7770n27arjPnAAw/UqvY1a9YoLS2tyvRu8+bNdd9992nv3r36/PPPa/cj1FJWVpYiIyN9r6+88kpJqvF3+bF169bpyJEjuv32233f/9ChQwoPD1d6erreffddSVJpaanee+893XPPPWrXrl2VMWoz3fzaa6/Jtm3deuutVc6TmJiorl27+s7zySef6ODBg3rggQeqfJ9Ro0YpLi6udj8GANQCq24BP/zud7/TsmXLNGvWLM2fP9+VMX/aUMTFxSk6OlqtW7eutv/w4cPVPt+1a9cqry3LUpcuXXzTkTt37pQkjRw58qw1HD16VC1btvS97tixY61q//rrr5Wenl5t/0UXXeR7383bz/z0tzpT83ffffeznzvzG1x99dU1vu/xeCT9u2Gsa807d+6U4zjV/kzOaNKkiaTTv4tU/c+uSZMm6tSpU53ODQA1odED/NCpUyfdddddWrhwoSZNmlTt/bOlPl6v96xjhoeH12qfJDmOU8tK/+1MWjd79mylpKTUeEzz5s2rvG7atKnf52kIdf1dzvwGy5YtU2JiYrX3IyLc+VehbduyLEtvvvlmjbX+9HcGgPpGowf46Xe/+52WL1+up59+utp7ZxKmn95I90yCUx/OpFVnOI6jr776Sr169ZIkde7cWdLp1CojI8PVc7dv3147duyotn/79u2+9xuDM79BfHz8z/4GZ9K0Tz/99GfHO1tD37lzZzmOo44dO6pbt25n/fyZ32Xnzp1VUsbvv/9ee/bsUe/evX/2/ABQW1yjB/ipc+fOuuuuu/TSSy+puLi4ynsej0etW7eutjr2hRdeqLd6/vSnP+nYsWO+1ytXrtSBAwc0dOhQSVKfPn3UuXNn/c///I+OHz9e7fOlpaV1Pvd//ud/auPGjdqwYYNvX3l5uRYuXKgOHTqoR48edR7bTYMHD5bH49FTTz2l77//vtr7Z36DNm3aqH///lq8eLEKCwurHPPj1LBZs2aSqjf0N910k8LDwzVjxoxqKaPjOL6p99TUVLVp00a5ubmqrKz0HZOXl8fTNgC4ikQPqIPHHntMy5Yt044dO9SzZ88q740ePVqzZs3S6NGjlZqaqvfee09ffvllvdXSqlUrXXHFFcrKylJJSYnmzZunLl26+G6LEhYWpj/+8Y8aOnSoevbsqaysLCUlJWnfvn1699135fF4fLcO8dekSZP05z//WUOHDtW4cePUqlUrLV26VHv27NH//u//Kiyscfx/SY/HoxdffFF33323LrvsMt12221q06aNCgsLtXr1avXr109/+MMfJEnPPfecrrjiCl122WW677771LFjR+3du1erV6/Wli1bJJ1unqXT/zu47bbb1KRJE/36179W586d9eSTT2ry5Mnau3evbrjhBsXGxmrPnj16/fXXdd9992nixIlq0qSJnnzySd1///26+uqrlZmZqT179mjJkiVcowfAVTR6QB106dJFd911l5YuXVrtvWnTpqm0tFQrV67Uq6++qqFDh+rNN99UfHx8vdQyZcoUbd26VTk5OTp27JiuueYavfDCC4qJifEdM3DgQG3YsEG///3v9Yc//EHHjx9XYmKi0tPTdf/999f53AkJCVq/fr0effRRPf/88zp58qR69eqlv/3tb7ruuuvc+HquueOOO3T++edr1qxZmj17tioqKpSUlKQrr7xSWVlZvuN69+6tjz76SFOnTtWLL76okydPqn379rr11lt9x/zHf/yHfv/73ys3N1dr166Vbdvas2ePmjVrpkmTJqlbt26aO3euZsyYIUlKTk7Wtddeq+uvv943xn333Sev16vZs2frt7/9rS655BKtWrVKU6dObbgfBYDxLKcuV3cDAACg0Wsc8yoAAABwHY0eAACAoWj0AAAADEWjBwAAYCgaPQAAAEPR6AEAABiq0d1Hz7Zt7d+/X7GxsWd9zBAAAKh/juPo2LFjOv/88xvNDdDhn0bX6O3fv1/JycmBLgMAAPygqKhIF1xwQaDLQB00ukYvNjZWknSF/lMRahLgagDU1XOfb/jlgwA0aseP2/pV2iHff5sRfBpdo3dmujZCTRRh0egBwSo2lmkewBRcShW8+DcxAACAoWj0AAAADEWjBwAAYCgaPQAAAEPR6AEAABiKRg8AAMBQNHoAAACGotEDAAAwFI0eAACAoWj0AAAADEWjBwAAYCgaPQAAAEPR6AEAABiKRg8AAMBQNHoAAACGotEDAAAwFI0eAACAoWj0AAAADEWjBwAAYCgaPQAAAEPR6AEAABgqItAFAAAA1KeTJ0+qsrLS1TEjIyMVHR3t6pj1gUYPAAAY6+TJk+rYvrmKD3pdHTcxMVF79uxp9M0ejR4AADBWZWWlig96tWdTe3li3blireyYrY59vlZlZSWNHgAAQKB5YsNca/SCCY0eAAAwntex5XXcGytYhF5rCwAAECJI9AAAgPFsObLlTqTn1jgNgUYPAAAYz5YttyZc3Rup/jF1CwAAYCgSPQAAYDyv48jruDPl6tY4DYFEDwAAwFAkegAAwHihuhiDRA8AAMBQJHoAAMB4thx5QzDRo9EDAADGY+oWAAAARiHRAwAAxuP2KgAAADAKiR4AADCe/cPm1ljBgkYPAAAYz+viqlu3xmkITN0CAAAYikQPAAAYz+uc3twaK1iQ6AEAABiKRA8AABiPxRgAAACGsmXJK8u1sYIFU7cAAACGItEDAADGs53Tm1tjBQsSPQAAAEOR6AEAAON5XbxGz61xGgKJHgAAgKFI9AAAgPFCNdGj0QMAAMazHUu249LtVVwapyEwdQsAAGAoEj0AAGC8UJ26JdEDAAAwFIkeAAAwnldh8rqUb3ldGaVh0OgBAADjOS4uxnBYjAEAAIBAI9EDAADGYzEGAAAAjEKiBwAAjOd1wuR1XFqM4bgyTIOg0QMAAMazZcl2aSLTVvB0ekzdAgAAGIpEDwAAGI/FGAAAADAKiR4AADCeu4sxuEYPAAAAAUaiBwAAjHd61a0719a5NU5DoNEDAADGsxUmL7dXAQAAgClI9AAAgPFYjAEAAACjkOgBAADj2QoLyUeg0egBAADjeR1LXselJ2O4NE5DYOoWAADAUCR6AADAeF4Xb6/iDaKpWxI9AAAAQ5HoAQAA49lOmGyXbq9ic3sVAAAABBqJHgAAMF6oXqNHowcAAIxny73botiujNIwmLoFAAAwFIkeAAAwnrtPxgienCx4KgUAAIBfSPQAAIDxvE6YvC7dXsWtcRoCjR4AADCeLUu23FqMwbNuAQAAEGAkegAAwHihOnUbPJUCAADALyR6AADAeO4+GSN4crLgqRQAAKCObMdydauLBQsWqEOHDoqOjlZ6ero2btx41mMHDhwoy7Kqbdddd51f56TRAwAAqGcrVqxQdna2pk+frs2bN6t3794aPHiwDh48WOPxr732mg4cOODbPv30U4WHh2v48OF+nZdGDwAAGM/+YerWja0uT8aYM2eOxowZo6ysLPXo0UO5ubmKiYnR4sWLazy+VatWSkxM9G3r1q1TTEwMjR4AAEBDKCsrq7JVVFTUeFxlZaU2bdqkjIwM376wsDBlZGRow4YNtTrXokWLdNttt6lZs2Z+1UijBwAAjGc7Ya5ukpScnKy4uDjflpOTU+O5Dx06JK/Xq4SEhCr7ExISVFxc/Iu1b9y4UZ9++qlGjx7t9/dm1S0AAEAdFBUVyePx+F5HRUXVy3kWLVqkSy65RGlpaX5/lkYPAAAYzytLXpceXXZmHI/HU6XRO5vWrVsrPDxcJSUlVfaXlJQoMTHxZz9bXl6uV155RU888USdamXqFgAAGK8+pm5rKzIyUn369FF+fv6/67Ft5efnq2/fvj/72b/85S+qqKjQXXfdVafvTaIHAABQz7KzszVy5EilpqYqLS1N8+bNU3l5ubKysiRJI0aMUFJSUrXr/BYtWqQbbrhB5513Xp3OS6MHAACM55VcnLr1X2ZmpkpLSzVt2jQVFxcrJSVFa9eu9S3QKCwsVFhY1aRwx44d+uCDD/T222/XuVYaPQAAgAYwduxYjR07tsb3CgoKqu3r3r27HMc5p3PS6AEAAOPV5dq6nxsrWNDoAQAA43mdMHldatDcGqchBE+lAAAA8AuJHgAAMJ4jS7ZLizEcl8ZpCCR6AAAAhiLRAwAAxuMaPQAAABiFRA8AABjPdizZjjvX1rk1TkOg0QMAAMbzKkxelyYy3RqnIQRPpQAAAPALiR4AADBeqE7dkugBAAAYyu9Gr7S0VA8++KDatWunqKgoJSYmavDgwfrwww8lSQsXLtTAgQPl8XhkWZaOHDnids0AAAB+sRXm6hYs/J66vfnmm1VZWamlS5eqU6dOKikpUX5+vg4fPixJOnHihIYMGaIhQ4Zo8uTJrhcMAADgL69jyevSlKtb4zQEvxq9I0eO6P3331dBQYEGDBggSWrfvr3S0tJ8x4wfP16SVFBQ4FqRAAAA8J9f2WPz5s3VvHlz/fWvf1VFRUV91QQAAOCqM4sx3NqChV+NXkREhPLy8rR06VK1aNFC/fr105QpU7R169Y6F1BRUaGysrIqGwAAAM6d31cT3nzzzdq/f79WrVqlIUOGqKCgQJdddpny8vLqVEBOTo7i4uJ8W3Jycp3GAQAAOBvHCZPt0uaY/qzb6OhoDRo0SFOnTtX69es1atQoTZ8+vU4FTJ48WUePHvVtRUVFdRoHAADgbLyyXN2ChSstaY8ePVReXl6nz0ZFRcnj8VTZAAAAcO78WnV7+PBhDR8+XPfcc4969eql2NhYffLJJ3rmmWc0bNgwSVJxcbGKi4v11VdfSZK2bdum2NhYtWvXTq1atXL/GwAAAPwC23HviRa248owDcKvRq958+ZKT0/X3LlztWvXLn3//fdKTk7WmDFjNGXKFElSbm6uZsyY4ftM//79JUlLlizRqFGj3KscAAAAP8tyHKdR9aVlZWWKi4vTQA1ThNUk0OUAqKOXCz8IdAkAztGxY7Z69Tioo0ePBu2lVWf6ipHv3qbI5pGujFl5vFJLr3olKH6X4Fk2AgAAAL/4/Qg0AACAYGPLku3Salm3xmkINHoAAMB4ofqsW6ZuAQAADEWiBwAAjHfmqRZujRUsgqdSAAAA+IVEDwAAGM+W5d4Nk1mMAQAA0Hg4Lq66dYKo0WPqFgAAwFAkegAAwHi24+LULbdXAQAAQKCR6AEAAOOF6u1VaPQAAIDxmLoFAACAUUj0AACA8WwXb68STPfRI9EDAAAwFIkeAAAwHtfoAQAAwCgkegAAwHihmujR6AEAAOOFaqPH1C0AAIChSPQAAIDxSPQAAABgFBI9AABgPEfu3ejYcWWUhkGjBwAAjMfULQAAAIxCogcAAIxHogcAAACjkOgBAADjkegBAADAKCR6AADAeKGa6NHoAQAA4zmOJcelBs2tcRoCU7cAAACGItEDAADGs2W59mQMt8ZpCCR6AAAAhiLRAwAAxmMxBgAAgKFYjAEAAACjkOgBAADjherULYkeAACAoWj0AACA8c5co+fWVhcLFixQhw4dFB0drfT0dG3cuPFnjz9y5IgeeughtW3bVlFRUerWrZvWrFnj1zmZugUAAMZzXJy6rUujt2LFCmVnZys3N1fp6emaN2+eBg8erB07dig+Pr7a8ZWVlRo0aJDi4+O1cuVKJSUl6euvv1aLFi38Oi+NHgAAQD2bM2eOxowZo6ysLElSbm6uVq9ercWLF2vSpEnVjl+8eLG+/fZbrV+/Xk2aNJEkdejQwe/zMnULAACM50hyHJe2H8YsKyurslVUVNR47srKSm3atEkZGRm+fWFhYcrIyNCGDRtq/MyqVavUt29fPfTQQ0pISNDFF1+sp556Sl6v16/vTaMHAABQB8nJyYqLi/NtOTk5NR536NAheb1eJSQkVNmfkJCg4uLiGj+ze/durVy5Ul6vV2vWrNHUqVP17LPP6sknn/SrRqZuAQCA8WxZslx+1m1RUZE8Ho9vf1RUlCvjS5Jt24qPj9fChQsVHh6uPn36aN++fZo9e7amT59e63Fo9AAAAOrA4/FUafTOpnXr1goPD1dJSUmV/SUlJUpMTKzxM23btlWTJk0UHh7u23fRRRepuLhYlZWVioyMrFWNTN0CAADjBfL2KpGRkerTp4/y8/N9+2zbVn5+vvr27VvjZ/r166evvvpKtm379n355Zdq27ZtrZs8iUYPAACEgDNPxnBr81d2drZefvllLV26VF988YUefPBBlZeX+1bhjhgxQpMnT/Yd/+CDD+rbb7/VI488oi+//FKrV6/WU089pYceesiv8zJ1CwAAUM8yMzNVWlqqadOmqbi4WCkpKVq7dq1vgUZhYaHCwv6dvyUnJ+utt97ShAkT1KtXLyUlJemRRx7Ro48+6td5afQAAIDxztwaxa2x6mLs2LEaO3Zsje8VFBRU29e3b1999NFHdTvZD5i6BQAAMBSJHgAAMN65PKO2prGCBY0eAAAwXqg2ekzdAgAAGIpEDwAAGM92LFkuJXF1ub1KoJDoAQAAGIpEDwAAGK8x3F4lEGj0AACA8U43em4txnBlmAbB1C0AAIChSPQAAIDxuL0KAAAAjEKiBwAAjOf8sLk1VrAg0QMAADAUiR4AADBeqF6jR6MHAADMF6Jzt0zdAgAAGIpEDwAAmM/FqVsF0dQtiR4AAIChSPQAAIDxeNYtAACAoUJ11S1TtwAAAIYi0QMAAOZzLPcWUZDoAQAAINBI9AAAgPFCdTEGiR4AAIChSPQAAID5QvQRaDR6AADAeNxeBQAAAEYh0QMAAKEhiKZc3UKiBwAAYCgSPQAAYLxQvUaPRg8AAJgvRFfdMnULAABgKBI9AAAQAqwfNrfGCg4kegAAAIYi0QMAAOYL0Wv0aPQAAID5QrTRY+oWAADAUCR6AADAfI51enNrrCBBogcAAGAoEj0AAGA8xzm9uTVWsCDRAwAAMBSJHgAAMF+Irrql0QMAAOZjMQYAAABMQqIHAACMZzmnN7fGChYkegAAAIYi0QMAAOZjMQYAAIChWIwBAAAAk5DoAQAA84Xo1C2JHgAAgKFI9AAAgPlI9AAAAGASEj0AAGC+EE30aPQAAID5uL0KAAAATEKiBwAAjMezbgEAAGAUEj0AAGC+EF2MQaIHAABgKBo9AAAAQ9HoAQAA41n694KMc97qWMOCBQvUoUMHRUdHKz09XRs3bjzrsXl5ebIsq8oWHR3t9zkb7TV64Z5YhVuRgS4DQB21i2ge6BIAnKOyCFvSwUCXYYQVK1YoOztbubm5Sk9P17x58zR48GDt2LFD8fHxNX7G4/Fox44dvteW5X+LSaIHAADMd+aGyW5tfpozZ47GjBmjrKws9ejRQ7m5uYqJidHixYvP+hnLspSYmOjbEhIS/D4vjR4AADCf4/Lmh8rKSm3atEkZGRm+fWFhYcrIyNCGDRvO+rnjx4+rffv2Sk5O1rBhw/TZZ5/5d2LR6AEAANRJWVlZla2ioqLG4w4dOiSv11stkUtISFBxcXGNn+nevbsWL16sN954Q8uXL5dt2/rVr36lb775xq8aafQAAID56iHRS05OVlxcnG/Lyclxrdy+fftqxIgRSklJ0YABA/Taa6+pTZs2eumll/wap9EuxgAAAGjMioqK5PF4fK+joqJqPK5169YKDw9XSUlJlf0lJSVKTEys1bmaNGmiSy+9VF999ZVfNZLoAQAA47l2a5UfPTPX4/FU2c7W6EVGRqpPnz7Kz8/37bNtW/n5+erbt2+t6vd6vdq2bZvatm3r1/cm0QMAAKhn2dnZGjlypFJTU5WWlqZ58+apvLxcWVlZkqQRI0YoKSnJN/37xBNP6PLLL1eXLl105MgRzZ49W19//bVGjx7t13lp9AAAgPkC/KzbzMxMlZaWatq0aSouLlZKSorWrl3rW6BRWFiosLB/T7R+9913GjNmjIqLi9WyZUv16dNH69evV48ePfw6r+U4TqN6NG9ZWZni4uJ0jecuRXDDZCBordn+XqBLAHCOyo7Zatltt44ePVrlWrRgcqav6PD7mQqrw5MlamKfPKm9Ux8Lit+Fa/QAAAAMxdQtAAAw3o8XUbgxVrAg0QMAADAUiR4AADBfHZ9Re9axggSNHgAAMF+AV90GClO3AAAAhiLRAwAAxmMxBgAAAIxCogcAAMwXotfo0egBAADzuTh1G0yNHlO3AAAAhiLRAwAA5gvRqVsSPQAAAEOR6AEAAPOR6AEAAMAkJHoAAMB43DAZAAAARqHRAwAAMBRTtwAAwHwsxgAAAIBJSPQAAIDxQnUxBo0eAAAIDUHUoLmFqVsAAABDkegBAADzsRgDAAAAJiHRAwAAxgvVxRgkegAAAIYi0QMAAOYL0Wv0aPQAAIDxmLoFAACAUUj0AACA+UJ06pZEDwAAwFAkegAAwHwhmujR6AEAAOOxGAMAAABGIdEDAADmC9GpWxI9AAAAQ5HoAQAA84VookejBwAAjMdiDAAAABiFRA8AAJgvRKduSfQAAAAMRaIHAACMxzV6AAAAMAqJHgAAMF+IXqNHowcAAMwXoo0eU7cAAACGItEDAADGs37Y3BorWJDoAQAAGIpEDwAAmC9Er9Gj0QMAAMbjPnoAAAAwCokeAAAwX4hO3ZLoAQAAGIpEDwAAhIYgSuLcQqIHAABgKBI9AABgvFBddUujBwAAzMdiDAAAAJiERA8AABgvVKduSfQAAAAMRaIHAADMF6LX6NHoAQAA4zF1CwAAgHqzYMECdejQQdHR0UpPT9fGjRtr9blXXnlFlmXphhtu8PucNHoAAMB8jsubn1asWKHs7GxNnz5dmzdvVu/evTV48GAdPHjwZz+3d+9eTZw4UVdeeaX/JxWNHgAAQL2bM2eOxowZo6ysLPXo0UO5ubmKiYnR4sWLz/oZr9erO++8UzNmzFCnTp3qdF4aPQAAYL4AJnqVlZXatGmTMjIyfPvCwsKUkZGhDRs2nPVzTzzxhOLj43Xvvff6d8IfYTEGAAAwXn0sxigrK6uyPyoqSlFRUdWOP3TokLxerxISEqrsT0hI0Pbt22s8xwcffKBFixZpy5Yt51QriR4AAEAdJCcnKy4uzrfl5OS4Mu6xY8d099136+WXX1br1q3PaSwSPQAAYL56uI9eUVGRPB6Pb3dNaZ4ktW7dWuHh4SopKamyv6SkRImJidWO37Vrl/bu3atf//rXvn22bUuSIiIitGPHDnXu3LlWpZLoAQAA1IHH46myna3Ri4yMVJ8+fZSfn+/bZ9u28vPz1bdv32rHX3jhhdq2bZu2bNni266//npdddVV2rJli5KTk2tdI4keAAAwnuU4shx3Ir26jJOdna2RI0cqNTVVaWlpmjdvnsrLy5WVlSVJGjFihJKSkpSTk6Po6GhdfPHFVT7fokULSaq2/5fQ6AEAANSzzMxMlZaWatq0aSouLlZKSorWrl3rW6BRWFiosDD3J1pp9AAAgPkawbNux44dq7Fjx9b4XkFBwc9+Ni8vr07npNEDAADG41m3AAAAMAqJHgAAMF8jmLoNBBI9AAAAQ5HoAQAA44XqNXo0egAAwHxM3QIAAMAkJHoAAMB4oTp1S6IHAABgKBI9AABgvhC9Ro9GDwAAhIRgmnJ1C1O3AAAAhiLRAwAA5nOc05tbYwUJEj0AAABDkegBAADjcXsVAAAAGIVEDwAAmI/bqwAAAJjJsk9vbo0VLJi6BQAAMBSJHgAAMF+ITt2S6AEAABiKRA8AABgvVG+vQqMHAADMx5MxAAAAYBISPQAAYLxQnbol0QMAADAUiR4AADAft1cBAACASUj0AACA8UL1Gj0aPQAAYD5urwIAAACTkOgBAADjherULYkeAACAofxu9EpLS/Xggw+qXbt2ioqKUmJiogYPHqyZM2fKsqyf3QoKCurhKwAAAPwCx+UtSPg9dXvzzTersrJSS5cuVadOnVRSUqL8/Hz17NlTBw4c8B33yCOPqKysTEuWLPHta9WqlTtVAwAA+CFUp279avSOHDmi999/XwUFBRowYIAkqX379kpLS6t2bNOmTVVRUaHExER3KgUAAIBf/Jq6bd68uZo3b66//vWvqqioqK+aAAAA3GU77m5Bwq9GLyIiQnl5eVq6dKlatGihfv36acqUKdq6dWudC6ioqFBZWVmVDQAAAOfO78UYN998s/bv369Vq1ZpyJAhKigo0GWXXaa8vLw6FZCTk6O4uDjflpycXKdxAAAAzipEF2PU6fYq0dHRGjRokKZOnar169dr1KhRmj59ep0KmDx5so4ePerbioqK6jQOAADA2Vj694KMc94C/WX84Mp99Hr06KHy8vI6fTYqKkoej6fKBgAAgHPn16rbw4cPa/jw4brnnnvUq1cvxcbG6pNPPtEzzzyjYcOG1VeNAAAA5yZEn3XrV6PXvHlzpaena+7cudq1a5e+//57JScna8yYMZoyZUp91QgAAIA68KvRi4qKUk5OjnJycn7x2LouzgAAAHBbqN4wmWfdAgAAGMrvR6ABAAAEHTdvixJEiR6NHgAAMJ7lOLJcWkTh1jgNgalbAAAAQ5HoAQAA89k/bG6NFSRI9AAAAAxFogcAAIwXqtfo0egBAADzheiqW6ZuAQAADEWiBwAAzBeiz7ol0QMAADAUiR4AADAez7oFAACAUUj0AACA+UL0Gj0aPQAAYDzLPr25NVawYOoWAADAUCR6AADAfCE6dUuiBwAAYCgSPQAAYD4egQYAAGAmy3Fc3epiwYIF6tChg6Kjo5Wenq6NGzee9djXXntNqampatGihZo1a6aUlBQtW7bM73PS6AEAANSzFStWKDs7W9OnT9fmzZvVu3dvDR48WAcPHqzx+FatWumxxx7Thg0btHXrVmVlZSkrK0tvvfWWX+el0QMAAOY7sxjDrc1Pc+bM0ZgxY5SVlaUePXooNzdXMTExWrx4cY3HDxw4UDfeeKMuuugide7cWY888oh69eqlDz74wK/z0ugBAADUo8rKSm3atEkZGRm+fWFhYcrIyNCGDRt+8fOO4yg/P187duxQ//79/To3izEAAID5HElu3ej4h0CvrKysyu6oqChFRUVVO/zQoUPyer1KSEiosj8hIUHbt28/62mOHj2qpKQkVVRUKDw8XC+88IIGDRrkV6kkegAAwHj1sRgjOTlZcXFxvi0nJ8fVmmNjY7Vlyxb94x//0MyZM5Wdna2CggK/xiDRAwAAqIOioiJ5PB7f65rSPElq3bq1wsPDVVJSUmV/SUmJEhMTzzp+WFiYunTpIklKSUnRF198oZycHA0cOLDWNZLoAQAA8zlycTHG6SE9Hk+V7WyNXmRkpPr06aP8/HzfPtu2lZ+fr759+9b6K9i2rYqKCr++NokeAABAPcvOztbIkSOVmpqqtLQ0zZs3T+Xl5crKypIkjRgxQklJSb7p35ycHKWmpqpz586qqKjQmjVrtGzZMr344ot+nZdGDwAAmC/Az7rNzMxUaWmppk2bpuLiYqWkpGjt2rW+BRqFhYUKC/v3RGt5ebl+85vf6JtvvlHTpk114YUXavny5crMzPTrvJbjNK4n85aVlSkuLk7XeO5ShBUZ6HIA1NGa7e8FugQA56jsmK2W3Xbr6NGjVa5FCyZn+oqrez+qiPCap1b9dcpboXf+7+mg+F1I9AAAgPlsSZaLYwUJGj0AAGC8c3lGbU1jBQtW3QIAABiKRA8AAJgvwIsxAoVEDwAAwFAkegAAwHwhmujR6AEAAPOFaKPH1C0AAIChSPQAAID5QvQ+eiR6AAAAhiLRAwAAxgvVGybT6AEAAPOxGAMAAAAmIdEDAADmsx3JcimJs0n0AAAAEGAkegAAwHxcowcAAACTkOgBAIAQ4GKip+BJ9Gj0AACA+Zi6BQAAgElI9AAAgPlsR65NuXJ7FQAAAAQaiR4AADCfY5/e3BorSNDoAQAA87EYAwAAACYh0QMAAOZjMQYAAABMQqIHAADMxzV6AAAAMAmJHgAAMJ8jFxM9d4ZpCDR6AADAfEzdAgAAwCQkegAAwHy2LcmlJ1rYwfNkDBI9AAAAQ5HoAQAA84XoNXo0egAAwHwh2ugxdQsAAGAoEj0AAGA+nnULAAAAk5DoAQAA4zmOLcdx57Yobo3TEGj0AACA+RzHvSlXFmMAAAAg0Ej0AACA+RwXF2OQ6AEAACDQSPQAAID5bFuyXFpEEUSLMUj0AAAADEWiBwAAzBei1+jR6AEAAOM5ti3HpanbYLqPHlO3AAAAhiLRAwAA5gvRqVsSPQAAAEOR6AEAAPPZjmSFXqJHowcAAMznOJLcuo9e8DR6TN0CAAAYikQPAAAYz7EdOS5N3TokegAAAAg0Ej0AAGA+x5Z71+hxw2QAAAAEGIkeAAAwXqheo0ejBwAAzBeiU7c0egAAwHin9L1rT0A7pe/dGagB0OgBAABjRUZGKjExUR8Ur3F13MTEREVGRro6Zn2g0QMAAMaKjo7Wnj17VFlZ6eq4kZGRio6OdnXM+kCjBwAAjBYdHR0UTVl94PYqAAAAhqLRAwAAMBSNHgAAgKFo9AAAAAxFowcAAGAoGj0AAABD0egBAAAYikYPAADAUDR6AAAAhqLRAwAAMBSNHgAAgKFo9AAAAAxFowcAAGAoGj0AAABD0egBAAAYikYPAADAUDR6AAAAhqLRAwAAMBSNHgAAgKFo9AAAAAwVEegCfspxHEnSKacywJUAOBdlx+xAlwDgHJUdP/33+Mx/mxF8Gl2jd+zYMUnS34+9GuBKAJyLlt0CXQEAtxw7dkxxcXGBLgN1YDmNrE23bVv79+9XbGysLMsKdDmoJ2VlZUpOTlZRUZE8Hk+gywFQB/w9Np/jODp27JjOP/98hYVxtVcwanSJXlhYmC644IJAl4EG4vF4+A8EEOT4e2w2krzgRnsOAABgKBo9AAAAQ9HoISCioqI0ffp0RUVFBboUAHXE32Og8Wt0izEAAADgDhI9AAAAQ9HoAQAAGIpGDwAAwFA0egAAAIai0QMAADAUjR4AAIChGt0j0GC+EydOqLCwUJWVlVX29+rVK0AVAfDHypUr9eqrr9b493jz5s0BqgpATUj00GBKS0v1X//1X4qNjVXPnj116aWXVtkANH7PPfecsrKylJCQoH/+859KS0vTeeedp927d2vo0KGBLg/AT9DoocGMHz9eR44c0ccff6ymTZtq7dq1Wrp0qbp27apVq1YFujwAtfDCCy9o4cKFev755xUZGan/9//+n9atW6dx48bp6NGjgS4PwE/wZAw0mLZt2+qNN95QWlqaPB6PPvnkE3Xr1k2rVq3SM888ow8++CDQJQL4BTExMfriiy/Uvn17xcfHa926derdu7d27typyy+/XIcPHw50iQB+hEQPDaa8vFzx8fGSpJYtW6q0tFSSdMkll3BdDxAkEhMT9e2330qS2rVrp48++kiStGfPHpEbAI0PjR4aTPfu3bVjxw5JUu/evfXSSy9p3759ys3NVdu2bQNcHYDauPrqq32XWmRlZWnChAkaNGiQMjMzdeONNwa4OgA/xdQtGszy5ct16tQpjRo1Sps2bdKQIUP07bffKjIyUnl5ecrMzAx0iQB+gW3bsm1bERGnb9rwyiuvaP369eratavuv/9+RUZGBrhCAD9Go4eAOXHihLZv36527dqpdevWgS4HAADj0OihwTzxxBOaOHGiYmJiquz/17/+pdmzZ2vatGkBqgyAP06ePKmtW7fq4MGDsm27ynvXX399gKoCUBMaPTSY8PBwHThwwLcg44zDhw8rPj5eXq83QJUBqK21a9dqxIgROnToULX3LMvi7zHQyLAYAw3GcRxZllVt///93/+pVatWAagIgL8efvhhDR8+XAcOHPBdr3dmo8kDGh8egYZ617JlS1mWJcuy1K1btyrNntfr1fHjx/XAAw8EsEIAtVVSUqLs7GwlJCQEuhQAtUCjh3o3b948OY6je+65RzNmzFBcXJzvvcjISHXo0EF9+/YNYIUAauuWW25RQUGBOnfuHOhSANQC1+ihwfz9739Xv379fLdlABB8Tpw4oeHDh6tNmza65JJL1KRJkyrvjxs3LkCVAagJjR4a1K5du7RkyRLt2rVL8+fPV3x8vN588021a9dOPXv2DHR5AH7BokWL9MADDyg6OlrnnXdelUsxLMvS7t27A1gdgJ+i0UOD+fvf/66hQ4eqX79+eu+99/TFF1+oU6dOmjVrlj755BOtXLky0CUC+AWJiYkaN26cJk2apLAw1vMBjR1/S9FgJk2apCeffFLr1q2rcvf8q6++2ve8TACNW2VlpTIzM2nygCDB31Q0mG3bttX4LMz4+Pga78kFoPEZOXKkVqxYEegyANQSV8WjwbRo0UIHDhxQx44dq+z/5z//qaSkpABVBcAfXq9XzzzzjN566y316tWr2mKMOXPmBKgyADWh0UODue222/Too4/qL3/5iyzLkm3b+vDDDzVx4kSNGDEi0OUBqIVt27bp0ksvlSR9+umnVd6r6YboAAKLxRhoMJWVlXrooYeUl5cnr9eriIgInTp1Snfeeafy8vIUHh4e6BIBADAKjR4aXFFRkbZt26by8nJdeuml6tKlS6BLAgDASEzdokEtWrRIc+fO1c6dOyVJXbt21fjx4zV69OgAVwagNk6ePKnnn39e7777rg4ePCjbtqu8v3nz5gBVBqAmNHpoMNOmTdOcOXP08MMP+x55tmHDBk2YMEGFhYV64oknAlwhgF9y77336u2339Ytt9yitLQ0rssDGjmmbtFg2rRpo+eee0633357lf1//vOf9fDDD3OLFSAIxMXFac2aNerXr1+gSwFQC9xHDw3m+++/V2pqarX9ffr00alTpwJQEQB/JSUlKTY2NtBlAKglGj00mLvvvlsvvvhitf0LFy7UnXfeGYCKAPjr2Wef1aOPPqqvv/460KUAqAWu0UO9ys7O9v2zZVn64x//qLfffluXX365JOnjjz9WYWEh99EDgkRqaqpOnjypTp06KSYmptoNk7/99tsAVQagJlyjh3p11VVX1eo4y7L0zjvv1HM1AM5VRkaGCgsLde+99yohIaHaYoyRI0cGqDIANaHRAwDUWkxMjDZs2KDevXsHuhQAtcA1egCAWrvwwgv1r3/9K9BlAKglGj0AQK3NmjVL//3f/62CggIdPnxYZWVlVTYAjQtTtwCAWgsLO50P/PTaPMdxZFmWvF5vIMoCcBasugUA1Nq7774b6BIA+IFEDwAAwFAkegAAv504cUKFhYWqrKyssr9Xr14BqghATWj0AAC1VlpaqqysLL355ps1vs81ekDjwqpbAECtjR8/XkeOHNHHH3+spk2bau3atVq6dKm6du2qVatWBbo8AD9BogcAqLV33nlHb7zxhlJTUxUWFqb27dtr0KBB8ng8ysnJ0XXXXRfoEgH8CIkeAKDWysvLFR8fL0lq2bKlSktLJUmXXHKJNm/eHMjSANSARg8AUGvdu3fXjh07JEm9e/fWSy+9pH379ik3N1dt27YNcHUAforbqwAAam358uU6deqURo0apU2bNmnIkCE6fPiwIiMjtXTpUmVmZga6RAA/QqMHAKizEydOaPv27WrXrp1at24d6HIA/ASLMQAAtZadnV3jfsuyFB0drS5dumjYsGFq1apVA1cGoCYkegCAWrvqqqu0efNmeb1ede/eXZL05ZdfKjw8XBdeeKF27Nghy7L0wQcfqEePHgGuFgCLMQAAtTZs2DBlZGRo//792rRpkzZt2qRvvvlGgwYN0u233659+/apf//+mjBhQqBLBSASPQCAH5KSkrRu3bpqad1nn32ma6+9Vvv27dPmzZt17bXX6tChQwGqEsAZJHoAgFo7evSoDh48WG1/aWmpysrKJEktWrSo9gxcAIFBowcAqLVhw4bpnnvu0euvv65vvvlG33zzjV5//XXde++9uuGGGyRJGzduVLdu3QJbKABJTN0CAPxw/PhxTZgwQX/605906tQpSVJERIRGjhypuXPnqlmzZtqyZYskKSUlJXCFApBEowcAqIPjx49r9+7dkqROnTqpefPmAa4IQE1o9AAAAAzFNXoAAACGotEDAAAwFI0eAACAoWj0AAAADEWjBwAAYCgaPQAAAEPR6AEAABiKRg8AAMBQ/x/Dbgj15JsN7AAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 107,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"calibrator.specify().execute().analyze()"
]
},
{
"cell_type": "markdown",
"id": "180bf22e-7ab4-492a-b336-93a9b239eabc",
"metadata": {},
"source": [
"According to the indices above, the variation in the end-of-simulation number of infected is more greatly influenced by *γ* over *β*. Hence, we should generally prioritise the former over the latter for calibration. \n",
"\n",
"## Sobol Sensitivity Analysis via Polynomial Chaos\n",
"\n",
"One advantage of `calisim` is the ability to compare different libraries, algorithms, and calibration methods. Let's construct Sobol sensitivity indices by fitting a polynomial chaos expansion surrogate model using the [OpenTurns engine](https://openturns.github.io/www/). \n",
"\n",
"We can reuse both the parameter specification and wrapper function defined above. But we will need to change the calibration specification."
]
},
{
"cell_type": "code",
"execution_count": 108,
"id": "a4bad659-ba2e-4c5d-a6e4-641038bbf6e8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 108,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA+kAAARTCAYAAAAXy2LiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbgUlEQVR4nO3deZyVdf3//+ewi7K4AQoopuCSILib+XEDya3Qj3342AKaaZmkhluWimZKmppmll8tkUrTXFJLVJDESjHccCnXlDAF1FyGJWGcmd8f/JhPE6gMDpw3zv1+u3GLc533dc5rZm7Yecx1netU1dfX1wcAAACouFaVHgAAAABYTKQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkA8Aq0KdPnxxwwAHN9ngzZsxIVVVVrr766ibve9hhh6VPnz7NNst7PWZVVVXOPPPMZn0eAPioE+kAsAxPPPFEDjnkkGy88cbp0KFDevbsmSFDhuTSSy+t9Gjv6cwzz2z2+AYAVq02lR4AAEpz//33Z88998xGG22UI488Mj169MhLL72UBx54IJdcckm+/vWvV3rED+XKK69MXV3dSn+ef/3rX2nTxksNAGgK/88JAP/hnHPOSZcuXfLggw+ma9euje579dVXKzNUM2rbtu0qeZ4OHTqskucBgI8Sp7sDwH/429/+lo9//ONLBXqSdOvWrdHtd999N2effXY23XTTtG/fPn369Mm3vvWtLFy4cJmPPXHixAwcODAdOnTIVlttlZtvvnmpNS+88EI++9nPZp111knHjh2z88475/bbb2+Wry1Z+v3jS97ffsEFF+SKK65o+Fp22GGHPPjgg0vtf8stt2TrrbdOhw4dsvXWW+c3v/nNMp9nWe9Jf/nll3PEEUdkww03TPv27bPJJpvk6KOPzqJFixrWvPXWWzn++OPTu3fvtG/fPptttlnOO++8pY7+X3fdddluu+3SqVOndO7cOf37988ll1yy4t8YACiAI+kA8B823njjTJ06NU8++WS23nrr91375S9/OePHj88hhxySE044IX/+858zduzYPPXUU0vF63PPPZfhw4fnq1/9akaOHJlx48bls5/9bO68884MGTIkSTJnzpx84hOfyIIFC3Lsscdm3XXXzfjx4/PpT386N954Yw466KCV9nVfe+21mTt3br7yla+kqqoq559/fg4++OC88MILDUffJ06cmP/+7//OVlttlbFjx+af//xnDj/88PTq1esDH/+VV17JjjvumLfeeitHHXVUtthii7z88su58cYbs2DBgrRr1y4LFizI7rvvnpdffjlf+cpXstFGG+X+++/PqaeemlmzZuXiiy9OkkyaNCmHHnpo9t5775x33nlJkqeeeir33XdfjjvuuJX2PQKAla4eAGhk4sSJ9a1bt65v3bp1/S677FJ/8skn19911131ixYtarRu+vTp9Unqv/zlLzfafuKJJ9Ynqf/973/fsG3jjTeuT1J/0003NWx7++236zfYYIP6QYMGNWw7/vjj65PU//GPf2zYNnfu3PpNNtmkvk+fPvW1tbX19fX19S+++GJ9kvpx48Y1+esbOXJk/cYbb9xwe8ljrbvuuvVvvPFGw/Zbb721Pkn9b3/724ZtAwcOrN9ggw3q33rrrYZtEydOrE/S6DHr6+vrk9SPGTOm4faIESPqW7VqVf/ggw8uNVNdXV19fX19/dlnn12/5ppr1j/77LON7v/mN79Z37p16/qZM2fW19fX1x933HH1nTt3rn/33Xeb/PUDQMmc7g4A/2HIkCGZOnVqPv3pT+exxx7L+eefn6FDh6Znz5657bbbGtZNmDAhSTJ69OhG+59wwglJstQp6htuuGGjI+GdO3fOiBEj8uijj2b27NkNj7njjjvmk5/8ZMO6tdZaK0cddVRmzJiRv/71r837xf6b4cOHZ+211264vdtuuyVZfPp9ksyaNSvTp0/PyJEj06VLl4Z1Q4YMyVZbbfW+j11XV5dbbrklBx54YLbffvul7q+qqkqS3HDDDdltt92y9tpr5/XXX2/4M3jw4NTW1uYPf/hDkqRr166ZP39+Jk2a9OG+aAAojEgHgGXYYYcdcvPNN+fNN9/MtGnTcuqpp2bu3Lk55JBDGkL573//e1q1apXNNtus0b49evRI165d8/e//73R9s0226whRpfo169fksXvC1/ymJtvvvlS82y55ZYN968sG220UaPbS4L9zTffbPTcffv2XWrfZc3871577bVUV1d/4NsHnnvuudx5551Zf/31G/0ZPHhwkv+7cN/Xvva19OvXL/vuu2969eqVL33pS7nzzjuX46sEgLJ5TzoAvI927dplhx12yA477JB+/frl8MMPzw033JAxY8Y0rPnP8F5dtW7depnb6+vrV9kMdXV1GTJkSE4++eRl3r/klxrdunXL9OnTc9ddd+WOO+7IHXfckXHjxmXEiBEZP378KpsXAJqbSAeA5bTkNO1Zs2YlWXyBubq6ujz33HMNR7qTxRd/e+utt7Lxxhs32v/5559PfX19o6h/9tlnk6Thausbb7xxnnnmmaWe++mnn264v1KWPPdzzz231H3Lmvnfrb/++uncuXOefPLJ91236aabZt68eQ1Hzt9Pu3btcuCBB+bAAw9MXV1dvva1r+X//b//l9NPP32psxsAYHXhdHcA+A/33HPPMo8eL3kP+pJTu/fbb78kabji+BIXXXRRkmT//fdvtP2VV15pdMX36urq/PznP8/AgQPTo0ePhsecNm1apk6d2rBu/vz5ueKKK9KnT58PfO/3yrTBBhtk4MCBGT9+fN5+++2G7ZMmTfrA98q3atUqw4YNy29/+9s89NBDS92/5Pv9P//zP5k6dWruuuuupda89dZbeffdd5Mk//znP5d6/AEDBiTJe378HQCsDhxJB4D/8PWvfz0LFizIQQcdlC222CKLFi3K/fffn+uvvz59+vTJ4YcfniTZZpttMnLkyFxxxRV56623svvuu2fatGkZP358hg0blj333LPR4/br1y9HHHFEHnzwwXTv3j1XXXVV5syZk3HjxjWs+eY3v5lf/epX2XfffXPsscdmnXXWyfjx4/Piiy/mpptuSqtWlf39+tixY7P//vvnk5/8ZL70pS/ljTfeyKWXXpqPf/zjmTdv3vvue+6552bixInZfffdc9RRR2XLLbfMrFmzcsMNN+RPf/pTunbtmpNOOim33XZbDjjggBx22GHZbrvtMn/+/DzxxBO58cYbM2PGjKy33nr58pe/nDfeeCN77bVXevXqlb///e+59NJLM3DgwEZnNQDA6kakA8B/uOCCC3LDDTdkwoQJueKKK7Jo0aJstNFG+drXvpbTTjstXbt2bVj705/+NB/72Mdy9dVX5ze/+U169OiRU089tdF71pfo27dvLr300px00kl55plnsskmm+T666/P0KFDG9Z07949999/f0455ZRceumleeeddzJgwID89re/XerIfCV86lOfyg033JDTTjstp556ajbddNOMGzcut956a6ZMmfK++/bs2TN//vOfc/rpp+eaa65JdXV1evbsmX333TcdO3ZMknTs2DH33ntvzj333Nxwww35+c9/ns6dO6dfv34566yzGq4q/4UvfCFXXHFFfvzjH+ett95Kjx49Mnz48Jx55pkV/0UGAHwYVfWr8mowAAAAwHvyq2YAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQrSp9ACrWl1dXV555ZV06tQpVVVVlR4HAACAj7j6+vrMnTs3G264YVq1ev9j5S0u0l955ZX07t270mMAAADQwrz00kvp1avX+65pcZHeqVOnJIu/OZ07d67wNDSHmpqaTJw4Mfvss0/atm1b6XEAACrCayIoV3V1dXr37t3Qo++nxUX6klPcO3fuLNI/ImpqatKxY8d07tzZ/yEBAC2W10RQvuV5y7ULxwEAAEAhRDoAAAAUQqQDAABAIVrce9IBAABWpfr6+tTU1OTdd9+t9CisRO3bt0/r1q0/9OOIdAAAgJVk4cKFmTFjRubNm1fpUVjJqqqqstlmm33oC5SLdAAAgJWgrq4uf/3rX9OmTZtssskmad++/XJd3ZvVT11dXWbNmpXnn38+/fr1y1prrbXCjyXSAQAAVoJ33nkndXV12WSTTT5UtLF62GCDDVJdXZ2bbropu+yyS/r167dCj+PCcQAAACtRq1ayqyVY8nP+17/+lbvuuisvvPDCij1Ocw4FAAAALdl6662XuXPnZsaMGSu0v0gHAABgue2xxx45/vjjKz1Gk8yYMSNVVVWZPn36Sn+uqqqqtG3bdoUvFijSAQAAaOSwww5LVVXVUn+ef/753HzzzTn77LM/1ONXVVXllltuaZ5hC/RhLhDownEAAAAs5VOf+lTGjRvXaNv666//gZ8FvmjRorRr125ljlaR515VX5cj6QAAAKWrrU2mTEl+9avF/1tbu9Kfsn379unRo0ejP61bt17qdPc+ffrk7LPPzogRI9K5c+ccddRRWbRoUUaNGpUNNtggHTp0yMYbb5yxY8c2rE+Sgw46KFVVVQ23l+WJJ57IXnvtlTXWWCPrrrtujjrqqEankR922GEZNmxYzjnnnGy44YbZfPPNkyTTpk3LoEGD0qFDh2y//fZ59NFHl3rsJ598Mvvuu2/WWmutdO/ePV/84hfz+uuvN9y/xx57ZNSoUTn++OOz3nrrZejQoR/iu7n8RDoAAEDJbr456dMn2XPP5HOfW/y/ffos3l6ICy64INtss00effTRnH766fnhD3+Y2267Lb/+9a/zzDPP5JprrmmI8QcffDBJMm7cuMyaNavh9n+aP39+hg4dmrXXXjsPPvhgbrjhhtx9990ZNWpUo3WTJ0/OM888k0mTJuV3v/td5s2blwMOOCBbbbVVHn744Zx55pk58cQTG+3z1ltvZa+99sqgQYPy0EMP5c4778ycOXPyP//zP43WjR8/Pu3atct9992Xyy+/vJm+W+/P6e4AAACluvnm5JBDkvr6xttffnnx9htvTA4+eKU89e9+97tGn+++77775oYbbljm2r322isnnHBCw+2ZM2emb9+++eQnP5mqqqpsvPHGDfetv/76SZKuXbumR48e7/n81157bd555538/Oc/z5prrpkk+dGPfpQDDzww5513Xrp3754kWXPNNfPTn/604VT0K664InV1dfnZz36WDh065OMf/3j+8Y9/5Oijj2547B/96EcZNGhQzj333IZtV111VXr37p1nn3224TPO+/btm/PPP3/5vmHNRKQDAACUqLY2Oe64pQM9Wbytqio5/vjkM59JPuB94itizz33zE9+8pOG20tCeVm23377RrcPO+ywDBkyJJtvvnk+9alP5YADDsg+++zTpOd/6qmnss022zR63l133TV1dXV55plnGiK9f//+jd4r/tRTT2XAgAHp0KFDw7Zddtml0WM/9thjueeeexr9EmKJv/3tbw2Rvt122zVp5uYg0gEAAEr0xz8m//jHe99fX5+89NLidXvs0exPv+aaa2azzTZb7rX/btttt82LL76YO+64I3fffXf+53/+J4MHD86NN964UuZsqnnz5jUckf9PG2ywwYd67A9LpAMAAJRo1qzmXbeKde7cOcOHD8/w4cNzyCGH5FOf+lTeeOONrLPOOmnbtm1qP+Did1tuuWWuvvrqzJ8/vyGW77vvvrRq1arhAnHvtd8vfvGLvPPOOw1H0x944IFGa7bddtvcdNNN6dOnT9q0KSuLXTgOAACgRP92RLdZ1q1CF110UX71q1/l6aefzrPPPpsbbrghPXr0SNeuXZMsvsL75MmTM3v27Lz55pvLfIzPf/7z6dChQ0aOHJknn3wy99xzT77+9a/ni1/8YsOp7svyuc99LlVVVTnyyCPz17/+NRMmTMgFF1zQaM0xxxyTN954I4ceemgefPDB/O1vf8tdd92Vww8//AN/ebCyiXQAAIAS7bZb0qvX4veeL0tVVdK79+J1henUqVPOP//8bL/99tlhhx0yY8aMTJgwIa1aLU7QCy+8MJMmTUrv3r0zaNCgZT5Gx44dc9ddd+WNN97IDjvskEMOOSR77713fvSjH73vc6+11lr57W9/myeeeCKDBg3Kt7/97aVOa99www1z3333pba2Nvvss0/69++f448/Pl27dm2YsVKq6uuXdRWCj67q6up06dIlb7/9djp37lzpcWgGNTU1mTBhQvbbb7+0bdu20uMAAFSE10TlWbBgQZ566qlsueWW6dix44o9yJKruyeNLyC3JNxX4tXdaZolP+8ZM2bkxRdfTL9+/fLpT386SdM61JF0AACAUh188OIQ79mz8fZevQT6R1RZ75AHAACgsYMPXvwxa3/84+KLxG2wweJT3FfCx65ReSIdAACgdK1br5SPWaM8TncHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACg4qqqqnLLLbeskufaY489cvzxx6+S52oqkQ4AAECDqqqq9/1z5plnvue+M2bMSFVVVaZPn77K5v2oaVPpAQAAACjHrFmzGv5+/fXX54wzzsgzzzzTsG2ttdaqxFhNVl9fn9ra2rRp0/zZW1tbm6qqqrRq1fzHvR1JBwAAKNRzzyWPPPLef557rvmfs0ePHg1/unTpkqqqqobb3bp1y0UXXZRevXqlffv2GThwYO68886GfTfZZJMkyaBBg1JVVZU99tgjSfLggw9myJAhWW+99dKlS5fsvvvueeSRR5o018KFC3PsscemW7du6dChQz75yU/mwQcfbLh/ypQpqaqqyh133JHtttsu7du3z5/+9KfMnz8/I0aMyFprrZUNNtggF1544TIf+8QTT0zPnj2z5pprZqeddsqUKVMa7r/66qvTtWvX3Hbbbdlqq63Svn37zJw5s0nzLy+RDgAAUKDnnkv69Uu22+69//Trt3JC/b1ccsklufDCC3PBBRfk8ccfz9ChQ/PpT386z/3/Q0ybNi1Jcvfdd2fWrFm5+eabkyRz587NyJEj86c//SkPPPBA+vbtm/322y9z585d7uc++eSTc9NNN2X8+PF55JFHstlmm2Xo0KF54403Gq375je/me9973t56qmnMmDAgJx00km59957c+utt2bixImZMmXKUr8gGDVqVKZOnZrrrrsujz/+eD772c/mU5/6VMPXlSQLFizIeeedl5/+9Kf5y1/+km7duq3Q9/CDON0dAACgQMvbr03o3A/tggsuyCmnnJL//d//TZKcd955ueeee3LxxRfnsssuy/rrr58kWXfdddOjR4+G/fbaa69Gj3PFFVeka9euuffee3PAAQd84PPOnz8/P/nJT3L11Vdn3333TZJceeWVmTRpUn72s5/lpJNOalj7ne98J0OGDEmSzJs3Lz/72c/yy1/+MnvvvXeSZPz48enVq1fD+pkzZ2bcuHGZOXNmNtxwwyTJiSeemDvvvDPjxo3LueeemySpqanJj3/842yzzTZN+6Y1kUgHAADgA1VXV+eVV17Jrrvu2mj7rrvumscee+x9950zZ05OO+20TJkyJa+++mpqa2uzYMGC5T5l/G9/+1tqamoaPXfbtm2z44475qmnnmq0dvvtt2+036JFi7LTTjs1bFtnnXWy+eabN9x+4oknUltbm379+jV6nIULF2bddddtuN2uXbsMGDBgueb9MEQ6AAAAK9XIkSPzz3/+M5dcckk23njjtG/fPrvssksWLVrU7M+15pprNmn9vHnz0rp16zz88MNp3bp1o/v+/SJ5a6yxRqqqqpplxvfjPekAAAB8oM6dO2fDDTfMfffd12j7fffdl6222irJ4qPNyeKrn//nmmOPPTb77bdfPv7xj6d9+/Z5/fXXl/u5N91007Rr167Rc9fU1OTBBx9seO732q9t27b585//3LDtzTffzLPPPttwe9CgQamtrc2rr76azTbbrNGffz9lf1VxJB0AAIDlctJJJ2XMmDHZdNNNM3DgwIwbNy7Tp0/PNddckyTp1q1b1lhjjdx5553p1atXOnTokC5duqRv3775xS9+ke233z7V1dU56aSTssYaayz386655po5+uijc9JJJ2WdddbJRhttlPPPPz8LFizIEUcc8Z77rbXWWjniiCNy0kknZd111023bt3y7W9/u9FHp/Xr1y+f//znM2LEiFx44YUZNGhQXnvttUyePDkDBgzI/vvvv+LfsBUg0gEAAFguxx57bN5+++2ccMIJefXVV7PVVlvltttuS9++fZMkbdq0yQ9/+MN85zvfyRlnnJHddtstU6ZMyc9+9rMcddRR2XbbbdO7d++ce+65OfHEE5v03N/73vdSV1eXL37xi5k7d26233773HXXXVl77bXfd7/vf//7mTdvXg488MB06tQpJ5xwQt5+++1Ga8aNG5fvfve7OeGEE/Lyyy9nvfXWy84777xcF7VrblX19fX1q/xZK6i6ujpdunTJ22+/nc6dO1d6HJpBTU1NJkyYkP322y9t27at9DgAABXhNVF5FixYkKeeeipbbrllOnbs2OT9H3lk8cesfZCHH0623XYFBqRZLfl5z5gxIy+++GL69euXT3/600ma1qHekw4AAFCgTp2adx2rB6e7AwAAFKhv3+TZZ9//c9A7dVq8jo8OkQ4AAFAoAd7yON0dAAAACiHSAQAAoBAiHQAAYCWqq6ur9AisAs31wWkiHQAAYCVo165dkmTevHkVnoRVYeHChUmSd99990M9jgvHAQAArARt2rTJeuutl5dffjlJstZaa6VVK8dJP4rq6ury0ksvZcGCBamtrf1QR9VFOgAAwEqy0UYbJUlDqPPRVVdXl9mzZ6e+vj6LFi1KpxX8AHuRDgAAsJJUVVVl4403Ttu2bTN58uRUV1enc+fOad26daVHoxnV19dn4cKFqa2tzVtvvZWuXbvmYx/72Ao9lkgHAABYyTbccMPsvffeufPOO/PGG28020XGKEtVVVW6du2awYMHZ+ONN16hxxDpAAAAq8CGG26YkSNHZt68eamtra30OKwErVq1ylprrZU2bVY8tUU6AADAKtK6det06dKl0mNQMJcWBAAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHQAAAAohEgHAACAQoh0AAAAKIRIBwAAgEKIdAAAAChERSP9D3/4Qw488MBsuOGGqaqqyi233PKB+0yZMiXbbrtt2rdvn8022yxXX331Sp8TAAAAVoWKRvr8+fOzzTbb5LLLLluu9S+++GL233//7Lnnnpk+fXqOP/74fPnLX85dd921kicFAACAla9NJZ983333zb777rvc6y+//PJssskmufDCC5MkW265Zf70pz/lBz/4QYYOHbqyxgQAAIBVYrV6T/rUqVMzePDgRtuGDh2aqVOnVmgiAAAAaD4VPZLeVLNnz0737t0bbevevXuqq6vzr3/9K2usscZS+yxcuDALFy5suF1dXZ0kqampSU1NzcodmFViyc/RzxMAaMm8JoJyNeXf5WoV6Sti7NixOeuss5baPnHixHTs2LECE7GyTJo0qdIjAABUnNdEUJ4FCxYs99rVKtJ79OiROXPmNNo2Z86cdO7ceZlH0ZPk1FNPzejRoxtuV1dXp3fv3tlnn33SuXPnlTovq0ZNTU0mTZqUIUOGpG3btpUeBwCgIrwmgnItOaN7eaxWkb7LLrtkwoQJjbZNmjQpu+yyy3vu0759+7Rv336p7W3btvUfr48YP1MAAK+JoERN+TdZ0QvHzZs3L9OnT8/06dOTLP6ItenTp2fmzJlJFh8FHzFiRMP6r371q3nhhRdy8skn5+mnn86Pf/zj/PrXv843vvGNSowPAAAAzaqikf7QQw9l0KBBGTRoUJJk9OjRGTRoUM4444wkyaxZsxqCPUk22WST3H777Zk0aVK22WabXHjhhfnpT3/q49cAAAD4SKjo6e577LFH6uvr3/P+q6++epn7PProoytxKgAAAKiM1epz0gEAAOCjTKQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABSiTaUHAAAAVh/PPZfMnfve93fqlPTtu+rmgY8akQ4AACyX555L+vX74HXPPivUYUU53R0AAFgu73cEfUXWAUsT6QAAAFAIkQ4AAACFEOkAAABQiIpH+mWXXZY+ffqkQ4cO2WmnnTJt2rT3XX/xxRdn8803zxprrJHevXvnG9/4Rt55551VNC0AAACsPBWN9Ouvvz6jR4/OmDFj8sgjj2SbbbbJ0KFD8+qrry5z/bXXXptvfvObGTNmTJ566qn87Gc/y/XXX59vfetbq3hyAAAAaH4VjfSLLrooRx55ZA4//PBstdVWufzyy9OxY8dcddVVy1x///33Z9ddd83nPve59OnTJ/vss08OPfTQDzz6DgAAAKuDikX6okWL8vDDD2fw4MH/N0yrVhk8eHCmTp26zH0+8YlP5OGHH26I8hdeeCETJkzIfvvtt0pmBgCAlqxTp+ZdByytTaWe+PXXX09tbW26d+/eaHv37t3z9NNPL3Ofz33uc3n99dfzyU9+MvX19Xn33Xfz1a9+9X1Pd1+4cGEWLlzYcLu6ujpJUlNTk5qammb4Sqi0JT9HP08AoCVbFa+J+vRJ/vKXZN68916z1lqL13lpBv+nKf8uKxbpK2LKlCk599xz8+Mf/zg77bRTnn/++Rx33HE5++yzc/rppy9zn7Fjx+ass85aavvEiRPTsWPHlT0yq9CkSZMqPQIAQMWV8JrouecqPQGUZcGCBcu9tqq+vr5+Jc7ynhYtWpSOHTvmxhtvzLBhwxq2jxw5Mm+99VZuvfXWpfbZbbfdsvPOO+f73/9+w7Zf/vKXOeqoozJv3ry0arX02fvLOpLeu3fvvP766+ncuXPzflFURE1NTSZNmpQhQ4akbdu2lR4HAKAivCaCclVXV2e99dbL22+//YEdWrEj6e3atct2222XyZMnN0R6XV1dJk+enFGjRi1znwULFiwV4q1bt06SvNfvGtq3b5/27dsvtb1t27b+4/UR42cKANCyXhP98pe/zHe+85106tQpn/zkJ/OXv/wl5513Xo499tgsWLAg9fX1Of/887PPPvskSdq0aZNvfetbufXWW9O6deuMGzcup59+ep566ql8+tOfzoUXXpgk6dOnTz73uc/lrrvuyrx58/LLX/4yF198ccMnUv3qV79KVVVVzj333PzmN7/JokWL0rt374wfPz7rrrtuJb8lFKop/yYrenX30aNH58orr8z48ePz1FNP5eijj878+fNz+OGHJ0lGjBiRU089tWH9gQcemJ/85Ce57rrr8uKLL2bSpEk5/fTTc+CBBzbEOgAA8NE3e/bsnHLKKfnDH/6Qhx56KK+88kqSpG/fvpkyZUoeffTR3H777Tn66KMb9qmtrU3//v3z2GOPZY899sh///d/5+qrr86TTz6ZG264ITNmzGhYu/baa+fhhx/OMccck3333Tdnn312/vrXv2bGjBn5wx/+kCT5yle+kgcffDCPPfZY9tprr1xwwQWr9HvAR1NF35M+fPjwvPbaaznjjDMye/bsDBw4MHfeeWfDxeRmzpzZ6Mj5aaedlqqqqpx22ml5+eWXs/766+fAAw/MOeecU6kvAQAAqIA///nP2XXXXdOjR48kyRe+8IVceumlmTdvXr785S/nr3/9a9q0aZOXXnopr7/+etZbb71UVVXloIMOSpIMHDgws2bNyjrrrJMk2WKLLTJjxoz06dMnSXLwwQc3rOvbt28+9rGPJUkGDBiQF198Mbvvvnvuv//+jB07NnPnzs2//vWvbLHFFqv4u8BHUcUvHDdq1Kj3PL19ypQpjW63adMmY8aMyZgxY1bBZAAAwOqiqqoqSfLtb3872267ba6//vpUVVVl3XXXzTvvvJNk8Uc+t2nTpuHv//622FatWuXdd99tuL3kvvdat3Dhwhx22GGZNm1aNt100/z2t7/NJZdcstK/Tj76Knq6OwAAwIrYcccdc//99+fVV19NfX19rrnmmiTJ22+/nZ49e6aqqio33nhj3njjjZXy/O+8807q6urSrVu31NbW5mc/+9lKeR5aHpEOAACsdjbYYIOce+65+eQnP5ntt98+66yzTrp06ZJTTz015557bgYOHJh77703G2200Up5/i5dumT06NEZMGBAdt555/Tr12+lPA8tT8U+gq1Sqqur06VLl+W69D2rh5qamkyYMCH77bdfi7mSKQDAf2qJr4nmzZuXtdZaK/X19fnqV7+aj33sYznllFMqPRYspSkd6kg6AACwWho7dmwGDRqUrbbaKvPmzXvPa13B6qTiF44DAABYEeecc45PeuIjx5F0AABY3dXWpuree9PzD39I1b33JrW1lZ4IWEEiHQAAVmc335z06ZM2Q4Zk+4suSpshQ5I+fRZvB1Y7Ih0AAFZXN9+cHHJI8o9/NN7+8suLtwt1WO2IdAAAWB3V1ibHHZcs68Oalmw7/ninvsNqRqQDAMDq6I9/XPoI+r+rr09eemnxOmC1IdIBAGB1NGtW864DiiDSAQBgdbTBBs27DiiCSAcAgNXRbrslvXolVVXLvr+qKunde/E6YLUh0gEAYHXUunVyySWL//6fob7k9sUXL14HrDZEOgAArK4OPji58cakZ8/G23v1Wrz94IMrMxewwtpUegAAAOBDOPjg5DOfybv33JPpd9yRgfvumzZ77ukIOqymRDoAAKzuWrdO/e675+X587PN7rsLdFiNOd0dAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgEBWP9Msuuyx9+vRJhw4dstNOO2XatGnvu/6tt97KMccckw022CDt27dPv379MmHChFU0LQAAAKw8bSr55Ndff31Gjx6dyy+/PDvttFMuvvjiDB06NM8880y6deu21PpFixZlyJAh6datW2688cb07Nkzf//739O1a9dVPzwAAAA0s4pG+kUXXZQjjzwyhx9+eJLk8ssvz+23356rrroq3/zmN5daf9VVV+WNN97I/fffn7Zt2yZJ+vTpsypHBgAAgJWmYqe7L1q0KA8//HAGDx78f8O0apXBgwdn6tSpy9zntttuyy677JJjjjkm3bt3z9Zbb51zzz03tbW1q2psAAAAWGkqdiT99ddfT21tbbp3795oe/fu3fP0008vc58XXnghv//97/P5z38+EyZMyPPPP5+vfe1rqampyZgxY5a5z8KFC7Nw4cKG29XV1UmSmpqa1NTUNNNXQyUt+Tn6eQIALZnXRFCupvy7rOjp7k1VV1eXbt265Yorrkjr1q2z3Xbb5eWXX873v//994z0sWPH5qyzzlpq+8SJE9OxY8eVPTKr0KRJkyo9AgBAxXlNBOVZsGDBcq+tWKSvt956ad26debMmdNo+5w5c9KjR49l7rPBBhukbdu2ad26dcO2LbfcMrNnz86iRYvSrl27pfY59dRTM3r06Ibb1dXV6d27d/bZZ5907ty5mb4aKqmmpiaTJk3KkCFDGq5VAADQ0nhNBOVackb38qhYpLdr1y7bbbddJk+enGHDhiVZfKR88uTJGTVq1DL32XXXXXPttdemrq4urVotfjv9s88+mw022GCZgZ4k7du3T/v27Zfa3rZtW//x+ojxMwUA8JoIStSUf5MV/Zz00aNH58orr8z48ePz1FNP5eijj878+fMbrvY+YsSInHrqqQ3rjz766Lzxxhs57rjj8uyzz+b222/Pueeem2OOOaZSXwIAAAA0m4q+J3348OF57bXXcsYZZ2T27NkZOHBg7rzzzoaLyc2cObPhiHmS9O7dO3fddVe+8Y1vZMCAAenZs2eOO+64nHLKKZX6EgAAAKDZVPzCcaNGjXrP09unTJmy1LZddtklDzzwwEqeCgAAAFa9ip7uDgAAAPwfkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFWKFI/9vf/pbTTjsthx56aF599dUkyR133JG//OUvzTocAAAAtCRNjvR77703/fv3z5///OfcfPPNmTdvXpLksccey5gxY5p9QAAAAGgpmhzp3/zmN/Pd7343kyZNSrt27Rq277XXXnnggQeadTgAAABoSZoc6U888UQOOuigpbZ369Ytr7/+erMMBQAAAC1RkyO9a9eumTVr1lLbH3300fTs2bNZhgIAAICWqMmR/r//+7855ZRTMnv27FRVVaWuri733XdfTjzxxIwYMWJlzAgAAAAtQpMj/dxzz80WW2yR3r17Z968edlqq63yX//1X/nEJz6R0047bWXMCAAAAC1Cm6bu0K5du1x55ZU5/fTT8+STT2bevHkZNGhQ+vbtuzLmAwAAgBajyZG+xEYbbZSNNtqoOWcBAACAFq3Jkf6lL33pfe+/6qqrVngYAAAAaMmaHOlvvvlmo9s1NTV58skn89Zbb2WvvfZqtsEAAACgpWlypP/mN79ZaltdXV2OPvrobLrpps0yFAAAALRETb66+zIfpFWrjB49Oj/4wQ+a4+EAAACgRWqWSE+Sv/3tb3n33Xeb6+EAAACgxWny6e6jR49udLu+vj6zZs3K7bffnpEjRzbbYAAAANDSNDnSH3300Ua3W7VqlfXXXz8XXnjhB175HQAAAHhvTY70e+65Z2XMAQAAAC1es70nHQAAAPhwlutI+qBBg1JVVbVcD/jII498qIEAAACgpVquSB82bNhKHgMAAABYrkgfM2bMyp4DAAAAWjzvSQcAAIBCNPnq7rW1tfnBD36QX//615k5c2YWLVrU6P433nij2YYDAACAlqTJR9LPOuusXHTRRRk+fHjefvvtjB49OgcffHBatWqVM888cyWMCAAAAC1DkyP9mmuuyZVXXpkTTjghbdq0yaGHHpqf/vSnOeOMM/LAAw+sjBkBAACgRWhypM+ePTv9+/dPkqy11lp5++23kyQHHHBAbr/99uadDgAAAFqQJkd6r169MmvWrCTJpptumokTJyZJHnzwwbRv3755pwMAAIAWpMmRftBBB2Xy5MlJkq9//es5/fTT07dv34wYMSJf+tKXmn1AAAAAaCmafHX3733vew1/Hz58eDbeeOPcf//96du3bw488MBmHQ4AAABakiZH+jvvvJMOHTo03N55552z8847N+tQAAAA0BI1+XT3bt26ZeTIkZk0aVLq6upWxkwAAADQIjU50sePH58FCxbkM5/5THr27Jnjjz8+Dz300MqYDQAAAFqUFbpw3A033JA5c+bk3HPPzV//+tfsvPPO6devX77zne+sjBkBAACgRWhypC/RqVOnHH744Zk4cWIef/zxrLnmmjnrrLOaczYAAABoUVY40t955538+te/zrBhw7LtttvmjTfeyEknndScswEAAECL0uSru99111259tprc8stt6RNmzY55JBDMnHixPzXf/3XypgPAAAAWowmR/pBBx2UAw44ID//+c+z3377pW3btitjLgAAAGhxmhzpc+bMSadOnVbGLAAAANCiNfk96QIdAAAAVo4VvnAcAAAA0LxEOgAAABRCpAMAAEAhmhzpr7322nve98QTT3yoYQAAAKAla3Kk9+/fP7fffvtS2y+44ILsuOOOzTIUAAAAtERNjvTRo0fnv//7v3P00UfnX//6V15++eXsvffeOf/883PttdeujBkBAACgRWhypJ988smZOnVq/vjHP2bAgAEZMGBA2rdvn8cffzwHHXTQypgRAAAAWoQVunDcZpttlq233jozZsxIdXV1hg8fnh49ejT3bAAAANCiNDnS77vvvgwYMCDPPfdcHn/88fzkJz/J17/+9QwfPjxvvvnmypgRAAAAWoQmR/pee+2V4cOH54EHHsiWW26ZL3/5y3n00Uczc+bM9O/ff2XMCAAAAC1Cm6buMHHixOy+++6Ntm266aa57777cs455zTbYAAAANDSNPlI+pJAf/7553PXXXflX//6V5Kkqqoqp59+evNOBwAAAC1IkyP9n//8Z/bee+/069cv++23X2bNmpUkOeKII3LiiSc2+4AAAADQUjQ50r/xjW+kbdu2mTlzZjp27Niwffjw4bnjjjuadTgAAABoSVboPel33XVXevXq1Wh737598/e//73ZBgMAAICWpslH0ufPn9/oCPoSb7zxRtq3b98sQwEAAEBL1ORI32233fLzn/+84XZVVVXq6upy/vnnZ88992zW4QAAAKAlafLp7ueff3723nvvPPTQQ1m0aFFOPvnk/OUvf8kbb7yR++67b2XMCAAAAC1Ck4+kb7311nn22WfzyU9+Mp/5zGcyf/78HHzwwXn00Uez6aabrowZAQAAoEVo8pH0JOnSpUu+/e1vN/csAAAA0KItV6Q//vjjy/2AAwYMWOFhAAAAoCVbrkgfOHBgqqqqUl9fn6qqqobt9fX1SdJoW21tbTOPCAAAAC3Dcr0n/cUXX8wLL7yQF198MTfddFM22WST/PjHP8706dMzffr0/PjHP86mm26am266aWXPCwAAAB9Zy3UkfeONN274+2c/+9n88Ic/zH777dewbcCAAendu3dOP/30DBs2rNmHBAAAgJagyVd3f+KJJ7LJJpsstX2TTTbJX//612YZCgAAAFqiJkf6lltumbFjx2bRokUN2xYtWpSxY8dmyy23bNbhAAAAoCVp8kewXX755TnwwAPTq1evhiu5P/7446mqqspvf/vbZh8QAAAAWoomR/qOO+6YF154Iddcc02efvrpJMnw4cPzuc99LmuuuWazDwgAAAAtRZMjPUnWXHPNHHXUUc09CwAAALRoKxTpzz33XO655568+uqrqaura3TfGWec0SyDAQAAQEvT5Ei/8sorc/TRR2e99dZLjx49UlVV1XBfVVWVSAcAAIAV1ORI/+53v5tzzjknp5xyysqYBwAAAFqsJn8E25tvvpnPfvazK2MWAAAAaNGaHOmf/exnM3HixJUxCwAAALRoTT7dfbPNNsvpp5+eBx54IP3790/btm0b3X/sscc223AAAADQkjQ50q+44oqstdZauffee3Pvvfc2uq+qqkqkAwAAwApqcqS/+OKLK2MOAAAAaPGa/J50AAAAYOVYriPpo0ePztlnn50111wzo0ePft+1F110UbMMBgAAAC3NckX6o48+mpqamoa/v5eqqqrmmQoAAABaoOWK9HvuuWeZfwcAAACaj/ekAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCFEOgAAABRCpAMAAEAhRDoAAAAUQqQDAABAIUQ6AAAAFEKkAwAAQCGKiPTLLrssffr0SYcOHbLTTjtl2rRpy7Xfddddl6qqqgwbNmzlDggAAACrQMUj/frrr8/o0aMzZsyYPPLII9lmm20ydOjQvPrqq++734wZM3LiiSdmt912W0WTAgAAwMpV8Ui/6KKLcuSRR+bwww/PVlttlcsvvzwdO3bMVVdd9Z771NbW5vOf/3zOOuusfOxjH1uF0wIAAMDKU9FIX7RoUR5++OEMHjy4YVurVq0yePDgTJ069T33+853vpNu3brliCOOWBVjAgAAwCrRppJP/vrrr6e2tjbdu3dvtL179+55+umnl7nPn/70p/zsZz/L9OnTl+s5Fi5cmIULFzbcrq6uTpLU1NSkpqZmxQanKEt+jn6eAEBL5jURlKsp/y4rGulNNXfu3Hzxi1/MlVdemfXWW2+59hk7dmzOOuuspbZPnDgxHTt2bO4RqaBJkyZVegQAgIrzmgjKs2DBguVeW9FIX2+99dK6devMmTOn0fY5c+akR48eS63/29/+lhkzZuTAAw9s2FZXV5ckadOmTZ555plsuummjfY59dRTM3r06Ibb1dXV6d27d/bZZ5907ty5Ob8cKqSmpiaTJk3KkCFD0rZt20qPAwBQEV4TQbmWnNG9PCoa6e3atct2222XyZMnN3yMWl1dXSZPnpxRo0YttX6LLbbIE0880Wjbaaedlrlz5+aSSy5J7969l9qnffv2ad++/VLb27Zt6z9eHzF+pgAAXhNBiZryb7Lip7uPHj06I0eOzPbbb58dd9wxF198cebPn5/DDz88STJixIj07NkzY8eOTYcOHbL11ls32r9r165JstR2AAAAWN1UPNKHDx+e1157LWeccUZmz56dgQMH5s4772y4mNzMmTPTqlXFPykOAAAAVrqKR3qSjBo1apmntyfJlClT3nffq6++uvkHAgAAgApwiBoAAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAK0abSA8Dq4rnnkrlz3/v+Tp2Svn1X3TwAAMBHj0iH5fDcc0m/fh+87tlnhToAALDinO4Oy+H9jqCvyDoAAIBlEekAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOmwHDp1at51AAAAy9Km0gPA6qBv3+TZZ9//c9A7dVq8DgAAYEWJdFhOpQR4mzZt8u677zZpn+985zs544wzVtJEAABAc3G6O7QA3/nOdyo9AgAAsBxEOqyGvvWtb6V///7Zdttt88QTTyRJ6urq8u1vfzs77rhjBgwYkG9961tJkm984xupra3NwIEDM3jw4CTJ17/+9eywww7p379/vvCFL2ThwoUV+1oAAID/I9JhNVNbW5tevXrliSeeyJgxY3L44YcnSa6++uokybRp0zJ9+vQ8+eSTueOOO/KDH/wgrVu3zvTp03P33XcnSc4888w8+OCDeeKJJ7LOOus07AsAAFSW96TDamjkyJFJks985jM5/PDDM3/+/EyYMCGPPfZYbr/99iTJ/Pnz89xzz2Xfffddav9bb701P/nJT/LOO+/k7bffTl1d3SqdHwAAWLYijqRfdtll6dOnTzp06JCddtop06ZNe8+1V155ZXbbbbesvfbaWXvttTN48OD3XQ8tRX19fb7//e9n+vTpmT59ep577rkce+yxS62bMWNGzjzzzEyYMCFPPPFETjrppLzzzjsVmBgAAPhPFY/066+/PqNHj86YMWPyyCOPZJtttsnQoUPz6quvLnP9lClTcuihh+aee+7J1KlT07t37+yzzz55+eWXV/HkUDm/+MUvkiS/+93v8rGPfSxrrrlm9t1334aj40nyyiuvZPbs2UmSjh07Zv78+UmS6urqrLHGGll77bWzYMGChscCAAAqr+Knu1900UU58sgjG95Xe/nll+f222/PVVddlW9+85tLrb/mmmsa3f7pT3+am266KZMnT86IESNWycxQSa1bt85LL72UAQMGpE2bNg3vJz/iiCPyj3/8I9tvv32qqqqy5ppr5uqrr06PHj0yatSobLfddunVq1fuvvvu7LHHHtliiy2y/vrrZ8cdd3QkHQAAClFVX19fX6knX7RoUTp27Jgbb7wxw4YNa9g+cuTIvPXWW7n11ls/8DHmzp2bbt265YYbbsgBBxzwgeurq6vTpUuXvP322+ncufOHGZ9C1NTUZMKECdlvv/3Stm3bSo8DAFARXhNBuZrSoRU9kv7666+ntrY23bt3b7S9e/fuefrpp5frMU455ZRsuOGGDR8t9Z8WLlzY6OOlqqurkyz+j1hNTc0KTk5Jlvwc/TwBgJbMayIoV1P+XVb8dPcP43vf+16uu+66TJkyJR06dFjmmrFjx+ass85aavvEiRPTsWPHlT0iK1ttbdb961/T880388gTT+SfW22VtG5d6akAACpm0qRJlR4B+A8LFixY7rUVjfT11lsvrVu3zpw5cxptnzNnTnr06PG++15wwQX53ve+l7vvvjsDBgx4z3WnnnpqRo8e3XC7urq64WJzTndfvVX95jdpPXp0qv7tooH1PXum9qKLUn/QQRWcDABg1aupqcmkSZMyZMgQp7tDYZac0b08Khrp7dq1y3bbbZfJkyc3vCe9rq4ukydPzqhRo95zv/PPPz/nnHNO7rrrrmy//fbv+xzt27dP+/btl9retm1b//Fand18c/K//5v8xyUVql55JW3+93+TG29MDj64QsMBAFSO17lQnqb8m6z4R7CNHj06V155ZcaPH5+nnnoqRx99dObPn99wtfcRI0bk1FNPbVh/3nnn5fTTT89VV12VPn36ZPbs2Zk9e3bmzZtXqS+BVa22NjnuuKUCPcn/bTv++MXrAAAAViMVf0/68OHD89prr+WMM87I7NmzM3DgwNx5550NF5ObOXNmWrX6v98l/OQnP8miRYtyyCGHNHqcMWPG5Mwzz1yVo1Mpf/xj8o9/vPf99fXJSy8tXrfHHqtsLAAAgA+r4pGeJKNGjXrP09unTJnS6PaMGTNW/kCUbdas5l0HAABQiIqf7g5NtsEGzbsOAACgECKd1c9uuyW9eiVVVcu+v6oq6d178ToAAIDViEhn9dO6dXLJJYv//p+hvuT2xRf7vHQAAGC1I9JZPR188OKPWevZs/H2Xr18/BoAALDaKuLCcbBCDj44+cxn8u4992T6HXdk4L77ps2eezqCDgAArLZEOqu31q1Tv/vueXn+/Gyz++4CHQAAWK053R0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAApRRKRfdtll6dOnTzp06JCddtop06ZNe9/1N9xwQ7bYYot06NAh/fv3z4QJE1bRpAAAALDyVDzSr7/++owePTpjxozJI488km222SZDhw7Nq6++usz1999/fw499NAcccQRefTRRzNs2LAMGzYsTz755CqeHAAAAJpXxSP9oosuypFHHpnDDz88W221VS6//PJ07NgxV1111TLXX3LJJfnUpz6Vk046KVtuuWXOPvvsbLvttvnRj360iicHAACA5lXRSF+0aFEefvjhDB48uGFbq1atMnjw4EydOnWZ+0ydOrXR+iQZOnToe64HAACA1UWbSj7566+/ntra2nTv3r3R9u7du+fpp59e5j6zZ89e5vrZs2cvc/3ChQuzcOHChtvV1dVJkpqamtTU1HyY8SnEkp+jnycA0JJ5TQTlasq/y4pG+qowduzYnHXWWUttnzhxYjp27FiBiVhZJk2aVOkRAAAqzmsiKM+CBQuWe21FI3299dZL69atM2fOnEbb58yZkx49eixznx49ejRp/amnnprRo0c33K6urk7v3r2zzz77pHPnzh/yK6AENTU1mTRpUoYMGZK2bdtWehwAgIrwmgjKteSM7uVR0Uhv165dtttuu0yePDnDhg1LktTV1WXy5MkZNWrUMvfZZZddMnny5Bx//PEN2yZNmpRddtllmevbt2+f9u3bL7W9bdu2/uP1EeNnCgDgNRGUqCn/Jit+uvvo0aMzcuTIbL/99tlxxx1z8cUXZ/78+Tn88MOTJCNGjEjPnj0zduzYJMlxxx2X3XffPRdeeGH233//XHfddXnooYdyxRVXVPLLAAAAgA+t4pE+fPjwvPbaaznjjDMye/bsDBw4MHfeeWfDxeFmzpyZVq3+7yL0n/jEJ3LttdfmtNNOy7e+9a307ds3t9xyS7beeutKfQkAAADQLKrq6+vrKz3EqlRdXZ0uXbrk7bff9p70j4iamppMmDAh++23n1O7AIAWy2siKFdTOrSin5MOAAAA/B+RDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFEOkAAABQCJEOAAAAhRDpAAAAUAiRDgAAAIUQ6QAAAFAIkQ4AAACFaFPpAVa1+vr6JEl1dXWFJ6G51NTUZMGCBamurk7btm0rPQ4AQEV4TQTlWtKfS3r0/bS4SJ87d26SpHfv3hWeBAAAgJZk7ty56dKly/uuqapfnpT/CKmrq8srr7ySTp06paqqqtLj0Ayqq6vTu3fvvPTSS+ncuXOlxwEAqAiviaBc9fX1mTt3bjbccMO0avX+7zpvcUfSW7VqlV69elV6DFaCzp07+z8kAKDF85oIyvRBR9CXcOE4AAAAKIRIBwAAgEKIdFZ77du3z5gxY9K+fftKjwIAUDFeE8FHQ4u7cBwAAACUypF0AAAAKIRIBwAAgEKIdAAAACiESAcAAIBCiHRWW3/4wx9y4IEHZsMNN0xVVVVuueWWSo8EALBKjR07NjvssEM6deqUbt26ZdiwYXnmmWcqPRbwIYh0Vlvz58/PNttsk8suu6zSowAAVMS9996bY445Jg888EAmTZqUmpqa7LPPPpk/f36lRwNWkI9g4yOhqqoqv/nNbzJs2LBKjwIAUDGvvfZaunXrlnvvvTf/9V//VelxgBXgSDoAAHxEvP3220mSddZZp8KTACtKpAMAwEdAXV1djj/++Oy6667ZeuutKz0OsILaVHoAAADgwzvmmGPy5JNP5k9/+lOlRwE+BJEOAACruVGjRuV3v/td/vCHP6RXr16VHgf4EEQ6AACspurr6/P1r389v/nNbzJlypRssskmlR4J+JBEOqutefPm5fnnn2+4/eKLL2b69OlZZ511stFGG1VwMgCAVeOYY47Jtddem1tvvTWdOnXK7NmzkyRdunTJGmusUeHpgBXhI9hYbU2ZMiV77rnnUttHjhyZq6++etUPBACwilVVVS1z+7hx43LYYYet2mGAZiHSAQAAoBA+gg0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBoML22GOPHH/88ZUeAwAogEgHgI+Aq6++Ol27dq30GADAhyTSAQAAoBAiHQAK8O6772bUqFHp0qVL1ltvvZx++umpr69vuH/hwoU58cQT07Nnz6y55prZaaedMmXKlCTJlClTcvjhh+ftt99OVVVVqqqqcuaZZyZJfvGLX2T77bdPp06d0qNHj3zuc5/Lq6+++r6zzJo1K/vvv3/WWGONbLLJJrn22mvTp0+fXHzxxQ1rLrroovTv3z9rrrlmevfuna997WuZN29ew/1Ljuz/7ne/y+abb56OHTvmkEMOyYIFCzJ+/Pj06dMna6+9do499tjU1tY27NenT59897vfzYgRI7LWWmtl4403zm233ZbXXnstn/nMZ7LWWmtlwIABeeihhxr2+ec//5lDDz00PXv2TMeOHdO/f//86le/+hA/DQCoHJEOAAUYP3582rRpk2nTpuWSSy7JRRddlJ/+9KcN948aNSpTp07Nddddl8cffzyf/exn86lPfSrPPfdcPvGJT+Tiiy9O586dM2vWrMyaNSsnnnhikqSmpiZnn312Hnvssdxyyy2ZMWNGDjvssPedZcSIEXnllVcyZcqU3HTTTbniiiuWCvtWrVrlhz/8Yf7yl79k/Pjx+f3vf5+TTz650ZoFCxbkhz/8Ya677rrceeedmTJlSg466KBMmDAhEyZMyC9+8Yv8v//3/3LjjTc22u8HP/hBdt111zz66KPZf//988UvfjEjRozIF77whTzyyCPZdNNNM2LEiIZfYrzzzjvZbrvtcvvtt+fJJ5/MUUcdlS9+8YuZNm3aiv44AKBy6gGAitp9993rt9xyy/q6urqGbaecckr9lltuWV9fX1//97//vb5169b1L7/8cqP99t577/pTTz21vr6+vn7cuHH1Xbp0+cDnevDBB+uT1M+dO3eZ9z/11FP1SeoffPDBhm3PPfdcfZL6H/zgB+/5uDfccEP9uuuu23B73Lhx9Unqn3/++YZtX/nKV+o7duzY6LmHDh1a/5WvfKXh9sYbb1z/hS98oeH2rFmz6pPUn3766Q3bpk6dWp+kftasWe85z/77719/wgknvOf9AFCqNpX8BQEAsNjOO++cqqqqhtu77LJLLrzwwtTW1uaJJ55IbW1t+vXr12ifhQsXZt11133fx3344Ydz5pln5rHHHsubb76Zurq6JMnMmTOz1VZbLbX+mWeeSZs2bbLttts2bNtss82y9tprN1p39913Z+zYsXn66adTXV2dd999N++8804WLFiQjh07Jkk6duyYTTfdtGGf7t27p0+fPllrrbUabfvPo/QDBgxodH+S9O/ff6ltr776anr06JHa2tqce+65+fWvf52XX345ixYtysKFCxvmAIDViUgHgMLNmzcvrVu3zsMPP5zWrVs3uu/fg/c/zZ8/P0OHDs3QoUNzzTXXZP3118/MmTMzdOjQLFq0aIXnmTFjRg444IAcffTROeecc7LOOuvkT3/6U4444ogsWrSoIY7btm3baL+qqqplblvyi4Ml/n3Nkl9cLGvbkv2+//3v55JLLsnFF1/c8D75448//kN9jQBQKSIdAArw5z//udHtBx54IH379k3r1q0zaNCg1NbW5tVXX81uu+22zP3btWvX6AJsSfL000/nn//8Z773ve+ld+/eSdLogmvLsvnmm+fdd9/No48+mu222y5J8vzzz+fNN99sWPPwww+nrq4uF154YVq1Wnx5m1//+tdN+4Kb0X333ZfPfOYz+cIXvpBkcbw/++yzyzxTAABK58JxAFCAmTNnZvTo0XnmmWfyq1/9KpdeemmOO+64JEm/fv3y+c9/PiNGjMjNN9+cF198MdOmTcvYsWNz++23J1l8VfR58+Zl8uTJef3117NgwYJstNFGadeuXS699NK88MILue2223L22We/7xxbbLFFBg8enKOOOirTpk3Lo48+mqOOOiprrLFGwxHszTbbLDU1NQ2P+4tf/CKXX375yv0GvY++fftm0qRJuf/++/PUU0/lK1/5SubMmVOxeQDgwxDpAFCAESNG5F//+ld23HHHHHPMMTnuuONy1FFHNdw/bty4jBgxIieccEI233zzDBs2LA8++GA22mijJMknPvGJfPWrX83w4cOz/vrr5/zzz8/666+fq6++OjfccEO22mqrfO9738sFF1zwgbP8/Oc/T/fu3fNf//VfOeigg3LkkUemU6dO6dChQ5Jkm222yUUXXZTzzjsvW2+9da655pqMHTt25XxjlsNpp52WbbfdNkOHDs0ee+yRHj16ZNiwYRWbBwA+jKr6+n/7EFYAgP/wj3/8I717987dd9+dvffeu9LjAMBHmkgHABr5/e9/n3nz5qV///6ZNWtWTj755Lz88st59tlnl7rwGwDQvFw4DgBopKamJt/61rfywgsvpFOnTvnEJz6Ra665RqADwCrgSDoAAAAUwoXjAAAAoBAiHQAAAAoh0gEAAKAQIh0AAAAKIdIBAACgECIdAAAACiHSAQAAoBAiHQAAAAoh0gEAAKAQ/x8c8Ad/Ikq5mwAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"specification = SensitivityAnalysisMethodModel(\n",
"\texperiment_name=\"openturns_sensitivity_analysis\",\n",
"\tparameter_spec=parameter_spec,\n",
"\tobserved_data=observed_data,\n",
"\tmethod=\"chaos_sobol\",\n",
"\torder=2,\n",
"\tn_samples=256,\n",
"\toutput_labels=[\"Number of Infected\"],\n",
"\tcalibration_func_kwargs=dict(t=observed_data.day),\n",
")\n",
"\n",
"calibrator = SensitivityAnalysisMethod(\n",
"\tcalibration_func=sensitivity_func, specification=specification, engine=\"openturns\"\n",
")\n",
"\n",
"calibrator.specify().execute().analyze()"
]
},
{
"cell_type": "markdown",
"id": "e70bee40-59c0-44c2-9fca-5d8adcb82fdf",
"metadata": {},
"source": [
"Our surrogate model produces comparable sensitivity indices to the Sobol method."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.19"
}
},
"nbformat": 4,
"nbformat_minor": 5
}