exoskeleton/code/analysis/metrics.py

"""Independent H1--H4 trial metrics.

These functions consume stored arrays only.  They deliberately do not import
the simulation or controller implementations, so paper endpoints can be
reconstructed independently from raw evidence.
"""

from __future__ import annotations

from typing import Any, Mapping

import numpy as np


METRIC_SCHEMA_VERSION = "1.3.0"


class MetricError(ValueError):
    pass


def _vector(value: Any, name: str, *, dtype=float) -> np.ndarray:
    array = np.asarray(value, dtype=dtype)
    if array.ndim != 1 or array.size == 0:
        raise MetricError(f"{name} must be a non-empty one-dimensional array")
    return array


def _boolean_vector(value: Any, name: str) -> np.ndarray:
    """Parse an evidence flag without allowing NaN/nonzero coercion to True."""
    raw = np.asarray(value)
    if raw.ndim != 1 or raw.size == 0:
        raise MetricError(f"{name} must be a non-empty one-dimensional array")
    if np.issubdtype(raw.dtype, np.bool_):
        return raw.astype(bool, copy=False)
    try:
        numeric = np.asarray(value, dtype=float)
    except (TypeError, ValueError) as error:
        raise MetricError(f"{name} must contain only 0/1 flags") from error
    if (
        not np.all(np.isfinite(numeric))
        or not np.all((numeric == 0.0) | (numeric == 1.0))
    ):
        raise MetricError(f"{name} must contain only finite 0/1 flags")
    return numeric.astype(bool)


def _matrix(value: Any, name: str, columns: int | None = None) -> np.ndarray:
    array = np.asarray(value, dtype=float)
    if array.ndim != 2 or array.shape[0] == 0:
        raise MetricError(f"{name} must be a non-empty two-dimensional array")
    if columns is not None and array.shape[1] != columns:
        raise MetricError(f"{name} must have {columns} columns")
    return array


def _same_rows(named: Mapping[str, np.ndarray]) -> int:
    counts = {name: value.shape[0] for name, value in named.items()}
    if len(set(counts.values())) != 1:
        raise MetricError(f"sample counts differ: {counts}")
    return next(iter(counts.values()))


def _finite(array: np.ndarray, name: str) -> None:
    if not np.all(np.isfinite(array)):
        raise MetricError(f"{name} contains non-finite values")


def _wrapped_delta(array: np.ndarray) -> np.ndarray:
    return (array + np.pi) % (2.0 * np.pi) - np.pi


def compute_h1_composite(
    *,
    mapping_valid: Any,
    position_error_m: Any,
    orientation_error_rad: Any,
    q_slave: Any,
    swivel_angle_rad: Any,
    master_step_norm: Any,
    position_threshold_m: float,
    orientation_threshold_rad: float,
    joint_step_threshold_rad: float,
    swivel_step_threshold_rad: float,
    input_step_threshold_rad: float,
    accepted: Any | None = None,
    commanded_reset: Any | None = None,
    degeneracy_transition: Any | None = None,
) -> dict[str, Any]:
    """Compute the trajectory-level H1 ``F_r``, ``D_r``, and ``C_r``."""
    valid = _boolean_vector(mapping_valid, "mapping_valid")
    e_position = _vector(position_error_m, "position_error_m")
    e_orientation = _vector(orientation_error_rad, "orientation_error_rad")
    slave = _matrix(q_slave, "q_slave")
    swivel = _vector(swivel_angle_rad, "swivel_angle_rad")
    input_step = _vector(master_step_norm, "master_step_norm")
    n = _same_rows(
        {
            "mapping_valid": valid,
            "position_error_m": e_position,
            "orientation_error_rad": e_orientation,
            "q_slave": slave,
            "swivel_angle_rad": swivel,
            "master_step_norm": input_step,
        }
    )
    accepted_array = (
        np.ones(n, dtype=bool)
        if accepted is None
        else _boolean_vector(accepted, "accepted")
    )
    reset = (
        np.zeros(n, dtype=bool)
        if commanded_reset is None
        else _boolean_vector(commanded_reset, "commanded_reset")
    )
    degeneracy = (
        np.zeros(n, dtype=bool)
        if degeneracy_transition is None
        else _boolean_vector(
            degeneracy_transition,
            "degeneracy_transition",
        )
    )
    _same_rows(
        {
            "accepted": accepted_array,
            "commanded_reset": reset,
            "degeneracy_transition": degeneracy,
            "mapping_valid": valid,
        }
    )
    for array, name in (
        (e_position, "position_error_m"),
        (e_orientation, "orientation_error_rad"),
        (slave, "q_slave"),
        (swivel, "swivel_angle_rad"),
        (input_step, "master_step_norm"),
    ):
        _finite(array, name)
    thresholds = {
        "position_threshold_m": position_threshold_m,
        "orientation_threshold_rad": orientation_threshold_rad,
        "joint_step_threshold_rad": joint_step_threshold_rad,
        "swivel_step_threshold_rad": swivel_step_threshold_rad,
        "input_step_threshold_rad": input_step_threshold_rad,
    }
    if any(not np.isfinite(value) or value < 0.0 for value in thresholds.values()):
        raise MetricError("H1 thresholds must be finite and non-negative")

    failure_mask = accepted_array & (
        ~valid
        | (e_position > position_threshold_m)
        | (e_orientation > orientation_threshold_rad)
    )
    if n > 1:
        joint_step = np.max(
            np.abs(_wrapped_delta(slave[1:] - slave[:-1])),
            axis=1,
        )
        swivel_step = np.abs(_wrapped_delta(swivel[1:] - swivel[:-1]))
        eligible = (
            accepted_array[1:]
            & accepted_array[:-1]
            & valid[1:]
            & valid[:-1]
            & (input_step[1:] <= input_step_threshold_rad)
            & ~reset[1:]
            & ~degeneracy[1:]
        )
        discontinuity_mask = eligible & (
            (joint_step > joint_step_threshold_rad)
            | (swivel_step > swivel_step_threshold_rad)
        )
    else:
        joint_step = np.empty(0, dtype=float)
        swivel_step = np.empty(0, dtype=float)
        eligible = np.empty(0, dtype=bool)
        discontinuity_mask = np.empty(0, dtype=bool)

    accepted_count = int(np.sum(accepted_array))
    eligible_count = int(np.sum(eligible))
    failure_metric_valid = accepted_count > 0
    discontinuity_metric_valid = eligible_count > 0
    composite_metric_valid = (
        failure_metric_valid and discontinuity_metric_valid
    )
    F_r = int(np.any(failure_mask)) if failure_metric_valid else None
    D_r = (
        int(np.any(discontinuity_mask))
        if discontinuity_metric_valid
        else None
    )
    return {
        "h1_F_r": F_r,
        "h1_D_r": D_r,
        "h1_C_r": (
            max(F_r, D_r) if composite_metric_valid else None
        ),
        "h1_failure_metric_valid": failure_metric_valid,
        "h1_discontinuity_metric_valid": discontinuity_metric_valid,
        "h1_composite_metric_valid": composite_metric_valid,
        "h1_accepted_sample_count": accepted_count,
        "h1_failure_sample_count": int(np.sum(failure_mask)),
        "h1_discontinuity_sample_count": int(np.sum(discontinuity_mask)),
        "h1_eligible_increment_count": eligible_count,
        "h1_mapping_valid_fraction": float(np.mean(valid[accepted_array]))
        if np.any(accepted_array)
        else 0.0,
        "h1_max_position_error_m": float(np.max(e_position[accepted_array]))
        if np.any(accepted_array)
        else 0.0,
        "h1_max_orientation_error_rad": float(np.max(e_orientation[accepted_array]))
        if np.any(accepted_array)
        else 0.0,
        "h1_max_joint_step_rad": float(np.max(joint_step))
        if joint_step.size
        else 0.0,
        "h1_max_swivel_step_rad": float(np.max(swivel_step))
        if swivel_step.size
        else 0.0,
    }


def compute_h1_timing_branch_metrics(
    *,
    pose_runtime_s: Any,
    warm_start: Any,
    phi_rad: Any,
    q_slave: Any,
    branch_smooth: Any,
    differential_applicable: Any,
    differential_valid: Any,
    differential_runtime_s: Any,
    differential_max_one_sided_consistency: Any,
    deadline_s: float = 0.020,
    minimum_phi_wrap_crossings: int = 0,
) -> dict[str, Any]:
    """Summarize H1 pose latency, actual differential, and phi-wrap stress."""
    pose_runtime = _vector(pose_runtime_s, "pose_runtime_s")
    warm = _boolean_vector(warm_start, "warm_start")
    phi = _vector(phi_rad, "phi_rad")
    slave = _matrix(q_slave, "q_slave")
    smooth = _boolean_vector(branch_smooth, "branch_smooth")
    applicable = _boolean_vector(
        differential_applicable, "differential_applicable"
    )
    differential = _boolean_vector(
        differential_valid, "differential_valid"
    )
    differential_runtime = _vector(
        differential_runtime_s, "differential_runtime_s"
    )
    differential_consistency = _vector(
        differential_max_one_sided_consistency,
        "differential_max_one_sided_consistency",
    )
    _same_rows(
        {
            "pose_runtime_s": pose_runtime,
            "warm_start": warm,
            "phi_rad": phi,
            "q_slave": slave,
            "branch_smooth": smooth,
            "differential_applicable": applicable,
            "differential_valid": differential,
            "differential_runtime_s": differential_runtime,
            "differential_max_one_sided_consistency": (
                differential_consistency
            ),
        }
    )
    for array, name in (
        (pose_runtime, "pose_runtime_s"),
        (phi, "phi_rad"),
        (slave, "q_slave"),
        (differential_runtime, "differential_runtime_s"),
    ):
        _finite(array, name)
    valid_differential_mask = applicable & differential
    if np.any(valid_differential_mask) and not np.all(
        np.isfinite(
            differential_consistency[valid_differential_mask]
        )
    ):
        raise MetricError(
            "valid differential samples need finite consistency evidence"
        )
    if np.any(pose_runtime < 0.0) or np.any(differential_runtime < 0.0):
        raise MetricError("H1 runtime samples must be non-negative")
    deadline_s = float(deadline_s)
    if not np.isfinite(deadline_s) or deadline_s <= 0.0:
        raise MetricError("H1 deadline_s must be finite and positive")
    minimum_phi_wrap_crossings = int(minimum_phi_wrap_crossings)
    if minimum_phi_wrap_crossings < 0:
        raise MetricError(
            "minimum_phi_wrap_crossings must be non-negative"
        )

    pose_ms = 1e3 * pose_runtime
    warm_pose_ms = pose_ms[warm]
    raw_phi_step = np.abs(np.diff(phi))
    wrap_crossing = raw_phi_step > np.pi
    slave_step = (
        np.max(np.abs(_wrapped_delta(slave[1:] - slave[:-1])), axis=1)
        if slave.shape[0] > 1
        else np.empty(0, dtype=float)
    )
    crossing_indices = np.flatnonzero(wrap_crossing) + 1
    crossing_count = int(np.sum(wrap_crossing))
    wrap_metric_valid = crossing_count >= minimum_phi_wrap_crossings
    result: dict[str, Any] = {
        "h1_pose_runtime_p50_ms": float(np.percentile(pose_ms, 50)),
        "h1_pose_runtime_p95_ms": float(np.percentile(pose_ms, 95)),
        "h1_pose_runtime_p99_ms": float(np.percentile(pose_ms, 99)),
        "h1_pose_runtime_max_ms": float(np.max(pose_ms)),
        "h1_pose_runtime_warm_p95_ms": (
            float(np.percentile(warm_pose_ms, 95))
            if warm_pose_ms.size
            else None
        ),
        "h1_pose_deadline_ms": 1e3 * deadline_s,
        "h1_pose_deadline_miss_count": int(
            np.sum(pose_runtime > deadline_s)
        ),
        "h1_pose_deadline_miss_fraction": float(
            np.mean(pose_runtime > deadline_s)
        ),
        "h1_branch_smooth_fraction": float(np.mean(smooth)),
        "h1_phi_raw_wrap_crossing_count": crossing_count,
        "h1_phi_raw_wrap_crossing_indices": crossing_indices.tolist(),
        "h1_phi_wrap_metric_valid": wrap_metric_valid,
        "h1_phi_wrap_minimum_crossing_count": (
            minimum_phi_wrap_crossings
        ),
        "h1_phi_wrap_crossing_max_slave_joint_step_rad": (
            float(np.max(slave_step[wrap_crossing]))
            if np.any(wrap_crossing)
            else None
        ),
        "h1_differential_applicable_fraction": float(
            np.mean(applicable)
        ),
    }
    if np.any(applicable):
        selected_runtime = differential_runtime[applicable]
        selected_consistency = differential_consistency[applicable]
        finite_consistency = selected_consistency[
            np.isfinite(selected_consistency)
        ]
        differential_ms = 1e3 * selected_runtime
        feedback_ready_mask = applicable & differential
        feedback_ready_ms = 1e3 * (
            pose_runtime[feedback_ready_mask]
            + differential_runtime[feedback_ready_mask]
        )
        result.update(
            {
                "h1_differential_valid_fraction": float(
                    np.mean(differential[applicable])
                ),
                "h1_differential_runtime_p50_ms": float(
                    np.percentile(differential_ms, 50)
                ),
                "h1_differential_runtime_p95_ms": float(
                    np.percentile(differential_ms, 95)
                ),
                "h1_differential_runtime_p99_ms": float(
                    np.percentile(differential_ms, 99)
                ),
                "h1_differential_runtime_max_ms": float(
                    np.max(differential_ms)
                ),
                "h1_differential_max_one_sided_consistency": float(
                    np.max(finite_consistency)
                )
                if finite_consistency.size
                else None,
                "h1_feedback_ready_valid_sample_count": int(
                    np.sum(feedback_ready_mask)
                ),
                "h1_feedback_ready_runtime_p95_ms": (
                    float(np.percentile(feedback_ready_ms, 95))
                    if feedback_ready_ms.size
                    else None
                ),
                "h1_feedback_ready_deadline_miss_count": (
                    int(np.sum(feedback_ready_ms > 1e3 * deadline_s))
                    if feedback_ready_ms.size
                    else None
                ),
                "h1_feedback_ready_deadline_miss_fraction": (
                    float(
                        np.mean(
                            feedback_ready_ms > 1e3 * deadline_s
                        )
                    )
                    if feedback_ready_ms.size
                    else None
                ),
            }
        )
    else:
        result.update(
            {
                "h1_differential_valid_fraction": None,
                "h1_differential_runtime_p50_ms": None,
                "h1_differential_runtime_p95_ms": None,
                "h1_differential_runtime_p99_ms": None,
                "h1_differential_runtime_max_ms": None,
                "h1_differential_max_one_sided_consistency": None,
                "h1_feedback_ready_valid_sample_count": 0,
                "h1_feedback_ready_runtime_p95_ms": None,
                "h1_feedback_ready_deadline_miss_count": None,
                "h1_feedback_ready_deadline_miss_fraction": None,
            }
        )
    return result


def compute_h1_audit_metrics(
    *,
    differential_valid: Any,
    validity_reason_code: Any,
    reach_clip_code: Any,
    joint_limit_active: Any,
    geometry_degenerate: Any,
    low_manipulability: Any,
    slave_min_singular_value: Any,
    slave_manipulability: Any,
    validity_reason_labels: Any | None = None,
) -> dict[str, Any]:
    """Summarize H1 stratum isolation and invalid-sample explanations."""
    differential = _boolean_vector(
        differential_valid, "differential_valid"
    )
    reason_raw = _vector(validity_reason_code, "validity_reason_code")
    reach_raw = _vector(reach_clip_code, "reach_clip_code")
    joint_limit = _boolean_vector(
        joint_limit_active, "joint_limit_active"
    )
    geometry = _boolean_vector(
        geometry_degenerate, "geometry_degenerate"
    )
    low_manip = _boolean_vector(
        low_manipulability, "low_manipulability"
    )
    minimum_singular = _vector(
        slave_min_singular_value, "slave_min_singular_value"
    )
    manipulability = _vector(slave_manipulability, "slave_manipulability")
    _same_rows(
        {
            "differential_valid": differential,
            "validity_reason_code": reason_raw,
            "reach_clip_code": reach_raw,
            "joint_limit_active": joint_limit,
            "geometry_degenerate": geometry,
            "low_manipulability": low_manip,
            "slave_min_singular_value": minimum_singular,
            "slave_manipulability": manipulability,
        }
    )
    _finite(reason_raw, "validity_reason_code")
    _finite(reach_raw, "reach_clip_code")
    _finite(minimum_singular, "slave_min_singular_value")
    _finite(manipulability, "slave_manipulability")
    if np.any(minimum_singular < 0.0) or np.any(manipulability < 0.0):
        raise MetricError(
            "H1 singular-value and manipulability diagnostics must be non-negative"
        )
    if not np.all(reason_raw == np.floor(reason_raw)) or np.any(reason_raw < 0):
        raise MetricError("validity_reason_code must contain non-negative integers")
    if not np.all(reach_raw == np.floor(reach_raw)) or not set(
        np.asarray(reach_raw, dtype=int).tolist()
    ).issubset({-1, 0, 1}):
        raise MetricError("reach_clip_code must contain only -1, 0, or 1")
    reason = np.asarray(reason_raw, dtype=int)
    reach = np.asarray(reach_raw, dtype=int)

    if validity_reason_labels is None:
        labels: list[str] = []
    elif (
        not isinstance(validity_reason_labels, list)
        or any(not isinstance(label, str) or not label for label in validity_reason_labels)
    ):
        raise MetricError("H1 validity_reason_labels must be a list of strings")
    else:
        labels = list(validity_reason_labels)
    maximum_code = int(np.max(reason))
    if labels and maximum_code >= len(labels):
        raise MetricError(
            "H1 validity_reason_labels does not cover every recorded code"
        )

    unique_codes, counts = np.unique(reason, return_counts=True)
    histogram = {
        (labels[int(code)] if labels else f"code_{int(code)}"): int(count)
        for code, count in zip(unique_codes, counts)
    }
    invalid = ~differential
    invalid_codes = reason[invalid]
    if invalid_codes.size:
        invalid_unique, invalid_counts = np.unique(
            invalid_codes, return_counts=True
        )
        highest = int(np.max(invalid_counts))
        # Ties use the lowest stable enum code.
        primary_code = int(
            np.min(invalid_unique[invalid_counts == highest])
        )
        primary_reason = (
            labels[primary_code] if labels else f"code_{primary_code}"
        )
    else:
        primary_reason = labels[0] if labels else "code_0"

    return {
        "h1_reach_clip_fraction": float(np.mean(reach != 0)),
        "h1_reach_clip_lower_fraction": float(np.mean(reach == -1)),
        "h1_reach_clip_upper_fraction": float(np.mean(reach == 1)),
        "h1_low_manipulability_fraction": float(np.mean(low_manip)),
        "h1_joint_limit_active_fraction": float(np.mean(joint_limit)),
        "h1_geometry_degenerate_fraction": float(np.mean(geometry)),
        "h1_min_slave_min_singular_value": float(np.min(minimum_singular)),
        "h1_min_slave_manipulability": float(np.min(manipulability)),
        "h1_validity_reason_histogram": histogram,
        "h1_primary_invalid_reason": primary_reason,
        "h1_invalid_reason_sample_count": int(invalid_codes.size),
        "h1_unexplained_invalid_fraction": float(
            np.mean(invalid & (reason == 0))
        ),
    }


def compute_h2_wrench_metrics(
    wrench_estimated: Any,
    wrench_reference: Any,
    *,
    sample_mask: Any | None = None,
    numerical_rank_deficient: Any | None = None,
    operationally_ill_conditioned: Any | None = None,
    numerical_rank_threshold: Any | None = None,
    operational_min_scaled_singular_threshold: Any | None = None,
    scaled_singular_values: Any | None = None,
    condition_number: Any | None = None,
) -> dict[str, Any]:
    estimate = _matrix(wrench_estimated, "wrench_estimated", columns=6)
    reference = _matrix(wrench_reference, "wrench_reference", columns=6)
    n = _same_rows({"wrench_estimated": estimate, "wrench_reference": reference})
    _finite(estimate, "wrench_estimated")
    _finite(reference, "wrench_reference")
    mask = (
        np.ones(n, dtype=bool)
        if sample_mask is None
        else _boolean_vector(sample_mask, "sample_mask")
    )
    if mask.shape[0] != n or not np.any(mask):
        raise MetricError("H2 sample_mask must select at least one aligned sample")
    error = estimate[mask] - reference[mask]
    force_norm = np.linalg.norm(error[:, :3], axis=1)
    moment_norm = np.linalg.norm(error[:, 3:], axis=1)
    result = {
        "h2_sample_count": int(error.shape[0]),
        "h2_force_rmse_N": float(np.sqrt(np.mean(force_norm**2))),
        "h2_moment_rmse_Nm": float(np.sqrt(np.mean(moment_norm**2))),
        "h2_force_mae_N": float(np.mean(force_norm)),
        "h2_moment_mae_Nm": float(np.mean(moment_norm)),
        "h2_force_bias_xyz_N": np.mean(error[:, :3], axis=0).tolist(),
        "h2_moment_bias_xyz_Nm": np.mean(error[:, 3:], axis=0).tolist(),
    }
    if numerical_rank_deficient is not None:
        numerical_flag = _boolean_vector(
            numerical_rank_deficient,
            "numerical_rank_deficient",
        )
        if numerical_flag.shape[0] != n:
            raise MetricError("H2 numerical-rank flag length mismatch")
        result["h2_numerical_rank_deficient_fraction"] = float(
            np.mean(numerical_flag[mask])
        )
    if operationally_ill_conditioned is not None:
        operational_flag = _boolean_vector(
            operationally_ill_conditioned,
            "operationally_ill_conditioned",
        )
        if operational_flag.shape[0] != n:
            raise MetricError("H2 operational-condition flag length mismatch")
        result["h2_operationally_ill_conditioned_fraction"] = float(
            np.mean(operational_flag[mask])
        )
    if numerical_rank_threshold is not None:
        numerical_threshold = _vector(
            numerical_rank_threshold,
            "numerical_rank_threshold",
        )
        if numerical_threshold.shape[0] != n:
            raise MetricError("H2 numerical-rank threshold length mismatch")
        _finite(numerical_threshold, "numerical_rank_threshold")
        result["h2_numerical_rank_threshold_min"] = float(
            np.min(numerical_threshold[mask])
        )
        result["h2_numerical_rank_threshold_max"] = float(
            np.max(numerical_threshold[mask])
        )
    if operational_min_scaled_singular_threshold is not None:
        operational_threshold = _vector(
            operational_min_scaled_singular_threshold,
            "operational_min_scaled_singular_threshold",
        )
        if operational_threshold.shape[0] != n:
            raise MetricError("H2 operational threshold length mismatch")
        _finite(
            operational_threshold,
            "operational_min_scaled_singular_threshold",
        )
        selected_threshold = operational_threshold[mask]
        if not np.allclose(
            selected_threshold,
            selected_threshold[0],
            rtol=0.0,
            atol=1e-15,
        ):
            raise MetricError(
                "H2 operational threshold must be constant within a trial"
            )
        result["h2_operational_min_scaled_singular_threshold"] = float(
            selected_threshold[0]
        )
    if scaled_singular_values is not None:
        singular_values = _matrix(
            scaled_singular_values,
            "scaled_singular_values",
        )
        if singular_values.shape[0] != n:
            raise MetricError("H2 singular-value sample count mismatch")
        _finite(singular_values, "scaled_singular_values")
        result["h2_min_scaled_singular_value"] = float(
            np.min(singular_values[mask, -1])
        )
    if condition_number is not None:
        condition = _vector(condition_number, "condition_number")
        if condition.shape[0] != n:
            raise MetricError("H2 condition-number sample count mismatch")
        if np.any(np.isnan(condition)) or np.any(condition < 0.0):
            raise MetricError(
                "H2 condition number must be non-negative and not NaN"
            )
        selected_condition = condition[mask]
        finite_condition = selected_condition[np.isfinite(selected_condition)]
        result["h2_infinite_condition_number_fraction"] = float(
            np.mean(~np.isfinite(selected_condition))
        )
        result["h2_max_finite_condition_number"] = (
            float(np.max(finite_condition))
            if finite_condition.size
            else None
        )
        # A rank-deficient sample has infinite condition number.  JSON cannot
        # represent infinity, so the explicit fraction above carries that
        # condition while this field is null in that case.
        result["h2_max_condition_number"] = (
            float(np.max(selected_condition))
            if np.all(np.isfinite(selected_condition))
            else None
        )
    return result


def _dt_array(dt: Any, count: int) -> np.ndarray:
    array = np.asarray(dt, dtype=float)
    if array.ndim == 0:
        array = np.full(count, float(array), dtype=float)
    if array.shape != (count,):
        raise MetricError(f"dt must be scalar or have shape ({count},)")
    if not np.all(np.isfinite(array)) or np.any(array <= 0.0):
        raise MetricError("dt must contain finite positive intervals")
    return array


def compute_h3_power_mismatch(
    *,
    tau_master_raw: Any,
    qd_master: Any,
    tau_slave_source: Any,
    qd_slave_source: Any,
    dt: Any,
    force_scale: float = 1.0,
    epsilon_energy_J: float = 1e-12,
    minimum_power_activity_J: float = 0.0,
    return_valid: Any | None = None,
) -> dict[str, Any]:
    tau_m = _matrix(tau_master_raw, "tau_master_raw")
    qd_m = _matrix(qd_master, "qd_master")
    tau_s = _matrix(tau_slave_source, "tau_slave_source")
    qd_s = _matrix(qd_slave_source, "qd_slave_source")
    n = _same_rows(
        {
            "tau_master_raw": tau_m,
            "qd_master": qd_m,
            "tau_slave_source": tau_s,
            "qd_slave_source": qd_s,
        }
    )
    if tau_m.shape != qd_m.shape or tau_s.shape != qd_s.shape:
        raise MetricError("torque and velocity shapes must match at each port")
    for array, name in (
        (tau_m, "tau_master_raw"),
        (qd_m, "qd_master"),
        (tau_s, "tau_slave_source"),
        (qd_s, "qd_slave_source"),
    ):
        _finite(array, name)
    intervals = _dt_array(dt, n)
    chi = (
        np.ones(n, dtype=bool)
        if return_valid is None
        else _boolean_vector(return_valid, "return_valid")
    )
    if chi.shape[0] != n:
        raise MetricError("return_valid length mismatch")
    force_scale = float(force_scale)
    epsilon_energy_J = float(epsilon_energy_J)
    minimum_power_activity_J = float(minimum_power_activity_J)
    if not np.isfinite(force_scale) or force_scale < 0.0:
        raise MetricError("force_scale must be finite and non-negative")
    if not np.isfinite(epsilon_energy_J) or epsilon_energy_J <= 0.0:
        raise MetricError("epsilon_energy_J must be finite and positive")
    if (
        not np.isfinite(minimum_power_activity_J)
        or minimum_power_activity_J < 0.0
    ):
        raise MetricError(
            "minimum_power_activity_J must be finite and non-negative"
        )

    power_master = np.einsum("ij,ij->i", tau_m, qd_m)
    power_slave = chi.astype(float) * np.einsum("ij,ij->i", tau_s, qd_s)
    scaled_slave = force_scale * power_slave
    numerator = float(np.sum(np.abs(power_master - scaled_slave) * intervals))
    power_activity = float(
        0.5
        * np.sum((np.abs(power_master) + np.abs(scaled_slave)) * intervals)
    )
    denominator = power_activity + epsilon_energy_J
    normalized = numerator / denominator
    normalized_valid = bool(power_activity >= minimum_power_activity_J)
    return {
        # Kept for backward-compatible diagnostics.  Confirmatory analysis must
        # use the validity flag/gated value when a nonzero activity gate is set.
        "h3_epsilon_P_act": normalized,
        "h3_epsilon_P_act_gated": normalized if normalized_valid else None,
        "h3_normalized_metric_valid": normalized_valid,
        "h3_minimum_power_activity_J": minimum_power_activity_J,
        "h3_power_activity_J": power_activity,
        "h3_absolute_power_mismatch_J": numerator,
        "h3_power_mismatch_numerator_J": numerator,
        "h3_power_normalizer_J": denominator,
        "h3_return_valid_fraction": float(np.mean(chi)),
        "h3_master_raw_work_J": float(np.sum(power_master * intervals)),
        "h3_scaled_slave_work_J": float(np.sum(scaled_slave * intervals)),
    }


def audit_h4_energy(
    *,
    energy_before_J: Any,
    energy_after_J: Any,
    tau_candidate: Any,
    tau_applied: Any,
    tau_accepted: Any | None = None,
    qd_master: Any,
    dt: Any,
    energy_min_J: float,
    energy_max_J: float,
    epsilon_torque_impulse_Nms: float = 1e-12,
    audit_tolerance_J: float = 1e-10,
    software_preclip_J: Any | None = None,
) -> dict[str, Any]:
    energy_before = _vector(energy_before_J, "energy_before_J")
    energy_after = _vector(energy_after_J, "energy_after_J")
    candidate = _matrix(tau_candidate, "tau_candidate")
    applied = _matrix(tau_applied, "tau_applied")
    accepted = (
        applied
        if tau_accepted is None
        else _matrix(tau_accepted, "tau_accepted")
    )
    velocity = _matrix(qd_master, "qd_master")
    n = _same_rows(
        {
            "energy_before_J": energy_before,
            "energy_after_J": energy_after,
            "tau_candidate": candidate,
            "tau_applied": applied,
            "tau_accepted": accepted,
            "qd_master": velocity,
        }
    )
    if not (
        candidate.shape
        == applied.shape
        == accepted.shape
        == velocity.shape
    ):
        raise MetricError(
            "candidate, projected, accepted torque, and velocity must align"
        )
    for array, name in (
        (energy_before, "energy_before_J"),
        (energy_after, "energy_after_J"),
        (candidate, "tau_candidate"),
        (applied, "tau_applied"),
        (accepted, "tau_accepted"),
        (velocity, "qd_master"),
    ):
        _finite(array, name)
    intervals = _dt_array(dt, n)
    energy_min_J = float(energy_min_J)
    energy_max_J = float(energy_max_J)
    tolerance = float(audit_tolerance_J)
    epsilon_tau = float(epsilon_torque_impulse_Nms)
    if not (
        np.isfinite(energy_min_J)
        and np.isfinite(energy_max_J)
        and 0.0 <= energy_min_J <= energy_max_J
    ):
        raise MetricError("invalid energy bounds")
    if not np.isfinite(tolerance) or tolerance < 0.0:
        raise MetricError("audit_tolerance_J must be finite and non-negative")
    if not np.isfinite(epsilon_tau) or epsilon_tau <= 0.0:
        raise MetricError(
            "epsilon_torque_impulse_Nms must be finite and positive"
        )

    # The deterministic gate is reconstructed at the actuator-accepted port.
    # If no drive readback exists, callers may omit tau_accepted and explicitly
    # declare that projected == accepted for that backend.
    applied_power = np.einsum("ij,ij->i", accepted, velocity)
    candidate_power = np.einsum("ij,ij->i", candidate, velocity)
    reconstructed_preclip = energy_before - applied_power * intervals
    reconstructed_after = np.clip(
        reconstructed_preclip,
        energy_min_J,
        energy_max_J,
    )
    deficit = np.maximum(0.0, energy_min_J - reconstructed_preclip)
    accounting_error = np.abs(energy_after - reconstructed_after)

    shadow_energy = float(energy_before[0])
    shadow_min = shadow_energy
    for power, interval in zip(candidate_power, intervals):
        shadow_energy = min(energy_max_J, shadow_energy - float(power * interval))
        shadow_min = min(shadow_min, shadow_energy)
    shadow_deficit = max(0.0, energy_min_J - shadow_min)
    projected_deficit = float(np.max(deficit))

    numerator = float(
        np.sum(np.linalg.norm(accepted - candidate, axis=1) * intervals)
    )
    denominator = float(
        np.sum(np.linalg.norm(candidate, axis=1) * intervals) + epsilon_tau
    )
    preclip_discrepancy = 0.0
    if software_preclip_J is not None:
        software = _vector(software_preclip_J, "software_preclip_J")
        if software.shape[0] != n:
            raise MetricError("software_preclip_J length mismatch")
        _finite(software, "software_preclip_J")
        preclip_discrepancy = float(
            np.max(np.abs(software - reconstructed_preclip))
        )
    audit_pass = (
        projected_deficit <= tolerance
        and float(np.max(accounting_error)) <= tolerance
        and preclip_discrepancy <= tolerance
    )
    return {
        "h4_energy_audit_pass": bool(audit_pass),
        "h4_preclip_floor_deficit_max_J": projected_deficit,
        "h4_energy_accounting_max_error_J": float(np.max(accounting_error)),
        "h4_software_preclip_max_error_J": preclip_discrepancy,
        "h4_shadow_floor_deficit_J": float(shadow_deficit),
        "h4_projected_floor_deficit_J": projected_deficit,
        "h4_delta_B_J": float(shadow_deficit - projected_deficit),
        "h4_D_proj": numerator / denominator,
        "h4_energy_min_J": energy_min_J,
        "h4_energy_max_J": energy_max_J,
        "h4_projection_distortion_numerator_Nms": numerator,
        "h4_projection_distortion_normalizer_Nms": denominator,
        "h4_shadow_energy_min_J": float(shadow_min),
        "h4_downstream_modification_max_Nm": float(
            np.max(np.linalg.norm(accepted - applied, axis=1))
        ),
    }


def compute_bilateral_diagnostics(
    *,
    master_tracking_error_rad: Any,
    slave_tracking_error_rad: Any,
    feedback_torque_Nm: Any,
    contact_force_N: Any,
    projection_factor: Any,
    wall_force_raw_N: Any | None = None,
    wall_force_applied_N: Any | None = None,
    wall_force_saturation_active: Any | None = None,
    wall_force_limit_N: Any | None = None,
    master_joint_limit_active: Any | None = None,
    slave_joint_limit_active: Any | None = None,
    master_velocity_limit_active: Any | None = None,
    slave_velocity_limit_active: Any | None = None,
    master_acceleration_limit_active: Any | None = None,
    slave_acceleration_limit_active: Any | None = None,
    master_torque_saturation_active: Any | None = None,
    slave_torque_saturation_active: Any | None = None,
    haptic_rate_limit_active: Any | None = None,
    haptic_torque_saturation_active: Any | None = None,
    energy_probe_raw_work_J: Any | None = None,
    projection_tolerance: float = 1e-12,
    contact_force_threshold_N: float = 1e-6,
    minimum_contact_fraction: float = 0.0,
    minimum_contact_rms_N: float = 0.0,
    maximum_force_limit_hit_fraction: float = 1.0,
    minimum_force_headroom_N: float = 0.0,
    maximum_master_tracking_rmse_rad: float | None = None,
    maximum_slave_tracking_rmse_rad: float | None = None,
    minimum_projection_intervention_fraction: float = 0.0,
    maximum_projection_intervention_fraction: float = 1.0,
    maximum_limit_active_fraction: float = 1.0,
    minimum_energy_probe_raw_work_J: float = 0.0,
) -> dict[str, Any]:
    """Compute secondary task/transparency diagnostics from stored samples."""
    master_error = _vector(
        master_tracking_error_rad, "master_tracking_error_rad"
    )
    slave_error = _vector(
        slave_tracking_error_rad, "slave_tracking_error_rad"
    )
    feedback = _matrix(feedback_torque_Nm, "feedback_torque_Nm")
    contact = _vector(contact_force_N, "contact_force_N")
    rho = _vector(projection_factor, "projection_factor")
    applied = (
        contact.copy()
        if wall_force_applied_N is None
        else _vector(wall_force_applied_N, "wall_force_applied_N")
    )
    raw = (
        applied.copy()
        if wall_force_raw_N is None
        else _vector(wall_force_raw_N, "wall_force_raw_N")
    )
    sample_count = contact.shape[0]
    if wall_force_limit_N is None:
        force_limit = None
    else:
        force_limit_array = np.asarray(wall_force_limit_N, dtype=float)
        if force_limit_array.ndim == 0:
            force_limit = np.full(sample_count, float(force_limit_array))
        else:
            force_limit = _vector(
                force_limit_array, "wall_force_limit_N"
            )
    if wall_force_saturation_active is None:
        saturation = (
            np.zeros(sample_count, dtype=bool)
            if force_limit is None
            else raw > force_limit
        )
    else:
        saturation = _boolean_vector(
            wall_force_saturation_active,
            "wall_force_saturation_active",
        )
    limit_inputs = {
        "master_joint_limit_active": master_joint_limit_active,
        "slave_joint_limit_active": slave_joint_limit_active,
        "master_velocity_limit_active": master_velocity_limit_active,
        "slave_velocity_limit_active": slave_velocity_limit_active,
        "master_acceleration_limit_active": master_acceleration_limit_active,
        "slave_acceleration_limit_active": slave_acceleration_limit_active,
        "master_torque_saturation_active": master_torque_saturation_active,
        "slave_torque_saturation_active": slave_torque_saturation_active,
        "haptic_rate_limit_active": haptic_rate_limit_active,
        "haptic_torque_saturation_active": (
            haptic_torque_saturation_active
        ),
    }
    limit_flags = {
        name: (
            np.zeros(sample_count, dtype=bool)
            if value is None
            else _boolean_vector(value, name)
        )
        for name, value in limit_inputs.items()
    }
    probe_work = (
        np.zeros(sample_count, dtype=float)
        if energy_probe_raw_work_J is None
        else _vector(
            energy_probe_raw_work_J,
            "energy_probe_raw_work_J",
        )
    )
    _same_rows(
        {
            "master_tracking_error_rad": master_error,
            "slave_tracking_error_rad": slave_error,
            "feedback_torque_Nm": feedback,
            "contact_force_N": contact,
            "projection_factor": rho,
            "wall_force_applied_N": applied,
            "wall_force_raw_N": raw,
            "wall_force_saturation_active": saturation,
            **(
                {}
                if force_limit is None
                else {"wall_force_limit_N": force_limit}
            ),
            **limit_flags,
            "energy_probe_raw_work_J": probe_work,
        }
    )
    for array, name in (
        (master_error, "master_tracking_error_rad"),
        (slave_error, "slave_tracking_error_rad"),
        (feedback, "feedback_torque_Nm"),
        (contact, "contact_force_N"),
        (rho, "projection_factor"),
        (applied, "wall_force_applied_N"),
        (raw, "wall_force_raw_N"),
        (probe_work, "energy_probe_raw_work_J"),
    ):
        _finite(array, name)
    nonnegative_arrays = (
        (contact, "contact_force_N"),
        (applied, "wall_force_applied_N"),
        (raw, "wall_force_raw_N"),
    )
    if any(np.any(array < 0.0) for array, _ in nonnegative_arrays):
        names = ", ".join(name for _, name in nonnegative_arrays)
        raise MetricError(f"{names} must be non-negative")
    audit_tolerance_N = 1e-12
    if np.any(raw + audit_tolerance_N < applied):
        raise MetricError(
            "wall_force_raw_N cannot be smaller than applied force"
        )
    if not np.allclose(
        contact,
        applied,
        rtol=0.0,
        atol=audit_tolerance_N,
    ):
        raise MetricError(
            "contact_force_N must equal the applied wall-force magnitude"
        )
    if np.any(np.diff(probe_work) < -1e-12):
        raise MetricError(
            "energy_probe_raw_work_J must be cumulative and non-decreasing"
        )
    if force_limit is not None:
        _finite(force_limit, "wall_force_limit_N")
        if np.any(force_limit <= 0.0):
            raise MetricError("wall_force_limit_N must be positive")
        expected_applied = np.minimum(raw, force_limit)
        expected_saturation = raw > force_limit
        if not np.allclose(
            applied,
            expected_applied,
            rtol=0.0,
            atol=audit_tolerance_N,
        ):
            raise MetricError(
                "applied wall force disagrees with raw force/force limit"
            )
        if not np.array_equal(saturation, expected_saturation):
            raise MetricError(
                "wall-force saturation flags disagree with raw force/limit"
            )
    projection_tolerance = float(projection_tolerance)
    contact_force_threshold_N = float(contact_force_threshold_N)
    if not np.isfinite(projection_tolerance) or projection_tolerance < 0.0:
        raise MetricError("projection_tolerance must be finite and non-negative")
    if (
        not np.isfinite(contact_force_threshold_N)
        or contact_force_threshold_N < 0.0
    ):
        raise MetricError(
            "contact_force_threshold_N must be finite and non-negative"
        )
    gate_values = {
        "minimum_contact_fraction": minimum_contact_fraction,
        "minimum_contact_rms_N": minimum_contact_rms_N,
        "maximum_force_limit_hit_fraction": (
            maximum_force_limit_hit_fraction
        ),
        "minimum_force_headroom_N": minimum_force_headroom_N,
        "maximum_projection_intervention_fraction": (
            maximum_projection_intervention_fraction
        ),
        "minimum_projection_intervention_fraction": (
            minimum_projection_intervention_fraction
        ),
        "maximum_limit_active_fraction": maximum_limit_active_fraction,
        "minimum_energy_probe_raw_work_J":
            minimum_energy_probe_raw_work_J,
    }
    if any(
        not np.isfinite(float(value)) or float(value) < 0.0
        for value in gate_values.values()
    ):
        raise MetricError("bilateral gate thresholds must be finite/non-negative")
    if (
        float(minimum_contact_fraction) > 1.0
        or float(maximum_force_limit_hit_fraction) > 1.0
        or float(maximum_projection_intervention_fraction) > 1.0
        or float(minimum_projection_intervention_fraction) > 1.0
        or float(maximum_limit_active_fraction) > 1.0
    ):
        raise MetricError("bilateral fraction thresholds must lie in [0, 1]")
    if (
        float(minimum_projection_intervention_fraction)
        > float(maximum_projection_intervention_fraction)
    ):
        raise MetricError(
            "minimum projection fraction cannot exceed its maximum"
        )
    for value, name in (
        (
            maximum_master_tracking_rmse_rad,
            "maximum_master_tracking_rmse_rad",
        ),
        (
            maximum_slave_tracking_rmse_rad,
            "maximum_slave_tracking_rmse_rad",
        ),
    ):
        if value is not None and (
            not np.isfinite(float(value)) or float(value) < 0.0
        ):
            raise MetricError(f"{name} must be finite and non-negative")
    feedback_norm = np.linalg.norm(feedback, axis=1)
    master_tracking_rmse = float(
        np.sqrt(np.mean(np.square(master_error)))
    )
    slave_tracking_rmse = float(
        np.sqrt(np.mean(np.square(slave_error)))
    )
    projection_intervention_fraction = float(
        np.mean(rho < (1.0 - projection_tolerance))
    )
    any_limit_active = np.logical_or.reduce(
        tuple(limit_flags.values())
    )
    limit_active_fraction = float(np.mean(any_limit_active))
    contact_mask = applied > contact_force_threshold_N
    contact_fraction = float(np.mean(contact_mask))
    contact_rms_all = float(np.sqrt(np.mean(np.square(applied))))
    contact_rms_active = (
        float(np.sqrt(np.mean(np.square(applied[contact_mask]))))
        if np.any(contact_mask)
        else 0.0
    )
    force_limit_hit_fraction = float(np.mean(saturation))
    force_limit_hit_contact_fraction = (
        float(np.mean(saturation[contact_mask]))
        if np.any(contact_mask)
        else 0.0
    )
    force_headroom_N = (
        None
        if force_limit is None
        else float(np.min(force_limit - applied))
    )
    contact_fraction_gate = (
        contact_fraction >= float(minimum_contact_fraction)
    )
    contact_rms_gate = (
        contact_rms_active >= float(minimum_contact_rms_N)
    )
    force_limit_gate = (
        force_limit_hit_fraction
        <= float(maximum_force_limit_hit_fraction)
    )
    force_headroom_gate = (
        True
        if force_headroom_N is None
        else force_headroom_N >= float(minimum_force_headroom_N)
    )
    master_tracking_gate = (
        True
        if maximum_master_tracking_rmse_rad is None
        else master_tracking_rmse
        <= float(maximum_master_tracking_rmse_rad)
    )
    slave_tracking_gate = (
        True
        if maximum_slave_tracking_rmse_rad is None
        else slave_tracking_rmse
        <= float(maximum_slave_tracking_rmse_rad)
    )
    projection_gate = (
        float(minimum_projection_intervention_fraction)
        <= projection_intervention_fraction
        <= float(maximum_projection_intervention_fraction)
    )
    limit_active_gate = (
        limit_active_fraction <= float(maximum_limit_active_fraction)
    )
    probe_work_final_J = float(probe_work[-1])
    probe_work_gate = (
        probe_work_final_J >= float(minimum_energy_probe_raw_work_J)
    )
    return {
        "bilateral_master_tracking_rmse_rad": master_tracking_rmse,
        "bilateral_slave_tracking_rmse_rad": slave_tracking_rmse,
        "bilateral_feedback_torque_rms_Nm": float(
            np.sqrt(np.mean(np.square(feedback_norm)))
        ),
        "bilateral_feedback_torque_peak_Nm": float(np.max(feedback_norm)),
        "bilateral_contact_force_rms_N": contact_rms_active,
        "bilateral_contact_force_rms_all_samples_N": contact_rms_all,
        "bilateral_contact_force_peak_N": float(np.max(applied)),
        "bilateral_wall_force_raw_peak_N": float(np.max(raw)),
        "bilateral_contact_fraction": contact_fraction,
        "bilateral_force_limit_hit_fraction": force_limit_hit_fraction,
        "bilateral_force_limit_hit_contact_fraction": (
            force_limit_hit_contact_fraction
        ),
        "bilateral_force_headroom_N": force_headroom_N,
        "bilateral_contact_fraction_gate_pass": bool(
            contact_fraction_gate
        ),
        "bilateral_contact_rms_gate_pass": bool(contact_rms_gate),
        "bilateral_force_limit_gate_pass": bool(force_limit_gate),
        "bilateral_force_headroom_gate_pass": bool(
            force_headroom_gate
        ),
        "bilateral_master_tracking_gate_pass": bool(
            master_tracking_gate
        ),
        "bilateral_slave_tracking_gate_pass": bool(
            slave_tracking_gate
        ),
        "bilateral_projection_gate_pass": bool(projection_gate),
        "bilateral_energy_probe_raw_work_J": probe_work_final_J,
        "bilateral_energy_probe_raw_work_gate_pass": bool(
            probe_work_gate
        ),
        "bilateral_limit_active_fraction": limit_active_fraction,
        "bilateral_limit_active_gate_pass": bool(limit_active_gate),
        **{
            f"bilateral_{name}_fraction": float(np.mean(flag))
            for name, flag in limit_flags.items()
        },
        "bilateral_stable_contact_gate_pass": bool(
            contact_fraction_gate
            and contact_rms_gate
            and force_limit_gate
            and force_headroom_gate
            and master_tracking_gate
            and slave_tracking_gate
            and projection_gate
            and limit_active_gate
            and probe_work_gate
        ),
        "bilateral_projection_intervention_fraction": float(
            projection_intervention_fraction
        ),
        "bilateral_projection_factor_min": float(np.min(rho)),
    }


_PACKET_EMPTY = 0
_PACKET_ACTIVE = 1
_PACKET_HELD = 2
_PACKET_TIMED_OUT = 3
_PACKET_RECOVERING = 4


def _packet_state_vector(value: Any, name: str) -> np.ndarray:
    raw = _vector(value, name)
    _finite(raw, name)
    if not np.all(raw == np.floor(raw)) or not np.all(
        (raw >= _PACKET_EMPTY) & (raw <= _PACKET_RECOVERING)
    ):
        raise MetricError(f"{name} must contain only integer states 0..4")
    return raw.astype(np.int64)


def _constant_sample_value(
    value: Any,
    name: str,
    sample_count: int,
) -> tuple[np.ndarray, float]:
    raw = np.asarray(value, dtype=float)
    if raw.ndim == 0:
        array = np.full(sample_count, float(raw))
    else:
        array = _vector(raw, name)
        if array.shape[0] != sample_count:
            raise MetricError(
                f"{name} has {array.shape[0]} samples; expected {sample_count}"
            )
    _finite(array, name)
    constant = float(array[0])
    if not np.allclose(array, constant, rtol=0.0, atol=1e-15):
        raise MetricError(f"{name} must be constant within one trial")
    return array, constant


def _packet_stream_diagnostics(
    *,
    direction: str,
    state: Any,
    active: Any,
    fresh: Any,
    age: Any,
    sequence: Any,
) -> dict[str, Any]:
    state_name = f"{direction}_packet_state"
    active_name = f"{direction}_packet_active"
    fresh_name = f"{direction}_packet_fresh"
    age_name = f"{direction}_packet_age"
    sequence_name = f"{direction}_packet_seq"
    states = _packet_state_vector(state, state_name)
    active_flags = _boolean_vector(active, active_name)
    fresh_flags = _boolean_vector(fresh, fresh_name)
    ages = _vector(age, age_name)
    sequences = _vector(sequence, sequence_name)
    _same_rows(
        {
            state_name: states,
            active_name: active_flags,
            fresh_name: fresh_flags,
            age_name: ages,
            sequence_name: sequences,
        }
    )

    expected_active = np.isin(
        states,
        (_PACKET_ACTIVE, _PACKET_HELD, _PACKET_RECOVERING),
    )
    expected_fresh = np.isin(
        states,
        (_PACKET_ACTIVE, _PACKET_RECOVERING),
    )
    if not np.array_equal(active_flags, expected_active):
        raise MetricError(
            f"{active_name} disagrees with {state_name}; active states are "
            "ACTIVE, HELD, and RECOVERING"
        )
    if not np.array_equal(fresh_flags, expected_fresh):
        raise MetricError(
            f"{fresh_name} disagrees with {state_name}; fresh states are "
            "ACTIVE and RECOVERING"
        )

    if np.any(
        active_flags
        & (~np.isfinite(ages) | (ages < 0.0))
    ):
        raise MetricError(
            f"{age_name} must be finite and non-negative while active"
        )
    if np.any(~active_flags & ~np.isnan(ages)):
        raise MetricError(f"{age_name} must be NaN while inactive")

    active_sequences = sequences[active_flags]
    if np.any(
        ~np.isfinite(active_sequences)
        | (active_sequences < 0.0)
        | (active_sequences != np.floor(active_sequences))
    ):
        raise MetricError(
            f"{sequence_name} must be a finite non-negative integer while active"
        )
    inactive_sequences = sequences[~active_flags]
    inactive_sequence_valid = np.isnan(inactive_sequences) | (
        inactive_sequences == -1.0
    )
    if not np.all(inactive_sequence_valid):
        raise MetricError(
            f"{sequence_name} must be -1 or NaN while inactive"
        )

    previous_sequence: int | None = None
    for index in np.flatnonzero(active_flags):
        current_sequence = int(sequences[index])
        if previous_sequence is None:
            if states[index] == _PACKET_HELD:
                raise MetricError(
                    f"{sequence_name} cannot start with a held packet"
                )
        elif fresh_flags[index]:
            if current_sequence <= previous_sequence:
                raise MetricError(
                    f"{sequence_name} must strictly increase on fresh packets"
                )
        elif current_sequence != previous_sequence:
            raise MetricError(
                f"{sequence_name} must remain constant while a packet is held"
            )
        previous_sequence = current_sequence

    fresh_sequences = sequences[fresh_flags].astype(np.int64)
    if fresh_sequences.size:
        internal_span = int(
            fresh_sequences[-1] - fresh_sequences[0] + 1
        )
        internal_missing_count = internal_span - int(fresh_sequences.size)
        internal_missing_fraction = (
            float(internal_missing_count / internal_span)
            if internal_span > 0
            else 0.0
        )
    else:
        internal_span = 0
        internal_missing_count = 0
        internal_missing_fraction = 0.0

    active_ages = ages[active_flags]
    age_statistics = {
        percentile: (
            float(np.percentile(active_ages, quantile))
            if active_ages.size
            else None
        )
        for percentile, quantile in (
            ("p50", 50),
            ("p95", 95),
            ("p99", 99),
            ("max", 100),
        )
    }
    prefix = f"network_{direction}"
    result = {
        f"{prefix}_active_fraction": float(np.mean(active_flags)),
        f"{prefix}_fresh_fraction": float(np.mean(fresh_flags)),
        f"{prefix}_empty_fraction": float(
            np.mean(states == _PACKET_EMPTY)
        ),
        f"{prefix}_held_fraction": float(
            np.mean(states == _PACKET_HELD)
        ),
        f"{prefix}_timeout_fraction": float(
            np.mean(states == _PACKET_TIMED_OUT)
        ),
        f"{prefix}_recovering_fraction": float(
            np.mean(states == _PACKET_RECOVERING)
        ),
        f"{prefix}_fresh_packet_count": int(fresh_sequences.size),
        f"{prefix}_internal_sequence_span": internal_span,
        f"{prefix}_internal_missing_count": internal_missing_count,
        f"{prefix}_internal_missing_fraction": (
            internal_missing_fraction
        ),
    }
    for percentile, statistic in age_statistics.items():
        result[f"{prefix}_packet_age_{percentile}_s"] = statistic
        # Keep the shorter spelling as an explicit alias for table consumers.
        result[f"{prefix}_age_{percentile}_s"] = statistic
    return result


def compute_network_diagnostics(
    *,
    forward_packet_state: Any,
    return_packet_state: Any,
    forward_packet_active: Any,
    return_packet_active: Any,
    forward_packet_fresh: Any,
    return_packet_fresh: Any,
    forward_packet_age: Any,
    return_packet_age: Any,
    forward_packet_seq: Any,
    return_packet_seq: Any,
    slave_zero_delay_tracking_error: Any,
    slave_reference_lag_error: Any,
    return_feedback_lag_error: Any,
    contact_force_norm: Any,
    contact_expected: Any,
    configured_forward_delay_s: Any,
    configured_return_delay_s: Any,
    configured_forward_jitter_s: Any,
    configured_return_jitter_s: Any,
    configured_forward_packet_loss: Any,
    configured_return_packet_loss: Any,
    configured_forward_timeout_s: Any,
    configured_return_timeout_s: Any,
    contact_force_threshold_N: float = 1e-6,
    maximum_zero_delay_tracking_rmse_rad: float | None = None,
    minimum_active_fraction: float = 0.0,
    minimum_forward_active_fraction: float | None = None,
    minimum_return_active_fraction: float | None = None,
    minimum_forward_fresh_fraction: float = 0.0,
    minimum_return_fresh_fraction: float = 0.0,
    maximum_forward_internal_missing_fraction: float = 1.0,
    maximum_return_internal_missing_fraction: float = 1.0,
    maximum_timeout_fraction: float = 1.0,
    maximum_forward_timeout_fraction: float | None = None,
    maximum_return_timeout_fraction: float | None = None,
    maximum_age_s: float | None = None,
    maximum_forward_age_s: float | None = None,
    maximum_return_age_s: float | None = None,
    minimum_contact_fraction: float = 0.0,
    minimum_contact_rms_N: float = 0.0,
    maximum_free_space_contact_peak_N: float | None = None,
) -> dict[str, Any]:
    """Audit packet-state evidence and compute independent Stage-B endpoints."""
    zero_delay_error = _vector(
        slave_zero_delay_tracking_error,
        "slave_zero_delay_tracking_error",
    )
    reference_lag_error = _vector(
        slave_reference_lag_error,
        "slave_reference_lag_error",
    )
    feedback_lag_error = _vector(
        return_feedback_lag_error,
        "return_feedback_lag_error",
    )
    contact_force = _vector(contact_force_norm, "contact_force_norm")
    expected_contact = _boolean_vector(
        contact_expected, "contact_expected"
    )
    sample_count = _same_rows(
        {
            "slave_zero_delay_tracking_error": zero_delay_error,
            "slave_reference_lag_error": reference_lag_error,
            "return_feedback_lag_error": feedback_lag_error,
            "contact_force_norm": contact_force,
            "contact_expected": expected_contact,
        }
    )
    if not np.all(expected_contact == expected_contact[0]):
        raise MetricError("contact_expected must be constant within one trial")
    for array, name in (
        (zero_delay_error, "slave_zero_delay_tracking_error"),
        (reference_lag_error, "slave_reference_lag_error"),
        (feedback_lag_error, "return_feedback_lag_error"),
        (contact_force, "contact_force_norm"),
    ):
        _finite(array, name)
        if np.any(array < 0.0):
            raise MetricError(f"{name} must be non-negative")

    forward_metrics = _packet_stream_diagnostics(
        direction="forward",
        state=forward_packet_state,
        active=forward_packet_active,
        fresh=forward_packet_fresh,
        age=forward_packet_age,
        sequence=forward_packet_seq,
    )
    return_metrics = _packet_stream_diagnostics(
        direction="return",
        state=return_packet_state,
        active=return_packet_active,
        fresh=return_packet_fresh,
        age=return_packet_age,
        sequence=return_packet_seq,
    )
    packet_lengths = {
        "forward_packet_state": np.asarray(forward_packet_state),
        "return_packet_state": np.asarray(return_packet_state),
        "forward_packet_active": np.asarray(forward_packet_active),
        "return_packet_active": np.asarray(return_packet_active),
    }
    if any(array.shape[0] != sample_count for array in packet_lengths.values()):
        raise MetricError(
            "network packet-state arrays must match tracking sample count"
        )

    configured_inputs = {
        "forward_delay_s": configured_forward_delay_s,
        "return_delay_s": configured_return_delay_s,
        "forward_jitter_s": configured_forward_jitter_s,
        "return_jitter_s": configured_return_jitter_s,
        "forward_packet_loss": configured_forward_packet_loss,
        "return_packet_loss": configured_return_packet_loss,
        "forward_timeout_s": configured_forward_timeout_s,
        "return_timeout_s": configured_return_timeout_s,
    }
    configured: dict[str, float] = {}
    for name, value in configured_inputs.items():
        _, configured[name] = _constant_sample_value(
            value,
            f"configured_{name}",
            sample_count,
        )
    for name in (
        "forward_delay_s",
        "return_delay_s",
        "forward_jitter_s",
        "return_jitter_s",
    ):
        if configured[name] < 0.0:
            raise MetricError(f"configured_{name} must be non-negative")
    for name in ("forward_packet_loss", "return_packet_loss"):
        if not 0.0 <= configured[name] <= 1.0:
            raise MetricError(f"configured_{name} must lie in [0, 1]")
    for name in ("forward_timeout_s", "return_timeout_s"):
        if configured[name] <= 0.0:
            raise MetricError(f"configured_{name} must be positive")

    contact_force_threshold_N = float(contact_force_threshold_N)
    if (
        not np.isfinite(contact_force_threshold_N)
        or contact_force_threshold_N < 0.0
    ):
        raise MetricError(
            "contact_force_threshold_N must be finite and non-negative"
        )

    minimum_forward_active_fraction = float(
        minimum_active_fraction
        if minimum_forward_active_fraction is None
        else minimum_forward_active_fraction
    )
    minimum_return_active_fraction = float(
        minimum_active_fraction
        if minimum_return_active_fraction is None
        else minimum_return_active_fraction
    )
    maximum_forward_timeout_fraction = float(
        maximum_timeout_fraction
        if maximum_forward_timeout_fraction is None
        else maximum_forward_timeout_fraction
    )
    maximum_return_timeout_fraction = float(
        maximum_timeout_fraction
        if maximum_return_timeout_fraction is None
        else maximum_return_timeout_fraction
    )
    maximum_forward_age_s = (
        maximum_age_s
        if maximum_forward_age_s is None
        else maximum_forward_age_s
    )
    maximum_return_age_s = (
        maximum_age_s
        if maximum_return_age_s is None
        else maximum_return_age_s
    )
    fraction_thresholds = {
        "minimum_forward_active_fraction": (
            minimum_forward_active_fraction
        ),
        "minimum_return_active_fraction": (
            minimum_return_active_fraction
        ),
        "minimum_forward_fresh_fraction": float(
            minimum_forward_fresh_fraction
        ),
        "minimum_return_fresh_fraction": float(
            minimum_return_fresh_fraction
        ),
        "maximum_forward_internal_missing_fraction": float(
            maximum_forward_internal_missing_fraction
        ),
        "maximum_return_internal_missing_fraction": float(
            maximum_return_internal_missing_fraction
        ),
        "maximum_forward_timeout_fraction": (
            maximum_forward_timeout_fraction
        ),
        "maximum_return_timeout_fraction": (
            maximum_return_timeout_fraction
        ),
        "minimum_contact_fraction": float(minimum_contact_fraction),
    }
    if any(
        not np.isfinite(value) or not 0.0 <= value <= 1.0
        for value in fraction_thresholds.values()
    ):
        raise MetricError("network fraction gates must lie in [0, 1]")
    nonnegative_optional_thresholds = {
        "maximum_zero_delay_tracking_rmse_rad": (
            maximum_zero_delay_tracking_rmse_rad
        ),
        "maximum_forward_age_s": maximum_forward_age_s,
        "maximum_return_age_s": maximum_return_age_s,
        "maximum_free_space_contact_peak_N": (
            maximum_free_space_contact_peak_N
        ),
    }
    for name, value in nonnegative_optional_thresholds.items():
        if value is not None and (
            not np.isfinite(float(value)) or float(value) < 0.0
        ):
            raise MetricError(f"{name} must be finite and non-negative")
    minimum_contact_rms_N = float(minimum_contact_rms_N)
    if (
        not np.isfinite(minimum_contact_rms_N)
        or minimum_contact_rms_N < 0.0
    ):
        raise MetricError(
            "minimum_contact_rms_N must be finite and non-negative"
        )

    zero_delay_rmse = float(
        np.sqrt(np.mean(np.square(zero_delay_error)))
    )
    reference_lag_rmse = float(
        np.sqrt(np.mean(np.square(reference_lag_error)))
    )
    feedback_lag_rmse = float(
        np.sqrt(np.mean(np.square(feedback_lag_error)))
    )
    contact_mask = contact_force > contact_force_threshold_N
    contact_fraction = float(np.mean(contact_mask))
    contact_rms_all = float(
        np.sqrt(np.mean(np.square(contact_force)))
    )
    contact_rms = (
        float(np.sqrt(np.mean(np.square(contact_force[contact_mask]))))
        if np.any(contact_mask)
        else 0.0
    )
    contact_peak = float(np.max(contact_force))
    contact_is_expected = bool(expected_contact[0])

    forward_active_gate = bool(
        forward_metrics["network_forward_active_fraction"]
        >= minimum_forward_active_fraction
    )
    return_active_gate = bool(
        return_metrics["network_return_active_fraction"]
        >= minimum_return_active_fraction
    )
    forward_fresh_gate = bool(
        forward_metrics["network_forward_fresh_fraction"]
        >= float(minimum_forward_fresh_fraction)
    )
    return_fresh_gate = bool(
        return_metrics["network_return_fresh_fraction"]
        >= float(minimum_return_fresh_fraction)
    )
    forward_internal_missing_gate = bool(
        forward_metrics["network_forward_internal_missing_fraction"]
        <= float(maximum_forward_internal_missing_fraction)
    )
    return_internal_missing_gate = bool(
        return_metrics["network_return_internal_missing_fraction"]
        <= float(maximum_return_internal_missing_fraction)
    )
    forward_timeout_gate = bool(
        forward_metrics["network_forward_timeout_fraction"]
        <= maximum_forward_timeout_fraction
    )
    return_timeout_gate = bool(
        return_metrics["network_return_timeout_fraction"]
        <= maximum_return_timeout_fraction
    )

    def age_gate(direction: str, maximum: float | None) -> bool:
        if maximum is None:
            return True
        age_max = forward_metrics[
            f"network_{direction}_packet_age_max_s"
        ] if direction == "forward" else return_metrics[
            f"network_{direction}_packet_age_max_s"
        ]
        return age_max is not None and age_max <= float(maximum)

    forward_age_gate = age_gate("forward", maximum_forward_age_s)
    return_age_gate = age_gate("return", maximum_return_age_s)
    zero_delay_gate = bool(
        maximum_zero_delay_tracking_rmse_rad is None
        or zero_delay_rmse
        <= float(maximum_zero_delay_tracking_rmse_rad)
    )
    contact_fraction_gate = bool(
        not contact_is_expected
        or contact_fraction >= float(minimum_contact_fraction)
    )
    contact_rms_gate = bool(
        not contact_is_expected
        or contact_rms >= minimum_contact_rms_N
    )
    free_space_peak_gate = bool(
        contact_is_expected
        or maximum_free_space_contact_peak_N is None
        or contact_peak <= float(maximum_free_space_contact_peak_N)
    )
    contact_condition_gate = bool(
        contact_fraction_gate
        and contact_rms_gate
        and free_space_peak_gate
    )
    local_gate = bool(
        zero_delay_gate
        and forward_active_gate
        and return_active_gate
        and forward_fresh_gate
        and return_fresh_gate
        and forward_internal_missing_gate
        and return_internal_missing_gate
        and forward_timeout_gate
        and return_timeout_gate
        and forward_age_gate
        and return_age_gate
        and contact_condition_gate
    )
    return {
        **forward_metrics,
        **return_metrics,
        **{
            f"network_configured_{name}": value
            for name, value in configured.items()
        },
        "network_contact_expected": contact_is_expected,
        "network_slave_zero_delay_tracking_rmse_rad": zero_delay_rmse,
        "network_slave_zero_delay_tracking_p95_rad": float(
            np.percentile(zero_delay_error, 95)
        ),
        "network_slave_zero_delay_tracking_max_rad": float(
            np.max(zero_delay_error)
        ),
        "network_zero_delay_tracking_rmse_rad": zero_delay_rmse,
        "network_zero_delay_tracking_p95_rad": float(
            np.percentile(zero_delay_error, 95)
        ),
        "network_zero_delay_tracking_max_rad": float(
            np.max(zero_delay_error)
        ),
        "network_slave_reference_lag_rmse_rad": reference_lag_rmse,
        "network_slave_reference_lag_p95_rad": float(
            np.percentile(reference_lag_error, 95)
        ),
        "network_slave_reference_lag_max_rad": float(
            np.max(reference_lag_error)
        ),
        "network_return_feedback_lag_rmse_Nm": feedback_lag_rmse,
        "network_return_feedback_lag_p95_Nm": float(
            np.percentile(feedback_lag_error, 95)
        ),
        "network_return_feedback_lag_max_Nm": float(
            np.max(feedback_lag_error)
        ),
        "network_contact_fraction": contact_fraction,
        "network_contact_force_rms_N": contact_rms,
        "network_contact_force_rms_all_samples_N": contact_rms_all,
        "network_contact_force_peak_N": contact_peak,
        "network_zero_delay_tracking_gate_pass": zero_delay_gate,
        "network_forward_active_gate_pass": forward_active_gate,
        "network_return_active_gate_pass": return_active_gate,
        "network_forward_fresh_gate_pass": forward_fresh_gate,
        "network_return_fresh_gate_pass": return_fresh_gate,
        "network_forward_internal_missing_gate_pass": (
            forward_internal_missing_gate
        ),
        "network_return_internal_missing_gate_pass": (
            return_internal_missing_gate
        ),
        "network_forward_timeout_gate_pass": forward_timeout_gate,
        "network_return_timeout_gate_pass": return_timeout_gate,
        "network_forward_age_gate_pass": forward_age_gate,
        "network_return_age_gate_pass": return_age_gate,
        "network_contact_fraction_gate_pass": contact_fraction_gate,
        "network_contact_rms_gate_pass": contact_rms_gate,
        "network_free_space_peak_gate_pass": free_space_peak_gate,
        "network_contact_condition_gate_pass": contact_condition_gate,
        "network_local_gate_pass": local_gate,
    }


def _field(
    samples: Mapping[str, Any],
    fields: Mapping[str, str],
    logical_name: str,
    default_name: str,
    *,
    optional: bool = False,
) -> Any:
    stored_name = fields.get(logical_name, default_name)
    if stored_name not in samples:
        if optional:
            return None
        raise MetricError(
            f"missing sample field {stored_name!r} for {logical_name!r}"
        )
    return samples[stored_name]


def _constant_trial_value(value: Any, name: str) -> float:
    """Resolve a scalar or constant non-empty sample vector."""
    array = np.asarray(value, dtype=float)
    if array.ndim == 0:
        result = float(array)
    elif array.ndim == 1 and array.size > 0:
        _finite(array, name)
        result = float(array[0])
        if not np.allclose(array, result, rtol=0.0, atol=1e-15):
            raise MetricError(f"{name} must be constant within one trial")
    else:
        raise MetricError(f"{name} must be a scalar or non-empty vector")
    if not np.isfinite(result):
        raise MetricError(f"{name} must be finite")
    return result


def _energy_bound_for_trial(
    samples: Mapping[str, Any],
    fields: Mapping[str, str],
    family_config: Mapping[str, Any],
    logical_name: str,
    default_field: str,
) -> float:
    """Prefer stored effective bounds and reject config/sample disagreement."""
    stored = _field(
        samples,
        fields,
        logical_name,
        default_field,
        optional=True,
    )
    configured = family_config.get(logical_name)
    if stored is None and configured is None:
        raise MetricError(
            f"H4 requires {logical_name} in samples or metric configuration"
        )
    if stored is None:
        return _constant_trial_value(configured, logical_name)
    stored_value = _constant_trial_value(stored, default_field)
    if configured is not None:
        configured_value = _constant_trial_value(configured, logical_name)
        if not np.isclose(
            stored_value, configured_value, rtol=0.0, atol=1e-15
        ):
            raise MetricError(
                f"{logical_name} disagrees with stored effective configuration"
            )
    return stored_value


def derive_trial_metrics(
    samples: Mapping[str, Any],
    configuration: Mapping[str, Any],
) -> dict[str, Any]:
    """Derive configured endpoint families from one stored trial."""
    enabled = configuration.get("enabled")
    if not isinstance(enabled, list) or not enabled:
        raise MetricError("metric configuration needs a non-empty enabled list")
    result: dict[str, Any] = {"metric_schema_version": METRIC_SCHEMA_VERSION}
    for family in enabled:
        family_config = configuration.get(family, {})
        fields = family_config.get("fields", {})
        if family == "h1":
            thresholds = family_config.get("thresholds", {})
            required_thresholds = (
                "position_threshold_m",
                "orientation_threshold_rad",
                "joint_step_threshold_rad",
                "swivel_step_threshold_rad",
                "input_step_threshold_rad",
            )
            missing = [name for name in required_thresholds if name not in thresholds]
            if missing:
                raise MetricError(f"missing H1 thresholds: {missing}")
            pose_success = _field(
                samples, fields, "pose_success", "map_pose_success"
            )
            differential_valid = _field(
                samples,
                fields,
                "differential_valid",
                "map_differential_valid",
            )
            branch_smooth = _field(
                samples,
                fields,
                "branch_smooth",
                "map_branch_smooth",
                optional=True,
            )
            if branch_smooth is None:
                # Pre-v3 evidence used map_differential_valid as the explicit
                # branch-smooth proxy.
                branch_smooth = differential_valid
            differential_applicable = _field(
                samples,
                fields,
                "differential_applicable",
                "map_differential_applicable",
                optional=True,
            )
            if differential_applicable is None:
                differential_applicable = np.zeros_like(
                    _boolean_vector(
                        differential_valid, "map_differential_valid"
                    ),
                    dtype=bool,
                )
            branch_array = _boolean_vector(
                branch_smooth, "map_branch_smooth"
            )
            applicable_array = _boolean_vector(
                differential_applicable, "map_differential_applicable"
            )
            mapping_valid = (
                _boolean_vector(pose_success, "map_pose_success")
                & branch_array
            )
            result.update(
                compute_h1_composite(
                    mapping_valid=mapping_valid,
                    position_error_m=_field(
                        samples,
                        fields,
                        "position_error_m",
                        "map_position_error_m",
                    ),
                    orientation_error_rad=_field(
                        samples,
                        fields,
                        "orientation_error_rad",
                        "map_orientation_error_rad",
                    ),
                    q_slave=_field(samples, fields, "q_slave", "map_q_slave"),
                    swivel_angle_rad=_field(
                        samples,
                        fields,
                        "swivel_angle_rad",
                        "map_swivel_angle_rad",
                    ),
                    master_step_norm=_field(
                        samples,
                        fields,
                        "master_step_norm",
                        "map_master_step_norm",
                    ),
                    accepted=_field(
                        samples,
                        fields,
                        "accepted",
                        "map_accepted",
                        optional=True,
                    ),
                    commanded_reset=_field(
                        samples,
                        fields,
                        "commanded_reset",
                        "map_commanded_reset",
                        optional=True,
                    ),
                    degeneracy_transition=_field(
                        samples,
                        fields,
                        "degeneracy_transition",
                        "map_degeneracy_transition",
                        optional=True,
                    ),
                    **thresholds,
                )
            )
            audit_fields = family_config.get("audit_fields", [])
            if audit_fields:
                required_audit_fields = {
                    "map_validity_reason_code",
                    "map_reach_clip_code",
                    "map_joint_limit_active",
                    "map_geometry_degenerate",
                    "map_slave_min_singular_value",
                    "map_slave_manipulability",
                    "map_low_manipulability",
                }
                if (
                    not isinstance(audit_fields, list)
                    or any(
                        not isinstance(name, str) or not name
                        for name in audit_fields
                    )
                ):
                    raise MetricError("H1 audit_fields must be a list of names")
                configured_audit_fields = set(audit_fields)
                if configured_audit_fields != required_audit_fields:
                    missing_audit = sorted(
                        required_audit_fields - configured_audit_fields
                    )
                    unknown_audit = sorted(
                        configured_audit_fields - required_audit_fields
                    )
                    raise MetricError(
                        "H1 audit_fields must declare the complete audit set; "
                        f"missing={missing_audit}, unknown={unknown_audit}"
                    )
                for stored_name in sorted(required_audit_fields):
                    if stored_name not in samples:
                        raise MetricError(
                            f"missing configured H1 audit field {stored_name!r}"
                        )
                result.update(
                    compute_h1_audit_metrics(
                        # H1 C_r is the common pose/branch endpoint. Actual
                        # differential validity is reported separately because
                        # the baseline methods do not yet expose an equivalent
                        # 7x7 differential implementation.
                        differential_valid=branch_array,
                        validity_reason_code=samples[
                            "map_validity_reason_code"
                        ],
                        reach_clip_code=samples["map_reach_clip_code"],
                        joint_limit_active=samples[
                            "map_joint_limit_active"
                        ],
                        geometry_degenerate=samples[
                            "map_geometry_degenerate"
                        ],
                        low_manipulability=samples[
                            "map_low_manipulability"
                        ],
                        slave_min_singular_value=samples[
                            "map_slave_min_singular_value"
                        ],
                        slave_manipulability=samples[
                            "map_slave_manipulability"
                        ],
                        validity_reason_labels=family_config.get(
                            "validity_reason_labels"
                        ),
                    )
                )
            branch_timing = family_config.get("branch_timing")
            if branch_timing is not None:
                if not isinstance(branch_timing, Mapping):
                    raise MetricError(
                        "H1 branch_timing configuration must be a mapping"
                    )
                strict_fields = {
                    "branch_smooth": "map_branch_smooth",
                    "differential_applicable":
                        "map_differential_applicable",
                    "differential_A": "map_differential_A",
                    "differential_runtime_s":
                        "map_differential_runtime_s",
                    "differential_max_one_sided_consistency": (
                        "map_differential_max_one_sided_consistency"
                    ),
                    "pose_runtime_s": "map_runtime_s",
                    "warm_start": "map_warm_start",
                    "phi_rad": "map_sew_phi_rad",
                }
                missing_strict = [
                    fields.get(logical_name, default_name)
                    for logical_name, default_name in strict_fields.items()
                    if fields.get(logical_name, default_name) not in samples
                ]
                if missing_strict:
                    raise MetricError(
                        "H1 branch_timing requires v3 evidence fields: "
                        f"{sorted(missing_strict)}"
                    )
                differential_A = np.asarray(
                    _field(
                        samples,
                        fields,
                        "differential_A",
                        "map_differential_A",
                    ),
                    dtype=float,
                )
                if differential_A.shape != (
                    applicable_array.shape[0],
                    7,
                    7,
                ):
                    raise MetricError(
                        "map_differential_A must have shape (samples, 7, 7)"
                    )
                differential_valid_array = _boolean_vector(
                    differential_valid, "map_differential_valid"
                )
                valid_differential_mask = (
                    applicable_array & differential_valid_array
                )
                if np.any(valid_differential_mask) and not np.all(
                    np.isfinite(
                        differential_A[valid_differential_mask]
                    )
                ):
                    raise MetricError(
                        "valid H1 differential matrices must be finite"
                    )
                result.update(
                    compute_h1_timing_branch_metrics(
                        pose_runtime_s=_field(
                            samples,
                            fields,
                            "pose_runtime_s",
                            "map_runtime_s",
                        ),
                        warm_start=_field(
                            samples,
                            fields,
                            "warm_start",
                            "map_warm_start",
                        ),
                        phi_rad=_field(
                            samples,
                            fields,
                            "phi_rad",
                            "map_sew_phi_rad",
                        ),
                        q_slave=_field(
                            samples, fields, "q_slave", "map_q_slave"
                        ),
                        branch_smooth=branch_smooth,
                        differential_applicable=(
                            differential_applicable
                        ),
                        differential_valid=differential_valid,
                        differential_runtime_s=_field(
                            samples,
                            fields,
                            "differential_runtime_s",
                            "map_differential_runtime_s",
                        ),
                        differential_max_one_sided_consistency=_field(
                            samples,
                            fields,
                            "differential_max_one_sided_consistency",
                            (
                                "map_differential_max_"
                                "one_sided_consistency"
                            ),
                        ),
                        deadline_s=branch_timing.get(
                            "deadline_s", 0.020
                        ),
                        minimum_phi_wrap_crossings=branch_timing.get(
                            "minimum_phi_wrap_crossings", 0
                        ),
                    )
                )
        elif family == "h2":
            result.update(
                compute_h2_wrench_metrics(
                    _field(
                        samples,
                        fields,
                        "wrench_estimated",
                        "wrench_estimated",
                    ),
                    _field(
                        samples,
                        fields,
                        "wrench_reference",
                        "wrench_reference",
                    ),
                    sample_mask=_field(
                        samples,
                        fields,
                        "sample_mask",
                        "wrench_sample_mask",
                        optional=True,
                    ),
                    numerical_rank_deficient=_field(
                        samples,
                        fields,
                        "numerical_rank_deficient",
                        "solver_numerical_rank_deficient",
                        optional=True,
                    ),
                    operationally_ill_conditioned=_field(
                        samples,
                        fields,
                        "operationally_ill_conditioned",
                        "solver_operationally_ill_conditioned",
                        optional=True,
                    ),
                    numerical_rank_threshold=_field(
                        samples,
                        fields,
                        "numerical_rank_threshold",
                        "solver_numerical_rank_threshold",
                        optional=True,
                    ),
                    operational_min_scaled_singular_threshold=_field(
                        samples,
                        fields,
                        "operational_min_scaled_singular_threshold",
                        "operational_min_scaled_singular_threshold",
                        optional=True,
                    ),
                    scaled_singular_values=_field(
                        samples,
                        fields,
                        "scaled_singular_values",
                        "scaled_singular_values",
                        optional=True,
                    ),
                    condition_number=_field(
                        samples,
                        fields,
                        "condition_number",
                        "solver_condition_number",
                        optional=True,
                    ),
                )
            )
        elif family == "h3":
            h3_eligible = _field(
                samples,
                fields,
                "h3_eligible",
                "h3_eligible",
                optional=True,
            )
            if h3_eligible is not None and not np.all(
                _boolean_vector(h3_eligible, "h3_eligible")
            ):
                raise MetricError(
                    "H3 is ineligible for samples containing an unpaired "
                    "synthetic energy probe"
                )
            probe_work = _field(
                samples,
                fields,
                "energy_probe_raw_work_J",
                "energy_probe_raw_work_J",
                optional=True,
            )
            if probe_work is not None:
                probe_work_array = _vector(
                    probe_work, "energy_probe_raw_work_J"
                )
                _finite(
                    probe_work_array, "energy_probe_raw_work_J"
                )
                if float(probe_work_array[-1]) > 1e-15:
                    raise MetricError(
                        "H3 is ineligible when synthetic energy-probe work "
                        "is nonzero"
                    )
            result.update(
                compute_h3_power_mismatch(
                    tau_master_raw=_field(
                        samples, fields, "tau_master_raw", "tau_master_raw"
                    ),
                    qd_master=_field(samples, fields, "qd_master", "qd_master"),
                    tau_slave_source=_field(
                        samples,
                        fields,
                        "tau_slave_source",
                        "tau_slave_source",
                    ),
                    qd_slave_source=_field(
                        samples,
                        fields,
                        "qd_slave_source",
                        "qd_slave_source",
                    ),
                    dt=_field(samples, fields, "dt", "dt"),
                    return_valid=_field(
                        samples,
                        fields,
                        "return_valid",
                        "return_valid",
                        optional=True,
                    ),
                    force_scale=family_config.get("force_scale", 1.0),
                    epsilon_energy_J=family_config.get(
                        "epsilon_energy_J", 1e-12
                    ),
                    minimum_power_activity_J=family_config.get(
                        "minimum_power_activity_J", 0.0
                    ),
                )
            )
        elif family == "h4":
            energy_min_J = _energy_bound_for_trial(
                samples,
                fields,
                family_config,
                "energy_min_J",
                "configured_energy_min_J",
            )
            energy_max_J = _energy_bound_for_trial(
                samples,
                fields,
                family_config,
                "energy_max_J",
                "configured_energy_max_J",
            )
            result.update(
                audit_h4_energy(
                    energy_before_J=_field(
                        samples, fields, "energy_before_J", "energy_before_J"
                    ),
                    energy_after_J=_field(
                        samples, fields, "energy_after_J", "energy_after_J"
                    ),
                    tau_candidate=_field(
                        samples, fields, "tau_candidate", "tau_master_candidate"
                    ),
                    tau_applied=_field(
                        samples, fields, "tau_applied", "tau_master_applied"
                    ),
                    tau_accepted=_field(
                        samples,
                        fields,
                        "tau_accepted",
                        "tau_master_accepted",
                        optional=True,
                    ),
                    qd_master=_field(samples, fields, "qd_master", "qd_master"),
                    dt=_field(samples, fields, "dt", "dt"),
                    software_preclip_J=_field(
                        samples,
                        fields,
                        "software_preclip_J",
                        "energy_preclip_J",
                        optional=True,
                    ),
                    energy_min_J=energy_min_J,
                    energy_max_J=energy_max_J,
                    epsilon_torque_impulse_Nms=family_config.get(
                        "epsilon_torque_impulse_Nms", 1e-12
                    ),
                    audit_tolerance_J=family_config.get(
                        "audit_tolerance_J", 1e-10
                    ),
                )
            )
        elif family == "bilateral":
            required_audit_fields = family_config.get(
                "required_audit_fields", []
            )
            if (
                not isinstance(required_audit_fields, list)
                or any(
                    not isinstance(name, str) or not name
                    for name in required_audit_fields
                )
            ):
                raise MetricError(
                    "bilateral required_audit_fields must be a list of names"
                )
            missing_audit_fields = sorted(
                name
                for name in required_audit_fields
                if name not in samples
            )
            if missing_audit_fields:
                raise MetricError(
                    "missing required bilateral audit fields: "
                    f"{missing_audit_fields}"
                )
            stored_force_limit = _field(
                samples,
                fields,
                "wall_force_limit_N",
                "configured_wall_force_limit_N",
                optional=True,
            )
            gate_config = family_config.get("gates", {})
            if not isinstance(gate_config, Mapping):
                raise MetricError("bilateral gates must be a mapping")
            result.update(
                compute_bilateral_diagnostics(
                    master_tracking_error_rad=_field(
                        samples,
                        fields,
                        "master_tracking_error_rad",
                        "master_tracking_error",
                    ),
                    slave_tracking_error_rad=_field(
                        samples,
                        fields,
                        "slave_tracking_error_rad",
                        "slave_tracking_error",
                    ),
                    feedback_torque_Nm=_field(
                        samples,
                        fields,
                        "feedback_torque_Nm",
                        "tau_master_applied",
                    ),
                    contact_force_N=_field(
                        samples,
                        fields,
                        "contact_force_N",
                        "contact_force_norm",
                    ),
                    projection_factor=_field(
                        samples,
                        fields,
                        "projection_factor",
                        "rho",
                    ),
                    wall_force_raw_N=_field(
                        samples,
                        fields,
                        "wall_force_raw_N",
                        "wall_force_raw_N",
                        optional=True,
                    ),
                    wall_force_applied_N=_field(
                        samples,
                        fields,
                        "wall_force_applied_N",
                        "wall_force_applied_N",
                        optional=True,
                    ),
                    wall_force_saturation_active=_field(
                        samples,
                        fields,
                        "wall_force_saturation_active",
                        "wall_force_saturation_active",
                        optional=True,
                    ),
                    wall_force_limit_N=(
                        stored_force_limit
                        if stored_force_limit is not None
                        else family_config.get("wall_force_limit_N")
                    ),
                    master_joint_limit_active=_field(
                        samples,
                        fields,
                        "master_joint_limit_active",
                        "master_joint_limit_active",
                        optional=True,
                    ),
                    slave_joint_limit_active=_field(
                        samples,
                        fields,
                        "slave_joint_limit_active",
                        "slave_joint_limit_active",
                        optional=True,
                    ),
                    master_velocity_limit_active=_field(
                        samples,
                        fields,
                        "master_velocity_limit_active",
                        "master_velocity_limit_active",
                        optional=True,
                    ),
                    slave_velocity_limit_active=_field(
                        samples,
                        fields,
                        "slave_velocity_limit_active",
                        "slave_velocity_limit_active",
                        optional=True,
                    ),
                    master_acceleration_limit_active=_field(
                        samples,
                        fields,
                        "master_acceleration_limit_active",
                        "master_acceleration_limit_active",
                        optional=True,
                    ),
                    slave_acceleration_limit_active=_field(
                        samples,
                        fields,
                        "slave_acceleration_limit_active",
                        "slave_acceleration_limit_active",
                        optional=True,
                    ),
                    master_torque_saturation_active=_field(
                        samples,
                        fields,
                        "master_torque_saturation_active",
                        "master_torque_saturation_active",
                        optional=True,
                    ),
                    slave_torque_saturation_active=_field(
                        samples,
                        fields,
                        "slave_torque_saturation_active",
                        "slave_torque_saturation_active",
                        optional=True,
                    ),
                    haptic_rate_limit_active=_field(
                        samples,
                        fields,
                        "haptic_rate_limit_active",
                        "haptic_rate_limit_active",
                        optional=True,
                    ),
                    haptic_torque_saturation_active=_field(
                        samples,
                        fields,
                        "haptic_torque_saturation_active",
                        "haptic_torque_saturation_active",
                        optional=True,
                    ),
                    energy_probe_raw_work_J=_field(
                        samples,
                        fields,
                        "energy_probe_raw_work_J",
                        "energy_probe_raw_work_J",
                        optional=True,
                    ),
                    projection_tolerance=family_config.get(
                        "projection_tolerance", 1e-12
                    ),
                    contact_force_threshold_N=family_config.get(
                        "contact_force_threshold_N", 1e-6
                    ),
                    minimum_contact_fraction=gate_config.get(
                        "minimum_contact_fraction", 0.0
                    ),
                    minimum_contact_rms_N=gate_config.get(
                        "minimum_contact_rms_N", 0.0
                    ),
                    maximum_force_limit_hit_fraction=gate_config.get(
                        "maximum_force_limit_hit_fraction", 1.0
                    ),
                    minimum_force_headroom_N=gate_config.get(
                        "minimum_force_headroom_N", 0.0
                    ),
                    maximum_master_tracking_rmse_rad=gate_config.get(
                        "maximum_master_tracking_rmse_rad"
                    ),
                    maximum_slave_tracking_rmse_rad=gate_config.get(
                        "maximum_slave_tracking_rmse_rad"
                    ),
                    minimum_projection_intervention_fraction=(
                        gate_config.get(
                            "minimum_projection_intervention_fraction",
                            0.0,
                        )
                    ),
                    maximum_projection_intervention_fraction=(
                        gate_config.get(
                            "maximum_projection_intervention_fraction",
                            1.0,
                        )
                    ),
                    maximum_limit_active_fraction=gate_config.get(
                        "maximum_limit_active_fraction", 1.0
                    ),
                    minimum_energy_probe_raw_work_J=gate_config.get(
                        "minimum_energy_probe_raw_work_J", 0.0
                    ),
                )
            )
        elif family == "network":
            required_upstream_metrics = (
                "h4_energy_audit_pass",
                "bilateral_force_limit_gate_pass",
                "bilateral_force_headroom_gate_pass",
                "bilateral_master_tracking_gate_pass",
                "bilateral_slave_tracking_gate_pass",
                "bilateral_projection_gate_pass",
                "bilateral_limit_active_gate_pass",
            )
            missing_upstream_metrics = [
                name
                for name in required_upstream_metrics
                if name not in result
            ]
            if missing_upstream_metrics:
                raise MetricError(
                    "network metrics must follow H4 and bilateral metrics; "
                    f"missing={missing_upstream_metrics}"
                )
            gate_config = family_config.get("gates", {})
            if not isinstance(gate_config, Mapping):
                raise MetricError("network gates must be a mapping")
            network_metrics = compute_network_diagnostics(
                forward_packet_state=_field(
                    samples,
                    fields,
                    "forward_packet_state",
                    "forward_packet_state",
                ),
                return_packet_state=_field(
                    samples,
                    fields,
                    "return_packet_state",
                    "return_packet_state",
                ),
                forward_packet_active=_field(
                    samples,
                    fields,
                    "forward_packet_active",
                    "forward_packet_active",
                ),
                return_packet_active=_field(
                    samples,
                    fields,
                    "return_packet_active",
                    "return_packet_active",
                ),
                forward_packet_fresh=_field(
                    samples,
                    fields,
                    "forward_packet_fresh",
                    "forward_packet_fresh",
                ),
                return_packet_fresh=_field(
                    samples,
                    fields,
                    "return_packet_fresh",
                    "return_packet_fresh",
                ),
                forward_packet_age=_field(
                    samples,
                    fields,
                    "forward_packet_age",
                    "forward_packet_age",
                ),
                return_packet_age=_field(
                    samples,
                    fields,
                    "return_packet_age",
                    "return_packet_age",
                ),
                forward_packet_seq=_field(
                    samples,
                    fields,
                    "forward_packet_seq",
                    "forward_packet_seq",
                ),
                return_packet_seq=_field(
                    samples,
                    fields,
                    "return_packet_seq",
                    "return_packet_seq",
                ),
                slave_zero_delay_tracking_error=_field(
                    samples,
                    fields,
                    "slave_zero_delay_tracking_error",
                    "slave_zero_delay_tracking_error",
                ),
                slave_reference_lag_error=_field(
                    samples,
                    fields,
                    "slave_reference_lag_error",
                    "slave_reference_lag_error",
                ),
                return_feedback_lag_error=_field(
                    samples,
                    fields,
                    "return_feedback_lag_error",
                    "return_feedback_lag_error",
                ),
                contact_force_norm=_field(
                    samples,
                    fields,
                    "contact_force_norm",
                    "contact_force_norm",
                ),
                contact_expected=_field(
                    samples,
                    fields,
                    "contact_expected",
                    "contact_expected",
                ),
                configured_forward_delay_s=_field(
                    samples,
                    fields,
                    "configured_forward_delay_s",
                    "configured_forward_delay_s",
                ),
                configured_return_delay_s=_field(
                    samples,
                    fields,
                    "configured_return_delay_s",
                    "configured_return_delay_s",
                ),
                configured_forward_jitter_s=_field(
                    samples,
                    fields,
                    "configured_forward_jitter_s",
                    "configured_forward_jitter_s",
                ),
                configured_return_jitter_s=_field(
                    samples,
                    fields,
                    "configured_return_jitter_s",
                    "configured_return_jitter_s",
                ),
                configured_forward_packet_loss=_field(
                    samples,
                    fields,
                    "configured_forward_packet_loss",
                    "configured_forward_packet_loss",
                ),
                configured_return_packet_loss=_field(
                    samples,
                    fields,
                    "configured_return_packet_loss",
                    "configured_return_packet_loss",
                ),
                configured_forward_timeout_s=_field(
                    samples,
                    fields,
                    "configured_forward_timeout_s",
                    "configured_forward_timeout_s",
                ),
                configured_return_timeout_s=_field(
                    samples,
                    fields,
                    "configured_return_timeout_s",
                    "configured_return_timeout_s",
                ),
                contact_force_threshold_N=family_config.get(
                    "contact_force_threshold_N", 1e-6
                ),
                maximum_zero_delay_tracking_rmse_rad=gate_config.get(
                    "maximum_slave_zero_delay_tracking_rmse_rad",
                    gate_config.get(
                        "maximum_zero_delay_tracking_rmse_rad"
                    ),
                ),
                minimum_active_fraction=gate_config.get(
                    "minimum_active_fraction", 0.0
                ),
                minimum_forward_active_fraction=gate_config.get(
                    "minimum_forward_active_fraction"
                ),
                minimum_return_active_fraction=gate_config.get(
                    "minimum_return_active_fraction"
                ),
                minimum_forward_fresh_fraction=gate_config.get(
                    "minimum_forward_fresh_fraction", 0.0
                ),
                minimum_return_fresh_fraction=gate_config.get(
                    "minimum_return_fresh_fraction", 0.0
                ),
                maximum_forward_internal_missing_fraction=gate_config.get(
                    "maximum_forward_internal_missing_fraction", 1.0
                ),
                maximum_return_internal_missing_fraction=gate_config.get(
                    "maximum_return_internal_missing_fraction", 1.0
                ),
                maximum_timeout_fraction=gate_config.get(
                    "maximum_timeout_fraction", 1.0
                ),
                maximum_forward_timeout_fraction=gate_config.get(
                    "maximum_forward_timeout_fraction"
                ),
                maximum_return_timeout_fraction=gate_config.get(
                    "maximum_return_timeout_fraction"
                ),
                maximum_age_s=gate_config.get("maximum_age_s"),
                maximum_forward_age_s=gate_config.get(
                    "maximum_forward_packet_age_s",
                    gate_config.get("maximum_forward_age_s"),
                ),
                maximum_return_age_s=gate_config.get(
                    "maximum_return_packet_age_s",
                    gate_config.get("maximum_return_age_s"),
                ),
                minimum_contact_fraction=gate_config.get(
                    "minimum_contact_fraction", 0.0
                ),
                minimum_contact_rms_N=gate_config.get(
                    "minimum_contact_rms_N", 0.0
                ),
                maximum_free_space_contact_peak_N=gate_config.get(
                    "maximum_free_space_contact_force_N",
                    gate_config.get(
                        "maximum_free_space_contact_peak_N"
                    ),
                ),
            )
            upstream_gate_pass = all(
                bool(result[name]) for name in required_upstream_metrics
            )
            network_metrics.update(
                {
                    "network_h4_bilateral_gate_pass": (
                        upstream_gate_pass
                    ),
                    "network_local_full_gate_pass": bool(
                        upstream_gate_pass
                        and network_metrics[
                            "network_local_gate_pass"
                        ]
                    ),
                }
            )
            result.update(network_metrics)
        elif family == "energy_challenge":
            required_metrics = (
                "h4_energy_audit_pass",
                "h4_shadow_floor_deficit_J",
                "h4_downstream_modification_max_Nm",
                "h4_D_proj",
                "bilateral_stable_contact_gate_pass",
            )
            missing_metrics = [
                name for name in required_metrics if name not in result
            ]
            if missing_metrics:
                raise MetricError(
                    "energy_challenge must follow H4 and bilateral metrics; "
                    f"missing={missing_metrics}"
                )
            minimum_shadow_deficit_J = float(
                family_config.get("minimum_shadow_deficit_J", 0.0)
            )
            maximum_downstream_modification_Nm = float(
                family_config.get(
                    "maximum_downstream_modification_Nm", 1e-10
                )
            )
            maximum_D_proj = float(
                family_config.get("maximum_D_proj", 1.0)
            )
            challenge_thresholds = (
                minimum_shadow_deficit_J,
                maximum_downstream_modification_Nm,
                maximum_D_proj,
            )
            if any(
                not np.isfinite(value) or value < 0.0
                for value in challenge_thresholds
            ):
                raise MetricError(
                    "energy challenge thresholds must be finite/non-negative"
                )
            shadow_gate = bool(
                result["h4_shadow_floor_deficit_J"]
                >= minimum_shadow_deficit_J
            )
            downstream_gate = bool(
                result["h4_downstream_modification_max_Nm"]
                <= maximum_downstream_modification_Nm
            )
            distortion_gate = bool(
                result["h4_D_proj"] <= maximum_D_proj
            )
            result.update(
                {
                    "energy_challenge_shadow_gate_pass": shadow_gate,
                    "energy_challenge_downstream_gate_pass":
                        downstream_gate,
                    "energy_challenge_distortion_gate_pass":
                        distortion_gate,
                    "energy_challenge_gate_pass": bool(
                        result["h4_energy_audit_pass"]
                        and result[
                            "bilateral_stable_contact_gate_pass"
                        ]
                        and shadow_gate
                        and downstream_gate
                        and distortion_gate
                    ),
                }
            )
        else:
            raise MetricError(f"unknown metric family {family!r}")
    return result