diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 39a12539e1..4c63887dfb 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -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}) diff --git a/cpp/examples/ide_initialization.cpp b/cpp/examples/ide_initialization.cpp index 29f7280e8b..9b5d14971e 100644 --- a/cpp/examples/ide_initialization.cpp +++ b/cpp/examples/ide_initialization.cpp @@ -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 total_population = + mio::CustomIndexArray total_population_init = mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 80 * 1e6); - mio::CustomIndexArray deaths = mio::CustomIndexArray( + mio::CustomIndexArray deaths_init = mio::CustomIndexArray( 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((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((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; - mio::TimeSeries init(num_agegroups * (size_t)mio::isecir::InfectionTransition::Count); - init.add_time_point>( + mio::TimeSeries transitions_init((size_t)mio::isecir::InfectionTransition::Count * num_agegroups); + transitions_init.add_time_point>( -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(). @@ -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"}, diff --git a/cpp/examples/ide_secir.cpp b/cpp/examples/ide_secir.cpp index 461e127153..8c75efe6d2 100644 --- a/cpp/examples/ide_secir.cpp +++ b/cpp/examples/ide_secir.cpp @@ -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 * @@ -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" @@ -32,96 +31,138 @@ int main() { + // This is a simple example to demonstrate how use the IDE-SECIR model. + using Vec = Eigen::VectorX; - size_t num_agegroups = 1; + // Define simulation parameters. + ScalarType t0 = 0.; + ScalarType tmax = 5.; + ScalarType dt = 0.01; // The step size will stay constant throughout the simulation. - ScalarType tmax = 10; - mio::CustomIndexArray N = - mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 10000.); - mio::CustomIndexArray deaths = - mio::CustomIndexArray(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 total_population_init = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 1000.); + mio::CustomIndexArray deaths_init = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 6.); - // Create TimeSeries with num_transitions * num_agegroups elements where transitions needed for simulation will be - // stored. - mio::TimeSeries 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 transitions_init(num_transitions * num_agegroups); - // Add time points for initialization of transitions. + // Define vector of transitions that will be added as values 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 smoothcos(2.0); - mio::StateAgeFunctionWrapper delaydistribution(smoothcos); - std::vector> 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::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 smoothcos1(3.0); + mio::StateAgeFunctionWrapper delaydistribution1(smoothcos1); + std::vector> 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::AgeGroup(0)] = vec_delaydistrib1; + // TransitionProbabilities std::vector 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::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()[group] = vec_prob; + } + // Contact patterns mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, num_agegroups); contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, 10.)); - model.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + model.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + // Furhter epidemiological parameters mio::ExponentialSurvivalFunction exponential(0.5); mio::StateAgeFunctionWrapper prob(exponential); - - model.parameters.get()[mio::AgeGroup(0)] = prob; - model.parameters.get()[mio::AgeGroup(0)] = prob; - model.parameters.get()[mio::AgeGroup(0)] = prob; - + for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) { + model.parameters.get()[group] = prob; + model.parameters.get()[group] = prob; + model.parameters.get()[group] = prob; + } model.parameters.set(0.1); - // Start the simulation on the 40th day of a year (i.e. in February). - model.parameters.set(40); + model.parameters.set( + 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); } diff --git a/cpp/examples/ide_secir_ageres.cpp b/cpp/examples/ide_secir_ageres.cpp deleted file mode 100644 index 66ce473850..0000000000 --- a/cpp/examples/ide_secir_ageres.cpp +++ /dev/null @@ -1,149 +0,0 @@ -/* -* Copyright (C) 2020-2025 MEmilio -* -* Authors: Hannah Tritzschak, Anna Wendler -* -* Contact: Martin J. Kuehn -* -* 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. -*/ -#include "ide_secir/model.h" -#include "ide_secir/infection_state.h" -#include "ide_secir/simulation.h" -#include "memilio/config.h" -#include "memilio/epidemiology/age_group.h" -#include "memilio/math/eigen.h" -#include "memilio/utils/custom_index_array.h" -#include "memilio/utils/time_series.h" -#include "memilio/epidemiology/uncertain_matrix.h" -#include "memilio/epidemiology/state_age_function.h" -#include "memilio/data/analyze_result.h" - -int main() -{ - using Vec = Eigen::VectorX; - - size_t num_agegroups = 2; - - ScalarType tmax = 5; - mio::CustomIndexArray N = - mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 10000.); - mio::CustomIndexArray deaths = - mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 6.); - ScalarType dt = 1.; - - int num_transitions = (int)mio::isecir::InfectionTransition::Count; - - // Create TimeSeries with num_transitions * num_agegroups elements where transitions needed for simulation will be - // stored. - mio::TimeSeries init(num_transitions * num_agegroups); - - // Add time points for initialization of transitions. - Vec vec_init(num_transitions * num_agegroups); - // Values for the Infectiontransitions are the same for all AgeGroups. - for (size_t group = 0; group < num_agegroups; ++group) { - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::SusceptibleToExposed] = 25.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 15.0; - vec_init[group * num_transitions + - (int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered] = 4.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere] = - 1.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSymptomsToRecovered] = 4.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical] = - 1.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedSevereToRecovered] = 1.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedCriticalToDead] = 1.0; - vec_init[group * num_transitions + (int)mio::isecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0; - } - - // 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); - } - - // Initialize model. - mio::isecir::Model model(std::move(init), N, deaths, num_agegroups); - - // Uncomment these 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; - // model.populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Count + - // (Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000; - // Initialization method with Recovered. - // model.populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0; - // model.populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Count + - // (Eigen::Index)mio::isecir::InfectionState::Recovered] = 0; - - // Set working parameters. - // First AgeGroup for Transition Distributions. - mio::SmootherCosine smoothcos1(2.0); - mio::StateAgeFunctionWrapper delaydistribution1(smoothcos1); - std::vector> vec_delaydistrib1(num_transitions, delaydistribution1); - // TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value. - vec_delaydistrib1[(int)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.); - - model.parameters.get()[mio::AgeGroup(0)] = vec_delaydistrib1; - - //Second AgeGroup for Transition Distributions. - mio::SmootherCosine smoothcos2(3.0); - mio::StateAgeFunctionWrapper delaydistribution2(smoothcos2); - std::vector> vec_delaydistrib2(num_transitions, delaydistribution2); - // TransitionDistribution is not used for SusceptibleToExposed. Therefore, the parameter can be set to any value. - vec_delaydistrib2[(int)mio::isecir::InfectionTransition::SusceptibleToExposed].set_distribution_parameter(-1.); - - model.parameters.get()[mio::AgeGroup(1)] = vec_delaydistrib2; - - std::vector 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; - for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) { - model.parameters.get()[group] = vec_prob; - } - - mio::ContactMatrixGroup contact_matrix = - mio::ContactMatrixGroup(1, static_cast(num_agegroups)); - contact_matrix[0] = - mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, 10.)); - model.parameters.get() = mio::UncertainContactMatrix(contact_matrix); - - mio::ExponentialSurvivalFunction exponential(0.5); - mio::StateAgeFunctionWrapper prob(exponential); - for (mio::AgeGroup group = mio::AgeGroup(0); group < mio::AgeGroup(num_agegroups); ++group) { - model.parameters.get()[group] = prob; - model.parameters.get()[group] = prob; - model.parameters.get()[group] = prob; - } - model.parameters.set(0.1); - // Start the simulation on the 40th day of a year (i.e. in February). - model.parameters.set(40); - - model.check_constraints(dt); - - // Carry out simulation. - mio::isecir::Simulation sim(model, dt); - sim.advance(tmax); - - auto interpolated_results = mio::interpolate_simulation_result(sim.get_result(), dt / 2.); - - 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", "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); -} diff --git a/cpp/examples/ide_secir_graph.cpp b/cpp/examples/ide_secir_graph.cpp index 7692a61d50..e3a5c822f5 100644 --- a/cpp/examples/ide_secir_graph.cpp +++ b/cpp/examples/ide_secir_graph.cpp @@ -37,9 +37,11 @@ int main() // Set time step size of IDE solver. Note that the time step size will be constant throughout the simulation. const ScalarType dt_ide_solver = 1.; + // Define number of age groups. size_t num_agegroups = 1; - mio::CustomIndexArray N_init = + // Define initial values for the total population and number of deaths. + mio::CustomIndexArray total_population_init = mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 10000.); mio::CustomIndexArray deaths_init = mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 13.10462213); @@ -77,7 +79,7 @@ int main() } // Initialize IDE model that will be used in in each node in the graph below. - mio::isecir::Model model(std::move(transitions_init), N_init, deaths_init, num_agegroups); + mio::isecir::Model model(std::move(transitions_init), total_population_init, deaths_init, num_agegroups); model.check_constraints(dt_ide_solver); // Set up graph with two nodes and no edges. To each node, we pass an id, the above constructed IDE model as well diff --git a/cpp/memilio/epidemiology/state_age_function.h b/cpp/memilio/epidemiology/state_age_function.h index 6387a60110..6b18545883 100644 --- a/cpp/memilio/epidemiology/state_age_function.h +++ b/cpp/memilio/epidemiology/state_age_function.h @@ -276,7 +276,8 @@ struct StateAgeFunction { * * This is a basic version to determine the mean value of a survival function * through numerical integration of the integral that describes the expected value. - * This basic implementation is only valid if the StateAgeFunction is of type a). Otherwise it should be overridden. + * This basic implementation is only valid if the StateAgeFunction is of type a) since we assume that the considered + * function converges to zero. Otherwise it should be overridden. * * For some specific derivations of StateAgeFunction%s there are more efficient ways to determine the * mean value which is why this member function is virtual and can be overridden (see, e.g., ExponentialSurvivalFunction). @@ -290,9 +291,13 @@ struct StateAgeFunction { { using std::ceil; if (!floating_point_equal(m_mean_tol, tol, 1e-14) || floating_point_equal(m_mean, -1., 1e-14)) { - // Integration using Trapezoidal rule. + // Integration using trapezoidal rule. + // Compute value for i=0. FP mean = 0.5 * dt * eval(FP(0 * dt)); FP supp_max_idx = ceil(get_support_max(dt, tol) / dt); + // We start with i=1 since the value for i=0 was already considered above when defining the variable mean. + // Note that we do not consider indices i>=supp_max_idx since it holds eval(i*dt)=supp_max_idx + // by definition of the support_max. for (int i = 1; i < supp_max_idx; i++) { mean += dt * eval(FP(i * dt)); } diff --git a/cpp/models/ide_secir/model.cpp b/cpp/models/ide_secir/model.cpp index 7c45c640ba..f7096e52e8 100644 --- a/cpp/models/ide_secir/model.cpp +++ b/cpp/models/ide_secir/model.cpp @@ -181,6 +181,14 @@ ScalarType Model::get_global_support_max(ScalarType dt) const return global_support_max; } +/* Disclaimer: The following discrete formulas are derived and explained in +Wendler A, Plötzke L, Tritzschak H, Kühn MJ: A nonstandard numerical scheme for a novel SECIR integro differential +equation-based model with nonexponentially distributed stay times. DOI:10.1016/j.amc.2025.129636. +Note that throughout this implementation, we assume that the initial values for the transitions contain information +about the number of individuals transitioning within a time interval (in contrast to at a time point), cf. Section 3.4 in +the paper. The formulas from the paper are adapted accordingly by scaling with the time step size dt. + */ + // ---- Functionality to calculate the sizes of the compartments for time t0. ---- void Model::compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_InfectionState, AgeGroup group, Eigen::Index idx_IncomingFlow, int idx_TransitionDistribution1,