Implement stats calculations for new calculated_stats_ container

src/HealthGPS/analysis_module.cpp

@@ -318,36 +318,95 @@ DALYsIndicator AnalysisModule::calculate_dalys(Population &population, unsigned
  .disability_adjusted_life_years = yll + yld};
 }
 
+void AnalysisModule::update_death_and_migration_stats(const Person &person, size_t index,
+ RuntimeContext &context) {
+
+ auto current_time = static_cast<unsigned int>(context.time_now());
+
+ if (!person.is_alive() && person.time_of_death() == current_time) {
+ calculated_stats_[index + channel_index_.at("deaths")]++;
+ float expected_life = definition_.life_expectancy().at(context.time_now(), person.gender);
+ double yll = std::max(expected_life - person.age, 0.0f) * DALY_UNITS;
+ calculated_stats_[index + channel_index_.at("mean_yll")] += yll;
+ calculated_stats_[index + channel_index_.at("mean_daly")] += yll;
+ }
+
+ if (person.has_emigrated() && person.time_of_migration() == current_time) {
+ calculated_stats_[index + channel_index_.at("emigrations")]++;
+ }
+}
+
+void AnalysisModule::update_calculated_stats_for_person(RuntimeContext &context,
+ const Person &person, size_t index) {
+ calculated_stats_[index + channel_index_.at("count")]++;
+
+ for (const auto &factor : context.mapping().entries()) {
+ double value = person.get_risk_factor_value(factor.key());
+ calculated_stats_[index + channel_index_.at("mean_" + factor.key().to_string())] += value;
+ }
+
+ for (const auto &[disease_name, disease_state] : person.diseases) {
+ if (disease_state.status == DiseaseStatus::active) {
+ calculated_stats_[index +
+ channel_index_.at("prevalence_" + disease_name.to_string())]++;
+ if (disease_state.start_time == context.time_now()) {
+ calculated_stats_[index +
+ channel_index_.at("incidence_" + disease_name.to_string())]++;
+ }
+ }
+ }
+}
+
 void AnalysisModule::calculate_population_statistics(RuntimeContext &context) {
- size_t num_factors_to_calculate =
- context.mapping().entries().size() - factors_to_calculate_.size();
 
  for (const auto &person : context.population()) {
- // Get the bin index for each factor
- std::vector<size_t> bin_indices;
- for (size_t i = 0; i < factors_to_calculate_.size(); i++) {
- double factor_value = person.get_risk_factor_value(factors_to_calculate_[i]);
- auto bin_index =
- static_cast<size_t>((factor_value - factor_min_values_[i]) / factor_bin_widths_[i]);
- bin_indices.push_back(bin_index);
- }
+ // First let's fetch the correct `calculated_stats_` bin index for this person
+ size_t index = calculate_index(person);
 
- // Calculate the index in the calculated_stats_ vector
- size_t index = 0;
- for (size_t i = 0; i < bin_indices.size() - 1; i++) {
- size_t accumulated_bins =
- std::accumulate(std::next(factor_bins_.cbegin(), i + 1), factor_bins_.cend(),
- size_t{1}, std::multiplies<>());
- index += bin_indices[i] * accumulated_bins * num_factors_to_calculate;
+ if (!person.is_active()) {
+ update_death_and_migration_stats(person, index, context);
+ continue;
  }
- index += bin_indices.back() * num_factors_to_calculate;
 
- // Now we can add the values of the factors that are not in factors_to_calculate_
+ update_calculated_stats_for_person(context, person, index);
+
+ double dw = calculate_disability_weight(person);
+ double yld = dw * DALY_UNITS;
+ calculated_stats_[index + channel_index_.at("mean_yld")] += yld;
+ calculated_stats_[index + channel_index_.at("mean_daly")] += yld;
+
+ classify_weight(person);
+ }
+
+ // For each bin in the calculated stats...
+ for (size_t i = 0; i < calculated_stats_.size(); i += channels_.size()) {
+ double count_F = calculated_stats_[i + channel_index_.at("count")];
+ double count_M = calculated_stats_[i + channel_index_.at("count")];
+ double deaths_F = calculated_stats_[i + channel_index_.at("deaths")];
+ double deaths_M = calculated_stats_[i + channel_index_.at("deaths")];
+
+ // Calculate in-place factor averages.
  for (const auto &factor : context.mapping().entries()) {
- if (std::find(factors_to_calculate_.cbegin(), factors_to_calculate_.cend(),
- factor.key()) == factors_to_calculate_.cend()) {
- calculated_stats_[index++] += person.get_risk_factor_value(factor.key());
- }
+ calculated_stats_[i + channel_index_.at("mean_" + factor.key().to_string())] /= count_F;
+ calculated_stats_[i + channel_index_.at("mean_" + factor.key().to_string())] /= count_M;
+ }
+
+ // Calculate in-place disease prevalence and incidence rates.
+ for (const auto &disease : context.diseases()) {
+ calculated_stats_[i + channel_index_.at("prevalence_" + disease.code.to_string())] /=
+ count_F;
+ calculated_stats_[i + channel_index_.at("prevalence_" + disease.code.to_string())] /=
+ count_M;
+ calculated_stats_[i + channel_index_.at("incidence_" + disease.code.to_string())] /=
+ count_F;
+ calculated_stats_[i + channel_index_.at("incidence_" + disease.code.to_string())] /=
+ count_M;
+ }
+
+ // Calculate in-place YLL/YLD/DALY averages.
+ for (const auto &column : {"mean_yll", "mean_yld", "mean_daly"}) {
+ calculated_stats_[i + channel_index_.at(column)] /= (count_F + deaths_F);
+ calculated_stats_[i + channel_index_.at(column)] /= (count_M + deaths_M);
  }
  }
 }
@@ -531,6 +590,25 @@ void AnalysisModule::classify_weight(DataSeries &series, const Person &entity) c
  }
 }
 
+void AnalysisModule::classify_weight(const Person &person) {
+ auto weight_class = weight_classifier_.classify_weight(person);
+ switch (weight_class) {
+ case WeightCategory::normal:
+ calculated_stats_[channel_index_.at("normal_weight")]++;
+ break;
+ case WeightCategory::overweight:
+ calculated_stats_[channel_index_.at("over_weight")]++;
+ calculated_stats_[channel_index_.at("above_weight")]++;
+ break;
+ case WeightCategory::obese:
+ calculated_stats_[channel_index_.at("obese_weight")]++;
+ calculated_stats_[channel_index_.at("above_weight")]++;
+ break;
+ default:
+ throw std::logic_error("Unknown weight classification category.");
+ }
+}
+
 void AnalysisModule::initialise_output_channels(RuntimeContext &context) {
  if (!channels_.empty()) {
  return;
@@ -560,6 +638,35 @@ void AnalysisModule::initialise_output_channels(RuntimeContext &context) {
  channels_.emplace_back("std_yld");
  channels_.emplace_back("mean_daly");
  channels_.emplace_back("std_daly");
+
+ // Since we will be performing frequent lookups, we will store the strings and indexes in a map
+ // for quick access.
+ for (size_t i = 0; i < channels_.size(); i++) {
+ channel_index_.emplace(channels_[i], i);
+ }
+}
+
+size_t AnalysisModule::calculate_index(const Person &person) const {
+ // Get the bin index for each factor
+ std::vector<size_t> bin_indices;
+ for (size_t i = 0; i < factors_to_calculate_.size(); i++) {
+ double factor_value = person.get_risk_factor_value(factors_to_calculate_[i]);
+ auto bin_index =
+ static_cast<size_t>((factor_value - factor_min_values_[i]) / factor_bin_widths_[i]);
+ bin_indices.push_back(bin_index);
+ }
+
+ // Calculate the index in the calculated_stats_ vector
+ size_t index = 0;
+ for (size_t i = 0; i < bin_indices.size() - 1; i++) {
+ size_t accumulated_bins =
+ std::accumulate(std::next(factor_bins_.cbegin(), i + 1), factor_bins_.cend(), size_t{1},
+ std::multiplies<>());
+ index += bin_indices[i] * accumulated_bins * channels_.size();
+ }
+ index += bin_indices.back() * channels_.size();
+
+ return index;
 }
 
 std::unique_ptr<AnalysisModule> build_analysis_module(Repository &repository,

src/HealthGPS/analysis_module.h

@@ -46,6 +46,7 @@ class AnalysisModule final : public UpdatableModule {
  WeightModel weight_classifier_;
  DoubleAgeGenderTable residual_disability_weight_;
  std::vector<std::string> channels_;
+ std::unordered_map<std::string, size_t> channel_index_;
  unsigned int comorbidities_;
  std::string name_{"Analysis"};
  std::vector<core::Identifier> factors_to_calculate_ = {"Gender"_id, "Age"_id};
@@ -65,13 +66,23 @@ class AnalysisModule final : public UpdatableModule {
  double calculate_disability_weight(const Person &entity) const;
  DALYsIndicator calculate_dalys(Population &population, unsigned int max_age,
  unsigned int death_year) const;
+ void update_death_and_migration_stats(const Person &person, size_t index,
+ RuntimeContext &context);
+ void update_calculated_stats_for_person(RuntimeContext &context, const Person &person,
+ size_t index);
 
  void calculate_population_statistics(RuntimeContext &context);
  void calculate_population_statistics(RuntimeContext &context, DataSeries &series) const;
 
  void classify_weight(hgps::DataSeries &series, const hgps::Person &entity) const;
+ void classify_weight(const Person &person);
  void initialise_output_channels(RuntimeContext &context);
 
+ /// @brief Calculates the bin index in `calculated_stats_` for a given person
+ /// @param person The person to calculate the index for
+ /// @return The index in `calculated_stats_`
+ size_t calculate_index(const Person &person) const;
+
  /// @brief Calculates the standard deviation of factors given data series containing means
  /// @param context The runtime context
  /// @param series The data series containing factor means