Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 MEmilio
3 : *
4 : * Authors: Nils Wassmuth, Rene Schmieding, Martin J. Kuehn
5 : *
6 : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
7 : *
8 : * Licensed under the Apache License, Version 2.0 (the "License");
9 : * you may not use this file except in compliance with the License.
10 : * You may obtain a copy of the License at
11 : *
12 : * http://www.apache.org/licenses/LICENSE-2.0
13 : *
14 : * Unless required by applicable law or agreed to in writing, software
15 : * distributed under the License is distributed on an "AS IS" BASIS,
16 : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 : * See the License for the specific language governing permissions and
18 : * limitations under the License.
19 : */
20 : #ifndef MIO_SDE_SIR_SIMULATION_H
21 : #define MIO_SDE_SIR_SIMULATION_H
22 :
23 : #include "memilio/compartments/flow_simulation.h"
24 : #include "memilio/compartments/simulation.h"
25 : #include "memilio/math/euler.h"
26 : #include "sde_sir/model.h"
27 :
28 : namespace mio
29 : {
30 : namespace ssir
31 : {
32 :
33 : /// @brief A specialized Simulation for mio::ssir::Model.
34 : class Simulation : public mio::Simulation<ScalarType, Model>
35 : {
36 : protected:
37 : using mio::Simulation<ScalarType, Model>::set_integrator;
38 :
39 : public:
40 : /**
41 : * @brief Set up the simulation with an ODE solver.
42 : * @param[in] model An instance of mio::ssir::Model.
43 : * @param[in] t0 Start time.
44 : * @param[in] dt Initial step size of integration.
45 : */
46 9 : Simulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 0.1)
47 9 : : mio::Simulation<ScalarType, Model>(model, t0, dt)
48 : {
49 9 : auto integrator = std::make_shared<mio::EulerIntegratorCore<ScalarType>>();
50 9 : set_integrator(integrator);
51 18 : }
52 :
53 : using mio::Simulation<ScalarType, Model>::Simulation;
54 : /**
55 : * @brief advance simulation to tmax
56 : * tmax must be greater than get_result().get_last_time_point()
57 : * @param tmax next stopping point of simulation
58 : */
59 18 : Eigen::Ref<Vector<>> advance(ScalarType tmax)
60 : {
61 18 : return get_ode_integrator().advance(
62 18 : [this](auto&& y, auto&& t, auto&& dydt) {
63 18 : get_model().step_size = get_dt();
64 18 : get_model().get_derivatives(y, y, t, dydt);
65 18 : },
66 18 : tmax, get_dt(), get_result());
67 : }
68 : };
69 :
70 : /**
71 : * @brief Run a Simulation of a mio::ssir::Model.
72 : * @param[in] t0 Start time.
73 : * @param[in] tmax End time.
74 : * @param[in] dt Initial step size of integration.
75 : * @param[in] model An instance of mio::ssir::Model.
76 : * @return A TimeSeries to represent the final simulation result
77 : */
78 0 : inline TimeSeries<ScalarType> simulate(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& model)
79 : {
80 0 : model.check_constraints();
81 0 : Simulation sim(model, t0, dt);
82 0 : sim.advance(tmax);
83 0 : return sim.get_result();
84 0 : }
85 :
86 : /// @brief A specialized FlowSimulation for mio::ssir::Model.
87 : class FlowSimulation : public mio::FlowSimulation<ScalarType, Model>
88 : {
89 : protected:
90 : using mio::FlowSimulation<ScalarType, Model>::set_integrator;
91 :
92 : public:
93 : /**
94 : * @brief Set up the simulation with an ODE solver.
95 : * @param[in] model An instance of mio::ssir::Model.
96 : * @param[in] t0 Start time.
97 : * @param[in] dt Initial step size of integration.
98 : */
99 9 : FlowSimulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 0.1)
100 9 : : mio::FlowSimulation<ScalarType, Model>(model, t0, dt)
101 : {
102 9 : auto integrator = std::make_shared<mio::EulerIntegratorCore<ScalarType>>();
103 9 : set_integrator(integrator);
104 18 : }
105 :
106 : /**
107 : * @brief Advance the simulation to tmax.
108 : * tmax must be greater than get_result().get_last_time_point().
109 : * @param[in] tmax Next stopping time of the simulation.
110 : */
111 18 : Eigen::Ref<Vector<>> advance(ScalarType tmax)
112 : {
113 18 : assert(get_flows().get_num_time_points() == get_result().get_num_time_points());
114 18 : auto result = this->get_ode_integrator().advance(
115 : // see the general mio::FlowSimulation for more details on this DerivFunction
116 126 : [this](auto&& flows, auto&& t, auto&& dflows_dt) {
117 18 : const auto& pop_result = get_result();
118 18 : auto& model = get_model();
119 : // compute current population
120 54 : model.get_derivatives(flows - get_flows().get_value(pop_result.get_num_time_points() - 1), m_pop);
121 36 : m_pop += pop_result.get_last_value();
122 : // compute the current change in flows with respect to the current population
123 18 : dflows_dt.setZero();
124 18 : model.step_size = get_dt(); // set the current step size
125 54 : model.get_flows(m_pop, m_pop, t, dflows_dt);
126 18 : },
127 18 : tmax, get_dt(), get_flows());
128 18 : compute_population_results();
129 36 : return result;
130 : }
131 : };
132 :
133 : /**
134 : * @brief Run a FlowSimulation of mio::ssir::Model.
135 : * @param[in] t0 Start time.
136 : * @param[in] tmax End time.
137 : * @param[in] dt Initial step size of integration.
138 : * @param[in] model An instance of mio::ssir::Model.
139 : * @return The simulation result as two TimeSeries. The first describes the compartments at each time point,
140 : * the second gives the corresponding flows that lead from t0 to each time point.
141 : */
142 : inline std::vector<TimeSeries<ScalarType>> simulate_flows(ScalarType t0, ScalarType tmax, ScalarType dt,
143 : Model const& model)
144 : {
145 : model.check_constraints();
146 : FlowSimulation sim(model, t0, dt);
147 : sim.advance(tmax);
148 : return {sim.get_result(), sim.get_flows()};
149 : }
150 :
151 : } // namespace ssir
152 : } // namespace mio
153 :
154 : #endif // MIO_SDE_SIR_SIMULATION_H
|
