From a871ba6aeb270d31e119840e813dad4f766709f0 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Wed, 3 Sep 2025 15:16:38 +0200 Subject: [PATCH 01/29] first working example --- cpp/examples/abm_minimal.cpp | 117 ++++++++++++++++++- cpp/memilio/compartments/parameter_studies.h | 4 +- cpp/models/abm/model.h | 2 + cpp/models/abm/simulation.h | 11 +- 4 files changed, 124 insertions(+), 10 deletions(-) diff --git a/cpp/examples/abm_minimal.cpp b/cpp/examples/abm_minimal.cpp index a78e0dddf7..d46618f5b8 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/cpp/examples/abm_minimal.cpp @@ -17,12 +17,72 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include "abm/simulation.h" +#include "abm/common_abm_loggers.h" +#include + +// shim class to make the abm simulation look like an ode simulation +class ResultSim : public mio::abm::Simulation +{ +public: + using Model = mio::abm::Model; + + ResultSim(const Model& m, double t, double) + : mio::abm::Simulation(mio::abm::TimePoint{mio::abm::days(t).seconds()}, Model(m)) + { + history.log(*this); // set initial results + } + + void advance(double tmax) + { + mio::abm::Simulation::advance(mio::abm::TimePoint{mio::abm::days(tmax).seconds()}, history); + } + + const mio::TimeSeries& get_result() const + { + return get<0>(history.get_log()); + } + + mio::History history{ + Eigen::Index(mio::abm::InfectionState::Count)}; +}; + +// graph ode specific function, overload with noop to avoid compiler errors +template +void calculate_mobility_returns(Eigen::Ref::Vector>, const ResultSim&, + Eigen::Ref::Vector>, double, double) +{ +} +#include "memilio/mobility/metapopulation_mobility_instant.h" + +namespace mio +{ + +// overload for ResultSim without mobility +auto make_mobility_sim(double t0, double dt, Graph, MobilityEdge>&& graph) +{ + using GraphT = GraphSimulation, MobilityEdge>, double, double, + void (*)(double, double, mio::MobilityEdge&, mio::SimulationNode&, + mio::SimulationNode&), + void (*)(double, double, mio::SimulationNode&)>; + return GraphT(t0, dt, std::move(graph), + static_cast&)>(&advance_model), + [](double, double, MobilityEdge&, SimulationNode&, SimulationNode&) {}); +} + +} // namespace mio + #include "abm/household.h" #include "abm/lockdown_rules.h" #include "abm/model.h" -#include "abm/common_abm_loggers.h" +#include "abm/time.h" +#include "memilio/compartments/parameter_studies.h" +#include "memilio/data/analyze_result.h" +#include "memilio/io/result_io.h" #include "memilio/utils/abstract_parameter_distribution.h" +#include "memilio/utils/time_series.h" +#include #include int main() @@ -154,14 +214,65 @@ int main() // Set start and end time for the simulation. auto t0 = mio::abm::TimePoint(0); auto tmax = t0 + mio::abm::days(10); - auto sim = mio::abm::Simulation(t0, std::move(model)); + // auto sim = mio::abm::Simulation(t0, std::move(model)); + const size_t num_runs = 10; // Create a history object to store the time series of the infection states. mio::History historyTimeSeries{ Eigen::Index(mio::abm::InfectionState::Count)}; + mio::ParameterStudy study(model, t0.days(), tmax.days(), num_runs); + + auto ensemble = study.run( + [](auto&& graph) { + // TODO: some actual sampling so that the percentiles show anything + return graph; + }, + [&](auto results_graph, auto&& /* run_idx */) { + auto interpolated_result = mio::interpolate_simulation_result(results_graph); + + auto params = std::vector{}; + params.reserve(results_graph.nodes().size()); + std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), std::back_inserter(params), + [](auto&& node) { + return node.property.get_simulation().get_model(); + }); + return std::make_pair(std::move(interpolated_result), std::move(params)); + }); + if (ensemble.size() > 0) { + auto ensemble_results = std::vector>>{}; + ensemble_results.reserve(ensemble.size()); + auto ensemble_params = std::vector>{}; + ensemble_params.reserve(ensemble.size()); + for (auto&& run : ensemble) { + ensemble_results.emplace_back(std::move(run.first)); + ensemble_params.emplace_back(std::move(run.second)); + } + + 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); + + const boost::filesystem::path save_dir = "/home/schm_r6/Code/memilio/memilio/results/"; + + if (!save_dir.empty()) { + + mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups, + (save_dir / ("Results_" + std::string("p05") + ".h5")).string())); + mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups, + (save_dir / ("Results_" + std::string("p25") + ".h5")).string())); + mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups, + (save_dir / ("Results_" + std::string("p50") + ".h5")).string())); + mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups, + (save_dir / ("Results_" + std::string("p75") + ".h5")).string())); + mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups, + (save_dir / ("Results_" + std::string("p95") + ".h5")).string())); + } + } // Run the simulation until tmax with the history object. - sim.advance(tmax, historyTimeSeries); + // sim.advance(tmax, historyTimeSeries); // The results are written into the file "abm_minimal.txt" as a table with 9 columns. // The first column is Time. The other columns correspond to the number of people with a certain infection state at this Time: diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 62cfe564a9..c4ae6d0132 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -43,7 +43,7 @@ namespace mio * 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. */ -template +template > class ParameterStudy { public: @@ -60,7 +60,7 @@ class ParameterStudy * The Graph type that stores simulations and their results of each run. * This is the output of ParameterStudies for each run. */ - using SimulationGraph = mio::Graph, mio::MobilityEdge>; + using SimulationGraph = mio::Graph, EdgeT>; /** * create study for graph of compartment models. 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/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. }; From 7f69da2b238f6c77619c31ad0854324eb3b9182e Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 9 Oct 2025 13:19:19 +0200 Subject: [PATCH 02/29] [wip] parameter studies v2 --- cpp/examples/abm_minimal.cpp | 117 ++++----- cpp/memilio/compartments/parameter_studies.h | 246 +++++++++++++++++++ 2 files changed, 297 insertions(+), 66 deletions(-) diff --git a/cpp/examples/abm_minimal.cpp b/cpp/examples/abm_minimal.cpp index 63c8b7776b..72ce3871ee 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/cpp/examples/abm_minimal.cpp @@ -20,6 +20,19 @@ #include "abm/simulation.h" #include "abm/common_abm_loggers.h" #include +#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/result_io.h" +#include "memilio/utils/abstract_parameter_distribution.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/time_series.h" + +#include +#include // shim class to make the abm simulation look like an ode simulation class ResultSim : public mio::abm::Simulation @@ -33,9 +46,15 @@ class ResultSim : public mio::abm::Simulation history.log(*this); // set initial results } - void advance(double tmax) + ResultSim(const Model& m, mio::abm::TimePoint t) + : mio::abm::Simulation(t, Model(m)) + { + history.log(*this); // set initial results + } + + void advance(mio::abm::TimePoint tmax) { - mio::abm::Simulation::advance(mio::abm::TimePoint{mio::abm::days(tmax).seconds()}, history); + mio::abm::Simulation::advance(tmax, history); } const mio::TimeSeries& get_result() const @@ -47,46 +66,9 @@ class ResultSim : public mio::abm::Simulation Eigen::Index(mio::abm::InfectionState::Count)}; }; -// graph ode specific function, overload with noop to avoid compiler errors -template -void calculate_mobility_returns(Eigen::Ref::Vector>, const ResultSim&, - Eigen::Ref::Vector>, double, double) -{ -} -#include "memilio/mobility/metapopulation_mobility_instant.h" - -namespace mio -{ - -// overload for ResultSim without mobility -auto make_mobility_sim(double t0, double dt, Graph, MobilityEdge>&& graph) -{ - using GraphT = GraphSimulation, MobilityEdge>, double, double, - void (*)(double, double, mio::MobilityEdge&, mio::SimulationNode&, - mio::SimulationNode&), - void (*)(double, double, mio::SimulationNode&)>; - return GraphT(t0, dt, std::move(graph), - static_cast&)>(&advance_model), - [](double, double, MobilityEdge&, SimulationNode&, SimulationNode&) {}); -} - -} // namespace mio - -#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/result_io.h" -#include "memilio/utils/abstract_parameter_distribution.h" -#include "memilio/utils/time_series.h" - -#include -#include - int main() { + mio::mpi::init(); // This is a minimal example with children and adults < 60 year old. // We divided them into 4 different age groups, which are defined as follows: mio::set_log_level(mio::LogLevel::warn); @@ -95,7 +77,6 @@ int main() 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); // Set same infection parameter for all age groups. For example, the incubation period is log normally distributed with parameters 4 and 1. @@ -111,7 +92,7 @@ int main() model.parameters.check_constraints(); // There are 10 households for each household group. - int n_households = 10; + int n_households = 100; // 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. @@ -213,40 +194,42 @@ int main() // Set start and end time for the simulation. auto t0 = mio::abm::TimePoint(0); - auto tmax = t0 + mio::abm::days(10); + auto tmax = t0 + mio::abm::days(100); // auto sim = mio::abm::Simulation(t0, std::move(model)); - const size_t num_runs = 10; + const size_t num_runs = 100; // Create a history object to store the time series of the infection states. mio::History historyTimeSeries{ Eigen::Index(mio::abm::InfectionState::Count)}; - mio::ParameterStudy study(model, t0.days(), tmax.days(), num_runs); + mio::ParameterStudy2 study( + model, t0, tmax, mio::abm::TimeSpan(0), num_runs); auto ensemble = study.run( - [](auto&& graph) { - // TODO: some actual sampling so that the percentiles show anything - return graph; + [](auto&& model_, auto t0_, auto) { + // TODO: change the parameters around a bit? + auto sim = ResultSim(model_, t0_); + sim.get_model().get_rng().generate_seeds(); + return sim; }, - [&](auto results_graph, auto&& /* run_idx */) { - auto interpolated_result = mio::interpolate_simulation_result(results_graph); + [](auto&& sim, auto&& /* run_idx */) { + auto interpolated_result = mio::interpolate_simulation_result(sim.get_result()); auto params = std::vector{}; - params.reserve(results_graph.nodes().size()); - std::transform(results_graph.nodes().begin(), results_graph.nodes().end(), std::back_inserter(params), - [](auto&& node) { - return node.property.get_simulation().get_model(); - }); - return std::make_pair(std::move(interpolated_result), std::move(params)); + return std::make_pair(std::move(interpolated_result), sim.get_model()); }); + if (ensemble.size() > 0) { auto ensemble_results = std::vector>>{}; ensemble_results.reserve(ensemble.size()); auto ensemble_params = std::vector>{}; ensemble_params.reserve(ensemble.size()); + int j = 0; + std::cout << "## " << ensemble[0].first.get_value(0).transpose() << "\n"; for (auto&& run : ensemble) { - ensemble_results.emplace_back(std::move(run.first)); - ensemble_params.emplace_back(std::move(run.second)); + std::cout << j++ << " " << run.first.get_last_value().transpose() << "\n"; + ensemble_results.emplace_back(std::vector{std::move(run.first)}); + ensemble_params.emplace_back(std::vector{std::move(run.second)}); } auto ensemble_results_p05 = ensemble_percentile(ensemble_results, 0.05); @@ -274,14 +257,16 @@ int main() // Run the simulation until tmax with the history object. // sim.advance(tmax, historyTimeSeries); - // The results are written into the file "abm_minimal.txt" as a table with 9 columns. - // The first column is Time. The other columns correspond to the number of people with a certain infection state at this Time: - // Time = Time in days, S = Susceptible, E = Exposed, I_NS = InfectedNoSymptoms, I_Sy = InfectedSymptoms, I_Sev = InfectedSevere, - // I_Crit = InfectedCritical, R = Recovered, D = Dead - std::ofstream outfile("abm_minimal.txt"); - std::get<0>(historyTimeSeries.get_log()) - .print_table(outfile, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4); - std::cout << "Results written to abm_minimal.txt" << std::endl; + // // The results are written into the file "abm_minimal.txt" as a table with 9 columns. + // // The first column is Time. The other columns correspond to the number of people with a certain infection state at this Time: + // // Time = Time in days, S = Susceptible, E = Exposed, I_NS = InfectedNoSymptoms, I_Sy = InfectedSymptoms, I_Sev = InfectedSevere, + // // I_Crit = InfectedCritical, R = Recovered, D = Dead + // std::ofstream outfile("abm_minimal.txt"); + // std::get<0>(historyTimeSeries.get_log()) + // .print_table(outfile, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4); + // std::cout << "Results written to abm_minimal.txt" << std::endl; + + mio::mpi::finalize(); return 0; } diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 98819b03d3..68dc1865fd 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -23,17 +23,22 @@ #include "memilio/io/binary_serializer.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 namespace mio { @@ -383,6 +388,247 @@ class ParameterStudy RandomNumberGenerator m_rng; }; +/** + * 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. + */ +template +class ParameterStudy2 +{ +public: + using Simulation = SimulationType; + using Parameters = ParameterType; + using Time = TimeType; + using Step = StepType; + + // 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 + */ + ParameterStudy2(const ParameterType& global_parameters, Time t0, Time tmax, Step dt, size_t num_runs) + : m_parameters(global_parameters) + , m_num_runs(num_runs) + , m_t0{t0} + , m_tmax{tmax} + , m_dt(dt) + { + } + + // /** + // * @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 + // */ + // template >> + // ParameterStudy2(typename Simulation::Model const& model, Time t0, Time tmax, Step dt, size_t num_runs) + // : ParameterStudy2(model, t0, tmax, dt, num_runs) + // { + // // TODO how is this supposed to work wrt. model? is this just a special case where ParameterType=Model? + // } + + /** + * @brief + * @param sample_simulation A function that accepts ParameterType and returns an instance of SimulationType. + * @param process_simulation_result A function that accepts S + */ + template + std::vector> + run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) + { + static_assert(std::is_invocable_r_v, + "Incorrect Type for create_simulation."); + static_assert(std::is_invocable_v, + "Incorrect Type for process_simulation_result."); + int num_procs, rank; +#ifdef MEMILIO_ENABLE_MPI + MPI_Comm_size(mpi::get_world(), &num_procs); + MPI_Comm_rank(mpi::get_world(), &rank); +#else + num_procs = 1; + rank = 0; +#endif + + //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 + m_rng.synchronize(); + thread_local_rng() = m_rng; + + auto run_distribution = distribute_runs(m_num_runs, num_procs); + auto 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_simulation(m_parameters, m_t0, m_dt); + 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)); + } + + //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))); + +#ifdef MEMILIO_ENABLE_MPI + //gather results + if (rank == 0) { + for (int src_rank = 1; src_rank < num_procs; ++src_rank) { + int bytes_size; + MPI_Recv(&bytes_size, 1, MPI_INT, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE); + 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{}); + if (!src_ensemble_results) { + log_error("Error receiving ensemble results from rank {}.", src_rank); + } + std::copy(src_ensemble_results.value().begin(), src_ensemble_results.value().end(), + std::back_inserter(ensemble_result)); + } + } + else { + auto bytes = serialize_binary(ensemble_result); + auto 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 + } +#endif + + 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(CreateSimulationFunction&& create_simulation) + { + return run(std::forward(create_simulation), + [](Simulation&& sim, size_t) -> Simulation&& { + return sim; + }); + } + + /** + * @brief returns the number of Monte Carlo runs + */ + size_t get_num_runs() const + { + return m_num_runs; + } + + /** + * @brief returns end point in simulation + */ + Time get_tmax() const + { + return m_tmax; + } + + /** + * @brief returns start point in simulation + */ + Time get_t0() const + { + return m_t0; + } + /** + * Get the input graph that the parameter study is run for. + * Use for graph simulations, use get_model for single node simulations. + * @{ + */ + const Parameters& get_parameters() const + { + return m_parameters; + } + /** @} */ + + RandomNumberGenerator& get_rng() + { + return m_rng; + } + +private: + 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; + + std::vector run_distribution(num_procs); + std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); + std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local); + + return run_distribution; + } + +private: + // Stores Graph with the names and ranges of all parameters + ParameterType m_parameters; + + size_t m_num_runs; + + // Start time (should be the same for all simulations) + Time m_t0; + // End time (should be the same for all simulations) + Time m_tmax; + // adaptive time step of the integrator (will be corrected if too large/small) + Step m_dt; + // + RandomNumberGenerator m_rng; +}; + +// //sample parameters and create simulation +// template +// GraphSimulation create_sampled_simulation(auto&& params_graph, +// SampleGraphFunction sample_graph) +// { +// SimulationGraph sim_graph; + +// auto sampled_graph = sample_graph(params_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); +// } + +// return make_mobility_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); +// } + } // namespace mio #endif // MIO_COMPARTMENTS_PARAMETER_STUDIES_H From eaaa5c2e52863e76575444dff0de15c3c9a6edee Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 9 Oct 2025 13:34:29 +0200 Subject: [PATCH 03/29] [wip] add run_idx to create_sim function --- cpp/examples/abm_minimal.cpp | 5 +++-- cpp/memilio/compartments/parameter_studies.h | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/cpp/examples/abm_minimal.cpp b/cpp/examples/abm_minimal.cpp index 72ce3871ee..5c90624b8b 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/cpp/examples/abm_minimal.cpp @@ -32,6 +32,7 @@ #include "memilio/utils/time_series.h" #include +#include #include // shim class to make the abm simulation look like an ode simulation @@ -206,10 +207,10 @@ int main() model, t0, tmax, mio::abm::TimeSpan(0), num_runs); auto ensemble = study.run( - [](auto&& model_, auto t0_, auto) { + [](auto&& model_, auto t0_, auto, size_t run_idx) { // TODO: change the parameters around a bit? auto sim = ResultSim(model_, t0_); - sim.get_model().get_rng().generate_seeds(); + sim.get_model().get_rng().set_counter(mio::Counter(run_idx)); return sim; }, [](auto&& sim, auto&& /* run_idx */) { diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 68dc1865fd..fcbc2ea8db 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -444,7 +444,7 @@ class ParameterStudy2 std::vector> run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) { - static_assert(std::is_invocable_r_v, + static_assert(std::is_invocable_r_v, "Incorrect Type for create_simulation."); static_assert(std::is_invocable_v, "Incorrect Type for process_simulation_result."); @@ -482,7 +482,7 @@ class ParameterStudy2 thread_local_rng().set_counter(run_rng_counter); //sample - auto sim = create_simulation(m_parameters, m_t0, m_dt); + auto sim = create_simulation(m_parameters, m_t0, m_dt, run_idx); log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", (thread_local_rng().get_counter() - run_rng_counter).get()); From 280afa5824257b630748799bfbb9eb48c0d89fff Mon Sep 17 00:00:00 2001 From: Sascha Korf <51127093+xsaschako@users.noreply.github.com> Date: Thu, 9 Oct 2025 15:03:31 +0200 Subject: [PATCH 04/29] mpi wokrks --- cpp/examples/abm_minimal.cpp | 116 ++++++++++++++++++++--------------- 1 file changed, 68 insertions(+), 48 deletions(-) diff --git a/cpp/examples/abm_minimal.cpp b/cpp/examples/abm_minimal.cpp index 5c90624b8b..df6392ec75 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/cpp/examples/abm_minimal.cpp @@ -35,6 +35,8 @@ #include #include +constexpr size_t num_age_groups = 4; + // shim class to make the abm simulation look like an ode simulation class ResultSim : public mio::abm::Simulation { @@ -67,19 +69,23 @@ class ResultSim : public mio::abm::Simulation Eigen::Index(mio::abm::InfectionState::Count)}; }; -int main() +mio::abm::Model make_model(uint32_t run_idx) { - mio::mpi::init(); - // This is a minimal example with children and adults < 60 year old. - // We divided them into 4 different age groups, which are defined as follows: - mio::set_log_level(mio::LogLevel::warn); - size_t num_age_groups = 4; + 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); + auto model = mio::abm::Model(num_age_groups); + std::initializer_list seeds = {14159265, 243141u + run_idx}; + auto rng = mio::RandomNumberGenerator(); + rng.seed(seeds); + auto run_rng_counter = + mio::rng_totalsequence_counter(static_cast(run_idx), mio::Counter(0)); + rng.set_counter(run_rng_counter); + 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.); @@ -164,9 +170,9 @@ int main() 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 rng = mio::abm::PersonalRandomNumberGenerator(person); + auto person_rng = mio::abm::PersonalRandomNumberGenerator(person); if (infection_state != mio::abm::InfectionState::Susceptible) { - person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(), + person.add_new_infection(mio::abm::Infection(person_rng, mio::abm::VirusVariant::Wildtype, person.get_age(), model.parameters, start_date, infection_state)); } } @@ -193,31 +199,45 @@ int main() 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(); + // This is a minimal example with children and adults < 60 year old. + // We divided them into 4 different age groups, which are defined as follows: + + 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(100); + auto tmax = t0 + mio::abm::days(5); // auto sim = mio::abm::Simulation(t0, std::move(model)); - const size_t num_runs = 100; + const size_t num_runs = 10; // Create a history object to store the time series of the infection states. mio::History historyTimeSeries{ Eigen::Index(mio::abm::InfectionState::Count)}; - mio::ParameterStudy2 study( - model, t0, tmax, mio::abm::TimeSpan(0), num_runs); + mio::ParameterStudy2 study( + 1, t0, tmax, mio::abm::TimeSpan(0), num_runs); auto ensemble = study.run( - [](auto&& model_, auto t0_, auto, size_t run_idx) { + [](auto&&, auto t0_, auto, size_t run_idx) { // TODO: change the parameters around a bit? - auto sim = ResultSim(model_, t0_); - sim.get_model().get_rng().set_counter(mio::Counter(run_idx)); + + auto sim = ResultSim(make_model(run_idx), t0_); return sim; }, - [](auto&& sim, auto&& /* run_idx */) { + [](auto&& sim, auto&& run_idx) { auto interpolated_result = mio::interpolate_simulation_result(sim.get_result()); - + std::string name_run = "abm_minimal_run_" + std::to_string(run_idx) + ".txt"; + std::ofstream outfile_run(name_run); + 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 " << name_run << std::endl; auto params = std::vector{}; - return std::make_pair(std::move(interpolated_result), sim.get_model()); + return std::move(interpolated_result); }); if (ensemble.size() > 0) { @@ -225,35 +245,35 @@ int main() ensemble_results.reserve(ensemble.size()); auto ensemble_params = std::vector>{}; ensemble_params.reserve(ensemble.size()); - int j = 0; - std::cout << "## " << ensemble[0].first.get_value(0).transpose() << "\n"; - for (auto&& run : ensemble) { - std::cout << j++ << " " << run.first.get_last_value().transpose() << "\n"; - ensemble_results.emplace_back(std::vector{std::move(run.first)}); - ensemble_params.emplace_back(std::vector{std::move(run.second)}); - } - - 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); - - const boost::filesystem::path save_dir = "/home/schm_r6/Code/memilio/memilio/results/"; - - if (!save_dir.empty()) { - - mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups, - (save_dir / ("Results_" + std::string("p05") + ".h5")).string())); - mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups, - (save_dir / ("Results_" + std::string("p25") + ".h5")).string())); - mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups, - (save_dir / ("Results_" + std::string("p50") + ".h5")).string())); - mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups, - (save_dir / ("Results_" + std::string("p75") + ".h5")).string())); - mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups, - (save_dir / ("Results_" + std::string("p95") + ".h5")).string())); - } + // int j = 0; + // std::cout << "## " << ensemble[0].first.get_value(0).transpose() << "\n"; + // for (auto&& run : ensemble) { + // std::cout << j++ << " " << run.first.get_last_value().transpose() << "\n"; + // ensemble_results.emplace_back(std::vector{std::move(run.first)}); + // ensemble_params.emplace_back(std::vector{std::move(run.second)}); + // } + + // 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); + + // const boost::filesystem::path save_dir = "/Users/saschakorf/Nosynch/Arbeit/memilio/cpp/examples/results"; + + // if (!save_dir.empty()) { + + // mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups, + // (save_dir / ("Results_" + std::string("p05") + ".h5")).string())); + // mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups, + // (save_dir / ("Results_" + std::string("p25") + ".h5")).string())); + // mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups, + // (save_dir / ("Results_" + std::string("p50") + ".h5")).string())); + // mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups, + // (save_dir / ("Results_" + std::string("p75") + ".h5")).string())); + // mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups, + // (save_dir / ("Results_" + std::string("p95") + ".h5")).string())); + // } } // Run the simulation until tmax with the history object. // sim.advance(tmax, historyTimeSeries); From f52035d0471244ea80c7db86cdeabf80e60ee660 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 16 Oct 2025 16:19:41 +0200 Subject: [PATCH 05/29] use ParameterStudy2 in ode_secir examples --- cpp/examples/ode_secir_parameter_study.cpp | 49 ++++++--------- .../ode_secir_parameter_study_graph.cpp | 25 ++++++-- cpp/memilio/compartments/parameter_studies.h | 59 +++++++++++++------ cpp/memilio/mobility/graph_simulation.h | 4 +- .../metapopulation_mobility_instant.h | 2 + cpp/memilio/utils/custom_index_array.h | 13 ++++ cpp/memilio/utils/index.h | 46 +++++++++++++++ 7 files changed, 142 insertions(+), 56 deletions(-) diff --git a/cpp/examples/ode_secir_parameter_study.cpp b/cpp/examples/ode_secir_parameter_study.cpp index 8ad8c1cd2b..751b0b6745 100644 --- a/cpp/examples/ode_secir_parameter_study.cpp +++ b/cpp/examples/ode_secir_parameter_study.cpp @@ -17,6 +17,9 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include "memilio/compartments/simulation.h" +#include "memilio/config.h" +#include "ode_secir/model.h" #include "ode_secir/parameters_io.h" #include "ode_secir/parameter_space.h" #include "memilio/compartments/parameter_studies.h" @@ -30,10 +33,7 @@ * @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)); @@ -46,32 +46,14 @@ 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, + BOOST_OUTCOME_TRY(mio::save_result({sim.get_result()}, {0}, sim.get_model().parameters.get_num_groups().get(), mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5")))); return mio::success(); @@ -83,6 +65,7 @@ int main() ScalarType t0 = 0; ScalarType tmax = 50; + ScalarType dt = 0.1; ScalarType cont_freq = 10; // see Polymod study @@ -141,14 +124,20 @@ int main() //create study auto num_runs = size_t(1); - mio::ParameterStudy> parameter_study(model, t0, tmax, num_runs); + // mio::ParameterStudy2, mio::osecir::Model, ScalarType> + // parameter_study(model, t0, tmax, dt, num_runs); + + auto parameter_study = + mio::make_parameter_study>(model, t0, tmax, dt, num_runs); //run study - auto sample_graph = [](auto&& graph) { - return mio::osecir::draw_sample(graph); + 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(); } diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index b2e8260cc3..c3cb5b7564 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -24,11 +24,13 @@ #include "memilio/io/epi_data.h" #include "memilio/io/result_io.h" #include "memilio/io/mobility_io.h" +#include "memilio/mobility/graph_simulation.h" #include "memilio/mobility/metapopulation_mobility_instant.h" #include "memilio/utils/logging.h" #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" @@ -266,8 +268,16 @@ int main() indices_save_edges); //run parameter study - auto parameter_study = mio::ParameterStudy>{ - params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)}; + // auto parameter_study = mio::ParameterStudy2< + // mio::GraphSimulation>, + // mio::MobilityEdge>, + // ScalarType, ScalarType>, + // mio::Graph, mio::MobilityParameters>, ScalarType>{ + // params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)}; + + auto parameter_study = mio::make_parameter_study_graph_ode>( + params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)); if (mio::mpi::is_root()) { printf("Seeds: "); @@ -279,10 +289,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 +332,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/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index fcbc2ea8db..69ec36ec80 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -412,7 +413,7 @@ class ParameterStudy2 * @param graph_sim_dt time step of graph simulation * @param num_runs number of runs */ - ParameterStudy2(const ParameterType& global_parameters, Time t0, Time tmax, Step dt, size_t num_runs) + ParameterStudy2(const Parameters& global_parameters, Time t0, Time tmax, Step dt, size_t num_runs) : m_parameters(global_parameters) , m_num_runs(num_runs) , m_t0{t0} @@ -482,7 +483,8 @@ class ParameterStudy2 thread_local_rng().set_counter(run_rng_counter); //sample - auto sim = create_simulation(m_parameters, m_t0, m_dt, run_idx); + Simulation 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()); @@ -611,23 +613,42 @@ class ParameterStudy2 RandomNumberGenerator m_rng; }; -// //sample parameters and create simulation -// template -// GraphSimulation create_sampled_simulation(auto&& params_graph, -// SampleGraphFunction sample_graph) -// { -// SimulationGraph sim_graph; - -// auto sampled_graph = sample_graph(params_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); -// } - -// return make_mobility_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); -// } +template +ParameterStudy2 make_parameter_study(const Parameters& global_parameters, Time t0, + Time tmax, Step dt, size_t num_runs) +{ + return {global_parameters, t0, tmax, dt, num_runs}; +} + +template +auto make_parameter_study_graph_ode(const Graph>& global_parameters, FP t0, + FP tmax, FP dt, size_t num_runs) +{ + using SimGraph = Graph, mio::MobilityEdge>; + using SimGraphSim = mio::GraphSimulation; + using Params = Graph>; + + return ParameterStudy2{global_parameters, t0, tmax, dt, num_runs}; +} + +//sample parameters and create simulation +template +GraphSim make_sampled_graph_simulation( + const Graph>& sampled_graph, + FP t0, FP dt_node_sim, FP dt_graph_sim) +{ + typename GraphSim::Graph 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 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 4a01e1392e..f0684557eb 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/custom_index_array.h b/cpp/memilio/utils/custom_index_array.h index d17d496275..e23991db03 100644 --- a/cpp/memilio/utils/custom_index_array.h +++ b/cpp/memilio/utils/custom_index_array.h @@ -23,6 +23,8 @@ #include "memilio/math/eigen_util.h" #include "memilio/utils/index.h" #include "memilio/utils/stl_util.h" +#include +#include namespace { @@ -90,6 +92,17 @@ std::enable_if_t<(I < (Index::size - 1)), std::pair> flatten_ind return {val + (size_t)mio::get(indices) * prod, prod * (size_t)mio::get(dimensions)}; } +template +void assign_from_vector(Index& index, const std::vector& values) +{ + assert(values.size() == sizeof...(Tags)); + + if constexpr (I < sizeof...(Tags)) { + get(index) = values[I]; + assign_from_vector(index, values); + } +} + template struct is_random_access_iterator : std::is_base_of::iterator_category, std::random_access_iterator_tag> { diff --git a/cpp/memilio/utils/index.h b/cpp/memilio/utils/index.h index 595db59f37..54e336c7c0 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 { From 47088f7873405f91261f09086e41ba95a0eaa616 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Tue, 21 Oct 2025 11:30:37 +0200 Subject: [PATCH 06/29] [wip] remove old ParameterStudies from cpp --- cpp/examples/abm_minimal.cpp | 2 +- cpp/examples/ode_secir_read_graph.cpp | 18 +- cpp/memilio/compartments/parameter_studies.h | 360 +------------------ cpp/tests/test_parameter_studies.cpp | 35 +- cpp/tests/test_save_results.cpp | 28 +- 5 files changed, 59 insertions(+), 384 deletions(-) diff --git a/cpp/examples/abm_minimal.cpp b/cpp/examples/abm_minimal.cpp index df6392ec75..4cace851ba 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/cpp/examples/abm_minimal.cpp @@ -237,7 +237,7 @@ int main() 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 " << name_run << std::endl; auto params = std::vector{}; - return std::move(interpolated_result); + return interpolated_result; }); if (ensemble.size() > 0) { diff --git a/cpp/examples/ode_secir_read_graph.cpp b/cpp/examples/ode_secir_read_graph.cpp index 8832aec976..09a88f46f5 100644 --- a/cpp/examples/ode_secir_read_graph.cpp +++ b/cpp/examples/ode_secir_read_graph.cpp @@ -112,8 +112,8 @@ int main(int argc, char** argv) std::cout << "Reading Mobility File..." << std::flush; auto read_mobility_result = mio::read_mobility_plain(filename); 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 +137,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 +146,20 @@ 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); + auto study = mio::make_parameter_study_graph_ode>(graph_read, t0, + tmax, 0.5, 2); + study.run([](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 69ec36ec80..cf266afab3 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -44,351 +44,6 @@ 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. - */ -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 = mio::Graph, EdgeT>; - - /** - * 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 - */ - ParameterStudy(const ParametersGraph& graph, FP t0, FP tmax, FP graph_sim_dt, size_t num_runs) - : m_graph(graph) - , m_num_runs(num_runs) - , m_t0{t0} - , m_tmax{tmax} - , m_dt_graph_sim(graph_sim_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 - */ - 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) - { - for (auto& params_node : m_graph.nodes()) { - set_params_distributions_normal(params_node, t0, tmax, dev_rel); - } - } - - /** - * @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) - { - m_graph.add_node(0, model); - } - - /** - * @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. - */ - template - std::vector> - run(SampleGraphFunction sample_graph, HandleSimulationResultFunction result_processing_function) - { - int num_procs, rank; -#ifdef MEMILIO_ENABLE_MPI - MPI_Comm_size(mpi::get_world(), &num_procs); - MPI_Comm_rank(mpi::get_world(), &rank); -#else - num_procs = 1; - rank = 0; -#endif - - //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 - m_rng.synchronize(); - thread_local_rng() = m_rng; - - auto run_distribution = distribute_runs(m_num_runs, num_procs); - auto 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); - - //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))); - -#ifdef MEMILIO_ENABLE_MPI - //gather results - if (rank == 0) { - for (int src_rank = 1; src_rank < num_procs; ++src_rank) { - int bytes_size; - MPI_Recv(&bytes_size, 1, MPI_INT, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE); - 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{}); - if (!src_ensemble_results) { - log_error("Error receiving ensemble results from rank {}.", src_rank); - } - std::copy(src_ensemble_results.value().begin(), src_ensemble_results.value().end(), - std::back_inserter(ensemble_result)); - } - } - else { - auto bytes = serialize_binary(ensemble_result); - auto 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 - } -#endif - - 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) - { - m_tmax = tmax; - } - - /** - * @brief returns end point in simulation - */ - FP get_tmax() const - { - return m_tmax; - } - - void set_t0(FP t0) - { - m_t0 = t0; - } - - /** - * @brief returns start point in simulation - */ - FP 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 - { - return m_graph.nodes()[0].property; - } - 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. - * @{ - */ - const ParametersGraph& get_model_graph() const - { - return m_graph; - } - ParametersGraph& get_model_graph() - { - return m_graph; - } - /** @} */ - - RandomNumberGenerator& get_rng() - { - return m_rng; - } - -private: - //sample parameters and create simulation - template - GraphSimulation create_sampled_simulation(SampleGraphFunction sample_graph) - { - SimulationGraph sim_graph; - - 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); - } - - return make_mobility_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); - } - - std::vector distribute_runs(size_t num_runs, int num_procs) - { - //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; - - std::vector run_distribution(num_procs); - std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); - std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local); - - 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; -}; - /** * Class that performs multiple simulation runs with randomly sampled parameters. * Can simulate mobility graphs with one simulation in each node or single simulations. @@ -442,7 +97,7 @@ class ParameterStudy2 * @param process_simulation_result A function that accepts S */ template - std::vector> + std::vector>> run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) { static_assert(std::is_invocable_r_v, @@ -469,7 +124,8 @@ class ParameterStudy2 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; + 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++) { @@ -536,10 +192,7 @@ class ParameterStudy2 template std::vector run(CreateSimulationFunction&& create_simulation) { - return run(std::forward(create_simulation), - [](Simulation&& sim, size_t) -> Simulation&& { - return sim; - }); + return run(std::forward(create_simulation), &result_forwarding_function); } /** @@ -597,6 +250,11 @@ class ParameterStudy2 return run_distribution; } + inline static Simulation&& result_forwarding_function(Simulation&& sim, size_t) + { + return std::move(sim); + } + private: // Stores Graph with the names and ranges of all parameters ParameterType m_parameters; diff --git a/cpp/tests/test_parameter_studies.cpp b/cpp/tests/test_parameter_studies.cpp index e0f46fbc5a..b7f63a6852 100644 --- a/cpp/tests/test_parameter_studies.cpp +++ b/cpp/tests/test_parameter_studies.cpp @@ -157,14 +157,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); + auto study = mio::make_parameter_study_graph_ode>(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([](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 +373,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); + auto parameter_study = mio::make_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([](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++) { diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index a15c62d419..0fd07553aa 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -169,8 +169,9 @@ 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; + auto parameter_study = + mio::make_parameter_study_graph_ode>(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 +184,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 +306,9 @@ 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; + auto parameter_study = + mio::make_parameter_study_graph_ode>(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 +322,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(); From 26460422b623c0bc9a609a7d598c1d41dc9e9181 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Tue, 21 Oct 2025 14:41:26 +0200 Subject: [PATCH 07/29] Fix up ode_secir_read_graph --- cpp/examples/ode_secir_read_graph.cpp | 51 +++++++++++---------------- cpp/memilio/io/cli.h | 3 ++ cpp/memilio/utils/base_dir.h | 8 +++-- cpp/memilio/utils/index.h | 16 ++++++--- cpp/memilio/utils/type_list.h | 6 ++++ 5 files changed, 47 insertions(+), 37 deletions(-) diff --git a/cpp/examples/ode_secir_read_graph.cpp b/cpp/examples/ode_secir_read_graph.cpp index 09a88f46f5..e2b9df4bd7 100644 --- a/cpp/examples/ode_secir_read_graph.cpp +++ b/cpp/examples/ode_secir_read_graph.cpp @@ -17,43 +17,34 @@ * 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/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/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,7 +101,7 @@ 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 << "\n" << read_mobility_result.error().formatted_message() << "\n"; std::cout << "Create the mobility file with MEmilio Epidata's getCommuterMobility.py file.\n"; 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/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 54e336c7c0..ac76fe3b38 100644 --- a/cpp/memilio/utils/index.h +++ b/cpp/memilio/utils/index.h @@ -227,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 { }; -template -struct index_of_type> : public index_of_type { +/// Specialization of index_of_type for Index. @see 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/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_ From 9b593f51f2b373a1b91acb42caf4ef3555613f20 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Wed, 22 Oct 2025 11:27:50 +0200 Subject: [PATCH 08/29] patch python bindings to mostly emulate old ParameterStudy behaviour --- cpp/memilio/compartments/parameter_studies.h | 27 ++-- .../bindings/compartments/parameter_studies.h | 103 +++++++++++++++ .../simulation/bindings/models/osecir.cpp | 121 +++++++++--------- .../simulation/bindings/models/osecirvvs.cpp | 107 +++++++++------- .../simulation_test/test_parameter_study.py | 6 +- 5 files changed, 240 insertions(+), 124 deletions(-) create mode 100644 pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index cf266afab3..73a5d596a9 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -23,23 +23,14 @@ #include "memilio/io/binary_serializer.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 { @@ -218,6 +209,11 @@ class ParameterStudy2 { return m_t0; } + + Time get_dt() const + { + return m_dt; + } /** * Get the input graph that the parameter study is run for. * Use for graph simulations, use get_model for single node simulations. @@ -227,6 +223,10 @@ class ParameterStudy2 { return m_parameters; } + Parameters& get_parameters() + { + return m_parameters; + } /** @} */ RandomNumberGenerator& get_rng() @@ -255,7 +255,6 @@ class ParameterStudy2 return std::move(sim); } -private: // Stores Graph with the names and ranges of all parameters ParameterType m_parameters; @@ -271,7 +270,7 @@ class ParameterStudy2 RandomNumberGenerator m_rng; }; -template +template ParameterStudy2 make_parameter_study(const Parameters& global_parameters, Time t0, Time tmax, Step dt, size_t num_runs) { @@ -282,11 +281,11 @@ template auto make_parameter_study_graph_ode(const Graph>& global_parameters, FP t0, FP tmax, FP dt, size_t num_runs) { - using SimGraph = Graph, mio::MobilityEdge>; - using SimGraphSim = mio::GraphSimulation; + using SimGraph = Graph, MobilityEdge>; + using SimGraphSim = GraphSimulation; using Params = Graph>; - return ParameterStudy2{global_parameters, t0, tmax, dt, num_runs}; + return ParameterStudy2{global_parameters, t0, tmax, dt, num_runs}; } //sample parameters and create simulation diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h b/pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h new file mode 100644 index 0000000000..8cb60ffe12 --- /dev/null +++ b/pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h @@ -0,0 +1,103 @@ +/* +* 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 "memilio/compartments/parameter_studies.h" +#include "memilio/compartments/compartmental_model.h" + +#include "pybind11/pybind11.h" +#include "pybind11/stl_bind.h" + +#include +#include +#include +#include +#include + +namespace py = pybind11; + +namespace pymio +{ + +/* + * @brief bind ParameterStudy for any model + */ +template +void bind_GrapParameterStudy( + py::module_& m, std::string const& name, std::vector argnames, + std::function create_simulation = + [](const mio::Graph>& g, double t0, double dt, + size_t) { + using GraphSim = mio::GraphSimulation< + double, + mio::Graph>, mio::MobilityEdge>, + double, double>; + auto copy = g; + return mio::make_sampled_graph_simulation(draw_sample(copy), t0, dt, dt) + }) +{ + assert(sizeof...(RunArgs) == argnames.size()); + using SimulationT = mio::GraphSimulation< + double, mio::Graph>, mio::MobilityEdge>, + double, double>; + using ParametersT = Graph < typename Sim::Model, MobilityParameters; + using StudyT = mio::ParameterStudy2; + using CreateSimulationFunctionT = std::function; + using ProcessSimulationResultFunctionT = std::function; + 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("parameters", py::overload_cast<>(&StudyT::get_parameters), + py::return_value_policy::reference_internal) + .def_property_readonly("rng", &StudyT::get_rng, py::return_value_policy::reference_internal) + .def( + "run", + [](StudyT& self, const ProcessSimulationResultFunctionT& handle_result, RunArgs args...) { + self.run(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), run_idx); + return 0; + }); + }, + py::arg("handle_result_func")) + .def("run", + [](StudyT& self) { //default argument doesn't seem to work with functions + return self.run(create_simulation); + }) + .def( + "run_single", + [](StudyT& self, ProcessSimulationResultFunctionT handle_result) { + self.run(create_simulation, [&handle_result](auto&& r, auto&& run_idx) { + handle_result(std::move(r.nodes()[0].property.get_simulation()), run_idx); + return 0; + }); + }, + py::arg("handle_result_func")) + .def("run_single", [](StudyT& self) { + return filter_graph_results(self.run(create_simulation)); + }); +} + +} // namespace pymio diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index 273da48f2d..5fc722ed8d 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp @@ -48,6 +48,7 @@ #include "pybind11/pybind11.h" #include "pybind11/stl_bind.h" #include "Eigen/Core" +#include #include namespace py = pybind11; @@ -58,88 +59,84 @@ 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, mio::MobilityEdge>, + double, double>>&& 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())); + results.emplace_back(std::move(graph_results[i].get_graph().nodes()[0].property.get_simulation())); } return std::move(results); } +/// Moves out graphs from GraphSimulation%s. Helps make the new ParameterStudy backwards compatible with older bindings. +template +std::vector extract_graph_from_graph_simulation(std::vector&& result_graph_sims) +{ + std::vector result_graphs; + result_graphs.reserve(result_graph_sims.size()); + for (const SimulationT& sim : result_graph_sims) { + result_graphs.emplace_back(std::move(sim.get_graph())); + } + return result_graphs; +} + /* * @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::ParameterStudy2; + + 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(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 extract_graph_from_graph_simulation(self.run(create_simulation)); }) .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(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 filter_graph_results(self.run(create_simulation)); }); } @@ -287,11 +284,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..7bb1210402 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -52,95 +52,110 @@ 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, mio::MobilityEdge>, + double, double>>&& 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())); + results.emplace_back(std::move(graph_results[i].get_graph().nodes()[0].property.get_simulation())); } return std::move(results); } +/// Moves out graphs from GraphSimulation%s. Helps make the new ParameterStudy backwards compatible with older bindings. +template +std::vector extract_graph_from_graph_simulation(std::vector&& result_graph_sims) +{ + std::vector result_graphs; + result_graphs.reserve(result_graph_sims.size()); + for (const SimulationT& sim : result_graph_sims) { + result_graphs.emplace_back(std::move(sim.get_graph())); + } + return result_graphs; +} + /* * @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::ParameterStudy2; + + 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); + [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 extract_graph_from_graph_simulation( + self.run([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); + })); }, py::arg("variant_high")) .def( "run_single", - [](mio::ParameterStudy& self, std::function handle_result, - bool variant_high) { + [](StudyT& self, std::function handle_result, bool variant_high) { self.run( - [variant_high](auto&& g) { - return draw_sample(g, variant_high); + [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 filter_graph_results( + self.run([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); + })); }, py::arg("variant_high")); } @@ -340,9 +355,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) From df4f000e6ac3151e1f55ad5f05ab91079ec62fe2 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Wed, 22 Oct 2025 11:46:10 +0200 Subject: [PATCH 09/29] explicit type conversion --- cpp/examples/abm_minimal.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/examples/abm_minimal.cpp b/cpp/examples/abm_minimal.cpp index 4cace851ba..23e8e9ad51 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/cpp/examples/abm_minimal.cpp @@ -33,6 +33,7 @@ #include #include +#include #include constexpr size_t num_age_groups = 4; @@ -227,7 +228,7 @@ int main() [](auto&&, auto t0_, auto, size_t run_idx) { // TODO: change the parameters around a bit? - auto sim = ResultSim(make_model(run_idx), t0_); + auto sim = ResultSim(make_model((uint32_t)run_idx), t0_); return sim; }, [](auto&& sim, auto&& run_idx) { From 2d278c0e389c2849f55b5e16b0fa79dd74bed14e Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Wed, 22 Oct 2025 12:22:20 +0200 Subject: [PATCH 10/29] another explicit type conversion --- cpp/examples/ode_secir_parameter_study.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/examples/ode_secir_parameter_study.cpp b/cpp/examples/ode_secir_parameter_study.cpp index 751b0b6745..60981be83e 100644 --- a/cpp/examples/ode_secir_parameter_study.cpp +++ b/cpp/examples/ode_secir_parameter_study.cpp @@ -53,7 +53,7 @@ mio::IOResult write_single_run_result(const size_t run, const mio::osecir: std::vector> all_results; std::vector ids; - BOOST_OUTCOME_TRY(mio::save_result({sim.get_result()}, {0}, sim.get_model().parameters.get_num_groups().get(), + 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(); From 5fa14f43451ebf81f482d5ccf6ef263e822e133c Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Wed, 22 Oct 2025 14:15:45 +0200 Subject: [PATCH 11/29] improve coverage --- cpp/tests/test_io_cli.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) 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\"."); From c6e292aaab4658f226be168d1be555bd3c9b79a4 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 23 Oct 2025 16:24:04 +0200 Subject: [PATCH 12/29] move abm study, restore minimal abm --- cpp/examples/CMakeLists.txt | 4 + cpp/examples/abm_minimal.cpp | 165 +++---------------- cpp/examples/abm_parameter_study.cpp | 230 +++++++++++++++++++++++++++ 3 files changed, 257 insertions(+), 142 deletions(-) create mode 100644 cpp/examples/abm_parameter_study.cpp diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 39a12539e1..d303ac241e 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -108,6 +108,10 @@ 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}) +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}) + 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 23e8e9ad51..643a17b1df 100644 --- a/cpp/examples/abm_minimal.cpp +++ b/cpp/examples/abm_minimal.cpp @@ -17,76 +17,26 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "abm/simulation.h" -#include "abm/common_abm_loggers.h" -#include #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/result_io.h" -#include "memilio/utils/abstract_parameter_distribution.h" -#include "memilio/utils/miompi.h" -#include "memilio/utils/time_series.h" +#include "abm/common_abm_loggers.h" -#include -#include -#include #include -constexpr size_t num_age_groups = 4; - -// shim class to make the abm simulation look like an ode simulation -class ResultSim : public mio::abm::Simulation -{ -public: - using Model = mio::abm::Model; - - ResultSim(const Model& m, double t, double) - : mio::abm::Simulation(mio::abm::TimePoint{mio::abm::days(t).seconds()}, Model(m)) - { - history.log(*this); // set initial results - } - - ResultSim(const Model& m, mio::abm::TimePoint t) - : mio::abm::Simulation(t, Model(m)) - { - history.log(*this); // set initial results - } - - void advance(mio::abm::TimePoint tmax) - { - mio::abm::Simulation::advance(tmax, history); - } - - const mio::TimeSeries& get_result() const - { - return get<0>(history.get_log()); - } - - mio::History history{ - Eigen::Index(mio::abm::InfectionState::Count)}; -}; - -mio::abm::Model make_model(uint32_t run_idx) +int main() { - + // This is a minimal example with children and adults < 60 year old. + // We divided them into 4 different age groups, which are defined as follows: + mio::set_log_level(mio::LogLevel::warn); + size_t num_age_groups = 4; 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); - std::initializer_list seeds = {14159265, 243141u + run_idx}; - auto rng = mio::RandomNumberGenerator(); - rng.seed(seeds); - auto run_rng_counter = - mio::rng_totalsequence_counter(static_cast(run_idx), mio::Counter(0)); - rng.set_counter(run_rng_counter); - model.get_rng() = rng; + // Create the model with 4 age groups. + auto model = mio::abm::Model(num_age_groups); // 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.); @@ -100,7 +50,7 @@ mio::abm::Model make_model(uint32_t run_idx) model.parameters.check_constraints(); // There are 10 households for each household group. - int n_households = 100; + 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. @@ -171,9 +121,9 @@ mio::abm::Model make_model(uint32_t run_idx) 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); + auto 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(), + person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(), model.parameters, start_date, infection_state)); } } @@ -200,95 +150,26 @@ mio::abm::Model make_model(uint32_t run_idx) 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(); - // This is a minimal example with children and adults < 60 year old. - // We divided them into 4 different age groups, which are defined as follows: - - 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 = 10; + auto tmax = t0 + mio::abm::days(10); + auto sim = mio::abm::Simulation(t0, std::move(model)); // Create a history object to store the time series of the infection states. mio::History historyTimeSeries{ Eigen::Index(mio::abm::InfectionState::Count)}; - mio::ParameterStudy2 study( - 1, t0, tmax, mio::abm::TimeSpan(0), num_runs); - - auto ensemble = study.run( - [](auto&&, auto t0_, auto, size_t run_idx) { - // TODO: change the parameters around a bit? - - auto sim = ResultSim(make_model((uint32_t)run_idx), t0_); - return sim; - }, - [](auto&& sim, auto&& run_idx) { - auto interpolated_result = mio::interpolate_simulation_result(sim.get_result()); - std::string name_run = "abm_minimal_run_" + std::to_string(run_idx) + ".txt"; - std::ofstream outfile_run(name_run); - 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 " << name_run << std::endl; - auto params = std::vector{}; - return interpolated_result; - }); - - if (ensemble.size() > 0) { - auto ensemble_results = std::vector>>{}; - ensemble_results.reserve(ensemble.size()); - auto ensemble_params = std::vector>{}; - ensemble_params.reserve(ensemble.size()); - // int j = 0; - // std::cout << "## " << ensemble[0].first.get_value(0).transpose() << "\n"; - // for (auto&& run : ensemble) { - // std::cout << j++ << " " << run.first.get_last_value().transpose() << "\n"; - // ensemble_results.emplace_back(std::vector{std::move(run.first)}); - // ensemble_params.emplace_back(std::vector{std::move(run.second)}); - // } - - // 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); - - // const boost::filesystem::path save_dir = "/Users/saschakorf/Nosynch/Arbeit/memilio/cpp/examples/results"; - - // if (!save_dir.empty()) { - - // mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups, - // (save_dir / ("Results_" + std::string("p05") + ".h5")).string())); - // mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups, - // (save_dir / ("Results_" + std::string("p25") + ".h5")).string())); - // mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups, - // (save_dir / ("Results_" + std::string("p50") + ".h5")).string())); - // mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups, - // (save_dir / ("Results_" + std::string("p75") + ".h5")).string())); - // mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups, - // (save_dir / ("Results_" + std::string("p95") + ".h5")).string())); - // } - } // Run the simulation until tmax with the history object. - // sim.advance(tmax, historyTimeSeries); - - // // The results are written into the file "abm_minimal.txt" as a table with 9 columns. - // // The first column is Time. The other columns correspond to the number of people with a certain infection state at this Time: - // // Time = Time in days, S = Susceptible, E = Exposed, I_NS = InfectedNoSymptoms, I_Sy = InfectedSymptoms, I_Sev = InfectedSevere, - // // I_Crit = InfectedCritical, R = Recovered, D = Dead - // std::ofstream outfile("abm_minimal.txt"); - // std::get<0>(historyTimeSeries.get_log()) - // .print_table(outfile, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4); - // std::cout << "Results written to abm_minimal.txt" << std::endl; - - mio::mpi::finalize(); + sim.advance(tmax, historyTimeSeries); + + // The results are written into the file "abm_minimal.txt" as a table with 9 columns. + // The first column is Time. The other columns correspond to the number of people with a certain infection state at this Time: + // Time = Time in days, S = Susceptible, E = Exposed, I_NS = InfectedNoSymptoms, I_Sy = InfectedSymptoms, I_Sev = InfectedSevere, + // I_Crit = InfectedCritical, R = Recovered, D = Dead + std::ofstream outfile("abm_minimal.txt"); + std::get<0>(historyTimeSeries.get_log()) + .print_table(outfile, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4); + std::cout << "Results written to abm_minimal.txt" << std::endl; return 0; } diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp new file mode 100644 index 0000000000..326cbdc12a --- /dev/null +++ b/cpp/examples/abm_parameter_study.cpp @@ -0,0 +1,230 @@ +/* +* 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 = 100; + + // 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 = 10; + + // Create a parameter study. We currently do not use Parameters or dt, so we set them to 0. + mio::ParameterStudy2, int, mio::abm::TimePoint, mio::abm::TimeSpan> + 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) { + auto sim = mio::abm::ResultSimulation(make_model(mio::thread_local_rng()), t0_); + return sim; + }, + [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; +} From 3af580cc798784582b9f163f99e607b1b6e75647 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 23 Oct 2025 18:52:35 +0200 Subject: [PATCH 13/29] clean up examples, improve mpi handling --- cpp/examples/abm_parameter_study.cpp | 4 +- cpp/examples/ode_secir_parameter_study.cpp | 26 ++-- cpp/examples/ode_secir_read_graph.cpp | 2 +- cpp/memilio/compartments/parameter_studies.h | 154 +++++++++++++------ cpp/memilio/io/io.h | 19 +++ cpp/memilio/utils/miompi.cpp | 3 + cpp/models/abm/CMakeLists.txt | 1 + cpp/models/abm/result_simulation.h | 65 ++++++++ 8 files changed, 212 insertions(+), 62 deletions(-) create mode 100644 cpp/models/abm/result_simulation.h diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 326cbdc12a..183c830c28 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -62,7 +62,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) model.parameters.check_constraints(); // There are 10 households for each household group. - int n_households = 100; + 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. @@ -175,7 +175,7 @@ int main() 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 = 10; + const size_t num_runs = 3; // Create a parameter study. We currently do not use Parameters or dt, so we set them to 0. mio::ParameterStudy2, int, mio::abm::TimePoint, mio::abm::TimeSpan> diff --git a/cpp/examples/ode_secir_parameter_study.cpp b/cpp/examples/ode_secir_parameter_study.cpp index 60981be83e..f8f05a2fc8 100644 --- a/cpp/examples/ode_secir_parameter_study.cpp +++ b/cpp/examples/ode_secir_parameter_study.cpp @@ -17,13 +17,14 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "memilio/compartments/simulation.h" #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" /** @@ -35,8 +36,8 @@ */ 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'; @@ -54,14 +55,15 @@ mio::IOResult write_single_run_result(const size_t run, const mio::osecir: std::vector ids; 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")))); + 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; @@ -123,7 +125,7 @@ int main() } //create study - auto num_runs = size_t(1); + auto num_runs = size_t(3); // mio::ParameterStudy2, mio::osecir::Model, ScalarType> // parameter_study(model, t0, tmax, dt, num_runs); @@ -141,9 +143,15 @@ int main() 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 sim.get_result(); //Result handler must return something, but only meaningful when using MPI. }; - 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}); + + auto result = parameter_study.run(sample_graph, handle_result); + + mio::mpi::finalize(); return 0; } diff --git a/cpp/examples/ode_secir_read_graph.cpp b/cpp/examples/ode_secir_read_graph.cpp index e2b9df4bd7..86fe53a405 100644 --- a/cpp/examples/ode_secir_read_graph.cpp +++ b/cpp/examples/ode_secir_read_graph.cpp @@ -146,7 +146,7 @@ int main(int argc, char** argv) std::cout << "Running Simulations..." << std::flush; auto study = mio::make_parameter_study_graph_ode>(graph_read, t0, tmax, 0.5, 2); - study.run([](auto&& g, auto t0_, auto dt_, auto) { + 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_); diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 73a5d596a9..69a89485b4 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -21,8 +21,10 @@ #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/mobility/metapopulation_mobility_instant.h" @@ -82,41 +84,24 @@ class ParameterStudy2 // // TODO how is this supposed to work wrt. model? is this just a special case where ParameterType=Model? // } - /** - * @brief - * @param sample_simulation A function that accepts ParameterType and returns an instance of SimulationType. - * @param process_simulation_result A function that accepts S - */ - template +private: + template std::vector>> - run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) + run_impl(size_t start_run_idx, size_t end_run_idx, CreateSimulationFunction&& create_simulation, + ProcessSimulationResultFunction&& process_simulation_result) { + assert(start_run_idx <= end_run_idx); static_assert(std::is_invocable_r_v, "Incorrect Type for create_simulation."); static_assert(std::is_invocable_v, "Incorrect Type for process_simulation_result."); - int num_procs, rank; -#ifdef MEMILIO_ENABLE_MPI - MPI_Comm_size(mpi::get_world(), &num_procs); - MPI_Comm_rank(mpi::get_world(), &rank); -#else - num_procs = 1; - rank = 0; -#endif - //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 - m_rng.synchronize(); - thread_local_rng() = m_rng; + using ProcessedResultT = + std::decay_t>; - auto run_distribution = distribute_runs(m_num_runs, num_procs); - auto 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)]; + thread_local_rng() = m_rng; - std::vector>> - ensemble_result; + 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++) { @@ -145,6 +130,76 @@ class ParameterStudy2 //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; + } + +public: + template + std::vector>> + run_serial(CreateSimulationFunction&& create_simulation, + ProcessSimulationResultFunction&& process_simulation_result) + { + return run_impl(0, m_num_runs, std::forward(create_simulation), + std::forward(process_simulation_result)); + } + + template + std::vector run_serial(CreateSimulationFunction&& create_simulation) + { + return run_serial(std::forward(create_simulation), + [](Simulation&& sim, size_t) -> Simulation&& { + return std::move(sim); + }); + } + + /** + * @brief run the sim + * @param sample_simulation A function that accepts ParameterType and returns an instance of SimulationType. + * @param process_simulation_result A function that accepts S + * side effect: sets seed and count of thread_local_rng + */ + template + std::vector>> + run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) + { + static_assert(std::is_invocable_r_v, + "Incorrect Type for create_simulation."); + static_assert(std::is_invocable_v, + "Incorrect Type for process_simulation_result."); + + using ProcessedResultT = + std::decay_t>; + int num_procs, rank; + +#ifdef MEMILIO_ENABLE_MPI + MPI_Comm_size(mpi::get_world(), &num_procs); + MPI_Comm_rank(mpi::get_world(), &rank); + + static_assert(is_serializable && + is_deserializable, + "The result type of the processing function must be serializable for usage with MPI. " + "Change the process_simulation_result argument of ParameterStudy::run, or use run_serial."); +#else + num_procs = 1; + rank = 0; +#endif + + //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 + m_rng.synchronize(); + // Note that this overwrites seed and counter of thread_local_rng, but it does not replace it. + thread_local_rng() = m_rng; + + auto run_distribution = distribute_runs(m_num_runs, num_procs); + auto 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 = + run_impl(start_run_idx, end_run_idx, std::forward(create_simulation), + std::forward(process_simulation_result)); + #ifdef MEMILIO_ENABLE_MPI //gather results if (rank == 0) { @@ -183,11 +238,14 @@ class ParameterStudy2 template std::vector run(CreateSimulationFunction&& create_simulation) { - return run(std::forward(create_simulation), &result_forwarding_function); + return run(std::forward(create_simulation), + [](Simulation&& sim, size_t) -> Simulation&& { + return std::move(sim); + }); } /** - * @brief returns the number of Monte Carlo runs + * @brief Return the number of total runs that the study will make. */ size_t get_num_runs() const { @@ -195,7 +253,7 @@ class ParameterStudy2 } /** - * @brief returns end point in simulation + * @brief Return the final time point for simulations. */ Time get_tmax() const { @@ -203,20 +261,23 @@ class ParameterStudy2 } /** - * @brief returns start point in simulation + * @brief Return the initial time point for simulations. */ Time get_t0() const { return m_t0; } + /** + * @brief Return the initial step sized used by simulations. + */ Time get_dt() const { return m_dt; } + /** - * 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 Parameters& get_parameters() const @@ -229,12 +290,18 @@ class ParameterStudy2 } /** @} */ + /** + * @brief Access the study's random number generator. + */ RandomNumberGenerator& get_rng() { return m_rng; } private: + /** + * @brief Create a vector with the number of runs each process should make. + */ std::vector distribute_runs(size_t num_runs, int num_procs) { assert(num_procs > 0); @@ -250,24 +317,11 @@ class ParameterStudy2 return run_distribution; } - inline static Simulation&& result_forwarding_function(Simulation&& sim, size_t) - { - return std::move(sim); - } - - // Stores Graph with the names and ranges of all parameters - ParameterType m_parameters; - - size_t m_num_runs; - - // Start time (should be the same for all simulations) - Time m_t0; - // End time (should be the same for all simulations) - Time m_tmax; - // adaptive time step of the integrator (will be corrected if too large/small) - Step m_dt; - // - 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. }; template diff --git a/cpp/memilio/io/io.h b/cpp/memilio/io/io.h index fdcf9424d4..13c2d30045 100644 --- a/cpp/memilio/io/io.h +++ b/cpp/memilio/io/io.h @@ -857,6 +857,15 @@ void serialize(IOContext& io, const T& t) serialize_internal(io, t); } +namespace details +{ +template +using is_serializable_expr_t = decltype(serialize(std::declval(), std::declval())); +} + +template +static constexpr bool is_serializable = is_expression_valid::value; + /** * Restores an object from the data stored in an IO context. * There must be provided for the type T either a free function `deserialize_internal(io, tag)` @@ -881,6 +890,16 @@ IOResult deserialize(IOContext& io, Tag tag) return deserialize_internal(io, tag); } +namespace details +{ +template +using is_deserializable_expr_t = + decltype(std::declval>() = deserialize(std::declval(), std::declval>())); +} + +template +static constexpr bool is_deserializable = is_expression_valid::value; + /** * @brief Returns the current working directory name */ 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/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/result_simulation.h b/cpp/models/abm/result_simulation.h new file mode 100644 index 0000000000..01e4b6713b --- /dev/null +++ b/cpp/models/abm/result_simulation.h @@ -0,0 +1,65 @@ +/* +* 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, copying the model. + ResultSimulation(const Model& m, TimePoint t) + : Simulation(t, Model(m)) + { + history.log(*this); // set initial results + } + + /** + * @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 From 437162200ce03de3f9b3f891770b569afd8441a2 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 23 Oct 2025 19:05:59 +0200 Subject: [PATCH 14/29] add hdf5 dep to abm study example --- cpp/examples/CMakeLists.txt | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index d303ac241e..fa12445425 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -108,9 +108,11 @@ 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}) -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}) +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) From fe9ef05ddbc5a72cd832f0bd0af174105cfaba14 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 24 Oct 2025 13:14:53 +0200 Subject: [PATCH 15/29] improve non-PS coverage --- cpp/tests/test_abm_simulation.cpp | 29 +++++++++++++++++++++++++++++ cpp/tests/test_utils.cpp | 13 +++++++++++++ 2 files changed, 42 insertions(+) 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_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")); +} From fcde395ae246b1331d270250349cfe2b48984bd9 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 24 Oct 2025 14:49:25 +0200 Subject: [PATCH 16/29] remove Simulation template from ParameterStudy2. this allows deduction of template arguments from the constructor, making the make_parameter_study functions obsolete. the simulation type is now deduced from the create_simulation argument, and checks on invocability have been moved to clean up deduced types. Add documentation. --- cpp/examples/abm_parameter_study.cpp | 8 +- cpp/examples/ode_secir_parameter_study.cpp | 14 +- .../ode_secir_parameter_study_graph.cpp | 19 +- cpp/examples/ode_secir_read_graph.cpp | 9 +- cpp/memilio/compartments/parameter_studies.h | 294 +++++++++--------- cpp/models/abm/result_simulation.h | 7 +- cpp/tests/test_parameter_studies.cpp | 17 +- cpp/tests/test_save_results.cpp | 14 +- .../bindings/compartments/parameter_studies.h | 103 ------ .../simulation/bindings/models/osecir.cpp | 44 +-- .../simulation/bindings/models/osecirvvs.cpp | 67 ++-- 11 files changed, 223 insertions(+), 373 deletions(-) delete mode 100644 pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 183c830c28..f0a2558d6f 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -177,9 +177,8 @@ int main() // auto sim = mio::abm::Simulation(t0, std::move(model)); const size_t num_runs = 3; - // Create a parameter study. We currently do not use Parameters or dt, so we set them to 0. - mio::ParameterStudy2, int, mio::abm::TimePoint, mio::abm::TimeSpan> - study(0, t0, tmax, mio::abm::TimeSpan(0), num_runs); + // Create a parameter study. The ABM currently does not use parameters or dt, so we set them both to 0. + mio::ParameterStudy2 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}); @@ -192,8 +191,7 @@ int main() auto ensemble_results = study.run( [](auto, auto t0_, auto, size_t) { - auto sim = mio::abm::ResultSimulation(make_model(mio::thread_local_rng()), t0_); - return sim; + 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()); diff --git a/cpp/examples/ode_secir_parameter_study.cpp b/cpp/examples/ode_secir_parameter_study.cpp index f8f05a2fc8..0d2ff9c7d7 100644 --- a/cpp/examples/ode_secir_parameter_study.cpp +++ b/cpp/examples/ode_secir_parameter_study.cpp @@ -69,6 +69,7 @@ int main() 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, @@ -124,15 +125,11 @@ int main() return -1; } - //create study + // create study auto num_runs = size_t(3); - // mio::ParameterStudy2, mio::osecir::Model, ScalarType> - // parameter_study(model, t0, tmax, dt, num_runs); + mio::ParameterStudy2 parameter_study(model, t0, tmax, dt, num_runs); - auto parameter_study = - mio::make_parameter_study>(model, t0, tmax, dt, num_runs); - - //run study + // 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); @@ -143,12 +140,13 @@ int main() if (!write_result_status) { std::cout << "Error writing result: " << write_result_status.error().formatted_message(); } - return sim.get_result(); //Result handler must return something, but only meaningful when using MPI. + return 0; // Result handler must return something. }; // 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(); diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index c3cb5b7564..b0f8202f0b 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -20,11 +20,9 @@ #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" -#include "memilio/mobility/graph_simulation.h" #include "memilio/mobility/metapopulation_mobility_instant.h" #include "memilio/utils/logging.h" #include "memilio/utils/miompi.h" @@ -33,11 +31,8 @@ #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. @@ -267,17 +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::ParameterStudy2< - // mio::GraphSimulation>, - // mio::MobilityEdge>, - // ScalarType, ScalarType>, - // mio::Graph, mio::MobilityParameters>, ScalarType>{ - // params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)}; - - auto parameter_study = mio::make_parameter_study_graph_ode>( - params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)); + mio::ParameterStudy2 parameter_study(params_graph, 0.0, num_days_sim, 0.5, size_t(num_runs)); if (mio::mpi::is_root()) { printf("Seeds: "); @@ -291,7 +276,7 @@ int main() auto ensemble = parameter_study.run( [](auto&& graph, ScalarType t0, ScalarType dt, size_t) { auto copy = graph; - return mio::make_sampled_graph_simulation( + return mio::make_sampled_graph_simulation>( mio::osecir::draw_sample(copy), t0, dt, dt); }, [&](auto&& results_sim, auto&& run_id) { diff --git a/cpp/examples/ode_secir_read_graph.cpp b/cpp/examples/ode_secir_read_graph.cpp index 86fe53a405..e282ef5468 100644 --- a/cpp/examples/ode_secir_read_graph.cpp +++ b/cpp/examples/ode_secir_read_graph.cpp @@ -18,12 +18,14 @@ * limitations under the License. */ #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" @@ -144,12 +146,11 @@ int main(int argc, char** argv) auto& graph_read = graph_read_result.value(); std::cout << "Running Simulations..." << std::flush; - auto study = mio::make_parameter_study_graph_ode>(graph_read, t0, - tmax, 0.5, 2); + mio::ParameterStudy2 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_); + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); }); std::cout << "Done" << std::endl; diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 69a89485b4..f3eb730b66 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -31,38 +31,53 @@ #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 ParameterStudy2 { public: - using Simulation = SimulationType; using Parameters = ParameterType; using Time = TimeType; using Step = StepType; - // TODO: replacement for "set_params_distributions_normal". Maybe a special ctor for UncertainParameterSet? +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. */ - ParameterStudy2(const Parameters& global_parameters, Time t0, Time tmax, Step dt, size_t num_runs) - : m_parameters(global_parameters) + ParameterStudy2(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} @@ -70,72 +85,25 @@ class ParameterStudy2 { } - // /** - // * @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 - // */ - // template >> - // ParameterStudy2(typename Simulation::Model const& model, Time t0, Time tmax, Step dt, size_t num_runs) - // : ParameterStudy2(model, t0, tmax, dt, num_runs) - // { - // // TODO how is this supposed to work wrt. model? is this just a special case where ParameterType=Model? - // } - -private: - template - std::vector>> - run_impl(size_t start_run_idx, size_t end_run_idx, CreateSimulationFunction&& create_simulation, - ProcessSimulationResultFunction&& process_simulation_result) - { - assert(start_run_idx <= end_run_idx); - static_assert(std::is_invocable_r_v, - "Incorrect Type for create_simulation."); - static_assert(std::is_invocable_v, - "Incorrect Type for process_simulation_result."); - - using ProcessedResultT = - std::decay_t>; - - 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 - 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 - Simulation 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)); - } - - //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; - } - -public: + /** + * @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. + * @{ + */ template - std::vector>> + std::vector> run_serial(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) { @@ -144,41 +112,45 @@ class ParameterStudy2 } template - std::vector run_serial(CreateSimulationFunction&& create_simulation) + std::vector> run_serial(CreateSimulationFunction&& create_simulation) { - return run_serial(std::forward(create_simulation), - [](Simulation&& sim, size_t) -> Simulation&& { - return std::move(sim); - }); + return run_serial( + std::forward(create_simulation), + [](SimulationT&& sim, size_t) -> SimulationT&& { + return std::move(sim); + }); } + /** @} */ /** - * @brief run the sim - * @param sample_simulation A function that accepts ParameterType and returns an instance of SimulationType. - * @param process_simulation_result A function that accepts S - * side effect: sets seed and count of thread_local_rng + * @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>> + std::vector> run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result) { - static_assert(std::is_invocable_r_v, - "Incorrect Type for create_simulation."); - static_assert(std::is_invocable_v, - "Incorrect Type for process_simulation_result."); - - using ProcessedResultT = - std::decay_t>; + 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); - - static_assert(is_serializable && - is_deserializable, - "The result type of the processing function must be serializable for usage with MPI. " - "Change the process_simulation_result argument of ParameterStudy::run, or use run_serial."); #else num_procs = 1; rank = 0; @@ -188,15 +160,13 @@ class ParameterStudy2 //So we set our own RNG to be used. //Assume that sampling uses the thread_local_rng() and isn't multithreaded m_rng.synchronize(); - // Note that this overwrites seed and counter of thread_local_rng, but it does not replace it. - 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)]; + size_t end_run_idx = start_run_idx + run_distribution[size_t(rank)]; - std::vector ensemble_result = + EnsembleResultT ensemble_result = run_impl(start_run_idx, end_run_idx, std::forward(create_simulation), std::forward(process_simulation_result)); @@ -209,7 +179,7 @@ class ParameterStudy2 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); } @@ -218,8 +188,8 @@ class ParameterStudy2 } } 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 @@ -230,16 +200,17 @@ class ParameterStudy2 } /** - * @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. + * @brief Carry out all simulations in the parameter study. @see run. + * Convenience function for a few number of runs, but can use more memory because it stores every single simulaiton. + * If this function causes errors, they may be caused by the simulation not being serializable. In that case, try + * using run_serial, provide a processing function, or make the simulation serializable. + * @return Vector of SimulationGraph for each run. */ template - std::vector run(CreateSimulationFunction&& create_simulation) + std::vector> run(CreateSimulationFunction&& create_simulation) { return run(std::forward(create_simulation), - [](Simulation&& sim, size_t) -> Simulation&& { + [](SimulationT&& sim, size_t) -> SimulationT&& { return std::move(sim); }); } @@ -299,16 +270,71 @@ class ParameterStudy2 } private: + /** + * @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) + { + assert(start_run_idx <= end_run_idx); + + // 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)); + } + + //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 Create a vector with the number of runs each process should make. + * @param num_runs The total number of runs. + * @param num_procs The total number of processes, i.e. the size of MPI_Comm. */ 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); @@ -324,31 +350,14 @@ class ParameterStudy2 RandomNumberGenerator m_rng; ///< The random number generator used by the study. }; -template -ParameterStudy2 make_parameter_study(const Parameters& global_parameters, Time t0, - Time tmax, Step dt, size_t num_runs) -{ - return {global_parameters, t0, tmax, dt, num_runs}; -} - -template -auto make_parameter_study_graph_ode(const Graph>& global_parameters, FP t0, - FP tmax, FP dt, size_t num_runs) -{ - using SimGraph = Graph, MobilityEdge>; - using SimGraphSim = GraphSimulation; - using Params = Graph>; - - return ParameterStudy2{global_parameters, t0, tmax, dt, num_runs}; -} - //sample parameters and create simulation -template -GraphSim make_sampled_graph_simulation( - const Graph>& sampled_graph, - FP t0, FP dt_node_sim, FP dt_graph_sim) +template +auto make_sampled_graph_simulation(const Graph>& sampled_graph, FP t0, + FP dt_node_sim, FP dt_graph_sim) { - typename GraphSim::Graph sim_graph; + 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); @@ -357,8 +366,7 @@ GraphSim make_sampled_graph_simulation( 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)); + return make_mobility_sim(t0, dt_graph_sim, std::move(sim_graph)); } } // namespace mio diff --git a/cpp/models/abm/result_simulation.h b/cpp/models/abm/result_simulation.h index 01e4b6713b..5a51a6e99e 100644 --- a/cpp/models/abm/result_simulation.h +++ b/cpp/models/abm/result_simulation.h @@ -33,11 +33,10 @@ class ResultSimulation : public Simulation public: using Model = M; - /// @brief Create a simulation, copying the model. - ResultSimulation(const Model& m, TimePoint t) - : Simulation(t, Model(m)) + /// @brief Create a simulation, moving the model. + ResultSimulation(Model&& m, TimePoint t) + : Simulation(t, std::move(m)) { - history.log(*this); // set initial results } /** diff --git a/cpp/tests/test_parameter_studies.cpp b/cpp/tests/test_parameter_studies.cpp index b7f63a6852..2d633d667a 100644 --- a/cpp/tests/test_parameter_studies.cpp +++ b/cpp/tests/test_parameter_studies.cpp @@ -17,6 +17,7 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include "memilio/config.h" #include "memilio/utils/parameter_distributions.h" #include "ode_secir/model.h" #include "ode_secir/parameter_space.h" @@ -157,12 +158,12 @@ 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::make_parameter_study_graph_ode>(graph, 0.0, 0.0, 0.5, 1); + mio::ParameterStudy2 study(graph, 0.0, 0.0, 0.5, 1); mio::log_rng_seeds(study.get_rng(), mio::LogLevel::warn); - auto ensemble_results = study.run([](auto&& g, auto t0_, auto dt_, auto) { + 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_); + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); }); auto& results = ensemble_results.at(0); @@ -373,11 +374,11 @@ 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); - auto parameter_study = mio::make_parameter_study>(model, t0, tmax, 0.1, 1); + mio::ParameterStudy2 parameter_study(model, t0, tmax, 0.1, 1); mio::log_rng_seeds(parameter_study.get_rng(), mio::LogLevel::warn); // Run parameter study - auto ensemble_results = parameter_study.run([](auto&& model_, auto t0_, auto dt_, auto) { + 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_); @@ -399,3 +400,7 @@ TEST(ParameterStudies, check_ensemble_run_result) } } } + +// TEST(ParameterStudies, run_mocks) { +// mio::ParameterStudy2 +// } diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index 0fd07553aa..4f87787fb9 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -170,8 +170,7 @@ TEST(TestSaveResult, save_result_with_params) mio::MobilityParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0))); auto num_runs = 3; - auto parameter_study = - mio::make_parameter_study_graph_ode>(graph, 0.0, 2.0, 0.5, num_runs); + mio::ParameterStudy2 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; @@ -186,8 +185,8 @@ TEST(TestSaveResult, save_result_with_params) parameter_study.run( [](auto&& g, auto t0_, auto dt_, auto) { auto copy = g; - return mio::make_sampled_graph_simulation(draw_sample(copy), - t0_, dt_, dt_); + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); }, [&](auto&& results, auto run_idx) { auto results_graph = results.get_graph(); @@ -307,8 +306,7 @@ TEST(TestSaveResult, save_percentiles_and_sums) indices_save_edges)); auto num_runs = 3; - auto parameter_study = - mio::make_parameter_study_graph_ode>(graph, 0.0, 2.0, 0.5, num_runs); + mio::ParameterStudy2 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; @@ -324,8 +322,8 @@ TEST(TestSaveResult, save_percentiles_and_sums) parameter_study.run( [](auto&& g, auto t0_, auto dt_, auto) { auto copy = g; - return mio::make_sampled_graph_simulation(draw_sample(copy), - t0_, dt_, dt_); + return mio::make_sampled_graph_simulation>(draw_sample(copy), t0_, + dt_, dt_); }, [&](auto&& results, auto /*run_idx*/) { auto results_graph = results.get_graph(); diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h b/pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h deleted file mode 100644 index 8cb60ffe12..0000000000 --- a/pycode/memilio-simulation/memilio/simulation/bindings/compartments/parameter_studies.h +++ /dev/null @@ -1,103 +0,0 @@ -/* -* 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 "memilio/compartments/parameter_studies.h" -#include "memilio/compartments/compartmental_model.h" - -#include "pybind11/pybind11.h" -#include "pybind11/stl_bind.h" - -#include -#include -#include -#include -#include - -namespace py = pybind11; - -namespace pymio -{ - -/* - * @brief bind ParameterStudy for any model - */ -template -void bind_GrapParameterStudy( - py::module_& m, std::string const& name, std::vector argnames, - std::function create_simulation = - [](const mio::Graph>& g, double t0, double dt, - size_t) { - using GraphSim = mio::GraphSimulation< - double, - mio::Graph>, mio::MobilityEdge>, - double, double>; - auto copy = g; - return mio::make_sampled_graph_simulation(draw_sample(copy), t0, dt, dt) - }) -{ - assert(sizeof...(RunArgs) == argnames.size()); - using SimulationT = mio::GraphSimulation< - double, mio::Graph>, mio::MobilityEdge>, - double, double>; - using ParametersT = Graph < typename Sim::Model, MobilityParameters; - using StudyT = mio::ParameterStudy2; - using CreateSimulationFunctionT = std::function; - using ProcessSimulationResultFunctionT = std::function; - 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("parameters", py::overload_cast<>(&StudyT::get_parameters), - py::return_value_policy::reference_internal) - .def_property_readonly("rng", &StudyT::get_rng, py::return_value_policy::reference_internal) - .def( - "run", - [](StudyT& self, const ProcessSimulationResultFunctionT& handle_result, RunArgs args...) { - self.run(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), run_idx); - return 0; - }); - }, - py::arg("handle_result_func")) - .def("run", - [](StudyT& self) { //default argument doesn't seem to work with functions - return self.run(create_simulation); - }) - .def( - "run_single", - [](StudyT& self, ProcessSimulationResultFunctionT handle_result) { - self.run(create_simulation, [&handle_result](auto&& r, auto&& run_idx) { - handle_result(std::move(r.nodes()[0].property.get_simulation()), run_idx); - return 0; - }); - }, - py::arg("handle_result_func")) - .def("run_single", [](StudyT& self) { - return filter_graph_results(self.run(create_simulation)); - }); -} - -} // namespace pymio diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index 5fc722ed8d..6d4e71dd9d 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp @@ -49,6 +49,7 @@ #include "pybind11/stl_bind.h" #include "Eigen/Core" #include +#include #include namespace py = pybind11; @@ -56,32 +57,6 @@ 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>, - double, double>>&& 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].get_graph().nodes()[0].property.get_simulation())); - } - return std::move(results); -} - -/// Moves out graphs from GraphSimulation%s. Helps make the new ParameterStudy backwards compatible with older bindings. -template -std::vector extract_graph_from_graph_simulation(std::vector&& result_graph_sims) -{ - std::vector result_graphs; - result_graphs.reserve(result_graph_sims.size()); - for (const SimulationT& sim : result_graph_sims) { - result_graphs.emplace_back(std::move(sim.get_graph())); - } - return result_graphs; -} - /* * @brief bind ParameterStudy for any model */ @@ -91,11 +66,11 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) using GraphT = mio::Graph, mio::MobilityEdge>; using SimulationT = mio::GraphSimulation; using ParametersT = mio::Graph>; - using StudyT = mio::ParameterStudy2; + using StudyT = mio::ParameterStudy2; 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); + return mio::make_sampled_graph_simulation(draw_sample(copy), t0, dt, dt); }; pymio::bind_class(m, name.c_str()) @@ -112,7 +87,7 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) .def( "run", [&create_simulation](StudyT& self, std::function handle_result) { - self.run(create_simulation, [&handle_result](auto&& g, auto&& run_idx) { + 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(). @@ -124,19 +99,24 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) py::arg("handle_result_func")) .def("run", [&create_simulation](StudyT& self) { //default argument doesn't seem to work with functions - return extract_graph_from_graph_simulation(self.run(create_simulation)); + return self.run_serial(create_simulation, [](SimulationT&& result, size_t) { + return std::move(result.get_graph()); + }); }) .def( "run_single", [&create_simulation](StudyT& self, std::function handle_result) { - self.run(create_simulation, [&handle_result](auto&& r, auto&& run_idx) { + 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", [&create_simulation](StudyT& self) { - return filter_graph_results(self.run(create_simulation)); + 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()); + }); }); } diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index 7bb1210402..1a86208697 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -49,31 +49,6 @@ 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>, - double, double>>&& 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].get_graph().nodes()[0].property.get_simulation())); - } - return std::move(results); -} - -/// Moves out graphs from GraphSimulation%s. Helps make the new ParameterStudy backwards compatible with older bindings. -template -std::vector extract_graph_from_graph_simulation(std::vector&& result_graph_sims) -{ - std::vector result_graphs; - result_graphs.reserve(result_graph_sims.size()); - for (const SimulationT& sim : result_graph_sims) { - result_graphs.emplace_back(std::move(sim.get_graph())); - } - return result_graphs; -} /* * @brief bind ParameterStudy for any model @@ -84,7 +59,7 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) using GraphT = mio::Graph, mio::MobilityEdge>; using SimulationT = mio::GraphSimulation; using ParametersT = mio::Graph>; - using StudyT = mio::ParameterStudy2; + using StudyT = mio::ParameterStudy2; pymio::bind_class(m, name.c_str()) .def(py::init(), py::arg("parameters"), py::arg("t0"), @@ -93,7 +68,6 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) .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<>(&StudyT::get_parameters, py::const_), @@ -104,11 +78,11 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) std::function, mio::MobilityEdge>, size_t)> handle_result, bool variant_high) { - self.run( + 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); + 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 @@ -124,22 +98,25 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) "run", [](StudyT& self, bool variant_high) { //default argument doesn't seem to work with functions - return extract_graph_from_graph_simulation( - self.run([variant_high](const ParametersT& g, double t0, double dt, size_t) { + 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); - })); + 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", [](StudyT& self, std::function handle_result, bool variant_high) { - self.run( + 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); + 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.get_graph().nodes()[0].property.get_simulation()), run_idx); @@ -150,12 +127,16 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) .def( "run_single", [](StudyT& self, bool variant_high) { - return filter_graph_results( - self.run([variant_high](const ParametersT& g, double t0, double dt, size_t) { + 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); - })); + 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")); } From 63b7a6d50650d9b50f77aa62a7a7b1a5908e91d4 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 24 Oct 2025 14:52:47 +0200 Subject: [PATCH 17/29] remove developement artifacts --- cpp/memilio/io/io.h | 19 ------------------- cpp/memilio/utils/custom_index_array.h | 13 ------------- 2 files changed, 32 deletions(-) diff --git a/cpp/memilio/io/io.h b/cpp/memilio/io/io.h index 13c2d30045..fdcf9424d4 100644 --- a/cpp/memilio/io/io.h +++ b/cpp/memilio/io/io.h @@ -857,15 +857,6 @@ void serialize(IOContext& io, const T& t) serialize_internal(io, t); } -namespace details -{ -template -using is_serializable_expr_t = decltype(serialize(std::declval(), std::declval())); -} - -template -static constexpr bool is_serializable = is_expression_valid::value; - /** * Restores an object from the data stored in an IO context. * There must be provided for the type T either a free function `deserialize_internal(io, tag)` @@ -890,16 +881,6 @@ IOResult deserialize(IOContext& io, Tag tag) return deserialize_internal(io, tag); } -namespace details -{ -template -using is_deserializable_expr_t = - decltype(std::declval>() = deserialize(std::declval(), std::declval>())); -} - -template -static constexpr bool is_deserializable = is_expression_valid::value; - /** * @brief Returns the current working directory name */ diff --git a/cpp/memilio/utils/custom_index_array.h b/cpp/memilio/utils/custom_index_array.h index e23991db03..d17d496275 100644 --- a/cpp/memilio/utils/custom_index_array.h +++ b/cpp/memilio/utils/custom_index_array.h @@ -23,8 +23,6 @@ #include "memilio/math/eigen_util.h" #include "memilio/utils/index.h" #include "memilio/utils/stl_util.h" -#include -#include namespace { @@ -92,17 +90,6 @@ std::enable_if_t<(I < (Index::size - 1)), std::pair> flatten_ind return {val + (size_t)mio::get(indices) * prod, prod * (size_t)mio::get(dimensions)}; } -template -void assign_from_vector(Index& index, const std::vector& values) -{ - assert(values.size() == sizeof...(Tags)); - - if constexpr (I < sizeof...(Tags)) { - get(index) = values[I]; - assign_from_vector(index, values); - } -} - template struct is_random_access_iterator : std::is_base_of::iterator_category, std::random_access_iterator_tag> { From 9bae572c2c12b9141468b9eed763bd73ac602ac5 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Mon, 27 Oct 2025 14:41:08 +0100 Subject: [PATCH 18/29] small cleanup, add test --- cpp/memilio/compartments/parameter_studies.h | 42 +++---------- cpp/tests/test_parameter_studies.cpp | 66 +++++++++++++++++++- 2 files changed, 72 insertions(+), 36 deletions(-) diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index f3eb730b66..f87610521e 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -199,49 +199,25 @@ class ParameterStudy2 return ensemble_result; } - /** - * @brief Carry out all simulations in the parameter study. @see run. - * Convenience function for a few number of runs, but can use more memory because it stores every single simulaiton. - * If this function causes errors, they may be caused by the simulation not being serializable. In that case, try - * using run_serial, provide a processing function, or make the simulation serializable. - * @return Vector of SimulationGraph for each run. - */ - template - std::vector> run(CreateSimulationFunction&& create_simulation) - { - return run(std::forward(create_simulation), - [](SimulationT&& sim, size_t) -> SimulationT&& { - return std::move(sim); - }); - } - - /** - * @brief Return the number of total runs that the study will make. - */ + /// @brief Return the number of total runs that the study will make. size_t get_num_runs() const { return m_num_runs; } - /** - * @brief Return the final time point for simulations. - */ + /// @brief Return the final time point for simulations. Time get_tmax() const { return m_tmax; } - /** - * @brief Return the initial time point for simulations. - */ + /// @brief Return the initial time point for simulations. Time get_t0() const { return m_t0; } - /** - * @brief Return the initial step sized used by simulations. - */ + /// @brief Return the initial step sized used by simulations. Time get_dt() const { return m_dt; @@ -261,9 +237,7 @@ class ParameterStudy2 } /** @} */ - /** - * @brief Access the study's random number generator. - */ + /// @brief Access the study's random number generator. RandomNumberGenerator& get_rng() { return m_rng; @@ -324,11 +298,13 @@ class ParameterStudy2 } /** - * @brief Create a vector with the number of runs each process should make. + * @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. */ - std::vector distribute_runs(size_t num_runs, int num_procs) + static std::vector distribute_runs(size_t num_runs, int num_procs) { assert(num_procs > 0); //evenly distribute runs diff --git a/cpp/tests/test_parameter_studies.cpp b/cpp/tests/test_parameter_studies.cpp index 2d633d667a..5140f3b54c 100644 --- a/cpp/tests/test_parameter_studies.cpp +++ b/cpp/tests/test_parameter_studies.cpp @@ -18,13 +18,16 @@ * 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) @@ -401,6 +404,63 @@ TEST(ParameterStudies, check_ensemble_run_result) } } -// TEST(ParameterStudies, run_mocks) { -// mio::ParameterStudy2 -// } +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::ParameterStudy2 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); + } +} From 651f1baea26ce4667480f66e46d81fbfb5f5c02a Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Tue, 28 Oct 2025 10:35:06 +0100 Subject: [PATCH 19/29] change name back to ParameterStudy --- cpp/examples/abm_parameter_study.cpp | 2 +- cpp/examples/ode_secir_parameter_study.cpp | 2 +- cpp/examples/ode_secir_parameter_study_graph.cpp | 2 +- cpp/examples/ode_secir_read_graph.cpp | 2 +- cpp/memilio/compartments/parameter_studies.h | 4 ++-- cpp/tests/test_parameter_studies.cpp | 6 +++--- cpp/tests/test_save_results.cpp | 4 ++-- .../memilio/simulation/bindings/models/osecir.cpp | 2 +- .../memilio/simulation/bindings/models/osecirvvs.cpp | 2 +- 9 files changed, 13 insertions(+), 13 deletions(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index f0a2558d6f..4d609f323c 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -178,7 +178,7 @@ int main() 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::ParameterStudy2 study(0, t0, tmax, mio::abm::TimeSpan(0), num_runs); + 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}); diff --git a/cpp/examples/ode_secir_parameter_study.cpp b/cpp/examples/ode_secir_parameter_study.cpp index 0d2ff9c7d7..7450226ef4 100644 --- a/cpp/examples/ode_secir_parameter_study.cpp +++ b/cpp/examples/ode_secir_parameter_study.cpp @@ -127,7 +127,7 @@ int main() // create study auto num_runs = size_t(3); - mio::ParameterStudy2 parameter_study(model, t0, tmax, dt, num_runs); + mio::ParameterStudy parameter_study(model, t0, tmax, dt, num_runs); // set up for run auto sample_graph = [](const auto& model_, ScalarType t0_, ScalarType dt_, size_t) { diff --git a/cpp/examples/ode_secir_parameter_study_graph.cpp b/cpp/examples/ode_secir_parameter_study_graph.cpp index b0f8202f0b..f17d019f71 100644 --- a/cpp/examples/ode_secir_parameter_study_graph.cpp +++ b/cpp/examples/ode_secir_parameter_study_graph.cpp @@ -262,7 +262,7 @@ int main() 2, 1, Eigen::VectorX::Constant(num_age_groups * (size_t)mio::osecir::InfectionState::Count, 0.2), indices_save_edges); - mio::ParameterStudy2 parameter_study(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: "); diff --git a/cpp/examples/ode_secir_read_graph.cpp b/cpp/examples/ode_secir_read_graph.cpp index e282ef5468..41cdf89b24 100644 --- a/cpp/examples/ode_secir_read_graph.cpp +++ b/cpp/examples/ode_secir_read_graph.cpp @@ -146,7 +146,7 @@ int main(int argc, char** argv) auto& graph_read = graph_read_result.value(); std::cout << "Running Simulations..." << std::flush; - mio::ParameterStudy2 study(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_, diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index f87610521e..359c30bedf 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -47,7 +47,7 @@ namespace mio * @tparam StepType The type used for time steps, e.g. double or TimeStep. */ template -class ParameterStudy2 +class ParameterStudy { public: using Parameters = ParameterType; @@ -76,7 +76,7 @@ class ParameterStudy2 * @param dt Initial time step of simulations. * @param num_runs Number of simulations that will be created and run. */ - ParameterStudy2(const Parameters& parameters, Time t0, Time tmax, Step dt, size_t num_runs) + 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} diff --git a/cpp/tests/test_parameter_studies.cpp b/cpp/tests/test_parameter_studies.cpp index 5140f3b54c..01e9f6e6a1 100644 --- a/cpp/tests/test_parameter_studies.cpp +++ b/cpp/tests/test_parameter_studies.cpp @@ -161,7 +161,7 @@ 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))); - mio::ParameterStudy2 study(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 ensemble_results = study.run_serial([](auto&& g, auto t0_, auto dt_, auto) { auto copy = g; @@ -377,7 +377,7 @@ 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::ParameterStudy2 parameter_study(model, t0, tmax, 0.1, 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 @@ -444,7 +444,7 @@ TEST(ParameterStudies, mocked_run) return s.tmax + i; }; const double process_sim_result = (num_runs * tmax) + num_runs * (num_runs - 1) / 2.; - mio::ParameterStudy2 study(params, t0, tmax, dt, num_runs); + 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); diff --git a/cpp/tests/test_save_results.cpp b/cpp/tests/test_save_results.cpp index 4f87787fb9..f706327ce6 100644 --- a/cpp/tests/test_save_results.cpp +++ b/cpp/tests/test_save_results.cpp @@ -170,7 +170,7 @@ TEST(TestSaveResult, save_result_with_params) mio::MobilityParameters(Eigen::VectorXd::Constant(Eigen::Index(num_groups * 10), 1.0))); auto num_runs = 3; - mio::ParameterStudy2 parameter_study(graph, 0.0, 2.0, 0.5, num_runs); + 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; @@ -306,7 +306,7 @@ TEST(TestSaveResult, save_percentiles_and_sums) indices_save_edges)); auto num_runs = 3; - mio::ParameterStudy2 parameter_study(graph, 0.0, 2.0, 0.5, num_runs); + 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; diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index 6d4e71dd9d..881c93b088 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp @@ -66,7 +66,7 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) using GraphT = mio::Graph, mio::MobilityEdge>; using SimulationT = mio::GraphSimulation; using ParametersT = mio::Graph>; - using StudyT = mio::ParameterStudy2; + using StudyT = mio::ParameterStudy; const auto create_simulation = [](const ParametersT& g, double t0, double dt, size_t) { auto copy = g; diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index 1a86208697..9c276310b8 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -59,7 +59,7 @@ void bind_ParameterStudy(py::module_& m, std::string const& name) using GraphT = mio::Graph, mio::MobilityEdge>; using SimulationT = mio::GraphSimulation; using ParametersT = mio::Graph>; - using StudyT = mio::ParameterStudy2; + using StudyT = mio::ParameterStudy; pymio::bind_class(m, name.c_str()) .def(py::init(), py::arg("parameters"), py::arg("t0"), From ea3acdea0b170f83cd0ff48f2766a18c9578f8cf Mon Sep 17 00:00:00 2001 From: Sascha <51127093+xsaschako@users.noreply.github.com> Date: Wed, 29 Oct 2025 11:33:05 +0100 Subject: [PATCH 20/29] update parameters to run --- cpp/CMakeLists.txt | 2 +- cpp/examples/abm_parameter_study.cpp | 10 +++++----- cpp/models/abm/common_abm_loggers.h | 17 +++++++++-------- 3 files changed, 15 insertions(+), 14 deletions(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 3ce345b3fd..8eb526845d 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -18,7 +18,7 @@ option(MEMILIO_SANITIZE_ADDRESS "Enable address sanitizer." OFF) option(MEMILIO_SANITIZE_UNDEFINED "Enable undefined behavior sanitizer." OFF) option(MEMILIO_BUILD_SHARED_LIBS "Build memilio as a shared library." ON) option(MEMILIO_BUILD_STATIC_LIBS "Build memilio as a static library." ON) -option(MEMILIO_ENABLE_MPI "Build memilio with MPI." OFF) +option(MEMILIO_ENABLE_MPI "Build memilio with MPI." ON) option(MEMILIO_ENABLE_OPENMP "Enable Multithreading with OpenMP." OFF) option(MEMILIO_ENABLE_WARNINGS "Build memilio with warnings." ON) option(MEMILIO_ENABLE_WARNINGS_AS_ERRORS "Build memilio with warnings as errors." ON) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 4d609f323c..cf822c8265 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -62,7 +62,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) model.parameters.check_constraints(); // There are 10 households for each household group. - int n_households = 10; + int n_households = 100000; // 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. @@ -119,7 +119,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) 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 end_date = mio::abm::TimePoint(0) + mio::abm::days(40); auto test_type = mio::abm::TestType::Antigen; auto test_parameters = model.parameters.get()[test_type]; auto testing_criteria_work = mio::abm::TestingCriteria(); @@ -159,7 +159,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) } // During the lockdown, social events are closed for 90% of people. - auto t_lockdown = mio::abm::TimePoint(0) + mio::abm::days(10); + auto t_lockdown = mio::abm::TimePoint(0) + mio::abm::days(40); mio::abm::close_social_events(t_lockdown, 0.9, model.parameters); return model; @@ -173,9 +173,9 @@ int main() // Set start and end time for the simulation. auto t0 = mio::abm::TimePoint(0); - auto tmax = t0 + mio::abm::days(5); + auto tmax = t0 + mio::abm::days(60); // auto sim = mio::abm::Simulation(t0, std::move(model)); - const size_t num_runs = 3; + const size_t num_runs = 128; // 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); diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index 0d490437f0..ca7db0f52e 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -166,7 +166,8 @@ struct LogDataForMobility : mio::LogAlways { * @brief Logger to log the TimeSeries of the number of Person%s in an #InfectionState. */ struct LogInfectionState : mio::LogAlways { - using Type = std::pair>; + using Type = std::pair; + /** * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. * @param[in] sim The simulation of the abm. @@ -174,15 +175,15 @@ struct LogInfectionState : mio::LogAlways { */ static Type log(const mio::abm::Simulation<>& sim) { + Eigen::VectorXd sum = Eigen::VectorXd::Zero(Eigen::Index(mio::abm::InfectionState::Count)); + auto curr_time = sim.get_time(); - Eigen::VectorX sum = - Eigen::VectorX::Zero(Eigen::Index(mio::abm::InfectionState::Count)); - auto curr_time = sim.get_time(); - PRAGMA_OMP(for) - for (auto& location : sim.get_model().get_locations()) { + // Otherwise log accordingly + for (auto&& person : sim.get_model().get_persons()) { for (uint32_t inf_state = 0; inf_state < (int)mio::abm::InfectionState::Count; inf_state++) { - sum[inf_state] += sim.get_model().get_subpopulation(location.get_id(), curr_time, - mio::abm::InfectionState(inf_state)); + if (person.get_infection_state(curr_time) == mio::abm::InfectionState(inf_state)) { + sum[inf_state] += 1; + } } } return std::make_pair(curr_time, sum); From f23f5ff5b3f0f06134a3cbb0589e61843f84aa2d Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Wed, 29 Oct 2025 13:03:37 +0100 Subject: [PATCH 21/29] Revert "update parameters to run" This reverts commit ea3acdea0b170f83cd0ff48f2766a18c9578f8cf. --- cpp/CMakeLists.txt | 2 +- cpp/examples/abm_parameter_study.cpp | 10 +++++----- cpp/models/abm/common_abm_loggers.h | 17 ++++++++--------- 3 files changed, 14 insertions(+), 15 deletions(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 8eb526845d..3ce345b3fd 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -18,7 +18,7 @@ option(MEMILIO_SANITIZE_ADDRESS "Enable address sanitizer." OFF) option(MEMILIO_SANITIZE_UNDEFINED "Enable undefined behavior sanitizer." OFF) option(MEMILIO_BUILD_SHARED_LIBS "Build memilio as a shared library." ON) option(MEMILIO_BUILD_STATIC_LIBS "Build memilio as a static library." ON) -option(MEMILIO_ENABLE_MPI "Build memilio with MPI." ON) +option(MEMILIO_ENABLE_MPI "Build memilio with MPI." OFF) option(MEMILIO_ENABLE_OPENMP "Enable Multithreading with OpenMP." OFF) option(MEMILIO_ENABLE_WARNINGS "Build memilio with warnings." ON) option(MEMILIO_ENABLE_WARNINGS_AS_ERRORS "Build memilio with warnings as errors." ON) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index cf822c8265..4d609f323c 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -62,7 +62,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) model.parameters.check_constraints(); // There are 10 households for each household group. - int n_households = 100000; + 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. @@ -119,7 +119,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) 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(40); + 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(); @@ -159,7 +159,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) } // During the lockdown, social events are closed for 90% of people. - auto t_lockdown = mio::abm::TimePoint(0) + mio::abm::days(40); + auto t_lockdown = mio::abm::TimePoint(0) + mio::abm::days(10); mio::abm::close_social_events(t_lockdown, 0.9, model.parameters); return model; @@ -173,9 +173,9 @@ int main() // Set start and end time for the simulation. auto t0 = mio::abm::TimePoint(0); - auto tmax = t0 + mio::abm::days(60); + auto tmax = t0 + mio::abm::days(5); // auto sim = mio::abm::Simulation(t0, std::move(model)); - const size_t num_runs = 128; + 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); diff --git a/cpp/models/abm/common_abm_loggers.h b/cpp/models/abm/common_abm_loggers.h index ca7db0f52e..0d490437f0 100644 --- a/cpp/models/abm/common_abm_loggers.h +++ b/cpp/models/abm/common_abm_loggers.h @@ -166,8 +166,7 @@ struct LogDataForMobility : mio::LogAlways { * @brief Logger to log the TimeSeries of the number of Person%s in an #InfectionState. */ struct LogInfectionState : mio::LogAlways { - using Type = std::pair; - + using Type = std::pair>; /** * @brief Log the TimeSeries of the number of Person%s in an #InfectionState. * @param[in] sim The simulation of the abm. @@ -175,15 +174,15 @@ struct LogInfectionState : mio::LogAlways { */ static Type log(const mio::abm::Simulation<>& sim) { - Eigen::VectorXd sum = Eigen::VectorXd::Zero(Eigen::Index(mio::abm::InfectionState::Count)); - auto curr_time = sim.get_time(); - // Otherwise log accordingly - for (auto&& person : sim.get_model().get_persons()) { + Eigen::VectorX sum = + Eigen::VectorX::Zero(Eigen::Index(mio::abm::InfectionState::Count)); + auto curr_time = sim.get_time(); + PRAGMA_OMP(for) + for (auto& location : sim.get_model().get_locations()) { for (uint32_t inf_state = 0; inf_state < (int)mio::abm::InfectionState::Count; inf_state++) { - if (person.get_infection_state(curr_time) == mio::abm::InfectionState(inf_state)) { - sum[inf_state] += 1; - } + sum[inf_state] += sim.get_model().get_subpopulation(location.get_id(), curr_time, + mio::abm::InfectionState(inf_state)); } } return std::make_pair(curr_time, sum); From 6b21559975faa9ca7f714af0fad696ace58d7005 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 6 Nov 2025 13:40:50 +0100 Subject: [PATCH 22/29] update rng log message --- cpp/memilio/compartments/parameter_studies.h | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 359c30bedf..9f81234230 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -281,12 +281,16 @@ class ParameterStudy //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()); + + [[maybe_unused]] const uint64_t create_counter = (thread_local_rng().get_counter() - run_rng_counter).get(); + log_debug("ParameterStudy: Generated {} random numbers creating simulation #{}.", create_counter, run_idx); //perform run sim.advance(m_tmax); + log_debug("ParameterStudy: Generated {} random numbers running simulation #{}.", + run_rng_counter.get() - create_counter, run_idx); + //handle result and store ensemble_result.emplace_back(process_simulation_result(std::move(sim), run_idx)); } @@ -326,7 +330,13 @@ class ParameterStudy RandomNumberGenerator m_rng; ///< The random number generator used by the study. }; -//sample parameters and create simulation +/** + * @brief Create a GraphSimulation from a parameter graph. + * @param[in] sampled_graph A graph of models as nodes and mobility parameters as edges, with pre-sampled values. + * @param[in] t0 Start time of the graph simulation. + * @param[in] dt_node_sim (Initial) time step used by each node in the GraphSimulation. + * @param[in] dt_graph_sim Time step used by the GraphSimulation itself. + */ template auto make_sampled_graph_simulation(const Graph>& sampled_graph, FP t0, FP dt_node_sim, FP dt_graph_sim) From 7dc86ac17f16eaf4c8bc4047569f93a3dfda3f5c Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 6 Nov 2025 16:16:06 +0100 Subject: [PATCH 23/29] [ci skip] Apply suggestions from code review Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com> --- cpp/examples/abm_parameter_study.cpp | 6 +++--- cpp/tests/test_abm_simulation.cpp | 2 +- cpp/tests/test_io_cli.cpp | 2 +- cpp/tests/test_parameter_studies.cpp | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 4d609f323c..aac074dfc8 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -52,13 +52,13 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& 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) + // Set the age groups that can go to school; here this 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. + // Check if the parameters satisfy their constraints. model.parameters.check_constraints(); // There are 10 households for each household group. @@ -90,7 +90,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) 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. + // Maximum contacts 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. diff --git a/cpp/tests/test_abm_simulation.cpp b/cpp/tests/test_abm_simulation.cpp index 856456baf8..8db16d72eb 100644 --- a/cpp/tests/test_abm_simulation.cpp +++ b/cpp/tests/test_abm_simulation.cpp @@ -161,7 +161,7 @@ TEST(TestSimulation, ResultSimulation) 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 + // run simulation. expect one timepoint 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); diff --git a/cpp/tests/test_io_cli.cpp b/cpp/tests/test_io_cli.cpp index 0721044ceb..4ecf31baba 100644 --- a/cpp/tests/test_io_cli.cpp +++ b/cpp/tests/test_io_cli.cpp @@ -315,7 +315,7 @@ TEST_F(TestCLI, test_set_param) EXPECT_THAT(print_wrap(result), IsFailure(mio::StatusCode::InvalidType)); 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")); + result = set.set_param(id_B, std::string("this 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()); diff --git a/cpp/tests/test_parameter_studies.cpp b/cpp/tests/test_parameter_studies.cpp index 01e9f6e6a1..91e51f738e 100644 --- a/cpp/tests/test_parameter_studies.cpp +++ b/cpp/tests/test_parameter_studies.cpp @@ -456,7 +456,7 @@ TEST(ParameterStudies, mocked_run) 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 + // Note: currently the tests do 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}) { From 3db2bcd0c25ab19ceed3b57e60715ca4dbd46f98 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 7 Nov 2025 09:31:46 +0100 Subject: [PATCH 24/29] remove out of scope index ctor --- cpp/memilio/utils/index.h | 46 --------------------------------------- 1 file changed, 46 deletions(-) diff --git a/cpp/memilio/utils/index.h b/cpp/memilio/utils/index.h index ac76fe3b38..cba56b1d57 100644 --- a/cpp/memilio/utils/index.h +++ b/cpp/memilio/utils/index.h @@ -21,11 +21,7 @@ #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 { @@ -33,41 +29,6 @@ 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 @@ -178,13 +139,6 @@ class Index } public: - template ::value>> - Index(IndexArgs... _index_args) - : indices(details::merge_indices(_index_args...)) - { - } - // comparison operators bool operator==(Index const& other) const { From 67dfb49e2d80c27cceb873fb71f9affc3b13f625 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 7 Nov 2025 09:47:01 +0100 Subject: [PATCH 25/29] fix comment --- cpp/examples/abm_parameter_study.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index aac074dfc8..7e400648f7 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -174,7 +174,7 @@ int main() // 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)); + // Set the number of simulations to run in the study 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. From 5c03d2969b8c100b5c996596101e8d383fac82e5 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 7 Nov 2025 09:49:26 +0100 Subject: [PATCH 26/29] Apply suggestions from code review Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com> --- cpp/memilio/compartments/parameter_studies.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index 9f81234230..ae85c3ff99 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -87,7 +87,7 @@ class ParameterStudy /** * @brief Run all simulations in serial. - * @param[in] create_simulation A callable sampling the study's parameters and return a simulation. + * @param[in] create_simulation A callable sampling the study's parameters and returning 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. * @@ -124,7 +124,7 @@ class ParameterStudy /** * @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] create_simulation A callable sampling the study's parameters and returning 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. * From 665f8c4738f662c0c723809f0e2a5d42d3c682dc Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 7 Nov 2025 18:40:44 +0100 Subject: [PATCH 27/29] [ci skip] clarify comment on abm study creation --- cpp/examples/abm_parameter_study.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 7e400648f7..9311348270 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -177,7 +177,10 @@ int main() // Set the number of simulations to run in the study 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. + // Create a parameter study. + // Note that the study for the ABM currently does not make use of the arguments "parameters" or "dt", as we create + // a new model for each simulation. Hence we set both arguments to 0. + // This is mostly due to https://github.com/SciCompMod/memilio/issues/1400 mio::ParameterStudy study(0, t0, tmax, mio::abm::TimeSpan(0), num_runs); // Optional: set seeds to get reproducable results From 35fb3a12d8d7fb522c90d0058d2e2e647fb8da83 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 7 Nov 2025 18:42:31 +0100 Subject: [PATCH 28/29] [ci skip] review suggestion --- cpp/examples/abm_parameter_study.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 9311348270..8446465c5f 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -55,7 +55,7 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng) // Set the age groups that can go to school; here this 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) + // Set the age groups that can go to work; here these are 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 constraints. From 7819b98c1bc8c672f9a79ff0e8829195da973ba2 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Mon, 10 Nov 2025 13:17:48 +0100 Subject: [PATCH 29/29] [ci skip] update parameter study class description --- cpp/memilio/compartments/parameter_studies.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index ae85c3ff99..516476fc22 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -41,10 +41,11 @@ namespace mio { /** - * @brief Class used to performs multiple simulation runs with randomly sampled parameters. + * @brief Class used to perform multiple simulation runs with randomly sampled parameters. + * Note that the type of simulation is not determined until calling one of the run functions. * @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. + * @tparam TimeType The time type used by the simulation, e.g. double or TimePoint. + * @tparam StepType The time step type used by the simulation, e.g. double or TimeStep. May be the same as TimeType. */ template class ParameterStudy