Skip to content

1392 update IDE-SECIR examples and improve some documentation #1402

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
annawendler merged 3 commits into from
Nov 4, 2025
Merged

1392 update IDE-SECIR examples and improve some documentation #1402

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ed40d2b
update IDE examples and improve some documentation
annawendler Oct 24, 2025
1d9a23d
small improvement
annawendler Oct 24, 2025
2352ea2
Apply suggestions from code review
annawendler Nov 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 0 additions & 4 deletions cpp/examples/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,10 +120,6 @@ add_executable(ide_secir_example ide_secir.cpp)
target_link_libraries(ide_secir_example PRIVATE memilio ide_secir)
target_compile_options(ide_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

add_executable(ide_secir_ageres_example ide_secir_ageres.cpp)
target_link_libraries(ide_secir_ageres_example PRIVATE memilio ide_secir)
target_compile_options(ide_secir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

add_executable(ide_secir_graph_example ide_secir_graph.cpp)
target_link_libraries(ide_secir_graph_example PRIVATE memilio ide_secir)
target_compile_options(ide_secir_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
Expand Down
44 changes: 27 additions & 17 deletions cpp/examples/ide_initialization.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,41 +54,51 @@ std::string setup(int argc, char** argv)

int main(int argc, char** argv)
{
// This is a simple example to demonstrate how to set initial data for the IDE-SECIR model using real data.
// This is a simple example to demonstrate how to set initial data for the IDE-SECIR model using reported data.
// A default initialization is used if no filename is provided in the command line.
// Have a look at the documentation of the set_initial_flows() function in models/ide_secir/parameters_io.h for a
// description of how to download suitable data.
// A valid filename could be for example "../../data/pydata/Germany/cases_all_germany_ma7.json" if the functionality to download real data is used.
// The default parameters of the IDE-SECIR model are used, so that the simulation results are not realistic and are for demonstration purpose only.
// A valid filename could be for example "../../data/pydata/Germany/cases_all_germany_ma7.json" if the functionality
// to download reported data is used.
// The default parameters of the IDE-SECIR model are used, so note that the simulation results are not realistic
// and are for demonstration purpose only.

// Define simulation parameters.
ScalarType t0 = 0.;
ScalarType tmax = 2.;
ScalarType dt = 0.01;

// Initialize model.
size_t num_agegroups = 1;
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population =
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population_init =
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 80 * 1e6);
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths = mio::CustomIndexArray<ScalarType, mio::AgeGroup>(
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths_init = mio::CustomIndexArray<ScalarType, mio::AgeGroup>(
mio::AgeGroup(num_agegroups),
0.); // The number of deaths will be overwritten if real data is used for initialization.
ScalarType dt = 0.5;
mio::isecir::Model model(mio::TimeSeries<ScalarType>((int)mio::isecir::InfectionTransition::Count),
total_population, deaths, num_agegroups);
0.); // The number of deaths will be overwritten if reported data is used for initialization.

mio::isecir::Model model(mio::TimeSeries<ScalarType>((size_t)mio::isecir::InfectionTransition::Count),
total_population_init, deaths_init, num_agegroups);

// Check provided parameters.
std::string filename = setup(argc, argv);
if (filename.empty()) {
std::cout << "You did not provide a valid filename. A default initialization is used." << std::endl;

// Create TimeSeries with num_transitions * num_agegroups elements where initial transitions needed for simulation
// will be stored. We require values for the transitions for a sufficient number of time points before the start of
// the simulation to initialize our model.
using Vec = Eigen::VectorX<ScalarType>;
mio::TimeSeries<ScalarType> init(num_agegroups * (size_t)mio::isecir::InfectionTransition::Count);
init.add_time_point<Eigen::VectorX<ScalarType>>(
mio::TimeSeries<ScalarType> transitions_init((size_t)mio::isecir::InfectionTransition::Count * num_agegroups);
transitions_init.add_time_point<Eigen::VectorX<ScalarType>>(
-7., Vec::Constant((size_t)mio::isecir::InfectionTransition::Count, 1. * dt));
while (init.get_last_time() < -dt / 2.) {
init.add_time_point(init.get_last_time() + dt,
Vec::Constant((int)mio::isecir::InfectionTransition::Count, 1. * dt));
while (transitions_init.get_last_time() < t0 - dt / 2.) {
transitions_init.add_time_point(transitions_init.get_last_time() + dt,
Vec::Constant((size_t)mio::isecir::InfectionTransition::Count, 1. * dt));
}
model.transitions = init;
model.transitions = transitions_init;
}
else {
// Use the real data for initialization.
// Use reported data for initialization.
// Here we assume that the file contains data without age resolution, hence we use read_confirmed_cases_noage()
// for reading the data and mio::ConfirmedCasesNoAgeEntry as EntryType in set_initial_flows().

Expand All @@ -112,7 +122,7 @@ int main(int argc, char** argv)

// Carry out simulation.
mio::isecir::Simulation sim(model, dt);
sim.advance(2.);
sim.advance(tmax);

// Print results.
sim.get_transitions().print_table({"S->E", "E->C", "C->I", "C->R", "I->H", "I->R", "H->U", "H->R", "U->D", "U->R"},
Expand Down
157 changes: 99 additions & 58 deletions cpp/examples/ide_secir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
/*
* Copyright (C) 2020-2025 MEmilio
*
* Authors: Anna Wendler, Lena Ploetzke
* Authors: Anna Wendler, Lena Ploetzke, Hannah Tritzschak
*
* Contact: Martin J. Kuehn <[email protected]>
*
Expand All @@ -17,7 +17,6 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/

#include "ide_secir/model.h"
#include "ide_secir/infection_state.h"
#include "ide_secir/simulation.h"
Expand All @@ -32,96 +31,138 @@

int main()
{
// This is a simple example to demonstrate how use the IDE-SECIR model.

using Vec = Eigen::VectorX<ScalarType>;

size_t num_agegroups = 1;
// Define simulation parameters.
ScalarType t0 = 0.;
ScalarType tmax = 5.;
ScalarType dt = 0.01; // The time step size will stay constant throughout the simulation.

ScalarType tmax = 10;
mio::CustomIndexArray<ScalarType, mio::AgeGroup> N =
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10000.);
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths =
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 13.10462213);
ScalarType dt = 1.;
// Define number of age groups.
size_t num_agegroups = 2;

int num_transitions = (int)mio::isecir::InfectionTransition::Count;
// Define initial values for the total population and number of deaths per age group.
mio::CustomIndexArray<ScalarType, mio::AgeGroup> total_population_init =
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 1000.);
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths_init =
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 6.);

// Create TimeSeries with num_transitions * num_agegroups elements where transitions needed for simulation will be
// stored.
mio::TimeSeries<ScalarType> init(num_transitions * num_agegroups);
// Create TimeSeries with num_transitions * num_agegroups elements where initial transitions needed for simulation
// will be stored. We require values for the transitions for a sufficient number of time points before the start of
// the simulation to initialize our model.
size_t num_transitions = (size_t)mio::isecir::InfectionTransition::Count;
mio::TimeSeries<ScalarType> transitions_init(num_transitions * num_agegroups);

// Add time points for initialization of transitions.
// Define vector of transitions that will be added to the time points of the TimeSeries transitions_init.
Vec vec_init(num_transitions * num_agegroups);
vec_init[(int)mio::isecir::InfectionTransition::SusceptibleToExposed] = 25.0;
vec_init[(int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 15.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] = 4.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] = 1.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = 4.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] = 1.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 1.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 1.0;
vec_init[(int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0;

for (size_t group = 0; group < num_agegroups; ++group) {
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::SusceptibleToExposed] = 25.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] =
15.0;
vec_init[group * num_transitions +
(size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] =
4.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] =
1.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = 4.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] =
1.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 1.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 1.0;
vec_init[group * num_transitions + (size_t)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0;
}
// Multiply vec_init with dt so that within a time interval of length 1, always the above number of
// individuals are transitioning from one compartment to another, irrespective of the chosen time step size.
vec_init = vec_init * dt;
// Add initial time point to time series.
init.add_time_point(-10, vec_init);
// Add further time points until time 0.
while (init.get_last_time() < -dt / 2) {
init.add_time_point(init.get_last_time() + dt, vec_init);

// In this example, we will set the TransitionDistributions below. For these distributions, setting the initial time
// point of the TimeSeries transitions_init at time -10 will give us a sufficient number of time points before t0=0.
// For more information on this, we refer to the documentation of TransitionDistributions in
// models/ide_secir/parameters.h.
transitions_init.add_time_point(-10, vec_init);
// Add further time points with distance dt until time t0.
while (transitions_init.get_last_time() < t0 - dt / 2) {
transitions_init.add_time_point(transitions_init.get_last_time() + dt, vec_init);
}

// Initialize model.
mio::isecir::Model model(std::move(init), N, deaths, num_agegroups);
mio::isecir::Model model(std::move(transitions_init), total_population_init, deaths_init, num_agegroups);

// Uncomment one of the code blocks below to use a different method to initialize the model, based on a
// given number of either Susceptibles or Recovered instead of using the TimeSeries transitions_init from above.

// Initialization method with given Susceptibles.
// size_t num_infstates = (size_t)mio::isecir::InfectionState::Count;
// for (size_t group = 0; group < num_agegroups; ++group) {
// model.populations.get_last_value()[group * num_infstates + (size_t)mio::isecir::InfectionState::Susceptible] =
// 900;
// }

// Uncomment one of the two lines to use a different method to initialize the model using the TimeSeries init.
// Initialization method with Susceptibles.
// model.populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000;
// Initialization method with Recovered.
// model.populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0;
// Initialization method with given Recovered.
// size_t num_infstates = (size_t)mio::isecir::InfectionState::Count;
// for (size_t group = 0; group < num_agegroups; ++group) {
// model.populations.get_last_value()[group * num_infstates + (size_t)mio::isecir::InfectionState::Recovered] = 10;
// }

// Set working parameters.
mio::SmootherCosine<ScalarType> smoothcos(2.0);
mio::StateAgeFunctionWrapper<ScalarType> delaydistribution(smoothcos);
std::vector<mio::StateAgeFunctionWrapper<ScalarType>> vec_delaydistrib(num_transitions, delaydistribution);
// TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value.
vec_delaydistrib[(int)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.);
vec_delaydistrib[(int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
.set_distribution_parameter(4.0);

model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(0)] = vec_delaydistrib;
// TransitionDistributions
// In the following, we explicitly set the TransitionDistributions for the first age group. If the model contains
// more age groups, the default distributions are used for these age groups.
mio::SmootherCosine<ScalarType> smoothcos1(3.0);
mio::StateAgeFunctionWrapper<ScalarType> delaydistribution1(smoothcos1);
std::vector<mio::StateAgeFunctionWrapper<ScalarType>> vec_delaydistrib1(num_transitions, delaydistribution1);
// TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value.
vec_delaydistrib1[(size_t)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.);
model.parameters.get<mio::isecir::TransitionDistributions>()[mio::AgeGroup(0)] = vec_delaydistrib1;

// TransitionProbabilities
std::vector<ScalarType> vec_prob(num_transitions, 0.5);
// The following probabilities must be 1, as there is no other way to go.
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] = 1;
vec_prob[Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] = 1;
model.parameters.get<mio::isecir::TransitionProbabilities>()[mio::AgeGroup(0)] = vec_prob;
vec_prob[(size_t)mio::isecir::InfectionTransition::SusceptibleToExposed] = 1;
vec_prob[(size_t)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 1;
for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) {
model.parameters.get<mio::isecir::TransitionProbabilities>()[group] = vec_prob;
}

// Contact patterns
mio::ContactMatrixGroup<ScalarType> contact_matrix = mio::ContactMatrixGroup<ScalarType>(1, num_agegroups);
contact_matrix[0] =
mio::ContactMatrix<ScalarType>(Eigen::MatrixX<ScalarType>::Constant(num_agegroups, num_agegroups, 10.));
model.parameters.get<mio::isecir::ContactPatterns>() = mio::UncertainContactMatrix<ScalarType>(contact_matrix);
model.parameters.get<mio::isecir::ContactPatterns>() = mio::UncertainContactMatrix(contact_matrix);

// Furhter epidemiological parameters
mio::ExponentialSurvivalFunction<ScalarType> exponential(0.5);
mio::StateAgeFunctionWrapper<ScalarType> prob(exponential);

model.parameters.get<mio::isecir::TransmissionProbabilityOnContact>()[mio::AgeGroup(0)] = prob;
model.parameters.get<mio::isecir::RelativeTransmissionNoSymptoms>()[mio::AgeGroup(0)] = prob;
model.parameters.get<mio::isecir::RiskOfInfectionFromSymptomatic>()[mio::AgeGroup(0)] = prob;

for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) {
model.parameters.get<mio::isecir::TransmissionProbabilityOnContact>()[group] = prob;
model.parameters.get<mio::isecir::RelativeTransmissionNoSymptoms>()[group] = prob;
model.parameters.get<mio::isecir::RiskOfInfectionFromSymptomatic>()[group] = prob;
}
model.parameters.set<mio::isecir::Seasonality>(0.1);
// Start the simulation on the 40th day of a year (i.e. in February).
model.parameters.set<mio::isecir::StartDay>(40);
model.parameters.set<mio::isecir::StartDay>(
40); // Start the simulation on the 40th day of a year (i.e. in February).

// Check if all model constraints regarding initial values and parameters are satisfied before simulating.
model.check_constraints(dt);

// Carry out simulation.
mio::isecir::Simulation sim(model, dt);
sim.advance(tmax);

// Interpolate results to days.
auto interpolated_results = mio::interpolate_simulation_result(sim.get_result(), dt / 2.);

interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
// Print results. Note that the column labels are suitable for a simulation with two age groups and may need to be
// adapted when the number of age groups is changed.
// interpolated_results.print_table(
// {"S1", "E1", "C1", "I1", "H1", "U1", "R1", "D1 ", "S2", "E2", "C2", "I2", "H2", "U2", "R2", "D2 "}, 16, 8);
// Uncomment this line to print the transitions.
// sim.get_transitions().print_table(
// {"S->E 1", "E->C 1", "C->I 1", "C->R 1", "I->H 1", "I->R 1", "H->U 1", "H->R 1", "U->D 1", "U->R 1"}, 16, 8);
// sim.get_transitions().print_table({"S->E 1", "E->C 1", "C->I 1", "C->R 1", "I->H 1", "I->R 1", "H->U 1",
// "H->R 1", "U->D 1", "U->R 1", "S->E 2", "E->C 2", "C->I 2", "C->R 2",
// "I->H 2", "I->R 2", "H->U 2", "H->R 2", "U->D 2", "U->R 2"},
// 16, 8);
}
Loading