add_executable(lct_secir_2_diseases_example lct_secir_2_diseases.cpp)

target_link_libraries(lct_secir_2_diseases_example PRIVATE memilio lct_secir_2_diseases)

target_compile_options(lct_secir_2_diseases_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

#include "lct_secir_2_diseases/model.h"

#include "lct_secir_2_diseases/infection_state.h"

#include "memilio/config.h"

#include "memilio/utils/time_series.h"

#include "memilio/epidemiology/lct_infection_state.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 LCT-SECIR-2-DISEASE model.

// One single AgeGroup/Category member is used here.

// Parameters, initial values and the number of subcompartments are not meant to represent a realistic scenario.

// The number of subcompartments can be chosen for most of the compartments:

constexpr size_t NumExposed_1a = 2, NumInfectedNoSymptoms_1a = 3, NumInfectedSymptoms_1a = 3,

NumInfectedSevere_1a = 3, NumInfectedCritical_1a = 2, NumExposed_2a = 1,

NumInfectedNoSymptoms_2a = 2, NumInfectedSymptoms_2a = 2, NumInfectedSevere_2a = 2,

NumInfectedCritical_2a = 1, NumExposed_1b = 2, NumInfectedNoSymptoms_1b = 3,

NumInfectedSymptoms_1b = 3, NumInfectedSevere_1b = 3, NumInfectedCritical_1b = 2,

NumExposed_2b = 1, NumInfectedNoSymptoms_2b = 2, NumInfectedSymptoms_2b = 2,

NumInfectedSevere_2b = 2, NumInfectedCritical_2b = 1;

// The compartment for Susceptible people and all compartments for Dead and Recovered people must have exactly one subcompartment:

constexpr size_t NumSusceptible = 1, NumDead_a = 1, NumDead_b = 1, NumRecovered_1a = 1, NumRecovered_1b = 1,

NumRecovered_2ab = 1;

using InfState = mio::lsecir2d::InfectionState;

using LctState =

mio::LctInfectionState<ScalarType, InfState, NumSusceptible, NumExposed_1a, NumInfectedNoSymptoms_1a,

NumInfectedSymptoms_1a, NumInfectedSevere_1a, NumInfectedCritical_1a, NumRecovered_1a,

NumDead_a, NumExposed_2a, NumInfectedNoSymptoms_2a, NumInfectedSymptoms_2a,

NumInfectedSevere_2a, NumInfectedCritical_2a, NumExposed_1b, NumInfectedNoSymptoms_1b,

NumInfectedSymptoms_1b, NumInfectedSevere_1b, NumInfectedCritical_1b, NumRecovered_1b,

NumDead_b, NumExposed_2b, NumInfectedNoSymptoms_2b, NumInfectedSymptoms_2b,

NumInfectedSevere_2b, NumInfectedCritical_2b, NumRecovered_2ab>;

using Model = mio::lsecir2d::Model<ScalarType, LctState>;

Model model;

ScalarType tmax = 5;

// Set Parameters.

model.parameters.get<mio::lsecir2d::TimeExposed_a<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeExposed_b<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedNoSymptoms_a<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedNoSymptoms_b<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedSymptoms_a<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedSymptoms_b<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedSevere_a<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedSevere_b<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedCritical_a<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TimeInfectedCritical_b<ScalarType>>()[0] = 3.;

model.parameters.get<mio::lsecir2d::TransmissionProbabilityOnContact_a<ScalarType>>()[0] = 0.1;

model.parameters.get<mio::lsecir2d::TransmissionProbabilityOnContact_b<ScalarType>>()[0] = 0.1;

mio::ContactMatrixGroup<ScalarType>& contact_matrix =

model.parameters.get<mio::lsecir2d::ContactPatterns<ScalarType>>();

contact_matrix[0] = mio::ContactMatrix<ScalarType>(Eigen::MatrixX<ScalarType>::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<ScalarType>(5.));

model.parameters.get<mio::lsecir2d::RelativeTransmissionNoSymptoms_a<ScalarType>>()[0] = 0.7;

model.parameters.get<mio::lsecir2d::RelativeTransmissionNoSymptoms_b<ScalarType>>()[0] = 0.7;

model.parameters.get<mio::lsecir2d::RiskOfInfectionFromSymptomatic_a<ScalarType>>()[0] = 0.25;

model.parameters.get<mio::lsecir2d::RiskOfInfectionFromSymptomatic_b<ScalarType>>()[0] = 0.25;

model.parameters.get<mio::lsecir2d::RecoveredPerInfectedNoSymptoms_a<ScalarType>>()[0] = 0.09;

model.parameters.get<mio::lsecir2d::RecoveredPerInfectedNoSymptoms_b<ScalarType>>()[0] = 0.09;

model.parameters.get<mio::lsecir2d::SeverePerInfectedSymptoms_a<ScalarType>>()[0] = 0.2;

model.parameters.get<mio::lsecir2d::SeverePerInfectedSymptoms_b<ScalarType>>()[0] = 0.2;

model.parameters.get<mio::lsecir2d::CriticalPerSevere_a<ScalarType>>()[0] = 0.25;

model.parameters.get<mio::lsecir2d::CriticalPerSevere_b<ScalarType>>()[0] = 0.25;

model.parameters.get<mio::lsecir2d::DeathsPerCritical_a<ScalarType>>()[0] = 0.3;

model.parameters.get<mio::lsecir2d::DeathsPerCritical_b<ScalarType>>()[0] = 0.3;

// Simple example how to initialize model without flows.

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

std::vector<std::vector<ScalarType>> initial_populations = {

{200}, {0, 0}, {30, 10, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0}, {0}, {0}, {0},

{0, 0}, {10, 0}, {0, 0}, {0}, {10, 0}, {30, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0},

{0}, {0}, {100}, {0, 0}, {0, 0}, {0, 0}, {0}, {0}};

// Assert that initial_populations has the right shape.

if (initial_populations.size() != (size_t)InfState::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)InfState::Susceptible].size() !=

LctState::get_num_subcompartments<InfState::Susceptible>()) ||

(initial_populations[(size_t)InfState::Exposed_1a].size() != NumExposed_1a) ||

(initial_populations[(size_t)InfState::InfectedNoSymptoms_1a].size() != NumInfectedNoSymptoms_1a) ||

(initial_populations[(size_t)InfState::InfectedSymptoms_1a].size() != NumInfectedSymptoms_1a) ||

(initial_populations[(size_t)InfState::InfectedSevere_1a].size() != NumInfectedSevere_1a) ||

(initial_populations[(size_t)InfState::InfectedCritical_1a].size() != NumInfectedCritical_1a) ||

(initial_populations[(size_t)InfState::Exposed_2a].size() != NumExposed_2a) ||

(initial_populations[(size_t)InfState::InfectedNoSymptoms_2a].size() != NumInfectedNoSymptoms_2a) ||

(initial_populations[(size_t)InfState::InfectedSymptoms_2a].size() != NumInfectedSymptoms_2a) ||

(initial_populations[(size_t)InfState::InfectedSevere_2a].size() != NumInfectedSevere_2a) ||

(initial_populations[(size_t)InfState::InfectedCritical_2a].size() != NumInfectedCritical_2a) ||

(initial_populations[(size_t)InfState::Recovered_1a].size() !=

LctState::get_num_subcompartments<InfState::Recovered_1a>()) ||

(initial_populations[(size_t)InfState::Dead_a].size() !=

LctState::get_num_subcompartments<InfState::Dead_a>()) ||

(initial_populations[(size_t)InfState::Exposed_1b].size() != NumExposed_1b) ||

(initial_populations[(size_t)InfState::InfectedNoSymptoms_1b].size() != NumInfectedNoSymptoms_1b) ||

(initial_populations[(size_t)InfState::InfectedSymptoms_1b].size() != NumInfectedSymptoms_1b) ||

(initial_populations[(size_t)InfState::InfectedSevere_1b].size() != NumInfectedSevere_1b) ||

(initial_populations[(size_t)InfState::InfectedCritical_1b].size() != NumInfectedCritical_1b) ||

(initial_populations[(size_t)InfState::Recovered_1b].size() !=

LctState::get_num_subcompartments<InfState::Recovered_1b>()) ||

(initial_populations[(size_t)InfState::Dead_b].size() !=

LctState::get_num_subcompartments<InfState::Dead_b>()) ||

(initial_populations[(size_t)InfState::Exposed_2b].size() != NumExposed_2b) ||

(initial_populations[(size_t)InfState::InfectedNoSymptoms_2b].size() != NumInfectedNoSymptoms_2b) ||

(initial_populations[(size_t)InfState::InfectedSymptoms_2b].size() != NumInfectedSymptoms_2b) ||

(initial_populations[(size_t)InfState::InfectedSevere_2b].size() != NumInfectedSevere_2b) ||

(initial_populations[(size_t)InfState::InfectedCritical_2b].size() != NumInfectedCritical_2b) ||

(initial_populations[(size_t)InfState::Recovered_2ab].size() !=

LctState::get_num_subcompartments<InfState::Recovered_2ab>())) {

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 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[i] = flat_initial_populations[i];

}

// Perform a simulation.

mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.1, model);

// The simulation result is divided by subcompartments.

// 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", " E1a", " C1a", " I1a", " H1a", " U1a", " Ra",

" Da", " E2a", " C2a", " I2a", " H2a", " U2a", " E1b",

" C1b", " I1b", " H1b", " U1b", " Rb", " Db", " E2b",

" C2b", " I2b", " H2b", " U2b", " Rab"},

6, 2);

}

for (int i = 0; i < (Eigen::Index)InfectionState::Count; i++) {

InfectionState State = static_cast<InfectionState>(i);

// first index of first subcompartment:

size_t index = 0;

for (size_t j = 0; j < (size_t)(State); j++) {

index = index + m_subcompartment_numbers[j];

}

// number of subcompartments:

size_t num_subcomp = m_subcompartment_numbers[(size_t)State];

compartments[i] = subcompartments.segment(index, num_subcomp).sum();

}

This directory contains a model implementation based on an ODE formulation using the Linear Chain Trick.

add_library(lct_secir_2_diseases

infection_state.h

model.h

model.cpp

parameters.h

)

target_link_libraries(lct_secir_2_diseases PUBLIC memilio)

target_include_directories(lct_secir_2_diseases PUBLIC

$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>

$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>

)

target_compile_options(lct_secir_2_diseases PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

This directory contains a model implementation for two different diseases or two variants of a disease based on an ODE formulation using the Linear Chain Trick.

To get started, check out the [official documentation](https://memilio.readthedocs.io/en/latest/cpp/models/lsecir2d.html)

or the [LCT2D minimal example](../../examples/lct_secir_2_diseases.cpp).

#ifndef LCT_SECIR_2_DISEASES_INFECTIONSTATE_H

#define LCT_SECIR_2_DISEASES_INFECTIONSTATE_H

namespace mio

{

namespace lsecir2d

{

enum class InfectionState

{

Susceptible = 0,

// Notation: State_[Infection number][disease]

// first infection with disease a

Exposed_1a = 1,

InfectedNoSymptoms_1a = 2,

InfectedSymptoms_1a = 3,

InfectedSevere_1a = 4,

InfectedCritical_1a = 5,

Recovered_1a = 6,

Dead_a = 7,

// second infection with disease a

Exposed_2a = 8,

InfectedNoSymptoms_2a = 9,

InfectedSymptoms_2a = 10,

InfectedSevere_2a = 11,

InfectedCritical_2a = 12,

// first infection with disease b

Exposed_1b = 13,

InfectedNoSymptoms_1b = 14,

InfectedSymptoms_1b = 15,

InfectedSevere_1b = 16,

InfectedCritical_1b = 17,

Recovered_1b = 18,

Dead_b = 19,

// second infection with disease b

Exposed_2b = 20,

InfectedNoSymptoms_2b = 21,

InfectedSymptoms_2b = 22,

InfectedSevere_2b = 23,

InfectedCritical_2b = 24,

// Recovered from both diseases

Recovered_2ab = 25,

Count = 26

};

} // namespace lsecir2d

} // namespace mio

#endif // LCT_SECIR_2_DISEASES_INFECTIONSTATE_H

#include "lct_secir_2_diseases/model.h"

namespace mio

{

namespace lsecir2d

{

} // namespace lsecir2d

} // namespace mio