Files
csaxs_bec/bin/debug/scan_artefact_check.py
T
menzel a198d2aaa5
CI for csaxs_bec / test (push) Successful in 1m48s
debug: compare each end of a line against its neighbours, not the line mean
Four scans agreed that the first point of every line loses a fixed
number of milliseconds rather than a fixed fraction -- 5.18 ms at 15 ms
per point, 5.36 at 20, 5.93 and 6.04 at 50 -- while the fraction itself
moved from 35% to 12%. Scan 473 was the exception at 16.09 ms, and its
own last/bulk of 1.207 said why: the intensity rises 20% along every
line in that scan, so the whole-line mean is a biased baseline and
manufactures a deficit at the low end.

Both ends are now measured against their own k nearest neighbours, and
the along-line trend is reported next to them so the confounder is
visible rather than folded in. On a synthetic file carrying a true 6.0 ms
deficit under a 20% ramp, the old baseline reported 14.46 ms and the new
one reports 6.09.

A channel only contributes to the summary if its deficit exceeds three
standard errors over the lines. mca8 sits at two counts and was donating
a meaningless 0.71 ms to the median.

Positioners are also now excluded by having a setpoint rather than by
their trace being monotonic. The shape test discarded any intensity that
drifts smoothly along a line -- which is precisely the case that
motivated this commit.
2026-08-25 20:45:14 +02:00

1216 lines
54 KiB
Python
Executable File

#!/usr/bin/env python3
"""Desk check for the July 2026 commissioning scan artefacts.
Read-only. Opens BEC master files and answers two questions without needing beam:
Bug 1 -- is the last point of every line systematically weaker?
Quantifies it per monitored signal, and checks whether the *first*
point is affected too (which would point at a different mechanism).
Bug 2 -- did the fast axis actually move on every line?
Compares the per-line travel of the fast-axis readback. A line where
motion silently failed shows up as a collapsed span while the scan
still reported success.
Usage
-----
# what is in the file (run this first if anything looks odd)
python3 scan_artefact_check.py --list S00473_master.h5
# analyse one or more scans; a clean scan alongside a bad one is ideal
python3 scan_artefact_check.py S00473_master.h5 S00254_master.h5
# or by scan number, resolving the standard cSAXS layout
python3 scan_artefact_check.py --root /sls/x12sa/data/p23080/raw/data \
--scan 473 --scan 450 --scan 254
# what the scan was actually asked to do, straight from the file
python3 scan_artefact_check.py --root ... --scan 450 --command
# one line, point by point -- the only thing that names the end-of-line
# mechanism, because the whole-line span is identical for all of them
python3 scan_artefact_check.py --root ... --scan 450 --trace-line 20
Needs only h5py and numpy.
Caveat worth knowing before trusting Bug 2's answer: if the fast axis is an
open-loop stepper with no encoder, its "readback" is really the commanded
position, and lost steps are invisible here. In that case a collapsed span
proves motion failed, but a healthy span does *not* prove it succeeded.
"""
from __future__ import annotations
import argparse
import glob
import os
import sys
import h5py
import numpy as np
DEVICES_PATH = "/entry/collection/devices"
READOUT_PATH = "/entry/collection/readout_groups"
# --------------------------------------------------------------------------- #
# file handling
# --------------------------------------------------------------------------- #
def resolve_scan(root: str, number: int) -> str | None:
"""Find S<nnnnn>_master.h5 under the standard S00000-00999/S<nnnnn>/ layout."""
pattern = os.path.join(root, "S*", f"S{number:05d}", f"S{number:05d}_master.h5")
hits = glob.glob(pattern)
if not hits: # some deployments omit the bucket directory
hits = glob.glob(os.path.join(root, f"S{number:05d}", f"S{number:05d}_master.h5"))
return hits[0] if hits else None
def walk_datasets(h5: h5py.File, top: str = "/") -> list[tuple[str, tuple]]:
"""Every dataset in the file, as (path, shape). Follows soft links only."""
found: list[tuple[str, tuple]] = []
def visit(name, obj):
if isinstance(obj, h5py.Dataset):
found.append((name, obj.shape))
h5[top].visititems(visit)
return found
def numeric_signals(h5: h5py.File) -> dict[str, np.ndarray]:
"""Collect 1-D numeric per-point signals, keyed by 'device/signal'."""
out: dict[str, np.ndarray] = {}
if DEVICES_PATH not in h5:
return out
for dev_name, dev in h5[DEVICES_PATH].items():
if not isinstance(dev, h5py.Group):
continue
for sig_name, sig in dev.items():
if not isinstance(sig, h5py.Group) or "value" not in sig:
continue
data = sig["value"]
if not isinstance(data, h5py.Dataset) or data.ndim != 1:
continue
if not np.issubdtype(data.dtype, np.number):
continue
arr = np.asarray(data[()], dtype=float)
if arr.size > 1:
out[f"{dev_name}/{sig_name}"] = arr
return out
WANTED_META = ("scan_name", "scan_number", "num_points", "num_lines", "frames_per_trigger")
def read_metadata(h5: h5py.File) -> dict:
"""Find the scan bookkeeping wherever it lives.
BEC nests this differently across versions, so search the whole tree for the
keys by *basename* rather than assuming a path. Attributes count too.
"""
meta: dict = {}
def take(key: str, raw):
if key not in WANTED_META or key in meta:
return
try:
val = raw[()] if isinstance(raw, h5py.Dataset) else raw
if isinstance(val, np.ndarray) and val.size == 1:
val = val.reshape(-1)[0]
meta[key] = val.item() if hasattr(val, "item") else val
except Exception: # pylint: disable=broad-except
pass
def visit(name, obj):
take(name.rsplit("/", 1)[-1], obj)
for akey, aval in getattr(obj, "attrs", {}).items():
take(akey, aval)
h5.visititems(visit)
for akey, aval in h5["/"].attrs.items():
take(akey, aval)
return {k: (v.decode() if isinstance(v, bytes) else v) for k, v in meta.items()}
def positioner_devices(h5: h5py.File) -> set[str]:
"""Devices that are motors, identified by having a `<dev>_user_setpoint`.
Needed because judging "which signal is the scanned axis" by shape alone
picks nonsense: the SLS ring current decays monotonically through a line and
scores better than the motor does. A setpoint sibling is what actually makes
a device a positioner.
"""
out: set[str] = set()
if DEVICES_PATH not in h5:
return out
for dev_name, dev in h5[DEVICES_PATH].items():
if isinstance(dev, h5py.Group) and any("user_setpoint" in k for k in dev):
out.add(dev_name)
return out
def monitored_names(h5: h5py.File) -> set[str]:
grp = f"{READOUT_PATH}/monitored"
return set(h5[grp].keys()) if grp in h5 else set()
# --------------------------------------------------------------------------- #
# geometry
# --------------------------------------------------------------------------- #
def infer_grid(n: int, meta: dict) -> tuple[int, int] | None:
"""(num_lines, points_per_line) for n samples, from metadata or factoring."""
lines = meta.get("num_lines")
frames = meta.get("frames_per_trigger")
if lines and frames and int(lines) * int(frames) == n:
return int(lines), int(frames)
if frames and n % int(frames) == 0:
return n // int(frames), int(frames)
if lines and n % int(lines) == 0:
return int(lines), n // int(lines)
return None
def line_breaks(arr: np.ndarray) -> np.ndarray | None:
"""Indices where a scanned axis starts a new line.
Returns the boundaries themselves rather than a count, because the monitored
stream is sampled opportunistically: the number of samples per line varies,
so the trace cannot be reshaped into a rectangle. Segmenting on the actual
boundaries handles ragged lines, which is what these files contain.
"""
d = np.diff(arr)
nz = np.abs(d[d != 0])
if nz.size < 4:
return None
typical = float(np.median(nz))
if typical <= 0:
return None
for cand in (np.where(np.abs(d) > 5 * typical)[0] + 1, # raster flyback
np.where(np.diff(np.sign(np.where(d == 0, 1, d))) != 0)[0] + 1): # snake
if cand.size >= 2:
sizes = np.diff(np.concatenate(([0], cand, [arr.size])))
if sizes.min() >= 2 and sizes.size >= 3:
return cand
return None
def detect_lines(arr: np.ndarray) -> int | None:
"""Recover the number of lines from a scanned axis, without metadata.
A raster fast axis ramps and then snaps back, so the line boundaries are the
handful of steps that are far larger than a typical one. A snake axis never
snaps back but reverses, so the boundaries are where the step changes sign.
Try both and accept whichever gives evenly-sized segments.
Returns the number of lines, or None if the signal does not look scanned.
"""
d = np.diff(arr)
nz = np.abs(d[d != 0])
if nz.size < 4:
return None
typical = float(np.median(nz))
if typical <= 0:
return None
def segments_ok(breaks: np.ndarray) -> int | None:
"""Number of lines implied by these boundaries.
Uses the *median* spacing rather than requiring every segment to match.
A line where motion failed produces no flyback and no reversal, so its
boundary is missing and its segment comes out as a double -- and those
are precisely the lines this script exists to find. Demanding even
segments would therefore fail on exactly the files that matter.
"""
if breaks.size < 2:
return None
sizes = np.diff(np.concatenate(([0], breaks, [arr.size])))
sizes = sizes[sizes >= 3]
if sizes.size == 0:
return None
per_line = int(np.median(sizes))
if per_line < 3 or arr.size % per_line:
return None
return arr.size // per_line
# raster: flyback jumps
flyback = segments_ok(np.where(np.abs(d) > 5 * typical)[0] + 1)
if flyback:
return flyback
# snake: direction reversals
sign = np.sign(d)
sign[sign == 0] = 1
return segments_ok(np.where(np.diff(sign) != 0)[0] + 1)
def monotonicity(block: np.ndarray) -> float:
"""How steadily a signal ramps within a line: 1.0 = perfect ramp, 0 = noise.
This is what separates a scanned axis from an intensity. Judging by
within-line *range* alone is not enough -- a noisy detector channel can have
a larger spread than the axis it is plotted against, and then gets mistaken
for the fast axis (which is exactly what the first version of this script
did). Correlation against the point index does not care about scale.
"""
idx = np.arange(block.shape[1], dtype=float)
scores = []
for row in block:
if np.ptp(row) <= 0 or not np.all(np.isfinite(row)):
continue # a stalled line tells us nothing about what kind of signal this is
scores.append(abs(np.corrcoef(idx, row)[0, 1]))
return float(np.mean(scores)) if scores else 0.0
def pick_fast_axis(signals: dict[str, np.ndarray], grid: tuple[int, int]) -> str | None:
"""The fast axis is the signal that ramps most cleanly *within* each line."""
lines, per_line = grid
best, best_score = None, 0.0
for name, arr in signals.items():
if arr.size != lines * per_line:
continue
block = arr.reshape(lines, per_line)
if float(np.median(np.ptp(block, axis=1))) <= 0:
continue # constant within lines -> slow axis or a setting
score = monotonicity(block)
if score > best_score:
best, best_score = name, score
return best if best_score > 0.9 else None
# --------------------------------------------------------------------------- #
# the two checks
# --------------------------------------------------------------------------- #
def check_end_of_line(signals, grid, monitored, fast_axis, req=None,
motors=frozenset()) -> None:
lines, per_line = grid
print(f"\n BUG 1 -- end-of-line intensity ({lines} lines x {per_line} points)")
if per_line < 6:
print(" too few points per line to judge; skipping")
return
exp = (req or {}).get("exp_time")
exp = float(exp) if isinstance(exp, (int, float)) and exp > 0 else None
k = max(1, min(3, (per_line - 2) // 2))
head = f" {'signal':<30} {'first/near':>11} {'last/near':>10} {'trend':>7}"
if exp:
head += f" {'lost @start':>12}"
print(head + " verdict")
print(" " + "-" * (len(head) + 22))
deficits = []
for name, arr in sorted(signals.items()):
if arr.size != lines * per_line:
continue
if fast_axis and name == fast_axis:
continue
if name.split("/")[0] in motors:
continue # a positioner, identified by having a setpoint
block = np.asarray(arr, dtype=float).reshape(lines, per_line)
if not motors and monotonicity(block) > 0.9:
continue # no positioner list to go on, so fall back to the shape
# Each end is compared against its own neighbours, not the whole-line
# mean. A line with an intensity trend along it makes the whole-line
# mean a biased baseline and manufactures a deficit at the low end that
# has nothing to do with the first point -- scan 473 rises 20% along
# each line and reported three times the true figure because of it.
with np.errstate(divide="ignore", invalid="ignore"):
r_first = block[:, 0] / block[:, 1:1 + k].mean(axis=1)
r_last = block[:, -1] / block[:, -1 - k:-1].mean(axis=1)
r_first = r_first[np.isfinite(r_first)]
r_last = r_last[np.isfinite(r_last)]
if r_first.size < 3 or r_last.size < 3:
continue
first, last = float(np.mean(r_first)), float(np.mean(r_last))
# Scatter over lines, so a channel with too few counts to say anything
# cannot contribute a number. mca8 sits at 2 counts and was polluting
# the summary with a deficit indistinguishable from Poisson noise.
sem = float(np.std(r_first)) / np.sqrt(r_first.size)
interior = block[:, 1:-1]
base = float(np.mean(interior[:, :k]))
trend = float(np.mean(interior[:, -k:])) / base if base else float("nan")
real = (1.0 - first) > max(3.0 * sem, 0.02)
if abs(last - 1) < 0.02 and not real:
verdict = "flat"
elif real and last < 0.98:
verdict = "both ends low (different mechanism)"
elif real:
verdict = "FIRST POINT LOW"
elif last < 0.98:
verdict = "LAST POINT LOW"
else:
verdict = ""
row = f" {name:<30} {first:>11.3f} {last:>10.3f} {trend:>7.3f}"
if exp:
if real:
lost = (1.0 - first) * exp
deficits.append(lost)
row += f" {lost * 1e3:>10.2f}ms"
else:
row += f" {'-':>12}"
print(row + f" {verdict}")
print(f" first/near = first point / mean of the next {k}; last/near likewise.")
print(" trend = interior level at the end of a line / at its start; far from")
print(" 1.00 means the line itself drifts, which is a separate effect.")
if deficits:
vals = np.array(deficits)
print(f"\n lost exposure at the first point of every line:"
f" {np.median(vals) * 1e3:.2f} ms"
f" (over {len(deficits)} channel(s), this scan ran at {exp * 1e3:.0f} ms/point)")
print(" A shutter opening late against the first gate costs a fixed number of")
print(" milliseconds, so this figure must repeat across scans taken at other")
print(" exposure times. A constant fraction instead would mean it scales with")
print(" dose and the shutter is innocent.")
# End-of-*scan* behaviour is a separate question from end-of-*line*.
print(f"\n Trend across the scan (last tenth of lines vs first tenth)")
for name, arr in sorted(signals.items()):
if arr.size != lines * per_line:
continue
if name.split("/")[0] in motors:
continue
block = np.asarray(arr, dtype=float).reshape(lines, per_line)
if not motors and monotonicity(block) > 0.9:
continue
j = max(1, lines // 10)
first_tenth, last_tenth = np.mean(block[:j]), np.mean(block[-j:])
if not np.isfinite(first_tenth) or abs(first_tenth) < 1e-30:
continue
rel = last_tenth / first_tenth
flag = " <-- degrades towards the end" if rel < 0.95 else ""
if abs(rel - 1) > 0.02:
print(f" {name:<44} {rel:>7.3f}{flag}")
def check_stationary_from_data(signals, grid, motors, req) -> None:
"""Where the stage stopped moving, from the data alone.
These stages have no encoders, so nothing can be asked where the sample
actually was. The data answers it anyway: two points taken at the same
physical position differ only by noise, while two taken at different
positions differ by the sample as well. Comparing every point of a line
against the last point, across all lines at once, shows where that stops
being true -- and the line-to-line pattern is the sample's own fingerprint,
so this needs no knowledge of what the sample is.
The discriminator is the *trend*, not the level: while the stage moves, the
residual grows with distance from the reference point, because the sample
decorrelates. While it stands still, the residual is flat at the noise
floor. A smooth sample has small residuals either way, which is why a plain
threshold on the residual would not do.
"""
lines, per_line = grid
print("\n WHERE THE STAGE STOPPED MOVING -- from the data, no encoder needed")
if lines < 8 or per_line < 8:
print(" too few lines or points per line to compare columns; skipping")
return
total = lines * per_line
exp = (req or {}).get("exp_time")
exp = float(exp) if isinstance(exp, (int, float)) and exp > 0 else None
reported = 0
for name, arr in sorted(signals.items()):
if arr.size != total or name.split("/")[0] in motors:
continue
block = np.asarray(arr, dtype=float).reshape(lines, per_line)
cols = block.T.copy()
# Remove each point's own level: an intensity drift along the line is a
# different effect, and would otherwise swamp the pattern comparison.
cols -= cols.mean(axis=1, keepdims=True)
for end, ref_i, label in ((True, -1, "last"), (False, 0, "first")):
contrast = float(np.std(cols[ref_i]))
if not np.isfinite(contrast) or contrast <= 0:
continue
resid = np.array([float(np.std(c - cols[ref_i])) for c in cols]) / contrast
base = float(resid[-2] if end else resid[1])
if not 0 < base < 0.45:
continue # noise-dominated: the sample leaves no fingerprint here
order = range(per_line - 2, -1, -1) if end else range(1, per_line)
run = 0
for i in order:
if resid[i] > 2 * base:
break
run += 1
if run < 3:
continue
if reported == 0:
print(" residual against the {} point of each line, in units of the"
" sample's".format(label))
print(" own line-to-line contrast. Flat = not moving; rising = moving.")
reported += 1
print(f"\n {name} contrast {contrast:.4g}, neighbour residual {base:.3f}")
k = min(per_line - 1, 18)
idx = np.unique(np.linspace(0, per_line - 2, k).astype(int))
print(" point " + " ".join(f"{i:>5d}" for i in idx))
print(" resid " + " ".join(f"{resid[i]:>5.2f}" for i in idx))
moving = per_line - run
print(f" the {label} {run} of {per_line} points are indistinguishable from"
f" the {label} point.")
print(f" The stage covered the line in ~{moving} points, so it ran about"
f" {(per_line - 1) / max(1, moving - 1):.2f}x")
print(" faster than the burst assumed"
+ (f", and {run * exp:.3g} s of the {per_line * exp:.3g} s burst"
" fired into a stage that was not moving." if exp else "."))
if not reported:
print(" no signal shows a run of points indistinguishable from the line end.")
print(" Either the stage moved throughout, or the sample has too little")
print(" line-to-line structure for this comparison to see anything.")
def report_distribution(flagged: np.ndarray, total: int, label: str) -> None:
"""Where in the scan the flagged lines fall.
Printing only the first N hides a cluster at the end -- and the reported
artefacts are said to appear towards the end of scans, which is exactly the
case a head-truncated list cannot show. So summarise position, then print
both ends of the list rather than just the top.
"""
if flagged.size == 0 or total <= 0:
return
tenths = np.zeros(10, dtype=int)
for i in flagged:
tenths[min(9, int(10 * i / total))] += 1
bar = " ".join(f"{c:>3}" for c in tenths)
print(f" {label} by tenth of scan (start -> end):")
print(f" {bar}")
last_tenth = tenths[-1] / max(1, flagged.size)
first_tenth = tenths[0] / max(1, flagged.size)
if last_tenth > 0.3:
print(f" -> {last_tenth:.0%} of them are in the final tenth: this clusters"
" at the END of the scan.")
elif first_tenth > 0.3:
print(f" -> {first_tenth:.0%} are in the first tenth: clusters at the START.")
else:
print(" -> spread through the scan, no obvious end-of-scan clustering.")
def numeric_timestamps(h5: h5py.File) -> dict[str, np.ndarray]:
"""Timestamps for the same keys numeric_signals returns, where present.
The monitored stream and the scan grid run on different clocks, so a
position and an intensity carrying the same index are not the same instant.
Timestamps are also what turns three position samples per line into a
velocity, which is the number that decides whether the axis ran too fast.
"""
out: dict[str, np.ndarray] = {}
if DEVICES_PATH not in h5:
return out
for dev_name, dev in h5[DEVICES_PATH].items():
if not isinstance(dev, h5py.Group):
continue
for sig_name, sig in dev.items():
if not isinstance(sig, h5py.Group) or "timestamp" not in sig:
continue
data = sig["timestamp"]
if not isinstance(data, h5py.Dataset) or data.ndim != 1:
continue
try:
arr = np.asarray(data[()], dtype=float)
except Exception: # pylint: disable=broad-except
continue
if arr.size > 1:
out[f"{dev_name}/{sig_name}"] = arr
return out
def read_request_inputs(h5: h5py.File) -> dict:
"""The scan request as a flat dict -- what the scan was actually asked to do.
Same source as --command, but returned instead of printed, so the checks can
measure the file against the command rather than against a guess. Written by
the thing that ran the scan, so it beats notebooks and memory.
"""
out: dict = {}
root = "/entry/collection/metadata/bec"
if root not in h5:
return out
for sub in ("request_inputs/inputs", "request_inputs/kwargs", "scan_parameters"):
path = f"{root}/{sub}"
if path not in h5:
continue
for key, node in h5[path].items():
if not isinstance(node, h5py.Dataset):
continue
try:
val = node[()]
except Exception: # pylint: disable=broad-except
continue
if isinstance(val, np.ndarray):
if val.size != 1:
continue
val = val.reshape(-1)[0]
if hasattr(val, "item"):
val = val.item()
if isinstance(val, bytes):
val = val.decode(errors="replace")
out.setdefault(key, val)
return out
# -- shared line segmentation, used by both Bug 2 and the line trace --------- #
def motor_candidates(signals, motors) -> dict:
"""Positioner readbacks only -- setpoints would answer the wrong question."""
return {n: a for n, a in signals.items()
if n.split("/")[0] in motors and "setpoint" not in n}
def slow_axis_breaks(cands: dict, expected: int | None):
"""(name, break indices) for the positioner that steps once per line.
The fast axis cannot delimit its own lines: a line where motion failed has
no flyback, so it merges into its neighbour and the resulting segment looks
perfectly healthy -- the stall erases the evidence of itself. The slow axis
steps once per line regardless of what the fast axis did.
"""
for name, arr in sorted(cands.items()):
changes = np.where(np.abs(np.diff(arr)) > 0)[0] + 1
if changes.size == 0:
continue
if expected and abs(changes.size - (expected - 1)) <= max(1, 0.05 * expected):
return name, changes
return None, None
def rank_fast_traces(cands: dict, slow_name, slow_breaks) -> list:
"""[(name, arr, bounds, segs, spans)] ranked by travel per line, most first.
The scanned axis is the one that travels furthest within a line; the slow
axis steps a little and everything else sits still.
"""
scored = []
for name, arr in sorted(cands.items()):
if name == slow_name:
continue
breaks = slow_breaks if slow_breaks is not None else line_breaks(arr)
if breaks is None:
continue
edges = np.concatenate(([0], np.asarray(breaks, dtype=int), [arr.size]))
bounds = [(int(a), int(b)) for a, b in zip(edges[:-1], edges[1:]) if b > a]
segs = [arr[a:b] for a, b in bounds]
spans = np.array([float(np.ptp(sg)) for sg in segs])
if spans.size == 0 or np.median(spans) <= 0:
continue
scored.append((name, arr, bounds, segs, spans))
scored.sort(key=lambda t: -float(np.median(t[4])))
return scored
def readback_coverage(segs, travel_cmd: float, direction: float):
"""(coverage, forward), each as a fraction of the commanded travel.
coverage -- how much of the line the samples span at all.
forward -- how far they get in the direction the line was commanded to run.
Repetition across lines cannot be the test for "is this a trace of the
motion": a healthy continuous scan repeats too, because every line does the
same thing. What separates a trace from a handful of loop checkpoints is
whether the samples follow the line at all.
"""
if travel_cmd <= 0 or not direction:
return None, None
cov, fwd = [], []
for sg in segs:
if sg.size < 2:
continue
cov.append(float(np.ptp(sg)) / travel_cmd)
cov_adv = direction * (sg - sg[0])
fwd.append(float(cov_adv.max()) / travel_cmd)
if not cov:
return None, None
return float(np.median(cov)), float(np.median(fwd))
def check_line_motion_from_motor(signals, motors, meta, req=None) -> None:
"""Bug 2, using a real positioner trace rather than whatever looks rampy."""
print("\n BUG 2 -- did the fast axis move on every line?")
cands = motor_candidates(signals, motors)
if not cands:
print(" no positioner readback in this file, so this cannot be answered here.")
print(" (cont_grid records motor position only in the monitored stream;")
print(" if that is absent the file cannot show whether a line moved.)")
return
expected = int(meta["num_lines"]) if meta.get("num_lines") else None
slow_name, slow_breaks = slow_axis_breaks(cands, expected)
if slow_breaks is not None:
print(f" lines delimited by the slow axis {slow_name}"
f" ({slow_breaks.size + 1} segments)")
scored = rank_fast_traces(cands, slow_name, slow_breaks)
if not scored:
print(f" positioner traces found ({', '.join(sorted(cands))}) but none shows")
print(" a line structure. For a continuous scan the position is sampled")
print(" while the line runs, so an irregular or too-sparse trace can hide it.")
return
if len(scored) > 1:
others = ", ".join(f"{n} (median span {np.median(sp):.4g}, {len(sg)} segments)"
for n, _, _bd, sg, sp in scored[1:4])
print(f" other positioner traces: {others}")
name, arr, _bounds, segs, spans = scored[0]
found = len(segs)
per_line = int(np.median([len(sg) for sg in segs]))
print(f" fast axis : {name} ({arr.size} samples, {found} line segments"
+ (f", metadata says {expected} lines" if expected else "") + ")")
if expected and abs(found - expected) > max(1, 0.02 * expected):
print(f" NOTE: segmented into {found} but metadata says {expected}. The extra"
" or missing")
print(" boundaries are themselves suspicious -- a stalled line has no flyback.")
median = float(np.median(spans))
if median <= 0:
print(" the axis does not move -- check that this is the scanned motor")
return
ratio = spans / median
print(f" samples per line : ~{per_line} median line span : {median:.6g}")
fs, fe = (req or {}).get("fast_start"), (req or {}).get("fast_end")
if isinstance(fs, (int, float)) and isinstance(fe, (int, float)) and fe != fs:
cov, fwd = readback_coverage(segs, abs(fe - fs), np.sign(fe - fs))
if cov is not None:
print(f" readback covers {cov:.0%} of the commanded travel,"
f" {fwd:.0%} of it in the commanded direction")
if fwd < 0.1:
print("\n STOP: the readback never advances along the commanded")
print(" direction. It is sampled at fixed points in the scan loop --")
print(" after the burst and during the return -- and not while the")
print(" line runs. The spans, stalls and step sizes below are")
print(" properties of the loop, not of the stage, and none of them is")
print(" evidence about the motion.\n")
elif cov < 0.6:
print(f"\n CAUTION: the readback spans only {cov:.0%} of the line, so"
" it misses")
print(" part of every traverse. The travel below is understated.\n")
# Enough geometry to check this file against the scan that was intended.
# A step size that disagrees with the scan command means either the wrong
# axis has been identified or the motion was not what was asked for -- and
# both make everything below meaningless.
grid_pts = int(meta.get("frames_per_trigger") or 0)
print(f" fast axis range : {arr.min():.6g} .. {arr.max():.6g}")
if slow_name is not None:
sarr = signals[slow_name]
print(f" slow axis range : {sarr.min():.6g} .. {sarr.max():.6g}"
f" ({found} lines -> {(sarr.max() - sarr.min()) / max(1, found - 1):.4g} per line)")
if grid_pts:
print(f" implied step : {median / grid_pts:.4g} per point"
f" ({median:.4g} over {grid_pts} points)")
print(" ^ compare with the step size in the scan command. A mismatch means")
print(" the wrong axis, the wrong units, or motion that was not commanded.")
# Whole-line spans only catch a line that never moved. An axis that stalls
# part-way and then continues still spans the full range, so look for runs of
# samples that did not advance.
stalled = []
for i, sg in enumerate(segs):
if sg.size < 4:
continue
step = np.abs(np.diff(sg))
typical = float(np.median(step))
if typical <= 0:
continue
frozen = step < 0.1 * typical
run, longest = 0, 0
for f in frozen:
run = run + 1 if f else 0
longest = max(longest, run)
if longest >= 3:
stalled.append((i, longest, sg.size))
if stalled:
print(f"\n within-line stalls: {len(stalled)} lines contain a run of >=3"
" consecutive samples that did not advance")
report_distribution(np.array([i for i, _, _ in stalled]), found, "within-line stalls")
for i, longest, size in stalled[:10]:
print(f" line {i:>5} {longest} frozen of {size} samples")
print(" A line can span its full range and still have stopped part-way;")
print(" this is the partial-stall case the span check cannot see.")
else:
print(" no within-line stalls either (no run of >=3 frozen samples).")
if per_line < 5:
print(f" WARNING: only ~{per_line} position samples per line. The motor trace")
print(" is far sparser than the data, so a partial stall can hide between")
print(" samples. Absence of evidence here is weak evidence of absence.")
dead = np.where(ratio < 0.1)[0]
low = np.where(ratio < 0.9)[0]
print(f" lines below 90% of median span : {len(low)}")
print(f" lines below 10% (motion absent) : {len(dead)}")
report_distribution(low, found, "lines below 90%")
report_distribution(dead, found, "lines with motion absent")
for i in list(dead[:10]) + ([-1] if len(dead) > 20 else []) + list(dead[-10:] if len(dead) > 20 else []):
if i == -1:
print(" ...")
continue
print(f" line {i:>5} span {spans[i]:.6g} ({ratio[i]:.2f}x)")
if len(dead):
print(" Collapsed spans are motion that did not happen while the scan still")
print(" reported success -- the signature we are looking for.")
elif len(low):
print(" No line is missing outright; the spread is a velocity/settling")
print(" variation rather than dropped motion.")
else:
print(" every line travelled a consistent distance.")
def check_line_motion(signals, grid, fast_axis) -> None:
lines, per_line = grid
print("\n BUG 2 -- did the fast axis move on every line?")
if fast_axis is None:
print(" could not identify a fast axis; re-run with --list and tell me the names")
return
block = signals[fast_axis].reshape(lines, per_line)
spans = np.ptp(block, axis=1)
median = float(np.median(spans))
print(f" fast axis : {fast_axis}")
print(f" median line span : {median:.6g}")
if median <= 0:
print(" the fast axis does not move at all in this file -- check the signal choice")
return
ratio = spans / median
bad = np.where(ratio < 0.9)[0]
dead = np.where(ratio < 0.1)[0]
print(f" lines below 90% of median span : {len(bad)}")
print(f" lines below 10% (motion absent) : {len(dead)}")
if len(bad):
print(f" {'line':>6} {'span':>14} {'vs median':>10}")
for i in bad[:25]:
print(f" {i:>6} {spans[i]:>14.6g} {ratio[i]:>9.2f}x")
if len(bad) > 25:
print(f" ... and {len(bad) - 25} more")
print("\n Lines with a collapsed span are motion that did not happen while the")
print(" scan still reported success -- the signature we are looking for.")
else:
print(" all lines travelled a consistent distance; no dropped motion visible here.")
starts, ends = block[:, 0], block[:, -1]
print(f" line start scatter : {np.std(starts):.3g} line end scatter : {np.std(ends):.3g}")
print(" (large scatter at one end suggests the return move races the next line)")
# --------------------------------------------------------------------------- #
def trace_line(h5, signals, times, grid, meta, motors, index) -> None:
"""Point-by-point dump of a single line.
Whole-line spans cannot separate the end-of-line mechanisms: a stage that
runs too fast and parks, one that overshoots and settles back, and one whose
flyback starts before the last trigger all travel the same distance. They
differ only in the *shape* of the position trace inside the line:
ramp then flat -> arrived early and parked
ramp past the target -> overshoot, then the servo pulls it back
ramp then reverses -> the return move began before the last trigger
flat while the setpoint -> the readback is not the real position
ramps
Each of those implies a different fix, which is why this exists.
"""
lines, per_line = grid
print(f"\n LINE TRACE -- line {index} of {lines}")
if not 0 <= index < lines:
print(f" outside 0..{lines - 1}")
return
req = read_request_inputs(h5)
start, end = req.get("fast_start"), req.get("fast_end")
step, exp, axis = req.get("fast_step_size"), req.get("exp_time"), req.get("fast_axis")
nominal, v_cmd = None, None
if isinstance(start, (int, float)) and isinstance(end, (int, float)):
nominal = np.linspace(float(start), float(end), per_line)
print(f" commanded : {axis or 'fast axis'} {start:.6f} -> {end:.6f}"
f" ({abs(end - start):.6g} over {per_line} points"
f", {'descending' if end < start else 'ascending'})")
if isinstance(exp, (int, float)) and exp > 0:
burst = per_line * float(exp)
line = f" timing : {exp:g} s/point -> burst {burst:.4g} s"
if isinstance(step, (int, float)) and step > 0:
v_cmd = float(step) / float(exp)
line += f", commanded velocity {v_cmd:.6g}/s"
print(line)
print(" (the traverse is meant to fill the whole burst)")
_trace_grid_data(signals, motors, grid, nominal, index)
_trace_position(signals, times, motors, meta, nominal, v_cmd, index, per_line)
def _trace_grid_data(signals, motors, grid, nominal, index) -> None:
"""The per-point data for this line, next to the position it was meant for."""
lines, per_line = grid
total = lines * per_line
sl = slice(index * per_line, (index + 1) * per_line)
ranked = []
for name, arr in signals.items():
if arr.size != total or name.split("/")[0] in motors:
continue
row = np.asarray(arr[sl], dtype=float)
bulk = float(np.median(np.abs(row)))
if not np.isfinite(bulk) or bulk <= 0:
continue
# Range in units of counting noise, not of level. A channel sitting at
# 2 counts has an enormous *relative* range and no information in it;
# ranking by ptp/median floats that Poisson noise to the top and buries
# the diode that actually formed the image.
ranked.append((float(np.ptp(row)) / max(np.sqrt(bulk), 1.0), name, row))
if not ranked:
print("\n no per-point data on the scan grid for this line.")
return
# Most structure first: a channel that is flat within the line says nothing
# about where the stage was.
ranked.sort(key=lambda t: -t[0])
show = ranked[:3]
lead = show[0][2]
lo, hi = float(np.min(lead)), float(np.max(lead))
print(f"\n data on the scan grid ({per_line} points, leading channel {show[0][1]}):")
header = f" {'idx':>4} {'nominal' if nominal is not None else 'point':>13}"
for _, name, _ in show:
header += f" {name.split('/')[-1][:14]:>15}"
print(header + " profile")
# A long line would bury the two ends, which are the only parts in question.
if per_line > 120:
rows = list(range(60)) + [-1] + list(range(per_line - 60, per_line))
else:
rows = list(range(per_line))
for i in rows:
if i == -1:
print(f" {'...':>4}")
continue
pos = f"{nominal[i]:13.6f}" if nominal is not None else f"{i:>13}"
cells = "".join(f" {row[i]:>15.6g}" for _, _, row in show)
frac = 0.0 if hi <= lo else (lead[i] - lo) / (hi - lo)
print(f" {i:>4} {pos}{cells} {'#' * int(round(frac * 24))}")
interior = lead[1:-1]
if interior.size:
bulk = float(np.median(interior))
if bulk:
print(f" first/bulk {lead[0] / bulk:.3f} last/bulk {lead[-1] / bulk:.3f}"
" (a shortfall at one end only is a truncated gate, not the sample)")
def _trace_position(signals, times, motors, meta, nominal, v_cmd, index, per_line) -> None:
"""The position trace for this line -- the part that names the mechanism."""
cands = motor_candidates(signals, motors)
if not cands:
print("\n no positioner readback in this file, so the position trace --")
print(" the decisive part of this check -- cannot be recovered.")
return
expected = int(meta["num_lines"]) if meta.get("num_lines") else None
slow_name, slow_breaks = slow_axis_breaks(cands, expected)
scored = rank_fast_traces(cands, slow_name, slow_breaks)
if not scored:
print("\n positioner traces present, but none shows a line structure.")
return
name, arr, bounds, _segs, _spans = scored[0]
if index >= len(bounds):
print(f"\n the position trace has only {len(bounds)} line segments.")
return
a, b = bounds[index]
seg = np.asarray(arr[a:b], dtype=float)
tarr = times.get(name)
t = None
if tarr is not None and tarr.size >= b:
t = np.asarray(tarr[a:b], dtype=float) - float(tarr[a])
print(f"\n position trace for this line, from {name} ({seg.size} samples)")
print(f" {'samp':>4} {'t [s]':>10} {'position':>15} {'velocity':>14}")
for i in range(seg.size):
ts = f"{t[i]:10.4f}" if t is not None else f"{'-':>10}"
if i and t is not None and t[i] > t[i - 1]:
vel = f"{(seg[i] - seg[i - 1]) / (t[i] - t[i - 1]):14.6g}"
else:
vel = f"{'-':>14}"
print(f" {i:>4} {ts} {seg[i]:>15.6f} {vel}")
print("\n what the shape says:")
if seg.size < 3:
print(f" only {seg.size} samples in this line -- too few to separate the")
print(" mechanisms. The monitored stream is sampled far more sparsely")
print(" than the data; that is a limit of the file, not a result.")
return
typical = float(np.median(np.abs(np.diff(seg))))
if typical <= 0:
print(" the axis does not move at all in this line.")
return
# Progress along the *commanded* direction. Taking the direction from the
# segment itself fails exactly where it matters: if the return move falls
# inside the segment, first and last sample are the same point and the
# apparent direction is zero -- which would report a flyback as healthy.
if nominal is not None:
direction = np.sign(nominal[-1] - nominal[0]) or 1.0
travel_cmd = float(abs(nominal[-1] - nominal[0]))
else:
head = seg[: max(2, seg.size // 4)]
direction = np.sign(head[-1] - head[0]) or 1.0
travel_cmd = float(np.ptp(seg))
adv = direction * (seg - seg[0]) # distance covered, signed forward
stop_i = int(np.argmax(adv)) # first sample at maximum travel
point_step = travel_cmd / max(1, per_line - 1)
# If nothing ever advances, these samples are not a trace of the traverse and
# nothing about the traverse can be read out of them. Saying so is the whole
# point: the alternative is to describe the return move as if it were the
# line, which is exactly the wrong answer confidently delivered.
if float(adv.max()) <= 0.1 * travel_cmd:
print(f" NOT A TRACE OF THE LINE. The samples span"
f" {float(seg.min()):.6f} .. {float(seg.max()):.6f}"
f" -- {float(np.ptp(seg)) / travel_cmd:.0%} of the commanded travel, and")
print(" none of it in the commanded direction. They are taken at fixed")
print(" points in the scan loop -- after the burst and during the return --")
print(" not while the line runs. No mechanism can be named from this file;")
print(" it needs a faster monitored readout or a measurement on the motor.")
return
if t is not None and stop_i > 0 and t[stop_i] > 0:
v_ach = adv[stop_i] / t[stop_i]
line = f" velocity while moving : {v_ach:.6g}/s"
if v_cmd:
line += f" commanded {v_cmd:.6g}/s ratio {v_ach / v_cmd:.2f}x"
print(line)
if v_cmd and v_ach > 1.15 * v_cmd:
print(" Faster than commanded. The axis will not run this slowly, so it")
print(" finishes the traverse early and the remaining triggers fire into")
print(" a stage that is no longer moving. Note that the scan validates")
print(" the requested velocity against base_velocity and that check")
print(" passed -- so the reported floor understates the real one.")
after = adv[stop_i:]
back = float(adv[stop_i] - after.min()) # how far it came back afterwards
overshoot = float(adv[stop_i] - travel_cmd)
frozen = 0
for gap in np.abs(np.diff(after)):
if gap > 0.1 * typical:
break
frozen += 1
if frozen:
print(f" stationary for {frozen} of {seg.size - 1} intervals after reaching"
f" full travel at sample {stop_i}"
f" ({frozen / max(1, seg.size - 1):.0%} of the line)")
if back > 0.25 * travel_cmd:
print(f" MECHANISM: the trace returns {back:.4g} toward the line start before")
print(" the line ends -- the flyback begins while triggers are still firing,")
print(" so the tail of the line re-samples it backwards, not standing still.")
elif overshoot > 2 * point_step:
print(f" MECHANISM: overshoot -- ran {overshoot:.4g} past the commanded end,")
print(" then settled back. The deceleration ramp is computed from the")
print(" commanded velocity, so a faster actual velocity overruns the target,")
print(" and the parked rows show a place outside the imaged range.")
elif frozen:
print(" MECHANISM: arrived early and parked. The seam in the image is then a")
print(" change of blur and scale rather than of position: moving points")
print(" average over their travel, parked points do not.")
else:
print(" the axis moves throughout this line -- whatever marks the image, it")
print(" is not this line stopping early.")
def dump_scan_command(h5: h5py.File) -> None:
"""Print what the scan was actually asked to do.
This is the authoritative record: the file writer stores the whole scan
metadata under /entry/collection/metadata/bec, including `request_inputs`,
which holds the arguments the scan was called with. It beats notebooks,
logbooks and memory because it is written by the thing that ran the scan.
"""
roots = [p for p in ("/entry/collection/metadata/bec",
"/entry/collection/metadata") if p in h5]
if not roots:
print(" no metadata group found -- cannot recover the command from this file")
return
for root in roots[:1]:
print(f" scan request, from {root}:")
shown = 0
def show(name, obj):
nonlocal shown
if not isinstance(obj, h5py.Dataset):
return
try:
val = obj[()]
except Exception: # pylint: disable=broad-except
return
if isinstance(val, bytes):
val = val.decode(errors="replace")
if isinstance(val, np.ndarray):
val = val.tolist()
text = str(val)
if len(text) > 300:
text = text[:300] + " ..."
print(f" {name:<46} {text}")
shown += 1
# request_inputs first: that is the command itself
for sub in ("request_inputs", "scan_args", "scan_parameters"):
if sub in h5[root]:
h5[f"{root}/{sub}"].visititems(lambda n, o, s=sub: show(f"{s}/{n}", o))
for key in sorted(h5[root]):
if key in ("request_inputs", "scan_args", "scan_parameters"):
continue
node = h5[f"{root}/{key}"]
if isinstance(node, h5py.Dataset):
show(key, node)
if not shown:
print(" (nothing readable -- run --list to see the raw structure)")
def analyse(path: str, do_list: bool, force_lines=None, force_points=None,
show_command: bool = False, trace_lines=()) -> None:
print("=" * 80)
print(os.path.basename(path))
print("=" * 80)
with h5py.File(path, "r") as h5:
meta = read_metadata(h5)
interesting = {k: meta[k] for k in
("scan_name", "scan_number", "num_points", "num_lines",
"frames_per_trigger") if k in meta}
print(" metadata:", interesting if interesting else "(none of the expected keys found)")
if show_command:
dump_scan_command(h5)
return
if do_list:
# Per-point arrays first: the file is dominated by scalar configuration
# entries, and a naive dump buries the data under hundreds of them.
arrays, scalars = [], 0
for name, shape in walk_datasets(h5):
size = int(np.prod(shape)) if shape else 1
if size > 1:
arrays.append((size, shape, name))
else:
scalars += 1
arrays.sort(key=lambda t: (-t[0], t[2]))
print(f"\n {len(arrays)} array datasets ({scalars} scalars hidden):")
for size, shape, name in arrays[:80]:
print(f" {str(shape):>14} {name}")
if len(arrays) > 80:
print(f" ... and {len(arrays) - 80} more")
if not arrays:
print(" none -- this file holds no per-point data at all,")
print(" which would itself explain the empty analysis.")
return
signals = numeric_signals(h5)
if not signals:
print(" no 1-D numeric signals found -- re-run with --list")
return
lengths: dict[int, list[str]] = {}
for name, arr in signals.items():
lengths.setdefault(arr.size, []).append(name)
print(f" {len(signals)} numeric signals, in {len(lengths)} different lengths:")
for size in sorted(lengths, reverse=True):
names = sorted(lengths[size])
shown = ", ".join(names[:6]) + (f", +{len(names) - 6} more" if len(names) > 6 else "")
print(f" {size:>8} points : {shown}")
# Prefer the length shared by the most signals: per-point data outnumbers
# whatever else is in the file. Taking the longest signal is wrong when a
# monitor is sampled on its own clock.
# Prefer the length the metadata says is the grid. These files hold two
# populations: the MCS channels, which are the per-point scan data, and the
# monitored devices, sampled opportunistically while the line runs. The
# monitored set is the larger of the two by signal count, so picking the
# most common length lands on the wrong data.
grid_n = None
if meta.get("num_lines") and meta.get("frames_per_trigger"):
want = int(meta["num_lines"]) * int(meta["frames_per_trigger"])
if want in lengths:
grid_n = want
n = grid_n if grid_n else max(lengths, key=lambda k: len(lengths[k]))
how = "metadata"
if force_lines and force_points:
grid, n, how = (force_lines, force_points), force_lines * force_points, "--lines/--points"
else:
grid = infer_grid(n, meta)
if grid is None:
# No usable metadata: recover the geometry from the data itself.
how = "detected from the scanned axis"
for cand in sorted(lengths[n]):
lines_found = detect_lines(signals[cand][:n])
if lines_found and n % lines_found == 0:
grid = (lines_found, n // lines_found)
how = f"detected from {cand}"
break
if grid:
print(f" geometry: {grid[0]} lines x {grid[1]} points ({how})")
if grid is None:
print(f"\n could not work out the line geometry for {n} points.")
print(" factor pairs that could plausibly be (lines x points-per-line):")
pairs = [(d, n // d) for d in range(2, int(n**0.5) + 1) if n % d == 0]
pairs += [(b, a) for a, b in pairs]
plausible = sorted({p for p in pairs if 2 <= p[0] <= 2000 and 4 <= p[1] <= 20000})
if plausible:
for a, b in plausible[:12]:
print(f" {a:>6} lines x {b:>6} points")
if len(plausible) > 12:
print(f" ... and {len(plausible) - 12} more")
else:
print(f" none -- {n} does not factor usefully, so the per-point arrays")
print(" may be padded or truncated, which is itself worth knowing.")
print("\n Re-run with --lines N --points M once you know the scan geometry")
print(" (the scan command that produced it gives both).")
return
mon = monitored_names(h5)
motors = positioner_devices(h5)
print(f" positioners in file: {', '.join(sorted(motors)) or '(none found)'}")
# The two checks read different arrays. Intensity lives on the scan grid;
# motor position does not -- it is sampled while the line runs, at whatever
# rate the readout loop managed, so it has its own length and its own
# line boundaries.
check_line_motion_from_motor(signals, motors, meta, read_request_inputs(h5))
req = read_request_inputs(h5)
check_end_of_line(signals, grid, mon, None, req, motors)
check_stationary_from_data(signals, grid, motors, req)
if trace_lines:
times = numeric_timestamps(h5)
for idx in trace_lines:
trace_line(h5, signals, times, grid, meta, motors, idx)
def main(argv=None) -> int:
ap = argparse.ArgumentParser(description=__doc__.splitlines()[0])
ap.add_argument("files", nargs="*", help="master .h5 files")
ap.add_argument("--scan", type=int, action="append", default=[], help="scan number")
ap.add_argument("--root", default="/sls/x12sa/data/p23080/raw/data",
help="data root used with --scan")
ap.add_argument("--list", action="store_true", help="dump the file structure and stop")
ap.add_argument("--command", action="store_true",
help="print the scan request recorded in the file, and stop")
ap.add_argument("--lines", type=int, help="override number of lines")
ap.add_argument("--points", type=int, help="override points per line")
ap.add_argument("--trace-line", type=int, action="append", default=[], metavar="N",
help="dump the point-by-point trace of line N (0-based); repeatable")
args = ap.parse_args(argv)
paths = list(args.files)
for num in args.scan:
found = resolve_scan(args.root, num)
if found:
paths.append(found)
else:
print(f"scan {num}: no master file under {args.root}", file=sys.stderr)
if not paths:
ap.error("give at least one file or --scan")
for path in paths:
try:
analyse(path, args.list, args.lines, args.points, args.command,
args.trace_line)
except Exception as exc: # pylint: disable=broad-except
print(f" FAILED on {path}: {type(exc).__name__}: {exc}")
print()
return 0
if __name__ == "__main__":
sys.exit(main())