From 421d728c053487c5ae98df8bcaeaa721f8dd26e4 Mon Sep 17 00:00:00 2001 From: Carlotta Gerstein <100771374+charlie0614@users.noreply.github.com> Date: Mon, 20 Apr 2026 12:55:43 +0200 Subject: [PATCH 01/12] switch to ScalarType in example docs --- cpp/examples/ode_secirvvs.cpp | 81 ++++++++++++++-------------- docs/source/cpp/models/omseirs4.rst | 6 +-- docs/source/cpp/models/osecir.rst | 18 +++---- docs/source/cpp/models/osecirts.rst | 16 +++--- docs/source/cpp/models/osecirvvs.rst | 26 ++++----- docs/source/cpp/models/oseir.rst | 16 +++--- docs/source/cpp/models/oseirdb.rst | 16 +++--- docs/source/cpp/models/oseirv.rst | 18 +++---- docs/source/cpp/models/osir.rst | 14 ++--- 9 files changed, 106 insertions(+), 105 deletions(-) diff --git a/cpp/examples/ode_secirvvs.cpp b/cpp/examples/ode_secirvvs.cpp index cfdd39b19d..8d9a1aba44 100644 --- a/cpp/examples/ode_secirvvs.cpp +++ b/cpp/examples/ode_secirvvs.cpp @@ -25,13 +25,13 @@ int main() { mio::set_log_level(mio::LogLevel::warn); - double t0 = 0; - double tmax = 30; - double dt = 0.1; + ScalarType t0 = 0; + ScalarType tmax = 30; + ScalarType dt = 0.1; mio::log_info("Simulating SECIRVVS; t={} ... {} with dt = {}.", t0, tmax, dt); - mio::osecirvvs::Model model(1); + mio::osecirvvs::Model model(1); for (mio::AgeGroup i = 0; i < model.parameters.get_num_groups(); i++) { model.populations[{i, mio::osecirvvs::InfectionState::ExposedNaive}] = 10; @@ -64,59 +64,60 @@ int main() {i, mio::osecirvvs::InfectionState::SusceptibleNaive}, 1000); } - model.parameters.get>() = 100; - model.parameters.get>() = 0.0143; - const size_t daily_vaccinations = 10; - model.parameters.get>().resize( + model.parameters.get>() = 100; + model.parameters.get>() = 0.0143; + const size_t daily_vaccinations = 10; + model.parameters.get>().resize( + mio::SimulationDay((size_t)tmax + 1)); + model.parameters.get>().resize( mio::SimulationDay((size_t)tmax + 1)); - model.parameters.get>().resize(mio::SimulationDay((size_t)tmax + 1)); for (size_t i = 0; i < tmax + 1; ++i) { - auto num_vaccinations = static_cast(i * daily_vaccinations); + auto num_vaccinations = static_cast(i * daily_vaccinations); model.parameters - .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = + .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = num_vaccinations; model.parameters - .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = + .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = num_vaccinations; } - model.parameters.get>() = 7; + model.parameters.get>() = 7; - auto& contacts = model.parameters.get>(); + auto& contacts = model.parameters.get>(); auto& contact_matrix = contacts.get_cont_freq_mat(); contact_matrix[0].get_baseline().setConstant(0.5); contact_matrix[0].get_baseline().diagonal().setConstant(5.0); - contact_matrix[0].add_damping(0.3, mio::SimulationTime(5.0)); + contact_matrix[0].add_damping(0.3, mio::SimulationTime(5.0)); //times - model.parameters.get>()[mio::AgeGroup(0)] = 3.33; - model.parameters.get>()[mio::AgeGroup(0)] = 1.87; - model.parameters.get>()[mio::AgeGroup(0)] = 7; - model.parameters.get>()[mio::AgeGroup(0)] = 6; - model.parameters.get>()[mio::AgeGroup(0)] = 7; + model.parameters.get>()[mio::AgeGroup(0)] = 3.33; + model.parameters.get>()[mio::AgeGroup(0)] = 1.87; + model.parameters.get>()[mio::AgeGroup(0)] = 7; + model.parameters.get>()[mio::AgeGroup(0)] = 6; + model.parameters.get>()[mio::AgeGroup(0)] = 7; //probabilities - model.parameters.get>()[mio::AgeGroup(0)] = 0.15; - model.parameters.get>()[mio::AgeGroup(0)] = 0.5; + model.parameters.get>()[mio::AgeGroup(0)] = 0.15; + model.parameters.get>()[mio::AgeGroup(0)] = 0.5; // The precise value between Risk* (situation under control) and MaxRisk* (situation not under control) // depends on incidence and test and trace capacity - model.parameters.get>()[mio::AgeGroup(0)] = 0.0; - model.parameters.get>()[mio::AgeGroup(0)] = 0.4; - model.parameters.get>()[mio::AgeGroup(0)] = 0.2; - model.parameters.get>()[mio::AgeGroup(0)] = 0.1; - model.parameters.get>()[mio::AgeGroup(0)] = 0.1; - model.parameters.get>()[mio::AgeGroup(0)] = 0.1; + model.parameters.get>()[mio::AgeGroup(0)] = 0.0; + model.parameters.get>()[mio::AgeGroup(0)] = 0.4; + model.parameters.get>()[mio::AgeGroup(0)] = 0.2; + model.parameters.get>()[mio::AgeGroup(0)] = 0.1; + model.parameters.get>()[mio::AgeGroup(0)] = 0.1; + model.parameters.get>()[mio::AgeGroup(0)] = 0.1; - model.parameters.get>()[mio::AgeGroup(0)] = 0.8; - model.parameters.get>()[mio::AgeGroup(0)] = 0.331; - model.parameters.get>()[mio::AgeGroup(0)] = 0.65; - model.parameters.get>()[mio::AgeGroup(0)] = 0.243; - model.parameters.get>()[mio::AgeGroup(0)] = - 0.1; - model.parameters.get>()[mio::AgeGroup(0)] = - 0.091; - model.parameters.get>()[mio::AgeGroup(0)] = 0.9; + model.parameters.get>()[mio::AgeGroup(0)] = 0.8; + model.parameters.get>()[mio::AgeGroup(0)] = 0.331; + model.parameters.get>()[mio::AgeGroup(0)] = 0.65; + model.parameters.get>()[mio::AgeGroup(0)] = 0.243; + model.parameters + .get>()[mio::AgeGroup(0)] = 0.1; + model.parameters + .get>()[mio::AgeGroup(0)] = 0.091; + model.parameters.get>()[mio::AgeGroup(0)] = 0.9; - model.parameters.get>() = 0.2; + model.parameters.get>() = 0.2; // The function apply_constraints() ensures that all parameters are within their defined bounds. // Note that negative values are set to zero instead of stopping the simulation. model.apply_constraints(); @@ -127,10 +128,10 @@ int main() // integrator->set_dt_max(1.0); // integrator->set_rel_tolerance(1e-4); // integrator->set_abs_tolerance(1e-1); - // mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model, std::move(integrator)); + // mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model, std::move(integrator)); // use default Cash-Karp adaptive integrator - mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model); + mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model); bool print_to_terminal = true; diff --git a/docs/source/cpp/models/omseirs4.rst b/docs/source/cpp/models/omseirs4.rst index 17b7cb216a..5da37b9ec7 100644 --- a/docs/source/cpp/models/omseirs4.rst +++ b/docs/source/cpp/models/omseirs4.rst @@ -123,9 +123,9 @@ Run a standard simulation via: .. code-block:: cpp - double t0 = 0.0; // days - double tmax = 3650; // 10 years - double dt = 1.0; // daily step + ScalarType t0 = 0.0; // days + ScalarType tmax = 3650; // 10 years + ScalarType dt = 1.0; // daily step auto timeseries = mio::simulate(t0, tmax, dt, model); Output diff --git a/docs/source/cpp/models/osecir.rst b/docs/source/cpp/models/osecir.rst index 50f651c698..5e1cd66293 100644 --- a/docs/source/cpp/models/osecir.rst +++ b/docs/source/cpp/models/osecir.rst @@ -181,7 +181,7 @@ Basic dampings can be added to the contact matrix as follows: .. code-block:: cpp // Create a contact matrix with constant contact rates between all groups - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); // Add a damping that reduces contacts by 70% starting at day 30 @@ -270,7 +270,7 @@ It will be active for at least 14 days and checked every 3 days. If the last che .. code-block:: cpp // Configure dynamic NPIs with thresholds - auto& dynamic_npis = params.get>(); + auto& dynamic_npis = params.get>(); dynamic_npis.set_interval(mio::SimulationTime(3.0)); // Check every 3 days dynamic_npis.set_duration(mio::SimulationTime(14.0)); // Apply for 14 days dynamic_npis.set_base_value(100'000); // Per 100,000 population @@ -289,12 +289,12 @@ Standard simulation: .. code-block:: cpp - double t0 = 0; // Start time - double tmax = 50; // End time - double dt = 0.1; // Time step + ScalarType t0 = 0; // Start time + ScalarType tmax = 50; // End time + ScalarType dt = 0.1; // Time step // Run a standard simulation - mio::TimeSeries secir = mio::osecir::simulate(t0, tmax, dt, model); + mio::TimeSeries secir = mio::osecir::simulate(t0, tmax, dt, model); Flow simulation for tracking transitions between compartments: @@ -314,7 +314,7 @@ For both simulation types, you can also specify a custom integrator: integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - mio::TimeSeries secir = mio::osecir::simulate(t0, tmax, dt, model, std::move(integrator)); + mio::TimeSeries secir = mio::osecir::simulate(t0, tmax, dt, model, std::move(integrator)); Output @@ -329,11 +329,11 @@ The output of the simulation is a `TimeSeries` object containing the sizes of ea // Access data at a specific time point Eigen::VectorXd value_at_time_i = secir.get_value(i); - double time_i = secir.get_time(i); + ScalarType time_i = secir.get_time(i); // Access the last time point Eigen::VectorXd last_value = secir.get_last_value(); - double last_time = secir.get_last_time(); + ScalarType last_time = secir.get_last_time(); For flow simulations, the result consists of two `mio::TimeSeries` objects, one for compartment sizes and one for flows: diff --git a/docs/source/cpp/models/osecirts.rst b/docs/source/cpp/models/osecirts.rst index d8fec1bfec..db54cefc7c 100644 --- a/docs/source/cpp/models/osecirts.rst +++ b/docs/source/cpp/models/osecirts.rst @@ -302,7 +302,7 @@ Basic dampings can be added to the contact matrix as follows: .. code-block:: cpp // Create a contact matrix with baseline contact rates - auto& contacts = model.parameters.get>(); + auto& contacts = model.parameters.get>(); auto& contact_matrix = contacts.get_cont_freq_mat(); contact_matrix[0].get_baseline().setConstant(0.5); contact_matrix[0].get_baseline().diagonal().setConstant(5.0); @@ -315,7 +315,7 @@ The model also supports dynamic NPIs based on epidemic thresholds: .. code-block:: cpp // Set threshold-based triggers for NPIs - auto& dynamic_npis = model.parameters.get>(); + auto& dynamic_npis = model.parameters.get>(); dynamic_npis.set_interval(mio::SimulationTime(3.0)); // Check every 3 days dynamic_npis.set_duration(mio::SimulationTime(14.0)); // Apply for 14 days dynamic_npis.set_base_value(100'000); // Per 100,000 population @@ -335,12 +335,12 @@ Standard simulation: .. code-block:: cpp - double t0 = 0; // Start time - double tmax = 50; // End time - double dt = 0.1; // Time step + ScalarType t0 = 0; // Start time + ScalarType tmax = 50; // End time + ScalarType dt = 0.1; // Time step // Run a standard simulation - mio::TimeSeries result = mio::osecirts::simulate(t0, tmax, dt, model); + mio::TimeSeries result = mio::osecirts::simulate(t0, tmax, dt, model); During simulation, the model handles several special processes: @@ -358,7 +358,7 @@ For both simulation types, you can also specify a custom integrator: integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - mio::TimeSeries result = mio::osecirts::simulate(t0, tmax, dt, model, std::move(integrator)); + mio::TimeSeries result = mio::osecirts::simulate(t0, tmax, dt, model, std::move(integrator)); Output ------ @@ -372,7 +372,7 @@ The output of the simulation is a `mio::TimeSeries` object containing the sizes // Access data at a specific time point Eigen::VectorXd value_at_time_i = result.get_value(i); - double time_i = result.get_time(i); + ScalarType time_i = result.get_time(i); // Access the last time point Eigen::VectorXd last_value = result.get_last_value(); diff --git a/docs/source/cpp/models/osecirvvs.rst b/docs/source/cpp/models/osecirvvs.rst index 6733210fe8..baa269f03b 100644 --- a/docs/source/cpp/models/osecirvvs.rst +++ b/docs/source/cpp/models/osecirvvs.rst @@ -265,19 +265,19 @@ After setting the initial populations, you also need to set the vaccination para // Prepare and resize vaccinations parameter for the entire simulation period const size_t daily_vaccinations = 10; - model.parameters.get>().resize( + model.parameters.get>().resize( mio::SimulationDay((size_t)tmax + 1)); - model.parameters.get>().resize( + model.parameters.get>().resize( mio::SimulationDay((size_t)tmax + 1)); // Set increasing number of vaccination over time for (size_t i = 0; i < tmax + 1; ++i) { - auto num_vaccinations = static_cast(i * daily_vaccinations); + auto num_vaccinations = static_cast(i * daily_vaccinations); model.parameters - .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = + .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = num_vaccinations; model.parameters - .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = + .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = num_vaccinations; } @@ -293,7 +293,7 @@ Basic dampings can be added to the contact matrix as follows: .. code-block:: cpp // Create a contact matrix with baseline contact rates - auto& contacts = model.parameters.get>(); + auto& contacts = model.parameters.get>(); auto& contact_matrix = contacts.get_cont_freq_mat(); contact_matrix[0].get_baseline().setConstant(0.5); contact_matrix[0].get_baseline().diagonal().setConstant(5.0); @@ -331,7 +331,7 @@ The model also supports dynamic NPIs based on epidemic thresholds: .. code-block:: cpp // Configure dynamic NPIs - auto& dynamic_npis = params.get>(); + auto& dynamic_npis = params.get>(); dynamic_npis.set_interval(mio::SimulationTime(3.0)); // Check NPI every 3 days dynamic_npis.set_duration(mio::SimulationTime(14.0)); // Apply NPI for 14 days dynamic_npis.set_base_value(100'000); // Base value to trigger NPI is population of 100,000 @@ -349,12 +349,12 @@ Basic simulation: .. code-block:: cpp - double t0 = 0; // Start time - double tmax = 30; // End time - double dt = 0.1; // Time step + ScalarType t0 = 0; // Start time + ScalarType tmax = 30; // End time + ScalarType dt = 0.1; // Time step // Run a standard simulation - mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model); + mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model); During simulation, the model handles several special processes: @@ -372,7 +372,7 @@ For both simulation types, you can also specify a custom integrator: integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model, std::move(integrator)); + mio::TimeSeries result = mio::osecirvvs::simulate(t0, tmax, dt, model, std::move(integrator)); Output ------ @@ -386,7 +386,7 @@ The output of the simulation is a `mio::TimeSeries` object containing the sizes // Access data at a specific time point Eigen::VectorXd value_at_time_i = result.get_value(i); - double time_i = result.get_time(i); + ScalarType time_i = result.get_time(i); // Access the last time point Eigen::VectorXd last_value = result.get_last_value(); diff --git a/docs/source/cpp/models/oseir.rst b/docs/source/cpp/models/oseir.rst index a22075ea7d..5b270a2edb 100644 --- a/docs/source/cpp/models/oseir.rst +++ b/docs/source/cpp/models/oseir.rst @@ -131,7 +131,7 @@ Basic dampings can be added to the ContactPatterns as follows: .. code-block:: cpp // Create a contact matrix with constant contact rates between all groups - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); // Add a damping that reduces contacts by 70% starting at day 30 @@ -160,12 +160,12 @@ Standard simulation: .. code-block:: cpp - double t0 = 0; // Start time - double tmax = 50; // End time - double dt = 0.1; // Initial step size + ScalarType t0 = 0; // Start time + ScalarType tmax = 50; // End time + ScalarType dt = 0.1; // Initial step size // Run a standard simulation - mio::TimeSeries result_sim = mio::oseir::simulate(t0, tmax, dt, model); + mio::TimeSeries result_sim = mio::oseir::simulate(t0, tmax, dt, model); Flow simulation for tracking transitions between compartments: @@ -185,7 +185,7 @@ For both simulation types, you can also specify a custom integrator: integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - mio::TimeSeries result_sim = mio::oseir::simulate(t0, tmax, dt, model, std::move(integrator)); + mio::TimeSeries result_sim = mio::oseir::simulate(t0, tmax, dt, model, std::move(integrator)); Output ------ @@ -200,11 +200,11 @@ standard simulation, you can access the results as follows: // Access data at specific time point Eigen::VectorXd value_at_time_point_i = result_sim.get_value(i); - double time_i = result_sim.get_time(i); + ScalarType time_i = result_sim.get_time(i); // Access the last time point Eigen::VectorXd last_value = result_sim.get_last_value(); - double last_time = result_sim.get_last_time(); + ScalarType last_time = result_sim.get_last_time(); For flow simulations, the result consists of two `mio::TimeSeries` objects, one for compartment sizes and one for flows: diff --git a/docs/source/cpp/models/oseirdb.rst b/docs/source/cpp/models/oseirdb.rst index 809dbda49c..527708d2d1 100644 --- a/docs/source/cpp/models/oseirdb.rst +++ b/docs/source/cpp/models/oseirdb.rst @@ -140,7 +140,7 @@ Basic dampings can be added to the ContactPatterns as follows: .. code-block:: cpp // Create a contact matrix with constant contact rates between all groups - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); // Add a damping that reduces contacts by 70% starting at day 30 @@ -169,12 +169,12 @@ Standard simulation: .. code-block:: cpp - double t0 = 0; // Start time - double tmax = 50; // End time - double dt = 0.1; // Initial step size + ScalarType t0 = 0; // Start time + ScalarType tmax = 50; // End time + ScalarType dt = 0.1; // Initial step size // Run a standard simulation - mio::TimeSeries result_sim = mio::oseirdb::simulate(t0, tmax, dt, model); + mio::TimeSeries result_sim = mio::oseirdb::simulate(t0, tmax, dt, model); Flow simulation for tracking transitions between compartments: @@ -194,7 +194,7 @@ For both simulation types, you can also specify a custom integrator: integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - mio::TimeSeries result_sim = mio::oseirdb::simulate(t0, tmax, dt, model, std::move(integrator)); + mio::TimeSeries result_sim = mio::oseirdb::simulate(t0, tmax, dt, model, std::move(integrator)); Output ------ @@ -209,11 +209,11 @@ standard simulation, you can access the results as follows: // Access data at specific time point Eigen::VectorXd value_at_time_point_i = result_sim.get_value(i); - double time_i = result_sim.get_time(i); + ScalarType time_i = result_sim.get_time(i); // Access the last time point Eigen::VectorXd last_value = result_sim.get_last_value(); - double last_time = result_sim.get_last_time(); + ScalarType last_time = result_sim.get_last_time(); For flow simulations, the result consists of two `mio::TimeSeries` objects, one for compartment sizes and one for flows: diff --git a/docs/source/cpp/models/oseirv.rst b/docs/source/cpp/models/oseirv.rst index b0bff91a3b..38ff17ecbc 100644 --- a/docs/source/cpp/models/oseirv.rst +++ b/docs/source/cpp/models/oseirv.rst @@ -52,7 +52,7 @@ broadly). The number of age groups is specified in the constructor: .. code-block:: cpp - mio::oseirv::Model model(num_agegroups); + mio::oseirv::Model model(num_agegroups); For stratifications with two or more dimensions, see :doc:`Model Creation <../ode_creation>`. @@ -127,12 +127,12 @@ Initial conditions are handled via the **Populations** class. Example for a sing .. code-block:: cpp - mio::oseirv::Model model(1); + mio::oseirv::Model model(1); // Set total population in age group 0 model.populations.set_total(total0); // Initialize vaccinated susceptibles (simple example) - double vc0 = 0.4; // vaccination coverage + ScalarType vc0 = 0.4; // vaccination coverage model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::SusceptibleVaccinated}] = vc0 * total0; model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = initial_infected; model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Exposed}] = initial_exposed; @@ -166,9 +166,9 @@ same pattern. Example with a Runge–Kutta integrator: .. code-block:: cpp - double t0 = 0.0; // start (weeks) - double tmax = 20.0; // end - double dt = 0.1; // initial step size + ScalarType t0 = 0.0; // start (weeks) + ScalarType tmax = 20.0; // end + ScalarType dt = 0.1; // initial step size auto integrator = std::make_unique(); integrator->set_dt_min(0.01); @@ -195,7 +195,7 @@ The result of a standard simulation is a ``mio::TimeSeries``: auto n_points = static_cast(sim.get_num_time_points()); Eigen::VectorXd val_i = sim.get_value(i); - double time_i = sim.get_time(i); + ScalarType time_i = sim.get_time(i); auto last_val = sim.get_last_value(); Printing and CSV export: @@ -217,8 +217,8 @@ Time-dependent changes of contact patterns (holidays, interventions) can be mode .. code-block:: cpp - mio::ContactMatrixGroup& cm_h = model.parameters.get>(); - mio::ContactMatrixGroup& cm_s = model.parameters.get>(); + mio::ContactMatrixGroup& cm_h = model.parameters.get>(); + mio::ContactMatrixGroup& cm_s = model.parameters.get>(); cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, base_contacts)); cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, base_contacts_sick)); diff --git a/docs/source/cpp/models/osir.rst b/docs/source/cpp/models/osir.rst index c7d0ce4614..e608c153ed 100644 --- a/docs/source/cpp/models/osir.rst +++ b/docs/source/cpp/models/osir.rst @@ -137,12 +137,12 @@ time: .. code-block:: cpp - double t0 = 0; // Start time - double tmax = 50; // End time - double dt = 0.1; // Time step + ScalarType t0 = 0; // Start time + ScalarType tmax = 50; // End time + ScalarType dt = 0.1; // Time step // Run a standard simulation - mio::TimeSeries result = mio::simulate(t0, tmax, dt, model); + mio::TimeSeries result = mio::simulate(t0, tmax, dt, model); You can also specify a custom integrator: @@ -154,7 +154,7 @@ You can also specify a custom integrator: integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - mio::TimeSeries result = mio::simulate(t0, tmax, dt, model, std::move(integrator)); + mio::TimeSeries result = mio::simulate(t0, tmax, dt, model, std::move(integrator)); Output @@ -170,11 +170,11 @@ a basic simulation, you can access the results as follows: // Access data at specific time point Eigen::VectorXd value_at_time_i = result.get_value(i); - double time_i = result.get_time(i); + ScalarType time_i = result.get_time(i); // Access the last time point Eigen::VectorXd last_value = result.get_last_value(); - double last_time = result.get_last_time(); + ScalarType last_time = result.get_last_time(); You can print the simulation results as a formatted table: From c90a58d4a865585b73f4f3799b8e6c2a51dc55fe Mon Sep 17 00:00:00 2001 From: Anna Wendler <106674756+annawendler@users.noreply.github.com> Date: Wed, 22 Apr 2026 14:45:38 +0200 Subject: [PATCH 02/12] update code examples in docu for lct and ide models --- docs/source/cpp/models/glsecir.rst | 139 ++++++++++++++-------------- docs/source/cpp/models/isecir.rst | 31 ++++--- docs/source/cpp/models/iseir.rst | 28 +++--- docs/source/cpp/models/lsecir.rst | 109 +++++++++++----------- docs/source/cpp/models/lsecir2d.rst | 103 +++++++++++---------- 5 files changed, 210 insertions(+), 200 deletions(-) diff --git a/docs/source/cpp/models/glsecir.rst b/docs/source/cpp/models/glsecir.rst index f7a1072fa8..fb53ec91d9 100644 --- a/docs/source/cpp/models/glsecir.rst +++ b/docs/source/cpp/models/glsecir.rst @@ -125,14 +125,14 @@ Note that in the GLCT model, we define two strains for the compartments `Infecte .. code-block:: cpp constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 6, NumInfectedSymptoms = 2, NumInfectedSevere = 2, - NumInfectedCritical = 10; - using Model = mio::glsecir::Model; - using LctState = Model::LctState; + NumInfectedCritical = 10; + using Model = mio::glsecir::Model; + using LctState = Model::LctState; using InfectionState = LctState::InfectionState; Model model; - + We continue by defining some epidemiological parameters needed throughout the model definition and initialization. .. code-block:: cpp @@ -164,13 +164,14 @@ We continue by defining some epidemiological parameters needed throughout the mo {750}, // Susceptible {30, 20}, // Exposed {20 * (1 - recoveredPerInfectedNoSymptoms), 10 * (1 - recoveredPerInfectedNoSymptoms), // InfectedNoSymptoms - 10 * (1 - recoveredPerInfectedNoSymptoms), 20 * recoveredPerInfectedNoSymptoms, - 10 * recoveredPerInfectedNoSymptoms, 10 * recoveredPerInfectedNoSymptoms}, + 10 * (1 - recoveredPerInfectedNoSymptoms), 20 * recoveredPerInfectedNoSymptoms, + 10 * recoveredPerInfectedNoSymptoms, 10 * recoveredPerInfectedNoSymptoms}, {50 * severePerInfectedSymptoms, 50 * (1 - severePerInfectedSymptoms)}, // InfectedSymptoms {50 * criticalPerSevere, 50 * (1 - criticalPerSevere)}, // InfectedSevere - {10 * deathsPerCritical, 10 * deathsPerCritical, 5 * deathsPerCritical, 3 * deathsPerCritical, // InfectedCritical - 2 * deathsPerCritical, 10 * (1 - deathsPerCritical), 10 * (1 - deathsPerCritical), 5 * (1 - deathsPerCritical), - 3 * (1 - deathsPerCritical), 2 * (1 - deathsPerCritical)}, + {10 * deathsPerCritical, 10 * deathsPerCritical, 5 * deathsPerCritical, + 3 * deathsPerCritical, // InfectedCritical + 2 * deathsPerCritical, 10 * (1 - deathsPerCritical), 10 * (1 - deathsPerCritical), 5 * (1 - deathsPerCritical), + 3 * (1 - deathsPerCritical), 2 * (1 - deathsPerCritical)}, {20}, // Recovered {10}}; // Dead @@ -183,18 +184,18 @@ Below, we assert that ``initial_populations`` has the right shape. return 1; } if ((initial_populations[(size_t)InfectionState::Susceptible].size() != - LctState::get_num_subcompartments()) || + LctState::get_num_subcompartments()) || (initial_populations[(size_t)InfectionState::Exposed].size() != NumExposed) || (initial_populations[(size_t)InfectionState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) || (initial_populations[(size_t)InfectionState::InfectedSymptoms].size() != NumInfectedSymptoms) || (initial_populations[(size_t)InfectionState::InfectedSevere].size() != NumInfectedSevere) || (initial_populations[(size_t)InfectionState::InfectedCritical].size() != NumInfectedCritical) || (initial_populations[(size_t)InfectionState::Recovered].size() != - LctState::get_num_subcompartments()) || + LctState::get_num_subcompartments()) || (initial_populations[(size_t)InfectionState::Dead].size() != - LctState::get_num_subcompartments())) { + LctState::get_num_subcompartments())) { mio::log_error("The length of at least one vector in initial_populations does not match the related number of " - "subcompartments."); + "subcompartments."); return 1; } @@ -210,7 +211,6 @@ Finally, we transfer the initial values in ``initial_populations`` to the model. model.populations[mio::Index(i)] = flat_initial_populations[i]; } - Since we want to recreate the LCT model as defined in the corresponding example, we set the parameters determining the transition distributions such that we obtain Erlang distributions. We will explain how to do this for the different compartments in the following. In general, we need to define a vector ``StartingProbabilities(...)`` and a matrix ``TransitionMatrix(...z)To(...*)`` for all compartments with subcompartments, i.e. `Exposed`, `InfectedNoSymptoms`, `InfectedSymptoms`, `InfectedSevere` and `InfectedCritical`. @@ -226,18 +226,17 @@ The get_default of the ``StartingProbabilities(...)`` returns the first unit vec .. code-block:: cpp - model.parameters.get() = - mio::glsecir::StartingProbabilitiesExposed().get_default( + model.parameters.get>() = + mio::glsecir::StartingProbabilitiesExposed().get_default( LctState::get_num_subcompartments()); The get_default function returns the ``TransitionMatrix`` that is required to have an Erlang-distributed stay time with an average of timeExposed. .. code-block:: cpp - model.parameters.get() = - mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default( - LctState::get_num_subcompartments(), timeExposed); - + model.parameters.get>() = + mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default( + LctState::get_num_subcompartments(), timeExposed); We continue with the compartment `InfectedNoSymptoms`. For InfectedNoSymptoms, two strains have to be defined, one for the transition `InfectedNoSymptomsToInfectedSymptoms` and one for the transition `InfectedNoSymptomsToRecovered`. The strains have a length of ``NumInfectedNoSymptoms/2`` each as we choose the same number of subcompartments for both strains. Note that the transition probability is included in the vector ``StartingProbabilitiesInfectedNoSymptoms``. @@ -249,19 +248,19 @@ The strains have a length of ``NumInfectedNoSymptoms/2`` each as we choose the s StartingProbabilitiesInfectedNoSymptoms[0] = 1 - recoveredPerInfectedNoSymptoms; StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = recoveredPerInfectedNoSymptoms; - model.parameters.get() = + model.parameters.get>() = StartingProbabilitiesInfectedNoSymptoms; Equal transition matrices for the strains have to be defined. They follow the same Erlang distribution such that we get the same result as with the LCT model that can only consider one strain. .. code-block:: cpp - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms().get_default( + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedNoSymptoms); - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default( + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedNoSymptoms); @@ -275,12 +274,13 @@ We proceed analogously for the remaining compartments `InfectedSymptoms`, `Infec StartingProbabilitiesInfectedSymptoms[0] = severePerInfectedSymptoms; StartingProbabilitiesInfectedSymptoms[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = 1 - severePerInfectedSymptoms; - model.parameters.get() = StartingProbabilitiesInfectedSymptoms; - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere().get_default( + model.parameters.get>() = + StartingProbabilitiesInfectedSymptoms; + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSymptoms); - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSymptoms); // InfectedSevere. @@ -289,12 +289,13 @@ We proceed analogously for the remaining compartments `InfectedSymptoms`, `Infec StartingProbabilitiesInfectedSevere[0] = criticalPerSevere; StartingProbabilitiesInfectedSevere[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = 1 - criticalPerSevere; - model.parameters.get() = StartingProbabilitiesInfectedSevere; - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical().get_default( + model.parameters.get>() = + StartingProbabilitiesInfectedSevere; + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSevere); - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSevere); // InfectedCritical. @@ -303,15 +304,15 @@ We proceed analogously for the remaining compartments `InfectedSymptoms`, `Infec StartingProbabilitiesInfectedCritical[0] = deathsPerCritical; StartingProbabilitiesInfectedCritical[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = 1 - deathsPerCritical; - model.parameters.get() = StartingProbabilitiesInfectedCritical; - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedCriticalToDead().get_default( + model.parameters.get>() = + StartingProbabilitiesInfectedCritical; + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedCriticalToDead().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedCritical); - model.parameters.get() = - mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default( + model.parameters.get>() = + mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedCritical); - .. _Nonpharmaceutical Interventions: Nonpharmaceutical Interventions ------------------------------- @@ -325,24 +326,24 @@ Basic dampings can be added to the contact matrix as follows: // Create a contact matrix with constant contact rates between all groups. ScalarType cont_freq = 10.; - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); - + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); + // Add a uniform damping across all age groups. - contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); For age-resolved models, you can apply different dampings to different groups: .. code-block:: cpp ScalarType cont_freq = 10.; - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, cont_freq)); - - // Add a damping that reduces contacts within the same age group by 70% starting at day 30. - contact_matrix.add_damping(Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(), - mio::SimulationTime(30.)); - + contact_matrix[0] = + mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, cont_freq)); + // Add a damping that reduces contacts within the same age group by 70% starting at day 30. + Eigen::MatrixX damping_matrix = Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(); + contact_matrix.add_damping(damping_matrix, mio::SimulationTime(30.)); For more complex scenarios, such as real-world venue closures or lockdown modeling, you can implement detailed NPIs with location-specific dampings. The GLCT-SECIR model supports contact matrices for different locations (e.g., home, school, work, other) and can apply different dampings to each location. @@ -400,22 +401,23 @@ Simulating the model from :math:`t_0` to :math:`t_{\max}` with initial step size .. code-block:: cpp - const ScalarType t0 = 0; - const ScalarType tmax = 10; - const ScalarType dt_init = 10; - mio::TimeSeries result = mio::simulate(t0, tmax, dt_init, model); + const ScalarType t0 = 0; + const ScalarType tmax = 10; + const ScalarType dt_init = 10; + mio::TimeSeries result = mio::simulate(t0, tmax, dt_init, model); You can also specify a custom integrator: .. code-block:: cpp - auto integrator = std::make_unique(); + auto integrator = std::make_unique>(); integrator->set_dt_min(0.3); integrator->set_dt_max(1.0); integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - - mio::TimeSeries result = mio::simulate(t0, tmax, dt, model, std::move(integrator)); + + mio::TimeSeries result = + mio::simulate(t0, tmax, dt_init, model, std::move(integrator)); Output ------ @@ -432,14 +434,15 @@ You can access the data in the `mio::TimeSeries` object as follows: // Get the number of time points. auto num_points = static_cast(result.get_num_time_points()); - - // Access data at a specific time point. - Eigen::VectorX value_at_time_i = result.get_value(i); - ScalarType time_i = result.get_time(i); - + + // Access data at a specific time point, e.g. i=0. + size_t i = 0; + Eigen::VectorX value_at_time_i = result.get_value(i); + ScalarType time_i = result.get_time(i); + // Access the last time point. - Eigen::VectorX last_value = result.get_last_value(); - ScalarType last_time = result.get_last_time(); + Eigen::VectorX last_value = result.get_last_value(); + ScalarType last_time = result.get_last_time(); You can print the simulation results as a formatted table: @@ -448,7 +451,7 @@ You can print the simulation results as a formatted table: // Print results to console with default formatting. result.print_table(); - + // Print with custom column labels. std::vector labels = {"S", "E", "C", "I", "H", "U", "R", "D"}; result.print_table(labels); @@ -458,7 +461,7 @@ Additionally, you can export the results to a CSV file: .. code-block:: cpp // Export results to CSV with default settings. - result.export_csv("simulation_results.csv"); + auto export_status = result.export_csv("simulation_results.csv"); Visualization diff --git a/docs/source/cpp/models/isecir.rst b/docs/source/cpp/models/isecir.rst index 9df8194094..ffae3b174e 100644 --- a/docs/source/cpp/models/isecir.rst +++ b/docs/source/cpp/models/isecir.rst @@ -138,14 +138,14 @@ Then we can define the initial flows as follows. .. code-block:: cpp - int num_transitions = (int)mio::isecir::InfectionTransition::Count; + size_t 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); - // Define vector with flows. - Vec vec_init(num_transitions * num_agegroups); + // Define vector with flows. + Eigen::VectorX 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; @@ -214,22 +214,24 @@ Basic dampings can be added to the contact matrix as follows: // Create a contact matrix with constant contact rates between all groups. ScalarType cont_freq = 10.; - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); // Add a uniform damping across all age groups. - contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.));; For age-resolved models, you can apply different dampings to different groups: .. code-block:: cpp ScalarType cont_freq = 10.; - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, cont_freq)); + contact_matrix[0] = + mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, cont_freq)); // Add a damping that reduces contacts within the same age group by 70% starting at day 30. - contact_matrix.add_damping(Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(), - mio::SimulationTime(30.)); + Eigen::MatrixX damping_matrix = Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(); + contact_matrix[0].add_damping(damping_matrix, mio::SimulationTime(30.)); + For more complex scenarios, such as real-world venue closures or lockdown modeling, you can implement detailed NPIs with location-specific dampings. The IDE-SECIR model supports contact matrices for different locations (e.g., home, school, work, other) and can apply different dampings to each location. @@ -325,12 +327,13 @@ You can access the data in the `TimeSeries` objects as follows: // Get the number of time points. auto num_points = static_cast(compartments.get_num_time_points()); - // Access data at a specific time point. - Eigen::VectorX value_at_time_i = compartments.get_value(i); - ScalarType time_i = compartments.get_time(i); + // Access data at a specific time point, e.g. i=0. + size_t i = 0; + Eigen::VectorX value_at_time_i = compartments.get_value(i); + ScalarType time_i // Access the last time point. - Eigen::VectorX last_value = compartments.get_last_value(); + Eigen::VectorX last_value = compartments.get_last_value(); ScalarType last_time = compartments.get_last_time(); @@ -350,7 +353,7 @@ Additionally, you can export the results to a CSV file: .. code-block:: cpp // Export results to CSV with default settings. - compartments.export_csv("simulation_results.csv"); + auto export_status = compartments.export_csv("simulation_results.csv"); Visualization diff --git a/docs/source/cpp/models/iseir.rst b/docs/source/cpp/models/iseir.rst index d1b6714848..8b8ecaf384 100644 --- a/docs/source/cpp/models/iseir.rst +++ b/docs/source/cpp/models/iseir.rst @@ -68,21 +68,20 @@ The initialization of the model can be done as follows where we set the `Suscept .. code-block:: cpp - using Vec = mio::TimeSeries::Vector; - + using Vector = Eigen::Matrix; int N = 810000; - double dt = 0.1; - mio::TimeSeries init(1); + ScalarType dt = 0.1; + mio::TimeSeries init(1); - init.add_time_point(-15.0, Vec::Constant(1, N * 0.95)); + init.add_time_point>(-15.0, Vector::Constant(1, N * 0.95)); while (init.get_last_time() < 0) { init.add_time_point(init.get_last_time() + dt, - Vec::Constant(1, (double)init.get_last_value()[0] + init.get_last_time())); + Vector::Constant(1, (ScalarType)init.get_last_value()[0] + init.get_last_time())); } // Initialize model. - mio::iseir::Model model(std::move(init), dt, N); + mio::iseir::Model model(std::move(init), dt, N); .. _Nonpharmaceutical Interventions: Nonpharmaceutical Interventions @@ -97,11 +96,11 @@ Basic dampings can be added to the contact matrix as follows: // Create a contact matrix with constant contact rates between all groups. ScalarType cont_freq = 10.; - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); // Add a damping that reduces contacts by 70% starting at day 30. - contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); Simulation @@ -133,12 +132,13 @@ points. You can access the results as follows: // Get the number of time points. auto num_points = static_cast(result.get_num_time_points()); - // Access data at a specific time point. - Eigen::VectorX value_at_time_i = result.get_value(i); + // Access data at a specific time point, e.g. i=0. + size_t i = 0; + Eigen::VectorX value_at_time_i = result.get_value(i); ScalarType time_i = result.get_time(i); // Access the last time point. - Eigen::VectorX last_value = result.get_last_value(); + Eigen::VectorX last_value = result.get_last_value(); ScalarType last_time = result.get_last_time(); The order of the compartments follows the definition in the `InfectionState` enum. @@ -159,7 +159,7 @@ Additionally, you can export the results to a CSV file: .. code-block:: cpp // Export results to CSV with default settings. - result.export_csv("simulation_results.csv"); + auto result_status = result.export_csv("simulation_results.csv"); Visualization diff --git a/docs/source/cpp/models/lsecir.rst b/docs/source/cpp/models/lsecir.rst index 3ccbc98152..9d43933a59 100644 --- a/docs/source/cpp/models/lsecir.rst +++ b/docs/source/cpp/models/lsecir.rst @@ -135,9 +135,9 @@ To initialize the model, we start by defining the number of subcompartments and constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 3, NumInfectedSymptoms = 1, NumInfectedSevere = 1, NumInfectedCritical = 5; using InfState = mio::lsecir::InfectionState; - using LctState = mio::LctInfectionState; - using Model = mio::lsecir::Model; + using LctState = mio::LctInfectionState; + using Model = mio::lsecir::Model; Model model; For the simulation, we need initial values for all (sub)compartments. If we do not set the initial values manually, these are set to :math:`0` by default. @@ -146,46 +146,42 @@ We start with constructing a vector ``initial_populations`` that we will pass on .. code-block:: cpp - std::vector> initial_populations = {{750}, {30, 20}, {20, 10, 10}, {50}, - {50}, {10, 10, 5, 3, 2}, {20}, {10}}; + std::vector> initial_populations = {{750}, {30, 20}, {20, 10, 10}, {50}, + {50}, {10, 10, 5, 3, 2}, {20}, {10}}; We assert that the vector has the correct size by checking that the number of `InfectionState`\s and the numbers of subcompartments are correct. .. code-block:: cpp - if (initial_populations.size() != (size_t)InfState::Count) { - mio::log_error( - "The number of vectors in initial_populations does not match the number of InfectionStates."); - return 1; - } - if ((initial_populations[(size_t)InfState::Susceptible].size() != - LctState::get_num_subcompartments()) || - (initial_populations[(size_t)InfState::Exposed].size() != NumExposed) || - (initial_populations[(size_t)InfState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) || - (initial_populations[(size_t)InfState::InfectedSymptoms].size() != NumInfectedSymptoms) || - (initial_populations[(size_t)InfState::InfectedSevere].size() != NumInfectedSevere) || - (initial_populations[(size_t)InfState::InfectedCritical].size() != NumInfectedCritical) || - (initial_populations[(size_t)InfState::Recovered].size() != - LctState::get_num_subcompartments()) || - (initial_populations[(size_t)InfState::Dead].size() != - LctState::get_num_subcompartments())) { - mio::log_error( - "The length of at least one vector in initial_populations does not match the related number of " - "subcompartments."); - return 1; - } + if (initial_populations.size() != (size_t)InfState::Count) { + mio::log_error("The number of vectors in initial_populations does not match the number of InfectionStates."); + return 1; + } + if ((initial_populations[(size_t)InfState::Susceptible].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(size_t)InfState::Exposed].size() != NumExposed) || + (initial_populations[(size_t)InfState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) || + (initial_populations[(size_t)InfState::InfectedSymptoms].size() != NumInfectedSymptoms) || + (initial_populations[(size_t)InfState::InfectedSevere].size() != NumInfectedSevere) || + (initial_populations[(size_t)InfState::InfectedCritical].size() != NumInfectedCritical) || + (initial_populations[(size_t)InfState::Recovered].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(size_t)InfState::Dead].size() != LctState::get_num_subcompartments())) { + mio::log_error("The length of at least one vector in initial_populations does not match the related number of " + "subcompartments."); + return 1; + } The initial populations in the model are set via: .. code-block:: cpp - std::vector flat_initial_populations; - for (auto&& vec : initial_populations) { - flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end()); - } - for (size_t i = 0; i < LctState::Count; i++) { - model.populations[i] = flat_initial_populations[i]; - } + std::vector flat_initial_populations; + for (auto&& vec : initial_populations) { + flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end()); + } + for (size_t i = 0; i < LctState::Count; i++) { + model.populations[i] = flat_initial_populations[i]; } @@ -207,22 +203,24 @@ Basic dampings can be added to the contact matrix as follows: // Create a contact matrix with constant contact rates between all groups. ScalarType cont_freq = 10.; - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); - + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); + // Add a uniform damping across all age groups. - contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); For age-resolved models, you can apply different dampings to different groups: .. code-block:: cpp ScalarType cont_freq = 10.; - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, cont_freq)); - + contact_matrix[0] = + mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, cont_freq)); + // Add a damping that reduces contacts within the same age group by 70% starting at day 30. - contact_matrix.add_damping(Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(), - mio::SimulationTime(30.)); + Eigen::MatrixX damping_matrix = Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(); + contact_matrix[0].add_damping(damping_matrix, mio::SimulationTime(30.)); For more complex scenarios, such as real-world venue closures or lockdown modeling, you can implement detailed NPIs with location-specific dampings. The LCT-SECIR model supports contact matrices for different locations (e.g., home, school, work, other) and can apply different dampings to each location. @@ -280,22 +278,22 @@ We can simulate the model from :math:`t_0` to :math:`t_{\max}` with initial step .. code-block:: cpp - ScalarType t0 = 0; - ScalarType tmax = 10; - ScalarType dt = 0.5; + ScalarType t0 = 0; + ScalarType tmax = 10; + ScalarType dt = 0.5; mio::TimeSeries result = mio::simulate(t0, tmax, dt, model); You can also specify a custom integrator: .. code-block:: cpp - auto integrator = std::make_unique(); + auto integrator = std::make_unique>(); integrator->set_dt_min(0.3); integrator->set_dt_max(1.0); integrator->set_rel_tolerance(1e-4); integrator->set_abs_tolerance(1e-1); - - mio::TimeSeries result = mio::simulate(t0, tmax, dt, model, std::move(integrator)); + + result = mio::simulate(t0, tmax, dt, model, std::move(integrator)); Output @@ -313,14 +311,15 @@ You can access the data in the `mio::TimeSeries` object as follows: // Get the number of time points. auto num_points = static_cast(result.get_num_time_points()); - - // Access data at a specific time point. - Eigen::VectorX value_at_time_i = result.get_value(i); - ScalarType time_i = result.get_time(i); - + + // Access data at a specific time point, e.g. i=0. + size_t i = 0; + Eigen::VectorX value_at_time_i = result.get_value(i); + ScalarType time_i = result.get_time(i); + // Access the last time point. - Eigen::VectorX last_value = result.get_last_value(); - ScalarType last_time = result.get_last_time(); + Eigen::VectorX last_value = result.get_last_value(); + ScalarType last_time = result.get_last_time(); You can print the simulation results as a formatted table: @@ -339,7 +338,7 @@ Additionally, you can export the results to a CSV file: .. code-block:: cpp // Export results to CSV with default settings. - result.export_csv("simulation_results.csv"); + auto result_status = result.export_csv("simulation_results.csv"); Visualization diff --git a/docs/source/cpp/models/lsecir2d.rst b/docs/source/cpp/models/lsecir2d.rst index ad96e8cda6..1229434356 100644 --- a/docs/source/cpp/models/lsecir2d.rst +++ b/docs/source/cpp/models/lsecir2d.rst @@ -274,21 +274,26 @@ To initialize the model, we start by defining the number of subcompartments for .. code-block:: cpp - constexpr size_t NumExposed_1a = 1, NumInfectedNoSymptoms_1a = 1, NumInfectedSymptoms_1a = 1, - NumInfectedSevere_1a = 1, NumInfectedCritical_1a = 1, NumExposed_2a = 1, - NumInfectedNoSymptoms_2a = 1, NumInfectedSymptoms_2a = 1, NumInfectedSevere_2a = 1, - NumInfectedCritical_2a = 1, NumExposed_1b = 1, NumInfectedNoSymptoms_1b = 1, - NumInfectedSymptoms_1b = 1, NumInfectedSevere_1b = 1, NumInfectedCritical_1b = 1, - NumExposed_2b = 1, NumInfectedNoSymptoms_2b = 1, NumInfectedSymptoms_2b = 1, - NumInfectedSevere_2b = 1, NumInfectedCritical_2b = 1; - using InfState = mio::lsecir2d::InfectionState; - using LctState = mio::LctInfectionState< - InfState, 1, NumExposed_1a, NumInfectedNoSymptoms_1a, NumInfectedSymptoms_1a, NumInfectedSevere_1a, - NumInfectedCritical_1a, 1, 1, NumExposed_2a, NumInfectedNoSymptoms_2a, NumInfectedSymptoms_2a, - NumInfectedSevere_2a, NumInfectedCritical_2a, NumExposed_1b, NumInfectedNoSymptoms_1b, NumInfectedSymptoms_1b, - NumInfectedSevere_1b, NumInfectedCritical_1b, 1, 1, NumExposed_2b, NumInfectedNoSymptoms_2b, - NumInfectedSymptoms_2b, NumInfectedSevere_2b, NumInfectedCritical_2b, 1>; - using Model = mio::lsecir2d::Model; + constexpr size_t NumExposed_1a = 2, NumInfectedNoSymptoms_1a = 3, NumInfectedSymptoms_1a = 3, + NumInfectedSevere_1a = 3, NumInfectedCritical_1a = 2, NumExposed_2a = 1, + NumInfectedNoSymptoms_2a = 2, NumInfectedSymptoms_2a = 2, NumInfectedSevere_2a = 2, + NumInfectedCritical_2a = 1, NumExposed_1b = 2, NumInfectedNoSymptoms_1b = 3, + NumInfectedSymptoms_1b = 3, NumInfectedSevere_1b = 3, NumInfectedCritical_1b = 2, + NumExposed_2b = 1, NumInfectedNoSymptoms_2b = 2, NumInfectedSymptoms_2b = 2, + NumInfectedSevere_2b = 2, NumInfectedCritical_2b = 1; + // The compartment for Susceptible people and all compartments for Dead and Recovered people must have exactly one subcompartment: + constexpr size_t NumSusceptible = 1, NumDead_a = 1, NumDead_b = 1, NumRecovered_1a = 1, NumRecovered_1b = 1, + NumRecovered_2ab = 1; + using InfState = mio::lsecir2d::InfectionState; + using LctState = + mio::LctInfectionState; + using Model = mio::lsecir2d::Model; Model model; For the simulation, we need initial values for all (sub)compartments. If we do not set the initial values manually, these are set to :math:`0` by default. @@ -298,7 +303,7 @@ that contains a vector with initial values for the respective subcompartments. .. code-block:: cpp - std::vector> initial_populations = { + std::vector> initial_populations = { {200}, {0, 0}, {30, 10, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0}, {0}, {0}, {0}, {0, 0}, {10, 0}, {0, 0}, {0}, {10, 0}, {30, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0}, {0}, {0}, {100}, {0, 0}, {0, 0}, {0, 0}, {0}, {0}}; @@ -324,8 +329,8 @@ We assert that the vector has the correct size by checking that the number of `I (initial_populations[(size_t)InfState::InfectedSymptoms_2a].size() != NumInfectedSymptoms_2a) || (initial_populations[(size_t)InfState::InfectedSevere_2a].size() != NumInfectedSevere_2a) || (initial_populations[(size_t)InfState::InfectedCritical_2a].size() != NumInfectedCritical_2a) || - (initial_populations[(size_t)InfState::Recovered_a].size() != - LctState::get_num_subcompartments()) || + (initial_populations[(size_t)InfState::Recovered_1a].size() != + LctState::get_num_subcompartments()) || (initial_populations[(size_t)InfState::Dead_a].size() != LctState::get_num_subcompartments()) || (initial_populations[(size_t)InfState::Exposed_1b].size() != NumExposed_1b) || @@ -338,12 +343,12 @@ We assert that the vector has the correct size by checking that the number of `I (initial_populations[(size_t)InfState::InfectedSymptoms_2b].size() != NumInfectedSymptoms_2b) || (initial_populations[(size_t)InfState::InfectedSevere_2b].size() != NumInfectedSevere_2b) || (initial_populations[(size_t)InfState::InfectedCritical_2b].size() != NumInfectedCritical_2b) || - (initial_populations[(size_t)InfState::Recovered_ab].size() != - LctState::get_num_subcompartments())) { - mio::log_error("The length of at least one vector in initial_populations does not match the related number of " - "subcompartments."); - return 1; - } + (initial_populations[(size_t)InfState::Recovered_2ab].size() != + LctState::get_num_subcompartments())) { + mio::log_error("The length of at least one vector in initial_populations does not match the related number of " + "subcompartments."); + return 1; + } The initial populations in the model are set via: @@ -369,22 +374,23 @@ Basic dampings can be added to the contact matrix as follows: .. code-block:: cpp // Create a contact matrix with constant contact rate 10 (one age group). - mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10)); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, 10)); // From SimulationTime 5, the contact pattern is reduced to 30% of the initial value. - contact_matrix[0].add_damping(0.7, mio::SimulationTime(5.)); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(5.)); For age-resolved models, you can apply different dampings to different groups: .. code-block:: cpp // Create a contact matrix with constant contact rate 10 between all age groups - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, 10)); - + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, 10)); + // Add a damping that reduces contacts within the same age group by 70% starting at day 5. - contact_matrix.add_damping(Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(), - mio::SimulationTime(5.)); + Eigen::MatrixX damping_matrix = Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(); + contact_matrix.add_damping(damping_matrix, mio::SimulationTime(5.)); Simulation ---------- @@ -393,12 +399,11 @@ We can simulate the model from :math:`t_0` to :math:`t_{\max}` with initial step .. code-block:: cpp - ScalarType t0 = 0; - ScalarType tmax = 10; - ScalarType dt = 0.5; + ScalarType t0 = 0; + ScalarType tmax = 10; + ScalarType dt = 0.5; mio::TimeSeries result = mio::simulate(t0, tmax, dt, model); - Output ------ @@ -414,14 +419,15 @@ You can access the data in the `mio::TimeSeries` object as follows: // Get the number of time points. auto num_points = static_cast(result.get_num_time_points()); - - // Access data at a specific time point. - Eigen::VectorX value_at_time_i = result.get_value(i); - ScalarType time_i = result.get_time(i); - + + // Access data at a specific time point, e.g. i=0. + size_t i = 0; + Eigen::VectorX value_at_time_i = result.get_value(i); + ScalarType time_i = result.get_time(i); + // Access the last time point. - Eigen::VectorX last_value = result.get_last_value(); - ScalarType last_time = result.get_last_time(); + Eigen::VectorX last_value = result.get_last_value(); + ScalarType last_time = result.get_last_time(); You can print the simulation results as a formatted table: @@ -430,20 +436,19 @@ You can print the simulation results as a formatted table: // Print results to console with default formatting. result.print_table(); - + // Print with custom column labels and width, and custom number of decimals. - results.print_table({" S", " E1a", " C1a", " I1a", " H1a", " U1a", " Ra", - " Da", " E2a", " C2a", " I2a", " H2a", " U2a", " E1b", - " C1b", " I1b", " H1b", " U1b", " Rb", " Db", " E2b", - " C2b", " I2b", " H2b", " U2b", " Rab"}, - 6, 2); + result.print_table({" S", " E1a", " C1a", " I1a", " H1a", " U1a", " R1a", " Da", " E2a", + " C2a", " I2a", " H2a", " U2a", " E1b", " C1b", " I1b", " H1b", " U1b", + " R1b", " Db", " E2b", " C2b", " I2b", " H2b", " U2b", " R2ab"}, + 6, 2); Additionally, you can export the results to a CSV file: .. code-block:: cpp // Export results to CSV with default settings. - result.export_csv("simulation_results.csv"); + auto export_status = result.export_csv("simulation_results.csv"); Visualization From 5785fdf3aa4581405e7667eaec70e67c425cd69f Mon Sep 17 00:00:00 2001 From: Anna Wendler <106674756+annawendler@users.noreply.github.com> Date: Thu, 7 May 2026 16:08:45 +0200 Subject: [PATCH 03/12] review suggestions for ide and lct --- docs/source/cpp/models/isecir.rst | 4 ++-- docs/source/cpp/models/lsecir.rst | 2 +- docs/source/cpp/models/lsecir2d.rst | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/cpp/models/isecir.rst b/docs/source/cpp/models/isecir.rst index b07ea48317..650152d328 100644 --- a/docs/source/cpp/models/isecir.rst +++ b/docs/source/cpp/models/isecir.rst @@ -220,7 +220,7 @@ Basic dampings can be added to the contact matrix as follows: contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); // Add a uniform damping across all age groups. - contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.));; + contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); For age-resolved models, you can apply different dampings to different groups: @@ -332,7 +332,7 @@ You can access the data in the ``TimeSeries`` objects as follows: // Access data at a specific time point, e.g. i=0. size_t i = 0; Eigen::VectorX value_at_time_i = compartments.get_value(i); - ScalarType time_i + ScalarType time_i = compartments.get_time(i); // Access the last time point. Eigen::VectorX last_value = compartments.get_last_value(); diff --git a/docs/source/cpp/models/lsecir.rst b/docs/source/cpp/models/lsecir.rst index cd57b4e067..766c264515 100644 --- a/docs/source/cpp/models/lsecir.rst +++ b/docs/source/cpp/models/lsecir.rst @@ -216,7 +216,7 @@ For age-resolved models, you can apply different dampings to different groups: ScalarType cont_freq = 10.; contact_matrix[0] = - mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, cont_freq)); + mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, cont_freq)); // Add a damping that reduces contacts within the same age group by 70% starting at day 30. Eigen::MatrixX damping_matrix = Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(); diff --git a/docs/source/cpp/models/lsecir2d.rst b/docs/source/cpp/models/lsecir2d.rst index 85fc325a09..852658aafd 100644 --- a/docs/source/cpp/models/lsecir2d.rst +++ b/docs/source/cpp/models/lsecir2d.rst @@ -386,7 +386,7 @@ For age-resolved models, you can apply different dampings to different groups: .. code-block:: cpp // Create a contact matrix with constant contact rate 10 between all age groups - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, 10)); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, 10)); // Add a damping that reduces contacts within the same age group by 70% starting at day 5. Eigen::MatrixX damping_matrix = Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(); From 5152663184e7f3838fdb509d243c79f852d33e7d Mon Sep 17 00:00:00 2001 From: Anna Wendler <106674756+annawendler@users.noreply.github.com> Date: Thu, 7 May 2026 16:40:02 +0200 Subject: [PATCH 04/12] update docu for installing python packages --- docs/source/python/m-epidata.rst | 5 ++++ docs/source/python/m-generation.rst | 4 +++ docs/source/python/m-plot.rst | 5 ++++ docs/source/python/m-simulation.rst | 30 ++++------------------ docs/source/python/m-surrogate.rst | 5 ++++ docs/source/python/python_packages.rst | 35 ++++++++++++++++++++++++-- 6 files changed, 57 insertions(+), 27 deletions(-) diff --git a/docs/source/python/m-epidata.rst b/docs/source/python/m-epidata.rst index 086e6e7632..63199936bd 100644 --- a/docs/source/python/m-epidata.rst +++ b/docs/source/python/m-epidata.rst @@ -4,6 +4,11 @@ MEmilio Epidata MEmilio Epidata provides modules and scripts to download epidemiological data from various different sources. The package as well as links to the sources can be found in the `pycode/memilio-epidata `_. +Installation +------------ + +See :doc:`python_packages/Installation` for a detailed installation guide. + Dependencies ------------ diff --git a/docs/source/python/m-generation.rst b/docs/source/python/m-generation.rst index bcabac83f1..c1a47e4732 100644 --- a/docs/source/python/m-generation.rst +++ b/docs/source/python/m-generation.rst @@ -11,6 +11,10 @@ The following figure from Chapter 5 outlines the workflow of the generator. Blue .. image:: https://github.com/SciCompMod/memilio/raw/main/pycode/memilio-generation/generator_workflow.png :alt: Workflow of the generator +Installation +------------ + +See :doc:`python_packages/Installation` for a detailed installation guide. Dependencies ---------- diff --git a/docs/source/python/m-plot.rst b/docs/source/python/m-plot.rst index 1effe1003d..99aec4a7c1 100644 --- a/docs/source/python/m-plot.rst +++ b/docs/source/python/m-plot.rst @@ -8,6 +8,11 @@ The package is contained inside the folder `pycode/memilio-plot `_. + Installation ------------ -The ``memilio-simulation`` package can be installed in two ways depending on your use case. - -**Option 1: Install from PyPI (Recommended - no C++ compiler required)** - -Pre-built wheels are provided for Linux and Windows on Python 3.8 to 3.13. - -.. code-block:: console - - pip install memilio-simulation - -This is the easiest way to get started. No C++ compiler or CMake is needed. - -**Option 2: Install from source (latest development version, or contributing)** - -If you need the latest unreleased code, or want to modify the bindings, build from source. - -Requirements: +The easiest way to get started is installing the memilio-simulation package via PyPI. +If you need the latest unreleased code, or want to modify the bindings, build from source. In that case, the following +tools are required: * A **C++20 compiler** * **CMake** >= 3.18 * **Ninja** build tool * All :doc:`C++ library dependencies <../getting_started>` -Install from the **root of the MEmilio repository** (the directory containing the top-level ``pyproject.toml``): - -.. code-block:: console - - pip install -e .[dev] +See :doc:`python_packages/Installation` for a detailed installation guide. -This is necessary because the C++ build requires access to the ``cpp/`` directory. -Note that this only installs ``memilio-simulation`` and not any other Python packages.The ``-e`` flag links the installation to your local source code so Python changes are reflected immediately. -C++ changes require re-running this command to recompile. Dependencies ------------ diff --git a/docs/source/python/m-surrogate.rst b/docs/source/python/m-surrogate.rst index 229ab79096..47794b34bb 100644 --- a/docs/source/python/m-surrogate.rst +++ b/docs/source/python/m-surrogate.rst @@ -14,6 +14,11 @@ For more details, we refer to: |Graph_Neural_Network_Surrogates| +Installation +------------ + +See :doc:`python_packages/Installation` for a detailed installation guide. + Dependencies ------------ diff --git a/docs/source/python/python_packages.rst b/docs/source/python/python_packages.rst index 2a58a46ae2..6d02515287 100644 --- a/docs/source/python/python_packages.rst +++ b/docs/source/python/python_packages.rst @@ -108,12 +108,37 @@ Please see the individual package documentation for more details on the function .. _Python_Installation: + + Installation ------------ +The Python packages can be installed in two ways depending on your use case. + +**Option 1: Install from PyPI (Recommended - no C++ compiler required; currently only supported for memilio-simulation)** + +Pre-built wheels are provided for Linux and Windows on Python 3.8 to 3.13. + +.. code-block:: console + + pip install memilio-simulation + +This is the easiest way to get started. No C++ compiler or CMake is needed. + +**Option 2: Install from source (latest development version, or contributing)** + Each package provides a ``pyproject.toml`` that installs the package and its dependencies with pip. + +If you want to install the memilio-simulation package: + +* The pyproject.toml is in the **root of the MEmilio repository**. + +If you want to install any of the other Python packages: + +* The pyproject.toml is in the respective folder ``pycode/memilio-*`` + The dependencies of the individual packages are denoted in their documentation. -The installation can be run with the following command (from the directory containing the ``pyproject.toml`` file) +The installation can be run with the following command from the directory containing the ``pyproject.toml`` file .. code-block:: console @@ -127,7 +152,13 @@ For development of code use this command instead python -m pip install -e .[dev] -This command allows you to work on the code without having to reinstall the package after a change. It also installs all additional dependencies required for development and maintenance. +This command allows you to work on the code without having . It also installs +all additional dependencies required for development and maintenance. + + +The ``-e`` flag links the installation to your local source code so Python changes are reflected immediately. Hence, +you do not need to reinstall the package after a changes in Python. C++ changes require re-running this command to +recompile. .. dropdown:: :fa:`gears` Build files for skbuild From 3b9e3d898414d1c2c267faa463ee89b41fa76e00 Mon Sep 17 00:00:00 2001 From: Anna Wendler <106674756+annawendler@users.noreply.github.com> Date: Mon, 11 May 2026 14:39:15 +0200 Subject: [PATCH 05/12] review suggestions --- docs/source/cpp/models/osecir.rst | 4 ++-- docs/source/cpp/models/oseir.rst | 4 ++-- docs/source/cpp/models/oseirdb.rst | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/source/cpp/models/osecir.rst b/docs/source/cpp/models/osecir.rst index af58ed95b5..3d27bf3cc9 100644 --- a/docs/source/cpp/models/osecir.rst +++ b/docs/source/cpp/models/osecir.rst @@ -181,8 +181,8 @@ Basic dampings can be added to the contact matrix as follows: .. code-block:: cpp // Create a contact matrix with constant contact rates between all groups - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); // Add a damping that reduces contacts by 70% starting at day 30 contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); diff --git a/docs/source/cpp/models/oseir.rst b/docs/source/cpp/models/oseir.rst index fb410e4c61..33d6dc2222 100644 --- a/docs/source/cpp/models/oseir.rst +++ b/docs/source/cpp/models/oseir.rst @@ -131,8 +131,8 @@ Basic dampings can be added to the ContactPatterns as follows: .. code-block:: cpp // Create a contact matrix with constant contact rates between all groups - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); // Add a damping that reduces contacts by 70% starting at day 30 contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); diff --git a/docs/source/cpp/models/oseirdb.rst b/docs/source/cpp/models/oseirdb.rst index c2cdec51d9..1dc1b59835 100644 --- a/docs/source/cpp/models/oseirdb.rst +++ b/docs/source/cpp/models/oseirdb.rst @@ -140,8 +140,8 @@ Basic dampings can be added to the ContactPatterns as follows: .. code-block:: cpp // Create a contact matrix with constant contact rates between all groups - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixX::Constant(1, 1, cont_freq)); // Add a damping that reduces contacts by 70% starting at day 30 contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); From 38b194cec8932a3e017870b7ec0e2b7b075c86de Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Mon, 11 May 2026 15:14:35 +0200 Subject: [PATCH 06/12] CHG: Change double to Scalartype --- docs/source/cpp/glct.rst | 4 ++-- docs/source/cpp/graph_metapop.rst | 4 ++-- docs/source/cpp/lct.rst | 4 ++-- docs/source/cpp/mobility_based_abm.rst | 2 +- docs/source/cpp/models/osecirts.rst | 14 +++++++------- docs/source/cpp/ode.rst | 4 ++-- docs/source/cpp/ode_creation.rst | 6 +++--- docs/source/cpp/ode_metapop.rst | 2 +- docs/source/cpp/smm.rst | 10 +++++----- .../source/python/m-simulation_common_patterns.rst | 2 +- docs/source/python/m-simulation_model_usage.rst | 2 +- 11 files changed, 27 insertions(+), 27 deletions(-) diff --git a/docs/source/cpp/glct.rst b/docs/source/cpp/glct.rst index eba11ca5a7..413bfa38dd 100644 --- a/docs/source/cpp/glct.rst +++ b/docs/source/cpp/glct.rst @@ -45,8 +45,8 @@ In the ``ContactPatterns``, each matrix element stores baseline contact rates :m group :math:`i` to group :math:`j`. The dimension of the matrix is automatically defined by the model initialization and it is reduced to one value if no stratification is used. The values can be adjusted during the simulation, e.g., through implementing nonpharmaceutical interventions, see the section on :ref:`Nonpharmaceutical Interventions GLCT`. -Parameters can be accessed via ``model.parameters.get>()`` and set via either -``model.parameters.set>(value)`` or ``model.parameters.get>() = value``. +Parameters can be accessed via ``model.parameters.get>()`` and set via either +``model.parameters.set>(value)`` or ``model.parameters.get>() = value``. Initial conditions diff --git a/docs/source/cpp/graph_metapop.rst b/docs/source/cpp/graph_metapop.rst index b8f2deb230..c2a8746bed 100644 --- a/docs/source/cpp/graph_metapop.rst +++ b/docs/source/cpp/graph_metapop.rst @@ -163,7 +163,7 @@ The following steps detail how to configure and execute a graph simulation: .. code-block:: cpp - mio::Graph>>, mio::MobilityEdgeStochastic> graph; + mio::Graph>>, mio::MobilityEdgeStochastic> graph; graph.add_node(1001, model_group1, t0); graph.add_node(1002, model_group2, t0); @@ -185,7 +185,7 @@ The following steps detail how to configure and execute a graph simulation: .. code-block:: cpp - mio::GraphBuilder>>, mio::MobilityEdgeStochastic> builder; + mio::GraphBuilder>>, mio::MobilityEdgeStochastic> builder; builder.add_node(1001, model_group1, t0); builder.add_node(1002, model_group2, t0); builder.add_edge(0, 1, std::move(transition_rates)); diff --git a/docs/source/cpp/lct.rst b/docs/source/cpp/lct.rst index 18d330915c..b005364b47 100644 --- a/docs/source/cpp/lct.rst +++ b/docs/source/cpp/lct.rst @@ -54,8 +54,8 @@ In the ``ContactPatterns``, each matrix element stores baseline contact rates :m The dimension of the matrix is automatically defined by the model initialization and is reduced to one value if no stratification is used. The values can be adjusted during the simulation, e.g., through implementing nonpharmaceutical interventions, see the section on :ref:`Nonpharmaceutical Interventions LCT`. -Parameters can be accessed via ``model.parameters.get>()`` and set via either -``model.parameters.get>() = value`` or ``model.parameters.set>(value)``. +Parameters can be accessed via ``model.parameters.get>()`` and set via either +``model.parameters.get>() = value`` or ``model.parameters.set>(value)``. Initial conditions diff --git a/docs/source/cpp/mobility_based_abm.rst b/docs/source/cpp/mobility_based_abm.rst index 74ac617933..8ca8c659dc 100644 --- a/docs/source/cpp/mobility_based_abm.rst +++ b/docs/source/cpp/mobility_based_abm.rst @@ -299,7 +299,7 @@ For infections to happen during the simulation, we have to initialize people wit .. code-block:: cpp // Assign infection state to each person randomly with specific distribution - std::vector infection_distribution{0.5, 0.3, 0.05, 0.05, 0.05, 0.05, 0.0, 0.0}; + std::vector infection_distribution{0.5, 0.3, 0.05, 0.05, 0.05, 0.05, 0.0, 0.0}; for (auto& person : model.get_persons()) { mio::abm::InfectionState infection_state = mio::abm::InfectionState( mio::DiscreteDistribution::get_instance()(mio::thread_local_rng(), infection_distribution)); diff --git a/docs/source/cpp/models/osecirts.rst b/docs/source/cpp/models/osecirts.rst index e7d6680be6..ac91d67775 100644 --- a/docs/source/cpp/models/osecirts.rst +++ b/docs/source/cpp/models/osecirts.rst @@ -277,17 +277,17 @@ After setting the initial populations, the daily vaccination parameters, which a const size_t daily_vaccinations = 10; const size_t num_days = 300; - model.parameters.get>().resize(mio::SimulationDay(num_days)); - model.parameters.get>().resize(mio::SimulationDay(num_days)); - model.parameters.get>().resize(mio::SimulationDay(num_days)); + model.parameters.get>().resize(mio::SimulationDay(num_days)); + model.parameters.get>().resize(mio::SimulationDay(num_days)); + model.parameters.get>().resize(mio::SimulationDay(num_days)); for (size_t i = 0; i < num_days; ++i) { for (mio::AgeGroup j = 0; j < nb_groups; ++j) { - auto num_vaccinations = static_cast(i * daily_vaccinations); - model.parameters.get>()[{j, mio::SimulationDay(i)}] = + auto num_vaccinations = static_cast(i * daily_vaccinations); + model.parameters.get>()[{j, mio::SimulationDay(i)}] = num_vaccinations; - model.parameters.get>()[{j, mio::SimulationDay(i)}] = + model.parameters.get>()[{j, mio::SimulationDay(i)}] = num_vaccinations; - model.parameters.get>()[{j, mio::SimulationDay(i)}] = + model.parameters.get>()[{j, mio::SimulationDay(i)}] = num_vaccinations; } } diff --git a/docs/source/cpp/ode.rst b/docs/source/cpp/ode.rst index a316eb911e..37e4eb87c1 100644 --- a/docs/source/cpp/ode.rst +++ b/docs/source/cpp/ode.rst @@ -56,8 +56,8 @@ In the ``ContactPatterns`` parameter, each matrix element stores baseline contac group :math:`i` and group :math:`j`. The dimension of the matrix is automatically defined by the model initiation and it is reduced to one value if no stratifcation is used. The values can be adjusted during the simulation, e.g., through implementing nonpharmaceutical interventions, see the section on :ref:`Nonpharmaceutical Interventions`. -Parameters can get accessed via ``model.parameters.get>()`` and set via either -``model.parameters.get>() = value`` or ``model.parameters.set>(value)``. +Parameters can get accessed via ``model.parameters.get>()`` and set via either +``model.parameters.get>() = value`` or ``model.parameters.set>(value)``. .. dropdown:: :fa:`gears` Expert's settings diff --git a/docs/source/cpp/ode_creation.rst b/docs/source/cpp/ode_creation.rst index 5fa3d6301a..ab803077fe 100644 --- a/docs/source/cpp/ode_creation.rst +++ b/docs/source/cpp/ode_creation.rst @@ -102,9 +102,9 @@ and for the contact rate :math:`\phi` a struct Avoid using the mathematical symbols of the constant as names for the struct. Their connection can be noted in the documentation of these structs. -The template :code:`FP` and the type :code:`UncertainValue` in these examples are commonly used throughout MEmilio. -:code:`FP` is a floating point type, usually :code:`double`. An :code:`UncertainValue` holds a value of type -:code:`FP` as well as (optionally) a distribution to sample new values from, e.g. for a parameter study. +The template :code:`ScalarType` and the type :code:`UncertainValue` in these examples are commonly used throughout MEmilio. +:code:`ScalarType` is a floating point type, usually :code:`double`. An :code:`UncertainValue` holds a value of type +:code:`ScalarType` as well as (optionally) a distribution to sample new values from, e.g. for a parameter study. Finally, define a type ``Parameters`` by listing all parameter structs as template arguments of a :code:`mio::ParameterSet`: diff --git a/docs/source/cpp/ode_metapop.rst b/docs/source/cpp/ode_metapop.rst index 3b276da4d8..ccbefab586 100644 --- a/docs/source/cpp/ode_metapop.rst +++ b/docs/source/cpp/ode_metapop.rst @@ -24,7 +24,7 @@ To set up a simulation of the ODE metapopulation model, you need to initialize t .. code-block:: cpp - mio::oseirmetapop::Model model(3, 1) + mio::oseirmetapop::Model model(3, 1) Set a population with the number of individuals in each region, age group and epidemiological compartment, e.g.: diff --git a/docs/source/cpp/smm.rst b/docs/source/cpp/smm.rst index 334ea4e819..38f65e30af 100644 --- a/docs/source/cpp/smm.rst +++ b/docs/source/cpp/smm.rst @@ -129,7 +129,7 @@ These populations have the class type ``Populations`` and can be set via: .. code-block:: cpp - double pop = 1000, numE = 0.001 * pop, numC = 0.0001 * pop, numI = 0.0001 * pop, numR = 0, numD = 0; + ScalarType pop = 1000, numE = 0.001 * pop, numC = 0.0001 * pop, numI = 0.0001 * pop, numR = 0, numD = 0; //Population is distributed equally to the regions for (size_t r = 0; r < num_regions; ++r) { @@ -184,9 +184,9 @@ until ``tmax``. The step size is only used to regularly save the system state du .. code-block:: cpp - double t0 = 0.0; - double dt = 0.1; - double tmax = 30.; + ScalarType t0 = 0.0; + ScalarType dt = 0.1; + ScalarType tmax = 30.; //Pass the model, t0 and dt to the Simulation auto sim = mio::smm::Simulation(model, t0, dt); @@ -250,7 +250,7 @@ We can define a model: Now, for accessing the population, all indices need to be given: --- code-block:: cpp +.. code-block:: cpp model.populations[{Region(r), InfectionState::S, Species(0)}] = 100; // ... diff --git a/docs/source/python/m-simulation_common_patterns.rst b/docs/source/python/m-simulation_common_patterns.rst index 8bb10cf38c..781508d69a 100644 --- a/docs/source/python/m-simulation_common_patterns.rst +++ b/docs/source/python/m-simulation_common_patterns.rst @@ -101,7 +101,7 @@ MEmilio uses ``mio::IOResult`` for handling return values and errors of io funct .. code-block:: c++ - mio::IOResult save_result(const std::vector>& results, const std::vector& ids, int num_groups, + mio::IOResult save_result(const std::vector>& results, const std::vector& ids, int num_groups, const std::string& filename) { ...; diff --git a/docs/source/python/m-simulation_model_usage.rst b/docs/source/python/m-simulation_model_usage.rst index 5daef8d840..feccd22b3c 100755 --- a/docs/source/python/m-simulation_model_usage.rst +++ b/docs/source/python/m-simulation_model_usage.rst @@ -40,7 +40,7 @@ Python and C++ to better understand the differences of both interfaces. #include "ode_seir/model.h" int num_groups = 1; - mio::oseir::Model model(num_groups); + mio::oseir::Model model(num_groups); Initialize parameters --------------------- From de8d8966c3e3840fd2916243df0a21037e62e4dd Mon Sep 17 00:00:00 2001 From: annawendler <106674756+annawendler@users.noreply.github.com> Date: Mon, 11 May 2026 15:22:19 +0200 Subject: [PATCH 07/12] Update docs/source/python/python_packages.rst Co-authored-by: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> --- docs/source/python/python_packages.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/python/python_packages.rst b/docs/source/python/python_packages.rst index 6d02515287..1d05d88615 100644 --- a/docs/source/python/python_packages.rst +++ b/docs/source/python/python_packages.rst @@ -152,7 +152,7 @@ For development of code use this command instead python -m pip install -e .[dev] -This command allows you to work on the code without having . It also installs +This command allows you to work on the code without having to reinstall the package after a change. It also installs all additional dependencies required for development and maintenance. From 49286575d497234948fbc753525ed3f40e5a6aa8 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Mon, 11 May 2026 15:24:45 +0200 Subject: [PATCH 08/12] CHG: Use FP in backend --- docs/source/cpp/ode_creation.rst | 6 +++--- docs/source/python/m-simulation_common_patterns.rst | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/cpp/ode_creation.rst b/docs/source/cpp/ode_creation.rst index ab803077fe..5fa3d6301a 100644 --- a/docs/source/cpp/ode_creation.rst +++ b/docs/source/cpp/ode_creation.rst @@ -102,9 +102,9 @@ and for the contact rate :math:`\phi` a struct Avoid using the mathematical symbols of the constant as names for the struct. Their connection can be noted in the documentation of these structs. -The template :code:`ScalarType` and the type :code:`UncertainValue` in these examples are commonly used throughout MEmilio. -:code:`ScalarType` is a floating point type, usually :code:`double`. An :code:`UncertainValue` holds a value of type -:code:`ScalarType` as well as (optionally) a distribution to sample new values from, e.g. for a parameter study. +The template :code:`FP` and the type :code:`UncertainValue` in these examples are commonly used throughout MEmilio. +:code:`FP` is a floating point type, usually :code:`double`. An :code:`UncertainValue` holds a value of type +:code:`FP` as well as (optionally) a distribution to sample new values from, e.g. for a parameter study. Finally, define a type ``Parameters`` by listing all parameter structs as template arguments of a :code:`mio::ParameterSet`: diff --git a/docs/source/python/m-simulation_common_patterns.rst b/docs/source/python/m-simulation_common_patterns.rst index 781508d69a..3e8998da6d 100644 --- a/docs/source/python/m-simulation_common_patterns.rst +++ b/docs/source/python/m-simulation_common_patterns.rst @@ -101,7 +101,7 @@ MEmilio uses ``mio::IOResult`` for handling return values and errors of io funct .. code-block:: c++ - mio::IOResult save_result(const std::vector>& results, const std::vector& ids, int num_groups, + mio::IOResult save_result(const std::vector>& results, const std::vector& ids, int num_groups, const std::string& filename) { ...; From bb5380f3f9532a7945e1cfb6e5b90d86aefb8ef5 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Mon, 11 May 2026 15:43:24 +0200 Subject: [PATCH 09/12] CHG: Use ScalarType instead of double in feedback examples --- cpp/examples/ode_secir_feedback.cpp | 65 +++++++++++---------- cpp/examples/ode_secir_feedback_graph.cpp | 69 ++++++++++++----------- 2 files changed, 69 insertions(+), 65 deletions(-) diff --git a/cpp/examples/ode_secir_feedback.cpp b/cpp/examples/ode_secir_feedback.cpp index 61876d6205..650a94dd60 100644 --- a/cpp/examples/ode_secir_feedback.cpp +++ b/cpp/examples/ode_secir_feedback.cpp @@ -21,32 +21,33 @@ #include "memilio/compartments/feedback_simulation.h" #include "memilio/utils/logging.h" -void initialize_model(mio::osecir::Model& model, int total_population, double cont_freq) +void initialize_model(mio::osecir::Model& model, int total_population, ScalarType cont_freq) { - model.parameters.set>(60); - model.parameters.set>(0.2); + model.parameters.set>(60); + model.parameters.set>(0.2); // time-related parameters - model.parameters.get>() = 3.2; - model.parameters.get>() = 2.0; - model.parameters.get>() = 5.8; - model.parameters.get>() = 9.5; - model.parameters.get>() = 7.1; + model.parameters.get>() = 3.2; + model.parameters.get>() = 2.0; + model.parameters.get>() = 5.8; + model.parameters.get>() = 9.5; + model.parameters.get>() = 7.1; // Set transmission and isolation parameters - model.parameters.get>() = 0.05; - model.parameters.get>() = 0.7; - model.parameters.get>() = 0.09; - model.parameters.get>() = 0.25; - model.parameters.get>() = 0.45; - model.parameters.get>() = 35; - model.parameters.get>() = 0.2; - model.parameters.get>() = 0.25; - model.parameters.get>() = 0.3; + model.parameters.get>() = 0.05; + model.parameters.get>() = 0.7; + model.parameters.get>() = 0.09; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.45; + model.parameters.get>() = 35; + model.parameters.get>() = 0.2; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.3; // contact matrix - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); // initial population model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = 40; @@ -65,14 +66,15 @@ void initialize_model(mio::osecir::Model& model, int total_population, d model.apply_constraints(); } -void initialize_feedback(mio::FeedbackSimulation>, - mio::osecir::ContactPatterns>& feedback_simulation) +void initialize_feedback( + mio::FeedbackSimulation>, + mio::osecir::ContactPatterns>& feedback_simulation) { // nominal ICU capacity - feedback_simulation.get_parameters().template get>() = 10; + feedback_simulation.get_parameters().template get>() = 10; // ICU occupancy in the past for memory kernel - auto& icu_occupancy = feedback_simulation.get_parameters().template get>(); + auto& icu_occupancy = feedback_simulation.get_parameters().template get>(); Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1); const auto cutoff = static_cast(feedback_simulation.get_parameters().template get()); for (int t = -cutoff; t <= 0; ++t) { @@ -80,8 +82,8 @@ void initialize_feedback(mio::FeedbackSimulation>() = {0.1}; - feedback_simulation.get_parameters().template get>() = {0.8}; + feedback_simulation.get_parameters().template get>() = {0.1}; + feedback_simulation.get_parameters().template get>() = {0.8}; } int main() @@ -92,12 +94,12 @@ int main() // risk which is calculated from the ICU occupancy using a memory kernel. mio::set_log_level(mio::LogLevel::warn); - const double tmax = 35; + const ScalarType tmax = 35; const int total_population = 1000; - const double cont_freq = 10; + const ScalarType cont_freq = 10; // create and initialize ODE model for a single age group - mio::osecir::Model model(1); + mio::osecir::Model model(1); initialize_model(model, total_population, cont_freq); // determine the index for the ICU state (InfectedCritical) for feedback mechanism @@ -105,10 +107,11 @@ int main() model.populations.get_flat_index({mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical})}; // create simulation objects: first a secir simulation, then a feedback simulation - auto simulation = mio::osecir::Simulation>>(model); + auto simulation = + mio::osecir::Simulation>>(model); auto feedback_simulation = - mio::FeedbackSimulation>, - mio::osecir::ContactPatterns>(std::move(simulation), icu_index); + mio::FeedbackSimulation>, + mio::osecir::ContactPatterns>(std::move(simulation), icu_index); // set up the parameters for the feedback simulation initialize_feedback(feedback_simulation); diff --git a/cpp/examples/ode_secir_feedback_graph.cpp b/cpp/examples/ode_secir_feedback_graph.cpp index 0b3572d0d4..a73a9917c6 100644 --- a/cpp/examples/ode_secir_feedback_graph.cpp +++ b/cpp/examples/ode_secir_feedback_graph.cpp @@ -26,46 +26,47 @@ #include // alias for the type of the simulation with feedback -using FeedbackSim = mio::FeedbackSimulation>, - mio::osecir::ContactPatterns>; +using FeedbackSim = mio::FeedbackSimulation>, + mio::osecir::ContactPatterns>; // helper function to initialize the model with population and parameters -void initialize_model(mio::osecir::Model& model, double cont_freq) +void initialize_model(mio::osecir::Model& model, ScalarType cont_freq) { - model.parameters.set>(60); - model.parameters.set>(0.2); + model.parameters.set>(60); + model.parameters.set>(0.2); // Mean stay times per compartment - model.parameters.get>() = 3.2; - model.parameters.get>() = 2.0; - model.parameters.get>() = 5.8; - model.parameters.get>() = 9.5; - model.parameters.get>() = 7.1; + model.parameters.get>() = 3.2; + model.parameters.get>() = 2.0; + model.parameters.get>() = 5.8; + model.parameters.get>() = 9.5; + model.parameters.get>() = 7.1; // Set transmission and isolation parameters - model.parameters.get>() = 0.05; - model.parameters.get>() = 0.7; - model.parameters.get>() = 0.09; - model.parameters.get>() = 0.25; - model.parameters.get>() = 0.45; - model.parameters.get>() = 35; - model.parameters.get>() = 0.2; - model.parameters.get>() = 0.25; - model.parameters.get>() = 0.3; + model.parameters.get>() = 0.05; + model.parameters.get>() = 0.7; + model.parameters.get>() = 0.09; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.45; + model.parameters.get>() = 35; + model.parameters.get>() = 0.2; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.3; // contact matrix - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq)); } // helper function to initialize the feedback mechanism parameters for a simulation void initialize_feedback(FeedbackSim& feedback_simulation) { // nominal ICU capacity - feedback_simulation.get_parameters().template get>() = 10; + feedback_simulation.get_parameters().template get>() = 10; // ICU occupancy in the past for memory kernel - auto& icu_occupancy = feedback_simulation.get_parameters().template get>(); + auto& icu_occupancy = feedback_simulation.get_parameters().template get>(); Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1); const auto cutoff = static_cast(feedback_simulation.get_parameters().template get()); for (int t = -cutoff; t <= 0; ++t) { @@ -73,24 +74,24 @@ void initialize_feedback(FeedbackSim& feedback_simulation) } // bounds for contact reduction measures - feedback_simulation.get_parameters().template get>() = {0.1}; - feedback_simulation.get_parameters().template get>() = {0.8}; + feedback_simulation.get_parameters().template get>() = {0.1}; + feedback_simulation.get_parameters().template get>() = {0.8}; // Set blending factors. The global blending factor is implicitly defined as 1 - local - regional. - feedback_simulation.get_parameters().template get>() = 0.5; - feedback_simulation.get_parameters().template get>() = 0.3; + feedback_simulation.get_parameters().template get>() = 0.5; + feedback_simulation.get_parameters().template get>() = 0.3; } // helper function to create the graph with nodes and edges -mio::Graph, mio::MobilityEdge> -create_graph(int num_nodes, int total_population, double cont_freq) +mio::Graph, mio::MobilityEdge> +create_graph(int num_nodes, int total_population, ScalarType cont_freq) { // Create a graph for the metapopulation simulation - mio::Graph, mio::MobilityEdge> g; + mio::Graph, mio::MobilityEdge> g; // Create models and add nodes to the graph for (int i = 0; i < num_nodes; ++i) { - mio::osecir::Model model(1); + mio::osecir::Model model(1); initialize_model(model, cont_freq); // Set initial populations (infection starts in the first node) @@ -120,7 +121,7 @@ create_graph(int num_nodes, int total_population, double cont_freq) model.populations.get_flat_index({mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical})}; // Create feedback simulation - auto feedback_sim = FeedbackSim(mio::Simulation>(model), icu_index); + auto feedback_sim = FeedbackSim(mio::Simulation>(model), icu_index); initialize_feedback(feedback_sim); // Node-ID-Logic: 1001-1005, 2001-2005, ... @@ -153,7 +154,7 @@ int main() const auto tmax = 10.; const auto dt = 0.5; const int total_population = 1000; - const double cont_freq = 2.7; + const ScalarType cont_freq = 2.7; const int num_nodes = 10; // Create the graph @@ -161,7 +162,7 @@ int main() // Create and run the simulation using Graph = decltype(g); - auto sim = mio::FeedbackGraphSimulation(std::move(g), t0, dt); + auto sim = mio::FeedbackGraphSimulation(std::move(g), t0, dt); sim.advance(tmax); // The output shows the compartments sizes for a node without any initial infections. From 42dac6823fe407a84617dfca8174881f28b345d1 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Mon, 11 May 2026 15:46:02 +0200 Subject: [PATCH 10/12] CHG: ScalarType to double --- cpp/examples/sde_sir.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/cpp/examples/sde_sir.cpp b/cpp/examples/sde_sir.cpp index 66746c4e2f..9866dafa47 100644 --- a/cpp/examples/sde_sir.cpp +++ b/cpp/examples/sde_sir.cpp @@ -25,15 +25,15 @@ int main() { mio::set_log_level(mio::LogLevel::debug); - double t0 = 0.; - double tmax = 5.; - double dt = 0.1; + ScalarType t0 = 0.; + ScalarType tmax = 5.; + ScalarType dt = 0.1; - double total_population = 10000; + ScalarType total_population = 10000; mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt); - mio::ssir::Model model; + mio::ssir::Model model; model.populations[{mio::Index(mio::ssir::InfectionState::Infected)}] = 100; model.populations[{mio::Index(mio::ssir::InfectionState::Recovered)}] = 1000; @@ -41,10 +41,11 @@ int main() total_population - model.populations[{mio::Index(mio::ssir::InfectionState::Infected)}] - model.populations[{mio::Index(mio::ssir::InfectionState::Recovered)}]; - model.parameters.set>(10); - model.parameters.set>(1); - model.parameters.get>().get_baseline()(0, 0) = 2.7; - model.parameters.get>().add_damping(0.6, mio::SimulationTime(12.5)); + model.parameters.set>(10); + model.parameters.set>(1); + model.parameters.get>().get_baseline()(0, 0) = 2.7; + model.parameters.get>().add_damping(0.6, + mio::SimulationTime(12.5)); model.check_constraints(); From 5098e3bd8d0b0ecb491091dae8f3405585a04492 Mon Sep 17 00:00:00 2001 From: jubicker <113909589+jubicker@users.noreply.github.com> Date: Mon, 11 May 2026 15:54:21 +0200 Subject: [PATCH 11/12] Replace double by ScalarType --- cpp/examples/d_abm.cpp | 13 +++++++------ docs/source/cpp/diffusive_abm.rst | 14 +++++++------- docs/source/cpp/temporal_hybrid.rst | 20 ++++++++++---------- 3 files changed, 24 insertions(+), 23 deletions(-) diff --git a/cpp/examples/d_abm.cpp b/cpp/examples/d_abm.cpp index 35ee9c16f2..68db33a561 100644 --- a/cpp/examples/d_abm.cpp +++ b/cpp/examples/d_abm.cpp @@ -21,6 +21,7 @@ #include "d_abm/quad_well.h" #include "d_abm/simulation.h" #include "d_abm/parameters.h" +#include "memilio/config.h" #include "memilio/utils/random_number_generator.h" #include "memilio/data/analyze_result.h" #include "memilio/epidemiology/adoption_rate.h" @@ -44,10 +45,10 @@ int main() using Model = mio::dabm::Model>; std::vector agents(1000); //Random variables for initialization of agents' position and infection state - auto& pos_rng = mio::UniformDistribution::get_instance(); + auto& pos_rng = mio::UniformDistribution::get_instance(); auto& sta_rng = mio::DiscreteDistribution::get_instance(); //Infection state distribution - std::vector pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.}; + std::vector pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.}; for (auto& a : agents) { //Agents are equally distributed in [-2,2]x[-2,2] at the beginning a.position = @@ -71,11 +72,11 @@ int main() } //Set interaction radius and noise term of the diffusion process - double interaction_radius = 0.5; - double noise = 0.4; + ScalarType interaction_radius = 0.5; + ScalarType noise = 0.4; - double dt = 0.1; - double tmax = 30.; + ScalarType dt = 0.1; + ScalarType tmax = 30.; Model model(agents, adoption_rates, interaction_radius, noise, {InfectionState::D}); auto sim = mio::dabm::Simulation(model, 0.0, dt); diff --git a/docs/source/cpp/diffusive_abm.rst b/docs/source/cpp/diffusive_abm.rst index 8863ce6bad..5e265f4109 100644 --- a/docs/source/cpp/diffusive_abm.rst +++ b/docs/source/cpp/diffusive_abm.rst @@ -98,11 +98,11 @@ The model has to be initialized with a vector of agents. Agents have two attribu std::vector agents(100); //Random variables for initialization of agents' position and infection state - auto& pos_rng = mio::UniformDistribution::get_instance(); + auto& pos_rng = mio::UniformDistribution::get_instance(); auto& sta_rng = mio::DiscreteDistribution::get_instance(); //Infection state distribution - std::vector pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.}; + std::vector pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.}; for (auto& a : agents) { //Agents are uniformly distributed in [-2,2]x[-2,2] @@ -114,8 +114,8 @@ Choosing an interaction radius of 0.5 and a noise term of 0.4, the model is init .. code-block:: cpp - double interaction_radius = 0.5; - double noise = 0.4; + ScalarType interaction_radius = 0.5; + ScalarType noise = 0.4; Model model(agents, adoption_rates, interaction_radius, noise); @@ -138,9 +138,9 @@ To simulate the model from `t0` to `tmax` with given step size `dt`, a ``Simulat .. code-block:: cpp - double t0 = 0.0; - double dt = 0.1; - double tmax = 30.; + ScalarType t0 = 0.0; + ScalarType dt = 0.1; + ScalarType tmax = 30.; //Pass the model, t0 and dt to the Simulation auto sim = mio::dabm::Simulation(model, t0, dt); diff --git a/docs/source/cpp/temporal_hybrid.rst b/docs/source/cpp/temporal_hybrid.rst index abf03ada6c..3020400b9c 100644 --- a/docs/source/cpp/temporal_hybrid.rst +++ b/docs/source/cpp/temporal_hybrid.rst @@ -32,7 +32,7 @@ The following example shows how to run a temporal-hybrid simulation of the diffu .. code-block:: cpp using ABM = mio::dabm::Model>; - using ODE = mio::osecir::Model; + using ODE = mio::osecir::Model; First, both models have to be initialized. Please see the documentation of the corresponding models for information about their parameters and how to initialize them. In the following example code block, we initialize a ODE-SECIR model with one age group and a diffusive ABM. We assume that the parameters of the diffusive ABM were created and initialized before and the parameters of the ODE-SECIR model are set such that they match the ABM parameters. @@ -50,19 +50,19 @@ After initializing the models, the corresponding simulations have to be created .. code-block:: cpp //Set t0 and internal dt for each model - double t0 = 0; - double dt = 0.1; + ScalarType t0 = 0; + ScalarType dt = 0.1; //Create simulations auto sim_abm = mio::dabm::Simulation(abm, t0, dt); auto sim_ode = mio::Simulation(ode, t0, dt); const auto result_fct_abm = [](const mio::dabm::Simulation>& sim, - double /*t*/) { + ScalarType /*t*/) { return sim.get_result(); }; - const auto result_fct_ode = [](const mio::Simulation& sim, double /*t*/) { + const auto result_fct_ode = [](const mio::Simulation& sim, ScalarType /*t*/) { return sim.get_result(); }; @@ -71,9 +71,9 @@ Initializing the temporal-hybrid simulation with a given step size `dt_switch` w .. code-block:: cpp //Create hybrid simulation - double dt_switch = 0.2; - mio::hybrid::TemporalHybridSimulation, - mio::TimeSeries> + ScalarType dt_switch = 0.2; + mio::hybrid::TemporalHybridSimulation, + mio::TimeSeries> hybrid_sim(sim_abm, sim_ode, result_fct_abm, result_fct_ode, true, t0, dt_switch); Before advancing the simulation until `tmax`, a switching condition has to be defined. In the example below, the temporal-hybrid model should switch from ABM to ODE if the number of infected individuals is bigger than 20 and it should switch back if the number is below 20. @@ -81,7 +81,7 @@ Before advancing the simulation until `tmax`, a switching condition has to be de .. code-block:: cpp //Define switching condition - const auto condition = [](const mio::TimeSeries& result_abm, const mio::TimeSeries& result_ode, + const auto condition = [](const mio::TimeSeries& result_abm, const mio::TimeSeries& result_ode, bool abm_used) { if (abm_used) { auto& last_value = result_abm.get_last_value().eval(); @@ -109,7 +109,7 @@ Before advancing the simulation until `tmax`, a switching condition has to be de }; //Simulate for 30 days - double tmax = 30.; + ScalarType tmax = 30.; hybrid_sim.advance(tmax, condition); The result ``TimeSeries`` objects of the two models used (which are returned by the above defined result functions) can be accessed and printed via From 30f044ec2ffd2621d014f1a1c578577979cf8b32 Mon Sep 17 00:00:00 2001 From: Anna Wendler <106674756+annawendler@users.noreply.github.com> Date: Mon, 11 May 2026 16:33:49 +0200 Subject: [PATCH 12/12] update docu regarding memilio-simulation wheels --- docs/source/python/python_packages.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/python/python_packages.rst b/docs/source/python/python_packages.rst index 1d05d88615..5471d08442 100644 --- a/docs/source/python/python_packages.rst +++ b/docs/source/python/python_packages.rst @@ -117,7 +117,7 @@ The Python packages can be installed in two ways depending on your use case. **Option 1: Install from PyPI (Recommended - no C++ compiler required; currently only supported for memilio-simulation)** -Pre-built wheels are provided for Linux and Windows on Python 3.8 to 3.13. +Pre-built wheels are provided for Linux and Windows on Python 3.9 to 3.13. .. code-block:: console