diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 90cd4121b7..91c7a394f5 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -185,6 +185,7 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/ode_seir) add_subdirectory(models/ode_seair) add_subdirectory(models/ode_sir) + add_subdirectory(models/ode_seir_metapop) add_subdirectory(models/sde_sir) add_subdirectory(models/sde_sirs) add_subdirectory(models/sde_seirvv) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index fea58a84cf..381559ea69 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -100,6 +100,10 @@ add_executable(graph_stochastic_mobility_example graph_stochastic_mobility.cpp) target_link_libraries(graph_stochastic_mobility_example PRIVATE memilio ode_secir) target_compile_options(graph_stochastic_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_metapop_seir_example ode_seir_metapop.cpp) +target_link_libraries(ode_metapop_seir_example PRIVATE memilio ode_seir_metapop) +target_compile_options(ode_metapop_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(abm_minimal_example abm_minimal.cpp) target_link_libraries(abm_minimal_example PRIVATE memilio abm) target_compile_options(abm_minimal_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/ode_seir_metapop.cpp b/cpp/examples/ode_seir_metapop.cpp new file mode 100644 index 0000000000..66b9f92c96 --- /dev/null +++ b/cpp/examples/ode_seir_metapop.cpp @@ -0,0 +1,44 @@ +#include "memilio/compartments/simulation.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_seir_metapop/model.h" +#include "models/ode_seir_metapop/parameters.h" +#include "memilio/geography/regions.h" + +int main() +{ + const ScalarType t0 = 0.; + const ScalarType tmax = 10; + ScalarType dt = 0.1; + + using namespace mio::oseirmetapop; + + Model model(3, 1); + + for (size_t i = 0; i < (size_t)model.parameters.get_num_regions(); i++) { + model.populations[{mio::regions::Region(i), mio::AgeGroup(0), InfectionState::Susceptible}] = + 10000; + } + + model.populations[{mio::regions::Region(0), mio::AgeGroup(0), InfectionState::Exposed}] += 100; + model.populations[{mio::regions::Region(0), mio::AgeGroup(0), InfectionState::Susceptible}] -= + 100; + + Eigen::MatrixXd mobility_data_commuter(3, 3); + mobility_data_commuter << 0.4, 0.3, 0.3, 0.2, 0.7, 0.1, 0.4, 0.1, 0.5; + + model.set_commuting_strengths(mobility_data_commuter); + + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() + .setConstant(2.7); + + model.parameters.set>(3.335); + model.parameters.set>(8.097612257); + model.parameters.set>(0.07333); + + auto result = simulate(t0, tmax, dt, model); + + result.print_table(); + return 0; +} diff --git a/cpp/memilio/epidemiology/contact_matrix.h b/cpp/memilio/epidemiology/contact_matrix.h index 83ce44517f..ce0b0704c0 100644 --- a/cpp/memilio/epidemiology/contact_matrix.h +++ b/cpp/memilio/epidemiology/contact_matrix.h @@ -31,7 +31,7 @@ namespace mio { /** - * represents the coefficient wise matrix (or vector) expression B - D * M + * Represents the coefficient wise matrix (or vector) expression B - D * (B - M) * where B is a baseline, M is a minimum and D is some time dependent complex damping factor. * Base class for e.g. time dependent contact matrices. * Coefficient wise expression, so B, D, M matrices must have the same shape. diff --git a/cpp/models/ode_seir_metapop/CMakeLists.txt b/cpp/models/ode_seir_metapop/CMakeLists.txt new file mode 100644 index 0000000000..108b9d7b87 --- /dev/null +++ b/cpp/models/ode_seir_metapop/CMakeLists.txt @@ -0,0 +1,11 @@ +add_library(ode_seir_metapop + model.h + model.cpp + parameters.h +) +target_link_libraries(ode_seir_metapop PUBLIC memilio) +target_include_directories(ode_seir_metapop PUBLIC + $ + $ +) +target_compile_options(ode_seir_metapop PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/ode_seir_metapop/infection_state.h b/cpp/models/ode_seir_metapop/infection_state.h new file mode 100644 index 0000000000..6127c8e868 --- /dev/null +++ b/cpp/models/ode_seir_metapop/infection_state.h @@ -0,0 +1,36 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Carlotta Gerstein +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#ifndef SEIRMETAPOP_INFECTIONSTATE_H +#define SEIRMETAPOP_INFECTIONSTATE_H + +#include "models/ode_seir/infection_state.h" + +namespace mio +{ +namespace oseirmetapop +{ + +using mio::oseir::InfectionState; + +} // namespace oseirmetapop +} // namespace mio + +#endif // SEIRMETAPOP_INFECTIONSTATE_H diff --git a/cpp/models/ode_seir_metapop/model.cpp b/cpp/models/ode_seir_metapop/model.cpp new file mode 100644 index 0000000000..bf0a8cd4eb --- /dev/null +++ b/cpp/models/ode_seir_metapop/model.cpp @@ -0,0 +1,28 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Carlotta Gerstein +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#include "ode_seir_metapop/model.h" + +namespace mio +{ +namespace oseirmetapop +{ + +} // namespace oseirmetapop +} // namespace mio diff --git a/cpp/models/ode_seir_metapop/model.h b/cpp/models/ode_seir_metapop/model.h new file mode 100644 index 0000000000..79a56ac627 --- /dev/null +++ b/cpp/models/ode_seir_metapop/model.h @@ -0,0 +1,289 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Carlotta Gerstein, Martin J. Kuehn, 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 ODESEIRMETAPOP_MODEL_H +#define ODESEIRMETAPOP_MODEL_H + +#include + +#include "memilio/compartments/flow_model.h" +#include "memilio/epidemiology/populations.h" +#include "models/ode_seir_metapop/parameters.h" +#include "models/ode_seir_metapop/infection_state.h" +#include "memilio/geography/regions.h" +#include "memilio/utils/time_series.h" + +namespace mio +{ +namespace oseirmetapop +{ + +/******************** +* define the model * +********************/ + +using Flows = TypeList, + Flow, + Flow>; + +/** + * @brief The Model holds the Parameters and the initial Populations and defines the ordinary differential equations. + */ +template +class Model : public FlowModel, + Parameters, Flows> +{ + + using Base = FlowModel, + Parameters, Flows>; + +public: + using typename Base::ParameterSet; + using typename Base::Populations; + + /** + * @brief Create a Model with the given number of Region%s and AgeGroups. + * @param[in] num_regions The number of Region%s. + * @param[in] num_agegroups The number of AgeGroup%s. + */ + Model(int num_regions, int num_agegroups) + : Base(Populations({Region(num_regions), AgeGroup(num_agegroups), InfectionState::Count}), + ParameterSet(Region(num_regions), AgeGroup(num_agegroups))) + { + } + + /** + * @brief Compute the values of the flows between compartments at a given time. + * @param[in] pop The current population by #InfectionState, AgeGroup and Region. + * @param[in] y The current population by #InfectionState, AgeGroup and Region. + * @param[in] t The current time. + * @param[in, out] flows The computed flows between compartments. + */ + void get_flows(Eigen::Ref> pop, Eigen::Ref> y, FP t, + Eigen::Ref> flows) const override + { + using enum InfectionState; + const auto& params = this->parameters; + const auto& population = this->populations; + const Eigen::MatrixXd commuting_strengths = + params.template get>().get_cont_freq_mat().get_matrix_at(t); + const size_t n_age_groups = (size_t)params.get_num_agegroups(); + const size_t n_regions = (size_t)params.get_num_regions(); + + Eigen::MatrixXd infected_pop(n_regions, n_age_groups); + for (size_t region_n = 0; region_n < n_regions; region_n++) { + for (size_t age_i = 0; age_i < n_age_groups; age_i++) { + infected_pop(region_n, age_i) = + pop[population.get_flat_index({Region(region_n), AgeGroup(age_i), Infected})]; + } + } + Eigen::MatrixXd infectious_share_per_region = commuting_strengths.transpose() * infected_pop; + for (size_t region_n = 0; region_n < n_regions; region_n++) { + for (size_t age_i = 0; age_i < n_age_groups; age_i++) { + infectious_share_per_region(region_n, age_i) /= + params.template get>()[{Region(region_n), AgeGroup(age_i)}]; + } + } + Eigen::MatrixXd infections_due_commuting = commuting_strengths * infectious_share_per_region; + for (size_t age_i = 0; age_i < n_age_groups; age_i++) { + for (size_t age_j = 0; age_j < n_age_groups; age_j++) { + for (size_t region_n = 0; region_n < n_regions; region_n++) { + const size_t Ejn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), Exposed}); + const size_t Ijn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), Infected}); + const size_t Rjn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), Recovered}); + const size_t Sjn = + population.get_flat_index({Region(region_n), AgeGroup(age_j), Susceptible}); + + const double Nj_inv = 1.0 / (pop[Sjn] + pop[Ejn] + pop[Ijn] + pop[Rjn]); + double coeffStoE = + 0.5 * + params.template get>().get_cont_freq_mat().get_matrix_at(t)(age_i, age_j) * + params.template get>()[AgeGroup(age_i)]; + + flows[Base::template get_flat_flow_index( + {Region(region_n), AgeGroup(age_i)})] += + (pop[Ijn] * Nj_inv + infections_due_commuting(region_n, age_j)) * coeffStoE * + y[population.get_flat_index({Region(region_n), AgeGroup(age_i), Susceptible})]; + } + } + for (size_t region = 0; region < n_regions; region++) { + flows[Base::template get_flat_flow_index( + {Region(region), AgeGroup(age_i)})] = + y[population.get_flat_index({Region(region), AgeGroup(age_i), Exposed})] / + params.template get>()[AgeGroup(age_i)]; + flows[Base::template get_flat_flow_index( + {Region(region), AgeGroup(age_i)})] = + y[population.get_flat_index({Region(region), AgeGroup(age_i), Infected})] / + params.template get>()[AgeGroup(age_i)]; + } + } + } + + /** + *@brief Computes the reproduction number at a given index time of the Model output obtained by the Simulation. + *@param[in] t_idx The index time at which the reproduction number is computed. + *@param[in] y The TimeSeries obtained from the Model Simulation. + *@returns The computed reproduction number at the provided index time. + */ + IOResult get_reproduction_number(size_t t_idx, const mio::TimeSeries& y) + { + if (!(t_idx < static_cast(y.get_num_time_points()))) { + return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries"); + } + + auto const& params = this->parameters; + auto const& pop = this->populations; + + const size_t num_age_groups = (size_t)params.get_num_agegroups(); + const size_t num_regions = (size_t)params.get_num_regions(); + constexpr size_t num_infected_compartments = 2; + const size_t total_infected_compartments = num_infected_compartments * num_age_groups * num_regions; + + ContactMatrixGroup const& contact_matrix = params.template get>(); + ContactMatrixGroup const& commuting_strengths = params.template get>(); + Populations const& population_after_commuting = params.template get>(); + + Eigen::MatrixXd F = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + Eigen::MatrixXd V = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + + for (auto i = AgeGroup(0); i < AgeGroup(num_age_groups); i++) { + for (auto n = Region(0); n < Region(num_regions); n++) { + size_t Si = pop.get_flat_index({n, i, InfectionState::Susceptible}); + for (auto j = AgeGroup(0); j < AgeGroup(num_age_groups); j++) { + for (auto m = Region(0); m < Region(num_regions); m++) { + auto const population_region = pop.template slice({m.get(), 1}); + auto const population_region_age = population_region.template slice({j.get(), 1}); + auto Njm = std::accumulate(population_region_age.begin(), population_region_age.end(), 0.); + + double coeffStoE = 0.5 * contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) * + params.template get>()[i]; + if (n == m) { + F(i.get() * num_regions + n.get(), num_age_groups * num_regions + j.get() * num_regions + + m.get()) += coeffStoE * y.get_value(t_idx)[Si] / Njm; + } + for (auto k = Region(0); k < Region(num_regions); k++) { + F(i.get() * num_regions + n.get(), + num_age_groups * num_regions + j.get() * num_regions + m.get()) += + coeffStoE * y.get_value(t_idx)[Si] * + commuting_strengths.get_matrix_at(y.get_time(t_idx))(n.get(), k.get()) * + commuting_strengths.get_matrix_at(y.get_time(t_idx))(m.get(), k.get()) / + population_after_commuting[{k, j}]; + } + } + } + + double T_Ei = params.template get>()[i]; + double T_Ii = params.template get>()[i]; + V(i.get() * num_regions + n.get(), i.get() * num_regions + n.get()) = 1.0 / T_Ei; + V(num_age_groups * num_regions + i.get() * num_regions + n.get(), i.get() * num_regions + n.get()) = + -1.0 / T_Ei; + V(num_age_groups * num_regions + i.get() * num_regions + n.get(), + num_age_groups * num_regions + i.get() * num_regions + n.get()) = 1.0 / T_Ii; + } + } + + V = V.inverse(); + + Eigen::MatrixXd NextGenMatrix = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + NextGenMatrix = F * V; + + //Compute the largest eigenvalue in absolute value + Eigen::ComplexEigenSolver ces; + + ces.compute(NextGenMatrix); + const Eigen::VectorXcd eigen_vals = ces.eigenvalues(); + + Eigen::VectorXd eigen_vals_abs; + eigen_vals_abs.resize(eigen_vals.size()); + + for (int i = 0; i < eigen_vals.size(); i++) { + eigen_vals_abs[i] = std::abs(eigen_vals[i]); + } + return mio::success(eigen_vals_abs.maxCoeff()); + } + + /** + *@brief Computes the reproduction number for all time points of the Model output obtained by the Simulation. + *@param[in] y The TimeSeries obtained from the Model Simulation. + *@returns Vector containing all reproduction numbers + */ + Eigen::VectorXd get_reproduction_numbers(const mio::TimeSeries& y) + { + auto num_time_points = y.get_num_time_points(); + Eigen::VectorXd temp(num_time_points); + for (size_t i = 0; i < static_cast(num_time_points); i++) { + temp[i] = get_reproduction_number(i, y).value(); + } + return temp; + } + + /** + * @brief Set the CommutingStrengths matrix and update the PopulationAfterCommuting. + * @param[in] commuting_strengths The matrix containing the fractions of the commuting population between Region%s. + */ + void set_commuting_strengths(const Eigen::MatrixXd& commuting_strengths) + { + auto& commuting_strengths_param = + this->parameters.template get>().get_cont_freq_mat()[0].get_baseline(); + commuting_strengths_param = commuting_strengths; + + auto number_regions = (size_t)this->parameters.get_num_regions(); + auto number_age_groups = (size_t)this->parameters.get_num_agegroups(); + auto& population = this->populations; + mio::Populations population_after_commuting( + {Region(number_regions), AgeGroup(number_age_groups)}, 0.); + + for (size_t region_n = 0; region_n < number_regions; ++region_n) { + for (size_t age = 0; age < number_age_groups; ++age) { + double population_n = 0; + for (size_t state = 0; state < (size_t)InfectionState::Count; state++) { + population_n += population[{Region(region_n), mio::AgeGroup(age), + InfectionState(state)}]; + } + population_after_commuting[{Region(region_n), mio::AgeGroup(age)}] += population_n; + for (size_t region_m = 0; region_m < number_regions; ++region_m) { + population_after_commuting[{Region(region_n), mio::AgeGroup(age)}] -= + commuting_strengths(region_n, region_m) * population_n; + population_after_commuting[{Region(region_m), mio::AgeGroup(age)}] += + commuting_strengths(region_n, region_m) * population_n; + } + } + } + this->parameters.template get>() = population_after_commuting; + } + + /** + * @brief Set the CommutingStrengths matrix without data. + * This function sets the CommutingStrengths matrix to the identity matrix, which corresponds to no commuting between Region%s. + * This prevents division by zero but does not produce meaningful results. + */ + void set_commuting_strengths() + { + auto number_regions = (size_t)this->parameters.get_num_regions(); + set_commuting_strengths(Eigen::MatrixXd::Identity(number_regions, number_regions)); + } +}; + +} // namespace oseirmetapop +} // namespace mio + +#endif // ODESEIRMOBILITY_MODEL_H diff --git a/cpp/models/ode_seir_metapop/parameters.h b/cpp/models/ode_seir_metapop/parameters.h new file mode 100644 index 0000000000..e293186888 --- /dev/null +++ b/cpp/models/ode_seir_metapop/parameters.h @@ -0,0 +1,324 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Carlotta Gerstein +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#ifndef ODESEIRMETAPOP_PARAMETERS_H +#define ODESEIRMETAPOP_PARAMETERS_H + +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/geography/regions.h" + +#include + +namespace mio +{ +namespace oseirmetapop +{ + +using Region = mio::regions::Region; + +/**************************************************** +* Define Parameters of the SEIR model with mobility * +****************************************************/ + +/** + * @brief Probability of getting infected from a contact. + */ +template +struct TransmissionProbabilityOnContact { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 1.0); + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/** + * @brief The latent time in day unit. + */ +template +struct TimeExposed { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 5.2); + } + static std::string name() + { + return "TimeExposed"; + } +}; + +/** + * @brief The infectious time in day unit. + */ +template +struct TimeInfected { + using Type = CustomIndexArray, AgeGroup>; + static Type get_default(Region, AgeGroup size) + { + return Type(size, 6.0); + } + static std::string name() + { + return "TimeInfected"; + } +}; + +/** + * @brief The contact patterns within the society are modelled using a + * ContactMatrix. + */ +template +struct ContactPatterns { + using Type = UncertainContactMatrix; + static Type get_default(Region, AgeGroup size) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/** + * @brief The commuting patterns between different Region%s are modelled using a ContactMatrix of size n_regions x n_regions. + * Each entry of the matrix represents the fraction of individuals commuting from one Region to another. The diagonal corresponds + * to the fraction of individuals staying in their Region. + */ +template +struct CommutingStrengths { + using Type = UncertainContactMatrix; + static Type get_default(Region size, AgeGroup) + { + return Type(1, static_cast((size_t)size)); + } + static std::string name() + { + return "CommutingStrengths"; + } +}; + +/** + * @brief The number of individuals in each Region and AgeGroup if commuting was applied. + * Computed as the sum of the number of individuals staying in their Region and the number of individuals commuting to this Region + * minus the number of individuals commuting from this Region. + */ +template +struct PopulationAfterCommuting { + using Type = Populations; + static Type get_default(Region size_regions, AgeGroup size_agegroups) + { + return Type({size_regions, size_agegroups}, 0.); + } + static std::string name() + { + return "PopulationAfterCommuting"; + } +}; + +template +using ParametersBase = ParameterSet, TimeExposed, TimeInfected, + ContactPatterns, CommutingStrengths, PopulationAfterCommuting>; + +/** + * @brief Parameters of the SEIR metapopulation model. + */ +template +class Parameters : public ParametersBase +{ +public: + Parameters(Region num_regions, AgeGroup num_agegroups) + : ParametersBase(num_regions, num_agegroups) + , m_num_regions{num_regions} + , m_num_agegroups(num_agegroups) + { + } + + Region get_num_regions() const + { + return m_num_regions; + } + + AgeGroup get_num_agegroups() const + { + return m_num_agegroups; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and applies them, if they do not. + * + * Attention: This function should be used with care. It is necessary for some test problems to run through quickly, + * but in a manual execution of an example, check_constraints() may be preferred. Note that the apply_constraints() + * function can and will not set Parameters to meaningful values in an epidemiological or virological context, + * as all models are designed to be transferable to multiple diseases. Consequently, only acceptable + * (like 0 or 1 for probabilities or small positive values for time spans) values are set here and a manual adaptation + * may often be necessary to have set meaningful values. + * + * @return Returns true if one ore more constraints were corrected, false otherwise. + */ + bool apply_constraints() + { + double tol_times = 1e-1; + int corrected = false; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + if (this->template get>()[i] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], tol_times); + this->template get>()[i] = tol_times; + corrected = true; + } + if (this->template get>()[i] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], tol_times); + this->template get>()[i] = tol_times; + corrected = true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_warning( + "Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", + this->template get>()[i], 0.0); + this->template get>() = 0.0; + corrected = true; + } + for (auto j = Region(0); j < Region(m_num_regions); j++) { + if (this->template get>()[{j, i}] <= 0.0) { + log_warning( + "Constraint check: Parameter PopulationAfterCommuting changed from {:.4f} to {:.4f}. Please " + "note that this only prevents division by zero. Consider to cancel and reset parameters.", + this->template get>()[{j, i}], 1.0); + this->template get>()[{j, i}] = 1.0; + corrected = true; + } + } + } + if ((this->template get>().get_cont_freq_mat().get_matrix_at(0).rowwise().sum() - + Eigen::VectorXd::Ones((size_t)this->get_num_regions())) + .cwiseAbs() + .maxCoeff() > 1e-10 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).minCoeff() < 0.0 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).maxCoeff() > 1.0) { + log_warning("Constraint check: Parameter CommutingStrengths does not ensure that the number of people " + "staying equals the complement of those leaving. Running without commuting."); + this->template get>().get_cont_freq_mat()[0].get_baseline() = + Eigen::MatrixXd::Identity((size_t)this->get_num_regions(), (size_t)this->get_num_regions()); + corrected = true; + } + return corrected; + } + + /** + * @brief Checks whether all Parameters satisfy their corresponding constraints and logs an error + * if constraints are not satisfied. + * @return Returns true if one constraint is not satisfied, otherwise false. + */ + bool check_constraints() const + { + const double tol_times = 1e-1; + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + if (this->template get>()[i] < tol_times) { + log_error( + "Constraint check: Parameter TimeExposed {:.4f} smaller or equal {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], 0.0); + return true; + } + if (this->template get>()[i] < tol_times) { + log_error( + "Constraint check: Parameter TimeInfected {:.4f} smaller or equal {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[i], 0.0); + return true; + } + if (this->template get>()[i] < 0.0 || + this->template get>()[i] > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact {:.4f} smaller {:.4f} or " + "greater {:.4f}", + this->template get>()[i], 0.0, 1.0); + return true; + } + for (auto j = Region(0); j < Region(m_num_regions); j++) { + if (this->template get>()[{j, i}] <= 0.0) { + log_error("Constraint check: Parameter PopulationAfterCommuting {:.4f} smaller or equal {:.4f}", + this->template get>()[{j, i}], 0.0); + return true; + } + } + } + if ((this->template get>().get_cont_freq_mat().get_matrix_at(0).rowwise().sum() - + Eigen::VectorXd::Ones((size_t)this->get_num_regions())) + .cwiseAbs() + .maxCoeff() > 1e-10 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).minCoeff() < 0.0 || + this->template get>().get_cont_freq_mat().get_matrix_at(0).maxCoeff() > 1.0) { + log_error("Constraint check: Parameter CommutingStrengths does not ensure that the number of people " + "staying equals the complement of those leaving."); + return true; + } + + return false; + } + +private: + Parameters(ParametersBase&& base) + : ParametersBase(std::move(base)) + , m_num_regions(base.get_num_regions()) + , m_num_agegroups(base.get_num_agegroups()) + { + } + +public: + /** + * @brief Deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); + return success(Parameters(std::move(base))); + } + +private: + Region m_num_regions; + AgeGroup m_num_agegroups; +}; + +} // namespace oseirmetapop +} // namespace mio + +#endif // SEIR_PARAMETERS_H diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 936dc46b89..cdf08e2bdf 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -86,6 +86,7 @@ set(TESTSOURCES temp_file_register.h test_graph_abm.cpp test_abstract_parameter_dist.cpp + test_odemetapop.cpp test_temporal_hybrid_model.cpp ) diff --git a/cpp/tests/data/ode-seir-compare-euler.csv b/cpp/tests/data/ode-seir-compare-euler.csv new file mode 100644 index 0000000000..66722c7ad1 --- /dev/null +++ b/cpp/tests/data/ode-seir-compare-euler.csv @@ -0,0 +1,502 @@ + # t S E I R +0.00000000000000 1049000.00000000000000 10000.00000000000000 1000.00000000000000 1000.00000000000000 +0.10000000000000 1048973.30537229031324 9834.38693540201530 1175.64102564102564 1016.66666666666663 +0.20000000000000 1048941.92287142318673 9676.64661058845013 1345.16983422781641 1036.26068376068383 +0.30000000000000 1048906.01605155738071 9526.46407255830127 1508.83969488661796 1058.68018099781420 +0.40000000000000 1048865.74174628220499 9383.53714566897543 1666.89359880309030 1083.82750924592460 +0.50000000000000 1048821.25033727660775 9247.57591725778730 1819.56467623974982 1111.60906922597610 +0.60000000000000 1048772.68601242476143 9118.30224370071664 1967.07659671122428 1141.93514716330537 +0.70000000000000 1048720.18701378209516 8995.44927611847015 2109.64395299105081 1174.71975710849256 +0.80000000000000 1048663.88587577571161 8878.76100496862091 2247.47262959732416 1209.88048965834332 +0.90000000000000 1048603.90965401218273 8767.99182279047091 2380.76015637907312 1247.33836681829871 +1.00000000000000 1048540.38014504511375 8662.90610439617012 2509.69604780077725 1287.01770275795002 +1.10000000000000 1048473.41409745207056 8563.27780382778474 2634.46212849889571 1328.84597022129628 +1.20000000000000 1048403.12341455032583 8468.89006742519268 2755.23284566162783 1372.75367236294460 +1.30000000000000 1048329.61534907401074 8379.53486237413017 2872.17556876134176 1418.67421979063829 +1.40000000000000 1048252.99269012128934 8295.01262012730876 2985.45087714815509 1466.54381260332730 +1.50000000000000 1048173.35394267027732 8215.13189411430903 3095.21283599300568 1516.30132722246321 +1.60000000000000 1048090.79349995055236 8139.70903117802754 3201.60926104916689 1567.88820782234666 +1.70000000000000 1048005.40180894767400 8068.56785619668608 3304.78197268254007 1621.24836217316602 +1.80000000000000 1047917.26552930800244 8001.53936937099570 3404.86703960315208 1676.32806171787502 +1.90000000000000 1047826.46768589981366 7938.46145567586791 3501.99501271305462 1733.07584571126085 +2.00000000000000 1047733.08781527902465 7879.17860599520827 3596.29114946929622 1791.44242925647836 +2.10000000000000 1047637.20210629759822 7823.54164947676873 3687.87562914471573 1851.38061508096666 +2.20000000000000 1047538.88353508408181 7771.40749666184274 3776.86375935403657 1912.84520890004524 +2.30000000000000 1047438.20199461758602 7722.63889296173602 3863.36617419804315 1975.79293822261252 +2.40000000000000 1047335.22441910672933 7677.10418206948179 3947.48902436451954 2040.18237445924660 +2.50000000000000 1047230.01490337902214 7634.67707891121063 4029.33415951106235 2105.97385819865531 +2.60000000000000 1047122.63481747731566 7595.23645175693855 4108.99930324186244 2173.12942752383969 +2.70000000000000 1047013.14291665202472 7558.66611312532677 4186.57822097802909 2241.61274924453755 +2.80000000000000 1046901.59544693224598 7524.85461913120889 4262.16088100901015 2311.38905292750451 +2.90000000000000 1046788.04624644946307 7493.69507693841297 4335.83360900112712 2382.42506761098821 +3.00000000000000 1046672.54684268392157 7465.08495999359275 4407.67923622812941 2454.68896109434036 +3.10000000000000 1046555.14654579432681 7438.92593072951149 4477.77724177804976 2528.15028169814241 +3.20000000000000 1046435.89253818674479 7415.12367043845279 4546.20388898039346 2602.77990239444307 +3.30000000000000 1046314.82996047323104 7393.58771602820434 4613.03235628787024 2678.54996721078305 +3.40000000000000 1046192.00199396267999 7374.23130338440114 4678.33286283746111 2755.43383981558100 +3.50000000000000 1046067.44993982347660 7356.97121707389397 4742.17278890653688 2833.40605419620533 +3.60000000000000 1045941.21329505008180 7341.72764613432537 4804.61679147105406 2912.44226734464746 +3.70000000000000 1045813.32982536125928 7328.42404570514464 4865.72691506450428 2992.51921386916501 +3.80000000000000 1045683.83563515322749 7316.98700426501910 4925.56269812827213 3073.61466245357360 +3.90000000000000 1045552.76523462473415 7307.34611624991157 4984.18127503635878 3155.70737408904461 +4.00000000000000 1045420.15160418860614 7299.43385983506141 5041.63747397004590 3238.77706200631746 +4.10000000000000 1045286.02625627722591 7293.18547967273389 5097.98391081096270 3322.80435323915162 +4.20000000000000 1045150.41929464729037 7288.53887438587572 5153.27107921422976 3407.77075175266782 +4.30000000000000 1045013.35947128455155 7285.43448862579680 5207.54743701679763 3493.65860307290495 +4.40000000000000 1044874.87424100411590 7283.81520950963841 5260.85948912983440 3580.45106035651816 +4.50000000000000 1044734.98981383931823 7283.62626726075268 5313.25186705798387 3668.13205184201524 +4.60000000000000 1044593.73120530904271 7284.81513988219194 5364.76740518254519 3756.68624962631475 +4.70000000000000 1044451.12228464707732 7287.33146170028613 5415.44721394005774 3846.09903971269068 +4.80000000000000 1044307.18582107743714 7291.12693562182631 5465.33075002247278 3936.35649327835836 +4.90000000000000 1044161.94352821342181 7296.15524895464205 5514.45588371995382 4027.44533911206645 +5.00000000000000 1044015.41610665619373 7302.37199264736955 5562.85896352246709 4119.35293717406603 +5.10000000000000 1043867.62328486575279 7309.73458381001910 5610.57487809159375 4212.06725323277351 +5.20000000000000 1043718.58385837392416 7318.20219138248194 5657.63711570949090 4305.57683453430036 +5.30000000000000 1043568.31572740629781 7327.73566482348087 5704.07782130758551 4399.87078646279224 +5.40000000000000 1043416.83593297796324 7338.29746569757208 5749.92785117342373 4494.93875015125195 +5.50000000000000 1043264.16069152322598 7349.85160204274962 5795.21682543010229 4590.77088100414221 +5.60000000000000 1043110.30542812077329 7362.36356540592442 5839.97317837888386 4687.35782809464399 +5.70000000000000 1042955.28480836923700 7375.80027043808968 5884.22420679191418 4784.69071440095831 +5.80000000000000 1042799.11276896891650 7390.12999694534756 5927.99611623842247 4882.76111784748991 +5.90000000000000 1042641.80254706146661 7405.32233429615735 5971.31406552442331 4981.56105311813008 +5.99999999999999 1042483.36670837807469 7421.34812808919378 6014.20220932266056 5081.08295421020375 +6.09999999999999 1042323.81717424478848 7438.17942899004629 6056.68373906643410 5181.31965769891485 +6.19999999999999 1042163.16524749016389 7455.78944364871586 6098.78092217795620 5282.26438668335504 +6.29999999999999 1042001.42163730144966 7474.15248761340627 6140.51513969900407 5383.91073538632099 +6.39999999999999 1041838.59648307040334 7493.24394015953840 6181.90692238889369 5486.25265438130464 +6.49999999999999 1041674.69937727064826 7513.04020095617489 6222.97598535214638 5589.28443642111961 +6.59999999999999 1041509.73938740557060 7533.51864849521189 6263.74126125569092 5693.00070284365574 +6.69999999999999 1041343.72507706412580 7554.65760021169444 6304.22093219300405 5797.39639053125029 +6.79999999999999 1041176.66452612215653 7576.43627422652935 6344.43246025026838 5902.46673940113396 +6.89999999999999 1041008.56535012181848 7598.83475264564095 6384.39261682737651 6008.20728040530503 +6.99999999999999 1040839.43471886427142 7621.83394635228979 6424.11751076446490 6114.61382401909486 +7.09999999999999 1040669.27937424718402 7645.41556123183727 6463.62261532260163 6221.68244919850258 +7.19999999999999 1040498.10564737697132 7669.56206577071134 6502.92279406527268 6329.40949278721291 +7.29999999999999 1040325.91947498614900 7694.25665997366468 6542.03232568541625 6437.79153935496743 +7.39999999999999 1040152.72641518386081 7719.48324554571263 6580.96492782092173 6546.82541144972402 +7.49999999999999 1039978.53166256635450 7745.22639728729064 6619.73377989978508 6656.50816024673895 +7.59999999999999 1039803.34006271406543 7771.47133565327385 6658.35154505441642 6766.83705657840164 +7.69999999999999 1039627.15612609940581 7798.20390042849249 6696.83039114299572 6877.80958232930880 +7.79999999999999 1039449.98404142993968 7825.41052547430809 6735.18201091423725 6989.42342218169233 +7.89999999999999 1039271.82768845011014 7853.07821450265783 6773.41764135042922 7101.67645569692922 +7.99999999999999 1039092.69065022270661 7881.19451783574004 6811.54808222220436 7214.56674971943630 +8.09999999999999 1038912.57622491195798 7909.74751011121953 6849.58371388713658 7328.09255108980597 +8.19999999999999 1038731.48743708815891 7938.72576889445736 6887.53451436295109 7442.25227965459180 +8.29999999999999 1038549.42704857385252 7968.11835416083159 6925.41007570487182 7557.04452156064144 +8.39999999999999 1038366.39756884961389 7997.91478861272299 6963.21961971544806 7672.46802282238968 +8.49999999999999 1038182.40126503794454 8028.10503879717271 7000.97201301402492 7788.52168315098061 +8.59999999999999 1037997.44017148204148 8058.67949699160363 7038.67578149194196 7905.20455003454754 +8.69999999999999 1037811.51609893597197 8089.62896382632698 7076.33912417845295 8022.51581305941363 +8.79999999999998 1037624.63064338255208 8120.94463161381373 7113.96992654136920 8140.45479846238777 +8.89999999999998 1037436.78519449394662 8152.61806835595416 7151.57577324543217 8259.02096390474435 +8.99999999999998 1037247.98094374907669 8184.64120240166994 7189.16396039049505 8378.21389345883472 +9.09999999999998 1037058.21889222227037 8217.00630772839759 7226.74150725068557 8498.03329279867648 +9.19999999999998 1036867.49985805689357 8249.70598982201773 7264.31516753487449 8618.47898458618874 +9.29999999999998 1036675.82448363606818 8282.73317213084374 7301.89144018792194 8739.55090404510338 +9.39999999999998 1036483.19324246328324 8316.08108307028670 7339.47657975140828 8861.24909471490173 +9.49999999999998 1036289.60644576488994 8349.74324355574026 7377.07660630177452 8983.57370437742429 +9.59999999999998 1036095.06424882542342 8383.71345504218698 7414.69731498307374 9106.52498114912123 +9.69999999999998 1035899.56665706692729 8417.98578804984027 7452.34428515083346 9230.10326973217161 +9.79999999999998 1035703.11353188287467 8452.55457115603895 7490.02288914286783 9354.30900781801938 +9.89999999999998 1035505.70459623611532 8487.41438043438211 7527.73830069220548 9479.14272263706698 +9.99999999999998 1035307.33944003155921 8522.56002932286719 7565.49550299671409 9604.60502764860394 +10.09999999999998 1035108.01752527162898 8557.98655890355622 7603.29929645938773 9730.69661936521516 +10.19999999999998 1034907.73819100391120 8593.68922857697908 7641.15430611269676 9857.41827430620469 +10.29999999999998 1034706.50065806997009 8629.66350711520136 7679.06498873986311 9984.77084607474899 +10.39999999999998 1034504.30403366254177 8665.90506407809880 7717.03563970538835 10112.75526255374643 +10.49999999999998 1034301.14731569960713 8702.40976157803925 7755.07039950667240 10241.37252321550295 +10.59999999999998 1034097.02939702232834 8739.17364637876926 7793.17326005806171 10370.62369654061331 +10.69999999999998 1033891.94906942511443 8776.19294231487220 7831.34807071822433 10500.50991754158167 +10.79999999999998 1033685.90502752293833 8813.46404301872462 7869.59854407128387 10631.03238538688493 +10.89999999999998 1033478.89587246428709 8850.98350494242368 7907.92826147173764 10762.19236112140607 +10.99999999999998 1033270.92011549521703 8888.74804066263641 7946.34067836276790 10893.99116547926860 +11.09999999999998 1033061.97618138056714 8926.75451245685144 7984.83912937715741 11026.43017678531396 +11.19999999999998 1032852.06241168873385 8964.99992613994982 8023.42683322965695 11159.51082894160027 +11.29999999999998 1032641.17706794477999 9003.48142515049040 8062.10689740928956 11293.23460949542823 +11.39999999999998 1032429.31833465816453 9042.19628487650698 8100.88232267972126 11427.60305778558359 +11.49999999999998 1032216.48432222986594 9081.14190721106388 8139.75600739550464 11562.61776316357827 +11.59999999999997 1032002.67306974332314 9120.31581532819291 8178.73075164168949 11698.28036328683629 +11.69999999999997 1031791.03632657043636 9156.56186974472439 8217.80926120397271 11834.59254248086472 +11.79999999999997 1031587.62644657120109 9183.88402147963097 8256.93350178156106 11971.55603016759778 +11.89999999999997 1031397.67450624494813 9197.22280754671010 8295.93109767775786 12109.17158853062392 +11.99999999999997 1031225.07032349496149 9192.95732092076469 8334.53524875902622 12247.43710682525307 +12.09999999999997 1031071.96818014676683 9169.27182348207680 8372.41396873331360 12386.34602763790281 +12.19999999999997 1030938.55839373019990 9126.34945944704668 8409.20588637267065 12525.88626045012461 +12.29999999999997 1030823.02987110440154 9066.37126169881230 8444.55917530710758 12666.03969188966948 +12.39999999999997 1030721.73071303451434 8993.31712627448724 8478.16981587952978 12806.78234481145410 +12.49999999999997 1030629.51430479332339 8912.58512824118043 8509.81539188938041 12948.08517507611214 +12.59999999999997 1030540.23905261140317 8830.46451257228910 8539.38100320868216 13089.91543160760193 +12.69999999999997 1030450.66139191738330 8750.22554802458217 8566.87461173031261 13232.23844832774739 +12.79999999999997 1030360.80313613510225 8671.81023557558910 8592.36693643271610 13375.01969185658527 +12.89999999999997 1030270.68534965766594 8595.16244059968813 8615.92640227888114 13518.22580746379754 +12.99999999999997 1030180.32837217091583 8520.22783269034517 8637.61921430371513 13661.82458083511301 +13.09999999999997 1030089.75184219516814 8446.95382742206675 8657.50942930936435 13805.78490107350808 +13.19999999999997 1029998.97471986932214 8375.28952998984460 8675.65902524565900 13950.07672489533070 +13.29999999999997 1029908.01530900120270 8305.18568066581793 8692.12796835034351 14094.67104198275774 +13.39999999999997 1029816.89127840858418 8236.59460201480942 8706.97427812141177 14239.53984145526374 +13.49999999999997 1029725.61968257243279 8169.47014781224516 8720.25409019146718 14384.65607942395400 +13.59999999999997 1029634.21698162471876 8103.76765360976515 8732.02171717184501 14529.99364759381206 +13.69999999999997 1029542.69906069256831 8039.44388889557194 8742.32970753198970 14675.52734288001011 +13.79999999999997 1029451.08124861877877 7976.45701079825085 8751.22890257752442 14821.23283800554418 +13.89999999999997 1029359.37833607883658 7914.76651928443425 8758.76849158837831 14967.08665304850365 +13.99999999999997 1029267.60459311318118 7854.33321380225789 8764.99606517635038 15113.06612790831059 +14.09999999999997 1029175.77378609508742 7795.11915132409922 8769.95766691960853 15259.14939566125031 +14.19999999999997 1029083.89919415011536 7737.08760574356620 8773.69784332974450 15405.31535677657666 +14.29999999999997 1028991.99362504621968 7680.20302858315245 8776.25969220521438 15551.54365416540531 +14.39999999999996 1028900.06943057093304 7624.43101097035560 8777.68490942326571 15697.81464903549204 +14.49999999999996 1028808.13852141203824 7569.73824684141709 8778.01383422077015 15844.10939752587910 +14.59999999999996 1028716.21238155802712 7516.09249733314300 8777.28549301275962 15990.40962809622579 +14.69999999999996 1028624.30208223336376 7463.46255632453358 8775.53764179587597 16136.69771964643769 +14.79999999999996 1028532.41829538356978 7411.81821709116957 8772.80680718244184 16282.95668034303526 +14.89999999999996 1028440.57130672491621 7361.13024003650116 8769.12832610935948 16429.17012712940777 +14.99999999999996 1028348.77102837234270 7311.37032146532511 8764.53638426464931 16575.32226589789570 +15.09999999999996 1028257.02701105899177 7262.51106336585144 8759.06405327303401 16721.39787230230650 +15.19999999999996 1028165.34845596156083 7214.52594416784177 8752.74332668064744 16867.38227319018915 +15.29999999999996 1028073.74422614241485 7167.38929044534143 8745.60515477765875 17013.26132863486782 +15.39999999999996 1027982.22285762254614 7121.07624953354389 8737.67947829633886 17159.02141454782759 +15.49999999999996 1027890.79257009632420 7075.56276303029063 8728.99526102089112 17304.64940585276781 +15.59999999999996 1027799.46127729932778 7030.82554115367839 8719.58052134420177 17450.13266020311494 +15.69999999999996 1027708.23659704113379 6986.84203792814151 8709.46236280552330 17595.45900222551791 +15.79999999999996 1027617.12586091377307 6943.59042717227931 8698.66700364199824 17740.61670827227499 +15.89999999999996 1027526.13612368551549 6901.04957926255156 8687.21980538589378 17885.59449166630657 +15.99999999999996 1027435.27417239139322 6859.19903864779553 8675.14530053835733 18030.38148842273949 +16.09999999999996 1027344.54653512872756 6818.01900209033647 8662.46721934953530 18174.96724343171081 +16.19999999999996 1027253.95948956860229 6777.49029761022393 8649.20851573390792 18319.34169708753689 +16.29999999999996 1027163.51907119178213 6737.59436410989656 8635.39139234879622 18463.49517234976884 +16.39999999999996 1027073.23108125769068 6698.31323165730373 8621.03732486304580 18607.41836222224811 +16.49999999999996 1026983.10109451541211 6659.62950240622013 8606.16708544207177 18751.10231763663251 +16.59999999999997 1026893.13446666521486 6621.52633213316858 8590.80076547456702 18894.53843572733240 +16.69999999999997 1026803.33634157862980 6583.98741237104878 8574.95779756537195 19037.71844848524051 +16.79999999999997 1026713.71165828395169 6546.99695312017957 8558.65697681821257 19180.63441177799541 +16.89999999999997 1026624.26515772601124 6510.53966611811575 8541.91648143124621 19323.27869472496604 +16.99999999999997 1026535.00138930708636 6474.60074865018305 8524.75389262761200 19465.64396941548694 +17.09999999999997 1026445.92471721535549 6439.16586788325912 8507.18621394247566 19607.72320095928080 +17.19999999999997 1026357.03932654880919 6404.22114570589747 8489.22988988734323 19749.50963785832209 +17.29999999999998 1026268.34922924044076 6369.75314405843164 8470.90082401177096 19890.99680268977681 +17.39999999999998 1026179.85826979123522 6335.74885073722544 8452.21439638192896 20032.17848308997418 +17.49999999999998 1026091.57013081805781 6302.19566565775040 8433.18548049486890 20173.04872302967124 +17.59999999999998 1026003.48833842074964 6269.08138756165772 8413.82845964670560 20313.60181437125357 +17.69999999999998 1025915.61626737576444 6236.39420115349640 8394.15724277236950 20453.83228869869708 +17.79999999999998 1025827.95714616158511 6204.12266465319226 8374.18527977398662 20593.73490941157070 +17.89999999999998 1025740.51406182057690 6172.25569775085205 8353.92557635441699 20733.30466407446875 +17.99999999999999 1025653.28996466379613 6140.78256995088032 8333.39070837192412 20872.53675701370958 +18.09999999999999 1025566.28767282282934 6109.69288929283175 8312.59283573144785 21011.42660215324213 +18.19999999999999 1025479.50987665401772 6078.97659143681904 8291.54371582745262 21149.96981608209899 +18.29999999999999 1025392.95914299995638 6048.62392910169092 8270.25471655283218 21288.16221134589068 +18.39999999999999 1025306.63791931280866 6018.62546184457460 8248.73682888788244 21425.99978995510537 +18.49999999999999 1025220.54853764350992 5988.97204617074931 8227.00067908291567 21563.47873710323620 +18.59999999999999 1025134.69321850163396 5959.65482596317088 8205.05654044763833 21700.59541508795155 +18.70000000000000 1025049.07407459034584 5930.66522322131277 8182.91434475998267 21837.34635742874525 +18.80000000000000 1024963.69311441958416 5901.99492909932178 8160.58369330670030 21973.72826317474392 +18.90000000000000 1024878.55224580236245 5873.63589523381324 8138.07386756760116 22109.73799139652328 +19.00000000000000 1024793.65327923744917 5845.58032535193888 8115.39383955494486 22245.37255585598177 +19.10000000000000 1024708.99793118191883 5817.82066715066867 8092.55228181913026 22380.62911984856328 +19.20000000000000 1024624.58782721811440 5790.34960443851378 8069.55757713145249 22515.50499121221583 +19.30000000000000 1024540.42450511700008 5763.16004953120864 8046.41782785436135 22649.99761749773825 +19.40000000000001 1024456.50941780256107 5736.24513589313483 8023.14086500931171 22784.10458129531253 +19.50000000000001 1024372.84393621969502 5709.59821101655052 7999.73425705197315 22917.82359571213601 +19.60000000000001 1024289.42935210885480 5683.21282953092395 7976.20531836424561 23051.15249999633670 +19.70000000000001 1024206.26688069081865 5657.08274653494482 7952.56111747223076 23184.08925530240595 +19.80000000000001 1024123.35766326379962 5631.20191114399677 7928.80848499900640 23316.63194059360831 +19.90000000000001 1024040.70276971661951 5605.56446024614161 7904.95402136076609 23448.77874867692663 +20.00000000000001 1023958.30320095969364 5580.16471245985758 7881.00410421461493 23580.52798236627132 +20.10000000000002 1023876.15989127755165 5554.99716228701982 7856.96489566603486 23711.87805076984660 +20.20000000000002 1023794.27371060429141 5530.05647445480008 7832.84234924378688 23842.82746569761366 +20.30000000000002 1023712.64546672534198 5505.33747844038680 7808.64221664975230 23973.37483818500914 +20.40000000000002 1023631.27590740774758 5480.83516317260455 7784.37005429098190 24103.51887512917165 +20.50000000000002 1023550.16572246071883 5456.54467190472133 7760.03122960098972 24233.25837603402033 +20.60000000000002 1023469.31554572971072 5432.46129725290120 7735.63092715708990 24362.59222986070381 +20.70000000000002 1023388.72595702507533 5408.58047639494816 7711.17415460036045 24491.51941197998894 +20.80000000000003 1023308.39748398831580 5384.89778642416695 7686.66574836461587 24620.03898122332976 +20.90000000000003 1023228.33060389722232 5361.40893985330786 7662.11037922054220 24748.15007702940784 +21.00000000000003 1023148.52574541268405 5338.10978026376051 7637.51255764096823 24875.85191668308471 +21.10000000000003 1023068.98329026834108 5314.99627809528829 7612.87663899304971 25003.14379264376839 +21.20000000000003 1022989.70357490540482 5292.06452657176578 7588.20682856294661 25130.02506996032025 +21.30000000000003 1022910.68689205381088 5269.31073775851746 7563.50718641840558 25256.49518376970445 +21.40000000000003 1022831.93349226226564 5246.73123874700559 7538.78163211448100 25382.55363687667705 +21.50000000000004 1022753.44358537835069 5224.32246796275012 7514.03394924745135 25508.19999741191714 +21.60000000000004 1022675.21734198008198 5202.08097159249337 7489.26778986184127 25633.43389656604268 +21.70000000000004 1022597.25489476136863 5180.00340012675588 7464.48667871528141 25758.25502639707338 +21.80000000000004 1022519.55633987160400 5158.08650501405464 7439.69401740579724 25882.66313770899433 +21.90000000000004 1022442.12173821218312 5136.32713542316924 7414.89308836597047 26006.65803799909190 +22.00000000000004 1022364.95111669029575 5114.72223510996901 7390.08705872826522 26130.23958947185747 +22.10000000000004 1022288.04446943197399 5093.26883938541141 7365.27898406567783 26253.40770711732694 +22.20000000000005 1022211.40175895544235 5071.96407218145123 7340.47181201173862 26376.16235685175343 +22.30000000000005 1022135.02291730628349 5050.80514321168994 7315.66838576374994 26498.50355371861588 +22.40000000000005 1022058.90784715558402 5029.78934522370218 7290.87144747304046 26620.43136014801348 +22.50000000000005 1021983.05642286187503 5008.91405134007437 7266.08364152586910 26741.94588427256531 +22.60000000000005 1021907.46849149861373 4988.17671248529496 7241.30751771851646 26863.04727829799594 +22.70000000000005 1021832.14387384813745 4967.57485489571354 7216.54553432997636 26983.73573692663922 +22.80000000000005 1021757.08236536290497 4947.10607770988463 7191.80006109554779 27104.01149583213919 +22.90000000000006 1021682.28373709553853 4926.76805063670326 7167.07338208452984 27223.87483018373314 +23.00000000000006 1021607.74773659813218 4906.55851169880953 7142.36769848510903 27343.32605321847586 +23.10000000000006 1021533.47408879233990 4886.47526504883626 7117.68513129943676 27462.36551485989548 +23.20000000000006 1021459.46249681105837 4866.51617885614087 7093.02772395179545 27580.99360038155282 +23.30000000000006 1021385.71264281205367 4846.67918326174913 7068.39744481265279 27699.21072911408191 +23.40000000000006 1021312.22418876562733 4826.96226839929386 7043.79618964132169 27817.01735319429281 +23.50000000000006 1021238.99677721585613 4807.36348247983096 7019.22578394984976 27934.41395635497975 +23.60000000000007 1021166.03003201726824 4787.88092993845839 6994.68798529068226 28051.40105275414317 +23.70000000000007 1021093.32355904695578 4768.51276964074168 6970.18448547055141 28167.97918584231957 +23.80000000000007 1021020.87694689375348 4749.25721314701786 6945.71691269297935 28284.14892726682956 +23.90000000000007 1020948.68976752448361 4730.11252303269976 6921.28683363169330 28399.91087581171087 +24.00000000000007 1020876.76157692843117 4711.07701126278062 6896.89575543717820 28515.26565637223757 +24.10000000000007 1020805.09191574051511 4692.14903761877667 6872.54512767853521 28630.21391896285786 +24.20000000000007 1020733.68030984408688 4673.32700817642308 6848.23634422271516 28744.75633775749884 +24.30000000000008 1020662.52627095382195 4654.60937383247528 6823.97074505316505 28858.89361016121256 +24.40000000000008 1020591.62929717975203 4635.99462887903337 6799.74961802982671 28972.62645591209730 +24.50000000000008 1020520.98887357185595 4617.48130962385585 6775.57420059238757 29085.95561621259549 +24.60000000000008 1020450.60447264777031 4599.06799305516870 6751.44568140861429 29198.88185288913519 +24.70000000000008 1020380.47555490233935 4580.75329554953896 6727.36520196953188 29311.40594757927829 +24.80000000000008 1020310.60156930063386 4562.53587162142321 6703.33385813317182 29423.52870094543687 +24.90000000000008 1020240.98195375478826 4544.41441271303393 6679.35270161854351 29535.25093191432461 +25.00000000000009 1020171.61613558477256 4526.38764602323317 6655.42274145143347 29646.57347694129930 +25.10000000000009 1020102.50353196414653 4508.45433337418672 6631.54494536358743 29757.49718929882147 +25.20000000000009 1020033.64355035114568 4490.61327011455796 6607.72024114677515 29868.02293838821424 +25.30000000000009 1019965.03558890544809 4472.86328405806125 6583.94951796319856 29978.15160907399331 +25.40000000000009 1019896.67903689073864 4455.20323445623762 6560.23362761364660 30087.88410104004652 +25.50000000000009 1019828.57327506458387 4437.63201100434071 6536.57338576475740 30197.22132816693920 +25.60000000000009 1019760.71767605491914 4420.14853287926690 6512.96957313671010 30306.16421792968686 +25.70000000000010 1019693.11160472419579 4402.75174780849193 6489.42293665262241 30414.71371081530015 +25.80000000000010 1019625.75441852118820 4385.44063116901725 6465.93419055088270 30522.87075975950938 +25.90000000000010 1019558.64546782162506 4368.21418511535103 6442.50401746161879 30630.63632960202449 +26.00000000000010 1019491.78409625682980 4351.07143773558710 6419.13306944845135 30738.01139655971929 +26.10000000000010 1019425.16964103200007 4334.01144223468236 6395.82196901666157 30844.99694771719442 +26.20000000000010 1019358.80143323354423 4317.03327614403952 6372.57131008884608 30951.59398053414043 +26.30000000000010 1019292.67879812594038 4300.13604055655833 6349.38165894910981 31057.80350236895538 +26.40000000000011 1019226.80105543928221 4283.31885938632604 6326.25355515681531 31163.62653001810759 +26.50000000000011 1019161.16751964681316 4266.58087865215566 6303.18751243086172 31269.06408927071971 +26.60000000000011 1019095.77750023303088 4249.92126578419538 6280.18401950545285 31374.11721447790114 +26.70000000000011 1019030.63030195317697 4233.33920895287247 6257.24354095826311 31478.78694813632683 +26.80000000000011 1018965.72522508364636 4216.83391641944399 6234.36651801189873 31583.07434048563300 +26.90000000000011 1018901.06156566448044 4200.40461590745781 6211.55336930951034 31686.98044911916440 +27.00000000000011 1018836.63861573312897 4184.05055399445155 6188.80449166539256 31790.50633860765447 +27.10000000000012 1018772.45566355064511 4167.77099552323216 6166.12026079137559 31893.65308013540925 +27.20000000000012 1018708.51199382019695 4151.56522303210750 6143.50103199978639 31996.42175114859856 +27.30000000000012 1018644.80688789824490 4135.43253620345331 6120.94714088374076 32098.81343501526135 +27.40000000000012 1018581.33962399850134 4119.37225133002994 6098.45890397548828 32200.82922069665801 +27.50000000000012 1018518.10947738913819 4103.38370079846936 6076.03661938352525 32302.47020242958388 +27.60000000000012 1018455.11572058289312 4087.46623258938280 6053.68056740915472 32403.73747941930924 +27.70000000000012 1018392.35762352123857 4071.61920979354818 6031.39101114315690 32504.63215554279668 +27.80000000000013 1018329.83445375203155 4055.84201014366272 6009.16819704321097 32605.15533906185010 +27.90000000000013 1018267.54547660099342 4040.13402556115580 5987.01235549268949 32705.30814234590434 +28.00000000000013 1018205.48995533760171 4024.49466171757422 5964.92370134142311 32805.09168160411355 +28.10000000000013 1018143.66715133516118 4008.92333761007330 5942.90243442901919 32904.50707662646892 +28.20000000000013 1018082.07632422528695 3993.41948515055583 5920.94874009129308 33003.55545053361857 +28.30000000000013 1018020.71673204726540 3977.98254876801866 5899.06278965035926 33102.23792953514203 +28.40000000000013 1017959.58763139217626 3962.61198502368143 5877.24474088890474 33200.55564269598108 +28.50000000000014 1017898.68827754235826 3947.30726223848387 5855.49473850916002 33298.50972171079775 +28.60000000000014 1017838.01792460528668 3932.06786013255305 5833.81291457705538 33396.10130068595026 +28.70000000000014 1017777.57582564360928 3916.89326947625341 5812.19938895203813 33493.33151592889772 +28.80000000000014 1017717.36123280052561 3901.78299175244774 5790.65426970302178 33590.20150574476429 +28.90000000000014 1017657.37339742039330 3886.73653882960571 5769.17765351090293 33686.71241023981565 +29.00000000000014 1017597.61157016560901 3871.75343264540788 5747.76962605808512 33782.86537113166560 +29.10000000000014 1017538.07500112883281 3856.83320490051210 5726.43026240542622 33878.66153156596556 +29.20000000000014 1017478.76293994218577 3841.97539676214819 5705.15962735701214 33974.10203593938786 +29.30000000000015 1017419.67463588167448 3827.17955857722700 5683.95777581315451 34069.18802972866979 +29.40000000000015 1017360.80933796847239 3812.44524959465753 5662.82475311198414 34163.92065932555852 +29.50000000000015 1017302.16629506670870 3797.77203769657262 5641.76059536001503 34258.30107187742396 +29.60000000000015 1017243.74475597706623 3783.15949913817622 5620.76532975202554 34352.33041513342323 +29.70000000000015 1017185.54396952816751 3768.60721829593604 5599.83897488061029 34446.00983729595464 +29.80000000000015 1017127.56318466376979 3754.11478742385043 5578.98154103572688 34539.34048687729955 +29.90000000000015 1017069.80165052728262 3739.68180641752861 5558.19303049456448 34632.32351256122638 +30.00000000000016 1017012.25861654325854 3725.30788258583425 5537.47343780204301 34724.96006306946947 +30.10000000000016 1016954.93333249562420 3710.99263042984603 5516.82275004224903 34817.25128703283553 +30.20000000000016 1016897.82504860369954 3696.73567142890124 5496.24094710109330 34909.19833286687208 +30.30000000000016 1016840.93301559472457 3682.53663383349249 5475.72800192047725 35000.80234865188686 +30.40000000000016 1016784.25648477429058 3668.39515246479641 5455.28388074424129 35092.06448201723106 +30.50000000000016 1016727.79470809409395 3654.31086852062026 5434.90854335616041 35182.98588002963515 +30.60000000000016 1016671.54693821712863 3640.28342938756214 5414.60194331023649 35273.56768908556842 +30.70000000000017 1016615.51242858031765 3626.31248845918071 5394.36402815354450 35363.81105480740371 +30.80000000000017 1016559.69043345528189 3612.39770495998209 5374.19473964186727 35453.71712194329302 +30.90000000000017 1016504.08020800643135 3598.53874377503553 5354.09401394834822 35543.28703427065921 +31.00000000000017 1016448.68100834696088 3584.73527528503882 5334.06178186539546 35632.52193450312916 +31.10000000000017 1016393.49209159298334 3570.98697520665519 5314.09796900004312 35721.42296420088678 +31.20000000000017 1016338.51271591545083 3557.29352443795187 5294.20249596299072 35809.99126368422003 +31.30000000000017 1016283.74214059009682 3543.65460890877648 5274.37527855151711 35898.22797195026942 +31.40000000000018 1016229.17962604551576 3530.06991943591220 5254.61622792646813 35986.13422659279604 +31.50000000000018 1016174.82443390937988 3516.53915158285645 5234.92525078351264 36073.71116372490360 +31.60000000000018 1016120.67582705314271 3503.06200552407608 5215.30224951884247 36160.95991790462722 +31.70000000000018 1016066.73306963429786 3489.63818591359177 5195.74712238950451 36247.88162206327252 +31.80000000000018 1016012.99542713793926 3476.26740175775467 5176.25976366853047 36334.47740743643226 +31.90000000000018 1015959.46216641599312 3462.94936629207859 5156.84006379503717 36420.74840349757142 +32.00000000000018 1015906.13255572505295 3449.68379686199842 5137.48790951945466 36506.69573789415881 +32.10000000000019 1015853.00586476305034 3436.47041480742564 5118.20318404404043 36592.32053638614889 +32.20000000000019 1015800.08136470394675 3423.30894535098150 5098.98576715883428 36677.62392278688640 +32.30000000000019 1015747.35832823149394 3410.19911748978757 5079.83553537319312 36762.60701890620112 +32.40000000000019 1015694.83602957113180 3397.14066389069831 5060.75236204305929 36847.27094449575088 +32.50000000000019 1015642.51374452118762 3384.13332078886651 5041.73611749408610 36931.61681719646731 +32.60000000000019 1015590.39075048232917 3371.17682788953425 5022.78666914076530 37015.64575248803885 +32.70000000000019 1015538.46632648562081 3358.27092827294200 5003.90388160167913 37099.35886364038743 +32.80000000000020 1015486.73975322023034 3345.41536830226096 4985.08761681100259 37182.75726166708046 +32.90000000000020 1015435.21031305915676 3332.60989753444483 4966.33773412637584 37265.84205528059829 +33.00000000000020 1015383.87729008402675 3319.85426863391285 4947.65409043326508 37348.61435084936966 +33.10000000000020 1015332.73997010907624 3307.14823728896863 4929.03654024592652 37431.07525235658977 +33.20000000000020 1015281.79764070396777 3294.49156213086872 4910.48493580507693 37513.22586136069003 +33.30000000000020 1015231.04959121532738 3281.88400465545374 4891.99912717238112 37595.06727695744485 +33.40000000000020 1015180.49511278781574 3269.32532914726289 4873.57896232185612 37676.60059574365005 +33.50000000000021 1015130.13349838391878 3256.81530260604859 4855.22428722829864 37757.82691178235109 +33.60000000000021 1015079.96404280269053 3244.35369467561304 4836.93494595281481 37838.74731656949007 +33.70000000000021 1015029.98604269814678 3231.94027757489675 4818.71078072556793 37919.36289900203701 +33.80000000000021 1014980.19879659614526 3219.57482603124208 4800.55163202581298 37999.67474534745998 +33.90000000000021 1014930.60160491103306 3207.25711721576317 4782.45733865931652 38079.68393921455572 +34.00000000000021 1014881.19376996112987 3194.98693068075363 4764.42773783324628 38159.39156152554642 +34.10000000000021 1014831.97459598351270 3182.76404829906869 4746.46266522860424 38238.79869048943510 +34.20000000000022 1014782.94338914833497 3170.58825420541552 4728.56195507028951 38317.90640157657617 +34.30000000000022 1014734.09945757186506 3158.45933473949117 4710.72544019486304 38396.71576749441738 +34.40000000000022 1014685.44211132929195 3146.37707839090945 4692.95295211609300 38475.22785816433316 +34.50000000000022 1014636.97066246683244 3134.34127574585773 4675.24432108834208 38553.44374069960031 +34.60000000000022 1014588.68442501290701 3122.35171943542991 4657.59937616788011 38631.36447938440688 +34.70000000000022 1014540.58271498896647 3110.40820408557875 4640.01794527217407 38708.99113565387233 +34.80000000000022 1014492.66485041962005 3098.51052626863702 4622.49985523723262 38786.32476807507192 +34.90000000000023 1014444.93015134206507 3086.65848445635902 4605.04493187305980 38863.36643232902861 +35.00000000000023 1014397.37793981516734 3074.85187897442893 4587.65300001728519 38940.11718119357829 +35.10000000000023 1014350.00753992784303 3063.09051195839311 4570.32388358701792 39016.57806452720251 +35.20000000000023 1014302.81827780685853 3051.37418731096659 4553.05740562899791 39092.75012925365445 +35.30000000000023 1014255.80948162428103 3039.70271066067380 4535.85338836808478 39168.63441934747243 +35.40000000000023 1014208.98048160434701 3028.07588932177669 4518.71165325414222 39244.23197582027205 +35.50000000000023 1014162.33061002986506 3016.49353225545201 4501.63202100737635 39319.54383670783864 +35.60000000000024 1014115.85920124826953 3004.95545003217421 4484.61431166216607 39394.57103705796180 +35.70000000000024 1014069.56559167685919 2993.46145479526831 4467.65834460944097 39469.31460891899769 +35.80000000000024 1014023.44911980815232 2982.01136022559285 4450.76393863765406 39543.77558132915146 +35.90000000000024 1013977.50912621442694 2970.60498150731746 4433.93091197239028 39617.95498030644376 +36.00000000000024 1013931.74495355179533 2959.24213529476174 4417.15908231465801 39691.85382883931743 +36.10000000000024 1013886.15594656451140 2947.92263968025554 4400.44826687790282 39765.47314687789185 +36.20000000000024 1013840.74145208788104 2936.64631416299517 4383.79828242378881 39838.81395132585749 +36.30000000000025 1013795.50081905198749 2925.41297961885857 4367.20894529678299 39911.87725603292347 +36.40000000000025 1013750.43339848390315 2914.22245827114966 4350.68007145758384 39984.66407178786903 +36.50000000000025 1013705.53854351071641 2903.07457366224253 4334.21147651542833 40057.17540631216252 +36.60000000000025 1013660.81560936104506 2891.96915062609605 4317.80297575931672 40129.41226425408968 +36.70000000000025 1013616.26395336736459 2880.90601526161163 4301.45438418818867 40201.37564718341309 +36.80000000000025 1013571.88293496717233 2869.88499490680670 4285.16551654008344 40273.06655358654825 +36.90000000000025 1013527.67191570426803 2858.90591811377681 4268.93618732031518 40344.48597886221978 +37.00000000000026 1013483.63025922991801 2847.96861462442166 4252.76621082870315 40415.63491531756154 +37.10000000000026 1013439.75733130308799 2837.07291534691103 4236.65540118587433 40486.51435216470418 +37.20000000000026 1013396.05249979125801 2826.21865233286644 4220.60357235867832 40557.12527551779931 +37.30000000000026 1013352.51513467018958 2815.40565875523544 4204.61053818474284 40627.46866839044378 +37.40000000000026 1013309.14460802404210 2804.63376888683570 4188.67611239618782 40697.54551069352601 +37.50000000000026 1013265.94029404502362 2793.90281807954761 4172.80010864253654 40767.35677923346520 +37.60000000000026 1013222.90156903269235 2783.21264274413716 4156.98234051284453 40836.90344771084347 +37.70000000000027 1013180.02781139337458 2772.56308033068399 4141.22262155706903 40906.18648671939445 +37.80000000000027 1013137.31840163888410 2761.95396930959851 4125.52076530670820 40975.20686374534853 +37.90000000000027 1013094.77272238547448 2751.38514915321002 4109.87658529472992 41043.96554316712718 +38.00000000000027 1013052.39015835244209 2740.85646031790384 4094.28989507481538 41112.46348625537212 +38.10000000000027 1013010.17009636049625 2730.36774422679264 4078.76050823993819 41180.70165117328725 +38.20000000000027 1012968.11192533001304 2719.91884325290448 4063.28823844030057 41248.68099297728622 +38.30000000000027 1012926.21503627905622 2709.50960070287101 4047.87289940064647 41316.40246361796017 +38.40000000000028 1012884.47882232116535 2699.13986080109726 4032.51430493697308 41383.86701194130728 +38.50000000000028 1012842.90267866326030 2688.80946867440252 4017.21226897266024 41451.07558369025355 +38.60000000000028 1012801.48600260296371 2678.51827033711197 4001.96660555403378 41518.02912150646443 +38.70000000000028 1012760.22819352627266 2668.26611267658882 3986.77712886538529 41584.72856493236759 +38.80000000000028 1012719.12865290453192 2658.05284343918947 3971.64365324346090 41651.17485041345935 +38.90000000000028 1012678.18678429175634 2647.87831121663066 3956.56599319143879 41717.36891130085132 +39.00000000000028 1012637.40199332148768 2637.74236543275401 3941.54396339241430 41783.31167785404250 +39.10000000000029 1012596.77368770365138 2627.64485633067579 3926.57737872240159 41849.00407724391698 +39.20000000000029 1012556.30127722152974 2617.58563496030956 3911.66605426287470 41914.44703355595993 +39.30000000000029 1012515.98417372792028 2607.56455316624897 3896.80980531285832 41979.64146779367729 +39.40000000000029 1012475.82179114187602 2597.58146357600071 3882.00844740058483 42044.58829788222647 +39.50000000000029 1012435.81354544521309 2587.63621958855401 3867.26179629472881 42109.28843867223623 +39.60000000000029 1012395.95885467843618 2577.72867536327931 3852.56966801523777 42173.74280194381572 +39.70000000000029 1012356.25713893712964 2567.85868580914257 3837.93187884376493 42237.95229641073820 +39.80000000000030 1012316.70782036799937 2558.02610657422611 3823.34824533372421 42301.91782772479928 +39.90000000000030 1012277.31032316491473 2548.23079403554675 3808.81858431997398 42365.64029848036444 +40.00000000000030 1012238.06407356448472 2538.47260528916013 3794.34271292814537 42429.12060821903287 +40.10000000000030 1012198.96849984209985 2528.75139814054546 3779.92044858362169 42492.35965343450516 +40.20000000000030 1012160.02303230774123 2519.06703109525733 3765.55160902018724 42555.35832757756725 +40.30000000000030 1012121.22710330132395 2509.41936334983848 3751.23601228834923 42618.11752106123458 +40.40000000000030 1012082.58014718838967 2499.80825478298675 3736.97347676334812 42680.63812126604171 +40.50000000000031 1012044.08160035556648 2490.23356594696361 3722.76382115286242 42742.92101254543377 +40.60000000000031 1012005.73090120579582 2480.69515805924084 3708.60686450442290 42804.96707623131806 +40.70000000000031 1011967.52749015414156 2471.19289299437651 3694.50242621253983 42866.77719063972472 +40.80000000000031 1011929.47080962255131 2461.72663327611099 3680.45032602555602 42928.35223107659840 +40.90000000000031 1011891.56030403519981 2452.29624206967810 3666.45038405223477 42989.69306984369177 +41.00000000000031 1011853.79541981383227 2442.90158317432588 3652.50242076808854 43050.80057624456094 +41.10000000000031 1011816.17560537264217 2433.54252101603515 3638.60625702145990 43111.67561659069906 +41.20000000000032 1011778.70031111326534 2424.21892064043550 3624.76171403935950 43172.31905420772091 +41.30000000000032 1011741.36898942012340 2414.93064770590809 3610.96861343307091 43232.73174944170751 +41.40000000000032 1011704.18109465483576 2405.67756847687178 3597.22677720353067 43292.91455966559442 +41.50000000000032 1011667.13608315144666 2396.45954981724663 3583.53602774648880 43352.86833928565466 +41.60000000000032 1011630.23341321118642 2387.27645918408643 3569.89618785745597 43412.59393974809791 +41.70000000000032 1011593.47254509723280 2378.12816462137926 3556.30708073644882 43472.09220954572083 +41.80000000000032 1011556.85294102958869 2369.01453475400604 3542.76852999253470 43531.36399422465911 +41.90000000000033 1011520.37406517949421 2359.93543878185665 3529.28035964818491 43590.41013639119774 +42.00000000000033 1011484.03538366453722 2350.89074647409552 3515.84239414344302 43649.23147571866866 +42.10000000000033 1011447.83636454283260 2341.88032816357190 3502.45445833991289 43707.82884895439202 +42.20000000000033 1011411.77647780801635 2332.90405474137287 3489.11637752457273 43766.20308992672653 +42.30000000000033 1011375.85519538365770 2323.96179765151101 3475.82797741342029 43824.35502955213451 +42.40000000000033 1011340.07199111767113 2315.05342888574569 3462.58908415495625 43882.28549584235589 +42.50000000000033 1011304.42634077707771 2306.17882097853226 3449.39952433350982 43939.99531391160417 +42.60000000000034 1011268.91772204241715 2297.33784700209435 3436.25912497241006 43997.48530598382786 +42.70000000000034 1011233.54561450204346 2288.53038056161677 3423.16771353701279 44054.75629140003730 +42.80000000000034 1011198.30949964688625 2279.75629579055521 3410.12511793758085 44111.80908662565344 +42.90000000000034 1011163.20886086462997 2271.01546734605699 3397.13116653202906 44168.64450525794382 +43.00000000000034 1011128.24318343412597 2262.30777040449357 3384.18568812853482 44225.26335803348047 +43.10000000000034 1011093.41195451992098 2253.63308065709543 3371.28851198801749 44281.66645283561957 +43.20000000000034 1011058.71466316643637 2244.99127430569115 3358.43946782649482 44337.85459470208298 +43.30000000000035 1011024.15080029226374 2236.38222805854548 3345.63838581731670 44393.82858583252528 +43.40000000000035 1010989.71985868492629 2227.80581912629350 3332.88509659328201 44449.58922559614439 +43.50000000000035 1010955.42133299470879 2219.26192521796656 3320.17943124864314 44505.13731053936499 +43.60000000000035 1010921.25471972906962 2210.75042453710967 3307.52122134099864 44560.47363439350738 +43.70000000000035 1010887.21951724705286 2202.27119577798476 3294.91029889308038 44615.59898808252183 +43.80000000000035 1010853.31522575358395 2193.82411812186047 3282.34649639443887 44670.51415973073745 +43.90000000000035 1010819.54134729353245 2185.40907123338138 3269.82964680302894 44725.21993467064749 +44.00000000000036 1010785.89738574612420 2177.02593525701832 3257.35958354669719 44779.71709545069461 +44.10000000000036 1010752.38284681923687 2168.67459081359357 3244.93614052457951 44834.00642184313620 +44.20000000000036 1010718.99723804334644 2160.35491899688350 3232.55915210840567 44888.08869085188053 +44.30000000000036 1010685.74006876617204 2152.06680137029070 3220.22845314371853 44941.96467672035214 +44.40000000000036 1010652.61085014650598 2143.81011996358575 3207.94387895100817 44995.63515093941533 +44.50000000000036 1010619.60909514874220 2135.58475726971983 3195.70526532676558 45049.10088225526852 +44.60000000000036 1010586.73431853693910 2127.39059624169931 3183.51244854445531 45102.36263667738240 +44.70000000000037 1010553.98603686911520 2119.22752028952618 3171.36526535541361 45155.42117748645978 +44.80000000000037 1010521.36376849131193 2111.09541327719853 3159.26355298967337 45208.27726524238096 +44.90000000000037 1010488.86703353188932 2102.99415951977244 3147.20714915671488 45260.93165779220726 +45.00000000000037 1010456.49535389582161 2094.92364378047978 3135.19589204614977 45313.38511027815548 +45.10000000000037 1010424.24825325876009 2086.88375126790379 3123.22962032833857 45365.63837514559418 +45.20000000000037 1010392.12525706132874 2078.87436763320829 3111.30817315494141 45417.69220215106907 +45.30000000000037 1010360.12589250341989 2070.89537896741922 3099.43139015940778 45469.54733837032109 +45.40000000000038 1010328.24968853814062 2062.94667179875705 3087.59911145740625 45521.20452820631181 +45.50000000000038 1010296.49617586610839 2055.02813309002067 3075.81117764719465 45572.66451339727064 +45.60000000000038 1010264.86488692986313 2047.13965023601622 3064.06742980993431 45623.92803302472748 +45.70000000000038 1010233.35535590804648 2039.28111106103484 3052.36770950994833 45674.99582352155994 +45.80000000000038 1010201.96711870923173 2031.45240381637450 3040.71185879493078 45725.86861868006235 +45.90000000000038 1010170.69971296656877 2023.65341717790579 3029.09972019609950 45776.54714965997846 +46.00000000000038 1010139.55267803196330 2015.88404024367969 3017.53113672830386 45827.03214499657770 +46.10000000000039 1010108.52555497013964 2008.14416253157719 3006.00595189008254 45877.32433060871699 +46.20000000000039 1010077.61788655293640 2000.43367397699808 2994.52400966367577 45927.42442980688793 +46.30000000000039 1010046.82921725360211 1992.75246493058899 2983.08515451499261 45977.33316330128582 +46.40000000000039 1010016.15909324109089 1985.10042615600855 2971.68923139353592 46027.05124920987146 +46.50000000000039 1009985.60706237400882 1977.47744882772918 2960.33608573228503 46076.57940306643286 +46.60000000000039 1009955.17267419537529 1969.88342452887355 2949.02556344753657 46125.91833782863978 +46.70000000000039 1009924.85547992656939 1962.31824524908643 2937.75751093870986 46175.06876388609817 +46.80000000000040 1009894.65503246150911 1954.78180338243851 2926.53177508811132 46224.03138906841195 +46.90000000000040 1009864.57088636117987 1947.27399172536366 2915.34820326066392 46272.80691865321569 +47.00000000000040 1009834.60259784793016 1939.79470347462620 2904.20664330360205 46321.39605537422904 +47.10000000000040 1009804.74972479965072 1932.34383222531937 2893.10694354613088 46369.79949942928943 +47.20000000000040 1009775.01182674407028 1924.92127196889305 2882.04895279905395 46418.01794848839199 +47.30000000000040 1009745.38846485316753 1917.52691709121018 2871.03252035436890 46466.05209770170768 +47.40000000000040 1009715.87920193735044 1910.16066237063046 2860.05749598483226 46513.90263970761589 +47.50000000000041 1009686.48360244010109 1902.82240297612157 2849.12372994349471 46561.57026464069349 +47.60000000000041 1009657.20123243203852 1895.51203446539625 2838.23107296320813 46609.05566013975476 +47.70000000000041 1009628.03165960544720 1888.22945278307475 2827.37937625610448 46656.35951135581126 +47.80000000000041 1009598.97445326845627 1880.97455425887188 2816.56849151304914 46703.48250096008269 +47.90000000000041 1009570.02918433956802 1873.74723560580810 2805.79827090306617 46750.42530915196403 +48.00000000000041 1009541.19542534218635 1866.54739391844441 2795.06856707274210 46797.18861366701458 +48.10000000000041 1009512.47275039879605 1859.37492667113861 2784.37923314560248 46843.77308978489600 +48.20000000000041 1009483.86073522537481 1852.22973171632430 2773.73012272146707 46890.17941033732495 +48.30000000000042 1009455.35895712592173 1845.11170728281149 2763.12108987578222 46936.40824571601843 +48.40000000000042 1009426.96699498686939 1838.02075197410727 2752.55198915893243 46982.46026388061728 +48.50000000000042 1009398.68442927161232 1830.95676476675681 2742.02267559552911 47028.33613036660245 +48.60000000000042 1009370.51084201491904 1823.91964500870381 2731.53300468368207 47074.03650829319668 +48.70000000000042 1009342.44581681734417 1816.90929241767094 2721.08283239424964 47119.56205837125890 +48.80000000000042 1009314.48893883975688 1809.92560707955704 2710.67201517006970 47164.91343891116412 +48.90000000000042 1009286.63979479786940 1802.96848944685371 2700.30040992517525 47210.09130583066872 +49.00000000000043 1009258.89797295676544 1796.03784033707893 2689.96787404398992 47255.09631266275392 +49.10000000000043 1009231.26306312531233 1789.13356093122729 2679.67426538050813 47299.92911056348385 +49.20000000000043 1009203.73465665103868 1782.25555277223748 2669.41944225745920 47344.59034831982717 +49.30000000000043 1009176.31234641419724 1775.40371776347524 2659.20326346545471 47389.08067235744966 +49.40000000000043 1009148.99572682264261 1768.57795816723274 2649.02558826212316 47433.40072674854309 +49.50000000000043 1009121.78439380647615 1761.77817660324195 2638.88627637122681 47477.55115321958147 +49.60000000000043 1009094.67794481245801 1755.00427604720403 2628.78518798176856 47521.53259115909896 +49.70000000000044 1009067.67597879865207 1748.25615982933255 2618.72218374708291 47565.34567762546067 +49.80000000000044 1009040.77809622907080 1741.53373163291076 2608.69712478391375 47608.99104735457513 +49.90000000000044 1009013.98389906843659 1734.83689549286328 2598.70987267148121 47652.46933276764321 +50.00000000000000 1008987.29299077682663 1728.16555579436863 2588.76028945058079 47695.78116397864505 \ No newline at end of file diff --git a/cpp/tests/data/ode-seir-metapop-compare.csv b/cpp/tests/data/ode-seir-metapop-compare.csv new file mode 100644 index 0000000000..6b6ed790b5 --- /dev/null +++ b/cpp/tests/data/ode-seir-metapop-compare.csv @@ -0,0 +1,26 @@ + # t S_0 E_0 I_0 R_0 S_1 E_1 I_1 R_1 S_2 E_2 I_2 R_2 S_3 E_3 I_3 R_3 +0.000000000000000000 1049000.000000000000000000 10000.000000000000000000 1000.000000000000000000 1000.000000000000000000 1049000.000000000000000000 10000.000000000000000000 1000.000000000000000000 1000.000000000000000000 1049000.000000000000000000 10000.000000000000000000 1000.000000000000000000 1000.000000000000000000 1049000.000000000000000000 10000.000000000000000000 1000.000000000000000000 1000.000000000000000000 +0.100000000000000006 1048970.988352040527388453 9838.271405749961559195 1172.626704245788005210 1018.113537963836847666 1048970.988352040527388453 9838.271405749961559195 1172.626704245788005210 1018.113537963836847666 1048970.988352040527388453 9838.271405749961559195 1172.626704245788005210 1018.113537963836847666 1048970.988352040527388453 9838.271405749961559195 1172.626704245788005210 1018.113537963836847666 +0.550000000000000044 1048786.562808100134134293 9199.938732493899806286 1880.226382551424649137 1133.272076854476154040 1048786.562808100134134293 9199.938732493899806286 1880.226382551424649137 1133.272076854476154040 1048786.562808100134134293 9199.938732493899806286 1880.226382551424649137 1133.272076854476154040 1048786.562808100134134293 9199.938732493899806286 1880.226382551424649137 1133.272076854476154040 +1.571335751995554197 1048092.401353218709118664 8194.697825347475372837 3146.000554222973278229 1566.900267210686934050 1048092.401353218709118664 8194.697825347475372837 3146.000554222973278229 1566.900267210686934050 1048092.401353218709118664 8194.697825347475372837 3146.000554222973278229 1566.900267210686934050 1048092.401353218709118664 8194.697825347475372837 3146.000554222973278229 1566.900267210686934050 +2.721576076432889124 1046961.386030971654690802 7588.350638500698551070 4176.225815231045999099 2274.037515296420224331 1046961.386030971654690802 7588.350638500698551070 4176.225815231045999099 2274.037515296420224331 1046961.386030971654690802 7588.350638500698551070 4176.225815231045999099 2274.037515296420224331 1046961.386030971654690802 7588.350638500698551070 4176.225815231045999099 2274.037515296420224331 +4.029943604597698403 1045349.619403100456111133 7330.778352213967991702 5036.529414663581519562 3283.072830021797926747 1045349.619403100456111133 7330.778352213967991702 5036.529414663581519562 3283.072830021797926747 1045349.619403100456111133 7330.778352213967991702 5036.529414663581519562 3283.072830021797926747 1045349.619403100456111133 7330.778352213967991702 5036.529414663581519562 3283.072830021797926747 +5.507034559643384952 1043220.440254993853159249 7377.465553596446625306 5783.675800650045857765 4618.418390759416979563 1043220.440254993853159249 7377.465553596446625306 5783.675800650045857765 4618.418390759416979563 1043220.440254993853159249 7377.465553596446625306 5783.675800650045857765 4618.418390759416979563 1043220.440254993853159249 7377.465553596446625306 5783.675800650045857765 4618.418390759416979563 +7.181045772330549859 1040493.505196049343794584 7685.356236648490266816 6488.497226337810388941 6332.641340964195478591 1040493.505196049343794584 7685.356236648490266816 6488.497226337810388941 6332.641340964195478591 1040493.505196049343794584 7685.356236648490266816 6488.497226337810388941 6332.641340964195478591 1040493.505196049343794584 7685.356236648490266816 6488.497226337810388941 6332.641340964195478591 +9.068602661997122283 1037073.405391793348826468 8222.868082855427928735 7214.758721446800336707 8488.967803904286483885 1037073.405391793348826468 8222.868082855427928735 7214.758721446800336707 8488.967803904286483885 1037073.405391793348826468 8222.868082855427928735 7214.758721446800336707 8488.967803904286483885 1037073.405391793348826468 8222.868082855427928735 7214.758721446800336707 8488.967803904286483885 +11.248067398941367756 1032694.429363684146665037 8997.711351056386774871 8047.601685562283819309 11260.257599697049954557 1032694.429363684146665037 8997.711351056386774871 8047.601685562283819309 11260.257599697049954557 1032694.429363684146665037 8997.711351056386774871 8047.601685562283819309 11260.257599697049954557 1032694.429363684146665037 8997.711351056386774871 8047.601685562283819309 11260.257599697049954557 +11.624942631692711359 1031892.499271133681759238 9142.254586236071190797 8194.879248038520017872 11770.366894591707023210 1031892.499271133681759238 9142.254586236071190797 8194.879248038520017872 11770.366894591707023210 1031892.499271133681759238 9142.254586236071190797 8194.879248038520017872 11770.366894591707023210 1031892.499271133681759238 9142.254586236071190797 8194.879248038520017872 11770.366894591707023210 +12.001817864444054962 1031192.262633586768060923 9177.270036856716615148 8340.744771348789072363 12289.722558207797192154 1031192.262633586768060923 9177.270036856716615148 8340.744771348789072363 12289.722558207797192154 1031192.262633586768060923 9177.270036856716615148 8340.744771348789072363 12289.722558207797192154 1031192.262633586768060923 9177.270036856716615148 8340.744771348789072363 12289.722558207797192154 +12.529611483968700725 1030606.423360764631070197 8845.787659644680388737 8516.279530909305321984 13031.509448681445064722 1030606.423360764631070197 8845.787659644680388737 8516.279530909305321984 13031.509448681445064722 1030606.423360764631070197 8845.787659644680388737 8516.279530909305321984 13031.509448681445064722 1030606.423360764631070197 8845.787659644680388737 8516.279530909305321984 13031.509448681445064722 +13.057405103493346488 1030131.377186079975217581 8443.798219163865724113 8638.425499267183113261 13786.399095489065075526 1030131.377186079975217581 8443.798219163865724113 8638.425499267183113261 13786.399095489065075526 1030131.377186079975217581 8443.798219163865724113 8638.425499267183113261 13786.399095489065075526 1030131.377186079975217581 8443.798219163865724113 8638.425499267183113261 13786.399095489065075526 +15.002878194169083415 1028352.653943004552274942 7292.666872994709592604 8738.639010714050527895 16616.040173286801291397 1028352.653943004552274942 7292.666872994709592604 8738.639010714050527895 16616.040173286801291397 1028352.653943004552274942 7292.666872994709592604 8738.639010714050527895 16616.040173286801291397 1028352.653943004552274942 7292.666872994709592604 8738.639010714050527895 16616.040173286801291397 +16.948351284844822118 1026594.450561258592642844 6481.962568733895750483 8505.735030830541290925 19417.851839176964858780 1026594.450561258592642844 6481.962568733895750483 8505.735030830541290925 19417.851839176964858780 1026594.450561258592642844 6481.962568733895750483 8505.735030830541290925 19417.851839176964858780 1026594.450561258592642844 6481.962568733895750483 8505.735030830541290925 19417.851839176964858780 +19.121289365575044883 1024712.823017425718717277 5803.026431943653733470 8062.481969843786828278 22421.668580786674283445 1024712.823017425718717277 5803.026431943653733470 8062.481969843786828278 22421.668580786674283445 1024712.823017425718717277 5803.026431943653733470 8062.481969843786828278 22421.668580786674283445 1024712.823017425718717277 5803.026431943653733470 8062.481969843786828278 22421.668580786674283445 +21.568220941483573938 1022730.819907458848319948 5201.036499887130958086 7476.448597236990281090 25591.694995416903111618 1022730.819907458848319948 5201.036499887130958086 7476.448597236990281090 25591.694995416903111618 1022730.819907458848319948 5201.036499887130958086 7476.448597236990281090 25591.694995416903111618 1022730.819907458848319948 5201.036499887130958086 7476.448597236990281090 25591.694995416903111618 +24.302867828470034794 1020699.985484831035137177 4646.677304804762570711 6807.136086702391366998 28846.201123661736346548 1020699.985484831035137177 4646.677304804762570711 6807.136086702391366998 28846.201123661736346548 1020699.985484831035137177 4646.677304804762570711 6807.136086702391366998 28846.201123661736346548 1020699.985484831035137177 4646.677304804762570711 6807.136086702391366998 28846.201123661736346548 +27.405783498438019308 1018625.438826579484157264 4112.061390983284582035 6085.048903788306233764 32177.450878648814978078 1018625.438826579484157264 4112.061390983284582035 6085.048903788306233764 32177.450878648814978078 1018625.438826579484157264 4112.061390983285491529 6085.048903788305324269 32177.450878648814978078 1018625.438826579484157264 4112.061390983284582035 6085.048903788306233764 32177.450878648814978078 +30.982097872048104392 1016514.052626653225161135 3582.158500640592137643 5328.958786543431415339 35574.830086162648512982 1016514.052626653225161135 3582.158500640592592390 5328.958786543431415339 35574.830086162648512982 1016514.052626653225161135 3582.158500640593047137 5328.958786543431415339 35574.830086162648512982 1016514.052626653225161135 3582.158500640592137643 5328.958786543431415339 35574.830086162648512982 +35.174157269526475034 1014377.523218245478346944 3050.869285028250033065 4551.750098480431915959 39019.857398245752847288 1014377.523218245478346944 3050.869285028250487812 4551.750098480431915959 39019.857398245752847288 1014377.523218245478346944 3050.869285028250487812 4551.750098480432825454 39019.857398245752847288 1014377.523218245478346944 3050.869285028250033065 4551.750098480431915959 39019.857398245752847288 +40.165160093022812760 1012242.354254877660423517 2520.413323501484228473 3767.284556028632778180 42469.947865592235757504 1012242.354254877660423517 2520.413323501484228473 3767.284556028634142422 42469.947865592235757504 1012242.354254877660423517 2520.413323501484683220 3767.284556028633687674 42469.947865592235757504 1012242.354254877660423517 2520.413323501484228473 3767.284556028632778180 42469.947865592235757504 +46.138354563962117538 1010170.613532148068770766 2004.483312102868012516 3000.391258669248145452 45824.511897079959453549 1010170.613532148068770766 2004.483312102867785143 3000.391258669249054947 45824.511897079959453549 1010170.613532148068770766 2004.483312102868694637 3000.391258669249054947 45824.511897079959453549 1010170.613532148068770766 2004.483312102868012516 3000.391258669248145452 45824.511897079959453549 +50.000000000000000000 1009063.586620887159369886 1728.076734493041385576 2588.504609791979419242 47619.832034827901225071 1009063.586620887159369886 1728.076734493041385576 2588.504609791979873989 47619.832034827901225071 1009063.586620887159369886 1728.076734493041840324 2588.504609791979873989 47619.832034827901225071 1009063.586620887159369886 1728.076734493041385576 2588.504609791979419242 47619.832034827901225071 \ No newline at end of file diff --git a/cpp/tests/test_odemetapop.cpp b/cpp/tests/test_odemetapop.cpp new file mode 100644 index 0000000000..deac64403a --- /dev/null +++ b/cpp/tests/test_odemetapop.cpp @@ -0,0 +1,309 @@ + +#include "load_test_data.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "ode_seir_metapop/model.h" +#include "ode_seir_metapop/parameters.h" +#include "memilio/math/euler.h" +#include "memilio/compartments/simulation.h" +#include +#include +#include + +class ModelTestOdeMetapop : public testing::Test +{ +public: + ModelTestOdeMetapop() + : model(4, 1) + { + } + ScalarType t0; + ScalarType tmax; + ScalarType dt; + ScalarType total_population_per_region; + mio::oseirmetapop::Model model; + +protected: + void SetUp() override + { + t0 = 0.; + tmax = 50.; + dt = 0.1; + + total_population_per_region = 1061000; + + for (size_t i = 0; i < (size_t)model.parameters.get_num_regions(); i++) { + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseir::InfectionState::Exposed)}] = 10000; + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseir::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseir::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible)}] = + total_population_per_region - + model.populations[{ + mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseir::InfectionState::Exposed)}] - + model.populations[{ + mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseir::InfectionState::Infected)}] - + model.populations[{ + mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseir::InfectionState::Recovered)}]; + } + model.set_commuting_strengths(); + } +}; + +TEST_F(ModelTestOdeMetapop, simulateDefault) +{ + mio::TimeSeries result = simulate(t0, tmax, dt, model); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); +} + +TEST_F(ModelTestOdeMetapop, checkPopulationConservation) +{ + auto result = simulate(t0, tmax, dt, model); + ScalarType num_persons = result.get_last_value().sum(); + EXPECT_NEAR(num_persons, total_population_per_region * (size_t)model.parameters.get_num_regions(), 1e-9); +} + +TEST_F(ModelTestOdeMetapop, compareWithPreviousRun) +{ + mio::oseir::Parameters local_params(1); + + local_params.set>(5.2); + local_params.set>(6); + local_params.set>(0.1); + mio::ContactMatrixGroup& contact_matrix = local_params.get>(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + model.set_local_parameters(local_params); + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + std::vector> refData = load_test_data_csv("ode-seir-metapop-compare.csv"); + auto result = mio::simulate>(t0, tmax, dt, model); + + ASSERT_EQ(refData.size(), static_cast(result.get_num_time_points())); + + for (Eigen::Index irow = 0; irow < result.get_num_time_points(); ++irow) { + double t = refData[static_cast(irow)][0]; + auto rel_tol = 1e-6; + + //test result diverges at damping because of changes, not worth fixing at the moment + if (t > 11.0 && t < 13.0) { + //strong divergence around damping + rel_tol = 0.5; + } + else if (t > 13.0) { + //minor divergence after damping + rel_tol = 1e-2; + } + mio::unused(rel_tol); + + ASSERT_NEAR(t, result.get_times()[irow], 1e-12) << "at row " << irow; + + for (size_t icol = 0; icol < 12; ++icol) { + double ref = refData[static_cast(irow)][icol + 1]; + double actual = result[irow][icol]; + + double tol = rel_tol * ref; + ASSERT_NEAR(ref, actual, tol) << "at row " << irow; + } + } +} + +TEST_F(ModelTestOdeMetapop, check_constraints_parameters) +{ + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + // model.check_constraints() combines the functions from population and parameters. + // We only want to test the functions for the parameters defined in parameters.h + ASSERT_EQ(model.parameters.check_constraints(), 0); + + mio::set_log_level(mio::LogLevel::off); + + mobility_data_commuter(0, 1) += 0.5; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 1) = -0.5; + mobility_data_commuter(0, 2) = 0.75; + mobility_data_commuter(0, 3) = 0.75; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 1) = 1.5; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 0.; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 0) = 0.; + mobility_data_commuter(0, 1) = 0.; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 1.; + model.set_commuting_strengths(mobility_data_commuter); + model.parameters.set>( + mio::Populations( + {mio::regions::Region(4), mio::AgeGroup(1)}, 0.)); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + // Nobody commutes to region 2 but everybody originating fron there commutes to other regions. + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0., 0.5, 0.5, 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mio::set_log_level(mio::LogLevel::warn); +} + +TEST_F(ModelTestOdeMetapop, apply_constraints_parameters) +{ + const ScalarType tol_times = 1e-1; + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + EXPECT_EQ(model.parameters.apply_constraints(), 0); + + mobility_data_commuter(0, 1) += 0.5; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>() + .get_cont_freq_mat()[0] + .get_baseline() + .isIdentity(), + true); + + mobility_data_commuter(0, 1) = -0.5; + mobility_data_commuter(0, 2) = 0.75; + mobility_data_commuter(0, 3) = 0.75; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>() + .get_cont_freq_mat()[0] + .get_baseline() + .isIdentity(), + true); + + mobility_data_commuter(0, 1) = 1.5; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 0.; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ(model.parameters.get>() + .get_cont_freq_mat()[0] + .get_baseline() + .isIdentity(), + true); + + mobility_data_commuter(0, 0) = 0.; + mobility_data_commuter(0, 1) = 0.; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 1.; + model.set_commuting_strengths(mobility_data_commuter); + model.parameters.set>( + mio::Populations( + {mio::regions::Region(4), mio::AgeGroup(1)}, 0.)); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR((model.parameters.get>()[{ + mio::regions::Region(3), mio::AgeGroup(0)}]), + 1.0, tol_times); + + // Nobody commutes to region 2 but everybody originating fron there commutes to other regions. + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0., 0.5, 0.5, 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR((model.parameters.get>()[{ + mio::regions::Region(1), mio::AgeGroup(0)}]), + 1.0, tol_times); + + mio::set_log_level(mio::LogLevel::warn); +} + +TEST(TestOdeMetapop, compareSEIR) +{ + ScalarType t0 = 0; + ScalarType tmax = 50.; + ScalarType dt = 0.1; + + ScalarType total_population = 1061000; + + mio::oseirmetapop::Model model(1, 1); + + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseir::InfectionState::Exposed)}] = 10000; + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseir::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseir::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible)}] = + total_population - + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseir::InfectionState::Exposed)}] - + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseir::InfectionState::Infected)}] - + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseir::InfectionState::Recovered)}]; + + mio::oseir::Parameters local_params(1); + + local_params.set>(5.2); + local_params.set>(6); + local_params.set>(0.1); + mio::ContactMatrixGroup& contact_matrix = local_params.get>(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + model.set_local_parameters(local_params); + + model.set_commuting_strengths(); + + model.check_constraints(); + + // Use the Euler integrator as adaptive methods make different time steps in this model due to restructured equations. + std::shared_ptr> integrator = std::make_shared>(); + + auto result = simulate(t0, tmax, dt, model, integrator); + std::vector> refData = load_test_data_csv("ode-seir-compare-euler.csv"); + + ASSERT_EQ(refData.size(), static_cast(result.get_num_time_points())); + + for (Eigen::Index irow = 0; irow < result.get_num_time_points(); ++irow) { + ScalarType t = refData[static_cast(irow)][0]; + auto rel_tol = 1e-10; + + //test result diverges at damping because of changes, not worth fixing at the moment + if (t > 11.0 && t < 13.0) { + //strong divergence around damping + rel_tol = 0.5; + } + else if (t > 13.0) { + //minor divergence after damping + rel_tol = 1e-2; + } + mio::unused(rel_tol); + + ASSERT_NEAR(t, result.get_times()[irow], 1e-12) << "at row " << irow; + + for (size_t icol = 0; icol < refData[static_cast(irow)].size() - 1; ++icol) { + ScalarType ref = refData[static_cast(irow)][icol + 1]; + ScalarType actual = result[irow][icol]; + + ScalarType tol = rel_tol * ref; + ASSERT_NEAR(ref, actual, tol) << "at row " << irow; + } + } +} diff --git a/docs/source/cpp/deterministic_metapop.rst b/docs/source/cpp/deterministic_metapop.rst new file mode 100644 index 0000000000..e8d40add30 --- /dev/null +++ b/docs/source/cpp/deterministic_metapop.rst @@ -0,0 +1,57 @@ +Deterministic metapopulation model +================================================= + +Introduction +---------------- + +This metapopulation model incorporates the effects of spatial heterogeneity and population dynamics into a system of ordinary differential equations (ODEs). A population is divided into several subpopulations, each representing spatially seperated regions which interact on some level. Each subpopulation is further divided into epidemiological compartments and eventually into age groups. +Commuting between regions is governed by a commuting matrix, which describes the fraction of individuals that commute from one region to another. Commuting is not performed explicitly, i.e., individuals are not exchanged between subpopulations, but rather their theoretical impact to transmission dynamics in other locations is considered, leading to a large system of ODEs which can be solved using standard numerical integration methods. + +Simulation +---------- + +The simulation runs in discrete time steps using standard numerical integration schemes. Several schemes including adaptive step size methods are available, see here. + +How to: Set up and run a simulation of the deterministic metapopulation model +---------------------------------------------------------------------------- + +To set up a simulation with the deterministic SEIR metapopulation model you need to initialize the model with the desired number of regions and age groups, e.g., 3 regions and 1 age group: +.. code-block:: cpp + + mio::oseirmetapop::Model model(3, 1) + +Set a population with the number of individuals in each region, age group and epidemiological compartment, e.g.: + +.. code-block:: cpp + + model.populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Susceptible}] = 900; + model.populations[{mio::oseirmetapop::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed}] = 100; + model.populations[{mio::oseirmetapop::Region(1), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Susceptible}] = 1000; + model.populations[{mio::oseirmetapop::Region(2), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Susceptible}] = 1000; + +and the epidemiological parameters, e.g.: + +.. code-block:: cpp + + model.parameters.template get>() + .get_cont_freq_mat()[0] + .get_baseline() + .setConstant(2.7); + + model.parameters.set>(3.335); + model.parameters.set>(8.097612257); + model.parameters.set>(0.07333); + +Construct an ``Eigen::MatrixXd`` of size :math:`n_{regions} \times n_{regions}` which describes the fraction of individuals commuting from one region to another. The matrix should satify the sum of each row equal to 1.0, e.g.: +.. code-block:: cpp + Eigen::MatrixXd mobility_data_commuter(3, 3); + mobility_data_commuter << 0.4, 0.3, 0.3, 0.2, 0.7, 0.1, 0.4, 0.1, 0.5; + +Set the commuting strengths matrix via the ``set_commuting_strengths`` method to ensure that the population after commuting is correctly updated: +.. code-block:: cpp + + model.set_commuting_strengths(mobility_data_commuter); + +Finally, to run the simulation from `t0` to `tmax` with a time step of `dt`, use the following command: +.. code-block:: cpp + simulate(t0, tmax, dt, model); \ No newline at end of file diff --git a/pycode/examples/plot/plotResultsMapGermany.py b/pycode/examples/plot/plotResultsMapGermany.py index 68c8940282..1633efe88c 100755 --- a/pycode/examples/plot/plotResultsMapGermany.py +++ b/pycode/examples/plot/plotResultsMapGermany.py @@ -68,13 +68,13 @@ try: population = pd.read_json( - 'data/pydata/Germany/county_current_population.json') + 'data/Germany/pydata/county_current_population.json') # pandas>1.5 raise FileNotFoundError instead of ValueError except (ValueError, FileNotFoundError): print("Population data was not found. Download it from the internet.") population = gpd.get_population_data( read_data=False, file_format=file_format, - out_folder='data/pydata/Germany/', no_raw=True, + out_folder='data/Germany/pydata/', no_raw=True, merge_eisenach=True) # For fitting of different age groups we need format ">X". diff --git a/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py b/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py new file mode 100644 index 0000000000..5232f2a3b5 --- /dev/null +++ b/pycode/memilio-epidata/memilio/epidata/getNRWCounties.py @@ -0,0 +1,56 @@ +import os + +import numpy as np +import pandas as pd + +from memilio.epidata import geoModificationGermany as geoger +from memilio.epidata import getDataIntoPandasDataFrame as gd + + +def main(): + """! Main program entry.""" + + arg_dict = gd.cli("commuter_official") + + directory = os.path.join( + arg_dict['out_folder'].split('/pydata/')[0], 'Germany/') + directory_mobility = os.path.join(directory, 'mobility/') + directory_population = os.path.join(directory, 'pydata/') + mobility_file = 'commuter_mobility_2022' + population_file = 'county_current_population' + + mobility_matrix = pd.read_csv( + os.path.join(directory_mobility + mobility_file + '.txt'), + sep=' ', header=None) + + # get county and state IDs + countyIDs = geoger.get_county_ids() + stateIDs = geoger.get_state_ids() + # get state ID to county ID map + stateID_to_countyID = geoger.get_stateid_to_countyids_map() + + # iterate over state_to_county map and replace IDs by numbering 0, ..., n + state_indices = [] + county_indices = [] + for state, counties in stateID_to_countyID.items(): + state_indices.append(stateIDs.index(state)) + county_indices.append( + np.array([countyIDs.index(county) for county in counties])) + + mobility_matrix_nrw = mobility_matrix.loc[county_indices[4], + county_indices[4]] + + gd.write_dataframe( + mobility_matrix_nrw, directory_mobility, mobility_file + '_nrw', 'txt', + param_dict={'sep': ' ', 'header': None, 'index': False}) + + population = pd.read_json(os.path.join( + directory_population + population_file + '.json')) + population_nrw = population.loc[county_indices[4]] + gd.write_dataframe(population_nrw, directory_population, + population_file + '_nrw', 'json') + + +if __name__ == "__main__": + + main() diff --git a/pycode/memilio-plot/memilio/plot/createGIF.py b/pycode/memilio-plot/memilio/plot/createGIF.py index 9ea23f9700..7806d30953 100644 --- a/pycode/memilio-plot/memilio/plot/createGIF.py +++ b/pycode/memilio-plot/memilio/plot/createGIF.py @@ -71,14 +71,14 @@ def create_plot_map(day, filename, files_input, output_path, compartments, file try: population = pd.read_json( - 'data/pydata/Germany/county_current_population.json') + 'data/Germany/pydata/county_current_population.json') # pandas>1.5 raise FileNotFoundError instead of ValueError except (ValueError, FileNotFoundError): print( "Population data was not found. Downloading it from the internet.") population = gpd.get_population_data( read_data=False, file_format=file_format, - out_folder='data/pydata/Germany/', no_raw=True, merge_eisenach=True) + out_folder='data/Germany/pydata/', no_raw=True, merge_eisenach=True) # For fitting of different age groups we need format ">X". age_group_values = list(age_groups.values())