From 95c6c4a53baff7da10ace02f0a04e301fed98286 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 08:49:01 +0100 Subject: [PATCH 01/16] bump dependencies --- cpp/thirdparty/CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cpp/thirdparty/CMakeLists.txt b/cpp/thirdparty/CMakeLists.txt index 5fb089603f..90adecd337 100644 --- a/cpp/thirdparty/CMakeLists.txt +++ b/cpp/thirdparty/CMakeLists.txt @@ -1,12 +1,12 @@ # Versions of the bundled libraries # If you like to upgrade, just change the number -set(MEMILIO_EIGEN_VERSION "3.3.9") -set(MEMILIO_SPDLOG_VERSION "1.11.0") +set(MEMILIO_EIGEN_VERSION "3.4.0") +set(MEMILIO_SPDLOG_VERSION "1.15.0") set(MEMILIO_BOOST_VERSION "1.84.0") set(MEMILIO_MINIMAL_BOOST_VERSION "1.76.0") -set(MEMILIO_JSONCPP_VERSION "1.9.5") +set(MEMILIO_JSONCPP_VERSION "1.9.6") set(MEMILIO_RANDOM123_VERSION "v1.14.0") -set(MEMILIO_IPOPT_VERSION "3.14.12") +set(MEMILIO_IPOPT_VERSION "3.14.17") # Gperftools for profiling; must be first, so that libraries are included in the profile if(MEMILIO_ENABLE_PROFILING) From d1e1bf141cf98d181be2d8f2f36f9bd7a2fed00d Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 09:00:13 +0100 Subject: [PATCH 02/16] make use of boost helpers for resizing eigen types with new version --- cpp/memilio/math/stepper_wrapper.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/cpp/memilio/math/stepper_wrapper.h b/cpp/memilio/math/stepper_wrapper.h index 428a7c2770..41a6d818f9 100644 --- a/cpp/memilio/math/stepper_wrapper.h +++ b/cpp/memilio/math/stepper_wrapper.h @@ -24,6 +24,7 @@ #include "memilio/utils/logging.h" #include "boost/numeric/odeint/external/eigen/eigen_algebra.hpp" +#include "boost/numeric/odeint/external/eigen/eigen_resize.hpp" #include "boost/numeric/odeint/stepper/controlled_runge_kutta.hpp" #include "boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp" #include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" @@ -44,7 +45,7 @@ class ControlledStepperWrapper : public mio::IntegratorCore using Stepper = boost::numeric::odeint::controlled_runge_kutta< ControlledStepper, FP, Vector, FP, boost::numeric::odeint::vector_space_algebra, typename boost::numeric::odeint::operations_dispatcher>::operations_type, - boost::numeric::odeint::never_resizer>>; + boost::numeric::odeint::initially_resizer>>; static constexpr bool is_fsal_stepper = std::is_same_v; static_assert(!is_fsal_stepper, @@ -109,7 +110,6 @@ class ControlledStepperWrapper : public mio::IntegratorCore step_result = m_stepper.try_step( // reorder arguments of the DerivFunction f for the stepper [&](const Vector& x, Vector& dxds, FP s) { - dxds.resizeLike(x); // boost resizers cannot resize Eigen::Vector, hence we need to do that here f(x, s, dxds); }, m_yt, t, m_ytp1, dt); @@ -186,7 +186,7 @@ class ExplicitStepperWrapper : public mio::IntegratorCore public: using Stepper = ExplicitStepper, FP, Vector, FP, boost::numeric::odeint::vector_space_algebra, typename boost::numeric::odeint::operations_dispatcher>::operations_type, - boost::numeric::odeint::never_resizer>; + boost::numeric::odeint::initially_resizer>; /** * @brief Set up the integrator. @@ -213,7 +213,6 @@ class ExplicitStepperWrapper : public mio::IntegratorCore m_stepper.do_step( // reorder arguments of the DerivFunction f for the stepper [&](const Vector& x, Vector& dxds, FP s) { - dxds.resizeLike(x); // do_step calls sys with a vector of size 0 for some reason f(x, s, dxds); }, ytp1, t, dt); From 638e34a7584ae3d19104dbdd70f8f86fdf9f684b Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 09:31:31 +0100 Subject: [PATCH 03/16] generalize and simplify logging for AD types --- cpp/memilio/utils/logging.h | 46 +++++++++++++------------------------ 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/cpp/memilio/utils/logging.h b/cpp/memilio/utils/logging.h index b20b504cb0..ac3f79376b 100644 --- a/cpp/memilio/utils/logging.h +++ b/cpp/memilio/utils/logging.h @@ -127,38 +127,24 @@ inline void log(LogLevel level, spdlog::string_view_t fmt, const Args&... args) } // namespace mio -/** - * @brief The fmt::formatter class provides formatting capabilities for the ad::gt1s::type - * to help SPD log to handle it. - */ -template <> -struct fmt::formatter::type> { - constexpr auto parse(format_parse_context& ctx) -> decltype(ctx.begin()) - { - return ctx.end(); - } - template - auto format(const ad::gt1s::type& input, FormatContext& ctx) -> decltype(ctx.out()) - { - return format_to(ctx.out(), "{}", ad::value(input)); - } -}; +namespace ad +{ +namespace internal +{ /** - * @brief The fmt::formatter class provides formatting capabilities for the ad::ga1s::type - * to help SPD log to handle it. + * @brief Format AD types (like ad::gt1s) using their value for logging with spdlog. + * + * If derivative information is needed as well, use `ad::derivative(...)` or define a `fmt::formatter<...>`. */ -template <> -struct fmt::formatter::type> { - constexpr auto parse(format_parse_context& ctx) -> decltype(ctx.begin()) - { - return ctx.end(); - } - template - auto format(const ad::ga1s::type& input, FormatContext& ctx) -> decltype(ctx.out()) - { - return format_to(ctx.out(), "{}", ad::value(input)); - } -}; +template +FP format_as(const active_type& ad_type) +{ + // Note: the format_as function needs to be in the same namespace as the value it takes + return value(ad_type); +} + +} // namespace internal +} // namespace ad #endif // LOGGING_H From 806c8c1be8a1dc96f4ee2e66e3890d4292a27825 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 09:31:55 +0100 Subject: [PATCH 04/16] specify logging for UncertainValue --- cpp/memilio/utils/uncertain_value.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/cpp/memilio/utils/uncertain_value.h b/cpp/memilio/utils/uncertain_value.h index f9af2f4163..27d660c72e 100644 --- a/cpp/memilio/utils/uncertain_value.h +++ b/cpp/memilio/utils/uncertain_value.h @@ -227,6 +227,15 @@ class UncertainValue std::unique_ptr m_dist; }; +/** + * @brief Format UncertainValues using their value for logging with spdlog. + */ +template +double format_as(const UncertainValue& uv) +{ + return uv; +} + // gtest printer // TODO: should be extended when UncertainValue gets operator== that compares distributions as well template From 889fe5b66eed7bfc58db8e0627982d4c277065f2 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 09:42:00 +0100 Subject: [PATCH 05/16] adapt ABM caches to new resize behaviour of Eigen types --- cpp/memilio/config.h | 1 - cpp/memilio/utils/custom_index_array.h | 18 ++++++++++++++++-- cpp/models/abm/model.cpp | 8 +++++--- cpp/tests/abm_helpers.cpp | 12 +++++++----- 4 files changed, 28 insertions(+), 11 deletions(-) diff --git a/cpp/memilio/config.h b/cpp/memilio/config.h index 98c6509f61..6f813d737d 100644 --- a/cpp/memilio/config.h +++ b/cpp/memilio/config.h @@ -26,7 +26,6 @@ #define MIO_CONFIG_H #include "memilio/config_internal.h" -#include using ScalarType = double; diff --git a/cpp/memilio/utils/custom_index_array.h b/cpp/memilio/utils/custom_index_array.h index 1738f6298f..4c141b68c0 100644 --- a/cpp/memilio/utils/custom_index_array.h +++ b/cpp/memilio/utils/custom_index_array.h @@ -235,8 +235,9 @@ class CustomIndexArray } /** - * Resize all dimensions. - * @param new_dims new dimensions. + * @brief Resize all dimensions. + ' Note that when increasing the overall size, new values may be uninitialized. + * @param new_dims New dimensions. */ void resize(Index new_dims) { @@ -245,6 +246,19 @@ class CustomIndexArray m_y.conservativeResize(m_numel); } + /** + * @brief Resize all dimensions, destroying all values. + * This Version of resize should only be used when the CustomIndexArray contains non-movable and non-copyable + * values, like atomics. New entries are all default initialized. + * @param new_dims New dimensions. + */ + void resize_destructive(Index new_dims) + { + m_dimensions = new_dims; + m_numel = product(m_dimensions); + m_y.resize(m_numel); + } + /** * Resize a single dimension. * @param new dimension. diff --git a/cpp/models/abm/model.cpp b/cpp/models/abm/model.cpp index 4c8ba1f7f5..54a1e8f98f 100755 --- a/cpp/models/abm/model.cpp +++ b/cpp/models/abm/model.cpp @@ -225,9 +225,11 @@ void Model::build_exposure_caches() m_contact_exposure_rates_cache.resize(num_locations); PRAGMA_OMP(taskloop) for (size_t i = 0; i < num_locations; i++) { - m_air_exposure_rates_cache[i].resize({CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count}); - m_contact_exposure_rates_cache[i].resize({CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count, - AgeGroup(parameters.get_num_groups())}); + m_air_exposure_rates_cache[i].resize_destructive( + {CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count}); + m_contact_exposure_rates_cache[i].resize_destructive({CellIndex(m_locations[i].get_cells().size()), + VirusVariant::Count, + AgeGroup(parameters.get_num_groups())}); } // implicit taskloop barrier m_are_exposure_caches_valid = false; m_exposure_caches_need_rebuild = false; diff --git a/cpp/tests/abm_helpers.cpp b/cpp/tests/abm_helpers.cpp index 06e3fdf63c..3ff3206841 100644 --- a/cpp/tests/abm_helpers.cpp +++ b/cpp/tests/abm_helpers.cpp @@ -35,8 +35,8 @@ mio::abm::Person make_test_person(mio::RandomNumberGenerator& rng, mio::abm::Loc return p; } -mio::abm::PersonId add_test_person(mio::abm::Model& model, mio::abm::LocationId loc_id, - mio::AgeGroup age, mio::abm::InfectionState infection_state, mio::abm::TimePoint t) +mio::abm::PersonId add_test_person(mio::abm::Model& model, mio::abm::LocationId loc_id, mio::AgeGroup age, + mio::abm::InfectionState infection_state, mio::abm::TimePoint t) { return model.add_person( make_test_person(model.get_rng(), model.get_location(loc_id), age, infection_state, t, model.parameters)); @@ -49,14 +49,16 @@ void interact_testing(mio::abm::PersonalRandomNumberGenerator& personal_rng, mio { // allocate and initialize air exposures with 0 mio::abm::AirExposureRates local_air_exposure; - local_air_exposure.resize({mio::abm::CellIndex(location.get_cells().size()), mio::abm::VirusVariant::Count}); + local_air_exposure.resize_destructive( + {mio::abm::CellIndex(location.get_cells().size()), mio::abm::VirusVariant::Count}); std::for_each(local_air_exposure.begin(), local_air_exposure.end(), [](auto& r) { r = 0.0; }); // allocate and initialize contact exposures with 0 mio::abm::ContactExposureRates local_contact_exposure; - local_contact_exposure.resize({mio::abm::CellIndex(location.get_cells().size()), mio::abm::VirusVariant::Count, - mio::AgeGroup(global_parameters.get_num_groups())}); + local_contact_exposure.resize_destructive({mio::abm::CellIndex(location.get_cells().size()), + mio::abm::VirusVariant::Count, + mio::AgeGroup(global_parameters.get_num_groups())}); std::for_each(local_contact_exposure.begin(), local_contact_exposure.end(), [](auto& r) { r = 0.0; }); From c97f85daa849d800848e0ef69419de8d16d045c7 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 11:40:09 +0100 Subject: [PATCH 06/16] work around Eigen now supporting STL iterators by making is_container more strict --- cpp/memilio/io/io.h | 6 +++++- cpp/tests/test_abm_serialization.cpp | 13 +++++++------ 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/cpp/memilio/io/io.h b/cpp/memilio/io/io.h index fe3e651288..6cc1575f05 100644 --- a/cpp/memilio/io/io.h +++ b/cpp/memilio/io/io.h @@ -787,7 +787,11 @@ using compare_iterators_t = decltype(std::declval().begin() != std::de * @tparam C any type. */ template -using is_container = is_expression_valid; +using is_container = + conjunction, + std::is_constructible().begin()), decltype(std::declval().end())>, + // Eigen types may pass as container, but we want to handle them separately + negation, C>>>; /** * serialize an STL compatible container. diff --git a/cpp/tests/test_abm_serialization.cpp b/cpp/tests/test_abm_serialization.cpp index f704b24288..834abc1059 100644 --- a/cpp/tests/test_abm_serialization.cpp +++ b/cpp/tests/test_abm_serialization.cpp @@ -32,6 +32,7 @@ #include "models/abm/person.h" #include "models/abm/trip_list.h" #include "models/abm/model.h" +#include #include #ifdef MEMILIO_HAS_JSONCPP @@ -126,12 +127,12 @@ TEST(TestAbmSerialization, TestingScheme) unsigned i = 1; // counter s.t. members have different values Json::Value testing_criteria; - std::array ages_bits{}; // initialize to false - ages_bits[i++] = true; - testing_criteria["ages"]["bitset"] = mio::serialize_json(ages_bits).value(); - std::array inf_st_bits{}; // initialize to false - inf_st_bits[i++] = true; - testing_criteria["infection_states"]["bitset"] = mio::serialize_json(inf_st_bits).value(); + std::bitset ages_bits{}; // initialize to false + ages_bits[i++] = true; + testing_criteria["ages"] = mio::serialize_json(ages_bits).value(); + std::bitset<(size_t)mio::abm::InfectionState::Count> inf_st_bits{}; // initialize to false + inf_st_bits[i++] = true; + testing_criteria["infection_states"] = mio::serialize_json(inf_st_bits).value(); Json::Value test_parameters; test_parameters["sensitivity"] = mio::serialize_json(mio::UncertainValue{(double)i++}).value(); From d3b94d5b355dc24edf9b938fdbe0a63e21d021c0 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 11:52:22 +0100 Subject: [PATCH 07/16] unbump ipopt --- cpp/thirdparty/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/thirdparty/CMakeLists.txt b/cpp/thirdparty/CMakeLists.txt index 90adecd337..c528354823 100644 --- a/cpp/thirdparty/CMakeLists.txt +++ b/cpp/thirdparty/CMakeLists.txt @@ -6,7 +6,7 @@ set(MEMILIO_BOOST_VERSION "1.84.0") set(MEMILIO_MINIMAL_BOOST_VERSION "1.76.0") set(MEMILIO_JSONCPP_VERSION "1.9.6") set(MEMILIO_RANDOM123_VERSION "v1.14.0") -set(MEMILIO_IPOPT_VERSION "3.14.17") +set(MEMILIO_IPOPT_VERSION "3.14.12") # Gperftools for profiling; must be first, so that libraries are included in the profile if(MEMILIO_ENABLE_PROFILING) From 1537bed3744a7d9143dafb5b60dd89c88b7b04a8 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 15:31:24 +0100 Subject: [PATCH 08/16] include Eigen as System library to avoid errors --- cpp/thirdparty/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/thirdparty/CMakeLists.txt b/cpp/thirdparty/CMakeLists.txt index c528354823..4835a4d606 100644 --- a/cpp/thirdparty/CMakeLists.txt +++ b/cpp/thirdparty/CMakeLists.txt @@ -58,7 +58,7 @@ if(MEMILIO_USE_BUNDLED_EIGEN) endif() add_library(eigen INTERFACE) - target_include_directories(eigen INTERFACE ${eigen_SOURCE_DIR}) + target_include_directories(eigen SYSTEM INTERFACE ${eigen_SOURCE_DIR}) add_library(Eigen3::Eigen ALIAS eigen) else() find_package(Eigen3 ${MEMILIO_EIGEN_VERSION} REQUIRED NO_MODULE) From e9cca2c917e47f10de3d9233932f74a86e57de6e Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Thu, 19 Dec 2024 15:32:00 +0100 Subject: [PATCH 09/16] bump cmake min version as we use install on a target from a subdirectory --- cpp/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 8becf01e39..bee7f41eaa 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.11) +cmake_minimum_required(VERSION 3.13) project(memilio VERSION 1.0.0) From ad52f9b3f8f095190420eb885791e7df366224db Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 3 Jan 2025 14:33:08 +0100 Subject: [PATCH 10/16] Update readme and comments according to review --- cpp/README.md | 6 +++--- cpp/memilio/utils/custom_index_array.h | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/README.md b/cpp/README.md index 11ba6e6b49..d6519b0c71 100644 --- a/cpp/README.md +++ b/cpp/README.md @@ -27,10 +27,10 @@ The following table lists the dependencies that are used. Most of them are requi | Library | Version | Required | Bundled | Notes | |---------|----------|----------|-----------------------|-------| -| spdlog | 1.11.0 | Yes | Yes (git repo) | https://github.com/gabime/spdlog | -| Eigen | 3.3.9 | Yes | Yes (git repo) | http://gitlab.com/libeigen/eigen | +| spdlog | 1.15.0 | Yes | Yes (git repo) | https://github.com/gabime/spdlog | +| Eigen | 3.4.0 | Yes | Yes (git repo) | http://gitlab.com/libeigen/eigen | | Boost | 1.84.0 | Yes | Yes (git repo) | https://github.com/boostorg/boost | -| JsonCpp | 1.9.5 | No | Yes (git repo) | https://github.com/open-source-parsers/jsoncpp | +| JsonCpp | 1.9.6 | No | Yes (git repo) | https://github.com/open-source-parsers/jsoncpp | | HDF5 | 1.12.0 | No | No | https://www.hdfgroup.org/, package libhdf5-dev on apt (Ubuntu) | | GoogleTest | 1.10 | For Tests only | Yes (git repo) | https://github.com/google/googletest | diff --git a/cpp/memilio/utils/custom_index_array.h b/cpp/memilio/utils/custom_index_array.h index 4c141b68c0..aa7b4a0a90 100644 --- a/cpp/memilio/utils/custom_index_array.h +++ b/cpp/memilio/utils/custom_index_array.h @@ -236,7 +236,7 @@ class CustomIndexArray /** * @brief Resize all dimensions. - ' Note that when increasing the overall size, new values may be uninitialized. + * Note that when increasing the overall size, new values may be uninitialized. * @param new_dims New dimensions. */ void resize(Index new_dims) @@ -248,7 +248,7 @@ class CustomIndexArray /** * @brief Resize all dimensions, destroying all values. - * This Version of resize should only be used when the CustomIndexArray contains non-movable and non-copyable + * This version of resize should only be used when the CustomIndexArray contains non-movable and non-copyable * values, like atomics. New entries are all default initialized. * @param new_dims New dimensions. */ From f527029fdcb3ec2d0715c92036ecb62725010a21 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 3 Jan 2025 14:34:30 +0100 Subject: [PATCH 11/16] remove unused cmake code (minimum version exceeds target of workaround) --- cpp/benchmarks/CMakeLists.txt | 35 +++++++++-------------------------- 1 file changed, 9 insertions(+), 26 deletions(-) diff --git a/cpp/benchmarks/CMakeLists.txt b/cpp/benchmarks/CMakeLists.txt index b6920b8e23..731fe8d566 100755 --- a/cpp/benchmarks/CMakeLists.txt +++ b/cpp/benchmarks/CMakeLists.txt @@ -5,34 +5,17 @@ set(BENCHMARK_ENABLE_INSTALL OFF CACHE BOOL "Don't install benchmark" FORCE) set(BENCHMARK_DOWNLOAD_DEPENDENCIES OFF CACHE BOOL "Don't download dependencies" FORCE) set(BENCHMARK_ENABLE_GTEST_TESTS OFF CACHE BOOL "Disable Google Test in benchmark" FORCE) -if(CMAKE_VERSION VERSION_LESS 3.11) - set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1") - - include(DownloadProject) - download_project(PROJ benchmark - GIT_REPOSITORY https://github.com/google/benchmark.git - GIT_TAG v1.6.1 - UPDATE_DISCONNECTED 1 - QUIET - ) - - # CMake warning suppression will not be needed in version 1.9 +include(FetchContent) +FetchContent_Declare(benchmark + GIT_REPOSITORY https://github.com/google/benchmark.git + GIT_TAG v1.6.1) +FetchContent_GetProperties(benchmark) + +if(NOT benchmark_POPULATED) + FetchContent_Populate(benchmark) set(CMAKE_SUPPRESS_DEVELOPER_WARNINGS 1 CACHE BOOL "") - add_subdirectory(${benchmark_SOURCE_DIR} ${benchmark_SOURCE_DIR} EXCLUDE_FROM_ALL) + add_subdirectory(${benchmark_SOURCE_DIR} ${benchmark_BINARY_DIR} EXCLUDE_FROM_ALL) unset(CMAKE_SUPPRESS_DEVELOPER_WARNINGS) -else() - include(FetchContent) - FetchContent_Declare(benchmark - GIT_REPOSITORY https://github.com/google/benchmark.git - GIT_TAG v1.6.1) - FetchContent_GetProperties(benchmark) - - if(NOT benchmark_POPULATED) - FetchContent_Populate(benchmark) - set(CMAKE_SUPPRESS_DEVELOPER_WARNINGS 1 CACHE BOOL "") - add_subdirectory(${benchmark_SOURCE_DIR} ${benchmark_BINARY_DIR} EXCLUDE_FROM_ALL) - unset(CMAKE_SUPPRESS_DEVELOPER_WARNINGS) - endif() endif() set_target_properties(benchmark PROPERTIES FOLDER "Extern") From a9a870100385b4d2cfcf23072fd6662dcdf88ed4 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 3 Jan 2025 14:35:00 +0100 Subject: [PATCH 12/16] improve format_as return types --- cpp/memilio/utils/logging.h | 2 +- cpp/memilio/utils/uncertain_value.h | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/memilio/utils/logging.h b/cpp/memilio/utils/logging.h index ac3f79376b..f6ef1dc69b 100644 --- a/cpp/memilio/utils/logging.h +++ b/cpp/memilio/utils/logging.h @@ -138,7 +138,7 @@ namespace internal * If derivative information is needed as well, use `ad::derivative(...)` or define a `fmt::formatter<...>`. */ template -FP format_as(const active_type& ad_type) +const FP& format_as(const active_type& ad_type) { // Note: the format_as function needs to be in the same namespace as the value it takes return value(ad_type); diff --git a/cpp/memilio/utils/uncertain_value.h b/cpp/memilio/utils/uncertain_value.h index 27d660c72e..90336d29ab 100644 --- a/cpp/memilio/utils/uncertain_value.h +++ b/cpp/memilio/utils/uncertain_value.h @@ -231,8 +231,9 @@ class UncertainValue * @brief Format UncertainValues using their value for logging with spdlog. */ template -double format_as(const UncertainValue& uv) +const FP& format_as(const UncertainValue& uv) { + // uses UncertainValue::operator const FP&() const return uv; } From 46d0ef9ca9cff999280700adf02f5b1402344ffb Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 3 Jan 2025 15:24:22 +0100 Subject: [PATCH 13/16] remove compiler diagnostics around eigen, replace mio::Vector<> by Eigen::VectorX<> --- cpp/examples/glct_secir.cpp | 16 +- cpp/examples/lct_secir.cpp | 8 +- cpp/memilio/compartments/compartmentalmodel.h | 28 ++-- cpp/memilio/compartments/flow_model.h | 29 ++-- cpp/memilio/compartments/flow_simulation.h | 4 +- cpp/memilio/compartments/simulation.h | 2 +- .../epidemiology/lct_infection_state.h | 4 +- cpp/memilio/epidemiology/lct_populations.h | 2 +- cpp/memilio/epidemiology/populations.h | 2 +- cpp/memilio/math/adapt_rk.h | 2 +- cpp/memilio/math/eigen.h | 26 +-- cpp/memilio/math/euler.h | 4 +- cpp/memilio/math/integrator.h | 9 +- cpp/memilio/math/stepper_wrapper.h | 25 +-- cpp/models/glct_secir/model.h | 28 ++-- cpp/models/glct_secir/parameters.h | 75 +++++---- cpp/models/ide_secir/model.cpp | 4 +- cpp/models/lct_secir/initializer_flows.h | 17 +- cpp/models/lct_secir/model.h | 21 ++- cpp/models/lct_secir/parameters.h | 24 +-- cpp/models/ode_seair/model.h | 4 +- cpp/models/ode_secir/model.h | 19 ++- cpp/models/ode_secirts/model.h | 14 +- cpp/models/ode_secirvvs/model.h | 14 +- cpp/models/ode_seir/model.h | 4 +- cpp/models/ode_sir/model.h | 4 +- cpp/models/sde_seirvv/model.h | 155 ++++++++++-------- cpp/models/sde_sir/model.h | 5 +- cpp/models/sde_sir/simulation.h | 4 +- cpp/models/sde_sirs/model.h | 4 +- cpp/models/sde_sirs/simulation.h | 4 +- cpp/tests/test_glct_secir.cpp | 32 ++-- cpp/tests/test_ide_parameters_io.cpp | 2 +- cpp/tests/test_ide_secir_ageres.cpp | 6 +- cpp/tests/test_lct_initializer_flows.cpp | 37 +++-- cpp/tests/test_lct_parameters_io.cpp | 8 +- cpp/tests/test_lct_secir.cpp | 34 ++-- 37 files changed, 347 insertions(+), 333 deletions(-) diff --git a/cpp/examples/glct_secir.cpp b/cpp/examples/glct_secir.cpp index ba1718f15c..59b08df184 100644 --- a/cpp/examples/glct_secir.cpp +++ b/cpp/examples/glct_secir.cpp @@ -133,8 +133,8 @@ int main() // InfectedNoSymptomsToInfectedSymptoms and one for InfectedNoSymptomsToRecovered. // The strains have a length of NumInfectedNoSymptoms/2. each. // The transition probability is included in the StartingProbability vector. - mio::Vector StartingProbabilitiesInfectedNoSymptoms = - mio::Vector::Zero(LctState::get_num_subcompartments()); + Eigen::VectorX StartingProbabilitiesInfectedNoSymptoms = + Eigen::VectorX::Zero(LctState::get_num_subcompartments()); StartingProbabilitiesInfectedNoSymptoms[0] = 1 - recoveredPerInfectedNoSymptoms; StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = recoveredPerInfectedNoSymptoms; @@ -152,8 +152,8 @@ int main() timeInfectedNoSymptoms); // Do the same for all compartments. // InfectedSymptoms. - mio::Vector StartingProbabilitiesInfectedSymptoms = - mio::Vector::Zero(LctState::get_num_subcompartments()); + Eigen::VectorX StartingProbabilitiesInfectedSymptoms = + Eigen::VectorX::Zero(LctState::get_num_subcompartments()); StartingProbabilitiesInfectedSymptoms[0] = severePerInfectedSymptoms; StartingProbabilitiesInfectedSymptoms[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = 1 - severePerInfectedSymptoms; @@ -165,8 +165,8 @@ int main() mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSymptoms); // InfectedSevere. - mio::Vector StartingProbabilitiesInfectedSevere = - mio::Vector::Zero(LctState::get_num_subcompartments()); + Eigen::VectorX StartingProbabilitiesInfectedSevere = + Eigen::VectorX::Zero(LctState::get_num_subcompartments()); StartingProbabilitiesInfectedSevere[0] = criticalPerSevere; StartingProbabilitiesInfectedSevere[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = 1 - criticalPerSevere; @@ -178,8 +178,8 @@ int main() mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSevere); // InfectedCritical. - mio::Vector StartingProbabilitiesInfectedCritical = - mio::Vector::Zero(LctState::get_num_subcompartments()); + Eigen::VectorX StartingProbabilitiesInfectedCritical = + Eigen::VectorX::Zero(LctState::get_num_subcompartments()); StartingProbabilitiesInfectedCritical[0] = deathsPerCritical; StartingProbabilitiesInfectedCritical[(Eigen::Index)( LctState::get_num_subcompartments() / 2.)] = 1 - deathsPerCritical; diff --git a/cpp/examples/lct_secir.cpp b/cpp/examples/lct_secir.cpp index 6678633171..9c9d501f84 100644 --- a/cpp/examples/lct_secir.cpp +++ b/cpp/examples/lct_secir.cpp @@ -75,10 +75,10 @@ int main() if (use_initializer_flows) { // Example how to use the class Initializer for the definition of an initial vector for the LCT model. - ScalarType dt = 0.001; - mio::Vector total_population = mio::Vector::Constant(1, 1000000.); - mio::Vector deaths = mio::Vector::Constant(1, 10.); - mio::Vector total_confirmed_cases = mio::Vector::Constant(1, 16000.); + ScalarType dt = 0.001; + Eigen::VectorX total_population = Eigen::VectorX::Constant(1, 1000000.); + Eigen::VectorX deaths = Eigen::VectorX::Constant(1, 10.); + Eigen::VectorX total_confirmed_cases = Eigen::VectorX::Constant(1, 16000.); // Create TimeSeries with num_transitions elements. int num_transitions = (int)mio::lsecir::InfectionTransition::Count; diff --git a/cpp/memilio/compartments/compartmentalmodel.h b/cpp/memilio/compartments/compartmentalmodel.h index 53e47e89c4..6b0aed725f 100644 --- a/cpp/memilio/compartments/compartmentalmodel.h +++ b/cpp/memilio/compartments/compartmentalmodel.h @@ -88,15 +88,15 @@ struct CompartmentalModel { { } - CompartmentalModel(const CompartmentalModel&) = default; - CompartmentalModel(CompartmentalModel&&) = default; + CompartmentalModel(const CompartmentalModel&) = default; + CompartmentalModel(CompartmentalModel&&) = default; CompartmentalModel& operator=(const CompartmentalModel&) = default; - CompartmentalModel& operator=(CompartmentalModel&&) = default; - virtual ~CompartmentalModel() = default; + CompartmentalModel& operator=(CompartmentalModel&&) = default; + virtual ~CompartmentalModel() = default; // REMARK: Not pure virtual for easier java/python bindings. - virtual void get_derivatives(Eigen::Ref>, Eigen::Ref> /*y*/, FP /*t*/, - Eigen::Ref> /*dydt*/) const {}; + virtual void get_derivatives(Eigen::Ref>, Eigen::Ref> /*y*/, + FP /*t*/, Eigen::Ref> /*dydt*/) const {}; /** * @brief This function evaluates the right-hand-side f of the ODE dydt = f(y, t). @@ -118,8 +118,8 @@ struct CompartmentalModel { * @param[in] t The current time. * @param[out] dydt A reference to the calculated output. */ - void eval_right_hand_side(Eigen::Ref> pop, Eigen::Ref> y, FP t, - Eigen::Ref> dydt) const + void eval_right_hand_side(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> dydt) const { dydt.setZero(); this->get_derivatives(pop, y, t, dydt); @@ -130,7 +130,7 @@ struct CompartmentalModel { * See eval_right_hand_side for more detail. * @return Current value of model populations as a flat vector. */ - Vector get_initial_values() const + Eigen::VectorX get_initial_values() const { return populations.get_compartments(); } @@ -183,10 +183,9 @@ struct CompartmentalModel { * @tparam M a type that has a eval_right_hand_side member function, e.g. a compartment model type. */ template -using eval_right_hand_side_expr_t = - decltype(std::declval().eval_right_hand_side(std::declval>>(), - std::declval>>(), - std::declval(), std::declval>>())); +using eval_right_hand_side_expr_t = decltype(std::declval().eval_right_hand_side( + std::declval>>(), std::declval>>(), + std::declval(), std::declval>>())); /** * Detect the get_initial_values member function of a compartment model. @@ -197,7 +196,8 @@ using eval_right_hand_side_expr_t = * @tparam M a type that has a get_initial_values member function, e.g. a compartment model type. */ template -using get_initial_values_expr_t = decltype(std::declval&>() = std::declval().get_initial_values()); +using get_initial_values_expr_t = + decltype(std::declval&>() = std::declval().get_initial_values()); /** * Template meta function to check if a type is a valid compartment model. diff --git a/cpp/memilio/compartments/flow_model.h b/cpp/memilio/compartments/flow_model.h index 41581633cd..9c1823359b 100644 --- a/cpp/memilio/compartments/flow_model.h +++ b/cpp/memilio/compartments/flow_model.h @@ -97,8 +97,8 @@ class FlowModel : public CompartmentalModel // Note: use get_flat_flow_index when accessing flows // Note: by convention, we compute incoming flows, thus entries in flows must be non-negative - virtual void get_flows(Eigen::Ref> /*pop*/, Eigen::Ref> /*y*/, FP /*t*/, - Eigen::Ref> /*flows*/) const = 0; + virtual void get_flows(Eigen::Ref> /*pop*/, Eigen::Ref> /*y*/, + FP /*t*/, Eigen::Ref> /*flows*/) const = 0; /** * @brief Compute the right-hand-side of the ODE dydt = f(y, t) from flow values. @@ -108,7 +108,7 @@ class FlowModel : public CompartmentalModel * @param[in] flows The current flow values (as calculated by get_flows) as a flat array. * @param[out] dydt A reference to the calculated output. */ - void get_derivatives(Eigen::Ref> flows, Eigen::Ref> dydt) const + void get_derivatives(Eigen::Ref> flows, Eigen::Ref> dydt) const { // set dydt to 0, then iteratively add all flow contributions dydt.setZero(); @@ -134,8 +134,8 @@ class FlowModel : public CompartmentalModel * @param[in] t The current time. * @param[out] dydt A reference to the calculated output. */ - void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, FP t, - Eigen::Ref> dydt) const override final + void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> dydt) const override final { m_flow_values.setZero(); get_flows(pop, y, t, m_flow_values); @@ -147,9 +147,9 @@ class FlowModel : public CompartmentalModel * This can be used as initial conditions in an ODE solver. By default, this is a zero vector. * @return The initial flows. */ - Vector get_initial_flows() const + Eigen::VectorX get_initial_flows() const { - return Vector::Zero((this->populations.numel() / static_cast(Comp::Count)) * Flows::size()); + return Eigen::VectorX::Zero((this->populations.numel() / static_cast(Comp::Count)) * Flows::size()); } /** @@ -204,7 +204,7 @@ class FlowModel : public CompartmentalModel } private: - mutable Vector m_flow_values; ///< Cache to avoid allocation in get_derivatives (using get_flows). + mutable Eigen::VectorX m_flow_values; ///< Cache to avoid allocation in get_derivatives (using get_flows). /** * @brief Compute the derivatives of the compartments. @@ -216,7 +216,7 @@ class FlowModel : public CompartmentalModel * @tparam I The index of a flow in FlowChart. */ template - inline void get_rhs_impl(Eigen::Ref> flows, Eigen::Ref> rhs, + inline void get_rhs_impl(Eigen::Ref> flows, Eigen::Ref> rhs, const FlowIndex& index) const { using Flow = type_at_index_t; @@ -244,16 +244,17 @@ class FlowModel : public CompartmentalModel */ template using get_derivatives_expr_t = decltype(std::declval().get_derivatives( - std::declval>>(), std::declval>>())); + std::declval>>(), std::declval>>())); template using get_flows_expr_t = - decltype(std::declval().get_flows(std::declval>>(), - std::declval>>(), std::declval(), - std::declval>>())); + decltype(std::declval().get_flows(std::declval>>(), + std::declval>>(), + std::declval(), std::declval>>())); template -using get_initial_flows_expr_t = decltype(std::declval>() = std::declval().get_initial_flows()); +using get_initial_flows_expr_t = + decltype(std::declval>() = std::declval().get_initial_flows()); /** @} */ /** diff --git a/cpp/memilio/compartments/flow_simulation.h b/cpp/memilio/compartments/flow_simulation.h index 4c44395a55..04e9113bd9 100644 --- a/cpp/memilio/compartments/flow_simulation.h +++ b/cpp/memilio/compartments/flow_simulation.h @@ -58,7 +58,7 @@ class FlowSimulation : public Simulation * tmax must be greater than get_result().get_last_time_point(). * @param[in] tmax Next stopping time of the simulation. */ - Eigen::Ref> advance(FP tmax) + Eigen::Ref> advance(FP tmax) { // the derivfunktion (i.e. the lambda passed to m_integrator.advance below) requires that there are at least // as many entries in m_flow_result as in Base::m_result @@ -132,7 +132,7 @@ class FlowSimulation : public Simulation } } - Vector m_pop; ///< pre-allocated temporary, used in right_hand_side() + Eigen::VectorX m_pop; ///< pre-allocated temporary, used in right_hand_side() private: mio::TimeSeries m_flow_result; ///< flow result of the simulation diff --git a/cpp/memilio/compartments/simulation.h b/cpp/memilio/compartments/simulation.h index 8b3307f87f..7594070672 100644 --- a/cpp/memilio/compartments/simulation.h +++ b/cpp/memilio/compartments/simulation.h @@ -95,7 +95,7 @@ class Simulation * tmax must be greater than get_result().get_last_time_point() * @param tmax next stopping point of simulation */ - Eigen::Ref> advance(FP tmax) + Eigen::Ref> advance(FP tmax) { return m_integrator.advance( [this](auto&& y, auto&& t, auto&& dydt) { diff --git a/cpp/memilio/epidemiology/lct_infection_state.h b/cpp/memilio/epidemiology/lct_infection_state.h index bf80a80079..c5a7344b11 100644 --- a/cpp/memilio/epidemiology/lct_infection_state.h +++ b/cpp/memilio/epidemiology/lct_infection_state.h @@ -90,11 +90,11 @@ class LctInfectionState * The size of the vector has to match the LctInfectionState. * @return Vector with accumulated values for the InfectionStates. */ - static Vector calculate_compartments(const Vector& subcompartments) + static Eigen::VectorX calculate_compartments(const Eigen::VectorX& subcompartments) { assert(subcompartments.rows() == Count); - Vector compartments((Eigen::Index)InfectionState::Count); + Eigen::VectorX compartments((Eigen::Index)InfectionState::Count); // Use segment of the vector subcompartments of each InfectionState and sum up the values of subcompartments. compartments[(Eigen::Index)InfectionState::Susceptible] = subcompartments[0]; compartments[(Eigen::Index)InfectionState::Exposed] = diff --git a/cpp/memilio/epidemiology/lct_populations.h b/cpp/memilio/epidemiology/lct_populations.h index 88e5d835d4..1f0ffe3aa0 100644 --- a/cpp/memilio/epidemiology/lct_populations.h +++ b/cpp/memilio/epidemiology/lct_populations.h @@ -123,7 +123,7 @@ class LctPopulations * This can be used as initial conditions for the ODE solver. * @return Eigen::VectorXd of populations. */ - inline Vector get_compartments() const + inline Eigen::VectorX get_compartments() const { return m_y.array().template cast(); } diff --git a/cpp/memilio/epidemiology/populations.h b/cpp/memilio/epidemiology/populations.h index 7b7852935c..e9739d2b10 100644 --- a/cpp/memilio/epidemiology/populations.h +++ b/cpp/memilio/epidemiology/populations.h @@ -79,7 +79,7 @@ class Populations : public CustomIndexArray, Categories...> * as initial conditions for the ODE solver * @return Eigen::VectorXd of populations */ - inline Vector get_compartments() const + inline Eigen::VectorX get_compartments() const { return this->array().template cast(); } diff --git a/cpp/memilio/math/adapt_rk.h b/cpp/memilio/math/adapt_rk.h index ab3c700d65..3cbb8563b9 100644 --- a/cpp/memilio/math/adapt_rk.h +++ b/cpp/memilio/math/adapt_rk.h @@ -292,7 +292,7 @@ class RKIntegratorCore : public IntegratorCore TableauFinal m_tab_final; FP m_abs_tol, m_rel_tol; mutable Eigen::Matrix m_kt_values; - mutable Vector m_yt_eval; + mutable Eigen::VectorX m_yt_eval; private: mutable Eigen::Array m_eps, diff --git a/cpp/memilio/math/eigen.h b/cpp/memilio/math/eigen.h index effd111876..8f3ea22223 100644 --- a/cpp/memilio/math/eigen.h +++ b/cpp/memilio/math/eigen.h @@ -20,31 +20,9 @@ #ifndef MIO_UTILS_EIGEN_H #define MIO_UTILS_EIGEN_H -#include "memilio/config.h" -#include "memilio/utils/compiler_diagnostics.h" - -/* this file wraps includes from eigen3 library to disable warnings. */ - -//C4996: some std functions that have been deprecated in c++17; maybe fixed in new eigen versions? -MSVC_WARNING_DISABLE_PUSH(4996) - -GCC_CLANG_DIAGNOSTIC(push) -GCC_CLANG_DIAGNOSTIC(ignored "-Wint-in-bool-context") -GCC_CLANG_DIAGNOSTIC(ignored "-Wshadow") +// this file wraps includes from eigen3 library to disable warnings +// Eigen is treated as a system library through cmake since #1168, so ALL warnings within the library are disabled #include -GCC_CLANG_DIAGNOSTIC(pop) - -MSVC_WARNING_POP() - -namespace mio -{ - -/// A vector of type FP from the Eigen library -template -using Vector = Eigen::Matrix; - -} // namespace mio - #endif // MIO_UTILS_EIGEN_H diff --git a/cpp/memilio/math/euler.h b/cpp/memilio/math/euler.h index 7c3c569141..45da43b6da 100644 --- a/cpp/memilio/math/euler.h +++ b/cpp/memilio/math/euler.h @@ -47,8 +47,8 @@ class EulerIntegratorCore : public IntegratorCore * @param[in,out] dt current time step h=dt * @param[out] ytp1 approximated value y(t+1) */ - bool step(const DerivFunction& f, Eigen::Ref> yt, FP& t, FP& dt, - Eigen::Ref> ytp1) const override + bool step(const DerivFunction& f, Eigen::Ref> yt, FP& t, FP& dt, + Eigen::Ref> ytp1) const override { // we are misusing the next step y as temporary space to store the derivative f(yt, t, ytp1); diff --git a/cpp/memilio/math/integrator.h b/cpp/memilio/math/integrator.h index 94214f529f..02d6d10cf5 100644 --- a/cpp/memilio/math/integrator.h +++ b/cpp/memilio/math/integrator.h @@ -34,7 +34,8 @@ namespace mio * Function template to be integrated. */ template -using DerivFunction = std::function> y, FP t, Eigen::Ref> dydt)>; +using DerivFunction = + std::function> y, FP t, Eigen::Ref> dydt)>; template class IntegratorCore @@ -79,8 +80,8 @@ class IntegratorCore * @return Always true for nonadaptive methods. * (If adaptive, returns whether the adaptive step sizing was successful.) */ - virtual bool step(const DerivFunction& f, Eigen::Ref> yt, FP& t, FP& dt, - Eigen::Ref> ytp1) const = 0; + virtual bool step(const DerivFunction& f, Eigen::Ref> yt, FP& t, FP& dt, + Eigen::Ref> ytp1) const = 0; /** * @brief Access lower bound to the step size dt. @@ -146,7 +147,7 @@ class OdeIntegrator * @return A reference to the last value in the results time series. */ - Eigen::Ref> advance(const DerivFunction& f, const FP tmax, FP& dt, TimeSeries& results) + Eigen::Ref> advance(const DerivFunction& f, const FP tmax, FP& dt, TimeSeries& results) { // hint at std functions for ADL using std::fabs; diff --git a/cpp/memilio/math/stepper_wrapper.h b/cpp/memilio/math/stepper_wrapper.h index 41a6d818f9..85b1ae40e9 100644 --- a/cpp/memilio/math/stepper_wrapper.h +++ b/cpp/memilio/math/stepper_wrapper.h @@ -43,8 +43,8 @@ template { using Stepper = boost::numeric::odeint::controlled_runge_kutta< - ControlledStepper, FP, Vector, FP, boost::numeric::odeint::vector_space_algebra, - typename boost::numeric::odeint::operations_dispatcher>::operations_type, + ControlledStepper, FP, Eigen::VectorX, FP, boost::numeric::odeint::vector_space_algebra, + typename boost::numeric::odeint::operations_dispatcher>::operations_type, boost::numeric::odeint::initially_resizer>>; static constexpr bool is_fsal_stepper = std::is_same_v; @@ -77,8 +77,8 @@ class ControlledStepperWrapper : public mio::IntegratorCore * @param[in,out] dt Current time step size h=dt. Overwritten by an estimated optimal step size for the next step. * @param[out] ytp1 The approximated value of y(t'). */ - bool step(const mio::DerivFunction& f, Eigen::Ref const> yt, FP& t, FP& dt, - Eigen::Ref> ytp1) const override + bool step(const mio::DerivFunction& f, Eigen::Ref const> yt, FP& t, FP& dt, + Eigen::Ref> ytp1) const override { using boost::numeric::odeint::fail; using std::max; @@ -109,7 +109,7 @@ class ControlledStepperWrapper : public mio::IntegratorCore if constexpr (!is_fsal_stepper) { // prevent compile time errors with fsal steppers step_result = m_stepper.try_step( // reorder arguments of the DerivFunction f for the stepper - [&](const Vector& x, Vector& dxds, FP s) { + [&](const Eigen::VectorX& x, Eigen::VectorX& dxds, FP s) { f(x, s, dxds); }, m_yt, t, m_ytp1, dt); @@ -170,7 +170,7 @@ class ControlledStepperWrapper : public mio::IntegratorCore } FP m_abs_tol, m_rel_tol; ///< Absolute and relative tolerances for integration. - mutable Vector m_ytp1, m_yt; ///< Temporary storage to avoid allocations in step function. + mutable Eigen::VectorX m_ytp1, m_yt; ///< Temporary storage to avoid allocations in step function. mutable Stepper m_stepper; ///< A stepper instance used for integration. }; @@ -184,9 +184,10 @@ template { public: - using Stepper = ExplicitStepper, FP, Vector, FP, boost::numeric::odeint::vector_space_algebra, - typename boost::numeric::odeint::operations_dispatcher>::operations_type, - boost::numeric::odeint::initially_resizer>; + using Stepper = + ExplicitStepper, FP, Eigen::VectorX, FP, boost::numeric::odeint::vector_space_algebra, + typename boost::numeric::odeint::operations_dispatcher>::operations_type, + boost::numeric::odeint::initially_resizer>; /** * @brief Set up the integrator. @@ -204,15 +205,15 @@ class ExplicitStepperWrapper : public mio::IntegratorCore * @param[in] dt Current time step size h=dt. * @param[out] ytp1 The approximated value of y(t+dt). */ - bool step(const mio::DerivFunction& f, Eigen::Ref const> yt, FP& t, FP& dt, - Eigen::Ref> ytp1) const override + bool step(const mio::DerivFunction& f, Eigen::Ref const> yt, FP& t, FP& dt, + Eigen::Ref> ytp1) const override { // copy the values from y(t) to ytp1, since we use the scheme do_step(sys, inout, t, dt) with // sys=f, inout=y(t) for in-place computation - also, this form is shared by several steppers in boost ytp1 = yt; m_stepper.do_step( // reorder arguments of the DerivFunction f for the stepper - [&](const Vector& x, Vector& dxds, FP s) { + [&](const Eigen::VectorX& x, Eigen::VectorX& dxds, FP s) { f(x, s, dxds); }, ytp1, t, dt); diff --git a/cpp/models/glct_secir/model.h b/cpp/models/glct_secir/model.h index 6501fe5d53..a90bea43a7 100644 --- a/cpp/models/glct_secir/model.h +++ b/cpp/models/glct_secir/model.h @@ -134,8 +134,9 @@ class Model * @param[in] t The current time. * @param[out] dydt A reference to the calculated output. */ - void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, ScalarType t, - Eigen::Ref> dydt) const override + void get_derivatives(Eigen::Ref> pop, + Eigen::Ref> y, ScalarType t, + Eigen::Ref> dydt) const override { dydt.setZero(); @@ -176,7 +177,7 @@ class Model dydt.segment(LctState::template get_first_index(), LctState::template get_num_subcompartments()) = -(params.template get() * - Vector::Ones(LctState::template get_num_subcompartments())) + Eigen::VectorX::Ones(LctState::template get_num_subcompartments())) .transpose() * y.segment(LctState::template get_first_index(), LctState::template get_num_subcompartments()) * @@ -202,7 +203,7 @@ class Model // Add flow directly to Recovered compartment. dydt[LctState::template get_first_index()] += -(params.template get() * - Vector::Ones(dimensionInfectedNoSymptomsToRecovered)) + Eigen::VectorX::Ones(dimensionInfectedNoSymptomsToRecovered)) .transpose() * y.segment(LctState::template get_first_index() + dimensionInfectedNoSymptomsToInfectedSymptoms, @@ -214,7 +215,7 @@ class Model LctState::template get_num_subcompartments()) = -params.template get() * (params.template get() * - Vector::Ones(dimensionInfectedNoSymptomsToInfectedSymptoms)) + Eigen::VectorX::Ones(dimensionInfectedNoSymptomsToInfectedSymptoms)) .transpose() * y.segment(LctState::template get_first_index(), dimensionInfectedNoSymptomsToInfectedSymptoms); @@ -239,7 +240,7 @@ class Model // Add flow directly to Recovered compartment. dydt[LctState::template get_first_index()] += -(params.template get() * - Vector::Ones(dimensionInfectedSymptomsToRecovered)) + Eigen::VectorX::Ones(dimensionInfectedSymptomsToRecovered)) .transpose() * y.segment(LctState::template get_first_index() + dimensionInfectedSymptomsToInfectedSevere, @@ -251,7 +252,7 @@ class Model LctState::template get_num_subcompartments()) = -params.template get() * (params.template get() * - Vector::Ones(dimensionInfectedSymptomsToInfectedSevere)) + Eigen::VectorX::Ones(dimensionInfectedSymptomsToInfectedSevere)) .transpose() * y.segment(LctState::template get_first_index(), dimensionInfectedSymptomsToInfectedSevere); @@ -276,7 +277,7 @@ class Model // Add flow directly to Recovered compartment. dydt[LctState::template get_first_index()] += -(params.template get() * - Vector::Ones(dimensionInfectedSevereToRecovered)) + Eigen::VectorX::Ones(dimensionInfectedSevereToRecovered)) .transpose() * y.segment(LctState::template get_first_index() + dimensionInfectedSevereToInfectedCritical, @@ -288,7 +289,7 @@ class Model LctState::template get_num_subcompartments()) = -params.template get() * (params.template get() * - Vector::Ones(dimensionInfectedSevereToInfectedCritical)) + Eigen::VectorX::Ones(dimensionInfectedSevereToInfectedCritical)) .transpose() * y.segment(LctState::template get_first_index(), dimensionInfectedSevereToInfectedCritical); @@ -312,7 +313,7 @@ class Model // Add flow directly to Recovered compartment. dydt[LctState::template get_first_index()] += -(params.template get() * - Vector::Ones(dimensionInfectedCriticalToRecovered)) + Eigen::VectorX::Ones(dimensionInfectedCriticalToRecovered)) .transpose() * y.segment(LctState::template get_first_index() + dimensionInfectedCriticalToDead, @@ -321,7 +322,7 @@ class Model // --- Dead. --- dydt[LctState::template get_first_index()] = -(params.template get() * - Vector::Ones(dimensionInfectedCriticalToDead)) + Eigen::VectorX::Ones(dimensionInfectedCriticalToDead)) .transpose() * y.segment(LctState::template get_first_index(), dimensionInfectedCriticalToDead); @@ -345,11 +346,12 @@ class Model if (!(LctState::Count == subcompartments_ts.get_num_elements())) { log_error("Result does not match InfectionStates of the model."); // Return a TimeSeries with values -1. - Vector error_output = Vector::Constant((Eigen::Index)InfectionState::Count, -1); + Eigen::VectorX error_output = + Eigen::VectorX::Constant((Eigen::Index)InfectionState::Count, -1); compartments_ts.add_time_point(-1, error_output); return compartments_ts; } - Vector compartments((Eigen::Index)InfectionState::Count); + Eigen::VectorX compartments((Eigen::Index)InfectionState::Count); for (Eigen::Index i = 0; i < subcompartments_ts.get_num_time_points(); ++i) { compartments_ts.add_time_point(subcompartments_ts.get_time(i), LctState::calculate_compartments(subcompartments_ts[i])); diff --git a/cpp/models/glct_secir/parameters.h b/cpp/models/glct_secir/parameters.h index e87b72b7ea..c6a9f45580 100644 --- a/cpp/models/glct_secir/parameters.h +++ b/cpp/models/glct_secir/parameters.h @@ -39,15 +39,15 @@ namespace glsecir /// @brief Vector with the probability to start in any of the subcompartments of the Exposed compartment. struct StartingProbabilitiesExposed { - using Type = Vector; + using Type = Eigen::VectorX; /** * @brief Default parameters can be used to get an Erlang distributed stay time in the Exposed compartment. * @param[in] numExposed Number of subcompartments of the Exposed compartment. */ static Type get_default(size_t numExposed) { - Vector def = Vector::Zero(numExposed); - def[0] = 1.; + Eigen::VectorX def = Eigen::VectorX::Zero(numExposed); + def[0] = 1.; return def; } static std::string name() @@ -67,7 +67,7 @@ struct TransitionMatrixExposedToInfectedNoSymptoms { static Type get_default(size_t numExposed, ScalarType timeExposed = 1.) { Eigen::MatrixXd def = - Vector::Constant(numExposed, -(ScalarType)numExposed / timeExposed).asDiagonal(); + Eigen::VectorX::Constant(numExposed, -(ScalarType)numExposed / timeExposed).asDiagonal(); def.diagonal(1).setConstant((ScalarType)numExposed / timeExposed); return def; } @@ -79,15 +79,15 @@ struct TransitionMatrixExposedToInfectedNoSymptoms { /// @brief Vector with the probability to start in any of the subcompartments of the InfectedNoSymptoms compartment. struct StartingProbabilitiesInfectedNoSymptoms { - using Type = Vector; + using Type = Eigen::VectorX; /** * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedNoSymptoms compartment. * @param[in] numInfectedNoSymptoms Number of subcompartments of the InfectedNoSymptoms compartment. */ static Type get_default(size_t numInfectedNoSymptoms) { - Vector def = Vector::Zero(numInfectedNoSymptoms); - def[0] = 1.; + Eigen::VectorX def = Eigen::VectorX::Zero(numInfectedNoSymptoms); + def[0] = 1.; return def; } static std::string name() @@ -110,7 +110,8 @@ struct TransitionMatrixInfectedNoSymptomsToInfectedSymptoms { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -134,7 +135,8 @@ struct TransitionMatrixInfectedNoSymptomsToRecovered { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -146,15 +148,15 @@ struct TransitionMatrixInfectedNoSymptomsToRecovered { /// @brief Vector with the probability to start in any of the subcompartments of the InfectedSymptoms compartment. struct StartingProbabilitiesInfectedSymptoms { - using Type = Vector; + using Type = Eigen::VectorX; /** * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedSymptoms compartment. * @param[in] numInfectedSymptoms Number of subcompartments of the InfectedSymptoms compartment. */ static Type get_default(size_t numInfectedSymptoms) { - Vector def = Vector::Zero(numInfectedSymptoms); - def[0] = 1.; + Eigen::VectorX def = Eigen::VectorX::Zero(numInfectedSymptoms); + def[0] = 1.; return def; } static std::string name() @@ -177,7 +179,8 @@ struct TransitionMatrixInfectedSymptomsToInfectedSevere { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -201,7 +204,8 @@ struct TransitionMatrixInfectedSymptomsToRecovered { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -213,15 +217,15 @@ struct TransitionMatrixInfectedSymptomsToRecovered { /// @brief Vector with the probability to start in any of the subcompartments of the InfectedSevere compartment. struct StartingProbabilitiesInfectedSevere { - using Type = Vector; + using Type = Eigen::VectorX; /** * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedSevere compartment. * @param[in] numInfectedSevere Number of subcompartments of the InfectedSevere compartment. */ static Type get_default(size_t numInfectedSevere) { - Vector def = Vector::Zero(numInfectedSevere); - def[0] = 1.; + Eigen::VectorX def = Eigen::VectorX::Zero(numInfectedSevere); + def[0] = 1.; return def; } static std::string name() @@ -244,7 +248,8 @@ struct TransitionMatrixInfectedSevereToInfectedCritical { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -268,7 +273,8 @@ struct TransitionMatrixInfectedSevereToRecovered { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -280,15 +286,15 @@ struct TransitionMatrixInfectedSevereToRecovered { /// @brief Vector with the probability to start in any of the subcompartments of the InfectedCritical compartment. struct StartingProbabilitiesInfectedCritical { - using Type = Vector; + using Type = Eigen::VectorX; /** * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedCritical compartment. * @param[in] numInfectedCritical Number of subcompartments of the InfectedCritical compartment. */ static Type get_default(size_t numInfectedCritical) { - Vector def = Vector::Zero(numInfectedCritical); - def[0] = 1.; + Eigen::VectorX def = Eigen::VectorX::Zero(numInfectedCritical); + def[0] = 1.; return def; } static std::string name() @@ -311,7 +317,8 @@ struct TransitionMatrixInfectedCriticalToDead { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -335,7 +342,8 @@ struct TransitionMatrixInfectedCriticalToRecovered { */ static Type get_default(size_t dimension, ScalarType time = 1.) { - Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + Eigen::MatrixXd def = + Eigen::VectorX::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); def.diagonal(1).setConstant((ScalarType)dimension / time); return def; } @@ -574,7 +582,7 @@ class Parameters : public ParametersBase // --- Check that we have no flows back from one compartment to the previous one // (only in between of the subcompartments). --- if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning( @@ -583,7 +591,8 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones( + this->get().rows())) .array() > 1e-10) .any()) { log_warning("Constraint check: The entries of TransitionMatrixInfectedNoSymptomsToInfectedSymptoms lead to " @@ -592,7 +601,7 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning( @@ -601,7 +610,7 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning( @@ -610,7 +619,7 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning( @@ -619,7 +628,7 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning( @@ -628,7 +637,7 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning("Constraint check: The entries of TransitionMatrixInfectedSevereToRecovered lead to a negative " @@ -636,7 +645,7 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning("Constraint check: The entries of TransitionMatrixInfectedCriticalToDead lead to a negative " @@ -644,7 +653,7 @@ class Parameters : public ParametersBase return true; } if (((this->get() * - Vector::Ones(this->get().rows())) + Eigen::VectorX::Ones(this->get().rows())) .array() > 1e-10) .any()) { log_warning( diff --git a/cpp/models/ide_secir/model.cpp b/cpp/models/ide_secir/model.cpp index 0319c313b1..4112997c7e 100644 --- a/cpp/models/ide_secir/model.cpp +++ b/cpp/models/ide_secir/model.cpp @@ -49,14 +49,14 @@ Model::Model(TimeSeries&& init, CustomIndexArray 0) { // Add first time point in m_populations according to last time point in m_transitions which is where we start // the simulation. - m_populations.add_time_point>( + m_populations.add_time_point>( m_transitions.get_last_time(), TimeSeries::Vector::Constant((size_t)InfectionState::Count * m_num_agegroups, 0.)); } else { // Initialize m_populations with zero as the first point of time if no data is provided for the transitions. // This can happen for example in the case of initialization with real data. - m_populations.add_time_point>( + m_populations.add_time_point>( 0, TimeSeries::Vector::Constant((size_t)InfectionState::Count * m_num_agegroups, 0.)); } diff --git a/cpp/models/lct_secir/initializer_flows.h b/cpp/models/lct_secir/initializer_flows.h index 28193a3fba..27db250f9e 100644 --- a/cpp/models/lct_secir/initializer_flows.h +++ b/cpp/models/lct_secir/initializer_flows.h @@ -85,11 +85,12 @@ class Initializer * @return Returns true if one (or more) constraint(s) of the model, the initial flows * or the computed initial value vector are not satisfied, otherwise false. */ - bool compute_initialization_vector(Vector const& total_population, Vector const& deaths, - Vector const& total_confirmed_cases) + bool compute_initialization_vector(Eigen::VectorX const& total_population, + Eigen::VectorX const& deaths, + Eigen::VectorX const& total_confirmed_cases) { - Vector init(m_model.populations.get_compartments()); + Eigen::VectorX init(m_model.populations.get_compartments()); bool error = compute_initialization_vector_impl(init, total_population, deaths, total_confirmed_cases); if (error) { @@ -131,9 +132,10 @@ class Initializer * otherwise false. */ template - bool compute_initialization_vector_impl(Vector& init, Vector const& total_population, - Vector const& deaths, - Vector const& total_confirmed_cases) + bool compute_initialization_vector_impl(Eigen::VectorX& init, + Eigen::VectorX const& total_population, + Eigen::VectorX const& deaths, + Eigen::VectorX const& total_confirmed_cases) { static_assert((Group < Model::num_groups) && (Group >= 0), "The template parameter Group should be valid."); using LctStateGroup = type_at_index_t; @@ -254,7 +256,8 @@ class Initializer * are not satisfied, otherwise false. */ template - bool compute_compartment(Vector& init, Eigen::Index idx_incoming_flow, ScalarType transition_rate) const + bool compute_compartment(Eigen::VectorX& init, Eigen::Index idx_incoming_flow, + ScalarType transition_rate) const { size_t first_index = m_model.populations.template get_first_index_of_group() + type_at_index_t::template get_first_index(); diff --git a/cpp/models/lct_secir/model.h b/cpp/models/lct_secir/model.h index 3cccdf93a2..9a3d5a765b 100644 --- a/cpp/models/lct_secir/model.h +++ b/cpp/models/lct_secir/model.h @@ -85,8 +85,9 @@ class Model * @param[in] t The current time. * @param[out] dydt A reference to the calculated output. */ - void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, ScalarType t, - Eigen::Ref> dydt) const override + void get_derivatives(Eigen::Ref> pop, + Eigen::Ref> y, ScalarType t, + Eigen::Ref> dydt) const override { // Vectors are sorted such that we first have all InfectionState%s for AgeGroup 0, // afterwards all for AgeGroup 1 and so on. @@ -113,11 +114,11 @@ class Model TimeSeries compartments_ts(num_compartments); if (!(this->populations.get_num_compartments() == (size_t)subcompartments_ts.get_num_elements())) { log_error("Result does not match InfectionState of the Model."); - Vector wrong_size = Vector::Constant(num_compartments, -1); + Eigen::VectorX wrong_size = Eigen::VectorX::Constant(num_compartments, -1); compartments_ts.add_time_point(-1, wrong_size); return compartments_ts; } - Vector compartments(num_compartments); + Eigen::VectorX compartments(num_compartments); for (Eigen::Index timepoint = 0; timepoint < subcompartments_ts.get_num_time_points(); ++timepoint) { compress_vector(subcompartments_ts[timepoint], compartments); compartments_ts.add_time_point(subcompartments_ts.get_time(timepoint), compartments); @@ -148,7 +149,8 @@ class Model * @param[out] compartments Reference to the vector where the output is stored. */ template - void compress_vector(const Vector& subcompartments, Vector& compartments) const + void compress_vector(const Eigen::VectorX& subcompartments, + Eigen::VectorX& compartments) const { static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid."); using LctStateGroup = type_at_index_t; @@ -182,8 +184,9 @@ class Model * @param[out] dydt A reference to the calculated output. */ template - void get_derivatives_impl(Eigen::Ref> pop, Eigen::Ref> y, - ScalarType t, Eigen::Ref> dydt) const + void get_derivatives_impl(Eigen::Ref> pop, + Eigen::Ref> y, ScalarType t, + Eigen::Ref> dydt) const { static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid."); using LctStateGroup = type_at_index_t; @@ -294,8 +297,8 @@ class Model * @param[out] dydt A reference to the calculated output. */ template - void interact(Eigen::Ref> pop, Eigen::Ref> y, ScalarType t, - Eigen::Ref> dydt) const + void interact(Eigen::Ref> pop, Eigen::Ref> y, + ScalarType t, Eigen::Ref> dydt) const { static_assert((Group1 < num_groups) && (Group1 >= 0) && (Group2 < num_groups) && (Group2 >= 0), "The template parameters Group1 & Group2 should be valid."); diff --git a/cpp/models/lct_secir/parameters.h b/cpp/models/lct_secir/parameters.h index 099b6c22ab..579ea5e9eb 100644 --- a/cpp/models/lct_secir/parameters.h +++ b/cpp/models/lct_secir/parameters.h @@ -41,7 +41,7 @@ namespace lsecir * @brief Average time spent in the Exposed compartment for each group. */ struct TimeExposed { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -57,7 +57,7 @@ struct TimeExposed { * symptoms or recover for each group in the SECIR model in day unit. */ struct TimeInfectedNoSymptoms { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -73,7 +73,7 @@ struct TimeInfectedNoSymptoms { * or recover for each group in the SECIR model in day unit. */ struct TimeInfectedSymptoms { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -89,7 +89,7 @@ struct TimeInfectedSymptoms { * SECIR model in day unit. */ struct TimeInfectedSevere { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -104,7 +104,7 @@ struct TimeInfectedSevere { * @brief Average time treated by ICU before dead or recover for each group in the SECIR model in day unit. */ struct TimeInfectedCritical { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -119,7 +119,7 @@ struct TimeInfectedCritical { * @brief Probability of getting infected from a contact for each group. */ struct TransmissionProbabilityOnContact { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -152,7 +152,7 @@ struct ContactPatterns { * @brief The relative InfectedNoSymptoms infectability for each group. */ struct RelativeTransmissionNoSymptoms { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -167,7 +167,7 @@ struct RelativeTransmissionNoSymptoms { * @brief The risk of infection from symptomatic cases for each group in the SECIR model. */ struct RiskOfInfectionFromSymptomatic { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 1.); @@ -182,7 +182,7 @@ struct RiskOfInfectionFromSymptomatic { * @brief The percentage of asymptomatic cases for each group in the SECIR model. */ struct RecoveredPerInfectedNoSymptoms { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 0.5); @@ -197,7 +197,7 @@ struct RecoveredPerInfectedNoSymptoms { * @brief The percentage of hospitalized patients per infected patients for each group in the SECIR model. */ struct SeverePerInfectedSymptoms { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 0.5); @@ -212,7 +212,7 @@ struct SeverePerInfectedSymptoms { * @brief The percentage of ICU patients per hospitalized patients for each group in the SECIR model. */ struct CriticalPerSevere { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 0.5); @@ -227,7 +227,7 @@ struct CriticalPerSevere { * @brief The percentage of dead patients per ICU patients for each group in the SECIR model. */ struct DeathsPerCritical { - using Type = Vector>; + using Type = Eigen::VectorX>; static Type get_default(size_t size) { return Type::Constant(size, 1, 0.1); diff --git a/cpp/models/ode_seair/model.h b/cpp/models/ode_seair/model.h index 6ab003481e..fb5744ab4a 100644 --- a/cpp/models/ode_seair/model.h +++ b/cpp/models/ode_seair/model.h @@ -46,8 +46,8 @@ class Model : public mio::CompartmentalModel> pop, Eigen::Ref> y, FP /* t */, - Eigen::Ref> dydt) const override + void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, FP /* t */, + Eigen::Ref> dydt) const override { auto& params = this->parameters; const auto pop_total = pop.sum(); diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index abef23c534..594d64dede 100755 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -76,8 +76,8 @@ class Model : public FlowModel> pop, Eigen::Ref> y, FP t, - Eigen::Ref> flows) const override + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override { auto const& params = this->parameters; AgeGroup n_agegroups = params.get_num_groups(); @@ -253,7 +253,7 @@ class Simulation; * @tparam Base simulation type that uses a secir compartment model. see Simulation. */ template >> -double get_infections_relative(const Simulation& model, FP t, const Eigen::Ref>& y); +double get_infections_relative(const Simulation& model, FP t, const Eigen::Ref>& y); /** * specialization of compartment model simulation for secir models. @@ -283,7 +283,7 @@ class Simulation : public BaseT * @param tmax next stopping point of simulation * @return value at tmax */ - Eigen::Ref> advance(FP tmax) + Eigen::Ref> advance(FP tmax) { auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis(); auto& dyn_npis = this->get_model().parameters.template get>(); @@ -382,7 +382,8 @@ inline auto simulate_flows(FP t0, FP tmax, FP dt, const Model& model, //see declaration above. template -double get_infections_relative(const Simulation& sim, FP /* t*/, const Eigen::Ref>& y) +double get_infections_relative(const Simulation& sim, FP /* t*/, + const Eigen::Ref>& y) { double sum_inf = 0; for (auto i = AgeGroup(0); i < sim.get_model().parameters.get_num_groups(); ++i) { @@ -608,9 +609,9 @@ IOResult get_reproduction_number(size_t t_idx, const Simulation& s */ template -Vector get_reproduction_numbers(const Simulation& sim) +Eigen::VectorX get_reproduction_numbers(const Simulation& sim) { - Vector temp(sim.get_result().get_num_time_points()); + Eigen::VectorX temp(sim.get_result().get_num_time_points()); for (int i = 0; i < sim.get_result().get_num_time_points(); i++) { temp[i] = get_reproduction_number((size_t)i, sim).value(); } @@ -660,7 +661,7 @@ IOResult get_reproduction_number(ScalarType t_value, const Simulatio * @tparam Base simulation type that uses a secir compartment model; see Simulation. */ template , FP>> -auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) +auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { auto& params = sim.get_model().parameters; //parameters as arrays @@ -690,7 +691,7 @@ auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Re } template , FP>> -auto test_commuters(Simulation& sim, Eigen::Ref> mobile_population, FP time) +auto test_commuters(Simulation& sim, Eigen::Ref> mobile_population, FP time) { auto& model = sim.get_model(); auto nondetection = 1.0; diff --git a/cpp/models/ode_secirts/model.h b/cpp/models/ode_secirts/model.h index a7f22c20e4..d71bad7933 100644 --- a/cpp/models/ode_secirts/model.h +++ b/cpp/models/ode_secirts/model.h @@ -115,8 +115,8 @@ class Model { } - void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, - Eigen::Ref> flows) const override + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override { auto const& params = this->parameters; AgeGroup n_agegroups = params.get_num_groups(); @@ -265,7 +265,7 @@ class Model // effective contact rate by contact rate between groups i and j and damping j FP season_val = (1 + params.template get>() * sin(3.141592653589793 * - (std::fmod((params.template get() + t), 365.0) / 182.5 + 0.5))); + (std::fmod((params.template get() + t), 365.0) / 182.5 + 0.5))); FP cont_freq_eff = season_val * contact_matrix.get_matrix_at(t)(static_cast((size_t)i), static_cast((size_t)j)); // without died people @@ -676,7 +676,7 @@ class Simulation; * @tparam Base simulation type that uses the SECIRS-type compartment model. see Simulation. */ template >> -FP get_infections_relative(const Simulation& model, FP t, const Eigen::Ref>& y); +FP get_infections_relative(const Simulation& model, FP t, const Eigen::Ref>& y); /** * specialization of compartment model simulation for the SECIRTS model. @@ -852,7 +852,7 @@ inline auto simulate_flows(FP t0, FP tmax, FP dt, const Model& model, //see declaration above. template -FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) +FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { FP sum_inf = 0; for (auto i = AgeGroup(0); i < sim.get_model().parameters.get_num_groups(); ++i) { @@ -881,7 +881,7 @@ FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eige * @tparam Base simulation type that uses a SECIRS-type compartment model. see Simulation. */ template , FP>> -auto get_migration_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) +auto get_migration_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { auto& params = sim.get_model().parameters; //parameters as arrays @@ -929,7 +929,7 @@ auto get_migration_factors(const Simulation& sim, FP /*t*/, const Eigen::R * @param[in] time Current simulation time, used to determine the commuter detection period. */ template , FP>> -auto test_commuters(Simulation& sim, Eigen::Ref> migrated, FP time) +auto test_commuters(Simulation& sim, Eigen::Ref> migrated, FP time) { auto& model = sim.get_model(); auto nondetection = 1.0; diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 73c2e8cb33..2b58a1d403 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -106,8 +106,8 @@ class Model { } - void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, - Eigen::Ref> flows) const override + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override { auto const& params = this->parameters; AgeGroup n_agegroups = params.get_num_groups(); @@ -554,7 +554,7 @@ class Simulation; * @tparam Base simulation type that uses a secir compartment model. see Simulation. */ template >> -FP get_infections_relative(const Simulation& model, FP t, const Eigen::Ref>& y); +FP get_infections_relative(const Simulation& model, FP t, const Eigen::Ref>& y); /** * specialization of compartment model simulation for the SECIRVVS model. @@ -665,7 +665,7 @@ class Simulation : public BaseT * @param tmax next stopping point of simulation * @return value at tmax */ - Eigen::Ref> advance(FP tmax) + Eigen::Ref> advance(FP tmax) { auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis(); auto& dyn_npis = this->get_model().parameters.template get>(); @@ -786,7 +786,7 @@ inline auto simulate_flows(FP t0, FP tmax, FP dt, const Model& model, //see declaration above. template -FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) +FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { FP sum_inf = 0; for (auto i = AgeGroup(0); i < sim.get_model().parameters.get_num_groups(); ++i) { @@ -815,7 +815,7 @@ FP get_infections_relative(const Simulation& sim, FP /*t*/, const Eige * @tparam Base simulation type that uses a secir compartment model. see Simulation. */ template , FP>> -auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) +auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Ref>& y) { auto& params = sim.get_model().parameters; @@ -857,7 +857,7 @@ auto get_mobility_factors(const Simulation& sim, FP /*t*/, const Eigen::Re } template , FP>> -auto test_commuters(Simulation& sim, Eigen::Ref> mobile_population, FP time) +auto test_commuters(Simulation& sim, Eigen::Ref> mobile_population, FP time) { auto& model = sim.get_model(); auto nondetection = 1.0; diff --git a/cpp/models/ode_seir/model.h b/cpp/models/ode_seir/model.h index f15e7f3759..2021674560 100644 --- a/cpp/models/ode_seir/model.h +++ b/cpp/models/ode_seir/model.h @@ -68,8 +68,8 @@ class Model { } - void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, - Eigen::Ref> flows) const override + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override { const Index age_groups = reduce_index>(this->populations.size()); const auto& params = this->parameters; diff --git a/cpp/models/ode_sir/model.h b/cpp/models/ode_sir/model.h index 5afd87dc33..1142aeed80 100644 --- a/cpp/models/ode_sir/model.h +++ b/cpp/models/ode_sir/model.h @@ -58,8 +58,8 @@ class Model { } - void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, FP t, - Eigen::Ref> dydt) const override + void get_derivatives(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> dydt) const override { auto params = this->parameters; AgeGroup n_agegroups = params.get_num_groups(); diff --git a/cpp/models/sde_seirvv/model.h b/cpp/models/sde_seirvv/model.h index 0bc9363dbb..25062edd0f 100644 --- a/cpp/models/sde_seirvv/model.h +++ b/cpp/models/sde_seirvv/model.h @@ -56,102 +56,117 @@ class Model : public FlowModel> pop, Eigen::Ref> y, ScalarType t, - Eigen::Ref> flows) const + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, + ScalarType t, Eigen::Ref> flows) const { - auto& params = this->parameters; + auto& params = this->parameters; params.get().get_matrix_at(t)(0, 0); ScalarType coeffStoIV1 = params.get().get_matrix_at(t)(0, 0) * - params.get() / populations.get_total(); + params.get() / populations.get_total(); ScalarType coeffStoIV2 = params.get().get_matrix_at(t)(0, 0) * - params.get() / populations.get_total(); - - // Normal distributed values for the stochastic part of the flows, variables are encoded - // in the following way: x_y is the stochastic part for the flow from x to y. Variant + params.get() / populations.get_total(); + + // Normal distributed values for the stochastic part of the flows, variables are encoded + // in the following way: x_y is the stochastic part for the flow from x to y. Variant // specific compartments also get an addendum v1 or v2 denoting the relevant variant. - ScalarType s_ev1 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType s_ev2 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType ev1_iv1 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType ev2_iv2 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType iv1_rv1 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType iv2_rv2 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType rv1_ev1v2 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType ev1v2_iv1v2 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); - ScalarType iv1v2_rv1v2 = mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType s_ev1 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType s_ev2 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType ev1_iv1 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType ev2_iv2 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType iv1_rv1 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType iv2_rv2 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType rv1_ev1v2 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType ev1v2_iv1v2 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); + ScalarType iv1v2_rv1v2 = + mio::DistributionAdapter>::get_instance()(rng, 0.0, 1.0); // Assuming that no person can change its InfectionState twice in a single time step, // take the minimum of the calculated flow and the source compartment, to ensure that // no compartment attains negative values. // Calculate inv_step_size and inv_sqrt_step_size for optimization. - ScalarType inv_step_size = 1.0 / step_size; + ScalarType inv_step_size = 1.0 / step_size; ScalarType inv_sqrt_step_size = 1.0 / sqrt(step_size); // Two outgoing flows from S so will clamp their sum to S * inv_step_size to ensure non-negative S. const ScalarType outflow1 = std::clamp( coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1] + - sqrt(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1]) - * inv_sqrt_step_size * s_ev1, + sqrt(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1]) * + inv_sqrt_step_size * s_ev1, 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size); - const ScalarType outflow2 = std::clamp( - coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * - (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) + - sqrt(coeffStoIV2 * y[(size_t)InfectionState::Susceptible] * (pop[(size_t)InfectionState::InfectedV1V2] - + pop[(size_t)InfectionState::InfectedV2])) * inv_sqrt_step_size * s_ev2, - 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size); + const ScalarType outflow2 = + std::clamp(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * + (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) + + sqrt(coeffStoIV2 * y[(size_t)InfectionState::Susceptible] * + (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2])) * + inv_sqrt_step_size * s_ev2, + 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size); const ScalarType outflow_sum = outflow1 + outflow2; - if (outflow_sum > 0) - { - const ScalarType scale = std::clamp(outflow_sum, 0.0 , y[(size_t)InfectionState::Susceptible] * inv_step_size) / outflow_sum; + if (outflow_sum > 0) { + const ScalarType scale = + std::clamp(outflow_sum, 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size) / outflow_sum; flows[get_flat_flow_index()] = outflow1 * scale; flows[get_flat_flow_index()] = outflow2 * scale; - } else - { + } + else { flows[get_flat_flow_index()] = 0; - flows[get_flat_flow_index()] = 0; + flows[get_flat_flow_index()] = 0; } - flows[get_flat_flow_index()] = std::clamp( - (1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1] + - sqrt((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1]) * inv_sqrt_step_size * ev1_iv1, - 0.0, y[(size_t)InfectionState::ExposedV1] * inv_step_size); - - flows[get_flat_flow_index()] = std::clamp( - (1.0 / params.get()) * y[(size_t)InfectionState::ExposedV2] + - sqrt((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV2]) * inv_sqrt_step_size * ev2_iv2, - 0.0, y[(size_t)InfectionState::ExposedV2] * inv_step_size); - - flows[get_flat_flow_index()] = std::clamp( - (1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1] + - sqrt((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1]) * inv_sqrt_step_size * iv1_rv1, - 0.0, y[(size_t)InfectionState::InfectedV1] * inv_step_size); - - flows[get_flat_flow_index()] = std::clamp( - (1.0 / params.get()) * y[(size_t)InfectionState::InfectedV2] + - sqrt((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV2]) * inv_sqrt_step_size * iv2_rv2, - 0.0, y[(size_t)InfectionState::InfectedV2] * inv_step_size); - - flows[get_flat_flow_index()] = std::clamp( - coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] - * (pop[(size_t)InfectionState::InfectedV1V2]+ pop[(size_t)InfectionState::InfectedV2]) + - sqrt(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] * (pop[(size_t)InfectionState::InfectedV1V2] - + pop[(size_t)InfectionState::InfectedV2])) * inv_sqrt_step_size * rv1_ev1v2, - 0.0, y[(size_t)InfectionState::RecoveredV1] * inv_step_size); - - flows[get_flat_flow_index()] = std::clamp( - (1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1V2] + - sqrt((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1V2]) / - sqrt(step_size) * ev1v2_iv1v2, - 0.0, y[(size_t)InfectionState::ExposedV1V2] * inv_step_size); - - flows[get_flat_flow_index()] = std::clamp( - (1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1V2] + - sqrt((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1V2]) / - sqrt(step_size) * iv1v2_rv1v2, - 0.0, y[(size_t)InfectionState::InfectedV1V2] * inv_step_size); + flows[get_flat_flow_index()] = + std::clamp((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1] + + sqrt((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1]) * + inv_sqrt_step_size * ev1_iv1, + 0.0, y[(size_t)InfectionState::ExposedV1] * inv_step_size); + + flows[get_flat_flow_index()] = + std::clamp((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV2] + + sqrt((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV2]) * + inv_sqrt_step_size * ev2_iv2, + 0.0, y[(size_t)InfectionState::ExposedV2] * inv_step_size); + + flows[get_flat_flow_index()] = + std::clamp((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1] + + sqrt((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1]) * + inv_sqrt_step_size * iv1_rv1, + 0.0, y[(size_t)InfectionState::InfectedV1] * inv_step_size); + + flows[get_flat_flow_index()] = + std::clamp((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV2] + + sqrt((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV2]) * + inv_sqrt_step_size * iv2_rv2, + 0.0, y[(size_t)InfectionState::InfectedV2] * inv_step_size); + + flows[get_flat_flow_index()] = + std::clamp(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] * + (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) + + sqrt(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] * + (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2])) * + inv_sqrt_step_size * rv1_ev1v2, + 0.0, y[(size_t)InfectionState::RecoveredV1] * inv_step_size); + + flows[get_flat_flow_index()] = + std::clamp((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1V2] + + sqrt((1.0 / params.get()) * y[(size_t)InfectionState::ExposedV1V2]) / + sqrt(step_size) * ev1v2_iv1v2, + 0.0, y[(size_t)InfectionState::ExposedV1V2] * inv_step_size); + + flows[get_flat_flow_index()] = + std::clamp((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1V2] + + sqrt((1.0 / params.get()) * y[(size_t)InfectionState::InfectedV1V2]) / + sqrt(step_size) * iv1v2_rv1v2, + 0.0, y[(size_t)InfectionState::InfectedV1V2] * inv_step_size); } ScalarType step_size; ///< A step size of the model with which the stochastic process is realized. diff --git a/cpp/models/sde_sir/model.h b/cpp/models/sde_sir/model.h index f35b65a5b2..dc564e6a7a 100644 --- a/cpp/models/sde_sir/model.h +++ b/cpp/models/sde_sir/model.h @@ -49,9 +49,8 @@ class Model : public FlowModel> pop, - Eigen::Ref> y, ScalarType t, - Eigen::Ref> flows) const + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, + ScalarType t, Eigen::Ref> flows) const { auto& params = this->parameters; ScalarType coeffStoI = params.get().get_matrix_at(t)(0, 0) * diff --git a/cpp/models/sde_sir/simulation.h b/cpp/models/sde_sir/simulation.h index 957103b86d..89eec3bff2 100644 --- a/cpp/models/sde_sir/simulation.h +++ b/cpp/models/sde_sir/simulation.h @@ -56,7 +56,7 @@ class Simulation : public mio::Simulation * tmax must be greater than get_result().get_last_time_point() * @param tmax next stopping point of simulation */ - Eigen::Ref> advance(ScalarType tmax) + Eigen::Ref> advance(ScalarType tmax) { return get_ode_integrator().advance( [this](auto&& y, auto&& t, auto&& dydt) { @@ -108,7 +108,7 @@ class FlowSimulation : public mio::FlowSimulation * tmax must be greater than get_result().get_last_time_point(). * @param[in] tmax Next stopping time of the simulation. */ - Eigen::Ref> advance(ScalarType tmax) + Eigen::Ref> advance(ScalarType tmax) { assert(get_flows().get_num_time_points() == get_result().get_num_time_points()); auto result = this->get_ode_integrator().advance( diff --git a/cpp/models/sde_sirs/model.h b/cpp/models/sde_sirs/model.h index 69e8d26d63..66c3b59d9b 100644 --- a/cpp/models/sde_sirs/model.h +++ b/cpp/models/sde_sirs/model.h @@ -50,8 +50,8 @@ class Model : public FlowModel> pop, Eigen::Ref> y, ScalarType t, - Eigen::Ref> flows) const + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, + ScalarType t, Eigen::Ref> flows) const { auto& params = this->parameters; ScalarType coeffStoI = params.get().get_matrix_at(t)(0, 0) * diff --git a/cpp/models/sde_sirs/simulation.h b/cpp/models/sde_sirs/simulation.h index 0b07b2243b..c26aaa44b9 100644 --- a/cpp/models/sde_sirs/simulation.h +++ b/cpp/models/sde_sirs/simulation.h @@ -56,7 +56,7 @@ class Simulation : public mio::Simulation * tmax must be greater than get_result().get_last_time_point() * @param tmax next stopping point of simulation */ - Eigen::Ref> advance(ScalarType tmax) + Eigen::Ref> advance(ScalarType tmax) { return get_ode_integrator().advance( [this](auto&& y, auto&& t, auto&& dydt) { @@ -108,7 +108,7 @@ class FlowSimulation : public mio::FlowSimulation * tmax must be greater than get_result().get_last_time_point(). * @param[in] tmax Next stopping time of the simulation. */ - Eigen::Ref> advance(ScalarType tmax) + Eigen::Ref> advance(ScalarType tmax) { assert(get_flows().get_num_time_points() == get_result().get_num_time_points()); auto result = this->get_ode_integrator().advance( diff --git a/cpp/tests/test_glct_secir.cpp b/cpp/tests/test_glct_secir.cpp index 6b7b9c4a60..33c52e5695 100644 --- a/cpp/tests/test_glct_secir.cpp +++ b/cpp/tests/test_glct_secir.cpp @@ -53,7 +53,7 @@ TEST(TestGLCTSecir, testEvalRightHandSide) mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default( LctState::get_num_subcompartments(), 3.2); // InfectedNoSymptoms. - mio::Vector StartingProbabilitiesInfectedNoSymptoms = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedNoSymptoms = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedNoSymptoms[0] = 1 - 0.09; StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)( @@ -67,7 +67,7 @@ TEST(TestGLCTSecir, testEvalRightHandSide) mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), 2.); // InfectedSymptoms. - mio::Vector StartingProbabilitiesInfectedSymptoms = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedSymptoms = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedSymptoms[0] = 0.2; StartingProbabilitiesInfectedSymptoms[(Eigen::Index)( @@ -80,7 +80,7 @@ TEST(TestGLCTSecir, testEvalRightHandSide) mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), 5.8); // InfectedSevere. - mio::Vector StartingProbabilitiesInfectedSevere = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedSevere = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedSevere[0] = 0.25; StartingProbabilitiesInfectedSevere[(Eigen::Index)( @@ -93,7 +93,7 @@ TEST(TestGLCTSecir, testEvalRightHandSide) mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), 9.5); // InfectedCritical. - mio::Vector StartingProbabilitiesInfectedCritical = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedCritical = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedCritical[0] = 0.3; StartingProbabilitiesInfectedCritical[(Eigen::Index)( @@ -135,16 +135,16 @@ TEST(TestGLCTSecir, testEvalRightHandSide) for (auto&& vec : initial_populations) { flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end()); } - mio::Vector pop(LctState::Count); + Eigen::VectorX pop(LctState::Count); for (size_t i = 0; i < LctState::Count; i++) { pop[i] = flat_initial_populations[i]; } // Compare the result of get_derivatives() with a hand calculated result. - mio::Vector dydt(LctState::Count); + Eigen::VectorX dydt(LctState::Count); model.get_derivatives(pop, pop, 0, dydt); // This vector is the equivalent of the result defined in the test suite testEvalRightHandSide of the LCT model. - mio::Vector compare(LctState::Count); + Eigen::VectorX compare(LctState::Count); compare << -15.3409, -3.4091, 6.25, -17.5 * 0.91, 15 * 0.91, 0 * 0.91, -17.5 * 0.09, 15 * 0.09, 0 * 0.09, 3.3052 * 0.2, 3.4483 * 0.2, 3.3052 * 0.8, 3.4483 * 0.8, -7.0417 * 0.25, 6.3158 * 0.25, -7.0417 * 0.75, 6.3158 * 0.75, -2.2906 * 0.3, -2.8169 * 0.3, -2.2906 * 0.7, -2.8169 * 0.7, 12.3899, 1.6901; @@ -175,7 +175,7 @@ class ModelTestGLCTSecir : public testing::Test mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default( LctState::get_num_subcompartments(), 3.2); // InfectedNoSymptoms. - mio::Vector StartingProbabilitiesInfectedNoSymptoms = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedNoSymptoms = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedNoSymptoms[0] = 1 - 0.09; StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)( @@ -189,7 +189,7 @@ class ModelTestGLCTSecir : public testing::Test mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), 2.); // InfectedSymptoms. - mio::Vector StartingProbabilitiesInfectedSymptoms = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedSymptoms = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedSymptoms[0] = 0.2; StartingProbabilitiesInfectedSymptoms[(Eigen::Index)( @@ -203,7 +203,7 @@ class ModelTestGLCTSecir : public testing::Test mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), 5.8); // InfectedSevere. - mio::Vector StartingProbabilitiesInfectedSevere = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedSevere = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedSevere[0] = 0.25; StartingProbabilitiesInfectedSevere[(Eigen::Index)( @@ -217,7 +217,7 @@ class ModelTestGLCTSecir : public testing::Test mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( (size_t)(LctState::get_num_subcompartments() / 2.), 9.5); // InfectedCritical. - mio::Vector StartingProbabilitiesInfectedCritical = mio::Vector::Zero( + Eigen::VectorX StartingProbabilitiesInfectedCritical = Eigen::VectorX::Zero( (Eigen::Index)LctState::get_num_subcompartments()); StartingProbabilitiesInfectedCritical[0] = 0.3; StartingProbabilitiesInfectedCritical[(Eigen::Index)( @@ -314,8 +314,8 @@ TEST_F(ModelTestGLCTSecir, testConstraintsModel) EXPECT_FALSE(constraint_check); // Check if the number of subcompartments does not match the dimension of the vector with StartingProbabilities. - mio::Vector wrong_size = mio::Vector::Zero(3); - wrong_size[0] = 1; + Eigen::VectorX wrong_size = Eigen::VectorX::Zero(3); + wrong_size[0] = 1; // Exposed. model->parameters.get() = wrong_size; constraint_check = model->check_constraints(); @@ -398,8 +398,8 @@ TEST_F(ModelTestGLCTSecir, testConstraintsParameters) (size_t)(LctState::get_num_subcompartments() / 2.), 7.1); // Check non matching dimensions of TransitionMatrix and vector with StartingProbabilities. - mio::Vector wrong_size = mio::Vector::Zero(3); - wrong_size[0] = 1; + Eigen::VectorX wrong_size = Eigen::VectorX::Zero(3); + wrong_size[0] = 1; // Exposed. model->parameters.get() = wrong_size; constraint_check = model->parameters.check_constraints(); @@ -512,7 +512,7 @@ TEST_F(ModelTestGLCTSecir, testCalculatePopWrongSize) size_t wrong_size = LctState::Count - 2; // Define TimeSeries with wrong_size elements. mio::TimeSeries wrong_num_elements(wrong_size); - mio::Vector vec_wrong_size = mio::Vector::Ones(wrong_size); + Eigen::VectorX vec_wrong_size = Eigen::VectorX::Ones(wrong_size); wrong_num_elements.add_time_point(-10, vec_wrong_size); wrong_num_elements.add_time_point(-9, vec_wrong_size); // Call the calculate_compartments function with the TimeSeries with a wrong number of elements. diff --git a/cpp/tests/test_ide_parameters_io.cpp b/cpp/tests/test_ide_parameters_io.cpp index 6da9bece05..dcad2031f9 100644 --- a/cpp/tests/test_ide_parameters_io.cpp +++ b/cpp/tests/test_ide_parameters_io.cpp @@ -106,7 +106,7 @@ TEST(TestIDEParametersIo, RKIcompareWithPreviousRun) } // Compare transitions at last time point with results from a previous run that are given here. - mio::Vector compare(num_transitions * num_agegroups); + Eigen::VectorX compare(num_transitions * num_agegroups); compare << 336.428571428600, 328.285714285701, 162.000000000000, 163.071428571425, 80.130989648839, 79.803571428575, 39.476374533415, 39.476374533415, 19.550404043081, 19.550404043081, 1105.714285714297, 1069.857142857200, 515.714285714250, 163.071428571425, 80.130989648839, 79.803571428575, 39.476374533415, 39.476374533415, diff --git a/cpp/tests/test_ide_secir_ageres.cpp b/cpp/tests/test_ide_secir_ageres.cpp index d05c26809d..e7048f4181 100644 --- a/cpp/tests/test_ide_secir_ageres.cpp +++ b/cpp/tests/test_ide_secir_ageres.cpp @@ -128,7 +128,7 @@ TEST(TestIdeAgeres, compareWithPreviousRun) // Compare compartments at last time point with results from a previous run that are given here. mio::TimeSeries compartments = sim.get_result(); - mio::Vector compare_compartments(num_compartments * num_agegroups); + Eigen::VectorX compare_compartments(num_compartments * num_agegroups); compare_compartments << 484.3056557672, 15.7685031055, 22.7020934123, 7.0615933479, 3.3491460693, 1.5803397070, 4454.5548070034, 10.6778615873, 484.3056557672, 31.0934010790, 21.1271954388, 23.6370809253, 3.9106794140, 1.7242153411, 4424.0110181177, 10.1907539167, 605.3820697090, 60.1973290710, 23.8046231705, 16.6085494134, @@ -143,7 +143,7 @@ TEST(TestIdeAgeres, compareWithPreviousRun) // Compare transitions at last time point with results from a previous run that are given here. mio::TimeSeries transitions = sim.get_transitions(); - mio::Vector compare_transitions(num_transitions * num_agegroups); + Eigen::VectorX compare_transitions(num_transitions * num_agegroups); compare_transitions << 31.5370062111, 30.6497959470, 14.1231866958, 14.7543908776, 6.6982921386, 6.6982921386, 3.1606794140, 3.1606794140, 1.4742153411, 1.4742153411, 31.5370062111, 29.5087817552, 14.7543908776, 14.1231866958, 6.3213588280, 6.3213588280, 2.9484306823, 2.9484306823, 1.3533839877, 1.3533839877, @@ -155,4 +155,4 @@ TEST(TestIdeAgeres, compareWithPreviousRun) for (int j = 0; j < compare_transitions.size(); j++) { ASSERT_NEAR(transitions.get_last_value()[j], compare_transitions[j], 1e-7); } -} \ No newline at end of file +} diff --git a/cpp/tests/test_lct_initializer_flows.cpp b/cpp/tests/test_lct_initializer_flows.cpp index 48cb71acf1..44357bca6c 100644 --- a/cpp/tests/test_lct_initializer_flows.cpp +++ b/cpp/tests/test_lct_initializer_flows.cpp @@ -39,7 +39,7 @@ TEST(TestInitializer, compareWithPrevious) using Model = mio::lsecir::Model; // Previous result. - mio::Vector compare(LctState::Count); + Eigen::VectorX compare(LctState::Count); compare << 82810889.00545, 850.70432, 970.04980, 315.32890, 391.51799, 391.39351, 565.45854, 580.79267, 85.97421, 86.02738, 80.26791, 189.53449, 167.57963, 329757.36512, 9710; @@ -66,9 +66,9 @@ TEST(TestInitializer, compareWithPrevious) model.parameters.get()[0] = 0.1; model.parameters.get()[0] = 0.1; - mio::Vector total_confirmed_cases = mio::Vector::Constant(1, 341223.); - mio::Vector deaths = mio::Vector::Constant(1, 9710.); - mio::Vector total_population = mio::Vector::Constant(1, 83155031.); + Eigen::VectorX total_confirmed_cases = Eigen::VectorX::Constant(1, 341223.); + Eigen::VectorX deaths = Eigen::VectorX::Constant(1, 9710.); + Eigen::VectorX total_population = Eigen::VectorX::Constant(1, 83155031.); // Add time points for initialization of transitions. mio::TimeSeries init((int)mio::lsecir::InfectionTransition::Count); @@ -110,7 +110,7 @@ TEST(TestInitializer, compareWithPreviousThreeGroups) using Model = mio::lsecir::Model; // Previous result. - mio::Vector compare(LctState::Count); + Eigen::VectorX compare(LctState::Count); compare << 82810889.00545, 850.70432, 970.04980, 315.32890, 391.51799, 391.39351, 565.45854, 580.79267, 85.97421, 86.02738, 80.26791, 189.53449, 167.57963, 329757.36512, 9710; @@ -139,11 +139,12 @@ TEST(TestInitializer, compareWithPreviousThreeGroups) mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, 10.)); - mio::Vector total_confirmed_cases = - mio::Vector::Constant(num_groups, 341223. / (ScalarType)num_groups); - mio::Vector deaths = mio::Vector::Constant(num_groups, 9710. / (ScalarType)num_groups); - mio::Vector total_population = - mio::Vector::Constant(num_groups, 83155031. / (ScalarType)num_groups); + Eigen::VectorX total_confirmed_cases = + Eigen::VectorX::Constant(num_groups, 341223. / (ScalarType)num_groups); + Eigen::VectorX deaths = + Eigen::VectorX::Constant(num_groups, 9710. / (ScalarType)num_groups); + Eigen::VectorX total_population = + Eigen::VectorX::Constant(num_groups, 83155031. / (ScalarType)num_groups); // Add time points for initialization of transitions. mio::TimeSeries init(num_groups * (size_t)mio::lsecir::InfectionTransition::Count); @@ -190,10 +191,10 @@ TEST(TestInitializer, testConstraints) // Deactivate temporarily log output for next tests. mio::set_log_level(mio::LogLevel::off); - ScalarType dt = 0.5; - mio::Vector total_confirmed_cases = mio::Vector::Constant(2, 341223.); - mio::Vector deaths = mio::Vector::Constant(2, 9710.); - mio::Vector total_population = mio::Vector::Constant(2, 83155031.); + ScalarType dt = 0.5; + Eigen::VectorX total_confirmed_cases = Eigen::VectorX::Constant(2, 341223.); + Eigen::VectorX deaths = Eigen::VectorX::Constant(2, 9710.); + Eigen::VectorX total_population = Eigen::VectorX::Constant(2, 83155031.); // Use a model with two groups. using InfState = mio::lsecir::InfectionState; using LctState = mio::LctInfectionState; @@ -205,7 +206,7 @@ TEST(TestInitializer, testConstraints) // Check wrong size of initial flows. mio::TimeSeries init_wrong_size(infectionTransition_count - 1); - mio::Vector vec_wrong_size = mio::Vector::Ones(infectionTransition_count - 1); + Eigen::VectorX vec_wrong_size = Eigen::VectorX::Ones(infectionTransition_count - 1); init_wrong_size.add_time_point(-50, vec_wrong_size); while (init_wrong_size.get_last_time() < 0) { init_wrong_size.add_time_point(init_wrong_size.get_last_time() + dt, vec_wrong_size); @@ -219,7 +220,7 @@ TEST(TestInitializer, testConstraints) // Check if last time of initial flows is not zero. mio::TimeSeries init_wrong(infectionTransition_count); - mio::Vector vec_init = mio::Vector::Ones(infectionTransition_count); + Eigen::VectorX vec_init = Eigen::VectorX::Ones(infectionTransition_count); init_wrong.add_time_point(-50, vec_init); while (init_wrong.get_last_time() < -5) { init_wrong.add_time_point(init_wrong.get_last_time() + dt, vec_init); @@ -334,7 +335,7 @@ TEST(TestInitializer, testConstraints) // Check with negative result for deaths. mio::TimeSeries init_negative_deaths(infectionTransition_count); vec_init[(int)mio::lsecir::InfectionTransition::InfectedSevereToInfectedCritical] = 1; - deaths = mio::Vector::Constant(2, -100.); + deaths = Eigen::VectorX::Constant(2, -100.); init_negative_deaths.add_time_point(-50., vec_init); while (init_negative_deaths.get_last_time() < 0) { init_negative_deaths.add_time_point(init_negative_deaths.get_last_time() + dt, vec_init); @@ -345,7 +346,7 @@ TEST(TestInitializer, testConstraints) EXPECT_TRUE(status); // Check with correct initialization values. - deaths = mio::Vector::Constant(2, 100.); + deaths = Eigen::VectorX::Constant(2, 100.); mio::TimeSeries init_right(infectionTransition_count); init_right.add_time_point(-50, vec_init); while (init_right.get_last_time() < 0) { diff --git a/cpp/tests/test_lct_parameters_io.cpp b/cpp/tests/test_lct_parameters_io.cpp index df94f36522..ab1873aa70 100644 --- a/cpp/tests/test_lct_parameters_io.cpp +++ b/cpp/tests/test_lct_parameters_io.cpp @@ -58,8 +58,8 @@ std::vector get_synthetic_rki_data_age() const int num_agegroups = 6; Json::Value js(Json::arrayValue); std::vector dates = {"2020-05-26", "2020-05-27", "2020-05-28", "2020-05-29", - "2020-05-30", "2020-05-31", "2020-06-01", "2020-06-02", - "2020-06-03", "2020-06-04", "2020-06-05"}; + "2020-05-30", "2020-05-31", "2020-06-01", "2020-06-02", + "2020-06-03", "2020-06-04", "2020-06-05"}; std::vector age_group_names = {"A00-A04", "A05-A14", "A15-A34", "A35-A59", "A60-A79", "A80+"}; for (int day = 0; day < 11; day++) { for (int age = 0; age < num_agegroups; age++) { @@ -106,7 +106,7 @@ TEST(TestLCTParametersIo, ReadPopulationDataRKI) ASSERT_THAT(print_wrap(read_result), IsSuccess()); // Result to compare the simulation result with. - mio::Vector compare((Eigen::Index)LctState::Count); + Eigen::VectorX compare((Eigen::Index)LctState::Count); // Calculate result using that the number of new confirmed cases at each day is one in the synthetic case data // and additionally the stay times, the number of subcompartments, and the transition probabilities. compare << 889.5, 1. / (2. * 0.8) * 2.3, 1. / (2. * 0.8) * 2.3, 1. / (3. * 0.8) * 1.3, 1. / (3. * 0.8) * 1.3, @@ -156,7 +156,7 @@ TEST(TestLCTParametersIo, ReadPopulationDataRKIAgeres) // in the synthetic case data and additionally the stay times, the number of subcompartments, // and the transition probabilities. size_t num_populations = (size_t)InfState::Count * num_agegroups; - mio::Vector compare(num_populations); + Eigen::VectorX compare(num_populations); for (size_t age = 0; age < num_agegroups; age++) { compare[(size_t)InfState::Count * age + (size_t)InfState::Exposed] = 1. / (1. - model.parameters.get()[age]) * diff --git a/cpp/tests/test_lct_secir.cpp b/cpp/tests/test_lct_secir.cpp index 9f0e0a999e..8b2d9e90ab 100644 --- a/cpp/tests/test_lct_secir.cpp +++ b/cpp/tests/test_lct_secir.cpp @@ -46,10 +46,10 @@ TEST(TestLCTSecir, simulateDefault) ScalarType tmax = 1; ScalarType dt = 0.1; - mio::Vector init = mio::Vector::Constant((Eigen::Index)InfState::Count, 15); - init[0] = 200; - init[3] = 50; - init[5] = 30; + Eigen::VectorX init = Eigen::VectorX::Constant((Eigen::Index)InfState::Count, 15); + init[0] = 200; + init[3] = 50; + init[5] = 30; Model model; for (size_t i = 0; i < LctState::Count; i++) { @@ -77,10 +77,10 @@ TEST(TestLCTSecir, compareWithOdeSecir) ScalarType dt = 0.1; // Initialization vector for both models. - mio::Vector init = mio::Vector::Constant((Eigen::Index)InfState::Count, 15); - init[0] = 200; - init[3] = 50; - init[5] = 30; + Eigen::VectorX init = Eigen::VectorX::Constant((Eigen::Index)InfState::Count, 15); + init[0] = 200; + init[3] = 50; + init[5] = 30; // Define LCT model. Model model_lct; @@ -239,10 +239,10 @@ TEST(TestLCTSecir, testEvalRightHandSide) // Compare the result of get_derivatives() with a hand calculated result. size_t num_subcompartments = LctState::Count; - mio::Vector dydt(num_subcompartments); + Eigen::VectorX dydt(num_subcompartments); model.get_derivatives(model.get_initial_values(), model.get_initial_values(), 0, dydt); - mio::Vector compare(num_subcompartments); + Eigen::VectorX compare(num_subcompartments); compare << -15.3409, -3.4091, 6.25, -17.5, 15, 0, 3.3052, 3.4483, -7.0417, 6.3158, -2.2906, -2.8169, 12.3899, 1.6901; @@ -316,10 +316,10 @@ TEST(TestLCTSecir, testEvalRightHandSideGroups) model.parameters.get()[1] = 0.5; // Compare the result of get_derivatives() with a hand calculated result. - mio::Vector dydt(num_subcompartments); + Eigen::VectorX dydt(num_subcompartments); model.get_derivatives(model.get_initial_values(), model.get_initial_values(), 0, dydt); - mio::Vector compare(num_subcompartments); + Eigen::VectorX compare(num_subcompartments); compare << -90.6818, 78.6818, 4, -4, 6, 0, -6.6, 4, -15.2, 12, -3.6, -4, 18.6, 0.8, -90.6818, 85.6818, -20, 12, -4.25, 2.25, 15, 0; @@ -333,7 +333,7 @@ TEST(TestLCTSecir, testEvalRightHandSideGroups) result.add_time_point(1, model.get_initial_values() + dydt); mio::TimeSeries population = model.calculate_compartments(result); // Sum of subcompartments in initial_population. - mio::Vector compare_population0(2 * (size_t)InfState::Count); + Eigen::VectorX compare_population0(2 * (size_t)InfState::Count); compare_population0 << 750, 50, 40, 50, 50, 30, 20, 10, 750, 10, 50, 1, 9, 0, 30, 100; ASSERT_EQ(compare_population0.size(), static_cast(population.get_num_elements())); EXPECT_NEAR(result.get_time(0), population.get_time(0), 1e-7); @@ -341,7 +341,7 @@ TEST(TestLCTSecir, testEvalRightHandSideGroups) EXPECT_NEAR(compare_population0[i], population.get_value(0)[i], 1e-3) << " at index " << i << ".\n"; } // Sum of subcompartments in compare vector. - mio::Vector compare_population1(2 * (size_t)InfState::Count); + Eigen::VectorX compare_population1(2 * (size_t)InfState::Count); compare_population1 << -90.6818, 82.6818, 2, -2.6, -3.2, -7.6, 18.6, 0.8, -90.6818, 85.6818, -20, 12, -4.25, 2.25, 15, 0; compare_population1 = compare_population0 + compare_population1; @@ -403,11 +403,11 @@ TEST(TestLCTSecir, testEvalRightHandSideThreeGroupsEqual) // Compare the result of get_derivatives() with a hand calculated result. - mio::Vector dydt(num_groups * num_subcompartments); + Eigen::VectorX dydt(num_groups * num_subcompartments); model.get_derivatives(model.get_initial_values(), model.get_initial_values(), 0, dydt); - mio::Vector compare(num_subcompartments); + Eigen::VectorX compare(num_subcompartments); compare << -15.3409, -3.4091, 6.25, -17.5, 15, 0, 3.3052, 3.4483, -7.0417, 6.3158, -2.2906, -2.8169, 12.3899, 1.6901; ScalarType sunum_groups = 0; @@ -601,7 +601,7 @@ TEST_F(ModelTestLCTSecir, testCalculatePopWrongSize) size_t wrong_size = LctState::Count - 2; // Define TimeSeries with wrong_size elements. mio::TimeSeries wrong_num_elements(wrong_size); - mio::Vector vec_wrong_size = mio::Vector::Ones(wrong_size); + Eigen::VectorX vec_wrong_size = Eigen::VectorX::Ones(wrong_size); wrong_num_elements.add_time_point(-10, vec_wrong_size); wrong_num_elements.add_time_point(-9, vec_wrong_size); // Call the calculate_compartments function with the TimeSeries with a wrong number of elements. From 00fd88564e90f1df6d85bf5ad7a2b4a86359b38c Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 3 Jan 2025 16:36:54 +0100 Subject: [PATCH 14/16] include config.h directly --- cpp/tests/matchers.h | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/tests/matchers.h b/cpp/tests/matchers.h index 4953625f75..3dcfd0adfd 100644 --- a/cpp/tests/matchers.h +++ b/cpp/tests/matchers.h @@ -20,6 +20,7 @@ #ifndef EPI_TESTS_MATCHERS_H #define EPI_TESTS_MATCHERS_H +#include "memilio/config.h" #include "memilio/utils/compiler_diagnostics.h" #include "memilio/math/floating_point.h" #include "memilio/io/io.h" From 0b5b77a344fbb4371c1179f9aff0fd352504fb18 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 3 Jan 2025 16:44:41 +0100 Subject: [PATCH 15/16] [ci skip] fix type in comment --- cpp/memilio/utils/logging.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/memilio/utils/logging.h b/cpp/memilio/utils/logging.h index f6ef1dc69b..b02649913e 100644 --- a/cpp/memilio/utils/logging.h +++ b/cpp/memilio/utils/logging.h @@ -133,7 +133,7 @@ namespace internal { /** - * @brief Format AD types (like ad::gt1s) using their value for logging with spdlog. + * @brief Format AD types (like ad::gt1s::type) using their value for logging with spdlog. * * If derivative information is needed as well, use `ad::derivative(...)` or define a `fmt::formatter<...>`. */ From af94264bc9839da8a69ea8462e33ecdaf1d292e2 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Mon, 6 Jan 2025 14:01:12 +0100 Subject: [PATCH 16/16] use resize_destructive as the new resize --- cpp/memilio/utils/custom_index_array.h | 27 +++++++++----------------- cpp/models/abm/model.cpp | 8 +++----- cpp/tests/abm_helpers.cpp | 8 +++----- 3 files changed, 15 insertions(+), 28 deletions(-) diff --git a/cpp/memilio/utils/custom_index_array.h b/cpp/memilio/utils/custom_index_array.h index aa7b4a0a90..fcd3ed51ed 100644 --- a/cpp/memilio/utils/custom_index_array.h +++ b/cpp/memilio/utils/custom_index_array.h @@ -235,24 +235,12 @@ class CustomIndexArray } /** - * @brief Resize all dimensions. - * Note that when increasing the overall size, new values may be uninitialized. + * @brief Resize all dimensions, invalidating entries. + * All entries of the CustomIndexArray should be reassigned a new value after a resize, as it may delete or reorder + * entries in an unexpected way. Newly added entries are default constructed. * @param new_dims New dimensions. */ void resize(Index new_dims) - { - m_dimensions = new_dims; - m_numel = product(m_dimensions); - m_y.conservativeResize(m_numel); - } - - /** - * @brief Resize all dimensions, destroying all values. - * This version of resize should only be used when the CustomIndexArray contains non-movable and non-copyable - * values, like atomics. New entries are all default initialized. - * @param new_dims New dimensions. - */ - void resize_destructive(Index new_dims) { m_dimensions = new_dims; m_numel = product(m_dimensions); @@ -260,15 +248,18 @@ class CustomIndexArray } /** - * Resize a single dimension. - * @param new dimension. + * Resize a single dimension, invalidating entries. + * All entries of the CustomIndexArray should be reassigned a new value after a resize, as it may delete or reorder + * entries in an unexpected way. Newly added entries are default constructed. + * @param new_dim New dimension size. + * @tparam Tag The dimension to resize. */ template void resize(mio::Index new_dim) { std::get>(m_dimensions.indices) = new_dim; m_numel = product(m_dimensions); - m_y.conservativeResize(m_numel); + m_y.resize(m_numel); } /** diff --git a/cpp/models/abm/model.cpp b/cpp/models/abm/model.cpp index 54a1e8f98f..4c8ba1f7f5 100755 --- a/cpp/models/abm/model.cpp +++ b/cpp/models/abm/model.cpp @@ -225,11 +225,9 @@ void Model::build_exposure_caches() m_contact_exposure_rates_cache.resize(num_locations); PRAGMA_OMP(taskloop) for (size_t i = 0; i < num_locations; i++) { - m_air_exposure_rates_cache[i].resize_destructive( - {CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count}); - m_contact_exposure_rates_cache[i].resize_destructive({CellIndex(m_locations[i].get_cells().size()), - VirusVariant::Count, - AgeGroup(parameters.get_num_groups())}); + m_air_exposure_rates_cache[i].resize({CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count}); + m_contact_exposure_rates_cache[i].resize({CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count, + AgeGroup(parameters.get_num_groups())}); } // implicit taskloop barrier m_are_exposure_caches_valid = false; m_exposure_caches_need_rebuild = false; diff --git a/cpp/tests/abm_helpers.cpp b/cpp/tests/abm_helpers.cpp index 3ff3206841..8ab1e79753 100644 --- a/cpp/tests/abm_helpers.cpp +++ b/cpp/tests/abm_helpers.cpp @@ -49,16 +49,14 @@ void interact_testing(mio::abm::PersonalRandomNumberGenerator& personal_rng, mio { // allocate and initialize air exposures with 0 mio::abm::AirExposureRates local_air_exposure; - local_air_exposure.resize_destructive( - {mio::abm::CellIndex(location.get_cells().size()), mio::abm::VirusVariant::Count}); + local_air_exposure.resize({mio::abm::CellIndex(location.get_cells().size()), mio::abm::VirusVariant::Count}); std::for_each(local_air_exposure.begin(), local_air_exposure.end(), [](auto& r) { r = 0.0; }); // allocate and initialize contact exposures with 0 mio::abm::ContactExposureRates local_contact_exposure; - local_contact_exposure.resize_destructive({mio::abm::CellIndex(location.get_cells().size()), - mio::abm::VirusVariant::Count, - mio::AgeGroup(global_parameters.get_num_groups())}); + local_contact_exposure.resize({mio::abm::CellIndex(location.get_cells().size()), mio::abm::VirusVariant::Count, + mio::AgeGroup(global_parameters.get_num_groups())}); std::for_each(local_contact_exposure.begin(), local_contact_exposure.end(), [](auto& r) { r = 0.0; });