From 3b8cedb9c7cdd661f9bfd61f11300cd3065d2562 Mon Sep 17 00:00:00 2001 From: jubicker Date: Thu, 12 Dec 2024 12:49:54 +0100 Subject: [PATCH 1/8] add diffusive ABM and tests --- cpp/CMakeLists.txt | 6 +- cpp/models/d_abm/CMakeLists.txt | 13 +++ cpp/models/d_abm/model.cpp | 29 +++++ cpp/models/d_abm/model.h | 52 +++++++++ cpp/models/d_abm/parameters.h | 64 ++++++++++ cpp/models/d_abm/quad_well.h | 201 ++++++++++++++++++++++++++++++++ cpp/models/d_abm/simulation.cpp | 24 ++++ cpp/models/d_abm/simulation.h | 162 +++++++++++++++++++++++++ cpp/tests/CMakeLists.txt | 3 +- cpp/tests/test_d_abm_model.cpp | 123 +++++++++++++++++++ 10 files changed, 674 insertions(+), 3 deletions(-) create mode 100644 cpp/models/d_abm/CMakeLists.txt create mode 100644 cpp/models/d_abm/model.cpp create mode 100644 cpp/models/d_abm/model.h create mode 100644 cpp/models/d_abm/parameters.h create mode 100644 cpp/models/d_abm/quad_well.h create mode 100644 cpp/models/d_abm/simulation.cpp create mode 100644 cpp/models/d_abm/simulation.h create mode 100644 cpp/tests/test_d_abm_model.cpp diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 8becf01e39..9228d92cc8 100644 --- a/cpp/CMakeLists.txt +++ b/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) diff --git a/cpp/models/d_abm/CMakeLists.txt b/cpp/models/d_abm/CMakeLists.txt new file mode 100644 index 0000000000..f69f3aecb8 --- /dev/null +++ b/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 + $ + $ +) +target_compile_options(dabm PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/d_abm/model.cpp b/cpp/models/d_abm/model.cpp new file mode 100644 index 0000000000..bcf658fd21 --- /dev/null +++ b/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 +* +* 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 diff --git a/cpp/models/d_abm/model.h b/cpp/models/d_abm/model.h new file mode 100644 index 0000000000..ec828b2360 --- /dev/null +++ b/cpp/models/d_abm/model.h @@ -0,0 +1,52 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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_D_ABM_MODEL_H +#define MIO_D_ABM_MODEL_H + +namespace mio +{ + +namespace dabm +{ +template +class Model : public Implementation +{ +public: + using Implementation::Implementation; + + using Implementation::adopt; + using Implementation::adoption_rate; + using Implementation::get_rng; + using Implementation::move; + using Implementation::time_point; + + using Status = typename Implementation::Status; + using Agent = typename Implementation::Agent; + + inline constexpr void check_constraints() const + { + } +}; + +} // namespace dabm +} // namespace mio + +#endif diff --git a/cpp/models/d_abm/parameters.h b/cpp/models/d_abm/parameters.h new file mode 100644 index 0000000000..b89f7800e6 --- /dev/null +++ b/cpp/models/d_abm/parameters.h @@ -0,0 +1,64 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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_d_ABM_PARAMETERS_H +#define MIO_d_ABM_PARAMETERS_H + +#include "memilio/config.h" +#include "memilio/utils/index.h" + +namespace mio +{ +namespace dabm +{ +// yet another "curiously recuring template pattern" +struct Region : public mio::Index { + Region(const size_t num_regions) + : mio::Index(num_regions) + { + } +}; + +// the AdoptionRate is considered to be of second order, if there are any "influences" with corresponding "factors". +// "from" is always an influence, scaled by "factor". +template +struct AdoptionRate { + Status from; // i + Status to; // j + Region region; // k + ScalarType factor; // gammahat_{ij}^k + std::vector influences; + std::vector factors; +}; + +template +struct AdoptionRates { + using Type = std::vector>; + const static std::string name() + { + return "AdoptionRates"; + } +}; + +} // namespace dabm + +} // namespace mio + +#endif diff --git a/cpp/models/d_abm/quad_well.h b/cpp/models/d_abm/quad_well.h new file mode 100644 index 0000000000..1f78ef3b95 --- /dev/null +++ b/cpp/models/d_abm/quad_well.h @@ -0,0 +1,201 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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 QUAD_WELL_H +#define QUAD_WELL_H + +#include "memilio/math/eigen.h" +#include "d_abm/parameters.h" +#include "memilio/utils/random_number_generator.h" +#include +#include + +using Position = Eigen::Vector2d; + +inline size_t well_index(const Position& p) +{ + // 0|1 + // -+- + // 2|3 + return (p.x() >= 0) + 2 * (p.y() < 0); +} + +template +class QuadWellModel +{ + +public: + using Status = InfectionState; + + struct Agent { + Position position; + Status status; + }; + + QuadWellModel(const std::vector& agents, const typename mio::dabm::AdoptionRates::Type& rates, + double contact_radius = 0.4, double sigma = 0.4, std::vector non_moving_states = {}) + : populations(agents) + , m_contact_radius(contact_radius) + , m_sigma(sigma) + , m_non_moving_states(non_moving_states) + , m_number_transitions(static_cast(Status::Count), Eigen::MatrixXd::Zero(4, 4)) + { + for (auto& agent : populations) { + mio::unused(agent); + assert(is_in_domain(agent.position)); + } + for (auto& r : rates) { + m_adoption_rates.emplace(std::forward_as_tuple(r.region, r.from, r.to), r); + //m_adoption_rates[{r.region, r.from, r.to}] = r; + } + } + + inline static constexpr void adopt(Agent& agent, const Status& new_status) + { + agent.status = new_status; + } + + double adoption_rate(const Agent& agent, const Status& new_status) const + { + double rate = 0; + // get the correct adoption rate + const size_t well = well_index(agent.position); + auto map_itr = m_adoption_rates.find({well, agent.status, new_status}); + if (map_itr != m_adoption_rates.end()) { + const auto& adoption_rate = map_itr->second; + // calculate the current rate, depending on order + if (adoption_rate.influences.size() == 0) { // first order adoption + // contact independant rate + rate = adoption_rate.factor; + } + else { // second order adoption + // accumulate rate per contact with a status in influences + size_t num_contacts = 0; + ScalarType influences = 0; + for (auto& contact : populations) { + // check if contact is indeed a contact + if (is_contact(agent, contact)) { + num_contacts++; + for (size_t i = 0; i < adoption_rate.influences.size(); i++) { + if (contact.status == adoption_rate.influences[i]) { + influences += adoption_rate.factors[i]; + } + } + } + } + // rate = factor * "concentration of contacts with status new_status" + if (num_contacts > 0) { + rate = adoption_rate.factor * (influences / num_contacts); + } + } + } + // else: no adoption from agent.status to new_status exist + return rate; + } + + void move(const double /*t*/, const double dt, Agent& agent) + { + const auto old_well = well_index(agent.position); + if (std::find(m_non_moving_states.begin(), m_non_moving_states.end(), agent.status) == + m_non_moving_states.end() && + std::find(m_non_moving_regions.begin(), m_non_moving_regions.end(), old_well) == + m_non_moving_regions.end()) { + Position p = {mio::DistributionAdapter>::get_instance()(m_rng, 0.0, 1.0), + mio::DistributionAdapter>::get_instance()(m_rng, 0.0, 1.0)}; + + agent.position = agent.position - dt * grad_U(agent.position) + (m_sigma * std::sqrt(dt)) * p; + const auto new_well = well_index(agent.position); + if (old_well != new_well) { + m_number_transitions[static_cast(agent.status)](old_well, new_well)++; + } + } + //else{agent has non-moving status or region} + } + + Eigen::VectorXd time_point() const + { + Eigen::VectorXd val = Eigen::VectorXd::Zero(4 * static_cast(Status::Count)); + for (auto& agent : populations) { + // split population into the wells given by grad_U + auto position = + static_cast(agent.status) + well_index(agent.position) * static_cast(Status::Count); + val[position] += 1; + } + return val; + } + + const std::vector& number_transitions() const + { + return m_number_transitions; + } + + std::vector& number_transitions() + { + return m_number_transitions; + } + + std::map, mio::dabm::AdoptionRate>& get_adoption_rates() + { + return m_adoption_rates; + } + + void set_non_moving_regions(std::vector non_moving_regions) + { + m_non_moving_regions = non_moving_regions; + } + + mio::RandomNumberGenerator& get_rng() + { + return m_rng; + } + + std::vector populations; + +private: + static Position grad_U(const Position x) + { + // U is a quad well potential + // U(x0,x1) = (x0^2 - 1)^2 + (x1^2 - 1)^2 + return {4 * x[0] * (x[0] * x[0] - 1), 4 * x[1] * (x[1] * x[1] - 1)}; + } + + bool is_contact(const Agent& agent, const Agent& contact) const + { + // test if contact is in the contact radius and test if agent and contact are different objects + return (agent.position - contact.position).norm() < m_contact_radius && (&agent != &contact) && + well_index(agent.position) == well_index(contact.position); + } + + bool is_in_domain(const Position& p) const + { + // restrict domain to [-2, 2]^2 where "escaping" is impossible, i.e. it holds x <= grad_U(x) for dt <= 0.1 + return -2 <= p[0] && p[0] <= 2 && -2 <= p[1] && p[1] <= 2; + } + + std::map, mio::dabm::AdoptionRate> m_adoption_rates; + double m_contact_radius; + double m_sigma; + std::vector m_non_moving_states; + std::vector m_non_moving_regions{}; + std::vector m_number_transitions; + mio::RandomNumberGenerator m_rng; +}; + +#endif diff --git a/cpp/models/d_abm/simulation.cpp b/cpp/models/d_abm/simulation.cpp new file mode 100644 index 0000000000..ddf6983917 --- /dev/null +++ b/cpp/models/d_abm/simulation.cpp @@ -0,0 +1,24 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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/simulation.h" + +namespace mio +{ +} // namespace mio diff --git a/cpp/models/d_abm/simulation.h b/cpp/models/d_abm/simulation.h new file mode 100644 index 0000000000..0e36f64907 --- /dev/null +++ b/cpp/models/d_abm/simulation.h @@ -0,0 +1,162 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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_D_ABM_SIMULATION_H +#define MIO_D_ABM_SIMULATION_H + +#include "d_abm/model.h" +#include "memilio/config.h" +#include "memilio/compartments/simulation.h" +#include "memilio/utils/random_number_generator.h" + +namespace mio +{ + +template +class Simulation> +{ + using Status = typename Implementation::Status; + +public: + using Model = dabm::Model; + + Simulation(const Model& model, ScalarType t0 = 0, ScalarType dt = 0.1) + : m_t0(t0) + , m_dt(dt) + , m_model(std::make_unique(model)) + , m_result(t0, m_model->time_point()) + { + assert(dt > 0); + m_current_rates.reserve(m_model->populations.size()); + m_current_events.reserve(m_model->populations.size()); + } + + void advance(const ScalarType t_max) + { + // draw time until an adoption takes place + ScalarType waiting_time = mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); + while (m_t0 < t_max) { + ScalarType dt = std::min({m_dt, t_max - m_t0}); + ScalarType remaining_time = dt; + while (remaining_time > 0) { + compute_current_rates_and_events(); // lambda_k (aka f-hat(N)) + ScalarType cumulative_adoption_rate = + std::accumulate(m_current_rates.begin(), m_current_rates.end(), 0.0); // Lambda + // status update + if (waiting_time < cumulative_adoption_rate * remaining_time) { + // draw which adoption takes place + const size_t event_id = + mio::DiscreteDistribution::get_instance()(m_model->get_rng(), m_current_rates); + Event& event = m_current_events[event_id]; + // perform adoption + m_model->adopt(event.agent, event.new_status); + // draw new waiting time + remaining_time -= waiting_time / cumulative_adoption_rate; + waiting_time = mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); + } + else { + // no event, decrease waiting time + waiting_time -= cumulative_adoption_rate * remaining_time; + break; + } + } + // position update + for (auto& agent : m_model->populations) { + m_model->move(m_t0, dt, agent); + } + m_t0 += dt; + // store result + m_result.add_time_point(m_t0, m_model->time_point()); + } + } + + void set_integrator(std::shared_ptr> /*integrator*/) + { + } + + /** + * @brief get_result returns the final simulation result + * @return a TimeSeries to represent the final simulation result + */ + TimeSeries& get_result() + { + return m_result; + } + + /** + * @brief get_result returns the final simulation result + * @return a TimeSeries to represent the final simulation result + */ + const TimeSeries& get_result() const + { + return m_result; + } + + /** + * @brief returns the simulation model used in simulation + */ + const Model& get_model() const + { + return *m_model; + } + + /** + * @brief returns the simulation model used in simulation + */ + Model& get_model() + { + return *m_model; + } + +private: + struct Event { + typename Model::Agent& agent; + Status new_status; + }; + /// @brief calculate values for m_current_rates and m_current_events + inline void compute_current_rates_and_events() + { + m_current_rates.clear(); + m_current_events.clear(); + // compute rate for each (agent, status) combination + for (auto& agent : m_model->populations) { + for (int s = 0; s < static_cast(Status::Count); s++) { + Status new_status = static_cast(s); + // check if an adoption from the agents status is possible + auto adoption_rate = m_model->adoption_rate(agent, new_status); + if (adoption_rate > 0) { + // add rate and corresponding event + m_current_rates.push_back(adoption_rate); + m_current_events.push_back({agent, new_status}); + } + } + } + } + + ScalarType m_t0, m_dt; + std::unique_ptr m_model; + std::vector m_current_rates; + std::vector m_current_events; // contains an event corresponding to each rate in m_current_rates + mio::TimeSeries m_result; +}; + +} // namespace mio + +#endif diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index e13b630b6e..32a960baec 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -47,6 +47,7 @@ set(TESTSOURCES test_contact_matrix.cpp test_type_safe.cpp test_custom_index_array.cpp + test_d_abm_model.cpp test_flows.cpp test_parameter_set.cpp test_matrix_shape.cpp @@ -99,7 +100,7 @@ endif() add_executable(memilio-test ${TESTSOURCES}) target_include_directories(memilio-test PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) -target_link_libraries(memilio-test PRIVATE memilio ode_secir ode_seir ode_secirvvs ode_secirts ode_seair ide_seir ide_secir lct_secir glct_secir abm gtest_main AD::AD) +target_link_libraries(memilio-test PRIVATE memilio ode_secir ode_seir ode_secirvvs ode_secirts ode_seair ide_seir ide_secir lct_secir glct_secir abm gtest_main AD::AD dabm) target_compile_options(memilio-test PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) # make unit tests find the test data files diff --git a/cpp/tests/test_d_abm_model.cpp b/cpp/tests/test_d_abm_model.cpp new file mode 100644 index 0000000000..c860526013 --- /dev/null +++ b/cpp/tests/test_d_abm_model.cpp @@ -0,0 +1,123 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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/parameters.h" +#include "d_abm/quad_well.h" +#include "d_abm/model.h" +#include "d_abm/simulation.h" +#include "memilio/utils/random_number_generator.h" +#include +#include +#include + +enum class InfectionState +{ + S, + E, + C, + I, + R, + D, + Count + +}; + +TEST(TestQuadWell, well) +{ + auto well = well_index(Eigen::Vector2d{-1, 1}); + EXPECT_EQ(well, size_t(0)); +} + +TEST(TestQuadWell, adopt) +{ + QuadWellModel::Agent agent{Eigen::Vector2d{-1, 1}, InfectionState::S}; + QuadWellModel qw({agent}, {}); + EXPECT_EQ(agent.status, InfectionState::S); + qw.adopt(agent, InfectionState::E); + EXPECT_EQ(agent.status, InfectionState::E); +} + +TEST(TestQuadWell, adoptionRate) +{ + QuadWellModel::Agent a1{Eigen::Vector2d{-1, 1}, InfectionState::S}; + QuadWellModel::Agent a2{Eigen::Vector2d{-1.2, 1}, InfectionState::I}; + QuadWellModel qw({a1, a2}, {{InfectionState::S, + InfectionState::E, + mio::dabm::Region(0), + 0.1, + {InfectionState::C, InfectionState::I}, + {1, 0.5}}}); + EXPECT_EQ(qw.adoption_rate(a1, InfectionState::E), 0.025); +} + +TEST(TestQuadWell, move) +{ + QuadWellModel::Agent a1{Eigen::Vector2d{-1.2, 1}, InfectionState::S}; + QuadWellModel::Agent a2{Eigen::Vector2d{-1.2, 1}, InfectionState::I}; + QuadWellModel qw({a1, a2}, {}, 0.1, 0., {InfectionState::I}); + qw.move(0, 0.1, a1); + qw.move(0, 0.1, a2); + EXPECT_EQ(a1.position[0], -0.9888); + EXPECT_EQ(a1.position[1], 1); + EXPECT_EQ(a2.position[0], -1.2); + EXPECT_EQ(a2.position[1], 1); +} + +TEST(TestQuadWell, setNonMovingRegions) +{ + QuadWellModel::Agent a1{Eigen::Vector2d{-1.2, 1}, InfectionState::S}; + QuadWellModel qw({a1}, {}); + qw.set_non_moving_regions({0}); + qw.move(0, 0.1, a1); + EXPECT_EQ(a1.position[0], -1.2); + EXPECT_EQ(a1.position[1], 1); +} + +TEST(TestDABMSimulation, Advance) +{ + using Model = mio::dabm::Model>; + auto& pos_rng = mio::UniformDistribution::get_instance(); + auto& sta_rng = mio::DiscreteDistribution::get_instance(); + std::vector pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.}; + std::vector agents(50); + for (auto& a : agents) { + a.position = + Eigen::Vector2d{pos_rng(mio::thread_local_rng(), -2., 2.), pos_rng(mio::thread_local_rng(), -2., 2.)}; + a.status = static_cast(sta_rng(mio::thread_local_rng(), pop_dist)); + } + std::vector> 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., {}, {}}); + } + Model model(agents, adoption_rates, 0.4, 0.4, {InfectionState::D}); + auto sim = mio::Simulation(model, 0.0, 0.1); + sim.advance(30.); + EXPECT_EQ(sim.get_result().get_last_time(), 30.); +} From 294b6014a6df73a87b1357ecddc2f415e16ca7ed Mon Sep 17 00:00:00 2001 From: jubicker Date: Thu, 12 Dec 2024 14:04:08 +0100 Subject: [PATCH 2/8] read me diffusive abm --- cpp/models/README.md | 14 +++++++++++++- cpp/models/d_abm/README.md | 14 ++++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) create mode 100644 cpp/models/d_abm/README.md diff --git a/cpp/models/README.md b/cpp/models/README.md index 7e49717f55..91fb146097 100644 --- a/cpp/models/README.md +++ b/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. diff --git a/cpp/models/d_abm/README.md b/cpp/models/d_abm/README.md new file mode 100644 index 0000000000..c78deec0fa --- /dev/null +++ b/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. From 9977cc8ce601f7bbaeffa9e8e139a6f5c08fe005 Mon Sep 17 00:00:00 2001 From: jubicker Date: Fri, 13 Dec 2024 11:16:28 +0100 Subject: [PATCH 3/8] add smm --- cpp/CMakeLists.txt | 1 + cpp/models/smm/CMakeLists.txt | 13 +++ cpp/models/smm/model.cpp | 29 +++++ cpp/models/smm/model.h | 89 ++++++++++++++ cpp/models/smm/parameters.h | 85 ++++++++++++++ cpp/models/smm/simulation.cpp | 24 ++++ cpp/models/smm/simulation.h | 205 +++++++++++++++++++++++++++++++++ cpp/tests/CMakeLists.txt | 3 +- cpp/tests/test_d_abm_model.cpp | 2 +- 9 files changed, 449 insertions(+), 2 deletions(-) create mode 100644 cpp/models/smm/CMakeLists.txt create mode 100644 cpp/models/smm/model.cpp create mode 100644 cpp/models/smm/model.h create mode 100644 cpp/models/smm/parameters.h create mode 100644 cpp/models/smm/simulation.cpp create mode 100644 cpp/models/smm/simulation.h diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 9228d92cc8..985c1d7d2c 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -164,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) diff --git a/cpp/models/smm/CMakeLists.txt b/cpp/models/smm/CMakeLists.txt new file mode 100644 index 0000000000..74996b370f --- /dev/null +++ b/cpp/models/smm/CMakeLists.txt @@ -0,0 +1,13 @@ +add_library(smm + parameters.h + model.h + model.cpp + simulation.h + simulation.cpp +) +target_link_libraries(smm PUBLIC memilio) +target_include_directories(smm PUBLIC + $ + $ +) +target_compile_options(smm PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/smm/model.cpp b/cpp/models/smm/model.cpp new file mode 100644 index 0000000000..6ffbf30d1f --- /dev/null +++ b/cpp/models/smm/model.cpp @@ -0,0 +1,29 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding +* +* Contact: Martin J. Kuehn +* +* 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/model.h" + +namespace mio +{ +namespace smm +{ + +} // namespace smm +} // namespace mio diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h new file mode 100644 index 0000000000..c302b469ee --- /dev/null +++ b/cpp/models/smm/model.h @@ -0,0 +1,89 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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_SMM_MODEL_H +#define MIO_SMM_MODEL_H + +#include "memilio/config.h" +#include "smm/parameters.h" +#include "memilio/compartments/compartmentalmodel.h" +#include "memilio/epidemiology/populations.h" + +namespace mio +{ +namespace smm +{ +template +class Model : public mio::CompartmentalModel, + ParametersBase> +{ + using Base = mio::CompartmentalModel, + ParametersBase>; + +public: + Model() + : Base(typename Base::Populations({static_cast(regions), Status::Count}, 0.), + typename Base::ParameterSet()) + { + } + + ScalarType evaluate(const AdoptionRate& rate, const Eigen::VectorXd& x) const + { + assert(rate.influences.size() == rate.factors.size()); + const auto& pop = this->populations; + const auto source = pop.get_flat_index({rate.region, rate.from}); + // determine order and calculate rate + if (rate.influences.size() == 0) { // first order adoption + return rate.factor * x[source]; + } + else { // second order adoption + ScalarType N = 0; + for (size_t s = 0; s < static_cast(Status::Count); ++s) { + N += x[pop.get_flat_index({rate.region, Status(s)})]; + } + // accumulate influences + ScalarType influences = 0.0; + for (size_t i = 0; i < rate.influences.size(); i++) { + influences += rate.factors[i] * x[pop.get_flat_index({rate.region, rate.influences[i]})]; + } + return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; + } + } + + ScalarType evaluate(const TransitionRate& rate, const Eigen::VectorXd& x) const + { + const auto source = this->populations.get_flat_index({rate.from, rate.status}); + return rate.factor * x[source]; + } + + mio::RandomNumberGenerator& get_rng() + { + return m_rng; + } + +private: + mio::RandomNumberGenerator m_rng; +}; + +} //namespace smm + +} // namespace mio + +#endif diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h new file mode 100644 index 0000000000..8103e8a310 --- /dev/null +++ b/cpp/models/smm/parameters.h @@ -0,0 +1,85 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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_SMM_PARAMETERS_H +#define MIO_SMM_PARAMETERS_H + +#include "memilio/config.h" +#include "memilio/utils/index.h" +#include "memilio/utils/parameter_set.h" + +namespace mio +{ +namespace smm +{ + +struct Region : public mio::Index { + Region(const size_t num_regions) + : mio::Index(num_regions) + { + } +}; + +// the AdoptionRate is considered to be of second order, if there are any "influences" with corresponding "factors". +// "from" is always an influence, scaled by "factor". +template +struct AdoptionRate { + Status from; // i + Status to; // j + Region region; // k + ScalarType factor; // gammahat_{ij}^k + std::vector influences; + std::vector factors; +}; + +template +struct AdoptionRates { + using Type = std::vector>; + const static std::string name() + { + return "AdoptionRates"; + } +}; + +template +struct TransitionRate { + Status status; // i + Region from; // k + Region to; // l + ScalarType factor; // lambda_i^{kl} +}; + +template +struct TransitionRates { + using Type = std::vector>; + const static std::string name() + { + return "TransitionRates"; + } +}; + +template +using ParametersBase = mio::ParameterSet, TransitionRates>; + +} // namespace smm + +} // namespace mio + +#endif diff --git a/cpp/models/smm/simulation.cpp b/cpp/models/smm/simulation.cpp new file mode 100644 index 0000000000..72280a646d --- /dev/null +++ b/cpp/models/smm/simulation.cpp @@ -0,0 +1,24 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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" + +namespace mio +{ +} // namespace mio diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h new file mode 100644 index 0000000000..dc5c1f5180 --- /dev/null +++ b/cpp/models/smm/simulation.h @@ -0,0 +1,205 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: René Schmieding, Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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_SMM_SIMULATION_H +#define MIO_SMM_SIMULATION_H + +#include "memilio/config.h" +#include "smm/model.h" +#include "smm/parameters.h" +#include "memilio/compartments/simulation.h" + +namespace mio +{ +template +class Simulation> +{ +public: +public: + using Model = smm::Model; + + /** + * @brief setup the simulation with an ODE solver + * @param[in] model: An instance of a compartmental model + * @param[in] t0 start time + * @param[in] dt initial step size of integration + */ + Simulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 1.) + : m_dt(dt) + , m_model(std::make_unique(model)) + , m_result(t0, m_model->get_initial_values()) + , m_internal_time(adoption_rates().size() + transition_rates().size(), t0) + , m_tp_next_event(adoption_rates().size() + transition_rates().size(), t0) + , m_waiting_times(adoption_rates().size() + transition_rates().size(), 0) + , m_current_rates(adoption_rates().size() + transition_rates().size(), 0) + { + assert(dt > 0); + assert(m_waiting_times.size() > 0); + assert(std::all_of(adoption_rates().begin(), adoption_rates().end(), [](auto&& r) { + return static_cast(r.region) < regions; + })); + assert(std::all_of(transition_rates().begin(), transition_rates().end(), [](auto&& r) { + return static_cast(r.from) < regions && static_cast(r.to) < regions; + })); + // initialize (internal) next event times by random values + for (size_t i = 0; i < m_tp_next_event.size(); i++) { + m_tp_next_event[i] += mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); + } + } + + /** + * @brief advance simulation to tmax + * tmax must be greater than get_result().get_last_time_point() + * @param tmax next stopping point of simulation + */ + Eigen::Ref advance(ScalarType tmax) + { + update_current_rates_and_waiting_times(); + size_t next_event = determine_next_event(); // index of the next event + ScalarType current_time = m_result.get_last_time(); + // set in the past to add a new time point immediately + ScalarType last_result_time = current_time - m_dt; + // iterate over time + while (current_time + m_waiting_times[next_event] < tmax) { + // update time + current_time += m_waiting_times[next_event]; + // regularily save current state in m_results + if (current_time > last_result_time + m_dt) { + last_result_time = current_time; + m_result.add_time_point(current_time); + // copy from the previous last value + m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; + } + // decide event type by index and perform it + if (next_event < adoption_rates().size()) { + // perform adoption event + const auto& rate = adoption_rates()[next_event]; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1; + m_model->populations[{rate.region, rate.from}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1; + m_model->populations[{rate.region, rate.to}] += 1; + } + else { + // perform transition event + const auto& rate = transition_rates()[next_event - adoption_rates().size()]; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1; + m_model->populations[{rate.from, rate.status}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1; + m_model->populations[{rate.to, rate.status}] += 1; + } + // update internal times + for (size_t i = 0; i < m_internal_time.size(); i++) { + m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; + } + // draw new "next event" time for the ocuured event + m_tp_next_event[next_event] += + mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); + // precalculate next event + update_current_rates_and_waiting_times(); + next_event = determine_next_event(); + } + // copy last result, if no event occurs between last_result_time and tmax + if (last_result_time < tmax) { + const auto& val = m_result.get_last_value(); + m_result.add_time_point(tmax, val); + } + return m_result.get_last_value(); + } + + /** + * @brief get_result returns the final simulation result + * @return a TimeSeries to represent the final simulation result + */ + TimeSeries& get_result() + { + return m_result; + } + + /** + * @brief get_result returns the final simulation result + * @return a TimeSeries to represent the final simulation result + */ + const TimeSeries& get_result() const + { + return m_result; + } + + /** + * @brief returns the simulation model used in simulation + */ + const Model& get_model() const + { + return *m_model; + } + + /** + * @brief returns the simulation model used in simulation + */ + Model& get_model() + { + return *m_model; + } + +private: + inline constexpr const typename smm::TransitionRates::Type& transition_rates() + { + return m_model->parameters.template get>(); + } + inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() + { + return m_model->parameters.template get>(); + } + + inline void update_current_rates_and_waiting_times() + { + size_t i = 0; // shared index for iterating both rates + for (const auto& rate : adoption_rates()) { + m_current_rates[i] = m_model->evaluate(rate, m_result.get_last_value()); + m_waiting_times[i] = (m_current_rates[i] > 0) + ? (m_tp_next_event[i] - m_internal_time[i]) / m_current_rates[i] + : std::numeric_limits::max(); + i++; + } + for (const auto& rate : transition_rates()) { + m_current_rates[i] = m_model->evaluate(rate, m_result.get_last_value()); + m_waiting_times[i] = (m_current_rates[i] > 0) + ? (m_tp_next_event[i] - m_internal_time[i]) / m_current_rates[i] + : std::numeric_limits::max(); + i++; + } + } + inline size_t determine_next_event() + { + return std::distance(m_waiting_times.begin(), std::min_element(m_waiting_times.begin(), m_waiting_times.end())); + } + + ScalarType m_dt; + std::unique_ptr m_model; + mio::TimeSeries m_result; + + std::vector m_internal_time; // internal times of all poisson processes (aka T_k) + std::vector m_tp_next_event; // internal time points of next event i after m_internal[i] (aka P_k) + std::vector m_waiting_times; // external times between m_internal_time and m_tp_next_event + std::vector m_current_rates; // current values of both types of rates +}; + +} // namespace mio + +#endif diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 32a960baec..a8a336c476 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -70,6 +70,7 @@ set(TESTSOURCES test_lct_initializer_flows.cpp test_glct_secir.cpp test_ad.cpp + test_smm_model.cpp abm_helpers.h abm_helpers.cpp actions.h @@ -100,7 +101,7 @@ endif() add_executable(memilio-test ${TESTSOURCES}) target_include_directories(memilio-test PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) -target_link_libraries(memilio-test PRIVATE memilio ode_secir ode_seir ode_secirvvs ode_secirts ode_seair ide_seir ide_secir lct_secir glct_secir abm gtest_main AD::AD dabm) +target_link_libraries(memilio-test PRIVATE memilio ode_secir ode_seir ode_secirvvs ode_secirts ode_seair ide_seir ide_secir lct_secir glct_secir abm gtest_main AD::AD dabm smm) target_compile_options(memilio-test PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) # make unit tests find the test data files diff --git a/cpp/tests/test_d_abm_model.cpp b/cpp/tests/test_d_abm_model.cpp index c860526013..356eb70bb4 100644 --- a/cpp/tests/test_d_abm_model.cpp +++ b/cpp/tests/test_d_abm_model.cpp @@ -90,7 +90,7 @@ TEST(TestQuadWell, setNonMovingRegions) EXPECT_EQ(a1.position[1], 1); } -TEST(TestDABMSimulation, Advance) +TEST(TestDABMSimulation, advance) { using Model = mio::dabm::Model>; auto& pos_rng = mio::UniformDistribution::get_instance(); From 0b0ebe3acbeaa7446e16b1d6cf77ae6cc0131f95 Mon Sep 17 00:00:00 2001 From: jubicker Date: Fri, 13 Dec 2024 11:16:40 +0100 Subject: [PATCH 4/8] add tests for smm --- cpp/tests/test_smm_model.cpp | 141 +++++++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 cpp/tests/test_smm_model.cpp diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp new file mode 100644 index 0000000000..3ed7226e88 --- /dev/null +++ b/cpp/tests/test_smm_model.cpp @@ -0,0 +1,141 @@ +/* +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* +* Authors: Julia Bicker +* +* Contact: Martin J. Kuehn +* +* 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 +#include +#include "smm/model.h" +#include "smm/parameters.h" +#include "smm/simulation.h" +#include + +enum class InfectionState +{ + S, + E, + C, + I, + R, + D, + Count + +}; + +TEST(TestSMM, evaluateAdoptionRate) +{ + using Model = mio::smm::Model<1, InfectionState>; + + Model model; + + std::vector> adoption_rates; + adoption_rates.push_back({InfectionState::S, + InfectionState::E, + mio::smm::Region(0), + 0.1, + {InfectionState::C, InfectionState::I}, + {1, 0.5}}); + adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::smm::Region(0), 0.2, {}, {}}); + + model.populations[{mio::smm::Region(0), InfectionState::S}] = 50; + model.populations[{mio::smm::Region(0), InfectionState::E}] = 10; + model.populations[{mio::smm::Region(0), InfectionState::C}] = 5; + model.populations[{mio::smm::Region(0), InfectionState::I}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::R}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::D}] = 0; + + EXPECT_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()), 5. / 13.); + EXPECT_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 2.); +} + +TEST(TestSMM, evaluateTransitionRate) +{ + using Model = mio::smm::Model<2, InfectionState>; + + Model model; + + model.populations[{mio::smm::Region(0), InfectionState::S}] = 50; + model.populations[{mio::smm::Region(0), InfectionState::E}] = 10; + model.populations[{mio::smm::Region(0), InfectionState::C}] = 5; + model.populations[{mio::smm::Region(0), InfectionState::I}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::R}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::D}] = 0; + + model.populations[{mio::smm::Region(1), InfectionState::S}] = 55; + model.populations[{mio::smm::Region(1), InfectionState::E}] = 10; + model.populations[{mio::smm::Region(1), InfectionState::C}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::I}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::R}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::D}] = 0; + + std::vector> transition_rates; + transition_rates.push_back({InfectionState::S, mio::smm::Region(0), mio::smm::Region(1), 0.01}); + transition_rates.push_back({InfectionState::E, mio::smm::Region(1), mio::smm::Region(0), 0.1}); + + EXPECT_EQ(model.evaluate(transition_rates[0], model.populations.get_compartments()), 0.5); + EXPECT_EQ(model.evaluate(transition_rates[1], model.populations.get_compartments()), 1.); +} + +TEST(TestSMMSimulation, advance) +{ + using Model = mio::smm::Model<2, InfectionState>; + + Model model; + + model.populations[{mio::smm::Region(0), InfectionState::S}] = 50; + model.populations[{mio::smm::Region(0), InfectionState::E}] = 10; + model.populations[{mio::smm::Region(0), InfectionState::C}] = 5; + model.populations[{mio::smm::Region(0), InfectionState::I}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::R}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::D}] = 0; + + model.populations[{mio::smm::Region(1), InfectionState::S}] = 55; + model.populations[{mio::smm::Region(1), InfectionState::E}] = 10; + model.populations[{mio::smm::Region(1), InfectionState::C}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::I}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::R}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::D}] = 0; + + std::vector> adoption_rates; + std::vector> transition_rates; + for (size_t region = 0; region < 2; ++region) { + adoption_rates.push_back({InfectionState::S, + InfectionState::E, + mio::smm::Region(region), + 0.1, + {InfectionState::C, InfectionState::I}, + {1, 0.5}}); + adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::smm::Region(region), 1.0 / 5., {}, {}}); + adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::smm::Region(region), 0.2 / 3., {}, {}}); + adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::smm::Region(region), 0.8 / 3., {}, {}}); + adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::smm::Region(region), 0.99 / 5., {}, {}}); + adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::smm::Region(region), 0.01 / 5., {}, {}}); + } + + for (size_t s = 0; s < static_cast(InfectionState::D); ++s) { + transition_rates.push_back({InfectionState(s), mio::smm::Region(0), mio::smm::Region(1), 0.01}); + transition_rates.push_back({InfectionState(s), mio::smm::Region(1), mio::smm::Region(0), 0.01}); + } + + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; + + auto sim = mio::Simulation(model, 0.0, 0.1); + sim.advance(30.); + EXPECT_EQ(sim.get_result().get_last_time(), 30.); +} From 2b95423e3b2403bae3803f67019c75161ea4e1c6 Mon Sep 17 00:00:00 2001 From: jubicker Date: Fri, 13 Dec 2024 11:41:35 +0100 Subject: [PATCH 5/8] smm readme --- cpp/models/smm/README.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 cpp/models/smm/README.md diff --git a/cpp/models/smm/README.md b/cpp/models/smm/README.md new file mode 100644 index 0000000000..73760775a0 --- /dev/null +++ b/cpp/models/smm/README.md @@ -0,0 +1,12 @@ +# Stochastic Metapopulation Model + +The stochastic metapopulation model uses a Markov process to simulate disease dynamics. Agents have an infection state and a position which is given by the region index the agent is in, hence agents in the same region have the same position. 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 as stochastic jumps with independent Poisson processes given by adoption rate functions. As agents' positions are given by their subregion index, there is no movement within one subregion. Movement between subregions is also modeled as stochastic jumps with independent Poisson processes given by transition rate functions. Gillespie's direct methode (stochastic simulation algorithm) is used for simulation. + +## Simulation + +At the beginning of the simulation, the waiting times for all events are drawn. Then the time is advanced until the time point of the next event - which can be a spatial transition or an infection state adoption - and the event takes places. The waiting times of the other events are updated and a new waiting time for the event that just happend is drawn. The simulation saves the system state in discrete time steps. From 8807c550d23bfd5c8bf398ad6fe289235543bc5b Mon Sep 17 00:00:00 2001 From: jubicker Date: Fri, 13 Dec 2024 12:36:49 +0100 Subject: [PATCH 6/8] msvc is a special snowflake --- cpp/tests/test_d_abm_model.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/tests/test_d_abm_model.cpp b/cpp/tests/test_d_abm_model.cpp index 356eb70bb4..1bcd8a2b8a 100644 --- a/cpp/tests/test_d_abm_model.cpp +++ b/cpp/tests/test_d_abm_model.cpp @@ -94,7 +94,7 @@ TEST(TestDABMSimulation, advance) { using Model = mio::dabm::Model>; auto& pos_rng = mio::UniformDistribution::get_instance(); - auto& sta_rng = mio::DiscreteDistribution::get_instance(); + auto& sta_rng = mio::DiscreteDistribution::get_instance(); std::vector pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.}; std::vector agents(50); for (auto& a : agents) { From 26f72588f955649f83f87c5d87f29a5894cd9e5f Mon Sep 17 00:00:00 2001 From: jubicker Date: Fri, 13 Dec 2024 15:06:52 +0100 Subject: [PATCH 7/8] examples smm and dabm --- cpp/examples/CMakeLists.txt | 27 +++++++---- cpp/examples/d_abm.cpp | 86 ++++++++++++++++++++++++++++++++ cpp/examples/smm.cpp | 97 +++++++++++++++++++++++++++++++++++++ 3 files changed, 200 insertions(+), 10 deletions(-) create mode 100644 cpp/examples/d_abm.cpp create mode 100644 cpp/examples/smm.cpp diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 8fd1f51c70..c8b63c0d9d 100644 --- a/cpp/examples/CMakeLists.txt +++ b/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}) @@ -127,6 +126,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) @@ -135,13 +142,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) @@ -163,7 +170,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() diff --git a/cpp/examples/d_abm.cpp b/cpp/examples/d_abm.cpp new file mode 100644 index 0000000000..65946685a8 --- /dev/null +++ b/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 +* +* 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 + +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>; + std::vector agents(1000); + //Random variables for initialization of agents' position and infection state + auto& pos_rng = mio::UniformDistribution::get_instance(); + auto& sta_rng = mio::DiscreteDistribution::get_instance(); + //Infection state distribution + std::vector 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(sta_rng(mio::thread_local_rng(), pop_dist)); + } + + //Set adoption rates + std::vector> 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(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 "}); +} diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp new file mode 100644 index 0000000000..78272fb306 --- /dev/null +++ b/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 +* +* 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; + + 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> adoption_rates; + std::vector> 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(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>() = adoption_rates; + model.parameters.get>() = transition_rates; + + double dt = 0.1; + double tmax = 30.; + + auto sim = mio::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 "}); +} From 428cffab99b30f092e692b3364ad0e39dc141dae Mon Sep 17 00:00:00 2001 From: jubicker Date: Wed, 18 Dec 2024 16:46:11 +0100 Subject: [PATCH 8/8] advance tests for smm and d_abm simulations --- cpp/models/d_abm/quad_well.h | 1 - cpp/models/smm/simulation.h | 2 +- cpp/tests/test_d_abm_model.cpp | 41 ++++++++++++----- cpp/tests/test_smm_model.cpp | 80 +++++++++++++++++++++++----------- 4 files changed, 86 insertions(+), 38 deletions(-) diff --git a/cpp/models/d_abm/quad_well.h b/cpp/models/d_abm/quad_well.h index 1f78ef3b95..5abfe2375d 100644 --- a/cpp/models/d_abm/quad_well.h +++ b/cpp/models/d_abm/quad_well.h @@ -63,7 +63,6 @@ class QuadWellModel } for (auto& r : rates) { m_adoption_rates.emplace(std::forward_as_tuple(r.region, r.from, r.to), r); - //m_adoption_rates[{r.region, r.from, r.to}] = r; } } diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index dc5c1f5180..0d30abd49f 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -108,7 +108,7 @@ class Simulation> for (size_t i = 0; i < m_internal_time.size(); i++) { m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; } - // draw new "next event" time for the ocuured event + // draw new "next event" time for the occured event m_tp_next_event[next_event] += mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); // precalculate next event diff --git a/cpp/tests/test_d_abm_model.cpp b/cpp/tests/test_d_abm_model.cpp index 1bcd8a2b8a..fc9105adbf 100644 --- a/cpp/tests/test_d_abm_model.cpp +++ b/cpp/tests/test_d_abm_model.cpp @@ -23,6 +23,7 @@ #include "d_abm/model.h" #include "d_abm/simulation.h" #include "memilio/utils/random_number_generator.h" +#include "abm_helpers.h" #include #include #include @@ -92,16 +93,11 @@ TEST(TestQuadWell, setNonMovingRegions) TEST(TestDABMSimulation, advance) { - using Model = mio::dabm::Model>; - auto& pos_rng = mio::UniformDistribution::get_instance(); - auto& sta_rng = mio::DiscreteDistribution::get_instance(); - std::vector pop_dist{0.98, 0.01, 0.005, 0.005, 0., 0.}; - std::vector agents(50); - for (auto& a : agents) { - a.position = - Eigen::Vector2d{pos_rng(mio::thread_local_rng(), -2., 2.), pos_rng(mio::thread_local_rng(), -2., 2.)}; - a.status = static_cast(sta_rng(mio::thread_local_rng(), pop_dist)); - } + using testing::Return; + using Model = mio::dabm::Model>; + QuadWellModel::Agent a1{Eigen::Vector2d{-1, 1}, InfectionState::S}; + QuadWellModel::Agent a2{Eigen::Vector2d{-1, 1}, InfectionState::R}; + QuadWellModel::Agent a3{Eigen::Vector2d{-1, 1}, InfectionState::I}; std::vector> adoption_rates; for (size_t region = 0; region < 4; ++region) { adoption_rates.push_back({InfectionState::S, @@ -116,8 +112,31 @@ TEST(TestDABMSimulation, advance) 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., {}, {}}); } - Model model(agents, adoption_rates, 0.4, 0.4, {InfectionState::D}); + Model model({a1, a2, a3}, adoption_rates, 0.4, 0.0, {InfectionState::D}); auto sim = mio::Simulation(model, 0.0, 0.1); + //Setup so first adoption event will be in second time step + ScopedMockDistribution>>> + mock_exponential_dist; + EXPECT_CALL(mock_exponential_dist.get_mock(), invoke) + .Times(testing::AtLeast(1)) + .WillOnce(Return(0.0226)) + .WillRepeatedly(testing::Return(1.0)); + //Setup so first adoption event will be the transition of a3 from I to R + ScopedMockDistribution>>> mock_discrete_dist; + EXPECT_CALL(mock_discrete_dist.get_mock(), invoke).Times(1).WillOnce(Return(size_t(1))); sim.advance(30.); + //Initial positions and infection state of all agents + EXPECT_EQ(sim.get_result().get_value(0)[static_cast(InfectionState::S)], 1); + EXPECT_EQ(sim.get_result().get_value(0)[static_cast(InfectionState::I)], 1); + EXPECT_EQ(sim.get_result().get_value(0)[static_cast(InfectionState::R)], 1); + //In the first time step no adoption event happens + EXPECT_EQ(sim.get_result().get_value(1)[static_cast(InfectionState::S)], 1); + EXPECT_EQ(sim.get_result().get_value(1)[static_cast(InfectionState::I)], 1); + EXPECT_EQ(sim.get_result().get_value(1)[static_cast(InfectionState::R)], 1); + //a3 transitions from I to R in the second time step + EXPECT_EQ(sim.get_result().get_value(2)[static_cast(InfectionState::S)], 1); + EXPECT_EQ(sim.get_result().get_value(2)[static_cast(InfectionState::I)], 0); + EXPECT_EQ(sim.get_result().get_value(2)[static_cast(InfectionState::R)], 2); + //check whether simulation advances until the end EXPECT_EQ(sim.get_result().get_last_time(), 30.); } diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 3ed7226e88..4c4ccaff11 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -23,6 +23,7 @@ #include "smm/model.h" #include "smm/parameters.h" #include "smm/simulation.h" +#include "abm_helpers.h" #include enum class InfectionState @@ -93,49 +94,78 @@ TEST(TestSMM, evaluateTransitionRate) TEST(TestSMMSimulation, advance) { + using testing::Return; using Model = mio::smm::Model<2, InfectionState>; Model model; - model.populations[{mio::smm::Region(0), InfectionState::S}] = 50; - model.populations[{mio::smm::Region(0), InfectionState::E}] = 10; - model.populations[{mio::smm::Region(0), InfectionState::C}] = 5; - model.populations[{mio::smm::Region(0), InfectionState::I}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::S}] = 1; + model.populations[{mio::smm::Region(0), InfectionState::E}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::C}] = 0; + model.populations[{mio::smm::Region(0), InfectionState::I}] = 1; model.populations[{mio::smm::Region(0), InfectionState::R}] = 0; model.populations[{mio::smm::Region(0), InfectionState::D}] = 0; - model.populations[{mio::smm::Region(1), InfectionState::S}] = 55; - model.populations[{mio::smm::Region(1), InfectionState::E}] = 10; + model.populations[{mio::smm::Region(1), InfectionState::S}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::E}] = 0; model.populations[{mio::smm::Region(1), InfectionState::C}] = 0; model.populations[{mio::smm::Region(1), InfectionState::I}] = 0; - model.populations[{mio::smm::Region(1), InfectionState::R}] = 0; + model.populations[{mio::smm::Region(1), InfectionState::R}] = 1; model.populations[{mio::smm::Region(1), InfectionState::D}] = 0; std::vector> adoption_rates; std::vector> transition_rates; - for (size_t region = 0; region < 2; ++region) { - adoption_rates.push_back({InfectionState::S, - InfectionState::E, - mio::smm::Region(region), - 0.1, - {InfectionState::C, InfectionState::I}, - {1, 0.5}}); - adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::smm::Region(region), 1.0 / 5., {}, {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::smm::Region(region), 0.2 / 3., {}, {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::smm::Region(region), 0.8 / 3., {}, {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::smm::Region(region), 0.99 / 5., {}, {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::smm::Region(region), 0.01 / 5., {}, {}}); - } - - for (size_t s = 0; s < static_cast(InfectionState::D); ++s) { - transition_rates.push_back({InfectionState(s), mio::smm::Region(0), mio::smm::Region(1), 0.01}); - transition_rates.push_back({InfectionState(s), mio::smm::Region(1), mio::smm::Region(0), 0.01}); - } + + adoption_rates.push_back({InfectionState::S, + InfectionState::E, + mio::smm::Region(0), + 0.1, + {InfectionState::C, InfectionState::I}, + {1, 0.5}}); + adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::smm::Region(0), 1.0 / 5., {}, {}}); + adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::smm::Region(0), 0.2 / 3., {}, {}}); + adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::smm::Region(0), 0.8 / 3., {}, {}}); + adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::smm::Region(0), 0.99 / 5., {}, {}}); + adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::smm::Region(0), 0.01 / 5., {}, {}}); + + transition_rates.push_back({InfectionState::R, mio::smm::Region(1), mio::smm::Region(0), 0.01}); model.parameters.get>() = adoption_rates; model.parameters.get>() = transition_rates; + ScopedMockDistribution>>> + mock_exponential_dist; + EXPECT_CALL(mock_exponential_dist.get_mock(), invoke) + .Times(testing::AtLeast(7)) + .WillOnce(Return(0.5)) + .WillOnce(Return(0.5)) + .WillOnce(Return(0.5)) + .WillOnce(Return(0.5)) + .WillOnce(Return(0.0397)) + .WillOnce(Return(0.5)) + .WillOnce(Return(0.0031)) + .WillRepeatedly(testing::Return(1.0)); + auto sim = mio::Simulation(model, 0.0, 0.1); sim.advance(30.); + //initial values + EXPECT_EQ(sim.get_result().get_value(0)[static_cast(InfectionState::S)], 1); + EXPECT_EQ(sim.get_result().get_value(0)[static_cast(InfectionState::I)], 1); + EXPECT_EQ(sim.get_result().get_value( + 0)[static_cast(InfectionState::Count) + static_cast(InfectionState::R)], + 1); + //no event happens in first time step + EXPECT_GE(sim.get_result().get_time(1), 0.2); + //adoption from I to R + EXPECT_EQ(sim.get_result().get_value(1)[static_cast(InfectionState::S)], 1); + EXPECT_EQ(sim.get_result().get_value(1)[static_cast(InfectionState::I)], 0); + EXPECT_EQ(sim.get_result().get_value(1)[static_cast(InfectionState::R)], 1); + EXPECT_EQ(sim.get_result().get_value( + 1)[static_cast(InfectionState::Count) + static_cast(InfectionState::R)], + 1); + //spatial transition + EXPECT_EQ(sim.get_result().get_value(2)[static_cast(InfectionState::S)], 1); + EXPECT_EQ(sim.get_result().get_value(2)[static_cast(InfectionState::I)], 0); + EXPECT_EQ(sim.get_result().get_value(2)[static_cast(InfectionState::R)], 2); EXPECT_EQ(sim.get_result().get_last_time(), 30.); }