Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Daniel Abele, Majid Abedi, Elisabeth Kluth, Carlotta Gerstein, Martin J. Kuehn, David Kerkmann, Khoa Nguyen
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 : #include "abm/model.h"
21 : #include "abm/location_id.h"
22 : #include "abm/location_type.h"
23 : #include "abm/intervention_type.h"
24 : #include "abm/person.h"
25 : #include "abm/location.h"
26 : #include "abm/mobility_rules.h"
27 : #include "abm/person_id.h"
28 : #include "memilio/epidemiology/age_group.h"
29 : #include "memilio/utils/logging.h"
30 : #include "memilio/utils/mioomp.h"
31 : #include "memilio/utils/stl_util.h"
32 : #include <cassert>
33 : #include <cstdint>
34 :
35 : namespace mio
36 : {
37 : namespace abm
38 : {
39 :
40 822 : LocationId Model::add_location(LocationType type, uint32_t num_cells)
41 : {
42 822 : LocationId id{static_cast<uint32_t>(m_locations.size())};
43 822 : m_locations.emplace_back(type, id, parameters.get_num_groups(), m_id, num_cells);
44 822 : m_has_locations[size_t(type)] = true;
45 :
46 : // mark caches for rebuild
47 822 : m_is_local_population_cache_valid = false;
48 822 : m_are_exposure_caches_valid = false;
49 822 : m_exposure_caches_need_rebuild = true;
50 :
51 1644 : return id;
52 : }
53 :
54 1045 : PersonId Model::add_person(const LocationId id, AgeGroup age)
55 : {
56 1045 : PersonId person_id = (static_cast<int64_t>(m_id)) << 32 | static_cast<uint32_t>(m_persons.size());
57 2090 : return add_person(Person(m_rng, get_location(id).get_type(), id, m_id, age, person_id));
58 : }
59 :
60 1327 : PersonId Model::add_person(Person&& person)
61 : {
62 1327 : assert(person.get_location() != LocationId::invalid_id() && "Added Person's location must be valid.");
63 1327 : assert(person.get_location() < LocationId((uint32_t)m_locations.size()) &&
64 : "Added Person's location is not in Model.");
65 1327 : assert(person.get_id() != PersonId::invalid_ID() && "Added Person's unique id must be valid.");
66 1327 : assert(person.get_age() < (AgeGroup)parameters.get_num_groups() && "Added Person's AgeGroup is too large.");
67 1327 : person.set_assigned_location(LocationType::Cemetery, m_cemetery_id, m_id);
68 1327 : m_persons.emplace_back(person);
69 1327 : m_activeness_statuses.push_back(true);
70 1327 : auto& new_person = m_persons.back();
71 1327 : if (m_is_local_population_cache_valid) {
72 3 : ++m_local_population_cache[new_person.get_location().get()];
73 : }
74 1327 : return new_person.get_id();
75 : }
76 :
77 774 : void Model::evolve(TimePoint t, TimeSpan dt)
78 : {
79 774 : begin_step(t, dt);
80 774 : log_info("ABM Model interaction.");
81 774 : interaction(t, dt);
82 774 : log_info("ABM Model mobility.");
83 774 : perform_mobility(t, dt);
84 774 : }
85 :
86 1063 : void Model::interaction(TimePoint t, TimeSpan dt)
87 : {
88 1063 : const uint32_t num_persons = static_cast<uint32_t>(m_persons.size());
89 : PRAGMA_OMP(parallel for)
90 3929 : for (uint32_t person_index = 0; person_index < num_persons; ++person_index) {
91 2866 : if (m_activeness_statuses[person_index]) {
92 2853 : assert(m_persons[person_index].get_location_model_id() == m_id &&
93 : "Person is not in this model but still active.");
94 2853 : interact(m_persons[person_index], t, dt);
95 : }
96 : }
97 1063 : }
98 :
99 774 : void Model::perform_mobility(TimePoint t, TimeSpan dt)
100 : {
101 774 : const uint32_t num_persons = static_cast<uint32_t>(m_persons.size());
102 : PRAGMA_OMP(parallel for)
103 3564 : for (uint32_t person_index = 0; person_index < num_persons; ++person_index) {
104 2790 : if (m_activeness_statuses[person_index]) {
105 2790 : Person& person = m_persons[person_index];
106 2790 : auto personal_rng = PersonalRandomNumberGenerator(person);
107 :
108 2790 : auto try_mobility_rule = [&](auto rule) -> bool {
109 : // run mobility rule and check if change of location can actually happen
110 82017 : auto target_type = rule(personal_rng, person, t, dt, parameters);
111 26676 : if (person.get_assigned_location_model_id(target_type) == m_id) {
112 35568 : const Location& target_location = get_location(find_location(target_type, person));
113 26676 : const LocationId current_location = person.get_location();
114 :
115 : // the Person cannot move if they do not wear mask as required at targeted location
116 8928 : if (target_location.is_mask_required() &&
117 72 : !person.is_compliant(personal_rng, InterventionType::Mask)) {
118 9 : return false;
119 : }
120 : // the Person cannot move if the capacity of targeted Location is reached
121 9045 : if (target_location.get_id() == current_location ||
122 324 : get_number_persons(target_location.get_id()) >= target_location.get_capacity().persons) {
123 8730 : return false;
124 : }
125 : // the Person cannot move if the performed TestingStrategy is positive
126 612 : if (!m_testing_strategy.run_strategy(personal_rng, person, target_location, t)) {
127 18 : return false;
128 : }
129 : // update worn mask to target location's requirements
130 135 : if (target_location.is_mask_required()) {
131 : // if the current MaskProtection level is lower than required, the Person changes mask
132 45 : if (parameters.get<MaskProtection>()[person.get_mask().get_type()] <
133 18 : parameters.get<MaskProtection>()[target_location.get_required_mask()]) {
134 36 : person.set_mask(target_location.get_required_mask(), t);
135 : }
136 : }
137 : else {
138 504 : person.set_mask(MaskType::None, t);
139 : }
140 : // all requirements are met, move to target location
141 405 : change_location(person, target_location.get_id());
142 135 : return true;
143 : }
144 0 : return false;
145 2790 : };
146 :
147 : // run mobility rules one after the other if the corresponding location type exists
148 : // shortcutting of bool operators ensures the rules stop after the first rule is applied
149 2790 : if (m_use_mobility_rules) {
150 10980 : (has_locations({LocationType::Cemetery}) && try_mobility_rule(&get_buried)) ||
151 5454 : (has_locations({LocationType::Home}) && try_mobility_rule(&return_home_when_recovered)) ||
152 5436 : (has_locations({LocationType::Hospital}) && try_mobility_rule(&go_to_hospital)) ||
153 5418 : (has_locations({LocationType::ICU}) && try_mobility_rule(&go_to_icu)) ||
154 5418 : (has_locations({LocationType::School, LocationType::Home}) && try_mobility_rule(&go_to_school)) ||
155 5391 : (has_locations({LocationType::Work, LocationType::Home}) && try_mobility_rule(&go_to_work)) ||
156 5373 : (has_locations({LocationType::BasicsShop, LocationType::Home}) && try_mobility_rule(&go_to_shop)) ||
157 5364 : (has_locations({LocationType::SocialEvent, LocationType::Home}) &&
158 8964 : try_mobility_rule(&go_to_event)) ||
159 5355 : (has_locations({LocationType::Home}) && try_mobility_rule(&go_to_quarantine));
160 : }
161 : else {
162 : // no daily routine mobility, just infection related
163 180 : (has_locations({LocationType::Cemetery}) && try_mobility_rule(&get_buried)) ||
164 90 : (has_locations({LocationType::Home}) && try_mobility_rule(&return_home_when_recovered)) ||
165 90 : (has_locations({LocationType::Hospital}) && try_mobility_rule(&go_to_hospital)) ||
166 225 : (has_locations({LocationType::ICU}) && try_mobility_rule(&go_to_icu)) ||
167 90 : (has_locations({LocationType::Home}) && try_mobility_rule(&go_to_quarantine));
168 : }
169 : }
170 : }
171 :
172 : // check if a person makes a trip
173 774 : bool weekend = t.is_weekend();
174 774 : size_t num_trips = m_trip_list.num_trips(weekend);
175 :
176 2205 : for (; m_trip_list.get_current_index() < num_trips &&
177 1521 : m_trip_list.get_next_trip_time(weekend).seconds() < (t + dt).time_since_midnight().seconds();
178 189 : m_trip_list.increase_index()) {
179 189 : auto& trip = m_trip_list.get_next_trip(weekend);
180 189 : auto& person = get_person(trip.person_id);
181 189 : auto personal_rng = PersonalRandomNumberGenerator(person);
182 : // skip the trip if the person is in quarantine or is dead
183 189 : if (person.is_in_quarantine(t, parameters) || person.get_infection_state(t) == InfectionState::Dead) {
184 0 : continue;
185 : }
186 189 : auto& target_location = get_location(trip.destination);
187 : // skip the trip if the Person wears mask as required at targeted location
188 189 : if (target_location.is_mask_required() && !person.is_compliant(personal_rng, InterventionType::Mask)) {
189 9 : continue;
190 : }
191 : // skip the trip if the performed TestingStrategy is positive
192 180 : if (!m_testing_strategy.run_strategy(personal_rng, person, target_location, t)) {
193 18 : continue;
194 : }
195 : // all requirements are met, move to target location
196 162 : change_location(person, target_location.get_id(), trip.trip_mode);
197 : // update worn mask to target location's requirements
198 162 : if (target_location.is_mask_required()) {
199 : // if the current MaskProtection level is lower than required, the Person changes mask
200 36 : if (parameters.get<MaskProtection>()[person.get_mask().get_type()] <
201 27 : parameters.get<MaskProtection>()[target_location.get_required_mask()]) {
202 9 : person.set_mask(target_location.get_required_mask(), t);
203 : }
204 : }
205 : else {
206 153 : person.set_mask(MaskType::None, t);
207 : }
208 : }
209 774 : if (((t).days() < std::floor((t + dt).days()))) {
210 63 : m_trip_list.reset_index();
211 : }
212 774 : }
213 :
214 100 : void Model::build_compute_local_population_cache() const
215 : {
216 : PRAGMA_OMP(single)
217 : {
218 100 : const size_t num_locations = m_locations.size();
219 100 : const size_t num_persons = m_persons.size();
220 100 : m_local_population_cache.resize(num_locations);
221 : PRAGMA_OMP(taskloop)
222 532 : for (size_t i = 0; i < num_locations; i++) {
223 432 : m_local_population_cache[i] = 0;
224 : } // implicit taskloop barrier
225 : PRAGMA_OMP(taskloop)
226 407 : for (size_t i = 0; i < num_persons; i++) {
227 307 : if (m_activeness_statuses[i]) {
228 301 : assert(m_persons[i].get_location_model_id() == m_id && "Person is not in this model but still active.");
229 301 : ++m_local_population_cache[m_persons[i].get_location().get()];
230 : }
231 : } // implicit taskloop barrier
232 : } // implicit single barrier
233 100 : }
234 :
235 97 : void Model::build_exposure_caches()
236 : {
237 : PRAGMA_OMP(single)
238 : {
239 97 : const size_t num_locations = m_locations.size();
240 97 : m_air_exposure_rates_cache.resize(num_locations);
241 97 : m_contact_exposure_rates_cache.resize(num_locations);
242 : PRAGMA_OMP(taskloop)
243 519 : for (size_t i = 0; i < num_locations; i++) {
244 422 : m_air_exposure_rates_cache[i].resize({CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count});
245 1266 : m_contact_exposure_rates_cache[i].resize({CellIndex(m_locations[i].get_cells().size()), VirusVariant::Count,
246 844 : AgeGroup(parameters.get_num_groups())});
247 : } // implicit taskloop barrier
248 97 : m_are_exposure_caches_valid = false;
249 97 : m_exposure_caches_need_rebuild = false;
250 : } // implicit single barrier
251 97 : }
252 :
253 1063 : void Model::compute_exposure_caches(TimePoint t, TimeSpan dt)
254 : {
255 : PRAGMA_OMP(single)
256 : {
257 : // if cache shape was changed (e.g. by add_location), rebuild it
258 1063 : if (m_exposure_caches_need_rebuild) {
259 97 : build_exposure_caches();
260 : }
261 : // use these const values to help omp recognize that the for loops are bounded
262 1063 : const auto num_locations = m_locations.size();
263 1063 : const auto num_persons = m_persons.size();
264 :
265 : // 1) reset all cached values
266 : // Note: we cannot easily reuse values, as they are time dependant (get_infection_state)
267 : PRAGMA_OMP(taskloop)
268 5214 : for (size_t i = 0; i < num_locations; ++i) {
269 4151 : const auto index = i;
270 4151 : auto& local_air_exposure = m_air_exposure_rates_cache[index];
271 4151 : std::for_each(local_air_exposure.begin(), local_air_exposure.end(), [](auto& r) {
272 4151 : r = 0.0;
273 4151 : });
274 4151 : auto& local_contact_exposure = m_contact_exposure_rates_cache[index];
275 4151 : std::for_each(local_contact_exposure.begin(), local_contact_exposure.end(), [](auto& r) {
276 23126 : r = 0.0;
277 23126 : });
278 : } // implicit taskloop barrier
279 : // here is an implicit (and needed) barrier from parallel for
280 :
281 : // 2) add all contributions from each person
282 : PRAGMA_OMP(taskloop)
283 3929 : for (size_t i = 0; i < num_persons; ++i) {
284 2866 : const Person& person = m_persons[i];
285 2866 : const auto location = person.get_location().get();
286 2866 : if (m_activeness_statuses[i]) {
287 2853 : assert(m_persons[i].get_location_model_id() == m_id && "Person is not in this model but still active.");
288 5706 : mio::abm::add_exposure_contribution(m_air_exposure_rates_cache[location],
289 2853 : m_contact_exposure_rates_cache[location], person,
290 2853 : get_location(person.get_location()), t, dt);
291 : }
292 : } // implicit taskloop barrier
293 : } // implicit single barrier
294 1063 : }
295 :
296 1063 : void Model::begin_step(TimePoint t, TimeSpan dt)
297 : {
298 1063 : m_testing_strategy.update_activity_status(t);
299 :
300 1063 : if (!m_is_local_population_cache_valid) {
301 100 : build_compute_local_population_cache();
302 100 : m_is_local_population_cache_valid = true;
303 : }
304 1063 : compute_exposure_caches(t, dt);
305 1063 : m_are_exposure_caches_valid = true;
306 1063 : }
307 :
308 720 : auto Model::get_locations() const -> Range<std::pair<ConstLocationIterator, ConstLocationIterator>>
309 : {
310 720 : return std::make_pair(m_locations.cbegin(), m_locations.cend());
311 : }
312 63 : auto Model::get_locations() -> Range<std::pair<LocationIterator, LocationIterator>>
313 : {
314 63 : return std::make_pair(m_locations.begin(), m_locations.end());
315 : }
316 :
317 486 : auto Model::get_persons() const -> Range<std::pair<ConstPersonIterator, ConstPersonIterator>>
318 : {
319 486 : return std::make_pair(m_persons.cbegin(), m_persons.cend());
320 : }
321 :
322 379 : auto Model::get_persons() -> Range<std::pair<PersonIterator, PersonIterator>>
323 : {
324 379 : return std::make_pair(m_persons.begin(), m_persons.end());
325 : }
326 :
327 0 : auto Model::get_activeness_statuses() const -> Range<std::pair<ConstActivenessIterator, ConstActivenessIterator>>
328 : {
329 0 : return std::make_pair(m_activeness_statuses.cbegin(), m_activeness_statuses.cend());
330 : }
331 :
332 5 : auto Model::get_activeness_statuses() -> Range<std::pair<ActivenessIterator, ActivenessIterator>>
333 : {
334 10 : return std::make_pair(m_activeness_statuses.begin(), m_activeness_statuses.end());
335 : }
336 :
337 81 : LocationId Model::find_location(LocationType type, const PersonId person) const
338 : {
339 162 : return find_location(type, get_person(person));
340 : }
341 :
342 9 : size_t Model::get_subpopulation_combined(TimePoint t, InfectionState s) const
343 : {
344 9 : return std::accumulate(m_locations.begin(), m_locations.end(), (size_t)0,
345 90 : [t, s, this](size_t running_sum, const Location& loc) {
346 90 : return running_sum + get_subpopulation(loc.get_id(), t, s);
347 9 : });
348 : }
349 :
350 18 : size_t Model::get_subpopulation_combined_per_location_type(TimePoint t, InfectionState s, LocationType type) const
351 : {
352 36 : return std::accumulate(
353 216 : m_locations.begin(), m_locations.end(), (size_t)0, [t, s, type, this](size_t running_sum, const Location& loc) {
354 378 : return loc.get_type() == type ? running_sum + get_subpopulation(loc.get_id(), t, s) : running_sum;
355 18 : });
356 : }
357 :
358 47 : TripList& Model::get_trip_list()
359 : {
360 47 : return m_trip_list;
361 : }
362 :
363 0 : const TripList& Model::get_trip_list() const
364 : {
365 0 : return m_trip_list;
366 : }
367 :
368 9 : void Model::use_mobility_rules(bool param)
369 : {
370 9 : m_use_mobility_rules = param;
371 9 : }
372 :
373 0 : bool Model::use_mobility_rules() const
374 : {
375 0 : return m_use_mobility_rules;
376 : }
377 :
378 81 : TestingStrategy& Model::get_testing_strategy()
379 : {
380 81 : return m_testing_strategy;
381 : }
382 :
383 0 : const TestingStrategy& Model::get_testing_strategy() const
384 : {
385 0 : return m_testing_strategy;
386 : }
387 :
388 : } // namespace abm
389 : } // namespace mio
|
