From 8fd732472f367f43cba2ecea5650d2d5cd62c16e Mon Sep 17 00:00:00 2001 From: Anna Wendler <106674756+annawendler@users.noreply.github.com> Date: Fri, 9 Jan 2026 16:42:26 +0100 Subject: [PATCH] correct indices and add test to check correctness of ageresolved simulation --- cpp/models/ide_secir/model.cpp | 10 +- cpp/tests/test_ide_parameters_io.cpp | 20 +-- cpp/tests/test_ide_secir.cpp | 37 ++-- cpp/tests/test_ide_secir_ageres.cpp | 244 ++++++++++++++++++++++++++- 4 files changed, 269 insertions(+), 42 deletions(-) diff --git a/cpp/models/ide_secir/model.cpp b/cpp/models/ide_secir/model.cpp index fb2f8c3aff..e17870f28a 100644 --- a/cpp/models/ide_secir/model.cpp +++ b/cpp/models/ide_secir/model.cpp @@ -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] - 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] - @@ -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()[group][idx_InfectionTransitions] * sum; } diff --git a/cpp/tests/test_ide_parameters_io.cpp b/cpp/tests/test_ide_parameters_io.cpp index 3834d38e69..e14ea00f0f 100644 --- a/cpp/tests/test_ide_parameters_io.cpp +++ b/cpp/tests/test_ide_parameters_io.cpp @@ -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 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()); diff --git a/cpp/tests/test_ide_secir.cpp b/cpp/tests/test_ide_secir.cpp index 09655e231b..9085aec22f 100755 --- a/cpp/tests/test_ide_secir.cpp +++ b/cpp/tests/test_ide_secir.cpp @@ -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::Vector; size_t num_agegroups = 1; ScalarType tmax = 0.5; + ScalarType dt = 0.5; + mio::CustomIndexArray N = mio::CustomIndexArray(mio::AgeGroup(num_agegroups), 10000.); mio::CustomIndexArray deaths = mio::CustomIndexArray(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 init(num_transitions * num_agegroups); // Add time points for initialization for transitions and death. @@ -274,26 +273,26 @@ TEST(IdeSecir, checkSimulationFunctions) mio::TimeSeries secihurd_simulated = sim.get_result(); mio::TimeSeries 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); } } diff --git a/cpp/tests/test_ide_secir_ageres.cpp b/cpp/tests/test_ide_secir_ageres.cpp index 34f9006699..c3f72b349f 100644 --- a/cpp/tests/test_ide_secir_ageres.cpp +++ b/cpp/tests/test_ide_secir_ageres.cpp @@ -149,10 +149,10 @@ TEST(TestIdeAgeres, compareWithPreviousRun) // Compare compartments at last time point with results from a previous run that are given here. mio::TimeSeries compartments = sim.get_result(); Eigen::VectorX compare_compartments(num_compartments * num_agegroups); - compare_compartments << 484.3056557672, 15.7685031055, 22.7020934123, 7.0615933479, 3.3491460693, 1.5803397070, - 4454.5548070034, 10.6778615873, 484.3056557672, 31.0934010790, 21.1271954388, 23.6370809253, 3.9106794140, - 1.7242153411, 4424.0110181177, 10.1907539167, 605.3820697090, 60.1973290710, 23.8046231705, 16.6085494134, - 3.6307172673, 1.6536810707, 4278.2949856871, 10.4280446109; + compare_compartments << 469.6949422507, 18.2547785569, 26.0031050672, 7.9621653957, 3.7350284051, 1.7375751388, + 4461.7342012617, 10.8782039239, 469.6949422507, 35.8652960540, 24.7473405359, 15.4322222059, 6.7061736912, + 2.9002347700, 4434.3772098068, 10.2765806855, 587.1186778134, 33.9201025102, 31.8752920696, 14.6856309847, + 6.6453275434, 2.9575629755, 4311.2552183628, 11.5421877404; ASSERT_EQ(compare_compartments.size(), static_cast(compartments.get_last_value().size())); @@ -164,11 +164,11 @@ TEST(TestIdeAgeres, compareWithPreviousRun) mio::TimeSeries transitions = sim.get_transitions(); Eigen::VectorX compare_transitions(num_transitions * num_agegroups); - compare_transitions << 31.5370062111, 30.6497959470, 14.1231866958, 14.7543908776, 6.6982921386, 6.6982921386, - 3.1606794140, 3.1606794140, 1.4742153411, 1.4742153411, 31.5370062111, 29.5087817552, 14.7543908776, - 14.1231866958, 6.3213588280, 6.3213588280, 2.9484306823, 2.9484306823, 1.3533839877, 1.3533839877, - 39.4212577639, 30.1092463410, 14.4459531059, 14.4459531059, 6.5114345346, 6.5114345346, 3.0573621415, - 3.0573621415, 1.4155355888, 1.4155355888; + compare_transitions << 36.5095571138, 35.2210349942, 15.9243307913, 16.7851751402, 7.4700568103, 7.4700568103, + 3.4751502776, 3.4751502776, 1.5959804568, 1.5959804568, 36.5095571138, 33.5703502804, 15.9243307913, + 14.9401136206, 6.9503005552, 6.9503005552, 3.0041079341, 3.0041079341, 1.3063243113, 1.3063243113, + 45.6369463923, 43.0480501661, 19.8862347266, 19.8862347266, 9.0259862757, 9.0259862757, 4.0260421953, + 4.0260421953, 1.7699583766, 1.7699583766; ASSERT_EQ(compare_transitions.size(), static_cast(transitions.get_last_value().size())); @@ -176,3 +176,229 @@ TEST(TestIdeAgeres, compareWithPreviousRun) ASSERT_NEAR(transitions.get_last_value()[j], compare_transitions[j], 1e-7); } } + +// Check results of a simulation with two age groups with an example calculated by hand, +// see https://doi.org/10.1016/j.amc.2025.129636 for the used formulas. +// The idea of this test is to mimic the checkSimulationFunctions test for the nonageresolved IDE-SECIR model. +// For this, we set up two models with two age groups, respectively. All parameters are set as in the non-ageresolved +// test. In the first model, we only have individuals in the first age group where the population size as well as the +// initial values for the deaths and transitions are chosen as in the non-ageresolved test. Hence, we can use the +// values from the other test to check for correctness. In the second model, we only have individuals in the second +// age group and the population size as well as the initial values for the deaths and transitions are scaled by some +// factor. This way, the expected results for this model are determined by the results from the other test and the +// scaling factor. With this, we can check that the correct indices are used for each age group. +TEST(IdeSecirAgeres, checkSimulationFunctions) +{ + using Vec = mio::TimeSeries::Vector; + size_t num_agegroups = 2; + ScalarType tmax = 0.5; + ScalarType dt = 0.5; + + // Define baseline values for the population size, the number of deaths and values for the transitions + // from SusceptibleToExposed and InfectedNoSymptomsToInfectedSymptoms for the considered age group. + // These values correspond to the values of the non-ageresolved checkSimulationFunctions test. + ScalarType population_baseline = 10000.; + ScalarType deaths_baseline = 10.; + ScalarType susceptibleToExposed_baseline = 1.; + ScalarType infectedNoSymptomsToInfectedSymptoms_baseline = 8.; + + // Scaling factor that determines by which factor the baseline values for the initialization will be scaled for + // the second age group. + ScalarType baseline_scaling = 2.; + + /*****************************************************************************************************************/ + + // Model with whole population in first age group with index 0. Use baseline values. + int considered_group = 0; + std::vector N_vec_group0 = {population_baseline, 0}; + std::vector deaths_vec_group0 = {deaths_baseline, 0}; + mio::CustomIndexArray N_group0 = mio::CustomIndexArray( + mio::AgeGroup(num_agegroups), N_vec_group0.begin(), N_vec_group0.end()); + mio::CustomIndexArray deaths_group0 = mio::CustomIndexArray( + mio::AgeGroup(num_agegroups), deaths_vec_group0.begin(), deaths_vec_group0.end()); + + // Create TimeSeries with num_transitions elements where transitions needed for simulation will be stored. + int num_transitions = (int)mio::isecir::InfectionTransition::Count; + mio::TimeSeries init_group0(num_transitions * num_agegroups); + + // Add time points for initialization for transitions and death. + Vec vec_init_group0 = Vec::Constant(num_transitions * num_agegroups, 0.); + // Set values for vec_init for group 0 as in nonageresolved test. + vec_init_group0[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::SusceptibleToExposed] = susceptibleToExposed_baseline; + vec_init_group0[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = + infectedNoSymptomsToInfectedSymptoms_baseline; + + // Add initial time point to TimeSeries. + init_group0.add_time_point(-0.5, vec_init_group0); + while (init_group0.get_last_time() < 0) { + init_group0.add_time_point(init_group0.get_last_time() + dt, vec_init_group0); + } + + // Initialize model. + mio::isecir::Model model_group0(std::move(init_group0), N_group0, deaths_group0, num_agegroups); + + // Set working parameters. + for (size_t group = 0; group < num_agegroups; group++) { + // In this example, SmootherCosine with parameter 1 (and thus with a maximum support of 1) + // is used for all TransitionDistribution%s for all groups. + mio::SmootherCosine smoothcos(1.0); + mio::StateAgeFunctionWrapper delaydistribution(smoothcos); + std::vector> vec_delaydistrib(num_transitions, delaydistribution); + model_group0.parameters.get()[mio::AgeGroup(group)] = vec_delaydistrib; + + std::vector vec_prob((int)mio::isecir::InfectionTransition::Count, 0.5); + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] = 1; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] = 1; + model_group0.parameters.get()[mio::AgeGroup(group)] = vec_prob; + + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, num_agegroups); + contact_matrix[0] = + mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, 4.)); + model_group0.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + mio::SmootherCosine smoothcos_prob(1.0); + mio::StateAgeFunctionWrapper prob(smoothcos_prob); + model_group0.parameters.get()[mio::AgeGroup(group)] = prob; + model_group0.parameters.get()[mio::AgeGroup(group)] = prob; + model_group0.parameters.get()[mio::AgeGroup(group)] = prob; + } + model_group0.parameters.set(0.); + model_group0.parameters.set(0); + + // Carry out simulation. + mio::isecir::Simulation sim_group0(model_group0, dt); + sim_group0.advance(tmax); + mio::TimeSeries secihurd_simulated_group0 = sim_group0.get_result(); + mio::TimeSeries transitions_simulated_group0 = sim_group0.get_transitions(); + + /*****************************************************************************************************************/ + + // Model with whole population in second age group with index 1. Use baseline values scaled by baseline_scaling. + considered_group = 1; + std::vector N_vec_group1 = {0., baseline_scaling * population_baseline}; + std::vector deaths_vec_group1 = {0., baseline_scaling * deaths_baseline}; + mio::CustomIndexArray N_group1 = mio::CustomIndexArray( + mio::AgeGroup(num_agegroups), N_vec_group1.begin(), N_vec_group1.end()); + mio::CustomIndexArray deaths_group1 = mio::CustomIndexArray( + mio::AgeGroup(num_agegroups), deaths_vec_group1.begin(), deaths_vec_group1.end()); + + // Create TimeSeries with num_transitions elements where transitions needed for simulation will be stored. + mio::TimeSeries init_group1(num_transitions * num_agegroups); + + // Add time points for initialization for transitions and death. + Vec vec_init_group1 = Vec::Constant(num_transitions * num_agegroups, 0.); + // Set values for vec_init for group 0 as in nonageresolved test. + vec_init_group1[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::SusceptibleToExposed] = + baseline_scaling * susceptibleToExposed_baseline; + vec_init_group1[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = + baseline_scaling * infectedNoSymptomsToInfectedSymptoms_baseline; + + // Add initial time point to TimeSeries. + init_group1.add_time_point(-0.5, vec_init_group1); + while (init_group1.get_last_time() < 0) { + init_group1.add_time_point(init_group1.get_last_time() + dt, vec_init_group1); + } + + // Initialize model. + mio::isecir::Model model_group1(std::move(init_group1), N_group1, deaths_group1, num_agegroups); + + // Set working parameters. + for (size_t group = 0; group < num_agegroups; group++) { + // In this example, SmootherCosine with parameter 1 (and thus with a maximum support of 1) + // is used for all TransitionDistribution%s for all groups. + mio::SmootherCosine smoothcos(1.0); + mio::StateAgeFunctionWrapper delaydistribution(smoothcos); + std::vector> vec_delaydistrib(num_transitions, delaydistribution); + model_group1.parameters.get()[mio::AgeGroup(group)] = vec_delaydistrib; + + std::vector vec_prob((int)mio::isecir::InfectionTransition::Count, 0.5); + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] = 1; + vec_prob[Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] = 1; + model_group1.parameters.get()[mio::AgeGroup(group)] = vec_prob; + + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, num_agegroups); + contact_matrix[0] = + mio::ContactMatrix(Eigen::MatrixX::Constant(num_agegroups, num_agegroups, 4.)); + model_group1.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + mio::SmootherCosine smoothcos_prob(1.0); + mio::StateAgeFunctionWrapper prob(smoothcos_prob); + model_group1.parameters.get()[mio::AgeGroup(group)] = prob; + model_group1.parameters.get()[mio::AgeGroup(group)] = prob; + model_group1.parameters.get()[mio::AgeGroup(group)] = prob; + } + model_group1.parameters.set(0.); + model_group1.parameters.set(0); + + // Carry out simulation. + mio::isecir::Simulation sim_group1(model_group1, dt); + sim_group1.advance(tmax); + mio::TimeSeries secihurd_simulated_group1 = sim_group1.get_result(); + mio::TimeSeries transitions_simulated_group1 = sim_group1.get_transitions(); + + /*****************************************************************************************************************/ + + // Define vectors for compartments and transitions with expected results from example + // (calculated by hand, see https://doi.org/10.1016/j.amc.2025.129636 for the used formulas). + std::vector secihurd_t0_baseline = {4995., 0.5, 0., 4., 0., 0., 4990.5, 10.}; + std::vector secihurd_t1_baseline = {4994.00020016, 0.49989992, 0.49994996, 0.12498749, + 1.03124687, 0.25781172, 4993.45699802, 10.12890586}; + std::vector transitions_t1_baseline = {0.99979984, 0.99989992, 0.24997498, 0.24997498, 2.06249374, + 2.06249374, 0.51562344, 0.51562344, 0.12890586, 0.12890586}; + + // Define vectors with expected results for first model. + int num_compartments = (int)mio::isecir::InfectionState::Count; + std::vector secihurd_t0_group0(num_compartments * num_agegroups, 0.); + std::vector secihurd_t1_group0(num_compartments * num_agegroups, 0.); + std::vector transitions_t1_group0(num_transitions * num_agegroups, 0.); + + // Copy expected results for first age group. + std::copy(secihurd_t0_baseline.begin(), secihurd_t0_baseline.end(), secihurd_t0_group0.begin()); + std::copy(secihurd_t1_baseline.begin(), secihurd_t1_baseline.end(), secihurd_t1_group0.begin()); + std::copy(transitions_t1_baseline.begin(), transitions_t1_baseline.end(), transitions_t1_group0.begin()); + + // Define vectors with expected results for second model. + std::vector secihurd_t0_group1(num_compartments * num_agegroups, 0.); + std::vector secihurd_t1_group1(num_compartments * num_agegroups, 0.); + std::vector transitions_t1_group1(num_transitions * num_agegroups, 0.); + + // Scale baseline vectors by baseline_scaling. + std::transform(secihurd_t0_baseline.begin(), secihurd_t0_baseline.end(), secihurd_t0_baseline.begin(), + [baseline_scaling](ScalarType x) { + return baseline_scaling * x; + }); + std::transform(secihurd_t1_baseline.begin(), secihurd_t1_baseline.end(), secihurd_t1_baseline.begin(), + [baseline_scaling](ScalarType x) { + return baseline_scaling * x; + }); + std::transform(transitions_t1_baseline.begin(), transitions_t1_baseline.end(), transitions_t1_baseline.begin(), + [baseline_scaling](ScalarType x) { + return baseline_scaling * x; + }); + + // Copy expected results for second age group. + std::copy(secihurd_t0_baseline.begin(), secihurd_t0_baseline.end(), secihurd_t0_group1.begin() + num_compartments); + std::copy(secihurd_t1_baseline.begin(), secihurd_t1_baseline.end(), secihurd_t1_group1.begin() + num_compartments); + std::copy(transitions_t1_baseline.begin(), transitions_t1_baseline.end(), + transitions_t1_group1.begin() + num_transitions); + + // Compare simulated compartments at time points t0 and t1. + for (Eigen::Index i = 0; i < (Eigen::Index)secihurd_simulated_group0.get_num_elements(); i++) { + EXPECT_NEAR(secihurd_simulated_group0[0][i], secihurd_t0_group0[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_group0[1][i], secihurd_t1_group0[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_group1[0][i], secihurd_t0_group1[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_group1[1][i], secihurd_t1_group1[i], 1e-8); + } + + // Compare simulated transitions with expected results at time point t1. + for (Eigen::Index i = 0; i < transitions_simulated_group0.get_num_elements(); i++) { + EXPECT_NEAR(transitions_simulated_group0[transitions_simulated_group0.get_num_time_points() - 1][i], + transitions_t1_group0[i], 1e-8); + EXPECT_NEAR(transitions_simulated_group1[transitions_simulated_group1.get_num_time_points() - 1][i], + transitions_t1_group1[i], 1e-8); + } +}