Skip to content

Commit 8f90ce7

Browse files
reneSchmreneSchm
authored
1400 The ABMs RNG is not used after initialization (#1534)
- Change PersonalRandomNumberGenerator to use the key from another RNG (usually the model's). - Remove the now unused key property and related methods from Person. - Modify the ABM parameter study example, such that it can return reproducible results with MPI or OpenMP. - Add a const version of abm::Model::get_rng. - Add new abm::Model::reset_rng function. - Fix interpolate_simulation_result returning a 0.0 time point with negative sign. - Add a stream operator for std::vector, allowing for easier printing. - Make git ignore the example_results directory.
1 parent 20cee90 commit 8f90ce7

24 files changed

Lines changed: 232 additions & 144 deletions

‎.gitignore‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -288,3 +288,6 @@ docs/source/cppapi
288288
.vscode/
289289

290290
# End of https://www.gitignore.io/api/c++,node,python
291+
292+
# Output directory from memilio examples
293+
example_results/

‎cpp/benchmarks/abm.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,7 @@ mio::abm::Simulation<> make_simulation(size_t num_persons, std::initializer_list
6868

6969
//infections and masks
7070
for (auto& person : model.get_persons()) {
71-
auto prng = mio::abm::PersonalRandomNumberGenerator(person);
71+
auto prng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
7272
//some % of people are infected, large enough to have some infection activity without everyone dying
7373
auto pct_infected = 0.05;
7474
if (mio::UniformDistribution<ScalarType>::get_instance()(prng, 0.0, 1.0) < pct_infected) {

‎cpp/examples/abm_history_object.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,7 @@ int main()
143143
// The infection states are chosen randomly.
144144
auto persons = model.get_persons();
145145
for (auto& person : persons) {
146-
auto rng = mio::abm::PersonalRandomNumberGenerator(person);
146+
auto rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
147147
mio::abm::InfectionState infection_state =
148148
(mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1));
149149
if (infection_state != mio::abm::InfectionState::Susceptible)

‎cpp/examples/abm_minimal.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -120,7 +120,7 @@ int main()
120120
for (auto& person : model.get_persons()) {
121121
mio::abm::InfectionState infection_state = mio::abm::InfectionState(
122122
mio::DiscreteDistribution<size_t>::get_instance()(mio::thread_local_rng(), infection_distribution));
123-
auto rng = mio::abm::PersonalRandomNumberGenerator(person);
123+
auto rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
124124
if (infection_state != mio::abm::InfectionState::Susceptible) {
125125
person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
126126
model.parameters, start_date, infection_state));

‎cpp/examples/abm_parameter_study.cpp‎

Lines changed: 43 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -46,8 +46,8 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng)
4646
const auto age_group_15_to_34 = mio::AgeGroup(2);
4747
const auto age_group_35_to_59 = mio::AgeGroup(3);
4848
// Create the model with 4 age groups.
49-
auto model = mio::abm::Model(num_age_groups);
50-
model.get_rng() = rng;
49+
auto model = mio::abm::Model(num_age_groups);
50+
model.reset_rng(rng.get_seeds());
5151

5252
// Set same infection parameter for all age groups. For example, the incubation period is log normally distributed with parameters 4 and 1.
5353
model.parameters.get<mio::abm::TimeExposedToNoSymptoms>() = mio::ParameterDistributionLogNormal(4., 1.);
@@ -131,8 +131,8 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng)
131131
std::vector<ScalarType> infection_distribution{0.5, 0.3, 0.05, 0.05, 0.05, 0.05, 0.0, 0.0};
132132
for (auto& person : model.get_persons()) {
133133
mio::abm::InfectionState infection_state = mio::abm::InfectionState(
134-
mio::DiscreteDistribution<size_t>::get_instance()(mio::thread_local_rng(), infection_distribution));
135-
auto person_rng = mio::abm::PersonalRandomNumberGenerator(person);
134+
mio::DiscreteDistribution<size_t>::get_instance()(model.get_rng(), infection_distribution));
135+
auto person_rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
136136
if (infection_state != mio::abm::InfectionState::Susceptible) {
137137
person.add_new_infection(mio::abm::Infection(person_rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
138138
model.parameters, start_date, infection_state));
@@ -175,53 +175,64 @@ int main()
175175
auto tmax = t0 + mio::abm::days(5);
176176
// Set the number of simulations to run in the study
177177
const size_t num_runs = 3;
178+
// Set up an RNG
179+
mio::RandomNumberGenerator rng;
178180

179-
// Create a parameter study.
180-
// Note that the study for the ABM currently does not make use of the arguments "parameters" or "dt", as we create
181-
// a new model for each simulation. Hence we set both arguments to 0.
182-
// This is mostly due to https://github.com/SciCompMod/memilio/issues/1400
183-
mio::ParameterStudy study(0, t0, tmax, mio::abm::TimeSpan(0), num_runs);
181+
// Optional: set seeds to get reproducible results
182+
// rng.seed({1234, 2345, 124324, 7567, 34534, 7});
183+
184+
// Always synchronise seeds before creating the model, otherwise different MPI ranks will create different models
185+
// If you do want to use a different model each run, move the model setup from the ParameterStudy creation into the
186+
// simulation setup (i.e. into the create_simulation function argument of ParameterStudy::run)
187+
rng.synchronize();
188+
auto model = make_model(rng);
184189

185-
// Optional: set seeds to get reproducable results
186-
// study.get_rng().seed({12341234, 53456, 63451, 5232576, 84586, 52345});
190+
// Create a parameter study.
191+
// Note that the ABM's step size is fixed, hence we set the effectively unused "dt" argument to 0
192+
mio::ParameterStudy study(std::move(model), t0, tmax, mio::abm::TimeSpan(0), num_runs);
193+
study.get_rng() = rng; // use the same RNG as the model
187194

188195
const std::string result_dir = mio::path_join(mio::base_dir(), "example_results");
189196
if (!mio::create_directory(result_dir)) {
190197
mio::log_error("Could not create result directory \"{}\".", result_dir);
191198
return 1;
192199
}
193200

201+
// Run the study
202+
// The first lambda ("create_simulation" argument) sets up the simulation, the second ("process_simulation_result")
203+
// allows us to process each simulations result. Be mindful of the memory used for storing these results!
194204
auto ensemble_results = study.run(
195-
[](auto, auto t0_, auto, size_t) {
196-
return mio::abm::ResultSimulation(make_model(mio::thread_local_rng()), t0_);
205+
[](auto&& model_, auto t0_, auto, size_t run_idx) {
206+
auto copy = model_;
207+
// using half of the RNG counter for the run index (soft-) limits both the number of runs and RNG draws
208+
// per person to 2^16 = 65536
209+
const auto ctr = mio::Counter<uint32_t>(static_cast<uint32_t>(run_idx) << 16);
210+
copy.reset_rng(ctr);
211+
return mio::abm::ResultSimulation(std::move(copy), t0_);
197212
},
198-
[result_dir](auto&& sim, auto&& run_idx) {
213+
[&result_dir](auto&& sim, auto&& run_idx) {
199214
auto interpolated_result = mio::interpolate_simulation_result(sim.get_result());
200215
std::string outpath = mio::path_join(result_dir, "abm_minimal_run_" + std::to_string(run_idx) + ".txt");
201216
std::ofstream outfile_run(outpath);
202217
sim.get_result().print_table(outfile_run, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4);
203218
std::cout << "Results written to " << outpath << std::endl;
204-
auto params = std::vector<mio::abm::Model>{};
205219
return std::vector{interpolated_result};
206220
});
207221

208-
if (ensemble_results.size() > 0) {
209-
auto ensemble_results_p05 = ensemble_percentile(ensemble_results, 0.05);
210-
auto ensemble_results_p25 = ensemble_percentile(ensemble_results, 0.25);
211-
auto ensemble_results_p50 = ensemble_percentile(ensemble_results, 0.50);
212-
auto ensemble_results_p75 = ensemble_percentile(ensemble_results, 0.75);
213-
auto ensemble_results_p95 = ensemble_percentile(ensemble_results, 0.95);
214-
215-
mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups,
216-
mio::path_join(result_dir, "Results_" + std::string("p05") + ".h5")));
217-
mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups,
218-
mio::path_join(result_dir, "Results_" + std::string("p25") + ".h5")));
219-
mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups,
220-
mio::path_join(result_dir, "Results_" + std::string("p50") + ".h5")));
221-
mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups,
222-
mio::path_join(result_dir, "Results_" + std::string("p75") + ".h5")));
223-
mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups,
224-
mio::path_join(result_dir, "Results_" + std::string("p95") + ".h5")));
222+
// The study collects all results on the root rank, so we only process the results there
223+
if (mio::mpi::is_root()) {
224+
const auto write_percentile = [&](double p) {
225+
std::ofstream out(mio::path_join(result_dir, fmt::format("Results_p{:0<4.2}.txt", p)));
226+
auto ensemble_percentiles = ensemble_percentile(ensemble_results, p);
227+
ensemble_percentiles.front().print_table(out, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7,
228+
4);
229+
};
230+
231+
write_percentile(0.05);
232+
write_percentile(0.25);
233+
write_percentile(0.50);
234+
write_percentile(0.75);
235+
write_percentile(0.95);
225236
}
226237

227238
mio::mpi::finalize();

‎cpp/examples/graph_abm.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -201,7 +201,7 @@ int main()
201201
for (auto& person : model1.get_persons()) {
202202
mio::abm::InfectionState infection_state = mio::abm::InfectionState(
203203
mio::DiscreteDistribution<size_t>::get_instance()(model1.get_rng(), infection_distribution_m1));
204-
auto rng = mio::abm::PersonalRandomNumberGenerator(person);
204+
auto rng = mio::abm::PersonalRandomNumberGenerator(model1.get_rng(), person);
205205
if (infection_state != mio::abm::InfectionState::Susceptible) {
206206
person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
207207
model1.parameters, start_date, infection_state));
@@ -231,7 +231,7 @@ int main()
231231
for (auto& person : model2.get_persons()) {
232232
mio::abm::InfectionState infection_state = mio::abm::InfectionState(
233233
mio::DiscreteDistribution<size_t>::get_instance()(model2.get_rng(), infection_distribution_m2));
234-
auto rng = mio::abm::PersonalRandomNumberGenerator(person);
234+
auto rng = mio::abm::PersonalRandomNumberGenerator(model2.get_rng(), person);
235235
if (infection_state != mio::abm::InfectionState::Susceptible) {
236236
person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
237237
model2.parameters, start_date, infection_state));

‎cpp/memilio/data/analyze_result.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -192,13 +192,13 @@ TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_re
192192
const auto t0 = simulation_result.get_time(0);
193193
const auto t_max = simulation_result.get_last_time();
194194
// add an additional day, if the first time point is within tolerance of floor(t0)
195-
const auto day0 = ceil(t0 - abs_tol);
195+
const auto day0 = static_cast<int>(ceil(t0 - abs_tol));
196196
// add an additional day, if the last time point is within tolerance of ceil(tmax)
197-
const auto day_max = floor(t_max + abs_tol);
197+
const auto day_max = static_cast<int>(floor(t_max + abs_tol));
198198

199199
// create interpolation_times vector with all days between day0 and day_max
200-
std::vector<FP> tps(static_cast<int>(day_max) - static_cast<int>(day0) + 1);
201-
std::iota(tps.begin(), tps.end(), day0);
200+
std::vector<FP> tps(day_max - day0 + 1);
201+
std::iota(tps.begin(), tps.end(), FP(day0));
202202

203203
return interpolate_simulation_result<FP>(simulation_result, tps);
204204
}

‎cpp/memilio/utils/stl_util.h‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,22 @@
3434
#include <cassert>
3535
#include <memory>
3636

37+
/// @brief Stream operator for writing std::vector to an std::ostream. Requires that T has a viable overload itself.
38+
template <class T>
39+
std::ostream& operator<<(std::ostream& out, const std::vector<T>& vec)
40+
{
41+
out << "{";
42+
// use size - 1 and handle the last entry separately, for nicer separator usage
43+
for (size_t i = 0; i < vec.size() - 1; i++) {
44+
out << vec[i] << ", ";
45+
}
46+
if (!vec.empty()) {
47+
out << vec.back();
48+
}
49+
out << "}";
50+
return out;
51+
}
52+
3753
namespace mio
3854
{
3955

‎cpp/models/abm/model.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -108,7 +108,7 @@ void Model::perform_mobility(TimePoint t, TimeSpan dt)
108108
for (uint32_t person_index = 0; person_index < num_persons; ++person_index) {
109109
if (m_activeness_statuses[person_index]) {
110110
Person& person = m_persons[person_index];
111-
auto personal_rng = PersonalRandomNumberGenerator(person);
111+
auto personal_rng = PersonalRandomNumberGenerator(m_rng, person);
112112

113113
auto try_mobility_rule = [&](auto rule) -> bool {
114114
// run mobility rule and check if change of location can actually happen
@@ -186,7 +186,7 @@ void Model::perform_mobility(TimePoint t, TimeSpan dt)
186186
m_trip_list.increase_index()) {
187187
auto& trip = m_trip_list.get_next_trip();
188188
auto& person = get_person(trip.person_id);
189-
auto personal_rng = PersonalRandomNumberGenerator(person);
189+
auto personal_rng = PersonalRandomNumberGenerator(m_rng, person);
190190
// skip the trip if the person is in quarantine or is dead
191191
if (person.is_in_quarantine(t, parameters) || person.get_infection_state(t) == InfectionState::Dead) {
192192
continue;

‎cpp/models/abm/model.h‎

Lines changed: 30 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -60,7 +60,7 @@ class Model
6060
using ActivenessIterator = std::vector<bool>::iterator;
6161
using ConstActivenessIterator = std::vector<bool>::const_iterator;
6262
using MobilityRuleType = LocationType (*)(PersonalRandomNumberGenerator&, const Person&, TimePoint, TimeSpan,
63-
const Parameters&);
63+
const Parameters&);
6464

6565
using Compartments = mio::abm::InfectionState;
6666
/**
@@ -333,11 +333,39 @@ class Model
333333
* Get the RandomNumberGenerator used by this Model for random events.
334334
* Persons use their own generators with the same key as the global one.
335335
* @return The random number generator.
336+
* @{
336337
*/
337338
RandomNumberGenerator& get_rng()
338339
{
339340
return m_rng;
340341
}
342+
const RandomNumberGenerator& get_rng() const
343+
{
344+
return m_rng;
345+
}
346+
/** @} */
347+
348+
/**
349+
* @brief Sets the RNG counters of the model and all persons to 0 (or to the optional counter argument).
350+
* @param[in] seeds Optional argument used to overwrite the current seed.
351+
* @param[in] counter Optional argument used as base value for each RNG counter.
352+
* Note: Both the model's and each person's RNG uses a 64 bit counter. Since the persons reserve half of the
353+
* counter for their ID, we only allow specifying the first 32 bits here, even for the model.
354+
* @{
355+
*/
356+
void reset_rng(Counter<uint32_t> counter = {})
357+
{
358+
m_rng.set_counter(Counter<uint64_t>(static_cast<uint64_t>(counter.get())));
359+
for (Person& person : get_persons()) {
360+
person.get_rng_counter() = counter;
361+
}
362+
}
363+
void reset_rng(const std::vector<uint32_t>& seeds, Counter<uint32_t> counter = {})
364+
{
365+
m_rng.seed(seeds);
366+
reset_rng(counter);
367+
}
368+
/** @} */
341369

342370
/**
343371
* Get the model id. Is only relevant for graph abm or hybrid model.
@@ -609,7 +637,7 @@ class Model
609637
compute_exposure_caches(t, dt);
610638
m_are_exposure_caches_valid = true;
611639
}
612-
auto personal_rng = PersonalRandomNumberGenerator(person);
640+
auto personal_rng = PersonalRandomNumberGenerator(m_rng, person);
613641
auto location = person.get_location();
614642
mio::abm::interact(personal_rng, person, get_location(location),
615643
m_local_population_by_age_cache[location.get()], m_air_exposure_rates_cache[location.get()],

0 commit comments

Comments
 (0)