abm/model.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Daniel Abele, Majid Abedi, Elisabeth Kluth, David Kerkmann, Sascha Korf, Martin J. Kuehn, Khoa Nguyen
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef MIO_ABM_MODEL_H
#define MIO_ABM_MODEL_H
 
#include "abm/infection_state.h"
#include "abm/model_functions.h"
#include "abm/location_type.h"
#include "abm/mobility_data.h"
#include "abm/parameters.h"
#include "abm/location.h"
#include "abm/person.h"
#include "abm/person_id.h"
#include "abm/time.h"
#include "abm/trip_list.h"
#include "abm/random_events.h"
#include "abm/testing_strategy.h"
#include "memilio/epidemiology/age_group.h"
#include "memilio/utils/logging.h"
#include "memilio/utils/random_number_generator.h"
#include "memilio/utils/stl_util.h"
 
#include <bitset>
#include <cstdint>
#include <vector>
 
namespace mio
{
namespace abm
{
 
class Model
{
public:
    using LocationIterator        = std::vector<Location>::iterator;
    using ConstLocationIterator   = std::vector<Location>::const_iterator;
    using PersonIterator          = std::vector<Person>::iterator;
    using ConstPersonIterator     = std::vector<Person>::const_iterator;
    using ActivenessIterator      = std::vector<char>::iterator;
    using ConstActivenessIterator = std::vector<char>::const_iterator;
    using MobilityRuleType        = LocationType (*)(PersonalRandomNumberGenerator&, const Person&, TimePoint, TimeSpan,
                                              const Parameters&);
 
    using Compartments = mio::abm::InfectionState;
    Model(size_t num_agegroups, int id = 0)
        : parameters(num_agegroups)
        , m_id(id)
        , m_trip_list()
        , m_use_mobility_rules(true)
        , m_cemetery_id(add_location(LocationType::Cemetery))
        , m_person_ids_equal_index(true)
    {
        assert(num_agegroups < MAX_NUM_AGE_GROUPS && "MAX_NUM_AGE_GROUPS exceeded.");
    }
 
    Model(const Parameters& params, int id = 0)
        : parameters(params.get_num_groups())
        , m_id(id)
        , m_trip_list()
        , m_use_mobility_rules(true)
        , m_cemetery_id(add_location(LocationType::Cemetery))
        , m_person_ids_equal_index(true)
    {
        parameters = params;
    }
 
    Model(const Model& other, int id = 0)
        : parameters(other.parameters)
        , m_local_population_cache()
        , m_local_population_by_age_cache()
        , m_air_exposure_rates_cache()
        , m_contact_exposure_rates_cache()
        , m_is_local_population_cache_valid(false)
        , m_are_exposure_caches_valid(false)
        , m_exposure_caches_need_rebuild(true)
        , m_id(id)
        , m_persons(other.m_persons)
        , m_locations(other.m_locations)
        , m_activeness_statuses(other.m_activeness_statuses)
        , m_has_locations(other.m_has_locations)
        , m_testing_strategy(other.m_testing_strategy)
        , m_trip_list(other.m_trip_list)
        , m_use_mobility_rules(other.m_use_mobility_rules)
        , m_mobility_rules(other.m_mobility_rules)
        , m_cemetery_id(other.m_cemetery_id)
        , m_rng(other.m_rng)
        , m_person_ids_equal_index(true)
    {
    }
    Model& operator=(const Model&) = default;
    Model(Model&&)                 = default;
    Model& operator=(Model&&)      = default;
 
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("Model");
        obj.add_element("parameters", parameters);
        // skip caches, they are rebuild by the deserialized model
        obj.add_list("persons", get_persons().begin(), get_persons().end());
        obj.add_list("locations", get_locations().begin(), get_locations().end());
        obj.add_list("activeness_statuses", m_activeness_statuses.begin(), m_activeness_statuses.end());
        obj.add_element("location_types", m_has_locations.to_ulong());
        obj.add_element("testing_strategy", m_testing_strategy);
        obj.add_element("trip_list", m_trip_list);
        obj.add_element("use_mobility_rules", m_use_mobility_rules);
        obj.add_element("cemetery_id", m_cemetery_id);
        obj.add_element("rng", m_rng);
    }
 
    template <class IOContext>
    static IOResult<Model> deserialize(IOContext& io)
    {
        auto obj                 = io.expect_object("Model");
        auto params              = obj.expect_element("parameters", Tag<Parameters>{});
        auto persons             = obj.expect_list("persons", Tag<Person>{});
        auto locations           = obj.expect_list("locations", Tag<Location>{});
        auto activeness_statuses = obj.expect_list("activeness_statuses", Tag<bool>{});
        auto location_types      = obj.expect_element("location_types", Tag<unsigned long>{});
        auto trip_list           = obj.expect_element("trip_list", Tag<TripList>{});
        auto use_mobility_rules  = obj.expect_element("use_mobility_rules", Tag<bool>{});
        auto cemetery_id         = obj.expect_element("cemetery_id", Tag<LocationId>{});
        auto rng                 = obj.expect_element("rng", Tag<RandomNumberGenerator>{});
        return apply(
            io,
            [](auto&& params_, auto&& persons_, auto&& locations_, auto&& activeness_statuses_, auto&& location_types_,
               auto&& trip_list_, auto&& use_mobility_rules_, auto&& cemetery_id_, auto&& rng_) {
                Model model{params_};
                model.m_persons.assign(persons_.cbegin(), persons_.cend());
                model.m_locations.assign(locations_.cbegin(), locations_.cend());
                model.m_activeness_statuses.assign(activeness_statuses_.cbegin(), activeness_statuses_.cend());
                model.m_has_locations      = location_types_;
                model.m_trip_list          = trip_list_;
                model.m_use_mobility_rules = use_mobility_rules_;
                model.m_cemetery_id        = cemetery_id_;
                model.m_rng                = rng_;
                return model;
            },
            params, persons, locations, activeness_statuses, location_types, trip_list, use_mobility_rules, cemetery_id,
            rng);
    }
 
    void begin_step(TimePoint t, TimeSpan dt);
 
    void evolve(TimePoint t, TimeSpan dt);
 
    LocationId add_location(LocationType type, uint32_t num_cells = 1);
 
    PersonId add_person(const LocationId id, AgeGroup age);
 
    PersonId add_person(Person&& person);
 
    Range<ConstLocationIterator> get_locations() const;
    Range<LocationIterator> get_locations();
    Range<ConstPersonIterator> get_persons() const;
    Range<PersonIterator> get_persons();
    Range<ConstActivenessIterator> get_activeness_statuses() const;
    Range<ActivenessIterator> get_activeness_statuses();
    LocationId find_location(LocationType type, const PersonId person) const;
 
    void assign_location(PersonId person, LocationId location)
    {
        assign_location(get_person(person), location);
    }
 
    size_t get_subpopulation_combined(TimePoint t, InfectionState s) const;
 
    size_t get_subpopulation_combined_per_location_type(TimePoint t, InfectionState s, LocationType type) const;
 
    TripList& get_trip_list();
 
    const TripList& get_trip_list() const;
 
    void use_mobility_rules(bool param);
    bool use_mobility_rules() const;
 
    bool has_location(LocationType type) const
    {
        return m_has_locations[size_t(type)];
    }
 
    template <class C = std::initializer_list<LocationType>>
    bool has_locations(const C& location_types) const
    {
        return std::all_of(location_types.begin(), location_types.end(), [&](auto loc) {
            return has_location(loc);
        });
    }
 
    TestingStrategy& get_testing_strategy();
 
    const TestingStrategy& get_testing_strategy() const;
 
    Parameters parameters;
 
    RandomNumberGenerator& get_rng()
    {
        return m_rng;
    }
    const RandomNumberGenerator& get_rng() const
    {
        return m_rng;
    }
    void reset_rng(Counter<uint32_t> counter = {})
    {
        m_rng.set_counter(Counter<uint64_t>(static_cast<uint64_t>(counter.get())));
        for (Person& person : get_persons()) {
            person.get_rng_counter() = counter;
        }
    }
    void reset_rng(const std::vector<uint32_t>& seeds, Counter<uint32_t> counter = {})
    {
        m_rng.seed(seeds);
        reset_rng(counter);
    }
    int get_id() const
    {
        return m_id;
    }
 
    const Person& get_person(PersonId person_id) const
    {
        return get_person_impl(*this, person_id);
    }
 
    Person& get_person(PersonId person_id)
    {
        return get_person_impl(*this, person_id);
    }
 
    size_t get_subpopulation(LocationId location, TimePoint t, InfectionState state) const
    {
        return std::count_if(m_persons.begin(), m_persons.end(), [&](auto&& p) {
            return p.get_location_model_id() == m_id && p.get_location() == location &&
                   p.get_infection_state(t) == state;
        });
    }
 
    size_t get_number_persons(LocationId location) const
    {
        if (!m_is_local_population_cache_valid) {
            build_compute_local_population_cache();
            m_is_local_population_cache_valid = true;
        }
        return m_local_population_cache[location.get()];
    }
 
    size_t get_number_persons_age(LocationId location, CellIndex cell_idx, AgeGroup age) const
    {
        if (!m_is_local_population_cache_valid) {
            build_compute_local_population_cache();
            m_is_local_population_cache_valid = true;
        }
        auto& m_local_population_by_age = m_local_population_by_age_cache[location.get()];
        return m_local_population_by_age[{cell_idx, age}];
    }
 
    // Change the Location of a Person. this requires that Location is part of this Model.
    inline void change_location(PersonId person, LocationId destination, TransportMode mode = TransportMode::Unknown,
                                const std::vector<uint32_t>& cells = {0})
    {
        change_location(get_person(person), destination, mode, cells);
    }
 
    inline void interact(PersonId person, TimePoint t, TimeSpan dt)
    {
        interact(get_person(person), t, dt);
    }
 
    const Location& get_location(LocationId id) const
    {
        assert(id != LocationId::invalid_id() && "Given LocationId must be valid.");
        assert(id < LocationId((uint32_t)m_locations.size()) && "Given LocationId is not in this Model.");
        return m_locations[id.get()];
    }
 
    void assign_location(Person& person, LocationId location)
    {
        person.set_assigned_location(get_location(location).get_type(), location, m_id);
    }
 
    Location& get_location(LocationId id)
    {
        assert(id != LocationId::invalid_id() && "Given LocationId must be valid.");
        assert(id < LocationId((uint32_t)m_locations.size()) && "Given LocationId is not in this Model.");
        return m_locations[id.get()];
    }
    inline Location& get_location(PersonId id)
    {
        return get_location(get_person(id).get_location());
    }
 
    inline const Location& get_location(PersonId id) const
    {
        return get_location(get_person(id).get_location());
    }
    uint32_t get_person_index(PersonId person_id) const
    {
        mio::log_debug("get_person_index is used leading to a search in m_persons.");
        auto it = std::find_if(m_persons.begin(), m_persons.end(), [person_id](auto& person) {
            return person.get_id() == person_id;
        });
        if (it == m_persons.end()) {
            log_error("Given PersonId is not in this Model.");
            return std::numeric_limits<uint32_t>::max();
        }
        else {
            return static_cast<uint32_t>(std::distance(m_persons.begin(), it));
        }
    }
 
protected:
    void interaction(TimePoint t, TimeSpan dt);
    void perform_mobility(TimePoint t, TimeSpan dt);
 
    void build_compute_local_population_cache() const;
 
    void build_exposure_caches();
 
    void compute_exposure_caches(TimePoint t, TimeSpan dt);
 
    // Change the Location of a Person. this requires that Location is part of this Model.
    inline void change_location(Person& person, LocationId destination, TransportMode mode = TransportMode::Unknown,
                                const std::vector<uint32_t>& cells = {0})
    {
        LocationId origin               = get_location(person).get_id();
        auto old_cells                  = person.get_cells();
        const bool has_changed_location = mio::abm::change_location(person, get_location(destination), mode, cells);
        // if the person has changed location, invalidate exposure caches but keep population caches valid
        if (has_changed_location) {
            m_are_exposure_caches_valid = false;
            if (m_is_local_population_cache_valid) {
                --m_local_population_cache[origin.get()];
                ++m_local_population_cache[destination.get()];
                for (CellIndex cell : old_cells) {
                    auto& local_population_by_age = m_local_population_by_age_cache[origin.get()];
                    --local_population_by_age[{cell, person.get_age()}];
                }
                for (CellIndex cell : cells) {
                    auto& local_population_by_age = m_local_population_by_age_cache[destination.get()];
                    ++local_population_by_age[{cell, person.get_age()}];
                }
            }
        }
    }
 
    inline Location& get_location(Person& person)
    {
        return get_location(person.get_location());
    }
 
    inline const Location& get_location(Person& person) const
    {
        return get_location(person.get_location());
    }
 
    LocationId find_location(LocationType type, const Person& person) const
    {
        auto location_id = person.get_assigned_location(type);
        assert(location_id != LocationId::invalid_id() && "The person has no assigned location of that type.");
        return location_id;
    }
 
    inline void interact(Person& person, TimePoint t, TimeSpan dt)
    {
        if (!m_are_exposure_caches_valid) {
            // checking caches is only needed for external calls
            // during simulation (i.e. in evolve()), the caches are computed in begin_step
            compute_exposure_caches(t, dt);
            m_are_exposure_caches_valid = true;
        }
        auto personal_rng = PersonalRandomNumberGenerator(m_rng, person);
        auto location     = person.get_location();
        mio::abm::interact(personal_rng, person, get_location(location),
                           m_local_population_by_age_cache[location.get()], m_air_exposure_rates_cache[location.get()],
                           m_contact_exposure_rates_cache[location.get()], t, dt, parameters);
    }
 
    template <class M>
    static std::conditional_t<std::is_const_v<M>, const Person&, Person&> get_person_impl(M& m, PersonId person_id)
    {
        if (m.m_person_ids_equal_index) {
            assert(static_cast<uint32_t>(person_id.get()) < m.m_persons.size() &&
                   "Given PersonId is not in this Model.");
            return m.m_persons[static_cast<uint32_t>(person_id.get())];
        }
        else {
            mio::log_warning("get_person is accessed by PersonId which does not align with the index of the person due "
                             "to former removal of persons. Therefore m_persons is searched.");
            auto it = std::find_if(m.m_persons.begin(), m.m_persons.end(), [person_id](auto& person) {
                return person.get_id() == person_id;
            });
            assert(it != m.m_persons.end() && "Given PersonId is not in this Model.");
            return *it;
        };
    }
 
    mutable Eigen::Matrix<std::atomic_int_fast32_t, Eigen::Dynamic, 1>
        m_local_population_cache; 
    mutable Eigen::Matrix<PopulationByAge, Eigen::Dynamic, 1>
        m_local_population_by_age_cache; 
    Eigen::Matrix<AirExposureRates, Eigen::Dynamic, 1>
        m_air_exposure_rates_cache; 
    Eigen::Matrix<ContactExposureRates, Eigen::Dynamic, 1>
        m_contact_exposure_rates_cache; 
    mutable bool m_is_local_population_cache_valid = false;
    bool m_are_exposure_caches_valid               = false;
    bool m_exposure_caches_need_rebuild            = true;
 
    int m_id; 
    std::vector<Person> m_persons; 
    std::vector<Location> m_locations; 
    std::vector<char> m_activeness_statuses; 
    std::bitset<size_t(LocationType::Count)>
        m_has_locations; 
    TestingStrategy m_testing_strategy; 
    TripList m_trip_list; 
    bool m_use_mobility_rules; 
    std::vector<MobilityRuleType> m_mobility_rules; 
    LocationId m_cemetery_id; // Central cemetery for all dead persons.
    RandomNumberGenerator m_rng; 
    bool m_person_ids_equal_index;
};
 
} // namespace abm
} // namespace mio
 
#endif