From 314c543392156ea299e9c511df6c9b12b1a3180a Mon Sep 17 00:00:00 2001 From: James Sweeney Date: Tue, 15 Sep 2026 18:07:59 -0700 Subject: [PATCH 1/9] Add synthetic damped-mode validation for XRay 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 --- 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 f94b0102..b825c837 100644 --- a/docs/xray/README.md +++ b/docs/xray/README.md @@ -30,6 +30,7 @@ uv run python examples/run_quickstarts.py --component xray --profile cpu | 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` | @@ -206,3 +207,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 3055a8a3..682b375d 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" @@ -4320,6 +4398,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 8a575e7c..1b1b568f 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, 144, 1), ("xscan", 42, 42, 171, 1), ("xrep", 6, 6, 103, 1), - ("xray", 33, 31, 390, 1), + ("xray", 34, 32, 398, 1), ] assert len(per_group) == 6 - assert sum(item[1] for item in per_group) == 92 - assert sum(item[1] + item[4] for item in per_group) == 97 - assert sum(item[2] for item in per_group) == 89 - assert sum(item[3] for item in per_group) == 839 + assert sum(item[1] for item in per_group) == 93 + assert sum(item[1] + item[4] for item in per_group) == 98 + assert sum(item[2] for item in per_group) == 90 + assert sum(item[3] for item in per_group) == 847 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 d4da117a4cb71fe561964e2e9aa465a885698e41 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Fri, 18 Sep 2026 20:29:11 -0700 Subject: [PATCH 2/9] 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 682b375d..987aa686 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -4450,7 +4450,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 d39bfb59db939461d43ad102a96c6baa1efcb0d2 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Wed, 23 Sep 2026 19:43:50 -0700 Subject: [PATCH 3/9] 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 | 14 +++++++- tests/xray/test_synthetic_validation.py | 44 +++++++++++++++++++++++ 3 files changed, 60 insertions(+), 3 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..fedcb9cd 100644 --- a/src/cuphoton/xray/synthetic_validation.py +++ b/src/cuphoton/xray/synthetic_validation.py @@ -18,7 +18,7 @@ import json import subprocess -from dataclasses import asdict, dataclass, field +from dataclasses import asdict, dataclass, field, replace from pathlib import Path from typing import Any, Callable, Sequence @@ -403,6 +403,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 +417,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) + ) 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 aa33a44eaf278f88a017cf90dd7cbb33b758a6bb Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Wed, 23 Sep 2026 19:44:07 -0700 Subject: [PATCH 4/9] 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 987aa686..853e37c0 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -26,6 +26,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, @@ -3746,8 +3753,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 @@ -4297,25 +4304,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, @@ -4398,31 +4386,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, @@ -4455,7 +4418,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'}") @@ -5489,9 +5452,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 @@ -5569,62 +5532,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 @@ -5692,7 +5599,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']}") @@ -5822,32 +5729,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 d0670d24400f81708d24f9453b139679901d1998 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Wed, 23 Sep 2026 20:19:37 -0700 Subject: [PATCH 5/9] 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 853e37c0..2438a52a 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -20,6 +20,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 @@ -4396,8 +4397,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 fedcb9cd..1a6f0e3d 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. @@ -107,6 +107,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 @@ -138,7 +142,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: @@ -147,6 +152,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): @@ -160,7 +169,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) @@ -189,8 +199,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)) @@ -202,8 +218,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], @@ -415,8 +439,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]), @@ -425,7 +455,7 @@ def validation_sweep( replace(m, amplitude=float(a), phase=float(p)) for m, a, p in zip(modes, amplitudes, phases) ) - 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 ) @@ -435,11 +465,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 ] @@ -523,19 +559,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( @@ -630,8 +685,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}), @@ -650,7 +706,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: @@ -716,8 +772,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"] @@ -782,8 +840,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 d2b4eceae0059a7f71d375b42704e607b42d99fb Mon Sep 17 00:00:00 2001 From: James Sweeney Date: Tue, 15 Sep 2026 18:15:15 -0700 Subject: [PATCH 6/9] Add per-trace and batched nonlinear mode refinement cuphoton.xray.mode_refinement refines the damped-mode model by nonlinear least squares from the linear-prediction result (analytic Jacobian, scipy.optimize.least_squares), seeding a mode the root filter dropped from the residual spectrum; RefinedModes carries Jacobian uncertainties and the residual rms. cuphoton.xray.mode_refinement_batched runs the same Levenberg-Marquardt for a batch of traces on NumPy or CuPy and matches the per-trace solver to about 1e-8. linear-prediction-validate gains --refine, which adds the refined estimator to the sweep and the summary; linear-prediction-refine-benchmark (lprb) times the SciPy loop against the batched solver on NumPy and CuPy (best of --repeat, transfers excluded, warm-up call). examples/xray_lp_modes_review.py writes a Bokeh review page of traces, reconstructions and modes. Documented in docs/xray/LINEAR-PREDICTION-REFINEMENT.md with RTX 4050 numbers. Overlaps the fit-diagnostics, iterative-fitting and GPU-batching work the maintainers describe in #1; posted so the overlap can be sorted out. Refs #1 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_014wyj5iPSHmtR3vSTvJWKfw Signed-off-by: James Sweeney --- docs/xray/LINEAR-PREDICTION-REFINEMENT.md | 124 ++++++++ docs/xray/README.md | 3 +- examples/xray_lp_modes_review.py | 175 +++++++++++ src/cuphoton/xray/commands.py | 106 +++++++ src/cuphoton/xray/mode_refinement.py | 226 +++++++++++++++ src/cuphoton/xray/mode_refinement_batched.py | 287 +++++++++++++++++++ tests/core/test_cli_contract.py | 10 +- tests/xray/test_cli.py | 22 ++ tests/xray/test_mode_refinement.py | 99 +++++++ tests/xray/test_mode_refinement_batched.py | 81 ++++++ 10 files changed, 1127 insertions(+), 6 deletions(-) create mode 100644 docs/xray/LINEAR-PREDICTION-REFINEMENT.md create mode 100644 examples/xray_lp_modes_review.py create mode 100644 src/cuphoton/xray/mode_refinement.py create mode 100644 src/cuphoton/xray/mode_refinement_batched.py create mode 100644 tests/xray/test_mode_refinement.py create mode 100644 tests/xray/test_mode_refinement_batched.py diff --git a/docs/xray/LINEAR-PREDICTION-REFINEMENT.md b/docs/xray/LINEAR-PREDICTION-REFINEMENT.md new file mode 100644 index 00000000..7c1ff59a --- /dev/null +++ b/docs/xray/LINEAR-PREDICTION-REFINEMENT.md @@ -0,0 +1,124 @@ +# Nonlinear refinement of linear-prediction modes + +Linear prediction gives frequencies and decays without an initial guess, but +it is not a maximum-likelihood estimator and its root filter can drop a +lightly damped mode under noise. `cuphoton.xray.mode_refinement` refines the +same damped-mode model by nonlinear least squares from the linear-prediction +result; `cuphoton.xray.mode_refinement_batched` does the same for a batch of +traces at once on NumPy or CuPy. Both fit + +``` +trace(t) = constant + sum_k amplitude_k exp(-decay_k t) cos(angular_frequency_k t + phase_k) +``` + +with an analytic Jacobian. Nothing here changes the linear-prediction +commands or their defaults. + +## Per-trace refinement + +`linear_prediction_refined(time, trace, n_components, n_modes)` runs +`linear_prediction_numpy`, keeps the `n_modes` strongest oscillating modes as +the starting point, and refines them and the constant with +`scipy.optimize.least_squares` (Levenberg-Marquardt, `xtol` and `ftol` +1e-12). When the linear-prediction step returned fewer than `n_modes` +oscillating modes, the missing ones are seeded one at a time from the +strongest peak of the residual spectrum (amplitude from the peak height, zero +decay and phase) and refined before the next seed. `refine_modes` is the +same solver from caller-supplied starting modes. + +`RefinedModes` carries the refined parameters ordered by amplitude, one-sigma +uncertainties from the Jacobian at the solution and the residual variance, +the reconstruction, the residual rms, the solver status, and how many modes +came from linear prediction (`initial_mode_count`) versus the residual +spectrum (`seeded_mode_count`). A residual rms well above the noise means +the model class, not the fit, is wrong. + +`linear-prediction-validate --refine` runs the validation sweep a second time +with this estimator (`estimator: cpu-refined` in the summary) so the two can +be compared level by level. Reference run on the default fixture, CPU, +default seed, 300 trials per level: + +```bash +uv run cuphoton xray lpv --refine --trials 300 --snr-db 40,30,20,10,5 +``` + +In that run the linear-prediction estimator loses the weak mode in 10 to 52 +percent of trials depending on the level, and its frequency scatter is 1.1 +to 3.8 times the bound. The refined estimator loses the weak mode in 1 of +300 trials at each of the two lowest levels of that run and in none above, +and its frequency and decay scatter are within 10 percent of the Cramer-Rao +bound at every level of that run (ratios 0.91 to 1.07; the sampling +uncertainty of a standard deviation from 300 trials is about 4 percent). +The refinement does not remove the bias a model mismatch produces, and the +residual ratio still flags that case. + +```bash +uv run cuphoton xray lpv --refine --snr-db 40,20,10,5 --trials 200 --output-dir /tmp/lpv-refined +uv run cuphoton xray lpv --refine --distortion chirp:0.05 --snr-db 30 +uv run python examples/xray_lp_modes_review.py --output /tmp/review.html +``` + +The example writes a Bokeh review page (requires the `viz` extra) with the +trace, the true model, the linear-prediction and refined reconstructions, +and the recovered modes against the true ones, for several noise levels and +an optional distortion. + +## Batched refinement + +`refine_modes_batched(xp, time, traces, theta0, n_modes)` takes traces shaped +`(B, N)` and starting parameters `(B, 4K+1)` in the same layout as the +per-trace solver (amplitude, decay, angular frequency, phase per mode, then +the constant) and runs Levenberg-Marquardt for the whole batch: the Jacobian +for all traces by broadcasting, the normal equations as a batch of small +`(4K+1)` systems solved with `xp.linalg.solve`, and per-trace damping (a +step that lowers the cost is accepted and the damping divided by 3, otherwise +the damping is multiplied by 5 and the trace keeps its parameters). Iteration +stops when every trace has converged (step below `tol` times the parameter +norm, default 1e-9) or after `max_iter` (default 60). `xp` is `numpy` or +`cupy`. On convergence the result matches the per-trace SciPy solver to about +1e-8 in every parameter (table below). + +The batch iterates until its slowest trace converges, so the NumPy batched +path is slower than the SciPy loop; the batched form pays off on the GPU. + +## Benchmark + +`linear-prediction-refine-benchmark` (`lprb`) builds `--traces` noisy copies +of the fixture (two modes, 96 samples, noise sigma 0.02), perturbs the true +parameters per trace to make identical starting points, and times three +paths: the SciPy solver in a Python loop over traces, the batched solver on +NumPy, and the batched solver on CuPy. Each path is timed `--repeat` times +(default 3) and the best is reported. For the batched paths the clock covers +the solver only: the arrays are on the device before the clock starts and +the device is synchronised before and after, so host-to-device transfer and +the copy of the result back are not included. One untimed CuPy call on 16 +traces runs first to absorb compilation and allocation. The command also +reports the largest parameter difference between each batched result and +the SciPy result. + +```bash +uv run cuphoton xray lprb --traces 256 --repeat 3 --json +uv run cuphoton xray lprb --traces 2048 --repeat 3 --json +uv run cuphoton xray lprb --traces 8192 --repeat 3 --json +``` + +### NVIDIA GeForce RTX 4050 Laptop GPU (Ada, compute capability 8.9, 6 GiB) + +Host: Intel Core i9-13900H, WSL2 Ubuntu 24.04 on Windows 11, driver 596.49, +CUDA runtime 13.2, CuPy 14.1.1, Python 3.12, NumPy and SciPy from the locked +`dev` and `gpu` extras. Best of 3 repeats; the three commands above. + +| traces | SciPy per trace, serial (s) | NumPy batched (s) | CuPy batched (s) | SciPy / CuPy | max abs parameter difference | +| --- | --- | --- | --- | --- | --- | +| 256 | 0.135 | 0.258 | 0.108 | 1.3 | 8.7e-9 | +| 2048 | 1.09 | 2.49 | 0.162 | 6.7 | 1.2e-8 | +| 8192 | 4.38 | 12.7 | 0.800 | 5.5 | 1.9e-8 | + +On this card the GPU path breaks even at a few hundred traces and is 5 to 7 +times faster than the SciPy loop from about two thousand traces upward. +Timings at 256 traces vary by tens of percent between runs on this laptop +(a second run of the 256-trace command gave 0.178 s for CuPy); the larger +batches are stable to a few percent. In a +`cupyx.profiler.benchmark` breakdown at 256 traces on the same card, one +iteration costs about 1.5 ms in the Jacobian (launch-bound elementwise +work), 0.2 ms in the normal equations and 0.2 ms in the batched solve. diff --git a/docs/xray/README.md b/docs/xray/README.md index b825c837..05316b13 100644 --- a/docs/xray/README.md +++ b/docs/xray/README.md @@ -30,7 +30,7 @@ uv run python examples/run_quickstarts.py --component xray --profile cpu | 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`) | +| Synthetic validation | `linear-prediction-validate` (`lpv`), `linear-prediction-refine-benchmark` (`lprb`) | | 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` | @@ -208,3 +208,4 @@ See [Data and artifact contracts](../data-artifacts.md#xray-hdf5-and-trace-produ - [Distributed detector artifacts](DISTRIBUTED-DETECTOR-ARTIFACTS.md) - [Validation visualization](VALIDATION-VIZ.md) - [Linear prediction validation](LINEAR-PREDICTION-VALIDATION.md) +- [Nonlinear refinement of linear-prediction modes](LINEAR-PREDICTION-REFINEMENT.md) diff --git a/examples/xray_lp_modes_review.py b/examples/xray_lp_modes_review.py new file mode 100644 index 00000000..a501d9f4 --- /dev/null +++ b/examples/xray_lp_modes_review.py @@ -0,0 +1,175 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +"""Review view for linear prediction on the synthetic damped-mode fixture. + +Writes one standalone Bokeh HTML file with, for a few noise levels and one +optional model mismatch, the trace, the true model, the linear-prediction +reconstruction, the refined reconstruction, and the recovered modes +against the true ones. No external data. Requires the ``viz`` extra. + + uv run python examples/xray_lp_modes_review.py --output review.html + uv run python examples/xray_lp_modes_review.py --distortion chirp:0.05 +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np + +from cuphoton.xray.linear_prediction import linear_prediction_numpy +from cuphoton.xray.mode_refinement import linear_prediction_refined +from cuphoton.xray.synthetic_validation import ( + distort_trace, + modes_model, + modes_to_theta, + synthetic_modes_trace, +) + + +def _mode_curve(a, d, w, p, t): + return a * np.exp(-d * t) * np.cos(w * t + p) + + +def build(output: Path, snr_db, distortion, seed: int, samples: int): + try: + from bokeh.layouts import gridplot + from bokeh.plotting import figure, output_file, save + except ImportError: # pragma: no cover - depends on the viz extra + sys.exit("bokeh is required: uv sync --extra viz") + + fx = synthetic_modes_trace(samples) + t = fx.time + clean = fx.clean + if distortion is not None: + clean = distort_trace(t, fx.modes, fx.constant, *distortion) + fine = np.linspace(t[0], t[-1], 2000) + true_fine = modes_model(modes_to_theta(fx.modes, fx.constant), fine, 2) + sig = float(np.std(clean - clean.mean())) + rng = np.random.default_rng(seed) + rows = [] + for snr in snr_db: + sigma = sig / 10 ** (snr / 20) + y = clean + rng.normal(0.0, sigma, size=t.shape) + lp = linear_prediction_numpy(t, y, 6) + rf = linear_prediction_refined(t, y, 6, len(fx.modes)) + lp_ratio = np.sqrt(np.mean((y - lp.reconstruction) ** 2)) / sigma + left = figure( + width=760, + height=280, + title=( + f"{snr:g} dB, sigma {sigma:.4g}: residual/noise " + f"lp {lp_ratio:.1f}, refined {rf.residual_rms / sigma:.1f}" + ), + x_axis_label="time", + y_axis_label="trace", + ) + left.scatter(t, y, size=4, color="black", legend_label="trace") + left.line(fine, true_fine, color="gray", legend_label="true model") + left.line( + t, + np.asarray(lp.reconstruction), + color="red", + legend_label="linear prediction", + ) + left.line( + t, + rf.reconstruction, + color="blue", + line_dash="dashed", + legend_label="refined", + ) + left.legend.location = "top_right" + left.legend.label_text_font_size = "8pt" + w = np.asarray(lp.angular_frequency) + lost = [ + m.angular_frequency + for m in fx.modes + if w.size == 0 or np.min(np.abs(w - m.angular_frequency)) > 0.3 + ] + right = figure( + width=520, + height=280, + title="modes: true (solid), linear prediction (dashed), " + "refined (dotted)" + (f"; lp lost w={lost[0]:g}" if lost else ""), + x_axis_label="time", + ) + colors = ["#1f77b4", "#2ca02c", "#ff7f0e", "#9467bd"] + for k, m in enumerate(fx.modes): + right.line( + fine, + _mode_curve( + m.amplitude, m.decay, m.angular_frequency, m.phase, fine + ), + color=colors[k % 4], + alpha=0.5, + legend_label=f"true w={m.angular_frequency:g}", + ) + for k in range(w.size): + if w[k] == 0: + continue + right.line( + fine, + _mode_curve( + float(lp.amplitude[k]), + float(lp.decay[k]), + float(w[k]), + float(lp.phase[k]), + fine, + ), + color="red", + line_dash="dashed", + legend_label=f"lp w={w[k]:.3f} d={lp.decay[k]:+.3f}", + ) + for k in range(rf.angular_frequency.size): + right.line( + fine, + _mode_curve( + float(rf.amplitude[k]), + float(rf.decay[k]), + float(rf.angular_frequency[k]), + float(rf.phase[k]), + fine, + ), + color="blue", + line_dash="dotted", + legend_label=( + f"refined w={rf.angular_frequency[k]:.3f} " + f"d={rf.decay[k]:+.3f}" + ), + ) + right.legend.label_text_font_size = "7pt" + rows.append([left, right]) + output_file(str(output), title="xray linear prediction on the fixture") + save(gridplot(rows)) + return output + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--output", default="xray_lp_modes_review.html") + parser.add_argument("--snr-db", default="40,20,10,5") + parser.add_argument( + "--distortion", + default="", + help="kind:amount, e.g. chirp:0.05 (see distort_trace)", + ) + parser.add_argument("--seed", type=int, default=11) + parser.add_argument("--samples", type=int, default=96) + args = parser.parse_args(argv) + snr = tuple(float(x) for x in args.snr_db.split(",") if x.strip()) + distortion = None + if args.distortion: + kind, _, amount = args.distortion.partition(":") + distortion = (kind.strip(), float(amount)) + out = build(Path(args.output), snr, distortion, args.seed, args.samples) + print(f"wrote {out}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/cuphoton/xray/commands.py b/src/cuphoton/xray/commands.py index 2438a52a..a881c7e7 100644 --- a/src/cuphoton/xray/commands.py +++ b/src/cuphoton/xray/commands.py @@ -664,6 +664,16 @@ class DistortionArg(StringInvariant): _required = False _default = None + refine = False + + class RefineArg(BoolInvariant): + _arg = "--refine" + _help = ( + "Also run the sweep with nonlinear least-squares refinement " + "of the linear-prediction modes." + ) + _required = False + output_dir = None class OutputDirArg(PathValueInvariant): @@ -682,6 +692,53 @@ class JsonArg(BoolInvariant): _required = False +class LinearPredictionRefineBenchmarkCommand(_XRayCommand): + _description_ = ( + "Benchmark per-trace SciPy refinement against the batched " + "Levenberg-Marquardt solver on NumPy and CuPy." + ) + _shortname_ = "lprb" + _handler_name_ = "_linear_prediction_refine_benchmark" + + samples = 96 + + class SamplesArg(IntegerInvariant): + _arg = "--samples" + _help = "Number of synthetic trace samples." + _required = False + _default = 96 + + traces = 256 + + class TracesArg(IntegerInvariant): + _arg = "--traces" + _help = "Number of traces in the batch." + _required = False + _default = 256 + + repeat = 3 + + class RepeatArg(IntegerInvariant): + _arg = "--repeat" + _help = "Timed repetitions per path; the best time is reported." + _required = False + _default = 3 + + no_gpu = False + + class NoGpuArg(BoolInvariant): + _arg = "--no-gpu" + _help = "Skip the CuPy comparison." + _required = False + + json = False + + class JsonArg(BoolInvariant): + _arg = "--json" + _help = "Emit machine-readable JSON." + _required = False + + class LinearPredictionBenchmarkCommand(_XRayCommand): _description_ = "Benchmark serial and batched linear-prediction P1." _shortname_ = "lpb" @@ -4417,6 +4474,24 @@ def _linear_prediction_validate(args): distortion=distortion, ) ] + if args.refine: + from .mode_refinement import linear_prediction_refined + + n_modes = len(sweeps[0].levels[0].modes) + sweeps.append( + validation_sweep( + samples=args.samples, + snr_db=snr, + trials=args.trials, + n_components=args.components, + seed=args.seed, + distortion=distortion, + estimator=lambda t, y, k: linear_prediction_refined( + t, y, k, n_modes + ), + backend="cpu-refined", + ) + ) if args.output_dir is not None: summary = write_validation_run(args.output_dir, sweeps) else: @@ -4434,6 +4509,37 @@ def _linear_prediction_validate(args): return 0 +def _linear_prediction_refine_benchmark(args): + from .mode_refinement_batched import benchmark_refinement + + result = benchmark_refinement( + samples=args.samples, + traces=args.traces, + run_gpu=not args.no_gpu, + repeat=args.repeat, + ) + payload = dataclass_asdict(result) + if args.json: + print(json.dumps(payload, indent=2, sort_keys=True)) + return 0 + print( + f"traces={result.traces} samples={result.samples} " + f"repeat={result.repeat}" + ) + print(f"cpu_serial_scipy_s={result.cpu_serial_scipy_s:.6g}") + print(f"numpy_batched_s={result.numpy_batched_s:.6g}") + print(f"max_abs_theta_diff_numpy={result.max_abs_theta_diff_numpy:.3g}") + if result.gpu_error: + print("gpu_status=unavailable") + print(f"gpu_error={result.gpu_error}") + elif result.cupy_batched_s is None: + print("gpu_status=skipped") + else: + print(f"cupy_batched_s={result.cupy_batched_s:.6g}") + print(f"max_abs_theta_diff_cupy={result.max_abs_theta_diff_cupy:.3g}") + return 0 + + def _linear_prediction_benchmark(args): from .linear_prediction import benchmark_linear_prediction_p1_batch diff --git a/src/cuphoton/xray/mode_refinement.py b/src/cuphoton/xray/mode_refinement.py new file mode 100644 index 00000000..be642730 --- /dev/null +++ b/src/cuphoton/xray/mode_refinement.py @@ -0,0 +1,226 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +"""Nonlinear least-squares refinement of damped modes, initialised by +linear prediction. + +Linear prediction gives frequencies and decays without an initial guess, +but it is not a maximum-likelihood estimator and its root filter can drop a +lightly damped mode under noise. Refining the same damped-mode model by +nonlinear least squares from the linear-prediction result recovers the +Cramer-Rao bound, and a mode the filter dropped can be seeded from the +peak of the residual spectrum before refinement. The model is + + trace(t) = c + sum_k A_k exp(-d_k t) cos(w_k t + p_k) + +with an analytic Jacobian; the solver is ``scipy.optimize.least_squares``. +NumPy only. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from scipy.optimize import least_squares + +from .linear_prediction import linear_prediction_numpy + + +@dataclass(frozen=True) +class RefinedModes: + """Refined damped modes and the fit diagnostics. + + Arrays are ordered by descending amplitude. ``sigma_*`` are one-sigma + uncertainties from the Jacobian at the solution and the residual + variance; ``initial_mode_count`` is how many oscillating modes the + linear-prediction step supplied and ``seeded_mode_count`` how many had + to be seeded from the residual spectrum. + """ + + time: np.ndarray + amplitude: np.ndarray + decay: np.ndarray + angular_frequency: np.ndarray + phase: np.ndarray + constant: float + sigma_amplitude: np.ndarray + sigma_decay: np.ndarray + sigma_angular_frequency: np.ndarray + sigma_phase: np.ndarray + reconstruction: np.ndarray + residual_rms: float + converged: bool + iterations: int + initial_mode_count: int + seeded_mode_count: int + + +def _model_and_jacobian(theta: np.ndarray, time: np.ndarray, n: int): + model = np.full(time.shape, float(theta[4 * n])) + jac = np.empty((time.size, 4 * n + 1)) + jac[:, 4 * n] = 1.0 + for k in range(n): + a, d, w, p = theta[4 * k : 4 * k + 4] + env = np.exp(-d * time) + arg = w * time + p + cos = np.cos(arg) + sin = np.sin(arg) + model += a * env * cos + jac[:, 4 * k] = env * cos + jac[:, 4 * k + 1] = -a * time * env * cos + jac[:, 4 * k + 2] = -a * time * env * sin + jac[:, 4 * k + 3] = -a * env * sin + return model, jac + + +def seed_from_residual( + time: np.ndarray, residual: np.ndarray, pad: int = 4096 +) -> tuple[float, float, float, float]: + """Seed one damped mode from the strongest peak of the residual + spectrum: amplitude from the peak height, zero decay, zero phase.""" + dt = float(time[1] - time[0]) + r = residual - residual.mean() + spec = np.abs(np.fft.rfft(r, pad)) + freq = 2.0 * np.pi * np.fft.rfftfreq(pad, d=dt) + k = int(np.argmax(spec[1:])) + 1 + return float(2.0 * spec[k] / time.size), 0.0, float(freq[k]), 0.0 + + +def refine_modes( + time: np.ndarray, + trace: np.ndarray, + initial_modes: list[tuple[float, float, float, float]], + constant: float, + *, + max_nfev: int = 2000, + initial_mode_count: int | None = None, + seeded_mode_count: int = 0, +) -> RefinedModes: + """Refine ``initial_modes`` (amplitude, decay, angular frequency, phase) + and the constant by nonlinear least squares.""" + time = np.asarray(time, dtype=np.float64) + trace = np.asarray(trace, dtype=np.float64) + n = len(initial_modes) + if n == 0: + raise ValueError("at least one initial mode is required") + theta0 = np.array( + [v for m in initial_modes for v in m] + [constant], dtype=np.float64 + ) + + def residual(theta): + return _model_and_jacobian(theta, time, n)[0] - trace + + def jacobian(theta): + return _model_and_jacobian(theta, time, n)[1] + + res = least_squares( + residual, + theta0, + jac=jacobian, + method="lm", + xtol=1e-12, + ftol=1e-12, + max_nfev=max_nfev, + ) + theta = res.x + model, jac = _model_and_jacobian(theta, time, n) + resid = trace - model + dof = max(time.size - theta.size, 1) + var = float(resid @ resid) / dof + try: + cov = np.linalg.inv(jac.T @ jac) * var + sigma = np.sqrt(np.clip(np.diag(cov), 0.0, None)) + except np.linalg.LinAlgError: + sigma = np.full(theta.size, np.nan) + amp = theta[0 : 4 * n : 4].copy() + dec = theta[1 : 4 * n : 4].copy() + freq = theta[2 : 4 * n : 4].copy() + ph = theta[3 : 4 * n : 4].copy() + # a negative amplitude or frequency is the same mode with a shifted phase + neg = amp < 0 + amp[neg] = -amp[neg] + ph[neg] = ph[neg] + np.pi + negf = freq < 0 + freq[negf] = -freq[negf] + ph[negf] = -ph[negf] + ph = (ph + np.pi) % (2.0 * np.pi) - np.pi + order = np.argsort(-amp) + sig = sigma.reshape(-1) + return RefinedModes( + time=time, + amplitude=amp[order], + decay=dec[order], + angular_frequency=freq[order], + phase=ph[order], + constant=float(theta[4 * n]), + sigma_amplitude=sig[0 : 4 * n : 4][order], + sigma_decay=sig[1 : 4 * n : 4][order], + sigma_angular_frequency=sig[2 : 4 * n : 4][order], + sigma_phase=sig[3 : 4 * n : 4][order], + reconstruction=model, + residual_rms=float(np.sqrt(np.mean(resid**2))), + converged=bool(res.success), + iterations=int(res.nfev), + initial_mode_count=( + n if initial_mode_count is None else initial_mode_count + ), + seeded_mode_count=seeded_mode_count, + ) + + +def linear_prediction_refined( + time, + trace, + n_components: int, + n_modes: int, + *, + roots_backend: str = "eigvals", +) -> RefinedModes: + """Linear prediction, then nonlinear refinement of the ``n_modes`` + strongest oscillating modes; modes the linear-prediction step did not + return are seeded one at a time from the residual spectrum.""" + time = np.asarray(time, dtype=np.float64) + trace = np.asarray(trace, dtype=np.float64) + lp = linear_prediction_numpy( + time, trace, n_components, roots_backend=roots_backend + ) + w = np.asarray(lp.angular_frequency, dtype=float) + d = np.asarray(lp.decay, dtype=float) + a = np.asarray(lp.amplitude, dtype=float) + p = np.asarray(lp.phase, dtype=float) + modes = [ + (float(a[i]), float(d[i]), float(w[i]), float(p[i])) + for i in range(w.size) + if w[i] > 0 + ] + modes.sort(key=lambda m: -m[0]) + modes = modes[:n_modes] + found = len(modes) + constant = float(np.mean(trace)) + residual = trace - np.asarray(lp.reconstruction, dtype=float) + seeded = 0 + while len(modes) < n_modes: + modes.append(seed_from_residual(time, residual)) + seeded += 1 + partial = refine_modes(time, trace, modes, constant) + residual = trace - partial.reconstruction + modes = [ + ( + float(partial.amplitude[i]), + float(partial.decay[i]), + float(partial.angular_frequency[i]), + float(partial.phase[i]), + ) + for i in range(len(modes)) + ] + constant = partial.constant + return refine_modes( + time, + trace, + modes, + constant, + initial_mode_count=found, + seeded_mode_count=seeded, + ) diff --git a/src/cuphoton/xray/mode_refinement_batched.py b/src/cuphoton/xray/mode_refinement_batched.py new file mode 100644 index 00000000..2dbc1d6b --- /dev/null +++ b/src/cuphoton/xray/mode_refinement_batched.py @@ -0,0 +1,287 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# +# SPDX-License-Identifier: Apache-2.0 + +"""Batched Levenberg-Marquardt refinement of damped modes on NumPy or CuPy. + +Many detector-row traces share one time axis and one model shape (K damped +modes plus a constant), so the refinement of :mod:`mode_refinement` can be +run for a whole batch at once: analytic Jacobians for all traces by +broadcasting, normal equations solved as a batch of small (4K+1) systems, +and per-trace damping with accept/reject. The array module is a parameter, +so the same code runs on NumPy and CuPy; the parameters are the same +(amplitude, decay, angular frequency, phase per mode, then the constant) +and, on convergence, the result matches the per-trace SciPy solver. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from time import perf_counter +from typing import Any + +import numpy as np + + +def _sync(xp) -> None: + dev = getattr(getattr(xp, "cuda", None), "Device", None) + if dev is not None: + dev().synchronize() + + +def model_and_jacobian_batched(xp, theta, time, n_modes: int): + """theta (B, 4K+1), time (N,) -> model (B, N), jacobian (B, N, 4K+1).""" + b = theta.shape[0] + n = time.shape[0] + p = 4 * n_modes + 1 + t = time[None, :] + model = xp.broadcast_to(theta[:, 4 * n_modes][:, None], (b, n)).copy() + jac = xp.empty((b, n, p), dtype=theta.dtype) + jac[:, :, 4 * n_modes] = 1.0 + for k in range(n_modes): + a = theta[:, 4 * k][:, None] + d = theta[:, 4 * k + 1][:, None] + w = theta[:, 4 * k + 2][:, None] + ph = theta[:, 4 * k + 3][:, None] + env = xp.exp(-d * t) + arg = w * t + ph + cos = xp.cos(arg) + sin = xp.sin(arg) + model += a * env * cos + jac[:, :, 4 * k] = env * cos + jac[:, :, 4 * k + 1] = -a * t * env * cos + jac[:, :, 4 * k + 2] = -a * t * env * sin + jac[:, :, 4 * k + 3] = -a * env * sin + return model, jac + + +@dataclass(frozen=True) +class BatchedRefinement: + """Refined parameters for a batch: ``theta`` (B, 4K+1) in the same + layout as the input, per-trace residual rms, convergence flags, the + number of iterations used and the wall time including any device + synchronisation.""" + + theta: Any + residual_rms: Any + converged: Any + iterations: int + elapsed_s: float + + +def refine_modes_batched( + xp, + time, + traces, + theta0, + n_modes: int, + *, + max_iter: int = 60, + tol: float = 1e-9, + lambda0: float = 1e-3, +) -> BatchedRefinement: + """Levenberg-Marquardt over a batch of traces (B, N) from starts + ``theta0`` (B, 4K+1). Damping is per trace: a step that lowers the + cost is accepted and its damping divided by 3, otherwise the damping + is multiplied by 5 and the trace keeps its parameters. Iteration + stops when every trace has converged or ``max_iter`` is reached.""" + t = xp.asarray(time, dtype=xp.float64) + y = xp.asarray(traces, dtype=xp.float64) + theta = xp.asarray(theta0, dtype=xp.float64).copy() + if y.ndim != 2 or theta.ndim != 2 or theta.shape[0] != y.shape[0]: + raise ValueError("traces must be (B, N) and theta0 (B, 4K+1)") + if theta.shape[1] != 4 * n_modes + 1: + raise ValueError("theta0 must have 4 * n_modes + 1 columns") + _sync(xp) + start = perf_counter() + b = y.shape[0] + lam = xp.full(b, lambda0, dtype=xp.float64) + model, jac = model_and_jacobian_batched(xp, theta, t, n_modes) + resid = model - y + cost = xp.sum(resid * resid, axis=1) + converged = xp.zeros(b, dtype=bool) + eye = xp.eye(theta.shape[1], dtype=xp.float64)[None, :, :] + iterations = 0 + for iterations in range(1, max_iter + 1): + jtj = xp.einsum("bnp,bnq->bpq", jac, jac) + jtr = xp.einsum("bnp,bn->bp", jac, resid) + diag = jtj * eye + step = -xp.linalg.solve( + jtj + lam[:, None, None] * diag, jtr[:, :, None] + )[:, :, 0] + trial = theta + step + model_t, jac_t = model_and_jacobian_batched(xp, trial, t, n_modes) + resid_t = model_t - y + cost_t = xp.sum(resid_t * resid_t, axis=1) + accept = (cost_t < cost) & ~converged + scale = xp.sqrt(xp.sum(theta * theta, axis=1)) + 1e-12 + small = xp.sqrt(xp.sum(step * step, axis=1)) < tol * scale + theta = xp.where(accept[:, None], trial, theta) + resid = xp.where(accept[:, None], resid_t, resid) + jac = xp.where(accept[:, None, None], jac_t, jac) + cost = xp.where(accept, cost_t, cost) + lam = xp.where(accept, lam / 3.0, lam * 5.0) + converged = converged | (accept & small) + if bool(xp.all(converged)): + break + _sync(xp) + rms = xp.sqrt(cost / y.shape[1]) + return BatchedRefinement( + theta=theta, + residual_rms=rms, + converged=converged, + iterations=iterations, + elapsed_s=perf_counter() - start, + ) + + +@dataclass(frozen=True) +class RefinementBenchmark: + traces: int + samples: int + n_modes: int + repeat: int + cpu_serial_scipy_s: float + numpy_batched_s: float + cupy_batched_s: float | None + max_abs_theta_diff_numpy: float + max_abs_theta_diff_cupy: float | None + gpu_error: str | None + + +def benchmark_refinement( + *, + samples: int = 96, + traces: int = 256, + n_modes: int = 2, + noise_sigma: float = 0.02, + seed: int = 0, + run_gpu: bool = True, + repeat: int = 3, +) -> RefinementBenchmark: + """Time the per-trace SciPy refinement against the batched solver on + NumPy and on CuPy, from identical perturbed starts, and report the + largest parameter difference between the batched and SciPy results. + + Every path is timed ``repeat`` times and the best time is reported. + The batched timings cover the solver only: the arrays are on the + device before the clock starts and the device is synchronised before + and after. The SciPy timing covers the Python loop over traces.""" + if repeat < 1: + raise ValueError("repeat must be positive") + from .mode_refinement import refine_modes + from .synthetic_validation import ( + DEFAULT_MODES, + modes_to_theta, + synthetic_modes_trace, + ) + + modes = DEFAULT_MODES[:n_modes] + fx = synthetic_modes_trace(samples, modes=modes) + rng = np.random.default_rng(seed) + y = fx.clean[None, :] + rng.normal(0.0, noise_sigma, (traces, samples)) + truth = modes_to_theta(modes, fx.constant) + theta0 = truth[None, :] * (1.0 + 0.05 * rng.standard_normal((traces, 1))) + theta0 = theta0 + 0.02 * rng.standard_normal(theta0.shape) + + # warm the device path once so compile and allocation costs are not + # charged to the timed run + if run_gpu: + try: + import cupy as cp + + refine_modes_batched(cp, fx.time, y[:16], theta0[:16], n_modes) + except Exception: # pragma: no cover - reported below + pass + + def serial_scipy(): + out = np.empty_like(theta0) + for i in range(traces): + init = [ + tuple(theta0[i, 4 * k : 4 * k + 4]) for k in range(n_modes) + ] + r = refine_modes(fx.time, y[i], init, float(theta0[i, -1])) + # refine_modes normalises signs and orders by amplitude; undo + # the ordering by matching frequencies to the start + for k in range(n_modes): + j = int( + np.argmin( + np.abs(r.angular_frequency - theta0[i, 4 * k + 2]) + ) + ) + out[i, 4 * k : 4 * k + 4] = ( + r.amplitude[j], + r.decay[j], + r.angular_frequency[j], + r.phase[j], + ) + out[i, -1] = r.constant + return out + + cpu_serial_s = float("inf") + for _ in range(repeat): + start = perf_counter() + serial = serial_scipy() + cpu_serial_s = min(cpu_serial_s, perf_counter() - start) + + res_np = min( + ( + refine_modes_batched(np, fx.time, y, theta0, n_modes) + for _ in range(repeat) + ), + key=lambda r: r.elapsed_s, + ) + diff_np = float( + np.max(np.abs(_canonical(res_np.theta, n_modes) - serial)) + ) + + cupy_s = None + diff_cp = None + gpu_error = None + if run_gpu: + try: + import cupy as cp + + res_cp = min( + ( + refine_modes_batched(cp, fx.time, y, theta0, n_modes) + for _ in range(repeat) + ), + key=lambda r: r.elapsed_s, + ) + theta_cp = cp.asnumpy(res_cp.theta) + cupy_s = res_cp.elapsed_s + diff_cp = float( + np.max(np.abs(_canonical(theta_cp, n_modes) - serial)) + ) + except Exception as exc: # pragma: no cover - depends on hardware + gpu_error = f"{type(exc).__name__}: {exc}" + return RefinementBenchmark( + traces=traces, + samples=samples, + n_modes=n_modes, + repeat=repeat, + cpu_serial_scipy_s=cpu_serial_s, + numpy_batched_s=res_np.elapsed_s, + cupy_batched_s=cupy_s, + max_abs_theta_diff_numpy=diff_np, + max_abs_theta_diff_cupy=diff_cp, + gpu_error=gpu_error, + ) + + +def _canonical(theta: np.ndarray, n_modes: int) -> np.ndarray: + """Same sign and phase conventions as :func:`refine_modes`.""" + out = np.array(theta, dtype=np.float64, copy=True) + for k in range(n_modes): + a = out[:, 4 * k] + w = out[:, 4 * k + 2] + ph = out[:, 4 * k + 3] + neg = a < 0 + a[neg] = -a[neg] + ph[neg] += np.pi + negf = w < 0 + w[negf] = -w[negf] + ph[negf] = -ph[negf] + out[:, 4 * k + 3] = (ph + np.pi) % (2.0 * np.pi) - np.pi + return out diff --git a/tests/core/test_cli_contract.py b/tests/core/test_cli_contract.py index 1b1b568f..8313b2e6 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, 144, 1), ("xscan", 42, 42, 171, 1), ("xrep", 6, 6, 103, 1), - ("xray", 34, 32, 398, 1), + ("xray", 35, 33, 404, 1), ] assert len(per_group) == 6 - assert sum(item[1] for item in per_group) == 93 - assert sum(item[1] + item[4] for item in per_group) == 98 - assert sum(item[2] for item in per_group) == 90 - assert sum(item[3] for item in per_group) == 847 + 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) == 853 def test_public_registry_order_and_component_derivations() -> None: diff --git a/tests/xray/test_cli.py b/tests/xray/test_cli.py index eec8d033..31af8a4b 100644 --- a/tests/xray/test_cli.py +++ b/tests/xray/test_cli.py @@ -65,6 +65,7 @@ def test_core_command_registry_loads_xray_commands(): "gpu-policy": "gp", "linear-prediction-benchmark": "lpb", "linear-prediction-validate": "lpv", + "linear-prediction-refine-benchmark": "lprb", "linear-prediction-fixed-stages-benchmark": "lpfsb", "linear-prediction-p2-benchmark": "lppb", "linear-prediction-profile-summary": "lpps", @@ -301,6 +302,27 @@ def broken(t, y, k): assert "any_mode_lost_rate=1.000" in output +def test_linear_prediction_validate_refine_adds_a_second_estimator(capsys): + assert ( + main(["lpv", "--trials", "3", "--snr-db", "30", "--refine", "--json"]) + == 0 + ) + payload = json.loads(capsys.readouterr().out) + assert payload["config"]["estimators"] == ["cpu", "cpu-refined"] + assert [r["estimator"] for r in payload["results"]] == [ + "cpu", + "cpu-refined", + ] + + +def test_linear_prediction_refine_benchmark_cli_cpu_only(capsys): + assert main(["lprb", "--traces", "4", "--no-gpu", "--json"]) == 0 + payload = json.loads(capsys.readouterr().out) + assert payload["traces"] == 4 + assert payload["cupy_batched_s"] is None + assert payload["max_abs_theta_diff_numpy"] < 1e-6 + + def test_gpu_policy(capsys): assert main(["gpu-policy"]) == 0 captured = capsys.readouterr() diff --git a/tests/xray/test_mode_refinement.py b/tests/xray/test_mode_refinement.py new file mode 100644 index 00000000..716bbc7d --- /dev/null +++ b/tests/xray/test_mode_refinement.py @@ -0,0 +1,99 @@ +# 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.mode_refinement import ( + _model_and_jacobian, + linear_prediction_refined, + refine_modes, + seed_from_residual, +) +from cuphoton.xray.synthetic_validation import ( + cramer_rao_bounds, + synthetic_modes_trace, + validation_sweep, +) + + +def test_analytic_jacobian_matches_finite_differences(): + fx = synthetic_modes_trace(96) + theta = np.array([1.25, 0.09, 2.4, 0.3, 0.45, 0.03, 0.9, -0.75, 0.15]) + _, jac = _model_and_jacobian(theta, fx.time, 2) + h = 1e-6 + for k in range(theta.size): + e = np.zeros_like(theta) + e[k] = h + num = ( + _model_and_jacobian(theta + e, fx.time, 2)[0] + - _model_and_jacobian(theta - e, fx.time, 2)[0] + ) / (2 * h) + np.testing.assert_allclose(jac[:, k], num, atol=1e-7) + + +def test_refinement_recovers_exact_modes_from_a_perturbed_start(): + fx = synthetic_modes_trace(96) + start = [(1.1, 0.05, 2.5, 0.0), (0.5, 0.0, 0.8, 0.0)] + r = refine_modes(fx.time, fx.trace, start, 0.0) + assert r.converged + np.testing.assert_allclose(r.angular_frequency, [2.4, 0.9], atol=1e-8) + np.testing.assert_allclose(r.decay, [0.09, 0.03], atol=1e-8) + np.testing.assert_allclose(r.amplitude, [1.25, 0.45], atol=1e-8) + assert abs(r.constant - 0.15) < 1e-8 + assert r.residual_rms < 1e-10 + + +def test_seed_from_residual_finds_the_missing_mode(): + fx = synthetic_modes_trace(256, duration=25.5) + weak = 0.45 * np.exp(-0.03 * fx.time) * np.cos(0.9 * fx.time - 0.75) + amp, dec, freq, _ = seed_from_residual(fx.time, weak) + assert abs(freq - 0.9) < 0.05 + assert dec == 0.0 and amp > 0 + + +def test_refined_estimator_keeps_the_weak_mode_under_noise(): + fx = synthetic_modes_trace(96) + sigma = float(np.std(fx.clean - fx.clean.mean())) / 100 + rng = np.random.default_rng(0) + lost = 0 + for _ in range(60): + y = fx.clean + rng.normal(0.0, sigma, size=fx.time.shape) + r = linear_prediction_refined(fx.time, y, 6, 2) + if np.min(np.abs(r.angular_frequency - 0.9)) > 0.05: + lost += 1 + assert lost == 0 + + +def test_refined_sweep_sits_near_the_cramer_rao_bound(): + sweep = validation_sweep( + snr_db=(30.0,), + trials=60, + seed=2, + estimator=lambda t, y, k: linear_prediction_refined(t, y, k, 2), + backend="cpu-refined", + ) + level = sweep.levels[0] + for mode in level.modes: + assert mode.loss_rate == 0.0 + assert 0.7 < mode.frequency_ratio < 1.4 + assert 0.7 < mode.decay_ratio < 1.4 + + +def test_refinement_rejects_empty_initial_modes(): + fx = synthetic_modes_trace(96) + with pytest.raises(ValueError): + refine_modes(fx.time, fx.trace, [], 0.0) + + +def test_bounds_are_the_reference_for_the_refined_uncertainties(): + # the per-fit one-sigma from the Jacobian should agree with the bound + # to within the scatter of a single residual-variance estimate + fx = synthetic_modes_trace(96) + sigma = float(np.std(fx.clean - fx.clean.mean())) / 1000 + b = cramer_rao_bounds(fx.modes, fx.constant, fx.time, sigma) + y = fx.clean + np.random.default_rng(5).normal(0.0, sigma, fx.time.shape) + r = linear_prediction_refined(fx.time, y, 6, 2) + ratio = r.sigma_angular_frequency / b["angular_frequency"] + assert np.all(ratio > 0.6) and np.all(ratio < 1.6) diff --git a/tests/xray/test_mode_refinement_batched.py b/tests/xray/test_mode_refinement_batched.py new file mode 100644 index 00000000..261dbb77 --- /dev/null +++ b/tests/xray/test_mode_refinement_batched.py @@ -0,0 +1,81 @@ +# 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.mode_refinement import _model_and_jacobian, refine_modes +from cuphoton.xray.mode_refinement_batched import ( + _canonical, + benchmark_refinement, + model_and_jacobian_batched, + refine_modes_batched, +) +from cuphoton.xray.synthetic_validation import ( + modes_to_theta, + synthetic_modes_trace, +) + + +def test_batched_model_and_jacobian_match_the_serial_ones(): + fx = synthetic_modes_trace(96) + theta = modes_to_theta(fx.modes, fx.constant) + batch = np.stack([theta, theta * 1.01, theta * 0.98]) + model, jac = model_and_jacobian_batched(np, batch, fx.time, 2) + for i in range(3): + m, j = _model_and_jacobian(batch[i], fx.time, 2) + np.testing.assert_allclose(model[i], m, rtol=0, atol=1e-12) + np.testing.assert_allclose(jac[i], j, rtol=0, atol=1e-12) + + +def test_batched_refinement_matches_scipy_per_trace(): + fx = synthetic_modes_trace(96) + rng = np.random.default_rng(1) + traces = 12 + y = fx.clean[None, :] + rng.normal(0.0, 0.02, (traces, fx.time.size)) + truth = modes_to_theta(fx.modes, fx.constant) + theta0 = truth[None, :] * (1.0 + 0.05 * rng.standard_normal((traces, 1))) + res = refine_modes_batched(np, fx.time, y, theta0, 2) + assert bool(np.all(res.converged)) + got = _canonical(res.theta, 2) + for i in range(traces): + init = [tuple(theta0[i, 4 * k : 4 * k + 4]) for k in range(2)] + r = refine_modes(fx.time, y[i], init, float(theta0[i, -1])) + j = int(np.argmin(np.abs(r.angular_frequency - 2.4))) + assert abs(got[i, 2] - r.angular_frequency[j]) < 1e-7 + assert abs(got[i, 1] - r.decay[j]) < 1e-7 + assert abs(got[i, -1] - r.constant) < 1e-7 + + +def test_batched_refinement_rejects_bad_shapes(): + fx = synthetic_modes_trace(96) + with pytest.raises(ValueError): + refine_modes_batched(np, fx.time, fx.trace, np.zeros((1, 9)), 2) + with pytest.raises(ValueError): + refine_modes_batched( + np, fx.time, fx.trace[None, :], np.zeros((1, 7)), 2 + ) + + +def test_benchmark_runs_on_cpu_and_reports_agreement(): + r = benchmark_refinement(traces=8, run_gpu=False) + assert r.traces == 8 + assert r.max_abs_theta_diff_numpy < 1e-6 + assert r.cupy_batched_s is None and r.gpu_error is None + + +def test_batched_refinement_on_cupy_matches_numpy(): + cp = pytest.importorskip("cupy") + try: + cp.cuda.runtime.getDeviceCount() + except Exception: # pragma: no cover - no device + pytest.skip("no CUDA device") + fx = synthetic_modes_trace(96) + rng = np.random.default_rng(2) + y = fx.clean[None, :] + rng.normal(0.0, 0.02, (32, fx.time.size)) + truth = modes_to_theta(fx.modes, fx.constant) + theta0 = truth[None, :] * (1.0 + 0.05 * rng.standard_normal((32, 1))) + a = refine_modes_batched(np, fx.time, y, theta0, 2) + b = refine_modes_batched(cp, fx.time, y, theta0, 2) + np.testing.assert_allclose(cp.asnumpy(b.theta), a.theta, atol=1e-7) From d21b1c1371137e2cd957140cc30e1aa46027cf30 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Fri, 18 Sep 2026 20:34:37 -0700 Subject: [PATCH 7/9] Correct refinement time origin and batched stopping Keep LP modal parameters relative to the first sample during refinement. Recognize stationary starts and use positive diagonal damping so a zero-amplitude mode does not make the batch solve singular. Signed-off-by: Trent Nelson --- src/cuphoton/xray/mode_refinement.py | 18 +++++--- src/cuphoton/xray/mode_refinement_batched.py | 18 +++++++- tests/xray/test_mode_refinement.py | 20 +++++++++ tests/xray/test_mode_refinement_batched.py | 45 ++++++++++++++++++++ 4 files changed, 93 insertions(+), 8 deletions(-) diff --git a/src/cuphoton/xray/mode_refinement.py b/src/cuphoton/xray/mode_refinement.py index be642730..549180ab 100644 --- a/src/cuphoton/xray/mode_refinement.py +++ b/src/cuphoton/xray/mode_refinement.py @@ -20,7 +20,7 @@ from __future__ import annotations -from dataclasses import dataclass +from dataclasses import dataclass, replace import numpy as np from scipy.optimize import least_squares @@ -180,12 +180,17 @@ def linear_prediction_refined( ) -> RefinedModes: """Linear prediction, then nonlinear refinement of the ``n_modes`` strongest oscillating modes; modes the linear-prediction step did not - return are seeded one at a time from the residual spectrum.""" + return are seeded one at a time from the residual spectrum. + + Amplitudes and phases use the first input sample as the time origin, + matching linear prediction; the returned time retains the input axis. + """ time = np.asarray(time, dtype=np.float64) trace = np.asarray(trace, dtype=np.float64) lp = linear_prediction_numpy( time, trace, n_components, roots_backend=roots_backend ) + fit_time = time - time[0] w = np.asarray(lp.angular_frequency, dtype=float) d = np.asarray(lp.decay, dtype=float) a = np.asarray(lp.amplitude, dtype=float) @@ -202,9 +207,9 @@ def linear_prediction_refined( residual = trace - np.asarray(lp.reconstruction, dtype=float) seeded = 0 while len(modes) < n_modes: - modes.append(seed_from_residual(time, residual)) + modes.append(seed_from_residual(fit_time, residual)) seeded += 1 - partial = refine_modes(time, trace, modes, constant) + partial = refine_modes(fit_time, trace, modes, constant) residual = trace - partial.reconstruction modes = [ ( @@ -216,11 +221,12 @@ def linear_prediction_refined( for i in range(len(modes)) ] constant = partial.constant - return refine_modes( - time, + result = refine_modes( + fit_time, trace, modes, constant, initial_mode_count=found, seeded_mode_count=seeded, ) + return replace(result, time=time) diff --git a/src/cuphoton/xray/mode_refinement_batched.py b/src/cuphoton/xray/mode_refinement_batched.py index 2dbc1d6b..949d7385 100644 --- a/src/cuphoton/xray/mode_refinement_batched.py +++ b/src/cuphoton/xray/mode_refinement_batched.py @@ -84,7 +84,11 @@ def refine_modes_batched( ``theta0`` (B, 4K+1). Damping is per trace: a step that lowers the cost is accepted and its damping divided by 3, otherwise the damping is multiplied by 5 and the trace keeps its parameters. Iteration - stops when every trace has converged or ``max_iter`` is reached.""" + stops when every trace has converged or ``max_iter`` is reached. + A scaled gradient check recognizes stationary starts independently of + step acceptance. Diagonal damping has a positive floor so an initially + zero-amplitude mode does not make the whole batch singular. + """ t = xp.asarray(time, dtype=xp.float64) y = xp.asarray(traces, dtype=xp.float64) theta = xp.asarray(theta0, dtype=xp.float64).copy() @@ -92,6 +96,8 @@ def refine_modes_batched( raise ValueError("traces must be (B, N) and theta0 (B, 4K+1)") if theta.shape[1] != 4 * n_modes + 1: raise ValueError("theta0 must have 4 * n_modes + 1 columns") + if not np.isfinite(lambda0) or lambda0 <= 0: + raise ValueError("lambda0 must be finite and positive") _sync(xp) start = perf_counter() b = y.shape[0] @@ -105,7 +111,15 @@ def refine_modes_batched( for iterations in range(1, max_iter + 1): jtj = xp.einsum("bnp,bnq->bpq", jac, jac) jtr = xp.einsum("bnp,bn->bp", jac, resid) - diag = jtj * eye + diagonal = xp.maximum(xp.diagonal(jtj, axis1=1, axis2=2), 1.0) + gradient_scale = ( + xp.sqrt(diagonal) * xp.maximum(xp.sqrt(cost), 1.0)[:, None] + ) + stationary = xp.max(xp.abs(jtr) / gradient_scale, axis=1) <= tol + converged = converged | stationary + if bool(xp.all(converged)): + break + diag = diagonal[:, :, None] * eye step = -xp.linalg.solve( jtj + lam[:, None, None] * diag, jtr[:, :, None] )[:, :, 0] diff --git a/tests/xray/test_mode_refinement.py b/tests/xray/test_mode_refinement.py index 716bbc7d..518b99f5 100644 --- a/tests/xray/test_mode_refinement.py +++ b/tests/xray/test_mode_refinement.py @@ -45,6 +45,26 @@ def test_refinement_recovers_exact_modes_from_a_perturbed_start(): assert r.residual_rms < 1e-10 +def test_lp_refinement_preserves_the_first_sample_time_origin(): + fx = synthetic_modes_trace(96) + base = linear_prediction_refined(fx.time, fx.trace, 6, 2) + shifted_time = fx.time + 1000.0 + shifted = linear_prediction_refined(shifted_time, fx.trace, 6, 2) + assert shifted.converged + assert shifted.residual_rms < 1e-10 + np.testing.assert_array_equal(shifted.time, shifted_time) + for field in ( + "amplitude", + "decay", + "angular_frequency", + "phase", + "reconstruction", + ): + np.testing.assert_allclose( + getattr(shifted, field), getattr(base, field), atol=1e-10 + ) + + def test_seed_from_residual_finds_the_missing_mode(): fx = synthetic_modes_trace(256, duration=25.5) weak = 0.45 * np.exp(-0.03 * fx.time) * np.cos(0.9 * fx.time - 0.75) diff --git a/tests/xray/test_mode_refinement_batched.py b/tests/xray/test_mode_refinement_batched.py index 261dbb77..30048fa3 100644 --- a/tests/xray/test_mode_refinement_batched.py +++ b/tests/xray/test_mode_refinement_batched.py @@ -58,6 +58,51 @@ def test_batched_refinement_rejects_bad_shapes(): ) +def test_batched_refinement_accepts_an_exact_start_without_improvement(): + fx = synthetic_modes_trace(96) + truth = modes_to_theta(fx.modes, fx.constant) + result = refine_modes_batched( + np, fx.time, fx.trace[None, :], truth[None, :], 2 + ) + assert result.converged.tolist() == [True] + assert result.iterations == 1 + assert result.residual_rms[0] < 1e-14 + + +def test_zero_amplitude_start_does_not_abort_other_rows(): + fx = synthetic_modes_trace(96) + truth = modes_to_theta(fx.modes, fx.constant) + starts = np.stack([truth, truth]) + starts[1, 4] = 0.0 + traces = np.stack([fx.trace, fx.trace]) + result = refine_modes_batched(np, fx.time, traces, starts, 2) + assert result.converged.tolist() == [True, True] + np.testing.assert_allclose( + result.theta, np.stack([truth, truth]), atol=1e-7 + ) + assert np.all(result.residual_rms < 1e-8) + + +def test_large_damping_does_not_make_a_nonstationary_start_converged(): + fx = synthetic_modes_trace(96) + start = modes_to_theta(fx.modes, fx.constant)[None, :] + start[0, -1] += 1.0 + result = refine_modes_batched( + np, fx.time, fx.trace[None, :], start, 2, lambda0=1e30, max_iter=1 + ) + assert result.converged.tolist() == [False] + + +@pytest.mark.parametrize("lambda0", [0.0, -1.0, np.nan, np.inf]) +def test_batched_refinement_requires_positive_finite_damping(lambda0): + fx = synthetic_modes_trace(96) + truth = modes_to_theta(fx.modes, fx.constant)[None, :] + with pytest.raises(ValueError, match="lambda0"): + refine_modes_batched( + np, fx.time, fx.trace[None, :], truth, 2, lambda0=lambda0 + ) + + def test_benchmark_runs_on_cpu_and_reports_agreement(): r = benchmark_refinement(traces=8, run_gpu=False) assert r.traces == 8 From 3c1049cf06a7acee6849844dd59ed296cd12ebea Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Thu, 24 Sep 2026 15:22:40 -0700 Subject: [PATCH 8/9] Make mode refinement convergence and uncertainty scale aware Normalize Jacobian columns for batched solves and require a small gradient to converge. Estimate local uncertainty with a rank-aware SVD and report unavailable uncertainties for rank-deficient fits. Signed-off-by: Trent Nelson --- src/cuphoton/xray/mode_refinement.py | 29 +++++++---- src/cuphoton/xray/mode_refinement_batched.py | 41 ++++++++------- tests/xray/test_mode_refinement.py | 45 ++++++++++++++++ tests/xray/test_mode_refinement_batched.py | 55 +++++++++++++++++--- 4 files changed, 134 insertions(+), 36 deletions(-) diff --git a/src/cuphoton/xray/mode_refinement.py b/src/cuphoton/xray/mode_refinement.py index 549180ab..632fd744 100644 --- a/src/cuphoton/xray/mode_refinement.py +++ b/src/cuphoton/xray/mode_refinement.py @@ -8,8 +8,8 @@ Linear prediction gives frequencies and decays without an initial guess, but it is not a maximum-likelihood estimator and its root filter can drop a lightly damped mode under noise. Refining the same damped-mode model by -nonlinear least squares from the linear-prediction result recovers the -Cramer-Rao bound, and a mode the filter dropped can be seeded from the +nonlinear least squares from the linear-prediction result can improve +parameter estimates, and a mode the filter dropped can be seeded from the peak of the residual spectrum before refinement. The model is trace(t) = c + sum_k A_k exp(-d_k t) cos(w_k t + p_k) @@ -33,9 +33,10 @@ class RefinedModes: """Refined damped modes and the fit diagnostics. Arrays are ordered by descending amplitude. ``sigma_*`` are one-sigma - uncertainties from the Jacobian at the solution and the residual - variance; ``initial_mode_count`` is how many oscillating modes the - linear-prediction step supplied and ``seeded_mode_count`` how many had + local uncertainties from the Jacobian and residual variance. They are + NaN when the scaled Jacobian is rank deficient or no residual degrees + of freedom remain. ``initial_mode_count`` is how many oscillating modes + the linear-prediction step supplied and ``seeded_mode_count`` how many had to be seeded from the residual spectrum. """ @@ -127,13 +128,21 @@ def jacobian(theta): theta = res.x model, jac = _model_and_jacobian(theta, time, n) resid = trace - model - dof = max(time.size - theta.size, 1) - var = float(resid @ resid) / dof + dof = time.size - theta.size + sigma = np.full(theta.size, np.nan) try: - cov = np.linalg.inv(jac.T @ jac) * var - sigma = np.sqrt(np.clip(np.diag(cov), 0.0, None)) + # Normalize columns before testing rank so parameter units do not + # determine identifiability. SVD avoids squaring the condition number. + scale = np.linalg.norm(jac, axis=0) + scale = np.where(scale > 0, scale, 1.0) + _, singular, vh = np.linalg.svd(jac / scale, full_matrices=False) + cutoff = np.finfo(float).eps * max(jac.shape) * singular[0] + if dof > 0 and np.all(singular > cutoff): + var = float(resid @ resid) / dof + sigma = np.sqrt(np.sum((vh / singular[:, None]) ** 2, axis=0)) + sigma *= np.sqrt(var) / scale except np.linalg.LinAlgError: - sigma = np.full(theta.size, np.nan) + pass amp = theta[0 : 4 * n : 4].copy() dec = theta[1 : 4 * n : 4].copy() freq = theta[2 : 4 * n : 4].copy() diff --git a/src/cuphoton/xray/mode_refinement_batched.py b/src/cuphoton/xray/mode_refinement_batched.py index 949d7385..2285bd42 100644 --- a/src/cuphoton/xray/mode_refinement_batched.py +++ b/src/cuphoton/xray/mode_refinement_batched.py @@ -10,8 +10,9 @@ broadcasting, normal equations solved as a batch of small (4K+1) systems, and per-trace damping with accept/reject. The array module is a parameter, so the same code runs on NumPy and CuPy; the parameters are the same -(amplitude, decay, angular frequency, phase per mode, then the constant) -and, on convergence, the result matches the per-trace SciPy solver. +(amplitude, decay, angular frequency, phase per mode, then the constant). +Both solvers find local optima; convergence alone does not establish that +they found the same optimum. """ from __future__ import annotations @@ -85,9 +86,11 @@ def refine_modes_batched( cost is accepted and its damping divided by 3, otherwise the damping is multiplied by 5 and the trace keeps its parameters. Iteration stops when every trace has converged or ``max_iter`` is reached. - A scaled gradient check recognizes stationary starts independently of - step acceptance. Diagonal damping has a positive floor so an initially - zero-amplitude mode does not make the whole batch singular. + Convergence requires a small gradient in Jacobian-scaled coordinates; + a step made tiny by damping does not imply stationarity. Scaling the + Jacobian columns makes the solve and stopping criterion independent of + parameter units. Positive damping also handles zero Jacobian columns + from an initial zero-amplitude mode. """ t = xp.asarray(time, dtype=xp.float64) y = xp.asarray(traces, dtype=xp.float64) @@ -109,35 +112,33 @@ def refine_modes_batched( eye = xp.eye(theta.shape[1], dtype=xp.float64)[None, :, :] iterations = 0 for iterations in range(1, max_iter + 1): - jtj = xp.einsum("bnp,bnq->bpq", jac, jac) - jtr = xp.einsum("bnp,bn->bp", jac, resid) - diagonal = xp.maximum(xp.diagonal(jtj, axis1=1, axis2=2), 1.0) - gradient_scale = ( - xp.sqrt(diagonal) * xp.maximum(xp.sqrt(cost), 1.0)[:, None] + scale = xp.sqrt(xp.sum(jac * jac, axis=1)) + scale = xp.where(scale > 0, scale, 1.0) + scaled_jac = jac / scale[:, None, :] + jtj = xp.einsum("bnp,bnq->bpq", scaled_jac, scaled_jac) + jtr = xp.einsum("bnp,bn->bp", scaled_jac, resid) + stationary = xp.max(xp.abs(jtr), axis=1) <= ( + tol * xp.maximum(xp.sqrt(cost), 1.0) ) - stationary = xp.max(xp.abs(jtr) / gradient_scale, axis=1) <= tol converged = converged | stationary if bool(xp.all(converged)): break - diag = diagonal[:, :, None] * eye - step = -xp.linalg.solve( - jtj + lam[:, None, None] * diag, jtr[:, :, None] - )[:, :, 0] + step = ( + -xp.linalg.solve(jtj + lam[:, None, None] * eye, jtr[:, :, None])[ + :, :, 0 + ] + / scale + ) trial = theta + step model_t, jac_t = model_and_jacobian_batched(xp, trial, t, n_modes) resid_t = model_t - y cost_t = xp.sum(resid_t * resid_t, axis=1) accept = (cost_t < cost) & ~converged - scale = xp.sqrt(xp.sum(theta * theta, axis=1)) + 1e-12 - small = xp.sqrt(xp.sum(step * step, axis=1)) < tol * scale theta = xp.where(accept[:, None], trial, theta) resid = xp.where(accept[:, None], resid_t, resid) jac = xp.where(accept[:, None, None], jac_t, jac) cost = xp.where(accept, cost_t, cost) lam = xp.where(accept, lam / 3.0, lam * 5.0) - converged = converged | (accept & small) - if bool(xp.all(converged)): - break _sync(xp) rms = xp.sqrt(cost / y.shape[1]) return BatchedRefinement( diff --git a/tests/xray/test_mode_refinement.py b/tests/xray/test_mode_refinement.py index 518b99f5..ca72de82 100644 --- a/tests/xray/test_mode_refinement.py +++ b/tests/xray/test_mode_refinement.py @@ -117,3 +117,48 @@ def test_bounds_are_the_reference_for_the_refined_uncertainties(): r = linear_prediction_refined(fx.time, y, 6, 2) ratio = r.sigma_angular_frequency / b["angular_frequency"] assert np.all(ratio > 0.6) and np.all(ratio < 1.6) + + +def test_rank_deficient_fit_reports_unavailable_uncertainty(): + fx = synthetic_modes_trace(96) + start = [(1.0, 0.09, 2.4, 0.3), (0.3, 0.09, 2.41, 0.3)] + theta = np.array([v for mode in start for v in mode] + [0.15]) + clean, _ = _model_and_jacobian(theta, fx.time, 2) + y = clean + np.random.default_rng(6).normal(0.0, 0.02, fx.time.size) + result = refine_modes(fx.time, y, start, 0.15) + assert result.converged + assert result.residual_rms < 0.025 + for field in ( + "sigma_amplitude", + "sigma_decay", + "sigma_angular_frequency", + "sigma_phase", + ): + assert np.all(np.isnan(getattr(result, field))) + + +@pytest.mark.parametrize("time_scale", [1e-9, 1e9]) +def test_identifiable_uncertainty_is_independent_of_time_units(time_scale): + fx = synthetic_modes_trace(96, noise_sigma=0.002, seed=5) + start = [ + (m.amplitude, m.decay, m.angular_frequency, m.phase) for m in fx.modes + ] + base = refine_modes(fx.time, fx.trace, start, fx.constant) + scaled_start = [ + (a, d / time_scale, w / time_scale, p) for a, d, w, p in start + ] + scaled = refine_modes( + fx.time * time_scale, fx.trace, scaled_start, fx.constant + ) + for field in ( + "sigma_amplitude", + "sigma_decay", + "sigma_angular_frequency", + "sigma_phase", + ): + actual = getattr(scaled, field) + if field in ("sigma_decay", "sigma_angular_frequency"): + actual = actual * time_scale + assert np.all(np.isfinite(actual)) + assert np.all(actual > 0) + np.testing.assert_allclose(actual, getattr(base, field), rtol=1e-4) diff --git a/tests/xray/test_mode_refinement_batched.py b/tests/xray/test_mode_refinement_batched.py index 30048fa3..8729f2f3 100644 --- a/tests/xray/test_mode_refinement_batched.py +++ b/tests/xray/test_mode_refinement_batched.py @@ -83,16 +83,46 @@ def test_zero_amplitude_start_does_not_abort_other_rows(): assert np.all(result.residual_rms < 1e-8) -def test_large_damping_does_not_make_a_nonstationary_start_converged(): +@pytest.mark.parametrize("lambda0", [1e10, 1e30]) +def test_large_damping_does_not_make_a_nonstationary_start_converged(lambda0): fx = synthetic_modes_trace(96) start = modes_to_theta(fx.modes, fx.constant)[None, :] start[0, -1] += 1.0 result = refine_modes_batched( - np, fx.time, fx.trace[None, :], start, 2, lambda0=1e30, max_iter=1 + np, fx.time, fx.trace[None, :], start, 2, lambda0=lambda0, max_iter=1 ) assert result.converged.tolist() == [False] +@pytest.mark.parametrize("time_scale", [1e-9, 1.0, 1e9]) +@pytest.mark.parametrize("backend", ["numpy", "cupy"]) +def test_batched_convergence_is_independent_of_time_units( + time_scale, backend +): + xp = np + if backend == "cupy": + xp = pytest.importorskip("cupy") + try: + xp.cuda.runtime.getDeviceCount() + except Exception: # pragma: no cover - no device + pytest.skip("no CUDA device") + fx = synthetic_modes_trace(96) + truth = modes_to_theta(fx.modes, fx.constant) + start = truth.copy() + start[[1, 2, 5, 6]] /= time_scale + start[[2, 6]] *= 1.05 + result = refine_modes_batched( + xp, fx.time * time_scale, fx.trace[None, :], start[None, :], 2 + ) + assert result.converged.tolist() == [True] + assert result.residual_rms[0] < 1e-8 + recovered = result.theta[0].copy() + if backend == "cupy": + recovered = xp.asnumpy(recovered) + recovered[[1, 2, 5, 6]] *= time_scale + np.testing.assert_allclose(recovered, truth, atol=1e-8) + + @pytest.mark.parametrize("lambda0", [0.0, -1.0, np.nan, np.inf]) def test_batched_refinement_requires_positive_finite_damping(lambda0): fx = synthetic_modes_trace(96) @@ -110,7 +140,8 @@ def test_benchmark_runs_on_cpu_and_reports_agreement(): assert r.cupy_batched_s is None and r.gpu_error is None -def test_batched_refinement_on_cupy_matches_numpy(): +@pytest.mark.parametrize("time_scale", [1e-9, 1.0, 1e9]) +def test_batched_refinement_on_cupy_matches_numpy(time_scale): cp = pytest.importorskip("cupy") try: cp.cuda.runtime.getDeviceCount() @@ -121,6 +152,18 @@ def test_batched_refinement_on_cupy_matches_numpy(): y = fx.clean[None, :] + rng.normal(0.0, 0.02, (32, fx.time.size)) truth = modes_to_theta(fx.modes, fx.constant) theta0 = truth[None, :] * (1.0 + 0.05 * rng.standard_normal((32, 1))) - a = refine_modes_batched(np, fx.time, y, theta0, 2) - b = refine_modes_batched(cp, fx.time, y, theta0, 2) - np.testing.assert_allclose(cp.asnumpy(b.theta), a.theta, atol=1e-7) + theta0[:, [1, 2, 5, 6]] /= time_scale + # Resolve noisy optima to a practical gradient tolerance; at 1e-9, + # rounding of cost differences can prevent further accepted steps. + a = refine_modes_batched(np, fx.time * time_scale, y, theta0, 2, tol=1e-8) + b = refine_modes_batched(cp, fx.time * time_scale, y, theta0, 2, tol=1e-8) + got = cp.asnumpy(b.theta) + expected = a.theta.copy() + got[:, [1, 2, 5, 6]] *= time_scale + expected[:, [1, 2, 5, 6]] *= time_scale + assert np.all(a.converged) + assert np.all(cp.asnumpy(b.converged)) + np.testing.assert_allclose(got, expected, rtol=0, atol=1e-7) + np.testing.assert_allclose( + cp.asnumpy(b.residual_rms), a.residual_rms, rtol=1e-10 + ) From 647fa994a0f0b3133bf64b5ff7265fc2d8371211 Mon Sep 17 00:00:00 2001 From: Trent Nelson Date: Thu, 24 Sep 2026 15:22:58 -0700 Subject: [PATCH 9/9] Document refinement methods without laptop performance claims Keep the reproducible harness and explain the experimental fit scope, local uncertainty assumptions, and convergence limits. Remove fixed-run accuracy statistics and laptop timing tables. Signed-off-by: Trent Nelson --- docs/xray/LINEAR-PREDICTION-REFINEMENT.md | 115 ++++++++++------------ docs/xray/README.md | 2 +- 2 files changed, 51 insertions(+), 66 deletions(-) diff --git a/docs/xray/LINEAR-PREDICTION-REFINEMENT.md b/docs/xray/LINEAR-PREDICTION-REFINEMENT.md index 7c1ff59a..9368ea26 100644 --- a/docs/xray/LINEAR-PREDICTION-REFINEMENT.md +++ b/docs/xray/LINEAR-PREDICTION-REFINEMENT.md @@ -1,18 +1,18 @@ # Nonlinear refinement of linear-prediction modes -Linear prediction gives frequencies and decays without an initial guess, but -it is not a maximum-likelihood estimator and its root filter can drop a -lightly damped mode under noise. `cuphoton.xray.mode_refinement` refines the -same damped-mode model by nonlinear least squares from the linear-prediction -result; `cuphoton.xray.mode_refinement_batched` does the same for a batch of -traces at once on NumPy or CuPy. Both fit +These experimental, opt-in tools refine linear-prediction modes by nonlinear +least squares. `cuphoton.xray.mode_refinement` fits one trace with SciPy; +`cuphoton.xray.mode_refinement_batched` fits a batch on NumPy or CuPy. Both use +an analytic Jacobian for the model ``` trace(t) = constant + sum_k amplitude_k exp(-decay_k t) cos(angular_frequency_k t + phase_k) ``` -with an analytic Jacobian. Nothing here changes the linear-prediction -commands or their defaults. +The fits are unconstrained: decay can become negative, and frequencies have +no bounds. They are separate from the existing production +`cuphoton.xray.iterative_fit` solver, which uses positive decay and bounded +frequencies. These tools have no detector or distributed pipeline integration. ## Per-trace refinement @@ -26,31 +26,26 @@ strongest peak of the residual spectrum (amplitude from the peak height, zero decay and phase) and refined before the next seed. `refine_modes` is the same solver from caller-supplied starting modes. -`RefinedModes` carries the refined parameters ordered by amplitude, one-sigma -uncertainties from the Jacobian at the solution and the residual variance, -the reconstruction, the residual rms, the solver status, and how many modes -came from linear prediction (`initial_mode_count`) versus the residual -spectrum (`seeded_mode_count`). A residual rms well above the noise means -the model class, not the fit, is wrong. +`RefinedModes` carries the refined parameters ordered by amplitude, local +one-sigma uncertainty estimates, the reconstruction, the residual rms, the +solver status, and how many modes came from linear prediction +(`initial_mode_count`) versus the residual spectrum (`seeded_mode_count`). +The uncertainty calculation uses an SVD of the column-normalized Jacobian +and the estimated residual variance. A rank-deficient Jacobian produces +`NaN` uncertainties. Finite estimates describe the local fit under the model +and independent, equal-variance noise assumptions; they do not establish +identifiability away from that solution or account for model mismatch. -`linear-prediction-validate --refine` runs the validation sweep a second time -with this estimator (`estimator: cpu-refined` in the summary) so the two can -be compared level by level. Reference run on the default fixture, CPU, -default seed, 300 trials per level: +A large residual can result from model mismatch, a poor local solution, or +failure to converge. Check the solver status and the reconstruction before +interpreting fitted parameters or uncertainties. -```bash -uv run cuphoton xray lpv --refine --trials 300 --snr-db 40,30,20,10,5 -``` - -In that run the linear-prediction estimator loses the weak mode in 10 to 52 -percent of trials depending on the level, and its frequency scatter is 1.1 -to 3.8 times the bound. The refined estimator loses the weak mode in 1 of -300 trials at each of the two lowest levels of that run and in none above, -and its frequency and decay scatter are within 10 percent of the Cramer-Rao -bound at every level of that run (ratios 0.91 to 1.07; the sampling -uncertainty of a standard deviation from 300 trials is about 4 percent). -The refinement does not remove the bias a model mismatch produces, and the -residual ratio still flags that case. +`linear-prediction-validate --refine` runs the validation sweep a second time +with this estimator (`estimator: cpu-refined` in the summary). Compare mode +loss, bias, scatter, and residuals across noise levels and distortions. +Results depend on the fixture, initialization, and noise realization; +refinement does not guarantee recovery of every mode or attainment of the +Cramer-Rao bound. ```bash uv run cuphoton xray lpv --refine --snr-db 40,20,10,5 --trials 200 --output-dir /tmp/lpv-refined @@ -68,18 +63,23 @@ an optional distortion. `refine_modes_batched(xp, time, traces, theta0, n_modes)` takes traces shaped `(B, N)` and starting parameters `(B, 4K+1)` in the same layout as the per-trace solver (amplitude, decay, angular frequency, phase per mode, then -the constant) and runs Levenberg-Marquardt for the whole batch: the Jacobian -for all traces by broadcasting, the normal equations as a batch of small -`(4K+1)` systems solved with `xp.linalg.solve`, and per-trace damping (a -step that lowers the cost is accepted and the damping divided by 3, otherwise -the damping is multiplied by 5 and the trace keeps its parameters). Iteration -stops when every trace has converged (step below `tol` times the parameter -norm, default 1e-9) or after `max_iter` (default 60). `xp` is `numpy` or -`cupy`. On convergence the result matches the per-trace SciPy solver to about -1e-8 in every parameter (table below). - -The batch iterates until its slowest trace converges, so the NumPy batched -path is slower than the SciPy loop; the batched form pays off on the GPU. +the constant). It normalizes the Jacobian columns and solves the damped +normal equations as a batch of small `(4K+1)` systems with `xp.linalg.solve`. +A step that lowers the cost is accepted and its damping divided by 3; +otherwise the damping is multiplied by 5 and that trace keeps its parameters. + +Convergence requires the largest absolute component of the +column-normalized Jacobian's residual gradient, divided by +`max(residual norm, 1)`, to be at most `tol` (default 1e-9). This criterion +is independent of the choice of time units. A tiny damped step alone does +not establish convergence. Iteration stops when every trace has converged +or after `max_iter` (default 60). Check each trace's convergence flag; +stationarity does not guarantee a good fit or agreement with another solver. + +At tight tolerances, rounding can prevent further cost reduction before the +gradient criterion is met. Such traces remain marked unconverged even when +their parameters agree closely with another fit. Choose a tolerance suited +to the required accuracy and inspect residuals as well as convergence flags. ## Benchmark @@ -91,8 +91,8 @@ NumPy, and the batched solver on CuPy. Each path is timed `--repeat` times (default 3) and the best is reported. For the batched paths the clock covers the solver only: the arrays are on the device before the clock starts and the device is synchronised before and after, so host-to-device transfer and -the copy of the result back are not included. One untimed CuPy call on 16 -traces runs first to absorb compilation and allocation. The command also +the copy of the result back are not included. One untimed CuPy call on up to +16 traces runs first to warm the device path. The command also reports the largest parameter difference between each batched result and the SciPy result. @@ -102,23 +102,8 @@ uv run cuphoton xray lprb --traces 2048 --repeat 3 --json uv run cuphoton xray lprb --traces 8192 --repeat 3 --json ``` -### NVIDIA GeForce RTX 4050 Laptop GPU (Ada, compute capability 8.9, 6 GiB) - -Host: Intel Core i9-13900H, WSL2 Ubuntu 24.04 on Windows 11, driver 596.49, -CUDA runtime 13.2, CuPy 14.1.1, Python 3.12, NumPy and SciPy from the locked -`dev` and `gpu` extras. Best of 3 repeats; the three commands above. - -| traces | SciPy per trace, serial (s) | NumPy batched (s) | CuPy batched (s) | SciPy / CuPy | max abs parameter difference | -| --- | --- | --- | --- | --- | --- | -| 256 | 0.135 | 0.258 | 0.108 | 1.3 | 8.7e-9 | -| 2048 | 1.09 | 2.49 | 0.162 | 6.7 | 1.2e-8 | -| 8192 | 4.38 | 12.7 | 0.800 | 5.5 | 1.9e-8 | - -On this card the GPU path breaks even at a few hundred traces and is 5 to 7 -times faster than the SciPy loop from about two thousand traces upward. -Timings at 256 traces vary by tens of percent between runs on this laptop -(a second run of the 256-trace command gave 0.178 s for CuPy); the larger -batches are stable to a few percent. In a -`cupyx.profiler.benchmark` breakdown at 256 traces on the same card, one -iteration costs about 1.5 ms in the Jacobian (launch-bound elementwise -work), 0.2 ms in the normal equations and 0.2 ms in the batched solve. +Use `--no-gpu` for the CPU comparison alone. These timings cover refinement +from supplied starts and exclude linear-prediction initialization, input I/O, +and detector processing. They cannot establish an end-to-end detector +speedup. Performance and parameter agreement need measurement on the target +hardware and representative traces. diff --git a/docs/xray/README.md b/docs/xray/README.md index 05316b13..39863c55 100644 --- a/docs/xray/README.md +++ b/docs/xray/README.md @@ -208,4 +208,4 @@ See [Data and artifact contracts](../data-artifacts.md#xray-hdf5-and-trace-produ - [Distributed detector artifacts](DISTRIBUTED-DETECTOR-ARTIFACTS.md) - [Validation visualization](VALIDATION-VIZ.md) - [Linear prediction validation](LINEAR-PREDICTION-VALIDATION.md) -- [Nonlinear refinement of linear-prediction modes](LINEAR-PREDICTION-REFINEMENT.md) +- [Experimental nonlinear refinement of linear-prediction modes](LINEAR-PREDICTION-REFINEMENT.md)