"""立方体キャビティと球キャビティを内面積で対応づけて比較する図を生成する（TEST_16）.

球（similar_radial・10 条件・レンダ済み）と立方体（サイズスイープ）は「同じ大きさ」の
定義が自明でない．多重反射の理論では放射輝度は総内面積と開口比で決まるため，
**内面積を横軸**に取って両者の指標カーブを重ねる（積分球としての比較）．

評価の作法は TEST_15 で確立したものを踏襲する:
    - 生の CoV では比較しない（MC ノイズ支配）．構造 CoV（decompose_cov）を使う
    - 形の比較は同心円（輪帯）プロファイルの方位平均で行い，ピークで正規化する
    - 輝度マップは 20 分割の輪帯平均の最大値を 1.0 に正規化する
    - angular_ratio は分母（動径分散）で跳ねるため，angular/radial 分散の内訳を併記する

解析関数は既存スクリプト（build_radial_comparison / build_radius_sweep_visual）から
import し，輪帯・正規化・CoV 分解の定義が図をまたいで厳密に一致するようにする．

実行方法:
    uv run python scripts/build_cube_sphere_comparison.py
"""

from __future__ import annotations

import argparse
import csv
import sys
from pathlib import Path

import matplotlib
import numpy as np

matplotlib.use("Agg")
import matplotlib.pyplot as plt  # noqa: E402

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

from build_radial_comparison import (  # noqa: E402
    normalize_to_peak,
    radial_stats,
    ring_min_max_ratio,
)
from build_radius_sweep_visual import decompose_cov, normalized_map  # noqa: E402
from sweep_cube_size import OUTPUT_BASE as CUBE_OUTPUT_BASE  # noqa: E402
from sweep_cube_size import RUN_DIR as CUBE_RUN_DIR  # noqa: E402
from sweep_cube_size import job_name as cube_job_name  # noqa: E402
from sweep_sphere_radius import OUTPUT_BASE as SPHERE_OUTPUT_BASE  # noqa: E402
from sweep_sphere_radius import job_name as sphere_job_name  # noqa: E402

plt.rcParams["font.family"] = ["Yu Gothic", "MS Gothic", "DejaVu Sans"]

# 球側は TEST_15 で最良だった similar_radial（リム水平出し）を比較対象にする
SPHERE_RUN_DIR = "board_baffle_similar_radial"

# 輪帯の分割数（評価パイプライン・既存比較図と同一）
N_BINS = 20

CUBE_COLOR = "tab:red"
SPHERE_COLOR = "tab:blue"


def read_summary(path: Path, key: str) -> list[dict[str, float]]:
    """summary.csv を読み，key 列の昇順の辞書リストで返す."""
    with path.open(encoding="utf-8") as fh:
        rows = [
            {k: (float(v) if v else float("nan")) for k, v in row.items()}
            for row in csv.DictReader(fh)
        ]
    return sorted(rows, key=lambda r: r[key])


def read_full_metrics(path: Path) -> dict[str, float]:
    """metrics.csv の 1 行目を全列 float で返す（angular/radial 分散の内訳併記用）."""
    with path.open(encoding="utf-8") as fh:
        row = next(csv.DictReader(fh))
    return {k: (float(v) if v else float("nan")) for k, v in row.items()}


def collect(
    rows: list[dict[str, float]], run_dir: Path, size_key: str, job_name_fn
) -> list[dict[str, float]]:
    """summary の各行に輪帯 min/max・構造 CoV・分散内訳を実測で付加して返す."""
    result = []
    for row in rows:
        job_dir = run_dir / job_name_fn(row[size_key])
        lum = job_dir / "luminance.npy"
        if not lum.exists():
            continue
        centers, means, _ = radial_stats(lum, N_BINS)
        _, structure, noise = decompose_cov(lum)
        enriched = dict(row)
        enriched["ring_min_max"] = ring_min_max_ratio(centers, means)
        enriched["structure_cov"] = structure
        enriched["noise_cov"] = noise
        enriched.update(
            {
                k: v
                for k, v in read_full_metrics(job_dir / "metrics.csv").items()
                if k in ("angular_variance", "radial_variance", "pca_primary_axis_angle_deg")
            }
        )
        result.append(enriched)
    return result


def nearest_by_area(area: float, rows: list[dict[str, float]]) -> dict[str, float]:
    """内面積が最も近い行を返す（対数距離で比較）."""
    return min(rows, key=lambda r: abs(np.log(r["cavity_area_mm2"] / area)))


def draw_cube_maps(fig, gs, cube_rows: list[dict[str, float]], run_dir: Path) -> None:
    """各サイズの輝度マップを「輪帯平均の最大値=1.0」で正規化して並べる（共通スケール）."""
    maps = [
        (row["size_mm"], normalized_map(run_dir / cube_job_name(row["size_mm"]) / "luminance.npy"))
        for row in cube_rows
    ]
    if not maps:
        return
    stacked = np.concatenate([m.ravel() for _, m in maps])
    vmin = float(np.percentile(stacked, 0.5))
    for i, (size, data) in enumerate(maps[:5]):
        ax = fig.add_subplot(gs[0, i])
        image = ax.imshow(data, cmap="hot", vmin=vmin, vmax=1.0)
        ax.set_title(f"L={size:g}mm", fontsize=10)
        ax.axis("off")
        if i == len(maps[:5]) - 1:
            fig.colorbar(image, ax=ax, fraction=0.046, label="最大輪帯（20 分割）=1.0")


def draw_cube_profiles(ax, cube_rows: list[dict[str, float]], run_dir: Path) -> None:
    """立方体全サイズの輪帯プロファイルを重ね描きする（ピーク=1 正規化・標準誤差帯付き）."""
    colors = plt.cm.viridis(np.linspace(0.0, 0.9, len(cube_rows)))
    for row, color in zip(cube_rows, colors, strict=True):
        lum = run_dir / cube_job_name(row["size_mm"]) / "luminance.npy"
        centers, means, sems = radial_stats(lum, N_BINS)
        means, sems = normalize_to_peak(centers, means, sems)
        ax.plot(centers, means, color=color, lw=1.6, label=f"L={row['size_mm']:g}")
        ax.fill_between(centers, means - sems, means + sems, color=color, alpha=0.35, lw=0)
    ax.axhline(1.0, color="gray", lw=0.8, ls="--")
    ax.set_xlabel("正規化楕円距離（0=中心, 1=ROI 四隅）")
    ax.set_ylabel("輪帯平均輝度（ピーク=1）")
    ax.set_title("立方体: 同心円プロファイルのサイズ依存（帯は ±1 標準誤差）", fontsize=11)
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3)


def draw_metric_vs_area(
    ax,
    cube_rows: list[dict[str, float]],
    sphere_rows: list[dict[str, float]],
    key: str,
    ylabel: str,
    title: str,
    ideal: float | None,
) -> None:
    """指標を内面積（対数軸）に対して球・立方体で重ね描きする."""
    for rows, color, marker, label in (
        (sphere_rows, SPHERE_COLOR, "o", "球 similar_radial（r=88〜250）"),
        (cube_rows, CUBE_COLOR, "s", "立方体（L スイープ）"),
    ):
        areas = [r["cavity_area_mm2"] for r in rows]
        ax.plot(areas, [r[key] for r in rows], marker=marker, color=color, label=label)
    if ideal is not None:
        ax.axhline(ideal, color="gray", lw=0.8, ls="--")
    ax.set_xscale("log")
    ax.set_xlabel("キャビティ内面積 (mm²・対数)")
    ax.set_ylabel(ylabel)
    ax.set_title(title, fontsize=11)
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3, which="both")


def draw_matched_pair(
    ax,
    cube_row: dict[str, float],
    sphere_row: dict[str, float],
    cube_run: Path,
    sphere_run: Path,
) -> None:
    """内面積が最も近い球と立方体の輪帯プロファイルを 1 枚に重ねる."""
    targets = [
        (
            cube_run / cube_job_name(cube_row["size_mm"]) / "luminance.npy",
            f"立方体 L={cube_row['size_mm']:g}（{cube_row['cavity_area_mm2'] / 1000:.0f}k mm²）",
            CUBE_COLOR,
            "-",
        ),
        (
            sphere_run / sphere_job_name(sphere_row["radius_mm"]) / "luminance.npy",
            f"球 r={sphere_row['radius_mm']:g}（{sphere_row['cavity_area_mm2'] / 1000:.0f}k mm²）",
            SPHERE_COLOR,
            "--",
        ),
    ]
    for path, label, color, ls in targets:
        centers, means, sems = radial_stats(path, N_BINS)
        means, sems = normalize_to_peak(centers, means, sems)
        ax.plot(centers, means, color=color, ls=ls, lw=2.0, label=label)
        ax.fill_between(centers, means - sems, means + sems, color=color, alpha=0.3, lw=0)
    ax.axhline(1.0, color="gray", lw=0.8, ls="--")
    ax.set_xlabel("正規化楕円距離")
    ax.set_ylabel("輪帯平均輝度（ピーク=1）")
    ax.set_title("内面積で対応づけた直接比較", fontsize=11)
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3)


def print_table(cube_rows: list[dict[str, float]], sphere_rows: list[dict[str, float]]) -> None:
    """比較表を標準出力に書く（angular/radial 分散の内訳を必ず併記する）."""
    print("\n=== 立方体（実測） ===")
    print(
        f"{'L':>6} {'内面積mm2':>10} {'構造CoV':>8} {'輪帯min/max':>11} "
        f"{'angular_ratio':>13} {'ang_var':>8} {'rad_var':>8} {'PCA_ecc':>8}"
    )
    for r in cube_rows:
        print(
            f"{r['size_mm']:6.0f} {r['cavity_area_mm2']:10.0f} {r['structure_cov']:8.4f} "
            f"{r['ring_min_max']:11.4f} {r['angular_ratio']:13.4f} "
            f"{r['angular_variance']:8.2f} {r['radial_variance']:8.2f} "
            f"{r['pca_eccentricity']:8.4f}"
        )
    print("\n=== 球 similar_radial（参照） ===")
    print(
        f"{'r':>6} {'内面積mm2':>10} {'構造CoV':>8} {'輪帯min/max':>11} "
        f"{'angular_ratio':>13} {'ang_var':>8} {'rad_var':>8} {'PCA_ecc':>8}"
    )
    for r in sphere_rows:
        print(
            f"{r['radius_mm']:6.0f} {r['cavity_area_mm2']:10.0f} {r['structure_cov']:8.4f} "
            f"{r['ring_min_max']:11.4f} {r['angular_ratio']:13.4f} "
            f"{r['angular_variance']:8.2f} {r['radial_variance']:8.2f} "
            f"{r['pca_eccentricity']:8.4f}"
        )

    # 同じ内面積での球の値を対数内挿して，立方体との比を出す（1 対 1 対応が無いため）
    log_area = np.log([r["cavity_area_mm2"] for r in sphere_rows])
    print("\n=== 同一内面積での比較（球は対数内挿） ===")
    print(
        f"{'L':>6} {'構造CoV 立方体':>13} {'球（内挿）':>10} {'比':>6}   "
        f"{'輪帯min/max 立方体':>17} {'球（内挿）':>10}"
    )
    for r in cube_rows:
        x = np.log(r["cavity_area_mm2"])
        sphere_cov = float(np.interp(x, log_area, [s["structure_cov"] for s in sphere_rows]))
        sphere_ring = float(np.interp(x, log_area, [s["ring_min_max"] for s in sphere_rows]))
        print(
            f"{r['size_mm']:6.0f} {r['structure_cov']:13.4f} {sphere_cov:10.4f} "
            f"{r['structure_cov'] / sphere_cov:6.2f}   "
            f"{r['ring_min_max']:17.4f} {sphere_ring:10.4f}"
        )


def main() -> int:
    """立方体 vs 球の比較図を生成する."""
    parser = argparse.ArgumentParser(description="立方体と球のキャビティ比較図を生成")
    parser.add_argument("--output", type=Path, default=None, help="出力 PNG パス")
    args = parser.parse_args()

    cube_run = CUBE_OUTPUT_BASE / CUBE_RUN_DIR
    sphere_run = SPHERE_OUTPUT_BASE / SPHERE_RUN_DIR
    cube_summary = cube_run / "summary.csv"
    sphere_summary = sphere_run / "summary.csv"
    for path in (cube_summary, sphere_summary):
        if not path.exists():
            print(f"エラー: {path} がありません（先にスイープを実行してください）")
            return 1

    cube_rows = collect(read_summary(cube_summary, "size_mm"), cube_run, "size_mm", cube_job_name)
    sphere_rows = collect(
        read_summary(sphere_summary, "radius_mm"), sphere_run, "radius_mm", sphere_job_name
    )
    output = args.output or cube_run / "cube_vs_sphere.png"

    fig = plt.figure(figsize=(20, 14))
    gs = fig.add_gridspec(3, 6, height_ratios=[1.2, 1, 1])
    draw_cube_maps(fig, gs, cube_rows, cube_run)
    draw_cube_profiles(fig.add_subplot(gs[1, 0:2]), cube_rows, cube_run)
    draw_metric_vs_area(
        fig.add_subplot(gs[1, 2:4]),
        cube_rows,
        sphere_rows,
        "ring_min_max",
        "輪帯 min/max（1 が理想）",
        "ドームの深さ vs 内面積",
        1.0,
    )
    draw_metric_vs_area(
        fig.add_subplot(gs[1, 4:6]),
        cube_rows,
        sphere_rows,
        "structure_cov",
        "構造 CoV（小さいほど均一）",
        "構造 CoV vs 内面積",
        None,
    )
    # 内面積の近い代表 3 ペア（小・中・大）で直接重ねる
    picks = [cube_rows[0], cube_rows[len(cube_rows) // 2], cube_rows[-1]]
    for i, cube_row in enumerate(picks):
        sphere_row = nearest_by_area(cube_row["cavity_area_mm2"], sphere_rows)
        draw_matched_pair(
            fig.add_subplot(gs[2, 2 * i : 2 * i + 2]), cube_row, sphere_row, cube_run, sphere_run
        )
    fig.suptitle(
        "立方体キャビティ vs 球キャビティ（board_baffle・内面積で対応づけ）"
        "　輪帯 20 分割・方位平均・ピーク正規化",
        fontsize=15,
    )
    fig.tight_layout()
    output.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(output, dpi=100)
    print(f"出力: {output}")
    print_table(cube_rows, sphere_rows)
    return 0


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