Skip to content

Better Atoms: a functional form for every atom, on a new physical-unit vocabulary - #83

Open
atrabattoni wants to merge 48 commits into
devfrom
feature/better-atoms
Open

Better Atoms: a functional form for every atom, on a new physical-unit vocabulary#83
atrabattoni wants to merge 48 commits into
devfrom
feature/better-atoms

Conversation

@atrabattoni

Copy link
Copy Markdown
Contributor

The atoms half of 0.2.9, on top of the tiles backend and the obspy engine
already merged by #82. docs/release-notes.md carries the 0.2.9 summary for a
user; this is the shape of the change.

One operation, two faces

The rule is now general rather than a property of a few routines: every atom
has a user-friendly functional form, and every function an atom form.
The
function is what one writes — xd.resample(da, 50.0) applies straight away —
and seeding it with ... yields the atom behind it, xd.resample(..., 50.0),
which composes with >>, pickles, and streams. The same code therefore serves
a slice in memory and an archive that does not fit in one; nothing has to be
rewritten to move between them. xdas.atoms.as_function generates the
functional form of any atom class, so the rule holds for atoms you write too.

Ordinary numpy expressions join in under the same seed: with ... in the data
slot, 20 * np.log10(np.abs(atom)) appends absolute → log10 → multiply to
the pipeline instead of computing. Fan-in (two atoms in one expression) raises
at the line that wrote it rather than silently computing.

A new default vocabulary, away from SciPy's parameters

The top-level processing functions no longer expose the arguments of the SciPy
routine underneath. Compare:

before (xdas.signal) now (xdas)
xs.resample(da, num=8192) xd.resample(da, 50.0) — a target rate in Hz
xs.decimate(da, q=2) xd.decimate(da, 50.0) — likewise
xs.filter(da, freq, btype, corners) xd.filter(da, (1.0, 10.0)) — a corner pair, None opening an end
xs.medfilt(da, kernel_dim) xd.medfilt(da, {"time": 0.5}) — seconds, not samples

An output sample count, a decimation factor and a normalised frequency all
describe the array; a rate, a corner and a window length describe the
measurement. Only the latter keep their meaning when the sampling rate
changes, which is what lets one pipeline be defined once and applied to
whatever it is given — and what makes the pipeline a faithful record of the
processing rather than of the file it was written against.

The new set, in its functional form:

xd.filter      xd.resample   xd.decimate    xd.integrate
xd.stft        xd.detrend    xd.taper       xd.differentiate
xd.hilbert     xd.medfilt    xd.rechunk     xd.sliding_mean_removal
xd.annotate    xd.trigger    xd.pick

each with an atom behind it (Filter, Resample, Decimate, … in
xdas.atoms.tasks, Annotate/Trigger/Picker in xdas.atoms.ml and
xdas.atoms.detect). The SciPy-shaped functions remain available in
xdas.signal.

The exact machine-parameter atoms — LFilter, SOSFilter, DownSample,
UpSample, and the new fused Polyphase — move down into an explicitly
expert layer, xdas.atoms.kernel. They stay public and importable from
xdas.atoms, but they are no longer the entry point: the vocabulary above
designs them from the data at the first call, and a user who never needs to
name filter coefficients never meets them. (Polyphase fuses
upsample/filter/downsample into one upfirdn pass, which is what makes
resampling 2.6–8.7× faster and keeps float32 from being promoted.)

The execution contract

call() now maps one input chunk to zero or more output chunks and
flush() drains what is buffered — which is what makes rechunking, windowing
and reductions expressible as atoms at all.

Every stateful atom judges the seams of its own input from the chunk
coordinates: state carries across a continuous seam, is flushed and restarted
at a gap or a rate change, and raises on an overlap. Eager calls split gappy
input into the same runs. "Chunked equals eager" is therefore a property to
test rather than to hope for, and xdas.testing.assert_chunk_invariant tests
it — quantified over re-chunkings at non-divisor sizes and over injected gaps.
It runs in CI for every stateful atom shipped. Atoms that cannot answer chunk
by chunk (the fft functions along the transformed dimension, the
whole-record routines) refuse instead of answering wrong.

process()

Now a method on every atom and the dispatch boundary for both ends of a run.
Sources: array, virtual array, path, directory, glob, DataCollection,
xdas.watch(dir), a ZeroMQ address, or any chunk iterable. Sinks: directory,
.csv, address, configured writer, or None to accumulate. Collections are
walked leaf by leaf with each result labelled by its tree path, and the
"memory_limit" config entry turns the two silent footguns — an eager call on
a huge virtual array, an accumulation that outgrows RAM — into errors naming
the streaming alternative. That limit defaults to a quarter of the memory the
process can use: the smaller of the machine's physical memory and any cgroup
limit (container, batch scheduler), so it is 4 GiB on a laptop and 250 GiB on a
terabyte node rather than one number that is wrong on both. The historical
process(atom, loader, writer) form is unchanged.

Picking

Annotate drives a SeisBench model from what its weight set declares —
window overlap, stacking, blinding, preprocessing — rather than from what the
architecture suggests, and overlaps device compute with transfers. Trigger
keys its thresholds on phase labels. Picker assembles the chain a weight set
describes, so xd.pick(dc, model) answers a whole network tree with one flat
pick table.

Also here

xdas.stack (collapse a collection level into a dimension, lazily, with
tolerance snapping and inner/outer joins), pool="processes" chunk IO through
Ray's shared-memory object store, DataFrame as a first-class collection leaf,
sel on string and categorical coordinates, and the numpy-dispatch fix where
registered defaults overrode explicit arguments.

Breaking

MLPicker is deprecated in favour of Annotate, kept as an alias until 0.4;
its output is now laid out sample-last. The re-exported Trigger defaults to
dim="time" rather than "last".

Label selection on a string coordinate raised UFuncTypeError before any
selection began: the overlap guard differenced the coordinate values,
and numpy has no subtract loop for string dtypes. The monotonicity check
now goes through the pandas index, which is dtype-generic; is_unique
keeps it strict, since pandas considers repeated values monotonic
increasing.

The guard itself now covers only ordered look-ups. A slice resolves its
bounds by searching and `method` searches for a neighbour, so both need
an axis whose values increase and stay guarded. Naming a label is a hash
look-up, well defined in any order, so it goes straight through — a
categorical axis such as ["P", "S", "N"] is selectable by label, and a
list returns its labels in the requested order. Ambiguous labels remain
the coordinate's business: to_index raises on them.
The dispatch wrapper applied its registration defaults after binding the
caller's arguments, so a registered `axis=-1` overwrote whatever the
caller passed: np.cumsum(da, 0) accumulated along the last axis,
positional or keyword alike. Defaults now fill in only the parameters
the caller did not bind, before numpy's own signature defaults.
DataArrayWriter concatenated with concat's "first" default, which is
only right when the chunked dimension leads the output: (distance, time)
chunks from a time-chunked source were stacked along distance. The
writer now takes the dimension as a dim argument; left unsaid it keeps
the "first" default.
The tree rebuild wrapped every non-DataArray leaf in DataArray(...), so
a collection holding tables raised on repr. A DataFrame now passes
through untouched, and the mapping repr asks the value what it is
rather than assuming everything non-array is a subtree.
The inverse of combine_by_coords, which concatenates along an existing
dimension: xd.stack(dc, "channel") turns each station's traces into one
(channel, time) array keyed by the level's own keys. The new dimension
is named after the level it collapsed (dim= chooses otherwise),
everything below the level is merged in lock-step, and tile-backed
leaves stay lazy. Leaves that disagree on their other coordinates raise
naming what disagreed; join="inner" trims to the common part and
join="outer" pads with NaN.

Coordinate agreement is judged on the sampling grid within a tolerance
rather than exactly: real networks record a channel a nanosecond off,
and exact comparison made outer joins silently interleave two grids into
mostly-NaN samples. The default tolerance is 1/100 sample where a
nominal interval is declared, tolerance=False keeps the strict
comparison, and true interleaving is refused rather than padded.

Also fixes concat opening a new dimension over VirtualSource-backed
arrays raising TypeError: whether the result can stay virtual is decided
before expand_dims, which no virtual source can follow.
Reading a chunk goes through one virtual layout, hence one h5py call
under the global HDF5 lock, so loader threads never decode concurrently
and extra threads only contend; compression on the write side contends
the same way. DataArrayLoader and DataArrayWriter now accept
pool="processes": each worker holds its own HDF5 lock, receives the
manifest of its chunk (a sliced virtual array, kilobytes) and reads its
own files.

The chunks come back through Ray's shared-memory object store rather
than a pickle pipe: written once by the worker, mapped zero-copy by the
parent. A pipe caps out well below memory bandwidth and spends parent
CPU on deserialization the signal chain needs — measured end to end on a
compressed ZFP archive at 16 workers, ingest 137 -> 1378 MiB/s and
egress 151 -> 1556 where a pickling pool lost end to end. max_workers is
enforced by parking submissions beyond it; Ray initializes lazily, an
already initialized runtime is left untouched, and shutdown leaves it up
because it is a session-wide resource.

The price of zero-copy is immutability: chunk data arrives read-only,
which atoms honor by allocating their outputs — pinned by a test running
a stateful pipeline over read-only chunks, ray installed or not. Ray is
the optional extra xdas[ray]; pool="threads" stays the default.
The core grows the vocabulary pipelines are built from:

- compose() and >>/>>=/da >> atom, with value semantics everywhere:
  composing never mutates an operand, so intermediate pipelines stay
  usable on their own. This also fixes the atomized Sequential-append
  aliasing bug (it mutated the input and returned None).
- Ufunc tracing under the ... seed: applying a numpy ufunc to an atom
  appends the operation to the pipeline instead of computing it, so
  20 * np.log10(np.abs(atom)) is a pipeline. The traced surface is
  ufuncs exactly; an expression involving two atoms (fan-in) raises at
  the line that wrote it, never bails silently into computation.
  Atom equality becomes identity so atoms can live in sets and traces
  cannot be confused by ==.
- as_function(cls) generates the function form of an atom class: a
  lowercase function taking the data first, eager on data, returning
  the configured atom on ..., extending a pipeline on an atom.
  atomized dispatches classes to it and keeps its function behaviour.
- fresh(): a stateless clone, config shared by reference (a model is
  never deep-copied), nested atoms recursed. initialized now recurses
  into nested atoms, so a pipeline holding an uninitialised filter no
  longer reports itself ready.
- A private whole-record guard: functions that need the whole record
  along their working dimension are marked at their definition site and
  refuse chunked execution along that dimension with a pointed error,
  resolving first/last aliases against the data before comparing.

DataArray defers to foreign __array_ufunc__ implementations so
da >> atom dispatches to the atom rather than being computed.
LFilter, SOSFilter, DownSample and UpSample move to xdas.atoms.kernel,
the expert layer of exact stateful primitives with machine parameters
(taps, integer factors) whose meaning depends on the sampling rate; they
stay importable from xdas.atoms.

The layer is born with Polyphase: upsample, FIR filter and downsample in
one scipy.signal.upfirdn pass, computing only the output samples that
survive the decimation and never materialising the zero-stuffed signal.
FIRFilter gains up=/down= and applies its taps through Polyphase
(designed at the upsampled rate, energy-compensated), and ResamplePoly
collapses its upsample/filter/downsample trio to that one child atom —
2.6x on a decimation by two along distance, 8.7x on 62.5 -> 50 Hz on a
254 MiB chunk. Chunked calls carry the filter memory and the output-grid
phase across chunks, so splitting the input does not change the result.

Taps are cast down to the data precision (float32 in, float32 out
instead of an lfilter promotion doubling downstream memory traffic), and
a target rate the coordinate resolution cannot represent exactly
declares its residual drift as jitter on the output coordinate. The
output-phase arithmetic is pinned against scipy's own upfirdn and
against the explicit kernel chain across up/down combinations and chunk
cuts.
xdas.atoms.tasks is the public processing vocabulary: Filter (a
(low, high) corner pair in Hz with None opening one end, ftype iir/fir,
zerophase), Decimate and Resample (target rate in Hz), Integrate and
Differentiate (stateful, carrying their seam state across chunks), plus
whole-record detrend, taper, hilbert, sliding_mean_removal and medfilt
with kernel lengths in physical units. Machine parameters (taps,
factors) are designed from the data at the first call and live in the
kernel layer, so a pipeline keeps its meaning when the sampling rate
changes.

Decimate applies its anti-alias taps through the polyphase kernel
(antialias.down carries the factor), never filtering at the full rate to
throw most of the result away. Zero-phase IIR filtering has no causal
streaming form, so that atom refuses chunked execution along its
dimension, as do the whole-record functions, each marked at its
definition site.

Every task atom gets a function form exported at the top level:
xd.decimate(da, 50.0) applies eagerly, xd.decimate(..., 50.0) returns
the atom, and passing an atom extends a pipeline.
Every stateful atom now judges the seams of its own input stream from
the chunk coordinates: a continuous chunk carries state across, a gap or
rate change flushes the previous run and reinitialises (redesigning
coefficients for the new rate), a backward overlap raises, and
on_discontinuity="reset"|"raise" makes strict runs opt-in. Incoming
chunks are split at internal discontinuities so state never crosses a
gap, and chunked stateful processing requires a regular coordinate along
the chunked dimension, raising with a pointer to to_regular() instead of
silently carrying state across unverifiable seams.

call() follows the transducer contract: one chunk in, zero or more
chunks out. The new flush() lifecycle drains buffered samples at the end
of the stream, at every seam, and at the end of every eager call;
Sequential.flush cascades codec-drain style, process() drains the
pipeline, and writers drop empty chunks. This fixes chunked DownSample
dropping its trailing samples when the length is not a multiple of the
factor. Reductions fall out of the contract (accumulate in call, emit at
flush), Atom.iter_chunks exposes the manual chunk loop as a plain
generator, and the new Rechunk kernel atom (function form xdas.rechunk)
restores a target chunk cadence without ever merging across a
discontinuity.

Eager calls auto-split gappy input into runs and re-join the outputs
with the gaps kept in the coordinates — announced, not silent: a warning
states how many discontinuities the source has and that state is flushed
and reset at each, named by the source's start so a collection walk
reports every leaf. The first/last dimension aliases are resolved
against the data before any comparison with the chunked dimension, so a
kernel built with its documented default no longer skips allocating its
seam state; UpSample survives one-sample chunks.

The commutation invariant — concat(atom(split(da, anywhere))) equals
atom(da) for arbitrary split points, tails included — is now in the test
suite for the stateful vocabulary.
Chunk-safety is a claim an atom makes about itself;
xdas.testing.assert_chunk_invariant is the evidence. It runs a pipeline
once eagerly and once streamed chunk by chunk and asserts the two agree
— shapes, values, coordinates, and pick tables compared as sets of rows
since eager and chunked walks order them differently. coord_atol admits
an explicitly declared sub-sample coordinate drift where rational
resampling reconstructs its grid segment by segment.

The invariant quantifies over cuts and gaps: the same stream is
re-chunked at derived non-divisor sizes so the boundaries land
elsewhere (cuts=, explicit dicts accepted), and inject_gaps places real
discontinuities in the input first so seam resets are exercised at
chunk boundaries that do not line up with them. A negative control is
part of the test suite: a pipeline that is not chunk-invariant fails
loudly rather than passing.
The spectral vocabulary joins the task-atom route: STFT takes its window
length and hop in physical units, both snapped — the window to the next
fast FFT size of the target so transforms stay efficient whatever the
sampling rate, the hop to a whole number of samples — so the actual grid
can differ from the request, and the docstring says so. scaling= chooses
"spectrum" (peak amplitudes) or "psd", so np.abs(stft)**2 composes to
an exact spectrogram; an expert nfft zero-pads the windowed frames.

Frames start at the first sample and advance by the hop; only fully
computable frames are ever emitted. Chunk by chunk, the unconsumed tail
is buffered across chunks and dropped at gaps and at the end of the
stream, so chunked processing emits exactly the frames of the eager
transform and no frame ever spans a discontinuity — held under
assert_chunk_invariant over cuts and gaps. One-sided spectra for real
data, centered two-sided for complex; scipy.signal.ShortTimeFFT does the
design and scaling internally. xd.stft is the function form.

The xdas.fft functions (fft, rfft, ifft, irfft) now declare whole-record
semantics: used as atoms in a chunked pipeline they raise along the
transformed dimension instead of silently computing one transform per
chunk, and Partial resolves the {input_dim: output_dim} mapping form of
dim, so transforming along another dimension than the chunked one keeps
working.
Tests for the corners the main suites walked past: in-place ufuncs trace
out of place, out= to a foreign atom and non-call ufunc methods raise,
right_shift with the atom on the left is an ordinary traced ufunc;
fresh() recurses into nested class atoms and the refusal helper is
conservative on unresolvable aliases and checks the keys of kernel
dicts; empty chunks are skipped, a one-sample chunk of a sampled
coordinate inherits the stream's rate at the seam, single-sample
streams have nothing to judge, first/last aliases resolve on eager
calls, dimensionless atoms map collections eagerly, folds work chunked,
unconcatenatable outputs fall back to a sequence; writers drop empty
chunks; the cut derivation stops when no new size exists; UpSample
survives a one-sample record and Polyphase an empty one.
process() becomes a method on every atom: pipeline.process(da,
out="results/") infers both ends. Sources dispatch on the input value
(get_source): an in-memory DataArray runs eagerly (or chunk by chunk
with chunks=), a virtual one streams through a DataArrayLoader whose
chunks="auto" aligns boundaries to the storage blocking (tile extents,
per-file extents) merged up to a byte budget, a path/directory/glob
opens with open_mfdataarray (multi-acquisition collections chain one
loader per run), tcp:// subscribes over ZeroMQ through a scheme
registry, and any iterable of chunks is consumed as is — the source
contract is iteration plus optional chunk_dim/nbytes/unbounded. Sinks
dispatch on the out spec crossed with the first output chunk
(get_writer), deferring writer creation to what the pipeline actually
emits: directories store DataArray chunks joined along the chunked
dimension (or SDS archives for Streams), *.csv appends DataFrames,
tcp:// publishes, out=None accumulates and returns the joined result,
writer instances pass through, and no empty outputs are ever created.
The historical process(atom, loader, writer) form keeps working.

Realtime is named: xd.watch(path, engine=...) wraps RealTimeLoader, and
unbounded sources get streaming semantics — no byte total on the
monitor, a clean KeyboardInterrupt that flushes the pipeline and
returns the writer result, and until= to stop at a coordinate value
(inclusive, truncating the last chunk).

Discontinuities are announced at the boundary too: a chunked source
warns once upfront with the count read off its coordinate (no data
touched), and a realtime source — which cannot be inspected upfront —
warns at each seam as it arrives.

The new memory_limit configuration entry (default 8 GiB) guards the two
footguns: an eager call on a huge virtual array and an out=None
accumulation that outgrows the limit both raise with the estimated size
and a pointer to .process(out=...).
Tests for the three corners of process() the main suites walked past: a
source whose chunked coordinate is dense gets no upfront discontinuity
scan, a realtime chunk that does not carry the chunked dimension leaves
the seam information untouched, and a realtime one-sample chunk of a
sampled coordinate inherits the stream's rate so the seam after it is
still judged correctly.
Everything a SeisBench picker does is a property of the weight set rather
than the architecture, so tests/fakemodel.py makes that an executable
contract instead of a warning in a document. FakeModel is a real
WaveformModel subclass with no weights and a closed-form forward pass;
WEIGHT_SETS holds five archetypes drawn from the cached PhaseNet
metadata, between them spanning ENZ/ZNE/Z12H, 3 and 4 input channels, the
NPS/PSN label flip, 50 and 100 Hz, blinding and overlap declared and
absent, per-phase thresholds declared and absent, and no filter, a flat
filter and a per-channel one. tests/conftest.py exposes the factory as
a fake_model fixture; import the module directly to parametrise on it.
Everything a SeisBench picker reads is a property of the weight set, so
the atom now reads it there: the overlap (in SeisBench's fraction-or-
samples form), the stacking rule, the blinding and every annotate
argument come off the model instance, overridable per call, with the
model's own annotate_batch_pre/post driving normalisation and blinding.
The component dimension is found by its labels ending with distinct
letters of component_order — never by its name — with the flexible
horizontal matching, and component_strategy spans SeisBench's range:
auto, clone, pad, a named slot, strict.

The output keeps the input's order among the batch dimensions but comes
out sample-last, (..., 'phase', dim): the characteristic function of one
phase of one channel is contiguous, which is the layout its consumers
reduce along. The end-aligned final window SeisBench appends is emitted
at flush(), so the output spans the input and stays chunk-invariant. A
model whose annotate_batch_post breaks the (batch, samples, classes)
stacking contract is named instead of surfacing as a broadcast error.

Chunked along its own dimension the sliding window carries across
chunks; chunked along another dimension the carry-over would leak one
chunk's time tail onto the next chunk's other lanes, so such a chunk is
a whole record run from a fresh state and settled on the spot — and the
component dimension is refused as a chunk axis, since the model reads
every component of a window at once.

MLPicker and xdas.mlpicker stay as DeprecationWarning aliases until 0.4.
The value pin holds: the numbers today's MLPicker produced on the DAS
layouts are unchanged, asserted through dimension names so they outlive
the layout change.
SeisBench's annotate_stream_pre filters the waveforms before resampling
them with whatever the weight set declares in filter_args/filter_kwargs;
skipping it feeds the network something it was not trained on. The
_model_filter builder understands both declared forms — flat, applied to
everything, and per channel, one glob pattern each, which SeisBench's
own DAS wrapper refuses — and translates them exactly: obspy filters
with corners=4, zerophase=False and a Butterworth in second-order
sections, which is what Filter defaults to. A zero-phase declaration
warns and doubles the order (exact zero-phase IIR has no causal
streaming form) and no corner may sit above half the Nyquist of the
model's own rate, both as SeisBench concedes.

The per-channel form is _ChannelFilter, a Filter subclass that filters
everything and restores the channels its pattern does not match: the
channel dimension is found by its labels exactly as Annotate finds it,
a pattern matching nothing is a silent no-op (stream.select semantics),
and ftype='fir' is refused since compensating the group delay on the
coordinate would leave the untouched channels on the wrong samples.

Both names are private: the stage exists for the picker assembly to
come, not as a user-facing bandpass — compose Filter yourself for that.
The agreement with obspy is pinned bitwise across every band that has a
Filter equivalent.
…ns, flush

Trigger gains the three pieces the picker needs from it, in its new home
xdas/atoms/detect.py alongside the rest of the detection vocabulary.

thresh accepts a mapping keyed on the phase coordinate Annotate
produces, as well as the scalar it took before. A label the mapping does
not name gets an infinite threshold and therefore never triggers, which
is what lets Annotate keep emitting the noise class: nothing downstream
has to slice it out of the characteristic function. Keying on the label
rather than on its position is a correctness requirement — the label
order belongs to the weight set and flips between them — so the numba
kernel now takes one threshold per lane instead of two scalars.

coords accepts scalar (0-d) coordinates and emits them as constant
columns, so a pick table carries its network/station/location identity.
Its default becomes 'auto' — identity first, measurement last: scalars,
then the other dimension coordinates, then the picked dimension — so
the columns do not depend on the input's dimension order. The columns
are resolved once at initialize, which is also what lets flush build a
table without the chunk it no longer has.

flush() closes the triggers still open at the end of a run, as
obspy.trigger_onset does at the end of an array; the last pick of a
record used to be dropped. Chunked along a dimension other than dim,
none of the state carries — such a chunk is a whole record of other
lanes, run from a fresh state and closed on the spot. Atom._join learns
to concatenate DataFrame chunks, since an eager call can now answer
with the picks of its call plus those of its flush and must still
return one table.

xdas/trigger.py becomes a compatibility module: find_picks stays
verbatim, Trigger is re-exported (its dim default is now 'time', not
'last'), and the module is imported eagerly so that the new lowercase
twin xd.trigger — which joins the top level with the other function
forms — is never shadowed by a later import of the module.
A collection walk now carries the tree path down to the leaves instead
of rebuilding provenance afterwards. Each leaf's result is labelled with
the path it was reached by as it is produced — one column per named
level, leading the table, filled with that level's key — so a pick found
under IA / DBNFM / -- comes out with its network, station and location
before its time and its value. Producing the labels at the leaf rather
than folding them in on the way back up is what will let a streaming
walk hand a leaf straight to a sink and still know whose it was.

Atoms may then declare a merge(results) hook folding those labelled
results. It is undefined on Atom, so an atom returning arrays sees
nothing change and its tree is rebuilt as before; Trigger.merge is a
plain concat, which is all it takes once the columns are there, so
xd.trigger(dc, ...) answers a whole network with one flat table; and
Sequential.merge delegates to the last stage declaring one. merge=False
opts out and returns the labelled tree. A column with two sources — the
tree key and a scalar coordinate of the same name — stays a single
column the tree path fills, warning on a genuine disagreement.
Positional levels contribute their index, and the flushed tail of a
folded sequence is attributed to the last element it came out of.

Annotating a collection no longer needs a prior xd.stack: Atom.gather
is a hook the walk consults on every mapping level before descending —
return the level collapsed to an array, or None (the default) to map
over its leaves. Annotate implements it, because what counts as a
component is a property of the model: Z12H and ENZ disagree about which
channels group together. It collapses through xd.stack, so the
structural checks, the grid snapping and the error messages are shared,
and tolerance= reaches them. Sequential.gather delegates to the first
stage that claims the level, so the gather happens once and before the
first stage runs — the per-channel filter a weight set ships selects
??H out of Z12H and can only do that once the channels are a dimension.

Recognition is deliberately conservative, since folding a station level
into a component axis would silently destroy the distinction between
stations: both the level's name (COMPONENT_LEVELS, or components=) and
its keys must resolve, a key must look like a whole channel code, and
the keys must agree on their length and band code. Keys that resolve to
nothing leave the level alone, so a DAS collection whose spatial axis
is called channel walks straight past; keys that resolve only in part
raise, naming the conflict. Pre-stacking keeps working: annotate(dc)
equals annotate(xd.stack(dc, 'channel')).
Picker(model) assembles the whole pipeline SeisBench's model.classify
runs — the weight set's own preprocessing filter when it ships one,
Resample to the weight set's own rate (not always 100 Hz: diting runs
at 50), Annotate, and Trigger with the weight set's own per-phase
thresholds — from the weight set and nothing else. Two pickers built on
one model class can differ in stage count, sampling rate and
thresholds, which is the point: everything a SeisBench picker does is a
property of the weights.

The thresholds come from _model_thresholds: one entry per picked label
as model_pick_labels resolves them — a model declaring a phases subset
(the EQTransformer family) picks exactly that subset, anything else
picks every label but the noise class — each looked up as SeisBench
does: the call wins, else the weight set's default_args, else the
model's documented default for that key, else the *_threshold
catch-all, else 0.3. Values pass through faithfully, including
iquique's P_threshold of 1.12, which simply never fires.

Annotate reads its labels through model_phases — labels=None falls back
to positional labels exactly as WaveformModel._predictions_to_stream
does, a callable is refused by name — and its class count now comes
from the labels rather than model.classes, which counts only the
picking head: EQTransformer sets it to 2 while labelling three outputs,
and sizing the buffers on it made the whole family unrunnable.

Being a Sequential rather than a factory function, a picker keeps
everything a pipeline can do — >> composes it, repr shows its stages,
it pickles, picker.process streams it — and inherits Annotate.gather
and Trigger.merge, so xd.pick(dc, model) answers a whole network tree
with one flat table. The dimension aliases first/last are refused: the
picks are annotated with coordinates, which are named.
atom.process(dc, ...) now walks a DataCollection exactly as atom(dc)
does, so atom.process(dc, out=None) == atom(dc). Anything else makes the
streaming form second-class precisely where streaming matters most: a
leaf too large to call eagerly at all is what the memory guard already
refuses, pointing at .process(da, out=...).

The walk mirrors Atom._walk step for step. Each mapping level is offered
to gather before anything is chunked — so a channel level becomes a
component dimension and the stacked array streams as one thing — then
recursed into; sequence levels fold element by element, the state
carrying across and the atom flushed once at the end, its tail
attributed to the last element; and each leaf streams through the
existing single-source path, with chunks= and until= applying per leaf.
One atom instance takes the leaves one at a time, reset between them,
because an atom holding a model either saturates the CPU or holds a lot
of device memory.

Sinks gain a rule per destination, and the produce-time path columns are
what make them work — every output chunk is labelled on its way to the
writer, never reconstructed afterwards:

- out=None accumulates per leaf and merges, giving the eager result;
- a *.csv, a URL or a ready writer instance is shared: every leaf
  appends to one table, the path columns keeping the rows apart, and
  the walk answers with that one result;
- a directory fans out, one subdirectory per leaf mirroring the tree
  path, since a directory of netcdf chunks describes one stream — a
  folded sequence is one stream and so writes to one directory.

merge= is accepted on process and popped walk-level, as __call__ does.
A DataSequence handed over as a collection now folds rather than
streaming as one concatenated result; a glob or directory that opens to
a sequence is still a single source, and get_source(sequence) asks for
the same of a collection in hand.

Verified while porting: xd.stack keeps a tile-backed leaf tile-backed,
so the gather of an ObsPy-style collection reads nothing.
The process loop stays serial and CPU atoms stay internally parallel;
only Annotate is asynchronous, behind a small bounded queue — the same
max_buffers pattern DataArrayLoader and DataArrayWriter already use.

Each completed window's reduction now allocates a fresh tensor (it
never views the circular stack), so on CUDA its device-to-host transfer
can be issued immediately — pinned staging buffer, non-blocking copy,
an event marking completion — while the sliding window moves on. call()
emits only the outputs whose transfer has completed, in order (the 0..n
contract already allows late emission); flush() drains whatever is
still in flight, so nothing survives the end of a run and chunk
invariance holds unchanged. At most max_buffers transfers are left
pending (default 2, 0 restores fully synchronous emission), which is
what bounds the staging memory. On the CPU there is nothing to wait
for: the queue completes on arrival and behavior is exactly the
synchronous reference.

Together with the H2D staging path that was already pinned and
non-blocking, the CPU keeps preparing and feeding windows while the GPU
computes and while results cross back, with no executor machinery: no
thread-per-stage, no stage budget, no second level of parallelism to
oversubscribe the first.
A new user-guide page next to the DAS pipeline material: xd.open then
xd.pick over a real 8-station network day, what the picker builds from a
weight set (three stages for original, four for obs), why thresholds and
classes are keyed by label rather than by position, and the same walk on
a DAS collection eagerly and through process(..., out=...). Its code is
not executed at build time — SeisBench and its weights are not
documentation dependencies — but every output shown was produced by
running it on feature/atoms, whose picking pipeline this branch
reproduces stage for stage.

The page states the one real SeisBench deviation as what it is: obspy's
Trace.resample defaults to window='hann' applied in the frequency
domain, which halves the amplitude at half the input Nyquist, where the
polyphase filter is flat. Fed the same resampled data, the two agree
pick for pick.

Annotate, Trigger and Picker get their own API section rather than
sitting under signal processing, the deprecated MLPicker drops out of
the listing, and the function forms annotate/pick/trigger join the
top-level list. Release notes gain the walk labelling/merge/gather
block, the process() collection walk and the async GPU queue.
The slow real-weight MLPicker tests took two small facts with them when
the fake-model suite replaced them: randn_wavefronts had no contract of
its own (shape, coordinates, seeded reproducibility — it is the DAS
synthetic the picking walkthrough builds on), and no test constructed
an Annotate without naming a device, so the CUDA-if-available default
went unexercised.
A stage that changes the number of samples along a dimension has to say
what became of the *other* coordinates attached to it. Polyphase and
UpSample copied them through untouched, so a decimated array kept its
full-length station codes: the dimension had 5 000 samples and the
coordinate naming them still had 10 000 values, and every lane was
labelled with the code of the lane at its own index — the first half of
the fibre, shifted by one channel each. Picking a DAS acquisition
answers with those codes, so the picks came out attributed to the wrong
channels; found by running the ABYSS pipeline against its archived
reference, whose stations sit one every 30.6 m where these sat one
every 15.3 m.

The labels now follow the samples: output k is drawn from input
k * down / up, and that is the input whose label it takes. They carry no
group-delay shift — a label names a source sample, where the dimension
coordinate names a position, which is what the delay compensation is
about — and the mapping depends on the sampling grid alone, so chunking
cannot move a label. DownSample was already right: isel subsamples every
coordinate attached to the dimension.

Note for the record: xd.concat drops non-dimensional coordinates along
the concatenated dimension, so re-joining chunked output loses them.
That is pre-existing (a plain split/concat round trip loses them too)
and left alone here; the streaming path is unaffected, since each chunk
reaches the sink already labelled.
A chunk of one sample declares no sampling interval of its own, so the
seam it stored had no rate to compare the next chunk against and the
judgement was skipped entirely: a gap right after such a chunk carried
state across silently. The seam is now judged on the rate of whichever
side knows one.

Two smaller repairs on the same layer: `Sequential.fresh` rebuilt the
clone by calling the constructor, which a pipeline that assembles its
own stages (a `Picker`) does not accept, and the `out=` guard of the
ufunc tracer compared tuples, so `out=` an array re-entered tracing
instead of refusing.

Along the way: the near-duplicate of the chunk joiner folded into the
one it duplicated, the one-call flush helper inlined into the cascade
that calls it, and the docstrings of the layer corrected.
`TileArray.root` is the path the whole archive relocates by, so it is
part of the class rather than an attribute the constructor happens to
set: declared and documented, it is also what the API page can point
at. Plus the typos the docs build walked past.
The atoms and processing pages still taught `Sequential([...])`,
`Partial` and hand-built loader/writer pairs — correct, but no longer
what one would write. They now open on the function forms and `>>`,
on physical parameters, and on `process(source, out=...)` with the
source and sink tables, keeping the explicit form as what to reach for
when the ends need configuring.

The picking page's annotation repr predated the sample-last layout, the
obspy page duplicated `xd.stack` by hand and called the legacy engine
an alias of the new one, and the API pages were missing
`trim_overlaps`, `DataCollection.select` and the collection hooks.
Rewritten from the development log it had become into what a 0.2.8 user
meets: a section each for the tiles backend and the atoms rework, the
rest of the new API in one paragraph apiece, and the breaking section
reserved for what one actually runs into. Everything else moved to
improvements or refactoring, and the details left to the docs.
The streaming page described the publisher and the subscriber but not
what one does with them now: naming an address on either end of
`process`, following a directory with `watch`, stopping on `until` or
on a keyboard interrupt. The FAQ still answered the chunked-filter
question with the kernel layer. A stray dangling reference and the one
ambiguous cross-reference left in the build go with them.
…king

Three ways the same promise was broken. Rounding an output sample's
source position could walk past the last input sample of the chunk, and
the clip that pulled it back landed on a different sample than the eager
call picked: flooring is in range by construction, so the label of an
output no longer depends on where the stream was cut.

The "first"/"last" alias was resolved against the *coordinates*, so a
dimension carrying none kept the literal alias, never matched the
chunked dimension, and the seam state was silently never allocated —
chunked output quietly differed from eager. Resolution now goes against
the dimensions, which a bare axis still has.

And a whole-record function whose `dim` defaults to the last dimension
refused *every* chunked dimension, since an unknown working dimension is
refused against all of them: the fft functions could not be used in a
pipeline chunked along another dimension, which is exactly what they
document. The default is now declared where the guard can read it.
The sink path was taken from the annotation path, which drops the levels
that have no name — the default for a plain `DataCollection({...})`. Every
leaf then wrote into the same directory, each writer restarting its
numbering, and reading one back gave another's data. The keys the walk
descended by now address the sink, whether or not they name a column.

Three more on the same path: `until=` was silently ignored on the eager
source (the truncation only existed in the chunked loop); the
accumulation guard sized chunks with `nbytes`, which tables and streams
do not have, so the one case that most needs a limit — a collection walk
accumulating into one pick table — was unbounded; and folding a sequence
asked for chunks larger than its last element, which raises where the
sibling path clamps.
A coordinate attached to the picked dimension — a tag, a code, anything
naming the samples rather than the lanes — was indexed with the absolute
sample number of the run against a single chunk's values, which raised
once the run was longer than a chunk. Those coordinates are now
accumulated over the run like the dimension coordinate they sit on, so a
trigger's onset still names something when the chunk it was found in is
gone. `coords=None` no longer raises on a dimension carrying no
coordinate: there is nothing to annotate with there.

`Annotate` shipped its companion coordinates at the length of the input
rather than of the emitted chunk, which the assembly then dropped, and
swallowed a stream shorter than one model window in silence where the
eager call says so. `MLPicker` gets its old signature back: the new
third positional argument is `components`, so the documented
`MLPicker(model, dim, device)` was setting the wrong one.
The point is not that a few routines gained physical parameters: it is
that the norm is now one operation with two faces — a user-friendly
functional form, and the atom behind it that composes and streams — and
that the default vocabulary has moved off SciPy's parameters onto the
quantities of the measurement. Those two lead, with the whole functional
roster named.

Which demotes the exact machine-parameter atoms to what they are: an
expert layer one should not have to meet, said at the end of the section
rather than as a headline about polyphase throughput. And the leftover
"New Features" heading, now holding nothing but the obspy engine, stack,
trim_overlaps and select, is named for what it holds.
A guard against footguns has to sit far enough above what one
legitimately loads at once that it only fires on a mistake, and no fixed
number does that across a laptop and a node with a terabyte: 8 GiB is a
quarter of the one and a hundredth of a percent of the other, refusing
ordinary work on the big machine while still being generous on the small.

The default is now a quarter of what the *process* can use — the smaller
of the machine's physical memory and any cgroup limit, so a container or
a batch allocation is respected rather than the hardware behind it — and
the error says what the limit is, since it is no longer a number one can
recite.
ZeroMQ drops what a publisher sends to a peer whose subscription has not
reached it yet, and being connected is not being subscribed. There is no
way to ask from the receiving end, so the tests slept and hoped — and hung
outright whenever the sleep turned out to be too short, which on a loaded
machine it eventually is.

A publisher now answers each new subscription with a greeting, in passing
as it streams. That is what the XPUB welcome message is for, but ZeroMQ
withholds it until the application reads the subscription it answers, which
nothing here ever did: submitting a packet now does, so the greeting has in
fact never been delivered before this. It is the ASN header, which is why a
subscriber joining a running interrogator can now be told the shape of the
stream — and skip forward to it, rather than choke on the first packet it
happens to land on. `wait_until_subscribed` returns on that greeting, and
returning is proof that nothing published from then on will be missed.

None of this asks anything of a real-time publisher, which streams whether
or not anyone listens and is never held up by its audience. Replaying a
recording is the one case that needs the other end to wait, since nothing a
subscriber does can hold back a stream already under way, and
`wait_for_subscribers` is for that alone.

Both subscribers take a timeout, so a stream gone quiet raises where it used
to hang. The ZMQ tests keep no sleep, cover joining a live flux, and no
longer end in a doctest whose publisher thread raced the reader for a
rebound global — that one hung the suite for the fifty minutes it took to
notice.
Ray earned its place by getting chunks across a process boundary without
pickling them, which is what caps a process pool: a loaded chunk crossing
back was serialized, transferred and deserialized, well below memory
bandwidth and on the parent's CPU. It brought a scheduler, a raylet, an
object store sized against the machine's memory and session state under
/tmp for a job whose chunks are bounded in size and whose flow is already
back-pressured -- everything that store is built to survive, and none of
what it is needed for.

An arena of shared memory does the same work in one module. The parent
cuts one /dev/shm file into fixed slots and hands them out; a worker writes
its chunk once into a slot and the parent maps the same pages, and a chunk
on its way to a writer is staged the same way. Reuse is the point:
allocating a block per chunk makes the kernel zero and fault every page
again, which is why the obvious spelling of this measures slower than the
pipe it replaces. What crosses is a ShmRef, the shape/dtype duck type
DataArray already accepts from a virtual array, so nothing else in the path
changes. Chunks arrive read-only, as they did from the object store.

Slots come back when the chunk that owns them is collected, so a streaming
consumer turns a handful of them forever. Anything that will not fit -- an
oversized chunk, an exhausted arena, a result that is not an array -- takes
the ordinary pickle path, which is slower and never wrong.

The pages die with the mappings that hold them, so any crash frees them.
What could outlive a run is the arena's name, unlinked at shutdown, again
by a finalizer at exit, and swept at the start of the next run if its owner
is gone.

Two things had to be taught. Loky recycles a worker whose resident memory
grows 300 MB past its baseline, and shared pages count toward it, so the
arena is discounted from that check or workers quit and respawn mid-run.
And the whole design rests on unlinking a mapped file and on signal zero as
a liveness probe, so it is POSIX-only; elsewhere the pool pickles as it did
before, which is what it did everywhere without ray installed.

pool="processes" is unchanged as a name, an interface and a contract. What
goes is the optional dependency.
The tasks these tests submit are defined beside them, and cloudpickle sends a
function defined in an importable module by reference -- so the worker was
being asked to import `tests`. It could, but only by accident: `python -m
pytest` puts the working directory on the path and loky passes that path on to
its workers. Run as `uv run pytest`, as the CI does, nothing puts the project
root on the path -- the import mode this suite uses names the module
`tests.test_pools` without making it importable -- and the worker failed to
unserialize the call.

Registering the module for pickle by value sends the tasks themselves instead,
which is what already happens for anything defined in a script. Nothing about
what is under test changes: the code the pool submits in earnest lives in
`xdas`, which a worker can always import.
Shutting a pool down dropped both the arena's name and this process's handle
on it, but the finalizer that stands in when nobody shuts it down only dropped
the name. The handle stayed in the module's table for good, and with it the
pages, so a process that opened pools in a loop and let them fall out of scope
accumulated one whole arena per pool -- twelve of them after twelve pools, in
the run that turned this up.

The mapping is still not closed, only forgotten: chunks the caller is holding
are views on it and keep it alive by themselves, which is the same promise
shutdown already made.

The module also claimed a crash of any kind frees the pages. That is true of
anything that runs Python on the way out, an interrupt included, but not of a
run killed outright: loky's workers never notice their parent has gone, and
they hold what they have mapped until somebody kills them too. It reproduces
with a bare loky pool and no shared memory anywhere near it, so the note now
says what actually happens rather than what the arena alone would do.
@codecov

codecov Bot commented Aug 12, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (e0f6e9f) to head (1a9d064).
⚠️ Report is 1 commits behind head on dev.

Additional details and impacted files
@@             Coverage Diff             @@
##               dev       #83     +/-   ##
===========================================
  Coverage   100.00%   100.00%             
===========================================
  Files           45        49      +4     
  Lines         6252      8450   +2198     
  Branches      1067      1524    +457     
===========================================
+ Hits          6252      8450   +2198     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

TdmsReader opens the file in __init__ and closes it in __exit__, but a
reader whose __init__ raises is never handed to the "with" that would
have closed it. Every auto-detection sniff of a file that is not TDMS
went down that path and leaked a descriptor.
None of the four ZMQ publishers and subscribers ever closed anything: a
socket kept its file descriptor, and its context an I/O thread, for the
life of the process, and the two ASN classes dropped their context on
the floor at construction. Closing is now part of the interface, shared
by all four through ZMQEndpoint: as a context manager where the endpoint
has a scope, with close() where it does not, and on garbage collection
for one that is merely dropped. A publisher process() opened itself from
a "tcp://" spec is closed by process(); one passed in stays the
caller's.
simplefilter("error") turned every warning into a failure, including the
ones the test has no opinion about, which is why the blocks it guarded
were emptied rather than kept. Naming the category instead says what the
test means -- no discontinuity was announced -- and leaves unrelated
warnings alone, so the four blocks that had been left standing with
nothing in them assert something again.

Also covers the last four paths the suite never reached: a coordinate
riding the picked dimension, a chunked record shorter than one model
window, an unchunked source truncated by "until", and a shared table no
leaf ever wrote to.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant