Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 MEmilio
3 : *
4 : * Authors: Daniel Abele, Majid Abedi, Elisabeth Kluth, David Kerkmann, Sascha Korf, Martin J. Kuehn, 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 : #ifndef MIO_ABM_MODEL_H
21 : #define MIO_ABM_MODEL_H
22 :
23 : #include "abm/model_functions.h"
24 : #include "abm/location_type.h"
25 : #include "abm/mobility_data.h"
26 : #include "abm/parameters.h"
27 : #include "abm/location.h"
28 : #include "abm/person.h"
29 : #include "abm/person_id.h"
30 : #include "abm/time.h"
31 : #include "abm/trip_list.h"
32 : #include "abm/random_events.h"
33 : #include "abm/testing_strategy.h"
34 : #include "memilio/epidemiology/age_group.h"
35 : #include "memilio/utils/random_number_generator.h"
36 : #include "memilio/utils/stl_util.h"
37 :
38 : #include <bitset>
39 : #include <vector>
40 :
41 : namespace mio
42 : {
43 : namespace abm
44 : {
45 :
46 : /**
47 : * @brief The Model of the Simulation.
48 : * It consists of Location%s and Person%s (Agents).
49 : */
50 : class Model
51 : {
52 : public:
53 : using LocationIterator = std::vector<Location>::iterator;
54 : using ConstLocationIterator = std::vector<Location>::const_iterator;
55 : using PersonIterator = std::vector<Person>::iterator;
56 : using ConstPersonIterator = std::vector<Person>::const_iterator;
57 :
58 : /**
59 : * @brief Create a Model.
60 : * @param[in] num_agegroups The number of AgeGroup%s in the simulated Model. Must be less than MAX_NUM_AGE_GROUPS.
61 : */
62 199 : Model(size_t num_agegroups)
63 199 : : parameters(num_agegroups)
64 199 : , m_trip_list()
65 199 : , m_use_mobility_rules(true)
66 398 : , m_cemetery_id(add_location(LocationType::Cemetery))
67 : {
68 199 : assert(num_agegroups < MAX_NUM_AGE_GROUPS && "MAX_NUM_AGE_GROUPS exceeded.");
69 199 : }
70 :
71 : /**
72 : * @brief Create a Model.
73 : * @param[in] params Initial simulation parameters.
74 : */
75 9 : Model(const Parameters& params)
76 9 : : parameters(params.get_num_groups())
77 9 : , m_trip_list()
78 9 : , m_use_mobility_rules(true)
79 18 : , m_cemetery_id(add_location(LocationType::Cemetery))
80 : {
81 9 : parameters = params;
82 9 : }
83 :
84 140 : Model(const Model& other)
85 140 : : parameters(other.parameters)
86 140 : , m_local_population_cache()
87 140 : , m_air_exposure_rates_cache()
88 140 : , m_contact_exposure_rates_cache()
89 140 : , m_is_local_population_cache_valid(false)
90 140 : , m_are_exposure_caches_valid(false)
91 140 : , m_exposure_caches_need_rebuild(true)
92 140 : , m_persons(other.m_persons)
93 140 : , m_locations(other.m_locations)
94 140 : , m_has_locations(other.m_has_locations)
95 140 : , m_testing_strategy(other.m_testing_strategy)
96 140 : , m_trip_list(other.m_trip_list)
97 140 : , m_use_mobility_rules(other.m_use_mobility_rules)
98 140 : , m_mobility_rules(other.m_mobility_rules)
99 140 : , m_cemetery_id(other.m_cemetery_id)
100 140 : , m_rng(other.m_rng)
101 : {
102 140 : }
103 : Model& operator=(const Model&) = default;
104 72 : Model(Model&&) = default;
105 : Model& operator=(Model&&) = default;
106 :
107 : /**
108 : * serialize this.
109 : * @see mio::serialize
110 : */
111 : template <class IOContext>
112 9 : void serialize(IOContext& io) const
113 : {
114 27 : auto obj = io.create_object("Model");
115 9 : obj.add_element("parameters", parameters);
116 : // skip caches, they are rebuild by the deserialized model
117 9 : obj.add_list("persons", get_persons().begin(), get_persons().end());
118 9 : obj.add_list("locations", get_locations().begin(), get_locations().end());
119 9 : obj.add_element("location_types", m_has_locations.to_ulong());
120 9 : obj.add_element("testing_strategy", m_testing_strategy);
121 9 : obj.add_element("trip_list", m_trip_list);
122 9 : obj.add_element("use_mobility_rules", m_use_mobility_rules);
123 9 : obj.add_element("cemetery_id", m_cemetery_id);
124 9 : obj.add_element("rng", m_rng);
125 18 : }
126 :
127 : /**
128 : * deserialize an object of this class.
129 : * @see mio::deserialize
130 : */
131 : template <class IOContext>
132 9 : static IOResult<Model> deserialize(IOContext& io)
133 : {
134 27 : auto obj = io.expect_object("Model");
135 27 : auto params = obj.expect_element("parameters", Tag<Parameters>{});
136 27 : auto persons = obj.expect_list("persons", Tag<Person>{});
137 27 : auto locations = obj.expect_list("locations", Tag<Location>{});
138 27 : auto location_types = obj.expect_element("location_types", Tag<unsigned long>{});
139 27 : auto trip_list = obj.expect_element("trip_list", Tag<TripList>{});
140 27 : auto use_mobility_rules = obj.expect_element("use_mobility_rules", Tag<bool>{});
141 27 : auto cemetery_id = obj.expect_element("cemetery_id", Tag<LocationId>{});
142 27 : auto rng = obj.expect_element("rng", Tag<RandomNumberGenerator>{});
143 : return apply(
144 : io,
145 9 : [](auto&& params_, auto&& persons_, auto&& locations_, auto&& location_types_, auto&& trip_list_,
146 : auto&& use_mobility_rules_, auto&& cemetery_id_, auto&& rng_) {
147 9 : Model model{params_};
148 9 : model.m_persons.assign(persons_.cbegin(), persons_.cend());
149 9 : model.m_locations.assign(locations_.cbegin(), locations_.cend());
150 9 : model.m_has_locations = location_types_;
151 9 : model.m_trip_list = trip_list_;
152 9 : model.m_use_mobility_rules = use_mobility_rules_;
153 9 : model.m_cemetery_id = cemetery_id_;
154 9 : model.m_rng = rng_;
155 18 : return model;
156 9 : },
157 18 : params, persons, locations, location_types, trip_list, use_mobility_rules, cemetery_id, rng);
158 9 : }
159 :
160 : /**
161 : * @brief Prepare the Model for the next Simulation step.
162 : * @param[in] t Current time.
163 : * @param[in] dt Length of the time step.
164 : */
165 : void begin_step(TimePoint t, TimeSpan dt);
166 :
167 : /**
168 : * @brief Evolve the Model one time step.
169 : * @param[in] t Current time.
170 : * @param[in] dt Length of the time step.
171 : */
172 : void evolve(TimePoint t, TimeSpan dt);
173 :
174 : /**
175 : * @brief Add a Location to the Model.
176 : * @param[in] type Type of Location to add.
177 : * @param[in] num_cells [Default: 1] Number of Cell%s that the Location is divided into.
178 : * @return ID of the newly created Location.
179 : */
180 : LocationId add_location(LocationType type, uint32_t num_cells = 1);
181 :
182 : /**
183 : * @brief Add a Person to the Model.
184 : * @param[in] id The LocationID of the initial Location of the Person.
185 : * @param[in] age AgeGroup of the person.
186 : * @return ID of the newly created Person.
187 : */
188 : PersonId add_person(const LocationId id, AgeGroup age);
189 :
190 : /**
191 : * @brief Adds a copy of a given Person to the Model.
192 : * @param[in] person The Person to copy from.
193 : * @return ID of the newly created Person.
194 : */
195 : PersonId add_person(Person&& person);
196 :
197 : /**
198 : * @brief Get a range of all Location%s in the Model.
199 : * @return A range of all Location%s.
200 : * @{
201 : */
202 : Range<std::pair<ConstLocationIterator, ConstLocationIterator>> get_locations() const;
203 : Range<std::pair<LocationIterator, LocationIterator>> get_locations();
204 : /** @} */
205 :
206 : /**
207 : * @brief Get a range of all Person%s in the Model.
208 : * @return A range of all Person%s.
209 : * @{
210 : */
211 : Range<std::pair<ConstPersonIterator, ConstPersonIterator>> get_persons() const;
212 : Range<std::pair<PersonIterator, PersonIterator>> get_persons();
213 : /** @} */
214 :
215 : /**
216 : * @brief Find an assigned Location of a Person.
217 : * @param[in] type The #LocationType that specifies the assigned Location.
218 : * @param[in] person PersonId of the Person.
219 : * @return ID of the Location of LocationType type assigend to person.
220 : */
221 : LocationId find_location(LocationType type, const PersonId person) const;
222 :
223 : /**
224 : * @brief Assign a Location to a Person.
225 : * A Person can have at most one assigned Location of a certain LocationType.
226 : * Assigning another Location of an already assigned LocationType will replace the prior assignment.
227 : * @param[in] person The PersonId of the person this location will be assigned to.
228 : * @param[in] location The LocationId of the Location.
229 : */
230 1098 : void assign_location(PersonId person, LocationId location)
231 : {
232 1098 : get_person(person).set_assigned_location(get_location(location).get_type(), location);
233 1098 : }
234 :
235 : /**
236 : * @brief Get the number of Persons in one #InfectionState at all Location%s.
237 : * @param[in] t Specified #TimePoint.
238 : * @param[in] s Specified #InfectionState.
239 : */
240 : size_t get_subpopulation_combined(TimePoint t, InfectionState s) const;
241 :
242 : /**
243 : * @brief Get the number of Persons in one #InfectionState at all Location%s of a type.
244 : * @param[in] t Specified #TimePoint.
245 : * @param[in] s Specified #InfectionState.
246 : * @param[in] type Specified #LocationType.
247 : */
248 : size_t get_subpopulation_combined_per_location_type(TimePoint t, InfectionState s, LocationType type) const;
249 :
250 : /**
251 : * @brief Get the mobility data.
252 : * @return Reference to the list of Trip%s that the Person%s make.
253 : */
254 : TripList& get_trip_list();
255 :
256 : const TripList& get_trip_list() const;
257 :
258 : /**
259 : * @brief Decide if mobility rules (like go to school/work) are used or not;
260 : * The mobility rules regarding hospitalization/ICU/quarantine are always used.
261 : * @param[in] param If true uses the mobility rules for changing location to school/work etc., else only the rules
262 : * regarding hospitalization/ICU/quarantine.
263 : */
264 : void use_mobility_rules(bool param);
265 : bool use_mobility_rules() const;
266 :
267 : /**
268 : * @brief Check if at least one Location with a specified LocationType exists.
269 : * @return True if there is at least one Location of LocationType `type`. False otherwise.
270 : */
271 28899 : bool has_location(LocationType type) const
272 : {
273 28899 : return m_has_locations[size_t(type)];
274 : }
275 :
276 : /**
277 : * @brief Check if at least one Location of every specified LocationType exists.
278 : * @tparam C A type of container of LocationType.
279 : * @param location_types A container of LocationType%s.
280 : * @return True if there is at least one Location of every LocationType in `location_types`. False otherwise.
281 : */
282 : template <class C = std::initializer_list<LocationType>>
283 24156 : bool has_locations(const C& location_types) const
284 : {
285 81954 : return std::all_of(location_types.begin(), location_types.end(), [&](auto loc) {
286 57798 : return has_location(loc);
287 24156 : });
288 : }
289 :
290 : /**
291 : * @brief Get the TestingStrategy.
292 : * @return Reference to the list of TestingScheme%s that are checked for testing.
293 : */
294 : TestingStrategy& get_testing_strategy();
295 :
296 : const TestingStrategy& get_testing_strategy() const;
297 :
298 : /**
299 : * @brief The simulation parameters of the Model.
300 : */
301 : Parameters parameters;
302 :
303 : /**
304 : * Get the RandomNumberGenerator used by this Model for random events.
305 : * Persons use their own generators with the same key as the global one.
306 : * @return The random number generator.
307 : */
308 936 : RandomNumberGenerator& get_rng()
309 : {
310 936 : return m_rng;
311 : }
312 :
313 : /**
314 : * @brief Add a TestingScheme to the set of schemes that are checked for testing at all Locations that have
315 : * the LocationType.
316 : * @param[in] loc_type LocationId key for TestingScheme to be added.
317 : * @param[in] scheme TestingScheme to be added.
318 : */
319 : void add_testing_scheme(const LocationType& loc_type, const TestingScheme& scheme);
320 :
321 : /**
322 : * @brief Remove a TestingScheme from the set of schemes that are checked for testing at all Locations that have
323 : * the LocationType.
324 : * @param[in] loc_type LocationId key for TestingScheme to be added.
325 : * @param[in] scheme TestingScheme to be added.
326 : */
327 : void remove_testing_scheme(const LocationType& loc_type, const TestingScheme& scheme);
328 :
329 : /**
330 : * @brief Get a reference to a Person from this Model.
331 : * @param[in] id A person's PersonId.
332 : * @return A reference to the Person.
333 : * @{
334 : */
335 18828 : Person& get_person(PersonId id)
336 : {
337 18828 : assert(id.get() < m_persons.size() && "Given PersonId is not in this Model.");
338 18828 : return m_persons[id.get()];
339 : }
340 :
341 8910 : const Person& get_person(PersonId id) const
342 : {
343 8910 : assert(id.get() < m_persons.size() && "Given PersonId is not in this Model.");
344 8910 : return m_persons[id.get()];
345 : }
346 : /** @} */
347 :
348 : /**
349 : * @brief Get the number of Person%s of a particular #InfectionState for all Cell%s.
350 : * @param[in] location A LocationId from the Model.
351 : * @param[in] t TimePoint of querry.
352 : * @param[in] state #InfectionState of interest.
353 : * @return Amount of Person%s of the #InfectionState in all Cell%s of the Location.
354 : */
355 27315 : size_t get_subpopulation(LocationId location, TimePoint t, InfectionState state) const
356 : {
357 27315 : return std::count_if(m_persons.begin(), m_persons.end(), [&](auto&& p) {
358 186516 : return p.get_location() == location && p.get_infection_state(t) == state;
359 27315 : });
360 : }
361 :
362 : /**
363 : * @brief Get the total number of Person%s at the Location.
364 : * @param[in] location A LocationId from the Model.
365 : * @return Number of Person%s in the location.
366 : */
367 270 : size_t get_number_persons(LocationId location) const
368 : {
369 270 : if (!m_is_local_population_cache_valid) {
370 0 : build_compute_local_population_cache();
371 : }
372 270 : return m_local_population_cache[location.get()];
373 : }
374 :
375 : // Change the Location of a Person. this requires that Location is part of this Model.
376 : /**
377 : * @brief Let a Person change to another Location.
378 : * @param[in] person PersonId of a Person from this Model.
379 : * @param[in] destination LocationId of the Location in this Model, which the Person should change to.
380 : * @param[in] mode The transport mode the person uses to change the Location.
381 : * @param[in] cells The cells within the destination the person should be in.
382 : */
383 315 : inline void change_location(PersonId person, LocationId destination, TransportMode mode = TransportMode::Unknown,
384 : const std::vector<uint32_t>& cells = {0})
385 : {
386 315 : LocationId origin = get_location(person).get_id();
387 : const bool has_changed_location =
388 315 : mio::abm::change_location(get_person(person), get_location(destination), mode, cells);
389 : // if the person has changed location, invalidate exposure caches but keep population caches valid
390 315 : if (has_changed_location) {
391 315 : m_are_exposure_caches_valid = false;
392 315 : if (m_is_local_population_cache_valid) {
393 315 : --m_local_population_cache[origin.get()];
394 315 : ++m_local_population_cache[destination.get()];
395 : }
396 : }
397 315 : }
398 :
399 : /**
400 : * @brief Let a person interact with the population at its current location.
401 : * @param[in] person PersonId of a person from this Model.
402 : * @param[in] t Time step of the simulation.
403 : * @param[in] dt Step size of the simulation.
404 : */
405 2772 : inline void interact(PersonId person, TimePoint t, TimeSpan dt)
406 : {
407 2772 : if (!m_are_exposure_caches_valid) {
408 : // checking caches is only needed for external calls
409 : // during simulation (i.e. in evolve()), the caches are computed in begin_step
410 0 : compute_exposure_caches(t, dt);
411 0 : m_are_exposure_caches_valid = true;
412 : }
413 2772 : auto personal_rng = PersonalRandomNumberGenerator(m_rng, get_person(person));
414 8316 : mio::abm::interact(personal_rng, get_person(person), get_location(person),
415 5544 : m_air_exposure_rates_cache[get_location(person).get_id().get()],
416 8316 : m_contact_exposure_rates_cache[get_location(person).get_id().get()], t, dt, parameters);
417 2772 : }
418 :
419 : /**
420 : * @brief Get a reference to a location in this Model.
421 : * @param[in] id LocationId of the Location.
422 : * @return Reference to the Location.
423 : * @{
424 : */
425 : const Location& get_location(LocationId id) const
426 : {
427 : assert(id != LocationId::invalid_id() && "Given LocationId must be valid.");
428 : assert(id < LocationId((uint32_t)m_locations.size()) && "Given LocationId is not in this Model.");
429 : return m_locations[id.get()];
430 : }
431 :
432 23391 : Location& get_location(LocationId id)
433 : {
434 23391 : assert(id != LocationId::invalid_id() && "Given LocationId must be valid.");
435 23391 : assert(id < LocationId((uint32_t)m_locations.size()) && "Given LocationId is not in this Model.");
436 23391 : return m_locations[id.get()];
437 : }
438 : /** @} */
439 :
440 : /**
441 : * @brief Get a reference to the location of a person.
442 : * @param[in] id PersonId of a person.
443 : * @return Reference to the Location.
444 : * @{
445 : */
446 11439 : inline Location& get_location(PersonId id)
447 : {
448 11439 : return get_location(get_person(id).get_location());
449 : }
450 :
451 : inline const Location& get_location(PersonId id) const
452 : {
453 : return get_location(get_person(id).get_location());
454 : }
455 : /** @} */
456 :
457 : private:
458 : /**
459 : * @brief Person%s interact at their Location and may become infected.
460 : * @param[in] t The current TimePoint.
461 : * @param[in] dt The length of the time step of the Simulation.
462 : */
463 : void interaction(TimePoint t, TimeSpan dt);
464 : /**
465 : * @brief Person%s change location in the Model according to rules.
466 : * @param[in] t The current TimePoint.
467 : * @param[in] dt The length of the time step of the Simulation.
468 : */
469 : void perform_mobility(TimePoint t, TimeSpan dt);
470 :
471 : /// @brief Shape the cache and store how many Person%s are at any Location. Use from single thread!
472 : void build_compute_local_population_cache() const;
473 :
474 : /// @brief Shape the air and contact exposure cache according to the current Location%s.
475 : void build_exposure_caches();
476 :
477 : /**
478 : * @brief Store all air/contact exposures for the current simulation step.
479 : * @param[in] t Current TimePoint of the simulation.
480 : * @param[in] dt The duration of the simulation step.
481 : */
482 : void compute_exposure_caches(TimePoint t, TimeSpan dt);
483 :
484 : mutable Eigen::Matrix<std::atomic_int_fast32_t, Eigen::Dynamic, 1>
485 : m_local_population_cache; ///< Current number of Persons in a given location.
486 : Eigen::Matrix<AirExposureRates, Eigen::Dynamic, 1>
487 : m_air_exposure_rates_cache; ///< Cache for local exposure through droplets in #transmissions/day.
488 : Eigen::Matrix<ContactExposureRates, Eigen::Dynamic, 1>
489 : m_contact_exposure_rates_cache; ///< Cache for local exposure through contacts in #transmissions/day.
490 : bool m_is_local_population_cache_valid = false;
491 : bool m_are_exposure_caches_valid = false;
492 : bool m_exposure_caches_need_rebuild = true;
493 :
494 : std::vector<Person> m_persons; ///< Vector of every Person.
495 : std::vector<Location> m_locations; ///< Vector of every Location.
496 : std::bitset<size_t(LocationType::Count)>
497 : m_has_locations; ///< Flags for each LocationType, set if a Location of that type exists.
498 : TestingStrategy m_testing_strategy; ///< List of TestingScheme%s that are checked for testing.
499 : TripList m_trip_list; ///< List of all Trip%s the Person%s do.
500 : bool m_use_mobility_rules; ///< Whether mobility rules are considered.
501 : std::vector<std::pair<LocationType (*)(PersonalRandomNumberGenerator&, const Person&, TimePoint, TimeSpan,
502 : const Parameters&),
503 : std::vector<LocationType>>>
504 : m_mobility_rules; ///< Rules that govern the mobility between Location%s.
505 : LocationId m_cemetery_id; // Central cemetery for all dead persons.
506 : RandomNumberGenerator m_rng; ///< Global random number generator
507 : };
508 :
509 : } // namespace abm
510 : } // namespace mio
511 :
512 : #endif
|
