Lab-part-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SIR epidemic model - [sensitivity analysis]\n",
    "\n",
    "In this notebook we will understand SIR epidemic model - its equations and parameters. \n",
    "Then we will perform sensitivity analysis of those parameters (with theoretical introduction)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SIR model\n",
    "\n",
    "`SIR` model describes dynamics of an epidemic for fixed numer of individuals (`N`) and given contact/infection rate together with mean recovery rate.  \n",
    "Our population is comprised of three groups: \n",
    "- `S` - Suspected\n",
    "- `I` - Infected\n",
    "- `R` - Recovered \n",
    "\n",
    "There's also a variant that includes those who had contact with infected person but haven't developed symptoms yet (`SEIR` with `E` - Exposed), and another variant with $\\alpha$ representing death rate caused by epidemic (coupled in $\\gamma+\\alpha$) and separate group `D` for deaths (`SIRD`).\n",
    "\n",
    "In `SIR` model where epidemic occurs in relatively short period of time we neglect birth rate and death rate (counting only non-epidemic deaths) in considered group. So we end up with two parameters: \n",
    "- $\\beta$ - contact rate (infected contacts $\\beta N$ other individuals per time step)\n",
    "- $\\gamma$ - mean recovery rate ($1/\\gamma$ days that recovery takes on average)\n",
    "\n",
    "<img src=\"https://i.ibb.co/XSgYy7w/Screenshot-2020-10-31-The-SIR-epidemic-model.png\" width=\"400\">\n",
    "\n",
    "At the end we will be interested in calculating reproduction ratio: \n",
    "$$R_0=\\frac{\\beta}{\\gamma}$$\n",
    "that will determine wheter the pandemic is in retreat ($R_0 < 1$) or in advance ($R_0 > 1$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model is described with given differential equations:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dS}{dt} & = -\\frac{\\beta SI}{N}\\\\\n",
    "\\frac{dI}{dt} & = \\frac{\\beta SI}{N} - \\gamma I \\\\\n",
    "\\frac{dR}{dt} & = \\gamma I\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert tf.__version__[0] == '2', 'Tensorflow 2.x version required'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2)\n",
    "tf.random.set_seed(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1: Complete equations according to formulas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deriv(y, N, beta, gamma):\n",
    "    S, I, R = y\n",
    "    dSdt = -beta * S * I / N\n",
    "    dIdt = beta * S * I / N - gamma * I\n",
    "    dRdt = gamma * I\n",
    "    return dSdt, dIdt, dRdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Euler scheme\n",
    "\n",
    "Euler scheme is a numerical procedure for integrating systems of differential equations which we will be using to solve SIR for our parameters.  \n",
    "We're looking for a solution $y(t)$ given $y'(t)=f(t,y(t))$ dynamics together with initial conditions $y(t_0)=y_0$.  \n",
    "\n",
    "With defined step size $h$ that $t_{n+1}=t_n+h$ we have formula:\n",
    "\n",
    "$$y_{n+1}=y_n+hf(t_n,y_n)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2: Complete the formula for each component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def euler(deriv, y0, dt, steps, N, beta, gamma):\n",
    "    S0, I0, R0 = y0\n",
    "    S_result, I_result, R_result = S0, I0, R0\n",
    "    for i in range(steps-1):\n",
    "        S1, I1, R1 = deriv((S0, I0, R0), N, beta, gamma)\n",
    "        next_s, next_i, next_r = S0 + dt*S1, I0 + dt*I1, R0 + dt*R1\n",
    "        S_result = tf.concat([S_result, next_s], axis=0)\n",
    "        I_result = tf.concat([I_result, next_i], axis=0)\n",
    "        R_result = tf.concat([R_result, next_r], axis=0)\n",
    "        S0, I0, R0 = next_s, next_i, next_r\n",
    "    return S_result, I_result, R_result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Total population\n",
    "N = tf.constant([1000.])\n",
    "\n",
    "# Initial infected and recovered individuals\n",
    "I0, R0 = tf.constant([1.]), tf.constant([0.])\n",
    "\n",
    "# Everyone else, S0, is susceptible to infection initially\n",
    "S0 = N - I0 - R0\n",
    "\n",
    "# Initial conditions tuple\n",
    "y0 = S0, I0, R0\n",
    "\n",
    "# Contact rate (beta), and mean recovery rate (gamma), (in 1/days)\n",
    "beta, gamma = 0.4, 1./10\n",
    "\n",
    "# A grid of time points (in days)\n",
    "days = 160\n",
    "t = np.linspace(0, days, days)\n",
    "dt = 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 3: Solve with euler scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "S, I, R = euler(deriv, y0, dt, days, N, beta, gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzmklEQVR4nO3deXwV1d348c/Jvdn3PUASEkjYBBQIiCwKghuiuC/to9j61Eef1q3VVuvP2r3a2ira1koVtZXHjaogailLQIgCskNYQwgJkIQASci+3Ht+f8wEE0iA5C5zk3zfr9d93Zk5c+d8M5DvnZw5c47SWiOEEKJ38LM6ACGEEN4jSV8IIXoRSfpCCNGLSNIXQoheRJK+EEL0InarAzibuLg4nZaWZnUYQgjRrWzcuPGY1jq+vTKfTvppaWls2LDB6jCEEKJbUUod7KhMmneEEKIXkaQvhBC9iCR9IYToRSTpCyFELyJJXwghepFzJn2l1Dyl1FGl1I5W22KUUkuVUvvM92hzu1JKvaSUylNKbVNKjW71mdnm/vuUUrM98+MIIYQ4m/O50n8TuPq0bU8Ay7XWmcBycx3gGiDTfN0HvALGlwTwDHAxMA54puWLQgghhPecs5++1voLpVTaaZtnAVPM5beAlcBPzO3/0MZ4zWuVUlFKqT7mvku11icAlFJLMb5I3nH9RzjT0ZP1/HPtQcIC7YQF2QkLtBMeZCc8yJ+IIH8SIwKJDPZHKeWJ6oUQwmd19eGsRK11sblcAiSay/2Aolb7HTK3dbT9DEqp+zD+SiA1NbVLwR2uqOPP2XmcbaqAQLsfiRFBpMQEMygxnMGJ4QxOCmdonwiC/G1dqlcIIXydy0/kaq21UsptM7ForecCcwGysrK6dNxRqdHk/3YGtY0Oqhuaqapvprqhmer6ZsprGzla1UDpyXpKT9ZTcLyW974uorbRAUCAzY+LUqOYOjiBa0f0ITU2xF0/mhBCWK6rSb9UKdVHa11sNt8cNbcfBlJa7ZdsbjvMN81BLdtXdrHu86KUIjTQTmigncSIs+/rdGoOldexq+QkGw+Wk5N3jOf+vZvn/r2bcWkx3DMxjSuHJWK3SWcn0XlaaxodjdQ119HkaKLR0XjGq8l55vaWfR3agcPpwKmdOLT57nS0WT69rGX99DKtNRpNy4x5LcsafSrWriy3OV4nP+/y+cU915zumkXQXfFMSJ7Aw+MfdsuxWutq0l8EzAaeNd8Xttr+A6XUuxg3bSvNL4YlwG9b3by9Eniy62G7l5+fIjU2hNTYEK66IAmAohO1LN5WzPx1B/nf+ZsYlBjGE9cMYergBLkX0EtorTnZcJKS6hJO1J2gor6izau8vpyK+goqGyqpaayhpqmG2qZaaptqqWlstdxUg1M7LfkZbMqGzc+Gn/I79QJQqFP/j1uWFeZ6F5dPP15nll3lrt9Jd8QC7omnX3i7LeAuU+f6dlNKvYNxlR4HlGL0wvkYeB9IBQ4Ct2mtTyjjJ/0zxk3aWuA7WusN5nG+C/zUPOxvtNZvnCu4rKwsbfWAaw6n5t87SvjDkt0UHK/l2hF9+M2Nw4kKCbA0LuEah9PBoZOHOFBxgAPlBzhYeZCS6hJKa0opqS459apvru/wGEH2IKKCoogMjCQsIIzQgFBC/EMI8Q8h1P/M5SB7EIH2QAJsAfj7+RNgCzj18redtm6W+9v8sfvZ8VN+2JSRvFuSeOuE3l6ZXJz0XkqpjVrrrHbLfHlidF9I+i2aHE7mfpHPC0v3EhcWyGuzsxjeL9LqsMQ51DbVkns0l+1Ht7Pj6A5yy3LZf2I/hZWFNDmb2uwbFxJHUljSN69Q4z0xLJH4kHiigqJOvSKDIgmyB1n0UwlxdpL03WjboQoeeHsTFbWNzL07i4kZcVaHJEwNzQ1sKt7El0VfsvbwWraUbGH/if2n2liD7EEMix9GZkwm6VHppEenMyB6AOlR6aREphBgk7/eRM8gSd/NSirrmT1vPfnHqnnrO+OYIInfEs3OZjYc2cDS/Uv5T/5/WH94PY2ORgDSo9IZ3Wc0IxJGMCJxBCMSRjAgegA2P+mOK3o+SfoeUFnbxK2vfklxZT0fPjCBzMRwq0PqFeqa6vg873MW7FzA53mfU1FfgUIxpu8YpqZN5ZLkS7gk5RKSwpKsDlUIy0jS95BD5bXc8JcvCfL349MHJxMZ4m91SD2Sw+lgyf4lvLX1LT7d+yk1TTXEBsdy/eDruTrjaqalTyM2JNbqMIXwGWdL+j49XaKvS44O4e93j+GWv33FLz7J5U+3X2R1SD3KkaojzNs8j79v+juFlYXEhcTxXyP/i1uG3cKUtCnY/eS/rxCdJb81LhqVGs0PpmYwZ/k+rrwgiauHS7OCq/Ye38vv1vyOt7e9TbOzmWnp03j+iueZNWSW3GwVwkWS9N3gB5dnsGxXKU99tJ0JGbFEBEkzT1fkHs3lN6t/w3u57xFgC+B/s/6XBy9+kIyYDKtDE6LHkHEF3MDf5sezN43keE0jc1flWx1Ot3O89jj3L76fEa+MYNGeRTx2yWMUPFzAnGvmSMIXws3kSt9NRiRHct2FfXltTT53X9KfhAh5cOdcHE4Hf9/0d55a8RSV9ZU8dPFDPH3p03JTVggPkit9N3rsykE0OzQvLNtndSg+b8+xPYx/fTwPfPoAIxNHsuX+Lbx49YuS8IXwMEn6btQ/NpQ7x6XywYYiSk92PGZLb6a15rVNrzF67mjyy/N55+Z3WHH3CoYnDLc6NCF6BUn6bnbvpHQcWjN/7UGrQ/E5J+pOcOsHt/K9T77H+OTxbLt/G3cMv0MGBhPCiyTpu1laXCiXD05g/rpCGpodVofjM7aXbmfM3DEs3LOQ30//PUvvWkq/CM8MHSuE6JgkfQ/4zsR0jtc08snW4nPv3Ass3ruYCfMm0NDcwJrvrOHxiY+fGtddCOFd8pvnARMzYslMCOMfXxVYHYrl5m2ex6x3ZzE4djBff+9rLk6+2OqQhOjVJOl7gFKK28emsO1QJfll1VaHY5k/fvlH7l10L9MHTGflPSulOUcIHyBJ30NmjuyLUrBo6xGrQ7HE818+z2NLH+PWYbfyyZ2fEBYQZnVIQggk6XtMUmQQF6fHsGjrEbdNuNxdzFk7h8eXPs5tF9zG/938fzJejhA+RJK+B11/YT/yy2rIPXLS6lC8Zv62+Tyy5BFuGnoT82+aLyNhCuFjJOl70DXDk/C3qV7TxLMsfxn3LLyHKWlTJOEL4aMk6XtQdGgAkzLiWJJbYnUoHrf3+F5u/eBWhsQN4ePbP5ZJw4XwUZL0PWzqkAQOHq/lwLEaq0PxmIr6Cq575zrsfnY+ufMTIoMirQ5JCNEBSfoeNmVQAgDZu49aHIlnaK255+N7yC/P51+3/Yu0qDSrQxJCnIUkfQ9LjQ1hQHwoK/eWWR2KR8xZN4eFexby3PTnuLT/pVaHI4Q4B0n6XjBlUAJr849T19izxuLZcGQDjy99nFmDZ/Ho+EetDkcIcR4k6XvBlMHxNDY7+Sr/mNWhuE1dUx13f3Q3iaGJvDHrDRkpU4huQpK+F4xLjyHY38bKPT2nieepFU+x69gu3pj1BtHB0VaHI4Q4T5L0vSDI38bY9BjW5h+3OhS3yCnM4cW1L/L9sd/nioFXWB2OEKITJOl7ybi0aPaWVlNe02h1KC5pdDRy3+L7SIlM4dnpz1odjhCikyTpe8nYtBgANhwstzgS1/wh5w/sLNvJX2b8RQZRE6IbcinpK6UeVUrlKqV2KKXeUUoFKaXSlVLrlFJ5Sqn3lFIB5r6B5nqeWZ7mlp+gm7gwJYoAmx9fF5ywOpQuO1B+gF998StuGXYLMwfNtDocIUQXdDnpK6X6AQ8BWVrr4YANuAN4DnhBa50BlAP3mh+5Fyg3t79g7tdrBPnbGJkcyfoD3Tfp/2TZT/BTfvzpyj9ZHYoQootcbd6xA8FKKTsQAhQDlwMLzPK3gBvM5VnmOmb5NNXL+vmNTY9hx+FKahubrQ6l01YfXM0HOz/gJxN/QkpkitXhCCG6qMtJX2t9GHgeKMRI9pXARqBCa92S1Q4BLdMl9QOKzM82m/vHdrX+7mhcWgzNTs2WwgqrQ+kUp3by6JJHSYlI4fGJj1sdjhDCBa4070RjXL2nA32BUOBqVwNSSt2nlNqglNpQVtZz+rUDjO4fjVKwvpu16y/YuYCNxRv57bTfEuIfYnU4QggXuNK8Mx04oLUu01o3AR8CE4Eos7kHIBk4bC4fBlIAzPJI4IyO61rruVrrLK11Vnx8vAvh+Z7IYH8GJYSztajC6lDOW7Ozmaezn2Z4wnC+NeJbVocjhHCRK0m/EBivlAox2+anATuBbOAWc5/ZwEJzeZG5jlm+Qve2eQSBEcmRbD9c2W2mUPzH1n+w9/hefj311/gp6eErRHfnSpv+OowbspuA7eax5gI/AX6olMrDaLN/3fzI60Csuf2HwBMuxN1tjUyO5Fh1I0cq660O5ZwaHY38YtUvGNdvHNcPvt7qcIQQbuDSfHZa62eAZ07bnA+Ma2ffeuBWV+rrCUYmRwGw/VAF/aKCrQ3mHN7e9jaFlYW8cu0rMqCaED2E/L3uZUOSwrH7KbYdqrQ6lLNyOB08u+ZZRiWN4pqMa6wORwjhJjJztZcF+dsYnBTO9sO+nfQX7FzAvhP7+ODWD+QqX4geRK70LTAyOZJth3z3Zq7WmmdznmVI3BBuGnqT1eEIIdxIkr4FRvSLorKuicITtVaH0q5VB1expWQLj13ymPTYEaKHkd9oC4xMjgTw2Xb9F9e+SFxInPTLF6IHkqRvgUGJ4fjbFDuO+F7S339iP4v2LOL+MfcT7O/bvYuEEJ0nSd8CAXY/BsaHsaekyupQzvDy+pex+9l5YOwDVocihPAASfoWGdongt3FvpX0axpreHPLm9x6wa30De9rdThCCA+QpG+RwUnhlJysp6LWd6ZPfC/3PSobKrl/zP1WhyKE8BBJ+hYZkhQOwG4fauJ5deOrDI0byqTUSVaHIoTwEEn6FhnaJwLAZ9r1NxdvZv3h9dyfdb88jCVEDyZJ3yIJ4YFEhfizu+Sk1aEAMHfjXILsQdw18i6rQxFCeJAkfYsopRiSFO4TzTv1zfW8m/suNw+9mejgaKvDEUJ4kCR9Cw1JimBPSRVOp7XDMSzas4iK+gpmXzj73DsLIbo1SfoWGpIUTm2jg0PldZbG8dbWt+gX3o/L0y+3NA4hhOdJ0rfQEPNm7i4L2/VLqktYkreEu0behc3PZlkcQgjvkKRvocyEMAD2lVrXrj9/23wc2sHsi6RpR4jeQJK+hUID7fSNDCLvaLVlMbyz4x2y+mYxJG6IZTEIIbxHkr7FBiaEkVdmTdLfd3wfG4s3cufwOy2pXwjhfZL0LZaREMb+ozWW9OB5L/c9AG4d1uunLhai15Ckb7GMhDDqmhwcqfR+D553d7zLpNRJpESmeL1uIYQ1JOlbLCPeuJnr7Xb9HUd3kFuWyx0X3OHVeoUQ1pKkb7GMBGuS/vu57+On/Lhl2C1erVcIYS1J+haLDQskOsSf/V6+mfvhrg+ZnDqZxLBEr9YrhLCWJH0fkJEQ5tUr/b3H95JblsuNQ270Wp1CCN8gSd8HeDvpf7TrIwBuHCpJX4jeRpK+DxgYH0Z5bRPHqxu8Ut9Huz9iTJ8xpEameqU+IYTvkKTvA7x5M/fwycOsO7yOm4be5PG6hBC+R5K+DxhodtvMP1bj8boW7VkEIO35QvRSkvR9QN+oYAJsfhR4Iekv3reYgdEDZawdIXopSfo+wOan6B8b4vEr/ZrGGpbnL2fmoJkyD64QvZRLSV8pFaWUWqCU2q2U2qWUukQpFaOUWqqU2me+R5v7KqXUS0qpPKXUNqXUaPf8CD1DWlyox6/0VxxYQYOjgZmDZnq0HiGE73L1Sn8O8G+t9RDgQmAX8ASwXGudCSw31wGuATLN133AKy7W3aMMiAvl4PFaHB4ceG3x3sWEBYRxaf9LPVaHEMK3dTnpK6UigUuB1wG01o1a6wpgFvCWudtbwA3m8izgH9qwFohSSvXpav09TVpcKI0OJ0cqPDPwmtaaT/d9ylUDryLAFuCROoQQvs+VK/10oAx4Qym1WSn1mlIqFEjUWheb+5QALc/59wOKWn3+kLmtDaXUfUqpDUqpDWVlZS6E172kx4UCcMBDTTxbS7dyuOqwNO0I0cu5kvTtwGjgFa31KKCGb5pyANBaa6BT7RVa67la6yytdVZ8fLwL4XUvLUm/4Lhnkv7n+z4H4OqMqz1yfCFE9+BK0j8EHNJarzPXF2B8CZS2NNuY70fN8sNA64Hbk81tAkgIDyQkwEZ+mWeS/pL9S7gw8UKSwpI8cnwhRPfQ5aSvtS4BipRSg81N04CdwCKgZZbt2cBCc3kRcLfZi2c8UNmqGajXU0qRFhvqkSv9qoYqcopyuGrgVW4/thCie7G7+PkHgflKqQAgH/gOxhfJ+0qpe4GDwG3mvp8BM4A8oNbcV7SSHh/KjsOVbj9udkE2zc5mrsqQpC9Eb+dS0tdabwGy2ima1s6+Gvi+K/X1dAPiQvn3jhIam50E2N333NySvCWE+IcwMWWi244phOie5IlcH5IWG4rDqSkqr3Xrcf+T/x+mpk0l0B7o1uMKIbofSfo+JD3e7MHjxm6b+eX55J3Ik/Z8IQQgSd+npMe6v6/+8vzlAEwfMN1txxRCdF+S9H1IdGgAUSH+bk362QXZJIUlyaiaQghAkr7PSYsNdVvS11qTXZDN1LSpMqqmEAKQpO9zBrhxtM09x/dQUl3C1LSpbjmeEKL7k6TvY9LiQjlSWU9do8PlY2UfyAZgarokfSGEQZK+j3HnGDzZBdkkRyQzMHqgy8cSQvQMkvR9zKmk72ITj9aalQUrpT1fCNGGJH0fk2YmfVenTswty6Wstkza84UQbUjS9zFhgXbiwwNdvtKX9nwhRHsk6fug9DjXu22uKFhBelQ6aVFp7glKCNEjSNL3QekuDrHs1E5WFaySph0hxBkk6fug9PhQjlU3crK+qUuf31qylfL6cmnaEUKcQZK+Dzo1X24XZ9HKLjDb8+VKXwhxGkn6PmiAi331swuyyYzJpF/EGfPOCyF6OUn6Pig1NgSl6NJ8uQ6ngy8OfiFX+UKIdknS90GBdhvJ0cFd6sGzrXQbJxtOclnaZR6ITAjR3UnS91HpcWFdSvo5RTkAMjWiEKJdkvR91IC4UPLLqjGmFj5/OUU5JEckkxqZ6qHIhBDdmSR9H5UeF0pNo4OyqoZOfW5N4RompkyU8XaEEO2SpO+j0rswBk9hZSGHTh6Sph0hRIck6fuoAfGdny83p9Bsz0+VpC+EaJ8kfR/VNzKYALtfp5L+msI1hAWEMTJxpAcjE0J0Z5L0fZSfnyI9NrRTffVzinIYnzweu5/dg5EJIbozSfo+zBhts/q89j3ZcJLtR7dLe74Q4qwk6fuw9PhQCk/U0uxwnnPftYfW4tROSfpCiLOSpO/D0uNCaXJoDlfUnXPfNYVr8FN+jE8e74XIhBDdlTT++rABrbpt9o8N/aZg82Z4913YtQsyM+Hpp8kpyuHCxAsJDwy3KFohRHcgV/o+rN0hlnNyYPJkePFFyMuDF1+kecgg1h38Upp2hBDnJEnfh8WEBhARZP+m2+amTTBjBvTtCwcPws6dsGEDW4dEUeOsZ5JOtjZgIYTPcznpK6VsSqnNSqnF5nq6UmqdUipPKfWeUirA3B5orueZ5Wmu1t3TKaVIjzcHXnM64d57ISICli+HpCRjp1GjWPPU3QBMfOwlOHrUwoiFEL7OHVf6DwO7Wq0/B7ygtc4AyoF7ze33AuXm9hfM/cQ5DDQHXuPDD2HLFvjd7yAlpc0+ORXbSA1OIrngBDzwgDWBCiG6BZeSvlIqGbgWeM1cV8DlwAJzl7eAG8zlWeY6Zvk0JaOCnVN6XCgl5TU4n/4ZDB0Kd97ZplxrTU5RDhMHToWnnza+HJYtsyhaIYSvc/VK/0Xgx0BLR/JYoEJr3WyuHwJa5uzrBxQBmOWV5v5tKKXuU0ptUEptKCsrczG87i89PpTrdn2B3+5d8Mtfgs3WprygooAjVUeYlDoJfvQjGDgQHnoImro2qboQomfrctJXSs0EjmqtN7oxHrTWc7XWWVrrrPj4eHceultKjwvlW1uXUN1/ANx00xnlbSZNCQw0evXs2gV/+5uXIxVCdAeuXOlPBK5XShUA72I068wBopRSLf3/k4HD5vJhIAXALI8EjrtQf6+QXn2Mi4t2sGPaLPA7858rpzCHiMAIhicMNzbMnAlTphht//X13g1WCOHzupz0tdZPaq2TtdZpwB3ACq31t4Fs4BZzt9nAQnN5kbmOWb5Cd3ZaqF4oZMF7ACwbNb3d8pZB1mx+rZp9nnkGiovhtde8EaIQohvxRD/9nwA/VErlYbTZv25ufx2INbf/EHjCA3X3LFrDP//JroyL2OQXeUZxRX0FO47uYFLKpLYFU6bApZfCs89CQ+dm3hJC9GxuSfpa65Va65nmcr7WepzWOkNrfavWusHcXm+uZ5jl+e6ou0fbtAl272b39OvbnUHrq6Kv0Oj2J0155hk4fBjefNPzcQohug15IteXLVgAdjvV199IRW0Tx6rbXrXnFOVgUzYu7nfxmZ+dOhXGjDFu7DrPPUqnEKJ3kKTvyz79FCZPJi3D6PW6r7Tt2PprCtdwUdJFhAaEnvlZpeDRR2H3bliyxBvRCiG6AUn6vqqoCLZvhxkzyEwwRs7cd7TqVHGTo4n1h9cb/fM7cuutxjg9L7zg6WiFEN2EJH1f9dlnxvu115IYEUh4kL3Nlf7mks3UNdedfWTNgAD4wQ9g6VLYscPDAQshugNJ+r7q008hLQ2GDEEpRWZCGHtLv7nSzyk0H8pq7yZua//zPxAcbLTtCyF6PUn6vqi+3hhJc8YMo20eGJQYTt7Rb6701xStIT0qnb7hfc9+rJgYmD0b3n5bRuAUQkjS90lffAG1tXDttac2ZSSEcbymkePVDcYga4U5577Kb/HII0Z/fRmaQYheT5K+L1q+HPz94bLLTm0alGjczN1bWk1+eT6lNaXnP1PW4MHGXw1/+YsMzSBELydJ3xetXAnjxkHoN10xW5J+3tGqtoOsna9HHzWad957z52RCiG6GUn6vubkSdi40Xi4qpXEiEDCA+3sLa1mTeEaIgMjuSDhgvM/7rRpMHy4cUNXhjwSoteSpO9r1qwBh8MYP6cVpRSZiWHsO1rFmsI1TEiZgJ/qxD+fUkbb/pYtxj0DIUSvJEnf12RnG/3rL7nkjKJBieHkFhex69guLu1/aeeP/a1vQVycdN8UoheTpO9rVq6Eiy+GkJAziob2iaC0YStA15J+cDDcfz8sXAj5Mt6dEL2RJH1fUllpjKx5Wnt+iwv6RlBvyyXAFkRW36yu1fHAA2C3w8svuxCoEKK7kqTvS1avNkbEPK09v8WQPhE0+OWSGjqSAFtA1+ro2xduuw1ef924aSyE6FUk6fuS7Gxjntt22vMBNHU0+uUTaRvpWj2PPAJVVfDGG64dRwjR7UjS9yUrV8L48RAU1G7xl0VfAk4aawe5Vk9WFkycCC+9ZPQUEkL0GpL0fUV5OWze3GF7PsDqwtX4YaOiIp2T9U2u1ffII8bN3MWLXTuOEKJbkaTvK1avNh6a6qA9H+CLg1+QGTMCP4LZXVzV4X7n5YYbIDVVum8K0ctI0vcV2dlGs87F7Ux9CDQ0N7D+8Hqmphvj8ew8UulafXY7PPig0aS0ZYtrxxJCdBuS9H3FypXGDdwO2vO/PvI1DY4GrsqYQlxYALlH3NDz5t57jfF95sxx/VhCiG5Bkr4vOHECtm49a3v+FweNoRMmp05maJ8Idrgj6UdHwz33wPz5cOiQ68cTQvg8Sfq+4Isvztmev7pwNRfEX0BsSCyjUqLYU3KS2sZm1+t+7DHj2YDnn3f9WEIInydJ3xdkZxtDJIwb126xw+kgpzCHyamTARiVGo1Tw9YiF9v1wZiS8a67YO5cmVlLiF5Akr4vWLkSJkwwHsxqx9bSrVQ1Vp0ab+eilCgANheVu6f+J54wJld54QX3HE8I4bMk6Vvt2DHYtu3s/fMPrgZgcn/jSj86NIABcaFsOljhnhgGDzaGZvjLX4znBYQQPZYkfau1jG1/lvb8FQUrGBg9kOSI5FPbRqVGs6WoHO2uCVF++lNjaAYZiE2IHk2SvtWys41hlMeObbe42dnMyoKVTEuf1mb7qNQojlU3UnSizj1xjBwJs2YZD2tVufjglxDCZ0nSt9rKlcY4OAHtj5q54cgGTjacZPqA6W22j06NBtzYrg/w1FNG884rr7jvmEIInyJJ30plZbBjx1mbdpblLwNganrbNv9BiWGEBNjYdNCNSX/sWLjqKvjDH+RqX4geqstJXymVopTKVkrtVErlKqUeNrfHKKWWKqX2me/R5nallHpJKZWnlNqmlBrtrh+i21q+3HifPr3jXQ4sZ1TSKOJC4tpst9v8GJUaxboDJ9wb069+Zdxclp48QvRIrlzpNwM/0loPA8YD31dKDQOeAJZrrTOB5eY6wDVApvm6D5A2hKVLISoKxoxpt7i2qZYvi748oz2/xYSBcewuqeJYdYP7Yho7Fm66yXhY69gx9x1XCOETupz0tdbFWutN5nIVsAvoB8wC3jJ3ewu4wVyeBfxDG9YCUUqpPl2tv9vT2kj6l18ONlu7u6wpXEOjo/GM9vwWEzOMq/+v9h93b2y//jXU1MBvf+ve4wohLOeWNn2lVBowClgHJGqti82iEiDRXO4HFLX62CFz2+nHuk8ptUEptaGsrMwd4fmmvDwoKjpr086SvCUE2AKYlDqp3fIR/SIJD7KTk+fmK/KhQ+E734E//xn27XPvsYUQlnI56SulwoB/AY9orduMAqaNTuSd6kiutZ6rtc7SWmfFx8e7Gp7vWrrUeL/iig53+SzvM6akTSE0ILTdcpuf4pIBsazed8x9/fVb/OpXRo+iH//YvccVQljKpaSvlPLHSPjztdYfmptLW5ptzPeWAV0OAymtPp5sbuudli2D/v1h4MB2i/PL89l9bDczMmac9TCTMuM4XFFH4Yla98bXpw88+SR8/LHxLIEQokdwpfeOAl4Hdmmt/9SqaBEw21yeDSxstf1usxfPeKCyVTNQ79LcDCtWGE07SrW7y+f7PgdgRubZk/6EgUa7fk6em9v1AX74Q2N2rYcegiYXp2cUQvgEV670JwJ3AZcrpbaYrxnAs8AVSql9wHRzHeAzIB/IA/4O/K8LdXdvX34JlZVw9dUd7vJZ3mdkxGSQGZt51kMNjA+lT2QQq/Z6YITM4GBjgpUdO4xJ1IUQ3Z69qx/UWq8B2r9MhTP6GJrt+9/van09yuLF4O8PV17ZbnFdUx0rDqzgvtH3nfNQSimuGJbI+xuKqG1sJiSgy/+k7Zs1C2bOhGeegdtvh+Tkc39GCOGz5IlcKyxeDJddBhER7RZnF2RT31x/zqadFldfkER9k5NVezzQ20kp4yrf6YTvf9/oaiqE6LYk6Xvb/v2wa5dx9dyBD3d9SHhAOJelXXZehxyXHkN0iD//zi1xV5RtpacbvXkWLYJ33/VMHUIIr5Ck722ffmq8d5D0m53NfLz7Y64bfB1B9vYnST+d3ebHlcOSWLHrKA3NDndF2tYjj8D48fDgg1Ba6pk6hBAeJ0nf2xYvhiFDOuyquapgFcfrjnPL0Fs6ddirRyRR1dDMl57oxQPGU8Pz5kF1NXzve9LMI0Q3JUnfm06cMIZSvu66DndZsHMBof6hXJ3Rcc+e9kwYGEt4kJ2FWzz46MPQofDcc/DJJ/C3v3muHiGEx0jS96aPPjL6u992W7vFDqeDj3Z/xIzMGQT7B3fq0IF2GzeO6sdnO0oor2l0R7Tte+gho6vpD39odOUUQnQrkvS96b33jGadDkbVXF24mtKaUm4eenOXDv+ti1NpbHby4WYPXu0rBW+8AZGRcPPNcPLkuT8jhPAZkvS9pazMeAr39ts7fAr3ra1vER4QznWDO27+OZshSRGMSo3infWF7h+Lp7WkJOMLbP9++O53pX1fiG5Ekr63/Otf4HAYSb8d1Y3VfJD7AbdfcDsh/iFdrubOcankHa3m6wI3zqjVnssug2efNX4uGYJZiG5Dkr63vPOOcSN0xIh2ixfsXEBNUw33XHSPS9XMHNmHyGB/5n6x36XjnJcf/Qi+/W34f/8PPvjA8/UJIVwmSd8b9u6FL76A//qvDpt23tzyJpkxmUxImeBSVSEBdu6dlM6yXUfZcbjSpWOdk1Lw2mswYQLcfTfk5Hi2PiGEyyTpe8PcuWC3G+3f7cg7kceqg6u456J7UB18KXTGPRPTiAiyM2e5FyZACQoyhl9OSTEeONu+3fN1CiG6TJK+p9XXw5tvGgOXJSW1u8tL617C38/f5aadFhFB/tw7aQBLd5Z6/mofID7emBQmNNQYRG7PHs/XKYToEkn6nvbhh3D8OPzP/7RbXFFfwbzN87hj+B30De/rtmrvmZhGTGgATy/cgdPphd41/fsbid/phKlTYfduz9cphOg0SfqepLUxz+yAATDtjNGmAfj7xr9T01TDo+MfdWvVkcH+/HTGUDYXVvDu10Xn/oA7DB1qzLLldMKUKbBli3fqFUKcN0n6nrRqFXz1FTz6KPideaqbHE28vP5lpqRNYVSfUW6v/ubR/bg4PYZnP9/F0ap6tx+/XcOGGUNNBATApZfKVItC+BhJ+p70m99AYiLce2+7xfM2z6PoZBGPXfKYR6pXSvGbG0fQ6HDy8DtbaHY4PVLPGYYMMWYHS0kx2vj/+ld5gEsIHyFJ31PWrjUmP3/8cWPawdPUNtXyi1W/YGLKxPOeLKUrMhLC+PUNI/gq/zh/XLrXY/WcITnZ6MJ51VXG5Cvf/S7U1XmvfiFEuyTpe4LW8PTTEBvb4Q3cl9e9THF1Mb+b9ju3dNM8m1vGJHPnuFReWbmfBRsPebSuNqKijIlXnnnG6ME0eTIcPOi9+oUQZ5Ck7wkffWRc5f/85xAWdkZxaXUpz+Y8y4zMGUzuP9krIf38+mFMzozjxwu2snjbEa/UCRj3Mn7+cyP579sHI0fC669Lc48QFpGk7261tcaN2xEj4P77293loX8/RG1TLc9f8bzXwgq023j1rjGM6R/Nw+9u4Z31hV6rGzDmENi8GUaPhv/+b6Otv6DAuzEIISTpu90vfwmFhUZXTbv9jOJFexbxfu77/OzSnzE0fqhXQwsJsPPGd8YxKSOOJz/czi8/2Uljs5du7oLRdXX5cnjlFeOex4gR8Mc/QkOD92IQopeTpO9Oy5bB739vXMleeukZxSXVJdy/+H5GJo7kxxN/bEGAEBZo5/XZWdwzIY15OQe48a857Cmp8l4Afn7GX0C5ucY5euwxGDwY5s83+vcLITxKkr67lJYaA6oNGQJz5pxR3Oho5Jb3b6GyoZK3b3wbf5u/BUEa7DY/fn79Bbx61xhKKuu59qXV/HxRLhW1Hpxx63SpqcYk8UuXQkyMce6ysoz7IQ4PTe4uhJCk7xZVVUabdWWlMblISNvx8LXW/OCzH5BTlMO86+cxIrH94ZW97aoLkvjPo5dy29gU/vFVAZOey+a3n+2iuNKLXSunT4cNG4wr/YoKuOkm44vzr3817o8IIdxKeXSGJRdlZWXpDRs2WB3G2TU0GKNLZmcbE4rMmtWmWGvNo0seZc66OTw56Ul+O803JxzZU1LFn7Pz+NTs2XPZoHhuGNWPKYMSiAzx0l8lzc3Glf7zz8P69cZfAN/+NsyebdwA9nDXViF6CqXURq11VrtlkvRdUF5uzBObnW30Q589u01xs7OZhz9/mL9u+CsPX/wwL1z1gsf75Luq6EQt731dxIKNhyg5WY/NTzGmfzTThiQwfkAsQ/tEEGD38B+IWhsPdr38MixcaHyxDh9ujNl/442QkeHZ+oXo5iTpe8LOnUbC37/fmEjk7rvbFJfVlHH7gtvJLsjmR5f8iD9c8QefT/itOZyarYcqWLHrKMt3H2VXsTEBeoDNj2F9I7goJYpBieEMjA9lYEIYsaEBnvn5ysuNJrO33jJ6/IAxsNusWXDttTBunDHOjxDiFEn67tTUBC+8YDxxGxEBCxYY88WatNa8u+NdHlnyCJX1lbw681VmXzT7LAfsHoor69hcWMGWIuO1/VAldU3f3HCNCLLTNyqYpMgg+kQGkRQRTJ/IIOIjAokM9icq2J9I82W3dfEvhQMH4JNPjKv/VauMG74hITBpkjGc86RJcNFF7T4QJ0RvIknfHZqa4P33jSEF9u83bji+8gokJABGsl9xYAW/WPULVheuZmzfsbx2/WuMTBxpceCe4XRqjlTWsb+shv1Hqyk4XkNxZT0llfUUV9ZzrLrjvvdhgfZTXwBhgXaCAmwE+/sR7G8jOMBGsL+d4ABjPch8Bdj8sNsUdpsf/n6K4OqTxG34kpivc4ham0PIPmP8fq0UTZmDaBx5EU0XjaL5otE4hg1FxcahFPgpZb6MAen8Wm1rKbf5fVMuRHckSb+rtDam//vgA2PogOJiuPBCY/TMGTNAKYqrivlg5we8vvl1tpVuo194P56a/BT3jbkPm5/Nutgt1tjspPRkPWXVDVTWNVFZ20RlXRMVtU1U1DWe2lbT2Exdk5P6Rgd1TcarvtFBbZMDRycmf4mtqWBkyT5GFu9jeGkeI0rySKo+caq8PCicAzF9ORDTj/zofhRE9+VIRDxHIuIoC43G2cG/1elfCq2/BxTGSsu21l8RLV8Y6pudOX3x1D5nKzu13jqqs9Xbfmy+xhfD8rUv+WlDE/jlrOFd+uzZkv6Zj4x6mFLqamAOYANe01o/6+0YOlRfD9u2GW3Ha9cawwMfPGj85lxzDY5XX6HgkmFsP5bLl8t+QnZBNhuOGF9KWX2zeHXmq9x94d0E2YMs/kGsF2D3IyUmhJSYkHPv3IHGZqfxJdDkoMnhpNmhaXY6aXJomhzGe7PDSbNT0+hw0uyYTrPDSZ1Tk9PsxL+slKjdOwg7uJ/wgwfoU5jP4KJdhO5Y0aYep81GbUw81Ql9qI5JoC4ymtqIKGojY6iNiKY2LIqaiGhqwiNpCI2gMSgEp834kmi5aGp97dSy2LJN803h6ddYrS+6Ovrc2Y5NO8dur15f4ovXmT4YEpmJ4R45rlev9JVSNmAvcAVwCPgauFNrvbO9/bt8pa+10eOjpsbo6936vaoKjh2DY8fQx8qoLy6iqjCPqiMFVFWUUhUAVYFQlhxNyaB+FA9M4EhSKHtrCtlzfA/1zcZkJAG2AMb2HcuMzBnMGjyLCxIu6OppEd5WU2M00R06BEVFbd+PHDGmtzx27OxPCIeEQHh4+6+wMGPC+JZXYODZl/39jSE7Wl42W8frHZXZbJzxp4jotXzpSn8ckKe1zgdQSr0LzALaTfpdtX3ZfG5feBdNftBkgyY/aG613GQz1h2BQJr5OkM5UE54czhJFUlkxmYyfcB0hsYNZVj8MEb1GSVX9N1VaKgx2ufIs9xvcTqNh8XMCwSOHYOyMuMBvKqq9l/FxbB3r7Hc0GD85Vhfb82lrVLGkBctXwTns96Zz5yt3s5s9+XPuLv+zrrmGmNsKjfzdtLvB7SesPUQcHHrHZRS9wH3AaSmpnapkpDUDC5IGI6/PQB//0D87YHY/QPxDwjCPyAY/9AI7KHh+AeHEuQfTHhAOOGB4W3e40LiSApLIjQgtIs/qujW/PyMh8NiYmDQoK4fR2vjobP6+rZfBK2Xm5qMfRwO4/305fMt0/rMl9PZufVz7XO2v346+nI725eer37G3fV3RUqKe49n8nqb/rlorecCc8Fo3unKMQYOHs8HP9vu1riE6BKljOYbf3+j6UcIi3l77J3DQOuvr2RzmxBCCC/wdtL/GshUSqUrpQKAO4BFXo5BCCF6La8272itm5VSPwCWYHTZnKe1zvVmDEII0Zt5vU1fa/0Z8Jm36xVCCCHj6QshRK8iSV8IIXoRSfpCCNGLSNIXQohexKdH2VRKlQEHXThEHHDMTeG4k8TVORJX50hcndMT4+qvtY5vr8Cnk76rlFIbOhp0yEoSV+dIXJ0jcXVOb4tLmneEEKIXkaQvhBC9SE9P+nOtDqADElfnSFydI3F1Tq+Kq0e36QshhGirp1/pCyGEaEWSvhBC9CI9Mukrpa5WSu1RSuUppZ6wMI4UpVS2UmqnUipXKfWwuT1GKbVUKbXPfI+2KD6bUmqzUmqxuZ6ulFpnnrf3zOGvvR1TlFJqgVJqt1Jql1LqEl84X0qpR81/wx1KqXeUUkFWnC+l1Dyl1FGl1I5W29o9P8rwkhnfNqXUaC/H9Qfz33GbUuojpVRUq7Inzbj2KKWu8lRcHcXWquxHSimtlIoz1y09Z+b2B83zlquU+n2r7e45Z1rrHvXCGLJ5PzAACAC2AsMsiqUPMNpcDseYFH4Y8HvgCXP7E8BzFsX3Q+D/gMXm+vvAHeby34AHLIjpLeC/zeUAIMrq84UxzecBILjVebrHivMFXAqMBna02tbu+QFmAJ8DChgPrPNyXFcCdnP5uVZxDTN/LwOBdPP31ebN2MztKRjDvB8E4nzknE0FlgGB5nqCu8+ZV35pvPkCLgGWtFp/EnjS6rjMWBYCVwB7gD7mtj7AHgtiSQaWA5cDi83/5Mda/ZK2OY9eiinSTK7qtO2Wni++mds5BmM48sXAVVadLyDttETR7vkBXgXubG8/b8R1WtmNwHxzuc3vpJl4L/HmOTO3LQAuBApaJX1LzxnGhcT0dvZz2znric077U2+3s+iWE5RSqUBo4B1QKLWutgsKgESLQjpReDHQMtM17FAhda62Vy34rylA2XAG2az02tKqVAsPl9a68PA80AhUAxUAhux/ny16Oj8+NLvwncxrqDBB+JSSs0CDmutt55WZHVsg4DJZrPhKqXUWHfH1ROTvs9RSoUB/wIe0VqfbF2mja9tr/abVUrNBI5qrTd6s97zYMf4c/cVrfUooAajueIUi85XNDAL40upLxAKXO3NGM6XFefnXJRSTwHNwHyrYwFQSoUAPwV+ZnUs7bBj/EU5HngceF8ppdxZQU9M+j41+bpSyh8j4c/XWn9obi5VSvUxy/sAR70c1kTgeqVUAfAuRhPPHCBKKdUym5oV5+0QcEhrvc5cX4DxJWD1+ZoOHNBal2mtm4APMc6h1eerRUfnx/LfBaXUPcBM4NvmF5IvxDUQ4wt8q/k7kAxsUkol+UBsh4APtWE9xl/ice6MqycmfZ+ZfN38hn4d2KW1/lOrokXAbHN5NkZbv9dorZ/UWidrrdMwzs8KrfW3gWzgFgvjKgGKlFKDzU3TgJ1YfL4wmnXGK6VCzH/TlrgsPV+tdHR+FgF3mz1SxgOVrZqBPE4pdTVGE+L1Wuva0+K9QykVqJRKBzKB9d6KS2u9XWudoLVOM38HDmF0uCjB4nMGfIxxMxel1CCMzgzHcOc58+TNE6teGHfg92Lc4X7KwjgmYfypvQ3YYr5mYLSfLwf2Ydypj7Ewxil803tngPkfKQ/4ALMHgZfjuQjYYJ6zj4FoXzhfwC+A3cAO4J8YvSi8fr6AdzDuKzRhJKt7Ozo/GDfn/2L+HmwHsrwcVx5GO3TL//2/tdr/KTOuPcA13j5np5UX8M2NXKvPWQDwtvn/bBNwubvPmQzDIIQQvUhPbN4RQgjRAUn6QgjRi0jSF0KIXkSSvhBC9CKS9IUQoheRpC+EEL2IJH0hhOhF/j/N4a02A8PRYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(t, S)\n",
    "ax.plot(t, I, color='r')\n",
    "ax.plot(t, R, color='g')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sensitivity analysis\n",
    "\n",
    "In sensitivity analysis we are interested in understanding how minor changes in input space contribute to output results for our black-box model defined as a continous scalar function $\\phi: \\mathbb{R}^p \\to \\mathbb{R}$.  \n",
    "\n",
    "Sensitivity analysis methods help us identify influential paremeters. There are global SA methods that compute sensitivity for whole input space, and local SA methods that perform it for given subarea. Both are widely used in physics, finance and engineering. \n",
    "\n",
    "### Global SA\n",
    "\n",
    "Global methods compute indices considering whole input spaces. Example of global method is Sobol Method (also referred as Variance-based sensitivity analysis) which decomposes output variance into orthogonal terms that are used to compute e.g. first-order sensitivity indices.\n",
    "\n",
    "(https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis)\n",
    "\n",
    "Estimates are computed with generated samples - whole process can be performed with SALib python library.\n",
    "\n",
    "### Local SA\n",
    "\n",
    "In local methods we do SA for given point in input space. To expand it for whole input space we heve e.g. Morris Method (global method) that performs a walk through input space and averages results.\n",
    "Here we will be implementing:\n",
    "- finite differences method (central difference estimator variant):\n",
    "\n",
    "$$\\Delta_C = \\frac{Y(\\theta + h) - Y(\\theta - h)}{2h}$$\n",
    "\n",
    "- adjoint method:\n",
    "\n",
    "$$\\Delta = \\frac{\\partial Y}{\\partial \\beta}(t, \\beta, \\gamma)$$\n",
    "\n",
    "We will perform forward solve and then backpropagate through solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_I(deriv, y0, dt, steps, N, beta, gamma):\n",
    "    _, I, _ = euler(deriv, y0, dt, steps, N, beta, gamma)\n",
    "    return I[:,None,None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = 0.001\n",
    "beta, gamma = tf.constant([0.4]), 1./10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4: Implement finite difference method (central variant)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=1912.5975>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res1 = model_I(deriv, y0, dt, days, N, beta + h, gamma)\n",
    "res2 = model_I(deriv, y0, dt, days, N, beta - h, gamma)\n",
    "\n",
    "(tf.reduce_sum(res1) - tf.reduce_sum(res2)) / (2.0*h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tf.GradientTape() as g:\n",
    "    g.watch(beta)\n",
    "    I = model_I(deriv, y0, dt, days, N, beta, gamma)\n",
    "dI_dbeta = g.gradient(I, beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=1911.5173>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dI_dbeta[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 5: Create plot for deltas for whole beta range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "deltas = []\n",
    "space = np.linspace(0.0, 0.4, 50)\n",
    "for beta in space:\n",
    "    beta = tf.constant([beta], dtype=tf.float32)\n",
    "    with tf.GradientTape() as g:\n",
    "        g.watch(beta)\n",
    "        I = model_I(deriv, y0, dt, days, N, beta, gamma)\n",
    "    dI_dβ = g.gradient(I, beta)[0].numpy()\n",
    "    deltas.append(dI_dβ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'beta')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwz0lEQVR4nO3deXxc1ZXg8d9Rad8lS7JleZGMBXhhtQEHQhIgAdNJx0w6JJAQnIQJnYR0emZ6OiHdPc1kYSbpyTQdZoCEBAdDaAxNZ4K7Y5q4AQMhGCzM4g1jIS+SrNWStVkqqVRn/qhbuCxU1lJVeiXpfD+f+tSr++57dVR26ei+e++7oqoYY4wxo0nxOgBjjDHJy5KEMcaYqCxJGGOMicqShDHGmKgsSRhjjIkq1esA4q2kpEQrKyu9DsMYY6aV1157rV1VS0eWz7gkUVlZSU1NjddhGGPMtCIih0crt8tNxhhjorIkYYwxJipLEsYYY6KyJGGMMSYqSxLGGGOisiRhjDEmKksSxhhjohozSYjIBhFpFZHdI8r/TETeFpE9IvJ3EeXfEZFaEdkvItdElK91ZbUicntEeZWIvOLKHxORdFee4V7Xuv2VcfmJzbT21K4mGo/3ex2GMbPGeFoSDwJrIwtE5ApgHXCeqq4AfuzKlwM3ACvcMfeKiE9EfMA9wLXAcuBGVxfgR8BdqroU6ARuceW3AJ2u/C5Xz8xird0DfO2Rndz7XK3XoRgza4yZJFT1BaBjRPHXgB+qqt/VaXXl64BNqupX1YNALXCxe9Sqap2qDgKbgHUiIsCVwBPu+I3AdRHn2ui2nwCucvXNLLXtnTYAXjvc6XEkxswek+2TOBO43F0Gel5ELnLlFUB9RL0GVxatfA5wXFUDI8pPOZfb3+Xqv4+I3CoiNSJS09bWNskfySS7bftDf4u809JDz8CQx9EYMztMNkmkAsXAGuAvgce9/CtfVe9X1dWqurq09H33pzIzQGA4yIsH2llYnEVQ4c36Lq9DMmZWmGySaAB+rSGvAkGgBGgEFkbUW+DKopUfAwpFJHVEOZHHuP0Frr6ZhXYeOU7PQIA/u6IaEdh5xC45GTMVJpskfgNcASAiZwLpQDuwGbjBjUyqAqqBV4EdQLUbyZROqHN7s6oq8BzwaXfe9cCTbnuze43b/6yrb2ah5/a3kpoirD1nHtVludYvYcwUGfNW4SLyKPARoEREGoA7gA3ABjcsdhBY736B7xGRx4G9QAC4TVWH3Xm+ATwN+IANqrrHvcW3gU0i8gPgdeABV/4A8LCI1BLqOL8hDj+vmaa27W9j1eIi8jPTWLW4iN++1UQwqKSk2FgGYxJpzCShqjdG2XVTlPp3AneOUr4F2DJKeR2h0U8jyweA68eKz8x8zV0D7Gvq5ttrzwbggkVFPPpqPe+29VI9N8/j6IyZ2WzGtUl6z78TGtX0kbNCgxJWLS4CrF/CmKlgScIkvW3725iXn8nZ80KthiUlORRmp7Hz8HFvAzNmFrAkYZLa0HCQ3x9o5yNnlRIeZS0iXLioiNesJWFMwlmSMEnttcOd9PgD711qCrtwUSG1rb10nbBJdcYkkiUJk9S27W8jNUW4bGnJKeUXhvsl6q01YUwiWZIwSW3b/lZWVxaRl5l2Svl5CwpJEXjd5ksYk1CWJEzSaurq5+3mHj5yVtn79uVkpLKsPN/6JYxJMEsSJmk9vz90s8YrRkkSABcuKuKNI8cZDtpEfGMSxZKESVrb9rdRXpDJmXNzR92/anERfYPD7G/umeLIjJk9LEmYpDQ0HOT3tacOfR3pwkU2qc6YRLMkYZJSzaFOev2BUfsjwhYWZ1GSm8FO67w2JmEsSZiktO2dVtJ87x/6Gik0qa7QWhLGJJAlCZOUnt/fxurFxeRmnP4elKsWF3Ho2Anae/1TFJkxs4slCZN0mrsGeLu5hyvOHnuVwfCkutePHE9wVMbMTpYkTNJ5t60XgHMqCsese05FAWk+sUWIjEmQMZOEiGwQkVa3wNDIfX8hIioiJe61iMjdIlIrIm+JyIURddeLyAH3WB9RvkpEdrlj7g6vlS0ixSKy1dXfKiJF8fmRTbJr6R4AYG5+xph1M9N8LJ9fYP0SxiTIeFoSDwJrRxaKyELgauBIRPG1hJYsrQZuBe5zdYsJrWh3CaEFhu6I+KV/H/CViOPC73U78IyqVgPPuNdmFmjtCfUvlOVnjqv+qkVFvNVwnKHhYCLDMmZWGjNJqOoLhJYPHeku4FtA5HTXdcBDGrIdKBSRcuAaYKuqdqhqJ7AVWOv25avqdrf86UPAdRHn2ui2N0aUmxmutdtPdrpvzE7rsAsXFzIwFGRfU3eCIzNm9plUn4SIrAMaVfXNEbsqgPqI1w2u7HTlDaOUA8xV1Sa33QzMPU08t4pIjYjUtLW1TfTHMUmmpWeAueNsRcDJleqsX8KY+JtwkhCRbOCvgL+Nfzijc62MqDfoUdX7VXW1qq4uLR17RIxJbm3dfkrzxu6PCCsvyKK8IJOdNsLJmLibTEviDKAKeFNEDgELgJ0iMg9oBBZG1F3gyk5XvmCUcoAWdzkK99w6iVjNNDTRlgTAivn5HGixezgZE28TThKquktVy1S1UlUrCV0iulBVm4HNwM1ulNMaoMtdMnoauFpEilyH9dXA025ft4iscaOabgaedG+1GQiPglofUW5mMFWltdtP2QRaEgBVJTkcbO8jaHeENSauxjME9lHgZeAsEWkQkVtOU30LUAfUAj8Hvg6gqh3A94Ed7vE9V4ar8wt3zLvAU678h8DHROQA8FH32sxwvf4A/UPDk0gSufgDQZrc8FljTHyMOXxEVW8cY39lxLYCt0WptwHYMEp5DbBylPJjwFVjxWdmlpbu0PDXiV5uqirJAeBgWx8VhVlxj8uY2cpmXJuk0toTaglMtCWxpNQlifbeuMdkzGxmScIklbb3JtJNLEmU5WWQne6jrr0vEWEZM2tZkjBJJXxLjvHOtg4Tkfc6r40x8WNJwiSV1m4/mWkp5I1ztnUkSxLGxJ8lCZNUWnr8zM3PjLpk6eksKcmhvuMEgwG7h5Mx8WJJwiSV1u6BCXdah1WV5hBUONJxIs5RGTN7WZIwSaWtx09Z3sT6I8Iq54RHONklJ2PixZKESSot3QMTHtkU9t5cCRsGa0zcWJIwSaPPH6BvcHjSLYnC7HSKc9KtJWFMHFmSMEkjvNjQeFaki6aqJIe6NksSxsSLJQmTNN6bIzHJlgTYMFhj4s2ShEkarZOcbR2pqiSH1h4/vf5AvMIyZlazJGGSRqtrScyNoSWxxHVeH7LWhDFxYUnCJI3WHj/pqSnkZ018tnVYlbvRn93DyZj4sCRhkkZ4It1kZluHvTdXwjqvjYkLSxImabR0+ye8jsRImWk+KgqzbK6EMXEynpXpNohIq4jsjij7XyLytoi8JSL/T0QKI/Z9R0RqRWS/iFwTUb7WldWKyO0R5VUi8oorf0xE0l15hntd6/ZXxuuHNsmptWfyt+SIZCOcjImf8bQkHgTWjijbCqxU1XOBd4DvAIjIcuAGYIU75l4R8YmID7gHuBZYDtzo6gL8CLhLVZcCnUB4edRbgE5XfperZ2aw1p7YWxLg5kq09xFaKNEYE4sxk4SqvgB0jCj7naqGxxhuBxa47XXAJlX1q+pBQutWX+wetapap6qDwCZgnYQuPl8JPOGO3whcF3GujW77CeAqieVitUlq/YPD9AwEKI1TS6JnIMCxvsE4RGbM7BaPPokvA0+57QqgPmJfgyuLVj4HOB6RcMLlp5zL7e9y9d9HRG4VkRoRqWlra4v5BzJTb7LLlo4mPMLJhsEaE7uYkoSI/DUQAB6JTziTo6r3q+pqVV1dWlrqZShmklq6w7fkiP1yU3iuhA2DNSZ2kx6QLiJfBD4BXKUnL/42Agsjqi1wZUQpPwYUikiqay1E1g+fq0FEUoECV9/MQO+1JGKYbR1WUZhFmk+s89qYOJhUS0JE1gLfAj6pqpErvGwGbnAjk6qAauBVYAdQ7UYypRPq3N7skstzwKfd8euBJyPOtd5tfxp4Vq0ncsZqdS2JWO7bFJbqS2FRcbbNlTAmDsZsSYjIo8BHgBIRaQDuIDSaKQPY6vqSt6vqV1V1j4g8DuwldBnqNlUdduf5BvA04AM2qOoe9xbfBjaJyA+A14EHXPkDwMMiUkuo4/yGOPy8Jkm19AyQ5hOKstPicr6qklxrSRgTB2MmCVW9cZTiB0YpC9e/E7hzlPItwJZRyusIjX4aWT4AXD9WfGZmaOsOrUgXrwFsS0pzeOFAG8GgkpJig+KMmSybcW2SQmuPPy79EWFVJTkMBoIc7eqP2zmNmY0sSZik0NIdn9nWYSeXMrVLTsbEwpKESQqtPf64dFqHLbEkYUxcWJIwnhsYGqarfyimZUtHKs3LICfdZ0uZGhMjSxLGc2098Rv+GiYiVJXajf6MiZUlCeO5eE6ki2TDYI2JnSUJ47mWOE6ki1RVkkND5wn8geG4nteY2cSShPFceG3reLcklpTkEFSo7zgxdmVjzKgsSRjPtfb4SU0RirPT43reyvCN/qzz2phJsyRhPNfS7ac0LyPuM6Or5tgwWGNiZUnCeC5ey5aOVJCdxpycdEsSxsTAkoTxXGu3n7I4rCMxmvBSpsaYybEkYTyXqJYEuCRhfRLGTJolCeMpf2CYzhNDcVmRbjRnzs2jvdfPsV5/Qs5vzExnScJ46uRs68S0JJaV5wOwr6knIec3ZqYbM0mIyAYRaRWR3RFlxSKyVUQOuOciVy4icreI1IrIWyJyYcQx6139AyKyPqJ8lYjscsfcLW5BgWjvYWaW1nCSiPMcibBl5XkA7GvqTsj5jZnpxtOSeBBYO6LsduAZVa0GnnGvAa4ltGRpNXArcB+EfuETWtHuEkILDN0R8Uv/PuArEcetHeM9zAwSz2VLRzMnN4OyvAxLEsZM0phJQlVfILR8aKR1wEa3vRG4LqL8IQ3ZDhSKSDlwDbBVVTtUtRPYCqx1+/JVdbtbv/qhEeca7T3MDJKo+zZFWj4/n72WJIyZlMn2ScxV1Sa33QzMddsVQH1EvQZXdrryhlHKT/ce7yMit4pIjYjUtLW1TeLHMV5p7fbjSxHm5CQuSSwrz+fdtl4GA8GEvYcxM1XMHdeuBaBxiGXS76Gq96vqalVdXVpamshQTJy1dA9QkpuOL4HrUC8rz2doWKlt7U3YexgzU002SbS4S0W451ZX3ggsjKi3wJWdrnzBKOWnew8zg8R7RbrRLLfOa2MmbbJJYjMQHqG0HngyovxmN8ppDdDlLhk9DVwtIkWuw/pq4Gm3r1tE1rhRTTePONdo72FmkNYef1xXpBtN5ZwcMlJTLEkYMwmpY1UQkUeBjwAlItJAaJTSD4HHReQW4DDwGVd9C/BHQC1wAvgSgKp2iMj3gR2u3vdUNdwZ/nVCI6iygKfcg9O8h5lBWrsHOH9hYULfI9WXwlnz8tjXbEnCmIkaM0mo6o1Rdl01Sl0Fbotyng3AhlHKa4CVo5QfG+09zMwxNBzkWN9gwibSRVo2L5+t+1pQVdxUHGPMONiMa+OZdnerjETdkiPSsvI8OvoG35u8Z4wZH0sSxjMnly2dgpaEuz2HzZcwZmIsSRjPhJctnYqWxNnv3cPJkoQxE2FJwnimJcH3bYpUkJVGRWGW3ejPmAmyJGE809zVjy9FKMlNfJKA0CUna0kYMzGWJIxnmroGmJuXkdDZ1pGWl+dR19bLwNDwlLyfMTOBJQnjmeauAcoLs6bs/ZaV5xNUeKfFLjkZM16WJIxnmroGmFeQ+E7rsGXWeW3MhFmSMJ5QVZq6+imfgpFNYYuKs8lJ91nntTETYEnCeKKrf4iBoeCUtiRSUoSz5uXZXAljJsCShPFEU1dojsT8KeyTgJMjnEJ3kDHGjMWShPFEU1c/wJS2JCCUJHoGAjQe75/S9zVmurIkYTwRbkmUT3GSWD4/3Hlt/RLGjIclCeOJ5q4BfCmS8AWHRjp7Xh4iNsLJmPGyJGE8cfT4AGVTOJEuLDs9lco5OZYkjBknSxLGE83d/VPeHxG2rDzPkoQx4xRTkhCR/ywie0Rkt4g8KiKZIlIlIq+ISK2IPCYi6a5uhntd6/ZXRpznO658v4hcE1G+1pXVisjtscRqkktT18CU90eELZuXz6FjJ+j1Bzx5f2Omk0knCRGpAL4JrFbVlYAPuAH4EXCXqi4FOoFb3CG3AJ2u/C5XDxFZ7o5bAawF7hURn4j4gHuAa4HlwI2urpnmVJWm4wOUF0zt8New8Mzr/bacqTFjivVyUyqQJSKpQDbQBFwJPOH2bwSuc9vr3Gvc/qsktI7kOmCTqvpV9SCh9bEvdo9aVa1T1UFgk6trprnu/gD9Q8PetSTmhxcgshFOxoxl0klCVRuBHwNHCCWHLuA14LiqhtvxDUCF264A6t2xAVd/TmT5iGOilb+PiNwqIjUiUtPW1jbZH8lMkaZub+ZIhM0vyCQ/M9X6JYwZh1guNxUR+su+CpgP5BC6XDTlVPV+VV2tqqtLS0u9CMFMwMk5Et5cbhIRW1vCmHGK5XLTR4GDqtqmqkPAr4HLgEJ3+QlgAdDothuBhQBufwFwLLJ8xDHRys0013Tcm4l0kZaV57O/uYdg0G7PYczpxJIkjgBrRCTb9S1cBewFngM+7eqsB55025vda9z+ZzV0A53NwA1u9FMVUA28CuwAqt1oqXRCndubY4jXJInmrn5SBMrypmZFutEsn5/PicFh6tp7PYvBmOkgdewqo1PVV0TkCWAnEABeB+4HfgtsEpEfuLIH3CEPAA+LSC3QQeiXPqq6R0QeJ5RgAsBtqjoMICLfAJ4mNHJqg6rumWy8Jnk0dQ1QlpdJqs+7aTprquYA8FLtMZaW5XkWhzHJbtJJAkBV7wDuGFFcR2hk0si6A8D1Uc5zJ3DnKOVbgC2xxGiSz1QvNjSaRXOyWVSczYsH2ll/aaWnsRiTzGzGtZlyTV39nvZHhF1eXcL2umMMDQe9DsWYpGVJwkyp0Ip03rckIJQkev0B3qg/7nUoxiQtSxJmSvX4A5wYHGa+R8NfI33gjBJSBF58x+bWGBONJQkzpcLDX5OhJVGQlcZ5Cwt5sbbd61CMSVqWJMyUCq9Ilwx9EgCXV5fyZv1xuvqHvA7FmKRkScJMqebwbOspXts6msurSwgqvPyutSaMGY0lCTOljnYNIB5PpIt0/sJCcjNSefGAJQljRmNJwkyp5q5+SnMzSPNwIl2kNF8Ka5bMsSRhTBTJ8U01s4aXiw1F86EzSzjScYLDx/q8DsWYpGNJwkyp5i7vFhuK5oNLSwCsNWHMKCxJmCmVLBPpIlWV5FBRmMXvLUkY8z6WJMyU6RkYotcfSLrLTSLC5dUlvPRuOwG7RYcxp7AkYaZMsg1/jfTB6hJ6BgK81djldSjGJBVLEmbKHO3yfrGhaC47owQRePEdu+RkTCRLEmbKNLvZ1vPyky9JFOWkc05FAb+vtfs4GRMppiQhIoUi8oSIvC0i+0TkAyJSLCJbReSAey5ydUVE7haRWhF5S0QujDjPelf/gIisjyhfJSK73DF3uxXwzDTV5CbSzU3CJAGh2dc7jxynZ8Bu0WFMWKwtiZ8A/6aqZwPnAfuA24FnVLUaeMa9BriW0NKk1cCtwH0AIlJMaOGiSwgtVnRHOLG4Ol+JOG5tjPEaDzUdH6AkN4P01ORswH5waSnDQWV7XYfXoRiTNCb9bRWRAuBDuOVJVXVQVY8D64CNrtpG4Dq3vQ54SEO2A4UiUg5cA2xV1Q5V7QS2AmvdvnxV3e7Wwn4o4lxmGmrqTr6JdJEuXFxIdrqPFw/YJSdjwmL5k64KaAN+KSKvi8gvRCQHmKuqTa5OMzDXbVcA9RHHN7iy05U3jFL+PiJyq4jUiEhNW5t9wZNVc1d/UvZHhGWk+rikqtjmSxgTIZYkkQpcCNynqhcAfZy8tASAawFoDO8xLqp6v6quVtXVpaWliX47M0lNXQPMT8Lhr5Eury6lrr2Phs4TXodiTFKIJUk0AA2q+op7/QShpNHiLhXhnlvd/kZgYcTxC1zZ6coXjFJupqFef4CegUDSzbYe6UNnhv7I+N2eFo8jMSY5TDpJqGozUC8iZ7miq4C9wGYgPEJpPfCk294M3OxGOa0ButxlqaeBq0WkyHVYXw087fZ1i8gaN6rp5ohzmWmmOckWG4pmaVku5y8s5B9fPUKoIWzM7JYa4/F/BjwiIulAHfAlQonncRG5BTgMfMbV3QL8EVALnHB1UdUOEfk+sMPV+56qhoeXfB14EMgCnnIPMw01uYl0ydwnEXbTmsX81396k+11HXzgjDleh2OMp2JKEqr6BrB6lF1XjVJXgduinGcDsGGU8hpgZSwxmuQQXts62fskAD5xbjk/+O1efrX9sCUJM+sl54B1M+OEWxJl+cmxIt3pZKb5uH7VAp7e00xr94DX4RjjKUsSZko0d/dTkptORqrP61DG5fOXLCYQVDbtqB+7sjEzmCUJMyWaknCxodOpLMnh8uoSHn31iN0+3MxqliTMlGg6nnyLDY3lC2sW09Q1wDNvt45d2ZgZypKEmRJNXf1JP/x1pCvPLqO8IJNfbT/sdSjGeMaShEm4Pn+A7oHAtLrcBJDqS+FzFy/ixQPtHGzv8zocYzxhScIkXFMSLzY0ls9evJDUFOEfX7HWhJmdLEmYhAsvWzrd+iQAyvIyuWbFPB6vaWBgaNjrcIyZcpYkTMI1TZNbckRz05rFdPUP8a9vNY1d2ZgZxpKESbhwSyJZV6Qby5olxSwty+Vh68A2s5AlCZNwR7v6Kc5JJzNtekykG0lEuOmSRbxZf5xdDV1eh2PMlLIkYRJuX1MPS8tyvQ4jJp9atYDsdB/3bqv1OhRjppQlCZNQgeEg+5q6OaeiwOtQYpKfmcZXP3wGT+1u5oV3bPVDM3tYkjAJVdvWiz8QnPZJAuDWDy2hck42d2zegz9gI53M7GBJwiRU+Br+yop8jyOJXWaaj++uW8nB9j5+/kKd1+EYMyUsSZiE2nO0m+x0H1Ul07tPIuzDZ5Zy7cp5/N/naqnvsHWwzcwXc5IQEZ+IvC4i/+peV4nIKyJSKyKPuVXrEJEM97rW7a+MOMd3XPl+EbkmonytK6sVkdtjjdVMvV2NXayYn48vRbwOJW7+2yeWIwjf+9e9XodiTMLFoyXx58C+iNc/Au5S1aVAJ3CLK78F6HTld7l6iMhy4AZgBbAWuNclHh9wD3AtsBy40dU108RwUNl7tJsV86d/f0Sk+YVZfPOqarbubeHZt1u8DseYhIopSYjIAuDjwC/cawGuBJ5wVTYC17ntde41bv9Vrv46YJOq+lX1IKE1sC92j1pVrVPVQWCTq2umibq2XvqHhmdEp/VIt3ywiqVludyxeY/drsPMaLG2JP4B+BYQXpVlDnBcVQPudQNQ4bYrgHoAt7/L1X+vfMQx0crfR0RuFZEaEalpa7PhicliV2Oo0/qcBTMvSaSnpvC9dSuo7+jn3m3veh2OMQkz6SQhIp8AWlX1tTjGMymqer+qrlbV1aWlpV6HY5zdjd1kpqWwpCTH61AS4tIzSvjkefP56fPvcshuJW5mqFhaEpcBnxSRQ4QuBV0J/AQoFJFUV2cB0Oi2G4GFAG5/AXAssnzEMdHKzTSxu7GL5eX5pPpm7iC6v/n4MtJ9Kdz+67dsmVMzI03626uq31HVBapaSajj+VlV/TzwHPBpV2098KTb3uxe4/Y/q6rqym9wo5+qgGrgVWAHUO1GS6W799g82XjN1AoGlT1Hu1g5A/sjIpXlZ3LHHy9ne10H3/0XG+1kZp7UsatM2LeBTSLyA+B14AFX/gDwsIjUAh2EfumjqntE5HFgLxAAblPVYQAR+QbwNOADNqjqngTEaxLg4LE++gaHZ3ySALh+9UIOtPZy/wt1LC3LZf2llV6HZEzcxCVJqOo2YJvbriM0MmlknQHg+ijH3wncOUr5FmBLPGI0U2t3uNN6FiQJgG+vPZu6tj6++y97WDwnm4+cVeZ1SMbExcy9WGw8tbuxi/TUlGl/99fx8qUIP7nhfM6al8+f/ePrvNPS43VIxsSFJQmTELsau1hWnk/aDO60HiknI5UH1q8mM93HLRt3cKzX73VIxsRs9nyDzZQJBpU9jd2cMwNu6jdR8wuz+PnNq2nt9vOnD79md4s1054lCRN3RzpO0OMPsHKG3Y5jvM5fWMj//sx51Bzu5NtPvMVwUL0OyZhJS8ToJjPLhWdaz4aRTdF84tz5HGrv48e/e4fB4SB//5nzp+3yrWZ2syRh4m730S7SfSmcOTfP61A89Y0rq8lI9XHnln0c632Vn69fTX5mmtdhGTMhdrnJxN3uxi7OmpdHeqr99/rKh5bwkxvOZ+eRTj7z05dp6R7wOiRjJsS+xSauVJXdjd2z+lLTSOvOr2DDFy+ivuMEn7r3D9S29nodkjHjZknCxFVDZz9d/UMzYrnSeLq8upRNt34Af2CY63/6B3Ye6fQ6JGPGxZKEiatds2ym9UScs6CAf/7apeRnpfG5n2/nsR1HCN2+zJjkZUnCxNXuxi5SU4Sz5s3uTutoFs/J4YmvXsoFC4v49j/v4uuP7OT4iUGvwzImKksSJq52NXZx5tw8MlJtuGc0pXkZ/Oo/XsLt157Nv+9rYe0/vMgfatu9DsuYUVmSMHGjquw52m2XmsbBlyJ89cNn8OuvXUZ2uo/PP/AK//OpfQwGbE0Kk1wsSZi4Odo1QEffICtn4HKliXLOggL+9Zsf5MaLF/Gz5+v41H0vsa+p2+uwjHmPJQkTN7sa3Ezr+TayaSKy01P5H//hHH72hVUcPT7Ax+9+kb/5zS46+6yvwngvljWuF4rIcyKyV0T2iMifu/JiEdkqIgfcc5ErFxG5W0RqReQtEbkw4lzrXf0DIrI+onyViOxyx9wtIhLLD2sSa8/RLnwpwrJySxKTcc2KeTz7Fx/mC2sW8+ir9Xzkx9vY+IdDtiyq8VQsLYkA8BequhxYA9wmIsuB24FnVLUaeMa9BriW0NKk1cCtwH0QSirAHcAlhBYruiOcWFydr0QctzaGeE2C7Wrsoros1+5RFIPC7HS+u24lW755OSvm53PH5j18/O7fW8e28Uwsa1w3qepOt90D7AMqgHXARldtI3Cd214HPKQh24FCESkHrgG2qmqHqnYCW4G1bl++qm53a2E/FHEuk2Q6+wZ5+d1jXFRZ7HUoM8JZ8/J45D9ewk9vWkXfYIDP/eIVvvzgDt6oP+51aGaWiUufhIhUAhcArwBzVbXJ7WoG5rrtCqA+4rAGV3a68oZRyk0SenTHEfyBIDetWex1KDOGiLB25Tz+/b98mL+85ix2Hunkunte4gsPvMKrBzu8Ds/MEjEnCRHJBf4Z+E+qesqwDNcCSPiUUhG5VURqRKSmra0t0W9nRhgaDvLwy4e5bOkcm0SXAJlpPm67YikvfftKvnPt2exr6uYzP3uZz/zsZX5/oN1mbZuEiilJiEgaoQTxiKr+2hW3uEtFuOdWV94ILIw4fIErO135glHK30dV71fV1aq6urS0NJYfyUzCv+1upqlrgC9fVuV1KDNaTkYqf/rhM3jxW1fyt59YzuFjfdz0wCusu+cl/qmmnoEhWwXPxF8so5sEeADYp6p/H7FrMxAeobQeeDKi/GY3ymkN0OUuSz0NXC0iRa7D+mrgabevW0TWuPe6OeJcJon88qWDVM7J5oqzyrwOZVbISvfx5Q9W8cK3ruAH162kzx/gL594i0v+xzPc+du9HGrv8zpEM4PEsujQZcAXgF0i8oYr+yvgh8DjInILcBj4jNu3BfgjoBY4AXwJQFU7ROT7wA5X73uqGr7g+nXgQSALeMo9TBJ5o/44O48c544/Xk5Kio1QnkoZqT5uWrOYz1+yiJfrjvGr7YfZ8NIhfv7iQT50Zimfv2QRV5xVZut6mJjITLueuXr1aq2pqfE6jFnjzze9zrP7Wnn5r64iN8MWOvRaS/cAj756hEdfPUJLt5/C7DQ+fk45111QwapFRZbITVQi8pqqrh5Zbt9qM2kt3QP89q0mbv5ApSWIJDE3P5P/9NEzue2KpbzwThu/eeMo/7yzgUdeOUJFYRbrzp/PdRdUzPqlZc342TfbTNqvth9mWJUvXlrpdShmhDRfClctm8tVy+bS6w/wuz3N/OaNo/z0+Xe5d9u7LCnJ4WMr5nL18rmcv7AIn7UwTBSWJMykDAwN88grR/josrksmpPtdTjmNHIzUvnUhQv41IULaOvx89TuJn63p4UHXjzIz56voyQ3navOnsvHls/l0qVzyE63XwvmJPvfYCZl8xtH6egb5EuXVXodipmA0rwMbv5AJTd/oJKu/iG27W9l694Wtuxq4rGaetJ9KVy4uJDLq0v54NISVlYUWCtjlrOOazNhqsq1P3kRgKf+/HLsvovT32AgyCsHj/HigXZePND+3u3KC7LSuPSMOVx6xhwuqirmzLI86/yeoazj2sTNy3XHeLu5h7/7k3MtQcwQ6akpXF5dyuXVocmobT1+/vBuKGH8/kA7T+1uBiA/M5XVlcVcVFnMRZVFnLOgwFYhnOEsSZgJ++VLhyjOSeeT58/3OhSTIKV5Gaw7v4J151egqtR39PPqoQ5qDnXw6qEOnn07dCOFdF8Ky8rzOHdBIecuKOD8hYUsKc21S1QziCUJMyH/+MoRtu5t4ZtXLrVbgs8SIsKiOdksmpPNp1eF7pRzrNfPjkOdvF7fyZv1x/l/rzfy8PbDAOSk+1hRUcDy8nyWl+ezrDyf6rl2C/npypKEGbeHXz7Ef3tyD1eeXcbXr1jqdTjGQ3NyM1i7ch5rV84DIBhU6tp7eaO+i7cajvNWQxeP7ain391PypcinFGaw7LyfM6cm8fSslyqy3JZVJxNqs9mhCcz67g24/LLlw7y3X/Zy8eWz+X/fu4Cuw5txjQcVA4f62NfUw/7mrrZ19TN2809NB7vf69Oui+FqpIcls7N5YySHCrdo2pODkU56R5GP/tYx7WZtJ+/UMedW/axdsU87r7xArsXkBkXX4qwpDSXJaW5fPzc8vfKe/0B3m3tpba1lwOtvdS29rC7sYundjURjPibtSArjcqSHBYXZ7OoOJuFxVksLMpmYXE25QWZ1gKZIpYkzGndt+1dfvRvb/Pxc8v5h8+eT5p9MU2McjNSOW9hIectLDylfDAQpL7zBIfa+zjY3sehY6Hn1+s7+e2uJoYjMogvRZiXn8n8wkzmF2ZRXpAV2i7Iorwwk3n5mRRlp9tw3TiwJGFGparc81wtP/7dO6w7fz7/+/rz7C83k1DpqSmcUZrLGaW579sXGA7S1DVAfecJ6jtOUN/RT0PnCY52DbDzSCfNXU0MDZ966TzNJ5TlZTI3P4N5BZmU5WVSmpcReuSGnktyM5iTm25//JyGJQlzioGhYZ58o5EH/3CYfU3dfOqCCv7X9efZkEbjqVRfCguLQ5eaOOP9+4NBpb3Xz9GuAZqO99PSPUBzt5/W7gGauwfY39zDC++00+sPjHr+wuw0inPSmZOTTnFOOsU5GczJSacoJ52i7DSKstMpdM9F2enkZabOmlaKJQkDQEPnCR7efpjHdtRz/MQQZ8/L44efOofrVy+0BGGSXkqKUJafSVl+JuePuIwVqX9wmPZeP229ftp6Tj46+gbp6BvkWJ+fg+191BzqpPPE4Cl9JKe8n0B+VhoFozzyMtPIz0oNPWemku9e52akkZuZSm56KjkZvmnTMrckMUsNDQc50NLL7qNdPLOvha17WwC4ZsU81l9aySVVxTab2sw4Wem+ky2SMQwHle7+ITpPDNJ5YojjEc/HTwzR1X/qo7Gzn67+IboHht536WvUWNJ85GSkkpsRes5xySO8nZ3hIzvdR3Z6qns+uZ2V5iPTPWel+chK95GZ5iM3IzXuf9QlfZIQkbXATwAf8AtV/aHHIU0bqkqvP0Brj5/W7tBfSLsau9hztIu3m3oYHA4CUJyTzp9++AxuWrOYisIsj6M2Jjn4UiR0uWmCQ3FVFX8gSPfAEN39AXoGhugeCNDnD9A7EKDHPff6h+j1B+j1D3PCH6BvMEB77yCHO07Q5w9wwj/MiaHhUzrsx/LLL17EFWfHdxnhpE4SIuID7gE+BjQAO0Rks6ru9Tay2Kgqw0FlWJVgEALBIMEgDAWDBIaVoeEggaASGA4yOBxkMBDEHzj57A8M4x8K0jcYoNd/8j9fr3+YXv8QHX2D7yWG8GSmsPzMVFZWFPClyypZUVHAyvn5VM7JmTXXV41JNBEhMy30l31ZjGs7qSqDw0H6B4fpGxymfzDAicFh+geH6R8aZmBoOPR6KFS2tOz9nf6xSuokAVwM1KpqHYCIbALWAXFPEnc/c4An32hEAVziVkL/SKFnCKoSnnuoqgQVFPes4bJQ/WAwVHfYlQX1ZHKYwB8G45KaIq7ZGmquFuekc96CQsryMijLz6AsL5OyvAwWFmezoCjLLiMZM02ICBmpPjJSfRR6tGxLsieJCqA+4nUDcMnISiJyK3ArwKJFiyb1RnPzMzh7Xj4ISOicoXMDIpAigrgCQUiRiHIJ1U9x+8LlKa7MlyLv7U8RwZcy4iFCSoqQ5hNSU1JI9Z3cTvMJGWk+MlJT3MP33nN2RugaZEZqiv3iN8YkRLIniXFR1fuB+yF0W47JnOOzFy3isxdNLsEYY8xMlexjsBqBhRGvF7gyY4wxUyDZk8QOoFpEqkQkHbgB2OxxTMYYM2sk9eUmVQ2IyDeApwkNgd2gqns8DssYY2aNpE4SAKq6BdjidRzGGDMbJfvlJmOMMR6yJGGMMSYqSxLGGGOisiRhjDEmqhm3xrWItAGHJ3l4CdAex3DixeKaGItrYiyuiUnWuCC22BaraunIwhmXJGIhIjWjLQTuNYtrYiyuibG4JiZZ44LExGaXm4wxxkRlScIYY0xUliROdb/XAURhcU2MxTUxFtfEJGtckIDYrE/CGGNMVNaSMMYYE5UlCWOMMVHNmiQhImtFZL+I1IrI7aPszxCRx9z+V0SkMmLfd1z5fhG5JhniEpFKEekXkTfc46dTHNeHRGSniARE5NMj9q0XkQPusT6J4hqO+Lziesv5ccT1X0Rkr4i8JSLPiMjiiH1efl6ni8vLz+urIrLLvffvRWR5xD4vv4+jxuX19zGi3p+IiIrI6oiy2D4vVZ3xD0K3GX8XWAKkA28Cy0fU+TrwU7d9A/CY217u6mcAVe48viSIqxLY7eHnVQmcCzwEfDqivBioc89FbrvI67jcvl4PP68rgGy3/bWIf0evP69R40qCzys/YvuTwL+5ba+/j9Hi8vT76OrlAS8A24HV8fq8ZktL4mKgVlXrVHUQ2ASsG1FnHbDRbT8BXCUi4so3qapfVQ8Cte58XseVSGPGpaqHVPUtIDji2GuAraraoaqdwFZgbRLElUjjies5VT3hXm4ntMoieP95RYsrkcYTV3fEyxwgPMLG0+/jaeJKpPH8ngD4PvAjYCCiLObPa7YkiQqgPuJ1gysbtY6qBoAuYM44j/UiLoAqEXldRJ4XkcvjFNN440rEsYk+d6aI1IjIdhG5Lk4xTSauW4CnJnnsVMUFHn9eInKbiLwL/B3wzYkc60Fc4OH3UUQuBBaq6m8neuxYkn7RIRNVE7BIVY+JyCrgNyKyYsRfOuZUi1W1UUSWAM+KyC5VfXcqAxCRm4DVwIen8n3HEiUuTz8vVb0HuEdEPgf8DRDX/prJihKXZ99HEUkB/h74YiLOP1taEo3AwojXC1zZqHVEJBUoAI6N89gpj8s1H48BqOprhK41njmFcSXi2ISeW1Ub3XMdsA24YCrjEpGPAn8NfFJV/RM51oO4PP+8ImwCrpvksVMSl8ffxzxgJbBNRA4Ba4DNrvM69s8rER0tyfYg1GKqI9RxE+74WTGizm2c2kH8uNtewakdP3XEr6MslrhKw3EQ6tBqBIqnKq6Iug/y/o7rg4Q6YYvcdjLEVQRkuO0S4ACjdP4l8N/xAkK/OKpHlHv6eZ0mLq8/r+qI7T8Gaty219/HaHElxffR1d/GyY7rmD+vmH+A6fIA/gh4x30h/tqVfY/QX08AmcA/EerYeRVYEnHsX7vj9gPXJkNcwJ8Ae4A3gJ3AH09xXBcRur7ZR6jFtSfi2C+7eGuBLyVDXMClwC73hdkF3DLFcf070OL+vd4ANifJ5zVqXEnwef0k4v/3c0T8UvT4+zhqXF5/H0fU3YZLEvH4vOy2HMYYY6KaLX0SxhhjJsGShDHGmKgsSRhjjInKkoQxxpioLEkYY4yJypKEMTFydwDdPYH6XxSR+YmMyZh4sSRhzNT7ImBJwkwLliSMiY9UEXlERPaJyBMiki0iq9zN3l4TkadFpNytcbEaeMStO5AlIn8rIjtEZLeI3D8Fd/k1ZtxsMp0xMZLQQlAHgQ+q6ksisgHYB/wHYJ2qtonIZ4FrVPXLIrIN+K+qWuOOL1bVDrf9MKFbr/yLFz+LMSPZXWCNiY96VX3Jbf8K+CtCN13b6hoGPkJ3Ch3NFSLyLSCb0L2c9gCWJExSsCRhTHyMbJL3ELpv1AdOd5CIZAL3ErrXTr2I/HdC9+syJilYn4Qx8bFIRMIJ4XOEVnkrDZeJSJqIrHD7ewjd3hlOJoR2EckFTlmX2xivWZIwJj72A7eJyD5Ct9n+P4R+4f9IRN4kdHfQS13dB4GfisgbgB/4ObAbeBrYMaVRGzMG67g2xhgTlbUkjDHGRGVJwhhjTFSWJIwxxkRlScIYY0xUliSMMcZEZUnCGGNMVJYkjDHGRPX/AedqbPESJtR1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(space, deltas)\n",
    "plt.xlabel('beta')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SALib\n",
    "\n",
    "`SALib` provides us with global sensitivity methods. Here we will use Sobol Analysis method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from SALib.sample import saltelli\n",
    "from SALib.analyze import sobol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = {\n",
    "    'num_vars': 2,\n",
    "    'names': ['beta', 'gamma'],\n",
    "    'bounds': [[0.0, 0.4], [1./18, 1./8]]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "param_values = saltelli.sample(problem, 100, calc_second_order=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x161120a60>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6q0lEQVR4nO2dfZBf5XXfP0fLGq/sRMJYmakXZOSGwEBlo1i8TGmcoNSG1GOkYhxhOylOPSVpwnQgRK089hAgaZGjSbBnwkzC2KSxcY0wSTRyIVFdCzetYhwtlmQijBwZMGidaRRgSWOWsJJO/9jfFXfv3pfn3vu83bv3O8Owur/7cp7nOffcc77Pec4jqsqAAQMGDOgvloUWYMCAAQMGuMVg6AcMGDCg5xgM/YABAwb0HIOhHzBgwICeYzD0AwYMGNBznBZagCze/OY36znnnBNajAEDBgzoFB577LG/U9VVeb9FZ+jPOeccpqamQosxYMCAAZ2CiHyv6LeBuhkwYMCAnmMw9AMGDBjQcwyGfsCAAQN6jsHQDxgwYEDPMRj6AQMGDOg5osu6GVAPO/dPs333Yb4/M8tbVk6w5crz2LRuMrRYvcTQ192H6zGMVUcGQ99h7Nw/zcf++HFm504AMD0zy8f++HGAKJSrTxj6uvtwPYYx68hA3XQY23cfPqVUCWbnTrB99+FAEvUXQ193H67HMGYdGTz6DuP7M7O1jseOWMNe6F5f2+zLmMclizJZXY9hzDoyGPrIUOelesvKCaZzlOgtKydci2kdMYe90K2+ttmXsY9LGlWyuh7DqvuH/GAO1E1ESBR1emYW5TVF3bl/Ovf8LVeex8T42IJjE+NjbLnyPA/S2kXMYS90q69t9mXs45JGlayux7Ds/nXfbdsYPPqIUKaoeV/+5FhXwuoyxBz2Qvu+znpzV5y/ikeePOZk3Gz2ZezjkkaVrK7fl7L7X75tT6132zYGQx8RmrxUm9ZNLlKULnGqCdqE1b7am9fXJsijFO579NlTv9umQ2xSFEX3UuDybXsW9HVovTNpd9MxNEXR/UN/MAfqJiIUvYh1XtDQIWJTNA2ru9DevEgtC5t0iE2KIu9eCdJ9HcM4xEyv2Xi322Aw9BHBhqJ2iVNNY9O6Se68Zi2TKycQYHLlBHdes7bS++pCe029NlveXdO+rLpXHpK+jmEcbLbbNkJ/hJY8dRM63EzDBofYJERs2we2+rBJWF3V3hjGt4hSyDvPFmxSFMm91mx9CM35vUy3mupd03FzTc00Rej5tCVt6GNMHWurqHX52bZ9ELoPy9obWrYEW648b4EceYiFYihDlW7Z0jsginGzjZAfISPqRkSuEpHDInJERLbm/P4uEfmmiBwXkWtTxy8Ska+LyCER+ZaIbLYpfFvEEG7aRt0QsW0fhO7DsvaGli1BHqXwc5etjpJiKENZX9fRu537p7nlgYOFYxPLuNnEzv3TXL5tD2u2PsTl2/Z4n0Oq9OhFZAy4G3g3cBTYJyK7VPWJ1GnPAh8Bfi1z+cvAv1HVvxaRtwCPichuVZ2xIXxb+J4J90Ej1A0R2/ZB0XkmVIUNFLW3TIYQqYGxUgp1YKJbVXr3iZ2P84VHn82lgKA5DRQzYogsTaibS4AjqvoUgIjcD2wEThl6VX1m9NvJ9IWq+p3U398Xkb8FVgEzbQW3AZ+rHX0Odh2j0rYPiq4X5tvsQ5Gz7U36uggxrmbtCsp0q0rvdu6fLjXy0IwGagKfczd118e4gAl1Mwk8l/r30dGxWhCRS4DXAd/N+e0GEZkSkaljx47VvXVj+JwJjzUcbdsHW648D8k5rhCsbWXpjF3gwvuK7bsPlxp5gCvOX+X8vfSdCho6hx48pVeKyD8BPg/8gqqezP6uqveo6npVXb9q1SpncmR5MsBbOlYMg52Htilpm9ZNNgrD66Auv1n23NBceFOuNjTHawMm+vDIk8ecp0n6drpC59CDGXUzDZyd+vdZo2NGEJEfBh4CPq6qj9YTzx6KqJM7r1nL3q0bnD8/5qJYbfnjSYdta0J5FfX15MqJ4Ea+CX0XA8drAyZppulyBa7a5tvpysu68h1Zmnj0+4BzRWSNiLwOuA7YZXLz0fl/AnxOVR9sLmZ7hKZOQi+YcAmXbWsybq7kaetVN9XB0LprC2WrbBP4cHyqPGzb0VMMC7kqPXpVPS4iNwK7gTHgXlU9JCJ3AFOquktELmbeoJ8BvE9EblfVC4GfBd4FnCkiHxnd8iOqesBBWwqxc/908AyM0AsmsmgzGZV37Z3XrG10vyo5mtb/Abt9bcOrbupJxkr7JTDVpfS4TM/MIrCA9vPl+JR52HnjfNOOA9y26xC3XX1hYx0KnXVltGBKVR8GHs4cuzX19z7mKZ3sdfcB97WUsRViysAIPdgJ2hgtmxSYiRxNKS/bfW0jc6JpW2Km/erqUnpcQq1arltlEmBmdq6TdFmC3te6GTIwFqMNFeC71nkslJcNr7ppW2Lpgzy00YdN6ybZu3UDT297L3u3bvBqQIueXTaeXaTLEvS+BILtDIwYaqe0RRuj5bvWeUjKKz3Wy0Q4oYvzi+p41U3bUnVdSJ2MnVYCO7u2JahqV6z2ofeG3mYGRt+zH0yMlo9a59l7haC8smOdZ+RhPu+7Dpq2pei60DoZM60E9funqi5RWbtCj0UZek/d2Ax7+5z9YNonrmudx0xJ5OGRJ/0t8MtDaJ2MeQyhfv8kGTJnLB9f9FtVu0KPRRl679HbDP27EKaaoE2f2OzP2DKR0vBZQ75NuB9aJ2MeQ2i3a1vdcQk9FmXovaEHe6F/7GGqCbLKe9fmi2r3jU0qJZZMpCx81ZAvSue7/cuH+PX3VafzrZgYZ2Z2Lve4L8Q6htDuna3brpjtQ++pG5uIOUw1WeThssZHH5bop2GyuMfG2BdRRC++PGc0NpJXaKjk+FKDz3c2ZvuwJDx6W4g1TDWdBDLJBW9CI8Q8CWWCsjanj19x/ioeefKY1bE3Secre8bMy4u9+bLjSw0+39lY7QMMhr42YgxTTRfzFFERyfGmBtv0+TGmnlW12bV8VRRRFX0UM13QBC50xOc7G6N9gIG66QVMJ4HGCuL55HjTrAGT5/suDWuK0JkSVRRR0ZiVXe+SLnBJ0VXpSN/oQZ8YPPoewNSrK8oFT443zRoweX4Mmy/kIXSmRNL2m3YcyP29aMyy19vygqs27HZJ0VV9dGOhB2OMTKuwZA29jcHK3sMFh2sC0zKoReWEJ0cGuSkNUPX8GIrKFSEG6mPTuslThb6ymHSQHVKEKkPu+mNd9tGNxVHo6nzUkqRubNAIefe479Fng1ATpmVQq8L8pjRA2fNjKiqXh1gyJWKQo8qjdh39lJUPDh15JQhN9TXFkvTobXgHJisnfXocJl5dVZjfdiFV3nmxF5WLJVMiBjmqjKnr6KcsMiyKeHw7CrF8cOpiSRp6G4PleuXkzv3T3P7lQ7w4SpNbOTHeqh52gqoPgu2sgZi39UsQS6aEDTnyKEkw+4BUGXLXOyVVfey2PHiQuROvzVmMj4l3R6Goj5aJsGbrQ9Fy9kvS0NvwTFyunNy5f3qRUs/MzrHlSweBuLnALGLd1q+PyOOPt3zpIAindKmMU64y5D6ijtKPXXZeumqncQcoKnqWTJrHytkvSY7eBh/qcuXk9t2HFxj5BHMnNXouMIsYuOelgjyabO6kLtKlIk7ZZK5n07owNeS37z7M3MmF7ch7H1ynYGb7KC/9NUbOfkl69DY8E5crJ8voDhdcoMt0MRt9bXvbw5Celkt5bFCPsdBYWdRZq+E6IybdR2u2PlRL3lBYkoYe7Ci0q5eijBayPfnk4+Vo008utj00ubbofm0/WC772pROTM7tEnyu1bCxUUls/bskqZvYseXK8xgfWxwSji+zP/kUe7pYLNse2kjJdd3XJnQidJM6M6EAbSRZ1B3nrlCTRoZeRK4SkcMickREtub8/i4R+aaIHBeRazO//ZmIzIjIf7cldN+xad0k2699x4LND1ZOjLP9A++wHkEUvQTTM7NRLDH3ve1hEcdrw0i7Ts1L+OOysgljIqfkjmF8TWEyf1CWh2+KOuOceP6zcydO9XnRGpbQqKRuRGQMuBt4N3AU2Cciu1T1idRpzwIfAX4t5xbbgeXAL7aWdgnBF1daFu7HkD3gc9vDMmrFhpH2EeYnY5XNDBlfJsbZN7Gi6p2wkf5pOs55W00mz4qxP008+kuAI6r6lKq+CtwPbEyfoKrPqOq3gJPZi1X1q8D/syHsAPsoC/dn505w044D3gtIpb3ql189Pm+kUnC17WGZN2fDW/QV5ud5v298/WnG2TddhekK8TKYjnPslGcWJpOxk8BzqX8fBS61KYSI3ADcALB69Wqbtx4wQtEEU1VRLfDr/WXXELz48hxjy4SVE+O8NDvndNvDMm/urs0XtfYWfa5+zXq/MWeH2MxESrc7ue/NOw4Y33fLleeVLsxK7hlr7aYiRJF1o6r3APcArF+/PsAyiH7DpOZ6mfKCv3IOt3/50CLP88RJRQSe3vbe2verQ4GVUSu2jHSo9MVYs0NcZSK1um/BwqzsPfMQuj+LYELdTANnp/591ujYgI7AJMw0ydjw4a28WLAzUtFxm6iiVkItFrKBWLNDXFEgTe9btjCrqr5VDP1ZBBOPfh9wroisYd7AXwd8yKlUDeBqIYqN3OnQC3ZMJpjSHquvHP7QyBubO69ZG3y8qtBEp2IompYHV5lITe/b9LrJyG1DpaFX1eMiciOwGxgD7lXVQyJyBzClqrtE5GLgT4AzgPeJyO2qeiGAiPxv4HzgjSJyFPioqu622Ygowz+HctWFadie0Ap5Iaovb2XlxDgzs4u995UT4zlnN0fR2Nx5zVr2bt1g9Vk20UanYlz16opSanrfquuK6jbV1RnftsEoj15VH1bVH1PVf6qq/3l07FZV3TX6e5+qnqWqb1DVMxMjP/rtJ1R1lapOjM6xauQhvvDPtVx1UTdst5G90BS3XX3hoiyb8WXCbVdfWHBFM8QyNnXRVbmL4IpSanrfsutsyup7HKOYjG2L2MI/W9fbQlFdnrKMhFDeny+KIZaxqQufcvugFlyNd9P7mlxnQ1bf+tcLQx9b+OdaribIpp3FQCkVwcdHJqaxqYOyeug7909b6zefOuJqvJvet+w603tWfSR9618vat3EFv65lqst+hb+N0GsY1OFouyoE6pWt64cdKQ5TOrl+Na/Xnj0IcI/k7B2qWU6dAmxjk0VEvlueeDgqc0uEvjaqLuvsEVVmVTR9K1/vTD04Df8Mwlrs0pz1+aLvBmR2MJG16j7goYcGxvYtG6SmwtWMtvcqLtPOlIFm1SV6UfS5zxYL6gb36gKa22UtG2KGMNGl6jb1yHHxiZs1N4pQ590xAQ2qSrXY9MEg6FvgKovdkh+0+TZttMnXW/fVoa6fd0X7tm1IQ6ZYgv+dcomVRXjR7I31I1PVIW1bZSmLU/oO2wMncFTt6/rlKGNmb/3wfGGSrENoVM2qarsKvNkD4Dbv3yI23Ydql2czwYGj74Bqr7YTUM3G7SC77AxtIdct70mx7tC73S59k4ZQuiUbS9807rJU/dMJs1ffHmOmdm5IDo1GPoGqAprmyqNDQX3HTaGzs6o216T80N/vJY6QuiUC6qqqgiaT53qJXXja0Vf2aIKqB9W21Bw32lbobMz6rbX5PymWxCW3TP7+xXnr+KRJ49FRw3FQFmF0inbVFWb7S5to3eGPjRnnKCJ0thScJ/cqo3t29qibnurzre5BWFekbjpmVnue/TZU9fHsjI5lnfnivNXLeif9PEuoWybzvQ5PtA76qbLYbct2sVnxkLo7AwXqDMOO/dPc8sDB0t1riqEz54fCrG8O488eazW8VhRtceDT4eodx59USg0PTNrtRaICxTRCgCXb9tjFE6H8MpcRRChaARTOijp6+wK1QSJLrYtgucLoedbYpOjLbJ6tHL5OKoEybrpnaEvC5diCI+rkDWadQ23yfLrLiA0jWDy8ary1JOw3CSEB1hhueZ+XYSeb4lNDhsIlaKaRe+om7JwaXbuBLc8cDC6NLky1A2n++INxUIjlKGsT9Nhuck2jQAilacUwgZdF8tCn1jk6BN659EnX8+bCmqBJFX+0ufGjLqGuy/eUBc+WEV9PSayYJ4iG8LnEz0w03Bf3LzoZ8uDB2svzglV6K1qS8cVE+OIwM07DrB992FjmWLIIIoFvfPoYV5hJ0sMW2yeYRnqLvzpizcUY72QLIr6+rd/9h25lUyTxU1Futm0bXnRz9wJbbQ4x/cirKLFaQB7t27grs0X8Y/HT/Liy/Xa0pVFb77QS0MP1eFyDJ6hSbjdpW0AbaILH6ymfW27bSa6HMq5qdLxKoquKYXXBerPJ3pH3SQoq9sN7jxD03AxL9y+eccBbtpxYMGO8k3CaRsTQKHD3hjrxRf1SZVMVdRE27aZTvb6dm5MJtSrKLqmFF7M1F+Id8vI0IvIVcCngTHgM6q6LfP7u4BPAW8HrlPVB1O/XQ98YvTP31TVP7QgtxGSzvO1oKdOpkiex5F8jrLX+Z65D53xkiCWjAVo3idF1915zVr2bt1gRba8RWt58E17mWSAVc0pNZ1zinWuKtS7VUndiMgYcDfwM8AFwAdF5ILMac8CHwH+W+baNwG/DlwKXAL8uoic0V5sc/ikMuqEi1WeRcgwcwh7FyNmCiGr42csH2d82cIUnhC0l4lXXUVjNaW5YqX+Qr1bJh79JcARVX0KQETuBzYCTyQnqOozo99OZq69EviKqr4w+v0rwFXAF1tLXgO+PMOi8DlP4U3CbZthZp1w0XbY23YHqNCUDcRPIeStvwjdhyZedRVF15TCi5H6g3CUkomhnwSeS/37KPMeugnyrl3U0yJyA3ADwOrVqw1vHRd27p9GIDd1Li9cNAm3bYWZuel3XzrI7V8+xMzLi9PvbIa9dUPVWGijLLpGIcQwT2NaB6lK1qZtiYn6SxBKH6LIulHVe1R1vaquX7WqW4WLEmzffTjXyAvkhovpcDs5Lw2bYWZu+t1JLUxZsxn29mUHqL5RCFWwkZ7YlwwwmwilDyYe/TRwdurfZ42OmWAa+KnMtV8zvLZTKAq9lGJPNO1xuAy166TfNc30qftsW8d9oW8UQhVsldJw4VXHQEs1RSh9MDH0+4BzRWQN84b7OuBDhvffDfyX1ATse4CP1ZayAygKycoWbqXhMsxskn7XVp7kZSxaBVoUqq6YGGdmdvEK0ZB1YLKG5a7NFwU3dq4R6wc3VmqvDkLoQyV1o6rHgRuZN9rfBh5Q1UMicoeIXA0gIheLyFHgA8Dvi8ih0bUvAL/B/MdiH3BHMjEbM5rUDYk5RDettWJ7TqDo41LWL0X1XtrUgWmDpbrCMtaVybFSe0Xwvcl5EYzy6FX1YeDhzLFbU3/vY56Wybv2XuDeFjJ6RVOPIeYQPSvbiolxfvDqceZOvOZvu54TSDBZ0S9F9V6a1oFpi75UA62LGDaUyUOskUYeYoo+ersyFppxeW1e7JhDdJ/pd0UvnUDlIqHYFrq0MSw2+jjWmvyh5IpNP8oQk5PQW0Pf9GvaJY+hDVx+lNrw7D49SRNj1bQtJtsLmm5sEltN/pByxRpp5CEmW9JbQ9/0a9oljyFWtOHZbVNgRQbV1Fg1bUsVl2zybB8eoe+oty1ipkiziMmW9NbQN/2aVnkMXU7tagvTtrfl2W1FG2XG3NRYNW1Lmf6ZPtu1R9jVqLdMP/J0FMJ8GGKKPnpr6Jt+Tcs8htChdEjUaXssnkyZQTU1Vi5WxLp+tin6FvUWrQBHOJV44POdjSn6iGJlrAu0SXfctC5/84WupXbZRJ22x5JqWmZQTdMHXayIdf1sU7SJemMY3yyKVoCns8vA7ztbZEt8o7cefZuvaRFFETpkDUkb1Wl7KE8m2z9FE6mJPKZ1WJq0peo6l882hYuoNyTqvIfpc7N6c8X5q3jkyWNRZRu1hWjOphwhsX79ep2amgr2/Gz4B/MvYbJRRNHqV1u1xZvI5UPRLt+2J1jbTZDXP+NjAjrv1SVI91nIlzYGgxFap2yjSEfzkOhtXh9kkfQJ5H+gY+kvEXlMVdfn/jYY+oUoM2hFXqCPgQ5taGM3CkX9c8bycZa/7rTOeWC+EMMHxxZyP/bLZAFHDwv11vTjkJQyidnZKTP0vaRu2ihvGUWRDlmnZ2YZE1nA97l8QULTRrGG6wmK+mHm5Tn23/oez9KEh+k7kLeQ7vJte6Ia4zptgcU6mncsOdf0/Sk7rwtrbHpn6NtmxlTxlnkcq4+Z/BgyHWJe+RtD/8QC29seVl3nEnVlKtLRIvlNC/69pcSj74KO9S7rpm1mjElGga3smzoFj2LNdIgFQ/+8hpi2PWxb1Mt1pptJwb9Ej7qsY73z6NtSHCYUhQ0apYmnUiVX0XNipVvaIK9dyYR5yEyfGPq3iX7u3D9daytME9iIEFxTlokcN+04UHhOdh4qtvE2Qe8MvY0QvoqisPGMJotV6lInMYbiNlDUrjuvWet1UizW/q2rn0k7yu7XBDZKJfig5DatmyzNqEvLGjN9WYbeUTc+wqs6zygKXX1MrvZ1gVcs7SqS45YHDgatV1/3HSgrK93m3bGh477oki7TMibonUdfVcLARthlSqOUeXxVnooNWUNn6riCr3ZVjUHR806oWvfsy2SporFWLh9HFW7ecYDtuw8btwMW0xZ1YCu6Bvd0SexZZW2xZPLoQ+SBN83JBzsLM0Ln3ruCj3aZ6EtVDrYtecpkgXJdadOOtvLHvvaibyjLo+8ddVOEEOF+VU7+ndesZXLlBML8S5W8ALZk7Ws46qNdJmNQlbFhK8Iok6VKzqbtsNGfZTo+wC96R90UIQSNYZKTn6f0TWUtC+F9L/ByCR9htskYJM+75YGDnMiJjG1NGDbRh+S3Ou2w3Z+xZCTFIkdILBlDH2JBTdN61E1kLctEycrhMjvE10vlOvvBdAzqFClzJUvZb3XaYbM/Y8lIKpJj6nsvFBYu6yOMqBsRuUpEDovIERHZmvP76SKyY/T7N0TknNHx14nIH4jI4yJyUER+yqr0NRCCxmgaujaRtU14bwvJSzU9M4vy2ksVMgOlKeqMgWuKokyWKjlD0XexZ0Z94dFne6Gnpqj06EVkDLgbeDdwFNgnIrtU9YnUaR8FXlTVHxWR64BPApuBfwegqmtF5EeAPxWRi1X1pO2GVCHUrHoTT6mJrG3Ce1uIaTPktqg7Bi4jDBNZin4LpfexZHwVPS9LtHVVT01hQt1cAhxR1acAROR+YCOQNvQbgdtGfz8I/K6ICHABsAdAVf9WRGaA9cBf2hA+DRPKoEuLHerK2ia8t4W2L7dr2qfu/WPSlzJZquRs244m4xJL7SHTWjbQ/bTjMphQN5PAc6l/Hx0dyz1HVY8DLwFnAgeBq0XkNBFZA7wTODv7ABG5QUSmRGTq2LFjtRvRJ8qgKdqE97ZgunNSHlyP4aAjzdC032LJ+MqTo2hf9y4UJ2sK1+mV9zL/YZgCPgX8BbBoCZ6q3qOq61V1/apVq2o/JBY+MCTKeGJfaW5tXm7XYzjoSDM07bdYUivz5PjwZauj+Aj5hAl1M81CL/ys0bG8c46KyGnACuB5nV+NdXNykoj8BfCdVhLnIBY+0BQ2KYqmdcddoIoPLpPV9RiW3X8ppt+ZtrnNuMRCfaXlSNo9O3eCMRFOqJ5awBiDrK5gYuj3AeeOqJdp4DrgQ5lzdgHXA18HrgX2qKqKyHLmV9/+QETeDRzPTOJaQSx8oAlspp3FksKWRtHLXSWr6zEsuv+KifHo+tA16uhNl96tKmTbfUL1lCff17FOUEndjDj3G4HdwLeBB1T1kIjcISJXj077LHCmiBwBfhVIUjB/BPimiHwb+E/Az9tuALTnA7OFxz6x8/FWNbTLYJNC6BIdUSWra0636P4idKYPbaGO3sTCtdtAl94X2zBaMKWqDwMPZ47dmvr7FeADOdc9AzjXiDYpZHnezX2PPnvqd9senk2Kouia6ZlZLt+2J6rd66va7ToNsOj+NxfUIY+V9rOBOjpY1G9AdFsOVqHNuxf6/WmL3qyMbcoHlpVoTWAzx9ZmKFyWOpb+QIH/rQ+zMGm3a0437/5Fdci7SE2Yoq4OZvstRsrQBE3fva62N40lU9SsCDY2B64Dm6FwVVEt3ytjyxArBRCrXC7Rts0x6FMTNG13V9ubRm88+qaouzlwW9ikKNL3arINXNXWck2yZ0xkjSn89SFXbGF/2zY33aqwbR+0vUfTdle1N7bxzcOSqUdfhLya2Vl0oYZ2WU1xyF8ZW1RvvE398wEL0cea7HXr19vog5D92HRfCd/jO9SjL0Hegoqfu2x18IUedWFzZWwMBdL6gj72l019MkXIfixrb1fGd8lTNxDPwo42aFP4KsHO/dPc/uVDvPjyXO4zfBZI8wmXobfvxXw+aIS6FIiNPgi5KLIs86gJZRoCg6HvEdoUvtq5f5otDx5k7kQxleezQJovuM6o8LngyGd2SB3nyEYfrJgYZ2Z2sQOyYmLc+B5tUJR5VITY3oclT90MmMf23YdLjfz4MvFWIC27gM1l4THXobfPrJ5YaQQbfSAFlciKjrtGWVp2jFlbg0efQtOwtwuz7lWoDDVHL5TrLBXfOcsuKIGsPrz/nZNedjOKteaTDZ2ZKaATi443QZ33uKxPY5zTGwz9CE0NTB8WU0B1muncCT21aMzlnIbvzUtsUyt5+vBHj017efljrkvTVmdct63ue1wkz+TKiSjf+95QN23D/aZhb6zhcl1sufI8xsfK42AfnqFJzrJNWsc2tRJrdkjX4bptdceta33dC4/ehlfdNOyNNVxOUKeMMVCadePDMyzz3PLG+aYdB7ht1yFuu/rC1ovObFArMWaHhPIwbVKabetZVV1Xd9xi6+sq9MLQ2wj3m4aGMYfLdT+ASXhdtDjFh7dStAClKGcZYGZ2rhVdZoOKSoxJ0XS2L32IJVW4rfNVZJxdlfJu8h7H0tcm6AV1Y8OLahqKxRzCdXF3oLJnl41nSLosvd1eHmLRB59oQ2HZ3PbRVI6Y32Mb6IVHb8OrbhqKVV0XMiOnC7sD1fHcqiaMQ9FlZal2TXcv6nomVxvdszkhbypHKCrG1zj3wtBvufK8RYt9xsek9te4qXErui50Rk7MtBLU7588WicNk3KzLl6qImMikFv7pQqh9cYG2uiezXmOOnL4pmJ8jnMvqBuAReRoBLXaQmfkxB6O1u2fhNY5Y/ni1ZBV7bJJB2RRZLyaflBD640NtNE9m/0Z8zvgc5x7Yei37z7M3MmFln3upAZ/MUJn5ITk2k3QpH82rZtk/63v4VObL6rVLpcvlW1jElpvbKCN7tnsz5jfAZ/j3AvqxmWHJeH+9Mxs7V3jY6BOYs4MaNM/ddvlUkds87sx6E0bZCmyuzZf5KVufNn9YnwHfI6zkaEXkauATwNjwGdUdVvm99OBzwHvBJ4HNqvqMyIyDnwG+PHRsz6nqndalB9w12F5u8aDOZd2xfmrFuw/mz4+oDyV0jZcv1Q2jYnPfqmLqnkOW7xzXn92fYI6C5/jXEndiMgYcDfwM8AFwAdF5ILMaR8FXlTVHwXuAj45Ov4B4HRVXcv8R+AXReQcS7KfgiseriybwiTsf+TJY7WOLzX4DKtj5mqziJVuMJnncEWRuZxjCQWf42zi0V8CHFHVpwBE5H5gI/BE6pyNwG2jvx8EfldEhPkp0TeIyGnABPAq8Pd2RH8NrlKjmq6Krfq9S1xrGi48Kl9hdddWMsZIN5ikPZrofBM9Mk25jNXrL0sj9iGfiaGfBJ5L/fsocGnROap6XEReAs5k3uhvBP4GWA7crKovtBU6Dy46rCpvu6oWdte51jT6kPIXo/HsEkyMeFXd+KZ6ZPoBiVFHY5DLddbNJcAJ4C3AGuAWEXlb9iQRuUFEpkRk6tixeGiNvHA/japa2LbpgrKCXq5ruFeF5D5ryBchBhn6DJO0x6q68U2pHZNnx5qWGoNcJh79NHB26t9njY7lnXN0RNOsYH5S9kPAn6nqHPC3IrIXWA88lb5YVe8B7oH5zcEbtMMJkq/tTTsO5P5eVQvbJl1Q5hUAzj2GMo8qBo8lBhnyZGo79tl7XHH+Ki+17fNgMnlYVTe+KZ1Z9eyd+6ej3dYvBgrXxNDvA84VkTXMG/TrmDfgaewCrge+DlwL7FFVFZFngQ3A50XkDcBlwKcsye4Fm9ZNnkqvzMJFGmARqrwC1zXcy2go3zXk8xCDDGnY+PDk3SOdxeX7Y2biuFTRlU3pzLJnx76tXwwUbiV1o6rHgRuB3cC3gQdU9ZCI3CEiV49O+yxwpogcAX4V2Do6fjfwRhE5xPwH4w9U9Vu2G+EaMWRslHkFPjyGsj6IwWOJQYY0bITrZVlfTe/ZFpvWTbJ36wae3vZe9m7dsOgDU/WutHmXip4d+7Z+MdgPozx6VX0YeDhz7NbU368wn0qZve4f8o53DbYomLxQ3vS+VV6Ba4+hrA/aRDy2UNQ/K5ePc/m2Pd6pDhsfHtNz23zMPrHzcb74jec4ocqYCB+89Gx+c9PaxvereldcZD/Fvq1fDBlfvVgZ6wNtKZi8MHzLlw6CcKoYW1koXsVR+lh4UdQHMSwMy+uf8THhH145fmojFZ9Uh41wvSrrq8k90/jEzscXjNsJ1VP/bmvsy/rXdvZTF7b1C53x1YtaN11AXng5d1IXVNyE4lC8bHFF6AU2MSwMy+uDN7zutEU1kLq0rV9V1leTe6bxxW88V+t4rIiBGokdg0efgsvFFjZC9jKvIKTHYEpTuF7Mku2DNVsfqiWvTdgI1/PuYTPrJinpYXq8DVyOvUtq1eQesS7SSmMw9CO4Ts8zDcOTc7sEE5oiRPpj6GwHGx9flx/wpEhf3nGb8DH2LqhVExljTOvNQ2+om7aLZVwvajAJw6GbIadJ6Gyrf+uM8xDSl+ODl55d63hTxLBgqApNZexC26AnHr2Nr6rr9LyqxVcw70mllSQmj6AMJqGzjf6tM85JOD07d6J2eemlgmTC1WbWTR6Kxnh6tNguhvFoqp9NrgtB9fTC0NtYLOMjzC9LRRTql0GOCVWhs43+rVPYKlteOvHku9KfvvCbm9ZaN+xZlNGWseh5U/2se10oqqcX1I0Nb9FXmJ/3nKTMZxoxhn9tYKN/Tce5K+H0UkEZbTk7d4KbdhzwXpsoSwFecf6qRvpZV69D6WYvPHob3qKvRQ15z4m1RgfYCzOz7V4xMY4I3LzjANt3H7a6Y1dsq2SXCspK8UI5bekzit25f5otDx5csH5lx77n2Hzx2bUzmurajVC6KeoglaoN1q9fr1NTU7WuyYZDMP9VjWFVnAku37ancMHH3q0bAkg0D1f92vS+2RcU5hdFbb/2HadqnhRRYxC+P/sMkzEt0vM0fIzRujv+x6lFdGmcsXyc/be+x+mzXb7rIvKYqq7P+60X1E3oBUNtEWt2iKsws9V9s37J6N/pHYjyEEN/9hkmY2qSeeYj6soz8mXHbSLUu94L6gbCLzFO0ITqiKEWRh5chZlN77t99+FFK13nTmphBc8EtrJturAwJhRMxjSt50Uf5K6tIamLUO96bwx9DGgzox7LhyoNV5lITe/b5AMh0CgkzqsD/0ePTTvJlmj7AYnhA2Q6pomeF1E9PqKulQW7YK2s2DGuCarmLXyhF9RNLOhbtoerMLPpfct2GTLZgcgUeRtRf+HRZ6Pc9DqWTbPrjmlIuvW2qy9kfNnC1b/jy4Tbrr7Q6nNiGRsYPHqr6Fu2h6sws+l9fVXwzPtgF6UstB3btmtAYtlwpcmYhopifdEnsYwNDIbeKorC12Ui1lcA+grXXb2MTe5r8oKa9ElV39Ux3m1prLbOQUzORd6YxkAr5cHHRyamsRkMvUXkeZwwvzLTZo5wVwopuUDbCp4mfVf0wc4ubLNBY7WdBwlduK0MS1lPIa6xGTh6i0h4x7zqfza5+r7NBfhE0zTAifExPnzZauucctt5kFhTc2HQ05jGZvDoLWPTukluLlgBaCtkiykk9AVbFEDdNEAf1FjZs8raHbJwm8l49ElP6+pfbEX1BkPvAK5DtphCQh+wSQHUTQP0gaJnlbUbFk4++yzcZjoefdHTuvoXY1E9I+pGRK4SkcMickREtub8frqI7Bj9/g0ROWd0/MMiciD130kRuchuE+KD65AtdEjYtvZ/XdikAEL3XR2UtTskLWL6bNt97VvvEtTt6xgpq0qPXkTGgLuBdwNHgX0isktVn0id9lHgRVX9URG5DvgksFlVvwB8YXSftcBOVT1guQ3RwXXoH3IlbYgJNpsUQHZ1ZrIHwC0PHOSmHQeiqlnfpN2mfdKGCjOVy6aehpzYrTsOsWytmYYJdXMJcERVnwIQkfuBjUDa0G8Ebhv9/SDwuyIiurBi2geB+1tL3BG4Dv1D5SCHyA22TQEkcmbDa4grM6Sq3U37pK3RrDMetvQ0ZE56Xf2LcWtNE+pmEkhvC390dCz3HFU9DrwEnJk5ZzPwxbwHiMgNIjIlIlPHjh0zkXtAIISYYHNBt+QZjgShw+wEZe1u0ydtqYUQ9FfIid267fW5taYpvEzGisilwMuq+ld5v6vqPcA9MF+m2Pbzq0KkvLomdetS+0AMi09CTLC5oKqabhGXwMdY2FoglkVboxmCOgw5sVu3vSbn+/5wmRj6aSC9W/BZo2N55xwVkdOAFcDzqd+vo8Cbd42qECnv9/seffbU9bGE8rEsPrni/FUL+id93CVsU1VlG74kvxfB51i0XSCWB1sb9fjUu6ryF65Rt71V5/v+cJlQN/uAc0VkjYi8jnmjvStzzi7g+tHf1wJ7En5eRJYBP0sgfr4qRCoL4fPOD4VYZvIfeTKfWis6HivKaqNXGZBYxqIpupR5lCBkETQX8D0GlR69qh4XkRuB3cAYcK+qHhKRO4ApVd0FfBb4vIgcAV5g/mOQ4F3Ac8lkrk/s3D9duU1f25oivhDL4pNY5GiLvOwb00UtRW2dnpm1XtPIBfKohSvOX8X23Ye5eceBWguCfNKIrqKIUG0Bf/SXEUevqg8DD2eO3Zr6+xXgAwXXfg24rLmIzZCE10VIQqSqED7BCge1qusglsUnschhA7apDyAKms8E6ba3XRAUC73ZBCHb4pP+6m2tmzJKJh0imWxvBpBTvsYIthZ5xBJuxyJHSJTpTJKTH6LmeFP0YUFQU/SpLWXobQmEMiohze1lQ6iilJ+ZBvtJFnkLU997wflu865QtOAoeTG65tE1QdLGmwpqGtmuVuoarhYEdQF9aksZemvoi8LryZUTi16+dAhVtEt7E2qiyFv4wqPPnvqgxL7dYBl/2SbkjSFVtA02rZss3fs01AYTTeBiQVBX0Ke2lKG31E1TisEmNVHkFWSjhlChYhWtVLYVWpuQN6Yt1tqgivaLwSs0oQ5dLAjqCvrUljL0xqPP8xDvvGZtba/RJkViOtEL/o2CySRUmTFvE/LGtMVaGySy3vLAwVMlFNJw5RWaRkOmE40uFgR1BX1qSxlEcxQ0JNavX69TU1O1rinaUd4kz9YlhZAnV3aXogSTKyfYu3WDleeaoIiiSsuxZutDubIK5dRYVTvK7vv0tvdWC+8QTfShjf41kc/0WUVjDERRvK3r9F1sEJHHVHV93m+9oG6a0giuKYS8RR4fvmx1FKGiiUdextG2CXnL7hsSTfXB52KeOrpeFl2Fpsv6Qt91Bb2gbprSCD4ohLwJ1PVvfVNwT8ZkEqps2XmbkDf0cvYitNEHXxPldXS9ijoMSZf1hb7rCnph6JvOnIdKrbJhFNqGvSbGtsqYN21HrLxoF1LtVkyMMzO7ONU3T9eLNqtPw2bb8nQS8sfZdl833eovJv1ziV4Y+qYeYldTq2ys5jM1tq481RCpolWIXR927p/mB68eX3R8fJnk6np2zUMebLUtTye3fOkgCMydWFzr32ZfL+WVvaboBUfflCPtamqVrdV8m9ZNsnfrBp7e9l72bt3QWyU3Rez6sH334VNGM403vv600gyZvVs38KnNFzltW55Ozp3URfImemqzr5fyyl5T9MKjh2YeYqwUQhViphi6HBLHrg9F42uyatt12+ro3vdnZq3Ik+haVeHCtsd9wPV70xtD3xQxUghViJVi6FpIXPRy+ZC1yYvddtxdtq3OmpFE3jby5KWZFj0ni6J5jlCFC328N72gbpYaYqUYuhQSh0zva/rsWMcd8mUbXyaMjy2sBmhL3qp9JMqeU1SgsGnhwrbw8d4sWY/eRqgUiqaoCntDyRVjSFwEW+l9Tfq66bNjppaKZMs7ZkPeMp2qWgxWRHU1KVxoAz7emyVp6G2ESqFpiqKwN6RcsVJKebDxcjXt6zbPjplqLJLNhbxF9MvKifHKldmx6akPKmlJUjdVoZJJIahYaYqQcsVMLWRhY3Vu076OdWVwl9CGfvGppya2xAeVtCQ9+jKPytRL8xFuNaEFQtInTSilsvNdwsbq3KZ9bfLsLmcvtYFpu9vQL7YpsCKZTW2JDyppSRr6stDNlD91Hf41pQVCh6V1KKWyBTWujZqNl71pX5t8ELuUvWQLddodSwZSmcyx2BJYotRNWehm6qW5Dv+a0gKx0id1F9T4QNsFY236uuzZsdKCrlGn3bHouY1S3j7aYuTRi8hVwKeBMeAzqrot8/vpwOeAdwLPA5tV9ZnRb28Hfh/4YeAkcPFoM/FgKPOoihZgZL+uoRagVNECsWZm1F1QA/WpHt90R5u+LpM1JP0WkjKq0+6Qep7uo6Ii74lMMdgSMKhHLyJjwHeAdwNHgX3AB1X1idQ5vwy8XVV/SUSuA/61qm4WkdOAbwI/r6oHReRMYEZVCxNgm9SjtwmftcXLYFIvvksoq42eRZIelx2H8WWygOqB18YGiGLcTFClY6HGPrTud0HnTRZqQbEOu+zPtvXoLwGOqOpTqvoqcD+wMXPORuAPR38/CPy0iAjwHuBbqnoQQFWfLzPyMaBp3RzbiCU0tYW6C2rqUj1dojuqZA019qH7sAs6X7VQCxaW8o7BloAZdTMJPJf691Hg0qJzVPW4iLwEnAn8GKAishtYBdyvqr+VfYCI3ADcALB69eq6bbCOGHKV88K5K85fxfbdh7l5x4FoKBlT1F1Qc/OOA8b3LqM0XNEdbSiOKooiW3VyTGSBwXU15qEXvMVKO6ZR1hfJzms2Snnbhuusm9OAfwFcDLwMfHUUXnw1fZKq3gPcA/PUjWOZgsPUSKSVJMZMjLrGrs6Cmia1U3xlG7UdCxPuNrmPzzEPnbEF8RjGIrTZQjMkTKibaeDs1L/PGh3LPWfEy69gflL2KPDnqvp3qvoy8DDw422F7jKa1jlxEVabLOYou9ZlrZi8MD4PSZjsM+xvOxamstoY8zpj3AXqJDS62kcmHv0+4FwRWcO8Qb8O+FDmnF3A9cDXgWuBPaqaUDb/UUSWA68CPwncZUv4KsS46KRJnZOd+6drl2KtQluv1PVWcFn6oghZztPHeLelOEwpirbPqTvGXaBOQqOrfVRp6Eec+43AbubTK+9V1UMicgcwpaq7gM8CnxeRI8ALzH8MUNUXReR3mP9YKPCwqj7kqC0LECPVAfVf3qQdRWgaVrc11D743CSML8vGSMvqK+y3QXGYyNr2OU3GuEkfxuhQ2UBZGeuutc+Io1fVh5mnXdLHbk39/QrwgYJr7wPuayFjIxQp+S0PHATCGfu6L2/ZLH+bkLGtofbJ58a2mbgveeo8J88o+SrTEaND1RZ9a1dvV8YWKfMJVW91x2ExR3rF+atqcXxlL2WbVK22hbV8cpUxpan5lMf0OUXzJSuX51c/tPkxDp2S6Qq+2tVmnqwOelvrpixrwyaXDPWKGv3RY9O8/52TPPLkMb4/M8vK5eOows07DpzaS9OkDkaWtqiLtl6pb64ytnC5TB6bVIZJu4uM0umnLWNifMxp8bTQKZmu0LdoqLcefVXWhq0BK8s+KXoBH3nyGHu3buCuzRfxytxJZmbnCjNXXHnONrzStrVi+ogQO1cV6fJLs3OlY2xD1r6WXPbRLp/RUG89+kSZb3ngICdyyjzYGrA2RY1MJstces6xecl9gOtspDyUzZeUjbENWWObP7EFH+3yGQ311tBD/qITsDtgZYNVNWFpOtAuDHIsmRIx1ai3gRBURlOj1ETWvPG685q1Xlfx+tBdH9Skz4SGXht6cD9gZYNV9QKGWokYS0ZBbDXqbSDEmDbV8bqyFunNndesXaTrrsbNp+66jnh9RkOV1St9I3T1yrqoqvhX5n2EqhYYS5XAuhUtY15iniB0Bcg6qCtrmd5AfgkK2+MWi+7ags3opKx6Ze89eteo8qbKvIJQq+xiyZRoUqM+dnRp5WRdWZvoje1xi0V3bcHXPNlg6C2gzWCFmBCNoXhVmRxF53YFXZrkriNrld740KlYdLdr6G165YBixFKYqW6N+gFhUaY3vnQqFt3tGgaPfgkiFnqhSI4YZBuwGCZ604dsmD5imIwdMGDAgB6g7VaCAwYMGDCgwxgM/YABAwb0HIOhHzBgwICeYzD0AwYMGNBzDIZ+wIABA3qO6LJuROQY8L2Gl78Z+DuL4thErLINctVHrLINctVDrHJBM9neqqqr8n6IztC3gYhMFaUXhUassg1y1Uessg1y1UOscoF92QbqZsCAAQN6jsHQDxgwYEDP0TdDf09oAUoQq2yDXPURq2yDXPUQq1xgWbZecfQDBgwYMGAx+ubRDxgwYMCADAZDP2DAgAE9R2cMvYhcJSKHReSIiGzN+f10Edkx+v0bInJO6rePjY4fFpErY5BLRM4RkVkROTD67/dsymUo27tE5JsiclxErs38dr2I/PXov+sjkutEqs92eZbrV0XkCRH5loh8VUTemvotZH+VyeWsvwxl+yUReXz0/P8jIhekfgv5XubK5fq9rJIrdd77RURFZH3qWPP+UtXo/wPGgO8CbwNeBxwELsic88vA743+vg7YMfr7gtH5pwNrRvcZi0Cuc4C/Ctxn5wBvBz4HXJs6/ibgqdH/zxj9fUZouUa//UPA/roCWD76+9+nxjJ0f+XK5bK/asj2w6m/rwb+bPR36PeySC5n76WJXKPzfgj4c+BRYL2N/uqKR38JcERVn1LVV4H7gY2ZczYCfzj6+0Hgp0VERsfvV9V/VNWngSOj+4WWyzUqZVPVZ1T1W8DJzLVXAl9R1RdU9UXgK8BVEcjlEiZyPaKqL4/++Shw1ujv0P1VJJdrmMj296l/vgFIsj+CvpclcrmEib0A+A3gk8ArqWOt+qsrhn4SeC7176OjY7nnqOpx4CXgTMNrQ8gFsEZE9ovI/xKRn7AkUx3ZXFzr+t6vF5EpEXlURDZZkqmJXB8F/rThtb7kAnf9ZSybiPyKiHwX+C3gP9S5NoBc4O69rJRLRH4cOFtVH6p7bRmGrQTD4W+A1ar6vIi8E9gpIhdmPI0Bi/FWVZ0WkbcBe0TkcVX9rk8BROTngPXAT/p8bhUK5AreX6p6N3C3iHwI+ARgdQ6jKQrkCvZeisgy4HeAj9i+d1c8+mng7NS/zxodyz1HRE4DVgDPG17rXa5RCPY8gKo+xjzn9mOW5DKVzcW1Tu+tqtOj/z8FfA1Y51MuEfmXwMeBq1X1H+tcG0Aul/1lLFsK9wObGl7rRS7H72WVXD8E/DPgayLyDHAZsGs0Iduuv1xMOjiYxDiN+QmuNbw2iXFh5pxfYeGk5wOjvy9k4STGU9ib9Gkj16pEDuYnZ6aBN/nss9S5/5XFk7FPMz+xeMbobyuytZTrDOD00d9vBv6anMksh2O5jvkX/9zM8aD9VSKXs/6qIdu5qb/fB0yN/g79XhbJ5ey9rKP7o/O/xmuTsa36y8qA+/gP+FfAd0YK/fHRsTuY92AAXg98iflJir8E3pa69uOj6w4DPxODXMD7gUPAAeCbwPsC9NnFzHN9P2A++jmUuvbfjmQ+AvxCDHIB/xx4fKTwjwMf9SzX/wT+72jMDgC7IumvXLlc95ehbJ9O6fkjpAxb4PcyVy7X72WVXJlzv8bI0Lftr6EEwoABAwb0HF3h6AcMGDBgQEMMhn7AgAEDeo7B0A8YMGBAzzEY+gEDBgzoOQZDP2DAgAE9x2DoBwwYMKDnGAz9gAEDBvQc/x9FpJiG0BqxKwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(param_values[:,0], param_values[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = np.zeros([param_values.shape[0]])\n",
    "\n",
    "for i, (beta, gamma) in enumerate(param_values):\n",
    "    Y[i] = tf.reduce_sum(model_I(deriv, y0, dt, days, N, beta, gamma))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "Si = sobol.analyze(problem, Y, calc_second_order=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'S1': array([0.7475557 , 0.16852292]),\n",
       " 'S1_conf': array([0.19277332, 0.14673407]),\n",
       " 'ST': array([0.89294114, 0.26704384]),\n",
       " 'ST_conf': array([0.22591018, 0.08714372])}"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Si"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}