"""ρ×ω×α スイープの輝度マップを「目で判断する」ための可視化.

指標（CoV / NCC 等）は分布を集約した量であり，形そのものの判断材料にはならない
（TEST_03 の主題）．本スクリプトは summary.csv の順位ではなく **輝度マップそのもの**
を並べ，どのパラメータで形がどう変わるか／変わらないかを目視できるようにする．

生成物:
  1. bsdf_grid_contact_sheet.png -- 全 60 条件 + 実機のコンタクトシート
     （行 = ρ，列 = (ω, α)．全セル共通カラースケールなので差がそのまま見える）
  2. bsdf_grid_rho_series.png    -- ω=0 純 diffuse の ρ 系列 + 実機．
     上段 = 輝度マップ，下段 = 実機との相対残差マップ [%]
  3. bsdf_grid_top_residual.png  -- スコア上位条件の残差マップ（発散カラーマップ）

残差の定義: 平均で正規化したうえで (sim - real) / mean(real) [%]．
正（赤）= シミュが明るすぎる領域，負（青）= 暗すぎる領域．

実行方法:
    uv run python scripts/build_bsdf_grid_visual.py [--top 4]
"""

from __future__ import annotations

import argparse
import csv
import sys
from pathlib import Path

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter, median_filter
from skimage.transform import resize

sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))

from sm_tias_lightsim.evaluation.roi_presets import SIM_ROI, extract_roi_xywh  # noqa: E402

matplotlib.rcParams["font.family"] = ["Meiryo", "Yu Gothic", "DejaVu Sans"]

ROOT = Path(__file__).resolve().parents[1]
SUMMARY_CSV = ROOT / "output" / "bsdf_grid_fix" / "summary.csv"
OUTPUT_BASE = ROOT / "output" / "bsdf_grid_fix"
REFERENCE_DIR = ROOT / "output" / "full_base_fix"
REAL_FULL = ROOT / "lightsim" / "luminance_maps" / "full_MiniTIAS_QM_20260601_124500.npy"
PANEL_DIR = ROOT / "output" / "_slide_panels"

RESIDUAL_LIMIT = 12.0  # 残差カラーバーの上下限 [%]


def load_rows() -> list[dict]:
    """summary.csv を数値化して読み込む."""
    if not SUMMARY_CSV.exists():
        raise SystemExit(f"エラー: {SUMMARY_CSV} がありません（スイープ未実行）")
    rows = []
    with SUMMARY_CSV.open(encoding="utf-8") as fh:
        for raw in csv.DictReader(fh):
            row: dict = {"job": raw["job"], "stage": raw["stage"]}
            for key, value in raw.items():
                if key not in ("job", "stage"):
                    row[key] = float(value) if value not in ("", None) else float("nan")
            rows.append(row)
    return rows


def _smooth(arr: np.ndarray, is_sim: bool) -> np.ndarray:
    """他パネルと同一の平滑化（median5 → gauss σ4）と ROI 抽出."""
    if arr.ndim == 3:
        arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2]
    arr = arr.astype(np.float64)
    if is_sim:
        arr = extract_roi_xywh(arr, SIM_ROI)
    return gaussian_filter(median_filter(arr, size=5), sigma=4)


def _job_map(job: str) -> np.ndarray:
    """ジョブ名から平滑化・mean 正規化済みの輝度マップを返す."""
    path = OUTPUT_BASE / job / "luminance.npy"
    if not path.exists():
        path = REFERENCE_DIR / "luminance.npy"
    img = _smooth(np.load(path), True)
    return img / img.mean()


def _real_map(shape: tuple[int, int]) -> np.ndarray:
    """実機マップを指定 shape に合わせ mean 正規化して返す."""
    real = _smooth(np.load(REAL_FULL), False)
    real = resize(real, shape, anti_aliasing=True)
    return real / real.mean()


def _residual(sim_norm: np.ndarray, real_norm: np.ndarray) -> np.ndarray:
    """mean 正規化済みどうしの相対残差 [%]（正＝シミュが明るすぎ）."""
    return 100.0 * (sim_norm - real_norm)


def contact_sheet(rows: list[dict]) -> Path:
    """全条件を 行=ρ × 列=(ω, α) で並べた一覧（共通カラースケール）."""
    rho_values = sorted({r["rho"] for r in rows})
    combos = sorted({(r["omega"], r["alpha"]) for r in rows})
    lookup = {(r["rho"], r["omega"], r["alpha"]): r for r in rows}

    ncols = len(combos) + 1  # 右端に実機列
    fig, axes = plt.subplots(len(rho_values), ncols,
                             figsize=(1.75 * ncols, 1.95 * len(rho_values)),
                             constrained_layout=True)
    sample = _job_map(rows[0]["job"])
    real = _real_map(sample.shape)
    vmin, vmax = np.percentile(real, [1, 99])

    for row_idx, rho in enumerate(rho_values):
        for col_idx, (omega, alpha) in enumerate(combos):
            ax = axes[row_idx][col_idx]
            entry = lookup.get((rho, omega, alpha))
            if entry is None:
                ax.axis("off")
                continue
            ax.imshow(_job_map(entry["job"]), cmap="hot", vmin=vmin, vmax=vmax)
            ax.set_title(f"NCC={entry['ncc_smooth']:.3f}", fontsize=6.5, pad=2)
            ax.set_xticks([])
            ax.set_yticks([])
            if col_idx == 0:
                ax.set_ylabel(f"ρ={rho:g}", fontsize=9)
            if row_idx == 0:
                label = "ω=0" if omega == 0 else f"ω={omega:g}\nα={alpha:g}"
                ax.text(0.5, 1.35, label, transform=ax.transAxes,
                        ha="center", va="bottom", fontsize=8)
        ax_real = axes[row_idx][-1]
        im = ax_real.imshow(real, cmap="hot", vmin=vmin, vmax=vmax)
        ax_real.set_xticks([])
        ax_real.set_yticks([])
        ax_real.set_title("実機", fontsize=7, pad=2)
        if row_idx == 0:
            ax_real.text(0.5, 1.35, "参照", transform=ax_real.transAxes,
                         ha="center", va="bottom", fontsize=8)

    fig.colorbar(im, ax=axes, shrink=0.5, label="mean 正規化輝度（全セル共通スケール）")
    fig.suptitle("全 60 条件の輝度マップ一覧 — 行=ρ，列=(ω, α)，右端=実機"
                 "（全セル共通カラースケール・mean 正規化）", fontsize=13)
    PANEL_DIR.mkdir(parents=True, exist_ok=True)
    path = PANEL_DIR / "bsdf_grid_contact_sheet.png"
    fig.savefig(path, dpi=140, bbox_inches="tight")
    plt.close(fig)
    return path


def rho_series(rows: list[dict]) -> Path:
    """ω=0 の ρ 系列を輝度マップ（上段）と残差マップ（下段）で並べる."""
    series = sorted([r for r in rows if r["omega"] == 0], key=lambda r: r["rho"])
    sample = _job_map(series[0]["job"])
    real = _real_map(sample.shape)
    vmin, vmax = np.percentile(real, [1, 99])

    ncols = len(series) + 1
    fig, axes = plt.subplots(2, ncols, figsize=(2.55 * ncols, 5.9),
                             constrained_layout=True)
    for col, entry in enumerate(series):
        sim = _job_map(entry["job"])
        axes[0][col].imshow(sim, cmap="hot", vmin=vmin, vmax=vmax)
        axes[0][col].set_title(f"ρ={entry['rho']:g}\nCoV={entry['cov']:.4f}  "
                               f"NCC={entry['ncc_smooth']:.3f}", fontsize=9)
        res = axes[1][col].imshow(_residual(sim, real), cmap="RdBu_r",
                                  vmin=-RESIDUAL_LIMIT, vmax=RESIDUAL_LIMIT)
        rms = float(np.sqrt(np.mean(_residual(sim, real) ** 2)))
        axes[1][col].set_title(f"残差 RMS={rms:.2f}%", fontsize=9)
    top = axes[0][-1].imshow(real, cmap="hot", vmin=vmin, vmax=vmax)
    axes[0][-1].set_title("実機 全灯\nCoV=0.0695", fontsize=9)
    axes[1][-1].axis("off")
    for ax in axes.ravel():
        if ax.axison:
            ax.set_xticks([])
            ax.set_yticks([])

    fig.colorbar(top, ax=axes[0], shrink=0.8, label="mean 正規化輝度")
    fig.colorbar(res, ax=axes[1], shrink=0.8, label="相対残差 [%]（赤=シミュが明るい）")
    fig.suptitle("純 diffuse（ω=0）の ρ 系列 — 上段: 輝度マップ / 下段: 実機との残差",
                 fontsize=13)
    PANEL_DIR.mkdir(parents=True, exist_ok=True)
    path = PANEL_DIR / "bsdf_grid_rho_series.png"
    fig.savefig(path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return path


def top_residual(rows: list[dict], top: int) -> Path:
    """スコア上位条件の残差マップを並べる（形の判断はここが主戦場）."""
    best = sorted(rows, key=lambda r: r["score"])[:top]
    sample = _job_map(best[0]["job"])
    real = _real_map(sample.shape)

    fig, axes = plt.subplots(1, top, figsize=(3.7 * top, 4.5), constrained_layout=True)
    for ax, entry in zip(np.atleast_1d(axes), best, strict=True):
        residual = _residual(_job_map(entry["job"]), real)
        im = ax.imshow(residual, cmap="RdBu_r", vmin=-RESIDUAL_LIMIT, vmax=RESIDUAL_LIMIT)
        rms = float(np.sqrt(np.mean(residual ** 2)))
        p95 = float(np.percentile(np.abs(residual), 95))
        ax.set_title(f"{entry['job']}\nRMS={rms:.2f}%  P95={p95:.2f}%  "
                     f"NCC={entry['ncc_smooth']:.3f}", fontsize=10)
        ax.set_xticks([])
        ax.set_yticks([])
    fig.colorbar(im, ax=axes, shrink=0.8, label="相対残差 [%]（赤=シミュが明るい / 青=暗い）")
    fig.suptitle("スコア上位条件の残差マップ — 赤青の構造が「材質では埋まらない形の差」",
                 fontsize=13)
    PANEL_DIR.mkdir(parents=True, exist_ok=True)
    path = PANEL_DIR / "bsdf_grid_top_residual.png"
    fig.savefig(path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return path


def main() -> int:
    """3 種の可視化を生成するメインエントリポイント."""
    parser = argparse.ArgumentParser(description="ρ×ω×α スイープの輝度マップ可視化")
    parser.add_argument("--top", type=int, default=4, help="残差マップに並べる上位条件数")
    args = parser.parse_args()

    rows = load_rows()
    print(f"読み込み: {len(rows)} 条件")
    print(f"一覧:     {contact_sheet(rows)}")
    print(f"ρ 系列:   {rho_series(rows)}")
    print(f"残差:     {top_residual(rows, args.top)}")
    return 0


if __name__ == "__main__":
    sys.exit(main())