diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 45a3d67773..4ebd70bc2e 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -161,6 +161,8 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/ode_seir) add_subdirectory(models/ode_seair) add_subdirectory(models/ode_sir) + add_subdirectory(models/ode_metapop) + add_subdirectory(models/ode_metapop_wang) 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 57e7033d5f..60b4ce4e9d 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -30,6 +30,41 @@ 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_metapop_wang_nrw examples_thesis/ode_metapop_wang_nrw.cpp) +target_link_libraries(ode_metapop_wang_nrw PRIVATE memilio ode_metapop_wang) +target_compile_options(ode_metapop_wang_nrw PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +if (MEMILIO_ENABLE_OPENMP) + add_executable(ode_metapop_timing examples_thesis/ode_metapop_timing.cpp) + target_link_libraries(ode_metapop_timing PRIVATE memilio ode_metapop likwid stdc++) + target_compile_options(ode_metapop_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + target_compile_definitions(ode_metapop_timing PRIVATE "-DLIKWID_PERFMON") +endif() + +add_executable(ode_metapop_steps examples_thesis/ode_metapop_steps.cpp) +target_link_libraries(ode_metapop_steps PRIVATE memilio ode_metapop) +target_compile_options(ode_metapop_steps PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(ode_metapop_seir_example ode_seir_metapop.cpp) +target_link_libraries(ode_metapop_seir_example PRIVATE memilio ode_metapop) +target_compile_options(ode_metapop_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(basic_reproduction_number_modela examples_thesis/basic_reproduction_number_modela.cpp) +target_link_libraries(basic_reproduction_number_modela PRIVATE memilio ode_seir) +target_compile_options(basic_reproduction_number_modela PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(basic_reproduction_number_modelb examples_thesis/basic_reproduction_number_modelb.cpp) +target_link_libraries(basic_reproduction_number_modelb PRIVATE memilio ode_metapop_wang) +target_compile_options(basic_reproduction_number_modelb PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(basic_reproduction_number_modelc examples_thesis/basic_reproduction_number_modelc.cpp) +target_link_libraries(basic_reproduction_number_modelc PRIVATE memilio ode_metapop) +target_compile_options(basic_reproduction_number_modelc PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(ode_metapop_nrw examples_thesis/ode_metapop_nrw.cpp) +target_link_libraries(ode_metapop_nrw PRIVATE memilio ode_metapop) +target_compile_options(ode_metapop_nrw 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}) @@ -92,6 +127,21 @@ add_executable(ode_secir_graph_example ode_secir_graph.cpp) target_link_libraries(ode_secir_graph_example PRIVATE memilio ode_secir) target_compile_options(ode_secir_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(graph_nrw examples_thesis/graph_nrw.cpp) +target_link_libraries(graph_nrw PRIVATE memilio ode_seir) +target_compile_options(graph_nrw PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(graph_steps examples_thesis/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 examples_thesis/graph_timing.cpp) + target_link_libraries(graph_timing PRIVATE memilio ode_seir likwid stdc++) + target_compile_options(graph_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + target_compile_definitions(graph_timing PRIVATE "-DLIKWID_PERFMON") +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/examples_thesis/basic_reproduction_number_modela.cpp b/cpp/examples/examples_thesis/basic_reproduction_number_modela.cpp new file mode 100644 index 0000000000..aa58d7cbea --- /dev/null +++ b/cpp/examples/examples_thesis/basic_reproduction_number_modela.cpp @@ -0,0 +1,64 @@ +#include "models/ode_seir/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 calculate_basic_reproduction_number(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; + } + + 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 << "\"Model A\": " << basic_reproduction_number << ", " << std::endl; +} + +int main() +{ + const ScalarType tmax = 1.; + size_t num_regions = 1; + + std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl; + + calculate_basic_reproduction_number(num_regions, tmax); + return 0; +} \ No newline at end of file diff --git a/cpp/examples/examples_thesis/basic_reproduction_number_modelb.cpp b/cpp/examples/examples_thesis/basic_reproduction_number_modelb.cpp new file mode 100644 index 0000000000..401e460b0e --- /dev/null +++ b/cpp/examples/examples_thesis/basic_reproduction_number_modelb.cpp @@ -0,0 +1,78 @@ + +#include "models/ode_metapop_wang/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 calculate_basic_reproduction_number(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::oseirmetapopwang::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::oseirmetapopwang::Region(i), mio::AgeGroup(j), + mio::oseirmetapopwang::InfectionState::Susceptible}] = 10000; + } + } + + 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 << "\"Model B\": " << basic_reproduction_number << "}" << std::endl; +} + +int main() +{ + const ScalarType tmax = 1.; + size_t num_regions = 1; + + std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl; + + calculate_basic_reproduction_number(num_regions, tmax); + return 0; +} \ No newline at end of file diff --git a/cpp/examples/examples_thesis/basic_reproduction_number_modelc.cpp b/cpp/examples/examples_thesis/basic_reproduction_number_modelc.cpp new file mode 100644 index 0000000000..683a8500c9 --- /dev/null +++ b/cpp/examples/examples_thesis/basic_reproduction_number_modelc.cpp @@ -0,0 +1,98 @@ +#include "models/ode_metapop/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 calculate_basic_reproduction_number(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::oseirmetapop::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::oseirmetapop::Region(i), mio::AgeGroup(j), + mio::oseirmetapop::InfectionState::Susceptible}] = 10000; + } + } + + 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::oseirmetapop::InfectionState::Count; state++) { + population_n += population[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age), + mio::oseirmetapop::InfectionState(state)}]; + } + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] += population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + population_after_commuting[{mio::oseirmetapop::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 = 1; + + std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl; + + calculate_basic_reproduction_number(num_regions, tmax); + return 0; +} \ No newline at end of file diff --git a/cpp/examples/examples_thesis/graph_nrw.cpp b/cpp/examples/examples_thesis/graph_nrw.cpp new file mode 100644 index 0000000000..ead6798cc5 --- /dev/null +++ b/cpp/examples/examples_thesis/graph_nrw.cpp @@ -0,0 +1,209 @@ + +#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) +{ + 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 / "Germany" / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + BOOST_OUTCOME_TRY( + auto&& minimum, + mio::read_mobility_plain( + (data_dir / "Germany" / "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(); +} + +/** + * 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) +{ + 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(); +} + +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 / "Germany" / "pydata" / "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 run(const fs::path& data_dir, double t0, double tmax, double dt) +{ + mio::set_log_level(mio::LogLevel::off); + // global parameters + const int num_age_groups = 6; + + mio::oseir::Parameters params(num_age_groups); + + BOOST_OUTCOME_TRY(set_covid_parameters(params)); + + // set contact matrix + BOOST_OUTCOME_TRY(set_contact_matrices(data_dir, params)); + + // 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 / "Germany" / "pydata" / "county_current_population_nrw.json").string(), true, + true)); + + 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 / "Germany" / "mobility" / "commuter_mobility_2022_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)); + } + } + + 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.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/examples_thesis/graph_steps.cpp b/cpp/examples/examples_thesis/graph_steps.cpp new file mode 100644 index 0000000000..48f0c6cfeb --- /dev/null +++ b/cpp/examples/examples_thesis/graph_steps.cpp @@ -0,0 +1,164 @@ + +#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(); + + std::cout << " \"Regions\": " << number_regions << "," << std::endl; + std::cout << " \"Steps Hotspot\": " << result_graph.nodes()[0].property.get_result().get_num_time_points() - 1 + << "," << std::endl; + std::cout << " \"Steps other Regions\": " + << result_graph.nodes()[99].property.get_result().get_num_time_points() - 1; +} + +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 << "}," << std::endl; + + return 0; +} diff --git a/cpp/examples/examples_thesis/graph_timing.cpp b/cpp/examples/examples_thesis/graph_timing.cpp new file mode 100644 index 0000000000..e32abde8b3 --- /dev/null +++ b/cpp/examples/examples_thesis/graph_timing.cpp @@ -0,0 +1,181 @@ + +#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 = false; + +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}] = 60000; + } + for (auto& node : nodes) { + node.parameters = params; + node.populations = population; + } + 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/examples_thesis/ode_metapop_nrw.cpp b/cpp/examples/examples_thesis/ode_metapop_nrw.cpp new file mode 100644 index 0000000000..7f4db879b3 --- /dev/null +++ b/cpp/examples/examples_thesis/ode_metapop_nrw.cpp @@ -0,0 +1,201 @@ + +#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_metapop/infection_state.h" +#include "models/ode_metapop/model.h" +#include "models/ode_metapop/parameters.h" +#include "models/ode_metapop/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::oseirmetapop::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::oseirmetapop::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 / "Germany" / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + BOOST_OUTCOME_TRY( + auto&& minimum, + mio::read_mobility_plain( + (data_dir / "Germany" / "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::oseirmetapop::Model& model, const fs::path& data_dir) +{ + BOOST_OUTCOME_TRY( + auto&& node_ids, + mio::get_node_ids((data_dir / "Germany" / "pydata" / "county_current_population_nrw.json").string(), true, + true)); + + BOOST_OUTCOME_TRY(const auto&& population_data, + mio::read_population_data( + (data_dir / "Germany" / "pydata" / "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::oseirmetapop::Region(region_idx), mio::AgeGroup(age), + mio::oseirmetapop::InfectionState::Susceptible}] = + entry.population[mio::AgeGroup(age)]; + } + } + } + + printf("Setting population data successful.\n"); + return mio::success(); +} + +template +mio::IOResult set_mobility_weights(mio::oseirmetapop::Model& model, const fs::path& data_dir) +{ + size_t number_regions = (size_t)model.parameters.get_num_regions(); + if (number_regions == 1) { + model.set_commuting_strengths(); + + return mio::success(); + } + else { + // mobility between nodes + BOOST_OUTCOME_TRY( + auto&& mobility_data_commuter, + mio::read_mobility_plain((data_dir / "Germany" / "mobility" / "commuter_mobility_2022_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::oseirmetapop::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.set_commuting_strengths(mobility_data_commuter); + + printf("Setting mobility weights successful.\n"); + return mio::success(); + } +} + +template +mio::IOResult set_parameters_and_population(mio::oseirmetapop::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::oseirmetapop::Region(27), mio::AgeGroup(4), mio::oseirmetapop::InfectionState::Susceptible}] -= + 100; + populations[{mio::oseirmetapop::Region(27), mio::AgeGroup(4), mio::oseirmetapop::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::oseirmetapop::Model model(number_regions, number_age_groups); + auto result_prepare_simulation = set_parameters_and_population(model, data_dir); + + 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.h5"); +} diff --git a/cpp/examples/examples_thesis/ode_metapop_steps.cpp b/cpp/examples/examples_thesis/ode_metapop_steps.cpp new file mode 100644 index 0000000000..7c6cb369a0 --- /dev/null +++ b/cpp/examples/examples_thesis/ode_metapop_steps.cpp @@ -0,0 +1,159 @@ + +#include "memilio/compartments/simulation.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_metapop/infection_state.h" +#include "models/ode_metapop/model.h" +#include "models/ode_metapop/parameters.h" +#include "models/ode_metapop/regions.h" + +#include + +bool age_groups = true; + +template +void set_contact_matrix(mio::oseirmetapop::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::oseirmetapop::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::oseirmetapop::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::oseirmetapop::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::oseirmetapop::Region(i), mio::AgeGroup(j), + mio::oseirmetapop::InfectionState::Susceptible}] = 10000; + } + } + model.populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed}] += + 100; + model.populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), + mio::oseirmetapop::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::oseirmetapop::InfectionState::Count; state++) { + population_n += populations[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age), + mio::oseirmetapop::InfectionState(state)}]; + } + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] += population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + population_after_commuting[{mio::oseirmetapop::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::oseirmetapop::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/examples_thesis/ode_metapop_timing.cpp b/cpp/examples/examples_thesis/ode_metapop_timing.cpp new file mode 100644 index 0000000000..7d09b72721 --- /dev/null +++ b/cpp/examples/examples_thesis/ode_metapop_timing.cpp @@ -0,0 +1,189 @@ +/* +* 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_metapop/infection_state.h" +#include "models/ode_metapop/model.h" +#include "models/ode_metapop/parameters.h" +#include "models/ode_metapop/regions.h" + +#include + +bool age_groups = false; + +template +void set_contact_matrix(mio::oseirmetapop::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::oseirmetapop::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::oseirmetapop::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::oseirmetapop::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::oseirmetapop::Region(i), mio::AgeGroup(j), + mio::oseirmetapop::InfectionState::Susceptible}] = 60000; + } + } + populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed}] += 100; + populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::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::oseirmetapop::InfectionState::Count; state++) { + population_n += populations[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age), + mio::oseirmetapop::InfectionState(state)}]; + } + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] += population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths[0].get_baseline()(region_n, region_m) * population_n; + population_after_commuting[{mio::oseirmetapop::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::oseirmetapop::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/examples_thesis/ode_metapop_wang_nrw.cpp b/cpp/examples/examples_thesis/ode_metapop_wang_nrw.cpp new file mode 100644 index 0000000000..62f200fcb8 --- /dev/null +++ b/cpp/examples/examples_thesis/ode_metapop_wang_nrw.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_metapop_wang/infection_state.h" +#include "models/ode_metapop_wang/model.h" +#include "models/ode_metapop_wang/parameters.h" +#include "models/ode_metapop_wang/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::oseirmetapopwang::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::oseirmetapopwang::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::oseirmetapopwang::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::oseirmetapopwang::Region(region_idx), mio::AgeGroup(age), + mio::oseirmetapopwang::InfectionState::Susceptible}] = + entry.population[mio::AgeGroup(age)]; + } + } + } + + printf("Setting population data successful.\n"); + return mio::success(); +} + +template +mio::IOResult set_mobility_weights(mio::oseirmetapopwang::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::oseirmetapopwang::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::oseirmetapopwang::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::oseirmetapopwang::Region(27), mio::AgeGroup(0), + mio::oseirmetapopwang::InfectionState::Susceptible}] -= 100; + populations[{mio::oseirmetapopwang::Region(27), mio::AgeGroup(0), + mio::oseirmetapopwang::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::oseirmetapopwang::Model model(number_regions, number_age_groups); + auto result_prepare_simulation = set_parameters_and_population(model, data_dir); + + 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_wang_nrw.h5"); +} diff --git a/cpp/examples/ode_seir_metapop.cpp b/cpp/examples/ode_seir_metapop.cpp new file mode 100644 index 0000000000..0c7519b1b0 --- /dev/null +++ b/cpp/examples/ode_seir_metapop.cpp @@ -0,0 +1,47 @@ +#include "memilio/compartments/simulation.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_metapop/infection_state.h" +#include "models/ode_metapop/model.h" +#include "models/ode_metapop/parameters.h" +#include "models/ode_metapop/regions.h" + +#include + +int main() +{ + const ScalarType t0 = 0.; + const ScalarType tmax = 10; + ScalarType dt = 0.1; + + size_t number_regions = 3; + size_t number_age_groups = 1; + + mio::oseirmetapop::Model model(number_regions, number_age_groups); + + for (size_t i = 0; i < number_regions; i++) { + model.populations[{mio::oseirmetapop::Region(i), mio::AgeGroup(0), + mio::oseirmetapop::InfectionState::Susceptible}] = 10000; + } + + model.populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed}] += + 100; + model.populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), + mio::oseirmetapop::InfectionState::Susceptible}] -= 100; + + Eigen::MatrixXd mobility_data_commuter(3, 3); + mobility_data_commuter << 0.4, 0.3, 0.3, 0.2, 0.7, 0.1, 0.4, 0.1, 0.5; + + model.set_commuting_strengths(mobility_data_commuter); + + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() + .setConstant(2.7); + + model.parameters.set>(3.335); + model.parameters.set>(8.097612257); + model.parameters.set>(0.07333); + + auto result = simulate(t0, tmax, dt, model); + return 0; +} diff --git a/cpp/memilio/compartments/simulation.h b/cpp/memilio/compartments/simulation.h index a72bee0386..f7b79c0b90 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_metapop/CMakeLists.txt b/cpp/models/ode_metapop/CMakeLists.txt new file mode 100644 index 0000000000..5b8a8c87f2 --- /dev/null +++ b/cpp/models/ode_metapop/CMakeLists.txt @@ -0,0 +1,13 @@ +add_library(ode_metapop + infection_state.h + model.h + model.cpp + parameters.h + regions.h +) +target_link_libraries(ode_metapop PUBLIC memilio) +target_include_directories(ode_metapop PUBLIC + $ + $ +) +target_compile_options(ode_metapop PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_metapop/infection_state.h b/cpp/models/ode_metapop/infection_state.h new file mode 100644 index 0000000000..03a2b70080 --- /dev/null +++ b/cpp/models/ode_metapop/infection_state.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITYIMPROVED_INFECTIONSTATE_H +#define ODESEIRMOBILITYIMPROVED_INFECTIONSTATE_H + +namespace mio +{ +namespace oseirmetapop +{ + +/** + * @brief The InfectionState enum describes the possible + * categories for the infectious state of persons + */ +enum class InfectionState +{ + Susceptible, + Exposed, + Infected, + Recovered, + Count +}; + +} // namespace oseirmetapop +} // namespace mio + +#endif // ODESEIR_INFECTIONSTATE_H diff --git a/cpp/models/ode_metapop/model.cpp b/cpp/models/ode_metapop/model.cpp new file mode 100644 index 0000000000..04ea73c6f3 --- /dev/null +++ b/cpp/models/ode_metapop/model.cpp @@ -0,0 +1,10 @@ + +#include "ode_metapop/model.h" + +namespace mio +{ +namespace oseirmetapop +{ + +} // namespace oseirmetapop +} // namespace mio diff --git a/cpp/models/ode_metapop/model.h b/cpp/models/ode_metapop/model.h new file mode 100644 index 0000000000..47906d1232 --- /dev/null +++ b/cpp/models/ode_metapop/model.h @@ -0,0 +1,248 @@ + +#ifndef ODESEIRMOBILITYIMPROVED_MODEL_H +#define ODESEIRMOBILITYIMPROVED_MODEL_H + +#include "memilio/compartments/flow_model.h" +#include "memilio/epidemiology/populations.h" +#include "models/ode_metapop/infection_state.h" +#include "models/ode_metapop/parameters.h" +#include "models/ode_metapop/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 oseirmetapop +{ + +/******************** + * 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 Eigen::MatrixXd commuting_strengths = + params.template get>().get_cont_freq_mat().get_matrix_at(t); + const size_t n_age_groups = (size_t)params.get_num_agegroups(); + const size_t n_regions = (size_t)params.get_num_regions(); + + Eigen::MatrixXd infected_pop(n_regions, n_age_groups); + for (size_t region_n = 0; region_n < n_regions; region_n++) { + for (size_t age_i = 0; age_i < n_age_groups; age_i++) { + infected_pop(region_n, age_i) = + pop[population.get_flat_index({Region(region_n), AgeGroup(age_i), InfectionState::Infected})]; + } + } + Eigen::MatrixXd infectious_share_per_region = commuting_strengths.transpose() * infected_pop; + for (size_t region_n = 0; region_n < n_regions; region_n++) { + for (size_t age_i = 0; age_i < n_age_groups; age_i++) { + infectious_share_per_region(region_n, age_i) /= + params.template get>()[{Region(region_n), AgeGroup(age_i)}]; + } + } + Eigen::MatrixXd infections_due_commuting = commuting_strengths * infectious_share_per_region; + for (size_t age_i = 0; age_i < n_age_groups; age_i++) { + for (size_t age_j = 0; age_j < n_age_groups; age_j++) { + for (size_t region_n = 0; region_n < 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 coeffStoE = + 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, age_j)) * coeffStoE * + y[population.get_flat_index({Region(region_n), AgeGroup(age_i), InfectionState::Susceptible})]; + } + } + for (size_t region = 0; region < 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>(); + Populations const& population_after_commuting = 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()) / + 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; + } + + void set_commuting_strengths(const Eigen::MatrixXd& commuting_strengths) + { + auto& commuting_strengths_param = + this->parameters.template get>().get_cont_freq_mat()[0].get_baseline(); + commuting_strengths_param = commuting_strengths; + + auto number_regions = (size_t)this->parameters.get_num_regions(); + auto number_age_groups = (size_t)this->parameters.get_num_agegroups(); + auto& population = this->populations; + auto& population_after_commuting = this->parameters.template get>(); + + 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::oseirmetapop::InfectionState::Count; state++) { + population_n += population[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age), + mio::oseirmetapop::InfectionState(state)}]; + } + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] += population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths(region_n, region_m) * population_n; + population_after_commuting[{mio::oseirmetapop::Region(region_m), mio::AgeGroup(age)}] += + commuting_strengths(region_n, region_m) * population_n; + } + } + } + } + + void set_commuting_strengths() + { + auto number_regions = (size_t)this->parameters.get_num_regions(); + set_commuting_strengths(Eigen::MatrixXd::Identity(number_regions, number_regions)); + } +}; // namespace oseirmetapop + +} // namespace oseirmetapop +} // namespace mio + +#endif // ODESEIRMOBILITY_MODEL_H diff --git a/cpp/models/ode_metapop/parameters.h b/cpp/models/ode_metapop/parameters.h new file mode 100644 index 0000000000..70bc13043b --- /dev/null +++ b/cpp/models/ode_metapop/parameters.h @@ -0,0 +1,300 @@ + +#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_metapop/regions.h" +#include "Eigen/Sparse" + +#include + +namespace mio +{ +namespace oseirmetapop +{ + +/**************************************************** + * 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 +struct PopulationAfterCommuting { + using Type = Populations; + static Type get_default(Region size_regions, AgeGroup size_agegroups) + { + return Type({size_regions, size_agegroups}, 0.); + } + static std::string name() + { + return "PopulationAfterCommuting"; + } +}; + +template +using ParametersBase = ParameterSet, TimeExposed, TimeInfected, + ContactPatterns, CommutingStrengths, PopulationAfterCommuting>; + +/** + * @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; + } + for (auto j = Region(0); j < Region(m_num_regions); j++) { + if (this->template get>()[{j, i}] <= 0.0) { + log_warning( + "Constraint check: Parameter PopulationAfterCommuting changed from {:.4f} to {:.4f}. Please " + "note that this only prevents division by zero. Consider to cancel and reset parameters.", + this->template get>()[{j, i}], 1.0); + this->template get>()[{j, i}] = 1.0; + corrected = true; + } + } + } + if ((this->template get>().get_cont_freq_mat().get_matrix_at(0).rowwise().sum() - + Eigen::VectorXd::Ones((size_t)this->get_num_regions())) + .cwiseAbs() + .maxCoeff() > 1e-10 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).minCoeff() < 0.0 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).maxCoeff() > 1.0) { + log_warning("Constraint check: Parameter CommutingStrengths does not ensure that the number of people " + "staying equals the complement of those leaving. Running without commuting."); + this->template get>().get_cont_freq_mat()[0].get_baseline() = + Eigen::MatrixXd::Identity((size_t)this->get_num_regions(), (size_t)this->get_num_regions()); + 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 + { + 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; + } + for (auto j = Region(0); j < Region(m_num_regions); j++) { + if (this->template get>()[{j, i}] <= 0.0) { + log_error("Constraint check: Parameter PopulationAfterCommuting {:.4f} smaller or equal {:.4f}", + this->template get>()[{j, i}], 0.0); + return true; + } + } + } + if ((this->template get>().get_cont_freq_mat().get_matrix_at(0).rowwise().sum() - + Eigen::VectorXd::Ones((size_t)this->get_num_regions())) + .cwiseAbs() + .maxCoeff() > 1e-10 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).minCoeff() < 0.0 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).maxCoeff() > 1.0) { + log_error("Constraint check: Parameter CommutingStrengths does not ensure that the number of people " + "staying equals the complement of those leaving."); + return true; + } + + return false; + } + +private: + Parameters(ParametersBase&& base) + : ParametersBase(std::move(base)) + , m_num_regions(base.get_num_regions()) + , m_num_agegroups(base.get_num_agegroups()) + { + } + +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 oseirmetapop +} // namespace mio + +#endif // SEIR_PARAMETERS_H diff --git a/cpp/models/ode_metapop/regions.h b/cpp/models/ode_metapop/regions.h new file mode 100644 index 0000000000..41f37337e2 --- /dev/null +++ b/cpp/models/ode_metapop/regions.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITYIMPROVED_REGIONS_H +#define ODESEIRMOBILITYIMPROVED_REGIONS_H + +#include "memilio/utils/index.h" + +namespace mio +{ +namespace oseirmetapop +{ + +/** + * @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 oseirmetapop +} // namespace mio + +#endif diff --git a/cpp/models/ode_metapop_wang/CMakeLists.txt b/cpp/models/ode_metapop_wang/CMakeLists.txt new file mode 100644 index 0000000000..fefa3afae1 --- /dev/null +++ b/cpp/models/ode_metapop_wang/CMakeLists.txt @@ -0,0 +1,13 @@ +add_library(ode_metapop_wang + infection_state.h + model.h + model.cpp + parameters.h + regions.h +) +target_link_libraries(ode_metapop_wang PUBLIC memilio) +target_include_directories(ode_metapop_wang PUBLIC + $ + $ +) +target_compile_options(ode_metapop_wang PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_metapop_wang/infection_state.h b/cpp/models/ode_metapop_wang/infection_state.h new file mode 100644 index 0000000000..9caf8c4d17 --- /dev/null +++ b/cpp/models/ode_metapop_wang/infection_state.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITY_INFECTIONSTATE_H +#define ODESEIRMOBILITY_INFECTIONSTATE_H + +namespace mio +{ +namespace oseirmetapopwang +{ + +/** + * @brief The InfectionState enum describes the possible + * categories for the infectious state of persons + */ +enum class InfectionState +{ + Susceptible, + Exposed, + Infected, + Recovered, + Count +}; + +} // namespace oseirmetapopwang +} // namespace mio + +#endif // ODESEIR_INFECTIONSTATE_H diff --git a/cpp/models/ode_metapop_wang/model.cpp b/cpp/models/ode_metapop_wang/model.cpp new file mode 100644 index 0000000000..3e93120de1 --- /dev/null +++ b/cpp/models/ode_metapop_wang/model.cpp @@ -0,0 +1,10 @@ + +#include "ode_metapop_wang/model.h" + +namespace mio +{ +namespace oseirmetapopwang +{ + +} // namespace oseirmetapopwang +} // namespace mio diff --git a/cpp/models/ode_metapop_wang/model.h b/cpp/models/ode_metapop_wang/model.h new file mode 100644 index 0000000000..a1bf45d44b --- /dev/null +++ b/cpp/models/ode_metapop_wang/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_metapop_wang/infection_state.h" +#include "models/ode_metapop_wang/parameters.h" +#include "models/ode_metapop_wang/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 oseirmetapopwang +{ + +/******************** + * 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 oseirmetapopwang +} // namespace mio + +#endif // ODESEIRMOBILITY_MODEL_H diff --git a/cpp/models/ode_metapop_wang/parameters.h b/cpp/models/ode_metapop_wang/parameters.h new file mode 100644 index 0000000000..6e1d7b45a1 --- /dev/null +++ b/cpp/models/ode_metapop_wang/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_metapop_wang/regions.h" + +#include + +namespace mio +{ +namespace oseirmetapopwang +{ + +/**************************************************** + * 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 oseirmetapopwang +} // namespace mio + +#endif // SEIR_PARAMETERS_H diff --git a/cpp/models/ode_metapop_wang/regions.h b/cpp/models/ode_metapop_wang/regions.h new file mode 100644 index 0000000000..361217436d --- /dev/null +++ b/cpp/models/ode_metapop_wang/regions.h @@ -0,0 +1,26 @@ + +#ifndef ODESEIRMOBILITY_REGIONS_H +#define ODESEIRMOBILITY_REGIONS_H + +#include "memilio/utils/index.h" + +namespace mio +{ +namespace oseirmetapopwang +{ + +/** + * @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 oseirmetapopwang +} // namespace mio + +#endif diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 74e59412fa..9a0f0ffbc3 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -84,6 +84,7 @@ set(TESTSOURCES sanitizers.cpp temp_file_register.h test_graph_abm.cpp + test_odemetapop.cpp ) if(MEMILIO_HAS_JSONCPP) diff --git a/cpp/tests/data/seir-compare.csv b/cpp/tests/data/seir-compare.csv index 239da95a90..677c4dcb36 100644 --- a/cpp/tests/data/seir-compare.csv +++ b/cpp/tests/data/seir-compare.csv @@ -1,38 +1,32 @@ # t S E I R -0.00000000000000 1049000.00000000000000 10000.00000000000000 1000.00000000000000 1000.00000000000000 -0.10020040080160 1048713.73075449210592 10092.73973116525303 1139.90293509658432 1053.62657924606287 -0.55110220440882 1046987.03773772297427 10913.21433861356854 1722.35062247182555 1377.39730119167007 -1.06524485978665 1044199.67649725102820 12546.50135443816725 2352.64073874483029 1901.18140956591151 -1.66118286922071 1039882.88403208740056 15277.40751934598848 3124.57561811346159 2715.13283045316348 -2.35595933715807 1033223.08161422202829 19620.39802993273770 4178.99925471840470 3977.52110112678702 -3.17804724413681 1022601.74967979907524 26616.01552748418908 5774.46915014757906 6007.76564256919391 -4.16892733448395 1004653.68742704729084 38433.08942357360502 8426.31628393779829 9486.90686544130585 -5.39212393400959 971794.34643305663485 59899.25324166747305 13285.70435961234216 16020.69596566357177 -6.75385596665597 915651.03697983513121 95955.08135336369742 21680.78975441652437 27713.09191238471249 -8.11558799930235 832006.67458390793763 148023.69851912744343 34434.60167287031800 46535.02522409435187 -9.20410138553930 742390.77543479565065 201243.67541724219336 48438.55588028273633 68926.99326767947059 -10.35945756851834 627065.38188297860324 264857.50646954111289 66979.08704119051981 102098.02460628982226 -11.34283729342433 518682.01091470837127 318045.68243877496570 84889.94336055630993 139382.36328596036765 -11.68513854174272 480717.31176111887908 334517.08517993631540 91305.17626443399058 154460.42679451091681 -12.02743979006111 449239.68215949274600 343612.28946372546488 97522.08154542848933 170625.94683135347441 -12.49960760584322 422867.08470962184947 338890.87290152284550 104702.56547239044448 194539.47691646512249 -13.03235808080705 398884.25140891782939 328683.93571044335840 110212.60941252954945 223219.20346810953924 -14.00378887908028 356681.26272533729207 311129.32350503187627 115039.57091181630676 278149.84285781485960 -15.03718720382851 315875.18265972472727 291988.74967131868470 115307.25494283462467 337828.81272612238536 -16.23194039567660 275050.90921587980120 268405.82756743824575 111736.34426414610061 405806.91895253618713 -17.60257755885122 236467.00561473876587 239969.97161763752229 104511.92403225359158 480051.09873537049862 -19.16760479273252 201836.59999166175839 207309.78507492708741 94020.34686743725615 557833.26806597434916 -20.88659483428036 173106.68153741111746 173218.67543084506178 81416.90213117960957 633257.74090056458954 -22.69375340601369 150852.65998401606339 140995.97202733816812 68301.84864338369516 700849.51934526243713 -24.50091197774702 134557.73041747722891 113246.00119738388457 56197.14653868280584 756999.12184645654634 -26.41330860711764 121973.81911217638117 88824.00341652805218 44973.30425799178920 805228.87321330432314 -28.32570523648826 112813.90882977400906 69081.93155133874097 35528.60144885115733 843575.55817003664561 -30.38229100048863 105671.11442069984332 52351.75372552395129 27272.90021063667882 875704.23164314008318 -32.53203962405642 100321.51816427615995 38955.94707539147203 20499.60398380133847 901222.93077653169166 -34.80205851628455 96313.32349675761361 28382.83233381239552 15052.79886064722450 921251.04530878353398 -37.21697173781426 93330.24045060748176 20191.58383072874494 10772.64150247535326 936705.53421618917491 -39.81010816461136 91131.38117136582150 13966.76619385667800 7485.05937048995293 948416.79326428822242 -42.63010344500196 89531.23135069153795 9332.43326583900853 5018.01799063735871 957118.31739283283241 -45.75944493788541 88384.30794850864913 5954.72807288235890 3209.52003529808644 963451.44394331169315 -49.04133237767380 87631.65751046490914 3712.10967585467961 2003.95478390116250 967652.27802978013642 -50.00000000000000 87471.72565843837219 3232.86999470437786 1745.83765843541028 968549.56668842269573 \ No newline at end of file +0.00000000000000 1049000.00000000000000 10000.00000000000000 1000.00000000000000 1000.00000000000000 +0.10000000000000 1048988.58104405831546 9820.84180564139024 1137.10656898469369 1053.47058131574704 +0.55000000000000 1048921.19260577531531 9071.41760821554271 1638.35319813403248 1369.03658787526888 +1.25368521533823 1048777.92434315709397 8057.65350011104692 2124.42279256048960 2039.99936417154072 +2.04101569383354 1048587.09683475596830 7102.84469418988101 2376.22157266856857 2933.83689838560213 +2.91269832810373 1048361.93679713841993 6213.99866413117434 2435.36705174795952 3988.69748698247349 +3.86108598819142 1048118.75420537753962 5400.19422411334563 2352.80445087728185 5128.24711963193658 +4.89524280895636 1047868.37748317280784 4653.03293158922588 2176.80497774428613 6301.78460749369606 +6.04258168918588 1047616.18589280080050 3957.58132139215195 1942.12075092638770 7484.11203488072078 +7.31202455451348 1047371.24218315689359 3317.02549236070581 1678.99995750832181 8632.73236697410903 +8.72217138862949 1047139.50138118292671 2731.33236444458089 1409.47894176413433 9719.68731260833556 +10.29866521143641 1046925.20438532845583 2200.92070243870512 1148.83682853681012 10725.03808369601029 +11.26330707694919 1046814.32675929483958 1929.26662121211734 1011.11757043176715 11245.28904906113530 +11.83968766077253 1046756.44987186952494 1781.52131508086131 936.25838507353467 11525.77042797587637 +12.41606824459587 1046723.19954850582872 1625.83279033276176 865.62759542813637 11785.34006573311126 +13.23077263638088 1046694.73595185985323 1416.35921567730816 770.52174862369714 12118.38308383898402 +14.29229910912589 1046662.49569856002927 1183.90972509897119 656.95470024970564 12496.63987609106516 +15.86173338430521 1046623.51623129635118 908.91131129897394 513.59281557679321 12953.97964182766191 +17.54650957365405 1046591.22194189496804 684.75600965202909 391.12584655378078 13332.89620189897141 +19.45489021353186 1046563.91490716743283 497.04591109170804 285.73393899509000 13653.30524274554955 +21.64765095553324 1046541.53458704100922 344.06106027643995 198.49179269541142 13915.91255998686938 +24.28135793516138 1046523.50970833166502 221.21786470759676 127.85456271077467 14127.41786424956445 +27.15441528731246 1046511.08352608617861 136.64984898643644 79.03686722492844 14273.22975770212179 +30.02747263946354 1046503.40388168359641 84.41354759733713 48.83710200550596 14363.34546871330349 +32.90052999161462 1046498.65905382949859 52.14589751457984 30.17173192072158 14419.02331673489971 +35.77358734376570 1046495.72778672503773 32.21290587549021 18.63911728165846 14453.42019011756929 +38.64664469591678 1046493.91696999303531 19.89941326560179 11.51439924071712 14474.66921750034635 +41.51970204806786 1046492.79833600565325 12.29280080570297 7.11301775236865 14487.79584543603960 +44.39275940021894 1046492.10730175010394 7.59384068675560 4.39405288810484 14495.90480467486668 +47.26581675237001 1046491.68041715247091 4.69107240917980 2.71441480176851 14500.91409563640809 +50.00000000000000 1046491.42674924037419 2.96615872057889 1.71632120750472 14503.89077083145639 \ No newline at end of file diff --git a/cpp/tests/test_odemetapop.cpp b/cpp/tests/test_odemetapop.cpp new file mode 100644 index 0000000000..bfedc322b5 --- /dev/null +++ b/cpp/tests/test_odemetapop.cpp @@ -0,0 +1,344 @@ + +#include "load_test_data.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "ode_metapop/model.h" +#include "ode_metapop/infection_state.h" +#include "ode_metapop/parameters.h" +#include "memilio/math/euler.h" +#include "memilio/compartments/simulation.h" +#include +#include +#include + +class ModelTestOdeMetapop : public testing::Test +{ +public: + ModelTestOdeMetapop() + : model(4, 1) + { + } + ScalarType t0; + ScalarType tmax; + ScalarType dt; + ScalarType total_population_per_region; + mio::oseirmetapop::Model model; + +protected: + void SetUp() override + { + t0 = 0.; + tmax = 50.; + dt = 0.1; + + total_population_per_region = 1061000; + + for (size_t i = 0; i < (size_t)model.parameters.get_num_regions(); i++) { + model.populations[{mio::Index( + mio::oseirmetapop::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed)}] = 10000; + model.populations[{mio::Index( + mio::oseirmetapop::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index( + mio::oseirmetapop::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index( + mio::oseirmetapop::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Susceptible)}] = + total_population_per_region - + model.populations[{ + mio::Index( + mio::oseirmetapop::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed)}] - + model.populations[{ + mio::Index( + mio::oseirmetapop::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] - + model.populations[{ + mio::Index( + mio::oseirmetapop::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}]; + } + model.set_commuting_strengths(); + } +}; + +TEST_F(ModelTestOdeMetapop, simulateDefault) +{ + mio::TimeSeries result = simulate(t0, tmax, dt, model); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); +} + +TEST_F(ModelTestOdeMetapop, checkPopulationConservation) +{ + auto result = simulate(t0, tmax, dt, model); + ScalarType num_persons = result.get_last_value().sum(); + EXPECT_NEAR(num_persons, total_population_per_region * (size_t)model.parameters.get_num_regions(), 1e-9); +} + +TEST_F(ModelTestOdeMetapop, check_constraints_parameters) +{ + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.04); + + model.parameters.get>().get_cont_freq_mat()[0].get_baseline()(0, 0) = + 10; + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + // 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>(-5.2); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.set>(5.2); + 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); + mobility_data_commuter(0, 1) += 0.5; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 1) = -0.5; + mobility_data_commuter(0, 2) = 0.75; + mobility_data_commuter(0, 3) = 0.75; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 1) = 1.5; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 0.; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 0) = 0.; + mobility_data_commuter(0, 1) = 0.; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 1.; + model.set_commuting_strengths(mobility_data_commuter); + model.parameters.set>( + mio::Populations( + {mio::oseirmetapop::Region(4), mio::AgeGroup(1)}, 0.)); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + // Nobody commutes to region 2 + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0., 0.5, 0.5, 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mio::set_log_level(mio::LogLevel::warn); +} + +TEST_F(ModelTestOdeMetapop, apply_constraints_parameters) +{ + const ScalarType tol_times = 1e-1; + model.parameters.set>(5.2); + model.parameters.set>(2); + model.parameters.set>(0.04); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + EXPECT_EQ(model.parameters.apply_constraints(), 0); + + mio::set_log_level(mio::LogLevel::off); + model.parameters.set>(-5.2); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>()[(mio::AgeGroup)0], tol_times); + + model.parameters.set>(1e-5); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>()[(mio::AgeGroup)0], tol_times); + + model.parameters.set>(10.); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR( + model.parameters.get>()[(mio::AgeGroup)0], 0.0, + 1e-14); + + model.parameters.set>(0.04); + mobility_data_commuter(0, 1) += 0.5; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>() + .get_cont_freq_mat()[0] + .get_baseline() + .isIdentity(), + true); + + mobility_data_commuter(0, 1) = -0.5; + mobility_data_commuter(0, 2) = 0.75; + mobility_data_commuter(0, 3) = 0.75; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>() + .get_cont_freq_mat()[0] + .get_baseline() + .isIdentity(), + true); + + mobility_data_commuter(0, 1) = 1.5; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 0.; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>() + .get_cont_freq_mat()[0] + .get_baseline() + .isIdentity(), + true); + + mobility_data_commuter(0, 0) = 0.; + mobility_data_commuter(0, 1) = 0.; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 1.; + model.set_commuting_strengths(mobility_data_commuter); + model.parameters.set>( + mio::Populations( + {mio::oseirmetapop::Region(4), mio::AgeGroup(1)}, 0.)); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR((model.parameters.get>()[{ + mio::oseirmetapop::Region(3), mio::AgeGroup(0)}]), + 1.0, tol_times); + + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0., 0.5, 0.5, 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR((model.parameters.get>()[{ + mio::oseirmetapop::Region(1), mio::AgeGroup(0)}]), + 1.0, tol_times); + + mio::set_log_level(mio::LogLevel::warn); +} + +TEST(TestOdeMetapop, compareSEIR) +{ + ScalarType t0 = 0; + ScalarType tmax = 50.; + ScalarType dt = 0.1; + + ScalarType total_population = 1061000; + + mio::oseirmetapop::Model model(1, 1); + + model.populations[{mio::Index( + mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed)}] = 10000; + model.populations[{mio::Index( + mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index( + mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index( + mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Susceptible)}] = + total_population - + model.populations[{mio::Index( + mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed)}] - + model.populations[{mio::Index( + mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] - + model.populations[{mio::Index( + mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}]; + + // The model with a single region should correspond to the SEIR model. + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.1); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + model.set_commuting_strengths(); + + model.check_constraints(); + + // Use the Euler integrator as adaptive methods make different time steps in this model due to restructured equations. + std::shared_ptr> integrator = std::make_shared>(); + + auto result = simulate(t0, tmax, dt, model, integrator); + std::vector> refData = load_test_data_csv("seir-compare-euler.csv"); + + ASSERT_EQ(refData.size(), static_cast(result.get_num_time_points())); + + for (Eigen::Index irow = 0; irow < result.get_num_time_points(); ++irow) { + ScalarType t = refData[static_cast(irow)][0]; + auto rel_tol = 1e-10; + + //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 < refData[static_cast(irow)].size() - 1; ++icol) { + ScalarType ref = refData[static_cast(irow)][icol + 1]; + ScalarType actual = result[irow][icol]; + + ScalarType tol = rel_tol * ref; + ASSERT_NEAR(ref, actual, tol) << "at row " << irow; + } + } +} + +TEST_F(ModelTestOdeMetapop, compareWithPreviousRun) +{ + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.1); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + std::vector> refData = load_test_data_csv("ode-seir-metapop-compare.csv"); + std::shared_ptr> 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; + } + } +} diff --git a/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py b/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py new file mode 100644 index 0000000000..5232f2a3b5 --- /dev/null +++ b/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py @@ -0,0 +1,56 @@ +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 = os.path.join( + arg_dict['out_folder'].split('/pydata/')[0], 'Germany/') + directory_mobility = os.path.join(directory, 'mobility/') + directory_population = os.path.join(directory, 'pydata/') + mobility_file = 'commuter_mobility_2022' + 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_equationbased.sh b/shellscripts/likwid_equationbased.sh new file mode 100644 index 0000000000..e4980ff289 --- /dev/null +++ b/shellscripts/likwid_equationbased.sh @@ -0,0 +1,28 @@ +#!/bin/bash +#SBATCH -N 1 +#SBATCH -n 1 +#SBATCH -c 1 +#SBATCH --exclusive +#SBATCH -t 5-0:00:00 +#SBATCH --output=likwid_equationbased-%A.out +#SBATCH --error=likwid_equationbased-%A.err +#SBATCH --exclude="be-cpu05, be-gpu01" +#SBATCH --job-name=likwid_mobilitymodels + +warm_up_runs=0 +runs=50 +regions=400 + +module purge +module load PrgEnv/gcc12-openmpi + +echo Running $1 on node $SLURM_JOB_NODELIST with $warm_up_runs warm up runs and $runs runs. +cd ../cpp/build +rm -rf CMakeCache.txt CMakeFiles/ +cmake -DCMAKE_BUILD_TYPE="Release" -DMEMILIO_ENABLE_OPENMP=ON -DMEMILIO_ENABLE_WARNINGS_AS_ERRORS=OFF .. +cmake --build . --target ode_metapop_timing + +perf record -C 0 likwid-pin ./bin/ode_metapop_timing $warm_up_runs $runs $regions +perf report +# srun --cpu-bind=cores --cpus-per-task=1 --cpu-freq=2200000-2200000 valgrind --tool=massif --detailed-freq=2 ./bin/ode_metapop_timing $warm_up_runs $runs $regions +# srun --cpu-bind=cores --cpus-per-task=1 --cpu-freq=2200000-2200000 likwid-perfctr -C 0 -g MEM_DP -m ./bin/ode_metapop_timing $warm_up_runs $runs $regions diff --git a/shellscripts/likwid_graphbased.sh b/shellscripts/likwid_graphbased.sh new file mode 100644 index 0000000000..4963e4e860 --- /dev/null +++ b/shellscripts/likwid_graphbased.sh @@ -0,0 +1,27 @@ +#!/bin/bash +#SBATCH -N 1 +#SBATCH -n 1 +#SBATCH -c 1 +#SBATCH --exclusive +#SBATCH -t 5-0:00:00 +#SBATCH --output=likwid_graphbased-%A.out +#SBATCH --error=likwid_graphbased-%A.err +#SBATCH --exclude="be-cpu05, be-gpu01" +#SBATCH --job-name=likwid_mobilitymodels + +warm_up_runs=0 +runs=1 +regions=80 + +module purge +module load PrgEnv/gcc13-openmpi + +echo Running $1 on node $SLURM_JOB_NODELIST with $warm_up_runs warm up runs and $runs runs. +cd ../cpp/build +rm -rf CMakeCache.txt CMakeFiles/ +CXX="g++ -fverbose-asm" cmake -DCMAKE_BUILD_TYPE="Release" -DMEMILIO_ENABLE_OPENMP=ON -DMEMILIO_ENABLE_WARNINGS_AS_ERRORS=OFF .. +# cmake --build . --target graph_timing + +# perf record ./bin/graph_timing $warm_up_runs $runs $regions +# perf report +srun --cpu-bind=cores --cpus-per-task=1 --cpu-freq=2200000-2200000 likwid-perfctr -C 0 -g MEM_DP -m ./bin/graph_timing $warm_up_runs $runs $regions 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/shellscripts/steps_mobilitymodels.sh b/shellscripts/steps_mobilitymodels.sh new file mode 100644 index 0000000000..c2c6d80d1f --- /dev/null +++ b/shellscripts/steps_mobilitymodels.sh @@ -0,0 +1,20 @@ +#!/bin/sh +#SBATCH -N 1 +#SBATCH -n 1 +#SBATCH -c 1 +#SBATCH --exclusive +#SBATCH -t 5-0:00:00 +#SBATCH --output=steps_equationbased-%A.out +#SBATCH --error=steps_equationbased-%A.err +#SBATCH --exclude="be-cpu05, be-gpu01" +#SBATCH --job-name=steps_mobilitymodels + +echo Running $1 on node $SLURM_JOB_NODELIST. +cd ../cpp/build +cmake -DCMAKE_BUILD_TYPE="Release" -DMEMILIO_ENABLE_OPENMP=ON .. +cmake --build . --target $1 + +for i in $(seq 2 12) +do + srun --cpu-bind=core --cpus-per-task=1 ./bin/$1 1e-$i +done \ No newline at end of file diff --git a/shellscripts/timing_mobilitymodels.sh b/shellscripts/timing_mobilitymodels.sh new file mode 100644 index 0000000000..a0d28bb84e --- /dev/null +++ b/shellscripts/timing_mobilitymodels.sh @@ -0,0 +1,27 @@ +#!/bin/bash +#SBATCH -N 1 +#SBATCH -n 1 +#SBATCH -c 1 +#SBATCH --exclusive +#SBATCH -t 5-0:00:00 +#SBATCH --output=timing_equationbased_noage_euler-%A.out +#SBATCH --error=timing_equationbased_noage_euler-%A.err +#SBATCH --exclude="be-cpu05, be-gpu01" +#SBATCH --job-name=timing_mobilitymodels + +warm_up_runs=10 +runs=100 + +echo Running $1 on node $SLURM_JOB_NODELIST with $warm_up_runs warm up runs and $runs runs. +cd ../cpp/build +rm -rf CMakeCache.txt CMakeFiles/ +cmake -DCMAKE_BUILD_TYPE="Release" -DMEMILIO_ENABLE_OPENMP=ON .. +cmake --build . --target $1 +for i in $(seq 1 1 4) +do + srun --cpu-bind=core --cpus-per-task=1 ./bin/$1 $warm_up_runs $runs $i +done +for i in $(seq 300 5 350) +do + srun --cpu-bind=core --cpus-per-task=1 ./bin/$1 $warm_up_runs $runs $i +done 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..dce00b1ff0 --- /dev/null +++ b/tools/plot_mobility_runtimes.py @@ -0,0 +1,236 @@ +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 = ['Model C (ODE)', 'Model D (Graph-ODE)'] + + +def plot_flops(name='number_flops'): + fig, ax = plt.subplots() + + def flops_equation_based(x, eta): + return (4*x**2+22*x)/eta + + def flops_graph_based(x, eta): + return (43*x**2+23*x/eta) + + 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, linestyle=linestyles[idx], + color=colors[0], + label=models[0] + ', $h=$' + str(eta)) + ax.plot( + x, flops_graph_based(x, eta), + linewidth=1.5, linestyle=linestyles[idx], + color=colors[1], + label=models[1] + ', $h=$' + 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 FLOP', 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() + plt.grid(True, linestyle='--') + + for idx, file in enumerate(files): + df = pd.read_json(file) + + ax1.plot(df["Regions"], df["Time"], + linestyle='--', marker='o', linewidth=1.2, label=models[idx]) + + 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)*200 + + def flops_graph_based(x): + return (43*x**2+230*x)*20 + + x = np.linspace(0, 400, 400) + + ax2.plot(x, flops_equation_based(x), linewidth=2, + color='darkblue', label='FLOP Model C') + ax2.plot(x, flops_graph_based(x), linewidth=2, + color='#9C180D', label='FLOP Model D') + ax2.set_ylabel('Number of FLOP', fontsize=fontsize_labels) + ax2.set_ylim(bottom=0.) + + handles1, labels1 = ax1.get_legend_handles_labels() + handles2, labels2 = ax2.get_legend_handles_labels() + + handles = handles1 + handles2 + labels = labels1 + labels2 + + plt.tight_layout() + fig.legend(handles, labels, loc='upper left', + bbox_to_anchor=(0.125, 0.925), ncols=2) + + 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.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.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(fontsize=fontsize_legends) + 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() + + +def plot_steps(files, name='', models=[]): + for idx, file in enumerate(files): + df = pd.read_json(file) + df.set_index('Absolute tolerance', inplace=True) + model_type = os.path.basename(file).split('_')[0] + if model_type == 'ode': + plt.plot(df, label=models[idx], color=colors[idx]) + else: + plt.plot(df['Steps Hotspot'], color=colors[idx], label=models[idx]) + plt.plot(df['Steps other Regions'], color=colors[idx]) + plt.fill_between( + df.index, df['Steps Hotspot'], + df['Steps other Regions'], + color=colors[idx], + alpha=0.15) + + plt.ylim(bottom=10.) + plt.xlim(left=df.index.min()/1.2, right=df.index.max()*1.2) + plt.yticks(fontsize=fontsize_legends) + plt.xticks(df.index, fontsize=fontsize_legends) + + plt.xscale('log') + plt.gca().invert_xaxis() + plt.ylabel('Number of steps', fontsize=fontsize_labels) + plt.xlabel('Absolute tolerance', fontsize=fontsize_labels) + plt.grid(True, linestyle='--') + plt.legend(fontsize=fontsize_legends) + + plot_dir = os.path.join(os.path.dirname(__file__), '../Plots') + plt.savefig( + os.path.join(plot_dir, 'compare_steps.png'), + bbox_inches='tight', dpi=500) + plt.close() + + +def plot_quotient(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.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') + plt.plot( + merged_df.index, merged_df[models[1]] / merged_df[models[0]], + label='Quotient', 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(fontsize=fontsize_legends) + 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') + + ode_timing_euler = os.path.join(result_dir, 'ode_timing_euler.json') + ode_timing_noage_euler = os.path.join( + result_dir, 'ode_timing_noage_euler.json') + graphbased_timing_euler = os.path.join( + result_dir, 'graphbased_timing_euler.json') + graphbased_timing_noage_euler = os.path.join( + result_dir, 'graphbased_timing_noage_euler.json') + + ode_steps = os.path.join(result_dir, 'ode_steps.json') + graphbased_steps = os.path.join(result_dir, 'graphbased_steps.json') + + timings_euler = [ode_timing_euler, graphbased_timing_euler] + timings_euler_noage = [ode_timing_noage_euler, + graphbased_timing_noage_euler] + + compare_runtimes( + timings_euler, name='compare_runtimes_euler', models=models) + compare_runtimes( + timings_euler_noage, name='compare_runtimes_euler_noage', + models=models) + plot_quotient(timings_euler, name='quotient', + title='Quotient', models=models) + plot_quotient(timings_euler_noage, name='quotient_noage', + title='Quotient', models=models) + compare_runtime_and_flops( + timings_euler_noage, 'compare_runtimes_and_flops') + plot_flops('number_flops') + plot_steps([ode_steps, graphbased_steps], models=models) diff --git a/tools/plot_results_mobility.py b/tools/plot_results_mobility.py new file mode 100644 index 0000000000..3e8cc74d7c --- /dev/null +++ b/tools/plot_results_mobility.py @@ -0,0 +1,475 @@ +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', + fontsize=10): + """! 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, edgecolor='black', + linewidth=0.1) + + 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, edgecolor='black', linewidth=0.1) + 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, edgecolor='black', linewidth=0.1) + + ax.set_title(legend[i], fontsize=fontsize) + 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=''): + + dfs = {} + min_vals = [] + max_vals = [] + + for date in range(10, 101, 20): + dfs[date] = extract_nrw_data_and_combine( + files=files, date=date, relative=True) + min_vals.append(dfs[date].drop(columns='Region').min().min()) + max_vals.append(dfs[date].drop(columns='Region').max().max()) + + min_val = min(min_vals) + max_val = max(max_vals) + + 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=12) + # cbar.set_ticks([min_val, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, max_val]) + # cbar.ax.xaxis.set_major_formatter(FormatStrFormatter('%.5f')) + plt.tight_layout() + 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): + plot_map_nrw( + dfs[date], 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, + fontsize=13) + + +def plot_difference(files, output_dir, name='difference2D'): + 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'] = ( + dfs_all[dfs_all.columns[1]] - dfs_all[dfs_all.columns[2]]).abs().sum() / total_population + + df_dif['difference'].plot(label='Relative difference') + df_dif['absolute value'].plot( + label='Relative difference in absolute value') + plt.xlim(left=0., right=101.) + plt.tight_layout() + plt.legend() + plt.grid(linestyle='--') + plt.savefig(os.path.join(output_dir, name)) + plt.close() + + +def plot_difference_maps(files, output_dir): + + df_dif1 = pd.DataFrame(columns=['Region']) + df_dif2 = pd.DataFrame(columns=['Region']) + + for date in range(10, 51, 10): + dfs_all = extract_nrw_data_and_combine( + files=files, date=date, relative=True) + df_dif1['Region'] = dfs_all['Region'] + df_dif1['Count (rel)' + str(date)] = dfs_all[dfs_all.columns[1] + ] - dfs_all[dfs_all.columns[2]] + + for date in range(60, 101, 10): + dfs_all = extract_nrw_data_and_combine( + files=files, date=date, relative=True) + df_dif2['Region'] = dfs_all['Region'] + + df_dif2['Count (rel)' + str(date)] = dfs_all[dfs_all.columns[1] + ] - dfs_all[dfs_all.columns[2]] + + min_val1 = df_dif1.drop(columns=['Region']).min().min() + max_val1 = df_dif1.drop(columns=['Region']).max().max() + min_val2 = df_dif2.drop(columns=['Region']).min().min() + max_val2 = df_dif2.drop(columns=['Region']).max().max() + min_val = min(min_val1, min_val2) + max_val = max(max_val1, max_val2) + maximum_abs = abs(max([min_val, max_val], key=abs)) + + cmap = 'seismic' + norm = mcolors.TwoSlopeNorm(vmin=-maximum_abs, vmax=maximum_abs, vcenter=0) + 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) + ticks = np.linspace(-maximum_abs, maximum_abs, 10) + cbar.ax.tick_params(labelsize=13) + cbar.set_ticks(ticks) + cbar.ax.xaxis.set_major_formatter(FormatStrFormatter('%.3f')) + plt.tight_layout() + plt.subplots_adjust(bottom=0.3) + plt.savefig(os.path.join(output_dir, 'colorbar_difference.png'), dpi=900) + plt.close() + + plot_map_nrw( + df_dif1, scale_colors=[-maximum_abs, maximum_abs], + legend=['Day ' + str(date) for date in range(10, 51, 10)], + title='', plot_colorbar=False, + output_path=output_dir, + fig_name="difference10-50", dpi=900, + outercolor='white', + log_scale=False, + cmap='seismic', + fontsize=17) + + plot_map_nrw( + df_dif2, scale_colors=[-maximum_abs, maximum_abs], + legend=['Day ' + str(date) for date in range(60, 101, 10)], + title='', plot_colorbar=False, + output_path=output_dir, + fig_name="difference60-100", 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': + 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(r'(\d+)') + df['Group'] = df['Group'].apply(pd.to_numeric, errors='coerce') + 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 = ['#2ca02c', '#ff7f0e', '#9C180D'] + 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_title(compartment, fontsize=8) + # ax.set_ylim(bottom=0.) + ax.set_xlim(left=0., right=dates.max()+1) + else: + df = pd.DataFrame() + regions = list(h5file.keys()) + population = 0 + for i in range(len(regions)): + for comp in compartments.keys(): + population += h5file[regions[i] + ]['Total'][:, compartments[comp]] + 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.5, color=colors[file_idx], + linestyle='--') + ax.set_title(compartment, fontsize=8) + ax.set_ylim(bottom=0.) + ax.set_xlim(left=0., right=df['Time'].max()+1) + + file_idx = file_idx+1 + ax.tick_params(labelsize=7, ) + 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 = {'Model B': 'results/ode_result_wang_nrw', + 'Model C': 'results/ode_result_nrw', + 'Model D': 'results/graph_result_nrw'} + + file_format = 'h5' + + models = ['Model B (ODE)', + 'Model C (ODE)', + 'Model D (Graph-ODE)'] + + plot_dir = os.path.join(os.path.dirname(__file__), '../Plots') + + # plot_maps(files=results, output_dir=plot_dir, + # legend=models, name='NRWAdaptiveDay') + # plot_difference_maps( + # files={key: value for key, value in results.items() + # if key in {'Model C', 'Model D'}}, + # output_dir=plot_dir) + # plot_difference( + # files={key: value for key, value in results.items() + # if key in {'Model C', 'Model D'}}, + # output_dir=plot_dir) + # compare_compartments(files=results, output_dir=plot_dir, legend=models) + plot_difference(files={'Old result': 'cpp/build/ode_result_nrw_test', 'New result': 'cpp/build/ode_result_nrw'}, output_dir=plot_dir, name='difference_ode_results_test') + # plot_difference(files={'Old result': 'results/graph_result_nrw', 'New result': 'cpp/build/graph_result_nrw'}, output_dir=plot_dir, name='difference_graph_results') \ No newline at end of file