From 337e5314a3c7342458c9b8c62aae0bf2a9535a72 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:04:53 -0700 Subject: [PATCH 01/10] refactor: remove obsolete compatibility paths --- tools/sync-gsa-results | 57 ------------------- workflow/rules/common.smk | 7 +-- .../prepare_faostat_animal_production.py | 31 +++------- 3 files changed, 9 insertions(+), 86 deletions(-) delete mode 100755 tools/sync-gsa-results diff --git a/tools/sync-gsa-results b/tools/sync-gsa-results deleted file mode 100755 index 23a8f476..00000000 --- a/tools/sync-gsa-results +++ /dev/null @@ -1,57 +0,0 @@ -#!/usr/bin/env bash -# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek -# -# SPDX-License-Identifier: GPL-3.0-or-later - -# Temporary script: rsync analysis results from shr and rename pce_* → gsa_* -# Usage: tools/sync-gsa-results -set -euo pipefail - -REMOTE="shr:~/GLADE-scratch" -REMOTE_NAME="pce_sensitivity" -LOCAL_NAME="gsa" - -LOCAL_ANALYSIS="results/${LOCAL_NAME}/analysis" -mkdir -p "$LOCAL_ANALYSIS" - -# Prefixes to rename: old → new -declare -A RENAMES=( - ["pce_"]="gsa_" - ["pce-l1-0p05_"]="gsa-l1-0p05_" - ["pce-l1-0p5_"]="gsa-l1-0p5_" -) - -# Rsync each prefix group, using --link-dest trick isn't needed since we -# rename on the fly. Instead: rsync to a temp dir, then rename in bulk. -TMPDIR=$(mktemp -d "${LOCAL_ANALYSIS}/.sync-XXXXXX") -trap 'rm -rf "$TMPDIR"' EXIT - -echo "==> Syncing analysis directories from ${REMOTE}..." -rsync -a --whole-file --compress --info=progress2 \ - --include='scen-pce*/' --include='scen-pce*/**' --exclude='*' \ - "${REMOTE}/results/${REMOTE_NAME}/analysis/" \ - "$TMPDIR/" - -echo "==> Renaming pce_* → gsa_*..." -renamed=0 -for old_prefix in "${!RENAMES[@]}"; do - new_prefix="${RENAMES[$old_prefix]}" - for dir in "$TMPDIR"/scen-${old_prefix}*; do - [ -d "$dir" ] || continue - base=$(basename "$dir") - newbase="${base/scen-${old_prefix}/scen-${new_prefix}}" - target="${LOCAL_ANALYSIS}/${newbase}" - if [ -d "$target" ]; then - # Merge into existing (overwrite) - cp -a "$dir"/* "$target"/ - else - mv "$dir" "$target" - fi - ((renamed++)) - done -done - -rm -rf "$TMPDIR" -trap - EXIT - -echo "==> Done. Renamed ${renamed} scenario directories into ${LOCAL_ANALYSIS}/" diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index aae65773..5ea2af32 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -70,15 +70,12 @@ def list_scenarios(): def get_effective_config(scenario_name): """Return the configuration with scenario overrides applied.""" - scenario_defs = load_scenario_defs() - # Start with a deep copy of the global config to avoid mutating it # We convert config to dict because it might be a Config object eff_config = copy.deepcopy(dict(config)) - if scenario_name and scenario_name in scenario_defs: - overrides = scenario_defs[scenario_name] - _recursive_update(eff_config, overrides) + if scenario_name: + _recursive_update(eff_config, load_scenario_defs()[scenario_name]) return eff_config diff --git a/workflow/scripts/prepare_faostat_animal_production.py b/workflow/scripts/prepare_faostat_animal_production.py index 73ea090d..8cffab5e 100644 --- a/workflow/scripts/prepare_faostat_animal_production.py +++ b/workflow/scripts/prepare_faostat_animal_production.py @@ -132,9 +132,7 @@ def main() -> None: "No FAOSTAT production records matched the configured animal products" ) - # Egg unit handling. Only hen eggs are mapped to the "eggs" model - # product (config/default.yaml). FAOSTAT may report egg production in - # tonnes ("t") or in thousands of eggs ("1000 No"). + # FAOSTAT QCL reports current egg production in tonnes. egg_mask_raw = df["product"] == "eggs" if egg_mask_raw.any(): if "Unit" not in df.columns: @@ -144,29 +142,14 @@ def main() -> None: egg_units = sorted(egg_unit_series.loc[egg_mask_raw].unique()) logger.info("Egg production units in source data: %s", ", ".join(egg_units)) - # Exact-match unit detection; avoid substring matches that could - # accept unexpected vintages silently. - thousand_eggs_mask = egg_mask_raw & egg_unit_series.isin({"1000 no"}) - tonnes_mask = egg_mask_raw & egg_unit_series.isin({"t"}) - unknown_mask = egg_mask_raw & ~(thousand_eggs_mask | tonnes_mask) - - if unknown_mask.any(): - unknown_units = sorted(egg_unit_series.loc[unknown_mask].unique()) + unexpected_units = sorted( + egg_unit_series.loc[egg_mask_raw & (egg_unit_series != "t")].unique() + ) + if unexpected_units: raise RuntimeError( "Unexpected FAOSTAT egg unit(s): " - + ", ".join(unknown_units) - + ". Expected 't' or '1000 No'." - ) - - if thousand_eggs_mask.any(): - # egg_thousands * 1000 eggs/thousand * 60g/egg / 1000 g/kg / 1000 kg/tonne - # = egg_thousands * 0.06 tonnes - df.loc[thousand_eggs_mask, "Value"] = ( - df.loc[thousand_eggs_mask, "Value"] * 0.06 - ) - logger.info( - "Converted %d egg records from '1000 No' to tonnes (60 g/egg)", - int(thousand_eggs_mask.sum()), + + ", ".join(unexpected_units) + + ". Expected 't'." ) result = ( From 4ccaeb8b001ce42a7908c4ba1b47f74166f0a37b Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:08:33 -0700 Subject: [PATCH 02/10] refactor: consolidate shared configuration helpers --- tests/test_ihme.py | 17 +++++++ tests/test_solve_namespace.py | 16 ++++++ workflow/rules/common.smk | 20 +------- workflow/scripts/ihme.py | 45 +++++++++++++++++ .../scripts/prepare_gbd_food_group_intake.py | 46 +---------------- workflow/scripts/prepare_gbd_mortality.py | 50 +------------------ workflow/scripts/solve_namespace.py | 3 +- 7 files changed, 85 insertions(+), 112 deletions(-) create mode 100644 tests/test_ihme.py create mode 100644 workflow/scripts/ihme.py diff --git a/tests/test_ihme.py b/tests/test_ihme.py new file mode 100644 index 00000000..73fb931a --- /dev/null +++ b/tests/test_ihme.py @@ -0,0 +1,17 @@ +# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from workflow.scripts.ihme import country_name_to_iso3 + + +def test_country_name_to_iso3_prefers_ihme_override() -> None: + assert country_name_to_iso3("Niger") == "NER" + + +def test_country_name_to_iso3_uses_pycountry() -> None: + assert country_name_to_iso3("Canada") == "CAN" + + +def test_country_name_to_iso3_returns_none_for_unknown_name() -> None: + assert country_name_to_iso3("Not a country") is None diff --git a/tests/test_solve_namespace.py b/tests/test_solve_namespace.py index 29d8f706..9be03129 100644 --- a/tests/test_solve_namespace.py +++ b/tests/test_solve_namespace.py @@ -9,6 +9,7 @@ import yaml from workflow.scripts.solve_namespace import ( + get_effective_config, validate_scenario_config_schemas, validate_scenario_overrides, ) @@ -84,3 +85,18 @@ def test_validates_one_representative_per_structure(self, base_config, monkeypat defs = {f"gsa_{i}": {"emissions": {"ghg_price": float(i)}} for i in range(50)} validate_scenario_config_schemas(base_config, defs, ".") assert len(calls) == 1 + + +def test_get_effective_config_applies_nested_overrides_without_mutating_base() -> None: + base = {"solving": {"solver": "highs", "threads": 1}} + scenarios = {"parallel": {"solving": {"threads": 4}}} + + effective = get_effective_config(base, "parallel", scenarios) + + assert effective == {"solving": {"solver": "highs", "threads": 4}} + assert base["solving"]["threads"] == 1 + + +def test_get_effective_config_rejects_unknown_scenario() -> None: + with pytest.raises(KeyError, match="missing"): + get_effective_config({}, "missing", {}) diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 5ea2af32..aa664dbc 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -7,7 +7,6 @@ Common configuration variables and helper functions shared across Snakemake rule This file should be included first in the main Snakefile, before any other rule files. """ -import copy import csv import hashlib import json @@ -22,6 +21,7 @@ from workflow.scripts.solve_namespace import ( _is_solve_time_key, _leaf_keys, deviation_penalty_uses_calibrated, + get_effective_config as _get_effective_config, health_input_paths, resolve_gbd_anchoring, resolve_pathvars, @@ -45,15 +45,6 @@ if ( ) -def _recursive_update(target, source): - for key, value in source.items(): - if isinstance(value, dict) and key in target and isinstance(target[key], dict): - _recursive_update(target[key], value) - else: - target[key] = value - return target - - def load_scenario_defs(): """Load scenario definitions from the config's `scenarios` key.""" global _SCENARIO_CACHE @@ -70,14 +61,7 @@ def list_scenarios(): def get_effective_config(scenario_name): """Return the configuration with scenario overrides applied.""" - # Start with a deep copy of the global config to avoid mutating it - # We convert config to dict because it might be a Config object - eff_config = copy.deepcopy(dict(config)) - - if scenario_name: - _recursive_update(eff_config, load_scenario_defs()[scenario_name]) - - return eff_config + return _get_effective_config(dict(config), scenario_name, load_scenario_defs()) def health_required(): diff --git a/workflow/scripts/ihme.py b/workflow/scripts/ihme.py new file mode 100644 index 00000000..fffdc286 --- /dev/null +++ b/workflow/scripts/ihme.py @@ -0,0 +1,45 @@ +# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek +# +# SPDX-License-Identifier: GPL-3.0-or-later + +"""Shared helpers for IHME source data.""" + +import pycountry + +COUNTRY_NAME_OVERRIDES = { + "Bolivia (Plurinational State of)": "BOL", + "Bonaire, Saint Eustatius and Saba": "BES", + "Cabo Verde": "CPV", + "Côte d'Ivoire": "CIV", + "Democratic People's Republic of Korea": "PRK", + "Democratic Republic of the Congo": "COD", + "French Guiana": "GUF", + "Iran (Islamic Republic of)": "IRN", + "Lao People's Democratic Republic": "LAO", + "Micronesia (Federated States of)": "FSM", + "Niger": "NER", + "Republic of Korea": "KOR", + "Republic of Moldova": "MDA", + "Republic of the Congo": "COG", + "Saint Barthélemy": "BLM", + "Saint Martin (French part)": "MAF", + "Sint Maarten (Dutch part)": "SXM", + "The former Yugoslav Republic of Macedonia": "MKD", + "Türkiye": "TUR", + "United Kingdom of Great Britain and Northern Ireland": "GBR", + "United Republic of Tanzania": "TZA", + "United States of America": "USA", + "United States Virgin Islands": "VIR", + "Venezuela (Bolivarian Republic of)": "VEN", + "Viet Nam": "VNM", +} + + +def country_name_to_iso3(name: str) -> str | None: + """Map an IHME location name to an ISO3 code.""" + if name in COUNTRY_NAME_OVERRIDES: + return COUNTRY_NAME_OVERRIDES[name] + try: + return pycountry.countries.search_fuzzy(name)[0].alpha_3 + except LookupError: + return None diff --git a/workflow/scripts/prepare_gbd_food_group_intake.py b/workflow/scripts/prepare_gbd_food_group_intake.py index 02102e71..d54920e6 100644 --- a/workflow/scripts/prepare_gbd_food_group_intake.py +++ b/workflow/scripts/prepare_gbd_food_group_intake.py @@ -55,41 +55,12 @@ from pathlib import Path import pandas as pd -import pycountry +from workflow.scripts.ihme import country_name_to_iso3 from workflow.scripts.logging_config import setup_script_logging logger = logging.getLogger(__name__) -# Reuse country name overrides from prepare_gbd_mortality.py -COUNTRY_NAME_OVERRIDES = { - "Bolivia (Plurinational State of)": "BOL", - "Bonaire, Saint Eustatius and Saba": "BES", - "Cabo Verde": "CPV", - "Côte d'Ivoire": "CIV", - "Democratic People's Republic of Korea": "PRK", - "Democratic Republic of the Congo": "COD", - "French Guiana": "GUF", - "Iran (Islamic Republic of)": "IRN", - "Lao People's Democratic Republic": "LAO", - "Micronesia (Federated States of)": "FSM", - "Niger": "NER", - "Republic of Korea": "KOR", - "Republic of Moldova": "MDA", - "Republic of the Congo": "COG", - "Saint Barthélemy": "BLM", - "Saint Martin (French part)": "MAF", - "Sint Maarten (Dutch part)": "SXM", - "The former Yugoslav Republic of Macedonia": "MKD", - "Türkiye": "TUR", - "United Kingdom of Great Britain and Northern Ireland": "GBR", - "United Republic of Tanzania": "TZA", - "United States of America": "USA", - "United States Virgin Islands": "VIR", - "Venezuela (Bolivarian Republic of)": "VEN", - "Viet Nam": "VNM", -} - # Map GBD 2023 risk factor file tokens to model food group names. # The token is the part of the filename between "..._DIET_" and the # "_Y2025..." date suffix (e.g. "LOW_IN_FRUITS", "HIGH_IN_RED_MEAT"). @@ -137,19 +108,6 @@ } -def map_country_name_to_iso3(name: str) -> str | None: - """Map GBD location name to ISO3 code using pycountry + manual overrides.""" - if name in COUNTRY_NAME_OVERRIDES: - return COUNTRY_NAME_OVERRIDES[name] - try: - matches = pycountry.countries.search_fuzzy(name) - if matches: - return matches[0].alpha_3 - except LookupError: - pass - return None - - def _wmean(values: pd.Series, weights: pd.Series) -> float: """Population-weighted mean, ignoring zero/NaN total weight.""" total = float(weights.sum()) @@ -186,7 +144,7 @@ def build_national_location_map(death_rates_path: str) -> dict[int, str]: ).drop_duplicates() loc_to_iso3: dict[int, str] = {} for location_id, location_name in ref.itertuples(index=False): - iso3 = map_country_name_to_iso3(location_name) + iso3 = country_name_to_iso3(location_name) if iso3 is not None: loc_to_iso3[int(location_id)] = iso3 logger.info( diff --git a/workflow/scripts/prepare_gbd_mortality.py b/workflow/scripts/prepare_gbd_mortality.py index 3133b460..1a7db115 100644 --- a/workflow/scripts/prepare_gbd_mortality.py +++ b/workflow/scripts/prepare_gbd_mortality.py @@ -17,42 +17,13 @@ from pathlib import Path import pandas as pd -import pycountry +from workflow.scripts.ihme import country_name_to_iso3 from workflow.scripts.logging_config import setup_script_logging # Logger will be configured in __main__ block logger = logging.getLogger(__name__) -# Manual overrides for country names that pycountry can't match -COUNTRY_NAME_OVERRIDES = { - "Bolivia (Plurinational State of)": "BOL", - "Bonaire, Saint Eustatius and Saba": "BES", - "Cabo Verde": "CPV", - "Côte d'Ivoire": "CIV", - "Democratic People's Republic of Korea": "PRK", - "Democratic Republic of the Congo": "COD", - "French Guiana": "GUF", # Use French data for French Guiana - "Iran (Islamic Republic of)": "IRN", - "Lao People's Democratic Republic": "LAO", - "Micronesia (Federated States of)": "FSM", - "Niger": "NER", # pycountry fuzzy search confuses with Nigeria (NGA) - "Republic of Korea": "KOR", - "Republic of Moldova": "MDA", - "Republic of the Congo": "COG", - "Saint Barthélemy": "BLM", - "Saint Martin (French part)": "MAF", - "Sint Maarten (Dutch part)": "SXM", - "The former Yugoslav Republic of Macedonia": "MKD", - "Türkiye": "TUR", - "United Kingdom of Great Britain and Northern Ireland": "GBR", - "United Republic of Tanzania": "TZA", - "United States of America": "USA", - "United States Virgin Islands": "VIR", - "Venezuela (Bolivarian Republic of)": "VEN", - "Viet Nam": "VNM", -} - # Map IHME cause names to model cause codes. The model's "Stroke" cause is # restricted to ischemic stroke (see prepare_relative_risks.CAUSE_MAP). CAUSE_MAP = { @@ -93,23 +64,6 @@ } -def map_country_name_to_iso3(name: str) -> str | None: - """Map IHME country name to ISO3 code using pycountry + manual overrides.""" - # Check manual overrides first - if name in COUNTRY_NAME_OVERRIDES: - return COUNTRY_NAME_OVERRIDES[name] - - # Try pycountry fuzzy search - try: - matches = pycountry.countries.search_fuzzy(name) - if matches: - return matches[0].alpha_3 - except LookupError: - pass - - return None - - def main() -> None: snakemake = globals().get("snakemake") # type: ignore if snakemake is None: @@ -158,7 +112,7 @@ def main() -> None: # Map country names to ISO3 (cache unique values to avoid repeated fuzzy searches) logger.info("Mapping country names to ISO3 codes...") unique_countries = df["location_name"].unique() - country_map = {name: map_country_name_to_iso3(name) for name in unique_countries} + country_map = {name: country_name_to_iso3(name) for name in unique_countries} df["country_iso3"] = df["location_name"].map(country_map) unmapped = df[df["country_iso3"].isna()]["location_name"].unique() diff --git a/workflow/scripts/solve_namespace.py b/workflow/scripts/solve_namespace.py index e6979d93..e48f0b8f 100644 --- a/workflow/scripts/solve_namespace.py +++ b/workflow/scripts/solve_namespace.py @@ -293,8 +293,7 @@ def get_effective_config( eff = copy.deepcopy(base_config) if not scenario_name: return eff - if scenario_name in scenario_defs: - _recursive_update(eff, scenario_defs[scenario_name]) + _recursive_update(eff, scenario_defs[scenario_name]) return eff From 47146e69aee9de269ed050f3a8ed99cb3a59b215 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:10:16 -0700 Subject: [PATCH 03/10] refactor: remove dead workflow and tooling paths --- tools/analyze_solve_profile.py | 141 --------- tools/export_model.py | 296 ------------------ tools/migrate-csv-to-parquet | 61 ---- tools/profile_solve_default.py | 82 ----- workflow/rules/analysis.smk | 13 +- workflow/rules/model.smk | 31 +- .../scripts/analysis/extract_water_metrics.py | 12 - workflow/scripts/build_model.py | 39 +-- workflow/scripts/build_model/biomass.py | 4 - workflow/scripts/solve_namespace.py | 6 + 10 files changed, 27 insertions(+), 658 deletions(-) delete mode 100644 tools/analyze_solve_profile.py delete mode 100755 tools/export_model.py delete mode 100755 tools/migrate-csv-to-parquet delete mode 100644 tools/profile_solve_default.py diff --git a/tools/analyze_solve_profile.py b/tools/analyze_solve_profile.py deleted file mode 100644 index 2a5479c5..00000000 --- a/tools/analyze_solve_profile.py +++ /dev/null @@ -1,141 +0,0 @@ -#!/usr/bin/env python -# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek -# -# SPDX-License-Identifier: GPL-3.0-or-later - -"""Analyze a cProfile dump from solve_model. - -Run:: - - pixi run python tools/analyze_solve_profile.py path/to/model_scen-X.prof - -Reports: - * Top functions by cumulative and total (self) time - * Aggregate time spent in linopy.* and pypsa.* - * Cumulative time of high-level checkpoints - (create_model, add_*_constraints, solve, etc.) -""" - -import argparse -from collections import defaultdict -from pathlib import Path -import pstats -import sys - - -def aggregate_by_package(stats: pstats.Stats) -> dict[str, tuple[float, float, int]]: - """Sum tottime and cumtime by package prefix.""" - buckets: dict[str, list[float]] = defaultdict(lambda: [0.0, 0.0, 0]) - for func, (_cc, _nc, tt, ct, _callers) in stats.stats.items(): - filename, _line, _name = func - if filename in ("~", ""): - bucket = "" - else: - parts = Path(filename).parts - # Map to top-level package directory under site-packages - try: - sp_idx = parts.index("site-packages") - bucket = parts[sp_idx + 1] if sp_idx + 1 < len(parts) else "" - except ValueError: - if "linopy" in filename: - bucket = "linopy" - elif "pypsa" in filename: - bucket = "pypsa" - elif "xarray" in filename: - bucket = "xarray" - elif "/workflow/" in filename or filename.startswith("workflow"): - bucket = "workflow" - else: - bucket = "" - b = buckets[bucket] - b[0] += tt - b[1] += ct # not strictly additive, but useful as upper bound - b[2] += 1 - return {k: (v[0], v[1], v[2]) for k, v in buckets.items()} - - -def find_checkpoint_times(stats: pstats.Stats, names: list[str]) -> dict[str, float]: - """Return cumtime for each function name (matches first occurrence).""" - found: dict[str, float] = {} - for func, (_cc, _nc, _tt, ct, _callers) in stats.stats.items(): - _filename, _line, fname = func - if fname in names and fname not in found: - found[fname] = ct - return found - - -def main() -> None: - ap = argparse.ArgumentParser() - ap.add_argument("profile", type=Path) - ap.add_argument("--top", type=int, default=30) - args = ap.parse_args() - - if not args.profile.is_file(): - print(f"profile not found: {args.profile}", file=sys.stderr) - sys.exit(1) - - stats = pstats.Stats(str(args.profile)) - stats.strip_dirs() - - total_runtime = stats.total_tt - print(f"Profile: {args.profile}") - print(f"Total cprofile runtime (tottime sum): {total_runtime:.2f} s\n") - - print("=" * 80) - print(f"Top {args.top} by cumulative time") - print("=" * 80) - stats.sort_stats("cumulative") - stats.print_stats(args.top) - - print("=" * 80) - print(f"Top {args.top} by self (tot) time") - print("=" * 80) - stats.sort_stats("tottime") - stats.print_stats(args.top) - - print("=" * 80) - print("Aggregate self time by package") - print("=" * 80) - by_pkg = aggregate_by_package(stats) - rows = sorted(by_pkg.items(), key=lambda kv: -kv[1][0]) - print(f"{'package':<30} {'tottime [s]':>12} {'n_funcs':>10}") - for pkg, (tt, _ct, n) in rows[:20]: - print(f"{pkg:<30} {tt:>12.3f} {n:>10}") - - print("\n" + "=" * 80) - print("Checkpoint cumulative times (high-level phases)") - print("=" * 80) - checkpoints = [ - "run_solve", - "_run_solve", - "create_model", - "_run", # linopy solver call wrapper - "solve", - "to_gurobipy", - "add_constraints", - "_extract_p_set_duals", - "assign_solution", - "assign_duals", - "post_processing", - "add_health_objective", - "add_production_stability_constraints", - "add_diet_stability_constraints", - "add_residue_feed_constraints", - "add_within_group_ratio_constraints", - "add_macronutrient_constraints", - "add_food_group_constraints", - "add_animal_growth_cap_constraints", - "add_crop_growth_cap_constraints", - "add_bounded_subsidy_constraints", - "_apply_forage_calibration", - ] - found = find_checkpoint_times(stats, checkpoints) - for name in checkpoints: - if name in found: - print(f" {name:<48} {found[name]:>8.2f} s") - else: - print(f" {name:<48} ") - - -if __name__ == "__main__": - main() diff --git a/tools/export_model.py b/tools/export_model.py deleted file mode 100755 index 3637f268..00000000 --- a/tools/export_model.py +++ /dev/null @@ -1,296 +0,0 @@ -#!/usr/bin/env python3 -# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek -# -# SPDX-License-Identifier: GPL-3.0-or-later - -"""Export GLADE model to MPS format using Gurobi's native export. - -This script builds the complete optimization model (including all constraints -added at solve time) and exports it to MPS format using Gurobi's native API, -which properly preserves integer/binary variable types. - -Usage: - pixi run -e gurobi python tools/export_model.py \ - --config config/yll.yaml \ - --scenario yll_20000 - -The exported MPS file can then be used for Gurobi parameter tuning: - pixi run -e gurobi python tools/tune_model.py results/yll/exported/model_scen-yll_20000.mps -""" - -import argparse -import logging -import os -from pathlib import Path - -import pandas as pd -import pypsa -import yaml - -from workflow.scripts.snakemake_utils import _recursive_update, apply_scenario_config -from workflow.scripts.solve_model.core import ( - add_food_group_constraints, - add_food_incentives_to_objective, - add_ghg_pricing_to_objective, - add_macronutrient_constraints, - add_residue_feed_constraints, - build_residue_feed_fraction_by_country, -) -from workflow.scripts.solve_model.health import add_health_objective -from workflow.scripts.solve_namespace import resolve_calibration_source_paths - -# Enable new PyPSA components API -pypsa.options.api.new_components_api = True - -logging.basicConfig( - level=logging.INFO, - format="%(asctime)s - %(levelname)s - %(message)s", -) -logger = logging.getLogger(__name__) - - -def resolve_path_root(raw_path: str, key: str) -> Path: - """Resolve environment variables and user-home markers in a path root.""" - resolved = os.path.expanduser(os.path.expandvars(raw_path)) - if "$" in resolved: - raise ValueError(f"Unresolved environment variable in config.paths.{key}") - return Path(resolved) - - -def load_config(config_path: str) -> dict: - """Load a config file and merge with default.yaml.""" - project_root = Path(__file__).parent.parent - - # Load default config - default_path = project_root / "config" / "default.yaml" - with open(default_path, encoding="utf-8") as f: - config = yaml.safe_load(f) - - # Load and merge user config - with open(config_path, encoding="utf-8") as f: - user_config = yaml.safe_load(f) - - _recursive_update(config, user_config) - return resolve_calibration_source_paths(config) - - -def build_and_export_model( - config_path: str, - scenario: str, - output_path: Path | None = None, -) -> Path: - """Build the complete model and export to MPS using Gurobi.""" - # Load and apply scenario config - config = load_config(config_path) - apply_scenario_config(config, scenario) - - config_name = config["name"] - paths_cfg = config["paths"] - results_dir = ( - resolve_path_root(paths_cfg["results_root"], "results_root") / config_name - ) - - # Determine output path - if output_path is None: - output_dir = results_dir / "exported" - output_dir.mkdir(parents=True, exist_ok=True) - output_path = output_dir / f"model_scen-{scenario}.mps" - - # Load the built network - network_path = results_dir / f"build/model_scen-{scenario}.nc" - if not network_path.exists(): - raise FileNotFoundError( - f"Built network not found: {network_path}\n" - f"First run: tools/smk -e gurobi -j4 --configfile {config_path} " - f"-- {network_path}" - ) - - logger.info("Loading network from %s", network_path) - n = pypsa.Network(network_path) - - # Add GHG pricing if enabled - if config["emissions"]["ghg_pricing_enabled"]: - ghg_price = float(config["emissions"]["ghg_price"]) - add_ghg_pricing_to_objective(n, ghg_price) - logger.info("Added GHG pricing: $%.2f/tCO2", ghg_price) - - # Add food incentives if enabled - if config["food_incentives"]["enabled"]: - sources = config["food_incentives"]["sources"] - incentive_paths = [str(Path(s.format(name=config_name))) for s in sources] - existing_paths = [p for p in incentive_paths if Path(p).exists()] - if existing_paths: - add_food_incentives_to_objective(n, existing_paths) - logger.info("Added food incentives from %d sources", len(existing_paths)) - - # Create linopy model - logger.info("Creating linopy model...") - n.optimize.create_model() - - # Load population data - processing_dir = ( - resolve_path_root(paths_cfg["processing_root"], "processing_root") / config_name - ) - population_path = processing_dir / "population.csv" - if not population_path.exists(): - raise FileNotFoundError(f"Population data not found: {population_path}") - population_df = pd.read_csv(population_path) - population_df["iso3"] = population_df["iso3"].astype(str).str.upper() - population_map = ( - population_df.set_index("iso3")["population"].astype(float).to_dict() - ) - - # Add macronutrient constraints - macronutrients = config.get("macronutrients") - if macronutrients: - add_macronutrient_constraints(n, macronutrients, population_map) - logger.info("Added macronutrient constraints") - - # Add food group constraints - food_group_cfg = config.get("food_groups", {}).get("constraints") - if food_group_cfg: - add_food_group_constraints(n, food_group_cfg, population_map, None) - logger.info("Added food group constraints") - - # Add residue feed constraints - max_feed_fraction = float(config["residues"]["max_feed_fraction"]) - m49_path = Path("data/curated/M49-codes.csv") - if m49_path.exists(): - max_feed_by_country = build_residue_feed_fraction_by_country( - config, str(m49_path) - ) - add_residue_feed_constraints(n, max_feed_fraction, max_feed_by_country) - logger.info("Added residue feed constraints") - - # Add health objective if enabled - if config["health"]["enabled"]: - health_data_dir = processing_dir / "health" - clusters_path = processing_dir / "health/country_clusters.csv" - - required_health_files = [ - health_data_dir / "risk_breakpoints.csv", - health_data_dir / "cluster_cause_baseline.csv", - health_data_dir / "cause_log_breakpoints.csv", - health_data_dir / "cluster_summary.csv", - clusters_path, - ] - - if all(f.exists() for f in required_health_files): - risk_factors = config["health"]["risk_factors"] - risk_cause_map = config["health"]["risk_cause_map"] - value_per_yll = float(config["health"]["value_per_yll"]) - add_health_objective( - n, - str(health_data_dir / "risk_breakpoints.csv"), - str(health_data_dir / "cluster_cause_baseline.csv"), - str(health_data_dir / "cause_log_breakpoints.csv"), - str(health_data_dir / "cluster_summary.csv"), - str(clusters_path), - str(population_path), - risk_factors, - risk_cause_map, - "gurobi", - value_per_yll, - ) - logger.info("Added health objective (value_per_yll=%.2f)", value_per_yll) - else: - missing = [f for f in required_health_files if not f.exists()] - logger.warning("Health data incomplete, skipping. Missing: %s", missing[:3]) - - logger.info( - "Model built: %d variables, %d constraints", - n.model.nvars, - n.model.ncons, - ) - - # Use linopy's to_gurobipy to build the native Gurobi model with SOS constraints - logger.info("Building Gurobi model via linopy.Model.to_gurobipy()...") - - # Convert linopy model to gurobipy model (includes SOS constraints) - gp_model = n.model.to_gurobipy() - gp_model.update() - - # Check model characteristics - num_sos = gp_model.NumSOS - logger.info( - "Gurobi model: %d vars (%d binary, %d integer), %d constraints, %d SOS constraints", - gp_model.NumVars, - gp_model.NumBinVars, - gp_model.NumIntVars, - gp_model.NumConstrs, - num_sos, - ) - - if gp_model.NumBinVars == 0 and gp_model.NumIntVars == 0: - if num_sos > 0: - logger.info("Model type: LP with SOS constraints (not a MIP)") - logger.info( - "Note: Gurobi handles SOS2 natively. The model will solve as " - "a continuous problem with SOS branching." - ) - else: - logger.info("Model type: Pure LP (no integer variables or SOS constraints)") - else: - logger.info("Model type: MIP") - - # Export using Gurobi's native MPS writer - logger.info("Exporting to %s", output_path) - gp_model.write(str(output_path)) - - logger.info("Export complete: %s", output_path) - return output_path - - -def main(): - parser = argparse.ArgumentParser( - description="Export GLADE model to MPS format for Gurobi tuning", - formatter_class=argparse.RawDescriptionHelpFormatter, - epilog=""" -Example: - # First build the network: - tools/smk -e gurobi -j4 --configfile config/yll.yaml -- results/yll/build/model_scen-yll_20000.nc - - # Then export for tuning: - pixi run -e gurobi python tools/export_model.py --config config/yll.yaml --scenario yll_20000 - - # Run tuning: - pixi run -e gurobi python tools/tune_model.py results/yll/exported/model_scen-yll_20000.mps - """, - ) - parser.add_argument( - "--config", - "-c", - required=True, - help="Path to config file (e.g., config/yll.yaml)", - ) - parser.add_argument( - "--scenario", - "-s", - required=True, - help="Scenario name (e.g., yll_20000)", - ) - parser.add_argument( - "--output", - "-o", - type=Path, - default=None, - help=( - "Output MPS file path " - "(default: {paths.results_root}/{name}/exported/model_scen-{scenario}.mps)" - ), - ) - - args = parser.parse_args() - - output_path = build_and_export_model( - args.config, - args.scenario, - args.output, - ) - - logger.info("\nNext steps:") - logger.info(" pixi run -e gurobi python tools/tune_model.py %s", output_path) - - -if __name__ == "__main__": - main() diff --git a/tools/migrate-csv-to-parquet b/tools/migrate-csv-to-parquet deleted file mode 100755 index 256021da..00000000 --- a/tools/migrate-csv-to-parquet +++ /dev/null @@ -1,61 +0,0 @@ -#!/usr/bin/env python3 -# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek -# -# SPDX-License-Identifier: GPL-3.0-or-later -"""One-off migration: convert analysis CSVs to Parquet in-place. - -Usage: - pixi run python tools/migrate-csv-to-parquet results/pce_sensitivity/analysis/ - -Converts all .csv files under the given directory tree to .parquet, -then removes the originals. Skips files that already have a .parquet -companion. Prints progress every 1000 files. -""" - -from pathlib import Path -import sys - -import pandas as pd - - -def main(): - if len(sys.argv) != 2: - print(f"Usage: {sys.argv[0]} ", file=sys.stderr) - sys.exit(1) - - root = Path(sys.argv[1]) - if not root.is_dir(): - print(f"Not a directory: {root}", file=sys.stderr) - sys.exit(1) - - csvs = sorted(root.rglob("*.csv")) - print(f"Found {len(csvs)} CSV files under {root}") - - converted = 0 - skipped = 0 - errors = 0 - - for i, csv_path in enumerate(csvs, 1): - parquet_path = csv_path.with_suffix(".parquet") - - if parquet_path.exists(): - skipped += 1 - continue - - try: - df = pd.read_csv(csv_path) - df.to_parquet(parquet_path) - csv_path.unlink() - converted += 1 - except Exception as e: - print(f" ERROR {csv_path}: {e}", file=sys.stderr) - errors += 1 - - if i % 1000 == 0: - print(f" ... {i}/{len(csvs)} processed ({converted} converted)") - - print(f"Done: {converted} converted, {skipped} skipped, {errors} errors") - - -if __name__ == "__main__": - main() diff --git a/tools/profile_solve_default.py b/tools/profile_solve_default.py deleted file mode 100644 index 53b92f4a..00000000 --- a/tools/profile_solve_default.py +++ /dev/null @@ -1,82 +0,0 @@ -#!/usr/bin/env python3 -# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek -# -# SPDX-License-Identifier: GPL-3.0-or-later -"""Profile-run the opt default scenario directly, skipping IIS computation. - -This harness reads the manifest produced by ``tools/export-solve-manifest``, -calls ``run_solve`` under cProfile, and saves the profile next to the output. - -Usage: - pixi run -e gurobi python tools/profile_solve_default.py [manifest.json] -""" - -import argparse -import cProfile -import json -import logging -from pathlib import Path -import sys -import time - -PROJECT_ROOT = Path(__file__).resolve().parent.parent -sys.path.insert(0, str(PROJECT_ROOT)) - - -def main() -> None: - ap = argparse.ArgumentParser() - ap.add_argument( - "manifest", - nargs="?", - default="/tmp/opt_manifest.json", - type=Path, - ) - ap.add_argument( - "--scenario-index", type=int, default=0, help="Index into manifest scenarios" - ) - args = ap.parse_args() - - manifest = json.loads(args.manifest.read_text()) - entry = manifest["scenarios"][args.scenario_index] - shared = manifest.get("shared_params") - - out_net = entry["outputs"]["network"] - Path(out_net).parent.mkdir(parents=True, exist_ok=True) - Path(entry["log"]).parent.mkdir(parents=True, exist_ok=True) - - from workflow.scripts.logging_config import setup_script_logging - from workflow.scripts.solve_model.core import _ShadowPriceLogFilter, run_solve - from workflow.scripts.solve_namespace import build_namespace - - logger = setup_script_logging(entry["log"]) - logging.getLogger("pypsa.optimization.optimize").addFilter(_ShadowPriceLogFilter()) - - smk = build_namespace(entry, shared) - - # Stub IIS so an infeasible solve does not burn minutes on diagnostics. - import linopy.model - - linopy.model.Model.compute_infeasibilities = lambda self: [] # type: ignore[method-assign] - - print(f"Scenario: {entry['scenario']}") - print(f"Output: {out_net}") - - prof_path = Path(out_net).with_suffix(".prof") - profiler = cProfile.Profile() - t0 = time.perf_counter() - profiler.enable() - try: - n = run_solve(smk, logger) - finally: - profiler.disable() - profiler.dump_stats(str(prof_path)) - print(f"Wallclock: {time.perf_counter() - t0:.2f} s") - print(f"Profile: {prof_path}") - - if n is not None: - # Skip netcdf export; not needed for timing. - Path(out_net).touch() - - -if __name__ == "__main__": - main() diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 8a8cfa10..89355190 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -23,17 +23,8 @@ rule prepare_faostat_emissions: "../scripts/prepare_faostat_emissions.py" -_ANALYSIS_SCRIPTS = expand( - "workflow/scripts/analysis/{script}", - script=[ - "extract_statistics.py", - "extract_net_emissions.py", - "extract_objective_breakdown.py", - "extract_ghg_attribution.py", - "extract_health_impacts.py", - "extract_baseline_deviation.py", - "extract_food_prices.py", - ], +_ANALYSIS_SCRIPTS = sorted( + str(path) for path in Path("workflow/scripts/analysis").glob("extract_*.py") ) from workflow.scripts.solve_namespace import ANALYSIS_OUTPUT_NAMES diff --git a/workflow/rules/model.smk b/workflow/rules/model.smk index 4ef4494d..c3dfd2b0 100644 --- a/workflow/rules/model.smk +++ b/workflow/rules/model.smk @@ -65,15 +65,6 @@ def build_model_fiber_baseline_input(wildcards): return {} -def build_model_grassland_calibration_input(wildcards): - """Grassland forage calibration is now applied at solve time. - - This stub remains to avoid breaking unpack() calls in build_model inputs - while the transition settles. - """ - return {} - - def build_model_fodder_yield_correction_input(wildcards): """Conditionally include fodder yield correction CSV.""" if config["fodder_decomposition"]["yield_corrections"]["enabled"]: @@ -138,7 +129,6 @@ rule build_model: unpack(yield_inputs), unpack(residue_yield_inputs), unpack(harvested_area_model_inputs), - unpack(build_model_grassland_calibration_input), unpack(build_model_fodder_yield_correction_input), unpack(build_model_yield_calibration_input), unpack(build_model_cost_calibration_input), @@ -182,22 +172,8 @@ rule build_model: faostat_pasture_area="/{name}/faostat_pasture_area.csv", current_grassland_area="/{name}/luc/current_grassland_area_by_class.csv", grazing_only_land="/{name}/land_grazing_only_by_class.csv", - build_scripts=expand( - "workflow/scripts/build_model/{script}", - script=[ - "animals.py", - "biomass.py", - "health.py", - "crops.py", - "food.py", - "grassland.py", - "infrastructure.py", - "land.py", - "nutrition.py", - "primary_resources.py", - "trade.py", - "utils.py", - ], + build_scripts=sorted( + str(path) for path in Path("workflow/scripts/build_model").glob("*.py") ), constants_script="workflow/scripts/constants.py", params: @@ -250,6 +226,9 @@ def solve_model_inputs(w): "m49": "data/curated/M49-codes.csv", "food_groups": "data/curated/food_groups.csv", "baseline_diet": f"/{w.name}/baseline_diet.csv", + "solve_scripts": sorted( + str(path) for path in Path("workflow/scripts/solve_model").glob("*.py") + ), } eff_cfg = get_effective_config(w.scenario) diff --git a/workflow/scripts/analysis/extract_water_metrics.py b/workflow/scripts/analysis/extract_water_metrics.py index 061f04f4..530cf235 100644 --- a/workflow/scripts/analysis/extract_water_metrics.py +++ b/workflow/scripts/analysis/extract_water_metrics.py @@ -25,9 +25,6 @@ - non-renewable groundwater depletion (Mm3 mined) from ``groundwater_nonrenewable`` tiers. -It also exposes the per-tier (CF, draw, source) table used to bin withdrawn -water by the scarcity of its source (the merit-order / withdrawal-by-CF figure). - Basis: the tier draw is *consumption* C -- crops need their net requirement E and the ``irrigation_delivery`` link draws ``C = E / eta_c`` from the pool. ``withdrawn_mm3`` therefore reports consumption; the estimated physical @@ -75,15 +72,6 @@ def _water_supply_draw(n: pypsa.Network) -> pd.DataFrame: ) -def extract_water_tiers(n: pypsa.Network) -> pd.DataFrame: - """Per-tier (region, cf, draw_mm3) table; the input to CF-binned views. - - One row per ``water_supply`` link. ``cf`` is the tier's marginal - characterisation factor; ``draw_mm3`` is the volume drawn from it. - """ - return _water_supply_draw(n).reset_index(names="link") - - def extract_water_by_region(n: pypsa.Network) -> pd.DataFrame: """Per-region water metrics. diff --git a/workflow/scripts/build_model.py b/workflow/scripts/build_model.py index 3997c944..1150fdaa 100644 --- a/workflow/scripts/build_model.py +++ b/workflow/scripts/build_model.py @@ -61,7 +61,6 @@ def filter(self, record: logging.LogRecord) -> bool: validation_cfg = snakemake.params.validation # type: ignore[attr-defined] use_actual_production = bool(validation_cfg["use_actual_production"]) - enforce_baseline = bool(validation_cfg["enforce_baseline_diet"]) enforce_baseline_feed = bool(validation_cfg["enforce_baseline_feed"]) # Enable land slack if explicitly requested or when using actual production enable_land_slack = bool(validation_cfg["land_slack"]) or use_actual_production @@ -300,6 +299,10 @@ def filter(self, record: logging.LogRecord) -> bool: "where they can be grown." ) + # Read regions once for correction mapping and model construction. + regions_df = gpd.read_file(snakemake.input.regions) + region_to_country = regions_df.set_index("region")["country"] + # Apply per-(country, crop) yield corrections. Two sources, applied # multiplicatively through the same per-cell rescaling: # * fodder_yield_corrections: Eurostat-anchored FDD-crop corrections, @@ -311,15 +314,12 @@ def filter(self, record: logging.LogRecord) -> bool: def _apply_yield_corrections(corr_df: pd.DataFrame, source_label: str) -> None: if corr_df.empty: return - _r2c_corr = gpd.read_file(snakemake.input.regions)[ - ["region", "country"] - ].set_index("region")["country"] n_adjusted = 0 for _, row in corr_df.iterrows(): country = str(row["country"]) crop = str(row["crop"]) factor = float(row["yield_correction_factor"]) - corr_regions = set(_r2c_corr[_r2c_corr == country].index) + corr_regions = set(region_to_country[region_to_country == country].index) for ws in ("r", "i"): key = f"{crop}_yield_{ws}" if key not in yields_data: @@ -345,9 +345,6 @@ def _apply_yield_corrections(corr_df: pd.DataFrame, source_label: str) -> None: if yield_cal_path: _apply_yield_corrections(read_csv(yield_cal_path), "FAOSTAT-target") - # Read regions - regions_df = gpd.read_file(snakemake.input.regions) - # Load class-level land areas land_class_df = read_csv(snakemake.input.land_area_by_class) # Expect columns: region, water_supply, resource_class, area_ha @@ -380,23 +377,16 @@ def _apply_yield_corrections(corr_df: pd.DataFrame, source_label: str) -> None: ) ch4_to_co2_factor = float(snakemake.params.emissions["ch4_to_co2_factor"]) n2o_to_co2_factor = float(snakemake.params.emissions["n2o_to_co2_factor"]) - try: - luc_coefficients_path = snakemake.input.luc_carbon_coefficients - luc_coeff_df = read_csv(luc_coefficients_path) - if not luc_coeff_df.empty: - luc_lef_lookup = utils._build_luc_lef_lookup(luc_coeff_df) - logger.info( - "Loaded LUC LEFs for %d (region, class, water, use) combinations", - len(luc_lef_lookup), - ) - else: - logger.warning( - "LUC carbon coefficients file is empty; skipping LUC emission adjustments" - ) - except (AttributeError, FileNotFoundError) as e: + luc_coeff_df = read_csv(snakemake.input.luc_carbon_coefficients) + if not luc_coeff_df.empty: + luc_lef_lookup = utils._build_luc_lef_lookup(luc_coeff_df) logger.info( - "LUC carbon coefficients not available (%s); skipping LUC emission adjustments", - type(e).__name__, + "Loaded LUC LEFs for %d (region, class, water, use) combinations", + len(luc_lef_lookup), + ) + else: + logger.warning( + "LUC carbon coefficients file is empty; skipping LUC emission adjustments" ) land_rainfed_df = land_class_df.xs("r", level="water_supply").copy() @@ -727,7 +717,6 @@ def _apply_yield_corrections(corr_df: pd.DataFrame, source_label: str) -> None: food_to_group = food_groups.set_index("food")["group"].to_dict() food_group_list = list(snakemake.params.food_groups) - macronutrient_cfg = snakemake.params.macronutrients nutrient_units = ( nutrition_data.reset_index() .drop_duplicates(subset=["nutrient"]) diff --git a/workflow/scripts/build_model/biomass.py b/workflow/scripts/build_model/biomass.py index ed8940a9..9a49a0e7 100644 --- a/workflow/scripts/build_model/biomass.py +++ b/workflow/scripts/build_model/biomass.py @@ -227,10 +227,6 @@ def add_biofuel_links( bus_index = n.buses.static.index - # Ensure bus_type column exists (default "food" for backward compatibility) - if "bus_type" not in biofuel_baseline.columns: - biofuel_baseline = biofuel_baseline.assign(bus_type="food") - # Aggregate baseline demand by (source_item, crop, country, bus_type). # Link names below carry only (source_item, country), so the groupby # must produce a unique row per such pair; otherwise n.links.add diff --git a/workflow/scripts/solve_namespace.py b/workflow/scripts/solve_namespace.py index e48f0b8f..0cd64d30 100644 --- a/workflow/scripts/solve_namespace.py +++ b/workflow/scripts/solve_namespace.py @@ -366,6 +366,9 @@ def rp(path: str) -> str: "m49": "data/curated/M49-codes.csv", "food_groups": "data/curated/food_groups.csv", "baseline_diet": rp("/{name}/baseline_diet.csv"), + "solve_scripts": sorted( + str(path) for path in Path("workflow/scripts/solve_model").glob("*.py") + ), } # Health processing inputs only when this scenario enables health (mirrors @@ -424,6 +427,9 @@ def rp(path: str) -> str: if inline_analysis: inputs["population"] = rp("/{name}/population.csv") + inputs["analysis_scripts"] = sorted( + str(path) for path in Path("workflow/scripts/analysis").glob("extract_*.py") + ) params: dict = { "health_enabled": eff["health"]["enabled"], From f945aab4a5c022c8cd690dccf316ad02ccdf6ba1 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:14:37 -0700 Subject: [PATCH 04/10] refactor: remove unused configuration branches --- CHANGELOG.md | 6 ++ config/default.yaml | 5 -- config/gsa.yaml | 4 +- config/gsa_fixed_diet.yaml | 3 +- config/schemas/config.schema.yaml | 18 +--- .../calibration/default/provenance.yaml | 2 - .../calibration/gbd-anchored/provenance.yaml | 2 - docs/cluster_execution.rst | 6 +- docs/sensitivity_analysis.rst | 9 +- docs/surrogate_modelling.rst | 8 +- tests/config/test_sensitivity.yaml | 86 ------------------- tests/test_gbd2019_rr_appendix.py | 50 ++--------- tests/test_sensitivity.py | 54 ------------ workflow/scripts/gbd2019_rr_appendix.py | 69 +++------------ .../scripts/generate_rr_age_attenuation.py | 4 +- workflow/scripts/solve_model/sensitivity.py | 14 +-- 16 files changed, 42 insertions(+), 298 deletions(-) delete mode 100644 tests/config/test_sensitivity.yaml diff --git a/CHANGELOG.md b/CHANGELOG.md index 864cb1e4..2bb6bcc2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -249,6 +249,12 @@ introduce breaking changes to configuration and outputs. - The MARS surrogate method; supported surrogates are now `pce`, `rf`, `xgb` and `mlp`. +- Unused configuration keys `health.ssb_sugar_g_per_100g`, + `data.gaez.climate_model_ensemble`, and + `sensitivity_analysis.default_surrogate`. Surrogate methods are selected + explicitly in target and bundle names. Sensitivity scenarios must use the + separate `food_loss` and `food_waste` factors instead of the removed + `food_loss_waste` convenience key. ### Fixed diff --git a/config/default.yaml b/config/default.yaml index f34259d4..7c852597 100644 --- a/config/default.yaml +++ b/config/default.yaml @@ -1396,7 +1396,6 @@ health: region_clusters: 30 breakpoint_rel_tol: 0.05 # Max Stage 1 PWL deviation as a fraction of each RR curve's amplitude (5%) log_rr_points: 15 - ssb_sugar_g_per_100g: 5.7 # ~=50 kcal per 226.8 g sugar-sweetened beverage (SSB) implies ~5.7 g sugar per 100 g value_per_yll: 50000 # USD_2024 per year of life lost intake_age_min: 11 # GDD adult band starts at 11; set to 11 to retain adult intake data. Note however that GDB chronic disease risk factors are for adults of >=25 years. # Dietary risk factors to consider (must match GBD risk-factor items) @@ -1746,7 +1745,6 @@ data: # GAEZ v5 parameters # Note: RES05 (yields/suitability) has ENSEMBLE, but RES02 (growing season) only has individual GCMs climate_model: "GFDL-ESM4" # Specific GCMs: "GFDL-ESM4", "IPSL-CM6A-LR", "MPI-ESM1-2-HR", "MRI-ESM2-0", "UKESM1-0-LL" - climate_model_ensemble: "ENSEMBLE" # Multi-model mean (only available for RES05, not RES02) period: "FP2140" # Future: "FP2140" (2021-2040), "FP4160" (2041-2060), "FP6180" (2061-2080), "FP8100" (2081-2100); Historical: "HP0120" (2001-2020), "HP8100" (1981-2000) climate_scenario: "SSP126" # "SSP126" (low emissions), "SSP370" (medium, ~RCP4.5), "SSP585" (high), "HIST" (historical) input_level: "H" # "H" (High), "L" (Low) @@ -1948,9 +1946,6 @@ remote_solve: sensitivity_analysis: holdout_fraction: 0.15 threads: 6 - # Method downstream consumers (uncertainty-band plots, notebooks) use when - # no explicit choice is given. Must match a key under ``methods``. - default_surrogate: mlp # When false (default), the surrogate-fit rule declares every Sobol scenario # the generator promises so Snakemake drives the full solve→analyse→ # surrogate chain in a single invocation (the canonical Snakemake idiom). diff --git a/config/gsa.yaml b/config/gsa.yaml index 6f86b7f0..8461a0db 100644 --- a/config/gsa.yaml +++ b/config/gsa.yaml @@ -125,7 +125,8 @@ scenarios: ch4: "{ch4_factor}" n2o: "{n2o_factor}" luc: "{luc_factor}" - food_loss_waste: "{flw_factor}" + food_loss: "{flw_factor}" + food_waste: "{flw_factor}" feed_conversion: "{fcr_factor}" health_relative_risk: protective: "{rr_protective}" @@ -147,7 +148,6 @@ scenarios: sensitivity_analysis: holdout_fraction: 0.15 threads: 6 - default_surrogate: mlp # GSA solves run outside Snakemake on the cluster (tools/batch-solve), so the # surrogate fit reads scenarios from disk rather than declaring every output. discover_scenarios_on_disk: true diff --git a/config/gsa_fixed_diet.yaml b/config/gsa_fixed_diet.yaml index 7c6cbcee..ba4b91e7 100644 --- a/config/gsa_fixed_diet.yaml +++ b/config/gsa_fixed_diet.yaml @@ -99,7 +99,8 @@ scenarios: ch4: "{ch4_factor}" n2o: "{n2o_factor}" luc: "{luc_factor}" - food_loss_waste: "{flw_factor}" + food_loss: "{flw_factor}" + food_waste: "{flw_factor}" feed_conversion: "{fcr_factor}" land: reforestation_cap: diff --git a/config/schemas/config.schema.yaml b/config/schemas/config.schema.yaml index 7c8b7992..1a7c35f8 100644 --- a/config/schemas/config.schema.yaml +++ b/config/schemas/config.schema.yaml @@ -1318,7 +1318,6 @@ properties: - region_clusters - breakpoint_rel_tol - log_rr_points - - ssb_sugar_g_per_100g - value_per_yll - risk_factors - causes @@ -1355,10 +1354,6 @@ properties: type: integer minimum: 0 description: "Lower bound (inclusive) of age buckets used for intake weighting" - ssb_sugar_g_per_100g: - type: number - minimum: 0 - description: "g sugar per 100g sugar-sweetened beverage" value_per_yll: type: number minimum: 0 @@ -1504,7 +1499,6 @@ properties: type: object required: - climate_model - - climate_model_ensemble - period - climate_scenario - input_level @@ -1516,9 +1510,6 @@ properties: climate_model: type: string description: "GAEZ v5 climate model (specific GCM)" - climate_model_ensemble: - type: string - description: "GAEZ v5 ensemble model" period: type: string pattern: "^(FP|HP)[0-9]{4}$" @@ -2209,7 +2200,7 @@ properties: sensitivity_analysis: type: object description: "Configuration for surrogate-based sensitivity analysis methods" - required: [outputs, holdout_fraction, threads, default_surrogate, methods, sobol, discover_scenarios_on_disk] + required: [outputs, holdout_fraction, threads, methods, sobol, discover_scenarios_on_disk] additionalProperties: false properties: outputs: @@ -2250,9 +2241,6 @@ properties: type: integer minimum: 1 description: "Number of threads for sensitivity analysis (controls sklearn n_jobs and BLAS)" - default_surrogate: - type: string - description: "Surrogate method downstream consumers (e.g. uncertainty-band plots, notebooks) pick when no method is specified explicitly. Must match a key under ``methods``." discover_scenarios_on_disk: type: boolean description: "If false (default), the surrogate-fit rule declares every scenario the generator promises so Snakemake drives the full solve→analyse→surrogate chain. If true, the rule scans the analysis directory and fits the surrogate on whatever scenarios have complete outputs on disk; intended for cluster sweeps where solves run outside Snakemake. Errors out if more than 50% of scenarios are missing." @@ -2423,10 +2411,6 @@ properties: type: number minimum: 0 description: "Factor for CO2 from land-use change" - food_loss_waste: - type: number - minimum: 0 - description: "Bundle convenience key: applies the same multiplicative factor to both food_loss and food_waste. Per-component keys take precedence when both are set. factor=1.5 raises the underlying loss/waste fractions by 50%; factor=1.0 is a no-op." food_loss: type: number minimum: 0 diff --git a/data/curated/calibration/default/provenance.yaml b/data/curated/calibration/default/provenance.yaml index 7512458e..ddba1a95 100644 --- a/data/curated/calibration/default/provenance.yaml +++ b/data/curated/calibration/default/provenance.yaml @@ -639,7 +639,6 @@ structural_config: data.faostat.fbs_production_element_code: 5511 data.faostat.qcl_production_element_code: 5510 data.gaez.climate_model: GFDL-ESM4 - data.gaez.climate_model_ensemble: ENSEMBLE data.gaez.climate_scenario: SSP126 data.gaez.input_level: H data.gaez.period: FP2140 @@ -969,7 +968,6 @@ structural_config: - legumes - red_meat - whole_grains - health.ssb_sugar_g_per_100g: 5.7 irrigation.irrigated_crops: all land.conversion_cost_forest_usd_per_ha: 8000 land.conversion_cost_nonforest_usd_per_ha: 2000 diff --git a/data/curated/calibration/gbd-anchored/provenance.yaml b/data/curated/calibration/gbd-anchored/provenance.yaml index a0925454..fafd0e77 100644 --- a/data/curated/calibration/gbd-anchored/provenance.yaml +++ b/data/curated/calibration/gbd-anchored/provenance.yaml @@ -639,7 +639,6 @@ structural_config: data.faostat.fbs_production_element_code: 5511 data.faostat.qcl_production_element_code: 5510 data.gaez.climate_model: GFDL-ESM4 - data.gaez.climate_model_ensemble: ENSEMBLE data.gaez.climate_scenario: SSP126 data.gaez.input_level: H data.gaez.period: FP2140 @@ -969,7 +968,6 @@ structural_config: - legumes - red_meat - whole_grains - health.ssb_sugar_g_per_100g: 5.7 irrigation.irrigated_crops: all land.conversion_cost_forest_usd_per_ha: 8000 land.conversion_cost_nonforest_usd_per_ha: 2000 diff --git a/docs/cluster_execution.rst b/docs/cluster_execution.rst index 2c059ba7..aebc8a44 100644 --- a/docs/cluster_execution.rst +++ b/docs/cluster_execution.rst @@ -188,10 +188,8 @@ flag at its default ``false``. --allowed-rules build_surrogate \ -- results/gsa/surrogates/surrogate_gsa_xgb.pkl -The method (``xgb``/``pce``/``rf``/``mlp``) is a wildcard of the rule; -``sensitivity_analysis.default_surrogate`` in the config picks the -default when downstream consumers (uncertainty plots, notebooks) load -a bundle without specifying a method. +The method (``xgb``/``pce``/``rf``/``mlp``) is a wildcard of the rule and is +always explicit in downstream bundle and plot targets. Repeat the command for each ``{group}`` needed (e.g., ``gsa-l1-low``, ``gsa-l1-high``) and for alternative surrogate types if you want to diff --git a/docs/sensitivity_analysis.rst b/docs/sensitivity_analysis.rst index 72abe5a6..14bdeecf 100644 --- a/docs/sensitivity_analysis.rst +++ b/docs/sensitivity_analysis.rst @@ -234,7 +234,6 @@ distributions rather than fixed value lists. # Surrogate fitting + Sobol settings (see surrogate_modelling for methods) sensitivity_analysis: holdout_fraction: 0.15 - default_surrogate: mlp sobol: outputs: [total_cost, co2, ch4, n2o, land_use, yll] grid_resolution: 15 # conditional-Sobol grid points per slice axis @@ -287,8 +286,8 @@ parameter per risk factor produces cause-specific adjustments automatically. is the entire value. **``sensitivity_analysis.sobol`` field reference** (the Sobol-index settings; -the surrogate-fitting fields -- ``methods``, ``outputs``, -``default_surrogate``, ``holdout_fraction``, ``discover_scenarios_on_disk`` -- +the surrogate-fitting fields -- ``methods``, ``outputs``, ``holdout_fraction``, +``discover_scenarios_on_disk`` -- are documented in :doc:`surrogate_modelling`): - ``outputs``: Allowlist of output names whose Sobol indices are computed and @@ -817,9 +816,7 @@ Sobol computation, policy sweeps, and uncertainty-band plots consume. Output paths use two wildcards: ``{group}`` identifies the scenario sampling group (e.g., ``gsa``, ``gsa-l1-low``) and ``{method}`` selects the surrogate type (``pce``, ``rf``, ``xgb``, ``mlp``). All methods consume the same -solved scenarios, and ``sensitivity_analysis.default_surrogate`` selects the -surrogate downstream consumers (notebooks, uncertainty plots) load by -default. +solved scenarios; downstream targets name the desired method explicitly. .. note:: diff --git a/docs/surrogate_modelling.rst b/docs/surrogate_modelling.rst index 05c5ee49..529495de 100644 --- a/docs/surrogate_modelling.rst +++ b/docs/surrogate_modelling.rst @@ -44,9 +44,8 @@ Surrogate methods ----------------- Four methods are supported: ``pce``, ``rf``, ``xgb``, and ``mlp``. The same -solved scenarios can be fitted by several methods independently; which one -downstream consumers load by default is set by -``sensitivity_analysis.default_surrogate`` (currently ``mlp``). +solved scenarios can be fitted by several methods independently. The method is +always explicit in the target or bundle filename. Polynomial Chaos Expansion (PCE) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -260,7 +259,6 @@ same solved scenarios without duplicating the design. sensitivity_analysis: holdout_fraction: 0.15 # fraction reserved for out-of-sample validation threads: 6 # sklearn n_jobs / BLAS threads - default_surrogate: mlp # method downstream consumers load by default discover_scenarios_on_disk: false sobol: # Sobol index settings (see sensitivity_analysis) outputs: [total_cost, co2, ch4, n2o, land_use, yll] @@ -301,8 +299,6 @@ same solved scenarios without duplicating the design. validation (e.g. 0.15). Set to 0 to disable holdout. - ``threads``: Threads for the fit (sklearn ``n_jobs`` and BLAS). Note that ``tools/smk -jN`` clamps rule threads to ``N``. -- ``default_surrogate``: Method downstream consumers (uncertainty-band plots, - notebooks) load when none is named explicitly. Must match a ``methods`` key. - ``discover_scenarios_on_disk``: When ``false`` (default), ``build_surrogate`` declares every Sobol scenario as an input so one ``tools/smk`` call drives the whole solve-analyse-surrogate chain. When ``true``, it instead scans the diff --git a/tests/config/test_sensitivity.yaml b/tests/config/test_sensitivity.yaml deleted file mode 100644 index c8847878..00000000 --- a/tests/config/test_sensitivity.yaml +++ /dev/null @@ -1,86 +0,0 @@ -# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek -# -# SPDX-License-Identifier: CC-BY-4.0 - -# Small test configuration for sensitivity analysis -# Uses test config base (reduced resolution) with minimal GSA samples - -name: "test_sensitivity" -scenarios: - # Small sensitivity test scenarios - # 3 parameters x 32 samples for fast testing - - # Default scenario for reference - default: {} - - _generators: - - name: gsa_{sample_id} - mode: sensitivity - samples: 32 - parameters: - yield_factor: - lower: 0.8 - upper: 1.2 - ch4_factor: - lower: 0.7 - upper: 1.3 - luc_factor: - lower: 0.5 - upper: 1.5 - template: - sensitivity: - crop_yields: - all: "{yield_factor}" - emission_factors: - ch4: "{ch4_factor}" - luc: "{luc_factor}" - -sensitivity_analysis: - holdout_fraction: 0.15 - threads: 2 - default_surrogate: pce - methods: - pce: - grid_resolution: 10 - method_options: - cross_truncation: 0.8 - -# Inherit from test config for fast execution -aggregation: - regions: - target_count: 200 - resource_class_quantiles: [0.5] - irrigated_area_source: "current" - -commodities: - hubs: 14 - -crops: - - wheat - - maize - - soybean - - white-potato - -# Keep this sensitivity fixture focused on its health and emissions parameters. -multiple_cropping: - wheat_maize: null - wheat_soybean: null - maize_soybean: null - -# Enable health for RR sensitivity testing -health: - enabled: true - value_per_yll: 50000 - -# Enable GHG pricing -emissions: - ghg_pricing_enabled: true - ghg_price: 100 - -# Use actual yields for more realistic results -validation: - use_actual_yields: true - -solving: - solver: "highs" - threads: 4 diff --git a/tests/test_gbd2019_rr_appendix.py b/tests/test_gbd2019_rr_appendix.py index 7eeff486..0a4aed79 100644 --- a/tests/test_gbd2019_rr_appendix.py +++ b/tests/test_gbd2019_rr_appendix.py @@ -62,22 +62,15 @@ class TestNormalizeExposure: """Tests for converting exposure text to g/day float.""" def test_basic_g_per_day(self): - assert _normalize_exposure("100 g/day", 1.0) == pytest.approx(100.0) - - def test_g_per_day_with_conversion(self): - assert _normalize_exposure("50 g/day", 2.0) == pytest.approx(100.0) + assert _normalize_exposure("100 g/day") == pytest.approx(100.0) def test_energy_based_raises(self): with pytest.raises(ValueError, match="Energy-based exposures"): - _normalize_exposure("5 %energy/day", 1.0) + _normalize_exposure("5 %energy/day") def test_no_unit_raises(self): with pytest.raises(ValueError, match="Unexpected exposure label"): - _normalize_exposure("100", 1.0) - - def test_none_conversion_raises(self): - with pytest.raises(ValueError, match="Missing conversion factor"): - _normalize_exposure("100 g/day", None) + _normalize_exposure("100") class TestExtractRiskBlocks: @@ -142,7 +135,7 @@ def test_valid_block_with_known_cause(self): for col in range(13, 28): row_data[col] = "0.80 (0.70, 0.90)" df = self._make_mock_df([{0: "Diet low in fruits"}, row_data]) - result = parse_gbd2019_rr_appendix(df, ssb_sugar_per_gram=0.1) + result = parse_gbd2019_rr_appendix(df) assert len(result) == 15 row = result.iloc[0] @@ -160,7 +153,7 @@ def test_unmapped_outcome_is_skipped(self): breast_row[col] = "0.95 (0.90, 1.00)" ihd_row[col] = "0.80 (0.70, 0.90)" df = self._make_mock_df([{0: "Diet low in fruits"}, breast_row, ihd_row]) - result = parse_gbd2019_rr_appendix(df, ssb_sugar_per_gram=0.1) + result = parse_gbd2019_rr_appendix(df) assert len(result) == 15 assert result.iloc[0]["cause"] == "CHD" @@ -169,7 +162,7 @@ def test_output_columns(self): for col in range(13, 28): row_data[col] = "0.80 (0.70, 0.90)" df = self._make_mock_df([{0: "Diet low in fruits"}, row_data]) - result = parse_gbd2019_rr_appendix(df, ssb_sugar_per_gram=0.1) + result = parse_gbd2019_rr_appendix(df) assert set(result.columns) == { "risk_factor", "cause", @@ -180,27 +173,13 @@ def test_output_columns(self): "rr_high", } - def test_sugar_conversion_applied(self): - ssb_sugar_per_gram = 0.11 - row_data = {0: "Diabetes mellitus type 2", 1: "200 g/day"} - for col in range(13, 28): - row_data[col] = "1.30 (1.20, 1.40)" - df = self._make_mock_df( - [{0: "Diet high in sugar-sweetened beverages"}, row_data] - ) - result = parse_gbd2019_rr_appendix(df, ssb_sugar_per_gram=ssb_sugar_per_gram) - assert result.iloc[0]["risk_factor"] == "sugar" - assert result.iloc[0]["exposure_g_per_day"] == pytest.approx( - 200.0 * ssb_sugar_per_gram - ) - def test_no_records_raises(self): row_data = {0: "Breast cancer", 1: "100 g/day"} for col in range(13, 28): row_data[col] = "0.95 (0.90, 1.00)" df = self._make_mock_df([{0: "Diet low in fruits"}, row_data]) with pytest.raises(ValueError, match="No dietary risk records"): - parse_gbd2019_rr_appendix(df, ssb_sugar_per_gram=0.1) + parse_gbd2019_rr_appendix(df) def test_duplicate_records_aggregated(self): row1 = {0: "Ischemic heart disease", 1: "100 g/day"} @@ -209,7 +188,7 @@ def test_duplicate_records_aggregated(self): row1[col] = "0.80 (0.70, 0.90)" row2[col] = "0.90 (0.85, 0.95)" df = self._make_mock_df([{0: "Diet low in fruits"}, row1, row2]) - result = parse_gbd2019_rr_appendix(df, ssb_sugar_per_gram=0.1) + result = parse_gbd2019_rr_appendix(df) assert result.iloc[0]["rr_mean"] == pytest.approx(0.85) assert result.iloc[0]["rr_low"] == pytest.approx(0.775) assert result.iloc[0]["rr_high"] == pytest.approx(0.925) @@ -230,26 +209,15 @@ def test_cause_map_maps_ihme_names(self): assert CAUSE_MAP["Colon and rectum cancer"] == "CRC" def test_risk_config_has_expected_risk_factors(self): - assert {v["risk_factor"] for v in RISK_CONFIG.values()} == { + assert set(RISK_CONFIG.values()) == { "fruits", "vegetables", "whole_grains", "legumes", "nuts_seeds", "red_meat", - "sugar", } def test_risk_config_keys_start_with_diet(self): for key in RISK_CONFIG: assert key.startswith("Diet"), f"Key '{key}' does not start with 'Diet'" - - def test_risk_config_ssb_has_none_conversion(self): - assert ( - RISK_CONFIG["Diet high in sugar-sweetened beverages"]["conversion"] is None - ) - - def test_risk_config_standard_factors_have_unit_conversion(self): - for name, config in RISK_CONFIG.items(): - if config["risk_factor"] != "sugar": - assert config["conversion"] == 1.0, f"Expected 1.0 for {name}" diff --git a/tests/test_sensitivity.py b/tests/test_sensitivity.py index 1f5c9280..fdb3bccb 100644 --- a/tests/test_sensitivity.py +++ b/tests/test_sensitivity.py @@ -579,60 +579,6 @@ def test_food_processing_is_invariant_to_waste(self, mock_network): ) -class TestFoodLossWasteBundle: - def test_bundle_only_config(self, mock_network): - """food_loss_waste applies the same factor to both loss and waste.""" - n = mock_network - orig_wheat = float( - n.links.static.loc["produce:wheat_rainfed:region1", "efficiency"] - ) - orig_consume_eff = float(n.links.static.loc["consume:flour:USA", "efficiency"]) - - cfg = {"food_loss_waste": 1.5} - apply_sensitivity_factors(n, cfg) - - # Loss side: wheat 10% loss -> 15%, mult 0.9 -> 0.85. - np.testing.assert_allclose( - n.links.static.loc["produce:wheat_rainfed:region1", "efficiency"], - _scaled_efficiency(orig_wheat, 0.9, 1.5), - ) - # Waste side: 20% waste -> 30%, mult 0.8 -> 0.7. - np.testing.assert_allclose( - n.links.static.loc["consume:flour:USA", "flw_multiplier"], 0.7 - ) - np.testing.assert_allclose( - n.links.static.loc["consume:flour:USA", "efficiency"], - orig_consume_eff * (0.7 / 0.8), - ) - - def test_component_keys_override_bundle(self, mock_network): - """food_loss / food_waste override food_loss_waste when both set.""" - n = mock_network - orig_wheat = float( - n.links.static.loc["produce:wheat_rainfed:region1", "efficiency"] - ) - orig_consume_eff = float(n.links.static.loc["consume:flour:USA", "efficiency"]) - - # bundle=1.5 would scale both; explicit food_loss=2.0 overrides loss, - # food_waste falls through to the bundle value (1.5). - cfg = {"food_loss_waste": 1.5, "food_loss": 2.0} - apply_sensitivity_factors(n, cfg) - - # Wheat: 10% loss * 2.0 -> 20%, mult 0.9 -> 0.8. - np.testing.assert_allclose( - n.links.static.loc["produce:wheat_rainfed:region1", "efficiency"], - _scaled_efficiency(orig_wheat, 0.9, 2.0), - ) - # Waste uses bundle factor 1.5. - np.testing.assert_allclose( - n.links.static.loc["consume:flour:USA", "flw_multiplier"], 0.7 - ) - np.testing.assert_allclose( - n.links.static.loc["consume:flour:USA", "efficiency"], - orig_consume_eff * (0.7 / 0.8), - ) - - class TestApplyFcrFactor: def test_scales_primary_and_coproduct_food_outputs(self, mock_network): """FCR scales every food-output bus on animal links (incl. co-products).""" diff --git a/workflow/scripts/gbd2019_rr_appendix.py b/workflow/scripts/gbd2019_rr_appendix.py index e1a1298a..0a3fb697 100644 --- a/workflow/scripts/gbd2019_rr_appendix.py +++ b/workflow/scripts/gbd2019_rr_appendix.py @@ -20,39 +20,14 @@ logger = logging.getLogger(__name__) -# Map IHME dietary risk names to model risk_factor identifiers and exposure conversion factors +# Map IHME dietary risk names to model risk-factor identifiers. RISK_CONFIG = { - "Diet low in fruits": {"risk_factor": "fruits", "unit": "g/day", "conversion": 1.0}, - "Diet low in vegetables": { - "risk_factor": "vegetables", - "unit": "g/day", - "conversion": 1.0, - }, - "Diet low in whole grains": { - "risk_factor": "whole_grains", - "unit": "g/day", - "conversion": 1.0, - }, - "Diet low in legumes": { - "risk_factor": "legumes", - "unit": "g/day", - "conversion": 1.0, - }, - "Diet low in nuts and seeds": { - "risk_factor": "nuts_seeds", - "unit": "g/day", - "conversion": 1.0, - }, - "Diet high in red meat": { - "risk_factor": "red_meat", - "unit": "g/day", - "conversion": 1.0, - }, - "Diet high in sugar-sweetened beverages": { - "risk_factor": "sugar", - "unit": "g/day", - "conversion": None, - }, + "Diet low in fruits": "fruits", + "Diet low in vegetables": "vegetables", + "Diet low in whole grains": "whole_grains", + "Diet low in legumes": "legumes", + "Diet low in nuts and seeds": "nuts_seeds", + "Diet high in red meat": "red_meat", } @@ -112,7 +87,7 @@ def _parse_rr_value(cell: object) -> tuple[float, float | None, float | None]: return mean, low, high -def _normalize_exposure(raw: str, conversion: float | None) -> float: +def _normalize_exposure(raw: str) -> float: """Convert exposure text like '100 g/day' into g/day as float.""" parts = raw.strip().split() if len(parts) < 2: @@ -129,10 +104,7 @@ def _normalize_exposure(raw: str, conversion: float | None) -> float: "Energy-based exposures are not supported in the current health module" ) - if conversion is None: - raise ValueError("Missing conversion factor for omega-3 exposure") - - return value * conversion + return value def _extract_risk_blocks(df: pd.DataFrame) -> dict[str, tuple[int, int]]: @@ -153,32 +125,15 @@ def _extract_risk_blocks(df: pd.DataFrame) -> dict[str, tuple[int, int]]: return bounds -def parse_gbd2019_rr_appendix( - df: pd.DataFrame, - ssb_sugar_per_gram: float, - basis_factor_by_risk: dict[str, float] | None = None, -) -> pd.DataFrame: +def parse_gbd2019_rr_appendix(df: pd.DataFrame) -> pd.DataFrame: """Parse the Excel sheet into tidy RR records (all 15 adult age buckets).""" - if basis_factor_by_risk is None: - basis_factor_by_risk = {} - records: list[dict[str, float | str]] = [] skipped_causes: dict[str, set[str]] = {} skipped_units: set[str] = set() blocks = _extract_risk_blocks(df) for risk_name, (start, end) in blocks.items(): - config = RISK_CONFIG[risk_name] - risk_id = config["risk_factor"] - conversion = config["conversion"] - - if risk_id == "sugar": - if ssb_sugar_per_gram <= 0: - raise ValueError("ssb_sugar_per_gram must be positive") - conversion = ssb_sugar_per_gram - - basis_factor = float(basis_factor_by_risk.get(risk_id, 1.0)) - effective_conversion = None if conversion is None else conversion * basis_factor + risk_id = RISK_CONFIG[risk_name] block = df.iloc[start:end] block = block[block[0].notna()] @@ -197,7 +152,7 @@ def parse_gbd2019_rr_appendix( cause = CAUSE_MAP[outcome] try: - exposure = _normalize_exposure(exposure_raw, effective_conversion) + exposure = _normalize_exposure(exposure_raw) except ValueError: skipped_units.add(exposure_raw) continue diff --git a/workflow/scripts/generate_rr_age_attenuation.py b/workflow/scripts/generate_rr_age_attenuation.py index 02005eaf..286bebb4 100644 --- a/workflow/scripts/generate_rr_age_attenuation.py +++ b/workflow/scripts/generate_rr_age_attenuation.py @@ -98,9 +98,7 @@ def _extract_shape(rr19: pd.DataFrame) -> dict[tuple[str, str, str], float]: def main() -> None: - rr19 = parse_gbd2019_rr_appendix( - pd.read_excel(GBD2019_RR_XLSX, header=None), ssb_sugar_per_gram=1.0 - ) + rr19 = parse_gbd2019_rr_appendix(pd.read_excel(GBD2019_RR_XLSX, header=None)) shape = _extract_shape(rr19) rows = [] diff --git a/workflow/scripts/solve_model/sensitivity.py b/workflow/scripts/solve_model/sensitivity.py index bdd7fcc9..03c18495 100644 --- a/workflow/scripts/solve_model/sensitivity.py +++ b/workflow/scripts/solve_model/sensitivity.py @@ -16,8 +16,6 @@ - Food loss (loss_fraction on crop_production / crop_production_multi / animal_production links) - Food waste (waste_fraction on food_consumption links) -- Food loss + waste bundle (food_loss_waste applies the same factor to - both unless an explicit food_loss / food_waste key overrides it) - Feed conversion ratios (efficiency on animal_production links) - Production costs (marginal_cost on crop_production, animal_production) @@ -105,10 +103,6 @@ def apply_sensitivity_factors(n: pypsa.Network, sensitivity_cfg: dict) -> None: - food_waste: float multiplier on the waste_fraction baked into food_consumption efficiencies and the ``flw_multiplier`` metadata column. - - food_loss_waste: float bundle convenience key. Applies the - same factor to both loss and waste, unless an explicit - food_loss or food_waste key is also set (per-component keys - override the bundle). - feed_conversion: float - costs: {crop: float, animal: float} """ @@ -123,12 +117,8 @@ def apply_sensitivity_factors(n: pypsa.Network, sensitivity_cfg: dict) -> None: if emission_cfg: _apply_emission_factors(n, emission_cfg) - # food_loss_waste is a bundle: it sets both food_loss and food_waste - # in one place. Per-component keys (food_loss, food_waste) take - # precedence when both are supplied. - bundle_factor = sensitivity_cfg.get("food_loss_waste", 1.0) - food_loss_factor = sensitivity_cfg.get("food_loss", bundle_factor) - food_waste_factor = sensitivity_cfg.get("food_waste", bundle_factor) + food_loss_factor = sensitivity_cfg.get("food_loss", 1.0) + food_waste_factor = sensitivity_cfg.get("food_waste", 1.0) if food_loss_factor != 1.0: _apply_food_loss_factor(n, food_loss_factor) From 4890c63662276fecb8fbe3c54e1c914e024f0d2c Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:17:14 -0700 Subject: [PATCH 05/10] refactor: simplify configuration validation --- config/schemas/config.schema.yaml | 10 +++- tests/test_config_validation.py | 57 +++++++++++++++++++++++ workflow/rules/common.smk | 1 - workflow/validation/__init__.py | 64 +++++++++++--------------- workflow/validation/consumer_values.py | 56 ++++++++++------------ 5 files changed, 117 insertions(+), 71 deletions(-) create mode 100644 tests/test_config_validation.py diff --git a/config/schemas/config.schema.yaml b/config/schemas/config.schema.yaml index 1a7c35f8..c7e0c284 100644 --- a/config/schemas/config.schema.yaml +++ b/config/schemas/config.schema.yaml @@ -9,6 +9,7 @@ description: | Based on config/default.yaml structure. type: object required: + - scenarios - planning_horizon - baseline_year - currency_base_year @@ -62,6 +63,8 @@ required: - deviation_penalty - solving - remote_solve + - sensitivity + - sensitivity_analysis - plotting additionalProperties: false @@ -171,6 +174,7 @@ properties: - min_crop_yield_t_per_ha - min_grassland_yield_t_per_ha - min_link_area_mha + - min_water_capacity_mm3 - min_water_requirement_m3_per_ha - min_co2_coefficient_tco2_per_ha - min_cost_correction_bnusd @@ -1141,6 +1145,8 @@ properties: - source - baseline_age - anchor_groups_to_gbd + - fbs_override_foods + - nhanes - source_basis - fbs - gdd_ia @@ -1318,12 +1324,14 @@ properties: - region_clusters - breakpoint_rel_tol - log_rr_points + - intake_age_min - value_per_yll - risk_factors - causes - risk_cause_map - gbd_rei_id - gbd_cause_id + - alternative_rr - clustering additionalProperties: false properties: @@ -2074,7 +2082,7 @@ properties: solving: type: object - required: [solver, io_api, threads, calculate_fixed_duals, options_gurobi, options_highs, time_limit, runtime, mem_mb, inline_analysis] + required: [solver, io_api, threads, calculate_fixed_duals, options_gurobi, options_highs, export_for_tuning, time_limit, runtime, mem_mb, inline_analysis] additionalProperties: false properties: solver: diff --git a/tests/test_config_validation.py b/tests/test_config_validation.py new file mode 100644 index 00000000..46bf6a47 --- /dev/null +++ b/tests/test_config_validation.py @@ -0,0 +1,57 @@ +# SPDX-FileCopyrightText: 2026 Koen van Greevenbroek +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from copy import deepcopy + +import pytest +import yaml + +from workflow.validation.consumer_values import validate_consumer_values + + +@pytest.fixture +def default_config() -> dict: + with open("config/default.yaml") as f: + return yaml.safe_load(f) + + +def test_schema_requires_fixed_default_sections(default_config) -> None: + with open("config/schemas/config.schema.yaml") as f: + schema = yaml.safe_load(f) + + assert set(default_config) <= set(schema["required"]) + for section in ("numerics", "diet", "health", "solving"): + section_schema = schema["properties"][section] + assert set(default_config[section]) <= set(section_schema["required"]) + + +def test_disabled_piecewise_utility_needs_no_baseline(default_config) -> None: + validate_consumer_values(default_config) + + +def test_piecewise_utility_requires_configured_baseline(default_config) -> None: + config = deepcopy(default_config) + config["food_utility_piecewise"]["enabled"] = True + + with pytest.raises(ValueError, match="baseline scenario 'baseline' is not defined"): + validate_consumer_values(config) + + +def test_piecewise_utility_requires_enforced_baseline(default_config) -> None: + config = deepcopy(default_config) + config["scenarios"]["baseline"] = {} + config["food_utility_piecewise"]["enabled"] = True + + with pytest.raises(ValueError, match="enforce_baseline_diet=true"): + validate_consumer_values(config) + + +def test_scenario_piecewise_utility_accepts_enforced_baseline(default_config) -> None: + config = deepcopy(default_config) + config["scenarios"] = { + "baseline": {"validation": {"enforce_baseline_diet": True}}, + "utility": {"food_utility_piecewise": {"enabled": True}}, + } + + validate_consumer_values(config) diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index aa664dbc..6049f946 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -174,7 +174,6 @@ def scenario_override_hash(scenario_name): # Extract configuration name and relevant config sections name = config["name"] gaez_cfg = config["data"]["gaez"] -grazing_cfg = config.get("grazing", {}) # Load GAEZ crop code mapping from CSV with open("data/curated/gaez_crop_code_mapping.csv", newline="") as _gaez_mapping_file: diff --git a/workflow/validation/__init__.py b/workflow/validation/__init__.py index 9c1bbc10..8f6fe22b 100644 --- a/workflow/validation/__init__.py +++ b/workflow/validation/__init__.py @@ -4,7 +4,6 @@ """Validation entry points for configuration and data consistency checks.""" -from collections.abc import Iterable from pathlib import Path from typing import Callable @@ -37,38 +36,35 @@ Validator = Callable[[dict, Path], None] -_CHECKS: dict[str, Validator] = { - "config_schema": validate_config_schema, - "calibration": validate_calibration, - "calibration_provenance": validate_calibration_provenance, - "commodities": validate_commodities, - "consumer_values": validate_consumer_values, - "optimal_taxes": validate_optimal_taxes, - "country_regions": validate_country_regions, - "food_groups": validate_food_groups, - "food_basis": validate_food_basis, - "diet_basis": validate_diet_basis, - "faostat_maps": validate_faostat_maps, - "crop_food_pathways": validate_crop_food_pathways, - "crop_groups": validate_crop_groups, - "crop_moisture_content": validate_crop_moisture_content, - "cropgrids_crops": validate_cropgrids_crops, - "gaez_crop_mapping": validate_gaez_crop_mapping, - "seed_rates": validate_seed_rates, - "health_map": validate_health_map, - "m49_codes": validate_m49_codes, - "multi_cropping": validate_multi_cropping, - "nutrition": validate_nutrition, - "sensitivity_generator": validate_sensitivity_generator, - "yield_unit_conversions": validate_yield_unit_conversions, -} +_SEMANTIC_CHECKS: tuple[tuple[str, Validator], ...] = ( + ("calibration", validate_calibration), + ("calibration_provenance", validate_calibration_provenance), + ("commodities", validate_commodities), + ("consumer_values", validate_consumer_values), + ("optimal_taxes", validate_optimal_taxes), + ("country_regions", validate_country_regions), + ("food_groups", validate_food_groups), + ("food_basis", validate_food_basis), + ("diet_basis", validate_diet_basis), + ("faostat_maps", validate_faostat_maps), + ("crop_food_pathways", validate_crop_food_pathways), + ("crop_groups", validate_crop_groups), + ("crop_moisture_content", validate_crop_moisture_content), + ("cropgrids_crops", validate_cropgrids_crops), + ("gaez_crop_mapping", validate_gaez_crop_mapping), + ("seed_rates", validate_seed_rates), + ("health_map", validate_health_map), + ("m49_codes", validate_m49_codes), + ("multi_cropping", validate_multi_cropping), + ("nutrition", validate_nutrition), + ("sensitivity_generator", validate_sensitivity_generator), + ("yield_unit_conversions", validate_yield_unit_conversions), +) def validate( config: dict, project_root: Path | None = None, - *, - enabled_checks: Iterable[str] | None = None, ) -> None: """Run configured validation checks against the active config and data. @@ -78,22 +74,14 @@ def validate( The merged Snakemake configuration dictionary. project_root: Root directory of the repository. Defaults to the current working directory. - enabled_checks: - Optional iterable of check names to run. When omitted, all registered checks - are executed. """ logger.info("Validating configuration and input datasets") root = Path(project_root) if project_root else Path.cwd() - check_names = tuple(enabled_checks) if enabled_checks else tuple(_CHECKS) + validate_config_schema(config, root) errors: list[str] = [] - for name in check_names: - try: - check = _CHECKS[name] - except KeyError as exc: - raise KeyError(f"Unknown validation check '{name}'") from exc - + for name, check in _SEMANTIC_CHECKS: try: check(config, root) except Exception as exc: diff --git a/workflow/validation/consumer_values.py b/workflow/validation/consumer_values.py index 2e4f5b59..52809844 100644 --- a/workflow/validation/consumer_values.py +++ b/workflow/validation/consumer_values.py @@ -4,37 +4,31 @@ """Validation checks for consumer values configuration.""" - -def _consumer_values_enabled(config: dict, scenario_defs: dict) -> bool: - def has_consumer_values_sources(cfg: dict) -> bool: - sources = [str(src) for src in cfg.get("sources", [])] - return any("consumer_values" in src for src in sources) - - base_cfg = config["food_incentives"] - if has_consumer_values_sources(base_cfg) and bool(base_cfg["enabled"]): - return True - - for overrides in scenario_defs.values(): - if not isinstance(overrides, dict): - continue - cv_cfg = overrides.get("food_incentives", {}) - if not isinstance(cv_cfg, dict) or not bool(cv_cfg.get("enabled", False)): - continue - merged_sources = cv_cfg.get("sources", base_cfg.get("sources", [])) - if any("consumer_values" in str(src) for src in merged_sources): - return True - - return False +from workflow.scenario_generators import expand_scenario_defs +from workflow.scripts.solve_namespace import get_effective_config def validate_consumer_values(config: dict, _project_root=None) -> None: - """Ensure consumer values runs have a baseline scenario defined.""" - scenario_defs = config.get("scenarios") or {} - - if not _consumer_values_enabled(config, scenario_defs): - return - - if "baseline" not in scenario_defs: - raise ValueError( - "consumer values incentives enabled but scenarios does not define a 'baseline' scenario" - ) + """Validate piecewise food utility against its configured baseline.""" + scenario_defs = expand_scenario_defs(config["scenarios"]) + effective_configs = [("base config", config)] + effective_configs.extend( + (f"scenario '{name}'", get_effective_config(config, name, scenario_defs)) + for name in scenario_defs + ) + + for label, effective in effective_configs: + if not effective["food_utility_piecewise"]["enabled"]: + continue + baseline = effective["consumer_values"]["baseline_scenario"] + if baseline not in scenario_defs: + raise ValueError( + f"{label} enables food_utility_piecewise but configured consumer " + f"values baseline scenario '{baseline}' is not defined" + ) + baseline_config = get_effective_config(config, baseline, scenario_defs) + if not baseline_config["validation"]["enforce_baseline_diet"]: + raise ValueError( + f"consumer values baseline scenario '{baseline}' must set " + "validation.enforce_baseline_diet=true" + ) From f0b89458ff07d726215c031531d854d47f3e79f8 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:19:20 -0700 Subject: [PATCH 06/10] docs: align interfaces with the current workflow --- docs/configuration.rst | 27 +++++++++++++------------- docs/current_diets.rst | 10 +++++----- docs/development.rst | 8 ++++---- docs/health.rst | 2 +- docs/land_use.rst | 5 +++-- docs/nutrition.rst | 7 +++---- docs/workflow.rst | 6 +++++- workflow/validation/cropgrids_crops.py | 2 +- 8 files changed, 36 insertions(+), 31 deletions(-) diff --git a/docs/configuration.rst b/docs/configuration.rst index 0a76d177..1dd59df2 100644 --- a/docs/configuration.rst +++ b/docs/configuration.rst @@ -335,16 +335,13 @@ Validation Options Set ``validation.enforce_baseline_diet`` to ``true`` to force the optimizer to match baseline consumption derived from the estimated baseline diet. When this flag is active, the ``diet.baseline_age`` and ``baseline_year`` settings determine which -cohort/year is enforced. Use ``validation.food_group_slack_marginal_cost`` to set the -penalty (USD\ :sub:`2024` per Mt) for the slack generators that backstop those fixed -food-group loads. Keep the value high so slack only activates when recorded production -cannot meet the enforced demand targets. +cohort/year is enforced. ``validation.slack_marginal_cost`` sets the common +penalty in bn USD per Mt or Mha for validation slack. Keep the value high so +slack only activates when recorded production cannot meet the enforced +demand targets. Set ``validation.enforce_baseline_feed`` to ``true`` to fix animal feed use to -GLEAM-derived baseline levels (see :ref:`gleam-feed-baseline`). The baseline is -scaled from GLEAM 2.0 (2010) to the reference year and calibrated against the -known GLEAM 3.0 global total using ``validation.gleam_calibration_year`` and -``validation.gleam_calibration_total_gt_dm``. +GLEAM-derived baseline levels (see :ref:`gleam-feed-baseline`). See :doc:`validation` for a detailed walkthrough of the validation workflow and diagnostic figures. @@ -361,7 +358,8 @@ the objective: curve per ``(food, country)`` pair. When ``food_utility_piecewise.enabled`` is ``true``, the workflow always reads -utility blocks from ``results/{name}/consumer_values/utility_blocks.csv``. +utility blocks from +``results/{name}/consumer_values/{baseline_scenario}/utility_blocks.csv``. These blocks are generated by ``calibrate_food_utility_blocks`` from: * baseline dual values extracted by ``extract_consumer_values``; and @@ -421,8 +419,11 @@ across components. * ``deviation_penalty.deviation_type``: ``absolute`` or ``relative``. * ``deviation_penalty.quadratic_cost``: shared coefficient for quadratic mode. * ``deviation_penalty.land.enabled`` plus per-component switches - ``land.crops.enabled``, ``land.grassland.enabled``, ``feed.enabled``, - ``diet.enabled``. ``land.land_conversion.enabled`` (default ``false``) + ``deviation_penalty.land.crops.enabled``, + ``deviation_penalty.land.grassland.enabled``, + ``deviation_penalty.feed.enabled``, and + ``deviation_penalty.diet.enabled``. + ``deviation_penalty.land.land_conversion.enabled`` (default ``false``) would additionally penalise land-use transitions, but is kept off because those carriers include sparing -- the penalty would tax reforestation from a zero baseline. @@ -434,8 +435,8 @@ across components. applied after sentinel resolution; lets scenarios scan around the calibrated central value without hard-coding absolute numbers. * ``deviation_penalty.land.crops.max_relative_deviation`` / - ``land.grassland.max_relative_deviation`` / - ``feed.max_relative_deviation``: hard-mode bounds. + ``deviation_penalty.land.grassland.max_relative_deviation`` / + ``deviation_penalty.feed.max_relative_deviation``: hard-mode bounds. **Behavior notes**: diff --git a/docs/current_diets.rst b/docs/current_diets.rst index 0cd0020d..0311610f 100644 --- a/docs/current_diets.rst +++ b/docs/current_diets.rst @@ -128,7 +128,7 @@ basis, so no further conversion is needed when reading before the conversion. GBD exposure is converted to the model basis at load time via -``diet.source_basis`` and ``diet.weight_conversion`` (cooked→dry for +``diet.source_basis`` and ``weight_conversion`` (cooked→dry for ``whole_grains`` and ``legumes`` at 0.45 and 0.40; cooked→fresh for ``red_meat`` at 1.43). NHANES values are intake-based and pass through unchanged. @@ -470,7 +470,7 @@ source (or NHANES for the USA). GBD exposure is converted to the model's basis at load time, per food-group, using ``diet.source_basis`` plus per-(source, country, food_group) overrides from ``data/curated/diet_source_basis_overrides.csv`` -and the conversion tables in ``diet.weight_conversion``. The script also +and the conversion tables in ``weight_conversion``. The script also logs cross-validation metrics: median and range of the source/GBD ratio across countries for every risk group, and GBD's milk exposure as a cross-check on the dairy total. @@ -766,12 +766,12 @@ Downstream Uses ``config.validation.enforce_baseline_diet`` is true, the solver adds per-food, per-country equality constraints on food consumption links. * **Within-group ratio fixing**: when - ``config.food_groups.fix_within_group_ratios`` is true, foods within + ``config.food_groups.fix_within_group_ratios.enabled`` is true, foods within each group are constrained to keep their baseline proportions while group totals may vary. * **Piecewise consumer utility calibration**: baseline per-food consumption and baseline food-equality duals together calibrate - ``results/{name}/consumer_values/utility_blocks.csv`` + ``results/{name}/consumer_values/{baseline_scenario}/utility_blocks.csv`` (:doc:`consumer_values`). * **Health impact assessment**: baseline consumption feeds the population-attributable fraction calculation (:doc:`health`). @@ -841,7 +841,7 @@ Workflow Integration (default ``"All ages"``). * ``config.diet.fbs_override_foods`` — foods anchored to FBS supply. See :ref:`Why animal products use FBS `. -* ``config.diet.source_basis`` and ``config.diet.weight_conversion`` — +* ``config.diet.source_basis`` and ``config.weight_conversion`` — per-source native bases and conversion tables. * ``config.diet.gdd_ia.cooked_to_raw`` — per-group cooked→raw inflation factors for GDD-IA (currently ``red_meat: 1.43``). diff --git a/docs/development.rst b/docs/development.rst index e85d5185..bab61d76 100644 --- a/docs/development.rst +++ b/docs/development.rst @@ -172,10 +172,10 @@ The project uses **pytest** for integration testing via the Snakemake Python API Test Configuration ~~~~~~~~~~~~~~~~~~ -Two dedicated config files drive the test suite: - -* **``tests/config/test.yaml``**: Minimal overrides on top of ``default.yaml`` — 200 regions, 2 resource classes, 9 crops, 14 trade hubs. Outputs to ``results/test/``. -* **``tests/config/test_scenarios.yaml``**: Two scenarios (``default`` and ``G``) to exercise the scenario mechanism and GHG pricing code path. +**``tests/config/test.yaml``** provides minimal overrides on top of +``default.yaml``: 200 regions, 2 resource classes, 11 crops, 14 trade hubs, and +two scenarios (``default`` and ``G``) that exercise scenario resolution and GHG +pricing. Outputs are written to ``results/test/``. Running Tests ~~~~~~~~~~~~~ diff --git a/docs/health.rst b/docs/health.rst index 283e47a8..467bddeb 100644 --- a/docs/health.rst +++ b/docs/health.rst @@ -341,7 +341,7 @@ Appendix 1, p. 171). - **Intake units**: Per-group basis matches the baseline-diet pipeline output (model basis; see :doc:`current_diets`). GBD exposure is converted to that basis at load time via ``diet.source_basis`` and - ``diet.weight_conversion``. + ``weight_conversion``. - **Alternative RR sources**: The ``health.alternative_rr`` config option allows substituting GBD dose-response curves with log-linear curves from literature meta-analyses on a per-risk-factor basis. By default, red meat uses literature diff --git a/docs/land_use.rst b/docs/land_use.rst index 01996152..5f0bddfd 100644 --- a/docs/land_use.rst +++ b/docs/land_use.rst @@ -339,8 +339,9 @@ Land Slack Validation runs that pin observed harvested area may encounter land-class mismatches. To maintain feasibility without globally loosening land limits, each land bus can receive a ``land_slack`` generator: - Controlled by ``validation.land_slack: true`` -- Marginal cost set by ``land.slack_marginal_cost`` (USD per Mha) -- Default ~5000 USD/ha ensures slack activates only as a last resort +- Marginal cost set by ``validation.slack_marginal_cost`` (bn USD per Mha) +- The default 50 bn USD/Mha (50,000 USD/ha) ensures slack activates only as a + last resort Multi-Cropping Land Correction ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/nutrition.rst b/docs/nutrition.rst index 859174a7..e3d0a4ac 100644 --- a/docs/nutrition.rst +++ b/docs/nutrition.rst @@ -146,14 +146,13 @@ to a per-country store instead of a nutrient mass: with operator chosen from ``min``/``max``/``equal`` in the config. As with macronutrients, an ``equal`` bound silences ``min``/``max`` for that group. -4. **Per-country equality from the baseline diet.** When the diet - module is configured to anchor a group to current per-country - consumption (``diet.enforce_baseline`` or an equality CSV), the +4. **Per-country equality from an explicit source.** When + ``food_groups.equal_by_country_source`` names an equality CSV, the solver builds a ``per_country_equal`` mapping ``{group: {country: g/person/day}}`` from the baseline diet and feeds it to ``add_food_group_constraints``. The equality RHS then uses the country-specific value instead of a global one — useful - when the goal is *to hold today's group mix fixed and let the model + when the goal is *to hold a country's group mix fixed and let the model choose within-group composition*. This setup keeps dietary diversity decoupled from macronutrient diff --git a/docs/workflow.rst b/docs/workflow.rst index 529206c5..bfe57942 100644 --- a/docs/workflow.rst +++ b/docs/workflow.rst @@ -24,7 +24,11 @@ Paths shown below use the default roots (``processing/``, ``results/``, Validation Hook --------------- -Before Snakemake resolves any rules, the ``workflow/Snakefile`` uses the ``onstart`` hook to run configuration/data validation powered by `Pydantic `__ and `Pandera `__. The checks live in ``workflow/validation/`` and currently ensure, for example, that every category in ``data/curated/food_groups.csv`` is listed under ``config.food_groups.included``. Add new validators by dropping another module in that package and registering it in ``workflow/validation/__init__.py``—the hook aggregates all errors and aborts the workflow if any check fails. +Before Snakemake includes any rule files, ``workflow/Snakefile`` directly +validates the merged configuration and input data. JSON Schema checks the +configuration structure first; semantic checks in ``workflow/validation/`` use +Pandera and focused Python validators. Semantic errors are aggregated before +the workflow aborts. The complete workflow dependency graph is shown below. Each node represents a Snakemake rule, and edges show dependencies between rules. diff --git a/workflow/validation/cropgrids_crops.py b/workflow/validation/cropgrids_crops.py index a340e57a..432882db 100644 --- a/workflow/validation/cropgrids_crops.py +++ b/workflow/validation/cropgrids_crops.py @@ -30,7 +30,7 @@ def validate_cropgrids_crops(config: dict, project_root: Path) -> None: 4. Every entry has a row in ``data/curated/cropgrids_crop_mapping.csv`` with a non-empty ``cropgrids_name`` (so the CROPGRIDS NetCDF can be extracted), a non-empty ``faostat_qcl_item_code`` (so FAOSTAT yield - can be looked up), and a non-empty ``faostat_qcl_yield_element_code``. + can be looked up). 5. No entry appears in ``data/curated/gaez_crop_code_mapping.csv``: it is reserved for GAEZ-backed crops, and a stray entry would otherwise feed the GAEZ download rules. From 7f4e6344a3c55a17ff9112b62f115992c1d6542f Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:20:58 -0700 Subject: [PATCH 07/10] refactor: unify analysis invocation paths --- tools/cluster-solve | 40 +-------- workflow/rules/analysis.smk | 4 +- workflow/scripts/analysis/analyze_model.py | 94 +++++++++------------ workflow/scripts/solve_and_analyze_model.py | 43 ++-------- 4 files changed, 51 insertions(+), 130 deletions(-) diff --git a/tools/cluster-solve b/tools/cluster-solve index 7a7f6b16..322661e5 100755 --- a/tools/cluster-solve +++ b/tools/cluster-solve @@ -64,43 +64,11 @@ def run_scenario( if inline_analysis: # Run analysis in-memory (same as solve_and_analyze_model.py) - from workflow.scripts.analysis.analyze_model import run_analysis - - output_paths = { - attr: getattr(smk.output, attr) - for attr in vars(smk.output) - if isinstance(getattr(smk.output, attr), str) - and getattr(smk.output, attr).endswith(".parquet") - } - - health_enabled = bool(smk.params.health_enabled) - run_analysis( - n, - output_paths=output_paths, - food_groups_path=smk.input.food_groups, - m49_codes_path=smk.input.m49, - population_path=smk.input.population, - ghg_price=float(smk.params.ghg_price), - ch4_gwp=float(smk.params.ch4_gwp), - n2o_gwp=float(smk.params.n2o_gwp), - value_per_yll=float(smk.params.health_value_per_yll), - health_risk_factors=list(smk.params.health_risk_factors), - logger=logger, - health_enabled=health_enabled, - risk_breakpoints_path=( - smk.input.health_risk_breakpoints if health_enabled else None - ), - health_cluster_cause_path=( - smk.input.health_cluster_cause if health_enabled else None - ), - health_cause_log_path=( - smk.input.health_cause_log if health_enabled else None - ), - health_clusters_path=( - smk.input.health_clusters if health_enabled else None - ), - tmrel_path=smk.input.health_tmrel if health_enabled else None, + from workflow.scripts.analysis.analyze_model import ( + run_analysis_from_namespace, ) + + run_analysis_from_namespace(smk, n, logger) else: # Write solved network to disk netcdf_config = smk.params.netcdf diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 89355190..a3bf5ae5 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -204,7 +204,7 @@ else: ), network="/{name}/solved/model_scen-{scenario}.nc", food_groups="data/curated/food_groups.csv", - m49_codes="data/curated/M49-codes.csv", + m49="data/curated/M49-codes.csv", population="/{name}/population.csv", analysis_scripts=_ANALYSIS_SCRIPTS, params: @@ -216,7 +216,7 @@ else: health_enabled=lambda w: get_effective_config(w.scenario)["health"][ "enabled" ], - value_per_yll=lambda w: get_effective_config(w.scenario)["health"][ + health_value_per_yll=lambda w: get_effective_config(w.scenario)["health"][ "value_per_yll" ], health_risk_factors=config["health"]["risk_factors"], diff --git a/workflow/scripts/analysis/analyze_model.py b/workflow/scripts/analysis/analyze_model.py index 000f6f54..d8c54416 100644 --- a/workflow/scripts/analysis/analyze_model.py +++ b/workflow/scripts/analysis/analyze_model.py @@ -57,23 +57,17 @@ from workflow.scripts.solve_namespace import ANALYSIS_OUTPUT_NAMES -def write_empty_outputs(output) -> None: - """Write empty Parquet files for all declared outputs. +def analysis_output_paths(output) -> dict[str, str]: + """Return the canonical analysis output paths from a namespace.""" + return {name: str(getattr(output, name)) for name in ANALYSIS_OUTPUT_NAMES} - Parameters - ---------- - output - Snakemake output object (or any object whose non-underscore string - attributes ending in ``.parquet`` should be created as empty files). - """ - for attr in dir(output): - if attr.startswith("_"): - continue - path = getattr(output, attr, None) - if isinstance(path, str) and path.endswith(".parquet"): - p = Path(path) - p.parent.mkdir(parents=True, exist_ok=True) - pd.DataFrame().to_parquet(p) + +def write_empty_outputs(output) -> None: + """Write empty Parquet files for all canonical analysis outputs.""" + for path in analysis_output_paths(output).values(): + p = Path(path) + p.parent.mkdir(parents=True, exist_ok=True) + pd.DataFrame().to_parquet(p) def run_analysis( @@ -239,6 +233,34 @@ def run_analysis( logger.info("Wrote all analysis outputs to %s", output_dir) +def run_analysis_from_namespace(smk, n: pypsa.Network, logger: logging.Logger) -> None: + """Run analysis using the shared Snakemake/manifest namespace contract.""" + health_enabled = bool(smk.params.health_enabled) + run_analysis( + n, + output_paths=analysis_output_paths(smk.output), + food_groups_path=smk.input.food_groups, + m49_codes_path=smk.input.m49, + population_path=smk.input.population, + ghg_price=float(smk.params.ghg_price), + ch4_gwp=float(smk.params.ch4_gwp), + n2o_gwp=float(smk.params.n2o_gwp), + value_per_yll=float(smk.params.health_value_per_yll), + health_risk_factors=list(smk.params.health_risk_factors), + logger=logger, + health_enabled=health_enabled, + risk_breakpoints_path=( + smk.input.health_risk_breakpoints if health_enabled else None + ), + health_cluster_cause_path=( + smk.input.health_cluster_cause if health_enabled else None + ), + health_cause_log_path=(smk.input.health_cause_log if health_enabled else None), + health_clusters_path=smk.input.health_clusters if health_enabled else None, + tmrel_path=smk.input.health_tmrel if health_enabled else None, + ) + + def main() -> None: """Snakemake entry point: load solved network and run analysis.""" logger = setup_script_logging(snakemake.log[0]) @@ -253,7 +275,7 @@ def main() -> None: try: n = load_solved_network(snakemake.input.network) - except (KeyError, Exception) as e: + except Exception as e: logger.warning( "Failed to load network (%s) — likely an unsolved model. " "Writing empty outputs.", @@ -271,43 +293,7 @@ def main() -> None: logger.info("Loaded network with %d links", len(n.links)) - # Build output_paths dict from snakemake.output - output_paths = { - attr: getattr(snakemake.output, attr) - for attr in dir(snakemake.output) - if not attr.startswith("_") - and isinstance(getattr(snakemake.output, attr), str) - and getattr(snakemake.output, attr).endswith(".parquet") - } - - health_enabled = bool(snakemake.params.health_enabled) - run_analysis( - n, - output_paths=output_paths, - food_groups_path=snakemake.input.food_groups, - m49_codes_path=snakemake.input.m49_codes, - population_path=snakemake.input.population, - ghg_price=float(snakemake.params.ghg_price), - ch4_gwp=float(snakemake.params.ch4_gwp), - n2o_gwp=float(snakemake.params.n2o_gwp), - value_per_yll=float(snakemake.params.value_per_yll), - health_risk_factors=list(snakemake.params.health_risk_factors), - logger=logger, - health_enabled=health_enabled, - risk_breakpoints_path=( - snakemake.input.health_risk_breakpoints if health_enabled else None - ), - health_cluster_cause_path=( - snakemake.input.health_cluster_cause if health_enabled else None - ), - health_cause_log_path=( - snakemake.input.health_cause_log if health_enabled else None - ), - health_clusters_path=( - snakemake.input.health_clusters if health_enabled else None - ), - tmrel_path=snakemake.input.health_tmrel if health_enabled else None, - ) + run_analysis_from_namespace(snakemake, n, logger) if __name__ == "__main__": diff --git a/workflow/scripts/solve_and_analyze_model.py b/workflow/scripts/solve_and_analyze_model.py index 41545063..019d6a36 100644 --- a/workflow/scripts/solve_and_analyze_model.py +++ b/workflow/scripts/solve_and_analyze_model.py @@ -13,7 +13,10 @@ import logging -from workflow.scripts.analysis.analyze_model import run_analysis, write_empty_outputs +from workflow.scripts.analysis.analyze_model import ( + run_analysis_from_namespace, + write_empty_outputs, +) from workflow.scripts.logging_config import setup_script_logging from workflow.scripts.solve_model.core import _ShadowPriceLogFilter, run_solve @@ -30,43 +33,7 @@ def main() -> None: write_empty_outputs(snakemake.output) return - # Phase 2: Analyze (using in-memory solved network) - output_paths = { - attr: getattr(snakemake.output, attr) - for attr in dir(snakemake.output) - if not attr.startswith("_") - and isinstance(getattr(snakemake.output, attr), str) - and getattr(snakemake.output, attr).endswith(".parquet") - } - - health_enabled = bool(snakemake.params.health_enabled) - run_analysis( - n, - output_paths=output_paths, - food_groups_path=snakemake.input.food_groups, - m49_codes_path=snakemake.input.m49, - population_path=snakemake.input.population, - ghg_price=float(snakemake.params.ghg_price), - ch4_gwp=float(snakemake.params.ch4_gwp), - n2o_gwp=float(snakemake.params.n2o_gwp), - value_per_yll=float(snakemake.params.health_value_per_yll), - health_risk_factors=list(snakemake.params.health_risk_factors), - logger=logger, - health_enabled=health_enabled, - risk_breakpoints_path=( - snakemake.input.health_risk_breakpoints if health_enabled else None - ), - health_cluster_cause_path=( - snakemake.input.health_cluster_cause if health_enabled else None - ), - health_cause_log_path=( - snakemake.input.health_cause_log if health_enabled else None - ), - health_clusters_path=( - snakemake.input.health_clusters if health_enabled else None - ), - tmrel_path=snakemake.input.health_tmrel if health_enabled else None, - ) + run_analysis_from_namespace(snakemake, n, logger) logger.info("Solve-and-analyze complete.") From 099425288c68c9f6cea35840ef1cb782419c44d7 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:22:17 -0700 Subject: [PATCH 08/10] fix: reject ignored scenario overrides --- tests/test_solve_namespace.py | 27 +++++++++++++++++++++++++++ workflow/rules/analysis.smk | 6 +++--- workflow/rules/model.smk | 6 +++--- workflow/scripts/solve_namespace.py | 16 +++++++++++----- 4 files changed, 44 insertions(+), 11 deletions(-) diff --git a/tests/test_solve_namespace.py b/tests/test_solve_namespace.py index 9be03129..b3e8572f 100644 --- a/tests/test_solve_namespace.py +++ b/tests/test_solve_namespace.py @@ -74,6 +74,33 @@ def test_rejects_unknown_sensitivity_key(self, base_config): with pytest.raises(ValueError, match="stale"): validate_scenario_config_schemas(base_config, defs, ".") + @pytest.mark.parametrize( + "override", + [ + {"solving": {"inline_analysis": True}}, + {"remote_solve": {"enabled": True}}, + {"plotting": {"comparison_scenarios": "all"}}, + {"food_groups": {"max_per_capita": {"grain": 100.0}}}, + {"biomass": {"marginal_values_usd_per_tonne": 100.0}}, + ], + ) + def test_rejects_parse_or_build_time_override(self, override): + with pytest.raises(ValueError, match="structural key"): + validate_scenario_overrides({"ignored": override}) + + def test_accepts_solver_resource_override(self): + validate_scenario_overrides( + { + "parallel": { + "solving": { + "threads": 4, + "runtime": 20, + "mem_mb": 8000, + } + } + } + ) + def test_validates_one_representative_per_structure(self, base_config, monkeypatch): """Thousands of same-template samples must cost one validation.""" import workflow.validation.config_schema as cs diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index a3bf5ae5..b38b8195 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -163,9 +163,9 @@ if config["solving"]["inline_analysis"]: "max_feed_fraction_by_region" ], countries=config["countries"], - export_for_tuning=lambda w: get_effective_config(w.scenario)[ - "solving" - ].get("export_for_tuning", False), + export_for_tuning=lambda w: get_effective_config(w.scenario)["solving"][ + "export_for_tuning" + ], netcdf=lambda w: get_effective_config(w.scenario)["netcdf"], scenario_hash=lambda w: scenario_override_hash(w.scenario), # --- analysis params --- diff --git a/workflow/rules/model.smk b/workflow/rules/model.smk index c3dfd2b0..67c85391 100644 --- a/workflow/rules/model.smk +++ b/workflow/rules/model.smk @@ -434,9 +434,9 @@ rule solve_model: "max_feed_fraction_by_region" ], countries=config["countries"], - export_for_tuning=lambda w: get_effective_config(w.scenario)["solving"].get( - "export_for_tuning", False - ), + export_for_tuning=lambda w: get_effective_config(w.scenario)["solving"][ + "export_for_tuning" + ], # Only used to force correct reruns when scenario definitions change. scenario_hash=lambda w: scenario_override_hash(w.scenario), output: diff --git a/workflow/scripts/solve_namespace.py b/workflow/scripts/solve_namespace.py index 0cd64d30..1708e62f 100644 --- a/workflow/scripts/solve_namespace.py +++ b/workflow/scripts/solve_namespace.py @@ -38,6 +38,7 @@ "health.relax_and_fix_max_gap", "validation.enforce_baseline_diet", "validation.animal_growth_cap", + "validation.crop_growth_cap", "deviation_penalty", "macronutrients", "food_utility_piecewise", @@ -45,8 +46,6 @@ "food_groups.constraints", "food_groups.fix_within_group_ratios", "food_groups.equal_by_country_source", - "food_groups.max_per_capita", - "biomass.marginal_values_usd_per_tonne", "biomass.biofuel_demand_scale", "land.regional_limit", "land.reforestation_cap", @@ -54,9 +53,16 @@ "exogenous_feed_calibration.enabled", "consumer_values", "sensitivity", - "solving", - "plotting", - "remote_solve", + "solving.solver", + "solving.io_api", + "solving.threads", + "solving.calculate_fixed_duals", + "solving.options_gurobi", + "solving.options_highs", + "solving.export_for_tuning", + "solving.time_limit", + "solving.runtime", + "solving.mem_mb", "netcdf", } From 5ec45e23f7148f8f77d26d1859c4fcb5060d0b99 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:24:51 -0700 Subject: [PATCH 09/10] refactor: require complete configuration --- CHANGELOG.md | 3 +++ config/default.yaml | 2 ++ config/schemas/config.schema.yaml | 3 ++- workflow/rules/analysis.smk | 5 ++--- workflow/rules/common.smk | 8 ++++---- workflow/rules/crops.smk | 4 ++-- workflow/rules/optimal_taxes.smk | 2 +- workflow/rules/plotting.smk | 11 +++++------ workflow/scripts/build_model.py | 2 +- .../doc_figures/validation_food_group_slack.py | 4 +--- workflow/scripts/snakemake_utils.py | 2 +- workflow/scripts/solve_namespace.py | 10 +++++----- workflow/validation/calibration.py | 2 +- workflow/validation/crop_food_pathways.py | 2 +- workflow/validation/crop_groups.py | 2 +- workflow/validation/cropgrids_crops.py | 2 +- workflow/validation/gaez_crop_mapping.py | 2 +- workflow/validation/health_map.py | 4 ++-- workflow/validation/optimal_taxes.py | 2 +- workflow/validation/sensitivity_generator.py | 2 +- 20 files changed, 38 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2bb6bcc2..146f619f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -95,6 +95,9 @@ introduce breaking changes to configuration and outputs. re-stamps a set without solving, for when a code change provably cannot move the artefacts. +- `config/default.yaml` now has the canonical name `default`. All shipped + configuration fields are validated as required; workflow code no longer + supplies hidden fallback values for missing keys. - `planning_horizon` now defaults to 2020, matching `baseline_year`, so an unmodified run solves the observed year the calibration artefacts are fit against. Configs that previously relied on the 2030 default (including diff --git a/config/default.yaml b/config/default.yaml index 7c852597..ff32fc59 100644 --- a/config/default.yaml +++ b/config/default.yaml @@ -3,6 +3,8 @@ # SPDX-License-Identifier: CC-BY-4.0 # Default configuration - override in project-specific configs +name: "default" + # --- section: scenario_metadata --- scenarios: # Each key represents a named scenario that can be activated via the diff --git a/config/schemas/config.schema.yaml b/config/schemas/config.schema.yaml index c7e0c284..8671db80 100644 --- a/config/schemas/config.schema.yaml +++ b/config/schemas/config.schema.yaml @@ -9,6 +9,7 @@ description: | Based on config/default.yaml structure. type: object required: + - name - scenarios - planning_horizon - baseline_year @@ -72,7 +73,7 @@ properties: name: type: string pattern: "^[a-zA-Z0-9_-]+$" - description: "Optional config name identifier (overrides default)" + description: "Config name used for processing and result directories" scenarios: type: object diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index b38b8195..20837933 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -245,7 +245,7 @@ def _sensitivity_generator_group(gen): def _sensitivity_generator(wildcards): """Return the sensitivity generator whose group matches the wildcard.""" - raw_defs = config.get("scenarios") or {} + raw_defs = config["scenarios"] group = wildcards.group generators = [ @@ -358,8 +358,7 @@ def _sensitivity_scenario_inputs(wildcards): def _sensitivity_method_config(wildcards): """Return the method-specific config dict from sensitivity_analysis.methods.""" method = wildcards.method - sa_cfg = config.get("sensitivity_analysis", {}) - methods = sa_cfg.get("methods", {}) + methods = config["sensitivity_analysis"]["methods"] if method not in methods: raise ValueError( f"Unknown sensitivity method '{method}'. " diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 6049f946..fe13ea06 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -49,7 +49,7 @@ def load_scenario_defs(): """Load scenario definitions from the config's `scenarios` key.""" global _SCENARIO_CACHE if _SCENARIO_CACHE is None: - raw_defs = config.get("scenarios") or {} + raw_defs = config["scenarios"] _SCENARIO_CACHE = expand_scenario_defs(raw_defs) return _SCENARIO_CACHE @@ -233,7 +233,7 @@ def gaez_crops(crops=None): crops bypass GAEZ entirely and must not appear in those input lists. """ base = list(crops) if crops is not None else list(config["crops"]) - cropgrids_set = set(config.get("cropgrids_crops") or []) + cropgrids_set = set(config["cropgrids_crops"]) return [c for c in base if c not in cropgrids_set] @@ -249,7 +249,7 @@ def irrigated_crops(): base = list(config["crops"]) else: base = list(irr_cfg) - cropgrids_set = set(config.get("cropgrids_crops") or []) + cropgrids_set = set(config["cropgrids_crops"]) return [c for c in base if c not in cropgrids_set] @@ -260,7 +260,7 @@ def gaez_path(kind: str, water_supply: str, crop: str) -> str: water_supply: "i" (irrigated) or "r" (rainfed) crop: crop name (e.g., "wheat") """ - cropgrids_set = set(config.get("cropgrids_crops") or []) + cropgrids_set = set(config["cropgrids_crops"]) if crop in cropgrids_set: raise ValueError( f"gaez_path() called for CROPGRIDS-backed crop '{crop}'; this " diff --git a/workflow/rules/crops.smk b/workflow/rules/crops.smk index bc28a366..3dfeff80 100644 --- a/workflow/rules/crops.smk +++ b/workflow/rules/crops.smk @@ -224,7 +224,7 @@ rule build_crop_yields_cropgrids: wildcard_constraints: crop=( "|".join(config["cropgrids_crops"]) - if config.get("cropgrids_crops") + if config["cropgrids_crops"] else "__never__" ), group: @@ -890,7 +890,7 @@ rule prepare_fiber_baseline: if config["cost_calibration"]["generate"]: _cal_scenario = config["cost_calibration"]["scenario"] - _cal_name = config.get("name", "default") + _cal_name = config["name"] rule extract_cost_calibration: input: diff --git a/workflow/rules/optimal_taxes.smk b/workflow/rules/optimal_taxes.smk index ad3d82a7..217ad2e4 100644 --- a/workflow/rules/optimal_taxes.smk +++ b/workflow/rules/optimal_taxes.smk @@ -12,7 +12,7 @@ This workflow: 5. Resolves the model with taxes/subsidies applied in the objective """ -plotting_cfg = config.get("plotting", {}) +plotting_cfg = config["plotting"] rule extract_optimal_consumption: diff --git a/workflow/rules/plotting.smk b/workflow/rules/plotting.smk index 6e114ccc..368db58f 100644 --- a/workflow/rules/plotting.smk +++ b/workflow/rules/plotting.smk @@ -4,9 +4,9 @@ gaez = config["data"]["gaez"] -plotting_cfg = config.get("plotting", {}) -food_group_colors = plotting_cfg.get("colors", {}).get("food_groups", {}) -_param_groups_cfg = plotting_cfg.get("colors", {}).get("parameter_groups", {}) +plotting_cfg = config["plotting"] +food_group_colors = plotting_cfg["colors"]["food_groups"] +_param_groups_cfg = plotting_cfg["colors"]["parameter_groups"] parameter_colors = { param: color for group in _param_groups_cfg.values() @@ -78,8 +78,7 @@ def _sobol_sensitivity_groups(): def _sobol_sensitivity_methods(): """Return all configured surrogate method names.""" - sa_cfg = config.get("sensitivity_analysis", {}) - return list(sa_cfg.get("methods", {}).keys()) + return list(config["sensitivity_analysis"]["methods"]) def _sobol_plot_targets(): @@ -574,7 +573,7 @@ rule plot_food_consumption_baseline_map: pdf="/{name}/plots/scen-{scenario}/food_consumption_baseline_map.pdf", csv="/{name}/plots/scen-{scenario}/food_consumption_baseline_map.csv", params: - age=config.get("diet", {}).get("baseline_age", "All ages"), + age=config["diet"]["baseline_age"], reference_year=config["baseline_year"], group_colors=food_group_colors, group: diff --git a/workflow/scripts/build_model.py b/workflow/scripts/build_model.py index 1150fdaa..6b56ac86 100644 --- a/workflow/scripts/build_model.py +++ b/workflow/scripts/build_model.py @@ -229,7 +229,7 @@ def filter(self, record: logging.LogRecord) -> bool: else: expected_irrigated_crops = set(map(str, irrigation_cfg)) # CROPGRIDS-backed crops are rainfed-only by construction. - cropgrids_crops = set(snakemake.config.get("cropgrids_crops") or []) # type: ignore[index] + cropgrids_crops = set(snakemake.config["cropgrids_crops"]) # type: ignore[index] expected_irrigated_crops -= cropgrids_crops # Read yields data and harvested area for each crop and water supply. diff --git a/workflow/scripts/doc_figures/validation_food_group_slack.py b/workflow/scripts/doc_figures/validation_food_group_slack.py index b4617634..c0d6c598 100644 --- a/workflow/scripts/doc_figures/validation_food_group_slack.py +++ b/workflow/scripts/doc_figures/validation_food_group_slack.py @@ -206,9 +206,7 @@ def main() -> None: demand = _aggregate_demand_by_group(network) consumption = _aggregate_consumption_by_group(network) - group_colors = ( - snakemake.config.get("plotting", {}).get("colors", {}).get("food_groups", {}) - ) + group_colors = snakemake.config["plotting"]["colors"]["food_groups"] _plot_two_panel( slack_df, diff --git a/workflow/scripts/snakemake_utils.py b/workflow/scripts/snakemake_utils.py index f8807fa9..8242d70a 100644 --- a/workflow/scripts/snakemake_utils.py +++ b/workflow/scripts/snakemake_utils.py @@ -91,7 +91,7 @@ def _recursive_update(target: dict, source: dict, _path: tuple[str, ...] = ()) - def load_scenarios(config: dict) -> dict: """Load scenario definitions from the config's `scenarios` key.""" - raw_defs = config.get("scenarios") or {} + raw_defs = config["scenarios"] return expand_scenario_defs(raw_defs) diff --git a/workflow/scripts/solve_namespace.py b/workflow/scripts/solve_namespace.py index 1708e62f..7faba875 100644 --- a/workflow/scripts/solve_namespace.py +++ b/workflow/scripts/solve_namespace.py @@ -320,12 +320,12 @@ def resolve_pathvars(path: str, path_roots: dict[str, str]) -> str: def default_path_roots(config: dict) -> dict[str, str]: """Resolve the standard four path roots from a config dict.""" - paths_cfg = config.get("paths", {}) or {} + paths_cfg = config["paths"] return { - "results": resolve_path_root(paths_cfg.get("results_root", "results")), - "processing": resolve_path_root(paths_cfg.get("processing_root", "processing")), - "logs": resolve_path_root(paths_cfg.get("logs_root", "logs")), - "benchmarks": resolve_path_root(paths_cfg.get("benchmarks_root", "benchmarks")), + "results": resolve_path_root(paths_cfg["results_root"]), + "processing": resolve_path_root(paths_cfg["processing_root"]), + "logs": resolve_path_root(paths_cfg["logs_root"]), + "benchmarks": resolve_path_root(paths_cfg["benchmarks_root"]), } diff --git a/workflow/validation/calibration.py b/workflow/validation/calibration.py index 9b74b819..5472c35b 100644 --- a/workflow/validation/calibration.py +++ b/workflow/validation/calibration.py @@ -84,7 +84,7 @@ def _resolve(config: dict, path: tuple) -> dict: def validate_calibration(config: dict, project_root: Path | None = None) -> None: """Ensure calibration sections use the canonical enabled/generate pattern.""" root = Path(project_root) if project_root else Path.cwd() - scenario_names = set((config.get("scenarios") or {}).keys()) + scenario_names = set(config["scenarios"]) errors: list[str] = [] for section in _CALIBRATION_SECTIONS: diff --git a/workflow/validation/crop_food_pathways.py b/workflow/validation/crop_food_pathways.py index 7a58f8ff..54ed309a 100644 --- a/workflow/validation/crop_food_pathways.py +++ b/workflow/validation/crop_food_pathways.py @@ -85,7 +85,7 @@ def validate_crop_food_pathways(config: dict, project_root: Path) -> None: ) # Check that byproducts listed in config appear as foods - byproducts_cfg = set(config.get("byproducts", [])) + byproducts_cfg = set(config["byproducts"]) foods_in_csv = set(df["food"].unique()) missing_byproducts = sorted(byproducts_cfg - foods_in_csv) diff --git a/workflow/validation/crop_groups.py b/workflow/validation/crop_groups.py index ca5390c2..6e456bfc 100644 --- a/workflow/validation/crop_groups.py +++ b/workflow/validation/crop_groups.py @@ -17,7 +17,7 @@ def validate_crop_groups(config: dict, project_root: Path) -> None: """ all_crops = set(config["crops"]) # Non-food crops and grassland can also appear on production maps - all_crops.update(config.get("non_food_crops", [])) + all_crops.update(config["non_food_crops"]) all_crops.add("grassland") group_cfg = config["plotting"]["crop_groups"] diff --git a/workflow/validation/cropgrids_crops.py b/workflow/validation/cropgrids_crops.py index 432882db..ad32176d 100644 --- a/workflow/validation/cropgrids_crops.py +++ b/workflow/validation/cropgrids_crops.py @@ -35,7 +35,7 @@ def validate_cropgrids_crops(config: dict, project_root: Path) -> None: reserved for GAEZ-backed crops, and a stray entry would otherwise feed the GAEZ download rules. """ - cropgrids_crops = list(config.get("cropgrids_crops") or []) + cropgrids_crops = list(config["cropgrids_crops"]) if not cropgrids_crops: return diff --git a/workflow/validation/gaez_crop_mapping.py b/workflow/validation/gaez_crop_mapping.py index 1fd621e1..d49735b5 100644 --- a/workflow/validation/gaez_crop_mapping.py +++ b/workflow/validation/gaez_crop_mapping.py @@ -40,7 +40,7 @@ def validate_gaez_crop_mapping(config: dict, project_root: Path) -> None: # Crops sourced from CROPGRIDS bypass GAEZ entirely and have no mapping # row by design (enforced by validate_cropgrids_crops). - cropgrids_crops = set(config.get("cropgrids_crops") or []) + cropgrids_crops = set(config["cropgrids_crops"]) config_crops = set(config["crops"]) - cropgrids_crops mapped_crops = set(df["crop_name"].unique()) diff --git a/workflow/validation/health_map.py b/workflow/validation/health_map.py index eafdc3d3..0599c7d2 100644 --- a/workflow/validation/health_map.py +++ b/workflow/validation/health_map.py @@ -10,7 +10,7 @@ def validate_health_map(config: dict, project_root: Path) -> None: - health = config.get("health", {}) + health = config["health"] risks = set(health.get("risk_factors", [])) causes = set(health.get("causes", [])) risk_cause_map: dict[str, list[str]] = health.get("risk_cause_map", {}) @@ -41,7 +41,7 @@ def validate_health_map(config: dict, project_root: Path) -> None: # Every risk factor needs a per-capita consumption cap: it both bounds the # food-group store (e_nom_max) and sets the upper end of the Stage 1 intake # breakpoint domain in prepare_health_costs. - max_per_capita = set(config.get("food_groups", {}).get("max_per_capita", {})) + max_per_capita = set(config["food_groups"]["max_per_capita"]) missing_caps = risks - max_per_capita if missing_caps: raise ValueError( diff --git a/workflow/validation/optimal_taxes.py b/workflow/validation/optimal_taxes.py index f402f6e0..8f321c64 100644 --- a/workflow/validation/optimal_taxes.py +++ b/workflow/validation/optimal_taxes.py @@ -21,7 +21,7 @@ def _optimal_taxes_enabled(config: dict, scenario_defs: dict) -> bool: def validate_optimal_taxes(config: dict, _project_root=None) -> None: """Ensure optimal taxes runs have required scenarios defined.""" - scenario_defs = config.get("scenarios") or {} + scenario_defs = config["scenarios"] if not _optimal_taxes_enabled(config, scenario_defs): return diff --git a/workflow/validation/sensitivity_generator.py b/workflow/validation/sensitivity_generator.py index af00f2ae..68c273f1 100644 --- a/workflow/validation/sensitivity_generator.py +++ b/workflow/validation/sensitivity_generator.py @@ -7,7 +7,7 @@ def validate_sensitivity_generator(config: dict, _project_root=None) -> None: """Ensure sensitivity generators have unique prefixes.""" - scenario_defs = config.get("scenarios") or {} + scenario_defs = config["scenarios"] sensitivity_generators = [ generator From 4c9f0230658b46b81b83efdc61df04f001887fa7 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Tue, 28 Jul 2026 14:34:26 -0700 Subject: [PATCH 10/10] fix: track structural calibration inputs --- data/curated/calibration/default/provenance.yaml | 16 ++++++++++++++++ .../calibration/gbd-anchored/provenance.yaml | 16 ++++++++++++++++ workflow/validation/calibration_provenance.py | 3 +++ 3 files changed, 35 insertions(+) diff --git a/data/curated/calibration/default/provenance.yaml b/data/curated/calibration/default/provenance.yaml index ddba1a95..0c5fe7e5 100644 --- a/data/curated/calibration/default/provenance.yaml +++ b/data/curated/calibration/default/provenance.yaml @@ -199,6 +199,7 @@ structural_config: - rendered-fat biomass.enforce_baseline_demand: true biomass.enforce_fiber_demand: true + biomass.marginal_values_usd_per_tonne: 0 byproducts: - beet-pulp - wheat-bran @@ -902,6 +903,21 @@ structural_config: - sugar - stimulants - animal_fat + food_groups.max_per_capita.animal_fat: 50 + food_groups.max_per_capita.dairy: 2865 + food_groups.max_per_capita.eggs: 213 + food_groups.max_per_capita.fruits: 658 + food_groups.max_per_capita.grain: 1403 + food_groups.max_per_capita.legumes: 300 + food_groups.max_per_capita.nuts_seeds: 79 + food_groups.max_per_capita.oil: 155 + food_groups.max_per_capita.poultry: 241 + food_groups.max_per_capita.red_meat: 285 + food_groups.max_per_capita.starchy_vegetable: 1221 + food_groups.max_per_capita.stimulants: 50 + food_groups.max_per_capita.sugar: 133 + food_groups.max_per_capita.vegetables: 785 + food_groups.max_per_capita.whole_grains: 300 food_loss_waste_calibration.food_groups: - vegetables - fruits diff --git a/data/curated/calibration/gbd-anchored/provenance.yaml b/data/curated/calibration/gbd-anchored/provenance.yaml index fafd0e77..95459eb8 100644 --- a/data/curated/calibration/gbd-anchored/provenance.yaml +++ b/data/curated/calibration/gbd-anchored/provenance.yaml @@ -199,6 +199,7 @@ structural_config: - rendered-fat biomass.enforce_baseline_demand: true biomass.enforce_fiber_demand: true + biomass.marginal_values_usd_per_tonne: 0 byproducts: - beet-pulp - wheat-bran @@ -902,6 +903,21 @@ structural_config: - sugar - stimulants - animal_fat + food_groups.max_per_capita.animal_fat: 50 + food_groups.max_per_capita.dairy: 2865 + food_groups.max_per_capita.eggs: 213 + food_groups.max_per_capita.fruits: 658 + food_groups.max_per_capita.grain: 1403 + food_groups.max_per_capita.legumes: 300 + food_groups.max_per_capita.nuts_seeds: 79 + food_groups.max_per_capita.oil: 155 + food_groups.max_per_capita.poultry: 241 + food_groups.max_per_capita.red_meat: 285 + food_groups.max_per_capita.starchy_vegetable: 1221 + food_groups.max_per_capita.stimulants: 50 + food_groups.max_per_capita.sugar: 133 + food_groups.max_per_capita.vegetables: 785 + food_groups.max_per_capita.whole_grains: 300 food_loss_waste_calibration.food_groups: - vegetables - fruits diff --git a/workflow/validation/calibration_provenance.py b/workflow/validation/calibration_provenance.py index 8f40872a..6b4ad4e5 100644 --- a/workflow/validation/calibration_provenance.py +++ b/workflow/validation/calibration_provenance.py @@ -47,6 +47,8 @@ "downloads", "credentials", "calibration", + "remote_solve", + "solving.inline_analysis", # Validation-mode switches: the artefacts are consumed by regular and # validation-mode solves alike (the calibration chain itself runs in # validation mode). @@ -55,6 +57,7 @@ # horizon by design. "planning_horizon", # Post-solve analysis only. + "plotting", "sensitivity_analysis", # Calibration application/generation machinery. Fit-relevant knobs in # these sections (e.g. food_loss_waste_calibration.food_groups,