Lab-part-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YWa2nHPPhVXP"
   },
   "source": [
    "# SIR epidemic model - [data assimilation]\n",
    "\n",
    "This notebook is the second part of inspecting SIR epidemic model. \n",
    "Here we will perform epidemic data assimilation with proposed numerical model using manually implemented 4D-Var method in TensorFlow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jxcqgv3UhVXb"
   },
   "source": [
    "## 4D-Var method\n",
    "\n",
    "Data assimilation is the process of optimizing the fit between gathered data and some prior (theoretical or based on previous measure).\n",
    "\n",
    "`4D-Var` is one of the methods for data assimilation - it widely used for weather forecast as a successor of `3D-Var` method. \n",
    "\n",
    "<img src=\"https://i.ibb.co/sWLCYSK/Andersson-NL115-figure-1-box-754px.png\" width=\"400\">\n",
    "\n",
    "`4D-Var` loss for $x$ model parameters is defined:\n",
    "\n",
    "$$J(x) = (x-x_b)^T \\boldsymbol{B}^{-1} (x-x_b) + \\sum_{i=0}^{n}(y_i-\\mathcal{H}_i(x_i))^T \\boldsymbol{R}_i^{-1} (y_i-\\mathcal{H}_i(x_i))$$\n",
    "\n",
    "where: \n",
    "- $x_b$ - our prior parameters\n",
    "- $\\boldsymbol{B}$ - a background error covariance\n",
    "- $\\boldsymbol{R}$ - a observational error covariance\n",
    "- $y$ - observation data\n",
    "- $\\mathcal{H}$ - operator that transforms input space vector into output space\n",
    "\n",
    "`3D-Var` method assumed that all observed data happened at one fixed point of time for given assimilation window. So when performing DA for weather conditions four times a day a six hours  time span is considered as one point in time. As we know weather conditions can change through time so `4D-Var` introduced time dimension to count in temporal changes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "8g1_KGGChVXh"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "Z5KGi_B4hVYO"
   },
   "outputs": [],
   "source": [
    "assert tf.__version__[0] == '2', 'Tensorflow 2.x version required'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "JUjsaU8xhVYq"
   },
   "outputs": [],
   "source": [
    "np.random.seed(2)\n",
    "tf.random.set_seed(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yjdyKip-hVZN"
   },
   "source": [
    "### Task 1: Complete equations according to formulas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "WMz1X2TnhVZR"
   },
   "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": {
    "id": "WI7wC9V_hVZv"
   },
   "source": [
    "### Task 2: Complete the formula for each component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "uTaglbljhVZ0"
   },
   "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": "code",
   "execution_count": 6,
   "metadata": {
    "id": "BuyvGHZuhVaV"
   },
   "outputs": [],
   "source": [
    "# Total population\n",
    "N = tf.constant([2500.])\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",
    "# A grid of time points (in days)\n",
    "days = 160\n",
    "t = np.linspace(0, days, days)\n",
    "dt = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "WSFqca1_hVbQ"
   },
   "outputs": [],
   "source": [
    "prior = tf.constant([[0.3],[0.1]])\n",
    "trainable_weights = tf.Variable([[0.2],[0.1]]) # posterior\n",
    "B_inv = tf.eye(2) * 100000.\n",
    "R_inv = tf.eye(1) * 0.00005"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "i0lRyt0qhVbj"
   },
   "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": 9,
   "metadata": {
    "id": "U0Hz7awYhVb5"
   },
   "outputs": [],
   "source": [
    "def model_x(x_pred, horizon=70):\n",
    "    return model_I(deriv, y0, dt, days, N, x_pred[0,0], x_pred[1,0])[:horizon]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 206
    },
    "id": "Mg0o1WlCjR6v",
    "outputId": "8528b4c0-4a09-48cf-f65e-394dc0d5adbb"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Date</th>\n",
       "      <th>Active cases</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2020-02-28</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2020-02-29</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2020-03-01</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2020-03-02</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2020-03-03</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Date  Active cases\n",
       "0  2020-02-28             1\n",
       "1  2020-02-29             1\n",
       "2  2020-03-01             1\n",
       "3  2020-03-02             1\n",
       "4  2020-03-03             2"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('New Zealand COVID active cases.csv', delimiter=';')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "VPO4pbYWhVcM"
   },
   "outputs": [],
   "source": [
    "raw = model_x(trainable_weights, horizon=days)\n",
    "y = data.loc[:, 'Active cases'].values\n",
    "y = y[:, None, None]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "br8UGDzwhVdX"
   },
   "source": [
    "### Task 3: Complete the loss function according to the formula\n",
    "\n",
    "Hint: You should use `tf.transpose` and `tf.reduce_sum`. In `tf.transpose` for observation part notice that first dimension is batch size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PAxRHqUmhVdc"
   },
   "source": [
    "`4D-Var` loss for $x$ model parameters is:\n",
    "\n",
    "$$J(x) = (x-x_b)^T \\boldsymbol{B}^{-1} (x-x_b) + \\sum_{i=0}^{n}(y_i-\\mathcal{H}_i(x_i))^T \\boldsymbol{R}_i^{-1} (y_i-\\mathcal{H}_i(x_i))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YxYBwW9AhVdg"
   },
   "source": [
    "With exemplary substitution it would be:\n",
    "\n",
    "\\begin{gather}\n",
    "J(x) = \n",
    "\\begin{bmatrix} x_{0} & x_{1} \\end{bmatrix} \n",
    "\\begin{bmatrix} 1 & 0 \\\\ 0 & 1 \\end{bmatrix}\n",
    "\\begin{bmatrix} x_{0} \\\\ x_{1} \\end{bmatrix}\n",
    "      +\n",
    "\\sum_{i=0}^{n}\n",
    "\\begin{bmatrix} y_{0}^i & y_{1}^i \\end{bmatrix} \n",
    "\\begin{bmatrix} 1 & 0 \\\\ 0 & 1 \\end{bmatrix}\n",
    "\\begin{bmatrix} y_{0}^i \\\\ y_{1}^i \\end{bmatrix}\n",
    "\\end{gather}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "id": "pJqhn7V_hVdj"
   },
   "outputs": [],
   "source": [
    "# (batch_size, rows, cols)\n",
    "@tf.function\n",
    "def custom_loss(prior, x_pred):\n",
    "    custom_loss = tf.transpose(x_pred-prior) @ B_inv @ (x_pred-prior) + \\\n",
    "      tf.reduce_sum(tf.transpose(y-model_x(x_pred),perm=[0,2,1]) @ R_inv @ (y-model_x(x_pred)))\n",
    "    return custom_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = tf.optimizers.Adam(learning_rate=0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loss for 0th iteration:  1591.4475\n",
      "Loss for 10th iteration:  1364.8013\n",
      "Loss for 20th iteration:  1175.9275\n",
      "Loss for 30th iteration:  1001.5325\n",
      "Loss for 40th iteration:  839.0573\n",
      "Loss for 50th iteration:  684.324\n",
      "Loss for 60th iteration:  541.3607\n",
      "Loss for 70th iteration:  415.45016\n",
      "Loss for 80th iteration:  311.19174\n",
      "Loss for 90th iteration:  229.98332\n",
      "Loss for 100th iteration:  170.2839\n",
      "Loss for 110th iteration:  128.70306\n",
      "Loss for 120th iteration:  101.1546\n",
      "Loss for 130th iteration:  83.74431\n",
      "Loss for 140th iteration:  73.22603\n",
      "Loss for 150th iteration:  67.14407\n",
      "Loss for 160th iteration:  63.777145\n",
      "Loss for 170th iteration:  61.99328\n",
      "Loss for 180th iteration:  61.089775\n",
      "Loss for 190th iteration:  60.652817\n"
     ]
    }
   ],
   "source": [
    "num_iters = 200\n",
    "penalty = np.zeros(num_iters, np.float64)\n",
    "\n",
    "for i in range(num_iters):\n",
    "    with tf.GradientTape() as tape:\n",
    "        loss = custom_loss(prior, trainable_weights)\n",
    "    grads = tape.gradient(loss, trainable_weights)\n",
    "    opt.apply_gradients([(grads, trainable_weights)])\n",
    "    penalty[i] = loss\n",
    "    if i % 10 == 0:\n",
    "        print(f'Loss for {i}th iteration: ', loss.numpy().squeeze())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x15f6e6a00>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAmRklEQVR4nO3deXgV9fn+8feTFZJAAiQsDfsOigKeIu5V3KvFpbUubWmrpdal+tXWutTqr6221lqrlVqxWnFF3GlLXapUbBUkQVYBCYuSsEUiYQ1keX5/nIlGJCxZzpzk3K/rynVmPjOTeTI5uc9kls+YuyMiIokhKewCREQkdhT6IiIJRKEvIpJAFPoiIglEoS8ikkBSwi5gb3Jzc713795hlyEi0qIUFhZ+7O55e5oW16Hfu3dvCgoKwi5DRKRFMbMP65umwzsiIglEoS8ikkAU+iIiCUShLyKSQPYZ+mb2sJltMLOFu7VfaWZLzGyRmf2uTvsNZlZkZkvN7JQ67acGbUVmdn3T/hgiIrI/9ufqnUeA+4BHaxvM7HhgLHCou+80s85B+1DgfOAg4EvAv81sYLDYBOAkoBiYbWZT3f39pvpBRERk3/YZ+u4+w8x679b8I+C37r4zmGdD0D4WmBy0rzSzImBUMK3I3VcAmNnkYF6FvohIDDX0mP5A4Bgzm2Vmb5rZl4P2fGB1nfmKg7b62r/AzMabWYGZFZSWljaouPIdldz16lKWl25t0PIiIq1VQ0M/BegIjAZ+CkwxM2uKgtx9ortH3D2Sl7fHG8r2qbK6hokzVvCX/yxvipJERFqNhoZ+MfC8R70L1AC5QAnQo8583YO2+tqbRW5WOheM6skL75VQsmlHc61GRKTFaWjovwgcDxCcqE0DPgamAuebWbqZ9QEGAO8Cs4EBZtbHzNKInuyd2sja9+oHx/YF4MEZK5pzNSIiLcr+XLL5FPAOMMjMis3sYuBhoG9wGedkYFyw178ImEL0BO3LwOXuXu3uVcAVwCvAYmBKMG+zyc9py9jh+Uye/REbt+5szlWJiLQYFs/PyI1EIt6YDteKNmzhpLtncMXx/bn25EFNWJmISPwys0J3j+xpWqu+I7d/53acMrQrk95exZaKyrDLEREJXasOfYDLju/H5ooqnpj1UdiliIiErtWH/iHdczi6fy4P/XclFZXVYZcjIhKqVh/6EN3bL92yk2cLi8MuRUQkVAkR+kf07cTwHjk8MGM5ldU1YZcjIhKahAh9M+OK4/uzumwHL8xptnvCRETiXkKEPsCYIZ05pHs2976xjF1V2tsXkcSUMKFvZvzfSQMp/mSHju2LSMJKmNAH+MrAPEb0zOG+N5axs0pX8ohI4kmo0DczrjlpIGvKK5gye/W+FxARaWUSKvQBju6fy5d7d2DC9OW6bl9EEk7Chb6Zce3Jg1i3uYJJb68KuxwRkZhKuNAHGN23E8cPymPC9CLKt6tPHhFJHAkZ+gDXnTqYLTur+PObRWGXIiISMwkb+kO6tefsEfn87X+rWKOna4lIgkjY0Ae45qSB4HD3ax+EXYqISEwkdOh375DBd47oxXNzinl/zeawyxERaXb787jEh81sQ/BoxN2nXWtmbma5wbiZ2b1mVmRm881sZJ15x5nZsuBrXNP+GA135QkDyG6byv/7+yLi+SliIiJNYX/29B8BTt290cx6ACcDdZ9OchrRh6EPAMYD9wfzdgRuAQ4HRgG3mFmHxhTeVLIzUrn25EHMWlnGtAXrwi5HRKRZ7TP03X0GULaHSXcD1wF1d4/HAo8GD0mfCeSYWTfgFOA1dy9z90+A19jDB0lYLhjVkyHd2nP7tMXs2KUbtkSk9WrQMX0zGwuUuPu83SblA3X7NygO2upr39P3Hm9mBWZWUFpa2pDyDlhyknHrmUMp2bSDB2Ysj8k6RUTCcMChb2YZwI3AL5q+HHD3ie4ecfdIXl5ec6xijw7v24mvHtKNv7y5nBJdwikirVRD9vT7AX2AeWa2CugOzDGzrkAJ0KPOvN2Dtvra48qNpw8B4PZpi0OuRESkeRxw6Lv7Anfv7O693b030UM1I919HTAV+E5wFc9ooNzd1wKvACebWYfgBO7JQVtcyc9py6XH9eOf89cyc8XGsMsREWly+3PJ5lPAO8AgMys2s4v3Mvs0YAVQBDwIXAbg7mXAr4DZwdcvg7a4c+lx/cjPacstLy3S83RFpNWxeL42PRKJeEFBQczX++qidYx/rJAbThvMD4/rF/P1i4g0hpkVuntkT9MS+o7c+px8UFdOGtqFP/57GcWfbA+7HBGRJqPQr8etXzso+jpVd+qKSOuh0K9Hfk5b/u+kAfx78QZefX992OWIiDQJhf5efO+oPgzu2o5bpy5i686qsMsREWk0hf5epCYncdvZw1i3uYI/qvtlEWkFFPr7cFivDlwwqid/e3sVi9aUh12OiEijKPT3w89OGUxO21RuemEh1TU6qSsiLZdCfz9kZ6Ty8zOGMHf1Jp5896N9LyAiEqcU+vvprOH5HNmvE797eQkbNleEXY6ISIMo9PeTmXHb2cPYWVXDrX9fFHY5IiINotA/AH1yM7lqzACmLVjHa7p2X0RaIIX+ARp/bF8Gd23HzS8uZEtFZdjliIgcEIX+AUpNTuI35wxj/ZYKfv/K0rDLERE5IAr9BhjRswPjjujNozM/ZM5Hn4RdjojIflPoN9BPThlE1/ZtuOG5BeyqUr/7ItIyKPQbKCs9hV+NPZil67fw4Fsrwi5HRGS/KPQb4cShXfjqsG7c8/oyVpRuDbscEZF92p/HJT5sZhvMbGGdtjvNbImZzTezF8wsp860G8ysyMyWmtkpddpPDdqKzOz6Jv9JQnLLmUNJT0nixhcWqN99EYl7+7On/whw6m5trwEHu/shwAfADQBmNhQ4HzgoWObPZpZsZsnABOA0YChwQTBvi9e5fRtuPH0IM1eU8UxBcdjliIjs1T5D391nAGW7tb3q7rUdzM8EugfDY4HJ7r7T3VcSfUD6qOCryN1XuPsuYHIwb6vwzUgPRvXuyG3TFlO6ZWfY5YiI1Kspjul/H/hXMJwPrK4zrThoq6/9C8xsvJkVmFlBaWlpE5TX/JKSjNvPGcaOXdX88h/vh12OiEi9GhX6ZnYTUAU80TTlgLtPdPeIu0fy8vKa6ts2u/6ds7j8+P78fd4api/ZEHY5IiJ71ODQN7PvAmcAF/lnZzBLgB51ZusetNXX3qpc+pW+9O+cxc9fXMg2PV5RROJQg0LfzE4FrgO+5u7b60yaCpxvZulm1gcYALwLzAYGmFkfM0sjerJ3auNKjz/pKcn89pxhlGzawV2v6vGKIhJ/9ueSzaeAd4BBZlZsZhcD9wHtgNfMbK6Z/QXA3RcBU4D3gZeBy929OjjpewXwCrAYmBLM2+pEenfkosN78sjbK5m3elPY5YiIfI7F87XlkUjECwoKwi7jgG2uqOSkP7xJTts0pl55FOkpyWGXJCIJxMwK3T2yp2m6I7cZtG+Tym/OGcbS9Vu4742isMsREfmUQr+ZnDC4C+eMzOfP/1nOwpLysMsREQEU+s3qljMOolNmGj95Zp564hSRuKDQb0bZGancfvYwlqzbwn3TdZhHRMKn0G9mJw7twjkj8vnz9CId5hGR0Cn0Y+AXZw6lgw7ziEgcUOjHQE5GGr8JDvNM0GEeEQmRQj9GThzahbNH5DNhehGL1ugwj4iEQ6EfQ7cEh3munTKPnVXVYZcjIglIoR9DORlp3HFu9DDPH9Q3j4iEQKEfYycM7sKFh/dk4lsrmLliY9jliEiCUeiH4KbTh9CrYwbXTpnH5orKsMsRkQSi0A9BZnoKf/jmcNaW7+D/TdWTtkQkdhT6IRnZswNXHN+f5+YU868Fa8MuR0QShEI/RFeOGcCw/GxufGEBGzZXhF2OiCQAhX6IUpOTuPubw9m+q5rrnptPPD/bQERah/15ctbDZrbBzBbWaetoZq+Z2bLgtUPQbmZ2r5kVmdl8MxtZZ5lxwfzLzGxc8/w4LU//zlncePoQ/rO0lMdnfRR2OSLSyu3Pnv4jwKm7tV0PvO7uA4DXg3GA04g+F3cAMB64H6IfEsAtwOHAKOCW2g8KgW+P7sUxA3K57Z/vs2z9lrDLEZFWbJ+h7+4zgLLdmscCk4LhScBZddof9aiZQI6ZdQNOAV5z9zJ3/wR4jS9+kCSspCTjrm8cSkZaCj+ePJeKSt2tKyLNo6HH9Lu4e+0lJ+uALsFwPrC6znzFQVt97RLo3L4Nv//GISxeu5nfvbw07HJEpJVq9Ilcj559bLIzkGY23swKzKygtLS0qb5ti3DC4C5898jePPy/lUxfuiHsckSkFWpo6K8PDtsQvNYmVAnQo8583YO2+tq/wN0nunvE3SN5eXkNLK/luv60wQzu2o6fPjOP0i07wy5HRFqZhob+VKD2CpxxwEt12r8TXMUzGigPDgO9ApxsZh2CE7gnB22ymzapyfzpghFsqaji2mfmUVOjyzhFpOnszyWbTwHvAIPMrNjMLgZ+C5xkZsuAE4NxgGnACqAIeBC4DMDdy4BfAbODr18GbbIHA7q04+YzhjLjg1Ie/t/KsMsRkVbE4vmGoEgk4gUFBWGXEQp354ePFTJ96QZeuOwoDs7PDrskEWkhzKzQ3SN7mqY7cuOUmXHHuYfQMTONH09+j+27qsIuSURaAYV+HOuQmcbd5w1n5cfb+NU/1BuniDSeQj/OHdk/l0uP68dT765mmnrjFJFGUui3ANecNJBDu2dz/XPzWbNpR9jliEgLptBvAVKTk7j3ghFU1zhXT55LVXVN2CWJSAul0G8henXK5NdnH8y7q8q4942isMsRkRZKod+CnD2iO+eO7M6f3ljG28s/DrscEWmBFPotzC/HHkSf3EyunjyXj7eqmwYROTAK/RYmMz2FCReOZNOOSq6dom4aROTAKPRboCHd2vOLM4by5gelPPjWirDLEZEWRKHfQl10eE9OH9aVO19ZypyPPgm7HBFpIRT6LZSZ8ZtzDqFrdhuufPI9yrdXhl2SiLQACv0WLLttKvddOJL1myv42XPziefO80QkPij0W7jhPXK47tRBvLxoHY/P/DDsckQkzin0W4FLju7LVwbl8at/LmbRmvKwyxGROKbQbwWSkoy7vnEoHTJSufLJ99i2U90wi8ieKfRbiU5Z6dxz/ghWbdzGzS8tDLscEYlTjQp9M/s/M1tkZgvN7Ckza2NmfcxslpkVmdnTZpYWzJsejBcF03s3yU8gnxrdtxM/HjOA5+eU8GxhcdjliEgcanDom1k+8GMg4u4HA8nA+cAdwN3u3h/4BLg4WORi4JOg/e5gPmliV54wgNF9O3Lziwsp2rA17HJEJM409vBOCtDWzFKADGAtcALwbDB9EnBWMDw2GCeYPsbMrJHrl90kJxn3nD+CtmnJXP7EHHbsqg67JBGJIw0OfXcvAX4PfEQ07MuBQmCTu9eeSSwG8oPhfGB1sGxVMH+n3b+vmY03swIzKygtLW1oeQmtS/s23P3N4XywYYuO74vI5zTm8E4HonvvfYAvAZnAqY0tyN0nunvE3SN5eXmN/XYJ67iBeVx5wgCeLSxmyuzVYZcjInGiMYd3TgRWunupu1cCzwNHATnB4R6A7kBJMFwC9AAIpmcDGxuxftmHq8YM4Oj+udz80kLeX7M57HJEJA40JvQ/AkabWUZwbH4M8D4wHfh6MM844KVgeGowTjD9DVe/Ac0qOcn44/nDyclI5bInCtlcof55RBJdY47pzyJ6QnYOsCD4XhOBnwHXmFkR0WP2DwWLPAR0CtqvAa5vRN2yn3Kz0rnvwpGs/mQHP3tW/fOIJDqL5xCIRCJeUFAQdhmtwsQZy7l92hJ+ccZQvn90n7DLEZFmZGaF7h7Z0zTdkZsgfnBMX04a2oXbpy1W//siCUyhnyDMjN9/41C65bThiifmULZtV9gliUgIFPoJJLttKvdfdBgfb9vF1U/P1fN1RRKQQj/BHJyfzS1nDmXGB6Xc+8aysMsRkRhT6CegC0f15NyR3fnjv5fx+uL1YZcjIjGk0E9AZsZtZx/Mwfntufrpuaz6eFvYJYlIjCj0E1Sb1GT+8q3DSEkyfvhYoR68IpIgFPoJrHuHDO69YATLNmzRg9VFEoRCP8EdMyCPn54ymH/MX8tD/10Zdjki0swU+sKlx/XltIO78pt/LeHt5R+HXY6INCOFvmBm3PmNQ+mTm8kVT77Hmk07wi5JRJqJQl8AyEpP4YFvH8auqhp+9HghFZV64pZIa6TQl0/1y8virvMOZV5xObe8tEgndkVaIYW+fM4pB3XliuP783TBaia9vSrsckSkiSn05QuuOWkgJw7pwq/+uZj/LtOJXZHWRKEvX5AUPHGrf14Wlz85h5W6Y1ek1VDoyx5lpafw13ERkgwumTRbj1oUaSUaFfpmlmNmz5rZEjNbbGZHmFlHM3vNzJYFrx2Cec3M7jWzIjObb2Yjm+ZHkObSo2MGf77oMD7cuJ0fP/Ue1eqKWaTFa+ye/j3Ay+4+GDgUWEz02bevu/sA4HU+exbuacCA4Gs8cH8j1y0xcES/Ttz6tYP4z9JS7nh5SdjliEgjNTj0zSwbOJbgwefuvsvdNwFjgUnBbJOAs4LhscCjHjUTyDGzbg1dv8TOt0b34tujezFxxgqeLSwOuxwRaYTG7On3AUqBv5nZe2b2VzPLBLq4+9pgnnVAl2A4H1hdZ/nioO1zzGy8mRWYWUFpaWkjypOm9Iszh3JE307c+PwCCj/UM3ZFWqrGhH4KMBK4391HANv47FAOAB69u+eADgS7+0R3j7h7JC8vrxHlSVNKTU7izxeNpGt2G374WKG6ahBpoRoT+sVAsbvPCsafJfohsL72sE3wuiGYXgL0qLN896BNWogOmWn8dVyEispqLplUwFb1wS/S4jQ49N19HbDazAYFTWOA94GpwLigbRzwUjA8FfhOcBXPaKC8zmEgaSEGdmnHny4cwdL1W7jyyTlUVdeEXZKIHIDGXr1zJfCEmc0HhgO3A78FTjKzZcCJwTjANGAFUAQ8CFzWyHVLSI4f1Jlfjj2I6UtLufXv6qNHpCVJaczC7j4XiOxh0pg9zOvA5Y1Zn8SPiw7vxUcbt/PAjBX06pjJD47tG3ZJIrIfGhX6kth+dupgVn+yndv/tZjuHdpy2jBdgSsS79QNgzRYUpLxh/OGM7xHDlc/PZf3PtKlnCLxTqEvjdImNZkHvxOhc/t0LplUwOqy7WGXJCJ7odCXRsvNSudv3x1FVY0z7uF3Kdu2K+ySRKQeCn1pEv07Z/HXcRFKNu3ge397l226hl8kLin0pcl8uXdHJlw4koVrNnPp44XsqtI1/CLxRqEvTerEoV34zTnDeGvZx/zkmXnUqDtmkbiiSzalyZ0X6UHZtl389l9L6JiZxi1nDsXMwi5LRFDoSzP54bF9+XjLTv7635XktUvn8uP7h12SiKDQl2ZiZtx4+hA2btvFna8spUNGGhce3jPsskQSnkJfmk1SkvG7rx/Cpu27uOnFBWSkJXPWiC88QkFEYkgncqVZpSYncf+3DmN0n05cM2Uu/5yvjlVFwqTQl2bXJjWZh74b4bBeHbhq8nu8umhd2CWJJCyFvsRERloKD3/3yxycn83lT85h+tIN+15IRJqcQl9ipl2bVCZ9fxSDurbjh48V8r+ij8MuSSThKPQlprLbpvLY9w+nb24mF0+azawVG8MuSSShKPQl5jpkpvH4JYfTvUMG33tkNu8sV/CLxEqjQ9/Mks3sPTP7RzDex8xmmVmRmT1tZmlBe3owXhRM793YdUvLlZuVzpM/OJz8nLZ892/vMuOD0rBLEkkITbGnfxWwuM74HcDd7t4f+AS4OGi/GPgkaL87mE8SWOd2bZg8fjT98rK4ZFIB/35/fdglibR6jQp9M+sOfBX4azBuwAnAs8Esk4CzguGxwTjB9DGmDlkSXqesdJ76wWiGdGvHpY8XMm2BruMXaU6N3dP/I3AdUNuHbidgk7vXdqZeDNTegpkPrAYIppcH83+OmY03swIzKygt1b/8iSA7I5XHLzmc4T1yuOLJOTxXWBx2SSKtVoND38zOADa4e2ET1oO7T3T3iLtH8vLymvJbSxyrvZzziH6duPaZeTw4Y0XYJYm0So3Z0z8K+JqZrQImEz2scw+QY2a1ffp0B0qC4RKgB0AwPRvQZRvyqcz06A1cXx3WjdumLeY30xbjrv74RZpSg0Pf3W9w9+7u3hs4H3jD3S8CpgNfD2YbB7wUDE8Nxgmmv+H6i5bdpKckc+8FI/j26F48MGMFP3lmPpXVegKXSFNpjl42fwZMNrNfA+8BDwXtDwGPmVkRUEb0g0LkC5KTjF+OPYjcrHTu/vcHlG3byX0XjiQzXZ3CijSWxfPOdiQS8YKCgrDLkBA9MetDbn5xIUO6teehcV+ma3absEsSiXtmVujukT1N0x25EtcuOrwXD437Mqs+3sbYCf9lYUl52CWJtGgKfYl7xw/uzLM/OpJkM8574B3dxCXSCAp9aRGGdGvPi5cfxYDOWfzgsQIenLFCV/aINIBCX1qMzu3bMHn8EZx2cFdum7aYqybPZfuuqn0vKCKfUuhLi9I2LZkJF47kulMH8ff5azjnz2/z4cZtYZcl0mIo9KXFMTMu+0p/Jn1vFGvLKzjzT//Vk7hE9pNCX1qsYwfm8Y8rj6Z7hwy+/8hs7np1KVW6kUtkrxT60qL16JjBcz86km8c1p0/vVHE+RNnUrJpR9hlicQthb60eG3Tkvnd1w/lnvOHs2TdFk6/5y1eXrgu7LJE4pJCX1qNscPz+ceVR9OzYwaXPl7Iz19coKt7RHaj0JdWpXduJs/96Eh+cEwfHp/5Eafd8xazV5WFXZZI3FDoS6uTlpLETV8dyuTxo6lx57wH3uHX/3ifisrqsEsTCZ1CX1qt0X078fJVx3LhqJ789b8rOf3etyj8UHv9ktgU+tKqZaancNvZw3js4lFU7Krm3Pvf4Ybn57Np+66wSxMJhUJfEsIxA/J47ZrjuOToPkwpKGbMXW/y/Jxi9d8jCUehLwkjMz2Fn58xlKlXHEWPjhlcM2UeFzw4k8VrN4ddmkjMKPQl4Rz0pWye/9GR/Pqsg1m8dgtfvfctrn9uPhu2VIRdmkiza3Dom1kPM5tuZu+b2SIzuypo72hmr5nZsuC1Q9BuZnavmRWZ2XwzG9lUP4TIgUpKMr41uhdv/vQrfO+oPjxbWMzxd/6HCdOLdJWPtGqN2dOvAq5196HAaOByMxsKXA+87u4DgNeDcYDTgAHB13jg/kasW6RJ5GSkcfMZQ3ntmuM4qn8ud76ylOPunM5j76xiZ5XCX1qfBoe+u6919znB8BZgMZAPjAUmBbNNAs4KhscCj3rUTCDHzLo1dP0iTalPbiYTvxNh8vjR9OyYwc0vLeKE37/J5Hc/olKduEkr0iTH9M2sNzACmAV0cfe1waR1QJdgOB9YXWex4qBt9+813swKzKygtLS0KcoT2W+j+3Ziyg+P4NHvjyK3XTrXP7+AMXe9yROzPtRhH2kVGh36ZpYFPAdc7e6fuwzCo9fDHdA1ce4+0d0j7h7Jy8trbHkiB8zMOHZgHi9ediQPjYuQk5HKTS8s5Og73mDC9CLKt1eGXaJIg6U0ZmEzSyUa+E+4+/NB83oz6+bua4PDN7VPtygBetRZvHvQJhKXzIwxQ7pwwuDOzFxRxgMzlnPnK0uZML2Ib365B98a3Yt+eVlhlylyQBoc+mZmwEPAYnf/Q51JU4FxwG+D15fqtF9hZpOBw4HyOoeBROKWmXFEv04c0a8Ti9du5sEZK3h85of87X+rOLp/Lt8a3YsTh3QmJVlXQEv8s4bekWhmRwNvAQuA2jNdNxI9rj8F6Al8CJzn7mXBh8R9wKnAduB77l6wt3VEIhEvKNjrLCKhKN2ykykFq3li5oesKa+gW3Ybzh3ZnXMP606f3Mywy5MEZ2aF7h7Z47R4vg1doS/xrqq6hulLS3l85oe8tayUGoeRPXM497DunDHsS2RnpIZdoiQghb5IDKzfXMGL75Xw3JxiPli/ldRk45gBeZx2cFdOHtpVHwASMwp9kRhydxaWbOaluSX8a+E6SjbtICXJOKp/Licf1IWvDOpMfk7bsMuUVkyhLxISd2d+cTnTFqxl2sK1rC6LPrR9QOcsvjIoj+MHdSbSuyNpKToJLE1HoS8SB9yd5aVb+c/SUqYv3cC7K8uorHYy05I5vG8nRvXpyKg+HRmWn02qrgSSRthb6DfqOn0R2X9mRv/O7ejfuR2XHNOXbTureHv5RqYv3cDMFRt5Y0n0lpa2qcmM7JXDqN6diPTuwMH52WS31fkAaRra0xeJE6VbdjJ7VRnvrixj1soylqzbTO2fZ+9OGQzrnsMh+dkM657NQV9qT7s2+iCQPdPhHZEWqHx7JfOKN7GgpJwFxeUsKCmnZNOOT6d/KbsNA7q0Y1DXdgzonMXALu3o3zmLzHT9A5/odHhHpAXKzkjl2IF5HDvwsz6oNm7dyYKSchat2cyy9Vv4YP1W3lmxkV1Vn/UEmpuVTq9OGfTqmEGPjhnR4U4Z9OyYSW5WGtH7JCVRaU9fpIWrrnE+KtvOB+u3ULRhKx9t3M6HZdtYXbaDNeU7qPsnnp6SRJf2bejavg1dstvQtX16dDw72tYpK52OGWm0b5uiD4cWTHv6Iq1YcpLRJzeTPrmZnHLQ56ftrKqm+JMd0Q+CjdtYU17BuvIK1m2uYH7xJl4tr2Bn1RefF5CcZHTISKNjZiodMtLolJUWjKeR3TaVrPQUstqkkJWeQrs2qbQLhrPapJCVlkJSkj4w4pVCX6QVS09Jpl9eVr29gbo75TsqWVtewfrNFZRt20XZtl18sn3XZ8PbKvlg/dZP2/fn4EBWegqZ6cm0TU2mzadfSdHXlGTapkXH01M+m9Y2NZnU5CRSk43U5CRS6g4n1bbZp/OkJCV9NhzMk5RkJBkkm2FmJAfj0XYL2gnao9MS7T8ahb5IAjMzcjLSyMlIY0i39vucv6bG2barii0VVWzdGX3dUlHJ1p1VbA3aNlfUDldSUVnDjspqKiqr2VlZQ9m2XezYVU1FVTUVlTWftu8K8elkFnxIJO3hAyEpyTA++2CwYP7oUHT4szYwrM7w5z9QzGrn333Zz743dZYd0q09913Y9I8SV+iLyH5LSrLgcE7TXi5aXeNUVFZTWV1DZbVTWV1DVbVTWRO8VtdE22qcyqoaKmucqqCtstqpCuZzh2p3atypqXFqPPq9a2rbgnF3p7qGOu3Rca8zXNte+5+NEx2u/Ucn2l5nukfn+WyYzy3L55b13aZH22qHcejZMaNJt3Ethb6IhC45yXSpaYzoXm8RkQSi0BcRSSAKfRGRBBLz0DezU81sqZkVmdn1sV6/iEgii2nom1kyMAE4DRgKXGBmQ2NZg4hIIov1nv4ooMjdV7j7LmAyMDbGNYiIJKxYh34+sLrOeHHQ9ikzG29mBWZWUFpaGtPiRERau7g7kevuE9094u6RvLy8fS8gIiL7LdZ3Q5QAPeqMdw/a9qiwsPBjM/uwEevLBT5uxPLNRXUdmHitC+K3NtV1YOK1LmhYbb3qmxDTrpXNLAX4ABhDNOxnAxe6+6JmWl9Bfd2Lhkl1HZh4rQvitzbVdWDitS5o+tpiuqfv7lVmdgXwCpAMPNxcgS8iIl8U884u3H0aMC3W6xURkTg8kdvEJoZdQD1U14GJ17ogfmtTXQcmXuuCJq4trh+XKCIiTau17+mLiEgdCn0RkQTSKkM/Xjp1M7MeZjbdzN43s0VmdlXQfquZlZjZ3ODr9JDqW2VmC4IaCoK2jmb2mpktC147xLimQXW2y1wz22xmV4exzczsYTPbYGYL67TtcftY1L3Be26+mTX9c+72XtedZrYkWPcLZpYTtPc2sx11tttfmquuvdRW7+/OzG4IttlSMzslxnU9XaemVWY2N2iP2TbbS0Y03/vM3VvVF9FLQZcDfYE0YB4wNKRaugEjg+F2RO9RGArcCvwkDrbVKiB3t7bfAdcHw9cDd4T8u1xH9EaTmG8z4FhgJLBwX9sHOB34F9HHm44GZsW4rpOBlGD4jjp19a47X0jbbI+/u+BvYR6QDvQJ/m6TY1XXbtPvAn4R6222l4xotvdZa9zTj5tO3dx9rbvPCYa3AIvZra+hODQWmBQMTwLOCq8UxgDL3b0xd2U3mLvPAMp2a65v+4wFHvWomUCOmXWLVV3u/qq7VwWjM4ne7R5z9Wyz+owFJrv7TndfCRQR/fuNaV1mZsB5wFPNse692UtGNNv7rDWG/j47dQuDmfUGRgCzgqYrgn/PHo71IZQ6HHjVzArNbHzQ1sXd1wbD64Au4ZQGwPl8/g8xHrZZfdsnnt533ye6N1irj5m9Z2ZvmtkxIdW0p99dvGyzY4D17r6sTlvMt9luGdFs77PWGPpxx8yygOeAq919M3A/0A8YDqwl+q9lGI5295FEn29wuZkdW3eiR/+fDOWaXjNLA74GPBM0xcs2+1SY26c+ZnYTUAU8ETStBXq6+wjgGuBJM2sf47Li7ne3mwv4/M5FzLfZHjLiU039PmuNoX9Anbo1NzNLJfrLfMLdnwdw9/XuXu3uNcCDNNO/tPvi7iXB6wbghaCO9bX/LgavG8KojegH0Rx3Xx/UGBfbjPq3T+jvOzP7LnAGcFEQFASHTjYGw4VEj5sPjGVde/ndxcM2SwHOAZ6ubYv1NttTRtCM77PWGPqzgQFm1ifYWzwfmBpGIcGxwoeAxe7+hzrtdY/BnQ0s3H3ZGNSWaWbtaoeJnghcSHRbjQtmGwe8FOvaAp/b+4qHbRaob/tMBb4TXF0xGiiv8+95szOzU4HrgK+5+/Y67XkWfWIdZtYXGACsiFVdwXrr+91NBc43s3Qz6xPU9m4sawNOBJa4e3FtQyy3WX0ZQXO+z2JxhjrWX0TPcH9A9BP6phDrOJrov2XzgbnB1+nAY8CCoH0q0C2E2voSvXJiHrCodjsBnYDXgWXAv4GOIdSWCWwEsuu0xXybEf3QWQtUEj12enF924fo1RQTgvfcAiAS47qKiB7rrX2f/SWY99zg9zsXmAOcGcI2q/d3B9wUbLOlwGmxrCtofwS4dLd5Y7bN9pIRzfY+UzcMIiIJpDUe3hERkXoo9EVEEohCX0QkgSj0RUQSiEJfRCSBKPRFRBKIQl9EJIH8fweuTDcx7S92AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(penalty)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(1, 1), dtype=float32, numpy=array([[21.093924]], dtype=float32)>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.transpose(trainable_weights-prior) @ B_inv @ (trainable_weights-prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=39.35813>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.reduce_sum(tf.transpose(y-model_x(trainable_weights),perm=[0,2,1]) @ R_inv @ (y-model_x(trainable_weights)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Variable 'Variable:0' shape=(2, 1) dtype=float32, numpy=\n",
       "array([[0.313344  ],\n",
       "       [0.09426614]], dtype=float32)>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trainable_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "custom_loss(prior, x_pred)\n",
      "  Args:\n",
      "    prior: float32 Tensor, shape=(2, 1)\n",
      "    x_pred: VariableSpec(shape=(2, 1), dtype=tf.float32, name='x_pred')\n",
      "  Returns:\n",
      "    float32 Tensor, shape=(1, 1)\n"
     ]
    }
   ],
   "source": [
    "print(custom_loss.pretty_printed_concrete_signatures())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NEDJBrBShVeS"
   },
   "source": [
    "## Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 283
    },
    "id": "QZ8onHV8hVeU",
    "outputId": "34033b88-315d-4b74-d338-16d392c5e432",
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1610ede80>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABH7klEQVR4nO3dd3hUVfrA8e+ZmVRIIQkJJASCQOhSDEUpgkgRaYoKigUb7q7dVRd297frVuuuZXUtrL2googgShcFRDTUAGnUkFBSSAKkZ+b8/riTGMgkmcAkkxnez/PMMzN37sy8c5O8OXPOe89RWmuEEEJ4F5O7AxBCCOF6ktyFEMILSXIXQggvJMldCCG8kCR3IYTwQhZ3BwAQERGh4+Li3B2GEEJ4lC1btuRqrds6eqxFJPe4uDgSExPdHYYQQngUpdShuh6TbhkhhPBCktyFEMILSXIXQggvJMldCCG8kCR3IYTwQi2iWkY0jcXbsnh2RSpHCkqIDg3gsfHdmTYgxt1hCSGagWoJs0ImJCRoKYV0rcXbspi3KImSCmv1Nh+TorW/hYLiCkn2QngBpdQWrXWCo8ekW8ZLPbsi9YzEDlBh0+QXV6CBrIISHvpkOwP+upLF27LcE6QQoslIcvdSRwpKnNovv7iCeYuSJMEL4WUkuXup6NAAp/ctqbDy7IrUJoxGCNHcJLl7qcfGdyfAx+z0/lkFJdJ6F8KLSLWMl6oaKH3ok+1OP2feoqQzniuE8FzScvdi0wbEECPdM0JckCS5ezlH3TM+JlXn/s4OxAohWjZJ7l5u2oAYnry2LzGhASggJjSAZ6/vV2eL3qSU9L0L4QWkz90LOTozdePcK2rtd/ZJTgBWraXvXQgvIC13L1N1ZmpWQUn1yUqO6tirWvRmVbuLRvrehfB8kty9jKMzU+tK1tMGxGCrY/oJ6XsXwrNJcvcydSXlurZ3DlQ8sHEBP718Cy8teYauuRlA406CEkK0PJLcvUxdSdnh9vx8lsz/DY9s+JC94bGM2fsTK966jwkHt/DY+O5NHKkQoilJcvcyjkofA3zMjpP1E0/Q+shhNr7yIY/9+nlG/OpN9kfF8cLql5kW69dMEQshmoIkdy/jqPTxyWv71q582b0bXnkF5sxh2G9uYuPcK/jT7Zfztxvmogry+Xb0tSzemumWzyCEOH9SCumFpg2IabiM8eGHISgI/vY3oMb87wHRPD/8ZuZ+9w53/Ot9ePRWKYkUwgNJy/1ClJoKq1bBvHkQEQGcWWXzdsIUcgJDufnHRVISKYSHkuR+IVqwAJSCm2+u3lSzmqbM4ssHAyZyxf5E/PeluSNCIcR5kuR+odHaSO6jRkF0dPXms6tpPhgwkTKzD/clfd3MAQohXEGS+4Vm2zZIS4OZM8/YfHaVTV6rUL7qM5rJO1ZBQUEzBymEOF+S3C80CxaAxQLTp5+x2VGVTfhv78dSWgqLF7slVCHEuVO6jtPPm1NCQoJOTEx0dxgez9GEYWdUumgNF10EvXvDV181/IJaQ5cu0L07fPNN0wUuhDgnSqktWusER49Jy91LODVh2IEDcPAgTJzo3IsqBTNmGJU1ublNEbYQoolIcvcSTk0Ytm6dcT1qlPMvPGMGWK2waNF5xyiEaD6S3L2EUxOGrVsHbdtCz57Ov3C/fka3zMcfn1+AQohmJcndSzQ4YZjWRnIfNcrobnFWVdfMunXSNSOEB3EquSulHlZK7VZK7VJKLVBK+SulOiulNiul9iqlPlFK+dr39bPf32t/PK5JP4EAnJgwbP9+OHwYRo+u93UWb8ti2FNr6Tx3GcOeWmv02U+ebPxzWL68qcIXQrhYg8ldKRUDPAAkaK37AGZgJvA08LzWuiuQD9xpf8qdQL59+/P2/UQTa3DCMCf62+sclFVREBUFX8sJTUJ4CmcnDrMAAUqpCiAQOApcAdxkf/xd4AngVWCq/TbAZ8DLSimlW0LNpZerd8KwdesgMhJ69Kjz+XUOyq5KZ9pVV8GXX0JlpVEnL4Ro0RpsuWuts4DngAyMpF4IbAEKtNaV9t0ygaqsEgMctj+30r5/+Nmvq5Sao5RKVEol5uTknO/nEA3ZsAFGjKi3v73eQdmJEyE/HzZvbqoIhRAu5Ey3TBuM1nhnIBpoBUw43zfWWr+htU7QWie0bdv2fF9O1Ccvz6hvHzy43t3qHZQdNw7MZli2rAkCFEK4mjMDqlcCB7TWOVrrCmARMAwIVUpVfT/vAFSdLZMFxALYHw8B8lwatWicrVuN60suqXc3R4Oy/haTMSgbEgLDh0tyF8JDOJPcM4ChSqlApZQCxgB7gG+B6+z73AZ8ab+9xH4f++Nrpb/dzaqmdhg4sN7dqgZlo0P8AbCYFOVWGxv35rI/5zRMmAA7d8Lx400dsRDiPDnT574ZY2B0K5Bkf84bwO+AR5RSezH61N+0P+VNINy+/RFgbhPELWpwWL5Y05Ytxhwxbdo0+FpT+0cztIsxRNK3QwjXXdKBpTuPcNWL69nbZ4ix09q1rv4IQggXc6rsQWv9Z+DPZ23eD9TqxNValwLXn39owhnVy+PZq1yqyheBXypnEhNh6FCnXu+DHw+xaGsWD47pxsNj4wH47bjuXPvfH7htRwXrQ0MxrVkDN97o+g8jhHAZOUPVwzU4p0xuLhw61GB/O0Da8VP8ZekexvSI5MEx3aq3RwX78/otl5BXaiUxrh96zRqXfgYhhOtJcvdwDc4ps2WLcZ3gcFbQM7y2bh++FhPPXt8Pk+nMksk+MSH8aVJvlrbtiTp40DjjVQjRYkly93ANzilTldwbGEw9UlDCkh1HmDmoI2GtfB3uc0NCBw4PuBQA28pV5xawEKJZSHL3cA3OKbNlC3TtapQy1uOtDQfQwJ0jOte5j8Vs4oZZV3K0dThZi5xY7EMI4TaS3D1cg3PK7NwJ/fvX+xonSytY8FMGU/pFE1PHN4EqE/q0J7lnAq02bcRmtbnmQwghXE4mCfECdc4pU1wM+/bBzTfX+/zVe45TVG7llks7NfheJpOizVVXEvbzCn5e+SODrrrsXMMWQjQhabl7s+RkY6revn3r3e3rpGNEh/gzIDa0elt9tfO9Z04CIOVT6ZoRoqWS5O7Ndu0yrvv0qXOXU6UVfJ+ew1V926Psk4o1tB6rb4/unG4TQdDPmzha6LhaRwjhXpLcvVlSEvj5GWen1mFtSjbllTYm9m1Xva3B2nmlUCMvZ3DGLj7enNEkoQshzo8kd2+2axf06mXM5liHr5OOEhXsx4DYX6YmcGY91lZjRxN9KpcNa7YgUwcJ0fJIcvdmu3bV2yVTUm5lXWoOE3q3O+OkpQZr5wFGjgSg0+5EtmYUuCRcIYTrSHL3YPVOGJafD1lZ9Q6m/nzwBGWVNq7oGXXG9gZr5wF690a3acOlmbtZuuOISz6PEMJ1JLl7qIYGPZ0ZTN24Lxcfs2JQ3JmzRTZYOw9gMqFGjGDksRS+2nkUq026ZoRoSaTO3UPVN+g5bUCMU8l90748BsS2IdC39q9BveuxVhk5kqglS1DHjvLj/jyGdY1o9OcQQjQNabl7qAYHPXftguBg6NDB4X6FxRUkZRVyaZday9s6z97vPuJosnTNCNHCSHL3UA0OeiYnG5UydSyIvWl/Hlpzfq3tAQOgVSumntrH6uTj0jUjRAsiyd1DNTjomZwMPXvW+fxN+3IJ8DHTv8ZZqY1mscCwYQw4mETu6XK2H84/99cSQriUJHcPVe+gZ0EBHDsGPXrU+fyN+/IY1DkMX8t5/gqMHElwejIR5adZuVvWVhWipZABVQ9W56Bnqv1M0jqS+4micvZmn+bagQ0MmDrD3u9+U0UGX+2JYt7Eur8tCCGaj7TcvVFKinFdR3LflmF0nwzs2PCC2Q0aNAj8/JiQm8r+3CL2Zp8+/9cUQpw3Se7eKCUFfHzgooscPrw1Ix+zSXFxh/oX8HCKvz8MGULXtO0ArNxz7PxfUwhx3iS5e6OUFOjWzRjwdGBbRgE92wc5rG8/J8OH47tjOwPDfViXkuOa1xRCnBdJ7t4oObnOLhmrTbPjcIFrumSqjBgBViszrFlsycjnZGmF615bCHFOJLl7m4oKY/WlOpJ72vFTFJVbGdAx1KmXq3f+miqXXgpKMexYKlabZmN67nl8ACGEK0hy9zb79kFlZZ017lsbMZja4Pw1VUJCoF8/opMSCfK38F2adM0I4W6S3L1NA5UyWw8VEN7Kl45hgQ2+VIOLdtQ0YgSmHzcxMi6U79JyZI53IdxMkru3qUru3bs7fHjb4XwGdAytXlKvPs4s2lFt+HAoLmaayuZoYSlpx6UkUgh3kuTubVJTIToagoJqPXSqtIL9OUX06xDq1Es5tWhHlREjABiSuQeA79KynYtXCNEkJLl7oHoHOdPSID7e4fN2HzkJQB8n69udWrSjSvv20KULwYk/0j0qiHWp0u8uhDtJcvcwDQ5y1pPcd2UVAtAn2rnkXnP+GgCzUtV97g6rZoYPhw0bGBUfwc8HT1BUVtnozyeEcA1J7h6m3kHOEycgN7fe5N4u2J+2QX5Ov9+0ATHVLXirfZC0zqqZESMgN5cJlkIqrJof9uU17sMJIVxGkruHqXeQMz3duFNHck/KKqRPTOOnHHC6amb4cAD6Hkwi0Ncs/e5CuJEkdw9T7yBnWppxx0FyP11Wyf7cIvqeQ3J3umomPh4iI7H8sJHLukSwLlVKIoVwF6eSu1IqVCn1mVIqRSmVrJS6VCkVppRapZRKt1+3se+rlFIvKaX2KqV2KqUGNu1HuLDUO8iZng5mM3TuXOt5e46cRGvo2yG40e/pdNWMUkbrff16RnVvS2Z+Cftzixr9fkKI8+dsy/1FYLnWugfQD0gG5gJrtNbdgDX2+wBXAd3slznAqy6N+AJX7yIdaWlGYvf1rfW8pKrB1HNouTeqamb4cDhwgCtaG/PLSNWMEO7R4LSASqkQYCQwG0BrXQ6UK6WmAqPsu70LrAN+B0wF3tPG9/Ef7a3+9lrroy6P/gJV5yId9ZVBZhUSGeRHZJD/Ob0fGH3vRwpKiA4N4LHx3R3HYK93j969hYsiotiQnsOdw2t/kxBCNC1n5nztDOQAbyul+gFbgAeBqBoJ+xgQZb8dAxyu8fxM+7YzkrtSag5Gy56OHTuea/yiitZGcrevjHS2XUfObTC1Sp3/UM7Wvz+0agXr1zNi7K/4NDGTskorfhZzg08VQriOM90yFmAg8KrWegBQxC9dMADYW+mNGjnTWr+htU7QWie0bdu2MU8Vjhw9CkVFDlvupRVW9uUU0at94/vbz9bgLJEWizFL5IYNjOjWlpIKK1sOycLZQjQ3Z5J7JpCptd5sv/8ZRrI/rpRqD2C/rqp7ywJiazy/g32baEr1VMrszT6N1abpeZ7J3elZIkeMgJ07GRpuxmJSrJcpgIVodg0md631MeCwUqpq9GwMsAdYAtxm33Yb8KX99hLgVnvVzFCgUPrbXafOlnM9yT3l2CkAerSvPd9MYzSq3l1rWm/5iYGd2vC9TAEsRLNzdp21+4EPlVK+wH7gdox/DJ8qpe4EDgE32Pf9GpgI7AWK7fsKF6hqOVcl2KqWM8C0tDRjPdMOHWo9L/noSfx9TMSFtzqv93e63n3oUKN7ZsMGRl55O8+tTCP3dBkRrZ0/M1YIcX6cKoXUWm+3949frLWeprXO11rnaa3HaK27aa2v1FqfsO+rtdb3aq27aK37aq0Tm/YjXDjqbTmnpRnrpppq/0hTjp2ke1QQZlPD0/zWx+l698BAuOQSWL+ekfHGeMrGvdI1I0RzkjNUPUiDUw846JLRWpN89BQ92p3/YGqj691/+oneYX60CfTh+zRJ7kI0J0nuHqSulnNskK+xvJ6D5J5zqowTReX0PM/+dmjgBKqzjRgB5eWYt25hWNcI1qfLVARCNCdn+9xFC/DY+O5n9LmD0XL+Y99AY2FsB8k9uXow9fxb7tCIevdhw4zr779n5Nhb+GrnUdKOn6Z7u/P/JyOEaJi03D1IXS3ncT7GIhx061brOclHjcd6uqBbpqYG690jIqBvX1i3juHdIgBYny5VM0I0F2m5exiHLecXPzOuHZVBHj1JdIg/IYE+Louh3qqdmrGNGgVvvkl0oJmuka35Pj2Xu0Zc5LI4hBB1k5a7N0hLg9BQo7V8luSjp1zWJVPF6Xr3UaOguBgSExnRLYLN+/MoPet5QoimIcndG1RNGKbOLHUsq7SyL+e0SwZTa3K63v3yy43rdesY2a0tZZU2Eg/KVARCNAdJ7h7iXBbF3pddRKVNu6QMsian693Dw+Hii+HbbxlyURi+ZhPfS7+7EM1CkrsHqHdOl5ISyMhwXClTNZjq4m6ZRtW7jxoFGzcSiI2EOJmKQIjmIsndA9Tbx71vn7HB4ZwyJ/GzmIgLD3RpPI2qdx81yvgH9PPPjOjWlpRjp8g+WerSeIQQtUm1jAeot4+7ngnDko+eIj4qCIvZ9f/Dna53HznSGAtYt44Rs+/j6eWwYW8u1w6sPQeOEMJ1pOXuAert4061V6g4qHFPOXbS5YOpZ2uw3r2q333dOnq1Dya8la9MASxEM5Dk7gHq7eNOSYGYGGjd+ozHc06VkXu63OWDqTU5Pb+7vd/dVFHO8G4RrE/PxWaTqQiEaEqS3D1AvX3cqanQvfZAZlMNptbUqHr3Gv3uuafLqueYF0I0Delz9xAO+7i1NpL7jTfW2j/lmJHcezThXC5O17vX7Hd/4FHAmIqgV3TT/eMR4kInLXdPlpMDBQV1tNxP0S7YnzatfJvs7Z2udw8Lg379YN06ooL96R4VJP3uQjQxSe6erGowtY5umaYeTG10vfsPP0BZGSPjI/jp4AlKymUqAiGaiiR3T1aV3Hv0OGNzeaWNfTmnXT6nzNnOqd5982ZGdGtLeaWNzQfymjQ+IS5k0ufuyVJSjHVTO3Y8Y/O+nNNUWHWTDqZWcbrefdQoMJth1SoG/2kYvhYT69NzGdU9ssljFOJCJC13T5aa6nDd1KrB1J4taWGMkBAYMgRWrsTfx8yQzmEyv7sQTUiSuyerswzyFL4WE50jWrkhqHqMGwc//wwnTjCiWwRpx0/XWXEjhDg/ktw9VXk57N9f52BqfFTrJpl24LyMG2eUb65Zw2h7d8y3qdluDkoI7yR97i3U4m1ZPLsilayCEsxKYdWamNAAHhvf3ejj3r8frFaHyT3l2Ckuj2/rlniPFJQQXTPOmgYNMrpnVq6k63XX0aFNAN+mZDNrSKdmjVWIC4Ek9xbo7GXsrNo4Vf+M5ewyHJdB5p4uI+dUWbMMptYVb53L7lksMGYMrFyJAq7oEcnCxExKK6z4n1VSKYQ4Py3se7sAx6f1V6k+vb+OGveUo8Zp/c05mOr0NARgdM1kZEBaGqN7RFJSYeXH/VISKYSrSXJvgRoaZDxSUGKUQbZrZ3Rz1FA1p0xT17jXisfZ7ePGGdcrV3LpReEE+Jj5NkX63YVwNUnuLVBdp/Wf8Xg9E4ZFBfsR1oTTDjiMx9ntnTtD166wahX+PmaGdQ1nbWo2WssskUK4kiT3FsjRaf1Vqk/vryO57zl6kt7RIQ6e2XQaNQ0BwNix8O23UF7O6B6RHD5Rwt7s080QqRAXDknuLVDN0/oBzEoBNU7v7+gPeXm1kntphZW92afp1YxdMnXFW9XnXmtudzC6Zk6fhh9/rC6JXCtdM0K4lFTLtFD1ntb/ww/G9VnJfW/2aSpt2i1T6VbF6lTVzOjRxlQEK1cSPXIkPdoFsTYlm3su79LscQvhraTl7onqqJTZc8QYTG3ulnsVp6tmQkJg6FBYuRIwSiITD+VTWFLRXKEK4fUkuXui1FTw8YG4uDM27zl6kla+ZjqGBbolrEZXzSQmQk4OV/SIxGrTMteMEC4kyd0TpaQYE4ZZzuxV23PkJD3bB2MyKbeE1aiqmauvNqYi+PprBnRsQ2igj/S7C+FCktw9kYNKGZtNs+foSbcuXdeoqpmBAyE6GpYswWxSjIpvy7rUHCqttmaKVgjv5nRyV0qZlVLblFJf2e93VkptVkrtVUp9opTytW/3s9/fa388rolivzBVVsK+fbWSe2Z+CafLKt3W3w6NrJpRCqZMgRUroLSU8b3bcaKonJ8P5rshciG8T2Na7g8CyTXuPw08r7XuCuQDd9q33wnk27c/b99PuMqBA1BRUXsw9WghgNsXnZ42IKa6BX/2nDi1EvyUKVBUBOvWcXn3tvhZTKzYfcwNUQvhfZxK7kqpDsDVwP/s9xVwBfCZfZd3gWn221Pt97E/Psa+v3CFPXuM6549z9x85CRmkyI+yv0LdDhdNTN6NLRqBUuWEOhr4fL4tizfdQybTc5WFeJ8OVvn/gLwOFCVOcKBAq11pf1+JlBVyBwDHAbQWlcqpQrt+5+x3L1Sag4wB6DjWcvEiXrs3m1c9+p1xuY9R0/SpW2rFjG7otNVM/7+RtXMkiXwyitM6NOOlXuOszOrkP6xoU0fqANaa9Ly0tiVvYt9+fsoqShBKUVMUAxdw7oyKGYQgT7uqUYSojEaTO5KqUlAttZ6i1JqlKveWGv9BvAGQEJCgjTVnLV7N8TGQtCZLfQ9R04yuHOYm4I6U3RoAFkOErzDqpkpU+CLL2DbNsb06IvFpFi+61izJ/e0vDRe/flVFqcu5mDBwTr38zX7Mix2GLf2u5Xrel1Ha9/WzRekEI3gTLfMMGCKUuog8DFGd8yLQKhSquqfQwegqkM1C4gFsD8eAsicrq6yezf07n3Gpvyico4Ulrq9v71Ko6pmrr7aGFxdsoSQQB8u7RLO8l1Hm20isfS8dK795Fp6vNyD/yb+l76RfXl90utsmbOFk3NPYv2TlfI/lnPgwQMsu2kZDw55kMyTmdz+5e3EPh/LX7/7K4Wlhc0SqxCN0WBy11rP01p30FrHATOBtVrrWcC3wHX23W4DvrTfXmK/j/3xtVqm/HMNq9WocT8ruVdN89urffNOGFaXRlXNtG0Ll11mdM0AE/q042BeManHTzVpjOXWcn6/5vf0/m9vVu9fzR9H/pFDDx1iyY1LmHPJHAa2H0iQXxAmZcLH7ENcaBwTu03kmbHPkHpfKutvX8/lnS7nz+v+TNf/dOXtbW9j01LGKVqO86lz/x3wiFJqL0af+pv27W8C4fbtjwBzzy9EUW3fPigrq5Xc99iTe8/27h9MrdLoqplt2yAzk7G9olAKlu9quqqZtLw0Ln3zUp7c8CQ39b2JtPvT+Ovov9KudTunnq+UYnjH4SyeuZgtc7YQHx7PHUvuYMx7Y8gozGiyuIVojEYld631Oq31JPvt/VrrwVrrrlrr67XWZfbtpfb7Xe2P72+KwC9IVZUyZyf3IydpF+xPeGs/NwRVN6erZqZMMa6XLiUyyJ+ETm2aLLmv2LuCQfMHcbDgIItnLOadae84ndQdGdh+IOtvX8/8yfNJPJLIxa9ezOd7PndhxEKcGzlD1ZPUUynTUvrba3K6aqZ7d2MBjy+Nnr0JfdqTcuwUB3OLXBrP/C3zufqjq4kLjWPbPduY2mOqS17XpEzcNfAutt+znR4RPbhu4XXMWz0Pq83xUolCNAdJ7p5k927o1Ala/1Kh4a453J3h9FwzSsHUqbB2LRQUML53FIBLT2h6ftPzzPlqDuO6jGPD7RvoGOL68tsuYV34bvZ3zBk4h6c2PsX0T6dTUlH/kolCNBVJ7p5k9+5arXZ3zuHekEZVzVx/vXHm7ZIldGgTSN+YEJYlHXVJHM/98ByPrHyE6T2ns3jmYoL8mm5sws/ix+uTX+c/V/2HJalLmPDhBKmmEW4hyd1TVFY6rJTZfcRIHD1bYMu9qmomNMCnepu/Tx2/coMHQ8eOsHAhAFP6RbMzs5D9Oee3/N7/tv6Px1Y9xg29b+Dj6z7G19w8a8veN/g+Ppr+EZsOb+Lydy7n2GmZVkE0L0nunmLvXigvr5Xcd2QWEuRvoZOb5nB3RlnlLyWC+cUVjitmlILrrjMW8CgsZEr/aJSCxduPnPP7Lk5ZzD1f3cOErhN4/5r3sZiad+GxmX1m8tVNX5F+Ip1hbw3jQP6BZn1/cWGT5N5CLN6WxbCn1hI3dxld5n1N3NxlDHtq7S9JcMcO47pfvzOel5RZyMUdQtw2h3tDnK6YAaNrprwcliwhKtify7qEs3hb1jmd0LTj2A5mLZrFoOhBfH7D583WYj/buC7jWHvrWvJL8rnivSs4XHjYLXGIC48k9xZg8bYs5i1Kqj5l32Fd+I4dxuIcNfrcyyqtpBw7Sd+YUHeE7ZRGrc40ZIgxtcKnnwIwrX8MGSeK2Xa4oFHvmVOUw9SPp9LGvw1fzPjC7XPBDOkwhJW3rOREyQnGvDeGo6dcM5YgRH0kubcAjlq3VapbuTt2QI8e4PdLLXvK0VNUWDX9OrSMM1MdqatixqSU466ZG24w5ng/cYIJfdrhZzHV3q8e5dZyrlt4HceLjvPFjC9oH9T+fMJ3mYToBL6Z9Q1HTh3hyvevJKdIlhQUTUuSewtQV+v2jMd37KjVJbMzswCAvi04uTuqmAHj24nDvvdZs4yqmYULCfL34cqeUXy18ygVTq7Q9OA3D/L9oe95c8qbDIoZ5IqP4DKXxV7GVzd9xf78/Yz7YBwFpQXuDkl4MUnuLUBdrdsqPXzKISvLQXIvJLyVb/UcLi1RVcWM2cGU/g773vv3N+aq//DD6uefKCp3avHst7a9xWtbXuPxyx7npr43uSJ8lxsVN4rFMxazO3s313xyDWWVZe4OSXgpSe4tQF2tWzDqwv8QW2HccZDcL+4QQktfC2XagBhsdQyKZhWUnNl6V8pova9fD4cOcXl8W0IDfVi8rf6qmeScZO77+j7GdB7DP8f805Xhu9z4ruN5e+rbrDu4jtlfzpYJx0STkOTeAjiaRREgJjSAJ6/ty/Bie2KrkdyLyytJzz5F3w6hzR3uOanv20mt7pmb7K3uBQvwtZi4um97Vu45xumySofPL60sZebnM2nl24r3r3kfs8n9C5Y0ZNbFs3hqzFN8vOtj5q6WufWE6zVv4a+o07QBMUwbEOP4wZd2QFSUcbHbfeQkNk2LHkyt6bHx3Zm3KMnhwHFV90z15+/c2ZgG+L334He/Y9qAGD7cnMHK3ce4dmCHWs9/dOWj7Dy+k2U3LWsxA6jOeHzY42QUZvDsD88SGxzL/UPud3dIwotIy90TOBhM3ZaRD0A/Ny1H11hV307qUmtQefZsSE6Gn37iko5tiA0L4LMtmbWe92XKl7zy8ys8PPRhJnab6OKom5ZSipeueomp3afy4PIHWZS8yN0hCS8iyb2lq6gw5pQ5K7lvPVRAp/BAIlrYNL/1mTYgps7B31rdNjNmQEAAvP02JpNiRkIsP+zL41DeLzNFZp3M4o4ldzCw/UCeHPNkU4beZMwmMx9N/4ghHYYwa9EsNmZsdHdIwktIcm/pdu82ztocOLB6k9aaLRn5DOzYxo2BnRunJxMLDobp0+Hjj6GkhOsuicWk4JOfjTM8tdbcvfRuSipKWDB9AX4Wz/knd7ZAn0CW3riU2OBYpnw8hfS8dHeHJLyAJPeWLjHRuE5IqN6UmV9CzqkyBnYMdU9M56Hm4LHil0Fjh+MNt98OhYXwxRe0C/FndPdIFm7JpNJq4+3tb/PN3m94+sqniQ+Pb/bP4WoRgREsv3k5ZmVm4kcTySuWZYfF+ZEB1ZYuMRFCQqBLl+pNW+397QM8sOUODQwe1zRqFMTFwfz5cNNNzBgUy5qUbD7eupWHVj/EqLhR3Dv43qYOt9lc1OYivpz5JaPfHc20T6ax+pbVHv2NRLiXtNxbusREo9Veo5Z9W0YBgb5merRrOWumnouqydI6nz1JWhWTCebMgXXrICWF0T0iiQzy5fE192LTNt6a8hYm5V2/wpfGXsq7095lQ8YG7lhyxzlNmiYESHJv2crKYOdOGHTmafRbM/K5uEMIFrPn/vhqTpamqWfx7DvuMCZMe+MNfMwm4jpt4mjZT8y77B90btPZLbE3tRl9ZvDPK/7JR0kf8cS6J9wdjvBQnpsdLgQ7dxrVMjX620srrOw5ctIjB1Nrcnoq4KgouPZaeOcdDhxNZlnGMwTY+mM7dWUzRtv85g6fyx397+Cv3/+V93a85+5whAeS5O5m9XZNOBhM3X64gEqb5pJOnp3cGzUV8K9+ha0gnzveuxazMnHdRX9n4ZZMSsq9dwFqpRSvTXqNMZ3HcNeSu1h3cJ27QxIeRpK7GzXYNZGYCBERxvJzdpv3n0ApSIgLc0/QLuL04tkAo0bx30lRrCtN4fnx/+beEUM5WVrJ4u3OTwXsiXzMPnx2w2d0DevKNZ9cQ0puirtDEh5EkrsbNdg1kZho9LfXGEzdfCCPnu2CCamxLqknaszi2Xvz9/G7S/K5Kh3uKO7B4M5h9GofzJsbDnj9gGOofyjLblqGr9mXqz+6WuaBF06T5O5G9XZNnDoFu3ad0SVTXmlja0Y+Qy7y7FY7OF/vbrVZmb14Nj5+Acz/LgT14osopbhrRGf2Zp/muzTvT3ad23RmycwlHDl1hKkfT6W0stTdIQkPIMndjertmti8GWw2GDasevvOzAJKK2wM6RzeXCE2qWkDYtg49wqen9EfgIc/2V5r3OGlzS+x8fBGXrrqJWJm/QoWLYJDh5h0cTRRwX68ueHCWHR6SIchfHDNB2zK3MTsxTJNsGiYJHc3qrdrYuNGoztm6NDqxzYfOAHA4M6e33KvUt+4Q2puKr9f+3smx0/mlotvgXvvNY7JCy/gazFx22VxrE/PZc+Rk+7+GM1ieq/pPHPlM3yy+xP+b+3/uTsc0cJJcnejersmNm6Evn2Ns1PtftyfR/eoIMJa+bovaBera9zhmeV7uG3xbQT6BPLG5DeMBUliY+HGG40zVvPymDW4E618zbz23T43Rd/8Hr3sUeYMnMM/N/yTt7a95e5wRAsm0w+4mcNT8a1W+PFHuPnm6k0VVhtbDuVz3SW15zP3ZHWNOyQXfURB+WY+uvYj2rVu98sDv/sdfPABvPIKIX/6Ezdf2on53+/nkbHxxEW0aqao3UcpxcsTX+Zg4UHu+eoeOoV0YsxFY9wdlmiBpOXeEiUlGQOqZ/W3F5dbGXqRd/S3V3E07lCuDlJg+YDWtmH4V4w488E+fWDyZHjpJSgq4s7hnbGYTRdU693H7MOn131Kj4geTP90Onty9rg7JNECSXJviTba5/Sukdy/T8vFpGBYlwg3BdU0zh530FSQ6/svTLQmtOw3PLpwZ+0pCebOhbw8eP11IoP8mTkols+3ZpJVx7cAbxTiH8Kym5YR4BPA1R9dzfHTx90dkmhhJLm3RBs3QnQ0dOpUven79Bz6xYYSEujZ9e1nO3uFpgLLAipMBwgvvx8zIVi1rj3nzGWXwZgx8PTTUFTEPZd3QaF4ee1eN3wC9+kY0pGlNy4luyibCR9OoLC00N0hiRZEkntLozWsX2+02u0nLxUWV7DjcAEjurV1c3BNo2qFpjJTMictn9GqciyBtiHVj5dUWPntpzvOTPB/+QtkZ8N//0tMaAAzB8eyMPEwGXnFbvgE7pMQncCiGxaxO3s3kxdMprjiwvr8om4NJnelVKxS6lul1B6l1G6l1IP27WFKqVVKqXT7dRv7dqWUekkptVcptVMpNbD+dxBnSE+HzEy44orqTRv35WLTMLKbd3XJ1HT/mFjyfJ/HrCMIq7i71uNWrXnok+1cNG8ZcXOXMWx9GceHXg7PPAOnT3Pv6K6YTYoX11x4qxiN7zqe9695nw0ZG7hh4Q1UWCvcHZJoAZxpuVcCv9Va9wKGAvcqpXoBc4E1WutuwBr7fYCrgG72yxzgVZdH7c3WrDGur/xl1sP16TkE+Vno7yGLYZ+L77Ofp1IdJbLiEUwE1rmfzT7bQFZBCQ/2mAq5ufD880QF+3PL0E58sS2T9OOnminqlmNGnxm8evWrLEtfxu1f3i4nOYmGk7vW+qjWeqv99ikgGYgBpgLv2nd7F5hmvz0VeE8bfgRClVLtXR2411q92pgozL7yktaa79NyuaxruEfP316fFXtX8Griqzxy6SP89/qba53YVZcfo+JZEX8pFU89DdnZ/GZ0V1r5Wnjqmwtzgq17Eu7hH1f8gw+TPuSh5Q95/bw7on6NyhZKqThgALAZiNJaH7U/dAyIst+OAQ7XeFqmfdvZrzVHKZWolErMyfH++UGcYrXCt98arXZ7f3va8dNkFZQwMt47+9tPlJzgjiV30KttL/5+xd+rB1jNNSZLq8/TI29DlZSw//7HCWvly29Gd2VNSjY/7Mtt4shbpnnD5/HI0Ef4z0//4a/f/dXd4Qg3cjq5K6VaA58DD2mtzzjfWxtNhEY1E7TWb2itE7TWCW3bemfiqo/Dedy3bYP8fKMSxG7l7mMAjO0ZVddLeSytNfd8dQ/ZRdm8f837+Fv8AWOA9V839HOqBb8/vAML+k+g42cfQGoqtw+LIyY0gH9+nYzNduG1XJVSPDvuWWb3n80T3z3BsxufdXdIwk2cSu5KKR+MxP6h1nqRffPxqu4W+3W2fXsWEFvj6R3s24RdXfOp7H7ffmhrJPdVycfpHxtKZLC/e4JtQm9seYPP9nzGP674BwPbnznuXtWCD3ViauMXh91IscUPHnoIf4uJxyd0Z1fWST5NPNzgc72RSZmYP3k+M/vM5PHVj/PcD8+5OyThBs5UyyjgTSBZa/3vGg8tAW6z374N+LLG9lvtVTNDgcIa3TeCuudTKV62wjgDM8popR8pKGFnZiHjentfqz3peBIPrXiI8V3G8+hljzrcZ9qAGLb/eRwvzOhPTB0zaALktmrDO2Nvg+XLYelSpvSLZlBcG55enkJBcXlTfYQWzWKy8P4173ND7xt4bNVj/HvTvxt+kvAqzrTchwG3AFcopbbbLxOBp4CxSql04Er7fYCvgf3AXmA+8BvXh+3ZHM2nElRWRL8DO2HixOptq5ONsw7H9fKu5F5UXsSMz2YQ6h/Ke9e8h0nV/2tYNTXwwaeu5oUZ/R3OpNnpT49Dr17w8MOo0lL+MqUPhSUVPLcytY5X9X4Wk4UPr/2Q63tdz29X/pbnNz3v7pBEM2pw4jCt9QagrtGtWjMW2fvf7z3PuLxadGhArVPlR+7fiq+tEqZMqd62as9xLopoRZe2rZs7xCb1wDcPkJKbwqpbVhHZKrJRz62aZO3ZFalkFZRgVsqYRXLtfsIf+DPDfzUD/v53ev3jH9x6aRzvbjrI9IEdGODhC4qfq6oEb9M2Hln5CCWVJcwbPs+YZVN4Ne+srWvhHM3jPn7/T5SFhlXP336iqJxN+/IY2zvKq/4Q39vxHm9tf4vfj/j9Oc9mOG1ATPUxtNrL/bIKSrg7M5iMydcbJzbt3Mlvx8XTLtifuZ8nUV554dZ9+5h9+Pi6j7n54pv5w9o/MHf1XCmTvABIcneDs+dx7xjkw4RDW/GbNgXMRtJflnSUSptmar9aVaQea8uRLdzz1T1c3ulynhj1xHm9Vl3jFvf0uxHatIG77ybIovj7tD6kHj91Qc0a6YjFZOHdae/y64Rf88wPz/CbZb+RE528nMzn7iZnzOO+bh2cKjyjS+bLbVl0i2xNz/ZB7gnQxXKKcrj202tpG9iWT6//FIvp/H716pwHvsLXmA74xhvh2WcZM28ek/tF85+16YzpGUnv6BCHz7sQmJSJVya+QrBfME9vfJrCskLenvo2fhY/d4cmmoC03FuCL78EPz8YOxaAwyeKSTyUz7QBMV7RJVNpq2TGZzM4fvo4i2YsanQ/uyN1rT+rgMXxw+H66+HPf4bt2/nrlN60CfTl4U+2U3pWa/9Co5TiqSuf4qkxT7Fg1wImfDiB/JJ8d4clmoAkd3ezWmHhQiOxtzYGTpfsOALAlH7R7ozMZX636nd8e/BbXpv0GgnRCS55zcfGd3c4yq+BZ1emwauvQng43HwzbUxWnrnuYtKOn+bZFRdu9UxNvxv+Oz645gM2Zmxk+NvDOVRwyN0hCReT5O5u338PWVkwaxZgnLW5eFsWl3RqQ2xY3RNoeYo3t77Jv3/8N/cNuo/Z/We77HWnDYip85TorIISFmeUwjvvwO7d8OCDjOoeyW2XduLNDQdYvUcWtgCYdfEsVty8gqyTWQx9cyhbjmxxd0jChSS5u9uHHxotdnt/+88H80nPPs0NCZ6/Vuryvcu556t7GNdlHP8e7/qTaOo7sWneoiQWR/YxVm2aPx8WLGDexJ70jg7mtwt3XFCrNtVndOfRbLxjI75mX4a/PZyPkj5yd0jCRSS5u1NpKXz2GVx7LQQarfQPNx8iyN/CZA/vktl2dBvXL7yePpF9WHj9QnzMrl9BylFJaZWSCqvRBfO3v8Hw4XD33fin7OGVmwZitWl+88GWC77/vUrvyN78fPfPDIoexKxFs3hs5WNYbXJsPJ0kd3f6+msoLKzuksk7XcY3SceYPrADgb6eW8iUUZjB1R9dTRv/Nnw962uC/YKb5H3OXqLvbEcKSsBigU8+geBgmDaNOFMZ/7qhHzsyC/nDF7uk3tsuslUkq29dzb2D7uW5Tc8x8aOJnCg54e6wxHmQ5O5Ob78N7dpVr7q0cEsm5VYbs4Z0dHNg5y6nKIerPryK4opivpn1DdFBTfsNpGqJPkeqK2qio2HRImOFq+uvZ3y3MB6+Mp7Pt2by5oYDTRqfJ/E1+/LyxJeZP3k+6w6uo/9r/fnh8A/uDkucI0nu7nLgACxbBnfdBRYLFVYb7286xODOYXSL8sza9rziPK58/0r25+/ny5lf0juyd7O8r6PumQAfM4+N7/7LhqFD4Y03YO1auPtu7h/dhYl92/GPr5P5JknmtavproF3sfGOjfiYfRj59kieXP+knPDkgSS5u8trr4HJBPfcA8DSHUfIKijhnpEXuTmwc5Nfks/Y98eSmpvK0huXcnnc5c323mef8RsTGsCT1/b95SSxKrfdZiys/d57mP7vj/z7hv4M7NiGBz/Zzk8HpAuipoToBLbO2cr0XtP5/drfM+GDCRw5dcTdYYlGUC2hzzEhIUEnJia6O4wmtXhbVvVkVwHWCja+chs7L+pHwfsLmNIvmvEvfI9JKZY/NMLjTlwqKC1g7Ptj2Xl8J1/O/JIJXSe4O6S6aW38Q50/H554ghOPzuO6137geGEp7905mEs6hbk7whZFa838rfN5cPmDBFgC+M9V/+Gmvjd53O+ot1JKbdFaOzx5RFruzaDm4hwAV+/5jrCSk7zedwLzFiXxj6+TSc8+za9HdfG4P5rc4lzGvT+OHcd28PkNn7fsxA7G8oWvvQazZ8MTTxD2/DMsuHsokcH+3PbWz2w5JGdr1qSUYs4lc9jxqx10j+jOzV/czPRPp5NdlN3wk4VbSXJvBjUnuTLZrPz6x89IbhvHpo4XU1Jh5b1NB4kNC2DSxZ61jviB/AMMe2sYSdlJfH7D50yKn+TukJxjMsH//ge33AL/939EvfI8C+4eSkRrX2576ye2ZUiCP1t8eDwbbt/A01c+zbL0ZfT+b2/e2f6OVBu1YJLcm0HNSa4mpWygy4lMXrpsZvUi2BVWzSNj47GYPefHsePYDi576zJyinJYfctqJnef7O6QGsdsNqqVZs2C3/+eds8/xYK7hxDe2pdb3/yJnw9KH/zZzCYzjw97nK1zthIfHs/tX97O5e9czu7s3e4OTTjgOdnEg1WV5JlsVu7/4WNSIzqyvPtl1Y/7mJVHTe27Zv8aRr4zEovJwvrb1zOs4zB3h3RuzGZjioLZs+Evf6H9A79iwa0DaBvkx6z5m41Fy0UtvSN7s/729cyfPJ/dObvp/3p/Hl/1OIWlhe4OTdQgyb0ZVJXqTUrZQLe8w/znspnoGkvL3T3iIkymlt/XrrXmuR+eY9wH44gNjmXTnZuardyxyVgs8NZbxpms779P9A1T+WJGdwZ0DOWhT7bzwuo06XpwwKRM3DXwLlLuTeGWi2/h2R+epet/uvLKT69QYa1wd3gCqZZpNkt/SGfQxOHk+bVm2uwXqFAmFNA7OpivHhjh7vAadKrsFHcuuZOFexYyved03p76NkF+nlmPX6cFC4xWfKdOlH+2iLkplSzamsU19lJL/zqmOhDGQiyPrnqUdQfX0T28O0+OeZJpPaZ5XIGAp5FqmRZg8jfv0a4wm96fv0vaU5MY0S2CVn4W5t/mmilwm1JyTjJD3xzK58mf88yVz7Dw+oXel9jBWOBjzRooKMB36GD+VfATv72yG19sy2LKyxtIPXbK3RG2WJdEX8LaW9ey9MalKKW49tNrGfjGQBanLJZvPm4iyb057NsHzz5rDN4NH87CxEzWp+fyuwndaR9S98yG7ma1Wfn3pn8z4PUBZBdls/LmlTw27DHvbo0NHw47dsCwYah77uH+1//Ah9PjOVFUwZSXN/D+j4ckWdVBKcWk+Ekk/TqJ96a9R1F5Edd8cg0D3xjIZ3s+k8nImpl0yzS1ykoYNQqSkmDPHtIswUx5eQMDYtvw4V1DWmxf+/78/cxePJv1GeuZ0n0Kb0x6g6jWUe4Oq/nYbPDcc/CHP0B0NIX/eokHTrbnu7QcxvaK4m9T+9AuxN/dUbZolbZKFiQt4G/f/430E+l0Du3Mw0Mf5vYBt9Pat7W7w/MK0i3jTk8/DRs3wquvUhQRxW8+3EprPx9enNm/RSb2cms5z/3wHBe/ejE7ju/gnanvsHjG4gsrsYNRC//448bPrlUrQq6fxjvLnuGfg0L5Pi2HMf9ax5sbDlBplTlX6mIxWbil3y0k35vMohsW0T6oPQ8sf4DY52P57YrfkpaX5u4QvZq03JvI4m1ZLH/9c/4z/xG+6TGcByY/hr/FRFmljQ/uGsKwrhHuDrGWb9K/4aEVD5GWl8ak+En8d+J/iQ2JdXdY7ldeDv/6l1FRYzaT/+g8Hm8/glX7T9KzfTB/m9qbhDiZtsAZmw5v4vkfn+eLlC+otFUyOm4091xyD1N7TMXfIt+EGqu+lrsk9yaweFsWr7y5ko/feohC/9ZMu/XfnPQ3vob6mBTPXt+v9qRWbrT92Hb+uPaPLEtfRnx4PC+Mf4Grul3l7rBangMH4P77YdkydIcO7LrjQX5l6UtWUSVjekTyyLh4ekeHuDtKj3Ds9DHe3vY2b2x9g4MFBwn2C2Z6z+nc1PcmRseNxmySyiRnSHJvZhP+tJiX/vsgkUUnuOaWf3Eg7MxEHhMawMa5V7gpul9sP7adv373V75I+YIQvxD+OPKPPDDkAXzNvu4OrWVbu9boi//xR2wXdWHt9XN43KcXJyoUV/dtzwNjutG9nRdWEzUBq83K2gNr+WjXR3y+53NOlZ+iXet2zOw9k5l9ZjIoZhAmJb3HdZHk3pxyc9ndewhd8w4z+/q/sKnTxbV2UcCBp65u/tgwTkT69uC3vLj5RZakLiHEL4SHhz7Mg0MfJNQ/1C0xeSStjfn4//hH2LEDW1Q7fpgwg7nhQ8n0CWJEtwjuGNaZy+PbtsixlZaopKKEZenL+CjpI5alL6PcWk5Uqyiu7nY1k+InMbbLWBmIPYsk9+Zy+DBcfTVlyancfc0f+P6iSxzu5o6W++ny07y/431e/vll9uTsITwgnPsG38dDQx+SpH4+tIZVq+D552H5crSfH6kjr+I/sZfydURPOrcN4sbBHZnaP5rIYOlTdlZ+ST7L0pexNG0py/cu52TZSXzNvoyOG82ErhMYHTeavlF9L/hWvST35vDDD8ZC18XFvPCbJ3mBOIe7BfiYHS8k0QQqbZWs3r+aD3Z+wBcpX1BcUczA9gO5f/D9zOwzUwawXC05GV56CT76CE6epLhdNF/3G8P/ogeTFhnHiPhIrh0Yw5ieUbT289w1cptbhbWCDRkb+CrtK5amLSX9RDoAYQFhjOw0ktFxoxkVN4rebXtfcH31ktybUkUFPPkk/O1v2Dp24ok7/8l7J1s53DUmNIDHxndv0sReWlnKtwe+ZWnaUj5P/pzsomxC/UO5odcNzO4/m6Edhnr3SUgtQUkJLFkC774LK1aAzUZBuw4s7zKELzomkNSxN0O6RzG2Vzuu7BkpLfpGyijM4LuD37Hu4DrWHVrH/vz9AAT5BpEQncDgmMHVl5igGK/+fZfk3lS++w4eegi2b+frvqOZe8U91VUxZ2uqrhitNal5qXx74FtW7l/Jqn2rKKooItAnkKu6XsXNF9/MVV2vws/i5/L3Fk44fhyWLoXFi9GrVqHKyykLCGRLp76sie7LD50uRvXpw9BuUQzrGs7gzmEE+fu4O2qPcqjgEN8f+p7NWZv5Kesnth/bToXNmLysXet29IvqR9/IvvSN6kvfyL70bNvTa761SnJ3tY0bjdb6smXkh7fj9yNu55vu9U9766pB1NPlp9l+bDuJRxLZnLWZdQfXcez0MQBig2OZFD+JyfGTGd15tNf8AnuNU6eM/vnVq9Fr1qDSjJN4Sv0C2BHVlW3t49nRPp7CiwcQ07sb/Tq2oX9sKN3bBeHjQXP9u1tZZRk7ju/gp6yf+PnIzyQdT2JPzh7KrGUAmJWZLmFd6BbWzbiEG9ddw7rSMaSjR3XtSHJ3hfx8tj/3Br7vvU2vzFROBATzv0HTeDNhKmU+DbeKG9tyt2kbhwsPk5KbQkpuCluPbSXxSCLJOclojJ9ZTFAMl8ddzqhOoxjdeTRd2njeMn0XtMOHjW9/mzdj+/FH2LEDU4XR4sxrFUpqeEfSI2I50LYTZfE98OnTm3bdYukWFUzXyNZ0DAvELJU4Tqm0VZKel05SdhJJx5NIzk0m/UQ6e0/spbiiuHo/X7MvF7W5iLjQOGKDY41LyC/XHYI7EOgT6MZPcqYLOrnXXJjarBRWrWtdK6DqKJgU2DT42yrpeXQvgzN3M+LAVoYc3oWPzUpKRCc+6j+BhX3HUuLrXMu4rkHU0+WnySjMOOOy98ReUnJTSMtLo6TylxWcIltFMih6EAnRCSREJ3BJ+0toH+RZy/KJBpSVwfbt8NNP6K1bKd+1B1NKMj6nf5mN8qRvIFkhkRwObceR0ChOR3egMrYTPrEdaN25A2GdY4luG0JkkB9tg/xkmuIGaK05evoo6XnppJ9Ir77OKMzg8MnDDteKDQsIo13rdrQNbEtkq8g6LxGBEYT4hTTpN4FmT+5KqQnAi4AZ+J/W+qn69j+X5N7YpF1nrNpGRFEhXU4cpkteJl3zDtMz+wD9j6bhX1kOwN6wDqyMH8o38cNIate1enm8mjSV2CjBporQFGNTxVgpBFMhV/YJIDy4lOyibLKLsskpzuHoqaPkl565VqdZmekY0pEeET3oEdGD7uHdjeuI7kS1ipJW+YVIazhyBPbsgT17KEtNpzRtL/rgQQKzMvAtLan1lBMBweS0CiWnVRsKgsMoCQ3HFtoG1SYUS3gYlvAwfMLa4NM2DP+IcAKiImgdGkxwoA/B/j4E+prld82utLKUrJNZHD55mMOFhzl88jCZJzOr/5arLmf/LdcU7BdMqH/oGZc2/m2qb0+On8wl0Y7LphvSrMldKWUG0oCxQCbwM3Cj1npPXc9pbHJfvC2LeYuSKCkvxbeyGL/KUnytpfhaS/C1luJTWYavrRSLtYyAiiJalRcRUHmawPIiAiqKCagspnVZIcGlJwkqOwnKRqkFSixQ4G8hKziU/WERHAptw5HgEIp9FJoKNGVoVY6mFJsqxkYxuvq6rN6YwwLCzvhPH9Uqio4hHc+4tA9qj8UkJXLCSVpDbi4cPAjHjlGWkcmpQ5mUHT6CPnYMc/Zx/PNyCCzIw6+s9j+BmspNFop8Ayjx8aPUx49yX3/K/AKo8POnws+fSv8AKv0DsPkHYAsIRAcEoPz8MPv5YvLzxeTnBz4+mPx8wdfYpnx9Mfn+8rjZzw+Tvy9mP1/MPhbMFgs+PmbMPj6YzWawmDBbzCiLBWUyYbZYUBYLJosZk8WM2WLGZLGgTAqzyYRJqRZzgli5tZzc4twzEn5ucS6FpYXkl+ZTUFpQ65Jfms/JspPMnzyfuwbedU7vW19yb4pMMhjYq7Xeb3/zj4GpQJ3JvbGeXZFKSYWVnkf/zvIuW1z1snaVQK5x0WYUvsZF+1bfNmk/TDoYi26HiUBMOhBlvzbRCpMOpG3rUO66rD8zL+lDRGAEPmapgBAuphS0bWtcAD/7xaGKCigogIICKnLyKMnJpTTnBOW5eVTk5WM7cQLryVPoomJ0cTEUFxNQWkJQaQk+BQX4lJXiU1aKb3kpfuVl+Lp5KT2rMmFVigplwqZMWE2/DDhrpdCo6m/Y2r6N6uuq26BR1Y+BQttvVm8/6zUUZz7HeJ9fXj8CRbhS9LQ/B1Xz9WtSQBusKpRs/1IY6KIDU0NTJPcY4HCN+5nAkLN3UkrNAeYAdOzYsVFvcKTAaIUU+Q1jwj6wKR+syse4NvkaF+VDpcmPcksAJT6BlPi0otwcAMqMwgyYUdp+jeMkbuznvOY8QUmIRvHxqf5H4NOtGz5A8Pm8ntVq/MMoL6eypJTykjIqS8uoLCnDWl6GtdS4byurwFpairWsDFt5BbbSMmzl5dgqKrFZrVjt11htaKvVeF2rFW2z2W/bwGZFW+33bTZUzfvahqq0GvtoQGu0thnfasAYQKtSc7s27quqh40n2zOyRmltbLPfr3pe9faq16jat+o1atxWNZ5X/RZgbK+hVWTk+fwk6uS2PgCt9RvAG2B0yzTmudGhAWQVlJARPg4Y59RzzECA/efnClUDr1X9+81xgpIQLYbZbFz8/bEEB7svkYg6NcXPJAuoOQl4B/s2l3lsfHejz72icct2nZ2Q66uWqWuQVpK4EMITNEVy/xnoppTqjJHUZwI3ufINqhJrQ9UykpCFEBcqlyd3rXWlUuo+YAVGb8hbWuvdrn6faQNiJFkLIUQdmqSrTGv9NfB1U7y2EEKIhsmEFUII4YUkuQshhBeS5C6EEF5IkrsQQnihFjErpFIqBzh0jk+PwJgvoKWRuBpH4mq8lhqbxNU45xNXJ611W0cPtIjkfj6UUol1TZzjThJX40hcjddSY5O4Gqep4pJuGSGE8EKS3IUQwgt5Q3J/w90B1EHiahyJq/FaamwSV+M0SVwe3+cuhBCiNm9ouQshhDiLJHchhPBCHp3clVITlFKpSqm9Sqm5bowjVin1rVJqj1Jqt1LqQfv2MKXUKqVUuv26jZviMyultimlvrLf76yU2mw/bp8opXzdEFOoUuozpVSKUipZKXVpSzheSqmH7T/DXUqpBUopf3ccL6XUW0qpbKXUrhrbHB4fZXjJHt9OpVQTLNpWb1zP2n+OO5VSXyilQms8Ns8eV6pSanxzxlXjsd8qpbRSKsJ+363Hy779fvsx262UeqbGdtcdL621R14wphPeB1wE+AI7gF5uiqU9MNB+OwhjgfBewDPAXPv2ucDTborvEeAj4Cv7/U+BmfbbrwG/dkNM7wJ32W/7AqHuPl4YS0QeAAJqHKfZ7jhewEiMlTV31djm8PgAE4FvMBbmHApsbua4xgEW++2na8TVy/536Qd0tv+9mpsrLvv2WIzpxw8BES3keI0GVgN+9vuRTXG8muWPpokO2qXAihr35wHz3B2XPZYvgbFAKtDevq09kOqGWDoAa4ArgK/sv9C5Nf4YzziOzRRTiD2JqrO2u/V48cv6v2EY02F/BYx31/EC4s5KCg6PD/A6cKOj/ZojrrMeuwb40H77jL9Je5K9tDnjAj4D+gEHayR3tx4vjMbClQ72c+nx8uRuGUcLcbt99Q6lVBwwANgMRGmtj9ofOgZEuSGkF4DHAZv9fjhQoLWutN93x3HrDOQAb9u7i/6nlGqFm4+X1joLeA7IAI4ChcAW3H+8qtR1fFrS38IdGK1icHNcSqmpQJbWesdZD7n7eMUDI+xdfd8ppQY1RVyenNxbHKVUa+Bz4CGt9cmaj2njX3Gz1p0qpSYB2VrrLc35vk6wYHxVfVVrPQAowuhmqOam49UGmIrxzycaaAVMaM4YnOWO49MQpdQfgErgwxYQSyDwe+BP7o7FAQvGt8OhwGPAp0op5eo38eTk3uQLcTeGUsoHI7F/qLVeZN98XCnV3v54eyC7mcMaBkxRSh0EPsbomnkRCFVKVa3C5Y7jlglkaq032+9/hpHs3X28rgQOaK1ztNYVwCKMY+ju41WlruPj9r8FpdRsYBIwy/6Px91xdcH4J73D/vvfAdiqlGrn5rjA+P1fpA0/YXyrjnB1XJ6c3KsX4rZXL8wElrgjEPt/3TeBZK31v2s8tAS4zX77Noy++GajtZ6nte6gtY7DOD5rtdazgG+B69wY1zHgsFKqu33TGGAPbj5eGN0xQ5VSgfafaVVcbj1eNdR1fJYAt9qrQIYChTW6b5qcUmoCRtffFK118VnxzlRK+SmlOgPdgJ+aIyatdZLWOlJrHWf//c/EKHo4hpuPF7AYY1AVpVQ8RkFBLq4+Xk01iNAcF4xR7zSMUeU/uDGO4RhfkXcC2+2XiRj922uAdIzR8TA3xjiKX6plLrL/0uwFFmIftW/mePoDifZjthho0xKOF/AXIAXYBbyPUbnQ7McLWIDR71+BkZjurOv4YAySv2L/O0gCEpo5rr0YfcVVv/uv1dj/D/a4UoGrmjOusx4/yC8Dqu4+Xr7AB/bfsa3AFU1xvGT6ASGE8EKe3C0jhBCiDpLchRDCC0lyF0IILyTJXQghvJAkdyGE8EKS3IUQwgtJchdCCC/0/xFbZVWjd46sAAAAAElFTkSuQmCC\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, tf.squeeze(model_x(prior, horizon=days)))\n",
    "ax.plot(t, tf.squeeze(model_x(trainable_weights, horizon=days)), color='r')\n",
    "ax.plot(t, tf.squeeze(raw), color='g')\n",
    "ax.scatter(t[:70], tf.squeeze(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4: What is the number of cases for day no. 80?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SMAoZ4BmhVer"
   },
   "source": [
    "### Task 5: Perform epidemic prediction, find params that won't exceed health system capacity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 283
    },
    "id": "SA1MqkNlhVet",
    "outputId": "1ea57282-96cf-47d9-8569-8c2741b033e8",
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x16ac9c940>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjA0lEQVR4nO3deZgU1b3/8feXZRQEBWWRCDjgjqhARrYRNeCCSESFuMS4JHiJSxaj+akkek2uZjFq1Ny44RK3/KIEN1SMCogbgoyAsjPjwqYoKOACynbuH6damn2G6e5TVf15PU8/1V1V0/2lmPlMzalT55hzDhERSZc6oQsQEZHcU7iLiKSQwl1EJIUU7iIiKaRwFxFJoXqhCwBo1qyZKy0tDV2GiEiivPXWW0udc823tC0W4V5aWkpFRUXoMkREEsXM5m1tm5plRERSSOEuIpJCCncRkRRSuIuIpJDCXUQkhRTuIiIppHAXEUkhhbvA8uVw990wZ07oSkQkRxTuxe6RR+DAA2HIEDj4YPiv/4KvvgpdlYjUksK9mI0fD2eeCXvvDaNHw89+BvfeC1dcEboyEamlWAw/IAE4B7/+NbRqBWPHwi67QJ8+fv3f/gannw69eoWuUkR2kM7ci9Vjj8Ebb8C11/pgz/jDH6C0FAYPhq+/DlaeiNSOwr0YrV8Pv/kNdOwI55238bZGjeD226Gy0rfHi0giKdyL0cSJPrwvvxzq1t18e9++/uLqbbcVvjYRyQmFezEaMQJKSuCkk7a83QwuuggqKuDNNwtbm4jkhMK92Djnw/2442C33ba+39lnQ+PG8Pe/F642EckZhXuxmTQJ5s+HQYO2vV/jxnDOOfDoo/DZZ4WpTURyRuFebEaMgPr1t94kk+3cc2H1anj66fzXJSI5pXAvNk884fuzN226/X3LyqBtW99tUkQSReFeTBYtgqoq395eHWZw6qnwwgvwxRf5rU1EckrhXkxefdUva3Ln6cCB8M038Oyz+alJRPJC4V5MXn3V36TUqVP1v6ZnT9hzTzXNiCSMwr2YvPYa9OgB9WowpFCdOnDyyTBqlIYjEEmQaoW7mf3KzGaY2XQz+5eZ7Wxm7cxsoplVmdmjZlYS7btT9Loq2l6a13+BVM/y5TBt2o4NBta/P6xc6X85iEgibDfczWwv4BdAmXOuI1AXOAO4HrjZObcvsAwYHH3JYGBZtP7maD8J7fXX/Q1MRxxR86896ijfffKFF3Jfl4jkRXWbZeoBDcysHtAQ+AjoDYyItj8AnBw9HxC9Jtrex8wsJ9XKjnv1VR/Q3brV/GsbNfK/FJ5/Pvd1iUhebDfcnXOLgBuB+fhQXwG8BSx3zq2NdlsI7BU93wtYEH3t2mj/PTZ9XzMbYmYVZlaxZMmS2v47ZHteew26dIGGDXfs6487Dt55Bz76KLd1iUheVKdZpin+bLwd8B1gF6BvbT/YOTfMOVfmnCtr3rx5bd9OtmXdOpg6Fbp23fH3OP54v3zxxZyUJCL5VZ1mmWOA951zS5xza4DHgXKgSdRMA9AaWBQ9XwS0AYi27wZ8mtOqpWaqqvy8qJ077/h7HHYYNG+uphmRhKhOuM8HuptZw6jtvA8wE3gJyIw+dS7wVPR8ZPSaaPtY55zLXclSY1Om+GVtwr1OHd808+KL/sKsiMRaddrcJ+IvjE4GpkVfMwy4ArjUzKrwber3Rl9yL7BHtP5S4Mo81C01MXWqv5jaoUPt3qd3b1iyBGbNyklZIpI/1bqbxTl3DXDNJqvfAzZrxHXOfQ38oPalSc5MmeJnViopqd37HHWUX44bV/tfFCKSV7pDNe2c8+FemyaZjPbtoXVrH+4iEmsK97T78EPflJKLcDeDo4+Gl19Wu7tIzCnc027qVL/MRbiDb5r55BOYPTs37ycieaFwT7tMT5lDD83N+x19tF+qaUYk1hTuaTd1Kuy7L+y6a27eb599YK+9fNOMiMSWwj3tZsyAQw7J3fuZ+aaZV15Ru7tIjCnc02z1aqiszH23xfJyP8bMBx/k9n1FJGcU7mlWWenHlclHuIMfRlhEYknhnmYzZ/plrsO9Y0do3FjhLhJjCvc0mznTt5EfcEBu37duXejeHcaPz+37ikjOKNzTbOZMf1dpgwa5f+/ycj9t34oVuX9vEak1hXuazZyZvzFgyst9b5kJE/Lz/iJSKwr3tFq7FubMyV+4d+vmhwFWu7tILCnc0+rdd2HNmvyFe+PG/q5XtbuLxJLCPa3y1VMmW3m5b5ZZu3b7+4pIQSnc0yoT7gcemL/PKC/30/e9807+PkNEdojCPa1mzYK2baFRo/x9Rs+efqmmGZHYUbin1dy5ue/fvqm2bf0gYrqoKhI7Cvc0cs6H+/775/dzzHzTjMJdJHYU7mm0ZIm/uSjf4Q6+aWbBAv8QkdhQuKfR3Ll+WYhwzwwipnZ3kVhRuKdRIcP9sMOgYUM1zYjEjMI9jebOhfr1Ye+98/9Z9etD1646cxeJGYV7Gs2d66fWq1u3MJ/Xs6efzu+rrwrzeSKyXQr3NJozpzBNMhk9evhJQSoqCveZIrJNCve0WbcOqqoKG+7du/vlG28U7jNFZJsU7mkzf76fO7WQ4d6smf88tbuLxIbCPW0yPWXyfXfqpnr08GfuzhX2c0VkixTuaVPIbpDZevaEpUv9UMMiEpzCPW0qK/1Y6y1aFPZze/TwSzXNiMSCwj1tqqp8N0izwn5uhw6w6666qCoSEwr3tMmEe6HVreun3tOZu0gsKNzTZO1aeP/9MOEOvmlm+nT4/PMwny8i31K4p8n8+T7gQ4V7z56wfj28+WaYzxeRbync06Sqyi9DhXu3bn6pdneR4BTuaRI63Js0gYMPVriLxEC1wt3MmpjZCDObbWazzKyHme1uZi+aWWW0bBrta2b2NzOrMrN3zKxLfv8J8q2qKmjQAFq1CldD5mam9evD1SAi1T5zvxX4j3PuQOAwYBZwJTDGObcfMCZ6DXACsF/0GALckdOKZeuqqmCffQrfDTJbjx6wfLkfvExEgqm3vR3MbDfgSOA8AOfcamC1mQ0Ajo52ewAYB1wBDAAedM45YEJ01t/KOfdRzqsHjj766M3WnXbaaVx00UWsXLmSfv36bbb9vPPO47zzzmPp0qUMGjRos+0XXnghp59+OgsWLODss8/ebPtll13G97//febMmcNPf/rTzbZfddVVHHPMMUydOpVLLrlks+1//OMf6dmzJ+PHj+c3v/nNZttvueUWOnXqxOjRo7nuuus2237XXXdxwAEH8PTTT3PTTTdt2DBpEjRowEMLFtCmTRseffRR7rhj89+tI0aMoFmzZtx///3cf//9m20fNWoUDRs25Pbbb2f48OGbbR83bhwAN954I88888xG2xqsX89zAOPHc+2IEYwZM2aj7XvssQePPfYYAEOHDuWNTZpwWrduzcMPPwzAJZdcwtSpUzfavv/++zNs2DAAhgwZwtzMHbmRTp06ccsttwDwox/9iIULF260vUePHvzpT38CYODAgXz66acbbe/Tpw9XX301ACeccAKrVq3aaHv//v359a9/Deh7b6PvvchDDz0U7nuvQQOee+45AK699trEfO9l/k25Vp0z93bAEuAfZjbFzO4xs12AllmBvRhoGT3fC8ieUHNhtG4jZjbEzCrMrGLJkiU7/i+QDVat8s0yITVsCE2bqt1dJDBz2xnoyczKgAlAuXNuopndCnwO/Nw51yRrv2XOuaZm9gzwZ+fca9H6McAVzrmtDvZdVlbmKjQWeO0sWABt28Kdd8IWzugK6sQT4YMPYMaMsHWIpJyZveWcK9vStuqcuS8EFjrnJkavRwBdgI/NrFX0Aa2AT6Lti4A2WV/fOlon+RS6p0y2Hj1g5kxYtix0JSJFa7vh7pxbDCwws8wYsn2AmcBI4Nxo3bnAU9HzkcA5Ua+Z7sCKfLW3S5a4hTvAxInb3k9E8ma7F1QjPwf+aWYlwHvAj/G/GIab2WBgHnBatO8ooB9QBayM9pV8q6qCkhJo3Tp0JX7C7Dp1/DgzffuGrkakKFUr3J1zU4Ettev02cK+Dri4dmVJjVVVQfv2hZsUe1saN4ZDDtFFVZGAdIdqWmT6uMdFz56+WWbdutCViBQlhXsaOBduqN+t6dEDvvhCPWZEAlG4p8HixbByZbzCvWdPv9T47iJBKNzTIDNvaZzCvX17aNkSXnstdCUiRUnhngZx6gaZYQa9esGrr4auRKQoKdzToKrK95LZe+/QlWysVy8/gcj8+aErESk6Cvc0qKqC0lKoXz90JRvr1csvdfYuUnAK9zSIW0+ZjEMP9X3e1e4uUnAK96TLdIOMUx/3jLp1fa8ZnbmLFJzCPek+/RRWrIjnmTv4ppkZM3ydIlIwCveki2NPmWyZdvfXXw9bh0iRUbgnXSbc49gsA34QsZIStbuLFJjCPekqK32f8vbtQ1eyZTvvDGVlancXKTCFe9JVVvr+7TvvHLqSrevVCyoq/BAJIlIQCvekq6yE/fYLXcW29eoFa9dq8g6RAlK4J5lzyQj38nLfdKSmGZGCUbgn2dKlvhtk3MO9SRM/eYcuqooUjMI9ySor/TLu4Q5wxBF+Zqa1a0NXIlIUFO5JlqRw79ULvvwSpkwJXYlIUVC4J1llpb/Fv1270JVs31FH+eW4cUHLECkWCvckq6yM52iQW9KqFXToAGPGhK5EpCgo3JMsCT1lsvXu7XvMrF4duhKR1FO4J1VSukFm693b38j05puhKxFJPYV7Ui1e7C9QJincjzrK93cfOzZ0JSKpp3BPqiT1lMnYfXfo3Fnt7iIFoHBPqky4779/2Dpqqk8f399d48yI5JXCPakqK30vmbZtQ1dSM717w5o1Gt9dJM8U7klVWemH+a1XL3QlNXPEEb5mtbuL5JXCPamS1lMmo1Ej6NZN4S6SZwr3JFq/3s/AlMRwB980U1EBy5eHrkQktRTuSfThh7BqVXLDvU8f/wvqlVdCVyKSWgr3JEpiN8hs3bv7maPUNCOSNwr3JEp6uO+0k7+wqnAXyRuFexJVVvqAbNMmdCU7rndvmDbN32krIjmncE+iykrYZx+ok+D/vr59/fL558PWIZJSCU6HIpbUbpDZOnXywwCPGhW6EpFUqna4m1ldM5tiZs9Er9uZ2UQzqzKzR82sJFq/U/S6Ktpemqfai9P69fDuu8kbdmBTZv7s/YUXNPWeSB7U5Mz9l8CsrNfXAzc75/YFlgGDo/WDgWXR+puj/SRXFiyAb75J/pk7QL9+vq/7hAmhKxFJnWqFu5m1Bk4E7oleG9AbGBHt8gBwcvR8QPSaaHufaH/JhTlz/DIN4X7ssX6aQDXNiORcdc/cbwEuB9ZHr/cAljvnMn9PLwT2ip7vBSwAiLaviPbfiJkNMbMKM6tYsmTJjlVfjGbP9suDDgpbRy7stpvvEqlwF8m57Ya7mfUHPnHOvZXLD3bODXPOlTnnypo3b57Lt063WbOgaVNo0SJ0Jblxwgnw9tuwaFHoSkRSpTpn7uXASWb2AfAIvjnmVqCJmWWGJGwNZH46FwFtAKLtuwGf5rDm4jZrlj9rT0tLV79+fvmf/4StQyRlthvuzrmhzrnWzrlS4AxgrHPuLOAlYFC027nAU9HzkdFrou1jnXMup1UXs9mz4cADQ1eROx07QuvWapoRybHa9HO/ArjUzKrwber3RuvvBfaI1l8KXFm7EuVby5bBxx+no709w8w3zbz4IqxeHboakdSoUbg758Y55/pHz99zznV1zu3rnPuBc+6baP3X0et9o+3v5aPwojQr6omapjN38E0zX3wB48eHrkQkNXSHapKkqadMtj59/JSBapoRyRmFe5LMmuUHDCstDV1JbjVuDEceCc8+G7oSkdRQuCfJrFl+2IG6dUNXknv9+8PMmX6GKRGpNYV7ksyenb4mmYyTT/bLJ58MWYVIaijck+Lrr+H999N3MTWjtNSPFPnEE6ErEUkFhXtSzJnjR4RM65k7wCmnwBtvaAIPkRxQuCfF9Ol+ecghYevIp1NOAedg5MjQlYgknsI9KaZN890Fkz6O+7Z07Ajt26tpRiQHFO5JMX26b2+vXz90JfljBgMHwujR/m5cEdlhCvekmDYt3U0yGaed5mdmUq8ZkVpRuCfB55/D/Pm+2SLtvvtdaNcOhg8PXYlIoinck2DGDL8shnA382fvo0fDpxopWmRHKdyTYNo0vyyGZhlQ04xIDijck2D6dGjUCNq2DV1JYXTuDPvsA488EroSkcRSuCfB9Om+SaZOkfx3mcEPfwhjx8KHH4auRiSRiiQtEsw53yxTDO3t2c46y9+Rq7N3kR2icI+7xYth6VI49NDQlRTWAQfA4YfDww+HrkQkkRTucTd5sl927hy2jhB+9COYMmVDbyERqTaFe9xNmeLboA87LHQlhXf66X7s+oceCl2JSOIo3ONu8mTYbz8/W1GxadkS+vb14b52behqRBJF4R53kycXZ5NMxuDBvsfM88+HrkQkURTucfbZZzBvHnTpErqScE48EZo3h/vuC12JSKIo3ONsyhS/LOYz95ISOOccP8b7J5+ErkYkMRTucaZw937yE9/mrgurItWmcI+zyZOhTRto1ix0JWF16ADl5XDXXf7GJhHZLoV7nE2ZUtzt7dkuvBAqK2HMmNCViCSCwj2uVqzwk2Ir3L1Bg/xfMLffHroSkURQuMfVpEl+XJnu3UNXEg877QTnn+8vrC5YELoakdhTuMfVhAl+2bVr2Dri5Kc/9b/w7rwzdCUisadwj6uJE+Ggg6BJk9CVxEdpKZx8sg/3lStDVyMSawr3OHLOn7mrSWZzl17qb+568MHQlYjEmsI9jt57zw/zq3DfXHm5Hwr45pvVLVJkGxTucZRpb1e4b87Mn73PnQvPPBO6GpHYUrjH0YQJsMsucPDBoSuJp4EDffv7n/7km7BEZDMK9ziaMME3PdStG7qSeKpfHy6/3B+nceNCVyMSSwr3uPn8c39nanl56Eri7cc/hj33hD/8IXQlIrG03XA3szZm9pKZzTSzGWb2y2j97mb2oplVRsum0Xozs7+ZWZWZvWNmusWyJl5/Hdatg6OPDl1JvO28M1x2mR+O4I03QlcjEjvVOXNfC1zmnOsAdAcuNrMOwJXAGOfcfsCY6DXACcB+0WMIcEfOq06zl1/2zQ49eoSuJP4uuMCP9f7f/x26EpHY2W64O+c+cs5Njp5/AcwC9gIGAA9Euz0AnBw9HwA86LwJQBMza5XrwlNr3Dh/V+ouu4SuJP4aNYKhQ2H0aLW9i2yiRm3uZlYKdAYmAi2dcx9FmxYDLaPnewHZg38sjNZt+l5DzKzCzCqWLFlS07rT6YsvoKJCTTI1ccEF8J3vwNVXq+eMSJZqh7uZNQIeAy5xzn2evc0554Aa/WQ554Y558qcc2XNmzevyZeml9rba65BA7jqKnjtNXj22dDViMRGtcLdzOrjg/2fzrnHo9UfZ5pbomVmDrRFQJusL28drZPtGTdO7e07YvBg2G8/3z1y7drQ1YjEQnV6yxhwLzDLOffXrE0jgXOj5+cCT2WtPyfqNdMdWJHVfCPbMnas79+u9vaaKSmBv/wFZs2Cu+8OXY1ILFTnzL0cOBvobWZTo0c/4M/AsWZWCRwTvQYYBbwHVAF3AxflvuwU+uQT397et2/oSpJpwAA48ki45ho/0YlIkau3vR2cc68BtpXNfbawvwMurmVdxef55/0FwX79QleSTGZ+MLGyMt818tZbQ1ckEpTuUI2LUaOgRQvo3Dl0JcnVpYufa/Xvf4e33w5djUhQCvc4WLfOn7mfcALU0X9JrVx3Hey+O1x0kYYElqKmJImDiRNh2TI1yeRC06Zwww0wfrym45OipnCPg1Gj/AiQxx4bupJ0OPdcfyyvuALmzQtdjUgQCvc4eOIJPwpk06ahK0kHMxg2zF+gzkyqLVJkFO6hzZgBM2fCaaeFriRdSkvh+uv9tYw7NHadFB+Fe2jDh/uLqAMHhq4kfS66yN83cNll/gYnkSKicA/JOR/uRx3lJ56Q3DKDf/zDjx75wx/C11+HrkikYBTuIU2bBrNnq0kmn/bcE+6/H6ZOhV/+MnQ1IgWjcA/p0Ud9k8ypp4auJN1OPNGP+z5sGDz4YOhqRApC4R7K2rX+jPK44/ydqZJf//M/fijlIUPgzTdDVyOSdwr3UJ57Dj780IeN5F+9ev76RqtWfpCxhQtDVySSVwr3UO6+27cH9+8fupLi0bw5PP00fPmlP+4aPVJSTOEewsKFftagn/zET84hhdOxI4wY4e8vOPVU+Oab0BWJ5IXCPYR77vGDWg0eHLqS4nT88XDffX5ylDPPhNWrQ1ckknMK90L76is/JO33vw/t24eupnidfbYf8/2JJ+AHP9AZvKSOwr3Q7rkHPv3Ud82TsH7xC/+LduRIf4ewbnKSFFG4F9Lq1XDTTdCrlybBjouLL/Zjzzz7LJxyCqxaFboikZzY7jR7kkMPPwwLFmic8bi54ALfVXLIEOjdG558Elq2DF2VSK3ozL1QvvoKrr4aunb1My5JvJx/Pvz73356vsMP1zR9kngK90L561/9TUs33eQHtJL4GTgQXnvN92QqL/dn8CIJpXAvhI8+8mOLn3oqHHFE6GpkW7p0gUmToEMH3wZ/+eXqSSOJpHDPN+fgZz+DNWt8wEv8tWoFL7/sZ3G64Qbo3l3jwUviKNzzbfhwePxxP3DVvvuGrkaqq0EDf+H7qaf8HcVduvhuk+vXh65MpFoU7vm0eLHvate1q58NSJLnpJP8uPvf+x78/Of+LH7SpNBViWyXwj1fVq/2dz6uXOlnA6qnXqeJteeevh98pitrt26+2+TSpaErE9kqhXu+XHqp73lx333+4pwkmxmcdRbMmQO/+pX/f91nH/jd72D58tDViWxG4Z4PN9wAt93mm2LOOCN0NZJLu+7qu7O+8w706QO//z20a+evqWgIYYkRhXuu3Xab7z53+unqHZNmHTr4C+VTpvgJzq+5Blq39u3yc+aErk5E4Z4zzsF11/lujwMGwEMPQd26oauSfOvUyd/sNHmy7xc/bBgceKCfPvGxxzQYmQSjcM+FlSv9xBtXX+2Hkh0+XJNwFJvOnf3k2wsWwLXXwsyZMGiQvxh7/vkwbhysWxe6SikiCvfaevttKCvzk13/7nfwwANQUhK6KgmlRQu46iqYNw9eeMH/Fffoo74r5Z57wo9/7JtzvvwydKWScuacC10DZWVlrqKiInQZNbNihb+Y9r//6+fmfOABOPbY0FVJHK1c6eduHTkSRo3yvWtKSnyf+e99z49E2a0b7LRT6EolYczsLedc2Ra3KdxraMUKf6fizTfDZ5/5P7n/+Edo1ix0ZZIEa9bA66/7fvMvveTb6p3zwd6pk/8r8PDD/eOAA3TdRrZJ4V5b69fDG2/4ppd//csP33viif7M/bvfDV2dJNmyZfDKK/Dqq1BRAW+9taHJplEjOPRQOOigjR977w111KIqCvcds2iRvwg2bhw895x/3bCh77d+8cV+rBGRXFu3znelnDTJP6ZN84OWLVmyYZ8GDXzf+r33htLSjZetWvmJRnbeOdS/QAqo4OFuZn2BW4G6wD3OuT9va/+g4b5yJXzwge/dMGOGf0yeDO++67fvtptvFx00yE9qveuuYeqU4vbppzB7tg/6WbPg/ff99+28eb55cFO77uov4LZsueHRogU0abL1R6NG+osgYQoa7mZWF5gLHAssBCYBZzrnZm7ta3Y43Neu9eG8atXWH59/7r/5M49ly/wPyqJFvtta9g+GGbRvD4ccAkce6W9OOewwtXtKvH3xhQ/5efP8YHUff7xhmf18e8Mk1Knjfyk0auT/Ss1+NGiw+bqGDf1fCCUlvuvvjizr1vWfW6fOhufbWpe9TZPebDPc8zGaVVegyjn3XvThjwADgK2G+w678UYYOrT6+++2G+y+OzRtCm3aQM+eftm2rW/LPPBA/w0rkiSNG0PHjv6xLWvW+JOd5cs3PFas2Pj18uX+mtLKlRs/li7d8HzVqg3LkMy2HPyZ0Dfb/Pm2tuXj+Za2ZdcP/u7m00+v/fHYRD7CfS9gQdbrhUC3TXcysyHAEIC2bdvu2Cf17u3HcWnQYMMjc5aReey6qw/0Jk00MqMUt/r1YY89/CMX1q/3d+CuWeMfq1dve7nputWr/XusX++vNWxpWd112UvwPZAyrRKZ55u+zvfzLW3LyH7dtGnt/h+2IljaOeeGAcPAN8vs0Jt07eofIlJ4deroL90Yy8fVk0VAm6zXraN1IiJSIPkI90nAfmbWzsxKgDOAkXn4HBER2YqcN8s459aa2c+A5/FdIe9zzs3I9eeIiMjW5aXN3Tk3ChiVj/cWEZHt0x0LIiIppHAXEUkhhbuISAop3EVEUigWo0Ka2RJg3g5+eTNgaQ7LyRXVVTOqq+biWpvqqpna1LW3c675ljbEItxrw8wqtjZwTkiqq2ZUV83FtTbVVTP5qkvNMiIiKaRwFxFJoTSE+7DQBWyF6qoZ1VVzca1NddVMXupKfJu7iIhsLg1n7iIisgmFu4hICiU63M2sr5nNMbMqM7syYB1tzOwlM5tpZjPM7JfR+t3N7EUzq4yW+ZlyZfv11TWzKWb2TPS6nZlNjI7bo9HQzIWuqYmZjTCz2WY2y8x6xOF4mdmvov/D6Wb2LzPbOcTxMrP7zOwTM5uetW6Lx8e8v0X1vWNmXQpc1w3R/+M7ZvaEmTXJ2jY0qmuOmR1fyLqytl1mZs7MmkWvgx6vaP3Po2M2w8z+krU+d8fLOZfIB3444XeB9kAJ8DbQIVAtrYAu0fPG+AnCOwB/Aa6M1l8JXB+ovkuB/w88E70eDpwRPb8TuDBATQ8A50fPS4AmoY8XforI94EGWcfpvBDHCzgS6AJMz1q3xeMD9AOeAwzoDkwscF3HAfWi59dn1dUh+rncCWgX/bzWLVRd0fo2+OHH5wHNYnK8vgeMBnaKXrfIx/EqyA9Nng5aD+D5rNdDgaGh64pqeQo4FpgDtIrWtQLmBKilNTAG6A08E31DL836YdzoOBaopt2iELVN1gc9XmyY/3d3/HDYzwDHhzpeQOkmobDF4wPcBZy5pf0KUdcm204B/hk93+hnMgrZHoWsCxgBHAZ8kBXuQY8X/mThmC3sl9PjleRmmS1NxL1XoFq+ZWalQGdgItDSOfdRtGkx0DJASbcAlwPro9d7AMudc2uj1yGOWztgCfCPqLnoHjPbhcDHyzm3CLgRmA98BKwA3iL88crY2vGJ08/CT/BnxRC4LjMbACxyzr29yabQx2t/oFfU1PeymR2ej7qSHO6xY2aNgMeAS5xzn2dvc/5XcUH7nZpZf+AT59xbhfzcaqiH/1P1DudcZ+ArfDPDtwIdr6bAAPwvn+8AuwB9C1lDdYU4PttjZr8F1gL/jEEtDYHfAP8dupYtqIf/67A78P+A4WZmuf6QJId7rCbiNrP6+GD/p3Pu8Wj1x2bWKtreCvikwGWVAyeZ2QfAI/immVuBJmaWmYUrxHFbCCx0zk2MXo/Ah33o43UM8L5zbolzbg3wOP4Yhj5eGVs7PsF/FszsPKA/cFb0iyd0Xfvgf0m/HX3/twYmm9megesC//3/uPPexP9V3SzXdSU53GMzEXf0W/deYJZz7q9Zm0YC50bPz8W3xReMc26oc661c64Uf3zGOufOAl4CBgWsazGwwMwOiFb1AWYS+Hjhm2O6m1nD6P80U1fQ45Vla8dnJHBO1AukO7Aiq/km78ysL77p7yTn3MpN6j3DzHYys3bAfsCbhajJOTfNOdfCOVcaff8vxHd6WEzg4wU8ib+oipntj+9QsJRcH698XUQoxAN/1Xsu/qrybwPWcQT+T+R3gKnRox++fXsMUIm/Or57wBqPZkNvmfbRN00V8G+iq/YFrqcTUBEdsyeBpnE4XsDvgdnAdOAhfM+Fgh8v4F/4dv81+GAavLXjg79Iflv0czANKCtwXVX4tuLM9/6dWfv/NqprDnBCIevaZPsHbLigGvp4lQAPR99jk4He+TheGn5ARCSFktwsIyIiW6FwFxFJIYW7iEgKKdxFRFJI4S4ikkIKdxGRFFK4i4ik0P8Bjcx+4j7QDMsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "health_sys_cap = 600\n",
    "ax.plot(t, tf.squeeze(model_x(trainable_weights, horizon=days)), color='r')\n",
    "ax.plot(t, [health_sys_cap for _ in t], '--k')"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "name": "Lab-part-2.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "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": 1
}