1210 abm paper logger

cpp/examples/abm_parameter_study.cpp

@@ -17,10 +17,11 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#include "abm/result_simulation.h"
+#include "abm/common_abm_loggers.h"
 #include "abm/household.h"
 #include "abm/lockdown_rules.h"
 #include "abm/model.h"
+#include "abm/simulation.h"
 #include "abm/time.h"
 
 #include "memilio/compartments/parameter_studies.h"
@@ -34,6 +35,7 @@
 #include "memilio/utils/stl_util.h"
 
 #include <string>
+#include <filesystem>
 
 constexpr size_t num_age_groups = 4;
 
@@ -164,6 +166,44 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng)
     return model;
 }
 
+/**
+ * @brief Simulation for this study, holding the Histories needed to evaluate a run.
+ * ParameterStudy only calls advance(tmax), so the Simulation has to own its Histories. Besides the result aggregated by
+ * #InfectionState, it logs the new Infections per LocationType and AgeGroup.
+ */
+class StudySimulation : public mio::abm::Simulation<>
+{
+public:
+    StudySimulation(mio::abm::Model&& model, mio::abm::TimePoint t0)
+        : Simulation(t0, std::move(model))
+        , history_detailed{Eigen::Index((size_t)mio::abm::LocationType::Count *
+                                        get_model().parameters.get_num_groups())}
+    {
+    }
+
+    void advance(mio::abm::TimePoint tmax)
+    {
+        Simulation::advance(tmax, history, history_detailed);
+    }
+
+    /// @brief The result aggregated by #InfectionState.
+    const mio::TimeSeries<ScalarType>& get_result() const
+    {
+        return std::get<0>(history.get_log());
+    }
+
+    /// @brief The new Infections per LocationType and AgeGroup.
+    const mio::TimeSeries<ScalarType>& get_result_detailed() const
+    {
+        return std::get<0>(history_detailed.get_log());
+    }
+
+private:
+    mio::History<mio::abm::TimeSeriesWriter, mio::abm::LogInfectionState> history{
+        Eigen::Index(mio::abm::InfectionState::Count)};
+    mio::History<mio::abm::TimeSeriesWriter, mio::abm::LogInfectionPerLocationTypePerAgeGroup> history_detailed;
+};
+
 int main()
 {
     mio::mpi::init();
@@ -172,9 +212,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 +234,11 @@ 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");
+
     // 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!
@@ -204,31 +249,76 @@ int main()
             // 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_);
+            return StudySimulation(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](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;
-            return std::vector{interpolated_result};
+
+            // 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;
+
+            // Both results are returned, so that the study gathers them from all MPI ranks on the root rank.
+            return std::vector{interpolated_result, interpolated_result_detailed};
         });
 
     // 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);
+        // Split the gathered results into one ensemble per result, because ensemble_percentile requires all entries
+        // of an ensemble to have the same number of elements, which the two results do not have.
+        std::vector<std::vector<mio::TimeSeries<ScalarType>>> ensemble_results_standard, ensemble_results_detailed;
+        for (auto& run_result : ensemble_results) {
+            ensemble_results_standard.push_back({std::move(run_result[0])});
+            ensemble_results_detailed.push_back({std::move(run_result[1])});
+        }
+
+        // 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_standard, 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.
+            const auto save_status = mio::save_result(
+                ensemble_percentiles, {0}, 1, (result_dir_standard / fmt::format("Results_p{:02d}.h5", p)).string());
+            if (!save_status) {
+                std::cout << "Error writing result: " << save_status.error().formatted_message() << std::endl;
+            }
+            const auto save_status_detailed =
+                mio::save_result(ensemble_percentiles_detailed, {0}, 1,
+                                 (result_dir_detailed / fmt::format("Results_p{:02d}.h5", p)).string());
+            if (!save_status_detailed) {
+                std::cout << "Error writing detailed result: " << save_status_detailed.error().formatted_message()
+                          << std::endl;
+            }
         };
 
-        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

@@ -30,17 +30,19 @@
 #include "abm/mobility_data.h"
 #include "memilio/utils/mioomp.h"
 
+#include <vector>
+
 namespace mio
 {
 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 +53,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 +160,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,19 +181,95 @@ 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,
-                                                                    mio::abm::InfectionState(inf_state));
+        const auto curr_time = sim.get_time();
+        const auto& model    = sim.get_model();
+        for (const Person& p : model.get_persons()) {
+            if (p.get_location_model_id() == model.get_id()) {
+                sum[(Eigen::Index)p.get_infection_state(curr_time)] += 1;
+            }
+        }
+        return std::make_pair(curr_time, sum);
+    }
+};
+
+/**
+* @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.
+* Unlike the other loggers, this one keeps a state: the time and the LocationType of each Person at its previous log.
+* Hence, it has to log after every time step, which is the case when used in a History passed to Simulation::advance.
+*/
+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 log. A Person is infected during the
+     * interaction at the beginning of a time step and may move to another Location afterwards. Hence, it is attributed
+     * to the LocationType it was at during the previous log, not to its current one.
+     * At the first log there is no previous log 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`.
+     */
+    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()));
+        const auto curr_time     = sim.get_time();
+        const auto persons       = sim.get_model().get_persons();
+        const size_t num_persons = persons.size();
+        // Only Persons that were already present at the previous log can be compared against it.
+        const size_t num_previous = m_prev_location_types.size();
+        m_prev_location_types.resize(num_persons);
+
+        size_t i = 0;
+        for (const Person& p : persons) {
+            if (i < num_previous && (p.get_infection_state(m_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)m_prev_location_types[i]);
+                sum[index] += 1;
             }
+            // Remember where the Person is now, which is where the interaction of the next time step happens.
+            m_prev_location_types[i] = p.get_location_type();
+            ++i;
         }
+        m_prev_time = curr_time;
         return std::make_pair(curr_time, sum);
     }
+
+private:
+    mio::abm::TimePoint m_prev_time{0}; ///< Time of the previous log.
+    std::vector<mio::abm::LocationType> m_prev_location_types; ///< LocationType of each Person at the previous log.
 };
 
 /**

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)),