diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index e982766dc3..52fb3dc63e 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -18,6 +18,10 @@ add_executable(ode_seir_example ode_seir.cpp) target_link_libraries(ode_seir_example PRIVATE memilio ode_seir) target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_seir_ageres_example ode_seir_ageres.cpp) +target_link_libraries(ode_seir_ageres_example PRIVATE memilio ode_seir) +target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(ode_sir_example ode_sir.cpp) target_link_libraries(ode_sir_example PRIVATE memilio ode_sir) target_compile_options(ode_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) @@ -30,6 +34,10 @@ add_executable(sde_sirs_example sde_sirs.cpp) target_link_libraries(sde_sirs_example PRIVATE memilio sde_sirs) target_compile_options(sde_sirs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(ode_sir_ageres_example ode_sir_ageres.cpp) +target_link_libraries(ode_sir_ageres_example PRIVATE memilio ode_sir) +target_compile_options(ode_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) + add_executable(seir_flows_example ode_seir_flows.cpp) target_link_libraries(seir_flows_example PRIVATE memilio ode_seir) target_compile_options(seir_flows_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/graph.cpp b/cpp/examples/graph.cpp index 49d0c48463..6ecbf8a5d3 100644 --- a/cpp/examples/graph.cpp +++ b/cpp/examples/graph.cpp @@ -31,20 +31,30 @@ int main() const auto tmax = 10.; const auto dt = 0.5; //time step of migration, daily migration every second step - mio::oseir::Model model; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 10000; + mio::oseir::Model model(1); + + // set population + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000; + + // set transition times model.parameters.set(1); - model.parameters.get().get_baseline()(0, 0) = 2.7; model.parameters.set(1); + // set contact matrix + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(2.7); + //two mostly identical groups auto model_group1 = model; auto model_group2 = model; + //some contact restrictions in group 1 - model_group1.parameters.get().add_damping(0.5, mio::SimulationTime(5)); + mio::ContactMatrixGroup& contact_matrix1 = model_group1.parameters.get().get_cont_freq_mat(); + contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5)); + //infection starts in group 1 - model_group1.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 9990; - model_group1.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 10; + model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990; + model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 10; mio::Graph>, mio::MigrationEdge> g; g.add_node(1001, model_group1, t0); diff --git a/cpp/examples/ode_seir.cpp b/cpp/examples/ode_seir.cpp index 0ba61eefa2..ad21bc76fb 100644 --- a/cpp/examples/ode_seir.cpp +++ b/cpp/examples/ode_seir.cpp @@ -23,6 +23,8 @@ #include "memilio/compartments/simulation.h" #include "memilio/utils/logging.h" +#include "memilio/utils/time_series.h" + int main() { mio::set_log_level(mio::LogLevel::debug); @@ -33,26 +35,25 @@ int main() mio::log_info("Simulating ODE SEIR; t={} ... {} with dt = {}.", t0, tmax, dt); - mio::oseir::Model model; + mio::oseir::Model model(1); - double total_population = 1061000; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 10000; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 1000; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 1000; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = - total_population - - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; + double total_population = 10000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; - model.parameters.set(0.04); model.parameters.set(5.2); - model.parameters.set(2); - - model.parameters.get().get_baseline()(0, 0) = 2.7; + model.parameters.set(6); + model.parameters.set(0.1); - // contacts increase by 100% after 12.5 days - model.parameters.get().add_damping(-1., mio::SimulationTime(12.5)); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.)); + model.check_constraints(); diff --git a/cpp/examples/ode_seir_ageres.cpp b/cpp/examples/ode_seir_ageres.cpp new file mode 100644 index 0000000000..03925413e5 --- /dev/null +++ b/cpp/examples/ode_seir_ageres.cpp @@ -0,0 +1,62 @@ +#include "ode_seir/model.h" +#include "ode_seir/infection_state.h" +#include "ode_seir/parameters.h" +#include "memilio/compartments/simulation.h" +#include "memilio/utils/logging.h" + +int main() +{ + mio::set_log_level(mio::LogLevel::debug); + + double t0 = 0; + double tmax = 50; + double dt = 0.001; + + mio::log_info("Simulating SEIR; t={} ... {} with dt = {}.", t0, tmax, dt); + + double cont_freq = 10; + + double nb_total_t0 = 10000, nb_exp_t0 = 100, nb_inf_t0 = 50, nb_rec_t0 = 10; + const size_t num_groups = 3; + + mio::oseir::Model model(num_groups); + double fact = 1.0 / num_groups; + + auto& params = model.parameters; + + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); i++) { + model.populations[{i, mio::oseir::InfectionState::Exposed}] = fact * nb_exp_t0; + model.populations[{i, mio::oseir::InfectionState::Infected}] = fact * nb_inf_t0; + model.populations[{i, mio::oseir::InfectionState::Recovered}] = fact * nb_rec_t0; + model.populations.set_difference_from_group_total({i, mio::oseir::InfectionState::Susceptible}, + fact * nb_total_t0); + + model.parameters.get()[i] = 5.2; + model.parameters.get()[i] = 6; + model.parameters.get()[i] = 0.04; + } + + mio::ContactMatrixGroup& contact_matrix = params.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, fact * cont_freq)); + contact_matrix.add_damping(Eigen::MatrixXd::Constant(num_groups, num_groups, 0.7), mio::SimulationTime(30.)); + + model.apply_constraints(); + + mio::TimeSeries seir = simulate(t0, tmax, dt, model); + + std::vector vars = {"S", "E", "I", "R"}; + printf("Number of time points :%d\n", static_cast(seir.get_num_time_points())); + printf("People in\n"); + + for (size_t k = 0; k < (size_t)mio::oseir::InfectionState::Count; k++) { + double dummy = 0; + + for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) { + printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i, + seir.get_last_value()[k + (size_t)mio::oseir::InfectionState::Count * (int)i]); + dummy += seir.get_last_value()[k + (size_t)mio::oseir::InfectionState::Count * (int)i]; + } + + printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy); + } +} diff --git a/cpp/examples/ode_seir_flows.cpp b/cpp/examples/ode_seir_flows.cpp index 97c71315bf..9e3adf5a28 100644 --- a/cpp/examples/ode_seir_flows.cpp +++ b/cpp/examples/ode_seir_flows.cpp @@ -38,25 +38,24 @@ int main() mio::log_info("Simulating SEIR; t={} ... {} with dt = {}.", t0, tmax, dt); - mio::oseir::Model model; + mio::oseir::Model model(1); + + constexpr double total_population = 10000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations.set_difference_from_group_total( + {mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}, total_population); - double total_population = 10000; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = - total_population - - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; - // suscetible now set with every other update - // params.nb_sus_t0 = params.nb_total_t0 - params.nb_exp_t0 - params.nb_inf_t0 - params.nb_rec_t0; model.parameters.set(5.2); model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); model.check_constraints(); + auto seir = simulate_flows(t0, tmax, dt, model); printf("Compartments: \n"); diff --git a/cpp/examples/ode_sir.cpp b/cpp/examples/ode_sir.cpp index 28097a4679..4b6769c5a5 100644 --- a/cpp/examples/ode_sir.cpp +++ b/cpp/examples/ode_sir.cpp @@ -24,6 +24,8 @@ #include "ode_sir/infection_state.h" #include "ode_sir/model.h" #include "ode_sir/parameters.h" +#include +#include int main() { @@ -37,24 +39,23 @@ int main() mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt); - mio::osir::Model model; + mio::osir::Model model(1); - model.populations[{mio::Index(mio::osir::InfectionState::Infected)}] = 1000; - model.populations[{mio::Index(mio::osir::InfectionState::Recovered)}] = 1000; - model.populations[{mio::Index(mio::osir::InfectionState::Susceptible)}] = - total_population - - model.populations[{mio::Index(mio::osir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::osir::InfectionState::Recovered)}]; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] = 1000; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}] = 1000; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}]; model.parameters.set(2); - model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 2.7; - model.parameters.get().add_damping(0.6, mio::SimulationTime(12.5)); - - auto integrator = std::make_shared(); + model.parameters.set(0.5); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); model.check_constraints(); - auto sir = simulate(t0, tmax, dt, model, integrator); + auto integrator = std::make_shared(); + auto sir = simulate(t0, tmax, dt, model, integrator); bool print_to_terminal = true; diff --git a/cpp/examples/ode_sir_ageres.cpp b/cpp/examples/ode_sir_ageres.cpp new file mode 100644 index 0000000000..60e6b88e24 --- /dev/null +++ b/cpp/examples/ode_sir_ageres.cpp @@ -0,0 +1,61 @@ +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/time_series.h" +#include "memilio/utils/logging.h" +#include "memilio/compartments/simulation.h" +#include "ode_sir/model.h" + +int main() +{ + mio::set_log_level(mio::LogLevel::debug); + + double t0 = 0; + double tmax = 50; + double dt = 0.1; + + mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt); + + double cont_freq = 10; // see Polymod study + + double nb_total_t0 = 10000, nb_inf_t0 = 50, nb_rec_t0 = 10; + + const size_t num_groups = 3; + mio::osir::Model model(num_groups); + + double fact = 1.0 / num_groups; + + auto& params = model.parameters; + + for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); i++) { + model.populations[{i, mio::osir::InfectionState::Infected}] = fact * nb_inf_t0; + model.populations[{i, mio::osir::InfectionState::Recovered}] = fact * nb_rec_t0; + model.populations.set_difference_from_group_total({i, mio::osir::InfectionState::Susceptible}, + fact * nb_total_t0); + + model.parameters.get()[i] = 2.0; + model.parameters.get()[i] = 0.3; + } + + mio::ContactMatrixGroup& contact_matrix = params.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, fact * cont_freq)); + contact_matrix.add_damping(Eigen::MatrixXd::Constant(num_groups, num_groups, 0.7), mio::SimulationTime(30.)); + + model.apply_constraints(); + + mio::TimeSeries sir = simulate(t0, tmax, dt, model); + + std::vector vars = {"S", "I", "R"}; + printf("Number of time points :%d\n", static_cast(sir.get_num_time_points())); + printf("People in\n"); + + for (size_t k = 0; k < (size_t)mio::osir::InfectionState::Count; k++) { + double dummy = 0; + + for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) { + printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i, + sir.get_last_value()[k + (size_t)mio::osir::InfectionState::Count * (int)i]); + dummy += sir.get_last_value()[k + (size_t)mio::osir::InfectionState::Count * (int)i]; + } + + printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy); + } +} \ No newline at end of file diff --git a/cpp/models/ode_seir/model.h b/cpp/models/ode_seir/model.h index 95e5459ca5..a0d68d4165 100644 --- a/cpp/models/ode_seir/model.h +++ b/cpp/models/ode_seir/model.h @@ -20,13 +20,10 @@ #ifndef SEIR_MODEL_H #define SEIR_MODEL_H -#include "memilio/compartments/flow_model.h" #include "memilio/epidemiology/populations.h" #include "memilio/epidemiology/contact_matrix.h" #include "memilio/utils/type_list.h" #include "memilio/compartments/compartmentalmodel.h" -#include "memilio/epidemiology/populations.h" -#include "memilio/epidemiology/contact_matrix.h" #include "memilio/io/io.h" #include "memilio/math/interpolation.h" #include "memilio/utils/time_series.h" @@ -34,6 +31,8 @@ #include "ode_seir/parameters.h" #include #include +#include "memilio/epidemiology/age_group.h" +#include "memilio/compartments/flow_simulation.h" namespace mio { @@ -50,29 +49,51 @@ using Flows = TypeList>; // clang-format on -class Model : public FlowModel, Parameters, Flows> +class Model : public FlowModel, Parameters, Flows> { - using Base = FlowModel, Parameters, Flows>; + using Base = FlowModel, Parameters, Flows>; public: - Model() - : Base(Populations({InfectionState::Count}, 0.), ParameterSet()) + Model(int num_agegroups) + : Base(Populations({AgeGroup(num_agegroups), InfectionState::Count}), ParameterSet(AgeGroup(num_agegroups))) { } void get_flows(Eigen::Ref pop, Eigen::Ref y, double t, Eigen::Ref flows) const override { - auto& params = this->parameters; - double coeffStoE = params.get().get_matrix_at(t)(0, 0) * - params.get() / populations.get_total(); - - flows[get_flat_flow_index()] = - coeffStoE * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected]; - flows[get_flat_flow_index()] = - (1.0 / params.get()) * y[(size_t)InfectionState::Exposed]; - flows[get_flat_flow_index()] = - (1.0 / params.get()) * y[(size_t)InfectionState::Infected]; + auto& params = this->parameters; + AgeGroup n_agegroups = params.get_num_groups(); + ContactMatrixGroup const& contact_matrix = params.get(); + + for (auto i = AgeGroup(0); i < n_agegroups; i++) { + size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible}); + size_t Ei = this->populations.get_flat_index({i, InfectionState::Exposed}); + size_t Ii = this->populations.get_flat_index({i, InfectionState::Infected}); + + for (auto j = AgeGroup(0); j < n_agegroups; j++) { + + size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible}); + size_t Ej = this->populations.get_flat_index({j, InfectionState::Exposed}); + size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected}); + size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered}); + + double Nj = pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj]; + double divNj = 1.0 / Nj; + + double coeffStoE = contact_matrix.get_matrix_at(t)(static_cast((size_t)i), + static_cast((size_t)j)) * + params.get()[i] * divNj; + + double dummy_S = y[Si] * y[Ij] * coeffStoE; + + flows[get_flat_flow_index({i})] += dummy_S; + } + flows[get_flat_flow_index({i})] = + (1.0 / params.get()[i]) * y[Ei]; + flows[get_flat_flow_index({i})] = + (1.0 / params.get()[i]) * y[Ii]; + } } /** @@ -87,18 +108,55 @@ class Model : public FlowModel, Para return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries"); } - ScalarType TimeInfected = this->parameters.get(); + auto const& params = this->parameters; + + const size_t num_groups = (size_t)params.get_num_groups(); + constexpr size_t num_infected_compartments = 2; + const size_t total_infected_compartments = num_infected_compartments * num_groups; + + ContactMatrixGroup const& contact_matrix = params.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_groups); i++) { + size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible}); + for (auto j = AgeGroup(0); j < AgeGroup(num_groups); j++) { + + double Nj = this->populations.get_group_total(j); + double divNj = 1.0 / Nj; + + double coeffStoE = contact_matrix.get_matrix_at(y.get_time(t_idx))( + static_cast((size_t)i), static_cast((size_t)j)) * + params.get()[i] * divNj; + F((size_t)i, (size_t)j + num_groups) = coeffStoE * y.get_value(t_idx)[Si]; + } + + double T_Ei = params.get()[i]; + double T_Ii = params.get()[i]; + V((size_t)i, (size_t)i) = 1.0 / T_Ei; + V((size_t)i + num_groups, (size_t)i) = -1.0 / T_Ei; + V((size_t)i + num_groups, (size_t)i + num_groups) = 1.0 / T_Ii; + } + + V = V.inverse(); + + Eigen::MatrixXd NextGenMatrix = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments); + NextGenMatrix = F * V; - ScalarType coeffStoE = this->parameters.get().get_matrix_at( - y.get_time(static_cast(t_idx)))(0, 0) * - this->parameters.get() / - this->populations.get_total(); + //Compute the largest eigenvalue in absolute value + Eigen::ComplexEigenSolver ces; - ScalarType result = - y.get_value(static_cast(t_idx))[(Eigen::Index)mio::oseir::InfectionState::Susceptible] * - TimeInfected * coeffStoE; + ces.compute(NextGenMatrix); + const Eigen::VectorXcd eigen_vals = ces.eigenvalues(); - return mio::success(result); + Eigen::VectorXd eigen_vals_abs; + eigen_vals_abs.resize(eigen_vals.size()); + + for (int i = 0; i < eigen_vals.size(); i++) { + eigen_vals_abs[i] = std::abs(eigen_vals[i]); + } + return mio::success(eigen_vals_abs.maxCoeff()); } /** diff --git a/cpp/models/ode_seir/parameters.h b/cpp/models/ode_seir/parameters.h index 1fa6d025f5..b4070f76ac 100644 --- a/cpp/models/ode_seir/parameters.h +++ b/cpp/models/ode_seir/parameters.h @@ -23,6 +23,9 @@ #include "memilio/utils/uncertain_value.h" #include "memilio/epidemiology/contact_matrix.h" #include "memilio/utils/parameter_set.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/custom_index_array.h" +#include "memilio/epidemiology/uncertain_matrix.h" #include @@ -39,10 +42,10 @@ namespace oseir * @brief probability of getting infected from a contact */ struct TransmissionProbabilityOnContact { - using Type = UncertainValue; - static Type get_default() + using Type = CustomIndexArray; + static Type get_default(AgeGroup size) { - return Type(1.0); + return Type(size, 1.); } static std::string name() { @@ -54,10 +57,10 @@ struct TransmissionProbabilityOnContact { * @brief the latent time in day unit */ struct TimeExposed { - using Type = UncertainValue; - static Type get_default() + using Type = CustomIndexArray; + static Type get_default(AgeGroup size) { - return Type(5.2); + return Type(size,5.2); } static std::string name() { @@ -69,10 +72,10 @@ struct TimeExposed { * @brief the infectious time in day unit */ struct TimeInfected { - using Type = UncertainValue; - static Type get_default() + using Type = CustomIndexArray; + static Type get_default(AgeGroup size) { - return Type(6.0); + return Type(size,6.0); } static std::string name() { @@ -84,10 +87,10 @@ struct TimeInfected { * @brief the contact patterns within the society are modelled using a ContactMatrix */ struct ContactPatterns { - using Type = ContactMatrix; - static Type get_default() + using Type = UncertainContactMatrix; + static Type get_default(AgeGroup size) { - return Type{1}; + return Type(1, static_cast((size_t)size)); } static std::string name() { @@ -103,11 +106,17 @@ using ParametersBase = ParameterSetget() < tol_times) { - log_warning("Constraint check: Parameter TimeExposed changed from {:.4f} to {:.4f}. Please note that " - "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " - "and reset parameters.", - this->get(), tol_times); - this->get() = tol_times; - corrected = true; - } - if (this->get() < tol_times) { - log_warning("Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " - "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " - "and reset parameters.", - this->get(), tol_times); - this->get() = tol_times; - corrected = true; - } - if (this->get() < 0.0 || - this->get() > 1.0) { - log_warning("Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", - this->get(), 0.0); - this->get() = 0.0; - corrected = true; + + for(auto i =AgeGroup(0); i < AgeGroup(m_num_groups);++i){ + if (this->get()[i] < tol_times) { + log_warning("Constraint check: Parameter TimeExposed changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->get()[i], tol_times); + this->get()[i] = tol_times; + corrected = true; + } + if (this->get()[i] < tol_times) { + log_warning("Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->get()[i], tol_times); + this->get()[i] = tol_times; + corrected = true; + } + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_warning("Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", + this->get()[i], 0.0); + this->get()[i] = 0.0; + corrected = true; + } } return corrected; } @@ -161,33 +173,36 @@ class Parameters : public ParametersBase { const double tol_times = 1e-1; - if (this->get() < tol_times) { + for(auto i = AgeGroup(0); i < m_num_groups; i++){ + if (this->get()[i] < tol_times) { log_error("Constraint check: Parameter TimeExposed {:.4f} smaller or equal {:.4f}. Please note that " "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " "and reset parameters.", - this->get(), 0.0); + this->get()[i], 0.0); return true; } - if (this->get() < tol_times) { + if (this->get()[i] < tol_times) { log_error("Constraint check: Parameter TimeInfected {:.4f} smaller or equal {:.4f}. Please note that " "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " "and reset parameters.", - this->get(), 0.0); + this->get()[i], 0.0); return true; } - if (this->get() < 0.0 || - this->get() > 1.0) { + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { log_error( "Constraint check: Parameter TransmissionProbabilityOnContact {:.4f} smaller {:.4f} or greater {:.4f}", - this->get(), 0.0, 1.0); + this->get()[i], 0.0, 1.0); return true; } - return false; + } + return false; } private: Parameters(ParametersBase&& base) : ParametersBase(std::move(base)) + , m_num_groups(get().get_cont_freq_mat().get_num_groups()) { } @@ -202,8 +217,9 @@ class Parameters : public ParametersBase BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); return success(Parameters(std::move(base))); } +private: + AgeGroup m_num_groups; }; - } // namespace oseir } // namespace mio diff --git a/cpp/models/ode_sir/model.h b/cpp/models/ode_sir/model.h index eafe9609cd..f3cf456b48 100644 --- a/cpp/models/ode_sir/model.h +++ b/cpp/models/ode_sir/model.h @@ -22,6 +22,7 @@ #define ODESIR_MODEL_H #include "memilio/compartments/compartmentalmodel.h" +#include "memilio/epidemiology/age_group.h" #include "memilio/epidemiology/populations.h" #include "memilio/epidemiology/contact_matrix.h" #include "ode_sir/infection_state.h" @@ -36,30 +37,48 @@ namespace osir * define the model * ********************/ -class Model : public CompartmentalModel, Parameters> +class Model : public CompartmentalModel, Parameters> { - using Base = CompartmentalModel, Parameters>; + using Base = CompartmentalModel, Parameters>; public: - Model() - : Base(Populations({InfectionState::Count}, 0.), ParameterSet()) + Model(int num_agegroups) + : Base(Populations({AgeGroup(num_agegroups), InfectionState::Count}), ParameterSet(AgeGroup(num_agegroups))) { } void get_derivatives(Eigen::Ref pop, Eigen::Ref y, double t, Eigen::Ref dydt) const override { - auto& params = this->parameters; - double coeffStoI = params.get().get_matrix_at(t)(0, 0) * - params.get() / populations.get_total(); - - dydt[(size_t)InfectionState::Susceptible] = - -coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected]; - dydt[(size_t)InfectionState::Infected] = - coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected] - - (1.0 / params.get()) * y[(size_t)InfectionState::Infected]; - dydt[(size_t)InfectionState::Recovered] = - (1.0 / params.get()) * y[(size_t)InfectionState::Infected]; + auto params = this->parameters; + AgeGroup n_agegroups = params.get_num_groups(); + ContactMatrixGroup const& contact_matrix = params.get(); + + for (auto i = AgeGroup(0); i < n_agegroups; i++) { + + size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible}); + size_t Ii = this->populations.get_flat_index({i, InfectionState::Infected}); + size_t Ri = this->populations.get_flat_index({i, InfectionState::Recovered}); + + for (auto j = AgeGroup(0); j < n_agegroups; j++) { + + size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible}); + size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected}); + size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered}); + + double Nj = pop[Sj] + pop[Ij] + pop[Rj]; + + double coeffStoI = contact_matrix.get_matrix_at(t)(static_cast((size_t)i), + static_cast((size_t)j)) * + params.get()[i] / Nj; + + dydt[Si] += -coeffStoI * y[Si] * pop[Ij]; + dydt[Ii] += coeffStoI * y[Si] * pop[Ij]; + } + dydt[Ii] -= (1.0 / params.get()[i]) * y[Ii]; + dydt[Ri] = (1.0 / params.get()[i]) * y[Ii]; + + } } }; diff --git a/cpp/models/ode_sir/parameters.h b/cpp/models/ode_sir/parameters.h index 229d1efd7e..8125c42842 100644 --- a/cpp/models/ode_sir/parameters.h +++ b/cpp/models/ode_sir/parameters.h @@ -21,9 +21,12 @@ #ifndef SIR_PARAMETERS_H #define SIR_PARAMETERS_H +#include "memilio/epidemiology/age_group.h" +#include "memilio/epidemiology/uncertain_matrix.h" #include "memilio/utils/uncertain_value.h" #include "memilio/epidemiology/contact_matrix.h" #include "memilio/utils/parameter_set.h" +#include "memilio/utils/custom_index_array.h" #include @@ -40,10 +43,10 @@ namespace osir * @brief probability of getting infected from a contact */ struct TransmissionProbabilityOnContact { - using Type = UncertainValue; - static Type get_default() + using Type = CustomIndexArray; + static Type get_default(AgeGroup size) { - return Type(1.0); + return Type(size,1.0); } static std::string name() { @@ -55,10 +58,10 @@ struct TransmissionProbabilityOnContact { * @brief the infectious time in day unit */ struct TimeInfected { - using Type = UncertainValue; - static Type get_default() + using Type = CustomIndexArray; + static Type get_default(AgeGroup size) { - return Type(6.0); + return Type(size,6.0); } static std::string name() { @@ -70,10 +73,10 @@ struct TimeInfected { * @brief the contact patterns within the society are modelled using a ContactMatrix */ struct ContactPatterns { - using Type = ContactMatrix; - static Type get_default() + using Type = UncertainContactMatrix; + static Type get_default(AgeGroup size) { - return Type{1}; + return Type(1, static_cast((size_t)size)); } static std::string name() { @@ -89,11 +92,17 @@ using ParametersBase = ParameterSetget() < tol_times) { - log_warning("Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " - "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " - "and reset parameters.", - this->get(), tol_times); - this->get() = tol_times; - corrected = true; - } - if (this->get() < 0.0 || - this->get() > 1.0) { - log_warning("Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", - this->get(), 0.0); - this->get() = 0.0; - corrected = true; + + for(auto i = AgeGroup(0); i < AgeGroup(m_num_groups);i++){ + if (this->get()[i] < tol_times) { + log_warning("Constraint check: Parameter TimeInfected changed from {:.4f} to {:.4f}. Please note that " + "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " + "and reset parameters.", + this->get()[i], tol_times); + this->get()[i] = tol_times; + corrected = true; + } + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_warning("Constraint check: Parameter TransmissionProbabilityOnContact changed from {:0.4f} to {:d} ", + this->get()[i], 0.0); + this->get() = 0.0; + corrected = true; + } } return corrected; } @@ -139,26 +151,30 @@ class Parameters : public ParametersBase { double tol_times = 1e-1; - if (this->get() < tol_times) { + for(auto i = AgeGroup(0); i < AgeGroup(m_num_groups);i++){ + + if (this->get()[i] < tol_times) { log_error("Constraint check: Parameter TimeInfected {:.4f} smaller or equal {:.4f}. Please note that " "unreasonably small compartment stays lead to massively increased run time. Consider to cancel " "and reset parameters.", - this->get(), 0.0); + this->get()[i], 0.0); return true; } - if (this->get() < 0.0 || - this->get() > 1.0) { + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { log_error( "Constraint check: Parameter TransmissionProbabilityOnContact {:.4f} smaller {:.4f} or greater {:.4f}", - this->get(), 0.0, 1.0); + this->get()[i], 0.0, 1.0); return true; } + } return false; } private: Parameters(ParametersBase&& base) : ParametersBase(std::move(base)) + , m_num_groups(get().get_cont_freq_mat().get_num_groups()) { } @@ -173,6 +189,8 @@ class Parameters : public ParametersBase BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); return success(Parameters(std::move(base))); } +private: + AgeGroup m_num_groups; }; } // namespace osir diff --git a/cpp/tests/sir_example_pre_agegroups.txt b/cpp/tests/sir_example_pre_agegroups.txt new file mode 100644 index 0000000000..53b1bc2c2b --- /dev/null +++ b/cpp/tests/sir_example_pre_agegroups.txt @@ -0,0 +1,501 @@ +Time S E I +0.00000 1059000.00000 1000.00000 1000.00000 +0.10020 1058729.96889 1219.93091 1050.10020 +0.20040 1058400.63359 1488.14742 1111.21898 +0.30060 1057999.01493 1815.20960 1185.77547 +0.40080 1057509.31518 2213.96699 1276.71783 +0.50100 1056912.31676 2700.04522 1387.63802 +0.60120 1056184.65788 3292.43129 1522.91083 +0.70140 1055297.96256 4014.17514 1687.86230 +0.80160 1054217.79966 4893.22706 1888.97328 +0.90180 1052902.44293 5963.43214 2134.12493 +1.00200 1051301.40267 7265.70326 2432.89408 +1.10220 1049353.69964 8849.39310 2796.90727 +1.20240 1046985.85489 10773.88148 3240.26363 +1.30261 1044109.57697 13110.38578 3780.03725 +1.40281 1040619.14330 15943.98649 4436.87021 +1.50301 1036388.49846 19375.83441 5235.66713 +1.60321 1031268.13444 23525.46525 6206.40031 +1.70341 1025081.88171 28533.08745 7385.03084 +1.80361 1017623.83388 34561.62188 8814.54424 +1.90381 1008655.75905 41798.15253 10546.08842 +2.00401 997905.52554 50454.29023 12640.18424 +2.10421 985067.28917 60764.75654 15167.95429 +2.20441 969804.44398 72983.27525 18212.28077 +2.30461 951756.59875 87374.64376 21868.75749 +2.40481 930552.04420 104201.71115 26246.24465 +2.50501 905827.21175 123706.01699 31466.77126 +2.60521 877254.33117 146081.20133 37664.46750 +2.70541 844577.66873 171439.16631 44983.16496 +2.80561 807657.17678 199770.52167 53572.30155 +2.90581 766516.03890 230903.11638 63580.84472 +3.00601 721385.67036 264465.19252 75149.13712 +3.10621 672738.91591 299862.18782 88398.89627 +3.20641 621300.67652 336277.27151 103422.05197 +3.30661 568026.44663 372703.94271 120269.61066 +3.40681 514044.27118 408013.57594 138942.15288 +3.50701 460564.07691 441052.20828 159383.71480 +3.60721 408767.87673 470751.60445 181480.51883 +3.70741 359701.22681 496233.50464 205065.26855 +3.80762 314187.13915 516886.19427 229926.66658 +3.90782 272777.52581 532399.70569 255822.76849 +4.00802 235746.62886 542757.27070 282496.10044 +4.10822 203120.24211 548191.40942 309688.34847 +4.20842 174727.76593 549119.38612 337152.84794 +4.30862 150262.69162 546073.46915 364663.83923 +4.40882 129339.87831 539637.89222 392022.22947 +4.50902 111542.63555 530399.16843 419058.19601 +4.60922 96457.07073 518911.62862 445631.30064 +4.70942 83694.28991 505676.83286 471628.87723 +4.80962 72902.66790 491133.94421 496963.38789 +4.90982 63772.86900 475657.83408 521569.29692 +5.01002 56038.08040 459562.06987 545399.84973 +5.11022 49471.40928 443104.58919 568424.00153 +5.21042 43881.83967 426494.53008 590623.63025 +5.31062 39109.66720 409899.24113 611991.09167 +5.41082 35021.96665 393450.90755 632527.12580 +5.51102 31508.39414 377252.51075 652239.09511 +5.61122 28477.46176 361383.01674 671139.52150 +5.71142 25853.32264 345901.79429 689244.88306 +5.81162 23573.04872 330852.31900 706574.63228 +5.91182 21584.35558 316265.24466 723150.39976 +6.01202 19843.71785 302160.93025 738995.35190 +6.11222 18314.81746 288551.50748 754133.67506 +6.21242 16967.27089 275442.56570 768590.16341 +6.31263 15775.58750 262834.52135 782389.89115 +6.41283 14718.31780 250723.72886 795557.95334 +6.51303 13777.35717 239103.38043 808119.26240 +6.61323 12937.37615 227964.23417 820098.38968 +6.71343 12185.35378 217295.20273 831519.44349 +6.81363 11510.19469 207083.82862 842405.97670 +6.91383 10902.41428 197316.66771 852780.91801 +7.01403 10353.87935 187979.59805 862666.52261 +7.11423 9857.59377 179058.06809 872084.33814 +7.21443 9407.52109 170537.29568 881055.18323 +7.31463 8998.43717 162402.42691 889599.13592 +7.41483 8625.80749 154638.66246 897735.53006 +7.51503 8285.68474 147231.35723 905482.95804 +7.61523 7974.62302 140166.09844 912859.27854 +7.71543 7689.60588 133428.76596 919881.62816 +7.81563 7427.98561 127005.57832 926566.43607 +7.91583 7187.43210 120883.12690 932929.44100 +8.01603 6965.88949 115048.40063 938985.70988 +8.11623 6761.53932 109488.80287 944749.65781 +8.21643 6572.76925 104192.16197 950235.06877 +8.31663 6398.14632 99146.73671 955455.11697 +8.41683 6236.39403 94341.21762 960422.38835 +8.51703 6086.37275 89764.72499 965148.90226 +8.61723 5947.06282 85406.80421 969646.13297 +8.71743 5817.54998 81257.41905 973925.03097 +8.81764 5697.01281 77306.94324 977996.04395 +8.91784 5584.71183 73546.15087 981869.13730 +9.01804 5479.98003 69966.20577 985553.81420 +9.11824 5382.21464 66558.65023 989059.13513 +9.21844 5290.86994 63315.39321 992393.73684 +9.31864 5205.45099 60228.69828 995565.85073 +9.41884 5125.50813 57291.17128 998583.32059 +9.51904 5050.63216 54495.74809 1001453.61975 +9.61924 4980.45009 51835.68226 1004183.86765 +9.71944 4914.62139 49304.53289 1006780.84572 +9.81964 4852.83473 46896.15258 1009251.01270 +9.91984 4794.80500 44604.67567 1011600.51934 +10.02004 4740.27077 42424.50670 1013835.22253 +10.12024 4688.99197 40350.30921 1015960.69882 +10.22044 4640.74788 38376.99473 1017982.25739 +10.32064 4595.33523 36499.71225 1019904.95252 +10.42084 4552.56669 34713.83789 1021733.59542 +10.52104 4512.26931 33014.96504 1023472.76565 +10.62124 4474.28330 31398.89468 1025126.82202 +10.72144 4438.46082 29861.62625 1026699.91293 +10.82164 4404.66494 28399.34866 1028195.98639 +10.92184 4372.76873 27008.43182 1029618.79945 +11.02204 4342.65437 25685.41833 1030971.92730 +11.12224 4314.21239 24427.01570 1032258.77191 +11.22244 4287.34102 23230.08869 1033482.57029 +11.32265 4261.94552 22091.65209 1034646.40239 +11.42285 4237.93764 21008.86378 1035753.19858 +11.52305 4215.23507 19979.01806 1036805.74687 +11.62325 4193.77791 18999.52241 1037806.69968 +11.72345 4173.91181 18067.50863 1038758.57956 +11.82365 4156.04662 17180.18802 1039663.76536 +11.92385 4140.42846 16335.07532 1040524.49622 +12.02405 4127.13038 15529.98285 1041342.88677 +12.12425 4116.06021 14762.99777 1042120.94202 +12.22445 4106.98283 14032.44600 1042860.57117 +12.32465 4099.55354 13336.84694 1043563.59953 +12.42485 4093.35823 12674.86354 1044231.77823 +12.52505 4087.95636 12045.25220 1044866.79143 +12.62525 4082.93409 11446.80492 1045470.26098 +12.72545 4078.16721 10878.08458 1046043.74820 +12.82565 4073.64245 10337.61512 1046588.74242 +12.92585 4069.34728 9823.99371 1047106.65901 +13.02605 4065.26981 9335.88713 1047598.84306 +13.12625 4061.39881 8872.02831 1048066.57288 +13.22645 4057.72365 8431.21307 1048511.06328 +13.32665 4054.23425 8012.29701 1048933.46874 +13.42685 4050.92108 7614.19249 1049334.88643 +13.52705 4047.77511 7235.86590 1049716.35900 +13.62725 4044.78777 6876.33491 1050078.87733 +13.72745 4041.95095 6534.66596 1050423.38309 +13.82766 4039.25699 6209.97185 1050750.77116 +13.92786 4036.69858 5901.40942 1051061.89199 +14.02806 4034.26884 5608.17737 1051357.55379 +14.12826 4031.96122 5329.51418 1051638.52460 +14.22846 4029.76952 5064.69615 1051905.53433 +14.32866 4027.68785 4813.03553 1052159.27662 +14.42886 4025.71064 4573.87869 1052400.41066 +14.52906 4023.83260 4346.60450 1052629.56290 +14.62926 4022.04871 4130.62263 1052847.32866 +14.72946 4020.35422 3925.37210 1053054.27368 +14.82966 4018.74460 3730.31979 1053250.93561 +14.92986 4017.21557 3544.95905 1053437.82538 +15.03006 4015.76308 3368.80838 1053615.42854 +15.13026 4014.38326 3201.41023 1053784.20651 +15.23046 4013.07246 3042.32974 1053944.59781 +15.33066 4011.82719 2891.15367 1054097.01914 +15.43086 4010.64418 2747.48931 1054241.86651 +15.53106 4009.52027 2610.96345 1054379.51628 +15.63126 4008.45252 2481.22141 1054510.32607 +15.73146 4007.43810 2357.92614 1054634.63576 +15.83166 4006.47432 2240.75734 1054752.76833 +15.93186 4005.55866 2129.41061 1054865.03073 +16.03206 4004.68870 2023.59668 1054971.71462 +16.13226 4003.86215 1923.04063 1055073.09722 +16.23246 4003.07683 1827.48123 1055169.44194 +16.33267 4002.33068 1736.67020 1055260.99912 +16.43287 4001.62174 1650.37161 1055348.00664 +16.53307 4000.94815 1568.36125 1055430.69059 +16.63327 4000.30814 1490.42605 1055509.26581 +16.73347 3999.70003 1416.36352 1055583.93645 +16.83367 3999.12223 1345.98122 1055654.89655 +16.93387 3998.57322 1279.09631 1055722.33048 +17.03407 3998.05156 1215.53498 1055786.41346 +17.13427 3997.55589 1155.13211 1055847.31200 +17.23447 3997.08491 1097.73074 1055905.18435 +17.33467 3996.63738 1043.18173 1055960.18088 +17.43487 3996.21214 991.34336 1056012.44450 +17.53507 3995.80808 942.08092 1056062.11100 +17.63527 3995.42413 895.26643 1056109.30944 +17.73547 3995.05930 850.77823 1056154.16247 +17.83567 3994.71263 808.50074 1056196.78663 +17.93587 3994.38321 768.32411 1056237.29268 +18.03607 3994.07019 730.14394 1056275.78587 +18.13627 3993.77275 693.86102 1056312.36623 +18.23647 3993.49011 659.38108 1056347.12880 +18.33667 3993.22154 626.61454 1056380.16393 +18.43687 3992.96632 595.47623 1056411.55744 +18.53707 3992.72381 565.88527 1056441.39092 +18.63727 3992.49336 537.76475 1056469.74189 +18.73747 3992.27438 511.04162 1056496.68401 +18.83768 3992.06628 485.64642 1056522.28729 +18.93788 3991.86854 461.51318 1056546.61828 +19.03808 3991.68064 438.57918 1056569.74018 +19.13828 3991.50208 416.78483 1056591.71309 +19.23848 3991.33240 396.07351 1056612.59409 +19.33868 3991.17116 376.39139 1056632.43745 +19.43888 3991.01794 357.68732 1056651.29474 +19.53908 3990.87234 339.91272 1056669.21494 +19.63928 3990.73398 323.02138 1056686.24464 +19.73948 3990.60250 306.96943 1056702.42807 +19.83968 3990.47756 291.71514 1056717.80730 +19.93988 3990.35883 277.21888 1056732.42229 +20.04008 3990.24600 263.44299 1056746.31101 +20.14028 3990.13878 250.35166 1056759.50956 +20.24048 3990.03689 237.91088 1056772.05223 +20.34068 3989.94007 226.08832 1056783.97161 +20.44088 3989.84807 214.85325 1056795.29868 +20.54108 3989.76063 204.17650 1056806.06287 +20.64128 3989.67755 194.03030 1056816.29215 +20.74148 3989.59859 184.38830 1056826.01311 +20.84168 3989.52356 175.22544 1056835.25100 +20.94188 3989.45226 166.51791 1056844.02983 +21.04208 3989.38450 158.24309 1056852.37241 +21.14228 3989.32011 150.37946 1056860.30042 +21.24248 3989.25893 142.90661 1056867.83446 +21.34269 3989.20078 135.80511 1056874.99411 +21.44289 3989.14552 129.05650 1056881.79798 +21.54309 3989.09301 122.64325 1056888.26373 +21.64329 3989.04311 116.54870 1056894.40818 +21.74349 3988.99570 110.75701 1056900.24730 +21.84369 3988.95063 105.25312 1056905.79625 +21.94389 3988.90781 100.02274 1056911.06945 +22.04409 3988.86712 95.05228 1056916.08061 +22.14429 3988.82845 90.32881 1056920.84275 +22.24449 3988.79170 85.84007 1056925.36824 +22.34469 3988.75677 81.57439 1056929.66884 +22.44489 3988.72359 77.52068 1056933.75573 +22.54509 3988.69205 73.66842 1056937.63954 +22.64529 3988.66208 70.00758 1056941.33034 +22.74549 3988.63360 66.52867 1056944.83773 +22.84569 3988.60653 63.22264 1056948.17083 +22.94589 3988.58081 60.08089 1056951.33830 +23.04609 3988.55637 57.09527 1056954.34836 +23.14629 3988.53314 54.25801 1056957.20885 +23.24649 3988.51107 51.56174 1056959.92718 +23.34669 3988.49010 48.99947 1056962.51044 +23.44689 3988.47016 46.56452 1056964.96532 +23.54709 3988.45122 44.25057 1056967.29821 +23.64729 3988.43322 42.05161 1056969.51518 +23.74749 3988.41611 39.96192 1056971.62197 +23.84770 3988.39986 37.97608 1056973.62407 +23.94790 3988.38441 36.08891 1056975.52668 +24.04810 3988.36973 34.29553 1056977.33474 +24.14830 3988.35577 32.59127 1056979.05295 +24.24850 3988.34252 30.97170 1056980.68578 +24.34870 3988.32992 29.43261 1056982.23747 +24.44890 3988.31795 27.97000 1056983.71205 +24.54910 3988.30657 26.58008 1056985.11335 +24.64930 3988.29575 25.25922 1056986.44502 +24.74950 3988.28548 24.00401 1056987.71051 +24.84970 3988.27572 22.81117 1056988.91312 +24.94990 3988.26644 21.67760 1056990.05596 +25.05010 3988.25762 20.60037 1056991.14201 +25.15030 3988.24924 19.57666 1056992.17410 +25.25050 3988.24127 18.60383 1056993.15489 +25.35070 3988.23371 17.67934 1056994.08695 +25.45090 3988.22652 16.80080 1056994.97269 +25.55110 3988.21968 15.96591 1056995.81441 +25.65130 3988.21319 15.17251 1056996.61431 +25.75150 3988.20702 14.41853 1056997.37445 +25.85170 3988.20115 13.70203 1056998.09682 +25.95190 3988.19558 13.02113 1056998.78330 +26.05210 3988.19028 12.37406 1056999.43566 +26.15230 3988.18525 11.75915 1057000.05560 +26.25251 3988.18046 11.17480 1057000.64474 +26.35271 3988.17592 10.61949 1057001.20460 +26.45291 3988.17160 10.09177 1057001.73663 +26.55311 3988.16749 9.59027 1057002.24223 +26.65331 3988.16359 9.11370 1057002.72271 +26.75351 3988.15988 8.66081 1057003.17931 +26.85371 3988.15636 8.23042 1057003.61322 +26.95391 3988.15301 7.82143 1057004.02556 +27.05411 3988.14983 7.43275 1057004.41742 +27.15431 3988.14681 7.06339 1057004.78980 +27.25451 3988.14394 6.71239 1057005.14368 +27.35471 3988.14120 6.37883 1057005.47997 +27.45491 3988.13861 6.06184 1057005.79955 +27.55511 3988.13614 5.76061 1057006.10325 +27.65531 3988.13380 5.47434 1057006.39186 +27.75551 3988.13157 5.20230 1057006.66612 +27.85571 3988.12946 4.94378 1057006.92676 +27.95591 3988.12745 4.69811 1057007.17444 +28.05611 3988.12554 4.46465 1057007.40982 +28.15631 3988.12372 4.24278 1057007.63350 +28.25651 3988.12199 4.03194 1057007.84606 +28.35671 3988.12035 3.83158 1057008.04806 +28.45691 3988.11880 3.64118 1057008.24003 +28.55711 3988.11731 3.46024 1057008.42245 +28.65731 3988.11591 3.28828 1057008.59581 +28.75752 3988.11457 3.12488 1057008.76055 +28.85772 3988.11330 2.96959 1057008.91711 +28.95792 3988.11209 2.82202 1057009.06589 +29.05812 3988.11094 2.68179 1057009.20727 +29.15832 3988.10985 2.54852 1057009.34163 +29.25852 3988.10882 2.42187 1057009.46931 +29.35872 3988.10783 2.30152 1057009.59065 +29.45892 3988.10689 2.18715 1057009.70595 +29.55912 3988.10600 2.07847 1057009.81553 +29.65932 3988.10516 1.97518 1057009.91966 +29.75952 3988.10436 1.87703 1057010.01862 +29.85972 3988.10359 1.78375 1057010.11266 +29.95992 3988.10287 1.69511 1057010.20202 +30.06012 3988.10218 1.61087 1057010.28695 +30.16032 3988.10152 1.53082 1057010.36765 +30.26052 3988.10090 1.45475 1057010.44435 +30.36072 3988.10031 1.38246 1057010.51723 +30.46092 3988.09974 1.31376 1057010.58649 +30.56112 3988.09921 1.24848 1057010.65231 +30.66132 3988.09870 1.18643 1057010.71486 +30.76152 3988.09822 1.12748 1057010.77430 +30.86172 3988.09776 1.07145 1057010.83079 +30.96192 3988.09733 1.01820 1057010.88447 +31.06212 3988.09691 0.96761 1057010.93548 +31.16232 3988.09652 0.91952 1057010.98396 +31.26253 3988.09614 0.87383 1057011.03003 +31.36273 3988.09579 0.83041 1057011.07381 +31.46293 3988.09545 0.78914 1057011.11541 +31.56313 3988.09513 0.74992 1057011.15495 +31.66333 3988.09482 0.71266 1057011.19252 +31.76353 3988.09453 0.67724 1057011.22822 +31.86373 3988.09426 0.64359 1057011.26215 +31.96393 3988.09400 0.61161 1057011.29440 +32.06413 3988.09375 0.58121 1057011.32504 +32.16433 3988.09351 0.55233 1057011.35416 +32.26453 3988.09329 0.52488 1057011.38183 +32.36473 3988.09307 0.49880 1057011.40813 +32.46493 3988.09287 0.47401 1057011.43312 +32.56513 3988.09268 0.45046 1057011.45686 +32.66533 3988.09249 0.42807 1057011.47943 +32.76553 3988.09232 0.40680 1057011.50088 +32.86573 3988.09216 0.38659 1057011.52126 +32.96593 3988.09200 0.36738 1057011.54063 +33.06613 3988.09185 0.34912 1057011.55903 +33.16633 3988.09171 0.33177 1057011.57652 +33.26653 3988.09157 0.31528 1057011.59315 +33.36673 3988.09144 0.29962 1057011.60894 +33.46693 3988.09132 0.28473 1057011.62395 +33.56713 3988.09121 0.27058 1057011.63822 +33.66733 3988.09110 0.25713 1057011.65177 +33.76754 3988.09099 0.24435 1057011.66465 +33.86774 3988.09089 0.23221 1057011.67690 +33.96794 3988.09080 0.22067 1057011.68853 +34.06814 3988.09071 0.20971 1057011.69959 +34.16834 3988.09062 0.19928 1057011.71009 +34.26854 3988.09054 0.18938 1057011.72008 +34.36874 3988.09046 0.17997 1057011.72957 +34.46894 3988.09039 0.17103 1057011.73858 +34.56914 3988.09032 0.16253 1057011.74715 +34.66934 3988.09026 0.15445 1057011.75529 +34.76954 3988.09019 0.14678 1057011.76303 +34.86974 3988.09013 0.13948 1057011.77038 +34.96994 3988.09008 0.13255 1057011.77737 +35.07014 3988.09002 0.12596 1057011.78401 +35.17034 3988.08997 0.11970 1057011.79032 +35.27054 3988.08992 0.11376 1057011.79632 +35.37074 3988.08988 0.10810 1057011.80202 +35.47094 3988.08983 0.10273 1057011.80744 +35.57114 3988.08979 0.09763 1057011.81258 +35.67134 3988.08975 0.09277 1057011.81747 +35.77154 3988.08971 0.08816 1057011.82212 +35.87174 3988.08968 0.08378 1057011.82654 +35.97194 3988.08964 0.07962 1057011.83074 +36.07214 3988.08961 0.07566 1057011.83473 +36.17234 3988.08958 0.07190 1057011.83852 +36.27255 3988.08955 0.06833 1057011.84212 +36.37275 3988.08952 0.06493 1057011.84554 +36.47295 3988.08950 0.06171 1057011.84880 +36.57315 3988.08947 0.05864 1057011.85189 +36.67335 3988.08945 0.05573 1057011.85483 +36.77355 3988.08942 0.05296 1057011.85762 +36.87375 3988.08940 0.05033 1057011.86027 +36.97395 3988.08938 0.04783 1057011.86279 +37.07415 3988.08936 0.04545 1057011.86519 +37.17435 3988.08934 0.04319 1057011.86747 +37.27455 3988.08933 0.04104 1057011.86963 +37.37475 3988.08931 0.03900 1057011.87169 +37.47495 3988.08929 0.03707 1057011.87364 +37.57515 3988.08928 0.03522 1057011.87550 +37.67535 3988.08926 0.03347 1057011.87726 +37.77555 3988.08925 0.03181 1057011.87894 +37.87575 3988.08924 0.03023 1057011.88053 +37.97595 3988.08923 0.02873 1057011.88205 +38.07615 3988.08921 0.02730 1057011.88349 +38.17635 3988.08920 0.02594 1057011.88485 +38.27655 3988.08919 0.02465 1057011.88615 +38.37675 3988.08918 0.02343 1057011.88739 +38.47695 3988.08917 0.02226 1057011.88856 +38.57715 3988.08916 0.02116 1057011.88968 +38.67735 3988.08916 0.02011 1057011.89074 +38.77756 3988.08915 0.01911 1057011.89175 +38.87776 3988.08914 0.01816 1057011.89270 +38.97796 3988.08913 0.01726 1057011.89361 +39.07816 3988.08913 0.01640 1057011.89448 +39.17836 3988.08912 0.01558 1057011.89530 +39.27856 3988.08911 0.01481 1057011.89608 +39.37876 3988.08911 0.01407 1057011.89682 +39.47896 3988.08910 0.01337 1057011.89753 +39.57916 3988.08909 0.01271 1057011.89820 +39.67936 3988.08909 0.01208 1057011.89883 +39.77956 3988.08908 0.01148 1057011.89944 +39.87976 3988.08908 0.01091 1057011.90001 +39.97996 3988.08908 0.01037 1057011.90056 +40.08016 3988.08907 0.00985 1057011.90108 +40.18036 3988.08907 0.00936 1057011.90157 +40.28056 3988.08906 0.00890 1057011.90204 +40.38076 3988.08906 0.00845 1057011.90249 +40.48096 3988.08906 0.00803 1057011.90291 +40.58116 3988.08905 0.00763 1057011.90331 +40.68136 3988.08905 0.00725 1057011.90370 +40.78156 3988.08905 0.00689 1057011.90406 +40.88176 3988.08904 0.00655 1057011.90440 +40.98196 3988.08904 0.00623 1057011.90473 +41.08216 3988.08904 0.00592 1057011.90504 +41.18236 3988.08904 0.00562 1057011.90534 +41.28257 3988.08903 0.00534 1057011.90562 +41.38277 3988.08903 0.00508 1057011.90589 +41.48297 3988.08903 0.00483 1057011.90614 +41.58317 3988.08903 0.00459 1057011.90639 +41.68337 3988.08903 0.00436 1057011.90662 +41.78357 3988.08902 0.00414 1057011.90683 +41.88377 3988.08902 0.00394 1057011.90704 +41.98397 3988.08902 0.00374 1057011.90724 +42.08417 3988.08902 0.00355 1057011.90743 +42.18437 3988.08902 0.00338 1057011.90760 +42.28457 3988.08902 0.00321 1057011.90777 +42.38477 3988.08902 0.00305 1057011.90793 +42.48497 3988.08901 0.00290 1057011.90809 +42.58517 3988.08901 0.00275 1057011.90823 +42.68537 3988.08901 0.00262 1057011.90837 +42.78557 3988.08901 0.00249 1057011.90850 +42.88577 3988.08901 0.00236 1057011.90863 +42.98597 3988.08901 0.00225 1057011.90874 +43.08617 3988.08901 0.00213 1057011.90886 +43.18637 3988.08901 0.00203 1057011.90896 +43.28657 3988.08901 0.00193 1057011.90907 +43.38677 3988.08901 0.00183 1057011.90916 +43.48697 3988.08901 0.00174 1057011.90925 +43.58717 3988.08900 0.00165 1057011.90934 +43.68737 3988.08900 0.00157 1057011.90942 +43.78758 3988.08900 0.00149 1057011.90950 +43.88778 3988.08900 0.00142 1057011.90958 +43.98798 3988.08900 0.00135 1057011.90965 +44.08818 3988.08900 0.00128 1057011.90972 +44.18838 3988.08900 0.00122 1057011.90978 +44.28858 3988.08900 0.00116 1057011.90984 +44.38878 3988.08900 0.00110 1057011.90990 +44.48898 3988.08900 0.00105 1057011.90995 +44.58918 3988.08900 0.00099 1057011.91001 +44.68938 3988.08900 0.00094 1057011.91006 +44.78958 3988.08900 0.00090 1057011.91010 +44.88978 3988.08900 0.00085 1057011.91015 +44.98998 3988.08900 0.00081 1057011.91019 +45.09018 3988.08900 0.00077 1057011.91023 +45.19038 3988.08900 0.00073 1057011.91027 +45.29058 3988.08900 0.00070 1057011.91031 +45.39078 3988.08900 0.00066 1057011.91034 +45.49098 3988.08900 0.00063 1057011.91038 +45.59118 3988.08900 0.00060 1057011.91041 +45.69138 3988.08900 0.00057 1057011.91044 +45.79158 3988.08900 0.00054 1057011.91047 +45.89178 3988.08900 0.00051 1057011.91049 +45.99198 3988.08899 0.00049 1057011.91052 +46.09218 3988.08899 0.00046 1057011.91054 +46.19238 3988.08899 0.00044 1057011.91057 +46.29259 3988.08899 0.00042 1057011.91059 +46.39279 3988.08899 0.00040 1057011.91061 +46.49299 3988.08899 0.00038 1057011.91063 +46.59319 3988.08899 0.00036 1057011.91065 +46.69339 3988.08899 0.00034 1057011.91067 +46.79359 3988.08899 0.00032 1057011.91068 +46.89379 3988.08899 0.00031 1057011.91070 +46.99399 3988.08899 0.00029 1057011.91071 +47.09419 3988.08899 0.00028 1057011.91073 +47.19439 3988.08899 0.00026 1057011.91074 +47.29459 3988.08899 0.00025 1057011.91076 +47.39479 3988.08899 0.00024 1057011.91077 +47.49499 3988.08899 0.00023 1057011.91078 +47.59519 3988.08899 0.00022 1057011.91079 +47.69539 3988.08899 0.00020 1057011.91080 +47.79559 3988.08899 0.00019 1057011.91081 +47.89579 3988.08899 0.00018 1057011.91082 +47.99599 3988.08899 0.00018 1057011.91083 +48.09619 3988.08899 0.00017 1057011.91084 +48.19639 3988.08899 0.00016 1057011.91085 +48.29659 3988.08899 0.00015 1057011.91086 +48.39679 3988.08899 0.00014 1057011.91086 +48.49699 3988.08899 0.00014 1057011.91087 +48.59719 3988.08899 0.00013 1057011.91088 +48.69739 3988.08899 0.00012 1057011.91089 +48.79760 3988.08899 0.00012 1057011.91089 +48.89780 3988.08899 0.00011 1057011.91090 +48.99800 3988.08899 0.00011 1057011.91090 +49.09820 3988.08899 0.00010 1057011.91091 +49.19840 3988.08899 0.00010 1057011.91091 +49.29860 3988.08899 0.00009 1057011.91092 +49.39880 3988.08899 0.00009 1057011.91092 +49.49900 3988.08899 0.00008 1057011.91093 +49.59920 3988.08899 0.00008 1057011.91093 +49.69940 3988.08899 0.00007 1057011.91093 +49.79960 3988.08899 0.00007 1057011.91094 +49.89980 3988.08899 0.00007 1057011.91094 +50.00000 3988.08899 0.00006 1057011.91095 \ No newline at end of file diff --git a/cpp/tests/test_flows.cpp b/cpp/tests/test_flows.cpp index 68eef35b7d..b3e1ead5ee 100644 --- a/cpp/tests/test_flows.cpp +++ b/cpp/tests/test_flows.cpp @@ -59,7 +59,7 @@ class TestModel : public mio::FlowModel public: TestModel(Populations::Index dimensions) - : Base(Populations(dimensions, 0.), mio::oseir::Parameters{}) + : Base(Populations(dimensions, 0.), mio::oseir::Parameters(mio::AgeGroup(1))) { } void get_flows(Eigen::Ref /*pop*/, Eigen::Ref /*y*/, double /*t*/, @@ -103,23 +103,23 @@ TEST(TestFlows, FlowSimulation) double tmax = 1; double dt = 0.001; - mio::oseir::Model model; + mio::oseir::Model model(1); double total_population = 10000; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = total_population - - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; // suscetible now set with every other update // params.nb_sus_t0 = params.nb_total_t0 - params.nb_exp_t0 - params.nb_inf_t0 - params.nb_rec_t0; model.parameters.set(5.2); model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + model.parameters.get().get_cont_freq_mat()[0].get_baseline().setConstant(10); model.check_constraints(); auto IC = std::make_shared(); @@ -145,23 +145,24 @@ TEST(TestFlows, CompareSimulations) double tmax = 1; double dt = 0.001; - mio::oseir::Model model; + mio::oseir::Model model(1); double total_population = 10000; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = total_population - - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; // suscetible now set with every other update // params.nb_sus_t0 = params.nb_total_t0 - params.nb_exp_t0 - params.nb_inf_t0 - params.nb_rec_t0; model.parameters.set(5.2); model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); model.check_constraints(); auto seir_sim_flows = simulate_flows(t0, tmax, dt, model); @@ -177,7 +178,7 @@ TEST(TestFlows, CompareSimulations) TEST(TestFlows, GetInitialFlows) { - mio::oseir::Model m; + mio::oseir::Model m(1); EXPECT_EQ(m.get_initial_flows().size(), 3); // 3 == Flows().size() EXPECT_EQ(m.get_initial_flows().norm(), 0); } diff --git a/cpp/tests/test_graph_simulation.cpp b/cpp/tests/test_graph_simulation.cpp index 6ae43c8864..97081da299 100644 --- a/cpp/tests/test_graph_simulation.cpp +++ b/cpp/tests/test_graph_simulation.cpp @@ -137,9 +137,9 @@ TEST(TestGraphSimulation, stopsAtTmaxStochastic) const auto tmax = 5.; const auto dt = 0.076; - mio::oseir::Model model; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 0.9; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 0.1; + mio::oseir::Model model(1); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 0.9; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 0.1; model.populations.set_total(1000); mio::Graph>, mio::MigrationEdgeStochastic> g; @@ -193,9 +193,9 @@ TEST(TestGraphSimulation, consistencyStochasticMobility) const auto tmax = 10.; const auto dt = 0.076; - mio::oseir::Model model; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 0.7; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 0.3; + mio::oseir::Model model(1); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 0.7; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 0.3; model.populations.set_total(1000); mio::Graph>, mio::MigrationEdgeStochastic> g; @@ -275,20 +275,21 @@ TEST(TestGraphSimulation, consistencyFlowMobility) double tmax = 1; double dt = 0.001; - mio::oseir::Model model; + mio::oseir::Model model(1); double total_population = 10000; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = total_population - - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; model.parameters.set(5.2); model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); model.check_constraints(); diff --git a/cpp/tests/test_mobility.cpp b/cpp/tests/test_mobility.cpp index 9a5895606c..9d4a6e98b9 100644 --- a/cpp/tests/test_mobility.cpp +++ b/cpp/tests/test_mobility.cpp @@ -38,17 +38,17 @@ TEST(TestMobility, compareNoMigrationWithSingleIntegration) auto tmax = 5; auto dt = 0.5; - mio::oseir::Model model1; - model1.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 0.9; - model1.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 0.1; + mio::oseir::Model model1(1); + model1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 0.9; + model1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 0.1; model1.populations.set_total(1000); - model1.parameters.get().get_baseline()(0, 0) = 10; + model1.parameters.get().get_cont_freq_mat()[0].get_baseline().setConstant(10); model1.parameters.set(0.4); model1.parameters.set(4); model1.parameters.set(10); auto model2 = model1; - model2.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = 1.; + model2.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 1.; model2.populations.set_total(500); auto graph_sim = mio::make_migration_sim( diff --git a/cpp/tests/test_odeseir.cpp b/cpp/tests/test_odeseir.cpp index 91070c305e..21ff5ac666 100644 --- a/cpp/tests/test_odeseir.cpp +++ b/cpp/tests/test_odeseir.cpp @@ -35,7 +35,7 @@ TEST(TestOdeSeir, simulateDefault) double tmax = 1; double dt = 0.1; - mio::oseir::Model model; + mio::oseir::Model model(1); mio::TimeSeries result = simulate(t0, tmax, dt, model); EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); @@ -43,6 +43,18 @@ TEST(TestOdeSeir, simulateDefault) class ModelTestOdeSeir : public testing::Test { + +public: + ModelTestOdeSeir() + : model(1) + { + } + double t0; + double tmax; + double dt; + double total_population; + mio::oseir::Model model; + protected: void SetUp() override { @@ -52,85 +64,33 @@ class ModelTestOdeSeir : public testing::Test total_population = 1061000; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 10000; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 1000; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 1000; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = - total_population - - this->model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - this->model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - this->model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; - // suscetible now set with every other update - // model.nb_sus_t0 = model.nb_total_t0 - model.nb_exp_t0 - model.nb_inf_t0 - model.nb_rec_t0; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 10000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 1000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 1000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; model.parameters.set(1.0); model.parameters.set(5.2); model.parameters.set(2); - model.parameters.get().get_baseline()(0, 0) = 2.7; - model.parameters.get().add_damping(0.6, mio::SimulationTime(12.5)); + mio::ContactMatrixGroup& contact_matrix = + model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(2.7); + contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5)); } - -public: - double t0; - double tmax; - double dt; - double total_population; - mio::oseir::Model model; }; -TEST_F(ModelTestOdeSeir, compareWithPreviousRun) -{ - /* - This test test the cpp model. The same test is implemented in python to compare the results of both simulations. - If this test is change the corresponding python test needs to be changed aswell (also updating the data file). - */ - std::vector> refData = load_test_data_csv("seir-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; - } - - ASSERT_NEAR(t, result.get_times()[irow], 1e-12) << "at row " << irow; - for (size_t icol = 0; icol < 4; ++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(ModelTestOdeSeir, checkPopulationConservation) { auto result = mio::simulate(t0, tmax, dt, model); - double num_persons = 0.0; - for (auto i = 0; i < result.get_last_value().size(); i++) { - num_persons += result.get_last_value()[i]; - } + double num_persons = result.get_last_value().sum(); EXPECT_NEAR(num_persons, total_population, 1e-8); } TEST_F(ModelTestOdeSeir, check_constraints_parameters) { - model.parameters.set(5.2); - model.parameters.set(6); - model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; - // 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); @@ -152,46 +112,47 @@ TEST_F(ModelTestOdeSeir, check_constraints_parameters) TEST(TestOdeSeir, apply_constraints_parameters) { const double tol_times = 1e-1; - mio::oseir::Model model; + mio::oseir::Model model(1); model.parameters.set(5.2); model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); EXPECT_EQ(model.parameters.apply_constraints(), 0); mio::set_log_level(mio::LogLevel::off); model.parameters.set(-5.2); EXPECT_EQ(model.parameters.apply_constraints(), 1); - EXPECT_EQ(model.parameters.get(), tol_times); + EXPECT_EQ(model.parameters.get()[(mio::AgeGroup)0], tol_times); model.parameters.set(1e-5); EXPECT_EQ(model.parameters.apply_constraints(), 1); - EXPECT_EQ(model.parameters.get(), tol_times); + EXPECT_EQ(model.parameters.get()[(mio::AgeGroup)0], tol_times); model.parameters.set(10.); EXPECT_EQ(model.parameters.apply_constraints(), 1); - EXPECT_NEAR(model.parameters.get(), 0.0, 1e-14); + EXPECT_NEAR(model.parameters.get()[(mio::AgeGroup)0], 0.0, 1e-14); mio::set_log_level(mio::LogLevel::warn); } TEST(TestOdeSeir, get_reproduction_numbers) { - mio::oseir::Model model; + mio::oseir::Model model(1); - double total_population = 10000; - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = - total_population - - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; + double total_population = 10000; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); model.apply_constraints(); @@ -216,6 +177,14 @@ TEST(TestOdeSeir, get_reproduction_numbers) mio::TimeSeries::Vector result_5(4); mio::TimeSeries::Vector result_6(4); + result_0.setZero(); + result_1.setZero(); + result_2.setZero(); + result_3.setZero(); + result_4.setZero(); + result_5.setZero(); + result_6.setZero(); + result_0[(Eigen::Index)mio::oseir::InfectionState::Susceptible] = 9700; result_1[(Eigen::Index)mio::oseir::InfectionState::Susceptible] = 9699.9611995799496071; result_2[(Eigen::Index)mio::oseir::InfectionState::Susceptible] = 9699.7865872644051706; @@ -238,7 +207,7 @@ TEST(TestOdeSeir, get_reproduction_numbers) EXPECT_NEAR(reproduction_numbers[i], checkReproductionNumbers[i], 1e-12); } - model.parameters.get().get_baseline()(0, 0) = 9; + contact_matrix[0].get_baseline().setConstant(9); auto reproduction_numbers2 = model.get_reproduction_numbers(result); @@ -246,7 +215,7 @@ TEST(TestOdeSeir, get_reproduction_numbers) EXPECT_NEAR(reproduction_numbers2[i], checkReproductionNumbers2[i], 1e-12); } - model.parameters.get().get_baseline()(0, 0) = 8; + contact_matrix[0].get_baseline().setConstant(8); auto reproduction_numbers3 = model.get_reproduction_numbers(result); @@ -260,21 +229,22 @@ TEST(TestOdeSeir, get_reproduction_numbers) TEST(TestOdeSeir, get_reproduction_number) { - mio::oseir::Model model; + mio::oseir::Model model(1); double total_population = 10000; //Initialize compartments to get total population of 10000 - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}] = 100; - model.populations[{mio::Index(mio::oseir::InfectionState::Susceptible)}] = - total_population - - model.populations[{mio::Index(mio::oseir::InfectionState::Exposed)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::oseir::InfectionState::Recovered)}]; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(10); model.apply_constraints(); @@ -317,11 +287,126 @@ TEST(TestOdeSeir, get_reproduction_number) EXPECT_NEAR(model.get_reproduction_number(0.7, result).value(), 2.3242860858116172196, 1e-12); EXPECT_NEAR(model.get_reproduction_number(0.0, result).value(), 2.3280000000000002913, 1e-12); - model.parameters.get().get_baseline()(0, 0) = 9; + contact_matrix[0].get_baseline().setConstant(9); EXPECT_NEAR(model.get_reproduction_number(0.1, result).value(), 2.0946073086586665113, 1e-12); EXPECT_NEAR(model.get_reproduction_number(0.3, result).value(), 2.0936545545126947765, 1e-12); - model.parameters.get().get_baseline()(0, 0) = 8; + contact_matrix[0].get_baseline().setConstant(8); EXPECT_NEAR(model.get_reproduction_number(0.2, result).value(), 1.8614409729718137676, 1e-12); EXPECT_NEAR(model.get_reproduction_number(0.9, result).value(), 1.858670429549998504, 1e-12); } + +TEST(TestSeir, get_flows) +{ + mio::oseir::Model model(1); + + constexpr double total_population = 400; + + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; + + model.parameters.set(2); + model.parameters.set(4); + model.parameters.set(1); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(1); + model.check_constraints(); + + auto dydt_default = Eigen::VectorXd(3); + dydt_default.setZero(); + auto y0 = model.get_initial_values(); + model.get_flows(y0, y0, 0, dydt_default); + + EXPECT_NEAR(dydt_default[0], 25, 1e-12); + EXPECT_NEAR(dydt_default[1], 50, 1e-12); + EXPECT_NEAR(dydt_default[2], 25, 1e-12); +} + +TEST(TestSeir, Simulation) +{ + double t0 = 0; + double tmax = 1; + double dt = 1; + + mio::oseir::Model model(1); + + constexpr double total_population = 400; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; + + model.parameters.set(2); + model.parameters.set(4); + model.parameters.set(1); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(1); + + model.check_constraints(); + + auto integrator = std::make_shared(); + + auto sim = simulate(t0, tmax, dt, model, integrator); + + EXPECT_EQ(sim.get_num_time_points(), 2); + + const auto& results_t1 = sim.get_last_value(); + EXPECT_NEAR(results_t1[0], 75, 1e-12); + EXPECT_NEAR(results_t1[1], 75, 1e-12); + EXPECT_NEAR(results_t1[2], 125, 1e-12); + EXPECT_NEAR(results_t1[3], 125, 1e-12); +} + +TEST(TestSeir, FlowSimulation) +{ + double t0 = 0; + double tmax = 1; + double dt = 1; + + mio::oseir::Model model(1); + + constexpr double total_population = 400; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}]; + + model.parameters.set(2); + model.parameters.set(4); + model.parameters.set(1); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(1); + + model.check_constraints(); + + auto integrator = std::make_shared(); + + auto sim = simulate_flows(t0, tmax, dt, model, integrator); + + // results + EXPECT_EQ(sim[0].get_num_time_points(), 2); + const auto& results_t1 = sim[0].get_last_value(); + EXPECT_NEAR(results_t1[0], 75, 1e-12); + EXPECT_NEAR(results_t1[1], 75, 1e-12); + EXPECT_NEAR(results_t1[2], 125, 1e-12); + EXPECT_NEAR(results_t1[3], 125, 1e-12); + + // flows + EXPECT_EQ(sim[1].get_num_time_points(), 2); + sim[1].print_table(); + const auto& flows_t1 = sim[1].get_last_value(); + EXPECT_NEAR(flows_t1[0], 25, 1e-12); + EXPECT_NEAR(flows_t1[1], 50, 1e-12); + EXPECT_NEAR(flows_t1[2], 25, 1e-12); +} diff --git a/cpp/tests/test_odesir.cpp b/cpp/tests/test_odesir.cpp index 1f7c9c014c..bd9233b0b1 100644 --- a/cpp/tests/test_odesir.cpp +++ b/cpp/tests/test_odesir.cpp @@ -28,6 +28,7 @@ #include #include #include +#include TEST(TestOdeSir, simulateDefault) { @@ -35,7 +36,7 @@ TEST(TestOdeSir, simulateDefault) double tmax = 1; double dt = 0.1; - mio::osir::Model model; + mio::osir::Model model(1); mio::TimeSeries result = simulate(t0, tmax, dt, model); EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); @@ -50,19 +51,25 @@ TEST(TestOdeSir, compareWithPreviousRun) double total_population = 1061000; - mio::osir::Model model; + mio::osir::Model model(1); - model.populations[{mio::Index(mio::osir::InfectionState::Infected)}] = 1000; - model.populations[{mio::Index(mio::osir::InfectionState::Recovered)}] = 1000; - model.populations[{mio::Index(mio::osir::InfectionState::Susceptible)}] = + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Susceptible)}] = total_population - - model.populations[{mio::Index(mio::osir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::osir::InfectionState::Recovered)}]; + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Infected)}] - + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Recovered)}]; model.parameters.set(1.0); model.parameters.set(2); - model.parameters.get().get_baseline()(0, 0) = 2.7; - model.parameters.get().add_damping(0.6, mio::SimulationTime(12.5)); + model.parameters.get().get_cont_freq_mat()[0].get_baseline().setConstant(2.7); + model.parameters.get().get_cont_freq_mat()[0].add_damping(0.6, + mio::SimulationTime(12.5)); std::vector> refData = load_test_data_csv("ode-sir-compare.csv"); auto integrator = std::make_shared(); @@ -105,19 +112,25 @@ TEST(TestOdeSir, checkPopulationConservation) double total_population = 1061000; - mio::osir::Model model; + mio::osir::Model model(1); - model.populations[{mio::Index(mio::osir::InfectionState::Infected)}] = 1000; - model.populations[{mio::Index(mio::osir::InfectionState::Recovered)}] = 1000; - model.populations[{mio::Index(mio::osir::InfectionState::Susceptible)}] = + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Infected)}] = 1000; + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Recovered)}] = 1000; + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Susceptible)}] = total_population - - model.populations[{mio::Index(mio::osir::InfectionState::Infected)}] - - model.populations[{mio::Index(mio::osir::InfectionState::Recovered)}]; + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Infected)}] - + model.populations[{mio::Index(0), + mio::Index(mio::osir::InfectionState::Recovered)}]; model.parameters.set(1.0); model.parameters.set(2); - model.parameters.get().get_baseline()(0, 0) = 2.7; - model.parameters.get().add_damping(0.6, mio::SimulationTime(12.5)); + model.parameters.get().get_cont_freq_mat()[0].get_baseline().setConstant(2.7); + model.parameters.get().get_cont_freq_mat()[0].add_damping(0.6, + mio::SimulationTime(12.5)); auto result = mio::simulate(t0, tmax, dt, model); double num_persons = 0.0; for (auto i = 0; i < result.get_last_value().size(); i++) { @@ -128,10 +141,10 @@ TEST(TestOdeSir, checkPopulationConservation) TEST(TestOdeSir, check_constraints_parameters) { - mio::osir::Model model; + mio::osir::Model model(1); model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + model.parameters.get().get_cont_freq_mat()[0].get_baseline().setConstant(10); // model.check_constraints() combines the functions from population and parameters. // We only want to test the functions for the parameters defined in parameters.h @@ -151,10 +164,10 @@ TEST(TestOdeSir, check_constraints_parameters) TEST(TestOdeSir, apply_constraints_parameters) { const double tol_times = 1e-1; - mio::osir::Model model; + mio::osir::Model model(1); model.parameters.set(6); model.parameters.set(0.04); - model.parameters.get().get_baseline()(0, 0) = 10; + model.parameters.get().get_cont_freq_mat()[0].get_baseline().setConstant(10); EXPECT_EQ(model.parameters.apply_constraints(), 0); @@ -162,10 +175,69 @@ TEST(TestOdeSir, apply_constraints_parameters) model.parameters.set(-2.5); EXPECT_EQ(model.parameters.apply_constraints(), 1); - EXPECT_EQ(model.parameters.get(), tol_times); + EXPECT_EQ(model.parameters.get()[(mio::AgeGroup)0], tol_times); model.parameters.set(10.); EXPECT_EQ(model.parameters.apply_constraints(), 1); - EXPECT_NEAR(model.parameters.get(), 0.0, 1e-14); + EXPECT_NEAR(model.parameters.get()[(mio::AgeGroup)0], 0.0, 1e-14); mio::set_log_level(mio::LogLevel::warn); +} + +TEST(Testsir, get_derivatives) +{ + mio::osir::Model model(1); + + constexpr auto total_population = 400; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}]; + + model.parameters.set(4); + model.parameters.set(1); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(1); + model.check_constraints(); + + auto dydt_default = Eigen::VectorXd(Eigen::Index(mio::osir::InfectionState::Count)); + dydt_default.setZero(); + auto y0 = model.get_initial_values(); + model.get_derivatives(y0, y0, 0, dydt_default); + + EXPECT_NEAR(dydt_default[0], -50, 1e-12); + EXPECT_NEAR(dydt_default[1], 25, 1e-12); + EXPECT_NEAR(dydt_default[2], 25, 1e-12); +} + +TEST(Testsir, Simulation) +{ + double t0 = 0; + double tmax = 1; + double dt = 1; + + mio::osir::Model model(1); + + constexpr auto total_population = 400; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] = 100; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}] = 100; + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Susceptible}] = + total_population - model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] - + model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}]; + + model.parameters.set(4); + model.parameters.set(1); + mio::ContactMatrixGroup& contact_matrix = model.parameters.get().get_cont_freq_mat(); + contact_matrix[0].get_baseline().setConstant(1); + model.check_constraints(); + + auto integrator = std::make_shared(); + auto sim = simulate(t0, tmax, dt, model, integrator); + + EXPECT_EQ(sim.get_num_time_points(), 2); + + const auto& results_t1 = sim.get_last_value(); + EXPECT_NEAR(results_t1[0], 150, 1e-12); + EXPECT_NEAR(results_t1[1], 125, 1e-12); + EXPECT_NEAR(results_t1[2], 125, 1e-12); } \ No newline at end of file diff --git a/pycode/examples/simulation/oseir_simple.py b/pycode/examples/simulation/oseir_simple.py index 96bfe779ab..4d01d40e29 100644 --- a/pycode/examples/simulation/oseir_simple.py +++ b/pycode/examples/simulation/oseir_simple.py @@ -21,7 +21,7 @@ import numpy as np -from memilio.simulation import Damping +from memilio.simulation import AgeGroup, Damping from memilio.simulation.oseir import Index_InfectionState from memilio.simulation.oseir import InfectionState as State from memilio.simulation.oseir import (Model, interpolate_simulation_result, @@ -40,26 +40,29 @@ def run_oseir_simulation(): dt = 0.1 # Initialize Parameters - model = Model() + num_groups = 1 + model = Model(num_groups) + A0 = AgeGroup(0) # Compartment transition duration - model.parameters.TimeExposed.value = 5.2 - model.parameters.TimeInfected.value = 6. + model.parameters.TimeExposed[A0] = 5.2 + model.parameters.TimeInfected[A0] = 6. # Compartment transition propabilities - model.parameters.TransmissionProbabilityOnContact.value = 1. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. # Initial number of people in each compartment - model.populations[Index_InfectionState(State.Exposed)] = 100 - model.populations[Index_InfectionState(State.Infected)] = 50 - model.populations[Index_InfectionState(State.Recovered)] = 10 + model.populations[A0, State.Exposed] = 100 + model.populations[A0, State.Infected] = 50 + model.populations[A0, State.Recovered] = 10 model.populations.set_difference_from_total( - (Index_InfectionState(State.Susceptible)), populations[0]) + (A0, State.Susceptible), populations[0]) - # model.parameters.ContactPatterns = ContactMatrix(np.r_[0.5]) - model.parameters.ContactPatterns.baseline = np.ones((1, 1)) - model.parameters.ContactPatterns.minimum = np.zeros((1, 1)) - model.parameters.ContactPatterns.add_damping( + model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones( + (num_groups, num_groups)) + model.parameters.ContactPatterns.cont_freq_mat[0].minimum = np.zeros( + (num_groups, num_groups)) + model.parameters.ContactPatterns.cont_freq_mat.add_damping( Damping(coeffs=np.r_[0.9], t=30.0, level=0, type=0)) # Check logical constraints to parameters diff --git a/pycode/examples/simulation/osir_simple.py b/pycode/examples/simulation/osir_simple.py index 336603e7cc..1f28c6f809 100644 --- a/pycode/examples/simulation/osir_simple.py +++ b/pycode/examples/simulation/osir_simple.py @@ -21,7 +21,7 @@ import numpy as np -from memilio.simulation import Damping +from memilio.simulation import AgeGroup, Damping from memilio.simulation.osir import Index_InfectionState from memilio.simulation.osir import InfectionState as State from memilio.simulation.osir import (Model, interpolate_simulation_result, @@ -40,24 +40,27 @@ def run_osir_simulation(): dt = 0.1 # Initialize Parameters - model = Model() + num_groups = 1 + model = Model(num_groups) + A0 = AgeGroup(0) # Compartment transition duration - model.parameters.TimeInfected.value = 6. + model.parameters.TimeInfected[A0] = 6. # Compartment transition propabilities - model.parameters.TransmissionProbabilityOnContact.value = 1. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. # Initial number of people in each compartment - model.populations[Index_InfectionState(State.Infected)] = 50 - model.populations[Index_InfectionState(State.Recovered)] = 10 + model.populations[A0, State.Infected] = 50 + model.populations[A0, State.Recovered] = 10 model.populations.set_difference_from_total( - (Index_InfectionState(State.Susceptible)), populations[0]) + (A0, State.Susceptible), populations[0]) - # model.parameters.ContactPatterns = ContactMatrix(np.r_[0.5]) - model.parameters.ContactPatterns.baseline = np.ones((1, 1)) - model.parameters.ContactPatterns.minimum = np.zeros((1, 1)) - model.parameters.ContactPatterns.add_damping( + model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones( + (num_groups, num_groups)) + model.parameters.ContactPatterns.cont_freq_mat[0].minimum = np.zeros( + (num_groups, num_groups)) + model.parameters.ContactPatterns.cont_freq_mat.add_damping( Damping(coeffs=np.r_[0.9], t=30.0, level=0, type=0)) # Check logical constraints to parameters diff --git a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt index 9d230ecb5f..432b7e149a 100644 --- a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt +++ b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt @@ -36,6 +36,12 @@ namespace py = pybind11; namespace pymio{ //specialization of pretty_name +template <> +std::string pretty_name() +{ + return "AgeGroup"; +} + template <> std::string pretty_name() { @@ -65,12 +71,12 @@ PYBIND11_MODULE(_simulation_test_oseir, m) .def("check_constraints", &mio::oseir::Parameters::check_constraints) .def("apply_constraints", &mio::oseir::Parameters::apply_constraints); - using Populations = mio::Populations; + using Populations = mio::Populations; pymio::bind_Population(m, "Population", mio::Tag{}); - pymio::bind_CompartmentalModel, mio::oseir::Parameters>(m, "ModelBase"); - py::class_, mio::oseir::Parameters, mio::TypeList, mio::Flow, mio::Flow>>>(m, "Model") - .def(py::init<>()); + pymio::bind_CompartmentalModel, mio::oseir::Parameters>(m, "ModelBase"); + py::class_, mio::oseir::Parameters, mio::TypeList, mio::Flow, mio::Flow>>>(m, "Model") + .def(py::init(), py::arg("num_agegroups")); diff --git a/pycode/memilio-simulation/memilio/simulation/oseir.cpp b/pycode/memilio-simulation/memilio/simulation/oseir.cpp index 683c69112c..631402deb4 100755 --- a/pycode/memilio-simulation/memilio/simulation/oseir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/oseir.cpp @@ -41,6 +41,12 @@ std::string pretty_name() return "InfectionState"; } +template <> +std::string pretty_name() +{ + return "AgeGroup"; +} + } // namespace pymio PYBIND11_MODULE(_simulation_oseir, m) @@ -66,15 +72,16 @@ PYBIND11_MODULE(_simulation_oseir, m) pymio::bind_ParameterSet(m, "ParametersBase"); py::class_(m, "Parameters") - .def(py::init<>()) + .def(py::init()) .def("check_constraints", &mio::oseir::Parameters::check_constraints); - using Populations = mio::Populations; - pymio::bind_Population(m, "Population", mio::Tag{}); - pymio::bind_CompartmentalModel(m, "ModelBase"); + using SeirPopulations = mio::Populations; + pymio::bind_Population(m, "SeirPopulations", mio::Tag{}); + + pymio::bind_CompartmentalModel(m, "ModelBase"); py::class_>(m, "Model") - .def(py::init<>()); + mio::CompartmentalModel>(m, "Model") + .def(py::init(), py::arg("num_agegroups")); m.def( "simulate", diff --git a/pycode/memilio-simulation/memilio/simulation/osir.cpp b/pycode/memilio-simulation/memilio/simulation/osir.cpp index 615196c954..1987f34563 100644 --- a/pycode/memilio-simulation/memilio/simulation/osir.cpp +++ b/pycode/memilio-simulation/memilio/simulation/osir.cpp @@ -40,6 +40,12 @@ std::string pretty_name() return "InfectionState"; } +template <> +std::string pretty_name() +{ + return "AgeGroup"; +} + } // namespace pymio PYBIND11_MODULE(_simulation_osir, m) @@ -64,15 +70,16 @@ PYBIND11_MODULE(_simulation_osir, m) pymio::bind_ParameterSet(m, "ParametersBase"); py::class_(m, "Parameters") - .def(py::init<>()) + .def(py::init()) .def("check_constraints", &mio::osir::Parameters::check_constraints); - using Populations = mio::Populations; - pymio::bind_Population(m, "Population", mio::Tag{}); - pymio::bind_CompartmentalModel(m, "ModelBase"); + using SirPopulations = mio::Populations; + pymio::bind_Population(m, "SirPopulations", mio::Tag{}); + + pymio::bind_CompartmentalModel(m, "ModelBase"); py::class_>(m, "Model") - .def(py::init<>()); + mio::CompartmentalModel>(m, "Model") + .def(py::init(), py::arg("num_agegroups")); m.def( "simulate", diff --git a/pycode/memilio-simulation/memilio/simulation_test/test_oseir.py b/pycode/memilio-simulation/memilio/simulation_test/test_oseir.py index 08928fab28..d97021fe7d 100644 --- a/pycode/memilio-simulation/memilio/simulation_test/test_oseir.py +++ b/pycode/memilio-simulation/memilio/simulation_test/test_oseir.py @@ -23,7 +23,7 @@ import numpy as np import pandas as pd -from memilio.simulation import Damping +from memilio.simulation import AgeGroup, Damping from memilio.simulation.oseir import Index_InfectionState from memilio.simulation.oseir import InfectionState as State from memilio.simulation.oseir import Model, simulate, simulate_flows @@ -35,27 +35,30 @@ class Test_oseir_integration(unittest.TestCase): def setUp(self): - model = Model() + model = Model(1) + A0 = AgeGroup(0) self.t0 = 0. self.tmax = 50. self.dt = 0.1002004008016032 total_population = 1061000 - model.populations[Index_InfectionState(State.Exposed)] = 10000 - model.populations[Index_InfectionState(State.Infected)] = 1000 - model.populations[Index_InfectionState(State.Recovered)] = 1000 + model.populations[A0, State.Exposed] = 10000 + model.populations[A0, State.Infected] = 1000 + model.populations[A0, State.Recovered] = 1000 model.populations.set_difference_from_total( - Index_InfectionState(State.Susceptible), total_population) + (A0, State.Susceptible), total_population) - model.parameters.TransmissionProbabilityOnContact.value = 1. - model.parameters.TimeExposed.value = 5.2 - model.parameters.TimeInfected.value = 2. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. + model.parameters.TimeExposed[A0] = 5.2 + model.parameters.TimeInfected[A0] = 2. - model.parameters.ContactPatterns.baseline = [[2.7]] - model.parameters.ContactPatterns.minimum = np.zeros((1, 1)) - model.parameters.ContactPatterns.add_damping( - Damping(coeffs=[[0.6]], t=12.5, level=0, type=0)) + model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones( + (1, 1)) * 2.7 + model.parameters.ContactPatterns.cont_freq_mat[0].minimum = np.zeros( + (1, 1)) + model.parameters.ContactPatterns.cont_freq_mat.add_damping( + Damping(coeffs=np.r_[0.6], t=12.5, level=0, type=0)) model.check_constraints() @@ -97,7 +100,7 @@ def test_compare_with_cpp(self): self.assertAlmostEqual( t, result.get_time(index_timestep), - delta=1e-12) + delta=1e-10) for index_compartment in range(0, 4): ref = timestep[index_compartment+1] @@ -121,26 +124,27 @@ def test_flow_simulation_simple(self): def test_check_constraints_parameters(self): - model = Model() + model = Model(1) + A0 = AgeGroup(0) - model.parameters.TimeExposed.value = 5.2 - model.parameters.TimeInfected.value = 6. - model.parameters.TransmissionProbabilityOnContact.value = 1. + model.parameters.TimeExposed[A0] = 5.2 + model.parameters.TimeInfected[A0] = 6. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. - model.parameters.TimeExposed.value = 5.2 - model.parameters.TimeInfected.value = 6. - model.parameters.TransmissionProbabilityOnContact.value = 1. + model.parameters.TimeExposed[A0] = 5.2 + model.parameters.TimeInfected[A0] = 6. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. self.assertEqual(model.parameters.check_constraints(), 0) - model.parameters.TimeExposed.value = -1. + model.parameters.TimeExposed[A0] = -1. self.assertEqual(model.parameters.check_constraints(), 1) - model.parameters.TimeExposed.value = 5.2 - model.parameters.TimeInfected.value = 0 + model.parameters.TimeExposed[A0] = 5.2 + model.parameters.TimeInfected[A0] = 0 self.assertEqual(model.parameters.check_constraints(), 1) - model.parameters.TimeInfected.value = 6. - model.parameters.TransmissionProbabilityOnContact.value = -1. + model.parameters.TimeInfected[A0] = 6. + model.parameters.TransmissionProbabilityOnContact[A0] = -1. self.assertEqual(model.parameters.check_constraints(), 1) diff --git a/pycode/memilio-simulation/memilio/simulation_test/test_osir.py b/pycode/memilio-simulation/memilio/simulation_test/test_osir.py index d777ac9b7d..1e1bd1cc1a 100644 --- a/pycode/memilio-simulation/memilio/simulation_test/test_osir.py +++ b/pycode/memilio-simulation/memilio/simulation_test/test_osir.py @@ -21,7 +21,7 @@ import numpy as np -from memilio.simulation import Damping +from memilio.simulation import AgeGroup, Damping from memilio.simulation.osir import Index_InfectionState from memilio.simulation.osir import InfectionState as State from memilio.simulation.osir import Model, simulate @@ -31,18 +31,21 @@ class Test_osir_integration(unittest.TestCase): def setUp(self): - model = Model() + model = Model(1) + A0 = AgeGroup(0) - model.parameters.TimeInfected.value = 6. - model.parameters.TransmissionProbabilityOnContact.value = 1. + model.parameters.TimeInfected[A0] = 6. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. - model.populations[Index_InfectionState(State.Susceptible)] = 4800 - model.populations[Index_InfectionState(State.Infected)] = 50 - model.populations[Index_InfectionState(State.Recovered)] = 50 + model.populations[A0, State.Susceptible] = 4800 + model.populations[A0, State.Infected] = 50 + model.populations[A0, State.Recovered] = 50 - model.parameters.ContactPatterns.baseline = np.ones((1, 1)) - model.parameters.ContactPatterns.minimum = np.zeros((1, 1)) - model.parameters.ContactPatterns.add_damping( + model.parameters.ContactPatterns.cont_freq_mat[0].baseline = np.ones( + (1, 1)) + model.parameters.ContactPatterns.cont_freq_mat[0].minimum = np.zeros( + (1, 1)) + model.parameters.ContactPatterns.cont_freq_mat.add_damping( Damping(coeffs=np.r_[0.9], t=30.0, level=0, type=0)) model.check_constraints() @@ -57,20 +60,21 @@ def test_simulate_simple(self): def test_check_constraints_parameters(self): - model = Model() + model = Model(1) + A0 = AgeGroup(0) - model.parameters.TimeInfected.value = 6. - model.parameters.TransmissionProbabilityOnContact.value = 1. + model.parameters.TimeInfected[A0] = 6. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. - model.parameters.TimeInfected.value = 6. - model.parameters.TransmissionProbabilityOnContact.value = 1. + model.parameters.TimeInfected[A0] = 6. + model.parameters.TransmissionProbabilityOnContact[A0] = 1. self.assertEqual(model.parameters.check_constraints(), 0) - model.parameters.TimeInfected.value = 0 + model.parameters.TimeInfected[A0] = 0 self.assertEqual(model.parameters.check_constraints(), 1) - model.parameters.TimeInfected.value = 6. - model.parameters.TransmissionProbabilityOnContact.value = -1. + model.parameters.TimeInfected[A0] = 6. + model.parameters.TransmissionProbabilityOnContact[A0] = -1. self.assertEqual(model.parameters.check_constraints(), 1)