Skip to content

769 new initialization for ide model #770

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
annawendler merged 9 commits into from
Feb 16, 2024
Merged

769 new initialization for ide model #770

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
d736a4d
added different initialization scheme
lenaploetzke Sep 15, 2023
e5a421a
added test for initialization and use deaths at 0
lenaploetzke Sep 16, 2023
4cd48a6
make initialization consistent
lenaploetzke Feb 1, 2024
52be71c
typo
lenaploetzke Feb 1, 2024
3c88ec4
Merge branch 'main' into 769-new-initialization-for-ide-model
lenaploetzke Feb 1, 2024
800f258
added review suggestions
lenaploetzke Feb 2, 2024
a1f25d6
added a comment
lenaploetzke Feb 13, 2024
1e6802d
Merge branch 'main' into 769-new-initialization-for-ide-model
lenaploetzke Feb 13, 2024
704f83a
changed test procedure
lenaploetzke Feb 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions cpp/examples/ide_secir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,10 @@ int main()
{
using Vec = mio::TimeSeries<ScalarType>::Vector;

ScalarType tmax = 10;
ScalarType N = 10000;
ScalarType Dead_before = 12;
ScalarType dt = 1;
ScalarType tmax = 10;
ScalarType N = 10000;
ScalarType deaths = 13.10462213;
ScalarType dt = 1;

int num_transitions = (int)mio::isecir::InfectionTransition::Count;

Expand Down Expand Up @@ -63,7 +63,7 @@ int main()
}

// Initialize model.
mio::isecir::Model model(std::move(init), N, Dead_before);
mio::isecir::Model model(std::move(init), N, deaths);

// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000;
// model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0;
Expand Down
124 changes: 76 additions & 48 deletions cpp/models/ide_secir/model.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,64 +29,33 @@ namespace mio
namespace isecir
{

Model::Model(TimeSeries<ScalarType>&& init, ScalarType N_init, ScalarType Dead_before,
Model::Model(TimeSeries<ScalarType>&& init, ScalarType N_init, ScalarType deaths, ScalarType total_confirmed_cases,
const ParameterSet& Parameterset_init)
: parameters{Parameterset_init}
, m_transitions{std::move(init)}
, m_populations{TimeSeries<ScalarType>(Eigen::Index(InfectionState::Count))}
, m_N{N_init}
, m_deaths_before{Dead_before}
, m_total_confirmed_cases{total_confirmed_cases}
{
m_deaths_before =
deaths - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)];
m_populations.add_time_point<Eigen::VectorXd>(
0, TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count, 0));
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] =
m_deaths_before + m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)];
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] = deaths;
}

void Model::initialize(ScalarType dt)
{
// compute Susceptibles at time 0 and m_forceofinfection at time -m_dt as initial values for discretization scheme
// use m_forceofinfection at -m_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);
if (m_forceofinfection > 0) {
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] /
(dt * m_forceofinfection);

//calculate other compartment sizes for t=0
if (m_total_confirmed_cases > 1e-12) {
m_initialization_method = 1;
other_compartments_current_timestep(dt);

//R; need an initial value for R, therefore do not calculate via compute_recovered()
// The scheme of the ODE model for initialization is applied here.
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)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
m_total_confirmed_cases - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
}
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
//calculate other compartment sizes for t=0
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)] =
m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
}
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, calculate Susceptibles via other compartments
//determining other compartment sizes is not dependent of Susceptibles(0), just of the transitions of the past.
//calculate other compartment sizes for t=0
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)] -
Expand All @@ -98,10 +67,69 @@ void Model::initialize(ScalarType dt)
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
}
else {
log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither "
"Susceptibles nor Recovered for time 0 were set. One of them should be larger 0.");
}

// compute Susceptibles at time 0 and m_forceofinfection at time -dt as initial values for discretization scheme
// use m_forceofinfection at -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);
if (m_forceofinfection > 1e-12) {
m_initialization_method = 2;
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] /
(dt * m_forceofinfection);

//calculate other compartment sizes for t=0
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)] =
m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
}
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 = 3;
//calculate other compartment sizes for t=0
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)] =
m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
}
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, calculate Susceptibles via other compartments
//determining other compartment sizes is not dependent of Susceptibles(0), just of the transitions of the past.
//calculate other compartment sizes for t=0
m_initialization_method = 4;
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)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
}
else {
m_initialization_method = -1;
log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither "
"Susceptibles nor Recovered or total confirmed cases for time 0 were set. One of them should be "
"larger 0.");
}
}
// compute m_forceofinfection at time 0 needed for further simulation
update_forceofinfection(dt);
}
Expand All @@ -118,11 +146,11 @@ void Model::compute_susceptibles(ScalarType dt)
void Model::compute_flow(int idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt)
{
ScalarType sum = 0;
/* In order to satisfy TransitionDistribution(m_dt*i) = 0 for all i >= k, k is determined by the maximum support of the distribution.
/* In order to satisfy TransitionDistribution(dt*i) = 0 for all i >= k, k is determined by the maximum support of the distribution.
Since we are using a backwards difference scheme to compute the derivative, we have that the
derivative of TransitionDistribution(m_dt*i) = 0 for all i >= k+1.
derivative of TransitionDistribution(dt*i) = 0 for all i >= k+1.

Hence calc_time_index goes until std::ceil(support_max/m_dt) since for std::ceil(support_max/m_dt)+1 all terms are already zero.
Hence calc_time_index goes until std::ceil(support_max/dt) since for std::ceil(support_max/dt)+1 all terms are already zero.
This needs to be adjusted if we are changing the finite difference scheme */

Eigen::Index calc_time_index = (Eigen::Index)std::ceil(
Expand All @@ -131,7 +159,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)* m_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;

Expand Down Expand Up @@ -214,7 +242,7 @@ void Model::update_forceofinfection(ScalarType dt, bool initialization)
ScalarType deaths;

if (initialization) {
// determine m_forceofinfection at time -m_dt which is the penultimate timepoint in m_transitions
// determine m_forceofinfection at time -dt which is the penultimate timepoint in m_transitions
num_time_points = m_transitions.get_num_time_points() - 1;
current_time = -dt;
deaths = m_deaths_before;
Expand Down
30 changes: 24 additions & 6 deletions cpp/models/ide_secir/model.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,17 +39,17 @@ class 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_init from one compartment to another.
* which transit within one timestep dt from one compartment to another.
* Possible transitions are specified in as #InfectionTransition%s.
* Considered points of times should have the distance dt_init and the last time point should be 0.
* Considered points of times should have the distance dt and the last time point should be 0.
* 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] dt_init The size of the time step used for numerical simulation.
* @param[in] N_init The population of the considered region.
* @param[in] Dead_before The total number of deaths at the time point - dt_init.
* @param[in] deaths The total number of deaths at the time zero.
* @param[in] total_confirmed_cases Total confirmed cases at time t0 can be set if it should be used for initialisation.
* @param[in, out] Parameterset_init Used Parameters for simulation.
*/
Model(TimeSeries<ScalarType>&& init, ScalarType N_init, ScalarType Dead_before,
Model(TimeSeries<ScalarType>&& init, ScalarType N_init, ScalarType deaths, ScalarType total_confirmed_cases = 0,
const ParameterSet& Parameterset_init = ParameterSet());

/**
Expand Down Expand Up @@ -205,6 +205,21 @@ class Model
m_tol = new_tol;
}

/**
* @brief Specifies a number associated with the method used for initialization.
*
* @returns 0 if the initialization method has not yet been selected,
* 1 if the method using the total number of confirmed cases at time 0 is used,
* 2 if the force of infection method is used,
* 3 if the initialization is calculated using a prior set value for S,
* 4 if the initialization is calculated using a prior set value for R and
* -1 if the initialization was not possible using any of the methods.
*/
int get_initialization_method()
{
return m_initialization_method;
}

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. */
TimeSeries<ScalarType>
Expand All @@ -215,8 +230,11 @@ class Model
private:
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}; ///< Deaths before start of simulation (at time -m_dt).
ScalarType m_deaths_before{0}; ///< Total number of deaths at the time point - 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.
};

} // namespace isecir
Expand Down
Loading