"""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