468 Abm parameter study

cpp/examples/CMakeLists.txt

@@ -108,6 +108,12 @@ add_executable(abm_minimal_example abm_minimal.cpp)
 target_link_libraries(abm_minimal_example PRIVATE memilio abm)
 target_compile_options(abm_minimal_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+if(MEMILIO_HAS_HDF5)
+    add_executable(abm_parameter_study_example abm_parameter_study.cpp)
+    target_link_libraries(abm_parameter_study_example PRIVATE memilio abm)
+    target_compile_options(abm_parameter_study_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+endif()
+
 add_executable(abm_history_example abm_history_object.cpp)
 target_link_libraries(abm_history_example PRIVATE memilio abm)
 target_compile_options(abm_history_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/examples/abm_minimal.cpp

@@ -21,7 +21,6 @@
 #include "abm/lockdown_rules.h"
 #include "abm/model.h"
 #include "abm/common_abm_loggers.h"
-#include "memilio/utils/abstract_parameter_distribution.h"
 
 #include <fstream>
 

cpp/examples/abm_parameter_study.cpp

@@ -0,0 +1,228 @@
+/*
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Rene Schmieding, Sascha Korf
+*
+* Contact: Martin J. Kuehn <[email protected]>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+#include "abm/result_simulation.h"
+#include "abm/household.h"
+#include "abm/lockdown_rules.h"
+#include "abm/model.h"
+#include "abm/time.h"
+
+#include "memilio/compartments/parameter_studies.h"
+#include "memilio/data/analyze_result.h"
+#include "memilio/io/io.h"
+#include "memilio/io/result_io.h"
+#include "memilio/utils/base_dir.h"
+#include "memilio/utils/logging.h"
+#include "memilio/utils/miompi.h"
+#include "memilio/utils/random_number_generator.h"
+#include "memilio/utils/stl_util.h"
+
+#include <string>
+
+constexpr size_t num_age_groups = 4;
+
+/// An ABM setup taken from abm_minimal.cpp.
+mio::abm::Model make_model(const mio::RandomNumberGenerator& rng)
+{
+
+    const auto age_group_0_to_4   = mio::AgeGroup(0);
+    const auto age_group_5_to_14  = mio::AgeGroup(1);
+    const auto age_group_15_to_34 = mio::AgeGroup(2);
+    const auto age_group_35_to_59 = mio::AgeGroup(3);
+    // Create the model with 4 age groups.
+    auto model      = mio::abm::Model(num_age_groups);
+    model.get_rng() = rng;
+
+    // Set same infection parameter for all age groups. For example, the incubation period is log normally distributed with parameters 4 and 1.
+    model.parameters.get<mio::abm::TimeExposedToNoSymptoms>() = mio::ParameterDistributionLogNormal(4., 1.);
+
+    // Set the age group the can go to school is AgeGroup(1) (i.e. 5-14)
+    model.parameters.get<mio::abm::AgeGroupGotoSchool>()                    = false;
+    model.parameters.get<mio::abm::AgeGroupGotoSchool>()[age_group_5_to_14] = true;
+    // Set the age group the can go to work is AgeGroup(2) and AgeGroup(3) (i.e. 15-34 and 35-59)
+    model.parameters.get<mio::abm::AgeGroupGotoWork>().set_multiple({age_group_15_to_34, age_group_35_to_59}, true);
+
+    // Check if the parameters satisfy their contraints.
+    model.parameters.check_constraints();
+
+    // There are 10 households for each household group.
+    int n_households = 10;
+
+    // For more than 1 family households we need families. These are parents and children and randoms (which are distributed like the data we have for these households).
+    auto child = mio::abm::HouseholdMember(num_age_groups); // A child is 50/50% 0-4 or 5-14.
+    child.set_age_weight(age_group_0_to_4, 1);
+    child.set_age_weight(age_group_5_to_14, 1);
+
+    auto parent = mio::abm::HouseholdMember(num_age_groups); // A parent is 50/50% 15-34 or 35-59.
+    parent.set_age_weight(age_group_15_to_34, 1);
+    parent.set_age_weight(age_group_35_to_59, 1);
+
+    // Two-person household with one parent and one child.
+    auto twoPersonHousehold_group = mio::abm::HouseholdGroup();
+    auto twoPersonHousehold_full  = mio::abm::Household();
+    twoPersonHousehold_full.add_members(child, 1);
+    twoPersonHousehold_full.add_members(parent, 1);
+    twoPersonHousehold_group.add_households(twoPersonHousehold_full, n_households);
+    add_household_group_to_model(model, twoPersonHousehold_group);
+
+    // Three-person household with two parent and one child.
+    auto threePersonHousehold_group = mio::abm::HouseholdGroup();
+    auto threePersonHousehold_full  = mio::abm::Household();
+    threePersonHousehold_full.add_members(child, 1);
+    threePersonHousehold_full.add_members(parent, 2);
+    threePersonHousehold_group.add_households(threePersonHousehold_full, n_households);
+    add_household_group_to_model(model, threePersonHousehold_group);
+
+    // Add one social event with 5 maximum contacts.
+    // Maximum contacs limit the number of people that a person can infect while being at this location.
+    auto event = model.add_location(mio::abm::LocationType::SocialEvent);
+    model.get_location(event).get_infection_parameters().set<mio::abm::MaximumContacts>(5);
+    // Add hospital and ICU with 5 maximum contacs.
+    auto hospital = model.add_location(mio::abm::LocationType::Hospital);
+    model.get_location(hospital).get_infection_parameters().set<mio::abm::MaximumContacts>(5);
+    auto icu = model.add_location(mio::abm::LocationType::ICU);
+    model.get_location(icu).get_infection_parameters().set<mio::abm::MaximumContacts>(5);
+    // Add one supermarket, maximum constacts are assumed to be 20.
+    auto shop = model.add_location(mio::abm::LocationType::BasicsShop);
+    model.get_location(shop).get_infection_parameters().set<mio::abm::MaximumContacts>(20);
+    // At every school, the maximum contacts are 20.
+    auto school = model.add_location(mio::abm::LocationType::School);
+    model.get_location(school).get_infection_parameters().set<mio::abm::MaximumContacts>(20);
+    // At every workplace, maximum contacts are 20.
+    auto work = model.add_location(mio::abm::LocationType::Work);
+    model.get_location(work).get_infection_parameters().set<mio::abm::MaximumContacts>(20);
+
+    // Increase aerosol transmission for all locations
+    model.parameters.get<mio::abm::AerosolTransmissionRates>() = 10.0;
+    // Increase contact rate for all people between 15 and 34 (i.e. people meet more often in the same location)
+    model.get_location(work)
+        .get_infection_parameters()
+        .get<mio::abm::ContactRates>()[{age_group_15_to_34, age_group_15_to_34}] = 10.0;
+
+    // People can get tested at work (and do this with 0.5 probability) from time point 0 to day 10.
+    auto validity_period       = mio::abm::days(1);
+    auto probability           = 0.5;
+    auto start_date            = mio::abm::TimePoint(0);
+    auto end_date              = mio::abm::TimePoint(0) + mio::abm::days(10);
+    auto test_type             = mio::abm::TestType::Antigen;
+    auto test_parameters       = model.parameters.get<mio::abm::TestData>()[test_type];
+    auto testing_criteria_work = mio::abm::TestingCriteria();
+    auto testing_scheme_work   = mio::abm::TestingScheme(testing_criteria_work, validity_period, start_date, end_date,
+                                                         test_parameters, probability);
+    model.get_testing_strategy().add_scheme(mio::abm::LocationType::Work, testing_scheme_work);
+
+    // Assign infection state to each person.
+    // The infection states are chosen randomly with the following distribution
+    std::vector<ScalarType> 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<size_t>::get_instance()(mio::thread_local_rng(), infection_distribution));
+        auto person_rng = mio::abm::PersonalRandomNumberGenerator(person);
+        if (infection_state != mio::abm::InfectionState::Susceptible) {
+            person.add_new_infection(mio::abm::Infection(person_rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
+                                                         model.parameters, start_date, infection_state));
+        }
+    }
+
+    // Assign locations to the people
+    for (auto& person : model.get_persons()) {
+        const auto id = person.get_id();
+        //assign shop and event
+        model.assign_location(id, event);
+        model.assign_location(id, shop);
+        //assign hospital and ICU
+        model.assign_location(id, hospital);
+        model.assign_location(id, icu);
+        //assign work/school to people depending on their age
+        if (person.get_age() == age_group_5_to_14) {
+            model.assign_location(id, school);
+        }
+        if (person.get_age() == age_group_15_to_34 || person.get_age() == age_group_35_to_59) {
+            model.assign_location(id, work);
+        }
+    }
+
+    // During the lockdown, social events are closed for 90% of people.
+    auto t_lockdown = mio::abm::TimePoint(0) + mio::abm::days(10);
+    mio::abm::close_social_events(t_lockdown, 0.9, model.parameters);
+
+    return model;
+}
+
+int main()
+{
+    mio::mpi::init();
+
+    mio::set_log_level(mio::LogLevel::warn);
+
+    // Set start and end time for the simulation.
+    auto t0   = mio::abm::TimePoint(0);
+    auto tmax = t0 + mio::abm::days(5);
+    // auto sim  = mio::abm::Simulation(t0, std::move(model));
+    const size_t num_runs = 3;
+
+    // Create a parameter study. The ABM currently does not use parameters or dt, so we set them both to 0.
+    mio::ParameterStudy study(0, t0, tmax, mio::abm::TimeSpan(0), num_runs);
+
+    // Optional: set seeds to get reproducable results
+    // study.get_rng().seed({12341234, 53456, 63451, 5232576, 84586, 52345});
+
+    const std::string result_dir = mio::path_join(mio::base_dir(), "example_results");
+    if (!mio::create_directory(result_dir)) {
+        mio::log_error("Could not create result directory \"{}\".", result_dir);
+        return 1;
+    }
+
+    auto ensemble_results = study.run(
+        [](auto, auto t0_, auto, size_t) {
+            return mio::abm::ResultSimulation(make_model(mio::thread_local_rng()), t0_);
+        },
+        [result_dir](auto&& sim, auto&& run_idx) {
+            auto interpolated_result = mio::interpolate_simulation_result(sim.get_result());
+            std::string outpath = mio::path_join(result_dir, "abm_minimal_run_" + std::to_string(run_idx) + ".txt");
+            std::ofstream outfile_run(outpath);
+            sim.get_result().print_table(outfile_run, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4);
+            std::cout << "Results written to " << outpath << std::endl;
+            auto params = std::vector<mio::abm::Model>{};
+            return std::vector{interpolated_result};
+        });
+
+    if (ensemble_results.size() > 0) {
+        auto ensemble_results_p05 = ensemble_percentile(ensemble_results, 0.05);
+        auto ensemble_results_p25 = ensemble_percentile(ensemble_results, 0.25);
+        auto ensemble_results_p50 = ensemble_percentile(ensemble_results, 0.50);
+        auto ensemble_results_p75 = ensemble_percentile(ensemble_results, 0.75);
+        auto ensemble_results_p95 = ensemble_percentile(ensemble_results, 0.95);
+
+        mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups,
+                                mio::path_join(result_dir, "Results_" + std::string("p05") + ".h5")));
+        mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups,
+                                mio::path_join(result_dir, "Results_" + std::string("p25") + ".h5")));
+        mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups,
+                                mio::path_join(result_dir, "Results_" + std::string("p50") + ".h5")));
+        mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups,
+                                mio::path_join(result_dir, "Results_" + std::string("p75") + ".h5")));
+        mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups,
+                                mio::path_join(result_dir, "Results_" + std::string("p95") + ".h5")));
+    }
+
+    mio::mpi::finalize();
+
+    return 0;
+}

cpp/examples/ode_secir_parameter_study.cpp

@@ -17,10 +17,14 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
+#include "memilio/config.h"
+#include "memilio/utils/base_dir.h"
+#include "memilio/utils/miompi.h"
+#include "memilio/utils/stl_util.h"
+#include "ode_secir/model.h"
 #include "ode_secir/parameters_io.h"
 #include "ode_secir/parameter_space.h"
 #include "memilio/compartments/parameter_studies.h"
-#include "memilio/mobility/metapopulation_mobility_instant.h"
 #include "memilio/io/result_io.h"
 
 /**
@@ -30,13 +34,10 @@
  * @param t0 starting point of simulation
  * @param tmax end point of simulation
  */
-mio::IOResult<void>
-write_single_run_result(const size_t run,
-                        const mio::Graph<mio::SimulationNode<ScalarType, mio::osecir::Simulation<ScalarType>>,
-                                         mio::MobilityEdge<ScalarType>>& graph)
+mio::IOResult<void> write_single_run_result(const size_t run, const mio::osecir::Simulation<ScalarType>& sim)
 {
-    std::string abs_path;
-    BOOST_OUTCOME_TRY(auto&& created, mio::create_directory("results", abs_path));
+    std::string abs_path = mio::path_join(mio::base_dir(), "example_results");
+    BOOST_OUTCOME_TRY(auto&& created, mio::create_directory(abs_path));
 
     if (run == 0) {
         std::cout << "Results are stored in " << abs_path << '\n';
@@ -46,44 +47,29 @@ write_single_run_result(const size_t run,
     }
 
     //write sampled parameters for this run
-    //omit edges to save space as they are not sampled
-    int inode = 0;
-    for (auto&& node : graph.nodes()) {
-        BOOST_OUTCOME_TRY(auto&& js_node_model, serialize_json(node.property.get_result(), mio::IOF_OmitDistributions));
-        Json::Value js_node(Json::objectValue);
-        js_node["NodeId"]  = node.id;
-        js_node["Model"]   = js_node_model;
-        auto node_filename = mio::path_join(abs_path, "Parameters_run" + std::to_string(run) + "_node" +
-                                                          std::to_string(inode++) + ".json");
-        BOOST_OUTCOME_TRY(mio::write_json(node_filename, js_node));
-    }
+    auto node_filename = mio::path_join(abs_path, "Parameters_run" + std::to_string(run) + ".json");
+    BOOST_OUTCOME_TRY(mio::write_json(node_filename, sim.get_result()));
 
     //write results for this run
     std::vector<mio::TimeSeries<ScalarType>> all_results;
     std::vector<int> ids;
 
-    ids.reserve(graph.nodes().size());
-    all_results.reserve(graph.nodes().size());
-    std::transform(graph.nodes().begin(), graph.nodes().end(), std::back_inserter(all_results), [](auto& node) {
-        return node.property.get_result();
-    });
-    std::transform(graph.nodes().begin(), graph.nodes().end(), std::back_inserter(ids), [](auto& node) {
-        return node.id;
-    });
-    auto num_groups = (int)(size_t)graph.nodes()[0].property.get_simulation().get_model().parameters.get_num_groups();
-    BOOST_OUTCOME_TRY(mio::save_result(all_results, ids, num_groups,
-                                       mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5"))));
+    BOOST_OUTCOME_TRY(mio::save_result({sim.get_result()}, {0}, (int)sim.get_model().parameters.get_num_groups().get(),
+                                       mio::path_join(abs_path, "Results_run" + std::to_string(run) + ".h5")));
 
     return mio::success();
 }
 
 int main()
 {
-    mio::set_log_level(mio::LogLevel::debug);
+    mio::mpi::init();
+    mio::set_log_level(mio::LogLevel::warn);
 
     ScalarType t0   = 0;
     ScalarType tmax = 50;
+    ScalarType dt   = 0.1;
 
+    // set up model with parameters
     ScalarType cont_freq = 10; // see Polymod study
 
     ScalarType num_total_t0 = 10000, num_exp_t0 = 100, num_inf_t0 = 50, num_car_t0 = 50, num_hosp_t0 = 20,
@@ -139,22 +125,31 @@ int main()
         return -1;
     }
 
-    //create study
-    auto num_runs = size_t(1);
-    mio::ParameterStudy<ScalarType, mio::osecir::Simulation<ScalarType>> parameter_study(model, t0, tmax, num_runs);
+    // create study
+    auto num_runs = size_t(3);
+    mio::ParameterStudy parameter_study(model, t0, tmax, dt, num_runs);
 
-    //run study
-    auto sample_graph = [](auto&& graph) {
-        return mio::osecir::draw_sample(graph);
+    // set up for run
+    auto sample_graph = [](const auto& model_, ScalarType t0_, ScalarType dt_, size_t) {
+        mio::osecir::Model<ScalarType> copy = model_;
+        mio::osecir::draw_sample(copy);
+        return mio::osecir::Simulation<ScalarType>(std::move(copy), t0_, dt_);
     };
-    auto handle_result = [](auto&& graph, auto&& run) {
-        auto write_result_status = write_single_run_result(run, graph);
+    auto handle_result = [](auto&& sim, auto&& run) {
+        auto write_result_status = write_single_run_result(run, sim);
         if (!write_result_status) {
             std::cout << "Error writing result: " << write_result_status.error().formatted_message();
         }
-        return 0; //Result handler must return something, but only meaningful when using MPI.
+        return 0; // Result handler must return something.
     };
-    parameter_study.run(sample_graph, handle_result);
+
+    // Optional: set seeds to get reproducable results
+    // parameter_study.get_rng().seed({1456, 157456, 521346, 35345, 6875, 6435});
+
+    // run study
+    auto result = parameter_study.run(sample_graph, handle_result);
+
+    mio::mpi::finalize();
 
     return 0;
 }