Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
3 : *
4 : * Authors: René Schmieding, Julia Bicker
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 :
21 : #ifndef MIO_D_ABM_QUAD_WELL_H
22 : #define MIO_D_ABM_QUAD_WELL_H
23 :
24 : #include "memilio/math/eigen.h"
25 : #include "d_abm/parameters.h"
26 : #include "memilio/utils/random_number_generator.h"
27 : #include "memilio/epidemiology/adoption_rate.h"
28 : #include <string>
29 : #include <vector>
30 :
31 : using Position = Eigen::Vector2d;
32 :
33 : /// @brief Get the index of the well containing the given position.
34 65443 : inline size_t well_index(const Position& p)
35 : {
36 : // 0|1
37 : // -+-
38 : // 2|3
39 65443 : return (p.x() >= 0) + 2 * (p.y() < 0);
40 : }
41 :
42 : /**
43 : * @brief Implementation of diffusive ABM, see dabm::Model.
44 : * This implementation defines a diffusion process for the potential F(x,y) = (x^2 -1)^2+(y^2-1)^2.
45 : * @tparam InfectionState An infection state enum.
46 : */
47 : template <class InfectionState>
48 : class QuadWellModel
49 : {
50 :
51 : public:
52 : using Status = InfectionState;
53 :
54 : /**
55 : * @brief Struct defining an agent for the diffusive ABM.
56 : */
57 : struct Agent {
58 : Position position;
59 : Status status;
60 : };
61 :
62 : /**
63 : * @brief Set up a diffusive ABM using the quadwell potential F(x,y) = (x^2 -1)^2+(y^2-1)^2.
64 : * @param[in] agents A vector of Agent%s representing the population.
65 : * @param[in] rates AdoptionRate%s defining InfectionState adoptions, see mio::AdoptionRate.
66 : * @param[in] contact_radius Contact radius for second-order adoptions.
67 : * @param[in] sigma Noise term for the diffusion process.
68 : * @param[in] non_moving_state InfectionStates that are excluded from movement e.g. Dead.
69 : */
70 49 : QuadWellModel(const std::vector<Agent>& agents, const std::vector<mio::AdoptionRate<Status>>& rates,
71 : double contact_radius = 0.4, double sigma = 0.4, std::vector<Status> non_moving_states = {})
72 49 : : populations(agents)
73 49 : , m_contact_radius(contact_radius)
74 49 : , m_sigma(sigma)
75 49 : , m_non_moving_states(non_moving_states)
76 98 : , m_number_transitions(static_cast<size_t>(Status::Count), Eigen::MatrixXd::Zero(4, 4))
77 : {
78 154 : for (auto& agent : populations) {
79 105 : mio::unused(agent);
80 105 : assert(is_in_domain(agent.position));
81 : }
82 231 : for (auto& r : rates) {
83 182 : m_adoption_rates.emplace(std::forward_as_tuple(r.region, r.from, r.to), r);
84 : }
85 49 : }
86 :
87 : /**
88 : * @brief Perform infection state adoption for an Agent.
89 : * @param[in, out] agent Agent whose infection state is changed.
90 : * @param[in] new_status Agent's new infection state.
91 : */
92 14 : inline static constexpr void adopt(Agent& agent, const Status& new_status)
93 : {
94 14 : agent.status = new_status;
95 14 : }
96 :
97 : /**
98 : * @brief Calculate adoption rate for an Agent.
99 : * @param[in] agent Agent for whom the adoption rate is calculated.
100 : * @param[in] new_status Target infection state of the adoption rate, see mio::AdoptionRate.to.
101 : * @return Value of agent-dependent AdoptionRate.
102 : */
103 37947 : double adoption_rate(const Agent& agent, const Status& new_status) const
104 : {
105 37947 : double rate = 0;
106 : // get the correct adoption rate
107 37947 : const size_t well = well_index(agent.position);
108 37947 : auto map_itr = m_adoption_rates.find({well, agent.status, new_status});
109 37947 : if (map_itr != m_adoption_rates.end()) {
110 2149 : const auto& adoption_rate = map_itr->second;
111 : // calculate the current rate, depending on order
112 2149 : if (adoption_rate.influences.size() == 0) { // first order adoption
113 : // contact independant rate
114 28 : rate = adoption_rate.factor;
115 : }
116 : else { // second order adoption
117 : // accumulate rate per contact with a status in influences
118 2121 : size_t num_contacts = 0;
119 2121 : ScalarType influences = 0;
120 8491 : for (auto& contact : populations) {
121 : // check if contact is indeed a contact
122 6370 : if (is_contact(agent, contact)) {
123 4249 : num_contacts++;
124 12747 : for (size_t i = 0; i < adoption_rate.influences.size(); i++) {
125 8498 : if (contact.status == adoption_rate.influences[i].status) {
126 35 : influences += adoption_rate.influences[i].factor;
127 : }
128 : }
129 : }
130 : }
131 : // rate = factor * "concentration of contacts with status new_status"
132 2121 : if (num_contacts > 0) {
133 2121 : rate = adoption_rate.factor * (influences / num_contacts);
134 : }
135 : }
136 : }
137 : // else: no adoption from agent.status to new_status exist
138 75894 : return rate;
139 : }
140 :
141 : /**
142 : * @brief Perform one integration step of the diffusion process for a given Agent using the Euler-Maruyama method.
143 : * @param[in] dt Step size.
144 : * @param[in] agent Agent to be moved.
145 : */
146 6328 : void move(const double /*t*/, const double dt, Agent& agent)
147 : {
148 6328 : const auto old_well = well_index(agent.position);
149 25312 : if (std::find(m_non_moving_states.begin(), m_non_moving_states.end(), agent.status) ==
150 31633 : m_non_moving_states.end() &&
151 12649 : std::find(m_non_moving_regions.begin(), m_non_moving_regions.end(), old_well) ==
152 12649 : m_non_moving_regions.end()) {
153 18942 : Position p = {mio::DistributionAdapter<std::normal_distribution<double>>::get_instance()(m_rng, 0.0, 1.0),
154 12628 : mio::DistributionAdapter<std::normal_distribution<double>>::get_instance()(m_rng, 0.0, 1.0)};
155 :
156 6314 : agent.position = agent.position - dt * grad_U(agent.position) + (m_sigma * std::sqrt(dt)) * p;
157 6314 : const auto new_well = well_index(agent.position);
158 6314 : if (old_well != new_well) {
159 7 : m_number_transitions[static_cast<size_t>(agent.status)](old_well, new_well)++;
160 : }
161 : }
162 : //else{agent has non-moving status or region}
163 6328 : }
164 :
165 : /**
166 : * @brief Calculate the current system state i.e. the populations for each region and infection state.
167 : * @return Vector containing the number of agents per infection state for each region.
168 : */
169 2107 : Eigen::VectorXd time_point() const
170 : {
171 2107 : Eigen::VectorXd val = Eigen::VectorXd::Zero(4 * static_cast<size_t>(Status::Count));
172 8428 : for (auto& agent : populations) {
173 : // split population into the wells given by grad_U
174 6321 : auto position =
175 6321 : static_cast<size_t>(agent.status) + well_index(agent.position) * static_cast<size_t>(Status::Count);
176 6321 : val[position] += 1;
177 : }
178 4214 : return val;
179 2107 : }
180 :
181 : /**
182 : * @brief Get the number of spatial transitions that happened until the current system state.
183 : * @return Matrix with entries (from_well, to_well) for every infection state.
184 : */
185 : const std::vector<Eigen::MatrixXd>& number_transitions() const
186 : {
187 : return m_number_transitions;
188 : }
189 :
190 7 : std::vector<Eigen::MatrixXd>& number_transitions()
191 : {
192 7 : return m_number_transitions;
193 : }
194 :
195 : /**
196 : * @brief Get AdoptionRate mapping.
197 : * @return Map of AdoptionRates based on their region index and source and target infection state.
198 : */
199 : std::map<std::tuple<mio::regions::Region, Status, Status>, mio::AdoptionRate<Status>>& get_adoption_rates()
200 : {
201 : return m_adoption_rates;
202 : }
203 :
204 : /**
205 : * @brief Set regions where no movement happens.
206 : */
207 7 : void set_non_moving_regions(std::vector<size_t> non_moving_regions)
208 : {
209 7 : m_non_moving_regions = non_moving_regions;
210 7 : }
211 :
212 : /**
213 : * Get the RandomNumberGenerator used by this Model for random events.
214 : * @return The random number generator.
215 : */
216 21 : mio::RandomNumberGenerator& get_rng()
217 : {
218 21 : return m_rng;
219 : }
220 :
221 : std::vector<Agent> populations; ///< Vector containing all Agent%s in the model.
222 :
223 : private:
224 : /**
225 : * @brief Claculate the gradient of the potential at a given position.
226 : * @param[in] x Position at which the gradient is evaluated.
227 : * @return Value of potential gradient.
228 : */
229 6314 : static Position grad_U(const Position x)
230 : {
231 : // U is a quad well potential
232 : // U(x0,x1) = (x0^2 - 1)^2 + (x1^2 - 1)^2
233 12628 : return {4 * x[0] * (x[0] * x[0] - 1), 4 * x[1] * (x[1] * x[1] - 1)};
234 : }
235 :
236 : /**
237 : * @brief Evaluate whether two agents are within their contact radius.
238 : * @param[in] agent Agent whose position specifies the contact area.
239 : * @param[in] contact Potential contact of agent.
240 : * @return Boolean specifying whether are within each others contact radius.
241 : */
242 6370 : bool is_contact(const Agent& agent, const Agent& contact) const
243 : {
244 : // test if contact is in the contact radius and test if agent and contact are different objects
245 16989 : return (agent.position - contact.position).norm() < m_contact_radius && (&agent != &contact) &&
246 16989 : well_index(agent.position) == well_index(contact.position);
247 : }
248 :
249 : /**
250 : * @brief Restrict domain to [-2, 2]^2 where "escaping" is impossible.
251 : * @param[in] p Position to check.
252 : * @return Boolean specifying whether p is in [-2, 2]^2.
253 : */
254 105 : bool is_in_domain(const Position& p) const
255 : {
256 105 : return -2 <= p[0] && p[0] <= 2 && -2 <= p[1] && p[1] <= 2;
257 : }
258 :
259 : std::map<std::tuple<mio::regions::Region, Status, Status>, mio::AdoptionRate<Status>>
260 : m_adoption_rates; ///< Map of AdoptionRates according to their region index and their from -> to infection states.
261 : double m_contact_radius; ///< Agents' interaction radius. Within this radius agents are considered as contacts.
262 : double m_sigma; ///< Noise term of the diffusion process.
263 : std::vector<Status> m_non_moving_states; ///< Infection states within which agents do not change their location.
264 : std::vector<size_t> m_non_moving_regions{}; ///< Regions without movement.
265 : std::vector<Eigen::MatrixXd>
266 : m_number_transitions; ///< Vector that contains for every infection state a matrix with entry (k,l) the number of spatial transitions from Region k to Regionl.
267 : mio::RandomNumberGenerator m_rng; ///< Model's random number generator.
268 : };
269 :
270 : #endif
|
