1159 add diffusive abm and smm

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)
@@ -162,6 +164,7 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/sde_sir)
     add_subdirectory(models/sde_sirs)
     add_subdirectory(models/sde_seirvv)
+    add_subdirectory(models/smm)
 endif()
 
 if(MEMILIO_BUILD_EXAMPLES)

cpp/examples/CMakeLists.txt

@@ -106,7 +106,6 @@ target_link_libraries(twitter_mobility_example PRIVATE memilio ode_secir)
 
 target_compile_options(twitter_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
-
 add_executable(ide_seir_example ide_seir.cpp)
 target_link_libraries(ide_seir_example PRIVATE memilio ide_seir)
 target_compile_options(ide_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
@@ -131,6 +130,14 @@ add_executable(history_example history.cpp)
 target_link_libraries(history_example PRIVATE memilio)
 target_compile_options(history_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(dabm_example d_abm.cpp)
+target_link_libraries(dabm_example PRIVATE memilio dabm)
+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)
@@ -139,13 +146,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)
@@ -167,7 +174,7 @@ if(MEMILIO_HAS_HDF5)
 endif()
 
 if(MEMILIO_HAS_JSONCPP)
-  add_executable(ide_initialization_example ide_initialization.cpp)
-  target_link_libraries(ide_initialization_example PRIVATE memilio ide_secir)
-  target_compile_options(ide_initialization_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+    add_executable(ide_initialization_example ide_initialization.cpp)
+    target_link_libraries(ide_initialization_example PRIVATE memilio ide_secir)
+    target_compile_options(ide_initialization_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 endif()

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 <[email protected]>
+*
+* 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 <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::dabm::AdoptionRate<InfectionState>> adoption_rates;
+    for (size_t region = 0; region < 4; ++region) {
+        adoption_rates.push_back({InfectionState::S,
+                                  InfectionState::E,
+                                  mio::dabm::Region(region),
+                                  0.1,
+                                  {InfectionState::C, InfectionState::I},
+                                  {1, 0.5}});
+        adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::dabm::Region(region), 1.0 / 5., {}, {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::dabm::Region(region), 0.2 / 3., {}, {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::dabm::Region(region), 0.8 / 3., {}, {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::dabm::Region(region), 0.99 / 5., {}, {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::dabm::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::Simulation<double, Model>(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,97 @@
+/* 
+* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+*
+* Authors: Julia Bicker, René Schmieding
+*
+* Contact: Martin J. Kuehn <[email protected]>
+*
+* 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"
+
+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::smm::Region(r), InfectionState::S}] =
+            (1000 - numE - numC - numI - numR - numD) / num_regions;
+        model.populations[{mio::smm::Region(r), InfectionState::E}] = numE / num_regions;
+        model.populations[{mio::smm::Region(r), InfectionState::C}] = numC / num_regions;
+        model.populations[{mio::smm::Region(r), InfectionState::I}] = numI / num_regions;
+        model.populations[{mio::smm::Region(r), InfectionState::R}] = numR / num_regions;
+        model.populations[{mio::smm::Region(r), InfectionState::D}] = numD / num_regions;
+    }
+
+    //Set infection state adoption and spatial transition rates
+    std::vector<mio::smm::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::smm::Region(r),
+                                  0.1,
+                                  {InfectionState::C, InfectionState::I},
+                                  {1, 0.5}});
+        adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::smm::Region(r), 1.0 / 5., {}, {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::smm::Region(r), 0.2 / 3., {}, {}});
+        adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::smm::Region(r), 0.8 / 3., {}, {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::smm::Region(r), 0.99 / 5., {}, {}});
+        adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::smm::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::smm::Region(i), mio::smm::Region(j), 0.01});
+                    transition_rates.push_back({InfectionState(s), mio::smm::Region(j), mio::smm::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::Simulation<double, Model>(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/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(dabm
+    model.h
+    model.cpp
+    simulation.h
+    simulation.cpp
+    parameters.h
+)
+target_link_libraries(dabm PUBLIC memilio)
+target_include_directories(dabm PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+target_compile_options(dabm PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/models/d_abm/README.md

@@ -0,0 +1,14 @@
+# Diffusive Agent-Based Model
+
+This agent-based model uses a Markov process to simulate disease dynamics. The features of an agent are position and infection state. The evolution of the system state is determined by the following master equation
+
+```math
+\partial_t p(X,Z;t) = G p(X,Z;t) + L p(X,Z;t)
+```
+The operator $G$ defines the infection state adoptions and only acts on $Z$, while $L$ defines location changes, only acting on $X$. Infection state adoptions are modeled with independent Poisson processes given by adoption rate functions. Movement is modeled with independent diffusion processes. A temporal Gillespie algorithm is used for simulation.
+
+The Model class needs an Implementation class as template argument which prvides the domain agents move and interact in. We here implemented a quadwell potential given in the class QuadWellModel, but any other suitable potential can be used as implementation. 
+
+## Simulation
+
+The simulation runs in discrete time steps. In every step we advance the model until the next infection state adoption event, then adopt the corresponding agent's infection state and draw a new waiting time until the next adoption event. If the waiting time until the next adoption event is bigger than the remaining time in the time step, we advance the model until the end of the time step.

cpp/models/d_abm/model.cpp

@@ -0,0 +1,29 @@
+/* 
+* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+*
+* Authors: René Schmieding
+*
+* Contact: Martin J. Kuehn <[email protected]>
+*
+* 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/model.h"
+
+namespace mio
+{
+namespace dabm
+{
+
+} // namespace dabm
+} // namespace mio