diff --git a/cpp/models/ide_secir/model.cpp b/cpp/models/ide_secir/model.cpp index 733bd1933d..ce24ae9c13 100644 --- a/cpp/models/ide_secir/model.cpp +++ b/cpp/models/ide_secir/model.cpp @@ -24,6 +24,8 @@ #include "memilio/utils/logging.h" #include "memilio/math/eigen.h" +#include "vector" + namespace mio { namespace isecir @@ -37,8 +39,6 @@ Model::Model(TimeSeries&& init, ScalarType N_init, ScalarType deaths , m_N{N_init} , m_total_confirmed_cases{total_confirmed_cases} { - m_deaths_before = - deaths - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)]; // Add first time point in m_populations according to last time point in m_transitions which is where we start the simulation. m_populations.add_time_point( m_transitions.get_last_time(), TimeSeries::Vector::Constant((int)InfectionState::Count, 0)); @@ -46,11 +46,72 @@ Model::Model(TimeSeries&& init, ScalarType N_init, ScalarType deaths m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] = deaths; } -void Model::initialize(ScalarType dt) +// ---- Functionality to calculate the sizes of the compartments for time t0. ---- +void Model::compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_InfectionState, + Eigen::Index idx_IncomingFlow, int idx_TransitionDistribution1, + int idx_TransitionDistribution2) +{ + ScalarType sum = 0; + + // Determine relevant calculation area and corresponding index. + ScalarType calc_time = + std::max(parameters.get()[idx_TransitionDistribution1].get_support_max(dt, m_tol), + parameters.get()[idx_TransitionDistribution2].get_support_max(dt, m_tol)); + + Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1; + + Eigen::Index num_time_points = m_transitions.get_num_time_points(); + + for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) { + + ScalarType state_age = (num_time_points - 1 - i) * dt; + + sum += (parameters.get()[idx_TransitionDistribution1] * + parameters.get()[idx_TransitionDistribution1].eval(state_age) + + (1 - parameters.get()[idx_TransitionDistribution1]) * + parameters.get()[idx_TransitionDistribution2].eval(state_age)) * + m_transitions[i + 1][idx_IncomingFlow]; + } + + m_populations.get_last_value()[idx_InfectionState] = sum; +} + +void Model::initial_compute_compartments_infection(ScalarType dt) { + // Exposed + compute_compartment_from_flows(dt, Eigen::Index(InfectionState::Exposed), + Eigen::Index(InfectionTransition::SusceptibleToExposed), + (int)InfectionTransition::ExposedToInfectedNoSymptoms); + // InfectedNoSymptoms + compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedNoSymptoms), + Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), + (int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms, + (int)InfectionTransition::InfectedNoSymptomsToRecovered); + // InfectedSymptoms + compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedSymptoms), + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), + (int)InfectionTransition::InfectedSymptomsToInfectedSevere, + (int)InfectionTransition::InfectedSymptomsToRecovered); + // InfectedSevere + compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedSevere), + Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + (int)InfectionTransition::InfectedSevereToInfectedCritical, + (int)InfectionTransition::InfectedSevereToRecovered); + // InfectedCritical + compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedCritical), + Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + (int)InfectionTransition::InfectedCriticalToDead, + (int)InfectionTransition::InfectedCriticalToRecovered); +} + +void Model::initial_compute_compartments(ScalarType dt) +{ + // The initialization method only affects the Susceptible and Recovered compartments. + // It is possible to calculate the sizes of the other compartments in advance because only the initial values of the flows are used. + initial_compute_compartments_infection(dt); + if (m_total_confirmed_cases > 1e-12) { m_initialization_method = 1; - other_compartments_current_timestep(dt); // The scheme of the ODE model for initialization is applied here. m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] = @@ -71,8 +132,6 @@ void Model::initialize(ScalarType dt) else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) { // Take initialized value for Susceptibles if value can't be calculated via the standard formula. m_initialization_method = 2; - // Calculate other compartment sizes for start time t0. - other_compartments_current_timestep(dt); // R; need an initial value for R, therefore do not calculate via compute_recovered() m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] = @@ -87,9 +146,7 @@ void Model::initialize(ScalarType dt) else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) { // If value for Recovered is initialized and standard method is not applicable (i.e., no value for total infections // or Susceptibles given directly), calculate Susceptibles via other compartments. - // The calculation of other compartments' values is not dependent on Susceptibles at time t0, i.e., S(t0), but only on the transitions of the past. m_initialization_method = 3; - other_compartments_current_timestep(dt); m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] = m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - @@ -104,7 +161,7 @@ void Model::initialize(ScalarType dt) // Compute Susceptibles at t0 and m_forceofinfection at time t0-dt as initial values for discretization scheme. // Use m_forceofinfection at t0-dt to be consistent with further calculations of S (see compute_susceptibles()), // where also the value of m_forceofinfection for the previous timestep is used. - update_forceofinfection(dt, true); + compute_forceofinfection(dt, true); if (m_forceofinfection > 1e-12) { m_initialization_method = 4; @@ -117,10 +174,7 @@ void Model::initialize(ScalarType dt) m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] / (dt * m_forceofinfection); - // Calculate other compartment sizes for t0. - other_compartments_current_timestep(dt); - - // R; need an initial value for R, therefore do not calculate via compute_recovered(). + // Recovered; calculated as the difference between the total population and the sum of the other compartment sizes. m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] = m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - @@ -151,14 +205,15 @@ void Model::initialize(ScalarType dt) } // Compute m_forceofinfection at time t0 needed for further simulation. - update_forceofinfection(dt); + compute_forceofinfection(dt); } +// ---- Functionality for the iterations of a simulation. ---- void Model::compute_susceptibles(ScalarType dt) { Eigen::Index num_time_points = m_populations.get_num_time_points(); // Using number of Susceptibles from previous time step and force of infection from previous time step: - // compute current number of Susceptibles and store Susceptibles in m_populations + // Compute current number of Susceptibles and store Susceptibles in m_populations. m_populations.get_last_value()[Eigen::Index(InfectionState::Susceptible)] = m_populations[num_time_points - 2][Eigen::Index(InfectionState::Susceptible)] / (1 + dt * m_forceofinfection); } @@ -179,7 +234,7 @@ void Model::compute_flow(int idx_InfectionTransitions, Eigen::Index idx_Incoming Eigen::Index num_time_points = m_transitions.get_num_time_points(); for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) { - // (num_time_points - 1 - i)* dt is the time, the individuals has already spent in this state. + // (num_time_points - 1 - i) * dt is the time, the individuals has already spent in this state. ScalarType state_age = (num_time_points - 1 - i) * dt; @@ -195,49 +250,94 @@ void Model::compute_flow(int idx_InfectionTransitions, Eigen::Index idx_Incoming void Model::flows_current_timestep(ScalarType dt) { - // Calculate flow from S to E with force of infection from previous time step und susceptibles from current time step. + // Calculate flow SusceptibleToExposed with the force of infection from previous time step und susceptibles from current time step. m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] = dt * m_forceofinfection * m_populations.get_last_value()[Eigen::Index(InfectionState::Susceptible)]; + // Calculate all other flows with compute_flow. - // Flow from E to C + // Flow ExposedToInfectedNoSymptoms. compute_flow((int)InfectionTransition::ExposedToInfectedNoSymptoms, Eigen::Index(InfectionTransition::SusceptibleToExposed), dt); - // Flow from C to I + // Flow InfectedNoSymptomsToInfectedSymptoms. compute_flow((int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms, Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt); - // Flow from C to R + // Flow InfectedNoSymptomsToRecovered. compute_flow((int)InfectionTransition::InfectedNoSymptomsToRecovered, Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt); - // Flow from I to H + // Flow InfectedSymptomsToInfectedSevere. compute_flow((int)InfectionTransition::InfectedSymptomsToInfectedSevere, Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt); - // Flow from I to R + // Flow InfectedSymptomsToRecovered. compute_flow((int)InfectionTransition::InfectedSymptomsToRecovered, Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt); - // Flow from H to U + // Flow InfectedSevereToInfectedCritical. compute_flow((int)InfectionTransition::InfectedSevereToInfectedCritical, Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt); - // Flow from to H to R + // Flow InfectedSevereToRecovered. compute_flow((int)InfectionTransition::InfectedSevereToRecovered, Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt); - // Flow from U to D + // Flow InfectedCriticalToDead. compute_flow((int)InfectionTransition::InfectedCriticalToDead, Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt); - // Flow from U to R + // Flow InfectedCriticalToRecovered. compute_flow((int)InfectionTransition::InfectedCriticalToRecovered, Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt); } -void Model::compute_deaths() +void Model::update_compartments() { - Eigen::Index num_time_points = m_populations.get_num_time_points(); + // Exposed + update_compartment_from_flow(InfectionState::Exposed, {InfectionTransition::SusceptibleToExposed}, + {InfectionTransition::ExposedToInfectedNoSymptoms}); + + // InfectedNoSymptoms + update_compartment_from_flow(InfectionState::InfectedNoSymptoms, {InfectionTransition::ExposedToInfectedNoSymptoms}, + {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms, + InfectionTransition::InfectedNoSymptomsToRecovered}); + + // InfectedSymptoms + update_compartment_from_flow( + InfectionState::InfectedSymptoms, {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms}, + {InfectionTransition::InfectedSymptomsToInfectedSevere, InfectionTransition::InfectedSymptomsToRecovered}); + + // InfectedSevere + update_compartment_from_flow( + InfectionState::InfectedSevere, {InfectionTransition::InfectedSymptomsToInfectedSevere}, + {InfectionTransition::InfectedSevereToInfectedCritical, InfectionTransition::InfectedSevereToRecovered}); + + // InfectedCritical + update_compartment_from_flow( + InfectionState::InfectedCritical, {InfectionTransition::InfectedSevereToInfectedCritical}, + {InfectionTransition::InfectedCriticalToDead, InfectionTransition::InfectedCriticalToRecovered}); + + // Recovered + update_compartment_from_flow( + InfectionState::Recovered, + {InfectionTransition::InfectedNoSymptomsToRecovered, InfectionTransition::InfectedSymptomsToRecovered, + InfectionTransition::InfectedSevereToRecovered, InfectionTransition::InfectedCriticalToRecovered}, + std::vector()); + + // Dead + update_compartment_from_flow(InfectionState::Dead, {InfectionTransition::InfectedCriticalToDead}, + std::vector()); +} - m_populations.get_last_value()[Eigen::Index(InfectionState::Dead)] = - m_populations[num_time_points - 2][Eigen::Index(InfectionState::Dead)] + - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)]; +void Model::update_compartment_from_flow(InfectionState infectionState, + std::vector const& IncomingFlows, + std::vector const& OutgoingFlows) +{ + Eigen::Index num_time_points = m_populations.get_num_time_points(); + ScalarType updated_compartment = m_populations[num_time_points - 2][Eigen::Index(infectionState)]; + for (const InfectionTransition& inflow : IncomingFlows) { + updated_compartment += m_transitions.get_last_value()[Eigen::Index(inflow)]; + } + for (const InfectionTransition& outflow : OutgoingFlows) { + updated_compartment -= m_transitions.get_last_value()[Eigen::Index(outflow)]; + } + m_populations.get_last_value()[Eigen::Index(infectionState)] = updated_compartment; } -void Model::update_forceofinfection(ScalarType dt, bool initialization) +void Model::compute_forceofinfection(ScalarType dt, bool initialization) { m_forceofinfection = 0; @@ -253,8 +353,8 @@ void Model::update_forceofinfection(ScalarType dt, bool initialization) .get_support_max(dt, m_tol)}); // Corresponding index. - /* need calc_time_index timesteps in sum, - subtract 1 because in the last summand all TransitionDistributions evaluate to 0 (by definition of support_max)*/ + // Need calc_time_index timesteps in sum, + // subtract 1 because in the last summand all TransitionDistributions evaluate to 0 (by definition of support_max). Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1; Eigen::Index num_time_points; @@ -266,7 +366,9 @@ void Model::update_forceofinfection(ScalarType dt, bool initialization) num_time_points = m_transitions.get_num_time_points() - 1; // Get time of penultimate timepoint in m_transitions. current_time = m_transitions.get_time(num_time_points - 1); - deaths = m_deaths_before; + // Determine the number of death at time t0-dt. + deaths = m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] - + m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)]; } else { // Determine m_forceofinfection for current last time in m_transitions. @@ -307,74 +409,5 @@ void Model::update_forceofinfection(ScalarType dt, bool initialization) m_forceofinfection = 1 / (m_N - deaths) * m_forceofinfection; } -void Model::compute_compartment(Eigen::Index idx_InfectionState, Eigen::Index idx_IncomingFlow, - int idx_TransitionDistribution1, int idx_TransitionDistribution2, ScalarType dt) -{ - ScalarType sum = 0; - - // Determine relevant calculation area and corresponding index. - ScalarType calc_time = - std::max(parameters.get()[idx_TransitionDistribution1].get_support_max(dt, m_tol), - parameters.get()[idx_TransitionDistribution2].get_support_max(dt, m_tol)); - - Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1; - - Eigen::Index num_time_points = m_transitions.get_num_time_points(); - - for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) { - - ScalarType state_age = (num_time_points - 1 - i) * dt; - - sum += (parameters.get()[idx_TransitionDistribution1] * - parameters.get()[idx_TransitionDistribution1].eval(state_age) + - (1 - parameters.get()[idx_TransitionDistribution1]) * - parameters.get()[idx_TransitionDistribution2].eval(state_age)) * - m_transitions[i + 1][idx_IncomingFlow]; - } - - m_populations.get_last_value()[idx_InfectionState] = sum; -} - -void Model::other_compartments_current_timestep(ScalarType dt) -{ - // E - // idx_TransitionDistribution2 is a dummy index as there is no transition from E to R in our model, - // write any transition here as probability of going from E to R is 0. - compute_compartment(Eigen::Index(InfectionState::Exposed), Eigen::Index(InfectionTransition::SusceptibleToExposed), - (int)InfectionTransition::ExposedToInfectedNoSymptoms, 0, dt); - // C - compute_compartment(Eigen::Index(InfectionState::InfectedNoSymptoms), - Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), - (int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms, - (int)InfectionTransition::InfectedNoSymptomsToRecovered, dt); - // I - compute_compartment(Eigen::Index(InfectionState::InfectedSymptoms), - Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), - (int)InfectionTransition::InfectedSymptomsToInfectedSevere, - (int)InfectionTransition::InfectedSymptomsToRecovered, dt); - // H - compute_compartment(Eigen::Index(InfectionState::InfectedSevere), - Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), - (int)InfectionTransition::InfectedSevereToInfectedCritical, - (int)InfectionTransition::InfectedSevereToRecovered, dt); - // U - compute_compartment(Eigen::Index(InfectionState::InfectedCritical), - Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), - (int)InfectionTransition::InfectedCriticalToDead, - (int)InfectionTransition::InfectedCriticalToRecovered, dt); -} - -void Model::compute_recovered() -{ - Eigen::Index num_time_points = m_populations.get_num_time_points(); - - m_populations.get_last_value()[Eigen::Index(InfectionState::Recovered)] = - m_populations[num_time_points - 2][Eigen::Index(InfectionState::Recovered)] + - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered)] + - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered)] + - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedSevereToRecovered)] + - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToRecovered)]; -} - } // namespace isecir } // namespace mio diff --git a/cpp/models/ide_secir/model.h b/cpp/models/ide_secir/model.h index 46eed025ca..2e3d3b1ea1 100644 --- a/cpp/models/ide_secir/model.h +++ b/cpp/models/ide_secir/model.h @@ -21,11 +21,14 @@ #ifndef IDESECIR_MODEL_H #define IDESECIR_MODEL_H +#include "abm/virus_variant.h" #include "ide_secir/parameters.h" #include "ide_secir/infection_state.h" #include "memilio/config.h" #include "memilio/utils/time_series.h" +#include "vector" + namespace mio { namespace isecir @@ -36,12 +39,13 @@ class Model public: /** - * @brief Constructor to create an IDE SECIR model. + * @brief Constructor to create an IDE-SECIR model. * * @param[in, out] init TimeSeries with the initial values of the number of individuals, * which transit within one timestep dt from one compartment to another. * Possible transitions are specified in #InfectionTransition%s. - * Considered time points should have the distance dt. The last time point determines the start time t0 of the simulation. + * Considered time points should have the distance dt. The last time point determines the start time t0 of the + * simulation. * The time history must reach a certain point in the past so that the simulation can be performed. * A warning is displayed if the condition is violated. * @param[in] N_init The population of the considered region. @@ -52,6 +56,118 @@ class Model Model(TimeSeries&& init, ScalarType N_init, ScalarType deaths, ScalarType total_confirmed_cases = 0, const ParameterSet& Parameterset_init = ParameterSet()); + // ---- Functionality to calculate the sizes of the compartments for time t0. ---- + /** + * @brief Compute the compartment specified in idx_InfectionState at the current time -- only using historic flow + * values and disrespecting potential, previous compartment value. + * + * The computation is meaningful for all compartments except Susceptible, Recovered and #Death + * and mostly needed for initialization. + * For Susceptible, Recovered and Dead, use corresponding alternative functions. + * + * @param[in] dt Time discretization step size. + * @param[in] idx_InfectionState Specifies the considered #InfectionState + * @param[in] idx_IncomingFlow Specifies the index of the infoming flow to #InfectionState in m_transitions. + * @param[in] idx_TransitionDistribution1 Specifies the index of the first relevant TransitionDistribution, + * related to a flow from the considered #InfectionState to any other #InfectionState. + * This index is also used for related probability. + * @param[in] idx_TransitionDistribution2 Specifies the index of the second relevant TransitionDistribution, + * related to a flow from the considered #InfectionState to any other #InfectionState (in most cases + * to Recovered). + * Related probability is calculated via 1-probability[idx_TransitionDistribution1]. + * Sometimes the second index is not needed, e.g., if probability[idx_TransitionDistribution1]=1. + */ + void compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_InfectionState, Eigen::Index idx_IncomingFlow, + int idx_TransitionDistribution1, int idx_TransitionDistribution2 = 0); + + /** + * @brief Computes the values of the infection compartments subset at initialization. + * + * The values for the compartments Exposed, InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and InfectedCritical + * for time t_0 are calculated using the initial data in form of flows. + * Calculated values are stored in m_populations. + * + * @param[in] dt Time discretization step size. + */ + void initial_compute_compartments_infection(ScalarType dt); + + /** + * @brief Computes the values of compartments at initialization. + * + * The initialization method is selected automatically based on the different values that need to be set beforehand. + * Infection compartments are always computed through historic flows. + * Initialization methods for Susceptible and Recovered are tested in the following order: + * 1.) If a positive number for the total number of confirmed cases is set, Recovered is set according to that + * value and Susceptible%s are derived. + * 2.) If Susceptible%s are set, Recovered will be derived. + * 3.) If Recovered are set directly, Susceptible%s are derived. + * 4.) If none of the above is set with positive value, the force of infection is used as in Messina et al (2021) + * to set the Susceptible%s. + * + * The function calculate_compartments_initialization() is used in every method for the compartments Exposed, + * InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and InfectedCritical. + * In addition, the force of infection is calculated for start time t_0. + * + * @param[in] dt Time discretization step size. + */ + void initial_compute_compartments(ScalarType dt); + + // ---- Functionality for the iterations of a simulation. ---- + /** + * @brief Computes number of Susceptible%s for the current last time in m_populations. + * + * Number is computed using previous number of Susceptible%s and the force of infection (also from previous timestep). + * Number is stored at the matching index in m_populations. + * @param[in] dt Time discretization step size. + */ + void compute_susceptibles(ScalarType dt); + + /** + * @brief Computes size of a flow. + * + * Computes size of one flow from #InfectionTransition, specified in idx_InfectionTransitions, for the current + * last time value in m_transitions. + * + * @param[in] idx_InfectionTransitions Specifies the considered flow from #InfectionTransition. + * @param[in] idx_IncomingFlow Index of the flow in #InfectionTransition, which goes to the considered starting + * compartment of the flow specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] dt Time step to compute flow for. + */ + void compute_flow(int idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt); + + /** + * @brief Sets all required flows for the current last timestep in m_transitions. + * + * New values are stored in m_transitions. Most values are computed via the function compute_flow(). + * + * @param[in] dt Time step. + */ + void flows_current_timestep(ScalarType dt); + + /** + * @brief Updates the values of all compartments except Susceptible at initialization. + * + * New values are stored in m_populations. The values are calculated using the compartment size in the previous + * time step and the related flows of the current time step. + * Therefore the flows of the current time step should be calculated before using this function. + * + */ + void update_compartments(); + + /** + * @brief Computes force of infection for the current last time in m_transitions. + * + * Computed value is stored in m_forceofinfection. + * + * @param[in] dt Time discretization step size. + * @param[in] initialization If true we are in the case of the initialization of the model. + * For this we need forceofinfection at time point t0-dt and not at the current last time + * (given by m_transitions) as in the other time steps. + */ + void compute_forceofinfection(ScalarType dt, bool initialization = false); + + // ---- Additional functionality such as constraint checking, setters and getters, etc. ---- /** * @brief Checks constraints on model parameters and initial data. * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. @@ -105,111 +221,6 @@ class Model return parameters.check_constraints(); } - /** - * @brief Initializes the model and sets compartment values at start time t0. - * - * The initialization method is selected automatically based on the different values that need to be set beforehand. - * Infection compartments are always computed through historic flow. - * Initialization methods for #Susceptible and #Recovered are tested in the following order: - * 1.) If a positive number for the total number of confirmed cases is set, #Recovered is set according to that value and #Susceptible%s are derived. - * 2.) If #Susceptible%s are set, #Recovered will be derived. - * 3.) If #Recovered are set directly, #Susceptible%s are derived. - * 4.) If none of the above is set with positive value, the force of infection is used as in Messina et al (2021) to set the #Susceptible%s. - * - * @param[in] dt Time discretization step size. - */ - void initialize(ScalarType dt); - - /** - * @brief Computes number of #Susceptible%s for the current last time in m_populations. - * - * Number is computed using previous number of #Susceptible%s and the force of infection (also from previous timestep). - * Number is stored at the matching index in m_populations. - * @param[in] dt Time discretization step size. - */ - void compute_susceptibles(ScalarType dt); - - /** - * @brief Computes size of a flow. - * - * Computes size of one flow from #InfectionTransition, specified in idx_InfectionTransitions, for the current - * last time value in m_transitions. - * - * @param[in] idx_InfectionTransitions Specifies the considered flow from #InfectionTransition. - * @param[in] idx_IncomingFlow Index of the flow in #InfectionTransition, which goes to the considered starting - * compartment of the flow specified in idx_InfectionTransitions. Size of considered flow is calculated via - * the value of this incoming flow. - * @param[in] dt Time step to compute flow for. - */ - void compute_flow(int idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt); - - /** - * @brief Sets all required flows for the current last timestep in m_transitions. - * - * New values are stored in m_transitions. Most values are computed via the function compute_flow(). - * - * @param[in] dt Time step. - */ - void flows_current_timestep(ScalarType dt); - - /** - * @brief Computes total number of people in the #Dead compartment for the current last time in m_populations. - * - * Number is stored in m_populations. - * - */ - void compute_deaths(); - - /** - * @brief Computes force of infection for the current last time in m_transitions. - * - * Computed value is stored in m_forceofinfection. - * - * @param[in] dt Time discretization step size. - * @param[in] initialization If true we are in the case of the initialization of the model. - * For this we need forceofinfection at time point t0-dt and not at the current last time (given by m_transitions) - * as in the other time steps. - */ - void update_forceofinfection(ScalarType dt, bool initialization = false); - - /** - * @brief Get the size of the compartment specified in idx_InfectionState at the current last time in m_populations. - * - * Calculation is reasonable for all compartments except #Susceptible, #Recovered and #Dead. - * Therefore, we have alternative functions for those compartments. - * - * @param[in] idx_InfectionState Specifies the considered #InfectionState - * @param[in] idx_IncomingFlow Specifies the index of the infoming flow to #InfectionState in m_transitions. - * @param[in] idx_TransitionDistribution1 Specifies the index of the first relevant TransitionDistribution, - * related to a flow from the considered #InfectionState to any other #InfectionState. - * This index is also used for related probability. - * @param[in] idx_TransitionDistribution2 Specifies the index of the second relevant TransitionDistribution, - * related to a flow from the considered #InfectionState to any other #InfectionState (in most cases to #Recovered). - * Necessary related probability is calculated via 1-probability[idx_TransitionDistribution1]. - * If the second index is not needed, eg if probability[idx_TransitionDistribution1]=1, - * just use an arbitrary legal index. - * @param[in] dt Time discretization step size. - */ - void compute_compartment(Eigen::Index idx_InfectionState, Eigen::Index idx_IncomingFlow, - int idx_TransitionDistribution1, int idx_TransitionDistribution2, ScalarType dt); - - /** - * @brief Sets all values of remaining compartments (compartments apart from #Susceptible, #Recovered, #Dead) for the current last timestep in m_populations. - * - * New values are stored in m_populations. Most values are computed via the function get_size_of_compartments(). - * - * @param[in] dt Time discretization step size. - */ - void other_compartments_current_timestep(ScalarType dt); - - /** - * @brief Computes total number of #Recovered for the current last time in m_populations. - * - * Number is stored in m_populations. - * - */ - void compute_recovered(); - /** * @brief Setter for the tolerance used to calculate the maximum support of the TransitionDistributions. * @@ -221,36 +232,53 @@ class Model } /** - * @brief Returns the number associated with the method selected automatically for initialization. - * - * @returns 0 if the model has not yet been initialized and the method has not yet been selected, - * 1 if the method using the total number of confirmed cases at time t0 is used, - * 2 if the initialization is calculated using a beforehand set value for #Susceptible, - * 3 if the initialization is calculated using a beforehand set value for #Recovered, - * 4 if the force of infection method is used, - * -1 if the initialization was not possible using any of the methods and - * -2 if the initialization was possible using one of the provided methods but the result is not appropriate. - */ - int get_initialization_method() + * @brief Returns the index of the automatically selected initialization method. + * + * The initialization method is selected automatically based on the different values that need to be set beforehand. + * Infection compartments are always computed through historic flow. + * Initialization methods for Susceptible and Recovered are tested in the following order: + * 1.) If a positive number for the total number of confirmed cases is set, Recovered is set according to that value + * and Susceptible%s are derived. + * 2.) If Susceptible%s are set, Recovered will be derived. + * 3.) If Recovered are set directly, Susceptible%s are derived. + * 4.) If none of the above is set with positive value, the force of infection is used as in Messina et al (2021) to + * set the Susceptible%s. + * + * @return Index representing the initialization method. + */ + int get_initialization_method_compartments() { return m_initialization_method; } + // ---- Public parameters. ---- ParameterSet parameters{}; ///< ParameterSet of Model Parameters. - /* Attention: m_populations and m_transitions do not necessarily have the same number of time points due to the initialization part. */ + // Attention: m_populations and m_transitions do not necessarily have the same number of time points due to the + // initialization part. TimeSeries m_transitions; ///< TimeSeries containing points of time and the corresponding number of transitions. - TimeSeries - m_populations; ///< TimeSeries containing points of time and the corresponding number of people in defined #InfectionState%s. + TimeSeries m_populations; ///< TimeSeries containing points of time and the corresponding number of + // people in defined #InfectionState%s. private: + /** + * @brief Updates the values of one compartment using flows. + * + * New value is stored in m_populations. The value is calculated using the compartment size in the previous + * time step and the related flows of the current time step. + * Therefore the flows of the current time step should be calculated before using this function. + */ + void update_compartment_from_flow(InfectionState infectionState, + std::vector const& IncomingFlows, + std::vector const& OutgoingFlows); + + // ---- Private parameters. ---- ScalarType m_forceofinfection{0}; ///< Force of infection term needed for numerical scheme. ScalarType m_N{0}; ///< Total population size of the considered region. - ScalarType m_deaths_before{0}; ///< Total number of deaths at the time point t0 - dt. ScalarType m_total_confirmed_cases{0}; ///< Total number of confirmed cases at time t0. ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the TransitionDistributions. - int m_initialization_method{ - 0}; ///< Gives the index of the method used for the initialization of the model. See also get_initialization_method() for the number code. + int m_initialization_method{0}; ///< Gives the index of the method used for the initialization of the model. + //See also get_initialization_method_compartments() for the number code. }; } // namespace isecir diff --git a/cpp/models/ide_secir/simulation.cpp b/cpp/models/ide_secir/simulation.cpp index 99711a9e9b..bc235ed251 100644 --- a/cpp/models/ide_secir/simulation.cpp +++ b/cpp/models/ide_secir/simulation.cpp @@ -34,9 +34,9 @@ void Simulation::advance(ScalarType tmax) { mio::log_info("Simulating IDE-SECIR from t0 = {} until tmax = {} with dt = {}.", m_model->m_transitions.get_last_time(), tmax, m_dt); - m_model->initialize(m_dt); + m_model->initial_compute_compartments(m_dt); - // for every time step: + // For every time step: while (m_model->m_transitions.get_last_time() < tmax - m_dt / 2) { m_model->m_transitions.add_time_point(m_model->m_transitions.get_last_time() + m_dt); @@ -45,19 +45,14 @@ void Simulation::advance(ScalarType tmax) // compute Susceptibles: m_model->compute_susceptibles(m_dt); - // compute flows: + // Compute flows: m_model->flows_current_timestep(m_dt); - // compute Dead compartment: - m_model->compute_deaths(); + // Update remaining compartments: + m_model->update_compartments(); - // compute m_forceofinfection - // (only used for calculation of Susceptibles and the InfectionTransition SusceptibleToExposed in the next timestep!): - m_model->update_forceofinfection(m_dt); - - // compute remaining compartments from flows: - m_model->other_compartments_current_timestep(m_dt); - m_model->compute_recovered(); + // Compute m_forceofinfection (only used for calculation of Susceptibles and flow SusceptibleToExposed in the next timestep!): + m_model->compute_forceofinfection(m_dt); } } diff --git a/cpp/tests/test_binary_serializer.cpp b/cpp/tests/test_binary_serializer.cpp index 16c5ef7d7d..dfc9f47251 100644 --- a/cpp/tests/test_binary_serializer.cpp +++ b/cpp/tests/test_binary_serializer.cpp @@ -192,7 +192,7 @@ TEST(BinarySerializer, type_check) { Foo foo; auto stream = mio::serialize_binary(foo, mio::IOF_IncludeTypeInfo); - auto r = mio::deserialize_binary(stream, mio::Tag{}, mio::IOF_IncludeTypeInfo); + auto r = mio::deserialize_binary(stream, mio::Tag{}, mio::IOF_IncludeTypeInfo); EXPECT_THAT(r, IsSuccess()); } @@ -200,7 +200,7 @@ TEST(BinarySerializer, fail_type_check) { Foo foo; auto stream = mio::serialize_binary(foo, mio::IOF_IncludeTypeInfo); - auto r = mio::deserialize_binary(stream, mio::Tag{}, mio::IOF_IncludeTypeInfo); + auto r = mio::deserialize_binary(stream, mio::Tag{}, mio::IOF_IncludeTypeInfo); EXPECT_THAT(r, IsFailure(mio::StatusCode::InvalidType)); } @@ -233,4 +233,4 @@ TEST(ByteStream, reset) auto p = (unsigned char*)&i; EXPECT_TRUE(stream.read(p, sizeof(i))); EXPECT_EQ(i, 0); -} +} \ No newline at end of file diff --git a/cpp/tests/test_ide_secir.cpp b/cpp/tests/test_ide_secir.cpp index 05602523b0..ea03fea553 100755 --- a/cpp/tests/test_ide_secir.cpp +++ b/cpp/tests/test_ide_secir.cpp @@ -285,7 +285,7 @@ TEST(IdeSecir, checkSimulationFunctions) } } -// Check if the model uses the correct method for initialization using the function get_initialization_method(). +// Check if the model uses the correct method for initialization using the function get_initialization_method_compartments(). TEST(IdeSecir, checkInitializations) { using Vec = mio::TimeSeries::Vector; @@ -311,14 +311,14 @@ TEST(IdeSecir, checkInitializations) mio::isecir::Model model1(std::move(init_copy1), N, deaths, 1000); // Check that the initialization method is not already set. - EXPECT_EQ(0, model1.get_initialization_method()); + EXPECT_EQ(0, model1.get_initialization_method_compartments()); // Carry out simulation. mio::isecir::Simulation sim1(model1, dt); sim1.advance(tmax); // Verify that the expected initialization method was used. - EXPECT_EQ(1, sim1.get_model().get_initialization_method()); + EXPECT_EQ(1, sim1.get_model().get_initialization_method_compartments()); // --- Case with S. /* !! For the other tests, the contact rate is set to 0 so that the force of infection is zero. @@ -339,7 +339,7 @@ TEST(IdeSecir, checkInitializations) sim2.advance(tmax); // Verify that the expected initialization method was used. - EXPECT_EQ(2, sim2.get_model().get_initialization_method()); + EXPECT_EQ(2, sim2.get_model().get_initialization_method_compartments()); // --- Case with R. mio::TimeSeries init_copy3(init); @@ -353,7 +353,7 @@ TEST(IdeSecir, checkInitializations) sim3.advance(tmax); // Verify that the expected initialization method was used. - EXPECT_EQ(3, sim3.get_model().get_initialization_method()); + EXPECT_EQ(3, sim3.get_model().get_initialization_method_compartments()); // --- Case with forceofinfection. mio::TimeSeries init_copy4(init); @@ -364,7 +364,7 @@ TEST(IdeSecir, checkInitializations) sim4.advance(tmax); // Verify that the expected initialization method was used. - EXPECT_EQ(4, sim4.get_model().get_initialization_method()); + EXPECT_EQ(4, sim4.get_model().get_initialization_method_compartments()); // --- Case without fitting initialization method. // Deactivate temporarily log output for next tests. Errors are expected here. @@ -381,7 +381,7 @@ TEST(IdeSecir, checkInitializations) sim5.advance(tmax); // Verify that initialization was not possible with one of the models methods. - EXPECT_EQ(-1, sim5.get_model().get_initialization_method()); + EXPECT_EQ(-1, sim5.get_model().get_initialization_method_compartments()); // --- Test with negative number of deaths. deaths = -10; @@ -397,7 +397,7 @@ TEST(IdeSecir, checkInitializations) sim6.advance(tmax); // Verify that initialization was possible but the result is not appropriate. - EXPECT_EQ(-2, sim6.get_model().get_initialization_method()); + EXPECT_EQ(-2, sim6.get_model().get_initialization_method_compartments()); // Reactive log output. mio::set_log_level(mio::LogLevel::warn);