diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index e6235e3b4f..72bf24afe2 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -151,6 +151,10 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/ode_seir) add_subdirectory(models/ode_seair) add_subdirectory(models/ode_sir) + add_subdirectory(models/ode_sir_mobility) + # add_subdirectory(models/ode_seir_mobility_massaction) + add_subdirectory(models/ode_seir_mobility_improved) + add_subdirectory(models/ode_seir_mobility) add_subdirectory(models/sde_sir) add_subdirectory(models/sde_sirs) add_subdirectory(models/sde_seirvv) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index f18ab284f7..1d5cc48c03 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -30,6 +30,32 @@ add_executable(sde_sir_example sde_sir.cpp) target_link_libraries(sde_sir_example PRIVATE memilio sde_sir) target_compile_options(sde_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_sir_mobility_example ode_sir_mobility.cpp) +target_link_libraries(ode_sir_mobility_example PRIVATE memilio ode_sir_mobility) +target_compile_options(ode_sir_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(ode_seir_mobility_example ode_seir_mobility.cpp) +target_link_libraries(ode_seir_mobility_example PRIVATE memilio ode_seir_mobility) +target_compile_options(ode_seir_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +if (MEMILIO_ENABLE_OPENMP) + add_executable(ode_seir_mobility_timing ode_seir_mobility_timing.cpp) + target_link_libraries(ode_seir_mobility_timing PRIVATE memilio ode_seir_mobility_improved) + target_compile_options(ode_seir_mobility_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +endif() + +add_executable(ode_seir_mobility_steps ode_seir_mobility_steps.cpp) +target_link_libraries(ode_seir_mobility_steps PRIVATE memilio ode_seir_mobility_improved) +target_compile_options(ode_seir_mobility_steps PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(basic_reproduction_numbers basic_reproduction_numbers.cpp) +target_link_libraries(basic_reproduction_numbers PRIVATE memilio ode_seir) +target_compile_options(basic_reproduction_numbers PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(ode_seir_mobility_example_improved ode_seir_mobility_improved.cpp) +target_link_libraries(ode_seir_mobility_example_improved PRIVATE memilio ode_seir_mobility_improved) +target_compile_options(ode_seir_mobility_example_improved PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(sde_sirs_example sde_sirs.cpp) target_link_libraries(sde_sirs_example PRIVATE memilio sde_sirs) target_compile_options(sde_sirs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) @@ -84,6 +110,20 @@ add_executable(graph_example graph.cpp) target_link_libraries(graph_example PRIVATE memilio ode_seir) target_compile_options(graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(graph_example_extended graph_extended.cpp) +target_link_libraries(graph_example_extended PRIVATE memilio ode_seir) +target_compile_options(graph_example_extended PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(graph_steps graph_steps.cpp) +target_link_libraries(graph_steps PRIVATE memilio ode_seir) +target_compile_options(graph_steps PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +if (MEMILIO_ENABLE_OPENMP) + add_executable(graph_timing graph_timing.cpp) + target_link_libraries(graph_timing PRIVATE memilio ode_seir) + target_compile_options(graph_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +endif() + add_executable(graph_stochastic_mobility_example graph_stochastic_mobility.cpp) target_link_libraries(graph_stochastic_mobility_example PRIVATE memilio ode_secir) target_compile_options(graph_stochastic_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/basic_reproduction_numbers.cpp b/cpp/examples/basic_reproduction_numbers.cpp new file mode 100644 index 0000000000..71cabd73f3 --- /dev/null +++ b/cpp/examples/basic_reproduction_numbers.cpp @@ -0,0 +1,196 @@ + +#include "models/ode_seir/model.h" +#include "models/ode_seir_mobility/model.h" +// #include "models/ode_seir_mobility_improved/model.h" + +#include "memilio/math/euler.h" +#include "memilio/compartments/simulation.h" +#include "memilio/utils/custom_index_array.h" + +Eigen::MatrixXd get_contact_matrix() +{ + Eigen::MatrixXd contact_matrix_eigen(6, 6); + contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, 0.0161, + 0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, 0.0458, 0.1939, + 0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724; + + return contact_matrix_eigen; +} + +const ScalarType TimeExposed[] = {3.335, 3.335, 3.335, 3.335, 3.335, 3.335}; +const ScalarType TimeInfected[] = {8.0096875, 8.0096875, 8.2182, 8.1158, 8.033, 7.985}; +const ScalarType TransmissionProbabilityOnContact[] = {0.03, 0.06, 0.06, 0.06, 0.09, 0.175}; + +void seir(size_t number_regions, ScalarType tmax) +{ + mio::set_log_level(mio::LogLevel::off); + ScalarType t0 = 0.; + ScalarType dt = 0.1; + ScalarType number_age_groups = 6; + + mio::oseir::Model model(number_age_groups); + auto& population = model.populations; + + for (size_t j = 0; j < number_age_groups; j++) { + + population[{mio::AgeGroup(j), mio::oseir::InfectionState::Susceptible}] = number_regions * 10000; + } + population[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] += 100; + population[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] -= 100; + + mio::ContactMatrixGroup& contact_matrix = + model.parameters.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline() = get_contact_matrix(); + + for (size_t j = 0; j < number_age_groups; j++) { + model.parameters.template get>()[mio::AgeGroup(j)] = TimeExposed[j]; + model.parameters.template get>()[mio::AgeGroup(j)] = TimeInfected[j]; + model.parameters.template get>()[mio::AgeGroup(j)] = + TransmissionProbabilityOnContact[j]; + } + + std::shared_ptr> integrator = std::make_shared>(); + + auto result = simulate(t0, tmax, dt, model, integrator); + + auto basic_reproduction_number = model.get_reproduction_number(t0, result).value(); + std::cout << "\"SEIR\": " << basic_reproduction_number << ", " << std::endl; +} + +void wang(size_t number_regions, ScalarType tmax) +{ + mio::set_log_level(mio::LogLevel::off); + ScalarType t0 = 0.; + ScalarType dt = 0.1; + ScalarType number_age_groups = 6; + + mio::oseirmobility::Model model(number_regions, number_age_groups); + auto& population = model.populations; + + for (size_t j = 0; j < number_age_groups; j++) { + for (size_t i = 0; i < number_regions; i++) { + population[{mio::oseirmobility::Region(i), mio::AgeGroup(j), + mio::oseirmobility::InfectionState::Susceptible}] = 10000; + } + } + population[{mio::oseirmobility::Region(0), mio::AgeGroup(0), mio::oseirmobility::InfectionState::Exposed}] += 100; + population[{mio::oseirmobility::Region(0), mio::AgeGroup(0), mio::oseirmobility::InfectionState::Susceptible}] -= + 100; + + double fraction_commuter = 1. / (2 * number_regions); + Eigen::MatrixXd mobility_data_commuter = + Eigen::MatrixXd::Constant(number_regions, number_regions, fraction_commuter) - + fraction_commuter * + Eigen::MatrixXd::Identity(number_regions, number_regions); // Ensure that the diagonal is zero + for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) { + mobility_data_commuter(county_idx_i, county_idx_i) = 1 - mobility_data_commuter.rowwise().sum()(county_idx_i); + } + model.parameters.template get>().get_cont_freq_mat()[0].get_baseline() = + mobility_data_commuter; + + mio::ContactMatrixGroup& contact_matrix = + model.parameters.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline() = get_contact_matrix(); + + for (size_t j = 0; j < number_age_groups; j++) { + model.parameters.template get>()[mio::AgeGroup(j)] = TimeExposed[j]; + model.parameters.template get>()[mio::AgeGroup(j)] = TimeInfected[j]; + model.parameters.template get>()[mio::AgeGroup(j)] = + TransmissionProbabilityOnContact[j]; + } + + std::shared_ptr> integrator = std::make_shared>(); + + auto result = simulate(t0, tmax, dt, model, integrator); + + auto basic_reproduction_number = model.get_reproduction_number(t0, result).value(); + std::cout << "\"Wang\": " << basic_reproduction_number << "}" << std::endl; +} + +// void metapopulation(size_t number_regions, ScalarType tmax) +// { +// mio::set_log_level(mio::LogLevel::off); +// ScalarType t0 = 0.; +// ScalarType dt = 0.1; +// ScalarType number_age_groups = 6; + +// mio::oseirmobilityimproved::Model model(number_regions, number_age_groups); +// auto& population = model.populations; + +// for (size_t j = 0; j < number_age_groups; j++) { +// for (size_t i = 0; i < number_regions; i++) { +// population[{mio::oseirmobilityimproved::Region(i), mio::AgeGroup(j), +// mio::oseirmobilityimproved::InfectionState::Susceptible}] = 10000; +// } +// } +// population[{mio::oseirmobilityimproved::Region(0), mio::AgeGroup(0), +// mio::oseirmobilityimproved::InfectionState::Exposed}] += 100; +// population[{mio::oseirmobilityimproved::Region(0), mio::AgeGroup(0), +// mio::oseirmobilityimproved::InfectionState::Susceptible}] -= 100; + +// double fraction_commuter = 1. / (2 * number_regions); +// Eigen::MatrixXd mobility_data_commuter = +// Eigen::MatrixXd::Constant(number_regions, number_regions, fraction_commuter) - +// fraction_commuter * +// Eigen::MatrixXd::Identity(number_regions, number_regions); // Ensure that the diagonal is zero +// for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) { +// mobility_data_commuter(county_idx_i, county_idx_i) = 1 - mobility_data_commuter.rowwise().sum()(county_idx_i); +// } +// model.parameters.template get>() +// .get_cont_freq_mat()[0] +// .get_baseline() = mobility_data_commuter; + +// mio::ContactMatrixGroup& contact_matrix = +// model.parameters.template get>().get_cont_freq_mat(); +// contact_matrix[0].get_baseline() = get_contact_matrix(); + +// for (size_t j = 0; j < number_age_groups; j++) { +// model.parameters.template get>()[mio::AgeGroup(j)] = TimeExposed[j]; +// model.parameters.template get>()[mio::AgeGroup(j)] = TimeInfected[j]; +// model.parameters +// .template get>()[mio::AgeGroup(j)] = +// TransmissionProbabilityOnContact[j]; +// } + +// mio::ContactMatrixGroup& commuting_strengths = +// model.parameters.template get>().get_cont_freq_mat(); + +// auto& population_after_commuting = model.m_population_after_commuting; +// for (size_t region_n = 0; region_n < number_regions; ++region_n) { +// for (size_t age = 0; age < number_age_groups; ++age) { +// double population_n = 0; +// for (size_t state = 0; state < (size_t)mio::oseirmobilityimproved::InfectionState::Count; state++) { +// population_n += population[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age), +// mio::oseirmobilityimproved::InfectionState(state)}]; +// } +// population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] += +// population_n; +// for (size_t region_m = 0; region_m < number_regions; ++region_m) { +// population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] -= +// commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; +// population_after_commuting[{mio::oseirmobilityimproved::Region(region_m), mio::AgeGroup(age)}] += +// commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; +// } +// } +// } + +// std::shared_ptr> integrator = std::make_shared>(); + +// auto result = simulate(t0, tmax, dt, model, integrator); + +// auto basic_reproduction_number = model.get_reproduction_number(t0, result).value(); +// std::cout << "\"Metapopulation\": " << basic_reproduction_number << "}" << std::endl; +// } + +int main() +{ + const ScalarType tmax = 1.; + size_t num_regions = 150; + + std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl; + + seir(num_regions, tmax); + wang(num_regions, tmax); + // metapopulation(num_regions, tmax); + return 0; +} \ No newline at end of file diff --git a/cpp/examples/graph.cpp b/cpp/examples/graph.cpp index ad944254af..b060e1bcb2 100644 --- a/cpp/examples/graph.cpp +++ b/cpp/examples/graph.cpp @@ -22,11 +22,12 @@ #include "ode_seir/parameters.h" #include "memilio/mobility/metapopulation_mobility_instant.h" #include "memilio/compartments/simulation.h" +#include "memilio/io/result_io.h" int main() { const auto t0 = 0.; - const auto tmax = 10.; + const auto tmax = 15.; const auto dt = 0.5; //time step of mobility, daily mobility every second step mio::oseir::Model<> model(1); @@ -34,9 +35,11 @@ int main() // set population model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000; + model.parameters.set>(1.); + // set transition times - model.parameters.set>(1); - model.parameters.set>(1); + model.parameters.set>(3.); + model.parameters.set>(5.); // set contact matrix mio::ContactMatrixGroup& contact_matrix = model.parameters.get>().get_cont_freq_mat(); @@ -47,9 +50,9 @@ int main() auto model_group2 = model; //some contact restrictions in group 1 - mio::ContactMatrixGroup& contact_matrix1 = - model_group1.parameters.get>().get_cont_freq_mat(); - contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5)); + // mio::ContactMatrixGroup& contact_matrix1 = + // model_group1.parameters.get>().get_cont_freq_mat(); + // contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5)); //infection starts in group 1 model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990; @@ -58,12 +61,29 @@ int main() mio::Graph>>, mio::MobilityEdge<>> g; g.add_node(1001, model_group1, t0); g.add_node(1002, model_group2, t0); - g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01)); + for (auto& node : g.nodes()) { + node.property.get_simulation().set_integrator(std::make_shared>()); + } + g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.05)); g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01)); auto sim = mio::make_mobility_sim(t0, dt, std::move(g)); sim.advance(tmax); + auto result_graph = std::move(sim).get_graph(); + auto result = mio::interpolate_simulation_result(result_graph); + + std::vector county_ids(result_graph.nodes().size()); + std::transform(result_graph.nodes().begin(), result_graph.nodes().end(), county_ids.begin(), [](auto& n) { + return n.id; + }); + + auto save_result_status = save_result(result, county_ids, 1, "graph_result.h5"); + + for (auto&& node : result_graph.nodes()) { + node.property.get_result().print_table(); + } + return 0; } diff --git a/cpp/examples/graph_extended.cpp b/cpp/examples/graph_extended.cpp new file mode 100644 index 0000000000..5fde967f6d --- /dev/null +++ b/cpp/examples/graph_extended.cpp @@ -0,0 +1,263 @@ + +#include "ode_seir/model.h" +#include "ode_seir/infection_state.h" +#include "ode_seir/parameters.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/compartments/simulation.h" +#include "memilio/io/result_io.h" +#include "memilio/io/epi_data.h" + +/** + * indices of contact matrix corresponding to locations where contacts occur. + */ +enum class ContactLocation +{ + Home = 0, + School, + Work, + Other, + Count, +}; + +static const std::map contact_locations = {{ContactLocation::Home, "home"}, + {ContactLocation::School, "school_pf_eig"}, + {ContactLocation::Work, "work"}, + {ContactLocation::Other, "other"}}; + +/** + * Set contact matrices. + * Reads contact matrices from files in the data directory. + * @param data_dir data directory. + * @param params Object that the contact matrices will be added to. + * @returns any io errors that happen during reading of the files. + */ +mio::IOResult set_contact_matrices(const fs::path& data_dir, mio::oseir::Parameters& params, + bool synthetic_population) +{ + if (!synthetic_population) { + //TODO: io error handling + auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), size_t(params.get_num_groups())); + for (auto&& contact_location : contact_locations) { + BOOST_OUTCOME_TRY(auto&& baseline, + mio::read_mobility_plain( + (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + BOOST_OUTCOME_TRY(auto&& minimum, + mio::read_mobility_plain( + (data_dir / "contacts" / ("minimum_" + contact_location.second + ".txt")).string())); + contact_matrices[size_t(contact_location.first)].get_baseline() = baseline; + contact_matrices[size_t(contact_location.first)].get_minimum() = minimum; + } + params.get>() = mio::UncertainContactMatrix(contact_matrices); + } + else { + mio::ContactMatrixGroup& contact_matrix = params.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(7.95 / (size_t)params.get_num_groups()); + } + + printf("Setting contact matrices successful.\n"); + return mio::success(); +} + +/** + * Set epidemiological parameters of Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +mio::IOResult set_covid_parameters(mio::oseir::Parameters& params, bool synthetic_population) +{ + params.template set>(3.335); + if (!synthetic_population) { + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + } + else { + params.template set>(8.097612257); + + params.template set>(0.07333); + } + + printf("Setting epidemiological parameters successful.\n"); + return mio::success(); +} + +mio::IOResult>> +set_population_data(const fs::path& data_dir, mio::oseir::Parameters& params, std::vector node_ids) +{ + size_t number_regions = node_ids.size(); + + std::vector> nodes(number_regions, + mio::oseir::Model(int(size_t(params.get_num_groups())))); + + for (auto& node : nodes) { + node.parameters = params; + } + + BOOST_OUTCOME_TRY(const auto&& population_data, + mio::read_population_data( + (data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true)); + + std::vector> vnum_population(node_ids.size(), + std::vector((size_t)params.get_num_groups(), 0.0)); + + for (auto&& entry : population_data) { + auto it = std::find_if(node_ids.begin(), node_ids.end(), [&entry](auto r) { + return r == 0 || + (entry.county_id && mio::regions::StateId(r) == mio::regions::get_state_id(int(*entry.county_id))) || + (entry.county_id && mio::regions::CountyId(r) == *entry.county_id) || + (entry.district_id && mio::regions::DistrictId(r) == *entry.district_id); + }); + if (it != node_ids.end()) { + auto region_idx = size_t(it - node_ids.begin()); + auto& num_population = vnum_population[region_idx]; + for (size_t age = 0; age < num_population.size(); age++) { + num_population[age] += entry.population[mio::AgeGroup(age)]; + } + } + } + + for (size_t region = 0; region < node_ids.size(); region++) { + auto num_groups = nodes[region].parameters.get_num_groups(); + for (auto i = mio::AgeGroup(0); i < num_groups; i++) { + nodes[region].populations.template set_difference_from_group_total( + {i, mio::oseir::InfectionState::Susceptible}, vnum_population[region][size_t(i)]); + } + } + nodes[27].populations[{mio::AgeGroup(4), mio::oseir::InfectionState::Susceptible}] -= 100; + nodes[27].populations[{mio::AgeGroup(4), mio::oseir::InfectionState::Exposed}] += 100; + + return mio::success(nodes); +} + +mio::IOResult>> +set_synthetic_population_data(mio::oseir::Parameters& params) +{ + size_t number_regions = 53; + + std::vector> nodes(number_regions, + mio::oseir::Model(int(size_t(params.get_num_groups())))); + + mio::Populations population( + {params.get_num_groups(), mio::oseir::InfectionState::Count}); + + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + population[{i, mio::oseir::InfectionState::Susceptible}] = 1000000; + } + for (auto& node : nodes) { + node.parameters = params; + node.populations = population; + } + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + nodes[0].populations[{i, mio::oseir::InfectionState::Exposed}] = 100; + nodes[0].populations[{i, mio::oseir::InfectionState::Susceptible}] = 999900; + } + + return mio::success(nodes); +} + +mio::IOResult run(const fs::path& data_dir, double t0, double tmax, double dt) +{ + mio::set_log_level(mio::LogLevel::off); + // global parameters + bool synthetic_population = false; + const int num_age_groups = 6; + if (num_age_groups != 6) { + synthetic_population = true; + } + mio::oseir::Parameters params(num_age_groups); + + BOOST_OUTCOME_TRY(set_covid_parameters(params, synthetic_population)); + + // set contact matrix + BOOST_OUTCOME_TRY(set_contact_matrices(data_dir, params, synthetic_population)); + + // graph of counties with populations and local parameters + // and mobility between counties + mio::Graph>>, mio::MobilityEdge<>> params_graph; + + BOOST_OUTCOME_TRY( + auto&& node_ids, + mio::get_node_ids((data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true, + true)); + + if (synthetic_population) { + BOOST_OUTCOME_TRY(auto&& nodes, set_synthetic_population_data(params)); + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + params_graph.add_node(node_ids[node_idx], nodes[node_idx]); + } + printf("Setting synthetic population successful.\n"); + } + else { + BOOST_OUTCOME_TRY(auto&& nodes, set_population_data(data_dir, params, node_ids)); + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + params_graph.add_node(node_ids[node_idx], nodes[node_idx]); + } + printf("Setting population from data successful.\n"); + } + + BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, + mio::read_mobility_plain((data_dir / "mobility" / "commuter_mobility_nrw.txt").string())); + if (mobility_data_commuter.rows() != Eigen::Index(params_graph.nodes().size()) || + mobility_data_commuter.cols() != Eigen::Index(params_graph.nodes().size())) { + return mio::failure(mio::StatusCode::InvalidValue, + "Mobility matrices do not have the correct size. You may need to run " + "transformMobilitydata.py from pycode memilio epidata package."); + } + + for (size_t county_idx_i = 0; county_idx_i < params_graph.nodes().size(); ++county_idx_i) { + for (size_t county_idx_j = 0; county_idx_j < params_graph.nodes().size(); ++county_idx_j) { + auto& populations = params_graph.nodes()[county_idx_i].property.get_simulation().get_model().populations; + + auto commuter_coeff_ij = mobility_data_commuter(county_idx_i, county_idx_j) / populations.get_total(); + params_graph.add_edge( + county_idx_i, county_idx_j, + Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count * size_t(params.get_num_groups()), + commuter_coeff_ij)); + } + } + + // for (auto& node : params_graph.nodes()) { + // node.property.get_simulation().set_integrator(std::make_shared>()); + // } + + auto sim = mio::make_mobility_sim(t0, dt, std::move(params_graph)); + + printf("Start Simulation\n"); + sim.advance(tmax); + + auto result_graph = std::move(sim).get_graph(); + auto result = mio::interpolate_simulation_result(result_graph); + + std::vector county_ids(result_graph.nodes().size()); + std::transform(result_graph.nodes().begin(), result_graph.nodes().end(), county_ids.begin(), [](auto& n) { + return n.id; + }); + + auto save_result_status = save_result(result, county_ids, num_age_groups, "graph_result_nrw_adaptive.h5"); + + return mio::success(); +} + +int main() +{ + const auto t0 = 0.; + const auto tmax = 100.; + const auto dt = 0.5; //time step of mobility, daily mobility every second step + + const std::string& data_dir = ""; + + auto result = run(data_dir, t0, tmax, dt); + + return 0; +} diff --git a/cpp/examples/graph_steps.cpp b/cpp/examples/graph_steps.cpp new file mode 100644 index 0000000000..d1b2479d64 --- /dev/null +++ b/cpp/examples/graph_steps.cpp @@ -0,0 +1,163 @@ + +#include "ode_seir/model.h" +#include "ode_seir/infection_state.h" +#include "ode_seir/parameters.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/compartments/simulation.h" + +#include + +bool age_groups = true; + +void set_contact_matrices(mio::oseir::Parameters& params) +{ + if (age_groups) { + Eigen::MatrixXd contact_matrix_eigen(6, 6); + contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, + 0.0161, 0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, + 0.0458, 0.1939, 0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724; + mio::ContactMatrixGroup& contact_matrix = + params.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline() = contact_matrix_eigen; + } + else { + mio::ContactMatrixGroup& contact_matrix = params.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(7.95); + } +} + +/** + * Set epidemiological parameters of Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +void set_covid_parameters(mio::oseir::Parameters& params) +{ + params.template set>(3.335); + + if (age_groups) { + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + } + else { + params.get>()[mio::AgeGroup(0)] = 0.07333; + params.get>()[mio::AgeGroup(0)] = 8.097612257; + } +} + +void set_population_data(mio::oseir::Parameters& params, + mio::Graph>>, + mio::MobilityEdge<>>& params_graph, + size_t number_regions) +{ + std::vector> nodes(number_regions, + mio::oseir::Model(int(size_t(params.get_num_groups())))); + + mio::Populations population( + {params.get_num_groups(), mio::oseir::InfectionState::Count}); + + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + population[{i, mio::oseir::InfectionState::Susceptible}] = 10000; + } + for (auto& node : nodes) { + node.parameters = params; + node.populations = population; + } + // for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + nodes[0].populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] += 100; + nodes[0].populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] -= 100; + // } + + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + params_graph.add_node(node_idx, nodes[node_idx]); + } +} + +void set_parameters_and_population(mio::Graph>>, + mio::MobilityEdge<>>& params_graph, + size_t number_regions) +{ + int num_age_groups = 1; + if (age_groups) { + num_age_groups = 6; + } + + mio::oseir::Parameters params(num_age_groups); + + set_covid_parameters(params); + + // set contact matrix + set_contact_matrices(params); + + set_population_data(params, params_graph, number_regions); + + for (size_t county_idx_i = 0; county_idx_i < params_graph.nodes().size(); ++county_idx_i) { + for (size_t county_idx_j = 0; county_idx_j < params_graph.nodes().size(); ++county_idx_j) { + double commuter_coeff_ij = 1. / (2 * number_regions); + if (county_idx_i == county_idx_j) { + commuter_coeff_ij = 0; + } + params_graph.add_edge( + county_idx_i, county_idx_j, + Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count * size_t(params.get_num_groups()), + commuter_coeff_ij)); + } + } +} + +void simulate(ScalarType tol, ScalarType tmax) +{ + ScalarType t0 = 0.; + ScalarType dt = 0.5; + size_t number_regions = 100; + + mio::Graph>>, mio::MobilityEdge<>> params_graph; + + set_parameters_and_population(params_graph, number_regions); + + using DefaultIntegratorCore = + mio::ControlledStepperWrapper; + + for (auto& node : params_graph.nodes()) { + node.property.get_simulation().set_integrator(std::make_shared(tol)); + } + + auto sim = mio::make_mobility_sim(t0, dt, std::move(params_graph)); + sim.advance(tmax); + + auto result_graph = std::move(sim).get_graph(); + + for (auto&& node : result_graph.nodes()) { + std::cout << " \"Steps Region " << node.id << "\": " << node.property.get_result().get_num_time_points() - 1 + << "," << std::endl; + } +} + +int main(int argc, char** argv) +{ + mio::set_log_level(mio::LogLevel::off); + const ScalarType tmax = 10; + ScalarType tol = 1e-2; + + if (argc > 1) { + tol = std::stod(argv[1]); + } + + std::cout << "{ \"Absolute tolerance\": " << tol << ", " << std::endl; + simulate(tol, tmax); + std::cout << "\n}," << std::endl; + + return 0; +} diff --git a/cpp/examples/graph_timing.cpp b/cpp/examples/graph_timing.cpp new file mode 100644 index 0000000000..7e7b9c2ba3 --- /dev/null +++ b/cpp/examples/graph_timing.cpp @@ -0,0 +1,182 @@ + +#include "ode_seir/model.h" +#include "ode_seir/infection_state.h" +#include "ode_seir/parameters.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/compartments/simulation.h" + +#include + +bool age_groups = true; + +void set_contact_matrices(mio::oseir::Parameters& params) +{ + if (age_groups) { + Eigen::MatrixXd contact_matrix_eigen(6, 6); + contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, + 0.0161, 0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, + 0.0458, 0.1939, 0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724; + mio::ContactMatrixGroup& contact_matrix = + params.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline() = contact_matrix_eigen; + } + else { + mio::ContactMatrixGroup& contact_matrix = params.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(7.95); + } +} + +/** + * Set epidemiological parameters of Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +void set_covid_parameters(mio::oseir::Parameters& params) +{ + params.template set>(3.335); + + if (age_groups) { + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + } + else { + params.get>()[mio::AgeGroup(0)] = 0.07333; + params.get>()[mio::AgeGroup(0)] = 8.097612257; + } +} + +void set_population_data(mio::oseir::Parameters& params, + mio::Graph>>, + mio::MobilityEdge<>>& params_graph, + size_t number_regions) +{ + std::vector> nodes(number_regions, + mio::oseir::Model(int(size_t(params.get_num_groups())))); + + mio::Populations population( + {params.get_num_groups(), mio::oseir::InfectionState::Count}); + + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + population[{i, mio::oseir::InfectionState::Susceptible}] = 10000; + } + for (auto& node : nodes) { + node.parameters = params; + node.populations = population; + } + // for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + nodes[0].populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] += 100; + nodes[0].populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] -= 100; + // } + + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + params_graph.add_node(node_idx, nodes[node_idx]); + } +} + +void set_parameters_and_population(mio::Graph>>, + mio::MobilityEdge<>>& params_graph, + size_t number_regions) +{ + int num_age_groups = 1; + if (age_groups) { + num_age_groups = 6; + } + + mio::oseir::Parameters params(num_age_groups); + + set_covid_parameters(params); + + // set contact matrix + set_contact_matrices(params); + + set_population_data(params, params_graph, number_regions); + + for (size_t county_idx_i = 0; county_idx_i < params_graph.nodes().size(); ++county_idx_i) { + for (size_t county_idx_j = 0; county_idx_j < params_graph.nodes().size(); ++county_idx_j) { + double commuter_coeff_ij = 1. / (2 * number_regions); + if (county_idx_i == county_idx_j) { + commuter_coeff_ij = 0; + } + params_graph.add_edge( + county_idx_i, county_idx_j, + Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count * size_t(params.get_num_groups()), + commuter_coeff_ij)); + } + } + + for (auto& node : params_graph.nodes()) { + node.property.get_simulation().set_integrator(std::make_shared>()); + } +} + +double simulate_runtime(size_t number_regions, ScalarType tmax) +{ + ScalarType t0 = 0.; + ScalarType dt = 0.5; + + mio::Graph>>, mio::MobilityEdge<>> params_graph; + + set_parameters_and_population(params_graph, number_regions); + + auto sim = mio::make_mobility_sim(t0, dt, std::move(params_graph)); + auto start_time = omp_get_wtime(); + sim.advance(tmax); + auto end_time = omp_get_wtime(); + + return end_time - start_time; +} + +void simulate(size_t number_regions, ScalarType tmax) +{ + ScalarType t0 = 0.; + ScalarType dt = 0.5; + + mio::Graph>>, mio::MobilityEdge<>> params_graph; + + set_parameters_and_population(params_graph, number_regions); + + auto sim = mio::make_mobility_sim(t0, dt, std::move(params_graph)); + sim.advance(tmax); +} + +int main(int argc, char** argv) +{ + mio::set_log_level(mio::LogLevel::off); + const ScalarType tmax = 20; + size_t warm_up = 10; + size_t num_runs = 100; + size_t num_regions = 10; + if (argc > 3) { + warm_up = std::stod(argv[1]); + num_runs = std::stod(argv[2]); + num_regions = std::stod(argv[3]); + } + + std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl; + // Warm up runs. + for (size_t i = 0; i < warm_up; i++) { + simulate(num_regions, tmax); + } + + // Runs with timing. + ScalarType total = 0; + for (size_t i = 0; i < num_runs; i++) { + double run_time = simulate_runtime(num_regions, tmax); + total += run_time; + } + std::cout << "\"Time\": " << total / num_runs << "\n}," << std::endl; + + return 0; +} diff --git a/cpp/examples/likwid_test.cpp b/cpp/examples/likwid_test.cpp new file mode 100644 index 0000000000..d2f26a0b56 --- /dev/null +++ b/cpp/examples/likwid_test.cpp @@ -0,0 +1,20 @@ +#include + +double work(double* a, size_t n) { + double s = 0; + for (size_t j=0; j set_covid_parameters(mio::oseir::Parameters& params, bool synthetic_population) +{ + params.template set>(3.335); + if (!synthetic_population) { + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + } + else { + params.template set>(8.097612257); + + params.template set>(0.07333); + } + + printf("Setting epidemiological parameters successful.\n"); + return mio::success(); +} + +/** + * indices of contact matrix corresponding to locations where contacts occur. + */ +enum class ContactLocation +{ + Home = 0, + School, + Work, + Other, + Count, +}; + +static const std::map contact_locations = {{ContactLocation::Home, "home"}, + {ContactLocation::School, "school_pf_eig"}, + {ContactLocation::Work, "work"}, + {ContactLocation::Other, "other"}}; + +/** + * Set contact matrices. + * Reads contact matrices from files in the data directory. + * @param data_dir data directory. + * @param params Object that the contact matrices will be added to. + * @returns any io errors that happen during reading of the files. + */ +mio::IOResult set_contact_matrices(const fs::path& data_dir, mio::oseir::Parameters& params, + bool synthetic_population) +{ + if (!synthetic_population) { + //TODO: io error handling + auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), size_t(params.get_num_groups())); + for (auto&& contact_location : contact_locations) { + BOOST_OUTCOME_TRY(auto&& baseline, + mio::read_mobility_plain( + (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + BOOST_OUTCOME_TRY(auto&& minimum, + mio::read_mobility_plain( + (data_dir / "contacts" / ("minimum_" + contact_location.second + ".txt")).string())); + contact_matrices[size_t(contact_location.first)].get_baseline() = baseline; + contact_matrices[size_t(contact_location.first)].get_minimum() = minimum; + } + params.get>() = mio::UncertainContactMatrix(contact_matrices); + } + else { + mio::ContactMatrixGroup& contact_matrix = params.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(7.95 / (size_t)params.get_num_groups()); + } + + printf("Setting contact matrices successful.\n"); + return mio::success(); +} + +template +mio::IOResult set_population_data(mio::oseir::Model& model, const fs::path& data_dir) +{ + BOOST_OUTCOME_TRY( + auto&& node_ids, + mio::get_node_ids((data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true, + true)); + + BOOST_OUTCOME_TRY(const auto&& population_data, + mio::read_population_data( + (data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true)); + + for (auto&& entry : population_data) { + auto it = std::find_if(node_ids.begin(), node_ids.end(), [&entry](auto r) { + return r == 0 || + (entry.county_id && mio::regions::StateId(r) == mio::regions::get_state_id(int(*entry.county_id))) || + (entry.county_id && mio::regions::CountyId(r) == *entry.county_id) || + (entry.district_id && mio::regions::DistrictId(r) == *entry.district_id); + }); + if (it != node_ids.end()) { + for (size_t age = 0; age < (size_t)model.parameters.get_num_groups(); age++) { + model.populations[{mio::AgeGroup(age), mio::oseir::InfectionState::Susceptible}] += + entry.population[mio::AgeGroup(age)]; + } + } + } + + printf("Setting population data successful.\n"); + return mio::success(); +} +template +mio::IOResult set_parameters_and_population(mio::oseir::Model& model, const fs::path& data_dir, + bool synthetic_population) +{ + auto& populations = model.populations; + auto& parameters = model.parameters; + + size_t number_age_groups = (size_t)parameters.get_num_groups(); + + if (synthetic_population) { + printf("Data is not compatible, using synthetic population instead.\n"); + for (size_t j = 0; j < number_age_groups; j++) { + model.populations[{mio::AgeGroup(j), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(j), mio::oseir::InfectionState::Susceptible}] = 999900; + } + } + else { + BOOST_OUTCOME_TRY(set_population_data(model, data_dir)); + populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] -= 100; + populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] += 100; + } + + BOOST_OUTCOME_TRY(set_contact_matrices(data_dir, parameters, synthetic_population)) + + BOOST_OUTCOME_TRY(set_covid_parameters(parameters, synthetic_population)); + + return mio::success(); +} int main() { mio::set_log_level(mio::LogLevel::debug); - ScalarType t0 = 0; + ScalarType t0 = 0.; ScalarType tmax = 50.; - ScalarType dt = 1.0; + ScalarType dt = 0.1; - mio::log_info("Simulating ODE SEIR; t={} ... {} with dt = {}.", t0, tmax, dt); + ScalarType number_age_groups = 6; + bool synthetic_population = false; + if (number_age_groups != 6) { + synthetic_population = true; + } - mio::oseir::Model model(1); + mio::log_info("Simulating ODE SEIR; t={} ... {} with dt = {}.", t0, tmax, dt); - ScalarType total_population = 10000; - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = - total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; + const std::string& data_dir = ""; - model.parameters.set>(5.2); - model.parameters.set>(6); - model.parameters.set>(0.1); + mio::oseir::Model model(number_age_groups); + auto result_prepare_simulation = set_parameters_and_population(model, data_dir, synthetic_population); - mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); - contact_matrix[0].get_baseline().setConstant(2.7); - contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); + std::shared_ptr> integrator = std::make_shared>(); - model.check_constraints(); + auto seir = simulate(t0, tmax, dt, model, integrator); - auto seir = simulate(t0, tmax, dt, model); + auto reproduction_numbers = model.get_reproduction_numbers(seir); + std::cout << "\nbasis reproduction number: " << reproduction_numbers[0] << "\n"; - seir.print_table({"S", "E", "I", "R"}); - std::cout << "\nnumber total: " << seir.get_last_value().sum() << "\n"; + // seir.print_table({"S", "E", "I", "R"}); + // std::cout << "\nnumber total: " << seir.get_last_value().sum() << "\n"; } diff --git a/cpp/examples/ode_seir_mobility.cpp b/cpp/examples/ode_seir_mobility.cpp new file mode 100644 index 0000000000..4b5600a215 --- /dev/null +++ b/cpp/examples/ode_seir_mobility.cpp @@ -0,0 +1,202 @@ + +#include "memilio/compartments/simulation.h" +#include "memilio/math/euler.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/io/mobility_io.h" +#include "models/ode_seir_mobility/infection_state.h" +#include "models/ode_seir_mobility/model.h" +#include "models/ode_seir_mobility/parameters.h" +#include "models/ode_seir_mobility/regions.h" +#include "memilio/io/io.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/result_io.h" + +/** + * Set epidemiological parameters of Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +mio::IOResult set_covid_parameters(mio::oseirmobility::Parameters& params) +{ + params.template set>(3.335); + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + + printf("Setting epidemiological parameters successful.\n"); + return mio::success(); +} + +/** + * indices of contact matrix corresponding to locations where contacts occur. + */ +enum class ContactLocation +{ + Home = 0, + School, + Work, + Other, + Count, +}; + +static const std::map contact_locations = {{ContactLocation::Home, "home"}, + {ContactLocation::School, "school_pf_eig"}, + {ContactLocation::Work, "work"}, + {ContactLocation::Other, "other"}}; + +/** + * Set contact matrices. + * Reads contact matrices from files in the data directory. + * @param data_dir data directory. + * @param params Object that the contact matrices will be added to. + * @returns any io errors that happen during reading of the files. + */ +mio::IOResult set_contact_matrices(const fs::path& data_dir, mio::oseirmobility::Parameters& params) +{ + //TODO: io error handling + auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), size_t(params.get_num_agegroups())); + for (auto&& contact_location : contact_locations) { + BOOST_OUTCOME_TRY(auto&& baseline, + mio::read_mobility_plain( + (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + BOOST_OUTCOME_TRY(auto&& minimum, + mio::read_mobility_plain( + (data_dir / "contacts" / ("minimum_" + contact_location.second + ".txt")).string())); + contact_matrices[size_t(contact_location.first)].get_baseline() = baseline; + contact_matrices[size_t(contact_location.first)].get_minimum() = minimum; + } + params.get>() = mio::UncertainContactMatrix(contact_matrices); + + printf("Setting contact matrices successful.\n"); + return mio::success(); +} + +template +mio::IOResult set_population_data(mio::oseirmobility::Model& model, const fs::path& data_dir) +{ + BOOST_OUTCOME_TRY( + auto&& node_ids, + mio::get_node_ids((data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true, + true)); + + BOOST_OUTCOME_TRY(const auto&& population_data, + mio::read_population_data( + (data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true)); + + for (auto&& entry : population_data) { + auto it = std::find_if(node_ids.begin(), node_ids.end(), [&entry](auto r) { + return r == 0 || + (entry.county_id && mio::regions::StateId(r) == mio::regions::get_state_id(int(*entry.county_id))) || + (entry.county_id && mio::regions::CountyId(r) == *entry.county_id) || + (entry.district_id && mio::regions::DistrictId(r) == *entry.district_id); + }); + if (it != node_ids.end()) { + auto region_idx = size_t(it - node_ids.begin()); + for (size_t age = 0; age < (size_t)model.parameters.get_num_agegroups(); age++) { + model.populations[{mio::oseirmobility::Region(region_idx), mio::AgeGroup(age), + mio::oseirmobility::InfectionState::Susceptible}] = + entry.population[mio::AgeGroup(age)]; + } + } + } + + printf("Setting population data successful.\n"); + return mio::success(); +} + +template +mio::IOResult set_mobility_weights(mio::oseirmobility::Model& model, const fs::path& data_dir) +{ + size_t number_regions = (size_t)model.parameters.get_num_regions(); + if (number_regions == 1) { + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() + .setConstant(1.0); + + return mio::success(); + } + else { + // mobility between nodes + BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, + mio::read_mobility_plain((data_dir / "mobility" / "commuter_mobility_nrw.txt").string())); + if (mobility_data_commuter.rows() != Eigen::Index(number_regions) || + mobility_data_commuter.cols() != Eigen::Index(number_regions)) { + return mio::failure(mio::StatusCode::InvalidValue, + "Mobility matrices do not have the correct size. You may need to run " + "transformMobilitydata.py from pycode memilio epidata package."); + } + + for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) { + auto population_i = model.populations.get_group_total(mio::oseirmobility::Region(county_idx_i)); + mobility_data_commuter.row(county_idx_i) /= population_i; + } + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() = mobility_data_commuter; + + printf("Setting mobility weights successful.\n"); + return mio::success(); + } +} + +template +mio::IOResult set_parameters_and_population(mio::oseirmobility::Model& model, const fs::path& data_dir) +{ + auto& populations = model.populations; + auto& parameters = model.parameters; + + BOOST_OUTCOME_TRY(set_population_data(model, data_dir)); + populations[{mio::oseirmobility::Region(27), mio::AgeGroup(0), mio::oseirmobility::InfectionState::Susceptible}] -= + 100; + populations[{mio::oseirmobility::Region(27), mio::AgeGroup(0), mio::oseirmobility::InfectionState::Exposed}] += 100; + BOOST_OUTCOME_TRY(set_mobility_weights(model, data_dir)); + + BOOST_OUTCOME_TRY(set_contact_matrices(data_dir, parameters)) + + BOOST_OUTCOME_TRY(set_covid_parameters(parameters)); + + return mio::success(); +} + +int main() +{ + mio::set_log_level(mio::LogLevel::debug); + + ScalarType t0 = 0.; + ScalarType tmax = 100.; + ScalarType dt = 0.1; + + ScalarType number_regions = 53; + ScalarType number_age_groups = 6; + + mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt); + + const std::string& data_dir = ""; + + mio::oseirmobility::Model model(number_regions, number_age_groups); + auto result_prepare_simulation = set_parameters_and_population(model, data_dir); + + // std::shared_ptr> integrator = std::make_shared>(); + + model.check_constraints(); + + auto result_from_sim = simulate(t0, tmax, dt, model); + + auto result = mio::interpolate_simulation_result(result_from_sim); + + auto save_result_status = + mio::save_result({result}, {1}, number_regions * number_age_groups, "ode_result_paper_nrw_adaptive.h5"); +} diff --git a/cpp/examples/ode_seir_mobility_improved.cpp b/cpp/examples/ode_seir_mobility_improved.cpp new file mode 100644 index 0000000000..92565b552d --- /dev/null +++ b/cpp/examples/ode_seir_mobility_improved.cpp @@ -0,0 +1,234 @@ + +#include "memilio/compartments/simulation.h" +#include "memilio/math/euler.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/io/mobility_io.h" +#include "models/ode_seir_mobility_improved/infection_state.h" +#include "models/ode_seir_mobility_improved/model.h" +#include "models/ode_seir_mobility_improved/parameters.h" +#include "models/ode_seir_mobility_improved/regions.h" +#include "memilio/io/io.h" +#include "memilio/io/result_io.h" +#include "memilio/io/epi_data.h" + +/** + * Set epidemiological parameters of Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +mio::IOResult set_covid_parameters(mio::oseirmobilityimproved::Parameters& params) +{ + params.template set>(3.335); + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + printf("Setting epidemiological parameters successful.\n"); + return mio::success(); +} + +/** + * indices of contact matrix corresponding to locations where contacts occur. + */ +enum class ContactLocation +{ + Home = 0, + School, + Work, + Other, + Count, +}; + +static const std::map contact_locations = {{ContactLocation::Home, "home"}, + {ContactLocation::School, "school_pf_eig"}, + {ContactLocation::Work, "work"}, + {ContactLocation::Other, "other"}}; + +/** + * Set contact matrices. + * Reads contact matrices from files in the data directory. + * @param data_dir data directory. + * @param params Object that the contact matrices will be added to. + * @returns any io errors that happen during reading of the files. + */ +mio::IOResult set_contact_matrices(const fs::path& data_dir, + mio::oseirmobilityimproved::Parameters& params) +{ + //TODO: io error handling + auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), size_t(params.get_num_agegroups())); + for (auto&& contact_location : contact_locations) { + BOOST_OUTCOME_TRY(auto&& baseline, + mio::read_mobility_plain( + (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + BOOST_OUTCOME_TRY(auto&& minimum, + mio::read_mobility_plain( + (data_dir / "contacts" / ("minimum_" + contact_location.second + ".txt")).string())); + contact_matrices[size_t(contact_location.first)].get_baseline() = baseline; + contact_matrices[size_t(contact_location.first)].get_minimum() = minimum; + } + params.get>() = + mio::UncertainContactMatrix(contact_matrices); + + printf("Setting contact matrices successful.\n"); + return mio::success(); +} + +template +mio::IOResult set_population_data(mio::oseirmobilityimproved::Model& model, const fs::path& data_dir) +{ + BOOST_OUTCOME_TRY( + auto&& node_ids, + mio::get_node_ids((data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true, + true)); + + BOOST_OUTCOME_TRY(const auto&& population_data, + mio::read_population_data( + (data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true)); + + for (auto&& entry : population_data) { + auto it = std::find_if(node_ids.begin(), node_ids.end(), [&entry](auto r) { + return r == 0 || + (entry.county_id && mio::regions::StateId(r) == mio::regions::get_state_id(int(*entry.county_id))) || + (entry.county_id && mio::regions::CountyId(r) == *entry.county_id) || + (entry.district_id && mio::regions::DistrictId(r) == *entry.district_id); + }); + if (it != node_ids.end()) { + auto region_idx = size_t(it - node_ids.begin()); + for (size_t age = 0; age < (size_t)model.parameters.get_num_agegroups(); age++) { + model.populations[{mio::oseirmobilityimproved::Region(region_idx), mio::AgeGroup(age), + mio::oseirmobilityimproved::InfectionState::Susceptible}] = + entry.population[mio::AgeGroup(age)]; + } + } + } + + printf("Setting population data successful.\n"); + return mio::success(); +} + +template +mio::IOResult set_mobility_weights(mio::oseirmobilityimproved::Model& model, const fs::path& data_dir) +{ + size_t number_regions = (size_t)model.parameters.get_num_regions(); + if (number_regions == 1) { + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() + .setConstant(1.0); + + return mio::success(); + } + else { + // mobility between nodes + BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, + mio::read_mobility_plain((data_dir / "mobility" / "commuter_mobility_nrw.txt").string())); + if (mobility_data_commuter.rows() != Eigen::Index(number_regions) || + mobility_data_commuter.cols() != Eigen::Index(number_regions)) { + return mio::failure(mio::StatusCode::InvalidValue, + "Mobility matrices do not have the correct size. You may need to run " + "transformMobilitydata.py from pycode memilio epidata package."); + } + + for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) { + auto population_i = model.populations.get_group_total(mio::oseirmobilityimproved::Region(county_idx_i)); + mobility_data_commuter.row(county_idx_i) /= population_i; + mobility_data_commuter(county_idx_i, county_idx_i) = + 1 - mobility_data_commuter.rowwise().sum()(county_idx_i); + } + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() = mobility_data_commuter; + + printf("Setting mobility weights successful.\n"); + return mio::success(); + } +} + +template +mio::IOResult set_parameters_and_population(mio::oseirmobilityimproved::Model& model, + const fs::path& data_dir) +{ + auto& populations = model.populations; + auto& parameters = model.parameters; + + size_t number_regions = (size_t)parameters.get_num_regions(); + size_t number_age_groups = (size_t)parameters.get_num_agegroups(); + + BOOST_OUTCOME_TRY(set_population_data(model, data_dir)); + populations[{mio::oseirmobilityimproved::Region(27), mio::AgeGroup(4), + mio::oseirmobilityimproved::InfectionState::Susceptible}] -= 100; + populations[{mio::oseirmobilityimproved::Region(27), mio::AgeGroup(4), + mio::oseirmobilityimproved::InfectionState::Exposed}] += 100; + + BOOST_OUTCOME_TRY(set_mobility_weights(model, data_dir)); + + BOOST_OUTCOME_TRY(set_contact_matrices(data_dir, parameters)) + + BOOST_OUTCOME_TRY(set_covid_parameters(parameters)); + + mio::ContactMatrixGroup& commuting_strengths = + parameters.template get>().get_cont_freq_mat(); + + auto& population_after_commuting = model.m_population_after_commuting; + for (size_t region_n = 0; region_n < number_regions; ++region_n) { + for (size_t age = 0; age < number_age_groups; ++age) { + double population_n = 0; + for (size_t state = 0; state < (size_t)mio::oseirmobilityimproved::InfectionState::Count; state++) { + population_n += populations[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age), + mio::oseirmobilityimproved::InfectionState(state)}]; + } + population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] += + population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + population_after_commuting[{mio::oseirmobilityimproved::Region(region_m), mio::AgeGroup(age)}] += + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + } + } + } + + return mio::success(); +} + +int main() +{ + mio::set_log_level(mio::LogLevel::debug); + + ScalarType t0 = 0.; + ScalarType tmax = 100.; + ScalarType dt = 0.1; + + ScalarType number_regions = 53; + ScalarType number_age_groups = 6; + + mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt); + + const std::string& data_dir = ""; + + mio::oseirmobilityimproved::Model model(number_regions, number_age_groups); + auto result_prepare_simulation = set_parameters_and_population(model, data_dir); + + // std::shared_ptr> integrator = std::make_shared>(); + + model.check_constraints(); + + printf("Start Simulation\n"); + auto result_from_sim = simulate(t0, tmax, dt, model); + + auto result = mio::interpolate_simulation_result(result_from_sim); + + auto save_result_status = + mio::save_result({result}, {1}, number_regions * number_age_groups, "ode_result_nrw_adaptive_test.h5"); +} diff --git a/cpp/examples/ode_seir_mobility_steps.cpp b/cpp/examples/ode_seir_mobility_steps.cpp new file mode 100644 index 0000000000..6068214be9 --- /dev/null +++ b/cpp/examples/ode_seir_mobility_steps.cpp @@ -0,0 +1,181 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Carlotta Gerstein +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "memilio/compartments/simulation.h" +#include "memilio/math/euler.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_seir_mobility_improved/infection_state.h" +#include "models/ode_seir_mobility_improved/model.h" +#include "models/ode_seir_mobility_improved/parameters.h" +#include "models/ode_seir_mobility_improved/regions.h" + +#include + +bool age_groups = true; + +template +void set_contact_matrix(mio::oseirmobilityimproved::Model& model) +{ + if (age_groups) { + Eigen::MatrixXd contact_matrix_eigen(6, 6); + contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, + 0.0161, 0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, + 0.0458, 0.1939, 0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724; + mio::ContactMatrixGroup& contact_matrix = + model.parameters.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline() = contact_matrix_eigen; + } + { + mio::ContactMatrixGroup& contact_matrix = + model.parameters.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(7.95); + } +} + +/** + * Set epidemiological parameters of Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +void set_covid_parameters(mio::oseirmobilityimproved::Parameters& params) +{ + params.template set>(3.335); + + if (age_groups) { + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + } + else { + params.get>()[mio::AgeGroup(0)] = 0.07333; + params.get>()[mio::AgeGroup(0)] = 8.097612257; + } +} + +template +void set_mobility_weights(mio::oseirmobilityimproved::Model& model) +{ + size_t number_regions = (size_t)model.parameters.get_num_regions(); + double fraction_commuter = 1. / (2 * number_regions); + Eigen::MatrixXd mobility_data_commuter = + Eigen::MatrixXd::Constant(number_regions, number_regions, fraction_commuter) - + fraction_commuter * + Eigen::MatrixXd::Identity(number_regions, number_regions); // Ensure that the diagonal is zero + for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) { + mobility_data_commuter(county_idx_i, county_idx_i) = 1 - mobility_data_commuter.rowwise().sum()(county_idx_i); + } + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() = mobility_data_commuter; +} + +template +void set_parameters_and_population(mio::oseirmobilityimproved::Model& model) +{ + auto& populations = model.populations; + auto& parameters = model.parameters; + + size_t number_regions = (size_t)parameters.get_num_regions(); + size_t number_age_groups = (size_t)parameters.get_num_agegroups(); + for (size_t j = 0; j < number_age_groups; j++) { + for (size_t i = 0; i < number_regions; i++) { + model.populations[{mio::oseirmobilityimproved::Region(i), mio::AgeGroup(j), + mio::oseirmobilityimproved::InfectionState::Susceptible}] = 10000; + } + } + model.populations[{mio::oseirmobilityimproved::Region(0), mio::AgeGroup(0), + mio::oseirmobilityimproved::InfectionState::Exposed}] += 100; + model.populations[{mio::oseirmobilityimproved::Region(0), mio::AgeGroup(0), + mio::oseirmobilityimproved::InfectionState::Susceptible}] -= 100; + set_mobility_weights(model); + + set_contact_matrix(model); + + set_covid_parameters(parameters); + + mio::ContactMatrixGroup& commuting_strengths = + parameters.template get>().get_cont_freq_mat(); + + auto& population_after_commuting = model.m_population_after_commuting; + for (size_t region_n = 0; region_n < number_regions; ++region_n) { + for (size_t age = 0; age < number_age_groups; ++age) { + double population_n = 0; + for (size_t state = 0; state < (size_t)mio::oseirmobilityimproved::InfectionState::Count; state++) { + population_n += populations[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age), + mio::oseirmobilityimproved::InfectionState(state)}]; + } + population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] += + population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + population_after_commuting[{mio::oseirmobilityimproved::Region(region_m), mio::AgeGroup(age)}] += + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + } + } + } +} + +void simulate(ScalarType tol, ScalarType tmax) +{ + mio::set_log_level(mio::LogLevel::off); + ScalarType t0 = 0.; + ScalarType dt = 0.1; + size_t number_regions = 100; + ScalarType number_age_groups = 1; + if (age_groups) { + number_age_groups = 6; + } + + mio::oseirmobilityimproved::Model model(number_regions, number_age_groups); + set_parameters_and_population(model); + using DefaultIntegratorCore = + mio::ControlledStepperWrapper; + + std::shared_ptr> integrator = std::make_shared(tol); + std::cout << "{ \"Absolute tolerance\": " << tol << ", " << std::endl; + + auto result = simulate(t0, tmax, dt, model, integrator); + std::cout << "\"Steps\": " << result.get_num_time_points() - 1 << "," << std::endl; +} + +int main(int argc, char** argv) +{ + const ScalarType tmax = 10; + ScalarType tol = 1e-12; + + if (argc > 1) { + tol = std::stod(argv[1]); + } + + simulate(tol, tmax); + return 0; +} diff --git a/cpp/examples/ode_seir_mobility_timing.cpp b/cpp/examples/ode_seir_mobility_timing.cpp new file mode 100644 index 0000000000..0b3b184941 --- /dev/null +++ b/cpp/examples/ode_seir_mobility_timing.cpp @@ -0,0 +1,192 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Carlotta Gerstein +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "memilio/compartments/simulation.h" +#include "memilio/math/euler.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_seir_mobility_improved/infection_state.h" +#include "models/ode_seir_mobility_improved/model.h" +#include "models/ode_seir_mobility_improved/parameters.h" +#include "models/ode_seir_mobility_improved/regions.h" + +#include + +bool age_groups = true; + +template +void set_contact_matrix(mio::oseirmobilityimproved::Model& model) +{ + if (age_groups) { + Eigen::MatrixXd contact_matrix_eigen(6, 6); + contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, + 0.0161, 0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, + 0.0458, 0.1939, 0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724; + mio::ContactMatrixGroup& contact_matrix = + model.parameters.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline() = contact_matrix_eigen; + } + { + mio::ContactMatrixGroup& contact_matrix = + model.parameters.template get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(7.95); + } +} + +/** + * Set epidemiological parameters of Sars-CoV-2 for a immune-naive + * population and wild type variant. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +void set_covid_parameters(mio::oseirmobilityimproved::Parameters& params) +{ + params.template set>(3.335); + + if (age_groups) { + params.get>()[mio::AgeGroup(0)] = 8.0096875; + params.get>()[mio::AgeGroup(1)] = 8.0096875; + params.get>()[mio::AgeGroup(2)] = 8.2182; + params.get>()[mio::AgeGroup(3)] = 8.1158; + params.get>()[mio::AgeGroup(4)] = 8.033; + params.get>()[mio::AgeGroup(5)] = 7.985; + + params.get>()[mio::AgeGroup(0)] = 0.03; + params.get>()[mio::AgeGroup(1)] = 0.06; + params.get>()[mio::AgeGroup(2)] = 0.06; + params.get>()[mio::AgeGroup(3)] = 0.06; + params.get>()[mio::AgeGroup(4)] = 0.09; + params.get>()[mio::AgeGroup(5)] = 0.175; + } + else { + params.get>()[mio::AgeGroup(0)] = 0.07333; + params.get>()[mio::AgeGroup(0)] = 8.097612257; + } +} + +template +void set_mobility_weights(mio::oseirmobilityimproved::Model& model) +{ + size_t number_regions = (size_t)model.parameters.get_num_regions(); + double fraction_commuter = 1. / (2 * number_regions); + Eigen::MatrixXd mobility_data_commuter = + Eigen::MatrixXd::Constant(number_regions, number_regions, fraction_commuter) - + fraction_commuter * + Eigen::MatrixXd::Identity(number_regions, number_regions); // Ensure that the diagonal is zero + for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) { + mobility_data_commuter(county_idx_i, county_idx_i) = 1 - mobility_data_commuter.rowwise().sum()(county_idx_i); + } + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() = mobility_data_commuter; +} + +template +void set_parameters_and_population(mio::oseirmobilityimproved::Model& model) +{ + auto& populations = model.populations; + auto& parameters = model.parameters; + + size_t number_regions = (size_t)parameters.get_num_regions(); + size_t number_age_groups = (size_t)parameters.get_num_agegroups(); + for (size_t j = 0; j < number_age_groups; j++) { + for (size_t i = 0; i < number_regions; i++) { + populations[{mio::oseirmobilityimproved::Region(i), mio::AgeGroup(j), + mio::oseirmobilityimproved::InfectionState::Susceptible}] = 10000; + } + } + populations[{mio::oseirmobilityimproved::Region(0), mio::AgeGroup(0), + mio::oseirmobilityimproved::InfectionState::Exposed}] += 100; + populations[{mio::oseirmobilityimproved::Region(0), mio::AgeGroup(0), + mio::oseirmobilityimproved::InfectionState::Susceptible}] -= 100; + set_mobility_weights(model); + + set_contact_matrix(model); + + set_covid_parameters(parameters); + + mio::ContactMatrixGroup& commuting_strengths = + parameters.template get>().get_cont_freq_mat(); + + auto& population_after_commuting = model.m_population_after_commuting; + for (size_t region_n = 0; region_n < number_regions; ++region_n) { + for (size_t age = 0; age < number_age_groups; ++age) { + double population_n = 0; + for (size_t state = 0; state < (size_t)mio::oseirmobilityimproved::InfectionState::Count; state++) { + population_n += populations[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age), + mio::oseirmobilityimproved::InfectionState(state)}]; + } + population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] += + population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{mio::oseirmobilityimproved::Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + population_after_commuting[{mio::oseirmobilityimproved::Region(region_m), mio::AgeGroup(age)}] += + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + } + } + } +} + +void simulate(size_t num_warm_up_runs, size_t num_runs, size_t number_regions, ScalarType tmax) +{ + mio::set_log_level(mio::LogLevel::off); + ScalarType t0 = 0.; + ScalarType dt = 0.1; + ScalarType number_age_groups = 1; + if (age_groups) { + number_age_groups = 6; + } + + mio::oseirmobilityimproved::Model model(number_regions, number_age_groups); + set_parameters_and_population(model); + + std::shared_ptr> integrator = std::make_shared>(); + + std::cout << "{ \"Regions\": " << number_regions << ", " << std::endl; + + // Warm up runs. + for (size_t i = 0; i < num_warm_up_runs; i++) { + simulate(t0, tmax, dt, model, integrator); + } + auto result = simulate(t0, tmax, dt, model, integrator); + + // Runs with timing. + ScalarType total = 0; + for (size_t i = 0; i < num_runs; i++) { + double runtime = simulate_runtime(t0, tmax, dt, model, integrator); + total += runtime; + } + std::cout << "\"Time\": " << total / num_runs << "\n}," << std::endl; +} + +int main(int argc, char** argv) +{ + const ScalarType tmax = 20; + size_t warm_up = 10; + size_t num_runs = 100; + size_t num_regions = 10; + if (argc > 3) { + warm_up = std::stod(argv[1]); + num_runs = std::stod(argv[2]); + num_regions = std::stod(argv[3]); + } + simulate(warm_up, num_runs, num_regions, tmax); + return 0; +} diff --git a/cpp/examples/ode_sir_mobility.cpp b/cpp/examples/ode_sir_mobility.cpp new file mode 100644 index 0000000000..7aaf725e45 --- /dev/null +++ b/cpp/examples/ode_sir_mobility.cpp @@ -0,0 +1,241 @@ + +#include "memilio/compartments/simulation.h" +#include "memilio/math/euler.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/io/mobility_io.h" +#include "ode_sir_mobility/infection_state.h" +#include "ode_sir_mobility/model.h" +#include "ode_sir_mobility/parameters.h" +#include "ode_sir_mobility/regions.h" +#include "memilio/io/io.h" + +mio::IOResult>>> read_path_mobility(const std::string& filename) +{ + BOOST_OUTCOME_TRY(auto&& num_lines, mio::count_lines(filename)); + + if (num_lines == 0) { + std::vector>> arr(0, std::vector>(0)); + return mio::success(arr); + } + + std::fstream file; + file.open(filename, std::ios::in); + if (!file.is_open()) { + return failure(mio::StatusCode::FileNotFound, filename); + } + + std::vector>> arr(std::sqrt(num_lines), + std::vector>(std::sqrt(num_lines))); + + try { + std::string tp; + while (getline(file, tp)) { + auto line = mio::split(tp, ' '); + int indx_x = std::stoi(line[0]); + int indx_y = std::stoi(line[1]); + if (indx_x != indx_y) { + auto path = std::accumulate(line.begin() + 2, line.end(), std::string("")); + + // string -> vector of integers + std::vector path_vec; + + // Remove the square brackets and \r + path = path.substr(1, path.size() - 3); + std::stringstream ss(path); + std::string token; + + // get numbers and save them in path_vec + while (std::getline(ss, token, ',')) { + path_vec.push_back(std::stoi(token)); + } + + // Sorted by end location + for (int number : path_vec) { + if (number != indx_x && number != indx_y) { + arr[indx_x][indx_y].push_back(number); + } + } + } + } + } + catch (std::runtime_error& ex) { + return failure(mio::StatusCode::InvalidFileFormat, filename + ": " + ex.what()); + } + + return mio::success(arr); +} + +template +mio::IOResult preprocess(const std::string& filename, mio::osirmobility::Model& model) +{ + BOOST_OUTCOME_TRY(auto&& mobility_paths, read_path_mobility(filename)); + size_t n_regions = (size_t)model.parameters.get_num_regions(); + for (size_t i = 0; i < n_regions; i++) { + for (size_t j = 0; j < n_regions; j++) { + if (j == i) { + continue; + } + std::sort(mobility_paths[i][j].begin(), mobility_paths[i][j].end()); + std::vector intersection_int; + std::vector intersection_region(intersection_int.size(), + mio::osirmobility::Region(0)); + for (size_t k = 0; k < n_regions; k++) { + if (k == i || k == j) { + continue; + } + std::sort(mobility_paths[k][j].begin(), mobility_paths[k][j].end()); + std::set_intersection(mobility_paths[i][j].begin(), mobility_paths[i][j].end(), + mobility_paths[k][j].begin(), mobility_paths[k][j].end(), + std::back_inserter(intersection_int)); + + if (intersection_int.begin() != intersection_int.end()) { + intersection_region.push_back(mio::osirmobility::Region(k)); + intersection_int.pop_back(); + } + } + if (intersection_region.begin() != intersection_region.end()) { + model.parameters.template get()[{ + mio::osirmobility::Region(i), mio::osirmobility::Region(j)}] = intersection_region; + } + } + } + return mio::success(); +} + +template +mio::IOResult set_mobility_weights(const std::string& mobility_data, const std::string& trip_chains, + mio::osirmobility::Model& model, size_t number_regions) +{ + BOOST_OUTCOME_TRY(preprocess(trip_chains, model)); + // mobility between nodes + BOOST_OUTCOME_TRY(auto&& mobility_data_commuter, + mio::read_mobility_plain(mobility_data + "mobility" + "commuter_migration_scaled.txt")); + if (mobility_data_commuter.rows() != Eigen::Index(number_regions) || + mobility_data_commuter.cols() != Eigen::Index(number_regions)) { + return mio::failure(mio::StatusCode::InvalidValue, + "Mobility matrices do not have the correct size. You may need to run " + "transformMobilitydata.py from pycode memilio epidata package."); + } + + for (auto age = mio::AgeGroup(0); age < model.parameters.get_num_agegroups(); age++) { + for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) { + for (size_t county_idx_j = 0; county_idx_j < number_regions; ++county_idx_j) { + //commuters + auto population_i = model.populations.get_group_total(mio::osirmobility::Region(county_idx_i)); + auto commuter_coeff_ij = mobility_data_commuter(county_idx_i, county_idx_j) / population_i; + if (commuter_coeff_ij > 4e-5) { + model.parameters.template get().push_back( + {mio::osirmobility::Region(county_idx_i), mio::osirmobility::Region(county_idx_j), + commuter_coeff_ij}); + } + } + } + } + return mio::success(); +} + +int main() +{ + mio::set_log_level(mio::LogLevel::debug); + + ScalarType t0 = 0.; + ScalarType tmax = 50.; + ScalarType dt = 1; + + ScalarType number_regions = 4; + ScalarType number_age_groups = 1; + ScalarType total_population_per_region = 10; + + mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt); + + const std::string& mobility_data = ""; + const std::string& trip_chain_data = ""; + + mio::osirmobility::Model model(number_regions, number_age_groups); + + for (size_t i = 0; i < number_regions; i++) { + model.populations[{mio::osirmobility::Region(i), mio::AgeGroup(0), + mio::osirmobility::InfectionState::Infected}] = 1; + model.populations[{mio::osirmobility::Region(i), mio::AgeGroup(0), + mio::osirmobility::InfectionState::Recovered}] = 0; + model.populations[{mio::osirmobility::Region(i), mio::AgeGroup(0), + mio::osirmobility::InfectionState::Susceptible}] = + total_population_per_region - + model.populations[{mio::osirmobility::Region(i), mio::AgeGroup(0), + mio::osirmobility::InfectionState::Infected}] - + model.populations[{mio::osirmobility::Region(i), mio::AgeGroup(0), + mio::osirmobility::InfectionState::Recovered}]; + } + + model.parameters.set>(2); + model.parameters.set>(0.04); + model.parameters.set>(1.); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(1.0); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(0), 0.2}); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(2), 0.6}); + model.parameters.get().push_back( + {mio::osirmobility::Region(2), mio::osirmobility::Region(0), 0.5}); + model.parameters.get().push_back( + {mio::osirmobility::Region(0), mio::osirmobility::Region(3), 1.0}); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(3), 0.2}); + + model.parameters.get()[{mio::osirmobility::Region(0), + mio::osirmobility::Region(1)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(0), + mio::osirmobility::Region(3)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(1), + mio::osirmobility::Region(0)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(1), + mio::osirmobility::Region(2)}] = {0}; + model.parameters.get()[{mio::osirmobility::Region(1), + mio::osirmobility::Region(3)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(2), + mio::osirmobility::Region(1)}] = {0}; + model.parameters.get()[{mio::osirmobility::Region(3), + mio::osirmobility::Region(0)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(3), + mio::osirmobility::Region(1)}] = {2}; + + // auto result_preprocess = set_mobility_weights(mobility_data, trip_chain_data, model, number_regions); + + std::shared_ptr> integrator = + std::make_shared>(); + + model.check_constraints(); + + auto sir = simulate(t0, tmax, dt, model, integrator); + + bool print_to_terminal = true; + + sir.print_table(); + + if (print_to_terminal) { + + std::vector vars = {"S", "I", "R"}; + printf("\n # t"); + for (size_t i = 0; i < (size_t)model.parameters.get_num_regions(); i++) { + for (size_t k = 0; k < (size_t)mio::osirmobility::InfectionState::Count; k++) { + printf(" %s_%d", vars[k].c_str(), (int)i); + } + } + + auto num_points = static_cast(sir.get_num_time_points()); + for (size_t i = 0; i < num_points; i++) { + printf("\n%.14f ", sir.get_time(i)); + for (size_t k = 0; k < (size_t)model.parameters.get_num_regions(); k++) { + for (size_t j = 0; j < (size_t)mio::osirmobility::InfectionState::Count; j++) { + printf(" %.14f", sir.get_value(i)[j + (size_t)mio::osirmobility::InfectionState::Count * (int)k]); + } + } + } + printf("\n"); + } +} diff --git a/cpp/memilio/compartments/simulation.h b/cpp/memilio/compartments/simulation.h index 8b3307f87f..b887744430 100644 --- a/cpp/memilio/compartments/simulation.h +++ b/cpp/memilio/compartments/simulation.h @@ -26,6 +26,8 @@ #include "memilio/math/stepper_wrapper.h" #include "memilio/utils/time_series.h" +#include + namespace mio { @@ -225,6 +227,22 @@ TimeSeries simulate(FP t0, FP tmax, FP dt, Model const& model, return sim.get_result(); } +/*Same function, used for timing*/ +template > +double simulate_runtime(FP t0, FP tmax, FP dt, Model const& model, + std::shared_ptr> integrator = nullptr) +{ + model.check_constraints(); + Sim sim(model, t0, dt); + if (integrator) { + sim.set_integrator(integrator); + } + double start_time = omp_get_wtime(); + sim.advance(tmax); + double end_time = omp_get_wtime(); + return end_time - start_time; +} + } // namespace mio #endif // SIMULATION_H diff --git a/cpp/models/ode_seir/model.h b/cpp/models/ode_seir/model.h index fa3bf4db93..3614ebc25e 100644 --- a/cpp/models/ode_seir/model.h +++ b/cpp/models/ode_seir/model.h @@ -80,7 +80,7 @@ class Model const size_t Ii = this->populations.get_flat_index({i, InfectionState::Infected}); for (auto j : make_index_range(age_groups)) { - const size_t Sj = this->populations.get_flat_index({i, InfectionState::Susceptible}); + const size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible}); const size_t Ej = this->populations.get_flat_index({j, InfectionState::Exposed}); const size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected}); const size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered}); diff --git a/cpp/models/ode_seir_mobility/CMakeLists.txt b/cpp/models/ode_seir_mobility/CMakeLists.txt new file mode 100644 index 0000000000..ecf6e3112b --- /dev/null +++ b/cpp/models/ode_seir_mobility/CMakeLists.txt @@ -0,0 +1,13 @@ +add_library(ode_seir_mobility + infection_state.h + model.h + model.cpp + parameters.h + regions.h +) +target_link_libraries(ode_seir_mobility PUBLIC memilio) +target_include_directories(ode_seir_mobility PUBLIC + $ + $ +) +target_compile_options(ode_seir_mobility PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_seir_mobility/infection_state.h b/cpp/models/ode_seir_mobility/infection_state.h new file mode 100644 index 0000000000..f4207708b3 --- /dev/null +++ b/cpp/models/ode_seir_mobility/infection_state.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITY_INFECTIONSTATE_H +#define ODESEIRMOBILITY_INFECTIONSTATE_H + +namespace mio +{ +namespace oseirmobility +{ + +/** + * @brief The InfectionState enum describes the possible + * categories for the infectious state of persons + */ +enum class InfectionState +{ + Susceptible, + Exposed, + Infected, + Recovered, + Count +}; + +} // namespace oseirmobility +} // namespace mio + +#endif // ODESEIR_INFECTIONSTATE_H diff --git a/cpp/models/ode_seir_mobility/model.cpp b/cpp/models/ode_seir_mobility/model.cpp new file mode 100644 index 0000000000..75494e52d6 --- /dev/null +++ b/cpp/models/ode_seir_mobility/model.cpp @@ -0,0 +1,10 @@ + +#include "ode_seir_mobility/model.h" + +namespace mio +{ +namespace oseirmobility +{ + +} // namespace oseirmobility +} // namespace mio diff --git a/cpp/models/ode_seir_mobility/model.h b/cpp/models/ode_seir_mobility/model.h new file mode 100644 index 0000000000..1f775495ae --- /dev/null +++ b/cpp/models/ode_seir_mobility/model.h @@ -0,0 +1,216 @@ + +#ifndef ODESEIRMOBILITY_MODEL_H +#define ODESEIRMOBILITY_MODEL_H + +#include "memilio/compartments/flow_model.h" +#include "memilio/epidemiology/populations.h" +#include "models/ode_seir_mobility/infection_state.h" +#include "models/ode_seir_mobility/parameters.h" +#include "models/ode_seir_mobility/regions.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/time_series.h" + +GCC_CLANG_DIAGNOSTIC(push) +GCC_CLANG_DIAGNOSTIC(ignored "-Wshadow") +#include +GCC_CLANG_DIAGNOSTIC(pop) + +namespace mio +{ +namespace oseirmobility +{ + +/******************** + * define the model * + ********************/ + +using Flows = TypeList, + Flow, + Flow>; + +template +class Model : public FlowModel, + Parameters, Flows> +{ + + using Base = + FlowModel, Parameters, Flows>; + +public: + using typename Base::ParameterSet; + using typename Base::Populations; + + Model(int num_regions, int num_agegroups) + : Base(Populations({Region(num_regions), AgeGroup(num_agegroups), InfectionState::Count}), + ParameterSet(Region(num_regions), AgeGroup(num_agegroups))) + { + } + + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override + { + const auto& params = this->parameters; + const auto& population = this->populations; + const auto& commuting_strengths = + params.template get>().get_cont_freq_mat().get_matrix_at(t); + const Index n_age_groups = reduce_index>(params.get_num_agegroups()); + const Index n_regions = reduce_index>(params.get_num_regions()); + + for (auto age_i : make_index_range(n_age_groups)) { + for (auto region_n : make_index_range(n_regions)) { + for (auto age_j : make_index_range(n_age_groups)) { + FP flow_SE_helper = 0; + const size_t Sjn = population.get_flat_index({region_n, age_j, InfectionState::Susceptible}); + const size_t Ejn = population.get_flat_index({region_n, age_j, InfectionState::Exposed}); + const size_t Ijn = population.get_flat_index({region_n, age_j, InfectionState::Infected}); + const size_t Rjn = population.get_flat_index({region_n, age_j, InfectionState::Recovered}); + + const double Njn_inv = 1.0 / (pop[Sjn] + pop[Ejn] + pop[Ijn] + pop[Rjn]); + + double coeffStoI = params.template get>().get_cont_freq_mat().get_matrix_at(t)( + age_i.get(), age_j.get()) * + params.template get>()[age_i]; + + for (auto region_m : make_index_range(n_regions)) { + const size_t Sjm = population.get_flat_index({region_m, age_j, InfectionState::Susceptible}); + const size_t Ejm = population.get_flat_index({region_m, age_j, InfectionState::Exposed}); + const size_t Ijm = population.get_flat_index({region_m, age_j, InfectionState::Infected}); + const size_t Rjm = population.get_flat_index({region_m, age_j, InfectionState::Recovered}); + + const double Njm_inv = 1.0 / (pop[Sjm] + pop[Ejm] + pop[Ijm] + pop[Rjm]); + if (region_n == region_m) { + flow_SE_helper += + pop[population.get_flat_index({region_n, age_j, InfectionState::Infected})] * Njn_inv; + continue; + } + flow_SE_helper += (commuting_strengths(region_n.get(), region_m.get()) * Njm_inv + + commuting_strengths(region_m.get(), region_n.get()) * Njn_inv) * + pop[population.get_flat_index({region_m, age_j, InfectionState::Infected})]; + } + flows[Base::template get_flat_flow_index( + {region_n, age_i})] += + flow_SE_helper * coeffStoI * + y[population.get_flat_index({region_n, age_i, InfectionState::Susceptible})]; + } + flows[Base::template get_flat_flow_index( + {region_n, age_i})] = (1.0 / params.template get>()[age_i]) * + y[population.get_flat_index({region_n, age_i, InfectionState::Exposed})]; + flows[Base::template get_flat_flow_index( + {region_n, age_i})] = (1.0 / params.template get>()[age_i]) * + y[population.get_flat_index({region_n, age_i, InfectionState::Infected})]; + } + } + } + + /** + *@brief Computes the reproduction number at a given index time of the Model output obtained by the Simulation. + *@param t_idx The index time at which the reproduction number is computed. + *@param y The TimeSeries obtained from the Model Simulation. + *@returns The computed reproduction number at the provided index time. + */ + IOResult get_reproduction_number(size_t t_idx, const mio::TimeSeries& y) + { + if (!(t_idx < static_cast(y.get_num_time_points()))) { + return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries"); + } + + auto const& params = this->parameters; + auto const& pop = this->populations; + + const size_t num_age_groups = (size_t)params.get_num_agegroups(); + const size_t num_regions = (size_t)params.get_num_regions(); + constexpr size_t num_infected_compartments = 2; + const size_t total_infected_compartments = num_infected_compartments * num_age_groups * num_regions; + + ContactMatrixGroup const& contact_matrix = params.template get>(); + ContactMatrixGroup const& commuting_strengths = params.template get>(); + + Eigen::MatrixXd F = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + Eigen::MatrixXd V = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + + for (auto i = AgeGroup(0); i < AgeGroup(num_age_groups); i++) { + for (auto n = Region(0); n < Region(num_regions); n++) { + size_t Si = pop.get_flat_index({n, i, InfectionState::Susceptible}); + for (auto j = AgeGroup(0); j < AgeGroup(num_age_groups); j++) { + auto const population_region_n = pop.template slice({(size_t)n, 1}); + auto const population_region_age_nj = population_region_n.template slice({(size_t)j, 1}); + auto Njn = std::accumulate(population_region_age_nj.begin(), population_region_age_nj.end(), 0.); + for (auto m = Region(0); m < Region(num_regions); m++) { + auto const population_region_m = pop.template slice({(size_t)m, 1}); + auto const population_region_age_mj = + population_region_m.template slice({(size_t)j, 1}); + auto Njm = + std::accumulate(population_region_age_mj.begin(), population_region_age_mj.end(), 0.); + + if (n == m) { + double coeffStoE = contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) * + params.template get>()[i] / + Njm; + F((size_t)i * num_regions + (size_t)n, + num_age_groups * num_regions + (size_t)j * num_regions + (size_t)m) = + coeffStoE * y.get_value(t_idx)[Si]; + } + else { + double coeffStoE = + contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) * + params.template get>()[i] * + (commuting_strengths.get_matrix_at(y.get_time(t_idx))(n.get(), m.get()) / Njm + + commuting_strengths.get_matrix_at(y.get_time(t_idx))(m.get(), n.get())) / + Njn; + F((size_t)i * num_regions + (size_t)n, + num_age_groups * num_regions + (size_t)j * num_regions + (size_t)m) = + coeffStoE * y.get_value(t_idx)[Si]; + } + } + } + + double T_Ei = params.template get>()[i]; + double T_Ii = params.template get>()[i]; + V((size_t)i * num_regions + (size_t)n, (size_t)i * num_regions + (size_t)n) = 1.0 / T_Ei; + V(num_age_groups * num_regions + (size_t)i * num_regions + (size_t)n, + (size_t)i * num_regions + (size_t)n) = -1.0 / T_Ei; + V(num_age_groups * num_regions + (size_t)i * num_regions + (size_t)n, + num_age_groups * num_regions + (size_t)i * num_regions + (size_t)n) = 1.0 / T_Ii; + } + } + + V = V.inverse(); + + Eigen::MatrixXd NextGenMatrix = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + NextGenMatrix = F * V; + + //Compute the largest eigenvalue in absolute value + Eigen::ComplexEigenSolver ces; + + ces.compute(NextGenMatrix); + const Eigen::VectorXcd eigen_vals = ces.eigenvalues(); + + Eigen::VectorXd eigen_vals_abs; + eigen_vals_abs.resize(eigen_vals.size()); + + for (int i = 0; i < eigen_vals.size(); i++) { + eigen_vals_abs[i] = std::abs(eigen_vals[i]); + } + return mio::success(eigen_vals_abs.maxCoeff()); + } + + /** + *@brief Computes the reproduction number for all time points of the Model output obtained by the Simulation. + *@param y The TimeSeries obtained from the Model Simulation. + *@returns vector containing all reproduction numbers + */ + Eigen::VectorXd get_reproduction_numbers(const mio::TimeSeries& y) + { + auto num_time_points = y.get_num_time_points(); + Eigen::VectorXd temp(num_time_points); + for (size_t i = 0; i < static_cast(num_time_points); i++) { + temp[i] = get_reproduction_number(i, y).value(); + } + return temp; + } +}; + +} // namespace oseirmobility +} // namespace mio + +#endif // ODESEIRMOBILITY_MODEL_H diff --git a/cpp/models/ode_seir_mobility/parameters.h b/cpp/models/ode_seir_mobility/parameters.h new file mode 100644 index 0000000000..ca6fd99f1f --- /dev/null +++ b/cpp/models/ode_seir_mobility/parameters.h @@ -0,0 +1,244 @@ + +#ifndef SEIRMOBILITY_PARAMETERS_H +#define SEIRMOBILITY_PARAMETERS_H + +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_seir_mobility/regions.h" + +#include + +namespace mio +{ +namespace oseirmobility +{ + +/**************************************************** + * Define Parameters of the SEIR model with mobility * + ****************************************************/ + +/** + * @brief Probability of getting infected from a contact. + */ +template +struct TransmissionProbabilityOnContact { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 1.0); + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/** + * @brief the latent time in day unit + */ +template +struct TimeExposed { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 5.2); + } + static std::string name() + { + return "TimeExposed"; + } +}; + +/** + * @brief The infectious time in day unit. + */ +template +struct TimeInfected { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 6.0); + } + static std::string name() + { + return "TimeInfected"; + } +}; + +/** + * @brief The contact patterns within the society are modelled using a ContactMatrix. + */ +template +struct ContactPatterns { + using Type = UncertainContactMatrix; + static Type get_default(Region, AgeGroup size) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/** + * @brief The contact patterns between different Region%s are modelled using a ContactMatrix. + */ +template +struct CommutingStrengths { + using Type = UncertainContactMatrix; + static Type get_default(Region size, AgeGroup) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "CommutingStrengths"; + } +}; + +template +using ParametersBase = ParameterSet, TimeExposed, TimeInfected, + ContactPatterns, CommutingStrengths>; + +/** + * @brief Parameters of SEIR model. + */ +template +class Parameters : public ParametersBase +{ +public: + Parameters(Region num_regions, AgeGroup num_agegroups) + : ParametersBase(num_regions, num_agegroups) + , m_num_regions{num_regions} + , m_num_agegroups(num_agegroups) + { + } + + Region get_num_regions() const + { + return m_num_regions; + } + + AgeGroup get_num_agegroups() const + { + return m_num_agegroups; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and applies them, if they do not. + * Time spans cannot be negative and probabilities can only take values between [0,1]. + * + * Attention: This function should be used with care. It is necessary for some test problems to run through quickly, + * but in a manual execution of an example, check_constraints() may be preferred. Note that the apply_constraints() + * function can and will not set Parameters to meaningful values in an epidemiological or virological context, + * as all models are designed to be transferable to multiple diseases. Consequently, only acceptable + * (like 0 or 1 for probabilities or small positive values for time spans) values are set here and a manual adaptation + * may often be necessary to have set meaningful values. + * + * @return Returns true if one ore more constraint were corrected, false otherwise. + */ + bool apply_constraints() + { + double tol_times = 1e-1; + + int corrected = false; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + if (this->template get>()[i] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], tol_times); + this->template get>()[i] = tol_times; + corrected = true; + } + if (this->template get>()[i] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], tol_times); + this->template get>()[i] = tol_times; + corrected = true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_warning( + "Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", + this->template get>()[i], 0.0); + this->template get>() = 0.0; + corrected = true; + } + } + return corrected; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and logs an error + * if constraints are not satisfied. + * @return Returns true if one constraint is not satisfied, otherwise false. + */ + bool check_constraints() const + { + double tol_times = 1e-1; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + if (this->template get>()[i] < tol_times) { + log_error( + "Constraint check: Parameter TimeExposed {:.4f} smaller or equal {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], 0.0); + return true; + } + if (this->template get>()[i] < tol_times) { + log_error( + "Constraint check: Parameter TimeInfected {:.4f} smaller or equal {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], 0.0); + return true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact {:.4f} smaller {:.4f} or " + "greater {:.4f}", + this->template get>()[i], 0.0, 1.0); + return true; + } + } + return false; + } + +private: + // Parameters(ParametersBase&& base) + // : ParametersBase(std::move(base)) //TODO: Adjust + // { + // } + +public: + /** + * deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); + return success(Parameters(std::move(base))); + } + +private: + Region m_num_regions; + AgeGroup m_num_agegroups; +}; + +} // namespace oseirmobility +} // namespace mio + +#endif // SEIR_PARAMETERS_H diff --git a/cpp/models/ode_seir_mobility/regions.h b/cpp/models/ode_seir_mobility/regions.h new file mode 100644 index 0000000000..4a8352b11b --- /dev/null +++ b/cpp/models/ode_seir_mobility/regions.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITY_REGIONS_H +#define ODESEIRMOBILITY_REGIONS_H + +#include "memilio/utils/index.h" + +namespace mio +{ +namespace oseirmobility +{ + +/** + * @brief The AgeGroup struct is used as a dynamically + * sized tag for all age dependent categories + */ +struct Region : public Index { + Region(size_t val) + : Index(val) + { + } +}; + +} // namespace oseirmobility +} // namespace mio + +#endif diff --git a/cpp/models/ode_seir_mobility_improved/CMakeLists.txt b/cpp/models/ode_seir_mobility_improved/CMakeLists.txt new file mode 100644 index 0000000000..82261701db --- /dev/null +++ b/cpp/models/ode_seir_mobility_improved/CMakeLists.txt @@ -0,0 +1,13 @@ +add_library(ode_seir_mobility_improved + infection_state.h + model.h + model.cpp + parameters.h + regions.h +) +target_link_libraries(ode_seir_mobility_improved PUBLIC memilio) +target_include_directories(ode_seir_mobility_improved PUBLIC + $ + $ +) +target_compile_options(ode_seir_mobility_improved PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_seir_mobility_improved/infection_state.h b/cpp/models/ode_seir_mobility_improved/infection_state.h new file mode 100644 index 0000000000..750a6415b7 --- /dev/null +++ b/cpp/models/ode_seir_mobility_improved/infection_state.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITYIMPROVED_INFECTIONSTATE_H +#define ODESEIRMOBILITYIMPROVED_INFECTIONSTATE_H + +namespace mio +{ +namespace oseirmobilityimproved +{ + +/** + * @brief The InfectionState enum describes the possible + * categories for the infectious state of persons + */ +enum class InfectionState +{ + Susceptible, + Exposed, + Infected, + Recovered, + Count +}; + +} // namespace oseirmobilityimproved +} // namespace mio + +#endif // ODESEIR_INFECTIONSTATE_H diff --git a/cpp/models/ode_seir_mobility_improved/model.cpp b/cpp/models/ode_seir_mobility_improved/model.cpp new file mode 100644 index 0000000000..567a8a1e86 --- /dev/null +++ b/cpp/models/ode_seir_mobility_improved/model.cpp @@ -0,0 +1,10 @@ + +#include "ode_seir_mobility_improved/model.h" + +namespace mio +{ +namespace oseirmobilityimproved +{ + +} // namespace oseirmobilityimproved +} // namespace mio diff --git a/cpp/models/ode_seir_mobility_improved/model.h b/cpp/models/ode_seir_mobility_improved/model.h new file mode 100644 index 0000000000..b8469169b2 --- /dev/null +++ b/cpp/models/ode_seir_mobility_improved/model.h @@ -0,0 +1,212 @@ + +#ifndef ODESEIRMOBILITYIMPROVED_MODEL_H +#define ODESEIRMOBILITYIMPROVED_MODEL_H + +#include "memilio/compartments/flow_model.h" +#include "memilio/epidemiology/populations.h" +#include "models/ode_seir_mobility_improved/infection_state.h" +#include "models/ode_seir_mobility_improved/parameters.h" +#include "models/ode_seir_mobility_improved/regions.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/time_series.h" + +GCC_CLANG_DIAGNOSTIC(push) +GCC_CLANG_DIAGNOSTIC(ignored "-Wshadow") +#include +GCC_CLANG_DIAGNOSTIC(pop) + +namespace mio +{ +namespace oseirmobilityimproved +{ + +/******************** + * define the model * + ********************/ + +using Flows = TypeList, + Flow, + Flow>; + +template +class Model : public FlowModel, + Parameters, Flows> +{ + + using Base = + FlowModel, Parameters, Flows>; + +public: + using typename Base::ParameterSet; + using typename Base::Populations; + + Model(int num_regions, int num_agegroups) + : Base(Populations({Region(num_regions), AgeGroup(num_agegroups), InfectionState::Count}), + ParameterSet(Region(num_regions), AgeGroup(num_agegroups))) + , m_population_after_commuting( + mio::Populations({Region(num_regions), AgeGroup(num_agegroups)})) + { + } + + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override + { + const auto& params = this->parameters; + const auto& population = this->populations; + const auto& commuting_strengths = + params.template get>().get_cont_freq_mat().get_matrix_at(t); + const Index n_age_groups = reduce_index>(params.get_num_agegroups()); + const Index n_regions = reduce_index>(params.get_num_regions()); + for (size_t age_i = 0; age_i < (size_t)n_age_groups; age_i++) { + for (size_t age_j = 0; age_j < (size_t)n_age_groups; age_j++) { + Eigen::VectorXd infectious_share_per_region = Eigen::VectorXd::Zero((size_t)n_regions); + for (size_t region_n = 0; region_n < (size_t)n_regions; region_n++) { + for (size_t region_m = 0; region_m < (size_t)n_regions; region_m++) { + infectious_share_per_region(region_n) += + commuting_strengths(region_m, region_n) * + pop[population.get_flat_index( + {Region(region_m), AgeGroup(age_j), InfectionState::Infected})]; + } + infectious_share_per_region(region_n) /= + m_population_after_commuting[{Region(region_n), AgeGroup(age_j)}]; + } + Eigen::VectorXd infections_due_commuting = commuting_strengths * infectious_share_per_region; + for (size_t region_n = 0; region_n < (size_t)n_regions; region_n++) { + const size_t Ejn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), InfectionState::Exposed}); + const size_t Ijn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), InfectionState::Infected}); + const size_t Rjn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), InfectionState::Recovered}); + const size_t Sjn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), InfectionState::Susceptible}); + + const double Nj_inv = 1.0 / (pop[Sjn] + pop[Ejn] + pop[Ijn] + pop[Rjn]); + double coeffStoI = + 0.5 * + params.template get>().get_cont_freq_mat().get_matrix_at(t)(age_i, age_j) * + params.template get>()[AgeGroup(age_i)]; + + flows[Base::template get_flat_flow_index( + {Region(region_n), AgeGroup(age_i)})] += + (pop[Ijn] * Nj_inv + infections_due_commuting(region_n)) * coeffStoI * + y[population.get_flat_index({Region(region_n), AgeGroup(age_i), InfectionState::Susceptible})]; + } + } + for (size_t region = 0; region < (size_t)n_regions; region++) { + flows[Base::template get_flat_flow_index( + {Region(region), AgeGroup(age_i)})] = + y[population.get_flat_index({Region(region), AgeGroup(age_i), InfectionState::Exposed})] / + params.template get>()[AgeGroup(age_i)]; + flows[Base::template get_flat_flow_index( + {Region(region), AgeGroup(age_i)})] = + y[population.get_flat_index({Region(region), AgeGroup(age_i), InfectionState::Infected})] / + params.template get>()[AgeGroup(age_i)]; + } + } + } + + /** + *@brief Computes the reproduction number at a given index time of the Model output obtained by the Simulation. + *@param t_idx The index time at which the reproduction number is computed. + *@param y The TimeSeries obtained from the Model Simulation. + *@returns The computed reproduction number at the provided index time. + */ + IOResult get_reproduction_number(size_t t_idx, const mio::TimeSeries& y) + { + if (!(t_idx < static_cast(y.get_num_time_points()))) { + return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries"); + } + + auto const& params = this->parameters; + auto const& pop = this->populations; + + const size_t num_age_groups = (size_t)params.get_num_agegroups(); + const size_t num_regions = (size_t)params.get_num_regions(); + constexpr size_t num_infected_compartments = 2; + const size_t total_infected_compartments = num_infected_compartments * num_age_groups * num_regions; + + ContactMatrixGroup const& contact_matrix = params.template get>(); + ContactMatrixGroup const& commuting_strengths = params.template get>(); + + Eigen::MatrixXd F = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + Eigen::MatrixXd V = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + + for (auto i = AgeGroup(0); i < AgeGroup(num_age_groups); i++) { + for (auto n = Region(0); n < Region(num_regions); n++) { + size_t Si = pop.get_flat_index({n, i, InfectionState::Susceptible}); + for (auto j = AgeGroup(0); j < AgeGroup(num_age_groups); j++) { + for (auto m = Region(0); m < Region(num_regions); m++) { + auto const population_region = pop.template slice({m.get(), 1}); + auto const population_region_age = population_region.template slice({j.get(), 1}); + auto Njm = std::accumulate(population_region_age.begin(), population_region_age.end(), 0.); + + double coeffStoE = 0.5 * contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) * + params.template get>()[i]; + if (n == m) { + F(i.get() * num_regions + n.get(), num_age_groups * num_regions + j.get() * num_regions + + m.get()) += coeffStoE * y.get_value(t_idx)[Si] / Njm; + } + for (auto k = Region(0); k < Region(num_regions); k++) { + F(i.get() * num_regions + n.get(), + num_age_groups * num_regions + j.get() * num_regions + m.get()) += + coeffStoE * y.get_value(t_idx)[Si] * + commuting_strengths.get_matrix_at(y.get_time(t_idx))(n.get(), k.get()) * + commuting_strengths.get_matrix_at(y.get_time(t_idx))(m.get(), k.get()) / + m_population_after_commuting[{k, j}]; + } + } + } + + double T_Ei = params.template get>()[i]; + double T_Ii = params.template get>()[i]; + V(i.get() * num_regions + n.get(), i.get() * num_regions + n.get()) = 1.0 / T_Ei; + V(num_age_groups * num_regions + i.get() * num_regions + n.get(), i.get() * num_regions + n.get()) = + -1.0 / T_Ei; + V(num_age_groups * num_regions + i.get() * num_regions + n.get(), + num_age_groups * num_regions + i.get() * num_regions + n.get()) = 1.0 / T_Ii; + } + } + + V = V.inverse(); + + Eigen::MatrixXd NextGenMatrix = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + NextGenMatrix = F * V; + + //Compute the largest eigenvalue in absolute value + Eigen::ComplexEigenSolver ces; + + ces.compute(NextGenMatrix); + const Eigen::VectorXcd eigen_vals = ces.eigenvalues(); + + Eigen::VectorXd eigen_vals_abs; + eigen_vals_abs.resize(eigen_vals.size()); + + for (int i = 0; i < eigen_vals.size(); i++) { + eigen_vals_abs[i] = std::abs(eigen_vals[i]); + } + return mio::success(eigen_vals_abs.maxCoeff()); + } + + /** + *@brief Computes the reproduction number for all time points of the Model output obtained by the Simulation. + *@param y The TimeSeries obtained from the Model Simulation. + *@returns vector containing all reproduction numbers + */ + Eigen::VectorXd get_reproduction_numbers(const mio::TimeSeries& y) + { + auto num_time_points = y.get_num_time_points(); + Eigen::VectorXd temp(num_time_points); + for (size_t i = 0; i < static_cast(num_time_points); i++) { + temp[i] = get_reproduction_number(i, y).value(); + } + return temp; + } + + mio::Populations m_population_after_commuting; +}; // namespace oseirmobilityimproved + +} // namespace oseirmobilityimproved +} // namespace mio + +#endif // ODESEIRMOBILITY_MODEL_H diff --git a/cpp/models/ode_seir_mobility_improved/parameters.h b/cpp/models/ode_seir_mobility_improved/parameters.h new file mode 100644 index 0000000000..49c548207f --- /dev/null +++ b/cpp/models/ode_seir_mobility_improved/parameters.h @@ -0,0 +1,261 @@ + +#ifndef SEIRMOBILITY_PARAMETERS_H +#define SEIRMOBILITY_PARAMETERS_H + +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_seir_mobility_improved/regions.h" +#include "Eigen/Sparse" + +#include + +namespace mio +{ +namespace oseirmobilityimproved +{ + +/**************************************************** + * Define Parameters of the SEIR model with mobility * + ****************************************************/ + +/** + * @brief Probability of getting infected from a contact. + */ +template +struct TransmissionProbabilityOnContact { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 1.0); + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/** + * @brief the latent time in day unit + */ +template +struct TimeExposed { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 5.2); + } + static std::string name() + { + return "TimeExposed"; + } +}; + +/** + * @brief The infectious time in day unit. + */ +template +struct TimeInfected { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 6.0); + } + static std::string name() + { + return "TimeInfected"; + } +}; + +/** + * @brief The contact patterns within the society are modelled using a ContactMatrix. + */ +template +struct ContactPatterns { + using Type = UncertainContactMatrix; + static Type get_default(Region, AgeGroup size) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/** + * @brief The contact patterns between different Region%s are modelled using a ContactMatrix. + */ +template +struct CommutingStrengths { + using Type = UncertainContactMatrix; + static Type get_default(Region size, AgeGroup) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "CommutingStrengths"; + } +}; + +/** + * @brief The sizes of the populations after commuting. + */ +template +struct PopulationSizes { + using Type = CustomIndexArray; + static Type get_default(Region size, AgeGroup) + { + return Type(size, 0.); + } + static std::string name() + { + return "PopulationSizes"; + } +}; + +template +using ParametersBase = ParameterSet, TimeExposed, TimeInfected, + ContactPatterns, CommutingStrengths, PopulationSizes>; + +/** + * @brief Parameters of SEIR model. + */ +template +class Parameters : public ParametersBase +{ +public: + Parameters(Region num_regions, AgeGroup num_agegroups) + : ParametersBase(num_regions, num_agegroups) + , m_num_regions{num_regions} + , m_num_agegroups(num_agegroups) + { + } + + Region get_num_regions() const + { + return m_num_regions; + } + + AgeGroup get_num_agegroups() const + { + return m_num_agegroups; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and applies them, if they do not. + * Time spans cannot be negative and probabilities can only take values between [0,1]. + * + * Attention: This function should be used with care. It is necessary for some test problems to run through quickly, + * but in a manual execution of an example, check_constraints() may be preferred. Note that the apply_constraints() + * function can and will not set Parameters to meaningful values in an epidemiological or virological context, + * as all models are designed to be transferable to multiple diseases. Consequently, only acceptable + * (like 0 or 1 for probabilities or small positive values for time spans) values are set here and a manual adaptation + * may often be necessary to have set meaningful values. + * + * @return Returns true if one ore more constraint were corrected, false otherwise. + */ + bool apply_constraints() + { + double tol_times = 1e-1; + + int corrected = false; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + if (this->template get>()[i] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], tol_times); + this->template get>()[i] = tol_times; + corrected = true; + } + if (this->template get>()[i] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], tol_times); + this->template get>()[i] = tol_times; + corrected = true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_warning( + "Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", + this->template get>()[i], 0.0); + this->template get>() = 0.0; + corrected = true; + } + } + return corrected; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and logs an error + * if constraints are not satisfied. + * @return Returns true if one constraint is not satisfied, otherwise false. + */ + bool check_constraints() const + { + double tol_times = 1e-1; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + if (this->template get>()[i] < tol_times) { + log_error( + "Constraint check: Parameter TimeExposed {:.4f} smaller or equal {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], 0.0); + return true; + } + if (this->template get>()[i] < tol_times) { + log_error( + "Constraint check: Parameter TimeInfected {:.4f} smaller or equal {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], 0.0); + return true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact {:.4f} smaller {:.4f} or " + "greater {:.4f}", + this->template get>()[i], 0.0, 1.0); + return true; + } + } + return false; + } + +private: + // Parameters(ParametersBase&& base) + // : ParametersBase(std::move(base)) //TODO: Adjust + // { + // } + +public: + /** + * deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); + return success(Parameters(std::move(base))); + } + +private: + Region m_num_regions; + AgeGroup m_num_agegroups; +}; + +} // namespace oseirmobilityimproved +} // namespace mio + +#endif // SEIR_PARAMETERS_H diff --git a/cpp/models/ode_seir_mobility_improved/regions.h b/cpp/models/ode_seir_mobility_improved/regions.h new file mode 100644 index 0000000000..d5f7657138 --- /dev/null +++ b/cpp/models/ode_seir_mobility_improved/regions.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITYIMPROVED_REGIONS_H +#define ODESEIRMOBILITYIMPROVED_REGIONS_H + +#include "memilio/utils/index.h" + +namespace mio +{ +namespace oseirmobilityimproved +{ + +/** + * @brief The Region struct is used as a dynamically + * sized tag for all region dependent categories + */ +struct Region : public Index { + Region(size_t val) + : Index(val) + { + } +}; + +} // namespace oseirmobilityimproved +} // namespace mio + +#endif diff --git a/cpp/models/ode_sir_mobility/CMakeLists.txt b/cpp/models/ode_sir_mobility/CMakeLists.txt new file mode 100644 index 0000000000..3a2f54adeb --- /dev/null +++ b/cpp/models/ode_sir_mobility/CMakeLists.txt @@ -0,0 +1,13 @@ +add_library(ode_sir_mobility + infection_state.h + model.h + model.cpp + parameters.h + regions.h +) +target_link_libraries(ode_sir_mobility PUBLIC memilio) +target_include_directories(ode_sir_mobility PUBLIC + $ + $ +) +target_compile_options(ode_sir_mobility PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_sir_mobility/README.md b/cpp/models/ode_sir_mobility/README.md new file mode 100644 index 0000000000..77f5698546 --- /dev/null +++ b/cpp/models/ode_sir_mobility/README.md @@ -0,0 +1,21 @@ + +# ODE SIR compartment model + +This model is a very simple ODE model with only three compartments and few parameters, mostly for demonstration of the MEmilio framework: +- Susceptible, may become infected at any time +- Infected, will be recovered after some time +- Recovered, recovered from infectious process (dead or recovered) + +We assume simulations over short periods of time, so that the population size can be considered constant and birth as well as (natural) mortality rates can be ignored. + +Below is an overview of the model architecture and its compartments. + +![SIR_model](https://github.com/SciCompMod/memilio/assets/69154294/01c9a2ae-2f5c-4bad-b7f0-34de651f2c73) +| Mathematical variable | C++ variable name | Description | +|---------------------------- | --------------- | -------------------------------------------------------------------------------------------------- | +| $\phi$ | `ContactPatterns` | Daily contact rate / Number of daily contacts. | +| $\rho$ | `TransmissionProbabilityOnContact` | Transmission risk for people located in the Susceptible compartment. | +| $N$ | `populations.get_total()` | Total population. | +| $T_{I}$ | `TimeInfected` | Time in days an individual stays in the Infected compartment. | + +An example can be found in [examples/ode_sir.cpp](../../examples/ode_sir.cpp) diff --git a/cpp/models/ode_sir_mobility/infection_state.h b/cpp/models/ode_sir_mobility/infection_state.h new file mode 100644 index 0000000000..dc98471c73 --- /dev/null +++ b/cpp/models/ode_sir_mobility/infection_state.h @@ -0,0 +1,25 @@ + +#ifndef ODESIRMOBILITY_INFECTIONSTATE_H +#define ODESIRMOBILITY_INFECTIONSTATE_H + +namespace mio +{ +namespace osirmobility +{ + +/** + * @brief The InfectionState enum describes the possible + * categories for the infectious state of persons + */ +enum class InfectionState +{ + Susceptible, + Infected, + Recovered, + Count +}; + +} // namespace osir +} // namespace mio + +#endif // ODESIR_INFECTIONSTATE_H diff --git a/cpp/models/ode_sir_mobility/model.cpp b/cpp/models/ode_sir_mobility/model.cpp new file mode 100644 index 0000000000..0f8bf7b573 --- /dev/null +++ b/cpp/models/ode_sir_mobility/model.cpp @@ -0,0 +1,10 @@ + +#include "ode_sir_mobility/model.h" + +namespace mio +{ +namespace osirmobility +{ + +} // namespace osir +} // namespace mio diff --git a/cpp/models/ode_sir_mobility/model.h b/cpp/models/ode_sir_mobility/model.h new file mode 100644 index 0000000000..30ef675728 --- /dev/null +++ b/cpp/models/ode_sir_mobility/model.h @@ -0,0 +1,114 @@ + +#ifndef ODESIRMOBILITY_MODEL_H +#define ODESIRMOBILITY_MODEL_H + +#include "memilio/compartments/flow_model.h" +#include "memilio/epidemiology/populations.h" +#include "ode_sir_mobility/infection_state.h" +#include "ode_sir_mobility/parameters.h" +#include "ode_sir_mobility/regions.h" +#include "memilio/epidemiology/age_group.h" + +namespace mio +{ +namespace osirmobility +{ + +/******************** + * define the model * + ********************/ + +using Flows = TypeList, + Flow>; + +template +class Model : public FlowModel, + Parameters, Flows> +{ + + using Base = + FlowModel, Parameters, Flows>; + +public: + using typename Base::ParameterSet; + using typename Base::Populations; + + Model(int num_regions, int num_agegroups) + : Base(Populations({Region(num_regions), AgeGroup(num_agegroups), InfectionState::Count}), + ParameterSet(Region(num_regions), AgeGroup(num_agegroups))) + { + } + // Einmal über den Vektor und später nochmal über die Regions + + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override + { + const auto& params = this->parameters; + const auto& population = this->populations; + + const Index n_age_groups = reduce_index>(params.get_num_agegroups()); + const Index n_regions = reduce_index>(params.get_num_regions()); + + for (auto age_i : make_index_range(n_age_groups)) { + for (auto age_j : make_index_range(n_age_groups)) { + double coeffStoI = params.template get>().get_cont_freq_mat().get_matrix_at(t)( + age_i.get(), age_j.get()) * + params.template get>()[age_i] / + population.get_group_total(age_j); + for (auto edge : params.template get()) { + auto start_region = get<0>(edge); + auto end_region = get<1>(edge); + auto strength = get(edge); + if (start_region == end_region) { + continue; + } + // s_n += h_mn/P_m * i_m + flows[Base::template get_flat_flow_index( + {start_region, age_i})] += + strength * pop[population.get_flat_index({end_region, age_j, InfectionState::Infected})]; + // s_m += h_mn/P_m * i_n + flows[Base::template get_flat_flow_index( + {end_region, age_i})] += + strength * pop[population.get_flat_index({start_region, age_j, InfectionState::Infected})]; + + // s_n += gamma * h_nm/P_n * sum(h_km/P_k * p_nm,k * i_k) + for (auto edge_commuter : params.template get()) { + auto start_region_commuter = get<0>(edge_commuter); + auto end_region_commuter = get<1>(edge_commuter); + auto strength_commuter = get(edge_commuter); + if (end_region_commuter != end_region || start_region_commuter == start_region || + ((std::find(params.template get()[{start_region, end_region}].begin(), + params.template get()[{start_region, end_region}].end(), + start_region_commuter)) == + params.template get()[{start_region, end_region}].end())) { + continue; + } + flows[Base::template get_flat_flow_index( + {start_region, age_i})] += + params.template get>() * strength * + strength_commuter * + pop[population.get_flat_index({start_region_commuter, age_j, InfectionState::Infected})]; + } + } + for (auto region : make_index_range(n_regions)) { + flows[Base::template get_flat_flow_index( + {region, age_i})] += pop[population.get_flat_index({region, age_j, InfectionState::Infected})]; + flows[Base::template get_flat_flow_index( + {region, age_i})] *= + coeffStoI * y[population.get_flat_index({region, age_j, InfectionState::Susceptible})]; + } + } + + for (auto region : make_index_range(n_regions)) { + flows[Base::template get_flat_flow_index( + {region, age_i})] = (1.0 / params.template get>()[age_i]) * + y[population.get_flat_index({region, age_i, InfectionState::Infected})]; + } + } + } +}; + +} // namespace osirmobility +} // namespace mio + +#endif // ODESIRMOBILITY_MODEL_H diff --git a/cpp/models/ode_sir_mobility/parameters.h b/cpp/models/ode_sir_mobility/parameters.h new file mode 100644 index 0000000000..18f7110543 --- /dev/null +++ b/cpp/models/ode_sir_mobility/parameters.h @@ -0,0 +1,286 @@ + +#ifndef SIRMOBILITY_PARAMETERS_H +#define SIRMOBILITY_PARAMETERS_H + +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/utils/custom_index_array.h" +#include "ode_sir_mobility/regions.h" + +#include + +namespace mio +{ +namespace osirmobility +{ + +/**************************************************** + * Define Parameters of the SIR model with mobility * + ****************************************************/ + +/** + * @brief Probability of getting infected from a contact. + */ +template +struct TransmissionProbabilityOnContact { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 1.0); + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/** + * @brief The infectious time in day unit. + */ +template +struct TimeInfected { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 6.0); + } + static std::string name() + { + return "TimeInfected"; + } +}; + +/** + * @brief The contact patterns within the society are modelled using a ContactMatrix. + */ +template +struct ContactPatterns { + using Type = UncertainContactMatrix; + static Type get_default(Region, AgeGroup size) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/** + * @brief The mean number of people migrating from one Region to another during a TimeStep. + */ +struct CommutingRatio { + using Type = std::vector>; + static Type get_default(Region, AgeGroup) + { + return Type({{Region(0), Region(0), 0.}}); + } + static std::string name() + { + return "CommutingRatio"; + } +}; + +/** + * @brief The ratio that regulates the infections during commuting. +*/ +template +struct ImpactTransmissionDuringCommuting { + using Type = UncertainValue; + static Type get_default(Region, AgeGroup) + { + return Type(0.); + } + static std::string name() + { + return "ImpactTransmissionDuringCommuting"; + } +}; + +/** + * @brief The Region%s that a person crosses when travelling from one Region to another. +*/ +struct PathIntersections { + using Type = CustomIndexArray, Region, Region>; + static Type get_default(Region size, AgeGroup) + { + return Type({size, size}); + } + static std::string name() + { + return "PathIntersections"; + } +}; + +template +using ParametersBase = ParameterSet, TimeInfected, ContactPatterns, + CommutingRatio, ImpactTransmissionDuringCommuting, PathIntersections>; + +/** + * @brief Parameters of SIR model. + */ +template +class Parameters : public ParametersBase +{ +public: + Parameters(Region num_regions, AgeGroup num_agegroups) + : ParametersBase(num_regions, num_agegroups) + , m_num_regions{num_regions} + , m_num_agegroups(num_agegroups) + { + } + + Region get_num_regions() const + { + return m_num_regions; + } + + AgeGroup get_num_agegroups() const + { + return m_num_agegroups; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and applies them, if they do not. + * Time spans cannot be negative and probabilities can only take values between [0,1]. + * + * Attention: This function should be used with care. It is necessary for some test problems to run through quickly, + * but in a manual execution of an example, check_constraints() may be preferred. Note that the apply_constraints() + * function can and will not set Parameters to meaningful values in an epidemiological or virological context, + * as all models are designed to be transferable to multiple diseases. Consequently, only acceptable + * (like 0 or 1 for probabilities or small positive values for time spans) values are set here and a manual adaptation + * may often be necessary to have set meaningful values. + * + * @return Returns true if one ore more constraint were corrected, false otherwise. + */ + bool apply_constraints() + { + double tol_times = 1e-1; + + int corrected = false; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + if (this->template get>()[i] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], tol_times); + this->template get>()[i] = tol_times; + corrected = true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_warning( + "Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", + this->template get>()[i], 0.0); + this->template get>() = 0.0; + corrected = true; + } + } + if (this->template get>() < 0.0 || + this->template get>() > 1.0) { + log_warning("Constraint check: Parameter ImpactTransmissionDuringCommuting changed from {:.4f} to {:.4f}.", + this->template get>(), 0.0); + this->template get>() = 0.0; + corrected = true; + } + for (auto& i : this->template get()) { + if (std::get(i) < 0.0 || std::get(i) > 1.0) { + log_warning("Constraint check: Parameter CommutingRatio changed from {:.4f} to {:.4f}.", + std::get(i), 0.0); + std::get(i) = 0.0; + corrected = true; + } + if (std::get<0>(i) < Region(0) || std::get<1>(i) < Region(0) || std::get<0>(i) >= m_num_regions || + std::get<1>(i) >= m_num_regions) { + log_warning( + "Constraint check: Removed entry of Parameter CommutingRatio because of non-existing Regions."); + auto it = std::find(this->template get().begin(), + this->template get().end(), i); + this->template get().erase(it); + corrected = true; + } + } + return corrected; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and logs an error + * if constraints are not satisfied. + * @return Returns true if one constraint is not satisfied, otherwise false. + */ + bool check_constraints() const + { + double tol_times = 1e-1; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + + if (this->template get>()[i] < tol_times) { + log_error( + "Constraint check: Parameter TimeInfected {:.4f} smaller or equal {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], 0.0); + return true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact {:.4f} smaller {:.4f} or " + "greater {:.4f}", + this->template get>()[i], 0.0, 1.0); + return true; + } + } + if (this->template get>() < 0.0 || + this->template get>() > 1.0) { + log_error( + "Constraint check: Parameter ImpactTransmissionDuringCommuting {:.4f} smaller {:.4f} or greater {:.4f}", + this->template get>(), 0.0, 1.0); + return true; + } + for (auto i : this->template get()) { + if (std::get(i) < 0.0 || std::get(i) > 1.0) { + log_error("Constraint check: Parameter CommutingRatio entry {:.4f} smaller {:.4f} or greater {:.4f}", + std::get(i), 0.0, 1.0); + return true; + } + if (std::get<0>(i) < Region(0) || std::get<1>(i) < Region(0) || std::get<0>(i) > m_num_regions || + std::get<1>(i) > m_num_regions) { + log_error("Constraint check: Parameter CommutingRatio has an entry with start or end Region " + "that does not appear in the model."); + return true; + } + } + return false; + } + +private: + // Parameters(ParametersBase&& base) + // : ParametersBase(std::move(base)) //TODO: Adjust + // { + // } + +public: + /** + * deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); + return success(Parameters(std::move(base))); + } + +private: + Region m_num_regions; + AgeGroup m_num_agegroups; +}; + +} // namespace osirmobility +} // namespace mio + +#endif // SIR_PARAMETERS_H diff --git a/cpp/models/ode_sir_mobility/regions.h b/cpp/models/ode_sir_mobility/regions.h new file mode 100644 index 0000000000..d3c8b528a7 --- /dev/null +++ b/cpp/models/ode_sir_mobility/regions.h @@ -0,0 +1,26 @@ + +#ifndef ODESIRMOBILITY_REGIONS_H +#define ODESIRMOBILITY_REGIONS_H + +#include "memilio/utils/index.h" + +namespace mio +{ +namespace osirmobility +{ + +/** + * @brief The AgeGroup struct is used as a dynamically + * sized tag for all age dependent categories + */ +struct Region : public Index { + Region(size_t val) + : Index(val) + { + } +}; + +} // namespace osirmobility +} // namespace mio + +#endif diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index fda62c5ef7..5375ff3c68 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -9,6 +9,7 @@ set(TESTSOURCES test_populations.cpp test_odeseir.cpp test_odesir.cpp + test_odesirmobility.cpp test_numericalIntegration.cpp test_smoother.cpp test_damping.cpp diff --git a/cpp/tests/test_odesirmobility.cpp b/cpp/tests/test_odesirmobility.cpp new file mode 100644 index 0000000000..247a8622ff --- /dev/null +++ b/cpp/tests/test_odesirmobility.cpp @@ -0,0 +1,248 @@ + +#include "load_test_data.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "ode_sir_mobility/model.h" +#include "ode_sir_mobility/infection_state.h" +#include "ode_sir_mobility/parameters.h" +#include "memilio/math/euler.h" +#include "memilio/compartments/simulation.h" +#include +#include +#include + +TEST(TestOdeSirMobility, simulateDefault) +{ + double t0 = 0; + double tmax = 1; + double dt = 0.1; + + size_t num_regions = 4; + + mio::osirmobility::Model model(num_regions); + mio::TimeSeries result = simulate(t0, tmax, dt, model); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); +} + +TEST(TestOdeSirMobility, compareWithPreviousRun) +{ + // initialization + double t0 = 0.; + double tmax = 3.; + double dt = 0.1; + + size_t number_regions = 4; + size_t number_age_groups = 1; + size_t total_population_per_region = 5000; + + mio::osirmobility::Model model(number_regions, number_age_groups); + + for (size_t i = 0; i < number_regions; i++) { + model.populations[{mio::Index( + mio::osirmobility::Region(i), mio::AgeGroup(0), mio::osirmobility::InfectionState::Infected)}] = 50; + model.populations[{mio::Index( + mio::osirmobility::Region(i), mio::AgeGroup(0), mio::osirmobility::InfectionState::Recovered)}] = 0; + model.populations[{mio::Index( + mio::osirmobility::Region(i), mio::AgeGroup(0), mio::osirmobility::InfectionState::Susceptible)}] = + total_population_per_region - + model.populations[{mio::Index( + mio::osirmobility::Region(i), mio::AgeGroup(0), mio::osirmobility::InfectionState::Infected)}] - + model.populations[{mio::Index( + mio::osirmobility::Region(i), mio::AgeGroup(0), mio::osirmobility::InfectionState::Recovered)}]; + } + model.parameters.set(1.0); + model.parameters.set(2); + + model.parameters.get().get_baseline()(0, 0) = 2.7; + model.parameters.get().add_damping(0.6, mio::SimulationTime(12.5)); + + model.parameters.set(1.); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(0), 0.2}); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(2), 0.6}); + model.parameters.get().push_back( + {mio::osirmobility::Region(2), mio::osirmobility::Region(0), 0.5}); + model.parameters.get().push_back( + {mio::osirmobility::Region(0), mio::osirmobility::Region(3), 1.0}); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(3), 0.2}); + + model.parameters.get()[{mio::osirmobility::Region(0), + mio::osirmobility::Region(1)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(0), + mio::osirmobility::Region(3)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(1), + mio::osirmobility::Region(0)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(1), + mio::osirmobility::Region(2)}] = {0}; + model.parameters.get()[{mio::osirmobility::Region(1), + mio::osirmobility::Region(3)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(2), + mio::osirmobility::Region(1)}] = {0}; + model.parameters.get()[{mio::osirmobility::Region(3), + mio::osirmobility::Region(0)}] = {2}; + model.parameters.get()[{mio::osirmobility::Region(3), + mio::osirmobility::Region(1)}] = {2}; + + std::vector> refData = load_test_data_csv("ode-sir-mobility-compare.csv"); + auto integrator = std::make_shared(); + auto result = mio::simulate(t0, tmax, dt, model, integrator); + + ASSERT_EQ(refData.size(), static_cast(result.get_num_time_points())); + + for (Eigen::Index irow = 0; irow < result.get_num_time_points(); ++irow) { + double t = refData[static_cast(irow)][0]; + auto rel_tol = 1e-6; + + //test result diverges at damping because of changes, not worth fixing at the moment + if (t > 11.0 && t < 13.0) { + //strong divergence around damping + rel_tol = 0.5; + } + else if (t > 13.0) { + //minor divergence after damping + rel_tol = 1e-2; + } + mio::unused(rel_tol); + + ASSERT_NEAR(t, result.get_times()[irow], 1e-12) << "at row " << irow; + + for (size_t icol = 0; icol < 12; ++icol) { + double ref = refData[static_cast(irow)][icol + 1]; + double actual = result[irow][icol]; + + double tol = rel_tol * ref; + ASSERT_NEAR(ref, actual, tol) << "at row " << irow; + } + } +} + +TEST(TestOdeSirMobility, checkPopulationConservation) +{ + // initialization + double t0 = 0.; + double tmax = 50.; + double dt = 0.1002004008016032; + + double population_per_region = 1061000; + mio::osirmobility::Region num_regions = 4; + + mio::osirmobility::Model model((size_t)num_regions); + + for (auto region = mio::osirmobility::Region(0); region < num_regions; ++region) { + model.populations[{region, mio::AgeGroup(0), mio::osirmobility::InfectionState::Infected}] = 1000; + model.populations[{region, mio::AgeGroup(0), mio::osirmobility::InfectionState::Recovered}] = 1000; + model.populations[{region, mio::AgeGroup(0), mio::osirmobility::InfectionState::Susceptible}] = + population_per_region - + model.populations[{region, mio::AgeGroup(0), mio::osirmobility::InfectionState::Infected}] - + model.populations[{region, mio::AgeGroup(0), mio::osirmobility::InfectionState::Recovered}]; + } + model.parameters.set(2); + model.parameters.set(0.04); + model.parameters.set(1.); + model.parameters.get().get_baseline()(0, 0) = 1.; + model.parameters.get().add_damping(0.6, mio::SimulationTime(12.5)); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(0), 0.5}); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(2), 0.8}); + model.parameters.get().push_back( + {mio::osirmobility::Region(2), mio::osirmobility::Region(0), 0.5}); + model.parameters.get().push_back( + {mio::osirmobility::Region(0), mio::osirmobility::Region(3), 1.0}); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(3), 0.8}); + + auto result = mio::simulate(t0, tmax, dt, model); + double num_persons = 0.0; + for (auto i = 0; i < result.get_last_value().size(); ++i) { + num_persons += result.get_last_value()[i]; + } + EXPECT_NEAR(num_persons, population_per_region * (size_t)num_regions, 1e-8); +} + +TEST(TestOdeSirMobility, check_constraints_parameters) +{ + mio::osirmobility::Region num_regions = 2; + + mio::osirmobility::Model model((size_t)num_regions); + model.parameters.set(6); + model.parameters.set(0.04); + model.parameters.set(1.); + model.parameters.get().get_baseline()(0, 0) = 10.; + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(0), 0.5}); + + // model.check_constraints() combines the functions from population and parameters. + // We only want to test the functions for the parameters defined in parameters.h + ASSERT_EQ(model.parameters.check_constraints(), 0); + + mio::set_log_level(mio::LogLevel::off); + + model.parameters.set(0); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.set(6); + model.parameters.set(10.); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.set(0.04); + model.parameters.set(10.); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.set(1.); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(0), 10.5}); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.get().pop_back(); + model.parameters.get().push_back( + {mio::osirmobility::Region(2), mio::osirmobility::Region(0), 0.5}); + mio::set_log_level(mio::LogLevel::warn); +} + +TEST(TestOdeSirMobility, apply_constraints_parameters) +{ + const double tol_times = 1e-1; + mio::osirmobility::Region num_regions = 2; + + mio::osirmobility::Model model((size_t)num_regions); + model.parameters.set(6); + model.parameters.set(0.04); + model.parameters.set(1.); + model.parameters.get().get_baseline()(0, 0) = 10.; + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(0), 0.5}); + + EXPECT_EQ(model.parameters.apply_constraints(), 0); + + mio::set_log_level(mio::LogLevel::off); + + model.parameters.set(-2.5); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get(), tol_times); + + model.parameters.set(10.); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR(model.parameters.get(), 0.0, 1e-14); + + model.parameters.set(0.04); + model.parameters.set(10.); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR(model.parameters.get(), 0.0, 1e-14); + + model.parameters.set(1.); + model.parameters.get().push_back( + {mio::osirmobility::Region(1), mio::osirmobility::Region(0), 10.5}); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR(std::get(model.parameters.get()[2]), 0.0, 1e-14); + + model.parameters.get().pop_back(); + model.parameters.get().push_back( + {mio::osirmobility::Region(2), mio::osirmobility::Region(0), 0.5}); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + // EXPECT_EQ(model.parameters.get().size(), 2); // 1 by default + 1 added + mio::set_log_level(mio::LogLevel::warn); +} \ No newline at end of file diff --git a/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py b/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py new file mode 100644 index 0000000000..a2c04550df --- /dev/null +++ b/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py @@ -0,0 +1,54 @@ +import os + +import numpy as np +import pandas as pd + +from memilio.epidata import geoModificationGermany as geoger +from memilio.epidata import getDataIntoPandasDataFrame as gd + +def main(): + """! Main program entry.""" + + arg_dict = gd.cli("commuter_official") + + directory = arg_dict['out_folder'].split('/pydata')[0] + directory_mobility = os.path.join(directory, 'mobility/') + directory_population = os.path.join(directory, 'pydata/Germany/') + mobility_file = 'commuter_mobility' + population_file = 'county_current_population' + + mobility_matrix = pd.read_csv( + os.path.join(directory_mobility + mobility_file + '.txt'), + sep=' ', header=None) + + # get county and state IDs + countyIDs = geoger.get_county_ids() + stateIDs = geoger.get_state_ids() + # get state ID to county ID map + stateID_to_countyID = geoger.get_stateid_to_countyids_map() + + # iterate over state_to_county map and replace IDs by numbering 0, ..., n + state_indices = [] + county_indices = [] + for state, counties in stateID_to_countyID.items(): + state_indices.append(stateIDs.index(state)) + county_indices.append( + np.array([countyIDs.index(county) for county in counties])) + + mobility_matrix_nrw = mobility_matrix.loc[county_indices[4], county_indices[4]] + + gd.write_dataframe( + mobility_matrix_nrw, directory_mobility, mobility_file + '_nrw', 'txt', + param_dict={'sep': ' ', 'header': None, 'index': False}) + + population = pd.read_json(os.path.join(directory_population + population_file + '.json')) + population_nrw = population.loc[county_indices[4]] + gd.write_dataframe(population_nrw, directory_population, population_file + '_nrw', 'json') + + + + + +if __name__ == "__main__": + + main() diff --git a/shellscripts/likwid_test.sh b/shellscripts/likwid_test.sh new file mode 100644 index 0000000000..d1e7a8c057 --- /dev/null +++ b/shellscripts/likwid_test.sh @@ -0,0 +1,16 @@ +#!/bin/sh +#SBATCH -N 1 +#SBATCH -n 1 +#SBATCH -c 1 +#SBATCH --exclusive +#SBATCH -t 5-0:00:00 +#SBATCH --output=likwid_test-%A.out +#SBATCH --error=likwid_test-%A.err +#SBATCH --exclude="be-cpu05, be-gpu01" +#SBATCH --job-name=likwid_test + + +cd ../cpp/examples +g++ -O3 -o likwid_test likwid_test.cpp + +srun --cpu-bind=cores --cpus-per-task=1 --cpu-freq=2200000-2200000 likwid-perfctr -C S0:1 -g MEM_DP ./likwid_test \ No newline at end of file diff --git a/tools/basic_reproduction_numbers.py b/tools/basic_reproduction_numbers.py new file mode 100644 index 0000000000..d1ddfb0fc2 --- /dev/null +++ b/tools/basic_reproduction_numbers.py @@ -0,0 +1,40 @@ +from sympy import * +init_printing() + +# 2 Age Groups: +def calculate_eigenvalues_explicit(): + rho=[0.07333, 0.07333, 0.07333] + S=[9990, 10000, 10000] + N=[10000, 10000, 10000] + T_I=[8.097612257, 8.097612257, 8.097612257] + phi=[[7.95/3, 7.95/3, 7.95/3], + [7.95/3, 7.95/3, 7.95/3], + [7.95/3, 7.95/3, 7.95/3]] + value1 = -(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))/(3*(sqrt(-4*(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))**3 + (27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]**2*N[1]**2*N[2]**2) + 2*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**2)/2 + 27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(2*N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(2*N[0]**2*N[1]**2*N[2]**2) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**(1/3)) - (sqrt(-4*(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))**3 + (27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]**2*N[1]**2*N[2]**2) + 2*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**2)/2 + 27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(2*N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(2*N[0]**2*N[1]**2*N[2]**2) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**(1/3)/3 - (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])/(3*N[0]*N[1]*N[2]) + value2 = -(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))/(3*(-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))**3 + (27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]**2*N[1]**2*N[2]**2) + 2*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**2)/2 + 27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(2*N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(2*N[0]**2*N[1]**2*N[2]**2) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**(1/3)) - (-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))**3 + (27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]**2*N[1]**2*N[2]**2) + 2*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**2)/2 + 27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(2*N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(2*N[0]**2*N[1]**2*N[2]**2) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**(1/3)/3 - (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])/(3*N[0]*N[1]*N[2]) + value3 = -(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))/(3*(-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))**3 + (27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]**2*N[1]**2*N[2]**2) + 2*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**2)/2 + 27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(2*N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(2*N[0]**2*N[1]**2*N[2]**2) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**(1/3)) - (-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]*N[1]*N[2]) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**2/(N[0]**2*N[1]**2*N[2]**2))**3 + (27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(N[0]**2*N[1]**2*N[2]**2) + 2*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**2)/2 + 27*(-S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][1]*phi[2][2]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][0]*phi[1][2]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][0]*phi[2][2]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][1]*phi[1][2]*phi[2][0]*rho[0]*rho[1]*rho[2] - S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][0]*phi[2][1]*rho[0]*rho[1]*rho[2] + S[0]*S[1]*S[2]*T_I[0]*T_I[1]*T_I[2]*phi[0][2]*phi[1][1]*phi[2][0]*rho[0]*rho[1]*rho[2])/(2*N[0]*N[1]*N[2]) - 9*(-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])*(N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][1]*phi[2][2]*rho[1]*rho[2] - N[0]*S[1]*S[2]*T_I[1]*T_I[2]*phi[1][2]*phi[2][1]*rho[1]*rho[2] + N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][0]*phi[2][2]*rho[0]*rho[2] - N[1]*S[0]*S[2]*T_I[0]*T_I[2]*phi[0][2]*phi[2][0]*rho[0]*rho[2] + N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][0]*phi[1][1]*rho[0]*rho[1] - N[2]*S[0]*S[1]*T_I[0]*T_I[1]*phi[0][1]*phi[1][0]*rho[0]*rho[1])/(2*N[0]**2*N[1]**2*N[2]**2) + (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])**3/(N[0]**3*N[1]**3*N[2]**3))**(1/3)/3 - (-N[0]*N[1]*S[2]*T_I[2]*phi[2][2]*rho[2] - N[0]*N[2]*S[1]*T_I[1]*phi[1][1]*rho[1] - N[1]*N[2]*S[0]*T_I[0]*phi[0][0]*rho[0])/(3*N[0]*N[1]*N[2]) + print(str(value1) + '\n' + str(value2) + '\n' + str(value3)) + +def calculate_eigenvalues(number_age_groups): + rho= symbols('rho[0]:%d'%number_age_groups) + S = symbols('S0:%d'%number_age_groups) + N = symbols('N0:%d'%number_age_groups) + T_I = symbols('T_I0:%d'%number_age_groups) + phi = [list(symbols('phi' + str(i) + '_0:%d'%number_age_groups)) for i in range(number_age_groups)] + # rho=[0.07333, 0.07333, 0.07333] + # S=[9990, 9990, 9990] + # N=[10000, 10000, 10000] + # T_I=[8.0096875, 8.0096875, 8.2182] + # phi=[[7.95/3, 7.95/3, 7.95/3], + # [7.95/3, 7.95/3, 7.95/3], + # [7.95/3, 7.95/3, 7.95/3]] + next_gen_matrix = zeros(number_age_groups, number_age_groups) + for i in range(number_age_groups): + for j in range(number_age_groups): + next_gen_matrix[i, j] = rho[i] * S[i] * phi[i][j] * T_I[j] / N[j] + p = next_gen_matrix.charpoly(lamda) + print(factor(p.as_expr())) + +if __name__ == '__main__': + calculate_eigenvalues(2) + # calculate_eigenvalues_explicit() \ No newline at end of file diff --git a/tools/plot_mobility_runtimes.py b/tools/plot_mobility_runtimes.py new file mode 100644 index 0000000000..ba2f0f466e --- /dev/null +++ b/tools/plot_mobility_runtimes.py @@ -0,0 +1,155 @@ +import matplotlib.pyplot as plt +import pandas as pd +import numpy as np + +import os + +colors = ['#1f77b4', '#2ca02c', '#ff7f0e'] +linestyles=['-', '--', '-.', ':'] +fontsize_labels = 16 +fontsize_legends = 12 + +models = ['Equation-based model', 'Graph-based model'] + +def plot_runtime(file, name=''): + fig = plt.figure() + df = pd.read_json(file) + + plt.plot(df["Regions"], df["Time"], + linestyle='--', marker='o', linewidth=1.2) + plt.ylim(bottom=df['Time'].min()) + plt.xlim(left=df["Regions"].min()-1, right=df["Regions"].max()+1) + plt.xlabel('Number of regions', fontsize=fontsize_labels) + plt.ylabel('Run time [seconds]', fontsize=fontsize_labels) + plt.yticks(fontsize=fontsize_legends) + plt.xticks(fontsize=fontsize_legends) + plt.grid(True, linestyle='--') + plt.tight_layout() + + plot_dir = os.path.join(os.path.dirname(__file__), '../Plots') + if name is None: + name = os.path.splitext(os.path.basename(file))[0] + plt.savefig(os.path.join(plot_dir, name), bbox_inches='tight', dpi=500) + plt.close() + +def plot_flops(name='number_flops'): + fig, ax = plt.subplots() + + def flops_equation_based(x, eta): + return (4*x**2+22*x+1)/eta + + def flops_graph_based(x, eta): + return (43*x**2+24*x/eta+2)*1 + + x = np.linspace(0, 400, 80) + + + for idx, eta in enumerate([0.05, 0.1, 0.2, 0.5]): + ax.plot(x, flops_equation_based(x, eta), linewidth=1.5, color=colors[0], linestyle=linestyles[idx], label='Model C, $\eta=$'+ str(eta)) + ax.plot(x, flops_graph_based(x, eta), linewidth=1.5, color=colors[1], linestyle=linestyles[idx], label='Model D, $\eta=$'+ str(eta)) + ax.set_ylim(bottom=0.) + ax.set_xlim(left=0., right=400.) + ax.set_xlabel('Number of regions', fontsize=fontsize_labels) + ax.set_ylabel('Number of FLOPs', fontsize=fontsize_labels) + + handles, labels = ax.get_legend_handles_labels() + sorted_handles_labels = sorted(zip(handles, labels), key=lambda x: x[1]) + + sorted_handles, sorted_labels = zip(*sorted_handles_labels) + + plt.tight_layout() + ax.legend(sorted_handles, sorted_labels, fontsize=fontsize_legends) + plt.grid(linestyle='--') + + plot_dir = os.path.join(os.path.dirname(__file__), '../Plots') + plt.savefig(os.path.join(plot_dir, name), bbox_inches='tight', dpi=500) + plt.close() + + + +def compare_runtime_and_flops(files, name=''): + fig, ax1 = plt.subplots() + + for file in files: + df = pd.read_json(file) + + ax1.plot(df["Regions"], df["Time"], + linestyle='--', marker='o', linewidth=1.2, label=file) + + ax1.set_ylim(bottom=0.) + ax1.set_xlim(left=0., right=400.) + ax1.set_xlabel('Number of regions', fontsize=fontsize_labels) + ax1.set_ylabel('Run time [seconds]', fontsize=fontsize_labels) + + ax2 = ax1.twinx() + + def flops_equation_based(x): + return (4*x**2+22*x+1)*200 + + def flops_graph_based(x): + return (43*x**2+240*x+2)*20 + + x = np.linspace(0, 400, 400) + + ax2.plot(x, flops_equation_based(x), linestyle='--', linewidth=1.2) + ax2.plot(x, flops_graph_based(x), linestyle='--', linewidth=1.2) + ax2.set_ylabel('Number of FLOPs', fontsize=fontsize_labels) + ax2.set_ylim(bottom=0.) + + plt.tight_layout() + + plot_dir = os.path.join(os.path.dirname(__file__), '../Plots') + plt.savefig(os.path.join(plot_dir, name), bbox_inches='tight', dpi=500) + plt.close() + +def compare_runtimes(files, name='', title='', models=[]): + merged_df = pd.DataFrame() + i = 0 + for file in files: + df = pd.read_json(file) + df = df.filter(items=['Regions', 'Time']) + # df.drop(thisFilter, inplace=True, axis=1) + df.rename(columns={'Time': models[i]}, inplace=True) + + if merged_df.empty: + merged_df = df + else: + merged_df = pd.merge(merged_df, df, on='Regions', how='outer') + i = i+1 + + merged_df = merged_df.set_index('Regions') + for column in merged_df.columns: + # plt.plot(merged_df['Regions'], column, + # linestyle='--', marker='o', linewidth=1.2) + plt.plot(merged_df.index, merged_df[column], label=column, + linestyle='--', marker='o', linewidth=1.2) + plt.ylim(bottom=0.) + plt.xlim(left=merged_df.index.min()-1, right=merged_df.index.max()+1) + plt.xlabel('Number of regions', fontsize=fontsize_labels) + plt.ylabel('Run time [seconds]', fontsize=fontsize_labels) + plt.yticks(fontsize=fontsize_legends) + plt.xticks(fontsize=fontsize_legends) + plt.grid(True, linestyle='--') + plt.legend() + plt.title(title) + plt.tight_layout() + + plot_dir = os.path.join(os.path.dirname(__file__), '../Plots') + plt.savefig(os.path.join(plot_dir, name), bbox_inches='tight', dpi=500) + plt.close() + +if __name__ == "__main__": + result_dir = os.path.join(os.path.dirname(__file__), '../results') + + result_equationbased_euler = os.path.join(result_dir, 'timing_equationbased_euler.json') + result_equationbased_noage_euler = os.path.join(result_dir, 'timing_equationbased_noage_euler.json') + result_graphbased_euler = os.path.join(result_dir, 'timing_graphbased_euler.json') + result_graphbased_noage_euler = os.path.join(result_dir, 'timing_graphbased_noage_euler.json') + + results_euler = [result_equationbased_euler, result_graphbased_euler] + results_euler_noage = [result_equationbased_noage_euler, result_graphbased_noage_euler] + + # compare_runtimes(results_euler, name='compare_runtimes_euler', models=models) + # compare_runtimes(results_euler_noage, name='compare_runtimes_euler_noage', models=models) + # compare_runtime_and_flops(results_euler_noage, 'compare_runtimes_and_flops') + plot_flops() diff --git a/tools/plot_results_mobility.py b/tools/plot_results_mobility.py new file mode 100644 index 0000000000..f8121c5676 --- /dev/null +++ b/tools/plot_results_mobility.py @@ -0,0 +1,385 @@ + +import datetime as dt +import os.path +import h5py + +import numpy as np +import pandas as pd + +import geopandas +from matplotlib.gridspec import GridSpec + +from memilio.epidata import geoModificationGermany as geoger + +import memilio.epidata.getPopulationData as gpd +from memilio.epidata import getDataIntoPandasDataFrame as gd +import memilio.plot.plotMap as pm + +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors +from matplotlib.ticker import FormatStrFormatter + +compartments = {'Susceptible': 0, + 'Exposed': 1, + 'Infected': 2, + 'Recovered': 3} + + +def plot_map_nrw(data: pd.DataFrame, + scale_colors: np.array([0, 1]), + legend: list = [], + title: str = '', + plot_colorbar: bool = True, + output_path: str = '', + fig_name: str = 'customPlot', + dpi: int = 300, + outercolor='white', + log_scale=False, + cmap='viridis'): + """! Plots region-specific information onto a interactive html map and + returning svg and png image. Allows the comparisons of a variable list of + data sets. + + @param[in] data Data to be plotted. First column must contain regional + specifier, following columns will be plotted for comparison. + @param[in] scale_colors Array of min-max-values to scale colorbar. + @param[in] legend List subtitles for different columns. Can be list of + empty strings. + @param[in] title Title of the plot. + @param[in] plot_colorbar Defines if a colorbar will be plotted. + @param[in] output_path Output path for the figure. + @param[in] fig_name Name of the figure created. + @param[in] dpi Dots-per-inch value for the exported figure. + @param[in] outercolor Background color of the plot image. + @param[in] log_scale Defines if the colorbar is plotted in log scale. + """ + region_classifier = data.columns[0] + region_data = data[region_classifier].to_numpy().astype(int) + + data_columns = data.columns[1:] + # Read and filter map data. + if np.isin(region_data, geoger.get_county_ids()).all(): + try: + map_data = geopandas.read_file( + os.path.join( + os.getcwd(), + 'tools/vg2500_12-31.utm32s.shape/vg2500/VG2500_KRS.shp')) + if '16056' in map_data.ARS.values: + map_data = pm.merge_eisenach(map_data) + # Remove information for plot. + map_data = map_data[['ARS', 'GEN', 'NUTS', 'geometry']] + # Use string values as in shape data file. + data[region_classifier] = data[region_classifier].astype( + 'str').str.zfill(5) + except FileNotFoundError: + pm.print_manual_download( + 'Georeferenzierung: UTM32s, Format: shape (ZIP, 5 MB)', + 'https://gdz.bkg.bund.de/index.php/default/verwaltungsgebiete-1-2-500-000-stand-31-12-vg2500-12-31.html') + else: + raise gd.DataError('Provide shape files regions to be plotted.') + + # Remove regions that are not input data table. + map_data = map_data[map_data.ARS.isin(data[region_classifier])] + + data['new_index'] = map_data.index.array + data = data.set_index('new_index') + + map_data[data_columns] = data.loc[:, data_columns] + + for i in range(len(data_columns)): + if legend[i] == '': + fname = 'data_column_' + str(i) + else: + fname = str(legend[i].replace(' ', '_')) + pm.save_interactive(data[data_columns[i]], os.path.join( + output_path, fname) + '.html', map_data, scale_colors) + + fig = plt.figure(figsize=(3.5 * len(data_columns), 3), facecolor=outercolor) + # Use n+2 many columns (1: legend + 2: empty space + 3-n: data sets) and + # n+2 rows where the top row is used for a potential title, the second row + # for the content and all other rows have height zero. + height_ratios = [0.05, 1] + # if len(data_columns) > 1: + # height_ratios = height_ratios + [ + # 0.0 for i in range(len(data_columns)-1)] + if plot_colorbar: + gs = GridSpec( + 2, len(data_columns)+2, figure=fig, + width_ratios=[1 for i in range(len(data_columns))]+[0.1, 0.2], + height_ratios=height_ratios) + else: + gs = GridSpec( + 2, len(data_columns), figure=fig, + width_ratios=[1 for i in range(len(data_columns))], + height_ratios=height_ratios) + + # Use top row for title. + tax = fig.add_subplot(gs[0, :]) + tax.set_axis_off() + tax.set_title(title, fontsize=16) + if plot_colorbar: + # Prepare colorbar. + cax = fig.add_subplot(gs[1, -2]) + + else: + cax = None + + if log_scale: + norm = mcolors.LogNorm(vmin=scale_colors[0], vmax=scale_colors[1]) + else: + norm = mcolors.TwoSlopeNorm(vmin=scale_colors[0], vmax=scale_colors[1], vcenter=0) + + for i in range(len(data_columns)): + + ax = fig.add_subplot(gs[1, i]) + if log_scale: + map_data.plot(data_columns[i], ax=ax, legend=False, + norm=norm, cmap=cmap) + + elif cax is not None: + map_data.plot(data_columns[i], ax=ax, cax=cax, legend=True, + vmin=scale_colors[0], vmax=scale_colors[1], cmap=cmap) + else: + # Do not plot colorbar. + map_data.plot(data_columns[i], ax=ax, legend=False, + vmin=scale_colors[0], vmax=scale_colors[1], cmap=cmap) + + ax.set_title(legend[i], fontsize=9) + ax.set_axis_off() + + if plot_colorbar: + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + cbar = fig.colorbar(sm, cax=cax) + # cbar.set_ticks([scale_colors[0], scale_colors[1]]) + # cbar.set_ticklabels([f'{scale_colors[0]:.3e}', f'{scale_colors[1]:.3e}'], fontsize=7) + + plt.subplots_adjust(bottom=0.1, left=0.1) + plt.tight_layout() + plt.savefig(os.path.join(output_path, fig_name + '.png'), dpi=dpi) + plt.close() + +def plot_maps(files, output_dir, legend, name=''): + + min_val = 10e-6 + max_val = 0.4 + + cmap = 'viridis' + norm = mcolors.LogNorm(vmin=min_val, vmax=max_val) + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + cbar_fig, cax = plt.subplots(figsize=(12, 1)) + cbar = plt.colorbar(sm, orientation='horizontal', cax=cax) + cbar.ax.tick_params(labelsize=10) + cbar.set_ticks([min_val, max_val]) + cbar.ax.xaxis.set_major_formatter(FormatStrFormatter('%.5f')) + plt.subplots_adjust(bottom=0.3) + plt.savefig(os.path.join(output_dir, 'colorbar.png'), dpi=300) + plt.close() + + for date in range(10, 101, 20): + dfs_all = extract_nrw_data_and_combine(files=files, date=date) + + plot_map_nrw( + dfs_all, scale_colors=[min_val, max_val], + legend=legend, + title='NRW - Simulation Day '+str(date), plot_colorbar=False, + output_path=output_dir, + fig_name=name+str(date), dpi=900, + outercolor='white', + log_scale=True) + +def plot_difference_2D(files, output_dir): + fig = plt.figure() + + df_dif = pd.DataFrame(columns=['Time', 'difference', 'absolute value']) + dates = [i for i in range(100)] + df_dif['Time'] = dates + df_dif.set_index('Time', inplace=True) + + total_population = 18190422. + + for date in range(100): + dfs_all = extract_nrw_data_and_combine(files=files, date=date, relative=False) + df_dif.loc[date,'difference'] = (dfs_all[dfs_all.columns[1]] - dfs_all[dfs_all.columns[2]]).sum() / total_population + df_dif.loc[date,'absolute value'] = abs(dfs_all[dfs_all.columns[1]] - dfs_all[dfs_all.columns[2]]).sum() / total_population + + # df_dif.set_index('Time', inplace=True) + df_dif['difference'].plot(label='Difference') + df_dif['absolute value'].plot(label='Difference in absolute value') + plt.legend() + plt.grid(linestyle='--') + plt.savefig(os.path.join(output_dir, 'difference2D.png')) + plt.close() + + +def plot_difference(files, output_dir): + fig = plt.figure() + + df_dif = pd.DataFrame(columns=['Region']) + + for date in range(60, 101, 10): + dfs_all = extract_nrw_data_and_combine(files=files, date=date, relative=True) + df_dif['Region'] = dfs_all['Region'] + + + df_dif['Count (rel)' + str(date)] = dfs_all[dfs_all.columns[1]] - dfs_all[dfs_all.columns[2]] + + # df_dif = df_dif[df_dif['Count (rel)' + str(date)] > 0] + + min_val = df_dif.drop(columns=['Region']).min().min() + max_val = df_dif.drop(columns=['Region']).max().max() + maximum_abs = abs(max([min_val, max_val], key=abs)) + + plot_map_nrw( + df_dif, scale_colors=[-maximum_abs, maximum_abs], + legend=['Day ' + str(date) for date in range(60, 101, 10)], + title='Difference between ODE and graph-based model', plot_colorbar=True, + output_path=output_dir, + fig_name="difference", dpi=900, + outercolor='white', + log_scale=False, + cmap='seismic') + +def extract_nrw_data_and_combine(files, date, relative=True): + age_groups = {0: '0-4', 1: '5-14', 2: '15-34', + 3: '35-59', 4: '60-79', 5: '80+'} + filter_age = None + + i = 0 + for file in files.values(): + model_type = os.path.basename(file).split('_')[0] + if model_type == 'ode': # result file of equation-based model has to be first + df = pm.extract_data( + file, region_spec=None, column=None, date=date, + filters={'Group': filter_age, 'InfectionState': [2]}, + output='matrix', + file_format=file_format) + df['Group'] = df.Group.str.extract('(\d+)') + df['Group'] = df['Group'].apply(pd.to_numeric, errors='coerce') + # df['Region'] = df['Group'] + df['Region'] = (df['Group']-1) // len(age_groups) + df = df.groupby(['Region'], as_index=False).agg({'Count': "sum"}) + + ids = geoger.get_county_ids() + ids = [id for id in ids if str(id).startswith('5')] + + if len(ids) != len(df): + raise gd.DataError("Data is not compatible with number of NRW counties.") + + df['Region'] = ids + else: + df = pm.extract_data( + file, region_spec=None, column=None, date=date, + filters={'Group': filter_age, 'InfectionState': [2]}, + file_format=file_format) + + df = df.apply(pd.to_numeric, errors='coerce') + + if relative: + + try: + population = pd.read_json( + 'data/pydata/Germany/county_current_population.json') + # pandas>1.5 raise FileNotFoundError instead of ValueError + except (ValueError, FileNotFoundError): + print("Population data was not found. Download it from the internet.") + population = gpd.get_population_data( + read_data=False, file_format=file_format, + out_folder='data/pydata/Germany/', no_raw=True, + merge_eisenach=True) + + # For fitting of different age groups we need format ">X". + age_group_values = list(age_groups.values()) + age_group_values[-1] = age_group_values[-1].replace('80+', '>79') + # scale data + df = pm.scale_dataframe_relative(df, age_group_values, population) + + if i == 0: + dfs_all = pd.DataFrame(df.iloc[:, 0]) + + dfs_all[df.columns[-1] + ' ' + str(i)] = df[df.columns[-1]] + i += 1 + + return dfs_all + + +def plot_total_compartment(files, output_dir, legend, compartment = 'Infected', name='', ax=None, print_legend=True): + + colors = ['#1f77b4', '#2ca02c', '#ff7f0e'] + file_idx = 0 + if ax is None: + fig, ax = plt.subplots() + ax.grid(True, linestyle='--') + for file in files.values(): + model_type = os.path.basename(file).split('_')[0] + # Load data. + h5file = h5py.File(file + '.h5', 'r') + if model_type=='ode': + dates = h5file['1']['Time'][:] + data = h5file['1']['Total'][:,compartments[compartment]] + ax.plot(dates, data, label=legend[file_idx], linewidth=2, color=colors[file_idx]) + ax.set_ylim(bottom=0.) + ax.set_xlim(left=0., right=dates.max()+1) + else: + df = pd.DataFrame() + regions = list(h5file.keys()) + for i in range(len(regions)): + df['Region'+str(i)] = h5file[regions[i]]['Total'][:, compartments[compartment]] + df['Total'] = df.sum(axis=1) + df['Time'] = h5file[regions[0]]['Time'][:] # hardcoded + ax.plot(df['Time'], df['Total'], label=legend[file_idx], linewidth=1.2, color=colors[file_idx]) + ax.set_ylim(bottom=0.) + ax.set_xlim(left=0., right=df['Time'].max()+1) + # ax.legend() + + file_idx = file_idx+1 + + plt.tight_layout() + if print_legend: + plt.legend() + plt.savefig(os.path.join(output_dir, name + '.png'), dpi=300) + + return ax + +def compare_compartments(files, output_dir, legend): + fig, axs = plt.subplots( + 2, 2, sharex='all') + axs = axs.flatten() + for i, compartment in enumerate(compartments.keys()): + plot_total_compartment(files=files, output_dir=output_dir, legend=legend, compartment=compartment, ax=axs[i], print_legend=False) + plt.tight_layout() + plt.subplots_adjust(bottom=0.15) + lines, labels = axs[0].get_legend_handles_labels() + fig.legend(lines, labels, ncol=len(models), loc='center', + fontsize=10, bbox_to_anchor=(0.5, 0.05)) + plt.savefig(os.path.join(output_dir, 'compare_all_compartments.png'), dpi=300) + plt.close() + + +if __name__ == '__main__': + + results_euler = {'Model B': 'cpp/build/ode_result_paper_nrw_euler', + 'Model C': 'cpp/build/ode_result_nrw_euler', + 'Model D': 'cpp/build/graph_result_nrw_euler'} + results_adaptive = {'Model B': 'cpp/build/ode_result_paper_nrw_adaptive', + 'Model C': 'cpp/build/ode_result_nrw_adaptive', + 'Model D': 'cpp/build/graph_result_nrw_adaptive'} + files_compare_solver = {'Data set 1': 'cpp/build/ode_result_nrw_euler', + 'Data set 2': 'cpp/build/ode_result_nrw_adaptive', + 'Data set 3': 'cpp/build/graph_result_nrw_euler', + 'Data set 4': 'cpp/build/graph_result_nrw_adaptive'} + file_format = 'h5' + + models = ['Model B', + 'Model C', + 'Model D'] + + plot_dir = os.path.join(os.path.dirname(__file__), '../Plots') + + plot_maps(files=results_adaptive, output_dir=plot_dir, legend=models, name='NRWAdaptiveDay') + # plot_difference_2D(files={key: value for key, value in results_adaptive.items() if key in { + # 'Model C', 'Model D'}}, output_dir=plot_dir) + compare_compartments(files=results_adaptive, output_dir=plot_dir, legend=models) + # plot_total_compartment(files=results_adaptive, output_dir=plot_dir, legend=models, + # compartment='Infected', name='infectives', title='Total infectives') diff --git a/tools/plot_results_mobilitymodels.py b/tools/plot_results_mobilitymodels.py new file mode 100644 index 0000000000..0426ca7747 --- /dev/null +++ b/tools/plot_results_mobilitymodels.py @@ -0,0 +1,205 @@ +import h5py +import os +import matplotlib.pyplot as plt + +import memilio.epidata.getDataIntoPandasDataFrame as gd + +# Define compartments. +secir_dict = {0: 'Susceptible', 1: 'Exposed', 2: 'Infected', 3: 'Recovered'} + +# Define color and style to be used while plotting for different models to +# make plots consistent. +color_dict = {0: '#1f77b4', + 1: '#2ca02c' + } +linestyle_dict = {"ODE SI": 'dashed', + "ODE Improved": 'dotted', + "Graph": 'dashdot' + } + + +def compare_all_compartments( + files, + legendplot, + filename_plot="compare_compartments"): + + fig, axs = plt.subplots( + 2, 2, sharex='all', num=filename_plot, tight_layout=False) + + # Add simulation results to plot. + for file in range(len(files)): + # Load data. + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + number_regions = len(list(h5file.keys())) + for region in range(number_regions): + if (len(list(h5file[list(h5file.keys())[region]].keys())) > 3): + data = h5file[list(h5file.keys())[region]] + dates = data['Time'][:] + + number_regions = len( + list(h5file[list(h5file.keys())[region]].keys())) - 2 + for region in range(number_regions): + total = data['Group' + str(region + 1)][:, :] + if (total.shape[1] != 4): + raise gd.DataError( + "Expected a different number of compartments.") + # Plot result. + if legendplot[file] in linestyle_dict: + for i in range(4): + axs[int(i / 2), + i % 2].plot(dates, + total[:, + i], + label=legendplot[file] + + " Region " + str(region), + linewidth=1.2, + linestyle=linestyle_dict[legendplot[file]], + color=color_dict[region]) + else: + for i in range(4): + axs[int(i / 2), i % 2].plot(dates, total[:, i], + label=legendplot[file], linewidth=1.2) + else: + data = h5file[list(h5file.keys())[region]] + dates = data['Time'][:] + # As there should be only one Group, total is the simulation + # result. + total = data['Total'][:, :] + if (total.shape[1] != 4): + raise gd.DataError( + "Expected a different number of compartments.") + # Plot result. + if legendplot[file] in linestyle_dict: + for i in range(4): + axs[int(i / 2), + i % 2].plot(dates, + total[:, + i], + label=legendplot[file] + + " Region " + str(region), + linewidth=1.2, + linestyle=linestyle_dict[legendplot[file]], + color=color_dict[region]) + else: + for i in range(4): + axs[int(i / 2), i % 2].plot(dates, total[:, i], + label=legendplot[file], linewidth=1.2) + h5file.close() + + # Define some characteristics of the plot. + for i in range(4): + axs[int(i / 2), i % 2].set_title(secir_dict[i], fontsize=8) + axs[int(i / 2), i % 2].set_xlim(left=0, right=dates[-1]) + axs[int(i / 2), i % 2].grid(True, linestyle='--') + axs[int(i / 2), i % 2].tick_params(axis='y', labelsize=7) + axs[int(i / 2), i % 2].tick_params(axis='x', labelsize=7) + # axs[int(i/2), i % 2].xaxis.set_ticks(np.arange(0, dates[-1]+1, 5)) + + fig.supxlabel('Time (in days)', fontsize=9) + + lines, labels = axs[0, 0].get_legend_handles_labels() + lgd = fig.legend(lines, labels, ncol=len(legendplot), loc='outside lower center', + fontsize=10, bbox_to_anchor=(0.5, - 0.06), bbox_transform=fig.transFigure) + + plt.tight_layout(pad=0, w_pad=0.5, h_pad=0.1) + plt.subplots_adjust(bottom=0.09) + + # Save result. + if not os.path.isdir('Plots'): + os.makedirs('Plots') + fig.savefig('Plots/' + filename_plot + '.png', + bbox_extra_artists=(lgd,), bbox_inches='tight', dpi=500) + + +def plot_new_infections( + files, + ylim, + legendplot, + filename_plot="compare_new_infections"): + + plt.figure(filename_plot) + + # Add simulation results to plot. + for file in range(len(files)): + # Load data. + h5file = h5py.File(str(files[file]) + '.h5', 'r') + + number_regions = len(list(h5file.keys())) + for region in range(number_regions): + if (len(list(h5file[list(h5file.keys())[region]].keys())) > 3): + data = h5file[list(h5file.keys())[region]] + dates = data['Time'][:] + + number_regions = len( + list(h5file[list(h5file.keys())[region]].keys())) - 2 + for region_ in range(number_regions): + total = data['Group' + str(region_ + 1)][:, :] + if (total.shape[1] != 4): + raise gd.DataError( + "Expected a different number of compartments.") + incidence = (total[:-1, 0] - total[1:, 0] + ) / (dates[1:] - dates[:-1]) + # Plot result. + if legendplot[file] in linestyle_dict: + plt.plot(dates[1:], + incidence, + linewidth=1.2, + linestyle=linestyle_dict[legendplot[file]], + color=color_dict[region_]) + else: + plt.plot(dates[1:], incidence, linewidth=1.2) + else: + data = h5file[list(h5file.keys())[region]] + dates = data['Time'][:] + # As there should be only one Group, total is the simulation + # result. + total = data['Total'][:, :] + if (total.shape[1] != 4): + raise gd.DataError( + "Expected a different number of compartments.") + incidence = (total[:-1, 0] - total[1:, 0]) / \ + (dates[1:] - dates[:-1]) + # Plot result. + if legendplot[file] in linestyle_dict: + plt.plot(dates[1:], + incidence, + linewidth=1.2, + linestyle=linestyle_dict[legendplot[file]], + color=color_dict[region]) + else: + plt.plot(dates[1:], incidence, linewidth=1.2) + + h5file.close() + + plt.xlabel('Time (in days)', fontsize=16) + # plt.xticks(np.arange(0, dates[-1]+1, 5)) + plt.yticks(fontsize=9) + plt.ylabel('New infections per day', fontsize=14) + plt.ylim(bottom=0, top=ylim) + plt.xlim(left=0, right=dates[-1]) + plt.legend(legendplot, fontsize=14, framealpha=0.5) + plt.grid(True, linestyle='--') + plt.tight_layout() + + # Save result. + if not os.path.isdir('Plots'): + os.makedirs('Plots') + plt.savefig( + 'Plots/' + + filename_plot + + '.png', + bbox_inches='tight', + dpi=500) + + +if __name__ == '__main__': + data_dir = os.path.join(os.path.dirname(__file__), "..", "cpp", "build") + plot_new_infections([os.path.join(data_dir, "ode_result_standard"), + os.path.join(data_dir, "ode_result_improved"), + os.path.join(data_dir, "graph_result")], + 2e3, legendplot=list(["ODE SI", "ODE Improved", "Graph"])) + compare_all_compartments([os.path.join(data_dir, "ode_result_standard"), + os.path.join(data_dir, "ode_result_improved"), + os.path.join(data_dir, "graph_result")], + legendplot=list(["ODE SI", "ODE Improved", "Graph"]))