#!/usr/bin/env python3
"""
FrESH temperature characterization - step analysis
Final paper-facing version.

Default behavior:
- Generates only the paper figure for the static step analysis
- Writes Table 1 and Appendix Table A1 summaries
- Prints manuscript-ready values to the console
- Saves a compact text/JSON report with the values used in the paper

Expected repository layout:
  .
  ├── analyze_step.py
  ├── data/
  │   ├── step_windows.json
  │   ├── RE1050_STEP_<run>_pt100.csv
  │   ├── RE1050_STEP_<run>_chiller.csv
  │   ├── RP1845_STEP_<run>_pt100.csv
  │   └── RP1845_STEP_<run>_chiller.csv
  └── output/
"""

from __future__ import annotations

import json
import re
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

CHILLERS = ["RE1050", "RP1845"]
SENSORS = ["Sensor1", "Sensor2", "Sensor3", "Sensor4", "Sensor5"]

CHILLER_COLORS = {
    "RE1050": "#2E6E4D",
    "RP1845": "#D2874C",
}

BASE_DIR = Path(__file__).resolve().parent
DATA_DIR = BASE_DIR / "data"
OUTPUT_DIR = BASE_DIR / "output"
STEP_WINDOWS_FILE = DATA_DIR / "step_windows.json"

SAVE_SUPPORTING_CSV = True
SAVE_LATEX_TABLES = True

PAIR_RE = re.compile(
    r"(?P<chiller>RE1050|RP1845)_STEP_(?P<run>.+)_(?P<dtype>pt100|chiller)\.csv$",
    re.IGNORECASE,
)


def ensure_output_dir() -> None:
    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)


def load_step_windows() -> dict:
    if not STEP_WINDOWS_FILE.exists():
        raise FileNotFoundError(f"Missing step windows file: {STEP_WINDOWS_FILE}")
    with open(STEP_WINDOWS_FILE, "r", encoding="utf-8") as f:
        return json.load(f)


def to_time_seconds(series: pd.Series) -> pd.Series:
    s = series.astype(str).str.strip()

    numeric = pd.to_numeric(s, errors="coerce")
    if numeric.notna().mean() > 0.8:
        return numeric.astype(float)

    t_only = pd.to_datetime(s, format="%H:%M:%S", errors="coerce")
    if t_only.notna().mean() > 0.8:
        return (
            t_only.dt.hour * 3600
            + t_only.dt.minute * 60
            + t_only.dt.second
        ).astype(float)

    dt = pd.to_datetime(s, errors="coerce")
    if dt.notna().mean() > 0.8:
        return (
            dt.dt.hour * 3600
            + dt.dt.minute * 60
            + dt.dt.second
            + dt.dt.microsecond / 1e6
        ).astype(float)

    td = pd.to_timedelta(s, errors="coerce")
    if td.notna().mean() > 0.8:
        return td.dt.total_seconds().astype(float)

    raise ValueError("Could not interpret time column as datetime, HH:MM:SS, or elapsed seconds.")


def parse_window_to_seconds(window_spec) -> tuple[float, float]:
    if isinstance(window_spec, dict):
        if "start_s" in window_spec and "end_s" in window_spec:
            return float(window_spec["start_s"]), float(window_spec["end_s"])
        raise ValueError(f"Unsupported dict window format: {window_spec}")

    if isinstance(window_spec, (list, tuple)) and len(window_spec) == 2:
        start_raw, end_raw = window_spec
        start_str = str(start_raw)
        end_str = str(end_raw)

        if len(start_str) == 5:
            start_str = f"{start_str}:00"
        if len(end_str) == 5:
            end_str = f"{end_str}:00"

        start_td = pd.to_timedelta(start_str)
        end_td = pd.to_timedelta(end_str)
        return float(start_td.total_seconds()), float(end_td.total_seconds())

    raise ValueError(f"Unsupported window format: {window_spec}")


def load_pt100_data(filepath: Path) -> pd.DataFrame:
    df = pd.read_csv(filepath, skiprows=1, header=None)

    if df.shape[1] < 6:
        raise ValueError(f"Unexpected PT100 file format: {filepath}")

    if df.shape[1] >= 7:
        df = df.iloc[:, :7].copy()
        df.columns = ["time_raw"] + SENSORS + ["Average"]
        df = df.drop(columns=["Average"], errors="ignore")
    else:
        df = df.iloc[:, :6].copy()
        df.columns = ["time_raw"] + SENSORS

    df["time_s"] = to_time_seconds(df["time_raw"])
    df = df.dropna(subset=["time_s"]).sort_values("time_s").reset_index(drop=True)

    for sensor in SENSORS:
        df[sensor] = pd.to_numeric(df[sensor], errors="coerce")

    return df


def load_chiller_data(filepath: Path) -> pd.DataFrame:
    df = pd.read_csv(filepath, skiprows=1, header=None)

    if df.shape[1] < 5:
        raise ValueError(f"Unexpected chiller file format: {filepath}")

    df = df.iloc[:, :5].copy()
    df.columns = ["time_raw", "SP", "BT", "RTD0", "RTD1"]

    df["time_s"] = to_time_seconds(df["time_raw"])
    df = df.dropna(subset=["time_s"]).sort_values("time_s").reset_index(drop=True)

    for col in ["SP", "BT", "RTD0", "RTD1"]:
        df[col] = pd.to_numeric(df[col], errors="coerce")

    return df


def merge_temperature_data(pt100_df: pd.DataFrame, chiller_df: pd.DataFrame) -> pd.DataFrame:
    left = pt100_df.copy()
    right = chiller_df.copy()

    left["time_s"] = pd.to_numeric(left["time_s"], errors="coerce").astype(float)
    right["time_s"] = pd.to_numeric(right["time_s"], errors="coerce").astype(float)

    left = left.dropna(subset=["time_s"]).sort_values("time_s")
    right = right.dropna(subset=["time_s"]).sort_values("time_s")

    merged = pd.merge_asof(
        left,
        right[["time_s", "SP", "BT", "RTD0", "RTD1"]],
        on="time_s",
        direction="nearest",
        tolerance=2.0,
    )
    return merged


def extract_step_data(data: pd.DataFrame, step_windows: dict) -> pd.DataFrame:
    data = data.copy()
    data["stepSP"] = np.nan

    for sp, window_spec in step_windows.items():
        sp_int = int(sp)
        start_sec, end_sec = parse_window_to_seconds(window_spec)
        start_extract = max(end_sec - 120.0, start_sec)

        mask = (data["time_s"] >= start_extract) & (data["time_s"] <= end_sec)
        data.loc[mask, "stepSP"] = sp_int

    step_data = data.dropna(subset=["stepSP"]).copy()
    step_data["stepSP"] = step_data["stepSP"].astype(int)
    return step_data


def fit_linear_model(x: np.ndarray, y: np.ndarray) -> tuple[float, float, float, float]:
    A = np.vstack([x, np.ones_like(x)]).T
    a, b = np.linalg.lstsq(A, y, rcond=None)[0]

    y_pred = a * x + b
    resid = y - y_pred

    ss_res = np.sum(resid ** 2)
    ss_tot = np.sum((y - y.mean()) ** 2)
    r2 = 1 - (ss_res / ss_tot) if ss_tot != 0 else np.nan
    rmse = float(np.sqrt(np.mean(resid ** 2)))

    return float(a), float(b), float(r2), rmse


def find_step_pairs(data_dir: Path) -> dict:
    grouped = {}
    for path in sorted(data_dir.rglob("*.csv")):
        m = PAIR_RE.search(path.name)
        if not m:
            continue
        key = (m.group("chiller").upper(), m.group("run"))
        grouped.setdefault(key, {})
        grouped[key][m.group("dtype").lower()] = path
    return grouped


def fit_models_per_run(
    chiller: str,
    run_key: str,
    step_windows: dict,
    pair: dict,
) -> pd.DataFrame | None:
    pt100_file = pair.get("pt100")
    chiller_file = pair.get("chiller")

    if pt100_file is None or chiller_file is None:
        return None

    pt100_df = load_pt100_data(pt100_file)
    chiller_df = load_chiller_data(chiller_file)
    data = merge_temperature_data(pt100_df, chiller_df)

    step_data = extract_step_data(data, step_windows)
    if step_data.empty:
        return None

    long_data = step_data.melt(
        id_vars=["time_s", "BT", "SP", "stepSP"],
        value_vars=SENSORS,
        var_name="sensor_id",
        value_name="T_sensor",
    )

    long_data["BT"] = pd.to_numeric(long_data["BT"], errors="coerce")
    long_data["T_sensor"] = pd.to_numeric(long_data["T_sensor"], errors="coerce")
    long_data = long_data.dropna(subset=["BT", "T_sensor"])

    models = []
    for sensor_id, sub in long_data.groupby("sensor_id"):
        step_mean = (
            sub.groupby("stepSP")[["BT", "T_sensor"]]
            .mean()
            .reset_index()
            .sort_values("stepSP")
        )

        if len(step_mean) < 2:
            continue

        x = step_mean["BT"].to_numpy()
        y = step_mean["T_sensor"].to_numpy()
        a, b, r2, rmse = fit_linear_model(x, y)

        models.append(
            {
                "chiller": chiller,
                "run": run_key,
                "sensor_id": sensor_id,
                "a_slope": a,
                "b_intercept": b,
                "R2": r2,
                "RMSE": rmse,
                "n_steps": len(step_mean),
            }
        )

    if not models:
        return None

    return pd.DataFrame(models)


def compute_step_offsets(step_windows_all: dict, step_pairs: dict) -> pd.DataFrame:
    rows = []

    for chiller in CHILLERS:
        run_map = step_windows_all.get(chiller, {})
        for run_key, step_windows in sorted(run_map.items()):
            pair = step_pairs.get((chiller, run_key), {})
            pt100_file = pair.get("pt100")
            chiller_file = pair.get("chiller")

            if pt100_file is None or chiller_file is None:
                continue

            pt100_df = load_pt100_data(pt100_file)
            chiller_df = load_chiller_data(chiller_file)
            data = merge_temperature_data(pt100_df, chiller_df)
            step_data = extract_step_data(data, step_windows)

            if step_data.empty:
                continue

            for step_value, sub in step_data.groupby("stepSP"):
                bt_mean = pd.to_numeric(sub["BT"], errors="coerce").mean()

                for sensor in SENSORS:
                    if sensor not in sub.columns:
                        continue

                    sensor_values = pd.to_numeric(sub[sensor], errors="coerce").dropna()
                    if sensor_values.empty or not np.isfinite(bt_mean):
                        continue

                    offset = (sensor_values - bt_mean).mean()

                    rows.append(
                        {
                            "Chiller": chiller,
                            "Run": run_key,
                            "Sensor": sensor,
                            "Step (°C)": int(step_value),
                            "Offset (°C)": float(offset),
                        }
                    )

    return pd.DataFrame(rows)


def build_chiller_summary(combined_df: pd.DataFrame) -> pd.DataFrame:
    summary = (
        combined_df.groupby("chiller")
        .agg(
            a_mean=("a_slope", "mean"),
            a_std=("a_slope", "std"),
            b_mean=("b_intercept", "mean"),
            b_std=("b_intercept", "std"),
            RMSE_mean=("RMSE", "mean"),
            RMSE_std=("RMSE", "std"),
            R2_mean=("R2", "mean"),
            R2_std=("R2", "std"),
            n_model_rows=("sensor_id", "size"),
            n_runs=("run", pd.Series.nunique),
        )
        .reset_index()
        .rename(columns={"chiller": "Chiller"})
    )
    return summary


def build_appendix_summary(combined_df: pd.DataFrame) -> pd.DataFrame:
    appendix = (
        combined_df.groupby(["chiller", "sensor_id"])
        .agg(
            a_mean=("a_slope", "mean"),
            a_std=("a_slope", "std"),
            b_mean=("b_intercept", "mean"),
            b_std=("b_intercept", "std"),
            RMSE_mean=("RMSE", "mean"),
            RMSE_std=("RMSE", "std"),
            R2_mean=("R2", "mean"),
            n_runs=("run", pd.Series.nunique),
        )
        .reset_index()
        .rename(columns={"chiller": "Chiller", "sensor_id": "Sensor"})
    )

    appendix["Sensor_num"] = appendix["Sensor"].str.extract(r"(\d+)").astype(int)
    appendix = appendix.sort_values(["Chiller", "Sensor_num"]).drop(columns=["Sensor_num"])
    return appendix


def build_step_summaries(offsets_df: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]:
    if offsets_df.empty:
        return pd.DataFrame(), pd.DataFrame()

    step_offsets_summary = (
        offsets_df.groupby(["Chiller", "Step (°C)"])["Offset (°C)"]
        .agg(["mean", "std", "count"])
        .reset_index()
        .rename(columns={"mean": "offset_mean", "std": "offset_std", "count": "n"})
        .sort_values(["Chiller", "Step (°C)"])
    )

    sensor_step_means = (
        offsets_df.groupby(["Chiller", "Sensor", "Step (°C)"])["Offset (°C)"]
        .mean()
        .reset_index()
        .sort_values(["Chiller", "Step (°C)", "Sensor"])
    )

    step_spread_summary = (
        sensor_step_means.groupby(["Chiller", "Step (°C)"])["Offset (°C)"]
        .agg(
            min_sensor_mean="min",
            max_sensor_mean="max",
            std_sensor_mean="std",
        )
        .reset_index()
        .sort_values(["Chiller", "Step (°C)"])
    )
    step_spread_summary["spread_sensor_mean"] = (
        step_spread_summary["max_sensor_mean"] - step_spread_summary["min_sensor_mean"]
    )

    return step_offsets_summary, step_spread_summary


def plot_offset_vs_bath_temperature_condensed(
    offsets_df: pd.DataFrame,
    chiller_summary: pd.DataFrame,
) -> None:
    if offsets_df.empty:
        raise RuntimeError("No step offsets available for plotting.")

    fig, ax = plt.subplots(figsize=(5.2, 4.2))

    for chiller in CHILLERS:
        color = CHILLER_COLORS[chiller]
        sub = offsets_df[offsets_df["Chiller"] == chiller].copy()
        if sub.empty:
            continue

        grouped = (
            sub.groupby("Step (°C)")["Offset (°C)"]
            .agg(["mean", "std"])
            .reset_index()
            .sort_values("Step (°C)")
        )

        steps = grouped["Step (°C)"].to_numpy()
        means = grouped["mean"].to_numpy()
        stds = grouped["std"].to_numpy()

        ax.errorbar(
            steps,
            means,
            yerr=stds,
            fmt="o",
            markersize=4.5,
            linewidth=1.4,
            capsize=3,
            color=color,
            ecolor=color,
            elinewidth=1.0,
            label=f"{chiller} (mean ±1σ)",
        )

        row = chiller_summary[chiller_summary["Chiller"] == chiller]
        if row.empty:
            continue

        a_mean = float(row["a_mean"].iloc[0])
        b_mean = float(row["b_mean"].iloc[0])

        bt_range = np.linspace(-35, 0, 200)
        offset_pred = (a_mean - 1.0) * bt_range + b_mean

        ax.plot(
            bt_range,
            offset_pred,
            linestyle="--",
            linewidth=1.3,
            color=color,
            alpha=0.90,
            label=f"{chiller} linear model",
        )

    ax.set_xlabel(r"Bath temperature $T_{BT}$ (°C)", fontsize=10)
    ax.set_ylabel(r"Offset $T_{\mathrm{well}} - T_{BT}$ (°C)", fontsize=10)
    ax.set_xlim(-36, 1)
    ax.set_xticks([-35, -30, -25, -20, -15, -10, -5, 0])
    ax.set_ylim(0, 2.2)
    ax.grid(True, alpha=0.25, linewidth=0.6)
    ax.tick_params(axis="both", labelsize=9)
    ax.legend(fontsize=8, frameon=False, loc="upper right")

    plt.tight_layout()
    plt.savefig(OUTPUT_DIR / "fig3_offset_condensed.png", dpi=300, bbox_inches="tight")
    plt.savefig(OUTPUT_DIR / "fig3_offset_condensed.pdf", bbox_inches="tight")

    plt.close()


def write_table_transfer_latex(chiller_summary: pd.DataFrame, outpath: Path) -> None:
    lines = [
        r"\begin{tabular}{lccc}",
        r"\hline",
        r"Chiller & $a$ & $b$ ($^{\circ}$C) & RMSE ($^{\circ}$C) \\",
        r"\hline",
    ]

    for _, row in chiller_summary.iterrows():
        lines.append(
            f"{row['Chiller']} & "
            f"{row['a_mean']:.3f} $\\pm$ {row['a_std']:.3f} & "
            f"{row['b_mean']:.3f} $\\pm$ {row['b_std']:.3f} & "
            f"{row['RMSE_mean']:.3f} \\\\"
        )

    lines.extend([r"\hline", r"\end{tabular}"])
    outpath.write_text("\n".join(lines) + "\n", encoding="utf-8")


def write_appendix_latex(appendix_summary: pd.DataFrame, outpath: Path) -> None:
    lines = [
        r"\begin{tabular}{llcccccccc}",
        r"\hline",
        r" &  & \multicolumn{2}{c}{$a$} & \multicolumn{2}{c}{$b$ ($^{\circ}$C)} & \multicolumn{2}{c}{RMSE ($^{\circ}$C)} & $R^2$ & $n$ runs \\",
        r"Chiller & Sensor & mean & std & mean & std & mean & std & mean & \\",
        r"\hline",
    ]

    for _, row in appendix_summary.iterrows():
        lines.append(
            f"{row['Chiller']} & {row['Sensor']} & "
            f"{row['a_mean']:.3f} & {row['a_std']:.3f} & "
            f"{row['b_mean']:.3f} & {row['b_std']:.3f} & "
            f"{row['RMSE_mean']:.3f} & {row['RMSE_std']:.3f} & "
            f"{row['R2_mean']:.3f} & {int(row['n_runs'])} \\\\"
        )

    lines.extend([r"\hline", r"\end{tabular}"])
    outpath.write_text("\n".join(lines) + "\n", encoding="utf-8")


def build_metrics_dict(
    combined_df: pd.DataFrame,
    chiller_summary: pd.DataFrame,
    appendix_summary: pd.DataFrame,
    step_spread_summary: pd.DataFrame,
) -> dict:
    metrics = {
        "valid_model_rows": int(len(combined_df)),
        "overall_mean_rmse": float(combined_df["RMSE"].mean()),
        "overall_mean_r2": float(combined_df["R2"].mean()),
        "overall_min_r2": float(combined_df["R2"].min()),
        "chillers": {},
    }

    for chiller in CHILLERS:
        csum = chiller_summary[chiller_summary["Chiller"] == chiller]
        if csum.empty:
            continue

        app = appendix_summary[appendix_summary["Chiller"] == chiller].copy()
        spread = step_spread_summary[step_spread_summary["Chiller"] == chiller].copy()

        metrics["chillers"][chiller] = {
            "a_mean": float(csum["a_mean"].iloc[0]),
            "a_std": float(csum["a_std"].iloc[0]),
            "b_mean": float(csum["b_mean"].iloc[0]),
            "b_std": float(csum["b_std"].iloc[0]),
            "rmse_mean": float(csum["RMSE_mean"].iloc[0]),
            "rmse_std": float(csum["RMSE_std"].iloc[0]),
            "r2_mean": float(csum["R2_mean"].iloc[0]),
            "r2_std": float(csum["R2_std"].iloc[0]),
            "n_runs": int(csum["n_runs"].iloc[0]),
            "appendix_rows": [
                {
                    "sensor": row["Sensor"],
                    "a_mean": float(row["a_mean"]),
                    "a_std": float(row["a_std"]),
                    "b_mean": float(row["b_mean"]),
                    "b_std": float(row["b_std"]),
                    "rmse_mean": float(row["RMSE_mean"]),
                    "rmse_std": float(row["RMSE_std"]),
                    "r2_mean": float(row["R2_mean"]),
                    "n_runs": int(row["n_runs"]),
                }
                for _, row in app.iterrows()
            ],
        }

        for step_target in [0, -35]:
            row = spread[spread["Step (°C)"] == step_target]
            if not row.empty:
                metrics["chillers"][chiller][f"step_{step_target}_sensor_range"] = {
                    "min_sensor_mean": float(row["min_sensor_mean"].iloc[0]),
                    "max_sensor_mean": float(row["max_sensor_mean"].iloc[0]),
                    "spread_sensor_mean": float(row["spread_sensor_mean"].iloc[0]),
                    "std_sensor_mean": float(row["std_sensor_mean"].iloc[0]),
                }

    return metrics


def write_report(
    combined_df: pd.DataFrame,
    chiller_summary: pd.DataFrame,
    appendix_summary: pd.DataFrame,
    step_spread_summary: pd.DataFrame,
    outpath: Path,
) -> None:
    lines = []
    lines.append("FrESH step-analysis report")
    lines.append("=" * 80)
    lines.append("")
    lines.append(f"Valid model rows: {len(combined_df)}")
    lines.append(f"Overall mean RMSE: {combined_df['RMSE'].mean():.3f} °C")
    lines.append(f"Overall mean R²: {combined_df['R2'].mean():.4f}")
    lines.append(f"Overall minimum R²: {combined_df['R2'].min():.4f}")
    lines.append("")

    lines.append("Table 1 values")
    lines.append("-" * 80)
    for _, row in chiller_summary.iterrows():
        lines.append(
            f"{row['Chiller']}: "
            f"a = {row['a_mean']:.3f} ± {row['a_std']:.3f}, "
            f"b = {row['b_mean']:.3f} ± {row['b_std']:.3f} °C, "
            f"RMSE = {row['RMSE_mean']:.3f} °C, "
            f"R² = {row['R2_mean']:.3f}, "
            f"n_runs = {int(row['n_runs'])}"
        )
    lines.append("")

    lines.append("Appendix A1 values")
    lines.append("-" * 80)
    for _, row in appendix_summary.iterrows():
        lines.append(
            f"{row['Chiller']} {row['Sensor']}: "
            f"a = {row['a_mean']:.3f} ± {row['a_std']:.3f}, "
            f"b = {row['b_mean']:.3f} ± {row['b_std']:.3f} °C, "
            f"RMSE = {row['RMSE_mean']:.3f} ± {row['RMSE_std']:.3f} °C, "
            f"R² = {row['R2_mean']:.3f}, "
            f"n_runs = {int(row['n_runs'])}"
        )
    lines.append("")

    lines.append("Step-offset ranges based on mean offset per sensor at each step")
    lines.append("-" * 80)
    for chiller in CHILLERS:
        lines.append(chiller)
        sub = step_spread_summary[step_spread_summary["Chiller"] == chiller].sort_values("Step (°C)")
        for _, row in sub.iterrows():
            lines.append(
                f"  Step {int(row['Step (°C)']):>3} °C: "
                f"min = {row['min_sensor_mean']:.3f} °C, "
                f"max = {row['max_sensor_mean']:.3f} °C, "
                f"spread = {row['spread_sensor_mean']:.3f} °C, "
                f"σ_sensors = {row['std_sensor_mean']:.3f} °C"
            )
        lines.append("")

    outpath.write_text("\n".join(lines) + "\n", encoding="utf-8")


def print_paper_summary(
    combined_df: pd.DataFrame,
    chiller_summary: pd.DataFrame,
    step_spread_summary: pd.DataFrame,
) -> None:
    print("\n" + "=" * 80)
    print("STEP ANALYSIS COMPLETE")
    print("=" * 80)

    print("\nPaper-facing outputs written to output/:")
    print("  - fig3_offset_condensed.png")
    print("  - table_transfer.csv")
    print("  - appendix_transfer.csv")
    print("  - step_report.txt")
    print("  - step_text_metrics.json")

    if SAVE_LATEX_TABLES:
        print("  - table_transfer.tex")
        print("  - appendix_transfer.tex")

    print("\nTable 1 values:")
    for _, row in chiller_summary.iterrows():
        print(
            f"  {row['Chiller']}: "
            f"a = {row['a_mean']:.3f} ± {row['a_std']:.3f}, "
            f"b = {row['b_mean']:.3f} ± {row['b_std']:.3f} °C, "
            f"RMSE = {row['RMSE_mean']:.3f} °C, "
            f"R² = {row['R2_mean']:.3f}"
        )

    print("\nResults-paragraph checks:")
    print(
        f"  Overall mean RMSE across all sensor-run fits: {combined_df['RMSE'].mean():.3f} °C"
    )
    print(
        f"  Overall mean R² across all sensor-run fits: {combined_df['R2'].mean():.4f}"
    )
    print(
        f"  Minimum R² across all sensor-run fits: {combined_df['R2'].min():.4f}"
    )

    for chiller in CHILLERS:
        sub = step_spread_summary[step_spread_summary["Chiller"] == chiller]

        row0 = sub[sub["Step (°C)"] == 0]
        row35 = sub[sub["Step (°C)"] == -35]

        if not row0.empty:
            print(
                f"  {chiller} at 0 °C: "
                f"{row0['min_sensor_mean'].iloc[0]:.2f} to {row0['max_sensor_mean'].iloc[0]:.2f} °C "
                f"(spread {row0['spread_sensor_mean'].iloc[0]:.2f} °C)"
            )

        if not row35.empty:
            print(
                f"  {chiller} at -35 °C: "
                f"{row35['min_sensor_mean'].iloc[0]:.2f} to {row35['max_sensor_mean'].iloc[0]:.2f} °C "
                f"(spread {row35['spread_sensor_mean'].iloc[0]:.2f} °C)"
            )


def main() -> None:
    ensure_output_dir()

    step_windows_all = load_step_windows()
    step_pairs = find_step_pairs(DATA_DIR)

    all_models = []
    for chiller in CHILLERS:
        run_map = step_windows_all.get(chiller, {})
        for run_key in sorted(run_map.keys()):
            pair = step_pairs.get((chiller, run_key), {})
            result = fit_models_per_run(chiller, run_key, run_map[run_key], pair)
            if result is not None:
                all_models.append(result)

    if not all_models:
        raise RuntimeError("No valid step runs were processed.")

    combined_df = pd.concat(all_models, ignore_index=True)
    chiller_summary = build_chiller_summary(combined_df)
    appendix_summary = build_appendix_summary(combined_df)

    offsets_df = compute_step_offsets(step_windows_all, step_pairs)
    step_offsets_summary, step_spread_summary = build_step_summaries(offsets_df)

    plot_offset_vs_bath_temperature_condensed(offsets_df, chiller_summary)

    chiller_summary.to_csv(OUTPUT_DIR / "table_transfer.csv", index=False)
    appendix_summary.to_csv(OUTPUT_DIR / "appendix_transfer.csv", index=False)

    if not step_offsets_summary.empty:
        step_offsets_summary.to_csv(OUTPUT_DIR / "step_offsets_summary.csv", index=False)
    if not step_spread_summary.empty:
        step_spread_summary.to_csv(OUTPUT_DIR / "step_sensor_offset_summary.csv", index=False)

    if SAVE_SUPPORTING_CSV:
        combined_df.to_csv(OUTPUT_DIR / "all_characterization_models.csv", index=False)
        if not offsets_df.empty:
            offsets_df.to_csv(OUTPUT_DIR / "step_offsets.csv", index=False)

    if SAVE_LATEX_TABLES:
        write_table_transfer_latex(chiller_summary, OUTPUT_DIR / "table_transfer.tex")
        write_appendix_latex(appendix_summary, OUTPUT_DIR / "appendix_transfer.tex")

    metrics = build_metrics_dict(
        combined_df=combined_df,
        chiller_summary=chiller_summary,
        appendix_summary=appendix_summary,
        step_spread_summary=step_spread_summary,
    )

    with open(OUTPUT_DIR / "step_text_metrics.json", "w", encoding="utf-8") as f:
        json.dump(metrics, f, indent=2)

    write_report(
        combined_df=combined_df,
        chiller_summary=chiller_summary,
        appendix_summary=appendix_summary,
        step_spread_summary=step_spread_summary,
        outpath=OUTPUT_DIR / "step_report.txt",
    )

    print_paper_summary(
        combined_df=combined_df,
        chiller_summary=chiller_summary,
        step_spread_summary=step_spread_summary,
    )


if __name__ == "__main__":
    main()