csu_core.py

"""Core utilities for the shuffled-vs-unshuffled comparison pipeline.

This file is auto-refactored out of the original `compare_shuffled_unshuffled.py` to improve readability.
It contains:
- Robust statistical helpers
- CSV probing helpers
- Plot saving utilities (Plotly HTML/PNG/EPS + an index page)
- Generic helpers used by feature extraction, pipeline, and mixed-model modules

Public API:
- compare_shuffled_unshuffled
- compare_participant_metrics
"""

from __future__ import annotations

import os
import re
import warnings
import math
from fnmatch import fnmatch
from typing import Optional, List, Tuple, Any, Dict

import numpy as np
from numpy.typing import NDArray
import pandas as pd
from scipy.stats import ttest_ind, linregress, pearsonr, mannwhitneyu

import plotly.graph_objects as go  # noqa:F401

import statsmodels.formula.api as smf  # noqa:F401
import statsmodels.api as sm  # noqa:F401

# Project helpers (used for plot output directories).
from helper import HMD_helper, apply_global_plotly_style
import common  # project module used by HMD_helper.save_plotly
from custom_logger import CustomLogger

# Shared yaw/quaternion column heuristics.
# These are used by file-selection helpers (e.g., `_score_trial_file`) and yaw feature extraction.
# Kept in a standalone module to avoid circular imports between `csu_core` and `csu_features`.
from csu_yaw_constants import _YAW_CANDIDATES, _QUAT_REGEX, _QUAT_LIST_COL_PAT  # noqa:F401

HAVE_SM = True
logger = CustomLogger(__name__)  # use custom logger


# ---------------------------------------------------------------------------
# Output precision helpers
# ---------------------------------------------------------------------------
# Keep raw numeric values in DataFrames/CSVs, but format console/report tables
# with at most three decimal places. Small non-zero values are shown in
# scientific notation so p values do not appear as exact zero.
DISPLAY_DECIMALS = 3
font_size = common.get_configs("font_size")


def _format_float_for_display(value: float, decimals: int = DISPLAY_DECIMALS) -> str:
    try:
        x = float(value)
    except Exception:
        return str(value)
    if np.isnan(x):
        return "nan"
    if np.isposinf(x):
        return "inf"
    if np.isneginf(x):
        return "-inf"
    if x != 0 and abs(x) < 10 ** (-decimals):
        return f"{x:.{decimals}e}"
    text = f"{x:.{decimals}f}".rstrip("0").rstrip(".")
    return "0" if text in {"-0", ""} else text


def _to_string_3dp(df: pd.DataFrame, *args, **kwargs) -> str:
    """Return a DataFrame string with floats shown to at most 3 decimals."""
    return df.to_string(*args, float_format=_format_float_for_display, **kwargs)


def _round_numeric_for_output(df: pd.DataFrame, decimals: int = DISPLAY_DECIMALS) -> pd.DataFrame:
    """Return a copy with numeric columns rounded for optional paper outputs."""
    out = df.copy()
    numeric_cols = out.select_dtypes(include=[np.number]).columns
    if len(numeric_cols):
        out.loc[:, numeric_cols] = out.loc[:, numeric_cols].round(decimals)
    return out


# Consistent dataset labels and colours used across all exported figures.
DATASET_LABEL_MAP = {
    "shuffled": "Randomised",
    "unshuffled": "Fixed order",
    "fixed-order": "Fixed order",
    "fixed order": "Fixed order",
    "randomized": "Randomised",
    "randomised": "Randomised",
}

DATASET_COLOUR_MAP = {
    "shuffled": "#1f77b4",
    "unshuffled": "#ff7f0e",
    "Randomised": "#1f77b4",
    "Fixed order": "#ff7f0e",
}

# Compatibility alias for earlier patches that used American spelling.
DATASET_COLOR_MAP = DATASET_COLOUR_MAP


# ---------------------------------------------------------------------------
# Per-figure text size controls for exported Plotly figures
# ---------------------------------------------------------------------------
# Edit these dictionaries when a specific figure needs larger or smaller text.
# The keys are the file names passed to _save_plot(..., name="...") without
# extension. Wildcards are supported through PLOT_TEXT_SIZE_BY_PATTERN.
#
# Available keys in a style dictionary:
#   x_label_size, y_label_size, x_tick_size, y_tick_size,
#   legend_size, legend_title_size, legend_title_text, legend_x, legend_y,
#   legend_xanchor, legend_yanchor, legend_orientation,
#   line_width, violin_line_width, marker_line_width, marker_size,
#   title_size, annotation_size, font_size, width, height
#
# Legend placement uses Plotly paper coordinates. For example,
# legend_x=0.98 and legend_y=0.98 places the legend near the top right,
# while legend_x=0.02 and legend_y=0.98 places it near the top left.
#
# These settings are applied centrally before HTML/PNG/EPS export, so all bar,
# line, scatter and violin plots can be controlled from one place.
PLOT_TEXT_SIZE_DEFAULT: Dict[str, Any] = {
    "font_size": font_size - 4,
    "x_label_size": font_size,
    "y_label_size": font_size,
    "x_tick_size": font_size - 4,
    "y_tick_size": font_size - 4,
    "legend_size": font_size - 4,
    "legend_title_size": font_size - 4,
    "title_size": font_size,
    "annotation_size": font_size - 4,
    "width": 1320,
    "height": 680,
}

PLOT_TEXT_SIZE_BY_KIND: Dict[str, Dict[str, Any]] = {
    "bar": {
        "x_label_size": 22,
        "y_label_size": 22,
        "x_tick_size": 18,
        "y_tick_size": 18,
        "legend_size": 18,
        "marker_line_width": 1,
        "height": 760,
    },
    "line": {
        "x_label_size": 24,
        "y_label_size": 24,
        "x_tick_size": 20,
        "y_tick_size": 20,
        "legend_size": 18,
        "line_width": 3,
        "height": 760,
    },
    "violin": {
        "x_label_size": 24,
        "y_label_size": 24,
        "x_tick_size": 20,
        "y_tick_size": 20,
        "legend_size": 18,
        "violin_line_width": 2,
        "height": 760,
    },
    "scatter": {
        "x_label_size": 22,
        "y_label_size": 22,
        "x_tick_size": 18,
        "y_tick_size": 18,
        "legend_size": 17,
        "height": 740,
    },
}

# Exact figure-specific overrides. Add any generated figure name here when you
# want that one plot to have different sizes from the default for its kind.
PLOT_TEXT_SIZE_BY_NAME: Dict[str, Dict[str, Any]] = {
    # Main manuscript line plots
    "curve_time_on_task_trigger_mean": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 34,
        "y_tick_size": 34,
        "legend_size": 34,
        "legend_x": 0.25,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "height": 820,
    },
    "curve_time_on_task_Q3": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 34,
        "y_tick_size": 34,
        "legend_size": 34,
        "legend_x": 0.25,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "height": 820,
    },
    "curve_time_on_task_dtrigger_sd": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 34,
        "y_tick_size": 34,
        "legend_size": 34,
        "legend_x": 0.25,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "height": 820,
    },

    # Appendix latency event absence line plots
    "missingness_press_over_trial": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 36,
        "legend_title_text": "",
        "legend_x": 0.98,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "height": 780,
    },
    "missingness_release_over_trial": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 36,
        "legend_title_text": "",
        "legend_x": 0.98,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "height": 780,
    },

    # AUC bar plot
    "F2_bar_auc_by_signal": {
        "x_label_size": 24,
        "y_label_size": 24,
        "x_tick_size": 18,
        "y_tick_size": 20,
        "legend_size": 18,
        "height": 820,
    },

    # Exposure based figure panels
    "MM5_curve_Q3_exposure_yielding": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 36,
        "legend_title_text": "",
        "legend_x": 0.25,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "height": 780,
    },
    "MM5_curve_Q3_exposure_eHMI": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 36,
        "legend_title_text": "",
        "legend_x": 0.85,
        "legend_y": 0.18,
        "legend_orientation": "v",
        "line_width": 4,
        "height": 780,
    },
    "MM5_forest_exposure_interactions": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "marker_size": 11,
        "width": 1500,
        "height": 720,
    },


    # Trial-position factor drift plots
    "factor_drift_yielding_over_trial_index": {
        "x_label_size": 28,
        "y_label_size": 28,
        "x_tick_size": 23,
        "y_tick_size": 23,
        "legend_size": 22,
        "legend_x": 0.85,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "marker_size": 8,
        "height": 820,
        "width": 1320,
    },
    "factor_drift_eHMIOn_over_trial_index": {
        "x_label_size": 28,
        "y_label_size": 28,
        "x_tick_size": 23,
        "y_tick_size": 23,
        "legend_size": 22,
        "legend_x": 0.8,
        "legend_y": 0.98,
        "legend_orientation": "v",
        "line_width": 4,
        "marker_size": 8,
        "height": 820,
        "width": 1320,
    },

    # Main manuscript participant-level learning and carryover violin plots
    "compare_participant_violin_E_carryover_Q3": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },
    "compare_participant_violin_E_drift_Q3": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },
    "compare_participant_violin_E_slope_Q3": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },
    "compare_participant_violin_E_carryover_dtrigger_sd": {
        "x_label_size": 28,
        "y_label_size": 28,
        "x_tick_size": 24,
        "y_tick_size": 24,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },
    "compare_participant_violin_E_carryover_trigger_mean": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },

    # Additional carryover factor violin plots
    "compare_participant_violin_E_carryover_prev_eHMIOn_Q3": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },
    "compare_participant_violin_E_carryover_prev_camera_Q3": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },
    "compare_participant_violin_E_carryover_prev_distPed_Q3": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "violin_line_width": 2.5,
        "marker_size": 8,
        "showlegend": False,
        "height": 820,
        "width": 1320,
    },

    # Reliability figures
    "reliability_trigger_mean_odd_even": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 30,
        "legend_title_text": "",
        "legend_x": 0.05,
        "legend_y": 0.95,
        "legend_xanchor": "left",
        "legend_yanchor": "top",
        "legend_orientation": "v",
        "marker_size": 12,
        "height": 760,
    },

    "reliability_trigger_mean_early_late": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 30,
        "legend_title_text": "",
        "legend_x": 0.05,
        "legend_y": 0.95,
        "legend_xanchor": "left",
        "legend_yanchor": "top",
        "legend_orientation": "v",
        "marker_size": 12,
        "height": 760,
    },

    "reliability_Q3_odd_even": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 30,
        "legend_title_text": "",
        "legend_x": 0.05,
        "legend_y": 0.95,
        "legend_xanchor": "left",
        "legend_yanchor": "top",
        "legend_orientation": "v",
        "marker_size": 12,
        "height": 760,
    },

    "reliability_Q3_early_late": {
        "x_label_size": 40,
        "y_label_size": 40,
        "x_tick_size": 36,
        "y_tick_size": 36,
        "legend_size": 30,
        "legend_title_text": "",
        "legend_x": 0.05,
        "legend_y": 0.95,
        "legend_xanchor": "left",
        "legend_yanchor": "top",
        "legend_orientation": "v",
        "marker_size": 12,
        "height": 760,
    },
}

# Wildcard overrides for groups of figures. Exact entries above still win.
PLOT_TEXT_SIZE_BY_PATTERN: List[Tuple[str, Dict[str, Any]]] = [
    ("compare_participant_violin_E_*", {"x_label_size": 26,
                                        "y_label_size": 26,
                                        "x_tick_size": 22,
                                        "y_tick_size": 22,
                                        "legend_size": 20,
                                        "legend_x": 0.98,
                                        "legend_y": 0.98,
                                        "violin_line_width": 2,
                                        "showlegend": False,
                                        "height": 800}),

    ("compare_participant_violin_breakmatched_*", {"x_label_size": 24,
                                                   "y_label_size": 24,
                                                   "x_tick_size": 20,
                                                   "y_tick_size": 20,
                                                   "showlegend": False,
                                                   "height": 760}),

    ("compare_participant_violin_*", {"x_label_size": 24,
                                      "y_label_size": 24,
                                      "x_tick_size": 20,
                                      "y_tick_size": 20,
                                      "showlegend": False,
                                      "height": 760}),

    ("compare_violin_*", {"x_label_size": 24,
                          "y_label_size": 24,
                          "x_tick_size": 20,
                          "y_tick_size": 20,
                          "showlegend": False,
                          "height": 760}),

    ("curve_time_on_task_*", {"x_label_size": 26,
                              "y_label_size": 26,
                              "x_tick_size": 22,
                              "y_tick_size": 22,
                              "legend_size": 20,
                              "legend_x": 0.98,
                              "legend_y": 0.98,
                              "line_width": 4,
                              "height": 800}),

    ("missingness_*_over_trial", {"x_label_size": 26,
                                  "y_label_size": 26,
                                  "x_tick_size": 22,
                                  "y_tick_size": 22,
                                  "legend_size": 20,
                                  "legend_x": 0.98,
                                  "legend_y": 0.98,
                                  "line_width": 4,
                                  "height": 780}),

    ("yaw_forward_fraction_by_context*", {"x_label_size": 22,
                                          "y_label_size": 22,
                                          "x_tick_size": 18,
                                          "y_tick_size": 18,
                                          "legend_size": 18,
                                          "height": 760}),

    ("F2_bar_*", {"x_label_size": 24,
                  "y_label_size": 24,
                  "x_tick_size": 18,
                  "y_tick_size": 20,
                  "height": 820}),

    ("F2_roc_*", {"x_label_size": 24,
                  "y_label_size": 24,
                  "x_tick_size": 20,
                  "y_tick_size": 20,
                  "legend_size": 18,
                  "height": 760}),

    ("F1_violin_*", {"x_label_size": 24,
                     "y_label_size": 24,
                     "x_tick_size": 20,
                     "y_tick_size": 20,
                     "showlegend": False,
                     "height": 760}),

    ("F1_scatter_*", {"x_label_size": 22,
                      "y_label_size": 22,
                      "x_tick_size": 18,
                      "y_tick_size": 18,
                      "height": 740}),
]


def _infer_plot_kind(fig: Any) -> str:
    """Infer a broad plot kind from Plotly traces."""
    try:
        trace_types = {str(getattr(tr, "type", "")).lower() for tr in fig.data}
    except Exception:
        return "default"
    if "violin" in trace_types or "box" in trace_types:
        return "violin"
    if "bar" in trace_types:
        return "bar"
    if "scatter" in trace_types:
        try:
            modes = [str(getattr(tr, "mode", "")) for tr in fig.data if str(getattr(tr,
                                                                                    "type", "")).lower() == "scatter"]
            if any("lines" in m for m in modes):
                return "line"
        except Exception:
            return "line"
        return "scatter"
    return "default"


def _plot_text_style_for_name(fig: Any, name: str) -> Dict[str, Any]:
    """Return the merged text/export-size style for one figure name."""
    style: Dict[str, Any] = dict(PLOT_TEXT_SIZE_DEFAULT)
    kind = _infer_plot_kind(fig)
    style.update(PLOT_TEXT_SIZE_BY_KIND.get(kind, {}))
    for pattern, overrides in PLOT_TEXT_SIZE_BY_PATTERN:
        if fnmatch(str(name), pattern):
            style.update(overrides)
    style.update(PLOT_TEXT_SIZE_BY_NAME.get(str(name), {}))
    return style


def _apply_plot_text_style(fig: Any, name: str) -> Any:
    """Apply figure-specific axis label and tick sizes before export."""
    if fig is None:
        return fig
    style = _plot_text_style_for_name(fig, name)

    try:
        fig.update_layout(
            font=dict(size=style.get("font_size", PLOT_TEXT_SIZE_DEFAULT["font_size"])),
            title_font=dict(size=style.get("title_size", PLOT_TEXT_SIZE_DEFAULT["title_size"])),
        )
        if not style.get("show_title", False):
            fig.update_layout(title_text=None)
    except Exception:
        pass

    try:
        fig.update_xaxes(
            title_font_size=style.get("x_label_size", PLOT_TEXT_SIZE_DEFAULT["x_label_size"]),
            tickfont_size=style.get("x_tick_size", PLOT_TEXT_SIZE_DEFAULT["x_tick_size"]),
            automargin=True,
        )
    except Exception:
        pass

    try:
        fig.update_yaxes(
            title_font_size=style.get("y_label_size", PLOT_TEXT_SIZE_DEFAULT["y_label_size"]),
            tickfont_size=style.get("y_tick_size", PLOT_TEXT_SIZE_DEFAULT["y_tick_size"]),
            automargin=True,
        )
    except Exception:
        pass

    try:
        legend_update = {
            "font": dict(size=style.get("legend_size", PLOT_TEXT_SIZE_DEFAULT["legend_size"])),
            "title": dict(
                text=style.get("legend_title_text", None),
                font=dict(size=style.get("legend_title_size", PLOT_TEXT_SIZE_DEFAULT["legend_title_size"])),
            ),
        }
        for style_key, legend_key in (
            ("legend_x", "x"),
            ("legend_y", "y"),
            ("legend_xanchor", "xanchor"),
            ("legend_yanchor", "yanchor"),
            ("legend_orientation", "orientation"),
            ("legend_traceorder", "traceorder"),
            ("legend_bgcolor", "bgcolor"),
            ("legend_bordercolor", "bordercolor"),
            ("legend_borderwidth", "borderwidth"),
        ):
            if style_key in style:
                legend_update[legend_key] = style[style_key]
        fig.update_layout(legend=legend_update)
    except Exception:
        pass

    try:
        if "showlegend" in style:
            showlegend = bool(style["showlegend"])
            fig.update_layout(showlegend=showlegend)
            for tr in fig.data:
                try:
                    tr.showlegend = showlegend
                except Exception:
                    pass
    except Exception:
        pass

    try:
        if style.get("line_width") is not None:
            lw = float(style["line_width"])
            for tr in fig.data:
                try:
                    if str(getattr(tr, "type", "")).lower() != "scatter":
                        continue
                    mode = str(getattr(tr, "mode", ""))
                    if "lines" not in mode:
                        continue
                    line = getattr(tr, "line", None)
                    current_width = getattr(line, "width", None) if line is not None else None
                    hoverinfo = str(getattr(tr, "hoverinfo", "") or "").lower()
                    fill = str(getattr(tr, "fill", "") or "").lower()
                    # Do not turn hidden confidence-band boundary traces into visible lines.
                    if current_width == 0 or hoverinfo == "skip" or fill not in {"", "none"}:
                        continue
                    tr.line.width = lw
                except Exception:
                    pass
    except Exception:
        pass

    try:
        if style.get("violin_line_width") is not None:
            fig.update_traces(line=dict(width=float(style["violin_line_width"])), selector=dict(type="violin"))
            fig.update_traces(line=dict(width=float(style["violin_line_width"])), selector=dict(type="box"))
    except Exception:
        pass

    try:
        if style.get("marker_line_width") is not None:
            fig.update_traces(marker_line_width=float(style["marker_line_width"]), selector=dict(type="bar"))
    except Exception:
        pass

    try:
        if style.get("marker_size") is not None:
            fig.update_traces(marker=dict(size=float(style["marker_size"])), selector=dict(type="scatter"))
    except Exception:
        pass

    try:
        fig.update_annotations(font_size=style.get("annotation_size", PLOT_TEXT_SIZE_DEFAULT["annotation_size"]))
    except Exception:
        pass

    # Mark figures styled by csu_core so helper.save_plotly does not reapply
    # its own per-name defaults and override values edited here.
    try:
        existing_meta = getattr(fig.layout, "meta", None)
        meta = dict(existing_meta) if isinstance(existing_meta, dict) else {}
        meta["csu_core_style_applied"] = True
        meta["csu_plot_name"] = str(name)
        fig.update_layout(meta=meta)
    except Exception:
        pass

    return fig


def _plot_export_size(name: str, fig: Any) -> Tuple[int, int]:
    """Return width and height for static exports of a named figure."""
    style = _plot_text_style_for_name(fig, name)
    return int(style.get("width", 1320)), int(style.get("height", 680))


def _normalise_dataset_token(value: Any) -> Any:
    if value is None:
        return value
    s = str(value)
    key = s.strip().lower().replace("_", " ")
    if key in DATASET_LABEL_MAP:
        return DATASET_LABEL_MAP[key]
    return value


def _rgba_from_hex(hex_colour: str, alpha: float) -> str:
    c = str(hex_colour).lstrip("#")
    if len(c) != 6:
        return f"rgba(0,0,0,{alpha})"
    r = int(c[0:2], 16)
    g = int(c[2:4], 16)
    b = int(c[4:6], 16)
    return f"rgba({r},{g},{b},{alpha})"


def _humanise_label(value: Any) -> Any:
    """Convert raw column names into publication-ready British English labels."""
    if value is None:
        return value
    if not isinstance(value, str):
        return value
    s = value.strip()
    if not s:
        return s

    # dataset names first
    norm = _normalise_dataset_token(s)
    if isinstance(norm, str) and norm != s:
        return norm

    direct = {
        'dataset': 'Dataset',
        'participant_id': 'Participant ID',
        'video_id': 'Video ID',
        'condition_name': 'Condition',
        'trial_index': 'Trial number',
        'trial_pos_norm': 'Normalised trial position',
        'bin_center': 'Normalised trial position',
        'y_mean': 'Mean value',
        'sem': 'Standard error',
        'q1': 'Q1 rating',
        'q2': 'Q2 rating',
        'q3': 'Q3 rating',
        'Q1': 'Q1 rating',
        'Q2': 'Q2 rating',
        'Q3': 'Q3 rating',
        'trigger_mean': 'Mean unsafety',
        'avg_trigger': 'Mean unsafety',
        'mean_trigger': 'Mean unsafety',
        'frac_time_unsafe': 'Fraction of trial time in unsafe zone',
        'unsafe_time_frac': 'Fraction of trial time in unsafe zone',
        'frac_unsafe': 'Fraction of trial time in unsafe zone',
        'dtrigger_sd': 'Unsafety volatility',
        'trigger_sd': 'Unsafety volatility',
        'volatility': 'Unsafety volatility',
        'n_transitions': 'Number of unsafety transitions',
        'transitions': 'Number of unsafety transitions',
        'num_transitions': 'Number of unsafety transitions',
        'latency_first_press_s': 'Latency to first button press (s)',
        'latency_first_release_s': 'Latency to first button release (s)',
        'press_release_hysteresis_s': 'Interval from first press to first release (s)',
        'yielding': 'Yielding',
        'eHMIOn': 'eHMI on',
        'prev_yielding': 'Previous trial yielding',
        'prev_eHMIOn': 'Previous trial eHMI on',
        'prev_camera': 'Previous trial camera condition',
        'prev_distPed': 'Previous trial pedestrian distance',
        'distPed': 'Pedestrian distance',
        'camera': 'Camera condition',
        'yaw_abs_mean': 'Mean absolute head yaw (deg)',
        'yaw_forward_frac_15': 'Fraction of time looking forwards',
        'yaw_sd': 'Standard deviation of head yaw (deg)',
        'yaw_entropy': 'Head yaw entropy',
        'yaw_speed_mean': 'Mean head yaw speed (deg/s)',
        'head_turn_count_15': 'Number of head turns',
        'head_turn_dwell_mean_s_15': 'Mean head turn dwell time (s)',
        'yaw_speed_pre_press_mean_1s': 'Mean head yaw speed before button press (deg/s)',
        'lag_turn_to_press_s_15': 'Lag from head turn to button press (s)',
        'n_trials_main': 'Number of main trials',
        'completion_frac': 'Completion fraction',
    }
    if s in direct:
        return direct[s]

    # compound metric names
    prefixes = [
        ('carryover_prev_yielding_', 'Carryover effect of previous yielding on {}'),
        ('carryover_prev_eHMIOn_', 'Carryover effect of previous eHMI on on {}'),
        ('carryover_prev_camera_', 'Carryover effect of previous camera condition on {}'),
        ('carryover_prev_distPed_', 'Carryover effect of previous pedestrian distance on {}'),
        ('drift_late_minus_early_', 'Late minus early change in {}'),
        ('slope_', 'Linear slope of {} over trial position'),
        ('post_break_reset_', 'Post break reset in {}'),
        ('break_reset_matched_', 'Break matched change in {}'),
        ('switch_cost_', 'Switch cost in {}'),
        ('missing_frac_', 'Missing fraction for {}'),
        ('baseline_mean_', 'Baseline mean of {}'),
        ('curve_time_on_task_', '{} over trial position'),
    ]
    for prefix, template in prefixes:
        if s.startswith(prefix):
            rest = s[len(prefix):]
            return template.format(_humanise_label(rest).lower())

    # general clean-up
    pretty = s.replace('_', ' ')
    pretty = re.sub(r'\bnormalized\b', 'normalised', pretty, flags=re.IGNORECASE)
    pretty = re.sub(r'\bunshuffled\b', 'Fixed order', pretty, flags=re.IGNORECASE)
    pretty = re.sub(r'\bshuffled\b', 'Randomised', pretty, flags=re.IGNORECASE)
    pretty = re.sub(r'\s+', ' ', pretty).strip()
    if pretty.islower():
        pretty = pretty.capitalize()
    return pretty


def _sanitise_axis_like(obj: Any) -> None:
    try:
        title = getattr(obj, 'title', None)
        if title is not None:
            text = getattr(title, 'text', None)
            if isinstance(text, str) and text:
                title.text = _humanise_label(text)
    except Exception:
        pass
    try:
        ticktext = getattr(obj, 'ticktext', None)
        if ticktext:
            obj.ticktext = tuple(_humanise_label(x) for x in ticktext)
    except Exception:
        pass


def _dataset_colour_from_any_label(label: Any) -> Optional[str]:
    if label is None:
        return None

    s = str(label).strip()
    if s in DATASET_COLOUR_MAP:
        return DATASET_COLOUR_MAP[s]

    key = s.lower().replace("_", " ")
    if key in {"shuffled", "randomised", "randomized"}:
        return DATASET_COLOUR_MAP["shuffled"]
    if key in {"unshuffled", "fixed order", "fixed-order"}:
        return DATASET_COLOUR_MAP["unshuffled"]
    return None


def _split_single_dataset_violin(fig: Any) -> None:
    """Turn a single dataset violin/box trace into one trace per dataset so colours can differ."""
    try:
        if len(fig.data) != 1:
            return

        tr = fig.data[0]
        if getattr(tr, "type", None) not in {"violin", "box"}:
            return

        x = list(getattr(tr, "x", []) or [])
        y = list(getattr(tr, "y", []) or [])
        if not x or len(x) != len(y):
            return

        x_norm = [_normalise_dataset_token(v) for v in x]
        wanted = ["Randomised", "Fixed order"]
        if not any(v in wanted for v in x_norm):
            return

        new_traces = []

        for ds in wanted:
            idx = [i for i, v in enumerate(x_norm) if v == ds]
            if not idx:
                continue

            xs = [ds] * len(idx)
            ys = [y[i] for i in idx]
            colour = DATASET_COLOUR_MAP.get(ds)

            common = dict(
                x=xs,
                y=ys,
                name=ds,
                legendgroup=ds,
                box_visible=bool(getattr(tr, "box", None) and getattr(tr.box, "visible", False)) or True,
                meanline_visible=bool(getattr(tr, "meanline", None) and getattr(tr.meanline, "visible", False)),
                points=getattr(tr, "points", "outliers"),
                jitter=getattr(tr, "jitter", None),
                pointpos=getattr(tr, "pointpos", None),
                opacity=getattr(tr, "opacity", None) or 0.75,
                showlegend=True,
            )

            if getattr(tr, "type", None) == "violin":
                nt = go.Violin(**common)
                nt.fillcolor = _rgba_from_hex(colour, 0.35)  # type: ignore[arg-type]
            else:
                nt = go.Box(**common)

            nt.line = dict(color=colour)
            nt.marker = dict(color=colour)
            new_traces.append(nt)

        if new_traces:
            fig.data = tuple(new_traces)

    except Exception:
        return


def _sanitise_figure_for_export(fig: Any) -> Any:
    if fig is None:
        return fig

    try:
        _split_single_dataset_violin(fig)
    except Exception:
        pass

    try:
        for tr in fig.data:
            try:
                if hasattr(tr, "name") and isinstance(tr.name, str):
                    tr.name = _humanise_label(tr.name)

                if hasattr(tr, "legendgroup") and isinstance(tr.legendgroup, str):
                    tr.legendgroup = _humanise_label(tr.legendgroup)

                ds_colour = None
                for cand in [getattr(tr, "name", None), getattr(tr, "legendgroup", None)]:
                    ds_colour = _dataset_colour_from_any_label(cand)
                    if ds_colour is not None:
                        break

                if ds_colour is not None:
                    if hasattr(tr, "line") and tr.line is not None:
                        tr.line.color = ds_colour

                    if hasattr(tr, "marker") and tr.marker is not None:
                        tr.marker.color = ds_colour
                        try:
                            if getattr(tr.marker, "line", None) is not None and getattr(tr.marker.line, "color", None) is None:  # noqa:E501
                                tr.marker.line.color = ds_colour
                        except Exception:
                            pass

                    if hasattr(tr, "fillcolor"):
                        # Lighter fill for violins and filled areas
                        tr.fillcolor = _rgba_from_hex(ds_colour, 0.20)

                if hasattr(tr, "x") and tr.x is not None:
                    tr.x = tuple(_normalise_dataset_token(v) if isinstance(v, str) else v for v in tr.x)

                if hasattr(tr, "y") and tr.y is not None and getattr(tr, "type", None) == "bar":
                    tr.y = tuple(_normalise_dataset_token(v) if isinstance(v, str) else v for v in tr.y)

                if hasattr(tr, "hovertemplate") and isinstance(tr.hovertemplate, str):
                    ht = tr.hovertemplate.replace("_", " ")
                    ht = ht.replace("unshuffled", "Fixed order").replace("shuffled", "Randomised")
                    ht = ht.replace("Fixed-order", "Fixed order")
                    tr.hovertemplate = ht

            except Exception:
                pass

        # Layout level palette so Plotly Express traces also follow the same order
        try:
            fig.update_layout(
                colorway=[
                    DATASET_COLOUR_MAP["shuffled"],
                    DATASET_COLOUR_MAP["unshuffled"],
                ],
                legend_title_text="",
            )
        except Exception:
            pass

        # layout titles
        try:
            if getattr(fig.layout, "title", None) and isinstance(fig.layout.title.text, str):
                fig.layout.title.text = _humanise_label(fig.layout.title.text)
        except Exception:
            pass

        for k in dir(fig.layout):
            if k.startswith("xaxis") or k.startswith("yaxis"):
                try:
                    _sanitise_axis_like(getattr(fig.layout, k))
                except Exception:
                    pass

        try:
            if getattr(fig.layout, "legend", None) and getattr(fig.layout.legend, "title", None):
                txt = getattr(fig.layout.legend.title, "text", None)
                if isinstance(txt, str):
                    fig.layout.legend.title.text = _humanise_label(txt)
        except Exception:
            pass

        try:
            anns = getattr(fig.layout, "annotations", None) or []
            for ann in anns:
                if hasattr(ann, "text") and isinstance(ann.text, str):
                    ann.text = _humanise_label(ann.text)
        except Exception:
            pass

    except Exception:
        pass

    try:
        fig = apply_global_plotly_style(fig)
    except Exception:
        pass

    return fig


# -----------------------
# Plot saving utilities
# -----------------------

_PLOTS_WRITTEN: List[Tuple[str, str]] = []  # (name, html_path) in primary output dir
_KALEIDO_WARNED = False


def _np_trapezoid(y: np.ndarray, x: np.ndarray) -> float:
    # Avoid numpy deprecation warnings (np.trapz -> np.trapezoid in newer numpy).
    if hasattr(np, "trapezoid"):
        return float(np.trapezoid(y, x))  # type: ignore[attr-defined]
    return float(np.trapz(y, x))  # type: ignore


def _safe_corr(x: pd.Series, y: pd.Series, min_n: int = 3) -> float:
    """Pearson correlation with guards for NaNs / constants."""
    mask = (~x.isna()) & (~y.isna())
    if mask.sum() < min_n:
        return np.nan
    xx = x.loc[mask]
    yy = y.loc[mask]
    if xx.nunique(dropna=True) < 2 or yy.nunique(dropna=True) < 2:
        return np.nan
    return float(xx.corr(yy))


def _is_near_constant(arr: NDArray[Any] | None, range_thr: float = 1e-6, sd_thr: float = 1e-6) -> bool:
    """True if array is (almost) constant."""
    if arr is None:
        return True

    a = np.asarray(arr, dtype=np.float64)
    a = a[np.isfinite(a)]
    if a.size < 2:
        return True

    rng: float = float(np.nanmax(a) - np.nanmin(a))
    sd: float = float(np.nanstd(a, ddof=1))

    return bool((rng <= range_thr) or (sd <= sd_thr))


def _perm_pvalue_mean_diff(xa: np.ndarray, xb: np.ndarray, n_resamples: int = 20000, seed: int = 0) -> float:
    """Two-sided permutation p-value for difference in means."""
    xa = np.asarray(xa, dtype=float)
    xb = np.asarray(xb, dtype=float)
    xa = xa[np.isfinite(xa)]
    xb = xb[np.isfinite(xb)]
    if xa.size == 0 or xb.size == 0:
        return np.nan
    obs = float(np.mean(xa) - np.mean(xb))
    pooled = np.concatenate([xa, xb])
    n_a = xa.size
    rng = np.random.default_rng(seed)
    count = 0
    for _ in range(int(n_resamples)):
        perm = rng.permutation(pooled)
        diff = float(np.mean(perm[:n_a]) - np.mean(perm[n_a:]))
        if abs(diff) >= abs(obs) - 1e-15:
            count += 1
    return float((count + 1) / (n_resamples + 1))


def _safe_two_sample_test(xa_in: Any, xb_in: Any, *, equal_var: bool = False, min_n: int = 2,
                          const_range_thr: float = 1e-6, const_sd_thr: float = 1e-6,
                          n_perm: int = 20000, seed: int = 0) -> Tuple[float, float, str, str]:
    """Robust 2-sample comparison returning (t_or_nan, p, test_name, note).

    - If both groups are (near) constant and equal -> p=1 (identical).
    - If (near) constant or precision-loss warning occurs -> permutation p-value.
    - Otherwise Welch t-test.
    Also tries Mann–Whitney U when appropriate (heavy ties).
    """
    xa = np.asarray(pd.to_numeric(pd.Series(xa_in), errors="coerce"), dtype=float)
    xb = np.asarray(pd.to_numeric(pd.Series(xb_in), errors="coerce"), dtype=float)
    xa = xa[np.isfinite(xa)]
    xb = xb[np.isfinite(xb)]
    if xa.size < min_n or xb.size < min_n:
        return (np.nan, np.nan, "insufficient_n", "n<2 in one group")

    mean_a = float(np.mean(xa))
    mean_b = float(np.mean(xb))
    sd_a = float(np.std(xa, ddof=1)) if xa.size > 1 else 0.0
    sd_b = float(np.std(xb, ddof=1)) if xb.size > 1 else 0.0

    a_const = _is_near_constant(xa, const_range_thr, const_sd_thr)
    b_const = _is_near_constant(xb, const_range_thr, const_sd_thr)

    if a_const and b_const:
        if np.isclose(mean_a, mean_b, atol=const_range_thr, rtol=0.0):
            return (0.0, 1.0, "constant_equal", "both groups (near) constant and equal")
        p_perm = _perm_pvalue_mean_diff(xa, xb, n_resamples=max(2000, n_perm // 5), seed=seed)
        return (np.nan, float(p_perm), "perm_mean_diff", "both groups (near) constant; Welch t-test unstable")

    # If one group constant or very small variance, Welch can be unstable; use Mann–Whitney or permutation
    if a_const or b_const or (sd_a < const_sd_thr * 10) or (sd_b < const_sd_thr * 10):
        # Try Mann–Whitney U (handles non-normal + ties); fall back to permutation.
        try:
            u_stat, p_u = mannwhitneyu(xa, xb, alternative="two-sided", method="auto")
            return (np.nan, float(p_u), "mannwhitneyu", "one group low-variance/constant; used MWU")
        except Exception:
            p_perm = _perm_pvalue_mean_diff(xa, xb, n_resamples=n_perm, seed=seed)
            return (np.nan, float(p_perm), "perm_mean_diff", "one group low-variance/constant; MWU failed")

    # Welch t-test, but capture precision-loss warning and fall back to permutation p-value if it happens
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always", RuntimeWarning)
        t_stat, p_val = ttest_ind(xa, xb, equal_var=equal_var)
    warn_text = "; ".join([str(ww.message) for ww in w]) if w else ""
    if "Precision loss occurred" in warn_text:
        p_perm = _perm_pvalue_mean_diff(xa, xb, n_resamples=n_perm, seed=seed)
        return (float(t_stat), float(p_perm), "welch+perm", "precision-loss warning; using permutation p-value")

    return (float(t_stat), float(p_val), "welch", "")


def _zscore(s: pd.Series) -> pd.Series:
    s = pd.to_numeric(s, errors="coerce")
    m = float(s.mean())
    sd = float(s.std(ddof=0))
    if sd == 0 or np.isnan(sd):
        return s * np.nan
    return (s - m) / sd


def _fisher_z(r: float) -> float:
    if r is None or np.isnan(r):
        return np.nan
    if abs(r) >= 1:
        return np.nan
    return float(np.arctanh(r))


def _xcorr_max_r_lag(x: np.ndarray, y: np.ndarray, dt_s: float, max_lag_s: float = 2.0, min_n: int = 10) -> tuple:
    """Max (absolute) Pearson cross-correlation and its lag.

    We compute r(k) = corr(x[t], y[t+k]). Positive lag means y occurs *after* x.
    In our usage: x=yaw_speed, y=trigger_derivative -> positive lag suggests head movement precedes risk updates.
    """
    if x is None or y is None:
        return (np.nan, np.nan)
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    n = min(x.size, y.size)
    if n < min_n or (not np.isfinite(dt_s)) or dt_s <= 0:
        return (np.nan, np.nan)
    x = x[:n]
    y = y[:n]

    max_k = int(round(max_lag_s / dt_s))
    if max_k < 1:
        return (np.nan, np.nan)

    best_r = np.nan
    best_k = 0

    for k in range(-max_k, max_k + 1):
        if k == 0:
            continue
        if k > 0:
            xx = x[:-k]
            yy = y[k:]
        else:
            kk = -k
            xx = x[kk:]
            yy = y[:-kk]

        if xx.size < min_n:
            continue
        mask = np.isfinite(xx) & np.isfinite(yy)
        if mask.sum() < min_n:
            continue

        xxm = xx[mask]
        yym = yy[mask]
        if np.unique(xxm).size < 2 or np.unique(yym).size < 2:
            continue

        r = np.corrcoef(xxm, yym)[0, 1]
        if not np.isfinite(r):
            continue

        if (not np.isfinite(best_r)) or (abs(r) > abs(best_r)):
            best_r = float(r)
            best_k = int(k)

    if not np.isfinite(best_r):
        return (np.nan, np.nan)
    return (best_r, float(best_k) * float(dt_s))


def _trial_num_display(x: pd.Series) -> pd.Series:
    """Convert 0-based trial_index to 1-based trial number for plotting.

    Many internal computations use 0-based indexing (0..39). For presentation,
    we display trials as 1..40.
    """
    xs = pd.to_numeric(x, errors="coerce")
    try:
        mn = float(np.nanmin(xs.to_numpy(dtype=float)))
    except Exception:
        return xs
    if np.isfinite(mn) and mn == 0.0:
        return xs + 1
    return xs


def _read_csv_flexible(fp: str) -> Optional[pd.DataFrame]:
    for sep in [",", ";", "\t"]:
        try:
            return pd.read_csv(fp, sep=sep)
        except Exception:
            continue
    return None


def _probe_csv_columns(fp: str) -> Optional[List[str]]:
    """Read only the header (and optionally 1 row) to get column names.

    This is used to avoid accidentally selecting processed/summary CSVs that share the same
    video_id in their filename but do not contain the raw HMD rotation columns.
    """
    for sep in [",", ";", "\t"]:
        try:
            df = pd.read_csv(fp, sep=sep, nrows=1)
            return [str(c) for c in df.columns]
        except Exception:
            continue
    return None


def _score_trial_file(fp: str) -> int:
    """Heuristic score for selecting the *raw* trial CSV for yaw extraction."""
    cols = _probe_csv_columns(fp)
    if not cols:
        return 0
    colset = set(cols)
    score = 0

    # Highest priority: the exact raw HMD quaternion columns used by the prior pipeline
    req = {"HMDRotationW", "HMDRotationX", "HMDRotationY", "HMDRotationZ"}
    if req.issubset(colset):
        score += 100

    # General quaternion columns (list-like or separate) – moderate priority
    if any(_QUAT_LIST_COL_PAT.search(c) for c in colset):
        score += 20
    # If we can match all 4 components via regex, add score
    comp_hits = 0
    for comp in ["w", "x", "y", "z"]:
        if any(_QUAT_REGEX[comp].search(c) for c in colset):
            comp_hits += 1
    if comp_hits == 4:
        score += 20

    # Yaw columns are helpful but lower priority (can be non-HMD yaw)
    if any(c in colset for c in _YAW_CANDIDATES):
        score += 10

    if "Timestamp" in colset:
        score += 5

    return score


def _choose_trial_file(csv_files: List[str], vid: str, folder: str) -> Optional[str]:
    """Pick the best candidate CSV for a given video_id within a participant folder.

    Replicates the behaviour of the previous pipeline (helper.export_participant_quaternion_matrix):
    prefer an exact file named '{video_id}.csv'. Otherwise, pick the file whose *columns* most
    resemble the raw trial log (contains HMDRotationW/X/Y/Z).
    """
    # 1) Exact match: <folder>/<vid>.csv
    exact = os.path.join(folder, f"{vid}.csv")
    if os.path.isfile(exact):
        return exact

    # 2) Candidates containing the vid in the basename
    candidates = [fp for fp in csv_files if vid in os.path.basename(fp)]
    if not candidates:
        return None

    # 3) If any basename exactly matches, prefer that
    exact2 = [fp for fp in candidates if os.path.basename(fp) == f"{vid}.csv"]
    if exact2:
        return exact2[0]

    # 4) Score by header contents
    scored = sorted(((_score_trial_file(fp), fp) for fp in candidates), reverse=True)
    best_score, best_fp = scored[0]

    # Tie-break: shorter basename tends to be the raw file (fewer prefixes/suffixes)
    top = [fp for sc, fp in scored if sc == best_score]
    if len(top) > 1:
        top = sorted(top, key=lambda x: (len(os.path.basename(x)), os.path.basename(x)))
        best_fp = top[0]

    return best_fp


def _bh_fdr(pvals: pd.Series) -> pd.Series:
    """Benjamini–Hochberg FDR correction (returns q-values aligned to input index)."""
    p = pd.to_numeric(pvals, errors="coerce")
    q = pd.Series(np.nan, index=p.index, dtype=float)

    p_non = p.dropna()
    n = int(p_non.shape[0])
    if n == 0:
        return q

    # Force a real numpy float array so np.argsort / arithmetic are well-typed
    pv: NDArray[np.float64] = p_non.to_numpy(dtype=np.float64, copy=False)

    order = np.argsort(pv)                # NDArray[int]
    ranked = pv[order]                    # NDArray[float]
    qs = ranked * (n / np.arange(1, n + 1, dtype=np.float64))
    qs = np.minimum.accumulate(qs[::-1])[::-1]

    q.loc[p_non.index[order]] = qs
    return q


def compare_participant_metrics(df: pd.DataFrame, metrics: List[str], group_col: str = "dataset", a: str = "shuffled",
                                b: str = "unshuffled", fdr: bool = True) -> pd.DataFrame:
    """Welch t-test + Cohen's d on participant-level metrics."""
    rows = []
    for m in metrics:
        if m not in df.columns:
            continue
        x = pd.to_numeric(df[df[group_col] == a][m], errors="coerce").dropna()
        y = pd.to_numeric(df[df[group_col] == b][m], errors="coerce").dropna()
        if len(x) < 2 or len(y) < 2:
            continue

        t, p, test_name, test_note = _safe_two_sample_test(x, y, equal_var=False)

        nx, ny = len(x), len(y)
        sx, sy = float(x.std(ddof=1)), float(y.std(ddof=1))
        pooled = np.sqrt(((nx - 1) * sx * sx + (ny - 1) * sy * sy) / (nx + ny - 2)) if (nx + ny - 2) > 0 else np.nan
        d = (float(x.mean()) - float(y.mean())) / pooled if pooled and not np.isnan(pooled) else np.nan

        rows.append({
            "metric": m,
            "n_shuffled": nx,
            "n_unshuffled": ny,
            "mean_shuffled": float(x.mean()),
            "mean_unshuffled": float(y.mean()),
            "t_stat": float(t),
            "p_value": float(p) if p is not None else np.nan,
            "test": test_name,
            "note": test_note,
            "cohens_d": float(d) if d is not None else np.nan,
        })

    out = pd.DataFrame(rows)
    if out.empty:
        return out
    out["q_value"] = _bh_fdr(out["p_value"]) if fdr else np.nan
    out = out.sort_values(["q_value", "p_value"]) if fdr else out.sort_values("p_value")
    return out


def _cohens_d(xa: np.ndarray, xb: np.ndarray) -> float:
    """Cohen's d using pooled SD (ignores NaNs; expects 1D arrays)."""
    xa = np.asarray(xa, dtype=float)
    xb = np.asarray(xb, dtype=float)
    xa = xa[~np.isnan(xa)]
    xb = xb[~np.isnan(xb)]
    if len(xa) < 2 or len(xb) < 2:
        return np.nan
    nx, ny = len(xa), len(xb)
    sx, sy = np.std(xa, ddof=1), np.std(xb, ddof=1)
    denom = (nx + ny - 2)
    if denom <= 0:
        return np.nan
    pooled = np.sqrt(((nx - 1) * sx * sx + (ny - 1) * sy * sy) / denom)
    if pooled == 0 or np.isnan(pooled):
        return np.nan
    return float((np.mean(xa) - np.mean(xb)) / pooled)


def _build_cond_key(df: pd.DataFrame) -> pd.Series:
    """Create a compact condition key from available mapping factors.

    Used when a feature table lacks a precomputed `condition_name` but comparisons
    request grouping by it.
    """
    cols: List[str] = []
    for c in ["yielding", "eHMIOn", "camera", "distPed"]:
        if c in df.columns:
            cols.append(c)

    if cols:
        parts = []
        for c in cols:
            parts.append(c + "=" + df[c].astype(str))
        return parts[0] if len(parts) == 1 else parts[0].str.cat(parts[1:], sep="|")

    if "condition_name" in df.columns:
        return df["condition_name"].astype(str)
    if "video_id" in df.columns:
        return df["video_id"].astype(str)
    return pd.Series(["cond"] * len(df), index=df.index)


def compare_shuffled_unshuffled(
    features_df: pd.DataFrame,
    dataset_col: str = "dataset",
    shuffled_label: str = "shuffled",
    unshuffled_label: str = "unshuffled",
    metrics: Optional[List[str]] = None,
    groupby: Optional[List[str]] = None,
    equal_var: bool = False,
    fdr: bool = True,
    *,
    unit: str = "trial",
    participant_col: str = "participant_id",
) -> pd.DataFrame:
    """Statistical comparison between shuffled vs unshuffled.

    IMPORTANT
    - If unit is "trial" (default), the function treats rows as independent observations.
      This is fine for between subject tables where each participant contributes at most
      one row per group (eg grouped by condition_name).
    - If unit is "participant", the function aggregates each metric to a participant level
      mean within each group (dataset + groupby + participant_id) before running the test.
      This avoids the clustering violation that occurs when many trials per participant
      are analysed as if they were independent.

    Output
    Produces a tidy table with group keys (if any), metric name, n mean sd per dataset,
    t test p values, Cohen's d, and optional BH FDR corrected p values.

    This is implemented locally so the script does not depend on helper.HMD_helper having
    the method (repo versions may differ).
    """
    if groupby is None:
        groupby = []
    df = features_df.copy()
    if dataset_col not in df.columns:
        return pd.DataFrame()

    # Robustness: some feature tables (e.g., yaw-only exports) may not carry the
    # mapping-derived condition_name. If requested, synthesize a stable key from
    # available factors before grouping.
    missing_gb = [c for c in groupby if c not in df.columns]
    if missing_gb:
        if (missing_gb == ["condition_name"]) and ("condition_name" not in df.columns):
            df["condition_name"] = _build_cond_key(df)
            missing_gb = []
        else:
            # Drop unavailable groupers rather than raising.
            keep = [c for c in groupby if c in df.columns]
            if keep != groupby:
                logger.warning(f"[compare] Dropping missing group-by columns: {missing_gb}. Using: {keep if keep else 'overall'}")  # noqa: E501
            groupby = keep

    df = df[df[dataset_col].isin([shuffled_label, unshuffled_label])]
    if df.empty:
        return pd.DataFrame()

    # Default: compare a sensible set of metrics (trigger + yaw), falling back to numeric columns.
    if metrics is None:
        default_metrics = [
            # trigger
            "frac_time_unsafe",
            "trigger_mean",
            "trigger_peak",
            "trigger_p95",
            "max_ramp_rate",
            "trigger_sd",
            "dtrigger_sd",
            "latency_first_press_s",
            "latency_first_release_s",
            "n_safe_to_unsafe",
            "n_unsafe_to_safe",
            "mean_unsafe_bout_s",
            "longest_unsafe_bout_s",
            "press_release_hysteresis",
            "time_to_brake_first_press_s",
            "anticipation_index",
            # yaw
            "yaw_forward_frac_15",
            "yaw_mean",
            "yaw_sd",
            "yaw_abs_mean",
            "yaw_range_deg",
            "yaw_turn_rate_mean_dps",
            "yaw_turn_rate_sd_dps",
        ]
        metrics = [m for m in default_metrics if m in df.columns]
        if not metrics:
            # Fall back to numeric columns excluding obvious IDs/factors.
            exclude = set([dataset_col, "participant_id", "video_id", "condition_name", "source_file"])
            metrics = [c for c in df.columns if c not in exclude and pd.api.types.is_numeric_dtype(df[c])]

    # Force numeric where possible
    for m in metrics:
        if m in df.columns:
            df[m] = pd.to_numeric(df[m], errors="coerce")

    # Optional: participant level aggregation to avoid clustering violations when many
    # trials per participant are present.
    unit_norm = str(unit).strip().lower()
    if unit_norm in {"participant", "participant_mean", "cluster", "clustered"}:
        if participant_col in df.columns:
            agg_cols: List[str] = []
            if groupby:
                agg_cols.extend(list(groupby))
            agg_cols.extend([dataset_col, participant_col])
            # Each metric becomes the participant mean within each group
            df = (
                df.groupby(agg_cols, dropna=False)[metrics]
                .mean()
                .reset_index()
            )
        else:
            logger.warning(f"[compare] unit=participant requested but {participant_col} missing; using trial level")
            unit_norm = "trial"
    else:
        unit_norm = "trial"

    if groupby:
        groups = df.groupby(groupby, dropna=False)
    else:
        groups = [((), df)]

    rows = []
    for gkey, gdf in groups:
        gdict = {}
        if groupby:
            if not isinstance(gkey, tuple):
                gkey = (gkey,)
            gdict = dict(zip(groupby, gkey))

        a = gdf[gdf[dataset_col] == shuffled_label]
        b = gdf[gdf[dataset_col] == unshuffled_label]

        for m in metrics:
            if m not in gdf.columns:
                continue
            xa = a[m].to_numpy(dtype=float)
            xb = b[m].to_numpy(dtype=float)
            xa = xa[~np.isnan(xa)]
            xb = xb[~np.isnan(xb)]

            if len(xa) < 2 or len(xb) < 2:
                t_stat, p_val, d = np.nan, np.nan, np.nan
                test_name, test_note = "insufficient_n", "n<2 in one group"
            else:
                t_stat, p_val, test_name, test_note = _safe_two_sample_test(xa, xb, equal_var=equal_var)
                d = _cohens_d(xa, xb)

            rows.append({
                **gdict,
                "unit": unit_norm,
                "metric": m,
                "n_shuffled": int(len(xa)),
                "mean_shuffled": float(np.mean(xa)) if len(xa) else np.nan,
                "sd_shuffled": float(np.std(xa, ddof=1)) if len(xa) > 1 else np.nan,
                "n_unshuffled": int(len(xb)),
                "mean_unshuffled": float(np.mean(xb)) if len(xb) else np.nan,
                "sd_unshuffled": float(np.std(xb, ddof=1)) if len(xb) > 1 else np.nan,
                "t": float(t_stat) if not np.isnan(t_stat) else np.nan,
                "p": float(p_val) if not np.isnan(p_val) else np.nan,
                "test": test_name,
                "note": test_note,
                "cohens_d": float(d) if not np.isnan(d) else np.nan,
            })

    out = pd.DataFrame(rows)
    if out.empty:
        return out

    if fdr:
        out["p_fdr"] = np.nan
        if groupby:
            for _, idxs in out.groupby(groupby, dropna=False).groups.items():
                out.loc[list(idxs), "p_fdr"] = _bh_fdr(out.loc[list(idxs), "p"]).values
        else:
            out["p_fdr"] = _bh_fdr(out["p"]).values

    return out


def _print_table(df: pd.DataFrame, title: str, max_rows: int = 12) -> None:
    logger.info("\n" + title)
    if df is None or df.empty:
        logger.info("(empty)")
        return
    with pd.option_context("display.max_rows", max_rows, "display.max_columns", 50, "display.width", 160):
        logger.info(_to_string_3dp(df.head(max_rows), index=False))


def _pick_col(df: pd.DataFrame, candidates: List[str]) -> Optional[str]:
    for c in candidates:
        if c in df.columns:
            return c
    return None


def _get_output_dir_for_logs() -> str:
    try:
        if common is not None:
            return str(common.get_configs("output"))
    except Exception:
        pass
    return "(see common.get_configs('output'))"


def _resolve_plot_dirs(h: HMD_helper, out_root: Optional[str] = None) -> List[str]:
    out_dirs: List[str] = []
    # 1) common output (usually "_output")
    try:
        if common is not None:
            out_dirs.append(str(common.get_configs("output")))
    except Exception:
        pass
    # 2) explicit out_root (usually "_compare_output")
    if out_root:
        out_dirs.append(str(out_root))
    # 3) figures folder (often where the user is looking)
    try:
        if getattr(h, "folder_figures", None):
            out_dirs.append(str(h.folder_figures))
    except Exception:
        pass

    # de-dup while keeping order
    seen = set()
    return [d for d in out_dirs if d and not (d in seen or seen.add(d))]


def _plotly_static_export_setup() -> None:
    """Make Plotly static export more robust before writing PNG/EPS/PDF."""
    try:
        import plotly.io as pio  # local import so the module still imports without Plotly IO extras
        try:
            pio.kaleido.scope.mathjax = None
        except Exception:
            pass
    except Exception:
        pass


def _convert_pdf_to_eps(pdf_path: str, eps_path: str) -> tuple[bool, str]:
    """Convert a Plotly PDF export to EPS using common command line tools."""
    import shutil
    import subprocess

    pdftops = shutil.which("pdftops")
    if pdftops:
        try:
            res = subprocess.run(
                [pdftops, "-eps", pdf_path, eps_path],
                stdout=subprocess.PIPE,
                stderr=subprocess.PIPE,
                text=True,
                check=False,
            )
            if res.returncode == 0 and os.path.exists(eps_path) and os.path.getsize(eps_path) > 0:
                return True, "pdftops"
            return False, f"pdftops failed: {res.stderr.strip() or res.stdout.strip()}"
        except Exception as e:
            return False, f"pdftops exception: {e}"

    return False, "pdftops not found"


def _write_plotly_static_exports(fig: Any, base_path: str, width: int, height: int) -> Dict[str, str]:
    """Write PNG and EPS for one Plotly figure.

    EPS is attempted with direct Plotly export first. If that fails, a PDF is
    exported and converted to EPS with pdftops. SVG export is intentionally not
    used, because this project does not need SVG files.

    Returns a mapping from extension to status text. The caller logs these statuses.
    """
    statuses: Dict[str, str] = {}
    _plotly_static_export_setup()

    png_path = f"{base_path}.png"
    try:
        fig.write_image(png_path, width=width, height=height)
        statuses["png"] = "ok"
    except Exception as e:
        statuses["png"] = f"failed: {e}"

    eps_path = f"{base_path}.eps"
    try:
        fig.write_image(eps_path, width=width, height=height)
        if os.path.exists(eps_path) and os.path.getsize(eps_path) > 0:
            statuses["eps"] = "ok"
            return statuses
        statuses["eps"] = "failed: empty EPS after direct export"
    except Exception as e:
        statuses["eps"] = f"direct failed: {e}"

    pdf_path = f"{base_path}.pdf"
    try:
        fig.write_image(pdf_path, width=width, height=height)
        ok, msg = _convert_pdf_to_eps(pdf_path, eps_path)
        if ok:
            statuses["eps"] = f"ok via {msg}"
            return statuses
        statuses["eps"] = statuses.get("eps", "") + f"; pdf conversion failed: {msg}"
    except Exception as e:
        statuses["eps"] = statuses.get("eps", "") + f"; pdf export failed: {e}"

    return statuses


def _save_plot(h: HMD_helper, fig, name: str, out_root: Optional[str] = None, record_index: bool = True,
               open_browser: bool = True, save_final: bool = True) -> None:
    """
    Save a Plotly figure.

    Preferred behaviour (when available): use `HMD_helper.save_plotly(...)` so figures can pop up
    one by one (interactive Plotly view).

    Fallback behaviour: write HTML via `fig.write_html(...)` and export PNG and EPS (best effort;
    requires Kaleido).

    Notes
    - Set environment variable `CSU_OPEN_BROWSER=1` to force opening each plot in your browser.
    - The function still writes plots to the relevant output folders.
    - Set ``save_final=True`` to mirror HTML, PNG and EPS into ``h.folder_figures``
      when the helper supports it and in the fallback path.
    """
    global _KALEIDO_WARNED
    if fig is None:
        return
    fig = _sanitise_figure_for_export(fig)
    fig = _apply_plot_text_style(fig, name)
    export_width, export_height = _plot_export_size(name, fig)

    # env var override (handy when you do not want to touch call sites)
    env_open = os.environ.get('CSU_OPEN_BROWSER', None)
    if env_open is not None and str(env_open).strip() != '':
        open_browser = str(env_open).strip().lower() in {'1', 'true', 'yes', 'y'}

    # If we have the project's helper, prefer it. It handles Kaleido MathJax quirks and browser opening.
    # NOTE: Some variants of HMD_helper.save_plotly also save PDF outputs. This project does not need PDF.
    # We therefore:
    # 1) Pass flags to disable PDF where supported.
    # 2) Propagate save_final so HTML, PNG and EPS can also be written to the figures folder.
    # 3) Always do an explicit EPS export ourselves as a backstop.
    try:
        if hasattr(h, "save_plotly"):
            # HMD_helper.save_plotly saves to `common.get_configs("output")` (and optionally to h.folder_figures).
            # Only use it when that matches the requested out_root, otherwise we fall back to multi dir saving.
            out_dirs = _resolve_plot_dirs(h, out_root=out_root)
            common_out = out_dirs[0] if out_dirs else None

            def _same_dir(a: Optional[str], b: Optional[str]) -> bool:
                if not a or not b:
                    return False
                try:
                    return os.path.abspath(str(a)) == os.path.abspath(str(b))
                except Exception:
                    return False

            if (out_root is None) or _same_dir(out_root, common_out):
                # Build kwargs defensively across helper versions.
                # Only pass parameters that exist in the helper signature.
                import inspect

                params = set()
                try:
                    params = set(inspect.signature(h.save_plotly).parameters.keys())
                except Exception:
                    params = set()

                kwargs = {}
                if 'name' in params:
                    kwargs['name'] = name
                if 'remove_margins' in params:
                    kwargs['remove_margins'] = False
                if 'width' in params:
                    kwargs['width'] = export_width
                if 'height' in params:
                    kwargs['height'] = export_height
                if 'save_eps' in params:
                    kwargs['save_eps'] = True
                if 'save_png' in params:
                    kwargs['save_png'] = True
                if 'save_html' in params:
                    kwargs['save_html'] = True
                if 'open_browser' in params:
                    kwargs['open_browser'] = bool(open_browser)
                if 'save_mp4' in params:
                    kwargs['save_mp4'] = False
                # Disable PDF like outputs where the helper supports it.
                if 'save_pdf' in params:
                    kwargs['save_pdf'] = False
                # Propagate save_final so helper based saves can mirror HTML, PNG and EPS
                # into the figures folder. Any PDF side effects are cleaned up below.
                if 'save_final' in params:
                    kwargs['save_final'] = bool(save_final)

                h.save_plotly(fig, **kwargs)

                # Even if the helper ran, ensure EPS exists (and remove any PDF created).
                out_dirs = _resolve_plot_dirs(h, out_root=out_root)
                for d in out_dirs:
                    try:
                        os.makedirs(d, exist_ok=True)
                    except Exception:
                        continue

                    # Static export backstop. This is intentionally independent of
                    # HMD_helper.save_plotly because some helper versions fail during
                    # PDF-to-EPS conversion even when HTML was written correctly.
                    try:
                        statuses = _write_plotly_static_exports(fig, os.path.join(d, name), export_width, export_height)  # noqa:E501
                        if statuses.get("eps", "").startswith("ok"):
                            logger.info(f"[plot] EPS wrote: {os.path.join(d, f'{name}.eps')} ({statuses.get('eps')})")
                        elif not _KALEIDO_WARNED:
                            _KALEIDO_WARNED = True
                            logger.warning(f"[plot] NOTE: EPS export failed for {name}: {statuses.get('eps')}")
                    except Exception as e:
                        if not _KALEIDO_WARNED:
                            _KALEIDO_WARNED = True
                            logger.warning(f"[plot] NOTE: static export failed for {name}: {e}")

                    # Remove PDF if the helper created one
                    try:
                        import glob
                        patterns = [
                            os.path.join(d, f"{name}*.pdf"),
                            os.path.join(d, "final", f"{name}*.pdf"),
                            os.path.join(d, "pdf", f"{name}*.pdf"),
                        ]
                        if save_final and getattr(h, "folder_figures", None):
                            patterns.extend([
                                os.path.join(str(h.folder_figures), f"{name}*.pdf"),
                                os.path.join(str(h.folder_figures), "final", f"{name}*.pdf"),
                                os.path.join(str(h.folder_figures), "pdf", f"{name}*.pdf"),
                            ])
                        for pat in patterns:
                            for pdf_path in glob.glob(pat):
                                try:
                                    if os.path.isfile(pdf_path):
                                        os.remove(pdf_path)
                                except Exception:
                                    pass
                    except Exception:
                        pass

                return
    except Exception as e:
        logger.error(f"[plot] NOTE: HMD_helper.save_plotly failed for {name}: {e}. Falling back to fig.write_html.")

    # ---- Fallback multi directory saving ----
    out_dirs = _resolve_plot_dirs(h, out_root=out_root)
    if not out_dirs:
        logger.warning(f"[plot] SKIP {name}: no output directories resolved")
        return

    wrote_html: List[str] = []
    primary_out = out_dirs[0]

    for d in out_dirs:
        try:
            os.makedirs(d, exist_ok=True)
        except Exception:
            continue

        # HTML
        html_path = os.path.join(d, f"{name}.html")
        try:
            fig.write_html(html_path, include_plotlyjs="cdn", full_html=True, auto_open=False)
            wrote_html.append(html_path)
            if record_index and (d == primary_out):
                _PLOTS_WRITTEN.append((name, html_path))
        except Exception as e:
            logger.error(f"[plot] FAILED html {html_path}: {e}")

        # PNG and EPS. EPS is tried directly and then through PDF/SVG conversion.
        try:
            statuses = _write_plotly_static_exports(fig, os.path.join(d, name), export_width, export_height)
            if statuses.get("eps", "").startswith("ok"):
                logger.info(f"[plot] EPS wrote: {os.path.join(d, f'{name}.eps')} ({statuses.get('eps')})")
            elif not _KALEIDO_WARNED:
                _KALEIDO_WARNED = True
                logger.warning(f"[plot] NOTE: EPS export failed for {name}: {statuses.get('eps')}")
        except Exception as e:
            if not _KALEIDO_WARNED:
                _KALEIDO_WARNED = True
                logger.warning(f"[plot] NOTE: static export failed for {name}: {e}")

    # Open the primary HTML for this figure if requested
    if wrote_html and open_browser:
        try:
            import webbrowser
            webbrowser.open("file://" + os.path.abspath(wrote_html[0]))
        except Exception:
            pass

    if wrote_html:
        logger.info(f"[plot] wrote: {name} -> " + ", ".join(wrote_html))
    else:
        logger.info(f"[plot] FAILED to write any HTML for {name}")


def _write_plot_index_and_open(h: HMD_helper) -> None:
    """Write a single HTML index linking all recorded plots, then auto-open it (one tab)."""
    if not _PLOTS_WRITTEN:
        return

    out_dirs = _resolve_plot_dirs(h, out_root=None)
    if not out_dirs:
        return
    out_dir = out_dirs[0]
    try:
        os.makedirs(out_dir, exist_ok=True)
    except Exception:
        return

    # build simple index
    rows = []
    for name, path in _PLOTS_WRITTEN:
        fn = os.path.basename(path)
        rows.append(f'<li><a href="{fn}">{name}</a></li>')
    html = f"""<!doctype html>
<html><head><meta charset="utf-8"><title>Multiped plots</title></head>
<body>
<h2>Multiped plots index</h2>
<p>Saved in: {out_dir}</p>
<ul>
{''.join(rows)}
</ul>
</body></html>"""

    index_path = os.path.join(out_dir, "index_plots.html")
    try:
        with open(index_path, "w", encoding="utf-8") as f:
            f.write(html)
        # auto-open (one tab)
        try:
            import webbrowser
            webbrowser.open("file://" + os.path.abspath(index_path))
        except Exception:
            pass
    except Exception as e:
        logger.error(f"[plot] FAILED to write index: {e}")


def _lin_slope(x: np.ndarray, y: np.ndarray) -> dict:
    """Return slope/intercept/r/p/n using scipy.stats.linregress; robust to NaNs."""
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    mask = np.isfinite(x) & np.isfinite(y)
    x = x[mask]
    y = y[mask]
    if x.size < 5 or np.unique(x).size < 2 or np.unique(y).size < 2:
        return {"n": int(x.size), "slope": np.nan, "r": np.nan, "p": np.nan}
    res = linregress(x, y)
    return {"n": int(x.size), "slope": float(res.slope), "r": float(res.rvalue), "p": float(res.pvalue)}


def _fisher_ci(r: float, n: int) -> Tuple[float, float]:
    """95% CI for Pearson r via Fisher z transform."""
    if n is None or n < 4 or r is None or np.isnan(r) or abs(r) >= 1.0:
        return (np.nan, np.nan)
    z = np.arctanh(r)
    se = 1.0 / math.sqrt(n - 3)
    lo = np.tanh(z - 1.96 * se)
    hi = np.tanh(z + 1.96 * se)
    return (float(lo), float(hi))


def _p_from_z(z: float) -> float:
    """Two-sided p-value from z using erfc (no extra deps)."""
    if z is None or np.isnan(z):
        return np.nan
    return float(math.erfc(abs(z) / math.sqrt(2.0)))


def _compare_independent_corr(r1: float, n1: int, r2: float, n2: int) -> Tuple[float, float]:
    """Fisher r-to-z test for difference between independent correlations."""
    if any(v is None or np.isnan(v) for v in [r1, r2]) or n1 < 4 or n2 < 4:
        return (np.nan, np.nan)
    z1 = np.arctanh(r1)
    z2 = np.arctanh(r2)
    se = math.sqrt(1.0 / (n1 - 3) + 1.0 / (n2 - 3))
    if se <= 0:
        return (np.nan, np.nan)
    z = (z1 - z2) / se
    return (float(z), _p_from_z(z))


def _read_csv_loose(path: str) -> Optional[pd.DataFrame]:
    if path is None or (not isinstance(path, str)) or (not os.path.exists(path)):
        return None
    try:
        return pd.read_csv(path)
    except Exception:
        try:
            return pd.read_csv(path, sep=";")
        except Exception as e:
            logger.error(f"[Q] Could not read CSV {path}: {e}")
            return None


def _roc_curve_and_auc(y: np.ndarray, score: np.ndarray):
    # Simple ROC curve + AUC (no sklearn dependency)
    y = np.asarray(y).astype(int)
    score = np.asarray(score).astype(float)
    mask = (~np.isnan(score)) & (~np.isnan(y))
    y = y[mask]
    score = score[mask]
    # Need both classes
    if y.size < 3 or len(np.unique(y)) < 2:
        return None, None, np.nan

    order = np.argsort(-score)  # descending score
    y_sorted = y[order]
    score_sorted = score[order]  # noqa: F841

    tp = np.cumsum(y_sorted == 1)
    fp = np.cumsum(y_sorted == 0)
    tp_total = tp[-1]
    fp_total = fp[-1]
    if tp_total == 0 or fp_total == 0:
        return None, None, np.nan

    tpr = tp / tp_total
    fpr = fp / fp_total

    # Add (0,0) at start
    tpr = np.concatenate([[0.0], tpr])
    fpr = np.concatenate([[0.0], fpr])

    auc = _np_trapezoid(tpr, fpr)
    return fpr, tpr, auc


def _corr_with_p(x: pd.Series, y: pd.Series, min_n: int = 5) -> Tuple[float, float, int]:
    mask = (~x.isna()) & (~y.isna())
    n0 = int(mask.sum())
    if n0 < min_n:
        return (np.nan, np.nan, n0)

    xx = pd.to_numeric(x.loc[mask], errors="coerce")
    yy = pd.to_numeric(y.loc[mask], errors="coerce")

    mask2 = (~xx.isna()) & (~yy.isna())
    xx = xx.loc[mask2]
    yy = yy.loc[mask2]

    n = int(xx.shape[0])
    if n < min_n or xx.nunique(dropna=True) < 2 or yy.nunique(dropna=True) < 2:
        return (np.nan, np.nan, n)

    # Force real numpy float arrays for SciPy typing + runtime consistency
    x_arr = xx.to_numpy(dtype=np.float64, copy=False)
    y_arr = yy.to_numpy(dtype=np.float64, copy=False)

    r, p = pearsonr(x_arr, y_arr)
    return (float(r), float(p), n)


def _ensure_dir(path: str) -> None:
    if path and not os.path.exists(path):
        os.makedirs(path, exist_ok=True)


# ---------------------------------------------------------------------------
# Star-import compatibility
# ---------------------------------------------------------------------------
# The refactor keeps many helpers with leading underscores to minimize changes
# to the original script bodies. Other modules import this file with
# `from csu_core import *`, but Python's default star-import behavior excludes
# underscore-prefixed names unless they are listed in __all__.
#
# To preserve the original code behavior without editing thousands of call
# sites, we export all public and underscore-prefixed symbols here.
__all__ = [k for k in globals().keys() if not k.startswith("__")]  # pyright: ignore[reportUnsupportedDunderAll]