"""Pure NumPy Fig2 primitives; no solver, phase selection, or probe DTW.

Adapted formula and inclusive ROI bounds from DynamisLab
src/drl_pinball/eval/wake_l2.py::mean_field_metric, roi_slices, inspected
2026-11-02 at HEAD 609116ec26f19766f5e9ce2b0df427d4e24faae8
(working-tree source, not asserted clean). The explicit common-fluid/finite
mask and motion summaries are additions. legacy_test/metrics.py::cycle_similarity
was inspected only: its six probe channels are not a field recovery metric.
"""
import math
import numpy as np


def roi_slices(shape, *, D, xs, y0=None, downstream_solid_extent):
    """Lab-fixed ROI: sensor-relative x/D=-6..14, y/D=-5..5, inclusive.

    Equivalent to pinball-front-relative x/D=4..24 when xs=front+10D.
    Supply actual geometry, e.g. D=40, xs=1600, y0=510.5.
    """
    ny, nx = shape
    y0 = (ny - 1) / 2 if y0 is None else float(y0)
    if D <= 0 or not all(np.isfinite(v) for v in (D, xs, y0, downstream_solid_extent)):
        raise ValueError("invalid geometry")
    x0, x1 = math.ceil(xs - 6 * D), math.floor(xs + 14 * D)
    lo, hi = math.ceil(y0 - 5 * D), math.floor(y0 + 5 * D)
    if x0 < 0 or lo < 0 or x1 >= nx or hi >= ny:
        raise ValueError("ROI outside grid")
    if x0 <= downstream_solid_extent:
        raise ValueError("ROI not downstream of solid extent")
    return slice(lo, hi + 1), slice(x0, x1 + 1)


def common_fluid_roi_mask(fluid_masks, fields, roi):
    """Intersect explicit boolean fluid masks, ROI, and finite 2-D fields.

    Pass every compared role/component/sample to freeze a common support over
    time. Zero velocity never implies solid. Masks must share the lab grid.
    """
    fields = [np.asarray(a, dtype=float) for a in fields]
    masks = [np.asarray(a) for a in fluid_masks]
    if not fields or not masks:
        raise ValueError("fields and explicit fluid masks required")
    shape = fields[0].shape
    if len(shape) != 2 or any(a.shape != shape for a in fields + masks):
        raise ValueError("equal-shape 2-D arrays required")
    if any(a.dtype != np.dtype(bool) for a in masks):
        raise ValueError("fluid masks must be boolean")
    support = np.zeros(shape, dtype=bool)
    ys, xs = roi
    for s, n in zip((ys, xs), shape):
        if not isinstance(s, slice) or s.step not in (None, 1) or not 0 <= (s.start or 0) < (n if s.stop is None else s.stop) <= n:
            raise ValueError("invalid ROI slices")
    support[ys, xs] = True
    for a in masks:
        support &= a
    for a in fields:
        support &= np.isfinite(a)
    if not support.any():
        raise ValueError("empty common finite fluid ROI")
    return support


def field_error(rux, ruy, tux, tuy, u0, support):
    """Eq1: sqrt(mean_support((dux**2 + duy**2)/U0**2)).

    For E(mean field), pass arithmetic mean fields from matched windows.
    Mean instantaneous E requires separate calls per snapshot and then a mean;
    these operations are not interchangeable. No target-energy denominator.
    """
    a = [np.asarray(v, dtype=float) for v in (rux, ruy, tux, tuy)]
    support = np.asarray(support)
    if a[0].ndim != 2 or any(v.shape != a[0].shape for v in a) or support.shape != a[0].shape or support.dtype != np.dtype(bool):
        raise ValueError("equal-shape 2-D fields and boolean support required")
    if not np.isfinite(u0) or u0 <= 0 or not support.any():
        raise ValueError("positive finite U0 and nonempty support required")
    if any(not np.isfinite(v[support]).all() for v in a):
        raise ValueError("nonfinite field on fixed support")
    du = (a[0][support] - a[2][support]) / u0
    dv = (a[1][support] - a[3][support]) / u0
    return float(np.sqrt(np.mean(du * du + dv * dv)))


def motion_summary(displacement):
    """Arithmetic mean, centered RMS, half peak-to-peak along time axis 0.

    Input is displacement in caller-declared units; no implicit D scaling.
    Use the same retained window as field evaluation. No sinusoid assumption.
    """
    a = np.asarray(displacement, dtype=float)
    if a.ndim == 0 or a.shape[0] == 0 or not np.isfinite(a).all():
        raise ValueError("nonempty finite time samples required")
    mean = a.mean(axis=0)
    return {"mean": mean, "centered_rms": np.sqrt(np.mean((a - mean) ** 2, axis=0)),
            "half_range": (a.max(axis=0) - a.min(axis=0)) / 2}