#!/usr/bin/env python3 from __future__ import annotations import argparse import html import json from pathlib import Path import sys from typing import Dict, List, Tuple import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd ROOT = Path(__file__).resolve().parents[1] sys.path.insert(0, str(ROOT)) from src.config_utils import prepare_output_dirs, read_config from src.models import get_model_params PRIMARY_CASES = [ ("pr_o", "RMAE_percent", True, "Precipitation RMAE (%)"), ("pr_o", "CORR", False, "Precipitation CORR"), ("t2_o", "MAE", True, "T2m MAE (°C)"), ("t2_o", "RMSE", True, "T2m RMSE (°C)"), ("t2_o", "CORR", False, "T2m CORR"), ] def _safe_wmean(g: pd.DataFrame) -> float: w = pd.to_numeric(g["n_grid"], errors="coerce").astype(float) x = pd.to_numeric(g["mean_value"], errors="coerce").astype(float) ok = np.isfinite(x) & np.isfinite(w) & (w > 0) if ok.sum() == 0: return np.nan return float(np.sum(x[ok] * w[ok]) / np.sum(w[ok])) def build_weighted_domain(summary_all: pd.DataFrame) -> pd.DataFrame: region_rows = summary_all[summary_all["region"].astype(str).str.match(r"^R[1-7]$", na=False)].copy() for col in ["stack_scope", "base_methods", "base_source_models"]: if col not in region_rows.columns: region_rows[col] = "" region_rows[col] = region_rows[col].fillna("") group_cols = [ "source_model", "lead", "target", "spatial_mode", "predictor_set", "method", "metric", "stack_scope", "base_methods", "base_source_models", ] return ( region_rows.groupby(group_cols, dropna=False)[["mean_value", "n_grid"]] .apply(_safe_wmean) .reset_index(name="weighted_mean") ) def best_vs_raw(weighted: pd.DataFrame) -> pd.DataFrame: rows = [] for target, metric, small_better, _title in PRIMARY_CASES: x = weighted[(weighted["target"] == target) & (weighted["metric"] == metric)].copy() for (source_model, lead), g in x.groupby(["source_model", "lead"]): raw = g[g["method"].astype(str).str.lower() == "raw"].copy() if raw.empty: continue best = g.sort_values("weighted_mean", ascending=small_better).iloc[0] raw_value = float(raw["weighted_mean"].iloc[0]) best_value = float(best["weighted_mean"]) if small_better: improvement = 100.0 * (raw_value - best_value) / raw_value absolute_change = raw_value - best_value else: improvement = 100.0 * (best_value - raw_value) / abs(raw_value) absolute_change = best_value - raw_value rows.append({ "source_model": source_model, "lead": int(lead), "target": target, "metric": metric, "raw_value": raw_value, "best_value": best_value, "absolute_change": absolute_change, "improvement_percent": improvement, "best_spatial_mode": best["spatial_mode"], "best_predictor_set": best["predictor_set"], "best_method": best["method"], "best_stack_scope": best.get("stack_scope", ""), }) return pd.DataFrame(rows).sort_values(["target", "metric", "source_model", "lead"]) if rows else pd.DataFrame() def best_per_source_lead(weighted: pd.DataFrame) -> pd.DataFrame: rows = [] for target, metric, small_better, _title in PRIMARY_CASES: x = weighted[(weighted["target"] == target) & (weighted["metric"] == metric)].copy() if x.empty: continue x = x.sort_values("weighted_mean", ascending=small_better) best = x.groupby(["source_model", "lead"], as_index=False).first() best["rank_metric"] = metric best["small_better"] = small_better rows.append(best) return pd.concat(rows, ignore_index=True) if rows else pd.DataFrame() def best_by_method(weighted: pd.DataFrame) -> pd.DataFrame: rows = [] for target, metric, small_better, _title in PRIMARY_CASES: x = weighted[(weighted["target"] == target) & (weighted["metric"] == metric)].copy() if x.empty: continue x = x.sort_values("weighted_mean", ascending=small_better) best = x.groupby(["source_model", "lead", "target", "method", "metric"], as_index=False).first() best["small_better"] = small_better rows.append(best) return pd.concat(rows, ignore_index=True) if rows else pd.DataFrame() def method_ranking(best_method_df: pd.DataFrame) -> pd.DataFrame: rows = [] for target, metric, small_better, _title in PRIMARY_CASES: x = best_method_df[(best_method_df["target"] == target) & (best_method_df["metric"] == metric)].copy() if x.empty: continue agg = x.groupby("method", as_index=False).agg( mean_score=("weighted_mean", "mean"), median_score=("weighted_mean", "median"), n_cases=("weighted_mean", "count"), ) agg = agg.sort_values("mean_score", ascending=small_better).reset_index(drop=True) agg.insert(0, "rank", np.arange(1, len(agg)+1)) agg.insert(1, "target", target) agg.insert(2, "metric", metric) agg["small_better"] = small_better rows.append(agg) return pd.concat(rows, ignore_index=True) if rows else pd.DataFrame() def stacking_vs_individual(weighted: pd.DataFrame) -> pd.DataFrame: rows = [] for target, metric, small_better, _title in PRIMARY_CASES: x = weighted[(weighted["target"] == target) & (weighted["metric"] == metric)].copy() if x.empty: continue for lead, g in x.groupby("lead"): indiv = g[g["stack_scope"].astype(str) == ""].copy() stack = g[g["stack_scope"].astype(str) != ""].copy() if indiv.empty or stack.empty: continue bi = indiv.sort_values("weighted_mean", ascending=small_better).iloc[0] bs = stack.sort_values("weighted_mean", ascending=small_better).iloc[0] vi = float(bi["weighted_mean"]); vs = float(bs["weighted_mean"]) gain = (vi - vs) if small_better else (vs - vi) rel = 100.0 * gain / abs(vi) if vi != 0 else np.nan rows.append({ "lead": int(lead), "target": target, "metric": metric, "best_individual_source": bi["source_model"], "best_individual_method": bi["method"], "best_individual_spatial_mode": bi["spatial_mode"], "best_individual_predictor_set": bi["predictor_set"], "best_individual_score": vi, "best_stack_source": bs["source_model"], "best_stack_method": bs["method"], "best_stack_scope": bs["stack_scope"], "best_stack_spatial_mode": bs["spatial_mode"], "best_stack_score": vs, "absolute_gain_stack_vs_individual": gain, "relative_gain_percent_stack_vs_individual": rel, "winner": "stacking" if gain > 0 else "individual", }) return pd.DataFrame(rows) def method_hyperparameter_catalogue(cfg: Dict) -> pd.DataFrame: methods = list(dict.fromkeys(list(cfg.get("methods", [])) + list(cfg.get("stacking", {}).get("base_methods", [])) + list(cfg.get("stacking", {}).get("meta_methods", [])))) rows = [] for m in methods: rows.append({"method": m, "hyperparameters_json": json.dumps(get_model_params(m, cfg), ensure_ascii=False, sort_keys=True)}) return pd.DataFrame(rows) def _fmt_df(df: pd.DataFrame, max_rows: int = 300) -> str: if df.empty: return "

No data.

" show = df.copy().head(max_rows) for c in show.columns: if pd.api.types.is_float_dtype(show[c]): show[c] = show[c].map(lambda v: "" if pd.isna(v) else f"{v:.4f}") return show.to_html(index=False, escape=False, classes="data-table") def _plot_method_ranking(rank: pd.DataFrame, report_fig_dir: Path) -> List[Path]: out_files = [] for target, metric, small_better, title in PRIMARY_CASES: d = rank[(rank["target"] == target) & (rank["metric"] == metric)].copy() if d.empty: continue d = d.sort_values("mean_score", ascending=small_better) fig, ax = plt.subplots(figsize=(9, max(4, 0.35*len(d))), constrained_layout=True) ax.barh(d["method"], d["mean_score"].astype(float)) ax.invert_yaxis() ax.set_xlabel(metric) ax.set_title(title + (": lower is better" if small_better else ": higher is better")) fname = report_fig_dir / f"method_ranking_{target}_{metric}.png" fig.savefig(fname, dpi=150) plt.close(fig) out_files.append(fname) return out_files def write_html_report(report_file: Path, cfg: Dict, weighted: pd.DataFrame, best_raw: pd.DataFrame, best_lead: pd.DataFrame, best_method_df: pd.DataFrame, rank: pd.DataFrame, stack_cmp: pd.DataFrame, hp: pd.DataFrame, stack_weights: pd.DataFrame, fig_files: List[Path]) -> None: report_file.parent.mkdir(parents=True, exist_ok=True) rel_figs = [p.relative_to(report_file.parent).as_posix() for p in fig_files] cfg_summary = {k: cfg.get(k) for k in ["data_dir", "mask_file", "source_models", "leads", "targets", "spatial_modes", "predictor_sets", "methods", "train_years", "test_years"]} cfg_summary["stacking"] = cfg.get("stacking", {}) top_tables = [] for target, metric, small_better, title in PRIMARY_CASES: d = weighted[(weighted["target"] == target) & (weighted["metric"] == metric)].copy().sort_values("weighted_mean", ascending=small_better).head(30) top_tables.append((title, d)) html_parts = ["""PA3 Stacking BC Report"""] html_parts.append("

PA3 ML/DL Bias Correction + Stacking Ensemble Report

") html_parts.append("
Bias correction uses monthly anomaly-to-anomaly learning. Climatology is computed from the training period only. Metrics are calculated on the independent test period only. RMAE (%) is reported only for precipitation; temperature uses MAE/RMSE/CORR.
") html_parts.append("

Table of contents

") for a,n in [("config","Experiment configuration"),("hyper","Model hyperparameters"),("raw","Best correction vs RAW"),("lead","Best configuration per source/lead"),("rank","Method ranking"),("stack","Stacking vs individual models"),("weights","Stacking weights"),("top","Top configurations"),("figs","Figures")]: html_parts.append(f"{html.escape(n)}") html_parts.append("
") html_parts.append("

Experiment configuration

" + html.escape(pd.Series(cfg_summary).to_string()) + "
") html_parts.append("

Model hyperparameters

This table records the actual hyperparameters used by each method variant.

" + _fmt_df(hp, 200)) html_parts.append("

Best correction vs RAW

") for target, metric, _small, title in PRIMARY_CASES: d = best_raw[(best_raw["target"]==target)&(best_raw["metric"]==metric)] html_parts.append(f"

{html.escape(title)}

" + _fmt_df(d, 200)) html_parts.append("

Best configuration per source model and lead

") for target, metric, _small, title in PRIMARY_CASES: d = best_lead[(best_lead["target"]==target)&(best_lead["rank_metric"]==metric)] html_parts.append(f"

{html.escape(title)}

" + _fmt_df(d, 300)) html_parts.append("

Method ranking

Each method is first allowed to choose its best spatial mode and predictor setup for each source/lead case. The scores are then averaged.
") for target, metric, _small, title in PRIMARY_CASES: d = rank[(rank["target"]==target)&(rank["metric"]==metric)] html_parts.append(f"

{html.escape(title)}

" + _fmt_df(d, 200)) html_parts.append("

Stacking vs individual models

Positive gain means the best stacking configuration improves upon the best individual RegCM/clWRF correction for that target/lead/metric.
" + _fmt_df(stack_cmp, 300)) html_parts.append("

Stacking weights

Weights are reported for linear meta-models when available. They are month/region-specific because PA3 trains separate monthly models.

" + _fmt_df(stack_weights, 500)) html_parts.append("

Top configurations

") for title, d in top_tables: html_parts.append(f"

{html.escape(title)}

" + _fmt_df(d, 30)) html_parts.append("

Figures

") for rel in rel_figs: html_parts.append(f"

{html.escape(Path(rel).stem)}

{html.escape(rel)}") html_parts.append("") report_file.write_text("\n".join(html_parts), encoding="utf-8") def main() -> None: parser = argparse.ArgumentParser(description="Analyze PA3 experiment summaries and build HTML report.") parser.add_argument("--config", default="config/experiment.yaml") args = parser.parse_args() cfg_path = Path(args.config).resolve() cfg = read_config(cfg_path) dirs = prepare_output_dirs(cfg, base_dir=cfg_path.parent.parent) report_dir = dirs["output"] / "report" fig_dir = report_dir / "figures"; fig_dir.mkdir(parents=True, exist_ok=True) summary_file = dirs["metrics"] / "summary_all.csv" if not summary_file.exists(): raise FileNotFoundError(f"Missing summary file: {summary_file}") summary = pd.read_csv(summary_file) weighted = build_weighted_domain(summary); weighted.to_csv(dirs["metrics"] / "summary_weighted_domain.csv", index=False) br = best_vs_raw(weighted); br.to_csv(dirs["metrics"] / "best_vs_raw.csv", index=False) bpl = best_per_source_lead(weighted); bpl.to_csv(dirs["metrics"] / "best_per_source_lead.csv", index=False) bbm = best_by_method(weighted); bbm.to_csv(dirs["metrics"] / "best_by_method.csv", index=False) rank = method_ranking(bbm); rank.to_csv(dirs["metrics"] / "method_ranking.csv", index=False) stack_cmp = stacking_vs_individual(weighted); stack_cmp.to_csv(dirs["metrics"] / "stacking_vs_individual.csv", index=False) hp = method_hyperparameter_catalogue(cfg); hp.to_csv(dirs["metrics"] / "method_hyperparameters.csv", index=False) weight_files = sorted(dirs["metrics"].glob("stack_weights_*.csv")) stack_weights = pd.concat([pd.read_csv(f) for f in weight_files], ignore_index=True) if weight_files else pd.DataFrame() if not stack_weights.empty: stack_weights.to_csv(dirs["metrics"] / "stack_weights_all.csv", index=False) figs = _plot_method_ranking(rank, fig_dir) write_html_report(report_dir / "index.html", cfg, weighted, br, bpl, bbm, rank, stack_cmp, hp, stack_weights, figs) print(f"Wrote HTML report: {report_dir / 'index.html'}") if __name__ == "__main__": main()