diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index f4b87f1a9a..ce6a28f40a 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -136,7 +136,7 @@ endif() # exclude some warnings we accept if(CMAKE_CXX_COMPILER_ID MATCHES "GNU") string(APPEND MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS - "-Wno-unknown-warning;-Wno-pragmas;-Wno-deprecated-copy;-Wno-expansion-to-defined;") + "-Wno-unknown-warning;-Wno-pragmas;-Wno-deprecated-copy;-Wno-expansion-to-defined;-Wno-shadow;") elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang") string(APPEND MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS "-Wno-unknown-warning-option;-Wno-deprecated;-Wno-gnu-zero-variadic-macro-arguments;") diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index c78b475ca3..aa28694c7d 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -112,6 +112,22 @@ 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}) +add_executable(asymmetric_graph_example asymmetric_graph.cpp) +target_link_libraries(asymmetric_graph_example PRIVATE memilio smm) +target_compile_options(asymmetric_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(asymmetric_graph_data asym_ex_data.cpp) +target_link_libraries(asymmetric_graph_data PRIVATE memilio smm) +target_compile_options(asymmetric_graph_data PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(asym_ben asym_graph_ben.cpp) +target_link_libraries(asym_ben PRIVATE memilio smm) +target_compile_options(asym_ben PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + +add_executable(asymmetric_params asymmetric_params.cpp) +target_link_libraries(asymmetric_params PRIVATE memilio smm) +target_compile_options(asymmetric_params PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(abm_minimal_example abm_minimal.cpp) target_link_libraries(abm_minimal_example PRIVATE memilio abm) target_compile_options(abm_minimal_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/abm_parameter_study.cpp b/cpp/examples/abm_parameter_study.cpp index 8446465c5f..cc13417410 100644 --- a/cpp/examples/abm_parameter_study.cpp +++ b/cpp/examples/abm_parameter_study.cpp @@ -228,4 +228,4 @@ int main() mio::mpi::finalize(); return 0; -} +} \ No newline at end of file diff --git a/cpp/examples/asym_ex_data.cpp b/cpp/examples/asym_ex_data.cpp new file mode 100644 index 0000000000..1149f2dfc5 --- /dev/null +++ b/cpp/examples/asym_ex_data.cpp @@ -0,0 +1,121 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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/config.h" +#include "memilio/geography/geolocation.h" +#include "memilio/geography/rtree.h" +#include "memilio/geography/distance.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_asymmetric.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_builder.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/index.h" +#include "memilio/utils/logging.h" +#include "memilio/timer/auto_timer.h" +#include "memilio/utils/parameter_distributions.h" +#include "memilio/utils/random_number_generator.h" +#include "smm/simulation.h" +#include "smm/parameters.h" +#include "fmd/infection_state.h" +#include "fmd/model.h" +#include "fmd/adoption_rates.h" +#include "thirdparty/csv.h" +#include +#include + +int main(int /*argc*/, char** /*argv*/) +{ + const auto t0 = 0.; + const auto tmax = 1000.; + const auto dt = 1.; //initial time step + using mio::fmd::InfectionState; + using Status = mio::Index; + using mio::regions::Region; + + //total compartment sizes + + using Model = mio::smm::Model; + + mio::fmd::Builder builder; + mio::log_info("Starting Graph generation"); + { + mio::timing::AutoTimer<"Graph Nodes Generation"> timer; + io::CSVReader<4> farms("/home/kilian/Documents/data/read_data/farms.csv"); + farms.read_header(io::ignore_extra_column, "id_dec", "x", "y", "farm_size"); + size_t farm_id, num_cows; + double latitude, longitude; + while (farms.read_row(farm_id, latitude, longitude, num_cows)) { + auto curr_model = mio::fmd::create_model(num_cows, mio::fmd::generic_adoption_rates()); + builder.add_node(farm_id, longitude, latitude, curr_model, t0); + } + } + mio::log_info("Nodes added to Graph"); + auto rng = mio::RandomNumberGenerator(); + + std::vector> interesting_indices; + interesting_indices.push_back( + {Model(Status{InfectionState::Count}, Region(1)).populations.get_flat_index({Region(0), InfectionState::I})}); + + { + mio::timing::AutoTimer<"Graph Edges Generation"> timer; + io::CSVReader<2> edges("/home/kilian/Documents/data/read_data/edges.csv"); + edges.read_header(io::ignore_extra_column, "from", "to"); + size_t from, to; + while (edges.read_row(from, to)) { + builder.add_edge(from, to, interesting_indices); + } + } + auto graph = std::move(builder).build(); + mio::log_info("Graph generated"); + + mio::fmd::generate_neighbours(graph, {mio::geo::kilometers(2.0)}); + + mio::log_info("Neighbors set"); + auto sim = mio::make_mobility_sim(t0, dt, std::move(graph)); + + io::CSVReader<4> exchanges("/home/kilian/Documents/data/read_data/exchanges.csv"); + exchanges.read_header(io::ignore_extra_column, "from_dec", "to_dec", "day", "weight"); + + size_t date, num_animals; + size_t from, to; + while (exchanges.read_row(from, to, date, num_animals)) { + sim.add_exchange(date, num_animals, from, to); + } + mio::log_info("Exchanges added"); + + // #ifdef MEMILIO_ENABLE_OPENMP + // #pragma omp parallel for + // for (size_t i = 0; i < 100; ++i) { + // #endif + + auto sim2(sim); + mio::log_info("new Simulation created"); + + sim2.advance(tmax); + mio::log_info("Simulation finished"); + + auto sth = sim2.exchanges_per_timestep().export_csv("Exchange_statistics.csv", {"Commuter Sick", "Commuter Total"}); + + sth = sim2.statistics_per_timestep().export_csv("Simulation_statistics.csv"); + + mio::log_info("Finished postprocessing"); + + return 0; +} diff --git a/cpp/examples/asym_graph_ben.cpp b/cpp/examples/asym_graph_ben.cpp new file mode 100644 index 0000000000..005fef597f --- /dev/null +++ b/cpp/examples/asym_graph_ben.cpp @@ -0,0 +1,168 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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/config.h" +#include "memilio/geography/geolocation.h" +#include "memilio/geography/rtree.h" +#include "memilio/geography/distance.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_asymmetric.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_builder.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/index.h" +#include "memilio/utils/logging.h" +#include "memilio/timer/auto_timer.h" +#include "memilio/utils/parameter_distributions.h" +#include "memilio/utils/random_number_generator.h" +#include "smm/simulation.h" +#include "smm/parameters.h" +#include "fmd/infection_state.h" +#include "fmd/model.h" +#include "fmd/adoption_rates.h" +#include "thirdparty/csv.h" +#include +#include +#include + +int main(int /*argc*/, char** /*argv*/) +{ + const auto t0 = 0.; + // const auto tmax = 100.; + const auto dt = 1.; //initial time step + + //total compartment sizes + using mio::fmd::InfectionState; + using Status = mio::Index; + using mio::regions::Region; + + using Model = mio::smm::Model; + auto home = Region(0); + + auto adoption_rates = mio::fmd::generic_adoption_rates(); + + mio::log_info("Reading CSVs"); + std::vector latitudes, longitudes; + std::vector farm_ids, num_cows_vec, dates, num_animals_exchanges; + std::vector froms, tos, from_exchanges, to_exchanges; + io::CSVReader<4> farms("../../../farms16mio.csv"); + farms.read_header(io::ignore_extra_column, "farms", "num_cows", "latitude", "longitude"); + int farm_id, num_cows; + double latitude, longitude; + while (farms.read_row(farm_id, num_cows, latitude, longitude)) { + farm_ids.push_back(farm_id); + num_cows_vec.push_back(num_cows); + latitudes.push_back(latitude); + longitudes.push_back(longitude); + } + io::CSVReader<2> edges("../../../edges16mio.csv"); + edges.read_header(io::ignore_extra_column, "from", "to"); + size_t from, to; + while (edges.read_row(from, to)) { + froms.push_back(from); + tos.push_back(to); + } + io::CSVReader<4> exchanges("../../../trade16mio.csv"); + exchanges.read_header(io::ignore_extra_column, "date", "num_animals", "from", "to"); + int date, num_animals; + while (exchanges.read_row(date, num_animals, from, to)) { + dates.push_back(date); + num_animals_exchanges.push_back(num_animals); + from_exchanges.push_back(from); + to_exchanges.push_back(to); + } + std::vector> interesting_indices; + interesting_indices.push_back( + {Model(Status{InfectionState::Count}, Region(1)).populations.get_flat_index({home, InfectionState::I})}); + + for (size_t num_edges_benchmark = 1024; num_edges_benchmark < 16000000; num_edges_benchmark *= 2) { + mio::log_info("Running with {} edges", num_edges_benchmark); + + mio::fmd::Graph graph; + mio::fmd::Builder builder; + + auto start = std::chrono::high_resolution_clock::now(); + + for (size_t i = 0; i < farm_ids.size(); ++i) { + Model curr_model(Status{InfectionState::Count}, Region(1)); + curr_model.populations[{home, InfectionState::S}] = num_cows_vec[i]; + curr_model.populations[{home, InfectionState::E}] = 0; + curr_model.populations[{home, InfectionState::I}] = 0; + curr_model.populations[{home, InfectionState::INS}] = 0; + curr_model.populations[{home, InfectionState::ICS}] = 0; + curr_model.populations[{home, InfectionState::R}] = 0; + curr_model.populations[{home, InfectionState::V}] = 0; + curr_model.populations[{home, InfectionState::D}] = 0; + curr_model.parameters.get>() = adoption_rates; + graph.add_node(farm_ids[i], longitudes[i], latitudes[i], curr_model, t0); + } + auto rng = mio::RandomNumberGenerator(); + + for (size_t i = 0; i < num_edges_benchmark; ++i) { + graph.add_edge(froms[i], tos[i], interesting_indices); + } + auto end = std::chrono::high_resolution_clock::now(); + auto duration = std::chrono::duration_cast(end - start); + + mio::log_info("Classic graph generation finished for {} edges in {} milliseconds", num_edges_benchmark, + duration.count()); + start = std::chrono::high_resolution_clock::now(); + for (size_t i = 0; i < farm_ids.size(); ++i) { + Model curr_model(Status{InfectionState::Count}, Region(1)); + curr_model.populations[{home, InfectionState::S}] = num_cows_vec[i]; + curr_model.populations[{home, InfectionState::E}] = 0; + curr_model.populations[{home, InfectionState::I}] = 0; + curr_model.populations[{home, InfectionState::INS}] = 0; + curr_model.populations[{home, InfectionState::ICS}] = 0; + curr_model.populations[{home, InfectionState::R}] = 0; + curr_model.populations[{home, InfectionState::V}] = 0; + curr_model.populations[{home, InfectionState::D}] = 0; + curr_model.parameters.get>() = adoption_rates; + builder.add_node(farm_ids[i], longitudes[i], latitudes[i], curr_model, t0); + } + for (size_t i = 0; i < num_edges_benchmark; ++i) { + builder.add_edge(from_exchanges[i], to_exchanges[i], interesting_indices); + } + auto graph2 = std::move(builder).build(true); + end = std::chrono::high_resolution_clock::now(); + duration = std::chrono::duration_cast(end - start); + + mio::log_info("New graph construction finished using {} edges in {} milliseconds", num_edges_benchmark, + duration.count()); + + // auto nodes = graph.nodes() | std::views::transform([](const auto& node) { + // return &node.property; + // }); + // auto tree = mio::geo::RTree(nodes.begin(), nodes.end()); + + // for (auto& node : graph.nodes()) { + // node.property.set_regional_neighbors( + // tree.in_range_indices_query(node.property.get_location(), {mio::geo::kilometers(2.0)})); + // } + + auto sim = mio::make_mobility_sim(t0, dt, std::move(graph)); + + for (size_t i = 0; i < dates.size(); ++i) { + sim.add_exchange(dates[i], num_animals_exchanges[i], from_exchanges[i], to_exchanges[i]); + } + + mio::log_info("Starting and immediately finishing simulation"); + } + return 0; +} diff --git a/cpp/examples/asymmetric_graph.cpp b/cpp/examples/asymmetric_graph.cpp new file mode 100644 index 0000000000..3094da5aa0 --- /dev/null +++ b/cpp/examples/asymmetric_graph.cpp @@ -0,0 +1,161 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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/config.h" +#include "memilio/geography/geolocation.h" +#include "memilio/geography/rtree.h" +#include "memilio/geography/distance.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_asymmetric.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_builder.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/index.h" +#include "memilio/utils/logging.h" +#include "memilio/timer/auto_timer.h" +#include "memilio/utils/parameter_distributions.h" +#include "memilio/utils/random_number_generator.h" +#include "smm/simulation.h" +#include "smm/parameters.h" +#include "fmd/infection_state.h" +#include "fmd/model.h" +#include "fmd/adoption_rates.h" +#include "thirdparty/csv.h" +#include +#include + +int main(int /*argc*/, char** /*argv*/) +{ + const auto t0 = 0.; + const auto tmax = 100.; + const auto dt = 1.; //initial time step + + using mio::fmd::InfectionState; + using Status = mio::Index; + using mio::regions::Region; + + //total compartment sizes + + using Model = mio::smm::Model; + auto home = Region(0); + + auto adoption_rates = mio::fmd::generic_adoption_rates(); + + mio::fmd::Builder builder; + mio::log_info("Starting Graph generation"); + { + mio::timing::AutoTimer<"Graph Nodes Generation"> timer; + io::CSVReader<4> farms("../../farms200000.csv"); + farms.read_header(io::ignore_extra_column, "farms", "num_cows", "latitude", "longitude"); + int farm_id, num_cows; + double latitude, longitude; + while (farms.read_row(farm_id, num_cows, latitude, longitude)) { + builder.add_node(farm_id, longitude, latitude, mio::fmd::create_model(num_cows, adoption_rates), t0); + } + } + mio::log_info("Nodes added to Graph"); + auto rng = mio::RandomNumberGenerator(); + + std::vector> interesting_indices; + interesting_indices.push_back( + {Model(Status{InfectionState::Count}, Region(1)).populations.get_flat_index({home, InfectionState::I})}); + io::CSVReader<2> edges("../../edges200000.csv"); + edges.read_header(io::ignore_extra_column, "from", "to"); + size_t from, to; + while (edges.read_row(from, to)) { + builder.add_edge(from, to, interesting_indices); + } + auto graph = std::move(builder).build(); + // // mio::log_info("Graph generated"); + + auto nodes = graph.nodes() | std::views::transform([](const auto& node) { + return &node.property; + }); + auto tree = mio::geo::RTree(nodes.begin(), nodes.end()); + // mio::log_info("RTree generated"); + + for (auto& node : graph.nodes()) { + mio::timing::AutoTimer<"neighbourhood search"> timer; + node.property.set_regional_neighbors( + tree.in_range_indices_query(node.property.get_location(), {mio::geo::kilometers(5.0)})); + } + + mio::log_info("Neighbors set"); + auto sim = mio::make_mobility_sim(t0, dt, std::move(graph)); + + io::CSVReader<4> exchanges("../../trade200000.csv"); + exchanges.read_header(io::ignore_extra_column, "date", "num_animals", "from", "to"); + + int date, num_animals; + while (exchanges.read_row(date, num_animals, from, to)) { + sim.add_exchange(date, num_animals, from, to); + } + // mio::log_info("Exchanges added"); + + // // #ifdef MEMILIO_ENABLE_OPENMP + // // #pragma omp parallel for + // // for (size_t i = 0; i < 100; ++i) { + // // #endif + + // // mio::log_info("new Simulation created"); + // sim2.infectionrisk = 0.1; + + sim.advance(tmax); + + mio::log_info("Simulation finished"); + + // #ifdef MEMILIO_ENABLE_OPENMP + // } + // #endif + + // sim2.get_graph().nodes()[28].property.get_result().print_table({"S", "E", "I", "R", "D"}); + // std::cout << "Second Table" << std::endl; + // sim2.get_graph().nodes()[1].property.get_result().print_table({"S", "E", "I", "R", "D"}); + + // auto& edge_1_0 = sim2.get_graph().edges()[1]; + // auto& results = edge_1_0.property.get_mobility_results(); + // results.print_table({"Commuter Sick", "Commuter Total"}); + + // auto exchange_results = sim2.sum_exchanges(); + // // mio::log_info("Sum of exchanged sick animals: {}", exchange_results[0]); + // // mio::log_info("Sum of exchanged animals: {}", exchange_results[1]); + + // auto sth = sim2.exchanges_per_timestep().export_csv("Exchange_statistics.csv", {"Commuter Sick", "Commuter Total"}); + + // // for (auto node : sim2.get_graph().nodes()) { + // // if (node.property.get_result().get_num_time_points() < num_time_points) { + // // mio::log_error("Node with inconsistent number of time points in results."); + // // } + // // } + + // // exchange_results = sim2.sum_nodes(); + // // // mio::log_info("{}, {}, {}, {}", exchange_results[0], exchange_results[1], exchange_results[2], exchange_results[3]); + + // sth = sim2.statistics_per_timestep().export_csv("Simulation_statistics.csv"); + // // // auto combined_results = sim2.combine_node_results(); + // // // combined_results.print_table({"S", "E", "I", "R", "D"}); + // // // auto ioresult = combined_results.export_csv("Simulation_results.csv"); + + // sim2.statistics_per_timestep({0, 1, 2, 3, 4}).print_table({"S", "E", "I", "R", "D"}); + // mio::log_info("Finished postprocessing"); + + mio::unused(sim.statistics_per_timestep().export_csv(mio::path_join("AsymmetricParams_single_run.csv"))); + + return 0; +} diff --git a/cpp/examples/asymmetric_params.cpp b/cpp/examples/asymmetric_params.cpp new file mode 100644 index 0000000000..2b1d2ccd12 --- /dev/null +++ b/cpp/examples/asymmetric_params.cpp @@ -0,0 +1,304 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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/config.h" +#include "memilio/geography/geolocation.h" +#include "memilio/compartments/parameter_studies.h" +#include "memilio/geography/rtree.h" +#include "memilio/geography/distance.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_asymmetric.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_builder.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/utils/index.h" +#include "memilio/utils/logging.h" +#include "memilio/timer/auto_timer.h" +#include "memilio/utils/parameter_distributions.h" +#include "memilio/utils/random_number_generator.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/base_dir.h" +#include "smm/simulation.h" +#include "smm/parameters.h" +#include "fmd/infection_state.h" +#include "fmd/model.h" +#include "fmd/adoption_rates.h" +#include "abm/time.h" +#include "thirdparty/csv.h" +#include +#include + +template +MPI_Datatype mpi_type(); + +template <> +inline MPI_Datatype mpi_type() +{ + return MPI_INT; +} +template <> +inline MPI_Datatype mpi_type() +{ + return MPI_DOUBLE; +} +template <> +inline MPI_Datatype mpi_type() +{ + return MPI_FLOAT; +} +template <> +inline MPI_Datatype mpi_type() +{ + return MPI_LONG_LONG; +} +// If your MPI supports it (MPI-3+), for fixed-width: +#include +template <> +inline MPI_Datatype mpi_type() +{ + return MPI_INT64_T; +} +template <> +inline MPI_Datatype mpi_type() +{ + return MPI_UINT64_T; +} + +template +void bcast_vector(std::vector& v, int root, MPI_Comm comm) +{ + int n = static_cast(v.size()); + MPI_Bcast(&n, 1, MPI_INT, root, comm); + if (n < 0) + std::abort(); + if (n == 0) { + v.clear(); + return; + } + if (v.size() != static_cast(n)) + v.resize(n); + MPI_Bcast(v.data(), n, mpi_type(), root, comm); +} + +int main(int /*argc*/, char** /*argv*/) +{ + mio::mpi::init(); + int size, rank; + MPI_Comm_size(MPI_COMM_WORLD, &size); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + const std::string result_dir = mio::path_join(mio::base_dir(), "example_results"); + + const auto t0 = 0.; + const auto tmax = 100.; + const auto dt = 1.; //initial time step + + const size_t num_runs = 100000; + + using mio::fmd::InfectionState; + using Status = mio::Index; + using mio::regions::Region; + + using Model = mio::smm::Model; + auto home = Region(0); + + auto adoption_rates = mio::fmd::generic_adoption_rates(); + + mio::fmd::Builder builder; + + // usage: + // if (rank==0) v = {1,2,3,4}; + // bcast_vector(v, 0, MPI_COMM_WORLD); + // now every rank has v + std::vector latitudes, longitudes; + std::vector farm_ids, num_cows_vec, dates, num_animals_exchanges; + std::vector froms, tos, from_exchanges, to_exchanges; + if (rank == 0) { + io::CSVReader<4> farms("../../farms200000.csv"); + farms.read_header(io::ignore_extra_column, "farms", "num_cows", "latitude", "longitude"); + int farm_id, num_cows; + double latitude, longitude; + while (farms.read_row(farm_id, num_cows, latitude, longitude)) { + farm_ids.push_back(farm_id); + num_cows_vec.push_back(num_cows); + latitudes.push_back(latitude); + longitudes.push_back(longitude); + } + io::CSVReader<2> edges("../../edges200000.csv"); + edges.read_header(io::ignore_extra_column, "from", "to"); + size_t from, to; + while (edges.read_row(from, to)) { + froms.push_back(from); + tos.push_back(to); + } + io::CSVReader<4> exchanges("../../trade200000.csv"); + exchanges.read_header(io::ignore_extra_column, "date", "num_animals", "from", "to"); + + int date, num_animals; + while (exchanges.read_row(date, num_animals, from, to)) { + dates.push_back(date); + num_animals_exchanges.push_back(num_animals); + from_exchanges.push_back(from); + to_exchanges.push_back(to); + } + if (!mio::create_directory(result_dir)) { + mio::log_error("Could not create result directory \"{}\".", result_dir); + return 1; + } + } + bcast_vector(farm_ids, 0, MPI_COMM_WORLD); + bcast_vector(num_cows_vec, 0, MPI_COMM_WORLD); + bcast_vector(latitudes, 0, MPI_COMM_WORLD); + bcast_vector(longitudes, 0, MPI_COMM_WORLD); + bcast_vector(froms, 0, MPI_COMM_WORLD); + bcast_vector(tos, 0, MPI_COMM_WORLD); + bcast_vector(dates, 0, MPI_COMM_WORLD); + bcast_vector(num_animals_exchanges, 0, MPI_COMM_WORLD); + bcast_vector(from_exchanges, 0, MPI_COMM_WORLD); + bcast_vector(to_exchanges, 0, MPI_COMM_WORLD); + + if (rank == 0) { + mio::log_info("Data broadcasted to all {} ranks.", size); + } + + for (size_t i = 0; i < farm_ids.size(); ++i) { + builder.add_node(farm_ids[i], longitudes[i], latitudes[i], + mio::fmd::create_model(num_cows_vec[i], adoption_rates), t0); + } + auto rng = mio::RandomNumberGenerator(); + + std::vector> interesting_indices; + interesting_indices.push_back( + {Model(Status{InfectionState::Count}, Region(1)).populations.get_flat_index({home, InfectionState::I})}); + for (size_t i = 0; i < froms.size(); ++i) { + builder.add_edge(froms[i], tos[i], interesting_indices); + } + auto graph = std::move(builder).build(); + auto nodes = graph.nodes() | std::views::transform([](const auto& node) { + return &node.property; + }); + + auto tree = mio::geo::RTree(nodes.begin(), nodes.end()); + + std::vector query_distances = {mio::geo::kilometers(5.0)}; + + std::vector>> locally_calculated_neighbors; + size_t num_calculations = farm_ids.size() / size; + // Perform your fair share of rtree queries. + for (size_t i = num_calculations * rank; i < num_calculations * (rank + 1); ++i) { + locally_calculated_neighbors.push_back( + tree.in_range_indices_query(graph.nodes()[i].property.get_location(), query_distances)); + } + size_t num_neighbour_queries = query_distances.size(); + // Look up how many results each query returned. + std::vector locally_saved_neighbour_list_sizes; + for (size_t outer_index = 0; outer_index < num_calculations; ++outer_index) { + for (size_t inner_index = 0; inner_index < num_neighbour_queries; ++inner_index) { + locally_saved_neighbour_list_sizes.push_back(locally_calculated_neighbors[outer_index][inner_index].size()); + } + } + int local_num_total_neighbours = + std::accumulate(locally_saved_neighbour_list_sizes.begin(), locally_saved_neighbour_list_sizes.end(), 0); + std::vector total_neighbours_per_rank(size, 0); + // Share with everybody the total number of all your query results. + MPI_Allgather(&local_num_total_neighbours, 1, MPI_INT, total_neighbours_per_rank.data(), 1, MPI_INT, + MPI_COMM_WORLD); + std::vector displacements(size, 0); + for (int i = 1; i < size; ++i) { + displacements[i] += displacements[i - 1]; + displacements[i] += total_neighbours_per_rank[i - 1]; + } + //Share with everybody the actual lists with query result counts. + std::vector all_neighbour_list_sizes(size * num_calculations * num_neighbour_queries); + MPI_Allgather(locally_saved_neighbour_list_sizes.data(), num_calculations * num_neighbour_queries, MPI_INT, + all_neighbour_list_sizes.data(), num_calculations * num_neighbour_queries, MPI_INT, MPI_COMM_WORLD); + std::vector flat_local; + for (const auto& nodes : locally_calculated_neighbors) { + for (const auto& neighbours : nodes) { + flat_local.insert(flat_local.end(), neighbours.begin(), neighbours.end()); + } + } + int total_elements = std::accumulate(all_neighbour_list_sizes.begin(), all_neighbour_list_sizes.end(), 0); + std::vector flat_global(total_elements); + //Share with everybody the actual lists with query results. + MPI_Allgatherv(flat_local.data(), flat_local.size(), MPI_INT, flat_global.data(), total_neighbours_per_rank.data(), + displacements.data(), MPI_INT, MPI_COMM_WORLD); + // Reconstruct the neighbor lists for all nodes. + size_t big_index = 0; + size_t pos = 0; + for (size_t index = 0; index < size * num_calculations; ++index) { + std::vector> neighbors; + for (size_t i = 0; i < num_neighbour_queries; ++i) { + std::vector current_neighbors(all_neighbour_list_sizes[big_index]); + for (int j = 0; j < all_neighbour_list_sizes[big_index]; ++j) { + current_neighbors[j] = flat_global[pos]; + pos++; + } + neighbors.push_back(current_neighbors); + big_index++; + } + graph.nodes()[index].property.set_regional_neighbors(neighbors); + } + // Calculate the nodes that were left out if the number of nodes is not divisible by the number of ranks. + for (auto index = num_calculations * size; index < farm_ids.size(); ++index) { + graph.nodes()[index].property.set_regional_neighbors( + tree.in_range_indices_query(graph.nodes()[index].property.get_location(), query_distances)); + } + + auto sim = mio::make_mobility_sim(t0, dt, std::move(graph)); + + for (size_t i = 0; i < dates.size(); ++i) { + sim.add_exchange(dates[i], num_animals_exchanges[i], from_exchanges[i], to_exchanges[i]); + } + if (rank == 0) { + mio::log_info("Starting the study:"); + } + mio::ParameterStudy study(sim, t0, tmax, dt, num_runs); + + auto get_simulation = [](const auto& sim, ScalarType, ScalarType, size_t run_id) { + auto sim2 = sim; + sim2.get_rng() = mio::thread_local_rng(); + int index_farm = mio::UniformIntDistribution::get_instance()(sim2.get_rng(), 0, + int(sim2.get_graph().nodes().size() - 1)); + auto E_index = sim2.get_graph().nodes()[0].property.get_simulation().get_model().populations.get_flat_index( + {Region(0), InfectionState::E}); + auto S_index = sim2.get_graph().nodes()[0].property.get_simulation().get_model().populations.get_flat_index( + {Region(0), InfectionState::S}); + auto num_sus = sim2.get_graph().nodes()[index_farm].property.get_result().get_last_value()[S_index]; + auto num_exposed = std::ceil(0.01 * num_sus); + sim2.get_graph().nodes()[index_farm].property.get_result().get_last_value()[E_index] = num_exposed; + sim2.get_graph().nodes()[index_farm].property.get_result().get_last_value()[S_index] = num_sus - num_exposed; + mio::log_info("Run {}: Infecting {} individuals at node {}.", run_id, num_exposed, index_farm); + return sim2; + }; + auto handle_result = [](auto&& sim, auto&& run) { + auto infection_numbers = sim.sum_nodes(); + auto dead_index = sim.get_graph().nodes()[0].property.get_simulation().get_model().populations.get_flat_index( + {Region(0), InfectionState::D}); + if (infection_numbers[dead_index] > 50) { + auto abs_path = mio::path_join(mio::base_dir(), "example_results"); + auto result = sim.statistics_per_timestep().export_csv( + mio::path_join(abs_path, "AsymmetricParams_run" + std::to_string(run) + ".csv")); + } + return 0; + }; + auto result = study.run(get_simulation, handle_result); + mio::mpi::finalize(); + return 0; +} diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 570ce17c6b..89dd7185b3 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: Julia Bicker, René Schmieding +* Authors: Julia Bicker, René Schmieding, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -18,11 +18,12 @@ * limitations under the License. */ +#include "memilio/config.h" +#include "memilio/epidemiology/age_group.h" #include "smm/simulation.h" #include "smm/model.h" #include "smm/parameters.h" #include "memilio/data/analyze_result.h" -#include "memilio/epidemiology/adoption_rate.h" enum class InfectionState { @@ -33,61 +34,73 @@ enum class InfectionState R, D, Count +}; +struct Species : public mio::Index { + Species(size_t val) + : Index(val) + { + } }; int main() { + using Age = mio::AgeGroup; + using Status = mio::Index; + using mio::regions::Region; + using enum InfectionState; //Example how to run the stochastic metapopulation models with four regions - const size_t num_regions = 4; - using Model = mio::smm::Model; + const size_t num_regions = 4; + const size_t num_age_groups = 1; + const size_t num_species = 1; + using Model = mio::smm::Model; ScalarType numE = 12, numC = 4, numI = 12, numR = 0, numD = 0; - Model model; + Model model(Status{Count, Age(num_age_groups), Species(num_species)}, Region(num_regions)); //Population are distributed uniformly to the four regions for (size_t r = 0; r < num_regions; ++r) { - model.populations[{mio::regions::Region(r), InfectionState::S}] = - (1000 - numE - numC - numI - numR - numD) / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::E}] = numE / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::C}] = numC / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::I}] = numI / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::R}] = numR / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::D}] = numD / num_regions; + model.populations[{Region(r), S, Age(0), Species(0)}] = (1000 - numE - numC - numI - numR - numD) / num_regions; + model.populations[{Region(r), E, Age(0), Species(0)}] = numE / num_regions; + model.populations[{Region(r), C, Age(0), Species(0)}] = numC / num_regions; + model.populations[{Region(r), I, Age(0), Species(0)}] = numI / num_regions; + model.populations[{Region(r), R, Age(0), Species(0)}] = numR / num_regions; + model.populations[{Region(r), D, Age(0), Species(0)}] = numD / num_regions; } + using AR = mio::smm::AdoptionRates; + using TR = mio::smm::TransitionRates; + //Set infection state adoption and spatial transition rates - std::vector> adoption_rates; - std::vector> transition_rates; + AR::Type adoption_rates; + TR::Type transition_rates; for (size_t r = 0; r < num_regions; ++r) { - adoption_rates.push_back({InfectionState::S, - InfectionState::E, - mio::regions::Region(r), + adoption_rates.push_back({{S, Age(0), Species(0)}, + {E, Age(0), Species(0)}, + Region(r), 0.1, - {{InfectionState::C, 1}, {InfectionState::I, 0.5}}}); - adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}}); + {{{C, Age(0), Species(0)}, 1}, {{I, Age(0), Species(0)}, 0.5}}}); + adoption_rates.push_back({{C, Age(0), Species(0)}, {R, Age(0), Species(0)}, Region(r), 0.2 / 3., {}}); + adoption_rates.push_back({{E, Age(0), Species(0)}, {C, Age(0), Species(0)}, Region(r), 1.0 / 5., {}}); + adoption_rates.push_back({{C, Age(0), Species(0)}, {I, Age(0), Species(0)}, Region(r), 0.8 / 3., {}}); + adoption_rates.push_back({{I, Age(0), Species(0)}, {R, Age(0), Species(0)}, Region(r), 0.99 / 5., {}}); + adoption_rates.push_back({{I, Age(0), Species(0)}, {D, Age(0), Species(0)}, Region(r), 0.01 / 5., {}}); } //Agents in infection state D do not transition - for (size_t s = 0; s < static_cast(InfectionState::D); ++s) { + for (size_t s = 0; s < static_cast(D); ++s) { for (size_t i = 0; i < num_regions; ++i) { for (size_t j = 0; j < num_regions; ++j) if (i != j) { - transition_rates.push_back( - {InfectionState(s), mio::regions::Region(i), mio::regions::Region(j), 0.01}); - transition_rates.push_back( - {InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01}); + transition_rates.push_back({{InfectionState(s), Age(0), Species(0)}, Region(i), Region(j), 0.01}); + transition_rates.push_back({{InfectionState(s), Age(0), Species(0)}, Region(j), Region(i), 0.01}); } } } - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get() = adoption_rates; + model.parameters.get() = transition_rates; ScalarType dt = 0.1; ScalarType tmax = 30.0; diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt index 658733a7f8..45722d358b 100644 --- a/cpp/memilio/CMakeLists.txt +++ b/cpp/memilio/CMakeLists.txt @@ -77,6 +77,8 @@ add_library(memilio mobility/metapopulation_mobility_instant.cpp mobility/metapopulation_mobility_stochastic.h mobility/metapopulation_mobility_stochastic.cpp + mobility/metapopulation_mobility_asymmetric.h + mobility/metapopulation_mobility_asymmetric.cpp mobility/graph_simulation.h mobility/graph_simulation.cpp mobility/graph.h diff --git a/cpp/memilio/compartments/parameter_studies.h b/cpp/memilio/compartments/parameter_studies.h index abe407d5ba..4cfad60a90 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -31,6 +31,7 @@ #include #include +#include #include #include #include diff --git a/cpp/memilio/data/analyze_result.h b/cpp/memilio/data/analyze_result.h index 5e80a05773..55964c37c2 100644 --- a/cpp/memilio/data/analyze_result.h +++ b/cpp/memilio/data/analyze_result.h @@ -335,6 +335,31 @@ std::vector> ensemble_percentile(const std::vector +TimeSeries ensemble_percentile(const std::vector>& ensemble_result, FP p) +{ + assert(p > 0.0 && p < 1.0 && "Invalid percentile value."); + + auto num_runs = ensemble_result.size(); + auto num_time_points = ensemble_result[0].get_num_time_points(); + auto num_elements = ensemble_result[0].get_num_elements(); + + TimeSeries percentile = TimeSeries::zero(num_time_points, num_elements); + + std::vector single_element(num_runs); //reused for each element + for (Eigen::Index time = 0; time < num_time_points; time++) { + percentile.get_time(time) = ensemble_result[0].get_time(time); + for (Eigen::Index elem = 0; elem < num_elements; elem++) { + std::transform(ensemble_result.begin(), ensemble_result.end(), single_element.begin(), [=](auto& run) { + return run[time][elem]; + }); + std::sort(single_element.begin(), single_element.end()); + percentile[time][elem] = single_element[static_cast(num_runs * p)]; + } + } + return percentile; +} + template FP result_distance_2norm(const std::vector>& result1, const std::vector>& result2) diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index 705f729025..8f72d650d2 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: René Schmieding, Julia Bicker +* Authors: René Schmieding, Julia Bicker, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -20,9 +20,8 @@ #ifndef MIO_EPI_ADOPTIONRATE_H #define MIO_EPI_ADOPTIONRATE_H -#include "memilio/utils/index.h" -#include "memilio/config.h" #include "memilio/geography/regions.h" +#include namespace mio { @@ -30,7 +29,7 @@ namespace mio /** * @brief Struct defining an influence for a second-order adoption. * The population having "status" is multiplied with "factor." - * @tparam Status An infection state enum. + * @tparam Status An infection state enum or MultiIndex. */ template struct Influence { @@ -43,15 +42,16 @@ struct Influence { * The AdoptionRate is considered to be of second-order if there are any "influences". * In the d_abm and smm simulations, "from" is implicitly an influence, scaled by "factor". This is multiplied by * the sum over all "influences", which scale their "status" with the respective "factor". - * @tparam Status An infection state enum. + * @tparam Status An infection state enum or MultiIndex. + * @tparam Region A MultiIndex. */ -template +template struct AdoptionRate { Status from; // i Status to; // j - mio::regions::Region region; // k + Region region; // k FP factor; // gammahat_{ij}^k - std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) + std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )}; }; } // namespace mio diff --git a/cpp/memilio/epidemiology/populations.h b/cpp/memilio/epidemiology/populations.h index 3c300bd79a..88a87bd1a0 100644 --- a/cpp/memilio/epidemiology/populations.h +++ b/cpp/memilio/epidemiology/populations.h @@ -307,6 +307,18 @@ class Populations : public CustomIndexArray, Categories...> } }; +/** + * @brief Population template specialization, forwarding categories from a MultiIndex to the Population. + * + * @tparam FP A floating point type, e.g., double. + * @tparam Categories Index categories. + */ +template +class Populations> : public Populations +{ + using Populations::Populations; +}; + } // namespace mio #endif // MIO_EPI_POPULATIONS_H diff --git a/cpp/memilio/geography/rtree.h b/cpp/memilio/geography/rtree.h index 740cf49878..211967d0f3 100644 --- a/cpp/memilio/geography/rtree.h +++ b/cpp/memilio/geography/rtree.h @@ -46,6 +46,12 @@ concept IsSphericalLocation = requires(const Location& loc) { std::is_floating_point_v && std::is_floating_point_v; }; +template +concept SphericalLocationIteratorType = std::indirectly_readable && requires(const Iter& loc) { + std::is_floating_point_vget_latitude())> && + std::is_floating_point_vget_longitude())>; +}; + /** * @brief R-tree for spatial queries of geographical locations on the sphere. * @@ -79,6 +85,26 @@ class RTree } } + /** + * @brief Construct a new RTree object with data given in a range + * + * @param first The beginning of the range + * @param last The end of the range + * The provided location data needs to provide get_latitude() and get_longitude(). + */ + template + RTree(Iter first, Iter last) + : rtree{} + { + size_t index = 0; + while (first != last) { + Point point((*first)->get_longitude(), (*first)->get_latitude()); + rtree.insert(Node(point, index)); + ++first; + ++index; + } + } + /** * @brief Return the number of data points stored in the RTree. * diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index a7a60330e0..87a6b16a5b 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -20,6 +20,7 @@ #ifndef GRAPH_H #define GRAPH_H +#include "memilio/timer/auto_timer.h" #include "memilio/utils/stl_util.h" #include "memilio/epidemiology/age_group.h" #include "memilio/utils/date.h" @@ -27,6 +28,7 @@ #include "memilio/utils/parameter_distributions.h" #include "memilio/epidemiology/damping.h" #include "memilio/geography/regions.h" +#include "memilio/utils/logging.h" #include #include #include @@ -237,6 +239,7 @@ class Graph requires std::constructible_from void add_edge(size_t start_node_idx, size_t end_node_idx, Args&&... args) { + mio::timing::AutoTimer<"Graph.add_edge()"> timer; assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); insert_sorted_replace(m_edges, Edge(start_node_idx, end_node_idx, std::forward(args)...), [](auto&& e1, auto&& e2) { @@ -245,6 +248,42 @@ class Graph }); } + template + Edge& lazy_add_edge(size_t start_node_idx, size_t end_node_idx, Args&&... args) + { + mio::timing::AutoTimer<"Graph.lazy_add_edge()"> timer; + assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); + m_edges.emplace_back(start_node_idx, end_node_idx, std::forward(args)...); + return m_edges.back(); + } + + Edge& sort_edges() + { + mio::timing::AutoTimer<"Graph.sort_edges()"> timer; + std::sort(m_edges.begin(), m_edges.end(), [](auto&& e1, auto&& e2) { + return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx + : e1.start_node_idx < e2.start_node_idx; + }); + return m_edges.back(); + } + + Edge& make_edges_unique() + { + mio::timing::AutoTimer<"Graph.make_edges_unique()"> timer; + std::vector> unique_edges; + unique_edges.reserve(m_edges.size()); + std::ranges::unique_copy(m_edges, std::back_inserter(unique_edges), [](auto&& e1, auto&& e2) { + return e1.start_node_idx == e2.start_node_idx && e1.end_node_idx == e2.end_node_idx; + }); + m_edges = std::move(unique_edges); + return m_edges.back(); + } + + void reserve_edges(size_t n) + { + m_edges.reserve(n); + } + /** * @brief range of nodes */ @@ -278,7 +317,9 @@ class Graph } /** - * @brief range of edges going out from a specific node + * @brief Range of edges going out from a specific node + * + * @param node_idx ID of node */ auto out_edges(size_t node_idx) { @@ -287,12 +328,26 @@ class Graph /** * @brief range of edges going out from a specific node + * + * @param node_idx ID of node */ auto out_edges(size_t node_idx) const { return out_edges(begin(m_edges), end(m_edges), node_idx); } + auto& get_edge(size_t node_a, size_t node_b) + { + auto edges = out_edges(node_a); + for (auto edge = edges.begin(); edge != edges.end(); ++edge) { + if (edge->end_node_idx == node_b) { + return *edge; + } + } + mio::log_error("Edge from {} to {} not found.", node_a, node_b); + throw std::out_of_range("Edge not found"); + } + private: template static auto out_edges(Iter b, Iter e, size_t idx) diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 0ceb33a821..1264773ec3 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -20,14 +20,150 @@ #ifndef MIO_MOBILITY_GRAPH_SIMULATION_H #define MIO_MOBILITY_GRAPH_SIMULATION_H +#include "memilio/config.h" #include "memilio/mobility/graph.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/logging.h" #include "memilio/utils/random_number_generator.h" +#include #include "memilio/compartments/feedback_simulation.h" #include "memilio/geography/regions.h" +#include "models/smm/parameters.h" namespace mio { +class MobilityParametersTimed +{ + +public: + /** + * @brief Construct a new Mobility Parameters Timed object + * + */ + MobilityParametersTimed() + : _exchanges{} {}; + MobilityParametersTimed(std::filebuf& input_data) + : _exchanges{} + { + insert_input_data(input_data); + }; + MobilityParametersTimed(double time, double number, size_t from, size_t to) + : _exchanges{} + { + _exchanges.push(ExchangeData(time, number, from, to)); + }; + + void add_exchange(double time, double number, size_t from, size_t to) + { + _exchanges.push(ExchangeData{time, number, from, to}); + } + + /** + * @brief Return the number of exchanged items in the next exchange event. + * + * @return auto + */ + auto next_event_number() + { + return _exchanges.top().number; + }; + /** + * @brief Return the timepoint of the next exchangeclass EdgePropertyT event. + * + * @return auto + */ + auto next_event_time() const + { + if (_exchanges.empty()) { + return std::numeric_limits::max(); + } + return _exchanges.top().time; + }; + /** + * @brief Return the destination node id of the next exchange + * + * @return auto + */ + auto next_event_edge_id() + { + return _exchanges.top().edge_id; + } + auto next_event_to_id() + { + return _exchanges.top().to; + } + auto next_event_from_id() + { + return _exchanges.top().from; + } + /** + * @brief Return a const reference to the next event + * + * @return auto + */ + auto next_event() + { + return _exchanges.top(); + } + /** + * @brief Delete the next event from the heap + * + * @return auto + */ + auto pop_next_event() + { + return _exchanges.pop(); + } + /** + * @brief Return the ExchangeData for the next exchange event and delete it from the list. + * + * @return auto + */ + auto process_next_event() + { + auto next_event = _exchanges.top(); + _exchanges.pop(); + return next_event; + }; + + auto size() const + { + return _exchanges.size(); + } + +private: + void insert_input_data(std::filebuf& input_data) + { + //... + mio::unused(input_data); + }; + + /** + * @brief Stores Timepoint and number of exchanged items for an exchange process. + * + * @param time Timepoint of the exchange process + * @param number Number of exchanged items + */ + struct ExchangeData { + double time; + double number; + size_t from; + size_t to; + int edge_id{}; + }; + + struct CompareExchangeData { + bool operator()(const ExchangeData& left, const ExchangeData& right) + { + return left.time > right.time; + }; + }; + +private: + std::priority_queue, CompareExchangeData> _exchanges; +}; + /** * @brief abstract simulation on a graph with alternating node and edge actions */ @@ -127,11 +263,10 @@ class GraphSimulationStochastic typename Graph::NodeProperty&, typename Graph::NodeProperty&)>, std::function> { - using Base = GraphSimulationBase, - std::function>; - + std::function>; using node_function = typename Base::node_function; using edge_function = typename Base::edge_function; @@ -256,6 +391,369 @@ class GraphSimulationStochastic RandomNumberGenerator m_rng; }; +template +class AsymmetricGraphSimulation : public GraphSimulationBase +{ + using Base = GraphSimulationBase; + using Base::GraphSimulationBase; + +public: + void advance(Timepoint t_max = 1.0) + { + mio::timing::AutoTimer<"Graph Simulation Advance"> timer; + auto dt = m_parameters.next_event_time() - Base::m_t; + while (Base::m_t < t_max) { + mio::log_debug("Time: {}", Base::m_t); + if (Base::m_t + dt > t_max) { + dt = t_max - Base::m_t; + } + + cull(dt); + vaccinate(dt); + + for (auto& n : Base::m_graph.nodes()) { + Base::m_node_func(Base::m_t, dt, n.property); + } + + Base::m_t += dt; + + while (m_parameters.next_event_time() == Base::m_t) { + auto next_event = m_parameters.process_next_event(); + if (Base::m_graph.nodes()[next_event.from].property.is_quarantined() || + Base::m_graph.nodes()[next_event.to].property.is_quarantined()) { + mio::log_debug("Mobility from node {} to node {} at time {} skipped due to quarantine.", + next_event.from, next_event.to, Base::m_t); + continue; + } + auto& e = Base::m_graph.get_edge(next_event.from, next_event.to); + mio::log_debug("{}, {}", e.start_node_idx, e.end_node_idx); + Base::m_edge_func(Base::m_t, next_event.number, e.property, + Base::m_graph.nodes()[e.start_node_idx].property, + Base::m_graph.nodes()[e.end_node_idx].property, m_rng); + } + dt = m_parameters.next_event_time() - Base::m_t; + apply_interventions(); + } + } + + /** + * @brief Go through the culling queue and cull as many animals as possible within dt. + * + * @param dt Time span for culling. + */ + void cull(ScalarType dt) + { + auto capacity = culling_capacity_per_day * dt; + while (!culling_queue.empty() && capacity > 0) { + auto [node_id, day] = culling_queue.front(); + auto animals = Base::m_graph.nodes()[node_id].property.get_result().get_last_value(); + auto num_animals = std::accumulate(animals.begin(), animals.end() - 1, 0); + if (num_animals <= capacity) { + mio::log_debug("Culling {} animals at node {} starting on day {} and going on for {} days.", + num_animals, node_id, Base::m_t, dt); + for (auto index = 0; index < (animals.size() - 1); index++) { + animals[index] = 0; + } + animals[animals.size() - 1] += num_animals; + culling_queue.pop(); + capacity -= num_animals; + } + else { + while (capacity > 0) { + mio::log_debug("Culling {} animals of node {} on day {} for {} days.", capacity, node_id, Base::m_t, + dt); + for (auto index = 0; index < (animals.size() - 1); index++) { + if (animals[index] < capacity) { + capacity -= animals[index]; + animals[animals.size() - 1] += animals[index]; + animals[index] = 0; + } + else { + animals[animals.size() - 1] += capacity; + animals[index] -= capacity; + capacity = 0; + break; + } + } + } + } + } + } + + /** + * @brief Go through the vaccination queue and vaccinate as many animals as possible within dt. + * + * @param dt Time span for vaccination. + */ + void vaccinate(ScalarType dt) + { + auto capacity = vaccination_capacity_per_day * dt; + while (!vaccination_queue.empty() && capacity > 0) { + auto [node_id, day] = vaccination_queue.front(); + auto animals = Base::m_graph.nodes()[node_id].property.get_result().get_last_value(); + auto num_animals = std::accumulate(animals.begin(), animals.end() - 1, 0); + if (num_animals <= capacity) { + mio::log_debug("Vaccinating {} animals at node {} starting on day {} and going on for {} days.", + num_animals, node_id, Base::m_t, dt); + for (auto index = 0; index < (animals.size() - 2); index++) { + animals[index] = 0; + } + animals[animals.size() - 2] += num_animals; + vaccination_queue.pop(); + capacity -= num_animals; + } + else { + while (capacity > 0) { + for (auto index = 0; index < (animals.size() - 1); index++) { + if (animals[index] < capacity) { + capacity -= animals[index]; + animals[animals.size() - 2] += animals[index]; + animals[index] = 0; + } + else { + animals[animals.size() - 2] += capacity; + animals[index] -= capacity; + capacity = 0; + break; + } + } + } + } + } + } + + void add_exchange(double time, double number, size_t from, size_t to) + { + m_parameters.add_exchange(time, number, from, to); + } + + /** + * get the mobility parameters. + */ + const MobilityParametersTimed& get_parameters() const + { + return m_parameters; + } + + auto sum_exchanges() + { + const auto size = Base::m_graph.edges()[0].property.get_mobility_results().get_num_elements(); + std::vector results(size, 0.0); + for (auto& n : Base::m_graph.edges()) { + assert(n.property.get_mobility_results().get_num_elements() == size); + for (auto result = n.property.get_mobility_results().begin(); + result != n.property.get_mobility_results().end(); ++result) { + for (int i = 0; i < size; i++) { + results[i] += (*result)[i]; + } + } + } + return results; + } + + auto exchanges_per_timestep() + { + const auto size = Base::m_graph.edges()[0].property.get_mobility_results().get_num_elements(); + std::vector timepoints; + // Collect all exchange timepoints + + for (auto& n : Base::m_graph.edges()) { + auto local_timepoints = n.property.get_mobility_results().get_time_points(); + for (auto t : local_timepoints) { + if (std::find(timepoints.begin(), timepoints.end(), t) == timepoints.end()) { + timepoints.push_back(t); + } + } + } + std::sort(timepoints.begin(), timepoints.end()); + auto results = TimeSeries::zero(timepoints.size(), size, timepoints); + for (auto& n : Base::m_graph.edges()) { + assert(n.property.get_mobility_results().get_num_elements() == size); + auto edge_result = n.property.get_mobility_results(); + // Add exchange data to big TimeSeries + for (Eigen::Index time_index = 0; time_index < edge_result.get_num_time_points(); ++time_index) { + auto time = edge_result.get_time(time_index); + auto index = results.get_index_of_time(time); + if (index == Eigen::Index(-1)) { + continue; + } + for (int i = 0; i < size; i++) { + results.get_value(index)[i] += edge_result.get_value(time_index)[i]; + } + } + } + return results; + } + + auto sum_nodes() + { + const auto size = Base::m_graph.nodes()[0].property.get_result().get_num_elements(); + std::vector results(size, 0.0); + for (auto& n : Base::m_graph.nodes()) { + assert(n.property.get_result().get_num_elements() == size); + for (int i = 0; i < size; i++) { + results[i] += n.property.get_result().get_last_value()[i]; + } + } + return results; + } + + auto statistics_per_timestep() + { + const auto size = Base::m_graph.nodes()[0].property.get_result().get_num_elements(); + std::vector timepoints; + // Collect all exchange timepoints => All simulations are stopped and write results at those + + for (auto& n : Base::m_graph.edges()) { + auto local_timepoints = n.property.get_mobility_results().get_time_points(); + for (auto t : local_timepoints) { + if (std::find(timepoints.begin(), timepoints.end(), t) == timepoints.end()) { + timepoints.push_back(t); + } + } + } + std::sort(timepoints.begin(), timepoints.end()); + auto results = TimeSeries::zero(timepoints.size(), size, timepoints); + for (auto& n : Base::m_graph.nodes()) { + assert(n.property.get_result().get_num_elements() == size); + auto node_timeseries = n.property.get_result(); + for (Eigen::Index time_index = 0; time_index < node_timeseries.get_num_time_points(); ++time_index) { + auto time = node_timeseries.get_time(time_index); + auto index = results.get_index_of_time(time); + if (index == Eigen::Index(-1)) { + continue; + } + for (int i = 0; i < size; i++) { + results.get_value(index)[i] += node_timeseries.get_value(time_index)[i]; + } + } + } + return results; + } + + auto return_all_time_series() + { + std::vector> all_time_series; + for (auto& n : Base::m_graph.nodes()) { + all_time_series.push_back(n.property.get_result()); + } + return all_time_series; + } + + auto statistics_per_timestep(std::vector node_indices) + { + assert(node_indices.size() > 0); + const auto size = Base::m_graph.nodes()[node_indices[0]].property.get_result().get_num_elements(); + std::vector timepoints; + // Collect all exchange timepoints => All simulations are stopped and write results at those + for (auto& n : Base::m_graph.edges()) { + auto local_timepoints = n.property.get_mobility_results().get_time_points(); + for (auto t : local_timepoints) { + if (std::find(timepoints.begin(), timepoints.end(), t) == timepoints.end()) { + timepoints.push_back(t); + } + } + } + std::sort(timepoints.begin(), timepoints.end()); + auto results = TimeSeries::zero(timepoints.size(), size, timepoints); + for (size_t node_index : node_indices) { + auto node_timeseries = Base::m_graph.nodes()[node_index].property.get_result(); + assert(node_timeseries.get_num_elements() == size); + for (Eigen::Index time_index = 0; time_index < node_timeseries.get_num_time_points(); ++time_index) { + auto time = node_timeseries.get_time(time_index); + auto index = results.get_index_of_time(time); + if (index == Eigen::Index(-1)) { + continue; + } + for (int i = 0; i < size; i++) { + results.get_value(index)[i] += node_timeseries.get_value(time_index)[i]; + } + } + } + return results; + } + + void apply_interventions() + { + // Set quarantine + // for (auto& n : Base::m_graph.nodes()) { + // if (n.property.get_result().get_last_value()[2] > 40) { + // mio::log_debug("Node {} is quarantined at time {}.", n.id, Base::m_t); + // n.property.set_quarantined(true); + // } + // else { + // mio::log_debug("Node {} is not quarantined at time {}.", n.id, Base::m_t); + // n.property.set_quarantined(false); + // } + // } + auto infectious_compartments = {1, 2, 3, 4}; + for (auto& n : Base::m_graph.nodes()) { + auto total_infections = 0; + for (auto comp : infectious_compartments) { + total_infections += n.property.get_result().get_last_value()[comp]; + } + if (total_infections > 1) { + mio::log_debug("Node {} spreads higher risk category at time {} because there are {} total infections.", + n.id, Base::m_t, total_infections); + mio::log_debug("This will affect {} other farms.", n.property.get_regional_neighbors()[0].size()); + for (size_t index = 0; index < n.property.get_regional_neighbors()[0].size(); ++index) { + const auto node_id = n.property.get_regional_neighbors()[0][index]; + mio::log_debug("Neighbor node id: {}", node_id); + auto neighbour_property = Base::m_graph.nodes()[node_id].property; + using Model = std::decay_t; + using AdoptionRate = + mio::smm::AdoptionRates; + neighbour_property.get_simulation().get_model().parameters.template get()[1].factor += + infectionrisk; + } + } + if (total_infections > 4) { + if (first_detection > Base::m_t) { + first_detection = Base::m_t; + } + mio::log_debug("Node {} is culled at time {} because there are {} total infections.", n.id, Base::m_t, + total_infections); + cull_node(n.id); + for (size_t index = 0; index < n.property.get_regional_neighbors()[0].size(); ++index) { + const auto node_id = n.property.get_regional_neighbors()[0][index]; + vaccinate(node_id); + } + } + } + } + + void cull_node(size_t node_id) + { + culling_queue.push(std::make_pair(node_id, Base::m_t)); + Base::m_graph.nodes()[node_id].property.set_quarantined(true); + mio::log_debug("It is day {} and we already want to cull {} farms.", Base::m_t, culling_queue.size()); + } + + void vaccinate_node(size_t node_id) + { + vaccination_queue.push(std::make_pair(node_id, Base::m_t)); + } + + RandomNumberGenerator& get_rng() + { + return m_rng; + } + + ScalarType infectionrisk = 0.0; + +private: + mio::MobilityParametersTimed m_parameters; + RandomNumberGenerator m_rng; + std::queue> culling_queue; + ScalarType culling_capacity_per_day = 2000; + std::queue> vaccination_queue; + ScalarType vaccination_capacity_per_day = 500; + ScalarType first_detection = std::numeric_limits::max(); +}; + template auto make_graph_sim(Timepoint t0, Timespan dt, Graph&& g, NodeF&& node_func, EdgeF&& edge_func) { diff --git a/cpp/memilio/mobility/metapopulation_mobility_asymmetric.cpp b/cpp/memilio/mobility/metapopulation_mobility_asymmetric.cpp new file mode 100644 index 0000000000..eb6bd917fc --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_asymmetric.cpp @@ -0,0 +1,24 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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/mobility/metapopulation_mobility_asymmetric.h" + +namespace mio +{ +} // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_asymmetric.h b/cpp/memilio/mobility/metapopulation_mobility_asymmetric.h new file mode 100644 index 0000000000..9a98d723f0 --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_asymmetric.h @@ -0,0 +1,263 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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. +*/ +#ifndef METAPOPULATION_MOBILITY_STOCHASTIC_H +#define METAPOPULATION_MOBILITY_STOCHASTIC_H + +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/random_number_generator.h" +#include "memilio/geography/geolocation.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" + +#include +#include +#include +#include +#include + +namespace mio +{ + +template +class LocationNode : public SimulationNode +{ + using Base = SimulationNode; + +public: + template ::value, void>> + LocationNode(FP latitude, FP longitude, Args&&... args) + : Base(std::forward(args)...) + , m_location(latitude, longitude) + , regional_neighbor_indices{} + { + } + + auto get_location() const + { + return m_location; + } + + void set_location(FP latitude, FP longitude) + { + m_location = mio::geo::GeographicalLocation(latitude, longitude); + } + + FP get_longitude() const + { + return m_location.get_longitude(); + } + + FP get_latitude() const + { + return m_location.get_latitude(); + } + + void set_regional_neighbors(const std::vector>& neighbors) + { + regional_neighbor_indices = neighbors; + } + + auto get_regional_neighbors() const + { + return regional_neighbor_indices; + } + + auto is_quarantined() const + { + return m_is_quarantined; + } + + void set_quarantined(bool quarantine) + { + m_is_quarantined = quarantine; + } + +private: + mio::geo::GeographicalLocation m_location; // location of the node + std::vector> regional_neighbor_indices; + bool m_is_quarantined{false}; +}; + +/** + * node functor for mobility-based simulation. + * @see SimulationNode::advance + */ +template +void advance_model(FP t, FP dt, LocationNode& node) +{ + node.advance(t, dt); +} + +/** + * represents the mobility between two nodes. + */ +template +class MobilityEdgeDirected +{ +public: + /** + * create edge with timed movement parameters. + * @param params mobility rate for each group and compartment + */ + + MobilityEdgeDirected(size_t size) + : m_mobility_results(size) + { + } + + MobilityEdgeDirected(size_t size, const std::vector>& saved_compartment_indices) + : m_mobility_results(size) + , m_saved_compartment_indices(saved_compartment_indices) + { + } + + MobilityEdgeDirected(const std::vector>& saved_compartment_indices) + : m_mobility_results(saved_compartment_indices.size() + 1) + , m_saved_compartment_indices(saved_compartment_indices) + { + } + + /** + * @brief Get the count of exchanges in selected compartments, along with the total number of exchanges. + * + * @return A reference to the TimeSeries object representing the mobility results. + */ + TimeSeries& get_mobility_results() + { + return m_mobility_results; + } + const TimeSeries& get_mobility_results() const + { + return m_mobility_results; + } + + /** + * compute mobility from node_from to node_to for a given event + * @param[in] event index specifying which compartment and age group change nodes + * @param node_from node that people changed from + * @param node_to node that people changed to + */ + template + void apply_mobility(const FP t, const FP num_moving, LocationNode& node_from, + LocationNode& node_to, mio::RandomNumberGenerator& rng); + +private: + // MobilityParametersTimed m_parameters; + TimeSeries m_mobility_results; + std::vector> m_saved_compartment_indices; + + void add_mobility_result_time_point(const FP t, std::vector& travellers) + { + const size_t save_indices_size = this->m_saved_compartment_indices.size(); + if (save_indices_size > 0) { + + Eigen::VectorXd condensed_values = Eigen::VectorXd::Zero(save_indices_size + 1); + + // sum up the values of m_saved_compartment_indices for each group (e.g. Age groups) + std::transform(this->m_saved_compartment_indices.begin(), this->m_saved_compartment_indices.end(), + condensed_values.data(), [&travellers](const auto& indices) { + return std::accumulate(indices.begin(), indices.end(), 0.0, + [&travellers](FP sum, auto i) { + return sum + travellers[i]; + }); + }); + + // the last value is the sum of commuters + condensed_values[save_indices_size] = std::accumulate(travellers.begin(), travellers.end(), 0.0); + + // Move the condensed values to the m_mobility_results time series + m_mobility_results.add_time_point(t, std::move(condensed_values)); + } + } +}; + +template +template +void MobilityEdgeDirected::apply_mobility(const FP t, const FP num_moving, LocationNode& node_from, + LocationNode& node_to, mio::RandomNumberGenerator& rng) +{ + // auto next_event = m_parameters.process_next_event(); + // auto num_moving = next_event.number; + // auto num_available = boost::numeric::ublas::sum(node_from.get_result().get_last_value()); + auto distribution = DiscreteDistributionInPlace(); + std::vector travellers(node_from.get_result().get_last_value().size(), 0); + if (num_moving > std::accumulate(node_from.get_result().get_last_value().begin(), + node_from.get_result().get_last_value().end(), 0.0)) { + mio::log_warning("Trying to move more individuals than available ({}) at time {}.", num_moving, t); + } + else { + for (int i = 0; i < num_moving; ++i) { + auto group = distribution(rng, {node_from.get_result().get_last_value()}); + node_from.get_result().get_last_value()[group] -= 1; + travellers[group] += 1; + node_to.get_result().get_last_value()[group] += 1; + } + } + add_mobility_result_time_point(t, travellers); +} + +template +void apply_timed_mobility(const FP t, const FP num_moving, MobilityEdgeDirected& edge, + LocationNode& node_from, LocationNode& node_to, + mio::RandomNumberGenerator& rng) +{ + edge.apply_mobility(t, num_moving, node_from, node_to, rng); +} + +/** + * create a mobility-based simulation. + * After every second time step, for each edge a portion of the population corresponding to the coefficients of the edge + * changes from one node to the other. In the next timestep, the mobile population returns to their "home" node. + * Returns are adjusted based on the development in the target node. + * @param t0 start time of the simulation + * @param dt time step between mobility + * @param graph set up for mobility-based simulation + * @{ + */ +template +AsymmetricGraphSimulation, MobilityEdgeDirected>> +make_mobility_sim(FP t0, FP dt, const Graph, MobilityEdgeDirected>& graph) +{ + using GraphSim = + AsymmetricGraphSimulation, MobilityEdgeDirected>, FP, FP, + void (*)(FP, FP, mio::MobilityEdgeDirected&, mio::LocationNode&, + mio::LocationNode&, mio::RandomNumberGenerator&), + void (*)(FP, FP, mio::LocationNode&)>; + return GraphSim(t0, dt, graph, &advance_model, &apply_timed_mobility); +} + +template +AsymmetricGraphSimulation, MobilityEdgeDirected>> +make_mobility_sim(FP t0, FP dt, Graph, MobilityEdgeDirected>&& graph) +{ + using GraphSim = + AsymmetricGraphSimulation, MobilityEdgeDirected>, FP, FP, + void (*)(FP, FP, mio::MobilityEdgeDirected&, mio::LocationNode&, + mio::LocationNode&, mio::RandomNumberGenerator&), + void (*)(FP, FP, mio::LocationNode&)>; + return GraphSim(t0, dt, std::move(graph), &advance_model, &apply_timed_mobility); +} + +/** @} */ + +} // namespace mio + +#endif //METAPOPULATION_MOBILITY_STOCHASTIC_H diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 6f8a0b4c0c..e6f88d0f52 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -100,7 +100,7 @@ class SimulationNode m_last_state = m_simulation.get_result().get_last_value(); } -private: +protected: Sim m_simulation; Eigen::VectorX m_last_state; FP m_t0; @@ -193,7 +193,7 @@ class MobilityParameters //@} /** - * Get/Setthe mobility coefficients. + * Get/Set the mobility coefficients. * The coefficients represent the (time-dependent) percentage of people moving * from one node to another by age and infection compartment. * @{ @@ -656,7 +656,7 @@ void mio::MobilityEdge::apply_mobility(FP t, FP dt, SimulationNode& } /** - * edge functor for mobility-based simulation. + * node functor for mobility-based simulation. * @see SimulationNode::advance */ template @@ -677,7 +677,8 @@ void apply_mobility(FP t, FP dt, MobilityEdge& mobilityEdge, SimulationNode< } /** - * create a mobility-based simulation. + * @brief create a mobility-based simulation. + * * After every second time step, for each edge a portion of the population corresponding to the coefficients of the edge * changes from one node to the other. In the next timestep, the mobile population returns to their "home" node. * Returns are adjusted based on the development in the target node. @@ -713,9 +714,9 @@ make_mobility_sim(FP t0, FP dt, Graph, MobilityEdge> /** @} */ /** - * Create a graph simulation without mobility. + * @brief Create a graph simulation without mobility. * - * Note that we set the time step of the graph simulation to infinity since we do not require any exchange between the + * Note that we set the time step of the graph simulation to infinity since we do not require any exchange between the * nodes. Hence, in each node, the simulation runs until tmax when advancing the simulation without interruption. * * @param t0 Start time of the simulation. diff --git a/cpp/memilio/utils/index.h b/cpp/memilio/utils/index.h index cba56b1d57..dcc84d267d 100644 --- a/cpp/memilio/utils/index.h +++ b/cpp/memilio/utils/index.h @@ -20,15 +20,75 @@ #ifndef INDEX_H #define INDEX_H +#include "memilio/io/io.h" #include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/metaprogramming.h" #include "memilio/utils/type_safe.h" +#include +#include +#include +#include + namespace mio { template class Index; +namespace details +{ + +/// @brief Function definition that accepts a MultiIndex, used for the definition of IsMultiIndex. +template +void is_multi_index_impl(Index); + +} // namespace details + +/// @brief A MultiIndex is an Index with any number of categories. Does accept empty or single category indices. +template +concept IsMultiIndex = requires(Index i) { details::is_multi_index_impl(i); }; + +namespace details +{ + +/// @brief Obtain a tuple of single-category indices from a Index or MultiIndex. +template +std::tuple...> get_tuple(const Index& i) +{ + if constexpr (sizeof...(Tags) == 1) { + return std::tuple(i); + } + else { + return i.indices; + } +} + +/// @brief Obtain a tuple of one single-category index from an enum value. +template +std::tuple> get_tuple(Enum i) + requires std::is_enum::value +{ + return std::tuple(Index(i)); +} + +/// @brief Merge a series of enums or MultIndex%s into a tuple of single-category indices. +template + requires((std::is_enum_v || IsMultiIndex) && ...) +decltype(auto) concatenate_indices_impl(IndexArgs&&... args) +{ + return std::tuple_cat(details::get_tuple(args)...); +} + +/** + * @brief Function declaration that allows type conversion from a tuple of single-category indices to MultiIndex. + * Used together with concatenate_indices_impl, this allows combining categories of multiple args into a single MultiIndex. + */ +template +Index tuple_to_index(std::tuple...>); + +} // namespace details + /** * @brief An Index with a single template parameter is a typesafe wrapper for size_t * that is associated with a Tag. It is used to index into a CustomIndexArray @@ -62,9 +122,10 @@ class MEMILIO_ENABLE_EBO Index : public TypeSafe>::TypeSafe; - static size_t constexpr size = 1; + static constexpr size_t size = 1; + static constexpr bool has_duplicates = false; - static Index constexpr Zero() + static constexpr Index Zero() { return Index((size_t)0); } @@ -72,8 +133,8 @@ class MEMILIO_ENABLE_EBO Index : public TypeSafe::value, void>* = nullptr> - Index(Dummy val) + Index(CategoryTag val) + requires std::is_enum_v : TypeSafe>((size_t)val) { } @@ -119,20 +180,31 @@ template class Index { public: - static size_t constexpr size = sizeof...(CategoryTag); + static constexpr size_t size = sizeof...(CategoryTag); + static constexpr bool has_duplicates = has_duplicates_v; + /// @brief Construct an Index filled with zeroes. static Index constexpr Zero() { return Index(Index::Zero()...); } - // constructor from Indices + /// @brief Constructor from individual Indices. Index(Index const&... _indices) : indices{_indices...} { } + /// @brief Constructor from mixed Indices and MultiIndices. + template + requires(sizeof...(IndexArgs) > 1) + Index(IndexArgs&&... subindices) + : indices(details::concatenate_indices_impl(std::forward(subindices)...)) + { + } + private: + /// @brief Internal constructor from a tuple. Index(const std::tuple...>& _indices) : indices(_indices) { @@ -150,6 +222,25 @@ class Index return !(this == other); } + bool operator<(Index const& other) const + { + // use apply to unfold both tuples, then use a fold expression to evaluate a pairwise less + return std::apply( + [&other](auto&&... indices_) { + return std::apply( + [&](auto&&... other_indices_) { + return ((indices_ < other_indices_) && ...); + }, + other.indices); + }, + indices); + } + + bool operator<=(Index const& other) const + { + return (*this == other) || (*this < other); + } + /** * serialize this. * @see mio::serialize @@ -197,70 +288,75 @@ struct index_of_type, Index> { static constexpr std::size_t value = 0; }; -// retrieves the Index at the Ith position for a Index with more than one Tag -template 1), void>* = nullptr> +/// @brief Retrieves the Index (by reference) at the Ith position of a MultiIndex. +template constexpr typename std::tuple_element...>>::type& get(Index& i) noexcept { - return std::get(i.indices); -} - -// retrieves the Index at the Ith position for a Index with one Tag, equals identity function -template * = nullptr> -constexpr typename std::tuple_element...>>::type& -get(Index& i) noexcept -{ - static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); - return i; -} - -// retrieves the Index at the Ith position for a Index with more than one Tag const version -template 1), void>* = nullptr> -constexpr typename std::tuple_element...>>::type const& -get(Index const& i) noexcept -{ - return std::get(i.indices); + if constexpr (sizeof...(CategoryTags) == 1) { + static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); + return i; + } + else { + return std::get(i.indices); + } } -// retrieves the Index at the Ith position for a Index with one Tag, equals identity function const version -template * = nullptr> +/// @brief Retrieves the Index (by const reference) at the Ith position of a MultiIndex. +template constexpr typename std::tuple_element...>>::type const& get(Index const& i) noexcept { - static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); - return i; -} - -// retrieves the Index for the tag Tag of a Index with more than one Tag -template 1), void>* = nullptr> -constexpr Index& get(Index& i) noexcept -{ - return std::get>(i.indices); + if constexpr (sizeof...(CategoryTags) == 1) { + static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); + return i; + } + else { + return std::get(i.indices); + } } -// retrieves the Index for the tag Tag of a Index with one Tag, equals identity function -template * = nullptr> +/// @brief Retrieves the Index (by reference) by its Tag in a MultiIndex. Requires unique tags. +template constexpr Index& get(Index& i) noexcept { - using IndexTag = std::tuple_element_t<0, std::tuple>; - static_assert(std::is_same::value, "Tags must match for an Index with just one template parameter"); - return i; + if constexpr (sizeof...(CategoryTags) == 1) { + using IndexTag = std::tuple_element_t<0, std::tuple>; + static_assert(std::is_same::value, + "Tags must match for an Index with just one template parameter"); + return i; + } + else { + return std::get>(i.indices); + } } -// retrieves the Index for the tag Tag for a Index with more than one Tag const version -template 1), void>* = nullptr> +/// @brief Retrieves the Index (by const reference) by its Tag in a MultiIndex. Requires unique tags. +template constexpr Index const& get(Index const& i) noexcept { - return std::get>(i.indices); + if constexpr (sizeof...(CategoryTags) == 1) { + using IndexTag = std::tuple_element_t<0, std::tuple>; + static_assert(std::is_same::value, + "Tags must match for an Index with just one template parameter"); + return i; + } + else { + return std::get>(i.indices); + } } -// retrieves the Index for the tag Tag for a Index with one Tag, equals identity function const version -template * = nullptr> -constexpr Index const& get(Index const& i) noexcept +/** + * @brief Combine several Index%s into one MultiIndex. + * @param args Either enum or MultiIndex values. + * @return A MultiIndex with all categories and values of the given Indices concatonated. + */ +template +decltype(auto) concatenate_indices(IndexArgs&&... args) { - using IndexTag = std::tuple_element_t<0, std::tuple>; - static_assert(std::is_same::value, "Tags must match for an Index with just one template parameter"); - return i; + using MergedIndex = + decltype(details::tuple_to_index(details::concatenate_indices_impl(std::declval()...))); + return MergedIndex(std::forward(args)...); } namespace details @@ -310,12 +406,27 @@ inline Index extend_index_impl(const Index& i, const * @tparam SubIndex An Index that contains a subset of the categories from SuperIndex. * @tparam SuperIndex Any Index. * @return A (sub)index with the given categories and values from index. + * @{ */ template -SubIndex reduce_index(const SuperIndex& index) +decltype(auto) reduce_index(const SuperIndex& index) +{ + if constexpr (SubIndex::size == 1 && std::is_base_of_v, SubIndex>) { + // this case handles reducing from e.g. Index directly to AgeGroup + // the default case would instead reduce to Index, which may cause conversion errors + return details::reduce_index_impl(index, mio::Tag>{}); + } + else { + return details::reduce_index_impl(index, mio::Tag{}); + } +} +template + requires std::is_enum_v +Index reduce_index(const SuperIndex& index) { - return details::reduce_index_impl(index, mio::Tag{}); + return details::reduce_index_impl(index, mio::Tag>{}); } +/** @} */ /** * @brief Create a SuperIndex by copying values from SubIndex, filling new categories with fill_value. diff --git a/cpp/memilio/utils/stl_util.h b/cpp/memilio/utils/stl_util.h index d663843c09..f5e2a82b56 100644 --- a/cpp/memilio/utils/stl_util.h +++ b/cpp/memilio/utils/stl_util.h @@ -21,6 +21,7 @@ #define STL_UTIL_H #include "memilio/utils/metaprogramming.h" +#include "memilio/timer/auto_timer.h" #include #include @@ -66,6 +67,7 @@ inline std::ostream& set_ostream_format(std::ostream& out, size_t width, size_t template void insert_sorted_replace(std::vector& vec, T const& item, Pred pred) { + mio::timing::AutoTimer<"insert_sorted_replace"> timer; auto bounds = std::equal_range(begin(vec), end(vec), item, pred); auto lb = bounds.first; auto ub = bounds.second; diff --git a/cpp/memilio/utils/time_series.h b/cpp/memilio/utils/time_series.h index e43c3318b2..730515a01a 100644 --- a/cpp/memilio/utils/time_series.h +++ b/cpp/memilio/utils/time_series.h @@ -158,6 +158,26 @@ class TimeSeries return value_matrix; } + /** + * @brief constructs TimeSeries instance and initializes it with zeros + * @param num_time_points number of time steps + * @param num_elements number of compartiments * number of groups + * @param timepoints vector of time points + * @return + */ + static TimeSeries zero(Eigen::Index num_time_points, Eigen::Index num_elements, std::vector timepoints) + { + assert(num_time_points == static_cast(timepoints.size())); + TimeSeries value_matrix(num_elements); + value_matrix.m_data = Matrix::Zero(num_elements + 1, num_time_points); + value_matrix.m_num_time_points = num_time_points; + for (Eigen::Index i = 0; i < num_time_points; ++i) { + value_matrix.m_data(0, i) = timepoints[i]; + } + + return value_matrix; + } + /** copy assignment */ TimeSeries& operator=(const TimeSeries& other) { @@ -280,6 +300,21 @@ class TimeSeries return m_data(0, i); } + auto get_time_points() + { + return m_data(0, Eigen::all); + } + + auto get_index_of_time(FP time) + { + for (Eigen::Index i = 0; i < m_num_time_points; ++i) { + if (m_data(0, i) == time) { + return i; + } + } + return Eigen::Index(-1); + } + /** * time of time point at index num_time_points - 1 */ diff --git a/cpp/models/fmd/adoption_rates.h b/cpp/models/fmd/adoption_rates.h new file mode 100644 index 0000000000..7a53797d09 --- /dev/null +++ b/cpp/models/fmd/adoption_rates.h @@ -0,0 +1,71 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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. +*/ +#ifndef FMD_ADOPTION_RATES_H +#define FMD_ADOPTION_RATES_H + +#include "memilio/config.h" +#include "memilio/geography/geolocation.h" +#include "memilio/geography/rtree.h" +#include "memilio/geography/distance.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_asymmetric.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_builder.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/logging.h" +#include "memilio/timer/auto_timer.h" +#include "memilio/utils/parameter_distributions.h" +#include "memilio/utils/random_number_generator.h" +#include "smm/simulation.h" +#include "smm/parameters.h" +#include "fmd/infection_state.h" +#include + +namespace mio +{ +namespace fmd +{ +using Status = mio::Index; +using mio::regions::Region; + +using AR = mio::AdoptionRate; + +std::vector generic_adoption_rates() +{ + auto home = Region(0); + std::vector adoption_rates; + adoption_rates.push_back({InfectionState::S, + InfectionState::E, + home, + 0.2, + {{InfectionState::I, 0.3}, {InfectionState::INS, 0.3}, {InfectionState::ICS, 0.8}}}); + adoption_rates.push_back({InfectionState::E, InfectionState::I, home, 0.2, {}}); + adoption_rates.push_back({InfectionState::I, InfectionState::INS, home, 0.1, {}}); + adoption_rates.push_back({InfectionState::I, InfectionState::ICS, home, 0.1, {}}); + adoption_rates.push_back({InfectionState::ICS, InfectionState::D, home, 0.6, {}}); + adoption_rates.push_back({InfectionState::ICS, InfectionState::R, home, 0.4, {}}); + adoption_rates.push_back({InfectionState::INS, InfectionState::R, home, 0.5, {}}); + return adoption_rates; +} + +} // namespace fmd +} // namespace mio + +#endif // FMD_ADOPTION_RATES_Hs \ No newline at end of file diff --git a/cpp/models/fmd/infection_state.h b/cpp/models/fmd/infection_state.h new file mode 100644 index 0000000000..c60acf75f9 --- /dev/null +++ b/cpp/models/fmd/infection_state.h @@ -0,0 +1,60 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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. +*/ +#ifndef FMD_INFECTION_STATE_H +#define FMD_INFECTION_STATE_H + +#include "memilio/config.h" +#include "memilio/geography/geolocation.h" +#include "memilio/geography/rtree.h" +#include "memilio/geography/distance.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_asymmetric.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_builder.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/logging.h" +#include "memilio/timer/auto_timer.h" +#include "memilio/utils/parameter_distributions.h" +#include "memilio/utils/random_number_generator.h" +#include "smm/simulation.h" +#include "smm/parameters.h" + +namespace mio +{ +namespace fmd +{ + +enum class InfectionState +{ + S, + E, + I, + INS, + ICS, + R, + V, + D, + Count +}; + +} // namespace fmd +} // namespace mio + +#endif // FMD_INFECTION_STATE_H \ No newline at end of file diff --git a/cpp/models/fmd/model.h b/cpp/models/fmd/model.h new file mode 100644 index 0000000000..681e749b0d --- /dev/null +++ b/cpp/models/fmd/model.h @@ -0,0 +1,94 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Kilian Volmer +* +* 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. +*/ +#ifndef FMD_MODEL_H +#define FMD_MODEL_H + +#include "memilio/config.h" +#include "memilio/geography/rtree.h" +#include "memilio/geography/distance.h" +#include "memilio/mobility/metapopulation_mobility_asymmetric.h" +#include "memilio/mobility/graph.h" +#include "memilio/mobility/graph_builder.h" +#include "smm/simulation.h" +#include "smm/parameters.h" +#include "fmd/infection_state.h" +#include "fmd/adoption_rates.h" + +namespace mio +{ +namespace fmd +{ + +using Model = mio::smm::Model; + +using Builder = mio::GraphBuilder< + mio::LocationNode>, + mio::MobilityEdgeDirected>; +using Graph = mio::Graph< + mio::LocationNode>, + mio::MobilityEdgeDirected>; + +/** + * @brief Create a model object with given number of susceptible individuals and adoption rates. + * + * @param num_sus Number of susceptible individuals. + * @param adoption_rates Vector of adoption rates. + * @return Model The created model. + */ +Model create_model(size_t num_sus, const std::vector& adoption_rates = generic_adoption_rates()) +{ + auto home = mio::regions::Region(0); + Model curr_model(Status{InfectionState::Count}, mio::regions::Region(1)); + curr_model.populations[{home, InfectionState::S}] = num_sus; + curr_model.populations[{home, InfectionState::E}] = 0; + curr_model.populations[{home, InfectionState::I}] = 0; + curr_model.populations[{home, InfectionState::INS}] = 0; + curr_model.populations[{home, InfectionState::ICS}] = 0; + curr_model.populations[{home, InfectionState::R}] = 0; + curr_model.populations[{home, InfectionState::V}] = 0; + curr_model.populations[{home, InfectionState::D}] = 0; + curr_model.parameters.get>() = adoption_rates; + return curr_model; +} + +/** + * @brief Generate and set regional neighbors for all nodes in the graph based on given distances. + * + * The function constructs an RTree from the node locations and queries it for neighbors within the specified distances. + * + * @param graph The graph for which to generate neighbors. + * @param distances The distances to consider for neighbor generation. + */ +void generate_neighbours(Graph& graph, std::vector distances) +{ + auto nodes = graph.nodes() | std::views::transform([](const auto& node) { + return &node.property; + }); + auto tree = mio::geo::RTree(nodes.begin(), nodes.end()); + + for (auto& node : graph.nodes()) { + node.property.set_regional_neighbors(tree.in_range_indices_query(node.property.get_location(), distances)); + } +} + +} // namespace fmd +} // namespace mio + +#endif // FMD_MODEL_H \ No newline at end of file diff --git a/cpp/models/hybrid/conversion_functions.cpp b/cpp/models/hybrid/conversion_functions.cpp index 8fa4b047b4..7fa8d4831a 100644 --- a/cpp/models/hybrid/conversion_functions.cpp +++ b/cpp/models/hybrid/conversion_functions.cpp @@ -36,7 +36,7 @@ namespace hybrid { template <> void convert_model(const dabm::Simulation>& current_model, - smm::Simulation& target_model) + smm::Simulation& target_model) { auto& current_result = current_model.get_result(); auto& target_result = target_model.get_result(); @@ -58,7 +58,7 @@ void convert_model(const dabm::Simulation -void convert_model(const smm::Simulation& current_model, +void convert_model(const smm::Simulation& current_model, dabm::Simulation>& target_model) { auto& current_result = current_model.get_result(); diff --git a/cpp/models/hybrid/conversion_functions.h b/cpp/models/hybrid/conversion_functions.h index 9d0e316543..2d15f6d638 100644 --- a/cpp/models/hybrid/conversion_functions.h +++ b/cpp/models/hybrid/conversion_functions.h @@ -37,11 +37,12 @@ namespace hybrid template <> void convert_model(const dabm::Simulation>& current_model, - smm::Simulation& target_model); + smm::Simulation>& target_model); template <> -void convert_model(const smm::Simulation& current_model, - dabm::Simulation>& target_model); +void convert_model( + const smm::Simulation>& current_model, + dabm::Simulation>& target_model); template <> void convert_model(const dabm::Simulation>& current_model, diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 56d9b56251..056878cfb7 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: René Schmieding, Julia Bicker +* Authors: René Schmieding, Julia Bicker, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -21,32 +21,46 @@ #ifndef MIO_SMM_MODEL_H #define MIO_SMM_MODEL_H -#include "memilio/config.h" +#include "memilio/utils/index.h" +#include "memilio/utils/index_range.h" #include "smm/parameters.h" #include "memilio/compartments/compartmental_model.h" #include "memilio/epidemiology/populations.h" #include "memilio/geography/regions.h" +#include namespace mio { namespace smm { +/// @brief The Index type used to define the SMM population. +template +using PopulationIndex = decltype(concatenate_indices(std::declval(), std::declval())); + /** * @brief Stochastic Metapopulation Model. - * @tparam regions Number of regions. - * @tparam Status An infection state enum. + * The stratification of the population of this model is split between "Status" and "Region". This split is mostly + * arbitrary, with the important distinction, that for second order rates the reference population + * (i.e., the N in S' = S * I / N) is calculated by accumulating subpopulations only over the Status. + * @tparam Comp An enum representing the infection states. Must also be contained in Status + * @tparam Status A MultiIndex allowing to further stratify infection state adoptions. + * @tparam Region A MultiIndex for "spatially" separate subpopulations, default is @ref mio::regions::Region. */ -template -class Model : public mio::CompartmentalModel, - ParametersBase> +template +class Model : public mio::CompartmentalModel>, + ParametersBase> { - using Base = mio::CompartmentalModel, - ParametersBase>; + using Base = mio::CompartmentalModel>, + ParametersBase>; + static_assert(!Base::Populations::Index::has_duplicates, "Do not use the same Index tag multiple times!"); public: - Model() - : Base(typename Base::Populations({static_cast(regions), Status::Count}, 0.0), + using Status = StatusT; + using Region = RegionT; + + Model(Status status_dimensions, Region region_dimensions) + : Base(typename Base::Populations(concatenate_indices(region_dimensions, status_dimensions)), typename Base::ParameterSet()) { } @@ -57,7 +71,7 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorX& x) const + FP evaluate(const AdoptionRate& rate, const Eigen::VectorXd& x) const { const auto& pop = this->populations; const auto source = pop.get_flat_index({rate.region, rate.from}); @@ -67,8 +81,8 @@ class Model : public mio::CompartmentalModel(Status::Count); ++s) { - N += x[pop.get_flat_index({rate.region, Status(s)})]; + for (auto status : make_index_range(reduce_index(this->populations.size()))) { + N += x[pop.get_flat_index({rate.region, status})]; } // accumulate influences FP influences = 0.0; @@ -86,7 +100,7 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorX& x) const + FP evaluate(const TransitionRate& rate, const Eigen::VectorXd& x) const { const auto source = this->populations.get_flat_index({rate.from, rate.status}); return rate.factor * x[source]; diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index d0a6011f6e..652d8891d8 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -33,11 +33,13 @@ namespace smm /** * @brief A vector of AdoptionRate%s, see mio::AdoptionRate - * @tparam Status An infection state enum. + * @tparam FP A floating point type, e.g., double. + * @tparam Status A MultiIndex, containing the infection state enum. + * @tparam Region A MultiIndex for spatial stratification. */ -template +template struct AdoptionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "AdoptionRates"; @@ -46,27 +48,35 @@ struct AdoptionRates { /** * @brief Struct defining a possible regional transition in a Model based on Poisson Processes. - * @tparam Status An infection state enum. + * @tparam FP A floating point type, e.g., double. + * @tparam Status A MultiIndex, containing the infection state enum. + * @tparam Region A MultiIndex for spatial stratification. */ -template +template struct TransitionRate { Status status; // i - mio::regions::Region from; // k - mio::regions::Region to; // l + Region from; // k + Region to; // l FP factor; // lambda_i^{kl} }; -template +/** + * @brief A vector of TransitionRate%s, see mio::TransitionRate + * @tparam FP A floating point type, e.g., double. + * @tparam Status A MultiIndex, containing the infection state enum. + * @tparam Region A MultiIndex for spatial stratification. + */ +template struct TransitionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "TransitionRates"; } }; -template -using ParametersBase = mio::ParameterSet, TransitionRates>; +template +using ParametersBase = mio::ParameterSet, TransitionRates>; } // namespace smm diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 49d4dcbb8f..fe0a83655f 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -21,10 +21,9 @@ #ifndef MIO_SMM_SIMULATION_H #define MIO_SMM_SIMULATION_H -#include "memilio/config.h" #include "smm/model.h" #include "smm/parameters.h" -#include "memilio/compartments/simulation.h" +#include "memilio/utils/time_series.h" namespace mio { @@ -37,12 +36,11 @@ namespace smm * @tparam regions The number of regions. * @tparam Status An infection state enum. */ -template +template class Simulation { public: -public: - using Model = smm::Model; + using Model = smm::Model; /** * @brief Set up the simulation for a Stochastic Metapopulation Model. @@ -61,18 +59,44 @@ class Simulation { assert(dt > 0); assert(m_waiting_times.size() > 0); - assert(std::all_of(adoption_rates().begin(), adoption_rates().end(), [](auto&& r) { - return static_cast(r.region) < regions; - })); - assert(std::all_of(transition_rates().begin(), transition_rates().end(), [](auto&& r) { - return static_cast(r.from) < regions && static_cast(r.to) < regions; - })); + assert(std::all_of(adoption_rates().begin(), adoption_rates().end(), + [regions = reduce_index(model.populations.size())](auto&& r) { + return r.region < regions; + })); + assert(std::all_of(transition_rates().begin(), transition_rates().end(), + [regions = reduce_index(model.populations.size())](auto&& r) { + return r.from < regions && r.to < regions; + })); // initialize (internal) next event times by random values for (size_t i = 0; i < m_tp_next_event.size(); i++) { m_tp_next_event[i] += mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); } } + Simulation(const Simulation& other) + : m_dt(other.m_dt) + , m_model(std::make_unique(*other.m_model)) + , m_result(other.m_result) + , m_internal_time(other.m_internal_time) + , m_tp_next_event(other.m_tp_next_event) + , m_waiting_times(other.m_waiting_times) + , m_current_rates(other.m_current_rates) + { + } + + Simulation& operator=(const Simulation& other) + { + if (this != &other) { + m_model = std::make_unique(*other.m_model); + m_result = other.m_result; + m_internal_time = other.m_internal_time; + m_tp_next_event = other.m_tp_next_event; + m_waiting_times = other.m_waiting_times; + m_current_rates = other.m_current_rates; + } + return *this; + } + /** * @brief Advance simulation to tmax. * This function performs a Gillespie algorithm. @@ -164,17 +188,17 @@ class Simulation /** * @brief Returns the model's transition rates. */ - inline constexpr const typename smm::TransitionRates::Type& transition_rates() + inline constexpr const typename smm::TransitionRates::Type& transition_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get>(); } /** * @brief Returns the model's adoption rates. */ - inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() + inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get>(); } /** @@ -217,7 +241,7 @@ class Simulation std::vector m_current_rates; ///< Current values of both types of rates i.e. adoption and transition rates. }; -} //namespace smm +} // namespace smm } // namespace mio #endif diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index b671cfa2d0..0e906dc13f 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: Julia Bicker +* Authors: Julia Bicker, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -18,18 +18,17 @@ * limitations under the License. */ -#include -#include -#include -#include -#include -#include "memilio/utils/compiler_diagnostics.h" +#include "abm_helpers.h" +#include "memilio/epidemiology/adoption_rate.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/geography/regions.h" #include "smm/model.h" #include "smm/parameters.h" #include "smm/simulation.h" -#include "abm_helpers.h" -#include "memilio/epidemiology/adoption_rate.h" -#include + +#include + +#include #include #include @@ -45,6 +44,14 @@ enum class InfectionState }; +enum class Infections +{ + S, + I, + R, + Count +}; + TEST(TestSMM, evaluateAdoptionRate) { //Test whether the adoption rates are evaluated correctly. @@ -53,9 +60,9 @@ TEST(TestSMM, evaluateAdoptionRate) // with N(from, region) the population in Region "region" having infection state "from" //Second-order adoption rates are given by: // rate.factor * N(rate.from, rate.region)/total_pop * sum (over all rate.influences) influence.factor * N(influence.status, rate.region) - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(1)); //Set adoption rates std::vector> adoption_rates; @@ -80,13 +87,52 @@ TEST(TestSMM, evaluateAdoptionRate) EXPECT_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 2.); } +TEST(TestSMM, evaluateinterregionalAdoptionRate) +{ + using Age = mio::AgeGroup; + using mio::regions::Region; + using Status = mio::Index; + using Model = mio::smm::Model>; + + std::vector>> adoption_rates; + //Second-order adoption + adoption_rates.push_back({{Infections::S, Age(0), Region(0)}, + {Infections::I, Age(0), Region(0)}, + {}, + 0.1, + {{{Infections::I, Age(1), Region(0)}, 0.1}, {{Infections::I, Age(1), Region(1)}, 0.2}}}); + //First-order adoption + adoption_rates.push_back({{Infections::I, Age(1), Region(0)}, {Infections::R, Age(1), Region(0)}, {}, 0.2, {}}); + Model model({Infections::Count, Age(2), Region(2)}, {}); + + model.populations[{Infections::S, Age(0), Region(0)}] = 50; + model.populations[{Infections::I, Age(0), Region(0)}] = 10; + model.populations[{Infections::R, Age(0), Region(0)}] = 0; + + model.populations[{Infections::S, Age(1), Region(0)}] = 100; + model.populations[{Infections::I, Age(1), Region(0)}] = 20; + model.populations[{Infections::R, Age(1), Region(0)}] = 0; + + model.populations[{Infections::S, Age(0), Region(1)}] = 40; + model.populations[{Infections::I, Age(0), Region(1)}] = 80; + model.populations[{Infections::R, Age(0), Region(1)}] = 0; + + model.populations[{Infections::S, Age(1), Region(1)}] = 80; + model.populations[{Infections::I, Age(1), Region(1)}] = 16; + model.populations[{Infections::R, Age(1), Region(1)}] = 0; + + EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()), + 0.1 * 50. * (0.1 * 20. * 1. + 0.2 * 16 * 1.) / (60 + 120 + 120 + 96)); + EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 4.); +} + TEST(TestSMM, evaluateTransitionRate) { //Same test as 'evaluateAdoptionRate' only for spatial transition rates. //Transition rates are given by: rate.factor * N(rate.status, rate.from) - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); //Initialize model populations model.populations[{mio::regions::Region(0), InfectionState::S}] = 50; model.populations[{mio::regions::Region(0), InfectionState::E}] = 10; @@ -114,9 +160,9 @@ TEST(TestSMMSimulation, advance) { //Test whether Gillespie algorithm calculates events in the correct order using testing::Return; - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); //Initialize model populations model.populations[{mio::regions::Region(0), InfectionState::S}] = 1; model.populations[{mio::regions::Region(0), InfectionState::E}] = 0; @@ -151,8 +197,8 @@ TEST(TestSMMSimulation, advance) //Spatial transition transition_rates.push_back({InfectionState::R, mio::regions::Region(1), mio::regions::Region(0), 0.01}); - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; //Mock exponential distribution to control the normalized waiting times that are drawn ScopedMockDistribution>>> @@ -195,9 +241,9 @@ TEST(TestSMMSimulation, stopsAtTmax) { //Test whether simulation stops at tmax and whether system state at tmax is saved if there are no adoptions/transitions using testing::Return; - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); //Set adoption and spatial transition rates std::vector> adoption_rates; @@ -211,8 +257,8 @@ TEST(TestSMMSimulation, stopsAtTmax) transition_rates.push_back({InfectionState::R, mio::regions::Region(1), mio::regions::Region(0), 0.01}); - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; //As populations are not set they have value 0 i.e. no events will happen //Simulate for 30 days @@ -228,7 +274,7 @@ TEST(TestSMMSimulation, convergence) //Test whether the mean number of transitions in one unit-time step corresponds to the expected number of transitions // given by rate * pop up to some tolerance using testing::Return; - using Model = mio::smm::Model; + using Model = mio::smm::Model; double pop = 1000; size_t num_runs = 100; @@ -243,7 +289,7 @@ TEST(TestSMMSimulation, convergence) //Do 100 unit-time step simulations with 1000 agents for (size_t n = 0; n < num_runs; ++n) { - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); model.populations[{mio::regions::Region(0), InfectionState::S}] = pop; model.populations[{mio::regions::Region(0), InfectionState::E}] = 0; @@ -252,13 +298,14 @@ TEST(TestSMMSimulation, convergence) model.populations[{mio::regions::Region(0), InfectionState::R}] = 0; model.populations[{mio::regions::Region(0), InfectionState::D}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::S}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::E}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::C}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::I}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::R}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::D}] = 0; - model.parameters.get>() = transition_rates; + model.populations[{mio::regions::Region(1), InfectionState::S}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::E}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::C}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::I}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::R}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::D}] = 0; + model.parameters.get>() = + transition_rates; auto sim = mio::smm::Simulation(model, 0.0, 1.0); sim.advance(1.); diff --git a/cpp/tests/test_temporal_hybrid_model.cpp b/cpp/tests/test_temporal_hybrid_model.cpp index 0ae31660ec..5f0ff64891 100644 --- a/cpp/tests/test_temporal_hybrid_model.cpp +++ b/cpp/tests/test_temporal_hybrid_model.cpp @@ -151,7 +151,7 @@ TEST(TestTemporalHybrid, test_advance) TEST(TestTemporalHybrid, test_conversion_dabm_smm) { using Model1 = mio::dabm::Model>; - using Model2 = mio::smm::Model; + using Model2 = mio::smm::Model; //Initialize agents for dabm SingleWell::Agent a1{Eigen::Vector2d{-0.5, 0}, @@ -162,13 +162,14 @@ TEST(TestTemporalHybrid, test_conversion_dabm_smm) mio::osecir::InfectionState::InfectedSymptoms}; Model1 model1({a1, a2, a3}, {}); - Model2 model2; - model2.parameters.get>().push_back( - {mio::osecir::InfectionState::Susceptible, - mio::osecir::InfectionState::Exposed, - mio::regions::Region(0), - 0.1, - {{mio::osecir::InfectionState::InfectedNoSymptoms, 1}, {mio::osecir::InfectionState::InfectedSymptoms, 0.5}}}); + Model2 model2(mio::osecir::InfectionState::Count, mio::regions::Region(1)); + model2.parameters.get>() + .push_back({mio::osecir::InfectionState::Susceptible, + mio::osecir::InfectionState::Exposed, + mio::regions::Region(0), + 0.1, + {{mio::osecir::InfectionState::InfectedNoSymptoms, 1}, + {mio::osecir::InfectionState::InfectedSymptoms, 0.5}}}); //Parameters for simulation double t0 = 0; diff --git a/cpp/thirdparty/csv.h b/cpp/thirdparty/csv.h new file mode 100644 index 0000000000..a4741f03ff --- /dev/null +++ b/cpp/thirdparty/csv.h @@ -0,0 +1,1371 @@ +// Copyright: (2012-2015) Ben Strasser +// License: BSD-3 +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// 3. Neither the name of the copyright holder nor the names of its contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#ifndef CSV_H +#define CSV_H + +#include +#include +#include +#include +#include +#include +#include +#ifndef CSV_IO_NO_THREAD +#include +#include +#include +#endif +#include +#include +#include +#include +#include + +namespace io +{ +//////////////////////////////////////////////////////////////////////////// +// LineReader // +//////////////////////////////////////////////////////////////////////////// + +namespace error +{ +struct base : std::exception { + virtual void format_error_message() const = 0; + + const char* what() const noexcept override + { + format_error_message(); + return error_message_buffer; + } + + mutable char error_message_buffer[2048]; +}; + +// this only affects the file name in the error message +const int max_file_name_length = 1024; + +struct with_file_name { + with_file_name() + { + std::memset(file_name, 0, sizeof(file_name)); + } + + void set_file_name(const char* file_name) + { + if (file_name != nullptr) { + // This call to strncpy has parenthesis around it + // to silence the GCC -Wstringop-truncation warning + (strncpy(this->file_name, file_name, sizeof(this->file_name))); + this->file_name[sizeof(this->file_name) - 1] = '\0'; + } + else { + this->file_name[0] = '\0'; + } + } + + char file_name[max_file_name_length + 1]; +}; + +struct with_file_line { + with_file_line() + { + file_line = -1; + } + + void set_file_line(int file_line) + { + this->file_line = file_line; + } + + int file_line; +}; + +struct with_errno { + with_errno() + { + errno_value = 0; + } + + void set_errno(int errno_value) + { + this->errno_value = errno_value; + } + + int errno_value; +}; + +struct can_not_open_file : base, with_file_name, with_errno { + void format_error_message() const override + { + if (errno_value != 0) + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Can not open file \"%s\" because \"%s\".", file_name, std::strerror(errno_value)); + else + std::snprintf(error_message_buffer, sizeof(error_message_buffer), "Can not open file \"%s\".", file_name); + } +}; + +struct line_length_limit_exceeded : base, with_file_name, with_file_line { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Line number %d in file \"%s\" exceeds the maximum length of 2^24-1.", file_line, file_name); + } +}; +} // namespace error + +class ByteSourceBase +{ +public: + virtual int read(char* buffer, int size) = 0; + virtual ~ByteSourceBase() + { + } +}; + +namespace detail +{ + +class OwningStdIOByteSourceBase : public ByteSourceBase +{ +public: + explicit OwningStdIOByteSourceBase(FILE* file) + : file(file) + { + // Tell the std library that we want to do the buffering ourself. + std::setvbuf(file, 0, _IONBF, 0); + } + + int read(char* buffer, int size) + { + return std::fread(buffer, 1, size, file); + } + + ~OwningStdIOByteSourceBase() + { + std::fclose(file); + } + +private: + FILE* file; +}; + +class NonOwningIStreamByteSource : public ByteSourceBase +{ +public: + explicit NonOwningIStreamByteSource(std::istream& in) + : in(in) + { + } + + int read(char* buffer, int size) + { + in.read(buffer, size); + return in.gcount(); + } + + ~NonOwningIStreamByteSource() + { + } + +private: + std::istream& in; +}; + +class NonOwningStringByteSource : public ByteSourceBase +{ +public: + NonOwningStringByteSource(const char* str, long long size) + : str(str) + , remaining_byte_count(size) + { + } + + int read(char* buffer, int desired_byte_count) + { + int to_copy_byte_count = desired_byte_count; + if (remaining_byte_count < to_copy_byte_count) + to_copy_byte_count = remaining_byte_count; + std::memcpy(buffer, str, to_copy_byte_count); + remaining_byte_count -= to_copy_byte_count; + str += to_copy_byte_count; + return to_copy_byte_count; + } + + ~NonOwningStringByteSource() + { + } + +private: + const char* str; + long long remaining_byte_count; +}; + +#ifndef CSV_IO_NO_THREAD +class AsynchronousReader +{ +public: + void init(std::unique_ptr arg_byte_source) + { + std::unique_lock guard(lock); + byte_source = std::move(arg_byte_source); + desired_byte_count = -1; + termination_requested = false; + worker = std::thread([&] { + std::unique_lock guard(lock); + try { + for (;;) { + read_requested_condition.wait(guard, [&] { + return desired_byte_count != -1 || termination_requested; + }); + if (termination_requested) + return; + + read_byte_count = byte_source->read(buffer, desired_byte_count); + desired_byte_count = -1; + if (read_byte_count == 0) + break; + read_finished_condition.notify_one(); + } + } + catch (...) { + read_error = std::current_exception(); + } + read_finished_condition.notify_one(); + }); + } + + bool is_valid() const + { + return byte_source != nullptr; + } + + void start_read(char* arg_buffer, int arg_desired_byte_count) + { + std::unique_lock guard(lock); + buffer = arg_buffer; + desired_byte_count = arg_desired_byte_count; + read_byte_count = -1; + read_requested_condition.notify_one(); + } + + int finish_read() + { + std::unique_lock guard(lock); + read_finished_condition.wait(guard, [&] { + return read_byte_count != -1 || read_error; + }); + if (read_error) + std::rethrow_exception(read_error); + else + return read_byte_count; + } + + ~AsynchronousReader() + { + if (byte_source != nullptr) { + { + std::unique_lock guard(lock); + termination_requested = true; + } + read_requested_condition.notify_one(); + worker.join(); + } + } + +private: + std::unique_ptr byte_source; + + std::thread worker; + + bool termination_requested; + std::exception_ptr read_error; + char* buffer; + int desired_byte_count; + int read_byte_count; + + std::mutex lock; + std::condition_variable read_finished_condition; + std::condition_variable read_requested_condition; +}; +#endif + +class SynchronousReader +{ +public: + void init(std::unique_ptr arg_byte_source) + { + byte_source = std::move(arg_byte_source); + } + + bool is_valid() const + { + return byte_source != nullptr; + } + + void start_read(char* arg_buffer, int arg_desired_byte_count) + { + buffer = arg_buffer; + desired_byte_count = arg_desired_byte_count; + } + + int finish_read() + { + return byte_source->read(buffer, desired_byte_count); + } + +private: + std::unique_ptr byte_source; + char* buffer; + int desired_byte_count; +}; +} // namespace detail + +class LineReader +{ +private: + static const int block_len = 1 << 20; + std::unique_ptr buffer; // must be constructed before (and thus + // destructed after) the reader! +#ifdef CSV_IO_NO_THREAD + detail::SynchronousReader reader; +#else + detail::AsynchronousReader reader; +#endif + int data_begin; + int data_end; + + char file_name[error::max_file_name_length + 1]; + unsigned file_line; + + static std::unique_ptr open_file(const char* file_name) + { + // We open the file in binary mode as it makes no difference under *nix + // and under Windows we handle \r\n newlines ourself. + FILE* file = std::fopen(file_name, "rb"); + if (file == 0) { + int x = errno; // store errno as soon as possible, doing it after + // constructor call can fail. + error::can_not_open_file err; + err.set_errno(x); + err.set_file_name(file_name); + throw err; + } + return std::unique_ptr(new detail::OwningStdIOByteSourceBase(file)); + } + + void init(std::unique_ptr byte_source) + { + file_line = 0; + + buffer = std::unique_ptr(new char[3 * block_len]); + data_begin = 0; + data_end = byte_source->read(buffer.get(), 2 * block_len); + + // Ignore UTF-8 BOM + if (data_end >= 3 && buffer[0] == '\xEF' && buffer[1] == '\xBB' && buffer[2] == '\xBF') + data_begin = 3; + + if (data_end == 2 * block_len) { + reader.init(std::move(byte_source)); + reader.start_read(buffer.get() + 2 * block_len, block_len); + } + } + +public: + LineReader() = delete; + LineReader(const LineReader&) = delete; + LineReader& operator=(const LineReader&) = delete; + + explicit LineReader(const char* file_name) + { + set_file_name(file_name); + init(open_file(file_name)); + } + + explicit LineReader(const std::string& file_name) + { + set_file_name(file_name.c_str()); + init(open_file(file_name.c_str())); + } + + LineReader(const char* file_name, std::unique_ptr byte_source) + { + set_file_name(file_name); + init(std::move(byte_source)); + } + + LineReader(const std::string& file_name, std::unique_ptr byte_source) + { + set_file_name(file_name.c_str()); + init(std::move(byte_source)); + } + + LineReader(const char* file_name, const char* data_begin, const char* data_end) + { + set_file_name(file_name); + init(std::unique_ptr(new detail::NonOwningStringByteSource(data_begin, data_end - data_begin))); + } + + LineReader(const std::string& file_name, const char* data_begin, const char* data_end) + { + set_file_name(file_name.c_str()); + init(std::unique_ptr(new detail::NonOwningStringByteSource(data_begin, data_end - data_begin))); + } + + LineReader(const char* file_name, FILE* file) + { + set_file_name(file_name); + init(std::unique_ptr(new detail::OwningStdIOByteSourceBase(file))); + } + + LineReader(const std::string& file_name, FILE* file) + { + set_file_name(file_name.c_str()); + init(std::unique_ptr(new detail::OwningStdIOByteSourceBase(file))); + } + + LineReader(const char* file_name, std::istream& in) + { + set_file_name(file_name); + init(std::unique_ptr(new detail::NonOwningIStreamByteSource(in))); + } + + LineReader(const std::string& file_name, std::istream& in) + { + set_file_name(file_name.c_str()); + init(std::unique_ptr(new detail::NonOwningIStreamByteSource(in))); + } + + void set_file_name(const std::string& file_name) + { + set_file_name(file_name.c_str()); + } + + void set_file_name(const char* file_name) + { + if (file_name != nullptr) { + strncpy(this->file_name, file_name, sizeof(this->file_name) - 1); + this->file_name[sizeof(this->file_name) - 1] = '\0'; + } + else { + this->file_name[0] = '\0'; + } + } + + const char* get_truncated_file_name() const + { + return file_name; + } + + void set_file_line(unsigned file_line) + { + this->file_line = file_line; + } + + unsigned get_file_line() const + { + return file_line; + } + + char* next_line() + { + if (data_begin == data_end) + return nullptr; + + ++file_line; + + assert(data_begin < data_end); + assert(data_end <= block_len * 2); + + if (data_begin >= block_len) { + std::memcpy(buffer.get(), buffer.get() + block_len, block_len); + data_begin -= block_len; + data_end -= block_len; + if (reader.is_valid()) { + data_end += reader.finish_read(); + std::memcpy(buffer.get() + block_len, buffer.get() + 2 * block_len, block_len); + reader.start_read(buffer.get() + 2 * block_len, block_len); + } + } + + int line_end = data_begin; + while (line_end != data_end && buffer[line_end] != '\n') { + ++line_end; + } + + if (line_end - data_begin + 1 > block_len) { + error::line_length_limit_exceeded err; + err.set_file_name(file_name); + err.set_file_line(file_line); + throw err; + } + + if (line_end != data_end && buffer[line_end] == '\n') { + buffer[line_end] = '\0'; + } + else { + // some files are missing the newline at the end of the + // last line + ++data_end; + buffer[line_end] = '\0'; + } + + // handle windows \r\n-line breaks + if (line_end != data_begin && buffer[line_end - 1] == '\r') + buffer[line_end - 1] = '\0'; + + char* ret = buffer.get() + data_begin; + data_begin = line_end + 1; + return ret; + } +}; + +//////////////////////////////////////////////////////////////////////////// +// CSV // +//////////////////////////////////////////////////////////////////////////// + +namespace error +{ +const int max_column_name_length = 63; +struct with_column_name { + with_column_name() + { + std::memset(column_name, 0, max_column_name_length + 1); + } + + void set_column_name(const char* column_name) + { + if (column_name != nullptr) { + std::strncpy(this->column_name, column_name, max_column_name_length); + this->column_name[max_column_name_length] = '\0'; + } + else { + this->column_name[0] = '\0'; + } + } + + char column_name[max_column_name_length + 1]; +}; + +const int max_column_content_length = 63; + +struct with_column_content { + with_column_content() + { + std::memset(column_content, 0, max_column_content_length + 1); + } + + void set_column_content(const char* column_content) + { + if (column_content != nullptr) { + std::strncpy(this->column_content, column_content, max_column_content_length); + this->column_content[max_column_content_length] = '\0'; + } + else { + this->column_content[0] = '\0'; + } + } + + char column_content[max_column_content_length + 1]; +}; + +struct extra_column_in_header : base, with_file_name, with_column_name { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(Extra column "%s" in header of file "%s".)", column_name, file_name); + } +}; + +struct missing_column_in_header : base, with_file_name, with_column_name { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(Missing column "%s" in header of file "%s".)", column_name, file_name); + } +}; + +struct duplicated_column_in_header : base, with_file_name, with_column_name { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(Duplicated column "%s" in header of file "%s".)", column_name, file_name); + } +}; + +struct header_missing : base, with_file_name { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), "Header missing in file \"%s\".", file_name); + } +}; + +struct too_few_columns : base, with_file_name, with_file_line { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), "Too few columns in line %d in file \"%s\".", + file_line, file_name); + } +}; + +struct too_many_columns : base, with_file_name, with_file_line { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), "Too many columns in line %d in file \"%s\".", + file_line, file_name); + } +}; + +struct escaped_string_not_closed : base, with_file_name, with_file_line { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Escaped string was not closed in line %d in file \"%s\".", file_line, file_name); + } +}; + +struct integer_must_be_positive : base, with_file_name, with_file_line, with_column_name, with_column_content { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" must be positive or 0 in column "%s" in file "%s" in line "%d".)", + column_content, column_name, file_name, file_line); + } +}; + +struct no_digit : base, with_file_name, with_file_line, with_column_name, with_column_content { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" contains an invalid digit in column "%s" in file "%s" in line "%d".)", + column_content, column_name, file_name, file_line); + } +}; + +struct integer_overflow : base, with_file_name, with_file_line, with_column_name, with_column_content { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" overflows in column "%s" in file "%s" in line "%d".)", column_content, + column_name, file_name, file_line); + } +}; + +struct integer_underflow : base, with_file_name, with_file_line, with_column_name, with_column_content { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" underflows in column "%s" in file "%s" in line "%d".)", column_content, + column_name, file_name, file_line); + } +}; + +struct invalid_single_character : base, with_file_name, with_file_line, with_column_name, with_column_content { + void format_error_message() const override + { + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The content "%s" of column "%s" in file "%s" in line "%d" is not a single character.)", + column_content, column_name, file_name, file_line); + } +}; +} // namespace error + +using ignore_column = unsigned int; +static const ignore_column ignore_no_column = 0; +static const ignore_column ignore_extra_column = 1; +static const ignore_column ignore_missing_column = 2; + +template +struct trim_chars { +private: + constexpr static bool is_trim_char(char) + { + return false; + } + + template + constexpr static bool is_trim_char(char c, char trim_char, OtherTrimChars... other_trim_chars) + { + return c == trim_char || is_trim_char(c, other_trim_chars...); + } + +public: + static void trim(char*& str_begin, char*& str_end) + { + while (str_begin != str_end && is_trim_char(*str_begin, trim_char_list...)) + ++str_begin; + while (str_begin != str_end && is_trim_char(*(str_end - 1), trim_char_list...)) + --str_end; + *str_end = '\0'; + } +}; + +struct no_comment { + static bool is_comment(const char*) + { + return false; + } +}; + +template +struct single_line_comment { +private: + constexpr static bool is_comment_start_char(char) + { + return false; + } + + template + constexpr static bool is_comment_start_char(char c, char comment_start_char, + OtherCommentStartChars... other_comment_start_chars) + { + return c == comment_start_char || is_comment_start_char(c, other_comment_start_chars...); + } + +public: + static bool is_comment(const char* line) + { + return is_comment_start_char(*line, comment_start_char_list...); + } +}; + +struct empty_line_comment { + static bool is_comment(const char* line) + { + if (*line == '\0') + return true; + while (*line == ' ' || *line == '\t') { + ++line; + if (*line == 0) + return true; + } + return false; + } +}; + +template +struct single_and_empty_line_comment { + static bool is_comment(const char* line) + { + return single_line_comment::is_comment(line) || + empty_line_comment::is_comment(line); + } +}; + +template +struct no_quote_escape { + static const char* find_next_column_end(const char* col_begin) + { + while (*col_begin != sep && *col_begin != '\0') + ++col_begin; + return col_begin; + } + + static void unescape(char*&, char*&) + { + } +}; + +template +struct double_quote_escape { + static const char* find_next_column_end(const char* col_begin) + { + while (*col_begin != sep && *col_begin != '\0') + if (*col_begin != quote) + ++col_begin; + else { + do { + ++col_begin; + while (*col_begin != quote) { + if (*col_begin == '\0') + throw error::escaped_string_not_closed(); + ++col_begin; + } + ++col_begin; + } while (*col_begin == quote); + } + return col_begin; + } + + static void unescape(char*& col_begin, char*& col_end) + { + if (col_end - col_begin >= 2) { + if (*col_begin == quote && *(col_end - 1) == quote) { + ++col_begin; + --col_end; + char* out = col_begin; + for (char* in = col_begin; in != col_end; ++in) { + if (*in == quote && (in + 1) != col_end && *(in + 1) == quote) { + ++in; + } + *out = *in; + ++out; + } + col_end = out; + *col_end = '\0'; + } + } + } +}; + +struct throw_on_overflow { + template + static void on_overflow(T&) + { + throw error::integer_overflow(); + } + + template + static void on_underflow(T&) + { + throw error::integer_underflow(); + } +}; + +struct ignore_overflow { + template + static void on_overflow(T&) + { + } + + template + static void on_underflow(T&) + { + } +}; + +struct set_to_max_on_overflow { + template + static void on_overflow(T& x) + { + // using (std::numeric_limits::max) instead of + // std::numeric_limits::max to make code including windows.h with its max + // macro happy + x = (std::numeric_limits::max)(); + } + + template + static void on_underflow(T& x) + { + x = (std::numeric_limits::min)(); + } +}; + +namespace detail +{ +template +void chop_next_column(char*& line, char*& col_begin, char*& col_end) +{ + assert(line != nullptr); + + col_begin = line; + // the col_begin + (... - col_begin) removes the constness + col_end = col_begin + (quote_policy::find_next_column_end(col_begin) - col_begin); + + if (*col_end == '\0') { + line = nullptr; + } + else { + *col_end = '\0'; + line = col_end + 1; + } +} + +template +void parse_line(char* line, char** sorted_col, const std::vector& col_order) +{ + for (int i : col_order) { + if (line == nullptr) + throw ::io::error::too_few_columns(); + char *col_begin, *col_end; + chop_next_column(line, col_begin, col_end); + + if (i != -1) { + trim_policy::trim(col_begin, col_end); + quote_policy::unescape(col_begin, col_end); + + sorted_col[i] = col_begin; + } + } + if (line != nullptr) + throw ::io::error::too_many_columns(); +} + +template +void parse_header_line(char* line, std::vector& col_order, const std::string* col_name, + ignore_column ignore_policy) +{ + col_order.clear(); + + bool found[column_count]; + std::fill(found, found + column_count, false); + while (line) { + char *col_begin, *col_end; + chop_next_column(line, col_begin, col_end); + + trim_policy::trim(col_begin, col_end); + quote_policy::unescape(col_begin, col_end); + + for (unsigned i = 0; i < column_count; ++i) + if (col_begin == col_name[i]) { + if (found[i]) { + error::duplicated_column_in_header err; + err.set_column_name(col_begin); + throw err; + } + found[i] = true; + col_order.push_back(i); + col_begin = 0; + break; + } + if (col_begin) { + if (ignore_policy & ::io::ignore_extra_column) + col_order.push_back(-1); + else { + error::extra_column_in_header err; + err.set_column_name(col_begin); + throw err; + } + } + } + if (!(ignore_policy & ::io::ignore_missing_column)) { + for (unsigned i = 0; i < column_count; ++i) { + if (!found[i]) { + error::missing_column_in_header err; + err.set_column_name(col_name[i].c_str()); + throw err; + } + } + } +} + +template +void parse(char* col, char& x) +{ + if (!*col) + throw error::invalid_single_character(); + x = *col; + ++col; + if (*col) + throw error::invalid_single_character(); +} + +template +void parse(char* col, std::string& x) +{ + x = col; +} + +template +void parse(char* col, const char*& x) +{ + x = col; +} + +template +void parse(char* col, char*& x) +{ + x = col; +} + +template +void parse_unsigned_integer(const char* col, T& x) +{ + x = 0; + while (*col != '\0') { + if ('0' <= *col && *col <= '9') { + T y = *col - '0'; + if (x > ((std::numeric_limits::max)() - y) / 10) { + overflow_policy::on_overflow(x); + return; + } + x = 10 * x + y; + } + else + throw error::no_digit(); + ++col; + } +} + +template +void parse(char* col, unsigned char& x) +{ + parse_unsigned_integer(col, x); +} +template +void parse(char* col, unsigned short& x) +{ + parse_unsigned_integer(col, x); +} +template +void parse(char* col, unsigned int& x) +{ + parse_unsigned_integer(col, x); +} +template +void parse(char* col, unsigned long& x) +{ + parse_unsigned_integer(col, x); +} +template +void parse(char* col, unsigned long long& x) +{ + parse_unsigned_integer(col, x); +} + +template +void parse_signed_integer(const char* col, T& x) +{ + if (*col == '-') { + ++col; + + x = 0; + while (*col != '\0') { + if ('0' <= *col && *col <= '9') { + T y = *col - '0'; + if (x < ((std::numeric_limits::min)() + y) / 10) { + overflow_policy::on_underflow(x); + return; + } + x = 10 * x - y; + } + else + throw error::no_digit(); + ++col; + } + return; + } + else if (*col == '+') + ++col; + parse_unsigned_integer(col, x); +} + +template +void parse(char* col, signed char& x) +{ + parse_signed_integer(col, x); +} +template +void parse(char* col, signed short& x) +{ + parse_signed_integer(col, x); +} +template +void parse(char* col, signed int& x) +{ + parse_signed_integer(col, x); +} +template +void parse(char* col, signed long& x) +{ + parse_signed_integer(col, x); +} +template +void parse(char* col, signed long long& x) +{ + parse_signed_integer(col, x); +} + +template +void parse_float(const char* col, T& x) +{ + bool is_neg = false; + if (*col == '-') { + is_neg = true; + ++col; + } + else if (*col == '+') + ++col; + + x = 0; + while ('0' <= *col && *col <= '9') { + int y = *col - '0'; + x *= 10; + x += y; + ++col; + } + + if (*col == '.' || *col == ',') { + ++col; + T pos = 1; + while ('0' <= *col && *col <= '9') { + pos /= 10; + int y = *col - '0'; + ++col; + x += y * pos; + } + } + + if (*col == 'e' || *col == 'E') { + ++col; + int e; + + parse_signed_integer(col, e); + + if (e != 0) { + T base; + if (e < 0) { + base = T(0.1); + e = -e; + } + else { + base = T(10); + } + + while (e != 1) { + if ((e & 1) == 0) { + base = base * base; + e >>= 1; + } + else { + x *= base; + --e; + } + } + x *= base; + } + } + else { + if (*col != '\0') + throw error::no_digit(); + } + + if (is_neg) + x = -x; +} + +template +void parse(char* col, float& x) +{ + parse_float(col, x); +} +template +void parse(char* col, double& x) +{ + parse_float(col, x); +} +template +void parse(char* col, long double& x) +{ + parse_float(col, x); +} + +template +void parse(char* col, T& x) +{ + // Mute unused variable compiler warning + (void)col; + (void)x; + // GCC evaluates "false" when reading the template and + // "sizeof(T)!=sizeof(T)" only when instantiating it. This is why + // this strange construct is used. + static_assert(sizeof(T) != sizeof(T), "Can not parse this type. Only builtin integrals, floats, " + "char, char*, const char* and std::string are supported"); +} + +} // namespace detail + +template , class quote_policy = no_quote_escape<','>, + class overflow_policy = throw_on_overflow, class comment_policy = no_comment> +class CSVReader +{ +private: + LineReader in; + + char* row[column_count]; + std::string column_names[column_count]; + + std::vector col_order; + + template + void set_column_names(std::string s, ColNames... cols) + { + column_names[column_count - sizeof...(ColNames) - 1] = std::move(s); + set_column_names(std::forward(cols)...); + } + + void set_column_names() + { + } + +public: + CSVReader() = delete; + CSVReader(const CSVReader&) = delete; + CSVReader& operator=(const CSVReader&); + + template + explicit CSVReader(Args&&... args) + : in(std::forward(args)...) + { + std::fill(row, row + column_count, nullptr); + col_order.resize(column_count); + for (unsigned i = 0; i < column_count; ++i) + col_order[i] = i; + for (unsigned i = 1; i <= column_count; ++i) + column_names[i - 1] = "col" + std::to_string(i); + } + + char* next_line() + { + return in.next_line(); + } + + template + void read_header(ignore_column ignore_policy, ColNames... cols) + { + static_assert(sizeof...(ColNames) >= column_count, "not enough column names specified"); + static_assert(sizeof...(ColNames) <= column_count, "too many column names specified"); + try { + set_column_names(std::forward(cols)...); + + char* line; + do { + line = in.next_line(); + if (!line) + throw error::header_missing(); + } while (comment_policy::is_comment(line)); + + detail::parse_header_line(line, col_order, column_names, + ignore_policy); + } + catch (error::with_file_name& err) { + err.set_file_name(in.get_truncated_file_name()); + throw; + } + } + + template + void set_header(ColNames... cols) + { + static_assert(sizeof...(ColNames) >= column_count, "not enough column names specified"); + static_assert(sizeof...(ColNames) <= column_count, "too many column names specified"); + set_column_names(std::forward(cols)...); + std::fill(row, row + column_count, nullptr); + col_order.resize(column_count); + for (unsigned i = 0; i < column_count; ++i) + col_order[i] = i; + } + + bool has_column(const std::string& name) const + { + return col_order.end() != + std::find(col_order.begin(), col_order.end(), + std::find(std::begin(column_names), std::end(column_names), name) - std::begin(column_names)); + } + + void set_file_name(const std::string& file_name) + { + in.set_file_name(file_name); + } + + void set_file_name(const char* file_name) + { + in.set_file_name(file_name); + } + + const char* get_truncated_file_name() const + { + return in.get_truncated_file_name(); + } + + void set_file_line(unsigned file_line) + { + in.set_file_line(file_line); + } + + unsigned get_file_line() const + { + return in.get_file_line(); + } + +private: + void parse_helper(std::size_t) + { + } + + template + void parse_helper(std::size_t r, T& t, ColType&... cols) + { + if (row[r]) { + try { + try { + ::io::detail::parse(row[r], t); + } + catch (error::with_column_content& err) { + err.set_column_content(row[r]); + throw; + } + } + catch (error::with_column_name& err) { + err.set_column_name(column_names[r].c_str()); + throw; + } + } + parse_helper(r + 1, cols...); + } + +public: + template + bool read_row(ColType&... cols) + { + static_assert(sizeof...(ColType) >= column_count, "not enough columns specified"); + static_assert(sizeof...(ColType) <= column_count, "too many columns specified"); + try { + try { + + char* line; + do { + line = in.next_line(); + if (!line) + return false; + } while (comment_policy::is_comment(line)); + + detail::parse_line(line, row, col_order); + + parse_helper(0, cols...); + } + catch (error::with_file_name& err) { + err.set_file_name(in.get_truncated_file_name()); + throw; + } + } + catch (error::with_file_line& err) { + err.set_file_line(in.get_file_line()); + throw; + } + + return true; + } +}; +} // namespace io +#endif \ No newline at end of file diff --git a/docs/source/cpp/smm.rst b/docs/source/cpp/smm.rst index e2fb413f75..b1a1901f3c 100644 --- a/docs/source/cpp/smm.rst +++ b/docs/source/cpp/smm.rst @@ -1,11 +1,20 @@ Stochastic metapopulation model =============================== -The stochastic metapopulation model uses a Markov process to simulate disease dynamics. Similar to the `Diffusive Agent-based Model` the Markov process is given by location and infection state changes. However, in contrast to the diffusive ABM, location changes are not given by agent-dependent diffusion processes, but by stochastic jumps between regions with the requirement that the domain is split into disjoint regions. Hence, there is no further spatial resolution within one region and locations or positions are only given by the region index. The evolution of the system state is determined by the following master equation +The stochastic metapopulation model uses a Markov process to simulate disease dynamics. Similar to the :doc:`diffusive_abm` the +Markov process is given by location and infection state changes. However, in contrast to the diffusive ABM, location changes are +not given by agent-dependent diffusion processes, but by stochastic jumps between regions with the requirement that the domain +is split into disjoint regions. Hence, there is no further spatial resolution within one region and locations or positions are +only given by the region index. The evolution of the system state is determined by the following master equation :math:`\partial_t p(X,Z;t) = G p(X,Z;t) + L p(X,Z;t)`. -The operator :math:`G` defines the infection state adoptions and only acts on :math:`Z`, the vector containing all subpopulations stratified by infection state. :math:`L` defines location changes, only acting on :math:`X`, the vector containing all subpopulations stratified by region. Infection state adoptions are modeled as stochastic jumps with independent Poisson processes given by adoption rate functions. Similar to the infection state dynamics, spatial transitions between regions are also modeled as stochastic jumps with independent Poisson processes given by transition rate functions. Gillespie's direct method (stochastic simulation algorithm) is used for simulation. +The operator :math:`G` defines the infection state adoptions and only acts on :math:`Z`, the vector containing all subpopulations +stratified by infection state. :math:`L` defines location changes, only acting on :math:`X`, the vector containing all subpopulations +stratified by region. Infection state adoptions are modeled as stochastic jumps with independent Poisson processes given by adoption +rate functions. Similar to the infection state dynamics, spatial transitions between regions are also modeled as stochastic jumps +with independent Poisson processes given by transition rate functions. Gillespie's direct method (stochastic simulation algorithm) +is used for the simulation. In the following, we present more details of the stochastic metapopulation model, including code examples. An overview of nonstandard but often used data types can be found under :doc:`data_types`. @@ -13,7 +22,8 @@ An overview of nonstandard but often used data types can be found under :doc:`da Infection states ---------------- -The model does not have fixed infection states, but gets an enum class of infection states as template argument. Thus it can be used with any set of infection states. +The model does not have fixed infection states, but gets an enum class of infection states as template argument. Thus it +can be used with any set of infection states. Using the infection states Susceptible (S), Exposed (E), Carrier (C), Infected (I), Recovered (R) and Dead (D), this can be done as follows: .. code-block:: cpp @@ -32,41 +42,44 @@ Using the infection states Susceptible (S), Exposed (E), Carrier (C), Infected ( const size_t num_regions = 2; - using Model = mio::smm::Model; + using Status = mio::Index; + using Model = mio::smm::Model; Model model; Infection state transitions --------------------------- -The infection state transitions are explicitly given by the adoption rates and are therefore subject to user input. Adoption rates always depend on their source infection state. If an adoption event requires interaction of agents (e.g. disease transmission), the corresponding rate depends not only on the source infection state, but also on other infection states, the **Influence**\s. An adoption rate that only depends on the source infection state, e.g. recovery or worsening of disease symptoms, is called `first-order` adoption rate and an adoption rate that has influences is called `second-order` adoption rate. Adoption rates are region-dependent; therefore it is possible to have different rates in two regions for the same infection state transition which can be useful when having e.g. region-dependent interventions or contact behavior. +The infection state transitions are explicitly given by the adoption rates and are therefore subject to user input. +Adoption rates always depend on their source infection state. If an adoption event requires interaction of agents (e.g. +disease transmission), the corresponding rate depends not only on the source infection state, but also on other infection +states, the **Influence**\s. An adoption rate that only depends on the source infection state, e.g. recovery or worsening +of disease symptoms, is called `first-order` adoption rate and an adoption rate that has influences is called `second-order` +adoption rate. Adoption rates are region-dependent; therefore it is possible to have different rates in two regions for +the same infection state transition which can be useful when having e.g. region-dependent interventions or contact behavior. -Using the infection states from above and two regions, there are five first-order adoption rates per region and one second-order adoption rate per region. In the example below, the second-order adoption rate (transition from S to E) differs between the regions: +Using the infection states from above and two regions, there are five first-order adoption rates per region and one second-order +adoption rate per region. In the example below, the second-order adoption rate (transition from S to E) differs between the regions: .. code-block:: cpp - std::vector> adoption_rates; + std::vector> adoption_rates; //Set first-order adoption rates for both regions for (size_t r = 0; r < num_regions; ++r) { - adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}}); + adoption_rates.push_back({{InfectionState::E}, {InfectionState::C}, mio::regions::Region(r), 1.0 / 5., {}}); + adoption_rates.push_back({{InfectionState::C}, {InfectionState::R}, mio::regions::Region(r), 0.2 / 3., {}}); + adoption_rates.push_back({{InfectionState::C}, {InfectionState::I}, mio::regions::Region(r), 0.8 / 3., {}}); + adoption_rates.push_back({{InfectionState::I}, {InfectionState::R}, mio::regions::Region(r), 0.99 / 5., {}}); + adoption_rates.push_back({{InfectionState::I}, {InfectionState::D}, mio::regions::Region(r), 0.01 / 5., {}}); } //Set second-order adoption rate different for the two regions // adoption rate has the form {i, j, k, c_{i,j}, {{tau1.state, tau1.factor}, {tau2.state, tau2.factor}}}, see the equation below - adoption_rates.push_back({InfectionState::S, InfectionState::E, mio::regions::Region(0), 0.1, {{InfectionState::C, 1}, {InfectionState::I, 0.5}}}); - adoption_rates.push_back({InfectionState::S, InfectionState::E, mio::regions::Region(1), 0.2, {{InfectionState::C, 1}, {InfectionState::I, 0.5}}}); + adoption_rates.push_back({{InfectionState::S}, {InfectionState::E}, mio::regions::Region(0), 0.1, {{{InfectionState::C}, 1}, {{InfectionState::I}, 0.5}}}); + adoption_rates.push_back({{InfectionState::S}, {InfectionState::E}, mio::regions::Region(1), 0.2, {{{InfectionState::C}, 1}, {{InfectionState::I}, 0.5}}}); //Initialize model parameter - model.parameters.get>() = adoption_rates; - -Sociodemographic stratification -------------------------------- - -Sociodemographic stratification e.g. by age, gender or immunity can be incorporated by stratifying the set of infection states passed as template to the model. + model.parameters.get>() = adoption_rates; Parameters ---------- @@ -110,8 +123,9 @@ with :math:`i^{(k)}` the population in region :math:`k` having infection state : Initial conditions ------------------ -Before running a simulation with the stochastic metapopulation model, the initial populations i.e. the number of agents per infection state for every region have to be set. -These populations have the class type **Populations** and can be set via: +Before running a simulation with the stochastic metapopulation model, the initial populations i.e. the number of agents +per infection state for every region have to be set. +These populations have the class type ``Populations`` and can be set via: .. code-block:: cpp @@ -128,7 +142,8 @@ These populations have the class type **Populations** and can be set via: } If individuals should transition between regions, the spatial transition rates of the model have to be initialized as well. -As the spatial transition rates are dependent on infection state, region changes for specific infection states can be prevented. Below, symmetric spatial transition rates are set for every region: +As the spatial transition rates are dependent on infection state, region changes for specific infection states can be +prevented. Below, symmetric spatial transition rates are set for every region: .. code-block:: cpp @@ -140,27 +155,32 @@ As the spatial transition rates are dependent on infection state, region changes if (i != j) { // transition rate has the form {i, k, l, \lambda^{(k,l)}_{i}} transition_rates.push_back( - {InfectionState(s), mio::regions::Region(i), mio::regions::Region(j), 0.01}); + {{InfectionState(s)}, mio::regions::Region(i), mio::regions::Region(j), 0.01}); transition_rates.push_back( - {InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01}); + {{InfectionState(s)}, mio::regions::Region(j), mio::regions::Region(i), 0.01}); } } } //Initialize model parameter - model.parameters.get>() = transition_rates; + model.parameters.get>() = transition_rates; Nonpharmaceutical interventions -------------------------------- -There are no nonpharmaceutical interventions (NPIs) explicitly implemented in the model. However, NPIs influencing the adoption or spatial transition rates can be realized by resetting the corresponding model parameters. +There are no nonpharmaceutical interventions (NPIs) explicitly implemented in the model. However, NPIs influencing the +adoption or spatial transition rates can be realized by resetting the corresponding model parameters. Simulation ---------- -At the beginning of the simulation, the waiting times for all events (infection state adoptions and spatial transitions) are drawn. Then the time is advanced until the time point of the next event - which can be a spatial transition or an infection state adoption - and the event takes places. The waiting times of the other events are updated and a new waiting time for the event that just happened is drawn. The simulation saves the system state in discrete time steps. +At the beginning of the simulation, the waiting times for all events (infection state adoptions and spatial transitions) +are drawn. Then the time is advanced until the time point of the next event - which can be a spatial transition or an +infection state adoption - and the event takes places. The waiting times of the other events are updated and a new waiting +time for the event that just happened is drawn. The simulation saves the system state in discrete time steps. -To simulate the model from `t0` to `tmax` with given step size `dt`, a **Simulation** has to be created and advanced until `tmax`. The step size is only used to regularly save the system state during the simulation. +To simulate the model from ``t0`` to ``tmax`` with given step size ``dt``, a **Simulation** has to be created and advanced +until ``tmax``. The step size is only used to regularly save the system state during the simulation. .. code-block:: cpp @@ -193,10 +213,58 @@ If one wants to interpolate the aggregated results to a ``mio::TimeSeries`` cont auto interpolated_results = mio::interpolate_simulation_result(sim.get_result()); + +Demographic Stratification +-------------------------- + +It is possible to add multiple indices to the ``Status`` to differentiate multiple groups on the same region. For example this +could represent the human and the mosquito populations in a specific region. To use this feature, one first of all has to +create a new index: + +.. code-block:: cpp + + struct Species : public mio::Index { + Species(size_t val): Index(val) + { + } + }; + +Then, one has to create the multiindex, where we reuse the ``InfectionState`` defined in the first example: + +.. code-block:: cpp + + using Status = mio::Index; + +Define the size for each index dimension that is not an enum class: + +.. code-block:: cpp + + const size_t num_species = 2; + +We can define a model: + +.. code-block:: cpp + + using Model = mio::smm::Model + Model model(Status{Count, Species(num_species)}, Region(num_regions)); + +Now, for accessing the population, all indexes need to be given: + +-- code-block:: cpp + + model.populations[{Region(r), InfectionState::S, Species(0)}] = 100; + // ... + adoption_rates.push_back({{InfectionState::S, Species(0)}, {InfectionState::E, Species(0)}, Region(r), 0.1, {{InfectionState::I, Species(1)}, 1}); + // ... + transition_rates.push_back({{InfectionState::S, Species(0)}, Region(i), Region(j), 0.01}); + + Examples -------- -An example of the stochastic metapopulation model with four regions can be found at: `examples/smm.cpp `_ +A full example with ``Status`` distributed according to ``InfectionState``, ``Age`` and ``Species`` can be found at +`examples/smm.cpp `_ + Overview of the ``smm`` namespace: