Skip to content

Commit f59a074

Browse files
MaxBetzDLRMaxBetzDLR
authored
1321 copy constructable simulation (#1345)
Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com> Main changes: - Move handling of IntegratorCore outside of Simulation class into member pointer integrator - Change IntegratorCore to be unique in integrator - Add Copy Constr/Assignment
1 parent e5a27fc commit f59a074

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

47 files changed

+370
-226
lines changed

‎cpp/benchmarks/flow_simulation_ode_secirvvs.cpp‎

Lines changed: 22 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -37,15 +37,18 @@ void flowless_sim(::benchmark::State& state)
3737
// create model
3838
Model model(cfg.num_agegroups);
3939
mio::benchmark::setup_model(model);
40-
// create simulation
41-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
42-
std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
43-
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
4440
// run benchmark
4541
for (auto _ : state) {
42+
// create simulation
43+
// exclude integrator creation from benchmark
44+
state.PauseTiming();
45+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> I =
46+
std::make_unique<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
47+
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
48+
state.ResumeTiming();
4649
// This code gets timed
4750
mio::benchmark::Simulation<mio::Simulation<ScalarType, Model>> sim(model, cfg.t0, cfg.dt);
48-
sim.set_integrator(I);
51+
sim.set_integrator_core(std::move(I));
4952
// run sim
5053
sim.advance(cfg.t_max);
5154
}
@@ -62,15 +65,17 @@ void flow_sim_comp_only(::benchmark::State& state)
6265
// create model
6366
Model model(cfg.num_agegroups);
6467
mio::benchmark::setup_model(model);
65-
// create simulation
66-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
67-
std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
68-
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
6968
// run benchmark
7069
for (auto _ : state) {
70+
// create simulation
71+
state.PauseTiming();
72+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> I =
73+
std::make_unique<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
74+
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
75+
state.ResumeTiming();
7176
// This code gets timed
7277
mio::osecirvvs::Simulation<ScalarType, mio::Simulation<ScalarType, Model>> sim(model, cfg.t0, cfg.dt);
73-
sim.set_integrator(I);
78+
sim.set_integrator_core(std::move(I));
7479
// run sim
7580
sim.advance(cfg.t_max);
7681
}
@@ -87,15 +92,17 @@ void flow_sim(::benchmark::State& state)
8792
// create model
8893
Model model(cfg.num_agegroups);
8994
mio::benchmark::setup_model(model);
90-
// create simulation
91-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
92-
std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
93-
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
9495
// run benchmark
9596
for (auto _ : state) {
97+
// create simulation
98+
state.PauseTiming();
99+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> I =
100+
std::make_unique<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
101+
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
102+
state.ResumeTiming();
96103
// This code gets timed
97104
mio::osecirvvs::Simulation<ScalarType, mio::FlowSimulation<ScalarType, Model>> sim(model, cfg.t0, cfg.dt);
98-
sim.set_integrator(I);
105+
sim.set_integrator_core(std::move(I));
99106
// run sim
100107
sim.advance(cfg.t_max);
101108
}

‎cpp/benchmarks/flow_simulation_ode_seir.cpp‎

Lines changed: 22 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -124,15 +124,18 @@ void flowless_sim(::benchmark::State& state)
124124
// create model
125125
Model model(cfg.num_agegroups);
126126
mio::benchmark::setup_model(model);
127-
// create simulation
128-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
129-
std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
130-
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
131127
mio::TimeSeries<ScalarType> results(static_cast<size_t>(Model::Compartments::Count));
132128
// run benchmark
133129
for (auto _ : state) {
130+
// create simulation
131+
// exclude integrator creation from benchmark
132+
state.PauseTiming();
133+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> I =
134+
std::make_unique<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
135+
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
136+
state.ResumeTiming();
134137
// This code gets timed
135-
results = mio::simulate(cfg.t0, cfg.t_max, cfg.dt, model, I);
138+
results = mio::simulate(cfg.t0, cfg.t_max, cfg.dt, model, std::move(I));
136139
}
137140
}
138141

@@ -147,15 +150,17 @@ void flow_sim_comp_only(::benchmark::State& state)
147150
// create model
148151
Model model(cfg.num_agegroups);
149152
mio::benchmark::setup_model(model);
150-
// create simulation
151-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
152-
std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
153-
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
154153
mio::TimeSeries<ScalarType> results(static_cast<size_t>(Model::Compartments::Count));
155154
// run benchmark
156155
for (auto _ : state) {
156+
// create simulation
157+
state.PauseTiming();
158+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> I =
159+
std::make_unique<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
160+
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
161+
state.ResumeTiming();
157162
// This code gets timed
158-
results = mio::simulate(cfg.t0, cfg.t_max, cfg.dt, model, I);
163+
results = mio::simulate(cfg.t0, cfg.t_max, cfg.dt, model, std::move(I));
159164
}
160165
}
161166

@@ -170,15 +175,17 @@ void flow_sim(::benchmark::State& state)
170175
// create model
171176
Model model(cfg.num_agegroups);
172177
mio::benchmark::setup_model(model);
173-
// create simulation
174-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
175-
std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
176-
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
177178
mio::TimeSeries<ScalarType> results(static_cast<size_t>(Model::Compartments::Count));
178179
// run benchmark
179180
for (auto _ : state) {
181+
// create simulation
182+
state.PauseTiming();
183+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> I =
184+
std::make_unique<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
185+
cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
186+
state.ResumeTiming();
180187
// This code gets timed
181-
results = mio::simulate_flows(cfg.t0, cfg.t_max, cfg.dt, model, I)[0];
188+
results = mio::simulate_flows(cfg.t0, cfg.t_max, cfg.dt, model, std::move(I))[0];
182189
}
183190
}
184191

‎cpp/benchmarks/graph_simulation.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ auto create_simulation()
122122
g;
123123
for (size_t county_id = 0; county_id < cfg.num_regions; county_id++) {
124124
g.add_node(county_id, model, cfg.t0);
125-
g.nodes()[county_id].property.get_simulation().set_integrator(std::make_shared<Integrator>());
125+
g.nodes()[county_id].property.get_simulation().set_integrator_core(std::make_unique<Integrator>());
126126
}
127127

128128
// Graph is always complete here

‎cpp/benchmarks/simulation.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -35,9 +35,9 @@ void simulation(::benchmark::State& state)
3535

3636
for (auto _ : state) {
3737
// This code gets timed
38-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
39-
std::make_shared<Integrator>(cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
40-
simulate(cfg.t0, cfg.t_max, cfg.dt, model, I);
38+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> I =
39+
std::make_unique<Integrator>(cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
40+
simulate(cfg.t0, cfg.t_max, cfg.dt, model, std::move(I));
4141
}
4242
}
4343

‎cpp/examples/ode_secir.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -83,12 +83,12 @@ int main()
8383
Example of using a different integrator
8484
All available integrators are listed in cpp/memilio/math/README.md
8585
86-
auto integrator = std::make_shared<mio::RKIntegratorCore>();
86+
auto integrator = std::make_unique<mio::RKIntegratorCore>();
8787
integrator->set_dt_min(0.3);
8888
integrator->set_dt_max(1.0);
8989
integrator->set_rel_tolerance(1e-4);
9090
integrator->set_abs_tolerance(1e-1);
91-
mio::TimeSeries<double> secir = simulate(t0, tmax, dt, model, integrator);
91+
mio::TimeSeries<double> secir = simulate<double>(t0, tmax, dt, model, std::move(integrator));
9292
*/
9393

9494
bool print_to_terminal = true;

‎cpp/examples/ode_secir_contact_changes.cpp‎

Lines changed: 9 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -40,8 +40,6 @@ int main()
4040

4141
double nb_total_t0 = 1000, nb_inf_t0 = 10;
4242

43-
auto integrator = std::make_shared<mio::EulerIntegratorCore<ScalarType>>();
44-
4543
// default model run to be compared against
4644
mio::osecir::Model model_a(1);
4745
const auto indx_flow_SE =
@@ -56,7 +54,9 @@ int main()
5654
// set probability of transmission and risk of infection to 1.
5755
model_a.parameters.get<mio::osecir::TransmissionProbabilityOnContact<ScalarType>>() = 1.0;
5856
model_a.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic<ScalarType>>() = 1.0;
59-
auto result_a = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_a, integrator);
57+
58+
mio::EulerIntegratorCore<ScalarType> integrator;
59+
auto result_a = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_a, integrator.clone());
6060
result_a[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
6161
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
6262
"I_Crit->D", "I_Crit->R"},
@@ -73,7 +73,8 @@ int main()
7373
mio::ContactMatrixGroup& contact_matrix_b = model_b.parameters.get<mio::osecir::ContactPatterns<ScalarType>>();
7474
contact_matrix_b[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
7575
contact_matrix_b[0].add_damping(0.5, mio::SimulationTime(0.)); // contact reduction happens here!
76-
auto result_b = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_b, integrator);
76+
77+
auto result_b = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_b, integrator.clone());
7778
result_b[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
7879
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
7980
"I_Crit->D", "I_Crit->R"},
@@ -91,7 +92,8 @@ int main()
9192
mio::ContactMatrixGroup& contact_matrix_c = model_c.parameters.get<mio::osecir::ContactPatterns<ScalarType>>();
9293
contact_matrix_c[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
9394
contact_matrix_c[0].add_damping(1., mio::SimulationTime(0.)); // contact reduction happens here!
94-
auto result_c = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_c, integrator);
95+
96+
auto result_c = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_c, integrator.clone());
9597
result_c[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
9698
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
9799
"I_Crit->D", "I_Crit->R"},
@@ -109,7 +111,8 @@ int main()
109111
mio::ContactMatrixGroup& contact_matrix_d = model_d.parameters.get<mio::osecir::ContactPatterns<ScalarType>>();
110112
contact_matrix_d[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
111113
contact_matrix_d[0].add_damping(-1., mio::SimulationTime(0.)); // contact increase happens here!
112-
auto result_d = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_d, integrator);
114+
115+
auto result_d = mio::simulate_flows<ScalarType>(t0, tmax, dt, model_d, integrator.clone());
113116
result_d[1].print_table({"S->E", "E->I_NS", "I_NS->I_Sy", "I_NS->R", "I_NSC->I_SyC", "I_NSC->R", "I_Sy->I_Sev",
114117
"I_Sy->R", "I_SyC->I_Sev", "I_SyC->R", "I_Sev->I_Crit", "I_Sev->R", "I_Sev->D",
115118
"I_Crit->D", "I_Crit->R"},

‎cpp/examples/ode_secirvvs.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -121,12 +121,12 @@ int main()
121121
model.apply_constraints();
122122

123123
// use adaptive Runge-Kutta-Fehlberg45 scheme as integrator
124-
// auto integrator = std::make_shared<mio::RKIntegratorCore>();
124+
// auto integrator = std::make_unique<mio::RKIntegratorCore>();
125125
// integrator->set_dt_min(0.3);
126126
// integrator->set_dt_max(1.0);
127127
// integrator->set_rel_tolerance(1e-4);
128128
// integrator->set_abs_tolerance(1e-1);
129-
// mio::TimeSeries<double> result = mio::osecirvvs::simulate(t0, tmax, dt, model, integrator);
129+
// mio::TimeSeries<double> result = mio::osecirvvs::simulate<double>(t0, tmax, dt, model, std::move(integrator));
130130

131131
// use default Cash-Karp adaptive integrator
132132
mio::TimeSeries<double> result = mio::osecirvvs::simulate<double>(t0, tmax, dt, model);

‎cpp/examples/ode_sir.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -56,9 +56,9 @@ int main()
5656
contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5));
5757
model.check_constraints();
5858

59-
std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> integrator =
60-
std::make_shared<mio::EulerIntegratorCore<ScalarType>>();
61-
auto sir = simulate(t0, tmax, dt, model, integrator);
59+
std::unique_ptr<mio::OdeIntegratorCore<ScalarType>> integrator =
60+
std::make_unique<mio::EulerIntegratorCore<ScalarType>>();
61+
auto sir = simulate(t0, tmax, dt, model, std::move(integrator));
6262

6363
// interpolate results
6464
auto interpolated_results = mio::interpolate_simulation_result(sir);

‎cpp/memilio/compartments/flow_simulation.h‎

100644100755
Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ class FlowSimulation : public details::FlowSimulationBase<FP, M, OdeIntegrator<F
4646
* @param[in] dt Initial step size of integration.
4747
*/
4848
FlowSimulation(Model const& model, FP t0 = 0., FP dt = 0.1)
49-
: Base(model, std::make_shared<DefaultIntegratorCore<FP>>(), t0, dt)
49+
: Base(model, std::make_unique<DefaultIntegratorCore<FP>>(), t0, dt)
5050
, m_pop(model.get_initial_values().size())
5151
{
5252
}
@@ -99,7 +99,7 @@ class FlowSimulation : public details::FlowSimulationBase<FP, M, OdeIntegrator<F
9999
* @param[in] tmax End time.
100100
* @param[in] dt Initial step size of integration.
101101
* @param[in] model An instance of a FlowModel.
102-
* @param[in] integrator Optionally override the IntegratorCore used by the FlowSimulation.
102+
* @param[in] integrator_core Optionally override the IntegratorCore used by the FlowSimulation.
103103
* @return The simulation result as two TimeSeries. The first describes the compartments at each time point,
104104
* the second gives the corresponding flows that lead from t0 to each time point.
105105
* @tparam FP a floating point type, e.g., double
@@ -108,12 +108,12 @@ class FlowSimulation : public details::FlowSimulationBase<FP, M, OdeIntegrator<F
108108
*/
109109
template <typename FP, class Model, class Sim = FlowSimulation<FP, Model>>
110110
std::vector<TimeSeries<FP>> simulate_flows(FP t0, FP tmax, FP dt, Model const& model,
111-
std::shared_ptr<OdeIntegratorCore<FP>> integrator = nullptr)
111+
std::unique_ptr<OdeIntegratorCore<FP>>&& integrator_core = nullptr)
112112
{
113113
model.check_constraints();
114114
Sim sim(model, t0, dt);
115-
if (integrator) {
116-
sim.set_integrator(integrator);
115+
if (integrator_core) {
116+
sim.set_integrator_core(std::move(integrator_core));
117117
}
118118
sim.advance(tmax);
119119
return {sim.get_result(), sim.get_flows()};

‎cpp/memilio/compartments/flow_simulation_base.h‎

100644100755
Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -50,11 +50,12 @@ class FlowSimulationBase : public SimulationBase<FP, M, Integrands...>
5050
/**
5151
* @brief Create a FlowSimulationBase.
5252
* @param[in] model An instance of a flow model.
53+
* @param[in] integrator_core A unique pointer to an object derived from IntegratorCore.
5354
* @param[in] t0 Start time.
5455
* @param[in] dt Initial step size of integration.
5556
*/
56-
FlowSimulationBase(Model const& model, std::shared_ptr<Core> integrator, FP t0, FP dt)
57-
: Base(model, integrator, t0, dt)
57+
FlowSimulationBase(Model const& model, std::unique_ptr<Core>&& integrator_core, FP t0, FP dt)
58+
: Base(model, std::move(integrator_core), t0, dt)
5859
, m_flow_result(t0, model.get_initial_flows())
5960
{
6061
}

0 commit comments

Comments
 (0)