Model/scripts/decompose_co2_pe_error.py

"""Decompose the API-path CO2/PE over-estimate against lodged EPC figures.

The corpus integration test surfaced a systematic +5-10% over-estimate on
CO2 and PE while cost/SAP is well-calibrated. This script profiles the
structure of that bias: multiplicative vs additive, per-component shares,
and segmentation by fuel / PV / heating type — to localise the cause.
"""
from __future__ import annotations

import json
import statistics as stats
from collections import defaultdict
from pathlib import Path
from typing import Any

from datatypes.epc.domain.mapper import EpcPropertyDataMapper
from domain.sap10_calculator.calculator import calculate_sap_from_inputs
from domain.sap10_calculator.rdsap.cert_to_inputs import (
    SAP_10_2_SPEC_PRICES,
    cert_to_inputs,
)

CORPUS = Path("backend/epc_api/json_samples/RdSAP-Schema-21.0.1/corpus.jsonl")


def main() -> None:
    docs = [
        json.loads(line)
        for line in CORPUS.read_text().splitlines()
        if line.strip()
    ]
    rows: list[dict[str, Any]] = []
    for doc in docs:
        lodged_sap = doc.get("energy_rating_current")
        lodged_co2 = doc.get("co2_emissions_current")  # t/yr
        lodged_pe = doc.get("energy_consumption_current")  # kWh/m2/yr
        if lodged_pe is None or lodged_co2 is None:
            continue
        try:
            epc = EpcPropertyDataMapper.from_api_response(doc)
            res = calculate_sap_from_inputs(
                cert_to_inputs(epc, prices=SAP_10_2_SPEC_PRICES)
            )
        except Exception:
            continue
        tfa = res.intermediate["tfa_m2"]
        im = res.intermediate
        rows.append({
            "cert": doc.get("rrn") or "",
            "sap_err": res.sap_score_continuous - (lodged_sap or 0),
            "our_pe": res.primary_energy_kwh_per_m2,
            "lodged_pe": lodged_pe,
            "pe_diff": res.primary_energy_kwh_per_m2 - lodged_pe,
            "pe_ratio": res.primary_energy_kwh_per_m2 / lodged_pe if lodged_pe else 0,
            "our_co2": res.co2_kg_per_yr / 1000.0,
            "lodged_co2": lodged_co2,
            "co2_diff": res.co2_kg_per_yr / 1000.0 - lodged_co2,
            "co2_ratio": (res.co2_kg_per_yr / 1000.0) / lodged_co2 if lodged_co2 else 0,
            # PE components per m2
            "space_pe": im["space_heating_pe_kwh_per_m2"],
            "hw_pe": im["hot_water_pe_kwh_per_m2"],
            "other_pe": im["other_pe_kwh_per_m2"],
            "pv_pe": im["pv_pe_offset_kwh_per_m2"],
            "mains_gas": (doc.get("sap_energy_source") or {}).get("mains_gas"),
            "has_pv": bool((doc.get("sap_energy_source") or {}).get("photovoltaic_supply")),
            "tfa": tfa,
        })

    n = len(rows)
    print(f"decomposed {n} certs\n" + "=" * 70)

    def summ(key: str) -> str:
        vals = [r[key] for r in rows]
        return (f"median={stats.median(vals):+.3f} mean={stats.mean(vals):+.3f} "
                f"p10={sorted(vals)[n//10]:+.3f} p90={sorted(vals)[n*9//10]:+.3f}")

    print(f"PE  diff  (kWh/m2): {summ('pe_diff')}")
    print(f"PE  ratio (our/lod): {summ('pe_ratio')}")
    print(f"CO2 diff  (t/yr)  : {summ('co2_diff')}")
    print(f"CO2 ratio (our/lod): {summ('co2_ratio')}")
    print("=" * 70)

    # multiplicative vs additive: correlate pe_diff with lodged_pe magnitude
    lod = [r["lodged_pe"] for r in rows]
    dif = [r["pe_diff"] for r in rows]
    mean_lod, mean_dif = stats.mean(lod), stats.mean(dif)
    cov = sum((l - mean_lod) * (d - mean_dif) for l, d in zip(lod, dif)) / n
    var = sum((l - mean_lod) ** 2 for l in lod) / n
    slope = cov / var
    print(f"PE: regress pe_diff ~ lodged_pe -> slope={slope:+.4f} intercept={mean_dif - slope*mean_lod:+.3f}")
    print("  (slope~0 => additive constant; slope>0 => multiplicative)")
    print("=" * 70)

    # segment
    def seg(name: str, pred: Any) -> None:
        s = [r for r in rows if pred(r)]
        if not s:
            return
        print(f"  {name:28s} n={len(s):4d}  PEdiff_med={stats.median(r['pe_diff'] for r in s):+6.2f}  "
              f"PEratio_med={stats.median(r['pe_ratio'] for r in s):.3f}  "
              f"CO2diff_med={stats.median(r['co2_diff'] for r in s):+.3f}  "
              f"SAPerr_med={stats.median(r['sap_err'] for r in s):+.2f}")

    print("SEGMENTS:")
    seg("mains_gas=Y", lambda r: r["mains_gas"] in (True, 1, "Y"))
    seg("mains_gas=N", lambda r: r["mains_gas"] not in (True, 1, "Y"))
    seg("ALL", lambda r: True)
    print("=" * 70)
    print("PE COMPONENT MEANS (kWh/m2):")
    for k in ("space_pe", "hw_pe", "other_pe", "pv_pe", "our_pe", "lodged_pe"):
        print(f"  {k:12s} mean={stats.mean(r[k] for r in rows):7.2f} "
              f"median={stats.median(r[k] for r in rows):7.2f}")
    print("=" * 70)
    # Well-calibrated-SAP certs: energy is right, so PE diff isolates factor/scope.
    cal = [r for r in rows if abs(r["sap_err"]) < 0.3]
    print(f"WELL-CALIBRATED-SAP (|sap_err|<0.3) n={len(cal)}:")
    print(f"  PE diff   median={stats.median(r['pe_diff'] for r in cal):+.2f} "
          f"ratio={stats.median(r['pe_ratio'] for r in cal):.3f}")
    print(f"  CO2 diff  median={stats.median(r['co2_diff'] for r in cal):+.3f} "
          f"ratio={stats.median(r['co2_ratio'] for r in cal):.3f}")


if __name__ == "__main__":
    main()