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;
 
@@ -173,9 +174,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(10);
     // Set the number of simulations to run in the study
-    const size_t num_runs = 3;
+    const size_t num_runs = 10;
 
     // 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
@@ -187,42 +188,74 @@ int main()
     // study.get_rng().seed({12341234, 53456, 63451, 5232576, 84586, 52345});
 
     const std::string result_dir = mio::path_join(mio::base_dir(), "example_results");
+    std::filesystem::remove_all(result_dir);
     if (!mio::create_directory(result_dir)) {
         mio::log_error("Could not create result directory \"{}\".", result_dir);
         return 1;
     }
 
+    const std::string result_directory_standard = mio::path_join(result_dir, "standart_results");
+    if (!mio::create_directory(result_directory_standard)) {
+        mio::log_error("Could not create result directory \"{}\".", result_directory_standard);
+        return 1;
+    }
+
+    const std::string result_dir_detailed = mio::path_join(result_dir, "detailed");
+    if (!mio::create_directory(result_dir_detailed)) {
+        mio::log_error("Could not create result directory \"{}\".", result_dir_detailed);
+        return 1;
+    }
+
     auto ensemble_results = study.run(
         [](auto, auto t0_, auto, size_t) {
             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());
+            auto interpolated_result          = mio::interpolate_simulation_result(sim.get_result());
+            auto interpolated_result_detailed = mio::interpolate_simulation_result(sim.get_result_detailed());
+
             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};
+
+            return std::vector<mio::TimeSeries<double>>{interpolated_result, interpolated_result_detailed};
         });
 
     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")));
+
+        std::vector<std::vector<mio::TimeSeries<double>>> ensemble_new, ensemble_detailed;
+        for (auto& run : ensemble_results) {
+            ensemble_new.push_back({std::move(run[0])});
+            ensemble_detailed.push_back({std::move(run[1])});
+        }
+
+        // Percentiles for aggregated results
+        auto new_p05 = ensemble_percentile(ensemble_new, 0.05);
+        auto new_p25 = ensemble_percentile(ensemble_new, 0.25);
+        auto new_p50 = ensemble_percentile(ensemble_new, 0.50);
+        auto new_p75 = ensemble_percentile(ensemble_new, 0.75);
+        auto new_p95 = ensemble_percentile(ensemble_new, 0.95);
+
+        mio::unused(save_result(new_p05, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p05.h5")));
+        mio::unused(save_result(new_p25, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p25.h5")));
+        mio::unused(save_result(new_p50, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p50.h5")));
+        mio::unused(save_result(new_p75, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p75.h5")));
+        mio::unused(save_result(new_p95, {0}, num_age_groups, mio::path_join(result_directory_standard, "Results_p95.h5")));
+
+        // Percentiles for detailed results
+        auto det_p05 = ensemble_percentile(ensemble_detailed, 0.05);
+        auto det_p25 = ensemble_percentile(ensemble_detailed, 0.25);
+        auto det_p50 = ensemble_percentile(ensemble_detailed, 0.50);
+        auto det_p75 = ensemble_percentile(ensemble_detailed, 0.75);
+        auto det_p95 = ensemble_percentile(ensemble_detailed, 0.95);
+
+        mio::unused(save_result(det_p05, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p05.h5")));
+        mio::unused(save_result(det_p25, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p25.h5")));
+        mio::unused(save_result(det_p50, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p50.h5")));
+        mio::unused(save_result(det_p75, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p75.h5")));
+        mio::unused(save_result(det_p95, {0}, num_age_groups, mio::path_join(result_dir_detailed, "Results_p95.h5")));
     }
 
     mio::mpi::finalize();

cpp/models/abm/common_abm_loggers.h

@@ -189,6 +189,67 @@ struct LogInfectionState : mio::LogAlways {
     }
 };
 
+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.
+     * @param[in] sim The simulation of the abm.
+     * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState.
+     */
+    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();
+        const auto persons   = sim.get_model().get_persons();
+
+
+        for (auto i = size_t(0); i < persons.size(); ++i) {
+            auto& p = persons[i];
+                auto index = (((size_t)(mio::abm::InfectionState::Count)) * ((uint32_t)p.get_age().get())) +
+                             ((uint32_t)p.get_infection_state(curr_time));
+                // PRAGMA_OMP(atomic)
+                sum[index] += 1;
+        }
+        return std::make_pair(curr_time, sum);
+    }
+};
+
+struct LogInfectionPerLocationTypePerAgeGroup : mio::LogAlways {
+    using Type = std::pair<mio::abm::TimePoint, Eigen::VectorXd>;
+    /** 
+     * @brief Log the TimeSeries of the number of Person%s in an #InfectionState.
+     * @param[in] sim The simulation of the abm.
+     * @return A pair of the TimePoint and the TimeSeries of the number of Person%s in an #InfectionState.
+     */
+    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();
+        const auto persons = sim.get_model().get_persons();
+
+
+        for (auto i = size_t(0); i < persons.size(); ++i) {
+            auto& p = persons[i];
+
+
+            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,33 @@ 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.
      */
-    const mio::TimeSeries<double>& get_result() 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 aggregated by infection states.
+     */
+    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.
+        Eigen::Index(InfectionState::Count)}; 
+
+    mio::History<TimeSeriesWriter, LogInfectionPerLocationTypePerAgeGroup> history_detailed{
+        Eigen::Index(LocationType::Count) * this->get_model().parameters.get_num_groups()};
+
+
 };
 
 } // 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-hours(1))
         , m_dt(hours(1))
     {
     }
@@ -85,6 +86,15 @@ class Simulation
         return m_t;
     }
 
+    /**
+     * @brief Get the previous time of the Simulation.
+     */
+    TimePoint get_prev_time() const
+    {
+        return m_t_prev;
+    }
+
+
     /**
      * @brief Get the Model that this Simulation evolves.
      */
@@ -103,11 +113,14 @@ 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 previous  TimePoint of the Simulation.
     TimeSpan m_dt; ///< The length of the time steps.
 };