ABM_hybrid_data_creation10.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "000644ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Covasim 3.1.4 (2022-10-22) — © 2020-2022 by IDM\n"
     ]
    }
   ],
   "source": [
    "import sciris as sc\n",
    "import pandas as pd\n",
    "import covasim as cv\n",
    "import covasim.parameters as cvp\n",
    "import pylab as pl\n",
    "import numpy as np\n",
    "import os\n",
    "import math\n",
    "import random\n",
    "import dask\n",
    "from dask.distributed import Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c9da2418",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set working directory\n",
    "new_directory = 'C:/Users/irisg/Documents/PhD/COVID_France/SEIR_vs_Rt_sims/ABM_2params_all_at_once10'\n",
    "os.chdir(new_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "64a87b34",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_hybrid(params):\n",
    "    beta, pop, initE, index = params\n",
    "    pars = dict(\n",
    "        pop_size=pop, pop_infected=initE, beta=beta,\n",
    "        pop_type='hybrid', location='France',\n",
    "        start_day='2020-02-01', end_day='2020-06-30',\n",
    "        use_waning=False, prog_by_age = True,\n",
    "        pop_scale = 10, rescale = True, \n",
    "        asymp_factor = 1\n",
    "    )\n",
    "    \n",
    "    # transmission reduction of NPI 1: exp(-1.45) = 0.2345703\n",
    "    NPI1 = cv.change_beta(days='2020-03-17', changes=0.2345703, layers = ['c', 'w', 's'])\n",
    "    # transmission reduction of NPI 2: exp(-0.8) = 0.449329\n",
    "    NPI2 = cv.change_beta(days='2020-05-11', changes=0.449329, layers = ['c', 'w', 's'])\n",
    "    \n",
    "    # no testing structure\n",
    "    tests = cv.test_prob(symp_prob = 0, asymp_prob = 0, start_day = '2020-02-15')\n",
    "\n",
    "    # run simulation\n",
    "    sim = cv.Sim(pars, interventions=[NPI1, NPI2, tests], label= f'Sim {index}')\n",
    "    sim.run()\n",
    "\n",
    "    # return  observations and Rt\n",
    "    return sim.results['r_eff'].values[1:151], sim.results['new_infectious'].values[1:151], sim.results['new_severe'].values[1:151], sim.results['new_deaths'].values[1:151], sim.results['n_severe'].values[1:151]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4d7b02af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulation settings: population size, number of regions and time frame, number of repetitions\n",
    "popsize_df = pd.read_csv('C:/Users/irisg/Documents/PhD/COVID_France/Dropbox_iris_covid/departement/simulations/Rt_trajectories/popsize_df.csv')\n",
    "popsizes = popsize_df['popsize']\n",
    "ps = popsizes.values\n",
    "\n",
    "depts = np.arange(1, 95)\n",
    "days = np.arange(1, 151)\n",
    "random_seed = np.arange(100, 201)\n",
    "reps = np.arange(1, 101)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "97422b8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "\n",
    "    # Run settings\n",
    "    n_workers = 12\n",
    "    \n",
    "    # Create and queue the Dask jobs\n",
    "    client = Client(n_workers=n_workers)\n",
    "    \n",
    "    for rep in reps:\n",
    "        \n",
    "        hybrid_sims = []\n",
    "       \n",
    "        # randomly sample starting parameters\n",
    "        np.random.seed(random_seed[rep])\n",
    "        betas = np.random.lognormal(math.log(0.016), 0.04, size=94)\n",
    "        initEs = np.random.normal(50, 10, size=94)\n",
    "        \n",
    "        indeces = np.arange(1, 95)\n",
    "   \n",
    "        # set simulation for all departments\n",
    "        for beta, pop, initE, index in zip(betas, popsizes, initEs, indeces):\n",
    "            params_list = [beta, pop, initE, index]\n",
    "            sim = dask.delayed(simulate_hybrid)(params_list)\n",
    "            hybrid_sims.append(sim)\n",
    "\n",
    "        # Run and process the simulations\n",
    "        res = list(dask.compute(*hybrid_sims))\n",
    "        \n",
    "        # Gather the simulated time series into a dataframe \n",
    "        Rt_list = []\n",
    "        IncI_list = []\n",
    "        IncH_list = []\n",
    "        IncD_list = []\n",
    "        PrevH_list = []\n",
    "        for item in res:\n",
    "            Rt_list.extend(item[0])\n",
    "            IncI_list.extend(item[1])\n",
    "            IncH_list.extend(item[2])\n",
    "            IncD_list.extend(item[3])\n",
    "            PrevH_list.extend(item[4])\n",
    "\n",
    "        # Create a DataFrame\n",
    "        hybrid_Inc_df = pd.DataFrame({'Rt': Rt_list, 'IncI': IncI_list, 'IncH': IncH_list, 'IncD': IncD_list, 'PrevH': PrevH_list})\n",
    "\n",
    "        \n",
    "        # create other rows necessary for Monolix estimation \n",
    "        hybrid_Inc_df['dept_id'] = np.repeat(depts, 150)\n",
    "        hybrid_Inc_df['day'] = np.tile(days, 94)\n",
    "        hybrid_Inc_df['lockdown1'] = np.where((hybrid_Inc_df['day'].between(45, 99)), 1, 0)\n",
    "        hybrid_Inc_df['BG1'] = np.where(hybrid_Inc_df['day'] >= 100, 1, 0)\n",
    "        hybrid_Inc_df['popsize'] = np.repeat(ps, 150)\n",
    "        \n",
    "        hybrid_Inc_df = hybrid_Inc_df[['dept_id', 'day', 'IncI', 'IncH', 'IncD', 'PrevH', 'Rt', 'lockdown1', 'BG1', 'popsize']]\n",
    "\n",
    "        # save data set\n",
    "        hybrid_Inc_df.to_csv(f'data_covasim_hybrid10_Rt_{rep}.csv', index=False, sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "d5faf432",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff212aaf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}