diff --git a/cpp/benchmarks/flow_simulation_ode_secirvvs.h b/cpp/benchmarks/flow_simulation_ode_secirvvs.h index 82a262cdf7..3df9d65c93 100644 --- a/cpp/benchmarks/flow_simulation_ode_secirvvs.h +++ b/cpp/benchmarks/flow_simulation_ode_secirvvs.h @@ -570,12 +570,17 @@ class Simulation : public Base 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); + const bool npi_expired = + floating_point_greater_equal(t, ScalarType(m_dynamic_npi.second), 1e-10); 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); + (exceeded_threshold->first > m_dynamic_npi.first || npi_expired)) { + // old npi was weaker or is expired + + const bool is_expiry_renewal = + npi_expired && !(exceeded_threshold->first > m_dynamic_npi.first); + ScalarType effective_delay = is_expiry_renewal ? ScalarType(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+delay and the end of the directive auto t_end = SimulationTime( std::min(direc_end, ScalarType(t_start + dyn_npis.get_duration()))); diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index dadef6739e..1343369af2 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -551,9 +551,10 @@ void mio::MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& 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); + const bool npi_expired = floating_point_greater_equal(t, FP(m_dynamic_npi.second), 1e-10); 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 + (exceeded_threshold->first > m_dynamic_npi.first || npi_expired)) { + //old NPI was weaker or is expired auto t_end = SimulationTime(t + dyn_npis.get_duration().get()); m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end); implement_dynamic_npis(m_parameters.get_coefficients(), exceeded_threshold->second, diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h index 154b6d0337..664d26f09c 100644 --- a/cpp/models/ode_secir/model.h +++ b/cpp/models/ode_secir/model.h @@ -315,15 +315,16 @@ class Simulation : public BaseT 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); + const bool npi_expired = floating_point_greater_equal(t, FP(m_dynamic_npi.second), 1e-10); 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 + (exceeded_threshold->first > m_dynamic_npi.first || npi_expired)) { + // 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); + npi_expired && !(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); @@ -341,7 +342,7 @@ class Simulation : public BaseT 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; + auto t_end_damping = SimulationTime(FP(t_end) + FP(1)); implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second, t_start_damping, t_end_damping, [](auto& g) { return make_contact_damping_matrix(g); diff --git a/cpp/models/ode_secirts/model.h b/cpp/models/ode_secirts/model.h index 9e05e543c9..78119b3c17 100644 --- a/cpp/models/ode_secirts/model.h +++ b/cpp/models/ode_secirts/model.h @@ -785,15 +785,16 @@ class Simulation : public BaseT 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); + const bool npi_expired = floating_point_greater_equal(t, FP(m_dynamic_npi.second), 1e-10); 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 + (exceeded_threshold->first > m_dynamic_npi.first || npi_expired)) { + //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); + npi_expired && !(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); @@ -813,7 +814,7 @@ class Simulation : public BaseT 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; + auto t_end_damping = SimulationTime(FP(t_end) + FP(1)); implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second, t_start_damping, t_end_damping, [](auto& g) { return make_contact_damping_matrix(g); diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h index 785756c0c3..9ed73d0abb 100644 --- a/cpp/models/ode_secirvvs/model.h +++ b/cpp/models/ode_secirvvs/model.h @@ -710,15 +710,16 @@ class Simulation : public BaseT 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); + const bool npi_expired = floating_point_greater_equal(t, FP(m_dynamic_npi.second), 1e-10); 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 + (exceeded_threshold->first > m_dynamic_npi.first || npi_expired)) { + //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); + npi_expired && !(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); @@ -738,7 +739,7 @@ class Simulation : public BaseT 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; + auto t_end_damping = SimulationTime(FP(t_end) + FP(1)); implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second, t_start_damping, t_end_damping, [](auto& g) { return make_contact_damping_matrix(g); diff --git a/cpp/tests/data/results_osecirvvs.h5 b/cpp/tests/data/results_osecirvvs.h5 index 34707a720b..432ace94dc 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 54baa7a091..aeeffeebd0 100644 --- a/cpp/tests/test_dynamic_npis.cpp +++ b/cpp/tests/test_dynamic_npis.cpp @@ -682,8 +682,8 @@ TEST(DynamicNPIs, secir_expiry_keep_alive) // // 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)] + // At t=8 (expiry): keep-alive acts, t_start=8, t_start_damping=8 (no +1!), + // t_end=13, t_end_damping=14 -> list collapses to [(t=4,0.5),(t=14,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); @@ -703,22 +703,22 @@ TEST(DynamicNPIs, secir_expiry_keep_alive) 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. + // Stop at t=12 so the 2nd keep-alive (which would act at t=13) does not trigger. mio::osecir::Simulation> sim(model, 1.0); - sim.advance(14.0); + sim.advance(12.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)]: + // Damping list after keep-alive collapses to [(t=4, 0.5), (t=14, 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); + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(13.0))(0, 0), 0.5); // still active + // Restore window [13, 14]: complete at t=14. + EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime(14.0))(0, 0), 1.0); } TEST(DynamicNPIs, secirvvs_threshold_safe)