Skip to content

897 add initialization method from flows for lct model #931

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

Closed
lenaploetzke wants to merge 22 commits into from
Closed

897 add initialization method from flows for lct model #931

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
43663a2
added lct model
lenaploetzke Feb 2, 2024
fde4506
Small update of README
lenaploetzke Feb 6, 2024
87333d7
rewrite description of example
lenaploetzke Feb 6, 2024
d0444e4
rewrite comment
lenaploetzke Feb 6, 2024
c7b10c7
included suggestions from code review
lenaploetzke Feb 6, 2024
83d9ef4
Apply suggestions from code review
lenaploetzke Feb 6, 2024
84153b3
resolved error
lenaploetzke Feb 6, 2024
aa1db0b
added parameters loc and scale
lenaploetzke Feb 7, 2024
298e628
added other SAF and extended boost archive
lenaploetzke Feb 7, 2024
09eddbf
Merge branch '892-implement-simple-lct-model' into 897-add-initializa…
lenaploetzke Feb 7, 2024
7ffa12d
added initializer from flows
lenaploetzke Feb 7, 2024
ba94faa
Merge branch 'main' into 897-add-initialization-method-from-flows-for…
lenaploetzke Feb 16, 2024
afb993d
added method for the tolerance
lenaploetzke Feb 16, 2024
b2b9691
added explanation to readme
lenaploetzke Feb 16, 2024
3073154
added Initializer to example
lenaploetzke Feb 16, 2024
ddef273
Merge branch 'main' into 897-add-initialization-method-from-flows-for…
lenaploetzke Feb 29, 2024
81c7fa2
Merge branch 'main' into 897-add-initialization-method-from-flows-for…
lenaploetzke Mar 7, 2024
27d5658
part of merge missed
lenaploetzke Mar 7, 2024
fee0583
add template to class Initializer
lenaploetzke Mar 8, 2024
5d000b0
fix error
lenaploetzke Mar 11, 2024
35509b8
maybe fix error now
lenaploetzke Mar 11, 2024
4cd50ce
Merge branch 'main' into 897-add-initialization-method-from-flows-for…
lenaploetzke Mar 22, 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
180 changes: 115 additions & 65 deletions cpp/examples/lct_secir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include "lct_secir/model.h"
#include "lct_secir/infection_state.h"
#include "lct_secir/simulation.h"
#include "lct_secir/initializer_flows.h"
#include "memilio/config.h"
#include "memilio/utils/time_series.h"
#include "memilio/epidemiology/uncertain_matrix.h"
Expand All @@ -37,73 +38,14 @@ int main()
using Model = mio::lsecir::Model<2, 3, 1, 1, 5>;
using LctState = Model::LctState;

ScalarType tmax = 20;
// Variable defines whether the class Initializer is used to define an initial vector from flows or whether a manually
// defined initial vector is used to initialize the LCT model.
bool use_initializer_flows = true;

// Define the initial value vector init with the distribution of the population into subcompartments.
// This method of defining the vector using a vector of vectors is a bit of overhead, but should remind you how
// the entries of the initial value vector relate to the defined template parameters of the model or the number of subcompartments.
// It is also possible to define the initial value vector directly.
std::vector<std::vector<ScalarType>> initial_populations = {{750}, {30, 20}, {20, 10, 10}, {50},
{50}, {10, 10, 5, 3, 2}, {20}, {10}};

// Assert that initial_populations has the right shape.
if (initial_populations.size() != (int)LctState::InfectionState::Count) {
mio::log_error("The number of vectors in initial_populations does not match the number of InfectionStates.");
return 1;
}
if ((initial_populations[(int)LctState::InfectionState::Susceptible].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Susceptible>()) ||
(initial_populations[(int)LctState::InfectionState::Exposed].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Exposed>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedNoSymptoms>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedSymptoms].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedSymptoms>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedSevere].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedSevere>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedCritical].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedCritical>()) ||
(initial_populations[(int)LctState::InfectionState::Recovered].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Recovered>()) ||
(initial_populations[(int)LctState::InfectionState::Dead].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Dead>())) {
mio::log_error("The length of at least one vector in initial_populations does not match the related number of "
"subcompartments.");
return 1;
}

// Transfer the initial values in initial_populations to the vector init.
Eigen::VectorXd init = Eigen::VectorXd::Zero(LctState::Count);
init[(int)LctState::get_first_index<LctState::InfectionState::Susceptible>()] =
initial_populations[(int)LctState::InfectionState::Susceptible][0];
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::Exposed>(); i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::Exposed>() + i] =
initial_populations[(int)LctState::InfectionState::Exposed][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedNoSymptoms>();
i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedNoSymptoms>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedSymptoms>(); i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedSymptoms>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedSymptoms][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedSevere>(); i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedSevere>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedSevere][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedCritical>(); i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedCritical>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedCritical][i];
}
init[(int)LctState::get_first_index<LctState::InfectionState::Recovered>()] =
initial_populations[(int)LctState::InfectionState::Recovered][0];
init[(int)LctState::get_first_index<LctState::InfectionState::Dead>()] =
initial_populations[(int)LctState::InfectionState::Dead][0];
ScalarType tmax = 20;

// Initialize model.
Model model(std::move(init));
// Initialize a model.
Model model(std::move(Eigen::VectorXd::Zero(LctState::Count)));

// Set Parameters.
model.parameters.get<mio::lsecir::TimeExposed>() = 3.2;
Expand All @@ -127,6 +69,114 @@ int main()
model.parameters.get<mio::lsecir::CriticalPerSevere>() = 0.25;
model.parameters.set<mio::lsecir::DeathsPerCritical>(0.3);

if (use_initializer_flows) {
// Example how to use the class Initializer for the definition of an initial vector for the LCT model.
// Create TimeSeries with num_transitions elements.
int num_transitions = (int)mio::lsecir::InfectionTransition::Count;
mio::TimeSeries<ScalarType> flows(num_transitions);

ScalarType dt = 0.1;

mio::TimeSeries<ScalarType>::Vector vec_flows(num_transitions);
vec_flows[(int)mio::lsecir::InfectionTransition::SusceptibleToExposed] = 2.0;
vec_flows[(int)mio::lsecir::InfectionTransition::ExposedToInfectedNoSymptoms] = 1.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] = 8.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedNoSymptomsToRecovered] = 4.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedSymptomsToInfectedSevere] = 1.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedSymptomsToRecovered] = 4.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedSevereToInfectedCritical] = 1.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedSevereToRecovered] = 1.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedCriticalToDead] = 1.0;
vec_flows[(int)mio::lsecir::InfectionTransition::InfectedCriticalToRecovered] = 1.0;
// Add initial time point to time series.
flows.add_time_point(-110, vec_flows);
// Add further time points until time 0.
while (flows.get_last_time() < -1e-10) {
flows.add_time_point(flows.get_last_time() + dt, vec_flows);
}

// Set initialization vector for the LCT model.
mio::lsecir::Initializer<Model> initializer(std::move(flows), model);
initializer.set_tol_for_support_max(1e-6);
auto status = initializer.compute_initialization_vector(1000000, 10, 16000);
if (status) {
return 1;
}
}
else {
// Simple example how to initialize model.
// Define the initial value vector init with the distribution of the population into subcompartments.
// This method of defining the vector using a vector of vectors is a bit of overhead, but should remind you how
// the entries of the initial value vector relate to the defined template parameters of the model or the number of subcompartments.
// It is also possible to define the initial value vector directly.
std::vector<std::vector<ScalarType>> initial_populations = {{750}, {30, 20}, {20, 10, 10}, {50},
{50}, {10, 10, 5, 3, 2}, {20}, {10}};

// Assert that initial_populations has the right shape.
if (initial_populations.size() != (int)LctState::InfectionState::Count) {
mio::log_error(
"The number of vectors in initial_populations does not match the number of InfectionStates.");
return 1;
}
if ((initial_populations[(int)LctState::InfectionState::Susceptible].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Susceptible>()) ||
(initial_populations[(int)LctState::InfectionState::Exposed].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Exposed>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedNoSymptoms>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedSymptoms].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedSymptoms>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedSevere].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedSevere>()) ||
(initial_populations[(int)LctState::InfectionState::InfectedCritical].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::InfectedCritical>()) ||
(initial_populations[(int)LctState::InfectionState::Recovered].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Recovered>()) ||
(initial_populations[(int)LctState::InfectionState::Dead].size() !=
LctState::get_num_subcompartments<LctState::InfectionState::Dead>())) {
mio::log_error(
"The length of at least one vector in initial_populations does not match the related number of "
"subcompartments.");
return 1;
}

// Define initial value vector with subcompartments and transfer the initial values in initial_populations to the vector.
Eigen::VectorXd init(LctState::Count);
init[(int)LctState::get_first_index<LctState::InfectionState::Susceptible>()] =
initial_populations[(int)LctState::InfectionState::Susceptible][0];
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::Exposed>(); i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::Exposed>() + i] =
initial_populations[(int)LctState::InfectionState::Exposed][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedNoSymptoms>();
i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedNoSymptoms>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedSymptoms>();
i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedSymptoms>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedSymptoms][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedSevere>();
i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedSevere>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedSevere][i];
}
for (unsigned int i = 0; i < LctState::get_num_subcompartments<LctState::InfectionState::InfectedCritical>();
i++) {
init[(int)LctState::get_first_index<LctState::InfectionState::InfectedCritical>() + i] =
initial_populations[(int)LctState::InfectionState::InfectedCritical][i];
}
init[(int)LctState::get_first_index<LctState::InfectionState::Recovered>()] =
initial_populations[(int)LctState::InfectionState::Recovered][0];
init[(int)LctState::get_first_index<LctState::InfectionState::Dead>()] =
initial_populations[(int)LctState::InfectionState::Dead][0];

// Set initial values for the model.
model.set_initial_values(std::move(init));
}

// Perform a simulation.
mio::TimeSeries<ScalarType> result = mio::lsecir::simulate(0, tmax, 0.5, model);
// Calculate the distribution in the InfectionState%s without subcompartments of the result and print it.
Expand Down
Loading