From 05f53e22274cdfe7e8b705247caaf481f9a8c108 Mon Sep 17 00:00:00 2001 From: James Sweeney Date: Tue, 15 Sep 2026 18:07:59 -0700 Subject: [PATCH 01/10] Add synthetic damped-mode and Cramer-Rao validation Adds cuphoton.xray.synthetic_validation: a fixture of known damped modes (the synthetic_trace modes plus white Gaussian noise), Cramer-Rao bounds from the exact Fisher information of the same model, a Monte Carlo sweep that reports per-mode bias, standard deviation and rmse of every parameter against the bound together with the mode-loss rate, optional model-mismatch traces with a residual-to-noise ratio, and a schema_version 1 summary (config, runtime via runtime_metadata, results per condition with attempted/successful/failed trials, artifacts). The linear-prediction-validate (lpv) command runs the sweep and, with --output-dir, writes summary.json and a Bokeh figure outside the checkout. Documented in docs/xray/LINEAR-PREDICTION-VALIDATION.md. The sweep is the reproducer for the mode loss discussed in #1: a lightly damped mode whose fitted decay crosses zero under noise is dropped by the root filter. No default behaviour changes. Refs #1 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_014wyj5iPSHmtR3vSTvJWKfw Signed-off-by: James Sweeney Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-VALIDATION.md | 131 ++++ docs/xray/README.md | 2 + src/cuphoton/xray/commands.py | 142 ++++ src/cuphoton/xray/synthetic_validation.py | 747 ++++++++++++++++++++++ tests/core/test_cli_contract.py | 10 +- tests/xray/test_cli.py | 51 ++ tests/xray/test_synthetic_validation.py | 207 ++++++ 7 files changed, 1285 insertions(+), 5 deletions(-) create mode 100644 docs/xray/LINEAR-PREDICTION-VALIDATION.md create mode 100644 src/cuphoton/xray/synthetic_validation.py create mode 100644 tests/xray/test_synthetic_validation.py diff --git a/docs/xray/LINEAR-PREDICTION-VALIDATION.md b/docs/xray/LINEAR-PREDICTION-VALIDATION.md new file mode 100644 index 00000000..510630e7 --- /dev/null +++ b/docs/xray/LINEAR-PREDICTION-VALIDATION.md @@ -0,0 +1,131 @@ +# Linear prediction validation on synthetic damped modes + +`linear-prediction-validate` (`lpv`) checks the linear-prediction fit against +a known answer. It builds traces from a small set of damped sinusoids, adds +white Gaussian noise at chosen signal-to-noise levels, runs the fit many +times, and compares the scatter and bias of every recovered parameter with +the Cramer-Rao lower bound for the same model. It also counts how often a +true mode is missing from the fit. No external data is needed. + +```bash +uv run cuphoton xray lpv --snr-db 40,30,20,10 --trials 200 --output-dir /tmp/lpv +uv run cuphoton xray lpv --distortion chirp:0.05 --snr-db 30 --output-dir /tmp/lpv-chirp +uv run cuphoton xray lpv --json +``` + +Keep `--output-dir` outside the checkout. It receives `summary.json` and, when +the `viz` extra is installed, `validation.html`. + +## Model and fixture + +The trace is + +``` +trace(t) = constant + sum_k amplitude_k exp(-decay_k t) cos(angular_frequency_k t + phase_k) +``` + +sampled at `samples` points from `t = 0` to `duration`. Time is in the input +time unit; angular frequency is in radians per time unit and decay in inverse +time units. The default fixture has two modes, `(1.25, 0.09, 2.4, 0.30)` and +`(0.45, 0.03, 0.9, -0.75)` as `(amplitude, decay, angular_frequency, phase)`, +constant `0.15`, 96 samples over 9.5 time units. With zero noise it is the +same trace as `synthetic_trace`, and the fit recovers the modes to better than +`1e-6`. `cuphoton.xray.synthetic_validation.synthetic_modes_trace` builds the +fixture; `DampedMode` describes a mode. + +Noise is white Gaussian, independent per sample. A level's `snr_db` is +`20 log10(rms of the mean-removed clean trace / sigma)`. One `numpy` generator +seeded by `--seed` draws every trial in order, so a run is reproducible from +the summary's `config`. + +## Bounds + +`cramer_rao_bounds` forms the Fisher information `J^T J / sigma^2` from the +numerical Jacobian of the model at the true parameters (all parameters free, +including the constant) and returns the square root of the diagonal of its +inverse per parameter. These are standard-deviation bounds; the summary +also carries the squared value as `crlb_variance`. For a single undamped +sinusoid the angular-frequency bound matches the closed form +`24 sigma^2 / (A^2 dt^2 N (N^2 - 1))` (Kay, Estimation Theory, 1993) to +within 3 percent, the difference being the freed decay and constant. + +## Mode matching and statistics + +A true mode is recovered in a trial when a fitted mode has an angular +frequency within `match_tolerance` (default 0.3 rad per time unit) of it; +fitted modes are assigned nearest-first and each is used at most once. The +per-mode `loss_rate` is the fraction of trials in which the mode was not +recovered. A trial is successful when the estimator returned and every true +mode was recovered; a trial in which the estimator raised counts in +`estimator_errors` and as a loss of every mode. + +For each mode and parameter (`amplitude`, `decay`, `angular_frequency`, +`phase`) the summary gives `bias`, `std` (ddof 1) and `rmse` computed over +the recovered trials only, the bound (`crlb_std`, `crlb_variance`) and +`std_over_crlb_std`. A negative fitted amplitude is folded into the phase +and the phase error is wrapped to `[-pi, pi)`. + +Read the decay statistics of a lightly damped mode together with its +`loss_rate`. The root filter keeps decaying roots only, so when noise pushes +the fitted decay of such a mode through zero the mode is dropped rather than +reported with a negative decay; the surviving decay estimates are truncated +at zero, and their scatter and bias are conditional on survival. The summary +labels these statistics `statistics_over: recovered_trials`. + +Every level also reports `residual_rms_over_sigma_median`, the median over +trials of the rms residual of the reconstruction divided by sigma. It is +close to 1 when the model class fits and well above 1 when it does not. + +## Model mismatch + +`--distortion kind:amount` replaces the clean trace with one outside the +model class: + +| kind | trace | +| --- | --- | +| `chirp` | every frequency rises by the fraction `amount` across the record | +| `gaussian_envelope` | the first mode decays as a Gaussian with 1/e time `amount` | +| `baseline_drift` | a linear ramp of `amount` across the record | +| `clip` | the trace is clipped above `amount` | +| `glitch` | one sample at mid-record is offset by `amount` | + +The truth used for bias and bounds stays the undistorted mode set, so the +statistics measure the estimator's response to the mismatch. A chirp or a +clipped trace can return modes with a confident, biased frequency and no +loss; the residual ratio is the number that says so. + +## Summary schema + +`summary.json` has `schema_version` 1 and five top-level keys. + +- `command`: `linear-prediction-validate`. +- `config`: the model string and units, the true modes and constant, the + sampling (`samples`, `duration`, `dt`, `start`), the noise definition and + levels, `seed`, `trials_per_level`, `n_components` (model order), the + mode-matching criterion and tolerance, how the statistics are computed, and + the estimators and distortions covered. +- `runtime`: `cuphoton.core.runtime.runtime_metadata` for the CPU backend + (package, Python, platform, NumPy versions, backend, device, dtype) plus + `scipy_version` and `source_revision` (the git revision of the checkout, + or `null`). +- `results`: one entry per sweep condition (estimator by distortion by + signal-to-noise level): `snr_db`, `noise_sigma`, `trials` with + `attempted`, `successful`, `failed` and `estimator_errors`, + `any_mode_lost_rate`, `residual_rms_over_sigma_median`, and `modes`, each + with `recovered_trials`, `loss_rate`, `statistics_over` and the four + parameter blocks described above. +- `artifacts`: filenames relative to the summary, `summary` and `figure` + (`null` when Bokeh is not installed). + +The figure shows, per estimator and per parameter, the sample standard +deviation against the bound versus signal-to-noise, and the loss rates. + +## Python entry points + +`cuphoton.xray.synthetic_validation` exposes `synthetic_modes_trace`, +`distort_trace`, `cramer_rao_bounds`, `match_modes`, `validation_sweep`, +`build_summary`, `write_validation_figure` and `write_validation_run`. +`validation_sweep` accepts another estimator callable with the signature +`estimator(time, trace, n_components)` returning an object with +`angular_frequency`, `decay`, `amplitude`, `phase` and `reconstruction` +(NumPy or CuPy arrays), so the same sweep can compare estimators. diff --git a/docs/xray/README.md b/docs/xray/README.md index 62b4b167..04ad1a4a 100644 --- a/docs/xray/README.md +++ b/docs/xray/README.md @@ -34,6 +34,7 @@ Choose a new output directory for each quickstart run. | Input inspection | `data-probe`, `roi-candidates` | | Trace work | `extract-trace`, `trace-smoke` | | Linear prediction | `linear-prediction-*`, `prediction-roots-benchmark`, `model-order-sweep` | +| Synthetic validation | `linear-prediction-validate` (`lpv`) | | Detector products | `detector-mask`, `detector-artifacts`, `detector-artifact-normalize`, `detector-artifact-compare` | | Distributed detector work | `detector-artifact-distributed`, `detector-artifact-merge` | | Review | `report`, `validation-viz`, `workflow-viz`, `phonon-viz` | @@ -295,3 +296,4 @@ See [Data and artifact contracts](../data-artifacts.md#xray-hdf5-and-trace-produ - [GPU-first behavior](GPU-FIRST.md) - [Distributed detector artifacts](DISTRIBUTED-DETECTOR-ARTIFACTS.md) - [Validation visualization](VALIDATION-VIZ.md) +- [Linear prediction validation](LINEAR-PREDICTION-VALIDATION.md) diff --git a/src/cuphoton/xray/commands.py b/src/cuphoton/xray/commands.py index d1c6861e..6929b47e 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -596,6 +596,84 @@ class JsonArg(BoolInvariant): _required = False +class LinearPredictionValidateCommand(_XRayCommand): + _description_ = ( + "Validate linear prediction on synthetic damped modes against " + "the Cramer-Rao bound." + ) + _shortname_ = "lpv" + _handler_name_ = "_linear_prediction_validate" + + samples = 96 + + class SamplesArg(IntegerInvariant): + _arg = "--samples" + _help = "Number of synthetic trace samples." + _required = False + _default = 96 + + components = 6 + + class ComponentsArg(IntegerInvariant): + _arg = "--components" + _help = "Number of SVD components to fit." + _required = False + _default = 6 + + trials = 200 + + class TrialsArg(IntegerInvariant): + _arg = "--trials" + _help = "Monte Carlo trials per signal-to-noise level." + _required = False + _default = 200 + + snr_db = "40,30,20,10" + + class SnrDbArg(StringInvariant): + _arg = "--snr-db" + _help = "Comma-separated signal-to-noise levels in dB." + _required = False + _default = "40,30,20,10" + + seed = 20260914 + + class SeedArg(IntegerInvariant): + _arg = "--seed" + _help = "Random seed for the noise draws." + _required = False + _default = 20260914 + + distortion = None + + class DistortionArg(StringInvariant): + _arg = "--distortion" + _help = ( + "Optional model mismatch as kind:amount, one of chirp, " + "gaussian_envelope, baseline_drift, clip, glitch " + "(for example chirp:0.05)." + ) + _required = False + _default = None + + output_dir = None + + class OutputDirArg(PathValueInvariant): + _arg = "--output-dir" + _help = ( + "Directory for summary.json and the validation figure " + "(created if needed; keep it outside the checkout)." + ) + _required = False + + json = False + + class JsonArg(BoolInvariant): + _arg = "--json" + _help = "Print the summary as JSON instead of the text report." + _required = False + + class LinearPredictionBenchmarkCommand(_XRayCommand): _description_ = "Benchmark serial and batched linear-prediction P1." _shortname_ = "lpb" @@ -4342,6 +4420,70 @@ def _linear_prediction_smoke(args): return 0 +def _print_validation_sweep(sweep): + print(f"estimator={sweep.backend}") + if sweep.distortion is not None: + print(f"distortion={sweep.distortion[0]}:{sweep.distortion[1]:g}") + for level in sweep.levels: + print( + f"snr_db={level.snr_db:g} sigma={level.noise_sigma:.4g} " + f"trials={level.trials_successful}/{level.trials_attempted} " + f"any_mode_lost_rate={level.any_mode_lost_rate:.3f} " + f"residual_ratio={level.residual_ratio:.2f}" + ) + for m in level.modes: + w = m.angular_frequency + d = m.decay + print( + f" mode w={w.truth:g} decay={d.truth:g}: " + f"freq std/crlb={w.std:.3g}/{w.crlb_std:.3g} " + f"({w.std_over_crlb_std:.2f}x) bias={w.bias:+.3g}; " + f"decay std/crlb={d.std:.3g}/{d.crlb_std:.3g} " + f"({d.std_over_crlb_std:.2f}x) bias={d.bias:+.3g}; " + f"recovered={m.recovered_trials} " + f"loss_rate={m.loss_rate:.3f}" + ) + + +def _linear_prediction_validate(args): + from .synthetic_validation import ( + build_summary, + validation_sweep, + write_validation_run, + ) + + snr = tuple(float(x) for x in str(args.snr_db).split(",") if x.strip()) + distortion = None + if args.distortion: + kind, _, amount = str(args.distortion).partition(":") + distortion = (kind.strip(), float(amount)) + sweeps = [ + validation_sweep( + samples=args.samples, + snr_db=snr, + trials=args.trials, + n_components=args.components, + seed=args.seed, + distortion=distortion, + ) + ] + if args.output_dir is not None: + summary = write_validation_run(args.output_dir, sweeps) + else: + summary = build_summary(sweeps) + if args.json: + print(json.dumps(summary, indent=2, sort_keys=True)) + return 0 + print(f"samples={sweeps[0].samples}") + print(f"signal_rms={sweeps[0].signal_rms:.6g}") + for sweep in sweeps: + _print_validation_sweep(sweep) + if args.output_dir is not None: + for key, name in summary["artifacts"].items(): + print(f"{key}={name if name else 'not written'}") + return 0 + + def _linear_prediction_benchmark(args): from .linear_prediction import benchmark_linear_prediction_p1_batch diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py new file mode 100644 index 00000000..1db5a48c --- /dev/null +++ b/src/cuphoton/xray/synthetic_validation.py @@ -0,0 +1,747 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +"""Synthetic damped-mode fixtures and Cramer-Rao validation for the +linear-prediction fit. + +The fixture generates traces with known modes so the workflow can be +validated without external datasets. The validation sweep compares the +estimator's scatter and bias with the Cramer-Rao lower bound computed from +the exact Fisher information of the same model, together with the rate at +which true modes are lost from the fit. Everything here runs on NumPy; the +estimator under test is :func:`linear_prediction_numpy` unless another +estimator callable is supplied. +""" + +from __future__ import annotations + +import json +import subprocess +from dataclasses import asdict, dataclass, field +from pathlib import Path +from typing import Any, Callable, Sequence + +import numpy as np + +from cuphoton.core.runtime import runtime_metadata + +from .linear_prediction import linear_prediction_numpy + +SUMMARY_SCHEMA_VERSION = 1 +PARAMETERS = ("amplitude", "decay", "angular_frequency", "phase") + + +@dataclass(frozen=True) +class DampedMode: + """One real damped sinusoid. + + ``amplitude * exp(-decay * t) * cos(angular_frequency * t + phase)`` + """ + + amplitude: float + decay: float + angular_frequency: float + phase: float = 0.0 + + +@dataclass(frozen=True) +class SyntheticModesTrace: + time: np.ndarray + trace: np.ndarray + clean: np.ndarray + modes: tuple[DampedMode, ...] + constant: float + noise_sigma: float + seed: int | None + + +DEFAULT_MODES = ( + DampedMode(1.25, 0.09, 2.4, 0.30), + DampedMode(0.45, 0.03, 0.9, -0.75), +) + + +def modes_model( + theta: np.ndarray, time: np.ndarray, n_modes: int +) -> np.ndarray: + """Evaluate the multi-mode model from a flat parameter vector. + + Four parameters per mode, then the constant. + """ + out = np.full(time.shape, float(theta[4 * n_modes])) + for k in range(n_modes): + a, d, w, p = theta[4 * k : 4 * k + 4] + out = out + a * np.exp(-d * time) * np.cos(w * time + p) + return out + + +def modes_to_theta( + modes: Sequence[DampedMode], constant: float +) -> np.ndarray: + return np.array( + [ + v + for m in modes + for v in (m.amplitude, m.decay, m.angular_frequency, m.phase) + ] + + [constant] + ) + + +def synthetic_modes_trace( + samples: int = 96, + *, + modes: Sequence[DampedMode] = DEFAULT_MODES, + constant: float = 0.15, + duration: float = 9.5, + noise_sigma: float = 0.0, + seed: int | None = None, +) -> SyntheticModesTrace: + """Build a trace with known modes plus optional white Gaussian noise. + + With the defaults and ``noise_sigma=0`` this reproduces + :func:`synthetic_trace` exactly. + """ + if samples < 16: + raise ValueError("samples must be at least 16") + time = np.linspace(0.0, duration, samples, dtype=np.float64) + clean = modes_model(modes_to_theta(modes, constant), time, len(modes)) + trace = clean + if noise_sigma > 0: + rng = np.random.default_rng(seed) + trace = clean + rng.normal(0.0, noise_sigma, size=time.shape) + return SyntheticModesTrace( + time, trace, clean, tuple(modes), constant, float(noise_sigma), seed + ) + + +DISTORTIONS = ( + "chirp", + "gaussian_envelope", + "baseline_drift", + "clip", + "glitch", +) + + +def distort_trace( + time: np.ndarray, + modes: Sequence[DampedMode], + constant: float, + kind: str, + amount: float, +) -> np.ndarray: + """Return a clean trace that departs from the damped-mode model. + + ``chirp``: every frequency rises by the fraction ``amount`` across the + record. ``gaussian_envelope``: the first mode decays as a Gaussian with + 1/e time ``amount`` instead of an exponential. ``baseline_drift``: a + linear ramp of ``amount`` across the record. ``clip``: the trace is + clipped above ``amount``. ``glitch``: one sample at mid-record is + offset by ``amount``. These are outside the model class on purpose: + the sweep reports how the fit responds and the residual-to-noise ratio + that should flag them. + """ + if kind not in DISTORTIONS: + raise ValueError(f"unknown distortion {kind!r}") + span = float(time[-1] - time[0]) + out = np.full(time.shape, float(constant)) + for k, m in enumerate(modes): + w = m.angular_frequency + if kind == "chirp": + w = w * (1.0 + amount * (time - time[0]) / span) + if kind == "gaussian_envelope" and k == 0: + env = np.exp(-((time / amount) ** 2)) + else: + env = np.exp(-m.decay * time) + out = out + m.amplitude * env * np.cos(w * time + m.phase) + if kind == "baseline_drift": + out = out + amount * (time - time[0]) / span + elif kind == "clip": + out = np.minimum(out, amount) + elif kind == "glitch": + out = out.copy() + out[time.size // 2] += amount + return out + + +def cramer_rao_bounds( + modes: Sequence[DampedMode], + constant: float, + time: np.ndarray, + noise_sigma: float, + step: float = 1e-6, +) -> dict[str, np.ndarray]: + """Cramer-Rao lower bound (standard deviation) per parameter for white + Gaussian noise of ``noise_sigma``. + + The Fisher information is J^T J / sigma^2 with J the numerical Jacobian + of :func:`modes_model`; the bound is the square root of the diagonal of + its inverse. Returned per mode as arrays ``amplitude``, ``decay``, + ``angular_frequency``, ``phase`` plus the scalar ``constant``. + """ + if noise_sigma <= 0: + raise ValueError("noise_sigma must be positive") + theta = modes_to_theta(modes, constant) + n = len(modes) + jac = np.empty((time.size, theta.size)) + for k in range(theta.size): + tp = theta.copy() + tm = theta.copy() + tp[k] += step + tm[k] -= step + jac[:, k] = (modes_model(tp, time, n) - modes_model(tm, time, n)) / ( + 2 * step + ) + fisher = jac.T @ jac / noise_sigma**2 + sd = np.sqrt(np.diag(np.linalg.inv(fisher))) + return { + "amplitude": sd[0 : 4 * n : 4], + "decay": sd[1 : 4 * n : 4], + "angular_frequency": sd[2 : 4 * n : 4], + "phase": sd[3 : 4 * n : 4], + "constant": sd[4 * n], + } + + +def _to_host(values) -> np.ndarray: + """Return a float NumPy array; CuPy arrays are copied to the host.""" + get = getattr(values, "get", None) + if get is not None and not isinstance(values, np.ndarray): + values = get() + return np.asarray(values, dtype=float) + + +def match_modes( + modes: Sequence[DampedMode], + angular_frequency: np.ndarray, + decay: np.ndarray, + *, + tolerance: float, +) -> list[int | None]: + """For each true mode, the index of the fitted mode within ``tolerance`` + (rad per time unit) or ``None``. + """ + out: list[int | None] = [] + used: set[int] = set() + for m in modes: + if angular_frequency.size == 0: + out.append(None) + continue + err = np.abs(angular_frequency - m.angular_frequency) + for idx in np.argsort(err): + if idx in used: + continue + if err[idx] <= tolerance: + out.append(int(idx)) + used.add(int(idx)) + break + else: + out.append(None) + return out + + +@dataclass(frozen=True) +class ParameterStats: + """Statistics of one estimated parameter over the trials in which the + mode was recovered, against the unconditional Cramer-Rao bound. + + ``crlb_std`` is the standard-deviation bound and ``crlb_variance`` its + square; ``std_over_crlb_std`` compares the sample standard deviation + with the former. + """ + + truth: float + bias: float + std: float + rmse: float + crlb_std: float + crlb_variance: float + std_over_crlb_std: float + + +@dataclass(frozen=True) +class ModeValidation: + """One true mode: recovery count, loss rate and per-parameter + statistics computed over the recovered trials only.""" + + angular_frequency: ParameterStats + decay: ParameterStats + amplitude: ParameterStats + phase: ParameterStats + recovered_trials: int + loss_rate: float + + @property + def frequency_ratio(self) -> float: + return self.angular_frequency.std_over_crlb_std + + @property + def decay_ratio(self) -> float: + return self.decay.std_over_crlb_std + + +@dataclass(frozen=True) +class ValidationLevel: + """One sweep condition. A trial is successful when the estimator + returned and every true mode was matched; ``estimator_errors`` counts + trials in which the estimator raised, which are also failed trials.""" + + snr_db: float + noise_sigma: float + trials_attempted: int + trials_successful: int + trials_failed: int + estimator_errors: int + n_components: int + modes: tuple[ModeValidation, ...] + any_mode_lost_rate: float + residual_ratio: float + + +@dataclass(frozen=True) +class ValidationSweep: + samples: int + duration: float + signal_rms: float + levels: tuple[ValidationLevel, ...] + match_tolerance: float + backend: str + modes: tuple[DampedMode, ...] = DEFAULT_MODES + constant: float = 0.15 + seed: int = 0 + distortion: tuple[str, float] | None = None + extra: dict = field(default_factory=dict) + + def to_dict(self) -> dict: + return asdict(self) + + +def _wrap_phase(values: np.ndarray) -> np.ndarray: + return (values + np.pi) % (2.0 * np.pi) - np.pi + + +def _canonical(amplitude: np.ndarray, phase: np.ndarray): + """A negative amplitude is the same mode with the phase shifted by pi.""" + amplitude = np.array(amplitude, dtype=float, copy=True) + phase = np.array(phase, dtype=float, copy=True) + neg = amplitude < 0 + amplitude[neg] = -amplitude[neg] + phase[neg] = phase[neg] + np.pi + return amplitude, _wrap_phase(phase) + + +def _stats( + samples: np.ndarray, truth: float, crlb_std: float, *, wrap: bool +) -> ParameterStats: + err = samples - truth + if wrap: + err = _wrap_phase(err) + n = err.size + bias = float(err.mean()) if n else float("nan") + std = float(err.std(ddof=1)) if n > 2 else float("nan") + rmse = float(np.sqrt(np.mean(err**2))) if n else float("nan") + return ParameterStats( + truth=float(truth), + bias=bias, + std=std, + rmse=rmse, + crlb_std=float(crlb_std), + crlb_variance=float(crlb_std) ** 2, + std_over_crlb_std=std / float(crlb_std), + ) + + +def validation_sweep( + *, + samples: int = 96, + modes: Sequence[DampedMode] = DEFAULT_MODES, + constant: float = 0.15, + duration: float = 9.5, + snr_db: Sequence[float] = (40.0, 30.0, 20.0, 10.0), + trials: int = 200, + n_components: int = 6, + match_tolerance: float = 0.3, + seed: int = 20260914, + estimator: Callable | None = None, + backend: str = "cpu", + distortion: tuple[str, float] | None = None, +) -> ValidationSweep: + """Monte Carlo sweep of the linear-prediction estimator against the + Cramer-Rao bound. + + ``snr_db`` is defined against the rms of the mean-removed clean trace. + For each level the per-mode scatter, bias and rmse of every parameter + are compared with the bound, and ``loss_rate`` is the fraction of + trials in which the true mode had no fitted mode within + ``match_tolerance`` (radians per time unit) of its angular frequency; + each fitted mode can match one true mode. + + The statistics are computed over the recovered trials only and the + bound is unconditional. The root filter keeps decaying roots, so the + surviving decay estimates are truncated at zero; for a lightly damped + mode the reported decay statistics are therefore conditional on + survival and can move non-monotonically with noise. Read them together + with ``loss_rate``. + + ``estimator`` may return CuPy arrays; they are moved to the host. A + trial in which the estimator raises is counted in ``estimator_errors`` + and as a loss of every mode. + + ``distortion`` = (kind, amount) replaces the clean trace with one from + :func:`distort_trace`, outside the model class; ``residual_ratio`` is + then the median rms residual of the reconstruction divided by the + noise sigma. It is close to 1 when the model fits and well above 1 + when the trace is not a sum of damped sinusoids. + """ + if trials < 1: + raise ValueError("trials must be at least 1") + if not snr_db: + raise ValueError("snr_db must contain at least one level") + est = estimator or (lambda t, y, k: linear_prediction_numpy(t, y, k)) + base = synthetic_modes_trace( + samples, modes=modes, constant=constant, duration=duration + ) + clean = base.clean + if distortion is not None: + clean = distort_trace( + base.time, modes, constant, distortion[0], distortion[1] + ) + signal_rms = float(np.std(clean - clean.mean())) + rng = np.random.default_rng(seed) + levels = [] + for snr in snr_db: + sigma = signal_rms / 10 ** (snr / 20) + bounds = cramer_rao_bounds(modes, constant, base.time, sigma) + found: list[dict[str, list[float]]] = [ + {name: [] for name in PARAMETERS} for _ in modes + ] + lost = np.zeros(len(modes), dtype=int) + any_lost = 0 + errors = 0 + ratios = [] + for _ in range(trials): + y = clean + rng.normal(0.0, sigma, size=base.time.shape) + try: + r = est(base.time, y, n_components) + except Exception: + errors += 1 + any_lost += 1 + lost += 1 + continue + w = _to_host(r.angular_frequency) + d = _to_host(r.decay) + a, p = _canonical(_to_host(r.amplitude), _to_host(r.phase)) + rec = _to_host(r.reconstruction) + ratios.append(float(np.sqrt(np.mean((y - rec) ** 2)) / sigma)) + hit = match_modes(modes, w, d, tolerance=match_tolerance) + if any(h is None for h in hit): + any_lost += 1 + for k, h in enumerate(hit): + if h is None: + lost[k] += 1 + continue + found[k]["amplitude"].append(a[h]) + found[k]["decay"].append(d[h]) + found[k]["angular_frequency"].append(w[h]) + found[k]["phase"].append(p[h]) + per_mode = [] + for k, m in enumerate(modes): + stats = { + name: _stats( + np.asarray(found[k][name]), + getattr(m, name), + bounds[name][k], + wrap=name == "phase", + ) + for name in PARAMETERS + } + per_mode.append( + ModeValidation( + angular_frequency=stats["angular_frequency"], + decay=stats["decay"], + amplitude=stats["amplitude"], + phase=stats["phase"], + recovered_trials=int(trials - lost[k]), + loss_rate=float(lost[k] / trials), + ) + ) + levels.append( + ValidationLevel( + snr_db=float(snr), + noise_sigma=float(sigma), + trials_attempted=trials, + trials_successful=trials - any_lost, + trials_failed=any_lost, + estimator_errors=errors, + n_components=n_components, + modes=tuple(per_mode), + any_mode_lost_rate=float(any_lost / trials), + residual_ratio=float(np.median(ratios)) + if ratios + else float("nan"), + ) + ) + return ValidationSweep( + samples=samples, + duration=duration, + signal_rms=signal_rms, + levels=tuple(levels), + match_tolerance=match_tolerance, + backend=backend, + modes=tuple(modes), + constant=constant, + seed=seed, + distortion=distortion, + ) + + +def source_revision() -> str | None: + """Git revision of the checkout containing this module, if any.""" + try: + out = subprocess.run( + ["git", "rev-parse", "HEAD"], + cwd=Path(__file__).resolve().parent, + capture_output=True, + text=True, + timeout=5, + ) + except (OSError, subprocess.SubprocessError): + return None + return out.stdout.strip() if out.returncode == 0 else None + + +def build_summary( + sweeps: Sequence[ValidationSweep], + *, + command: str = "linear-prediction-validate", + artifacts: dict[str, str | None] | None = None, +) -> dict[str, Any]: + """Assemble the ``schema_version`` 1 summary of one or more sweeps. + + ``config`` records the truth, the sampling, the noise definition, the + seed, the trial count, the model order and the mode-matching rule; + ``runtime`` the package, Python, NumPy and SciPy versions, the source + revision and the backend; ``results`` one entry per sweep condition + (estimator by distortion by signal-to-noise level); ``artifacts`` the + relative filenames written next to the summary. + """ + if not sweeps: + raise ValueError("at least one sweep is required") + first = sweeps[0] + dt = first.duration / (first.samples - 1) + config = { + "model": ( + "trace(t) = constant + sum_k amplitude_k exp(-decay_k t) " + "cos(angular_frequency_k t + phase_k)" + ), + "units": { + "time": "input time unit", + "amplitude": "trace units", + "decay": "1 / time unit", + "angular_frequency": "rad / time unit", + "phase": "rad", + }, + "truth": { + "modes": [asdict(m) for m in first.modes], + "constant": first.constant, + }, + "sampling": { + "samples": first.samples, + "duration": first.duration, + "dt": dt, + "start": 0.0, + }, + "noise": { + "kind": "white Gaussian, independent per sample", + "snr_db_definition": ( + "20 log10(rms of the mean-removed clean trace / sigma)" + ), + "signal_rms": first.signal_rms, + "snr_db": [lvl.snr_db for lvl in first.levels], + }, + "seed": first.seed, + "trials_per_level": first.levels[0].trials_attempted, + "n_components": first.levels[0].n_components, + "mode_matching": { + "criterion": ( + "a true mode is recovered when a fitted mode lies within " + "the tolerance of its angular frequency; fitted modes are " + "assigned nearest-first and each is used at most once" + ), + "tolerance_rad_per_time": first.match_tolerance, + }, + "statistics": ( + "bias, std (ddof 1) and rmse are computed over the recovered " + "trials of each mode; crlb_std and crlb_variance are the " + "unconditional Cramer-Rao bounds from the exact Fisher " + "information of the model at the true parameters; the phase " + "error is wrapped to [-pi, pi)" + ), + "estimators": sorted({s.backend for s in sweeps}), + "distortions": sorted( + { + f"{s.distortion[0]}:{s.distortion[1]:g}" + if s.distortion + else "none" + for s in sweeps + } + ), + } + runtime = runtime_metadata(backend="cpu", dtype="float64") + try: + from scipy import __version__ as scipy_version + except ImportError: # pragma: no cover - scipy is a dependency + scipy_version = None + runtime["scipy_version"] = scipy_version + runtime["source_revision"] = source_revision() + results = [] + for sweep in sweeps: + for lvl in sweep.levels: + results.append( + { + "estimator": sweep.backend, + "distortion": ( + { + "kind": sweep.distortion[0], + "amount": sweep.distortion[1], + } + if sweep.distortion + else None + ), + "snr_db": lvl.snr_db, + "noise_sigma": lvl.noise_sigma, + "trials": { + "attempted": lvl.trials_attempted, + "successful": lvl.trials_successful, + "failed": lvl.trials_failed, + "estimator_errors": lvl.estimator_errors, + }, + "any_mode_lost_rate": lvl.any_mode_lost_rate, + "residual_rms_over_sigma_median": lvl.residual_ratio, + "modes": [ + { + "recovered_trials": m.recovered_trials, + "loss_rate": m.loss_rate, + "statistics_over": "recovered_trials", + **{ + name: asdict(getattr(m, name)) + for name in PARAMETERS + }, + } + for m in lvl.modes + ], + } + ) + return { + "schema_version": SUMMARY_SCHEMA_VERSION, + "command": command, + "config": config, + "runtime": runtime, + "results": results, + "artifacts": dict(artifacts or {}), + } + + +def write_validation_figure( + sweeps: Sequence[ValidationSweep], path: Path, *, title: str +) -> bool: + """Standalone Bokeh HTML: per estimator and parameter the sample + standard deviation against the bound versus signal-to-noise, and the + loss rate. Returns False when the ``viz`` extra is not installed.""" + try: + from bokeh.layouts import gridplot + from bokeh.plotting import figure, output_file, save + except ImportError: + return False + palette = ["#1f77b4", "#d62728", "#2ca02c", "#9467bd", "#ff7f0e"] + rows = [] + for sweep in sweeps: + snr = [lvl.snr_db for lvl in sweep.levels] + label = sweep.backend + ( + f" {sweep.distortion[0]}:{sweep.distortion[1]:g}" + if sweep.distortion + else "" + ) + row = [] + for name in ("angular_frequency", "decay"): + fig = figure( + width=360, + height=260, + y_axis_type="log", + title=f"{label}: {name} std (points) vs bound (dashed)", + x_axis_label="snr (dB)", + y_axis_label="std", + ) + for k, mode in enumerate(sweep.modes): + std = [ + getattr(lvl.modes[k], name).std for lvl in sweep.levels + ] + crlb = [ + getattr(lvl.modes[k], name).crlb_std + for lvl in sweep.levels + ] + color = palette[k % len(palette)] + fig.line(snr, crlb, color=color, line_dash="dashed") + fig.scatter( + snr, + std, + color=color, + size=7, + legend_label=f"mode w={mode.angular_frequency:g}", + ) + fig.legend.label_text_font_size = "8pt" + row.append(fig) + fig = figure( + width=360, + height=260, + title=f"{label}: loss rate", + x_axis_label="snr (dB)", + y_axis_label="fraction of trials", + ) + for k, mode in enumerate(sweep.modes): + fig.line( + snr, + [lvl.modes[k].loss_rate for lvl in sweep.levels], + color=palette[k % len(palette)], + legend_label=f"mode w={mode.angular_frequency:g}", + ) + fig.line( + snr, + [lvl.any_mode_lost_rate for lvl in sweep.levels], + color="black", + line_dash="dotted", + legend_label="any mode", + ) + fig.legend.label_text_font_size = "8pt" + row.append(fig) + rows.append(row) + output_file(str(path), title=title) + save(gridplot(rows)) + return True + + +def write_validation_run( + output_dir: Path | str, + sweeps: Sequence[ValidationSweep], + *, + command: str = "linear-prediction-validate", + figure_name: str = "validation.html", +) -> dict[str, Any]: + """Write ``summary.json`` and, when Bokeh is available, the figure + into ``output_dir`` (created if needed). Returns the summary.""" + out = Path(output_dir) + out.mkdir(parents=True, exist_ok=True) + have_figure = write_validation_figure( + sweeps, out / figure_name, title="xray linear prediction validation" + ) + artifacts = { + "summary": "summary.json", + "figure": figure_name if have_figure else None, + } + summary = build_summary(sweeps, command=command, artifacts=artifacts) + with open(out / "summary.json", "w", encoding="utf-8", newline="\n") as f: + json.dump(summary, f, indent=2, sort_keys=True) + f.write("\n") + return summary diff --git a/tests/core/test_cli_contract.py b/tests/core/test_cli_contract.py index ad4fdec1..3ef79c25 100644 --- a/tests/core/test_cli_contract.py +++ b/tests/core/test_cli_contract.py @@ -75,13 +75,13 @@ def test_public_command_surface_counts_are_exact() -> None: ("xpois", 7, 7, 162, 1), ("xscan", 44, 44, 229, 1), ("xrep", 6, 6, 123, 1), - ("xray", 33, 31, 392, 1), + ("xray", 34, 32, 400, 1), ] assert len(per_group) == 6 - assert sum(item[1] for item in per_group) == 94 - assert sum(item[1] + item[4] for item in per_group) == 99 - assert sum(item[2] for item in per_group) == 91 - assert sum(item[3] for item in per_group) == 961 + assert sum(item[1] for item in per_group) == 95 + assert sum(item[1] + item[4] for item in per_group) == 100 + assert sum(item[2] for item in per_group) == 92 + assert sum(item[3] for item in per_group) == 969 def test_public_registry_order_and_component_derivations() -> None: diff --git a/tests/xray/test_cli.py b/tests/xray/test_cli.py index 10b78bdf..7ea68a11 100644 --- a/tests/xray/test_cli.py +++ b/tests/xray/test_cli.py @@ -64,6 +64,7 @@ def test_core_command_registry_loads_xray_commands(): "extract-trace": "et", "gpu-policy": "gp", "linear-prediction-benchmark": "lpb", + "linear-prediction-validate": "lpv", "linear-prediction-fixed-stages-benchmark": "lpfsb", "linear-prediction-p2-benchmark": "lppb", "linear-prediction-profile-summary": "lpps", @@ -205,6 +206,56 @@ def test_detector_artifacts_cli_rejects_unknown_diagnostic_level( assert "invalid choice: 'verbose'" in captured.err +def test_linear_prediction_validate_text_json_and_output_dir( + tmp_path, capsys +): + # the optional --distortion string must be allowed to be absent + assert main(["lpv", "--trials", "3", "--snr-db", "30"]) == 0 + assert "loss_rate=" in capsys.readouterr().out + assert ( + main( + [ + "lpv", + "--trials", + "3", + "--snr-db", + "30", + "--distortion", + "chirp:0.05", + "--json", + ] + ) + == 0 + ) + payload = json.loads(capsys.readouterr().out) + assert payload["schema_version"] == 1 + assert payload["results"][0]["distortion"] == { + "kind": "chirp", + "amount": 0.05, + } + assert payload["artifacts"] == {} + out = tmp_path / "lpv" + assert ( + main( + [ + "lpv", + "--trials", + "3", + "--snr-db", + "30", + "--output-dir", + str(out), + ] + ) + == 0 + ) + text = capsys.readouterr().out + assert "summary=summary.json" in text + written = json.loads((out / "summary.json").read_text()) + assert written["results"][0]["trials"]["attempted"] == 3 + assert written["artifacts"]["summary"] == "summary.json" + + def test_gpu_policy(capsys): assert main(["gpu-policy"]) == 0 captured = capsys.readouterr() diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py new file mode 100644 index 00000000..59fca9d6 --- /dev/null +++ b/tests/xray/test_synthetic_validation.py @@ -0,0 +1,207 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +import numpy as np +import pytest + +from cuphoton.xray.linear_prediction import ( + linear_prediction_numpy, + synthetic_trace, +) +from cuphoton.xray.synthetic_validation import ( + PARAMETERS, + DampedMode, + build_summary, + cramer_rao_bounds, + distort_trace, + match_modes, + synthetic_modes_trace, + validation_sweep, + write_validation_run, +) + + +def test_fixture_reproduces_synthetic_trace_at_zero_noise(): + time, trace = synthetic_trace(96) + fx = synthetic_modes_trace(96) + np.testing.assert_allclose(fx.time, time) + np.testing.assert_allclose(fx.trace, trace) + assert fx.noise_sigma == 0.0 + + +def test_fixture_noise_is_reproducible_and_has_the_requested_scale(): + a = synthetic_modes_trace(256, noise_sigma=0.05, seed=7) + b = synthetic_modes_trace(256, noise_sigma=0.05, seed=7) + np.testing.assert_array_equal(a.trace, b.trace) + resid = a.trace - a.clean + assert abs(resid.std() - 0.05) < 0.01 + + +def test_estimator_recovers_known_modes_exactly_at_zero_noise(): + fx = synthetic_modes_trace(96) + r = linear_prediction_numpy(fx.time, fx.trace, 6) + hit = match_modes( + fx.modes, + np.asarray(r.angular_frequency), + np.asarray(r.decay), + tolerance=0.05, + ) + assert all(h is not None for h in hit) + for m, h in zip(fx.modes, hit): + assert abs(float(r.angular_frequency[h]) - m.angular_frequency) < 1e-6 + assert abs(float(r.decay[h]) - m.decay) < 1e-6 + + +def test_cramer_rao_bound_scales_with_noise_and_shrinks_with_trace_length(): + fx = synthetic_modes_trace(96) + b1 = cramer_rao_bounds(fx.modes, fx.constant, fx.time, 0.01) + b2 = cramer_rao_bounds(fx.modes, fx.constant, fx.time, 0.02) + np.testing.assert_allclose( + b2["angular_frequency"], 2 * b1["angular_frequency"], rtol=1e-6 + ) + longer = synthetic_modes_trace(192, duration=19.0) + b3 = cramer_rao_bounds(longer.modes, longer.constant, longer.time, 0.01) + assert np.all(b3["angular_frequency"] < b1["angular_frequency"]) + + +def test_single_mode_bound_matches_the_undamped_closed_form(): + # Real sinusoid A cos(w t + phi) in white Gaussian noise, unknown A, w and + # phi, samples at t = n dt for n = 0..N-1 (Kay, Fundamentals of + # Statistical Signal Processing: Estimation Theory, 1993, sinusoidal + # parameter estimation): var(w) >= 24 sigma^2 / (A^2 dt^2 N (N^2 - 1)). + # The complex-exponential form has 12 in place of 24. The numerical + # bound also frees the decay and the constant, so it may sit slightly + # above the closed form, never below it. + mode = (DampedMode(1.0, 0.0, 2.0, 0.3),) + fx = synthetic_modes_trace(256, modes=mode, constant=0.0, duration=25.5) + dt = fx.time[1] - fx.time[0] + n = fx.time.size + sigma = 0.01 + closed = np.sqrt(24 * sigma**2 / (1.0**2 * dt**2 * n * (n**2 - 1))) + numerical = cramer_rao_bounds(mode, 0.0, fx.time, sigma)[ + "angular_frequency" + ][0] + assert numerical >= closed * 0.999 + assert numerical < closed * 1.03 + + +def test_validation_sweep_reports_bound_ratios_and_loss_rates(): + sweep = validation_sweep(snr_db=(40.0, 20.0), trials=40, seed=1) + assert len(sweep.levels) == 2 + strong = sweep.levels[0].modes[0] + w = strong.angular_frequency + assert w.crlb_std > 0 and w.std > 0 + assert w.crlb_variance == pytest.approx(w.crlb_std**2) + assert w.std_over_crlb_std >= 1.0 # an estimator cannot beat the bound + assert w.std_over_crlb_std < 10.0 + assert abs(w.bias) < 0.01 + assert w.rmse >= abs(w.bias) + for name in PARAMETERS: + assert getattr(strong, name).truth == getattr(sweep.modes[0], name) + assert abs(strong.phase.bias) < 0.05 + for level in sweep.levels: + assert level.trials_attempted == 40 + assert level.trials_successful + level.trials_failed == 40 + assert level.estimator_errors == 0 + for m in level.modes: + assert 0.0 <= m.loss_rate <= 1.0 + assert m.recovered_trials == round(40 * (1 - m.loss_rate)) + d = sweep.to_dict() + assert d["levels"][0]["snr_db"] == 40.0 + + +def test_estimator_errors_count_as_failed_trials(): + def broken(t, y, k): + raise RuntimeError("no fit") + + sweep = validation_sweep(snr_db=(30.0,), trials=3, estimator=broken) + level = sweep.levels[0] + assert level.estimator_errors == 3 + assert level.trials_failed == 3 and level.trials_successful == 0 + assert all(m.loss_rate == 1.0 for m in level.modes) + assert np.isnan(level.residual_ratio) + + +def test_summary_schema_and_run_artifacts(tmp_path): + sweep = validation_sweep(snr_db=(30.0,), trials=6, seed=3) + summary = build_summary([sweep]) + assert summary["schema_version"] == 1 + assert summary["command"] == "linear-prediction-validate" + config = summary["config"] + assert config["truth"]["modes"][0]["angular_frequency"] == 2.4 + assert config["sampling"]["samples"] == 96 + assert config["seed"] == 3 and config["trials_per_level"] == 6 + assert config["n_components"] == 6 + assert config["mode_matching"]["tolerance_rad_per_time"] == 0.3 + assert summary["runtime"]["backend"] == "cpu" + assert summary["runtime"]["numpy_version"] + result = summary["results"][0] + assert result["trials"]["attempted"] == 6 + mode = result["modes"][0] + assert mode["statistics_over"] == "recovered_trials" + for name in PARAMETERS: + assert set(mode[name]) == { + "truth", + "bias", + "std", + "rmse", + "crlb_std", + "crlb_variance", + "std_over_crlb_std", + } + out = tmp_path / "run" + written = write_validation_run(out, [sweep]) + assert (out / "summary.json").exists() + assert written["artifacts"]["summary"] == "summary.json" + figure = written["artifacts"]["figure"] + assert figure is None or (out / figure).exists() + with pytest.raises(ValueError): + build_summary([]) + + +def test_invalid_inputs_are_rejected(): + fx = synthetic_modes_trace(96) + with pytest.raises(ValueError): + cramer_rao_bounds(fx.modes, fx.constant, fx.time, 0.0) + with pytest.raises(ValueError): + validation_sweep(trials=0, snr_db=(40.0,)) + with pytest.raises(ValueError): + validation_sweep(trials=5, snr_db=()) + with pytest.raises(ValueError): + synthetic_modes_trace(8) + + +def test_undistorted_residual_ratio_is_near_one(): + sweep = validation_sweep(snr_db=(30.0,), trials=30, seed=9) + assert 0.7 < sweep.levels[0].residual_ratio < 1.6 + + +def test_distortions_raise_the_residual_ratio(): + for kind, amount in [ + ("chirp", 0.2), + ("gaussian_envelope", 6.0), + ("baseline_drift", 0.3), + ("clip", 1.0), + ]: + sweep = validation_sweep( + snr_db=(30.0,), trials=20, seed=9, distortion=(kind, amount) + ) + assert sweep.levels[0].residual_ratio > 4.0, kind + assert sweep.distortion == (kind, amount) + + +def test_distort_trace_kinds_and_rejection(): + fx = synthetic_modes_trace(96) + for kind, amount in [ + ("chirp", 0.1), + ("gaussian_envelope", 5.0), + ("baseline_drift", 0.2), + ("clip", 1.0), + ("glitch", 0.5), + ]: + y = distort_trace(fx.time, fx.modes, fx.constant, kind, amount) + assert y.shape == fx.time.shape + assert not np.allclose(y, fx.clean), kind + with pytest.raises(ValueError): + distort_trace(fx.time, fx.modes, fx.constant, "square", 1.0) From f8037ebd010ad88a5d964a947cf8be72e19805e9 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Fri, 18 Sep 2026 20:29:11 -0700 Subject: [PATCH 02/10] Correct synthetic validation statistics and provenance Keep small-trial summaries valid JSON, integrate the requested chirp, and retain floating parameters for Fisher derivatives. Match nearby modes globally and reject incompatible combined sweep settings. Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-VALIDATION.md | 19 +++- src/cuphoton/xray/commands.py | 2 +- src/cuphoton/xray/synthetic_validation.py | 109 ++++++++++++++++------ tests/xray/test_synthetic_validation.py | 84 +++++++++++++++++ 4 files changed, 181 insertions(+), 33 deletions(-) diff --git a/docs/xray/LINEAR-PREDICTION-VALIDATION.md b/docs/xray/LINEAR-PREDICTION-VALIDATION.md index 510630e7..6470528b 100644 --- a/docs/xray/LINEAR-PREDICTION-VALIDATION.md +++ b/docs/xray/LINEAR-PREDICTION-VALIDATION.md @@ -53,7 +53,8 @@ within 3 percent, the difference being the freed decay and constant. A true mode is recovered in a trial when a fitted mode has an angular frequency within `match_tolerance` (default 0.3 rad per time unit) of it; -fitted modes are assigned nearest-first and each is used at most once. The +the assignment first maximizes the number of recovered modes, then minimizes +their total frequency error. Each fitted mode is used at most once. The per-mode `loss_rate` is the fraction of trials in which the mode was not recovered. A trial is successful when the estimator returned and every true mode was recovered; a trial in which the estimator raised counts in @@ -64,6 +65,8 @@ For each mode and parameter (`amplitude`, `decay`, `angular_frequency`, the recovered trials only, the bound (`crlb_std`, `crlb_variance`) and `std_over_crlb_std`. A negative fitted amplitude is folded into the phase and the phase error is wrapped to `[-pi, pi)`. +Undefined statistics are written as JSON `null`: standard deviation requires +at least two recovered trials, and bias and RMSE require at least one. Read the decay statistics of a lightly damped mode together with its `loss_rate`. The root filter keeps decaying roots only, so when noise pushes @@ -74,7 +77,9 @@ labels these statistics `statistics_over: recovered_trials`. Every level also reports `residual_rms_over_sigma_median`, the median over trials of the rms residual of the reconstruction divided by sigma. It is -close to 1 when the model class fits and well above 1 when it does not. +near 1 for an accurate reconstruction at the specified noise level. A large +value can indicate model mismatch, missed modes, or estimator error; it does +not by itself distinguish these causes. ## Model mismatch @@ -92,7 +97,7 @@ model class: The truth used for bias and bounds stays the undistorted mode set, so the statistics measure the estimator's response to the mismatch. A chirp or a clipped trace can return modes with a confident, biased frequency and no -loss; the residual ratio is the number that says so. +loss; the residual ratio supplies an additional reconstruction diagnostic. ## Summary schema @@ -109,7 +114,7 @@ loss; the residual ratio is the number that says so. `scipy_version` and `source_revision` (the git revision of the checkout, or `null`). - `results`: one entry per sweep condition (estimator by distortion by - signal-to-noise level): `snr_db`, `noise_sigma`, `trials` with + signal-to-noise level): `snr_db`, `signal_rms`, `noise_sigma`, `trials` with `attempted`, `successful`, `failed` and `estimator_errors`, `any_mode_lost_rate`, `residual_rms_over_sigma_median`, and `modes`, each with `recovered_trials`, `loss_rate`, `statistics_over` and the four @@ -119,6 +124,8 @@ loss; the residual ratio is the number that says so. The figure shows, per estimator and per parameter, the sample standard deviation against the bound versus signal-to-noise, and the loss rates. +The signal RMS in `config.noise` describes the first sweep; each result +records its own signal RMS because distortions can change it. ## Python entry points @@ -129,3 +136,7 @@ deviation against the bound versus signal-to-noise, and the loss rates. `estimator(time, trace, n_components)` returning an object with `angular_frequency`, `decay`, `amplitude`, `phase` and `reconstruction` (NumPy or CuPy arrays), so the same sweep can compare estimators. +`build_summary` combines sweeps with shared truth, sampling, seed, matching, +SNR levels, trial counts and model order. It rejects incompatible settings +instead of describing later sweeps with the first sweep's configuration. +Write separate summaries for changes to those settings. diff --git a/src/cuphoton/xray/commands.py b/src/cuphoton/xray/commands.py index 6929b47e..a09273cb 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -4472,7 +4472,7 @@ def _linear_prediction_validate(args): else: summary = build_summary(sweeps) if args.json: - print(json.dumps(summary, indent=2, sort_keys=True)) + print(json.dumps(summary, indent=2, sort_keys=True, allow_nan=False)) return 0 print(f"samples={sweeps[0].samples}") print(f"signal_rms={sweeps[0].signal_rms:.6g}") diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index 1db5a48c..5fb62b96 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -23,6 +23,7 @@ from typing import Any, Callable, Sequence import numpy as np +from scipy.optimize import linear_sum_assignment from cuphoton.core.runtime import runtime_metadata @@ -85,7 +86,8 @@ def modes_to_theta( for m in modes for v in (m.amplitude, m.decay, m.angular_frequency, m.phase) ] - + [constant] + + [constant], + dtype=np.float64, ) @@ -141,21 +143,27 @@ def distort_trace( clipped above ``amount``. ``glitch``: one sample at mid-record is offset by ``amount``. These are outside the model class on purpose: the sweep reports how the fit responds and the residual-to-noise ratio - that should flag them. + that can help diagnose them. """ if kind not in DISTORTIONS: raise ValueError(f"unknown distortion {kind!r}") span = float(time[-1] - time[0]) out = np.full(time.shape, float(constant)) for k, m in enumerate(modes): - w = m.angular_frequency + phase = m.angular_frequency * time + m.phase if kind == "chirp": - w = w * (1.0 + amount * (time - time[0]) / span) + phase = phase + ( + 0.5 + * m.angular_frequency + * amount + * (time - time[0]) ** 2 + / span + ) if kind == "gaussian_envelope" and k == 0: env = np.exp(-((time / amount) ** 2)) else: env = np.exp(-m.decay * time) - out = out + m.amplitude * env * np.cos(w * time + m.phase) + out = out + m.amplitude * env * np.cos(phase) if kind == "baseline_drift": out = out + amount * (time - time[0]) / span elif kind == "clip": @@ -223,22 +231,28 @@ def match_modes( """For each true mode, the index of the fitted mode within ``tolerance`` (rad per time unit) or ``None``. """ - out: list[int | None] = [] - used: set[int] = set() - for m in modes: - if angular_frequency.size == 0: - out.append(None) - continue - err = np.abs(angular_frequency - m.angular_frequency) - for idx in np.argsort(err): - if idx in used: - continue - if err[idx] <= tolerance: - out.append(int(idx)) - used.add(int(idx)) - break - else: - out.append(None) + if not np.isfinite(tolerance) or tolerance < 0: + raise ValueError("tolerance must be finite and nonnegative") + out: list[int | None] = [None] * len(modes) + if not modes or angular_frequency.size == 0: + return out + err = np.abs( + np.array([m.angular_frequency for m in modes])[:, None] + - angular_frequency[None, :] + ) + # A missing match costs more than every valid distance combined: + # maximize recovered modes first, then minimize total frequency error. + missing_cost = len(modes) + 1 + cost = np.full( + (len(modes), angular_frequency.size + len(modes)), float(missing_cost) + ) + cost[:, : angular_frequency.size] = np.where( + err <= tolerance, err / (tolerance or 1.0), np.inf + ) + rows, columns = linear_sum_assignment(cost) + for row, column in zip(rows, columns): + if column < angular_frequency.size: + out[int(row)] = int(column) return out @@ -340,7 +354,7 @@ def _stats( err = _wrap_phase(err) n = err.size bias = float(err.mean()) if n else float("nan") - std = float(err.std(ddof=1)) if n > 2 else float("nan") + std = float(err.std(ddof=1)) if n > 1 else float("nan") rmse = float(np.sqrt(np.mean(err**2))) if n else float("nan") return ParameterStats( truth=float(truth), @@ -392,8 +406,8 @@ def validation_sweep( ``distortion`` = (kind, amount) replaces the clean trace with one from :func:`distort_trace`, outside the model class; ``residual_ratio`` is then the median rms residual of the reconstruction divided by the - noise sigma. It is close to 1 when the model fits and well above 1 - when the trace is not a sum of damped sinusoids. + noise sigma. Large values can reflect model mismatch, missed modes or + estimator error; this ratio alone does not distinguish those causes. """ if trials < 1: raise ValueError("trials must be at least 1") @@ -530,6 +544,34 @@ def build_summary( if not sweeps: raise ValueError("at least one sweep is required") first = sweeps[0] + shared_fields = ( + "samples", + "duration", + "modes", + "constant", + "seed", + "match_tolerance", + ) + first_levels = [ + (level.snr_db, level.trials_attempted, level.n_components) + for level in first.levels + ] + for sweep in sweeps[1:]: + if ( + any( + getattr(sweep, name) != getattr(first, name) + for name in shared_fields + ) + or [ + (level.snr_db, level.trials_attempted, level.n_components) + for level in sweep.levels + ] + != first_levels + ): + raise ValueError( + "summary sweeps must share truth, sampling, seed, " + "matching and trial settings" + ) dt = first.duration / (first.samples - 1) config = { "model": ( @@ -568,7 +610,8 @@ def build_summary( "criterion": ( "a true mode is recovered when a fitted mode lies within " "the tolerance of its angular frequency; fitted modes are " - "assigned nearest-first and each is used at most once" + "assigned to maximize matches, then minimize total frequency " + "error, and each is used at most once" ), "tolerance_rad_per_time": first.match_tolerance, }, @@ -611,6 +654,7 @@ def build_summary( else None ), "snr_db": lvl.snr_db, + "signal_rms": sweep.signal_rms, "noise_sigma": lvl.noise_sigma, "trials": { "attempted": lvl.trials_attempted, @@ -619,14 +663,23 @@ def build_summary( "estimator_errors": lvl.estimator_errors, }, "any_mode_lost_rate": lvl.any_mode_lost_rate, - "residual_rms_over_sigma_median": lvl.residual_ratio, + "residual_rms_over_sigma_median": ( + lvl.residual_ratio + if np.isfinite(lvl.residual_ratio) + else None + ), "modes": [ { "recovered_trials": m.recovered_trials, "loss_rate": m.loss_rate, "statistics_over": "recovered_trials", **{ - name: asdict(getattr(m, name)) + name: { + key: value if np.isfinite(value) else None + for key, value in asdict( + getattr(m, name) + ).items() + } for name in PARAMETERS }, } @@ -742,6 +795,6 @@ def write_validation_run( } summary = build_summary(sweeps, command=command, artifacts=artifacts) with open(out / "summary.json", "w", encoding="utf-8", newline="\n") as f: - json.dump(summary, f, indent=2, sort_keys=True) + json.dump(summary, f, indent=2, sort_keys=True, allow_nan=False) f.write("\n") return summary diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py index 59fca9d6..44eef03e 100644 --- a/tests/xray/test_synthetic_validation.py +++ b/tests/xray/test_synthetic_validation.py @@ -2,6 +2,9 @@ # # SPDX-License-Identifier: Apache-2.0 +import json +from dataclasses import replace + import numpy as np import pytest @@ -86,6 +89,41 @@ def test_single_mode_bound_matches_the_undamped_closed_form(): assert numerical < closed * 1.03 +def test_integer_mode_parameters_have_the_same_bounds_as_floats(): + time = np.linspace(0, 25.5, 256) + actual = cramer_rao_bounds((DampedMode(1, 0, 2, 0),), 0, time, 0.01) + expected = cramer_rao_bounds( + (DampedMode(1.0, 0.0, 2.0, 0.0),), 0.0, time, 0.01 + ) + for name in (*PARAMETERS, "constant"): + np.testing.assert_allclose(actual[name], expected[name]) + assert np.all(np.isfinite(actual[name])) + + +@pytest.mark.parametrize("start", [0.0, 3.0]) +def test_chirp_integrates_the_requested_frequency_ramp(start): + time = np.linspace(start, start + 10.0, 256) + # Integrating omega(t) = 2 + .04 * (t - start) gives this phase. + expected = np.cos(2 * time + 0.02 * (time - start) ** 2 + 0.3) + actual = distort_trace(time, (DampedMode(1, 0, 2, 0.3),), 0, "chirp", 0.2) + np.testing.assert_allclose(actual, expected) + + +def test_matching_recovers_close_modes_independent_of_truth_order(): + modes = (DampedMode(1, 0, 1.0), DampedMode(1, 0, 1.15)) + fitted = np.array([1.05, 0.85]) + decay = np.zeros(2) + assert match_modes(modes, fitted, decay, tolerance=0.2) == [1, 0] + assert match_modes(modes[::-1], fitted, decay, tolerance=0.2) == [0, 1] + assert match_modes(modes, fitted[:1], decay[:1], tolerance=0.2) == [ + 0, + None, + ] + assert match_modes( + modes, np.array([np.nan]), decay[:1], tolerance=0.2 + ) == [None, None] + + def test_validation_sweep_reports_bound_ratios_and_loss_rates(): sweep = validation_sweep(snr_db=(40.0, 20.0), trials=40, seed=1) assert len(sweep.levels) == 2 @@ -160,6 +198,52 @@ def test_summary_schema_and_run_artifacts(tmp_path): build_summary([]) +@pytest.mark.parametrize("trials", [1, 2]) +def test_small_sweeps_write_strict_json_and_two_trials_have_std( + tmp_path, trials +): + sweep = validation_sweep(snr_db=(30.0,), trials=trials, seed=3) + summary = write_validation_run(tmp_path, [sweep]) + encoded = json.dumps(summary, allow_nan=False) + assert json.loads(encoded) == json.loads( + (tmp_path / "summary.json").read_text() + ) + strong = summary["results"][0]["modes"][0] + assert strong["recovered_trials"] == trials + assert (strong["angular_frequency"]["std"] is None) == (trials == 1) + + +def test_failed_sweep_statistics_are_null_in_summary(): + def broken(t, y, k): + raise RuntimeError("no fit") + + sweep = validation_sweep(snr_db=(30.0,), trials=2, estimator=broken) + summary = build_summary([sweep]) + json.dumps(summary, allow_nan=False) + result = summary["results"][0] + assert result["residual_rms_over_sigma_median"] is None + assert result["modes"][0]["angular_frequency"]["bias"] is None + + +@pytest.mark.parametrize( + "change", [{"samples": 192}, {"seed": 77}, {"match_tolerance": 0.1}] +) +def test_summary_rejects_incompatible_sweep_provenance(change): + sweep = validation_sweep(snr_db=(30.0,), trials=2) + with pytest.raises(ValueError, match="summary sweeps must share"): + build_summary([sweep, replace(sweep, **change)]) + + +def test_summary_preserves_distortion_signal_scale(): + sweep = validation_sweep(snr_db=(30.0,), trials=2) + distorted = replace( + sweep, backend="other", distortion=("clip", 0.5), signal_rms=0.1 + ) + summary = build_summary([sweep, distorted]) + assert summary["results"][0]["signal_rms"] == sweep.signal_rms + assert summary["results"][1]["signal_rms"] == 0.1 + + def test_invalid_inputs_are_rejected(): fx = synthetic_modes_trace(96) with pytest.raises(ValueError): From ace536d33779252fc7850cbe399e2a1abfbb3df1 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Wed, 23 Sep 2026 19:43:50 -0700 Subject: [PATCH 03/10] Normalize validation truth before comparing fitted modes Use the same amplitude and phase convention for truth, fitted modes, and summary metadata so equivalent modes have zero parameter error. Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-VALIDATION.md | 5 +-- src/cuphoton/xray/synthetic_validation.py | 19 ++++++++-- tests/xray/test_synthetic_validation.py | 44 +++++++++++++++++++++++ 3 files changed, 63 insertions(+), 5 deletions(-) diff --git a/docs/xray/LINEAR-PREDICTION-VALIDATION.md b/docs/xray/LINEAR-PREDICTION-VALIDATION.md index 6470528b..6cbb876c 100644 --- a/docs/xray/LINEAR-PREDICTION-VALIDATION.md +++ b/docs/xray/LINEAR-PREDICTION-VALIDATION.md @@ -63,8 +63,9 @@ mode was recovered; a trial in which the estimator raised counts in For each mode and parameter (`amplitude`, `decay`, `angular_frequency`, `phase`) the summary gives `bias`, `std` (ddof 1) and `rmse` computed over the recovered trials only, the bound (`crlb_std`, `crlb_variance`) and -`std_over_crlb_std`. A negative fitted amplitude is folded into the phase -and the phase error is wrapped to `[-pi, pi)`. +`std_over_crlb_std`. Truth and fitted modes use nonnegative amplitudes, +folding a negative amplitude into the phase. Both phases and phase errors +are wrapped to `[-pi, pi)`; the summary records truth in this convention. Undefined statistics are written as JSON `null`: standard deviation requires at least two recovered trials, and bias and RMSE require at least one. diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index 5fb62b96..0e41aab0 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -18,9 +18,10 @@ import json import subprocess -from dataclasses import asdict, dataclass, field +from collections.abc import Callable, Sequence +from dataclasses import asdict, dataclass, field, replace from pathlib import Path -from typing import Any, Callable, Sequence +from typing import Any import numpy as np from scipy.optimize import linear_sum_assignment @@ -250,7 +251,7 @@ def match_modes( err <= tolerance, err / (tolerance or 1.0), np.inf ) rows, columns = linear_sum_assignment(cost) - for row, column in zip(rows, columns): + for row, column in zip(rows, columns, strict=True): if column < angular_frequency.size: out[int(row)] = int(column) return out @@ -403,6 +404,10 @@ def validation_sweep( trial in which the estimator raises is counted in ``estimator_errors`` and as a loss of every mode. + Truth and fitted modes use nonnegative amplitudes and phases in + ``[-pi, pi)``. The returned truth uses the same convention as the + parameter statistics. + ``distortion`` = (kind, amount) replaces the clean trace with one from :func:`distort_trace`, outside the model class; ``residual_ratio`` is then the median rms residual of the reconstruction divided by the @@ -413,6 +418,14 @@ def validation_sweep( raise ValueError("trials must be at least 1") if not snr_db: raise ValueError("snr_db must contain at least one level") + amplitudes, phases = _canonical( + np.array([m.amplitude for m in modes]), + np.array([m.phase for m in modes]), + ) + modes = tuple( + replace(m, amplitude=float(a), phase=float(p)) + for m, a, p in zip(modes, amplitudes, phases, strict=True) + ) est = estimator or (lambda t, y, k: linear_prediction_numpy(t, y, k)) base = synthetic_modes_trace( samples, modes=modes, constant=constant, duration=duration diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py index 44eef03e..02943fbf 100644 --- a/tests/xray/test_synthetic_validation.py +++ b/tests/xray/test_synthetic_validation.py @@ -4,6 +4,7 @@ import json from dataclasses import replace +from types import SimpleNamespace import numpy as np import pytest @@ -161,6 +162,49 @@ def broken(t, y, k): assert np.isnan(level.residual_ratio) +@pytest.mark.parametrize("amplitude", [-1.0, 1.0]) +@pytest.mark.parametrize("phase", [0.3, 0.3 + 4 * np.pi]) +def test_sweep_uses_the_same_amplitude_phase_convention_for_truth_and_fit( + amplitude, phase +): + mode = DampedMode(amplitude, 0.1, 2.0, phase) + fixture = synthetic_modes_trace(modes=(mode,)) + + def exact_estimator(t, y, k): + return SimpleNamespace( + amplitude=np.array([mode.amplitude]), + decay=np.array([mode.decay]), + angular_frequency=np.array([mode.angular_frequency]), + phase=np.array([mode.phase]), + reconstruction=fixture.clean, + ) + + sweep = validation_sweep( + modes=(mode,), + snr_db=(30.0,), + trials=2, + n_components=1, + estimator=exact_estimator, + ) + result = sweep.levels[0].modes[0] + assert result.recovered_trials == 2 + for name in PARAMETERS: + stats = getattr(result, name) + assert stats.bias == pytest.approx(0.0, abs=1e-14) + assert stats.rmse == pytest.approx(0.0, abs=1e-14) + assert stats.truth == getattr(sweep.modes[0], name) + + canonical = synthetic_modes_trace(modes=sweep.modes) + np.testing.assert_allclose(canonical.clean, fixture.clean, atol=1e-14) + summary = build_summary([sweep]) + truth = summary["config"]["truth"]["modes"][0] + assert truth["amplitude"] == 1.0 + expected_phase = 0.3 - np.pi if amplitude < 0 else 0.3 + assert truth["phase"] == pytest.approx(expected_phase) + for name in PARAMETERS: + assert summary["results"][0]["modes"][0][name]["truth"] == truth[name] + + def test_summary_schema_and_run_artifacts(tmp_path): sweep = validation_sweep(snr_db=(30.0,), trials=6, seed=3) summary = build_summary([sweep]) From a39a14f7811bf7c989b1f2bdca65503ab864e3e0 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Wed, 23 Sep 2026 19:44:07 -0700 Subject: [PATCH 04/10] Move XRay command print helpers into a separate module Keep text formatting separate from command setup and execution. Signed-off-by: Trent Nelson --- src/cuphoton/xray/_cli_output.py | 133 ++++++++++++++++++++++++++++ src/cuphoton/xray/commands.py | 145 +++---------------------------- 2 files changed, 146 insertions(+), 132 deletions(-) create mode 100644 src/cuphoton/xray/_cli_output.py diff --git a/src/cuphoton/xray/_cli_output.py b/src/cuphoton/xray/_cli_output.py new file mode 100644 index 00000000..90f0276e --- /dev/null +++ b/src/cuphoton/xray/_cli_output.py @@ -0,0 +1,133 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +"""Text output helpers for XRay commands.""" + +from __future__ import annotations + + +def print_file_probe(label, probe): + print(f"{label}_path={probe.path}") + print(f"{label}_size_bytes={probe.size_bytes}") + print(f"{label}_schema={probe.schema}") + print(f"{label}_ipm_pairs={','.join(probe.ipm_pairs) or '-'}") + print(f"{label}_keys={','.join(probe.keys)}") + for dataset in probe.datasets: + shape = "x".join(str(dim) for dim in dataset.shape) + chunks = ( + "-" + if dataset.chunks is None + else "x".join(str(dim) for dim in dataset.chunks) + ) + print( + f"{label}_dataset={dataset.name} " + f"shape={shape} dtype={dataset.dtype} chunks={chunks}" + ) + + +def print_validation_sweep(sweep): + print(f"estimator={sweep.backend}") + if sweep.distortion is not None: + print(f"distortion={sweep.distortion[0]}:{sweep.distortion[1]:g}") + for level in sweep.levels: + print( + f"snr_db={level.snr_db:g} sigma={level.noise_sigma:.4g} " + f"trials={level.trials_successful}/{level.trials_attempted} " + f"any_mode_lost_rate={level.any_mode_lost_rate:.3f} " + f"residual_ratio={level.residual_ratio:.2f}" + ) + for m in level.modes: + w = m.angular_frequency + d = m.decay + print( + f" mode w={w.truth:g} decay={d.truth:g}: " + f"freq std/crlb={w.std:.3g}/{w.crlb_std:.3g} " + f"({w.std_over_crlb_std:.2f}x) bias={w.bias:+.3g}; " + f"decay std/crlb={d.std:.3g}/{d.crlb_std:.3g} " + f"({d.std_over_crlb_std:.2f}x) bias={d.bias:+.3g}; " + f"recovered={m.recovered_trials} " + f"loss_rate={m.loss_rate:.3f}" + ) + + +def print_model_order_sweep_payload(payload): + print(f"source={payload['source']['kind']}") + print(f"samples={payload['samples']}") + print(f"roots_backend={payload['roots_backend']}") + print(f"best_components={payload['best_components']}") + print(f"best_selected_model_order={payload['best_selected_model_order']}") + print(f"best_rms_residual={payload['best_rms_residual']:.6g}") + print( + "best_reconstruction_rms_error=" + f"{payload['best_reconstruction_rms_error']:.6g}" + ) + for entry in payload["entries"]: + print( + f"components={entry['components']} " + f"selected_model_order={entry['selected_model_order']} " + f"rms_residual={entry['rms_residual']:.6g} " + "reconstruction_rms_error=" + f"{entry['reconstruction_rms_error']:.6g} " + f"chi2={entry['chi2']:.6g} " + f"selected_roots={entry['selected_root_count']} " + f"decaying_roots={entry['decaying_root_count']}" + ) + + +def print_model_order_sweep_batch(payload): + print(f"source={payload['source']['kind']}") + print(f"trace_count={payload['trace_count']}") + print(f"samples={payload['samples']}") + print(f"roots_backend={payload['roots_backend']}") + component_counts = ",".join( + str(item) for item in payload["component_counts"] + ) + print(f"component_counts={component_counts}") + print( + "best_components_unique=" + + ",".join(str(item) for item in payload["best_components_unique"]) + ) + print( + "best_selected_model_orders_unique=" + + ",".join( + str(item) for item in payload["best_selected_model_orders_unique"] + ) + ) + for trace_payload in payload["traces"]: + print(f"trace_index={trace_payload['trace_index']}") + print(f" best_components={trace_payload['best_components']}") + print( + " best_selected_model_order=" + f"{trace_payload['best_selected_model_order']}" + ) + print( + " best_reconstruction_rms_error=" + f"{trace_payload['best_reconstruction_rms_error']:.6g}" + ) + + +def print_subspace_benchmark_batch(payload): + print(f"source={payload['source']['kind']}") + print(f"trace_count={payload['trace_count']}") + print(f"samples={payload['samples']}") + print(f"model_order={payload['model_order']}") + print(f"components={payload['components']}") + print( + "baseline_rms_residual_range=" + f"{payload['baseline_rms_residual_min']:.6g}:" + f"{payload['baseline_rms_residual_max']:.6g}" + ) + for summary in payload["method_summary"]: + print(f"method={summary['method']}") + print(f" svd_backend={summary['svd_backend']}") + print(f" trace_count={summary['trace_count']}") + print( + " rms_residual_range=" + f"{summary['rms_residual_min']:.6g}:" + f"{summary['rms_residual_max']:.6g}" + ) + print( + " max_abs_reconstruction_diff_max=" + f"{summary['max_abs_reconstruction_diff_max']:.6g}" + ) diff --git a/src/cuphoton/xray/commands.py b/src/cuphoton/xray/commands.py index a09273cb..534e01b1 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -27,6 +27,13 @@ VariablePositionalInvariant, ) +from ._cli_output import ( + print_file_probe, + print_model_order_sweep_batch, + print_model_order_sweep_payload, + print_subspace_benchmark_batch, + print_validation_sweep, +) from ._types import FIT_DIAGNOSTICS_LEVELS from .doctor import ( collect_doctor_report, @@ -3766,8 +3773,8 @@ def _data_probe(args): if args.json: print(json.dumps(payload, indent=2, sort_keys=True)) else: - _print_file_probe("on", result.on) - _print_file_probe("off", result.off) + print_file_probe("on", result.on) + print_file_probe("off", result.off) return 0 @@ -4319,25 +4326,6 @@ def _detector_mask(args): return 0 -def _print_file_probe(label, probe): - print(f"{label}_path={probe.path}") - print(f"{label}_size_bytes={probe.size_bytes}") - print(f"{label}_schema={probe.schema}") - print(f"{label}_ipm_pairs={','.join(probe.ipm_pairs) or '-'}") - print(f"{label}_keys={','.join(probe.keys)}") - for dataset in probe.datasets: - shape = "x".join(str(dim) for dim in dataset.shape) - chunks = ( - "-" - if dataset.chunks is None - else "x".join(str(dim) for dim in dataset.chunks) - ) - print( - f"{label}_dataset={dataset.name} " - f"shape={shape} dtype={dataset.dtype} chunks={chunks}" - ) - - def _linear_prediction_smoke(args): from .linear_prediction import ( LinearPredictionComparison, @@ -4420,31 +4408,6 @@ def _linear_prediction_smoke(args): return 0 -def _print_validation_sweep(sweep): - print(f"estimator={sweep.backend}") - if sweep.distortion is not None: - print(f"distortion={sweep.distortion[0]}:{sweep.distortion[1]:g}") - for level in sweep.levels: - print( - f"snr_db={level.snr_db:g} sigma={level.noise_sigma:.4g} " - f"trials={level.trials_successful}/{level.trials_attempted} " - f"any_mode_lost_rate={level.any_mode_lost_rate:.3f} " - f"residual_ratio={level.residual_ratio:.2f}" - ) - for m in level.modes: - w = m.angular_frequency - d = m.decay - print( - f" mode w={w.truth:g} decay={d.truth:g}: " - f"freq std/crlb={w.std:.3g}/{w.crlb_std:.3g} " - f"({w.std_over_crlb_std:.2f}x) bias={w.bias:+.3g}; " - f"decay std/crlb={d.std:.3g}/{d.crlb_std:.3g} " - f"({d.std_over_crlb_std:.2f}x) bias={d.bias:+.3g}; " - f"recovered={m.recovered_trials} " - f"loss_rate={m.loss_rate:.3f}" - ) - - def _linear_prediction_validate(args): from .synthetic_validation import ( build_summary, @@ -4477,7 +4440,7 @@ def _linear_prediction_validate(args): print(f"samples={sweeps[0].samples}") print(f"signal_rms={sweeps[0].signal_rms:.6g}") for sweep in sweeps: - _print_validation_sweep(sweep) + print_validation_sweep(sweep) if args.output_dir is not None: for key, name in summary["artifacts"].items(): print(f"{key}={name if name else 'not written'}") @@ -5511,9 +5474,9 @@ def _model_order_sweep(args): if args.json: print(json.dumps(payload, indent=2, sort_keys=True)) elif payload["source"]["kind"] == "trace-npz-batch": - _print_model_order_sweep_batch(payload) + print_model_order_sweep_batch(payload) else: - _print_model_order_sweep_payload(payload) + print_model_order_sweep_payload(payload) return 0 @@ -5591,62 +5554,6 @@ def _model_order_sweep_batch_payload( } -def _print_model_order_sweep_payload(payload): - print(f"source={payload['source']['kind']}") - print(f"samples={payload['samples']}") - print(f"roots_backend={payload['roots_backend']}") - print(f"best_components={payload['best_components']}") - print(f"best_selected_model_order={payload['best_selected_model_order']}") - print(f"best_rms_residual={payload['best_rms_residual']:.6g}") - print( - "best_reconstruction_rms_error=" - f"{payload['best_reconstruction_rms_error']:.6g}" - ) - for entry in payload["entries"]: - print( - f"components={entry['components']} " - f"selected_model_order={entry['selected_model_order']} " - f"rms_residual={entry['rms_residual']:.6g} " - "reconstruction_rms_error=" - f"{entry['reconstruction_rms_error']:.6g} " - f"chi2={entry['chi2']:.6g} " - f"selected_roots={entry['selected_root_count']} " - f"decaying_roots={entry['decaying_root_count']}" - ) - - -def _print_model_order_sweep_batch(payload): - print(f"source={payload['source']['kind']}") - print(f"trace_count={payload['trace_count']}") - print(f"samples={payload['samples']}") - print(f"roots_backend={payload['roots_backend']}") - component_counts = ",".join( - str(item) for item in payload["component_counts"] - ) - print(f"component_counts={component_counts}") - print( - "best_components_unique=" - + ",".join(str(item) for item in payload["best_components_unique"]) - ) - print( - "best_selected_model_orders_unique=" - + ",".join( - str(item) for item in payload["best_selected_model_orders_unique"] - ) - ) - for trace_payload in payload["traces"]: - print(f"trace_index={trace_payload['trace_index']}") - print(f" best_components={trace_payload['best_components']}") - print( - " best_selected_model_order=" - f"{trace_payload['best_selected_model_order']}" - ) - print( - " best_reconstruction_rms_error=" - f"{trace_payload['best_reconstruction_rms_error']:.6g}" - ) - - def _subspace_benchmark(args): from .subspace import compare_subspace_methods @@ -5714,7 +5621,7 @@ def _subspace_benchmark(args): if args.json: print(json.dumps(payload, indent=2, sort_keys=True)) elif payload["source"]["kind"] == "trace-npz-batch": - _print_subspace_benchmark_batch(payload) + print_subspace_benchmark_batch(payload) else: print(f"source={payload['source']['kind']}") print(f"samples={payload['samples']}") @@ -5844,32 +5751,6 @@ def _subspace_benchmark_batch_payload( } -def _print_subspace_benchmark_batch(payload): - print(f"source={payload['source']['kind']}") - print(f"trace_count={payload['trace_count']}") - print(f"samples={payload['samples']}") - print(f"model_order={payload['model_order']}") - print(f"components={payload['components']}") - print( - "baseline_rms_residual_range=" - f"{payload['baseline_rms_residual_min']:.6g}:" - f"{payload['baseline_rms_residual_max']:.6g}" - ) - for summary in payload["method_summary"]: - print(f"method={summary['method']}") - print(f" svd_backend={summary['svd_backend']}") - print(f" trace_count={summary['trace_count']}") - print( - " rms_residual_range=" - f"{summary['rms_residual_min']:.6g}:" - f"{summary['rms_residual_max']:.6g}" - ) - print( - " max_abs_reconstruction_diff_max=" - f"{summary['max_abs_reconstruction_diff_max']:.6g}" - ) - - def _subspace_acceptance(args): from .subspace import compare_subspace_methods From 29a03df798175ab8bb0519691fcf88b734bd618e Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Wed, 23 Sep 2026 20:19:37 -0700 Subject: [PATCH 05/10] Harden synthetic validation inputs and reporting Reject invalid numerical inputs before trials and distinguish estimator exceptions from mode loss in the text report. Record revisions only for tracked source checkouts, together with tracked- file dirtiness. Embed plot resources for offline use without changing Bokeh output state. Clarify numerical Fisher bounds and conditional scatter, and anchor the Gaussian distortion envelope at the first sample. Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-VALIDATION.md | 22 +++- src/cuphoton/xray/_cli_output.py | 1 + src/cuphoton/xray/commands.py | 17 ++- src/cuphoton/xray/synthetic_validation.py | 115 +++++++++++++----- tests/xray/test_cli.py | 45 +++++++ tests/xray/test_synthetic_validation.py | 65 +++++++++- .../test_synthetic_validation_provenance.py | 115 ++++++++++++++++++ tests/xray/test_synthetic_validation_viz.py | 79 ++++++++++++ 8 files changed, 420 insertions(+), 39 deletions(-) create mode 100644 tests/xray/test_synthetic_validation_provenance.py create mode 100644 tests/xray/test_synthetic_validation_viz.py diff --git a/docs/xray/LINEAR-PREDICTION-VALIDATION.md b/docs/xray/LINEAR-PREDICTION-VALIDATION.md index 6cbb876c..e81f3554 100644 --- a/docs/xray/LINEAR-PREDICTION-VALIDATION.md +++ b/docs/xray/LINEAR-PREDICTION-VALIDATION.md @@ -14,7 +14,8 @@ uv run cuphoton xray lpv --json ``` Keep `--output-dir` outside the checkout. It receives `summary.json` and, when -the `viz` extra is installed, `validation.html`. +the `viz` extra is installed, `validation.html`. The HTML includes its Bokeh +resources and can be viewed offline. ## Model and fixture @@ -37,6 +38,8 @@ Noise is white Gaussian, independent per sample. A level's `snr_db` is `20 log10(rms of the mean-removed clean trace / sigma)`. One `numpy` generator seeded by `--seed` draws every trial in order, so a run is reproducible from the summary's `config`. +SNR levels must be finite and produce finite, positive noise scales; +`--components` and `--trials` must be positive integers. ## Bounds @@ -58,7 +61,8 @@ their total frequency error. Each fitted mode is used at most once. The per-mode `loss_rate` is the fraction of trials in which the mode was not recovered. A trial is successful when the estimator returned and every true mode was recovered; a trial in which the estimator raised counts in -`estimator_errors` and as a loss of every mode. +`estimator_errors` and as a loss of every mode. The text report includes +the error count so estimator failures can be distinguished from mode loss. For each mode and parameter (`amplitude`, `decay`, `angular_frequency`, `phase`) the summary gives `bias`, `std` (ddof 1) and `rmse` computed over @@ -74,7 +78,8 @@ Read the decay statistics of a lightly damped mode together with its the fitted decay of such a mode through zero the mode is dropped rather than reported with a negative decay; the surviving decay estimates are truncated at zero, and their scatter and bias are conditional on survival. The summary -labels these statistics `statistics_over: recovered_trials`. +labels these statistics `statistics_over: recovered_trials`. Conditional +statistics and finite-sample scatter can fall below the unconditional bound. Every level also reports `residual_rms_over_sigma_median`, the median over trials of the rms residual of the reconstruction divided by sigma. It is @@ -95,6 +100,9 @@ model class: | `clip` | the trace is clipped above `amount` | | `glitch` | one sample at mid-record is offset by `amount` | +Amounts must be finite. The Gaussian 1/e time must be positive and is measured +from the first sample. + The truth used for bias and bounds stays the undistorted mode set, so the statistics measure the estimator's response to the mismatch. A chirp or a clipped trace can return modes with a confident, biased frequency and no @@ -112,8 +120,12 @@ loss; the residual ratio supplies an additional reconstruction diagnostic. the estimators and distortions covered. - `runtime`: `cuphoton.core.runtime.runtime_metadata` for the CPU backend (package, Python, platform, NumPy versions, backend, device, dtype) plus - `scipy_version` and `source_revision` (the git revision of the checkout, - or `null`). + `scipy_version`, `source_revision`, and `source_dirty`. A revision is + recorded only when this module is tracked by its Git checkout; + `source_dirty` reports tracked changes anywhere in that checkout. + Both are `null` for an installed or untracked module, or when Git metadata + is unavailable. Ignored and untracked run artifacts do not mark the source + dirty. - `results`: one entry per sweep condition (estimator by distortion by signal-to-noise level): `snr_db`, `signal_rms`, `noise_sigma`, `trials` with `attempted`, `successful`, `failed` and `estimator_errors`, diff --git a/src/cuphoton/xray/_cli_output.py b/src/cuphoton/xray/_cli_output.py index 90f0276e..7bfaa973 100644 --- a/src/cuphoton/xray/_cli_output.py +++ b/src/cuphoton/xray/_cli_output.py @@ -34,6 +34,7 @@ def print_validation_sweep(sweep): print( f"snr_db={level.snr_db:g} sigma={level.noise_sigma:.4g} " f"trials={level.trials_successful}/{level.trials_attempted} " + f"estimator_errors={level.estimator_errors} " f"any_mode_lost_rate={level.any_mode_lost_rate:.3f} " f"residual_ratio={level.residual_ratio:.2f}" ) diff --git a/src/cuphoton/xray/commands.py b/src/cuphoton/xray/commands.py index 534e01b1..c97852b0 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -21,6 +21,7 @@ NonNegativeIntegerInvariant, PairInvariant, PathValueInvariant, + PositiveIntegerInvariant, SequenceInvariant, SetInvariant, StringInvariant, @@ -621,7 +622,7 @@ class SamplesArg(IntegerInvariant): components = 6 - class ComponentsArg(IntegerInvariant): + class ComponentsArg(PositiveIntegerInvariant): _arg = "--components" _help = "Number of SVD components to fit." _required = False @@ -629,7 +630,7 @@ class ComponentsArg(IntegerInvariant): trials = 200 - class TrialsArg(IntegerInvariant): + class TrialsArg(PositiveIntegerInvariant): _arg = "--trials" _help = "Monte Carlo trials per signal-to-noise level." _required = False @@ -4418,8 +4419,16 @@ def _linear_prediction_validate(args): snr = tuple(float(x) for x in str(args.snr_db).split(",") if x.strip()) distortion = None if args.distortion: - kind, _, amount = str(args.distortion).partition(":") - distortion = (kind.strip(), float(amount)) + try: + kind, amount = str(args.distortion).split(":", maxsplit=1) + if not kind.strip(): + raise ValueError("missing distortion kind") + distortion = (kind.strip(), float(amount)) + except ValueError as exc: + raise ValueError( + "--distortion must be kind:amount with a numeric amount " + "(for example chirp:0.05)" + ) from exc sweeps = [ validation_sweep( samples=args.samples, diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index 0e41aab0..b4a45d4d 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -8,7 +8,7 @@ The fixture generates traces with known modes so the workflow can be validated without external datasets. The validation sweep compares the estimator's scatter and bias with the Cramer-Rao lower bound computed from -the exact Fisher information of the same model, together with the rate at +the Fisher information of the same model, together with the rate at which true modes are lost from the fit. Everything here runs on NumPy; the estimator under test is :func:`linear_prediction_numpy` unless another estimator callable is supplied. @@ -108,6 +108,10 @@ def synthetic_modes_trace( """ if samples < 16: raise ValueError("samples must be at least 16") + if not np.isfinite(duration) or duration <= 0: + raise ValueError("duration must be finite and positive") + if not np.isfinite(noise_sigma) or noise_sigma < 0: + raise ValueError("noise_sigma must be finite and nonnegative") time = np.linspace(0.0, duration, samples, dtype=np.float64) clean = modes_model(modes_to_theta(modes, constant), time, len(modes)) trace = clean @@ -139,7 +143,8 @@ def distort_trace( ``chirp``: every frequency rises by the fraction ``amount`` across the record. ``gaussian_envelope``: the first mode decays as a Gaussian with - 1/e time ``amount`` instead of an exponential. ``baseline_drift``: a + 1/e time ``amount`` from the first sample instead of an exponential. + ``baseline_drift``: a linear ramp of ``amount`` across the record. ``clip``: the trace is clipped above ``amount``. ``glitch``: one sample at mid-record is offset by ``amount``. These are outside the model class on purpose: @@ -148,6 +153,10 @@ def distort_trace( """ if kind not in DISTORTIONS: raise ValueError(f"unknown distortion {kind!r}") + if not np.isfinite(amount): + raise ValueError("distortion amount must be finite") + if kind == "gaussian_envelope" and amount <= 0: + raise ValueError("gaussian_envelope amount must be positive") span = float(time[-1] - time[0]) out = np.full(time.shape, float(constant)) for k, m in enumerate(modes): @@ -161,7 +170,8 @@ def distort_trace( / span ) if kind == "gaussian_envelope" and k == 0: - env = np.exp(-((time / amount) ** 2)) + with np.errstate(over="ignore", under="ignore"): + env = np.exp(-(((time - time[0]) / amount) ** 2)) else: env = np.exp(-m.decay * time) out = out + m.amplitude * env * np.cos(phase) @@ -190,8 +200,14 @@ def cramer_rao_bounds( its inverse. Returned per mode as arrays ``amplitude``, ``decay``, ``angular_frequency``, ``phase`` plus the scalar ``constant``. """ - if noise_sigma <= 0: - raise ValueError("noise_sigma must be positive") + if not np.isfinite(noise_sigma) or noise_sigma <= 0: + raise ValueError("noise_sigma must be finite and positive") + with np.errstate(over="ignore", under="ignore"): + noise_variance = np.square(np.float64(noise_sigma)) + if not np.isfinite(noise_variance) or noise_variance <= 0: + raise ValueError("noise_sigma must have a finite, positive variance") + if not np.isfinite(step) or step <= 0: + raise ValueError("step must be finite and positive") theta = modes_to_theta(modes, constant) n = len(modes) jac = np.empty((time.size, theta.size)) @@ -203,8 +219,16 @@ def cramer_rao_bounds( jac[:, k] = (modes_model(tp, time, n) - modes_model(tm, time, n)) / ( 2 * step ) - fisher = jac.T @ jac / noise_sigma**2 - sd = np.sqrt(np.diag(np.linalg.inv(fisher))) + try: + with np.errstate(over="raise", divide="raise", invalid="raise"): + fisher = jac.T @ jac / noise_variance + sd = np.sqrt(np.diag(np.linalg.inv(fisher))) + except FloatingPointError as exc: + raise ValueError( + "noise_sigma produces nonfinite Cramer-Rao bounds" + ) from exc + if not np.all(np.isfinite(sd)) or np.any(sd <= 0): + raise ValueError("Cramer-Rao bounds must be finite and positive") return { "amplitude": sd[0 : 4 * n : 4], "decay": sd[1 : 4 * n : 4], @@ -416,8 +440,14 @@ def validation_sweep( """ if trials < 1: raise ValueError("trials must be at least 1") - if not snr_db: + if n_components <= 0: + raise ValueError("n_components must be positive") + if not np.isfinite(match_tolerance) or match_tolerance < 0: + raise ValueError("match_tolerance must be finite and nonnegative") + if len(snr_db) == 0: raise ValueError("snr_db must contain at least one level") + if not np.all(np.isfinite(snr_db)): + raise ValueError("snr_db levels must be finite") amplitudes, phases = _canonical( np.array([m.amplitude for m in modes]), np.array([m.phase for m in modes]), @@ -426,7 +456,7 @@ def validation_sweep( replace(m, amplitude=float(a), phase=float(p)) for m, a, p in zip(modes, amplitudes, phases, strict=True) ) - est = estimator or (lambda t, y, k: linear_prediction_numpy(t, y, k)) + est = estimator or linear_prediction_numpy base = synthetic_modes_trace( samples, modes=modes, constant=constant, duration=duration ) @@ -436,11 +466,17 @@ def validation_sweep( base.time, modes, constant, distortion[0], distortion[1] ) signal_rms = float(np.std(clean - clean.mean())) + with np.errstate(over="ignore", under="ignore", divide="ignore"): + noise_sigmas = signal_rms / np.power(10.0, np.asarray(snr_db) / 20) + # Validate every level before trials, including underflow/overflow in + # the SNR conversion, so a later invalid level cannot yield partial work. + bounds_per_level = [ + cramer_rao_bounds(modes, constant, base.time, float(sigma)) + for sigma in noise_sigmas + ] rng = np.random.default_rng(seed) levels = [] - for snr in snr_db: - sigma = signal_rms / 10 ** (snr / 20) - bounds = cramer_rao_bounds(modes, constant, base.time, sigma) + for snr, sigma, bounds in zip(snr_db, noise_sigmas, bounds_per_level): found: list[dict[str, list[float]]] = [ {name: [] for name in PARAMETERS} for _ in modes ] @@ -524,19 +560,38 @@ def validation_sweep( ) -def source_revision() -> str | None: - """Git revision of the checkout containing this module, if any.""" +def _source_provenance() -> tuple[str | None, bool | None]: + """Revision and tracked-file dirtiness for this module's checkout.""" try: - out = subprocess.run( - ["git", "rev-parse", "HEAD"], - cwd=Path(__file__).resolve().parent, - capture_output=True, - text=True, - timeout=5, + source = Path(__file__).resolve() + options = { + "cwd": source.parent, + "capture_output": True, + "text": True, + "timeout": 5, + "check": True, + } + # A wheel installed in another project's .venv can have an enclosing + # Git checkout without belonging to that checkout's source tree. + subprocess.run( + ["git", "ls-files", "--error-unmatch", "--", source.name], + **options, + ) + revision = subprocess.run( + ["git", "rev-parse", "HEAD"], **options + ).stdout.strip() + status = subprocess.run( + ["git", "status", "--porcelain", "--untracked-files=no"], + **options, ) except (OSError, subprocess.SubprocessError): - return None - return out.stdout.strip() if out.returncode == 0 else None + return None, None + return revision, bool(status.stdout.strip()) + + +def source_revision() -> str | None: + """Git revision when this module belongs to a tracked source checkout.""" + return _source_provenance()[0] def build_summary( @@ -631,8 +686,9 @@ def build_summary( "statistics": ( "bias, std (ddof 1) and rmse are computed over the recovered " "trials of each mode; crlb_std and crlb_variance are the " - "unconditional Cramer-Rao bounds from the exact Fisher " - "information of the model at the true parameters; the phase " + "unconditional Cramer-Rao bounds from the model's Fisher " + "information using a numerical Jacobian at the true " + "parameters; the phase " "error is wrapped to [-pi, pi)" ), "estimators": sorted({s.backend for s in sweeps}), @@ -651,7 +707,7 @@ def build_summary( except ImportError: # pragma: no cover - scipy is a dependency scipy_version = None runtime["scipy_version"] = scipy_version - runtime["source_revision"] = source_revision() + runtime["source_revision"], runtime["source_dirty"] = _source_provenance() results = [] for sweep in sweeps: for lvl in sweep.levels: @@ -717,8 +773,10 @@ def write_validation_figure( standard deviation against the bound versus signal-to-noise, and the loss rate. Returns False when the ``viz`` extra is not installed.""" try: + from bokeh.embed import file_html from bokeh.layouts import gridplot - from bokeh.plotting import figure, output_file, save + from bokeh.plotting import figure + from bokeh.resources import INLINE except ImportError: return False palette = ["#1f77b4", "#d62728", "#2ca02c", "#9467bd", "#ff7f0e"] @@ -783,8 +841,9 @@ def write_validation_figure( fig.legend.label_text_font_size = "8pt" row.append(fig) rows.append(row) - output_file(str(path), title=title) - save(gridplot(rows)) + path.write_text( + file_html(gridplot(rows), INLINE, title), encoding="utf-8" + ) return True diff --git a/tests/xray/test_cli.py b/tests/xray/test_cli.py index 7ea68a11..eec8d033 100644 --- a/tests/xray/test_cli.py +++ b/tests/xray/test_cli.py @@ -256,6 +256,51 @@ def test_linear_prediction_validate_text_json_and_output_dir( assert written["artifacts"]["summary"] == "summary.json" +@pytest.mark.parametrize("option", ["--components", "--trials"]) +@pytest.mark.parametrize("value", ["0", "-1"]) +def test_linear_prediction_validate_requires_positive_counts( + capsys, option, value +): + assert main(["lpv", option, value]) == 2 + captured = capsys.readouterr() + assert f"argument {option}: must be at least 1" in captured.err + assert not captured.out + + +@pytest.mark.parametrize("json_args", [[], ["--json"]]) +@pytest.mark.parametrize( + ("args", "message"), + [ + (["--snr-db", "nan"], "snr_db"), + (["--snr-db", "inf"], "snr_db"), + (["--distortion", "gaussian_envelope:0"], "gaussian_envelope"), + (["--distortion", "chirp"], "--distortion must be kind:amount"), + (["--distortion", "chirp:bad"], "--distortion must be kind:amount"), + ], +) +def test_linear_prediction_validate_rejects_invalid_reports( + capsys, args, message, json_args +): + _assert_cli_error(capsys, ["lpv", *args, *json_args], message) + + +def test_linear_prediction_validate_text_distinguishes_estimator_errors( + capsys, monkeypatch +): + from cuphoton.xray import synthetic_validation + + def broken(t, y, k): + raise RuntimeError("estimator failed") + + monkeypatch.setattr( + synthetic_validation, "linear_prediction_numpy", broken + ) + assert main(["lpv", "--trials", "2", "--snr-db", "30"]) == 0 + output = capsys.readouterr().out + assert "trials=0/2 estimator_errors=2" in output + assert "any_mode_lost_rate=1.000" in output + + def test_gpu_policy(capsys): assert main(["gpu-policy"]) == 0 captured = capsys.readouterr() diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py index 02943fbf..73d00647 100644 --- a/tests/xray/test_synthetic_validation.py +++ b/tests/xray/test_synthetic_validation.py @@ -110,6 +110,16 @@ def test_chirp_integrates_the_requested_frequency_ramp(start): np.testing.assert_allclose(actual, expected) +@pytest.mark.parametrize("start", [0.0, 3.0]) +def test_gaussian_envelope_starts_at_the_first_sample(start): + time = np.linspace(start, start + 10.0, 256) + expected = np.exp(-(((time - start) / 5.0) ** 2)) * np.cos(2 * time + 0.3) + actual = distort_trace( + time, (DampedMode(1, 0.1, 2, 0.3),), 0, "gaussian_envelope", 5.0 + ) + np.testing.assert_allclose(actual, expected) + + def test_matching_recovers_close_modes_independent_of_truth_order(): modes = (DampedMode(1, 0, 1.0), DampedMode(1, 0, 1.15)) fitted = np.array([1.05, 0.85]) @@ -132,8 +142,9 @@ def test_validation_sweep_reports_bound_ratios_and_loss_rates(): w = strong.angular_frequency assert w.crlb_std > 0 and w.std > 0 assert w.crlb_variance == pytest.approx(w.crlb_std**2) - assert w.std_over_crlb_std >= 1.0 # an estimator cannot beat the bound - assert w.std_over_crlb_std < 10.0 + # A coarse plausibility check: finite-sample, recovered-trial scatter + # from a biased estimator need not exceed the unconditional CRLB. + assert 0.2 < w.std_over_crlb_std < 10.0 assert abs(w.bias) < 0.01 assert w.rmse >= abs(w.bias) for name in PARAMETERS: @@ -300,6 +311,56 @@ def test_invalid_inputs_are_rejected(): synthetic_modes_trace(8) +@pytest.mark.parametrize("sigma", [-0.1, np.nan, np.inf, -np.inf]) +def test_fixture_rejects_invalid_noise_sigma(sigma): + with pytest.raises(ValueError, match="noise_sigma"): + synthetic_modes_trace(noise_sigma=sigma) + + +@pytest.mark.parametrize( + "sigma", [np.nan, np.inf, -0.1, 1e-200, 1e-160, 1e200] +) +def test_bounds_reject_invalid_noise_scales_without_runtime_warnings(sigma): + fx = synthetic_modes_trace() + with ( + np.errstate(all="raise"), + pytest.raises(ValueError, match="noise_sigma"), + ): + cramer_rao_bounds(fx.modes, fx.constant, fx.time, sigma) + + +@pytest.mark.parametrize( + "kwargs, message", + [ + ({"snr_db": (30.0, np.nan)}, "snr_db"), + ({"snr_db": (30.0, np.inf)}, "snr_db"), + ({"snr_db": (30.0, -np.inf)}, "snr_db"), + ({"snr_db": (30.0, 1e308)}, "noise_sigma"), + ({"snr_db": (30.0, -1e308)}, "noise_sigma"), + ({"snr_db": (30.0, 4000.0)}, "noise_sigma"), + ({"snr_db": (30.0, -4000.0)}, "noise_sigma"), + ({"n_components": 0}, "n_components"), + ({"distortion": ("gaussian_envelope", 0.0)}, "must be positive"), + ({"distortion": ("gaussian_envelope", -1.0)}, "must be positive"), + ({"distortion": ("glitch", np.inf)}, "must be finite"), + ({"match_tolerance": np.nan}, "match_tolerance"), + ({"duration": 0.0}, "duration"), + ], +) +def test_invalid_sweep_input_is_rejected_before_any_estimator_calls( + kwargs, message +): + calls = [] + + def estimator(*args): + calls.append(args) + raise RuntimeError("invalid input must not reach the estimator") + + with np.errstate(all="raise"), pytest.raises(ValueError, match=message): + validation_sweep(trials=1, estimator=estimator, **kwargs) + assert not calls + + def test_undistorted_residual_ratio_is_near_one(): sweep = validation_sweep(snr_db=(30.0,), trials=30, seed=9) assert 0.7 < sweep.levels[0].residual_ratio < 1.6 diff --git a/tests/xray/test_synthetic_validation_provenance.py b/tests/xray/test_synthetic_validation_provenance.py new file mode 100644 index 00000000..958416c6 --- /dev/null +++ b/tests/xray/test_synthetic_validation_provenance.py @@ -0,0 +1,115 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +import subprocess + +import pytest + +from cuphoton.xray import synthetic_validation + + +@pytest.fixture +def source_checkout(tmp_path, monkeypatch): + def git(*args): + return subprocess.run( + [ + "git", + "-c", + "user.name=cuPhoton tests", + "-c", + "user.email=test@nvidia.com", + "-c", + "commit.gpgsign=false", + "-c", + "core.hooksPath=/dev/null", + *args, + ], + cwd=tmp_path, + capture_output=True, + text=True, + check=True, + timeout=5, + ).stdout.strip() + + git("init", "-q") + source = tmp_path / "src/cuphoton/xray/synthetic_validation.py" + source.parent.mkdir(parents=True) + source.write_text("# source checkout\n") + (tmp_path / ".gitignore").write_text(".venv/\n") + git("add", "src", ".gitignore") + git("commit", "-qm", "Initial test source") + monkeypatch.setattr(synthetic_validation, "__file__", str(source)) + return source, git + + +def test_clean_tracked_source_reports_its_revision(source_checkout): + source, git = source_checkout + revision = git("rev-parse", "HEAD") + assert synthetic_validation._source_provenance() == (revision, False) + assert synthetic_validation.source_revision() == revision + # Generated run files do not change the tracked source provenance. + (source.parent / "summary.json").write_text("{}\n") + assert synthetic_validation._source_provenance() == (revision, False) + + +@pytest.mark.parametrize("staged", [False, True]) +def test_changed_tracked_source_is_dirty(source_checkout, staged): + source, git = source_checkout + source.write_text("# changed source\n") + if staged: + git("add", "src") + assert synthetic_validation._source_provenance() == ( + git("rev-parse", "HEAD"), + True, + ) + + +def test_dirty_tracks_changes_outside_the_module_directory(source_checkout): + source, git = source_checkout + checkout = source.parents[3] + (checkout / ".gitignore").write_text(".venv/\n*.json\n") + assert synthetic_validation._source_provenance() == ( + git("rev-parse", "HEAD"), + True, + ) + + +@pytest.mark.parametrize("directory", [".venv", "untracked"]) +def test_untracked_install_does_not_report_enclosing_checkout( + source_checkout, monkeypatch, directory +): + source, _ = source_checkout + installed = ( + source.parents[3] + / directory + / "lib/python3.12/site-packages/cuphoton/xray/synthetic_validation.py" + ) + installed.parent.mkdir(parents=True) + installed.write_text("# installed package\n") + monkeypatch.setattr(synthetic_validation, "__file__", str(installed)) + assert synthetic_validation._source_provenance() == (None, None) + assert synthetic_validation.source_revision() is None + + +def test_non_git_install_has_unknown_provenance(tmp_path, monkeypatch): + source = tmp_path / "synthetic_validation.py" + source.write_text("# installed package\n") + monkeypatch.setattr(synthetic_validation, "__file__", str(source)) + assert synthetic_validation._source_provenance() == (None, None) + + +def test_git_unavailable_has_unknown_provenance(source_checkout, monkeypatch): + monkeypatch.setenv("PATH", "") + assert synthetic_validation._source_provenance() == (None, None) + + +@pytest.mark.parametrize("dirty", [False, True]) +def test_summary_includes_tracked_source_provenance(source_checkout, dirty): + source, git = source_checkout + if dirty: + source.write_text("# changed source\n") + sweep = synthetic_validation.validation_sweep(snr_db=(30,), trials=2) + runtime = synthetic_validation.build_summary([sweep])["runtime"] + assert runtime["source_revision"] == git("rev-parse", "HEAD") + assert runtime["source_dirty"] is dirty diff --git a/tests/xray/test_synthetic_validation_viz.py b/tests/xray/test_synthetic_validation_viz.py new file mode 100644 index 00000000..fe7d0ca1 --- /dev/null +++ b/tests/xray/test_synthetic_validation_viz.py @@ -0,0 +1,79 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +import sys +from html.parser import HTMLParser + +import pytest + +from cuphoton.xray.synthetic_validation import ( + validation_sweep, + write_validation_figure, +) + + +class _ResourceParser(HTMLParser): + def __init__(self): + super().__init__() + self.resource_urls = [] + + def handle_starttag(self, tag, attrs): + resource_attr = {"script": "src", "link": "href"}.get(tag) + if resource_attr is not None: + url = dict(attrs).get(resource_attr) + if url is not None: + self.resource_urls.append(url) + + +@pytest.fixture +def sweep(): + return validation_sweep(snr_db=(30.0,), trials=2) + + +def test_validation_figure_embeds_resources_for_offline_use(tmp_path, sweep): + pytest.importorskip("bokeh") + path = tmp_path / "validation.html" + + assert write_validation_figure([sweep], path, title="Validation figure") + + html = path.read_text(encoding="utf-8") + resources = _ResourceParser() + resources.feed(html) + assert resources.resource_urls == [] + assert "/* BEGIN bokeh.min.js */" in html + assert "Validation figure" in html + assert "angular_frequency std (points) vs bound (dashed)" in html + assert "decay std (points) vs bound (dashed)" in html + assert "loss rate" in html + + +def test_validation_figure_preserves_bokeh_output_state( + tmp_path, sweep, monkeypatch +): + pytest.importorskip("bokeh") + from bokeh.io import state + + current = state.State() + monkeypatch.setattr(state, "_STATE", current) + current.output_file(tmp_path / "existing.html", title="Existing output") + existing_file = current.file + existing_document = current.document + + assert write_validation_figure( + [sweep], tmp_path / "validation.html", title="Validation figure" + ) + + assert current.file is existing_file + assert current.document is existing_document + assert current.document.roots == [] + + +def test_validation_figure_returns_false_without_bokeh(tmp_path, monkeypatch): + monkeypatch.setitem(sys.modules, "bokeh", None) + monkeypatch.setitem(sys.modules, "bokeh.embed", None) + monkeypatch.setitem(sys.modules, "bokeh.layouts", None) + path = tmp_path / "validation.html" + + assert not write_validation_figure([], path, title="Validation figure") + assert not path.exists() From 4fe90282f0f9d33147010e674d4eec207c04346c Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Thu, 24 Sep 2026 22:20:55 -0700 Subject: [PATCH 06/10] Reject singular Fisher information in XRay validation Report unidentifiable modes with the validation error context. Clarify that the frequency standard-deviation bound uses the square root of the closed-form variance. Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-VALIDATION.md | 9 +++++---- src/cuphoton/xray/synthetic_validation.py | 5 +++-- tests/xray/test_synthetic_validation.py | 13 +++++++++++++ 3 files changed, 21 insertions(+), 6 deletions(-) diff --git a/docs/xray/LINEAR-PREDICTION-VALIDATION.md b/docs/xray/LINEAR-PREDICTION-VALIDATION.md index e81f3554..947a5aed 100644 --- a/docs/xray/LINEAR-PREDICTION-VALIDATION.md +++ b/docs/xray/LINEAR-PREDICTION-VALIDATION.md @@ -47,10 +47,11 @@ SNR levels must be finite and produce finite, positive noise scales; numerical Jacobian of the model at the true parameters (all parameters free, including the constant) and returns the square root of the diagonal of its inverse per parameter. These are standard-deviation bounds; the summary -also carries the squared value as `crlb_variance`. For a single undamped -sinusoid the angular-frequency bound matches the closed form -`24 sigma^2 / (A^2 dt^2 N (N^2 - 1))` (Kay, Estimation Theory, 1993) to -within 3 percent, the difference being the freed decay and constant. +also carries the squared value as `crlb_variance`. For the single undamped +sinusoid test fixture, the angular-frequency bound matches the square root +of the closed-form variance `24 sigma^2 / (A^2 dt^2 N (N^2 - 1))` +(Kay, Estimation Theory, 1993) to within 3 percent. The numerical calculation +also frees the decay and constant. ## Mode matching and statistics diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index b4a45d4d..b149ebf7 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -223,9 +223,10 @@ def cramer_rao_bounds( with np.errstate(over="raise", divide="raise", invalid="raise"): fisher = jac.T @ jac / noise_variance sd = np.sqrt(np.diag(np.linalg.inv(fisher))) - except FloatingPointError as exc: + except (FloatingPointError, np.linalg.LinAlgError) as exc: raise ValueError( - "noise_sigma produces nonfinite Cramer-Rao bounds" + "noise_sigma or modes produce singular or nonfinite " + "Cramer-Rao bounds" ) from exc if not np.all(np.isfinite(sd)) or np.any(sd <= 0): raise ValueError("Cramer-Rao bounds must be finite and positive") diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py index 73d00647..ac647da2 100644 --- a/tests/xray/test_synthetic_validation.py +++ b/tests/xray/test_synthetic_validation.py @@ -101,6 +101,19 @@ def test_integer_mode_parameters_have_the_same_bounds_as_floats(): assert np.all(np.isfinite(actual[name])) +@pytest.mark.parametrize( + "modes", + [ + (DampedMode(0.0, 0.1, 2.0, 0.3),), + (DampedMode(1.0, 0.1, 2.0, 0.3),) * 2, + ], + ids=["zero-amplitude", "duplicate-modes"], +) +def test_bounds_reject_unidentifiable_modes(modes): + with pytest.raises(ValueError, match="Cramer-Rao bounds"): + cramer_rao_bounds(modes, 0.0, np.linspace(0.0, 9.5, 96), 0.01) + + @pytest.mark.parametrize("start", [0.0, 3.0]) def test_chirp_integrates_the_requested_frequency_ramp(start): time = np.linspace(start, start + 10.0, 256) From e9d9f6c9e845110e30161b875700b03f73bf852b Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Mon, 28 Sep 2026 14:48:36 -0700 Subject: [PATCH 07/10] Apply current lint rules to XRay validation Signed-off-by: Trent Nelson --- src/cuphoton/xray/synthetic_validation.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index b149ebf7..471eaf4e 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -477,7 +477,9 @@ def validation_sweep( ] rng = np.random.default_rng(seed) levels = [] - for snr, sigma, bounds in zip(snr_db, noise_sigmas, bounds_per_level): + for snr, sigma, bounds in zip( + snr_db, noise_sigmas, bounds_per_level, strict=True + ): found: list[dict[str, list[float]]] = [ {name: [] for name in PARAMETERS} for _ in modes ] @@ -565,7 +567,7 @@ def _source_provenance() -> tuple[str | None, bool | None]: """Revision and tracked-file dirtiness for this module's checkout.""" try: source = Path(__file__).resolve() - options = { + options: dict[str, Any] = { "cwd": source.parent, "capture_output": True, "text": True, From a4f74be5ea8ad1f428a35a6fc728f5ed1d5a8806 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Wed, 30 Sep 2026 19:10:51 -0700 Subject: [PATCH 08/10] Move XRay reports and validation into domain modules Give probe, validation, and benchmark results their own text formatting. Keep batch aggregation and validation orchestration in the owning modules, with commands routing reports through their injected output. Signed-off-by: Trent Nelson --- src/cuphoton/xray/_cli_output.py | 134 ------- src/cuphoton/xray/commands.py | 407 +++++----------------- src/cuphoton/xray/hdf5.py | 35 ++ src/cuphoton/xray/linear_prediction.py | 156 ++++++++- src/cuphoton/xray/subspace.py | 220 +++++++++++- src/cuphoton/xray/synthetic_validation.py | 91 +++++ tests/xray/test_cli.py | 47 +++ tests/xray/test_linear_prediction.py | 50 +++ tests/xray/test_subspace.py | 99 ++++++ tests/xray/test_synthetic_validation.py | 36 +- 10 files changed, 817 insertions(+), 458 deletions(-) delete mode 100644 src/cuphoton/xray/_cli_output.py diff --git a/src/cuphoton/xray/_cli_output.py b/src/cuphoton/xray/_cli_output.py deleted file mode 100644 index 7bfaa973..00000000 --- a/src/cuphoton/xray/_cli_output.py +++ /dev/null @@ -1,134 +0,0 @@ -# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. -# -# SPDX-License-Identifier: Apache-2.0 - -"""Text output helpers for XRay commands.""" - -from __future__ import annotations - - -def print_file_probe(label, probe): - print(f"{label}_path={probe.path}") - print(f"{label}_size_bytes={probe.size_bytes}") - print(f"{label}_schema={probe.schema}") - print(f"{label}_ipm_pairs={','.join(probe.ipm_pairs) or '-'}") - print(f"{label}_keys={','.join(probe.keys)}") - for dataset in probe.datasets: - shape = "x".join(str(dim) for dim in dataset.shape) - chunks = ( - "-" - if dataset.chunks is None - else "x".join(str(dim) for dim in dataset.chunks) - ) - print( - f"{label}_dataset={dataset.name} " - f"shape={shape} dtype={dataset.dtype} chunks={chunks}" - ) - - -def print_validation_sweep(sweep): - print(f"estimator={sweep.backend}") - if sweep.distortion is not None: - print(f"distortion={sweep.distortion[0]}:{sweep.distortion[1]:g}") - for level in sweep.levels: - print( - f"snr_db={level.snr_db:g} sigma={level.noise_sigma:.4g} " - f"trials={level.trials_successful}/{level.trials_attempted} " - f"estimator_errors={level.estimator_errors} " - f"any_mode_lost_rate={level.any_mode_lost_rate:.3f} " - f"residual_ratio={level.residual_ratio:.2f}" - ) - for m in level.modes: - w = m.angular_frequency - d = m.decay - print( - f" mode w={w.truth:g} decay={d.truth:g}: " - f"freq std/crlb={w.std:.3g}/{w.crlb_std:.3g} " - f"({w.std_over_crlb_std:.2f}x) bias={w.bias:+.3g}; " - f"decay std/crlb={d.std:.3g}/{d.crlb_std:.3g} " - f"({d.std_over_crlb_std:.2f}x) bias={d.bias:+.3g}; " - f"recovered={m.recovered_trials} " - f"loss_rate={m.loss_rate:.3f}" - ) - - -def print_model_order_sweep_payload(payload): - print(f"source={payload['source']['kind']}") - print(f"samples={payload['samples']}") - print(f"roots_backend={payload['roots_backend']}") - print(f"best_components={payload['best_components']}") - print(f"best_selected_model_order={payload['best_selected_model_order']}") - print(f"best_rms_residual={payload['best_rms_residual']:.6g}") - print( - "best_reconstruction_rms_error=" - f"{payload['best_reconstruction_rms_error']:.6g}" - ) - for entry in payload["entries"]: - print( - f"components={entry['components']} " - f"selected_model_order={entry['selected_model_order']} " - f"rms_residual={entry['rms_residual']:.6g} " - "reconstruction_rms_error=" - f"{entry['reconstruction_rms_error']:.6g} " - f"chi2={entry['chi2']:.6g} " - f"selected_roots={entry['selected_root_count']} " - f"decaying_roots={entry['decaying_root_count']}" - ) - - -def print_model_order_sweep_batch(payload): - print(f"source={payload['source']['kind']}") - print(f"trace_count={payload['trace_count']}") - print(f"samples={payload['samples']}") - print(f"roots_backend={payload['roots_backend']}") - component_counts = ",".join( - str(item) for item in payload["component_counts"] - ) - print(f"component_counts={component_counts}") - print( - "best_components_unique=" - + ",".join(str(item) for item in payload["best_components_unique"]) - ) - print( - "best_selected_model_orders_unique=" - + ",".join( - str(item) for item in payload["best_selected_model_orders_unique"] - ) - ) - for trace_payload in payload["traces"]: - print(f"trace_index={trace_payload['trace_index']}") - print(f" best_components={trace_payload['best_components']}") - print( - " best_selected_model_order=" - f"{trace_payload['best_selected_model_order']}" - ) - print( - " best_reconstruction_rms_error=" - f"{trace_payload['best_reconstruction_rms_error']:.6g}" - ) - - -def print_subspace_benchmark_batch(payload): - print(f"source={payload['source']['kind']}") - print(f"trace_count={payload['trace_count']}") - print(f"samples={payload['samples']}") - print(f"model_order={payload['model_order']}") - print(f"components={payload['components']}") - print( - "baseline_rms_residual_range=" - f"{payload['baseline_rms_residual_min']:.6g}:" - f"{payload['baseline_rms_residual_max']:.6g}" - ) - for summary in payload["method_summary"]: - print(f"method={summary['method']}") - print(f" svd_backend={summary['svd_backend']}") - print(f" trace_count={summary['trace_count']}") - print( - " rms_residual_range=" - f"{summary['rms_residual_min']:.6g}:" - f"{summary['rms_residual_max']:.6g}" - ) - print( - " max_abs_reconstruction_diff_max=" - f"{summary['max_abs_reconstruction_diff_max']:.6g}" - ) diff --git a/src/cuphoton/xray/commands.py b/src/cuphoton/xray/commands.py index c97852b0..e73fe23b 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -28,13 +28,6 @@ VariablePositionalInvariant, ) -from ._cli_output import ( - print_file_probe, - print_model_order_sweep_batch, - print_model_order_sweep_payload, - print_subspace_benchmark_batch, - print_validation_sweep, -) from ._types import FIT_DIAGNOSTICS_LEVELS from .doctor import ( collect_doctor_report, @@ -3772,10 +3765,9 @@ def _data_probe(args): result = probe_hdf5_pair(h5dir=args.h5dir, fon=args.fon, foff=args.foff) payload = result.to_dict() if args.json: - print(json.dumps(payload, indent=2, sort_keys=True)) + args._out(json.dumps(payload, indent=2, sort_keys=True)) else: - print_file_probe("on", result.on) - print_file_probe("off", result.off) + args._out(result) return 0 @@ -4410,11 +4402,7 @@ def _linear_prediction_smoke(args): def _linear_prediction_validate(args): - from .synthetic_validation import ( - build_summary, - validation_sweep, - write_validation_run, - ) + from .synthetic_validation import validate_linear_prediction snr = tuple(float(x) for x in str(args.snr_db).split(",") if x.strip()) distortion = None @@ -4429,30 +4417,22 @@ def _linear_prediction_validate(args): "--distortion must be kind:amount with a numeric amount " "(for example chirp:0.05)" ) from exc - sweeps = [ - validation_sweep( - samples=args.samples, - snr_db=snr, - trials=args.trials, - n_components=args.components, - seed=args.seed, - distortion=distortion, + report = validate_linear_prediction( + samples=args.samples, + snr_db=snr, + trials=args.trials, + n_components=args.components, + seed=args.seed, + distortion=distortion, + output_dir=args.output_dir, + ) + args._out( + json.dumps( + report.to_dict(), indent=2, sort_keys=True, allow_nan=False ) - ] - if args.output_dir is not None: - summary = write_validation_run(args.output_dir, sweeps) - else: - summary = build_summary(sweeps) - if args.json: - print(json.dumps(summary, indent=2, sort_keys=True, allow_nan=False)) - return 0 - print(f"samples={sweeps[0].samples}") - print(f"signal_rms={sweeps[0].signal_rms:.6g}") - for sweep in sweeps: - print_validation_sweep(sweep) - if args.output_dir is not None: - for key, name in summary["artifacts"].items(): - print(f"{key}={name if name else 'not written'}") + if args.json + else report + ) return 0 @@ -5439,132 +5419,54 @@ def _prediction_roots_benchmark(args): def _model_order_sweep(args): - from .linear_prediction import model_order_sweep + from .linear_prediction import model_order_sweep, model_order_sweep_batch components = _component_sweep_values(args) trace_paths = _trace_npz_batch_paths_from_args(args, required=False) + source = None if trace_paths: _reject_trace_npz_hdf5_options(args) - if len(trace_paths) == 1: - time, trace, source = _load_trace_npz(trace_paths[0]) - payload = _single_model_order_sweep_payload( - model_order_sweep( - time, - trace, - components, - roots_backend=args.roots_backend, - relative_tolerance=args.relative_tolerance, - ), - source=source, - ) - else: - time, trace_rows, sources = _load_trace_npz_batch(trace_paths) - payload = _model_order_sweep_batch_payload( - time, - trace_rows, - sources, - components=components, - roots_backend=args.roots_backend, - relative_tolerance=args.relative_tolerance, - ) - else: - time, trace, source = _trace_for_workbench(args) - payload = _single_model_order_sweep_payload( - model_order_sweep( - time, - trace, - components, - roots_backend=args.roots_backend, - relative_tolerance=args.relative_tolerance, - ), - source=source, + if len(trace_paths) > 1: + time, trace_rows, sources = _load_trace_npz_batch(trace_paths) + result = model_order_sweep_batch( + time, + trace_rows, + sources, + components=components, + roots_backend=args.roots_backend, + relative_tolerance=args.relative_tolerance, ) - - if args.json: - print(json.dumps(payload, indent=2, sort_keys=True)) - elif payload["source"]["kind"] == "trace-npz-batch": - print_model_order_sweep_batch(payload) else: - print_model_order_sweep_payload(payload) - return 0 - - -def _single_model_order_sweep_payload(result, *, source): - payload = dataclass_asdict(result) - payload["source"] = source - return payload - - -def _model_order_sweep_batch_payload( - time, - trace_rows, - sources, - *, - components, - roots_backend: str, - relative_tolerance: float, -): - from .linear_prediction import model_order_sweep - - trace_payloads = [] - best_components = [] - best_selected_model_orders = [] - best_rms_residuals = [] - best_reconstruction_errors = [] - for trace_index, (trace, source) in enumerate( - zip(trace_rows, sources, strict=True) - ): + time, trace, source = ( + _load_trace_npz(trace_paths[0]) + if trace_paths + else _trace_for_workbench(args) + ) result = model_order_sweep( time, trace, components, - roots_backend=roots_backend, - relative_tolerance=relative_tolerance, - ) - trace_payload = _single_model_order_sweep_payload( - result, - source=source, - ) - trace_payload["trace_index"] = trace_index - trace_payloads.append(trace_payload) - best_components.append(int(result.best_components)) - best_selected_model_orders.append( - int(result.best_selected_model_order) - ) - best_rms_residuals.append(float(result.best_rms_residual)) - best_reconstruction_errors.append( - float(result.best_reconstruction_rms_error) + roots_backend=args.roots_backend, + relative_tolerance=args.relative_tolerance, ) - return { - "source": { - "kind": "trace-npz-batch", - "trace_count": len(sources), - "traces": sources, - }, - "trace_count": len(sources), - "samples": int(time.shape[0]), - "roots_backend": roots_backend, - "relative_tolerance": float(relative_tolerance), - "component_counts": tuple(int(item) for item in components), - "best_components_by_trace": tuple(best_components), - "best_components_unique": tuple(sorted(set(best_components))), - "best_selected_model_orders_by_trace": tuple( - best_selected_model_orders - ), - "best_selected_model_orders_unique": tuple( - sorted(set(best_selected_model_orders)) - ), - "best_rms_residual_min": min(best_rms_residuals), - "best_rms_residual_max": max(best_rms_residuals), - "best_reconstruction_rms_error_min": min(best_reconstruction_errors), - "best_reconstruction_rms_error_max": max(best_reconstruction_errors), - "traces": tuple(trace_payloads), - } + if args.json: + payload = result.to_dict() + if source is not None: + payload["source"] = source + args._out(json.dumps(payload, indent=2, sort_keys=True)) + else: + if source is not None: + args._out(f"source={source['kind']}") + args._out(result) + return 0 def _subspace_benchmark(args): - from .subspace import compare_subspace_methods + from .subspace import ( + compare_subspace_methods, + compare_subspace_methods_batch, + ) trace_paths = _trace_npz_batch_paths_from_args(args, required=False) svd_backends = _subspace_svd_backends(args) @@ -5573,191 +5475,54 @@ def _subspace_benchmark(args): if args.method is None else tuple(args.method) ) + source = None if trace_paths: _reject_trace_npz_hdf5_options(args) - if len(trace_paths) == 1: - time, trace, source = _load_trace_npz(trace_paths[0]) - payload = _single_subspace_benchmark_payload( - compare_subspace_methods( - time, - trace, - model_order=args.model_order, - components=args.components, - methods=methods, - svd_backends=svd_backends, - pencil_rows=args.pencil_rows, - randomized_oversamples=args.randomized_oversamples, - power_iterations=args.power_iterations, - random_seed=args.random_seed, - ), - source=source, - svd_backends=svd_backends, - ) - else: - time, trace_rows, sources = _load_trace_npz_batch(trace_paths) - payload = _subspace_benchmark_batch_payload( - time, - trace_rows, - sources, - model_order=args.model_order, - components=args.components, - methods=methods, - svd_backends=svd_backends, - pencil_rows=args.pencil_rows, - randomized_oversamples=args.randomized_oversamples, - power_iterations=args.power_iterations, - random_seed=args.random_seed, - ) - else: - time, trace, source = _trace_for_workbench(args) - payload = _single_subspace_benchmark_payload( - compare_subspace_methods( - time, - trace, - model_order=args.model_order, - components=args.components, - methods=methods, - svd_backends=svd_backends, - pencil_rows=args.pencil_rows, - randomized_oversamples=args.randomized_oversamples, - power_iterations=args.power_iterations, - random_seed=args.random_seed, - ), - source=source, + if len(trace_paths) > 1: + time, trace_rows, sources = _load_trace_npz_batch(trace_paths) + result = compare_subspace_methods_batch( + time, + trace_rows, + sources, + model_order=args.model_order, + components=args.components, + methods=methods, svd_backends=svd_backends, + pencil_rows=args.pencil_rows, + randomized_oversamples=args.randomized_oversamples, + power_iterations=args.power_iterations, + random_seed=args.random_seed, ) - - if args.json: - print(json.dumps(payload, indent=2, sort_keys=True)) - elif payload["source"]["kind"] == "trace-npz-batch": - print_subspace_benchmark_batch(payload) else: - print(f"source={payload['source']['kind']}") - print(f"samples={payload['samples']}") - print(f"model_order={payload['model_order']}") - print(f"components={payload['components']}") - print(f"baseline_chi2={payload['baseline_chi2']:.6g}") - print(f"baseline_rms_residual={payload['baseline_rms_residual']:.6g}") - for method in payload["methods"]: - print(f"method={method['method']}") - print(f" svd_backend={method['svd_backend']}") - print(f" svd_rank={method['svd_rank']}") - print(f" rms_residual={method['rms_residual']:.6g}") - print( - " max_abs_reconstruction_diff=" - f"{method['max_abs_reconstruction_diff']:.6g}" - ) - print(f" elapsed_s={method['elapsed_s']:.6g}") - return 0 - - -def _single_subspace_benchmark_payload(result, *, source, svd_backends): - payload = dataclass_asdict(result) - payload["source"] = source - payload["svd_backends"] = svd_backends - return payload - - -def _subspace_benchmark_batch_payload( - time, - trace_rows, - sources, - *, - model_order: int, - components: int, - methods, - svd_backends, - pencil_rows: int | None, - randomized_oversamples: int, - power_iterations: int, - random_seed: int, -): - from .subspace import compare_subspace_methods - - method_summary = { - (method, svd_backend): { - "method": method, - "svd_backend": svd_backend, - "trace_count": 0, - "rms_residual_min": float("inf"), - "rms_residual_max": 0.0, - "max_abs_reconstruction_diff_max": 0.0, - "elapsed_s_total": 0.0, - } - for method in methods - for svd_backend in svd_backends - } - trace_payloads = [] - baseline_chi2_values = [] - baseline_rms_values = [] - - for trace_index, (trace, source) in enumerate( - zip(trace_rows, sources, strict=True) - ): + time, trace, source = ( + _load_trace_npz(trace_paths[0]) + if trace_paths + else _trace_for_workbench(args) + ) result = compare_subspace_methods( time, trace, - model_order=model_order, - components=components, + model_order=args.model_order, + components=args.components, methods=methods, svd_backends=svd_backends, - pencil_rows=pencil_rows, - randomized_oversamples=randomized_oversamples, - power_iterations=power_iterations, - random_seed=random_seed, - ) - trace_payload = _single_subspace_benchmark_payload( - result, - source=source, - svd_backends=svd_backends, + pencil_rows=args.pencil_rows, + randomized_oversamples=args.randomized_oversamples, + power_iterations=args.power_iterations, + random_seed=args.random_seed, ) - trace_payload["trace_index"] = trace_index - trace_payloads.append(trace_payload) - baseline_chi2_values.append(float(result.baseline_chi2)) - baseline_rms_values.append(float(result.baseline_rms_residual)) - for method_result in result.methods: - summary = method_summary[ - (method_result.method, method_result.svd_backend) - ] - summary["trace_count"] += 1 - summary["rms_residual_min"] = min( - float(summary["rms_residual_min"]), - float(method_result.rms_residual), - ) - summary["rms_residual_max"] = max( - float(summary["rms_residual_max"]), - float(method_result.rms_residual), - ) - summary["max_abs_reconstruction_diff_max"] = max( - float(summary["max_abs_reconstruction_diff_max"]), - float(method_result.max_abs_reconstruction_diff), - ) - summary["elapsed_s_total"] += float(method_result.elapsed_s) - - summaries = tuple( - method_summary[(method, svd_backend)] - for method in methods - for svd_backend in svd_backends - ) - return { - "source": { - "kind": "trace-npz-batch", - "trace_count": len(sources), - "traces": sources, - }, - "trace_count": len(sources), - "samples": int(time.shape[0]), - "model_order": int(model_order), - "components": int(components), - "svd_backends": svd_backends, - "baseline_chi2_min": min(baseline_chi2_values), - "baseline_chi2_max": max(baseline_chi2_values), - "baseline_rms_residual_min": min(baseline_rms_values), - "baseline_rms_residual_max": max(baseline_rms_values), - "method_summary": summaries, - "traces": tuple(trace_payloads), - } + if args.json: + payload = result.to_dict() + if source is not None: + payload["source"] = source + payload["svd_backends"] = svd_backends + args._out(json.dumps(payload, indent=2, sort_keys=True)) + else: + if source is not None: + args._out(f"source={source['kind']}") + args._out(result) + return 0 def _subspace_acceptance(args): diff --git a/src/cuphoton/xray/hdf5.py b/src/cuphoton/xray/hdf5.py index bbae346f..3e0df47e 100644 --- a/src/cuphoton/xray/hdf5.py +++ b/src/cuphoton/xray/hdf5.py @@ -39,6 +39,15 @@ class DatasetProbe: dtype: str chunks: tuple[int, ...] | None + def __str__(self) -> str: + shape = "x".join(str(dim) for dim in self.shape) + chunks = ( + "-" + if self.chunks is None + else "x".join(str(dim) for dim in self.chunks) + ) + return f"{self.name} shape={shape} dtype={self.dtype} chunks={chunks}" + @dataclass(frozen=True) class FileProbe: @@ -50,12 +59,38 @@ class FileProbe: ipm_pairs: tuple[str, ...] load_plan: Hdf5LoadPlan | None + def __str__(self) -> str: + return self.format_text() + + def format_text(self, *, label: str | None = None) -> str: + """Summarize the file, optionally prefixing fields with a label.""" + prefix = "" if label is None else f"{label}_" + lines = [ + f"{prefix}path={self.path}", + f"{prefix}size_bytes={self.size_bytes}", + f"{prefix}schema={self.schema}", + f"{prefix}ipm_pairs={','.join(self.ipm_pairs) or '-'}", + f"{prefix}keys={','.join(self.keys)}", + ] + lines.extend( + f"{prefix}dataset={dataset}" for dataset in self.datasets + ) + return "\n".join(lines) + @dataclass(frozen=True) class Hdf5PairProbe: on: FileProbe off: FileProbe + def __str__(self) -> str: + return "\n".join( + ( + self.on.format_text(label="on"), + self.off.format_text(label="off"), + ) + ) + def to_dict(self) -> dict[str, Any]: return { "on": _file_probe_to_dict(self.on), diff --git a/src/cuphoton/xray/linear_prediction.py b/src/cuphoton/xray/linear_prediction.py index 8d522a3b..6c0f56b4 100644 --- a/src/cuphoton/xray/linear_prediction.py +++ b/src/cuphoton/xray/linear_prediction.py @@ -4,7 +4,8 @@ from __future__ import annotations -from dataclasses import dataclass +from collections.abc import Sequence +from dataclasses import asdict, dataclass from math import pi from time import perf_counter from typing import Any @@ -431,6 +432,17 @@ class ModelOrderSweepEntry: singular_value_ratio: float | None elapsed_s: float + def __str__(self) -> str: + return ( + f"components={self.components} " + f"selected_model_order={self.selected_model_order} " + f"rms_residual={self.rms_residual:.6g} " + f"reconstruction_rms_error={self.reconstruction_rms_error:.6g} " + f"chi2={self.chi2:.6g} " + f"selected_roots={self.selected_root_count} " + f"decaying_roots={self.decaying_root_count}" + ) + @dataclass(frozen=True) class ModelOrderSweep: @@ -443,6 +455,115 @@ class ModelOrderSweep: best_reconstruction_rms_error: float entries: tuple[ModelOrderSweepEntry, ...] + def to_dict(self) -> dict[str, Any]: + return asdict(self) + + def __str__(self) -> str: + return "\n".join( + ( + f"samples={self.samples}", + f"roots_backend={self.roots_backend}", + f"best_components={self.best_components}", + f"best_selected_model_order={self.best_selected_model_order}", + f"best_rms_residual={self.best_rms_residual:.6g}", + "best_reconstruction_rms_error=" + f"{self.best_reconstruction_rms_error:.6g}", + *(str(entry) for entry in self.entries), + ) + ) + + +@dataclass(frozen=True) +class ModelOrderSweepBatch: + """Per-trace model-order sweeps and their input provenance.""" + + samples: int + roots_backend: str + relative_tolerance: float + component_counts: tuple[int, ...] + sources: tuple[dict[str, Any], ...] + traces: tuple[ModelOrderSweep, ...] + + def to_dict(self) -> dict[str, Any]: + best_components = tuple( + trace.best_components for trace in self.traces + ) + best_orders = tuple( + trace.best_selected_model_order for trace in self.traces + ) + best_residuals = tuple( + trace.best_rms_residual for trace in self.traces + ) + best_errors = tuple( + trace.best_reconstruction_rms_error for trace in self.traces + ) + trace_payloads = [] + for index, (trace, source) in enumerate( + zip(self.traces, self.sources, strict=True) + ): + payload = trace.to_dict() + payload["source"] = source + payload["trace_index"] = index + trace_payloads.append(payload) + return { + "source": { + "kind": "trace-npz-batch", + "trace_count": len(self.traces), + "traces": self.sources, + }, + "trace_count": len(self.traces), + "samples": self.samples, + "roots_backend": self.roots_backend, + "relative_tolerance": self.relative_tolerance, + "component_counts": self.component_counts, + "best_components_by_trace": best_components, + "best_components_unique": tuple(sorted(set(best_components))), + "best_selected_model_orders_by_trace": best_orders, + "best_selected_model_orders_unique": tuple( + sorted(set(best_orders)) + ), + "best_rms_residual_min": min(best_residuals), + "best_rms_residual_max": max(best_residuals), + "best_reconstruction_rms_error_min": min(best_errors), + "best_reconstruction_rms_error_max": max(best_errors), + "traces": tuple(trace_payloads), + } + + def __str__(self) -> str: + lines = [ + "source=trace-npz-batch", + f"trace_count={len(self.traces)}", + f"samples={self.samples}", + f"roots_backend={self.roots_backend}", + "component_counts=" + ",".join(map(str, self.component_counts)), + "best_components_unique=" + + ",".join( + str(item) + for item in sorted( + {trace.best_components for trace in self.traces} + ) + ), + "best_selected_model_orders_unique=" + + ",".join( + str(item) + for item in sorted( + {trace.best_selected_model_order for trace in self.traces} + ) + ), + ] + for index, trace in enumerate(self.traces): + lines.extend( + ( + f"trace_index={index}", + f" best_components={trace.best_components}", + " best_selected_model_order=" + f"{trace.best_selected_model_order}", + " best_reconstruction_rms_error=" + f"{trace.best_reconstruction_rms_error:.6g}", + ) + ) + return "\n".join(lines) + @dataclass(frozen=True) class _PredictionRootsResult: @@ -1268,6 +1389,39 @@ def model_order_sweep( ) +def model_order_sweep_batch( + time, + trace_rows, + sources: Sequence[dict[str, Any]], + *, + components: tuple[int, ...], + roots_backend: str = "eigvals", + relative_tolerance: float = 0.01, +) -> ModelOrderSweepBatch: + """Sweep component counts for rows sharing one time axis.""" + + results = tuple( + model_order_sweep( + time, + trace, + components, + roots_backend=roots_backend, + relative_tolerance=relative_tolerance, + ) + for trace, _source in zip(trace_rows, sources, strict=True) + ) + if not results: + raise ValueError("trace_rows must contain at least one trace") + return ModelOrderSweepBatch( + samples=results[0].samples, + roots_backend=roots_backend, + relative_tolerance=float(relative_tolerance), + component_counts=tuple(int(item) for item in components), + sources=tuple(sources), + traces=results, + ) + + def benchmark_linear_prediction_p1_batch( *, samples: int = 96, diff --git a/src/cuphoton/xray/subspace.py b/src/cuphoton/xray/subspace.py index 3088aa1d..98afb635 100644 --- a/src/cuphoton/xray/subspace.py +++ b/src/cuphoton/xray/subspace.py @@ -4,7 +4,7 @@ from __future__ import annotations -from dataclasses import dataclass +from dataclasses import asdict, dataclass from time import perf_counter import numpy as np @@ -27,6 +27,22 @@ class SubspaceMethodBenchmark: roots_real: tuple[float, ...] roots_imag: tuple[float, ...] + def __str__(self) -> str: + return "\n".join( + ( + f"method={self.method}", + f" svd_backend={self.svd_backend}", + f" svd_rank={self.svd_rank}", + f" rms_residual={self.rms_residual:.6g}", + " max_abs_reconstruction_diff=" + f"{self.max_abs_reconstruction_diff:.6g}", + f" elapsed_s={self.elapsed_s:.6g}", + ) + ) + + def to_dict(self) -> dict[str, object]: + return asdict(self) + @dataclass(frozen=True) class SubspaceBenchmark: @@ -37,6 +53,159 @@ class SubspaceBenchmark: baseline_rms_residual: float methods: tuple[SubspaceMethodBenchmark, ...] + def __str__(self) -> str: + return "\n".join( + ( + f"samples={self.samples}", + f"model_order={self.model_order}", + f"components={self.components}", + f"baseline_chi2={self.baseline_chi2:.6g}", + f"baseline_rms_residual={self.baseline_rms_residual:.6g}", + *(str(method) for method in self.methods), + ) + ) + + def to_dict(self) -> dict[str, object]: + return asdict(self) + + +@dataclass(frozen=True) +class _SubspaceMethodSummary: + method: str + svd_backend: str + trace_count: int + rms_residual_min: float + rms_residual_max: float + max_abs_reconstruction_diff_max: float + elapsed_s_total: float + + def __str__(self) -> str: + return "\n".join( + ( + f"method={self.method}", + f" svd_backend={self.svd_backend}", + f" trace_count={self.trace_count}", + " rms_residual_range=" + f"{self.rms_residual_min:.6g}:{self.rms_residual_max:.6g}", + " max_abs_reconstruction_diff_max=" + f"{self.max_abs_reconstruction_diff_max:.6g}", + ) + ) + + +@dataclass(frozen=True) +class SubspaceBenchmarkBatch: + samples: int + model_order: int + components: int + methods: tuple[str, ...] + svd_backends: tuple[str, ...] + sources: tuple[dict[str, object], ...] + traces: tuple[SubspaceBenchmark, ...] + + @property + def method_summary(self) -> tuple[_SubspaceMethodSummary, ...]: + summaries = [] + for method in self.methods: + for svd_backend in self.svd_backends: + results = tuple( + result + for trace in self.traces + for result in trace.methods + if (result.method, result.svd_backend) + == (method, svd_backend) + ) + summaries.append( + _SubspaceMethodSummary( + method=method, + svd_backend=svd_backend, + trace_count=len(results), + rms_residual_min=min( + [ + float("inf"), + *(result.rms_residual for result in results), + ] + ), + rms_residual_max=max( + [ + 0.0, + *(result.rms_residual for result in results), + ] + ), + max_abs_reconstruction_diff_max=max( + [ + 0.0, + *( + result.max_abs_reconstruction_diff + for result in results + ), + ] + ), + elapsed_s_total=sum( + (result.elapsed_s for result in results), 0.0 + ), + ) + ) + return tuple(summaries) + + def __str__(self) -> str: + baseline_rms = tuple( + trace.baseline_rms_residual for trace in self.traces + ) + return "\n".join( + ( + "source=trace-npz-batch", + f"trace_count={len(self.traces)}", + f"samples={self.samples}", + f"model_order={self.model_order}", + f"components={self.components}", + "baseline_rms_residual_range=" + f"{min(baseline_rms):.6g}:{max(baseline_rms):.6g}", + *(str(summary) for summary in self.method_summary), + ) + ) + + def to_dict(self) -> dict[str, object]: + traces = tuple( + { + **trace.to_dict(), + "source": source, + "svd_backends": self.svd_backends, + "trace_index": index, + } + for index, (trace, source) in enumerate( + zip(self.traces, self.sources, strict=True) + ) + ) + return { + "source": { + "kind": "trace-npz-batch", + "trace_count": len(self.sources), + "traces": self.sources, + }, + "trace_count": len(self.sources), + "samples": self.samples, + "model_order": self.model_order, + "components": self.components, + "svd_backends": self.svd_backends, + "baseline_chi2_min": min( + trace.baseline_chi2 for trace in self.traces + ), + "baseline_chi2_max": max( + trace.baseline_chi2 for trace in self.traces + ), + "baseline_rms_residual_min": min( + trace.baseline_rms_residual for trace in self.traces + ), + "baseline_rms_residual_max": max( + trace.baseline_rms_residual for trace in self.traces + ), + "method_summary": tuple( + asdict(summary) for summary in self.method_summary + ), + "traces": traces, + } + @dataclass(frozen=True) class _SubspaceSvd: @@ -137,6 +306,53 @@ def compare_subspace_methods( ) +def compare_subspace_methods_batch( + time, + trace_rows, + sources, + *, + model_order: int, + components: int, + methods: tuple[str, ...] = ("matrix-pencil", "esprit"), + svd_backends: tuple[str, ...] = ("full",), + pencil_rows: int | None = None, + randomized_oversamples: int = 8, + power_iterations: int = 1, + random_seed: int = 0, +) -> SubspaceBenchmarkBatch: + """Compare trace rows while retaining their individual source metadata.""" + + results = [] + trace_sources = [] + for trace, source in zip(trace_rows, sources, strict=True): + results.append( + compare_subspace_methods( + time, + trace, + model_order=model_order, + components=components, + methods=methods, + svd_backends=svd_backends, + pencil_rows=pencil_rows, + randomized_oversamples=randomized_oversamples, + power_iterations=power_iterations, + random_seed=random_seed, + ) + ) + trace_sources.append(source) + if not results: + raise ValueError("at least one trace is required") + return SubspaceBenchmarkBatch( + samples=int(np.asarray(time).shape[0]), + model_order=int(model_order), + components=int(components), + methods=tuple(methods), + svd_backends=tuple(svd_backends), + sources=tuple(trace_sources), + traces=tuple(results), + ) + + def synthetic_subspace_benchmark( *, samples: int = 96, @@ -458,8 +674,10 @@ def _sort_roots(roots: np.ndarray): __all__ = [ "SubspaceBenchmark", + "SubspaceBenchmarkBatch", "SubspaceMethodBenchmark", "compare_subspace_methods", + "compare_subspace_methods_batch", "esprit_roots", "matrix_pencil_roots", "synthetic_subspace_benchmark", diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index 471eaf4e..8d43e1e4 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -313,6 +313,19 @@ class ModeValidation: recovered_trials: int loss_rate: float + def __str__(self) -> str: + w = self.angular_frequency + d = self.decay + return ( + f"mode w={w.truth:g} decay={d.truth:g}: " + f"freq std/crlb={w.std:.3g}/{w.crlb_std:.3g} " + f"({w.std_over_crlb_std:.2f}x) bias={w.bias:+.3g}; " + f"decay std/crlb={d.std:.3g}/{d.crlb_std:.3g} " + f"({d.std_over_crlb_std:.2f}x) bias={d.bias:+.3g}; " + f"recovered={self.recovered_trials} " + f"loss_rate={self.loss_rate:.3f}" + ) + @property def frequency_ratio(self) -> float: return self.angular_frequency.std_over_crlb_std @@ -339,6 +352,17 @@ class ValidationLevel: any_mode_lost_rate: float residual_ratio: float + def __str__(self) -> str: + lines = [ + f"snr_db={self.snr_db:g} sigma={self.noise_sigma:.4g} " + f"trials={self.trials_successful}/{self.trials_attempted} " + f"estimator_errors={self.estimator_errors} " + f"any_mode_lost_rate={self.any_mode_lost_rate:.3f} " + f"residual_ratio={self.residual_ratio:.2f}" + ] + lines.extend(f" {mode}" for mode in self.modes) + return "\n".join(lines) + @dataclass(frozen=True) class ValidationSweep: @@ -354,10 +378,43 @@ class ValidationSweep: distortion: tuple[str, float] | None = None extra: dict = field(default_factory=dict) + def __str__(self) -> str: + lines = [f"estimator={self.backend}"] + if self.distortion is not None: + lines.append( + f"distortion={self.distortion[0]}:{self.distortion[1]:g}" + ) + lines.extend(str(level) for level in self.levels) + return "\n".join(lines) + def to_dict(self) -> dict: return asdict(self) +@dataclass(frozen=True) +class ValidationReport: + """Validation sweeps with their serializable summary and artifacts.""" + + sweeps: tuple[ValidationSweep, ...] + summary: dict[str, Any] + + def __str__(self) -> str: + first = self.sweeps[0] + lines = [ + f"samples={first.samples}", + f"signal_rms={first.signal_rms:.6g}", + ] + lines.extend(str(sweep) for sweep in self.sweeps) + lines.extend( + f"{key}={name if name else 'not written'}" + for key, name in self.summary["artifacts"].items() + ) + return "\n".join(lines) + + def to_dict(self) -> dict[str, Any]: + return self.summary + + def _wrap_phase(values: np.ndarray) -> np.ndarray: return (values + np.pi) % (2.0 * np.pi) - np.pi @@ -563,6 +620,40 @@ def validation_sweep( ) +def validate_linear_prediction( + *, + samples: int = 96, + snr_db: Sequence[float] = (40.0, 30.0, 20.0, 10.0), + trials: int = 200, + n_components: int = 6, + seed: int = 20260914, + distortion: tuple[str, float] | None = None, + output_dir: Path | str | None = None, +) -> ValidationReport: + """Validate linear prediction and optionally write its run artifacts. + + The returned report provides text and strict-JSON-compatible summaries + of the same sweep. Set ``output_dir`` to write ``summary.json`` and, + when Bokeh is available, a validation figure. + """ + sweeps = ( + validation_sweep( + samples=samples, + snr_db=snr_db, + trials=trials, + n_components=n_components, + seed=seed, + distortion=distortion, + ), + ) + summary = ( + build_summary(sweeps) + if output_dir is None + else write_validation_run(output_dir, sweeps) + ) + return ValidationReport(sweeps=sweeps, summary=summary) + + def _source_provenance() -> tuple[str | None, bool | None]: """Revision and tracked-file dirtiness for this module's checkout.""" try: diff --git a/tests/xray/test_cli.py b/tests/xray/test_cli.py index eec8d033..eed6adbe 100644 --- a/tests/xray/test_cli.py +++ b/tests/xray/test_cli.py @@ -4,6 +4,7 @@ from __future__ import annotations +import io import json import numpy as np @@ -37,6 +38,52 @@ def test_cli_help(capsys): assert "GPU-accelerated X-ray detector analysis tools" in captured.out +@pytest.mark.parametrize("as_json", [False, True]) +@pytest.mark.parametrize( + ("command_name", "marker"), + [ + ("DataProbeCommand", "on"), + ("LinearPredictionValidateCommand", "samples"), + ("ModelOrderSweepCommand", "trace_count"), + ("SubspaceBenchmarkCommand", "trace_count"), + ], +) +def test_report_commands_use_injected_output( + command_name, marker, as_json, tmp_path, capsys +): + from cuphoton.xray import commands + + options = {"json": as_json} + if command_name == "DataProbeCommand": + h5py = pytest.importorskip("h5py") + with h5py.File(tmp_path / "cube.h5", "w") as handle: + handle.create_dataset("imgs", data=np.zeros((2, 3, 4))) + options.update(h5dir=tmp_path, fon="cube.h5", foff="cube.h5") + elif command_name == "LinearPredictionValidateCommand": + options.update(samples=48, snr_db="30", trials=2) + else: + time, rows = synthetic_trace_batch(samples=48, traces=2) + paths = [tmp_path / f"trace-{index}.npz" for index in range(2)] + for path, trace in zip(paths, rows, strict=True): + np.savez(path, time=time, trace=trace) + options["trace_npz"] = paths + if command_name == "ModelOrderSweepCommand": + options["component_values"] = [2, 4] + + output, error = io.StringIO(), io.StringIO() + command = getattr(commands, command_name)(output=output, error=error) + command.start(options) + + captured = capsys.readouterr() + assert captured.out == captured.err == error.getvalue() == "" + text = output.getvalue() + assert text.endswith("\n") + if as_json: + assert isinstance(json.loads(text), dict) + else: + assert marker in text + + def test_app_spec_uses_canonical_namespace(): spec = get_component("xray") diff --git a/tests/xray/test_linear_prediction.py b/tests/xray/test_linear_prediction.py index 0a9d9072..d78ab3c2 100644 --- a/tests/xray/test_linear_prediction.py +++ b/tests/xray/test_linear_prediction.py @@ -35,6 +35,7 @@ linear_prediction_p2_artifacts_numpy, linear_prediction_variable_artifacts_cupy_batched, model_order_sweep, + model_order_sweep_batch, synthetic_prediction_coefficients, synthetic_trace, synthetic_trace_batch, @@ -1910,6 +1911,55 @@ def test_model_order_sweep_chooses_smallest_qualified_order(): assert result.best_selected_model_order == 2 +def test_model_order_sweep_batch_retains_results_and_provenance(): + time, rows = synthetic_trace_batch(samples=48, traces=2) + sources = [ + {"kind": "trace-npz", "path": "first.npz"}, + {"kind": "trace-npz", "path": "second.npz"}, + ] + batch = model_order_sweep_batch( + time.tolist(), + rows, + sources, + components=(6, 2, 4), + relative_tolerance=0.05, + ) + expected = [ + model_order_sweep(time, row, (6, 2, 4), relative_tolerance=0.05) + for row in rows + ] + + payload = batch.to_dict() + assert batch.sources == tuple(sources) + assert payload["source"]["traces"] == tuple(sources) + assert payload["trace_count"] == payload["source"]["trace_count"] == 2 + assert payload["component_counts"] == (6, 2, 4) + assert payload["best_components_by_trace"] == tuple( + result.best_components for result in expected + ) + assert payload["best_selected_model_orders_unique"] == tuple( + sorted({result.best_selected_model_order for result in expected}) + ) + residuals = [result.best_rms_residual for result in expected] + assert payload["best_rms_residual_min"] == pytest.approx(min(residuals)) + assert payload["best_rms_residual_max"] == pytest.approx(max(residuals)) + for index, trace in enumerate(payload["traces"]): + assert trace["trace_index"] == index + assert trace["source"] == sources[index] + assert trace["best_components"] == expected[index].best_components + assert trace["best_rms_residual"] == pytest.approx(residuals[index]) + assert trace["entries"] == batch.traces[index].to_dict()["entries"] + + +@pytest.mark.parametrize("source_count", [1, 3]) +def test_model_order_sweep_batch_rejects_unmatched_provenance(source_count): + time, rows = synthetic_trace_batch(samples=32, traces=2) + sources = [{"kind": "trace-npz"}] * source_count + + with pytest.raises(ValueError, match="zip"): + model_order_sweep_batch(time, rows, sources, components=(2,)) + + def test_compare_cpu_gpu_can_skip_gpu(): comparison = compare_cpu_gpu( samples=64, diff --git a/tests/xray/test_subspace.py b/tests/xray/test_subspace.py index 4445a0a3..56ce18ee 100644 --- a/tests/xray/test_subspace.py +++ b/tests/xray/test_subspace.py @@ -12,6 +12,7 @@ from cuphoton.xray.linear_prediction import synthetic_trace from cuphoton.xray.subspace import ( compare_subspace_methods, + compare_subspace_methods_batch, esprit_roots, matrix_pencil_roots, ) @@ -120,6 +121,104 @@ def test_compare_subspace_methods_can_use_partial_svd(): assert len(method.roots_real) == 5 +def test_subspace_batch_preserves_trace_results_and_source_order(): + time, trace = synthetic_trace(samples=64) + trace_rows = np.stack((trace, trace + 0.1 * np.cos(3.0 * time))) + sources = ( + {"kind": "trace-npz", "path": "trace-y2.npz"}, + { + "kind": "trace-npz", + "path": "trace-y10.npz", + "summary": {"detector_row": 10}, + "time_matches_first": True, + }, + ) + options = { + "model_order": 4, + "components": 6, + "methods": ("esprit", "matrix-pencil"), + "svd_backends": ("randomized", "full"), + "pencil_rows": 24, + "randomized_oversamples": 3, + "power_iterations": 2, + "random_seed": 9, + } + batch = compare_subspace_methods_batch( + time, trace_rows, sources, **options + ) + payload = batch.to_dict() + + assert payload["source"]["traces"] == sources + assert payload["trace_count"] == 2 + for index, (trace_result, trace_payload) in enumerate( + zip(batch.traces, payload["traces"], strict=True) + ): + expected = compare_subspace_methods( + time, trace_rows[index], **options + ) + assert trace_payload["trace_index"] == index + assert trace_payload["source"] == sources[index] + assert trace_payload["svd_backends"] == options["svd_backends"] + assert trace_result.baseline_chi2 == expected.baseline_chi2 + for result, expected_method in zip( + trace_result.methods, expected.methods, strict=True + ): + assert result.pencil_rows == 24 + assert result.rms_residual == pytest.approx( + expected_method.rms_residual + ) + np.testing.assert_allclose( + result.roots_real, expected_method.roots_real + ) + np.testing.assert_allclose( + result.roots_imag, expected_method.roots_imag + ) + + assert [ + (summary.method, summary.svd_backend) + for summary in batch.method_summary + ] == [ + ("esprit", "randomized"), + ("esprit", "full"), + ("matrix-pencil", "randomized"), + ("matrix-pencil", "full"), + ] + for summary in batch.method_summary: + methods = [ + method + for result in batch.traces + for method in result.methods + if (method.method, method.svd_backend) + == (summary.method, summary.svd_backend) + ] + assert summary.trace_count == 2 + assert summary.rms_residual_min == min( + method.rms_residual for method in methods + ) + assert summary.rms_residual_max == max( + method.rms_residual for method in methods + ) + assert summary.max_abs_reconstruction_diff_max == max( + method.max_abs_reconstruction_diff for method in methods + ) + assert summary.elapsed_s_total == pytest.approx( + sum(method.elapsed_s for method in methods) + ) + + assert payload["baseline_chi2_min"] == min( + result.baseline_chi2 for result in batch.traces + ) + assert payload["baseline_chi2_max"] == max( + result.baseline_chi2 for result in batch.traces + ) + assert payload["baseline_rms_residual_min"] == min( + result.baseline_rms_residual for result in batch.traces + ) + assert payload["baseline_rms_residual_max"] == max( + result.baseline_rms_residual for result in batch.traces + ) + + def test_matrix_pencil_rejects_rankless_trace(): with pytest.raises(ValueError, match="rank"): matrix_pencil_roots(np.zeros(32), model_order=5) diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py index ac647da2..6d5033c8 100644 --- a/tests/xray/test_synthetic_validation.py +++ b/tests/xray/test_synthetic_validation.py @@ -21,6 +21,7 @@ distort_trace, match_modes, synthetic_modes_trace, + validate_linear_prediction, validation_sweep, write_validation_run, ) @@ -52,7 +53,7 @@ def test_estimator_recovers_known_modes_exactly_at_zero_noise(): tolerance=0.05, ) assert all(h is not None for h in hit) - for m, h in zip(fx.modes, hit): + for m, h in zip(fx.modes, hit, strict=True): assert abs(float(r.angular_frequency[h]) - m.angular_frequency) < 1e-6 assert abs(float(r.decay[h]) - m.decay) < 1e-6 @@ -266,6 +267,39 @@ def test_summary_schema_and_run_artifacts(tmp_path): build_summary([]) +@pytest.mark.parametrize("write_artifacts", [False, True]) +def test_validation_report_summarizes_the_run(tmp_path, write_artifacts): + out = tmp_path / "run" if write_artifacts else None + report = validate_linear_prediction( + samples=64, + snr_db=(30.0,), + trials=1, + n_components=4, + seed=7, + distortion=("chirp", 0.05), + output_dir=out, + ) + summary = report.to_dict() + assert json.loads(json.dumps(summary, allow_nan=False)) == summary + config = summary["config"] + assert config["sampling"]["samples"] == 64 + assert config["seed"] == 7 and config["trials_per_level"] == 1 + assert config["n_components"] == 4 + assert config["distortions"] == ["chirp:0.05"] + assert report.sweeps[0].levels[0].snr_db == 30.0 + text = str(report) + assert text.startswith("samples=64\nsignal_rms=") + assert "estimator=cpu\ndistortion=chirp:0.05\nsnr_db=30 " in text + if out is not None: + assert json.loads((out / "summary.json").read_text()) == summary + assert "summary=summary.json" in text + for name in summary["artifacts"].values(): + assert name is None or (out / name).is_file() + else: + assert summary["artifacts"] == {} + assert not (tmp_path / "run").exists() + + @pytest.mark.parametrize("trials", [1, 2]) def test_small_sweeps_write_strict_json_and_two_trials_have_std( tmp_path, trials From 832652136c3d1057d057bee3a9bc98414e6d0498 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Thu, 8 Oct 2026 12:57:07 -0700 Subject: [PATCH 09/10] Correct XRay bounds and equivalent mode representations Use analytic derivatives and scaled singular values to reject unidentifiable model parameters before reporting Cramer-Rao bounds. Normalize frequency signs with phases so equivalent traces report the same recovery and parameter statistics. Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-VALIDATION.md | 21 ++-- src/cuphoton/xray/linear_prediction.py | 2 + src/cuphoton/xray/synthetic_validation.py | 109 ++++++++++++------ tests/xray/test_synthetic_validation.py | 128 +++++++++++++++++++--- 4 files changed, 205 insertions(+), 55 deletions(-) diff --git a/docs/xray/LINEAR-PREDICTION-VALIDATION.md b/docs/xray/LINEAR-PREDICTION-VALIDATION.md index 947a5aed..614c6fb8 100644 --- a/docs/xray/LINEAR-PREDICTION-VALIDATION.md +++ b/docs/xray/LINEAR-PREDICTION-VALIDATION.md @@ -43,15 +43,19 @@ SNR levels must be finite and produce finite, positive noise scales; ## Bounds -`cramer_rao_bounds` forms the Fisher information `J^T J / sigma^2` from the -numerical Jacobian of the model at the true parameters (all parameters free, -including the constant) and returns the square root of the diagonal of its -inverse per parameter. These are standard-deviation bounds; the summary +`cramer_rao_bounds` uses the analytic Jacobian of the model at the true +parameters (all parameters free, including the constant). A column-scaled +singular-value decomposition computes the bounds associated with the Fisher +information `J^T J / sigma^2` without forming or inverting that matrix. +Numerically rank-deficient models, including modes with identical frequency +and decay but different phases, raise `ValueError`: the individual parameters +cannot be identified. These are standard-deviation bounds; the summary also carries the squared value as `crlb_variance`. For the single undamped sinusoid test fixture, the angular-frequency bound matches the square root of the closed-form variance `24 sigma^2 / (A^2 dt^2 N (N^2 - 1))` -(Kay, Estimation Theory, 1993) to within 3 percent. The numerical calculation -also frees the decay and constant. +(Kay, Estimation Theory, 1993) to within 3 percent. The calculation also frees +the decay and constant. The legacy `step` argument remains accepted for +compatibility; analytic derivatives do not use a finite-difference step. ## Mode matching and statistics @@ -68,8 +72,9 @@ the error count so estimator failures can be distinguished from mode loss. For each mode and parameter (`amplitude`, `decay`, `angular_frequency`, `phase`) the summary gives `bias`, `std` (ddof 1) and `rmse` computed over the recovered trials only, the bound (`crlb_std`, `crlb_variance`) and -`std_over_crlb_std`. Truth and fitted modes use nonnegative amplitudes, -folding a negative amplitude into the phase. Both phases and phase errors +`std_over_crlb_std`. Truth and fitted modes use nonnegative amplitudes and +frequencies. A negative frequency is reflected along with its phase; a +negative amplitude is folded into the phase. Both phases and phase errors are wrapped to `[-pi, pi)`; the summary records truth in this convention. Undefined statistics are written as JSON `null`: standard deviation requires at least two recovered trials, and bias and RMSE require at least one. diff --git a/src/cuphoton/xray/linear_prediction.py b/src/cuphoton/xray/linear_prediction.py index 6c0f56b4..d93bd85e 100644 --- a/src/cuphoton/xray/linear_prediction.py +++ b/src/cuphoton/xray/linear_prediction.py @@ -712,6 +712,7 @@ def amplitudes_and_phases_cupy(a0: Any, a1: Any): "LinearPredictionVariableArtifactsBenchmark", "LinearPredictionRuntimeBridgeBenchmark", "ModelOrderSweep", + "ModelOrderSweepBatch", "ModelOrderSweepEntry", "PredictionRootsBackendBenchmark", "PredictionRootsBenchmark", @@ -740,6 +741,7 @@ def amplitudes_and_phases_cupy(a0: Any, a1: Any): "linear_prediction_cupy", "linear_prediction_numpy", "model_order_sweep", + "model_order_sweep_batch", "synthetic_prediction_coefficients", "synthetic_trace", "synthetic_trace_batch", diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index 8d43e1e4..b8d9f26d 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -92,6 +92,24 @@ def modes_to_theta( ) +def _modes_jacobian( + theta: np.ndarray, time: np.ndarray, n_modes: int +) -> np.ndarray: + """Analytic model derivatives in the order used by ``modes_to_theta``.""" + jac = np.ones((time.size, theta.size)) + for k in range(n_modes): + amplitude, decay, frequency, phase = theta[4 * k : 4 * k + 4] + envelope = np.exp(-decay * time) + angle = frequency * time + phase + cosine = envelope * np.cos(angle) + sine = envelope * np.sin(angle) + jac[:, 4 * k] = cosine + jac[:, 4 * k + 1] = -amplitude * time * cosine + jac[:, 4 * k + 2] = -amplitude * time * sine + jac[:, 4 * k + 3] = -amplitude * sine + return jac + + def synthetic_modes_trace( samples: int = 96, *, @@ -195,10 +213,13 @@ def cramer_rao_bounds( """Cramer-Rao lower bound (standard deviation) per parameter for white Gaussian noise of ``noise_sigma``. - The Fisher information is J^T J / sigma^2 with J the numerical Jacobian - of :func:`modes_model`; the bound is the square root of the diagonal of - its inverse. Returned per mode as arrays ``amplitude``, ``decay``, - ``angular_frequency``, ``phase`` plus the scalar ``constant``. + The Fisher information is J^T J / sigma^2 with J the analytic Jacobian + of :func:`modes_model`. A column-scaled SVD gives the bounds without + forming or inverting the normal matrix. Models with numerically + unidentifiable parameters raise ``ValueError``. Returned per mode as + arrays ``amplitude``, ``decay``, ``angular_frequency``, ``phase`` plus + the scalar ``constant``. ``step`` remains accepted for compatibility + but is unused by the analytic calculation. """ if not np.isfinite(noise_sigma) or noise_sigma <= 0: raise ValueError("noise_sigma must be finite and positive") @@ -209,27 +230,39 @@ def cramer_rao_bounds( if not np.isfinite(step) or step <= 0: raise ValueError("step must be finite and positive") theta = modes_to_theta(modes, constant) + time = np.asarray(time, dtype=np.float64) + if time.ndim != 1 or not np.all(np.isfinite(time)): + raise ValueError("time must be a finite one-dimensional array") + if not np.all(np.isfinite(theta)): + raise ValueError("modes and constant must be finite") n = len(modes) - jac = np.empty((time.size, theta.size)) - for k in range(theta.size): - tp = theta.copy() - tm = theta.copy() - tp[k] += step - tm[k] -= step - jac[:, k] = (modes_model(tp, time, n) - modes_model(tm, time, n)) / ( - 2 * step - ) try: with np.errstate(over="raise", divide="raise", invalid="raise"): - fisher = jac.T @ jac / noise_variance - sd = np.sqrt(np.diag(np.linalg.inv(fisher))) + jac = _modes_jacobian(theta, time, n) + scales = np.linalg.norm(jac, axis=0) + if time.size < theta.size or np.any(scales == 0): + raise ValueError( + "Cramer-Rao bounds require identifiable modes" + ) + _, singular, vh = np.linalg.svd(jac / scales, full_matrices=False) + cutoff = np.finfo(float).eps * max(jac.shape) * singular[0] + if singular[-1] <= cutoff: + raise ValueError( + "Cramer-Rao bounds require identifiable modes" + ) + sd = noise_sigma * np.linalg.norm(vh / singular[:, None], axis=0) + sd /= scales except (FloatingPointError, np.linalg.LinAlgError) as exc: raise ValueError( "noise_sigma or modes produce singular or nonfinite " "Cramer-Rao bounds" ) from exc - if not np.all(np.isfinite(sd)) or np.any(sd <= 0): - raise ValueError("Cramer-Rao bounds must be finite and positive") + with np.errstate(over="ignore", under="ignore"): + variance = np.square(sd) + if not np.all(np.isfinite(variance)) or np.any(variance <= 0): + raise ValueError( + "noise_sigma or modes produce nonfinite or zero Cramer-Rao bounds" + ) return { "amplitude": sd[0 : 4 * n : 4], "decay": sd[1 : 4 * n : 4], @@ -263,8 +296,8 @@ def match_modes( if not modes or angular_frequency.size == 0: return out err = np.abs( - np.array([m.angular_frequency for m in modes])[:, None] - - angular_frequency[None, :] + np.abs([m.angular_frequency for m in modes])[:, None] + - np.abs(angular_frequency)[None, :] ) # A missing match costs more than every valid distance combined: # maximize recovered modes first, then minimize total frequency error. @@ -419,14 +452,20 @@ def _wrap_phase(values: np.ndarray) -> np.ndarray: return (values + np.pi) % (2.0 * np.pi) - np.pi -def _canonical(amplitude: np.ndarray, phase: np.ndarray): - """A negative amplitude is the same mode with the phase shifted by pi.""" +def _canonical( + amplitude: np.ndarray, angular_frequency: np.ndarray, phase: np.ndarray +): + """Normalize amplitude/frequency signs and wrap equivalent phases.""" amplitude = np.array(amplitude, dtype=float, copy=True) + angular_frequency = np.array(angular_frequency, dtype=float, copy=True) phase = np.array(phase, dtype=float, copy=True) + negative_frequency = angular_frequency < 0 + angular_frequency[negative_frequency] *= -1 + phase[negative_frequency] *= -1 neg = amplitude < 0 amplitude[neg] = -amplitude[neg] phase[neg] = phase[neg] + np.pi - return amplitude, _wrap_phase(phase) + return amplitude, angular_frequency, _wrap_phase(phase) def _stats( @@ -486,9 +525,9 @@ def validation_sweep( trial in which the estimator raises is counted in ``estimator_errors`` and as a loss of every mode. - Truth and fitted modes use nonnegative amplitudes and phases in - ``[-pi, pi)``. The returned truth uses the same convention as the - parameter statistics. + Truth and fitted modes use nonnegative amplitudes and frequencies, + with equivalent phases in ``[-pi, pi)``. The returned truth uses the + same convention as the parameter statistics. ``distortion`` = (kind, amount) replaces the clean trace with one from :func:`distort_trace`, outside the model class; ``residual_ratio`` is @@ -506,13 +545,18 @@ def validation_sweep( raise ValueError("snr_db must contain at least one level") if not np.all(np.isfinite(snr_db)): raise ValueError("snr_db levels must be finite") - amplitudes, phases = _canonical( + amplitudes, frequencies, phases = _canonical( np.array([m.amplitude for m in modes]), + np.array([m.angular_frequency for m in modes]), np.array([m.phase for m in modes]), ) modes = tuple( - replace(m, amplitude=float(a), phase=float(p)) - for m, a, p in zip(modes, amplitudes, phases, strict=True) + replace( + m, amplitude=float(a), angular_frequency=float(w), phase=float(p) + ) + for m, a, w, p in zip( + modes, amplitudes, frequencies, phases, strict=True + ) ) est = estimator or linear_prediction_numpy base = synthetic_modes_trace( @@ -553,9 +597,12 @@ def validation_sweep( any_lost += 1 lost += 1 continue - w = _to_host(r.angular_frequency) d = _to_host(r.decay) - a, p = _canonical(_to_host(r.amplitude), _to_host(r.phase)) + a, w, p = _canonical( + _to_host(r.amplitude), + _to_host(r.angular_frequency), + _to_host(r.phase), + ) rec = _to_host(r.reconstruction) ratios.append(float(np.sqrt(np.mean((y - rec) ** 2)) / sigma)) hit = match_modes(modes, w, d, tolerance=match_tolerance) @@ -781,7 +828,7 @@ def build_summary( "bias, std (ddof 1) and rmse are computed over the recovered " "trials of each mode; crlb_std and crlb_variance are the " "unconditional Cramer-Rao bounds from the model's Fisher " - "information using a numerical Jacobian at the true " + "information using an analytic Jacobian at the true " "parameters; the phase " "error is wrapped to [-pi, pi)" ), diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py index 6d5033c8..e18df026 100644 --- a/tests/xray/test_synthetic_validation.py +++ b/tests/xray/test_synthetic_validation.py @@ -16,10 +16,13 @@ from cuphoton.xray.synthetic_validation import ( PARAMETERS, DampedMode, + _modes_jacobian, build_summary, cramer_rao_bounds, distort_trace, match_modes, + modes_model, + modes_to_theta, synthetic_modes_trace, validate_linear_prediction, validation_sweep, @@ -58,15 +61,21 @@ def test_estimator_recovers_known_modes_exactly_at_zero_noise(): assert abs(float(r.decay[h]) - m.decay) < 1e-6 -def test_cramer_rao_bound_scales_with_noise_and_shrinks_with_trace_length(): +@pytest.mark.parametrize("sigma", [0.01, 1e-160]) +def test_cramer_rao_bound_scales_with_noise_and_shrinks_with_trace_length( + sigma, +): fx = synthetic_modes_trace(96) - b1 = cramer_rao_bounds(fx.modes, fx.constant, fx.time, 0.01) - b2 = cramer_rao_bounds(fx.modes, fx.constant, fx.time, 0.02) + b1 = cramer_rao_bounds(fx.modes, fx.constant, fx.time, sigma) + b2 = cramer_rao_bounds(fx.modes, fx.constant, fx.time, 2 * sigma) np.testing.assert_allclose( - b2["angular_frequency"], 2 * b1["angular_frequency"], rtol=1e-6 + b2["angular_frequency"], + 2 * b1["angular_frequency"], + rtol=1e-6, + atol=0, ) longer = synthetic_modes_trace(192, duration=19.0) - b3 = cramer_rao_bounds(longer.modes, longer.constant, longer.time, 0.01) + b3 = cramer_rao_bounds(longer.modes, longer.constant, longer.time, sigma) assert np.all(b3["angular_frequency"] < b1["angular_frequency"]) @@ -102,19 +111,87 @@ def test_integer_mode_parameters_have_the_same_bounds_as_floats(): assert np.all(np.isfinite(actual[name])) +def test_analytic_jacobian_and_bounds_match_finite_differences(): + fixture = synthetic_modes_trace() + theta = modes_to_theta(fixture.modes, fixture.constant) + n_modes = len(fixture.modes) + step = 1e-6 + numerical = np.column_stack( + [ + ( + modes_model(theta + delta, fixture.time, n_modes) + - modes_model(theta - delta, fixture.time, n_modes) + ) + / (2 * step) + for delta in step * np.eye(theta.size) + ] + ) + np.testing.assert_allclose( + _modes_jacobian(theta, fixture.time, n_modes), + numerical, + rtol=1e-7, + atol=3e-10, + ) + expected = 0.01 * np.sqrt(np.diag(np.linalg.inv(numerical.T @ numerical))) + actual = cramer_rao_bounds( + fixture.modes, fixture.constant, fixture.time, 0.01 + ) + for index, name in enumerate(PARAMETERS): + np.testing.assert_allclose( + actual[name], expected[index:-1:4], rtol=1e-8 + ) + np.testing.assert_allclose(actual["constant"], expected[-1], rtol=1e-8) + + +@pytest.mark.parametrize("time_scale", [1e-9, 1.0, 1e9]) +def test_bounds_resolve_nearby_modes_in_different_time_units(time_scale): + modes = (DampedMode(1, 0.01, 2, 0.3), DampedMode(0.5, 0.01, 2.02, 1.5)) + time = np.linspace(0, 100, 512) + expected = cramer_rao_bounds(modes, 0.15, time, 0.01) + scaled = tuple( + replace( + m, + decay=m.decay / time_scale, + angular_frequency=m.angular_frequency / time_scale, + ) + for m in modes + ) + actual = cramer_rao_bounds(scaled, 0.15, time * time_scale, 0.01) + for name in (*PARAMETERS, "constant"): + factor = time_scale if name in ("decay", "angular_frequency") else 1 + np.testing.assert_allclose( + actual[name] * factor, expected[name], rtol=1e-8 + ) + assert np.all(np.isfinite(actual[name])) + assert np.all(actual[name] > 0) + + @pytest.mark.parametrize( "modes", [ (DampedMode(0.0, 0.1, 2.0, 0.3),), (DampedMode(1.0, 0.1, 2.0, 0.3),) * 2, + (DampedMode(1.0, 0.1, 2.0, 0.3), DampedMode(0.5, 0.1, 2.0, 1.5)), + ], + ids=[ + "zero-amplitude", + "duplicate-modes", + "coincident-phase-distinct-modes", ], - ids=["zero-amplitude", "duplicate-modes"], ) def test_bounds_reject_unidentifiable_modes(modes): with pytest.raises(ValueError, match="Cramer-Rao bounds"): cramer_rao_bounds(modes, 0.0, np.linspace(0.0, 9.5, 96), 0.01) +def test_bounds_reject_more_parameters_than_samples(): + fixture = synthetic_modes_trace() + with pytest.raises(ValueError, match="Cramer-Rao bounds"): + cramer_rao_bounds( + fixture.modes, fixture.constant, fixture.time[:8], 0.01 + ) + + @pytest.mark.parametrize("start", [0.0, 3.0]) def test_chirp_integrates_the_requested_frequency_ramp(start): time = np.linspace(start, start + 10.0, 256) @@ -134,9 +211,12 @@ def test_gaussian_envelope_starts_at_the_first_sample(start): np.testing.assert_allclose(actual, expected) -def test_matching_recovers_close_modes_independent_of_truth_order(): - modes = (DampedMode(1, 0, 1.0), DampedMode(1, 0, 1.15)) - fitted = np.array([1.05, 0.85]) +@pytest.mark.parametrize("frequency_sign", [1, -1]) +def test_matching_recovers_close_modes_independent_of_truth_order( + frequency_sign, +): + modes = (DampedMode(1, 0, frequency_sign * 1.0), DampedMode(1, 0, 1.15)) + fitted = np.array([1.05, frequency_sign * 0.85]) decay = np.zeros(2) assert match_modes(modes, fitted, decay, tolerance=0.2) == [1, 0] assert match_modes(modes[::-1], fitted, decay, tolerance=0.2) == [0, 1] @@ -189,18 +269,21 @@ def broken(t, y, k): @pytest.mark.parametrize("amplitude", [-1.0, 1.0]) @pytest.mark.parametrize("phase", [0.3, 0.3 + 4 * np.pi]) +@pytest.mark.parametrize("frequency_sign", [1, -1]) def test_sweep_uses_the_same_amplitude_phase_convention_for_truth_and_fit( - amplitude, phase + amplitude, phase, frequency_sign ): - mode = DampedMode(amplitude, 0.1, 2.0, phase) + mode = DampedMode( + amplitude, 0.1, frequency_sign * 2.0, frequency_sign * phase + ) fixture = synthetic_modes_trace(modes=(mode,)) def exact_estimator(t, y, k): return SimpleNamespace( amplitude=np.array([mode.amplitude]), decay=np.array([mode.decay]), - angular_frequency=np.array([mode.angular_frequency]), - phase=np.array([mode.phase]), + angular_frequency=np.array([-mode.angular_frequency]), + phase=np.array([-mode.phase]), reconstruction=fixture.clean, ) @@ -224,12 +307,27 @@ def exact_estimator(t, y, k): summary = build_summary([sweep]) truth = summary["config"]["truth"]["modes"][0] assert truth["amplitude"] == 1.0 + assert truth["angular_frequency"] == 2.0 expected_phase = 0.3 - np.pi if amplitude < 0 else 0.3 assert truth["phase"] == pytest.approx(expected_phase) for name in PARAMETERS: assert summary["results"][0]["modes"][0][name]["truth"] == truth[name] +def test_equivalent_frequency_signs_have_the_same_mode_recovery(): + positive = DampedMode(1, 0.1, 2, 0.3) + negative = replace(positive, angular_frequency=-2, phase=-0.3) + sweeps = [ + validation_sweep( + modes=(mode,), snr_db=(80,), trials=2, n_components=4 + ) + for mode in (positive, negative) + ] + assert sweeps[0].levels == sweeps[1].levels + assert sweeps[0].modes == sweeps[1].modes + assert sweeps[0].levels[0].modes[0].recovered_trials == 2 + + def test_summary_schema_and_run_artifacts(tmp_path): sweep = validation_sweep(snr_db=(30.0,), trials=6, seed=3) summary = build_summary([sweep]) @@ -364,9 +462,7 @@ def test_fixture_rejects_invalid_noise_sigma(sigma): synthetic_modes_trace(noise_sigma=sigma) -@pytest.mark.parametrize( - "sigma", [np.nan, np.inf, -0.1, 1e-200, 1e-160, 1e200] -) +@pytest.mark.parametrize("sigma", [np.nan, np.inf, -0.1, 1e-200, 1e200]) def test_bounds_reject_invalid_noise_scales_without_runtime_warnings(sigma): fx = synthetic_modes_trace() with ( From b7282ede7a9ef93e24c9c0894183c0b2a730b450 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Mon, 12 Oct 2026 13:09:18 -0700 Subject: [PATCH 10/10] Use decay to break XRay mode-matching ties Keep recovered-mode count and frequency error ahead of decay error. Compare assignment costs exactly so decay only resolves frequency ties, including reversed modes and extra fitted candidates. Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-VALIDATION.md | 3 +- src/cuphoton/xray/synthetic_validation.py | 119 +++++++++++++++++++--- tests/xray/test_synthetic_validation.py | 98 ++++++++++++++++++ 3 files changed, 203 insertions(+), 17 deletions(-) diff --git a/docs/xray/LINEAR-PREDICTION-VALIDATION.md b/docs/xray/LINEAR-PREDICTION-VALIDATION.md index 614c6fb8..a3a8e369 100644 --- a/docs/xray/LINEAR-PREDICTION-VALIDATION.md +++ b/docs/xray/LINEAR-PREDICTION-VALIDATION.md @@ -62,7 +62,8 @@ compatibility; analytic derivatives do not use a finite-difference step. A true mode is recovered in a trial when a fitted mode has an angular frequency within `match_tolerance` (default 0.3 rad per time unit) of it; the assignment first maximizes the number of recovered modes, then minimizes -their total frequency error. Each fitted mode is used at most once. The +their total frequency error, using total decay error to break ties. +Each fitted mode is used at most once. The per-mode `loss_rate` is the fraction of trials in which the mode was not recovered. A trial is successful when the estimator returned and every true mode was recovered; a trial in which the estimator raised counts in diff --git a/src/cuphoton/xray/synthetic_validation.py b/src/cuphoton/xray/synthetic_validation.py index b8d9f26d..778de29b 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -24,7 +24,6 @@ from typing import Any import numpy as np -from scipy.optimize import linear_sum_assignment from cuphoton.core.runtime import runtime_metadata @@ -280,6 +279,75 @@ def _to_host(values) -> np.ndarray: return np.asarray(values, dtype=float) +def _integer_errors(errors: np.ndarray, n_modes: int) -> list[int]: + """Scale binary floats exactly, ranking nonfinite errors last.""" + ratios = [ + float(error).as_integer_ratio() if np.isfinite(error) else None + for error in errors + ] + denominator = max((r[1] for r in ratios if r is not None), default=1) + scaled = [ + r[0] * (denominator // r[1]) if r is not None else None + for r in ratios + ] + unavailable = ( + n_modes * max((x for x in scaled if x is not None), default=0) + 1 + ) + return [unavailable if x is None else x for x in scaled] + + +def _minimum_mode_assignment(cost: list[list[int]]) -> list[int]: + """Rectangular Hungarian assignment without rounding the integer costs.""" + n_rows, n_columns = len(cost), len(cost[0]) + row_price = [0] * (n_rows + 1) + column_price = [0] * (n_columns + 1) + owner = [0] * (n_columns + 1) + previous = [0] * (n_columns + 1) + for row in range(1, n_rows + 1): + owner[0] = row + column = 0 + slack: list[int | None] = [None] * (n_columns + 1) + visited = [False] * (n_columns + 1) + while True: + visited[column] = True + current_row = owner[column] + step = None + next_column = 0 + for candidate in range(1, n_columns + 1): + if visited[candidate]: + continue + reduced = ( + cost[current_row - 1][candidate - 1] + - row_price[current_row] + - column_price[candidate] + ) + best = slack[candidate] + if best is None or reduced < best: + best = slack[candidate] = reduced + previous[candidate] = column + if step is None or best < step: + step, next_column = best, candidate + assert step is not None + for candidate, best in enumerate(slack): + if visited[candidate]: + row_price[owner[candidate]] += step + column_price[candidate] -= step + elif best is not None: + slack[candidate] = best - step + column = next_column + if owner[column] == 0: + break + while column: + predecessor = previous[column] + owner[column] = owner[predecessor] + column = predecessor + assignment = [0] * n_rows + for column, row in enumerate(owner[1:]): + if row: + assignment[row - 1] = column + return assignment + + def match_modes( modes: Sequence[DampedMode], angular_frequency: np.ndarray, @@ -287,8 +355,11 @@ def match_modes( *, tolerance: float, ) -> list[int | None]: - """For each true mode, the index of the fitted mode within ``tolerance`` - (rad per time unit) or ``None``. + """Match within the frequency ``tolerance`` (rad per time unit). + + Maximize recovered modes, then minimize total frequency error, then + total decay error. Return fitted indices or ``None`` for missing modes. + Nonfinite decay errors rank last without preventing frequency recovery. """ if not np.isfinite(tolerance) or tolerance < 0: raise ValueError("tolerance must be finite and nonnegative") @@ -299,19 +370,35 @@ def match_modes( np.abs([m.angular_frequency for m in modes])[:, None] - np.abs(angular_frequency)[None, :] ) - # A missing match costs more than every valid distance combined: - # maximize recovered modes first, then minimize total frequency error. - missing_cost = len(modes) + 1 - cost = np.full( - (len(modes), angular_frequency.size + len(modes)), float(missing_cost) - ) - cost[:, : angular_frequency.size] = np.where( - err <= tolerance, err / (tolerance or 1.0), np.inf - ) - rows, columns = linear_sum_assignment(cost) - for row, column in zip(rows, columns, strict=True): + valid = err <= tolerance + if not np.any(valid): + return out + with np.errstate(over="ignore", invalid="ignore"): + decay_err = np.abs( + np.array([m.decay for m in modes])[:, None] - decay[None, :] + ) + frequency_cost = _integer_errors(err[valid], len(modes)) + decay_cost = _integer_errors(decay_err[valid], len(modes)) + # Integer bounds make each priority dominate every possible total of + # the lower priorities. No epsilon lets decay outweigh frequency. + frequency_scale = len(modes) * max(decay_cost) + 1 + real_cost = [ + frequency * frequency_scale + damping + for frequency, damping in zip(frequency_cost, decay_cost, strict=True) + ] + missing_cost = len(modes) * max(real_cost) + 1 + invalid_cost = (len(modes) + 1) * missing_cost + cost = [ + [invalid_cost] * angular_frequency.size + [missing_cost] * len(modes) + for _ in modes + ] + for (row, column), value in zip( + np.argwhere(valid), real_cost, strict=True + ): + cost[row][column] = value + for row, column in enumerate(_minimum_mode_assignment(cost)): if column < angular_frequency.size: - out[int(row)] = int(column) + out[row] = column return out @@ -820,7 +907,7 @@ def build_summary( "a true mode is recovered when a fitted mode lies within " "the tolerance of its angular frequency; fitted modes are " "assigned to maximize matches, then minimize total frequency " - "error, and each is used at most once" + "error, then total decay error, and each is used at most once" ), "tolerance_rad_per_time": first.match_tolerance, }, diff --git a/tests/xray/test_synthetic_validation.py b/tests/xray/test_synthetic_validation.py index e18df026..13331132 100644 --- a/tests/xray/test_synthetic_validation.py +++ b/tests/xray/test_synthetic_validation.py @@ -229,6 +229,104 @@ def test_matching_recovers_close_modes_independent_of_truth_order( ) == [None, None] +@pytest.mark.parametrize( + "truth_frequency,truth_decay,fitted_frequency,fitted_decay,tolerance,expected", + [ + # Different frequencies can also have equal total frequency error. + ([1, 2], [0.1, 0.8], [3, 4], [0.8, 0.1], 3, [1, 0]), + # Include fitted modes left unused by the primary assignment. + ([2], [0.1], [1, 3], [0.8, 0.1], 1, [1]), + # A partial recovery goes to the true mode with the matching decay. + ([2, 2], [0.1, 0.8], [2], [0.8], 0, [None, 0]), + # Better decay must not sacrifice an otherwise recoverable mode. + ([1, 1.15], [0, 100], [1.05, 0.85], [0, 100], 0.2, [1, 0]), + # Nor may it outweigh even one representable step in frequency error. + ([1], [0], [np.nextafter(1.0, 2.0), 1], [0, 1e100], 1, [1]), + # Rounded floating-point totals hide this genuine frequency advantage. + ( + [1, 2, 10], + [0, 1, 2], + [1.5, np.nextafter(1.5, 2.0), 1000], + [1, 0, 2], + 1000, + [0, 1, 2], + ), + ], +) +def test_matching_uses_decay_only_after_recovery_and_frequency( + truth_frequency, + truth_decay, + fitted_frequency, + fitted_decay, + tolerance, + expected, +): + modes = tuple( + DampedMode(1, d, w) + for w, d in zip(truth_frequency, truth_decay, strict=True) + ) + for reverse in (False, True): + order = slice(None, None, -1 if reverse else 1) + actual = match_modes( + modes, + np.array(fitted_frequency)[order], + np.array(fitted_decay)[order], + tolerance=tolerance, + ) + assert actual == [ + len(fitted_frequency) - 1 - h if reverse and h is not None else h + for h in expected + ] + + +@pytest.mark.parametrize("unavailable_decay", [np.nan, np.inf, -np.inf]) +def test_matching_nonfinite_decay_does_not_prevent_frequency_recovery( + unavailable_decay, +): + modes = (DampedMode(1, 0.1, 2),) + assert match_modes( + modes, np.array([2]), np.array([unavailable_decay]), tolerance=0 + ) == [0] + assert match_modes( + modes, + np.array([2, 2]), + np.array([unavailable_decay, 0.1]), + tolerance=0, + ) == [1] + + +@pytest.mark.parametrize("reverse", [False, True]) +def test_sweep_matches_exact_modes_with_shared_frequency_by_decay(reverse): + modes = (DampedMode(1, 0.05, 2, 0.3), DampedMode(0.4, 0.5, 2, 1)) + fitted_modes = modes[::-1] if reverse else modes + fixture = synthetic_modes_trace(modes=modes) + + def exact_estimator(t, y, k): + return SimpleNamespace( + **{ + name: np.array([getattr(m, name) for m in fitted_modes]) + for name in PARAMETERS + }, + reconstruction=fixture.clean, + ) + + level = validation_sweep( + modes=modes, + snr_db=(80,), + trials=4, + seed=314, + estimator=exact_estimator, + ).levels[0] + assert level.trials_successful == 4 + for mode in level.modes: + assert mode.recovered_trials == 4 + assert mode.loss_rate == 0 + for name in PARAMETERS: + stats = getattr(mode, name) + assert stats.bias == pytest.approx(0, abs=1e-14) + assert stats.rmse == pytest.approx(0, abs=1e-14) + + def test_validation_sweep_reports_bound_ratios_and_loss_rates(): sweep = validation_sweep(snr_db=(40.0, 20.0), trials=40, seed=1) assert len(sweep.levels) == 2