diff --git a/cpp/models/ide_secir/model.cpp b/cpp/models/ide_secir/model.cpp index fb2f8c3aff..6f7664d45c 100644 --- a/cpp/models/ide_secir/model.cpp +++ b/cpp/models/ide_secir/model.cpp @@ -267,11 +267,14 @@ 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; } @@ -279,7 +282,7 @@ void Model::initial_compute_compartments(ScalarType dt) 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; } @@ -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]; } } @@ -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) { @@ -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 { @@ -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 { @@ -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()[group][idx_InfectionTransitions] * sum; } diff --git a/cpp/models/ide_secir/simulation.cpp b/cpp/models/ide_secir/simulation.cpp index 366ff98315..db14a921f6 100644 --- a/cpp/models/ide_secir/simulation.cpp +++ b/cpp/models/ide_secir/simulation.cpp @@ -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: 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..67408945e5 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,166 @@ 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, where the first group is as in the nonageresolved +// test and the second group is a multiple of the first group, but the proportion of infected and dead individuals +// etc. is the same as in the first group. Both models have the same initial conditions but differ in their contact +// matrices. In the first model, individuals are only interacting with individuals in their own group, i.e. we only have +// intra-group contacts. In the second model, individuals are only interacting with individuals of the other group, +// i.e. we only have inter-group contacts. All remaining parameters are set as in the non-ageresolved test. +// This way, the expected results for both models are determined by the results from the nonageresolved test and the +// scaling factor. With this, we can check that the correct indices are used for each age group, both when individuals +// are interacting only with their own group or with another group. +TEST(TestIdeAgeres, 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.; + + // ***Set up models*** + std::vector N_vec = {population_baseline, baseline_scaling * population_baseline}; + std::vector deaths_vec = {deaths_baseline, baseline_scaling * deaths_baseline}; + mio::CustomIndexArray N = + mio::CustomIndexArray(mio::AgeGroup(num_agegroups), N_vec.begin(), N_vec.end()); + mio::CustomIndexArray deaths = mio::CustomIndexArray( + mio::AgeGroup(num_agegroups), deaths_vec.begin(), deaths_vec.end()); + + // Create TimeSeries with num_transitions elements where transitions needed for simulation will be stored. + size_t num_transitions = (size_t)mio::isecir::InfectionTransition::Count; + mio::TimeSeries init_intra(num_transitions * num_agegroups); + + // Add time points for transitions for initialization. + Vec vec_init = Vec::Constant(num_transitions * num_agegroups, 0.); + // Set values for vec_init for group 0 with baseline values. + int considered_group = 0; + vec_init[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::SusceptibleToExposed] = susceptibleToExposed_baseline; + vec_init[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = + infectedNoSymptomsToInfectedSymptoms_baseline; + // Set values for vec_init for group 1 with baseline values scaled by baseline_scaling. + considered_group = 1; + vec_init[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::SusceptibleToExposed] = + baseline_scaling * susceptibleToExposed_baseline; + vec_init[considered_group * (int)mio::isecir::InfectionTransition::Count + + (int)mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = + baseline_scaling * infectedNoSymptomsToInfectedSymptoms_baseline; + + // Add time points to TimeSeries. + init_intra.add_time_point(-0.5, vec_init); + while (init_intra.get_last_time() < 0) { + init_intra.add_time_point(init_intra.get_last_time() + dt, vec_init); + } + + // Initialize model. + mio::isecir::Model model_intra(std::move(init_intra), N, deaths, num_agegroups); + + // Set working parameters for both models. + 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_intra.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_intra.parameters.get()[mio::AgeGroup(group)] = vec_prob; + + mio::SmootherCosine smoothcos_prob(1.0); + mio::StateAgeFunctionWrapper prob(smoothcos_prob); + model_intra.parameters.get()[mio::AgeGroup(group)] = prob; + model_intra.parameters.get()[mio::AgeGroup(group)] = prob; + model_intra.parameters.get()[mio::AgeGroup(group)] = prob; + } + model_intra.parameters.set(0.); + model_intra.parameters.set(0); + + // Define model_inter as copy of model_intra, i.e. all initial conditions and parameters are the same in both models. + mio::isecir::Model model_inter = model_intra; + + // Here we set the contact matrices for both models which is where they differ from each other. + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, num_agegroups); + // With this contact matrix, individuals are only interacting with individuals in their own group. + Eigen::MatrixX matrix_intra{{4., 0.}, {0., 4.}}; + contact_matrix[0] = mio::ContactMatrix(matrix_intra); + model_intra.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + // With this contact matrix, individuals are only interacting with individuals of the other group. + Eigen::MatrixX matrix_inter{{0., 4.}, {4., 0.}}; + contact_matrix[0] = mio::ContactMatrix(matrix_inter); + model_inter.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + // ***Simulate*** + mio::isecir::Simulation sim_intra(model_intra, dt); + sim_intra.advance(tmax); + mio::TimeSeries secihurd_simulated_intra = sim_intra.get_result(); + mio::TimeSeries transitions_simulated_intra = sim_intra.get_transitions(); + + mio::isecir::Simulation sim_inter(model_inter, dt); + sim_inter.advance(tmax); + mio::TimeSeries secihurd_simulated_inter = sim_inter.get_result(); + mio::TimeSeries transitions_simulated_inter = sim_inter.get_transitions(); + + // ***Test for correctness*** + // 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}; + + ASSERT_EQ(secihurd_t0_baseline.size() * num_agegroups, secihurd_simulated_intra.get_num_elements()); + // Compare simulated compartments at time points t0 and t1. + // In both models, the two age groups have the same number of contacts where the only difference lies in whom they + // have contact with. In both age groups, the probability of meeting an infected individualis the same, hence we + // expect the same results in both models. + for (size_t i = 0; i < secihurd_t0_baseline.size(); i++) { + const size_t age1 = i; + const size_t age2 = i + secihurd_t0_baseline.size(); + // Test model_intra. + EXPECT_NEAR(secihurd_simulated_intra[0][age1], secihurd_t0_baseline[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_intra[0][age2], baseline_scaling * secihurd_t0_baseline[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_intra[1][age1], secihurd_t1_baseline[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_intra[1][age2], baseline_scaling * secihurd_t1_baseline[i], 1e-8); + // Test model_inter. + EXPECT_NEAR(secihurd_simulated_inter[0][age1], secihurd_t0_baseline[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_inter[0][age2], baseline_scaling * secihurd_t0_baseline[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_inter[1][age1], secihurd_t1_baseline[i], 1e-8); + EXPECT_NEAR(secihurd_simulated_inter[1][age2], baseline_scaling * secihurd_t1_baseline[i], 1e-8); + } + + ASSERT_EQ(transitions_t1_baseline.size() * num_agegroups, transitions_simulated_intra.get_num_elements()); + // Compare simulated transitions with expected results at time point t1. + for (size_t i = 0; i < transitions_t1_baseline.size(); i++) { + const size_t age1 = i; + const size_t age2 = i + transitions_t1_baseline.size(); + // Test model_intra. + EXPECT_NEAR(transitions_simulated_intra.get_last_value()[age1], transitions_t1_baseline[i], 1e-8); + EXPECT_NEAR(transitions_simulated_intra.get_last_value()[age2], baseline_scaling * transitions_t1_baseline[i], + 1e-8); + // Test model_inter. + EXPECT_NEAR(transitions_simulated_inter.get_last_value()[age1], transitions_t1_baseline[i], 1e-8); + EXPECT_NEAR(transitions_simulated_inter.get_last_value()[age2], baseline_scaling * transitions_t1_baseline[i], + 1e-8); + } +}