// Add first time point in m_populations according to last time point in m_transitions which is where we start

// the simulation.

calc_time = std::max(m_transitiondistributions_support_max[idx_TransitionDistribution1],

m_transitiondistributions_support_max[idx_TransitionDistribution2]);

calc_time = m_transitiondistributions_support_max[idx_TransitionDistribution1];

// It is possible to calculate the sizes of the other compartments in advance because only the initial values of

// the flows are used.

// 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.

// Recovered; calculated as the difference between the total population and the sum of the other

// compartment sizes.

// Another check would be if the sum of the compartments is equal to N, but in all initialization methods one of

// the compartments is initialized via N - the others.

// This also means that if a compartment is greater than N, we will always have one or more compartments less than

// zero.

/* In order to satisfy TransitionDistribution(dt*i) = 0 for all i >= k, k is determined by the maximum support of

Eigen::Index calc_time_index =

(Eigen::Index)std::ceil(m_transitiondistributions_support_max[idx_InfectionTransitions] / dt);

// (current_time_index - i) is the index corresponding to time the individuals have already spent in this state.

sum += m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index - i] *

m_transitions[i + 1][idx_IncomingFlow];

// Calculate flow SusceptibleToExposed with force of infection from previous time step and Susceptibles from

// current time step.

ScalarType calc_time =

std::max({m_transitiondistributions_support_max[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms],

m_transitiondistributions_support_max[(int)InfectionTransition::InfectedNoSymptomsToRecovered],

m_transitiondistributions_support_max[(int)InfectionTransition::InfectedSymptomsToInfectedSevere],

m_transitiondistributions_support_max[(int)InfectionTransition::InfectedSymptomsToRecovered]});

ScalarType state_age = state_age_index * dt;

(m_transitiondistributions_in_forceofinfection[0][num_time_points - i - 1] *

m_transitiondistributions_in_forceofinfection[1][num_time_points - i - 1] *

void Model::set_transitiondistributions_support_max(ScalarType dt)

{

// We do not consider the transition SusceptibleToExposed as it is not needed in the computations.

for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) {

m_transitiondistributions_support_max[transition] =

parameters.get<TransitionDistributions>()[transition].get_support_max(dt, m_tol);

}

}

void Model::set_transitiondistributions_derivative(ScalarType dt)

{

// We do not consider the transition SusceptibleToExposed as it is not needed in the computations.

for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) {

Eigen::Index support_max_index =

(Eigen::Index)std::ceil(m_transitiondistributions_support_max[transition] / dt);

// Create vec_derivative that contains the value of the approximated derivative for all necessary time points.

// Here, we evaluate the derivative at time points t_0, ..., t_{support_max_index}.

std::vector<ScalarType> vec_derivative(support_max_index + 1, 0.);

for (Eigen::Index i = 0; i <= support_max_index; i++) {

// Compute state_age for considered index.

ScalarType state_age = (ScalarType)i * dt;

// Compute derivative.

vec_derivative[i] = (parameters.get<TransitionDistributions>()[transition].eval(state_age) -

parameters.get<TransitionDistributions>()[transition].eval(state_age - dt)) /

dt;

}

m_transitiondistributions_derivative[transition] = vec_derivative;

}

}

void Model::set_transitiondistributions_in_forceofinfection(ScalarType dt)

{

// Relevant transitions for force of infection term.

std::vector<std::vector<int>> relevant_transitions = {

{(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,

(int)InfectionTransition::InfectedNoSymptomsToRecovered},

{(int)InfectionTransition::InfectedSymptomsToInfectedSevere,

(int)InfectionTransition::InfectedSymptomsToRecovered}};

// Determine the relevant calculation area = union of the supports of the relevant transition distributions.

ScalarType calc_time = std::max({m_transitiondistributions_support_max[relevant_transitions[0][0]],

m_transitiondistributions_support_max[relevant_transitions[0][1]],

m_transitiondistributions_support_max[relevant_transitions[1][0]],

m_transitiondistributions_support_max[relevant_transitions[1][1]]});

// Corresponding index.

// Need to evaluate survival functions at t_0, ..., t_{calc_time_index} for computation of force of infection,

// 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;

// Compute contributions from survival functions and transition probabilities starting in InfectedNoSymptoms and

// InfectedSymptoms, respectively, on force of infection term.

for (int contribution = 0; contribution < 2; contribution++) {

std::vector<ScalarType> vec_contribution_to_foi(calc_time_index + 1, 0.);

for (Eigen::Index i = 0; i <= calc_time_index; i++) {

// Compute state_age for all indices from t_0, ..., t_{calc_time_index}.

ScalarType state_age = i * dt;

vec_contribution_to_foi[i] =

parameters.get<TransitionProbabilities>()[relevant_transitions[contribution][0]] *

parameters.get<TransitionDistributions>()[relevant_transitions[contribution][0]].eval(state_age) +

parameters.get<TransitionProbabilities>()[relevant_transitions[contribution][1]] *

parameters.get<TransitionDistributions>()[relevant_transitions[contribution][1]].eval(state_age);

}

m_transitiondistributions_in_forceofinfection[contribution] = vec_contribution_to_foi;

}

}

// Note that this function computes the global_support_max via the get_support_max() method and does not make use

// of the vector m_transitiondistributions_support_max. This is because the global_support_max is already used in

// check_constraints and we cannot ensure that the vector has already been computed when checking for constraints

// (which usually happens before setting the initial flows and simulating).