4 add files for lct paper

2024_Ploetzke_et_al_Revisiting_Linear_Chain_Trick/CMakeLists.txt

@@ -27,6 +27,7 @@ FetchContent_MakeAvailable(memilio)
 set(MEMILIO_BUILD_TESTS OFF)
 set(MEMILIO_BUILD_EXAMPLES OFF)
 set(MEMILIO_BUILD_SIMULATIONS OFF)
+set(MEMILIO_ENABLE_OPENM OFF)
 
 # Add the subdirectory for MEmilio build.
 add_subdirectory(${memilio_SOURCE_DIR}/cpp ${memilio_BINARY_DIR})
@@ -36,5 +37,19 @@ 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)
+
+add_executable(lct_impact_age_resolution lct_impact_age_resolution.cpp)
+target_link_libraries(lct_impact_age_resolution PRIVATE memilio lct_secir)
+
+add_executable(lct_covid19_inspired_scenario lct_covid19_inspired_scenario.cpp)
+target_link_libraries(lct_covid19_inspired_scenario PRIVATE memilio lct_secir)
+
+# Just build run time measurement target if MEMILIO_ENABLE_OPENMP is explicitly set to on.
+if (MEMILIO_ENABLE_OPENMP)
+  add_executable(lct_runtime lct_runtime.cpp)
+  target_link_libraries(lct_runtime PRIVATE memilio lct_secir)
+  # target_compile_options(lct_runtime PRIVATE -O0)
+endif()

2024_Ploetzke_et_al_Revisiting_Linear_Chain_Trick/README.md

@@ -5,7 +5,26 @@ This directory contains all files related to the paper
 - _Lena Plötzke, Anna Wendler, René Schmieding, Martin J. Kühn. Revisiting the Linear Chain Trick in epidemiological models: Implications of underlying assumptions for numerical solutions,  (2024)._ 
 https://doi.org/10.48550/arXiv.2412.09140.
 
-Below is an overview of the files and the paper sections they belong to.
-
 ## Requirements
 We require the libraries `JsonCpp` and `HDF5` for running the simulations (these are optional for the MEmilio project, see [MEmilio cpp README](https://github.com/SciCompMod/memilio/blob/main/cpp/README.md)).
+
+The memilio.epidata package needs to be installed for the python plot scripts and the data download. 
+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.
+
+## Overview
+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 target `lct_runtime` is only built if the option `MEMILIO_ENABLE_OPENMP` is enabled. 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. Please execute the file [download_infection_data](download_infection_data.py) for this purpose.
+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. 

2024_Ploetzke_et_al_Revisiting_Linear_Chain_Trick/download_data.py

@@ -0,0 +1,74 @@
+#############################################################################
+# Copyright (C) 2020-2025 MEmilio
+#
+# Authors: Lena Ploetzke
+#
+# Contact: Martin J. Kuehn <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+############################################################################
+from datetime import date
+import pandas as pd
+import os
+from memilio.epidata import defaultDict as dd
+from memilio.epidata import getCaseData
+from memilio.epidata import getDIVIData
+
+
+def main():
+
+    data_folder= os.path.join(os.path.dirname(
+        __file__), "data")
+
+    if not os.path.isdir(data_folder):
+        os.makedirs(data_folder)
+
+    # Download data in format required.
+    getCaseData.get_case_data(read_data=False,
+                              file_format=dd.defaultDict['file_format'],
+                              out_folder=data_folder,
+                              no_raw=False,
+                              start_date=date(2020, 1, 1),
+                              end_date=date(2020, 12, 31),
+                              impute_dates=True,
+                              moving_average=0,
+                              make_plot=dd.defaultDict['make_plot'],
+                              split_berlin=dd.defaultDict['split_berlin'],
+                              rep_date=dd.defaultDict['rep_date'],
+                              files='All'
+                              )
+    getCaseData.get_case_data(read_data=False,
+                              file_format=dd.defaultDict['file_format'],
+                              out_folder=data_folder,
+                              no_raw=False,
+                              start_date=date(2020, 1, 1),
+                              end_date=date(2020, 12, 31),
+                              impute_dates=True,
+                              moving_average=7,
+                              make_plot=dd.defaultDict['make_plot'],
+                              split_berlin=dd.defaultDict['split_berlin'],
+                              rep_date=dd.defaultDict['rep_date'],
+                              files='All'
+                              )
+    getDIVIData.get_divi_data(read_data=False,
+                              file_format=dd.defaultDict['file_format'],
+                              out_folder=data_folder,
+                              start_date=date(2020, 1, 1),
+                              end_date=date(2020, 12, 31),
+                              impute_dates=True,
+                              moving_average=0
+                              )
+
+
+if __name__ == "__main__":
+    main()

2024_Ploetzke_et_al_Revisiting_Linear_Chain_Trick/get_data_covid19_inspired.sh

@@ -0,0 +1,48 @@
+#!/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.
+if [ ! -d "build/" ]; then
+    mkdir "build/"
+fi
+cd build/
+cmake ..
+
+infection_data_dir="../data"
+contact_data_dir="_deps/memilio-src/data/contacts/"
+result_dir="../simulation_results/"
+if [ ! -d "$result_dir" ]; then
+    mkdir "$result_dir"
+fi
+result_dir="../simulation_results/simulation_lct_covid19/"
+if [ ! -d "$result_dir" ]; then
+    mkdir "$result_dir"
+fi
+
+# 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 $infection_data_dir $contact_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 $infection_data_dir $contact_data_dir $result_dir $year $month_oct $day_oct $RelativeTransmissionNoSymptoms $RiskOfInfectionFromSymptomatic $scale_contacts_oct $scale_confirmed_cases_oct $npi_size_oct 

2024_Ploetzke_et_al_Revisiting_Linear_Chain_Trick/get_data_numerical_experiments.sh

@@ -0,0 +1,110 @@
+#!/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.
+if [ ! -d "build/" ]; then
+    mkdir "build/"
+fi
+cd build/
+cmake ..
+
+result_dir="../simulation_results/"
+if [ ! -d "$result_dir" ]; then
+    mkdir "$result_dir"
+fi
+result_dir="../simulation_results/simulation_lct_numerical_experiments/"
+if [ ! -d "$result_dir" ]; then
+    mkdir "$result_dir"
+fi
+
+
+subdir_dropReff="$result_dir/dropReff/"
+if [ ! -d "$subdir_dropReff" ]; then
+    mkdir "$subdir_dropReff"
+fi
+subdir_riseReffTo2short="$result_dir/riseReffTo2short/"
+if [ ! -d "$subdir_riseReffTo2short" ]; then
+    mkdir "$subdir_riseReffTo2short"
+fi
+subdir_riseReffTo2shortTEhalved="$result_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="$result_dir/riseRefflong/"
+if [ ! -d "$subdir_riseRefflong" ]; then
+    mkdir "$subdir_riseRefflong"
+fi
+subdir_riseRefflongTEhalved="$result_dir/riseRefflongTEhalved/"
+if [ ! -d "$subdir_riseRefflongTEhalved" ]; then
+    mkdir "$subdir_riseRefflongTEhalved"
+fi
+subdir_riseRefflongTEdoubled="$result_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="$result_dir/age_resolution/"
+if [ ! -d "$subdir_age_resolution" ]; then
+    mkdir "$subdir_age_resolution"
+fi
+cmake --build . --target lct_impact_age_resolution 
+contact_data_dir="_deps/memilio-src/data/contacts/"
+./bin/lct_impact_age_resolution $contact_data_dir $subdir_age_resolution 

2024_Ploetzke_et_al_Revisiting_Linear_Chain_Trick/get_runtimes_lct.sh

@@ -0,0 +1,33 @@
+#!/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.
+
+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

2024_Ploetzke_et_al_Revisiting_Linear_Chain_Trick/git_tag.cmake

@@ -1,3 +1,2 @@
 # GIT_TAG_MEMILIO can be set to a commit hash or a branch name.
-#set(GIT_TAG_MEMILIO "8958e3159c7cdb64265f93f21d025f76d33a91c0")
-set(GIT_TAG_MEMILIO "1121-lct-paper")
+set(GIT_TAG_MEMILIO "09ecf790c5bcef30dc95b21eee8cffa99f8bddd6")