Merge branch 'main' of github.com:SciCompMod/memilio-tutorials

Author: annawendler

.github/workflows/run-tutorials.yml

@@ -36,7 +36,7 @@ jobs:
           # Install memilio-simulation from the cloned repo (C++ extension, needs cmake/ninja)
           uv add --dev ./memilio/pycode/memilio-simulation
           # Install tutorial dependencies
-          uv add marimo matplotlib pandas numpy pyabc[pyarrow] bayesflow tensorflow nbformat
+          uv add -r requirements.txt -c requirements.txt
 
       - name: Run tutorials and export to HTML
         run: |

cpp-tutorials/abm/parameter_setter.cpp

@@ -24,8 +24,8 @@
  *
  * MEmilio provides an agent-based model (ABM) where each individual ("agent")
  * has an age, a home, and assigned locations they visit during the day.
- * Agents transition between infection states (Susceptible, Exposed, ..., Dead)
- * based on stochastic contact events at their current location.
+ * Agents' transitions between infection states (Susceptible, Exposed, ..., Dead)
+ * are based on stochastic contact events at their current location.
  *
  * This tutorial demonstrates how to:
  *  1. Define age groups and configure model parameters.
@@ -37,9 +37,9 @@
  *
  * Key concept -- HouseholdMember age weights:
  *   Each HouseholdMember carries an integer weight for every AgeGroup.
- *   When a person is created for that member slot, the model draws its age
+ *   When a person is created as a HouseholdMember, the model draws its age
  *   from a discrete distribution proportional to these weights.
- *   Example: weights {1, 1, 0, 0} give a 50/50 chance of AgeGroup 0 or 1.
+ *   Example: weights {1, 1, 0, 0} gives a 50/50 chance of AgeGroup 0 or 1.
  */
 
 #include "abm/household.h"
@@ -56,7 +56,7 @@ int main(int argc, char* argv[])
 {
     // Usage: tutorial_abm_household [n_households] [infected_frac] [sim_days]
     //   n_households  : number of each household type           (default: 125)
-    //   infected_frac : fraction initially infected             (default: 0.2)
+    //   infected_frac : fraction of initially infected             (default: 0.2)
     //   sim_days      : simulation duration in days             (default: 30)
     int arg_n_households     = (argc > 1) ? std::atoi(argv[1]) : 125;
     double arg_infected_frac = (argc > 2) ? std::atof(argv[2]) : 0.2;
@@ -77,14 +77,10 @@ int main(int argc, char* argv[])
     const auto age_group_60_to_79 = mio::AgeGroup(4); // seniors
     const auto age_group_80_plus  = mio::AgeGroup(5); // elderly
 
-    // *** Create the model and set infection parameters. ***
-    // The Model holds all persons, locations, and parameters. We pass in the
-    // number of age groups so that all parameter arrays are sized correctly.
-    // `set_local_parameters` and `set_world_parameters` fill in realistic
-    // epidemiological values (see parameter_setter.h).
+    // *** Create the model. ***
+    // The Model holds all persons, locations, and parameters. We need to add the locations
+    // before we set the parameters as some parameters depend on the locations.
     auto model = mio::abm::Model(num_age_groups);
-    set_local_parameters(model);
-    set_world_parameters(model.parameters);
 
     // Define which age groups are eligible to go to school and to work.
     // The AgeGroupGotoSchool / AgeGroupGotoWork arrays default to false for
@@ -121,7 +117,7 @@ int main(int argc, char* argv[])
 
     // *** Compose households and add them to the model. ***
     // A Household collects (member_type, count) pairs. A HouseholdGroup bundles
-    // many copies of a Household template. `add_household_group_to_model`
+    // many copies of a Household template. `add_household_group_to_model`,
     // creates the actual persons and their home locations.
     //
     // CLI parameters (see usage at top of main):
@@ -187,10 +183,16 @@ int main(int argc, char* argv[])
     // One workplace for all working adults.
     auto work = model.add_location(mio::abm::LocationType::Work);
 
+    // *** Set paramters for all locations. ***
+    // `set_local_parameters` and `set_world_parameters` fill in realistic
+    // epidemiological values (see parameter_setter.h).
+    set_local_parameters(model);
+    set_world_parameters(model.parameters);
+
     // *** Assign initial infection states. ***
     // Each person draws a random infection state from the distribution below.
     // Persons who are not Susceptible receive a full Infection object so their
-    // viral-load course and state transitions are properly initialised.
+    // viral-load course and state transitions are properly initialized.
     //
     //  Index | InfectionState                 | Probability
     //  ------|--------------------------------|------------
@@ -202,7 +204,7 @@ int main(int argc, char* argv[])
     //    5   | InfectedCritical   (I_Crit)    | 0.0
     //    6   | Recovered                      | 0.0
     //    7   | Dead                           | 0.0
-    auto start_date = mio::abm::TimePoint(0); // t = 0 s from the simulation epoch
+    auto start_date = mio::abm::TimePoint(0); // t = 0
 
     // Build infection distribution from the infected fraction.
     // The non-susceptible portion is split: 25% Exposed, 50% I_NS, 25% I_Sy.
@@ -264,4 +266,4 @@ int main(int argc, char* argv[])
     std::cout << "Results written to abm_household.txt\n";
 
     return 0;
-}
+}

cpp-tutorials/abm/tutorial_abm_households.cpp

@@ -64,7 +64,7 @@ const auto age_group_80_plus  = mio::AgeGroup(5);
 //   - Is active for the full 30-day simulation window.
 //   - Uses a PCR test (sensitivity/specificity come from TestData defaults).
 //   - Tests are valid for `validity_days`; a negative result exempts retesting for validity_days days (Default: 3).
-//   - Targets persons in any infected state (Exposed through Critical).
+//   - Targets persons in any infected state (Exposed to Critical).
 //   - Applies to all age groups except school children (5-14).
 //   - Tests with `testing_probability` probability upon entry (Default: 100%).
 //   - Applies to all public locations (SocialEvent, School, Work, BasicsShop)
@@ -90,7 +90,7 @@ void add_npi_testing_strategies_to_world(mio::abm::Model& model, double testing_
 
     auto testing_criteria = mio::abm::TestingCriteria(ages_to_test, states_to_test);
     auto testing_scheme   = mio::abm::TestingScheme(testing_criteria, validity, start_date_test, end_date_test,
-                                                  pcr_test_parameters, testing_probability);
+                                                    pcr_test_parameters, testing_probability);
 
     // Attach the scheme to all public location types.
     model.get_testing_strategy().add_scheme(mio::abm::LocationType::SocialEvent, testing_scheme);
@@ -117,8 +117,6 @@ int main(int argc, char* argv[])
 
     // *** Create the model and set infection parameters. ***
     auto model = mio::abm::Model(num_age_groups);
-    set_local_parameters(model);
-    set_world_parameters(model.parameters);
 
     // Define which age groups are eligible to go to school and to work.
     // The AgeGroupGotoSchool / AgeGroupGotoWork arrays default to false for
@@ -218,6 +216,10 @@ int main(int argc, char* argv[])
     // One workplace for all working adults.
     auto work = model.add_location(mio::abm::LocationType::Work);
 
+    // *** Set paramters for all locations (same as Tutorial 1). ***
+    set_local_parameters(model);
+    set_world_parameters(model.parameters);
+
     // *** Assign initial infection states. ***
     //
     //  Index | InfectionState          | Probability
@@ -291,4 +293,4 @@ int main(int argc, char* argv[])
     std::cout << "Results written to abm_tests.txt\n";
 
     return 0;
-}
+}

cpp-tutorials/abm/tutorial_abm_testing.cpp

@@ -186,7 +186,7 @@ int main(int argc, char* argv[])
     //   use_data_vacc     : 1 = apply JSON data-driven vaccination, 0 = skip (default: 0)
     //                         Note that no vaccination takes place if json file should be used but cannot be found!
     double arg_vacc_rate       = (argc > 1) ? std::atof(argv[1]) : 0.0;
-    int arg_n_households       = (argc > 2) ? std::atoi(argv[2]) : 1000;
+    int arg_n_households       = (argc > 2) ? std::atoi(argv[2]) : 125;
     double arg_protection_peak = (argc > 3) ? std::atof(argv[3]) : 0.67;
     int arg_use_data_vacc      = (argc > 4) ? std::atoi(argv[4]) : 0;
 
@@ -202,10 +202,8 @@ int main(int argc, char* argv[])
     const auto age_group_60_to_79 = mio::AgeGroup(4); // seniors
     const auto age_group_80_plus  = mio::AgeGroup(5); // elderly
 
-    // *** Create the model and set infection parameters (same as Tutorial 1). ***
+    // *** Create the model (same as Tutorial 1). ***
     auto model = mio::abm::Model(num_age_groups);
-    set_local_parameters(model);
-    set_world_parameters(model.parameters);
 
     // Define which age groups are eligible to go to school and to work.
     // The AgeGroupGotoSchool / AgeGroupGotoWork arrays default to false for
@@ -238,7 +236,7 @@ int main(int argc, char* argv[])
         model.parameters.get<mio::abm::SeverityProtectionFactor>()[{mio::abm::ProtectionType::GenericVaccine, age,
                                                                     mio::abm::VirusVariant::Wildtype}] =
             mio::TimeSeriesFunctor<ScalarType>{mio::TimeSeriesFunctorType::LinearInterpolation,
-                                               {{0, 0.0}, {1, 0.85}, {180, 0.70}}};
+                                               {{0, 0.0}, {14, 0.85}, {180, 0.70}}};
     }
 
     std::cout << "Vaccination protection factors configured.\n";
@@ -321,6 +319,10 @@ int main(int argc, char* argv[])
     // One workplace for all working adults.
     auto work = model.add_location(mio::abm::LocationType::Work);
 
+    // *** Set paramters for all locations (same as Tutorial 1). ***
+    set_local_parameters(model);
+    set_world_parameters(model.parameters);
+
     // *** Assign initial infection states. ***
     //
     //  Index | InfectionState          | Probability
@@ -419,10 +421,9 @@ int main(int argc, char* argv[])
         // ── Data-driven mode ─────────────────────────────────────────
         // Read real-world vaccination counts from a JSON file and apply
         // them day-by-day to the model's persons.
-        const std::string vacc_json_path =
-            "/Users/saschakorf/Documents/Promotion/memilio-tutorials/cpp-tutorials/abm/vacc_county_ageinf_ma7.json";
-        const int county_id       = 1002;
-        const mio::Date sim_start = mio::Date(2020, 10, 1);
+        const std::string vacc_json_path = "./vacc_county_ageinf_ma7.json";
+        const int county_id              = 1002;
+        const mio::Date sim_start        = mio::Date(2020, 10, 1);
 
         auto vacc_map = prepare_vaccination_data(vacc_json_path, num_age_groups, county_id);
 
@@ -452,4 +453,4 @@ int main(int argc, char* argv[])
     std::cout << "Results written to abm_vaccination.txt\n";
 
     return 0;
-}
+}

cpp-tutorials/abm/tutorial_abm_vaccination.cpp

@@ -48,12 +48,12 @@ int main()
     // all currently active dampings of one level and different types are summed up. If two dampings have different
     // levels (independent of the type) they are combined multiplicatively. In the following we apply a `Damping`
     // of 0.9 after 10 days and another damping of 0.6 after 20 days which means that the contacts are reduced
-    // by 10% and 40%, respectively. To always retain a minimum level of contacts, a minimum contact frequency can
+    // by 90% and 60%, respectively. To always retain a minimum level of contacts, a minimum contact frequency can
     // be set that is never deceeded. In our example we set this minimum contact rate to 0.
     contact_matrix[0].add_damping(0.9, mio::SimulationTime<ScalarType>(10.));
     contact_matrix[0].add_damping(0.6, mio::SimulationTime<ScalarType>(20.));
 
-    //EXERCISE: Please create a damping that replaces the 10% reduction after 10 days by a 20% reduction. Additionally, add a damping after 40 days that increases the contact rate by 50%.
+    //EXERCISE: Please create a damping that replaces the 90% reduction after 10 days by a 20% reduction. Additionally, add a damping after 40 days that increases the contact rate by 50%.
 
     // Again, we start with 0.5% of the population initially in `Exposed` and 0.5% initially in `InfectedNoSymptoms`
     // while the remaining 99% is `Susceptible`.

cpp-tutorials/exercises/exercise03.cpp

@@ -1 +1 @@
-set(GIT_TAG_MEMILIO "70648381af0d37ff0f6ab6be8995bcae1a110fc6")
+set(GIT_TAG_MEMILIO "20cee90e7ce55b0107f0607df36e50b9e4a49889")

cpp-tutorials/git_tag.cmake

@@ -46,7 +46,7 @@ int main()
     // all currently active dampings of one level and different types are summed up. If two dampings have different
     // levels (independent of the type) they are combined multiplicatively. In the following we apply a `Damping`
     // of 0.9 after 10 days and another damping of 0.6 after 20 days which means that the contacts are reduced
-    // by 10% and 40%, respectively. To always retain a minimum level of contacts, a minimum contact frequency can
+    // by 90% and 60%, respectively. To always retain a minimum level of contacts, a minimum contact frequency can
     // be set that is never deceeded. In our example we set this minimum contact rate to 0.
     contact_matrix[0].add_damping(0.9, mio::SimulationTime<ScalarType>(10.));
     contact_matrix[0].add_damping(0.6, mio::SimulationTime<ScalarType>(20.));

cpp-tutorials/tutorial03.cpp

@@ -250,7 +250,7 @@ def _(mo):
 
 @app.cell
 def _(result):
-    result.print_table()
+    print(result.print_table(return_string=True))
     return
 
 
@@ -280,7 +280,7 @@ def _(mo):
 
 @app.cell
 def _(interpolated_result):
-    interpolated_result.print_table()
+    print(interpolated_result.print_table(return_string=True))
     return
 
 

exercises/exercise01.py

@@ -47,7 +47,7 @@ def _(mo):
 def _(np):
     import memilio.simulation.osecir as osecir
     from memilio.simulation import AgeGroup, LogLevel, set_log_level
-    set_log_level(LogLevel.Off)
+    set_log_level(LogLevel.Error)
 
     # Initialize total population, simulation start time, timeframe, and step size
     total_population = 100000
@@ -76,11 +76,14 @@ def _(np):
     model.parameters.DeathsPerCritical[group] = 0.3
 
     # Set contact frequency
-    model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones((1, 1)) * 10
+    model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones(
+        (1, 1)) * 10
 
     # Initialize compartments
-    model.populations[group, osecir.InfectionState.Exposed] = 0.005 * total_population
-    model.populations[group, osecir.InfectionState.InfectedNoSymptoms] = 0.005 * total_population
+    model.populations[group,
+                      osecir.InfectionState.Exposed] = 0.005 * total_population
+    model.populations[group,
+                      osecir.InfectionState.InfectedNoSymptoms] = 0.005 * total_population
     model.populations.set_difference_from_total(
         (group, osecir.InfectionState.Susceptible), total_population)
     return dt, model, osecir, t0, tmax
@@ -122,7 +125,8 @@ def _(dt, model, osecir, t0, tmax):
 
 @app.cell
 def _(results):
-    compartments = results[0]  # The first element of results contains the compartment data
+    # The first element of results contains the compartment data
+    compartments = results[0]
     flows = results[1]  # The second element of results contains the flow data
     return compartments, flows
 
@@ -162,7 +166,6 @@ def _(flows, np, osecir):
     plot_time = time_days[1:]
 
     print(flow_array.shape)
-    print("ji")
     return daily_flows, plot_time
 
 
@@ -196,12 +199,13 @@ def _(mo):
 def _(compartments, daily_flows, osecir, plot_time, plt):
     # 1. Prepare compartment data for comparison
     # Interpolate to match the time grid of our flows
-    comp_array = osecir.interpolate_simulation_result(compartments).as_ndarray()
+    comp_array = osecir.interpolate_simulation_result(
+        compartments).as_ndarray()
 
     # Extract the "InfectedSevere" compartment (Currently Hospitalized)
     # Adding 1 because row 0 is the time axis
     state_severe_idx = 1 + int(osecir.InfectionState.InfectedSevere)
-    current_severe = comp_array[state_severe_idx, 1:] 
+    current_severe = comp_array[state_severe_idx, 1:]
 
     # 2. Flow indices
     new_symptomatic_idx = 2
@@ -216,8 +220,10 @@ def _(compartments, daily_flows, osecir, plot_time, plt):
     fig, ax = plt.subplots(1, 2, figsize=(14, 6))
 
     # --- Subplot 1: Disease Progression (The Time Lag) ---
-    ax[0].bar(plot_time, daily_flows[new_symptomatic_idx, :], color='coral', alpha=0.5, label='New Symptomatic Cases')
-    ax[0].bar(plot_time, daily_flows[new_hospitalized_idx, :], color='darkred', alpha=0.8, label='New Hospitalizations')
+    ax[0].bar(plot_time, daily_flows[new_symptomatic_idx, :],
+              color='coral', alpha=0.5, label='New Symptomatic Cases')
+    ax[0].bar(plot_time, daily_flows[new_hospitalized_idx, :],
+              color='darkred', alpha=0.8, label='New Hospitalizations')
     # EXERCISE: Extend subplot 1 to also show new ICU admissions.
     # Hint: add a third bar for daily_flows[new_icu_idx, :] with a suitable color and label,
     # e.g. color='maroon', alpha=0.9, label='New ICU Admissions'.
@@ -228,14 +234,16 @@ def _(compartments, daily_flows, osecir, plot_time, plt):
 
     # --- Subplot 2: Incidence vs. Bed Occupancy ---
     # Axis 1: Daily New Cases (Incidence)
-    ax[1].bar(plot_time, daily_flows[new_hospitalized_idx, :], color='steelblue', alpha=0.5, label='Daily Admissions (Incidence)')
+    ax[1].bar(plot_time, daily_flows[new_hospitalized_idx, :],
+              color='steelblue', alpha=0.5, label='Daily Admissions (Incidence)')
     ax[1].set_xlabel('Time [days]')
     ax[1].set_ylabel('Daily New Hospitalizations [#]', color='steelblue')
     ax[1].tick_params(axis='y', labelcolor='steelblue')
 
     # Axis 2: Current Occupancy
     ax2 = ax[1].twinx()
-    ax2.plot(plot_time, current_severe, color='black', linewidth=3, label='Currently Hospitalized')
+    ax2.plot(plot_time, current_severe, color='black',
+             linewidth=3, label='Currently Hospitalized')
     ax2.set_ylabel('Total Hospital Bed Occupancy [#]', color='black')
     ax2.tick_params(axis='y', labelcolor='black')