Model/scripts/validate_epc_prediction.py

"""Leave-one-out accuracy harness for EPC Prediction (ADR-0029).

Runs entirely against the frozen postcode-clustered corpus
(`fetch_epc_prediction_corpus.py`). For every cert that has neighbours, it
drops that cert from its postcode cohort, predicts it from the rest using only
its *guaranteed* inputs (property type + built form), and compares the predicted
`EpcPropertyData` to the actual one.

Reports the ADR-0029 metrics:
  - classification rate: main wall construction (extend as coverage grows);
  - geometry residuals: floor area, window count + total window area, building
    parts (mean signed + mean absolute);
  - SAP reported three ways — predicted-then-calculated vs (a) the actual lodged
    SAP, (b) the calculator on the actual components, (c) the neighbour-mean SAP
    baseline (the number predict-then-calculate must beat).

USAGE
-----
    PYTHONPATH=. python scripts/validate_epc_prediction.py

Corpus dir: $EPC_PREDICTION_CORPUS (default /tmp/epc_prediction_corpus).
"""

from __future__ import annotations

import os
import statistics
from pathlib import Path
from typing import Optional

from datatypes.epc.domain.epc_property_data import EpcPropertyData
from domain.epc_prediction.comparable_properties import ComparableProperty
from domain.epc_prediction.validation import (
    evaluate_component_accuracy,
    iter_predictions,
)
from domain.sap10_calculator.calculator import Sap10Calculator, SapResult
from harness.epc_prediction_corpus import load_corpus

_KG_PER_TONNE: float = 1000.0

CORPUS = Path(os.environ.get("EPC_PREDICTION_CORPUS", "/tmp/epc_prediction_corpus"))


def _result(
    calculator: Sap10Calculator, epc: EpcPropertyData
) -> Optional[SapResult]:
    try:
        return calculator.calculate(epc)
    except Exception:  # noqa: BLE001 — some pictures don't score; count as misses
        return None


def _co2_tonnes(result: SapResult) -> float:
    """Calculated annual CO2 in tonnes, matching the lodged `co2_emissions_current`
    scale (see domain/property_baseline/performance.py)."""
    return result.co2_kg_per_yr / _KG_PER_TONNE


def main() -> None:
    cohorts = load_corpus(CORPUS)
    calculator = Sap10Calculator()

    # PRIMARY signal — Component Accuracy, calculator-free (the shared scorer).
    accuracy = evaluate_component_accuracy(cohorts)

    print(f"corpus: {CORPUS}")
    print(f"predicted {accuracy.targets} SAP-10.2 held-out targets\n")
    print("--- Component Accuracy (PRIMARY, calculator-independent) ---")
    for name, (hits, total) in accuracy.classification.items():
        if total:
            print(f"CLASSIFICATION  {name}: {hits}/{total} = {hits / total:.1%}")
    print()
    _floor_area_error(cohorts)
    _residual("floor_area (m2)", accuracy.residuals.get("floor_area", []))
    _residual("window_count", accuracy.residuals.get("window_count", []))
    _residual(
        "total_window_area (m2)", accuracy.residuals.get("total_window_area", [])
    )
    _residual("building_parts", accuracy.residuals.get("building_parts", []))
    _residual("door_count", accuracy.residuals.get("door_count", []))

    # SECONDARY guard — end-to-end vs API-lodged, calculator-FLOORED. Re-walks the
    # same held-out targets (one orchestration via iter_predictions).
    sap_vs_lodged: list[float] = []
    co2_vs_lodged: list[float] = []
    pei_vs_lodged: list[float] = []
    # Calculator floors — calc(actual) vs lodged — per metric. Each is the error
    # the end-to-end cannot beat (the API-path mapper/calculator residual, a
    # separate workstream), so it attributes how much of a metric's pred-vs-lodged
    # gap is the calculator vs the prediction. PEI carries a far larger floor than
    # SAP (~16 vs ~1.6 kWh/m2 / pts), so the headline PEI MAE must not be read as
    # pure prediction error (issue #1228).
    sap_floor: list[float] = []
    co2_floor: list[float] = []
    pei_floor: list[float] = []
    for predicted, actual in iter_predictions(cohorts):
        pred_result = _result(calculator, predicted)
        actual_result = _result(calculator, actual)
        lodged_sap = actual.energy_rating_current
        lodged_co2 = actual.co2_emissions_current
        lodged_pei = actual.energy_consumption_current
        if pred_result is not None:
            if lodged_sap is not None:
                sap_vs_lodged.append(
                    abs(pred_result.sap_score_continuous - lodged_sap)
                )
            if lodged_co2 is not None:
                co2_vs_lodged.append(abs(_co2_tonnes(pred_result) - lodged_co2))
            if lodged_pei is not None:
                pei_vs_lodged.append(
                    abs(pred_result.primary_energy_kwh_per_m2 - lodged_pei)
                )
        if actual_result is not None:
            if lodged_sap is not None:
                sap_floor.append(
                    abs(actual_result.sap_score_continuous - lodged_sap)
                )
            if lodged_co2 is not None:
                co2_floor.append(abs(_co2_tonnes(actual_result) - lodged_co2))
            if lodged_pei is not None:
                pei_floor.append(
                    abs(actual_result.primary_energy_kwh_per_m2 - lodged_pei)
                )

    print()
    print("--- End-to-end vs API-lodged (SECONDARY, calculator-FLOORED) ---")
    _sap_line("SAP |pred − lodged|", sap_vs_lodged)
    _sap_line("CO2 (t) |pred − lodged|", co2_vs_lodged)
    _sap_line("PEI (kWh/m2) |pred − lodged|", pei_vs_lodged)
    _sap_line("  floor: SAP |calc(actual) − lodged|", sap_floor)
    _sap_line("  floor: CO2 |calc(actual) − lodged|", co2_floor)
    _sap_line("  floor: PEI |calc(actual) − lodged|", pei_floor)


def _floor_area_error(cohorts: list[list[ComparableProperty]]) -> None:
    """Floor-area accuracy as MAE (m²) and MAPE (% of the actual), plus the
    typical (median actual) size — so the absolute error can be read relative to
    how big dwellings are. The predicted area is the cohort median, set
    independently of the geo/similarity weighting that drives the categoricals."""
    pairs = [
        (predicted.total_floor_area_m2, actual.total_floor_area_m2)
        for predicted, actual in iter_predictions(cohorts)
    ]
    valid = [(p, a) for p, a in pairs if a]
    if not valid:
        print("RESIDUAL  floor_area: (none)")
        return
    mae = statistics.mean(abs(p - a) for p, a in valid)
    mape = statistics.mean(abs(p - a) / a for p, a in valid)
    typical = statistics.median(a for _, a in valid)
    print(
        f"RESIDUAL  floor_area: MAE {mae:.2f} m2 | MAPE {mape:.1%} | "
        f"typical (median actual) {typical:.0f} m2 (n={len(valid)})"
    )


def _residual(label: str, values: list[float]) -> None:
    if not values:
        print(f"RESIDUAL  {label}: (none)")
        return
    mean_signed = statistics.mean(values)
    mean_abs = statistics.mean(abs(v) for v in values)
    print(f"RESIDUAL  {label}: mean {mean_signed:+.2f} | mean|·| {mean_abs:.2f} "
          f"(n={len(values)})")


def _sap_line(label: str, values: list[float]) -> None:
    if not values:
        print(f"{label}: (none)")
        return
    print(f"{label}: MAE {statistics.mean(values):.2f} | "
          f"median {statistics.median(values):.2f} (n={len(values)})")


if __name__ == "__main__":
    main()