Skip to content
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.

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

946 implement generalized linear chain trick #1058

Merged
merged 34 commits into from
Nov 26, 2024
Merged
Show file tree
Hide file tree
Changes from 33 commits
Commits
Show all changes
34 commits
Select commit Hold shift + click to select a range
913de3a
Not working: first ideas to make derived of CompartmentalModel
lenaploetzke May 10, 2024
2606932
implemented glct model
lenaploetzke May 17, 2024
6d8f910
Merge branch 'main' into 946-implement-generalized-linear-chain-trick
lenaploetzke May 17, 2024
427157a
some adaptions to the glct model
lenaploetzke May 29, 2024
0bf1529
add constraints to check_constraints functions
lenaploetzke Jun 4, 2024
f38904f
correct signs
lenaploetzke Jun 4, 2024
8e66b49
Merge branch 'main' into 946-implement-generalized-linear-chain-trick
lenaploetzke Jun 21, 2024
8d2011f
Merge branch 'main' into 914-make-the-lct-secir-model-a-derived-class…
lenaploetzke Jun 24, 2024
be02fef
add ScalarType as template argument
lenaploetzke Jun 24, 2024
f9a3447
make LCT model a derived class of CompartmentalModel
lenaploetzke Jun 25, 2024
deaaeed
added test; changed int to size_t in LctInfectionState
lenaploetzke Jun 26, 2024
dd6bb88
fix bug
lenaploetzke Jun 26, 2024
99a6867
rly fix bug
lenaploetzke Jun 26, 2024
1edf777
add explicit conversion
lenaploetzke Jun 26, 2024
ff06d05
Merge branch '914-make-the-lct-secir-model-a-derived-class-of-flowmod…
lenaploetzke Jun 27, 2024
ff78759
fix error
lenaploetzke Jun 27, 2024
3614f55
make dimensions in parameters of size_t
lenaploetzke Jun 27, 2024
3e4b68a
added test
lenaploetzke Jun 28, 2024
7470dd1
added README
lenaploetzke Jun 28, 2024
34d9e6e
update README
lenaploetzke Jul 1, 2024
f8c9f89
Merge branch 'main' into 946-implement-generalized-linear-chain-trick
lenaploetzke Jul 2, 2024
6c06904
cleanup
lenaploetzke Jul 2, 2024
6e65611
Merge branch 'main' into 946-implement-generalized-linear-chain-trick
lenaploetzke Jul 22, 2024
a6674fc
[ci skip] adapt comments
lenaploetzke Jul 22, 2024
411e1f7
Merge branch 'main' into 946-implement-generalized-linear-chain-trick
lenaploetzke Aug 7, 2024
a2c3938
Merge branch 'main' into 946-implement-generalized-linear-chain-trick
lenaploetzke Oct 7, 2024
17ff131
use calculate_compartments function of ict_infection_state for LCT an…
lenaploetzke Oct 7, 2024
60da5ce
Apply suggestions from code review
lenaploetzke Oct 14, 2024
74dbc56
[ci skip] reviewer suggestions
lenaploetzke Oct 15, 2024
5262748
Include tikz picture in README
lenaploetzke Oct 18, 2024
cfcb46a
mainly documentation
lenaploetzke Oct 21, 2024
c09e34e
use version that is rendered correctly in github
lenaploetzke Oct 21, 2024
5828c41
Update README.md
lenaploetzke Oct 21, 2024
6a612dc
Typos
lenaploetzke Nov 26, 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
1 change: 1 addition & 0 deletions cpp/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,7 @@ if(MEMILIO_BUILD_MODELS)
add_subdirectory(models/ode_secir)
add_subdirectory(models/ode_secirvvs)
add_subdirectory(models/lct_secir)
add_subdirectory(models/glct_secir)
add_subdirectory(models/ide_secir)
add_subdirectory(models/ide_seir)
add_subdirectory(models/ode_seir)
Expand Down
4 changes: 4 additions & 0 deletions cpp/examples/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,10 @@ add_executable(lct_secir_example lct_secir.cpp)
target_link_libraries(lct_secir_example PRIVATE memilio lct_secir)
target_compile_options(lct_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

add_executable(glct_secir_example glct_secir.cpp)
target_link_libraries(glct_secir_example PRIVATE memilio glct_secir)
target_compile_options(glct_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

add_executable(history_example history.cpp)
target_link_libraries(history_example PRIVATE memilio)
target_compile_options(history_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
Expand Down
211 changes: 211 additions & 0 deletions cpp/examples/glct_secir.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,211 @@
/*
* Copyright (C) 2020-2024 MEmilio
*
* Authors: Lena Ploetzke
*
* Contact: Martin J. Kuehn <[email protected]>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

#include "glct_secir/model.h"
#include "glct_secir/parameters.h"
#include "memilio/config.h"
#include "memilio/utils/time_series.h"
#include "memilio/epidemiology/uncertain_matrix.h"
#include "memilio/math/eigen.h"
#include "memilio/utils/logging.h"
#include "memilio/compartments/simulation.h"
#include "memilio/data/analyze_result.h"

#include <vector>

int main()
{
// Simple example to demonstrate how to run a simulation using an GLCT SECIR model.
// This example recreates the lct_secir.cpp example using the GLCT model.
// This means, that we use the corresponding initial numbers, parameters and numbers of subcompartments that are
// equivalent to the choices in the LCT example.
// Parameters, initial values and the number of subcompartments are not meant to represent a realistic scenario.
// We need to double the choices of the number of subcompartments for some compartments,
// as we define different strains for different transition possibilities. For both strains, the same number of
// subcompartments is chosen. The transition probabilities are defined in the StartingProbabilities.
constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 6, NumInfectedSymptoms = 2, NumInfectedSevere = 2,
NumInfectedCritical = 10;
using Model = mio::glsecir::Model<NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms, NumInfectedSevere,
NumInfectedCritical>;
using LctState = Model::LctState;
using InfectionState = LctState::InfectionState;
lenaploetzke marked this conversation as resolved.
Show resolved Hide resolved

Model model;

const ScalarType tmax = 10;
const ScalarType t0 = 0;
const ScalarType dt_init = 10; ///< Initially used step size for adaptive method.

// Define some epidemiological parameters needed throughput the model definition and initialization.
const ScalarType TimeExposed = 3.2;
const ScalarType TimeInfectedNoSymptoms = 2.;
const ScalarType TimeInfectedSymptoms = 5.8;
const ScalarType TimeInfectedSevere = 9.5;
const ScalarType TimeInfectedCritical = 7.1;
const ScalarType RecoveredPerInfectedNoSymptoms = 0.09;
const ScalarType SeverePerInfectedSymptoms = 0.2;
const ScalarType CriticalPerSevere = 0.25;
const ScalarType DeathsPerCritical = 0.3;
lenaploetzke marked this conversation as resolved.
Show resolved Hide resolved

// Define the initial values with the distribution of the population into subcompartments.
// This method of defining the initial values using a vector of vectors is not necessary, 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 values directly.
// Split the initial population defined in the lct_secir.cpp example according to the transition probabilities in
// the two strains per compartment.
std::vector<std::vector<ScalarType>> initial_populations = {
{750},
{30, 20},
{20 * (1 - RecoveredPerInfectedNoSymptoms), 10 * (1 - RecoveredPerInfectedNoSymptoms),
10 * (1 - RecoveredPerInfectedNoSymptoms), 20 * RecoveredPerInfectedNoSymptoms,
10 * RecoveredPerInfectedNoSymptoms, 10 * RecoveredPerInfectedNoSymptoms},
{50 * SeverePerInfectedSymptoms, 50 * (1 - SeverePerInfectedSymptoms)},
{50 * CriticalPerSevere, 50 * (1 - CriticalPerSevere)},
{10 * DeathsPerCritical, 10 * DeathsPerCritical, 5 * DeathsPerCritical, 3 * DeathsPerCritical,
2 * DeathsPerCritical, 10 * (1 - DeathsPerCritical), 10 * (1 - DeathsPerCritical), 5 * (1 - DeathsPerCritical),
3 * (1 - DeathsPerCritical), 2 * (1 - DeathsPerCritical)},
{20},
{10}};
annawendler marked this conversation as resolved.
Show resolved Hide resolved

// Assert that initial_populations has the right shape.
if (initial_populations.size() != (size_t)InfectionState::Count) {
mio::log_error("The number of vectors in initial_populations does not match the number of InfectionStates.");
return 1;
}
if ((initial_populations[(size_t)InfectionState::Susceptible].size() !=
LctState::get_num_subcompartments<InfectionState::Susceptible>()) ||
(initial_populations[(size_t)InfectionState::Exposed].size() != NumExposed) ||
(initial_populations[(size_t)InfectionState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) ||
(initial_populations[(size_t)InfectionState::InfectedSymptoms].size() != NumInfectedSymptoms) ||
(initial_populations[(size_t)InfectionState::InfectedSevere].size() != NumInfectedSevere) ||
(initial_populations[(size_t)InfectionState::InfectedCritical].size() != NumInfectedCritical) ||
(initial_populations[(size_t)InfectionState::Recovered].size() !=
LctState::get_num_subcompartments<InfectionState::Recovered>()) ||
(initial_populations[(size_t)InfectionState::Dead].size() !=
LctState::get_num_subcompartments<InfectionState::Dead>())) {
mio::log_error("The length of at least one vector in initial_populations does not match the related number of "
"subcompartments.");
annawendler marked this conversation as resolved.
Show resolved Hide resolved
return 1;
}

// Transfer the initial values in initial_populations to the model.
std::vector<ScalarType> flat_initial_populations;
for (auto&& vec : initial_populations) {
flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end());
}
for (size_t i = 0; i < LctState::Count; i++) {
model.populations[mio::Index<LctState>(i)] = flat_initial_populations[i];
}

// Set Parameters according to the LCT model definitions, e.g. use Erlang-distributions.
// Exposed.
// The get_default of the StartingProbabilities returns the first unit vector of the defined size.
// It is necessary to set it although the default method is used to define the length of the vector.
model.parameters.get<mio::glsecir::StartingProbabilitiesExposed>() =
mio::glsecir::StartingProbabilitiesExposed().get_default(
LctState::get_num_subcompartments<InfectionState::Exposed>());
// The get_default function returns the TransitionMatrix that is required to have an Erlang-distributed
// stay time with an average of TimeExposed.
model.parameters.get<mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms>() =
mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default(
LctState::get_num_subcompartments<InfectionState::Exposed>(), TimeExposed);
// This definition of the StartingProbability and the TransitionMatrix lead to an Erlang-distributed latent stage.

// InfectedNoSymptoms.
// For InfectedNoSymptoms, two strains has to be defined, one for the Transition
// InfectedNoSymptomsToInfectedSymptoms and one for InfectedNoSymptomsToRecovered.
// The strains have a length of NumInfectedNoSymptoms/2. each.
// The transition probability is included in the StartingProbability vector.
mio::Vector<ScalarType> StartingProbabilitiesInfectedNoSymptoms =
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>());
StartingProbabilitiesInfectedNoSymptoms[0] = 1 - RecoveredPerInfectedNoSymptoms;
StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)(
LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.)] = RecoveredPerInfectedNoSymptoms;
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedNoSymptoms>() =
StartingProbabilitiesInfectedNoSymptoms;
// Define equal TransitionMatrices for the strains.
// They follow the same Erlang-distribution such that we get the same result as with one strain in the LCT model.
model.parameters.get<mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms>() =
mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.),
TimeInfectedNoSymptoms);
model.parameters.get<mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered>() =
mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.),
TimeInfectedNoSymptoms);
// Do the same for all compartments.
// InfectedSymptoms.
mio::Vector<ScalarType> StartingProbabilitiesInfectedSymptoms =
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>());
StartingProbabilitiesInfectedSymptoms[0] = SeverePerInfectedSymptoms;
StartingProbabilitiesInfectedSymptoms[(Eigen::Index)(
LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.)] = 1 - SeverePerInfectedSymptoms;
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedSymptoms>() = StartingProbabilitiesInfectedSymptoms;
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere>() =
mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.), TimeInfectedSymptoms);
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered>() =
mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.), TimeInfectedSymptoms);
// InfectedSevere.
mio::Vector<ScalarType> StartingProbabilitiesInfectedSevere =
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedSevere>());
StartingProbabilitiesInfectedSevere[0] = CriticalPerSevere;
StartingProbabilitiesInfectedSevere[(Eigen::Index)(
LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.)] = 1 - CriticalPerSevere;
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedSevere>() = StartingProbabilitiesInfectedSevere;
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical>() =
mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.), TimeInfectedSevere);
model.parameters.get<mio::glsecir::TransitionMatrixInfectedSevereToRecovered>() =
mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.), TimeInfectedSevere);
// InfectedCritical.
mio::Vector<ScalarType> StartingProbabilitiesInfectedCritical =
mio::Vector<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedCritical>());
StartingProbabilitiesInfectedCritical[0] = DeathsPerCritical;
StartingProbabilitiesInfectedCritical[(Eigen::Index)(
LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.)] = 1 - DeathsPerCritical;
model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedCritical>() = StartingProbabilitiesInfectedCritical;
model.parameters.get<mio::glsecir::TransitionMatrixInfectedCriticalToDead>() =
mio::glsecir::TransitionMatrixInfectedCriticalToDead().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.), TimeInfectedCritical);
model.parameters.get<mio::glsecir::TransitionMatrixInfectedCriticalToRecovered>() =
mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default(
(size_t)(LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.), TimeInfectedCritical);

model.parameters.get<mio::glsecir::TransmissionProbabilityOnContact>() = 0.05;

mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::glsecir::ContactPatterns>();
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10));
// From SimulationTime 5, the contact pattern is reduced to 30% of the initial value.
contact_matrix[0].add_damping(0.7, mio::SimulationTime(5.));

model.parameters.get<mio::glsecir::RelativeTransmissionNoSymptoms>() = 0.7;
model.parameters.get<mio::glsecir::RiskOfInfectionFromSymptomatic>() = 0.25;

// Perform a simulation.
mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(t0, tmax, dt_init, model);
// The simulation result is divided by subcompartments as in the LCT model.
// We call the function calculate_compartments to get a result according to the InfectionStates.
mio::TimeSeries<ScalarType> population_no_subcompartments = model.calculate_compartments(result);
auto interpolated_result = mio::interpolate_simulation_result(population_no_subcompartments, 0.1);
interpolated_result.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 10, 4);
}
9 changes: 4 additions & 5 deletions cpp/examples/lct_secir.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,9 @@ int main()

// 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;
bool use_initializer_flows = false;

ScalarType tmax = 20;
ScalarType tmax = 10;

// Set Parameters.
model.parameters.get<mio::lsecir::TimeExposed>()[0] = 3.2;
Expand Down Expand Up @@ -155,10 +155,9 @@ int main()

// Perform a simulation.
mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.5, model);
//result.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
// The simulation result is divided by subcompartments.
// We call the function calculate_comparttments to get a result according to the InfectionStates.
// We call the function calculate_compartments to get a result according to the InfectionStates.
mio::TimeSeries<ScalarType> population_no_subcompartments = model.calculate_compartments(result);
auto interpolated_results = mio::interpolate_simulation_result(population_no_subcompartments);
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 12, 4);
}
50 changes: 50 additions & 0 deletions cpp/memilio/epidemiology/lct_infection_state.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,9 @@
#ifndef MIO_EPI_LCT_INFECTION_STATE_H
#define MIO_EPI_LCT_INFECTION_STATE_H

#include "memilio/config.h"
#include "memilio/math/eigen.h"

#include <array>

namespace mio
Expand Down Expand Up @@ -78,6 +81,53 @@ class LctInfectionState
return index;
}

/**
* @brief Cumulates a vector with the number of individuals in each subcompartment (with subcompartments
* according to the LctInfectionState) to produce a Vector that divides the population only into the infection
* states defined in InfectionStates.
*
* @param[in] subcompartments Vector with number of individuals in each subcompartment.
* The size of the vector has to match the LctInfectionState.
* @return Vector with accumulated values for the InfectionStates.
*/
static Vector<ScalarType> calculate_compartments(const Vector<ScalarType>& subcompartments)
{
assert(subcompartments.rows() == Count);
annawendler marked this conversation as resolved.
Show resolved Hide resolved

Vector<ScalarType> compartments((Eigen::Index)InfectionState::Count);
// Use segment of the vector subcompartments of each InfectionState and sum up the values of subcompartments.
compartments[(Eigen::Index)InfectionState::Susceptible] = subcompartments[0];
compartments[(Eigen::Index)InfectionState::Exposed] =
subcompartments
.segment(get_first_index<InfectionState::Exposed>(), get_num_subcompartments<InfectionState::Exposed>())
.sum();
compartments[(Eigen::Index)InfectionState::InfectedNoSymptoms] =
subcompartments
.segment(get_first_index<InfectionState::InfectedNoSymptoms>(),
get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
.sum();
compartments[(Eigen::Index)InfectionState::InfectedSymptoms] =
subcompartments
.segment(get_first_index<InfectionState::InfectedSymptoms>(),
get_num_subcompartments<InfectionState::InfectedSymptoms>())
.sum();
compartments[(Eigen::Index)InfectionState::InfectedSevere] =
subcompartments
.segment(get_first_index<InfectionState::InfectedSevere>(),
get_num_subcompartments<InfectionState::InfectedSevere>())
.sum();
compartments[(Eigen::Index)InfectionState::InfectedCritical] =
subcompartments
.segment(get_first_index<InfectionState::InfectedCritical>(),
get_num_subcompartments<InfectionState::InfectedCritical>())
.sum();
compartments[(Eigen::Index)InfectionState::Recovered] =
subcompartments[get_first_index<InfectionState::Recovered>()];
compartments[(Eigen::Index)InfectionState::Dead] = subcompartments[get_first_index<InfectionState::Dead>()];

return compartments;
}

static constexpr size_t Count{(... + Ns)};

private:
Expand Down
12 changes: 12 additions & 0 deletions cpp/models/glct_secir/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
add_library(glct_secir
infection_state.h
model.h
model.cpp
parameters.h
)
target_link_libraries(glct_secir PUBLIC memilio)
target_include_directories(glct_secir PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
)
target_compile_options(glct_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
Loading
Loading