Newer
Older
SmTIAS-LightSim / scripts / build_cavity_comparison.py
"""キャビティ形状(ロフト接合の設計)どうしを半径軸で比較する図を生成する.

同じ LED 構成(board_baffle)のまま,ロフト接合の設計だけを変えた系列を重ねて描く.
比較は CoV ではなく**同心円(輪帯)プロファイル**を基本にする.本スイープは絶対輝度が
半径で 1 桁変わり生の CoV が MC ノイズに支配されるため([TEST_15](../docs/06_TEST/)),
方位平均でノイズを 1/sqrt(N) に落とした量でないと半径間の比較ができない.

出す量:
    内側の落ち込み (ピーク→d=0.5) : 球のサイズが支配する成分
    外側の落ち込み (d=0.5→0.95)   : ロフトが作る「井戸」が支配する成分
    輪帯 min/max                  : ドーム全体の深さ(1 が理想)
    構造 CoV                      : ノイズを除いた真の不均一(decompose_cov)

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

from __future__ import annotations

import argparse
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  # noqa: E402
from sweep_sphere_radius import CAVITIES, OUTPUT_BASE, RADII, job_name, run_dir_name  # noqa: E402

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

# 表示名と色(キャビティ形状ごと)
STYLES = {
    "fixed_cut": ("接続円固定・3 次(喉あり)", "tab:red", "--"),
    "similar_g2": ("相似・5 次 G2(喉あり)", "tab:green", "-."),
    "similar_radial": ("相似・リム水平出し(喉なし)", "tab:blue", "-"),
}
# 輪帯プロファイルを重ねて見せる代表半径
PROFILE_RADII = (100.0, 250.0)


def collect(variant: str, cavity: str, bins: int) -> dict[str, list[float]]:
    """1 つのキャビティ形状について,半径ごとの各指標を集める(欠けている半径は飛ばす)."""
    out: dict[str, list[float]] = {k: [] for k in
                                   ("radius", "inner", "outer", "ratio", "structure")}
    for radius in RADII:
        path = OUTPUT_BASE / run_dir_name(variant, cavity) / job_name(radius) / "luminance.npy"
        if not path.exists():
            continue
        centers, means, sems = radial_stats(path, bins)
        ratio = ring_min_max_ratio(centers, means)
        _, structure, _ = decompose_cov(path)
        means, _ = normalize_to_peak(centers, means, sems)
        valid = ~np.isnan(means)
        mid = float(np.interp(0.5, centers[valid], means[valid]))
        edge = float(np.interp(0.95, centers[valid], means[valid]))
        out["radius"].append(radius)
        out["inner"].append((1.0 - mid) * 100.0)
        out["outer"].append((mid - edge) * 100.0)
        out["ratio"].append(ratio)
        out["structure"].append(structure)
    return out


def draw_profiles(ax, variant: str, bins: int, radius: float) -> None:
    """1 つの半径について,各キャビティ形状の輪帯プロファイルを重ねる."""
    for cavity, (label, color, style) in STYLES.items():
        path = OUTPUT_BASE / run_dir_name(variant, cavity) / job_name(radius) / "luminance.npy"
        if not path.exists():
            continue
        centers, means, sems = radial_stats(path, bins)
        means, sems = normalize_to_peak(centers, means, sems)
        ax.plot(centers, means, style, color=color, lw=2.0, label=label)
        ax.fill_between(centers, means - sems, means + sems, color=color, alpha=0.25, lw=0)
    ax.axhline(1.0, color="gray", lw=0.8, ls=":")
    ax.set_xlabel("正規化楕円距離")
    ax.set_ylabel("輪帯平均輝度(ピーク=1)")
    ax.set_title(f"同心円プロファイル r={radius:g}mm", fontsize=12)
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)


def draw_series(ax, series: dict[str, dict], key: str, title: str, ylabel: str,
                ideal: float | None = None) -> None:
    """指標 1 つについて,各キャビティ形状の半径依存カーブを描く."""
    for cavity, data in series.items():
        if not data["radius"]:
            continue
        label, color, style = STYLES[cavity]
        ax.plot(data["radius"], data[key], style, color=color, marker="o", lw=2.0, label=label)
    if ideal is not None:
        ax.axhline(ideal, color="gray", lw=0.8, ls="--")
    ax.set_xlabel("球内半径 (mm)")
    ax.set_ylabel(ylabel)
    ax.set_title(title, fontsize=12)
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)


def main() -> int:
    """キャビティ形状の比較図を生成する."""
    parser = argparse.ArgumentParser(description="キャビティ形状の比較図を生成")
    parser.add_argument("--variant", default="board_baffle", help="LED 構成")
    parser.add_argument("--bins", type=int, default=40, help="輪帯の分割数")
    parser.add_argument("--output", type=Path, default=None, help="出力 PNG パス")
    args = parser.parse_args()

    series = {c: collect(args.variant, c, args.bins) for c in CAVITIES}
    output = args.output or OUTPUT_BASE / "cavity_comparison.png"

    fig, axes = plt.subplots(2, 3, figsize=(20, 12))
    for ax, radius in zip(axes[0][:2], PROFILE_RADII, strict=True):
        draw_profiles(ax, args.variant, args.bins, radius)
    draw_series(axes[0][2], series, "ratio", "ドームの深さ(1 が理想)", "輪帯 min/max", 1.0)
    draw_series(axes[1][0], series, "inner", "内側の落ち込み(ピーク→d=0.5)", "落ち込み [%]")
    draw_series(axes[1][1], series, "outer", "外側の落ち込み(d=0.5→0.95)", "落ち込み [%]")
    draw_series(axes[1][2], series, "structure", "構造 CoV(MC ノイズ除去後)", "構造 CoV")

    fig.suptitle(
        f"ロフト接合の設計比較 [{args.variant}] 同心円平均でノイズを除いた比較", fontsize=15
    )
    fig.tight_layout()
    output.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(output, dpi=100)
    print(f"出力: {output}")
    for cavity, data in series.items():
        if data["radius"]:
            best = int(np.argmax(data["ratio"]))
            print(f"  {cavity:16s} {len(data['radius']):2d} 条件  "
                  f"最良 r={data['radius'][best]:.0f} で 輪帯 min/max={data['ratio'][best]:.4f} "
                  f"構造 CoV={data['structure'][best]:.4f}")
    return 0


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