diff --git a/cpp/benchmarks/CMakeLists.txt b/cpp/benchmarks/CMakeLists.txt index 731fe8d566..cb2d6a8bbf 100755 --- a/cpp/benchmarks/CMakeLists.txt +++ b/cpp/benchmarks/CMakeLists.txt @@ -34,5 +34,8 @@ target_link_libraries(simulation_benchmark PRIVATE memilio ode_secir benchmark:: add_executable(graph_simulation_benchmark graph_simulation.cpp) target_link_libraries(graph_simulation_benchmark PRIVATE memilio ode_secirvvs benchmark::benchmark) +add_executable(secirvvs_advance_benchmark secirvvs_advance.cpp) +target_link_libraries(secirvvs_advance_benchmark PRIVATE memilio ode_secirvvs benchmark::benchmark) + add_executable(abm_benchmark abm.cpp) target_link_libraries(abm_benchmark PRIVATE abm benchmark::benchmark) diff --git a/cpp/benchmarks/flow_simulation_ode_secirvvs.h b/cpp/benchmarks/flow_simulation_ode_secirvvs.h index c1587283a0..be893af609 100644 --- a/cpp/benchmarks/flow_simulation_ode_secirvvs.h +++ b/cpp/benchmarks/flow_simulation_ode_secirvvs.h @@ -42,8 +42,8 @@ class FlowlessModel : public CompartmentalModel, - osecirvvs::Parameters>; + mio::Populations, + osecirvvs::Parameters>; public: FlowlessModel(const Populations& pop, const ParameterSet& params) @@ -539,7 +539,6 @@ class Simulation : public Base */ Eigen::Ref> advance(ScalarType tmax) { - auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis(); auto& dyn_npis = this->get_model().parameters.template get>(); auto& contact_patterns = this->get_model().parameters.template get>(); @@ -549,45 +548,36 @@ class Simulation : public Base this->get_model().parameters.template get>(); ScalarType delay_npi_implementation; - auto t = Base::get_result().get_last_time(); - const auto dt = dyn_npis.get_interval().get(); + auto t = Base::get_result().get_last_time(); while (t < tmax) { - auto dt_eff = std::min({dt, tmax - t, m_t_last_npi_check + dt - t}); - if (dt_eff >= 1.0) { - dt_eff = 1.0; + if (t > 0) { + delay_npi_implementation = ScalarType(dyn_npis.get_implementation_delay()); + } + else { // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay. + delay_npi_implementation = 0; } - if (t == 0) { //this->apply_vaccination(t); // done in init now? this->apply_variant(t, base_infectiousness); } - Base::advance(t + dt_eff); - if (t + 0.5 + dt_eff - std::floor(t + 0.5) >= 1) { - this->apply_vaccination(t + 0.5 + dt_eff); - this->apply_variant(t, base_infectiousness); - } - - if (t > 0) { - delay_npi_implementation = 7; - } - else { - delay_npi_implementation = 0; - } - t = t + dt_eff; if (dyn_npis.get_thresholds().size() > 0) { - if (floating_point_greater_equal(t, m_t_last_npi_check + dt)) { - if (t < t_end_dyn_npis) { - auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * - dyn_npis.get_base_value(); - auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); - if (exceeded_threshold != dyn_npis.get_thresholds().end() && - (exceeded_threshold->first > m_dynamic_npi.first || - t > ScalarType(m_dynamic_npi.second))) { //old npi was weaker or is expired - + ScalarType direc_begin = ScalarType(dyn_npis.get_directive_begin()); + ScalarType direc_end = ScalarType(dyn_npis.get_directive_end()); + if (floating_point_greater_equal(t, direc_begin, 1e-10) && t < direc_end) { + auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * + dyn_npis.get_base_value(); + auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); + if (exceeded_threshold != dyn_npis.get_thresholds().end() && + (exceeded_threshold->first > m_dynamic_npi.first || + t > ScalarType(m_dynamic_npi.second))) { // old npi was weaker or is expired + + if (t + delay_npi_implementation < direc_end) { auto t_start = SimulationTime(t + delay_npi_implementation); - auto t_end = t_start + SimulationTime(dyn_npis.get_duration()); + // set the end to the minimum of start+delay and the end of the directive + auto t_end = SimulationTime( + min(direc_end, ScalarType(t_start + dyn_npis.get_duration()))); this->get_model().parameters.get_start_commuter_detection() = t_start.get(); this->get_model().parameters.get_end_commuter_detection() = t_end.get(); m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end); @@ -597,12 +587,16 @@ class Simulation : public Base }); } } - m_t_last_npi_check = t; } } - else { - m_t_last_npi_check = t; + + auto dt_eff = min(1.0, tmax - t); + Base::advance(t + dt_eff); + if (t + 0.5 + dt_eff - std::floor(t + 0.5) >= 1) { + this->apply_vaccination(t + 0.5 + dt_eff); + this->apply_variant(t, base_infectiousness); } + t = t + dt_eff; } this->get_model().parameters.template get>() = @@ -719,6 +713,15 @@ void setup_model(Model& model) model.parameters.template get>()[AgeGroup(0)] = 0.9; model.parameters.template get>() = 0.2; + + auto& npis = model.parameters.template get>(); + auto npi_groups = Eigen::VectorXd::Ones(contact_matrix[0].get_num_groups()); + npis.set_threshold(0.01 * 100'000, {DampingSampling(0.5, DampingLevel(0), DampingType(0), + SimulationTime(0), {0}, npi_groups)}); + npis.set_base_value(100'000); + npis.set_implementation_delay(SimulationTime(7.0)); + npis.set_duration(SimulationTime(14.0)); + // The function apply_constraints() ensures that all parameters are within their defined bounds. // Note that negative values are set to zero instead of stopping the simulation. model.apply_constraints(); diff --git a/cpp/benchmarks/secirvvs_advance.cpp b/cpp/benchmarks/secirvvs_advance.cpp new file mode 100644 index 0000000000..e4a432371b --- /dev/null +++ b/cpp/benchmarks/secirvvs_advance.cpp @@ -0,0 +1,93 @@ +/* +* Copyright (C) 2020-2026 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. +*/ + +/** + * Benchmark comparing the overhead of the model-specific osecirvvs::Simulation::advance() + * (which uses substeps and calls apply_vaccination / apply_variant / dynamic NPI checks) + * versus the generic mio::Simulation::advance() (a single integrator call for the entire range). + */ + +#include "ode_secirvvs/model.h" +#include "memilio/compartments/simulation.h" +#include "memilio/utils/logging.h" +#include "benchmark/benchmark.h" + +using FP = double; +using Model = mio::osecirvvs::Model; + +static Model make_model(bool with_npis = false) +{ + constexpr int tmax = 10; + Model model(1); + model.populations[{mio::AgeGroup(0), mio::osecirvvs::InfectionState::InfectedSymptomsNaive}] = 100.0; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecirvvs::InfectionState::SusceptibleNaive}, + 10000.0); + model.parameters.get>().resize(mio::SimulationDay(size_t(tmax + 1))); + model.parameters.get>().resize(mio::SimulationDay(size_t(tmax + 1))); + if (with_npis) { + auto& npis = model.parameters.get>(); + npis.set_threshold(0.01 * 100'000, {mio::DampingSampling{1.0, + mio::DampingLevel(0), + mio::DampingType(0), + mio::SimulationTime(0), + {0}, + Eigen::VectorXd::Ones(1)}}); + npis.set_duration(mio::SimulationTime(14.0)); + npis.set_base_value(100'000); + } + return model; +} + +// Generic advance: single integrator call for the full [t0, tmax] range. +static void BM_generic(benchmark::State& state) +{ + mio::set_log_level(mio::LogLevel::off); + auto model = make_model(); + for (auto _ : state) { + mio::simulate(0., 10., 0.1, model); + } +} + +// Model-specific advance without dynamic NPIs: 1-day loop with apply_vaccination + apply_variant, +// dynamic NPI threshold check is skipped (thresholds empty). +static void BM_secirvvs_no_npis(benchmark::State& state) +{ + mio::set_log_level(mio::LogLevel::off); + auto model = make_model(/*with_npis=*/false); + for (auto _ : state) { + mio::osecirvvs::simulate(0., 10., 0.1, model); + } +} + +// Model-specific advance with dynamic NPIs: same as above plus get_infections_relative + +// threshold comparison on every day step. +static void BM_secirvvs_with_npis(benchmark::State& state) +{ + mio::set_log_level(mio::LogLevel::off); + auto model = make_model(/*with_npis=*/true); + for (auto _ : state) { + mio::osecirvvs::simulate(0., 10., 0.1, model); + } +} + +BENCHMARK(BM_generic)->Name("SECIRVVS generic advance"); +BENCHMARK(BM_secirvvs_no_npis)->Name("SECIRVVS model-specific advance (no dynamic NPIs)"); +BENCHMARK(BM_secirvvs_with_npis)->Name("SECIRVVS model-specific advance (with dynamic NPIs)"); +BENCHMARK_MAIN(); diff --git a/cpp/examples/ode_secirvvs.cpp b/cpp/examples/ode_secirvvs.cpp index cfdd39b19d..0dc318cbb3 100644 --- a/cpp/examples/ode_secirvvs.cpp +++ b/cpp/examples/ode_secirvvs.cpp @@ -79,7 +79,6 @@ int main() .get>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] = num_vaccinations; } - model.parameters.get>() = 7; auto& contacts = model.parameters.get>(); auto& contact_matrix = contacts.get_cont_freq_mat(); diff --git a/cpp/memilio/epidemiology/dynamic_npis.h b/cpp/memilio/epidemiology/dynamic_npis.h index 02467fa3c7..73f2a1e968 100644 --- a/cpp/memilio/epidemiology/dynamic_npis.h +++ b/cpp/memilio/epidemiology/dynamic_npis.h @@ -23,6 +23,8 @@ #include "memilio/epidemiology/damping_sampling.h" #include "memilio/utils/stl_util.h" +#include + namespace mio { @@ -54,7 +56,7 @@ class DynamicNPIs */ auto get_max_exceeded_threshold(FP value) { - //thresholds are sorted by value descending, so upper_bound returns the first threshold that is smaller using binary search + // thresholds are sorted by value descending, so upper_bound returns the first threshold that is smaller using binary search auto iter_max_exceeded_threshold = std::upper_bound(m_thresholds.begin(), m_thresholds.end(), value, [](auto& val, auto& t2) { return val > t2.first; @@ -97,22 +99,24 @@ class DynamicNPIs } /** - * Get/Set the interval at which the NPIs are checked. + * Get/Set the implementation delay at which the NPIs are implemented after threshold exceedance. + * This parameter imitates delayed reaction times when automatic implementations should be realized. * @{ */ /** - * @return the interval at which the NPIs are checked. + * @return the implementation delay after which the NPIs are implemented upon threshold exceedance. */ - SimulationTime get_interval() const + SimulationTime get_implementation_delay() const { - return m_interval; + return m_delay; } /** - * @param interval The interval at which the NPIs are checked. + * @param delay The implementation delay after which the NPIs are implemented upon threshold exceedance. */ - void set_interval(SimulationTime interval) + void set_implementation_delay(SimulationTime delay) { - m_interval = interval; + assert(delay >= SimulationTime(0.0) && "Implementation delay must be non-negative."); + m_delay = delay; } /**@}*/ @@ -130,7 +134,7 @@ class DynamicNPIs return m_base; } /** - * @return The base value of the thresholds. + * @param v Sets the base value of the thresholds. */ void set_base_value(FP v) { @@ -138,6 +142,50 @@ class DynamicNPIs } /**@}*/ + /** + * Get/Set the first day of the simulation for which a DynamicNPI *can* be activated. + * This parameter imitates the beginning date of a legal directive. + * @{ + */ + /** + * @return the first day of a legal directive prescribing the DynamicNPI. + */ + SimulationTime get_directive_begin() const + { + return m_directive_begin; + } + /** + * @param begin The first day of a legal directive prescribing the DynamicNPI. + */ + void set_directive_begin(SimulationTime begin) + { + assert(begin < m_directive_end && "Directive begin must be before directive end."); + m_directive_begin = begin; + } + /**@}*/ + + /** + * Get/Set the last day of the simulation for which a DynamicNPI *can* be active. + * This parameter imitates the last date of a legal directive and ends all active DynamicNPIs. + * @{ + */ + /** + * @return the last day of a legal directive prescribing the DynamicNPI. + */ + SimulationTime get_directive_end() const + { + return m_directive_end; + } + /** + * @param end The last day of a legal directive prescribing the DynamicNPI. + */ + void set_directive_end(SimulationTime end) + { + assert(m_directive_begin < end && "Directive end must be strictly after directive begin."); + m_directive_end = end; + } + /**@}*/ + /** * draw a random sample from the damping distributions */ @@ -160,8 +208,10 @@ class DynamicNPIs auto obj = io.create_object("DynamicNPIs"); obj.add_list("Thresholds", get_thresholds().begin(), get_thresholds().end()); obj.add_element("Duration", get_duration()); - obj.add_element("Interval", get_interval()); + obj.add_element("Delay", get_implementation_delay()); obj.add_element("BaseValue", get_base_value()); + obj.add_element("DirectiveBegin", get_directive_begin()); + obj.add_element("DirectiveEnd", get_directive_end()); } /** @@ -174,27 +224,33 @@ class DynamicNPIs auto obj = io.expect_object("DynamicNPIs"); auto t = obj.expect_list("Thresholds", Tag>>>{}); auto d = obj.expect_element("Duration", Tag>{}); - auto i = obj.expect_element("Interval", Tag>{}); + auto i = obj.expect_element("Delay", Tag>{}); auto b = obj.expect_element("BaseValue", Tag{}); + auto f = obj.expect_element("DirectiveBegin", Tag>{}); + auto l = obj.expect_element("DirectiveEnd", Tag>{}); return apply( io, - [](auto&& t_, auto&& d_, auto&& i_, auto&& b_) { + [](auto&& t_, auto&& d_, auto&& i_, auto&& b_, auto&& f_, auto&& l_) { auto npis = DynamicNPIs(); npis.set_duration(d_); - npis.set_interval(i_); + npis.set_implementation_delay(i_); npis.set_base_value(b_); for (auto&& e : t_) { npis.set_threshold(e.first, e.second); } + npis.set_directive_begin(f_); + npis.set_directive_end(l_); return npis; }, - t, d, i, b); + t, d, i, b, f, l); } private: std::vector>>> m_thresholds; SimulationTime m_duration{14.0}; - SimulationTime m_interval{3.0}; + SimulationTime m_delay{0.0}; + SimulationTime m_directive_begin{SimulationTime(std::numeric_limits::lowest())}; + SimulationTime m_directive_end{SimulationTime(std::numeric_limits::max())}; FP m_base{1.0}; }; @@ -287,18 +343,18 @@ void implement_dynamic_npis(DampingExprGroup& damping_expr_group, const std::vec auto active = get_active_damping(damping_expr, level, type, begin); auto active_end = get_active_damping(damping_expr, level, type, end) - .eval(); //copy because it may be removed or changed + .eval(); // copy because it may be removed or changed auto value = make_matrix(npi.get_value().value() * npi.get_group_weights()); auto npi_implemented = false; - //add begin of npi if not already bigger + // add begin of npi if not already bigger if ((active.array() < value.array()).any()) { damping_expr.add_damping(max(value, active), level, type, begin); npi_implemented = true; } - //replace dampings during the new npi + // replace dampings during the new npi auto damping_indices = get_damping_indices(damping_expr, level, type, begin, end); for (auto& i : damping_indices) { auto& d = damping_expr.get_dampings()[i]; @@ -306,27 +362,27 @@ void implement_dynamic_npis(DampingExprGroup& damping_expr_group, const std::vec npi_implemented = true; } - //add end of npi to restore active dampings if any change was made + // add end of npi to restore active dampings if any change was made if (npi_implemented) { damping_expr.add_damping(active_end, level, type, end); } } } - //remove duplicates that accumulated because of dampings that become active during the time span - //a damping is obsolete if the most recent damping of the same type and level has the same value + // remove duplicates that accumulated because of dampings that become active during the time span + // a damping is obsolete if the most recent damping of the same type and level has the same value for (auto& damping_expr : damping_expr_group) { - //go from the back so indices aren't invalidated when dampings are removed - //use indices to loop instead of reverse iterators because removing invalidates the current iterator + // go from the back so indices aren't invalidated when dampings are removed + // use indices to loop instead of reverse iterators because removing invalidates the current iterator for (auto i = int(0); i < int(damping_expr.get_dampings().size()) - 1; ++i) { auto it = damping_expr.get_dampings().rbegin() + i; - //look for previous damping of the same type/level + // look for previous damping of the same type/level auto it_prev = std::find_if(it + 1, damping_expr.get_dampings().rend(), [&di = *it](auto& dj) { return di.get_level() == dj.get_level() && di.get_type() == dj.get_type(); }); - //remove if match is found and has same value + // remove if match is found and has same value if (it_prev != damping_expr.get_dampings().rend() && it->get_coeffs() == it_prev->get_coeffs()) { damping_expr.remove_damping(damping_expr.get_dampings().size() - 1 - i); } diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 8e44714e78..dadef6739e 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -43,8 +43,7 @@ class SimulationNode using Simulation = Sim; template - requires std::is_constructible_v - SimulationNode(Args&&... args) + requires std::is_constructible_v SimulationNode(Args&&... args) : m_simulation(std::forward(args)...) , m_last_state(m_simulation.get_result().get_last_value()) , m_t0(m_simulation.get_result().get_last_time()) @@ -388,10 +387,10 @@ class MobilityEdge /** @} */ /** - * compute mobility from node_from to node_to. - * mobility is based on coefficients. + * @brief Compute mobility from node_from to node_to. + * Mobility is based on coefficients. * The mobile population is added to the current state of node_to, subtracted from node_from. - * on return, the mobile population (adjusted for infections) is subtracted from node_to, added to node_from. + * On return, the mobile population (adjusted for infections) is subtracted from node_to, added to node_from. * @param t current time * @param dt last time step (fixed to 0.5 for mobility model) * @param node_from node that people changed from, return to @@ -448,7 +447,7 @@ void MobilityEdge::add_mobility_result_time_point(const FP t) } /** - * adjust number of people that changed node when they return according to the model. + * @brief Adjust number of people that changed node when they return according to the model. * E.g. during the time in the other node, some people who left as susceptible will return exposed. * Implemented for general compartmentmodel simulations, overload for your custom model if necessary * so that it can be found with argument-dependent lookup, i.e. in the same namespace as the model. @@ -472,7 +471,7 @@ void calculate_mobility_returns(Eigen::Ref::Vector> mobi } /** - * get the percantage of infected people of the total population in the node + * @brief Get the percentage of infected people of the total population in the node. * If dynamic NPIs are enabled, there needs to be an overload of get_infections_relative(model, y) * for the Model type that can be found with argument-dependent lookup. Ideally define get_infections_relative * in the same namespace as the Model type. @@ -495,7 +494,7 @@ FP get_infections_relative(const SimulationNode& node, FP t, const Eige } /** - * Get an additional mobility factor. + * @brief Get an additional mobility factor. * The absolute mobility for each compartment is computed by c_i * y_i * f_i, wher c_i is the coefficient set in * MobilityParameters, y_i is the current compartment population, f_i is the factor returned by this function. * This factor is optional, default 1.0. If you need to adjust mobility in that way, overload get_mobility_factors(model, t, y) @@ -520,7 +519,7 @@ auto get_mobility_factors(const SimulationNode& node, FP t, const Eigen } /** - * Test persons when moving from their source node. + * @brief Test persons when moving from their source node. * May transfer persons between compartments, e.g., if an infection was detected. * This feature is optional, default implementation does nothing. * In order to support this feature for your model, implement a test_commuters overload @@ -547,14 +546,8 @@ template void mio::MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& node_from, SimulationNode& node_to) { - //check dynamic npis - if (m_t_last_dynamic_npi_check == -std::numeric_limits::infinity()) { - m_t_last_dynamic_npi_check = node_from.get_t0(); - } - auto& dyn_npis = m_parameters.get_dynamic_npis_infected(); - if (dyn_npis.get_thresholds().size() > 0 && - floating_point_greater_equal(t, m_t_last_dynamic_npi_check + dyn_npis.get_interval().get())) { + if (dyn_npis.get_thresholds().size() > 0) { auto inf_rel = get_infections_relative(node_from, t, node_from.get_last_state()) * dyn_npis.get_base_value(); auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); @@ -569,7 +562,6 @@ void mio::MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& m_parameters.get_coefficients().get_shape(), g); }); } - m_t_last_dynamic_npi_check = t; } //returns @@ -651,7 +643,7 @@ void apply_mobility(FP t, FP dt, MobilityEdge& mobilityEdge, SimulationNode< } /** - * create a mobility-based simulation. + * @brief Create a mobility-based simulation. * After every second time step, for each edge a portion of the population corresponding to the coefficients of the edge * changes from one node to the other. In the next timestep, the mobile population returns to their "home" node. * Returns are adjusted based on the development in the target node. @@ -687,7 +679,7 @@ make_mobility_sim(FP t0, FP dt, Graph, MobilityEdge> /** @} */ /** - * Create a graph simulation without mobility. + * @brief Create a graph simulation without mobility. * * Note that we set the time step of the graph simulation to infinity since we do not require any exchange between the * nodes. Hence, in each node, the simulation runs until tmax when advancing the simulation without interruption. diff --git a/cpp/models/ode_secir/analyze_result.h b/cpp/models/ode_secir/analyze_result.h index 9e035d2b4c..3c8acfcc40 100644 --- a/cpp/models/ode_secir/analyze_result.h +++ b/cpp/models/ode_secir/analyze_result.h @@ -63,62 +63,65 @@ std::vector ensemble_params_percentile(const std::vector auto& { - return model.populations[{i, (InfectionState)compart}]; - }); + param_percentil( + node, [ compart, i ](auto&& model) -> auto& { + return model.populations[{i, (InfectionState)compart}]; + }); } // times - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { return model.parameters.template get>()[i]; }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); //probs - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { return model.parameters.template get>()[i]; }); + param_percentil( + node, [i](auto&& model) -> auto& { return model.parameters.template get>()[i]; }); } // group independent params - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); for (size_t run = 0; run < num_runs; run++) { auto const& params = ensemble_params[run][node]; diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index 512d890565..0185dae3d7 100644 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -290,53 +290,71 @@ class Simulation : public BaseT Eigen::Ref> advance(FP tmax) { using std::min; - auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis(); - auto& dyn_npis = this->get_model().parameters.template get>(); - auto& contact_patterns = this->get_model().parameters.template get>(); + auto& dyn_npis = this->get_model().parameters.template get>(); - FP delay_npi_implementation; // delay which is needed to implement a NPI that is criterion-dependent - FP t = BaseT::get_result().get_last_time(); - const FP dt = dyn_npis.get_thresholds().size() > 0 ? FP(dyn_npis.get_interval().get()) : FP(tmax); + // If no dynamic NPI thresholds are set, directly use the base advance + if (dyn_npis.get_thresholds().size() == 0) { + return BaseT::advance(tmax); + } - while (t < tmax) { - FP dt_eff = min(dt, tmax - t); - dt_eff = min(dt_eff, m_t_last_npi_check + dt - t); + auto& contact_patterns = this->get_model().parameters.template get>(); - BaseT::advance(t + dt_eff); + FP delay_npi_implementation; + FP t = BaseT::get_result().get_last_time(); + while (t < tmax) { if (t > 0) { - delay_npi_implementation = - this->get_model().parameters.template get>(); + delay_npi_implementation = FP(dyn_npis.get_implementation_delay()); } - else { // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay. + else { // DynamicNPIs for t=0 are treated as 'from-start NPIs'. I.e., do not enforce delay. delay_npi_implementation = 0; } - t = t + dt_eff; if (dyn_npis.get_thresholds().size() > 0) { - if (floating_point_greater_equal(t, m_t_last_npi_check + dt)) { - if (t < t_end_dyn_npis) { - auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * - dyn_npis.get_base_value(); - auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); - if (exceeded_threshold != dyn_npis.get_thresholds().end() && - (exceeded_threshold->first > m_dynamic_npi.first || - t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired - - auto t_start = SimulationTime(t + delay_npi_implementation); - auto t_end = t_start + SimulationTime(FP(dyn_npis.get_duration())); + FP direc_begin = FP(dyn_npis.get_directive_begin()); + FP direc_end = FP(dyn_npis.get_directive_end()); + if (floating_point_greater_equal(t, direc_begin, 1e-10) && t < direc_end) { + auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * + dyn_npis.get_base_value(); + auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); + if (exceeded_threshold != dyn_npis.get_thresholds().end() && + (exceeded_threshold->first > m_dynamic_npi.first || + t > FP(m_dynamic_npi.second))) { // old npi was weaker or is expired + + // Keep-alive: if the NPI expired but the threshold is still exceeded at the + // same level, renew immediately without delay to avoid a gap. + // Apply implementation delay only if stronger NPI needed. + bool is_expiry_renewal = + (t > FP(m_dynamic_npi.second)) && !(exceeded_threshold->first > m_dynamic_npi.first); + FP effective_delay = is_expiry_renewal ? FP(0) : delay_npi_implementation; + if (t + effective_delay < direc_end) { + auto t_start = SimulationTime(t + effective_delay); + // set the end to the minimum of start+duration and the end of the directive + auto t_end = SimulationTime(min(direc_end, FP(t_start + dyn_npis.get_duration()))); m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end); + // For new NPIs (t_start > 0): shift t_start_damping by +1 so the smooth + // transition window [t_start, t_start+1] lies in the future, consistent with + // predefined dampings. For t_start = 0 (global t0): the window [-1, 0] + // is in the past and no shift is needed. + // For keep-alive renewals (is_expiry_renewal): do not shift t_start_damping. + // Since t_start == t_end_damping_old, the new start damping has the same + // (time, level, type) as the previous entry. + // Therefore, the contact matrix stays constant at the NPI level with no dip. + auto t_start_damping = (!is_expiry_renewal && FP(t_start) > FP(0)) + ? SimulationTime(FP(t_start) + FP(1)) + : t_start; + auto t_end_damping = (FP(t_start) > FP(0)) ? SimulationTime(FP(t_end) + FP(1)) : t_end; implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second, - t_start, t_end, [](auto& g) { + t_start_damping, t_end_damping, [](auto& g) { return make_contact_damping_matrix(g); }); } } - m_t_last_npi_check = t; } } - else { - m_t_last_npi_check = t; - } + + FP dt_eff = min(1.0, tmax - t); + BaseT::advance(t + dt_eff); + t = t + dt_eff; } return this->get_result().get_last_value(); diff --git a/cpp/models/ode_secir/parameter_space.h b/cpp/models/ode_secir/parameter_space.h index 2ac6ddab5c..bbaf8ba657 100644 --- a/cpp/models/ode_secir/parameter_space.h +++ b/cpp/models/ode_secir/parameter_space.h @@ -70,7 +70,6 @@ void set_params_distributions_normal(Model& model, FP t0, FP tmax, FP dev_re set_distribution(model.parameters.template get>(), 0.0); set_distribution(model.parameters.template get>()); set_distribution(model.parameters.template get>()); - set_distribution(model.parameters.template get>()); // populations for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) { @@ -216,8 +215,6 @@ Graph, MobilityParameters> draw_sample(Graph, MobilityPa shared_contacts.draw_sample_dampings(); auto& shared_dynamic_npis = shared_params_model.parameters.template get>(); shared_dynamic_npis.draw_sample(); - auto& shared_dynamic_npis_delay = shared_params_model.parameters.template get>(); - shared_dynamic_npis_delay.draw_sample(); for (auto& params_node : graph.nodes()) { auto& node_model = params_node.property; diff --git a/cpp/models/ode_secir/parameters.h b/cpp/models/ode_secir/parameters.h index 0de7e86b1d..fe73d87e39 100644 --- a/cpp/models/ode_secir/parameters.h +++ b/cpp/models/ode_secir/parameters.h @@ -351,22 +351,6 @@ struct DynamicNPIsInfectedSymptoms { } }; -/** - * @brief The delay with which DynamicNPIs are implemented and enforced after exceedance of threshold. - */ -template -struct DynamicNPIsImplementationDelay { - using Type = UncertainValue; - static Type get_default(AgeGroup /*size*/) - { - return Type(0.0); - } - static std::string name() - { - return "DynamicNPIsImplementationDelay"; - } -}; - /** * @brief capacity to test and trace contacts of infected for quarantine per day. */ @@ -402,12 +386,12 @@ struct TestAndTraceCapacityMaxRisk { template using ParametersBase = ParameterSet, Seasonality, ICUCapacity, TestAndTraceCapacity, - TestAndTraceCapacityMaxRisk, ContactPatterns, DynamicNPIsImplementationDelay, - DynamicNPIsInfectedSymptoms, TimeExposed, TimeInfectedNoSymptoms, TimeInfectedSymptoms, - TimeInfectedSevere, TimeInfectedCritical, TransmissionProbabilityOnContact, - RelativeTransmissionNoSymptoms, RecoveredPerInfectedNoSymptoms, - RiskOfInfectionFromSymptomatic, MaxRiskOfInfectionFromSymptomatic, - SeverePerInfectedSymptoms, CriticalPerSevere, DeathsPerSevere, DeathsPerCritical>; + TestAndTraceCapacityMaxRisk, ContactPatterns, DynamicNPIsInfectedSymptoms, TimeExposed, + TimeInfectedNoSymptoms, TimeInfectedSymptoms, TimeInfectedSevere, TimeInfectedCritical, + TransmissionProbabilityOnContact, RelativeTransmissionNoSymptoms, + RecoveredPerInfectedNoSymptoms, RiskOfInfectionFromSymptomatic, + MaxRiskOfInfectionFromSymptomatic, SeverePerInfectedSymptoms, CriticalPerSevere, + DeathsPerSevere, DeathsPerCritical>; /** * @brief Parameters of an age-resolved SECIR/SECIHURD model. @@ -465,18 +449,6 @@ class Parameters : public ParametersBase return m_end_commuter_detection; } - /** - * Time in simulation after which no dynamic NPIs are applied. - */ - FP& get_end_dynamic_npis() - { - return m_end_dynamic_npis; - } - FP get_end_dynamic_npis() const - { - return m_end_dynamic_npis; - } - /** * @brief Checks whether all Parameters satisfy their corresponding constraints and applies them, if they do not. * Time spans cannot be negative and probabilities can only take values between [0,1]. @@ -509,13 +481,6 @@ class Parameters : public ParametersBase corrected = true; } - if (this->template get>() < 0.0) { - log_warning("Constraint check: Parameter DynamicNPIsImplementationDelay changed from {} to {}", - this->template get>(), 0); - this->template set>(0); - corrected = true; - } - if (this->template get>() < 0.0) { log_warning("Constraint check: Parameter TestAndTraceCapacity changed from {} to {}", this->template get>(), 0); @@ -677,11 +642,6 @@ class Parameters : public ParametersBase return true; } - if (this->template get>() < 0.0) { - log_error("Constraint check: Parameter DynamicNPIsImplementationDelay smaller {}", 0); - return true; - } - const FP tol_times = 1e-1; // accepted tolerance for compartment stays for (auto i = AgeGroup(0); i < AgeGroup(m_num_groups); ++i) { @@ -805,7 +765,6 @@ class Parameters : public ParametersBase FP m_commuter_nondetection = 0.0; FP m_start_commuter_detection = 0.0; FP m_end_commuter_detection = 0.0; - FP m_end_dynamic_npis = std::numeric_limits::max(); }; /** diff --git a/cpp/models/ode_secirts/analyze_result.h b/cpp/models/ode_secirts/analyze_result.h index 4970ae932f..496e51e014 100644 --- a/cpp/models/ode_secirts/analyze_result.h +++ b/cpp/models/ode_secirts/analyze_result.h @@ -72,109 +72,124 @@ std::vector ensemble_params_percentile(const std::vector(0); compart < InfectionState::Count; ++compart) { - param_percentil(node, [compart, i](auto&& model) -> auto& { - return model.populations[{i, compart}]; - }); + param_percentil( + node, [ compart, i ](auto&& model) -> auto& { + return model.populations[{i, compart}]; + }); } // times - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { return model.parameters.template get>()[i]; }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); //probs - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { return model.parameters.template get>()[i]; }); + param_percentil( + node, [i](auto&& model) -> auto& { return model.parameters.template get>()[i]; }); //vaccinations - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - - for (auto day = SimulationDay(0); day < num_days; ++day) { - param_percentil(node, [i, day](auto&& model) -> auto& { - return model.parameters.template get>()[{i, day}]; + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; }); - param_percentil(node, [i, day](auto&& model) -> auto& { - return model.parameters.template get>()[{i, day}]; + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; }); - param_percentil(node, [i, day](auto&& model) -> auto& { - return model.parameters.template get>()[{i, day}]; + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + + for (auto day = SimulationDay(0); day < num_days; ++day) { + param_percentil( + node, [ i, day ](auto&& model) -> auto& { + return model.parameters.template get>()[{i, day}]; + }); + param_percentil( + node, [ i, day ](auto&& model) -> auto& { + return model.parameters.template get>()[{i, day}]; + }); + param_percentil( + node, [ i, day ](auto&& model) -> auto& { + return model.parameters.template get>()[{i, day}]; + }); } //virus variants - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); } // group independent params - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); for (size_t run = 0; run < num_runs; run++) { diff --git a/cpp/models/ode_secirts/model.h b/cpp/models/ode_secirts/model.h index 29ad4a2ce7..b58bb28f4d 100644 --- a/cpp/models/ode_secirts/model.h +++ b/cpp/models/ode_secirts/model.h @@ -759,7 +759,6 @@ class Simulation : public BaseT using std::floor; using std::min; - auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis(); auto& dyn_npis = this->get_model().parameters.template get>(); auto& contact_patterns = this->get_model().parameters.template get>(); // const size_t num_groups = (size_t)this->get_model().parameters.get_num_groups(); @@ -769,57 +768,68 @@ class Simulation : public BaseT auto base_infectiousness = this->get_model().parameters.template get>(); FP delay_npi_implementation; - auto t = BaseT::get_result().get_last_time(); - const auto dt = dyn_npis.get_interval().get(); + FP t = BaseT::get_result().get_last_time(); while (t < tmax) { - auto dt_eff = min({dt, tmax - t, m_t_last_npi_check + dt - t}); - if (dt_eff >= 1.0) { - dt_eff = 1.0; - } - - BaseT::advance(t + dt_eff); - if (t + 0.5 + dt_eff - floor(t + 0.5) >= 1) { - this->apply_variant(t, base_infectiousness); - } - if (t > 0) { - delay_npi_implementation = - this->get_model().parameters.template get>(); + delay_npi_implementation = FP(dyn_npis.get_implementation_delay()); } else { - // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay. + // DynamicNPIs for t=0 are treated as 'from-start NPIs'. I.e., do not enforce delay. delay_npi_implementation = 0; } - t = t + dt_eff; if (dyn_npis.get_thresholds().size() > 0) { - if (floating_point_greater_equal(t, m_t_last_npi_check + dt)) { - if (t < t_end_dyn_npis) { - auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * - dyn_npis.get_base_value(); - auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); - if (exceeded_threshold != dyn_npis.get_thresholds().end() && - (exceeded_threshold->first > m_dynamic_npi.first || - t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired - - auto t_start = SimulationTime(t + delay_npi_implementation); - auto t_end = t_start + SimulationTime(dyn_npis.get_duration()); + FP direc_begin = FP(dyn_npis.get_directive_begin()); + FP direc_end = FP(dyn_npis.get_directive_end()); + if (floating_point_greater_equal(t, direc_begin, 1e-10) && t < direc_end) { + auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * + dyn_npis.get_base_value(); + auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); + if (exceeded_threshold != dyn_npis.get_thresholds().end() && + (exceeded_threshold->first > m_dynamic_npi.first || + t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired + + // Keep-alive: if the NPI expired but the threshold is still exceeded at the + // same level, renew immediately without delay to avoid a gap. + // Apply implementation delay only if stronger NPI needed. + bool is_expiry_renewal = + (t > FP(m_dynamic_npi.second)) && !(exceeded_threshold->first > m_dynamic_npi.first); + FP effective_delay = is_expiry_renewal ? FP(0) : delay_npi_implementation; + if (t + effective_delay < direc_end) { + auto t_start = SimulationTime(t + effective_delay); + // set the end to the minimum of start+duration and the end of the directive + auto t_end = SimulationTime(min(direc_end, FP(t_start + dyn_npis.get_duration()))); this->get_model().parameters.get_start_commuter_detection() = (FP)t_start; this->get_model().parameters.get_end_commuter_detection() = (FP)t_end; m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end); + // For new NPIs (t_start > 0): shift t_start_damping by +1 so the smooth + // transition window [t_start, t_start+1] lies in the future, consistent with + // predefined dampings. For t_start = 0 (global t0): the window [-1, 0] + // is in the past and no shift is needed. + // For keep-alive renewals (is_expiry_renewal): do not shift t_start_damping. + // Since t_start == t_end_damping_old, the new start damping has the same + // (time, level, type) as the previous entry. + // Therefore, the contact matrix stays constant at the NPI level with no dip. + auto t_start_damping = (!is_expiry_renewal && FP(t_start) > FP(0)) + ? SimulationTime(FP(t_start) + FP(1)) + : t_start; + auto t_end_damping = (FP(t_start) > FP(0)) ? SimulationTime(FP(t_end) + FP(1)) : t_end; implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second, - t_start, t_end, [](auto& g) { + t_start_damping, t_end_damping, [](auto& g) { return make_contact_damping_matrix(g); }); } } - m_t_last_npi_check = t; } } - else { - m_t_last_npi_check = t; + + auto dt_eff = min(1.0, tmax - t); + BaseT::advance(t + dt_eff); + if (t + 0.5 + dt_eff - floor(t + 0.5) >= 1) { + this->apply_variant(t, base_infectiousness); } + t = t + dt_eff; } // reset TransmissionProbabilityOnContact. This is important for the graph simulation where the advance // function is called multiple times for the same model. diff --git a/cpp/models/ode_secirts/parameters.h b/cpp/models/ode_secirts/parameters.h index 370f1ed602..df26815127 100644 --- a/cpp/models/ode_secirts/parameters.h +++ b/cpp/models/ode_secirts/parameters.h @@ -746,22 +746,6 @@ struct InfectiousnessNewVariant { } }; -/** - * @brief The delay with which DynamicNPIs are implemented and enforced after exceedance of threshold. - */ -template -struct DynamicNPIsImplementationDelay { - using Type = UncertainValue; - static Type get_default(AgeGroup /*size*/) - { - return Type(0.0); - } - static std::string name() - { - return "DynamicNPIsImplementationDelay"; - } -}; - template using ParametersBase = ParameterSet< StartDay, Seasonality, ICUCapacity, TestAndTraceCapacity, TestAndTraceCapacityMaxRiskNoSymptoms, @@ -776,8 +760,7 @@ using ParametersBase = ParameterSet< DailyBoosterVaccinations, ReducExposedPartialImmunity, ReducExposedImprovedImmunity, ReducInfectedSymptomsPartialImmunity, ReducInfectedSymptomsImprovedImmunity, ReducInfectedSevereCriticalDeadPartialImmunity, ReducInfectedSevereCriticalDeadImprovedImmunity, - ReducTimeInfectedMild, InfectiousnessNewVariant, DynamicNPIsImplementationDelay, - StartDayNewVariant>; + ReducTimeInfectedMild, InfectiousnessNewVariant, StartDayNewVariant>; /** * @brief Parameters of the age-resolved SECIRS-type model with high temporary immunity upon immunization and waning immunity over @@ -836,18 +819,6 @@ class Parameters : public ParametersBase return m_end_commuter_detection; } - /** - * Time in simulation after which no dynamic NPIs are applied. - */ - FP& get_end_dynamic_npis() - { - return m_end_dynamic_npis; - } - FP get_end_dynamic_npis() const - { - return m_end_dynamic_npis; - } - /** * @brief Checks whether all Parameters satisfy their corresponding constraints and applies them, if they do not. * Time spans cannot be negative and probabilities can only take values between [0,1]. @@ -899,13 +870,6 @@ class Parameters : public ParametersBase corrected = true; } - if (this->template get>() < 0.0) { - log_warning("Constraint check: Parameter DynamicNPIsImplementationDelay changed from {} to {}", - this->template get>(), 0); - this->template set>(0); - corrected = true; - } - const FP tol_times = 1e-1; // accepted tolerance for compartment stays for (auto i = AgeGroup(0); i < AgeGroup(m_num_groups); ++i) { @@ -1177,11 +1141,6 @@ class Parameters : public ParametersBase return true; } - if (this->template get>() < 0.0) { - log_error("Constraint check: Parameter DynamicNPIsImplementationDelay smaller {}", 0); - return true; - } - for (auto i = AgeGroup(0); i < AgeGroup(m_num_groups); ++i) { if (this->template get>()[i] < tol_times) { @@ -1397,7 +1356,6 @@ class Parameters : public ParametersBase FP m_commuter_nondetection = 0.0; FP m_start_commuter_detection = 0.0; FP m_end_commuter_detection = 0.0; - FP m_end_dynamic_npis = std::numeric_limits::max(); }; } // namespace osecirts diff --git a/cpp/models/ode_secirvvs/analyze_result.h b/cpp/models/ode_secirvvs/analyze_result.h index 8c9cab6a42..2544bb9a98 100644 --- a/cpp/models/ode_secirvvs/analyze_result.h +++ b/cpp/models/ode_secirvvs/analyze_result.h @@ -67,109 +67,128 @@ std::vector ensemble_params_percentile(const std::vector(0); compart < InfectionState::Count; ++compart) { - param_percentil(node, [compart, i](auto&& model) -> auto& { - return model.populations[{i, compart}]; - }); + param_percentil( + node, [ compart, i ](auto&& model) -> auto& { + return model.populations[{i, compart}]; + }); } // times - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { return model.parameters.template get>()[i]; }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); //probs - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); //vaccinations - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); - - for (auto day = SimulationDay(0); day < num_days; ++day) { - param_percentil(node, [i, day](auto&& model) -> auto& { - return model.parameters.template get>()[{i, day}]; + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; }); - param_percentil(node, [i, day](auto&& model) -> auto& { - return model.parameters.template get>()[{i, day}]; + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); + + for (auto day = SimulationDay(0); day < num_days; ++day) { + param_percentil( + node, [ i, day ](auto&& model) -> auto& { + return model.parameters.template get>()[{i, day}]; + }); + param_percentil( + node, [ i, day ](auto&& model) -> auto& { + return model.parameters.template get>()[{i, day}]; + }); } //virus variants - param_percentil(node, [i](auto&& model) -> auto& { - return model.parameters.template get>()[i]; - }); + param_percentil( + node, [i](auto&& model) -> auto& { + return model.parameters.template get>()[i]; + }); } // group independent params - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); - param_percentil(node, [](auto&& model) -> auto& { - return model.parameters.template get>(); - }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); + param_percentil( + node, [](auto&& model) -> auto& { return model.parameters.template get>(); }); for (size_t run = 0; run < num_runs; run++) { auto const& params = ensemble_params[run][node]; diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 9a9c842d2e..5184846c88 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -681,7 +681,6 @@ class Simulation : public BaseT using std::floor; using std::min; - auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis(); auto& dyn_npis = this->get_model().parameters.template get>(); auto& contact_patterns = this->get_model().parameters.template get>(); // const size_t num_groups = (size_t)this->get_model().parameters.get_num_groups(); @@ -691,63 +690,72 @@ class Simulation : public BaseT auto base_infectiousness = this->get_model().parameters.template get>(); FP delay_npi_implementation; - FP t = BaseT::get_result().get_last_time(); - const FP dt = dyn_npis.get_thresholds().size() > 0 ? dyn_npis.get_interval().get() : tmax; + FP t = BaseT::get_result().get_last_time(); while (t < tmax) { - FP dt_eff = min(dt, tmax - t); - dt_eff = min(dt_eff, m_t_last_npi_check + dt - t); - - if (dt_eff >= 1.0) { - dt_eff = 1.0; + if (t > 0) { + delay_npi_implementation = FP(dyn_npis.get_implementation_delay()); + } + else { // DynamicNPIs for t=0 are treated as 'from-start NPIs'. I.e., do not enforce delay. + delay_npi_implementation = 0; } - if (t == 0) { //this->apply_vaccination(t); // done in init now? this->apply_variant(t, base_infectiousness); } - BaseT::advance(t + dt_eff); - if (t + 0.5 + dt_eff - floor(t + 0.5) >= 1) { - this->apply_vaccination(t + 0.5 + dt_eff); - this->apply_variant(t, base_infectiousness); - } - - if (t > 0) { - delay_npi_implementation = - this->get_model().parameters.template get>(); - } - else { // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay. - delay_npi_implementation = 0; - } - t = t + dt_eff; if (dyn_npis.get_thresholds().size() > 0) { - if (floating_point_greater_equal(t, m_t_last_npi_check + dt)) { - if (t < t_end_dyn_npis) { - auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * - dyn_npis.get_base_value(); - auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); - if (exceeded_threshold != dyn_npis.get_thresholds().end() && - (exceeded_threshold->first > m_dynamic_npi.first || - t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired - - auto t_start = SimulationTime(t + delay_npi_implementation); - auto t_end = t_start + SimulationTime(dyn_npis.get_duration()); - this->get_model().parameters.get_start_commuter_detection() = t_start.get(); - this->get_model().parameters.get_end_commuter_detection() = t_end.get(); + FP direc_begin = FP(dyn_npis.get_directive_begin()); + FP direc_end = FP(dyn_npis.get_directive_end()); + if (floating_point_greater_equal(t, direc_begin, 1e-10) && t < direc_end) { + auto inf_rel = get_infections_relative(*this, t, this->get_result().get_last_value()) * + dyn_npis.get_base_value(); + auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel); + if (exceeded_threshold != dyn_npis.get_thresholds().end() && + (exceeded_threshold->first > m_dynamic_npi.first || + t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired + + // Keep-alive: if the NPI expired but the threshold is still exceeded at the + // same level, renew immediately without delay to avoid a gap. + // Apply implementation delay only if stronger NPI needed. + bool is_expiry_renewal = + (t > FP(m_dynamic_npi.second)) && !(exceeded_threshold->first > m_dynamic_npi.first); + FP effective_delay = is_expiry_renewal ? FP(0) : delay_npi_implementation; + if (t + effective_delay < direc_end) { + auto t_start = SimulationTime(t + effective_delay); + // set the end to the minimum of start+duration and the end of the directive + auto t_end = SimulationTime(min(direc_end, FP(t_start + dyn_npis.get_duration()))); + this->get_model().parameters.get_start_commuter_detection() = (FP)t_start; + this->get_model().parameters.get_end_commuter_detection() = (FP)t_end; m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end); + // For new NPIs (t_start > 0): shift t_start_damping by +1 so the smooth + // transition window [t_start, t_start+1] lies in the future, consistent with + // predefined dampings. For t_start = 0 (global t0): the window [-1, 0] + // is in the past and no shift is needed. + // For keep-alive renewals (is_expiry_renewal): do not shift t_start_damping. + // Since t_start == t_end_damping_old, the new start damping has the same + // (time, level, type) as the previous entry. + // Therefore, the contact matrix stays constant at the NPI level with no dip. + auto t_start_damping = (!is_expiry_renewal && FP(t_start) > FP(0)) + ? SimulationTime(FP(t_start) + FP(1)) + : t_start; + auto t_end_damping = (FP(t_start) > FP(0)) ? SimulationTime(FP(t_end) + FP(1)) : t_end; implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second, - t_start, t_end, [](auto& g) { + t_start_damping, t_end_damping, [](auto& g) { return make_contact_damping_matrix(g); }); } } - m_t_last_npi_check = t; } } - else { - m_t_last_npi_check = t; + + FP dt_eff = min(1.0, tmax - t); + BaseT::advance(t + dt_eff); + if (t + 0.5 + dt_eff - floor(t + 0.5) >= 1) { + this->apply_vaccination(t + 0.5 + dt_eff); + this->apply_variant(t, base_infectiousness); } + t = t + dt_eff; } // reset TransmissionProbabilityOnContact. This is important for the graph simulation where the advance // function is called multiple times for the same model. diff --git a/cpp/models/ode_secirvvs/parameter_space.h b/cpp/models/ode_secirvvs/parameter_space.h index 875c08f5ef..a9078fadcb 100644 --- a/cpp/models/ode_secirvvs/parameter_space.h +++ b/cpp/models/ode_secirvvs/parameter_space.h @@ -176,8 +176,6 @@ Graph, MobilityParameters> draw_sample(Graph, MobilityPa shared_contacts.draw_sample_dampings(); auto& shared_dynamic_npis = shared_params_model.parameters.template get>(); shared_dynamic_npis.draw_sample(); - auto& shared_dynamic_npis_delay = shared_params_model.parameters.template get>(); - shared_dynamic_npis_delay.draw_sample(); for (auto& params_node : graph.nodes()) { auto& node_model = params_node.property; diff --git a/cpp/models/ode_secirvvs/parameters.h b/cpp/models/ode_secirvvs/parameters.h index 3a04a3ca40..bc724c52c1 100644 --- a/cpp/models/ode_secirvvs/parameters.h +++ b/cpp/models/ode_secirvvs/parameters.h @@ -185,22 +185,6 @@ struct DynamicNPIsInfectedSymptoms { } }; -/** - * @brief The delay with which DynamicNPIs are implemented and enforced after exceedance of threshold. - */ -template -struct DynamicNPIsImplementationDelay { - using Type = UncertainValue; - static Type get_default(AgeGroup /*size*/) - { - return Type(0.0); - } - static std::string name() - { - return "DynamicNPIsImplementationDelay"; - } -}; - /** * @brief the (mean) latent time in day unit */ @@ -643,17 +627,16 @@ struct InfectiousnessNewVariant { template using ParametersBase = ParameterSet< StartDay, Seasonality, ICUCapacity, TestAndTraceCapacity, TestAndTraceCapacityMaxRiskNoSymptoms, - TestAndTraceCapacityMaxRiskSymptoms, ContactPatterns, DynamicNPIsImplementationDelay, - DynamicNPIsInfectedSymptoms, TimeExposed, TimeInfectedNoSymptoms, TimeInfectedSymptoms, - TimeInfectedSevere, TimeInfectedCritical, TransmissionProbabilityOnContact, - RelativeTransmissionNoSymptoms, RecoveredPerInfectedNoSymptoms, RiskOfInfectionFromSymptomatic, - MaxRiskOfInfectionFromSymptomatic, SeverePerInfectedSymptoms, CriticalPerSevere, DeathsPerSevere, - DeathsPerCritical, VaccinationGap, DaysUntilEffectivePartialImmunity, - DaysUntilEffectiveImprovedImmunity, DailyFullVaccinations, DailyPartialVaccinations, - ReducExposedPartialImmunity, ReducExposedImprovedImmunity, ReducInfectedSymptomsPartialImmunity, - ReducInfectedSymptomsImprovedImmunity, ReducInfectedSevereCriticalDeadPartialImmunity, - ReducInfectedSevereCriticalDeadImprovedImmunity, ReducTimeInfectedMild, InfectiousnessNewVariant, - StartDayNewVariant>; + TestAndTraceCapacityMaxRiskSymptoms, ContactPatterns, DynamicNPIsInfectedSymptoms, TimeExposed, + TimeInfectedNoSymptoms, TimeInfectedSymptoms, TimeInfectedSevere, TimeInfectedCritical, + TransmissionProbabilityOnContact, RelativeTransmissionNoSymptoms, RecoveredPerInfectedNoSymptoms, + RiskOfInfectionFromSymptomatic, MaxRiskOfInfectionFromSymptomatic, SeverePerInfectedSymptoms, + CriticalPerSevere, DeathsPerSevere, DeathsPerCritical, VaccinationGap, + DaysUntilEffectivePartialImmunity, DaysUntilEffectiveImprovedImmunity, DailyFullVaccinations, + DailyPartialVaccinations, ReducExposedPartialImmunity, ReducExposedImprovedImmunity, + ReducInfectedSymptomsPartialImmunity, ReducInfectedSymptomsImprovedImmunity, + ReducInfectedSevereCriticalDeadPartialImmunity, ReducInfectedSevereCriticalDeadImprovedImmunity, + ReducTimeInfectedMild, InfectiousnessNewVariant, StartDayNewVariant>; /** * @brief Parameters of an age-resolved SECIR/SECIHURD model with paths for partial and improved immunity through vaccination. @@ -711,18 +694,6 @@ class Parameters : public ParametersBase return m_end_commuter_detection; } - /** - * Time in simulation after which no dynamic NPIs are applied. - */ - FP& get_end_dynamic_npis() - { - return m_end_dynamic_npis; - } - FP get_end_dynamic_npis() const - { - return m_end_dynamic_npis; - } - /** * @brief Checks whether all Parameters satisfy their corresponding constraints and applies them, if they do not. * Time spans cannot be negative and probabilities can only take values between [0,1]. @@ -774,13 +745,6 @@ class Parameters : public ParametersBase corrected = true; } - if (this->template get>() < 0.0) { - log_warning("Constraint check: Parameter DynamicNPIsImplementationDelay changed from {} to {}", - this->template get>(), 0); - this->template set>(0); - corrected = true; - } - const FP tol_times = 1e-1; // accepted tolerance for compartment stays for (auto i = AgeGroup(0); i < AgeGroup(m_num_groups); ++i) { @@ -1013,11 +977,6 @@ class Parameters : public ParametersBase return true; } - if (this->template get>() < 0.0) { - log_error("Constraint check: Parameter DynamicNPIsImplementationDelay smaller {}", 0); - return true; - } - for (auto i = AgeGroup(0); i < AgeGroup(m_num_groups); ++i) { if (this->template get>()[i] < tol_times) { @@ -1200,7 +1159,6 @@ class Parameters : public ParametersBase FP m_commuter_nondetection = 0.0; FP m_start_commuter_detection = 0.0; FP m_end_commuter_detection = 0.0; - FP m_end_dynamic_npis = std::numeric_limits::max(); }; } // namespace osecirvvs diff --git a/cpp/tests/data/results_osecirvvs.h5 b/cpp/tests/data/results_osecirvvs.h5 index 0e3e79575e..34707a720b 100644 Binary files a/cpp/tests/data/results_osecirvvs.h5 and b/cpp/tests/data/results_osecirvvs.h5 differ diff --git a/cpp/tests/test_dynamic_npis.cpp b/cpp/tests/test_dynamic_npis.cpp index b545878abc..0b494999a7 100644 --- a/cpp/tests/test_dynamic_npis.cpp +++ b/cpp/tests/test_dynamic_npis.cpp @@ -22,6 +22,7 @@ #include "memilio/utils/compiler_diagnostics.h" #include "ode_secir/model.h" #include "ode_secirvvs/model.h" +#include "ode_secirts/model.h" #include "matchers.h" #include @@ -63,6 +64,25 @@ TEST(DynamicNPIs, get_threshold) EXPECT_EQ(npis.get_max_exceeded_threshold(0.5), npis.get_thresholds().end()); } +TEST(DynamicNPIs, get_and_set_delay) +{ + mio::DynamicNPIs npis; + EXPECT_EQ(npis.get_implementation_delay(), mio::SimulationTime(0.0)); + npis.set_implementation_delay(mio::SimulationTime(7.0)); + EXPECT_EQ(npis.get_implementation_delay(), mio::SimulationTime(7.0)); +} + +TEST(DynamicNPIs, get_and_set_directive_begin_end) +{ + mio::DynamicNPIs npis; + EXPECT_LE(npis.get_directive_begin(), mio::SimulationTime(-1e4)); + EXPECT_GE(npis.get_directive_end(), mio::SimulationTime(1e4)); + npis.set_directive_begin(mio::SimulationTime(7.0)); + EXPECT_EQ(npis.get_directive_begin(), mio::SimulationTime(7.0)); + npis.set_directive_end(mio::SimulationTime(14.0)); + EXPECT_EQ(npis.get_directive_end(), mio::SimulationTime(14.0)); +} + TEST(DynamicNPIs, get_damping_indices) { using Damping = mio::Damping>; @@ -315,7 +335,6 @@ TEST(DynamicNPIs, mobility) Eigen::VectorXd::Ones(2)}}); npis.set_duration(mio::SimulationTime(5.0)); npis.set_base_value(100'000); - npis.set_interval(mio::SimulationTime(3.0)); mio::MobilityCoefficientGroup coeffs(1, 2); mio::MobilityParameters parameters(coeffs); @@ -323,11 +342,11 @@ TEST(DynamicNPIs, mobility) mio::MobilityEdge edge(parameters); - EXPECT_EQ(edge.get_parameters().get_coefficients()[0].get_dampings().size(), 0); //initial + EXPECT_EQ(edge.get_parameters().get_coefficients()[0].get_dampings().size(), 0); // initial edge.apply_mobility(0.5, 0.5, node_from, node_to); - EXPECT_EQ(edge.get_parameters().get_coefficients()[0].get_dampings().size(), 0); //not check at the beginning + EXPECT_EQ(edge.get_parameters().get_coefficients()[0].get_dampings().size(), 0); // no check at the beginning EXPECT_CALL(node_from.get_simulation(), advance).Times(1).WillOnce([&](auto t) { node_from.get_simulation().result.add_time_point(t, last_state_safe); @@ -336,7 +355,7 @@ TEST(DynamicNPIs, mobility) node_to.advance(3.0, 2.5); edge.apply_mobility(3.0, 2.5, node_from, node_to); - EXPECT_EQ(edge.get_parameters().get_coefficients()[0].get_dampings().size(), 0); //threshold not exceeded + EXPECT_EQ(edge.get_parameters().get_coefficients()[0].get_dampings().size(), 0); // threshold not exceeded EXPECT_CALL(node_from.get_simulation(), advance).Times(1).WillOnce([&](auto t) { node_from.get_simulation().result.add_time_point(t, last_state_crit); @@ -346,7 +365,7 @@ TEST(DynamicNPIs, mobility) edge.apply_mobility(4.5, 1.5, node_from, node_to); EXPECT_EQ(edge.get_parameters().get_coefficients()[0].get_dampings().size(), - 0); //threshold exceeded, but only check every 3 days + 2); // threshold exceeded, NPI immediately applied EXPECT_CALL(node_from.get_simulation(), advance).Times(1).WillOnce([&](auto t) { node_from.get_simulation().result.add_time_point(t, last_state_crit); @@ -445,8 +464,8 @@ TEST(DynamicNPIs, secir_threshold_exceeded) Eigen::VectorXd::Ones(1)}}); npis.set_duration(mio::SimulationTime(5.0)); npis.set_base_value(50'000); - model.parameters.get>() = npis; - model.parameters.get>() = 0.0; + npis.set_implementation_delay(mio::SimulationTime(0.0)); + model.parameters.get>() = npis; EXPECT_EQ(model.parameters.get>().get_cont_freq_mat()[0].get_dampings().size(), 0); @@ -485,33 +504,225 @@ TEST(DynamicNPIs, secir_delayed_implementation) EXPECT_EQ(model.parameters.get>().get_cont_freq_mat()[0].get_dampings().size(), 0); - // start with t0 = 0.0 - model.parameters.get>() = 3.0; + // start with t0 = 0.0 so dynamicNPIs are active from the start + // t0=0 is special: delay is not enforced, t_start=0, smooth window [-1, 0] -> full effect at t=0 + npis.set_implementation_delay(mio::SimulationTime(3.0)); + model.parameters.get>() = npis; mio::osecir::Simulation> sim(model, 0.0); sim.advance(3.0); mio::ContactMatrixGroup const& contact_matrix = sim.get_model().parameters.template get>(); - EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(2.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 0.5); EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 0.5); - // second simulation with t0 = 1.0, so the NPIs are implemented at tmax + delay = 6.0 - const auto tmax = 4.0; - model.parameters.get>() = 2.0; + // second simulation with t0 = 1.0: NPI check at t=1, delay=2 -> t_start=3, damping at t_start+1=4 + // smooth window [3, 4]: at t=3 still 1.0, at t=4 full effect 0.5 + npis.set_implementation_delay(mio::SimulationTime(2.0)); + model.parameters.get>() = npis; mio::osecir::Simulation> sim_2(model, 1.0); - sim_2.advance(tmax); + sim_2.advance(5.0); mio::ContactMatrixGroup const& contact_matrix_sim_2 = sim_2.get_model().parameters.template get>(); - EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(5.0))(0, 0), 1.0); - EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(6.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(2.0))(0, 0), 1.0); // before ramp + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); // full effect - // third simulation with t0 = 1.0, so the NPIs are implemented at tmax + delay = 14.0 - model.parameters.get>() = 10.0; + // third simulation: delay=10 -> t_start=11, damping at t_start+1=12 + // smooth window [11, 12]: at t=11 still 1.0, at t=12 full effect 0.5 + npis.set_implementation_delay(mio::SimulationTime(10.0)); + model.parameters.get>() = npis; mio::osecir::Simulation> sim_3(model, 1.0); sim_3.advance(4.0); mio::ContactMatrixGroup const& contact_matrix_sim_3 = sim_3.get_model().parameters.template get>(); - EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(13.0))(0, 0), 1.0); - EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(14.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(11.0))(0, 0), 1.0); // before ramp + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(12.0))(0, 0), 0.5); // full effect +} + +TEST(DynamicNPIs, secir_implementation_with_directives) +{ + mio::osecir::Model model(1); + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = 10; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, 100); + + mio::ContactMatrixGroup& cm = model.parameters.get>(); + cm[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1.0)); + + mio::DynamicNPIs npis; + npis.set_threshold(0.05 * 50'000, {mio::DampingSampling{0.5, + mio::DampingLevel(0), + mio::DampingType(0), + mio::SimulationTime(0), + {0}, + Eigen::VectorXd::Ones(1)}}); + npis.set_duration(mio::SimulationTime(5.0)); + npis.set_base_value(50'000); + model.parameters.get>() = npis; + + EXPECT_EQ(model.parameters.get>().get_cont_freq_mat()[0].get_dampings().size(), + 0); + + // directive begin is after the simulation, so no NPI is implemented + npis.set_directive_begin(mio::SimulationTime(5.0)); + model.parameters.get>() = npis; + mio::osecir::Simulation> sim(model, 0.0); + sim.advance(3.0); + mio::ContactMatrixGroup const& contact_matrix = + sim.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); + + // directive begin is satisfied + // t0=0: delay not enforced, t_start=0, smooth window [-1,0] -> full effect immediately at t=0 + // duration=5 -> t_end=5, t_end_damping=6: smooth window [5,6] -> back to 1.0 at t=6 + npis.set_implementation_delay(mio::SimulationTime(2.0)); // not used as t0=0 + npis.set_directive_begin(mio::SimulationTime(0.0)); + model.parameters.get>() = npis; + mio::osecir::Simulation> sim_2(model, 0.0); + sim_2.advance(3.0); + mio::ContactMatrixGroup const& contact_matrix_sim_2 = + sim_2.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(0.0))(0, 0), + 0.5); // full effect at t=0 (window in past) + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(6.0))(0, 0), + 1.0); // lifted after duration+1 + + // directive begin is satisfied but directive end ends the NPI earlier + // t_start=0, t_end=min(3,5)=3, t_end_damping=4: window [3,4] -> 1.0 at t=4 + npis.set_directive_end(mio::SimulationTime(3.0)); + model.parameters.get>() = npis; + mio::osecir::Simulation> sim_3(model, 0.0); + sim_3.advance(5.0); + mio::ContactMatrixGroup const& contact_matrix_sim_3 = + sim_3.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(2.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 1.0); // lifted at t_end+1=4 + + // directive begin is satisfied (now with delay>0 as t0=1) + // t_start=3, t_start_damping=4, smooth window [3,4] -> 0.5 at t=4 + // duration=5 -> t_end=8, t_end_damping=9, smooth window [8,9] -> 1.0 at t=9 + npis.set_implementation_delay(mio::SimulationTime(2.0)); + npis.set_directive_begin(mio::SimulationTime(0.0)); + npis.set_directive_end(mio::SimulationTime(1000000.)); + model.parameters.get>() = npis; + mio::osecir::Simulation> sim_4(model, 1.0); + sim_4.advance(5.0); + mio::ContactMatrixGroup const& contact_matrix_sim_4 = + sim_4.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); // ramp not yet done + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); // full effect + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(8.0))(0, 0), 0.5); // still active + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(9.0))(0, 0), 1.0); // lifted at t_end+1=9 + + // directive begin is satisfied but directive end ends the NPI earlier (now with delay>0 as t0=1) + // t_start=3, t_end=min(4,8)=4, t_end_damping=5: 1.0 at t=5 + npis.set_directive_end(mio::SimulationTime(4.0)); + model.parameters.get>() = npis; + mio::osecir::Simulation> sim_5(model, 1.0); + sim_5.advance(6.0); + mio::ContactMatrixGroup const& contact_matrix_sim_5 = + sim_5.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(4.0))(0, 0), + 0.5); // full effect at t_start+1=4 + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(5.0))(0, 0), 1.0); // lifted at t_end+1=5 +} + +TEST(DynamicNPIs, secir_backward_consistency_with_predefined) +{ + // Verify the core reproducibility guarantee of the +1 offset: + // A simulation with dynamic NPIs (threshold exceeded at t=1, delay=2 -> t_start=3) + // must produce the same contact matrix as a simulation with predefined dampings + // placed at t_start_damping = t_start+1 = 4 and t_end_damping = t_end+1 = 9. + mio::osecir::Model model(1); + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = 10; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, 100); + + mio::ContactMatrixGroup& cm = model.parameters.get>(); + cm[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1.0)); + + mio::DynamicNPIs npis; + npis.set_threshold(0.05 * 50'000, {mio::DampingSampling{0.5, + mio::DampingLevel(0), + mio::DampingType(0), + mio::SimulationTime(0), + {0}, + Eigen::VectorXd::Ones(1)}}); + npis.set_duration(mio::SimulationTime(5.0)); + npis.set_base_value(50'000); + npis.set_implementation_delay(mio::SimulationTime(2.0)); + model.parameters.get>() = npis; + + // t0=1: threshold exceeded at t=1, delay=2 -> t_start=3, t_start_damping=4, + // duration=5 -> t_end=8, t_end_damping=9. + // Advance only to t=8 so that the keep-alive renewal (which acts at t=9 when t > t_end=8) + // never triggers. The damping list then stays identical to the predefined dampings for all + // time points, including t=9 (restore window [8,9] fully captured in the list). + mio::osecir::Simulation> sim(model, 1.0); + sim.advance(8.0); + + // Equivalent predefined dampings: place damping at t_start_damping=4 and restore at t_end_damping=9. + // smoother_cosine uses window [t-1, t], so damping at t=4 gives window [3,4]. + mio::ContactMatrixGroup cm_expected(1, 1); + cm_expected[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1.0)); + cm_expected[0].add_damping(0.5, mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime(4.0)); + cm_expected[0].add_damping(0.0, mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime(9.0)); + + auto const& cm_dynamic = sim.get_model().parameters.get>().get_cont_freq_mat(); + for (double t : {1.0, 2.0, 3.0, 3.5, 4.0, 5.0, 7.0, 8.0, 8.5, 9.0}) { + EXPECT_DOUBLE_EQ(cm_dynamic.get_matrix_at(mio::SimulationTime(t))(0, 0), + cm_expected.get_matrix_at(mio::SimulationTime(t))(0, 0)) + << "Mismatch at t=" << t; + } +} + +TEST(DynamicNPIs, secir_expiry_keep_alive) +{ + // Verify keep-alive: when the NPI expires but the threshold is still exceeded, + // the NPI is renewed immediately with no delay gap and no dip between expiry and renewal. + // + // Timeline (delay=2, duration=5, t0=1): + // 1st NPI: t_start=3, t_start_damping=4, t_end=8, t_end_damping=9 + // At t=9 (expiry): keep-alive acts, t_start=9, t_start_damping=9 (no +1!), + // t_end=14, t_end_damping=15 -> list collapses to [(t=4,0.5),(t=15,0.0)] + mio::osecir::Model model(1); + model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = 10; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, 100); + + mio::ContactMatrixGroup& cm = model.parameters.get>(); + cm[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1.0)); + + mio::DynamicNPIs npis; + npis.set_threshold(0.05 * 50'000, {mio::DampingSampling{0.5, + mio::DampingLevel(0), + mio::DampingType(0), + mio::SimulationTime(0), + {0}, + Eigen::VectorXd::Ones(1)}}); + npis.set_duration(mio::SimulationTime(5.0)); + npis.set_base_value(50'000); + npis.set_implementation_delay(mio::SimulationTime(2.0)); + model.parameters.get>() = npis; + + // Stop at t=14 so the 2nd keep-alive (which would act at t=15) does not trigger. + mio::osecir::Simulation> sim(model, 1.0); + sim.advance(14.0); + + auto const& cm_dyn = sim.get_model().parameters.get>().get_cont_freq_mat(); + + // Damping list after keep-alive collapses to [(t=4, 0.5), (t=15, 0.0)]: + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); // 1st NPI active + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(8.0))(0, 0), 0.5); // still active + // no dip at t=9 as keep-alive overwrote the restore entry, contact stays at 0.5. + // (Without this fix the restore would restore to 1.0 at t=9.) + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(9.0))(0, 0), 0.5); + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(10.0))(0, 0), 0.5); // still constant + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(14.0))(0, 0), 0.5); // still active + // Restore window [14, 15]: complete at t=15. + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(15.0))(0, 0), 1.0); } TEST(DynamicNPIs, secirvvs_threshold_safe) @@ -573,8 +784,8 @@ TEST(DynamicNPIs, secirvvs_threshold_exceeded) Eigen::VectorXd::Ones(1)}}); npis.set_duration(mio::SimulationTime(5.0)); npis.set_base_value(50'000); - model.parameters.get>() = npis; - model.parameters.get>() = 0.0; + npis.set_implementation_delay(mio::SimulationTime(0.0)); + model.parameters.get>() = npis; EXPECT_EQ( model.parameters.get>().get_cont_freq_mat()[0].get_dampings().size(), @@ -621,30 +832,280 @@ TEST(DynamicNPIs, secirvvs_delayed_implementation) model.parameters.get>().get_cont_freq_mat()[0].get_dampings().size(), 0); - // start with t0 = 0.0 + // start with t0 = 0.0 so dynamicNPIs are active from the start + // t0=0: delay not enforced, t_start=0, smooth window [-1,0] -> full effect at t=0 mio::osecirvvs::Simulation> sim(model, 0.0); sim.advance(3.0); mio::ContactMatrixGroup const& contact_matrix = sim.get_model().parameters.template get>(); - EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(2.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 0.5); EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 0.5); - // second simulation with t0 = 1.0, so the NPIs are implemented at tmax + delay = 6.0 - const auto tmax = 4.0; - model.parameters.get>() = 2.0; + // second simulation with t0 = 1.0: t_start=3, t_start_damping=4, smooth window [3,4] + // -> at t=3 still 1.0, at t=4 full effect 0.5 + npis.set_implementation_delay(mio::SimulationTime(2.0)); + model.parameters.get>() = npis; mio::osecirvvs::Simulation> sim_2(model, 1.0); - sim_2.advance(tmax); + sim_2.advance(5.0); mio::ContactMatrixGroup const& contact_matrix_sim_2 = sim_2.get_model().parameters.template get>(); - EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(5.0))(0, 0), 1.0); - EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(6.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); // ramp not yet done + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); // full effect - // third simulation with t0 = 1.0, so the NPIs are implemented at tmax + delay = 14.0 - model.parameters.get>() = 10.0; + // third simulation: delay=10 -> t_start=11, t_start_damping=12 + npis.set_implementation_delay(mio::SimulationTime(10.0)); + model.parameters.get>() = npis; mio::osecirvvs::Simulation> sim_3(model, 1.0); sim_3.advance(4.0); mio::ContactMatrixGroup const& contact_matrix_sim_3 = sim_3.get_model().parameters.template get>(); - EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(13.0))(0, 0), 1.0); - EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(14.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(11.0))(0, 0), 1.0); // ramp start + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(12.0))(0, 0), 0.5); // full effect +} + +TEST(DynamicNPIs, secirvvs_implementation_with_directives) +{ + mio::osecirvvs::Model model(1); + model.populations[{mio::AgeGroup(0), mio::osecirvvs::InfectionState::InfectedSymptomsNaive}] = 10; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecirvvs::InfectionState::SusceptibleNaive}, + 100); + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + + mio::ContactMatrixGroup& cm = model.parameters.get>(); + cm[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1.0)); + + mio::DynamicNPIs npis; + npis.set_threshold(0.05 * 50'000, {mio::DampingSampling{0.5, + mio::DampingLevel(0), + mio::DampingType(0), + mio::SimulationTime(0), + {0}, + Eigen::VectorXd::Ones(1)}}); + npis.set_duration(mio::SimulationTime(5.0)); + npis.set_base_value(50'000); + + // directive begin is after the simulation, so no NPI is implemented + npis.set_directive_begin(mio::SimulationTime(5.0)); + model.parameters.get>() = npis; + mio::osecirvvs::Simulation> sim(model, 0.0); + sim.advance(3.0); + mio::ContactMatrixGroup const& contact_matrix = + sim.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); + + // directive begin is satisfied + // t0=0: delay not enforced, t_start=0, smooth window [-1,0] -> full effect at t=0 + // duration=5 -> t_end=5, t_end_damping=6: back to 1.0 at t=6 + npis.set_implementation_delay(mio::SimulationTime(2.0)); // not used as t0=0 + npis.set_directive_begin(mio::SimulationTime(0.0)); + model.parameters.get>() = npis; + mio::osecirvvs::Simulation> sim_2(model, 0.0); + sim_2.advance(3.0); + mio::ContactMatrixGroup const& contact_matrix_sim_2 = + sim_2.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(6.0))(0, 0), 1.0); // lifted at t_end+1=6 + + // directive begin satisfied but directive end cuts NPI earlier + // t_start=0, t_end=min(3,5)=3, t_end_damping=4: 1.0 at t=4 + npis.set_directive_end(mio::SimulationTime(3.0)); + model.parameters.get>() = npis; + mio::osecirvvs::Simulation> sim_3(model, 0.0); + sim_3.advance(5.0); + mio::ContactMatrixGroup const& contact_matrix_sim_3 = + sim_3.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(2.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 1.0); // lifted at t_end+1=4 + + // directive begin satisfied (with delay > 0 as t0 = 1) + // t_start=3, t_start_damping=4: 0.5 at t=4; t_end=8, t_end_damping=9: 1.0 at t=9 + npis.set_implementation_delay(mio::SimulationTime(2.0)); + npis.set_directive_begin(mio::SimulationTime(0.0)); + npis.set_directive_end(mio::SimulationTime(1000000.)); + model.parameters.get>() = npis; + mio::osecirvvs::Simulation> sim_4(model, 1.0); + sim_4.advance(5.0); + mio::ContactMatrixGroup const& contact_matrix_sim_4 = + sim_4.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); // ramp start + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); // full effect + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(8.0))(0, 0), 0.5); // still active + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(9.0))(0, 0), 1.0); // lifted at t_end+1=9 + + // directive begin satisfied but directive end cuts NPI earlier (with delay > 0 as t0 = 1) + // t_end=min(4,8)=4, t_end_damping=5: 1.0 at t=5 + npis.set_directive_end(mio::SimulationTime(4.0)); + model.parameters.get>() = npis; + mio::osecirvvs::Simulation> sim_5(model, 1.0); + sim_5.advance(6.0); + mio::ContactMatrixGroup const& contact_matrix_sim_5 = + sim_5.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(4.0))(0, 0), + 0.5); // full effect at t_start+1=4 + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(5.0))(0, 0), 1.0); // lifted at t_end+1=5 +} + +TEST(DynamicNPIs, osecirts_delayed_implementation) +{ + mio::osecirts::Model model(1); + model.populations[{mio::AgeGroup(0), mio::osecirts::InfectionState::InfectedSymptomsNaive}] = 10; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecirts::InfectionState::SusceptibleNaive}, + 100); + + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + + mio::ContactMatrixGroup& cm = model.parameters.get>(); + cm[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1.0)); + + mio::DynamicNPIs npis; + npis.set_threshold(0.05 * 50'000, {mio::DampingSampling{0.5, + mio::DampingLevel(0), + mio::DampingType(0), + mio::SimulationTime(0), + {0}, + Eigen::VectorXd::Ones(1)}}); + npis.set_duration(mio::SimulationTime(5.0)); + npis.set_base_value(50'000); + model.parameters.get>() = npis; + + EXPECT_EQ( + model.parameters.get>().get_cont_freq_mat()[0].get_dampings().size(), 0); + + // start with t0 = 0.0: delay not enforced, t_start=0, smooth window [-1,0] -> full effect at t=0 + mio::osecirts::Simulation> sim(model, 0.0); + sim.advance(3.0); + mio::ContactMatrixGroup const& contact_matrix = + sim.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 0.5); + + // second simulation with t0 = 1.0: t_start=3, t_start_damping=4, smooth window [3,4] + // -> at t=3 still 1.0, at t=4 full effect 0.5 + npis.set_implementation_delay(mio::SimulationTime(2.0)); + model.parameters.get>() = npis; + mio::osecirts::Simulation> sim_2(model, 1.0); + sim_2.advance(5.0); + mio::ContactMatrixGroup const& contact_matrix_sim_2 = + sim_2.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); // ramp start + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); // full effect + + // third simulation: delay=10 -> t_start=11, t_start_damping=12 + npis.set_implementation_delay(mio::SimulationTime(10.0)); + model.parameters.get>() = npis; + mio::osecirts::Simulation> sim_3_secirts(model, 1.0); + sim_3_secirts.advance(4.0); + mio::ContactMatrixGroup const& contact_matrix_sim_3_secirts = + sim_3_secirts.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_3_secirts.get_matrix_at(mio::SimulationTime(11.0))(0, 0), 1.0); // ramp start + EXPECT_EQ(contact_matrix_sim_3_secirts.get_matrix_at(mio::SimulationTime(12.0))(0, 0), 0.5); // full effect +} + +TEST(DynamicNPIs, osecirts_implementation_with_directives) +{ + mio::osecirts::Model model(1); + model.populations[{mio::AgeGroup(0), mio::osecirts::InfectionState::InfectedSymptomsNaive}] = 10; + model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecirts::InfectionState::SusceptibleNaive}, + 100); + + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + model.parameters.get>().resize(mio::SimulationDay(size_t(1000))); + model.parameters.get>().array().setConstant(0); + + mio::ContactMatrixGroup& cm = model.parameters.get>(); + cm[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 1.0)); + + mio::DynamicNPIs npis; + npis.set_threshold(0.05 * 50'000, {mio::DampingSampling{0.5, + mio::DampingLevel(0), + mio::DampingType(0), + mio::SimulationTime(0), + {0}, + Eigen::VectorXd::Ones(1)}}); + npis.set_duration(mio::SimulationTime(5.0)); + npis.set_base_value(50'000); + model.parameters.get>() = npis; + + // directive begin is after the simulation, so no NPI is implemented + npis.set_directive_begin(mio::SimulationTime(5.0)); + model.parameters.get>() = npis; + mio::osecirts::Simulation> sim(model, 0.0); + sim.advance(3.0); + mio::ContactMatrixGroup const& contact_matrix = + sim.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); + + // directive begin is satisfied (t0=0, so delay is not enforced) + // t_start=0, smooth window [-1,0] -> full effect at t=0 + // duration=5 -> t_end=5, t_end_damping=6: back to 1.0 at t=6 + npis.set_implementation_delay(mio::SimulationTime(2.0)); // not used as t0=0 + npis.set_directive_begin(mio::SimulationTime(0.0)); + model.parameters.get>() = npis; + mio::osecirts::Simulation> sim_2(model, 0.0); + sim_2.advance(3.0); + mio::ContactMatrixGroup const& contact_matrix_sim_2 = + sim_2.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_2.get_matrix_at(mio::SimulationTime(6.0))(0, 0), 1.0); // lifted at t_end+1=6 + + // directive begin is satisfied but directive end ends the NPI earlier + // t_start=0, t_end=min(3,5)=3, t_end_damping=4: 1.0 at t=4 + npis.set_directive_end(mio::SimulationTime(3.0)); + model.parameters.get>() = npis; + mio::osecirts::Simulation> sim_3(model, 0.0); + sim_3.advance(5.0); + mio::ContactMatrixGroup const& contact_matrix_sim_3 = + sim_3.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(0.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(2.0))(0, 0), 0.5); + EXPECT_EQ(contact_matrix_sim_3.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 1.0); // lifted at t_end+1=4 + + // directive begin is satisfied (now with delay>0 as t0=1) + // t_start=3, t_start_damping=4: 0.5 at t=4; t_end=8, t_end_damping=9: 1.0 at t=9 + npis.set_implementation_delay(mio::SimulationTime(2.0)); + npis.set_directive_begin(mio::SimulationTime(0.0)); + npis.set_directive_end(mio::SimulationTime(1000000.)); + model.parameters.get>() = npis; + mio::osecirts::Simulation> sim_4(model, 1.0); + sim_4.advance(5.0); + mio::ContactMatrixGroup const& contact_matrix_sim_4 = + sim_4.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(3.0))(0, 0), 1.0); // ramp start + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(4.0))(0, 0), 0.5); // full effect + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(8.0))(0, 0), 0.5); // still active + EXPECT_EQ(contact_matrix_sim_4.get_matrix_at(mio::SimulationTime(9.0))(0, 0), 1.0); // lifted at t_end+1=9 + + // directive end ends the NPI earlier (now with delay>0 as t0=1) + // t_end=min(4,8)=4, t_end_damping=5: 1.0 at t=5 + npis.set_directive_end(mio::SimulationTime(4.0)); + model.parameters.get>() = npis; + mio::osecirts::Simulation> sim_5(model, 1.0); + sim_5.advance(6.0); + mio::ContactMatrixGroup const& contact_matrix_sim_5 = + sim_5.get_model().parameters.template get>(); + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(1.0))(0, 0), 1.0); + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(4.0))(0, 0), + 0.5); // full effect at t_start+1=4 + EXPECT_EQ(contact_matrix_sim_5.get_matrix_at(mio::SimulationTime(5.0))(0, 0), 1.0); // lifted at t_end+1=5 } diff --git a/cpp/tests/test_odesecir.cpp b/cpp/tests/test_odesecir.cpp index 71025ec26f..e1524e5b70 100755 --- a/cpp/tests/test_odesecir.cpp +++ b/cpp/tests/test_odesecir.cpp @@ -1197,13 +1197,6 @@ TEST(TestOdeSecir, check_constraints_parameters) model.parameters.get>()[(mio::AgeGroup)0] = 0.3; model.parameters.set>(1.1); ASSERT_EQ(model.parameters.check_constraints(), 1); - - model.parameters.set>(1.0); - model.parameters.set>(-4); - ASSERT_EQ(model.parameters.check_constraints(), 1); - - model.parameters.set>(3); - EXPECT_EQ(model.parameters.check_constraints(), 0); } TEST(TestOdeSecir, apply_constraints_parameters) @@ -1292,12 +1285,6 @@ TEST(TestOdeSecir, apply_constraints_parameters) model.parameters.set>(1.1); EXPECT_EQ(model.parameters.apply_constraints(), 1); EXPECT_EQ(model.parameters.get>()[indx_agegroup], 0); - - model.parameters.set>(-4); - EXPECT_EQ(model.parameters.apply_constraints(), 1); - EXPECT_EQ(model.parameters.get>(), 0); - - EXPECT_EQ(model.parameters.apply_constraints(), 0); } #if defined(MEMILIO_HAS_JSONCPP) diff --git a/cpp/tests/test_odesecirts.cpp b/cpp/tests/test_odesecirts.cpp index efd416e962..c798e9154a 100755 --- a/cpp/tests/test_odesecirts.cpp +++ b/cpp/tests/test_odesecirts.cpp @@ -464,9 +464,8 @@ void set_contact_parameters(mio::osecirts::Model::ParameterSet& paramete npis.set_threshold(10.0, {mio::DampingSampling(npi_value, mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime(0), {0}, npi_groups)}); npis.set_base_value(100'000); - npis.set_interval(mio::SimulationTime(3.0)); npis.set_duration(mio::SimulationTime(14.0)); - parameters.get_end_dynamic_npis() = 10.0; //required for dynamic NPIs to have effect in this model + // npis.set_directive_end(mio::SimulationTime(10.0)); // --> can probably be removed??? } void set_covid_parameters(mio::osecirts::Model::ParameterSet& params, bool set_invalid_initial_value) diff --git a/cpp/tests/test_odesecirvvs.cpp b/cpp/tests/test_odesecirvvs.cpp index 99ba435e3a..3273cfd6fa 100755 --- a/cpp/tests/test_odesecirvvs.cpp +++ b/cpp/tests/test_odesecirvvs.cpp @@ -292,10 +292,8 @@ void set_contact_parameters(mio::osecirvvs::Model::ParameterSet& paramet npis.set_threshold(10.0, {mio::DampingSampling(npi_value, mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime(0), {0}, npi_groups)}); npis.set_base_value(100'000); - npis.set_interval(mio::SimulationTime(3.0)); + npis.set_implementation_delay(mio::SimulationTime(7.0)); npis.set_duration(mio::SimulationTime(14.0)); - parameters.get_end_dynamic_npis() = 10.0; //required for dynamic NPIs to have effect in this model - parameters.template get>() = 7; } void set_covid_parameters(mio::osecirvvs::Model::ParameterSet& params, bool set_invalid_initial_value) @@ -1532,10 +1530,6 @@ TEST(TestOdeSECIRVVS, check_constraints_parameters) model.parameters.set>(1); EXPECT_EQ(model.parameters.check_constraints(), 0); - model.parameters.set>(1); - model.parameters.set>(-4); - ASSERT_EQ(model.parameters.check_constraints(), 1); - mio::set_log_level(mio::LogLevel::warn); } @@ -1677,10 +1671,6 @@ TEST(TestOdeSECIRVVS, apply_constraints_parameters) EXPECT_EQ(model.parameters.apply_constraints(), 1); EXPECT_EQ(model.parameters.get>()[indx_agegroup], 1); - model.parameters.set>(-4); - EXPECT_EQ(model.parameters.apply_constraints(), 1); - EXPECT_EQ(model.parameters.get>(), 0); - EXPECT_EQ(model.parameters.apply_constraints(), 0); mio::set_log_level(mio::LogLevel::warn); } diff --git a/cpp/tests/test_uncertain_ad_compatibility.cpp b/cpp/tests/test_uncertain_ad_compatibility.cpp index ed14eb2139..d928b2c8f5 100644 --- a/cpp/tests/test_uncertain_ad_compatibility.cpp +++ b/cpp/tests/test_uncertain_ad_compatibility.cpp @@ -138,10 +138,9 @@ TEST(TestUncertainADCompatibility, create_model) mio::osecirvvs::Model model(num_age_groups); auto& params = model.parameters; - params.template get>() = 100; - params.template get>() = 0.0143; - params.template get>() = 0.2; - params.template get>() = 7; + params.template get>() = 100; + params.template get>() = 0.0143; + params.template get>() = 0.2; using std::floor; size_t tmax_days = static_cast(floor(tmax)); diff --git a/docs/source/cpp/models/osecir.rst b/docs/source/cpp/models/osecir.rst index 50f651c698..f07577a57f 100644 --- a/docs/source/cpp/models/osecir.rst +++ b/docs/source/cpp/models/osecir.rst @@ -262,16 +262,17 @@ A complex lockdown scenario with multiple interventions starting on a specific d A more advanced structure to automatically activate interventions based on threshold criteria is given by **DynamicNPIs**. Dynamic NPIs can be configured to trigger when the number of symptomatic infected individuals exceeds a certain relative threshold in the population. In contrast to static NPIs which are active as long as no other NPI gets implemented, dynamic NPIs are checked at regular intervals and get -activated for a defined duration when the threshold is exceeded. As above, different dampings `contact_dampings` can be assigned to different contact locations -and are then triggered all at once the threshold is exceeded. -The following example shows how to set up dynamic NPIs based on the number of 200 symptomatic infected individuals per 100,000 population. -It will be active for at least 14 days and checked every 3 days. If the last check after day 14 is negative, the NPI will be deactivated. +activated for a defined duration when the threshold is exceeded. For most realistic studies, an additional delay parameter for automated implementation +can be set. This parameter imitates a delayed reaction to exceedance of the considered threshold. +As above, different dampings `contact_dampings` can be assigned to different contact locations and are then triggered all at once the threshold +is exceeded. The following example shows how to set up dynamic NPIs based on the number of 200 symptomatic infected individuals per 100,000 population. +It will be active for at least 14 days. If the last check after day 14 is negative, the NPI will be deactivated. .. code-block:: cpp // Configure dynamic NPIs with thresholds auto& dynamic_npis = params.get>(); - dynamic_npis.set_interval(mio::SimulationTime(3.0)); // Check every 3 days + dynamic_npis.set_implementation_delay(mio::SimulationTime(0.0)); // Simulate no implementation delay dynamic_npis.set_duration(mio::SimulationTime(14.0)); // Apply for 14 days dynamic_npis.set_base_value(100'000); // Per 100,000 population dynamic_npis.set_threshold(200.0, contact_dampings); // Trigger at 200 cases per 100,000 diff --git a/docs/source/cpp/models/osecirts.rst b/docs/source/cpp/models/osecirts.rst index d8fec1bfec..62ae66d5ce 100644 --- a/docs/source/cpp/models/osecirts.rst +++ b/docs/source/cpp/models/osecirts.rst @@ -316,7 +316,7 @@ The model also supports dynamic NPIs based on epidemic thresholds: // Set threshold-based triggers for NPIs auto& dynamic_npis = model.parameters.get>(); - dynamic_npis.set_interval(mio::SimulationTime(3.0)); // Check every 3 days + dynamic_npis.set_implementation_delay(mio::SimulationTime(0.0)); // Simulate no implementation delay dynamic_npis.set_duration(mio::SimulationTime(14.0)); // Apply for 14 days dynamic_npis.set_base_value(100'000); // Per 100,000 population dynamic_npis.set_threshold(200.0, dampings); // Trigger at 200 cases per 100,000 diff --git a/docs/source/cpp/models/osecirvvs.rst b/docs/source/cpp/models/osecirvvs.rst index 6733210fe8..aceb6c999a 100644 --- a/docs/source/cpp/models/osecirvvs.rst +++ b/docs/source/cpp/models/osecirvvs.rst @@ -162,9 +162,6 @@ The model includes all parameters from the basic ODE-SECIR model plus additional * - :math:`TTC_{maxSym}` - ``TestAndTraceCapacityMaxRiskSymptoms`` - Multiplier for test and trace capacity for symptomatic cases. - * - :math:`T_{dyndelay}` - - ``DynamicNPIsImplementationDelay`` - - Delay in days for implementing dynamic NPIs after threshold exceedance. * - :math:`\lambda_{N,i}` - ``ext_inf_force_dummy`` - Force of infection for susceptibles with naive immunity. @@ -332,7 +329,7 @@ The model also supports dynamic NPIs based on epidemic thresholds: // Configure dynamic NPIs auto& dynamic_npis = params.get>(); - dynamic_npis.set_interval(mio::SimulationTime(3.0)); // Check NPI every 3 days + dynamic_npis.set_implementation_delay(mio::SimulationTime(0.0)); // Simulate no implementation delay dynamic_npis.set_duration(mio::SimulationTime(14.0)); // Apply NPI for 14 days dynamic_npis.set_base_value(100'000); // Base value to trigger NPI is population of 100,000 dynamic_npis.set_threshold(200.0, dampings); // Trigger at 200 cases per 100,000 diff --git a/pycode/examples/simulation/2020_npis_sarscov2_wildtype_germany.py b/pycode/examples/simulation/2020_npis_sarscov2_wildtype_germany.py index efb7d6963f..e0cc43c9d8 100644 --- a/pycode/examples/simulation/2020_npis_sarscov2_wildtype_germany.py +++ b/pycode/examples/simulation/2020_npis_sarscov2_wildtype_germany.py @@ -502,7 +502,7 @@ def senior_awareness(t, min, max): local_npis.append(physical_distancing_work_other(0, 0.6, 0.8)) local_npis.append(senior_awareness(0, 0.0, 0.0)) - dynamic_npis.interval = 3.0 + dynamic_npis.implementation_delay = 0.0 dynamic_npis.duration = 14.0 dynamic_npis.base_value = 100000 dynamic_npis.set_threshold(200.0, local_npis) diff --git a/pycode/examples/simulation/2021_vaccination_sarscov2_delta_germany.py b/pycode/examples/simulation/2021_vaccination_sarscov2_delta_germany.py index 533d2676bd..a055e9ecff 100644 --- a/pycode/examples/simulation/2021_vaccination_sarscov2_delta_germany.py +++ b/pycode/examples/simulation/2021_vaccination_sarscov2_delta_germany.py @@ -506,10 +506,8 @@ def senior_awareness(t, min, max): masks_high = 0.0 masks_narrow = 0.0 - start_open = mio.Date( - 2021, month_open, 1) + start_open = mio.Date(2021, month_open, 1) start_summer = start_open - self.start_date - params.end_dynamic_npis = start_summer if start_open < end_date: dampings.append(contacts_at_home(start_summer, 0.0, 0.0)) @@ -586,7 +584,8 @@ def senior_awareness(t, min, max): physical_distancing_other(0, 0.2 + narrow, 0.4 - narrow)) dynamic_npi_dampings2.append(senior_awareness(0, 0.0, 0.0)) - dynamic_npis.interval = 1.0 + dynamic_npis.implementation_delay = 7.0 + dynamic_npis.directive_end = start_summer dynamic_npis.duration = 14.0 dynamic_npis.base_value = 100000 dynamic_npis.set_threshold(35.0, dynamic_npi_dampings) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/dynamic_npis.h b/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/dynamic_npis.h index 9970dd3b1b..caae6d6110 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/dynamic_npis.h +++ b/pycode/memilio-simulation/memilio/simulation/bindings/epidemiology/dynamic_npis.h @@ -39,12 +39,28 @@ void bind_dynamicNPI_members(pybind11::module_& m, std::string const& name) bind_class(m, name.c_str()) .def(pybind11::init<>()) .def_property( - "interval", + "implementation_delay", [](mio::DynamicNPIs& self) { - return static_cast(self.get_interval()); + return static_cast(self.get_implementation_delay()); }, [](mio::DynamicNPIs& self, double v) { - self.set_interval(mio::SimulationTime(v)); + self.set_implementation_delay(mio::SimulationTime(v)); + }) + .def_property( + "directive_begin", + [](mio::DynamicNPIs& self) { + return static_cast(self.get_directive_begin()); + }, + [](mio::DynamicNPIs& self, double v) { + self.set_directive_begin(mio::SimulationTime(v)); + }) + .def_property( + "directive_end", + [](mio::DynamicNPIs& self) { + return static_cast(self.get_directive_end()); + }, + [](mio::DynamicNPIs& self, double v) { + self.set_directive_end(mio::SimulationTime(v)); }) .def_property( "duration", diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index 288c95fd95..a07991f1ba 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp @@ -145,14 +145,6 @@ PYBIND11_MODULE(_simulation_osecirvvs, m) [](mio::osecirvvs::Parameters& self, double v) { self.get_end_commuter_detection() = v; }) - .def_property( - "end_dynamic_npis", - [](const mio::osecirvvs::Parameters& self) { - return self.get_end_dynamic_npis(); - }, - [](mio::osecirvvs::Parameters& self, double v) { - self.get_end_dynamic_npis() = v; - }) .def("check_constraints", &mio::osecirvvs::Parameters::check_constraints) .def("apply_constraints", &mio::osecirvvs::Parameters::apply_constraints); @@ -239,9 +231,9 @@ PYBIND11_MODULE(_simulation_osecirvvs, m) mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity}; auto weights = std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}; auto result = mio::set_edges, - mio::MobilityParameters, mio::MobilityCoefficientGroup, - mio::osecirvvs::InfectionState, decltype(mio::read_mobility_plain)>( + ContactLocation, mio::osecirvvs::Model, + mio::MobilityParameters, mio::MobilityCoefficientGroup, + mio::osecirvvs::InfectionState, decltype(mio::read_mobility_plain)>( mobility_data_file, params_graph, mobile_comp, contact_locations_size, mio::read_mobility_plain, weights); return pymio::check_and_throw(result); diff --git a/pycode/memilio-simulation/tests/test_dynamic_npis.py b/pycode/memilio-simulation/tests/test_dynamic_npis.py index de8fc83ec4..61dd1971d0 100644 --- a/pycode/memilio-simulation/tests/test_dynamic_npis.py +++ b/pycode/memilio-simulation/tests/test_dynamic_npis.py @@ -32,10 +32,14 @@ def test_dynamic_npis(self): """ """ model = Model(0) dynamic_npis = model.parameters.DynamicNPIsInfectedSymptoms - dynamic_npis.interval = 3.0 + dynamic_npis.implementation_delay = 3.0 + dynamic_npis.directive_begin = 1.0 + dynamic_npis.directive_end = 10.0 dynamic_npis.duration = 14.0 dynamic_npis.base_value = 100000 - self.assertEqual(dynamic_npis.interval, 3.0) + self.assertEqual(dynamic_npis.implementation_delay, 3.0) + self.assertEqual(dynamic_npis.directive_begin, 1.0) + self.assertEqual(dynamic_npis.directive_end, 10.0) self.assertEqual(dynamic_npis.duration, 14.0) self.assertEqual(dynamic_npis.base_value, 100000) self.assertEqual(len(dynamic_npis.threshold), 0)