155 add vaccination model from paper

cpp/CMakeLists.txt

@@ -74,6 +74,7 @@ add_subdirectory(memilio)
 if (MEMILIO_BUILD_MODELS)
     add_subdirectory(models/abm)
     add_subdirectory(models/secir)
+    add_subdirectory(models/ode_secirvvs)
     add_subdirectory(models/seir)
 endif()
 if (MEMILIO_BUILD_EXAMPLES)

cpp/examples/parameter_study_secir.cpp

@@ -19,8 +19,9 @@
 */
 #include "secir/secir_parameters_io.h"
 #include "secir/parameter_space.h"
-#include "secir/parameter_studies.h"
+#include "memilio/compartments/parameter_studies.h"
 #include "memilio/mobility/mobility.h"
+#include "memilio/io/result_io.h"
 
 /**
  * @brief creates xml file with a single run parameter study with std 0 (used to save parameters of individual runs)
@@ -68,8 +69,10 @@ write_single_run_result(const int run,
     std::transform(graph.nodes().begin(), graph.nodes().end(), std::back_inserter(ids), [](auto& node) {
         return node.id;
     });
+    auto num_groups = (int)(size_t)graph.nodes()[0].property.get_simulation().get_model().parameters.get_num_groups();
     BOOST_OUTCOME_TRY(
-        mio::save_result(all_results, ids, mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5"))));
+        mio::save_result(all_results, ids, num_groups,
+                    mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5"))));
 
     return mio::success();
 }
@@ -163,17 +166,20 @@ int main()
 
     //create study
     auto num_runs = size_t(1);
-    mio::ParameterStudy<mio::SecirSimulation<>> parameter_study(model, t0, tmax, 0.2, num_runs);
+    mio::ParameterStudy<mio::SecirSimulation<>> parameter_study(model, t0, tmax, num_runs);
 
     //run study
-    int run     = 0;
-    auto lambda = [&run](auto&& graph) {
+    int run           = 0;
+    auto sample_graph = [](auto&& graph) {
+        return mio::draw_sample(graph);
+    };
+    auto handle_result = [&run](auto&& graph) {
         auto write_result_status = write_single_run_result(run++, graph);
         if (!write_result_status) {
             std::cout << "Error writing result: " << write_result_status.error().formatted_message();
         }
     };
-    parameter_study.run(lambda);
+    parameter_study.run(sample_graph, handle_result);
 
     return 0;
 }

cpp/examples/read_graph.cpp

@@ -18,8 +18,10 @@
 * limitations under the License.
 */
 #include "memilio/mobility/mobility.h"
-#include "secir/secir_parameters_io.h"
 #include "memilio/io/mobility_io.h"
+#include "memilio/compartments/parameter_studies.h"
+#include "secir/parameter_space.h"
+#include "secir/secir_parameters_io.h"
 #include <data_dir.h>
 #include <iostream>
 

cpp/memilio/CMakeLists.txt

@@ -2,6 +2,9 @@ configure_file(config_internal.h.in memilio/config_internal.h)
 
 add_library(memilio
     config.h
+    data/analyze_result.h
+    data/analyze_result.cpp
+    epidemiology/age_group.h
     epidemiology/populations.h
     epidemiology/damping.cpp
     epidemiology/damping.h
@@ -18,13 +21,18 @@ add_library(memilio
     epidemiology/holiday_data_de.ipp
     compartments/compartmentalmodel.h
     compartments/simulation.h
+    compartments/parameter_studies.h
     io/io.h
     io/io.cpp
     io/hdf5_cpp.h
     io/json_serializer.h
     io/json_serializer.cpp
     io/mobility_io.h
     io/mobility_io.cpp
+    io/result_io.h
+    io/result_io.cpp
+    io/epi_data.h
+    io/epi_data.cpp
     math/euler.cpp
     math/euler.h
     math/smoother.h

cpp/memilio/compartments/compartmentalmodel.h

@@ -71,9 +71,10 @@ using has_apply_constraints_member_function = is_expression_valid<details::apply
  * studies
  *
  */
-template <class Pop, class Params>
+template <class Comp, class Pop, class Params>
 struct CompartmentalModel {
 public:
+    using Compartments = Comp;
     using Populations = Pop;
     using ParameterSet = Params;
 

cpp/models/secir/parameter_studies.h → cpp/memilio/compartments/parameter_studies.h

@@ -20,8 +20,6 @@
 #ifndef PARAMETER_STUDIES_H
 #define PARAMETER_STUDIES_H
 
-#include "secir/secir.h"
-#include "secir/parameter_space.h"
 #include "memilio/utils/time_series.h"
 #include "memilio/mobility/mobility.h"
 #include "memilio/compartments/simulation.h"
@@ -89,33 +87,18 @@ class ParameterStudy
         m_graph.add_node(0, model);
     }
 
-    /**
-     * @brief create study for single compartment model with normal distributions.
-     * Sets all parameters to normal distribution with specified relative deviation.
-     * @param params SecirParams object 
-     * @param t0 start time of simulations
-     * @param tmax end time of simulations
-     * @param dev_rel relative deviation of parameters distributions
-     * @param num_runs number of runs in ensemble run
-     */
-    ParameterStudy(typename Simulation::Model const& model, double t0, double tmax, double dev_rel, size_t num_runs)
-        : ParameterStudy(model, t0, tmax, num_runs)
-    {
-        set_params_distributions_normal(m_graph.nodes()[0].property, t0, tmax, dev_rel);
-    }
-
     /*
      * @brief Carry out all simulations in the parameter study.
      * Save memory and enable more runs by immediately processing and/or discarding the result.
      * @param result_processing_function Processing function for simulation results, e.g., output function.
      *                                   Receives the result after each run is completed.
      */
-    template<class HandleSimulationResultFunction>
-    void run(HandleSimulationResultFunction result_processing_function)
+    template <class SampleGraphFunction, class HandleSimulationResultFunction>
+    void run(SampleGraphFunction sample_graph, HandleSimulationResultFunction result_processing_function)
     {
         // Iterate over all parameters in the parameter space
         for (size_t i = 0; i < m_num_runs; i++) {
-            auto sim = create_sampled_simulation();
+            auto sim = create_sampled_simulation(sample_graph);
             sim.advance(m_tmax);
 
             result_processing_function(std::move(sim).get_graph());
@@ -127,12 +110,13 @@ class ParameterStudy
      * Convinience function for a few number of runs, but uses a lot of memory.
      * @return vector of results of each run.
      */
-    std::vector<mio::Graph<mio::SimulationNode<Simulation>, mio::MigrationEdge>> run()
+    template <class SampleGraphFunction>
+    std::vector<mio::Graph<mio::SimulationNode<Simulation>, mio::MigrationEdge>> run(SampleGraphFunction sample_graph)
     {
         std::vector<mio::Graph<mio::SimulationNode<Simulation>, mio::MigrationEdge>> ensemble_result;
         ensemble_result.reserve(m_num_runs);
 
-        run([&ensemble_result](auto&& r) {
+        run(sample_graph, [&ensemble_result](auto&& r) {
             ensemble_result.emplace_back(std::move(r));
         });
 
@@ -219,47 +203,18 @@ class ParameterStudy
 
 private:
     //sample parameters and create simulation
-    mio::GraphSimulation<mio::Graph<mio::SimulationNode<Simulation>, mio::MigrationEdge>> create_sampled_simulation()
+    template <class SampleGraphFunction>
+    mio::GraphSimulation<mio::Graph<mio::SimulationNode<Simulation>, mio::MigrationEdge>>
+    create_sampled_simulation(SampleGraphFunction sample_graph)
     {
         mio::Graph<mio::SimulationNode<Simulation>, mio::MigrationEdge> sim_graph;
 
-        //sample global parameters
-        auto& shared_params_model = m_graph.nodes()[0].property;
-        draw_sample_infection(shared_params_model);
-        auto& shared_contacts = shared_params_model.parameters.template get<mio::ContactPatterns>();
-        shared_contacts.draw_sample_dampings();
-        auto& shared_dynamic_npis = shared_params_model.parameters.template get<DynamicNPIsInfected>();
-        shared_dynamic_npis.draw_sample();
-
-        for (auto& params_node : m_graph.nodes()) {
-            auto& node_model = params_node.property;
-
-            //sample local parameters
-            draw_sample_demographics(params_node.property);
-
-            //copy global parameters
-            //save demographic parameters so they aren't overwritten
-            auto local_icu_capacity = node_model.parameters.template get<ICUCapacity>();
-            auto local_tnt_capacity = node_model.parameters.template get<TestAndTraceCapacity>();
-            auto local_holidays = node_model.parameters.template get<ContactPatterns>().get_school_holidays();
-            node_model.parameters = shared_params_model.parameters;
-            node_model.parameters.template get<ICUCapacity>() = local_icu_capacity;
-            node_model.parameters.template get<TestAndTraceCapacity>() = local_tnt_capacity;
-            node_model.parameters.template get<ContactPatterns>().get_school_holidays() = local_holidays;
-
-            node_model.parameters.template get<ContactPatterns>().make_matrix();
-            node_model.apply_constraints();
-
-            sim_graph.add_node(params_node.id, node_model, m_t0, m_dt_integration);
+        auto sampled_graph = sample_graph(m_graph);
+        for (auto&& node : sampled_graph.nodes()) {
+            sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration);
         }
-
-        for (auto& edge : m_graph.edges()) {
-            auto edge_params = edge.property;
-            apply_dampings(edge_params.get_coefficients(), shared_contacts.get_dampings(), [&edge_params](auto& v) {
-                return make_migration_damping_vector(edge_params.get_coefficients().get_shape(), v);
-            });
-            edge_params.set_dynamic_npis_infected(shared_dynamic_npis);
-            sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge_params);
+        for (auto&& edge : sampled_graph.edges()) {
+            sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property);
         }
 
         return make_migration_sim(m_t0, m_dt_graph_sim, std::move(sim_graph));

cpp/memilio/data/analyze_result.cpp

@@ -0,0 +1,165 @@
+/* 
+* Copyright (C) 2020-2021 German Aerospace Center (DLR-SC)
+*
+* Authors: Wadim Koslow, Daniel Abele
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+#include "memilio/data/analyze_result.h"
+
+#include <algorithm>
+#include <cassert>
+
+namespace mio
+{
+
+TimeSeries<double> interpolate_simulation_result(const TimeSeries<double>& simulation_result)
+{
+    assert(simulation_result.get_num_time_points() > 0 && "TimeSeries must not be empty.");
+
+    const auto t0      = simulation_result.get_time(0);
+    const auto tmax    = simulation_result.get_last_time();
+    const auto day0    = static_cast<int>(floor(t0));
+    const auto day_max = static_cast<int>(ceil(tmax));
+
+    auto day = day0;
+    TimeSeries<double> interpolated(simulation_result.get_num_elements());
+    interpolated.reserve(day_max - day0 + 1);
+    interpolated.add_time_point(day, simulation_result.get_value(0));
+    day++;
+
+    //interpolate between pair of time points that lie on either side of each integer day
+    //TODO: extrapolate first and last point
+    for (int i = 0; i < simulation_result.get_num_time_points() - 1;) {
+        //only go to next pair of time points if no time point is added.
+        //otherwise check the same time points again
+        //in case there is more than one day between the two time points
+        if (simulation_result.get_time(i) < day && simulation_result.get_time(i + 1) >= day) {
+            auto weight = (day - simulation_result.get_time(i)) /
+                          (simulation_result.get_time(i + 1) - simulation_result.get_time(i));
+            interpolated.add_time_point(day, simulation_result[i] +
+                                                 (simulation_result[i + 1] - simulation_result[i]) * weight);
+            ++day;
+        }
+        else {
+            ++i;
+        }
+    }
+
+    if (day_max > tmax) {
+        interpolated.add_time_point(day, simulation_result.get_last_value());
+    }
+
+    return interpolated;
+}
+
+std::vector<std::vector<TimeSeries<double>>>
+sum_nodes(const std::vector<std::vector<TimeSeries<double>>>& ensemble_result)
+{
+    auto num_runs        = ensemble_result.size();
+    auto num_nodes       = ensemble_result[0].size();
+    auto num_time_points = ensemble_result[0][0].get_num_time_points();
+    auto num_elements    = ensemble_result[0][0].get_num_elements();
+
+    std::vector<std::vector<TimeSeries<double>>> sum_result(
+        num_runs, std::vector<TimeSeries<double>>(1, TimeSeries<double>::zero(num_time_points, num_elements)));
+
+    for (size_t run = 0; run < num_runs; run++) {
+        for (Eigen::Index time = 0; time < num_time_points; time++) {
+            sum_result[run][0].get_time(time) = ensemble_result[run][0].get_time(time);
+            for (size_t node = 0; node < num_nodes; node++) {
+                sum_result[run][0][time] += ensemble_result[run][node][time];
+            }
+        }
+    }
+    return sum_result;
+}
+
+std::vector<TimeSeries<double>> ensemble_mean(const std::vector<std::vector<TimeSeries<double>>>& ensemble_result)
+{
+    auto num_runs        = ensemble_result.size();
+    auto num_nodes       = ensemble_result[0].size();
+    auto num_time_points = ensemble_result[0][0].get_num_time_points();
+    auto num_elements    = ensemble_result[0][0].get_num_elements();
+
+    std::vector<TimeSeries<double>> mean(num_nodes, TimeSeries<double>::zero(num_time_points, num_elements));
+
+    for (size_t run = 0; run < num_runs; run++) {
+        assert(ensemble_result[run].size() == num_nodes && "ensemble results not uniform.");
+        for (size_t node = 0; node < num_nodes; node++) {
+            assert(ensemble_result[run][node].get_num_time_points() == num_time_points &&
+                   "ensemble results not uniform.");
+            for (Eigen::Index time = 0; time < num_time_points; time++) {
+                assert(ensemble_result[run][node].get_num_elements() == num_elements &&
+                       "ensemble results not uniform.");
+                mean[node].get_time(time) = ensemble_result[run][node].get_time(time);
+                mean[node][time] += ensemble_result[run][node][time] / num_runs;
+            }
+        }
+    }
+
+    return mean;
+}
+
+std::vector<TimeSeries<double>> ensemble_percentile(const std::vector<std::vector<TimeSeries<double>>>& ensemble_result,
+                                                    double p)
+{
+    assert(p > 0.0 && p < 1.0 && "Invalid percentile value.");
+
+    auto num_runs        = ensemble_result.size();
+    auto num_nodes       = ensemble_result[0].size();
+    auto num_time_points = ensemble_result[0][0].get_num_time_points();
+    auto num_elements    = ensemble_result[0][0].get_num_elements();
+
+    std::vector<TimeSeries<double>> percentile(num_nodes, TimeSeries<double>::zero(num_time_points, num_elements));
+
+    std::vector<double> single_element_ensemble(num_runs); //reused for each element
+    for (size_t node = 0; node < num_nodes; node++) {
+        for (Eigen::Index time = 0; time < num_time_points; time++) {
+            percentile[node].get_time(time) = ensemble_result[0][node].get_time(time);
+            for (Eigen::Index elem = 0; elem < num_elements; elem++) {
+                std::transform(ensemble_result.begin(), ensemble_result.end(), single_element_ensemble.begin(),
+                               [=](auto& run) {
+                                   return run[node][time][elem];
+                               });
+                std::sort(single_element_ensemble.begin(), single_element_ensemble.end());
+                percentile[node][time][elem] = single_element_ensemble[static_cast<size_t>(num_runs * p)];
+            }
+        }
+    }
+    return percentile;
+}
+
+double result_distance_2norm(const std::vector<mio::TimeSeries<double>>& result1,
+                             const std::vector<mio::TimeSeries<double>>& result2)
+{
+    assert(result1.size() == result2.size());
+    assert(result1.size() > 0);
+    assert(result1[0].get_num_time_points() > 0);
+    assert(result1[0].get_num_elements() > 0);
+
+    auto norm_sqr = 0.0;
+    for (auto iter_node1 = result1.begin(), iter_node2 = result2.begin(); iter_node1 < result1.end();
+         ++iter_node1, ++iter_node2) {
+        for (Eigen::Index time_idx = 0; time_idx < iter_node1->get_num_time_points(); ++time_idx) {
+            auto v1 = (*iter_node1)[time_idx];
+            auto v2 = (*iter_node2)[time_idx];
+            norm_sqr += ((v1 - v2).array() * (v1 - v2).array()).sum();
+        }
+    }
+    return std::sqrt(norm_sqr);
+}
+
+} // namespace mio