Demo_Vaccination.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vaccination (Demo)\n",
    "In order to demonstrate the simulations, which were performed to evaluate the different vaccination-strategies in our [preprint](https://www.medrxiv.org/content/10.1101/2020.12.16.20248301v1), we generated a jupyter notebook. We suggest to get familiar with the `Demo.ipynb` notebook first, in order to get an intuition for our model and its basic application and characteristics.\n",
    "\n",
    "This notebook subdivied the following code blocks:\n",
    "\n",
    "1. [\"Initiate world\"](#buildW) - initialize small or large world for Gangelt (small = 10% of population and buildings, large = 100% of population and buildings)\n",
    "\n",
    "2. [Simulating an initial infection-wave](#runIni) to be used in defining one of the tested strategies\n",
    "\t—> Be aware that here exists a checkpoint for the user \n",
    "    \n",
    "3. [Define vaccination orders](#defOrder) for different vaccination strategies.\n",
    "\n",
    "4. [Running simulations with different vaccination-fractions](#runVacScreen) of the population, for the different strategies.\n",
    "\n",
    "5. [Generating a plot](#plots), resembling figure 3 in our manuscript.\n",
    "\n",
    "Please be aware, that we set the default world to be used to the small version, in order to avoid excessive runtimes for this demonstration. \n",
    "Since this world is not fully representative to the large world, we used for the manuscript; the infection-dynamics differ and thus the results of the vaccintion-screens do too.\n",
    "Furthermore we have set the vaccination-fraction increments to 20% (where we used 5%-steps in the manuscript), also to reduce runtime.\n",
    "\n",
    "#### Be aware: Excessive Runtime!\n",
    "The (default) reduced version of the vaccination screens has a runtime of around **1-2 hours**, but there exists the possibility to use the non-reduced version (as we did in the manuscript); however be aware, that we expect a runtime of **12-36 hours** for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gerda.core.virusPropagationModel import *\n",
    "import gerda.utilities.VPM_plotting as vpm_plot\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "import glob\n",
    "import os\n",
    "import copy\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.cm as cm\n",
    "import numpy as np\n",
    "#%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "t0=time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build world <a id='buildW'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geopath = 'input_data/geo/'\n",
    "geofiles = {\n",
    "            0:'Buildings_Gangelt_MA_1.csv',\n",
    "            1:'Buildings_Gangelt_MA_3.csv',\n",
    "}         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choose which world to initate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we choose, which version of the world to use for the vaccination-screens in this notebook.\n",
    "By default we use the reduced version of the world, in order to avoid excessive runtimes (reduced version approx. 2h, complete version approx. 12h).\n",
    "However one can comment out the second line and uncomment the first line, to run the screen on the large world, we based the results in the manuscript on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#world_to_pick = 0 ## choosing complete version\n",
    "world_to_pick = 1 ## choosing reduced version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initiate the selected world without infections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_initially_infected=0\n",
    "modeledWorld = ModeledPopulatedWorld(initial_infections=n_initially_infected, geofile_name=geopath+geofiles[world_to_pick], input_schedules='schedules_v2')\n",
    "n_people = modeledWorld.number_of_people\n",
    "print('Number of agents: '+ str(n_people))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of the function to initially infect specified susceptible agents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All agents with IDs in the list 'IDs' will be initially infected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def infect_world(world, IDs=[1]):\n",
    "    ID_list = world.get_remaining_possible_initial_infections(IDs)\n",
    "    world.initialize_infection(specific_people_ids=ID_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run initial simulation<a id='runIni'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we run an intitial (unperturbed) infection wave, to be used in the vaccination strategy which is based on the forecasted infection wave."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define Simulation characteristics\n",
    "\n",
    "timesteps_initial_sim: Number of simulated timesteps in hours for the initial simulation <br>\n",
    "general_infectivity: Scaling factor for the general infectivity (default value = 0.3) <br>\n",
    "general_interaction_frequency: Expected interaction-frequency over all individuals per hour (default-value=2/h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timesteps_initial_sim=3000\n",
    "general_infectivity=0.3\n",
    "general_interaction_frequency=2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Introduce infected agents\n",
    "\n",
    "Before starting the simulation the generated world is copied to ensure the same origin for each simulation. The number of initial infections is defined by 'n_initially_infected' (default=4)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world_to_use = copy.deepcopy(modeledWorld)\n",
    "\n",
    "n_initially_infected=4\n",
    "\n",
    "infect_world(world_to_use, IDs=[i+1 for i in range(n_initially_infected)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set up simulation with specified characteristics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_initial_wave = Simulation(world_to_use, timesteps_initial_sim, run_immediately=False)\n",
    "\n",
    "simulation_initial_wave.change_agent_attributes({'all':{'behaviour_as_infected':{'value':general_infectivity,'type':'replacement'}}})\n",
    "simulation_initial_wave.interaction_frequency=general_interaction_frequency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_initial_wave.simulate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CHECKPOINT FOR USER:\n",
    "\n",
    "If the plot below does not show an infection outbreak (due to the stochastic nature of our model, some infection-waves may die out); please re-run the two cells above and continue if the plot shows an outbreak."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot status-trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_initial_wave.plot_status_timecourse()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define vaccination orders for different vaccination strategies <a id='defOrder'></a>\n",
    "\n",
    "Vaccination orders are defined by ordered lists of agent IDs following specified sorting rules.\n",
    "These are:\n",
    "1. sorted by [forecasted infection wave](#sorted)\n",
    "2. sorted by [overrepresentation](#overrep) of agent type in infection spreading\n",
    "3. sorted by [agent interactivity](#activity)\n",
    "4. sorted by [households](#house)\n",
    "5. sorted by [age](#age) (descending)\n",
    "6. sorted [randomly](#random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sorting strategies, based on previous simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following two vaccination strategies are based on results derived from predceeding simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Sorted by forecasted infection wave <a id='sorted'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we return the order, in which agents got infected in the previous simulation (simulation_initial_wave).\n",
    "This gives us a vaccination-order, based on the simulated outbreak."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to sort agents according to forecasted infection wave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_previous_infections(previous_simulation):\n",
    "    df_inf = previous_simulation.get_infection_event_information()\n",
    "    infected = list(df_inf['h_ID'].values)\n",
    "    all_agents = infected+[a.ID for a in previous_simulation.people if a.ID not in infected]\n",
    "    return all_agents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create list of agents sorted by forecasted infection wave, by executing defined function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agents_forecasted = get_previous_infections(simulation_initial_wave)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Sorted by over-representation (agent types) <a id='overrep'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we sort agents by the the number of infections caused, relative to the number of individuals of the respective type (termed over-representation). Order of individual agents of the repective is not specified."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to sort agents by type. The sequence of types is defined in sorted list \"agent_types\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ids_by_agent_type(world,agent_types_sorted_by_overrepresentation):\n",
    "    ids_by_type = {s_type:[] for s_type in agent_types_sorted_by_overrepresentation}\n",
    "    for p in world.people:\n",
    "        ids_by_type[p.schedule['type']].append(p.ID)\n",
    "    ordered_ids = []\n",
    "    for s_type in agent_types_sorted_by_overrepresentation:\n",
    "        ordered_ids.extend(ids_by_type[s_type])\n",
    "    return ordered_ids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create list of agents sorted by type, by executing defined function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that the order of agent types, acording to over-representation, is not determined automatically in this notebook. Instead we used the order, applied in the manuscript and hard-coded it in the variable \"agent_types_sorted_by_overrepresentation\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agent_types_sorted_by_overrepresentation = ['under_age', 'adult', 'teacher', 'medical_professional', 'public_worker', 'pensioner']\n",
    "agents_by_type = get_ids_by_agent_type(modeledWorld,agent_types_sorted_by_overrepresentation=agent_types_sorted_by_overrepresentation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Sorted by Interactivity <a id='activity'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we sort agents, based on their \"interactivity\". This is determined by performing a preceding simulation of 168 hours (1 week), without any introduced infection (so the determining the interactivity is not affected by the infection-spread); and subsequntly sorting agents according to the number of interactions, they had."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to sort agents with decreasing number of total interactions from a preceding simulation without infection spreading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ids_by_interactions(simulation_no_infections):\n",
    "    contacts = simulation_no_infections.get_contact_distributions(min_t=0, max_t=168)[1]\n",
    "    contacts_mean = contacts.groupby('ID').mean()\n",
    "    contacts_mean.reset_index(inplace=True)\n",
    "    contacts_sorted = contacts_mean.sort_values('interactions', axis=0, ascending=False)\n",
    "    to_recover_list = list(contacts_sorted['ID'].values)\n",
    "    return to_recover_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the preceding simulation without infection spreading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timesteps_blank_sim=168\n",
    "uninfectedWorld = copy.deepcopy(modeledWorld)\n",
    "simulation_no_infections = Simulation(uninfectedWorld, timesteps_blank_sim, run_immediately=False)\n",
    "simulation_no_infections.interaction_frequency=general_interaction_frequency\n",
    "simulation_no_infections.simulate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create the list of agents sorted by decreasing number of total interactions, by executing defined function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agents_by_interactions = get_ids_by_interactions(simulation_no_infections)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sorting-strategies, based on population-structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following strategies are based on the population structure of the simulated world and do not require preceeding simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Sorted by household <a id='house'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define function to sort agents by their households. First one member of each household is chosen. When there is one member for each household a second member is chosen for each household of size 2 and bigger. After that a third member is chosen for each household of size 3 and bigger. ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ids_by_households(world):\n",
    "    ai = world.get_agent_info()\n",
    "    home_count_dict = dict(zip(list(ai.groupby('Home').count().index), list(ai.groupby('Home').count()['ID'])))\n",
    "    ai['Home_size']=ai['Home'].map(home_count_dict)\n",
    "\n",
    "    ai_sorted = ai.sort_values('Home')\n",
    "    home_list = ai_sorted['Home'].values\n",
    "\n",
    "    count_list=[]\n",
    "    predecessor = 0\n",
    "    count = 1\n",
    "    for x in home_list:\n",
    "        if x == predecessor:\n",
    "            count +=1\n",
    "        else:\n",
    "            count = 1\n",
    "        count_list.append(count)\n",
    "        predecessor = x\n",
    "\n",
    "    ai_sorted['Home_position'] = count_list\n",
    "\n",
    "    ai_sorted.sort_values('Home_position') ## dataframe\n",
    "    households_sorted = ai_sorted.sort_values(by=['Home_position'])\n",
    "    to_recover_list = list(households_sorted['ID']) ### id list\n",
    "    return to_recover_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create the list of agents sorted by households, by executing defined function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agents_by_households = get_ids_by_households(modeledWorld)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Sorted by age <a id='age'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we sort agents by age (in descending order) and thus implicitly by likelyhood for severe consequences of infection, to vaccinate accordingly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to sort agents by age, descending"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ids_by_age(world):\n",
    "    ages_and_ids = [(p.age, p.ID) for p in world.people]\n",
    "    ages_and_ids_df = pd.DataFrame(ages_and_ids, columns=['age','id'])\n",
    "    ages_and_ids_df.sort_values('age', ascending=False, inplace=True)\n",
    "    to_recover_list = list(ages_and_ids_df['id'])\n",
    "    return to_recover_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create the list of agents sorted by age, descending"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agents_by_age = get_ids_by_age(modeledWorld)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Random sorting <a id='random'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last vaccination-strategy is to vaccinate individuals in random order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to randomly shuffle the order of agents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_random_id_list(world):\n",
    "    agent_ids = [p.ID for p in world.people]\n",
    "    random.shuffle(agent_ids)\n",
    "    return agent_ids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a list of agents with random order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agents_random = get_random_id_list(modeledWorld)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Old then active sorting <a id='old_active'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bonus vaccination-strategy is to vaccinate old people first and then switch to interactive individuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "old_age_limit=55"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to vaccinate old people first and at a given age limit switch to vaccinate people by decreasing interactivity. It reuses the list of agents sorted by interactivity from step 3 above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_old_active_id_list(world):\n",
    "    old_ages_and_ids = [(p.age, p.ID) for p in world.people if p.age >= old_age_limit]\n",
    "    old_ages_and_ids_df = pd.DataFrame(old_ages_and_ids, columns=['age','id'])\n",
    "    old_ages_and_ids_df.sort_values('age', ascending=False, inplace=True)\n",
    "    old_agents_list = list(old_ages_and_ids_df['id'])\n",
    "    young_agents_by_interaction = [id for id in agents_by_interactions if id not in old_agents_list]\n",
    "    agent_list = old_agents_list + young_agents_by_interaction\n",
    "    return agent_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a list of agents starting with oldest and below a certain age sorted by interactivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agents_old_active = get_old_active_id_list(modeledWorld)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run vaccination-fraction screens for different strategies<a id='runVacScreen'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With ordered list of individuals, by different vaccination-strategies; we now perform the screen of vaccination-fractions among the population, in specified increments. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to vaccinate a certain fraction of agents in a specified order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function vaccinates (setting initial agent-state to recovered) a given fraction of individuals in the world to be simulated (argument \"world\") (specified by argument \"frac\"), according to the strategy-specific agent-order (argument \"ordered_agents\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define function to vaccinate list of agents\n",
    "def vaccinate_world(world, frac, ordered_agents):\n",
    "    agents_to_vaccinate = ordered_agents[:int(frac*len(ordered_agents))]\n",
    "    for p in world.people:\n",
    "        if p.ID in agents_to_vaccinate:\n",
    "            p.set_initially_recovered()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a dictionary to access the previously created ordered lists by keywords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agent_orders = {'forecasted':agents_forecasted, 'overrepresentation':agents_by_type,\n",
    "                'interactions':agents_by_interactions, 'households':agents_by_households,\n",
    "                'age':agents_by_age, 'random':agents_random, 'old_active':agents_old_active}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define the fractions of agents to vaccinate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define function to simulate the vaccination scenario with the defined fractions and agent orders"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This cell defines the function, to run a vaccination-fraction screen for one strategy.\n",
    "It takes a world to simulate (\"modeledWorld\"), a list of ordered agents to be vaccinated (\"ordered_agents\"), a list of vaccination-fractions to be tested (\"fractions\"), the number of seed-infections (\"n_initially_infected\" by default 4) and the number of time-steps to be simulated (\"timesteps_vaccination_sim\" by default 3000) as arguments. \n",
    "For each specified fractions it vaccinates the first n agents in the list of ordered agents, sets the specified number of seed-infections (making sure none of these are vaccinated), sets the general infectivity (\"behaviour_as_infected\") to the default value (specified by perviously defined variable \"general_infectivity\") and finally  does the same for the interactivity (set to default value \"general_interaction_frequency\"). Then the system is simulated and relevant results (total number of total infections, total number of deceased and the maximum demand for ICU-treatment) are returned. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_fractions(modeledWorld, ordered_agents, fractions, n_initially_infected=4, timesteps_vaccination_sim=3000):\n",
    "    results=pd.DataFrame(columns=['fraction','infected','icued','dead'])\n",
    "    for frac in fractions:\n",
    "        world_to_use = copy.deepcopy(modeledWorld)\n",
    "        vaccinate_world(world_to_use, frac, ordered_agents)\n",
    "        infect_world(world_to_use, IDs=[i+1 for i in range(n_initially_infected)])\n",
    "        simulation_vaccinated_fraction = Simulation(world_to_use, timesteps_vaccination_sim, run_immediately=False)\n",
    "        simulation_vaccinated_fraction.change_agent_attributes({'all':{'behaviour_as_infected':{'value':general_infectivity,'type':'replacement'}}})\n",
    "        simulation_vaccinated_fraction.interaction_frequency=general_interaction_frequency\n",
    "        simulation_vaccinated_fraction.simulate()\n",
    "        cummulative_flags = simulation_vaccinated_fraction.get_flag_sums_over_time(specific_flags=['ICUed'])\n",
    "        dead_agents = simulation_vaccinated_fraction.get_status_trajectories(specific_statuses=['D'])['D']['D'].iloc[-1]\n",
    "        infected_agents = simulation_vaccinated_fraction.get_infection_event_information().shape[0]\n",
    "        icued_agents = cummulative_flags['ICUed'].max()\n",
    "        results.loc[frac,'fraction']=frac\n",
    "        results.loc[frac,'infected']=infected_agents\n",
    "        results.loc[frac,'icued']=icued_agents\n",
    "        results.loc[frac,'dead']=dead_agents\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a list of strategies to test by keywords, run the simulations for all fractions and all strategies and save relevant results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below defines the list of strategies, one wants to test (by default all), sets the number of simulated timesteps for each run (recommended default is 3000) and finally defines the range and increment sizes of vaccination-fractions.\n",
    "Then it runs the simulate_fractions-function for each strategy and writes te resulting pandas.DataFrames in a dictionary (results_dict) with the strategy-names as keys."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#strategies_to_test = ['forecasted', 'overrepresentation', 'interactions', 'households', 'age', 'random']\n",
    "strategies_to_test = ['interactions', 'age', 'random', 'old_active']\n",
    "timesteps_vaccination_sim=3000\n",
    "\n",
    "lowest_vaccination_fraction=0.1\n",
    "highest_vaccination_fraction=0.91\n",
    "increment_size=0.2\n",
    "\n",
    "fractions = np.arange(lowest_vaccination_fraction, highest_vaccination_fraction, increment_size)\n",
    "\n",
    "results_dict = {strategy: simulate_fractions(modeledWorld, agent_orders[strategy], fractions,timesteps_vaccination_sim=timesteps_vaccination_sim) for strategy in strategies_to_test}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot results of vaccination screens<a id='plots'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we generate a plot, representing figure 3 in the manuscript. The very left panel corresponds to Panel A in figure 3 and shows the number of total infections, relative to the number of non-vaccinated individuals. The middle panel corresponds to panel 3B, showing the same for deceased individuals. The very right panel corresponds to panel 3D and shows the maximum ICU-demand per 100000. Please note the Panel 3C from the manuscript can not be obtained in this notebook, since it relies on several replicates (100 in our study) of each simulation; and is thus infeasible for demonstration-purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,3, figsize=(15,5))\n",
    "\n",
    "colors = {'random':'lightseagreen', 'interactions':'violet', 'age':'gold', 'old_active':'black',\n",
    "          'overrepresentation':'royalblue', 'households':'lightgreen', 'forecasted':'orange'}\n",
    "\n",
    "for strategy in results_dict:\n",
    "    fractions = results_dict[strategy]['fraction']\n",
    "    axs[0].plot(fractions, results_dict[strategy]['infected']/(modeledWorld.number_of_people*(1-fractions)),\n",
    "                label=strategy, lw=2, alpha=0.7, color=colors[strategy])\n",
    "    axs[1].plot(fractions, results_dict[strategy]['dead']/(modeledWorld.number_of_people*(1-fractions)),\n",
    "                label=strategy, lw=2, alpha=0.7, color=colors[strategy])\n",
    "    axs[2].plot(fractions, (results_dict[strategy]['icued']/modeledWorld.number_of_people)*100000.0,\n",
    "                label=strategy, lw=2, alpha=0.7, color=colors[strategy])\n",
    "    \n",
    "for ax in axs:\n",
    "    ax.set_xlabel('fraction of immunized population')\n",
    "    \n",
    "axs[0].set_ylabel('fraction of non immunized population')\n",
    "axs[1].set_ylabel('fraction of non immunized population')\n",
    "axs[2].set_ylabel('max ICU demand per 100.000')\n",
    "\n",
    "axs[0].set_title('Infections')\n",
    "axs[1].set_title('Death')\n",
    "axs[2].set_title('ICU Demand')\n",
    "\n",
    "axs[0].legend()\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('\\nage limit: ', old_age_limit)\n",
    "print('\\nmax vacc: ', highest_vaccination_fraction)\n",
    "for strategy in results_dict:\n",
    "    print('\\n', strategy, '\\n', results_dict[strategy], '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(modeledWorld.number_of_people)\n",
    "# print(len(agents_by_age))\n",
    "# print(len(agents_by_interactions))\n",
    "print(len(agents_old_active))\n",
    "\n",
    "print(time.time()-t0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}