Sector and income

example_new/KevinHall.json

@@ -67,6 +67,43 @@
             "Description": "string"
         }
     },
+    "RuralPrevalence": [
+        {
+            "AgeRange": [0, 17],
+            "Female": 0.65,
+            "Male": 0.65
+        },
+        {
+            "AgeRange": [18, 100],
+            "Female": 0.6,
+            "Male": 0.6
+        }
+    ],
+    "IncomeModels": [
+        {
+            "Name": "Category_1",
+            "Intercept": 1.0,
+            "Coefficients": {}
+        },
+        {
+            "Name": "Category_2",
+            "Intercept": 0.75498278113169,
+            "Coefficients": {
+                "Sex": 0.0204338261883629,
+                "Over18": 0.0486699373325097,
+                "Sector": -0.706477682886734
+            }
+        },
+        {
+            "Name": "Category_3",
+            "Intercept": 1.48856946873517,
+            "Coefficients": {
+                "Sex": 0.000641255498596025,
+                "Over18": 0.206622749570132,
+                "Sector": -1.71313287940798
+            }
+        }
+    ],
     "AgeMeanHeight": {
         "Male": [
             0.498, 0.757, 0.868, 0.952, 1.023, 1.092, 1.155, 1.219, 1.28, 1.333,

src/HealthGPS.Console/model_parser.cpp

@@ -129,31 +129,50 @@ std::unique_ptr<hgps::StaticLinearModelDefinition>
 load_staticlinear_risk_model_definition(const poco::json &opt, const host::Configuration &config) {
     MEASURE_FUNCTION();
 
-    // Risk factor linear models.
-    hgps::LinearModelParams models;
-    for (const auto &factor : opt["RiskFactorModels"]) {
-        auto factor_key = factor["Name"].get<hgps::core::Identifier>();
-        models.intercepts[factor_key] = factor["Intercept"].get<double>();
-        models.coefficients[factor_key] =
-            factor["Coefficients"].get<std::unordered_map<hgps::core::Identifier, double>>();
-    }
-
-    // Risk factor names and correlation matrix.
-    std::vector<hgps::core::Identifier> names;
+    // Risk factor linear models and correlation matrix.
+    std::vector<hgps::LinearModelParams> risk_factor_models;
     const auto correlations_file_info =
         host::get_file_info(opt["RiskFactorCorrelationFile"], config.root_path);
     const auto correlations_table = load_datatable_from_csv(correlations_file_info);
     Eigen::MatrixXd correlations{correlations_table.num_rows(), correlations_table.num_columns()};
-    for (size_t col = 0; col < correlations_table.num_columns(); col++) {
-        names.emplace_back(correlations_table.column(col).name());
-        for (size_t row = 0; row < correlations_table.num_rows(); row++) {
-            correlations(row, col) =
-                std::any_cast<double>(correlations_table.column(col).value(row));
+
+    for (size_t i = 0; i < opt["RiskFactorModels"].size(); i++) {
+        // Risk factor model.
+        const auto &factor = opt["RiskFactorModels"][i];
+        hgps::LinearModelParams model;
+        model.name = factor["Name"].get<hgps::core::Identifier>();
+        model.intercept = factor["Intercept"].get<double>();
+        model.coefficients =
+            factor["Coefficients"].get<std::unordered_map<hgps::core::Identifier, double>>();
+        risk_factor_models.emplace_back(std::move(model));
+
+        // Correlation matrix column.
+        for (size_t j = 0; j < correlations_table.num_rows(); j++) {
+            correlations(i, j) = std::any_cast<double>(correlations_table.column(i).value(j));
         }
+
+        // Check correlation matrix column name matches risk factor name.
+        auto column_name = hgps::core::Identifier{correlations_table.column(i).name()};
+        if (model.name != column_name) {
+            throw hgps::core::HgpsException{
+                fmt::format("Risk factor {} name ({}) does not match correlation matrix "
+                            "column {} name ({})",
+                            i, model.name.to_string(), i, column_name.to_string())};
+        }
+    }
+
+    // Check correlation matrix column count matches risk factor count.
+    if (opt["RiskFactorModels"].size() != correlations_table.num_columns()) {
+        throw hgps::core::HgpsException{
+            fmt::format("Risk factor count ({}) does not match correlation "
+                        "matrix column count ({})",
+                        opt["RiskFactorModels"].size(), correlations_table.num_columns())};
     }
+
+    // Compute Cholesky decomposition of correlation matrix.
     auto cholesky = Eigen::MatrixXd{Eigen::LLT<Eigen::MatrixXd>{correlations}.matrixL()};
 
-    return std::make_unique<hgps::StaticLinearModelDefinition>(std::move(names), std::move(models),
+    return std::make_unique<hgps::StaticLinearModelDefinition>(std::move(risk_factor_models),
                                                                std::move(cholesky));
 }
 
@@ -253,14 +272,10 @@ load_ebhlm_risk_model_definition(const poco::json &opt) {
 std::unique_ptr<hgps::KevinHallModelDefinition>
 load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configuration &config) {
     MEASURE_FUNCTION();
-    std::unordered_map<hgps::core::Identifier, double> energy_equation;
-    std::unordered_map<hgps::core::Identifier, std::pair<double, double>> nutrient_ranges;
-    std::unordered_map<hgps::core::Identifier, std::map<hgps::core::Identifier, double>>
-        nutrient_equations;
-    std::unordered_map<hgps::core::Identifier, std::optional<double>> food_prices;
-    std::unordered_map<hgps::core::Gender, std::vector<double>> age_mean_height;
 
     // Nutrient groups.
+    std::unordered_map<hgps::core::Identifier, double> energy_equation;
+    std::unordered_map<hgps::core::Identifier, std::pair<double, double>> nutrient_ranges;
     for (const auto &nutrient : opt["Nutrients"]) {
         auto nutrient_key = nutrient["Name"].get<hgps::core::Identifier>();
         nutrient_ranges[nutrient_key] = nutrient["Range"].get<std::pair<double, double>>();
@@ -272,6 +287,9 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur
     }
 
     // Food groups.
+    std::unordered_map<hgps::core::Identifier, std::map<hgps::core::Identifier, double>>
+        nutrient_equations;
+    std::unordered_map<hgps::core::Identifier, std::optional<double>> food_prices;
     for (const auto &food : opt["Foods"]) {
         auto food_key = food["Name"].get<hgps::core::Identifier>();
         food_prices[food_key] = food["Price"].get<std::optional<double>>();
@@ -291,7 +309,28 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur
     const auto food_data_file_info = host::get_file_info(opt["FoodsDataFile"], config.root_path);
     const auto food_data_table = load_datatable_from_csv(food_data_file_info);
 
+    // Rural sector prevalence for age groups and sex.
+    std::map<hgps::core::IntegerInterval, std::unordered_map<hgps::core::Gender, double>>
+        rural_prevalence;
+    for (const auto &age_group : opt["RuralPrevalence"]) {
+        auto age_group_range = age_group["AgeRange"].get<hgps::core::IntegerInterval>();
+        rural_prevalence[age_group_range] = {{hgps::core::Gender::female, age_group["Female"]},
+                                             {hgps::core::Gender::male, age_group["Male"]}};
+    }
+
+    // Income models for different income classifications.
+    std::vector<hgps::LinearModelParams> income_models;
+    for (const auto &factor : opt["IncomeModels"]) {
+        hgps::LinearModelParams model;
+        model.name = factor["Name"].get<hgps::core::Identifier>();
+        model.intercept = factor["Intercept"].get<double>();
+        model.coefficients =
+            factor["Coefficients"].get<std::unordered_map<hgps::core::Identifier, double>>();
+        income_models.emplace_back(std::move(model));
+    }
+
     // Load M/F average heights for age.
+    std::unordered_map<hgps::core::Gender, std::vector<double>> age_mean_height;
     const auto max_age = static_cast<size_t>(config.settings.age_range.upper());
     auto male_height = opt["AgeMeanHeight"]["Male"].get<std::vector<double>>();
     auto female_height = opt["AgeMeanHeight"]["Female"].get<std::vector<double>>();
@@ -306,7 +345,8 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur
 
     return std::make_unique<hgps::KevinHallModelDefinition>(
         std::move(energy_equation), std::move(nutrient_ranges), std::move(nutrient_equations),
-        std::move(food_prices), std::move(age_mean_height));
+        std::move(food_prices), std::move(rural_prevalence), std::move(income_models),
+        std::move(age_mean_height));
 }
 
 std::pair<hgps::RiskFactorModelType, std::unique_ptr<hgps::RiskFactorModelDefinition>>

src/HealthGPS.Core/forward_type.h

@@ -35,6 +35,18 @@ enum class DiseaseGroup : uint8_t {
     cancer
 };
 
+/// @brief Enumerates sector types
+enum class Sector : uint8_t {
+    /// @brief Unknown sector
+    unknown,
+
+    /// @brief Urban sector
+    urban,
+
+    /// @brief Rural sector
+    rural
+};
+
 /// @brief C++20 concept for numeric columns types
 template <typename T>
 concept Numerical = std::is_arithmetic_v<T>;
@@ -60,4 +72,5 @@ class DoubleDataTableColumn;
 class IntegerDataTableColumn;
 
 class DataTableColumnVisitor;
+
 } // namespace hgps::core

src/HealthGPS/dynamic_hierarchical_linear_model.cpp

@@ -26,10 +26,7 @@ RiskFactorModelType DynamicHierarchicalLinearModel::type() const noexcept {
 std::string DynamicHierarchicalLinearModel::name() const noexcept { return "Dynamic"; }
 
 void DynamicHierarchicalLinearModel::generate_risk_factors(
-    [[maybe_unused]] RuntimeContext &context) {
-    throw core::HgpsException(
-        "DynamicHierarchicalLinearModel::generate_risk_factors not yet implemented.");
-}
+    [[maybe_unused]] RuntimeContext &context) {}
 
 void DynamicHierarchicalLinearModel::update_risk_factors(RuntimeContext &context) {
     auto age_key = core::Identifier{"age"};

src/HealthGPS/dynamic_hierarchical_linear_model.h

@@ -48,7 +48,6 @@ class DynamicHierarchicalLinearModel final : public RiskFactorModel {
 
     std::string name() const noexcept override;
 
-    /// @throws std::logic_error the dynamic model does not generate risk factors.
     void generate_risk_factors(RuntimeContext &context) override;
 
     void update_risk_factors(RuntimeContext &context) override;

src/HealthGPS/kevin_hall_model.cpp

@@ -33,9 +33,13 @@ KevinHallModel::KevinHallModel(
     const std::unordered_map<core::Identifier, std::map<core::Identifier, double>>
         &nutrient_equations,
     const std::unordered_map<core::Identifier, std::optional<double>> &food_prices,
+    const std::map<hgps::core::IntegerInterval, std::unordered_map<hgps::core::Gender, double>>
+        &rural_prevalence,
+    const std::vector<LinearModelParams> &income_models,
     const std::unordered_map<core::Gender, std::vector<double>> &age_mean_height)
     : energy_equation_{energy_equation}, nutrient_ranges_{nutrient_ranges},
       nutrient_equations_{nutrient_equations}, food_prices_{food_prices},
+      rural_prevalence_{rural_prevalence}, income_models_{income_models},
       age_mean_height_{age_mean_height} {
 
     if (energy_equation_.empty()) {
@@ -50,6 +54,12 @@ KevinHallModel::KevinHallModel(
     if (food_prices_.empty()) {
         throw core::HgpsException("Food price mapping is empty");
     }
+    if (rural_prevalence_.empty()) {
+        throw core::HgpsException("Rural prevalence mapping is empty");
+    }
+    if (income_models_.empty()) {
+        throw core::HgpsException("Income models list is empty");
+    }
     if (age_mean_height_.empty()) {
         throw core::HgpsException("Age mean height mapping is empty");
     }
@@ -59,19 +69,34 @@ RiskFactorModelType KevinHallModel::type() const noexcept { return RiskFactorMod
 
 std::string KevinHallModel::name() const noexcept { return "Dynamic"; }
 
-void KevinHallModel::generate_risk_factors([[maybe_unused]] RuntimeContext &context) {
-    throw core::HgpsException("KevinHallModel::generate_risk_factors not yet implemented.");
+void KevinHallModel::generate_risk_factors(RuntimeContext &context) {
+
+    // Initialise sector for everyone.
+    for (auto &person : context.population()) {
+        initialise_sector(context, person);
+    }
 }
 
 void KevinHallModel::update_risk_factors(RuntimeContext &context) {
     hgps::Population &population = context.population();
-    double mean_sim_body_weight = 0.0;
-    double mean_adjustment_coefficient = 0.0;
+
+    // Initialise sector for newborns.
+    for (auto &person : population) {
+
+        // Do not initialise non-newborns.
+        if (person.age > 0) {
+            continue;
+        }
+
+        initialise_sector(context, person);
+    }
 
     // TODO: Compute target body weight.
     const float target_BW = 100.0;
 
     // Trial run.
+    double mean_sim_body_weight = 0.0;
+    double mean_adjustment_coefficient = 0.0;
     for (auto &person : population) {
         // Ignore if inactive.
         if (!person.is_active()) {
@@ -114,6 +139,61 @@ void KevinHallModel::update_risk_factors(RuntimeContext &context) {
     }
 }
 
+void KevinHallModel::initialise_sector(RuntimeContext &context, Person &person) const {
+
+    // Get prevalence for age and sex (default is zero).
+    double rural_prevalence = 0.0;
+    for (const auto &[age_range, prevalence_by_sex] : rural_prevalence_) {
+        if (age_range.contains(person.age)) {
+            rural_prevalence = prevalence_by_sex.at(person.gender);
+            break;
+        }
+    }
+
+    // Sample the person's sector.
+    double rand = context.random().next_double();
+    auto sector = rand < rural_prevalence ? core::Sector::rural : core::Sector::urban;
+    person.sector = sector;
+}
+
+void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) const {
+
+    // Compute logits for each income category.
+    auto logits = std::vector<double>(income_models_.size());
+    for (size_t i = 0; i < income_models_.size(); i++) {
+        logits[i] = income_models_[i].intercept;
+        for (const auto &[factor_name, coefficient] : income_models_[i].coefficients) {
+            logits[i] += coefficient * person.get_risk_factor_value(factor_name);
+        }
+    }
+
+    // Compute softmax probabilities for each income category.
+    auto e_logits = std::vector<double>(income_models_.size());
+    double e_logits_sum = 0.0;
+    for (size_t i = 0; i < income_models_.size(); i++) {
+        e_logits[i] = exp(logits[i]);
+        e_logits_sum += e_logits[i];
+    }
+
+    // Compute income category probabilities.
+    auto probabilities = std::vector<double>(income_models_.size());
+    for (size_t i = 0; i < income_models_.size(); i++) {
+        probabilities[i] = e_logits[i] / e_logits_sum;
+    }
+
+    // Compute income category.
+    double rand = context.random().next_double();
+    for (size_t i = 0; i < probabilities.size(); i++) {
+        if (rand < probabilities[i]) {
+            person.income = income_models_[i].name;
+            return;
+        }
+        rand -= probabilities[i];
+    }
+
+    throw core::HgpsException("Logic Error: failed to initialise income category");
+}
+
 SimulatePersonState KevinHallModel::simulate_person(Person &person, double shift) const {
     // Initial simulated person state.
     const double H_0 = person.get_risk_factor_value(H_key);
@@ -280,9 +360,13 @@ KevinHallModelDefinition::KevinHallModelDefinition(
     std::unordered_map<core::Identifier, std::pair<double, double>> nutrient_ranges,
     std::unordered_map<core::Identifier, std::map<core::Identifier, double>> nutrient_equations,
     std::unordered_map<core::Identifier, std::optional<double>> food_prices,
+    std::map<hgps::core::IntegerInterval, std::unordered_map<hgps::core::Gender, double>>
+        rural_prevalence,
+    std::vector<LinearModelParams> income_models,
     std::unordered_map<core::Gender, std::vector<double>> age_mean_height)
     : energy_equation_{std::move(energy_equation)}, nutrient_ranges_{std::move(nutrient_ranges)},
       nutrient_equations_{std::move(nutrient_equations)}, food_prices_{std::move(food_prices)},
+      rural_prevalence_{std::move(rural_prevalence)}, income_models_{std::move(income_models)},
       age_mean_height_{std::move(age_mean_height)} {
 
     if (energy_equation_.empty()) {
@@ -297,14 +381,21 @@ KevinHallModelDefinition::KevinHallModelDefinition(
     if (food_prices_.empty()) {
         throw core::HgpsException("Food prices mapping is empty");
     }
+    if (rural_prevalence_.empty()) {
+        throw core::HgpsException("Rural prevalence mapping is empty");
+    }
+    if (income_models_.empty()) {
+        throw core::HgpsException("Income models list is empty");
+    }
     if (age_mean_height_.empty()) {
         throw core::HgpsException("Age mean height mapping is empty");
     }
 }
 
 std::unique_ptr<RiskFactorModel> KevinHallModelDefinition::create_model() const {
     return std::make_unique<KevinHallModel>(energy_equation_, nutrient_ranges_, nutrient_equations_,
-                                            food_prices_, age_mean_height_);
+                                            food_prices_, rural_prevalence_, income_models_,
+                                            age_mean_height_);
 }
 
 } // namespace hgps