From c111ba3666dec3ec963606b38d3788b334906494 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 29 Sep 2025 17:03:35 +0200 Subject: [PATCH 01/16] first draft influenza model --- cpp/CMakeLists.txt | 1 + cpp/examples/CMakeLists.txt | 4 + cpp/models/ode_seirv/CMakeLists.txt | 12 + cpp/models/ode_seirv/infection_state.h | 42 +++ cpp/models/ode_seirv/model.cpp | 20 ++ cpp/models/ode_seirv/model.h | 197 ++++++++++++++ cpp/models/ode_seirv/parameters.h | 358 +++++++++++++++++++++++++ 7 files changed, 634 insertions(+) create mode 100644 cpp/models/ode_seirv/CMakeLists.txt create mode 100644 cpp/models/ode_seirv/infection_state.h create mode 100644 cpp/models/ode_seirv/model.cpp create mode 100644 cpp/models/ode_seirv/model.h create mode 100644 cpp/models/ode_seirv/parameters.h diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 3ce345b3fd..6af0d5d74e 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -183,6 +183,7 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/ide_secir) add_subdirectory(models/ide_seir) add_subdirectory(models/ode_seir) + add_subdirectory(models/ode_seirv) add_subdirectory(models/ode_seair) add_subdirectory(models/ode_sir) add_subdirectory(models/sde_sir) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 39a12539e1..7b87f40bd2 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -68,6 +68,10 @@ add_executable(seir_flows_example ode_seir_flows.cpp) target_link_libraries(seir_flows_example PRIVATE memilio ode_seir) target_compile_options(seir_flows_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_seirv_example ode_seirv_example.cpp) +target_link_libraries(ode_seirv_example PRIVATE memilio ode_seirv) +target_compile_options(ode_seirv_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(ode_secir_example ode_secir.cpp) target_link_libraries(ode_secir_example PRIVATE memilio ode_secir) target_compile_options(ode_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_seirv/CMakeLists.txt b/cpp/models/ode_seirv/CMakeLists.txt new file mode 100644 index 0000000000..b6068e1b9d --- /dev/null +++ b/cpp/models/ode_seirv/CMakeLists.txt @@ -0,0 +1,12 @@ +add_library(ode_seirv + infection_state.h + model.h + model.cpp + parameters.h +) +target_link_libraries(ode_seirv PUBLIC memilio) +target_include_directories(ode_seirv PUBLIC + $ + $ +) +target_compile_options(ode_seirv PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_seirv/infection_state.h b/cpp/models/ode_seirv/infection_state.h new file mode 100644 index 0000000000..7eb0078eb9 --- /dev/null +++ b/cpp/models/ode_seirv/infection_state.h @@ -0,0 +1,42 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#ifndef SEIRV_INFECTIONSTATE_H +#define SEIRV_INFECTIONSTATE_H +namespace mio +{ +namespace oseirv +{ + +enum class InfectionState +{ + Susceptible, + SusceptibleVaccinated, + Exposed, + ExposedVaccinated, + Infected, + InfectedVaccinated, + Recovered, + RecoveredVaccinated, + Count +}; + +} // namespace oseirv +} // namespace mio +#endif diff --git a/cpp/models/ode_seirv/model.cpp b/cpp/models/ode_seirv/model.cpp new file mode 100644 index 0000000000..5c62bf43cf --- /dev/null +++ b/cpp/models/ode_seirv/model.cpp @@ -0,0 +1,20 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#include "ode_seirv/model.h" diff --git a/cpp/models/ode_seirv/model.h b/cpp/models/ode_seirv/model.h new file mode 100644 index 0000000000..f77fe64a95 --- /dev/null +++ b/cpp/models/ode_seirv/model.h @@ -0,0 +1,197 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#ifndef SEIRV_MODEL_H +#define SEIRV_MODEL_H + +#include "memilio/compartments/flow_model.h" +#include "memilio/config.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/epidemiology/populations.h" +#include "memilio/utils/time_series.h" +#include "ode_seirv/infection_state.h" +#include "ode_seirv/parameters.h" + +GCC_CLANG_DIAGNOSTIC(push) +GCC_CLANG_DIAGNOSTIC(ignored "-Wshadow") +#include +#include +GCC_CLANG_DIAGNOSTIC(pop) + +#include +#include + +namespace mio +{ +namespace oseirv +{ + +// clang-format off +using Flows = TypeList, + Flow, + Flow, + Flow, + Flow, + Flow>; +// clang-format on + +template +class Model + : public FlowModel, Parameters, Flows> +{ + using Base = FlowModel, Parameters, Flows>; + +public: + using ParameterSet = typename Base::ParameterSet; + using PopulationsType = typename Base::Populations; + + Model(const PopulationsType& pop, const ParameterSet& params) + : Base(pop, params) + { + } + Model(int num_agegroups) + : Base(PopulationsType({AgeGroup(num_agegroups), InfectionState::Count}), ParameterSet(AgeGroup(num_agegroups))) + { + } + + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override + { + const auto& P = this->parameters; + const Index AG = reduce_index>(this->populations.size()); + const size_t G = (size_t)P.get_num_groups(); + + // Contact matrices (may be time-dependent via damping): + const auto cm_h_expr = + P.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); + const auto cm_s_expr = + P.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); + + // Evaluate expressions to concrete matrices and extract mixing parameter as plain FP + Eigen::MatrixX H = cm_h_expr; + Eigen::MatrixX S = cm_s_expr; + const FP m = P.template get>(); + Eigen::MatrixX B = H + m * S; + + // Normalization ν(m) (spectral radius of Σ*B), Σ = φ * diag(σ_i) + Eigen::VectorX sigma(G); + for (size_t i = 0; i < G; ++i) + sigma[(Eigen::Index)i] = P.template get>()[AgeGroup(i)] * P.template get>(); + Eigen::MatrixX Sigma = sigma.asDiagonal(); + Eigen::MatrixX NG = Sigma * B; + + Eigen::ComplexEigenSolver> ces; + ces.compute(NG); + FP nu = ces.eigenvalues().cwiseAbs().maxCoeff(); + if (nu > FP(0)) { + B /= nu; + } + + const FP Re = P.template get>(); + const FP gamma = P.template get>(); + const FP delta = P.template get>(); + const FP tz = P.template get>(); + const FP ts = P.template get>(); + const FP lambda0 = P.template get>(); + const FP rho = P.template get>(); + + const FP season = std::exp(delta * std::sin(FP(2) * std::numbers::pi_v * (t / FP(52.0) - tz + ts))); + + for (auto i : make_index_range(AG)) { + // Flat indices + const size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible}); + const size_t SVi = this->populations.get_flat_index({i, InfectionState::SusceptibleVaccinated}); + const size_t Ei = this->populations.get_flat_index({i, InfectionState::Exposed}); + const size_t EVi = this->populations.get_flat_index({i, InfectionState::ExposedVaccinated}); + const size_t Ii = this->populations.get_flat_index({i, InfectionState::Infected}); + const size_t IVi = this->populations.get_flat_index({i, InfectionState::InfectedVaccinated}); + + // λ_i(t) = Re*γ*season * sum_j B_{ij} * ((I_j+I^V_j)/N_j)^ρ + λ0 + FP sum = FP(0); + for (auto j : make_index_range(AG)) { + const size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible}); + const size_t SVj = this->populations.get_flat_index({j, InfectionState::SusceptibleVaccinated}); + const size_t Ej = this->populations.get_flat_index({j, InfectionState::Exposed}); + const size_t EVj = this->populations.get_flat_index({j, InfectionState::ExposedVaccinated}); + const size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected}); + const size_t IVj = this->populations.get_flat_index({j, InfectionState::InfectedVaccinated}); + const size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered}); + const size_t RVj = this->populations.get_flat_index({j, InfectionState::RecoveredVaccinated}); + + const FP Nj = pop[Sj] + pop[SVj] + pop[Ej] + pop[EVj] + pop[Ij] + pop[IVj] + pop[Rj] + pop[RVj]; + const FP inf_frac = (Nj > FP(0)) ? ((pop[Ij] + pop[IVj]) / Nj) : FP(0); + const FP term = (rho == FP(1)) ? inf_frac : std::pow(inf_frac, rho); + sum += B(i.get(), j.get()) * term; + } + const FP lambda_i = Re * gamma * season * sum + lambda0; + + // Flows: S->E, S^V->E^V, E->I, E^V->I^V, I->R, I^V->R^V + flows[Base::template get_flat_flow_index(i)] = + lambda_i * y[Si]; + flows[Base::template get_flat_flow_index(i)] = lambda_i * y[SVi]; + + flows[Base::template get_flat_flow_index(i)] = + gamma * y[Ei]; + flows[Base::template get_flat_flow_index(i)] = gamma * y[EVi]; + + flows[Base::template get_flat_flow_index(i)] = + gamma * y[Ii]; + flows[Base::template get_flat_flow_index(i)] = gamma * y[IVi]; + } + } + + // Convenience: initialize according to paper: doi.org/10.1186/s12879-017-2344-6 + void set_initial_conditions_from_sigma_phi_vacc() + { + const auto& P = this->parameters; + const size_t G = (size_t)P.get_num_groups(); + for (size_t i = 0; i < G; ++i) { + using IndexT = typename PopulationsType::Index; + + // Read total N_i from the currently stored group total: + FP Ni = this->populations.get_group_total(gi); + + FP sigma = P.template get>()[gi]; + FP phi = P.template get>(); + FP VC = P.template get>()[gi]; + FP VE = P.template get>()[gi]; + + this->populations[IndexT{AgeGroup(i), InfectionState::Exposed}] = 0; + this->populations[IndexT{AgeGroup(i), InfectionState::ExposedVaccinated}] = 0; + this->populations[IndexT{AgeGroup(i), InfectionState::Infected}] = 0; + this->populations[IndexT{AgeGroup(i), InfectionState::InfectedVaccinated}] = 0; + + this->populations[IndexT{AgeGroup(i), InfectionState::Susceptible}] = phi * sigma * (FP(1) - VC) * Ni; + this->populations[IndexT{AgeGroup(i), InfectionState::SusceptibleVaccinated}] = + phi * sigma * (FP(1) - VE) * VC * Ni; + + this->populations[IndexT{AgeGroup(i), InfectionState::Recovered}] = + (FP(1) - phi * sigma) * (FP(1) - VC) * Ni; + this->populations[IndexT{AgeGroup(i), InfectionState::RecoveredVaccinated}] = + (FP(1) - phi * sigma * (FP(1) - VE)) * VC * Ni; + } + } +}; + +} // namespace oseirv +} // namespace mio +#endif diff --git a/cpp/models/ode_seirv/parameters.h b/cpp/models/ode_seirv/parameters.h new file mode 100644 index 0000000000..08771b3df5 --- /dev/null +++ b/cpp/models/ode_seirv/parameters.h @@ -0,0 +1,358 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#ifndef SEIRV_PARAMETERS_H +#define SEIRV_PARAMETERS_H + +#include "memilio/epidemiology/age_group.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/utils/parameter_set.h" + +namespace mio +{ +namespace oseirv +{ + +/** + * @brief Baseline transmissibility R_e (dimensionless). + * + * Multiplies the force of infection λ(t) after contact-matrix normalization. Main scalar for transmission + * intensity. Default: 1.0. + */ +template +struct BaselineTransmissibility { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(1.0); + } + static std::string name() + { + return "BaselineTransmissibility"; + } +}; + +/** + * @brief Transition/recovery rate γ in 1/week. + * + * Controls E→I and I→Rand is included in λ(t) to decouple from R_e. + * Use the model time unit (weeks). Default: 1/2.0 + */ +template +struct Gamma { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(1.0 / 2.0); + } + static std::string name() + { + return "Gamma"; + } +}; + +/** + * @brief Seasonality amplitude δ (dimensionless). + * + * Amplitude of the seasonal modulation exp(δ·sin(2π(t/52 − t_z + t_s))) in λ(t). + * Set to 0.0 to disable seasonality. Default: 0.0. + */ +template +struct SeasonalityAmplitude { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(0.0); + } + static std::string name() + { + return "SeasonalityAmplitude"; + } +}; + +/** + * @brief Subtype-specific seasonal shift t_z. + * + * Phase shift for the sinusoid in λ(t); may be adjusted by a seasonal-shift correction if implemented. Default: 0.0. + */ +template +struct ShiftTZ { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(0.0); + } + static std::string name() + { + return "ShiftTZ"; + } +}; + +/** + * @brief Season-specific fine shift t_s. + * + * Additional small seasonal timing correction per season. Default: 0.0. + */ +template +struct ShiftTS { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(0.0); + } + static std::string name() + { + return "ShiftTS"; + } +}; + +/** + * @brief Outside force of infection λ₀ in 1/week. + * + * Additive external hazard that can seed infections when no infectives are present. + * Default: 0.0. + */ +template +struct OutsideFoI { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(0.0); + } + static std::string name() + { + return "OutsideFoI"; + } +}; + +/** + * @brief Clustering/concavity parameter ρ (dimensionless). + * + * Exponent on the infectious fraction in λ(t). ρ=1 recovers mass-action; 0<ρ<1 yields concave, + * clustered transmission. Must be > 0. Default: 1.0. + */ +template +struct ClusteringRho { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(1.0); + } + static std::string name() + { + return "ClusteringRho"; + } +}; + +/** + * @brief Mixing parameter m for “sick” contacts (dimensionless). + * + * Weight for the symptomatically sick contact matrix in β_eff = ν(m)^{-1}(β_healthy + m·β_sick). + * Combines share of symptomatic cases and their higher infectiousness. Default: 1.0. + */ +template +struct SickMixingM { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(1.0); + } + static std::string name() + { + return "SickMixingM"; + } +}; + +/** + * @brief Contact patterns of healthy people (age-structured contact frequencies). + * + * Used in β_eff, can be time-varying via damping. Unit: contacts per model time step (week). + */ +template +struct ContactPatternsHealthy { + using Type = UncertainContactMatrix; + static Type get_default(AgeGroup size) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "ContactPatternsHealthy"; + } +}; + +/** + * @brief Contact patterns of symptomatically sick people (age-structured contact frequencies). + * + * Used in β_eff together with m, can be time-varying via damping. + * Unit: contacts per model time step (week). + */ +template +struct ContactPatternsSick { + using Type = UncertainContactMatrix; + static Type get_default(AgeGroup size) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "ContactPatternsSick"; + } +}; + +/** + * @brief Age-specific baseline susceptibility σ_i. + * + * Represents pre-existing subtype/age-specific immunity, used with Φ at t0 to form φ·σ_i. + * Default: 1.0 (fully susceptible). + */ +template +struct SigmaByAge { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(AgeGroup size) + { + return Type(size, 1.0); + } + static std::string name() + { + return "SigmaByAge"; + } +}; + +/** + * @brief Vaccination coverage VC_i at t0 (dimensionless, in [0,1]). + * + * Age-specific share vaccinated at season start; used in initial conditions. + * Default: 0.0. + */ +template +struct VaccineCoverage { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(AgeGroup size) + { + return Type(size, 0.0); + } + static std::string name() + { + return "VaccineCoverage"; + } +}; + +/** + * @brief Vaccine effectiveness VE_i (dimensionless, in [0,1]). + * + * Reduces effective susceptibility of vaccinated individuals in the initial state via (1 − VE_i). Default: 0.0. + */ +template +struct VaccineEffectiveness { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(AgeGroup size) + { + return Type(size, 0.0); + } + static std::string name() + { + return "VaccineEffectiveness"; + } +}; + +/** + * @brief Φ (phi): fraction of the population that remains susceptible at t0 (dimensionless, typically in [0,1]). + * + * Scales the susceptible pool at season start: φ·σ_i). Default: 1.0. + */ +template +struct Phi { + using Type = UncertainValue; + static Type get_default(AgeGroup) + { + return Type(1.0); + } + static std::string name() + { + return "Phi"; + } +}; + +template +using ParametersBase = + ParameterSet, Gamma, SeasonalityAmplitude, ShiftTZ, ShiftTS, + OutsideFoI, ClusteringRho, SickMixingM, ContactPatternsHealthy, + ContactPatternsSick, SigmaByAge, VaccineCoverage, VaccineEffectiveness, Phi>; + +/** + * @brief Parameter set for the age-resolved SEIRV model (S,E,I,R plus vaccinated states) as per the appendix. + * + * Contains seasonality, contact-pattern and immunity/vaccination parameters. + * Model time unit is week; contact matrices may be time-dependent (damping). + */ +template +class Parameters : public ParametersBase +{ +public: + /** + * @brief Construct with the number of age groups. + * @param ng Number of age groups. + */ + Parameters(AgeGroup ng) + : ParametersBase(ng) + , m_num_groups{ng} + { + } + + /** + * @brief Returns the number of age groups. + */ + AgeGroup get_num_groups() const + { + return m_num_groups; + } + + /** + * @brief Apply minimal constraints to keep parameters valid. + * + * @return true if one or more parameters were corrected; false otherwise. + */ + bool apply_constraints() + { + bool corrected = false; + if (this->template get>() <= 0.0) { + this->template get>() = 1.0; + corrected = true; + } + if (this->template get>() < 0.0) { + this->template get>() = 0.0; + corrected = true; + } + if (this->template get>() < 0.0) { + this->template get>() = 0.0; + corrected = true; + } + return corrected; + } + +private: + AgeGroup m_num_groups; +}; + +} // namespace oseirv +} // namespace mio +#endif From 1e8a70f9522d24d366dc96b9972d46544d8377c2 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 6 Oct 2025 12:25:14 +0200 Subject: [PATCH 02/16] add example and readme --- cpp/examples/ode_seirv_example.cpp | 117 +++++++++++++++++++++++++++++ cpp/models/ode_seirv/README.md | 17 +++++ cpp/models/ode_seirv/model.h | 38 +--------- 3 files changed, 137 insertions(+), 35 deletions(-) create mode 100644 cpp/examples/ode_seirv_example.cpp create mode 100644 cpp/models/ode_seirv/README.md diff --git a/cpp/examples/ode_seirv_example.cpp b/cpp/examples/ode_seirv_example.cpp new file mode 100644 index 0000000000..71ad2b5429 --- /dev/null +++ b/cpp/examples/ode_seirv_example.cpp @@ -0,0 +1,117 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "ode_seirv/model.h" +#include "memilio/compartments/simulation.h" +#include "memilio/utils/logging.h" + +/** + * @brief set_initial_population sets the initial population of the model + * + * We assume that no one is initially infected or exposed at the beginning of the season. + * Infections are seeded later via the OutsideFoI parameter. + * The destinction into the 2 layers of susceptibility is done via the parameters. + */ +template +void set_initial_population(mio::oseirv::Model& model, const size_t total_pop) +{ + auto& params = model.parameters; + auto& pop = model.populations; + const size_t num_groups = (size_t)params.get_num_groups(); + + for (size_t i = 0; i < num_groups; ++i) { + pop.template set_difference_from_group_total( + {mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}, total_pop / num_groups); + + // Total population N_i as currently stored in group totals + FP Ni = pop.get_group_total(mio::AgeGroup(i)); + + FP sigma = params.template get>()[mio::AgeGroup(i)]; + FP phi = params.template get>(); + FP VC = params.template get>()[mio::AgeGroup(i)]; + FP VE = params.template get>()[mio::AgeGroup(i)]; + + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Exposed}] = 0; + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::ExposedVaccinated}] = 0; + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Infected}] = 0; + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::InfectedVaccinated}] = 0; + + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}] = phi * sigma * (FP(1) - VC) * Ni; + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::SusceptibleVaccinated}] = + phi * sigma * (FP(1) - VE) * VC * Ni; + + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Recovered}] = (FP(1) - phi * sigma) * (FP(1) - VC) * Ni; + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::RecoveredVaccinated}] = + (FP(1) - phi * sigma * (FP(1) - VE)) * VC * Ni; + } +} + +// Example usage of the SEIRV ODE model presented in https://doi.org/10.1186/s12879-017-2344-6 +int main() +{ + using FP = double; + mio::set_log_level(mio::LogLevel::debug); + + const FP t0 = 0., tmax = 42.; + const FP dt = 0.1; + + const size_t num_groups = 1; + const auto total_pop = 1e5; + mio::oseirv::Model model((int)num_groups); + auto& parameters = model.parameters; + + FP cont_freq = 10.0; + FP fact = FP(1) / FP(num_groups); + + // Healthy contacts + mio::ContactMatrixGroup& contacts_healthy = parameters.get>(); + contacts_healthy[0] = mio::ContactMatrix( + Eigen::MatrixXd::Constant((Eigen::Index)num_groups, (Eigen::Index)num_groups, fact * cont_freq)); + + // Sick contacts (here 20% reduced) + mio::ContactMatrixGroup& contacts_sick = parameters.get>(); + contacts_sick[0] = mio::ContactMatrix( + Eigen::MatrixXd::Constant((Eigen::Index)num_groups, (Eigen::Index)num_groups, fact * cont_freq * 0.8)); + + // Parameters + parameters.get>() = 1.2; + parameters.get>() = 1.0 / 2.0; + parameters.get>() = 1.0; + parameters.get>() = 0.0; + parameters.get>() = 0.0; + parameters.get>() = 1e-6; + parameters.get>() = 0.9; + parameters.get>() = 2.0; + + for (size_t i = 0; i < num_groups; ++i) { + parameters.get>()[mio::AgeGroup(i)] = 1.0; + parameters.get>()[mio::AgeGroup(i)] = 0.3; + parameters.get>()[mio::AgeGroup(i)] = 0.5; + model.populations.set_difference_from_group_total( + {mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}, 1e5 / num_groups); + } + + set_initial_population(model, total_pop); + + auto seirv = mio::simulate(t0, tmax, dt, model); + + std::vector vars = {"S", "SV", "E", "EV", "I", "IV", "R", "RV"}; + seirv.print_table(vars); +} diff --git a/cpp/models/ode_seirv/README.md b/cpp/models/ode_seirv/README.md new file mode 100644 index 0000000000..ced80a5a1a --- /dev/null +++ b/cpp/models/ode_seirv/README.md @@ -0,0 +1,17 @@ +# ODE SEIRV (seasonal influenza) + +A age‑structured **SEIRV** model implemented into the MEmilio ODE framework. +It extends the classic SEIR by a vaccinated layer (Sᵛ, Eᵛ, Iᵛ, Rᵛ) and follows the seasonal influenza +model structure from: + +> **Weidemann, F., Remschmidt, C., Buda, S. et al.** +> *Is the impact of childhood influenza vaccination less than expected: a transmission modelling study.* +> **BMC Infectious Diseases** 17, 258 (2017). +> https://doi.org/10.1186/s12879-017-2344-6 + +## Reference + +Please cite the original study if you use this model: + +> Weidemann, F., Remschmidt, C., Buda, S. et al. *Is the impact of childhood influenza vaccination less than expected: a transmission modelling study.* **BMC Infect Dis** 17, 258 (2017). +> https://doi.org/10.1186/s12879-017-2344-6 diff --git a/cpp/models/ode_seirv/model.h b/cpp/models/ode_seirv/model.h index f77fe64a95..1430664bb6 100644 --- a/cpp/models/ode_seirv/model.h +++ b/cpp/models/ode_seirv/model.h @@ -75,7 +75,7 @@ class Model { const auto& P = this->parameters; const Index AG = reduce_index>(this->populations.size()); - const size_t G = (size_t)P.get_num_groups(); + const size_t num_groups = (size_t)P.get_num_groups(); // Contact matrices (may be time-dependent via damping): const auto cm_h_expr = @@ -90,8 +90,8 @@ class Model Eigen::MatrixX B = H + m * S; // Normalization ν(m) (spectral radius of Σ*B), Σ = φ * diag(σ_i) - Eigen::VectorX sigma(G); - for (size_t i = 0; i < G; ++i) + Eigen::VectorX sigma(num_groups); + for (size_t i = 0; i < num_groups; ++i) sigma[(Eigen::Index)i] = P.template get>()[AgeGroup(i)] * P.template get>(); Eigen::MatrixX Sigma = sigma.asDiagonal(); Eigen::MatrixX NG = Sigma * B; @@ -158,38 +158,6 @@ class Model InfectionState::RecoveredVaccinated>(i)] = gamma * y[IVi]; } } - - // Convenience: initialize according to paper: doi.org/10.1186/s12879-017-2344-6 - void set_initial_conditions_from_sigma_phi_vacc() - { - const auto& P = this->parameters; - const size_t G = (size_t)P.get_num_groups(); - for (size_t i = 0; i < G; ++i) { - using IndexT = typename PopulationsType::Index; - - // Read total N_i from the currently stored group total: - FP Ni = this->populations.get_group_total(gi); - - FP sigma = P.template get>()[gi]; - FP phi = P.template get>(); - FP VC = P.template get>()[gi]; - FP VE = P.template get>()[gi]; - - this->populations[IndexT{AgeGroup(i), InfectionState::Exposed}] = 0; - this->populations[IndexT{AgeGroup(i), InfectionState::ExposedVaccinated}] = 0; - this->populations[IndexT{AgeGroup(i), InfectionState::Infected}] = 0; - this->populations[IndexT{AgeGroup(i), InfectionState::InfectedVaccinated}] = 0; - - this->populations[IndexT{AgeGroup(i), InfectionState::Susceptible}] = phi * sigma * (FP(1) - VC) * Ni; - this->populations[IndexT{AgeGroup(i), InfectionState::SusceptibleVaccinated}] = - phi * sigma * (FP(1) - VE) * VC * Ni; - - this->populations[IndexT{AgeGroup(i), InfectionState::Recovered}] = - (FP(1) - phi * sigma) * (FP(1) - VC) * Ni; - this->populations[IndexT{AgeGroup(i), InfectionState::RecoveredVaccinated}] = - (FP(1) - phi * sigma * (FP(1) - VE)) * VC * Ni; - } - } }; } // namespace oseirv From 35262a67d88ce53ab23626942523e91402e59eef Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 6 Oct 2025 14:55:48 +0200 Subject: [PATCH 03/16] adjust params and cast in example --- cpp/examples/ode_seirv_example.cpp | 2 +- cpp/models/ode_seirv/parameters.h | 24 +++++++++++------------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/cpp/examples/ode_seirv_example.cpp b/cpp/examples/ode_seirv_example.cpp index 71ad2b5429..eec8b0c13c 100644 --- a/cpp/examples/ode_seirv_example.cpp +++ b/cpp/examples/ode_seirv_example.cpp @@ -30,7 +30,7 @@ * The destinction into the 2 layers of susceptibility is done via the parameters. */ template -void set_initial_population(mio::oseirv::Model& model, const size_t total_pop) +void set_initial_population(mio::oseirv::Model& model, const FP total_pop) { auto& params = model.parameters; auto& pop = model.populations; diff --git a/cpp/models/ode_seirv/parameters.h b/cpp/models/ode_seirv/parameters.h index 08771b3df5..26aecc0178 100644 --- a/cpp/models/ode_seirv/parameters.h +++ b/cpp/models/ode_seirv/parameters.h @@ -34,7 +34,7 @@ namespace oseirv /** * @brief Baseline transmissibility R_e (dimensionless). * - * Multiplies the force of infection λ(t) after contact-matrix normalization. Main scalar for transmission + * Multiplies the force of infection after contact-matrix normalization. Main scalar for transmission * intensity. Default: 1.0. */ template @@ -51,9 +51,9 @@ struct BaselineTransmissibility { }; /** - * @brief Transition/recovery rate γ in 1/week. + * @brief Transition/recovery rate in 1/week. * - * Controls E→I and I→Rand is included in λ(t) to decouple from R_e. + * Controls E→I and I→Rand is included in the force of infection to decouple from R_e. * Use the model time unit (weeks). Default: 1/2.0 */ template @@ -70,9 +70,9 @@ struct Gamma { }; /** - * @brief Seasonality amplitude δ (dimensionless). + * @brief Seasonality amplitude (dimensionless). * - * Amplitude of the seasonal modulation exp(δ·sin(2π(t/52 − t_z + t_s))) in λ(t). + * Amplitude of the seasonal modulation. * Set to 0.0 to disable seasonality. Default: 0.0. */ template @@ -91,7 +91,7 @@ struct SeasonalityAmplitude { /** * @brief Subtype-specific seasonal shift t_z. * - * Phase shift for the sinusoid in λ(t); may be adjusted by a seasonal-shift correction if implemented. Default: 0.0. + * Phase shift for the sinusoid in the force of incetion.; may be adjusted by a seasonal-shift correction if implemented. Default: 0.0. */ template struct ShiftTZ { @@ -125,9 +125,10 @@ struct ShiftTS { }; /** - * @brief Outside force of infection λ₀ in 1/week. + * @brief Outside force of infection in 1/week. * * Additive external hazard that can seed infections when no infectives are present. + * Used for the beginning of the season (to import infections) or to model travel/imported cases. * Default: 0.0. */ template @@ -146,8 +147,7 @@ struct OutsideFoI { /** * @brief Clustering/concavity parameter ρ (dimensionless). * - * Exponent on the infectious fraction in λ(t). ρ=1 recovers mass-action; 0<ρ<1 yields concave, - * clustered transmission. Must be > 0. Default: 1.0. + * Exponent on the infectious fraction in force of infection. Must be > 0. Default: 1.0. */ template struct ClusteringRho { @@ -221,7 +221,7 @@ struct ContactPatternsSick { /** * @brief Age-specific baseline susceptibility σ_i. * - * Represents pre-existing subtype/age-specific immunity, used with Φ at t0 to form φ·σ_i. + * Represents pre-existing subtype/age-specific immunity. * Default: 1.0 (fully susceptible). */ template @@ -275,9 +275,7 @@ struct VaccineEffectiveness { }; /** - * @brief Φ (phi): fraction of the population that remains susceptible at t0 (dimensionless, typically in [0,1]). - * - * Scales the susceptible pool at season start: φ·σ_i). Default: 1.0. + * @brief phi: fraction of the population that remains susceptible at t0 (dimensionless, typically in [0,1]). */ template struct Phi { From 7019e1d091f4cf0b47abda6f6eb801f59261217e Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 6 Oct 2025 15:02:50 +0200 Subject: [PATCH 04/16] [ci skip] rtd doc --- docs/source/cpp/models/oseirv.rst | 242 ++++++++++++++++++++++++++++++ 1 file changed, 242 insertions(+) create mode 100644 docs/source/cpp/models/oseirv.rst diff --git a/docs/source/cpp/models/oseirv.rst b/docs/source/cpp/models/oseirv.rst new file mode 100644 index 0000000000..539f05d676 --- /dev/null +++ b/docs/source/cpp/models/oseirv.rst @@ -0,0 +1,242 @@ +ODE-based SEIRV model +====================== + +The ODE-SEIRV module extends the classic SEIR model by an explicit vaccinated layer (``S^V, E^V, I^V, R^V``) and targets +seasonal influenza–type applications (cf. Weidemann et al. 2017). It is age-structured and uses a normalized contact +matrix formulation with explicit seasonality. Immunity after recovery is assumed to last for the considered season. Also, there is no transition between the unvaccinated and vaccinated compartments during the simulation horizon (i.e., no vaccinations after season start). + +The infection states and transitions are illustrated in the following figure. + +.. image:: https://martinkuehn.eu/research/images/seirv.png + :alt: SEIRV_model + + +Infection States +---------------- + +The model contains the following **InfectionState**\s: + +.. code-block:: RST + + `Susceptible` + `SusceptibleVaccinated` + `Exposed` + `ExposedVaccinated` + `Infected` + `InfectedVaccinated` + `Recovered` + `RecoveredVaccinated` + + +Infection State Transitions +--------------------------- + +The SEIRV model is implemented as a **FlowModel**. Thus, in each time step, the flows (new infections, progressions, +recoveries) are computed explicitly in addition to compartment values. The defined transitions `FromState, ToState` are: + +.. code-block:: RST + + `Susceptible, Exposed` + `SusceptibleVaccinated, ExposedVaccinated` + `Exposed, Infected` + `ExposedVaccinated, InfectedVaccinated` + `Infected, Recovered` + `InfectedVaccinated, RecoveredVaccinated` + + +Sociodemographic Stratification +-------------------------------- + +The population can be stratified by one sociodemographic dimension denoted **AgeGroup** (can be interpreted more +broadly). The number of age groups is specified in the constructor: + +.. code-block:: cpp + + mio::oseirv::Model model(num_agegroups); + +For stratifications with two or more dimensions, see :doc:`Model Creation <../ode_creation>`. + + +Parameters +---------- + +The model uses the following parameters (time unit: week): + +.. list-table:: + :header-rows: 1 + :widths: 20 25 55 + + * - Mathematical Symbol + - C++ Name / Type + - Description + * - :math:`R_e` + - ``BaselineTransmissibility`` + - Baseline transmissibility (dimensionless); scales the normalized force of infection. + * - :math:`\gamma` + - ``Gamma`` + - Transition / recovery rate (1/week) for E -> I and I -> R. + * - :math:`\delta` + - ``SeasonalityAmplitude`` + - Amplitude of the seasonal modulation :math:`\exp(\delta\,\sin(2\pi(t/52 - t_z + t_s)))`. + * - :math:`t_z` + - ``ShiftTZ`` + - Coarse (subtype-specific) seasonal phase shift. + * - :math:`t_s` + - ``ShiftTS`` + - Fine seasonal phase adjustment per season. + * - :math:`\lambda_0` + - ``OutsideFoI`` + - External (additive) force of infection, can seed infections. + * - :math:`\rho` + - ``ClusteringRho`` + - Clustering exponent on the infectious fraction. + * - :math:`m` + - ``SickMixingM`` + - Mixing weight for symptomatic (“sick”) contacts in the blended contact matrix. + * - :math:`C^{H}` + - ``ContactPatternsHealthy`` + - Age-structured contact matrix (healthy). Can be time-dependent via damping. + * - :math:`C^{S}` + - ``ContactPatternsSick`` + - Age-structured contact matrix (symptomatic), combined using :math:`m`. + * - :math:`\sigma_i` + - ``SigmaByAge`` + - Age-specific baseline susceptibility (pre-existing immunity modifier). + * - :math:`VC_i` + - ``VaccineCoverage`` + - Vaccination coverage per age group at season start (share vaccinated). + * - :math:`VE_i` + - ``VaccineEffectiveness`` + - Vaccine effectiveness (reducing effective susceptibility). + * - :math:`\phi_0` + - ``Phi`` + - Fraction of the total population forming the effectively susceptible pool at :math:`t_0`. + +Note: ``VaccineCoverage`` and ``VaccineEffectiveness`` are only used for initialization. Transitions presently +apply identical hazards to vaccinated and unvaccinated susceptible compartments. Future extensions may introduce +differential infection hazards. + + +Initial Conditions +------------------ + +Initial conditions are handled via the **Populations** class. Example for a single age group: + +.. code-block:: cpp + + mio::oseirv::Model model(1); + // Set total population in age group 0 + model.populations.set_total(total0); + + // Initialize vaccinated susceptibles (simple example) + double vc0 = 0.4; // vaccination coverage + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::SusceptibleVaccinated}] = vc0 * total0; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = initial_infected; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Exposed}] = initial_exposed; + + // Other states (Recovered / RecoveredVaccinated) often 0 at season start + + // Set remaining susceptibles as difference + model.populations.set_difference_from_total( + {mio::AgeGroup(0), mio::oseirv::InfectionState::Susceptible}, total0); + +For age-resolved simulations, repeat for each age group; ``set_difference_from_group_total`` ensures correct +susceptible counts: + +.. code-block:: cpp + + for (auto a = mio::AgeGroup(0); a < num_agegroups; ++a) { + model.populations[{a, mio::oseirv::InfectionState::Exposed}] = exposed0 / num_agegroups; + model.populations[{a, mio::oseirv::InfectionState::Infected}] = infected0 / num_agegroups; + model.populations[{a, mio::oseirv::InfectionState::SusceptibleVaccinated}] = vc[a.get()] * group_size[a.get()]; + model.populations.set_difference_from_group_total( + {a, mio::oseirv::InfectionState::Susceptible}, group_size[a.get()]); + } + + +Simulation +---------- + +Like other ODE models in MEmilio, the SEIRV model can be simulated with standard compartment output or with explicit +flows. Once integrated with utility wrappers (analogous to ``oseir::simulate`` / ``simulate_flows``) usage follows the +same pattern. Example with a Runge–Kutta integrator: + +.. code-block:: cpp + + double t0 = 0.0; // start (weeks) + double tmax = 20.0; // end + double dt = 0.1; // initial step size + + auto integrator = std::make_unique(); + integrator->set_dt_min(0.01); + integrator->set_dt_max(0.5); + integrator->set_rel_tolerance(1e-4); + integrator->set_abs_tolerance(1e-6); + + auto sim = mio::simulate(t0, tmax, dt, model, std::move(integrator)); + +Flow simulation (when explicit flows are required): + +.. code-block:: cpp + + auto flowsim = mio::simulate_flows(t0, tmax, dt, model); + // flowsim[0] = compartment sizes, flowsim[1] = flows + + +Output +------ + +The result of a standard simulation is a ``mio::TimeSeries``: + +.. code-block:: cpp + + auto n_points = static_cast(sim.get_num_time_points()); + Eigen::VectorXd val_i = sim.get_value(i); + double time_i = sim.get_time(i); + auto last_val = sim.get_last_value(); + +Printing and CSV export: + +.. code-block:: cpp + + sim.print_table(); + std::vector labels = {"S","S_V","E","E_V","I","I_V","R","R_V"}; + sim.print_table(labels); + sim.export_csv("seirv_results.csv"); + + + +Contact Changes / Interventions +-------------------------------- + +Time-dependent changes of contact patterns (holidays, interventions) can be modeled via dampings (``add_damping``) on +``ContactPatternsHealthy`` and/or ``ContactPatternsSick``: + +.. code-block:: cpp + + mio::ContactMatrixGroup& cm_h = model.parameters.get>(); + mio::ContactMatrixGroup& cm_s = model.parameters.get>(); + cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, base_contacts)); + cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, base_contacts_sick)); + + // Reduce healthy contacts by 40% starting at week 5 + cm_h[0].add_damping(0.6, mio::SimulationTime(5.0)); + + +Visualization +------------- + +For visualization you can use the Python package :doc:`m-plot <../../python/m-plot>` as in the other models. + + +Literature +---------- + +* Weidemann, F., Remschmidt, C., Buda, S. et al. *Is the impact of childhood influenza vaccination less than expected: a transmission modelling study.* BMC Infectious Diseases 17, 258 (2017). https://doi.org/10.1186/s12879-017-2344-6 + + +Overview of the ``oseirv`` namespace: +------------------------------------- + +.. doxygennamespace:: mio::oseirv + From c3f4136a2fcdd85de1ca611843e4c4bb27d34712 Mon Sep 17 00:00:00 2001 From: HenrZu <69154294+HenrZu@users.noreply.github.com> Date: Mon, 6 Oct 2025 21:51:44 +0200 Subject: [PATCH 05/16] add tests and rework model.h --- cpp/models/ode_seirv/model.h | 50 +++--- cpp/tests/CMakeLists.txt | 1 + cpp/tests/test_odeseirv.cpp | 284 +++++++++++++++++++++++++++++++++++ 3 files changed, 314 insertions(+), 21 deletions(-) create mode 100644 cpp/tests/test_odeseirv.cpp diff --git a/cpp/models/ode_seirv/model.h b/cpp/models/ode_seirv/model.h index 1430664bb6..d047e4e96f 100644 --- a/cpp/models/ode_seirv/model.h +++ b/cpp/models/ode_seirv/model.h @@ -73,47 +73,55 @@ class Model void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, Eigen::Ref> flows) const override { - const auto& P = this->parameters; - const Index AG = reduce_index>(this->populations.size()); - const size_t num_groups = (size_t)P.get_num_groups(); + const auto& params = this->parameters; + const size_t num_groups = (size_t)params.get_num_groups(); - // Contact matrices (may be time-dependent via damping): + // Contact matrices const auto cm_h_expr = - P.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); + params.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); const auto cm_s_expr = - P.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); + params.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); - // Evaluate expressions to concrete matrices and extract mixing parameter as plain FP + // Get effective contact matrix B = H + m*S Eigen::MatrixX H = cm_h_expr; Eigen::MatrixX S = cm_s_expr; - const FP m = P.template get>(); + const FP m = params.template get>(); Eigen::MatrixX B = H + m * S; - // Normalization ν(m) (spectral radius of Σ*B), Σ = φ * diag(σ_i) + // Normalization ν(m): prepare susceptibility scaling and next-generation matrix (NG) + // sigma_i = Phi * SigmaByAge[i] (age-specific susceptibility scaled by a global factor) Eigen::VectorX sigma(num_groups); for (size_t i = 0; i < num_groups; ++i) - sigma[(Eigen::Index)i] = P.template get>()[AgeGroup(i)] * P.template get>(); + sigma[(Eigen::Index)i] = + params.template get>()[AgeGroup(i)] * params.template get>(); + // Sigma is a diagonal matrix applying susceptibility on the receiving (row) side of contacts. Eigen::MatrixX Sigma = sigma.asDiagonal(); - Eigen::MatrixX NG = Sigma * B; + + // NG = Sigma * B is the (unnormalized) next-generation matrix. Its spectral radius will be used + // below to scale contacts B so that structure (mixing pattern) and transmissibility magnitude (BaselineTransmissibility) + // are cleanly separated. + Eigen::MatrixX NG = Sigma * B; Eigen::ComplexEigenSolver> ces; ces.compute(NG); FP nu = ces.eigenvalues().cwiseAbs().maxCoeff(); - if (nu > FP(0)) { + if (nu > FP(0.)) { B /= nu; } - const FP Re = P.template get>(); - const FP gamma = P.template get>(); - const FP delta = P.template get>(); - const FP tz = P.template get>(); - const FP ts = P.template get>(); - const FP lambda0 = P.template get>(); - const FP rho = P.template get>(); + const FP Re = + params.template get>(); // baseline transmissibility scaling (R-like factor) + const FP gamma = params.template get>(); // progression/recovery rate (1 / mean duration) + const FP delta = params.template get>(); // amplitude of seasonal modulation + const FP tz = params.template get>(); // base phase shift for seasonality + const FP ts = params.template get>(); // additional phase shift (e.g. scenario-specific) + const FP lambda0 = params.template get>(); // constant external force of infection + const FP rho = params.template get< + ClusteringRho>(); // clustering exponent; rho<1 dampens, rho>1 amplifies prevalence effect const FP season = std::exp(delta * std::sin(FP(2) * std::numbers::pi_v * (t / FP(52.0) - tz + ts))); - for (auto i : make_index_range(AG)) { + for (auto i : make_index_range((mio::AgeGroup)num_groups)) { // Flat indices const size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible}); const size_t SVi = this->populations.get_flat_index({i, InfectionState::SusceptibleVaccinated}); @@ -124,7 +132,7 @@ class Model // λ_i(t) = Re*γ*season * sum_j B_{ij} * ((I_j+I^V_j)/N_j)^ρ + λ0 FP sum = FP(0); - for (auto j : make_index_range(AG)) { + for (auto j : make_index_range((mio::AgeGroup)num_groups)) { const size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible}); const size_t SVj = this->populations.get_flat_index({j, InfectionState::SusceptibleVaccinated}); const size_t Ej = this->populations.get_flat_index({j, InfectionState::Exposed}); diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 01db381a7e..f2933e0855 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -8,6 +8,7 @@ set(TESTSOURCES testmain.cpp test_populations.cpp test_odeseir.cpp + test_odeseirv.cpp test_odesir.cpp test_numericalIntegration.cpp test_smoother.cpp diff --git a/cpp/tests/test_odeseirv.cpp b/cpp/tests/test_odeseirv.cpp new file mode 100644 index 0000000000..38985ec5e4 --- /dev/null +++ b/cpp/tests/test_odeseirv.cpp @@ -0,0 +1,284 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "memilio/config.h" +#include "memilio/math/integrator.h" +#include "memilio/math/euler.h" +#include "memilio/utils/time_series.h" +#include "memilio/compartments/simulation.h" +#include "memilio/compartments/flow_simulation.h" + +#include "ode_seirv/model.h" +#include "ode_seirv/infection_state.h" +#include "ode_seirv/parameters.h" + +#include + +#include +#include + +TEST(TestOdeSeirv, simulateDefault) +{ + double t0 = 0.0; + double tmax = 1.0; + double dt = 0.1; + + mio::oseirv::Model model(1); + auto result = simulate(t0, tmax, dt, model); // generic template simulate + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-12); + EXPECT_EQ(result.get_num_elements(), (Eigen::Index)mio::oseirv::InfectionState::Count); +} + +TEST(TestOdeSeirv, populationConservation) +{ + mio::oseirv::Model model(1); + const double S = 300; + const double SV = 200; + const double E = 50; + const double EV = 25; + const double I = 40; + const double IV = 10; + const double R = 0; + const double RV = 0; + double total = S + SV + E + EV + I + IV + R + RV; + + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Susceptible}] = S; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::SusceptibleVaccinated}] = SV; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Exposed}] = E; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::ExposedVaccinated}] = EV; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = I; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::InfectedVaccinated}] = IV; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Recovered}] = R; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::RecoveredVaccinated}] = RV; + + mio::ContactMatrixGroup& cm_h = model.parameters.get>(); + cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); + mio::ContactMatrixGroup& cm_s = model.parameters.get>(); + cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); + + model.parameters.set>(1.0); + model.parameters.set>(1.0); + + double t0 = 0.0, tmax = 5.0, dt = 0.5; + auto sim = simulate(t0, tmax, dt, model); + auto last = sim.get_last_value(); + EXPECT_NEAR(last.sum(), total, 1e-8); +} + +TEST(TestOdeSeirv, applyConstraints) +{ + mio::oseirv::Model model(1); + // First: defaults should need no correction + EXPECT_FALSE(model.parameters.apply_constraints()); + + model.parameters.set>(0.0); // invalid, must become >0 (1.0) + model.parameters.set>(-2.); // invalid -> 0 + model.parameters.set>(-0.5); // invalid -> 0 + EXPECT_TRUE(model.parameters.apply_constraints()); + EXPECT_GT(model.parameters.get>(), 0.0); + EXPECT_EQ(model.parameters.get>(), 0.0); + EXPECT_EQ(model.parameters.get>(), 0.0); +} + +TEST(TestOdeSeirv, flowsSingleAgeGroup) +{ + mio::oseirv::Model model(1); + + // Populations + const double S = 300, SV = 200, E = 50, EV = 25, I = 40, IV = 10; // R, RV = 0 + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Susceptible}] = S; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::SusceptibleVaccinated}] = SV; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Exposed}] = E; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::ExposedVaccinated}] = EV; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = I; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::InfectedVaccinated}] = IV; + + // Parameters: identity contacts, gamma=1, Re=1, others zero ⇒ λ = (I+IV)/N + mio::ContactMatrixGroup& cm_h = model.parameters.get>(); + cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); + mio::ContactMatrixGroup& cm_s = model.parameters.get>(); + cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); + model.parameters.set>(1.0); + model.parameters.set>(1.0); + model.parameters.set>(1.0); + model.parameters.set>(0.0); + model.parameters.set>(0.0); + + auto y0 = model.get_initial_values(); + auto pop = y0; // same vector for signature + Eigen::VectorXd flows(6); // number of flows for one age group + flows.setZero(); + model.get_flows(pop, y0, 0.0, flows); + + const double N = S + SV + E + EV + I + IV; // (R,RV = 0) + const double lambda = (I + IV) / N; // expected force of infection + const double f_SE = S * lambda; + const double f_SV_EV = SV * lambda; + const double f_E_I = E; // gamma * E + const double f_EV_IV = EV; // gamma * EV + const double f_I_R = I; // gamma * I + const double f_IV_RV = IV; // gamma * IV + + auto idx_SE = model.template get_flat_flow_index(mio::AgeGroup(0)); + auto idx_SV = model.template get_flat_flow_index(mio::AgeGroup(0)); + auto idx_EI = + model.template get_flat_flow_index( + mio::AgeGroup(0)); + auto idx_EVI = + model.template get_flat_flow_index(mio::AgeGroup(0)); + auto idx_IR = model.template get_flat_flow_index(mio::AgeGroup(0)); + auto idx_IVR = + model.template get_flat_flow_index(mio::AgeGroup(0)); + + EXPECT_NEAR(flows[idx_SE], f_SE, 1e-12); + EXPECT_NEAR(flows[idx_SV], f_SV_EV, 1e-12); + EXPECT_NEAR(flows[idx_EI], f_E_I, 1e-12); + EXPECT_NEAR(flows[idx_EVI], f_EV_IV, 1e-12); + EXPECT_NEAR(flows[idx_IR], f_I_R, 1e-12); + EXPECT_NEAR(flows[idx_IVR], f_IV_RV, 1e-12); +} + +TEST(TestOdeSeirv, flowsTwoAgeGroupsIdentityContacts) +{ + mio::oseirv::Model model(2); + mio::ContactMatrixGroup& cm_h = model.parameters.get>(); + // Each group only contacts itself so that λ_i = I_i / N_i + Eigen::MatrixXd Id2 = Eigen::MatrixXd::Identity(2, 2); + cm_h[0] = mio::ContactMatrix(Id2); + mio::ContactMatrixGroup& cm_s = model.parameters.get>(); + cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Zero(2, 2)); + model.parameters.set>(1.0); + model.parameters.set>(1.0); + model.parameters.set>(1.0); + model.parameters.set>(0.0); + model.parameters.set>(0.0); + + // Group 0 + double S0 = 100, I0 = 20; // others zero + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Susceptible}] = S0; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = I0; + // Group 1 + double S1 = 80, I1 = 10; + model.populations[{mio::AgeGroup(1), mio::oseirv::InfectionState::Susceptible}] = S1; + model.populations[{mio::AgeGroup(1), mio::oseirv::InfectionState::Infected}] = I1; + + auto y0 = model.get_initial_values(); + auto pop = y0; + Eigen::VectorXd flows(12); // 6 flows * 2 age groups + flows.setZero(); + model.get_flows(pop, y0, 0.0, flows); + + double N0 = S0 + I0; + double N1 = S1 + I1; + double lambda0 = I0 / N0; // identity contacts => only own group contributes + double lambda1 = I1 / N1; + + auto idx_SE_0 = model.template get_flat_flow_index(mio::AgeGroup(0)); + auto idx_SE_1 = model.template get_flat_flow_index(mio::AgeGroup(1)); + auto idx_EI_0 = + model.template get_flat_flow_index( + mio::AgeGroup(0)); + auto idx_EI_1 = + model.template get_flat_flow_index( + mio::AgeGroup(1)); + auto idx_IR_0 = model.template get_flat_flow_index(mio::AgeGroup(0)); + auto idx_IR_1 = model.template get_flat_flow_index(mio::AgeGroup(1)); + + EXPECT_NEAR(flows[idx_SE_0], S0 * lambda0, 1e-12); + EXPECT_NEAR(flows[idx_SE_1], S1 * lambda1, 1e-12); + // Exposed are zero => progressions must be zero + EXPECT_NEAR(flows[idx_EI_0], 0.0, 1e-12); + EXPECT_NEAR(flows[idx_EI_1], 0.0, 1e-12); + EXPECT_NEAR(flows[idx_IR_0], I0, 1e-12); + EXPECT_NEAR(flows[idx_IR_1], I1, 1e-12); +} + +TEST(TestOdeSeirv, zeroPopulationNoNan) +{ + mio::oseirv::Model model(1); + model.populations.set_total(0.0); + auto y0 = model.get_initial_values(); + auto pop = y0; + Eigen::VectorXd flows(6); + flows.setZero(); + model.get_flows(pop, y0, 0.0, flows); + for (int i = 0; i < flows.size(); ++i) { + EXPECT_FALSE(std::isnan(flows[i])); + EXPECT_EQ(flows[i], 0.0); + } +} + +TEST(TestOdeSeirv, simulationEuler) +{ + mio::oseirv::Model model(1); + // Simple initial values + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Exposed}] = 10; + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = 5; + model.populations.set_difference_from_group_total( + {mio::AgeGroup(0), mio::oseirv::InfectionState::Susceptible}, 100.0); + + // Identity contacts + mio::ContactMatrixGroup& cm_h = model.parameters.get>(); + cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); + mio::ContactMatrixGroup& cm_s = model.parameters.get>(); + cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 0.)); + model.parameters.set>(1.0); + model.parameters.set>(1.0); + + double t0 = 0.0, tmax = 1.0, dt = 0.5; + std::unique_ptr> integrator = std::make_unique>(); + auto sim = simulate(t0, tmax, dt, model, std::move(integrator)); + EXPECT_EQ(sim.get_num_time_points(), 3); // t=0,0.5,1.0 + // Sanity: all values non-negative + for (Eigen::Index i = 0; i < sim.get_last_value().size(); ++i) { + EXPECT_GE(sim.get_last_value()[i], 0.0); + } +} + +TEST(TestOdeSeirv, flowSimulation) +{ + mio::oseirv::Model model(1); + model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = 5; + model.populations.set_difference_from_group_total( + {mio::AgeGroup(0), mio::oseirv::InfectionState::Susceptible}, 100.0); + + mio::ContactMatrixGroup& cm_h = model.parameters.get>(); + cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); + model.parameters.set>(1.0); + model.parameters.set>(1.0); + + double t0 = 0.0, tmax = 0.5, dt = 0.5; + std::unique_ptr> integrator = std::make_unique>(); + auto sim = simulate_flows(t0, tmax, dt, model, std::move(integrator)); + EXPECT_EQ(sim[0].get_num_time_points(), 2); + EXPECT_EQ(sim[1].get_num_time_points(), 2); + // Flow vector size should be 6 for one age group + EXPECT_EQ(sim[1].get_last_value().size(), 6); +} From 8be7ebd5b86056484b56a2c26d50588160333c91 Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 09:50:42 +0100 Subject: [PATCH 06/16] [ci skip] example review suggestions --- cpp/examples/ode_seirv_example.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/cpp/examples/ode_seirv_example.cpp b/cpp/examples/ode_seirv_example.cpp index eec8b0c13c..60965f568f 100644 --- a/cpp/examples/ode_seirv_example.cpp +++ b/cpp/examples/ode_seirv_example.cpp @@ -28,6 +28,9 @@ * We assume that no one is initially infected or exposed at the beginning of the season. * Infections are seeded later via the OutsideFoI parameter. * The destinction into the 2 layers of susceptibility is done via the parameters. + * @tparam FP floating point type + * @param[in, out] model SEIRV model + * @param[in] total_pop total population size */ template void set_initial_population(mio::oseirv::Model& model, const FP total_pop) @@ -77,18 +80,18 @@ int main() mio::oseirv::Model model((int)num_groups); auto& parameters = model.parameters; - FP cont_freq = 10.0; - FP fact = FP(1) / FP(num_groups); + FP cont_freq = 10.0; + FP cont_freq_group = FP(1) / FP(num_groups); // Healthy contacts mio::ContactMatrixGroup& contacts_healthy = parameters.get>(); contacts_healthy[0] = mio::ContactMatrix( - Eigen::MatrixXd::Constant((Eigen::Index)num_groups, (Eigen::Index)num_groups, fact * cont_freq)); + Eigen::MatrixXd::Constant((Eigen::Index)num_groups, (Eigen::Index)num_groups, cont_freq_group * cont_freq)); // Sick contacts (here 20% reduced) mio::ContactMatrixGroup& contacts_sick = parameters.get>(); - contacts_sick[0] = mio::ContactMatrix( - Eigen::MatrixXd::Constant((Eigen::Index)num_groups, (Eigen::Index)num_groups, fact * cont_freq * 0.8)); + contacts_sick[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant( + (Eigen::Index)num_groups, (Eigen::Index)num_groups, cont_freq_group * cont_freq * 0.8)); // Parameters parameters.get>() = 1.2; From dbaba7345334ab7a4b7395a14f87b194593024bd Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 10:05:25 +0100 Subject: [PATCH 07/16] [ci skip] readme and parameters model.h --- cpp/models/ode_seirv/README.md | 6 +++ cpp/models/ode_seirv/model.h | 74 +++++++++++++++++----------------- 2 files changed, 42 insertions(+), 38 deletions(-) diff --git a/cpp/models/ode_seirv/README.md b/cpp/models/ode_seirv/README.md index ced80a5a1a..0d2e2fc65e 100644 --- a/cpp/models/ode_seirv/README.md +++ b/cpp/models/ode_seirv/README.md @@ -9,6 +9,12 @@ model structure from: > **BMC Infectious Diseases** 17, 258 (2017). > https://doi.org/10.1186/s12879-017-2344-6 +## Get started + +To get started, check out the +[official documentation](https://memilio.readthedocs.io/en/latest/cpp/models/oseirv.html) +or the [code example](../../examples/ode_seirv_example.cpp). + ## Reference Please cite the original study if you use this model: diff --git a/cpp/models/ode_seirv/model.h b/cpp/models/ode_seirv/model.h index d047e4e96f..3f372ce06b 100644 --- a/cpp/models/ode_seirv/model.h +++ b/cpp/models/ode_seirv/model.h @@ -76,50 +76,49 @@ class Model const auto& params = this->parameters; const size_t num_groups = (size_t)params.get_num_groups(); - // Contact matrices - const auto cm_h_expr = + // Get effective contact matrix effective_contacts = contacts_healthy + m*contacts_sick + Eigen::MatrixX contacts_healthy = params.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); - const auto cm_s_expr = + Eigen::MatrixX contacts_sick = params.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); - - // Get effective contact matrix B = H + m*S - Eigen::MatrixX H = cm_h_expr; - Eigen::MatrixX S = cm_s_expr; - const FP m = params.template get>(); - Eigen::MatrixX B = H + m * S; + const FP m = params.template get>(); + Eigen::MatrixX effective_contacts = contacts_healthy + m * contacts_sick; // Normalization ν(m): prepare susceptibility scaling and next-generation matrix (NG) // sigma_i = Phi * SigmaByAge[i] (age-specific susceptibility scaled by a global factor) - Eigen::VectorX sigma(num_groups); + Eigen::VectorX susceptibility(num_groups); for (size_t i = 0; i < num_groups; ++i) - sigma[(Eigen::Index)i] = + susceptibility[(Eigen::Index)i] = params.template get>()[AgeGroup(i)] * params.template get>(); // Sigma is a diagonal matrix applying susceptibility on the receiving (row) side of contacts. - Eigen::MatrixX Sigma = sigma.asDiagonal(); + Eigen::MatrixX Sigma = susceptibility.asDiagonal(); - // NG = Sigma * B is the (unnormalized) next-generation matrix. Its spectral radius will be used - // below to scale contacts B so that structure (mixing pattern) and transmissibility magnitude (BaselineTransmissibility) + // NG = Sigma * effective_contacts is the (unnormalized) next-generation matrix. Its spectral radius will be used + // below to scale contacts effective_contacts so that structure (mixing pattern) and transmissibility magnitude (BaselineTransmissibility) // are cleanly separated. - Eigen::MatrixX NG = Sigma * B; + Eigen::MatrixX NG = Sigma * effective_contacts; Eigen::ComplexEigenSolver> ces; ces.compute(NG); - FP nu = ces.eigenvalues().cwiseAbs().maxCoeff(); - if (nu > FP(0.)) { - B /= nu; + FP scale_transmissibility = ces.eigenvalues().cwiseAbs().maxCoeff(); + if (scale_transmissibility > FP(0.)) { + effective_contacts /= scale_transmissibility; } - const FP Re = + const FP transmissibility_baseline = params.template get>(); // baseline transmissibility scaling (R-like factor) - const FP gamma = params.template get>(); // progression/recovery rate (1 / mean duration) - const FP delta = params.template get>(); // amplitude of seasonal modulation - const FP tz = params.template get>(); // base phase shift for seasonality - const FP ts = params.template get>(); // additional phase shift (e.g. scenario-specific) - const FP lambda0 = params.template get>(); // constant external force of infection - const FP rho = params.template get< - ClusteringRho>(); // clustering exponent; rho<1 dampens, rho>1 amplifies prevalence effect - - const FP season = std::exp(delta * std::sin(FP(2) * std::numbers::pi_v * (t / FP(52.0) - tz + ts))); + const FP recovery_rate = params.template get>(); // progression/recovery rate (1 / mean duration) + const FP season_amplitude = params.template get>(); // amplitude of seasonal modulation + const FP season_shift_per_subtype = params.template get>(); // base phase shift for seasonality + const FP season_shift_per_season = + params.template get>(); // additional phase shift (e.g. season-specific) + const FP outside_foi = params.template get>(); // constant external force of infection + const FP clustering_exp = params.template get>(); // clustering exponent; clustering_exp<1 dampens, clustering_exp>1 amplifies prevalence effect + + const FP season = + std::exp(season_amplitude * std::sin(FP(2) * std::numbers::pi_v * + (t / FP(52.0) - season_shift_per_subtype + season_shift_per_season))); for (auto i : make_index_range((mio::AgeGroup)num_groups)) { // Flat indices @@ -130,7 +129,6 @@ class Model const size_t Ii = this->populations.get_flat_index({i, InfectionState::Infected}); const size_t IVi = this->populations.get_flat_index({i, InfectionState::InfectedVaccinated}); - // λ_i(t) = Re*γ*season * sum_j B_{ij} * ((I_j+I^V_j)/N_j)^ρ + λ0 FP sum = FP(0); for (auto j : make_index_range((mio::AgeGroup)num_groups)) { const size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible}); @@ -144,26 +142,26 @@ class Model const FP Nj = pop[Sj] + pop[SVj] + pop[Ej] + pop[EVj] + pop[Ij] + pop[IVj] + pop[Rj] + pop[RVj]; const FP inf_frac = (Nj > FP(0)) ? ((pop[Ij] + pop[IVj]) / Nj) : FP(0); - const FP term = (rho == FP(1)) ? inf_frac : std::pow(inf_frac, rho); - sum += B(i.get(), j.get()) * term; + const FP transmission_rate = (clustering_exp == FP(1)) ? inf_frac : std::pow(inf_frac, clustering_exp); + sum += effective_contacts(i.get(), j.get()) * transmission_rate; } - const FP lambda_i = Re * gamma * season * sum + lambda0; + const FP foi_i = transmissibility_baseline * recovery_rate * season * sum + lambda0; // Flows: S->E, S^V->E^V, E->I, E^V->I^V, I->R, I^V->R^V flows[Base::template get_flat_flow_index(i)] = - lambda_i * y[Si]; + foi_i * y[Si]; flows[Base::template get_flat_flow_index(i)] = lambda_i * y[SVi]; + InfectionState::ExposedVaccinated>(i)] = foi_i * y[SVi]; flows[Base::template get_flat_flow_index(i)] = - gamma * y[Ei]; + recovery_rate * y[Ei]; flows[Base::template get_flat_flow_index(i)] = gamma * y[EVi]; + InfectionState::InfectedVaccinated>(i)] = recovery_rate * y[EVi]; flows[Base::template get_flat_flow_index(i)] = - gamma * y[Ii]; + recovery_rate * y[Ii]; flows[Base::template get_flat_flow_index(i)] = gamma * y[IVi]; + InfectionState::RecoveredVaccinated>(i)] = recovery_rate * y[IVi]; } } }; From 9b868e459931b378092ebc284b6c59209d3569e1 Mon Sep 17 00:00:00 2001 From: Henrik Zunker <69154294+HenrZu@users.noreply.github.com> Date: Tue, 28 Oct 2025 10:08:00 +0100 Subject: [PATCH 08/16] Apply suggestions from code review Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com> --- cpp/models/ode_seirv/parameters.h | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/cpp/models/ode_seirv/parameters.h b/cpp/models/ode_seirv/parameters.h index 26aecc0178..ecbbca6eca 100644 --- a/cpp/models/ode_seirv/parameters.h +++ b/cpp/models/ode_seirv/parameters.h @@ -53,11 +53,11 @@ struct BaselineTransmissibility { /** * @brief Transition/recovery rate in 1/week. * - * Controls E→I and I→Rand is included in the force of infection to decouple from R_e. + * Controls E→I and I→R and is included in the force of infection to decouple the recovery rate from R_e. * Use the model time unit (weeks). Default: 1/2.0 */ template -struct Gamma { +struct RecoveryRate { using Type = UncertainValue; static Type get_default(AgeGroup) { @@ -91,10 +91,10 @@ struct SeasonalityAmplitude { /** * @brief Subtype-specific seasonal shift t_z. * - * Phase shift for the sinusoid in the force of incetion.; may be adjusted by a seasonal-shift correction if implemented. Default: 0.0. + * Phase shift for the sinusoid in the force of infection; may be adjusted by a seasonal-shift correction if implemented. Default: 0.0. */ template -struct ShiftTZ { +struct SeasonalityShiftPerSubtype { using Type = UncertainValue; static Type get_default(AgeGroup) { @@ -112,7 +112,7 @@ struct ShiftTZ { * Additional small seasonal timing correction per season. Default: 0.0. */ template -struct ShiftTS { +struct SeasonalityShiftPerSeason{ using Type = UncertainValue; static Type get_default(AgeGroup) { @@ -128,7 +128,7 @@ struct ShiftTS { * @brief Outside force of infection in 1/week. * * Additive external hazard that can seed infections when no infectives are present. - * Used for the beginning of the season (to import infections) or to model travel/imported cases. + * Used in the beginning of the season (to import infections) or to model travel/imported cases. * Default: 0.0. */ template @@ -150,7 +150,7 @@ struct OutsideFoI { * Exponent on the infectious fraction in force of infection. Must be > 0. Default: 1.0. */ template -struct ClusteringRho { +struct ClusteringExponent{ using Type = UncertainValue; static Type get_default(AgeGroup) { @@ -169,7 +169,7 @@ struct ClusteringRho { * Combines share of symptomatic cases and their higher infectiousness. Default: 1.0. */ template -struct SickMixingM { +struct SickMixing { using Type = UncertainValue; static Type get_default(AgeGroup) { @@ -225,7 +225,7 @@ struct ContactPatternsSick { * Default: 1.0 (fully susceptible). */ template -struct SigmaByAge { +struct SusceptibilityByAge { using Type = CustomIndexArray, AgeGroup>; static Type get_default(AgeGroup size) { @@ -275,10 +275,10 @@ struct VaccineEffectiveness { }; /** - * @brief phi: fraction of the population that remains susceptible at t0 (dimensionless, typically in [0,1]). + * @brief Fraction of the population that remains susceptible at t0 phi (dimensionless, typically in [0,1]). */ template -struct Phi { +struct SusceptibleFraction { using Type = UncertainValue; static Type get_default(AgeGroup) { From 394bf7540227b128ed5df27b9b6fd9ab7c7da86c Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 11:02:14 +0100 Subject: [PATCH 09/16] adjust params --- cpp/examples/ode_seirv_example.cpp | 22 ++++++++++----------- cpp/models/ode_seirv/model.h | 17 ++++++++-------- cpp/models/ode_seirv/parameters.h | 31 +++++++++++++++--------------- cpp/tests/test_odeseirv.cpp | 18 ++++++++--------- docs/source/cpp/models/oseirv.rst | 14 +++++++------- 5 files changed, 52 insertions(+), 50 deletions(-) diff --git a/cpp/examples/ode_seirv_example.cpp b/cpp/examples/ode_seirv_example.cpp index 60965f568f..166b90bc7d 100644 --- a/cpp/examples/ode_seirv_example.cpp +++ b/cpp/examples/ode_seirv_example.cpp @@ -46,8 +46,8 @@ void set_initial_population(mio::oseirv::Model& model, const FP total_pop) // Total population N_i as currently stored in group totals FP Ni = pop.get_group_total(mio::AgeGroup(i)); - FP sigma = params.template get>()[mio::AgeGroup(i)]; - FP phi = params.template get>(); + FP sigma = params.template get>()[mio::AgeGroup(i)]; + FP phi = params.template get>(); FP VC = params.template get>()[mio::AgeGroup(i)]; FP VE = params.template get>()[mio::AgeGroup(i)]; @@ -94,17 +94,17 @@ int main() (Eigen::Index)num_groups, (Eigen::Index)num_groups, cont_freq_group * cont_freq * 0.8)); // Parameters - parameters.get>() = 1.2; - parameters.get>() = 1.0 / 2.0; - parameters.get>() = 1.0; - parameters.get>() = 0.0; - parameters.get>() = 0.0; - parameters.get>() = 1e-6; - parameters.get>() = 0.9; - parameters.get>() = 2.0; + parameters.get>() = 1.2; + parameters.get>() = 1.0 / 2.0; + parameters.get>() = 1.0; + parameters.get>() = 0.0; + parameters.get>() = 0.0; + parameters.get>() = 1e-6; + parameters.get>() = 0.9; + parameters.get>() = 2.0; for (size_t i = 0; i < num_groups; ++i) { - parameters.get>()[mio::AgeGroup(i)] = 1.0; + parameters.get>()[mio::AgeGroup(i)] = 1.0; parameters.get>()[mio::AgeGroup(i)] = 0.3; parameters.get>()[mio::AgeGroup(i)] = 0.5; model.populations.set_difference_from_group_total( diff --git a/cpp/models/ode_seirv/model.h b/cpp/models/ode_seirv/model.h index 3f372ce06b..a9c2b6a730 100644 --- a/cpp/models/ode_seirv/model.h +++ b/cpp/models/ode_seirv/model.h @@ -81,15 +81,14 @@ class Model params.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); Eigen::MatrixX contacts_sick = params.template get>().get_cont_freq_mat().get_matrix_at(SimulationTime(t)); - const FP m = params.template get>(); + const FP m = params.template get>(); Eigen::MatrixX effective_contacts = contacts_healthy + m * contacts_sick; // Normalization ν(m): prepare susceptibility scaling and next-generation matrix (NG) - // sigma_i = Phi * SigmaByAge[i] (age-specific susceptibility scaled by a global factor) Eigen::VectorX susceptibility(num_groups); for (size_t i = 0; i < num_groups; ++i) - susceptibility[(Eigen::Index)i] = - params.template get>()[AgeGroup(i)] * params.template get>(); + susceptibility[(Eigen::Index)i] = params.template get>()[AgeGroup(i)] * + params.template get>(); // Sigma is a diagonal matrix applying susceptibility on the receiving (row) side of contacts. Eigen::MatrixX Sigma = susceptibility.asDiagonal(); @@ -107,13 +106,15 @@ class Model const FP transmissibility_baseline = params.template get>(); // baseline transmissibility scaling (R-like factor) - const FP recovery_rate = params.template get>(); // progression/recovery rate (1 / mean duration) + const FP recovery_rate = + params.template get>(); // progression/recovery rate (1 / mean duration) const FP season_amplitude = params.template get>(); // amplitude of seasonal modulation - const FP season_shift_per_subtype = params.template get>(); // base phase shift for seasonality + const FP season_shift_per_subtype = + params.template get>(); // base phase shift for seasonality const FP season_shift_per_season = - params.template get>(); // additional phase shift (e.g. season-specific) + params.template get>(); // additional phase shift (e.g. season-specific) const FP outside_foi = params.template get>(); // constant external force of infection - const FP clustering_exp = params.template get>(); // clustering exponent; clustering_exp<1 dampens, clustering_exp>1 amplifies prevalence effect const FP season = diff --git a/cpp/models/ode_seirv/parameters.h b/cpp/models/ode_seirv/parameters.h index ecbbca6eca..81e34ae876 100644 --- a/cpp/models/ode_seirv/parameters.h +++ b/cpp/models/ode_seirv/parameters.h @@ -65,7 +65,7 @@ struct RecoveryRate { } static std::string name() { - return "Gamma"; + return "RecoveryRate"; } }; @@ -102,7 +102,7 @@ struct SeasonalityShiftPerSubtype { } static std::string name() { - return "ShiftTZ"; + return "SeasonalityShiftPerSubtype"; } }; @@ -112,7 +112,7 @@ struct SeasonalityShiftPerSubtype { * Additional small seasonal timing correction per season. Default: 0.0. */ template -struct SeasonalityShiftPerSeason{ +struct SeasonalityShiftPerSeason { using Type = UncertainValue; static Type get_default(AgeGroup) { @@ -120,7 +120,7 @@ struct SeasonalityShiftPerSeason{ } static std::string name() { - return "ShiftTS"; + return "SeasonalityShiftPerSeason"; } }; @@ -150,7 +150,7 @@ struct OutsideFoI { * Exponent on the infectious fraction in force of infection. Must be > 0. Default: 1.0. */ template -struct ClusteringExponent{ +struct ClusteringExponent { using Type = UncertainValue; static Type get_default(AgeGroup) { @@ -158,7 +158,7 @@ struct ClusteringExponent{ } static std::string name() { - return "ClusteringRho"; + return "ClusteringExponent"; } }; @@ -177,7 +177,7 @@ struct SickMixing { } static std::string name() { - return "SickMixingM"; + return "SickMixing"; } }; @@ -233,7 +233,7 @@ struct SusceptibilityByAge { } static std::string name() { - return "SigmaByAge"; + return "CustomIndexArray"; } }; @@ -286,15 +286,16 @@ struct SusceptibleFraction { } static std::string name() { - return "Phi"; + return "SusceptibleFraction"; } }; template using ParametersBase = - ParameterSet, Gamma, SeasonalityAmplitude, ShiftTZ, ShiftTS, - OutsideFoI, ClusteringRho, SickMixingM, ContactPatternsHealthy, - ContactPatternsSick, SigmaByAge, VaccineCoverage, VaccineEffectiveness, Phi>; + ParameterSet, RecoveryRate, SeasonalityAmplitude, + SeasonalityShiftPerSubtype, SeasonalityShiftPerSeason, OutsideFoI, ClusteringExponent, + SickMixing, ContactPatternsHealthy, ContactPatternsSick, CustomIndexArray, + VaccineCoverage, VaccineEffectiveness, SusceptibleFraction>; /** * @brief Parameter set for the age-resolved SEIRV model (S,E,I,R plus vaccinated states) as per the appendix. @@ -332,9 +333,9 @@ class Parameters : public ParametersBase bool apply_constraints() { bool corrected = false; - if (this->template get>() <= 0.0) { - this->template get>() = 1.0; - corrected = true; + if (this->template get>() <= 0.0) { + this->template get>() = 1.0; + corrected = true; } if (this->template get>() < 0.0) { this->template get>() = 0.0; diff --git a/cpp/tests/test_odeseirv.cpp b/cpp/tests/test_odeseirv.cpp index 38985ec5e4..a91b07a3e2 100644 --- a/cpp/tests/test_odeseirv.cpp +++ b/cpp/tests/test_odeseirv.cpp @@ -75,7 +75,7 @@ TEST(TestOdeSeirv, populationConservation) cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); model.parameters.set>(1.0); - model.parameters.set>(1.0); + model.parameters.set>(1.0); double t0 = 0.0, tmax = 5.0, dt = 0.5; auto sim = simulate(t0, tmax, dt, model); @@ -89,11 +89,11 @@ TEST(TestOdeSeirv, applyConstraints) // First: defaults should need no correction EXPECT_FALSE(model.parameters.apply_constraints()); - model.parameters.set>(0.0); // invalid, must become >0 (1.0) + model.parameters.set>(0.0); // invalid, must become >0 (1.0) model.parameters.set>(-2.); // invalid -> 0 model.parameters.set>(-0.5); // invalid -> 0 EXPECT_TRUE(model.parameters.apply_constraints()); - EXPECT_GT(model.parameters.get>(), 0.0); + EXPECT_GT(model.parameters.get>(), 0.0); EXPECT_EQ(model.parameters.get>(), 0.0); EXPECT_EQ(model.parameters.get>(), 0.0); } @@ -117,8 +117,8 @@ TEST(TestOdeSeirv, flowsSingleAgeGroup) mio::ContactMatrixGroup& cm_s = model.parameters.get>(); cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); model.parameters.set>(1.0); - model.parameters.set>(1.0); - model.parameters.set>(1.0); + model.parameters.set>(1.0); + model.parameters.set>(1.0); model.parameters.set>(0.0); model.parameters.set>(0.0); @@ -171,8 +171,8 @@ TEST(TestOdeSeirv, flowsTwoAgeGroupsIdentityContacts) mio::ContactMatrixGroup& cm_s = model.parameters.get>(); cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Zero(2, 2)); model.parameters.set>(1.0); - model.parameters.set>(1.0); - model.parameters.set>(1.0); + model.parameters.set>(1.0); + model.parameters.set>(1.0); model.parameters.set>(0.0); model.parameters.set>(0.0); @@ -250,7 +250,7 @@ TEST(TestOdeSeirv, simulationEuler) mio::ContactMatrixGroup& cm_s = model.parameters.get>(); cm_s[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 0.)); model.parameters.set>(1.0); - model.parameters.set>(1.0); + model.parameters.set>(1.0); double t0 = 0.0, tmax = 1.0, dt = 0.5; std::unique_ptr> integrator = std::make_unique>(); @@ -272,7 +272,7 @@ TEST(TestOdeSeirv, flowSimulation) mio::ContactMatrixGroup& cm_h = model.parameters.get>(); cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); model.parameters.set>(1.0); - model.parameters.set>(1.0); + model.parameters.set>(1.0); double t0 = 0.0, tmax = 0.5, dt = 0.5; std::unique_ptr> integrator = std::make_unique>(); diff --git a/docs/source/cpp/models/oseirv.rst b/docs/source/cpp/models/oseirv.rst index 539f05d676..cd0591a013 100644 --- a/docs/source/cpp/models/oseirv.rst +++ b/docs/source/cpp/models/oseirv.rst @@ -73,25 +73,25 @@ The model uses the following parameters (time unit: week): - ``BaselineTransmissibility`` - Baseline transmissibility (dimensionless); scales the normalized force of infection. * - :math:`\gamma` - - ``Gamma`` + - ``RecoveryRate`` - Transition / recovery rate (1/week) for E -> I and I -> R. * - :math:`\delta` - ``SeasonalityAmplitude`` - Amplitude of the seasonal modulation :math:`\exp(\delta\,\sin(2\pi(t/52 - t_z + t_s)))`. * - :math:`t_z` - - ``ShiftTZ`` + - ``SeasonalityShiftPerSubtype`` - Coarse (subtype-specific) seasonal phase shift. * - :math:`t_s` - - ``ShiftTS`` + - ``SeasonalityShiftPerSeason`` - Fine seasonal phase adjustment per season. * - :math:`\lambda_0` - ``OutsideFoI`` - External (additive) force of infection, can seed infections. * - :math:`\rho` - - ``ClusteringRho`` + - ``ClusteringExponent`` - Clustering exponent on the infectious fraction. * - :math:`m` - - ``SickMixingM`` + - ``SickMixing`` - Mixing weight for symptomatic (“sick”) contacts in the blended contact matrix. * - :math:`C^{H}` - ``ContactPatternsHealthy`` @@ -100,7 +100,7 @@ The model uses the following parameters (time unit: week): - ``ContactPatternsSick`` - Age-structured contact matrix (symptomatic), combined using :math:`m`. * - :math:`\sigma_i` - - ``SigmaByAge`` + - ``CustomIndexArray`` - Age-specific baseline susceptibility (pre-existing immunity modifier). * - :math:`VC_i` - ``VaccineCoverage`` @@ -109,7 +109,7 @@ The model uses the following parameters (time unit: week): - ``VaccineEffectiveness`` - Vaccine effectiveness (reducing effective susceptibility). * - :math:`\phi_0` - - ``Phi`` + - ``SusceptibleFraction`` - Fraction of the total population forming the effectively susceptible pool at :math:`t_0`. Note: ``VaccineCoverage`` and ``VaccineEffectiveness`` are only used for initialization. Transitions presently From e497f1e04d958da2296a5c9745256996690db8de Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 11:13:33 +0100 Subject: [PATCH 10/16] [ci skip] update tests --- cpp/tests/test_odeseirv.cpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/cpp/tests/test_odeseirv.cpp b/cpp/tests/test_odeseirv.cpp index a91b07a3e2..5d0d997ea8 100644 --- a/cpp/tests/test_odeseirv.cpp +++ b/cpp/tests/test_odeseirv.cpp @@ -93,7 +93,7 @@ TEST(TestOdeSeirv, applyConstraints) model.parameters.set>(-2.); // invalid -> 0 model.parameters.set>(-0.5); // invalid -> 0 EXPECT_TRUE(model.parameters.apply_constraints()); - EXPECT_GT(model.parameters.get>(), 0.0); + EXPECT_EQ(model.parameters.get>(), 0.0); EXPECT_EQ(model.parameters.get>(), 0.0); EXPECT_EQ(model.parameters.get>(), 0.0); } @@ -111,7 +111,7 @@ TEST(TestOdeSeirv, flowsSingleAgeGroup) model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = I; model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::InfectedVaccinated}] = IV; - // Parameters: identity contacts, gamma=1, Re=1, others zero ⇒ λ = (I+IV)/N + // Parameters: identity contacts, RecoveryRate=1, ClusteringExponent=1, BaselineTransmissibility=1, others zero ⇒ λ = (I+IV)/N mio::ContactMatrixGroup& cm_h = model.parameters.get>(); cm_h[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1)); mio::ContactMatrixGroup& cm_s = model.parameters.get>(); @@ -132,10 +132,10 @@ TEST(TestOdeSeirv, flowsSingleAgeGroup) const double lambda = (I + IV) / N; // expected force of infection const double f_SE = S * lambda; const double f_SV_EV = SV * lambda; - const double f_E_I = E; // gamma * E - const double f_EV_IV = EV; // gamma * EV - const double f_I_R = I; // gamma * I - const double f_IV_RV = IV; // gamma * IV + const double f_E_I = RecoveryRate * E; + const double f_EV_IV = RecoveryRate * EV; + const double f_I_R = RecoveryRate * I; + const double f_IV_RV = RecoveryRate * IV; auto idx_SE = model.template get_flat_flow_index(mio::AgeGroup(0)); @@ -144,7 +144,7 @@ TEST(TestOdeSeirv, flowsSingleAgeGroup) auto idx_EI = model.template get_flat_flow_index( mio::AgeGroup(0)); - auto idx_EVI = + auto idx_EIV = model.template get_flat_flow_index(mio::AgeGroup(0)); auto idx_IR = model.template get_flat_flow_index model(2); mio::ContactMatrixGroup& cm_h = model.parameters.get>(); - // Each group only contacts itself so that λ_i = I_i / N_i + // Let each group have only contacts with itself and set other parameters, so that λ_i = I_i / N_i Eigen::MatrixXd Id2 = Eigen::MatrixXd::Identity(2, 2); cm_h[0] = mio::ContactMatrix(Id2); mio::ContactMatrixGroup& cm_s = model.parameters.get>(); @@ -176,6 +176,7 @@ TEST(TestOdeSeirv, flowsTwoAgeGroupsIdentityContacts) model.parameters.set>(0.0); model.parameters.set>(0.0); + // Only population in the non-vaccinated susceptible and infected compartments // Group 0 double S0 = 100, I0 = 20; // others zero model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Susceptible}] = S0; @@ -262,7 +263,7 @@ TEST(TestOdeSeirv, simulationEuler) } } -TEST(TestOdeSeirv, flowSimulation) +TEST(TestOdeSeirv, flowSimulationEuler) { mio::oseirv::Model model(1); model.populations[{mio::AgeGroup(0), mio::oseirv::InfectionState::Infected}] = 5; From ac047e93f8fcb140d571d22c061d045cb9a8b972 Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 11:18:22 +0100 Subject: [PATCH 11/16] fix typos --- cpp/models/ode_seirv/model.h | 4 ++-- cpp/models/ode_seirv/parameters.h | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/cpp/models/ode_seirv/model.h b/cpp/models/ode_seirv/model.h index a9c2b6a730..54a0addc7c 100644 --- a/cpp/models/ode_seirv/model.h +++ b/cpp/models/ode_seirv/model.h @@ -87,7 +87,7 @@ class Model // Normalization ν(m): prepare susceptibility scaling and next-generation matrix (NG) Eigen::VectorX susceptibility(num_groups); for (size_t i = 0; i < num_groups; ++i) - susceptibility[(Eigen::Index)i] = params.template get>()[AgeGroup(i)] * + susceptibility[(Eigen::Index)i] = params.template get>()[AgeGroup(i)] * params.template get>(); // Sigma is a diagonal matrix applying susceptibility on the receiving (row) side of contacts. Eigen::MatrixX Sigma = susceptibility.asDiagonal(); @@ -146,7 +146,7 @@ class Model const FP transmission_rate = (clustering_exp == FP(1)) ? inf_frac : std::pow(inf_frac, clustering_exp); sum += effective_contacts(i.get(), j.get()) * transmission_rate; } - const FP foi_i = transmissibility_baseline * recovery_rate * season * sum + lambda0; + const FP foi_i = transmissibility_baseline * recovery_rate * season * sum + outside_foi; // Flows: S->E, S^V->E^V, E->I, E^V->I^V, I->R, I^V->R^V flows[Base::template get_flat_flow_index(i)] = diff --git a/cpp/models/ode_seirv/parameters.h b/cpp/models/ode_seirv/parameters.h index 81e34ae876..e8149230fc 100644 --- a/cpp/models/ode_seirv/parameters.h +++ b/cpp/models/ode_seirv/parameters.h @@ -233,7 +233,7 @@ struct SusceptibilityByAge { } static std::string name() { - return "CustomIndexArray"; + return "SusceptibilityByAge"; } }; @@ -294,7 +294,7 @@ template using ParametersBase = ParameterSet, RecoveryRate, SeasonalityAmplitude, SeasonalityShiftPerSubtype, SeasonalityShiftPerSeason, OutsideFoI, ClusteringExponent, - SickMixing, ContactPatternsHealthy, ContactPatternsSick, CustomIndexArray, + SickMixing, ContactPatternsHealthy, ContactPatternsSick, SusceptibilityByAge, VaccineCoverage, VaccineEffectiveness, SusceptibleFraction>; /** From 72a6832aa309fa41e77b4ec913e0b57183bb5850 Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 13:11:55 +0100 Subject: [PATCH 12/16] fix notation in example --- cpp/examples/ode_seirv_example.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/cpp/examples/ode_seirv_example.cpp b/cpp/examples/ode_seirv_example.cpp index 166b90bc7d..d9f1adbfca 100644 --- a/cpp/examples/ode_seirv_example.cpp +++ b/cpp/examples/ode_seirv_example.cpp @@ -46,23 +46,23 @@ void set_initial_population(mio::oseirv::Model& model, const FP total_pop) // Total population N_i as currently stored in group totals FP Ni = pop.get_group_total(mio::AgeGroup(i)); - FP sigma = params.template get>()[mio::AgeGroup(i)]; - FP phi = params.template get>(); - FP VC = params.template get>()[mio::AgeGroup(i)]; - FP VE = params.template get>()[mio::AgeGroup(i)]; + FP s_age = params.template get>()[mio::AgeGroup(i)]; + FP s_frac = params.template get>(); + FP vc = params.template get>()[mio::AgeGroup(i)]; + FP ve = params.template get>()[mio::AgeGroup(i)]; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Exposed}] = 0; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::ExposedVaccinated}] = 0; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Infected}] = 0; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::InfectedVaccinated}] = 0; - pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}] = phi * sigma * (FP(1) - VC) * Ni; + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}] = s_frac * s_age * (FP(1) - vc) * Ni; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::SusceptibleVaccinated}] = - phi * sigma * (FP(1) - VE) * VC * Ni; + s_frac * s_age * (FP(1) - ve) * vc * Ni; - pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Recovered}] = (FP(1) - phi * sigma) * (FP(1) - VC) * Ni; + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Recovered}] = (FP(1) - s_frac * s_age) * (FP(1) - vc) * Ni; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::RecoveredVaccinated}] = - (FP(1) - phi * sigma * (FP(1) - VE)) * VC * Ni; + (FP(1) - s_frac * s_age * (FP(1) - ve)) * vc * Ni; } } @@ -104,7 +104,7 @@ int main() parameters.get>() = 2.0; for (size_t i = 0; i < num_groups; ++i) { - parameters.get>()[mio::AgeGroup(i)] = 1.0; + parameters.get>()[mio::AgeGroup(i)] = 1.0; parameters.get>()[mio::AgeGroup(i)] = 0.3; parameters.get>()[mio::AgeGroup(i)] = 0.5; model.populations.set_difference_from_group_total( From c12cd0a5766e60f91f4ede8ce4a038e9b2704eae Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 13:12:32 +0100 Subject: [PATCH 13/16] [ci skip] format --- cpp/examples/ode_seirv_example.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/examples/ode_seirv_example.cpp b/cpp/examples/ode_seirv_example.cpp index d9f1adbfca..46b5928369 100644 --- a/cpp/examples/ode_seirv_example.cpp +++ b/cpp/examples/ode_seirv_example.cpp @@ -46,10 +46,10 @@ void set_initial_population(mio::oseirv::Model& model, const FP total_pop) // Total population N_i as currently stored in group totals FP Ni = pop.get_group_total(mio::AgeGroup(i)); - FP s_age = params.template get>()[mio::AgeGroup(i)]; - FP s_frac = params.template get>(); - FP vc = params.template get>()[mio::AgeGroup(i)]; - FP ve = params.template get>()[mio::AgeGroup(i)]; + FP s_age = params.template get>()[mio::AgeGroup(i)]; + FP s_frac = params.template get>(); + FP vc = params.template get>()[mio::AgeGroup(i)]; + FP ve = params.template get>()[mio::AgeGroup(i)]; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Exposed}] = 0; pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::ExposedVaccinated}] = 0; @@ -104,7 +104,7 @@ int main() parameters.get>() = 2.0; for (size_t i = 0; i < num_groups; ++i) { - parameters.get>()[mio::AgeGroup(i)] = 1.0; + parameters.get>()[mio::AgeGroup(i)] = 1.0; parameters.get>()[mio::AgeGroup(i)] = 0.3; parameters.get>()[mio::AgeGroup(i)] = 0.5; model.populations.set_difference_from_group_total( From 82505c825798d2035f942f126c0892b749313c98 Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 14:03:30 +0100 Subject: [PATCH 14/16] . --- cpp/tests/test_odeseirv.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/cpp/tests/test_odeseirv.cpp b/cpp/tests/test_odeseirv.cpp index 5d0d997ea8..3901032a95 100644 --- a/cpp/tests/test_odeseirv.cpp +++ b/cpp/tests/test_odeseirv.cpp @@ -128,14 +128,15 @@ TEST(TestOdeSeirv, flowsSingleAgeGroup) flows.setZero(); model.get_flows(pop, y0, 0.0, flows); + const auto recovery_rate = model.parameters.get>(); const double N = S + SV + E + EV + I + IV; // (R,RV = 0) const double lambda = (I + IV) / N; // expected force of infection const double f_SE = S * lambda; const double f_SV_EV = SV * lambda; - const double f_E_I = RecoveryRate * E; - const double f_EV_IV = RecoveryRate * EV; - const double f_I_R = RecoveryRate * I; - const double f_IV_RV = RecoveryRate * IV; + const double f_E_I = recovery_rate * E; + const double f_EV_IV = recovery_rate * EV; + const double f_I_R = recovery_rate * I; + const double f_IV_RV = recovery_rate * IV; auto idx_SE = model.template get_flat_flow_index(mio::AgeGroup(0)); From f979fecb0662f5aba022b6a7efb079dddbb68613 Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 16:38:48 +0100 Subject: [PATCH 15/16] . --- cpp/tests/test_odeseirv.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/tests/test_odeseirv.cpp b/cpp/tests/test_odeseirv.cpp index 3901032a95..349ab0ff0e 100644 --- a/cpp/tests/test_odeseirv.cpp +++ b/cpp/tests/test_odeseirv.cpp @@ -89,7 +89,7 @@ TEST(TestOdeSeirv, applyConstraints) // First: defaults should need no correction EXPECT_FALSE(model.parameters.apply_constraints()); - model.parameters.set>(0.0); // invalid, must become >0 (1.0) + model.parameters.set>(-0.5); // invalid, must become >0 (1.0) model.parameters.set>(-2.); // invalid -> 0 model.parameters.set>(-0.5); // invalid -> 0 EXPECT_TRUE(model.parameters.apply_constraints()); From 9fd5e08d015b1e0e5e15d2035283038dc5e513b1 Mon Sep 17 00:00:00 2001 From: Henrik Zunker Date: Tue, 28 Oct 2025 16:47:36 +0100 Subject: [PATCH 16/16] . --- cpp/tests/test_odeseirv.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/tests/test_odeseirv.cpp b/cpp/tests/test_odeseirv.cpp index 349ab0ff0e..5919da987b 100644 --- a/cpp/tests/test_odeseirv.cpp +++ b/cpp/tests/test_odeseirv.cpp @@ -93,7 +93,7 @@ TEST(TestOdeSeirv, applyConstraints) model.parameters.set>(-2.); // invalid -> 0 model.parameters.set>(-0.5); // invalid -> 0 EXPECT_TRUE(model.parameters.apply_constraints()); - EXPECT_EQ(model.parameters.get>(), 0.0); + EXPECT_EQ(model.parameters.get>(), 1.0); EXPECT_EQ(model.parameters.get>(), 0.0); EXPECT_EQ(model.parameters.get>(), 0.0); }