diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 45a3d67773..4855a17a68 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -156,6 +156,8 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/ode_secirvvs) add_subdirectory(models/lct_secir) add_subdirectory(models/glct_secir) + add_subdirectory(models/ide_endemic_secir) + add_subdirectory(models/ide_endemic_sird) add_subdirectory(models/ide_secir) add_subdirectory(models/ide_seir) add_subdirectory(models/ode_seir) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index e7c5cb99f4..a255394e3a 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -104,6 +104,14 @@ add_executable(twitter_mobility_example twitter_mobility.cpp) target_link_libraries(twitter_mobility_example PRIVATE memilio ode_secir) target_compile_options(twitter_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ide_endemic_secir_example ide_endemic_secir.cpp) +target_link_libraries(ide_endemic_secir_example PRIVATE memilio ide_endemic_secir) +target_compile_options(ide_endemic_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(ide_endemic_sird_example ide_endemic_sird.cpp) +target_link_libraries(ide_endemic_sird_example PRIVATE memilio ide_endemic_sird) +target_compile_options(ide_endemic_sird_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(ide_seir_example ide_seir.cpp) target_link_libraries(ide_seir_example PRIVATE memilio ide_seir) target_compile_options(ide_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) @@ -175,8 +183,19 @@ if(MEMILIO_HAS_HDF5) target_compile_options(ode_secir_save_results_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) endif() + if(MEMILIO_HAS_JSONCPP) add_executable(ide_initialization_example ide_initialization.cpp) target_link_libraries(ide_initialization_example PRIVATE memilio ide_secir) target_compile_options(ide_initialization_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) endif() + +if(MEMILIO_HAS_HDF5) + add_executable(ide_endemic_secir_save_results_example ide_endemic_secir_save_results.cpp) + target_link_libraries(ide_endemic_secir_save_results_example PRIVATE memilio ide_endemic_secir) + target_compile_options(ide_endemic_secir_save_results_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + + add_executable(ide_endemic_sird_save_results_example ide_endemic_sird_save_results.cpp) + target_link_libraries(ide_endemic_sird_save_results_example PRIVATE memilio ide_endemic_sird) + target_compile_options(ide_endemic_sird_save_results_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +endif() \ No newline at end of file diff --git a/cpp/examples/ide_endemic_secir.cpp b/cpp/examples/ide_endemic_secir.cpp new file mode 100644 index 0000000000..13113dc3cc --- /dev/null +++ b/cpp/examples/ide_endemic_secir.cpp @@ -0,0 +1,192 @@ +#include "ide_endemic_secir/computed_parameters.h" +#include "ide_endemic_secir/model.h" +#include "ide_endemic_secir/infection_state.h" +#include "ide_endemic_secir/normalized_model.h" +#include "ide_endemic_secir/parameters.h" +#include "ide_endemic_secir/simulation.h" +#include "memilio/config.h" +#include "memilio/math/eigen.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/time_series.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/data/analyze_result.h" +#include +#include +#include + +int main() +{ + using Vec = mio::TimeSeries::Vector; + + ScalarType tmax = 1; + ScalarType dt = 0.00001; + + int num_states = static_cast(mio::endisecir::InfectionState::Count); + int num_transitions = static_cast(mio::endisecir::InfectionTransition::Count); + + // Create TimeSeries with num_states elements where states needed for simulation will be stored. + mio::TimeSeries init(num_states); + + Vec vec_init(num_states); + + vec_init[static_cast(mio::endisecir::InfectionState::Susceptible)] = 100000.; + vec_init[static_cast(mio::endisecir::InfectionState::Exposed)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedNoSymptoms)] = 10.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedSymptoms)] = 20.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedSevere)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedCritical)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::Recovered)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::Dead)] = 0.; + + init.add_time_point(0, vec_init); + + mio::endisecir::CompParameters computed_parameters(std::move(init)); + + //Set working parameters + + // mio::ExponentialSurvivalFunction exp(3.0); + // mio::StateAgeFunctionWrapper delaydistribution(exp); + // std::vector vec_delaydistribution(num_transitions, delaydistribution); + + mio::SmootherCosine smoothcos(2.0); + mio::StateAgeFunctionWrapper delaydistribution(smoothcos); + std::vector vec_delaydistribution(num_transitions, delaydistribution); + + // Uncomment for Lognorm. + // mio::ConstantFunction initialfunc(0); + // mio::StateAgeFunctionWrapper delaydistributioninit(initialfunc); + // std::vector vec_delaydistribution(num_transitions, delaydistributioninit); + // // ExposedToInfectedNoSymptoms + // mio::LognormSurvivalFunction survivalExposedToInfectedNoSymptoms(0.3, 0, 4.2); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::ExposedToInfectedNoSymptoms].set_state_age_function( + // survivalExposedToInfectedNoSymptoms); + // // InfectedNoSymptomsToInfectedSymptoms + // mio::LognormSurvivalFunction survivalInfectedNoSymptomsToInfectedSymptoms(0.7, 0, 0.8); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] + // .set_state_age_function(survivalInfectedNoSymptomsToInfectedSymptoms); + // // InfectedNoSymptomsToRecovered + // mio::LognormSurvivalFunction survivalInfectedNoSymptomsToRecovered(0.2, 0, 7.7); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedNoSymptomsToRecovered] + // .set_state_age_function(survivalInfectedNoSymptomsToRecovered); + // // InfectedSymptomsToInfectedSevere + // mio::LognormSurvivalFunction survivalInfectedSymptomsToInfectedSevere(0.7, 0, 5.3); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSymptomsToInfectedSevere] + // .set_state_age_function(survivalInfectedSymptomsToInfectedSevere); + // // InfectedSymptomsToRecovered + // mio::LognormSurvivalFunction survivalInfectedSymptomsToRecovered(0.2, 0, 7.8); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSymptomsToRecovered].set_state_age_function( + // survivalInfectedSymptomsToRecovered); + // // InfectedSevereToInfectedCritical + // mio::LognormSurvivalFunction survivalInfectedSevereToInfectedCritical(1.0, 0, 0.9); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSevereToInfectedCritical] + // .set_state_age_function(survivalInfectedSevereToInfectedCritical); + // // InfectedSevereToRecovered + // mio::LognormSurvivalFunction survivalInfectedSevereToRecovered(0.3, 0, 17.1); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSevereToRecovered].set_state_age_function( + // survivalInfectedSevereToRecovered); + // // InfectedCriticalToDead + // mio::LognormSurvivalFunction survivalInfectedCriticalToDead(0.4, 0, 9.8); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedCriticalToDead].set_state_age_function( + // survivalInfectedCriticalToDead); + // // InfectedCriticalToRecovered + // mio::LognormSurvivalFunction survivalInfectedCriticalToRecovered(0.3, 0, 17.1); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedCriticalToRecovered].set_state_age_function( + // survivalInfectedCriticalToRecovered); + + computed_parameters.parameters.get() = vec_delaydistribution; + + std::vector vec_prob((int)mio::endisecir::InfectionTransition::Count, 1.); + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] = 0.8; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedNoSymptomsToRecovered)] = 1 - 0.8; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSymptomsToInfectedSevere)] = 0.1; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSymptomsToRecovered)] = 1 - 0.1; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSevereToInfectedCritical)] = 0.2; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSevereToRecovered)] = 1 - 0.2; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedCriticalToDead)] = 0.4; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedCriticalToRecovered)] = 1 - 0.4; + + computed_parameters.parameters.set(vec_prob); + + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.)); + computed_parameters.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + mio::ConstantFunction constant(0.1); + mio::StateAgeFunctionWrapper constant_prob(constant); + + computed_parameters.parameters.get() = constant_prob; + + mio::ExponentialSurvivalFunction exponential(0.5); + mio::StateAgeFunctionWrapper exponential_prob(exponential); + + computed_parameters.parameters.get() = exponential_prob; + computed_parameters.parameters.get() = exponential_prob; + + computed_parameters.parameters.set(4e-1); + computed_parameters.parameters.set(3e-1); + + computed_parameters.set_tol_for_support_max(1e-6); + + mio::endisecir::Model model(computed_parameters); + + mio::endisecir::NormModel normmodel(computed_parameters); + + // start the simulation. + mio::endisecir::Simulation sim(computed_parameters, model, normmodel, dt); + sim.advance(tmax); + + //Get the compartments of model and print them. + auto interpolated_results = mio::interpolate_simulation_result(sim.get_compartments(), dt / 2.); + interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8); + // Uncomment to print the compartments computed with the update scheme. + // auto interpolated_results_update = mio::interpolate_simulation_result(sim.get_compartments_update(), dt / 2.); + // interpolated_results_update.print_table({"US", "UE", "UC", "UI", "UH", "UU", "UR", "UD"}, 16, 8); + + //Get the commpartments of normmodel and print them. + // auto interpolated_normresults = mio::interpolate_simulation_result(sim.get_normmodel_compartments(), dt / 2.); + // interpolated_normresults.print_table({"s", "e", "c", "i", "h", "u", "r"}, 16, 8); + + // Uncomment to print the normalized compartments of model. + // sim.get_normalizedcompartments().print_table({"S/N", "E/N", "C/N", "I/N", "H/N", "U/N", "R/N"}, 16, 8); + + // Uncomment to print the transitions of model. + // sim.get_transitions().print_table({"S->E", "E->C", "C->I", "C->R", "I->H", "I->R", "H->U", "H->R", "U->D", "U->R"}, + // 16, 8); + // sim.get_transitions_update().print_table( + // {"US->UE", "UE->UC", "UC->UI", "UC->UR", "UI->UH", "UI->UR", "UH->UU", "UH->UR", "UU->UD", "UU->UR"}, 16, 8); + + // Uncomment to print the transitions of normmodel. + // sim.get_normmodel_transitions().print_table( + // {"s->e", "e->c", "c->i", "c->r", "i->h", "i->r , "h->u", "h->r", "u->d", "u->r"}, 16, 8); + + // Uncomment to print the force of infection of model. + auto interpolated_FoI = mio::interpolate_simulation_result(sim.get_forceofinfections(), dt / 2.); + interpolated_FoI.print_table({"FoI"}, 16, 8); + // sim.get_forceofinfections_update().print_table({"FoIUpdate"}, 16, 8); + + // Uncomment to print the force of infection of normmodel. + // sim.get_normmodel_forceofinfections().print_table({"norm FoI"}, 16, 8); + + //Uncomment to print the reproduction number + std::cout << "The reproduction number Rc = " << sim.get_reproductionnumber_c() << "\n"; + + // Uncomment to print the the values T_z1^z2 + //for (int i = 0; i < (int)sim.get_T().size(); i++) { + // std::cout << "T_" << i << " = " << sim.get_T()[i] << "\n"; + //} + + // Uncomment to print the the values V^z + //for (int i = 0; i < (int)sim.get_V().size(); i++) { + // std::cout << "V_" << i << " = " << sim.get_V()[i] << "\n"; + //} + + // Uncomment to print the values W_z + //for (int i = 0; i < (int)sim.get_W().size(); i++) { + // std::cout << "W_" << i << " = " << sim.get_W()[i] << "\n"; + //} + + // Uncomment to print the total population size. + // sim.get_totalpopulations().print_table({"N"}, 16, 9); + // sim.get_totalpopulations_update().print_table({"UN"}, 16, 9); +} \ No newline at end of file diff --git a/cpp/examples/ide_endemic_secir_save_results.cpp b/cpp/examples/ide_endemic_secir_save_results.cpp new file mode 100644 index 0000000000..9ccc42d053 --- /dev/null +++ b/cpp/examples/ide_endemic_secir_save_results.cpp @@ -0,0 +1,229 @@ +#include "memilio/io/result_io.h" +#include "memilio/utils/time_series.h" +#include "memilio/config.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/math/floating_point.h" + +#include "ide_endemic_secir/infection_state.h" +#include "ide_endemic_secir/model.h" +#include "ide_endemic_secir/simulation.h" +#include "ide_endemic_secir/parameters.h" + +#include "boost/filesystem.hpp" +#include +#include +#include + +mio::IOResult simulate_endidemodel(ScalarType tmax, std::string save_dir = "") +{ + using Vec = mio::TimeSeries::Vector; + + ScalarType dt = 1.0; + + int num_states = static_cast(mio::endisecir::InfectionState::Count); + int num_transitions = static_cast(mio::endisecir::InfectionTransition::Count); + + // Create TimeSeries with num_states elements where states needed for simulation will be stored. + mio::TimeSeries init(num_states); + + Vec vec_init(num_states); + + vec_init[static_cast(mio::endisecir::InfectionState::Susceptible)] = 100000.; + vec_init[static_cast(mio::endisecir::InfectionState::Exposed)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedNoSymptoms)] = 10.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedSymptoms)] = 20.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedSevere)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::InfectedCritical)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::Recovered)] = 0.; + vec_init[static_cast(mio::endisecir::InfectionState::Dead)] = 0.; + + init.add_time_point(0, vec_init); + + mio::endisecir::CompParameters computed_parameters(std::move(init)); + + //Set working parameters + + // mio::ExponentialSurvivalFunction exp(3.0); + // mio::StateAgeFunctionWrapper delaydistribution(exp); + // std::vector vec_delaydistribution(num_transitions, delaydistribution); + + mio::SmootherCosine smoothcos(6.0); + mio::StateAgeFunctionWrapper delaydistribution(smoothcos); + std::vector vec_delaydistribution(num_transitions, delaydistribution); + + //Uncomment for Lognorm. + // mio::ConstantFunction initialfunc(0); + // mio::StateAgeFunctionWrapper delaydistributioninit(initialfunc); + // std::vector vec_delaydistribution(num_transitions, delaydistributioninit); + // // ExposedToInfectedNoSymptoms + // mio::LognormSurvivalFunction survivalExposedToInfectedNoSymptoms(0.3, 0, 4.2); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::ExposedToInfectedNoSymptoms].set_state_age_function( + // survivalExposedToInfectedNoSymptoms); + // // InfectedNoSymptomsToInfectedSymptoms + // mio::LognormSurvivalFunction survivalInfectedNoSymptomsToInfectedSymptoms(0.7, 0, 0.8); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] + // .set_state_age_function(survivalInfectedNoSymptomsToInfectedSymptoms); + // // InfectedNoSymptomsToRecovered + // mio::LognormSurvivalFunction survivalInfectedNoSymptomsToRecovered(0.2, 0, 7.7); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedNoSymptomsToRecovered] + // .set_state_age_function(survivalInfectedNoSymptomsToRecovered); + // // InfectedSymptomsToInfectedSevere + // mio::LognormSurvivalFunction survivalInfectedSymptomsToInfectedSevere(0.7, 0, 5.3); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSymptomsToInfectedSevere] + // .set_state_age_function(survivalInfectedSymptomsToInfectedSevere); + // // InfectedSymptomsToRecovered + // mio::LognormSurvivalFunction survivalInfectedSymptomsToRecovered(0.2, 0, 7.8); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSymptomsToRecovered].set_state_age_function( + // survivalInfectedSymptomsToRecovered); + // // InfectedSevereToInfectedCritical + // mio::LognormSurvivalFunction survivalInfectedSevereToInfectedCritical(1.0, 0, 0.9); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSevereToInfectedCritical] + // .set_state_age_function(survivalInfectedSevereToInfectedCritical); + // // InfectedSevereToRecovered + // mio::LognormSurvivalFunction survivalInfectedSevereToRecovered(0.3, 0, 17.1); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedSevereToRecovered].set_state_age_function( + // survivalInfectedSevereToRecovered); + // // InfectedCriticalToDead + // mio::LognormSurvivalFunction survivalInfectedCriticalToDead(0.4, 0, 9.8); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedCriticalToDead].set_state_age_function( + // survivalInfectedCriticalToDead); + // // InfectedCriticalToRecovered + // mio::LognormSurvivalFunction survivalInfectedCriticalToRecovered(0.3, 0, 17.1); + // vec_delaydistribution[(int)mio::endisecir::InfectionTransition::InfectedCriticalToRecovered].set_state_age_function( + // survivalInfectedCriticalToRecovered); + + computed_parameters.parameters.get() = vec_delaydistribution; + + std::vector vec_prob((int)mio::endisecir::InfectionTransition::Count, 1.); + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] = 0.8; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedNoSymptomsToRecovered)] = 1 - 0.8; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSymptomsToInfectedSevere)] = 0.1; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSymptomsToRecovered)] = 1 - 0.1; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSevereToInfectedCritical)] = 0.2; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSevereToRecovered)] = 1 - 0.2; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedCriticalToDead)] = 0.4; + vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedCriticalToRecovered)] = 1 - 0.4; + + //Different values: + // std::vector vec_prob((int)mio::endisecir::InfectionTransition::Count, 1.); + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] = 0.5; + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedNoSymptomsToRecovered)] = 1 - 0.5; + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSymptomsToInfectedSevere)] = 0.5; + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSymptomsToRecovered)] = 1 - 0.5; + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSevereToInfectedCritical)] = 0.5; + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedSevereToRecovered)] = 1 - 0.5; + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedCriticalToDead)] = 0.8; + // vec_prob[Eigen::Index(mio::endisecir::InfectionTransition::InfectedCriticalToRecovered)] = 1 - 0.8; + + computed_parameters.parameters.set(vec_prob); + + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.)); + computed_parameters.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + mio::ConstantFunction constant(0.1); + mio::StateAgeFunctionWrapper constant_prob(constant); + + computed_parameters.parameters.get() = constant_prob; + + mio::ExponentialSurvivalFunction exponential(0.5); + mio::StateAgeFunctionWrapper exponential_prob(exponential); + + computed_parameters.parameters.get() = exponential_prob; + computed_parameters.parameters.get() = exponential_prob; + + computed_parameters.parameters.set(9e-2); + computed_parameters.parameters.set(7e-3); + + //computed_parameters.set_tol_for_support_max(1e-6); + + mio::endisecir::Model model(computed_parameters); + + mio::endisecir::NormModel normmodel(computed_parameters); + + // start the simulation. + mio::endisecir::Simulation sim(computed_parameters, model, normmodel, dt); + sim.advance(tmax); + + if (!save_dir.empty()) { + std::string tmax_string = std::to_string(tmax); + std::string dt_string = std::to_string(dt); + + std::string filename_ide = save_dir + "analysis_endide_SECIR_" + tmax_string.substr(0, tmax_string.find(".")) + + "_" + dt_string.substr(0, dt_string.find(".") + 5); + + //Save files of Model. + //Save total population. + std::string filename_ide_totalpopulation1 = filename_ide + "_totalpopulation1.h5"; + mio::IOResult save_result_status_tp1 = + mio::save_result({sim.get_totalpopulations()}, {0}, 1, filename_ide_totalpopulation1); + //Save derivative of total population. + std::string filename_ide_dertotalpopulation1 = filename_ide + "_dertotalpopulation1.h5"; + mio::IOResult save_result_status_dtp1 = + mio::save_result({sim.get_totalpopulations_derivative()}, {0}, 1, filename_ide_dertotalpopulation1); + + //Save compartments. + std::string filename_ide_compartments1 = filename_ide + "_compartments1.h5"; + mio::IOResult save_result_status_c1 = + mio::save_result({sim.get_compartments()}, {0}, 1, filename_ide_compartments1); + + //Save normalized compartments. + std::string filename_ide_normcompartments = filename_ide + "_normcompartments.h5"; + mio::IOResult save_result_status_nc = + mio::save_result({sim.get_normalizedcompartments()}, {0}, 1, filename_ide_normcompartments); + + //Save the force of infection. + std::string filename_ide_forceofinfection = filename_ide + "_forceofinfection.h5"; + mio::IOResult save_result_status_foi = + mio::save_result({sim.get_forceofinfections()}, {0}, 1, filename_ide_forceofinfection); + + //Save files of NormModel. + //Save compartments. + std::string filename_ide_normmod_compartments = filename_ide + "_normmod_compartments1.h5"; + mio::IOResult save_result_status_nmc1 = + mio::save_result({sim.get_normmodel_compartments()}, {0}, 1, filename_ide_normmod_compartments); + + //Save the force of infection. + std::string filename_ide_normmod_forceofinfection = filename_ide + "_normmod_forceofinfection1.h5"; + mio::IOResult save_result_status_nmfoi = + mio::save_result({sim.get_normmodel_forceofinfections()}, {0}, 1, filename_ide_normmod_forceofinfection); + + //Safe difference between normalized compartments. + std::string filename_ide_difference_normcomp = filename_ide + "_difference_normcomp.h5"; + mio::IOResult save_result_status_dnm = + mio::save_result({sim.get_difference_normalizationcomp()}, {0}, 1, filename_ide_difference_normcomp); + + //Safe difference between Force of Infections. + std::string filename_ide_difference_normFoI = filename_ide + "_difference_normfoi.h5"; + mio::IOResult save_result_status_dnf = + mio::save_result({sim.get_difference_normalizationFoi()}, {0}, 1, filename_ide_difference_normFoI); + + if (!save_result_status_tp1 || !save_result_status_c1 || !save_result_status_nc || !save_result_status_foi || + !save_result_status_nmc1 || !save_result_status_nmfoi || !save_result_status_dnm || + !save_result_status_dnf || !save_result_status_dtp1) { + return mio::failure(mio::StatusCode::UnknownError, "Error while saving results."); + } + } + // Print the reproduction number. + std::cout << "The reproduction number Rc for Birth rate > Death rate is " << sim.get_reproductionnumber_c() + << "\n"; + + return mio::success(); +} + +int main() +{ + std::string result_dir = "/localdata1/trit_ha/code/memilio-1/PythonPlotsEndIDESECIR/simulation_results/"; + + // Define tmax. + ScalarType tmax = 500; + + auto result_ide = simulate_endidemodel(tmax, result_dir); + if (!result_ide) { + printf("%s\n", result_ide.error().formatted_message().c_str()); + return -1; + } + + return 0; +} \ No newline at end of file diff --git a/cpp/examples/ide_endemic_sird.cpp b/cpp/examples/ide_endemic_sird.cpp new file mode 100644 index 0000000000..6e43b9e475 --- /dev/null +++ b/cpp/examples/ide_endemic_sird.cpp @@ -0,0 +1,178 @@ +#include "ide_endemic_sird/computed_parameters.h" +#include "ide_endemic_sird/model.h" +#include "ide_endemic_sird/infection_state.h" +#include "ide_endemic_sird/normalized_model.h" +#include "ide_endemic_sird/parameters.h" +#include "ide_endemic_sird/simulation.h" +#include "memilio/config.h" +#include "memilio/math/eigen.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/time_series.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/data/analyze_result.h" +#include +#include +#include + +int main() +{ + using Vec = mio::TimeSeries::Vector; + + ScalarType tmax = 100; + ScalarType dt = 1.0; + + int num_states = static_cast(mio::endisird::InfectionState::Count); + int num_transitions = static_cast(mio::endisird::InfectionTransition::Count); + + // Create TimeSeries with num_states elements where states needed for simulation will be stored. + mio::TimeSeries init(num_states); + + Vec vec_init(num_states); + + vec_init[static_cast(mio::endisird::InfectionState::Susceptible)] = 100000.; + vec_init[static_cast(mio::endisird::InfectionState::Infected)] = 40.; + vec_init[static_cast(mio::endisird::InfectionState::Recovered)] = 0.; + vec_init[static_cast(mio::endisird::InfectionState::Dead)] = 0.; + + init.add_time_point(0, vec_init); + + mio::endisird::CompParameters computed_parameters(std::move(init)); + + //Set working parameters + + // mio::ExponentialSurvivalFunction exp(3.0); + // mio::StateAgeFunctionWrapper delaydistribution(exp); + // std::vector vec_delaydistribution(num_transitions, delaydistribution); + + mio::SmootherCosine smoothcos(8.0); + mio::StateAgeFunctionWrapper delaydistribution(smoothcos); + std::vector vec_delaydistribution(num_transitions, delaydistribution); + + // mio::ConstantFunction initialfunc(0); + // mio::StateAgeFunctionWrapper delaydistributioninit(initialfunc); + // std::vector vec_delaydistribution(num_transitions, delaydistributioninit); + // // InfectedToDead + // mio::SmootherCosine survivalExposedToInfectedNoSymptoms(6.0); + // vec_delaydistribution[(int)mio::endisird::InfectionTransition::InfectedToDead].set_state_age_function( + // survivalExposedToInfectedNoSymptoms); + // // InfectedToRecovered + // mio::SmootherCosine survivalInfectedNoSymptomsToInfectedSymptoms(8.0); + // vec_delaydistribution[(int)mio::endisird::InfectionTransition::InfectedToRecovered].set_state_age_function( + // survivalInfectedNoSymptomsToInfectedSymptoms); + + // mio::ConstantFunction initialfunc(0); + // mio::StateAgeFunctionWrapper delaydistributioninit(initialfunc); + // std::vector vec_delaydistribution(num_transitions, delaydistributioninit); + // // InfectedToDead + // mio::LognormSurvivalFunction survivalExposedToInfectedNoSymptoms(0.8, 0, 4.2); + // vec_delaydistribution[(int)mio::endisird::InfectionTransition::InfectedToDead].set_state_age_function( + // survivalExposedToInfectedNoSymptoms); + // // InfectedToRecovered + // mio::LognormSurvivalFunction survivalInfectedNoSymptomsToInfectedSymptoms(0.2, 0, 0.8); + // vec_delaydistribution[(int)mio::endisird::InfectionTransition::InfectedToRecovered].set_state_age_function( + // survivalInfectedNoSymptomsToInfectedSymptoms); + + computed_parameters.parameters.get() = vec_delaydistribution; + + std::vector vec_prob((int)mio::endisird::InfectionTransition::Count, 1.); + vec_prob[Eigen::Index(mio::endisird::InfectionTransition::InfectedToDead)] = 0.1; + vec_prob[Eigen::Index(mio::endisird::InfectionTransition::InfectedToRecovered)] = 1 - 0.1; + + computed_parameters.parameters.set(vec_prob); + + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.)); + computed_parameters.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + mio::ConstantFunction constant(0.1); + mio::StateAgeFunctionWrapper constant_prob(constant); + + computed_parameters.parameters.get() = constant_prob; + + mio::ExponentialSurvivalFunction exponential(0.5); + mio::StateAgeFunctionWrapper exponential_prob(exponential); + + computed_parameters.parameters.get() = exponential_prob; + + computed_parameters.parameters.set(4e-3); + computed_parameters.parameters.set(3e-3); + + computed_parameters.set_tol_for_support_max(1e-6); + + mio::endisird::Model model(computed_parameters); + + mio::endisird::NormModel normmodel(computed_parameters); + + // start the simulation. + mio::endisird::Simulation sim(computed_parameters, model, normmodel, dt); + sim.advance(tmax); + + //Get the compartments of model and print them. + // auto interpolated_results = mio::interpolate_simulation_result(sim.get_compartments(), dt / 2.); + // interpolated_results.print_table({"S", "I", "R", "D "}, 16, 8); + + //Get the commpartments of normmodel and print them. + auto interpolated_normresults = mio::interpolate_simulation_result(sim.get_normmodel_compartments(), dt / 2.); + interpolated_normresults.print_table({"s", "i", "r"}, 16, 8); + + // Uncomment to print the normalized compartments of model. + // sim.get_normalizedcompartments().print_table({"S/N", "I/N", "R/N"}, 16, 8); + + // Uncomment to print the transitions of model. + // sim.get_transitions().print_table({"S->I", "I->D", "I->R"}, + // 16, 8); + + // Uncomment to print the transitions of normmodel. + // sim.get_normmodel_transitions().print_table( + // {"s->i", "i->d", "i->r }, 16, 8); + + // Uncomment to print the force of infection of model. + // auto interpolated_FoI = mio::interpolate_simulation_result(sim.get_forceofinfections(), dt / 2.); + // interpolated_FoI.print_table({"FoI"}, 16, 8); + + // Uncomment to print the force of infection of normmodel. + // sim.get_normmodel_forceofinfections().print_table({"norm FoI"}, 16, 8); + + sim.get_size().print_table({"sum"}, 16, 8); + //Uncomment to print the reproduction number + // std::cout << "The reproduction number Rc = " << sim.get_reproductionnumber_c() << "\n"; + + // Uncomment to print the the values T_z1^z2 + // for (int i = 0; i < (int)sim.get_T().size(); i++) { + // std::cout << "T_" << i << " = " << sim.get_T()[i] << "\n"; + // } + + // Uncomment to print W_i + + // std::cout << "W_i = " << sim.get_W() << "\n"; + + // Uncomment to print the equilibria + // std::cout << "l_1 = " << sim.get_Equilibirum_FoI()[0] << "\n"; + // std::cout << "l_2 = " << sim.get_Equilibirum_FoI()[1] << "\n"; + // std::cout << "s_1 = " << sim.get_Equilibirum_compartments()[0][(int)mio::endisird::InfectionState::Susceptible] + // << "\n"; + // std::cout << "s_2 = " << sim.get_Equilibirum_compartments()[1][(int)mio::endisird::InfectionState::Susceptible] + // << "\n"; + // std::cout << "i_1 = " << sim.get_Equilibirum_compartments()[0][(int)mio::endisird::InfectionState::Infected] + // << "\n"; + // std::cout << "i_2 = " << sim.get_Equilibirum_compartments()[1][(int)mio::endisird::InfectionState::Infected] + // << "\n"; + // std::cout << "r_1 = " << sim.get_Equilibirum_compartments()[0][(int)mio::endisird::InfectionState::Recovered] + // << "\n"; + // std::cout << "r_2 = " << sim.get_Equilibirum_compartments()[1][(int)mio::endisird::InfectionState::Recovered] + // << "\n"; + // std::cout << "sigmaid_1 = " + // << sim.get_Equilibirum_transitions()[0][(int)mio::endisird::InfectionTransition::InfectedToDead] << "\n"; + // std::cout << "sigmaid_2 = " + // << sim.get_Equilibirum_transitions()[1][(int)mio::endisird::InfectionTransition::InfectedToDead] << "\n"; + // std::cout << "sigmair_1 = " + // << sim.get_Equilibirum_transitions()[0][(int)mio::endisird::InfectionTransition::InfectedToRecovered] + // << "\n"; + // std::cout << "sigmair_2 = " + // << sim.get_Equilibirum_transitions()[1][(int)mio::endisird::InfectionTransition::InfectedToRecovered] + // << "\n"; + + // Uncomment to print the total population size. + // sim.get_totalpopulations().print_table({"N"}, 16, 9); +} \ No newline at end of file diff --git a/cpp/examples/ide_endemic_sird_save_results.cpp b/cpp/examples/ide_endemic_sird_save_results.cpp new file mode 100644 index 0000000000..39f37f5f91 --- /dev/null +++ b/cpp/examples/ide_endemic_sird_save_results.cpp @@ -0,0 +1,213 @@ +#include "memilio/io/result_io.h" +#include "memilio/utils/time_series.h" +#include "memilio/config.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/math/floating_point.h" + +#include "ide_endemic_sird/infection_state.h" +#include "ide_endemic_sird/model.h" +#include "ide_endemic_sird/simulation.h" +#include "ide_endemic_sird/parameters.h" + +#include "boost/filesystem.hpp" +#include +#include +#include + +mio::IOResult simulate_endidemodel(ScalarType tmax, std::string save_dir = "") +{ + using Vec = mio::TimeSeries::Vector; + + ScalarType dt = 1.0; + + int num_states = static_cast(mio::endisird::InfectionState::Count); + int num_transitions = static_cast(mio::endisird::InfectionTransition::Count); + + // Create TimeSeries with num_states elements where states needed for simulation will be stored. + mio::TimeSeries init(num_states); + + Vec vec_init(num_states); + + vec_init[static_cast(mio::endisird::InfectionState::Susceptible)] = 87500.; + vec_init[static_cast(mio::endisird::InfectionState::Infected)] = 190.; + vec_init[static_cast(mio::endisird::InfectionState::Recovered)] = 12310.; + vec_init[static_cast(mio::endisird::InfectionState::Dead)] = 0.; + + init.add_time_point(0, vec_init); + + mio::endisird::CompParameters computed_parameters(std::move(init)); + + //Set working parameters + + // mio::ExponentialSurvivalFunction exp(3.0); + // mio::StateAgeFunctionWrapper delaydistribution(exp); + // std::vector vec_delaydistribution(num_transitions, delaydistribution); + + mio::SmootherCosine smoothcos(8.0); + mio::StateAgeFunctionWrapper delaydistribution(smoothcos); + std::vector vec_delaydistribution(num_transitions, delaydistribution); + + // mio::ConstantFunction initialfunc(0); + // mio::StateAgeFunctionWrapper delaydistributioninit(initialfunc); + // std::vector vec_delaydistribution(num_transitions, delaydistributioninit); + // // InfectedToDead + // mio::SmootherCosine survivalExposedToInfectedNoSymptoms(14.0); + // vec_delaydistribution[(int)mio::endisird::InfectionTransition::InfectedToDead].set_state_age_function( + // survivalExposedToInfectedNoSymptoms); + // // InfectedToRecovered + // mio::SmootherCosine survivalInfectedNoSymptomsToInfectedSymptoms(8.0); + // vec_delaydistribution[(int)mio::endisird::InfectionTransition::InfectedToRecovered].set_state_age_function( + // survivalInfectedNoSymptomsToInfectedSymptoms); + + computed_parameters.parameters.get() = vec_delaydistribution; + + std::vector vec_prob((int)mio::endisird::InfectionTransition::Count, 1.); + vec_prob[Eigen::Index(mio::endisird::InfectionTransition::InfectedToDead)] = 0.1; + vec_prob[Eigen::Index(mio::endisird::InfectionTransition::InfectedToRecovered)] = 1 - 0.1; + + computed_parameters.parameters.set(vec_prob); + + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.)); + computed_parameters.parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + mio::ConstantFunction constant(0.1); + mio::StateAgeFunctionWrapper constant_prob(constant); + + computed_parameters.parameters.get() = constant_prob; + + mio::ExponentialSurvivalFunction exponential(0.5); + mio::StateAgeFunctionWrapper exponential_prob(exponential); + + computed_parameters.parameters.get() = exponential_prob; + + computed_parameters.parameters.set(4e-3); + computed_parameters.parameters.set(3e-3); + + //computed_parameters.set_tol_for_support_max(1e-6); + + mio::endisird::Model model(computed_parameters); + + mio::endisird::NormModel normmodel(computed_parameters); + + // start the simulation. + mio::endisird::Simulation sim(computed_parameters, model, normmodel, dt); + sim.advance(tmax); + + if (!save_dir.empty()) { + std::string tmax_string = std::to_string(tmax); + std::string dt_string = std::to_string(dt); + + std::string filename_ide = save_dir + "analysis_endide_SIRD_" + tmax_string.substr(0, tmax_string.find(".")) + + "_" + dt_string.substr(0, dt_string.find(".") + 5); + + //Save files of Model. + //Save total population. + std::string filename_ide_totalpopulation1 = filename_ide + "_totalpopulation.h5"; + mio::IOResult save_result_status_tp1 = + mio::save_result({sim.get_totalpopulations()}, {0}, 1, filename_ide_totalpopulation1); + //Save derivative of total population. + std::string filename_ide_dertotalpopulation1 = filename_ide + "_dertotalpopulation.h5"; + mio::IOResult save_result_status_dtp1 = + mio::save_result({sim.get_totalpopulations_derivative()}, {0}, 1, filename_ide_dertotalpopulation1); + + //Save compartments. + std::string filename_ide_compartments1 = filename_ide + "_compartments.h5"; + mio::IOResult save_result_status_c1 = + mio::save_result({sim.get_compartments()}, {0}, 1, filename_ide_compartments1); + + //Save normalized compartments. + std::string filename_ide_normcompartments = filename_ide + "_normcompartments.h5"; + mio::IOResult save_result_status_nc = + mio::save_result({sim.get_normalizedcompartments()}, {0}, 1, filename_ide_normcompartments); + + //Save the force of infection. + std::string filename_ide_forceofinfection = filename_ide + "_forceofinfection2.h5"; + mio::IOResult save_result_status_foi = + mio::save_result({sim.get_forceofinfections()}, {0}, 1, filename_ide_forceofinfection); + + //Save files of NormModel. + //Save compartments. + std::string filename_ide_normmod_compartments = filename_ide + "_normmod_compartments.h5"; + mio::IOResult save_result_status_nmc1 = + mio::save_result({sim.get_normmodel_compartments()}, {0}, 1, filename_ide_normmod_compartments); + + //Save the force of infection. + std::string filename_ide_normmod_forceofinfection = filename_ide + "_normmod_forceofinfection.h5"; + mio::IOResult save_result_status_nmfoi = + mio::save_result({sim.get_normmodel_forceofinfections()}, {0}, 1, filename_ide_normmod_forceofinfection); + + //Safe difference between normalized compartments. + std::string filename_ide_difference_normcomp = filename_ide + "_difference_normcomp.h5"; + mio::IOResult save_result_status_dnm = + mio::save_result({sim.get_difference_normalizationcomp()}, {0}, 1, filename_ide_difference_normcomp); + + //Safe difference between Force of Infections. + std::string filename_ide_difference_normFoI = filename_ide + "_difference_normfoi.h5"; + mio::IOResult save_result_status_dnf = + mio::save_result({sim.get_difference_normalizationFoi()}, {0}, 1, filename_ide_difference_normFoI); + + if (!save_result_status_tp1 || !save_result_status_c1 || !save_result_status_nc || !save_result_status_foi || + !save_result_status_nmc1 || !save_result_status_nmfoi || !save_result_status_dnm || + !save_result_status_dnf || !save_result_status_dtp1) { + return mio::failure(mio::StatusCode::UnknownError, "Error while saving results."); + } + } + // Print the reproduction number. + //Uncomment to print the reproduction number + std::cout << "The reproduction number Rc = " << sim.get_reproductionnumber_c() << "\n"; + + // Uncomment to print the the values T_z1^z2 + for (int i = 0; i < (int)sim.get_T().size(); i++) { + std::cout << "T_" << i << " = " << sim.get_T()[i] << "\n"; + } + + // Uncomment to print W_i + + std::cout << "W_i = " << sim.get_W() << "\n"; + + // Uncomment to print the equilibria + std::cout << "l_1 = " << sim.get_Equilibirum_FoI()[0] << "\n"; + std::cout << "l_2 = " << sim.get_Equilibirum_FoI()[1] << "\n"; + std::cout << "s_1 = " << sim.get_Equilibirum_compartments()[0][(int)mio::endisird::InfectionState::Susceptible] + << "\n"; + std::cout << "s_2 = " << sim.get_Equilibirum_compartments()[1][(int)mio::endisird::InfectionState::Susceptible] + << "\n"; + std::cout << "i_1 = " << sim.get_Equilibirum_compartments()[0][(int)mio::endisird::InfectionState::Infected] + << "\n"; + std::cout << "i_2 = " << sim.get_Equilibirum_compartments()[1][(int)mio::endisird::InfectionState::Infected] + << "\n"; + std::cout << "r_1 = " << sim.get_Equilibirum_compartments()[0][(int)mio::endisird::InfectionState::Recovered] + << "\n"; + std::cout << "r_2 = " << sim.get_Equilibirum_compartments()[1][(int)mio::endisird::InfectionState::Recovered] + << "\n"; + std::cout << "sigmaid_1 = " + << sim.get_Equilibirum_transitions()[0][(int)mio::endisird::InfectionTransition::InfectedToDead] << "\n"; + std::cout << "sigmaid_2 = " + << sim.get_Equilibirum_transitions()[1][(int)mio::endisird::InfectionTransition::InfectedToDead] << "\n"; + std::cout << "sigmair_1 = " + << sim.get_Equilibirum_transitions()[0][(int)mio::endisird::InfectionTransition::InfectedToRecovered] + << "\n"; + std::cout << "sigmair_2 = " + << sim.get_Equilibirum_transitions()[1][(int)mio::endisird::InfectionTransition::InfectedToRecovered] + << "\n"; + + return mio::success(); +} + +int main() +{ + std::string result_dir = "/localdata1/trit_ha/code/memilio-1/PythonPlotsEndIDESIRD/simulation_results/"; + + // Define tmax. + ScalarType tmax = 4000; + + auto result_ide = simulate_endidemodel(tmax, result_dir); + if (!result_ide) { + printf("%s\n", result_ide.error().formatted_message().c_str()); + return -1; + } + + return 0; +} \ No newline at end of file diff --git a/cpp/models/ide_endemic_secir/CMakeLists.txt b/cpp/models/ide_endemic_secir/CMakeLists.txt new file mode 100644 index 0000000000..b335f4c23c --- /dev/null +++ b/cpp/models/ide_endemic_secir/CMakeLists.txt @@ -0,0 +1,17 @@ +add_library(ide_endemic_secir + infection_state.h + model.h + model.cpp + normalized_model.h + normalized_model.cpp + simulation.h + simulation.cpp + parameters.h + computed_parameters.h +) +target_link_libraries(ide_endemic_secir PUBLIC memilio) +target_include_directories(ide_endemic_secir PUBLIC + $ + $ +) +target_compile_options(ide_endemic_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ide_endemic_secir/computed_parameters.h b/cpp/models/ide_endemic_secir/computed_parameters.h new file mode 100644 index 0000000000..a2bc681ec0 --- /dev/null +++ b/cpp/models/ide_endemic_secir/computed_parameters.h @@ -0,0 +1,493 @@ +#ifndef IDE_END_SECIR_COMPPARAMS_H +#define IDE_END_SECIR_COMPPARAMS_H + +#include "ide_endemic_secir/infection_state.h" +#include "ide_endemic_secir/parameters.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" + +#include "vector" +#include +#include +#include +#include +#include + +namespace mio +{ +namespace endisecir +{ +class CompParameters; +class Model; +class Simulation; + +class CompParameters +{ + using ParameterSet = Parameters; + +public: + CompParameters(TimeSeries&& states_init) + : parameters{Parameters()} + , m_statesinit{std::move(states_init)} + { + m_totalpopulationinit = std::accumulate(m_statesinit[0].begin(), m_statesinit[0].end(), 0) - + m_statesinit[0][(int)InfectionState::Dead]; + } + + bool check_constraints() const + { + if (!(static_cast(m_statesinit.get_num_elements()) == static_cast(InfectionState::Count))) { + log_error( + " A variable given for model construction is not valid. Number of elements in vector of populations " + "does not match the required number."); + return true; + } + for (int i = 0; i < static_cast(InfectionState::Count); i++) { + if (m_statesinit[0][i] < 0) { + log_error("Initialization failed. Initial values for populations are less than zero."); + return true; + } + } + return parameters.check_constraints(); + } + + /** + * @brief Setter for the tolerance used to calculate the maximum support of the TransitionDistributions. + * + * @param[in] new_tol New tolerance. + */ + void set_tol_for_support_max(ScalarType new_tol) + { + m_tol = new_tol; + } + + // ---- Public parameters. ---- + ParameterSet parameters; ///< ParameterSet of Model Parameters. + +private: + // ---- Functionality to set vectors with necessary information regarding TransitionDistributions. ---- + /** + * @brief Setter for the vector m_transitiondistributions_support_max that contains the support_max for all + * TransitionDistributions. + * + * This determines how many summands are required when calculating flows, the force of infection or compartments. + * + * @param[in] dt Time step size. + */ + void set_transitiondistributions_support_max(ScalarType dt) + { + // The transition Susceptible to Exposed is not needed in the computations. + for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) { + m_transitiondistributions_support_max[transition] = + parameters.get()[transition].get_support_max(dt, m_tol); + } + } + + /** + * @brief Setter for the vector m_transitiondistributions. + * + * Here we compute the weighted transition distributions for every compartment. Meaning for every compartment we weight + * the distributions of the outgoing transitions with their probability and compute the sum of the two weighed distributions. + * + * @param[in] dt Time step size. + */ + void set_transitiondistributions(ScalarType dt) + { + //Exposed state: + Eigen::Index support_max_index = (Eigen::Index)std::ceil( + m_transitiondistributions_support_max[(int)InfectionTransition::ExposedToInfectedNoSymptoms] / dt); + + std::vector vec_contribution_to_foi_1(support_max_index + 1, 0.); + for (Eigen::Index i = 0; i <= support_max_index; i++) { + ///Compute the state_age. + ScalarType state_age = (ScalarType)i * dt; + vec_contribution_to_foi_1[i] = + parameters.get()[(int)InfectionTransition::ExposedToInfectedNoSymptoms].eval( + state_age); + } + m_transitiondistributions[(int)InfectionState::Exposed] = vec_contribution_to_foi_1; + + // Vector containing the transitions for the states with two outgoing flows. + // We will compute the Transition Distribution for InfectedNoSymptoms and InfectedSymptoms until m_caltime, as we need + // that many values to compute the mean_infectivity. + + std::vector> vector_transitions = { + {(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms, + (int)InfectionTransition::InfectedNoSymptomsToRecovered}, + {(int)InfectionTransition::InfectedSymptomsToInfectedSevere, + (int)InfectionTransition::InfectedSymptomsToRecovered}, + {(int)InfectionTransition::InfectedSevereToInfectedCritical, + (int)InfectionTransition::InfectedSevereToRecovered}, + {(int)InfectionTransition::InfectedCriticalToDead, (int)InfectionTransition::InfectedCriticalToRecovered}}; + + for (int state = 2; state < (int)InfectionState::Count - 2; state++) { + Eigen::Index calc_time_index = std::ceil( + std::max(m_transitiondistributions_support_max[(Eigen::Index)vector_transitions[state - 2][0]], + m_transitiondistributions_support_max[(Eigen::Index)vector_transitions[state - 2][1]]) / + dt); + + std::vector vec_contribution_to_foi_2(calc_time_index + 1, 0.); + for (Eigen::Index i = 0; i <= calc_time_index; i++) { + ///Compute the state_age. + ScalarType state_age = (ScalarType)i * dt; + + vec_contribution_to_foi_2[i] += + parameters.get()[vector_transitions[state - 2][0]] * + parameters.get()[vector_transitions[state - 2][0]].eval(state_age) + + parameters.get()[vector_transitions[state - 2][1]] * + parameters.get()[vector_transitions[state - 2][1]].eval(state_age); + } + m_transitiondistributions[state] = vec_contribution_to_foi_2; + } + } + + /** + * @brief Setter for the vector m_transitiondistributions_derivative that contains the approximated derivative for + * all TransitionDistributions for all necessary time points. + * + * The derivative is approximated using a backwards difference scheme. + * The number of necessary time points for each TransitionDistribution is determined using + * m_transitiondistributions_support_max. + * + * @param[in] dt Time step size. + */ + void set_transitiondistributions_derivative(ScalarType dt) + { + // The transition SusceptibleToExposed 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 stores the value of the approximated derivative for all necessary time points. + std::vector vec_derivative(support_max_index + 1, 0.); + + for (Eigen::Index i = 0; i <= support_max_index; i++) { + // Compute state_age. + ScalarType state_age = (ScalarType)i * dt; + //Compute the apprximate derivative by a backwards difference scheme. + vec_derivative[i] = (parameters.get()[transition].eval(state_age) - + parameters.get()[transition].eval(state_age - dt)) / + dt; + } + m_transitiondistributions_derivative[transition] = vec_derivative; + } + } + + /** + *@brief Setter for the vector m_B that contains the approximated value of B for all for all necessary time points. + * + * The value m_B is needed for the compuation of the force of infection term and also for the compuation of + * m_infectivity. + * + * @param[in] dt Time step size. + */ + void set_B(ScalarType dt) + { + // Determine the calc_time_index that is the maximum of the support max of the TransitionDistributions that are + // used in this computation. + ScalarType calc_time = + std::max(m_transitiondistributions_support_max[(int)InfectionTransition::InfectedSymptomsToInfectedSevere], + m_transitiondistributions_support_max[(int)InfectionTransition::InfectedSymptomsToRecovered]); + Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1; + + m_B = std::vector(calc_time_index + 1, 0.); + + // We start the for loop for time_point_index = 1, as the next for loop where we compute the sum, goes up to + // time_point_index-1. Therefore, for time_point_index = 0, the for loop would start at 0 and go up to -1, which + // is not possible. The value m_B(0) is just set 0 zero. + for (Eigen::Index time_point_index = 1; time_point_index <= calc_time_index; time_point_index++) { + + ScalarType sum = 0; + Eigen::Index max_support_index = std::ceil( + m_transitiondistributions_support_max[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] / + dt); + Eigen::Index starting_point = std::max(0, (int)time_point_index - (int)max_support_index); + + for (int i = starting_point; i <= time_point_index - 1; i++) { + ScalarType state_age_i = static_cast(i) * dt; + + sum += + parameters.get().eval(state_age_i) * + m_transitiondistributions[static_cast(InfectionState::InfectedSymptoms)][i] * + m_transitiondistributions_derivative[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] + [time_point_index - i]; + } + m_B[time_point_index] = + parameters + .get()[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] * + dt * sum; + } + } + + /** + *@brief Setter for the vector m_meaninfectivity that contains the approximated value of the mean infectivity for all + * for all necessary time points. + * + * This values is needed the compute the reproduction numer and the force of infection term. + * + * @param[in] dt Time step size. + */ + void set_infectivity(ScalarType dt) + { + // Compute the calc_time_indedx corresponding to a_1: + ScalarType calc_time_1 = std::max( + m_transitiondistributions_support_max[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms], + m_transitiondistributions_support_max[(int)InfectionTransition::InfectedNoSymptomsToRecovered]); + Eigen::Index calc_time_index_1 = (Eigen::Index)std::ceil(calc_time_1 / dt) - 1; + + // Compute the calc_time_indedx corresponding to a_2: + Eigen::Index calc_time_index_2 = m_B.size() - 1; + + // All in all calc_time_index for the infectivity: + Eigen::Index calc_time_index = std::max(calc_time_index_1 + 1, calc_time_index_2 + 1); + + m_infectivity = std::vector(calc_time_index, 0.); + + for (Eigen::Index time_point_index = 1; time_point_index < calc_time_index; time_point_index++) { + ScalarType a_1 = 0; + ScalarType a_2 = 0; + + ScalarType time_point_age = (ScalarType)time_point_index * dt; + //We compute a_1 and a_2 at time_point. + Eigen::Index max_support_index = std::ceil( + m_transitiondistributions_support_max[(int)InfectionTransition::ExposedToInfectedNoSymptoms] / dt); + Eigen::Index starting_point = std::max(0, (int)time_point_index - (int)max_support_index); + + //compute a_1: + + for (int i = starting_point; i < std::min(calc_time_index_1, time_point_index); i++) { + + ScalarType state_age_i = static_cast(i) * dt; + a_1 += parameters.get().eval(state_age_i) * + m_transitiondistributions[static_cast(InfectionState::InfectedNoSymptoms)][i] * + m_transitiondistributions_derivative[(int)InfectionTransition::ExposedToInfectedNoSymptoms] + [time_point_index - i]; + } + //compute a_2: + + for (int i = starting_point; i < std::min(calc_time_index_2, time_point_index); i++) { + + a_2 += m_transitiondistributions_derivative[(int)InfectionTransition::ExposedToInfectedNoSymptoms] + [time_point_index - i] * + m_B[i]; + } + m_infectivity[time_point_index] = + dt * std::exp(-parameters.get() * time_point_age) * (a_2 - a_1); + } + } + + /** + *@brief Setter for the vectors m_FoI_0 and m_NormFoI_0 that contain the approximated values of the function FoI_0, + * that is a part of the force of infection term for all necessary time points. + */ + void set_FoI_0(ScalarType dt) + { + Eigen::Index calc_time_index = m_transitiondistributions[(int)InfectionState::InfectedNoSymptoms].size(); + m_FoI_0 = std::vector(calc_time_index, 0.); + m_NormFoI_0 = std::vector(calc_time_index, 0.); + for (Eigen::Index time_point_index = 0; time_point_index < calc_time_index; time_point_index++) { + ScalarType time_point_age = (ScalarType)time_point_index * dt; + m_FoI_0[time_point_index] = + (parameters.get().eval(time_point_age) * + parameters.get().get_cont_freq_mat().get_matrix_at(time_point_age)(0, 0)) * + (m_statesinit[0][(int)InfectionState::InfectedNoSymptoms] * + m_transitiondistributions[(int)InfectionState::InfectedNoSymptoms][time_point_index] * + std::exp(-parameters.get() * time_point_age) * + parameters.get().eval(time_point_age) + + m_statesinit[0][(int)InfectionState::InfectedSymptoms] * + m_transitiondistributions[(int)InfectionState::InfectedSymptoms][time_point_index] * + std::exp(-parameters.get() * time_point_age) * + parameters.get().eval(time_point_age)); + m_NormFoI_0[time_point_index] = + (parameters.get().eval(time_point_age) * + parameters.get().get_cont_freq_mat().get_matrix_at(time_point_age)(0, 0)) * + ((m_statesinit[0][(int)InfectionState::InfectedNoSymptoms] * + m_transitiondistributions[(int)InfectionState::InfectedNoSymptoms][time_point_index] * + std::exp(-parameters.get() * time_point_age) * + parameters.get().eval(time_point_age) + + m_statesinit[0][(int)InfectionState::InfectedSymptoms] * + m_transitiondistributions[(int)InfectionState::InfectedSymptoms][time_point_index] * + std::exp(-parameters.get() * time_point_age) * + parameters.get().eval(time_point_age)) / + m_totalpopulationinit); + } + } + + /** + * @brief Setter for the vectors m_InitFoI and m_NormInitFoI that contain the approximated values of the functions + * f(standard model) and g(normalized model) that is a part of the force of infection term for all necessary time + * points. + * + * @param[in] dt Time step size. + */ + void set_InitFoI(ScalarType dt) + { + Eigen::Index calc_time_index = std::max(m_infectivity.size(), m_B.size()); + m_InitFoI.resize(calc_time_index, 0.); + m_NormInitFoI.resize(calc_time_index, 0.); + for (Eigen::Index time_point_index = 0; time_point_index < calc_time_index; time_point_index++) { + ScalarType sum = 0; + ScalarType sum_norm = 0; + ScalarType time_point_age = (ScalarType)time_point_index * dt; + if (time_point_index < (int)m_B.size()) { + sum -= m_statesinit[0][(int)InfectionState::InfectedNoSymptoms] * + std::exp(-parameters.get() * time_point_age) * m_B[time_point_index]; + sum_norm -= m_statesinit[0][(int)InfectionState::InfectedNoSymptoms] / m_totalpopulationinit * + std::exp(-parameters.get() * time_point_age) * m_B[time_point_index]; + } + if (time_point_index < (int)m_infectivity.size()) { + sum += m_statesinit[0][(int)InfectionState::Exposed] * m_infectivity[time_point_index]; + sum_norm += m_statesinit[0][(int)InfectionState::Exposed] / m_totalpopulationinit * + m_infectivity[time_point_index]; + } + m_InitFoI[time_point_index] = + parameters.get().eval(time_point_age) * + parameters.get().get_cont_freq_mat().get_matrix_at(time_point_age)(0, 0) * sum; + m_NormInitFoI[time_point_index] = + parameters.get().eval(time_point_age) * + parameters.get().get_cont_freq_mat().get_matrix_at(time_point_age)(0, 0) * sum_norm; + } + } + + // ---- Parameters needed for the analysis of the model ---- + /** + * @brief Setter for the Reproduction number m_reproductionnumber_c. + * + */ + void set_reproductionnumber_c(ScalarType dt) + { + // Determine the corresponding time index. + + Eigen::Index calc_time_index = m_infectivity.size() - 1; + ScalarType sum = 0; + for (int i = 0; i <= calc_time_index; i++) { + sum += m_infectivity[i]; + } + m_reproductionnumber_c = parameters.get().eval(0) * + parameters.get().get_cont_freq_mat().get_matrix_at(0)(0, 0) * dt * + sum; + } + + /** + * @brief Setter for m_T. + * + * @param[in] dt step size. + */ + void set_T(ScalarType dt) + { + // The value T_z1^z2 is not defined for the transition SusceptiblesToExposed. + for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) { + Eigen::Index support_max_index = + (Eigen::Index)std::ceil(m_transitiondistributions_support_max[transition] / dt); + ScalarType sum = 0; + + for (Eigen::Index i = 0; i <= support_max_index; i++) { + // Compute state_age. + ScalarType state_age = (ScalarType)i * dt; + //Compute the apprximate derivative by a backwards difference scheme. + sum -= std::exp(-parameters.get() * state_age) * + m_transitiondistributions_derivative[transition][i]; + } + m_T[transition] = parameters.get()[transition] * sum; + } + } + + /** + * @brief Setter for m_V. + * + * @param[in] dt step size. + */ + void set_V() + { + // The value V^z is not defined for the compartments Susceptibles and Exposed + // InfectedNoSymptoms: + m_V[(int)InfectionState::InfectedNoSymptoms] = m_T[(int)InfectionTransition::ExposedToInfectedNoSymptoms]; + // InfectedSymptoms: + m_V[(int)InfectionState::InfectedSymptoms] = + m_T[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] * + m_V[(int)InfectionState::InfectedNoSymptoms]; + // InfectedSevere: + m_V[(int)InfectionState::InfectedSevere] = m_T[(int)InfectionTransition::InfectedSymptomsToInfectedSevere] * + m_V[(int)InfectionState::InfectedSymptoms]; + // InfectedCritical: + m_V[(int)InfectionState::InfectedCritical] = + m_T[(int)InfectionTransition::InfectedSevereToInfectedCritical] * m_V[(int)InfectionState::InfectedSevere]; + // Dead: + m_V[(int)InfectionState::Dead] = + m_T[(int)InfectionTransition::InfectedCriticalToDead] * m_V[(int)InfectionState::InfectedCritical]; + // Recovered: + m_V[(int)InfectionState::Recovered] = + m_T[(int)InfectionTransition::InfectedNoSymptomsToRecovered] * + m_V[(int)InfectionState::InfectedNoSymptoms] + + m_T[(int)InfectionTransition::InfectedSymptomsToRecovered] * m_V[(int)InfectionState::InfectedSymptoms] + + m_T[(int)InfectionTransition::InfectedSevereToRecovered] * m_V[(int)InfectionState::InfectedSevere] + + m_T[(int)InfectionTransition::InfectedCriticalToRecovered] * m_V[(int)InfectionState::InfectedCritical]; + } + + /** + * @brief Setter for m_W. + * + * m_W contains the values W_z for compartments z. As W_z is only defined for the compartments Exposed, InfectedNoSymptoms, + * InfectedSymptoms, InfectedSevere and InfectedCrititcal we will set W_z for the other compartments to zero. + * + * @param[in] dt step size. + */ + void set_W(ScalarType dt) + { + for (int state = 1; state < (int)InfectionState::Count - 2; state++) { + Eigen::Index calc_time_index = m_transitiondistributions[state].size(); + ScalarType sum = 0; + + for (Eigen::Index i = 0; i < calc_time_index; i++) { + // Compute state_age. + ScalarType state_age = (ScalarType)i * dt; + //Compute the apprximate derivative by a backwards difference scheme. + sum += std::exp(-parameters.get() * state_age) * m_transitiondistributions[state][i]; + } + m_W[state] = sum; + } + } + + // ---- Private parameters. ---- + TimeSeries m_statesinit; ///< TimeSeries containing the initial values for the compartments. + ScalarType m_totalpopulationinit; + ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the TransitionDistributions. + std::vector m_transitiondistributions_support_max{ + std::vector((int)InfectionTransition::Count, 0.)}; ///< A vector containing the support_max + // for all TransitionDistributions. + std::vector> m_transitiondistributions{std::vector>( + (int)InfectionState::Count, + std::vector(1, 0.))}; ///< A vector containing the weighted TransitionDistributions. + std::vector> m_transitiondistributions_derivative{std::vector>( + (int)InfectionTransition::Count, std::vector(1, 0.))}; ///< A Vector containing + // the approximated derivative for all TransitionDistributions for all necessary time points. + std::vector m_B{std::vector(1, 0.)}; ///< A Vector contaiing the appriximated + // values for B. + std::vector m_infectivity{ + std::vector(1, 0.)}; ///< A vector containing the approximated mean infectivity for all time points. + ScalarType m_reproductionnumber_c; ///< The control Reproduction number + std::vector m_FoI_0{std::vector(1, 0.)}; ///< A vector containing the approcimated + // value of the function FoI_0 used for the computation of the force of infection in the standard model + std::vector m_NormFoI_0{std::vector(1, 0.)}; ///< A vector containing the approcimated + // value of the function FoI_0 used for the computation of the force of infection in the normalized model + std::vector m_InitFoI = std::vector(1, 0.); ///< A vector containing the approcimated + // value of the function FoI_0 used for the computation of the force of infection in the standard model + std::vector m_NormInitFoI{std::vector(1, 0.)}; ///< A vector containing the approcimated + // value of the function FoI_0 used for the computation of the force of infection in the normalized model + std::vector m_T{std::vector((int)InfectionTransition::Count, 0.)}; ///< A vector + // containing the approximated value for T_z1^z2 for every Flow z1 to z2. + std::vector m_V{std::vector((int)InfectionTransition::Count, 0.)}; ///< A vector + // containing the approximated value for V^z for every compartment z. + std::vector m_W{std::vector((int)InfectionState::Count, 0.)}; ///< A vector containing+ + // the approximated value for W_z for every compartment z. + // ---- Friend classes/functions. ---- + friend class Model; + friend class NormModel; + friend class Simulation; +}; // namespace endisecir + +} // namespace endisecir + +} // namespace mio + +#endif //IDE_END_SECIR_COMPPARAMS_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_secir/infection_state.h b/cpp/models/ide_endemic_secir/infection_state.h new file mode 100644 index 0000000000..e2e7563d40 --- /dev/null +++ b/cpp/models/ide_endemic_secir/infection_state.h @@ -0,0 +1,49 @@ +#ifndef IDE_END_SECIR_INFECTIONSTATE_H +#define IDE_END_SECIR_INFECTIONSTATE_H + +namespace mio +{ + +namespace endisecir +{ + +/** + * @brief The #InfectionState enum describes the possible + * categories for the infectious state of persons + */ +enum class InfectionState +{ + Susceptible = 0, + Exposed = 1, + InfectedNoSymptoms = 2, + InfectedSymptoms = 3, + InfectedSevere = 4, + InfectedCritical = 5, + Recovered = 6, + Dead = 7, + Count = 8 +}; + +/** + * @brief The #InfectionTransition enum describes the possible + * transitions of the infectious state of persons. + */ +enum class InfectionTransition +{ + SusceptibleToExposed = 0, + ExposedToInfectedNoSymptoms = 1, + InfectedNoSymptomsToInfectedSymptoms = 2, + InfectedNoSymptomsToRecovered = 3, + InfectedSymptomsToInfectedSevere = 4, + InfectedSymptomsToRecovered = 5, + InfectedSevereToInfectedCritical = 6, + InfectedSevereToRecovered = 7, + InfectedCriticalToDead = 8, + InfectedCriticalToRecovered = 9, + Count = 10 +}; + +} // namespace endisecir +} // namespace mio + +#endif //IDE_END_SECIR_INFECTIONSTATE_H diff --git a/cpp/models/ide_endemic_secir/model.cpp b/cpp/models/ide_endemic_secir/model.cpp new file mode 100644 index 0000000000..b5031ae015 --- /dev/null +++ b/cpp/models/ide_endemic_secir/model.cpp @@ -0,0 +1,412 @@ +#include "ide_endemic_secir/model.h" +#include "ide_endemic_secir/computed_parameters.h" +#include "ide_endemic_secir/parameters.h" +#include "ide_endemic_secir/infection_state.h" +#include "memilio/config.h" +#include "memilio/utils/logging.h" + +#include "memilio/utils/time_series.h" +#include "vector" +#include +#include +#include +#include +#include +#include +#include + +namespace mio +{ +namespace endisecir +{ +Model::Model(CompParameters const& compparams) + : compparameters{std::make_shared(compparams)} + , transitions{TimeSeries(Eigen::Index(InfectionTransition::Count))} + , transitions_update{TimeSeries(Eigen::Index(InfectionTransition::Count))} + , populations{compparameters->m_statesinit} + , populations_update{compparameters->m_statesinit} + +{ + + // Set flows at start time t0. + // As we assume that all individuals have infectio age 0 at time t0, the flows at t0 are set to 0. + transitions.add_time_point( + 0, TimeSeries::Vector::Constant(static_cast(InfectionTransition::Count), 0.)); + transitions_update.add_time_point( + 0, TimeSeries::Vector::Constant(static_cast(InfectionTransition::Count), 0.)); + + // Set population size at start timt t0. + ScalarType init_populationsize = + std::accumulate(populations[0].begin(), populations[0].end(), 0) - populations[0][(int)InfectionState::Dead]; + m_totalpopulation.add_time_point(0, TimeSeries::Vector::Constant(1, init_populationsize)); + m_totalpopulationupdate.add_time_point(0, TimeSeries::Vector::Constant(1, init_populationsize)); + m_totalpopulation_derivative.add_time_point(0, TimeSeries::Vector::Constant(1, 0.)); + + // Set the force of infection term at time t0. + m_forceofinfection.add_time_point(0, TimeSeries::Vector::Constant(1, compparameters->m_FoI_0[0])); + m_forceofinfectionupdate.add_time_point(0, TimeSeries::Vector::Constant(1, compparameters->m_FoI_0[0])); + + //Set normalized_populations at start time t0. + TimeSeries::Vector vec_normalizedpopulations = + TimeSeries::Vector(Eigen::Index(InfectionState::Count) - 1); + for (int infection_state = 0; infection_state < Eigen::Index(InfectionState::Count) - 1; infection_state++) { + vec_normalizedpopulations[infection_state] = populations[0][infection_state] / m_totalpopulation[0][0]; + } + m_normalizedpopulations.add_time_point(0, vec_normalizedpopulations); +} + +bool Model::check_constraints() const +{ + if (!(static_cast(populations.get_num_elements()) == static_cast(InfectionState::Count))) { + log_error(" A variable given for model construction is not valid. Number of elements in vector of populations " + "does not match the required number."); + return true; + } + + for (int i = 0; i < static_cast(InfectionState::Count); i++) { + if (populations[0][i] < 0) { + log_error("Initialization failed. Initial values for populations are less than zero."); + return true; + } + } + return compparameters->check_constraints(); +} + +// ----Functionality for the iterations of a simulation. ---- +void Model::compute_susceptibles(ScalarType dt) +{ + Eigen::Index num_time_points = populations.get_num_time_points(); + populations.get_last_value()[static_cast(InfectionState::Susceptible)] = + (populations[num_time_points - 2][static_cast(InfectionState::Susceptible)] + + dt * m_totalpopulation[num_time_points - 2][0] * compparameters->parameters.get()) / + (1 + dt * (m_forceofinfection[num_time_points - 2][0] + compparameters->parameters.get())); + + //Computation of susceptibles using the update formula for the other compartments. + //Formula for Susceptibles stays the same, but we use m_totalpopulationupdate and m_forceofinfectionupdate + populations_update.get_last_value()[static_cast(InfectionState::Susceptible)] = + (populations_update[num_time_points - 2][static_cast(InfectionState::Susceptible)] + + dt * (m_totalpopulationupdate[num_time_points - 2][0] * compparameters->parameters.get())) / + (1 + + dt * (m_forceofinfectionupdate[num_time_points - 2][0] + compparameters->parameters.get())); +} + +void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt) +{ + Eigen::Index calc_time_index = + (Eigen::Index)std::ceil(compparameters->m_transitiondistributions_support_max[idx_InfectionTransitions] / dt); + ScalarType current_time_age = static_cast(current_time_index) * dt; + + ScalarType sum1 = 0; + ScalarType sum2 = 0; + //Determine the starting point of the for loop. + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + + for (Eigen::Index i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = static_cast(i) * dt; + sum1 += transitions[i + 1][idx_IncomingFlow] * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)) * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index - i]; + //For the update formula version of the model: + sum2 += transitions_update[i + 1][idx_IncomingFlow] * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)) * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index - i]; + } + if (current_time_index <= calc_time_index) { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum1 - + std::exp((-compparameters->parameters.get()) * (current_time_age)) * + populations[0][idx_CurrentCompartment] * + compparameters->parameters.get()[idx_InfectionTransitions] * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index]; + + //For the update formula version of the model: + transitions_update.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum2 - + std::exp((-compparameters->parameters.get()) * (current_time_age)) * + populations_update[0][idx_CurrentCompartment] * + compparameters->parameters.get()[idx_InfectionTransitions] * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index]; + } + else { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum1; + + //For the update formula version of the model: + transitions_update.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum2; + } +} + +void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt) +{ + Eigen::Index current_time_index = transitions.get_num_time_points() - 1; + compute_flow(idx_InfectionTransitions, idx_IncomingFlow, idx_CurrentCompartment, current_time_index, dt); +} + +void Model::flows_currents_timestep(ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + // Calculate the transition SusceptibleToExposed with force of infection from previous time step and Susceptibles from + // current time step. + transitions.get_last_value()[static_cast(InfectionTransition::SusceptibleToExposed)] = + m_forceofinfection[current_time_index - 1][0] * + populations.get_last_value()[static_cast(InfectionState::Susceptible)]; + + ///For the update formula version of the model: + transitions_update.get_last_value()[static_cast(InfectionTransition::SusceptibleToExposed)] = + m_forceofinfectionupdate[current_time_index - 1][0] * + populations_update.get_last_value()[static_cast(InfectionState::Susceptible)]; + + // Calculate the other Transitions with compute_flow. + // Exposed to InfectedNoSymptoms + compute_flow(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), + Eigen::Index(InfectionTransition::SusceptibleToExposed), Eigen::Index(InfectionState::Exposed), dt); + // InfectedNoSymptoms to InfectedSymptoms + compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), + Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), + Eigen::Index(InfectionState::InfectedNoSymptoms), dt); + // InfectedNoSymptoms to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered), + Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), + Eigen::Index(InfectionState::InfectedNoSymptoms), dt); + // InfectedSymptoms to InfectedSevere + compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), + Eigen::Index(InfectionState::InfectedSymptoms), dt); + // InfectedSymptoms to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered), + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), + Eigen::Index(InfectionState::InfectedSymptoms), dt); + // InfectedSevere to InfectedCritical + compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + Eigen::Index(InfectionState::InfectedSevere), dt); + // InfectedCritical to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered), + Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + Eigen::Index(InfectionState::InfectedSevere), dt); + // InfectedCritical to Dead + compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead), + Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + Eigen::Index(InfectionState::InfectedCritical), dt); + // InfectedCritical to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered), + Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + Eigen::Index(InfectionState::InfectedCritical), dt); +} + +void Model::update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, + bool NaturalDeathispossible, bool Transitionispossible, ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + ScalarType current_time_age = (ScalarType)current_time_index * dt; + Eigen::Index calc_time_index = current_time_index; + if (Transitionispossible) { + calc_time_index = compparameters->m_transitiondistributions[(int)infectionState].size() - 1; + } + + ScalarType sum = 0; + + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + for (int i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = (ScalarType)i * dt; + ScalarType sum_inflows = 0; + for (const InfectionTransition& inflow : IncomingFlows) { + sum_inflows += transitions[i + 1][(int)inflow]; + } + if (NaturalDeathispossible && Transitionispossible) { + sum += compparameters->m_transitiondistributions[(int)infectionState][current_time_index - i] * + sum_inflows * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + else if (NaturalDeathispossible && !Transitionispossible) { + sum += sum_inflows * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + // The case !NaturalDeathispossible && Transitionispossible is not possible, as if NaturalDeath is not possible + // this means you are in the Death compartment and then Transition is also not possible. + else { + sum += sum_inflows; + } + } + if (NaturalDeathispossible && Transitionispossible) { + if (current_time_index <= calc_time_index) { + populations.get_last_value()[(int)infectionState] = + dt * sum + compparameters->m_transitiondistributions[(int)infectionState][current_time_index] * + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } + else { + populations.get_last_value()[(int)infectionState] = dt * sum; + } + } + else if (NaturalDeathispossible && !Transitionispossible) { + populations.get_last_value()[(int)infectionState] = + dt * sum + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } + else { + populations.get_last_value()[(int)infectionState] = dt * sum + populations[0][(int)infectionState]; + } +} +void Model::update_compartment_from_flow(InfectionState infectionState, + std::vector const& IncomingFlows, + std::vector const& OutgoingFlows, + bool NaturalDeathispossible, ScalarType dt) +{ + Eigen::Index num_time_points = populations.get_num_time_points(); + + ScalarType updated_compartment = populations_update[num_time_points - 2][static_cast(infectionState)]; + + for (const InfectionTransition& inflow : IncomingFlows) { + updated_compartment += dt * transitions_update.get_last_value()[(int)inflow]; + } + for (const InfectionTransition& outflow : OutgoingFlows) { + updated_compartment -= dt * transitions_update.get_last_value()[(int)outflow]; + } + + if (NaturalDeathispossible) { + updated_compartment = updated_compartment / (1 + dt * compparameters->parameters.get()); + } + populations_update.get_last_value()[(int)infectionState] = updated_compartment; +} + +void Model::update_compartments(ScalarType dt) +{ + + // Exposed + update_compartment_from_flow(InfectionState::Exposed, {InfectionTransition::SusceptibleToExposed}, + {InfectionTransition::ExposedToInfectedNoSymptoms}, true, dt); + update_compartment_with_sum(InfectionState::Exposed, {InfectionTransition::SusceptibleToExposed}, true, true, dt); + // InfectedNoSymptoms + update_compartment_from_flow( + InfectionState::InfectedNoSymptoms, {InfectionTransition::ExposedToInfectedNoSymptoms}, + {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms, InfectionTransition::InfectedNoSymptomsToRecovered}, + true, dt); + update_compartment_with_sum(InfectionState::InfectedNoSymptoms, {InfectionTransition::ExposedToInfectedNoSymptoms}, + true, true, dt); + + // InfectedSymptoms + update_compartment_from_flow( + InfectionState::InfectedSymptoms, {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms}, + {InfectionTransition::InfectedSymptomsToInfectedSevere, InfectionTransition::InfectedSymptomsToRecovered}, true, + dt); + update_compartment_with_sum(InfectionState::InfectedSymptoms, + {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms}, true, true, dt); + + // InfectedSevere + update_compartment_from_flow( + InfectionState::InfectedSevere, {InfectionTransition::InfectedSymptomsToInfectedSevere}, + {InfectionTransition::InfectedSevereToInfectedCritical, InfectionTransition::InfectedSevereToRecovered}, true, + dt); + update_compartment_with_sum(InfectionState::InfectedSevere, {InfectionTransition::InfectedSymptomsToInfectedSevere}, + true, true, dt); + + // InfectedCritical + update_compartment_from_flow( + InfectionState::InfectedCritical, {InfectionTransition::InfectedSevereToInfectedCritical}, + {InfectionTransition::InfectedCriticalToDead, InfectionTransition::InfectedCriticalToRecovered}, true, dt); + update_compartment_with_sum(InfectionState::InfectedCritical, + {InfectionTransition::InfectedSevereToInfectedCritical}, true, true, dt); + // Recovered + update_compartment_from_flow(InfectionState::Recovered, + { + InfectionTransition::InfectedNoSymptomsToRecovered, + InfectionTransition::InfectedSymptomsToRecovered, + InfectionTransition::InfectedSevereToRecovered, + InfectionTransition::InfectedCriticalToRecovered, + }, + std::vector(), true, dt); + update_compartment_with_sum(InfectionState::Recovered, + { + InfectionTransition::InfectedNoSymptomsToRecovered, + InfectionTransition::InfectedSymptomsToRecovered, + InfectionTransition::InfectedSevereToRecovered, + InfectionTransition::InfectedCriticalToRecovered, + }, + true, false, dt); + + // Dead + update_compartment_from_flow(InfectionState::Dead, {InfectionTransition::InfectedCriticalToDead}, + std::vector(), false, dt); + update_compartment_with_sum(InfectionState::Dead, {InfectionTransition::InfectedCriticalToDead}, false, false, dt); +} + +void Model::compute_populationsize() +{ + ScalarType sum1 = 0; + ScalarType sum2 = 0; + for (int state = 0; state < Eigen::Index(InfectionState::Count) - 1; state++) { + sum1 += populations.get_last_value()[state]; + sum2 += populations_update.get_last_value()[state]; + } + m_totalpopulation.get_last_value()[0] = sum1; + m_totalpopulationupdate.get_last_value()[0] = sum2; + // Here we comoute the derivative of the total population that is given by + // BirthRate * N(t) - DeathRate * N(t) - Transition[InfectedCritical>ToDeath](t). + m_totalpopulation_derivative.get_last_value()[0] = + (compparameters->parameters.get() - compparameters->parameters.get()) * + m_totalpopulation.get_last_value()[0] - + transitions.get_last_value()[(int)InfectionTransition::InfectedCriticalToDead]; +} + +void Model::compute_normalizedcompartments() +{ + for (int infection_state = 0; infection_state < Eigen::Index(InfectionState::Count) - 1; infection_state++) { + m_normalizedpopulations.get_last_value()[infection_state] = + populations.get_last_value()[infection_state] / m_totalpopulation.get_last_value()[0]; + } +} + +void Model::compute_forceofinfection(ScalarType dt) +{ + + Eigen::Index num_time_points = populations.get_num_time_points(); + ScalarType current_time = populations.get_last_time(); + + // Determine the starting point of the for loop. + Eigen::Index starting_point = std::max(0, (int)num_time_points - 1 - (int)compparameters->m_infectivity.size()); + + ScalarType sum = 0.; + ScalarType sum_update = 0.; + // Compute the sum in the force of infection term. + for (Eigen::Index i = starting_point; i < num_time_points - 1; i++) { + sum += compparameters->m_infectivity[num_time_points - 1 - i] * + populations[i + 1][(int)InfectionState::Susceptible] * m_forceofinfection[i][0]; + sum_update += compparameters->m_infectivity[num_time_points - 1 - i] * + populations_update[i + 1][(int)InfectionState::Susceptible] * m_forceofinfectionupdate[i][0]; + } + + ScalarType Foi_0 = 0; + ScalarType Foi_0_update = 0; + ScalarType f = 0; + ScalarType f_update = 0; + + // Add inital functions for the force of infection in case they still have an influence. + if (num_time_points <= (int)compparameters->m_FoI_0.size()) { + Foi_0 = compparameters->m_FoI_0[num_time_points - 1] / m_totalpopulation[num_time_points - 1][0]; + Foi_0_update = compparameters->m_FoI_0[num_time_points - 1] / m_totalpopulationupdate[num_time_points - 1][0]; + } + if (num_time_points <= (int)compparameters->m_InitFoI.size()) { + f = compparameters->m_InitFoI[num_time_points - 1] / m_totalpopulation[num_time_points - 1][0]; + f_update = compparameters->m_InitFoI[num_time_points - 1] / m_totalpopulationupdate[num_time_points - 1][0]; + } + m_forceofinfection.get_last_value()[0] = + (dt * sum * compparameters->parameters.get().eval(current_time) * + compparameters->parameters.get().get_cont_freq_mat().get_matrix_at(current_time)(0, 0)) / + m_totalpopulation[num_time_points - 1][0] + + Foi_0 + f; + m_forceofinfectionupdate.get_last_value()[0] = + (dt * sum_update * + (compparameters->parameters.get().eval(current_time) * + compparameters->parameters.get().get_cont_freq_mat().get_matrix_at(current_time)(0, 0))) / + m_totalpopulationupdate[num_time_points - 1][0] + + Foi_0_update + f_update; +} + +}; // namespace endisecir + +} // namespace mio diff --git a/cpp/models/ide_endemic_secir/model.h b/cpp/models/ide_endemic_secir/model.h new file mode 100644 index 0000000000..82e4e61043 --- /dev/null +++ b/cpp/models/ide_endemic_secir/model.h @@ -0,0 +1,212 @@ +#ifndef IDE_END_SECIR_MODEL_H +#define IDE_END_SECIR_MODEL_H + +#include "ide_endemic_secir/infection_state.h" +#include "ide_endemic_secir/parameters.h" +#include "ide_endemic_secir/computed_parameters.h" +#include "memilio/config.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/time_series.h" + +#include "vector" +#include +#include +#include + +namespace mio +{ +namespace endisecir +{ +// Forward declaration of friend classes/functions of Model. +class Model; +class Simulation; + +class Model +{ + using ParameterSet = Parameters; + +public: + /** + * @brief Constructor to create an endemic IDE_SECIR model. + * + * @param[in] TODO!!!!!! + */ + + Model(CompParameters const& compparams); + + /** + * @brief Checks constraints on model parameters and initial data. + * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. + */ + bool check_constraints() const; + + /** + * @brief Setter for the tolerance used to calculate the maximum support of the TransitionDistributions. + * + * @param[in] new_tol New tolerance. + */ + + // ---- Public parameters. ---- + std::shared_ptr compparameters; + TimeSeries + transitions; ///< TimesSeries containing points of time and the corresponding number of individuals transitioning + // from one #InfectionState to another as defined in #Infection%s. + TimeSeries + transitions_update; ///< TimesSeries containing points of time and the corresponding number of individuals transitioning + // from one #InfectionState to another as defined in #Infection%s. In this case we use the update formula version. + TimeSeries populations; ///< TimeSeries containing points of time and the corresponding number of people + // in defined #InfectionState%s. In this case we compute them by a sum. + TimeSeries + populations_update; ///< TimeSeries containing points of time and the corresponding number of people + // in defined #InfectionState%s. We compute them by an update formula. + +private: + // ---- Functionality for the iterations of a simulation. + + /** + * @brief Computes number of Susceptibles for the current last time in populations. + * + * Number is computed by using the previous number of Susceptibles, total population of the last time point and + * the force of infection (also from the last time point). + * Number is stored at the matching index in populations and populations_update. + * + * @param[in] dt Time discretization step size. + */ + void compute_susceptibles(ScalarType dt); + + /** + * @brief Computes sizeof a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the time index + * current_time_index. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] current_time_index The time index the transition should be computed for. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt); + + /** + * @brief Computes size of a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the current + * last time value in transitions. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt); + + /** + * @brief Sets all required transitions for the current last timestep in transitions. + * + * New values are stored in transitions. Most values are computed via the function compute_flow() + * + * @param[in] dt time discretization step size. + */ + void flows_currents_timestep(ScalarType dt); + + /** + * @brief Updates the values of one compartment, specified in infectionState, using all past transitions. + * + * New value is stored in populations. The values are calculated using all past values for the incoming flows + * including the current time step. + * Therefore the flows of the current time step should be calculated before using this function. + * @param[in] infectionState Specifies the #InfectionState we want to update. + * @param[in] IncomingFlows + * @param[in] NaturalDeathispossible Boolian that determines if there is the possibility of Natural Death in infectionState. + * @param[in] Transitionispossible Boolian that determines if it is possible to transition from the current InfectionState + * into another + * @param[in] dt + */ + void update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, bool NaturalDeathispossible, + bool Transitionispossible, ScalarType dt); + /** + * @brief Updates the values of one compartment, specified in infectionState, using the transitions. + * + * New value is stored in populations_update. The values are calculated using the compartment size in the previous + * time step and the related flows of the current time step. + * Therefore the flows of the current time step should be calculated before using this function. + * @param[in] infectionState Specifies the #InfectionState we want to update. + * @param[in] IncomingFlows + * @param[in] OutgoingFlows + * @param[in] NaturalDeathispossible Boolian that determines if there is the possibility of Natural Death in infectionState. + */ + void update_compartment_from_flow(InfectionState infectionState, + std::vector const& IncomingFlows, + std::vector const& OutgoingFlowsm, + bool NaturalDeathispossible, ScalarType dt); + + /** + * @brief Updates the values of all compartments except Susceptible + * + * New value is stored in populations. Values are computed via the function update_compartment_from_flow + */ + void update_compartments(ScalarType dt); + + /** + * @brief Compute the population size for the current last time in populations. + * + * The population size is computed as the sum of all compartments. + * The population size is stored in the vector m_populationsize. + */ + void compute_populationsize(); + + /** + * @brief Compute the normalized compartments for the current last time in m_normalizedpopulations. + * + * The normalized compartments are computed as populations / m_populationsize. + */ + void compute_normalizedcompartments(); + + /** + * @brief Computes the force of infection for the current last time in transitions. + * + * Computed value is stored in m_forceofinfection. + * + * @param[in] dt Time discretization step size. + */ + void compute_forceofinfection(ScalarType dt); + + // ---- Private parameters. ---- + + TimeSeries m_forceofinfection{ + TimeSeries(1)}; ///< TimeSeries containing the Force of infection term for every time point, + // needed for the numerical scheme. + TimeSeries m_forceofinfectionupdate{ + TimeSeries(1)}; ///< TimeSeries containing the Force of infection term for every time point, + // needed for the numerical scheme. For the numerical scheme using the update formula for the compartments. + TimeSeries m_totalpopulation{TimeSeries( + 1)}; ///< TimeSeries containing the total population size of the considered region for each time point. + TimeSeries m_totalpopulationupdate{TimeSeries( + 1)}; ///< TimeSeries containing the total population size of the considered region for each time point. + //In this case we use the compartments from populations_update. + TimeSeries m_totalpopulation_derivative{TimeSeries( + 1)}; ///< TimeSeries containing the derivative of the total population size of the considered + // region for each time point. + TimeSeries m_normalizedpopulations{ + TimeSeries(Eigen::Index(InfectionState::Count) - + 1)}; ///< TimeSeries containing points of time and the corresponding portion + // of people in defined #IndectionState%s. + + // ---- Friend classes/functions. ---- + // In the Simulation class, the actual simulation is performed which is why it needs access to the here + // defined (and private) functions to solve the model equations. + friend class Simulation; +}; + +} // namespace endisecir +} // namespace mio + +#endif //IDE_END_SECIR_MODEL_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_secir/normalized_model.cpp b/cpp/models/ide_endemic_secir/normalized_model.cpp new file mode 100644 index 0000000000..f570aff0c7 --- /dev/null +++ b/cpp/models/ide_endemic_secir/normalized_model.cpp @@ -0,0 +1,339 @@ +#include "ide_endemic_secir/normalized_model.h" +#include "ide_endemic_secir/computed_parameters.h" +#include "ide_endemic_secir/model.h" +#include "ide_endemic_secir/parameters.h" +#include "ide_endemic_secir/infection_state.h" +#include "memilio/config.h" +#include "memilio/utils/logging.h" + +#include "memilio/utils/time_series.h" +#include "vector" +#include +#include +#include +#include +#include +#include +#include + +namespace mio +{ +namespace endisecir +{ + +NormModel::NormModel(CompParameters const& compparams) + : compparameters{std::make_shared(compparams)} + , transitions{TimeSeries(Eigen::Index(InfectionTransition::Count))} + , populations{TimeSeries(Eigen::Index(InfectionState::Count) - 1)} +{ + //Set populations at start time t0. + TimeSeries::Vector vec_initnormalizedpopulations = + TimeSeries::Vector(Eigen::Index(InfectionState::Count) - 1); + for (int infection_state = 0; infection_state < Eigen::Index(InfectionState::Count) - 1; infection_state++) { + vec_initnormalizedpopulations[infection_state] = + compparameters->m_statesinit[0][infection_state] / compparameters->m_totalpopulationinit; + } + populations.add_time_point(0, vec_initnormalizedpopulations); + + // Set flows at start time t0. + // As we assume that all individuals have infectio age 0 at time t0, the flows at t0 are set to 0. + transitions.add_time_point( + 0, TimeSeries::Vector::Constant(static_cast(InfectionTransition::Count), 0.)); + + m_forceofinfection.add_time_point(0, TimeSeries::Vector::Constant(1, compparameters->m_NormFoI_0[0])); +} + +bool NormModel::check_constraints() const +{ + if (!(static_cast(populations.get_num_elements()) == static_cast(InfectionState::Count) - 1)) { + log_error(" A variable given for model construction is not valid. Number of elements in vector of populations " + "does not match the required number."); + return true; + } + + for (int i = 0; i < static_cast(InfectionState::Count) - 1; i++) { + if (populations[0][i] < 0) { + log_error("Initialization failed. Initial values for populations are less than zero."); + return true; + } + } + return compparameters->check_constraints(); +} +// ----Functionality for Initialization. ---- + +// ----Functionality for the iterations of a simulation. ---- +void NormModel::compute_susceptibles(ScalarType dt) +{ + Eigen::Index num_time_points = populations.get_num_time_points(); + populations.get_last_value()[static_cast(InfectionState::Susceptible)] = + (populations[num_time_points - 2][static_cast(InfectionState::Susceptible)] + + dt * compparameters->parameters.get()) / + (1 + dt * (m_forceofinfection[num_time_points - 2][0] + compparameters->parameters.get()) - + transitions[num_time_points - 2][(Eigen::Index)InfectionTransition::InfectedCriticalToDead]); +} + +void NormModel::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt) +{ + Eigen::Index calc_time_index = + (Eigen::Index)std::ceil(compparameters->m_transitiondistributions_support_max[idx_InfectionTransitions] / dt); + ScalarType current_time_age = static_cast(current_time_index) * dt; + ScalarType sum = 0; + //Determine the starting point of the for loop. + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + + for (Eigen::Index i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = static_cast(i) * dt; + + sum += (transitions[i + 1][idx_IncomingFlow] + + (compparameters->parameters.get() + + transitions[i][(Eigen::Index)InfectionTransition::InfectedCriticalToDead] - + compparameters->parameters.get()) * + populations[i][idx_CurrentCompartment]) * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)) * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index - i]; + } + if (current_time_index <= calc_time_index) { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum - + std::exp((-compparameters->parameters.get()) * (current_time_age)) * + populations[0][idx_CurrentCompartment] * + compparameters->parameters.get()[idx_InfectionTransitions] * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index]; + } + else { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum; + } +} + +void NormModel::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt) +{ + Eigen::Index current_time_index = transitions.get_num_time_points() - 1; + compute_flow(idx_InfectionTransitions, idx_IncomingFlow, idx_CurrentCompartment, current_time_index, dt); +} + +void NormModel::flows_currents_timestep(ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + // Calculate the transition SusceptibleToExposed with force of infection from previous time step and Susceptibles from + // current time step. + transitions.get_last_value()[static_cast(InfectionTransition::SusceptibleToExposed)] = + m_forceofinfection[current_time_index - 1][0] * + populations.get_last_value()[static_cast(InfectionState::Susceptible)]; + + // Calculate the other Transitions with compute_flow. + // Exposed to InfectedNoSymptoms + compute_flow(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), + Eigen::Index(InfectionTransition::SusceptibleToExposed), Eigen::Index(InfectionState::Exposed), dt); + // InfectedNoSymptoms to InfectedSymptoms + compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), + Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), + Eigen::Index(InfectionState::InfectedNoSymptoms), dt); + // InfectedNoSymptoms to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered), + Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), + Eigen::Index(InfectionState::InfectedNoSymptoms), dt); + // InfectedSymptoms to InfectedSevere + compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), + Eigen::Index(InfectionState::InfectedSymptoms), dt); + // InfectedSymptoms to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered), + Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), + Eigen::Index(InfectionState::InfectedSymptoms), dt); + // InfectedSevere to InfectedCritical + compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + Eigen::Index(InfectionState::InfectedSevere), dt); + // InfectedCritical to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered), + Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), + Eigen::Index(InfectionState::InfectedSevere), dt); + // InfectedCritical to Dead + compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead), + Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + Eigen::Index(InfectionState::InfectedCritical), dt); + // InfectedCritical to Recovered + compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered), + Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), + Eigen::Index(InfectionState::InfectedCritical), dt); +} + +void NormModel::update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, + bool Transitionispossible, ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + ScalarType current_time_age = (ScalarType)current_time_index * dt; + Eigen::Index calc_time_index = current_time_index; + if (Transitionispossible) { + calc_time_index = compparameters->m_transitiondistributions[(int)infectionState].size() - 1; + } + + ScalarType sum = 0; + + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + for (int i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = (ScalarType)i * dt; + ScalarType sum_inflows = 0; + for (const InfectionTransition& inflow : IncomingFlows) { + sum_inflows += transitions[i + 1][(int)inflow]; + } + if (Transitionispossible) { + sum += compparameters->m_transitiondistributions[(int)infectionState][current_time_index - i] * + (sum_inflows + (compparameters->parameters.get() + + transitions[i][(Eigen::Index)InfectionTransition::InfectedCriticalToDead] - + compparameters->parameters.get()) * + populations[i][(int)infectionState]) * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + else { + sum += (sum_inflows + (compparameters->parameters.get() + + transitions[i][(Eigen::Index)InfectionTransition::InfectedCriticalToDead] - + compparameters->parameters.get()) * + populations[i][(int)infectionState]) * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + } + if (Transitionispossible) { + if (current_time_index <= calc_time_index) { + populations.get_last_value()[(int)infectionState] = + dt * sum + compparameters->m_transitiondistributions[(int)infectionState][current_time_index] * + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } + else { + populations.get_last_value()[(int)infectionState] = dt * sum; + } + } + else { + populations.get_last_value()[(int)infectionState] = + dt * sum + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } +} + +void NormModel::update_compartments(ScalarType dt) +{ + + // Exposed + update_compartment_with_sum(InfectionState::Exposed, {InfectionTransition::SusceptibleToExposed}, true, dt); + // InfectedNoSymptoms + update_compartment_with_sum(InfectionState::InfectedNoSymptoms, {InfectionTransition::ExposedToInfectedNoSymptoms}, + true, dt); + + // InfectedSymptoms + update_compartment_with_sum(InfectionState::InfectedSymptoms, + {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms}, true, dt); + + // InfectedSevere + update_compartment_with_sum(InfectionState::InfectedSevere, {InfectionTransition::InfectedSymptomsToInfectedSevere}, + true, dt); + + // InfectedCritical + update_compartment_with_sum(InfectionState::InfectedCritical, + {InfectionTransition::InfectedSevereToInfectedCritical}, true, dt); + // Recovered + update_compartment_with_sum(InfectionState::Recovered, + { + InfectionTransition::InfectedNoSymptomsToRecovered, + InfectionTransition::InfectedSymptomsToRecovered, + InfectionTransition::InfectedSevereToRecovered, + InfectionTransition::InfectedCriticalToRecovered, + }, + false, dt); +} + +void NormModel::compute_forceofinfection(ScalarType dt) +{ + + Eigen::Index num_time_points = populations.get_num_time_points(); + ScalarType current_time = populations.get_last_time(); + + // Determine the starting point of the for loop. + Eigen::Index starting_point1 = std::max(0, (int)num_time_points - 1 - (int)compparameters->m_infectivity.size()); + + ScalarType sum1 = 0.; + // Compute the first sum in the force of infection term. + for (Eigen::Index i = starting_point1; i < num_time_points - 1; i++) { + sum1 += compparameters->m_infectivity[num_time_points - 1 - i] * + (populations[i + 1][(int)InfectionState::Susceptible] * m_forceofinfection[i][0] + + (compparameters->parameters.get() + + transitions[i + 1][(int)InfectionTransition::InfectedCriticalToDead] - + compparameters->parameters.get()) * + populations[i + 1][(int)InfectionState::Exposed]); + } + + // Determine the starting point of the for loop. + Eigen::Index starting_point2 = std::max(0, (int)num_time_points - 1 - (int)compparameters->m_B.size()); + + ScalarType sum2 = 0.; + // Compute the second sum in the force of infection term. + for (Eigen::Index i = starting_point2; i < num_time_points - 1; i++) { + ScalarType state_age = static_cast(i) * dt; + sum2 -= (compparameters->parameters.get() + + transitions[i + 1][(int)InfectionTransition::InfectedCriticalToDead] - + compparameters->parameters.get()) * + std::exp(-compparameters->parameters.get() * (current_time - state_age)) * + populations[i + 1][(int)InfectionState::InfectedNoSymptoms] * + compparameters->m_B[num_time_points - 1 - i]; + } + + // Determine the starting point of the for loop. + Eigen::Index starting_point3 = + std::max(0, (int)num_time_points - 1 - + (int)compparameters->m_transitiondistributions[(int)InfectionState::InfectedNoSymptoms].size()); + + ScalarType sum3 = 0.; + // Compute the third sum in the force of infection term. + for (Eigen::Index i = starting_point3; i < num_time_points - 1; i++) { + ScalarType state_age = static_cast(i) * dt; + sum3 += + (compparameters->parameters.get() + + transitions[i + 1][(int)InfectionTransition::InfectedCriticalToDead] - + compparameters->parameters.get()) * + std::exp(-compparameters->parameters.get() * (current_time - state_age)) * + populations[i + 1][(int)InfectionState::InfectedNoSymptoms] * + compparameters->parameters.get().eval(current_time - state_age) * + compparameters->m_transitiondistributions[(int)InfectionState::InfectedNoSymptoms][num_time_points - 1 - i]; + } + + // Determine the starting point of the for loop. + Eigen::Index starting_point4 = + std::max(0, (int)num_time_points - 1 - + (int)compparameters->m_transitiondistributions[(int)InfectionState::InfectedSymptoms].size()); + + ScalarType sum4 = 0.; + // Compute the third sum in the force of infection term. + for (Eigen::Index i = starting_point4; i < num_time_points - 1; i++) { + ScalarType state_age = static_cast(i) * dt; + sum4 += + (compparameters->parameters.get() + + transitions[i + 1][(int)InfectionTransition::InfectedCriticalToDead] - + compparameters->parameters.get()) * + std::exp(-compparameters->parameters.get() * (current_time - state_age)) * + populations[i + 1][(int)InfectionState::InfectedSymptoms] * + compparameters->parameters.get().eval(current_time - state_age) * + compparameters->m_transitiondistributions[(int)InfectionState::InfectedSymptoms][num_time_points - 1 - i]; + } + ScalarType sum = sum1 + sum2 + sum3 + sum4; + + m_forceofinfection.get_last_value()[0] = + dt * sum * + (compparameters->parameters.get().eval(current_time) * + compparameters->parameters.get().get_cont_freq_mat().get_matrix_at(current_time)(0, 0)); + + // Add inital functions for the force of infection in case they still have an influence. + if (num_time_points <= (int)compparameters->m_NormFoI_0.size()) { + m_forceofinfection.get_last_value()[0] += compparameters->m_NormFoI_0[num_time_points - 1]; + } + if (num_time_points <= (int)compparameters->m_NormInitFoI.size()) { + m_forceofinfection.get_last_value()[0] += compparameters->m_NormInitFoI[num_time_points - 1]; + } +} + +} // namespace endisecir + +} // namespace mio \ No newline at end of file diff --git a/cpp/models/ide_endemic_secir/normalized_model.h b/cpp/models/ide_endemic_secir/normalized_model.h new file mode 100644 index 0000000000..a2809f12e2 --- /dev/null +++ b/cpp/models/ide_endemic_secir/normalized_model.h @@ -0,0 +1,155 @@ +#ifndef IDE_END_SECIR_NORMMODEL_H +#define IDE_END_SECIR_NORMMODEL_H + +#include "ide_endemic_secir/infection_state.h" +#include "ide_endemic_secir/parameters.h" +#include "ide_endemic_secir/computed_parameters.h" +#include "ide_endemic_secir/model.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" + +#include "vector" +#include +#include +#include + +namespace mio +{ +namespace endisecir +{ + +// Forward declaration of friend classes/functions of Model. +class NormModel; +class Simulation; + +class NormModel +{ + using ParameterSet = Parameters; + +public: + /** + * @brief Constructor to create an endemic IDE_SECIR model. + * + * @param[in] TODO!!!!!! + */ + + NormModel(CompParameters const& compparams); + + /** + * @brief Checks constraints on model parameters and initial data. + * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. + */ + bool check_constraints() const; + + // ---- Public parameters. ---- + std::shared_ptr compparameters; + TimeSeries + transitions; ///< TimesSeries containing points of time and the corresponding number of individuals transitioning + // from one #InfectionState to another as defined in #Infection%s. + TimeSeries populations; ///< TimeSeries containing points of time and the corresponding number of people + // in defined #InfectionState%s. + +private: + // ---- Functionality for the iterations of a simulation. + + /** + * @brief Computes number of Susceptibles for the current last time in populations. + * + * Number is computed by using the previous number of Susceptibles, total population of the last time point and + * the force of infection (also from the last time point). + * Number is stored at the matching index in populations. + * + * @param[in] dt Time discretization step size. + */ + void compute_susceptibles(ScalarType dt); + + /** + * @brief Computes sizeof a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the time index + * current_time_index. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] current_time_index The time index the transition should be computed for. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt); + + /** + * @brief Computes size of a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the current + * last time value in transitions. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt); + + /** + * @brief Sets all required transitions for the current last timestep in transitions. + * + * New values are stored in transitions. Most values are computed via the function compute_flow() + * + * @param[in] dt time discretization step size. + */ + void flows_currents_timestep(ScalarType dt); + + /** + * @brief Updates the values of one compartment, specified in infectionState, using all past transitions. + * + * New value is stored in populations. The values are calculated using all past values for the incoming flows + * including the current time step. + * Therefore the flows of the current time step should be calculated before using this function. + * @param[in] infectionState Specifies the #InfectionState we want to update. + * @param[in] IncomingFlows + * @param[in] NaturalDeathispossible Boolian that determines if there is the possibility of Natural Death in infectionState. + * @param[in] Transitionispossible Boolian that determines if it is possible to transition from the current InfectionState + * into another + * @param[in] dt + */ + void update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, bool Transitionispossible, + ScalarType dt); + + /** + * @brief Updates the values of all compartments except Susceptible + * + * New value is stored in populations. Values are computed via the function update_compartment_from_flow + */ + void update_compartments(ScalarType dt); + + /** + * @brief Computes the force of infection for the current last time in transitions. + * + * Computed value is stored in m_forceofinfection. + * + * @param[in] dt Time discretization step size. + */ + void compute_forceofinfection(ScalarType dt); + + // ---- Private parameters. ---- + + TimeSeries m_forceofinfection{ + TimeSeries(1)}; ///< TimeSeries containing the Force of infection term for every time point, + + // ---- Friend classes/functions. ---- + // In the Simulation class, the actual simulation is performed which is why it needs access to the here + // defined (and private) functions to solve the model equations. + friend class Simulation; +}; + +} // namespace endisecir +} // namespace mio + +#endif // IDE_END_SECIR_NORMMODEL_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_secir/parameters.h b/cpp/models/ide_endemic_secir/parameters.h new file mode 100644 index 0000000000..e09cf6e8f9 --- /dev/null +++ b/cpp/models/ide_endemic_secir/parameters.h @@ -0,0 +1,314 @@ +#ifndef IDE_END_SECIR_PARAMS_H +#define IDE_END_SECIR_PARAMS_H + +#include "memilio/config.h" +#include "memilio/epidemiology/contact_matrix.h" +#include "memilio/math/floating_point.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/parameter_set.h" +#include "ide_endemic_secir/infection_state.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/epidemiology/uncertain_matrix.h" + +#include +#include +#include + +namespace mio +{ +namespace endisecir +{ + +/************************************************** +* Define Parameters of the IDE-END-SECIHURD model * +**************************************************/ + +/** + * @brief Transitions distribution for each transition in #Infection Transition. + * + * we use as a default SmootherCosine functions for all transitions with m_paramter=2. + */ + +struct TransitionDistributions { + using Type = std::vector; + static Type get_default() + { + SmootherCosine smoothcos(2.0); + StateAgeFunctionWrapper delaydistribution(smoothcos); + return Type((int)InfectionTransition::Count, delaydistribution); + } + + static std::string mame() + { + return "TransitionDistributions"; + } +}; + +/** + * @brief Defines the probability for each possible transitions to take this transition- + */ + +struct TransitionProbabilities { + /*For consistency, also define TransitionProbabilities for each transition in #InfectionTransition. + Transition Probabilities should be set to 1 if there is no possible other flow from starting compartment.*/ + using Type = std::vector; + static Type get_default() + { + std::vector probs((int)InfectionTransition::Count, 0.5); + // Set the following probablities to 1 as there is no other option to go anywhere else. + probs[Eigen::Index(InfectionTransition::SusceptibleToExposed)] = 1; + probs[Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)] = 1; + return Type(probs); + } + static std::string name() + { + return "TransitionsProbabilities"; + } +}; + +/** + * @brief The contact patterns within the society are modelled using an UncertainContactMatrix. + */ +struct ContactPatterns { + using Type = UncertainContactMatrix; + static Type get_default() + { + ContactMatrixGroup contact_matrix = ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.)); + return Type(contact_matrix); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/** +* @brief Probability of getting infected from a contact. +*/ +struct TransmissionProbabilityOnContact { + using Type = StateAgeFunctionWrapper; + static Type get_default() + { + ConstantFunction constfunc(1.0); + return Type(constfunc); + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/** +* @brief The relative InfectedNoSymptoms infectability. +*/ +struct RelativeTransmissionNoSymptoms { + using Type = StateAgeFunctionWrapper; + static Type get_default() + { + ConstantFunction constfunc(1.0); + return Type(constfunc); + } + static std::string name() + { + return "RelativeTrabsnissionNoSymptoms"; + } +}; + +/** +* @brief The risk of infection from symptomatic cases in the SECIR model. +*/ +struct RiskOfInfectionFromSymptomatic { + using Type = StateAgeFunctionWrapper; + static Type get_default() + { + ConstantFunction constfunc(1.0); + return Type(constfunc); + } + static std::string name() + { + return "RiskOfInfectionFromSymptomatic"; + } +}; + +/** + * @brief The natural birth rate. + */ +struct NaturalBirthRate { + using Type = ScalarType; + static Type get_default() + { + return Type(5e-5); + } + static std::string name() + { + return "NaturalBirthRate"; + } +}; + +/** + * @brief The natural death rate. + */ +struct NaturalDeathRate { + using Type = ScalarType; + static Type get_default() + { + return Type(3e-5); + } + static std::string name() + { + return "NaturalDeathRate"; + } +}; + +// Define Parameterset for IDE-END-SECIR model. +using ParametersBase = + ParameterSet; + +/** + * @brief Parameters of an endemic SECIR/SECIHURD model. + */ + +class Parameters : public ParametersBase +{ +public: + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and logs an error. + */ + bool check_constraints() const + { + // For paramerers potentially depending on the infectious age, the values are checked equidistantly on a + // realistic maximum window. + size_t infectious_window_check = 50; // parameter defining minmal window on x-axis. + + //Check if all the probabilitys are between 0 and 1. + for (size_t i = 0; i < infectious_window_check; i++) { + if (this->get().eval((static_cast(i))) < 0.0 || + this->get().eval(static_cast(i)) > 1.0) { + log_error( + "Constraint check: TransmissionProbabilityOnContact smaller {:d} or larger {:d} at some time {:d}", + 0, 1, i); + return true; + } + } + + for (size_t i = 0; i < infectious_window_check; i++) { + if (this->get().eval((static_cast(i))) < 0.0 || + this->get().eval(static_cast(i)) > 1.0) { + log_error( + "Constraint check: RelativeTransmissionNoSymptoms smaller {:d} or larger {:d} at some time {:d}", 0, + 1, i); + return true; + } + } + + for (size_t i = 0; i < infectious_window_check; i++) { + if (this->get().eval((static_cast(i))) < 0.0 || + this->get().eval(static_cast(i)) > 1.0) { + log_error( + "Constraint check: RiskOfInfectionFromSymptomatic smaller {:d} or larger {:d} at some time {:d}", 0, + 1, i); + return true; + } + } + + for (size_t i = 0; i < static_cast(InfectionTransition::Count); i++) { + if (this->get()[i] < 0.0 || this->get()[i] > 1.0) { + log_error("Constraint check: One parameter in TransitionProbabilities smaller {:d} or larger {:d}", 0, + 1); + return true; + } + } + + // The TransitionProbabilities SusceptibleToExposed and ExposedToInfectedNoSymptoms should be 1. + if (!floating_point_equal( + this->get()[static_cast(InfectionTransition::SusceptibleToExposed)], 1.0, + 1e-14)) { + log_error("Constraint check: Parameter transition probability for SusceptibleToExposed unequal to {:d}", 1); + return true; + } + + if (!floating_point_equal(this->get()[static_cast( + InfectionTransition::ExposedToInfectedNoSymptoms)], + 1.0, 1e-14)) { + log_error("Constraint check: Parameter transition probability for ExposedToInfectedNoSymptoms unequal " + "to {:d}", + 1); + return true; + } + + // All TransitionProbabilities where the Transition starts in the same compartment should sum up to 1. + + if (!floating_point_equal(this->get()[static_cast( + InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] + + this->get()[static_cast( + InfectionTransition::InfectedNoSymptomsToRecovered)], + 1.0, 1e-14)) { + log_error("Constraint check: Sum of transition probability for InfectedNoSymptomsToInfectedSymptoms and " + "InfectedNoSymptomsToRecovered not equal to {:d}", + 1); + return true; + } + + if (!floating_point_equal(this->get()[static_cast( + InfectionTransition::InfectedSymptomsToInfectedSevere)] + + this->get()[static_cast( + InfectionTransition::InfectedSymptomsToRecovered)], + 1.0, 1e-14)) { + log_error("Constraint check: Sum of transition probability for InfectedSymptomsToInfectedSevere and " + "InfectedSymptomsToRecovered not equal to {:d}", + 1); + return true; + } + + if (!floating_point_equal(this->get()[static_cast( + InfectionTransition::InfectedSevereToInfectedCritical)] + + this->get()[static_cast( + InfectionTransition::InfectedSevereToRecovered)], + 1.0, 1e-14)) { + log_error("Constraint check: Sum of transition probability for InfectedSevereToInfectedCritical and " + "InfectedSevereToRecovered not equal to {:d}", + 1); + return true; + } + + if (!floating_point_equal( + this->get()[static_cast(InfectionTransition::InfectedCriticalToDead)] + + this->get()[static_cast( + InfectionTransition::InfectedCriticalToRecovered)], + 1.0, 1e-14)) { + log_error("Constraint check: Sum of transition probability for InfectedCriticalToDead and " + "InfectedCriticalToRecovered not equal to {:d}", + 1); + return true; + } + + /* The first entry of TransitionDistributions is not checked because the distribution S->E is never used + (and it makes no sense to use the distribution). The support does not need to be valid.*/ + for (size_t i = 1; i < static_cast(InfectionTransition::Count); i++) { + if (floating_point_less(this->get()[i].get_support_max(10), 0.0, 1e-14)) { + log_error("Constraint check: One function in TransitionDistributions has invalid support."); + return true; + } + } + + if (this->get() < 0.0) { + log_warning("Constraint check: Parameter NaturalBirthRate should be greater than {:d}", 0); + return true; + } + + if (this->get() < 0.0) { + log_warning("Constraint check: Parameter NaturalDeathRate should be greater than {:d}", 0); + return true; + } + + return false; + } +}; + +} // namespace endisecir + +} // namespace mio + +#endif //IDE_END_SECIR_PARAMS_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_secir/simulation.cpp b/cpp/models/ide_endemic_secir/simulation.cpp new file mode 100644 index 0000000000..1c41aad2d9 --- /dev/null +++ b/cpp/models/ide_endemic_secir/simulation.cpp @@ -0,0 +1,117 @@ + +#include "ide_endemic_secir/simulation.h" +#include "ide_endemic_secir/computed_parameters.h" +#include "ide_endemic_secir/model.h" +#include "ide_endemic_secir/normalized_model.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include + +namespace mio +{ +namespace endisecir +{ + +void Simulation::advance(ScalarType tmax) +{ + mio::log_info("Simulating IDE-END-SECIR from t0 = {} until tmax = {} with dt = {}.", + m_model->transitions.get_last_time(), tmax, m_dt); + + m_model->compparameters->set_transitiondistributions_support_max(m_dt); + m_model->compparameters->set_transitiondistributions(m_dt); + m_model->compparameters->set_transitiondistributions_derivative(m_dt); + m_model->compparameters->set_B(m_dt); + m_model->compparameters->set_infectivity(m_dt); + m_model->compparameters->set_FoI_0(m_dt); + m_model->compparameters->set_InitFoI(m_dt); + m_model->compparameters->set_reproductionnumber_c(m_dt); + m_model->compparameters->set_T(m_dt); + m_model->compparameters->set_V(); + m_model->compparameters->set_W(m_dt); + + m_normmodel->compparameters->set_transitiondistributions_support_max(m_dt); + m_normmodel->compparameters->set_transitiondistributions(m_dt); + m_normmodel->compparameters->set_transitiondistributions_derivative(m_dt); + m_normmodel->compparameters->set_B(m_dt); + m_normmodel->compparameters->set_infectivity(m_dt); + m_normmodel->compparameters->set_FoI_0(m_dt); + m_normmodel->compparameters->set_InitFoI(m_dt); + m_normmodel->compparameters->set_reproductionnumber_c(m_dt); + m_normmodel->compparameters->set_T(m_dt); + m_normmodel->compparameters->set_V(); + m_normmodel->compparameters->set_W(m_dt); + + m_difference_normalizedcompartments.add_time_point(0); + m_difference_normalizedFoI.add_time_point(0); + compute_difference_normalizations(); + + // For every time step: + while (m_model->transitions.get_last_time() < tmax - m_dt / 2) { + + m_difference_normalizedcompartments.add_time_point(m_difference_normalizedcompartments.get_last_time() + m_dt); + m_difference_normalizedFoI.add_time_point(m_difference_normalizedFoI.get_last_time() + m_dt); + //standard model: + m_model->transitions.add_time_point(m_model->transitions.get_last_time() + m_dt); + m_model->transitions_update.add_time_point(m_model->transitions_update.get_last_time() + m_dt); + m_model->populations.add_time_point(m_model->populations.get_last_time() + m_dt); + m_model->populations_update.add_time_point(m_model->populations_update.get_last_time() + m_dt); + m_model->m_forceofinfection.add_time_point(m_model->m_forceofinfection.get_last_time() + m_dt); + m_model->m_forceofinfectionupdate.add_time_point(m_model->m_forceofinfectionupdate.get_last_time() + m_dt); + m_model->m_totalpopulation.add_time_point(m_model->m_totalpopulation.get_last_time() + m_dt); + m_model->m_totalpopulationupdate.add_time_point(m_model->m_totalpopulationupdate.get_last_time() + m_dt); + m_model->m_totalpopulation_derivative.add_time_point(m_model->m_totalpopulation_derivative.get_last_time() + + m_dt); + m_model->m_normalizedpopulations.add_time_point(m_model->m_normalizedpopulations.get_last_time() + m_dt); + + // Compute susceptibles: + m_model->compute_susceptibles(m_dt); + + // Compute flows: + m_model->flows_currents_timestep(m_dt); + + // Compute remaining compartments: + m_model->update_compartments(m_dt); + + // Compute m_populationsize: + m_model->compute_populationsize(); + + // Compute normalized compartments: + m_model->compute_normalizedcompartments(); + + // Compute m_forceofinfection; + m_model->compute_forceofinfection(m_dt); + + // normalized model: + m_normmodel->transitions.add_time_point(m_normmodel->transitions.get_last_time() + m_dt); + m_normmodel->populations.add_time_point(m_normmodel->populations.get_last_time() + m_dt); + m_normmodel->m_forceofinfection.add_time_point(m_normmodel->m_forceofinfection.get_last_time() + m_dt); + + // Compute susceptibles: + m_normmodel->compute_susceptibles(m_dt); + + // Compute flows: + m_normmodel->flows_currents_timestep(m_dt); + + // Compute remaining compartments: + m_normmodel->update_compartments(m_dt); + + // Compute m_forceofinfection; + m_normmodel->compute_forceofinfection(m_dt); + + //Compute the difference of the two normalized compartments. + compute_difference_normalizations(); + } +} + +TimeSeries simulate(ScalarType tmax, ScalarType dt, CompParameters const& m_compparams, + Model const& m_model, NormModel const& m_normmodel) +{ + m_model.check_constraints(); + Simulation sim(m_compparams, m_model, m_normmodel, dt); + sim.advance(tmax); + return sim.get_compartments(); +} + +} // namespace endisecir + +} // namespace mio diff --git a/cpp/models/ide_endemic_secir/simulation.h b/cpp/models/ide_endemic_secir/simulation.h new file mode 100644 index 0000000000..e1bb5640ed --- /dev/null +++ b/cpp/models/ide_endemic_secir/simulation.h @@ -0,0 +1,287 @@ +#ifndef IDE_END_SECIR_SIMULATION_H +#define IDE_END_SECIR_SIMULATION_H + +#include "ide_endemic_secir/computed_parameters.h" +#include "ide_endemic_secir/model.h" +#include "ide_endemic_secir/normalized_model.h" +#include "memilio/config.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/time_series.h" +#include +#include + +namespace mio +{ +namespace endisecir +{ + +class Simulation +{ +public: + Simulation(CompParameters const& compparams, Model const& model, NormModel const& normmodel, ScalarType dt = 0.1) + : m_compparams(std::make_unique(compparams)) + , m_model(std::make_unique(model)) + , m_normmodel(std::make_unique(normmodel)) + , m_dt(dt) + { + } + + /** + * @brief Run the simulation until time tmax. + */ + void advance(ScalarType tmax); + + /** + * @brief Get the result of the simulation for the compartments of m_model + * Return the number of persons in all #InfectionState%s + */ + TimeSeries get_compartments() + { + return m_model->populations; + } + /** + * @brief Get the result of the simulation for the compartments of m_model, where we use the update scheme. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries get_compartments_update() + { + return m_model->populations_update; + } + + /** + * @brief Get the result of the simulation for the normalized compartments, where we use m_model and the compartments + * computed using the sum scheme. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries get_normalizedcompartments() + { + return m_model->m_normalizedpopulations; + } + /** + * @brief Get the result of the simulation for the compartments of the normalized model. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries get_normmodel_compartments() + { + return m_normmodel->populations; + } + + /** + * @brief Get the result of the simulation for the compartments. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries& get_compartments() const + { + return m_model->populations; + } + + /** + * @brief Get the result of the simulation for the compartments of m_model, where we use the update scheme. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries& get_compartments_update() const + { + return m_model->populations_update; + } + + /** + * @brief Get the result of the simulation for the normalized compartments, where we use m_model. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries& get_normalizedcompartments() const + { + return m_model->m_normalizedpopulations; + } + + /** + * @brief Get the result of the simulation for the compartments of the normalized model. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries& get_normmodel_compartments() const + { + return m_normmodel->populations; + } + + /** + * @brief Get the transitions between the different #InfectionState%s for m_model. + */ + TimeSeries const& get_transitions() + { + return m_model->transitions; + } + /** + * @brief Get the transitions between the different #InfectionState%s for m_model using the update scheme. + */ + TimeSeries const& get_transitions_update() + { + return m_model->transitions_update; + } + + /** + * @brief Get the transitions between the different #InfectionState%s of m_normmodel. + */ + TimeSeries const& get_normmodel_transitions() + { + return m_normmodel->transitions; + } + + /** + * @brief Get the force of infection term of m_model. + */ + TimeSeries const& get_forceofinfections() + { + return m_model->m_forceofinfection; + } + + /** + * @brief Get the force of infection term of m_model using the update scheme. + */ + TimeSeries const& get_forceofinfections_update() + { + return m_model->m_forceofinfectionupdate; + } + + /** + * @brief Get the force of infection term of m_normmodel. + */ + TimeSeries const& get_normmodel_forceofinfections() + { + return m_normmodel->m_forceofinfection; + } + + /** + * @brief Get the total population of m_model. + */ + TimeSeries const& get_totalpopulations() + { + return m_model->m_totalpopulation; + } + + /** + * @brief Get the total population of m_model using the update scheme. + */ + TimeSeries const& get_totalpopulations_update() + { + return m_model->m_totalpopulationupdate; + } + + /** + * @brief Get the derivative of the total population of m_model. + */ + TimeSeries const& get_totalpopulations_derivative() + { + return m_model->m_totalpopulation_derivative; + } + + /** + * @brief Get the reproduction numer. + */ + ScalarType const& get_reproductionnumber_c() + { + return m_model->compparameters->m_reproductionnumber_c; + } + + /** + * @brief Get T. + */ + std::vector const& get_T() + { + return m_model->compparameters->m_T; + } + + /** + * @brief Get V. + */ + std::vector const& get_V() + { + return m_model->compparameters->m_V; + } + + /** + * @brief Get W. + */ + std::vector const& get_W() + { + return m_model->compparameters->m_W; + } + + /** + * @brief returns the simulation model used in simulation. + */ + const Model& get_model() const + { + return *m_model; + } + + /** + * @brief returns the simulation model used in simulation. + */ + Model& get_model() + { + return *m_model; + } + + /** + * @brief returns the simulation normmodel used in simulation. + */ + const NormModel& get_normmodel() const + { + return *m_normmodel; + } + + /** + * @brief returns the simulation normmodel used in simulation. + */ + NormModel& get_normmodel() + { + return *m_normmodel; + } + + /** + * @brief get the size of the time step of the simulation. + */ + ScalarType get_stepsize() + { + return m_dt; + } + + void compute_difference_normalizations() + { + for (int state = 0; state < (int)InfectionState::Count - 1; state++) { + m_difference_normalizedcompartments.get_last_value()[state] = + std::abs(m_normmodel->populations.get_last_value()[state] - + m_model->m_normalizedpopulations.get_last_value()[state]); + } + m_difference_normalizedFoI.get_last_value()[0] = std::abs(m_normmodel->m_forceofinfection.get_last_value()[0] - + m_model->m_forceofinfection.get_last_value()[0]); + } + + TimeSeries const& get_difference_normalizationcomp() + { + return m_difference_normalizedcompartments; + } + + TimeSeries const& get_difference_normalizationFoi() + { + return m_difference_normalizedFoI; + } + +private: + std::unique_ptr m_compparams; ///< Unique pointer to the computed Parameters. + std::unique_ptr m_model; ///< Unique pointer to the simulated Model. + std::unique_ptr m_normmodel; ///< Unique pointer to the simulated normalized Model. + ScalarType m_dt; ///< Time step used for numerical computations in simulation. + TimeSeries m_difference_normalizedcompartments{TimeSeries( + (int)InfectionState::Count - 1)}; ///< TimeSeries containing the difference of the compartments + // computed by NormModel and the normalized compartments computed in Model. + TimeSeries m_difference_normalizedFoI{ + TimeSeries(1)}; ///< TimeSeries containing the difference of the force of infection terms + // computed by NormModel and Model. +}; + +TimeSeries simulate(ScalarType tmax, ScalarType dt, Model const& model); + +} // namespace endisecir +} // namespace mio + +#endif //IDE_END_SECIR_SIMULATION_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_sird/CMakeLists.txt b/cpp/models/ide_endemic_sird/CMakeLists.txt new file mode 100644 index 0000000000..62d95776f5 --- /dev/null +++ b/cpp/models/ide_endemic_sird/CMakeLists.txt @@ -0,0 +1,17 @@ +add_library(ide_endemic_sird + infection_state.h + model.h + model.cpp + normalized_model.h + normalized_model.cpp + simulation.h + simulation.cpp + parameters.h + computed_parameters.h +) +target_link_libraries(ide_endemic_sird PUBLIC memilio) +target_include_directories(ide_endemic_sird PUBLIC + $ + $ +) +target_compile_options(ide_endemic_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ide_endemic_sird/computed_parameters.h b/cpp/models/ide_endemic_sird/computed_parameters.h new file mode 100644 index 0000000000..b1cd3950b6 --- /dev/null +++ b/cpp/models/ide_endemic_sird/computed_parameters.h @@ -0,0 +1,370 @@ +#ifndef IDE_END_SIRD_COMPPARAMS_H +#define IDE_END_SIRD_COMPPARAMS_H + +#include "ide_endemic_sird/infection_state.h" +#include "ide_endemic_sird/parameters.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" + +#include "vector" +#include +#include +#include +#include +#include + +namespace mio +{ +namespace endisird +{ +class CompParameters; +class Model; +class Simulation; + +class CompParameters +{ + using ParameterSet = Parameters; + +public: + CompParameters(TimeSeries&& states_init) + : parameters{Parameters()} + , m_statesinit{std::move(states_init)} + { + m_totalpopulationinit = std::accumulate(m_statesinit[0].begin(), m_statesinit[0].end(), 0) - + m_statesinit[0][(int)InfectionState::Dead]; + } + + bool check_constraints() const + { + if (!(static_cast(m_statesinit.get_num_elements()) == static_cast(InfectionState::Count))) { + log_error( + " A variable given for model construction is not valid. Number of elements in vector of populations " + "does not match the required number."); + return true; + } + for (int i = 0; i < static_cast(InfectionState::Count); i++) { + if (m_statesinit[0][i] < 0) { + log_error("Initialization failed. Initial values for populations are less than zero."); + return true; + } + } + return parameters.check_constraints(); + } + + /** + * @brief Setter for the tolerance used to calculate the maximum support of the TransitionDistributions. + * + * @param[in] new_tol New tolerance. + */ + void set_tol_for_support_max(ScalarType new_tol) + { + m_tol = new_tol; + } + + // ---- Public parameters. ---- + ParameterSet parameters; ///< ParameterSet of Model Parameters. + +private: + // ---- Functionality to set vectors with necessary information regarding TransitionDistributions. ---- + /** + * @brief Setter for the vector m_transitiondistributions_support_max that contains the support_max for all + * TransitionDistributions. + * + * This determines how many summands are required when calculating flows, the force of infection or compartments. + * + * @param[in] dt Time step size. + */ + void set_transitiondistributions_support_max(ScalarType dt) + { + // The transition Susceptible to Exposed is not needed in the computations. + for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) { + m_transitiondistributions_support_max[transition] = + parameters.get()[transition].get_support_max(dt, m_tol); + } + } + + /** + * @brief Setter for the vector m_transitiondistributions. + * + * Here we compute the weighted transition distributions for every compartment. Meaning for every compartment we weight + * the distributions of the outgoing transitions with their probability and compute the sum of the two weighed distributions. + * + * @param[in] dt Time step size. + */ + void set_transitiondistributions(ScalarType dt) + { + + // Vector containing the outgoing transitions for the Infected state. + + std::vector vector_transitions = {(int)InfectionTransition::InfectedToDead, + (int)InfectionTransition::InfectedToRecovered}; + + Eigen::Index calc_time_index = + std::ceil(std::max(m_transitiondistributions_support_max[(Eigen::Index)vector_transitions[0]], + m_transitiondistributions_support_max[(Eigen::Index)vector_transitions[1]]) / + dt); + + std::vector vec_contribution_to_foi_2(calc_time_index + 1, 0.); + for (Eigen::Index i = 0; i <= calc_time_index; i++) { + ///Compute the state_age. + ScalarType state_age = (ScalarType)i * dt; + + vec_contribution_to_foi_2[i] += + parameters.get()[vector_transitions[0]] * + parameters.get()[vector_transitions[0]].eval(state_age) + + parameters.get()[vector_transitions[1]] * + parameters.get()[vector_transitions[1]].eval(state_age); + m_transitiondistributions = vec_contribution_to_foi_2; + } + } + + /** + * @brief Setter for the vector m_transitiondistributions_derivative that contains the approximated derivative for + * all TransitionDistributions for all necessary time points. + * + * The derivative is approximated using a backwards difference scheme. + * The number of necessary time points for each TransitionDistribution is determined using + * m_transitiondistributions_support_max. + * + * @param[in] dt Time step size. + */ + void set_transitiondistributions_derivative(ScalarType dt) + { + // The transition SusceptibleToExposed 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 stores the value of the approximated derivative for all necessary time points. + std::vector vec_derivative(support_max_index + 1, 0.); + + for (Eigen::Index i = 0; i <= support_max_index; i++) { + // Compute state_age. + ScalarType state_age = (ScalarType)i * dt; + //Compute the apprximate derivative by a backwards difference scheme. + vec_derivative[i] = (parameters.get()[transition].eval(state_age) - + parameters.get()[transition].eval(state_age - dt)) / + dt; + } + m_transitiondistributions_derivative[transition] = vec_derivative; + } + } + + /** + *@brief Setter for the vector m_meaninfectivity that contains the approximated value of the mean infectivity for all + * for all necessary time points. + * + * This values is needed the compute the reproduction numer and the force of infection term. + * + * @param[in] dt Time step size. + */ + void set_infectivity(ScalarType dt) + { + // Compute the calc_time_indedx corresponding to a_1: + ScalarType calc_time_index = m_transitiondistributions.size(); + + m_infectivity = std::vector(calc_time_index, 0.); + + for (Eigen::Index time_point_index = 1; time_point_index < calc_time_index; time_point_index++) { + ScalarType time_point_age = (ScalarType)time_point_index * dt; + m_infectivity[time_point_index] = parameters.get().eval(time_point_age) * + m_transitiondistributions[time_point_index] * + std::exp(-parameters.get() * time_point_age); + } + } + + /** + *@brief Setter for the vectors m_FoI_0 and m_NormFoI_0 that contain the approximated values of the function FoI_0, + * that is a part of the force of infection term for all necessary time points. + */ + void set_FoI_0(ScalarType dt) + { + Eigen::Index calc_time_index = m_transitiondistributions.size(); + m_FoI_0 = std::vector(calc_time_index, 0.); + m_NormFoI_0 = std::vector(calc_time_index, 0.); + for (Eigen::Index time_point_index = 0; time_point_index < calc_time_index; time_point_index++) { + ScalarType time_point_age = (ScalarType)time_point_index * dt; + m_FoI_0[time_point_index] = + (parameters.get().eval(time_point_age) * + parameters.get().get_cont_freq_mat().get_matrix_at(time_point_age)(0, 0)) * + m_statesinit[0][(int)InfectionState::Infected] * m_infectivity[time_point_index]; + m_NormFoI_0[time_point_index] = + (parameters.get().eval(time_point_age) * + parameters.get().get_cont_freq_mat().get_matrix_at(time_point_age)(0, 0)) * + (m_statesinit[0][(int)InfectionState::Infected] / m_totalpopulationinit) * + m_infectivity[time_point_index]; + } + } + + // ---- Parameters needed for the analysis of the model ---- + /** + * @brief Setter for the Reproduction number m_reproductionnumber_c. + * + */ + void set_reproductionnumber_c(ScalarType dt) + { + // Determine the corresponding time index. + + Eigen::Index calc_time_index = m_infectivity.size() - 1; + ScalarType sum = 0; + for (int i = 0; i <= calc_time_index; i++) { + sum += m_infectivity[i]; + } + m_reproductionnumber_c = parameters.get().eval(0) * + parameters.get().get_cont_freq_mat().get_matrix_at(0)(0, 0) * dt * + sum; + } + + /** + * @brief Setter for m_T. + * + * @param[in] dt step size. + */ + void set_T(ScalarType dt) + { + // The value T_z1^z2 is not defined for the transition SusceptiblesToInfected. + for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) { + Eigen::Index support_max_index = + (Eigen::Index)std::ceil(m_transitiondistributions_support_max[transition] / dt); + ScalarType sum = 0; + + for (Eigen::Index i = 0; i <= support_max_index; i++) { + // Compute state_age. + ScalarType state_age = (ScalarType)i * dt; + //Compute the apprximate derivative by a backwards difference scheme. + sum -= std::exp(-parameters.get() * state_age) * + m_transitiondistributions_derivative[transition][i]; + } + m_T[transition] = parameters.get()[transition] * sum; + } + } + + /** + * @brief Setter for m_W. + * + * m_W contains the values W_i + * + * @param[in] dt step size. + */ + void set_W(ScalarType dt) + { + + Eigen::Index calc_time_index = m_transitiondistributions.size(); + ScalarType sum = 0; + + for (Eigen::Index i = 0; i < calc_time_index; i++) { + // Compute state_age. + ScalarType state_age = (ScalarType)i * dt; + //Compute the apprximate derivative by a backwards difference scheme. + sum += std::exp(-parameters.get() * state_age) * m_transitiondistributions[i]; + } + m_W = sum; + } + + void set_equilibria() + { + // We start by setting the equilibrium for the infected compartments. + ScalarType p = (m_W / m_T[(int)InfectionTransition::InfectedToDead]) * + (parameters.get() - 1 / m_W - parameters.get()) - + parameters.get() / + (m_reproductionnumber_c - m_T[(int)InfectionTransition::InfectedToDead]); + ScalarType q = parameters.get() * m_W / + ((m_reproductionnumber_c - m_T[(int)InfectionTransition::InfectedToDead]) * + m_T[(int)InfectionTransition::InfectedToDead]) * + ((1 / m_W) - parameters.get() + parameters.get() - + m_reproductionnumber_c); + m_compartments_equilibrium[0][(int)InfectionState::Infected] = -p / 2 + std::sqrt((p / 2) * (p / 2) - q); + m_compartments_equilibrium[1][(int)InfectionState::Infected] = -p / 2 - std::sqrt((p / 2) * (p / 2) - q); + + // From this we set the other equilibria points. + // Susceptibles: + m_compartments_equilibrium[0][(int)InfectionState::Susceptible] = + parameters.get() / + (m_compartments_equilibrium[0][(int)InfectionState::Infected] * + (m_reproductionnumber_c / m_W - m_T[(int)InfectionTransition::InfectedToDead] / m_W) + + parameters.get()); + m_compartments_equilibrium[1][(int)InfectionState::Susceptible] = + parameters.get() / + (m_compartments_equilibrium[1][(int)InfectionState::Infected] * + (m_reproductionnumber_c / m_W - m_T[(int)InfectionTransition::InfectedToDead] / m_W) + + parameters.get()); + // Force of Infection: + m_FoI_equilibrium[0] = + m_compartments_equilibrium[0][(int)InfectionState::Infected] * (m_reproductionnumber_c / m_W); + m_FoI_equilibrium[1] = + m_compartments_equilibrium[1][(int)InfectionState::Infected] * (m_reproductionnumber_c / m_W); + // InfectedToDead: + m_transitions_equilibrium[0][(int)InfectionTransition::InfectedToDead] = + m_compartments_equilibrium[0][(int)InfectionState::Infected] * + (m_T[(int)InfectionTransition::InfectedToDead] / m_W); + m_transitions_equilibrium[1][(int)InfectionTransition::InfectedToDead] = + m_compartments_equilibrium[1][(int)InfectionState::Infected] * + (m_T[(int)InfectionTransition::InfectedToDead] / m_W); + // InfectedToRecovered: + m_transitions_equilibrium[0][(int)InfectionTransition::InfectedToRecovered] = + (m_compartments_equilibrium[0][(int)InfectionState::Susceptible] * m_FoI_equilibrium[0] + + (parameters.get() + + m_transitions_equilibrium[0][(int)InfectionTransition::InfectedToDead] - + parameters.get()) * + m_compartments_equilibrium[0][(int)InfectionState::Infected]) * + m_T[(int)InfectionTransition::InfectedToDead]; + m_transitions_equilibrium[1][(int)InfectionTransition::InfectedToRecovered] = + (m_compartments_equilibrium[1][(int)InfectionState::Susceptible] * m_FoI_equilibrium[0] + + (parameters.get() + + m_transitions_equilibrium[1][(int)InfectionTransition::InfectedToDead] - + parameters.get()) * + m_compartments_equilibrium[1][(int)InfectionState::Infected]) * + m_T[(int)InfectionTransition::InfectedToDead]; + // Recovered: + m_compartments_equilibrium[0][(int)InfectionState::Recovered] = + m_transitions_equilibrium[0][(int)InfectionTransition::InfectedToRecovered] / + (parameters.get() - + m_transitions_equilibrium[0][(int)InfectionTransition::InfectedToDead]); + m_compartments_equilibrium[1][(int)InfectionState::Recovered] = + m_transitions_equilibrium[1][(int)InfectionTransition::InfectedToRecovered] / + (parameters.get() - + m_transitions_equilibrium[1][(int)InfectionTransition::InfectedToDead]); + } + + // ---- Private parameters. ---- + TimeSeries m_statesinit; ///< TimeSeries containing the initial values for the compartments. + ScalarType m_totalpopulationinit; + ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the TransitionDistributions. + std::vector m_transitiondistributions_support_max{ + std::vector((int)InfectionTransition::Count, 0.)}; ///< A vector containing the support_max + // for all TransitionDistributions. + std::vector m_transitiondistributions{ + std::vector(1, 0.)}; ///< A vector containing the weighted TransitionDistributions + // for the Infected state. + std::vector> m_transitiondistributions_derivative{std::vector>( + (int)InfectionTransition::Count, std::vector(1, 0.))}; ///< A Vector containing + // the approximated derivative for all TransitionDistributions for all necessary time points. + std::vector m_infectivity{ + std::vector(1, 0.)}; ///< A vector containing the approximated mean infectivity for all time points. + ScalarType m_reproductionnumber_c; ///< The control Reproduction number + std::vector m_FoI_0{std::vector(1, 0.)}; ///< A vector containing the approcimated + // value of the function FoI_0 used for the computation of the force of infection in the standard model + std::vector m_NormFoI_0{std::vector(1, 0.)}; ///< A vector containing the approcimated + // value of the function FoI_0 used for the computation of the force of infection in the normalized model + std::vector m_T{std::vector((int)InfectionTransition::Count, 0.)}; ///< A vector + // containing the approximated value for T_z1^z2 for every Flow z1 to z2. + ScalarType m_W{0}; ///< ScalarType of the value W_i. + std::vector> m_compartments_equilibrium{std::vector>( + 2, std::vector((int)InfectionState::Count - 1, 0.))}; ///< Vector containing the + // two computed equilibria points for the compartments of NormModel. + std::vector> m_transitions_equilibrium{std::vector>( + 2, std::vector((int)InfectionState::Count - 1, 0.))}; ///< Vector containing the + // two computed equilibria points for the transitions of NormModel. + std::vector m_FoI_equilibrium{std::vector(2, 0.)}; ///< A Vector containing the two + // computed equilibria points for the force of infection of NormModel. + // ---- Friend classes/functions. ---- + friend class Model; + friend class NormModel; + friend class Simulation; +}; + +} // namespace endisird + +} // namespace mio + +#endif //IDE_END_SIRD_COMPPARAMS_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_sird/infection_state.h b/cpp/models/ide_endemic_sird/infection_state.h new file mode 100644 index 0000000000..03b46edcc1 --- /dev/null +++ b/cpp/models/ide_endemic_sird/infection_state.h @@ -0,0 +1,38 @@ +#ifndef IDE_END_SIRD_INFECTIONSTATE_H +#define IDE_END_SIRD_INFECTIONSTATE_H + +namespace mio +{ + +namespace endisird +{ + +/** + * @brief The #InfectionState enum describes the possible + * categories for the infectious state of persons + */ +enum class InfectionState +{ + Susceptible = 0, + Infected = 1, + Recovered = 2, + Dead = 3, + Count = 4 +}; + +/** + * @brief The #InfectionTransition enum describes the possible + * transitions of the infectious state of persons. + */ +enum class InfectionTransition +{ + SusceptibleToInfected = 0, + InfectedToDead = 1, + InfectedToRecovered = 2, + Count = 3, +}; + +} // namespace endisird +} // namespace mio + +#endif //IDE_END_SIRD_INFECTIONSTATE_H diff --git a/cpp/models/ide_endemic_sird/model.cpp b/cpp/models/ide_endemic_sird/model.cpp new file mode 100644 index 0000000000..15b6456682 --- /dev/null +++ b/cpp/models/ide_endemic_sird/model.cpp @@ -0,0 +1,259 @@ +#include "ide_endemic_sird/model.h" +#include "ide_endemic_sird/computed_parameters.h" +#include "ide_endemic_sird/parameters.h" +#include "ide_endemic_sird/infection_state.h" +#include "memilio/config.h" +#include "memilio/utils/logging.h" + +#include "memilio/utils/time_series.h" +#include "vector" +#include +#include +#include +#include +#include +#include +#include + +namespace mio +{ +namespace endisird +{ +Model::Model(CompParameters const& compparams) + : compparameters{std::make_shared(compparams)} + , transitions{TimeSeries(Eigen::Index(InfectionTransition::Count))} + , populations{compparameters->m_statesinit} + +{ + + // Set flows at start time t0. + // As we assume that all individuals have infectio age 0 at time t0, the flows at t0 are set to 0. + transitions.add_time_point( + 0, TimeSeries::Vector::Constant(static_cast(InfectionTransition::Count), 0.)); + // Set population size at start timt t0. + ScalarType init_populationsize = + std::accumulate(populations[0].begin(), populations[0].end(), 0) - populations[0][(int)InfectionState::Dead]; + m_totalpopulation.add_time_point(0, TimeSeries::Vector::Constant(1, init_populationsize)); + m_totalpopulation_derivative.add_time_point(0, TimeSeries::Vector::Constant(1, 0.)); + + // Set the force of infection term at time t0. + m_forceofinfection.add_time_point(0, TimeSeries::Vector::Constant(1, compparameters->m_FoI_0[0])); + + //Set normalized_populations at start time t0. + TimeSeries::Vector vec_normalizedpopulations = + TimeSeries::Vector(Eigen::Index(InfectionState::Count) - 1); + for (int infection_state = 0; infection_state < Eigen::Index(InfectionState::Count) - 1; infection_state++) { + vec_normalizedpopulations[infection_state] = populations[0][infection_state] / m_totalpopulation[0][0]; + } + m_normalizedpopulations.add_time_point(0, vec_normalizedpopulations); +} + +bool Model::check_constraints() const +{ + if (!(static_cast(populations.get_num_elements()) == static_cast(InfectionState::Count))) { + log_error(" A variable given for model construction is not valid. Number of elements in vector of populations " + "does not match the required number."); + return true; + } + + for (int i = 0; i < static_cast(InfectionState::Count); i++) { + if (populations[0][i] < 0) { + log_error("Initialization failed. Initial values for populations are less than zero."); + return true; + } + } + return compparameters->check_constraints(); +} + +// ----Functionality for the iterations of a simulation. ---- +void Model::compute_susceptibles(ScalarType dt) +{ + Eigen::Index num_time_points = populations.get_num_time_points(); + populations.get_last_value()[static_cast(InfectionState::Susceptible)] = + (populations[num_time_points - 2][static_cast(InfectionState::Susceptible)] + + dt * m_totalpopulation[num_time_points - 2][0] * compparameters->parameters.get()) / + (1 + dt * (m_forceofinfection[num_time_points - 2][0] + compparameters->parameters.get())); +} + +void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt) +{ + Eigen::Index calc_time_index = + (Eigen::Index)std::ceil(compparameters->m_transitiondistributions_support_max[idx_InfectionTransitions] / dt); + ScalarType current_time_age = static_cast(current_time_index) * dt; + + ScalarType sum = 0; + //Determine the starting point of the for loop. + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + + for (Eigen::Index i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = static_cast(i) * dt; + sum += transitions[i + 1][idx_IncomingFlow] * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)) * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index - i]; + } + if (current_time_index <= calc_time_index) { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum - + std::exp((-compparameters->parameters.get()) * (current_time_age)) * + populations[0][idx_CurrentCompartment] * + compparameters->parameters.get()[idx_InfectionTransitions] * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index]; + } + else { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum; + } +} + +void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt) +{ + Eigen::Index current_time_index = transitions.get_num_time_points() - 1; + compute_flow(idx_InfectionTransitions, idx_IncomingFlow, idx_CurrentCompartment, current_time_index, dt); +} + +void Model::flows_currents_timestep(ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + // Calculate the transition SusceptibleToExposed with force of infection from previous time step and Susceptibles from + // current time step. + transitions.get_last_value()[static_cast(InfectionTransition::SusceptibleToInfected)] = + m_forceofinfection[current_time_index - 1][0] * + populations.get_last_value()[static_cast(InfectionState::Susceptible)]; + + // Calculate the other Transitions with compute_flow. + // Infected To Dead: + compute_flow(Eigen::Index(InfectionTransition::InfectedToDead), + Eigen::Index(InfectionTransition::SusceptibleToInfected), Eigen::Index(InfectionState::Infected), dt); + // Infected To Recovered: + compute_flow(Eigen::Index(InfectionTransition::InfectedToRecovered), + Eigen::Index(InfectionTransition::SusceptibleToInfected), Eigen::Index(InfectionState::Infected), dt); +} + +void Model::update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, + bool NaturalDeathispossible, bool Transitionispossible, ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + ScalarType current_time_age = (ScalarType)current_time_index * dt; + Eigen::Index calc_time_index = current_time_index; + if (Transitionispossible) { + calc_time_index = compparameters->m_transitiondistributions.size() - 1; + } + + ScalarType sum = 0; + + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + for (int i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = (ScalarType)i * dt; + ScalarType sum_inflows = 0; + for (const InfectionTransition& inflow : IncomingFlows) { + sum_inflows += transitions[i + 1][(int)inflow]; + } + if (NaturalDeathispossible && Transitionispossible) { + sum += compparameters->m_transitiondistributions[current_time_index - i] * sum_inflows * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + else if (NaturalDeathispossible && !Transitionispossible) { + sum += sum_inflows * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + // The case !NaturalDeathispossible && Transitionispossible is not possible, as if NaturalDeath is not possible + // this means you are in the Death compartment and then Transition is also not possible. + else { + sum += sum_inflows; + } + } + if (NaturalDeathispossible && Transitionispossible) { + if (current_time_index <= calc_time_index) { + populations.get_last_value()[(int)infectionState] = + dt * sum + compparameters->m_transitiondistributions[current_time_index] * + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } + else { + populations.get_last_value()[(int)infectionState] = dt * sum; + } + } + else if (NaturalDeathispossible && !Transitionispossible) { + populations.get_last_value()[(int)infectionState] = + dt * sum + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } + else { + populations.get_last_value()[(int)infectionState] = dt * sum + populations[0][(int)infectionState]; + } +} + +void Model::update_compartments(ScalarType dt) +{ + + // Infected + update_compartment_with_sum(InfectionState::Infected, {InfectionTransition::SusceptibleToInfected}, true, true, dt); + // Recovered + update_compartment_with_sum(InfectionState::Recovered, + { + InfectionTransition::InfectedToRecovered, + }, + true, false, dt); + + // Dead + update_compartment_with_sum(InfectionState::Dead, {InfectionTransition::InfectedToDead}, false, false, dt); +} + +void Model::compute_populationsize() +{ + ScalarType sum = 0; + for (int state = 0; state < Eigen::Index(InfectionState::Count) - 1; state++) { + sum += populations.get_last_value()[state]; + } + m_totalpopulation.get_last_value()[0] = sum; + // Here we comoute the derivative of the total population that is given by + // BirthRate * N(t) - DeathRate * N(t) - Transition[InfectedCritical>ToDeath](t). + m_totalpopulation_derivative.get_last_value()[0] = + (compparameters->parameters.get() - compparameters->parameters.get()) * + m_totalpopulation.get_last_value()[0] - + transitions.get_last_value()[(int)InfectionTransition::InfectedToDead]; +} + +void Model::compute_normalizedcompartments() +{ + for (int infection_state = 0; infection_state < Eigen::Index(InfectionState::Count) - 1; infection_state++) { + m_normalizedpopulations.get_last_value()[infection_state] = + populations.get_last_value()[infection_state] / m_totalpopulation.get_last_value()[0]; + } +} + +void Model::compute_forceofinfection(ScalarType dt) +{ + + Eigen::Index num_time_points = populations.get_num_time_points(); + ScalarType current_time = populations.get_last_time(); + + // Determine the starting point of the for loop. + Eigen::Index starting_point = std::max(0, (int)num_time_points - 1 - (int)compparameters->m_infectivity.size()); + + ScalarType sum = 0.; + // Compute the sum in the force of infection term. + for (Eigen::Index i = starting_point; i < num_time_points - 1; i++) { + sum += compparameters->m_infectivity[num_time_points - 1 - i] * + populations[i + 1][(int)InfectionState::Susceptible] * m_forceofinfection[i][0]; + } + + ScalarType Foi_0 = 0; + + // Add inital functions for the force of infection in case they still have an influence. + if (num_time_points <= (int)compparameters->m_FoI_0.size()) { + Foi_0 = compparameters->m_FoI_0[num_time_points - 1] / m_totalpopulation[num_time_points - 1][0]; + } + m_forceofinfection.get_last_value()[0] = + (dt * sum * compparameters->parameters.get().eval(current_time) * + compparameters->parameters.get().get_cont_freq_mat().get_matrix_at(current_time)(0, 0)) / + m_totalpopulation[num_time_points - 1][0] + + Foi_0; +} + +}; // namespace endisird + +} // namespace mio diff --git a/cpp/models/ide_endemic_sird/model.h b/cpp/models/ide_endemic_sird/model.h new file mode 100644 index 0000000000..017e98527e --- /dev/null +++ b/cpp/models/ide_endemic_sird/model.h @@ -0,0 +1,185 @@ +#ifndef IDE_END_SIRD_MODEL_H +#define IDE_END_SIRD_MODEL_H + +#include "ide_endemic_sird/infection_state.h" +#include "ide_endemic_sird/parameters.h" +#include "ide_endemic_sird/computed_parameters.h" +#include "memilio/config.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/time_series.h" + +#include "vector" +#include +#include +#include + +namespace mio +{ +namespace endisird +{ +// Forward declaration of friend classes/functions of Model. +class Model; +class Simulation; + +class Model +{ + using ParameterSet = Parameters; + +public: + /** + * @brief Constructor to create an endemic IDE_SIRD model. + * + * @param[in] TODO!!!!!! + */ + + Model(CompParameters const& compparams); + + /** + * @brief Checks constraints on model parameters and initial data. + * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. + */ + bool check_constraints() const; + + /** + * @brief Setter for the tolerance used to calculate the maximum support of the TransitionDistributions. + * + * @param[in] new_tol New tolerance. + */ + + // ---- Public parameters. ---- + std::shared_ptr compparameters; + TimeSeries + transitions; ///< TimesSeries containing points of time and the corresponding number of individuals transitioning + // from one #InfectionState to another as defined in #Infection%s. + TimeSeries populations; ///< TimeSeries containing points of time and the corresponding number of people + // in defined #InfectionState%s. In this case we compute them by a sum. + +private: + // ---- Functionality for the iterations of a simulation. + + /** + * @brief Computes number of Susceptibles for the current last time in populations. + * + * Number is computed by using the previous number of Susceptibles, total population of the last time point and + * the force of infection (also from the last time point). + * Number is stored at the matching index in populations and populations_update. + * + * @param[in] dt Time discretization step size. + */ + void compute_susceptibles(ScalarType dt); + + /** + * @brief Computes sizeof a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the time index + * current_time_index. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] current_time_index The time index the transition should be computed for. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt); + + /** + * @brief Computes size of a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the current + * last time value in transitions. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt); + + /** + * @brief Sets all required transitions for the current last timestep in transitions. + * + * New values are stored in transitions. Most values are computed via the function compute_flow() + * + * @param[in] dt time discretization step size. + */ + void flows_currents_timestep(ScalarType dt); + + /** + * @brief Updates the values of one compartment, specified in infectionState, using all past transitions. + * + * New value is stored in populations. The values are calculated using all past values for the incoming flows + * including the current time step. + * Therefore the flows of the current time step should be calculated before using this function. + * @param[in] infectionState Specifies the #InfectionState we want to update. + * @param[in] IncomingFlows + * @param[in] NaturalDeathispossible Boolian that determines if there is the possibility of Natural Death in infectionState. + * @param[in] Transitionispossible Boolian that determines if it is possible to transition from the current InfectionState + * into another + * @param[in] dt + */ + void update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, bool NaturalDeathispossible, + bool Transitionispossible, ScalarType dt); + + /** + * @brief Updates the values of all compartments except Susceptible + * + * New value is stored in populations. Values are computed via the function update_compartment_from_flow + */ + void update_compartments(ScalarType dt); + + /** + * @brief Compute the population size for the current last time in populations. + * + * The population size is computed as the sum of all compartments. + * The population size is stored in the vector m_populationsize. + */ + void compute_populationsize(); + + /** + * @brief Compute the normalized compartments for the current last time in m_normalizedpopulations. + * + * The normalized compartments are computed as populations / m_populationsize. + */ + void compute_normalizedcompartments(); + + /** + * @brief Computes the force of infection for the current last time in transitions. + * + * Computed value is stored in m_forceofinfection. + * + * @param[in] dt Time discretization step size. + */ + void compute_forceofinfection(ScalarType dt); + + // ---- Private parameters. ---- + + TimeSeries m_forceofinfection{ + TimeSeries(1)}; ///< TimeSeries containing the Force of infection term for every time point, + // needed for the numerical scheme. + TimeSeries m_totalpopulation{TimeSeries( + 1)}; ///< TimeSeries containing the total population size of the considered region for each time point. + TimeSeries m_totalpopulation_derivative{TimeSeries( + 1)}; ///< TimeSeries containing the derivative of the total population size of the considered + // region for each time point. + TimeSeries m_normalizedpopulations{ + TimeSeries(Eigen::Index(InfectionState::Count) - + 1)}; ///< TimeSeries containing points of time and the corresponding portion + // of people in defined #IndectionState%s. + + // ---- Friend classes/functions. ---- + // In the Simulation class, the actual simulation is performed which is why it needs access to the here + // defined (and private) functions to solve the model equations. + friend class Simulation; +}; + +} // namespace endisird +} // namespace mio + +#endif //IDE_END_SIRD_MODEL_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_sird/normalized_model.cpp b/cpp/models/ide_endemic_sird/normalized_model.cpp new file mode 100644 index 0000000000..e924d9e345 --- /dev/null +++ b/cpp/models/ide_endemic_sird/normalized_model.cpp @@ -0,0 +1,242 @@ +#include "ide_endemic_sird/normalized_model.h" +#include "ide_endemic_sird/computed_parameters.h" +#include "ide_endemic_sird/parameters.h" +#include "ide_endemic_sird/infection_state.h" +#include "memilio/config.h" +#include "memilio/utils/logging.h" + +#include "memilio/utils/time_series.h" +#include "vector" +#include +#include +#include +#include +#include +#include +#include + +namespace mio +{ +namespace endisird +{ + +NormModel::NormModel(CompParameters const& compparams) + : compparameters{std::make_shared(compparams)} + , transitions{TimeSeries(Eigen::Index(InfectionTransition::Count))} + , populations{TimeSeries(Eigen::Index(InfectionState::Count) - 1)} +{ + //Set populations at start time t0. + TimeSeries::Vector vec_initnormalizedpopulations = + TimeSeries::Vector(Eigen::Index(InfectionState::Count) - 1); + for (int infection_state = 0; infection_state < Eigen::Index(InfectionState::Count) - 1; infection_state++) { + vec_initnormalizedpopulations[infection_state] = + compparameters->m_statesinit[0][infection_state] / compparameters->m_totalpopulationinit; + } + populations.add_time_point(0, vec_initnormalizedpopulations); + + // Set flows at start time t0. + // As we assume that all individuals have infectio age 0 at time t0, the flows at t0 are set to 0. + transitions.add_time_point( + 0, TimeSeries::Vector::Constant(static_cast(InfectionTransition::Count), 0.)); + + m_forceofinfection.add_time_point(0, TimeSeries::Vector::Constant(1, compparameters->m_NormFoI_0[0])); + m_size.add_time_point(0, TimeSeries::Vector::Constant(1, 1.0)); +} + +bool NormModel::check_constraints() const +{ + if (!(static_cast(populations.get_num_elements()) == static_cast(InfectionState::Count) - 1)) { + log_error(" A variable given for model construction is not valid. Number of elements in vector of populations " + "does not match the required number."); + return true; + } + + for (int i = 0; i < static_cast(InfectionState::Count) - 1; i++) { + if (populations[0][i] < 0) { + log_error("Initialization failed. Initial values for populations are less than zero."); + return true; + } + } + return compparameters->check_constraints(); +} + +// ----Functionality for the iterations of a simulation. ---- +void NormModel::compute_susceptibles(ScalarType dt) +{ + Eigen::Index num_time_points = populations.get_num_time_points(); + populations.get_last_value()[static_cast(InfectionState::Susceptible)] = + (populations[num_time_points - 2][static_cast(InfectionState::Susceptible)] + + dt * compparameters->parameters.get()) / + (1 + dt * (m_forceofinfection[num_time_points - 2][0] + compparameters->parameters.get()) - + transitions[num_time_points - 2][(Eigen::Index)InfectionTransition::InfectedToDead]); +} + +void NormModel::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt) +{ + Eigen::Index calc_time_index = + (Eigen::Index)std::ceil(compparameters->m_transitiondistributions_support_max[idx_InfectionTransitions] / dt); + ScalarType current_time_age = static_cast(current_time_index) * dt; + ScalarType sum = 0; + //Determine the starting point of the for loop. + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + + for (Eigen::Index i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = static_cast(i) * dt; + + sum += + (transitions[i + 1][idx_IncomingFlow] + (compparameters->parameters.get() + + transitions[i][(Eigen::Index)InfectionTransition::InfectedToDead] - + compparameters->parameters.get()) * + populations[i][idx_CurrentCompartment]) * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)) * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index - i]; + } + if (current_time_index <= calc_time_index) { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum - + std::exp((-compparameters->parameters.get()) * (current_time_age)) * + populations[0][idx_CurrentCompartment] * + compparameters->parameters.get()[idx_InfectionTransitions] * + compparameters->m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index]; + } + else { + transitions.get_value(current_time_index)[idx_InfectionTransitions] = + (-dt) * compparameters->parameters.get()[idx_InfectionTransitions] * sum; + } +} + +void NormModel::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt) +{ + Eigen::Index current_time_index = transitions.get_num_time_points() - 1; + compute_flow(idx_InfectionTransitions, idx_IncomingFlow, idx_CurrentCompartment, current_time_index, dt); +} + +void NormModel::flows_currents_timestep(ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + // Calculate the transition SusceptibleToExposed with force of infection from previous time step and Susceptibles from + // current time step. + transitions.get_last_value()[static_cast(InfectionTransition::SusceptibleToInfected)] = + m_forceofinfection[current_time_index - 1][0] * + populations.get_last_value()[static_cast(InfectionState::Susceptible)]; + + // Calculate the other Transitions with compute_flow. + // Infected To Dead: + compute_flow(Eigen::Index(InfectionTransition::InfectedToDead), + Eigen::Index(InfectionTransition::SusceptibleToInfected), Eigen::Index(InfectionState::Infected), dt); + // Infected To Recovered: + compute_flow(Eigen::Index(InfectionTransition::InfectedToRecovered), + Eigen::Index(InfectionTransition::SusceptibleToInfected), Eigen::Index(InfectionState::Infected), dt); +} + +void NormModel::update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, + bool Transitionispossible, ScalarType dt) +{ + Eigen::Index current_time_index = populations.get_num_time_points() - 1; + ScalarType current_time_age = (ScalarType)current_time_index * dt; + Eigen::Index calc_time_index = current_time_index; + if (Transitionispossible) { + calc_time_index = compparameters->m_transitiondistributions.size() - 1; + } + + ScalarType sum = 0; + + Eigen::Index starting_point = std::max(0, (int)current_time_index - (int)calc_time_index); + for (int i = starting_point; i < current_time_index; i++) { + ScalarType state_age_i = (ScalarType)i * dt; + ScalarType sum_inflows = 0; + for (const InfectionTransition& inflow : IncomingFlows) { + sum_inflows += transitions[i + 1][(int)inflow]; + } + if (Transitionispossible) { + sum += compparameters->m_transitiondistributions[current_time_index - i] * + (sum_inflows + (compparameters->parameters.get() + + transitions[i][(Eigen::Index)InfectionTransition::InfectedToDead] - + compparameters->parameters.get()) * + populations[i][(int)infectionState]) * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + else { + sum += (sum_inflows + (compparameters->parameters.get() + + transitions[i][(Eigen::Index)InfectionTransition::InfectedToDead] - + compparameters->parameters.get()) * + populations[i][(int)infectionState]) * + std::exp(-compparameters->parameters.get() * (current_time_age - state_age_i)); + } + } + if (Transitionispossible) { + if (current_time_index <= calc_time_index) { + populations.get_last_value()[(int)infectionState] = + dt * sum + compparameters->m_transitiondistributions[current_time_index] * + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } + else { + populations.get_last_value()[(int)infectionState] = dt * sum; + } + } + else { + populations.get_last_value()[(int)infectionState] = + dt * sum + populations[0][(int)infectionState] * + std::exp(-compparameters->parameters.get() * (current_time_age)); + } +} + +void NormModel::update_compartments(ScalarType dt) +{ + + // Infected + update_compartment_with_sum(InfectionState::Infected, {InfectionTransition::SusceptibleToInfected}, true, dt); + // Recovered + update_compartment_with_sum(InfectionState::Recovered, + { + InfectionTransition::InfectedToRecovered, + }, + false, dt); +} +void NormModel::compute_size() +{ + ScalarType sum = 0; + for (int state = 0; state < Eigen::Index(InfectionState::Count) - 1; state++) { + sum += populations.get_last_value()[state]; + } + m_size.get_last_value()[0] = sum; +} + +void NormModel::compute_forceofinfection(ScalarType dt) +{ + + Eigen::Index num_time_points = populations.get_num_time_points(); + ScalarType current_time = populations.get_last_time(); + + // Determine the starting point of the for loop. + Eigen::Index starting_point1 = std::max(0, (int)num_time_points - 1 - (int)compparameters->m_infectivity.size()); + + ScalarType sum = 0.; + // Compute the first sum in the force of infection term. + for (Eigen::Index i = starting_point1; i < num_time_points - 1; i++) { + sum += compparameters->m_infectivity[num_time_points - 1 - i] * + (populations[i + 1][(int)InfectionState::Susceptible] * m_forceofinfection[i][0] + + (compparameters->parameters.get() + + transitions[i + 1][(int)InfectionTransition::InfectedToDead] - + compparameters->parameters.get()) * + populations[i + 1][(int)InfectionState::Infected]); + } + + m_forceofinfection.get_last_value()[0] = + dt * sum * + (compparameters->parameters.get().eval(current_time) * + compparameters->parameters.get().get_cont_freq_mat().get_matrix_at(current_time)(0, 0)); + + // Add inital functions for the force of infection in case they still have an influence. + if (num_time_points <= (int)compparameters->m_NormFoI_0.size()) { + m_forceofinfection.get_last_value()[0] += compparameters->m_NormFoI_0[num_time_points - 1]; + } +} + +} // namespace endisird + +} // namespace mio \ No newline at end of file diff --git a/cpp/models/ide_endemic_sird/normalized_model.h b/cpp/models/ide_endemic_sird/normalized_model.h new file mode 100644 index 0000000000..e6cca2560d --- /dev/null +++ b/cpp/models/ide_endemic_sird/normalized_model.h @@ -0,0 +1,158 @@ +#ifndef IDE_END_SIRD_NORMMODEL_H +#define IDE_END_SIRD_NORMMODEL_H + +#include "ide_endemic_sird/infection_state.h" +#include "ide_endemic_sird/parameters.h" +#include "ide_endemic_sird/computed_parameters.h" +#include "ide_endemic_sird/model.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" + +#include "vector" +#include +#include +#include + +namespace mio +{ +namespace endisird +{ + +// Forward declaration of friend classes/functions of Model. +class NormModel; +class Simulation; + +class NormModel +{ + using ParameterSet = Parameters; + +public: + /** + * @brief Constructor to create an endemic IDE_SECIR model. + * + * @param[in] TODO!!!!!! + */ + + NormModel(CompParameters const& compparams); + + /** + * @brief Checks constraints on model parameters and initial data. + * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. + */ + bool check_constraints() const; + + // ---- Public parameters. ---- + std::shared_ptr compparameters; + TimeSeries + transitions; ///< TimesSeries containing points of time and the corresponding number of individuals transitioning + // from one #InfectionState to another as defined in #Infection%s. + TimeSeries populations; ///< TimeSeries containing points of time and the corresponding number of people + // in defined #InfectionState%s. + +private: + // ---- Functionality for the iterations of a simulation. + + /** + * @brief Computes number of Susceptibles for the current last time in populations. + * + * Number is computed by using the previous number of Susceptibles, total population of the last time point and + * the force of infection (also from the last time point). + * Number is stored at the matching index in populations. + * + * @param[in] dt Time discretization step size. + */ + void compute_susceptibles(ScalarType dt); + + /** + * @brief Computes sizeof a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the time index + * current_time_index. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] current_time_index The time index the transition should be computed for. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, Eigen::Index current_time_index, ScalarType dt); + + /** + * @brief Computes size of a flow + * + * Computes size of one Transition from #InfectionTransition, specified in idx_InfectionTransitions, for the current + * last time value in transitions. + * + * @param[in] idx_InfectionTransitions Specifies the considered transition from #InfectionTransition + * @param[in] idx_IncomingFlow Index of the transition in #InfectionTransition, which goes to the considered starting + * compartment of the transition specified in idx_InfectionTransitions. Size of considered flow is calculated via + * the value of this incoming flow. + * @param[in] idx_CurrentCompartment Index of the Compartment we flow out. + * @param[in] dt Time discretization step size. + */ + void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, + Eigen::Index idx_CurrentCompartment, ScalarType dt); + + /** + * @brief Sets all required transitions for the current last timestep in transitions. + * + * New values are stored in transitions. Most values are computed via the function compute_flow() + * + * @param[in] dt time discretization step size. + */ + void flows_currents_timestep(ScalarType dt); + + /** + * @brief Updates the values of one compartment, specified in infectionState, using all past transitions. + * + * New value is stored in populations. The values are calculated using all past values for the incoming flows + * including the current time step. + * Therefore the flows of the current time step should be calculated before using this function. + * @param[in] infectionState Specifies the #InfectionState we want to update. + * @param[in] IncomingFlows + * @param[in] NaturalDeathispossible Boolian that determines if there is the possibility of Natural Death in infectionState. + * @param[in] Transitionispossible Boolian that determines if it is possible to transition from the current InfectionState + * into another + * @param[in] dt + */ + void update_compartment_with_sum(InfectionState infectionState, + std::vector const& IncomingFlows, bool Transitionispossible, + ScalarType dt); + + /** + * @brief Updates the values of all compartments except Susceptible + * + * New value is stored in populations. Values are computed via the function update_compartment_from_flow + */ + void update_compartments(ScalarType dt); + + void compute_size(); + + /** + * @brief Computes the force of infection for the current last time in transitions. + * + * Computed value is stored in m_forceofinfection. + * + * @param[in] dt Time discretization step size. + */ + void compute_forceofinfection(ScalarType dt); + + // ---- Private parameters. ---- + + TimeSeries m_forceofinfection{ + TimeSeries(1)}; ///< TimeSeries containing the Force of infection term for every time point, + TimeSeries m_size{TimeSeries(1)}; + + // ---- Friend classes/functions. ---- + // In the Simulation class, the actual simulation is performed which is why it needs access to the here + // defined (and private) functions to solve the model equations. + friend class Simulation; +}; + +} // namespace endisird +} // namespace mio + +#endif // IDE_END_SIRD_NORMMODEL_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_sird/parameters.h b/cpp/models/ide_endemic_sird/parameters.h new file mode 100644 index 0000000000..e9b64fc3a3 --- /dev/null +++ b/cpp/models/ide_endemic_sird/parameters.h @@ -0,0 +1,244 @@ +#ifndef IDE_END_SIRD_PARAMS_H +#define IDE_END_SIRD_PARAMS_H + +#include "memilio/config.h" +#include "memilio/epidemiology/contact_matrix.h" +#include "memilio/math/floating_point.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/parameter_set.h" +#include "ide_endemic_sird/infection_state.h" +#include "memilio/epidemiology/state_age_function.h" +#include "memilio/epidemiology/uncertain_matrix.h" + +#include +#include +#include + +namespace mio +{ +namespace endisird +{ + +/************************************************** +* Define Parameters of the IDE-END-SIRD model * +**************************************************/ + +/** + * @brief Transitions distribution for each transition in #Infection Transition. + * + * we use as a default SmootherCosine functions for all transitions with m_paramter=2. + */ + +struct TransitionDistributions { + using Type = std::vector; + static Type get_default() + { + SmootherCosine smoothcos(2.0); + StateAgeFunctionWrapper delaydistribution(smoothcos); + return Type((int)InfectionTransition::Count, delaydistribution); + } + + static std::string mame() + { + return "TransitionDistributions"; + } +}; + +/** + * @brief Defines the probability for each possible transitions to take this transition- + */ + +struct TransitionProbabilities { + /*For consistency, also define TransitionProbabilities for each transition in #InfectionTransition. + Transition Probabilities should be set to 1 if there is no possible other flow from starting compartment.*/ + using Type = std::vector; + static Type get_default() + { + std::vector probs((int)InfectionTransition::Count, 0.5); + // Set the following probablities to 1 as there is no other option to go anywhere else. + probs[Eigen::Index(InfectionTransition::SusceptibleToInfected)] = 1; + return Type(probs); + } + static std::string name() + { + return "TransitionsProbabilities"; + } +}; + +/** + * @brief The contact patterns within the society are modelled using an UncertainContactMatrix. + */ +struct ContactPatterns { + using Type = UncertainContactMatrix; + static Type get_default() + { + ContactMatrixGroup contact_matrix = ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.)); + return Type(contact_matrix); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/** +* @brief Probability of getting infected from a contact. +*/ +struct TransmissionProbabilityOnContact { + using Type = StateAgeFunctionWrapper; + static Type get_default() + { + ConstantFunction constfunc(1.0); + return Type(constfunc); + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/** +* @brief The risk of infection from symptomatic cases in the SECIR model. +*/ +struct RiskOfInfectionFromSymptomatic { + using Type = StateAgeFunctionWrapper; + static Type get_default() + { + ConstantFunction constfunc(1.0); + return Type(constfunc); + } + static std::string name() + { + return "RiskOfInfectionFromSymptomatic"; + } +}; + +/** + * @brief The natural birth rate. + */ +struct NaturalBirthRate { + using Type = ScalarType; + static Type get_default() + { + return Type(5e-5); + } + static std::string name() + { + return "NaturalBirthRate"; + } +}; + +/** + * @brief The natural death rate. + */ +struct NaturalDeathRate { + using Type = ScalarType; + static Type get_default() + { + return Type(3e-5); + } + static std::string name() + { + return "NaturalDeathRate"; + } +}; + +// Define Parameterset for IDE-END-SECIR model. +using ParametersBase = + ParameterSet; + +/** + * @brief Parameters of an endemic SECIR/SECIHURD model. + */ + +class Parameters : public ParametersBase +{ +public: + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and logs an error. + */ + bool check_constraints() const + { + // For paramerers potentially depending on the infectious age, the values are checked equidistantly on a + // realistic maximum window. + size_t infectious_window_check = 50; // parameter defining minmal window on x-axis. + + //Check if all the probabilitys are between 0 and 1. + for (size_t i = 0; i < infectious_window_check; i++) { + if (this->get().eval((static_cast(i))) < 0.0 || + this->get().eval(static_cast(i)) > 1.0) { + log_error( + "Constraint check: TransmissionProbabilityOnContact smaller {:d} or larger {:d} at some time {:d}", + 0, 1, i); + return true; + } + } + + for (size_t i = 0; i < infectious_window_check; i++) { + if (this->get().eval((static_cast(i))) < 0.0 || + this->get().eval(static_cast(i)) > 1.0) { + log_error( + "Constraint check: RiskOfInfectionFromSymptomatic smaller {:d} or larger {:d} at some time {:d}", 0, + 1, i); + return true; + } + } + + for (size_t i = 0; i < static_cast(InfectionTransition::Count); i++) { + if (this->get()[i] < 0.0 || this->get()[i] > 1.0) { + log_error("Constraint check: One parameter in TransitionProbabilities smaller {:d} or larger {:d}", 0, + 1); + return true; + } + } + + // The TransitionProbabilitie SusceptibleToInfected should be 1. + if (!floating_point_equal( + this->get()[static_cast(InfectionTransition::SusceptibleToInfected)], 1.0, + 1e-14)) { + log_error("Constraint check: Parameter transition probability for SusceptibleToExposed unequal to {:d}", 1); + return true; + } + + // All TransitionProbabilities where the Transition starts in the same compartment should sum up to 1. + + if (!floating_point_equal( + this->get()[static_cast(InfectionTransition::InfectedToDead)] + + this->get()[static_cast(InfectionTransition::InfectedToRecovered)], + 1.0, 1e-14)) { + log_error("Constraint check: Sum of transition probability for InfectedToDead and " + "InfectedToRecovered not equal to {:d}", + 1); + return true; + } + + /* The first entry of TransitionDistributions is not checked because the distribution S->E is never used + (and it makes no sense to use the distribution). The support does not need to be valid.*/ + for (size_t i = 1; i < static_cast(InfectionTransition::Count); i++) { + if (floating_point_less(this->get()[i].get_support_max(10), 0.0, 1e-14)) { + log_error("Constraint check: One function in TransitionDistributions has invalid support."); + return true; + } + } + + if (this->get() < 0.0) { + log_warning("Constraint check: Parameter NaturalBirthRate should be greater than {:d}", 0); + return true; + } + + if (this->get() < 0.0) { + log_warning("Constraint check: Parameter NaturalDeathRate should be greater than {:d}", 0); + return true; + } + + return false; + } +}; + +} // namespace endisird + +} // namespace mio + +#endif //IDE_END_SIRD_PARAMS_H \ No newline at end of file diff --git a/cpp/models/ide_endemic_sird/simulation.cpp b/cpp/models/ide_endemic_sird/simulation.cpp new file mode 100644 index 0000000000..921e9e6bb4 --- /dev/null +++ b/cpp/models/ide_endemic_sird/simulation.cpp @@ -0,0 +1,111 @@ + +#include "ide_endemic_sird/simulation.h" +#include "ide_endemic_sird/computed_parameters.h" +#include "ide_endemic_sird/model.h" +#include "ide_endemic_sird/normalized_model.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include + +namespace mio +{ +namespace endisird +{ + +void Simulation::advance(ScalarType tmax) +{ + mio::log_info("Simulating IDE-END-SECIR from t0 = {} until tmax = {} with dt = {}.", + m_model->transitions.get_last_time(), tmax, m_dt); + + m_model->compparameters->set_transitiondistributions_support_max(m_dt); + m_model->compparameters->set_transitiondistributions(m_dt); + m_model->compparameters->set_transitiondistributions_derivative(m_dt); + m_model->compparameters->set_infectivity(m_dt); + m_model->compparameters->set_FoI_0(m_dt); + m_model->compparameters->set_reproductionnumber_c(m_dt); + m_model->compparameters->set_T(m_dt); + m_model->compparameters->set_W(m_dt); + + m_normmodel->compparameters->set_transitiondistributions_support_max(m_dt); + m_normmodel->compparameters->set_transitiondistributions(m_dt); + m_normmodel->compparameters->set_transitiondistributions_derivative(m_dt); + m_normmodel->compparameters->set_infectivity(m_dt); + m_normmodel->compparameters->set_FoI_0(m_dt); + m_normmodel->compparameters->set_reproductionnumber_c(m_dt); + m_normmodel->compparameters->set_T(m_dt); + m_normmodel->compparameters->set_W(m_dt); + m_normmodel->compparameters->set_equilibria(); + + m_difference_normalizedcompartments.add_time_point(0); + m_difference_normalizedFoI.add_time_point(0); + compute_difference_normalizations(); + + // For every time step: + while (m_model->transitions.get_last_time() < tmax - m_dt / 2) { + + m_difference_normalizedcompartments.add_time_point(m_difference_normalizedcompartments.get_last_time() + m_dt); + m_difference_normalizedFoI.add_time_point(m_difference_normalizedFoI.get_last_time() + m_dt); + //standard model : + m_model->transitions.add_time_point(m_model->transitions.get_last_time() + m_dt); + m_model->populations.add_time_point(m_model->populations.get_last_time() + m_dt); + m_model->m_forceofinfection.add_time_point(m_model->m_forceofinfection.get_last_time() + m_dt); + m_model->m_totalpopulation.add_time_point(m_model->m_totalpopulation.get_last_time() + m_dt); + m_model->m_totalpopulation_derivative.add_time_point(m_model->m_totalpopulation_derivative.get_last_time() + + m_dt); + m_model->m_normalizedpopulations.add_time_point(m_model->m_normalizedpopulations.get_last_time() + m_dt); + + // Compute susceptibles: + m_model->compute_susceptibles(m_dt); + + // Compute flows: + m_model->flows_currents_timestep(m_dt); + + // Compute remaining compartments: + m_model->update_compartments(m_dt); + + // Compute m_populationsize: + m_model->compute_populationsize(); + + // Compute normalized compartments: + m_model->compute_normalizedcompartments(); + + // Compute m_forceofinfection; + m_model->compute_forceofinfection(m_dt); + + // normalized model: + m_normmodel->transitions.add_time_point(m_normmodel->transitions.get_last_time() + m_dt); + m_normmodel->populations.add_time_point(m_normmodel->populations.get_last_time() + m_dt); + m_normmodel->m_forceofinfection.add_time_point(m_normmodel->m_forceofinfection.get_last_time() + m_dt); + m_normmodel->m_size.add_time_point(m_normmodel->m_size.get_last_time() + m_dt); + + // Compute susceptibles: + m_normmodel->compute_susceptibles(m_dt); + + // Compute flows: + m_normmodel->flows_currents_timestep(m_dt); + + // Compute remaining compartments: + m_normmodel->update_compartments(m_dt); + + m_normmodel->compute_size(); + + // Compute m_forceofinfection; + m_normmodel->compute_forceofinfection(m_dt); + + //Compute the difference of the two normalized compartments. + compute_difference_normalizations(); + } +} + +TimeSeries simulate(ScalarType tmax, ScalarType dt, CompParameters const& m_compparams, + Model const& m_model, NormModel const& m_normmodel) +{ + m_model.check_constraints(); + Simulation sim(m_compparams, m_model, m_normmodel, dt); + sim.advance(tmax); + return sim.get_compartments(); +} + +} // namespace endisird + +} // namespace mio diff --git a/cpp/models/ide_endemic_sird/simulation.h b/cpp/models/ide_endemic_sird/simulation.h new file mode 100644 index 0000000000..7888303d1f --- /dev/null +++ b/cpp/models/ide_endemic_sird/simulation.h @@ -0,0 +1,257 @@ +#ifndef IDE_END_SIRD_SIMULATION_H +#define IDE_END_SIRD_SIMULATION_H + +#include "ide_endemic_sird/computed_parameters.h" +#include "ide_endemic_sird/model.h" +#include "ide_endemic_sird/normalized_model.h" +#include "memilio/config.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/time_series.h" +#include +#include + +namespace mio +{ +namespace endisird +{ +class Simulation +{ +public: + Simulation(CompParameters const& compparams, Model const& model, NormModel const& normmodel, ScalarType dt = 0.1) + : m_compparams(std::make_unique(compparams)) + , m_model(std::make_unique(model)) + , m_normmodel(std::make_unique(normmodel)) + , m_dt(dt) + { + } + + /** + * @brief Run the simulation until time tmax. + */ + void advance(ScalarType tmax); + + /** + * @brief Get the result of the simulation for the compartments of m_model + * Return the number of persons in all #InfectionState%s + */ + TimeSeries get_compartments() + { + return m_model->populations; + } + + /** + * @brief Get the result of the simulation for the normalized compartments, where we use m_model and the compartments + * computed using the sum scheme. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries get_normalizedcompartments() + { + return m_model->m_normalizedpopulations; + } + /** + * @brief Get the result of the simulation for the compartments of the normalized model. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries get_normmodel_compartments() + { + return m_normmodel->populations; + } + + /** + * @brief Get the result of the simulation for the compartments. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries& get_compartments() const + { + return m_model->populations; + } + + /** + * @brief Get the result of the simulation for the normalized compartments, where we use m_model. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries& get_normalizedcompartments() const + { + return m_model->m_normalizedpopulations; + } + + /** + * @brief Get the result of the simulation for the compartments of the normalized model. + * Return the number of persons in all #InfectionState%s + */ + TimeSeries& get_normmodel_compartments() const + { + return m_normmodel->populations; + } + + /** + * @brief Get the transitions between the different #InfectionState%s for m_model. + */ + TimeSeries const& get_transitions() + { + return m_model->transitions; + } + /** + * @brief Get the transitions between the different #InfectionState%s of m_normmodel. + */ + TimeSeries const& get_normmodel_transitions() + { + return m_normmodel->transitions; + } + + /** + * @brief Get the force of infection term of m_model. + */ + TimeSeries const& get_forceofinfections() + { + return m_model->m_forceofinfection; + } + + /** + * @brief Get the force of infection term of m_normmodel. + */ + TimeSeries const& get_normmodel_forceofinfections() + { + return m_normmodel->m_forceofinfection; + } + + /** + * @brief Get the total population of m_model. + */ + TimeSeries const& get_totalpopulations() + { + return m_model->m_totalpopulation; + } + + TimeSeries const& get_size() + { + return m_normmodel->m_size; + } + + /** + * @brief Get the derivative of the total population of m_model. + */ + TimeSeries const& get_totalpopulations_derivative() + { + return m_model->m_totalpopulation_derivative; + } + + /** + * @brief Get the reproduction numer. + */ + ScalarType const& get_reproductionnumber_c() + { + return m_model->compparameters->m_reproductionnumber_c; + } + + /** + * @brief Get T. + */ + std::vector const& get_T() + { + return m_model->compparameters->m_T; + } + + /** + * @brief Get W. + */ + ScalarType const& get_W() + { + return m_model->compparameters->m_W; + } + + std::vector> const& get_Equilibirum_compartments() + { + return m_normmodel->compparameters->m_compartments_equilibrium; + } + + std::vector> const& get_Equilibirum_transitions() + { + return m_normmodel->compparameters->m_transitions_equilibrium; + } + + std::vector const& get_Equilibirum_FoI() + { + return m_normmodel->compparameters->m_FoI_equilibrium; + } + + /** + * @brief returns the simulation model used in simulation. + */ + const Model& get_model() const + { + return *m_model; + } + + /** + * @brief returns the simulation model used in simulation. + */ + Model& get_model() + { + return *m_model; + } + + /** + * @brief returns the simulation normmodel used in simulation. + */ + const NormModel& get_normmodel() const + { + return *m_normmodel; + } + + /** + * @brief returns the simulation normmodel used in simulation. + */ + NormModel& get_normmodel() + { + return *m_normmodel; + } + + /** + * @brief get the size of the time step of the simulation. + */ + ScalarType get_stepsize() + { + return m_dt; + } + + void compute_difference_normalizations() + { + for (int state = 0; state < (int)InfectionState::Count - 1; state++) { + m_difference_normalizedcompartments.get_last_value()[state] = + std::abs(m_normmodel->populations.get_last_value()[state] - + m_model->m_normalizedpopulations.get_last_value()[state]); + } + m_difference_normalizedFoI.get_last_value()[0] = std::abs(m_normmodel->m_forceofinfection.get_last_value()[0] - + m_model->m_forceofinfection.get_last_value()[0]); + } + + TimeSeries const& get_difference_normalizationcomp() + { + return m_difference_normalizedcompartments; + } + + TimeSeries const& get_difference_normalizationFoi() + { + return m_difference_normalizedFoI; + } + +private: + std::unique_ptr m_compparams; ///< Unique pointer to the computed Parameters. + std::unique_ptr m_model; ///< Unique pointer to the simulated Model. + std::unique_ptr m_normmodel; ///< Unique pointer to the simulated normalized Model. + ScalarType m_dt; ///< Time step used for numerical computations in simulation. + TimeSeries m_difference_normalizedcompartments{TimeSeries( + (int)InfectionState::Count - 1)}; ///< TimeSeries containing the difference of the compartments + // computed by NormModel and the normalized compartments computed in Model. + TimeSeries m_difference_normalizedFoI{ + TimeSeries(1)}; ///< TimeSeries containing the difference of the force of infection terms + // computed by NormModel and Model. +}; + +TimeSeries simulate(ScalarType tmax, ScalarType dt, Model const& model); + +} // namespace endisird +} // namespace mio + +#endif //IDE_END_SIRD_SIMULATION_H \ No newline at end of file diff --git a/cpp/models/ide_secir/model.h b/cpp/models/ide_secir/model.h index 98cff9ea7b..bb96e033e0 100644 --- a/cpp/models/ide_secir/model.h +++ b/cpp/models/ide_secir/model.h @@ -68,7 +68,7 @@ class Model CustomIndexArray total_confirmed_cases_init = CustomIndexArray()); // ---- Additional functionality such as constraint checking, setters and getters, etc. ---- - /** + /**c * @brief Checks constraints on model parameters and initial data. * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. */ diff --git a/cpp/models/ide_secir/simulation.h b/cpp/models/ide_secir/simulation.h index ce5cce5f28..688347fce8 100644 --- a/cpp/models/ide_secir/simulation.h +++ b/cpp/models/ide_secir/simulation.h @@ -1,3 +1,4 @@ + /* * Copyright (C) 2020-2025 MEmilio *