multichannel_interface.py

"""TNFR field observations on graphs reconstructed from multichannel signals.

Unlike the delay-coordinate graph in :mod:`tnfr.validation.temporal_interface`,
this adapter supplies one node per simultaneous signal channel. Hilbert phases,
amplitude-envelope pressure and a phase-locking graph are observation choices;
they do not establish physical wiring or a closed nodal evolution law.

The shared field owners compute local wrapped phase separation ``|∇φ|``,
circular curvature ``K_φ``, pressure aggregation ``Φ_s`` and ``ξ_C``. The last
is a static pressure-coherence product fit or a separate spectral fallback.
``MultichannelWindowSeries`` currently retains only its numeric value, so this
adapter cannot identify its fit/fallback provenance or certify a measured
correlation length.

Pipeline
--------
``multi-channel signals -> per-channel Hilbert phase/amplitude -> phase-locking
coupling graph (nodes = channels) -> TNFR spatial tetrad per window`` compared
against the Kuramoto order parameter ``R``, mean phase-locking value and
spatial phase dispersion.

Scope and honesty
-----------------
- Local phase separation and global phase order can correlate; neither
  determines the other on an arbitrary supplied graph.
- Amplitude and phase are different transforms of the same signals. An
  amplitude-only pressure proxy does not prove their statistical independence.
- Different field formulas and finite univariate AUC rankings do not establish
  additional predictive information, physical emergence or a measured law.
- All functions are read-only telemetry; the only graph mutation is the
  construction of new coupling graphs and the setting of ``phase``/``theta`` and
  ``dnfr``/``delta_nfr`` node attributes at build time.

References
----------
- ``src/tnfr/validation/temporal_interface.py`` — single-series counterpart
- ``src/tnfr/physics/canonical.py`` — canonical tetrad field functions
- AGENTS.md §"Telemetry & Structural Field Tetrad"
- Kuramoto, "Chemical Oscillations, Waves, and Turbulence" (1984);
  Acebrón et al., Rev. Mod. Phys. 77, 137 (2005) — order parameter.
"""

from __future__ import annotations

from dataclasses import dataclass, field
from typing import Any, Mapping, Sequence

try:
    import numpy as np
except ImportError:  # pragma: no cover - numpy is a core dependency
    np = None  # type: ignore[assignment]

try:
    import networkx as nx
except ImportError:  # pragma: no cover - optional dependency guard
    nx = None  # type: ignore[assignment]

from ..physics.canonical import (
    CoherenceLengthEstimate,
    compute_phase_curvature,
    compute_phase_gradient,
    compute_structural_potential,
    estimate_coherence_length_with_provenance,
)
from .temporal_interface import _coherence_comparison, _coherence_window_metadata

__all__ = [
    "MultichannelConfig",
    "MultichannelWindowSeries",
    "SynchronyDiscrimination",
    "fft_bandpass",
    "analytic_phase_amplitude",
    "phase_amplitude_matrices",
    "mean_field_order",
    "kuramoto_order_parameter",
    "phase_locking_matrix",
    "phase_offsets",
    "amplitude_pressure",
    "build_coupling_graph",
    "multichannel_window_series",
    "evaluate_synchrony_discrimination",
]

_EPS = 1e-12


def _require_numpy() -> None:
    if np is None:  # pragma: no cover - core dependency guard
        raise RuntimeError("numpy is required for multi-channel analysis")


def _require_networkx() -> None:
    if nx is None:  # pragma: no cover - optional dependency guard
        raise RuntimeError("networkx is required for multi-channel analysis")


# ---------------------------------------------------------------------------
# Configuration and result containers
# ---------------------------------------------------------------------------


@dataclass(frozen=True)
class MultichannelConfig:
    """Configuration for the multichannel signal-observation pipeline.

    Parameters
    ----------
    window:
        Number of samples per analysis window.
    step:
        Stride (in samples) between successive windows.
    k_neighbours:
        Neighbours per channel in the phase-locking coupling graph.
    sampling_rate:
        Sampling rate (Hz) of the signals; only used when ``bandpass`` is set.
    bandpass:
        Optional ``(low, high)`` band (Hz) applied to each channel before the
        Hilbert transform.  ``None`` uses the broadband signal.
    """

    window: int = 512
    step: int = 128
    k_neighbours: int = 4
    sampling_rate: float = 1.0
    bandpass: tuple[float, float] | None = None

    def __post_init__(self) -> None:
        if self.window < 16:
            raise ValueError("window must be >= 16 samples")
        if self.step < 1:
            raise ValueError("step must be >= 1")
        if self.k_neighbours < 1:
            raise ValueError("k_neighbours must be >= 1")
        if self.sampling_rate <= 0.0:
            raise ValueError("sampling_rate must be > 0")
        if self.bandpass is not None:
            low, high = self.bandpass
            if not (0.0 <= low < high):
                raise ValueError("bandpass must satisfy 0 <= low < high")


@dataclass(frozen=True)
class MultichannelWindowSeries:
    """Per-window TNFR spatial tetrad plus the matched synchronisation baselines.

    All arrays are aligned by index; entry ``i`` corresponds to the window whose
    most-recent sample is ``window_end[i]``.  TNFR channels are ``grad_phi``,
    ``k_phi``, ``xi_c`` and ``phi_s``; baselines are ``order_parameter``
    (Kuramoto ``R``), ``mean_plv`` and ``phase_dispersion``.
    Raw xi values retain aligned canonical estimates with method and units;
    numeric-only legacy construction cannot establish a comparable xi series.
    """

    window_end: "np.ndarray"
    grad_phi: "np.ndarray"
    k_phi: "np.ndarray"
    xi_c: "np.ndarray"
    phi_s: "np.ndarray"
    order_parameter: "np.ndarray"
    mean_plv: "np.ndarray"
    phase_dispersion: "np.ndarray"
    xi_c_provenance: tuple[CoherenceLengthEstimate, ...] = ()

    def as_dict(self) -> dict[str, Any]:
        return {
            "window_end": [int(v) for v in self.window_end],
            "grad_phi": [float(v) for v in self.grad_phi],
            "k_phi": [float(v) for v in self.k_phi],
            "xi_c": [float(v) for v in self.xi_c],
            "phi_s": [float(v) for v in self.phi_s],
            "order_parameter": [float(v) for v in self.order_parameter],
            "mean_plv": [float(v) for v in self.mean_plv],
            "phase_dispersion": [float(v) for v in self.phase_dispersion],
            "xi_c_provenance": [
                _coherence_window_metadata(item) for item in self.xi_c_provenance
            ],
        }


@dataclass(frozen=True)
class SynchronyDiscrimination:
    """Honest comparison of TNFR tetrad vs synchronisation baselines.

    ``auc`` maps each indicator to its discriminative ROC-AUC against a binary
    regime label (``max(auc, 1 - auc)`` so the arbitrary regime direction does
    not matter).  An AUC near 0.5 means the indicator does not separate the two
    regimes; near 1.0 means clean separation. If the labelled windows lack
    either class, ``metadata['auc_available']`` is false and the retained 0.5
    values are compatibility placeholders, not performance measurements.
    Per-channel availability separately excludes an incomparable xi series
    from best-channel selection while preserving the other comparisons.
    """

    indicators: tuple[str, ...]
    auc: Mapping[str, float]
    tnfr_indicators: tuple[str, ...]
    baseline_indicators: tuple[str, ...]
    best_tnfr: tuple[str, float]
    best_baseline: tuple[str, float]
    n_windows: int
    n_positive_windows: int
    interpretation: str
    metadata: Mapping[str, Any] = field(default_factory=dict)

    def summary(self) -> str:
        if self.metadata.get("auc_available") is False:
            return self.interpretation
        bt, bv = self.best_tnfr
        bb, bbv = self.best_baseline
        return (
            f"best TNFR: {bt}={bv:.3f} | best baseline: {bb}={bbv:.3f} | "
            f"windows={self.n_windows} (positive={self.n_positive_windows})"
        )


# ---------------------------------------------------------------------------
# Signal -> phase/amplitude
# ---------------------------------------------------------------------------


def fft_bandpass(
    signal: Sequence[float], *, sampling_rate: float, low: float, high: float
) -> "np.ndarray":
    """Zero-phase FFT band-pass (pure numpy; no SciPy dependency).

    A brick-wall mask keeps Fourier components with frequency in ``[low, high]``.
    Suitable for isolating an oscillatory band (e.g. EEG alpha) before phase
    extraction; ringing from the sharp cut-off is acceptable for telemetry.
    """
    _require_numpy()
    x = np.asarray(signal, dtype=float)
    n = x.size
    if n < 2:
        return x.copy()
    x = x - float(np.mean(x))
    freqs = np.fft.rfftfreq(n, d=1.0 / float(sampling_rate))
    spectrum = np.fft.rfft(x)
    mask = (freqs >= low) & (freqs <= high)
    spectrum = spectrum * mask
    return np.fft.irfft(spectrum, n=n)


def _analytic_signal(x: "np.ndarray") -> "np.ndarray":
    """Discrete analytic signal via FFT (identical to ``scipy.signal.hilbert``)."""
    n = x.size
    spectrum = np.fft.fft(x)
    h = np.zeros(n, dtype=float)
    if n % 2 == 0:
        h[0] = 1.0
        h[n // 2] = 1.0
        h[1 : n // 2] = 2.0
    else:
        h[0] = 1.0
        h[1 : (n + 1) // 2] = 2.0
    return np.fft.ifft(spectrum * h)


def analytic_phase_amplitude(
    signal: Sequence[float],
    *,
    sampling_rate: float = 1.0,
    bandpass: tuple[float, float] | None = None,
) -> tuple["np.ndarray", "np.ndarray"]:
    """Return measured instantaneous ``(phase, amplitude)`` of one channel.

    The signal mean is removed first; if ``bandpass`` is given the signal is
    band-limited before the Hilbert transform so the analytic phase tracks a
    well-defined oscillatory component.
    """
    _require_numpy()
    x = np.asarray(signal, dtype=float)
    if x.ndim != 1:
        raise ValueError("signal must be one-dimensional")
    n = x.size
    if n < 2:
        return np.zeros(n, dtype=float), np.zeros(n, dtype=float)
    x = x - float(np.mean(x))
    if bandpass is not None:
        x = fft_bandpass(
            x, sampling_rate=sampling_rate, low=bandpass[0], high=bandpass[1]
        )
    analytic = _analytic_signal(x)
    return np.angle(analytic), np.abs(analytic)


def phase_amplitude_matrices(
    signals: "np.ndarray",
    *,
    sampling_rate: float = 1.0,
    bandpass: tuple[float, float] | None = None,
) -> tuple["np.ndarray", "np.ndarray"]:
    """Compute phase and amplitude matrices for all channels.

    ``signals`` has shape ``(n_channels, n_samples)``.  Returns ``(phase,
    amplitude)`` of the same shape.
    """
    _require_numpy()
    data = np.asarray(signals, dtype=float)
    if data.ndim != 2:
        raise ValueError("signals must have shape (n_channels, n_samples)")
    n_channels, n_samples = data.shape
    phase = np.empty_like(data)
    amplitude = np.empty_like(data)
    for j in range(n_channels):
        phase[j], amplitude[j] = analytic_phase_amplitude(
            data[j], sampling_rate=sampling_rate, bandpass=bandpass
        )
    return phase, amplitude


# ---------------------------------------------------------------------------
# Synchronisation descriptors
# ---------------------------------------------------------------------------


def mean_field_order(
    phase_matrix: "np.ndarray",
) -> tuple["np.ndarray", "np.ndarray"]:
    """Return the Kuramoto order time series ``(R(t), Ψ(t))``.

    ``R(t) = |⟨e^{iφ_j(t)}⟩_j|`` is the instantaneous global synchrony and
    ``Ψ(t)`` the collective phase.
    """
    _require_numpy()
    phases = np.asarray(phase_matrix, dtype=float)
    z = np.mean(np.exp(1j * phases), axis=0)
    return np.abs(z), np.angle(z)


def kuramoto_order_parameter(phase_matrix: "np.ndarray") -> float:
    """Mean Kuramoto order parameter ``⟨R(t)⟩`` over the window (baseline)."""
    r, _ = mean_field_order(phase_matrix)
    return float(np.mean(r))


def phase_locking_matrix(phase_matrix: "np.ndarray") -> "np.ndarray":
    """Pairwise phase-locking value matrix ``PLV_{jk} = |⟨e^{i(φ_j-φ_k)}⟩_t|``.

    Returns an ``(n_channels, n_channels)`` matrix with unit diagonal.
    """
    _require_numpy()
    phases = np.asarray(phase_matrix, dtype=float)
    n_samples = phases.shape[1]
    e = np.exp(1j * phases)
    plv = np.abs(e @ e.conj().T) / float(max(1, n_samples))
    np.fill_diagonal(plv, 1.0)
    return np.clip(plv, 0.0, 1.0)


def phase_offsets(phase_matrix: "np.ndarray") -> "np.ndarray":
    """Mean phase offset of each channel from the collective rhythm.

    ``θ_j = angle(⟨e^{i(φ_j(t) - Ψ(t))}⟩_t)`` where ``Ψ(t)`` is the mean-field
    phase.  Channels locked to the collective have ``θ_j ≈ 0``; leading/lagging
    channels have ``θ_j ≠ 0``.  The result is a bounded scalar per channel,
    suitable as the TNFR node phase.
    """
    _require_numpy()
    phases = np.asarray(phase_matrix, dtype=float)
    _, psi = mean_field_order(phases)
    relative = phases - psi[np.newaxis, :]
    return np.angle(np.mean(np.exp(1j * relative), axis=1))


def amplitude_pressure(amplitude_matrix: "np.ndarray") -> "np.ndarray":
    """Supply an amplitude-envelope pressure proxy without a phase argument.

    ``ΔNFR_j = |⟨A_j⟩_t - ⟨A⟩| / (⟨A⟩ + ε)`` — the relative deviation of each
    channel's mean envelope from the network mean. It does not consume phase
    explicitly; that does not establish statistical independence or an
    independently measured physical pressure.
    """
    _require_numpy()
    amp = np.asarray(amplitude_matrix, dtype=float)
    mean_amp = np.mean(amp, axis=1)
    global_mean = float(np.mean(mean_amp))
    return np.abs(mean_amp - global_mean) / (global_mean + _EPS)


# ---------------------------------------------------------------------------
# Coupling graph
# ---------------------------------------------------------------------------


def build_coupling_graph(
    phase_matrix: "np.ndarray",
    amplitude_matrix: "np.ndarray",
    *,
    k_neighbours: int,
) -> Any:
    """Build the phase-locking coupling graph for one window.

    Nodes are channels.  Edges connect each channel to its ``k_neighbours`` most
    strongly phase-locked partners (coupling distance ``1 − PLV``).  Each node
    carries the measured phase offset (``phase``/``theta``) and the amplitude
    pressure (``dnfr``/``delta_nfr``) so the canonical tetrad field functions
    operate on real telemetry.
    """
    _require_numpy()
    _require_networkx()
    phases = np.asarray(phase_matrix, dtype=float)
    amps = np.asarray(amplitude_matrix, dtype=float)
    n_channels = phases.shape[0]

    plv = phase_locking_matrix(phases)
    distance = 1.0 - plv
    np.fill_diagonal(distance, 0.0)
    theta = phase_offsets(phases)
    pressure = amplitude_pressure(amps)

    G = nx.Graph()
    for j in range(n_channels):
        G.add_node(
            j,
            phase=float(theta[j]),
            theta=float(theta[j]),
            dnfr=float(pressure[j]),
            delta_nfr=float(pressure[j]),
        )

    k = max(1, min(int(k_neighbours), n_channels - 1))
    for j in range(n_channels):
        order = np.argsort(distance[j])
        added = 0
        for idx in order:
            neighbour = int(idx)
            if neighbour == j:
                continue
            G.add_edge(j, neighbour, distance=float(distance[j, neighbour]))
            added += 1
            if added >= k:
                break
    return G


def _mean_abs_dict(values: Mapping[Any, float]) -> float:
    if not values:
        return float("nan")
    return float(np.mean([abs(v) for v in values.values()]))


# ---------------------------------------------------------------------------
# Windowed tetrad series
# ---------------------------------------------------------------------------


def _as_channel_matrix(signals: "np.ndarray") -> "np.ndarray":
    data = np.asarray(signals, dtype=float)
    if data.ndim != 2:
        raise ValueError("signals must have shape (n_channels, n_samples)")
    return data


def multichannel_window_series(
    signals: "np.ndarray", *, config: MultichannelConfig | None = None
) -> MultichannelWindowSeries:
    """Compute the per-window spatial tetrad and synchronisation baselines.

    Parameters
    ----------
    signals:
        Array of shape ``(n_channels, n_samples)``.
    config:
        Pipeline configuration (window, step, neighbours, optional band-pass).
    """
    _require_numpy()
    cfg = config or MultichannelConfig()
    data = _as_channel_matrix(signals)
    n_channels, n_samples = data.shape
    if n_channels < 3:
        raise ValueError("multi-channel analysis needs >= 3 channels")
    if n_samples < cfg.window:
        raise ValueError(
            f"signal length ({n_samples}) shorter than window ({cfg.window})"
        )

    ends: list[int] = []
    grad: list[float] = []
    kphi: list[float] = []
    xi: list[float] = []
    xi_provenance: list[CoherenceLengthEstimate] = []
    phis: list[float] = []
    order: list[float] = []
    plv_mean: list[float] = []
    dispersion: list[float] = []

    upper = np.triu_indices(n_channels, k=1)
    start = 0
    while start + cfg.window <= n_samples:
        seg = data[:, start : start + cfg.window]
        phase, amplitude = phase_amplitude_matrices(
            seg, sampling_rate=cfg.sampling_rate, bandpass=cfg.bandpass
        )
        graph = build_coupling_graph(phase, amplitude, k_neighbours=cfg.k_neighbours)
        grad.append(_mean_abs_dict(compute_phase_gradient(graph)))
        kphi.append(_mean_abs_dict(compute_phase_curvature(graph)))
        estimate = estimate_coherence_length_with_provenance(graph)
        xi.append(float(estimate.value))
        xi_provenance.append(estimate)
        phis.append(_mean_abs_dict(compute_structural_potential(graph)))

        order.append(kuramoto_order_parameter(phase))
        plv = phase_locking_matrix(phase)
        plv_mean.append(float(np.mean(plv[upper])))
        theta = phase_offsets(phase)
        dispersion.append(float(1.0 - np.abs(np.mean(np.exp(1j * theta)))))

        ends.append(start + cfg.window - 1)
        start += cfg.step

    return MultichannelWindowSeries(
        window_end=np.asarray(ends, dtype=int),
        grad_phi=np.asarray(grad, dtype=float),
        k_phi=np.asarray(kphi, dtype=float),
        xi_c=np.asarray(xi, dtype=float),
        phi_s=np.asarray(phis, dtype=float),
        order_parameter=np.asarray(order, dtype=float),
        mean_plv=np.asarray(plv_mean, dtype=float),
        phase_dispersion=np.asarray(dispersion, dtype=float),
        xi_c_provenance=tuple(xi_provenance),
    )


# ---------------------------------------------------------------------------
# Regime discrimination (AUC)
# ---------------------------------------------------------------------------


def _average_ranks(sorted_values: "np.ndarray") -> "np.ndarray":
    n = sorted_values.size
    ranks = np.empty(n, dtype=float)
    i = 0
    while i < n:
        j = i
        while j + 1 < n and sorted_values[j + 1] == sorted_values[i]:
            j += 1
        ranks[i : j + 1] = (i + j) / 2.0 + 1.0
        i = j + 1
    return ranks


def _rank_auc(values: "np.ndarray", labels: "np.ndarray") -> float:
    """Mann–Whitney ROC-AUC (tie-aware); ``labels`` is a boolean array."""
    v = np.asarray(values, dtype=float)
    y = np.asarray(labels, dtype=bool)
    mask = np.isfinite(v)
    v = v[mask]
    y = y[mask]
    n_pos = int(np.count_nonzero(y))
    n_neg = int(y.size - n_pos)
    if n_pos == 0 or n_neg == 0:
        return 0.5
    order = np.argsort(v, kind="mergesort")
    ranks = _average_ranks(v[order])
    rank_of = np.empty_like(ranks)
    rank_of[order] = ranks
    auc = (rank_of[y].sum() - n_pos * (n_pos + 1) / 2.0) / (n_pos * n_neg)
    return float(auc)


def _discriminative_auc(values: "np.ndarray", labels: "np.ndarray") -> float:
    auc = _rank_auc(values, labels)
    return max(auc, 1.0 - auc)


_TNFR_CHANNELS = ("grad_phi", "k_phi", "xi_c", "phi_s")
_BASELINE_CHANNELS = ("order_parameter", "mean_plv", "phase_dispersion")


def _binary_sample_labels(
    labels: Sequence[float], *, n_samples: int | None = None
) -> "np.ndarray":
    """Admit a finite binary vector without coercing labels or losing alignment."""
    arr = np.asarray(labels)
    if arr.ndim != 1 or arr.dtype.kind not in "biuf":
        raise ValueError("labels must be a one-dimensional numeric binary vector")
    if not np.all((arr == 0) | (arr == 1)):
        raise ValueError("labels must contain only finite binary values 0 or 1")
    if n_samples is not None and arr.size != n_samples:
        raise ValueError("labels must match the number of signal samples")
    return arr


def _window_labels(
    labels: Sequence[float], series: MultichannelWindowSeries, window: int
) -> "np.ndarray":
    """Majority binary label per window (label at each window's sample span)."""
    arr = _binary_sample_labels(labels)
    ends = np.asarray(series.window_end)
    if (
        ends.ndim != 1
        or ends.dtype.kind not in "iu"
        or np.any(ends < window - 1)
        or np.any(ends >= arr.size)
    ):
        raise ValueError("window endpoints must index complete labelled sample spans")
    out = np.empty(series.window_end.size, dtype=bool)
    for i, end in enumerate(ends):
        start = int(end) - window + 1
        segment = arr[start : int(end) + 1]
        out[i] = bool(np.mean(segment) >= 0.5)
    return out


def evaluate_synchrony_discrimination(
    signals: "np.ndarray",
    labels: Sequence[float],
    *,
    config: MultichannelConfig | None = None,
) -> SynchronyDiscrimination:
    """Compare TNFR tetrad vs synchrony baselines at separating two regimes.

    ``labels`` is a per-sample binary regime label (0/1) aligned with the
    signal samples.  Each window receives its majority label; for every
    indicator the discriminative ROC-AUC against the window labels is reported.

    Scores describe this supplied labelled sample. Distinct formulas do not
    prove independent information or out-of-sample prediction. The report's
    comparison bands are descriptive policies, not uncertainty intervals.
    Raw ``xi_c`` values retain per-window estimator provenance. Mixed estimators
    or distance conventions are unavailable for a pooled xi comparison.
    """
    _require_numpy()
    cfg = config or MultichannelConfig()
    data = _as_channel_matrix(signals)
    labels = _binary_sample_labels(labels, n_samples=data.shape[1])
    series = multichannel_window_series(data, config=cfg)
    return _evaluate_synchrony_series(
        series, labels, config=cfg, n_channels=data.shape[0], n_samples=data.shape[1]
    )


def _evaluate_synchrony_series(
    series: MultichannelWindowSeries,
    labels: Sequence[float],
    *,
    config: MultichannelConfig,
    n_channels: int,
    n_samples: int,
) -> SynchronyDiscrimination:
    """Score an already computed series; benchmark summaries reuse the same data."""
    cfg = config
    labels = _binary_sample_labels(labels, n_samples=n_samples)
    window_labels = _window_labels(labels, series, cfg.window)
    xi_values, xi_comparison = _coherence_comparison(
        series.xi_c, series.xi_c_provenance
    )

    channels = {
        "grad_phi": series.grad_phi,
        "k_phi": series.k_phi,
        "xi_c": xi_values,
        "phi_s": series.phi_s,
        "order_parameter": series.order_parameter,
        "mean_plv": series.mean_plv,
        "phase_dispersion": series.phase_dispersion,
    }
    auc = {}
    channel_availability = {}
    channel_reasons = {}
    for name, values in channels.items():
        finite_labels = window_labels[np.isfinite(values)]
        available = bool(np.any(finite_labels) and np.any(~finite_labels))
        channel_availability[name] = available
        channel_reasons[name] = (
            xi_comparison["reason"]
            if name == "xi_c" and not xi_comparison["available"]
            else (None if available else "finite_windows_lack_both_label_classes")
        )
        auc[name] = _discriminative_auc(values, window_labels) if available else 0.5

    def _best(group: tuple[str, ...]) -> tuple[str, float]:
        candidates = tuple(name for name in group if channel_availability[name])
        best_name = max(candidates or group, key=lambda n: auc[n])
        return best_name, auc[best_name]

    best_tnfr = _best(_TNFR_CHANNELS)
    best_baseline = _best(_BASELINE_CHANNELS)
    n_windows = int(series.window_end.size)
    n_pos = int(np.count_nonzero(window_labels))
    unavailable_reason = None
    if not n_windows:
        unavailable_reason = "no_windows"
    elif not n_pos:
        unavailable_reason = "missing_positive_class"
    elif n_pos == n_windows:
        unavailable_reason = "missing_negative_class"
    elif not any(channel_availability[name] for name in _TNFR_CHANNELS) or not any(
        channel_availability[name] for name in _BASELINE_CHANNELS
    ):
        unavailable_reason = "no_available_channel_in_required_family"

    gap = best_tnfr[1] - best_baseline[1]
    if unavailable_reason is not None:
        interpretation = (
            "AUC comparison unavailable: finite comparable channels in both "
            "families must cover both label classes. Retained 0.5 scores and best-channel entries are "
            "compatibility placeholders, not evidence of discrimination or "
            "chance performance."
        )
    elif best_tnfr[1] < 0.6 and best_baseline[1] < 0.6:
        interpretation = (
            "Neither the TNFR tetrad nor the synchrony baselines separate the "
            "two regimes well (all AUC < 0.6); the regimes may not differ in "
            "phase-coupling structure."
        )
    elif abs(gap) < 0.03:
        interpretation = (
            f"TNFR '{best_tnfr[0]}' (AUC={best_tnfr[1]:.3f}) and baseline "
            f"'{best_baseline[0]}' (AUC={best_baseline[1]:.3f}) are comparable "
            "regime discriminators."
        )
    elif gap > 0:
        interpretation = (
            f"TNFR '{best_tnfr[0]}' (AUC={best_tnfr[1]:.3f}) separates the "
            f"regimes more strongly than the best baseline '{best_baseline[0]}' "
            f"(AUC={best_baseline[1]:.3f}); Δ={gap:+.3f}."
        )
    else:
        interpretation = (
            f"Baseline '{best_baseline[0]}' (AUC={best_baseline[1]:.3f}) "
            f"separates the regimes more strongly than the best TNFR channel "
            f"'{best_tnfr[0]}' (AUC={best_tnfr[1]:.3f}); Δ={gap:+.3f}."
        )

    if unavailable_reason is None and not channel_availability["xi_c"]:
        interpretation += (
            " xi_c is excluded from channel selection: "
            f"{channel_reasons['xi_c']}; its 0.5 score is an unavailable placeholder."
        )

    return SynchronyDiscrimination(
        indicators=tuple(channels),
        auc=auc,
        tnfr_indicators=_TNFR_CHANNELS,
        baseline_indicators=_BASELINE_CHANNELS,
        best_tnfr=best_tnfr,
        best_baseline=best_baseline,
        n_windows=n_windows,
        n_positive_windows=n_pos,
        interpretation=interpretation,
        metadata={
            "auc_available": unavailable_reason is None,
            "auc_unavailable_reason": unavailable_reason,
            "channel_auc_available": channel_availability,
            "channel_auc_unavailable_reason": channel_reasons,
            "xi_c_comparison": xi_comparison,
            "n_channels": int(n_channels),
            "config": {
                "window": cfg.window,
                "step": cfg.step,
                "k_neighbours": cfg.k_neighbours,
                "sampling_rate": cfg.sampling_rate,
                "bandpass": list(cfg.bandpass) if cfg.bandpass else None,
            },
        },
    )