Skip to content

Commit 67e70b9

Browse files
HenrZuHenrZu
committed
add is_expiry_renewal to avoid dips if dynamic npis is extended
1 parent 470a7e7 commit 67e70b9

4 files changed

Lines changed: 112 additions & 27 deletions

File tree

‎cpp/models/ode_secir/model.h‎

Lines changed: 20 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -305,7 +305,7 @@ class Simulation : public BaseT
305305
if (t > 0) {
306306
delay_npi_implementation = FP(dyn_npis.get_implementation_delay());
307307
}
308-
else { // DynamicNPIs for t=0 are 'misused' to be 'from-start NPIs'. I.e., do not enforce delay.
308+
else { // DynamicNPIs for t=0 are treated as 'from-start NPIs'. I.e., do not enforce delay.
309309
delay_npi_implementation = 0;
310310
}
311311

@@ -320,16 +320,28 @@ class Simulation : public BaseT
320320
(exceeded_threshold->first > m_dynamic_npi.first ||
321321
t > FP(m_dynamic_npi.second))) { // old npi was weaker or is expired
322322

323-
if (t + delay_npi_implementation < direc_end) {
324-
auto t_start = SimulationTime<FP>(t + delay_npi_implementation);
323+
// Keep-alive: if the NPI expired but the threshold is still exceeded at the
324+
// same level, renew immediately without delay to avoid a gap.
325+
// Apply implementation delay only if stronger NPI needed.
326+
bool is_expiry_renewal =
327+
(t > FP(m_dynamic_npi.second)) && !(exceeded_threshold->first > m_dynamic_npi.first);
328+
FP effective_delay = is_expiry_renewal ? FP(0) : delay_npi_implementation;
329+
if (t + effective_delay < direc_end) {
330+
auto t_start = SimulationTime<FP>(t + effective_delay);
325331
// set the end to the minimum of start+duration and the end of the directive
326332
auto t_end = SimulationTime<FP>(min<FP>(direc_end, FP(t_start + dyn_npis.get_duration())));
327333
m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
328-
// For t_start > 0: shift dampings by +1 so the smooth transition window
329-
// [t_start, t_start+1] lies in the future, consistent with predefined dampings.
330-
// For t_start = 0: window [-1, 0] is in the past, so keep as is.
331-
auto t_start_damping =
332-
(FP(t_start) > FP(0)) ? SimulationTime<FP>(FP(t_start) + FP(1)) : t_start;
334+
// For new NPIs (t_start > 0): shift t_start_damping by +1 so the smooth
335+
// transition window [t_start, t_start+1] lies in the future, consistent with
336+
// predefined dampings. For t_start = 0 (global t0): the window [-1, 0]
337+
// is in the past and no shift is needed.
338+
// For keep-alive renewals (is_expiry_renewal): do not shift t_start_damping.
339+
// Since t_start == t_end_damping_old, the new start damping has the same
340+
// (time, level, type) as the previous entry.
341+
// Therefore, the contact matrix stays constant at the NPI level with no dip.
342+
auto t_start_damping = (!is_expiry_renewal && FP(t_start) > FP(0))
343+
? SimulationTime<FP>(FP(t_start) + FP(1))
344+
: t_start;
333345
auto t_end_damping = (FP(t_start) > FP(0)) ? SimulationTime<FP>(FP(t_end) + FP(1)) : t_end;
334346
implement_dynamic_npis<FP>(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second,
335347
t_start_damping, t_end_damping, [](auto& g) {

‎cpp/models/ode_secirts/model.h‎

Lines changed: 20 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -775,7 +775,7 @@ class Simulation : public BaseT
775775
delay_npi_implementation = FP(dyn_npis.get_implementation_delay());
776776
}
777777
else {
778-
// DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay.
778+
// DynamicNPIs for t=0 are treated as 'from-start NPIs'. I.e., do not enforce delay.
779779
delay_npi_implementation = 0;
780780
}
781781

@@ -790,18 +790,30 @@ class Simulation : public BaseT
790790
(exceeded_threshold->first > m_dynamic_npi.first ||
791791
t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired
792792

793-
if (t + delay_npi_implementation < direc_end) {
794-
auto t_start = SimulationTime<FP>(t + delay_npi_implementation);
793+
// Keep-alive: if the NPI expired but the threshold is still exceeded at the
794+
// same level, renew immediately without delay to avoid a gap.
795+
// Apply implementation delay only if stronger NPI needed.
796+
bool is_expiry_renewal =
797+
(t > FP(m_dynamic_npi.second)) && !(exceeded_threshold->first > m_dynamic_npi.first);
798+
FP effective_delay = is_expiry_renewal ? FP(0) : delay_npi_implementation;
799+
if (t + effective_delay < direc_end) {
800+
auto t_start = SimulationTime<FP>(t + effective_delay);
795801
// set the end to the minimum of start+duration and the end of the directive
796802
auto t_end = SimulationTime<FP>(min<FP>(direc_end, FP(t_start + dyn_npis.get_duration())));
797803
this->get_model().parameters.get_start_commuter_detection() = (FP)t_start;
798804
this->get_model().parameters.get_end_commuter_detection() = (FP)t_end;
799805
m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
800-
// For t_start > 0: shift dampings by +1 so the smooth transition window
801-
// [t_start, t_start+1] lies in the future, consistent with predefined dampings.
802-
// For t_start = 0: window [-1, 0] is in the past, so keep as is.
803-
auto t_start_damping =
804-
(FP(t_start) > FP(0)) ? SimulationTime<FP>(FP(t_start) + FP(1)) : t_start;
806+
// For new NPIs (t_start > 0): shift t_start_damping by +1 so the smooth
807+
// transition window [t_start, t_start+1] lies in the future, consistent with
808+
// predefined dampings. For t_start = 0 (global t0): the window [-1, 0]
809+
// is in the past and no shift is needed.
810+
// For keep-alive renewals (is_expiry_renewal): do not shift t_start_damping.
811+
// Since t_start == t_end_damping_old, the new start damping has the same
812+
// (time, level, type) as the previous entry.
813+
// Therefore, the contact matrix stays constant at the NPI level with no dip.
814+
auto t_start_damping = (!is_expiry_renewal && FP(t_start) > FP(0))
815+
? SimulationTime<FP>(FP(t_start) + FP(1))
816+
: t_start;
805817
auto t_end_damping = (FP(t_start) > FP(0)) ? SimulationTime<FP>(FP(t_end) + FP(1)) : t_end;
806818
implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second,
807819
t_start_damping, t_end_damping, [](auto& g) {

‎cpp/models/ode_secirvvs/model.h‎

Lines changed: 20 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -696,7 +696,7 @@ class Simulation : public BaseT
696696
if (t > 0) {
697697
delay_npi_implementation = FP(dyn_npis.get_implementation_delay());
698698
}
699-
else { // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay.
699+
else { // DynamicNPIs for t=0 are treated as 'from-start NPIs'. I.e., do not enforce delay.
700700
delay_npi_implementation = 0;
701701
}
702702
if (t == 0) {
@@ -715,18 +715,30 @@ class Simulation : public BaseT
715715
(exceeded_threshold->first > m_dynamic_npi.first ||
716716
t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired
717717

718-
if (t + delay_npi_implementation < direc_end) {
719-
auto t_start = SimulationTime<FP>(t + delay_npi_implementation);
718+
// Keep-alive: if the NPI expired but the threshold is still exceeded at the
719+
// same level, renew immediately without delay to avoid a gap.
720+
// Apply implementation delay only if stronger NPI needed.
721+
bool is_expiry_renewal =
722+
(t > FP(m_dynamic_npi.second)) && !(exceeded_threshold->first > m_dynamic_npi.first);
723+
FP effective_delay = is_expiry_renewal ? FP(0) : delay_npi_implementation;
724+
if (t + effective_delay < direc_end) {
725+
auto t_start = SimulationTime<FP>(t + effective_delay);
720726
// set the end to the minimum of start+duration and the end of the directive
721727
auto t_end = SimulationTime<FP>(min<FP>(direc_end, FP(t_start + dyn_npis.get_duration())));
722728
this->get_model().parameters.get_start_commuter_detection() = (FP)t_start;
723729
this->get_model().parameters.get_end_commuter_detection() = (FP)t_end;
724730
m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
725-
// For t_start > 0: shift dampings by +1 so the smooth transition window
726-
// [t_start, t_start+1] lies in the future, consistent with predefined dampings.
727-
// For t_start = 0: window [-1, 0] is in the past, so keep as is.
728-
auto t_start_damping =
729-
(FP(t_start) > FP(0)) ? SimulationTime<FP>(FP(t_start) + FP(1)) : t_start;
731+
// For new NPIs (t_start > 0): shift t_start_damping by +1 so the smooth
732+
// transition window [t_start, t_start+1] lies in the future, consistent with
733+
// predefined dampings. For t_start = 0 (global t0): the window [-1, 0]
734+
// is in the past and no shift is needed.
735+
// For keep-alive renewals (is_expiry_renewal): do not shift t_start_damping.
736+
// Since t_start == t_end_damping_old, the new start damping has the same
737+
// (time, level, type) as the previous entry.
738+
// Therefore, the contact matrix stays constant at the NPI level with no dip.
739+
auto t_start_damping = (!is_expiry_renewal && FP(t_start) > FP(0))
740+
? SimulationTime<FP>(FP(t_start) + FP(1))
741+
: t_start;
730742
auto t_end_damping = (FP(t_start) > FP(0)) ? SimulationTime<FP>(FP(t_end) + FP(1)) : t_end;
731743
implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second,
732744
t_start_damping, t_end_damping, [](auto& g) {

‎cpp/tests/test_dynamic_npis.cpp‎

Lines changed: 52 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -657,9 +657,12 @@ TEST(DynamicNPIs, secir_backward_consistency_with_predefined)
657657
model.parameters.get<mio::osecir::DynamicNPIsInfectedSymptoms<double>>() = npis;
658658

659659
// t0=1: threshold exceeded at t=1, delay=2 -> t_start=3, t_start_damping=4,
660-
// duration=5 -> t_end=8, t_end_damping=9
660+
// duration=5 -> t_end=8, t_end_damping=9.
661+
// Advance only to t=8 so that the keep-alive renewal (which acts at t=9 when t > t_end=8)
662+
// never triggers. The damping list then stays identical to the predefined dampings for all
663+
// time points, including t=9 (restore window [8,9] fully captured in the list).
661664
mio::osecir::Simulation<double, mio_test::MockSimulation<mio::osecir::Model>> sim(model, 1.0);
662-
sim.advance(10.0);
665+
sim.advance(8.0);
663666

664667
// Equivalent predefined dampings: place damping at t_start_damping=4 and restore at t_end_damping=9.
665668
// smoother_cosine uses window [t-1, t], so damping at t=4 gives window [3,4].
@@ -669,13 +672,59 @@ TEST(DynamicNPIs, secir_backward_consistency_with_predefined)
669672
cm_expected[0].add_damping(0.0, mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime<double>(9.0));
670673

671674
auto const& cm_dynamic = sim.get_model().parameters.get<mio::osecir::ContactPatterns<double>>().get_cont_freq_mat();
672-
for (double t : {1.0, 2.0, 3.0, 3.5, 4.0, 5.0, 7.0, 8.0, 8.5, 9.0, 10.0}) {
675+
for (double t : {1.0, 2.0, 3.0, 3.5, 4.0, 5.0, 7.0, 8.0, 8.5, 9.0}) {
673676
EXPECT_DOUBLE_EQ(cm_dynamic.get_matrix_at(mio::SimulationTime<double>(t))(0, 0),
674677
cm_expected.get_matrix_at(mio::SimulationTime<double>(t))(0, 0))
675678
<< "Mismatch at t=" << t;
676679
}
677680
}
678681

682+
TEST(DynamicNPIs, secir_expiry_keep_alive)
683+
{
684+
// Verify keep-alive: when the NPI expires but the threshold is still exceeded,
685+
// the NPI is renewed immediately with no delay gap and no dip between expiry and renewal.
686+
//
687+
// Timeline (delay=2, duration=5, t0=1):
688+
// 1st NPI: t_start=3, t_start_damping=4, t_end=8, t_end_damping=9
689+
// At t=9 (expiry): keep-alive acts, t_start=9, t_start_damping=9 (no +1!),
690+
// t_end=14, t_end_damping=15 -> list collapses to [(t=4,0.5),(t=15,0.0)]
691+
mio::osecir::Model<double> model(1);
692+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = 10;
693+
model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, 100);
694+
695+
mio::ContactMatrixGroup<double>& cm = model.parameters.get<mio::osecir::ContactPatterns<double>>();
696+
cm[0] = mio::ContactMatrix<double>(Eigen::MatrixXd::Constant(1, 1, 1.0));
697+
698+
mio::DynamicNPIs<double> npis;
699+
npis.set_threshold(0.05 * 50'000, {mio::DampingSampling<double>{0.5,
700+
mio::DampingLevel(0),
701+
mio::DampingType(0),
702+
mio::SimulationTime<double>(0),
703+
{0},
704+
Eigen::VectorXd::Ones(1)}});
705+
npis.set_duration(mio::SimulationTime<double>(5.0));
706+
npis.set_base_value(50'000);
707+
npis.set_implementation_delay(mio::SimulationTime<double>(2.0));
708+
model.parameters.get<mio::osecir::DynamicNPIsInfectedSymptoms<double>>() = npis;
709+
710+
// Stop at t=14 so the 2nd keep-alive (which would act at t=15) does not trigger.
711+
mio::osecir::Simulation<double, mio_test::MockSimulation<mio::osecir::Model>> sim(model, 1.0);
712+
sim.advance(14.0);
713+
714+
auto const& cm_dyn = sim.get_model().parameters.get<mio::osecir::ContactPatterns<double>>().get_cont_freq_mat();
715+
716+
// Damping list after keep-alive collapses to [(t=4, 0.5), (t=15, 0.0)]:
717+
EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime<double>(4.0))(0, 0), 0.5); // 1st NPI active
718+
EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime<double>(8.0))(0, 0), 0.5); // still active
719+
// no dip at t=9 as keep-alive overwrote the restore entry, contact stays at 0.5.
720+
// (Without this fix the restore would restore to 1.0 at t=9.)
721+
EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime<double>(9.0))(0, 0), 0.5);
722+
EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime<double>(10.0))(0, 0), 0.5); // still constant
723+
EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime<double>(14.0))(0, 0), 0.5); // still active
724+
// Restore window [14, 15]: complete at t=15.
725+
EXPECT_DOUBLE_EQ(cm_dyn.get_matrix_at(mio::SimulationTime<double>(15.0))(0, 0), 1.0);
726+
}
727+
679728
TEST(DynamicNPIs, secirvvs_threshold_safe)
680729
{
681730
mio::osecirvvs::Model<double> model(1);

0 commit comments

Comments
 (0)