From f3a011ce8efb1f1e76c7ec8ee6e84f2b817120ce Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Thu, 4 Jun 2026 22:35:15 -0400 Subject: [PATCH 1/9] Add new function to scale V using authoritative x --- .../test_scale_V_from_authoritative_x.py | 30 ++++++++++ bedrock/transform/eeio/derived_cornerstone.py | 58 +++++++++++++++++-- 2 files changed, 82 insertions(+), 6 deletions(-) create mode 100644 bedrock/transform/eeio/__tests__/test_scale_V_from_authoritative_x.py diff --git a/bedrock/transform/eeio/__tests__/test_scale_V_from_authoritative_x.py b/bedrock/transform/eeio/__tests__/test_scale_V_from_authoritative_x.py new file mode 100644 index 00000000..81b910fc --- /dev/null +++ b/bedrock/transform/eeio/__tests__/test_scale_V_from_authoritative_x.py @@ -0,0 +1,30 @@ +from __future__ import annotations + +import numpy as np + +from bedrock.transform.eeio.derived_cornerstone import ( + derive_cornerstone_x_after_redefinition, + scale_cornerstone_V_with_authoritative_x, +) +from bedrock.utils.config.usa_config import get_usa_config + + +def test_row_sums_equal_authoritative_x() -> None: + """V_new row sums must match the authoritative model-year x for all non-zero industries.""" + scale_cornerstone_V_with_authoritative_x.cache_clear() + derive_cornerstone_x_after_redefinition.cache_clear() + + cfg = get_usa_config() + V_new = scale_cornerstone_V_with_authoritative_x() + x_new = derive_cornerstone_x_after_redefinition(year=cfg.model_base_year) + + row_sums = V_new.sum(axis=1) + x_aligned = x_new.reindex(row_sums.index).fillna(0.0) + + mask = x_aligned.values != 0 + np.testing.assert_allclose( + row_sums.values[mask], + x_aligned.values[mask], + rtol=1e-6, + err_msg="V_new row sums do not match authoritative x_new", + ) diff --git a/bedrock/transform/eeio/derived_cornerstone.py b/bedrock/transform/eeio/derived_cornerstone.py index d082134d..08cf0a36 100644 --- a/bedrock/transform/eeio/derived_cornerstone.py +++ b/bedrock/transform/eeio/derived_cornerstone.py @@ -326,26 +326,28 @@ def _distribute_waste_parent_x_using_v_row_shares( @functools.cache @pa.check_output(CornerstoneXVectorSchema) -def derive_cornerstone_x_after_redefinition() -> pd.Series[float]: +def derive_cornerstone_x_after_redefinition(year: int = 0) -> pd.Series[float]: """Gross industry output in Cornerstone schema, after BEA redefinitions. - Uses gross-output time series for the configured GHG data year, selecting - before/after-redefinition source from config, then expands it to - Cornerstone industries via the BEA→Cornerstone industry correspondence. + Uses gross-output time series for *year* (defaults to + ``usa_ghg_data_year`` when *year* is 0), selecting before/after-redefinition + source from config, then expands it to Cornerstone industries via the + BEA→Cornerstone industry correspondence. For one-to-many splits (e.g. waste 562000), ``expand_vector`` first duplicates the parent scalar to each child. When waste disaggregation is on, those waste rows are then replaced so each child gets a share of the parent total consistent with row sums of disaggregated ``V`` (same nominal - level as the GHG-year BEA gross output, split from 2017 Make structure). + level as the BEA gross output for *year*, split from 2017 Make structure). This is the industry ``x`` in ``derive_cornerstone_B_via_vnorm`` when ``use_E_data_year_for_x_in_B`` is True; otherwise that path uses ``derive_cornerstone_x()``. """ cfg = get_usa_config() + effective_year = cfg.usa_ghg_data_year if year == 0 else year x_bea = derive_gross_output( - target_year=cfg.usa_ghg_data_year, + target_year=effective_year, iot_before_or_after_redefinition=cfg.iot_before_or_after_redefinition, ) x_cs = expand_vector(x_bea, CS_INDUSTRY_LIST, cs_industry_to_bea_map()) @@ -398,6 +400,50 @@ def derive_cornerstone_Vnorm_scrap_corrected( return V_scrap_corrected +@functools.cache +def scale_cornerstone_V_with_authoritative_x() -> pd.DataFrame: + """Estimate V rescaled to match model-year gross industry output. + + Derives x_new from the BEA gross-output time series at ``model_base_year`` + (after redefinitions) via ``derive_cornerstone_x_after_redefinition``. + Inflates the base 2017 V to model-year dollars, computes its row sums + (x_model_year), then scales each row i proportionally by + ``x_new[i] / x_model_year[i]``. Industries with zero model-year output + receive a zero row. + + Returns + ------- + pd.DataFrame + New V with ``V_new.sum(axis=1) ≈ x_new`` for all industries that have + non-zero model-year output. + """ + cfg = get_usa_config() + x_new = derive_cornerstone_x_after_redefinition(year=cfg.model_base_year) + + V_model_year = derive_cornerstone_V( + apply_inflation=True, target_year=cfg.model_base_year + ) + x_model_year = V_model_year.sum(axis=1) + x_new_aligned = x_new.reindex(x_model_year.index).fillna(0.0) + + scale = pd.Series( + np.where( + x_model_year.values != 0, + x_new_aligned.values / x_model_year.values, + 0.0, + ), + index=x_model_year.index, + ) + V_new = V_model_year.multiply(scale, axis=0) + + mask = x_model_year.values != 0 + assert np.allclose( + V_new.sum(axis=1).values[mask], x_new_aligned.values[mask], rtol=1e-6 + ), "Row sums of V_new do not match x_new" + + return V_new + + # --------------------------------------------------------------------------- # Base 2017 IO matrices — U # --------------------------------------------------------------------------- From 1cc599d7b8b59e51bd51c4429c19d78207ede71d Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Thu, 4 Jun 2026 22:44:11 -0400 Subject: [PATCH 2/9] add derive_q_from_scaled_cornerstone_V_from_authoritative_x --- bedrock/transform/eeio/derived_cornerstone.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/bedrock/transform/eeio/derived_cornerstone.py b/bedrock/transform/eeio/derived_cornerstone.py index 08cf0a36..5c44bcfa 100644 --- a/bedrock/transform/eeio/derived_cornerstone.py +++ b/bedrock/transform/eeio/derived_cornerstone.py @@ -138,7 +138,7 @@ _BEDROCK_PKG_ROOT = pathlib.Path(__file__).resolve().parents[2] -_WASTE_ORIGINAL_CODE = "562000" +_WASTE_ORIGINAL_CODE = '562000' _WASTE_NEW_CODES: list[str] = list(WASTE_DISAGG_COMMODITIES[_WASTE_ORIGINAL_CODE]) # --------------------------------------------------------------------------- @@ -351,7 +351,7 @@ def derive_cornerstone_x_after_redefinition(year: int = 0) -> pd.Series[float]: iot_before_or_after_redefinition=cfg.iot_before_or_after_redefinition, ) x_cs = expand_vector(x_bea, CS_INDUSTRY_LIST, cs_industry_to_bea_map()) - x_cs.index.name = "sector" + x_cs.index.name = 'sector' return _distribute_waste_parent_x_using_v_row_shares(x_cs) @@ -439,11 +439,16 @@ def scale_cornerstone_V_with_authoritative_x() -> pd.DataFrame: mask = x_model_year.values != 0 assert np.allclose( V_new.sum(axis=1).values[mask], x_new_aligned.values[mask], rtol=1e-6 - ), "Row sums of V_new do not match x_new" + ), 'Row sums of V_new do not match x_new' return V_new +def derive_q_from_scaled_cornerstone_V_from_authoritative_x() -> pd.Series[float]: + V = scale_cornerstone_V_with_authoritative_x() + return compute_q(V=V) + + # --------------------------------------------------------------------------- # Base 2017 IO matrices — U # --------------------------------------------------------------------------- From 6ea638b726dc859bd2d566681a48a5a9e066fdc3 Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Fri, 5 Jun 2026 12:53:38 -0400 Subject: [PATCH 3/9] Add compare_q script --- bedrock/analysis/q/compare_alt_q.py | 525 ++++++++++++++++++++++++++++ 1 file changed, 525 insertions(+) create mode 100644 bedrock/analysis/q/compare_alt_q.py diff --git a/bedrock/analysis/q/compare_alt_q.py b/bedrock/analysis/q/compare_alt_q.py new file mode 100644 index 00000000..5913dd12 --- /dev/null +++ b/bedrock/analysis/q/compare_alt_q.py @@ -0,0 +1,525 @@ +"""Compare three q-estimation pathways at Cornerstone detail and BEA summary level. + +Pathways +-------- +1. q_V_scaled + derive_q_from_scaled_cornerstone_V_from_authoritative_x() + V is inflated to model_base_year then rescaled so each industry row sum equals + the authoritative BEA gross-output x. q = V.sum(axis=0). + +2. q_Aq_scaled + scale_cornerstone_q (with dollar-year adjustment) + inflate_with_commodity_pi. + Replicates the derive_cornerstone_Aq_scaled() branch that is active when + scale_a_matrix_with_summary_tables=True AND adjust_summary_A_and_q_dollar_year=True. + The dollar-year flag is hardwired here so the comparison is independent of the + live config. + +3. q_inflated + inflate_cornerstone_q_or_y_with_commodity_pi applied directly to the base 2017 q + (no structural summary-table scaling). This is the commodity-price-index-only path. + +Reference vector +---------------- +x_authoritative + derive_cornerstone_x_after_redefinition(year=model_base_year): BEA annual gross + industry output at model_base_year, after redefinitions, expanded to the 405-sector + Cornerstone schema. Already in model_base_year current dollars — no additional + inflation applied. Industry codes and commodity codes share the same 405-sector + namespace, so each q pathway is compared directly to this x. + +At the summary level each pathway is summed via the Cornerstone→BEA-summary mapping +and compared both to derive_summary_q_usa(model_base_year) (the "Make Table" q) and +to the aggregated x_authoritative. + +Usage +----- + python -m bedrock.analysis.compare_alt_q +""" + +from __future__ import annotations + +import logging +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from bedrock.transform.eeio.derived_2017 import derive_summary_q_usa +from bedrock.transform.eeio.derived_cornerstone import ( + derive_cornerstone_Aq, + derive_cornerstone_V, + derive_cornerstone_x, + derive_cornerstone_x_after_redefinition, + derive_q_from_scaled_cornerstone_V_from_authoritative_x, +) +from bedrock.utils.config.usa_config import get_usa_config +from bedrock.utils.economic.inflation_helpers_cornerstone import ( + adjust_summary_q_dollar_year, + inflate_cornerstone_q_or_y_with_commodity_pi, +) +from bedrock.utils.math.formulas import compute_q +from bedrock.utils.taxonomy.bea_v2017_to_ceda_v7_helpers import ( + load_bea_v2017_summary_to_cornerstone, +) + +logger = logging.getLogger(__name__) + +OUTPUT_DIR = Path(__file__).parent / "output" +PLOTS_DIR = OUTPUT_DIR / "plots" + +_PATHWAY_STYLE: list[tuple[str, str, str]] = [ + ("q_V_scaled", "V_scaled", "#1f77b4"), + ("q_Aq_scaled", "Aq_scaled", "#ff7f0e"), + ("q_inflated", "inflated", "#2ca02c"), +] + + +# --------------------------------------------------------------------------- +# Pathway builders +# --------------------------------------------------------------------------- + + +def _build_q_V_scaled() -> pd.Series: + """Pathway 1: q from V rescaled to authoritative gross-output x.""" + return derive_q_from_scaled_cornerstone_V_from_authoritative_x() + + +def _build_q_Aq_scaled( + base_q: pd.Series, + detail_year: int, + model_year: int, +) -> pd.Series: + """Pathway 2: scale (dollar-year-adjusted) + inflate with commodity PI. + + Replicates derive_cornerstone_Aq_scaled() with + scale_a_matrix_with_summary_tables=True and + adjust_summary_A_and_q_dollar_year=True, regardless of the active config. + """ + # --- scale step (dollar-year-adjusted) --- + q_summary_target = derive_summary_q_usa(model_year) + q_summary_target_adj = adjust_summary_q_dollar_year( + q_summary=q_summary_target, + from_year=model_year, + to_year=detail_year, + ) + q_summary_base = derive_summary_q_usa(detail_year) + ratio = (q_summary_target_adj / q_summary_base).fillna(1.0) + ratio = ratio.replace([np.inf, -np.inf], 1.0) + + summary_to_cornerstone = load_bea_v2017_summary_to_cornerstone() + q_scaled = base_q.copy() + for summary_sector, val in ratio.items(): + if summary_sector not in summary_to_cornerstone: + continue + sectors = [ + s for s in summary_to_cornerstone[summary_sector] if s in q_scaled.index + ] + if sectors: + q_scaled.loc[sectors] *= float(val) + + # --- inflate step --- + return inflate_cornerstone_q_or_y_with_commodity_pi( + q_scaled, original_year=detail_year, target_year=model_year + ) + + +def _build_q_inflated( + base_q: pd.Series, + detail_year: int, + model_year: int, +) -> pd.Series: + """Pathway 3: inflate base 2017 q directly with commodity PI (no scaling).""" + return inflate_cornerstone_q_or_y_with_commodity_pi( + base_q, original_year=detail_year, target_year=model_year + ) + + +# --------------------------------------------------------------------------- +# Aggregation to summary level +# --------------------------------------------------------------------------- + + +def _cornerstone_to_summary_map() -> dict[str, str]: + """cornerstone_code → bea_summary_code lookup.""" + summary_to_cornerstone = load_bea_v2017_summary_to_cornerstone() + return { + code: str(summary) + for summary, codes in summary_to_cornerstone.items() + for code in codes + } + + +def aggregate_q_to_summary(q: pd.Series, cs_to_summary: dict[str, str]) -> pd.Series: + """Sum q within each BEA summary parent. + + Any sector in q that is absent from cs_to_summary is silently dropped + by the valid_summaries filter. A warning is emitted with the count and + dollar sum so that unexpected losses surface immediately. + """ + unmapped = [c for c in q.index if c not in cs_to_summary] + if unmapped: + lost = q.loc[unmapped].sum() + logger.warning( + "aggregate_q_to_summary: %d sector(s) not in cs_to_summary " + "and will be dropped (total $M %.0f): %s", + len(unmapped), + lost, + unmapped, + ) + mapped = q.rename(index=cs_to_summary) + agg = mapped.groupby(level=0).sum() + valid_summaries = set(cs_to_summary.values()) + agg = agg.loc[[c for c in agg.index if c in valid_summaries]] + return agg + + +# --------------------------------------------------------------------------- +# Comparison tables +# --------------------------------------------------------------------------- + + +def _pct_diff(a: pd.Series, b: pd.Series, label_a: str, label_b: str) -> pd.Series: + """(a − b) / b × 100, aligned on common index.""" + common = a.index.intersection(b.index) + denom = b.reindex(common).replace(0, np.nan) + return ((a.reindex(common) - denom) / denom * 100).rename( + f"pct_diff_{label_a}_vs_{label_b}" + ) + + +def _build_detail_comparison( + q_V_scaled: pd.Series, + q_Aq_scaled: pd.Series, + q_inflated: pd.Series, + x_authoritative: pd.Series, +) -> pd.DataFrame: + """Side-by-side detail table with pairwise pct differences. + + x_authoritative is industry gross output; q series are commodity gross + output. Both share the 405-sector Cornerstone code space so direct + element-wise comparison is valid. + """ + common = ( + q_V_scaled.index.intersection(q_Aq_scaled.index) + .intersection(q_inflated.index) + .intersection(x_authoritative.index) + ) + df = pd.DataFrame( + { + "x_authoritative": x_authoritative.reindex(common), + "q_V_scaled": q_V_scaled.reindex(common), + "q_Aq_scaled": q_Aq_scaled.reindex(common), + "q_inflated": q_inflated.reindex(common), + } + ) + df["pct_Aq_vs_inflated"] = _pct_diff( + q_Aq_scaled, q_inflated, "Aq_scaled", "inflated" + ).reindex(common) + df["pct_V_vs_Aq"] = _pct_diff( + q_V_scaled, q_Aq_scaled, "V_scaled", "Aq_scaled" + ).reindex(common) + df["pct_V_vs_inflated"] = _pct_diff( + q_V_scaled, q_inflated, "V_scaled", "inflated" + ).reindex(common) + df["pct_V_vs_x"] = _pct_diff(q_V_scaled, x_authoritative, "V_scaled", "x").reindex( + common + ) + df["pct_Aq_vs_x"] = _pct_diff( + q_Aq_scaled, x_authoritative, "Aq_scaled", "x" + ).reindex(common) + df["pct_inflated_vs_x"] = _pct_diff( + q_inflated, x_authoritative, "inflated", "x" + ).reindex(common) + return df.sort_values("x_authoritative", ascending=False) + + +def _build_summary_comparison( + q_V_scaled: pd.Series, + q_Aq_scaled: pd.Series, + q_inflated: pd.Series, + x_authoritative: pd.Series, + q_make_table: pd.Series, + cs_to_summary: dict[str, str], +) -> pd.DataFrame: + """Summary-level table: each pathway vs Make Table q and vs x_authoritative.""" + agg_V = aggregate_q_to_summary(q_V_scaled, cs_to_summary) + agg_Aq = aggregate_q_to_summary(q_Aq_scaled, cs_to_summary) + agg_inf = aggregate_q_to_summary(q_inflated, cs_to_summary) + agg_x = aggregate_q_to_summary(x_authoritative, cs_to_summary) + + common = ( + agg_V.index.intersection(agg_Aq.index) + .intersection(agg_inf.index) + .intersection(agg_x.index) + .intersection(q_make_table.index) + ) + df = pd.DataFrame( + { + "x_authoritative_agg": agg_x.reindex(common), + "q_make_table": q_make_table.reindex(common), + "q_V_scaled_agg": agg_V.reindex(common), + "q_Aq_scaled_agg": agg_Aq.reindex(common), + "q_inflated_agg": agg_inf.reindex(common), + } + ) + denom_mt = df["q_make_table"].replace(0, np.nan) + denom_x = df["x_authoritative_agg"].replace(0, np.nan) + for col, label in [ + ("q_V_scaled_agg", "V_scaled"), + ("q_Aq_scaled_agg", "Aq_scaled"), + ("q_inflated_agg", "inflated"), + ]: + df[f"pct_{label}_vs_make_table"] = (df[col] - denom_mt) / denom_mt * 100 + df[f"pct_{label}_vs_x"] = (df[col] - denom_x) / denom_x * 100 + + return df.sort_values("x_authoritative_agg", ascending=False) + + +# --------------------------------------------------------------------------- +# Plots +# --------------------------------------------------------------------------- + + +def _scatter_panel( + ax: plt.Axes, + x: np.ndarray, + y: np.ndarray, + title: str, + xlabel: str, + ylabel: str, + color: str, +) -> None: + """Shared helper: log-log scatter with identity line on a single Axes.""" + mask = (x > 0) & (y > 0) + ax.scatter(x[mask], y[mask], s=5, alpha=0.5, color=color, zorder=2) + lo = min(x[mask].min(), y[mask].min()) * 0.9 + hi = max(x[mask].max(), y[mask].max()) * 1.1 + ax.plot([lo, hi], [lo, hi], "k--", lw=0.8, zorder=3) + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_title(title, fontsize=10) + ax.set_xlabel(xlabel, fontsize=8) + ax.set_ylabel(ylabel, fontsize=8) + ax.grid(True, alpha=0.25, which="both") + + +def _plot_detail_scatter( + detail_df: pd.DataFrame, + q_2017: pd.Series, + x_2017: pd.Series, + model_year: int, +) -> Path: + """1×4 log-log scatter at detail level. + + Panels 1–3: each q pathway vs x_authoritative (model_year). + Panel 4 (reference): 2017 q vs 2017 x, both from the Cornerstone Make + matrix after redefinitions — shows the baseline q/x relationship before + any scaling or inflation is applied. + """ + fig, axes = plt.subplots(1, 4, figsize=(20, 5)) + fig.suptitle( + f"Detail-level q vs x — model_year {model_year} | panel 4: 2017 reference", + fontsize=11, + ) + + x_auth = detail_df["x_authoritative"].to_numpy(dtype=float) + for ax, (col, label, color) in zip(axes[:3], _PATHWAY_STYLE): + _scatter_panel( + ax, + x=x_auth, + y=detail_df[col].to_numpy(dtype=float), + title=label, + xlabel="x_authoritative ($M)", + ylabel=f"q {label} ($M)", + color=color, + ) + + # Reference panel: 2017 q vs 2017 x + common_2017 = q_2017.index.intersection(x_2017.index) + _scatter_panel( + axes[3], + x=x_2017.reindex(common_2017).to_numpy(dtype=float), + y=q_2017.reindex(common_2017).to_numpy(dtype=float), + title="2017 reference (q vs x, same V)", + xlabel="x 2017 Cornerstone Make ($M)", + ylabel="q 2017 Cornerstone Make ($M)", + color="#9467bd", + ) + + fig.tight_layout() + path = PLOTS_DIR / f"detail_scatter_{model_year}.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +def _plot_summary_levels(summary_df: pd.DataFrame, model_year: int) -> Path: + """Ranked line plot of summary-level absolute q values (log y-axis). + + Sectors sorted by x_authoritative_agg descending (left = largest). + """ + df = summary_df.sort_values("x_authoritative_agg", ascending=False) + ranks = np.arange(len(df)) + + series = [ + ("x_authoritative_agg", "x_authoritative", "black", "o", 2.0), + ("q_make_table", "q_make_table", "#7f7f7f", "s", 1.5), + ] + [(f"{col}_agg", label, color, "^", 1.0) for col, label, color in _PATHWAY_STYLE] + + fig, ax = plt.subplots(figsize=(14, 5)) + for colname, label, color, marker, lw in series: + vals = df[colname].to_numpy(dtype=float) + ax.plot( + ranks, + vals, + marker=marker, + markersize=3, + linewidth=lw, + color=color, + label=label, + alpha=0.85, + ) + + ax.set_yscale("log") + ax.set_xticks(ranks[::4]) + ax.set_xticklabels(df.index[::4], rotation=45, ha="right", fontsize=7) + ax.set_xlabel("BEA summary sector (ranked by x_authoritative)", fontsize=9) + ax.set_ylabel("$M (log scale)", fontsize=9) + ax.set_title(f"Summary-level q pathways vs references — {model_year}", fontsize=11) + ax.legend(fontsize=8, frameon=False) + ax.grid(True, axis="y", alpha=0.3) + fig.tight_layout() + path = PLOTS_DIR / f"summary_levels_{model_year}.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +def _plot_summary_pct_deviations(summary_df: pd.DataFrame, model_year: int) -> Path: + """1×2 grouped horizontal bar chart: % deviation from x_authoritative (left) + and from Make Table q (right), sectors sorted by x_authoritative_agg descending. + """ + df = summary_df.sort_values( + "x_authoritative_agg", ascending=True + ) # ascending → biggest at top + sectors = df.index.tolist() + n = len(sectors) + y = np.arange(n) + bar_h = 0.25 + + fig, axes = plt.subplots(1, 2, figsize=(16, max(6, n * 0.22)), sharey=True) + fig.suptitle(f"Summary-level q % deviation — model_year {model_year}", fontsize=11) + + for ax, ref_suffix, ref_label in [ + (axes[0], "x", "x_authoritative"), + (axes[1], "make_table", "Make Table q"), + ]: + for i, (_, label, color) in enumerate(_PATHWAY_STYLE): + pct_col = f"pct_{label}_vs_{ref_suffix}" + vals = df[pct_col].to_numpy(dtype=float) + ax.barh( + y + (i - 1) * bar_h, + vals, + height=bar_h, + color=color, + label=label, + alpha=0.8, + ) + ax.axvline(0, color="black", linewidth=0.8) + ax.set_xlabel("% deviation", fontsize=9) + ax.set_title(f"vs {ref_label}", fontsize=10) + ax.grid(True, axis="x", alpha=0.3) + ax.legend(fontsize=8, frameon=False) + + axes[0].set_yticks(y) + axes[0].set_yticklabels(sectors, fontsize=7) + fig.tight_layout() + path = PLOTS_DIR / f"summary_pct_deviations_{model_year}.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + + +def main() -> None: + cfg = get_usa_config() + detail_year: int = cfg.usa_detail_original_year + model_year: int = cfg.model_base_year + + print(f"compare_alt_q detail_year={detail_year} model_year={model_year}") + + base_q = derive_cornerstone_Aq().scaled_q + x_2017 = derive_cornerstone_x() + q_2017 = compute_q(V=derive_cornerstone_V()) + q_V_scaled = _build_q_V_scaled() + q_Aq_scaled = _build_q_Aq_scaled(base_q, detail_year, model_year) + q_inflated = _build_q_inflated(base_q, detail_year, model_year) + x_authoritative = derive_cornerstone_x_after_redefinition(year=model_year) + q_make_table = derive_summary_q_usa(model_year) + + detail_df = _build_detail_comparison( + q_V_scaled, q_Aq_scaled, q_inflated, x_authoritative + ) + cs_to_summary = _cornerstone_to_summary_map() + summary_df = _build_summary_comparison( + q_V_scaled, + q_Aq_scaled, + q_inflated, + x_authoritative, + q_make_table, + cs_to_summary, + ) + + # ------------------------------------------------------------------ + # Printed summary stats + # ------------------------------------------------------------------ + print("\nAggregate totals at summary level ($M):") + for col, label in [ + ("x_authoritative_agg", "x_authoritative"), + ("q_make_table", "q_make_table"), + ("q_V_scaled_agg", "V_scaled"), + ("q_Aq_scaled_agg", "Aq_scaled"), + ("q_inflated_agg", "inflated"), + ]: + print(f" {label:<20} {summary_df[col].sum():>14,.0f}") + + print("\nOverall pct deviation from Make Table q:") + denom_mt = summary_df["q_make_table"].sum() + for col, label, _ in _PATHWAY_STYLE: + agg_col = f"{col}_agg" + print( + f" {label:<18} {(summary_df[agg_col].sum() - denom_mt) / denom_mt * 100:+.3f}%" + ) + + print("\nOverall pct deviation from x_authoritative:") + denom_x = summary_df["x_authoritative_agg"].sum() + for col, label, _ in _PATHWAY_STYLE: + agg_col = f"{col}_agg" + print( + f" {label:<18} {(summary_df[agg_col].sum() - denom_x) / denom_x * 100:+.3f}%" + ) + + # ------------------------------------------------------------------ + # Plots + # ------------------------------------------------------------------ + PLOTS_DIR.mkdir(parents=True, exist_ok=True) + p1 = _plot_detail_scatter( + detail_df, q_2017=q_2017, x_2017=x_2017, model_year=model_year + ) + p2 = _plot_summary_levels(summary_df, model_year) + p3 = _plot_summary_pct_deviations(summary_df, model_year) + print(f"\nPlots saved to {PLOTS_DIR}:") + for p in (p1, p2, p3): + print(f" {p.name}") + + +if __name__ == "__main__": + logging.basicConfig( + level=logging.INFO, format="%(levelname)s %(name)s: %(message)s" + ) + main() From b38c5ee853aa3513456a7e8e2043bdac0f95560b Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Fri, 5 Jun 2026 13:10:07 -0400 Subject: [PATCH 4/9] Fix type errors --- bedrock/analysis/q/compare_alt_q.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/bedrock/analysis/q/compare_alt_q.py b/bedrock/analysis/q/compare_alt_q.py index 5913dd12..d17d0686 100644 --- a/bedrock/analysis/q/compare_alt_q.py +++ b/bedrock/analysis/q/compare_alt_q.py @@ -109,12 +109,11 @@ def _build_q_Aq_scaled( summary_to_cornerstone = load_bea_v2017_summary_to_cornerstone() q_scaled = base_q.copy() - for summary_sector, val in ratio.items(): - if summary_sector not in summary_to_cornerstone: + for summary_sector, cs_sectors in summary_to_cornerstone.items(): + if summary_sector not in ratio.index: continue - sectors = [ - s for s in summary_to_cornerstone[summary_sector] if s in q_scaled.index - ] + val = ratio.loc[summary_sector] + sectors = [s for s in cs_sectors if s in q_scaled.index] if sectors: q_scaled.loc[sectors] *= float(val) From 54cfbb8675549f54d3551f5c97e04332f7696b62 Mon Sep 17 00:00:00 2001 From: "Dr. Wesley Ingwersen" Date: Fri, 5 Jun 2026 13:32:51 -0400 Subject: [PATCH 5/9] Update bedrock/analysis/q/compare_alt_q.py Co-authored-by: Ben Young (ERG) <44471635+bl-young@users.noreply.github.com> --- bedrock/analysis/q/compare_alt_q.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bedrock/analysis/q/compare_alt_q.py b/bedrock/analysis/q/compare_alt_q.py index d17d0686..16e3e285 100644 --- a/bedrock/analysis/q/compare_alt_q.py +++ b/bedrock/analysis/q/compare_alt_q.py @@ -33,7 +33,7 @@ Usage ----- - python -m bedrock.analysis.compare_alt_q + python -m bedrock.analysis.q.compare_alt_q """ from __future__ import annotations From 610b8845830e5b31758924ffefc92795947fdfc2 Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Fri, 5 Jun 2026 14:59:22 -0400 Subject: [PATCH 6/9] fixes mapping issue where output was being lost --- bedrock/analysis/q/compare_alt_q.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/bedrock/analysis/q/compare_alt_q.py b/bedrock/analysis/q/compare_alt_q.py index d17d0686..2f00619b 100644 --- a/bedrock/analysis/q/compare_alt_q.py +++ b/bedrock/analysis/q/compare_alt_q.py @@ -142,11 +142,16 @@ def _build_q_inflated( def _cornerstone_to_summary_map() -> dict[str, str]: """cornerstone_code → bea_summary_code lookup.""" summary_to_cornerstone = load_bea_v2017_summary_to_cornerstone() - return { + mapping = { code: str(summary) for summary, codes in summary_to_cornerstone.items() for code in codes } + # BEA industry 331314 (secondary aluminum smelting) is not a commodity code; + # after redefinitions its output is attributed to commodity 331313. + if "331313" in mapping: + mapping["331314"] = mapping["331313"] + return mapping def aggregate_q_to_summary(q: pd.Series, cs_to_summary: dict[str, str]) -> pd.Series: From cb9c990929f813d65788a429b6bcba451dd64792 Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Fri, 5 Jun 2026 15:08:57 -0400 Subject: [PATCH 7/9] fix typing errors in derived_cornerstone --- bedrock/transform/eeio/derived_cornerstone.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/bedrock/transform/eeio/derived_cornerstone.py b/bedrock/transform/eeio/derived_cornerstone.py index 5c44bcfa..2c9138ad 100644 --- a/bedrock/transform/eeio/derived_cornerstone.py +++ b/bedrock/transform/eeio/derived_cornerstone.py @@ -27,6 +27,7 @@ import functools import pathlib from dataclasses import dataclass +from typing import cast import numpy as np import pandas as pd @@ -122,6 +123,7 @@ SingleRegionYtotAndTradeVectorSet, SingleRegionYVectorSet, ) +from bedrock.utils.taxonomy.bea.matrix_mappings import USA_GROSS_INDUSTRY_OUTPUT_YEARS from bedrock.utils.taxonomy.bea.v2017_final_demand import ( USA_2017_FINAL_DEMAND_EXPORT_CODE, USA_2017_FINAL_DEMAND_IMPORT_CODE, @@ -345,7 +347,7 @@ def derive_cornerstone_x_after_redefinition(year: int = 0) -> pd.Series[float]: ``derive_cornerstone_x()``. """ cfg = get_usa_config() - effective_year = cfg.usa_ghg_data_year if year == 0 else year + effective_year = cfg.usa_ghg_data_year if year == 0 else cast("USA_GROSS_INDUSTRY_OUTPUT_YEARS", year) x_bea = derive_gross_output( target_year=effective_year, iot_before_or_after_redefinition=cfg.iot_before_or_after_redefinition, @@ -426,19 +428,21 @@ def scale_cornerstone_V_with_authoritative_x() -> pd.DataFrame: x_model_year = V_model_year.sum(axis=1) x_new_aligned = x_new.reindex(x_model_year.index).fillna(0.0) + x_model_year_np = x_model_year.to_numpy(dtype=float) + x_new_aligned_np = x_new_aligned.to_numpy(dtype=float) scale = pd.Series( np.where( - x_model_year.values != 0, - x_new_aligned.values / x_model_year.values, + x_model_year_np != 0, + x_new_aligned_np / x_model_year_np, 0.0, ), index=x_model_year.index, ) V_new = V_model_year.multiply(scale, axis=0) - mask = x_model_year.values != 0 + mask = x_model_year_np != 0 assert np.allclose( - V_new.sum(axis=1).values[mask], x_new_aligned.values[mask], rtol=1e-6 + V_new.sum(axis=1).to_numpy(dtype=float)[mask], x_new_aligned_np[mask], rtol=1e-6 ), 'Row sums of V_new do not match x_new' return V_new From c83024716f2602f34b4554f05bb7eb78f9389f43 Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Fri, 5 Jun 2026 15:11:57 -0400 Subject: [PATCH 8/9] black and ruff derived_cornerstone.py --- bedrock/transform/eeio/derived_cornerstone.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/bedrock/transform/eeio/derived_cornerstone.py b/bedrock/transform/eeio/derived_cornerstone.py index 2c9138ad..f84732b0 100644 --- a/bedrock/transform/eeio/derived_cornerstone.py +++ b/bedrock/transform/eeio/derived_cornerstone.py @@ -347,7 +347,11 @@ def derive_cornerstone_x_after_redefinition(year: int = 0) -> pd.Series[float]: ``derive_cornerstone_x()``. """ cfg = get_usa_config() - effective_year = cfg.usa_ghg_data_year if year == 0 else cast("USA_GROSS_INDUSTRY_OUTPUT_YEARS", year) + effective_year = ( + cfg.usa_ghg_data_year + if year == 0 + else cast('USA_GROSS_INDUSTRY_OUTPUT_YEARS', year) + ) x_bea = derive_gross_output( target_year=effective_year, iot_before_or_after_redefinition=cfg.iot_before_or_after_redefinition, From 7a0365e3cbc7e5624668b9c2958781f7ff0af40c Mon Sep 17 00:00:00 2001 From: WesIngwersen Date: Mon, 15 Jun 2026 15:58:29 -0400 Subject: [PATCH 9/9] fix type errors --- bedrock/transform/eeio/derived_cornerstone.py | 1 + 1 file changed, 1 insertion(+) diff --git a/bedrock/transform/eeio/derived_cornerstone.py b/bedrock/transform/eeio/derived_cornerstone.py index 3c793114..8c8aaeca 100644 --- a/bedrock/transform/eeio/derived_cornerstone.py +++ b/bedrock/transform/eeio/derived_cornerstone.py @@ -26,6 +26,7 @@ from __future__ import annotations import functools +from typing import cast import numpy as np import pandas as pd