#!/usr/bin/env python3
"""
概要:
    熱伝導率フィッティングのための過渡温度スペクトルを読み込み、検査します。
詳細説明:
    読み込み処理はプロット処理から意図的に分離されています。
    これにより、将来的に別の機器やエクスポート形式をサポートする際に、
    新しい読み込み関数として追加しやすくなっています。
"""
from __future__ import annotations

import argparse
import csv
import re
import warnings
from pathlib import Path
from typing import Any

import matplotlib.pyplot as plt
import numpy as np


KAPPA_RE = re.compile(r"^\s*(?:k|kappa|κ)\s*=\s*([-+0-9.eE]+)\s*$", re.I)


def _read_k_sweep_csv(path: str | Path) -> tuple[list[np.ndarray], list[np.ndarray], dict[str, Any]]:
    """
    概要:
        横持ちCSV形式を読み込みます。
    詳細説明:
        time_s, experiment, k=値 などの列を含むワイドCSV形式に対応しています。
    引数:
        :param path: 読み込むファイルのパス
        :type path: str | Path
    戻り値:
        :returns: 時間配列のリスト、信号配列のリスト、メタデータの辞書のタプル
        :rtype: tuple[list[np.ndarray], list[np.ndarray], dict[str, Any]]
    例外:
        :raises ValueError: ヘッダーやデータ行が不足している場合に発生
    """
    path = Path(path)
    with path.open("r", encoding="utf-8-sig", newline="") as fh:
        # Sniffer permits comma, tab, or semicolon exports while retaining CSV quoting.
        sample = fh.read(8192)
        fh.seek(0)
        try:
            dialect = csv.Sniffer().sniff(sample, delimiters=",\t;")
        except csv.Error:
            dialect = csv.excel
        rows = list(csv.reader(fh, dialect))

    if len(rows) < 2:
        raise ValueError(f"{path}: header and at least one numeric data row are required.")
    header = [cell.strip() for cell in rows[0]]
    if len(header) < 2:
        raise ValueError(f"{path}: at least a time column and one spectrum are required.")

    values = np.full((len(rows) - 1, len(header)), np.nan, dtype=float)
    for i, row in enumerate(rows[1:]):
        for j, cell in enumerate(row[: len(header)]):
            try:
                values[i, j] = float(cell.strip())
            except (ValueError, AttributeError):
                pass
    good_time = np.isfinite(values[:, 0])
    if not np.any(good_time):
        raise ValueError(f"{path}: no numeric values were found in the first (time) column.")

    x = values[good_time, 0]
    y_list = [values[good_time, j] for j in range(1, len(header))]
    labels = header[1:]
    kappa = []
    experiment_indices = []
    for i, label in enumerate(labels):
        match = KAPPA_RE.match(label)
        kappa.append(float(match.group(1)) if match else None)
        if label.lower() in {"experiment", "experimental", "measured", "measurement"}:
            experiment_indices.append(i)

    inf: dict[str, Any] = {
        "filename": str(path),
        "format": "wide k-sweep CSV",
        "x_label": header[0],
        "y_labels": labels,
        "kappa": kappa,
        "experiment_indices": experiment_indices,
    }
    # x_list is intentionally a list: future formats may have a different time grid per spectrum.
    return [x.copy() for _ in y_list], y_list, inf


def read_data_files(filename: str | Path) -> tuple[list[np.ndarray], list[np.ndarray], dict[str, Any]]:
    """
    概要:
        サポートされているデータファイルからデータとメタデータを返します。
    詳細説明:
        別の読み込み関数を定義し、ここで選択することでフォーマットを追加できます。
        現在の読み込み処理は、timeという名前の最初の列と、k=数値またはkappa=数値
        という名前の計算列を検出します。
    引数:
        :param filename: 入力ファイル名
        :type filename: str | Path
    戻り値:
        :returns: x_list、y_list、メタデータinfのタプル
        :rtype: tuple[list[np.ndarray], list[np.ndarray], dict[str, Any]]
    例外:
        :raises ValueError: サポートされていない入力形式の場合に発生
    """
    path = Path(filename)
    with path.open("r", encoding="utf-8-sig", newline="") as fh:
        first_line = fh.readline().lower()
    if "time" in first_line and ("k=" in first_line or "kappa=" in first_line or "κ=" in first_line):
        return _read_k_sweep_csv(path)
    raise ValueError(
        f"Unsupported input format: {path}. Add its reader to read_data_files()."
    )


def _plot_split_time_axis(ax_linear, ax_log, x_list, y_list, inf: dict[str, Any]) -> None:
    """
    概要:
        隣り合う軸に t <= 0 を線形で、t > 0 を対数でプロットします。
    引数:
        :param ax_linear: 線形スケールの軸
        :type ax_linear: matplotlib.axes.Axes
        :param ax_log: 対数スケールの軸
        :type ax_log: matplotlib.axes.Axes
        :param x_list: xデータ配列のリスト
        :type x_list: list[np.ndarray]
        :param y_list: yデータ配列のリスト
        :type y_list: list[np.ndarray]
        :param inf: メタデータの辞書
        :type inf: dict[str, Any]
    例外:
        :raises ValueError: 対数軸にするための正の時間値が存在しない場合に発生
    """
    labels = inf["y_labels"]
    kappas = inf["kappa"]
    experiment = set(inf["experiment_indices"])
    positive = np.concatenate([x[np.isfinite(x) & (x > 0)] for x in x_list])
    if positive.size == 0:
        raise ValueError("xscale='log' requires at least one positive time value.")
    tmin = positive.min()

    for i, (x, y) in enumerate(zip(x_list, y_list)):
        style = {"lw": 1.25, "alpha": 0.8}
        if i in experiment:
            style.update(color="black", lw=1.0, alpha=0.65, label=labels[i], zorder=3)
        else:
            style["label"] = labels[i]
        neg = x <= 0
        pos = x > 0
        ax_linear.plot(x[neg], y[neg], **style)
        ax_log.plot(x[pos], y[pos], **style)

    ax_log.set_xscale("log")
    ax_linear.set_xlim(min(np.nanmin(x) for x in x_list), 0.0)
    ax_log.set_xlim(tmin, max(np.nanmax(x) for x in x_list))
    ax_linear.spines["right"].set_visible(False)
    ax_log.spines["left"].set_visible(False)
    ax_log.tick_params(left=False, labelleft=False)
    ax_linear.grid(True, alpha=0.25)
    ax_log.grid(True, which="both", alpha=0.25)
    # small diagonal marks make the scale change explicit
    d = 0.012
    kwargs = dict(color="k", clip_on=False, lw=0.8)
    ax_linear.plot((1 - d, 1 + d), (-d, +d), transform=ax_linear.transAxes, **kwargs)
    ax_linear.plot((1 - d, 1 + d), (1 - d, 1 + d), transform=ax_linear.transAxes, **kwargs)
    ax_log.plot((-d, +d), (-d, +d), transform=ax_log.transAxes, **kwargs)
    ax_log.plot((-d, +d), (1 - d, 1 + d), transform=ax_log.transAxes, **kwargs)


def plot_data(x_list, y_list, inf: dict[str, Any], xscale: str = "linear"):
    """
    概要:
        測定されたスペクトルとkスイープスペクトルを可視化します。
    引数:
        :param x_list: xデータ配列のリスト
        :type x_list: list[np.ndarray]
        :param y_list: yデータ配列のリスト
        :type y_list: list[np.ndarray]
        :param inf: メタデータの辞書
        :type inf: dict[str, Any]
        :param xscale: x軸のスケール。linearまたはlogを指定
        :type xscale: str
    戻り値:
        :returns: 作成されたMatplotlibのFigureオブジェクト
        :rtype: matplotlib.figure.Figure
    例外:
        :raises ValueError: xscaleが無効な値の場合に発生
    """
    if xscale not in {"linear", "log"}:
        raise ValueError("xscale must be 'linear' or 'log'.")
    if xscale == "log":
        fig, (ax0, ax1) = plt.subplots(1, 2, sharey=True, figsize=(11, 6),
                                       gridspec_kw={"width_ratios": [1.0, 2.1], "wspace": 0.05})
        _plot_split_time_axis(ax0, ax1, x_list, y_list, inf)
        ax0.set_xlabel("time (s), linear: t ≤ 0")
        ax1.set_xlabel("time (s), logarithmic: t > 0")
        axes = (ax0, ax1)
    else:
        fig, ax = plt.subplots(figsize=(10, 6))
        for i, (x, y) in enumerate(zip(x_list, y_list)):
            style = {"lw": 1.25, "alpha": 0.8, "label": inf["y_labels"][i]}
            if i in inf["experiment_indices"]:
                style.update(color="black", lw=1.0, alpha=0.65, zorder=3)
            ax.plot(x, y, **style)
        ax.grid(True, alpha=0.25)
        ax.set_xlabel("time (s)")
        axes = (ax,)

    axes[0].set_ylabel("temperature signal")
    axes[-1].legend(loc="best", fontsize=8, ncol=2)
    fig.suptitle(Path(inf["filename"]).name)
    fig.tight_layout()
    return fig


def smooth_series(x: np.ndarray, y: np.ndarray, nsmooth: int, method: str) -> np.ndarray:
    """
    概要:
        中心移動平均または局所3次多項式を使用して系列を平滑化します。
    引数:
        :param x: xデータ配列
        :type x: np.ndarray
        :param y: yデータ配列
        :type y: np.ndarray
        :param nsmooth: 平滑化に使用するデータ点数
        :type nsmooth: int
        :param method: 平滑化手法。meanまたはpoly3を指定
        :type method: str
    戻り値:
        :returns: 平滑化されたyデータ配列
        :rtype: np.ndarray
    例外:
        :raises ValueError: nsmoothが偶数、またはデータ点数を超えている場合などに発生
    """
    if nsmooth <= 1:
        return y.copy()
    if nsmooth % 2 == 0:
        raise ValueError("nsmooth must be odd so that smoothing is centred at each data point.")
    if nsmooth > len(y):
        raise ValueError("nsmooth cannot exceed the number of data points.")
    if method not in {"mean", "poly3"}:
        raise ValueError("smooth method must be 'mean' or 'poly3'.")
    half = nsmooth // 2
    out = np.empty_like(y, dtype=float)
    for i in range(len(y)):
        start = max(0, min(i - half, len(y) - nsmooth))
        stop = start + nsmooth
        xx, yy = x[start:stop], y[start:stop]
        valid = np.isfinite(xx) & np.isfinite(yy)
        if method == "mean":
            out[i] = np.mean(yy[valid]) if np.any(valid) else np.nan
        else:
            degree = min(3, int(np.count_nonzero(valid)) - 1)
            out[i] = (np.polyval(np.polyfit(xx[valid] - x[i], yy[valid], degree), 0.0)
                      if degree >= 0 else np.nan)
    return out


def plot_difference(x_list, y_list, inf: dict[str, Any], xscale: str = "linear",
                    nsmooth: int = 1, smooth: str = "poly3"):
    """
    概要:
        利用可能なすべてのkappaについて残差 y_experiment - y_simulation をプロットします。
    引数:
        :param x_list: xデータ配列のリスト
        :type x_list: list[np.ndarray]
        :param y_list: yデータ配列のリスト
        :type y_list: list[np.ndarray]
        :param inf: メタデータの辞書
        :type inf: dict[str, Any]
        :param xscale: x軸のスケール
        :type xscale: str
        :param nsmooth: 平滑化に使用する点数
        :type nsmooth: int
        :param smooth: 平滑化手法
        :type smooth: str
    戻り値:
        :returns: 作成されたMatplotlibのFigureオブジェクト
        :rtype: matplotlib.figure.Figure
    例外:
        :raises ValueError: 実験列が1つではない場合や対数スケールに正の時間がない場合に発生
    """
    exp_indices = inf["experiment_indices"]
    if len(exp_indices) != 1:
        raise ValueError("Exactly one experimental column is required for a difference plot.")
    i_exp = exp_indices[0]
    y_exp = np.asarray(y_list[i_exp], dtype=float)

    if xscale == "log":
        fig, (ax0, ax1) = plt.subplots(1, 2, sharey=True, figsize=(11, 6),
                                       gridspec_kw={"width_ratios": [1.0, 2.1], "wspace": 0.05})
        positive = np.concatenate([x[np.isfinite(x) & (x > 0)] for x in x_list])
        if positive.size == 0:
            raise ValueError("xscale='log' requires at least one positive time value.")
        tmin = positive.min()
        for i, (x, y, k) in enumerate(zip(x_list, y_list, inf["kappa"])):
            if k is None:
                continue
            residual = smooth_series(x, y_exp - np.asarray(y, dtype=float), nsmooth, smooth)
            ax0.plot(x[x <= 0], residual[x <= 0], lw=1.1, alpha=0.8, label=f"κ={k:g}")
            ax1.plot(x[x > 0], residual[x > 0], lw=1.1, alpha=0.8, label=f"κ={k:g}")
        ax1.set_xscale("log")
        ax0.set_xlim(min(np.nanmin(x) for x in x_list), 0.0)
        ax1.set_xlim(tmin, max(np.nanmax(x) for x in x_list))
        ax0.spines["right"].set_visible(False)
        ax1.spines["left"].set_visible(False)
        ax1.tick_params(left=False, labelleft=False)
        ax0.set_xlabel("time (s), linear: t ≤ 0")
        ax1.set_xlabel("time (s), logarithmic: t > 0")
        axes = (ax0, ax1)
    else:
        fig, ax = plt.subplots(figsize=(10, 6))
        for x, y, k in zip(x_list, y_list, inf["kappa"]):
            if k is not None:
                residual = smooth_series(x, y_exp - np.asarray(y, dtype=float), nsmooth, smooth)
                ax.plot(x, residual, lw=1.1, alpha=0.8, label=f"κ={k:g}")
        ax.set_xlabel("time (s)")
        axes = (ax,)

    for ax in axes:
        ax.axhline(0.0, color="black", lw=0.8)
        ax.grid(True, alpha=0.25)
    axes[0].set_ylabel("experiment − simulation")
    axes[-1].legend(loc="best", fontsize=8, ncol=2)
    smooth_label = "unsmoothed" if nsmooth <= 1 else f"{smooth}, n={nsmooth}"
    fig.suptitle(f"Residual spectra ({smooth_label}): {Path(inf['filename']).name}")
    fig.tight_layout()
    return fig


def calculate_sse(
    x_list, y_list, inf: dict[str, Any], sse_kind: str = "raw",
    tmin: float | None = None, tmax: float | None = None, rel_error: float = 0.0,
) -> tuple[np.ndarray, np.ndarray, int]:
    """
    概要:
        実験スペクトルと各シミュレーション間のkappa値とSSEを返します。
    詳細説明:
        sse_kindがrawの場合は通常のSSEを使用します。sse_kindがlogの場合は
        残差をy_scale(t)で割ったものの二乗和を使用します。正のy_expに対しては
        y_scale=y_expとし、これは対数信号残差の1次の形になります。
        正でない値では、y_scaleはそれまでに観測された最小の正の実験値となり、
        不安定な対数変換を回避します。rel_errorを指定した場合、raw SSEに限り
        残差は 1 / max(s_noise, rel_error * abs(y_exp)) で重み付けされます。
        ここでs_noiseは選択された時間範囲内の最良の非重み付けフィットから推定されます。
    引数:
        :param x_list: xデータ配列のリスト
        :type x_list: list[np.ndarray]
        :param y_list: yデータ配列のリスト
        :type y_list: list[np.ndarray]
        :param inf: メタデータの辞書
        :type inf: dict[str, Any]
        :param sse_kind: SSEの種類。rawまたはlogを指定
        :type sse_kind: str
        :param tmin: 対象とする最小時間
        :type tmin: float | None
        :param tmax: 対象とする最大時間
        :type tmax: float | None
        :param rel_error: 相対誤差の重み付け係数
        :type rel_error: float
    戻り値:
        :returns: kappaの配列、SSEの配列、有効データ点数のタプル
        :rtype: tuple[np.ndarray, np.ndarray, int]
    例外:
        :raises ValueError: 引数の組み合わせが不正な場合や対象データがない場合に発生
    """
    exp_indices = inf["experiment_indices"]
    if len(exp_indices) != 1:
        raise ValueError("Exactly one experimental column is required for SSE calculation.")
    y_exp = np.asarray(y_list[exp_indices[0]], dtype=float)
    time = np.asarray(x_list[exp_indices[0]], dtype=float)
    if sse_kind not in {"raw", "log"}:
        raise ValueError("sse_kind must be 'raw' or 'log'.")
    if rel_error < 0:
        raise ValueError("rel_error must be non-negative.")
    if sse_kind == "log" and rel_error:
        raise ValueError("--rel-error cannot be combined with sse-log; both specify relative weighting.")
    time_mask = np.isfinite(time)
    if tmin is not None:
        time_mask &= time >= tmin
    if tmax is not None:
        time_mask &= time <= tmax
    if not np.any(time_mask):
        raise ValueError("The selected --tmin/--tmax range contains no data points.")
    if sse_kind == "log":
        positive = np.where(y_exp > 0, y_exp, np.inf)
        running_floor = np.minimum.accumulate(positive)
        if not np.isfinite(running_floor).any():
            raise ValueError("log SSE requires at least one positive experimental signal.")
        # If the first samples are non-positive, use the first available floor.
        first_floor = running_floor[np.isfinite(running_floor)][0]
        scale = np.where(np.isfinite(running_floor), running_floor, first_floor)
        scale = np.where(y_exp > 0, y_exp, scale)
    kappa, raw_sse, residuals = [], [], []
    n_points: int | None = None
    for y_sim, k in zip(y_list, inf["kappa"]):
        if k is None:
            continue
        valid = time_mask & np.isfinite(y_exp) & np.isfinite(y_sim)
        if not np.any(valid):
            continue
        if n_points is None:
            n_points = int(np.count_nonzero(valid))
        elif n_points != int(np.count_nonzero(valid)):
            raise ValueError("Simulation columns have different valid-point counts; add explicit weighting first.")
        residual = y_exp[valid] - np.asarray(y_sim)[valid]
        if sse_kind == "log":
            residual = residual / scale[valid]
        kappa.append(k)
        raw_sse.append(np.sum(residual**2))
        residuals.append(residual)
    if not kappa:
        raise ValueError("No simulation columns named k=<number> were found.")
    kappa = np.asarray(kappa)
    raw_sse = np.asarray(raw_sse)
    if rel_error:
        i_best = int(np.argmin(raw_sse))
        noise_floor = np.sqrt(raw_sse[i_best] / max((n_points or 1) - 1, 1))
        # This max() form prevents the low-signal tail from receiving unlimited weight.
        sigma = np.maximum(noise_floor, rel_error * np.abs(y_exp))
        sse = []
        # Rebuild directly from the κ columns; this keeps masks identical to raw SSE.
        for j, k in enumerate(inf["kappa"]):
            if k is None:
                continue
            y_sim = np.asarray(y_list[j], dtype=float)
            valid = time_mask & np.isfinite(y_exp) & np.isfinite(y_sim)
            sse.append(np.sum(((y_exp[valid] - y_sim[valid]) / sigma[valid]) ** 2))
        sse = np.asarray(sse)
    else:
        sse = raw_sse
    order = np.argsort(kappa)
    return kappa[order], sse[order], n_points or 0


def plot_sse(kappa: np.ndarray, sse: np.ndarray, sse_kind: str = "raw"):
    """
    概要:
        kappaに対するSSEの依存性をプロットします。
    引数:
        :param kappa: kappa値の配列
        :type kappa: np.ndarray
        :param sse: SSEの配列
        :type sse: np.ndarray
        :param sse_kind: SSEの種類
        :type sse_kind: str
    戻り値:
        :returns: 作成されたMatplotlibのFigureオブジェクト
        :rtype: matplotlib.figure.Figure
    """
    fig, ax = plt.subplots(figsize=(8, 5))
    ax.plot(kappa, sse, "o-", color="tab:blue")
    best = np.argmin(sse)
    ax.axvline(kappa[best], color="tab:red", ls="--", lw=1, label=f"minimum: κ={kappa[best]:g}")
    ylabel = "log-weighted SSE" if sse_kind == "log" else "SSE"
    ax.set(xlabel="thermal conductivity, κ", ylabel=ylabel, title="SSE of experiment vs. simulation")
    ax.grid(True, alpha=0.3)
    ax.legend()
    fig.tight_layout()
    return fig


def likelihood_kernel_from_sse(sse: np.ndarray, noise_std: float) -> tuple[np.ndarray, np.ndarray]:
    """
    概要:
        delta SSEと相対的なガウス誤差尤度カーネルを返します。
    詳細説明:
        共通の点ごとの標準偏差 noise_std を持つ独立した正規残差について、
        相対カーネルは exp(-(SSE - SSE_min) / (2 * noise_std * noise_std)) となります。
    引数:
        :param sse: SSEの配列
        :type sse: np.ndarray
        :param noise_std: ノイズの標準偏差
        :type noise_std: float
    戻り値:
        :returns: delta SSEの配列と尤度カーネルの配列のタプル
        :rtype: tuple[np.ndarray, np.ndarray]
    例外:
        :raises ValueError: noise_stdが正ではない場合に発生
    """
    if noise_std <= 0:
        raise ValueError("noise_std must be positive.")
    delta_chi2 = (sse - np.min(sse)) / noise_std**2
    return delta_chi2, np.exp(-0.5 * delta_chi2)


def estimate_kappa_std(kappa: np.ndarray, delta_sse: np.ndarray) -> tuple[float, float, float]:
    """
    概要:
        delta SSEが1となるkappaを見つけ、k_best、sigma_left、sigma_rightを報告します。
    詳細説明:
        delta SSEが1になることは、カーネルが exp(-1/2) に等しいことと同義です。
        最小値の両側で、隣接する計算済みkappa点の間に線形補間が使用されます。
        交差が見つからない場合はNaNが返されます。
    引数:
        :param kappa: kappa値の配列
        :type kappa: np.ndarray
        :param delta_sse: delta SSEの配列
        :type delta_sse: np.ndarray
    戻り値:
        :returns: 最良のkappa、左側の不確かさ、右側の不確かさのタプル
        :rtype: tuple[float, float, float]
    """
    i0 = int(np.argmin(delta_sse))
    k0 = float(kappa[i0])

    def crossing(indices: np.ndarray) -> float:
        for ia, ib in zip(indices[:-1], indices[1:]):
            ya, yb = delta_sse[ia], delta_sse[ib]
            if (ya - 1.0) * (yb - 1.0) <= 0 and ya != yb:
                return float(kappa[ia] + (1.0 - ya) * (kappa[ib] - kappa[ia]) / (yb - ya))
        return float("nan")

    # On the left enumerate from far-left toward the minimum; right is natural.
    left_cross = crossing(np.arange(0, i0 + 1)) if i0 else float("nan")
    right_cross = crossing(np.arange(i0, len(kappa))) if i0 < len(kappa) - 1 else float("nan")
    return k0, k0 - left_cross, right_cross - k0


def likelihood_widths_from_reference(kappa: np.ndarray, sse: np.ndarray, k_ref: float,
                                     sse_ref: float, sse_scale: float) -> tuple[float, float]:
    """
    概要:
        連続的な最適値の周りにおける、線形補間されたSSE尤度の2つの交差幅を見つけます。
    引数:
        :param kappa: kappa値の配列
        :type kappa: np.ndarray
        :param sse: SSEの配列
        :type sse: np.ndarray
        :param k_ref: 基準となるkappa値
        :type k_ref: float
        :param sse_ref: 基準となるSSE値
        :type sse_ref: float
        :param sse_scale: 交差を判定するためのSSEのスケール
        :type sse_scale: float
    戻り値:
        :returns: 左側の不確かさと右側の不確かさのタプル
        :rtype: tuple[float, float]
    """
    threshold = sse_ref + sse_scale
    crossings = []
    for ka, kb, sa, sb in zip(kappa[:-1], kappa[1:], sse[:-1], sse[1:]):
        if (sa - threshold) * (sb - threshold) <= 0 and sa != sb:
            crossings.append(float(ka + (threshold - sa) * (kb - ka) / (sb - sa)))
    left = [v for v in crossings if v < k_ref]
    right = [v for v in crossings if v > k_ref]
    return (k_ref - max(left) if left else float("nan"),
            min(right) - k_ref if right else float("nan"))


def fit_local_sse_polynomial(kappa: np.ndarray, sse: np.ndarray) -> dict[str, Any]:
    """
    概要:
        局所的なSSEの最小値をフィッティングします。
    詳細説明:
        中心の5点を用いて3次多項式でフィッティングしますが、点が少ない場合は2次多項式を
        使用します。フィッティングは意図的に局所的になっています。大域的な多項式を使うと、
        最小値から遠く離れて急増するSSEが、最適値付近の形状に不自然な影響を与えるためです。
    引数:
        :param kappa: kappa値の配列
        :type kappa: np.ndarray
        :param sse: SSEの配列
        :type sse: np.ndarray
    戻り値:
        :returns: 多項式フィッティングの結果と状態を含む辞書
        :rtype: dict[str, Any]
    """
    i0 = int(np.argmin(sse))
    n = len(kappa)
    if i0 >= 2 and i0 <= n - 3:
        indices = np.arange(i0 - 2, i0 + 3)
        degree = 3
        message = None
    elif i0 >= 1 and i0 <= n - 2:
        indices = np.arange(i0 - 1, i0 + 2)
        degree = 2
        message = "Only three bracketing κ points are available; quadratic interpolation used."
        warnings.warn(message, RuntimeWarning, stacklevel=2)
    else:
        message = "Minimum is at the edge of the κ sweep; polynomial interpolation is unavailable."
        warnings.warn(message, RuntimeWarning, stacklevel=2)
        return {"ok": False, "message": message, "kappa": float(kappa[i0])}

    k_ref = float(kappa[i0])
    x = kappa[indices] - k_ref
    coeff = np.polyfit(x, sse[indices], degree)
    stationary = np.roots(np.polyder(coeff))
    xmin, xmax = x.min(), x.max()
    candidates = []
    for root in stationary:
        if abs(root.imag) < 1e-8 and xmin <= root.real <= xmax:
            curvature = np.polyval(np.polyder(coeff, 2), root.real)
            if curvature > 0:
                candidates.append(float(root.real))
    if not candidates:
        message = "Local polynomial has no valid minimum inside its fitting range; using discrete minimum."
        warnings.warn(message, RuntimeWarning, stacklevel=2)
        return {"ok": False, "message": message, "kappa": float(kappa[i0])}
    x_best = min(candidates, key=abs)
    return {
        "ok": True, "degree": degree, "indices": indices, "coeff": coeff,
        "k_ref": k_ref, "kappa": k_ref + x_best,
        "sse_min": float(np.polyval(coeff, x_best)), "x_best": x_best,
        "x_limits": (float(xmin), float(xmax)),
        "message": message,
    }


def polynomial_sigma(polyfit: dict[str, Any], sse_scale: float) -> tuple[float, float]:
    """
    概要:
        フィットされたkappaの範囲内で、delta SSEがsse_scaleとなる多項式SSEの幅を見つけます。
    引数:
        :param polyfit: 多項式フィッティングの結果を格納した辞書
        :type polyfit: dict[str, Any]
        :param sse_scale: 交差を判定するためのSSEのスケール
        :type sse_scale: float
    戻り値:
        :returns: 左側の不確かさと右側の不確かさのタプル
        :rtype: tuple[float, float]
    """
    if not polyfit.get("ok"):
        return float("nan"), float("nan")
    coeff = polyfit["coeff"].copy()
    coeff[-1] -= polyfit["sse_min"] + sse_scale
    roots = np.roots(coeff)
    k_best = polyfit["kappa"]
    xmin, xmax = polyfit["x_limits"]
    k_left_limit = polyfit["k_ref"] + xmin
    k_right_limit = polyfit["k_ref"] + xmax
    real_k = [polyfit["k_ref"] + r.real for r in roots if abs(r.imag) < 1e-8]
    tolerance = 1e-10 * max(1.0, abs(k_best))
    left = [v for v in real_k if k_left_limit - tolerance <= v < k_best]
    right = [v for v in real_k if k_best < v <= k_right_limit + tolerance]
    return (k_best - max(left) if left else float("nan"),
            min(right) - k_best if right else float("nan"))


def plot_likelihood_kernel(kappa: np.ndarray, delta_chi2: np.ndarray, k0: float, sigma_left: float,
                           sigma_right: float, polyfit: dict[str, Any] | None = None,
                           sse_scale: float | None = None,
                           poly_sigma_left: float = float("nan"), poly_sigma_right: float = float("nan")):
    """
    概要:
        尤度カーネルをプロットして可視化します。
    引数:
        :param kappa: kappa値の配列
        :type kappa: np.ndarray
        :param delta_chi2: delta chi2の配列
        :type delta_chi2: np.ndarray
        :param k0: 最良のkappa値
        :type k0: float
        :param sigma_left: 左側の不確かさ
        :type sigma_left: float
        :param sigma_right: 右側の不確かさ
        :type sigma_right: float
        :param polyfit: 多項式フィッティングの結果の辞書
        :type polyfit: dict[str, Any] | None
        :param sse_scale: SSEのスケール
        :type sse_scale: float | None
        :param poly_sigma_left: 多項式による左側の不確かさ
        :type poly_sigma_left: float
        :param poly_sigma_right: 多項式による右側の不確かさ
        :type poly_sigma_right: float
    戻り値:
        :returns: 作成されたMatplotlibのFigureオブジェクト
        :rtype: matplotlib.figure.Figure
    """
    fig, ax = plt.subplots(figsize=(8, 5))
    level = np.exp(-0.5)
    # Interpolate SSE (equivalently Δχ²), then apply the exponential.  Thus the
    # displayed curve and the Δχ²=1 / exp(-1/2) uncertainty criterion agree.
    k_curve, delta_curve = kappa, delta_chi2
    if polyfit and polyfit.get("ok"):
        # Insert the continuous regression minimum into the likelihood curve.
        # This avoids treating the closest sampled κ as its maximum.
        k_curve = np.append(kappa, polyfit["kappa"])
        delta_curve = np.append(delta_chi2, 0.0)
        order = np.argsort(k_curve)
        k_curve, delta_curve = k_curve[order], delta_curve[order]
    k_dense = np.linspace(kappa.min(), kappa.max(), 2000)
    kernel_dense = np.exp(-0.5 * np.interp(k_dense, k_curve, delta_curve))
    kernel = np.exp(-0.5 * delta_chi2)
    ax.plot(k_dense, kernel_dense, "-", color="tab:purple", label="relative likelihood kernel")
    ax.plot(kappa, kernel, "o", color="tab:purple", label="calculated κ points")
    ax.axhline(level, color="0.35", ls=":", label="exp(−1/2)")
    ax.axvline(k0, color="tab:red", ls="--", lw=1, label=f"best κ={k0:g}")
    if np.isfinite(sigma_left):
        ax.axvline(k0 - sigma_left, color="tab:green", ls="--", lw=1)
    if np.isfinite(sigma_right):
        ax.axvline(k0 + sigma_right, color="tab:green", ls="--", lw=1)
    if polyfit and polyfit.get("ok") and sse_scale is not None:
        xmin, xmax = polyfit["x_limits"]
        kp = np.linspace(polyfit["k_ref"] + xmin, polyfit["k_ref"] + xmax, 400)
        sse_p = np.polyval(polyfit["coeff"], kp - polyfit["k_ref"])
        kernel_p = np.exp(-(sse_p - polyfit["sse_min"]) / (2 * sse_scale))
        ax.plot(kp, kernel_p, color="tab:orange", lw=1.7,
                label=f"local SSE polynomial (degree {polyfit['degree']})")
        ax.axvline(polyfit["kappa"], color="tab:orange", ls="--", lw=1,
                   label=f"polynomial best κ={polyfit['kappa']:.4g}")
        if np.isfinite(poly_sigma_left):
            ax.axvline(polyfit["kappa"] - poly_sigma_left, color="tab:orange", ls=":", lw=1)
        if np.isfinite(poly_sigma_right):
            ax.axvline(polyfit["kappa"] + poly_sigma_right, color="tab:orange", ls=":", lw=1)
    ax.set(xlabel="thermal conductivity, κ", ylabel="relative likelihood", ylim=(-0.03, 1.05),
           title="Likelihood kernel from SSE")
    ax.grid(True, alpha=0.3)
    ax.legend()
    fig.tight_layout()
    return fig


def main() -> None:
    """
    概要:
        コマンドライン引数を解析し、データの読み込み、SSE計算、およびプロットを実行します。
    """
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("input", help="CSV containing experiment and k-sweep spectra")
    parser.add_argument("--mode", choices=["read", "diff", "sse", "std", "sse-log", "std-log"], default="read",
                        help="read: spectra; diff: experiment−simulation; sse/std: ordinary SSE; sse-log/std-log: log-weighted SSE")
    parser.add_argument("--xscale", choices=["linear", "log"], default="linear")
    parser.add_argument("--tmin", type=float, help="first time included in SSE/likelihood calculation (s)")
    parser.add_argument("--tmax", type=float, help="last time included in SSE/likelihood calculation (s)")
    parser.add_argument("--rel-error", type=float, default=0.0, metavar="A",
                        help="relative error a: use sigma(t)=max(noise_floor, a*abs(y_exp)); raw SSE modes only")
    parser.add_argument("--nsmooth", type=int, default=1,
                        help="odd number of points for mode=diff smoothing; 1 disables smoothing")
    parser.add_argument("--smooth", choices=["poly3", "mean"], default="poly3",
                        help="mode=diff smoothing: local cubic polynomial (poly3) or moving mean")
    parser.add_argument("--output", help="save plot (e.g. spectra.png); omit to open a window")
    args = parser.parse_args()
    x_list, y_list, inf = read_data_files(args.input)
    if args.mode == "read":
        fig = plot_data(x_list, y_list, inf, args.xscale)
    elif args.mode == "diff":
        fig = plot_difference(x_list, y_list, inf, args.xscale, args.nsmooth, args.smooth)
    else:
        sse_kind = "log" if args.mode.endswith("-log") else "raw"
        kappa, sse, n_points = calculate_sse(
            x_list, y_list, inf, sse_kind,
            tmin=args.tmin, tmax=args.tmax, rel_error=args.rel_error,
        )
        if args.mode.startswith("sse"):
            fig = plot_sse(kappa, sse, sse_kind)
        else:
            # The residual standard deviation is estimated at the best κ.
            # One fitted parameter gives N−1 residual degrees of freedom.
            noise_std = np.sqrt(np.min(sse) / max(n_points - 1, 1))
            polyfit = fit_local_sse_polynomial(kappa, sse)
            if polyfit.get("ok"):
                # The regression supplies the continuous likelihood maximum.
                k0 = polyfit["kappa"]
                sse_ref = polyfit["sse_min"]
            else:
                k0 = float(kappa[np.argmin(sse)])
                sse_ref = float(np.min(sse))
            delta_chi2 = (sse - sse_ref) / noise_std**2
            kernel = np.exp(-0.5 * delta_chi2)
            sigma_left, sigma_right = likelihood_widths_from_reference(
                kappa, sse, k0, sse_ref, noise_std**2,
            )
            poly_sigma_left, poly_sigma_right = polynomial_sigma(polyfit, noise_std**2)
            if polyfit.get("ok") and (not np.isfinite(poly_sigma_left) or not np.isfinite(poly_sigma_right)):
                warnings.warn(
                    "The local polynomial fitting range does not contain both exp(-1/2) crossings; "
                    "use the likelihood-interpolation widths for uncertainty.",
                    RuntimeWarning,
                    stacklevel=1,
                )
            fig = plot_likelihood_kernel(
                kappa, delta_chi2, k0, sigma_left, sigma_right,
                polyfit=polyfit, sse_scale=noise_std**2,
                poly_sigma_left=poly_sigma_left, poly_sigma_right=poly_sigma_right,
            )
            print(f"best kappa = {k0:.8g}")
            print(f"estimated pointwise residual standard deviation = {noise_std:.8g}")
            print(f"s_eff from likelihood (kernel=exp(-1/2)): -{sigma_left:.8g}, +{sigma_right:.8g}")
            if polyfit.get("ok"):
                print(f"polynomial degree = {polyfit['degree']}  (local κ range: "
                      f"{polyfit['k_ref'] + polyfit['x_limits'][0]:.8g} to "
                      f"{polyfit['k_ref'] + polyfit['x_limits'][1]:.8g})")
                print(f"s_eff from polynomial (kernel=exp(-1/2)): "
                      f"-{poly_sigma_left:.8g}, +{poly_sigma_right:.8g}")
            else:
                print(f"WARNING: {polyfit['message']}")
    if args.output:
        fig.savefig(args.output, dpi=180, bbox_inches="tight")
    else:
        plt.show()


if __name__ == "__main__":
    main()