diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 3c5e4e7369..e996d0a39c 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -21,6 +21,8 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI add_executable(ode_sir_example ode_sir.cpp) target_link_libraries(ode_sir_example PRIVATE memilio ode_sir) +add_executable(sim_mobility_detailed ode_secirvvs_mobility_detailed.cpp) +target_link_libraries(sim_mobility_detailed PRIVATE memilio ode_secirvvs) add_executable(seir_flows_example ode_seir_flows.cpp) target_link_libraries(seir_flows_example PRIVATE memilio ode_seir) diff --git a/cpp/examples/ode_secirvvs_mobility_detailed.cpp b/cpp/examples/ode_secirvvs_mobility_detailed.cpp new file mode 100644 index 0000000000..ad63b1c309 --- /dev/null +++ b/cpp/examples/ode_secirvvs_mobility_detailed.cpp @@ -0,0 +1,632 @@ +/* +* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC) +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#include "memilio/compartments/flow_simulation.h" +#include "memilio/config.h" +#include "memilio/data/analyze_result.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/logging.h" +#include "ode_secirvvs/model.h" +#include "ode_secirvvs/infection_state.h" +#include "ode_secirvvs/parameters.h" +#include "ode_secirvvs/parameters_io.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" +#include "memilio/compartments/simulation.h" +#include "memilio/io/result_io.h" +#include "memilio/compartments/parameter_studies.h" +#include "memilio/utils/miompi.h" +#include "memilio/utils/random_number_generator.h" + +#include "memilio/geography/regions.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/mobility_io.h" +#include "memilio/utils/date.h" +#include "ode_secirvvs/parameter_space.h" +#include "memilio/utils/stl_util.h" +#include "boost/filesystem.hpp" +#include +#include +#include + +#include +#include +#include + +/** + * Set a value and distribution of an UncertainValue. + * Assigns average of min and max as a value and UNIFORM(min, max) as a distribution. + * @param p uncertain value to set. + * @param min minimum of distribution. + * @param max minimum of distribution. + */ +void assign_uniform_distribution(mio::UncertainValue& p, double min, double max) +{ + p = mio::UncertainValue(0.5 * (max + min)); + p.set_distribution(mio::ParameterDistributionUniform(min, max)); +} + +/** + * Set a value and distribution of an array of UncertainValues. + * Assigns average of min[i] and max[i] as a value and UNIFORM(min[i], max[i]) as a distribution for + * each element i of the array. + * @param array array of UncertainValues to set. + * @param min minimum of distribution for each element of array. + * @param max minimum of distribution for each element of array. + */ +template +void array_assign_uniform_distribution(mio::CustomIndexArray& array, + const double (&min)[N], const double (&max)[N]) +{ + assert(N == array.numel()); + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(N); ++i) { + assign_uniform_distribution(array[i], min[size_t(i)], max[size_t(i)]); + } +} + +/** + * Set a value and distribution of an array of UncertainValues. + * Assigns average of min and max as a value and UNIFORM(min, max) as a distribution to every element of the array. + * @param array array of UncertainValues to set. + * @param min minimum of distribution. + * @param max minimum of distribution. + */ +void array_assign_uniform_distribution(mio::CustomIndexArray& array, double min, + double max) +{ + for (auto i = mio::AgeGroup(0); i < array.size(); ++i) { + assign_uniform_distribution(array[i], min, max); + } +} + +/** + * Set epidemiological parameters of Covid19. + * @param params Object that the parameters will be added to. + * @returns Currently generates no errors. + */ +mio::IOResult set_covid_parameters(mio::osecirvvs::Parameters& params) +{ + //times + const double incubationTime = 3.1; // doi.org/10.3201/eid2806.220158 + const double serialIntervalMin = 2.38; // doi.org/10.1016/j.lanepe.2022.100446 + const double serialIntervalMax = 2.38; // doi.org/10.1016/j.lanepe.2022.100446 + const double timeInfectedSymptomsMin = 6.58; //https://doi.org/10.1016/S0140-6736(22)00327-0 + const double timeInfectedSymptomsMax = 7.16; //https://doi.org/10.1016/S0140-6736(22)00327-0 + const double timeInfectedSevereMin[] = {1.8, 1.8, 1.8, 2.5, 3.5, 4.91}; // doi.org/10.1186/s12879-022-07971-6 + const double timeInfectedSevereMax[] = {2.3, 2.3, 2.3, 3.67, 5, 7.01}; // doi.org/10.1186/s12879-022-07971-6 + + const ScalarType fact_icu = 1.; + const double timeInfectedCriticalMin[] = {fact_icu * 4.95, fact_icu * 4.95, fact_icu * 4.86, + fact_icu * 14.14, fact_icu * 14.4, fact_icu * 10.}; + const double timeInfectedCriticalMax[] = {fact_icu * 8.95, fact_icu * 8.95, fact_icu * 8.86, + fact_icu * 20.58, fact_icu * 19.8, fact_icu * 13.2}; + + array_assign_uniform_distribution(params.get(), incubationTime, incubationTime); + array_assign_uniform_distribution(params.get(), serialIntervalMin, + serialIntervalMax); + array_assign_uniform_distribution(params.get(), timeInfectedSymptomsMin, + timeInfectedSymptomsMax); + array_assign_uniform_distribution(params.get(), timeInfectedSevereMin, + timeInfectedSevereMax); + array_assign_uniform_distribution(params.get(), timeInfectedCriticalMin, + timeInfectedCriticalMax); + + double fac_variant = 1.94; //https://doi.org/10.7554/eLife.78933 + const double transmissionProbabilityOnContactMin[] = {0.02 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, + 0.05 * fac_variant, 0.08 * fac_variant, 0.1 * fac_variant}; + + const double transmissionProbabilityOnContactMax[] = {0.04 * fac_variant, 0.07 * fac_variant, 0.07 * fac_variant, + 0.07 * fac_variant, 0.10 * fac_variant, 0.15 * fac_variant}; + + const double relativeTransmissionNoSymptomsMin = 0.5; + const double relativeTransmissionNoSymptomsMax = 0.5; + const double riskOfInfectionFromSymptomaticMin = 0.0; // beta (depends on incidence and test and trace capacity) + const double riskOfInfectionFromSymptomaticMax = 0.2; + const double maxRiskOfInfectionFromSymptomaticMin = 0.4; + const double maxRiskOfInfectionFromSymptomaticMax = 0.5; + + const double recoveredPerInfectedNoSymptomsMin[] = {0.2, 0.25, 0.2, + 0.2, 0.175, 0.1}; // doi.org/10.1101/2022.05.05.22274697 + const double recoveredPerInfectedNoSymptomsMax[] = {0.4, 0.45, 0.35, + 0.3, 0.25, 0.15}; // doi.org/10.1101/2022.05.05.22274697 + + // doi:10.1136/bmjgh-2023-0123 + + // Factors from https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(22)00462-7/fulltext + const double severePerInfectedSymptomsMin[] = {1 * 0.006, 0.8 * 0.006, 0.4 * 0.015, + 0.3 * 0.049, 0.25 * 0.15, 0.35 * 0.2}; + const double severePerInfectedSymptomsMax[] = {1 * 0.009, 0.8 * 0.009, 0.4 * 0.023, 0.3 * 0.074, + 0.25 * 0.18, 0.35 * 0.25}; // quelle 2021 paper + factors + + // delta paper + // risk of icu admission https://doi.org/10.1177/14034948221108548 + const double fac_icu = 0.52; + const double criticalPerSevereMin[] = {0.05 * fac_icu, 0.05 * fac_icu, 0.05 * fac_icu, + 0.10 * fac_icu, 0.25 * fac_icu, 0.35 * fac_icu}; + const double criticalPerSevereMax[] = {0.10 * fac_icu, 0.10 * fac_icu, 0.10 * fac_icu, + 0.20 * fac_icu, 0.35 * fac_icu, 0.45 * fac_icu}; + + // 61% less risk doi:10.1136/bmjgh-2023-0123 + // https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10347449/pdf/bmjgh-2023-012328.pdf + const double fac_dead = 0.39; + const double deathsPerCriticalMin[] = {fac_dead * 0.00, fac_dead * 0.00, fac_dead * 0.10, + fac_dead * 0.10, fac_dead * 0.30, fac_dead * 0.5}; // 2021 paper + const double deathsPerCriticalMax[] = {fac_dead * 0.10, fac_dead * 0.10, fac_dead * 0.18, + fac_dead * 0.18, fac_dead * 0.50, fac_dead * 0.7}; + + const double reducExposedPartialImmunityMin = 0.569; // doi.org/10.1136/bmj-2022-071502 + const double reducExposedPartialImmunityMax = 0.637; // doi.org/10.1136/bmj-2022-071502 + const double reducExposedImprovedImmunityMin = + 0.34 * reducExposedPartialImmunityMin; // https://jamanetwork.com/journals/jama/fullarticle/2788487 0.19346 + const double reducExposedImprovedImmunityMax = + 0.34 * reducExposedPartialImmunityMax; // https://jamanetwork.com/journals/jama/fullarticle/2788487 0.21658 + + const double reducInfectedSymptomsPartialImmunityMin = 0.746; // doi.org/10.1056/NEJMoa2119451 + const double reducInfectedSymptomsPartialImmunityMax = 0.961; // doi.org/10.1056/NEJMoa2119451 + const double reducInfectedSymptomsImprovedImmunityMin = 0.295; // doi.org/10.1056/NEJMoa2119451 + const double reducInfectedSymptomsImprovedImmunityMax = 0.344; // doi.org/10.1056/NEJMoa2119451 + + const double reducInfectedSevereCriticalDeadPartialImmunityMin = + 0.52; // www.assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/1050721/Vaccine-surveillance-report-week-4.pdf + const double reducInfectedSevereCriticalDeadPartialImmunityMax = + 0.82; // www.assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/1050721/Vaccine-surveillance-report-week-4.pdf + const double reducInfectedSevereCriticalDeadImprovedImmunityMin = 0.1; // doi.org/10.1136/bmj-2022-071502 + const double reducInfectedSevereCriticalDeadImprovedImmunityMax = 0.19; // doi.org/10.1136/bmj-2022-071502 + const double temp_reducTimeInfectedMild = 0.5; // doi.org/10.1101/2021.09.24.21263978 + + array_assign_uniform_distribution(params.get(), + transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); + array_assign_uniform_distribution(params.get(), + relativeTransmissionNoSymptomsMin, relativeTransmissionNoSymptomsMax); + array_assign_uniform_distribution(params.get(), + riskOfInfectionFromSymptomaticMin, riskOfInfectionFromSymptomaticMax); + array_assign_uniform_distribution(params.get(), + maxRiskOfInfectionFromSymptomaticMin, maxRiskOfInfectionFromSymptomaticMax); + array_assign_uniform_distribution(params.get(), + recoveredPerInfectedNoSymptomsMin, recoveredPerInfectedNoSymptomsMax); + array_assign_uniform_distribution(params.get(), + severePerInfectedSymptomsMin, severePerInfectedSymptomsMax); + array_assign_uniform_distribution(params.get(), criticalPerSevereMin, + criticalPerSevereMax); + array_assign_uniform_distribution(params.get(), deathsPerCriticalMin, + deathsPerCriticalMax); + + array_assign_uniform_distribution(params.get(), + reducExposedPartialImmunityMin, reducExposedPartialImmunityMax); + array_assign_uniform_distribution(params.get(), + reducExposedImprovedImmunityMin, reducExposedImprovedImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSymptomsPartialImmunityMin, reducInfectedSymptomsPartialImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSymptomsImprovedImmunityMin, + reducInfectedSymptomsImprovedImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSevereCriticalDeadPartialImmunityMin, + reducInfectedSevereCriticalDeadPartialImmunityMax); + array_assign_uniform_distribution(params.get(), + reducInfectedSevereCriticalDeadImprovedImmunityMin, + reducInfectedSevereCriticalDeadImprovedImmunityMax); + array_assign_uniform_distribution(params.get(), temp_reducTimeInfectedMild, + temp_reducTimeInfectedMild); + + //sasonality + const double seasonality_min = 0.1; + const double seasonality_max = 0.3; + + assign_uniform_distribution(params.get(), seasonality_min, seasonality_max); + return mio::success(); +} + +enum class ContactLocation +{ + Home = 0, + School, + Work, + Other, + Count, +}; + +/** + * different types of NPI, used as DampingType. + */ +enum class Intervention +{ + Home, + SchoolClosure, + HomeOffice, + GatheringBanFacilitiesClosure, + PhysicalDistanceAndMasks, + SeniorAwareness, +}; + +/** + * different level of NPI, used as DampingLevel. + */ +enum class InterventionLevel +{ + Main, + PhysicalDistanceAndMasks, + SeniorAwareness, + Holidays, +}; + +static const std::map contact_locations = {{ContactLocation::Home, "home"}, + {ContactLocation::School, "school_pf_eig"}, + {ContactLocation::Work, "work"}, + {ContactLocation::Other, "other"}}; + +/** + * Set contact matrices. + * Reads contact matrices from files in the data directory. + * @param data_dir data directory. + * @param params Object that the contact matrices will be added to. + * @returns any io errors that happen during reading of the files. + */ +mio::IOResult set_contact_matrices(const fs::path& data_dir, mio::osecirvvs::Parameters& params, + ScalarType scale_contacts = 1.0, ScalarType share_staying = 1.0) +{ + auto contact_transport_status = mio::read_mobility_plain(data_dir.string() + "//contacts//contacts_transport.txt"); + auto contact_matrix_transport = contact_transport_status.value(); + auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), size_t(params.get_num_groups())); + for (auto&& contact_location : contact_locations) { + BOOST_OUTCOME_TRY(baseline, + mio::read_mobility_plain( + (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string())); + + if (contact_location.second == "other") { + contact_matrices[size_t(contact_location.first)].get_baseline() = + (1 - share_staying) * abs(baseline - contact_matrix_transport) / scale_contacts + + share_staying * abs(baseline - contact_matrix_transport); + contact_matrices[size_t(contact_location.first)].get_minimum() = Eigen::MatrixXd::Zero(6, 6); + } + else { + contact_matrices[size_t(contact_location.first)].get_baseline() = + (1 - share_staying) * baseline / scale_contacts + share_staying * baseline; + contact_matrices[size_t(contact_location.first)].get_minimum() = Eigen::MatrixXd::Zero(6, 6); + } + } + params.get() = mio::UncertainContactMatrix(contact_matrices); + + return mio::success(); +} + +mio::IOResult set_contact_matrices_transport(const fs::path& data_dir, mio::osecirvvs::Parameters& params, + ScalarType scale_contacts = 1.0) +{ + auto contact_transport_status = mio::read_mobility_plain(data_dir.string() + "//contacts//contacts_transport.txt"); + auto contact_matrix_transport = contact_transport_status.value(); + auto contact_matrices = mio::ContactMatrixGroup(1, size_t(params.get_num_groups())); + // ScalarType const polymod_share_contacts_transport = 1 / 0.2770885028949545; + + contact_matrices[0].get_baseline() = contact_matrix_transport / scale_contacts; + contact_matrices[0].get_minimum() = Eigen::MatrixXd::Zero(6, 6); + + params.get() = mio::UncertainContactMatrix(contact_matrices); + + return mio::success(); +} + +// reset population in graph +void init_pop_cologne_szenario(mio::Graph& graph) +{ + std::vector> immunity = {{0.04, 0.61, 0.35}, {0.04, 0.61, 0.35}, {0.075, 0.62, 0.305}, + {0.08, 0.62, 0.3}, {0.035, 0.58, 0.385}, {0.01, 0.41, 0.58}}; + + for (auto& node : graph.nodes()) { + for (auto age = mio::AgeGroup(0); age < mio::AgeGroup(6); age++) { + auto pop_age = 0.0; + for (auto inf_state = mio::Index(0); + inf_state < mio::osecirvvs::InfectionState::Count; ++inf_state) { + pop_age += node.property.populations[{age, inf_state}]; + node.property.populations[{age, inf_state}] = 0.0; + } + size_t immunity_index = static_cast(age); + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleNaive}] = + pop_age * immunity[immunity_index][0]; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptiblePartialImmunity}] = + pop_age * immunity[immunity_index][1]; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity}] = + pop_age * immunity[immunity_index][2]; + } + if (node.id == 5315) { + // infect p% of population + ScalarType p = 0.05; + for (auto age = mio::AgeGroup(0); age < graph.nodes()[0].property.parameters.get_num_groups(); age++) { + node.property.populations[{mio::AgeGroup(age), mio::osecirvvs::InfectionState::InfectedSymptomsNaive}] = + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleNaive}] * p; + node.property.populations[{age, mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity}] = + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptiblePartialImmunity}] * p; + node.property.populations[{mio::AgeGroup(age), + mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity}] = + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity}] * p; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleNaive}] *= 0.99; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptiblePartialImmunity}] *= 0.99; + node.property.populations[{age, mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity}] *= 0.99; + } + } + } +} + +/** + * Create the input graph for the parameter study. + * Reads files from the data directory. + * @param start_date start date of the simulation. + * @param end_date end date of the simulation. + * @param data_dir data directory. + * @returns created graph or any io errors that happen during reading of the files. + */ +mio::IOResult> +get_graph(const std::string data_dir, const int num_days, bool masks, bool ffp2, bool szenario_cologne, bool edges) +{ + auto mobility_data_commuter = + mio::read_mobility_plain(("/localdata1/code/memilio/test/commuter_migration_with_locals.txt")); + auto mob_data = mobility_data_commuter.value(); + auto start_date = mio::Date(2021, 8, 1); + auto end_date = mio::Date(2021, 11, 1); + + // global parameters + const int num_age_groups = 6; + mio::GraphDetailed params_graph; + + // Nodes + mio::osecirvvs::Parameters params(num_age_groups); + + params.get() = mio::get_day_in_year(start_date); + auto params_status = set_covid_parameters(params); + auto contacts_status = set_contact_matrices(data_dir, params); + + params.get() = mio::get_day_in_year(start_date); + + // Set nodes + auto scaling_factor_infected = std::vector(size_t(params.get_num_groups()), 1.0); + auto scaling_factor_icu = 1.0; + auto tnt_capacity_factor = 1.43 / 100000.; + + const auto& read_function_nodes = mio::osecirvvs::read_input_data_county; + const auto& read_function_edges = mio::read_mobility_plain; + const auto& node_id_function = mio::get_node_ids; + + const auto& set_node_function = + mio::set_nodes_detailed; + BOOST_OUTCOME_TRY(set_node_function( + params, start_date, end_date, data_dir, data_dir + "//pydata//Germany//county_current_population.json", + mio::path_join(data_dir, "mobility", "activity_duration_work.txt"), true, params_graph, read_function_nodes, + node_id_function, scaling_factor_infected, scaling_factor_icu, tnt_capacity_factor, + mio::get_offset_in_days(end_date, start_date), false)); + + for (auto& node : params_graph.nodes()) { + if (masks) { + + const double fact_surgical_mask = 0.1; + const double fact_ffp2 = 0.001; + + double factor_mask[] = {1, fact_surgical_mask, fact_ffp2, fact_ffp2, fact_ffp2, fact_ffp2}; + if (!ffp2) { + std::cout << "surgical masks" + << "\n"; + // multipliziere alle außer den ersten EIntrag mit 40 + for (size_t j = 1; j < 6; j++) { + factor_mask[j] = factor_mask[j] * fact_surgical_mask / fact_ffp2; + } + } + + double fac_variant = 1.94; //https://doi.org/10.7554/eLife.78933 + double transmissionProbabilityOnContactMin[] = {0.02 * fac_variant, 0.05 * fac_variant, 0.05 * fac_variant, + 0.05 * fac_variant, 0.08 * fac_variant, 0.1 * fac_variant}; + + double transmissionProbabilityOnContactMax[] = {0.04 * fac_variant, 0.07 * fac_variant, 0.07 * fac_variant, + 0.07 * fac_variant, 0.10 * fac_variant, 0.15 * fac_variant}; + for (int i = 0; i < 6; i++) { + transmissionProbabilityOnContactMin[i] = transmissionProbabilityOnContactMin[i] * factor_mask[i]; + transmissionProbabilityOnContactMax[i] = transmissionProbabilityOnContactMax[i] * factor_mask[i]; + } + array_assign_uniform_distribution( + node.mobility.parameters.get(), + transmissionProbabilityOnContactMin, transmissionProbabilityOnContactMax); + } + + for (size_t t_idx = 0; t_idx < num_days; ++t_idx) { + auto t = mio::SimulationDay((size_t)t_idx); + for (auto j = mio::AgeGroup(0); j < params.get_num_groups(); j++) { + node.mobility.parameters.template get()[{j, t}] = 0; + node.mobility.parameters.template get()[{j, t}] = 0; + } + } + } + + if (szenario_cologne) { + init_pop_cologne_szenario(params_graph); + } + + // Edges + if (edges) { + auto migrating_compartments = { + mio::osecirvvs::InfectionState::SusceptibleNaive, + mio::osecirvvs::InfectionState::ExposedNaive, + mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive, + mio::osecirvvs::InfectionState::InfectedSymptomsNaive, + mio::osecirvvs::InfectionState::SusceptibleImprovedImmunity, + mio::osecirvvs::InfectionState::SusceptiblePartialImmunity, + mio::osecirvvs::InfectionState::ExposedPartialImmunity, + mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity, + mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity, + mio::osecirvvs::InfectionState::ExposedImprovedImmunity, + mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity, + mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity, + }; + + const auto& read_function_edges = mio::read_mobility_plain; + const auto& set_edge_function = + mio::set_edges_detailed; + BOOST_OUTCOME_TRY(set_edge_function(mio::path_join(data_dir, "mobility", "travel_times_pathes.txt"), + data_dir + "//mobility//commuter_migration_with_locals.txt", + mio::path_join(data_dir, "mobility", "wegketten_ohne_komma.txt"), + params_graph, migrating_compartments, contact_locations.size(), + std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}, false)); + + // average travel time + const ScalarType avg_traveltime = std::accumulate(params_graph.edges().begin(), params_graph.edges().end(), 0.0, + [](double sum, const auto& e) { + return sum + e.traveltime; + }) / + params_graph.edges().size(); + + // average share of commuters for all counties relative to the total population + const ScalarType total_population = std::accumulate(params_graph.nodes().begin(), params_graph.nodes().end(), + 0.0, [](double sum, const auto& n) { + return sum + n.property.populations.get_total(); + }); + const ScalarType num_commuters = std::accumulate(params_graph.edges().begin(), params_graph.edges().end(), 0.0, + [mob_data](double sum, const auto& e) { + auto start_node = e.start_node_idx; + auto end_node = e.end_node_idx; + return sum + mob_data(start_node, end_node); + }); + const ScalarType avg_commuter_share = num_commuters / total_population; + + if (mio::mpi::is_root()) { + std::cout << "avg commuter share = " << avg_commuter_share << "\n"; + std::cout << "avg travel time = " << avg_traveltime << "\n"; + } + + // scale all contact matrices + for (auto& node : params_graph.nodes()) { + // mobility node + auto contacts_status_mobility = + set_contact_matrices_transport(data_dir, node.mobility.parameters, avg_traveltime * avg_commuter_share); + + // local node + auto contacts_status = + set_contact_matrices(data_dir, node.property.parameters, 1 - avg_traveltime, 1 - avg_commuter_share); + } + mio::ContactMatrixGroup const& contact_matrix3 = + params_graph.nodes()[0].mobility.parameters.get(); + std::cout << "Mobility node: contact matrix at t = 0\n"; + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(6); i++) { + for (auto j = mio::AgeGroup(0); j < mio::AgeGroup(6); j++) { + std::cout << contact_matrix3.get_matrix_at(0)(static_cast((size_t)i), + static_cast((size_t)j)) + << " "; + } + std::cout << "\n"; + } + } + + if (mio::mpi::is_root()) { + std::cout << "Anzahl kanten = " << params_graph.edges().size() << "\n"; + } + return params_graph; +} +mio::IOResult run() +{ + // mio::set_log_level(mio::LogLevel::critical); + const auto num_days = 70; + const int num_runs = 150; + const bool masks = true; + const bool ffp2 = true; + const bool edges = true; + + // wenn masks false und ffp2 true, dann error ausgeben + if (!masks && ffp2) { + mio::log_error("ffp2 only possible with masks"); + } + std::string data_dir = "/localdata1/test/memilio/data"; + + bool szenario_cologne = false; + + BOOST_OUTCOME_TRY(created, get_graph(data_dir, num_days, masks, ffp2, szenario_cologne, edges)); + auto params_graph = created; + + std::string res_dir = "/localdata1/code/memilio/results_paper/mask_" + std::to_string(masks); + + if (ffp2) + res_dir = res_dir + "_ffp2"; + if (szenario_cologne) + res_dir = res_dir + "_cologne"; + if (!edges) + res_dir = res_dir + "_no_edges"; + if (!edges && (ffp2 || masks)) + mio::log_error("no edges only possible without masks and ffp2"); + + res_dir += std::string(ffp2 ? "_ffp2" : "") + std::string(szenario_cologne ? "_cologne" : "") + + std::string(!edges ? "_no_edges" : ""); + + // check if boths dir exist, otherwise create them + if (!fs::exists(res_dir)) { + fs::create_directory(res_dir); + } + + std::vector county_ids(params_graph.nodes().size()); + std::transform(params_graph.nodes().begin(), params_graph.nodes().end(), county_ids.begin(), [](auto& n) { + return n.id; + }); + + // parameter study + auto parameter_study = mio::ParameterStudyDetailed>( + params_graph, 0.0, num_days, 0.01, num_runs); + if (mio::mpi::is_root()) { + printf("Seeds: "); + for (auto s : parameter_study.get_rng().get_seeds()) { + printf("%u, ", s); + } + printf("\n"); + } + auto save_single_run_result = mio::IOResult(mio::success()); + auto ensemble = parameter_study.run( + [&](auto&& graph) { + return draw_sample(graph); + }, + [&](auto results_graph, auto&& run_idx) { + auto interpolated_result = mio::interpolate_simulation_result(results_graph); + auto params = std::vector(); + + return std::make_pair(interpolated_result, params); + }); + + if (ensemble.size() > 0) { + auto ensemble_results = std::vector>>{}; + ensemble_results.reserve(ensemble.size()); + auto ensemble_params = std::vector>{}; + ensemble_params.reserve(ensemble.size()); + for (auto&& run : ensemble) { + ensemble_results.emplace_back(std::move(run.first)); + ensemble_params.emplace_back(std::move(run.second)); + } + BOOST_OUTCOME_TRY(save_results(ensemble_results, ensemble_params, county_ids, res_dir, false)); + } + + return mio::success(); +} + +int main(int, char**) +{ + mio::set_log_level(mio::LogLevel::warn); + mio::mpi::init(); + + auto result = run(); + if (!result) { + printf("%s\n", result.error().formatted_message().c_str()); + mio::mpi::finalize(); + return -1; + } + mio::mpi::finalize(); + return 0; +} \ No newline at end of file diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt index be65a209c4..775da2bf7f 100644 --- a/cpp/memilio/CMakeLists.txt +++ b/cpp/memilio/CMakeLists.txt @@ -59,6 +59,8 @@ add_library(memilio mobility/metapopulation_mobility_instant.cpp mobility/metapopulation_mobility_stochastic.h mobility/metapopulation_mobility_stochastic.cpp + mobility/metapopulation_mobility_detailed.cpp + mobility/metapopulation_mobility_detailed.h 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 73ec8fa461..893dadd14f 100644 --- a/cpp/memilio/compartments/parameter_studies.h +++ b/cpp/memilio/compartments/parameter_studies.h @@ -38,6 +38,22 @@ namespace mio { +template +struct has_stay_duration : std::false_type { +}; + +template +struct has_stay_duration().stay_duration)>> : std::true_type { +}; + +template +struct has_traveltime : std::false_type { +}; + +template +struct has_traveltime().traveltime)>> : std::true_type { +}; + /** * Class that performs multiple simulation runs with randomly sampled parameters. * Can simulate migration graphs with one simulation in each node or single simulations. @@ -140,8 +156,8 @@ class ParameterStudy #else num_procs = 1; rank = 0; -#endif - +#endif + //The ParameterDistributions used for sampling parameters use thread_local_rng() //So we set our own RNG to be used. //Assume that sampling uses the thread_local_rng() and isn't multithreaded @@ -149,7 +165,8 @@ class ParameterStudy thread_local_rng() = m_rng; auto run_distribution = distribute_runs(m_num_runs, num_procs); - auto start_run_idx = std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0)); + auto start_run_idx = + std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0)); auto end_run_idx = start_run_idx + run_distribution[size_t(rank)]; std::vector> ensemble_result; @@ -167,7 +184,8 @@ class ParameterStudy //sample auto sim = create_sampled_simulation(sample_graph); - log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", (thread_local_rng().get_counter() - run_rng_counter).get()); + log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", + (thread_local_rng().get_counter() - run_rng_counter).get()); //perform run sim.advance(m_tmax); @@ -328,10 +346,31 @@ class ParameterStudy } /** @} */ - RandomNumberGenerator& get_rng() { + RandomNumberGenerator& get_rng() + { return m_rng; } +protected: + // adaptive time step of the integrator (will be corrected if too large/small) + double m_dt_integration = 0.1; + // + RandomNumberGenerator m_rng; + + std::vector distribute_runs(size_t num_runs, int num_procs) + { + //evenly distribute runs + //lower processes do one more run if runs are not evenly distributable + auto num_runs_local = num_runs / num_procs; //integer division! + auto remainder = num_runs % num_procs; + + std::vector run_distribution(num_procs); + std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); + std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local); + + return run_distribution; + } + private: //sample parameters and create simulation template @@ -341,29 +380,25 @@ class ParameterStudy auto sampled_graph = sample_graph(m_graph); for (auto&& node : sampled_graph.nodes()) { - sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration); + if constexpr (has_stay_duration::value) { + sim_graph.add_node(node.id, node.stay_duration, node.property, m_t0, m_dt_integration); + } + else { + sim_graph.add_node(node.id, node.property, m_t0, m_dt_integration); + } } for (auto&& edge : sampled_graph.edges()) { - sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + if constexpr (has_traveltime::value) { + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.property); + } + else { + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property); + } } return make_migration_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); } - std::vector distribute_runs(size_t num_runs, int num_procs) - { - //evenly distribute runs - //lower processes do one more run if runs are not evenly distributable - auto num_runs_local = num_runs / num_procs; //integer division! - auto remainder = num_runs % num_procs; - - std::vector run_distribution(num_procs); - std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1); - std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local); - - return run_distribution; - } - private: // Stores Graph with the names and ranges of all parameters ParametersGraph m_graph; @@ -376,10 +411,6 @@ class ParameterStudy double m_tmax; // time step of the graph double m_dt_graph_sim; - // adaptive time step of the integrator (will be corrected if too large/small) - double m_dt_integration = 0.1; - // - RandomNumberGenerator m_rng; }; } // namespace mio diff --git a/cpp/memilio/compartments/simulation.h b/cpp/memilio/compartments/simulation.h index ff2b4fff70..e2179de30b 100644 --- a/cpp/memilio/compartments/simulation.h +++ b/cpp/memilio/compartments/simulation.h @@ -59,6 +59,20 @@ class Simulation { } + /** + * @brief Copy constructor. + * @param[in] other The simulation to copy. + */ + Simulation(const Simulation& other) + : m_integratorCore(other.m_integratorCore) + , m_model(std::make_unique(*other.m_model)) + , m_integrator(m_integratorCore) + , m_result(other.m_result) + , m_dt(other.m_dt) + { + m_integrator.set_integrator(m_integratorCore); + } + /** * @brief Set the integrator core used in the simulation. * @param[in] integrator A shared pointer to an object derived from IntegratorCore. diff --git a/cpp/memilio/data/analyze_result.h b/cpp/memilio/data/analyze_result.h index ecb501247a..ab9577cbdc 100644 --- a/cpp/memilio/data/analyze_result.h +++ b/cpp/memilio/data/analyze_result.h @@ -22,6 +22,7 @@ #include "memilio/utils/time_series.h" #include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" #include #include @@ -120,6 +121,18 @@ interpolate_simulation_result(const Graph, MigrationE return interpolated; } +template +std::vector> +interpolate_simulation_result(const GraphDetailed, MigrationEdge>& graph_result) +{ + std::vector> interpolated; + interpolated.reserve(graph_result.nodes().size()); + std::transform(graph_result.nodes().begin(), graph_result.nodes().end(), std::back_inserter(interpolated), + [](auto& n) { + return interpolate_simulation_result(n.property.get_result()); + }); + return interpolated; +} /** * Compute the distance between two SECIR simulation results. * The distance is the 2-norm of the element-wise difference of the two results. diff --git a/cpp/memilio/io/mobility_io.cpp b/cpp/memilio/io/mobility_io.cpp index 5b0169fa04..d60933447f 100644 --- a/cpp/memilio/io/mobility_io.cpp +++ b/cpp/memilio/io/mobility_io.cpp @@ -135,4 +135,93 @@ IOResult read_mobility_plain(const std::string& filename) return success(migration); } +IOResult read_duration_stay(const std::string& filename) +{ + BOOST_OUTCOME_TRY(num_lines, count_lines(filename)); + + if (num_lines == 0) { + return success(Eigen::MatrixXd(0, 0)); + } + + std::fstream file; + file.open(filename, std::ios::in); + if (!file.is_open()) { + return failure(StatusCode::FileNotFound, filename); + } + + Eigen::VectorXd duration(num_lines); + + try { + std::string tp; + int linenumber = 0; + while (getline(file, tp)) { + auto line = split(tp, ' '); + Eigen::Index i = static_cast(linenumber); + duration(i) = std::stod(line[0]); + linenumber++; + } + } + catch (std::runtime_error& ex) { + return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); + } + + return success(duration); +} + +IOResult>>> read_path_mobility(const std::string& filename) +{ + BOOST_OUTCOME_TRY(num_lines, count_lines(filename)); + + if (num_lines == 0) { + std::vector>> arr(0, std::vector>(0)); + return success(arr); + } + + std::fstream file; + file.open(filename, std::ios::in); + if (!file.is_open()) { + return failure(StatusCode::FileNotFound, filename); + } + + const int num_nodes = static_cast(std::sqrt(num_lines)); + std::vector>> arr(num_nodes, std::vector>(num_nodes)); + + try { + std::string tp; + while (getline(file, tp)) { + auto line = split(tp, ' '); + auto indx_x = std::stoi(line[0]); + auto indx_y = std::stoi(line[1]); + if (indx_x != indx_y) { + auto path = std::accumulate(line.begin() + 2, line.end(), std::string("")); + + // string -> vector of integers + std::vector path_vec; + + // Remove the square brackets and \r + path = path.substr(1, path.size() - 3); + std::stringstream ss(path); + std::string token; + + // get numbers and save them in path_vec + while (std::getline(ss, token, ',')) { + path_vec.push_back(std::stoi(token)); + } + + for (int number : path_vec) { + arr[indx_x][indx_y].push_back(number); + } + } + else { + arr[indx_x][indx_y].push_back(static_cast(indx_x)); + } + } + } + catch (std::runtime_error& ex) { + return failure(StatusCode::InvalidFileFormat, filename + ": " + ex.what()); + } + + return success(arr); +} + } // namespace mio diff --git a/cpp/memilio/io/mobility_io.h b/cpp/memilio/io/mobility_io.h index 7059951600..90df68aaa3 100644 --- a/cpp/memilio/io/mobility_io.h +++ b/cpp/memilio/io/mobility_io.h @@ -17,8 +17,8 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#ifndef READ_TWITTER_H -#define READ_TWITTER_H +#ifndef MIO_MOBILITY_IO_H +#define MIO_MOBILITY_IO_H #include "memilio/config.h" #include "memilio/math/eigen.h" @@ -63,6 +63,22 @@ IOResult read_mobility_formatted(const std::string& filename); */ IOResult read_mobility_plain(const std::string& filename); +/** + * @brief Reads txt file containing the duration of stay in each county. + Writes it into a Eigen vector of size N, where N is the number of local entites. + * @param filename name of file to be read + * @return IOResult containing the duration of stay in each local entity + */ +IOResult read_duration_stay(const std::string& filename); + +/** + * @brief For each edge we have the path from the start node to the end node. This functions reads the file and returns the path for each edge. + * + * @param filename Filename of the file containing the paths. + * @return IOResult>>> contains the paths for each edge. + */ +IOResult>>> read_path_mobility(const std::string& filename); + #ifdef MEMILIO_HAS_JSONCPP /** @@ -133,7 +149,8 @@ IOResult write_graph(const Graph& graph, const * @param read_edges boolean value that decides whether the edges of the graph should also be read in. */ template -IOResult> read_graph(const std::string& directory, int ioflags = IOF_None, bool read_edges = true) +IOResult> read_graph(const std::string& directory, int ioflags = IOF_None, + bool read_edges = true) { std::string abs_path; if (!file_exists(directory, abs_path)) { @@ -160,7 +177,7 @@ IOResult> read_graph(const std::string& direct } //read edges; nodes must already be available for that) - if(read_edges){ + if (read_edges) { for (auto inode = size_t(0); inode < graph.nodes().size(); ++inode) { //list of edges auto edge_filename = path_join(abs_path, "GraphEdges_node" + std::to_string(inode) + ".json"); @@ -177,7 +194,7 @@ IOResult> read_graph(const std::string& direct if (end_node_idx >= graph.nodes().size()) { log_error("EndNodeIndex not in range of number of graph nodes."); return failure(StatusCode::OutOfRange, - edge_filename + ", EndNodeIndex not in range of number of graph nodes."); + edge_filename + ", EndNodeIndex not in range of number of graph nodes."); } BOOST_OUTCOME_TRY(parameters, deserialize_json(e["Parameters"], Tag{}, ioflags)); graph.add_edge(start_node_idx, end_node_idx, parameters); @@ -192,4 +209,4 @@ IOResult> read_graph(const std::string& direct } // namespace mio -#endif // READ_TWITTER_H +#endif // MIO_MOBILITY_IO_H diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index 03975b8c2f..7a8328fb67 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -20,8 +20,10 @@ #ifndef MIO_MOBILITY_GRAPH_SIMULATION_H #define MIO_MOBILITY_GRAPH_SIMULATION_H +#include "memilio/math/euler.h" #include "memilio/mobility/graph.h" #include "memilio/utils/random_number_generator.h" +#include namespace mio { @@ -181,8 +183,8 @@ class GraphSimulationStochastic //perform transition Base::m_edge_func(event_edge.property, flat_index, - Base::m_graph.nodes()[event_edge.start_node_idx].property, - Base::m_graph.nodes()[event_edge.end_node_idx].property); + Base::m_graph.nodes()[event_edge.start_node_idx].property, + Base::m_graph.nodes()[event_edge.end_node_idx].property); //calculate new cumulative rate cumulative_rate = get_cumulative_transition_rate(); @@ -247,6 +249,484 @@ class GraphSimulationStochastic RandomNumberGenerator m_rng; }; +template > +class GraphSimulationDetailed +{ + +public: + using edge_function = edge_f; + using node_function = std::function; + + GraphSimulationDetailed(double t0, double dt, const GraphDetailed& g, const node_function& node_func, + const edge_function&& edge_func) + : m_dt(dt) + , m_t(double(t0)) + , m_graph(g) + , m_node_func(node_func) + , m_edge_func(edge_func) + { + } + + GraphSimulationDetailed(double t0, double dt, GraphDetailed&& g, const node_function& node_func, + const edge_function&& edge_func) + : m_dt(dt) + , m_t(double(t0)) + , m_graph(std::move(g)) + , m_node_func(node_func) + , m_edge_func(edge_func) + { + } + + inline double round_second_decimal(double x) + { + return std::round(x * 100.0) / 100.0; + } + + double get_t() const + { + return m_t; + } + + GraphDetailed& get_graph() & + { + return m_graph; + } + + const GraphDetailed& get_graph() const& + { + return m_graph; + } + + GraphDetailed&& get_graph() && + { + return std::move(m_graph); + } + + void advance(double t_max = 1.0) + { + ScalarType dt_first_mobility = + std::accumulate(m_graph.edges().begin(), m_graph.edges().end(), std::numeric_limits::max(), + [&](ScalarType current_min, const auto& e) { + auto traveltime_per_region = round_second_decimal(e.traveltime / e.path.size()); + if (traveltime_per_region < 0.01) + traveltime_per_region = 0.01; + auto start_mobility = + round_second_decimal(1 - 2 * (traveltime_per_region * e.path.size()) - + m_graph.nodes()[e.end_node_idx].stay_duration); + if (start_mobility < 0) { + start_mobility = 0.; + } + return std::min(current_min, start_mobility); + }); + + // set population to zero in mobility nodes before starting + for (auto& n : m_graph.nodes()) { + n.mobility.get_result().get_last_value().setZero(); + } + + auto min_dt = 0.01; + double t_begin = m_t - 1.; + + // calc schedule for each edge + precompute_schedule(); + while (m_t - epsilon < t_max) { + // auto start = std::chrono::system_clock::now(); + + t_begin += 1; + if (m_t + dt_first_mobility > t_max) { + dt_first_mobility = t_max - m_t; + for (auto& n : m_graph.nodes()) { + m_node_func(m_t, dt_first_mobility, n.property); + } + break; + } + + for (auto& node : m_graph.nodes()) { + node.mobility.get_simulation().set_integrator(std::make_shared()); + } + + size_t indx_schedule = 0; + while (t_begin + 1 > m_t + 1e-10) { + advance_edges(indx_schedule); + + // first we integrate the nodes in time. Afterwards the update on the edges is done. + // We start with the edges since the values for t0 are given. + advance_local_nodes(indx_schedule); + advance_mobility_nodes(indx_schedule); + + indx_schedule++; + m_t += min_dt; + } + // At the end of the day. we set each compartment zero for all mobility nodes since we have to estimate + // the state of the indivuals moving between different counties. + // Therefore there can be differences with the states given by the integrator used for the mobility node. + for (auto& n : m_graph.nodes()) { + n.mobility.get_result().get_last_value().setZero(); + } + } + } + +private: + GraphDetailed m_graph; + double m_t; + double m_dt; + node_function m_node_func; + edge_function m_edge_func; + + // describes the schedule for each edge, i.e. which node is visited at which time step + std::vector> schedule_edges; + + // defines if people from a edge are currently in a mobility node. This is necessary to determine the correct + // position in the schedule. Otherwise we could add a specific identifier in schedule_edges. + std::vector> mobility_schedule_edges; + + // describes the syncronization steps which are necessary for the integrator in the edges and nodes. + std::vector> mb_int_schedule; + std::vector> ln_int_schedule; + + // defines the edges on which mobility takes place for each time step + std::vector> edges_mobility; + + // same for the nodes + std::vector> nodes_mobility; + std::vector> nodes_mobility_m; + + // first mobility activites per edge + std::vector first_mobility; + + const double epsilon = std::numeric_limits::epsilon(); + + void precompute_schedule() + { + // print transi + const size_t timesteps = 100; + schedule_edges.reserve(m_graph.edges().size()); + mobility_schedule_edges.reserve(m_graph.edges().size()); + + for (auto& e : m_graph.edges()) { + // 100 since we round to second decimal + std::vector schedule(timesteps, 0.); + std::vector is_mobility_node(timesteps, false); + + double traveltime_per_region = round_second_decimal(e.traveltime / e.path.size()); + + if (traveltime_per_region < 0.01) + traveltime_per_region = 0.01; + + double start_mobility = round_second_decimal(0 + 1 - 2 * (traveltime_per_region * e.path.size()) - + m_graph.nodes()[e.end_node_idx].stay_duration); + if (start_mobility < 0.0) { + start_mobility = 0.; + } + + double arrive_at = start_mobility + traveltime_per_region * e.path.size(); + + std::fill(schedule.begin(), schedule.begin() + static_cast((start_mobility + epsilon) * 100), + e.start_node_idx); + + // all values true between start mob und arrive + std::fill(is_mobility_node.begin() + static_cast((start_mobility + epsilon) * 100), + is_mobility_node.begin() + static_cast((arrive_at + epsilon) * 100), true); + int count_true = std::count(is_mobility_node.begin(), is_mobility_node.end(), true); + size_t current_index = static_cast((start_mobility + epsilon) * 100); + for (size_t county : e.path) { + std::fill(schedule.begin() + current_index, + schedule.begin() + current_index + + static_cast((traveltime_per_region + epsilon) * 100), + county); + current_index += static_cast((traveltime_per_region + epsilon) * 100); + } + + // stay in destination county + std::fill(schedule.begin() + current_index, schedule.end() - count_true, e.path[e.path.size() - 1]); + + // revert trip for return. + std::fill(is_mobility_node.end() - count_true, is_mobility_node.end(), true); + + auto first_true = std::find(is_mobility_node.begin(), is_mobility_node.end(), true); + auto last_true = std::find(is_mobility_node.rbegin(), is_mobility_node.rend(), true); + + // If all values are false, we get an error + if (first_true != is_mobility_node.end() && last_true != is_mobility_node.rend()) { + int first_index = std::distance(is_mobility_node.begin(), first_true); + std::vector path_reversed(schedule.begin() + first_index, + schedule.begin() + first_index + count_true); + std::reverse(path_reversed.begin(), path_reversed.end()); + std::copy(path_reversed.begin(), path_reversed.end(), schedule.end() - count_true); + schedule_edges.push_back(schedule); + mobility_schedule_edges.push_back(is_mobility_node); + } + else { + log_error("Error in creating schedule."); + } + } + + // iterate over nodes + size_t indx_node = 0; + for (auto& n : m_graph.nodes()) { + // local node and automatical add zero, since we want to start at t=0 and start with integrating all nodes to the next dt + std::vector order_local_node = {0}; + std::vector indx_edges; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][0] == indx_node && !mobility_schedule_edges[indx_edge][0]) { + indx_edges.push_back(indx_edge); + } + } + + for (size_t indx_t = 1; indx_t < timesteps; ++indx_t) { + std::vector indx_edges_new; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][indx_t] == indx_node && !mobility_schedule_edges[indx_edge][indx_t]) { + indx_edges_new.push_back(indx_edge); + } + } + + if (indx_edges_new.size() != indx_edges.size() || + !std::equal(indx_edges.begin(), indx_edges.end(), indx_edges_new.begin())) { + order_local_node.push_back(indx_t); + indx_edges = indx_edges_new; + } + } + if (order_local_node[order_local_node.size() - 1] != timesteps - 1) + order_local_node.push_back(timesteps - 1); + ln_int_schedule.push_back(order_local_node); + + // mobility node + std::vector order_mobility_node; + std::vector indx_edges_mobility; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][0] == indx_node && mobility_schedule_edges[indx_edge][0]) { + indx_edges_mobility.push_back(indx_edge); + } + } + for (size_t indx_t = 1; indx_t < timesteps; ++indx_t) { + std::vector indx_edges_mobility_new; + for (size_t indx_edge = 0; indx_edge < schedule_edges.size(); ++indx_edge) { + if (schedule_edges[indx_edge][indx_t] == indx_node && mobility_schedule_edges[indx_edge][indx_t]) { + indx_edges_mobility_new.push_back(indx_edge); + } + } + + if (indx_edges_mobility_new.size() != indx_edges_mobility.size() || + !std::equal(indx_edges_mobility.begin(), indx_edges_mobility.end(), + indx_edges_mobility_new.begin())) { + order_mobility_node.push_back(indx_t); + indx_edges_mobility = indx_edges_mobility_new; + } + } + mb_int_schedule.push_back(order_mobility_node); + indx_node++; + } + + auto indx_edge = 0; + first_mobility.reserve(m_graph.edges().size()); + for (auto& e : m_graph.edges()) { + this->first_mobility[indx_edge] = std::distance( + mobility_schedule_edges[indx_edge].begin(), + std::find(mobility_schedule_edges[indx_edge].begin(), mobility_schedule_edges[indx_edge].end(), true)); + indx_edge++; + } + + edges_mobility.reserve(timesteps); + nodes_mobility.reserve(timesteps); + nodes_mobility_m.reserve(timesteps); + + // we handle indx_current = 0 separately since we want have added them always in the + // intregration schedule so this would lead to wrong results + // here we iterate over first mobility activities and add the edges indx when there is a start at t = 0 + std::vector temp_edges_mobility; + indx_edge = 0; + for (auto& start_time : first_mobility) { + if (start_time == 0) { + temp_edges_mobility.push_back(indx_edge); + } + indx_edge++; + } + edges_mobility.push_back(std::move(temp_edges_mobility)); + + // same for nodes + std::vector temp_nodes_mobility(m_graph.nodes().size()); + std::iota(temp_nodes_mobility.begin(), temp_nodes_mobility.end(), 0); + nodes_mobility.emplace_back(std::move(temp_nodes_mobility)); + + std::vector temp_nodes_mobility_m; + size_t node_indx = 0; + for (auto& n : m_graph.nodes()) { + if (std::binary_search(mb_int_schedule[node_indx].begin(), mb_int_schedule[node_indx].end(), 0)) { + temp_nodes_mobility_m.push_back(node_indx); + node_indx++; + } + } + nodes_mobility_m.push_back(temp_nodes_mobility_m); + + for (size_t indx_current = 1; indx_current < timesteps; ++indx_current) { + std::vector temp_edge_mobility; + indx_edge = 0; + for (auto& e : m_graph.edges()) { + auto current_node_indx = schedule_edges[indx_edge][indx_current]; + if (indx_current >= first_mobility[indx_edge]) { + if (mobility_schedule_edges[indx_edge][indx_current] && + std::binary_search(mb_int_schedule[current_node_indx].begin(), + mb_int_schedule[current_node_indx].end(), indx_current)) + temp_edge_mobility.push_back(indx_edge); + else if (!mobility_schedule_edges[indx_edge][indx_current] && + std::binary_search(ln_int_schedule[current_node_indx].begin(), + ln_int_schedule[current_node_indx].end(), indx_current)) + temp_edge_mobility.push_back(indx_edge); + } + indx_edge++; + } + edges_mobility.push_back(temp_edge_mobility); + + // reset temp_nodes_mobility + temp_nodes_mobility.clear(); + temp_nodes_mobility_m.clear(); + node_indx = 0; + for (auto& n : m_graph.nodes()) { + if (std::binary_search(ln_int_schedule[node_indx].begin(), ln_int_schedule[node_indx].end(), + indx_current)) { + temp_nodes_mobility.push_back(node_indx); + } + + if (std::binary_search(mb_int_schedule[node_indx].begin(), mb_int_schedule[node_indx].end(), + indx_current)) { + temp_nodes_mobility_m.push_back(node_indx); + } + node_indx++; + } + nodes_mobility.push_back(temp_nodes_mobility); + nodes_mobility_m.push_back(temp_nodes_mobility_m); + } + } + + void advance_edges(size_t indx_schedule) + { + for (const auto& edge_indx : edges_mobility[indx_schedule]) { + auto& e = m_graph.edges()[edge_indx]; + // first mobility activity + if (indx_schedule == first_mobility[edge_indx]) { + auto& node_from = m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + auto& node_to = m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility; + // m_edge_func(m_t, 0.0, e.property, node_from, node_to, 0); + } + // next mobility activity + else if (indx_schedule > first_mobility[edge_indx]) { + auto current_node_indx = schedule_edges[edge_indx][indx_schedule]; + bool in_mobility_node = mobility_schedule_edges[edge_indx][indx_schedule]; + + // determine dt, which is equal to the last integration/syncronization point in the current node + auto integrator_schedule_row = ln_int_schedule[current_node_indx]; + if (in_mobility_node) + integrator_schedule_row = mb_int_schedule[current_node_indx]; + // search the index of indx_schedule in the integrator schedule + const size_t indx_current = std::distance( + integrator_schedule_row.begin(), + std::lower_bound(integrator_schedule_row.begin(), integrator_schedule_row.end(), indx_schedule)); + + if (integrator_schedule_row[indx_current] != indx_schedule) + std::cout << "Error in schedule." + << "\n"; + + ScalarType dt_mobility; + if (indx_current == 0) { + dt_mobility = round_second_decimal(e.traveltime / e.path.size()); + if (dt_mobility < 0.01) + dt_mobility = 0.01; + } + else { + dt_mobility = + round_second_decimal((static_cast(integrator_schedule_row[indx_current]) - + static_cast(integrator_schedule_row[indx_current - 1])) / + 100 + + epsilon); + } + + // We have two cases. Either, we want to send the individuals to the next node, or we just want + // to update their state since a syncronization step is necessary in the current node. + if ((schedule_edges[edge_indx][indx_schedule] != schedule_edges[edge_indx][indx_schedule - 1]) || + (mobility_schedule_edges[edge_indx][indx_schedule] != + mobility_schedule_edges[edge_indx][indx_schedule - 1])) { + auto& node_from = mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + + auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + + if (indx_schedule < mobility_schedule_edges[edge_indx].size() - 1) { + // m_edge_func(m_t, dt_mobility, e.property, node_from, node_to, 1); + } + // else + // the last time step is handled differently since we have to close the timeseries + // m_edge_func(m_t, dt_mobility, e.property, node_from, + // m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property, 2); + } + else { + auto& node_from = mobility_schedule_edges[edge_indx][indx_schedule - 1] + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule - 1]].property; + + auto& node_to = mobility_schedule_edges[edge_indx][indx_schedule] + ? m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].mobility + : m_graph.nodes()[schedule_edges[edge_indx][indx_schedule]].property; + + assert(node_from.get_result().get_last_value() == node_to.get_result().get_last_value()); + // m_edge_func(m_t, dt_mobility, e.property, node_from, node_to, 3); + } + } + } + } + + void advance_local_nodes(size_t indx_schedule) + { + for (const auto& n_indx : nodes_mobility[indx_schedule]) { + auto& n = m_graph.nodes()[n_indx]; + const size_t indx_current = std::distance( + ln_int_schedule[n_indx].begin(), + std::lower_bound(ln_int_schedule[n_indx].begin(), ln_int_schedule[n_indx].end(), indx_schedule)); + + const size_t val_next = + (indx_current == ln_int_schedule[n_indx].size() - 1) ? 100 : ln_int_schedule[n_indx][indx_current + 1]; + const ScalarType next_dt = + round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); + m_node_func(m_t, next_dt, n.property); + } + } + + void advance_mobility_nodes(size_t indx_schedule) + { + for (const size_t& n_indx : nodes_mobility_m[indx_schedule]) { + auto& n = m_graph.nodes()[n_indx]; + const size_t indx_current = std::distance( + mb_int_schedule[n_indx].begin(), + std::lower_bound(mb_int_schedule[n_indx].begin(), mb_int_schedule[n_indx].end(), indx_schedule)); + const size_t val_next = + (indx_current == mb_int_schedule[n_indx].size() - 1) ? 100 : mb_int_schedule[n_indx][indx_current + 1]; + const ScalarType next_dt = + round_second_decimal((static_cast(val_next) - indx_schedule) / 100 + epsilon); + + // get all time points from the last integration step + auto& last_time_point = n.mobility.get_result().get_time(n.mobility.get_result().get_num_time_points() - 1); + // wenn last_time_point nicht innerhalb eines intervalls von 1-e10 von t liegt, dann setzte den letzten Zeitpunkt auf m_t + if (std::fabs(last_time_point - m_t) > 1e-10) { + n.mobility.get_result().get_last_time() = m_t; + } + m_node_func(m_t, next_dt, n.mobility); + Eigen::Index indx_min; + while (n.mobility.get_result().get_last_value().minCoeff(&indx_min) < -1e-7) { + Eigen::Index indx_max; + n.mobility.get_result().get_last_value().maxCoeff(&indx_max); + n.mobility.get_result().get_last_value()[indx_max] -= + n.mobility.get_result().get_last_value()[indx_min]; + n.mobility.get_result().get_last_value()[indx_min] = 0; + } + } + } +}; + template auto make_graph_sim(double t0, double dt, Graph&& g, NodeF&& node_func, EdgeF&& edge_func) { @@ -261,5 +741,12 @@ auto make_graph_sim_stochastic(double t0, double dt, Graph&& g, NodeF&& node_fun t0, dt, std::forward(g), std::forward(node_func), std::forward(edge_func)); } +template +auto make_graph_sim_detailed(double t0, double dt, GraphDetailed&& g, NodeF&& node_func, EdgeF&& edge_func) +{ + return GraphSimulationDetailed>( + t0, dt, std::forward(g), std::forward(node_func), std::forward(edge_func)); +} + } // namespace mio #endif //MIO_MOBILITY_GRAPH_SIMULATION_H diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp b/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp new file mode 100644 index 0000000000..b5c7d1a178 --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.cpp @@ -0,0 +1,25 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Henrik Zunker +* +* 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_detailed.h" + +namespace mio +{ + +} // namespace mio diff --git a/cpp/memilio/mobility/metapopulation_mobility_detailed.h b/cpp/memilio/mobility/metapopulation_mobility_detailed.h new file mode 100644 index 0000000000..af21ebcb9a --- /dev/null +++ b/cpp/memilio/mobility/metapopulation_mobility_detailed.h @@ -0,0 +1,688 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Henrik Zunker +* +* 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_DETAILED_H +#define METAPOPULATION_MOBILITY_DETAILED_H + +#include "memilio/compartments/parameter_studies.h" +#include "memilio/epidemiology/simulation_day.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/time_series.h" +#include "memilio/math/eigen.h" +#include "memilio/math/eigen_util.h" +#include "memilio/utils/metaprogramming.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "memilio/math/euler.h" +#include "memilio/epidemiology/contact_matrix.h" +#include "memilio/epidemiology/dynamic_npis.h" +#include "memilio/compartments/simulation.h" +#include "memilio/utils/date.h" +#include "memilio/mobility/graph.h" +#include "memilio/io/mobility_io.h" + +#include "boost/filesystem.hpp" + +#include +#include + +namespace mio +{ + +template +struct NodeDetailed : Node { + template + NodeDetailed(int node_id, Args&&... args) + : Node(node_id, std::forward(args)...) + , stay_duration(0.5) + , mobility(std::forward(args)...) + { + } + + template + NodeDetailed(int node_id, double duration, Model property_arg, Model mobility_arg, double m_t0, + double m_dt_integration) + : Node(node_id, property_arg, m_t0, m_dt_integration) + , stay_duration(duration) + , mobility(mobility_arg, m_t0, m_dt_integration) + { + } + + template + NodeDetailed(int node_id, double duration, Args&&... args) + : Node(node_id, std::forward(args)...) + , stay_duration(duration) + , mobility(std::forward(args)...) + { + } + + NodeDetailed(int node_id, double duration, NodePropertyT property_arg, NodePropertyT mobility_pt_arg) + : Node(node_id, property_arg) + , stay_duration(duration) + , mobility(mobility_pt_arg) + { + } + + double stay_duration; + NodePropertyT mobility; +}; + +template +struct EdgeDetailed : Edge { + template + EdgeDetailed(size_t start, size_t end, Args&&... args) + : Edge(start, end, std::forward(args)...) + , traveltime(0.) + , path{static_cast(start), static_cast(end)} + { + } + + template + EdgeDetailed(size_t start, size_t end, double t_travel, Args&&... args) + : Edge(start, end, std::forward(args)...) + , traveltime(t_travel) + , path{static_cast(start), static_cast(end)} + { + } + + template + EdgeDetailed(size_t start, size_t end, double t_travel, std::vector path_mobility, Args&&... args) + : Edge(start, end, std::forward(args)...) + , traveltime(t_travel) + , path(path_mobility) + { + } + double traveltime; + std::vector path; +}; + +template +class GraphDetailed : public Graph +{ +public: + using Graph::Graph; + + /** + * @brief range of nodes + */ + auto nodes() + { + return make_range(begin(m_nodes), end(m_nodes)); + } + + /** + * @brief range of nodes + */ + auto nodes() const + { + return make_range(begin(m_nodes), end(m_nodes)); + }; + + /** + * @brief range of edges + */ + auto edges() + { + return make_range(begin(m_edges), end(m_edges)); + } + + /** + * @brief range of edges + */ + auto edges() const + { + return make_range(begin(m_edges), end(m_edges)); + } + + /** + * @brief range of edges going out from a specific node + */ + auto out_edges(size_t node_idx) + { + return out_edges(begin(m_edges), end(m_edges), node_idx); + } + + /** + * @brief range of edges going out from a specific node + */ + auto out_edges(size_t node_idx) const + { + return out_edges(begin(m_edges), end(m_edges), node_idx); + } + + template + NodeDetailed& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2) + { + m_nodes.emplace_back(id, duration_stay, model1, model2); + return m_nodes.back(); + } + + template + NodeDetailed& add_node(int id, double duration_stay, ModelType& model1, ModelType& model2, + double m_t0, double m_dt_integration) + { + m_nodes.emplace_back(id, duration_stay, model1, model2, m_t0, m_dt_integration); + return m_nodes.back(); + } + + EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, + mio::MigrationParameters& args) + { + assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); + return *insert_sorted_replace(m_edges, + EdgeDetailed(start_node_idx, end_node_idx, traveltime, args), + [](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; + }); + } + + template + EdgeDetailed& add_edge(size_t start_node_idx, size_t end_node_idx, double traveltime, + std::vector path, Args&&... args) + { + assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx); + return *insert_sorted_replace( + m_edges, + EdgeDetailed(start_node_idx, end_node_idx, traveltime, path, std::forward(args)...), + [](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; + }); + } + +private: + std::vector> m_nodes; + std::vector> m_edges; +}; + +/** + * @brief Sets the graph nodes for counties or districts. + * Reads the node ids which could refer to districts or counties and the epidemiological + * data from json files and creates one node for each id. Every node contains a model. + * @param[in] params Model Parameters that are used for every node. + * @param[in] start_date Start date for which the data should be read. + * @param[in] end_data End date for which the data should be read. + * @param[in] data_dir Directory that contains the data files. + * @param[in] population_data_path Path to json file containing the population data. + * @param[in] stay_times_data_path Path to txt file containing the stay times for the considered local entities. + * @param[in] is_node_for_county Specifies whether the node ids should be county ids (true) or district ids (false). + * @param[in, out] params_graph Graph whose nodes are set by the function. + * @param[in] read_func Function that reads input data for german counties and sets Model compartments. + * @param[in] node_func Function that returns the county ids. + * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county. + * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting. + * @param[in] tnt_capacity_factor Factor for test and trace capacity. + * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation. + * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files. + */ +template +IOResult set_nodes_detailed(const Parameters& params, Date start_date, Date end_date, const fs::path& data_dir, + const std::string& population_data_path, const std::string& stay_times_data_path, + bool is_node_for_county, GraphDetailed& params_graph, + ReadFunction&& read_func, NodeIdFunction&& node_func, + const std::vector& scaling_factor_inf, double scaling_factor_icu, + double tnt_capacity_factor, int num_days = 0, bool export_time_series = false) +{ + + BOOST_OUTCOME_TRY(duration_stay, mio::read_duration_stay(stay_times_data_path)); + BOOST_OUTCOME_TRY(node_ids, node_func(population_data_path, is_node_for_county)); + + std::vector nodes(node_ids.size(), Model(int(size_t(params.get_num_groups())))); + + for (auto& node : nodes) { + node.parameters = params; + } + BOOST_OUTCOME_TRY(read_func(nodes, start_date, node_ids, scaling_factor_inf, scaling_factor_icu, data_dir.string(), + num_days, export_time_series)); + + for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) { + + auto tnt_capacity = nodes[node_idx].populations.get_total() * tnt_capacity_factor; + + //local parameters + auto& tnt_value = nodes[node_idx].parameters.template get(); + tnt_value = mio::UncertainValue(0.5 * (1.2 * tnt_capacity + 0.8 * tnt_capacity)); + tnt_value.set_distribution(mio::ParameterDistributionUniform(0.8 * tnt_capacity, 1.2 * tnt_capacity)); + + //holiday periods + auto id = int(mio::regions::CountyId(node_ids[node_idx])); + auto holiday_periods = mio::regions::get_holidays(mio::regions::get_state_id(id), start_date, end_date); + auto& contacts = nodes[node_idx].parameters.template get(); + contacts.get_school_holidays() = + std::vector>(holiday_periods.size()); + std::transform( + holiday_periods.begin(), holiday_periods.end(), contacts.get_school_holidays().begin(), [=](auto& period) { + return std::make_pair(mio::SimulationTime(mio::get_offset_in_days(period.first, start_date)), + mio::SimulationTime(mio::get_offset_in_days(period.second, start_date))); + }); + + //uncertainty in populations + for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + for (auto j = mio::Index(0); j < Model::Compartments::Count; ++j) { + auto& compartment_value = nodes[node_idx].populations[{i, j}]; + compartment_value = + mio::UncertainValue(0.5 * (1.1 * double(compartment_value) + 0.9 * double(compartment_value))); + compartment_value.set_distribution(mio::ParameterDistributionUniform(0.9 * double(compartment_value), + 1.1 * double(compartment_value))); + } + } + + // Add mobility node + auto mobility = nodes[node_idx]; + mobility.populations.set_total(0); + + params_graph.add_node(node_ids[node_idx], duration_stay((Eigen::Index)node_idx), nodes[node_idx], mobility); + } + return success(); +} + +/** + * @brief Sets the graph edges. + * Reads the commuting matrices, travel times and paths from data and creates one edge for each pair of nodes. + * @param[in] travel_times_path Path to txt file containing the travel times between counties. + * @param[in] mobility_data_path Path to txt file containing the commuting matrices. + * @param[in] travelpath_path Path to txt file containing the paths between counties. + * @param[in, out] params_graph Graph whose nodes are set by the function. + * @param[in] migrating_compartments Compartments that commute. + * @param[in] contact_locations_size Number of contact locations. + * @param[in] commuting_weights Vector with a commuting weight for every AgeGroup. + */ +template +IOResult +set_edges_detailed(const std::string& travel_times_path, const std::string mobility_data_path, + const std::string& travelpath_path, GraphDetailed& params_graph, + std::initializer_list& migrating_compartments, size_t contact_locations_size, + std::vector commuting_weights = std::vector{}, + ScalarType theshold_edges = 4e-5) +{ + BOOST_OUTCOME_TRY(mobility_data_commuter, mio::read_mobility_plain(mobility_data_path)); + BOOST_OUTCOME_TRY(travel_times, mio::read_mobility_plain(travel_times_path)); + BOOST_OUTCOME_TRY(path_mobility, mio::read_path_mobility(travelpath_path)); + + for (size_t county_idx_i = 0; county_idx_i < params_graph.nodes().size(); ++county_idx_i) { + for (size_t county_idx_j = 0; county_idx_j < params_graph.nodes().size(); ++county_idx_j) { + auto& populations = params_graph.nodes()[county_idx_i].property.populations; + + // mobility coefficients have the same number of components as the contact matrices. + // so that the same NPIs/dampings can be used for both (e.g. more home office => fewer commuters) + auto mobility_coeffs = MigrationCoefficientGroup(contact_locations_size, populations.numel()); + auto num_age_groups = (size_t)params_graph.nodes()[county_idx_i].property.parameters.get_num_groups(); + commuting_weights = + (commuting_weights.size() == 0 ? std::vector(num_age_groups, 1.0) : commuting_weights); + + //commuters + auto working_population = 0.0; + auto min_commuter_age = mio::AgeGroup(2); + auto max_commuter_age = mio::AgeGroup(4); //this group is partially retired, only partially commutes + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + working_population += populations.get_group_total(age) * commuting_weights[size_t(age)]; + } + auto commuter_coeff_ij = mobility_data_commuter(county_idx_i, county_idx_j) / working_population; + for (auto age = min_commuter_age; age <= max_commuter_age; ++age) { + for (auto compartment : migrating_compartments) { + auto coeff_index = populations.get_flat_index({age, compartment}); + mobility_coeffs[size_t(ContactLocation::Work)].get_baseline()[coeff_index] = + commuter_coeff_ij * commuting_weights[size_t(age)]; + } + } + + auto path = path_mobility[county_idx_i][county_idx_j]; + if (static_cast(path[0]) != county_idx_i || + static_cast(path[path.size() - 1]) != county_idx_j) + std::cout << "Wrong Path for edge " << county_idx_i << " " << county_idx_j << "\n"; + + //only add edges with mobility above thresholds for performance + if (commuter_coeff_ij > theshold_edges) { + params_graph.add_edge(county_idx_i, county_idx_j, travel_times(county_idx_i, county_idx_j), + path_mobility[county_idx_i][county_idx_j], std::move(mobility_coeffs)); + } + } + } + + return success(); +} + +// class MigrationEdgeDetailed : public MigrationEdge +// { +// public: +// template +// void apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to, int mode); +// }; + +// /** +// * edge functor for migration simulation. +// * @see MigrationEdge::apply_migration +// */ +// template +// void apply_migration(double t, double dt, MigrationEdgeDetailed& migrationEdgeDetailed, SimulationNode& node_from, +// SimulationNode& node_to, int mode) +// { +// migrationEdgeDetailed.apply_migration(t, dt, node_from, node_to, mode); +// } + +template +class ParameterStudyDetailed : public ParameterStudy +{ +public: + /** + * The type of simulation of a single node of the graph. + */ + using Simulation = S; + /** + * The Graph type that stores the parametes of the simulation. + * This is the input of ParameterStudies. + */ + using ParametersGraphDetailed = mio::GraphDetailed; + /** + * The Graph type that stores simulations and their results of each run. + * This is the output of ParameterStudies for each run. + */ + using SimulationGraphDetailed = mio::GraphDetailed, mio::MigrationEdge>; + + /** + * create study for graph of compartment models. + * @param graph graph of parameters + * @param t0 start time of simulations + * @param tmax end time of simulations + * @param graph_sim_dt time step of graph simulation + * @param num_runs number of runs + */ + ParameterStudyDetailed(const ParametersGraphDetailed& graph, double t0, double tmax, double graph_sim_dt, + size_t num_runs) + : ParameterStudy(graph, t0, tmax, graph_sim_dt, num_runs) + , m_graph(graph) + , m_num_runs(num_runs) + , m_t0{t0} + , m_tmax{tmax} + , m_dt_graph_sim(graph_sim_dt) + { + } + + /* + * @brief Carry out all simulations in the parameter study. + * Save memory and enable more runs by immediately processing and/or discarding the result. + * The result processing function is called when a run is finished. It receives the result of the run + * (a SimulationGraphDetailed object) and an ordered index. The values returned by the result processing function + * are gathered and returned as a list. + * This function is parallelized if memilio is configured with MEMILIO_ENABLE_MPI. + * The MPI processes each contribute a share of the runs. The sample function and result processing function + * are called in the same process that performs the run. The results returned by the result processing function are + * gathered at the root process and returned as a list by the root in the same order as if the programm + * were running sequentially. Processes other than the root return an empty list. + * @param sample_graph Function that receives the ParametersGraph and returns a sampled copy. + * @param result_processing_function Processing function for simulation results, e.g., output function. + * @returns At the root process, a list of values per run that have been returned from the result processing function. + * At all other processes, an empty list. + * @tparam SampleGraphFunction Callable type, accepts instance of ParametersGraph. + * @tparam HandleSimulationResultFunction Callable type, accepts instance of SimulationGraphDetailed and an index of type size_t. + */ + template + std::vector> + run(SampleGraphFunction sample_graph, HandleSimulationResultFunction result_processing_function) + { + int num_procs, rank; +#ifdef MEMILIO_ENABLE_MPI + MPI_Comm_size(mpi::get_world(), &num_procs); + MPI_Comm_rank(mpi::get_world(), &rank); +#else + num_procs = 1; + rank = 0; +#endif + + //The ParameterDistributions used for sampling parameters use thread_local_rng() + //So we set our own RNG to be used. + //Assume that sampling uses the thread_local_rng() and isn't multithreaded + this->m_rng.synchronize(); + thread_local_rng() = this->m_rng; + + auto run_distribution = this->distribute_runs(m_num_runs, num_procs); + auto start_run_idx = + std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0)); + auto end_run_idx = start_run_idx + run_distribution[size_t(rank)]; + + std::vector> + ensemble_result; + ensemble_result.reserve(m_num_runs); + + for (size_t run_idx = start_run_idx; run_idx < end_run_idx; run_idx++) { + log(LogLevel::info, "ParameterStudies: run {}", run_idx); + + //prepare rng for this run by setting the counter to the right offset + //Add the old counter so that this call of run() produces different results + //from the previous call + auto run_rng_counter = + this->m_rng.get_counter() + + rng_totalsequence_counter(static_cast(run_idx), Counter(0)); + thread_local_rng().set_counter(run_rng_counter); + + //sample + auto sim = create_sampled_sim(sample_graph); + log(LogLevel::info, "ParameterStudies: Generated {} random numbers.", + (thread_local_rng().get_counter() - run_rng_counter).get()); + + //perform run + sim.advance(m_tmax); + + //handle result and store + ensemble_result.emplace_back(result_processing_function(std::move(sim).get_graph(), run_idx)); + } + + //Set the counter of our RNG so that future calls of run() produce different parameters. + this->m_rng.set_counter(this->m_rng.get_counter() + + rng_totalsequence_counter(m_num_runs, Counter(0))); + +#ifdef MEMILIO_ENABLE_MPI + //gather results + if (rank == 0) { + for (int src_rank = 1; src_rank < num_procs; ++src_rank) { + int bytes_size; + MPI_Recv(&bytes_size, 1, MPI_INT, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE); + ByteStream bytes(bytes_size); + MPI_Recv(bytes.data(), bytes.data_size(), MPI_BYTE, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE); + + auto src_ensemble_results = deserialize_binary(bytes, Tag{}); + if (!src_ensemble_results) { + log_error("Error receiving ensemble results from rank {}.", src_rank); + } + std::copy(src_ensemble_results.value().begin(), src_ensemble_results.value().end(), + std::back_inserter(ensemble_result)); + } + } + else { + auto bytes = serialize_binary(ensemble_result); + auto bytes_size = int(bytes.data_size()); + MPI_Send(&bytes_size, 1, MPI_INT, 0, 0, mpi::get_world()); + MPI_Send(bytes.data(), bytes.data_size(), MPI_BYTE, 0, 0, mpi::get_world()); + ensemble_result.clear(); //only return root process + } +#endif + + return ensemble_result; + } + +private: + //sample parameters and create simulation + template + mio::GraphSimulationDetailed create_sampled_sim(SampleGraphFunction sample_graph) + { + SimulationGraphDetailed sim_graph; + + auto sampled_graph = sample_graph(m_graph); + for (auto&& node : sampled_graph.nodes()) { + sim_graph.add_node(node.id, node.stay_duration, node.property, node.mobility, m_t0, this->m_dt_integration); + } + for (auto&& edge : sampled_graph.edges()) { + sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.property); + } + return make_migration_sim(m_t0, m_dt_graph_sim, std::move(sim_graph)); + } + +private: + // Stores Graph with the names and ranges of all parameters + ParametersGraphDetailed m_graph; + size_t m_num_runs; + double m_t0; + double m_tmax; + double m_dt_graph_sim; +}; + +/** + * @brief number of migrated people when they return according to the model. + * E.g. during the time in the other node, some people who left as susceptible will return exposed. + * Implemented for general compartmentmodel simulations, overload for your custom model if necessary + * so that it can be found with argument-dependent lookup, i.e. in the same namespace as the model. + * @param migrated number of people that migrated as input, number of people that return as output + * @param sim Simulation that is used for the migration + * @param integrator Integrator that is used for the estimation. Has to be a one-step scheme. + * @param total total population in the node that the people migrated to. + * @param t time of migration + * @param dt timestep + */ +template ::value>> +void update_status_migrated(Eigen::Ref::Vector> migrated, const Sim& sim, IntegratorCore& integrator, + Eigen::Ref::Vector> total, double t, double dt) +{ + auto y0 = migrated.eval(); + auto y1 = migrated; + integrator.step( + [&](auto&& y, auto&& t_, auto&& dydt) { + sim.get_model().get_derivatives(total, y, t_, dydt); + }, + y0, t, dt, y1); +} + +template +using Vector = Eigen::Matrix; + +/* + * move migrated people from one node to another. + * @param migrated number of people that migrated + * @param results_from current state of node that people migrated from + * @param results_to current state of node that people migrated to +*/ +template +void move_migrated(Eigen::Ref> migrated, Eigen::Ref> results_from, + Eigen::Ref> results_to) +{ + // The performance of the low order integrator often fails for small compartments i.e. compartments are sometimes negative. + // We handel this by checking for values below 0 in m_migrated and add them to the highest value + Eigen::Index max_index_migrated; + for (size_t i = 0; i < migrated.size(); ++i) { + migrated.maxCoeff(&max_index_migrated); + if (migrated(i) < 0) { + migrated(max_index_migrated) -= migrated(i); + migrated(i) = 0; + } + if (results_from(i) < migrated(i)) { + migrated(max_index_migrated) -= (migrated(i) - results_from(i)); + migrated(i) = results_from(i); + } + results_from(i) -= migrated(i); + results_to(i) += migrated(i); + } +} + +// template +// void MigrationEdge::apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to, +// int mode) // +// { + +// if (mode == 0) { +// //normal daily migration +// m_migrated.add_time_point( +// t, (node_from.get_last_state().array() * m_parameters.get_coefficients().get_matrix_at(t).array() * +// get_migration_factors(node_from, t, node_from.get_last_state()).array()) +// .matrix()); +// m_return_times.add_time_point(t + dt); +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_to.get_result().get_last_value()); +// } +// // change county of migrated +// else if (mode == 1) { +// // update status of migrated before moving to next county +// IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); +// update_status_migrated(m_migrated.get_last_value(), node_from.get_simulation(), integrator_node, +// node_from.get_result().get_last_value(), t, dt); +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_to.get_result().get_last_value()); +// } +// // option for last time point to remove time points +// else if (mode == 2) { +// Eigen::Index idx = m_return_times.get_num_time_points() - 1; +// IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); +// update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, +// node_from.get_result().get_last_value(), t, dt); + +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_to.get_result().get_last_value()); + +// for (Eigen::Index i = m_return_times.get_num_time_points() - 1; i >= 0; --i) { +// if (m_return_times.get_time(i) <= t) { +// m_migrated.remove_time_point(i); +// m_return_times.remove_time_point(i); +// } +// } +// } +// // just update status of migrated +// else if (mode == 3) { +// Eigen::Index idx = m_return_times.get_num_time_points() - 1; +// IntegratorCore& integrator_node = node_from.get_simulation().get_integrator(); +// update_status_migrated(m_migrated[idx], node_from.get_simulation(), integrator_node, +// node_from.get_result().get_last_value(), t, dt); + +// move_migrated(m_migrated.get_last_value(), node_from.get_result().get_last_value(), +// node_from.get_result().get_last_value()); +// } +// else { +// std::cout << "Invalid input mode. Should be 0 or 1." +// << "\n"; +// } +// } + +/** + * create a migration simulation. + * After every second time step, for each edge a portion of the population corresponding to the coefficients of the edge + * moves from one node to the other. In the next timestep, the migrated population return 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 migrations + * @param graph set up for migration simulation + * @{ + */ +template +GraphSimulationDetailed, MigrationEdge>> +make_migration_sim(double t0, double dt, const GraphDetailed, MigrationEdge>& graph) +{ + return make_graph_sim_detailed(t0, dt, graph, &evolve_model, &apply_migration); +} + +template +GraphSimulationDetailed, MigrationEdge>> +make_migration_sim(double t0, double dt, GraphDetailed, MigrationEdge>&& graph) +{ + return make_graph_sim_detailed(t0, dt, std::move(graph), &evolve_model, &apply_migration); +} + +} // namespace mio + +#endif //METAPOPULATION_MOBILITY_DETAILED_H diff --git a/cpp/memilio/mobility/metapopulation_mobility_instant.h b/cpp/memilio/mobility/metapopulation_mobility_instant.h index 5178315036..6d7e6ecf5d 100644 --- a/cpp/memilio/mobility/metapopulation_mobility_instant.h +++ b/cpp/memilio/mobility/metapopulation_mobility_instant.h @@ -54,6 +54,20 @@ class SimulationNode { } + SimulationNode(const Sim& property_arg, double t0, double dt_integration) + : m_simulation(property_arg, t0, dt_integration) + , m_last_state(m_simulation.get_result().get_last_value()) + , m_t0(m_simulation.get_result().get_last_time()) + { + } + + SimulationNode(const SimulationNode& other) + : m_simulation(other.m_simulation) + , m_last_state(other.m_last_state) + , m_t0(other.m_t0) + { + } + /** * get the result of the simulation in this node. * @{ @@ -292,10 +306,12 @@ class MigrationEdge template void apply_migration(double t, double dt, SimulationNode& node_from, SimulationNode& node_to); -private: +protected: MigrationParameters m_parameters; TimeSeries m_migrated; TimeSeries m_return_times; + +private: bool m_return_migrated; double m_t_last_dynamic_npi_check = -std::numeric_limits::infinity(); std::pair m_dynamic_npi = {-std::numeric_limits::max(), SimulationTime(0)}; @@ -389,8 +405,8 @@ auto get_migration_factors(const SimulationNode& node, double t, const Eige * detect a get_migration_factors function for the Model type. */ template -using test_commuters_expr_t = decltype(test_commuters( - std::declval(), std::declval&>(), std::declval())); +using test_commuters_expr_t = decltype( + test_commuters(std::declval(), std::declval&>(), std::declval())); /** * Test persons when migrating from their source node. diff --git a/cpp/models/ode_secirvvs/parameter_space.cpp b/cpp/models/ode_secirvvs/parameter_space.cpp index 1748c307ae..9e26a39715 100644 --- a/cpp/models/ode_secirvvs/parameter_space.cpp +++ b/cpp/models/ode_secirvvs/parameter_space.cpp @@ -23,6 +23,8 @@ #include "ode_secirvvs/infection_state.h" #include "ode_secirvvs/model.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" + namespace mio { namespace osecirvvs @@ -191,5 +193,53 @@ Graph draw_sample(Graph& return sampled_graph; } +GraphDetailed draw_sample(GraphDetailed& graph) +{ + GraphDetailed sampled_graph; + + //sample global parameters + auto& shared_params_model = graph.nodes()[0].property; + draw_sample_infection(shared_params_model); + auto& shared_contacts = shared_params_model.parameters.template get(); + shared_contacts.draw_sample_dampings(); + auto& shared_dynamic_npis = shared_params_model.parameters.template get(); + shared_dynamic_npis.draw_sample(); + + for (auto& params_node : graph.nodes()) { + auto& node_model = params_node.property; + + //sample local parameters + draw_sample_demographics(params_node.property); + + //copy global parameters + //save demographic parameters so they aren't overwritten + auto local_icu_capacity = node_model.parameters.template get(); + auto local_tnt_capacity = node_model.parameters.template get(); + auto local_holidays = node_model.parameters.template get().get_school_holidays(); + auto local_daily_v1 = node_model.parameters.template get(); + auto local_daily_v2 = node_model.parameters.template get(); + node_model.parameters = shared_params_model.parameters; + node_model.parameters.template get() = local_icu_capacity; + node_model.parameters.template get() = local_tnt_capacity; + node_model.parameters.template get().get_school_holidays() = local_holidays; + node_model.parameters.template get() = local_daily_v1; + node_model.parameters.template get() = local_daily_v2; + + node_model.parameters.template get().make_matrix(); + node_model.apply_constraints(); + + sampled_graph.add_node(params_node.id, params_node.stay_duration, node_model, params_node.mobility); + } + + for (auto& edge : graph.edges()) { + auto edge_params = edge.property; + //no dynamic NPIs + //TODO: add switch to optionally enable dynamic NPIs to edges + sampled_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.traveltime, edge.path, edge_params); + } + + return sampled_graph; +} + } // namespace osecirvvs } // namespace mio diff --git a/cpp/models/ode_secirvvs/parameter_space.h b/cpp/models/ode_secirvvs/parameter_space.h index a883cc482a..1943806f45 100644 --- a/cpp/models/ode_secirvvs/parameter_space.h +++ b/cpp/models/ode_secirvvs/parameter_space.h @@ -21,6 +21,8 @@ #define ODESECIRVVS_PARAMETER_SPACE_H #include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" + #include "memilio/utils/memory.h" #include "memilio/utils/logging.h" #include "memilio/utils/parameter_distributions.h" @@ -63,6 +65,8 @@ void draw_sample(Model& model); */ Graph draw_sample(Graph& graph, bool variant_high); +GraphDetailed draw_sample(GraphDetailed& graph); + } // namespace osecirvvs } // namespace mio diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index c19b9fcec2..f31153406b 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -49,6 +49,7 @@ set(TESTSOURCES test_io_framework.cpp test_binary_serializer.cpp test_compartmentsimulation.cpp + test_mobility_detailed.cpp test_mobility_io.cpp test_transform_iterator.cpp test_metaprogramming.cpp @@ -68,17 +69,17 @@ set(TESTSOURCES ) if(MEMILIO_HAS_JSONCPP) -set(TESTSOURCES ${TESTSOURCES} -test_json_serializer.cpp -test_epi_data_io.cpp -) + set(TESTSOURCES ${TESTSOURCES} + test_json_serializer.cpp + test_epi_data_io.cpp + ) endif() if(MEMILIO_HAS_JSONCPP AND MEMILIO_HAS_HDF5) -set(TESTSOURCES ${TESTSOURCES} -test_save_parameters.cpp -test_save_results.cpp -) + set(TESTSOURCES ${TESTSOURCES} + test_save_parameters.cpp + test_save_results.cpp + ) endif() add_executable(memilio-test ${TESTSOURCES}) diff --git a/cpp/tests/test_mobility_detailed.cpp b/cpp/tests/test_mobility_detailed.cpp new file mode 100644 index 0000000000..b7fa254fa9 --- /dev/null +++ b/cpp/tests/test_mobility_detailed.cpp @@ -0,0 +1,47 @@ + +#include "gtest/gtest.h" +#include "gmock/gmock.h" +#include "memilio/mobility/graph_simulation.h" +#include "memilio/mobility/metapopulation_mobility_detailed.h" +#include "memilio/compartments/simulation.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "memilio/utils/compiler_diagnostics.h" +#include "ode_secir/model.h" +#include "ode_secir/parameters.h" + +namespace mio +{ +mio::IOResult read_duration_stay(const std::string& filename) +{ + mio::unused(filename); + return mio::success(Eigen::MatrixXd(0, 0)); +} + +TEST(TestMobilityDetailed, AddNodeTest) +{ + mio::osecir::Model model(1); + GraphDetailed>, int> sim_graph; + + auto& node = sim_graph.add_node(1, 0.5, model, model, 0.0, 0.1); + EXPECT_EQ(node.stay_duration, 0.5); +} + +TEST(TestMobilityDetailed, AddNodeWithTimeIntegrationTest) +{ + mio::osecir::Model model(1); + GraphDetailed model_graph; + auto& node = model_graph.add_node(1, 0.5, model, model); + EXPECT_EQ(node.stay_duration, 0.5); +} + +TEST(TestMobilityDetailed, AddEdgeTest) +{ + mio::osecir::Model model(1); + GraphDetailed model_graph; + model_graph.add_node(0, 0.5, model, model); + model_graph.add_node(1, 0.5, model, model); + auto& edge = model_graph.add_edge(0, 1, 0.5, {1, 2}); + EXPECT_EQ(edge.traveltime, 0.5); +} + +} // namespace mio \ No newline at end of file