1121 LCT paper

cpp/CMakeLists.txt

@@ -174,6 +174,7 @@ endif()
 
 if(MEMILIO_BUILD_SIMULATIONS)
     add_subdirectory(simulations)
+    add_subdirectory(simulations/2024_Ploetzke_Linear_Chain_Trick)
 endif()
 
 if(MEMILIO_BUILD_BENCHMARKS)

cpp/simulations/2024_Ploetzke_Linear_Chain_Trick/CMakeLists.txt

@@ -0,0 +1,35 @@
+if(MEMILIO_HAS_HDF5)
+    if(NOT DEFINED NUM_SUBCOMPARTMENTS)
+            set(NUM_SUBCOMPARTMENTS "42")
+        endif()
+    add_definitions(-DNUM_SUBCOMPARTMENTS=${NUM_SUBCOMPARTMENTS})
+    add_executable(lct_impact_distribution_assumption lct_impact_distribution_assumption.cpp)
+    target_link_libraries(lct_impact_distribution_assumption PRIVATE memilio lct_secir)
+    target_compile_options(lct_impact_distribution_assumption PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+endif()
+
+if(MEMILIO_HAS_JSONCPP AND MEMILIO_HAS_HDF5)
+    add_executable(lct_impact_age_resolution lct_impact_age_resolution.cpp)
+    target_link_libraries(lct_impact_age_resolution PRIVATE memilio lct_secir Boost::filesystem)
+    target_compile_options(lct_impact_age_resolution PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+    if(NOT DEFINED NUM_SUBCOMPARTMENTS)
+            set(NUM_SUBCOMPARTMENTS "42")
+    endif()
+    add_definitions(-DNUM_SUBCOMPARTMENTS=${NUM_SUBCOMPARTMENTS})
+    add_executable(lct_covid19_inspired_scenario lct_covid19_inspired_scenario.cpp)
+    target_link_libraries(lct_covid19_inspired_scenario PRIVATE memilio lct_secir Boost::filesystem)
+    target_compile_options(lct_covid19_inspired_scenario PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+endif()
+
+if (MEMILIO_ENABLE_OPENMP)
+    if(NOT DEFINED NUM_SUBCOMPARTMENTS)
+            set(NUM_SUBCOMPARTMENTS "42")
+        endif()
+    add_definitions(-DNUM_SUBCOMPARTMENTS=${NUM_SUBCOMPARTMENTS})
+
+    add_executable(lct_runtime lct_runtime.cpp)
+    target_link_libraries(lct_runtime PRIVATE memilio lct_secir)
+    target_compile_options(lct_runtime PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+    # target_compile_options(lct_runtime PRIVATE -O0)
+endif()

cpp/simulations/2024_Ploetzke_Linear_Chain_Trick/README.md

@@ -0,0 +1,35 @@
+# Revisiting the Linear Chain Trick in epidemiological models: Implications of underlying assumptions for numerical solutions #
+
+In this directory you will find all files related to the Paper 
+
+- _Lena Plötzke , Anna Wendler , René Schmieding , Martin J. Kühn (2024). Revisiting the Linear Chain Trick in epidemiological models: Implications of underlying assumptions for numerical solutions._ 
+https://doi.org/10.48550/arXiv.2412.09140.
+
+Below is an overview of the files and the paper sections they belong to.
+
+- Section 4.2: The file [lct_impact_distribution_assumption](lct_impact_distribution_assumption.cpp) provides the functionality to run simulations to assess the impact of the distribution assumption or different numbers of subcompartments. The dynamics at change points and epidemic peaks can be examined using this script. The population is not divided into age groups for these experiments. All simulation results are created and saved in `.h5` files when the shell script [get_data_numerical_experiments](get_data_numerical_experiments.sh) is executed. The visualizations of these simulation results in the paper were created using the python script [plot_numerical_experiments](plot_numerical_experiments.py).
+
+- Section 4.3: With the file [lct_impact_age_resolution](lct_impact_age_resolution.cpp), one can run simulations to assess the impact of including an age resolution. The simulation results are created and saved together with the results for section 4.2 with the shellscript [get_data_numerical_experiments](get_data_numerical_experiments.sh). The visualizations are also created with [plot_numerical_experiments](plot_numerical_experiments.py).
+
+- Section 4.4: Run time measurements are possible with the file [lct_runtime](lct_runtime.cpp). `OpenMP` is used to measure the run times. The Slurm script [get_runtimes_lct](get_runtimes_lct.sh) can be used to define a job to measure the run time for different numbers of subcompartments. This script can be easily adapted to use an adaptive solver. To use the optimization flag `-O0`, uncomment the suitable line in the [CMakeLists file](CMakeLists.txt). Visualizations of the run times in the paper were created using the python script [plot_runtimes_lct](plot_runtimes_lct.py).
+
+- Section 4.5: A COVID-19 inspired scenario in Germany in 2020 is defined in the file [lct_covid19_inspired_scenario](lct_covid19_inspired_scenario.cpp). The simulation results are created and saved with the shell script [get_data_covid19_inspired](get_data_covid19_inspired.sh). 
+    The simulation is initialized using real data, which has to be downloaded beforehand:
+    1. Reported case data from the RKI can be downloaded via [getCaseData](../../../pycode/memilio-epidata/memilio/epidata/getCaseData.py).
+
+    2. Reported number of patients in intensive care units (DIVI data) can be downloaded via [getDIVIData](../../../pycode/memilio-epidata/memilio/epidata/getDIVIData.py).
+
+    - The option impute_dates should be set to True and moving_average=0 while downloading both datasets. For the other function parameters, one can use default parameters from [defaultDict](../../../pycode/memilio-epidata/memilio/epidata/defaultDict.py) (and matching dates).
+
+    The visualizations of the simulation results in the paper were created using the python script [plot_covid19_inspired](plot_covid19_inspired.py).
+
+- Figure 2 and Figure 12: These figures are not based on simulation results. Figure 2 contains a visualization of the density and the survival function of Erlang distributions with different parameter choices. Figure 12 shows the age-resolved contact pattern for Germany. Both plots are created using [plot_details](plot_details.py).
+
+
+For most of the above `.cpp` files, the number of subcompartments used in the LCT models for all compartments can be controlled via the preprocessor macro NUM_SUBCOMPARTMENTS. Have a look at the files for further documentation or the shell scripts for the usage. 
+
+# Requirements
+We require the libraries `JsonCpp` and `HDF5` for running the scripts (these are optional for the whole project, see [README](../../README.md)).
+
+The memilio.epidata package needs to be installed for the python plot scripts. 
+Have a look at the [pycode README](../../../pycode/README.rst) and the [memilio-epidata README](../../../pycode/memilio-epidata/README.rst) for instructions how to install the package.

cpp/simulations/2024_Ploetzke_Linear_Chain_Trick/get_data_covid19_inspired.sh

@@ -0,0 +1,43 @@
+#!/bin/bash
+
+## This script can be used to get all simulation data using the lct_covid19_inspired_scenario.cpp file.
+## It is necessary to download RKI and DIVI data beforehand. For more information see the README.
+
+# Define and construct relevant folders.
+cd ../../
+if [ ! -d "build/" ]; then
+    mkdir "build/"
+fi
+cd build/
+data_dir="../../data"
+result_dir="$data_dir/simulation_lct_covid19/"
+if [ ! -d "$result_dir" ]; then
+    mkdir "$result_dir"
+fi
+cmake ..
+
+# Define parameters used as command line arguments.
+year=2020
+RelativeTransmissionNoSymptoms=1.
+RiskOfInfectionFromSymptomatic=0.3
+month_oct=10
+day_oct=1
+scale_contacts_oct=0.6537
+npi_size_oct=0.3
+scale_confirmed_cases_oct=1.2
+
+# Compile with different numbers of subcompartments and run simulations.
+for num_subcomp in 1 3 10 50
+do
+    cmake -DNUM_SUBCOMPARTMENTS=$num_subcomp -DCMAKE_BUILD_TYPE="Release" .
+    cmake --build . --target lct_covid19_inspired_scenario
+
+    # Simulation for 01/10/2020.
+    ./bin/lct_covid19_inspired_scenario $data_dir $result_dir $year $month_oct $day_oct $RelativeTransmissionNoSymptoms $RiskOfInfectionFromSymptomatic $scale_contacts_oct $scale_confirmed_cases_oct $npi_size_oct 
+done
+
+# Setup with numbers of subcompartments so that each corresponds to the approximate stay time in the compartment.
+# This is done by setting the makro NUM_SUBCOMPARTMENTS to zero.
+cmake -DNUM_SUBCOMPARTMENTS=0 -DCMAKE_BUILD_TYPE="Release" .
+cmake --build . --target lct_covid19_inspired_scenario
+./bin/lct_covid19_inspired_scenario $data_dir $result_dir $year $month_oct $day_oct $RelativeTransmissionNoSymptoms $RiskOfInfectionFromSymptomatic $scale_contacts_oct $scale_confirmed_cases_oct $npi_size_oct 

cpp/simulations/2024_Ploetzke_Linear_Chain_Trick/get_data_numerical_experiments.sh

@@ -0,0 +1,106 @@
+#!/bin/bash
+
+## This script can be used to get all simulation data for the numerical experiments regarding the impact of the 
+## distribution assumption and of the age resolution.
+## The files lct_impact_distribution_assumption.cpp and lct_impact_age_resolution.cpp are used.
+
+# Define and construct relevant folders.
+cd ../../
+if [ ! -d "build/" ]; then
+    mkdir "build/"
+fi
+cd build/
+cmake ..
+
+dir="../../data/simulation_lct_numerical_experiments"
+if [ ! -d "$dir" ]; then
+    mkdir "$dir"
+fi
+
+subdir_dropReff="$dir/dropReff/"
+if [ ! -d "$subdir_dropReff" ]; then
+    mkdir "$subdir_dropReff"
+fi
+subdir_riseReffTo2short="$dir/riseReffTo2short/"
+if [ ! -d "$subdir_riseReffTo2short" ]; then
+    mkdir "$subdir_riseReffTo2short"
+fi
+subdir_riseReffTo2shortTEhalved="$dir/riseReffTo2shortTEhalved/"
+if [ ! -d "$subdir_riseReffTo2shortTEhalved" ]; then
+    mkdir "$subdir_riseReffTo2shortTEhalved"
+fi
+
+# Compile with the different numbers of subcompartments and run with different setups.
+for num_subcomp in 1 3 10 50
+do
+    cmake -DNUM_SUBCOMPARTMENTS=$num_subcomp -DCMAKE_BUILD_TYPE="Release" .
+    cmake --build . --target lct_impact_distribution_assumption
+
+    # First case: Decrease the effective reproduction number at simulation day 2 to 0.5 and simulate for 12 days.
+    Reff2=0.5
+    simulation_days=12
+    ./bin/lct_impact_distribution_assumption $Reff2 $simulation_days $subdir_dropReff
+
+    # Second case: Increase the effective reproduction number at simulation day 2 to 2 and simulate for 12 days.
+    Reff2=2.
+    simulation_days=12
+    # Additionally save result with division in subcompartments.
+    if [ "$num_subcomp" -eq 10 ] || [ "$num_subcomp" -eq 50 ]; then
+        ./bin/lct_impact_distribution_assumption $Reff2 $simulation_days $subdir_riseReffTo2short 1
+    fi
+    ./bin/lct_impact_distribution_assumption $Reff2 $simulation_days $subdir_riseReffTo2short
+
+    # Third case: Second case but with TimeExposed scaled by 0.5.
+    scale_TimeExposed=0.5
+    # Additionally save result with division in subcompartments.
+    if [ "$num_subcomp" -eq 50 ]; then
+        ./bin/lct_impact_distribution_assumption $Reff2 $simulation_days $subdir_riseReffTo2shortTEhalved 1 $scale_TimeExposed
+    fi
+    ./bin/lct_impact_distribution_assumption $Reff2 $simulation_days $subdir_riseReffTo2shortTEhalved 0 $scale_TimeExposed
+
+    # Fourth case: Print final sizes without saving results.
+    simulation_days=500
+    ./bin/lct_impact_distribution_assumption 2 $simulation_days "" 0 1.0 1
+    ./bin/lct_impact_distribution_assumption 4 $simulation_days "" 0 1.0 1
+    ./bin/lct_impact_distribution_assumption 10 $simulation_days "" 0 1.0 1
+done
+
+
+# Fifth case: Increase the effective reproduction number at simulation day 2 to different values and 
+# simulate for 200 days to compare epidemic peaks.
+# Also perform simulations with TimeExposed scaled by 0.5 or 2.
+# Define and construct relevant folders.
+subdir_riseRefflong="$dir/riseRefflong/"
+if [ ! -d "$subdir_riseRefflong" ]; then
+    mkdir "$subdir_riseRefflong"
+fi
+subdir_riseRefflongTEhalved="$dir/riseRefflongTEhalved/"
+if [ ! -d "$subdir_riseRefflongTEhalved" ]; then
+    mkdir "$subdir_riseRefflongTEhalved"
+fi
+subdir_riseRefflongTEdoubled="$dir/riseRefflongTEdoubled/"
+if [ ! -d "$subdir_riseRefflongTEdoubled" ]; then
+    mkdir "$subdir_riseRefflongTEdoubled"
+fi
+simulationdays=200
+Reff2s=(2 3 4 5 6 7 8 9 10)
+for num_subcomp in 1 2 3 4 5 6 7 8 9 10 50
+do
+    cmake -DNUM_SUBCOMPARTMENTS=$num_subcomp -DCMAKE_BUILD_TYPE="Release" .
+    cmake --build . --target lct_impact_distribution_assumption
+    for index in {0..8}
+    do
+        ./bin/lct_impact_distribution_assumption ${Reff2s[index]} ${simulationdays} $subdir_riseRefflong
+        ./bin/lct_impact_distribution_assumption ${Reff2s[index]} ${simulationdays} $subdir_riseRefflongTEhalved 0 0.5
+        ./bin/lct_impact_distribution_assumption ${Reff2s[index]} ${simulationdays} $subdir_riseRefflongTEdoubled 0 2.0
+    done
+done
+
+# Sixth case: Simulation for the impact of age resolution.
+subdir_age_resolution="$dir/age_resolution/"
+if [ ! -d "$subdir_age_resolution" ]; then
+    mkdir "$subdir_age_resolution"
+fi
+cmake --build . --target lct_impact_age_resolution 
+contact_data_dir="../../data/contacts/"
+./bin/lct_impact_age_resolution $contact_data_dir $subdir_age_resolution 

cpp/simulations/2024_Ploetzke_Linear_Chain_Trick/get_runtimes_lct.sh

@@ -0,0 +1,34 @@
+#!/bin/bash
+#SBATCH --job-name=lct-performance
+#SBATCH --output=lct-%A.out
+#SBATCH --error=lct-%A.err
+#SBATCH --nodes=1
+#SBATCH --ntasks=1
+#SBATCH --cpus-per-task=1
+#SBATCH --exclusive
+#SBATCH --exclude="be-cpu05, be-gpu01"
+#SBATCH --time=5-0:00:00
+
+## This script can be used to monitor runtimes for the lct model using the file lct_runtime.cpp.
+## The command "sbatch get_runtimes_lct.sh ./bin/lct_runtime" should be run in this folder 
+## to have valid folder instructions.
+
+num_runs=100
+num_warm_up_runs=10
+# Use 1 to measure run times for an adaptive solver, 0 for fixed step sizes.
+use_adaptive_solver=0
+echo Running $1 on node $SLURM_JOB_NODELIST with $warm_up_runs warm up runs and $runs runs.
+
+cd ../../
+if [ ! -d "build/" ]; then
+    mkdir "build/"
+fi
+cd build/
+cmake -DCMAKE_BUILD_TYPE="Release" -DMEMILIO_ENABLE_OPENMP=ON ..
+
+for i in {1..200}
+do  
+    cmake -DNUM_SUBCOMPARTMENTS=$i -DCMAKE_BUILD_TYPE="Release" -DMEMILIO_ENABLE_OPENMP=ON .
+    cmake --build . --target lct_runtime
+    srun --cpus-per-task=1 --cpu-bind=cores ./$1 $num_runs $num_warm_up_runs $use_adaptive_solver
+done