diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 39a12539e1..fa12445425 100644 --- a/cpp/examples/CMakeLists.txt +++ b/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}) diff --git a/cpp/examples/abm_minimal.cpp b/cpp/examples/abm_minimal.cpp index 156b212b4f..643a17b1df 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/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 diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp new file mode 100644 index 0000000000..4d609f323c --- /dev/null +++ b/cpp/examples/abm_parameter_study.cpp @@ -0,0 +1,228 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Rene Schmieding, Sascha Korf +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#include "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 + +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::ParameterDistributionLogNormal(4., 1.); + + // Set the age group the can go to school is AgeGroup(1) (i.e. 5-14) + model.parameters.get() = false; + model.parameters.get()[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().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(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(5); + auto icu = model.add_location(mio::abm::LocationType::ICU); + model.get_location(icu).get_infection_parameters().set(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(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(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(20); + + // Increase aerosol transmission for all locations + model.parameters.get() = 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()[{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()[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 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)); + 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{}; + 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; +} diff --git a/cpp/examples/ode_secir_parameter_study.cpp b/cpp/examples/ode_secir_parameter_study.cpp index 8ad8c1cd2b..7450226ef4 100644 --- a/cpp/examples/ode_secir_parameter_study.cpp +++ b/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 -write_single_run_result(const size_t run, - const mio::Graph>, - mio::MobilityEdge>& graph) +mio::IOResult write_single_run_result(const size_t run, const mio::osecir::Simulation& 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> all_results; std::vector 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> 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 copy = model_; + mio::osecir::draw_sample(copy); + return mio::osecir::Simulation(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; } diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index b2e8260cc3..f17d019f71 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -20,7 +20,6 @@ #include "memilio/compartments/parameter_studies.h" #include "memilio/config.h" -#include "memilio/geography/regions.h" #include "memilio/io/epi_data.h" #include "memilio/io/result_io.h" #include "memilio/io/mobility_io.h" @@ -29,13 +28,11 @@ #include "memilio/utils/miompi.h" #include "memilio/utils/random_number_generator.h" #include "memilio/utils/time_series.h" +#include "ode_secir/model.h" #include "ode_secir/parameters_io.h" #include "ode_secir/parameter_space.h" -#include "boost/filesystem.hpp" -#include "memilio/utils/stl_util.h" #include #include -#include /** * Set a value and distribution of an UncertainValue. @@ -265,9 +262,7 @@ int main() 2, 1, Eigen::VectorX::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.2), indices_save_edges); - //run parameter study - auto parameter_study = mio::ParameterStudy>{ - params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)}; + mio::ParameterStudy parameter_study(params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)); if (mio::mpi::is_root()) { printf("Seeds: "); @@ -279,10 +274,13 @@ int main() auto save_single_run_result = mio::IOResult(mio::success()); auto ensemble = parameter_study.run( - [](auto&& graph) { - return draw_sample(graph); + [](auto&& graph, ScalarType t0, ScalarType dt, size_t) { + auto copy = graph; + return mio::make_sampled_graph_simulation>( + mio::osecir::draw_sample(copy), t0, dt, dt); }, - [&](auto results_graph, auto&& run_id) { + [&](auto&& results_sim, auto&& run_id) { + auto results_graph = results_sim.get_graph(); auto interpolated_result = mio::interpolate_simulation_result(results_graph); auto params = std::vector>{}; @@ -319,7 +317,7 @@ int main() boost::filesystem::path results_dir("test_results"); bool created = boost::filesystem::create_directories(results_dir); if (created) { - mio::log_info("Directory '{:s}' was created.", results_dir.string()); + mio::log_info("Directory '{}' was created.", results_dir.string()); } auto county_ids = std::vector{1001, 1002, 1003}; diff --git a/cpp/examples/ode_secir_read_graph.cpp b/cpp/examples/ode_secir_read_graph.cpp index 8832aec976..41cdf89b24 100644 --- a/cpp/examples/ode_secir_read_graph.cpp +++ b/cpp/examples/ode_secir_read_graph.cpp @@ -17,43 +17,36 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "memilio/mobility/metapopulation_mobility_instant.h" -#include "memilio/io/mobility_io.h" #include "memilio/compartments/parameter_studies.h" +#include "memilio/config.h" +#include "memilio/io/cli.h" +#include "memilio/io/mobility_io.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/base_dir.h" + +#include "memilio/utils/stl_util.h" +#include "ode_secir/model.h" #include "ode_secir/parameter_space.h" #include "ode_secir/parameters_io.h" -#include -#include -std::string setup(int argc, char** argv, const std::string data_dir) -{ - if (argc == 2) { - std::cout << "Using file " << argv[1] << " in data/Germany/mobility." << std::endl; - return mio::path_join(data_dir, "Germany", "mobility", (std::string)argv[1]); - } - else { - if (argc > 2) { - mio::log_error("Too many arguments given."); - } - else { - mio::log_warning("No arguments given."); - } - auto mobility_file = "commuter_mobility_2022.txt"; - std::cout << "Using file " << mobility_file << " in data/Germany/mobility." << std::endl; - std::cout << "Usage: read_graph MOBILITY_FILE" - << "\n\n"; - std::cout - << "This example performs a simulation based on mobility data from the German Federal Employment Agency." - << std::endl; - return mio::path_join(data_dir, "Germany", "mobility", mobility_file); - } -} +#include int main(int argc, char** argv) { mio::set_log_level(mio::LogLevel::critical); - std::string data_dir = DATA_DIR; - std::string filename = setup(argc, argv, data_dir); + + auto parameters = + mio::cli::ParameterSetBuilder() + .add<"MobilityFile">( + mio::path_join(mio::base_dir(), "data", "Germany", "mobility", "commuter_mobility_2022.txt"), + {.description = "Create the mobility file with MEmilio Epidata's getCommuterMobility.py file."}) + .build(); + + auto result = mio::command_line_interface(argv[0], argc, argv, parameters, {"MobilityFile"}); + if (!result) { + std::cout << result.error().message(); + return result.error().code().value(); + } const auto t0 = 0.; const auto tmax = 10.; @@ -110,10 +103,10 @@ int main(int argc, char** argv) mio::osecir::set_params_distributions_normal(model, t0, tmax, 0.2); std::cout << "Reading Mobility File..." << std::flush; - auto read_mobility_result = mio::read_mobility_plain(filename); + auto read_mobility_result = mio::read_mobility_plain(parameters.get<"MobilityFile">()); if (!read_mobility_result) { - std::cout << read_mobility_result.error().formatted_message() << '\n'; - std::cout << "Create the mobility file with MEmilio Epidata's getCommuterMobility.py file." << '\n'; + std::cout << "\n" << read_mobility_result.error().formatted_message() << "\n"; + std::cout << "Create the mobility file with MEmilio Epidata's getCommuterMobility.py file.\n"; return 0; } auto& commuter_mobility = read_mobility_result.value(); @@ -137,7 +130,8 @@ int main(int argc, char** argv) std::cout << "Writing Json Files..." << std::flush; auto write_status = mio::write_graph(graph, "graph_parameters"); if (!write_status) { - std::cout << "Error: " << write_status.error().formatted_message(); + std::cout << "\n" << write_status.error().formatted_message(); + return 0; } std::cout << "Done" << std::endl; @@ -145,13 +139,19 @@ int main(int argc, char** argv) auto graph_read_result = mio::read_graph>("graph_parameters"); if (!graph_read_result) { - std::cout << "Error: " << graph_read_result.error().formatted_message(); + std::cout << "\n" << graph_read_result.error().formatted_message(); + return 0; } std::cout << "Done" << std::endl; auto& graph_read = graph_read_result.value(); std::cout << "Running Simulations..." << std::flush; - auto study = mio::ParameterStudy>(graph_read, t0, tmax, 0.5, 2); + mio::ParameterStudy study(graph_read, t0, tmax, 0.5, 2); + study.run_serial([](auto&& g, auto t0_, auto dt_, auto) { + auto copy = g; + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); + }); std::cout << "Done" << std::endl; return 0; diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 6de772c360..359c30bedf 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -21,118 +21,133 @@ #define MIO_COMPARTMENTS_PARAMETER_STUDIES_H #include "memilio/io/binary_serializer.h" +#include "memilio/io/io.h" #include "memilio/mobility/graph_simulation.h" #include "memilio/utils/logging.h" +#include "memilio/utils/metaprogramming.h" #include "memilio/utils/miompi.h" #include "memilio/utils/random_number_generator.h" -#include "memilio/utils/time_series.h" #include "memilio/mobility/metapopulation_mobility_instant.h" -#include "memilio/compartments/simulation.h" +#include #include +#include #include -#include -#include -#include +#include +#include +#include namespace mio { /** - * Class that performs multiple simulation runs with randomly sampled parameters. - * Can simulate mobility graphs with one simulation in each node or single simulations. - * @tparam S type of simulation that runs in one node of the graph. + * @brief Class used to performs multiple simulation runs with randomly sampled parameters. + * @tparam ParameterType The parameters used to create simulations. + * @tparam TimeType The type used for time, e.g. double or TimePoint. + * @tparam StepType The type used for time steps, e.g. double or TimeStep. */ -template +template class ParameterStudy { public: - /** - * The type of simulation of a single node of the graph. - */ - using Simulation = S; - /** - * The Graph type that stores the parametes of the simulation. - * This is the input of ParameterStudies. - */ - using ParametersGraph = Graph>; - /** - * The Graph type that stores simulations and their results of each run. - * This is the output of ParameterStudies for each run. - */ - using SimulationGraph = Graph, MobilityEdge>; + using Parameters = ParameterType; + using Time = TimeType; + using Step = StepType; + +private: + template + requires std::is_invocable_v + using SimulationT = std::decay_t>; + + template + requires std::is_invocable_v, size_t> + using ProcessedResultT = std::decay_t< + std::invoke_result_t, size_t>>; + +public: + // TODO: replacement for "set_params_distributions_normal"? Maybe a special ctor for UncertainParameterSet? /** - * create study for graph of compartment models. - * @param graph graph of parameters - * @param t0 start time of simulations - * @param tmax end time of simulations - * @param graph_sim_dt time step of graph simulation - * @param num_runs number of runs + * @brief Create a parameter study with some parameters. + * The simulation type is determined when calling any "run" member function. + * @param parameters The parameters used to create simulations. + * @param t0 Start time of simulations. + * @param tmax End time of simulations. + * @param dt Initial time step of simulations. + * @param num_runs Number of simulations that will be created and run. */ - ParameterStudy(const ParametersGraph& graph, FP t0, FP tmax, FP graph_sim_dt, size_t num_runs) - : m_graph(graph) + ParameterStudy(const Parameters& parameters, Time t0, Time tmax, Step dt, size_t num_runs) + : m_parameters(parameters) , m_num_runs(num_runs) , m_t0{t0} , m_tmax{tmax} - , m_dt_graph_sim(graph_sim_dt) + , m_dt(dt) { } /** - * create study for graph of compartment models. - * Creates distributions for all parameters of the models in the graph. - * @param graph graph of parameters - * @param t0 start time of simulations - * @param tmax end time of simulations - * @param dev_rel relative deviation of the created distributions from the initial value. - * @param graph_sim_dt time step of graph simulation - * @param num_runs number of runs + * @brief Run all simulations in serial. + * @param[in] create_simulation A callable sampling the study's parameters and return a simulation. + * @param[in] process_simulation_result (Optional) A callable that takes the simulation and processes its result. + * @return A vector that contains (processed) simulation results for each run. + * + * Important side effect: Calling this function overwrites seed and counter of thread_local_rng(). + * Use this RNG when sampling parameters in create_simulation. + * + * The function signature for create_simulation is + * `SimulationT(const Parameters& study_parameters, Time t0, Step dt, size_t run_idx)`, + * where SimulationT is some kind of simulation. + * The function signature for process_simulation_result is + * `ProcessedResultT(SimulationT&&, size_t run_index)`, + * where ProcessedResultT is a (de)serializable result. + * @{ */ - ParameterStudy(const ParametersGraph& graph, FP t0, FP tmax, FP dev_rel, FP graph_sim_dt, size_t num_runs) - : ParameterStudy(graph, t0, tmax, graph_sim_dt, num_runs) + template + std::vector> + run_serial(CreateSimulationFunction&& create_simulation, + ProcessSimulationResultFunction&& process_simulation_result) { - for (auto& params_node : m_graph.nodes()) { - set_params_distributions_normal(params_node, t0, tmax, dev_rel); - } + return run_impl(0, m_num_runs, std::forward(create_simulation), + std::forward(process_simulation_result)); } - /** - * @brief Create study for single compartment model. - * @param model compartment model with initial values - * @param t0 start time of simulations - * @param tmax end time of simulations - * @param num_runs number of runs in ensemble run - */ - ParameterStudy(typename Simulation::Model const& model, FP t0, FP tmax, size_t num_runs) - : ParameterStudy({}, t0, tmax, tmax - t0, num_runs) + template + std::vector> run_serial(CreateSimulationFunction&& create_simulation) { - m_graph.add_node(0, model); + return run_serial( + std::forward(create_simulation), + [](SimulationT&& sim, size_t) -> SimulationT&& { + return std::move(sim); + }); } + /** @} */ /** - * @brief Carry out all simulations in the parameter study. - * Save memory and enable more runs by immediately processing and/or discarding the result. - * The result processing function is called when a run is finished. It receives the result of the run - * (a SimulationGraph object) and an ordered index. The values returned by the result processing function - * are gathered and returned as a list. - * This function is parallelized if memilio is configured with MEMILIO_ENABLE_MPI. - * The MPI processes each contribute a share of the runs. The sample function and result processing function - * are called in the same process that performs the run. The results returned by the result processing function are - * gathered at the root process and returned as a list by the root in the same order as if the programm - * were running sequentially. Processes other than the root return an empty list. - * @param sample_graph Function that receives the ParametersGraph and returns a sampled copy. - * @param result_processing_function Processing function for simulation results, e.g., output function. - * @returns At the root process, a list of values per run that have been returned from the result processing function. - * At all other processes, an empty list. - * @tparam SampleGraphFunction Callable type, accepts instance of ParametersGraph. - * @tparam HandleSimulationResultFunction Callable type, accepts instance of SimulationGraph and an index of type size_t. + * @brief Run all simulations distributed over multiple MPI ranks. + * @param[in] create_simulation A callable sampling the study's parameters and return a simulation. + * @param[in] process_simulation_result A callable that takes the simulation and processes its result. + * @return A vector that contains processed simulation results for each run. + * + * Important: Do not forget to use mio::mpi::init and finalize when using this function! + * + * Important side effect: Calling this function overwrites seed and counter of thread_local_rng(). + * Use this RNG when sampling parameters in create_simulation. + * + * The function signature for create_simulation is + * `SimulationT(const Parameters& study_parameters, Time t0, Step dt, size_t run_idx)`, + * where SimulationT is some kind of simulation. + * The function signature for process_simulation_result is + * `ProcessedResultT(SimulationT&&, size_t run_index)`, + * where ProcessedResultT is a (de)serializable result. */ - template - std::vector> - run(SampleGraphFunction sample_graph, HandleSimulationResultFunction result_processing_function) + template + std::vector> + run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) { + using EnsembleResultT = + std::vector>; int num_procs, rank; + #ifdef MEMILIO_ENABLE_MPI MPI_Comm_size(mpi::get_world(), &num_procs); MPI_Comm_rank(mpi::get_world(), &rank); @@ -145,40 +160,15 @@ class ParameterStudy //So we set our own RNG to be used. //Assume that sampling uses the thread_local_rng() and isn't multithreaded m_rng.synchronize(); - thread_local_rng() = m_rng; - auto run_distribution = distribute_runs(m_num_runs, num_procs); - auto start_run_idx = + std::vector run_distribution = distribute_runs(m_num_runs, num_procs); + size_t start_run_idx = std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0)); - auto end_run_idx = start_run_idx + run_distribution[size_t(rank)]; - - std::vector> ensemble_result; - ensemble_result.reserve(m_num_runs); - - for (size_t run_idx = start_run_idx; run_idx < end_run_idx; run_idx++) { - log(LogLevel::info, "ParameterStudies: run {}", run_idx); - - //prepare rng for this run by setting the counter to the right offset - //Add the old counter so that this call of run() produces different results - //from the previous call - auto run_rng_counter = m_rng.get_counter() + rng_totalsequence_counter( - static_cast(run_idx), Counter(0)); - thread_local_rng().set_counter(run_rng_counter); - - //sample - auto sim = create_sampled_simulation(sample_graph); - log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", - (thread_local_rng().get_counter() - run_rng_counter).get()); - - //perform run - sim.advance(m_tmax); + size_t end_run_idx = start_run_idx + run_distribution[size_t(rank)]; - //handle result and store - ensemble_result.emplace_back(result_processing_function(std::move(sim).get_graph(), run_idx)); - } - - //Set the counter of our RNG so that future calls of run() produce different parameters. - m_rng.set_counter(m_rng.get_counter() + rng_totalsequence_counter(m_num_runs, Counter(0))); + EnsembleResultT ensemble_result = + run_impl(start_run_idx, end_run_idx, std::forward(create_simulation), + std::forward(process_simulation_result)); #ifdef MEMILIO_ENABLE_MPI //gather results @@ -189,7 +179,7 @@ class ParameterStudy ByteStream bytes(bytes_size); MPI_Recv(bytes.data(), bytes.data_size(), MPI_BYTE, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE); - auto src_ensemble_results = deserialize_binary(bytes, Tag{}); + IOResult src_ensemble_results = deserialize_binary(bytes, Tag{}); if (!src_ensemble_results) { log_error("Error receiving ensemble results from rank {}.", src_rank); } @@ -198,8 +188,8 @@ class ParameterStudy } } else { - auto bytes = serialize_binary(ensemble_result); - auto bytes_size = int(bytes.data_size()); + ByteStream bytes = serialize_binary(ensemble_result); + int bytes_size = int(bytes.data_size()); MPI_Send(&bytes_size, 1, MPI_INT, 0, 0, mpi::get_world()); MPI_Send(bytes.data(), bytes.data_size(), MPI_BYTE, 0, 0, mpi::get_world()); ensemble_result.clear(); //only return root process @@ -209,154 +199,118 @@ class ParameterStudy return ensemble_result; } - /** - * @brief Carry out all simulations in the parameter study. - * Convenience function for a few number of runs, but can use more memory because it stores all runs until the end. - * Unlike the other overload, this function is not MPI-parallel. - * @return vector of SimulationGraph for each run. - */ - template - std::vector run(SampleGraphFunction sample_graph) - { - std::vector ensemble_result; - ensemble_result.reserve(m_num_runs); - - //The ParameterDistributions used for sampling parameters use thread_local_rng() - //So we set our own RNG to be used. - //Assume that sampling uses the thread_local_rng() and isn't multithreaded - thread_local_rng() = m_rng; - - for (size_t i = 0; i < m_num_runs; i++) { - log(LogLevel::info, "ParameterStudies: run {}", i); - - //prepare rng for this run by setting the counter to the right offset - //Add the old counter so that this call of run() produces different results - //from the previous call - auto run_rng_counter = m_rng.get_counter() + - rng_totalsequence_counter(static_cast(i), Counter(0)); - thread_local_rng().set_counter(run_rng_counter); - - auto sim = create_sampled_simulation(sample_graph); - log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", - (thread_local_rng().get_counter() - run_rng_counter).get()); - - sim.advance(m_tmax); - - ensemble_result.emplace_back(std::move(sim).get_graph()); - } - - //Set the counter of our RNG so that future calls of run() produce different parameters. - m_rng.set_counter(m_rng.get_counter() + rng_totalsequence_counter(m_num_runs, Counter(0))); - - return ensemble_result; - } - - /** - * @brief sets the number of Monte Carlo runs - * @param[in] num_runs number of runs - */ - void set_num_runs(size_t num_runs) - { - m_num_runs = num_runs; - } - - /** - * @brief returns the number of Monte Carlo runs - */ - int get_num_runs() const - { - return static_cast(m_num_runs); - } - - /** - * @brief sets end point in simulation - * @param[in] tmax end point in simulation - */ - void set_tmax(FP tmax) + /// @brief Return the number of total runs that the study will make. + size_t get_num_runs() const { - m_tmax = tmax; + return m_num_runs; } - /** - * @brief returns end point in simulation - */ - FP get_tmax() const + /// @brief Return the final time point for simulations. + Time get_tmax() const { return m_tmax; } - void set_t0(FP t0) - { - m_t0 = t0; - } - - /** - * @brief returns start point in simulation - */ - FP get_t0() const + /// @brief Return the initial time point for simulations. + Time get_t0() const { return m_t0; } - /** - * Get the input model that the parameter study is run for. - * Use for single node simulations, use get_model_graph for graph simulations. - * @{ - */ - const typename Simulation::Model& get_model() const + /// @brief Return the initial step sized used by simulations. + Time get_dt() const { - return m_graph.nodes()[0].property; + return m_dt; } - typename Simulation::Model& get_model() - { - return m_graph.nodes()[0].property; - } - /** @} */ /** - * Get the input graph that the parameter study is run for. - * Use for graph simulations, use get_model for single node simulations. + * @brief Get the input parameters that each simulation in the study is created from. * @{ */ - const ParametersGraph& get_model_graph() const + const Parameters& get_parameters() const { - return m_graph; + return m_parameters; } - ParametersGraph& get_model_graph() + Parameters& get_parameters() { - return m_graph; + return m_parameters; } /** @} */ + /// @brief Access the study's random number generator. RandomNumberGenerator& get_rng() { return m_rng; } private: - //sample parameters and create simulation - template - GraphSimulation create_sampled_simulation(SampleGraphFunction sample_graph) + /** + * @brief Main loop creating and running simulations. + * @param[in] start_run_idx, end_run_idx Range of indices. Performs one run for each index. + * @param[in] create_simulation A callable sampling the study's parameters and return a simulation. + * @param[in] process_simulation_result A callable that takes the simulation and processes its result. + * @return A vector that contains processed simulation results for each run. + * + * Important side effect: Calling this function overwrites seed and counter of thread_local_rng(). + * Use this RNG when sampling parameters in create_simulation. + */ + template + std::vector> + run_impl(size_t start_run_idx, size_t end_run_idx, CreateSimulationFunction&& create_simulation, + ProcessSimulationResultFunction&& process_simulation_result) { - SimulationGraph sim_graph; + assert(start_run_idx <= end_run_idx); - auto sampled_graph = sample_graph(m_graph); - for (auto&& node : sampled_graph.nodes()) { - sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration); - } - for (auto&& edge : sampled_graph.edges()) { - sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + // Note that this overwrites seed and counter of thread_local_rng, but it does not replace it. + thread_local_rng() = m_rng; + + std::vector> ensemble_result; + ensemble_result.reserve(m_num_runs); + + for (size_t run_idx = start_run_idx; run_idx < end_run_idx; run_idx++) { + log(LogLevel::info, "ParameterStudies: run {}", run_idx); + + //prepare rng for this run by setting the counter to the right offset + //Add the old counter so that this call of run() produces different results + //from the previous call + Counter run_rng_counter = + m_rng.get_counter() + + rng_totalsequence_counter(static_cast(run_idx), Counter(0)); + thread_local_rng().set_counter(run_rng_counter); + + //sample + SimulationT sim = + create_simulation(std::as_const(m_parameters), std::as_const(m_t0), std::as_const(m_dt), run_idx); + log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", + (thread_local_rng().get_counter() - run_rng_counter).get()); + + //perform run + sim.advance(m_tmax); + + //handle result and store + ensemble_result.emplace_back(process_simulation_result(std::move(sim), run_idx)); } - return make_mobility_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); + //Set the counter of our RNG so that future calls of run() produce different parameters. + m_rng.set_counter(m_rng.get_counter() + rng_totalsequence_counter(m_num_runs, Counter(0))); + + return ensemble_result; } - std::vector distribute_runs(size_t num_runs, int num_procs) + /** + * @brief Distribute a number of runs over a number of processes. + * Processes with low ranks get additional runs, if the number is not evenly divisible. + * @param num_runs The total number of runs. + * @param num_procs The total number of processes, i.e. the size of MPI_Comm. + * @return A vector of size num_procs with the number of runs each process should make. + */ + static std::vector distribute_runs(size_t num_runs, int num_procs) { + assert(num_procs > 0); //evenly distribute runs //lower processes do one more run if runs are not evenly distributable - auto num_runs_local = num_runs / num_procs; //integer division! - auto remainder = num_runs % num_procs; + size_t num_runs_local = num_runs / num_procs; //integer division! + size_t remainder = num_runs % num_procs; std::vector run_distribution(num_procs); std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); @@ -365,24 +319,32 @@ class ParameterStudy return run_distribution; } -private: - // Stores Graph with the names and ranges of all parameters - ParametersGraph m_graph; - - size_t m_num_runs; - - // Start time (should be the same for all simulations) - FP m_t0; - // End time (should be the same for all simulations) - FP m_tmax; - // time step of the graph - FP m_dt_graph_sim; - // adaptive time step of the integrator (will be corrected if too large/small) - FP m_dt_integration = 0.1; - // - RandomNumberGenerator m_rng; + ParameterType m_parameters; ///< Stores parameters used to create a simulation for each run. + size_t m_num_runs; ///< Total number of runs (i.e. simulations) to do when calling "run". + Time m_t0, m_tmax; ///< Start and end time for the simulations. + Step m_dt; ///< Initial step size of the simulation. Some integrators may adapt their step size during simulation. + RandomNumberGenerator m_rng; ///< The random number generator used by the study. }; +//sample parameters and create simulation +template +auto make_sampled_graph_simulation(const Graph>& sampled_graph, FP t0, + FP dt_node_sim, FP dt_graph_sim) +{ + using SimGraph = Graph, MobilityEdge>; + + SimGraph sim_graph; + + for (auto&& node : sampled_graph.nodes()) { + sim_graph.add_node(node.id, node.property, t0, dt_node_sim); + } + for (auto&& edge : sampled_graph.edges()) { + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + } + + return make_mobility_sim(t0, dt_graph_sim, std::move(sim_graph)); +} + } // namespace mio #endif // MIO_COMPARTMENTS_PARAMETER_STUDIES_H diff --git a/cpp/memilio/io/cli.h b/cpp/memilio/io/cli.h index 1f0fefcf69..395f8b1ce7 100644 --- a/cpp/memilio/io/cli.h +++ b/cpp/memilio/io/cli.h @@ -293,6 +293,9 @@ class AbstractParameter : public DatalessParameter else { // deserialize failed // insert more information to the error message std::string msg = "While setting \"" + name + "\": " + param_result.error().message(); + if (param_result.error().message() == "Json value is not a string.") { + msg += " Try using escaped quotes (\\\") around your input strings."; + } return mio::failure(param_result.error().code(), msg); } }) diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 33770f8a47..0ceb33a821 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -31,10 +31,12 @@ namespace mio /** * @brief abstract simulation on a graph with alternating node and edge actions */ -template +template class GraphSimulationBase { public: + using Graph = GraphT; + using node_function = node_f; using edge_function = edge_f; diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 7046cb63ac..6f8a0b4c0c 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -46,6 +46,8 @@ template class SimulationNode { public: + using Simulation = Sim; + template ::value, void>> SimulationNode(Args&&... args) : m_simulation(std::forward(args)...) diff --git a/cpp/memilio/utils/base_dir.h b/cpp/memilio/utils/base_dir.h index 0b396ceb27..d390dcd6ce 100644 --- a/cpp/memilio/utils/base_dir.h +++ b/cpp/memilio/utils/base_dir.h @@ -22,13 +22,15 @@ #include "memilio/config.h" +#include + namespace mio { /** - * @brief Returns path to the repo directory. -*/ -constexpr std::string mio::base_dir() + * @brief Returns the absolute path to the project directory. + */ +const static std::string base_dir() { return MEMILIO_BASE_DIR; } diff --git a/cpp/memilio/utils/index.h b/cpp/memilio/utils/index.h index 595db59f37..ac76fe3b38 100644 --- a/cpp/memilio/utils/index.h +++ b/cpp/memilio/utils/index.h @@ -21,7 +21,11 @@ #define INDEX_H #include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/metaprogramming.h" #include "memilio/utils/type_safe.h" +#include +#include +#include namespace mio { @@ -29,6 +33,41 @@ namespace mio template class Index; +namespace details +{ + +template +std::tuple...> get_tuple(Index i) +{ + if constexpr (sizeof...(Tags) == 1) { + return std::tuple(i); + } + else { + return i.indices; + } +} + +template +std::tuple> get_tuple(Enum i) + requires std::is_enum::value +{ + return std::tuple(Index(i)); +} + +// template +// Index merge_index_from_tuple(std::tuple...>); + +template +auto merge_indices(IndexArgs... args) +{ + return std::tuple_cat(get_tuple(args)...); +} + +template +using merge_indices_expr = decltype(std::tuple(get_tuple(std::declval())...)); + +} // namespace details + /** * @brief An Index with a single template parameter is a typesafe wrapper for size_t * that is associated with a Tag. It is used to index into a CustomIndexArray @@ -139,6 +178,13 @@ class Index } public: + template ::value>> + Index(IndexArgs... _index_args) + : indices(details::merge_indices(_index_args...)) + { + } + // comparison operators bool operator==(Index const& other) const { @@ -181,12 +227,20 @@ class Index std::tuple...> indices; }; -template -struct type_at_index> : public type_at_index { +/// Specialization of type_at_index for Index. @see type_at_index. +template +struct type_at_index> : public type_at_index { +}; + +/// Specialization of index_of_type for Index. @see index_of_type. +template +struct index_of_type> : public index_of_type { }; -template -struct index_of_type> : public index_of_type { +/// Specialization of index_of_type for Index. Resolves ambiguity when using Index%s as items. @see index_of_type. +template +struct index_of_type, Index> { + static constexpr std::size_t value = 0; }; // retrieves the Index at the Ith position for a Index with more than one Tag diff --git a/cpp/memilio/utils/miompi.cpp b/cpp/memilio/utils/miompi.cpp index afbfcff117..8c1079d538 100644 --- a/cpp/memilio/utils/miompi.cpp +++ b/cpp/memilio/utils/miompi.cpp @@ -18,6 +18,7 @@ * limitations under the License. */ #include "memilio/utils/miompi.h" +#include "memilio/utils/logging.h" #ifdef MEMILIO_ENABLE_MPI #include @@ -41,6 +42,8 @@ void init() { #ifdef MEMILIO_ENABLE_MPI MPI_Init(nullptr, nullptr); +#else + mio::log_debug("Using mio::mpi::init without MPI being enabled."); #endif } diff --git a/cpp/memilio/utils/type_list.h b/cpp/memilio/utils/type_list.h index 6ff50a720a..43d02aea12 100644 --- a/cpp/memilio/utils/type_list.h +++ b/cpp/memilio/utils/type_list.h @@ -61,6 +61,12 @@ template struct index_of_type> : public index_of_type { }; +/// Specialization of index_of_type for TypeList. Resolves ambiguity when using TypeLists as items. @see index_of_type. +template +struct index_of_type, TypeList> { + static constexpr std::size_t value = 0; +}; + } // namespace mio #endif // MIO_UTILS_TYPE_LIST_H_ diff --git a/cpp/models/abm/CMakeLists.txt b/cpp/models/abm/CMakeLists.txt index 50ce4258f6..ae04b68b2d 100644 --- a/cpp/models/abm/CMakeLists.txt +++ b/cpp/models/abm/CMakeLists.txt @@ -4,6 +4,7 @@ add_library(abm location_id.h household.cpp household.h + result_simulation.h simulation.cpp simulation.h person.cpp diff --git a/cpp/models/abm/model.h b/cpp/models/abm/model.h index 0c6b320e25..acbda1fa3d 100644 --- a/cpp/models/abm/model.h +++ b/cpp/models/abm/model.h @@ -20,6 +20,7 @@ #ifndef MIO_ABM_MODEL_H #define MIO_ABM_MODEL_H +#include "abm/infection_state.h" #include "abm/model_functions.h" #include "abm/location_type.h" #include "abm/mobility_data.h" @@ -61,6 +62,7 @@ class Model using MobilityRuleType = LocationType (*)(PersonalRandomNumberGenerator&, const Person&, TimePoint, TimeSpan, const Parameters&); + using Compartments = mio::abm::InfectionState; /** * @brief Create a Model. * @param[in] num_agegroups The number of AgeGroup%s in the simulated Model. Must be less than MAX_NUM_AGE_GROUPS. diff --git a/cpp/models/abm/result_simulation.h b/cpp/models/abm/result_simulation.h new file mode 100644 index 0000000000..5a51a6e99e --- /dev/null +++ b/cpp/models/abm/result_simulation.h @@ -0,0 +1,64 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Rene Schmieding +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "abm/common_abm_loggers.h" +#include "abm/simulation.h" + +namespace mio +{ +namespace abm +{ + +/// @brief Simulation holding its own History to provide a get_result member. Can be used for a ParameterStudy. +template +class ResultSimulation : public Simulation +{ +public: + using Model = M; + + /// @brief Create a simulation, moving the model. + ResultSimulation(Model&& m, TimePoint t) + : Simulation(t, std::move(m)) + { + } + + /** + * @brief Run the simulation until the given time point. + * @param tmax Final time point for the simualtion. + */ + void advance(TimePoint tmax) + { + Simulation::advance(tmax, history); + } + + /** + * @brief Return the simulation result aggregated by infection states. + */ + const mio::TimeSeries& get_result() const + { + return get<0>(history.get_log()); + } + + mio::History history{ + Eigen::Index(InfectionState::Count)}; ///< History used to create the result TimeSeries. +}; + +} // namespace abm +} // namespace mio diff --git a/cpp/models/abm/simulation.h b/cpp/models/abm/simulation.h index 2f973cf092..2778e7cea8 100644 --- a/cpp/models/abm/simulation.h +++ b/cpp/models/abm/simulation.h @@ -35,14 +35,15 @@ namespace abm template class Simulation { - public: + using Model = M; + /** * @brief Create a simulation. * @param[in] t0 The starting time of the Simulation. * @param[in] model The Model to simulate. */ - Simulation(TimePoint t0, M&& model) + Simulation(TimePoint t0, Model&& model) : m_model(std::move(model)) , m_t(t0) , m_dt(hours(1)) @@ -87,11 +88,11 @@ class Simulation /** * @brief Get the Model that this Simulation evolves. */ - M& get_model() + Model& get_model() { return m_model; } - const M& get_model() const + const Model& get_model() const { return m_model; } @@ -105,7 +106,7 @@ class Simulation m_t += m_dt; } - M m_model; ///< The Model to simulate. + Model m_model; ///< The Model to simulate. TimePoint m_t; ///< The current TimePoint of the Simulation. TimeSpan m_dt; ///< The length of the time steps. }; diff --git a/cpp/tests/test_abm_simulation.cpp b/cpp/tests/test_abm_simulation.cpp index 53c372fe7c..856456baf8 100644 --- a/cpp/tests/test_abm_simulation.cpp +++ b/cpp/tests/test_abm_simulation.cpp @@ -18,10 +18,14 @@ * limitations under the License. */ #include "abm/location_type.h" +#include "abm/simulation.h" +#include "abm/result_simulation.h" +#include "abm/time.h" #include "abm_helpers.h" #include "abm/common_abm_loggers.h" #include "matchers.h" #include "memilio/io/history.h" +#include #include TEST(TestSimulation, advance_random) @@ -143,3 +147,28 @@ TEST(TestSimulation, advanceWithCommonHistory) 3); // Check if all persons are in the delta-logger Mobility helper entry 0, 3 persons EXPECT_EQ(logMobilityInfoDelta[1].size(), 3); // Check if all persons are in the delta-log first entry, 3 persons } + +TEST(TestSimulation, ResultSimulation) +{ + // run a ResultSimulation on a minimal setup + auto model = mio::abm::Model(num_age_groups); + auto location = model.add_location(mio::abm::LocationType::Home); + auto person = model.add_person(location, age_group_15_to_34); + + model.assign_location(person, location); + + const auto t0 = mio::abm::TimePoint(0) + mio::abm::hours(100); + const auto tmax = t0 + mio::abm::hours(50); + auto sim = mio::abm::ResultSimulation(std::move(model), t0); + + // run simulation. expect one timepopint per day, but nothing to change in the results + sim.advance(tmax); + const size_t N = (size_t)(tmax - t0).hours() + 1; + ASSERT_EQ(sim.get_result().get_num_time_points(), N); + EXPECT_THAT(sim.get_result().get_times(), ElementsAreLinspace(t0.days(), tmax.days(), N)); + for (const auto& tp : sim.get_result()) { + EXPECT_EQ(tp.sum(), 1.0); + } + EXPECT_EQ(sim.get_result().get_value(0)[(Eigen::Index)mio::abm::InfectionState::Susceptible], 1.0); + EXPECT_EQ(sim.get_result().get_value(N - 1)[(Eigen::Index)mio::abm::InfectionState::Susceptible], 1.0); +} diff --git a/cpp/tests/test_io_cli.cpp b/cpp/tests/test_io_cli.cpp index 1b70f67681..0721044ceb 100644 --- a/cpp/tests/test_io_cli.cpp +++ b/cpp/tests/test_io_cli.cpp @@ -301,18 +301,24 @@ TEST_F(TestCLI, test_get_param) TEST_F(TestCLI, test_set_param) { Params p; - auto set = mio::cli::details::AbstractSet::build(p).value(); - auto id = mio::cli::details::Identifier::make_raw("A"); + auto set = mio::cli::details::AbstractSet::build(p).value(); + auto id_A = mio::cli::details::Identifier::make_raw("A"); + auto id_B = mio::cli::details::Identifier::make_raw("B"); + auto wrong_id = mio::cli::details::Identifier::make_raw("NotAnOption"); // use normally - EXPECT_THAT(print_wrap(set.set_param(id, std::string("5.2"))), IsSuccess()); + EXPECT_THAT(print_wrap(set.set_param(id_A, std::string("5.2"))), IsSuccess()); EXPECT_EQ(p.get(), 5.2); + EXPECT_THAT(print_wrap(set.set_param(id_B, std::string("\"5.2\""))), IsSuccess()); + EXPECT_EQ(p.get(), "5.2"); // cause errors - auto result = set.set_param(id, std::string("definitely a double")); + auto result = set.set_param(id_A, std::string("definitely a double")); EXPECT_THAT(print_wrap(result), IsFailure(mio::StatusCode::InvalidType)); - EXPECT_EQ(result.error().formatted_message(), - "Invalid type: While setting \"A\": Json value is not a double.\n" - "* Line 1, Column 1\n Syntax error: value, object or array expected.\n"); - result = set.set_param(mio::cli::details::Identifier::make_raw("NotAnOption"), std::string()); + EXPECT_EQ(result.error().message(), "While setting \"A\": Json value is not a double.\n" + "* Line 1, Column 1\n Syntax error: value, object or array expected.\n"); + result = set.set_param(id_B, std::string("this is string is missing quotes")); + EXPECT_THAT(print_wrap(result), IsFailure(mio::StatusCode::InvalidType)); + EXPECT_THAT(result.error().message(), testing::HasSubstr("\\\"")); // check only for additional hint + result = set.set_param(wrong_id, std::string()); EXPECT_THAT(print_wrap(result), IsFailure(mio::StatusCode::KeyNotFound)); EXPECT_EQ(result.error().formatted_message(), "Key not found: Could not set parameter: No such option \"NotAnOption\"."); diff --git a/cpp/tests/test_parameter_studies.cpp b/cpp/tests/test_parameter_studies.cpp index e0f46fbc5a..01e9f6e6a1 100644 --- a/cpp/tests/test_parameter_studies.cpp +++ b/cpp/tests/test_parameter_studies.cpp @@ -17,13 +17,17 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include "memilio/config.h" +#include "memilio/utils/miompi.h" #include "memilio/utils/parameter_distributions.h" #include "ode_secir/model.h" #include "ode_secir/parameter_space.h" #include "memilio/compartments/parameter_studies.h" #include "memilio/mobility/metapopulation_mobility_instant.h" #include "memilio/utils/random_number_generator.h" +#include #include +#include #include TEST(ParameterStudies, sample_from_secir_params) @@ -157,14 +161,17 @@ TEST(ParameterStudies, sample_graph) graph.add_node(1, model); graph.add_edge(0, 1, mio::MobilityParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 8), 1.0))); - auto study = mio::ParameterStudy>(graph, 0.0, 0.0, 0.5, 1); + mio::ParameterStudy study(graph, 0.0, 0.0, 0.5, 1); mio::log_rng_seeds(study.get_rng(), mio::LogLevel::warn); - auto results = study.run([](auto&& g) { - return draw_sample(g); + auto ensemble_results = study.run_serial([](auto&& g, auto t0_, auto dt_, auto) { + auto copy = g; + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); }); - EXPECT_EQ(results[0].edges()[0].property.get_parameters().get_coefficients()[0].get_dampings().size(), 1); - for (auto& node : results[0].nodes()) { + auto& results = ensemble_results.at(0); + EXPECT_EQ(results.get_graph().edges()[0].property.get_parameters().get_coefficients()[0].get_dampings().size(), 1); + for (auto& node : results.get_graph().nodes()) { auto& result_model = node.property.get_simulation().get_model(); EXPECT_EQ(result_model.parameters.get>() .get_cont_freq_mat()[0] @@ -370,26 +377,24 @@ TEST(ParameterStudies, check_ensemble_run_result) mio::ContactMatrix(Eigen::MatrixXd::Constant((size_t)num_groups, (size_t)num_groups, fact * cont_freq)); mio::osecir::set_params_distributions_normal(model, t0, tmax, 0.2); - mio::ParameterStudy> parameter_study(model, t0, tmax, 1); + mio::ParameterStudy parameter_study(model, t0, tmax, 0.1, 1); mio::log_rng_seeds(parameter_study.get_rng(), mio::LogLevel::warn); // Run parameter study - parameter_study.set_num_runs(1); - auto graph_results = parameter_study.run([](auto&& g) { - return draw_sample(g); + auto ensemble_results = parameter_study.run_serial([](auto&& model_, auto t0_, auto dt_, auto) { + auto copy = model_; + draw_sample(copy); + return mio::osecir::Simulation(copy, t0_, dt_); }); - std::vector> results; - for (size_t i = 0; i < graph_results.size(); i++) { - results.push_back(std::move(graph_results[i].nodes()[0].property.get_result())); - } + const mio::TimeSeries& results = ensemble_results.at(0).get_result(); - for (Eigen::Index i = 0; i < results[0].get_num_time_points(); i++) { + for (Eigen::Index i = 0; i < results.get_num_time_points(); i++) { std::vector total_at_ti((size_t)mio::osecir::InfectionState::Count, 0); - for (Eigen::Index j = 0; j < results[0][i].size(); j++) { // number of compartments per time step - EXPECT_GE(results[0][i][j], 0.0) << " day " << results[0].get_time(i) << " group " << j; - total_at_ti[static_cast(j) / (size_t)mio::osecir::InfectionState::Count] += results[0][i][j]; + for (Eigen::Index j = 0; j < results[i].size(); j++) { // number of compartments per time step + EXPECT_GE(results[i][j], 0.0) << " day " << results.get_time(i) << " group " << j; + total_at_ti[static_cast(j) / (size_t)mio::osecir::InfectionState::Count] += results[i][j]; } for (auto j = mio::AgeGroup(0); j < params.get_num_groups(); j++) { @@ -398,3 +403,64 @@ TEST(ParameterStudies, check_ensemble_run_result) } } } + +namespace +{ + +struct MockStudyParams { + const int init, run; +}; + +struct MockStudySim { + MockStudySim(const MockStudyParams& p_, double t0_, double dt_) + : p(p_) + , t0(t0_) + , dt(dt_) + { + } + void advance(double t) + { + tmax = t; + } + + MockStudyParams p; + double t0, dt; + double tmax = 0; +}; + +} // namespace + +TEST(ParameterStudies, mocked_run) +{ + // run a very simple study, that works with mpi + const double t0 = 20, tmax = 21, dt = 22; + const MockStudyParams params{23, -1}; + const size_t num_runs = 5; // enough to notice MPI effects + const auto make_sim = [&](auto&& params_, auto t0_, auto dt_, auto i_) { + MockStudyParams cp{params_.init, (int)i_}; + return MockStudySim(cp, t0_, dt_); + }; + const auto process_sim = [&](MockStudySim&& s, size_t i) { + return s.tmax + i; + }; + const double process_sim_result = (num_runs * tmax) + num_runs * (num_runs - 1) / 2.; + mio::ParameterStudy study(params, t0, tmax, dt, num_runs); + // case: run_serial without processing; expect created simulations in order + auto result_serial = study.run_serial(make_sim); + EXPECT_EQ(result_serial.size(), num_runs); + for (int i = 0; const auto& sim : result_serial) { + EXPECT_EQ(sim.t0, t0); + EXPECT_EQ(sim.dt, dt); + EXPECT_EQ(sim.tmax, tmax); + EXPECT_EQ(sim.p.init, params.init); + EXPECT_EQ(sim.p.run, i++); + } + // case: run and run_serial with processing; expect the same (unordered) result for both, on all ranks + // Note: currently the tests are not make use of MPI, so we expect the same result from each rank + auto result_serial_processed = study.run_serial(make_sim, process_sim); + auto result_parallel = study.run(make_sim, process_sim); + for (const auto& result : {result_serial_processed, result_parallel}) { + EXPECT_EQ(result.size(), num_runs); + EXPECT_EQ(std::accumulate(result.begin(), result.end(), 0.0), process_sim_result); + } +} diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index a15c62d419..f706327ce6 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -169,8 +169,8 @@ TEST(TestSaveResult, save_result_with_params) graph.add_edge(0, 1, mio::MobilityParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0))); - auto num_runs = 3; - auto parameter_study = mio::ParameterStudy>(graph, 0.0, 2.0, 0.5, num_runs); + auto num_runs = 3; + mio::ParameterStudy parameter_study(graph, 0.0, 2.0, 0.5, num_runs); mio::log_rng_seeds(parameter_study.get_rng(), mio::LogLevel::warn); TempFileRegister tmp_file_register; @@ -183,10 +183,13 @@ TEST(TestSaveResult, save_result_with_params) ensemble_params.reserve(size_t(num_runs)); auto save_result_status = mio::IOResult(mio::success()); parameter_study.run( - [](auto&& g) { - return draw_sample(g); + [](auto&& g, auto t0_, auto dt_, auto) { + auto copy = g; + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); }, - [&](auto results_graph, auto run_idx) { + [&](auto&& results, auto run_idx) { + auto results_graph = results.get_graph(); ensemble_results.push_back(mio::interpolate_simulation_result(results_graph)); ensemble_params.emplace_back(); @@ -302,8 +305,8 @@ TEST(TestSaveResult, save_percentiles_and_sums) mio::MobilityParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0), indices_save_edges)); - auto num_runs = 3; - auto parameter_study = mio::ParameterStudy>(graph, 0.0, 2.0, 0.5, num_runs); + auto num_runs = 3; + mio::ParameterStudy parameter_study(graph, 0.0, 2.0, 0.5, num_runs); mio::log_rng_seeds(parameter_study.get_rng(), mio::LogLevel::warn); TempFileRegister tmp_file_register; @@ -317,10 +320,13 @@ TEST(TestSaveResult, save_percentiles_and_sums) auto ensemble_edges = std::vector>>{}; ensemble_edges.reserve(size_t(num_runs)); parameter_study.run( - [](auto&& g) { - return draw_sample(g); + [](auto&& g, auto t0_, auto dt_, auto) { + auto copy = g; + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); }, - [&](auto results_graph, auto /*run_idx*/) { + [&](auto&& results, auto /*run_idx*/) { + auto results_graph = results.get_graph(); ensemble_results.push_back(mio::interpolate_simulation_result(results_graph)); ensemble_params.emplace_back(); diff --git a/cpp/tests/test_utils.cpp b/cpp/tests/test_utils.cpp index cbffefdfff..f5283a81b8 100644 --- a/cpp/tests/test_utils.cpp +++ b/cpp/tests/test_utils.cpp @@ -17,6 +17,7 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include "memilio/utils/base_dir.h" #include "memilio/utils/index.h" #include "memilio/utils/index_range.h" #include "memilio/utils/logging.h" @@ -27,6 +28,8 @@ #include #include +#include + template struct CategoryTag : public mio::Index> { CategoryTag(size_t value) @@ -160,3 +163,13 @@ TEST(TestUtils, RedirectLogger) logger.release(); } + +TEST(TestUtils, base_dir) +{ + auto base_dir = boost::filesystem::path(mio::base_dir()); + // check that the path exists + EXPECT_TRUE(boost::filesystem::exists(base_dir)); + // check that the path is correct, by sampling some fixed paths from project files + EXPECT_TRUE(boost::filesystem::exists(base_dir / "cpp" / "memilio")); + EXPECT_TRUE(boost::filesystem::exists(base_dir / "pycode" / "memilio-epidata")); +} diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index 273da48f2d..881c93b088 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp @@ -48,6 +48,8 @@ #include "pybind11/pybind11.h" #include "pybind11/stl_bind.h" #include "Eigen/Core" +#include +#include #include namespace py = pybind11; @@ -55,91 +57,66 @@ namespace py = pybind11; namespace { -//select only the first node of the graph of each run, used for parameterstudy with single nodes -template -std::vector filter_graph_results( - std::vector, mio::MobilityEdge>>&& graph_results) -{ - std::vector results; - results.reserve(graph_results.size()); - for (auto i = size_t(0); i < graph_results.size(); ++i) { - results.emplace_back(std::move(graph_results[i].nodes()[0].property.get_simulation())); - } - return std::move(results); -} - /* * @brief bind ParameterStudy for any model */ -template +template void bind_ParameterStudy(py::module_& m, std::string const& name) { - pymio::bind_class, pymio::EnablePickling::Never>(m, name.c_str()) - .def(py::init(), py::arg("model"), py::arg("t0"), - py::arg("tmax"), py::arg("num_runs")) - .def(py::init>&, double, double, - double, size_t>(), - py::arg("model_graph"), py::arg("t0"), py::arg("tmax"), py::arg("dt"), py::arg("num_runs")) - .def_property("num_runs", &mio::ParameterStudy::get_num_runs, - &mio::ParameterStudy::set_num_runs) - .def_property("tmax", &mio::ParameterStudy::get_tmax, - &mio::ParameterStudy::set_tmax) - .def_property("t0", &mio::ParameterStudy::get_t0, - &mio::ParameterStudy::set_t0) - .def_property_readonly("model", py::overload_cast<>(&mio::ParameterStudy::get_model), - py::return_value_policy::reference_internal) - .def_property_readonly("model", - py::overload_cast<>(&mio::ParameterStudy::get_model, py::const_), + using GraphT = mio::Graph, mio::MobilityEdge>; + using SimulationT = mio::GraphSimulation; + using ParametersT = mio::Graph>; + using StudyT = mio::ParameterStudy; + + const auto create_simulation = [](const ParametersT& g, double t0, double dt, size_t) { + auto copy = g; + return mio::make_sampled_graph_simulation(draw_sample(copy), t0, dt, dt); + }; + + pymio::bind_class(m, name.c_str()) + .def(py::init(), py::arg("parameters"), py::arg("t0"), + py::arg("tmax"), py::arg("dt"), py::arg("num_runs")) + .def_property_readonly("num_runs", &StudyT::get_num_runs) + .def_property_readonly("tmax", &StudyT::get_tmax) + .def_property_readonly("t0", &StudyT::get_t0) + .def_property_readonly("dt", &StudyT::get_dt) + .def_property_readonly("model_graph", py::overload_cast<>(&StudyT::get_parameters), py::return_value_policy::reference_internal) - .def_property_readonly("model_graph", - py::overload_cast<>(&mio::ParameterStudy::get_model_graph), + .def_property_readonly("model_graph", py::overload_cast<>(&StudyT::get_parameters, py::const_), py::return_value_policy::reference_internal) - .def_property_readonly( - "model_graph", py::overload_cast<>(&mio::ParameterStudy::get_model_graph, py::const_), - py::return_value_policy::reference_internal) .def( "run", - [](mio::ParameterStudy& self, - std::function, mio::MobilityEdge>, - size_t)> - handle_result) { - self.run( - [](auto&& g) { - return draw_sample(g); - }, - [&handle_result](auto&& g, auto&& run_idx) { - //handle_result_function needs to return something - //we don't want to run an unknown python object through parameterstudies, so - //we just return 0 and ignore the list returned by run(). - //So python will behave slightly different than c++ - handle_result(std::move(g), run_idx); - return 0; - }); + [&create_simulation](StudyT& self, std::function handle_result) { + self.run_serial(create_simulation, [&handle_result](auto&& g, auto&& run_idx) { + //handle_result_function needs to return something + //we don't want to run an unknown python object through parameterstudies, so + //we just return 0 and ignore the list returned by run(). + //So python will behave slightly different than c++ + handle_result(std::move(g.get_graph()), run_idx); + return 0; + }); }, py::arg("handle_result_func")) .def("run", - [](mio::ParameterStudy& self) { //default argument doesn't seem to work with functions - return self.run([](auto&& g) { - return draw_sample(g); + [&create_simulation](StudyT& self) { //default argument doesn't seem to work with functions + return self.run_serial(create_simulation, [](SimulationT&& result, size_t) { + return std::move(result.get_graph()); }); }) .def( "run_single", - [](mio::ParameterStudy& self, std::function handle_result) { - self.run( - [](auto&& g) { - return draw_sample(g); - }, - [&handle_result](auto&& r, auto&& run_idx) { - handle_result(std::move(r.nodes()[0].property.get_simulation()), run_idx); - return 0; - }); + [&create_simulation](StudyT& self, std::function handle_result) { + self.run_serial(create_simulation, [&handle_result](auto&& r, auto&& run_idx) { + handle_result(std::move(r.get_graph().nodes()[0].property.get_simulation()), run_idx); + return 0; + }); }, py::arg("handle_result_func")) - .def("run_single", [](mio::ParameterStudy& self) { - return filter_graph_results(self.run([](auto&& g) { - return draw_sample(g); - })); + .def("run_single", [&create_simulation](StudyT& self) { + return self.run_serial(create_simulation, [](SimulationT&& result, size_t) { + //select only the first node of the graph of each run, used for parameterstudy with single nodes + return std::move(result.get_graph().nodes()[0].property.get_simulation()); + }); }); } @@ -287,11 +264,11 @@ PYBIND11_MODULE(_simulation_osecir, m) mio::osecir::InfectionState::InfectedSymptoms, mio::osecir::InfectionState::Recovered}; auto weights = std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}; auto result = mio::set_edges, mio::MobilityParameters, - mio::MobilityCoefficientGroup, mio::osecir::InfectionState, - decltype(mio::read_mobility_plain)>(mobility_data_file, params_graph, - mobile_comp, contact_locations_size, - mio::read_mobility_plain, weights); + ContactLocation, mio::osecir::Model, mio::MobilityParameters, + mio::MobilityCoefficientGroup, mio::osecir::InfectionState, + decltype(mio::read_mobility_plain)>(mobility_data_file, params_graph, + mobile_comp, contact_locations_size, + mio::read_mobility_plain, weights); return pymio::check_and_throw(result); }, py::return_value_policy::move); diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index 8e1e922404..9c276310b8 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -49,98 +49,94 @@ namespace py = pybind11; namespace { -//select only the first node of the graph of each run, used for parameterstudy with single nodes -template -std::vector filter_graph_results( - std::vector, mio::MobilityEdge>>&& graph_results) -{ - std::vector results; - results.reserve(graph_results.size()); - for (auto i = size_t(0); i < graph_results.size(); ++i) { - results.emplace_back(std::move(graph_results[i].nodes()[0].property.get_simulation())); - } - return std::move(results); -} /* * @brief bind ParameterStudy for any model */ -template +template void bind_ParameterStudy(py::module_& m, std::string const& name) { - pymio::bind_class, pymio::EnablePickling::Never>(m, name.c_str()) - .def(py::init(), py::arg("model"), py::arg("t0"), - py::arg("tmax"), py::arg("num_runs")) - .def(py::init>&, double, double, - double, size_t>(), - py::arg("model_graph"), py::arg("t0"), py::arg("tmax"), py::arg("dt"), py::arg("num_runs")) - .def_property("num_runs", &mio::ParameterStudy::get_num_runs, - &mio::ParameterStudy::set_num_runs) - .def_property("tmax", &mio::ParameterStudy::get_tmax, - &mio::ParameterStudy::set_tmax) - .def_property("t0", &mio::ParameterStudy::get_t0, - &mio::ParameterStudy::set_t0) - .def_property_readonly("model", py::overload_cast<>(&mio::ParameterStudy::get_model), - py::return_value_policy::reference_internal) - .def_property_readonly("model", - py::overload_cast<>(&mio::ParameterStudy::get_model, py::const_), + using GraphT = mio::Graph, mio::MobilityEdge>; + using SimulationT = mio::GraphSimulation; + using ParametersT = mio::Graph>; + using StudyT = mio::ParameterStudy; + + pymio::bind_class(m, name.c_str()) + .def(py::init(), py::arg("parameters"), py::arg("t0"), + py::arg("tmax"), py::arg("dt"), py::arg("num_runs")) + .def_property_readonly("num_runs", &StudyT::get_num_runs) + .def_property_readonly("tmax", &StudyT::get_tmax) + .def_property_readonly("t0", &StudyT::get_t0) + .def_property_readonly("dt", &StudyT::get_dt) + .def_property_readonly("model_graph", py::overload_cast<>(&StudyT::get_parameters), py::return_value_policy::reference_internal) - .def_property_readonly("model_graph", - py::overload_cast<>(&mio::ParameterStudy::get_model_graph), + .def_property_readonly("model_graph", py::overload_cast<>(&StudyT::get_parameters, py::const_), py::return_value_policy::reference_internal) - .def_property_readonly( - "model_graph", py::overload_cast<>(&mio::ParameterStudy::get_model_graph, py::const_), - py::return_value_policy::reference_internal) .def( "run", - [](mio::ParameterStudy& self, - std::function, mio::MobilityEdge>, - size_t)> + [](StudyT& self, + std::function, mio::MobilityEdge>, size_t)> handle_result, bool variant_high) { - self.run( - [variant_high](auto&& g) { - return draw_sample(g, variant_high); + self.run_serial( + [variant_high](const ParametersT& g, double t0, double dt, size_t) { + auto copy = g; + return mio::make_sampled_graph_simulation(draw_sample(copy, variant_high), t0, dt, + dt); }, [&handle_result](auto&& g, auto&& run_idx) { //handle_result_function needs to return something //we don't want to run an unknown python object through parameterstudies, so //we just return 0 and ignore the list returned by run(). //So python will behave slightly different than c++ - handle_result(std::move(g), run_idx); + handle_result(std::move(g.get_graph()), run_idx); return 0; }); }, py::arg("handle_result_func"), py::arg("variant_high")) .def( "run", - [](mio::ParameterStudy& self, + [](StudyT& self, bool variant_high) { //default argument doesn't seem to work with functions - return self.run([variant_high](auto&& g) { - return draw_sample(g, variant_high); - }); + return self.run_serial( + [variant_high](const ParametersT& g, double t0, double dt, size_t) { + auto copy = g; + return mio::make_sampled_graph_simulation(draw_sample(copy, variant_high), t0, dt, + dt); + }, + [](SimulationT&& result, size_t) { + return std::move(result.get_graph()); + }); }, py::arg("variant_high")) .def( "run_single", - [](mio::ParameterStudy& self, std::function handle_result, - bool variant_high) { - self.run( - [variant_high](auto&& g) { - return draw_sample(g, variant_high); + [](StudyT& self, std::function handle_result, bool variant_high) { + self.run_serial( + [variant_high](const ParametersT& g, double t0, double dt, size_t) { + auto copy = g; + return mio::make_sampled_graph_simulation(draw_sample(copy, variant_high), t0, dt, + dt); }, [&handle_result](auto&& r, auto&& run_idx) { - handle_result(std::move(r.nodes()[0].property.get_simulation()), run_idx); + handle_result(std::move(r.get_graph().nodes()[0].property.get_simulation()), run_idx); return 0; }); }, py::arg("handle_result_func"), py::arg("variant_high")) .def( "run_single", - [](mio::ParameterStudy& self, bool variant_high) { - return filter_graph_results(self.run([variant_high](auto&& g) { - return draw_sample(g, variant_high); - })); + [](StudyT& self, bool variant_high) { + return self.run_serial( + [variant_high](const ParametersT& g, double t0, double dt, size_t) { + auto copy = g; + return mio::make_sampled_graph_simulation(draw_sample(copy, variant_high), t0, dt, + dt); + }, + [](SimulationT&& result, size_t) { + //select only the first node of the graph of each run, used for parameterstudy with single nodes + return std::move(result.get_graph().nodes()[0].property.get_simulation()); + }); }, py::arg("variant_high")); } @@ -340,9 +336,9 @@ PYBIND11_MODULE(_simulation_osecirvvs, m) mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity}; auto weights = std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}; auto result = mio::set_edges, - mio::MobilityParameters, mio::MobilityCoefficientGroup, - mio::osecirvvs::InfectionState, decltype(mio::read_mobility_plain)>( + ContactLocation, mio::osecirvvs::Model, + mio::MobilityParameters, mio::MobilityCoefficientGroup, + mio::osecirvvs::InfectionState, decltype(mio::read_mobility_plain)>( mobility_data_file, params_graph, mobile_comp, contact_locations_size, mio::read_mobility_plain, weights); return pymio::check_and_throw(result); diff --git a/pycode/memilio-simulation/memilio/simulation_test/test_parameter_study.py b/pycode/memilio-simulation/memilio/simulation_test/test_parameter_study.py index c636e77a7c..2c5a2c1414 100644 --- a/pycode/memilio-simulation/memilio/simulation_test/test_parameter_study.py +++ b/pycode/memilio-simulation/memilio/simulation_test/test_parameter_study.py @@ -85,12 +85,14 @@ def test_graph(self): def test_run(self): """ """ - model = self._get_model() + graph = osecir.ModelGraph() + graph.add_node(0, self._get_model()) t0 = 1 tmax = 10 + dt = 0.1 num_runs = 3 - study = osecir.ParameterStudy(model, t0, tmax, num_runs) + study = osecir.ParameterStudy(graph, t0, tmax, dt, num_runs) self.assertEqual(study.t0, t0) self.assertEqual(study.tmax, tmax)