diff --git a/cpp/memilio/geography/regions.h b/cpp/memilio/geography/regions.h index c8dbeb5978..a359d5bce9 100644 --- a/cpp/memilio/geography/regions.h +++ b/cpp/memilio/geography/regions.h @@ -77,6 +77,8 @@ DECL_TYPESAFE(int, CountyId); DECL_TYPESAFE(int, DistrictId); +DECL_TYPESAFE(int, ProvinciaId); + /** * get the id of the state that the specified county is in. * @param[in, out] county a county id. diff --git a/cpp/memilio/io/epi_data.cpp b/cpp/memilio/io/epi_data.cpp index 24e86a00e0..e95c6ef809 100644 --- a/cpp/memilio/io/epi_data.cpp +++ b/cpp/memilio/io/epi_data.cpp @@ -27,10 +27,10 @@ namespace mio std::vector ConfirmedCasesDataEntry::age_group_names = {"A00-A04", "A05-A14", "A15-A34", "A35-A59", "A60-A79", "A80+"}; - -std::vector PopulationDataEntry::age_group_names = { +std::vector PopulationDataEntry::age_group_names = { "<3 years", "3-5 years", "6-14 years", "15-17 years", "18-24 years", "25-29 years", "30-39 years", "40-49 years", "50-64 years", "65-74 years", ">74 years"}; +std::vector PopulationDataEntrySpain::age_group_names = {"Population"}; std::vector VaccinationDataEntry::age_group_names = {"0-4", "5-14", "15-34", "35-59", "60-79", "80-99"}; @@ -49,11 +49,14 @@ IOResult> get_node_ids(const std::string& path, bool is_node_fo } } else { - if (entry.district_id) { + if (entry.state_id) { + id.push_back(entry.state_id->get()); + } + else if (entry.district_id) { id.push_back(entry.district_id->get()); } else { - return failure(StatusCode::InvalidValue, "Population data file is missing district ids."); + return failure(StatusCode::InvalidValue, "Population data file is missing district and state ids."); } } } @@ -62,6 +65,13 @@ IOResult> get_node_ids(const std::string& path, bool is_node_fo id.erase(std::unique(id.begin(), id.end()), id.end()); return success(id); } + +IOResult> get_country_id(const std::string& /*path*/, bool /*is_node_for_county*/, + bool /*rki_age_groups*/) +{ + std::vector id = {0}; + return success(id); +} } // namespace mio #endif //MEMILIO_HAS_JSONCPP diff --git a/cpp/memilio/io/epi_data.h b/cpp/memilio/io/epi_data.h index 3e4572b870..048eb33c12 100644 --- a/cpp/memilio/io/epi_data.h +++ b/cpp/memilio/io/epi_data.h @@ -77,6 +77,9 @@ class ConfirmedCasesNoAgeEntry double num_recovered; double num_deaths; Date date; + boost::optional state_id; + boost::optional county_id; + boost::optional district_id; template static IOResult deserialize(IOContext& io) @@ -86,12 +89,15 @@ class ConfirmedCasesNoAgeEntry auto num_recovered = obj.expect_element("Recovered", Tag{}); auto num_deaths = obj.expect_element("Deaths", Tag{}); auto date = obj.expect_element("Date", Tag{}); + auto state_id = obj.expect_optional("ID_State", Tag{}); + auto county_id = obj.expect_optional("ID_County", Tag{}); + auto district_id = obj.expect_optional("ID_District", Tag{}); return apply( io, - [](auto&& nc, auto&& nr, auto&& nd, auto&& d) { - return ConfirmedCasesNoAgeEntry{nc, nr, nd, d}; + [](auto&& nc, auto&& nr, auto&& nd, auto&& d, auto&& sid, auto&& cid, auto&& did) { + return ConfirmedCasesNoAgeEntry{nc, nr, nd, d, sid, cid, did}; }, - num_confirmed, num_recovered, num_deaths, date); + num_confirmed, num_recovered, num_deaths, date, state_id, county_id, district_id); } }; @@ -142,29 +148,42 @@ class ConfirmedCasesDataEntry { auto obj = io.expect_object("ConfirmedCasesDataEntry"); auto num_confirmed = obj.expect_element("Confirmed", Tag{}); - auto num_recovered = obj.expect_element("Recovered", Tag{}); - auto num_deaths = obj.expect_element("Deaths", Tag{}); + auto num_recovered = obj.expect_optional("Recovered", Tag{}); + auto num_deaths = obj.expect_optional("Deaths", Tag{}); auto date = obj.expect_element("Date", Tag{}); - auto age_group_str = obj.expect_element("Age_RKI", Tag{}); + auto age_group_str = obj.expect_optional("Age_RKI", Tag{}); auto state_id = obj.expect_optional("ID_State", Tag{}); auto county_id = obj.expect_optional("ID_County", Tag{}); auto district_id = obj.expect_optional("ID_District", Tag{}); return apply( io, - [](auto&& nc, auto&& nr, auto&& nd, auto&& d, auto&& a_str, auto&& sid, auto&& cid, + [](auto&& nc, auto&& nr, auto&& nd, auto&& d, auto&& a_str_opt, auto&& sid, auto&& cid, auto&& did) -> IOResult { - auto a = AgeGroup(0); - auto it = std::find(age_group_names.begin(), age_group_names.end(), a_str); - if (it != age_group_names.end()) { - a = AgeGroup(size_t(it - age_group_names.begin())); - } - else if (a_str == "unknown") { - a = AgeGroup(age_group_names.size()); + auto a = AgeGroup(0); // default if Age_RKI missing + + // if no age group is given, use "total" as name + if (!a_str_opt) { + if (age_group_names.size() != 1) { + age_group_names.clear(); + age_group_names.push_back("total"); + } } else { - return failure(StatusCode::InvalidValue, "Invalid confirmed cases data age group."); + const auto& a_str = *a_str_opt; + auto it = std::find(age_group_names.begin(), age_group_names.end(), a_str); + if (it != age_group_names.end()) { + a = AgeGroup(size_t(it - age_group_names.begin())); + } + else if (a_str == "unknown") { + a = AgeGroup(age_group_names.size()); + } + else { + return failure(StatusCode::InvalidValue, "Invalid confirmed cases data age group."); + } } - return success(ConfirmedCasesDataEntry{nc, nr, nd, d, a, sid, cid, did}); + double nrec = nr ? *nr : 0.0; + double ndead = nd ? *nd : 0.0; + return success(ConfirmedCasesDataEntry{nc, nrec, ndead, d, a, sid, cid, did}); }, num_confirmed, num_recovered, num_deaths, date, age_group_str, state_id, county_id, district_id); } @@ -331,6 +350,35 @@ class PopulationDataEntry } }; +class PopulationDataEntrySpain +{ +public: + static std::vector age_group_names; + + CustomIndexArray population; + boost::optional provincia_id; + + template + static IOResult deserialize(IoContext& io) + { + auto obj = io.expect_object("PopulationDataEntrySpain"); + auto provincia = obj.expect_optional("ID_Provincia", Tag{}); + std::vector> age_groups; + age_groups.reserve(age_group_names.size()); + std::transform(age_group_names.begin(), age_group_names.end(), std::back_inserter(age_groups), + [&obj](auto&& age_name) { + return obj.expect_element(age_name, Tag{}); + }); + return apply( + io, + [](auto&& ag, auto&& pid) { + return PopulationDataEntrySpain{ + CustomIndexArray(AgeGroup(ag.size()), ag.begin(), ag.end()), pid}; + }, + details::unpack_all(age_groups), provincia); + } +}; + namespace details { inline void get_rki_age_interpolation_coefficients(const std::vector& age_ranges, @@ -435,6 +483,19 @@ inline IOResult> deserialize_population_data(co } } +/** + * Deserialize population data from a JSON value. + * Age groups are interpolated to RKI age groups. + * @param jsvalue JSON value that contains the population data. + * @param rki_age_groups Specifies whether population data should be interpolated to rki age groups. + * @return list of population data. + */ +inline IOResult> deserialize_population_data_spain(const Json::Value& jsvalue) +{ + BOOST_OUTCOME_TRY(auto&& population_data, deserialize_json(jsvalue, Tag>{})); + return success(population_data); +} + /** * Deserialize population data from a JSON file. * Age groups are interpolated to RKI age groups. @@ -448,6 +509,18 @@ inline IOResult> read_population_data(const std return deserialize_population_data(jsvalue, rki_age_group); } +/** + * Deserialize population data from a JSON file. + * Age groups are interpolated to RKI age groups. + * @param filename JSON file that contains the population data. + * @return list of population data. + */ +inline IOResult> read_population_data_spain(const std::string& filename) +{ + BOOST_OUTCOME_TRY(auto&& jsvalue, read_json(filename)); + return deserialize_population_data_spain(jsvalue); +} + /** * @brief Sets the age groups' names for the ConfirmedCasesDataEntry%s. * @param[in] names age group names @@ -474,6 +547,8 @@ IOResult set_vaccination_data_age_group_names(std::vector nam * @return list of node ids. */ IOResult> get_node_ids(const std::string& path, bool is_node_for_county, bool rki_age_groups = true); +IOResult> get_country_id(const std::string& /*path*/, bool /*is_node_for_county*/, + bool /*rki_age_groups*/ = true); /** * Represents an entry in a vaccination data file. diff --git a/cpp/memilio/io/parameters_io.cpp b/cpp/memilio/io/parameters_io.cpp index f9f9f8b3c0..450a38f902 100644 --- a/cpp/memilio/io/parameters_io.cpp +++ b/cpp/memilio/io/parameters_io.cpp @@ -62,11 +62,43 @@ IOResult>> read_population_data(const std::vecto return success(vnum_population); } +IOResult>> +read_population_data_spain(const std::vector& population_data, + const std::vector& vregion) +{ + std::vector> vnum_population(vregion.size(), std::vector(1, 0.0)); + + for (auto&& provincia_entry : population_data) { + if (provincia_entry.provincia_id) { + for (size_t idx = 0; idx < vregion.size(); ++idx) { + if (vregion[idx] == provincia_entry.provincia_id->get()) { + vnum_population[idx][0] += provincia_entry.population[AgeGroup(0)]; + } + } + } + // id 0 means the whole country + for (size_t idx = 0; idx < vregion.size(); ++idx) { + if (vregion[idx] == 0) { + vnum_population[idx][0] += provincia_entry.population[AgeGroup(0)]; + } + } + } + + return success(vnum_population); +} + IOResult>> read_population_data(const std::string& path, const std::vector& vregion) { BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data(path)); return read_population_data(population_data, vregion); } + +IOResult>> read_population_data_spain(const std::string& path, + const std::vector& vregion) +{ + BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data_spain(path)); + return read_population_data_spain(population_data, vregion); +} } // namespace mio #endif //MEMILIO_HAS_JSONCPP diff --git a/cpp/memilio/io/parameters_io.h b/cpp/memilio/io/parameters_io.h index 3ea00fc73b..bbb0769d51 100644 --- a/cpp/memilio/io/parameters_io.h +++ b/cpp/memilio/io/parameters_io.h @@ -136,6 +136,24 @@ IOResult>> read_population_data(const std::vecto IOResult>> read_population_data(const std::string& path, const std::vector& vregion); +/** + * @brief Reads (Spain) population data from a vector of Spain population data entries. + * Uses provincias IDs and a single aggregated age group. + * @param[in] population_data Vector of Spain population data entries. + * @param[in] vregion Vector of keys representing the provincias (or country = 0) of interest. + */ +IOResult>> +read_population_data_spain(const std::vector& population_data, + const std::vector& vregion); + +/** + * @brief Reads (Spain) population data from census data file with provincias. + * @param[in] path Path to the provincias population data file. + * @param[in] vregion Vector of keys representing the provincias (or country = 0) of interest. + */ +IOResult>> read_population_data_spain(const std::string& path, + const std::vector& vregion); + } // namespace mio #endif //MEMILIO_HAS_JSONCPP diff --git a/cpp/memilio/mobility/graph.h b/cpp/memilio/mobility/graph.h index eb8bd35990..1a616be241 100644 --- a/cpp/memilio/mobility/graph.h +++ b/cpp/memilio/mobility/graph.h @@ -307,15 +307,15 @@ IOResult set_nodes(const Parameters& params, Date start_date, Date end_dat }); //uncertainty in populations - for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { - for (auto j = Index(0); j < Model::Compartments::Count; ++j) { - auto& compartment_value = nodes[node_idx].populations[{i, j}]; - compartment_value = - 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))); - } - } + // for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) { + // for (auto j = Index(0); j < Model::Compartments::Count; ++j) { + // auto& compartment_value = nodes[node_idx].populations[{i, j}]; + // compartment_value = + // 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))); + // } + // } params_graph.add_node(node_ids[node_idx], nodes[node_idx]); } @@ -368,7 +368,7 @@ IOResult set_edges(const fs::path& mobility_data_file, Graph(); ++age) { for (auto compartment : mobile_compartments) { auto coeff_index = populations.get_flat_index({age, compartment}); - mobility_coeffs[size_t(ContactLocation::Work)].get_baseline()[coeff_index] = + mobility_coeffs[size_t(ContactLocation::Home)].get_baseline()[coeff_index] = commuter_coeff_ij * commuting_weights[size_t(age)]; } } diff --git a/cpp/models/ode_secir/parameters_io.h b/cpp/models/ode_secir/parameters_io.h index ab334d922d..6f9a89acda 100644 --- a/cpp/models/ode_secir/parameters_io.h +++ b/cpp/models/ode_secir/parameters_io.h @@ -75,7 +75,17 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect const std::vector& scaling_factor_inf) { const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size(); - assert(scaling_factor_inf.size() == num_age_groups); + // allow single scalar scaling that is broadcast to all age groups + assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups); + + // Set scaling factors to match num age groups + std::vector scaling_factor_inf_full; + if (scaling_factor_inf.size() == 1) { + scaling_factor_inf_full.assign(num_age_groups, scaling_factor_inf[0]); + } + else { + scaling_factor_inf_full = scaling_factor_inf; + } std::vector> t_InfectedNoSymptoms{model.size()}; std::vector> t_Exposed{model.size()}; @@ -89,7 +99,12 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect std::vector> mu_U_D{model.size()}; for (size_t node = 0; node < model.size(); ++node) { - for (size_t group = 0; group < num_age_groups; group++) { + const size_t model_groups = (size_t)model[node].parameters.get_num_groups(); + assert(model_groups == 1 || model_groups == num_age_groups); + for (size_t ag = 0; ag < num_age_groups; ag++) { + // If the model has fewer groups than casedata entries available, + // reuse group 0 parameters for all RKI age groups + const size_t group = (model_groups == num_age_groups) ? ag : 0; t_Exposed[node].push_back( static_cast(std::round(model[node].parameters.template get>()[(AgeGroup)group]))); @@ -121,26 +136,49 @@ IOResult set_confirmed_cases_data(std::vector>& model, std::vect BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec, t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere, - t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf)); + t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf_full)); for (size_t node = 0; node < model.size(); node++) { if (std::accumulate(num_InfectedSymptoms[node].begin(), num_InfectedSymptoms[node].end(), 0.0) > 0) { size_t num_groups = (size_t)model[node].parameters.get_num_groups(); - for (size_t i = 0; i < num_groups; i++) { - model[node].populations[{AgeGroup(i), InfectionState::Exposed}] = num_Exposed[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptoms}] = - num_InfectedNoSymptoms[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsConfirmed}] = 0; - model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptoms}] = - num_InfectedSymptoms[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptomsConfirmed}] = 0; - model[node].populations[{AgeGroup(i), InfectionState::InfectedSevere}] = num_InfectedSevere[node][i]; - // Only set the number of ICU patients here, if the date is not available in the data. + if (num_groups == num_age_groups) { + for (size_t i = 0; i < num_groups; i++) { + model[node].populations[{AgeGroup(i), InfectionState::Exposed}] = num_Exposed[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptoms}] = + num_InfectedNoSymptoms[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptoms}] = + num_InfectedSymptoms[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i), InfectionState::InfectedSevere}] = + num_InfectedSevere[node][i]; + // Only set the number of ICU patients here, if the date is not available in the data. + if (!is_divi_data_available(date)) { + model[node].populations[{AgeGroup(i), InfectionState::InfectedCritical}] = num_icu[node][i]; + } + model[node].populations[{AgeGroup(i), InfectionState::Dead}] = num_death[node][i]; + model[node].populations[{AgeGroup(i), InfectionState::Recovered}] = num_rec[node][i]; + } + } + else { + const auto sum_vec = [](const std::vector& v) { + return std::accumulate(v.begin(), v.end(), 0.0); + }; + const size_t i0 = 0; + model[node].populations[{AgeGroup(i0), InfectionState::Exposed}] = sum_vec(num_Exposed[node]); + model[node].populations[{AgeGroup(i0), InfectionState::InfectedNoSymptoms}] = + sum_vec(num_InfectedNoSymptoms[node]); + model[node].populations[{AgeGroup(i0), InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i0), InfectionState::InfectedSymptoms}] = + sum_vec(num_InfectedSymptoms[node]); + model[node].populations[{AgeGroup(i0), InfectionState::InfectedSymptomsConfirmed}] = 0; + model[node].populations[{AgeGroup(i0), InfectionState::InfectedSevere}] = + sum_vec(num_InfectedSevere[node]); if (!is_divi_data_available(date)) { - model[node].populations[{AgeGroup(i), InfectionState::InfectedCritical}] = num_icu[node][i]; + model[node].populations[{AgeGroup(i0), InfectionState::InfectedCritical}] = sum_vec(num_icu[node]); } - model[node].populations[{AgeGroup(i), InfectionState::Dead}] = num_death[node][i]; - model[node].populations[{AgeGroup(i), InfectionState::Recovered}] = num_rec[node][i]; + model[node].populations[{AgeGroup(i0), InfectionState::Dead}] = sum_vec(num_death[node]); + model[node].populations[{AgeGroup(i0), InfectionState::Recovered}] = sum_vec(num_rec[node]); } } else { @@ -231,10 +269,20 @@ IOResult set_population_data(std::vector>& model, assert(num_population.size() == vregion.size()); assert(model.size() == vregion.size()); for (size_t region = 0; region < vregion.size(); region++) { - auto num_groups = model[region].parameters.get_num_groups(); - for (auto i = AgeGroup(0); i < num_groups; i++) { + const auto model_groups = (size_t)model[region].parameters.get_num_groups(); + const auto data_groups = num_population[region].size(); + assert(data_groups == model_groups || (model_groups == 1 && data_groups >= 1)); + + if (data_groups == model_groups) { + for (auto i = AgeGroup(0); i < model[region].parameters.get_num_groups(); i++) { + model[region].populations.template set_difference_from_group_total( + {i, InfectionState::Susceptible}, num_population[region][(size_t)i]); + } + } + else if (model_groups == 1 && data_groups >= 1) { + const double total = std::accumulate(num_population[region].begin(), num_population[region].end(), 0.0); model[region].populations.template set_difference_from_group_total( - {i, InfectionState::Susceptible}, num_population[region][size_t(i)]); + {AgeGroup(0), InfectionState::Susceptible}, total); } } return success(); @@ -256,6 +304,15 @@ IOResult set_population_data(std::vector>& model, const std::str return success(); } +template +IOResult set_population_data_provincias(std::vector>& model, const std::string& path, + const std::vector& vregion) +{ + BOOST_OUTCOME_TRY(const auto&& num_population, read_population_data_spain(path, vregion)); + BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion)); + return success(); +} + } //namespace details #ifdef MEMILIO_HAS_HDF5 @@ -283,8 +340,8 @@ IOResult export_input_data_county_timeseries( const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path) { const auto num_age_groups = (size_t)models[0].parameters.get_num_groups(); - assert(scaling_factor_inf.size() == num_age_groups); - assert(num_age_groups == ConfirmedCasesDataEntry::age_group_names.size()); + // allow scalar scaling factor as convenience for 1-group models + assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups); assert(models.size() == region.size()); std::vector> extrapolated_data( region.size(), TimeSeries::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups)); @@ -333,9 +390,10 @@ IOResult export_input_data_county_timeseries( * @param[in] pydata_dir Directory of files. */ template -IOResult read_input_data_germany(std::vector& model, Date date, +IOResult read_input_data_germany(std::vector& model, Date date, std::vector& /*state*/, const std::vector& scaling_factor_inf, double scaling_factor_icu, - const std::string& pydata_dir) + const std::string& pydata_dir, int /*num_days*/ = 0, + bool /*export_time_series*/ = false) { BOOST_OUTCOME_TRY( details::set_divi_data(model, path_join(pydata_dir, "germany_divi.json"), {0}, date, scaling_factor_icu)); @@ -358,7 +416,8 @@ IOResult read_input_data_germany(std::vector& model, Date date, template IOResult read_input_data_state(std::vector& model, Date date, std::vector& state, const std::vector& scaling_factor_inf, double scaling_factor_icu, - const std::string& pydata_dir) + const std::string& pydata_dir, int /*num_days*/ = 0, + bool /*export_time_series*/ = false) { BOOST_OUTCOME_TRY( @@ -407,6 +466,33 @@ IOResult read_input_data_county(std::vector& model, Date date, cons return success(); } +/** + * @brief Reads population data from population files for the specefied county. + * @param[in, out] model Vector of model in which the data is set. + * @param[in] date Date for which the data should be read. + * @param[in] county Vector of region keys of counties of interest. + * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data. + * @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data. + * @param[in] pydata_dir Directory of files. + * @param[in] num_days [Default: 0] Number of days to be simulated; required to extrapolate real data. + * @param[in] export_time_series [Default: false] If true, reads data for each day of simulation and writes it in the same directory as the input files. + */ +template +IOResult read_input_data_provincias(std::vector& model, Date date, const std::vector& provincias, + const std::vector& scaling_factor_inf, double scaling_factor_icu, + const std::string& pydata_dir, int /*num_days*/ = 0, + bool /*export_time_series*/ = false) +{ + BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "provincia_icu.json"), provincias, date, + scaling_factor_icu)); + + BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_pronvincias.json"), + provincias, date, scaling_factor_inf)); + BOOST_OUTCOME_TRY(details::set_population_data_provincias( + model, path_join(pydata_dir, "provincias_current_population.json"), provincias)); + return success(); +} + /** * @brief reads population data from population files for the specified nodes * @param[in, out] model vector of model in which the data is set diff --git a/cpp/tests/test_odesecir.cpp b/cpp/tests/test_odesecir.cpp index cb4ad98a2f..b59da559b1 100755 --- a/cpp/tests/test_odesecir.cpp +++ b/cpp/tests/test_odesecir.cpp @@ -29,6 +29,7 @@ #include "memilio/io/parameters_io.h" #include "memilio/data/analyze_result.h" #include "memilio/math/adapt_rk.h" +#include "memilio/geography/regions.h" #include @@ -1516,5 +1517,184 @@ TEST(TestOdeSecir, read_population_data_failure) EXPECT_EQ(result.error().message(), "File with county population expected."); } +TEST(TestOdeSecirIO, read_input_data_county_aggregates_one_group) +{ + // Set up two models with different number of age groups. + const size_t num_age_groups = 6; + std::vector> models6{mio::osecir::Model((int)num_age_groups)}; + std::vector> models1{mio::osecir::Model(1)}; + + // Relevant parameters for model with 6 age groups + for (auto i = mio::AgeGroup(0); i < (mio::AgeGroup)num_age_groups; ++i) { + models6[0].parameters.get>()[i] = 0.2; + models6[0].parameters.get>()[i] = 0.25; + } + + // Relevant parameters for model with 1 age group + models1[0].parameters.get>()[mio::AgeGroup(0)] = 0.2; + models1[0].parameters.get>()[mio::AgeGroup(0)] = 0.25; + + const auto pydata_dir_Germany = mio::path_join(TEST_DATA_DIR, "Germany", "pydata"); + const std::vector counties{1002}; + const auto date = mio::Date(2020, 12, 1); + + std::vector scale6(num_age_groups, 1.0); + std::vector scale1{1.0}; + + // Initialize both models + ASSERT_THAT(mio::osecir::read_input_data_county(models6, date, counties, scale6, 1.0, pydata_dir_Germany), + IsSuccess()); + ASSERT_THAT(mio::osecir::read_input_data_county(models1, date, counties, scale1, 1.0, pydata_dir_Germany), + IsSuccess()); + + // Aggreagate the results from the model with 6 age groups and compare with the model with 1 age group + const auto& m6 = models6[0]; + const auto& m1 = models1[0]; + const double tol = 1e-10; + for (int s = 0; s < (int)mio::osecir::InfectionState::Count; ++s) { + double sum6 = 0.0; + for (size_t ag = 0; ag < num_age_groups; ++ag) { + sum6 += m6.populations[{mio::AgeGroup(ag), (mio::osecir::InfectionState)s}].value(); + } + const double v1 = m1.populations[{mio::AgeGroup(0), (mio::osecir::InfectionState)s}].value(); + EXPECT_NEAR(sum6, v1, tol); + } + + // Total population + EXPECT_NEAR(m6.populations.get_total(), m1.populations.get_total(), tol); +} + +TEST(TestOdeSecirIO, set_population_data_single_age_group) +{ + const size_t num_age_groups = 6; + + // Create two models: one with 6 age groups, one with 1 age group + std::vector> models6{mio::osecir::Model((int)num_age_groups)}; + std::vector> models1{mio::osecir::Model(1)}; + + // Test population data with 6 different values for age groups + std::vector> population_data6 = {{10000.0, 20000.0, 30000.0, 25000.0, 15000.0, 8000.0}}; + std::vector> population_data1 = {{108000.0}}; // sum of all age groups + std::vector regions = {1002}; + + // Set population data for both models + EXPECT_THAT(mio::osecir::details::set_population_data(models6, population_data6, regions), IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_population_data(models1, population_data1, regions), IsSuccess()); + + // Sum all compartments across age groups in 6-group model and compare 1-group model + const double tol = 1e-10; + for (int s = 0; s < (int)mio::osecir::InfectionState::Count; ++s) { + double sum6 = 0.0; + for (size_t ag = 0; ag < num_age_groups; ++ag) { + sum6 += models6[0].populations[{mio::AgeGroup(ag), (mio::osecir::InfectionState)s}].value(); + } + double val1 = models1[0].populations[{mio::AgeGroup(0), (mio::osecir::InfectionState)s}].value(); + + EXPECT_NEAR(sum6, val1, tol); + } + + // Total population should also match + EXPECT_NEAR(models6[0].populations.get_total(), models1[0].populations.get_total(), tol); +} + +TEST(TestOdeSecirIO, set_confirmed_cases_data_single_age_group) +{ + const size_t num_age_groups = 6; + + // Create two models: one with 6 age groups, one with 1 age group + std::vector> models6{mio::osecir::Model((int)num_age_groups)}; + std::vector> models1{mio::osecir::Model(1)}; + + // Create case data for all 6 age groups over multiple days (current day + 6 days back) + std::vector case_data; + + for (int day_offset = -6; day_offset <= 0; ++day_offset) { + mio::Date current_date = mio::offset_date_by_days(mio::Date(2020, 12, 1), day_offset); + + for (int age_group = 0; age_group < 6; ++age_group) { + double base_confirmed = 80.0 + age_group * 8.0 + (day_offset + 6) * 5.0; + double base_recovered = 40.0 + age_group * 4.0 + (day_offset + 6) * 3.0; + double base_deaths = 3.0 + age_group * 0.5 + (day_offset + 6) * 0.5; + + mio::ConfirmedCasesDataEntry entry{base_confirmed, + base_recovered, + base_deaths, + current_date, + mio::AgeGroup(age_group), + {}, + mio::regions::CountyId(1002), + {}}; + case_data.push_back(entry); + } + } + + std::vector regions = {1002}; + std::vector scaling_factors = {1.0}; + + // Set confirmed cases data for both models + EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models6, case_data, regions, mio::Date(2020, 12, 1), + scaling_factors), + IsSuccess()); + EXPECT_THAT(mio::osecir::details::set_confirmed_cases_data(models1, case_data, regions, mio::Date(2020, 12, 1), + scaling_factors), + IsSuccess()); + + // Sum all compartments across age groups in 6-group model should be equal to 1-group model + for (int s = 0; s < (int)mio::osecir::InfectionState::Count; ++s) { + double sum6 = 0.0; + for (size_t ag = 0; ag < num_age_groups; ++ag) { + sum6 += models6[0].populations[{mio::AgeGroup(ag), (mio::osecir::InfectionState)s}].value(); + } + + double val1 = models1[0].populations[{mio::AgeGroup(0), (mio::osecir::InfectionState)s}].value(); + + EXPECT_NEAR(sum6, val1, 1e-10); + } + + // Total population + EXPECT_NEAR(models6[0].populations.get_total(), models1[0].populations.get_total(), 1e-10); +} + +TEST(TestOdeSecirIO, set_divi_data_single_age_group) +{ + // Create models with 6 age groups and 1 age group + std::vector> models_6_groups{mio::osecir::Model(6)}; + std::vector> models_1_group{mio::osecir::Model(1)}; + + // Set relevant parameters for all age groups + for (int i = 0; i < 6; i++) { + models_6_groups[0].parameters.get>()[mio::AgeGroup(i)] = 0.2; + models_6_groups[0].parameters.get>()[mio::AgeGroup(i)] = 0.25; + } + + // Set relevant parameters for 1 age group model + models_1_group[0].parameters.get>()[mio::AgeGroup(0)] = 0.2; + models_1_group[0].parameters.get>()[mio::AgeGroup(0)] = 0.25; + + // Apply DIVI data to both models + std::vector regions = {1002}; + double scaling_factor_icu = 1.0; + mio::Date date(2020, 12, 1); + std::string divi_data_path = mio::path_join(TEST_DATA_DIR, "Germany", "pydata", "county_divi_ma7.json"); + auto result_6_groups = + mio::osecir::details::set_divi_data(models_6_groups, divi_data_path, regions, date, scaling_factor_icu); + auto result_1_group = + mio::osecir::details::set_divi_data(models_1_group, divi_data_path, regions, date, scaling_factor_icu); + + EXPECT_THAT(result_6_groups, IsSuccess()); + EXPECT_THAT(result_1_group, IsSuccess()); + + // Calculate totals after applying DIVI data + double total_icu_6_groups_after = 0.0; + for (int i = 0; i < 6; i++) { + total_icu_6_groups_after += + models_6_groups[0].populations[{mio::AgeGroup(i), mio::osecir::InfectionState::InfectedCritical}].value(); + } + double icu_1_group_after = + models_1_group[0].populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}].value(); + + EXPECT_NEAR(total_icu_6_groups_after, icu_1_group_after, 1e-10); +} + #endif #endif diff --git a/pycode/examples/simulation/graph_germany_nuts0.py b/pycode/examples/simulation/graph_germany_nuts0.py new file mode 100644 index 0000000000..b9d4843430 --- /dev/null +++ b/pycode/examples/simulation/graph_germany_nuts0.py @@ -0,0 +1,250 @@ +############################################################################# +# Copyright (C) 2020-2025 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. +############################################################################# +import numpy as np +import datetime +import os +import memilio.simulation as mio +import memilio.simulation.osecir as osecir +import matplotlib.pyplot as plt + +from enum import Enum +from memilio.simulation.osecir import (Model, Simulation, + interpolate_simulation_result) + +import pickle + + +class Simulation: + """ """ + + def __init__(self, data_dir, start_date, results_dir): + self.num_groups = 1 + self.data_dir = data_dir + self.start_date = start_date + self.results_dir = results_dir + if not os.path.exists(self.results_dir): + os.makedirs(self.results_dir) + + def set_covid_parameters(self, model): + """ + + :param model: + + """ + model.parameters.TimeExposed[mio.AgeGroup(0)] = 3.335 + model.parameters.TimeInfectedNoSymptoms[mio.AgeGroup(0)] = 2.58916 + model.parameters.TimeInfectedSymptoms[mio.AgeGroup(0)] = 6.94547 + model.parameters.TimeInfectedSevere[mio.AgeGroup(0)] = 7.28196 + model.parameters.TimeInfectedCritical[mio.AgeGroup(0)] = 13.066 + + # probabilities + model.parameters.TransmissionProbabilityOnContact[mio.AgeGroup( + 0)] = 0.07333 + model.parameters.RelativeTransmissionNoSymptoms[mio.AgeGroup(0)] = 1 + + model.parameters.RecoveredPerInfectedNoSymptoms[mio.AgeGroup( + 0)] = 0.2069 + model.parameters.SeverePerInfectedSymptoms[mio.AgeGroup(0)] = 0.07864 + model.parameters.CriticalPerSevere[mio.AgeGroup(0)] = 0.17318 + model.parameters.DeathsPerCritical[mio.AgeGroup(0)] = 0.21718 + + # start day is set to the n-th day of the year + model.parameters.StartDay = self.start_date.timetuple().tm_yday + + model.parameters.Seasonality = mio.UncertainValue(0.2) + + def set_contact_matrices(self, model): + """ + + :param model: + + """ + contact_matrices = mio.ContactMatrixGroup(1, self.num_groups) + + baseline = np.ones((self.num_groups, self.num_groups)) * 7.95 + minimum = np.ones((self.num_groups, self.num_groups)) * 0 + contact_matrices[0] = mio.ContactMatrix(baseline, minimum) + model.parameters.ContactPatterns.cont_freq_mat = contact_matrices + + def set_npis(self, params, end_date, damping_value): + """ + + :param params: + :param end_date: + + """ + + start_damping = datetime.date( + 2020, 12, 18) + + if start_damping < end_date: + start_date = (start_damping - self.start_date).days + params.ContactPatterns.cont_freq_mat[0].add_damping( + mio.Damping(np.r_[damping_value], t=start_date)) + + def get_graph(self, end_date, damping_value): + """ + + :param end_date: + + """ + print("Initializing model...") + model = Model(self.num_groups) + self.set_covid_parameters(model) + self.set_contact_matrices(model) + self.set_npis(model.parameters, end_date, damping_value) + print("Model initialized.") + + graph = osecir.ModelGraph() + + scaling_factor_infected = [2.5] + scaling_factor_icu = 1.0 + + data_dir_Germany = os.path.join(self.data_dir, "Germany") + pydata_dir = os.path.join(data_dir_Germany, "pydata") + + path_population_data = os.path.join( + pydata_dir, "county_current_population_germany.json") + + print("Setting nodes...") + mio.osecir.set_node_germany( + model.parameters, + mio.Date(self.start_date.year, + self.start_date.month, self.start_date.day), + mio.Date(end_date.year, + end_date.month, end_date.day), pydata_dir, + path_population_data, False, graph, scaling_factor_infected, + scaling_factor_icu, 0.9, 0, False) + + print("Graph created.") + + return graph + + def run(self, num_days_sim, damping_value, save_graph=True): + """ + + :param num_days_sim: + :param num_runs: (Default value = 10) + :param save_graph: (Default value = True) + :param create_gif: (Default value = True) + + """ + mio.set_log_level(mio.LogLevel.Warning) + end_date = self.start_date + datetime.timedelta(days=num_days_sim) + num_runs = 10 + + graph = self.get_graph(end_date, damping_value) + + if save_graph: + path_graph = os.path.join(self.results_dir, "graph") + if not os.path.exists(path_graph): + os.makedirs(path_graph) + osecir.write_graph(graph, path_graph) + + mobility_graph = osecir.MobilityGraph() + for node_idx in range(graph.num_nodes): + mobility_graph.add_node(graph.get_node( + node_idx).id, graph.get_node(node_idx).property) + for edge_idx in range(graph.num_edges): + mobility_graph.add_edge( + graph.get_edge(edge_idx).start_node_idx, + graph.get_edge(edge_idx).end_node_idx, + graph.get_edge(edge_idx).property) + + mobility_sim = osecir.MobilitySimulation(mobility_graph, t0=0, dt=0.5) + mobility_sim.advance(num_days_sim) + + results = [] + for node_idx in range(graph.num_nodes): + results.append(osecir.interpolate_simulation_result( + mobility_sim.graph.get_node(node_idx).property.result)) + + osecir.interpolate_simulation_result( + mobility_sim.graph.get_node(0).property.result).export_csv( + 'test.csv') + + return results + + +def run_germany_nuts0_simulation(damping_value): + mio.set_log_level(mio.LogLevel.Warning) + file_path = os.path.dirname(os.path.abspath(__file__)) + + sim = Simulation( + data_dir=os.path.join(file_path, "../../../data"), + start_date=datetime.date(year=2020, month=12, day=12), + results_dir=os.path.join(file_path, "../../../results_osecir")) + num_days_sim = 50 + + results = sim.run(num_days_sim, damping_value) + + return {"region" + str(region): results[region] for region in range(len(results))} + + +def prior(): + damping_value = np.random.uniform(0.0, 1.0) + return {"damping_value": damping_value} + + +if __name__ == "__main__": + + run_germany_nuts0_simulation(0.5) + # import os + # os.environ["KERAS_BACKEND"] = "jax" + + # import bayesflow as bf + + # simulator = bf.simulators.make_simulator([prior, run_germany_nuts3_simulation]) + # # trainings_data = simulator.sample(5) + + # with open('trainings_data.pickle', 'wb') as f: + # pickle.dump(trainings_data, f, pickle.HIGHEST_PROTOCOL) + + # with open('trainings_data.pickle', 'rb') as f: + # trainings_data = pickle.load(f) + + # # trainings_data = {k:v for k, v in trainings_data.items() if k in ('damping_value', 'region0', 'region1')} + # print("Loaded training data:", trainings_data) + + # trainings_data = simulator.sample(2) + # validation_data = simulator.sample(2) + + # adapter = ( + # bf.Adapter() + # .to_array() + # .convert_dtype("float64", "float32") + # .constrain("damping_value", lower=0.0, upper=1.0) + # .concatenate(["region"+str(region) for region in range(len(trainings_data)-1)], into="summary_variables") + # .rename("damping_value", "inference_variables") + # #.standardize("summary_variables") + # ) + + # summary_network = bf.networks.TimeSeriesNetwork(summary_dim=4) + # inference_network = bf.networks.CouplingFlow() + + # workflow = bf.BasicWorkflow( + # simulator=simulator, + # adapter=adapter, + # summary_network=summary_network, + # inference_network=inference_network + # ) + + # history = workflow.fit_offline(data=trainings_data, epochs=2, batch_size=2, validation_data=trainings_data) + # f = bf.diagnostics.plots.loss(history) diff --git a/pycode/examples/simulation/graph_germany_nuts1.py b/pycode/examples/simulation/graph_germany_nuts1.py new file mode 100644 index 0000000000..40389ee00f --- /dev/null +++ b/pycode/examples/simulation/graph_germany_nuts1.py @@ -0,0 +1,274 @@ +############################################################################# +# Copyright (C) 2020-2025 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. +############################################################################# +import numpy as np +import datetime +import os +import memilio.simulation as mio +import memilio.simulation.osecir as osecir +import matplotlib.pyplot as plt + +from enum import Enum +from memilio.simulation.osecir import (Model, Simulation, + interpolate_simulation_result) + +import pickle + + +class Simulation: + """ """ + + def __init__(self, data_dir, start_date, results_dir): + self.num_groups = 1 + self.data_dir = data_dir + self.start_date = start_date + self.results_dir = results_dir + if not os.path.exists(self.results_dir): + os.makedirs(self.results_dir) + + def set_covid_parameters(self, model): + """ + + :param model: + + """ + model.parameters.TimeExposed[mio.AgeGroup(0)] = 3.335 + model.parameters.TimeInfectedNoSymptoms[mio.AgeGroup(0)] = 2.58916 + model.parameters.TimeInfectedSymptoms[mio.AgeGroup(0)] = 6.94547 + model.parameters.TimeInfectedSevere[mio.AgeGroup(0)] = 7.28196 + model.parameters.TimeInfectedCritical[mio.AgeGroup(0)] = 13.066 + + # probabilities + model.parameters.TransmissionProbabilityOnContact[mio.AgeGroup( + 0)] = 0.07333 + model.parameters.RelativeTransmissionNoSymptoms[mio.AgeGroup(0)] = 1 + + model.parameters.RecoveredPerInfectedNoSymptoms[mio.AgeGroup( + 0)] = 0.2069 + model.parameters.SeverePerInfectedSymptoms[mio.AgeGroup(0)] = 0.07864 + model.parameters.CriticalPerSevere[mio.AgeGroup(0)] = 0.17318 + model.parameters.DeathsPerCritical[mio.AgeGroup(0)] = 0.21718 + + # start day is set to the n-th day of the year + model.parameters.StartDay = self.start_date.timetuple().tm_yday + + model.parameters.Seasonality = mio.UncertainValue(0.2) + + def set_contact_matrices(self, model): + """ + + :param model: + + """ + contact_matrices = mio.ContactMatrixGroup(1, self.num_groups) + + baseline = np.ones((self.num_groups, self.num_groups)) * 7.95 + minimum = np.ones((self.num_groups, self.num_groups)) * 0 + contact_matrices[0] = mio.ContactMatrix(baseline, minimum) + model.parameters.ContactPatterns.cont_freq_mat = contact_matrices + + def set_npis(self, params, end_date, damping_value): + """ + + :param params: + :param end_date: + + """ + + start_damping = datetime.date( + 2020, 12, 18) + + if start_damping < end_date: + start_date = (start_damping - self.start_date).days + params.ContactPatterns.cont_freq_mat[0].add_damping( + mio.Damping(np.r_[damping_value], t=start_date)) + + def get_graph(self, end_date): + """ + + :param end_date: + + """ + print("Initializing model...") + model = Model(self.num_groups) + self.set_covid_parameters(model) + self.set_contact_matrices(model) + print("Model initialized.") + + graph = osecir.ModelGraph() + + scaling_factor_infected = [2.5] + scaling_factor_icu = 1.0 + tnt_capacity_factor = 7.5 / 100000. + + data_dir_Germany = os.path.join(self.data_dir, "Germany") + mobility_data_file = os.path.join( + data_dir_Germany, "mobility", "commuter_mobility_2022_states.txt") + pydata_dir = os.path.join(data_dir_Germany, "pydata") + + path_population_data = os.path.join( + pydata_dir, "county_current_population_states.json") + + print("Setting nodes...") + mio.osecir.set_nodes_states( + model.parameters, + mio.Date(self.start_date.year, + self.start_date.month, self.start_date.day), + mio.Date(end_date.year, + end_date.month, end_date.day), pydata_dir, + path_population_data, False, graph, scaling_factor_infected, + scaling_factor_icu, 0, 0, False) + + print("Setting edges...") + mio.osecir.set_edges( + mobility_data_file, graph, 1) + + print("Graph created.") + + return graph + + def run(self, num_days_sim, damping_values, save_graph=True): + """ + + :param num_days_sim: + :param num_runs: (Default value = 10) + :param save_graph: (Default value = True) + :param create_gif: (Default value = True) + + """ + mio.set_log_level(mio.LogLevel.Warning) + end_date = self.start_date + datetime.timedelta(days=num_days_sim) + + graph = self.get_graph(end_date) + + if save_graph: + path_graph = os.path.join(self.results_dir, "graph") + if not os.path.exists(path_graph): + os.makedirs(path_graph) + osecir.write_graph(graph, path_graph) + + mobility_graph = osecir.MobilityGraph() + for node_idx in range(graph.num_nodes): + # if node_idx < 5: + node = graph.get_node(node_idx) + self.set_npis(node.property.parameters, + end_date, damping_values[node_idx]) + mobility_graph.add_node(node.id, node.property) + # else: + # node = graph.get_node(node_idx) + # self.set_npis(node.property.parameters, end_date, 0.5) + # mobility_graph.add_node(node.id, node.property) + for edge_idx in range(graph.num_edges): + mobility_graph.add_edge( + graph.get_edge(edge_idx).start_node_idx, + graph.get_edge(edge_idx).end_node_idx, + graph.get_edge(edge_idx).property) + mobility_sim = osecir.MobilitySimulation(mobility_graph, t0=0, dt=0.5) + mobility_sim.advance(num_days_sim) + + results = [] + for node_idx in range(graph.num_nodes): + results.append(osecir.interpolate_simulation_result( + mobility_sim.graph.get_node(node_idx).property.result)) + + return results + + +def run_germany_nuts1_simulation(damping_values): + mio.set_log_level(mio.LogLevel.Warning) + file_path = os.path.dirname(os.path.abspath(__file__)) + + sim = Simulation( + data_dir=os.path.join(file_path, "../../../data"), + start_date=datetime.date(year=2020, month=12, day=12), + results_dir=os.path.join(file_path, "../../../results_osecir")) + num_days_sim = 50 + + results = sim.run(num_days_sim, damping_values) + results[0].export_csv('test.csv') + + return {f'region{region}': results[region] for region in range(len(results))} + + +def prior(): + damping_values = np.random.uniform(0.0, 1.0, 16) + return {'damping_values': damping_values} + + +if __name__ == "__main__": + test = prior() + run_germany_nuts1_simulation(test['damping_values']) + + # import os + # os.environ["KERAS_BACKEND"] = "tensorflow" + + # import bayesflow as bf + + # simulator = bf.simulators.make_simulator([prior, run_germany_nuts1_simulation]) + # trainings_data = simulator.sample(1) + + # for region in range(16): + # trainings_data[f'region{region}'] = trainings_data[f'region{region}'][:,:, 8][..., np.newaxis] + + # with open('validation_data_16param.pickle', 'wb') as f: + # pickle.dump(trainings_data, f, pickle.HIGHEST_PROTOCOL) + + # with open('trainings_data1_16param.pickle', 'rb') as f: + # trainings_data = pickle.load(f) + # trainings_data['damping_values'] = trainings_data['damping_values'][:, :16] + # for i in range(9): + # with open(f'trainings_data{i+2}_16param.pickle', 'rb') as f: + # data = pickle.load(f) + # data['damping_values'] = data['damping_values'][:, :16] + # trainings_data = {k: np.concatenate([trainings_data[k], data[k]]) for k in trainings_data.keys()} + + # with open('validation_data_16param.pickle', 'rb') as f: + # validation_data = pickle.load(f) + # validation_data['damping_values'] = validation_data['damping_values'][:, :16] + + # adapter = ( + # bf.Adapter() + # .to_array() + # .convert_dtype("float64", "float32") + # .constrain("damping_values", lower=0.0, upper=1.0) + # .rename("damping_values", "inference_variables") + # .concatenate([f'region{i}' for i in range(16)], into="summary_variables", axis=-1) + # .log("summary_variables", p1=True) + # ) + + # summary_network = bf.networks.TimeSeriesNetwork(summary_dim=32) + # inference_network = bf.networks.CouplingFlow() + + # workflow = bf.BasicWorkflow( + # simulator=simulator, + # adapter=adapter, + # summary_network=summary_network, + # inference_network=inference_network, + # standardize='all' + # ) + + # history = workflow.fit_offline(data=trainings_data, epochs=100, batch_size=32, validation_data=validation_data) + + # # workflow.approximator.save(filepath=os.path.join(os.path.dirname(__file__), "model_1params.keras")) + + # plots = workflow.plot_default_diagnostics(test_data=validation_data, calibration_ecdf_kwargs={'difference': True, 'stacked': True}) + # plots['losses'].savefig('losses_couplingflow_16param.png') + # plots['recovery'].savefig('recovery_couplingflow_16param.png') + # plots['calibration_ecdf'].savefig('calibration_ecdf_couplingflow_16param.png') + # plots['z_score_contraction'].savefig('z_score_contraction_couplingflow_16param.png') diff --git a/pycode/examples/simulation/graph_germany_nuts3.py b/pycode/examples/simulation/graph_germany_nuts3.py new file mode 100644 index 0000000000..5153d0f7a7 --- /dev/null +++ b/pycode/examples/simulation/graph_germany_nuts3.py @@ -0,0 +1,412 @@ +############################################################################# +# Copyright (C) 2020-2025 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. +############################################################################# +import numpy as np +import datetime +import os +import memilio.simulation as mio +import memilio.simulation.osecir as osecir +import matplotlib.pyplot as plt + +from enum import Enum +from memilio.simulation.osecir import (Model, Simulation, + interpolate_simulation_result) +import pandas as pd +from memilio.epidata import defaultDict as dd +import pickle + +excluded_ids = [11001, 11002, 11003, 11004, 11005, 11006, 11007, 11008, 11009, 11010, 11011, 11012, 16056] +no_icu_ids = [7338, 9374, 9473, 9573] +region_ids = [region_id for region_id in dd.County.keys() if region_id not in excluded_ids] + +inference_params = ['damping_values', 't_E', 't_ISy', 't_ISev', 'transmission_prob'] + +def plot_region_median_mad( + data: np.ndarray, + region: int, + true_data = None, + ax = None, + label = None, + color = "red" +): + if data.ndim != 3: + raise ValueError("Array not of shape (samples, time_points, regions)") + if true_data is not None: + if true_data.shape != data.shape[1:]: + raise ValueError("True data shape does not match data shape") + n_samples, n_time, n_regions = data.shape + if not (0 <= region < n_regions): + raise IndexError + + x = np.arange(n_time) + vals = data[:, :, region] # (samples, time_points) + med = np.median(vals, axis=0) + mad = np.median(np.abs(vals - med), axis=0) + + if ax is None: + fig, ax = plt.subplots() + + line, = ax.plot(x, med, lw=2, label=label or f"Region {region}", color=color) + band = ax.fill_between(x, med - mad, med + mad, alpha=0.25, color=color) + if true_data is not None: + true_vals = true_data[:, region] # (time_points,) + ax.plot(x, true_vals, lw=2, color="black", label="True data") + + ax.set_xlabel("Time") + ax.set_ylabel("ICU") + ax.set_title(f"Region {region}") + if label is not None: + ax.legend() + return line, band + + +def plot_aggregated_over_regions( + data: np.ndarray, + region_agg = np.sum, + true_data = None, + ax = None, + label = None, + color = 'red' +): + if data.ndim != 3: + raise ValueError("Array not of shape (samples, time_points, regions)") + if true_data is not None: + if true_data.shape != data.shape[1:]: + raise ValueError("True data shape does not match data shape") + + # Aggregate over regions + agg_over_regions = region_agg(data, axis=-1) # (samples, time_points) + + # Aggregate over samples + agg_median = np.median(agg_over_regions, axis=0) # (time_points, ) + agg_mad = np.median( + np.abs(agg_over_regions - agg_median[None]), + axis=0 + ) + + x = np.arange(agg_median.shape[0]) + if ax is None: + fig, ax = plt.subplots() + + line, = ax.plot(x, agg_median, lw=2, label=label or "Aggregated over regions", color=color) + band = ax.fill_between( + x, + agg_median - agg_mad, + agg_median + agg_mad, + alpha=0.25, + color=color + ) + if true_data is not None: + true_vals = region_agg(true_data, axis=-1) # (time_points,) + ax.plot(x, true_vals, lw=2, color="black", label="True data") + + ax.set_xlabel("Time") + ax.set_ylabel("ICU") + if label is not None: + ax.legend() + + return line, band + +class Simulation: + """ """ + + def __init__(self, data_dir, start_date, results_dir): + self.num_groups = 1 + self.data_dir = data_dir + self.start_date = start_date + self.results_dir = results_dir + if not os.path.exists(self.results_dir): + os.makedirs(self.results_dir) + + def set_covid_parameters(self, model, t_E, t_ISy, t_ISev, t_Cr, transmission_prob): + """ + + :param model: + + """ + model.parameters.TimeExposed[mio.AgeGroup(0)] = t_E + model.parameters.TimeInfectedNoSymptoms[mio.AgeGroup(0)] = 5.2 - t_E + model.parameters.TimeInfectedSymptoms[mio.AgeGroup(0)] = t_ISy + model.parameters.TimeInfectedSevere[mio.AgeGroup(0)] = t_ISev + model.parameters.TimeInfectedCritical[mio.AgeGroup(0)] = t_Cr + + # probabilities + model.parameters.TransmissionProbabilityOnContact[mio.AgeGroup(0)] = transmission_prob + model.parameters.RelativeTransmissionNoSymptoms[mio.AgeGroup(0)] = 1 + + model.parameters.RecoveredPerInfectedNoSymptoms[mio.AgeGroup(0)] = 0.2069 + model.parameters.SeverePerInfectedSymptoms[mio.AgeGroup(0)] = 0.07864 + model.parameters.CriticalPerSevere[mio.AgeGroup(0)] = 0.17318 + model.parameters.DeathsPerCritical[mio.AgeGroup(0)] = 0.21718 + + # start day is set to the n-th day of the year + model.parameters.StartDay = self.start_date.timetuple().tm_yday + + model.parameters.Seasonality = mio.UncertainValue(0.2) + + def set_contact_matrices(self, model): + """ + + :param model: + + """ + contact_matrices = mio.ContactMatrixGroup(1, self.num_groups) + + baseline = np.ones((self.num_groups, self.num_groups)) * 7.95 + minimum = np.zeros((self.num_groups, self.num_groups)) + contact_matrices[0] = mio.ContactMatrix(baseline, minimum) + model.parameters.ContactPatterns.cont_freq_mat = contact_matrices + + def set_npis(self, params, end_date, damping_value): + """ + + :param params: + :param end_date: + + """ + + start_damping = datetime.date( + year=2020, month=10, day=8) + + if start_damping < end_date: + start_date = (start_damping - self.start_date).days + params.ContactPatterns.cont_freq_mat[0].add_damping(mio.Damping(np.r_[damping_value], t=start_date)) + + def get_graph(self, end_date, t_E, t_ISy, t_ISev, t_Cr, transmission_prob): + """ + + :param end_date: + + """ + print("Initializing model...") + model = Model(self.num_groups) + self.set_covid_parameters(model, t_E, t_ISy, t_ISev, t_Cr, transmission_prob) + self.set_contact_matrices(model) + print("Model initialized.") + + graph = osecir.ModelGraph() + + scaling_factor_infected = [2.5] + scaling_factor_icu = 1.0 + tnt_capacity_factor = 7.5 / 100000. + + data_dir_Germany = os.path.join(self.data_dir, "Germany") + mobility_data_file = os.path.join( + data_dir_Germany, "mobility", "commuter_mobility_2022.txt") + pydata_dir = os.path.join(data_dir_Germany, "pydata") + + path_population_data = os.path.join(pydata_dir, + "county_current_population.json") + + print("Setting nodes...") + mio.osecir.set_nodes( + model.parameters, + mio.Date(self.start_date.year, + self.start_date.month, self.start_date.day), + mio.Date(end_date.year, + end_date.month, end_date.day), pydata_dir, + path_population_data, True, graph, scaling_factor_infected, + scaling_factor_icu, 0, 0, False) + + print("Setting edges...") + mio.osecir.set_edges( + mobility_data_file, graph, 1) + + print("Graph created.") + + return graph + + def run(self, num_days_sim, damping_values, t_E, t_ISy, t_ISev, t_Cr, transmission_prob, save_graph=True): + """ + + :param num_days_sim: + :param num_runs: (Default value = 10) + :param save_graph: (Default value = True) + :param create_gif: (Default value = True) + + """ + mio.set_log_level(mio.LogLevel.Warning) + end_date = self.start_date + datetime.timedelta(days=num_days_sim) + + graph = self.get_graph(end_date, t_E, t_ISy, t_ISev, t_Cr, transmission_prob) + + mobility_graph = osecir.MobilityGraph() + for node_idx in range(graph.num_nodes): + node = graph.get_node(node_idx) + self.set_npis(node.property.parameters, end_date, damping_values[node.id // 1000 -1]) + mobility_graph.add_node(node.id, node.property) + for edge_idx in range(graph.num_edges): + mobility_graph.add_edge( + graph.get_edge(edge_idx).start_node_idx, + graph.get_edge(edge_idx).end_node_idx, + graph.get_edge(edge_idx).property) + mobility_sim = osecir.MobilitySimulation(mobility_graph, t0=0, dt=0.5) + mobility_sim.advance(num_days_sim) + + results = {} + for node_idx in range(mobility_sim.graph.num_nodes): + node = mobility_sim.graph.get_node(node_idx) + if node.id in no_icu_ids: + results[f'no_icu_region{node_idx}'] = osecir.interpolate_simulation_result(node.property.result) + else: + results[f'region{node_idx}'] = osecir.interpolate_simulation_result(node.property.result) + + return results + +def run_germany_nuts3_simulation(damping_values, t_E, t_ISy, t_ISev, t_Cr, transmission_prob): + mio.set_log_level(mio.LogLevel.Warning) + file_path = os.path.dirname(os.path.abspath(__file__)) + + sim = Simulation( + data_dir=os.path.join(file_path, "../../../data"), + start_date=datetime.date(year=2020, month=10, day=1), + results_dir=os.path.join(file_path, "../../../results_osecir")) + num_days_sim = 60 + + results = sim.run(num_days_sim, damping_values, t_E, t_ISy, t_ISev, t_Cr, transmission_prob) + + return results + +def prior(): + damping_values = np.random.uniform(0.0, 1.0, 16) + t_E = np.random.uniform(1., 5.2) + t_ISy = np.random.uniform(4., 10.) + t_ISev = np.random.uniform(5., 10.) + t_Cr = np.random.uniform(9., 17.) + transmission_prob = np.random.uniform(0., 0.2) + return {'damping_values': damping_values, + 't_E': t_E, + 't_ISy': t_ISy, + 't_ISev': t_ISev, + 't_Cr': t_Cr, + 'transmission_prob': transmission_prob} + +def load_divi_data(): + file_path = os.path.dirname(os.path.abspath(__file__)) + divi_path = os.path.join(file_path, "../../../data/Germany/pydata") + + data = pd.read_json(os.path.join(divi_path, "county_divi_all_dates.json")) + data = data[data['Date']>= np.datetime64(datetime.date(2020, 10, 1))] + data = data[data['Date'] <= np.datetime64(datetime.date(2020, 10, 1) + datetime.timedelta(days=50))] + data = data.sort_values(by=['ID_County', 'Date']) + divi_data = data.pivot(index='Date', columns='ID_County', values='ICU') + divi_dict = {f"region{i}": divi_data[region_id].to_numpy()[None, :, None] for i, region_id in enumerate(region_ids) if region_id not in no_icu_ids} + + return divi_data.to_numpy(), divi_dict + + +if __name__ == "__main__": + + import os + os.environ["KERAS_BACKEND"] = "tensorflow" + + import bayesflow as bf + from tensorflow import keras + + simulator = bf.simulators.make_simulator([prior, run_germany_nuts3_simulation]) + trainings_data = simulator.sample(1000) + + for key in trainings_data.keys(): + if key not in inference_params: + trainings_data[key] = trainings_data[key][:, :, 7][..., np.newaxis] + + with open('trainings_data10_counties_wcovidparams_oct.pickle', 'wb') as f: + pickle.dump(trainings_data, f, pickle.HIGHEST_PROTOCOL) + + # with open('trainings_data1_counties_wcovidparams_oct.pickle', 'rb') as f: + # trainings_data = pickle.load(f) + # trainings_data = {k: np.round(v) if ('region' in k) else v for k, v in trainings_data.items()} + # for i in range(19): + # with open(f'trainings_data{i+2}_counties_wcovidparams_oct.pickle', 'rb') as f: + # data = pickle.load(f) + # trainings_data = {k: np.concatenate([trainings_data[k], np.round(data[k])]) if ('region' in k) else np.concatenate([trainings_data[k], data[k]]) for k in trainings_data.keys()} + + # with open('validation_data_counties_wcovidparams_oct.pickle', 'rb') as f: + # validation_data = pickle.load(f) + # divi_dict = {k: np.round(v) if ('region' in k) else v for k, v in validation_data.items()} + + # adapter = ( + # bf.Adapter() + # .to_array() + # .convert_dtype("float64", "float32") + # .constrain("damping_values", lower=0.0, upper=1.0) + # .constrain("t_E", lower=1.0, upper=6.0) + # .constrain("t_ISy", lower=5.0, upper=10.0) + # .constrain("t_ISev", lower=2.0, upper=8.0) + # .concatenate(["damping_values", "t_E", "t_ISy", "t_ISev"], into="inference_variables", axis=-1) + # .concatenate([f'region{i}' for i in range(len(region_ids)) if region_ids[i] not in no_icu_ids], into="summary_variables", axis=-1) + # .log("summary_variables", p1=True) + # ) + + # print("summary_variables shape:", adapter(trainings_data)["summary_variables"].shape) + # print("inference_variables shape:", adapter(trainings_data)["inference_variables"].shape) + + # summary_network = bf.networks.TimeSeriesNetwork(summary_dim=38) + # inference_network = bf.networks.CouplingFlow(depth=7, transform='spline') + + # workflow = bf.BasicWorkflow( + # simulator=simulator, + # adapter=adapter, + # summary_network=summary_network, + # inference_network=inference_network, + # standardize='all' + # ) + + # history = workflow.fit_offline(data=trainings_data, epochs=100, batch_size=32, validation_data=validation_data) + + # workflow.approximator.save(filepath=os.path.join(os.path.dirname(__file__), "model_countylvl_wcovidparams_oct.keras")) + + # plots = workflow.plot_default_diagnostics(test_data=validation_data, calibration_ecdf_kwargs={'difference': True, 'stacked': True}) + # plots['losses'].savefig('losses_countylvl_wcovidparams2_oct.png') + # plots['recovery'].savefig('recovery_countylvl_wcovidparams2_oct.png') + # plots['calibration_ecdf'].savefig('calibration_ecdf_countylvl_wcovidparams2_oct.png') + # plots['z_score_contraction'].savefig('z_score_contraction_countylvl_wcovidparams2_oct.png') + + # divi_data, divi_dict = load_divi_data() + # # divi_data = np.concatenate( + # # [validation_data[f'region{i}'] for i in range(len(region_ids)) if region_ids[i] not in no_icu_ids], + # # axis=-1 + # # ) + # workflow.approximator = keras.models.load_model(os.path.join(os.path.dirname(__file__), "model_countylvl_wcovidparams_oct.keras")) + + # samples = workflow.sample(conditions=divi_dict, num_samples=1000) + # samples = np.concatenate([samples[key] for key in inference_params], axis=-1) + # samples = np.squeeze(samples) + # sims = [] + # for i in range(samples.shape[0]): + # result = run_germany_nuts3_simulation(samples[i][:16], *samples[i][16:]) + # for key in result.keys(): + # result[key] = np.array(result[key])[:, 7, None] + # sims.append(np.concatenate([result[key] for key in result.keys() if key.startswith('region')], axis=-1)) + # sims = np.array(sims) + # sims = np.floor(sims) + + # np.random.seed(42) + # fig, ax = plt.subplots(nrows=2, ncols=5, figsize=(12, 5), layout="constrained") + # ax = ax.flatten() + # rand_index = np.random.choice(sims.shape[-1], replace=False, size=len(ax)) + # for i, a in enumerate(ax): + # plot_region_median_mad(sims, region=rand_index[i], true_data=divi_data, label=r"Median $\pm$ Mad", ax=a) + # plt.savefig('random_regions_wcovidparams_oct.png') + # # plt.show() + # #%% + # plot_aggregated_over_regions(sims, true_data=divi_data, label="Region Aggregated Median $\pm$ Mad") + # plt.savefig('region_aggregated_wcovidparams_oct.png') + # # plt.show() + # # %% \ No newline at end of file diff --git a/pycode/examples/simulation/graph_spain_nuts3.py b/pycode/examples/simulation/graph_spain_nuts3.py new file mode 100644 index 0000000000..081fb1ea22 --- /dev/null +++ b/pycode/examples/simulation/graph_spain_nuts3.py @@ -0,0 +1,266 @@ +############################################################################# +# Copyright (C) 2020-2025 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. +############################################################################# +import numpy as np +import datetime +import os +import memilio.simulation as mio +import memilio.simulation.osecir as osecir +import matplotlib.pyplot as plt + +from enum import Enum +from memilio.simulation.osecir import (Model, Simulation, + interpolate_simulation_result) + +import pickle + + +class Simulation: + """ """ + + def __init__(self, data_dir, start_date, results_dir): + self.num_groups = 1 + self.data_dir = data_dir + self.start_date = start_date + self.results_dir = results_dir + if not os.path.exists(self.results_dir): + os.makedirs(self.results_dir) + + def set_covid_parameters(self, model): + """ + + :param model: + + """ + model.parameters.TimeExposed[mio.AgeGroup(0)] = 3.335 + model.parameters.TimeInfectedNoSymptoms[mio.AgeGroup(0)] = 2.58916 + model.parameters.TimeInfectedSymptoms[mio.AgeGroup(0)] = 6.94547 + model.parameters.TimeInfectedSevere[mio.AgeGroup(0)] = 7.28196 + model.parameters.TimeInfectedCritical[mio.AgeGroup(0)] = 13.066 + + # probabilities + model.parameters.TransmissionProbabilityOnContact[mio.AgeGroup( + 0)] = 0.07333 + model.parameters.RelativeTransmissionNoSymptoms[mio.AgeGroup(0)] = 1 + + model.parameters.RecoveredPerInfectedNoSymptoms[mio.AgeGroup( + 0)] = 0.2069 + model.parameters.SeverePerInfectedSymptoms[mio.AgeGroup(0)] = 0.07864 + model.parameters.CriticalPerSevere[mio.AgeGroup(0)] = 0.17318 + model.parameters.DeathsPerCritical[mio.AgeGroup(0)] = 0.21718 + + # start day is set to the n-th day of the year + model.parameters.StartDay = self.start_date.timetuple().tm_yday + + model.parameters.Seasonality = mio.UncertainValue(0.2) + + def set_contact_matrices(self, model): + """ + + :param model: + + """ + contact_matrices = mio.ContactMatrixGroup(1, self.num_groups) + + baseline = np.ones((self.num_groups, self.num_groups)) * 12.32 + minimum = np.ones((self.num_groups, self.num_groups)) * 0 + contact_matrices[0] = mio.ContactMatrix(baseline, minimum) + model.parameters.ContactPatterns.cont_freq_mat = contact_matrices + + def set_npis(self, params, end_date, damping_value): + """ + + :param params: + :param end_date: + + """ + + start_damping = datetime.date( + 2020, 12, 18) + + if start_damping < end_date: + start_date = (start_damping - self.start_date).days + params.ContactPatterns.cont_freq_mat[0].add_damping( + mio.Damping(np.r_[damping_value], t=start_date)) + + def get_graph(self, end_date): + """ + + :param end_date: + + """ + print("Initializing model...") + model = Model(self.num_groups) + self.set_covid_parameters(model) + self.set_contact_matrices(model) + print("Model initialized.") + + graph = osecir.ModelGraph() + + scaling_factor_infected = [2.5] + scaling_factor_icu = 1.0 + tnt_capacity_factor = 7.5 / 100000. + + data_dir_Spain = os.path.join(self.data_dir, "Spain") + mobility_data_file = os.path.join( + data_dir_Spain, "mobility", "commuter_mobility_2022.txt") + pydata_dir = os.path.join(data_dir_Spain, "pydata") + + path_population_data = os.path.join( + pydata_dir, "provincias_current_population.json") + + print("Setting nodes...") + mio.osecir.set_nodes_provincias( + model.parameters, + mio.Date(self.start_date.year, + self.start_date.month, self.start_date.day), + mio.Date(end_date.year, + end_date.month, end_date.day), pydata_dir, + path_population_data, True, graph, scaling_factor_infected, + scaling_factor_icu, 0, 0, False, False) + + print("Setting edges...") + mio.osecir.set_edges( + mobility_data_file, graph, 1) + + print("Graph created.") + + return graph + + def run(self, num_days_sim, damping_values, save_graph=True): + """ + + :param num_days_sim: + :param num_runs: (Default value = 10) + :param save_graph: (Default value = True) + :param create_gif: (Default value = True) + + """ + mio.set_log_level(mio.LogLevel.Warning) + end_date = self.start_date + datetime.timedelta(days=num_days_sim) + + graph = self.get_graph(end_date) + + if save_graph: + path_graph = os.path.join(self.results_dir, "graph") + if not os.path.exists(path_graph): + os.makedirs(path_graph) + osecir.write_graph(graph, path_graph) + + mobility_graph = osecir.MobilityGraph() + for node_idx in range(graph.num_nodes): + node = graph.get_node(node_idx) + self.set_npis(node.property.parameters, + end_date, damping_values[node_idx]) + mobility_graph.add_node(node.id, node.property) + for edge_idx in range(graph.num_edges): + mobility_graph.add_edge( + graph.get_edge(edge_idx).start_node_idx, + graph.get_edge(edge_idx).end_node_idx, + graph.get_edge(edge_idx).property) + mobility_sim = osecir.MobilitySimulation(mobility_graph, t0=0, dt=0.5) + mobility_sim.advance(num_days_sim) + + results = [] + for node_idx in range(graph.num_nodes): + results.append(osecir.interpolate_simulation_result( + mobility_sim.graph.get_node(node_idx).property.result)) + + return results + + +def run_spain_nuts3_simulation(damping_values): + mio.set_log_level(mio.LogLevel.Warning) + file_path = os.path.dirname(os.path.abspath(__file__)) + + sim = Simulation( + data_dir=os.path.join(file_path, "../../../data"), + start_date=datetime.date(year=2020, month=12, day=12), + results_dir=os.path.join(file_path, "../../../results_osecir")) + num_days_sim = 50 + + results = sim.run(num_days_sim, damping_values) + + return {f'region{region}': results[region] for region in range(len(results))} + + +def prior(): + damping_values = np.random.uniform(0.0, 1.0, 47) + return {'damping_values': damping_values} + + +if __name__ == "__main__": + test = prior() + run_spain_nuts3_simulation(test['damping_values']) + # import os + # os.environ["KERAS_BACKEND"] = "tensorflow" + + # import bayesflow as bf + + # simulator = bf.simulators.make_simulator([prior, run_germany_nuts3_simulation]) + # # trainings_data = simulator.sample(1000) + + # # for region in range(400): + # # trainings_data[f'region{region}'] = trainings_data[f'region{region}'][:,:, 8][..., np.newaxis] + + # # with open('validation_data_400params.pickle', 'wb') as f: + # # pickle.dump(trainings_data, f, pickle.HIGHEST_PROTOCOL) + + # with open('trainings_data1_400params.pickle', 'rb') as f: + # trainings_data = pickle.load(f) + # for i in range(9): + # with open(f'trainings_data{i+2}_400params.pickle', 'rb') as f: + # data = pickle.load(f) + # trainings_data = {k: np.concatenate([trainings_data[k], data[k]]) for k in trainings_data.keys()} + + # with open('validation_data_400params.pickle', 'rb') as f: + # validation_data = pickle.load(f) + + # adapter = ( + # bf.Adapter() + # .to_array() + # .convert_dtype("float64", "float32") + # .constrain("damping_values", lower=0.0, upper=1.0) + # .rename("damping_values", "inference_variables") + # .concatenate([f'region{i}' for i in range(400)], into="summary_variables", axis=-1) + # .log("summary_variables", p1=True) + # ) + + # print("summary_variables shape:", adapter(trainings_data)["summary_variables"].shape) + + # summary_network = bf.networks.TimeSeriesNetwork(summary_dim=700, recurrent_dim=256) + # inference_network = bf.networks.DiffusionModel(subnet_kwargs={'widths': {512, 512, 512, 512, 512}}) + + # workflow = bf.BasicWorkflow( + # simulator=simulator, + # adapter=adapter, + # summary_network=summary_network, + # inference_network=inference_network, + # standardize='all' + # ) + + # history = workflow.fit_offline(data=trainings_data, epochs=1000, batch_size=32, validation_data=validation_data) + + # # workflow.approximator.save(filepath=os.path.join(os.path.dirname(__file__), "model_10params.keras")) + + # plots = workflow.plot_default_diagnostics(test_data=validation_data, calibration_ecdf_kwargs={'difference': True, 'stacked': True}) + # plots['losses'].savefig('losses_diffusionmodel_400params.png') + # plots['recovery'].savefig('recovery_diffusionmodel_400params.png') + # plots['calibration_ecdf'].savefig('calibration_ecdf_diffusionmodel_400params.png') + # plots['z_score_contraction'].savefig('z_score_contraction_diffusionmodel_400params.png') diff --git a/pycode/memilio-epidata/memilio/epidata/defaultDict.py b/pycode/memilio-epidata/memilio/epidata/defaultDict.py index 73389b82d2..281f88c130 100644 --- a/pycode/memilio-epidata/memilio/epidata/defaultDict.py +++ b/pycode/memilio-epidata/memilio/epidata/defaultDict.py @@ -45,9 +45,9 @@ 'out_folder': default_file_path, 'update_data': False, 'start_date': date(2020, 1, 1), - 'end_date': date.today(), + 'end_date': date(2021, 1, 1), 'split_berlin': False, - 'impute_dates': False, + 'impute_dates': True, 'moving_average': 0, 'file_format': 'json_timeasstring', 'no_raw': False, @@ -704,3 +704,111 @@ def invert_dict(dict_to_invert): """ return {val: key for key, val in dict_to_invert.items()} + + +Provincias = {2: 'Araba/Álava', + 3: 'Albacete', + 4: 'Alicante/Alacant', + 5: 'Almería', + 6: 'Ávila', + 7: 'Badajoz', + 8: 'Balears, Illes', + 9: 'Barcelona', + 10: 'Burgos', + 11: 'Cáceres', + 12: 'Cádiz', + 13: 'Castellón/Castelló', + 14: 'Ciudad Real', + 15: 'Córdoba', + 16: 'Coruña, A', + 17: 'Cuenca', + 18: 'Girona', + 19: 'Granada', + 20: 'Guadalajara', + 21: 'Gipuzkoa', + 22: 'Huelva', + 23: 'Huesca', + 24: 'Jaén', + 25: 'León', + 26: 'Lleida', + 27: 'Rioja, La', + 28: 'Lugo', + 29: 'Madrid', + 30: 'Málaga', + 31: 'Murcia', + 32: 'Navarra', + 33: 'Asturias', + 34: 'Palencia', + 35: 'Palmas, Las', + 36: 'Pontevedra', + 37: 'Salamanca', + 38: 'Santa Cruz de Tenerife', + 39: 'Cantabria', + 40: 'Segovia', + 41: 'Sevilla', + 42: 'Soria', + 43: 'Tarragona', + 44: 'Teruel', + 45: 'Toledo', + 46: 'Valencia/València', + 47: 'Valladolid', + 48: 'Bizkaia', + 49: 'Zamora', + 50: 'Zaragoza', + 51: 'Ceuta', + 52: 'Melilla', + 53: 'Ourense'} + + +Provincia_ISO_to_ID = {'A': 4, + 'AB': 3, + 'AL': 5, + 'AV': 6, + 'B': 9, + 'BA': 7, + 'BI': 48, + 'BU': 10, + 'C': 16, + 'CA': 12, + 'CC': 11, + 'CE': 51, + 'CO': 15, + 'CR': 14, + 'CS': 13, + 'CU': 17, + 'GC': 35, + 'GI': 18, + 'GR': 19, + 'GU': 20, + 'H': 22, + 'HU': 23, + 'J': 24, + 'L': 26, + 'LE': 25, + 'LO': 27, + 'LU': 28, + 'M': 29, + 'MA': 30, + 'ML': 52, + 'MU': 31, + 'NA': 32, + 'O': 33, + 'OR': 53, + 'P': 34, + 'PM': 8, + 'PO': 36, + 'SA': 37, + 'S': 39, + 'SE': 41, + 'SG': 40, + 'SO': 42, + 'SS': 21, + 'T': 43, + 'TE': 44, + 'TF': 38, + 'TO': 45, + 'V': 46, + 'VA': 47, + 'VI': 2, + 'ZA': 48, + 'Z': 50} diff --git a/pycode/memilio-epidata/memilio/epidata/getPopulationData.py b/pycode/memilio-epidata/memilio/epidata/getPopulationData.py index c7d0b401de..07e8ab7415 100644 --- a/pycode/memilio-epidata/memilio/epidata/getPopulationData.py +++ b/pycode/memilio-epidata/memilio/epidata/getPopulationData.py @@ -55,8 +55,9 @@ def read_population_data(ref_year): req = requests.get(download_url) df_pop_raw = pd.read_csv(io.StringIO(req.text), sep=';', header=5) except pd.errors.ParserError: - gd.default_print('Warning', 'Data for year '+str(ref_year) + - ' is not available; downloading newest data instead.') + gd.default_print( + 'Warning', 'Data for year ' + str(ref_year) + + ' is not available; downloading newest data instead.') ref_year = None if ref_year is None: download_url = 'https://www.regionalstatistik.de/genesis/online?operation=download&code=12411-02-03-4&option=csv' @@ -66,7 +67,9 @@ def read_population_data(ref_year): return df_pop_raw, ref_year -def export_population_dataframe(df_pop: pd.DataFrame, directory: str, file_format: str, merge_eisenach: bool, ref_year): +def export_population_dataframe( + df_pop: pd.DataFrame, directory: str, file_format: str, + merge_eisenach: bool, ref_year): """ Writes population dataframe into directory with new column names and age groups :param df_pop: Population data DataFrame to be exported pd.DataFrame @@ -140,6 +143,10 @@ def export_population_dataframe(df_pop: pd.DataFrame, directory: str, file_forma columns=dd.EngEng["idCounty"]) gd.write_dataframe(df_pop_export, directory, filename, file_format) + gd.write_dataframe(aggregate_to_state_level(df_pop_export), + directory, filename + '_states', file_format) + gd.write_dataframe(aggregate_to_country_level(df_pop_export), + directory, filename + '_germany', file_format) return df_pop_export @@ -190,8 +197,8 @@ def assign_population_data(df_pop_raw, counties, age_cols, idCounty_idx): # direct assignment of population data found df_pop.loc[df_pop[dd.EngEng['idCounty']] == df_pop_raw.loc [start_idx, dd.EngEng['idCounty']], - age_cols] = df_pop_raw.loc[start_idx: start_idx + num_age_groups - 1, dd.EngEng - ['number']].values.astype(int) + age_cols] = df_pop_raw.loc[start_idx: start_idx + + num_age_groups - 1, dd.EngEng['number']].values.astype(int) # Berlin and Hamburg elif county_id + '000' in counties[:, 1]: # direct assignment of population data found @@ -262,7 +269,8 @@ def fetch_population_data(read_data: bool = dd.defaultDict['read_data'], if read_data == True: gd.default_print( - 'Warning', 'Read_data is not supportet for getPopulationData.py. Setting read_data = False') + 'Warning', + 'Read_data is not supportet for getPopulationData.py. Setting read_data = False') read_data = False directory = os.path.join(out_folder, 'Germany', 'pydata') @@ -444,6 +452,28 @@ def get_population_data(read_data: bool = dd.defaultDict['read_data'], return df_pop_export +def aggregate_to_state_level(df_pop: pd.DataFrame): + + countyIDtostateID = geoger.get_countyid_to_stateid_map() + + df_pop['ID_State'] = df_pop[dd.EngEng['idCounty']].map(countyIDtostateID) + df_pop = df_pop.drop( + columns='ID_County').groupby( + 'ID_State', as_index=True).sum() + df_pop['ID_State'] = df_pop.index + return df_pop + + +def aggregate_to_country_level(df_pop: pd.DataFrame): + + df_pop['ID_Country'] = 0 + df_pop = df_pop.drop( + columns=['ID_County', 'ID_State']).groupby( + 'ID_Country', as_index=True).sum() + df_pop['ID_Country'] = df_pop.index + return df_pop + + def main(): """ Main program entry.""" diff --git a/pycode/memilio-epidata/memilio/epidata/getSimulationDataSpain.py b/pycode/memilio-epidata/memilio/epidata/getSimulationDataSpain.py new file mode 100644 index 0000000000..8ccde09fc5 --- /dev/null +++ b/pycode/memilio-epidata/memilio/epidata/getSimulationDataSpain.py @@ -0,0 +1,202 @@ +import pandas as pd +import os +import io +import requests + +from pyspainmobility import Mobility, Zones +import defaultDict as dd + +from memilio.epidata import getDataIntoPandasDataFrame as gd + + +def fetch_population_data(): + download_url = 'https://servicios.ine.es/wstempus/js/es/DATOS_TABLA/67988?tip=AM&' + req = requests.get(download_url) + req.encoding = 'ISO-8859-1' + + df = pd.read_json(io.StringIO(req.text)) + df = df[['MetaData', 'Data']] + df = df[df['MetaData'].apply( + lambda x: x[0]['T3_Variable'] == 'Provincias')] + df = df[df['MetaData'].apply( + lambda x: x[1]['Nombre'] == 'Total')] + df['ID_Provincia'] = df['MetaData'].apply(lambda x: x[0]['Id']) + df['Population'] = df['Data'].apply(lambda x: x[0]['Valor']) + return df[['ID_Provincia', 'Population']] + + +def remove_islands(df, column_labels=['ID_Provincia']): + for label in column_labels: + df = df[~df[label].isin([51, 52, 8, 35, 38])] + return df + + +def get_population_data(): + df = fetch_population_data() + df = remove_islands(df) + + return df + + +def fetch_icu_data(): + download_url = 'https://www.sanidad.gob.es/areas/alertasEmergenciasSanitarias/alertasActuales/nCov/documentos/Datos_Capacidad_Asistencial_Historico_14072023.csv' + req = requests.get(download_url) + req.encoding = 'ISO-8859-1' + + df = pd.read_csv(io.StringIO(req.text), sep=';') + return df + + +def preprocess_icu_data(df): + df_icu = df[df["Unidad"] == "U. Críticas SIN respirador"] + df_icu_vent = df[df["Unidad"] == "U. Críticas CON respirador"] + + df_icu = df_icu[['Fecha', 'ID_Provincia', 'OCUPADAS_COVID19']].rename( + columns={'OCUPADAS_COVID19': 'ICU'}) + df_icu_vent = df_icu_vent[['Fecha', 'ID_Provincia', 'OCUPADAS_COVID19']].rename( + columns={'OCUPADAS_COVID19': 'ICU_ventilated'}) + + df_merged = pd.merge(df_icu, df_icu_vent, on=[ + 'Fecha', 'ID_Provincia'], how='outer') + df_merged['Fecha'] = pd.to_datetime( + df_merged['Fecha'], format='%d/%m/%Y').dt.strftime('%Y-%m-%d') + df_merged.rename(columns={'Fecha': 'Date'}, inplace=True) + + return df_merged + + +def get_icu_data(): + df = fetch_icu_data() + df.rename(columns={'Cod_Provincia': 'ID_Provincia'}, inplace=True) + df = remove_islands(df) + df = preprocess_icu_data(df) + + return df + + +def fetch_case_data(): + download_url = 'https://cnecovid.isciii.es/covid19/resources/casos_diagnostico_provincia.csv' + req = requests.get(download_url) + + df = pd.read_csv(io.StringIO(req.text), sep=',') + + return df + + +def preprocess_case_data(df): + df['provincia_iso'] = df['provincia_iso'].map(dd.Provincia_ISO_to_ID) + df = df.rename( + columns={'provincia_iso': 'ID_County', 'fecha': 'Date', + 'num_casos': 'Confirmed'})[ + ['ID_County', 'Date', 'Confirmed']] + # ensure correct types + df['ID_County'] = pd.to_numeric( + df['ID_County'], errors='coerce').astype('Int64') + df = df[df['ID_County'].notna()].copy() + df['ID_County'] = df['ID_County'].astype(int) + return df + + +def get_case_data(): + df_raw = fetch_case_data() + df = preprocess_case_data(df_raw) + df = df.groupby(['ID_County', 'Date'], as_index=False)['Confirmed'].sum() + df = df.sort_values(['ID_County', 'Date']) + df['Confirmed'] = df.groupby('ID_County')['Confirmed'].cumsum() + return df + + +def download_mobility_data(data_dir, start_date, end_date, level): + mobility_data = Mobility(version=2, zones=level, + start_date=start_date, end_date=end_date, output_directory=data_dir) + mobility_data.get_od_data() + + +def preprocess_mobility_data(df, data_dir): + df.drop(columns=['hour', 'trips_total_length_km'], inplace=True) + if not os.path.exists(os.path.join(data_dir, 'poblacion.csv')): + download_url = 'https://movilidad-opendata.mitma.es/zonificacion/poblacion.csv' + req = requests.get(download_url) + with open(os.path.join(data_dir, 'poblacion.csv'), 'wb') as f: + f.write(req.content) + + zonification_df = pd.read_csv(os.path.join(data_dir, 'poblacion.csv'), sep='|')[ + ['municipio', 'provincia', 'poblacion']] + poblacion = zonification_df.groupby('provincia')[ + 'poblacion'].sum() + zonification_df.drop_duplicates( + subset=['municipio', 'provincia'], inplace=True) + municipio_to_provincia = dict( + zip(zonification_df['municipio'], zonification_df['provincia'])) + df['id_origin'] = df['id_origin'].map(municipio_to_provincia) + df['id_destination'] = df['id_destination'].map(municipio_to_provincia) + + df.query('id_origin != id_destination', inplace=True) + df = df.groupby(['date', 'id_origin', 'id_destination'], + as_index=False)['n_trips'].sum() + + df['n_trips'] = df.apply( + lambda row: row['n_trips'] / poblacion[row['id_origin']], axis=1) + + df = df.groupby(['id_origin', 'id_destination'], + as_index=False)['n_trips'].mean() + + df = remove_islands(df, ['id_origin', 'id_destination']) + + return df + + +def get_mobility_data(data_dir, start_date='2022-08-01', end_date='2022-08-31', level='municipios'): + filename = f'Viajes_{level}_{start_date}_{end_date}_v2.parquet' + if not os.path.exists(os.path.join(data_dir, filename)): + print( + f"File {os.path.exists(os.path.join(data_dir, filename))} does not exist. Downloading mobility data...") + download_mobility_data(data_dir, start_date, end_date, level) + + df = pd.read_parquet(os.path.join(data_dir, filename)) + df = preprocess_mobility_data(df, data_dir) + + return df + + +if __name__ == "__main__": + + data_dir = os.path.join(os.path.dirname( + os.path.abspath(__file__)), "../../../../data/Spain") + pydata_dir = os.path.join(data_dir, 'pydata') + os.makedirs(pydata_dir, exist_ok=True) + mobility_dir = os.path.join(data_dir, 'mobility') + os.makedirs(mobility_dir, exist_ok=True) + + df = get_population_data() + df.to_json(os.path.join( + pydata_dir, 'provincias_current_population.json'), orient='records') + + df = get_icu_data() + # rename to match expected column name in parameters_io.h + if 'ID_Provincia' in df.columns: + df.rename(columns={'ID_Provincia': 'ID_County'}, inplace=True) + # and should also be int + if 'ID_County' in df.columns: + df['ID_County'] = pd.to_numeric(df['ID_County'], errors='coerce') + df = df[df['ID_County'].notna()].copy() + df['ID_County'] = df['ID_County'].astype(int) + df.to_json(os.path.join(pydata_dir, 'provincia_icu.json'), orient='records') + + df = get_case_data() + # same for case data + if 'ID_Provincia' in df.columns: + df.rename(columns={'ID_Provincia': 'ID_County'}, inplace=True) + if 'ID_County' in df.columns: + df['ID_County'] = pd.to_numeric(df['ID_County'], errors='coerce') + df = df[df['ID_County'].notna()].copy() + df['ID_County'] = df['ID_County'].astype(int) + df.to_json(os.path.join( + pydata_dir, 'cases_all_pronvincias.json'), orient='records') + + df = get_mobility_data(mobility_dir) + matrix = df.pivot(index='id_origin', + columns='id_destination', values='n_trips').fillna(0) + + gd.write_dataframe(matrix, data_dir, 'commuter_mobility', 'txt', { + 'sep': ' ', 'index': False, 'header': False}) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index 50ea9b65b8..b0d82a94e9 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp @@ -200,7 +200,13 @@ PYBIND11_MODULE(_simulation_osecir, m) mio::osecir::ParametersBase>(m, "Parameters") .def(py::init()) .def("check_constraints", &mio::osecir::Parameters::check_constraints) - .def("apply_constraints", &mio::osecir::Parameters::apply_constraints); + .def("apply_constraints", &mio::osecir::Parameters::apply_constraints) + .def_property( + "end_commuter_detection", + [](const mio::osecir::Parameters& self) -> auto { return self.get_end_commuter_detection(); }, + [](mio::osecir::Parameters& self, double p) { + self.get_end_commuter_detection() = p; + }); using Populations = mio::Populations; pymio::bind_Population(m, "Populations", mio::Tag::Populations>{}); @@ -264,6 +270,43 @@ PYBIND11_MODULE(_simulation_osecir, m) }, py::return_value_policy::move); + m.def( + "set_nodes_states", + [](const mio::osecir::Parameters& params, mio::Date start_date, mio::Date end_date, + const std::string& data_dir, const std::string& population_data_path, bool is_node_for_county, + mio::Graph, mio::MobilityParameters>& params_graph, + const std::vector& scaling_factor_inf, double scaling_factor_icu, double tnt_capacity_factor, + int num_days = 0, bool export_time_series = false) { + auto result = mio::set_nodes< + mio::osecir::TestAndTraceCapacity, mio::osecir::ContactPatterns, + mio::osecir::Model, mio::MobilityParameters, mio::osecir::Parameters, + decltype(mio::osecir::read_input_data_state>), decltype(mio::get_node_ids)>( + params, start_date, end_date, data_dir, population_data_path, is_node_for_county, params_graph, + mio::osecir::read_input_data_state>, mio::get_node_ids, scaling_factor_inf, + scaling_factor_icu, tnt_capacity_factor, num_days, export_time_series); + return pymio::check_and_throw(result); + }, + py::return_value_policy::move); + + m.def( + "set_node_germany", + [](const mio::osecir::Parameters& params, mio::Date start_date, mio::Date end_date, + const std::string& data_dir, const std::string& population_data_path, bool is_node_for_county, + mio::Graph, mio::MobilityParameters>& params_graph, + const std::vector& scaling_factor_inf, double scaling_factor_icu, double tnt_capacity_factor, + int num_days = 0, bool export_time_series = false) { + auto result = mio::set_nodes, + mio::osecir::ContactPatterns, mio::osecir::Model, + mio::MobilityParameters, mio::osecir::Parameters, + decltype(mio::osecir::read_input_data_germany>), + decltype(mio::get_country_id)>( + params, start_date, end_date, data_dir, population_data_path, is_node_for_county, params_graph, + mio::osecir::read_input_data_germany>, mio::get_country_id, + scaling_factor_inf, scaling_factor_icu, tnt_capacity_factor, num_days, export_time_series); + return pymio::check_and_throw(result); + }, + py::return_value_policy::move); + pymio::iterable_enum(m, "ContactLocation") .value("Home", ContactLocation::Home) .value("School", ContactLocation::School) @@ -278,12 +321,11 @@ PYBIND11_MODULE(_simulation_osecir, m) auto mobile_comp = {mio::osecir::InfectionState::Susceptible, mio::osecir::InfectionState::Exposed, mio::osecir::InfectionState::InfectedNoSymptoms, mio::osecir::InfectionState::InfectedSymptoms, mio::osecir::InfectionState::Recovered}; - auto weights = std::vector{0., 0., 1.0, 1.0, 0.33, 0., 0.}; + // auto weights = std::vector{1.0}; auto result = mio::set_edges, mio::MobilityParameters, mio::MobilityCoefficientGroup, mio::osecir::InfectionState, - decltype(mio::read_mobility_plain)>(mobility_data_file, params_graph, - mobile_comp, contact_locations_size, - mio::read_mobility_plain, weights); + decltype(mio::read_mobility_plain)>( + mobility_data_file, params_graph, mobile_comp, contact_locations_size, mio::read_mobility_plain, {1.0}); return pymio::check_and_throw(result); }, py::return_value_policy::move); @@ -304,6 +346,16 @@ PYBIND11_MODULE(_simulation_osecir, m) return pymio::check_and_throw(result); }, py::return_value_policy::move); + m.def( + "read_input_data_provincias", + [](std::vector>& model, mio::Date date, const std::vector& provincias, + const std::vector& scaling_factor_inf, double scaling_factor_icu, const std::string& dir, + int num_days = 0, bool export_time_series = false) { + auto result = mio::osecir::read_input_data_provincias>( + model, date, provincias, scaling_factor_inf, scaling_factor_icu, dir, num_days, export_time_series); + return pymio::check_and_throw(result); + }, + py::return_value_policy::move); #endif // MEMILIO_HAS_JSONCPP m.def("interpolate_simulation_result", diff --git a/shellscripts/fitting_graphmodel.sh b/shellscripts/fitting_graphmodel.sh new file mode 100644 index 0000000000..fa53bc085b --- /dev/null +++ b/shellscripts/fitting_graphmodel.sh @@ -0,0 +1,15 @@ +#!/bin/bash +#SBATCH -N 1 +#SBATCH -n 1 +#SBATCH -c 1 +#SBATCH -t 5-0:00:00 +#SBATCH --output=shellscripts/train_countylvl-%A.out +#SBATCH --error=shellscripts/train_countylvl-%A.err +#SBATCH --job-name=train_countylvl +#SBATCH --partition=gpu +#SBATCH --gpus=1 + +module load PrgEnv/gcc13-openmpi-python +module load cuda/12.9.0-none-none-6wnenm2 +source venv/bin/activate +srun --cpu-bind=core python pycode/examples/simulation/graph_germany_nuts3.py \ No newline at end of file