1159 add diffusive abm and smm (#1162)

Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com>

add diffusive ABM (model, simulation, parameters and quadwell potential)
add SMM (model, simulation, parameters)
add tests for both models
add examples for both models

Author: jubicker

cpp/CMakeLists.txt

@@ -24,8 +24,8 @@ option(MEMILIO_ENABLE_PROFILING "Enable runtime performance profiling of memilio
 
 mark_as_advanced(MEMILIO_USE_BUNDLED_SPDLOG MEMILIO_SANITIZE_ADDRESS MEMILIO_SANITIZE_UNDEFINED)
 
-# try to treat AppleClang as Clang, but warn about missing support 
-if (CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
+# try to treat AppleClang as Clang, but warn about missing support
+if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
     message(WARNING "The compiler ID \"AppleClang\" is not supported, trying to compile with \"Clang\" options.")
     set(CMAKE_CXX_COMPILER_ID "Clang")
 endif()
@@ -82,6 +82,7 @@ if(CMAKE_BUILD_TYPE STREQUAL "Debug" OR CMAKE_BUILD_TYPE STREQUAL "DEBUG")
         message(STATUS "Coverage enabled")
         include(CodeCoverage)
         append_coverage_compiler_flags()
+
         # In addition to standard flags, disable elision and inlining to prevent e.g. closing brackets being marked as
         # uncovered.
         set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fno-elide-constructors -fno-default-inline")
@@ -149,6 +150,7 @@ add_subdirectory(memilio)
 
 if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/abm)
+    add_subdirectory(models/d_abm)
     add_subdirectory(models/ode_secir)
     add_subdirectory(models/ode_secirts)
     add_subdirectory(models/ode_secirvvs)
@@ -163,6 +165,7 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/sde_sirs)
     add_subdirectory(models/sde_seirvv)
     add_subdirectory(models/graph_abm)
+    add_subdirectory(models/smm)
 endif()
 
 if(MEMILIO_BUILD_EXAMPLES)

cpp/examples/CMakeLists.txt

@@ -132,6 +132,14 @@ add_executable(graph_abm_example graph_abm.cpp)
 target_link_libraries(graph_abm_example PRIVATE memilio graph_abm abm)
 target_compile_options(graph_abm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(dabm_example d_abm.cpp)
+target_link_libraries(dabm_example PRIVATE memilio d_abm)
+target_compile_options(dabm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(smm_example smm.cpp)
+target_link_libraries(smm_example PRIVATE memilio smm)
+target_compile_options(smm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 if(MEMILIO_HAS_JSONCPP)
     add_executable(ode_secir_read_graph_example ode_secir_read_graph.cpp)
     target_link_libraries(ode_secir_read_graph_example PRIVATE memilio ode_secir)
@@ -140,13 +148,13 @@ if(MEMILIO_HAS_JSONCPP)
 endif()
 
 if(MEMILIO_HAS_HDF5 AND MEMILIO_HAS_JSONCPP)
-  add_executable(ode_secir_parameter_study_example ode_secir_parameter_study.cpp)
-  target_link_libraries(ode_secir_parameter_study_example PRIVATE memilio ode_secir)
-  target_compile_options(ode_secir_parameter_study_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+    add_executable(ode_secir_parameter_study_example ode_secir_parameter_study.cpp)
+    target_link_libraries(ode_secir_parameter_study_example PRIVATE memilio ode_secir)
+    target_compile_options(ode_secir_parameter_study_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
-  add_executable(ode_secir_parameter_study_graph ode_secir_parameter_study_graph.cpp)
-  target_link_libraries(ode_secir_parameter_study_graph PRIVATE memilio ode_secir)
-  target_compile_options(ode_secir_parameter_study_graph PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+    add_executable(ode_secir_parameter_study_graph ode_secir_parameter_study_graph.cpp)
+    target_link_libraries(ode_secir_parameter_study_graph PRIVATE memilio ode_secir)
+    target_compile_options(ode_secir_parameter_study_graph PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 endif()
 
 if(MEMILIO_HAS_JSONCPP)

cpp/examples/d_abm.cpp

@@ -0,0 +1,86 @@
+/* 
+* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+*
+* Authors: Julia Bicker, René Schmieding
+*
+* 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 "d_abm/quad_well.h"
+#include "d_abm/simulation.h"
+#include "d_abm/parameters.h"
+#include "memilio/utils/random_number_generator.h"
+#include "memilio/data/analyze_result.h"
+#include "memilio/epidemiology/adoption_rate.h"
+#include <vector>
+
+enum class InfectionState
+{
+    S,
+    E,
+    C,
+    I,
+    R,
+    D,
+    Count
+
+};
+
+int main()
+{
+    //Example how to run a simulation of the diffusive ABM using the quadwell potential
+    using Model = mio::dabm::Model<QuadWellModel<InfectionState>>;
+    std::vector<Model::Agent> agents(1000);
+    //Random variables for initialization of agents' position and infection state
+    auto& pos_rng = mio::UniformDistribution<double>::get_instance();
+    auto& sta_rng = mio::DiscreteDistribution<size_t>::get_instance();
+    //Infection state distribution
+    std::vector<double> pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.};
+    for (auto& a : agents) {
+        //Agents are equally distributed in [-2,2]x[-2,2] at the beginning
+        a.position =
+            Eigen::Vector2d{pos_rng(mio::thread_local_rng(), -2., 2.), pos_rng(mio::thread_local_rng(), -2., 2.)};
+        a.status = static_cast<InfectionState>(sta_rng(mio::thread_local_rng(), pop_dist));
+    }
+
+    //Set adoption rates
+    std::vector<mio::AdoptionRate<InfectionState>> adoption_rates;
+    for (size_t region = 0; region < 4; ++region) {
+        adoption_rates.push_back({InfectionState::S,
+                                  InfectionState::E,
+                                  mio::regions::Region(region),
+                                  0.1,
+                                  {{InfectionState::C, 1}, {InfectionState::I, 0.5}}});
+        adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::regions::Region(region), 1.0 / 5., {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::regions::Region(region), 0.2 / 3., {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::regions::Region(region), 0.8 / 3., {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::regions::Region(region), 0.99 / 5., {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::regions::Region(region), 0.01 / 5., {}});
+    }
+
+    //Set interaction radius and noise term of the diffusion process
+    double interaction_radius = 0.5;
+    double noise              = 0.4;
+
+    double dt   = 0.1;
+    double tmax = 30.;
+
+    Model model(agents, adoption_rates, interaction_radius, noise, {InfectionState::D});
+    auto sim = mio::dabm::Simulation(model, 0.0, dt);
+    sim.advance(tmax);
+
+    auto interpolated_results = mio::interpolate_simulation_result(sim.get_result());
+    interpolated_results.print_table({"S", "E", "C", "I", "R", "D "});
+}

cpp/examples/smm.cpp

@@ -0,0 +1,99 @@
+/* 
+* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+*
+* Authors: Julia Bicker, René Schmieding
+*
+* 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 "smm/simulation.h"
+#include "smm/parameters.h"
+#include "memilio/data/analyze_result.h"
+#include "memilio/epidemiology/adoption_rate.h"
+
+enum class InfectionState
+{
+    S,
+    E,
+    C,
+    I,
+    R,
+    D,
+    Count
+
+};
+
+int main()
+{
+
+    //Example how to run the stochastic metapopulation models with four regions
+    const size_t num_regions = 4;
+    using Model              = mio::smm::Model<num_regions, InfectionState>;
+
+    double numE = 12, numC = 4, numI = 12, numR = 0, numD = 0;
+
+    Model model;
+    //Population are distributed uniformly to the four regions
+    for (size_t r = 0; r < num_regions; ++r) {
+        model.populations[{mio::regions::Region(r), InfectionState::S}] =
+            (1000 - numE - numC - numI - numR - numD) / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::E}] = numE / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::C}] = numC / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::I}] = numI / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::R}] = numR / num_regions;
+        model.populations[{mio::regions::Region(r), InfectionState::D}] = numD / num_regions;
+    }
+
+    //Set infection state adoption and spatial transition rates
+    std::vector<mio::AdoptionRate<InfectionState>> adoption_rates;
+    std::vector<mio::smm::TransitionRate<InfectionState>> transition_rates;
+    for (size_t r = 0; r < num_regions; ++r) {
+        adoption_rates.push_back({InfectionState::S,
+                                  InfectionState::E,
+                                  mio::regions::Region(r),
+                                  0.1,
+                                  {{InfectionState::C, 1}, {InfectionState::I, 0.5}}});
+        adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}});
+    }
+
+    //Agents in infection state D do not transition
+    for (size_t s = 0; s < static_cast<size_t>(InfectionState::D); ++s) {
+        for (size_t i = 0; i < num_regions; ++i) {
+            for (size_t j = 0; j < num_regions; ++j)
+                if (i != j) {
+                    transition_rates.push_back(
+                        {InfectionState(s), mio::regions::Region(i), mio::regions::Region(j), 0.01});
+                    transition_rates.push_back(
+                        {InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01});
+                }
+        }
+    }
+
+    model.parameters.get<mio::smm::AdoptionRates<InfectionState>>()   = adoption_rates;
+    model.parameters.get<mio::smm::TransitionRates<InfectionState>>() = transition_rates;
+
+    double dt   = 0.1;
+    double tmax = 30.;
+
+    auto sim = mio::smm::Simulation(model, 0.0, dt);
+    sim.advance(tmax);
+
+    auto interpolated_results = mio::interpolate_simulation_result(sim.get_result());
+    interpolated_results.print_table({"S", "E", "C", "I", "R", "D "});
+}

cpp/memilio/CMakeLists.txt

@@ -18,6 +18,7 @@ add_library(memilio
     epidemiology/dynamic_npis.cpp
     epidemiology/lct_infection_state.h
     epidemiology/lct_populations.h
+    epidemiology/adoption_rate.h
     geography/regions.h
     geography/regions.cpp
     epidemiology/simulation_day.h

cpp/memilio/epidemiology/adoption_rate.h

@@ -0,0 +1,59 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: René Schmieding, Julia Bicker
+*
+* 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.
+*/
+#ifndef MIO_EPI_ADOPTIONRATE_H
+#define MIO_EPI_ADOPTIONRATE_H
+
+#include "memilio/utils/index.h"
+#include "memilio/config.h"
+#include "memilio/geography/regions.h"
+
+namespace mio
+{
+
+/**
+ * @brief Struct defining an influence for a second-order adoption.
+ * The population having "status" is multiplied with "factor."
+ * @tparam Status An infection state enum.
+ */
+template <class Status>
+struct Influence {
+    Status status;
+    ScalarType factor;
+};
+
+/**
+ * @brief Struct defining a possible status adoption in a Model based on Poisson Processes.
+ * The AdoptionRate is considered to be of second-order if there are any "influences".
+ * In the d_abm and smm simulations, "from" is implicitly an influence, scaled by "factor". This is multiplied by
+ * the sum over all "influences", which scale their "status" with the respective "factor".
+ * @tparam Status An infection state enum.
+ */
+template <class Status>
+struct AdoptionRate {
+    Status from; // i
+    Status to; // j
+    mio::regions::Region region; // k
+    ScalarType factor; // gammahat_{ij}^k
+    std::vector<Influence<Status>> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )
+};
+
+} // namespace mio
+
+#endif

cpp/memilio/geography/regions.h

@@ -23,6 +23,7 @@
 #include "memilio/utils/date.h"
 #include "memilio/utils/stl_util.h"
 #include "memilio/utils/type_safe.h"
+#include "memilio/utils/index.h"
 
 #include "boost/filesystem.hpp"
 
@@ -37,6 +38,14 @@ namespace mio
 namespace regions
 {
 
+/// @biref Index for enumerating subregions (cities, counties, etc.) of the modelled area.
+struct Region : public mio::Index<Region> {
+    Region(const size_t num_regions)
+        : mio::Index<Region>(num_regions)
+    {
+    }
+};
+
 /**
          * Id of a state.
          * For Germany the Ids are:

cpp/memilio/utils/random_number_generator.h

@@ -662,6 +662,13 @@ using DiscreteDistribution = DistributionAdapter<DiscreteDistributionInPlace<Int
 template <class Real>
 using ExponentialDistribution = DistributionAdapter<std::exponential_distribution<Real>>;
 
+/**
+ * adapted std::normal_distribution.
+ * @see DistributionAdapter
+ */
+template <class Real>
+using NormalDistribution = DistributionAdapter<std::normal_distribution<Real>>;
+
 /**
  * adapted std::uniform_int_distribution.
  * @see DistributionAdapter

cpp/models/README.md

@@ -1,3 +1,15 @@
 # MEmilio Models #
 
-Contains different concrete models that are built using the MEmilio C++ library. See the corresponding directory for details about each model.
+Contains different concrete models that are built using the MEmilio C++ library. Some models contain a spatial resolution and some do not contain spatial resolution, hence cannot capture spatially heterogenous dynamics.
+The models with spatial resolution are:
+  - abm
+  - d_abm
+
+The models without spatial resolution are:
+  - ode_*
+  - ide_*
+  - lct_*
+  - glct_*
+  - sde_*
+  
+See the corresponding directory for details about each model.

cpp/models/d_abm/CMakeLists.txt

@@ -0,0 +1,13 @@
+add_library(d_abm
+    model.h
+    model.cpp
+    simulation.h
+    simulation.cpp
+    parameters.h
+)
+target_link_libraries(d_abm PUBLIC memilio)
+target_include_directories(d_abm PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+target_compile_options(d_abm PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})