shared_utils.py

"""
shared_utils.py
===============
Single shared module for the NIH Catalyze Smart Catheter project.

Covers three analysis pipelines:
  - Ammonia/CO₂ extraction: Master sheet loaders, ratio computation
  - Spline interpolation: NaOH @ pH=11 (PCHIP) and citric acid @ pH=4 (linear)
  - Biomarker analysis: statistics, by-pig plots, summary table

Pipeline order (CSV/PNG artifacts default to ``catalyze/output/``):
  ammonia_analysis.ipynb       → output/ammonia_*.csv, …
  spline_interpolation.ipynb   → output/naoh_*.csv, output/citric_*.csv, …
  biomarker_analysis.ipynb     → output/biomarker_summary.csv + PNGs

Dependencies: pandas, numpy, scipy, matplotlib, openpyxl
"""

import re
import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from scipy.interpolate import PchipInterpolator, interp1d
from scipy.stats import mannwhitneyu, ttest_ind

# Avoid globally suppressing warnings in a shared library.
# If a specific warning stream becomes noisy, handle it at the call site.


# ============================================================================
# CONSTANTS
# ============================================================================

HEALTHY_PIGS    = ["PIG-04", "PIG-06", "PIG-07", "PIG-08"]
NONHEALTHY_PIGS = ["PIG-01", "PIG-02", "PIG-03", "PIG-05", "PIG-09", "PIG-10"]
ALL_PIGS_ORDER  = [f"PIG-{i:02d}" for i in range(1, 11)]

# Master-sheet column names (p1 … p10) grouped for time-series plots
GROUPS_PIG_COLS = {
    "Healthy Controls":     ["p4", "p6", "p7", "p8"],
    "Non-Healthy Controls": ["p1", "p2", "p3", "p5", "p9", "p10"],
}

HEALTHY_COLOR        = "#90EE90"   # light green
NONHEALTHY_COLOR     = "#FFB6C1"   # light pink
MEAN_LINE_HEALTHY    = "#228B22"   # dark green dashed mean line
MEAN_LINE_NONHEALTHY = "#C71585"   # dark pink dashed mean line

# Repository root (folder containing this file) and default artifact directory
CATALYZE_DIR         = Path(__file__).resolve().parent
DEFAULT_OUTPUT_DIR   = CATALYZE_DIR / "output"


def ensure_parent_dir(file_path):
    """
    Create the parent directory for a file path if it is missing.

    Parameters
    ----------
    file_path : str or pathlib.Path
        Path whose parent directory should exist before writing a file.
    """
    Path(file_path).parent.mkdir(parents=True, exist_ok=True)


# ============================================================================
# STATISTICS
# ============================================================================

def calculate_statistics(healthy_data, nonhealthy_data):
    """
    Compare two groups with t-test, Mann-Whitney U, and Cohen's d.

    Parameters
    ----------
    healthy_data : array-like
        Observations for the healthy group.
    nonhealthy_data : array-like
        Observations for the non-healthy group.

    Returns
    -------
    dict
        Keys: healthy_mean, healthy_std, nonhealthy_mean, nonhealthy_std,
        p_value_ttest, p_value_mannwhitney, cohens_d, diff_percent,
        significance_status, effect_size.
    """
    healthy_mean    = np.mean(healthy_data)
    healthy_std     = np.std(healthy_data, ddof=1)
    nonhealthy_mean = np.mean(nonhealthy_data)
    nonhealthy_std  = np.std(nonhealthy_data, ddof=1)

    t_stat, p_val_ttest = ttest_ind(healthy_data, nonhealthy_data)
    u_stat, p_val_mw    = mannwhitneyu(healthy_data, nonhealthy_data, alternative="two-sided")

    pooled_std = np.sqrt((healthy_std**2 + nonhealthy_std**2) / 2)
    cohens_d   = (healthy_mean - nonhealthy_mean) / pooled_std
    diff_pct   = 100 * (healthy_mean - nonhealthy_mean) / nonhealthy_mean

    if p_val_ttest < 0.001:
        sig_status = "Highly Significant"
    elif p_val_ttest < 0.01:
        sig_status = "Very Significant"
    elif p_val_ttest < 0.05:
        sig_status = "Significant"
    else:
        sig_status = "Not Significant"

    d_abs = abs(cohens_d)
    if d_abs < 0.2:
        effect_size = "None"
    elif d_abs < 0.5:
        effect_size = "Small"
    elif d_abs < 0.8:
        effect_size = "Medium"
    else:
        effect_size = "Large"

    return {
        "healthy_mean":          healthy_mean,
        "healthy_std":           healthy_std,
        "nonhealthy_mean":       nonhealthy_mean,
        "nonhealthy_std":        nonhealthy_std,
        "p_value_ttest":         p_val_ttest,
        "p_value_mannwhitney":   p_val_mw,
        "cohens_d":              cohens_d,
        "diff_percent":          diff_pct,
        "significance_status":   sig_status,
        "effect_size":           effect_size,
    }


# ============================================================================
# PLOTTING — by-pig boxplots
# ============================================================================

def pig_facecolor(pig_id):
    """Return box face color for a pig (green = healthy, pink = non-healthy)."""
    return HEALTHY_COLOR if pig_id in HEALTHY_PIGS else NONHEALTHY_COLOR


def format_mean_for_legend(x):
    """Format a pooled mean for legend text, adapting to its magnitude."""
    if not np.isfinite(x):
        return "n/a"
    axv = abs(x)
    if axv >= 1000 or (axv > 0 and axv < 1e-3):
        return f"{x:.3e}"
    if axv >= 100:
        return f"{x:.1f}"
    if axv >= 1:
        return f"{x:.2f}"
    return f"{x:.3g}"


def plot_biomarkers_by_pig_figure(panels, suptitle, subtitle, output_path, y_axis_scales=None):
    """
    One row of subplots: per-pig boxplot (PIG-01 … PIG-10), green vs pink by group.

    Parameters
    ----------
    panels : list of tuple
        Each entry: (panel_letter, distribution_title, y_axis_label, pig_to_values)
        where pig_to_values maps 'PIG-XX' to a 1-D array of measurements.
    suptitle : str
        Figure title.
    subtitle : str
        Subtitle line (e.g. legend description).
    output_path : str or pathlib.Path
        Output PNG path (300 dpi).
    y_axis_scales : list of str, optional
        Same length as panels; each entry 'linear' or 'log'.
    """
    n = len(panels)
    if n == 0:
        return

    scales = (["linear"] * n) if y_axis_scales is None else list(y_axis_scales)
    if len(scales) != n:
        raise ValueError(f"y_axis_scales length ({len(scales)}) must match panels ({n}).")

    fig, axes = plt.subplots(1, n, figsize=(6.2 * n, 6.0))
    if n == 1:
        axes = [axes]

    for ax, (letter, dist_title, ylabel, pig_value_dict), yscale in zip(axes, panels, scales):
        data       = []
        facecolors = []
        for pig in ALL_PIGS_ORDER:
            vals = np.asarray(pig_value_dict.get(pig, []), dtype=float)
            vals = vals[np.isfinite(vals)]
            if yscale == "log":
                vals = vals[vals > 0]
            data.append(vals.tolist() if vals.size else [])
            facecolors.append(pig_facecolor(pig))

        bp = ax.boxplot(
            data,
            labels=ALL_PIGS_ORDER,
            patch_artist=True,
            widths=0.55,
            showfliers=True,
            medianprops=dict(color="black", linewidth=2),
            boxprops=dict(linewidth=1.5),
        )
        for patch, fc in zip(bp["boxes"], facecolors):
            patch.set_facecolor(fc)

        if yscale == "log":
            ax.set_yscale("log")

        h_parts  = [np.asarray(data[i]) for i, pig in enumerate(ALL_PIGS_ORDER)
                    if pig in HEALTHY_PIGS    and len(data[i]) > 0]
        nh_parts = [np.asarray(data[i]) for i, pig in enumerate(ALL_PIGS_ORDER)
                    if pig in NONHEALTHY_PIGS and len(data[i]) > 0]
        h_pooled  = np.concatenate(h_parts)  if h_parts  else np.array([])
        nh_pooled = np.concatenate(nh_parts) if nh_parts else np.array([])

        hm  = float(np.mean(h_pooled))  if h_pooled.size  else float("nan")
        nhm = float(np.mean(nh_pooled)) if nh_pooled.size else float("nan")
        if h_pooled.size:
            ax.axhline(hm,  color=MEAN_LINE_HEALTHY,    linestyle="--", linewidth=2.2, zorder=1)
        if nh_pooled.size:
            ax.axhline(nhm, color=MEAN_LINE_NONHEALTHY, linestyle="--", linewidth=2.2, zorder=1)

        leg_lines = [
            Line2D([0], [0], color=MEAN_LINE_HEALTHY,    linestyle="--", linewidth=2.2),
            Line2D([0], [0], color=MEAN_LINE_NONHEALTHY, linestyle="--", linewidth=2.2),
        ]
        ax.legend(
            leg_lines,
            [f"Healthy mean\n{format_mean_for_legend(hm)}",
             f"Non-Healthy mean\n{format_mean_for_legend(nhm)}"],
            loc="upper left", fontsize=9, framealpha=0.95,
        )

        ax.set_title(f"{letter}. {dist_title}", fontsize=12, fontweight="bold", pad=8)
        ax.set_xlabel("Pig ID",  fontsize=11, fontweight="bold")
        ax.set_ylabel(ylabel,    fontsize=11, fontweight="bold")
        ax.tick_params(axis="x", rotation=45)
        ax.grid(True, alpha=0.3, axis="y", which="both" if yscale == "log" else "major")

    fig.suptitle(suptitle, fontsize=14, fontweight="bold", y=1.02)
    fig.text(0.5, 0.965, subtitle, ha="center", fontsize=11, style="italic")
    plt.tight_layout(rect=[0, 0, 1, 0.88])
    out = Path(output_path)
    out.parent.mkdir(parents=True, exist_ok=True)
    plt.savefig(out, dpi=300, bbox_inches="tight")
    plt.close()
    print(f"✓ Saved: {out}")


def plot_timeseries_by_group(df, ylabel, title, log_y=False, groups=None):
    """
    Plot time-series concentration data side-by-side for each pig group.

    Parameters
    ----------
    df : pandas.DataFrame
        Wide format with a 'Cumulative time (min)' column and pig columns (p1, …).
    ylabel : str
        Y-axis label.
    title : str
        Figure super-title.
    log_y : bool
        Use log scale on the y-axis.
    groups : dict, optional
        Maps group label → list of pig column names. Defaults to GROUPS_PIG_COLS.
    """
    if groups is None:
        groups = GROUPS_PIG_COLS

    fig, axes = plt.subplots(1, len(groups), figsize=(7 * len(groups), 6))
    if len(groups) == 1:
        axes = [axes]

    for ax, (group_label, pig_cols) in zip(axes, groups.items()):
        for pig_col in pig_cols:
            if pig_col not in df.columns:
                continue
            mask = df[pig_col].notna()
            ax.plot(
                df.loc[mask, "Cumulative time (min)"],
                df.loc[mask, pig_col],
                marker="o", markersize=4, linewidth=1.5, label=pig_col,
            )
        if log_y:
            ax.set_yscale("log")
        ax.set_title(group_label)
        ax.set_xlabel("Cumulative time (min)")
        ax.set_ylabel(ylabel)
        ax.legend(title="Pig", bbox_to_anchor=(1.01, 1), loc="upper left")
        ax.grid(True, alpha=0.3)

    all_ylims = [ax.get_ylim() for ax in axes]
    for ax in axes:
        ax.set_ylim(min(lim[0] for lim in all_ylims), max(lim[1] for lim in all_ylims))

    plt.suptitle(title)
    plt.tight_layout()
    plt.show()


def plot_combined_timeseries(ammonia_df, co2_df, ratio_df):
    """
    Three-panel figure: ammonia, CO₂, and ratio vs cumulative time (all pigs).

    Parameters
    ----------
    ammonia_df : pandas.DataFrame
    co2_df     : pandas.DataFrame
    ratio_df   : pandas.DataFrame
        All three must have a 'Cumulative time (min)' column and pig columns.
    """
    fig, (ax_amm, ax_co2, ax_ratio) = plt.subplots(1, 3, figsize=(21, 6))

    pig_cols = [c for c in ammonia_df.columns if c.startswith("p")]
    for pig_col in pig_cols:
        for ax, df, log_y in [
            (ax_amm,   ammonia_df, False),
            (ax_co2,   co2_df,    False),
            (ax_ratio, ratio_df,  True),
        ]:
            if pig_col not in df.columns:
                continue
            mask = df[pig_col].notna()
            ax.plot(
                df.loc[mask, "Cumulative time (min)"],
                df.loc[mask, pig_col],
                marker="o", markersize=4, linewidth=1.5, label=pig_col,
            )

    panel_cfg = [
        (ax_amm,   "Ammonia concentration (mM)",      "Ammonia",             False),
        (ax_co2,   "Carbon dioxide concentration (mM)", "Carbon Dioxide",     False),
        (ax_ratio, "Ammonia / CO\u2082 ratio (log scale)", "Ammonia / CO\u2082", True),
    ]
    for ax, ylabel, panel_title, log_y in panel_cfg:
        ax.set_title(panel_title)
        ax.set_xlabel("Cumulative time (min)")
        ax.set_ylabel(ylabel)
        if log_y:
            ax.set_yscale("log")
        ax.legend(title="Pig", bbox_to_anchor=(1.01, 1), loc="upper left")
        ax.grid(True, alpha=0.3)

    plt.suptitle("Ammonia, CO\u2082, and Ammonia/CO\u2082 ratio vs cumulative time (Control Pigs)")
    plt.tight_layout()
    plt.show()


# ============================================================================
# DATA HELPERS
# ============================================================================

def pooled_from_pig_series_dict(pig_value_dict, pig_ids):
    """
    Concatenate all finite observations for pigs listed in pig_ids.

    Parameters
    ----------
    pig_value_dict : dict
        Maps 'PIG-XX' to a 1-D array-like of values.
    pig_ids : iterable of str
        Pig identifiers to include.

    Returns
    -------
    numpy.ndarray
        Pooled finite values, shape (0,) if none available.
    """
    parts = []
    for p in pig_ids:
        v = np.asarray(pig_value_dict.get(p, []), dtype=float)
        v = v[np.isfinite(v)]
        if v.size:
            parts.append(v)
    return np.concatenate(parts) if parts else np.array([])


def master_wide_df_to_pig_dict(df):
    """
    Convert a wide Master-sheet DataFrame (columns p1 … p10) to a pig dict.

    Parameters
    ----------
    df : pandas.DataFrame
        Must include p1 … p10 columns (or a subset).

    Returns
    -------
    dict
        'PIG-XX' → 1-D numpy.ndarray of finite values.
    """
    out = {}
    for pig in ALL_PIGS_ORDER:
        col = f"p{int(pig.split('-')[1])}"
        if col not in df.columns:
            out[pig] = np.array([])
            continue
        v = df[col].dropna().to_numpy(dtype=float)
        out[pig] = v[np.isfinite(v)]
    return out


def pig_ids_to_master_sheet_columns(pig_ids):
    """
    Map 'PIG-04' style identifiers to Master-sheet column names ('p4').

    Parameters
    ----------
    pig_ids : list of str

    Returns
    -------
    list of str
    """
    return [f"p{int(p.split('-')[1])}" for p in pig_ids]


def pooled_biomarker_values_for_pigs(df, pig_ids):
    """
    Pool all non-null, finite values for the given pigs from a wide DataFrame.

    Parameters
    ----------
    df : pandas.DataFrame
        Wide table with one column per pig (p1, …).
    pig_ids : list of str
        'PIG-' identifiers.

    Returns
    -------
    numpy.ndarray
        1-D pooled array (may be empty).
    """
    cols  = pig_ids_to_master_sheet_columns(pig_ids)
    parts = []
    for col in cols:
        if col not in df.columns:
            continue
        series = df[col].dropna().to_numpy(dtype=float)
        series = series[np.isfinite(series)]
        if series.size:
            parts.append(series)
    return np.concatenate(parts) if parts else np.array([])


def pig_dict_from_spline_long_csv(csv_path, value_column_prefix):
    """
    Build a per-pig value dict from a long-format spline results CSV.

    Parameters
    ----------
    csv_path : str or pathlib.Path
        CSV with 'Pig' column and a value column.
    value_column_prefix : str
        The first column whose stripped name starts with this prefix is used
        (e.g. 'NaOH @ pH=11', 'Citric @ pH=4').

    Returns
    -------
    dict
        'PIG-01' … 'PIG-10' → 1-D float arrays (empty if data missing).
    """
    out  = {p: np.array([], dtype=float) for p in ALL_PIGS_ORDER}
    path = Path(csv_path)
    if not path.is_file():
        return out
    df = pd.read_csv(path)
    if "Pig" not in df.columns:
        return out
    val_col = next(
        (c for c in df.columns if str(c).strip().startswith(value_column_prefix)),
        None,
    )
    if val_col is None:
        return out
    for pig, grp in df.groupby("Pig"):
        pig = str(pig).strip()
        if pig not in out:
            continue
        v = grp[val_col].dropna().to_numpy(dtype=float)
        out[pig] = v[np.isfinite(v)]
    return out


def scale_pig_series_dict(pig_dict, factor):
    """
    Multiply each pig's observation vector by a scalar (e.g. mol/L → mM).

    Parameters
    ----------
    pig_dict : dict
        Maps pig ID to 1-D array-like.
    factor : float

    Returns
    -------
    dict
        New dict with scaled numpy.ndarray values.
    """
    return {pig: np.asarray(vals, dtype=float) * factor for pig, vals in pig_dict.items()}


# ============================================================================
# EXCEL LOADERS — phosphate, creatinine, urine pH
# ============================================================================

def read_phosphate_creatinine_by_pig(input_file):
    """
    Load per-pig phosphate and creatinine series from the titration workbook.

    Parameters
    ----------
    input_file : str or pathlib.Path
        Excel workbook with a 'PHOSPHATE and CREATININE' sheet.

    Returns
    -------
    tuple of dict
        (phosphate_data, creatinine_data), each mapping 'PIG-XX' to values.
    """
    df = pd.read_excel(input_file, sheet_name="PHOSPHATE and CREATININE", nrows=13)

    phosphate_cols = [3, 6, 9, 12, 15, 18, 21, 24, 28, 31]
    pig_names      = [f"PIG-{i:02d}" for i in range(1, 11)]
    phosphate_data = {}
    for pig, col_idx in zip(pig_names, phosphate_cols):
        raw = df.iloc[:, col_idx].dropna().to_numpy(dtype=float)
        phosphate_data[pig] = raw[np.isfinite(raw)]

    creatinine_cols = [38, 41, 44, 47, 50, 53, 56, 59, 62, 65]
    creatinine_pig_order = [
        "PIG-03", "PIG-01", "PIG-02", "PIG-04", "PIG-05",
        "PIG-06", "PIG-07", "PIG-08", "PIG-09", "PIG-10",
    ]
    creatinine_data = {}
    for pig, col_idx in zip(creatinine_pig_order, creatinine_cols):
        raw = df.iloc[:, col_idx].dropna().to_numpy(dtype=float)
        creatinine_data[pig] = raw[np.isfinite(raw)]

    return phosphate_data, creatinine_data


def load_urine_ph_by_pig(input_file):
    """
    Load direct urine pH per pig from the 'pHi, OSM, Vol rate' sheet.

    Parameters
    ----------
    input_file : str or pathlib.Path

    Returns
    -------
    dict
        'PIG-XX' → 1-D array of pH values.
    """
    df  = pd.read_excel(input_file, sheet_name="pHi, OSM, Vol rate", nrows=18)
    out = {}
    for pig in ALL_PIGS_ORDER:
        if pig in df.columns:
            v = df[pig].dropna().to_numpy(dtype=float)
            out[pig] = v[np.isfinite(v)]
        else:
            out[pig] = np.array([])
    return out


# ============================================================================
# AMMONIA / CO₂ — Master sheet extraction
# ============================================================================

_SAMPLE_LABEL_RE = re.compile(r"^UB?\d{2,}$", re.IGNORECASE)


def _is_sample_label(value):
    """Return True only if value matches expected sample ID pattern (U01, UB03, …)."""
    return bool(_SAMPLE_LABEL_RE.match(str(value).strip()))


def _label_sort_key(label):
    """Sort UB labels before U labels, both in ascending numeric order."""
    m = re.match(r"^(UB?)(\d+)$", label, re.IGNORECASE)
    if m:
        return (0 if m.group(1).upper() == "UB" else 1, int(m.group(2)))
    return (2, label)


def _cumulative_time(label, max_ub):
    """
    Convert a sample label to cumulative time in minutes.

    The last UB label maps to time 0; earlier UBs step back by -20 each.
    U labels step forward by 20 (U01 = 20 min, U02 = 40 min, …).
    """
    m   = re.match(r"^(UB?)(\d+)$", label, re.IGNORECASE)
    num = int(m.group(2))
    if m.group(1).upper() == "UB":
        return (num - max_ub) * 20
    return num * 20


def _load_sheet(filepath, sheet_name, label_col):
    """
    Load an Excel sheet, auto-detecting the header row by locating label_col.

    Parameters
    ----------
    filepath : str or pathlib.Path
    sheet_name : str
    label_col : str
        Column name used to identify the header row.

    Returns
    -------
    pandas.DataFrame
        Data below the header row, with proper column names.
    """
    df_raw     = pd.read_excel(filepath, sheet_name=sheet_name, header=None)
    header_row = 0
    label_lower = label_col.strip().lower()
    for i in range(min(10, len(df_raw))):
        if (df_raw.iloc[i].astype(str).str.strip().str.lower() == label_lower).any():
            header_row = i
            break
    df_raw.columns = df_raw.iloc[header_row]
    return df_raw.iloc[header_row + 1:].reset_index(drop=True)


def _find_column(df, patterns, match_type="exact"):
    """
    Find a DataFrame column by name matching.

    Parameters
    ----------
    df : pandas.DataFrame
    patterns : list of str
        For 'exact': normalized column name must equal one of the patterns.
        For 'contains': column name must contain all substrings.
    match_type : {'exact', 'contains'}

    Returns
    -------
    Column label or None.
    """
    for col in df.columns:
        c_lower = re.sub(r"\s+", " ", str(col).strip()).lower()
        if match_type == "exact":
            normalized = [re.sub(r"\s+", " ", p).lower() for p in patterns]
            if c_lower in normalized:
                return col
        else:
            if all(s.lower() in c_lower for s in patterns):
                return col
    return None


def _extract_pig_series(df, pig_col_name, label_col, conc_col):
    """
    Extract label and concentration from a loaded sheet for one pig.

    Parameters
    ----------
    df : pandas.DataFrame
        Output of _load_sheet.
    pig_col_name : str
        Column name to assign to concentration values (e.g. 'p4').
    label_col : str
        Name of the sample label column.
    conc_col : str
        Name of the concentration column (exact match tried first, then contains).

    Returns
    -------
    pandas.DataFrame or None
        Two columns ['Label', pig_col_name], filtered to valid sample IDs.
    """
    lc = _find_column(df, [label_col], match_type="exact")
    cc = _find_column(df, [conc_col],  match_type="exact")
    if cc is None:
        keywords = [w for w in re.split(r"[\s()]+", conc_col) if len(w) > 1]
        cc = _find_column(df, keywords, match_type="contains")
    if lc is None or cc is None:
        return None
    result          = df[[lc, cc]].copy()
    result.columns  = ["Label", pig_col_name]
    result["Label"] = result["Label"].astype(str).str.strip()
    return result[result["Label"].apply(_is_sample_label)].copy()


def _get_pig_number(filename):
    """Extract pig number from a Master sheet filename (PFI01 → 1, PFI10 → 10)."""
    m = re.search(r"PFI(\d+)", filename, re.IGNORECASE)
    return int(m.group(1)) if m else None


def _load_analyte_from_master_sheets(filepaths, sheet_name, label_col, conc_col, verbose=True):
    """
    Load one analyte from all Master sheet workbooks.

    Returns
    -------
    tuple
        (data_dict, pig_cols, all_labels) where data_dict maps (label, pig_col) → value.
    """
    data     = {}
    pig_cols = []
    for filepath in filepaths:
        pig_num = _get_pig_number(filepath.name)
        if verbose:
            print(f"  {filepath.name} ...", end=" ", flush=True)
        if pig_num is None:
            if verbose:
                print("skipped (no pig number)")
            continue
        pig_col = f"p{pig_num}"
        df_raw  = _load_sheet(filepath, sheet_name, label_col)
        df_pig  = _extract_pig_series(df_raw, pig_col, label_col, conc_col)
        if df_pig is None:
            if verbose:
                print("skipped (missing columns)")
            continue
        pig_cols.append(pig_col)
        for _, row in df_pig.iterrows():
            key = (row["Label"], pig_col)
            if key not in data:
                data[key] = row[pig_col]
        if verbose:
            print("OK")

    all_labels = sorted({k[0] for k in data.keys()}, key=_label_sort_key)
    return data, pig_cols, all_labels


def _build_analyte_df(data, pig_cols, all_labels):
    """
    Build a wide time-series DataFrame from an extracted analyte data dict.

    Parameters
    ----------
    data : dict
        Maps (label, pig_col) → concentration value.
    pig_cols : list of str
    all_labels : list of str

    Returns
    -------
    pandas.DataFrame
        Columns: Label, Cumulative time (min), p1, p2, …
    """
    df = pd.DataFrame({"Label": all_labels})
    for pig_col in pig_cols:
        df[pig_col] = df["Label"].map(lambda L, pc=pig_col: data.get((L, pc), np.nan))

    ub_nums = [
        int(re.match(r"^UB(\d+)$", L, re.IGNORECASE).group(1))
        for L in df["Label"] if re.match(r"^UB\d+$", L, re.IGNORECASE)
    ]
    max_ub = max(ub_nums) if ub_nums else 0
    df.insert(
        1, "Cumulative time (min)",
        df["Label"].apply(lambda L: _cumulative_time(L, max_ub)),
    )
    return df


def load_ammonia_co2_dataframes(master_dir="Master sheets", verbose=True):
    """
    Load all Master sheet Excel files and return ammonia, CO₂, and ratio DataFrames.

    Parameters
    ----------
    master_dir : str or pathlib.Path
        Folder containing PFI##-Control-Master sheet_*.xlsx workbooks.
    verbose : bool
        Print per-file progress.

    Returns
    -------
    tuple of pandas.DataFrame
        (ammonia_df, co2_df, ratio_df) in wide format:
        columns = [Label, Cumulative time (min), p1, p2, …, p10].
    """
    master_path = Path(master_dir)
    if not master_path.exists():
        raise FileNotFoundError(f"Master sheets folder not found: {master_path}")

    filepaths = [
        f for f in sorted(master_path.glob("*.xlsx")) + sorted(master_path.glob("*.xls"))
        if not f.name.startswith("~")
    ]

    if verbose:
        print("Loading Ammonia sheets …")
    amm_data, amm_cols, amm_labels = _load_analyte_from_master_sheets(
        filepaths, "Ammonia", "Label", "Ammonia concentration (mM)", verbose=verbose,
    )
    ammonia_df = _build_analyte_df(amm_data, amm_cols, amm_labels)

    if verbose:
        print("\nLoading Carbon dioxide sheets …")
    co2_data, co2_cols, co2_labels = _load_analyte_from_master_sheets(
        filepaths, "Carbon dioxide", "Sample", "Concentration (mM)", verbose=verbose,
    )
    co2_df = _build_analyte_df(co2_data, co2_cols, co2_labels)

    ratio_df  = ammonia_df[["Label", "Cumulative time (min)"]].copy()
    pig_cols  = [c for c in ammonia_df.columns if c.startswith("p")]
    for pig_col in pig_cols:
        if pig_col in co2_df.columns:
            ratio_df[pig_col] = ammonia_df[pig_col].values / co2_df[pig_col].values

    if verbose:
        print(f"\nAmmonia shape: {ammonia_df.shape}, CO₂ shape: {co2_df.shape}, "
              f"Ratio shape: {ratio_df.shape}")

    return ammonia_df, co2_df, ratio_df


# ============================================================================
# SPLINE INTERPOLATION
# ============================================================================

def extract_titration_data(filepath, sheet_name, verbose=True):
    """
    Extract NaOH titration data from an Excel ACIDITY sheet (columns A–K).

    Column B = NaOH concentration; columns C–K = pH per time point.

    Parameters
    ----------
    filepath : str or pathlib.Path
    sheet_name : str
        Sheet name (typically 'ACIDITY').
    verbose : bool

    Returns
    -------
    dict
        'PIG-XX' → {'time_points', 'naoh_conc', 'ph_matrix'}.
    """
    if verbose:
        print("=" * 80)
        print("EXTRACTING DATA FROM EXCEL (ACIDITY)")
        print("=" * 80)

    df       = pd.read_excel(filepath, sheet_name=sheet_name)
    df_cols  = df.iloc[:, :11]

    pig_starts = [
        (str(val), idx)
        for idx, val in enumerate(df_cols.iloc[:, 0])
        if pd.notna(val) and "PIG-" in str(val)
    ]
    if verbose:
        print(f"\nFound {len(pig_starts)} pigs")

    all_pig_data = {}
    for pig_name, start_row in pig_starts:
        time_points = [
            float(val)
            for val in df_cols.iloc[start_row, 2:11]
            if pd.notna(val) and str(val) != "min" and _can_float(val)
        ]
        naoh_conc = []
        ph_matrix = []
        row_idx   = start_row + 2

        while row_idx < len(df_cols):
            naoh_val = df_cols.iloc[row_idx, 1]
            if pd.isna(naoh_val):
                break
            try:
                naoh_conc.append(float(naoh_val))
                ph_row = [
                    float(df_cols.iloc[row_idx, c]) if pd.notna(df_cols.iloc[row_idx, c]) else np.nan
                    for c in range(2, 11)
                ]
                ph_matrix.append(ph_row)
                row_idx += 1
            except (ValueError, TypeError):
                break

        all_pig_data[pig_name] = {
            "time_points": np.array(time_points),
            "naoh_conc":   np.array(naoh_conc),
            "ph_matrix":   np.array(ph_matrix),
        }
        if verbose and naoh_conc:
            print(f"  {pig_name}: {len(time_points)} times, {len(naoh_conc)} conc × pH")

    if verbose:
        print(f"✓ Data extraction complete for {len(all_pig_data)} pigs\n")
    return all_pig_data


def extract_alkalinity_data(filepath, sheet_name, verbose=True):
    """
    Extract citric acid titration data from an Excel ALKALINITY sheet (columns A–K).

    Same layout as ACIDITY: column B = citric concentration, columns C–K = pH per time.

    Parameters
    ----------
    filepath : str or pathlib.Path
    sheet_name : str
        Sheet name (typically 'ALKALINITY').
    verbose : bool

    Returns
    -------
    dict
        'PIG-XX' → {'time_points', 'citric_conc', 'ph_matrix'}.
    """
    if verbose:
        print("=" * 80)
        print("EXTRACTING CITRIC ACID DATA FROM ALKALINITY SHEET")
        print("=" * 80)

    df      = pd.read_excel(filepath, sheet_name=sheet_name)
    df_cols = df.iloc[:, :11]

    pig_starts = [
        (str(val), idx)
        for idx, val in enumerate(df_cols.iloc[:, 0])
        if pd.notna(val) and "PIG-" in str(val)
    ]
    if verbose:
        print(f"\nFound {len(pig_starts)} pigs in ALKALINITY sheet")

    all_pig_data = {}
    for pig_name, start_row in pig_starts:
        time_points = [
            float(val)
            for val in df_cols.iloc[start_row, 2:11]
            if pd.notna(val) and str(val) != "min" and _can_float(val)
        ]
        citric_conc = []
        ph_matrix   = []
        row_idx     = start_row + 2

        while row_idx < len(df_cols):
            citric_val = df_cols.iloc[row_idx, 1]
            if pd.isna(citric_val):
                break
            try:
                citric_conc.append(float(citric_val))
                ph_row = [
                    float(df_cols.iloc[row_idx, c]) if pd.notna(df_cols.iloc[row_idx, c]) else np.nan
                    for c in range(2, 11)
                ]
                ph_matrix.append(ph_row)
                row_idx += 1
            except (ValueError, TypeError):
                break

        all_pig_data[pig_name] = {
            "time_points": np.array(time_points),
            "citric_conc": np.array(citric_conc),
            "ph_matrix":   np.array(ph_matrix),
        }

    if verbose:
        print(f"✓ Citric data extraction complete for {len(all_pig_data)} pigs\n")
    return all_pig_data


def _can_float(val):
    """Return True if val can be converted to float."""
    try:
        float(val)
        return True
    except (ValueError, TypeError):
        return False


def interpolate_naoh_at_pH(naoh_conc, ph_values, target_pH):
    """
    Use PCHIP monotonic cubic spline to find NaOH concentration at target pH.

    Parameters
    ----------
    naoh_conc : array-like
    ph_values : array-like
    target_pH : float

    Returns
    -------
    tuple
        (naoh_at_target: float, success: bool). Returns (nan, False) on failure.
    """
    valid = ~np.isnan(ph_values)
    naoh  = naoh_conc[valid]
    ph    = ph_values[valid]

    if len(naoh) < 2 or target_pH < ph.min() or target_pH > ph.max():
        return np.nan, False

    try:
        interp = PchipInterpolator(ph, naoh)
        return float(interp(target_pH)), True
    except Exception:
        return np.nan, False


def interpolate_citric_at_pH4(citric_conc, ph_values, target_pH=4.0):
    """
    Linear interpolation to find citric acid concentration at target pH.

    pH decreases as citric acid is added, so the series is sorted descending.

    Parameters
    ----------
    citric_conc : array-like
    ph_values : array-like
    target_pH : float

    Returns
    -------
    tuple
        (citric_at_target: float, success: bool). Returns (nan, False) on failure.
    """
    valid  = ~np.isnan(ph_values)
    citric = np.asarray(citric_conc)[valid]
    ph     = np.asarray(ph_values)[valid]

    if len(citric) < 2 or target_pH < ph.min() or target_pH > ph.max():
        return np.nan, False

    try:
        sort_idx     = np.argsort(ph)[::-1]
        ph_sorted    = ph[sort_idx]
        citric_sorted = citric[sort_idx]
        unique_mask  = np.concatenate([[True], ph_sorted[1:] != ph_sorted[:-1]])
        ph_unique    = ph_sorted[unique_mask]
        citric_unique = citric_sorted[unique_mask]
        if len(ph_unique) < 2:
            return np.nan, False
        interp = interp1d(ph_unique, citric_unique, kind="linear")
        return float(interp(target_pH)), True
    except Exception:
        return np.nan, False


def process_all_pigs(pig_data, target_pH, verbose=True):
    """
    Apply PCHIP NaOH interpolation at target pH for every pig and time point.

    Parameters
    ----------
    pig_data : dict
        Output of extract_titration_data.
    target_pH : float
    verbose : bool

    Returns
    -------
    tuple
        (results, params_list) where results maps pig_name → {'time', 'naoh_pH11'}.
    """
    if verbose:
        print("=" * 80)
        print("NaOH SPLINE INTERPOLATION (PCHIP)")
        print("=" * 80)
        print(f"Target pH: {target_pH}\n")

    results       = {}
    params_list   = []
    success_count = 0
    total_count   = 0

    for pig_name in sorted(pig_data.keys()):
        data      = pig_data[pig_name]
        times     = data["time_points"]
        naoh_conc = data["naoh_conc"]
        ph_matrix = data["ph_matrix"]

        naoh_series = []
        for t_idx, time in enumerate(times):
            ph_vals = ph_matrix[:, t_idx]
            total_count += 1
            naoh_val, success = interpolate_naoh_at_pH(naoh_conc, ph_vals, target_pH)
            naoh_series.append(naoh_val)

            valid   = ~np.isnan(ph_vals)
            n_pts   = int(np.sum(valid))
            ph_rng  = (f"{ph_vals[valid].min():.2f}-{ph_vals[valid].max():.2f}"
                       if n_pts > 0 else "N/A")

            if success:
                success_count += 1
                if verbose:
                    print(f"  {pig_name} t={time:.0f} min: NaOH={naoh_val:.6e} mol/L")
            params_list.append({
                "Pig":                  pig_name,
                "Time (min)":           time,
                "N data points":        n_pts,
                "pH range":             ph_rng,
                "NaOH @ pH=11 (mol/L)": naoh_val if success else np.nan,
                "Status":               "Success" if success else "Failed - pH out of range",
            })

        results[pig_name] = {"time": times, "naoh_pH11": np.array(naoh_series)}

    if verbose:
        pct = 100 * success_count / max(1, total_count)
        print(f"\nSUMMARY: {success_count}/{total_count} successful ({pct:.1f}%)\n")
    return results, params_list


def process_citric_all_pigs(pig_data, target_pH, verbose=True):
    """
    Apply linear citric acid interpolation at target pH for every pig and time point.

    Parameters
    ----------
    pig_data : dict
        Output of extract_alkalinity_data.
    target_pH : float
    verbose : bool

    Returns
    -------
    dict
        pig_name → {'time', 'citric_pH4'}.
    """
    if verbose:
        print("=" * 80)
        print("CITRIC ACID INTERPOLATION (Linear)")
        print("=" * 80)
        print(f"Target pH: {target_pH}\n")

    results       = {}
    success_count = 0
    total_count   = 0

    for pig_name in sorted(pig_data.keys()):
        data        = pig_data[pig_name]
        times       = data["time_points"]
        citric_conc = data["citric_conc"]
        ph_matrix   = data["ph_matrix"]

        citric_series = []
        for t_idx in range(len(times)):
            ph_vals = ph_matrix[:, t_idx]
            total_count += 1
            citric_val, success = interpolate_citric_at_pH4(citric_conc, ph_vals, target_pH)
            citric_series.append(citric_val)
            if success:
                success_count += 1

        results[pig_name] = {"time": times, "citric_pH4": np.array(citric_series)}

    if verbose:
        pct = 100 * success_count / max(1, total_count)
        print(f"SUMMARY: {success_count}/{total_count} successful ({pct:.1f}%)\n")
    return results


def save_results(results, params_list, output_naoh, output_full, verbose=True):
    """
    Save NaOH interpolation results to two CSV files.

    Parameters
    ----------
    results : dict
        Output of process_all_pigs.
    params_list : list of dict
    output_naoh : str
        Simplified CSV (Pig, Time, NaOH value).
    output_full : str
        Full parameter table CSV.
    verbose : bool

    Returns
    -------
    tuple
        (simple_df, params_df).
    """
    if verbose:
        print("=" * 80)
        print("SAVING NaOH RESULTS")
        print("=" * 80)

    simple_data = [
        {"Pig": pig_name, "Time (min)": t, "NaOH @ pH=11 (mol/L)": naoh}
        for pig_name in sorted(results.keys())
        for t, naoh in zip(results[pig_name]["time"], results[pig_name]["naoh_pH11"])
        if not np.isnan(naoh)
    ]
    simple_df = pd.DataFrame(simple_data)
    ensure_parent_dir(output_naoh)
    simple_df.to_csv(output_naoh, index=False)

    params_df = pd.DataFrame(params_list)
    ensure_parent_dir(output_full)
    params_df.to_csv(output_full, index=False)

    if verbose:
        print(f"✓ {output_naoh}: {len(simple_df)} rows")
        print(f"✓ {output_full}: {len(params_df)} rows")
        success_df = params_df[params_df["Status"] == "Success"]
        if len(success_df) > 0:
            col = "NaOH @ pH=11 (mol/L)"
            print(f"  Mean: {success_df[col].mean():.6e} mol/L\n")
    return simple_df, params_df


def save_citric_results(results, output_file, verbose=True):
    """
    Save citric acid @ pH=4 interpolation results to CSV.

    Parameters
    ----------
    results : dict
        Output of process_citric_all_pigs.
    output_file : str
    verbose : bool

    Returns
    -------
    pandas.DataFrame
    """
    simple_data = [
        {"Pig": pig_name, "Time (min)": t, "Citric @ pH=4": citric}
        for pig_name in sorted(results.keys())
        for t, citric in zip(results[pig_name]["time"], results[pig_name]["citric_pH4"])
        if not np.isnan(citric)
    ]
    df = pd.DataFrame(simple_data)
    ensure_parent_dir(output_file)
    df.to_csv(output_file, index=False)
    if verbose:
        print(f"✓ {output_file}: {len(df)} rows\n")
    return df


def create_group_plots(results, healthy_pigs, nonhealthy_pigs, output_file=None, display=False):
    """
    Side-by-side plots of NaOH @ pH=11 for healthy vs non-healthy pigs over time.

    Parameters
    ----------
    results : dict
        Output of process_all_pigs.
    healthy_pigs : list of str
    nonhealthy_pigs : list of str
    output_file : str or None
        Save path; if None the figure is not saved.
    display : bool
        Call plt.show() for inline notebook display.

    Returns
    -------
    matplotlib.figure.Figure
    """
    healthy_colors    = ["#FFD700", "#4682B4", "#4682B4", "#4682B4"]
    nonhealthy_colors = ["#808080", "#FF6347", "#90EE90", "#8A2BE2", "#FF6347", "#FFD700"]

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

    for ax, pig_list, colors, panel_title in [
        (ax1, healthy_pigs,    healthy_colors,    "HEALTHY CONTROLS"),
        (ax2, nonhealthy_pigs, nonhealthy_colors, "NON-HEALTHY CONTROLS"),
    ]:
        for idx, pig in enumerate(pig_list):
            if pig not in results:
                continue
            data   = results[pig]
            times  = data["time"]
            naoh   = data["naoh_pH11"]
            valid  = ~np.isnan(naoh)
            t_v, n_v = times[valid], naoh[valid]
            color  = colors[idx % len(colors)]
            ax.plot(t_v, n_v, "o",  color=color, markersize=6, alpha=0.7)
            ax.plot(t_v, n_v, "--", color=color, linewidth=1.5, alpha=0.5)
            if len(t_v) > 0:
                ax.text(t_v[-1] + 5, n_v[-1], pig, fontsize=10, fontweight="bold", va="center")

        ax.set_title(panel_title, fontsize=14, fontweight="bold", pad=15,
                     bbox=dict(boxstyle="round", facecolor="white", edgecolor="black", linewidth=2))
        ax.set_xlabel("Time (min)", fontsize=12)
        ax.set_ylabel("NaOH (mol to pH=11)", fontsize=12)
        ax.set_xlim([0, 400])
        ax.set_ylim([0, 1.6e-4])
        ax.ticklabel_format(style="scientific", axis="y", scilimits=(0, 0))
        ax.grid(True, alpha=0.3)

    plt.suptitle("Acidity Titration – 10 Control Pigs", fontsize=16, fontweight="bold", y=0.98)
    plt.tight_layout()

    if output_file:
        ensure_parent_dir(output_file)
        fig.savefig(output_file, dpi=300, bbox_inches="tight")
        print(f"✓ Plot saved to: {output_file}")
    if display:
        plt.show()
    else:
        plt.close()
    return fig


# ============================================================================
# ANALYSIS PIPELINES
# ============================================================================

def analyze_phosphate_creatinine(input_file, phosphate_creatinine_figure_path):
    """
    Plot by-pig phosphate/creatinine boxplots and print group statistics.

    Parameters
    ----------
    input_file : str or pathlib.Path
        Titration workbook with 'PHOSPHATE and CREATININE' sheet.
    phosphate_creatinine_figure_path : str or pathlib.Path
        Output PNG path.

    Returns
    -------
    dict
        {'phosphate': stats_dict, 'creatinine': stats_dict}.
    """
    print("\n" + "=" * 100)
    print("PHOSPHATE AND CREATININE ANALYSIS")
    print("=" * 100)

    phosphate_data, creatinine_data = read_phosphate_creatinine_by_pig(input_file)

    plot_biomarkers_by_pig_figure(
        [
            ("A", "Phosphate Distribution by Pig",  "Phosphate concentration (mg/dL)",  phosphate_data),
            ("B", "Creatinine Distribution by Pig", "Creatinine concentration (mg/dL)", creatinine_data),
        ],
        "Phosphate and Creatinine by Pig: Healthy vs Non-Healthy",
        "Green = Healthy | Pink = Non-Healthy",
        phosphate_creatinine_figure_path,
    )

    healthy_phosphate    = np.concatenate([phosphate_data[p]    for p in HEALTHY_PIGS])
    nonhealthy_phosphate = np.concatenate([phosphate_data[p]    for p in NONHEALTHY_PIGS])
    healthy_creatinine   = np.concatenate([creatinine_data[p]   for p in HEALTHY_PIGS])
    nonhealthy_creatinine = np.concatenate([creatinine_data[p]  for p in NONHEALTHY_PIGS])

    phosphate_stats  = calculate_statistics(healthy_phosphate,   nonhealthy_phosphate)
    creatinine_stats = calculate_statistics(healthy_creatinine,  nonhealthy_creatinine)

    for name, stats in [("Phosphate", phosphate_stats), ("Creatinine", creatinine_stats)]:
        print(f"\n{name}:")
        print(f"  Healthy:     {stats['healthy_mean']:.2f} ± {stats['healthy_std']:.2f}")
        print(f"  Non-Healthy: {stats['nonhealthy_mean']:.2f} ± {stats['nonhealthy_std']:.2f}")
        print(f"  p-value:     {stats['p_value_ttest']:.6f} ({stats['significance_status']})")
        print(f"  Cohen's d:   {stats['cohens_d']:.3f} ({stats['effect_size']})")

    return {"phosphate": phosphate_stats, "creatinine": creatinine_stats}


def load_naoh_from_csv(csv_path):
    """
    Load NaOH @ pH=11 statistics from naoh_pH11_results.csv.

    The CSV is produced by spline_interpolation; this function does not repeat
    the interpolation logic.

    Parameters
    ----------
    csv_path : str or pathlib.Path

    Returns
    -------
    dict
        {'naoh_pH11': stats_dict} or {} if file is missing or malformed.
    """
    print("\n" + "=" * 100)
    print("NaOH @ pH=11 (from spline_interpolation output)")
    print("=" * 100)

    try:
        df = pd.read_csv(csv_path)
    except FileNotFoundError:
        print(f"  WARNING: {csv_path} not found. Run spline_interpolation first.")
        return {}

    col_naoh = "NaOH @ pH=11 (mol/L)"
    if col_naoh not in df.columns or "Pig" not in df.columns:
        print(f"  WARNING: Expected columns 'Pig', '{col_naoh}' in {csv_path}")
        return {}

    healthy_naoh    = df[df["Pig"].isin(HEALTHY_PIGS)][col_naoh].dropna().values
    nonhealthy_naoh = df[df["Pig"].isin(NONHEALTHY_PIGS)][col_naoh].dropna().values

    if len(healthy_naoh) < 2 or len(nonhealthy_naoh) < 2:
        print("  WARNING: Insufficient data for statistical comparison")
        return {}

    stats = calculate_statistics(healthy_naoh, nonhealthy_naoh)
    print(f"\nNaOH @ pH=11:")
    print(f"  Healthy:     {stats['healthy_mean']:.6e} ± {stats['healthy_std']:.6e} mol/L")
    print(f"  Non-Healthy: {stats['nonhealthy_mean']:.6e} ± {stats['nonhealthy_std']:.6e} mol/L")
    print(f"  p-value:     {stats['p_value_ttest']:.6f} ({stats['significance_status']})")
    print(f"  Cohen's d:   {stats['cohens_d']:.3f} ({stats['effect_size']})")
    return {"naoh_pH11": stats}


def load_citric_from_csv(csv_path):
    """
    Load Citric @ pH=4 statistics from a citric spline CSV.

    The CSV is produced by spline_interpolation; this function does not repeat
    the interpolation logic.

    Parameters
    ----------
    csv_path : str or pathlib.Path

    Returns
    -------
    dict
        {'citric_pH4': stats_dict} or {} if file is missing or malformed.
    """
    print("\n" + "=" * 100)
    print("Citric @ pH=4 (from spline_interpolation output)")
    print("=" * 100)

    try:
        df = pd.read_csv(csv_path)
    except FileNotFoundError:
        print(f"  WARNING: {csv_path} not found. Run spline_interpolation first.")
        return {}

    if "Pig" not in df.columns:
        print(f"  WARNING: Expected column 'Pig' in {csv_path}")
        return {}

    citric_col = next(
        (c for c in df.columns if str(c).strip().startswith("Citric @ pH=4")),
        None,
    )
    if citric_col is None:
        print(f"  WARNING: No 'Citric @ pH=4' column found in {csv_path}")
        return {}

    healthy_citric    = df[df["Pig"].isin(HEALTHY_PIGS)][citric_col].dropna().values
    nonhealthy_citric = df[df["Pig"].isin(NONHEALTHY_PIGS)][citric_col].dropna().values

    if len(healthy_citric) < 2 or len(nonhealthy_citric) < 2:
        print("  WARNING: Insufficient data for Citric @ pH=4 statistical comparison")
        return {}

    stats = calculate_statistics(healthy_citric, nonhealthy_citric)
    print(f"\nCitric @ pH=4:")
    print(f"  Healthy:     {stats['healthy_mean']:.4f} ± {stats['healthy_std']:.4f}")
    print(f"  Non-Healthy: {stats['nonhealthy_mean']:.4f} ± {stats['nonhealthy_std']:.4f}")
    print(f"  p-value:     {stats['p_value_ttest']:.6f} ({stats['significance_status']})")
    print(f"  Cohen's d:   {stats['cohens_d']:.3f} ({stats['effect_size']})")
    return {"citric_pH4": stats}


def analyze_direct_pH(input_file):
    """
    Analyze direct urine pH measurements from the titration workbook.

    Parameters
    ----------
    input_file : str or pathlib.Path

    Returns
    -------
    dict
        {'urine_pH': stats_dict}.
    """
    print("\n" + "=" * 100)
    print("DIRECT URINE pH ANALYSIS")
    print("=" * 100)

    df = pd.read_excel(input_file, sheet_name="pHi, OSM, Vol rate", nrows=18)

    healthy_pH    = np.array([v for p in HEALTHY_PIGS    if p in df.columns
                               for v in df[p].dropna().values])
    nonhealthy_pH = np.array([v for p in NONHEALTHY_PIGS if p in df.columns
                               for v in df[p].dropna().values])

    stats = calculate_statistics(healthy_pH, nonhealthy_pH)
    print(f"\nUrine pH:")
    print(f"  Healthy:     {stats['healthy_mean']:.2f} ± {stats['healthy_std']:.2f}")
    print(f"  Non-Healthy: {stats['nonhealthy_mean']:.2f} ± {stats['nonhealthy_std']:.2f}")
    print(f"  p-value:     {stats['p_value_ttest']:.6f} ({stats['significance_status']})")
    print(f"  Cohen's d:   {stats['cohens_d']:.3f} ({stats['effect_size']})")
    return {"urine_pH": stats}


def analyze_ammonia_co2_ratio_boxplots(ammonia_csv, co2_csv, ratio_csv, boxplot_output_path):
    """
    Healthy vs non-healthy comparison for ammonia, CO₂, and ammonia/CO₂ ratio.

    Loads CSVs from ammonia_analysis, pools all time points per group, computes
    statistics, and saves a pooled 1×3 boxplot figure.

    Parameters
    ----------
    ammonia_csv : str or pathlib.Path
    co2_csv     : str or pathlib.Path
    ratio_csv   : str or pathlib.Path
    boxplot_output_path : str or pathlib.Path
        Output PNG for the pooled group comparison.

    Returns
    -------
    dict
        Analyte name → stats_dict for analytes with sufficient data.
    """
    print("\n" + "=" * 100)
    print("AMMONIA, CO₂, AND AMMONIA/CO₂ RATIO (from ammonia_analysis CSVs)")
    print("=" * 100)

    required = [("ammonia", ammonia_csv), ("CO2", co2_csv), ("ammonia/CO2 ratio", ratio_csv)]
    missing  = [label for label, p in required if not Path(p).exists()]
    if missing:
        print("  WARNING: Missing CSV file(s):")
        for label, p in required:
            if not Path(p).exists():
                print(f"    - {p} ({label})")
        print("  Run ammonia_analysis.ipynb through the save step first.")
        return {}

    ammonia_df = pd.read_csv(ammonia_csv)
    co2_df     = pd.read_csv(co2_csv)
    ratio_df   = pd.read_csv(ratio_csv)

    analytes = [
        ("Ammonia (mM)",     ammonia_df, "Ammonia (mM)"),
        ("CO2 (mM)",         co2_df,     "Carbon dioxide (mM)"),
        ("Ammonia/CO2 ratio", ratio_df,  "Ammonia / CO2 ratio"),
    ]

    stats_by_name = {}
    plot_series   = []

    for display_name, df, y_axis_short in analytes:
        h_vals  = pooled_biomarker_values_for_pigs(df, HEALTHY_PIGS)
        nh_vals = pooled_biomarker_values_for_pigs(df, NONHEALTHY_PIGS)
        if h_vals.size < 2 or nh_vals.size < 2:
            print(f"\n{display_name}: insufficient pooled data for tests; skipping.")
            continue
        st = calculate_statistics(h_vals, nh_vals)
        stats_by_name[display_name] = st
        plot_series.append((display_name, h_vals, nh_vals, y_axis_short))

        print(f"\n{display_name}:")
        print(f"  Healthy:     {st['healthy_mean']:.4g} ± {st['healthy_std']:.4g} (n={h_vals.size})")
        print(f"  Non-Healthy: {st['nonhealthy_mean']:.4g} ± {st['nonhealthy_std']:.4g} (n={nh_vals.size})")
        print(f"  p-value:     {st['p_value_ttest']:.6f} ({st['significance_status']})")
        print(f"  Cohen's d:   {st['cohens_d']:.3f} ({st['effect_size']})")

    if not plot_series:
        print("  No figure written (no analyte passed data checks).")
        return stats_by_name

    n_panels = len(plot_series)
    fig, axes = plt.subplots(1, n_panels, figsize=(5 * n_panels, 5.5))
    if n_panels == 1:
        axes = [axes]

    for ax, (display_name, h_vals, nh_vals, y_axis_short) in zip(axes, plot_series):
        bp = ax.boxplot(
            [h_vals, nh_vals],
            labels=["Healthy", "Non-Healthy"],
            patch_artist=True,
            widths=0.55,
            showfliers=True,
            medianprops=dict(color="black", linewidth=2),
            boxprops=dict(linewidth=1.5),
        )
        bp["boxes"][0].set_facecolor(HEALTHY_COLOR)
        bp["boxes"][1].set_facecolor(NONHEALTHY_COLOR)
        ax.set_ylabel(y_axis_short, fontsize=11, fontweight="bold")
        ax.set_title(display_name,  fontsize=12, fontweight="bold", pad=8)
        ax.grid(True, alpha=0.3, axis="y")

    fig.suptitle(
        "Healthy vs Non-Healthy: Ammonia, CO₂, and Ammonia/CO₂ ratio (all time points pooled)",
        fontsize=13, fontweight="bold", y=1.02,
    )
    plt.tight_layout()
    ensure_parent_dir(boxplot_output_path)
    plt.savefig(boxplot_output_path, dpi=300, bbox_inches="tight")
    plt.close()
    print(f"\n✓ Saved: {boxplot_output_path}")

    return stats_by_name


def plot_ammonia_co2_ratio_by_pig_figure(ammonia_csv, co2_csv, ratio_csv, output_path):
    """
    Save a 1×3 by-pig boxplot figure for ammonia, CO₂, and ratio.

    Parameters
    ----------
    ammonia_csv : str or pathlib.Path
    co2_csv     : str or pathlib.Path
    ratio_csv   : str or pathlib.Path
    output_path : str or pathlib.Path
    """
    paths = [Path(ammonia_csv), Path(co2_csv), Path(ratio_csv)]
    if any(not p.is_file() for p in paths):
        return

    ammonia_df = pd.read_csv(ammonia_csv)
    co2_df     = pd.read_csv(co2_csv)
    ratio_df   = pd.read_csv(ratio_csv)

    panels = [
        ("A", "Ammonia Distribution by Pig",          "Ammonia concentration (mM)",       master_wide_df_to_pig_dict(ammonia_df)),
        ("B", "Carbon Dioxide Distribution by Pig",   "Carbon dioxide concentration (mM)", master_wide_df_to_pig_dict(co2_df)),
        ("C", "Ammonia / CO2 Ratio Distribution by Pig", "Ammonia / CO2 ratio (log scale)", master_wide_df_to_pig_dict(ratio_df)),
    ]
    plot_biomarkers_by_pig_figure(
        panels,
        "Ammonia, CO₂, and Ammonia/CO₂ Ratio by Pig: Healthy vs Non-Healthy",
        "Green = Healthy | Pink = Non-Healthy",
        output_path,
        y_axis_scales=["linear", "linear", "log"],
    )


def plot_naoh_citric_by_pig_figure(naoh_csv, citric_csv, output_path):
    """
    Save a 1×2 by-pig boxplot figure for NaOH @ pH=11 and Citric @ pH=4.

    CSV values are in mol/L; displayed values are multiplied by 1000 (mM).

    Parameters
    ----------
    naoh_csv   : str or pathlib.Path
    citric_csv : str or pathlib.Path
    output_path : str or pathlib.Path
    """
    naoh_path   = Path(naoh_csv)
    citric_path = Path(citric_csv)
    if not naoh_path.is_file():
        print(f"  WARNING: NaOH CSV missing: {naoh_path}")
        return
    if not citric_path.is_file():
        print(f"  WARNING: Citric CSV missing: {citric_path}")
        return

    naoh_dict   = scale_pig_series_dict(pig_dict_from_spline_long_csv(naoh_path,   "NaOH @ pH=11"),  1000.0)
    citric_dict = scale_pig_series_dict(pig_dict_from_spline_long_csv(citric_path, "Citric @ pH=4"), 1000.0)

    plot_biomarkers_by_pig_figure(
        [
            ("A", "NaOH @ pH=11 Distribution by Pig",  "NaOH @ pH=11 (mM/L)",  naoh_dict),
            ("B", "Citric @ pH=4 Distribution by Pig", "Citric @ pH=4 (mM/L)", citric_dict),
        ],
        "NaOH @ pH=11 and Citric @ pH=4 by Pig: Healthy vs Non-Healthy",
        "Green = Healthy | Pink = Non-Healthy",
        output_path,
    )


def plot_urine_ph_by_pig_figure(input_file, output_path):
    """
    Save a by-pig boxplot figure for direct urine pH.

    Parameters
    ----------
    input_file  : str or pathlib.Path
    output_path : str or pathlib.Path
    """
    pig_dict = load_urine_ph_by_pig(input_file)
    plot_biomarkers_by_pig_figure(
        [("A", "Urine pH Distribution by Pig", "Urine pH", pig_dict)],
        "Urine pH by Pig: Healthy vs Non-Healthy",
        "Green = Healthy | Pink = Non-Healthy",
        output_path,
    )


def create_biomarker_summary(all_stats, output_file):
    """
    Build and save the biomarker summary CSV table.

    Parameters
    ----------
    all_stats : dict
        Biomarker name → stats dict from calculate_statistics.
    output_file : str

    Returns
    -------
    pandas.DataFrame
    """
    print("\n" + "=" * 100)
    print("CREATING BIOMARKER SUMMARY TABLE")
    print("=" * 100)

    rows = []
    for biomarker_name, stats in all_stats.items():
        if not stats:
            continue
        p = stats["p_value_ttest"]
        rows.append({
            "Biomarker":           biomarker_name,
            "Healthy Mean":        f"{stats['healthy_mean']:.2f}",
            "Non-Healthy Mean":    f"{stats['nonhealthy_mean']:.2f}",
            "Difference (%)":      f"{stats['diff_percent']:+.1f}%",
            "p-value":             f"{p:.3f}" if p >= 0.001 else "<0.001",
            "Significance Status": stats["significance_status"],
            "Cohen's d":           f"{stats['cohens_d']:.2f}",
            "Effect Size":         stats["effect_size"],
        })

    df = pd.DataFrame(rows)
    ensure_parent_dir(output_file)
    df.to_csv(output_file, index=False)
    print(f"\n✓ Biomarker summary saved to: {output_file}")
    print(f"\n{df.to_string(index=False)}")
    return df


def resolve_citric_csv_path(base_dir, primary_filename="citric_acid_pH4_results.csv",
                             alternate_filenames=("citric_pH4_results.csv",)):
    """
    Return the path to a citric spline CSV, trying primary then alternates.

    Parameters
    ----------
    base_dir : str or pathlib.Path
    primary_filename : str
    alternate_filenames : tuple of str

    Returns
    -------
    str
        First existing path found, otherwise the primary path (even if absent).
    """
    base = Path(base_dir)
    for name in (primary_filename,) + tuple(alternate_filenames):
        p = base / name
        if p.is_file():
            return str(p)
    return str(base / primary_filename)


def run_full_biomarker_analysis(
    output_dir=None,
    titration_xlsx="Titration_Data_02052026.xlsx",
    project_root=None,
    naoh_csv_name="naoh_pH11_results.csv",
    citric_primary="citric_acid_pH4_results.csv",
    citric_alternates=("citric_pH4_results.csv",),
    ammonia_csv="ammonia_concentration.csv",
    co2_csv="co2_concentration.csv",
    ratio_csv="ammonia_co2_ratio.csv",
    ammonia_pooled_boxplot_png="ammonia_co2_ratio_boxplots.png",
    phosphate_creatinine_by_pig_png="phosphate_creatinine_boxplots.png",
    ammonia_by_pig_png="ammonia_co2_ratio_by_pig_boxplots.png",
    naoh_citric_by_pig_png="naoh_citric_by_pig_boxplots.png",
    urine_ph_by_pig_png="urine_ph_by_pig_boxplots.png",
    biomarker_summary_csv="biomarker_summary.csv",
):
    """
    Run the full biomarker pipeline: statistics, figures, and summary CSV.

    Steps:
      1. Phosphate & creatinine (from titration workbook)
      2. NaOH @ pH=11 and Citric @ pH=4 (from spline CSV outputs)
      3. Urine pH (from titration workbook)
      4. Ammonia, CO₂, and ratio (from ammonia_analysis CSV outputs)
      5. Save biomarker_summary.csv

    Parameters
    ----------
    output_dir : str or pathlib.Path, optional
        Directory for CSV/PNG outputs. Defaults to ``DEFAULT_OUTPUT_DIR``
        (``catalyze/output``).
    titration_xlsx : str or pathlib.Path
        Excel workbook with phosphate, creatinine, and pH sheets. If relative,
        resolved under ``project_root`` (not under ``output_dir``).
    project_root : str or pathlib.Path, optional
        Catalyze project folder for resolving relative ``titration_xlsx``.
        Defaults to ``CATALYZE_DIR``.
    Remaining parameters : str
        Output filenames resolved under output_dir.

    Returns
    -------
    pandas.DataFrame
        Biomarker summary table.
    """
    print("\n" + "=" * 100)
    print("NIH CATALYZE SMART CATHETER PROJECT - BIOMARKER ANALYSIS")
    print("=" * 100)

    if output_dir is None:
        output_dir = DEFAULT_OUTPUT_DIR
    if project_root is None:
        project_root = CATALYZE_DIR

    base = Path(output_dir)
    base.mkdir(parents=True, exist_ok=True)

    titration_path = Path(titration_xlsx)
    if not titration_path.is_absolute():
        titration_path = Path(project_root) / titration_xlsx

    all_stats = {}

    pc_stats = analyze_phosphate_creatinine(titration_path, base / phosphate_creatinine_by_pig_png)
    all_stats["Phosphate"] = pc_stats["phosphate"]
    all_stats["Creatinine"] = pc_stats["creatinine"]

    naoh_csv = str(base / naoh_csv_name)
    naoh_stats = load_naoh_from_csv(naoh_csv)
    if naoh_stats:
        all_stats["NaOH @ pH=11"] = naoh_stats["naoh_pH11"]

    citric_csv = resolve_citric_csv_path(base, citric_primary, citric_alternates)
    citric_stats = load_citric_from_csv(citric_csv)
    if citric_stats:
        all_stats["Citric @ pH=4"] = citric_stats["citric_pH4"]
    plot_naoh_citric_by_pig_figure(naoh_csv, citric_csv, base / naoh_citric_by_pig_png)

    ph_stats = analyze_direct_pH(titration_path)
    all_stats["Urine pH"] = ph_stats["urine_pH"]
    plot_urine_ph_by_pig_figure(titration_path, base / urine_ph_by_pig_png)

    ammo_ratio_stats = analyze_ammonia_co2_ratio_boxplots(
        str(base / ammonia_csv), str(base / co2_csv), str(base / ratio_csv),
        str(base / ammonia_pooled_boxplot_png),
    )
    all_stats.update(ammo_ratio_stats)
    plot_ammonia_co2_ratio_by_pig_figure(
        str(base / ammonia_csv), str(base / co2_csv), str(base / ratio_csv),
        base / ammonia_by_pig_png,
    )

    summary_df = create_biomarker_summary(all_stats, str(base / biomarker_summary_csv))

    print("\n" + "=" * 100)
    print("ANALYSIS COMPLETE")
    print("=" * 100)
    print("\nKEY FINDINGS:")
    print("  1. Phosphate shows HIGHLY SIGNIFICANT difference (primary biomarker)")
    print("  2. Urine pH shows NO difference (baseline preserved)")
    print("  3. Conclusion: Specific tubular dysfunction, not generalized failure")

    return summary_df