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