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_SEIRVV_SIMULATION_H
21 : #define MIO_SDE_SEIRVV_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_seirvv/model.h"
27 :
28 : namespace mio
29 : {
30 : namespace sseirvv
31 : {
32 :
33 : /// @brief A specialized Simulation for mio::sseirvv::Model.
34 : class Simulation : public mio::Simulation<ScalarType, Model>
35 : {
36 : protected:
37 : using mio::Simulation<ScalarType, Model>::set_integrator;
38 : public:
39 : /**
40 : * @brief Set up the simulation with an SDE solver.
41 : * @param[in] model An instance of mio::sseirvv::Model.
42 : * @param[in] t0 Start time.
43 : * @param[in] dt Initial step size of integration.
44 : */
45 9 : Simulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 0.1)
46 9 : : mio::Simulation<ScalarType, Model>(model, t0, dt)
47 : {
48 9 : auto integrator = std::make_shared<mio::EulerIntegratorCore<ScalarType>>();
49 9 : set_integrator(integrator);
50 18 : }
51 :
52 : using mio::Simulation<ScalarType, Model>::Simulation;
53 : /**
54 : * @brief Advance simulation to tmax.
55 : * tmax must be greater than get_result().get_last_time_point().
56 : * @param tmax Next stopping point of simulation.
57 : */
58 18 : Eigen::Ref<Eigen::VectorXd> advance(ScalarType tmax)
59 : {
60 18 : return get_ode_integrator().advance(
61 18 : [this](auto&& y, auto&& t, auto&& dydt) {
62 18 : get_model().step_size = get_dt();
63 18 : get_model().get_derivatives(y, y, t, dydt);
64 18 : },
65 18 : tmax, get_dt(), get_result());
66 : }
67 : };
68 :
69 : /**
70 : * @brief Run a Simulation of a mio::sseirvv::Model.
71 : * @param[in] t0 Start time.
72 : * @param[in] tmax End time.
73 : * @param[in] dt Initial step size of integration.
74 : * @param[in] model An instance of mio::sseirvv::Model.
75 : * @return A TimeSeries to represent the final simulation result
76 : */
77 0 : inline TimeSeries<ScalarType> simulate(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& model)
78 : {
79 0 : model.check_constraints();
80 0 : Simulation sim(model, t0, dt);
81 0 : sim.advance(tmax);
82 0 : return sim.get_result();
83 0 : }
84 :
85 : /// @brief A specialized FlowSimulation for mio::sseirvv::Model.
86 : class FlowSimulation : public mio::FlowSimulation<ScalarType, Model>
87 : {
88 : protected:
89 : using mio::FlowSimulation<ScalarType, Model>::set_integrator;
90 :
91 : public:
92 : /**
93 : * @brief Set up the simulation with an ODE solver.
94 : * @param[in] model An instance of mio::sseirvv::Model.
95 : * @param[in] t0 Start time.
96 : * @param[in] dt Initial step size of integration.
97 : */
98 9 : FlowSimulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 0.1)
99 9 : : mio::FlowSimulation<ScalarType, Model>(model, t0, dt)
100 : {
101 9 : auto integrator = std::make_shared<mio::EulerIntegratorCore<ScalarType>>();
102 9 : set_integrator(integrator);
103 18 : }
104 :
105 : /**
106 : * @brief Advance the simulation to tmax.
107 : * tmax must be greater than get_result().get_last_time_point().
108 : * @param[in] tmax Next stopping time of the simulation.
109 : */
110 18 : Eigen::Ref<Eigen::VectorXd> advance(ScalarType tmax)
111 : {
112 18 : assert(get_flows().get_num_time_points() == get_result().get_num_time_points());
113 18 : auto result = this->get_ode_integrator().advance(
114 : // See the general mio::FlowSimulation for more details on this DerivFunction.
115 126 : [this](auto&& flows, auto&& t, auto&& dflows_dt) {
116 18 : const auto& pop_result = get_result();
117 18 : auto& model = get_model();
118 : // Compute current population.
119 54 : model.get_derivatives(flows - get_flows().get_value(pop_result.get_num_time_points() - 1), m_pop);
120 36 : m_pop += pop_result.get_last_value();
121 : // Compute the current change in flows with respect to the current population.
122 18 : dflows_dt.setZero();
123 18 : model.step_size = get_dt(); // set the current step size
124 54 : model.get_flows(m_pop, m_pop, t, dflows_dt);
125 18 : },
126 18 : tmax, get_dt(), get_flows());
127 18 : compute_population_results();
128 36 : return result;
129 : }
130 : };
131 :
132 : /**
133 : * @brief Run a FlowSimulation of mio::sseirvv::Model.
134 : * @param[in] t0 Start time.
135 : * @param[in] tmax End time.
136 : * @param[in] dt Initial step size of integration.
137 : * @param[in] model An instance of mio::sseirvv::Model.
138 : * @return The simulation result as two TimeSeries. The first describes the compartments at each time point,
139 : * the second gives the corresponding flows that lead from t0 to each time point.
140 : */
141 : inline std::vector<TimeSeries<ScalarType>> simulate_flows(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& model)
142 : {
143 : model.check_constraints();
144 : FlowSimulation sim(model, t0, dt);
145 : sim.advance(tmax);
146 : return {sim.get_result(), sim.get_flows()};
147 : }
148 :
149 : } // namespace sseirvv
150 : } // namespace mio
151 :
152 : #endif // MIO_SDE_SEIRVV_SIMULATION_H
|
