ABM_rm_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_rm(params):\n",
    "    beta, pop, initE, index = params\n",
    "    pars = dict(\n",
    "        pop_size=pop, pop_infected=initE, beta=beta,\n",
    "        pop_type='random', location='France',\n",
    "        start_day='2020-02-01', end_day='2020-06-30',\n",
    "        use_waning=False, prog_by_age = False,\n",
    "        pop_scale = 10, rescale = True, \n",
    "        asymp_factor = 1,\n",
    "        prognoses = dict(age_cutoffs = np.array([0, 50]),\n",
    "                         sus_ORs = np.array([1, 1]),\n",
    "                         trans_ORs = np.array([1, 1]),\n",
    "                         comorbidities = np.array([1, 1]),\n",
    "                         symp_probs = np.array([1, 1]),\n",
    "                         severe_probs = np.array([0.1, 0.1]),\n",
    "                         crit_probs = np.array([1, 1]),\n",
    "                         death_probs = np.array([0.1, 0.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)\n",
    "    # transmission reduction of NPI 2: exp(-0.8) = 0.449329\n",
    "    NPI2 = cv.change_beta(days='2020-05-11', changes=0.449329)\n",
    "    \n",
    "    # no testing structure\n",
    "    tests = cv.test_prob(symp_prob = 0, asymp_prob = 0, start_day = '2020-02-01')\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": 14,
   "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": 15,
   "id": "97422b8f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\irisg\\anaconda3\\lib\\site-packages\\distributed\\node.py:179: UserWarning: Port 8787 is already in use.\n",
      "Perhaps you already have a cluster running?\n",
      "Hosting the HTTP server on port 52591 instead\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "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",
    "        rm_sims = []\n",
    "       \n",
    "        # randomly sample starting parameters\n",
    "        np.random.seed(random_seed[rep])\n",
    "        betas = np.random.lognormal(math.log(0.014), 0.04, size=94)\n",
    "        initEs = np.random.normal(50, 10, size=94)\n",
    "        \n",
    "        index = np.arange(1, 95)\n",
    "   \n",
    "        # set simulation for all departments\n",
    "        for beta, pop, initE, i in zip(betas, popsizes, initEs, index):\n",
    "            params_list = [beta, pop, initE, i]\n",
    "            sim = dask.delayed(simulate_rm)(params_list)\n",
    "            rm_sims.append(sim)\n",
    "\n",
    "        # Run and process the simulations\n",
    "        res = list(dask.compute(*rm_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",
    "        rm_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 and scale observations to 10,000 population\n",
    "        rm_Inc_df['dept_id'] = np.repeat(depts, 150)\n",
    "        rm_Inc_df['day'] = np.tile(days, 94)\n",
    "        rm_Inc_df['lockdown1'] = np.where((rm_Inc_df['day'].between(45, 99)), 1, 0)\n",
    "        rm_Inc_df['BG1'] = np.where(rm_Inc_df['day'] >= 100, 1, 0)\n",
    "        rm_Inc_df['popsize'] = np.repeat(ps, 150)\n",
    "        \n",
    "        rm_Inc_df = rm_Inc_df[['dept_id', 'day', 'IncI', 'IncH', 'IncD', 'PrevH', 'Rt', 'lockdown1', 'BG1', 'popsize']]\n",
    "\n",
    "        # save data set\n",
    "        rm_Inc_df.to_csv(f'data_covasim_rm10_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
}