From de623fd1959e5b78d920bf08eceee3ecb924e117 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 12:13:50 +0100 Subject: [PATCH 01/29] Add SectorPrevalence to KevinHall.json config. --- example_new/KevinHall.json | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/example_new/KevinHall.json b/example_new/KevinHall.json index 0ae394a1f..69ec3416c 100644 --- a/example_new/KevinHall.json +++ b/example_new/KevinHall.json @@ -67,6 +67,18 @@ "Description": "string" } }, + "SectorPrevalence": { + "Rural": { + "Male": { + "18": 0.65, + "100": 0.6 + }, + "Female": { + "18": 0.65, + "100": 0.6 + } + } + }, "AgeMeanHeight": { "Male": [ 0.498, 0.757, 0.868, 0.952, 1.023, 1.092, 1.155, 1.219, 1.28, 1.333, From d638cef24444fc59231498cc6d8a21a889734fd0 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 14:50:33 +0100 Subject: [PATCH 02/29] Load Rural Prevalence into model parser from JSON. --- example_new/KevinHall.json | 22 +++++++++++----------- src/HealthGPS.Console/model_parser.cpp | 10 ++++++++++ 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/example_new/KevinHall.json b/example_new/KevinHall.json index 69ec3416c..ba860f21c 100644 --- a/example_new/KevinHall.json +++ b/example_new/KevinHall.json @@ -67,18 +67,18 @@ "Description": "string" } }, - "SectorPrevalence": { - "Rural": { - "Male": { - "18": 0.65, - "100": 0.6 - }, - "Female": { - "18": 0.65, - "100": 0.6 - } + "RuralPrevalence": [ + { + "AgeRange": [0, 17], + "Female": 0.65, + "Male": 0.65 + }, + { + "AgeRange": [18, 100], + "Female": 0.6, + "Male": 0.6 } - }, + ], "AgeMeanHeight": { "Male": [ 0.498, 0.757, 0.868, 0.952, 1.023, 1.092, 1.155, 1.219, 1.28, 1.333, diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index ca65fad00..4a3523ce5 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -291,6 +291,16 @@ 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> + rural_prevalence; + for (const auto &age_group : opt["RuralPrevalence"]) { + auto age_group_range = age_group["AgeRange"].get(); + rural_prevalence[std::move(age_group_range)] = { + {hgps::core::Gender::female, age_group["Female"]}, + {hgps::core::Gender::male, age_group["Male"]}}; + } + // Load M/F average heights for age. const auto max_age = static_cast(config.settings.age_range.upper()); auto male_height = opt["AgeMeanHeight"]["Male"].get>(); From 585e93630266fb5567a81fb53986e28b857fd936 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 15:39:08 +0100 Subject: [PATCH 03/29] Load Rural Prevalence into KH model. --- src/HealthGPS.Console/model_parser.cpp | 14 +++++++------- src/HealthGPS/kevin_hall_model.cpp | 13 ++++++++++--- src/HealthGPS/kevin_hall_model.h | 10 ++++++++++ 3 files changed, 27 insertions(+), 10 deletions(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index 4a3523ce5..1d2a03245 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -253,14 +253,10 @@ load_ebhlm_risk_model_definition(const poco::json &opt) { std::unique_ptr load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configuration &config) { MEASURE_FUNCTION(); - std::unordered_map energy_equation; - std::unordered_map> nutrient_ranges; - std::unordered_map> - nutrient_equations; - std::unordered_map> food_prices; - std::unordered_map> age_mean_height; // Nutrient groups. + std::unordered_map energy_equation; + std::unordered_map> nutrient_ranges; for (const auto &nutrient : opt["Nutrients"]) { auto nutrient_key = nutrient["Name"].get(); nutrient_ranges[nutrient_key] = nutrient["Range"].get>(); @@ -272,6 +268,9 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur } // Food groups. + std::unordered_map> + nutrient_equations; + std::unordered_map> food_prices; for (const auto &food : opt["Foods"]) { auto food_key = food["Name"].get(); food_prices[food_key] = food["Price"].get>(); @@ -302,6 +301,7 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur } // Load M/F average heights for age. + std::unordered_map> age_mean_height; const auto max_age = static_cast(config.settings.age_range.upper()); auto male_height = opt["AgeMeanHeight"]["Male"].get>(); auto female_height = opt["AgeMeanHeight"]["Female"].get>(); @@ -316,7 +316,7 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur return std::make_unique( 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(age_mean_height)); } std::pair> diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index b36bea1b9..7eec4ae36 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -33,10 +33,12 @@ KevinHallModel::KevinHallModel( const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, + const std::map> + &rural_prevalence, const std::unordered_map> &age_mean_height) : energy_equation_{energy_equation}, nutrient_ranges_{nutrient_ranges}, nutrient_equations_{nutrient_equations}, food_prices_{food_prices}, - age_mean_height_{age_mean_height} { + rural_prevalence_{rural_prevalence}, age_mean_height_{age_mean_height} { if (energy_equation_.empty()) { throw core::HgpsException("Energy equation mapping is empty"); @@ -280,10 +282,12 @@ KevinHallModelDefinition::KevinHallModelDefinition( std::unordered_map> nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, + std::map> + rural_prevalence, std::unordered_map> 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)}, - age_mean_height_{std::move(age_mean_height)} { + rural_prevalence_{std::move(rural_prevalence)}, age_mean_height_{std::move(age_mean_height)} { if (energy_equation_.empty()) { throw core::HgpsException("Energy equation mapping is empty"); @@ -297,6 +301,9 @@ 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 (age_mean_height_.empty()) { throw core::HgpsException("Age mean height mapping is empty"); } @@ -304,7 +311,7 @@ KevinHallModelDefinition::KevinHallModelDefinition( std::unique_ptr KevinHallModelDefinition::create_model() const { return std::make_unique(energy_equation_, nutrient_ranges_, nutrient_equations_, - food_prices_, age_mean_height_); + food_prices_, rural_prevalence_, age_mean_height_); } } // namespace hgps diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index 816f12e53..8316a9170 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -53,6 +53,7 @@ class KevinHallModel final : public RiskFactorModel { /// @param nutrient_ranges The minimum and maximum nutrient values /// @param nutrient_equations The nutrient coefficients for each food group /// @param food_prices The unit price for each food group + /// @param rural_prevalence Rural sector prevalence for age groups and sex /// @param age_mean_height The mean height at all ages (male and female) KevinHallModel( const std::unordered_map &energy_equation, @@ -60,6 +61,8 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, + const std::map> + &rural_prevalence, const std::unordered_map> &age_mean_height); RiskFactorModelType type() const noexcept override; @@ -77,6 +80,8 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &nutrient_equations_; const std::unordered_map> &food_prices_; + const std::map> + &rural_prevalence_; const std::unordered_map> &age_mean_height_; // Model parameters. @@ -164,6 +169,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { /// @param nutrient_ranges The minimum and maximum nutrient values /// @param nutrient_equations The nutrient coefficients for each food group /// @param food_prices The unit price for each food group + /// @param rural_prevalence Rural sector prevalence for age groups and sex /// @param age_mean_height The mean height at all ages (male and female) /// @throws std::invalid_argument for empty arguments KevinHallModelDefinition( @@ -171,6 +177,8 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, + std::map> + rural_prevalence, std::unordered_map> age_mean_height); /// @brief Construct a new KevinHallModel from this definition @@ -182,6 +190,8 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> nutrient_ranges_; std::unordered_map> nutrient_equations_; std::unordered_map> food_prices_; + std::map> + rural_prevalence_; std::unordered_map> age_mean_height_; }; From 192290661ff181222faf44de4a00b3147f0c8a1f Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 16:16:33 +0100 Subject: [PATCH 04/29] We want all our RF models' generate_ methods to run. --- src/HealthGPS/dynamic_hierarchical_linear_model.cpp | 5 +---- src/HealthGPS/dynamic_hierarchical_linear_model.h | 1 - src/HealthGPS/kevin_hall_model.cpp | 4 +--- src/HealthGPS/kevin_hall_model.h | 1 - src/HealthGPS/riskfactor.cpp | 3 +++ 5 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/HealthGPS/dynamic_hierarchical_linear_model.cpp b/src/HealthGPS/dynamic_hierarchical_linear_model.cpp index f30d0f097..b49ac10d5 100644 --- a/src/HealthGPS/dynamic_hierarchical_linear_model.cpp +++ b/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"}; diff --git a/src/HealthGPS/dynamic_hierarchical_linear_model.h b/src/HealthGPS/dynamic_hierarchical_linear_model.h index c8d7faa94..7871504c0 100644 --- a/src/HealthGPS/dynamic_hierarchical_linear_model.h +++ b/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; diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 7eec4ae36..ed6362a28 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -61,9 +61,7 @@ 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([[maybe_unused]] RuntimeContext &context) {} void KevinHallModel::update_risk_factors(RuntimeContext &context) { hgps::Population &population = context.population(); diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index 8316a9170..a765e53a7 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -69,7 +69,6 @@ class KevinHallModel 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; diff --git a/src/HealthGPS/riskfactor.cpp b/src/HealthGPS/riskfactor.cpp index 84b020f3d..3871067af 100644 --- a/src/HealthGPS/riskfactor.cpp +++ b/src/HealthGPS/riskfactor.cpp @@ -45,6 +45,9 @@ RiskFactorModel &RiskFactorModule::at(const RiskFactorModelType &model_type) con void RiskFactorModule::initialise_population(RuntimeContext &context) { auto &static_model = models_.at(RiskFactorModelType::Static); static_model->generate_risk_factors(context); + + auto &dynamic_model = models_.at(RiskFactorModelType::Dynamic); + dynamic_model->generate_risk_factors(context); } void RiskFactorModule::update_population(RuntimeContext &context) { From 3a6534eb18d7e46dbae1aafef0ba88c0d8ff4071 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 17:27:39 +0100 Subject: [PATCH 05/29] Add stub initialise_sector method to KH model class. --- src/HealthGPS/kevin_hall_model.cpp | 25 ++++++++++++++++++++++--- src/HealthGPS/kevin_hall_model.h | 4 ++++ 2 files changed, 26 insertions(+), 3 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index ed6362a28..c51eaf6cf 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -61,17 +61,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) {} +void KevinHallModel::generate_risk_factors([[maybe_unused]] RuntimeContext &context) { + + // Initialise sector for everyone. + for (auto &person : context.population()) { + initialise_sector(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(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 +131,8 @@ void KevinHallModel::update_risk_factors(RuntimeContext &context) { } } +void KevinHallModel::initialise_sector([[maybe_unused]] Person &person) const {} + SimulatePersonState KevinHallModel::simulate_person(Person &person, double shift) const { // Initial simulated person state. const double H_0 = person.get_risk_factor_value(H_key); diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index a765e53a7..ee4ffa333 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -95,6 +95,10 @@ class KevinHallModel final : public RiskFactorModel { static constexpr double xi_Na = 3000.0; // Na from ECF changes (mg/L/day). static constexpr double xi_CI = 4000.0; // Na from carbohydrate changes (mg/day). + /// @brief Initialise the sector of a person + /// @param person The person to initialise sector for + void initialise_sector(Person &person) const; + /// @brief Simulates the energy balance model for a given person /// @param person The person to simulate /// @param shift Model adjustment term From e5a312c851f32cc0f43dd0569755d1a03b1b0555 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 18:07:20 +0100 Subject: [PATCH 06/29] Add enum type for sector (region). --- src/HealthGPS.Core/forward_type.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/HealthGPS.Core/forward_type.h b/src/HealthGPS.Core/forward_type.h index 4f7852919..d0f16dc5f 100644 --- a/src/HealthGPS.Core/forward_type.h +++ b/src/HealthGPS.Core/forward_type.h @@ -35,6 +35,15 @@ enum class DiseaseGroup : uint8_t { cancer }; +/// @brief Enumerates sector types +enum class Sector : uint8_t { + /// @brief Urban sector + urban, + + /// @brief Rural sector + rural +}; + /// @brief C++20 concept for numeric columns types template concept Numerical = std::is_arithmetic_v; @@ -60,4 +69,5 @@ class DoubleDataTableColumn; class IntegerDataTableColumn; class DataTableColumnVisitor; + } // namespace hgps::core From 87c4eac5abe73b8a11b7dc84b110242d5d366cc8 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 18:21:26 +0100 Subject: [PATCH 07/29] Add 'unknown' sector enum type, and inititalise people to unknown sector. --- src/HealthGPS.Core/forward_type.h | 3 +++ src/HealthGPS/person.h | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/HealthGPS.Core/forward_type.h b/src/HealthGPS.Core/forward_type.h index d0f16dc5f..8f175fb10 100644 --- a/src/HealthGPS.Core/forward_type.h +++ b/src/HealthGPS.Core/forward_type.h @@ -37,6 +37,9 @@ enum class DiseaseGroup : uint8_t { /// @brief Enumerates sector types enum class Sector : uint8_t { + /// @brief Unknown sector + unknown, + /// @brief Urban sector urban, diff --git a/src/HealthGPS/person.h b/src/HealthGPS/person.h index 6c176cc8d..8cc1fd5f9 100644 --- a/src/HealthGPS/person.h +++ b/src/HealthGPS/person.h @@ -56,6 +56,9 @@ struct Person { /// @brief Current age in years unsigned int age{}; + /// @brief Sector (region) assigned value + core::Sector sector{core::Sector::unknown}; + /// @brief Social-economic status (SES) assigned value double ses{}; From 5f9a7fee670d3493a9583a2ac4451b5c1aede3d8 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 18:23:31 +0100 Subject: [PATCH 08/29] Finish initialise_sector code in KevinHall. --- src/HealthGPS/kevin_hall_model.cpp | 23 +++++++++++++++++++---- src/HealthGPS/kevin_hall_model.h | 3 ++- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index c51eaf6cf..22cebfd7a 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -61,11 +61,11 @@ RiskFactorModelType KevinHallModel::type() const noexcept { return RiskFactorMod std::string KevinHallModel::name() const noexcept { return "Dynamic"; } -void KevinHallModel::generate_risk_factors([[maybe_unused]] RuntimeContext &context) { +void KevinHallModel::generate_risk_factors(RuntimeContext &context) { // Initialise sector for everyone. for (auto &person : context.population()) { - initialise_sector(person); + initialise_sector(context, person); } } @@ -80,7 +80,7 @@ void KevinHallModel::update_risk_factors(RuntimeContext &context) { continue; } - initialise_sector(person); + initialise_sector(context, person); } // TODO: Compute target body weight. @@ -131,7 +131,22 @@ void KevinHallModel::update_risk_factors(RuntimeContext &context) { } } -void KevinHallModel::initialise_sector([[maybe_unused]] Person &person) const {} +void KevinHallModel::initialise_sector(RuntimeContext &context, Person &person) const { + + // Get prevalence for age and sex (default is zero). + double prevalence = 0.0; + for (const auto &[age_range, prevalence_by_sex] : rural_prevalence_) { + if (age_range.contains(person.age)) { + prevalence = prevalence_by_sex.at(person.gender); + break; + } + } + + // Sample the person's sector. + double rand = context.random().next_double(); + auto sector = rand < prevalence ? core::Sector::rural : core::Sector::urban; + person.sector = sector; +} SimulatePersonState KevinHallModel::simulate_person(Person &person, double shift) const { // Initial simulated person state. diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index ee4ffa333..e8dfd99ec 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -96,8 +96,9 @@ class KevinHallModel final : public RiskFactorModel { static constexpr double xi_CI = 4000.0; // Na from carbohydrate changes (mg/day). /// @brief Initialise the sector of a person + /// @param context The runtime context /// @param person The person to initialise sector for - void initialise_sector(Person &person) const; + void initialise_sector(RuntimeContext &context, Person &person) const; /// @brief Simulates the energy balance model for a given person /// @param person The person to simulate From eb3557267d0c3cd4e7b60913c725465c5e05097e Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 18:26:26 +0100 Subject: [PATCH 09/29] Clarify that we are using rural prevalence here. --- src/HealthGPS/kevin_hall_model.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 22cebfd7a..f211f78fa 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -134,17 +134,17 @@ 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 prevalence = 0.0; + double rural_prevalence = 0.0; for (const auto &[age_range, prevalence_by_sex] : rural_prevalence_) { if (age_range.contains(person.age)) { - prevalence = prevalence_by_sex.at(person.gender); + rural_prevalence = prevalence_by_sex.at(person.gender); break; } } // Sample the person's sector. double rand = context.random().next_double(); - auto sector = rand < prevalence ? core::Sector::rural : core::Sector::urban; + auto sector = rand < rural_prevalence ? core::Sector::rural : core::Sector::urban; person.sector = sector; } From c404c0a3413cb367cc249a05995d14929d1d216b Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Fri, 20 Oct 2023 19:10:53 +0100 Subject: [PATCH 10/29] Don't std move trivial type Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/HealthGPS.Console/model_parser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index 1d2a03245..3ba634ff1 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -295,7 +295,7 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur rural_prevalence; for (const auto &age_group : opt["RuralPrevalence"]) { auto age_group_range = age_group["AgeRange"].get(); - rural_prevalence[std::move(age_group_range)] = { + rural_prevalence[age_group_range] = { {hgps::core::Gender::female, age_group["Female"]}, {hgps::core::Gender::male, age_group["Male"]}}; } From f4d6638a331f43e93dc6f43526a3b2b4c2faa2ea Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Mon, 23 Oct 2023 14:56:47 +0100 Subject: [PATCH 11/29] model parser clang-tidy. --- src/HealthGPS.Console/model_parser.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index 3ba634ff1..4f899f714 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -295,9 +295,8 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur rural_prevalence; for (const auto &age_group : opt["RuralPrevalence"]) { auto age_group_range = age_group["AgeRange"].get(); - rural_prevalence[age_group_range] = { - {hgps::core::Gender::female, age_group["Female"]}, - {hgps::core::Gender::male, age_group["Male"]}}; + rural_prevalence[age_group_range] = {{hgps::core::Gender::female, age_group["Female"]}, + {hgps::core::Gender::male, age_group["Male"]}}; } // Load M/F average heights for age. From b4d72ab1d9e109fdaf9bf1c3870d8203be6c0d94 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Mon, 23 Oct 2023 16:35:09 +0100 Subject: [PATCH 12/29] Add Income model parameters to KevinHall config JSON. --- example_new/KevinHall.json | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/example_new/KevinHall.json b/example_new/KevinHall.json index ba860f21c..78989ab85 100644 --- a/example_new/KevinHall.json +++ b/example_new/KevinHall.json @@ -79,6 +79,26 @@ "Male": 0.6 } ], + "IncomeModel": [ + { + "Category": 2, + "Intercept": 0.75498278113169, + "Coefficients": { + "Sex": 0.0204338261883629, + "Over18": 0.0486699373325097, + "Sector": -0.706477682886734 + } + }, + { + "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, From d7f5b61453c9e2edf1ba3a32377782bf06291c1d Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Mon, 23 Oct 2023 17:09:52 +0100 Subject: [PATCH 13/29] Load income models from JSON in mdoel parser. --- example_new/KevinHall.json | 6 +++--- src/HealthGPS.Console/model_parser.cpp | 9 +++++++++ 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/example_new/KevinHall.json b/example_new/KevinHall.json index 78989ab85..6dce89474 100644 --- a/example_new/KevinHall.json +++ b/example_new/KevinHall.json @@ -79,9 +79,9 @@ "Male": 0.6 } ], - "IncomeModel": [ + "IncomeModels": [ { - "Category": 2, + "Name": "Category_2", "Intercept": 0.75498278113169, "Coefficients": { "Sex": 0.0204338261883629, @@ -90,7 +90,7 @@ } }, { - "Category": 3, + "Name": "Category_3", "Intercept": 1.48856946873517, "Coefficients": { "Sex": 0.000641255498596025, diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index 4f899f714..eac7e6b1b 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -299,6 +299,15 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur {hgps::core::Gender::male, age_group["Male"]}}; } + // Income models for different income classifications. + hgps::LinearModelParams income_models; + for (const auto &factor : opt["IncomeModels"]) { + auto income_class_key = factor["Name"].get(); + income_models.intercepts[income_class_key] = factor["Intercept"].get(); + income_models.coefficients[income_class_key] = + factor["Coefficients"].get>(); + } + // Load M/F average heights for age. std::unordered_map> age_mean_height; const auto max_age = static_cast(config.settings.age_range.upper()); From 7c0846fcb929364a0ee3ddbb77e74b80b363023f Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Mon, 23 Oct 2023 17:34:21 +0100 Subject: [PATCH 14/29] Load income models from mdoel parser into kevin hall model. --- src/HealthGPS.Console/model_parser.cpp | 3 ++- src/HealthGPS/kevin_hall_model.cpp | 20 +++++++++++++++++--- src/HealthGPS/kevin_hall_model.h | 9 +++++++++ 3 files changed, 28 insertions(+), 4 deletions(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index eac7e6b1b..a1996d465 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -324,7 +324,8 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur return std::make_unique( std::move(energy_equation), std::move(nutrient_ranges), std::move(nutrient_equations), - std::move(food_prices), std::move(rural_prevalence), std::move(age_mean_height)); + std::move(food_prices), std::move(rural_prevalence), std::move(income_models), + std::move(age_mean_height)); } std::pair> diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index f211f78fa..f6f285913 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -35,10 +35,12 @@ KevinHallModel::KevinHallModel( const std::unordered_map> &food_prices, const std::map> &rural_prevalence, + const LinearModelParams &income_models, const std::unordered_map> &age_mean_height) : energy_equation_{energy_equation}, nutrient_ranges_{nutrient_ranges}, nutrient_equations_{nutrient_equations}, food_prices_{food_prices}, - rural_prevalence_{rural_prevalence}, age_mean_height_{age_mean_height} { + rural_prevalence_{rural_prevalence}, income_models_{income_models}, + age_mean_height_{age_mean_height} { if (energy_equation_.empty()) { throw core::HgpsException("Energy equation mapping is empty"); @@ -52,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_.intercepts.empty() || income_models_.coefficients.empty()) { + throw core::HgpsException("Income models mapping is incomplete"); + } if (age_mean_height_.empty()) { throw core::HgpsException("Age mean height mapping is empty"); } @@ -316,10 +324,12 @@ KevinHallModelDefinition::KevinHallModelDefinition( std::unordered_map> food_prices, std::map> rural_prevalence, + LinearModelParams income_models, std::unordered_map> 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)}, age_mean_height_{std::move(age_mean_height)} { + rural_prevalence_{std::move(rural_prevalence)}, income_models_{std::move(income_models)}, + age_mean_height_{std::move(age_mean_height)} { if (energy_equation_.empty()) { throw core::HgpsException("Energy equation mapping is empty"); @@ -336,6 +346,9 @@ KevinHallModelDefinition::KevinHallModelDefinition( if (rural_prevalence_.empty()) { throw core::HgpsException("Rural prevalence mapping is empty"); } + if (income_models_.intercepts.empty() || income_models_.coefficients.empty()) { + throw core::HgpsException("Income models mapping is incomplete"); + } if (age_mean_height_.empty()) { throw core::HgpsException("Age mean height mapping is empty"); } @@ -343,7 +356,8 @@ KevinHallModelDefinition::KevinHallModelDefinition( std::unique_ptr KevinHallModelDefinition::create_model() const { return std::make_unique(energy_equation_, nutrient_ranges_, nutrient_equations_, - food_prices_, rural_prevalence_, age_mean_height_); + food_prices_, rural_prevalence_, income_models_, + age_mean_height_); } } // namespace hgps diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index e8dfd99ec..dc1d2ce70 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -1,5 +1,8 @@ #pragma once +// TODO: LinearModelParams (in static_linear_model.h) should be moved somewhere better. +#include "static_linear_model.h" + #include "interfaces.h" #include "mapping.h" @@ -54,6 +57,7 @@ class KevinHallModel final : public RiskFactorModel { /// @param nutrient_equations The nutrient coefficients for each food group /// @param food_prices The unit price for each food group /// @param rural_prevalence Rural sector prevalence for age groups and sex + /// @param income_models The income models for each income category /// @param age_mean_height The mean height at all ages (male and female) KevinHallModel( const std::unordered_map &energy_equation, @@ -63,6 +67,7 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &food_prices, const std::map> &rural_prevalence, + const LinearModelParams &income_models, const std::unordered_map> &age_mean_height); RiskFactorModelType type() const noexcept override; @@ -81,6 +86,7 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &food_prices_; const std::map> &rural_prevalence_; + const LinearModelParams &income_models_; const std::unordered_map> &age_mean_height_; // Model parameters. @@ -174,6 +180,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { /// @param nutrient_equations The nutrient coefficients for each food group /// @param food_prices The unit price for each food group /// @param rural_prevalence Rural sector prevalence for age groups and sex + /// @param income_models The income models for each income category /// @param age_mean_height The mean height at all ages (male and female) /// @throws std::invalid_argument for empty arguments KevinHallModelDefinition( @@ -183,6 +190,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> food_prices, std::map> rural_prevalence, + LinearModelParams income_models, std::unordered_map> age_mean_height); /// @brief Construct a new KevinHallModel from this definition @@ -196,6 +204,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> food_prices_; std::map> rural_prevalence_; + LinearModelParams income_models_; std::unordered_map> age_mean_height_; }; From 1aea169fdb350122355fbbf7b5b27cffa9ff05d2 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Tue, 24 Oct 2023 13:40:22 +0100 Subject: [PATCH 15/29] Add missing category_1 income model. --- example_new/KevinHall.json | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/example_new/KevinHall.json b/example_new/KevinHall.json index 6dce89474..9bea36af3 100644 --- a/example_new/KevinHall.json +++ b/example_new/KevinHall.json @@ -80,6 +80,11 @@ } ], "IncomeModels": [ + { + "Name": "Category_1", + "Intercept": 1.0, + "Coefficients": {} + }, { "Name": "Category_2", "Intercept": 0.75498278113169, From 3c6ad90a8741f03815612dfae8ed4b51f7f61e00 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Tue, 24 Oct 2023 16:38:10 +0100 Subject: [PATCH 16/29] LinearModelParams should be a vector. --- src/HealthGPS.Console/model_parser.cpp | 63 +++++++++++++++++--------- src/HealthGPS/kevin_hall_model.cpp | 26 ++++++++--- src/HealthGPS/kevin_hall_model.h | 12 +++-- src/HealthGPS/static_linear_model.cpp | 42 ++++++----------- src/HealthGPS/static_linear_model.h | 20 ++++---- 5 files changed, 92 insertions(+), 71 deletions(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index a1996d465..f5d34e6d5 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -129,31 +129,50 @@ std::unique_ptr 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(); - models.intercepts[factor_key] = factor["Intercept"].get(); - models.coefficients[factor_key] = - factor["Coefficients"].get>(); - } - - // Risk factor names and correlation matrix. - std::vector names; + // Risk factor linear models and correlation matrix. + std::vector 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(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(); + model.intercept = factor["Intercept"].get(); + model.coefficients = + factor["Coefficients"].get>(); + 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(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{correlations}.matrixL()}; - return std::make_unique(std::move(names), std::move(models), + return std::make_unique(std::move(risk_factor_models), std::move(cholesky)); } @@ -300,12 +319,14 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur } // Income models for different income classifications. - hgps::LinearModelParams income_models; + std::vector income_models; for (const auto &factor : opt["IncomeModels"]) { - auto income_class_key = factor["Name"].get(); - income_models.intercepts[income_class_key] = factor["Intercept"].get(); - income_models.coefficients[income_class_key] = + hgps::LinearModelParams model; + model.name = factor["Name"].get(); + model.intercept = factor["Intercept"].get(); + model.coefficients = factor["Coefficients"].get>(); + income_models.emplace_back(std::move(model)); } // Load M/F average heights for age. diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index f6f285913..ea88a9b6c 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -35,7 +35,7 @@ KevinHallModel::KevinHallModel( const std::unordered_map> &food_prices, const std::map> &rural_prevalence, - const LinearModelParams &income_models, + const std::vector &income_models, const std::unordered_map> &age_mean_height) : energy_equation_{energy_equation}, nutrient_ranges_{nutrient_ranges}, nutrient_equations_{nutrient_equations}, food_prices_{food_prices}, @@ -57,8 +57,8 @@ KevinHallModel::KevinHallModel( if (rural_prevalence_.empty()) { throw core::HgpsException("Rural prevalence mapping is empty"); } - if (income_models_.intercepts.empty() || income_models_.coefficients.empty()) { - throw core::HgpsException("Income models mapping is incomplete"); + 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"); @@ -156,6 +156,20 @@ void KevinHallModel::initialise_sector(RuntimeContext &context, Person &person) person.sector = sector; } +// void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) const { +// auto logits = std::vector(income_models_.intercepts.size(), 0.0); + +// // Get income model for person's sector. +// const auto &income_model = income_models_.coefficients.at(person.sector); + +// // Compute income category. +// double income = income_models_.intercepts.at(person.sector); +// for (const auto &[factor_name, coefficient] : income_model) { +// income += coefficient * person.get_risk_factor_value(factor_name); +// } + +// } + SimulatePersonState KevinHallModel::simulate_person(Person &person, double shift) const { // Initial simulated person state. const double H_0 = person.get_risk_factor_value(H_key); @@ -324,7 +338,7 @@ KevinHallModelDefinition::KevinHallModelDefinition( std::unordered_map> food_prices, std::map> rural_prevalence, - LinearModelParams income_models, + std::vector income_models, std::unordered_map> 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)}, @@ -346,8 +360,8 @@ KevinHallModelDefinition::KevinHallModelDefinition( if (rural_prevalence_.empty()) { throw core::HgpsException("Rural prevalence mapping is empty"); } - if (income_models_.intercepts.empty() || income_models_.coefficients.empty()) { - throw core::HgpsException("Income models mapping is incomplete"); + 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"); diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index dc1d2ce70..d6f92665d 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -67,7 +67,7 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &food_prices, const std::map> &rural_prevalence, - const LinearModelParams &income_models, + const std::vector &income_models, const std::unordered_map> &age_mean_height); RiskFactorModelType type() const noexcept override; @@ -86,7 +86,7 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &food_prices_; const std::map> &rural_prevalence_; - const LinearModelParams &income_models_; + const std::vector &income_models_; const std::unordered_map> &age_mean_height_; // Model parameters. @@ -106,6 +106,10 @@ class KevinHallModel final : public RiskFactorModel { /// @param person The person to initialise sector for void initialise_sector(RuntimeContext &context, Person &person) const; + // /// @brief Initialise the income category of a person + // /// @param person The person to initialise sector for + // void initialise_income(RuntimeContext &context, Person &person) const; + /// @brief Simulates the energy balance model for a given person /// @param person The person to simulate /// @param shift Model adjustment term @@ -190,7 +194,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> food_prices, std::map> rural_prevalence, - LinearModelParams income_models, + std::vector income_models, std::unordered_map> age_mean_height); /// @brief Construct a new KevinHallModel from this definition @@ -204,7 +208,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> food_prices_; std::map> rural_prevalence_; - LinearModelParams income_models_; + std::vector income_models_; std::unordered_map> age_mean_height_; }; diff --git a/src/HealthGPS/static_linear_model.cpp b/src/HealthGPS/static_linear_model.cpp index 51563481a..e6b7ff031 100644 --- a/src/HealthGPS/static_linear_model.cpp +++ b/src/HealthGPS/static_linear_model.cpp @@ -6,18 +6,13 @@ namespace hgps { -StaticLinearModel::StaticLinearModel(std::vector risk_factor_names, - LinearModelParams risk_factor_models, +StaticLinearModel::StaticLinearModel(std::vector risk_factor_models, Eigen::MatrixXd risk_factor_cholesky) - : risk_factor_names_{std::move(risk_factor_names)}, - risk_factor_models_{std::move(risk_factor_models)}, + : risk_factor_models_{std::move(risk_factor_models)}, risk_factor_cholesky_{std::move(risk_factor_cholesky)} { - if (risk_factor_names_.empty()) { - throw core::HgpsException("Risk factor names list is empty"); - } - if (risk_factor_models_.intercepts.empty() || risk_factor_models_.coefficients.empty()) { - throw core::HgpsException("Risk factor models mapping is incomplete"); + if (risk_factor_models_.empty()) { + throw core::HgpsException("Risk factor model list is empty"); } if (!risk_factor_cholesky_.allFinite()) { throw core::HgpsException("Risk factor Cholesky matrix contains non-finite values"); @@ -64,7 +59,7 @@ void StaticLinearModel::update_risk_factors(RuntimeContext &context) { Eigen::VectorXd StaticLinearModel::correlated_samples(RuntimeContext &context) { // Correlated samples using Cholesky decomposition. - Eigen::VectorXd samples{risk_factor_names_.size()}; + Eigen::VectorXd samples{risk_factor_models_.size()}; std::ranges::generate(samples, [&context] { return context.random().next_normal(0.0, 1.0); }); samples = risk_factor_cholesky_ * samples; @@ -80,30 +75,22 @@ Eigen::VectorXd StaticLinearModel::correlated_samples(RuntimeContext &context) { void StaticLinearModel::linear_approximation(Person &person) { // Approximate risk factor values for person with linear models. - for (const auto &factor_name : risk_factor_names_) { - double factor = risk_factor_models_.intercepts.at(factor_name); - const auto &coefficients = risk_factor_models_.coefficients.at(factor_name); - - for (const auto &[coefficient_name, coefficient_value] : coefficients) { + for (const auto &model : risk_factor_models_) { + double factor = model.intercept; + for (const auto &[coefficient_name, coefficient_value] : model.coefficients) { factor += coefficient_value * person.get_risk_factor_value(coefficient_name); } - - person.risk_factors[factor_name] = factor; + person.risk_factors[model.name] = factor; } } StaticLinearModelDefinition::StaticLinearModelDefinition( - std::vector risk_factor_names, LinearModelParams risk_factor_models, - Eigen::MatrixXd risk_factor_cholesky) - : risk_factor_names_{std::move(risk_factor_names)}, - risk_factor_models_{std::move(risk_factor_models)}, + std::vector risk_factor_models, Eigen::MatrixXd risk_factor_cholesky) + : risk_factor_models_{std::move(risk_factor_models)}, risk_factor_cholesky_{std::move(risk_factor_cholesky)} { - if (risk_factor_names_.empty()) { - throw core::HgpsException("Risk factor names list is empty"); - } - if (risk_factor_models_.intercepts.empty() || risk_factor_models_.coefficients.empty()) { - throw core::HgpsException("Risk factor models mapping is incomplete"); + if (risk_factor_models_.empty()) { + throw core::HgpsException("Risk factor model list is empty"); } if (!risk_factor_cholesky_.allFinite()) { throw core::HgpsException("Risk factor Cholesky matrix contains non-finite values"); @@ -111,8 +98,7 @@ StaticLinearModelDefinition::StaticLinearModelDefinition( } std::unique_ptr StaticLinearModelDefinition::create_model() const { - return std::make_unique(risk_factor_names_, risk_factor_models_, - risk_factor_cholesky_); + return std::make_unique(risk_factor_models_, risk_factor_cholesky_); } } // namespace hgps diff --git a/src/HealthGPS/static_linear_model.h b/src/HealthGPS/static_linear_model.h index 63dc4667f..3928dba86 100644 --- a/src/HealthGPS/static_linear_model.h +++ b/src/HealthGPS/static_linear_model.h @@ -9,8 +9,9 @@ namespace hgps { /// @brief Defines the linear model parameters used to initialise risk factors struct LinearModelParams { - std::unordered_map intercepts; - std::unordered_map> coefficients; + core::Identifier name; + double intercept; + std::unordered_map coefficients; }; /// @brief Implements the static linear model type @@ -19,12 +20,11 @@ struct LinearModelParams { class StaticLinearModel final : public RiskFactorModel { public: /// @brief Initialises a new instance of the StaticLinearModel class - /// @param risk_factor_names An ordered list of risk factor names /// @param risk_factor_models The linear models used to initialise a person's risk factor values /// @param risk_factor_cholesky The Cholesky decomposition of the risk factor correlation matrix /// @throws HgpsException for invalid arguments - StaticLinearModel(std::vector risk_factor_names, - LinearModelParams risk_factor_models, Eigen::MatrixXd risk_factor_cholesky); + StaticLinearModel(std::vector risk_factor_models, + Eigen::MatrixXd risk_factor_cholesky); RiskFactorModelType type() const noexcept override; @@ -39,8 +39,7 @@ class StaticLinearModel final : public RiskFactorModel { void linear_approximation(Person &person); private: - const std::vector risk_factor_names_; - const LinearModelParams risk_factor_models_; + const std::vector risk_factor_models_; const Eigen::MatrixXd risk_factor_cholesky_; }; @@ -48,12 +47,10 @@ class StaticLinearModel final : public RiskFactorModel { class StaticLinearModelDefinition final : public RiskFactorModelDefinition { public: /// @brief Initialises a new instance of the StaticLinearModelDefinition class - /// @param risk_factor_names An ordered list of risk factor names /// @param risk_factor_models The linear models used to initialise a person's risk factor values /// @param risk_factor_cholesky The Cholesky decomposition of the risk factor correlation matrix /// @throws HgpsException for invalid arguments - StaticLinearModelDefinition(std::vector risk_factor_names, - LinearModelParams risk_factor_models, + StaticLinearModelDefinition(std::vector risk_factor_models, Eigen::MatrixXd risk_factor_cholesky); /// @brief Construct a new StaticLinearModel from this definition @@ -61,8 +58,7 @@ class StaticLinearModelDefinition final : public RiskFactorModelDefinition { std::unique_ptr create_model() const override; private: - std::vector risk_factor_names_; - LinearModelParams risk_factor_models_; + std::vector risk_factor_models_; Eigen::MatrixXd risk_factor_cholesky_; }; From c72c2bc5116c698b3ae0ed8b7840b4719ab05678 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Tue, 24 Oct 2023 22:56:28 +0100 Subject: [PATCH 17/29] Add initialise_income to kevin hall model. --- src/HealthGPS/kevin_hall_model.cpp | 44 +++++++++++++++++++++++------- src/HealthGPS/kevin_hall_model.h | 7 +++-- src/HealthGPS/person.h | 3 ++ 3 files changed, 41 insertions(+), 13 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index ea88a9b6c..17ee01b3e 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -156,19 +156,43 @@ void KevinHallModel::initialise_sector(RuntimeContext &context, Person &person) person.sector = sector; } -// void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) const { -// auto logits = std::vector(income_models_.intercepts.size(), 0.0); +void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) const { + + // Compute logits for each income category. + auto logits = std::vector(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); + } + } -// // Get income model for person's sector. -// const auto &income_model = income_models_.coefficients.at(person.sector); + // Compute softmax probabilities for each income category. + auto e_logits = std::vector(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. -// double income = income_models_.intercepts.at(person.sector); -// for (const auto &[factor_name, coefficient] : income_model) { -// income += coefficient * person.get_risk_factor_value(factor_name); -// } + // Compute income category probabilities. + auto probabilities = std::vector(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. diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index d6f92665d..57f8cdd36 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -106,9 +106,10 @@ class KevinHallModel final : public RiskFactorModel { /// @param person The person to initialise sector for void initialise_sector(RuntimeContext &context, Person &person) const; - // /// @brief Initialise the income category of a person - // /// @param person The person to initialise sector for - // void initialise_income(RuntimeContext &context, Person &person) const; + /// @brief Initialise the income category of a person + /// @param context The runtime context + /// @param person The person to initialise sector for + void initialise_income(RuntimeContext &context, Person &person) const; /// @brief Simulates the energy balance model for a given person /// @param person The person to simulate diff --git a/src/HealthGPS/person.h b/src/HealthGPS/person.h index 8cc1fd5f9..9287c44c4 100644 --- a/src/HealthGPS/person.h +++ b/src/HealthGPS/person.h @@ -62,6 +62,9 @@ struct Person { /// @brief Social-economic status (SES) assigned value double ses{}; + /// @brief Income category + core::Identifier income{}; + /// @brief Current risk factors values std::map risk_factors; From 61ab5b78b8cce4afc70752ad1d10d373a8500abf Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Wed, 25 Oct 2023 11:32:55 +0100 Subject: [PATCH 18/29] Fix read after std::move. --- src/HealthGPS.Console/model_parser.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index f5d34e6d5..c75308116 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -144,12 +144,6 @@ load_staticlinear_risk_model_definition(const poco::json &opt, const host::Confi model.intercept = factor["Intercept"].get(); model.coefficients = factor["Coefficients"].get>(); - 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(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()}; @@ -159,6 +153,12 @@ load_staticlinear_risk_model_definition(const poco::json &opt, const host::Confi "column {} name ({})", i, model.name.to_string(), i, column_name.to_string())}; } + + // Write data structures. + for (size_t j = 0; j < correlations_table.num_rows(); j++) { + correlations(i, j) = std::any_cast(correlations_table.column(i).value(j)); + } + risk_factor_models.emplace_back(std::move(model)); } // Check correlation matrix column count matches risk factor count. From 06e553e955fe97ffcf94bc247c51fd9ecabeb0af Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Wed, 25 Oct 2023 12:47:28 +0100 Subject: [PATCH 19/29] Call initialise_income in KH model. --- src/HealthGPS/kevin_hall_model.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 17ee01b3e..281a9c52a 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -71,17 +71,17 @@ std::string KevinHallModel::name() const noexcept { return "Dynamic"; } void KevinHallModel::generate_risk_factors(RuntimeContext &context) { - // Initialise sector for everyone. + // Initialise sector and income for everyone. for (auto &person : context.population()) { initialise_sector(context, person); + initialise_income(context, person); } } void KevinHallModel::update_risk_factors(RuntimeContext &context) { - hgps::Population &population = context.population(); - // Initialise sector for newborns. - for (auto &person : population) { + // Initialise sector and income for newborns. + for (auto &person : context.population()) { // Do not initialise non-newborns. if (person.age > 0) { @@ -89,6 +89,7 @@ void KevinHallModel::update_risk_factors(RuntimeContext &context) { } initialise_sector(context, person); + initialise_income(context, person); } // TODO: Compute target body weight. @@ -97,7 +98,7 @@ void KevinHallModel::update_risk_factors(RuntimeContext &context) { // Trial run. double mean_sim_body_weight = 0.0; double mean_adjustment_coefficient = 0.0; - for (auto &person : population) { + for (auto &person : context.population()) { // Ignore if inactive. if (!person.is_active()) { continue; @@ -110,13 +111,13 @@ void KevinHallModel::update_risk_factors(RuntimeContext &context) { } // Compute model adjustment term. - const size_t population_size = population.current_active_size(); + const size_t population_size = context.population().current_active_size(); mean_sim_body_weight /= population_size; mean_adjustment_coefficient /= population_size; double shift = (target_BW - mean_sim_body_weight) / mean_adjustment_coefficient; // Final run. - for (auto &person : population) { + for (auto &person : context.population()) { // Ignore if inactive. if (!person.is_active()) { continue; From e373e53b49118d075ca1c1372f266dbbf71f3e6e Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Wed, 25 Oct 2023 14:40:49 +0100 Subject: [PATCH 20/29] Add sector update code to kevinhall. --- example_new/KevinHall.json | 8 +++--- src/HealthGPS.Console/model_parser.cpp | 8 +++--- src/HealthGPS/kevin_hall_model.cpp | 38 +++++++++++++++++++------- src/HealthGPS/kevin_hall_model.h | 13 ++++++--- 4 files changed, 45 insertions(+), 22 deletions(-) diff --git a/example_new/KevinHall.json b/example_new/KevinHall.json index 9bea36af3..a0f1a8a9d 100644 --- a/example_new/KevinHall.json +++ b/example_new/KevinHall.json @@ -69,12 +69,12 @@ }, "RuralPrevalence": [ { - "AgeRange": [0, 17], + "Name": "Under18", "Female": 0.65, "Male": 0.65 }, { - "AgeRange": [18, 100], + "Name": "From18", "Female": 0.6, "Male": 0.6 } @@ -89,7 +89,7 @@ "Name": "Category_2", "Intercept": 0.75498278113169, "Coefficients": { - "Sex": 0.0204338261883629, + "Gender": 0.0204338261883629, "Over18": 0.0486699373325097, "Sector": -0.706477682886734 } @@ -98,7 +98,7 @@ "Name": "Category_3", "Intercept": 1.48856946873517, "Coefficients": { - "Sex": 0.000641255498596025, + "Gender": 0.000641255498596025, "Over18": 0.206622749570132, "Sector": -1.71313287940798 } diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index c75308116..9f402839e 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -310,12 +310,12 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur const auto food_data_table = load_datatable_from_csv(food_data_file_info); // Rural sector prevalence for age groups and sex. - std::map> + std::map> rural_prevalence; for (const auto &age_group : opt["RuralPrevalence"]) { - auto age_group_range = age_group["AgeRange"].get(); - rural_prevalence[age_group_range] = {{hgps::core::Gender::female, age_group["Female"]}, - {hgps::core::Gender::male, age_group["Male"]}}; + auto age_group_name = age_group["Name"].get(); + rural_prevalence[age_group_name] = {{hgps::core::Gender::female, age_group["Female"]}, + {hgps::core::Gender::male, age_group["Male"]}}; } // Income models for different income classifications. diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 281a9c52a..08dafd64c 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -33,7 +33,7 @@ KevinHallModel::KevinHallModel( const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, - const std::map> + const std::map> &rural_prevalence, const std::vector &income_models, const std::unordered_map> &age_mean_height) @@ -142,21 +142,39 @@ 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; - } + // Get rural prevalence for age group and sex. + double prevalence; + if (person.age < 18) { + prevalence = rural_prevalence_.at("Under18"_id).at(person.gender); + } else { + prevalence = rural_prevalence_.at("From18"_id).at(person.gender); } // Sample the person's sector. double rand = context.random().next_double(); - auto sector = rand < rural_prevalence ? core::Sector::rural : core::Sector::urban; + auto sector = rand < prevalence ? core::Sector::rural : core::Sector::urban; person.sector = sector; } +void KevinHallModel::update_sector(RuntimeContext &context, Person &person) const { + + // Only update rural sector 18 year olds. + if ((person.age != 18) || (person.sector != core::Sector::rural)) { + return; + } + + // Get rural prevalence for age group and sex. + double prevalence_under18 = rural_prevalence_.at("Under18"_id).at(person.gender); + double prevalence_from18 = rural_prevalence_.at("From18"_id).at(person.gender); + + // Compute random rural to urban transition. + double rand = context.random().next_double(); + double p_rural_to_urban = 1.0 - prevalence_from18 / prevalence_under18; + if (rand < p_rural_to_urban) { + person.sector = core::Sector::urban; + } +} + void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) const { // Compute logits for each income category. @@ -361,7 +379,7 @@ KevinHallModelDefinition::KevinHallModelDefinition( std::unordered_map> nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, - std::map> + std::map> rural_prevalence, std::vector income_models, std::unordered_map> age_mean_height) diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index 57f8cdd36..986613acb 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -65,7 +65,7 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, - const std::map> + const std::map> &rural_prevalence, const std::vector &income_models, const std::unordered_map> &age_mean_height); @@ -84,7 +84,7 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &nutrient_equations_; const std::unordered_map> &food_prices_; - const std::map> + const std::map> &rural_prevalence_; const std::vector &income_models_; const std::unordered_map> &age_mean_height_; @@ -106,6 +106,11 @@ class KevinHallModel final : public RiskFactorModel { /// @param person The person to initialise sector for void initialise_sector(RuntimeContext &context, Person &person) const; + /// @brief Update the sector of a person + /// @param context The runtime context + /// @param person The person to update sector for + void update_sector(RuntimeContext &context, Person &person) const; + /// @brief Initialise the income category of a person /// @param context The runtime context /// @param person The person to initialise sector for @@ -193,7 +198,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, - std::map> + std::map> rural_prevalence, std::vector income_models, std::unordered_map> age_mean_height); @@ -207,7 +212,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map> nutrient_ranges_; std::unordered_map> nutrient_equations_; std::unordered_map> food_prices_; - std::map> + std::map> rural_prevalence_; std::vector income_models_; std::unordered_map> age_mean_height_; From f5a22d29fdbbf3469702d9434c87b715524a486b Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Wed, 25 Oct 2023 15:40:37 +0100 Subject: [PATCH 21/29] Code to convert sector to a real, for use as a coefficient. Add error checking for unknown gender and sector. --- src/HealthGPS/person.cpp | 20 ++++++++++++++++++-- src/HealthGPS/person.h | 11 +++++++++-- 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/src/HealthGPS/person.cpp b/src/HealthGPS/person.cpp index 4a0690195..e77f7482e 100644 --- a/src/HealthGPS/person.cpp +++ b/src/HealthGPS/person.cpp @@ -1,5 +1,7 @@ #include "person.h" +#include "HealthGPS.Core/exception.h" + namespace hgps { std::atomic Person::newUID{0}; @@ -10,6 +12,7 @@ std::map> Person::curren {"Age"_id, [](const Person &p) { return static_cast(p.age); }}, {"Age2"_id, [](const Person &p) { return pow(p.age, 2); }}, {"Age3"_id, [](const Person &p) { return pow(p.age, 3); }}, + {"Sector"_id, [](const Person &p) { return p.sector_to_value(); }}, {"SES"_id, [](const Person &p) { return p.ses; }}, // HACK: ew, gross... allows us to mock risk factors we don't have data for yet @@ -56,14 +59,27 @@ double Person::get_risk_factor_value(const core::Identifier &key) const { throw std::out_of_range("Risk factor not found: " + key.to_string()); } -float Person::gender_to_value() const noexcept { +float Person::gender_to_value() const { + if (gender == core::Gender::unknown) { + throw core::HgpsException("Gender is unknown."); + } return gender == core::Gender::male ? 1.0f : 0.0f; } -std::string Person::gender_to_string() const noexcept { +std::string Person::gender_to_string() const { + if (gender == core::Gender::unknown) { + throw core::HgpsException("Gender is unknown."); + } return gender == core::Gender::male ? "male" : "female"; } +float Person::sector_to_value() const { + if (sector == core::Sector::unknown) { + throw core::HgpsException("Sector is unknown."); + } + return sector == core::Sector::urban ? 0.0f : 1.0f; +} + void Person::emigrate(const unsigned int time) { if (!is_active()) { throw std::logic_error("Entity must be active prior to emigrate."); diff --git a/src/HealthGPS/person.h b/src/HealthGPS/person.h index 9287c44c4..2c2b75677 100644 --- a/src/HealthGPS/person.h +++ b/src/HealthGPS/person.h @@ -100,11 +100,18 @@ struct Person { /// @brief Gets the gender enumeration as a number for analysis /// @return The gender associated value - float gender_to_value() const noexcept; + /// @throws HgpsException if gender is unknown + float gender_to_value() const; /// @brief Gets the gender enumeration name string /// @return The gender name - std::string gender_to_string() const noexcept; + /// @throws HgpsException if gender is unknown + std::string gender_to_string() const; + + /// @brief Gets the sector enumeration as a number + /// @return The sector value (0 for urban, 1 for rural) + /// @throws HgpsException if sector is unknown + float sector_to_value() const; /// @brief Emigrate this instance from the virtual population /// @param time Migration time From 2326ca9898de9d70130eae7d6e0549bc89fb6058 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Wed, 25 Oct 2023 16:02:59 +0100 Subject: [PATCH 22/29] Age group name consistency. --- example_new/KevinHall.json | 2 +- src/HealthGPS/kevin_hall_model.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/example_new/KevinHall.json b/example_new/KevinHall.json index a0f1a8a9d..c49ecc9d5 100644 --- a/example_new/KevinHall.json +++ b/example_new/KevinHall.json @@ -74,7 +74,7 @@ "Male": 0.65 }, { - "Name": "From18", + "Name": "Over18", "Female": 0.6, "Male": 0.6 } diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 08dafd64c..2dacb164e 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -147,7 +147,7 @@ void KevinHallModel::initialise_sector(RuntimeContext &context, Person &person) if (person.age < 18) { prevalence = rural_prevalence_.at("Under18"_id).at(person.gender); } else { - prevalence = rural_prevalence_.at("From18"_id).at(person.gender); + prevalence = rural_prevalence_.at("Over18"_id).at(person.gender); } // Sample the person's sector. @@ -165,11 +165,11 @@ void KevinHallModel::update_sector(RuntimeContext &context, Person &person) cons // Get rural prevalence for age group and sex. double prevalence_under18 = rural_prevalence_.at("Under18"_id).at(person.gender); - double prevalence_from18 = rural_prevalence_.at("From18"_id).at(person.gender); + double prevalence_over18 = rural_prevalence_.at("Over18"_id).at(person.gender); // Compute random rural to urban transition. double rand = context.random().next_double(); - double p_rural_to_urban = 1.0 - prevalence_from18 / prevalence_under18; + double p_rural_to_urban = 1.0 - prevalence_over18 / prevalence_under18; if (rand < p_rural_to_urban) { person.sector = core::Sector::urban; } From 3985efa42b14af1d860047bffcc8ad52ee6f61ac Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Wed, 25 Oct 2023 16:16:32 +0100 Subject: [PATCH 23/29] Add person method to check over 18. --- src/HealthGPS/person.cpp | 3 +++ src/HealthGPS/person.h | 4 ++++ 2 files changed, 7 insertions(+) diff --git a/src/HealthGPS/person.cpp b/src/HealthGPS/person.cpp index e77f7482e..ee511357c 100644 --- a/src/HealthGPS/person.cpp +++ b/src/HealthGPS/person.cpp @@ -12,6 +12,7 @@ std::map> Person::curren {"Age"_id, [](const Person &p) { return static_cast(p.age); }}, {"Age2"_id, [](const Person &p) { return pow(p.age, 2); }}, {"Age3"_id, [](const Person &p) { return pow(p.age, 3); }}, + {"Over18"_id, [](const Person &p) { return static_cast(p.over_18()); }}, {"Sector"_id, [](const Person &p) { return p.sector_to_value(); }}, {"SES"_id, [](const Person &p) { return p.ses; }}, @@ -80,6 +81,8 @@ float Person::sector_to_value() const { return sector == core::Sector::urban ? 0.0f : 1.0f; } +bool Person::over_18() const noexcept { return age >= 18; } + void Person::emigrate(const unsigned int time) { if (!is_active()) { throw std::logic_error("Entity must be active prior to emigrate."); diff --git a/src/HealthGPS/person.h b/src/HealthGPS/person.h index 2c2b75677..d7d89db96 100644 --- a/src/HealthGPS/person.h +++ b/src/HealthGPS/person.h @@ -113,6 +113,10 @@ struct Person { /// @throws HgpsException if sector is unknown float sector_to_value() const; + /// @brief Check if person is an adult (18 or over) + /// @return true if person is 18 or over; else false + bool over_18() const noexcept; + /// @brief Emigrate this instance from the virtual population /// @param time Migration time /// @throws std::logic_error for attempting to set to non-active individuals. From 905223abfe17c1777372acf7678f09b3dfda1fe4 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Wed, 25 Oct 2023 16:37:11 +0100 Subject: [PATCH 24/29] Add kevinhall code to update income. --- src/HealthGPS/kevin_hall_model.cpp | 27 ++++++++++++++++++--------- src/HealthGPS/kevin_hall_model.h | 5 +++++ 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 2dacb164e..84d0de412 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -71,7 +71,7 @@ std::string KevinHallModel::name() const noexcept { return "Dynamic"; } void KevinHallModel::generate_risk_factors(RuntimeContext &context) { - // Initialise sector and income for everyone. + // Initialise everyone. for (auto &person : context.population()) { initialise_sector(context, person); initialise_income(context, person); @@ -80,16 +80,15 @@ void KevinHallModel::generate_risk_factors(RuntimeContext &context) { void KevinHallModel::update_risk_factors(RuntimeContext &context) { - // Initialise sector and income for newborns. + // Initialise newborns and update others. for (auto &person : context.population()) { - - // Do not initialise non-newborns. - if (person.age > 0) { - continue; + if (person.age == 0) { + initialise_sector(context, person); + initialise_income(context, person); + } else { + update_sector(context, person); + update_income(context, person); } - - initialise_sector(context, person); - initialise_income(context, person); } // TODO: Compute target body weight. @@ -213,6 +212,16 @@ void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) throw core::HgpsException("Logic Error: failed to initialise income category"); } +void KevinHallModel::update_income(RuntimeContext &context, Person &person) const { + + // Only update 18 year olds. + if (person.age != 18) { + return; + } + + initialise_income(context, person); +} + SimulatePersonState KevinHallModel::simulate_person(Person &person, double shift) const { // Initial simulated person state. const double H_0 = person.get_risk_factor_value(H_key); diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index 986613acb..ddf666c0a 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -116,6 +116,11 @@ class KevinHallModel final : public RiskFactorModel { /// @param person The person to initialise sector for void initialise_income(RuntimeContext &context, Person &person) const; + /// @brief Update the income category of a person + /// @param context The runtime context + /// @param person The person to update sector for + void update_income(RuntimeContext &context, Person &person) const; + /// @brief Simulates the energy balance model for a given person /// @param person The person to simulate /// @param shift Model adjustment term From 71926bdbaa933736cb7051f89758089c26324431 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Thu, 26 Oct 2023 13:48:54 +0100 Subject: [PATCH 25/29] Change weird if Co-authored-by: Alex Dewar --- src/HealthGPS/kevin_hall_model.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 84d0de412..4fbfed3c2 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -215,11 +215,9 @@ void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) void KevinHallModel::update_income(RuntimeContext &context, Person &person) const { // Only update 18 year olds. - if (person.age != 18) { - return; + if (person.age == 18) { + initialise_income(context, person); } - - initialise_income(context, person); } SimulatePersonState KevinHallModel::simulate_person(Person &person, double shift) const { From d32fd44216d144e001cc4e6e574dda45e65839b6 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Thu, 26 Oct 2023 14:39:14 +0100 Subject: [PATCH 26/29] nutrient_ranges is now a Double interval, interval now has its own clamp method, and disallowed nonsensical intervals at construction. --- src/HealthGPS.Console/model_parser.cpp | 8 ++----- src/HealthGPS.Core/interval.h | 23 ++++++++++++++------- src/HealthGPS.Tests/AgeGenderTable.Test.cpp | 2 +- src/HealthGPS.Tests/Interval.Test.cpp | 18 ---------------- src/HealthGPS/gender_table.h | 4 ++-- src/HealthGPS/kevin_hall_model.cpp | 11 ++-------- src/HealthGPS/kevin_hall_model.h | 18 ++++++---------- 7 files changed, 28 insertions(+), 56 deletions(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index 9f402839e..ffeae8991 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -275,14 +275,10 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur // Nutrient groups. std::unordered_map energy_equation; - std::unordered_map> nutrient_ranges; + std::unordered_map nutrient_ranges; for (const auto &nutrient : opt["Nutrients"]) { auto nutrient_key = nutrient["Name"].get(); - nutrient_ranges[nutrient_key] = nutrient["Range"].get>(); - if (nutrient_ranges[nutrient_key].first > nutrient_ranges[nutrient_key].second) { - throw hgps::core::HgpsException{ - fmt::format("Nutrient range is invalid: {}", nutrient_key.to_string())}; - } + nutrient_ranges[nutrient_key] = nutrient["Range"].get(); energy_equation[nutrient_key] = nutrient["Energy"].get(); } diff --git a/src/HealthGPS.Core/interval.h b/src/HealthGPS.Core/interval.h index ef749b806..680ea315e 100644 --- a/src/HealthGPS.Core/interval.h +++ b/src/HealthGPS.Core/interval.h @@ -1,6 +1,10 @@ #pragma once + +#include "HealthGPS.Core/exception.h" #include "forward_type.h" #include "string_util.h" + +#include #include namespace hgps::core { @@ -16,7 +20,11 @@ template class Interval { /// @param lower_value Lower bound value /// @param upper_value Upper bound value explicit Interval(TYPE lower_value, TYPE upper_value) - : lower_{lower_value}, upper_{upper_value} {} + : lower_{lower_value}, upper_{upper_value} { + if (lower_ > upper_) { + throw HgpsException(fmt::format("Invalid interval: {}-{}", lower_, upper_)); + } + } /// @brief Gets the interval lower bound /// @return The lower bound value @@ -33,13 +41,7 @@ template class Interval { /// @brief Determines whether a value is in the Interval. /// @param value The value to check /// @return true if the value is in the interval; otherwise, false. - bool contains(TYPE value) const noexcept { - if (lower_ < upper_) { - return lower_ <= value && value <= upper_; - } - - return lower_ >= value && value >= upper_; - } + bool contains(TYPE value) const noexcept { return lower_ <= value && value <= upper_; } /// @brief Determines whether an Interval is inside this instance interval. /// @param other The other Interval to check @@ -48,6 +50,11 @@ template class Interval { return contains(other.lower_) && contains(other.upper_); } + /// @brief Clamp a given value to the interval boundaries + /// @param value The value to clamp + /// @return The clamped value + TYPE clamp(TYPE value) const noexcept { return std::clamp(value, lower_, upper_); } + /// @brief Convert this instance to a string representation /// @return The equivalent string representation std::string to_string() const noexcept { return fmt::format("{}-{}", lower_, upper_); } diff --git a/src/HealthGPS.Tests/AgeGenderTable.Test.cpp b/src/HealthGPS.Tests/AgeGenderTable.Test.cpp index a5114b9e7..da739652d 100644 --- a/src/HealthGPS.Tests/AgeGenderTable.Test.cpp +++ b/src/HealthGPS.Tests/AgeGenderTable.Test.cpp @@ -159,7 +159,7 @@ TEST(TestHealthGPS_AgeGenderTable, CreateWithWrongRangerThrows) { using namespace hgps; auto negative_range = core::IntegerInterval(-1, 10); - auto inverted_range = core::IntegerInterval(10, 1); + auto inverted_range = core::IntegerInterval(1, 1); ASSERT_THROW(create_age_gender_table(negative_range), std::invalid_argument); ASSERT_THROW(create_age_gender_table(inverted_range), std::invalid_argument); diff --git a/src/HealthGPS.Tests/Interval.Test.cpp b/src/HealthGPS.Tests/Interval.Test.cpp index 7f25d4488..72e7318c1 100644 --- a/src/HealthGPS.Tests/Interval.Test.cpp +++ b/src/HealthGPS.Tests/Interval.Test.cpp @@ -30,24 +30,6 @@ TEST(TestCore_Interval, CreatePositive) { ASSERT_TRUE(animal.contains(dog)); } -TEST(TestCore_Interval, CreateNegative) { - using namespace hgps::core; - auto lower = 0; - auto upper = -10; - auto len = upper - lower; - auto mid = len / 2; - auto animal = IntegerInterval{lower, upper}; - auto cat = IntegerInterval{mid, upper}; - auto dog = IntegerInterval{lower, mid}; - - ASSERT_EQ(lower, animal.lower()); - ASSERT_EQ(upper, animal.upper()); - ASSERT_EQ(len, animal.length()); - ASSERT_TRUE(animal.contains(mid)); - ASSERT_TRUE(animal.contains(cat)); - ASSERT_TRUE(animal.contains(dog)); -} - TEST(TestCore_Interval, Comparable) { using namespace hgps::core; diff --git a/src/HealthGPS/gender_table.h b/src/HealthGPS/gender_table.h index f8d9b4a57..f2046ed87 100644 --- a/src/HealthGPS/gender_table.h +++ b/src/HealthGPS/gender_table.h @@ -157,10 +157,10 @@ GenderTable create_integer_gender_table(const core::IntegerInterval & /// @tparam TYPE The values data type /// @param age_range The age breakpoints range /// @return A new instance of the AgeGenderTable class -/// @throws std::out_of_range for age range 'lower' of negative value or less than the 'upper' value +/// @throws std::out_of_range for age range 'lower' of negative value or equal to the 'upper' value template AgeGenderTable create_age_gender_table(const core::IntegerInterval &age_range) { - if (age_range.lower() < 0 || age_range.lower() >= age_range.upper()) { + if (age_range.lower() < 0 || age_range.lower() == age_range.upper()) { throw std::invalid_argument( "The 'age lower' value must be greater than zero and less than the 'age upper' value."); } diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 4fbfed3c2..43e5ef286 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -3,7 +3,6 @@ #include "HealthGPS.Core/exception.h" -#include #include /* @@ -29,7 +28,7 @@ const core::Identifier CI_key{"Carbohydrate"}; KevinHallModel::KevinHallModel( const std::unordered_map &energy_equation, - const std::unordered_map> &nutrient_ranges, + const std::unordered_map &nutrient_ranges, const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, @@ -375,15 +374,9 @@ double KevinHallModel::compute_AT(double EI, double EI_0) const { return beta_AT * delta_EI; } -double KevinHallModel::bounded_nutrient_value(const core::Identifier &nutrient, - double value) const { - const auto &range = nutrient_ranges_.at(nutrient); - return std::clamp(range.first, range.second, value); -} - KevinHallModelDefinition::KevinHallModelDefinition( std::unordered_map energy_equation, - std::unordered_map> nutrient_ranges, + std::unordered_map nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, std::map> diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index ddf666c0a..f0c683666 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -53,7 +53,7 @@ class KevinHallModel final : public RiskFactorModel { public: /// @brief Initialises a new instance of the KevinHallModel class /// @param energy_equation The energy coefficients for each nutrient - /// @param nutrient_ranges The minimum and maximum nutrient values + /// @param nutrient_ranges The interval boundaries for nutrient values /// @param nutrient_equations The nutrient coefficients for each food group /// @param food_prices The unit price for each food group /// @param rural_prevalence Rural sector prevalence for age groups and sex @@ -61,7 +61,7 @@ class KevinHallModel final : public RiskFactorModel { /// @param age_mean_height The mean height at all ages (male and female) KevinHallModel( const std::unordered_map &energy_equation, - const std::unordered_map> &nutrient_ranges, + const std::unordered_map &nutrient_ranges, const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, @@ -80,7 +80,7 @@ class KevinHallModel final : public RiskFactorModel { private: const std::unordered_map &energy_equation_; - const std::unordered_map> &nutrient_ranges_; + const std::unordered_map &nutrient_ranges_; const std::unordered_map> &nutrient_equations_; const std::unordered_map> &food_prices_; @@ -178,12 +178,6 @@ class KevinHallModel final : public RiskFactorModel { /// @param EI_0 The initial energy intake /// @return The computed adaptive thermogenesis double compute_AT(double EI, double EI_0) const; - - /// @brief Return the nutrient value bounded within its range - /// @param nutrient The nutrient Identifier - /// @param value The nutrient value to bound - /// @return The bounded nutrient value - double bounded_nutrient_value(const core::Identifier &nutrient, double value) const; }; /// @brief Defines the energy balance model data type @@ -191,7 +185,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { public: /// @brief Initialises a new instance of the KevinHallModelDefinition class /// @param energy_equation The energy coefficients for each nutrient - /// @param nutrient_ranges The minimum and maximum nutrient values + /// @param nutrient_ranges The interval boundaries for nutrient values /// @param nutrient_equations The nutrient coefficients for each food group /// @param food_prices The unit price for each food group /// @param rural_prevalence Rural sector prevalence for age groups and sex @@ -200,7 +194,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { /// @throws std::invalid_argument for empty arguments KevinHallModelDefinition( std::unordered_map energy_equation, - std::unordered_map> nutrient_ranges, + std::unordered_map nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, std::map> @@ -214,7 +208,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { private: std::unordered_map energy_equation_; - std::unordered_map> nutrient_ranges_; + std::unordered_map nutrient_ranges_; std::unordered_map> nutrient_equations_; std::unordered_map> food_prices_; std::map> From a44f19cd803d165bf38efb96e2f5e89b7c12b26c Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Thu, 26 Oct 2023 14:48:27 +0100 Subject: [PATCH 27/29] KH: rural_prevalence should be an unordered_map. --- src/HealthGPS.Console/model_parser.cpp | 2 +- src/HealthGPS/kevin_hall_model.cpp | 4 ++-- src/HealthGPS/kevin_hall_model.h | 10 +++++----- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/HealthGPS.Console/model_parser.cpp b/src/HealthGPS.Console/model_parser.cpp index ffeae8991..43bc12272 100644 --- a/src/HealthGPS.Console/model_parser.cpp +++ b/src/HealthGPS.Console/model_parser.cpp @@ -306,7 +306,7 @@ load_kevinhall_risk_model_definition(const poco::json &opt, const host::Configur const auto food_data_table = load_datatable_from_csv(food_data_file_info); // Rural sector prevalence for age groups and sex. - std::map> + std::unordered_map> rural_prevalence; for (const auto &age_group : opt["RuralPrevalence"]) { auto age_group_name = age_group["Name"].get(); diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index 43e5ef286..ba4e496f3 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -32,7 +32,7 @@ KevinHallModel::KevinHallModel( const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, - const std::map> + const std::unordered_map> &rural_prevalence, const std::vector &income_models, const std::unordered_map> &age_mean_height) @@ -379,7 +379,7 @@ KevinHallModelDefinition::KevinHallModelDefinition( std::unordered_map nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, - std::map> + std::unordered_map> rural_prevalence, std::vector income_models, std::unordered_map> age_mean_height) diff --git a/src/HealthGPS/kevin_hall_model.h b/src/HealthGPS/kevin_hall_model.h index f0c683666..182ea9ad4 100644 --- a/src/HealthGPS/kevin_hall_model.h +++ b/src/HealthGPS/kevin_hall_model.h @@ -65,8 +65,8 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &nutrient_equations, const std::unordered_map> &food_prices, - const std::map> - &rural_prevalence, + const std::unordered_map> &rural_prevalence, const std::vector &income_models, const std::unordered_map> &age_mean_height); @@ -84,7 +84,7 @@ class KevinHallModel final : public RiskFactorModel { const std::unordered_map> &nutrient_equations_; const std::unordered_map> &food_prices_; - const std::map> + const std::unordered_map> &rural_prevalence_; const std::vector &income_models_; const std::unordered_map> &age_mean_height_; @@ -197,7 +197,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map nutrient_ranges, std::unordered_map> nutrient_equations, std::unordered_map> food_prices, - std::map> + std::unordered_map> rural_prevalence, std::vector income_models, std::unordered_map> age_mean_height); @@ -211,7 +211,7 @@ class KevinHallModelDefinition final : public RiskFactorModelDefinition { std::unordered_map nutrient_ranges_; std::unordered_map> nutrient_equations_; std::unordered_map> food_prices_; - std::map> + std::unordered_map> rural_prevalence_; std::vector income_models_; std::unordered_map> age_mean_height_; From 0ab8ef60bf676c2342be371c9584b83144f31812 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Thu, 26 Oct 2023 15:23:14 +0100 Subject: [PATCH 28/29] KH: use vector reserve, rather than letting it zero init. --- src/HealthGPS/kevin_hall_model.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index ba4e496f3..ef868ac5f 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -176,7 +176,8 @@ void KevinHallModel::update_sector(RuntimeContext &context, Person &person) cons void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) const { // Compute logits for each income category. - auto logits = std::vector(income_models_.size()); + auto logits = std::vector{}; + logits.reserve(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) { @@ -185,7 +186,8 @@ void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) } // Compute softmax probabilities for each income category. - auto e_logits = std::vector(income_models_.size()); + auto e_logits = std::vector{}; + e_logits.reserve(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]); @@ -193,7 +195,8 @@ void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) } // Compute income category probabilities. - auto probabilities = std::vector(income_models_.size()); + auto probabilities = std::vector{}; + probabilities.reserve(income_models_.size()); for (size_t i = 0; i < income_models_.size(); i++) { probabilities[i] = e_logits[i] / e_logits_sum; } From 1ebaa1acbb4a8707985e2486b39ad40ca8f1cc36 Mon Sep 17 00:00:00 2001 From: James Paul Turner Date: Thu, 26 Oct 2023 16:23:35 +0100 Subject: [PATCH 29/29] KH: push_back, don't assign, into reserved vector. --- src/HealthGPS/kevin_hall_model.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/HealthGPS/kevin_hall_model.cpp b/src/HealthGPS/kevin_hall_model.cpp index ef868ac5f..8635ce215 100644 --- a/src/HealthGPS/kevin_hall_model.cpp +++ b/src/HealthGPS/kevin_hall_model.cpp @@ -178,10 +178,10 @@ void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) // Compute logits for each income category. auto logits = std::vector{}; logits.reserve(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); + for (const auto &income_model : income_models_) { + logits.push_back(income_model.intercept); + for (const auto &[factor_name, coefficient] : income_model.coefficients) { + logits.back() += coefficient * person.get_risk_factor_value(factor_name); } } @@ -189,21 +189,21 @@ void KevinHallModel::initialise_income(RuntimeContext &context, Person &person) auto e_logits = std::vector{}; e_logits.reserve(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]; + for (const auto &logit : logits) { + e_logits.push_back(exp(logit)); + e_logits_sum += e_logits.back(); } // Compute income category probabilities. auto probabilities = std::vector{}; probabilities.reserve(income_models_.size()); - for (size_t i = 0; i < income_models_.size(); i++) { - probabilities[i] = e_logits[i] / e_logits_sum; + for (const auto &e_logit : e_logits) { + probabilities.push_back(e_logit / e_logits_sum); } // Compute income category. double rand = context.random().next_double(); - for (size_t i = 0; i < probabilities.size(); i++) { + for (size_t i = 0; i < income_models_.size(); i++) { if (rand < probabilities[i]) { person.income = income_models_[i].name; return;