1274 normalizing of exposure rate (#1322)

DavidKerkmann · reneSchm · web-flow · commit cff9ccb8709f · 2025-10-28T07:56:48.000+01:00
- Add a cache for the number of persons per cell and age group.
- Use this cache to normalize the exposure rate according to the number of persons inside a cell.

Co-authored-by: reneSchm &lt;49305466+reneSchm@users.noreply.github.com&gt;
diff --git a/cpp/models/abm/model.cpp b/cpp/models/abm/model.cpp
@@ -71,6 +71,10 @@ PersonId Model::add_person(Person&& person)
     auto& new_person = m_persons.back();
     if (m_is_local_population_cache_valid) {
         ++m_local_population_cache[new_person.get_location().get()];
+        for (CellIndex cell : new_person.get_cells()) {
+            auto& local_population_by_age = m_local_population_by_age_cache[new_person.get_location().get()];
+            ++local_population_by_age[{cell, new_person.get_age()}];
+        }
     }
     return new_person.get_id();
 }
@@ -221,18 +225,33 @@ void Model::build_compute_local_population_cache() const
 {
     PRAGMA_OMP(single)
     {
+        log_debug("Building local population cache for ABM.");
         const size_t num_locations = m_locations.size();
         const size_t num_persons   = m_persons.size();
         m_local_population_cache.resize(num_locations);
+        m_local_population_by_age_cache.resize(num_locations);
+
         PRAGMA_OMP(taskloop)
         for (size_t i = 0; i < num_locations; i++) {
             m_local_population_cache[i] = 0;
+            m_local_population_by_age_cache[i].resize(
+                {CellIndex(m_locations[i].get_cells().size()), AgeGroup(parameters.get_num_groups())});
+            for (CellIndex cell = 0; cell < CellIndex(m_locations[i].get_cells().size()); ++cell) {
+                for (AgeGroup group = AgeGroup(0); group < AgeGroup(parameters.get_num_groups()); ++group) {
+                    auto& local_population_by_age          = m_local_population_by_age_cache[i];
+                    local_population_by_age[{cell, group}] = 0;
+                }
+            }
         } // implicit taskloop barrier
         PRAGMA_OMP(taskloop)
         for (size_t i = 0; i < num_persons; i++) {
             if (m_activeness_statuses[i]) {
                 assert(m_persons[i].get_location_model_id() == m_id && "Person is not in this model but still active.");
                 ++m_local_population_cache[m_persons[i].get_location().get()];
+                for (CellIndex cell : m_persons[i].get_cells()) {
+                    auto& local_population_by_age = m_local_population_by_age_cache[m_persons[i].get_location().get()];
+                    ++local_population_by_age[{cell, m_persons[i].get_age()}];
+                }
             }
         } // implicit taskloop barrier
     } // implicit single barrier
@@ -242,6 +261,7 @@ void Model::build_exposure_caches()
 {
     PRAGMA_OMP(single)
     {
+        mio::log_debug("Building exposure caches for ABM.");
         const size_t num_locations = m_locations.size();
         m_air_exposure_rates_cache.resize(num_locations);
         m_contact_exposure_rates_cache.resize(num_locations);
@@ -297,6 +317,20 @@ void Model::compute_exposure_caches(TimePoint t, TimeSpan dt)
                                                     get_location(person.get_location()), parameters, t, dt);
             }
         } // implicit taskloop barrier
+
+        // 3) normalize contributions for each location
+
+        // make sure that the caches are computed
+        if (!m_is_local_population_cache_valid) {
+            build_compute_local_population_cache();
+            m_is_local_population_cache_valid = true;
+        }
+        PRAGMA_OMP(taskloop)
+        for (size_t location = 0; location < num_locations; ++location) {
+            mio::abm::normalize_exposure_contribution(m_contact_exposure_rates_cache[location],
+                                                      m_local_population_by_age_cache[location]);
+        } // implicit taskloop barrier
+
     } // implicit single barrier
 }
 
diff --git a/cpp/models/abm/model.h b/cpp/models/abm/model.h
@@ -94,6 +94,7 @@ class Model
     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)
@@ -374,16 +375,34 @@ class Model
     /**
      * @brief Get the total number of Person%s at the Location.
      * @param[in] location A LocationId from the Model.
-     * @return Number of Person%s in the location.
+     * @return Number of Person%s at the location.
      */
     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()];
     }
 
+    /**
+     * @brief Get the number of Person%s of a specific AgeGroup in a specific Cell at the Location.
+     * @param[in] location A LocationId from the Model.
+     * @param[in] cell_idx Index of the Cell.
+     * @param[in] AgeGroup An AgeGroup from the Model.
+     * @return Number of Person%s of the AgeGroup in the Cell at the Location.
+     */
+    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.
     /**
      * @brief Let a Person change to another Location.
@@ -503,6 +522,7 @@ class Model
 
     /**
      * @brief Store all air/contact exposures for the current simulation step.
+     * This will also compute the local population cache if it is not valid, as it is required for the computation of the exposure rates.
      * @param[in] t Current TimePoint of the simulation.
      * @param[in] dt The duration of the simulation step.
      */
@@ -520,13 +540,22 @@ class Model
                                 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()}];
+                }
             }
         }
     }
@@ -575,9 +604,10 @@ class Model
             m_are_exposure_caches_valid = true;
         }
         auto personal_rng = PersonalRandomNumberGenerator(person);
-        mio::abm::interact(personal_rng, person, get_location(person.get_location()),
-                           m_air_exposure_rates_cache[person.get_location().get()],
-                           m_contact_exposure_rates_cache[person.get_location().get()], t, dt, parameters);
+        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);
     }
 
     /**
@@ -608,13 +638,15 @@ class Model
 
     mutable Eigen::Matrix<std::atomic_int_fast32_t, Eigen::Dynamic, 1>
         m_local_population_cache; ///< Current number of Persons in a given location.
+    mutable Eigen::Matrix<PopulationByAge, Eigen::Dynamic, 1>
+        m_local_population_by_age_cache; ///<Current number of Persons per AgeGroup in a given location.
     Eigen::Matrix<AirExposureRates, Eigen::Dynamic, 1>
         m_air_exposure_rates_cache; ///< Cache for local exposure through droplets in #transmissions/day.
     Eigen::Matrix<ContactExposureRates, Eigen::Dynamic, 1>
         m_contact_exposure_rates_cache; ///< Cache for local exposure through contacts in #transmissions/day.
-    bool m_is_local_population_cache_valid = false;
-    bool m_are_exposure_caches_valid       = false;
-    bool m_exposure_caches_need_rebuild    = true;
+    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; ///< Model id. Is only used for abm graph model or hybrid model.
     std::vector<Person> m_persons; ///< Vector of every Person.
diff --git a/cpp/models/abm/model_functions.cpp b/cpp/models/abm/model_functions.cpp
@@ -34,13 +34,22 @@ namespace abm
 
 ScalarType daily_transmissions_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
                                            const VirusVariant virus, const AgeGroup age_receiver,
-                                           const LocalInfectionParameters& params)
+                                           size_t age_receiver_group_size, const LocalInfectionParameters& params)
 {
     assert(age_receiver < rates.size<AgeGroup>());
     ScalarType prob = 0;
     for (AgeGroup age_transmitter(0); age_transmitter < rates.size<AgeGroup>(); ++age_transmitter) {
-        prob +=
-            rates[{cell_index, virus, age_transmitter}] * params.get<ContactRates>()[{age_receiver, age_transmitter}];
+        if (age_receiver == age_transmitter &&
+            age_receiver_group_size > 1) // adjust for the person not meeting themself
+        {
+            prob += rates[{cell_index, virus, age_transmitter}] *
+                    params.get<ContactRates>()[{age_receiver, age_transmitter}] * age_receiver_group_size /
+                    (age_receiver_group_size - 1);
+        }
+        else {
+            prob += rates[{cell_index, virus, age_transmitter}] *
+                    params.get<ContactRates>()[{age_receiver, age_transmitter}];
+        }
     }
     return prob;
 }
@@ -52,8 +61,9 @@ ScalarType daily_transmissions_by_air(const AirExposureRates& rates, const CellI
 }
 
 void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const Location& location,
-              const AirExposureRates& local_air_exposure, const ContactExposureRates& local_contact_exposure,
-              const TimePoint t, const TimeSpan dt, const Parameters& global_parameters)
+              const PopulationByAge& local_population_by_age, const AirExposureRates& local_air_exposure,
+              const ContactExposureRates& local_contact_exposure, const TimePoint t, const TimeSpan dt,
+              const Parameters& global_parameters)
 {
     // make sure all dimensions are set correctly and all indices are valid
     assert(location.get_cells().size() == local_air_exposure.size<CellIndex>().get());
@@ -72,13 +82,14 @@ void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const
         auto age_receiver          = person.get_age();
         ScalarType mask_protection = person.get_mask_protective_factor(global_parameters);
         assert(person.get_cells().size() && "Person is in multiple cells. Interact logic is incorrect at the moment.");
-        for (auto cell_index : person.get_cells()) {
+        for (CellIndex cell_index : person.get_cells()) {
             std::pair<VirusVariant, ScalarType> local_indiv_trans_prob[static_cast<uint32_t>(VirusVariant::Count)];
             for (uint32_t v = 0; v != static_cast<uint32_t>(VirusVariant::Count); ++v) {
-                VirusVariant virus = static_cast<VirusVariant>(v);
+                VirusVariant virus                      = static_cast<VirusVariant>(v);
+                size_t local_population_by_age_receiver = local_population_by_age[{cell_index, age_receiver}];
                 ScalarType local_indiv_trans_prob_v =
                     (daily_transmissions_by_contacts(local_contact_exposure, cell_index, virus, age_receiver,
-                                                     local_parameters) +
+                                                     local_population_by_age_receiver, local_parameters) +
                      daily_transmissions_by_air(local_air_exposure, cell_index, virus, global_parameters)) *
                     (1 - mask_protection) * (1 - person.get_protection_factor(t, virus, global_parameters));
 
@@ -129,6 +140,23 @@ void add_exposure_contribution(AirExposureRates& local_air_exposure, ContactExpo
     }
 }
 
+void normalize_exposure_contribution(ContactExposureRates& local_contact_exposure,
+                                     const PopulationByAge& local_population_by_age)
+{
+    // make sure all dimensions are set correctly and all indices are valid
+    assert(local_population_by_age.size<AgeGroup>() == local_contact_exposure.size<AgeGroup>());
+    assert(local_population_by_age.size<CellIndex>() == local_contact_exposure.size<CellIndex>());
+    assert(local_contact_exposure.size<VirusVariant>() == VirusVariant::Count);
+
+    for (auto index : make_index_range(local_contact_exposure.size())) {
+        auto age_index = reduce_index<Index<CellIndex, AgeGroup>>(index);
+        if (local_population_by_age[age_index] > 0) {
+            // this instruction is not and does not need to be atomic
+            local_contact_exposure[index] = local_contact_exposure[index] / local_population_by_age[age_index];
+        }
+    }
+}
+
 bool change_location(Person& person, const Location& destination, const TransportMode mode,
                      const std::vector<uint32_t>& cells)
 {
diff --git a/cpp/models/abm/model_functions.h b/cpp/models/abm/model_functions.h
@@ -37,12 +37,13 @@ namespace abm
  * @param[in] cell_index Cell index of the Cell.
  * @param[in] virus VirusVariant of interest.
  * @param[in] age_receiver AgeGroup of the receiving Person.
+ * @param[in] age_receiver_group_size Number of persons in the AgeGroup of the receiving Person.
  * @param[in] params The local infection parameters.
  * @return Average amount of Infection%s with the virus from the AgeGroup of the transmitter per day.
  */
 ScalarType daily_transmissions_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
                                            const VirusVariant virus, const AgeGroup age_receiver,
-                                           const LocalInfectionParameters& params);
+                                           size_t age_receiver_group_size, const LocalInfectionParameters& params);
 
 /**
  * @brief Compute the number of daily transmissions for aerosol transmission of a virus in a cell.
@@ -69,20 +70,32 @@ void add_exposure_contribution(AirExposureRates& local_air_exposure, ContactExpo
                                const Person& person, const Location& location, const Parameters& params,
                                const TimePoint t, const TimeSpan dt);
 
+using PopulationByAge = CustomIndexArray<std::atomic_int_fast32_t, CellIndex, AgeGroup>;
+
+/**
+ * @brief Normalize contact exposure rate to average exposure per contact per time (from total exposure per time).
+ * @param[in, out] local_contact_exposure Exposure rates through contacts for the local population.
+ * @param[in] local_population_by_age Local population by AgeGroup%s.
+ */
+void normalize_exposure_contribution(ContactExposureRates& local_contact_exposure,
+                                     const PopulationByAge& local_population_by_age);
+
 /**
  * @brief Let a Person interact with the population at its current Location, possibly getting infected.
  * @param[in, out] rng PersonalRandomNumberGenerator for this Person.
  * @param[in, out] person The person to interact with the local population.
  * @param[in] location The person's current location.
+ * @param[in] local_population_by_age Local population by AgeGroup%s at the given location.
  * @param[in] local_air_exposure Precomputed exposure rates by aerosols for the local population.
  * @param[in] local_contact_exposure Precomputed exposure by rates contacts for the local population.
  * @param[in] t Current Simulation time.
  * @param[in] dt Length of the current Simulation time step.
  * @param[in] global_parameters Parameters of the Model.
  */
 void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const Location& location,
-              const AirExposureRates& local_air_exposure, const ContactExposureRates& local_contact_exposure,
-              const TimePoint t, const TimeSpan dt, const Parameters& global_parameters);
+              const PopulationByAge& local_population_by_age, const AirExposureRates& local_air_exposure,
+              const ContactExposureRates& local_contact_exposure, const TimePoint t, const TimeSpan dt,
+              const Parameters& global_parameters);
 /**
  * @brief Change a persons location to another location.
  * If the person already is at the destination, neither mode nor cells are set.
diff --git a/cpp/tests/abm_helpers.cpp b/cpp/tests/abm_helpers.cpp
@@ -65,7 +65,19 @@ void interact_testing(mio::abm::PersonalRandomNumberGenerator& personal_rng, mio
     for (const mio::abm::Person& p : local_population) {
         add_exposure_contribution(local_air_exposure, local_contact_exposure, p, location, global_parameters, t, dt);
     }
+    // count local population by age group and cell
+    mio::abm::PopulationByAge local_population_by_age;
+    local_population_by_age.resize(
+        {mio::abm::CellIndex(location.get_cells().size()), mio::AgeGroup(global_parameters.get_num_groups())});
+    std::for_each(local_population_by_age.begin(), local_population_by_age.end(), [](auto& r) {
+        r = 0; // initialize with 0
+    });
+    for (const mio::abm::Person& p : local_population) {
+        for (mio::abm::CellIndex cell_index : p.get_cells()) {
+            local_population_by_age[{cell_index, p.get_age()}]++;
+        }
+    }
     // run interaction
-    mio::abm::interact(personal_rng, person, location, local_air_exposure, local_contact_exposure, t, dt,
-                       global_parameters);
+    mio::abm::interact(personal_rng, person, location, local_population_by_age, local_air_exposure,
+                       local_contact_exposure, t, dt, global_parameters);
 }
diff --git a/cpp/tests/test_abm_model.cpp b/cpp/tests/test_abm_model.cpp
diff --git a/cpp/tests/utils.h b/cpp/tests/utils.h
