Skip to content

1469 Correct nonageresolved indices in IDE-SECIR model #1472

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
annawendler wants to merge 2 commits into
base: main
Choose a base branch
from
Open

1469 Correct nonageresolved indices in IDE-SECIR model #1472

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8fd7324
correct indices and add test to check correctness of ageresolved simu…
annawendler Jan 9, 2026
4aa8008
review suggestions
annawendler Feb 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 35 additions & 30 deletions cpp/models/ide_secir/model.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -267,19 +267,22 @@ void Model::initial_compute_compartments(ScalarType dt)
// the transitions are used.
initial_compute_compartments_infection(dt);

// Define index for initial time point of populations which defines the start of the simulation.
const Eigen::Index start_time_point = 0;

// We store in two Booleans if there are Susceptibles or Recovered given for every age group.
bool susceptibles_given = true;
for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
if (populations[Eigen::Index(0)][Si] < 1e-12) {
if (populations[start_time_point][Si] < 1e-12) {
susceptibles_given = false;
break;
}
}
bool recovered_given = true;
for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
if (populations[Eigen::Index(0)][Ri] < 1e-12) {
if (populations[start_time_point][Ri] < 1e-12) {
recovered_given = false;
break;
}
Expand All @@ -301,16 +304,15 @@ void Model::initial_compute_compartments(ScalarType dt)
int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
// The scheme of the ODE model for initialization is applied here.
populations[Eigen::Index(0)][Ri] = total_confirmed_cases[group] - populations[Eigen::Index(0)][ISyi] -
populations[Eigen::Index(0)][ISevi] -
populations[Eigen::Index(0)][ICri] -
populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];

populations[Eigen::Index(0)][Si] = m_N[group] - populations[Eigen::Index(0)][Ei] -
populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
populations[Eigen::Index(0)][ISevi] -
populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Ri] -
populations[Eigen::Index(0)][Di];
populations[start_time_point][Ri] = total_confirmed_cases[group] - populations[start_time_point][ISyi] -
populations[start_time_point][ISevi] -
populations[start_time_point][ICri] - populations[start_time_point][Di];

populations[start_time_point][Si] =
m_N[group] - populations[start_time_point][Ei] - populations[start_time_point][INSi] -
populations[start_time_point][ISyi] - populations[start_time_point][ISevi] -
populations[start_time_point][ICri] - populations[start_time_point][Ri] -
populations[start_time_point][Di];
}
}

Expand All @@ -327,11 +329,11 @@ void Model::initial_compute_compartments(ScalarType dt)
int ICri = get_state_flat_index(Eigen::Index(InfectionState::InfectedCritical), group);
int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
populations[Eigen::Index(0)][Ri] = m_N[group] - populations[Eigen::Index(0)][Si] -
populations[Eigen::Index(0)][Ei] - populations[Eigen::Index(0)][INSi] -
populations[Eigen::Index(0)][ISyi] -
populations[Eigen::Index(0)][ISevi] -
populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Di];
populations[start_time_point][Ri] =
m_N[group] - populations[start_time_point][Si] - populations[start_time_point][Ei] -
populations[start_time_point][INSi] - populations[start_time_point][ISyi] -
populations[start_time_point][ISevi] - populations[start_time_point][ICri] -
populations[start_time_point][Di];
}
}
else if (recovered_given) {
Expand All @@ -347,11 +349,11 @@ void Model::initial_compute_compartments(ScalarType dt)
int ICri = get_state_flat_index(Eigen::Index(InfectionState::InfectedCritical), group);
int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
populations[Eigen::Index(0)][Si] = m_N[group] - populations[Eigen::Index(0)][Ei] -
populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
populations[Eigen::Index(0)][ISevi] -
populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Ri] -
populations[Eigen::Index(0)][Di];
populations[start_time_point][Si] =
m_N[group] - populations[start_time_point][Ei] - populations[start_time_point][INSi] -
populations[start_time_point][ISyi] - populations[start_time_point][ISevi] -
populations[start_time_point][ICri] - populations[start_time_point][Ri] -
populations[start_time_point][Di];
}
}
else {
Expand Down Expand Up @@ -379,16 +381,16 @@ void Model::initial_compute_compartments(ScalarType dt)
with a number less than zero, so that a log_error is thrown.
However, this initialization method is consistent with the numerical solver of the model equations,
so it may sometimes make sense to use this method. */
populations[Eigen::Index(0)][Si] =
populations[start_time_point][Si] =
transitions.get_last_value()[StEi] / (dt * m_forceofinfection[group]);

// Recovered; calculated as the difference between the total population and the sum of the other
// compartment sizes.
populations[Eigen::Index(0)][Ri] =
m_N[group] - populations[Eigen::Index(0)][Si] - populations[Eigen::Index(0)][Ei] -
populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
populations[Eigen::Index(0)][ISevi] - populations[Eigen::Index(0)][ICri] -
populations[Eigen::Index(0)][Di];
populations[start_time_point][Ri] =
m_N[group] - populations[start_time_point][Si] - populations[start_time_point][Ei] -
populations[start_time_point][INSi] - populations[start_time_point][ISyi] -
populations[start_time_point][ISevi] - populations[start_time_point][ICri] -
populations[start_time_point][Di];
}
}
else {
Expand Down Expand Up @@ -448,13 +450,16 @@ void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx
Eigen::Index calc_time_index =
(Eigen::Index)std::ceil(m_transitiondistributions_support_max[group][idx_InfectionTransitions] / dt);

int transition_idx = get_transition_flat_index(idx_InfectionTransitions, size_t(group));
// Index referring to the incoming flow in TimeSeries transitions.
int inflow_index = get_transition_flat_index(idx_IncomingFlow, size_t(group));
for (Eigen::Index i = current_time_index - calc_time_index; i < current_time_index; i++) {
// (current_time_index - i) is the index corresponding to time the individuals have already spent in this state.
sum += m_transitiondistributions_derivative[group][idx_InfectionTransitions][current_time_index - i] *
transitions[i + 1][idx_IncomingFlow];
transitions[i + 1][inflow_index];
}

// Index referring to the here computed transition in TimeSeries transitions.
int transition_idx = get_transition_flat_index(idx_InfectionTransitions, size_t(group));
transitions.get_value(current_time_index)[transition_idx] =
(-dt) * parameters.get<TransitionProbabilities>()[group][idx_InfectionTransitions] * sum;
}
Expand Down
2 changes: 1 addition & 1 deletion cpp/models/ide_secir/simulation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ void Simulation::advance(ScalarType tmax)
m_model->transitions.add_time_point(m_model->transitions.get_last_time() + m_dt);
m_model->populations.add_time_point(m_model->populations.get_last_time() + m_dt);

// compute Susceptibles:
// Compute Susceptibles:
m_model->compute_susceptibles(m_dt);

// Compute flows:
Expand Down
20 changes: 10 additions & 10 deletions cpp/tests/test_ide_parameters_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,16 +203,16 @@ TEST(TestIDEParametersIo, RKIcompareWithPreviousRunAgeRes)

// Compare transitions at last time point with results from a previous run that are given here.
Eigen::VectorX<ScalarType> compare(num_transitions * num_agegroups);
compare << 336.428571428600, 328.285714285701, 162.000000000000, 163.071428571425, 80.130989648839, 79.803571428575,
39.476374533415, 39.476374533415, 19.550404043081, 19.550404043081, 1105.714285714297, 1069.857142857200,
515.714285714250, 163.071428571425, 80.130989648839, 79.803571428575, 39.476374533415, 39.476374533415,
19.550404043081, 19.550404043081, 5819.000000000000, 5744.000000000000, 2806.428571428551, 163.071428571425,
80.130989648839, 79.803571428575, 39.476374533415, 39.476374533415, 19.550404043081, 19.550404043081,
6685.142857142899, 6572.857142857101, 3200.714285714304, 163.071428571425, 80.130989648839, 79.803571428575,
39.476374533415, 39.476374533415, 19.550404043081, 19.550404043081, 2376.000000000000, 2342.285714285696,
1142.571428571406, 163.071428571425, 80.130989648839, 79.803571428575, 39.476374533415, 39.476374533415,
19.550404043081, 19.550404043081, 966.714285714304, 946.142857142797, 457.214285714301, 163.071428571425,
80.130989648839, 79.803571428575, 39.476374533415, 39.476374533415, 19.550404043081, 19.550404043081;
compare << 336.4285714286, 328.2857142857, 162.0000000000, 163.0714285714, 80.1309896488, 79.8035714286,
39.4763745334, 39.4763745334, 19.5504040431, 19.5504040431, 1105.7142857143, 1069.8571428572, 515.7142857142,
525.3214285714, 254.4848638825, 253.2142857143, 124.6903391854, 124.6903391854, 61.1148381577, 61.1148381577,
5819.0000000000, 5744.0000000000, 2806.4285714286, 2839.2142857143, 1394.5112118989, 1391.2321428571,
689.6518098426, 689.6518098426, 340.5829657792, 340.5829657792, 6685.1428571429, 6572.8571428571,
3200.7142857143, 3243.5714285714, 1591.9783266355, 1588.8214285714, 787.5403666800, 787.5403666800,
388.5439635160, 388.5439635160, 2376.0000000000, 2342.2857142857, 1142.5714285714, 1156.8571428571,
566.0586818750, 564.0892857143, 277.9264675819, 277.9264675819, 136.1445518771, 136.1445518771, 966.7142857143,
946.1428571428, 457.2142857143, 465.1428571428, 226.4021912199, 225.5714285714, 110.8246575673, 110.8246575673,
54.0704051215, 54.0704051215;

mio::isecir::Simulation sim(model, dt);
ASSERT_EQ(compare.size(), model.transitions.get_last_value().size());
Expand Down
37 changes: 18 additions & 19 deletions cpp/tests/test_ide_secir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,22 +211,21 @@ TEST(IdeSecir, checkStartTime)
}

// Check results of our simulation with an example calculated by hand,
// for calculations see internal Overleaf document.
// TODO: Add link to material when published.
// see https://doi.org/10.1016/j.amc.2025.129636 for the used formulas.
TEST(IdeSecir, checkSimulationFunctions)
{
using Vec = mio::TimeSeries<ScalarType>::Vector;
size_t num_agegroups = 1;
ScalarType tmax = 0.5;
ScalarType dt = 0.5;

mio::CustomIndexArray<ScalarType, mio::AgeGroup> N =
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10000.);
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths =
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10.);
ScalarType dt = 0.5;

int num_transitions = (int)mio::isecir::InfectionTransition::Count;

// Create TimeSeries with num_transitions elements where transitions needed for simulation will be stored.
int num_transitions = (int)mio::isecir::InfectionTransition::Count;
mio::TimeSeries<ScalarType> init(num_transitions * num_agegroups);

// Add time points for initialization for transitions and death.
Expand Down Expand Up @@ -274,26 +273,26 @@ TEST(IdeSecir, checkSimulationFunctions)
mio::TimeSeries<ScalarType> secihurd_simulated = sim.get_result();
mio::TimeSeries<ScalarType> transitions_simulated = sim.get_transitions();

// Define vectors for compartments and transitions with values from example
// (calculated by hand, see internal Overleaf document).
// TODO: Add link to material when published.
Vec secihurd0((int)mio::isecir::InfectionState::Count);
Vec secihurd1((int)mio::isecir::InfectionState::Count);
Vec transitions1(num_transitions);
secihurd0 << 4995, 0.5, 0, 4, 0, 0, 4990.5, 10;
secihurd1 << 4994.00020016, 0.49989992, 0.49994996, 0.12498749, 1.03124687, 0.25781172, 4993.45699802, 10.12890586;
transitions1 << 0.99979984, 0.99989992, 0.24997498, 0.24997498, 2.06249374, 2.06249374, 0.51562344, 0.51562344,
// Define vectors for compartments and transitions at t0 and t1 with values from example
// (calculated by hand, see https://doi.org/10.1016/j.amc.2025.129636 for the used formulas).
Vec secihurd_t0((int)mio::isecir::InfectionState::Count);
Vec secihurd_t1((int)mio::isecir::InfectionState::Count);
Vec transitions_t1(num_transitions);
secihurd_t0 << 4995, 0.5, 0, 4, 0, 0, 4990.5, 10;
secihurd_t1 << 4994.00020016, 0.49989992, 0.49994996, 0.12498749, 1.03124687, 0.25781172, 4993.45699802,
10.12890586;
transitions_t1 << 0.99979984, 0.99989992, 0.24997498, 0.24997498, 2.06249374, 2.06249374, 0.51562344, 0.51562344,
0.12890586, 0.12890586;

// Compare SECIHURD compartments at times 0 and 1.
// Compare simulated compartments at time points t0 and t1.
for (Eigen::Index i = 0; i < (Eigen::Index)mio::isecir::InfectionState::Count; i++) {
EXPECT_NEAR(secihurd_simulated[0][i], secihurd0[i], 1e-8);
EXPECT_NEAR(secihurd_simulated[1][i], secihurd1[i], 1e-8);
EXPECT_NEAR(secihurd_simulated[0][i], secihurd_t0[i], 1e-8);
EXPECT_NEAR(secihurd_simulated[1][i], secihurd_t1[i], 1e-8);
}

// Compare transitions at time 1.
// Compare simulated transitions with expected results at time point t1.
for (Eigen::Index i = 0; i < num_transitions; i++) {
EXPECT_NEAR(transitions_simulated[transitions_simulated.get_num_time_points() - 1][i], transitions1[i], 1e-8);
EXPECT_NEAR(transitions_simulated[transitions_simulated.get_num_time_points() - 1][i], transitions_t1[i], 1e-8);
}
}

Expand Down
Loading