Changes according to review

Author: DavidKerkmann

cpp/models/abm/model.cpp

@@ -72,8 +72,8 @@ PersonId Model::add_person(Person&& person)
     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_per_age = m_local_population_per_age_cache[new_person.get_location().get()];
-            ++local_population_per_age[{cell, new_person.get_age()}];
+            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();
@@ -224,17 +224,17 @@ void Model::build_compute_local_population_cache() const
         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_per_age_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_per_age_cache[i].resize(
+            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_per_age          = m_local_population_per_age_cache[i];
-                    local_population_per_age[{cell, group}] = 0;
+                    auto& local_population_by_age          = m_local_population_by_age_cache[i];
+                    local_population_by_age[{cell, group}] = 0;
                 }
             }
         } // implicit taskloop barrier
@@ -244,9 +244,8 @@ void Model::build_compute_local_population_cache() const
                 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_per_age =
-                        m_local_population_per_age_cache[m_persons[i].get_location().get()];
-                    ++local_population_per_age[{cell, m_persons[i].get_age()}];
+                    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
@@ -315,11 +314,16 @@ void Model::compute_exposure_caches(TimePoint t, TimeSpan 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_per_age_cache[location],
-                                                      parameters.get_num_groups());
+                                                      m_local_population_by_age_cache[location]);
         } // implicit taskloop barrier
 
     } // implicit single barrier

cpp/models/abm/model.h

@@ -94,7 +94,7 @@ class Model
     Model(const Model& other, int id = 0)
         : parameters(other.parameters)
         , m_local_population_cache()
-        , m_local_population_per_age_cache()
+        , m_local_population_by_age_cache()
         , m_air_exposure_rates_cache()
         , m_contact_exposure_rates_cache()
         , m_is_local_population_cache_valid(false)
@@ -415,8 +415,8 @@ class Model
             build_compute_local_population_cache();
             m_is_local_population_cache_valid = true;
         }
-        auto& m_local_population_per_age = m_local_population_per_age_cache[location.get()];
-        return m_local_population_per_age[{cell_idx, age}];
+        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.
@@ -564,12 +564,12 @@ class Model
                 --m_local_population_cache[origin.get()];
                 ++m_local_population_cache[destination.get()];
                 for (CellIndex cell : old_cells) {
-                    auto& local_population_per_age = m_local_population_per_age_cache[origin.get()];
-                    --local_population_per_age[{cell, person.get_age()}];
+                    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_per_age = m_local_population_per_age_cache[destination.get()];
-                    ++local_population_per_age[{cell, person.get_age()}];
+                    auto& local_population_by_age = m_local_population_by_age_cache[destination.get()];
+                    ++local_population_by_age[{cell, person.get_age()}];
                 }
             }
         }
@@ -621,7 +621,7 @@ class Model
         auto personal_rng = PersonalRandomNumberGenerator(person);
         auto location     = person.get_location();
         mio::abm::interact(personal_rng, person, get_location(location),
-                           m_local_population_per_age_cache[location.get()], m_air_exposure_rates_cache[location.get()],
+                           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);
     }
 
@@ -653,8 +653,8 @@ 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<CustomIndexArray<std::atomic_int_fast32_t, CellIndex, AgeGroup>, Eigen::Dynamic, 1>
-        m_local_population_per_age_cache; ///<Current number of Persons per AgeGroup 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>

cpp/models/abm/model_functions.cpp

@@ -34,19 +34,17 @@ namespace abm
 
 ScalarType daily_transmissions_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
                                            const VirusVariant virus, const AgeGroup age_receiver,
-                                           size_t num_persons_in_age_receiverer_group,
-                                           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) {
         if (age_receiver == age_transmitter &&
-            num_persons_in_age_receiverer_group > 1) // adjust for the person not meeting themself
+            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}] * num_persons_in_age_receiverer_group /
-                    (num_persons_in_age_receiverer_group - 1) *
-                    params.get<ContactRates>()[{age_receiver, 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}] *
@@ -63,9 +61,9 @@ ScalarType daily_transmissions_by_air(const AirExposureRates& rates, const CellI
 }
 
 void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const Location& location,
-              const CustomIndexArray<std::atomic_int_fast32_t, CellIndex, AgeGroup>& local_population_per_age,
-              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());
@@ -87,11 +85,11 @@ void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const
         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);
-                size_t local_population_per_age_receiver = local_population_per_age[{cell_index, age_receiver}];
+                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_population_per_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));
 
@@ -142,19 +140,18 @@ void add_exposure_contribution(AirExposureRates& local_air_exposure, ContactExpo
     }
 }
 
-void normalize_exposure_contribution(
-    ContactExposureRates& local_contact_exposure,
-    const CustomIndexArray<std::atomic_int_fast32_t, CellIndex, AgeGroup>& local_population_per_age,
-    size_t num_agegroups)
+void normalize_exposure_contribution(ContactExposureRates& local_contact_exposure,
+                                     const PopulationByAge& local_population_by_age)
 {
-    for (CellIndex cell = 0; cell < local_contact_exposure.size<CellIndex>(); ++cell) {
-        for (auto&& virus : enum_members<VirusVariant>()) {
-            for (auto age_group = AgeGroup(0); age_group < AgeGroup(num_agegroups); age_group++) {
-                if (local_population_per_age[{cell, age_group}] > 0) {
-                    local_contact_exposure[{cell, virus, age_group}] =
-                        local_contact_exposure[{cell, virus, age_group}] / local_population_per_age[{cell, age_group}];
-                }
-            }
+    // 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) {
+            local_contact_exposure[index] = local_contact_exposure[index] / local_population_by_age[age_index];
         }
     }
 }

cpp/models/abm/model_functions.h

@@ -37,14 +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] num_persons_in_age_receiverer_group Number of persons in the 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,
-                                           size_t num_persons_in_age_receiverer_group,
-                                           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.
@@ -71,33 +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 by number of people in age groups
- * @param[in, out] local_contact_exposure Exposure by rates contacts for the local population.
- * @param[in] local_population_per_age Local population per AgeGroup.
- * @param[in] num_agegroup The number of AgeGroups in the model.
+ * @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 CustomIndexArray<std::atomic_int_fast32_t, CellIndex, AgeGroup>& local_population_per_age,
-    size_t num_agegroups);
+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_per_age Local population per AgeGroup at the 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 CustomIndexArray<std::atomic_int_fast32_t, CellIndex, AgeGroup>& local_population_per_age,
-              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.

cpp/tests/abm_helpers.cpp

@@ -65,19 +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 per age group and cell
-    mio::CustomIndexArray<std::atomic_int_fast32_t, mio::abm::CellIndex, mio::AgeGroup> local_population_per_age;
-    local_population_per_age.resize(
+    // 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_per_age.begin(), local_population_per_age.end(), [](auto& r) {
+    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_per_age[{cell_index, p.get_age()}]++;
+            local_population_by_age[{cell_index, p.get_age()}]++;
         }
     }
     // run interaction
-    mio::abm::interact(personal_rng, person, location, local_population_per_age, local_air_exposure,
+    mio::abm::interact(personal_rng, person, location, local_population_by_age, local_air_exposure,
                        local_contact_exposure, t, dt, global_parameters);
 }

cpp/tests/test_abm_model.cpp

@@ -216,12 +216,12 @@ TEST_F(TestModel, exposureContributionNormalization)
     contact_exposure_rates.resize({mio::abm::CellIndex{1}, mio::abm::VirusVariant{1}, mio::AgeGroup{1}});
     contact_exposure_rates[{mio::abm::CellIndex{0}, mio::abm::VirusVariant::Wildtype, mio::AgeGroup{0}}] = 10.0;
 
-    mio::CustomIndexArray<std::atomic_int_fast32_t, mio::abm::CellIndex, mio::AgeGroup> local_population_per_age;
-    local_population_per_age.resize({mio::abm::CellIndex{1}, mio::AgeGroup{1}});
-    local_population_per_age[{mio::abm::CellIndex{0}, mio::AgeGroup{0}}] =
+    mio::abm::PopulationByAge local_population_by_age;
+    local_population_by_age.resize({mio::abm::CellIndex{1}, mio::AgeGroup{1}});
+    local_population_by_age[{mio::abm::CellIndex{0}, mio::AgeGroup{0}}] =
         2; // Set population for cell 0 and age group 0 to 2
 
-    mio::abm::normalize_exposure_contribution(contact_exposure_rates, local_population_per_age, 1);
+    mio::abm::normalize_exposure_contribution(contact_exposure_rates, local_population_by_age);
 
     auto& rate = contact_exposure_rates[{mio::abm::CellIndex{0}, mio::abm::VirusVariant::Wildtype, mio::AgeGroup{0}}];
     EXPECT_EQ(rate, 5.0);