cruiseship/simulation/figure3c_efficacy_scenario.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mtick\n",
    "import pandas as pd\n",
    "from joblib import Parallel, delayed\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "cwd = os.getcwd()\n",
    "model_path = os.path.abspath(os.path.join(cwd, os.pardir, 'models'))\n",
    "sys.path.append(model_path)\n",
    "from cumul_diamond_ODE_model import *\n",
    "from cumul_treated_diamond_ODE_model import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "untreated_model = cumulDiamondODEModel()\n",
    "treated_model = cumulTreatedDiamondODEModel()\n",
    "times = np.linspace(0,31,100)\n",
    "n_simulations = 100\n",
    "nu = 2 / 4.3\n",
    "inferred_parameters = pd.read_csv('../data/inferred_parameters.csv')\n",
    "data = pd.read_csv('../data/data_cases_symp.csv')\n",
    "data['onset_date'] = pd.to_datetime(data['onset_date'], format='%d-%b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pre-sample parameter sets to ensure both models use the same parameters for each simulation\n",
    "sampled_indices = np.random.randint(0, inferred_parameters.shape[0], size=n_simulations)\n",
    "sampled_parameters = [inferred_parameters.iloc[idx].to_list() for idx in sampled_indices]\n",
    "\n",
    "coverage_array = [(1, 1), (1, 0), (0, 1), (0.5, 0.5)]\n",
    "efficacy_array = np.linspace(0, 1, 100)\n",
    "\n",
    "def simulate_models(params, times):\n",
    "    untreated_values = untreated_model.simulate(params[:8], times)\n",
    "    cumul_infected_untreated_passengers = untreated_values[-1,22]\n",
    "    cumul_infected_untreated_crew = untreated_values[-1,23]\n",
    "    total_untreated_infections = cumul_infected_untreated_passengers + cumul_infected_untreated_crew\n",
    "\n",
    "    averted_results = {}\n",
    "    for eff in efficacy_array:\n",
    "        for coverage in coverage_array:\n",
    "            treated_values = treated_model.simulate(params[:8] + [eff, eff, coverage[0], coverage[1]], times)\n",
    "            cumul_infected_treated_passengers = treated_values[-1,22]\n",
    "            cumul_infected_treated_crew = treated_values[-1,23]\n",
    "            total_treated_infections = cumul_infected_treated_passengers + cumul_infected_treated_crew\n",
    "\n",
    "            infections_averted_percentage = 100 * (total_untreated_infections - total_treated_infections) / total_untreated_infections if total_untreated_infections > 0 else 0\n",
    "            \n",
    "            averted_results.setdefault((eff, coverage[0], coverage[1]), []).append(infections_averted_percentage)\n",
    "\n",
    "    return averted_results\n",
    "\n",
    "# Execute simulations in parallel\n",
    "all_results = Parallel(n_jobs=-1)(delayed(simulate_models)(params, times) for params in sampled_parameters)\n",
    "\n",
    "# Combine results from all simulations\n",
    "combined_results = {}\n",
    "for result in all_results:\n",
    "    for key, values in result.items():\n",
    "        combined_results.setdefault(key, []).extend(values)\n",
    "\n",
    "# Calculate percentiles for each combination of efficacy and coverage\n",
    "percentile_results = []\n",
    "for key, values in combined_results.items():\n",
    "    efficacy, passenger_coverage, crew_coverage = key\n",
    "    median_averted = np.percentile(values, 50)\n",
    "    p5_averted = np.percentile(values, 5)\n",
    "    p95_averted = np.percentile(values, 95)\n",
    "    percentile_results.append((efficacy, passenger_coverage, crew_coverage, median_averted, p5_averted, p95_averted))\n",
    "\n",
    "# Convert results to a DataFrame\n",
    "columns = [\"Efficacy\", \"Passenger Coverage\", \"Crew Coverage\", \"Median Percentage Infections Averted\", \"5th Percentile Percentage Infections Averted\", \"95th Percentile Percentage Infections Averted\"]\n",
    "results_df = pd.DataFrame(percentile_results, columns=columns)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting the color palette to viridis\n",
    "colors = plt.cm.viridis(np.linspace(0, 1, len(coverage_array)))\n",
    "\n",
    "# Creating the plot\n",
    "plt.figure(figsize=(6, 5))\n",
    "plt.rcParams['font.sans-serif']=='Arial'\n",
    "\n",
    "for i, coverage in enumerate(coverage_array):\n",
    "    # Filtering the data for the current coverage scenario\n",
    "    scenario_data = results_df[(results_df['Passenger Coverage'] == coverage[0]) & (results_df['Crew Coverage'] == coverage[1])]\n",
    "\n",
    "    # Sorting by efficacy for consistent plotting\n",
    "    scenario_data = scenario_data.sort_values('Efficacy')\n",
    "\n",
    "    # Plotting the median percentage of infections averted\n",
    "    plt.plot(scenario_data['Efficacy'] * 100, scenario_data['Median Percentage Infections Averted'], color=colors[i], label=f'Passenger: {int(coverage[0]*100)}%, Crew: {int(coverage[1]*100)}%')\n",
    "\n",
    "    # Adding the shaded region for the 5th-95th percentile\n",
    "    plt.fill_between(scenario_data['Efficacy'] * 100, scenario_data['5th Percentile Percentage Infections Averted'], scenario_data['95th Percentile Percentage Infections Averted'], color=colors[i], alpha=0.3)\n",
    "\n",
    "\n",
    "plt.plot(np.nan, np.nan, color='black', label='Median')\n",
    "plt.fill_between([np.nan,np.nan], [np.nan, np.nan], alpha=0.3, color='black', label='90% confidence interval')\n",
    "\n",
    "# Formatting the y-axis to show percentage\n",
    "plt.gca().yaxis.set_major_formatter(mtick.PercentFormatter())\n",
    "\n",
    "# Formatting the x-axis to show percentage\n",
    "plt.gca().xaxis.set_major_formatter(mtick.PercentFormatter())\n",
    "\n",
    "# Adjusting y-axis limits\n",
    "plt.ylim(0, 100)\n",
    "plt.xlim(0, 100)\n",
    "\n",
    "# Adding labels and title\n",
    "plt.xlabel('Efficacy', weight='bold', fontsize=12)\n",
    "plt.ylabel('Infections averted', weight='bold', fontsize=12)\n",
    "plt.legend(title='Coverage', loc='lower right', fontsize=8)\n",
    "\n",
    "# Show the plot\n",
    "plt.tight_layout()\n",
    "\n",
    "# Define the directory path and the filename\n",
    "directory = '../figures'\n",
    "filename = 'figure3c.svg'\n",
    "\n",
    "# Check if the directory exists\n",
    "if not os.path.exists(directory):\n",
    "    # If the directory does not exist, create it\n",
    "    os.makedirs(directory)\n",
    "\n",
    "# Save the figure\n",
    "plt.savefig(os.path.join(directory, filename), bbox_inches='tight')\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "diamondenv",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}