385 traveltime aware mobility scheme

cpp/examples/CMakeLists.txt

@@ -76,6 +76,10 @@ add_executable(ode_secirvvs_example ode_secirvvs.cpp)
 target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs)
 target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ode_seir_detailed_mobility ode_seir_detailed_mobility.cpp)
+target_link_libraries(ode_seir_detailed_mobility PRIVATE memilio ode_seir)
+target_compile_options(ode_seir_detailed_mobility PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 add_executable(ode_secirts_example ode_secirts.cpp)
 target_link_libraries(ode_secirts_example PRIVATE memilio ode_secirts)
 target_compile_options(ode_secirts_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/examples/ode_seir_detailed_mobility.cpp

@@ -0,0 +1,101 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Henrik Zunker
+*
+* 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 "ode_seir/model.h"
+#include "ode_seir/infection_state.h"
+#include "ode_seir/parameters.h"
+#include "memilio/mobility/metapopulation_mobility_instant.h"
+#include "memilio/mobility/metapopulation_mobility_detailed.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/compartments/flow_simulation.h"
+#include "memilio/data/analyze_result.h"
+
+int main()
+{
+    const auto t0   = 0.;
+    const auto tmax = 10.;
+    const auto dt   = 0.5; //time step of migration, daily migration every second step
+
+    mio::oseir::Model<> model(1);
+
+    // set population
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000;
+
+    // set transition times
+    model.parameters.set<mio::oseir::TimeExposed<>>(1);
+    model.parameters.set<mio::oseir::TimeInfected<>>(1);
+
+    // set contact matrix
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns<>>().get_cont_freq_mat();
+    contact_matrix[0].get_baseline().setConstant(2.7);
+
+    //two mostly identical groups
+    auto model_group1 = model;
+    auto model_group2 = model;
+
+    //some contact restrictions in group 1
+    mio::ContactMatrixGroup& contact_matrix1 =
+        model_group1.parameters.get<mio::oseir::ContactPatterns<>>().get_cont_freq_mat();
+    contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5));
+
+    //infection starts in group 1
+    model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990;
+    model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}]     = 10;
+
+    auto sim1          = mio::FlowSimulation<ScalarType, mio::oseir::Model<>>(model_group1, t0, dt);
+    auto sim2          = mio::FlowSimulation<ScalarType, mio::oseir::Model<>>(model_group2, t0, dt);
+    double stay_time_1 = 0.3;
+    double stay_time_2 = 0.3;
+    // mio::ExtendedGraph<mio::Simulation<ScalarType, mio::oseir::Model<>>> g;
+    mio::ExtendedGraph<mio::FlowSimulation<ScalarType, mio::oseir::Model<double>>> g;
+    g.add_node(1, sim1, sim2, stay_time_1);
+    g.add_node(2, sim1, sim2, stay_time_2);
+
+    double traveltime      = 0.1;
+    std::vector<int> path1 = {0, 1};
+    std::vector<int> path2 = {1, 0};
+    g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01), traveltime, path1);
+    g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01), traveltime, path2);
+
+    auto sim = mio::make_mobility_sim(t0, dt, std::move(g));
+
+    sim.advance(tmax);
+    // results node 1
+    std::cout << "Results node 1" << std::endl;
+    auto interpolated_sim1 =
+        mio::interpolate_simulation_result(sim.get_graph().nodes()[0].property.base_sim.get_result());
+    interpolated_sim1.print_table({"S", "E", "I", "R"});
+
+    // results node 1 mobility_sim
+    std::cout << "Mobility results node 1" << std::endl;
+    auto interpolated_sim1_mobility = sim.get_graph().nodes()[0].property.mobility_sim.get_result();
+    interpolated_sim1_mobility.print_table({"S", "E", "I", "R"});
+
+    // results node 2
+    std::cout << "Results node 2" << std::endl;
+    auto interpolated_sim2 = sim.get_graph().nodes()[1].property.base_sim.get_result();
+    interpolated_sim2.print_table({"S", "E", "I", "R"});
+
+    // results node 2 mobility_sim
+    std::cout << "Mobility results node 2" << std::endl;
+    auto interpolated_sim2_mobility = sim.get_graph().nodes()[1].property.mobility_sim.get_result();
+    interpolated_sim2_mobility.print_table({"S", "E", "I", "R"});
+
+    return 0;
+}

cpp/memilio/CMakeLists.txt

@@ -66,6 +66,8 @@ add_library(memilio
     mobility/metapopulation_mobility_instant.cpp
     mobility/metapopulation_mobility_stochastic.h
     mobility/metapopulation_mobility_stochastic.cpp
+    mobility/metapopulation_mobility_detailed.cpp
+    mobility/metapopulation_mobility_detailed.h
     mobility/graph_simulation.h
     mobility/graph_simulation.cpp
     mobility/graph.h

cpp/memilio/compartments/flow_model.h

@@ -55,8 +55,8 @@ using filtered_tuple_t = decltype(filter_tuple<OmittedTag>(std::declval<Tuple>()
 
 // Remove all occurrences of OmittedTag from the types in an Index = IndexTemplate<types...>.
 template <class OmittedTag, template <class...> class IndexTemplate, class Index>
-using filtered_index_t = decltype(as_index<IndexTemplate>(
-    std::declval<filtered_tuple_t<OmittedTag, decltype(as_tuple(std::declval<Index>()))>>()));
+using filtered_index_t = decltype(
+    as_index<IndexTemplate>(std::declval<filtered_tuple_t<OmittedTag, decltype(as_tuple(std::declval<Index>()))>>()));
 
 } //namespace details
 
@@ -203,6 +203,26 @@ class FlowModel : public CompartmentalModel<FP, Comp, Pop, Params>
         return index_of_type_v<Flow<Source, Target>, Flows>;
     }
 
+    /**
+     * @brief Returns the current flow values.
+     *
+     * @return A constant reference to an Eigen::VectorX containing the current flow values.
+     */
+    Eigen::VectorX<FP>& get_flow_values() const
+    {
+        return m_flow_values;
+    }
+
+    /**
+     * @brief Sets the flow values.
+     *
+     * @param flows A constant reference to an Eigen::VectorX containing flow values.
+     */
+    void set_flow_values(const Eigen::VectorX<FP> flows)
+    {
+        m_flow_values = flows;
+    }
+
 private:
     mutable Eigen::VectorX<FP> m_flow_values; ///< Cache to avoid allocation in get_derivatives (using get_flows).
 

cpp/memilio/compartments/parameter_studies.h

@@ -35,32 +35,51 @@
 #include <limits>
 #include <numeric>
 
+#include <type_traits>
+
+// Check if a type has a member called 'stay_duration'
+template <typename T, typename = void>
+struct has_stay_duration : std::false_type {
+};
+
+template <typename T>
+struct has_stay_duration<T, std::void_t<decltype(std::declval<T>().stay_duration)>> : std::true_type {
+};
+
+// Check if a type has a member called 'travel_time'
+template <typename T, typename = void>
+struct has_travel_time : std::false_type {
+};
+
+template <typename T>
+struct has_travel_time<T, std::void_t<decltype(std::declval<T>().travel_time)>> : std::true_type {
+};
+
+// Check if a type has a member called 'path'
+template <typename T, typename = void>
+struct has_path : std::false_type {
+};
+
+template <typename T>
+struct has_path<T, std::void_t<decltype(std::declval<T>().path)>> : std::true_type {
+};
+
 namespace mio
 {
 
 /**
  * Class that performs multiple simulation runs with randomly sampled parameters.
  * Can simulate mobility graphs with one simulation in each node or single simulations.
  * @tparam S type of simulation that runs in one node of the graph.
+ * @tparam ParametersGraph stores the parameters of the simulation. This is the input of ParameterStudies.
+ * @tparam SimulationGraph stores simulations and their results of each run. This is the output of ParameterStudies for each run.
  */
-template <class S>
+template <class S, class ParametersGraph = Graph<typename S::Model, MobilityParameters<double>>,
+          class SimulationGraph = Graph<SimulationNode<S>, MobilityEdge<double>>>
 class ParameterStudy
 {
 public:
-    /**
-    * The type of simulation of a single node of the graph.
-    */
     using Simulation = S;
-    /**
-    * The Graph type that stores the parametes of the simulation.
-    * This is the input of ParameterStudies.
-    */
-    using ParametersGraph = mio::Graph<typename Simulation::Model, mio::MobilityParameters<double>>;
-    /**
-    * The Graph type that stores simulations and their results of each run.
-    * This is the output of ParameterStudies for each run.
-    */
-    using SimulationGraph = mio::Graph<mio::SimulationNode<Simulation>, mio::MobilityEdge<double>>;
 
     /**
      * create study for graph of compartment models.
@@ -336,18 +355,84 @@ class ParameterStudy
     }
 
 private:
-    //sample parameters and create simulation
+    /**
+     * @brief Adds a node to the graph using the overload based on the node's property.
+     *
+     * This function checks whether the node property type has a 'stay_duration' member.
+     * If so, it uses the specialized version that extracts additional simulation parameters.
+     * Otherwise, we just use the simulation time (m_t0) and integration step size (m_dt_integration).
+     *
+     * @tparam GraphType The type of the graph.
+     * @tparam NodeType The type of the node.
+     * @param graph The graph to which the node will be added.
+     * @param node The node to add.
+     */
+    template <typename GraphType, typename NodeType>
+    void add_node(GraphType& graph, const NodeType& node)
+    {
+        // Deduce the property type of the node.
+        using PropertyType = std::decay_t<decltype(node.property)>;
+
+        // If the property type has a member called 'stay_duration', use the specialized version.
+        if constexpr (has_stay_duration<PropertyType>::value) {
+            graph.add_node(node.id, node.property.base_sim, node.property.mobility_sim, node.property.stay_duration);
+        }
+        else {
+            graph.add_node(node.id, node.property, m_t0, m_dt_integration);
+        }
+    }
+
+    /**
+     * @brief Adds an edge to the graph using the appropriate overload based on the edge's property.
+     *
+     * This function checks whether the edge property type has both 'travel_time' and 'path' members.
+     * If so, it uses the specialized version that uses additional parameters.
+     * Otherwise, we just use the edge property.
+     *
+     * @tparam GraphType The type of the graph.
+     * @tparam EdgeType The type of the edge.
+     * @param graph The graph to which the edge will be added.
+     * @param edge The edge to add.
+     */
+    template <typename GraphType, typename EdgeType>
+    void add_edge(GraphType& graph, const EdgeType& edge)
+    {
+        // Deduce the property type of the edge.
+        using PropertyType = std::decay_t<decltype(edge.property)>;
+
+        // Check if the property type has both 'travel_time' and 'path' members.
+        if constexpr (has_travel_time<PropertyType>::value && has_path<PropertyType>::value) {
+            graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property.get_parameters(),
+                           edge.property.travel_time, edge.property.path);
+        }
+        else {
+            graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property);
+        }
+    }
+
+    /**
+     * @brief Creates a mobility simulation based on a sampled graph.
+     * @tparam SampleGraphFunction The type of the function that samples the graph.
+     * @param sample_graph A function that takes the internal graph and returns a sampled graph.
+     * @return A mobility simulation created from the sampled simulation graph.
+     */
     template <class SampleGraphFunction>
-    mio::GraphSimulation<SimulationGraph> create_sampled_simulation(SampleGraphFunction sample_graph)
+    auto create_sampled_simulation(SampleGraphFunction sample_graph)
     {
+        // Create an empty simulation graph.
         SimulationGraph sim_graph;
 
+        // Sample the original graph using the provided sampling function.
         auto sampled_graph = sample_graph(m_graph);
+
+        // Iterate over each node in the sampled graph and add it to the simulation graph.
         for (auto&& node : sampled_graph.nodes()) {
-            sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration);
+            add_node(sim_graph, node);
         }
+
+        // Iterate over each edge in the sampled graph and add it to the simulation graph.
         for (auto&& edge : sampled_graph.edges()) {
-            sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property);
+            add_edge(sim_graph, edge);
         }
 
         return make_mobility_sim(m_t0, m_dt_graph_sim, std::move(sim_graph));

cpp/memilio/compartments/simulation.h

@@ -60,6 +60,19 @@ class Simulation
     {
     }
 
+    /**
+     * @brief Copy constructor
+     * @param[in] other Another simulation to copy
+     */
+    Simulation(const Simulation& other)
+        : m_integratorCore(other.m_integratorCore)
+        , m_model(std::make_unique<Model>(*other.m_model))
+        , m_integrator(m_integratorCore)
+        , m_result(other.m_result)
+        , m_dt(other.m_dt)
+    {
+    }
+
     /**
      * @brief Set the integrator core used in the simulation.
      * @param[in] integrator A shared pointer to an object derived from IntegratorCore.