Skip to content

Commit 8fd7324

Browse files
annawendlerannawendler
committed
correct indices and add test to check correctness of ageresolved simulation
1 parent 79088be commit 8fd7324

File tree

4 files changed

+269
-42
lines changed

4 files changed

+269
-42
lines changed

‎cpp/models/ide_secir/model.cpp‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -303,8 +303,7 @@ void Model::initial_compute_compartments(ScalarType dt)
303303
// The scheme of the ODE model for initialization is applied here.
304304
populations[Eigen::Index(0)][Ri] = total_confirmed_cases[group] - populations[Eigen::Index(0)][ISyi] -
305305
populations[Eigen::Index(0)][ISevi] -
306-
populations[Eigen::Index(0)][ICri] -
307-
populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
306+
populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Di];
308307

309308
populations[Eigen::Index(0)][Si] = m_N[group] - populations[Eigen::Index(0)][Ei] -
310309
populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
@@ -448,13 +447,16 @@ void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx
448447
Eigen::Index calc_time_index =
449448
(Eigen::Index)std::ceil(m_transitiondistributions_support_max[group][idx_InfectionTransitions] / dt);
450449

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

458+
// Index referring to the here computed transition in TimeSeries transitions.
459+
int transition_idx = get_transition_flat_index(idx_InfectionTransitions, size_t(group));
458460
transitions.get_value(current_time_index)[transition_idx] =
459461
(-dt) * parameters.get<TransitionProbabilities>()[group][idx_InfectionTransitions] * sum;
460462
}

‎cpp/tests/test_ide_parameters_io.cpp‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -203,16 +203,16 @@ TEST(TestIDEParametersIo, RKIcompareWithPreviousRunAgeRes)
203203

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

217217
mio::isecir::Simulation sim(model, dt);
218218
ASSERT_EQ(compare.size(), model.transitions.get_last_value().size());

‎cpp/tests/test_ide_secir.cpp‎

Lines changed: 18 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -211,22 +211,21 @@ TEST(IdeSecir, checkStartTime)
211211
}
212212

213213
// Check results of our simulation with an example calculated by hand,
214-
// for calculations see internal Overleaf document.
215-
// TODO: Add link to material when published.
214+
// see https://doi.org/10.1016/j.amc.2025.129636 for the used formulas.
216215
TEST(IdeSecir, checkSimulationFunctions)
217216
{
218217
using Vec = mio::TimeSeries<ScalarType>::Vector;
219218
size_t num_agegroups = 1;
220219
ScalarType tmax = 0.5;
220+
ScalarType dt = 0.5;
221+
221222
mio::CustomIndexArray<ScalarType, mio::AgeGroup> N =
222223
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10000.);
223224
mio::CustomIndexArray<ScalarType, mio::AgeGroup> deaths =
224225
mio::CustomIndexArray<ScalarType, mio::AgeGroup>(mio::AgeGroup(num_agegroups), 10.);
225-
ScalarType dt = 0.5;
226-
227-
int num_transitions = (int)mio::isecir::InfectionTransition::Count;
228226

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

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

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

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

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

0 commit comments

Comments
 (0)