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 1 commit into
base: main
Choose a base branch
from
+269 −42
Open

1469 Correct nonageresolved indices in IDE-SECIR model #1472

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
10 changes: 6 additions & 4 deletions cpp/models/ide_secir/model.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -303,8 +303,7 @@ void Model::initial_compute_compartments(ScalarType dt)
// 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] -
Copy link
Member

@reneSchm reneSchm Jan 13, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Off topic, but if a value is used as much as this Eigen::Index(0), it helps readability to give it a name. E.g. const Eigen::Index init_time_point = 0, if I read the context correctly.

populations[Eigen::Index(0)][ISevi] -
populations[Eigen::Index(0)][ICri] -
populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Di];

populations[Eigen::Index(0)][Si] = m_N[group] - populations[Eigen::Index(0)][Ei] -
populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
Expand Down Expand Up @@ -448,13 +447,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
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