1210 abm paper logger

cpp/examples/abm_parameter_study.cpp

@@ -34,6 +34,7 @@
 #include "memilio/utils/stl_util.h"
 
 #include <string>
+#include <filesystem>
 
 constexpr size_t num_age_groups = 4;
 
@@ -172,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(30);
     // Set the number of simulations to run in the study
-    const size_t num_runs = 3;
+    const size_t num_runs = 11;
     // Set up an RNG
     mio::RandomNumberGenerator rng;
 
@@ -194,6 +195,14 @@ int main()
 
     const auto result_dir = mio::create_directories_or_exit(mio::example_results_dir("abm_parameter_study"));
 
+    // Each run writes its aggregated and its detailed (per location type and age group) result into its own
+    // subdirectory of the result directory.
+    const auto result_dir_standard = mio::create_directories_or_exit(result_dir / "standard_results");
+    const auto result_dir_detailed = mio::create_directories_or_exit(result_dir / "detailed_results");
+
+    // Collects the detailed result of each run, see the process_simulation_result lambda below.
+    std::vector<std::vector<mio::TimeSeries<ScalarType>>> ensemble_results_detailed;
+
     // 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!
@@ -206,29 +215,61 @@ int main()
             copy.reset_rng(ctr);
             return mio::abm::ResultSimulation(std::move(copy), t0_);
         },
-        [&result_dir](auto&& sim, auto&& run_idx) {
-            auto interpolated_result = mio::interpolate_simulation_result(sim.get_result());
-            auto outpath             = result_dir / ("abm_minimal_run_" + std::to_string(run_idx) + ".txt");
+        [&result_dir_standard, &result_dir_detailed, &ensemble_results_detailed](auto&& sim, auto&& run_idx) {
+            auto interpolated_result          = mio::interpolate_simulation_result(sim.get_result());
+            auto interpolated_result_detailed = mio::interpolate_simulation_result(sim.get_result_detailed());
+
+            const auto outpath = result_dir_standard / ("abm_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.string() << std::endl;
+
+            // The detailed result has one column per (location type, age group) pair, so we let print_table
+            // generate the column names itself.
+            const auto outpath_detailed = result_dir_detailed / ("abm_run_" + std::to_string(run_idx) + ".txt");
+            std::ofstream outfile_run_detailed(outpath_detailed);
+            sim.get_result_detailed().print_table(outfile_run_detailed, {}, 7, 4);
+            std::cout << "Detailed results written to " << outpath_detailed.string() << std::endl;
+
+            // The detailed result is collected separately, because ensemble_percentile requires all entries of an
+            // ensemble to have the same number of elements, which the two results do not have.
+            ensemble_results_detailed.push_back({interpolated_result_detailed});
+
             return std::vector{interpolated_result};
         });
 
     // 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(result_dir / fmt::format("Results_p{:0<4.2}.txt", p));
-            auto ensemble_percentiles = ensemble_percentile(ensemble_results, p);
+        // The percentiles are written in the "Results_p05.h5" layout expected by
+        // pycode/memilio-plot/memilio/plot/plotAbmInfectionStates.py, so that the results of this example can be
+        // plotted directly. The percentile is given in whole percent to match that file name.
+        const auto write_percentile = [&](int p) {
+            auto ensemble_percentiles          = ensemble_percentile(ensemble_results, p / 100.0);
+            auto ensemble_percentiles_detailed = ensemble_percentile(ensemble_results_detailed, p / 100.0);
+
+            std::ofstream out(result_dir_standard / fmt::format("Results_p{:02d}.txt", p));
             ensemble_percentiles.front().print_table(out, {"S", "E", "I_NS", "I_Sy", "I_Sev", "I_Crit", "R", "D"}, 7,
                                                      4);
+
+            // save_result splits each row into num_groups groups. Both results are a single vector per time point,
+            // so they are written as one group each, whose "Total" holds all entries of that vector. Passing the
+            // number of age groups here would instead split those entries across that many groups.
+            mio::unused(mio::save_result(ensemble_percentiles, {0}, 1,
+                                         (result_dir_standard / fmt::format("Results_p{:02d}.h5", p)).string()));
+            mio::unused(mio::save_result(ensemble_percentiles_detailed, {0}, 1,
+                                         (result_dir_detailed / fmt::format("Results_p{:02d}.h5", p)).string()));
         };
 
-        write_percentile(0.05);
-        write_percentile(0.25);
-        write_percentile(0.50);
-        write_percentile(0.75);
-        write_percentile(0.95);
+        write_percentile(5);
+        write_percentile(25);
+        write_percentile(50);
+        write_percentile(75);
+        write_percentile(95);
+
+        std::cout << "\nPercentiles written. Plot them with (in folder memilio-plot) &&:\n"
+                  << "   python -m memilio.plot.plotAbmInfectionStates \\\n"
+                  << "      --path-to-infection-states " << result_dir_standard.string() << " \\\n"
+                  << "      --path-to-loc-types " << result_dir_detailed.string() << std::endl;
     }
 
     mio::mpi::finalize();

cpp/models/abm/common_abm_loggers.h

@@ -37,10 +37,10 @@ namespace abm
 
 /**
  * @brief Struct to save specific mobility data of an agent.
- * The data consists of:
- * 
+ * The data consists of the agent's id, the Location%s it moved between, the times it left and arrived,
+ * the TransportMode and ActivityType of the trip, and the agent's InfectionState.
  */
-struct mobility_data {
+struct MobilityData {
     uint32_t agent_id;
     uint32_t from_id;
     uint32_t to_id;
@@ -51,6 +51,11 @@ struct mobility_data {
     mio::abm::InfectionState infection_state;
 };
 
+/**
+ * @brief Deduce the ActivityType an agent is most likely pursuing from the LocationType it is at.
+ * @param[in] current_location The type of the Location the agent is currently at.
+ * @return The ActivityType associated with the given LocationType, or ActivityType::UnknownActivity if there is none.
+ */
 constexpr mio::abm::ActivityType guess_activity_type(mio::abm::LocationType current_location)
 {
     switch (current_location) {
@@ -153,7 +158,7 @@ struct LogDataForMobility : mio::LogAlways {
     static Type log(const mio::abm::Simulation<>& sim)
     {
         Type mobility_data{};
-        for (Person p : sim.get_model().get_persons()) {
+        for (const Person& p : sim.get_model().get_persons()) {
             mobility_data.push_back(
                 std::make_tuple(p.get_id(), p.get_location(), sim.get_time(), p.get_last_transport_mode(),
                                 guess_activity_type(p.get_location_type()), p.get_infection_state(sim.get_time())));
@@ -174,11 +179,9 @@ struct LogInfectionState : mio::LogAlways {
      */
     static Type log(const mio::abm::Simulation<>& sim)
     {
-
         Eigen::VectorX<ScalarType> sum =
             Eigen::VectorX<ScalarType>::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++) {
                 sum[inf_state] += sim.get_model().get_subpopulation(location.get_id(), curr_time,
@@ -189,6 +192,75 @@ struct LogInfectionState : mio::LogAlways {
     }
 };
 
+/**
+* @brief Logger to log the TimeSeries of the number of Person%s in an #InfectionState per AgeGroup.
+* This is a finer grained variant of LogInfectionState. It is not used by ResultSimulation, but provided here so that
+* a Simulation can be given an age resolved History where that resolution is needed.
+*/
+struct LogInfectionStatePerAgeGroup : mio::LogAlways {
+    using Type = std::pair<mio::abm::TimePoint, Eigen::VectorXd>;
+    /**
+     * @brief Log the TimeSeries of the number of Person%s in an #InfectionState per AgeGroup.
+     * @param[in] sim The simulation of the abm.
+     * @return A pair of the TimePoint and a vector counting the Person%s per AgeGroup and #InfectionState,
+     * indexed by `age_group * InfectionState::Count + infection_state`.
+     */
+    static Type log(const mio::abm::Simulation<>& sim)
+    {
+        Eigen::VectorXd sum = Eigen::VectorXd::Zero(
+            Eigen::Index((size_t)mio::abm::InfectionState::Count * sim.get_model().parameters.get_num_groups()));
+        const auto curr_time = sim.get_time();
+
+        for (const Person& p : sim.get_model().get_persons()) {
+            auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) +
+                         ((uint32_t)p.get_infection_state(curr_time));
+            sum[index] += 1;
+        }
+        return std::make_pair(curr_time, sum);
+    }
+};
+
+/**
+* @brief Logger to log the TimeSeries of new #Infection%s per LocationType and AgeGroup.
+*/
+struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways {
+    using Type = std::pair<mio::abm::TimePoint, Eigen::VectorXd>;
+    /**
+     * @brief Log the TimeSeries of new #Infection%s per LocationType and AgeGroup.
+     * A Person is counted if it became #InfectionState::Exposed since the previous time step. It is attributed to the
+     * LocationType it is at when the transition is observed. At the first time step there is no previous time step to
+     * compare against, so no Person is counted. This avoids miscounting the initially infected Person%s of the
+     * Model as new #Infection%s.
+     * @param[in] sim The simulation of the abm.
+     * @return A pair of the TimePoint and a vector counting the newly exposed Person%s per AgeGroup and LocationType,
+     * indexed by `age_group * LocationType::Count + location_type`.
+     */
+    static Type log(const mio::abm::Simulation<>& sim)
+    {
+        Eigen::VectorXd sum = Eigen::VectorXd::Zero(
+            Eigen::Index((size_t)mio::abm::LocationType::Count * sim.get_model().parameters.get_num_groups()));
+        auto curr_time = sim.get_time();
+        auto prev_time = sim.get_prev_time();
+
+        // Before the first time step is evolved there is no previous state to compare against. Persons that are
+        // infected at initialization must not be counted as new Infections.
+        if (prev_time >= curr_time) {
+            return std::make_pair(curr_time, sum);
+        }
+
+        for (const Person& p : sim.get_model().get_persons()) {
+            if ((p.get_infection_state(prev_time) != mio::abm::InfectionState::Exposed) &&
+                (p.get_infection_state(curr_time) == mio::abm::InfectionState::Exposed)) {
+                auto index = (((size_t)(mio::abm::LocationType::Count)) * ((uint32_t)p.get_age().get())) +
+                             ((uint32_t)p.get_location_type());
+                sum[index] += 1;
+            }
+        }
+        return std::make_pair(curr_time, sum);
+    }
+};
+
+
 /**
 * @brief This is like the DataWriterToMemory, but it only logs time series data.
 * @tparam Loggers The loggers that are used to log data. The loggers must return a touple with a TimePoint and a value.

cpp/models/abm/result_simulation.h

@@ -45,19 +45,31 @@ class ResultSimulation : public Simulation<M>
      */
     void advance(TimePoint tmax)
     {
-        Simulation<Model>::advance(tmax, history);
+        Simulation<Model>::advance(tmax, history, history_detailed);
     }
 
     /**
-     * @brief Return the simulation result aggregated by infection states.
+     * @brief Return the simulation result aggregated by #InfectionState.
      */
     const mio::TimeSeries<double>& get_result() const
     {
-        return get<0>(history.get_log());
+        return std::get<0>(history.get_log());
+    }
+
+    /**
+     * @brief Return the detailed simulation result, i.e. the new #Infection%s aggregated by LocationType and AgeGroup.
+     */
+    const mio::TimeSeries<double>& get_result_detailed() const
+    {
+        return std::get<0>(history_detailed.get_log());
     }
 
     mio::History<TimeSeriesWriter, LogInfectionState> history{
         Eigen::Index(InfectionState::Count)}; ///< History used to create the result TimeSeries.
+
+    mio::History<TimeSeriesWriter, LogInfectionPerLocationTypePerAgeGroup> history_detailed{
+        Eigen::Index(LocationType::Count) *
+        this->get_model().parameters.get_num_groups()}; ///< History used to create the detailed result TimeSeries.
 };
 
 } // namespace abm

cpp/models/abm/simulation.h

@@ -46,6 +46,7 @@ class Simulation
     Simulation(TimePoint t0, Model&& model)
         : m_model(std::move(model))
         , m_t(t0)
+        , m_t_prev(t0)
         , m_dt(hours(1))
     {
     }
@@ -85,6 +86,15 @@ class Simulation
         return m_t;
     }
 
+    /**
+     * @brief Get the time of the Simulation before the last evolve step.
+     * Equal to get_time() as long as no time step has been evolved yet.
+     */
+    TimePoint get_prev_time() const
+    {
+        return m_t_prev;
+    }
+
     /**
      * @brief Get the Model that this Simulation evolves.
      */
@@ -103,11 +113,13 @@ class Simulation
     {
         auto dt = std::min(m_dt, tmax - m_t);
         m_model.evolve(m_t, dt);
+        m_t_prev = m_t;
         m_t += m_dt;
     }
 
     Model m_model; ///< The Model to simulate.
     TimePoint m_t; ///< The current TimePoint of the Simulation.
+    TimePoint m_t_prev; ///< The TimePoint of the Simulation before the last evolve step.
     TimeSpan m_dt; ///< The length of the time steps.
 };
 

cpp/tests/test_abm_mobility_rules.cpp

@@ -495,8 +495,9 @@ TEST_F(TestMobilityRules, hospital)
     EXPECT_EQ(mio::abm::go_to_hospital(rng_inf, p_inf, t, dt, mio::abm::Parameters(num_age_groups)),
               mio::abm::LocationType::Hospital);
 
+    // The infection has to start at t, otherwise it may already have progressed to InfectedSevere by then.
     auto p_car =
-        make_test_person(this->get_rng(), home, age_group_15_to_34, mio::abm::InfectionState::InfectedSymptoms);
+        make_test_person(this->get_rng(), home, age_group_15_to_34, mio::abm::InfectionState::InfectedSymptoms, t);
     auto rng_car = mio::abm::PersonalRandomNumberGenerator(this->get_rng(), p_car);
     // Ensure person has infection symptoms still stay at home
     EXPECT_EQ(mio::abm::go_to_hospital(rng_car, p_car, t, dt, mio::abm::Parameters(num_age_groups)),