385 mobility scheme with traveltimes and transmission during changing nodes

cpp/examples/CMakeLists.txt

@@ -21,6 +21,8 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI
 add_executable(ode_sir_example ode_sir.cpp)
 target_link_libraries(ode_sir_example PRIVATE memilio ode_sir)
 
+add_executable(sim_mobility_detailed ode_secirvvs_mobility_detailed.cpp)
+target_link_libraries(sim_mobility_detailed PRIVATE memilio ode_secirvvs)
 
 add_executable(seir_flows_example ode_seir_flows.cpp)
 target_link_libraries(seir_flows_example PRIVATE memilio ode_seir)

cpp/memilio/CMakeLists.txt

@@ -59,6 +59,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/parameter_studies.h

@@ -38,6 +38,22 @@
 namespace mio
 {
 
+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 {
+};
+
+template <typename T, typename = void>
+struct has_traveltime : std::false_type {
+};
+
+template <typename T>
+struct has_traveltime<T, std::void_t<decltype(std::declval<T>().traveltime)>> : std::true_type {
+};
+
 /**
  * Class that performs multiple simulation runs with randomly sampled parameters.
  * Can simulate migration graphs with one simulation in each node or single simulations.
@@ -140,16 +156,17 @@ class ParameterStudy
 #else
         num_procs = 1;
         rank      = 0;
-#endif  
-        
+#endif
+
         //The ParameterDistributions used for sampling parameters use thread_local_rng()
         //So we set our own RNG to be used.
         //Assume that sampling uses the thread_local_rng() and isn't multithreaded
         m_rng.synchronize();
         thread_local_rng() = m_rng;
 
         auto run_distribution = distribute_runs(m_num_runs, num_procs);
-        auto start_run_idx = std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0));
+        auto start_run_idx =
+            std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0));
         auto end_run_idx = start_run_idx + run_distribution[size_t(rank)];
 
         std::vector<std::invoke_result_t<HandleSimulationResultFunction, SimulationGraph, size_t>> ensemble_result;
@@ -167,7 +184,8 @@ class ParameterStudy
 
             //sample
             auto sim = create_sampled_simulation(sample_graph);
-            log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", (thread_local_rng().get_counter() - run_rng_counter).get());
+            log(LogLevel::info, "ParameterStudies: Generated {} random numbers.",
+                (thread_local_rng().get_counter() - run_rng_counter).get());
 
             //perform run
             sim.advance(m_tmax);
@@ -328,10 +346,31 @@ class ParameterStudy
     }
     /** @} */
 
-    RandomNumberGenerator& get_rng() {
+    RandomNumberGenerator& get_rng()
+    {
         return m_rng;
     }
 
+protected:
+    // adaptive time step of the integrator (will be corrected if too large/small)
+    double m_dt_integration = 0.1;
+    //
+    RandomNumberGenerator m_rng;
+
+    std::vector<size_t> distribute_runs(size_t num_runs, int num_procs)
+    {
+        //evenly distribute runs
+        //lower processes do one more run if runs are not evenly distributable
+        auto num_runs_local = num_runs / num_procs; //integer division!
+        auto remainder      = num_runs % num_procs;
+
+        std::vector<size_t> run_distribution(num_procs);
+        std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1);
+        std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local);
+
+        return run_distribution;
+    }
+
 private:
     //sample parameters and create simulation
     template <class SampleGraphFunction>
@@ -341,29 +380,25 @@ class ParameterStudy
 
         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);
+            if constexpr (has_stay_duration<decltype(node)>::value) {
+                sim_graph.add_node(node.id, node.stay_duration, node.property, m_t0, m_dt_integration);
+            }
+            else {
+                sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration);
+            }
         }
         for (auto&& edge : sampled_graph.edges()) {
-            sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property);
+            if constexpr (has_traveltime<decltype(edge)>::value) {
+                sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.property);
+            }
+            else {
+                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));
     }
 
-    std::vector<size_t> distribute_runs(size_t num_runs, int num_procs)
-    {
-        //evenly distribute runs
-        //lower processes do one more run if runs are not evenly distributable
-        auto num_runs_local = num_runs / num_procs; //integer division!
-        auto remainder = num_runs % num_procs;
-
-        std::vector<size_t> run_distribution(num_procs);
-        std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1);
-        std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local);
-
-        return run_distribution;
-    }
-
 private:
     // Stores Graph with the names and ranges of all parameters
     ParametersGraph m_graph;
@@ -376,10 +411,6 @@ class ParameterStudy
     double m_tmax;
     // time step of the graph
     double m_dt_graph_sim;
-    // adaptive time step of the integrator (will be corrected if too large/small)
-    double m_dt_integration = 0.1;
-    //
-    RandomNumberGenerator m_rng;
 };
 
 } // namespace mio

cpp/memilio/compartments/simulation.h

@@ -59,6 +59,20 @@ class Simulation
     {
     }
 
+    /**
+     * @brief Copy constructor.
+     * @param[in] other The 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)
+    {
+        m_integrator.set_integrator(m_integratorCore);
+    }
+
     /**
      * @brief Set the integrator core used in the simulation.
      * @param[in] integrator A shared pointer to an object derived from IntegratorCore.

cpp/memilio/data/analyze_result.h

@@ -22,6 +22,7 @@
 
 #include "memilio/utils/time_series.h"
 #include "memilio/mobility/metapopulation_mobility_instant.h"
+#include "memilio/mobility/metapopulation_mobility_detailed.h"
 
 #include <functional>
 #include <vector>
@@ -120,6 +121,18 @@ interpolate_simulation_result(const Graph<SimulationNode<Simulation>, MigrationE
     return interpolated;
 }
 
+template <class Simulation>
+std::vector<TimeSeries<double>>
+interpolate_simulation_result(const GraphDetailed<SimulationNode<Simulation>, MigrationEdge>& graph_result)
+{
+    std::vector<TimeSeries<double>> interpolated;
+    interpolated.reserve(graph_result.nodes().size());
+    std::transform(graph_result.nodes().begin(), graph_result.nodes().end(), std::back_inserter(interpolated),
+                   [](auto& n) {
+                       return interpolate_simulation_result(n.property.get_result());
+                   });
+    return interpolated;
+}
 /**
  * Compute the distance between two SECIR simulation results.
  * The distance is the 2-norm of the element-wise difference of the two results.

cpp/memilio/io/mobility_io.cpp

@@ -135,4 +135,93 @@ IOResult<Eigen::MatrixXd> read_mobility_plain(const std::string& filename)
     return success(migration);
 }
 
+IOResult<Eigen::MatrixXd> read_duration_stay(const std::string& filename)
+{
+    BOOST_OUTCOME_TRY(num_lines, count_lines(filename));
+
+    if (num_lines == 0) {
+        return success(Eigen::MatrixXd(0, 0));
+    }
+
+    std::fstream file;
+    file.open(filename, std::ios::in);
+    if (!file.is_open()) {
+        return failure(StatusCode::FileNotFound, filename);
+    }
+
+    Eigen::VectorXd duration(num_lines);
+
+    try {
+        std::string tp;
+        int linenumber = 0;
+        while (getline(file, tp)) {
+            auto line      = split(tp, ' ');
+            Eigen::Index i = static_cast<Eigen::Index>(linenumber);
+            duration(i)    = std::stod(line[0]);
+            linenumber++;
+        }
+    }
+    catch (std::runtime_error& ex) {
+        return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what());
+    }
+
+    return success(duration);
+}
+
+IOResult<std::vector<std::vector<std::vector<int>>>> read_path_mobility(const std::string& filename)
+{
+    BOOST_OUTCOME_TRY(num_lines, count_lines(filename));
+
+    if (num_lines == 0) {
+        std::vector<std::vector<std::vector<int>>> arr(0, std::vector<std::vector<int>>(0));
+        return success(arr);
+    }
+
+    std::fstream file;
+    file.open(filename, std::ios::in);
+    if (!file.is_open()) {
+        return failure(StatusCode::FileNotFound, filename);
+    }
+
+    const int num_nodes = static_cast<int>(std::sqrt(num_lines));
+    std::vector<std::vector<std::vector<int>>> arr(num_nodes, std::vector<std::vector<int>>(num_nodes));
+
+    try {
+        std::string tp;
+        while (getline(file, tp)) {
+            auto line   = split(tp, ' ');
+            auto indx_x = std::stoi(line[0]);
+            auto indx_y = std::stoi(line[1]);
+            if (indx_x != indx_y) {
+                auto path = std::accumulate(line.begin() + 2, line.end(), std::string(""));
+
+                // string -> vector of integers
+                std::vector<int> path_vec;
+
+                // Remove the square brackets and \r
+                path = path.substr(1, path.size() - 3);
+                std::stringstream ss(path);
+                std::string token;
+
+                // get numbers and save them in path_vec
+                while (std::getline(ss, token, ',')) {
+                    path_vec.push_back(std::stoi(token));
+                }
+
+                for (int number : path_vec) {
+                    arr[indx_x][indx_y].push_back(number);
+                }
+            }
+            else {
+                arr[indx_x][indx_y].push_back(static_cast<int>(indx_x));
+            }
+        }
+    }
+    catch (std::runtime_error& ex) {
+        return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what());
+    }
+
+    return success(arr);
+}
+
 } // namespace mio

cpp/memilio/io/mobility_io.h

@@ -17,8 +17,8 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#ifndef READ_TWITTER_H
-#define READ_TWITTER_H
+#ifndef MIO_MOBILITY_IO_H
+#define MIO_MOBILITY_IO_H
 
 #include "memilio/config.h"
 #include "memilio/math/eigen.h"
@@ -63,6 +63,22 @@ IOResult<Eigen::MatrixXd> read_mobility_formatted(const std::string& filename);
  */
 IOResult<Eigen::MatrixXd> read_mobility_plain(const std::string& filename);
 
+/**
+ * @brief Reads txt file containing the duration of stay in each county.
+          Writes it into a Eigen vector of size N, where N is the number of local entites.
+ * @param filename name of file to be read
+ * @return IOResult<Eigen::MatrixXd> containing the duration of stay in each local entity
+ */
+IOResult<Eigen::MatrixXd> read_duration_stay(const std::string& filename);
+
+/**
+ * @brief For each edge we have the path from the start node to the end node. This functions reads the file and returns the path for each edge.
+ * 
+ * @param filename Filename of the file containing the paths.
+ * @return IOResult<std::vector<std::vector<std::vector<int>>>> contains the paths for each edge. 
+ */
+IOResult<std::vector<std::vector<std::vector<int>>>> read_path_mobility(const std::string& filename);
+
 #ifdef MEMILIO_HAS_JSONCPP
 
 /**
@@ -133,7 +149,8 @@ IOResult<void> write_graph(const Graph<Model, MigrationParameters>& graph, const
  * @param read_edges boolean value that decides whether the edges of the graph should also be read in.
  */
 template <class Model>
-IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& directory, int ioflags = IOF_None, bool read_edges = true)
+IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& directory, int ioflags = IOF_None,
+                                                       bool read_edges = true)
 {
     std::string abs_path;
     if (!file_exists(directory, abs_path)) {
@@ -160,7 +177,7 @@ IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& direct
     }
 
     //read edges; nodes must already be available for that)
-    if(read_edges){
+    if (read_edges) {
         for (auto inode = size_t(0); inode < graph.nodes().size(); ++inode) {
             //list of edges
             auto edge_filename = path_join(abs_path, "GraphEdges_node" + std::to_string(inode) + ".json");
@@ -177,7 +194,7 @@ IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& direct
                 if (end_node_idx >= graph.nodes().size()) {
                     log_error("EndNodeIndex not in range of number of graph nodes.");
                     return failure(StatusCode::OutOfRange,
-                                edge_filename + ", EndNodeIndex not in range of number of graph nodes.");
+                                   edge_filename + ", EndNodeIndex not in range of number of graph nodes.");
                 }
                 BOOST_OUTCOME_TRY(parameters, deserialize_json(e["Parameters"], Tag<MigrationParameters>{}, ioflags));
                 graph.add_edge(start_node_idx, end_node_idx, parameters);
@@ -192,4 +209,4 @@ IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& direct
 
 } // namespace mio
 
-#endif // READ_TWITTER_H
+#endif // MIO_MOBILITY_IO_H