1400 The ABMs RNG is not used after initialization

.gitignore

@@ -288,3 +288,6 @@ docs/source/cppapi
 .vscode/
 
 # End of https://www.gitignore.io/api/c++,node,python
+
+# Output directory from memilio examples
+example_results/

cpp/benchmarks/abm.cpp

@@ -68,7 +68,7 @@ mio::abm::Simulation<> make_simulation(size_t num_persons, std::initializer_list
 
     //infections and masks
     for (auto& person : model.get_persons()) {
-        auto prng = mio::abm::PersonalRandomNumberGenerator(person);
+        auto prng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
         //some % of people are infected, large enough to have some infection activity without everyone dying
         auto pct_infected = 0.05;
         if (mio::UniformDistribution<ScalarType>::get_instance()(prng, 0.0, 1.0) < pct_infected) {

cpp/examples/abm_history_object.cpp

@@ -143,7 +143,7 @@ int main()
     // The infection states are chosen randomly.
     auto persons = model.get_persons();
     for (auto& person : persons) {
-        auto rng = mio::abm::PersonalRandomNumberGenerator(person);
+        auto rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
         mio::abm::InfectionState infection_state =
             (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1));
         if (infection_state != mio::abm::InfectionState::Susceptible)

cpp/examples/abm_minimal.cpp

@@ -120,7 +120,7 @@ int main()
     for (auto& person : model.get_persons()) {
         mio::abm::InfectionState infection_state = mio::abm::InfectionState(
             mio::DiscreteDistribution<size_t>::get_instance()(mio::thread_local_rng(), infection_distribution));
-        auto rng = mio::abm::PersonalRandomNumberGenerator(person);
+        auto rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
         if (infection_state != mio::abm::InfectionState::Susceptible) {
             person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
                                                          model.parameters, start_date, infection_state));

cpp/examples/abm_parameter_study.cpp

@@ -46,8 +46,8 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng)
     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;
+    auto model = mio::abm::Model(num_age_groups);
+    model.reset_rng(rng.get_seeds());
 
     // Set same infection parameter for all age groups. For example, the incubation period is log normally distributed with parameters 4 and 1.
     model.parameters.get<mio::abm::TimeExposedToNoSymptoms>() = mio::ParameterDistributionLogNormal(4., 1.);
@@ -131,8 +131,8 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng)
     std::vector<ScalarType> infection_distribution{0.5, 0.3, 0.05, 0.05, 0.05, 0.05, 0.0, 0.0};
     for (auto& person : model.get_persons()) {
         mio::abm::InfectionState infection_state = mio::abm::InfectionState(
-            mio::DiscreteDistribution<size_t>::get_instance()(mio::thread_local_rng(), infection_distribution));
-        auto person_rng = mio::abm::PersonalRandomNumberGenerator(person);
+            mio::DiscreteDistribution<size_t>::get_instance()(model.get_rng(), infection_distribution));
+        auto person_rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), 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));
@@ -175,53 +175,64 @@ int main()
     auto tmax = t0 + mio::abm::days(5);
     // Set the number of simulations to run in the study
     const size_t num_runs = 3;
+    // Set up an RNG
+    mio::RandomNumberGenerator rng;
 
-    // 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 reproducible results
+    // rng.seed({1234, 2345, 124324, 7567, 34534, 7});
+
+    // Always synchronise seeds before creating the model, otherwise different MPI ranks will create different models
+    // If you do want to use a different model each run, move the model setup from the ParameterStudy creation into the
+    // simulation setup (i.e. into the create_simulation function argument of ParameterStudy::run)
+    rng.synchronize();
+    auto model = make_model(rng);
 
-    // Optional: set seeds to get reproducable results
-    // study.get_rng().seed({12341234, 53456, 63451, 5232576, 84586, 52345});
+    // Create a parameter study.
+    // Note that the ABM's step size is fixed, hence we set the effectively unused "dt" argument to 0
+    mio::ParameterStudy study(std::move(model), t0, tmax, mio::abm::TimeSpan(0), num_runs);
+    study.get_rng() = rng; // use the same RNG as the model
 
     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;
     }
 
+    // Run the study
+    // The first lambda ("create_simulation" argument) sets up the simulation, the second ("process_simulation_result")
+    // allows us to process each simulations result. Be mindful of the memory used for storing these results!
     auto ensemble_results = study.run(
-        [](auto, auto t0_, auto, size_t) {
-            return mio::abm::ResultSimulation(make_model(mio::thread_local_rng()), t0_);
+        [](auto&& model_, auto t0_, auto, size_t run_idx) {
+            auto copy = model_;
+            // using half of the RNG counter for the run index (soft-) limits both the number of runs and RNG draws
+            // per person to 2^16 = 65536
+            const auto ctr = mio::Counter<uint32_t>(static_cast<uint32_t>(run_idx) << 16);
+            copy.reset_rng(ctr);
+            return mio::abm::ResultSimulation(std::move(copy), t0_);
         },
-        [result_dir](auto&& sim, auto&& run_idx) {
+        [&result_dir](auto&& sim, auto&& run_idx) {
             auto interpolated_result = mio::interpolate_simulation_result(sim.get_result());
             std::string outpath = mio::path_join(result_dir, "abm_minimal_run_" + std::to_string(run_idx) + ".txt");
             std::ofstream outfile_run(outpath);
             sim.get_result().print_table(outfile_run, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7, 4);
             std::cout << "Results written to " << outpath << std::endl;
-            auto params = std::vector<mio::abm::Model>{};
             return std::vector{interpolated_result};
         });
 
-    if (ensemble_results.size() > 0) {
-        auto ensemble_results_p05 = ensemble_percentile(ensemble_results, 0.05);
-        auto ensemble_results_p25 = ensemble_percentile(ensemble_results, 0.25);
-        auto ensemble_results_p50 = ensemble_percentile(ensemble_results, 0.50);
-        auto ensemble_results_p75 = ensemble_percentile(ensemble_results, 0.75);
-        auto ensemble_results_p95 = ensemble_percentile(ensemble_results, 0.95);
-
-        mio::unused(save_result(ensemble_results_p05, {0}, num_age_groups,
-                                mio::path_join(result_dir, "Results_" + std::string("p05") + ".h5")));
-        mio::unused(save_result(ensemble_results_p25, {0}, num_age_groups,
-                                mio::path_join(result_dir, "Results_" + std::string("p25") + ".h5")));
-        mio::unused(save_result(ensemble_results_p50, {0}, num_age_groups,
-                                mio::path_join(result_dir, "Results_" + std::string("p50") + ".h5")));
-        mio::unused(save_result(ensemble_results_p75, {0}, num_age_groups,
-                                mio::path_join(result_dir, "Results_" + std::string("p75") + ".h5")));
-        mio::unused(save_result(ensemble_results_p95, {0}, num_age_groups,
-                                mio::path_join(result_dir, "Results_" + std::string("p95") + ".h5")));
+    // The study collects all results on the root rank, so we only process the results there
+    if (mio::mpi::is_root()) {
+        const auto write_percentile = [&](double p) {
+            std::ofstream out(mio::path_join(result_dir, fmt::format("Results_p{:0<4.2}.txt", p)));
+            auto ensemble_percentiles = ensemble_percentile(ensemble_results, p);
+            ensemble_percentiles.front().print_table(out, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7,
+                                                     4);
+        };
+
+        write_percentile(0.05);
+        write_percentile(0.25);
+        write_percentile(0.50);
+        write_percentile(0.75);
+        write_percentile(0.95);
     }
 
     mio::mpi::finalize();

cpp/examples/graph_abm.cpp

@@ -201,7 +201,7 @@ int main()
     for (auto& person : model1.get_persons()) {
         mio::abm::InfectionState infection_state = mio::abm::InfectionState(
             mio::DiscreteDistribution<size_t>::get_instance()(model1.get_rng(), infection_distribution_m1));
-        auto rng = mio::abm::PersonalRandomNumberGenerator(person);
+        auto rng = mio::abm::PersonalRandomNumberGenerator(model1.get_rng(), person);
         if (infection_state != mio::abm::InfectionState::Susceptible) {
             person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
                                                          model1.parameters, start_date, infection_state));
@@ -231,7 +231,7 @@ int main()
     for (auto& person : model2.get_persons()) {
         mio::abm::InfectionState infection_state = mio::abm::InfectionState(
             mio::DiscreteDistribution<size_t>::get_instance()(model2.get_rng(), infection_distribution_m2));
-        auto rng = mio::abm::PersonalRandomNumberGenerator(person);
+        auto rng = mio::abm::PersonalRandomNumberGenerator(model2.get_rng(), person);
         if (infection_state != mio::abm::InfectionState::Susceptible) {
             person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
                                                          model2.parameters, start_date, infection_state));

cpp/memilio/data/analyze_result.h

@@ -192,13 +192,13 @@ TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_re
     const auto t0    = simulation_result.get_time(0);
     const auto t_max = simulation_result.get_last_time();
     // add an additional day, if the first time point is within tolerance of floor(t0)
-    const auto day0 = ceil(t0 - abs_tol);
+    const auto day0 = static_cast<int>(ceil(t0 - abs_tol));
     // add an additional day, if the last time point is within tolerance of ceil(tmax)
-    const auto day_max = floor(t_max + abs_tol);
+    const auto day_max = static_cast<int>(floor(t_max + abs_tol));
 
     // create interpolation_times vector with all days between day0 and day_max
-    std::vector<FP> tps(static_cast<int>(day_max) - static_cast<int>(day0) + 1);
-    std::iota(tps.begin(), tps.end(), day0);
+    std::vector<FP> tps(day_max - day0 + 1);
+    std::iota(tps.begin(), tps.end(), FP(day0));
 
     return interpolate_simulation_result<FP>(simulation_result, tps);
 }

cpp/memilio/utils/stl_util.h

@@ -34,6 +34,22 @@
 #include <cassert>
 #include <memory>
 
+/// @brief Stream operator for writing std::vector to an std::ostream. Requires that T has a viable overload itself.
+template <class T>
+std::ostream& operator<<(std::ostream& out, const std::vector<T>& vec)
+{
+    out << "{";
+    // use size - 1 and handle the last entry separately, for nicer separator usage
+    for (size_t i = 0; i < vec.size() - 1; i++) {
+        out << vec[i] << ", ";
+    }
+    if (!vec.empty()) {
+        out << vec.back();
+    }
+    out << "}";
+    return out;
+}
+
 namespace mio
 {
 

cpp/models/abm/model.cpp

@@ -108,7 +108,7 @@ void Model::perform_mobility(TimePoint t, TimeSpan dt)
     for (uint32_t person_index = 0; person_index < num_persons; ++person_index) {
         if (m_activeness_statuses[person_index]) {
             Person& person    = m_persons[person_index];
-            auto personal_rng = PersonalRandomNumberGenerator(person);
+            auto personal_rng = PersonalRandomNumberGenerator(m_rng, person);
 
             auto try_mobility_rule = [&](auto rule) -> bool {
                 // run mobility rule and check if change of location can actually happen
@@ -186,7 +186,7 @@ void Model::perform_mobility(TimePoint t, TimeSpan dt)
          m_trip_list.increase_index()) {
         auto& trip        = m_trip_list.get_next_trip();
         auto& person      = get_person(trip.person_id);
-        auto personal_rng = PersonalRandomNumberGenerator(person);
+        auto personal_rng = PersonalRandomNumberGenerator(m_rng, person);
         // skip the trip if the person is in quarantine or is dead
         if (person.is_in_quarantine(t, parameters) || person.get_infection_state(t) == InfectionState::Dead) {
             continue;

cpp/models/abm/model.h

@@ -60,7 +60,7 @@ class Model
     using ActivenessIterator      = std::vector<bool>::iterator;
     using ConstActivenessIterator = std::vector<bool>::const_iterator;
     using MobilityRuleType        = LocationType (*)(PersonalRandomNumberGenerator&, const Person&, TimePoint, TimeSpan,
-                                              const Parameters&);
+                                                     const Parameters&);
 
     using Compartments = mio::abm::InfectionState;
     /**
@@ -333,11 +333,39 @@ class Model
     * Get the RandomNumberGenerator used by this Model for random events.
     * Persons use their own generators with the same key as the global one. 
     * @return The random number generator.
+    * @{
     */
     RandomNumberGenerator& get_rng()
     {
         return m_rng;
     }
+    const RandomNumberGenerator& get_rng() const
+    {
+        return m_rng;
+    }
+    /** @} */
+
+    /**
+     * @brief Sets the RNG counters of the model and all persons to 0 (or to the optional counter argument).
+     * @param[in] seeds Optional argument used to overwrite the current seed.
+     * @param[in] counter Optional argument used as base value for each RNG counter.
+     * Note: Both the model's and each person's RNG uses a 64 bit counter. Since the persons reserve half of the
+     * counter for their ID, we only allow specifying the first 32 bits here, even for the model.
+     * @{
+     */
+    void reset_rng(Counter<uint32_t> counter = {})
+    {
+        m_rng.set_counter(Counter<uint64_t>(static_cast<uint64_t>(counter.get())));
+        for (Person& person : get_persons()) {
+            person.get_rng_counter() = counter;
+        }
+    }
+    void reset_rng(const std::vector<uint32_t>& seeds, Counter<uint32_t> counter = {})
+    {
+        m_rng.seed(seeds);
+        reset_rng(counter);
+    }
+    /** @} */
 
     /**
      * Get the model id. Is only relevant for graph abm or hybrid model.
@@ -609,7 +637,7 @@ class Model
             compute_exposure_caches(t, dt);
             m_are_exposure_caches_valid = true;
         }
-        auto personal_rng = PersonalRandomNumberGenerator(person);
+        auto personal_rng = PersonalRandomNumberGenerator(m_rng, person);
         auto location     = person.get_location();
         mio::abm::interact(personal_rng, person, get_location(location),
                            m_local_population_by_age_cache[location.get()], m_air_exposure_rates_cache[location.get()],

cpp/models/abm/person.cpp

@@ -49,8 +49,6 @@ Person::Person(mio::RandomNumberGenerator& rng, LocationType location_type, Loca
     , m_test_results({TestType::Count}, TestResult())
     , m_assigned_location_model_ids((int)LocationType::Count)
     , m_person_id(person_id)
-    , m_rng_key(rng.get_key())
-    , m_rng_index(static_cast<uint32_t>(person_id.get()))
 {
     m_random_workgroup        = UniformDistribution<ScalarType>::get_instance()(rng);
     m_random_schoolgroup      = UniformDistribution<ScalarType>::get_instance()(rng);
@@ -62,7 +60,6 @@ Person::Person(const Person& other, PersonId person_id)
     : Person(other)
 {
     m_person_id = person_id;
-    m_rng_index = static_cast<uint32_t>(person_id.get());
 }
 
 bool Person::is_infected(TimePoint t) const

cpp/models/abm/person.h

@@ -23,7 +23,6 @@
 #include "abm/infection.h"
 #include "abm/infection_state.h"
 #include "abm/location_id.h"
-#include "abm/location.h"
 #include "abm/location_type.h"
 #include "abm/parameters.h"
 #include "abm/person_id.h"
@@ -398,24 +397,6 @@ class Person
         return m_rng_counter;
     }
 
-    /**
-     * @brief Get this Person's index that is used for the RandomNumberGenerator.
-     * @see mio::abm::PersonalRandomNumberGenerator.
-     */
-    uint32_t get_rng_index()
-    {
-        return m_rng_index;
-    }
-
-    /**
-     * @brief Get this Person's key that is used for the RandomNumberGenerator.
-     * @see mio::abm::PersonalRandomNumberGenerator.
-     */
-    mio::Key<uint64_t> get_rng_key()
-    {
-        return m_rng_key;
-    }
-
     /**
      * @brief Get the latest #ProtectionType and its initial TimePoint of the Person.
      * @param[in] t TimePoint to check.
@@ -444,8 +425,7 @@ class Person
             .add("last_transport_mode", m_last_transport_mode)
             .add("rng_counter", m_rng_counter)
             .add("test_results", m_test_results)
-            .add("id", m_person_id)
-            .add("rng_index", m_rng_index);
+            .add("id", m_person_id);
     }
 
     /**
@@ -490,8 +470,6 @@ class Person
     std::vector<int>
         m_assigned_location_model_ids; ///< Vector with model ids of the assigned locations. Only used in graph abm.
     PersonId m_person_id; ///< Unique identifier of a person.
-    mio::Key<uint64_t> m_rng_key; ///< Key for PersonalRandomNumberGenerator
-    uint32_t m_rng_index; ///< Index for PersonalRandomNumberGenerator.
     Counter<uint32_t> m_rng_counter{0}; ///< counter for RandomNumberGenerator.
 };
 

cpp/models/abm/personal_rng.cpp

@@ -35,8 +35,9 @@ PersonalRandomNumberGenerator::PersonalRandomNumberGenerator(mio::Key<uint64_t>
 {
 }
 
-PersonalRandomNumberGenerator::PersonalRandomNumberGenerator(Person& person)
-    : PersonalRandomNumberGenerator(person.get_rng_key(), person.get_rng_index(), person.get_rng_counter())
+PersonalRandomNumberGenerator::PersonalRandomNumberGenerator(const RandomNumberGenerator& model_rng, Person& person)
+    : PersonalRandomNumberGenerator(model_rng.get_key(), static_cast<uint32_t>(person.get_id().get()),
+                                    person.get_rng_counter())
 {
 }
 

cpp/models/abm/personal_rng.h

@@ -22,7 +22,6 @@
 #define MIO_ABM_PERSONAL_RNG_H
 
 #include "memilio/utils/random_number_generator.h"
-#include "models/abm/person_id.h"
 #include <cstdint>
 
 namespace mio
@@ -61,7 +60,7 @@ class PersonalRandomNumberGenerator : public mio::RandomNumberGeneratorBase<Pers
      * Creates a RandomNumberGenerator for a person.
      * @param person Reference to the Person who's counter will be used.
      */
-    PersonalRandomNumberGenerator(Person& person);
+    PersonalRandomNumberGenerator(const RandomNumberGenerator& model_rng, Person& person);
 
     /**
      * @return Get the key.