Skip to content

Commit 104dc58

Browse files
annawendlerannawendler
committed
add exercises for lct and ide
1 parent 71ffcef commit 104dc58

4 files changed

Lines changed: 305 additions & 4 deletions

File tree

‎cpp-tutorials/CMakeLists.txt‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -39,10 +39,16 @@ set(MEMILIO_BUILD_SIMULATIONS OFF)
3939
add_subdirectory(${memilio_SOURCE_DIR}/cpp ${memilio_BINARY_DIR})
4040

4141
add_executable(tutorial_ide tutorial_ide.cpp)
42-
target_link_libraries(tutorial_ide PRIVATE memilio ide_secir Boost::filesystem)
42+
target_link_libraries(tutorial_ide PRIVATE memilio ide_secir)
43+
44+
add_executable(exercise_ide exercise_ide.cpp)
45+
target_link_libraries(exercise_ide PRIVATE memilio ide_secir)
4346

4447
add_executable(tutorial_lct tutorial_lct.cpp)
45-
target_link_libraries(tutorial_lct PRIVATE memilio lct_secir Boost::filesystem)
48+
target_link_libraries(tutorial_lct PRIVATE memilio lct_secir)
49+
50+
add_executable(exercise_lct exercise_lct.cpp)
51+
target_link_libraries(exercise_lct PRIVATE memilio lct_secir)
4652

4753
# add_executable(tutorial_abm tutorial_abm_household1.cpp)
4854
# target_link_libraries(tutorial_abm PRIVATE memilio abm Boost::filesystem)

‎cpp-tutorials/exercise_ide.cpp‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
#include "ide_secir/model.h"
2+
#include "ide_secir/infection_state.h"
3+
#include "ide_secir/simulation.h"
4+
#include "memilio/config.h"
5+
#include "memilio/epidemiology/age_group.h"
6+
#include "memilio/math/eigen.h"
7+
#include "memilio/utils/custom_index_array.h"
8+
#include "memilio/utils/time_series.h"
9+
#include "memilio/epidemiology/uncertain_matrix.h"
10+
#include "memilio/epidemiology/state_age_function.h"
11+
#include "memilio/data/analyze_result.h"
12+
13+
int main()
14+
{
15+
// MEmilio implements two models based on integro-differential equations (IDEs) with different infection states.
16+
// IDE-based models are a generalization of ODE-based models. Whereas ODE-based models assume an exponential
17+
// distribution regarding the time spent in an infection state, IDE-based models allow for arbitrary distributions.
18+
19+
// The following example shows how to set up and run a simple IDE-SECIR model without any further stratification.
20+
21+
using Vec = Eigen::VectorX<ScalarType>;
22+
23+
/*** Model setup ***/
24+
// Define simulation parameters.
25+
ScalarType t0 = 0.;
26+
ScalarType tmax = 20.;
27+
ScalarType dt = 0.01; // The step size will stay constant throughout the simulation.
28+
29+
// We define the number of age groups used in the model.
30+
size_t num_agegroups = 1;
31+
32+
// Next, we define initial values for the total population and number of deaths per age group.
33+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population_init =
34+
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 1000.);
35+
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths_init =
36+
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 0.);
37+
38+
// Now, we create a TimeSeries object with num_transitions * num_agegroups elements where initial transitions needed for simulation
39+
// will be stored. We require values for the transitions for a sufficient number of time points before the start of
40+
// the simulation to initialize our model.
41+
size_t num_transitions = (size_t)mio::isecir::InfectionTransition::Count;
42+
mio::TimeSeries<ScalarType> transitions_init(num_transitions);
43+
44+
// Then, we define vector of transitions that will be added as values to the time points of the TimeSeries transitions_init.
45+
Vec vec_init(num_transitions);
46+
vec_init[(size_t)mio::isecir::InfectionTransition::SusceptibleToExposed] = 3.0;
47+
vec_init[(size_t)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 0.0;
48+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 0.0;
49+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] = 0.0;
50+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] = 0.0;
51+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = 0.0;
52+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] = 0.0;
53+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 0.0;
54+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 0.0;
55+
vec_init[(size_t)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 0.0;
56+
57+
// We multiply vec_init with dt so that within a time interval of length 1, always the above number of
58+
// individuals are transitioning from one compartment to another, irrespective of the chosen time step size.
59+
vec_init = vec_init * dt;
60+
61+
// In this example, we will set the TransitionDistributions below. For these distributions, setting the initial time
62+
// point of the TimeSeries transitions_init at time -10 will give us a sufficient number of time points before t0.
63+
// For more information on this, we refer to the documentation of TransitionDistributions in
64+
// models/ide_secir/parameters.h.
65+
66+
// TODO: Set the initial time point of the TimeSeries transitions_init where the initial time is -10 and the value
67+
// is given by vec_init.
68+
69+
// We add further time points with distance dt until time t0. Here, we already define the start time of our simulation
70+
// as we will start the simulation at the last time point in transitions_init.
71+
while (transitions_init.get_last_time() < t0 - dt / 2) {
72+
transitions_init.add_time_point(transitions_init.get_last_time() + dt, vec_init);
73+
}
74+
75+
// With this, we can initialize the model.
76+
mio::isecir::Model model(std::move(transitions_init), total_population_init, deaths_init, num_agegroups);
77+
78+
// We set the number of initially Recovered individuals to 10.
79+
model.populations[0][(size_t)mio::isecir::InfectionState::Recovered] = 10.;
80+
81+
// After having initialized the model, we now set the epidemiological parameters.
82+
// We start with setting the transition distributions between the InfectionStates. With this, we define the times
83+
// individuals spend on average in the respective InfectionStates.
84+
// We start by defining a SmootherCosine object that is defined by an initial distribtion parameter.
85+
mio::SmootherCosine<ScalarType> smoothcos(4.0);
86+
// This is passed to a StateAgeFunctionWrapper object which is the type of oject used within the model to allow
87+
// for variable transition distributions.
88+
mio::StateAgeFunctionWrapper<ScalarType> transition_distribution(smoothcos);
89+
// We define a vector of StateAgeFunctionsWrappers where we set each transition distribution of our model to the
90+
// above defined transition_distribution.
91+
std::vector<mio::StateAgeFunctionWrapper<ScalarType>> vec_transition_distribution(num_transitions,
92+
transition_distribution);
93+
// Each transition can be set individually and is not necessary that all transition follow the same distribution
94+
// with the same parameter. For this, MEmilio already implements several possible distributions such as exponential,
95+
// gamma and lognormal distributions but arbitrary distributions can be implemented, see
96+
// cpp/memilio/epidemiology/state_age_function.h.
97+
//
98+
// TODO: Set the distribution for the transition from InfectedNoSymptoms to InfectedSymptoms with a lognormal
99+
// distribution with an initial distribution parameter of 0.3.
100+
101+
// Finally, the TransitionDistributions of the IDE-SECIR model are set according to vec_transition_distribution.
102+
model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(0)] = vec_transition_distribution;
103+
104+
// The following parameters define the relevant transition probabilities between InfectionStates.
105+
std::vector<ScalarType> vec_prob(num_transitions, 0.5);
106+
// The following probabilities must be 1, as there is no other way to go.
107+
vec_prob[(size_t)mio::isecir::InfectionTransition::SusceptibleToExposed] = 1;
108+
vec_prob[(size_t)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 1;
109+
// The TransitionProbabilities of the IDE-SECIR model are set according to vec_prob.
110+
model.parameters.get<mio::isecir::TransitionProbabilities>()[mio::AgeGroup(0)] = vec_prob;
111+
112+
// Further epidemiological parameters define the transmission probability of Susceptibles on contact with
113+
// infectious individuals, the relative risk of transmission from individuals that are infectious but not
114+
// symptomatic, and the risk of infection from indidivuals that are infectious and symptomatic. These parameters
115+
// can be set as functions of the time since infection and are set as the survival function of the exponential distribution here.
116+
mio::ExponentialSurvivalFunction<ScalarType> exponential(0.3);
117+
mio::StateAgeFunctionWrapper<ScalarType> prob(exponential);
118+
model.parameters.get<mio::isecir::TransmissionProbabilityOnContact>()[mio::AgeGroup(0)] = prob;
119+
model.parameters.get<mio::isecir::RelativeTransmissionNoSymptoms>()[mio::AgeGroup(0)] = prob;
120+
model.parameters.get<mio::isecir::RiskOfInfectionFromSymptomatic>()[mio::AgeGroup(0)] = prob;
121+
122+
// In order to include seasonality, we can set the seasonality's impact and the start day of
123+
// simulation. Seasonality is modeled by a sinoidal function. With a Seasonality value of 0.2, the risk of
124+
// transmission on January 1st is 50 % higher than on July 1st.
125+
model.parameters.get<mio::isecir::Seasonality>() = 0.1;
126+
model.parameters.get<mio::isecir::StartDay>() =
127+
40.; // Start the simulation on the 40th day of a year (i.e. February 9).
128+
129+
// Transmission is driven by the risk of transmission per contact and the contact matrix that defines the
130+
// average daily number of contacts between individuals. For a model with a single age group, the contact matrix
131+
// reduces to a simple scalar value (here, set to 10).
132+
mio::ContactMatrixGroup<ScalarType>& contact_matrix = model.parameters.get<mio::isecir::ContactPatterns>();
133+
contact_matrix[0] =
134+
mio::ContactMatrix<ScalarType>(Eigen::MatrixX<ScalarType>::Constant(num_agegroups, num_agegroups, 10.));
135+
136+
// We call the global_support_max method to get the maximum support_max over all TransitionDistributions. This
137+
// determines how many historic time points we need when initializing the model.
138+
std::cout << "Global support max: " << model.get_global_support_max(dt) << std::endl;
139+
140+
// *** Model simulation ***
141+
// We check if all model constraints regarding initial values and parameters are satisfied before simulating.
142+
// Note: MEmilio's check_constraints() returns True if a constraint is violated, and False if everything is fine.
143+
model.check_constraints(dt);
144+
145+
// We then perform a simulation.
146+
mio::isecir::Simulation sim(model, dt);
147+
sim.advance(tmax);
148+
149+
// We interpolate the simulation results for compartments and flows to days and print the results.
150+
auto interpolated_results = mio::interpolate_simulation_result(sim.get_result());
151+
auto interpolated_results_flows = mio::interpolate_simulation_result(sim.get_transitions());
152+
153+
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D"}, 12, 4);
154+
interpolated_results_flows.print_table(
155+
{"S->E", "E->C", "C->I", "C->R", "I->H", "I->R", "H->U", "H->R", "U->D", "U->R"}, 12, 4);
156+
157+
// We export the results as csv which is saved in the current folder. Then we can plot the results using plot_secir_results.py.
158+
auto export_status = sim.get_result().export_csv("../../cpp-tutorials/results_ide.csv");
159+
}

0 commit comments

Comments
 (0)