sim/experiments.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Code setup for Experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import sys\n",
    "if '..' not in sys.path:\n",
    "    sys.path.append('..')\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import copy\n",
    "import scipy as sp\n",
    "import math\n",
    "import seaborn\n",
    "import pickle\n",
    "import warnings\n",
    "import matplotlib\n",
    "import re\n",
    "import multiprocessing\n",
    "\n",
    "from lib.mobilitysim import MobilitySimulator\n",
    "from lib.dynamics import DiseaseModel\n",
    "from lib.inference import * \n",
    "from bayes_opt import BayesianOptimization\n",
    "from lib.parallel import *\n",
    "from lib.distributions import CovidDistributions\n",
    "from lib.plot import Plotter\n",
    "from lib.data import collect_data_from_df\n",
    "from lib.measures import (\n",
    "    MeasureList, \n",
    "    BetaMultiplierMeasure, \n",
    "    BetaMultiplierMeasureByType,\n",
    "    SocialDistancingForAllMeasure, \n",
    "    SocialDistancingForKGroups,\n",
    "    SocialDistancingByAgeMeasure,\n",
    "    SocialDistancingForPositiveMeasure, \n",
    "    ComplianceForAllMeasure,\n",
    "    Interval)\n",
    "\n",
    "from lib.mobilitysim import MobilitySimulator\n",
    "from lib.town_data import generate_population, generate_sites, compute_distances\n",
    "from lib.town_maps import MapIllustrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_settings = 'lib/tu_settings_20_10.pk'\n",
    "random_repeats = 1\n",
    "cpus_used = multiprocessing.cpu_count() \n",
    "c = 0 # seed\n",
    "FIGSIZE = (8, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import town settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import downsampled version of Tübingen (population downsampled by 20x, sites by 10x) that was used used to generate mobility traces on the fly for inference (sites, home locations, etc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# See town-generator.ipynb for an example on how to create own settings\n",
    "with open(mob_settings, 'rb') as fp:\n",
    "    obj = pickle.load(fp)\n",
    "mob = MobilitySimulator(**obj)\n",
    "np.random.seed(c)\n",
    "runstr = f'run{c}_'\n",
    "\n",
    "# General note for plotting: `errorevery` has to be set proportionally to `acc`, and best to keep `acc` as it is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import COVID-19 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "days_until_lockdown = 13 # March 10 - March 23\n",
    "days_after_lockdown = 20 # March 24 - April 12\n",
    "days_present = days_until_lockdown + days_after_lockdown + 2 # simulate 2 more days due to test lag\n",
    "days_future = 12 * 7 # projecting 12 weeks into the future"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_downsample = 10\n",
    "new_cases_ = collect_data_from_df('LK Tübingen', 'new', until=days_present)\n",
    "resistant_cases_ = collect_data_from_df('LK Tübingen', 'recovered', until=days_present)\n",
    "fatality_cases_ = collect_data_from_df('LK Tübingen', 'fatality', until=days_present)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fatality rate per age group\n",
    "num_age_groups = fatality_cases_.shape[1] \n",
    "fatality_rates_by_age = (fatality_cases_[-1, :] / \\\n",
    "    (new_cases_[-1, :] + fatality_cases_[-1, :] + resistant_cases_[-1, :]))\n",
    "\n",
    "print('Empirical fatality rates per age group:  ', fatality_rates_by_age.tolist())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scale down cases based on number of people in simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_cases, resistant_cases, fatality_cases = (\n",
    "    1/case_downsample * new_cases_, \n",
    "    1/case_downsample * resistant_cases_, \n",
    "    1/case_downsample * fatality_cases_)\n",
    "new_cases, resistant_cases, fatality_cases = np.ceil(new_cases), np.ceil(resistant_cases), np.ceil(fatality_cases)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define initial seed count (based on infection counts on March 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_seeds = {\n",
    "    'expo' : 1,\n",
    "    'ipre' : 1,\n",
    "    'isym' : 3,\n",
    "    'iasy' : 3,\n",
    "}\n",
    "present_seeds = {\n",
    "    'expo': 3,\n",
    "    'ipre': 1,\n",
    "    'iasy': 0,\n",
    "    'isym_notposi': 8,\n",
    "    'isym_posi': 4,\n",
    "    'resi_notposi': 78,\n",
    "    'resi_posi': 110,\n",
    "}\n",
    "max_time_present = 24.0 * (days_present)\n",
    "max_time_future = 24.0 * (days_future) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define standard testing parameters, same used for inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standard_testing(max_time):\n",
    "    standard_testing_params = {\n",
    "        'testing_t_window'    : [0.0, max_time], # in hours\n",
    "        'testing_frequency'   : 24.0,     # in hours\n",
    "        'test_reporting_lag'  : 48.0,     # in hours (actual and self-report delay)\n",
    "        'tests_per_batch'     : 10,       # assume 300 tests/day in LK Tübingen\n",
    "        'test_smart_delta'    : 24.0 * 3, # in hours\n",
    "        'test_smart_duration' : 24.0 * 7, # in hours\n",
    "        'test_smart_action'   : 'isolate', \n",
    "        'test_smart_num_contacts'   : 10, \n",
    "        'test_targets'        : 'isym',\n",
    "        'test_queue_policy'   : 'fifo',\n",
    "        'smart_tracing'       : None, \n",
    "    }\n",
    "    return standard_testing_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define distributions as estimated by literature \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantiate correct distributions\n",
    "distributions = CovidDistributions(fatality_rates_by_age=fatality_rates_by_age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set epidemic parameters as inferred using Bayesian optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inferred_params = {\n",
    "    'betas': [1.1383] * 5,  # site infectivity by type\n",
    "    'alpha': 0.3224,\n",
    "    'mu': 0.2072\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define function to run general type of experiment, fixing the above settings.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run(tparam, measure_list, t, local_seeds):\n",
    "\n",
    "    # add standard measure of positives staying isolated\n",
    "    measure_list +=  [\n",
    "        SocialDistancingForPositiveMeasure(\n",
    "            t_window=Interval(0.0, t), p_stay_home=1.0)\n",
    "    ]\n",
    "    measure_list = MeasureList(measure_list)\n",
    "\n",
    "    # run simulations\n",
    "    summary = launch_parallel_simulations(\n",
    "        mob_settings, \n",
    "        distributions, \n",
    "        random_repeats, multiprocessing.cpu_count(), \n",
    "        inferred_params, local_seeds, tparam, measure_list, \n",
    "        max_time=t, \n",
    "        num_people=mob.num_people, \n",
    "        num_sites=mob.num_sites, \n",
    "        site_loc=mob.site_loc, \n",
    "        home_loc=mob.home_loc, verbose=False)\n",
    "    return summary\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_summary(summary, filename):\n",
    "    with open('summaries/' + filename, 'wb') as fp:\n",
    "        pickle.dump(summary, fp)\n",
    "    \n",
    "def load_summary(filename):\n",
    "    with open('summaries/' + filename, 'rb') as fp:\n",
    "        summary = pickle.load(fp)\n",
    "    return summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiments for \"A spatiotemporal model for Sars-CoV-2\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Model inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2.1 Model projections in uncontrolled scenario until present day (April 13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_params_opt = standard_testing(max_time_present)\n",
    "measure_list = []\n",
    "\n",
    "summary_INF = run(testing_params_opt, measure_list, max_time_present, initial_seeds)\n",
    "\n",
    "save_summary(summary_INF, 'summary_INF.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_INF = load_summary('summary_INF.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "plotter.plot_positives_vs_target(\n",
    "    summary_INF, new_cases.sum(axis=1), \n",
    "    test_lag=2,\n",
    "    title=('Predictions for uncontrolled scenario after March 23'\n",
    "           'for inferred parameters'), \n",
    "    filename=runstr + 'opt_00',\n",
    "    figsize=FIGSIZE,\n",
    "    errorevery=1, acc=500, \n",
    "    lockdown_at=days_until_lockdown,\n",
    "    ymax=250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = new_cases.sum(axis=1)\n",
    "targets = np.hstack((targets[0], np.diff(targets)))\n",
    "\n",
    "plotter = Plotter()\n",
    "plotter.plot_daily_infected(\n",
    "    summary_INF,\n",
    "    title='', \n",
    "    filename=runstr + 'opt_001',\n",
    "    figsize=(6, 4),\n",
    "    legend_loc=0,\n",
    "    errorevery=100, acc=1000, \n",
    "    lockdown_at=days_until_lockdown,\n",
    "    lockdown_label_y=500,\n",
    "    start_date='2020-03-10',\n",
    "    show_target=targets,\n",
    "    ymax=1800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "plotter.plot_daily_rts(\n",
    "    summary_INF,\n",
    "    filename=runstr + 'opt_003',\n",
    "    start_date='2020-03-10',\n",
    "    sigma=0.37, # prior found by MLE, set to None to recompute\n",
    "    figsize=(6, 4), \n",
    "    subplot_adjust=None,\n",
    "    lockdown_label='Lockdown', \n",
    "    lockdown_at=days_until_lockdown,\n",
    "    lockdown_label_y=2.6,\n",
    "    titles=None,\n",
    "    ymax=4.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3. Social distancing and Business restrictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3.1. Reproducing the observed course of Sars-CoV-2 under the taken interventions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# {0: 'education', 1: 'social', 2: 'bus_stop', 3: 'office', 4: 'supermarket'}\n",
    "tuned_p_stay_home = 0.2\n",
    "tuned_site_multipliers = [0.0, 0.0, 0.5, 0.5, 0.5]\n",
    "\n",
    "observed_measures_SD = [\n",
    "    # education close, social closes, rest reduced by 50%\n",
    "    BetaMultiplierMeasureByType(\n",
    "        t_window=Interval(24.0 * days_until_lockdown, max_time_present), \n",
    "        beta_multiplier=tuned_site_multipliers),\n",
    "    \n",
    "    # 40 % less activite of all due to contact constraints\n",
    "    SocialDistancingForAllMeasure(\n",
    "     t_window=Interval(24.0 * days_until_lockdown, max_time_present), \n",
    "        p_stay_home=tuned_p_stay_home)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_params_SD_0 = standard_testing(max_time_present)\n",
    "summary_SD_0 = run(testing_params_SD_0, observed_measures_SD, max_time_present, initial_seeds)\n",
    "\n",
    "save_summary(summary_SD_0, 'summary_SD_0.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_SD_0 = load_summary('summary_SD_0.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract average states at end of simulation\n",
    "SD_0_expo = summary_SD_0.state['expo'].sum(axis=1).mean()\n",
    "SD_0_iasy = summary_SD_0.state['iasy'].sum(axis=1).mean()\n",
    "SD_0_ipre = summary_SD_0.state['ipre'].sum(axis=1).mean()\n",
    "SD_0_isym = summary_SD_0.state['isym'].sum(axis=1).mean()\n",
    "\n",
    "ave_present_seeds = {\n",
    "    'expo' : SD_0_expo,\n",
    "    'iasy' : SD_0_iasy,\n",
    "    'ipre' : SD_0_ipre,\n",
    "    'isym' : SD_0_isym,\n",
    "}\n",
    "\n",
    "with open(f'ave_present_seeds_{c}.txt', 'w+') as logfile:\n",
    "    logfile.write(str(ave_present_seeds) + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "plotter.plot_positives_vs_target(\n",
    "    summary_SD_0, new_cases.sum(axis=1), \n",
    "    test_lag=2,\n",
    "    title='Predictions for observed controlled scenario after March 23', \n",
    "    filename=runstr + 'SD_000',\n",
    "    figsize=FIGSIZE,\n",
    "    errorevery=1, acc=500, \n",
    "    lockdown_at=days_until_lockdown,\n",
    "    ymax=250)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = new_cases.sum(axis=1)\n",
    "targets = np.hstack((targets[0], np.diff(targets)))\n",
    "\n",
    "plotter = Plotter()\n",
    "plotter.plot_daily_infected(\n",
    "    summary_SD_0,\n",
    "    title='', \n",
    "    filename=runstr + 'SD_001',\n",
    "    errorevery=100, acc=1000, \n",
    "    figsize=(6, 4),\n",
    "    start_date='2020-03-10',\n",
    "    show_target=targets,\n",
    "    lockdown_at=days_until_lockdown,\n",
    "    lockdown_label_y=140,\n",
    "    lockdown_label='Lockdown',\n",
    "    ymax=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "plotter.plot_cumulative_infected(\n",
    "    summary_SD_0,\n",
    "    title='', \n",
    "    filename=runstr + 'SD_002',\n",
    "    figsize=(6, 4),\n",
    "    legend_loc='upper left',\n",
    "    errorevery=100, acc=1000, \n",
    "    start_date='2020-03-10',\n",
    "    lockdown_at=days_until_lockdown,\n",
    "    lockdown_label_y=140,\n",
    "    show_target=new_cases.sum(axis=1),\n",
    "    ymax=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "plotter.plot_daily_rts(\n",
    "    summary_SD_0,\n",
    "    filename=runstr + 'SD_003',\n",
    "    start_date='2020-03-10',\n",
    "    sigma=0.37, \n",
    "    figsize=(6, 4), \n",
    "    subplot_adjust=None,\n",
    "    lockdown_label='Lockdown', \n",
    "    lockdown_at=days_until_lockdown,\n",
    "    lockdown_label_y=2.6,\n",
    "    ymax=4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Only simulate the future from here onwards \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# baseline (no measure for the future starting April 12)\n",
    "future_baseline = run(standard_testing(max_time_future), [], max_time_future, present_seeds) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3.2. How long should the current most restrictive measures have to continue, as of April 12?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use settings as above and simulate in the future, with additional measures of the same length\n",
    "\n",
    "continue_for_weeks_SD_2 = [1, 2, 4, 8]\n",
    "testing_params_SD_2 = standard_testing(max_time_future)\n",
    "summaries_SD_2 = []\n",
    "\n",
    "for j, w in enumerate(continue_for_weeks_SD_2):\n",
    "    \n",
    "    # same as tuned before but continued for longer\n",
    "    m = [BetaMultiplierMeasureByType(\n",
    "            t_window=Interval(0.0, (24.0 * 7 * w)), \n",
    "            beta_multiplier=tuned_site_multipliers),\n",
    "\n",
    "         SocialDistancingForAllMeasure(\n",
    "         t_window=Interval(0.0, (24.0 * 7 * w)), \n",
    "            p_stay_home=tuned_p_stay_home)\n",
    "    ] \n",
    "\n",
    "    res = run(testing_params_SD_2, m, max_time_future, present_seeds) \n",
    "    summaries_SD_2.append(res)\n",
    "    print(w, ' done.')\n",
    "\n",
    "save_summary(summaries_SD_2, 'summaries_SD_2.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "continue_for_weeks_SD_2 = [1, 2, 4, 8]\n",
    "summaries_SD_2 = load_summary('summaries_SD_2.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "for j, (w, s) in enumerate(zip(continue_for_weeks_SD_2, summaries_SD_2)):\n",
    "    \n",
    "    plotter.plot_daily_infected(\n",
    "        s,\n",
    "        title=f'Infections with restrictive measures continued for {w} weeks', \n",
    "        filename=runstr + f'SD_01{j}',\n",
    "        figsize=FIGSIZE,\n",
    "        errorevery=14, acc=500, \n",
    "        lockdown_end=w * 7,\n",
    "        show_target=new_cases.sum(axis=1),\n",
    "        ymax=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "plotter.compare_total_infections(\n",
    "    summaries_SD_2, \n",
    "    titles=[f'Extended by {w} weeks' for w in continue_for_weeks_SD_2], \n",
    "    figtitle='',\n",
    "    start_date='2020-04-12',\n",
    "    filename=runstr + 'SD_030', \n",
    "    figsize=(6,4), acc=500, \n",
    "    legend_is_left=True,\n",
    "    ymax=2000, errorevery=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Best values found by MLE for the variance parameter sigma of the Prior\n",
    "best_sigmas = [0.22, 0.18, 0.05, 0.05]\n",
    "\n",
    "plotter = Plotter()\n",
    "plotter.plot_daily_rts(\n",
    "    summaries_SD_2,\n",
    "    filename=runstr + f'SD2_003',\n",
    "    titles=[f'Extended by {w} weeks' for w in continue_for_weeks_SD_2], \n",
    "    colors=plotter.color_different_scenarios,\n",
    "    start_date='2020-04-12',\n",
    "    sigma=best_sigmas, \n",
    "    fill_between=False,\n",
    "    draw_dots=False,\n",
    "    figsize=(6, 4),\n",
    "    errorevery=2,\n",
    "    show_legend=True,\n",
    "    xtick_interval=2,\n",
    "    ymax=4.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3.3. What is the benefit of restricting demographic or random groups?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Experiment (i): Businesses and People below 60 years stop restrictive measures after 2 weeks, people above 60 years continue for X more weeks after that"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use settings as above and simulate in the future, with additional measures of the same length\n",
    "\n",
    "continue_for_weeks_SD_3 = [0, 4, 8, 12]\n",
    "testing_params_SD_3 = standard_testing(max_time_future)\n",
    "summaries_SD_3 = []\n",
    "cutoff_SD_3 = 2\n",
    "for j, w in enumerate(continue_for_weeks_SD_3):\n",
    "    \n",
    "    # same as tuned before but continued for longer\n",
    "    m = [BetaMultiplierMeasureByType(\n",
    "            t_window=Interval(0.0, (24.0 * 7 * cutoff_SD_3)), \n",
    "            beta_multiplier=tuned_site_multipliers),\n",
    "\n",
    "         SocialDistancingForAllMeasure(\n",
    "         t_window=Interval(0.0, (24.0 * 7 * cutoff_SD_3)), \n",
    "            p_stay_home=tuned_p_stay_home),\n",
    "         \n",
    "         SocialDistancingByAgeMeasure(\n",
    "         t_window=Interval((24.0 * 7 * cutoff_SD_3), \n",
    "                           (24.0 * 7 * (cutoff_SD_3  + w))), \n",
    "         p_stay_home=[0.0, 0.0, 0.0, 0.0, 0.9, 0.9])\n",
    "    ] \n",
    "\n",
    "    res = run(testing_params_SD_3, m, max_time_future, present_seeds) \n",
    "    summaries_SD_3.append(res)\n",
    "    print(w, ' done.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_summary(summaries_SD_3, 'summaries_SD_3.pk')\n",
    "summaries_SD_3 = load_summary('summaries_SD_3.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "titles_SD_3 = [f'Extended by {w} weeks' for w in continue_for_weeks_SD_3]\n",
    "plotter.compare_total_fatalities_and_hospitalizations(\n",
    "    summaries_SD_3, \n",
    "    titles=titles_SD_3, \n",
    "    figtitle='Hospitalizations and Fatalities when social distancing of ages 60+',\n",
    "    filename=runstr + 'SD_040', \n",
    "    figsize=FIGSIZE,\n",
    "    acc=500,\n",
    "    ymax=50, errorevery=14)\n",
    "\n",
    "plotter.compare_total_infections(\n",
    "    summaries_SD_3, \n",
    "    titles=titles_SD_3, \n",
    "    figtitle=r'Total infections when social distancing of ages 60+',\n",
    "    filename=runstr + 'SD_041', \n",
    "    figsize=FIGSIZE, acc=500, \n",
    "    ymax=500, errorevery=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Map plot of change in activity patterns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bbox = (48.4900, 48.5485, 9.0224, 9.1061) # coordinate bounding box\n",
    "# ill = MapIllustrator()\n",
    "# map_name = 'Map title'\n",
    "# map_time = 24.0 * 30\n",
    "# map_summary = summaries_44_1['advanced'][0]\n",
    "# sitemap = ill.checkin_rate_map(\n",
    "#     bbox=bbox, site_loc=mob.site_loc, map_name=map_name, \n",
    "#     sim=map_summary, t=map_time, r=0)\n",
    "# sitemap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Experiment (ii): Alternating quarantine of K random subsets of the population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K_groups = [2, 3, 4]\n",
    "testing_params_SD_4 = standard_testing(max_time_future)\n",
    "summaries_SD_4 = []\n",
    "K_groups_weeks = 6\n",
    "\n",
    "# no quarantine baseline (standard measures)\n",
    "summaries_SD_4.append(future_baseline)\n",
    "\n",
    "# alternating quarantines\n",
    "for j, k in enumerate(K_groups):\n",
    "    \n",
    "    # same as tuned plus containing K groups alternatinglye in 4 week period afterwards\n",
    "    m = [SocialDistancingForKGroups(\n",
    "         t_window=Interval(0.0, (24.0 * 7 * K_groups_weeks)), \n",
    "         K=k)] \n",
    "\n",
    "    res = run(testing_params_SD_4, m, max_time_future, present_seeds) \n",
    "    summaries_SD_4.append(res)\n",
    "    print(k, ' done.')\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_summary(summaries_SD_4, 'summaries_SD_4.pk')\n",
    "summaries_SD_4 = load_summary('summaries_SD_4.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "titles_SD_4 = ['No curfew'] + [f'K = {k}' for k in K_groups]\n",
    "plotter.compare_total_infections(\n",
    "    summaries_SD_4, \n",
    "    titles=titles_SD_4, \n",
    "    figtitle=r'Curfew of K random groups for staggered K - 1 day windows',\n",
    "    filename=runstr + 'SD_050', \n",
    "    figsize=FIGSIZE, acc=500, \n",
    "    ymax=500, errorevery=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3.4. Can we control the outbreak using only contact tracing and isolation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_params_SD_5 = standard_testing(max_time_future)\n",
    "testing_params_SD_5['test_smart_delta'] = 24.0 * 3 # time window considered for inspecting contacts\n",
    "testing_params_SD_5['test_smart_action'] = 'isolate'\n",
    "testing_params_SD_5['test_targets'] = 'isym'\n",
    "\n",
    "isolation_days = [3, 7] # how many days selected people have to stay in isolation\n",
    "contacts_isolated = [10, 25] # how many contacts are isolated in the `test_smart_delta` window\n",
    "duration_weeks_SD_5 = 12 # strategies tested for 12 weeks starting today\n",
    "\n",
    "summaries_SD_5 = dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# baseline taken from above\n",
    "for days in isolation_days:\n",
    "    for contacts in contacts_isolated:\n",
    "        summaries_SD_5[(days, contacts)] = [future_baseline]\n",
    "\n",
    "\n",
    "# Isolation strategies\n",
    "for days in isolation_days:\n",
    "    for contacts in contacts_isolated:\n",
    "        for policy in ['basic', 'advanced']:\n",
    "\n",
    "            testing_params_SD_5['smart_tracing'] = policy\n",
    "            testing_params_SD_5['test_smart_num_contacts'] = contacts\n",
    "\n",
    "\n",
    "            # same as tuned plus different isolation strategies for contact tracing\n",
    "            m = [SocialDistancingForSmartTracing(\n",
    "                    t_window=Interval(*testing_params_SD_5['testing_t_window']), \n",
    "                    p_stay_home=1.0, \n",
    "                    test_smart_duration=24.0 * days)\n",
    "            ]\n",
    "\n",
    "            res = run(testing_params_SD_5, m, max_time_future, present_seeds) \n",
    "            summaries_SD_5[(days, contacts)].append(res)\n",
    "            \n",
    "            print(days, contacts, policy, ' done.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_summary(summaries_SD_5, 'summaries_SD_5.pk')\n",
    "summaries_SD_5 = load_summary('summaries_SD_5.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "\n",
    "for i, days in enumerate(isolation_days):\n",
    "    for j, contacts in enumerate(contacts_isolated):\n",
    "        summaries_ = summaries_SD_5[(days, contacts)]\n",
    "        titles_SD_5_ = ['No contact tracing', 'Basic contact tracing', 'Advanced contact tracing']\n",
    "        plotter.compare_total_infections(\n",
    "            summaries_, \n",
    "            titles=titles_SD_5_, \n",
    "            figtitle=(f'Strategies for isolating {contacts} contacts.\\n'\n",
    "                      f'Selected people are isolated for {days} days.'),\n",
    "            filename=runstr + f'SD_5{i}{j}', \n",
    "            figsize=FIGSIZE, acc=500, \n",
    "            ymax=300, errorevery=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "\n",
    "for i, days in enumerate(isolation_days):\n",
    "    for j, contacts in enumerate(contacts_isolated):\n",
    "        summaries_ = summaries_SD_5[(days, contacts)]\n",
    "        titles_SD_5_ = ['No contact tracing', 'Basic contact tracing', 'Advanced contact tracing']\n",
    "        for k, titl in enumerate(titles_SD_5_):\n",
    "            plotter.plot_daily_at_home(\n",
    "                summaries_[k],\n",
    "                title='Home isolation for ' + titl, \n",
    "                filename=runstr + f'SD_6{i}{j}{k}', figsize=FIGSIZE, \n",
    "                errorevery=14, acc=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3.5. Effects  of compliance on the efficacy of isolation for smart  tracing strategies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_params_SD_6 = standard_testing(max_time_future)\n",
    "testing_params_SD_6['test_smart_delta'] = 24.0 * 3 # time window considered for inspecting contacts\n",
    "testing_params_SD_6['test_smart_action'] = 'isolate'\n",
    "testing_params_SD_6['test_targets'] = 'isym'\n",
    "testing_params_SD_6['test_smart_num_contacts'] = 25\n",
    "isolation_days_SD_6 = 7 # how many days selected people have to stay in isolation\n",
    "duration_weeks_SD_6 = 12 # strategies tested for 12 weeks starting today\n",
    "\n",
    "summaries_SD_6 = dict()\n",
    "\n",
    "p_compliance = [0.1, 0.25, 0.5, 1.0] \n",
    "\n",
    "for policy in ['basic', 'advanced']:\n",
    "    summaries_ = []\n",
    "    testing_params_SD_6['smart_tracing'] = policy\n",
    "    \n",
    "    for p in p_compliance:\n",
    "\n",
    "        m = [SocialDistancingForSmartTracing(\n",
    "                t_window=Interval(*testing_params_SD_6['testing_t_window']), \n",
    "                p_stay_home=1.0, \n",
    "                test_smart_duration=24.0 * isolation_days_SD_6),\n",
    "             ComplianceForAllMeasure(\n",
    "                t_window=Interval(*testing_params_SD_6['testing_t_window']), \n",
    "                p_compliance=p)\n",
    "        ]\n",
    "        res = run(testing_params_SD_6, m, max_time_future, present_seeds) \n",
    "        summaries_.append(res)\n",
    "\n",
    "        print(policy, p, ' done.')\n",
    "        \n",
    "    summaries_SD_6[policy] = summaries_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_summary(summaries_SD_6, 'summaries_SD_6.pk')\n",
    "summaries_SD_6 = load_summary('summaries_SD_6.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "titles_SD_6_ = list(['Tracking compliance '+ str(int(p*100.0)) + ' \\%' for p in p_compliance])\n",
    "for j, policy in enumerate(['basic', 'advanced']):\n",
    "    summaries_ = summaries_SD_6[policy]\n",
    "    plotter.compare_total_infections(\n",
    "        summaries_, \n",
    "        titles=titles_SD_6_, \n",
    "        figtitle=(f'Infections for compliance levels for ' + policy + ' tracing isolation'),\n",
    "        filename=runstr + f'SD_6{j}', \n",
    "        figsize=FIGSIZE, acc=500, \n",
    "        ymax=200, errorevery=14)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.4 Testing and Tracing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.4.1. Effects of smart testing strategies using contact tracing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_params_TT_0 = standard_testing(max_time_future)\n",
    "testing_params_TT_0['test_smart_delta'] = 24.0 * 3 # time window considered for inspecting contacts\n",
    "testing_params_TT_0['test_smart_action'] = 'test'\n",
    "testing_params_TT_0['test_targets'] = 'isym'\n",
    "\n",
    "contacts_inspected = [10, 25] # how many contacts are inspected in the `test_smart_delta` window\n",
    "duration_weeks_TT_0 = 12 # strategies tested for 6 weeks starting today\n",
    "\n",
    "summaries_TT_0 = dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# baseline taken from above\n",
    "for contacts in contacts_inspected:\n",
    "    summaries_TT_0[contacts] = [future_baseline]\n",
    "\n",
    "# Testing strategies\n",
    "for contacts in contacts_inspected:\n",
    "    for policy in ['basic', 'advanced']:\n",
    "        testing_params_TT_0['smart_tracing'] = policy\n",
    "        testing_params_TT_0['test_smart_num_contacts'] = contacts\n",
    "        m = []\n",
    "        res = run(testing_params_TT_0, m, max_time_future, present_seeds) \n",
    "        summaries_TT_0[contacts].append(res)\n",
    "\n",
    "        print(contacts, policy, ' done.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_summary(summaries_TT_0, 'summaries_TT_0.pk')\n",
    "summaries_TT_0 = load_summary('summaries_TT_0.pk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = Plotter()\n",
    "for j, contacts in enumerate(contacts_inspected):\n",
    "    summaries_ = summaries_TT_0[contacts]\n",
    "    titles_TT_0_ = ['No contact tracing', 'Basic contact tracing', 'Advanced contact tracing']\n",
    "    plotter.compare_total_infections(\n",
    "        summaries_, \n",
    "        titles=titles_TT_0_, \n",
    "        figtitle=(f'Tracing strategies when testing {contacts} contacts'),\n",
    "        filename=runstr + f'TT_00{j}', \n",
    "        show_positives=True,\n",
    "        figsize=FIGSIZE, acc=500, \n",
    "        ymax=200, errorevery=14)\n",
    "    \n",
    "    for l, titl in enumerate(titles_TT_0_):\n",
    "        plotter.plot_daily_tested(summaries_[l], \n",
    "            title='Test results for ' + titl, \n",
    "            filename=runstr + f'TT_1{j}{k}', \n",
    "            figsize=FIGSIZE, errorevery=14, acc=500, ymax=200)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}