diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 9b6715b3fa..4882dc262b 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -190,6 +190,7 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/ode_seirdb) 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 e63da90ac1..502d109c85 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -112,6 +112,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..1bfdf88d9c --- /dev/null +++ b/cpp/examples/ode_seir_metapop.cpp @@ -0,0 +1,80 @@ +/* +* Copyright (C) 2020-2026 MEmilio +* +* Authors: Carlotta Gerstein +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "memilio/compartments/simulation.h" +#include "memilio/utils/custom_index_array.h" +#include "models/ode_seir_metapop/model.h" +#include "models/ode_seir_metapop/parameters.h" +#include "memilio/geography/regions.h" +#include "memilio/data/analyze_result.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::Infected}] += 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.get>()[{mio::regions::Region(0), mio::AgeGroup(0)}] = 3.; + model.parameters.get>()[{mio::regions::Region(1), mio::AgeGroup(0)}] = 4.; + model.parameters.get>()[{mio::regions::Region(2), mio::AgeGroup(0)}] = 5.; + model.parameters.get>()[{mio::regions::Region(0), mio::AgeGroup(0)}] = 7.; + model.parameters.get>()[{mio::regions::Region(1), mio::AgeGroup(0)}] = 8.; + model.parameters.get>()[{mio::regions::Region(2), mio::AgeGroup(0)}] = 9.; + model.parameters.set>(0.07333); + + auto result = simulate(t0, tmax, dt, model); + auto interpolated_result = mio::interpolate_simulation_result(result); + + std::vector vars = {"S", "E", "I", "R"}; + printf("Infected individuals per Region over time [%%]:\n"); + for (size_t i = 0; i < (size_t)model.parameters.get_num_regions(); i++) { + printf("\t Region %zd", i); + } + for (ScalarType t : interpolated_result.get_times()) { + printf("\n %f", t); + for (size_t i = 0; i < (size_t)model.parameters.get_num_regions(); i++) { + + printf("\t %.5f ", + interpolated_result.get_value((size_t)t)[(size_t)model.parameters.get_num_regions() * + ((size_t)mio::oseirmetapop::InfectionState::Count - 1) + + i] / + model.populations.get_group_total(mio::regions::Region(i)) * 100); + } + } + return 0; +} diff --git a/cpp/memilio/epidemiology/contact_matrix.h b/cpp/memilio/epidemiology/contact_matrix.h index 47e34e3a8a..e7adf17324 100644 --- a/cpp/memilio/epidemiology/contact_matrix.h +++ b/cpp/memilio/epidemiology/contact_matrix.h @@ -32,7 +32,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..66a792cff9 --- /dev/null +++ b/cpp/models/ode_seir_metapop/CMakeLists.txt @@ -0,0 +1,12 @@ +add_library(ode_seir_metapop + infection_state.h + 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/README.md b/cpp/models/ode_seir_metapop/README.md new file mode 100644 index 0000000000..c7bd1f6974 --- /dev/null +++ b/cpp/models/ode_seir_metapop/README.md @@ -0,0 +1,9 @@ +# ODE-based SEIR model + +This directory contains a metapopulation model implementation based on an ODE formulation, modeling the mobility implicitly by including the impact of other regions into the force of infection. +To get started, check out the [official documentation](https://memilio.readthedocs.io/en/latest/cpp/ode_metapop.html) +or the code example: + +- [examples/ode_seir_metapop.cpp](../../examples/ode_seir_metapop.cpp) + +For the metapopulation model approach that models the mobility explicitly, have a look at the [corresponding documentation])(https://memilio.readthedocs.io/en/latest/cpp/graph_metapop.html). \ No newline at end of file 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..65a570d9ad --- /dev/null +++ b/cpp/models/ode_seir_metapop/infection_state.h @@ -0,0 +1,36 @@ +/* +* Copyright (C) 2020-2026 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..43dff29446 --- /dev/null +++ b/cpp/models/ode_seir_metapop/model.cpp @@ -0,0 +1,28 @@ +/* +* Copyright (C) 2020-2026 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..96906464b2 --- /dev/null +++ b/cpp/models/ode_seir_metapop/model.h @@ -0,0 +1,310 @@ +/* +* Copyright (C) 2020-2026 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 for every region 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 AgeGroup%s. + * @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(SimulationTime(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 region = 0; region < n_regions; 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++) { + const size_t Ejn = population.get_flat_index({Region(region), AgeGroup(age_j), Exposed}); + const size_t Ijn = population.get_flat_index({Region(region), AgeGroup(age_j), Infected}); + const size_t Rjn = population.get_flat_index({Region(region), AgeGroup(age_j), Recovered}); + const size_t Sjn = population.get_flat_index({Region(region), AgeGroup(age_j), Susceptible}); + + const FP Njn = pop[Sjn] + pop[Ejn] + pop[Ijn] + pop[Rjn]; + const FP divNjn = (Njn < Limits::zero_tolerance()) ? FP(0.0) : FP(1.0 / Njn); + FP coeffStoE = + 0.5 * + params.template get>().get_cont_freq_mat().get_matrix_at( + SimulationTime(t))(age_i, age_j) * + params.template get>()[{Region(region), AgeGroup(age_i)}]; + + flows[Base::template get_flat_flow_index( + {Region(region), AgeGroup(age_i)})] += + (pop[Ijn] * divNjn + infections_due_commuting(region, age_j)) * coeffStoE * + y[population.get_flat_index({Region(region), AgeGroup(age_i), Susceptible})]; + } + 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>()[{Region(region), 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>()[{Region(region), 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>(); + auto 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(SimulationTime(y.get_time(t_idx)))(i.get(), j.get()) * + params.template get>()[{n, 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(SimulationTime(y.get_time(t_idx)))(n.get(), + k.get()) * + commuting_strengths.get_matrix_at(SimulationTime(y.get_time(t_idx)))(m.get(), + k.get()) / + population_after_commuting[{k, j}]; + } + } + } + + double T_Ei = params.template get>()[{n, i}]; + double T_Ii = params.template get>()[{n, 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); + FP rho = ces.eigenvalues().cwiseAbs().maxCoeff(); + + return mio::success(rho); + } + + /** + *@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)); + } + + /** + * serialize this. + * @see mio::serialize + */ + template + void serialize(IOContext& io) const + { + auto obj = io.create_object("Model"); + obj.add_element("Parameters", this->parameters); + obj.add_element("Populations", this->populations); + } + + /** + * deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + auto obj = io.expect_object("Model"); + auto par = obj.expect_element("Parameters", Tag{}); + auto pop = obj.expect_element("Populations", Tag{}); + return apply( + io, + [](auto&& par_, auto&& pop_) { + return Model{pop_, par_}; + }, + par, pop); + } +}; + +} // namespace oseirmetapop +} // namespace mio + +#endif // SEIRMETAPOP_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..d3f8368a43 --- /dev/null +++ b/cpp/models/ode_seir_metapop/parameters.h @@ -0,0 +1,345 @@ +/* +* Copyright (C) 2020-2026 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, Region, AgeGroup>; + static Type get_default(Region size_region, AgeGroup size_age) + { + return Type({size_region, size_age}, 1.0); + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/** + * @brief The latent time in day unit. + */ +template +struct TimeExposed { + using Type = CustomIndexArray, Region, AgeGroup>; + static Type get_default(Region size_region, AgeGroup size_age) + { + return Type({size_region, size_age}, 5.2); + } + static std::string name() + { + return "TimeExposed"; + } +}; + +/** + * @brief The infectious time in day unit. + */ +template +struct TimeInfected { + using Type = CustomIndexArray, Region, AgeGroup>; + static Type get_default(Region size_region, AgeGroup size_age) + { + return Type({size_region, size_age}, 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() + { + const FP tol_times = 1e-1; + bool corrected = false; + + for (auto i = AgeGroup(0); i < AgeGroup(m_num_agegroups); i++) { + for (auto j = Region(0); j < Region(m_num_regions); j++) { + if (this->template get>()[{j, i}] < tol_times) { + log_warning( + "Constraint check: Parameter TimeExposed changed from {} to {}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[{j, i}], tol_times); + this->template get>()[{j, i}] = tol_times; + corrected = true; + } + if (this->template get>()[{j, i}] < tol_times) { + log_warning( + "Constraint check: Parameter TimeInfected changed from {} to {}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[{j, i}], tol_times); + this->template get>()[{j, i}] = tol_times; + corrected = true; + } + if (this->template get>()[{j, i}] < 0.0 || + this->template get>()[{j, i}] > 1.0) { + log_warning("Constraint check: Parameter TransmissionProbabilityOnContact changed from {} to {} ", + this->template get>()[{j, i}], 0.0); + this->template get>() = 0.0; + corrected = true; + } + + if (this->template get>()[{j, i}] <= 0.0) { + log_warning( + "Constraint check: Parameter PopulationAfterCommuting changed from {} to {}. 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(SimulationTime(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(SimulationTime(0)) + .minCoeff() < 0.0 || + this->template get>() + .get_cont_freq_mat() + .get_matrix_at(SimulationTime(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++) { + for (auto j = Region(0); j < Region(m_num_regions); j++) { + if (this->template get>()[{j, i}] < tol_times) { + log_error( + "Constraint check: Parameter TimeExposed {} smaller or equal {}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[{j, i}], 0.0); + return true; + } + if (this->template get>()[{j, i}] < tol_times) { + log_error( + "Constraint check: Parameter TimeInfected {} smaller or equal {}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->template get>()[{j, i}], 0.0); + return true; + } + if (this->template get>()[{j, i}] < 0.0 || + this->template get>()[{j, i}] > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact {} smaller {} or " + "greater {:.4f}", + this->template get>()[{j, i}], 0.0, 1.0); + return true; + } + + if (this->template get>()[{j, i}] <= 0.0) { + log_error("Constraint check: Parameter PopulationAfterCommuting {} smaller or equal {}", + this->template get>()[{j, i}], 0.0); + return true; + } + } + } + if ((this->template get>() + .get_cont_freq_mat() + .get_matrix_at(SimulationTime(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(SimulationTime(0)) + .minCoeff() < 0.0 || + this->template get>() + .get_cont_freq_mat() + .get_matrix_at(SimulationTime(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 // SEIRMETAPOP_PARAMETERS_H diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 6447c6eabc..7c2c89cfff 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -90,6 +90,7 @@ set(TESTSOURCES temp_file_register.h test_graph_abm.cpp test_abstract_parameter_dist.cpp + test_odemetapop.cpp test_temporal_hybrid_model.cpp test_odemseirs4.cpp test_uncertain_ad_compatibility.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..1e3f0cfac8 --- /dev/null +++ b/cpp/tests/test_odemetapop.cpp @@ -0,0 +1,481 @@ + +#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::oseirmetapop::InfectionState::Exposed)}] = 10000; + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Susceptible)}] = + total_population_per_region - + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed)}] - + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] - + model.populations[{mio::Index( + mio::regions::Region(i), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}]; + } + model.set_commuting_strengths(); + } +}; + +TEST_F(ModelTestOdeMetapop, simulateDefault) +{ + mio::TimeSeries result = simulate(t0, tmax, dt, model); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); +} + +TEST_F(ModelTestOdeMetapop, checkPopulationConservation) +{ + auto result = simulate(t0, tmax, dt, model); + ScalarType num_persons = result.get_last_value().sum(); + EXPECT_NEAR(num_persons, total_population_per_region * (size_t)model.parameters.get_num_regions(), 1e-9); +} + +TEST_F(ModelTestOdeMetapop, compareWithPreviousRun) +{ + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.1); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + std::vector> refData = load_test_data_csv("ode-seir-metapop-compare.csv"); + 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) +{ + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.04); + model.parameters.get>().get_cont_freq_mat()[0].get_baseline()(0, 0) = + 10; + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + // model.check_constraints() combines the functions from population and parameters. + // We only want to test the functions for the parameters defined in parameters.h + ASSERT_EQ(model.parameters.check_constraints(), 0); + + mio::set_log_level(mio::LogLevel::off); + + model.parameters.set>(-5.2); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.set>(5.2); + model.parameters.set>(0); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.set>(6); + model.parameters.set>(10.); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + model.parameters.set>(0.04); + mobility_data_commuter(0, 1) += 0.5; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 1) = -0.5; + mobility_data_commuter(0, 2) = 0.75; + mobility_data_commuter(0, 3) = 0.75; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 1) = 1.5; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 0.; + model.set_commuting_strengths(mobility_data_commuter); + ASSERT_EQ(model.parameters.check_constraints(), 1); + + mobility_data_commuter(0, 0) = 0.; + mobility_data_commuter(0, 1) = 0.; + mobility_data_commuter(0, 2) = 0.; + mobility_data_commuter(0, 3) = 1.; + model.set_commuting_strengths(mobility_data_commuter); + model.parameters.set>( + mio::Populations({mio::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; + model.parameters.set>(5.2); + model.parameters.set>(2); + model.parameters.set>(0.04); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + mio::set_log_level(mio::LogLevel::off); + model.parameters.set>(-5.2); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ((model.parameters + .get>()[{mio::regions::Region(0), mio::AgeGroup(0)}]), + tol_times); + + model.parameters.set>(1e-5); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_EQ((model.parameters + .get>()[{(mio::regions::Region)0, (mio::AgeGroup)0}]), + tol_times); + + model.parameters.set>(10.); + EXPECT_EQ(model.parameters.apply_constraints(), 1); + EXPECT_NEAR((model.parameters.get>()[{ + (mio::regions::Region)0, (mio::AgeGroup)0}]), + 0.0, 1e-14); + + model.parameters.set>(0.04); + + 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_F(ModelTestOdeMetapop, get_reproduction_numbers) +{ + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.04); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0].get_baseline().setConstant(10); + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + model.apply_constraints(); + + tmax = 1; + + Eigen::VectorXd checkReproductionNumbers(4); + checkReproductionNumbers << 2.3728557964184716, 2.3726126970254491, 2.3710580092591909, 2.3687971820697493; + + Eigen::VectorXd checkReproductionNumbers2(4); + checkReproductionNumbers2 << 2.1355702167766251, 2.1353514273229033, 2.1339522083332727, 2.1319174638627727; + + Eigen::VectorXd checkReproductionNumbers3(4); + checkReproductionNumbers3 << 1.8982846371347806, 1.8980901576203593, 1.8968464074073521, 1.8950377456557974; + + mio::TimeSeries result((size_t)model.parameters.get_num_regions() * + (size_t)mio::oseir::InfectionState::Count); + mio::TimeSeries::Vector result_0(16); + mio::TimeSeries::Vector result_1(16); + mio::TimeSeries::Vector result_2(16); + mio::TimeSeries::Vector result_3(16); + + result_0.setZero(); + result_1.setZero(); + result_2.setZero(); + result_3.setZero(); + + for (auto region = mio::regions::Region(0); region < mio::regions::Region(4); region++) { + result_0[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1049000.00000; + result_1[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1048892.52981; + result_2[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1048205.22826; + result_3[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1047205.75424; + } + + result.add_time_point(0, result_0); + result.add_time_point(0.10000, result_1); + result.add_time_point(0.55000, result_2); + result.add_time_point(1.00000, result_3); + + auto reproduction_numbers = model.get_reproduction_numbers(result); + + for (int i = 0; i < reproduction_numbers.size(); i++) { + EXPECT_NEAR(reproduction_numbers[i], checkReproductionNumbers[i], 1e-12); + } + + contact_matrix[0].get_baseline().setConstant(9); + + auto reproduction_numbers2 = model.get_reproduction_numbers(result); + + for (int i = 0; i < reproduction_numbers2.size(); i++) { + EXPECT_NEAR(reproduction_numbers2[i], checkReproductionNumbers2[i], 1e-12); + } + + contact_matrix[0].get_baseline().setConstant(8); + + auto reproduction_numbers3 = model.get_reproduction_numbers(result); + + for (int i = 0; i < reproduction_numbers2.size(); i++) { + EXPECT_NEAR(reproduction_numbers3[i], checkReproductionNumbers3[i], 1e-12); + } +} + +TEST_F(ModelTestOdeMetapop, get_reproduction_number) +{ + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.04); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0].get_baseline().setConstant(10); + + Eigen::MatrixXd mobility_data_commuter((size_t)model.parameters.get_num_regions(), + (size_t)model.parameters.get_num_regions()); + mobility_data_commuter << 0., 0., 0., 1., 0.2, 0., 0.6, 0.2, 0., 0.5, 0.5, 0., 0., 0., 0., 1.; + model.set_commuting_strengths(mobility_data_commuter); + + model.apply_constraints(); + + tmax = 1; + + mio::TimeSeries result((size_t)model.parameters.get_num_regions() * + (size_t)mio::oseir::InfectionState::Count); + mio::TimeSeries::Vector result_0(16); + mio::TimeSeries::Vector result_1(16); + mio::TimeSeries::Vector result_2(16); + mio::TimeSeries::Vector result_3(16); + + result_0.setZero(); + result_1.setZero(); + result_2.setZero(); + result_3.setZero(); + + for (auto region = mio::regions::Region(0); region < mio::regions::Region(4); region++) { + result_0[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1049000.00000; + result_1[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1048892.52981; + result_2[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1048205.22826; + result_3[model.populations.get_flat_index( + {region, mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible})] = 1047205.75424; + } + + result.add_time_point(0, result_0); + result.add_time_point(0.10000, result_1); + result.add_time_point(0.55000, result_2); + result.add_time_point(1.00000, result_3); + + EXPECT_NEAR(model.get_reproduction_number((size_t)0, result).value(), 2.3728557964184716, 1e-12); + EXPECT_NEAR(model.get_reproduction_number((size_t)1, result).value(), 2.3726126970254491, 1e-12); + EXPECT_NEAR(model.get_reproduction_number((size_t)2, result).value(), 2.3710580092591909, 1e-12); + + contact_matrix[0].get_baseline().setConstant(9); + EXPECT_NEAR(model.get_reproduction_number((size_t)0, result).value(), 2.1355702167766251, 1e-12); + EXPECT_NEAR(model.get_reproduction_number((size_t)1, result).value(), 2.1353514273229033, 1e-12); + + contact_matrix[0].get_baseline().setConstant(8); + EXPECT_NEAR(model.get_reproduction_number((size_t)0, result).value(), 1.8982846371347806, 1e-12); + EXPECT_NEAR(model.get_reproduction_number((size_t)1, result).value(), 1.8980901576203593, 1e-12); +} + +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::oseirmetapop::InfectionState::Exposed)}] = 10000; + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Susceptible)}] = + total_population - + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Exposed)}] - + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Infected)}] - + model.populations[{mio::Index( + mio::regions::Region(0), mio::AgeGroup(0), mio::oseirmetapop::InfectionState::Recovered)}]; + + model.parameters.set>(5.2); + model.parameters.set>(6); + model.parameters.set>(0.1); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get>(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); + + model.set_commuting_strengths(); + + model.check_constraints(); + + // Use the Euler integrator as adaptive methods make different time steps in this model due to restructured equations. + auto result = simulate(t0, tmax, dt, model, std::make_unique>()); + 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/metapop.rst b/docs/source/cpp/metapop.rst index 77d80eb44e..f1e1fedb41 100644 --- a/docs/source/cpp/metapop.rst +++ b/docs/source/cpp/metapop.rst @@ -8,3 +8,4 @@ MEmilio implements metapopulation models through a generalization of aggregated graph_metapop smm + ode_metapop diff --git a/docs/source/cpp/ode_metapop.rst b/docs/source/cpp/ode_metapop.rst new file mode 100644 index 0000000000..1320ac6e80 --- /dev/null +++ b/docs/source/cpp/ode_metapop.rst @@ -0,0 +1,73 @@ +ODE-based metapopulation model +============================== + +Introduction +---------------- + +This model incorporates spatial resolution into an ODE-based metapopulation model by including the effects of other regions into the force of infection term of the 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. + +.. note:: + + In comparison to the :doc:`Graph-based metapopulation model`, 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. For a dense network, this approach scales better with the number of regions than the graph-based approach. Currently, the implicit mobility approach is only available for the SEIR model. + +Simulation +---------- + +As the model is given as a large system of ODEs, the simulation can be performed by solving the system using a standard ODE solver. The model is build on the same setup as the simpler ODE-models, so we refer to :doc:`ode` for more details. + + +How to: Set up and run the ODE-based metapopulation model +--------------------------------------------------------- + +To set up a simulation of the ODE 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); + +Currently, the contact matrix is set globally, so it is the same for every region. The other parameters can be set for every region individually by + +.. code-block:: cpp + + model.parameters.get>()[{mio::regions::Region(0), mio::AgeGroup(0)}] = 3. + +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 satifsy the sum of each row equal to 1.0, e.g.: + +.. code-block:: cpp + + Eigen::MatrixXd commuting_strengths(3, 3); + commuting_strengths << 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(commuting_strengths); + +If no mobility data is provided, you should still run ``model.set_commuting_strengths()`` which sets the matrix to an identity matrix, reflecting no mobility between regions. Without setting this, you will encounter errors during the simulation. + +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 439ae9f5f8..2972713a0a 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-plot/memilio/plot/createGIF.py b/pycode/memilio-plot/memilio/plot/createGIF.py index f257a19e01..06d9d9db1b 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())