From 36959f399789a425ce13368ceb4c91ad7c64482b Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 21:21:19 +0200 Subject: [PATCH 1/8] Subdivide a relation into pieces rather than at cuts Channel-level selection needs the same split conform_to_inventory does, but it keeps runs of channels rather than partitioning a span: a row it matches nothing of leaves the spool, and one it matches entirely passes through untouched. Cuts cannot say either. So build_subdivision_plan now takes the pieces themselves, and the snapping which turns cuts into pieces is subdivision_pieces beside it. A piece is modified when it differs from its own row rather than when its row had cuts, which says the same thing about a cut row and the right thing about a kept one. Divide the native types when snapping. Converting each operand to float seconds first rounds three times where this rounds once, and a comment claiming CoordRange._get_index's fudge factor was too weak was simply wrong -- it divides natively and is fine. One rounding is still a rounding, though: measured over 400k random float triples, the ratio lands a hair high on 5% of on-grid cuts and a hair low on 6% of off-grid ones, so both correction loops earn their place along a distance axis even though a datetime one cannot reach them. --- dascore/core/spool.py | 116 +++++++++++++++---------- dascore/utils/chunk_plan.py | 108 +++++++++++++++-------- tests/test_proc/test_proc_inventory.py | 2 +- tests/test_utils/test_chunk.py | 108 +++++++++++++++++++---- 4 files changed, 237 insertions(+), 97 deletions(-) diff --git a/dascore/core/spool.py b/dascore/core/spool.py index 827e692d..35a44f48 100644 --- a/dascore/core/spool.py +++ b/dascore/core/spool.py @@ -49,6 +49,7 @@ build_chunk_plan, build_subdivision_plan, samples_adjusted_envelopes, + subdivision_pieces, ) from dascore.utils.display import get_dascore_text, get_nice_text from dascore.utils.docs import compose_docstring, get_docstring @@ -154,56 +155,57 @@ def _report_unconformed(rows: pd.DataFrame, on_unresolved: str) -> None: warnings.warn(msg, UserWarning, stacklevel=3) -def _check_one_acquisition(source_rows: pd.DataFrame, epochs) -> None: +def _refuse_rows(source_rows: pd.DataFrame, reasons, summary: str) -> None: """ - Refuse the patches whose acquisition changes partway through. + Raise naming the patches an inventory verb cannot handle, and why. + + Every refusal these verbs make has the same shape — a few patches out + of an archive, each with its own particular; naming them is the whole + value, since the fix is nearly always to one file or one inventory + entry. `reasons` holds a particular per row and None where the row is + fine, so a caller judges rows without also formatting them. + """ + named = [ + f"{path} ({reason})" + for path, reason in zip(source_rows["source_path"], reasons, strict=True) + if reason is not None + ] + if not named: + return + msg = f"{len(named)} patch(es) {summary}: {_first_few(named)}." + raise PatchError(msg) + + +def _acquisition_conflicts(epochs) -> list: + """ + Name the patches whose acquisition changes partway through. Subdividing cannot reconcile these the way it reconciles a change of optical path — see `RowEpochs.conflict` for why. `on_unresolved` does not wave them through either: the inventory describes such a patch twice rather than not at all. """ - conflicted = [ - f"{path} at {row.conflict}" - for path, row in zip(source_rows["source_path"], epochs, strict=True) - if row.conflict is not None - ] - if not conflicted: - return - msg = ( - f"{len(conflicted)} patch(es) span a change of acquisition, which " - f"subdividing cannot reconcile: {_first_few(conflicted)}. Select the " - "side you want, or correct the inventory." - ) - raise PatchError(msg) + return [None if x.conflict is None else f"at {x.conflict}" for x in epochs] -def _check_subdividable(source_rows: pd.DataFrame, rows: pd.DataFrame, cuts) -> None: +def _unsubdividable(rows: pd.DataFrame, pieces, name: str) -> list: """ - Refuse a patch which must be split but states no sampling interval. - - The pieces are found on the patch's own sample grid, which its time - step is the only description of: a piece ends one step short of - where the next begins. Without a step both pieces would have to - claim the boundary instant, and envelopes are inclusive, so a sample - landing there would appear in both. Duplicating a sample is a worse - answer than saying so — the caller asked for metadata - reconciliation, not for the data to be restructured. + Name the patches which must be split but state no sampling interval. + + The pieces are found on the patch's own sample grid, which its step + is the only description of: a piece ends one step short of where the + next begins. Without a step both pieces would have to claim the + boundary value, and envelopes are inclusive, so a sample landing + there would appear in both. Duplicating a sample is a worse answer + than saying so — the caller asked for the metadata to be reconciled, + not for the data to be restructured. """ - bad = [ - f"{path} at {row_cuts[0]}" - for path, step, row_cuts in zip( - source_rows["source_path"], rows["time_step"], cuts, strict=True - ) - if row_cuts and (pd.isnull(step) or not step) + return [ + f"at {row_pieces[0]}" + if row_pieces and (pd.isnull(step) or not step) + else None + for step, row_pieces in zip(rows[f"{name}_step"], pieces, strict=True) ] - if not bad: - return - msg = ( - f"{len(bad)} patch(es) must be subdivided at an epoch boundary but " - f"state no time step to find their samples with: {_first_few(bad)}." - ) - raise PatchError(msg) def _unstated(values) -> np.ndarray: @@ -1364,7 +1366,6 @@ def conform_to_inventory( >>> mixed = dc.spool([patch, other]).attach_inventory(inventory) >>> assert len(mixed.conform_to_inventory(on_unresolved="drop")) == 1 """ - from dascore.io.index.planned import derived_catalog # noqa: PLC0415 from dascore.proc.inventory import ( # noqa: PLC0415 _NO_EPOCHS, resolve_row_epochs, @@ -1387,7 +1388,12 @@ def conform_to_inventory( if columns is None else resolve_row_epochs(new._inventory, *columns) ) - _check_one_acquisition(source_rows, epochs) + _refuse_rows( + source_rows, + _acquisition_conflicts(epochs), + "span a change of acquisition, which subdividing cannot " + "reconcile; select the side you want, or correct the inventory", + ) described = np.array([x.described for x in epochs], dtype=bool) if not described.all(): _report_unconformed(source_rows[~described], on_unresolved) @@ -1396,16 +1402,36 @@ def conform_to_inventory( if not any(cuts): # nothing to subdivide: a filter is the whole job return new._restrict_to_rows(kept["_patch_id"].to_numpy()) sources = source_rows[described].reset_index(drop=True) - _check_subdividable(sources, kept, cuts) + _refuse_rows( + sources, + _unsubdividable(kept, cuts, "time"), + "must be subdivided at an epoch boundary but state no time step " + "to find their samples with", + ) + return new._subdivided( + sources, kept, subdivision_pieces(kept, cuts, "time"), "time" + ) + + def _subdivided(self, sources, rows, pieces, name: str) -> Self: + """ + Return the spool whose patches are the given pieces of these rows. + + The pieces are a plan rather than a rewritten relation, so the + outputs *are* the contents rows — `len` and `get_contents` stay + exact — and loading goes through the machinery which already + trims a member on extraction. + """ + from dascore.io.index.planned import derived_catalog # noqa: PLC0415 + catalog = derived_catalog( source_rows=sources, - plan=build_subdivision_plan(kept, cuts, "time"), - parent=new._catalog, + plan=build_subdivision_plan(rows, pieces, name), + parent=self._catalog, merge_kwargs={}, mode="chunk", - origin_path=new.spool_path, + origin_path=self.spool_path, ) - return new._new_from_catalog(catalog) + return self._new_from_catalog(catalog) def _with_inventory(self, inventory, on_unresolved, valid, method) -> Self: """ diff --git a/dascore/utils/chunk_plan.py b/dascore/utils/chunk_plan.py index 92936a74..efa7025f 100644 --- a/dascore/utils/chunk_plan.py +++ b/dascore/utils/chunk_plan.py @@ -1044,18 +1044,15 @@ def _snapped_cuts(cuts, start, step) -> list: exactly on opens the piece that cut opens, and a cut falling between two samples leaves the earlier one behind. """ - span = to_float(step) out = [] for cut in cuts: - # The ratio is inexact — `to_float` rounds both operands and then - # the quotient — so it only starts the search, and the grid it - # describes settles it. `CoordRange._get_index`'s fudge factor is - # not enough here: at a million samples the error outgrows any - # fixed tolerance, and an index one too high puts the boundary - # sample in the epoch *before* the boundary, which is the one - # place it must not go. Comparing the values themselves is exact, - # and each loop steps at most once for any plausible error. - index = math.ceil(to_float(cut - start) / span) + # Dividing the native types rounds once, where converting each to + # float seconds first would round three times — enough to lift an + # on-grid cut above its own index, putting the boundary sample in + # the piece before the boundary. One rounding is still a rounding, + # so the ratio only starts the search and comparing values on the + # grid settles it; each loop steps at most once. + index = math.ceil((cut - start) / step) while start + (index - 1) * step >= cut: index -= 1 while start + index * step < cut: @@ -1073,28 +1070,32 @@ def _snapped_cuts(cuts, start, step) -> list: return sorted(out) -def build_subdivision_plan(df: pd.DataFrame, cuts, name: str) -> ChunkPlan: +def _dim_columns(df: pd.DataFrame, name: str) -> tuple[str, str, str]: + """Return the envelope columns of one dimension, checked present.""" + columns = (f"{name}_min", f"{name}_max", f"{name}_step") + assert set(columns).issubset(df.columns) + return columns + + +def subdivision_pieces(df: pd.DataFrame, cuts, name: str) -> list[list[tuple]]: """ - Build a plan splitting each row of a relation at its own cut values. + Return the inclusive pieces each row's cuts divide it into. - Unlike a chunk plan, no row ever meets another: each output is one - contiguous piece of exactly one source row, so the pieces of a row - partition its samples — none is dropped, none is duplicated, and a - row with no cuts passes through as a single unmodified output. This - is what an operation which must re-describe patches rather than - restructure them (`Spool.conform_to_inventory`) needs. + The pieces of a row partition its samples: none is dropped and none + is duplicated, whatever instant a cut falls on, because each piece + ends one step short of where the next begins. Parameters ---------- df - The relation to subdivide; one row per patch. + The relation being subdivided; one row per patch. cuts One sequence of cut values per row, in the row's own units, each above that row's minimum and no greater than its maximum — a cut on the maximum yields a one-sample final piece. Order does not matter. A cut opens a new piece at the first sample at or after it, so a row with `n` distinct cuts becomes at most `n + 1` - outputs. + pieces. name The dimension being subdivided. @@ -1104,28 +1105,67 @@ def build_subdivision_plan(df: pd.DataFrame, cuts, name: str) -> ChunkPlan: callers must reject a null or zero step themselves — they are the ones which can name the file it came from. """ - min_name, max_name = f"{name}_min", f"{name}_max" - step_name = f"{name}_step" - assert {min_name, max_name, step_name}.issubset(df.columns) - # One entry per row, even where it is empty: a short sequence would - # drop the rows past its end from the plan, and so from the spool. + min_name, max_name, step_name = _dim_columns(df, name) assert len(cuts) == len(df) - df = _ensure_patch_id(df).reset_index(drop=True) - positions, lows, highs, modified = [], [], [], [] + df = df.reset_index(drop=True) + out = [] for position, row_cuts in enumerate(cuts): start, stop = df.at[position, min_name], df.at[position, max_name] # Envelopes are value-ordered whatever the coordinate's # orientation, so the grid is walked by the step's magnitude. step = abs(df.at[position, step_name]) assert not len(row_cuts) or (not pd.isnull(step) and step) - grid = _snapped_cuts(row_cuts, start, step) - bounds = [start, *grid] - for index, low in enumerate(bounds): - high = bounds[index + 1] - step if index + 1 < len(bounds) else stop + bounds = [start, *_snapped_cuts(row_cuts, start, step)] + out.append( + [ + (low, bounds[index + 1] - step if index + 1 < len(bounds) else stop) + for index, low in enumerate(bounds) + ] + ) + return out + + +def build_subdivision_plan(df: pd.DataFrame, pieces, name: str) -> ChunkPlan: + """ + Build a plan cutting each row of a relation into its own pieces. + + Unlike a chunk plan, no row ever meets another: each output is one + contiguous piece of exactly one source row. A row handed the whole of + its own envelope passes through as a single unmodified output, and + one handed no pieces at all leaves the relation — which is how a + selection along the dimension drops a patch none of whose samples it + keeps. This is what operations which re-describe patches rather than + restructure them (`Spool.conform_to_inventory`, inventory-backed + channel selection) need. + + Parameters + ---------- + df + The relation to subdivide; one row per patch. + pieces + One sequence of inclusive `(low, high)` envelopes per row, in the + row's own units and already on its sample grid. They must be + disjoint and ascending within a row; `subdivision_pieces` builds + them from cut values, and a mask over the row's samples gives + them directly. + name + The dimension being subdivided. + """ + min_name, max_name, step_name = _dim_columns(df, name) + # One entry per row, even where it is empty: a short sequence would + # drop the rows past its end from the plan, and so from the spool. + assert len(pieces) == len(df) + df = _ensure_patch_id(df).reset_index(drop=True) + positions, lows, highs, modified = [], [], [], [] + for position, row_pieces in enumerate(pieces): + whole = (df.at[position, min_name], df.at[position, max_name]) + for piece in row_pieces: positions.append(position) - lows.append(low) - highs.append(high) - modified.append(bool(grid)) + lows.append(piece[0]) + highs.append(piece[1]) + # A piece covering its whole row *is* the row, and saying so + # is what lets it load without a trim. + modified.append(tuple(piece) != whole) ids = np.arange(len(positions), dtype=np.int64) # Outputs are not file rows: source bookkeeping stays on the members, # and the dimension's structural identity described the whole row. diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index 3562ab21..71244465 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -2358,7 +2358,7 @@ def test_an_acquisition_change_raises(self, patch, inventory, off_grid_boundary) named = re.escape(spool.get_contents()["source_path"].iloc[0]) for policy in ("raise", "warn", "drop"): with pytest.raises( - PatchError, match=f"change of acquisition.*{named} at .*" + PatchError, match=f"change of acquisition.*{named} \\(at .*" ): spool.conform_to_inventory(on_unresolved=policy) diff --git a/tests/test_utils/test_chunk.py b/tests/test_utils/test_chunk.py index 32b384ec..487723b8 100644 --- a/tests/test_utils/test_chunk.py +++ b/tests/test_utils/test_chunk.py @@ -15,6 +15,7 @@ _normalize_chunk_units, build_chunk_plan, build_subdivision_plan, + subdivision_pieces, ) from dascore.utils.time import to_timedelta64 @@ -539,25 +540,30 @@ def test_mixed_dtypes_upcast(self, sized_df): assert plan.params["size"]["partitions"][0]["dtype"] == "float64" -def _one_row_df(start, step, samples): - """A one-row relation over an evenly sampled time dimension.""" +def _one_row_df(start, step, samples, name="time"): + """A one-row relation over one evenly sampled dimension.""" return pd.DataFrame( { - "time_min": [start], - "time_max": [start + step * (samples - 1)], - "time_step": [step], + f"{name}_min": [start], + f"{name}_max": [start + step * (samples - 1)], + f"{name}_step": [step], "_patch_id": [0], } ) -class TestBuildSubdivisionPlan: - """Splitting each row of a relation at its own cut values.""" +def _cut_plan(df, cuts, name="time"): + """The plan a relation's per-row cuts produce, through both helpers.""" + return build_subdivision_plan(df, subdivision_pieces(df, cuts, name), name) + + +class TestSubdivisionPieces: + """Turning each row's cut values into pieces of its sample grid.""" def test_uncut_rows_pass_through(self, contiguous_df): """A row with no cuts is one unmodified output of the same span.""" df = contiguous_df.assign(_patch_id=np.arange(len(contiguous_df))) - plan = build_subdivision_plan(df, [()] * len(df), "time") + plan = _cut_plan(df, [()] * len(df)) assert len(plan.outputs) == len(df) assert not plan.members["_modified"].any() assert (plan.outputs["time_min"].values == df["time_min"].values).all() @@ -567,7 +573,7 @@ def test_pieces_partition_the_samples(self): start, step = STARTTIME, np.timedelta64(10, "ms") df = _one_row_df(start, step, 100) cut = start + step * 40 - plan = build_subdivision_plan(df, [(cut,)], "time") + plan = _cut_plan(df, [(cut,)]) assert len(plan.outputs) == 2 first, second = plan.outputs.iloc[0], plan.outputs.iloc[1] assert second["time_min"] == cut @@ -579,15 +585,16 @@ def test_an_on_grid_cut_lands_on_its_own_sample(self, index): """ A cut sitting exactly on a sample opens the piece at that sample. - The index is found by a float division, which is inexact — at - these step and index combinations the naive quotient rounds just - above the integer — so an unchecked `ceil` puts the boundary - sample one piece too early. Only a grid comparison settles it. + Dividing the native types keeps this exact for a time + coordinate — datetime64 arithmetic is integer nanoseconds — where + converting each operand to float seconds first rounds three + times and puts the boundary sample one piece too early at these + step and index combinations. """ start, step = STARTTIME, np.timedelta64(17060962, "ns") df = _one_row_df(start, step, index + 10) cut = start + step * index - plan = build_subdivision_plan(df, [(cut,)], "time") + plan = _cut_plan(df, [(cut,)]) assert plan.outputs["time_min"].iloc[1] == cut def test_a_cut_a_hair_past_a_sample_opens_at_the_next(self): @@ -602,15 +609,48 @@ def test_a_cut_a_hair_past_a_sample_opens_at_the_next(self): start, step = STARTTIME, np.timedelta64(1, "D").astype("timedelta64[ns]") df = _one_row_df(start, step, 1100) cut = start + step * 1000 + np.timedelta64(1, "ns") - plan = build_subdivision_plan(df, [(cut,)], "time") + plan = _cut_plan(df, [(cut,)]) assert plan.outputs["time_min"].iloc[1] == start + step * 1001 + def test_an_on_grid_float_cut_is_not_pushed_past_its_sample(self): + """ + A float grid rounds where a datetime one cannot, and both ways. + + `1.0 + 0.1` is the textbook case: the difference back out is + `0.10000000000000009`, so the ratio is a hair above 1 and `ceil` + answers 2 — the sample the cut lands exactly on would open the + piece *before* it. Comparing against the grid is what settles + it, and a distance axis makes this the ordinary case rather than + the exotic one. + """ + start, step = 1.0, 0.1 + df = _one_row_df(start, step, 10, name="distance") + cut = start + step + assert (cut - start) / step > 1 # the misleading ratio itself + plan = _cut_plan(df, [(cut,)], name="distance") + assert plan.outputs["distance_min"].iloc[1] == cut + + def test_a_float_cut_past_a_sample_still_opens_at_the_next(self): + """ + The increment, on a float grid: a cut one ulp past a sample. + + The subtraction cannot represent the gap, so the ratio says the + cut sits exactly on the sample. It does not — the sample is + behind it, and belongs to the piece before the cut. + """ + start, step = 0.1, 0.25 + df = _one_row_df(start, step, 10, name="distance") + cut = np.nextafter(start + step, np.inf) + assert (cut - start) / step == 1 # the ratio cannot see the gap + plan = _cut_plan(df, [(cut,)], name="distance") + assert plan.outputs["distance_min"].iloc[1] == start + 2 * step + def test_a_cut_on_the_maximum_yields_one_sample(self): """A row ending exactly on a boundary keeps its last sample apart.""" start, step = STARTTIME, np.timedelta64(10, "ms") df = _one_row_df(start, step, 50) cut = df["time_max"].iloc[0] - plan = build_subdivision_plan(df, [(cut,)], "time") + plan = _cut_plan(df, [(cut,)]) assert len(plan.outputs) == 2 assert plan.outputs["time_min"].iloc[1] == plan.outputs["time_max"].iloc[1] @@ -618,4 +658,38 @@ def test_one_cut_sequence_per_row(self): """A short sequence would silently drop rows from the plan.""" df = _one_row_df(STARTTIME, np.timedelta64(10, "ms"), 10) with pytest.raises(AssertionError): - build_subdivision_plan(pd.concat([df, df]), [()], "time") + subdivision_pieces(pd.concat([df, df]), [()], "time") + + +class TestBuildSubdivisionPlan: + """Building a plan from pieces given directly, as a selection does.""" + + def test_a_row_given_no_pieces_leaves(self): + """This is how a selection drops a patch it keeps no sample of.""" + df = _one_row_df(0.0, 1.0, 10, name="distance") + plan = build_subdivision_plan(df, [[]], "distance") + assert not len(plan.outputs) + assert not len(plan.members) + + def test_a_whole_row_piece_is_unmodified(self): + """A piece which *is* its row loads without a trim.""" + df = _one_row_df(0.0, 1.0, 10, name="distance") + whole = (df["distance_min"].iloc[0], df["distance_max"].iloc[0]) + plan = build_subdivision_plan(df, [[whole]], "distance") + assert len(plan.outputs) == 1 + assert not plan.members["_modified"].any() + + def test_disjoint_pieces_become_separate_outputs(self): + """Two runs of kept channels are two patches, with a hole between.""" + df = _one_row_df(0.0, 1.0, 10, name="distance") + plan = build_subdivision_plan(df, [[(0.0, 2.0), (7.0, 9.0)]], "distance") + assert len(plan.outputs) == 2 + assert plan.outputs["distance_max"].iloc[0] == 2.0 + assert plan.outputs["distance_min"].iloc[1] == 7.0 + assert plan.members["_modified"].all() + + def test_one_piece_sequence_per_row(self): + """A short sequence would silently drop rows from the plan.""" + df = _one_row_df(0.0, 1.0, 10, name="distance") + with pytest.raises(AssertionError): + build_subdivision_plan(pd.concat([df, df]), [[]], "distance") From c765c2fe3a7041c2ecd2b3e9aff3a56efe7c303d Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 21:23:37 +0200 Subject: [PATCH 2/8] Add Patch.unselect, the complement of Patch.select Channel-level trimming needs to keep what a selection would have removed, and a patch is the one place a range complement makes sense: it can have samples taken out of its middle, where a spool would have to cut every patch into the pieces on either side. That is the property Spool.unselect refuses coordinates for, seen from the other side. It asks select itself which samples it would keep and inverts the answer, so a selector cannot come to mean one thing in select and another in its complement. Naming several coordinates complements each on its own: the true complement of a block is an L, which no array can hold. --- dascore/core/patch.py | 1 + dascore/proc/coords.py | 74 +++++++++++++++++++++++++ docs/tutorial/patch.qmd | 26 +++++++++ tests/test_proc/test_proc_coords.py | 84 +++++++++++++++++++++++++++++ 4 files changed, 185 insertions(+) diff --git a/dascore/core/patch.py b/dascore/core/patch.py index 32bbb042..1ccc8f3c 100644 --- a/dascore/core/patch.py +++ b/dascore/core/patch.py @@ -403,6 +403,7 @@ def get_patch_name(self, *args, **kwargs) -> str: # --- processing funcs select = dascore.proc.select + unselect = dascore.proc.unselect order = dascore.proc.order correlate = dascore.proc.correlate diff --git a/dascore/proc/coords.py b/dascore/proc/coords.py index 4cbb2573..09f4fecd 100644 --- a/dascore/proc/coords.py +++ b/dascore/proc/coords.py @@ -566,6 +566,80 @@ def select( return patch.new(data=data, coords=new_coords) +@patch_function(history=None) +def unselect( + patch: PatchType, *, copy=False, relative=False, samples=False, **kwargs +) -> PatchType: + """ + Return the patch outside a selection. + + The complement of [`Patch.select`](`dascore.Patch.select`): it takes + the same selectors and removes exactly the samples that selection + would have kept. + + Parameters + ---------- + patch + The patch object. + copy + If True, copy the resulting data. This is needed so the old + array can get gc'ed and memory freed. + relative + If True, unselect ranges are relative to the start of coordinate, if + positive, or the end of the coordinate, if negative. + samples + If True, the query meaning is in samples. + **kwargs + Used to specify the coordinate on which data are unselected. + + Examples + -------- + >>> import dascore as dc + >>> from dascore.examples import get_example_patch + >>> patch = get_example_patch() + >>> + >>> # Drop meters 50 to 300, keeping what lies outside them. + >>> outside = patch.unselect(distance=(50, 300)) + >>> + >>> # Drop the first ten distance samples. + >>> trimmed = patch.unselect(distance=(..., 10), samples=True) + + Notes + ----- + - Removing a range from the middle of a coordinate leaves a hole in + it, so the result is no longer evenly sampled and the coordinate + becomes a monotonic array. That is exactly what + [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses + coordinates for: at spool level the complement of a range is a hole + in every patch rather than a choice between patches. + + - Each named coordinate is complemented on its own. Selecting on two + coordinates keeps the samples in both ranges, and everything + outside that is an L rather than a block, which no array can hold — + so `unselect` removes each named range instead, which is the part + of the complement that is expressible. + """ + if invalid_coords := set(kwargs) - set(patch.coords.coord_map): + invalid_list = sorted(invalid_coords) + valid_list = sorted(patch.coords.coord_map) + msg = ( + f"Coordinate(s) {invalid_list} not found in patch coordinates: " + f"{valid_list}" + ) + raise PatchCoordinateError(msg) + complements = {} + for name, value in kwargs.items(): + coord = patch.coords.coord_map[name] + # Asking select itself which samples it would keep is what stops + # the two from drifting: one selector cannot come to mean + # different things in select and its complement. + _, indexer = coord.select(value, relative=relative, samples=samples) + kept = np.zeros(len(coord), dtype=bool) + kept[indexer] = True + complements[name] = ~kept + return patch.select(**complements, copy=copy) + + @patch_function(history=None) def order( patch: PatchType, *, copy=False, relative=False, samples=False, **kwargs diff --git a/docs/tutorial/patch.qmd b/docs/tutorial/patch.qmd index cf9eea52..dc969689 100644 --- a/docs/tutorial/patch.qmd +++ b/docs/tutorial/patch.qmd @@ -438,6 +438,32 @@ sub_patch = patch.select(distance=np.array([0, 12, 10, 9]), samples=True) assert len(sub_patch.get_array('distance')) == 4 ``` +## Unselect + +[`Patch.unselect`](`dascore.Patch.unselect`) is the complement of `select`: it takes the same selectors and removes exactly the samples `select` would have kept. + +```{python} +import dascore as dc + +patch = dc.get_example_patch() + +# Everything outside meters 50 to 200. +outside = patch.unselect(distance=(50, 200)) + +# Together the two account for every channel, and share none. +inside = patch.select(distance=(50, 200)) +assert len(outside.get_array("distance")) + len(inside.get_array("distance")) == 300 +``` + +Removing a range from the middle leaves a hole, so the coordinate is no longer evenly sampled: + +```{python} +assert patch.get_coord("distance").step is not None +assert outside.get_coord("distance").step is None +``` + +That is why [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses coordinates. A patch can have samples removed from its middle; at spool level the complement of a range would be a hole in every patch rather than a choice between patches. When several coordinates are named, each is complemented on its own — the true complement of a block is an L, which no array can hold. + ## Order Order is similar to [`Patch.select`](`dascore.Patch.select`), but will re-arrange data to the order specified by a value array. This may also cause parts of the patch to be duplicated. diff --git a/tests/test_proc/test_proc_coords.py b/tests/test_proc/test_proc_coords.py index deff814e..4f26bccb 100644 --- a/tests/test_proc/test_proc_coords.py +++ b/tests/test_proc/test_proc_coords.py @@ -644,6 +644,90 @@ def test_single_sample_retains_step(self, random_patch): assert time.step == orig.step +class TestUnselect: + """Keeping what a selection would have removed.""" + + def test_complements_select(self, random_patch): + """Together the two account for every sample, and share none.""" + selector = (50, 200) + kept = random_patch.select(distance=selector).get_array("distance") + dropped = random_patch.unselect(distance=selector).get_array("distance") + whole = random_patch.get_array("distance") + assert not set(kept) & set(dropped) + assert sorted([*kept, *dropped]) == sorted(whole) + + def test_interior_range_leaves_a_hole(self, random_patch): + """ + The property which makes a range complement wrong for a spool. + + A patch can have samples removed from its middle; a spool would + have to cut every patch into the pieces on either side. + """ + out = random_patch.unselect(distance=(50, 200)) + coord = out.get_coord("distance") + assert coord.step is None + values = out.get_array("distance") + assert not ((values >= 50) & (values <= 200)).any() + assert values.min() < 50 < 200 < values.max() + + def test_data_follows_the_coordinate(self, random_patch): + """The rows removed are the rows the selection would have kept.""" + axis = random_patch.dims.index("distance") + values = random_patch.get_array("distance") + out = random_patch.unselect(distance=(50, 200)) + wanted = random_patch.data.take( + np.flatnonzero(~((values >= 50) & (values <= 200))), axis=axis + ) + assert np.array_equal(out.data, wanted) + + def test_samples(self, random_patch): + """A sample range is complemented in samples too.""" + out = random_patch.unselect(distance=(..., 10), samples=True) + whole = random_patch.get_array("distance") + assert np.array_equal(out.get_array("distance"), whole[10:]) + + def test_relative(self, random_patch): + """Relative selectors mean what they mean in select.""" + kept = random_patch.select(time=(1, None), relative=True) + dropped = random_patch.unselect(time=(1, None), relative=True) + assert len(kept.get_array("time")) + len(dropped.get_array("time")) == len( + random_patch.get_array("time") + ) + + def test_unselecting_everything_empties_the_dimension(self, random_patch): + """Removing the whole span is legal, and says so with a shape.""" + coord = random_patch.get_coord("distance") + out = random_patch.unselect(distance=(coord.min(), coord.max())) + assert out.shape[random_patch.dims.index("distance")] == 0 + + def test_each_named_coordinate_is_complemented(self, random_patch): + """ + Two names remove two ranges rather than the one intersection. + + The complement of a block is an L, which no array can hold, so + unselect removes the part of it which is expressible. + """ + time = random_patch.get_array("time") + window = (time[0], time[4]) + out = random_patch.unselect(distance=(50, 60), time=window) + assert len(out.get_array("distance")) == len( + random_patch.get_array("distance") + ) - len(random_patch.select(distance=(50, 60)).get_array("distance")) + assert len(out.get_array("time")) == len(time) - 5 + + def test_unknown_coordinate_raises(self, random_patch): + """A misspelled name is the error it is in select.""" + with pytest.raises(PatchCoordinateError, match="not found in patch"): + random_patch.unselect(not_a_coord=(1, 2)) + + def test_non_dimensional_coordinate(self, random_patch): + """A coordinate along a dimension trims that dimension.""" + size = random_patch.coord_shapes["distance"][0] + patch = random_patch.update_coords(quality=("distance", np.arange(size))) + out = patch.unselect(quality=(0, 9)) + assert len(out.get_array("quality")) == size - 10 + + class TestOrder: """Tests for ordering Patches.""" From 00a137c94dce8856b669cf91c799e755575ad94c Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 21:33:51 +0200 Subject: [PATCH 3/8] Select and unselect on the coordinates along the fiber A track, an annotation group, or a geometry axis describes channels rather than whole patches, so selecting on one trims each patch to the channels which match and subdivides it where the matching region is disjoint. len can grow, which is why this builds a plan rather than filtering rows. None of the semantics are new. The values are projected onto the channels by the function Patch.enrich projects them with, and judged by the predicate the index applies to a stated attr, so globs, sequences, ranges and booleans mean here exactly what they mean there -- and selection cannot disagree with enrichment about which channel belongs to what. Only None needed its own reading: it is how a query spells the undefined marker _fill_from_intervals writes. The mask runs along one dimension, so unlike a patch's rectangle its complement is exact, and unselect takes it. Naming attrs as well stays one complement rather than two: a patch the attrs never matched keeps every channel, since the selection never held it. Only a real dimension can be trimmed, since what a non-dimensional coordinate says about the one it runs along is not in the index, and a name the index already uses keeps its own meaning -- distance is the patch's axis whether or not the inventory could also place it. --- dascore/core/spool.py | 194 +++++++++++++++++----- dascore/proc/inventory.py | 215 +++++++++++++++++++++++++ tests/test_core/test_spool.py | 6 +- tests/test_proc/test_proc_inventory.py | 55 ++++--- 4 files changed, 411 insertions(+), 59 deletions(-) diff --git a/dascore/core/spool.py b/dascore/core/spool.py index 35a44f48..ca538c1c 100644 --- a/dascore/core/spool.py +++ b/dascore/core/spool.py @@ -234,6 +234,25 @@ def _requested_names(_attrs, _coords, kwargs) -> set[str]: return out +def _without_names(spec: namespace_select_type, names) -> namespace_select_type: + """ + Return an `_attrs`/`_coords` argument with some names taken out. + + A channel-level name is answered by the inventory rather than the + index, so it has to leave the spec before the index sees it — and the + tag form as well as the mapping one, since a tag left behind + designates a bare keyword which went with it. + """ + if not names or spec is None: + return spec + if isinstance(spec, Mapping): + return {name: value for name, value in spec.items() if name not in names} + if isinstance(spec, str): + return None if spec in names else spec + kept = [x for x in spec if x not in names] + return kept or None + + def _namespace_names(spec: namespace_select_type) -> set[str]: """ Return the names an `_attrs`/`_coords` argument designates. @@ -921,8 +940,8 @@ def select( **kwargs, ) -> Self: """{doc}.""" - inventory_query, _attrs, _coords, kwargs = self._split_inventory_query( - _attrs, _coords, kwargs, samples + attr_query, channel_query, _attrs, _coords, kwargs = ( + self._split_inventory_query(_attrs, _coords, kwargs, samples) ) catalog = self._catalog.select( _attrs=_attrs, @@ -932,8 +951,10 @@ def select( **kwargs, ) out = self._new_from_catalog(catalog) - if inventory_query: - out = out._select_from_inventory(inventory_query) + if attr_query: + out = out._select_from_inventory(attr_query) + if channel_query: + out = out._select_channels(channel_query) return out @compose_docstring(doc=get_docstring(BaseSpool.unselect)) @@ -947,14 +968,11 @@ def unselect( """{doc}.""" requested = _requested_names(_attrs, _coords, kwargs) known_attrs, known_coords = self._index_names() - selectable = set() + selectable, channels = set(), {} if self._inventory is not None: names = self._inventory.get_names() - self._check_channel_level( - requested, - set(names.coords), - known_attrs | known_coords | set(names.attrs), - _namespace_names(_coords), + channels = self._channel_query( + requested, names, known_attrs, known_coords, _coords, kwargs ) selectable = set(names.attrs) - known_coords attrs, coords = resolve_selector_namespaces( @@ -962,21 +980,23 @@ def unselect( known_coords, _attrs=_attrs, _coords=_coords, - kwargs=kwargs, + kwargs={k: v for k, v in kwargs.items() if k not in channels}, ) if coords: msg = ( - f"{sorted(coords)} name coordinates, which unselect cannot " - "take yet: selecting on one trims each patch, so removing " - "a range means cutting every patch into the pieces outside " - "it rather than choosing between patches. Select the ranges " - "to keep instead." + f"{sorted(coords)} name coordinates of the patches " + "themselves, which unselect cannot take: removing a range " + "means cutting every patch into the pieces outside it " + "rather than choosing between patches. Select the ranges " + "to keep instead, or use Patch.unselect on each patch. The " + "coordinates an inventory defines along the fiber are " + "different -- removing one of those chooses channels." ) raise InvalidSpoolQueryError(msg) # A no-op selector selects everything, so its complement is an # empty spool -- and "remove nothing" reads just as naturally. stated = {k: v for k, v in attrs.items() if v is not None and v is not Ellipsis} - if not stated: + if not stated and not channels: msg = ( "unselect needs something to remove; " f"{sorted(requested) or 'nothing'} names no selection. " @@ -986,31 +1006,58 @@ def unselect( raise ParameterError(msg) # The complement is taken against select itself rather than by # negating each predicate, so the two can never drift apart. - removed = self.select(_attrs=stated)._catalog._ordered_ids() - return self._restrict_to_rows(removed, keep=False) + if not channels: + removed = self.select(_attrs=stated)._catalog._ordered_ids() + return self._restrict_to_rows(removed, keep=False) + # With both, the complement is still one set: a patch keeps every + # channel unless the attrs matched it, and the channels the fiber + # query did not match when they did. Complementing the two halves + # apart would drop a patch the whole selection never held. + matched = self if not stated else self.select(_attrs=stated) + return self._select_channels( + channels, complement=True, applies_to=matched._catalog._ordered_ids() + ) def _index_names(self) -> tuple[set[str], set[str]]: """The attr and coord names the index knows, read once per call.""" backend = self._catalog.backend return set(backend.attr_names()), set(backend.coord_names()) - def _check_channel_level(self, requested, coords, known, wanted_coords) -> None: + def _channel_query( + self, requested, names, known_attrs, known_coords, _coords, kwargs + ) -> dict: """ - Raise for a name the inventory defines along the fiber. + Return the selectors naming coordinates the inventory runs along + the fiber. A name which is also an attr resolves to the attr, as bare names always do, so only one the caller put in `_coords` is read as the coordinate — an annotation group may share an acquisition field's - name, and selecting on the field must keep working. - """ - candidates = (requested - known) | (wanted_coords & coords) - if channel_level := sorted(candidates & coords): + name, and selecting on the field must keep working. A name the + index already uses for a coordinate keeps that meaning outright: + `distance` is the patch's own axis whether or not an inventory + could also place it on the fiber, and an inventory must not + quietly move a name out of the namespace it has always been in. + """ + coords = set(names.coords) - known_coords + known = known_attrs | known_coords | set(names.attrs) + candidates = (requested - known) | (_namespace_names(_coords) & coords) + wanted = candidates & coords + if not wanted: + return {} + # The tag form of `_coords` designates bare kwargs, so the selector + # is among them; only the mapping form carries one itself. + mapped = _coords if isinstance(_coords, Mapping) else {} + stated = {**mapped, **kwargs} + if named := sorted(wanted - set(stated)): msg = ( - f"{channel_level} name coordinates the attached inventory " - "defines along the fiber, which selection cannot trim to " - "yet. Enrich the patches and select on each one instead." + f"{named} name coordinates the attached inventory defines " + "along the fiber. They are selected as ordinary keywords or " + "through _coords, not through _attrs, since they describe " + "channels rather than whole patches." ) raise InvalidSpoolQueryError(msg) + return {name: stated[name] for name in wanted} def _split_inventory_query(self, _attrs, _coords, kwargs, samples): """ @@ -1021,16 +1068,23 @@ def _split_inventory_query(self, _attrs, _coords, kwargs, samples): and the inventory only fills in the others. """ if self._inventory is None: - return {}, _attrs, _coords, kwargs + return {}, {}, _attrs, _coords, kwargs requested = _requested_names(_attrs, _coords, kwargs) known_attrs, known_coords = self._index_names() names = self._inventory.get_names() - self._check_channel_level( - requested, - set(names.coords), - known_attrs | known_coords | set(names.attrs), - _namespace_names(_coords), + channels = self._channel_query( + requested, names, known_attrs, known_coords, _coords, kwargs ) + if channels and samples: + msg = ( + f"{sorted(channels)} name coordinates the inventory defines " + "along the fiber, which have no sample numbering of their " + "own: the channels they describe are the patch's, and " + "samples=True asks about that axis instead." + ) + raise InvalidSpoolQueryError(msg) + kwargs = {k: v for k, v in kwargs.items() if k not in channels} + _coords = _without_names(_coords, channels) # A name the index already uses for a coordinate keeps its meaning; # bare names resolve to attrs first, and an inventory must not # quietly move one out of the namespace it has always been in. @@ -1038,7 +1092,7 @@ def _split_inventory_query(self, _attrs, _coords, kwargs, samples): # samples=True selections are coordinate-only, so an attr among # them is an error the index states better than this can. if samples or not requested & selectable: - return {}, _attrs, _coords, kwargs + return {}, channels, _attrs, _coords, kwargs attrs, coords = resolve_selector_namespaces( known_attrs | selectable, known_coords, @@ -1054,7 +1108,75 @@ def _split_inventory_query(self, _attrs, _coords, kwargs, samples): # A bare None or ... selects everything, here as everywhere. if value is not None and value is not Ellipsis: query[name] = value - return query, attrs, coords, {} + return query, channels, attrs, coords, {} + + def _select_channels(self, query: dict, *, complement=False, applies_to=None): + """ + Trim each patch to the channels the inventory says match. + + Where the matching region is disjoint along the fiber — a track + which passes in and out of the selection, or an uncovered zone in + the middle of a path — the patch is subdivided so each contiguous + run becomes its own patch. Selection therefore changes both the + shape and the number of patches, which is why it builds a plan + rather than filtering rows. + + Parameters + ---------- + query + The channel-level selectors, by inventory name. + complement + Keep the channels the query does *not* match. The mask is one + dimensional, so unlike a patch's rectangle its complement is + exactly expressible. + applies_to + The rows the query judges; any other row keeps every channel. + `unselect` uses it to leave a patch its attrs never matched + whole, which is what makes the two halves one complement. + """ + from dascore.proc.inventory import ( # noqa: PLC0415 + resolve_channel_pieces, + resolve_contexts, + ) + + source_rows, working = self._plan_frames() + if not len(working): + return self + columns = _resolution_columns(working) + contexts = ( + np.full(len(working), None, dtype=object) + if columns is None + else resolve_contexts(self._inventory, *columns) + ) + if applies_to is not None: + judged = np.isin(working["_patch_id"].to_numpy(), np.asarray(applies_to)) + contexts = np.where(judged, contexts, None) + name, pieces, reasons = resolve_channel_pieces( + self._inventory, contexts, working, query, complement=complement + ) + _refuse_rows( + source_rows, + reasons, + "are described by the inventory but cannot have their channels " + "placed along the fiber", + ) + if name is None: + # No row has a fiber to be judged along, so the query matched + # nothing: an empty spool, or the whole of it complemented. + return self if complement else self._restrict_to_rows([]) + bounds = list(zip(working[f"{name}_min"], working[f"{name}_max"], strict=True)) + whole = [ + len(row) == 1 and tuple(row[0]) == pair + for row, pair in zip(pieces, bounds, strict=True) + ] + if all(whole): # every patch kept entire: nothing to plan + return self + if all(keep or not row for keep, row in zip(whole, pieces, strict=True)): + # Every patch is kept whole or dropped, so this is a filter and + # the relation it presents need not be rebuilt. + kept = working["_patch_id"].to_numpy()[[bool(x) for x in pieces]] + return self._restrict_to_rows(kept) + return self._subdivided(source_rows, working, pieces, name) def _select_from_inventory(self, query: dict) -> Self: """ diff --git a/dascore/proc/inventory.py b/dascore/proc/inventory.py index 2cd8ccae..be9a9b3d 100644 --- a/dascore/proc/inventory.py +++ b/dascore/proc/inventory.py @@ -30,6 +30,7 @@ ) from dascore.exceptions import ( InvalidInventoryError, + InvalidSpoolQueryError, ParameterError, PatchError, UnresolvedPatchError, @@ -770,6 +771,220 @@ def _get_coord_values(inventory, path, name, distances): return _get_annotation_coord(path, name, distances) +# --- channel selection over index rows -------------------------------- + + +def _channel_placement(dims: set[str], acquisition) -> tuple: + """ + Return the dimension a row's channels are placed along, and its axis. + + The patch-level twin, `_get_channel_axes`, reads the map's axes off + whichever of the patch's coordinates state them, dimensional or not. + Here only a dimension will do: a plan trims a dimension, and what a + non-dimensional coordinate says about the one it runs along is not in + the index. Refusing rather than guessing is what keeps selection from + answering differently than enrichment would. + + Returns + ------- + A `(name, axis)` pair, or `(None, reason)` naming why there is none. + """ + dist_map = acquisition.distance_map + if dist_map is None: + return None, ( + f"{acquisition.code!r} defines no distance_map, so its channels " + "cannot be placed on the optical path" + ) + found = {} + for axis in dist_map.axes: + # One coordinate per axis, the same preference `_get_channel_axes` + # applies; two axes landing on one dimension is that dimension + # stating both, which the map's own validator has already agreed. + for name in _AXIS_COORDS[axis]: + if name in dims: + found.setdefault(name, axis) + break + if not found: + wanted = sorted({x for axis in dist_map.axes for x in _AXIS_COORDS[axis]}) + return None, ( + f"has dimensions {sorted(dims)}, and {acquisition.code!r} places " + f"channels by one of {wanted}" + ) + if len(found) > 1: + return None, ( + f"carries {sorted(found)} as separate dimensions, so which of " + f"them {acquisition.code!r} places its channels by is ambiguous" + ) + name, axis = next(iter(found.items())) + return name, axis + + +def _undefined_mask(values) -> np.ndarray: + """ + Return the channels an interval track states nothing about. + + `None` is how a query spells absence, and `_fill_from_intervals` + decides how absence is stored: a string array has no null, so it is + the empty string there, NaN in a numeric one, and a membership group + is simply False where nothing includes it. The two must agree, which + is why this reads as that function's mirror. + """ + array = np.asarray(values) + if array.dtype == bool: + return ~array + if np.issubdtype(array.dtype, np.number): + return np.isnan(array) + return array == "" + + +def _mask_pieces(mask: np.ndarray, grid: np.ndarray, low, high) -> list[tuple]: + """ + Return the inclusive envelope of each run of kept channels. + + The row's own bounds are used at its ends rather than the grid's + rebuilt ones. They are the same channel, but a float grid can land an + ulp off, and a piece which does not compare equal to the row it + covers would be trimmed on load instead of passing through. + """ + edges = np.flatnonzero(np.diff(np.concatenate([[0], mask.view(np.int8), [0]]))) + last = len(grid) - 1 + return [ + (low if start == 0 else grid[start], high if stop - 1 == last else grid[stop - 1]) + for start, stop in zip(edges[::2], edges[1::2], strict=True) + ] + + +def resolve_channel_pieces( + inventory, contexts, frame, query, *, complement: bool = False +) -> tuple: + """ + Return the channel dimension, each row's kept pieces, and refusals. + + The pieces are what a selection along the fiber keeps: one per + contiguous run of matching channels, so a query a path answers in two + places subdivides the row into two. A row with no context keeps + nothing and says nothing about it — an inventory-backed selector + silently does not match a patch the inventory is silent about, just + as a patch lacking an attr is not selected on it. + + Every value is projected onto the channels the way `Patch.enrich` + projects it and judged by the predicate the index applies to a stated + attr, so a selector cannot mean one thing here and another in either. + + Parameters + ---------- + inventory + The inventory to resolve against. + contexts + Each row's resolved context, or None where it has none. + frame + The relation being selected; one row per patch. + query + The channel-level selectors, by inventory name. + complement + Keep the channels the query does *not* match. The mask runs along + one dimension, so unlike a patch's rectangle its complement is + exactly expressible; a row with no context keeps everything, + being a row the selection never held. + + Returns + ------- + A `(name, pieces, reasons)` triple. `reasons` holds a refusal per row + and None elsewhere; when any row is refused the other two are None, + since the caller raises rather than selecting. + """ + placements, reasons = [], [] + for context, dims in zip(contexts, frame["dims"], strict=True): + placement = (None, None) + if context is not None: + placement = _channel_placement(set(dims.split(",")), context.acquisition) + placements.append(placement) + # The second half of a placement is the axis when there is one and + # the reason there is not, so only a nameless one carries a reason. + reasons.append(placement[1] if context is not None and not placement[0] else None) + named = {name for name, _ in placements if name is not None} + if len(named) > 1: + joined = ", ".join(sorted(named)) + msg = ( + "The patches of this spool place their channels along different " + f"dimensions ({joined}), so a selection along the fiber has no " + "one dimension to trim. Select them apart first." + ) + raise InvalidSpoolQueryError(msg) + name = next(iter(named), None) + if any(x is not None for x in reasons) or name is None: + # No dimension to trim along, so nothing was judged and there are + # no pieces to report; the caller reads that off `name`. + return name, None, reasons + pieces, unusable = _row_pieces( + inventory, contexts, placements, frame, name, query, complement + ) + return name, pieces, [x or y for x, y in zip(reasons, unusable, strict=True)] + + +def _row_bounds(frame, name: str): + """Return each row's inclusive envelope along one dimension.""" + return zip( + frame[f"{name}_min"].to_numpy(), frame[f"{name}_max"].to_numpy(), strict=True + ) + + +def _row_pieces( + inventory, contexts, placements, frame, name, query, complement +) -> tuple[list[list], list]: + """Return the kept pieces of each row, and any row with no usable grid.""" + steps = frame[f"{name}_step"].to_numpy() + # Keyed by identity because that is what is cheap: sibling epochs + # resolve to one shared context object, and `__eq__` on an inventory + # model dumps the whole subtree. A miss only recomputes. + cache: dict[tuple, list] = {} + out, unusable = [], [] + rows = zip(contexts, placements, _row_bounds(frame, name), steps, strict=True) + for context, (dim, axis), (low, high), step in rows: + unusable.append(None) + if context is None or dim is None: + # A row the inventory is silent about is not selected, and so + # its complement keeps it whole. + out.append([(low, high)] if complement else []) + continue + if pd.isnull(step) or not step: + # Which channels match is decided on the sample grid, and the + # step is its only description. Guessing would trim the wrong + # channels silently, which is worse than saying so. + unusable[-1] = "states no channel spacing to place its channels on" + out.append([]) + continue + key = (id(context), axis, low, high, step, complement) + if (pieces := cache.get(key)) is None: + grid = np.arange(round((high - low) / step) + 1) * step + low + distances = context.acquisition.channel_to_distance(grid, axis=axis) + mask = np.ones(len(grid), dtype=bool) + for wanted, selector in query.items(): + mask &= _channel_matches( + inventory, context.optical_path, wanted, selector, distances + ) + pieces = _mask_pieces(~mask if complement else mask, grid, low, high) + cache[key] = pieces + out.append(pieces) + return out, unusable + + +def _channel_matches(inventory, path, name, selector, distances) -> np.ndarray: + """Return the channels one selector matches, as the index would.""" + from dascore.io.index.query import evaluate_attr_predicate # noqa: PLC0415 + + values = _get_coord_values(inventory, path, name, distances) + if values is None: + # A name this path defines nowhere states nothing about any of its + # channels, so it matches none of them -- listing a name is not + # promising a value for it. + return np.zeros(len(distances), dtype=bool) + values = values.values if isinstance(values, BaseCoord) else values + if selector is None: + return _undefined_mask(values) + return evaluate_attr_predicate(list(values), name, selector) + + def _coords_equal(existing, values) -> bool: """Return True when a patch coordinate already holds these values.""" other = values if isinstance(values, BaseCoord) else get_coord(data=values) diff --git a/tests/test_core/test_spool.py b/tests/test_core/test_spool.py index c05503ca..1f857789 100644 --- a/tests/test_core/test_spool.py +++ b/tests/test_core/test_spool.py @@ -591,9 +591,9 @@ def test_coordinates_raise(self, diverse_spool): cuts every patch into the pieces outside it -- one patch becoming two -- rather than choosing between patches. """ - with pytest.raises(InvalidSpoolQueryError, match="cannot take yet"): + with pytest.raises(InvalidSpoolQueryError, match="unselect cannot take"): diverse_spool.unselect(time=("2020-01-03", None)) - with pytest.raises(InvalidSpoolQueryError, match="cannot take yet"): + with pytest.raises(InvalidSpoolQueryError, match="unselect cannot take"): diverse_spool.unselect(_coords={"time": ("2020-01-03", None)}) def test_unknown_name_raises(self, diverse_spool): @@ -653,7 +653,7 @@ def test_naming_only_a_no_op_raises(self, diverse_spool): def test_coords_tag_form_raises(self, diverse_spool): """The tag form names bare kwargs, and is refused the same way.""" - with pytest.raises(InvalidSpoolQueryError, match="cannot take yet"): + with pytest.raises(InvalidSpoolQueryError, match="unselect cannot take"): diverse_spool.unselect(_coords="time", time=("2020-01-03", None)) def test_base_spool_unselect_raises(self, random_spool): diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index 71244465..2c92c0da 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -722,11 +722,11 @@ def test_policy_survives_a_union(self, patch, inventory): class TestInventoryQueryError: """An attached inventory changes what an unknown query name means.""" - def test_inventory_field_names_the_inventory(self, patch, inventory): + def test_inventory_field_is_answered_by_the_inventory(self, patch, inventory): """The index's 'no such attribute' would deny a field which exists.""" spool = dc.spool(patch).attach_inventory(inventory) - with pytest.raises(InvalidSpoolQueryError, match="along the fiber"): - spool.select(coupling="cement") + # No channel is coupled in cement, so the field answers with none. + assert len(spool.select(coupling="cement")) == 0 def test_plain_spool_keeps_its_message(self, patch): """With no inventory there is nothing to add.""" @@ -1456,18 +1456,24 @@ def test_original_is_unchanged(self, two_patch_spool, inventory): spool.select(gauge_length=99.0) assert len(spool) == 2 - def test_channel_level_name_says_so(self, two_patch_spool, inventory): + def test_channel_level_name_trims_channels(self, two_patch_spool, inventory): """ - A track name is a real inventory name which select cannot use yet. + A track name selects channels rather than whole patches. - Saying which of the two it is takes the names accessor, and this - is what it was added for. + Telling one from an attr takes the names accessor, and this is + what it was added for. """ spool = two_patch_spool.attach_inventory(inventory) - with pytest.raises(InvalidSpoolQueryError, match="along the fiber"): - spool.select(coupling="trench") - with pytest.raises(InvalidSpoolQueryError, match="along the fiber"): - spool.select(zone="east") + out = spool.select(coupling="trench") + assert len(out) == 2 + assert set(out.get_contents()["distance_max"]) == {150.0} + + def test_a_channel_level_name_nothing_matches_keeps_nothing( + self, two_patch_spool, inventory + ): + """A value no channel holds selects no channel, so no patch.""" + spool = two_patch_spool.attach_inventory(inventory) + assert len(spool.select(zone="east")) == 0 def test_unknown_name_still_raises(self, two_patch_spool, inventory): """A name neither side knows is a misspelling, and says so.""" @@ -1556,11 +1562,18 @@ def test_stated_value_wins(self, patch, inventory): "random" ] - def test_channel_level_name_says_so(self, two_patch_spool, inventory): - """A track name is as unsupported here as it is in select.""" + def test_channel_level_name_removes_channels(self, two_patch_spool, inventory): + """ + A track name removes channels rather than whole patches. + + Which is the one thing `unselect` refuses a real coordinate for: + a coordinate range would leave a hole in every patch, while a + channel query chooses which channels a patch holds. + """ spool = two_patch_spool.attach_inventory(inventory) - with pytest.raises(InvalidSpoolQueryError, match="along the fiber"): - spool.unselect(coupling="trench") + out = spool.unselect(coupling="trench") + assert len(out) == 2 + assert set(out.get_contents()["distance_min"]) == {151.0} class TestInventorySelectCoverage: @@ -1858,7 +1871,7 @@ def test_a_group_may_share_an_attrs_name(self, patch, inventory): Bare names resolve to attrs first, so selecting on the field has to keep working; only a caller who asked for `_coords` means the - group, which is the channel-level half and not supported yet. + group, which trims channels rather than choosing patches. """ path = inventory.networks[0].fiber_arrays[0].optical_paths[0] clash = inventory.replace( @@ -1880,11 +1893,13 @@ def test_a_group_may_share_an_attrs_name(self, patch, inventory): spool = dc.spool([patch]).attach_inventory(clash) assert len(spool.select(gauge_length=10.0)) == 1 assert len(spool.select(_attrs={"gauge_length": 10.0})) == 1 - for form in ({"gauge_length": 10.0}, "gauge_length", ["gauge_length"]): + # The group covers the first metre of the path, which is the + # patch's first channel and nothing else. + for form in ({"gauge_length": "odd"}, "gauge_length", ["gauge_length"]): # The mapping form and both tag forms all name the coordinate. - kwargs = {} if isinstance(form, Mapping) else {"gauge_length": 10.0} - with pytest.raises(InvalidSpoolQueryError, match="along the fiber"): - spool.select(_coords=form, **kwargs) + kwargs = {} if isinstance(form, Mapping) else {"gauge_length": "odd"} + out = spool.select(_coords=form, **kwargs) + assert out.get_contents()["distance_max"].tolist() == [0.0] def test_a_field_scalar_on_another_path_is_kept(self): """ From 4c1cb2231ebe3967742e3c39557244551eadadac Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 21:42:32 +0200 Subject: [PATCH 4/8] Add Spool.split_by, one patch per value along the fiber It is the same subdivision channel selection performs, asked a different way: instead of one query deciding which channels to keep, each value a group takes decides an output of its own. Groups may overlap, so a channel can land in more than one output and the pieces of a row need not be disjoint. Every kind splits -- strings by value, a membership group into the channels it includes and those it does not, a numeric one by each distinct measurement -- and include/exclude are globs over the value written as a string, so one vocabulary covers all three. The value is stamped on each output so overlapping siblings stay apart and later operations can select on it. That needed the plan resolver to know an output may state attrs of its own rather than inheriting its members': assembling one says nothing about why it was cut out. Naming the columns rather than serializing them keeps each value's own type, which a numeric group needs. --- dascore/core/spool.py | 134 +++++++++++++++++++++++++- dascore/io/index/planned.py | 29 +++++- dascore/proc/inventory.py | 185 +++++++++++++++++++++++++++--------- 3 files changed, 296 insertions(+), 52 deletions(-) diff --git a/dascore/core/spool.py b/dascore/core/spool.py index ca538c1c..67b0224c 100644 --- a/dascore/core/spool.py +++ b/dascore/core/spool.py @@ -6,6 +6,7 @@ import inspect import warnings from collections.abc import Callable, Generator, Iterator, Mapping, Sequence +from dataclasses import replace from functools import singledispatch from pathlib import Path from typing import TYPE_CHECKING, ClassVar, Literal, TypeVar, overload @@ -56,6 +57,7 @@ from dascore.utils.misc import ( _spool_map, deep_equality_check, + iterate, ) from dascore.utils.namespace import NamespaceOwner from dascore.utils.patch import ( @@ -234,6 +236,35 @@ def _requested_names(_attrs, _coords, kwargs) -> set[str]: return out +def _glob_filter(include, exclude): + """ + Return a predicate deciding which split values to keep. + + The patterns are matched against each value written as a string, so + one vocabulary covers every kind a group can hold: `"hole_*"` reads a + categorical group, and a membership group is `"True"` and `"False"`. + Globs mean what they mean everywhere else here, which is what SQLite + means by them rather than what `fnmatch` does. + """ + from dascore.io.index.query import glob_to_regex # noqa: PLC0415 + + patterns = tuple( + None if spec is None else [glob_to_regex(str(x)) for x in iterate(spec)] + for spec in (include, exclude) + ) + + def keep(value) -> bool: + wanted, unwanted = patterns + text = str(value) + if unwanted is not None and any(x.match(text) for x in unwanted): + # Excluding wins, so naming a family and carving one out of it + # reads in either order. + return False + return wanted is None or any(x.match(text) for x in wanted) + + return keep + + def _without_names(spec: namespace_select_type, names) -> namespace_select_type: """ Return an `_attrs`/`_coords` argument with some names taken out. @@ -1421,6 +1452,95 @@ def enrich( new._on_unresolved = on_unresolved return new + def split_by( + self, + name: str, + *, + include: str | Sequence[str] | None = None, + exclude: str | Sequence[str] | None = None, + stamp: bool = True, + ) -> Self: + """ + Expand the spool into one patch per value of an inventory coordinate. + + Most often an annotation group. Every kind of group splits: a + categorical one by each of its strings, a membership group into + the channels it includes and those it does not, and a numeric one + by each distinct measurement. Groups are allowed to overlap, so a + channel may appear in more than one output patch, and a patch + whose channels take several values becomes several patches — this + can greatly expand the spool. + + Parameters + ---------- + name + The inventory-derived coordinate to split on. + include, exclude + Glob patterns matched against each value written as a string, + so `"trench_?a"` and `"hole_*"` both work and `True` is + matched as `"True"`. With `include`, only values matching one + of them are kept; `exclude` drops the values it matches, and + wins where both match. + stamp + Whether to record the value on each output patch as an attr + named after the coordinate, so overlapping siblings stay + distinguishable and later operations can select on it. Pass + False for a nested split, where the second should not + overwrite the first. + + Examples + -------- + >>> import dascore as dc + >>> from dascore.examples import inventory_patch_pair + >>> + >>> patch, inventory = inventory_patch_pair() + >>> spool = dc.spool(patch).attach_inventory(inventory) + >>> + >>> # The example path annotates two zones along the fiber. + >>> zones = spool.split_by("zone") + >>> assert len(zones) == 2 + >>> assert set(zones.get_contents()["zone"]) == {"north", "south"} + >>> + >>> # Which can be narrowed by a glob over the values. + >>> assert len(spool.split_by("zone", include="nor*")) == 1 + """ + from dascore.proc.inventory import resolve_split_pieces # noqa: PLC0415 + + if self._inventory is None: + msg = ( + "Spool.split_by needs an inventory to split on: the values " + "it expands into are the ones an inventory states along the " + "fiber. Attach one with Spool.attach_inventory." + ) + raise ParameterError(msg) + source_rows, working = self._plan_frames() + contexts = self._plan_contexts(working) + dim, rows, reasons = resolve_split_pieces( + self._inventory, contexts, working, name, _glob_filter(include, exclude) + ) + _refuse_rows( + source_rows, + reasons, + "are described by the inventory but cannot have their channels " + "placed along the fiber", + ) + if dim is None: # nothing to split: no row has a fiber to split on + return self._restrict_to_rows([]) + pieces = [[piece for _, piece in row] for row in rows] + marks = None + if stamp: + marks = (name, [value for row in rows for value, _ in row]) + return self._subdivided(source_rows, working, pieces, dim, marks) + + def _plan_contexts(self, working) -> np.ndarray: + """Resolve each row of a planning frame to its inventory context.""" + from dascore.proc.inventory import resolve_contexts # noqa: PLC0415 + + columns = _resolution_columns(working) + if columns is None: + return np.full(len(working), None, dtype=object) + return resolve_contexts(self._inventory, *columns) + def conform_to_inventory( self, inventory=None, @@ -1534,24 +1654,32 @@ def conform_to_inventory( sources, kept, subdivision_pieces(kept, cuts, "time"), "time" ) - def _subdivided(self, sources, rows, pieces, name: str) -> Self: + def _subdivided(self, sources, rows, pieces, name: str, stamp=None) -> Self: """ Return the spool whose patches are the given pieces of these rows. The pieces are a plan rather than a rewritten relation, so the outputs *are* the contents rows — `len` and `get_contents` stay exact — and loading goes through the machinery which already - trims a member on extraction. + trims a member on extraction. `stamp` names an attr to record on + each output, in the order the pieces were given. """ from dascore.io.index.planned import derived_catalog # noqa: PLC0415 + plan = build_subdivision_plan(rows, pieces, name) + stamped = () + if stamp is not None: + stamp_name, values = stamp + plan = replace(plan, outputs=plan.outputs.assign(**{stamp_name: values})) + stamped = (stamp_name,) catalog = derived_catalog( source_rows=sources, - plan=build_subdivision_plan(rows, pieces, name), + plan=plan, parent=self._catalog, merge_kwargs={}, mode="chunk", origin_path=self.spool_path, + stamped=stamped, ) return self._new_from_catalog(catalog) diff --git a/dascore/io/index/planned.py b/dascore/io/index/planned.py index 52ebfb74..a27b8715 100644 --- a/dascore/io/index/planned.py +++ b/dascore/io/index/planned.py @@ -367,6 +367,7 @@ def __init__( mode: str = "chunk", check_behavior: WARN_LEVELS = "warn", origin_path=None, + stamped: tuple[str, ...] = (), ): if "output_id" not in member_rows.columns: msg = "member_rows must carry an output_id column." @@ -382,6 +383,8 @@ def __init__( self.check_behavior = check_behavior # informational only: the directory/file the plan derived from self.origin_path = origin_path + # attrs the outputs state about themselves rather than inherit + self.stamped = tuple(stamped) def live_entries(self) -> dict[str, dc.Patch]: """Expose the loader's live registry (for absorption/transfer).""" @@ -460,7 +463,7 @@ def resolve(self, row: Mapping, **trim) -> dc.Patch: if self.mode == "identity": # one untouched member per output; residuals apply at load assert len(members) == 1 - return self._load_member(members.iloc[0].to_dict()) + return self._stamp(self._load_member(members.iloc[0].to_dict()), row) if self.mode == "concat": patches = [ self._load_member(kwargs) for kwargs in members.to_dict("records") @@ -469,11 +472,29 @@ def resolve(self, row: Mapping, **trim) -> dc.Patch: patches, check_behavior=self.check_behavior, **{self.dim: None} ) assert len(out) == 1 - return out[0] + return self._stamp(out[0], row) joined = members.assign(current_index=output_id) patches = self._assembler()._patch_from_instruction_df(joined) assert len(patches) == 1 - return patches[0] + return self._stamp(patches[0], row) + + def _stamp(self, patch: dc.Patch, row: Mapping) -> dc.Patch: + """ + Apply the attrs the outputs state about themselves, if any. + + An output is assembled from its members, so it carries their + attrs and knows nothing of why it was cut out. `stamped` is how + an operation which does know says so -- `Spool.split_by` + recording which value each patch was split on -- and it keeps + the patch which comes out agreeing with the row `get_contents` + shows for it. + """ + if not self.stamped: + return patch + return patch.update_attrs(**{x: row[x] for x in self.stamped}) + + + def _residual_ranges(residuals) -> dict: @@ -506,6 +527,7 @@ def derived_catalog( mode: str = "chunk", check_behavior: WARN_LEVELS = "warn", origin_path=None, + stamped: tuple[str, ...] = (), ) -> PatchCatalog: """ Materialize a plan into a fresh in-memory catalog. @@ -574,6 +596,7 @@ def derived_catalog( mode=mode, check_behavior=check_behavior, origin_path=origin_path, + stamped=stamped, ) backend = get_backend(":memory:") coord_dims_map = {} if parent is None else parent.backend.coord_dims_map() diff --git a/dascore/proc/inventory.py b/dascore/proc/inventory.py index be9a9b3d..961cb956 100644 --- a/dascore/proc/inventory.py +++ b/dascore/proc/inventory.py @@ -893,6 +893,34 @@ def resolve_channel_pieces( and None elsewhere; when any row is refused the other two are None, since the caller raises rather than selecting. """ + name, placements, reasons = channel_placements(contexts, frame) + if any(x is not None for x in reasons) or name is None: + # No dimension to trim along, so nothing was judged and there are + # no pieces to report; the caller reads that off `name`. + return name, None, reasons + out, unusable = [], [] + for row in _placed_rows(contexts, placements, frame, name): + unusable.append(row.reason) + if row.grid is None: + # A row the inventory is silent about is not selected, and so + # its complement keeps it whole; one it cannot place is refused + # above rather than answered for. + out.append([(row.low, row.high)] if complement and not row.reason else []) + continue + mask = np.ones(len(row.grid), dtype=bool) + for wanted, selector in query.items(): + mask &= _channel_matches( + inventory, row.context.optical_path, wanted, selector, row.distances + ) + out.append(_mask_pieces(~mask if complement else mask, *row.grid_bounds)) + return name, out, [x or y for x, y in zip(reasons, unusable, strict=True)] + + +def channel_placements(contexts, frame) -> tuple: + """ + Return the one dimension a relation's channels run along, and per row + how each is placed on it, or why it cannot be. + """ placements, reasons = [], [] for context, dims in zip(contexts, frame["dims"], strict=True): placement = (None, None) @@ -902,71 +930,136 @@ def resolve_channel_pieces( # The second half of a placement is the axis when there is one and # the reason there is not, so only a nameless one carries a reason. reasons.append(placement[1] if context is not None and not placement[0] else None) - named = {name for name, _ in placements if name is not None} + named = {x for x, _ in placements if x is not None} if len(named) > 1: joined = ", ".join(sorted(named)) msg = ( "The patches of this spool place their channels along different " - f"dimensions ({joined}), so a selection along the fiber has no " - "one dimension to trim. Select them apart first." + f"dimensions ({joined}), so an operation along the fiber has no " + "one dimension to work on. Select them apart first." ) raise InvalidSpoolQueryError(msg) - name = next(iter(named), None) - if any(x is not None for x in reasons) or name is None: - # No dimension to trim along, so nothing was judged and there are - # no pieces to report; the caller reads that off `name`. - return name, None, reasons - pieces, unusable = _row_pieces( - inventory, contexts, placements, frame, name, query, complement - ) - return name, pieces, [x or y for x, y in zip(reasons, unusable, strict=True)] + return next(iter(named), None), placements, reasons -def _row_bounds(frame, name: str): - """Return each row's inclusive envelope along one dimension.""" - return zip( - frame[f"{name}_min"].to_numpy(), frame[f"{name}_max"].to_numpy(), strict=True - ) +class PlacedRow(NamedTuple): + """One index row's channels, placed on its optical path.""" + + context: Any + low: Any + high: Any + grid: np.ndarray | None + distances: np.ndarray | None + reason: str | None + + @property + def grid_bounds(self) -> tuple: + """The arguments `_mask_pieces` reads a mask against.""" + return self.grid, self.low, self.high -def _row_pieces( - inventory, contexts, placements, frame, name, query, complement -) -> tuple[list[list], list]: - """Return the kept pieces of each row, and any row with no usable grid.""" +def _placed_rows(contexts, placements, frame, name: str): + """ + Yield each row's channel grid and where it lands on the optical path. + + A row with no context, or no usable step to rebuild its grid with, + yields no grid; the two are different answers, so only the second + carries a reason to refuse it with. + """ steps = frame[f"{name}_step"].to_numpy() + bounds = zip(frame[f"{name}_min"], frame[f"{name}_max"], strict=True) # Keyed by identity because that is what is cheap: sibling epochs # resolve to one shared context object, and `__eq__` on an inventory # model dumps the whole subtree. A miss only recomputes. - cache: dict[tuple, list] = {} - out, unusable = [], [] - rows = zip(contexts, placements, _row_bounds(frame, name), steps, strict=True) - for context, (dim, axis), (low, high), step in rows: - unusable.append(None) + cache: dict[tuple, tuple] = {} + for context, (dim, axis), (low, high), step in zip( + contexts, placements, bounds, steps, strict=True + ): if context is None or dim is None: - # A row the inventory is silent about is not selected, and so - # its complement keeps it whole. - out.append([(low, high)] if complement else []) + yield PlacedRow(context, low, high, None, None, None) continue if pd.isnull(step) or not step: - # Which channels match is decided on the sample grid, and the - # step is its only description. Guessing would trim the wrong - # channels silently, which is worse than saying so. - unusable[-1] = "states no channel spacing to place its channels on" - out.append([]) + # Which channels are which is decided on the sample grid, and + # the step is its only description. Guessing would trim the + # wrong channels silently, which is worse than saying so. + reason = "states no channel spacing to place its channels on" + yield PlacedRow(context, low, high, None, None, reason) continue - key = (id(context), axis, low, high, step, complement) - if (pieces := cache.get(key)) is None: + key = (id(context), axis, low, high, step) + if (placed := cache.get(key)) is None: grid = np.arange(round((high - low) / step) + 1) * step + low - distances = context.acquisition.channel_to_distance(grid, axis=axis) - mask = np.ones(len(grid), dtype=bool) - for wanted, selector in query.items(): - mask &= _channel_matches( - inventory, context.optical_path, wanted, selector, distances - ) - pieces = _mask_pieces(~mask if complement else mask, grid, low, high) - cache[key] = pieces - out.append(pieces) - return out, unusable + placed = (grid, context.acquisition.channel_to_distance(grid, axis=axis)) + cache[key] = placed + yield PlacedRow(context, low, high, *placed, None) + + +def resolve_split_pieces(inventory, contexts, frame, name, keep) -> tuple: + """ + Return the channel dimension, each row's pieces by value, and refusals. + + Splitting expands the spool into one patch per value a group takes + along each row, so a row is answered with a list of `(value, piece)` + pairs rather than pieces alone. A group whose values overlap gives + some channels to more than one output, which is the whole point of a + membership group and why the pieces of a row need not be disjoint. + + Parameters + ---------- + inventory + The inventory to resolve against. + contexts + Each row's resolved context, or None where it has none. + frame + The relation being split; one row per patch. + name + The inventory-derived coordinate to split on. + keep + Decides which values to emit, by the value itself. + + Returns + ------- + A `(dim, rows, reasons)` triple, `rows` holding the `(value, piece)` + pairs of each row in value order. + """ + dim, placements, reasons = channel_placements(contexts, frame) + if any(x is not None for x in reasons) or dim is None: + return dim, None, reasons + out, unusable = [], [] + for row in _placed_rows(contexts, placements, frame, dim): + unusable.append(row.reason) + out.append([]) + if row.grid is None: + continue + path = row.context.optical_path + values = _get_coord_values(inventory, path, name, row.distances) + if values is None: + # A group this path defines nowhere puts none of its channels + # anywhere, so the row contributes no output at all. + continue + values = values.values if isinstance(values, BaseCoord) else values + for value in _split_values(values): + if not keep(value): + continue + pieces = _mask_pieces(np.asarray(values) == value, *row.grid_bounds) + out[-1].extend((value, piece) for piece in pieces) + return dim, out, [x or y for x, y in zip(reasons, unusable, strict=True)] + + +def _split_values(values) -> list: + """ + Return the distinct values a group takes, in a stable order. + + Every kind splits, one output per value: a categorical group by its + strings, a membership group into the channels it includes and those + it does not, and a numeric one by each distinct measurement. Sorted + so the outputs of a spool do not depend on which channel came first. + """ + array = np.asarray(values) + if np.issubdtype(array.dtype, np.floating): + # NaN is how a numeric group spells "no value here", and it is + # never equal to itself, so it could not name an output anyway. + array = array[~np.isnan(array)] + return sorted(set(array.tolist())) def _channel_matches(inventory, path, name, selector, distances) -> np.ndarray: From faca4222f8761780d0ad14efb8179d78d4a9877d Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 21:55:06 +0200 Subject: [PATCH 5/8] Test the fiber operations, and lint them Covers every added line: the trims themselves and the data behind them, disjoint matches, composition with a time split, selecting twice, the complement, and the rows a fiber query cannot answer for -- a lag-time patch, an unevenly sampled one, an acquisition with no map, and a patch carrying two channel axes at once. Four tests which pinned the "not supported yet" error now pin what happens instead. One of them found a real bug: a name the index already uses for a coordinate has to keep its own meaning, and _coords={"distance": ...} was being read as optical distance. Absence is not a value to split on, which the first draft got wrong for strings -- the empty string a string coordinate carries instead of a null was making an output of its own. A membership group is the exception, since False there says something about every channel. --- dascore/core/spool.py | 22 +- dascore/io/index/planned.py | 3 - dascore/proc/coords.py | 3 +- dascore/proc/inventory.py | 25 +- tests/test_proc/test_proc_inventory.py | 573 +++++++++++++++++++++++++ tests/test_utils/test_chunk.py | 2 +- 6 files changed, 599 insertions(+), 29 deletions(-) diff --git a/dascore/core/spool.py b/dascore/core/spool.py index 67b0224c..a6cab27e 100644 --- a/dascore/core/spool.py +++ b/dascore/core/spool.py @@ -203,9 +203,7 @@ def _unsubdividable(rows: pd.DataFrame, pieces, name: str) -> list: not for the data to be restructured. """ return [ - f"at {row_pieces[0]}" - if row_pieces and (pd.isnull(step) or not step) - else None + f"at {row_pieces[0]}" if row_pieces and (pd.isnull(step) or not step) else None for step, row_pieces in zip(rows[f"{name}_step"], pieces, strict=True) ] @@ -277,7 +275,7 @@ def _without_names(spec: namespace_select_type, names) -> namespace_select_type: if not names or spec is None: return spec if isinstance(spec, Mapping): - return {name: value for name, value in spec.items() if name not in names} + return {str(k): v for k, v in spec.items() if k not in names} if isinstance(spec, str): return None if spec in names else spec kept = [x for x in spec if x not in names] @@ -1165,23 +1163,17 @@ def _select_channels(self, query: dict, *, complement=False, applies_to=None): `unselect` uses it to leave a patch its attrs never matched whole, which is what makes the two halves one complement. """ - from dascore.proc.inventory import ( # noqa: PLC0415 - resolve_channel_pieces, - resolve_contexts, - ) + from dascore.proc.inventory import resolve_channel_pieces # noqa: PLC0415 source_rows, working = self._plan_frames() if not len(working): return self - columns = _resolution_columns(working) - contexts = ( - np.full(len(working), None, dtype=object) - if columns is None - else resolve_contexts(self._inventory, *columns) - ) + contexts = self._plan_contexts(working) if applies_to is not None: + # A row the attrs did not match is a row the selection never + # held, so it is left unjudged rather than judged and kept. judged = np.isin(working["_patch_id"].to_numpy(), np.asarray(applies_to)) - contexts = np.where(judged, contexts, None) + contexts[~judged] = None name, pieces, reasons = resolve_channel_pieces( self._inventory, contexts, working, query, complement=complement ) diff --git a/dascore/io/index/planned.py b/dascore/io/index/planned.py index a27b8715..d98dadfe 100644 --- a/dascore/io/index/planned.py +++ b/dascore/io/index/planned.py @@ -494,9 +494,6 @@ def _stamp(self, patch: dc.Patch, row: Mapping) -> dc.Patch: return patch.update_attrs(**{x: row[x] for x in self.stamped}) - - - def _residual_ranges(residuals) -> dict: """Envelope-applicable value ranges from a residual tuple. diff --git a/dascore/proc/coords.py b/dascore/proc/coords.py index 09f4fecd..e60de253 100644 --- a/dascore/proc/coords.py +++ b/dascore/proc/coords.py @@ -623,8 +623,7 @@ def unselect( invalid_list = sorted(invalid_coords) valid_list = sorted(patch.coords.coord_map) msg = ( - f"Coordinate(s) {invalid_list} not found in patch coordinates: " - f"{valid_list}" + f"Coordinate(s) {invalid_list} not found in patch coordinates: {valid_list}" ) raise PatchCoordinateError(msg) complements = {} diff --git a/dascore/proc/inventory.py b/dascore/proc/inventory.py index 961cb956..e61bf3a6 100644 --- a/dascore/proc/inventory.py +++ b/dascore/proc/inventory.py @@ -849,7 +849,10 @@ def _mask_pieces(mask: np.ndarray, grid: np.ndarray, low, high) -> list[tuple]: edges = np.flatnonzero(np.diff(np.concatenate([[0], mask.view(np.int8), [0]]))) last = len(grid) - 1 return [ - (low if start == 0 else grid[start], high if stop - 1 == last else grid[stop - 1]) + ( + low if start == 0 else grid[start], + high if stop - 1 == last else grid[stop - 1], + ) for start, stop in zip(edges[::2], edges[1::2], strict=True) ] @@ -929,7 +932,8 @@ def channel_placements(contexts, frame) -> tuple: placements.append(placement) # The second half of a placement is the axis when there is one and # the reason there is not, so only a nameless one carries a reason. - reasons.append(placement[1] if context is not None and not placement[0] else None) + placed = context is not None and not placement[0] + reasons.append(placement[1] if placed else None) named = {x for x, _ in placements if x is not None} if len(named) > 1: joined = ", ".join(sorted(named)) @@ -990,7 +994,8 @@ def _placed_rows(contexts, placements, frame, name: str): grid = np.arange(round((high - low) / step) + 1) * step + low placed = (grid, context.acquisition.channel_to_distance(grid, axis=axis)) cache[key] = placed - yield PlacedRow(context, low, high, *placed, None) + grid, distances = placed + yield PlacedRow(context, low, high, grid, distances, None) def resolve_split_pieces(inventory, contexts, frame, name, keep) -> tuple: @@ -1053,13 +1058,17 @@ def _split_values(values) -> list: strings, a membership group into the channels it includes and those it does not, and a numeric one by each distinct measurement. Sorted so the outputs of a spool do not depend on which channel came first. + + Absence is not a value, so the channels a group says nothing about + make no output of their own — with the one exception a membership + group is: `False` there means "not in this group", which is a + statement about every channel rather than the absence of one. """ array = np.asarray(values) - if np.issubdtype(array.dtype, np.floating): - # NaN is how a numeric group spells "no value here", and it is - # never equal to itself, so it could not name an output anyway. - array = array[~np.isnan(array)] - return sorted(set(array.tolist())) + if array.dtype == bool: + return sorted(set(array.tolist())) + keep = ~np.isnan(array) if np.issubdtype(array.dtype, np.number) else array != "" + return sorted(set(array[keep].tolist())) def _channel_matches(inventory, path, name, selector, distances) -> np.ndarray: diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index 2c92c0da..ff0a05df 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -2542,3 +2542,576 @@ def test_passing_an_inventory_clears_enrichment(self, patch, path_epochs): out = spool.conform_to_inventory(path_epochs) assert out._enrich_kwargs is None assert "gauge_length" not in dict(out[0].attrs) + + +@pytest.fixture(scope="module") +def two_zones(inventory): + """The example inventory, with a group covering two separate stretches.""" + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + return inventory.replace( + path, + path.new( + annotations=( + *path.annotations, + OpticalPathAnnotation( + start_distance=110.0, end_distance=150.0, group="hole", value="a" + ), + OpticalPathAnnotation( + start_distance=300.0, end_distance=340.0, group="hole", value="a" + ), + ) + ), + ) + + +def _channels(spool): + """Every channel a spool holds, gathered from its patches.""" + return np.concatenate([x.get_array("distance") for x in spool]) + + +class TestChannelSelect: + """Selecting on the coordinates an inventory defines along the fiber.""" + + def test_trims_to_the_matching_channels(self, patch, inventory): + """ + A track name keeps the channels it covers and no others. + + The example path is coupled from 100 to 250 m along the fiber and + the patch's channels start at 100, so the trench runs from the + patch's channel 0 to its channel 150. + """ + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.select(coupling="trench") + assert len(out) == 1 + assert out[0].get_coord("distance").min() == 0 + assert out[0].get_coord("distance").max() == 150 + + def test_len_is_exact_before_anything_loads(self, patch, inventory): + """Selection is metadata work: the count is final without a read.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.select(coupling="trench") + contents = out.get_contents() + assert len(out) == len(contents) == 1 + assert contents["distance_max"].iloc[0] == 150.0 + + def test_the_data_is_the_channels_it_names(self, patch, inventory): + """ + The rows kept are the rows the coordinate names. + + Trimming the envelope without trimming the same rows out of the + array is the one failure which would not show up in the contents. + """ + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.select(coupling="trench")[0] + whole = patch.get_array("distance") + rows = np.searchsorted(whole, out.get_array("distance")) + assert np.array_equal(out.data, patch.data[rows]) + + def test_a_disjoint_match_subdivides(self, patch, two_zones): + """ + A group covering two stretches gives two patches, not one span. + + This is why selection builds a plan: `len` grows, and the hole + between the two runs belongs to neither. + """ + spool = dc.spool(patch).attach_inventory(two_zones) + out = spool.select(hole="a") + assert len(out) == 2 + contents = out.get_contents() + assert contents["distance_min"].tolist() == [10.0, 200.0] + assert contents["distance_max"].tolist() == [50.0, 240.0] + + def test_a_disjoint_match_keeps_each_piece_whole(self, patch, two_zones): + """Each piece holds exactly the rows its own envelope names.""" + spool = dc.spool(patch).attach_inventory(two_zones) + whole = patch.get_array("distance") + for piece in spool.select(hole="a"): + rows = np.searchsorted(whole, piece.get_array("distance")) + assert np.array_equal(piece.data, patch.data[rows]) + + def test_selection_agrees_with_enrichment(self, patch, two_zones): + """ + Every channel kept really does hold the value asked for. + + The two run on one projection, so this is the property which + makes that worth doing rather than a coincidence to maintain. + """ + spool = dc.spool(patch).attach_inventory(two_zones) + for piece in spool.select(hole="a").enrich(): + assert set(piece.get_array("hole")) == {"a"} + + def test_a_value_nothing_holds_keeps_nothing(self, patch, inventory): + """ + A patch with no matching channel is no more selected than one + which lacks the attr entirely. + """ + spool = dc.spool(patch).attach_inventory(inventory) + assert len(spool.select(coupling="cement")) == 0 + + def test_an_undescribed_patch_is_silently_not_selected(self, patch, inventory): + """The one agreed exception to loud-by-default.""" + other = patch.update_attrs(acquisition_key="DAS.R2D1..OTHER") + spool = dc.spool([patch, other]).attach_inventory(inventory) + with warnings.catch_warnings(): + warnings.simplefilter("error") + assert len(spool.select(coupling="trench")) == 1 + + def test_a_straddling_patch_is_not_selected(self, patch, inventory): + """ + Described twice is not one answer, so it is not selected. + + `conform_to_inventory` is what turns such a row into rows which + each resolve; select is a filter and says nothing about it. + """ + coord = patch.get_coord("time") + middle = coord.min() + (coord.max() - coord.min()) / 2 + assert coord.min() < middle <= coord.max() + split = _split_epochs(inventory, middle, second={"name": "moved"}) + spool = dc.spool(patch).attach_inventory(split) + assert len(spool.select(coupling="trench")) == 0 + + def test_conforming_first_makes_it_selectable(self, patch, inventory): + """The pieces each resolve, so each is judged on its own.""" + coord = patch.get_coord("time") + middle = coord.min() + (coord.max() - coord.min()) / 2 + split = _split_epochs(inventory, middle, second={"name": "moved"}) + spool = dc.spool(patch).attach_inventory(split) + out = spool.conform_to_inventory().select(coupling="trench") + assert len(out) == 2 + assert set(out.get_contents()["distance_max"]) == {150.0} + + def test_composes_with_a_time_split(self, patch, inventory): + """Both subdivisions hold, and the pieces still hold their own data.""" + coord = patch.get_coord("time") + middle = coord.min() + (coord.max() - coord.min()) / 2 + split = _split_epochs(inventory, middle, second={"name": "moved"}) + spool = dc.spool(patch).attach_inventory(split) + out = spool.conform_to_inventory().select(coupling="trench") + times, dists = patch.get_array("time"), patch.get_array("distance") + for piece in out: + rows = np.searchsorted(dists, piece.get_array("distance")) + cols = np.searchsorted(times, piece.get_array("time")) + assert np.array_equal(piece.data, patch.data[np.ix_(rows, cols)]) + + def test_two_selections_are_both_applied(self, patch, two_zones): + """ + Selecting twice narrows; it does not re-plan from the source. + + Chunking the same dimension twice replaces the first plan, which + is right for chunk and would silently undo a selection here. + """ + spool = dc.spool(patch).attach_inventory(two_zones) + out = spool.select(hole="a").select(zone="north") + assert len(out) == 1 + assert out.get_contents()["distance_max"].tolist() == [50.0] + + def test_kwargs_are_and(self, patch, two_zones): + """Two names in one call are judged together, as the index does.""" + spool = dc.spool(patch).attach_inventory(two_zones) + both = spool.select(hole="a", zone="north") + chained = spool.select(hole="a").select(zone="north") + assert both.get_contents()["distance_max"].tolist() == ( + chained.get_contents()["distance_max"].tolist() + ) + + def test_a_whole_match_changes_nothing(self, patch, inventory): + """Keeping every channel needs no plan, so the spool is unchanged.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.select(geometry="trench") + assert len(out) == 1 + assert out[0].shape == patch.shape + + def test_selector_shapes_match_the_attr_side(self, patch, inventory): + """A glob, a sequence, and a scalar all mean what they mean there.""" + spool = dc.spool(patch).attach_inventory(inventory) + scalar = spool.select(coupling="trench").get_contents()["distance_max"] + for selector in ("tren*", "trenc?", ["trench", "conduit"]): + got = spool.select(coupling=selector).get_contents()["distance_max"] + assert got.tolist() == scalar.tolist() + + def test_a_numeric_range_selects_a_stretch(self, patch, inventory): + """A qualified numeric field takes a range, as a coordinate does.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.select(**{"optical_components.optical_length": (400.0, 600.0)}) + assert len(out) == 1 + + def test_none_matches_the_undefined_channels(self, patch, inventory): + """ + `None` is how a query spells the marker absence is stored as. + + The path is coupled only to 250 m, so the channels past it hold + the empty string a string coordinate has instead of a null. + """ + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.select(coupling=None) + assert len(out) == 1 + assert out.get_contents()["distance_min"].tolist() == [151.0] + + def test_a_membership_group_takes_a_boolean(self, patch, inventory): + """A membership group is False where nothing includes it.""" + spool = dc.spool(patch).attach_inventory(inventory) + noisy = spool.select(noisy=True).get_contents() + assert noisy["distance_min"].tolist() == [50.0] + assert len(spool.select(noisy=False)) == 2 + + def test_a_name_this_path_defines_nowhere_matches_nothing(self, patch, inventory): + """ + Listing a name is not promising a value for it. + + One path recording a track makes the name selectable for the + whole inventory, so a patch whose own path records none must + answer with no channel rather than with an error. + """ + array = inventory.networks[0].fiber_arrays[0] + acquisition, path = array.acquisitions[0], array.optical_paths[0] + both = inventory.replace( + array, + array.new( + acquisitions=(acquisition, acquisition.new(location_code="01")), + optical_paths=(path, path.new(location_code="01", coupling=())), + ), + ) + assert "coupling" in both.get_names().coords + other = patch.update_attrs(acquisition_key="DAS.R2D1.01.RAW", tag="second") + spool = dc.spool([patch, other]).attach_inventory(both) + # The first patch's path is coupled; the second's records none. + out = spool.select(coupling="trench") + assert out.get_contents()["tag"].tolist() == ["random"] + + def test_the_patch_axis_keeps_its_own_meaning(self, patch, inventory): + """ + `distance` is the patch's axis, inventory attached or not. + + The inventory could also place it on the fiber — its channels + start at 100 m along the path — so a name the index already uses + for a coordinate has to keep it, or attaching would move a name + out of the namespace it has always been in. + """ + spool = dc.spool(patch).attach_inventory(inventory) + for form in ({}, {"_coords": {"distance": (0, 100)}}): + kwargs = {"distance": (0, 100)} if not form else {} + out = spool.select(**form, **kwargs) + assert out.get_contents()["distance_max"].tolist() == [100.0] + + def test_samples_with_a_channel_name_raises(self, patch, inventory): + """A fiber coordinate has no sample numbering of its own.""" + spool = dc.spool(patch).attach_inventory(inventory) + with pytest.raises(InvalidSpoolQueryError, match="no sample numbering"): + spool.select(coupling="trench", samples=True) + + def test_a_channel_name_through_attrs_raises(self, patch, inventory): + """It describes channels, so the attrs namespace is the wrong one.""" + spool = dc.spool(patch).attach_inventory(inventory) + with pytest.raises(InvalidSpoolQueryError, match="along the fiber"): + spool.select(_attrs={"coupling": "trench"}) + + def test_an_acquisition_with_no_map_refuses(self, patch, inventory): + """Its channels cannot be placed, and guessing would trim wrongly.""" + without = _replace_acquisition(inventory, distance_map=None) + spool = dc.spool(patch).attach_inventory(without) + with pytest.raises(PatchError, match="no distance_map"): + spool.select(coupling="trench") + + def test_a_patch_without_the_axis_refuses(self, inventory, patch): + """A patch carrying no dimension the map places channels by.""" + lag = patch.rename_coords(distance="offset") + spool = dc.spool(lag).attach_inventory(inventory) + with pytest.raises(PatchError, match="places channels by"): + spool.select(coupling="trench") + + +class TestChannelUnselect: + """Removing the channels a selection would have kept.""" + + def test_complements_the_selection(self, patch, two_zones): + """Together the two hold every channel, and share none.""" + spool = dc.spool(patch).attach_inventory(two_zones) + kept = _channels(spool.select(hole="a")) + dropped = _channels(spool.unselect(hole="a")) + assert not set(kept) & set(dropped) + assert sorted([*kept, *dropped]) == sorted(patch.get_array("distance")) + + def test_removes_the_undefined_channels(self, patch, inventory): + """The spec's own example: drop the channels with no coupling.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.unselect(coupling=None) + assert out.get_contents()["distance_max"].tolist() == [150.0] + + def test_an_undescribed_patch_is_kept_whole(self, patch, inventory): + """The selection never held it, so its complement keeps all of it.""" + other = patch.update_attrs(acquisition_key="DAS.R2D1..OTHER") + spool = dc.spool([patch, other]).attach_inventory(inventory) + out = spool.unselect(coupling="trench") + assert len(out) == 2 + assert sorted(out.get_contents()["distance_min"]) == [0.0, 151.0] + + def test_attrs_and_channels_stay_one_complement(self, patch, inventory): + """ + A patch the attrs never matched keeps every channel. + + Complementing the two halves apart would drop it, since the + selection it is the complement of never held it. + """ + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.unselect(gauge_length=99.0, coupling="trench") + assert len(out) == 1 + assert out[0].shape == patch.shape + + def test_attrs_and_channels_trim_the_matched_patch(self, patch, inventory): + """And one the attrs did match loses the channels which matched.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.unselect(gauge_length=10.0, coupling="trench") + assert out.get_contents()["distance_min"].tolist() == [151.0] + + def test_a_patch_coordinate_still_raises(self, patch, inventory): + """A real coordinate range is a trim, and unselect is not one.""" + spool = dc.spool(patch).attach_inventory(inventory) + with pytest.raises(InvalidSpoolQueryError, match="unselect cannot take"): + spool.unselect(_coords={"distance": (0, 100)}) + + +class TestSplitBy: + """Expanding a spool into one patch per value along the fiber.""" + + def test_one_patch_per_value(self, patch, inventory): + """The example path annotates two zones, so two patches come out.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("zone") + assert len(out) == 2 + assert out.get_contents()["zone"].tolist() == ["north", "south"] + + def test_the_pieces_hold_their_own_channels(self, patch, inventory): + """Each output holds the rows its value covers, and no others.""" + spool = dc.spool(patch).attach_inventory(inventory) + whole = patch.get_array("distance") + for piece in spool.split_by("zone"): + rows = np.searchsorted(whole, piece.get_array("distance")) + assert np.array_equal(piece.data, patch.data[rows]) + + def test_the_value_is_stamped_on_the_patch(self, patch, inventory): + """So overlapping siblings stay apart once they are patches.""" + spool = dc.spool(patch).attach_inventory(inventory) + assert [x.attrs.zone for x in spool.split_by("zone")] == ["north", "south"] + + def test_stamp_false_leaves_the_attrs_alone(self, patch, inventory): + """Which is what a nested split wants of the second one.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("zone", stamp=False) + assert "zone" not in out.get_contents().columns + assert not any(dict(x.attrs).get("zone") for x in out) + + def test_the_stamp_is_selectable(self, patch, inventory): + """A stamped value is an ordinary attr, so it filters like one.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("zone") + assert len(out.select(zone="north")) == 1 + + def test_a_membership_group_splits_in_two(self, patch, inventory): + """Both sides come out: the channels included and those not.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("noisy") + assert sorted(out.get_contents()["noisy"].tolist()) == [False, False, True] + + def test_overlapping_values_share_channels(self, patch, inventory): + """ + A group may overlap another, so a channel can land in two outputs. + + `noisy` runs from 150 to 300 m and the zones meet at 200, so the + two split differently over the same fiber. + """ + spool = dc.spool(patch).attach_inventory(inventory) + zones = _channels(spool.split_by("zone")) + noisy = _channels(spool.split_by("noisy")) + assert sorted(zones) == sorted(patch.get_array("distance")) + assert sorted(noisy) == sorted(patch.get_array("distance")) + + def test_a_disjoint_value_becomes_several_patches(self, patch, two_zones): + """One value covering two stretches keeps them apart.""" + spool = dc.spool(patch).attach_inventory(two_zones) + out = spool.split_by("hole") + assert len(out) == 2 + assert out.get_contents()["hole"].tolist() == ["a", "a"] + + def test_include_keeps_only_what_it_names(self, patch, inventory): + """The globs read the value written as a string.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("zone", include="nor*") + assert out.get_contents()["zone"].tolist() == ["north"] + + def test_exclude_wins_over_include(self, patch, inventory): + """So naming a family and carving one out of it reads either way.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("zone", include=("nor*", "sou*"), exclude="north") + assert out.get_contents()["zone"].tolist() == ["south"] + + def test_a_group_no_path_defines_yields_nothing(self, patch, inventory): + """There is no value to split on, so there is no output.""" + spool = dc.spool(patch).attach_inventory(inventory) + assert len(spool.split_by("not_a_group")) == 0 + + def test_needs_an_inventory(self, patch): + """The values it expands into are ones an inventory states.""" + with pytest.raises(ParameterError, match="needs an inventory"): + dc.spool(patch).split_by("zone") + + def test_a_nested_split_keeps_the_first_stamp(self, patch, inventory): + """Which is what `stamp=False` is for.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("zone").split_by("noisy", stamp=False) + assert set(out.get_contents()["zone"]) == {"north", "south"} + + +class TestChannelSelectEdges: + """Rows a fiber query cannot answer for, and the ways it says so.""" + + def test_a_lag_time_patch_has_no_epoch_to_resolve_at(self, patch, inventory): + """ + Without instants there is no context, so nothing is selected. + + The same thing `Patch.enrich` refuses to guess at: a correlation's + lag times are not the moments the fiber was in some state. + """ + lags = patch.get_coord("time").values - patch.get_coord("time").min() + spool = dc.spool(patch.update_coords(time=lags)).attach_inventory(inventory) + assert len(spool.select(coupling="trench")) == 0 + + def test_an_empty_spool_stays_empty(self, patch, inventory): + """With no rows there is nothing to resolve, and no work to do.""" + spool = dc.spool(patch).attach_inventory(inventory).select(tag="nope") + assert len(spool) == 0 + assert len(spool.select(coupling="trench")) == 0 + + def test_an_unevenly_sampled_patch_refuses(self, patch, inventory): + """ + Which channels match is decided on the sample grid. + + A patch with a hole in its distance coordinate records no + spacing, so the grid cannot be rebuilt; trimming the wrong + channels quietly is worse than saying so. + """ + holed = patch.unselect(distance=(50, 200)) + assert holed.get_coord("distance").step is None + spool = dc.spool(holed).attach_inventory(inventory) + with pytest.raises(PatchError, match="no channel spacing"): + spool.select(coupling="trench") + + def test_patches_on_different_channel_dimensions_refuse(self, patch, inventory): + """ + One spool, one dimension to trim: two is no answer at all. + + The patch-level twin refuses the same shape, where one patch + carries two coordinates the map could be read on. + """ + distance = patch.get_coord("distance") + both_axes = _replace_acquisition( + inventory, + distance_map=DistanceMap( + channel=(float(distance.min()), float(distance.max())), + instrument_distance=(float(distance.min()), float(distance.max())), + distance=(100.0, 100.0 + float(distance.max() - distance.min())), + ), + ) + channels = dc.Patch( + data=patch.data, + coords={ + "channel": patch.get_array("distance"), + "time": patch.get_array("time"), + }, + dims=("channel", "time"), + attrs=patch.attrs, + ) + spool = dc.spool([patch, channels]).attach_inventory(both_axes) + with pytest.raises(InvalidSpoolQueryError, match="different"): + spool.select(coupling="trench") + + def test_none_on_a_membership_group(self, patch, inventory): + """A membership group says something about every channel: False.""" + spool = dc.spool(patch).attach_inventory(inventory) + undefined = spool.select(noisy=None).get_contents() + assert undefined["distance_min"].tolist() == ( + spool.select(noisy=False).get_contents()["distance_min"].tolist() + ) + + def test_none_on_a_numeric_group(self, patch, inventory): + """A numeric group spells absence NaN, which no range matches.""" + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + numeric = inventory.replace( + path, + path.new( + annotations=( + *path.annotations, + OpticalPathAnnotation( + start_distance=100.0, + end_distance=200.0, + group="frost_depth", + value=1.5, + ), + ) + ), + ) + spool = dc.spool(patch).attach_inventory(numeric) + # The group covers the first hundred channels, so the rest are NaN. + assert spool.select(frost_depth=None).get_contents()[ + "distance_min" + ].tolist() == [101.0] + assert spool.select(frost_depth=(1.0, 2.0)).get_contents()[ + "distance_max" + ].tolist() == [100.0] + + def test_splitting_an_undescribed_spool_yields_nothing(self, patch, inventory): + """No fiber to split on means no output, not an error.""" + other = patch.update_attrs(acquisition_key="DAS.R2D1..OTHER") + spool = dc.spool(other).attach_inventory(inventory) + assert len(spool.split_by("zone")) == 0 + + def test_splitting_skips_a_row_it_cannot_place(self, patch, inventory): + """A described patch splits; one the inventory is silent about does not.""" + other = patch.update_attrs(acquisition_key="DAS.R2D1..OTHER", tag="second") + spool = dc.spool([patch, other]).attach_inventory(inventory) + out = spool.split_by("zone") + assert out.get_contents()["zone"].tolist() == ["north", "south"] + + def test_splitting_an_unevenly_sampled_patch_refuses(self, patch, inventory): + """The grid is what values are read on here too.""" + holed = patch.unselect(distance=(50, 200)) + spool = dc.spool(holed).attach_inventory(inventory) + with pytest.raises(PatchError, match="no channel spacing"): + spool.split_by("zone") + + def test_one_patch_carrying_two_channel_axes_refuses(self, patch, inventory): + """ + Both are dimensions, so neither is obviously the channel one. + + The patch-level resolver settles this by projecting each and + checking they agree; two *dimensions* are different axes rather + than two spellings of one, so there is nothing to check. + """ + distance = patch.get_coord("distance") + both_axes = _replace_acquisition( + inventory, + distance_map=DistanceMap( + channel=(float(distance.min()), float(distance.max())), + instrument_distance=(float(distance.min()), float(distance.max())), + distance=(100.0, 100.0 + float(distance.max() - distance.min())), + ), + ) + stacked = patch.data[:3][..., np.newaxis] + cube = dc.Patch( + data=np.broadcast_to(stacked, (3, patch.shape[1], 2)).copy(), + coords={ + "channel": np.arange(3.0), + "time": patch.get_array("time"), + "distance": np.arange(2.0), + }, + dims=("channel", "time", "distance"), + attrs=patch.attrs, + ) + spool = dc.spool(cube).attach_inventory(both_axes) + with pytest.raises(PatchError, match="ambiguous"): + spool.select(coupling="trench") + + def test_splitting_a_lag_time_patch_yields_nothing(self, patch, inventory): + """Without instants there is no context to read a group from.""" + lags = patch.get_coord("time").values - patch.get_coord("time").min() + spool = dc.spool(patch.update_coords(time=lags)).attach_inventory(inventory) + assert len(spool.split_by("zone")) == 0 diff --git a/tests/test_utils/test_chunk.py b/tests/test_utils/test_chunk.py index 487723b8..e5c867ff 100644 --- a/tests/test_utils/test_chunk.py +++ b/tests/test_utils/test_chunk.py @@ -665,7 +665,7 @@ class TestBuildSubdivisionPlan: """Building a plan from pieces given directly, as a selection does.""" def test_a_row_given_no_pieces_leaves(self): - """This is how a selection drops a patch it keeps no sample of.""" + """A selection drops a patch this way when it keeps no sample.""" df = _one_row_df(0.0, 1.0, 10, name="distance") plan = build_subdivision_plan(df, [[]], "distance") assert not len(plan.outputs) From c83531eac3fc041bcc2c5d24498de6bce733377d Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 22:33:10 +0200 Subject: [PATCH 6/8] Address the adversarial review of the fiber operations Three reviewers independently found the same defect, which is the worst of these: the channel grid was rebuilt with the signed step, so a reverse-sorted patch counted to a negative number of samples, got an empty grid, and was silently dropped by every fiber query. The sibling in chunk_plan.py takes the magnitude for exactly this reason. Worse in kind, though only one reviewer got near it: re-planning the same dimension collapses onto the sources, which is sound only while a plan's pieces cover them. A selection's do not, so select(...).chunk(distance=...) loaded back the channels it had removed while the contents went on describing the ones it kept. A plan now records whether it drops samples, read off the pieces rather than taken on trust, and a lossy one never collapses. The rest, each with a test: a pathless acquisition is valid and projects nothing rather than being dereferenced; the units the inventory documents for a field reach the predicate, so a metre selector no longer meets a unitless value; unselect strips channel names from _coords as select does, so both spellings work; a bare ... selects everything here as everywhere; samples and relative are refused beside a fiber name rather than ignored; split_by refuses a name the inventory could not contribute instead of returning an empty spool; and Patch.unselect handles two coordinates on one dimension and a dimension carrying no values of its own. Prose: split_by claimed a channel could land in two outputs of one call. It cannot -- a channel resolves to one value of a group, the same one enrichment projects -- and the test named for it asserted the partition it actually is. Both now say so, here and in the spec. BaseSpool.unselect still said coordinates were refused outright. Tests: five asserted only inside a loop over a spool whose length they never pinned, so an implementation selecting nothing passed them; the complement tests were satisfied by keeping nothing and dropping everything; and no test was red without the grid arithmetic this branch changed. That last one is now a property over 900 random float grids, verified red with the correction loops removed. --- dascore/core/spool.py | 168 +++++++++---- dascore/io/index/planned.py | 37 ++- dascore/proc/coords.py | 62 +++-- dascore/proc/inventory.py | 132 ++++++---- dascore/utils/chunk_plan.py | 7 +- docs/tutorial/patch.qmd | 2 +- tests/test_proc/test_proc_coords.py | 17 +- tests/test_proc/test_proc_inventory.py | 323 ++++++++++++++++++++++--- tests/test_utils/test_chunk.py | 64 ++++- 9 files changed, 643 insertions(+), 169 deletions(-) diff --git a/dascore/core/spool.py b/dascore/core/spool.py index a6cab27e..17194e2d 100644 --- a/dascore/core/spool.py +++ b/dascore/core/spool.py @@ -8,6 +8,7 @@ from collections.abc import Callable, Generator, Iterator, Mapping, Sequence from dataclasses import replace from functools import singledispatch +from itertools import pairwise from pathlib import Path from typing import TYPE_CHECKING, ClassVar, Literal, TypeVar, overload @@ -98,6 +99,13 @@ def _copy_public_dataframe(frame: pd.DataFrame) -> pd.DataFrame: # policy here removes the patch, and only the noise differs. _VALID_ON_UNRESOLVED_CONFORM = ("raise", "warn", "drop") +# Written once because both fiber verbs refuse for it, and two spellings +# would let them start explaining the same refusal differently. +_UNPLACEABLE = ( + "are described by the inventory but cannot have their channels placed " + "along the fiber" +) + _UNRESOLVED_WARNING = ( "The attached inventory does not describe every patch in this spool, and " "those it does not describe were not enriched. Use on_unresolved='raise' " @@ -234,6 +242,47 @@ def _requested_names(_attrs, _coords, kwargs) -> set[str]: return out +def _drops_samples(rows: pd.DataFrame, pieces, name: str) -> bool: + """ + Return whether any row keeps less than the whole of itself. + + Read off the pieces rather than taken on trust, because getting it + wrong is silent either way: a plan wrongly called lossless loads back + the samples it dropped, and one wrongly called lossy only forgoes a + collapse. Conform's pieces always cover their row, so it answers + False here without the caller having to say so. + """ + lows = rows[f"{name}_min"].to_numpy() + highs = rows[f"{name}_max"].to_numpy() + steps = rows[f"{name}_step"].to_numpy() + for row_pieces, low, high, step in zip(pieces, lows, highs, steps, strict=True): + if not row_pieces: + continue # a row kept out entirely leaves no outputs to collapse + if row_pieces[0][0] != low or row_pieces[-1][1] != high: + return True + # Conform's pieces meet exactly one step apart, so they cover the + # row between them; a selection's have the channels it dropped in + # the gaps, which is the whole difference. + gap = abs(step) + if any(b[0] != a[1] + gap for a, b in pairwise(row_pieces)): + return True + return False + + +def _stated_channels(channels: dict) -> dict: + """ + Return the channel selectors which actually select something. + + A bare `...` selects everything here as everywhere, so it asks for no + trimming at all — but the name still had to be recognized as the + inventory's, or the index would be left to complain that it has never + heard of it. `None` is not dropped beside it: on a coordinate the + fiber defines it spells the undefined marker, which is a statement + about which channels to keep rather than the absence of one. + """ + return {name: value for name, value in channels.items() if value is not Ellipsis} + + def _glob_filter(include, exclude): """ Return a predicate deciding which split values to keep. @@ -263,6 +312,11 @@ def keep(value) -> bool: return keep +def _without_keys(mapping: Mapping, names) -> dict: + """Return the mapping without the named keys.""" + return {str(k): v for k, v in mapping.items() if k not in names} + + def _without_names(spec: namespace_select_type, names) -> namespace_select_type: """ Return an `_attrs`/`_coords` argument with some names taken out. @@ -275,7 +329,7 @@ def _without_names(spec: namespace_select_type, names) -> namespace_select_type: if not names or spec is None: return spec if isinstance(spec, Mapping): - return {str(k): v for k, v in spec.items() if k not in names} + return _without_keys(spec, names) if isinstance(spec, str): return None if spec in names else spec kept = [x for x in spec if x not in names] @@ -657,12 +711,17 @@ def unselect( one bad tag, an instrument being serviced — without spelling the rest of the archive as a selection. - Coordinates are not accepted yet. Selecting on one trims each - patch to the range rather than choosing between patches, so the - complement is every patch cut into the pieces outside it — one - patch becoming two. That is subdivision rather than filtering, - and it needs machinery this does not have; until it does, select - the ranges to keep. + The patches' own coordinates are not accepted. Selecting on one + trims each patch to the range rather than choosing between + patches, so the complement is every patch cut into the pieces + outside it — one patch becoming two. Select the ranges to keep + instead, or use [`Patch.unselect`](`dascore.Patch.unselect`) on + each patch, which can take samples out of its middle. + + The coordinates an attached inventory defines along the fiber are + different, and are accepted: removing one of those chooses which + channels a patch holds. A patch may then be cut into the pieces + the query did not match, so `len` can grow here as well. Naming nothing raises, and so does naming only no-op selectors (`None`, `...`). `select()` with no selection is the whole spool, @@ -676,7 +735,9 @@ def unselect( Attribute selections, in the forms [`select`](`dascore.core.spool.Spool.select`) accepts. _coords - Accepted only to say so; see above. + The patches' own coordinates only to say they are refused; a + coordinate an attached inventory defines along the fiber is + accepted and chooses channels. See above. **kwargs The selection whose matches are removed. @@ -970,7 +1031,7 @@ def select( ) -> Self: """{doc}.""" attr_query, channel_query, _attrs, _coords, kwargs = ( - self._split_inventory_query(_attrs, _coords, kwargs, samples) + self._split_inventory_query(_attrs, _coords, kwargs, samples, relative) ) catalog = self._catalog.select( _attrs=_attrs, @@ -982,7 +1043,7 @@ def select( out = self._new_from_catalog(catalog) if attr_query: out = out._select_from_inventory(attr_query) - if channel_query: + if channel_query := _stated_channels(channel_query): out = out._select_channels(channel_query) return out @@ -1008,8 +1069,8 @@ def unselect( known_attrs | selectable, known_coords, _attrs=_attrs, - _coords=_coords, - kwargs={k: v for k, v in kwargs.items() if k not in channels}, + _coords=_without_names(_coords, channels), + kwargs=_without_keys(kwargs, channels), ) if coords: msg = ( @@ -1025,7 +1086,7 @@ def unselect( # A no-op selector selects everything, so its complement is an # empty spool -- and "remove nothing" reads just as naturally. stated = {k: v for k, v in attrs.items() if v is not None and v is not Ellipsis} - if not stated and not channels: + if not stated and not _stated_channels(channels): msg = ( "unselect needs something to remove; " f"{sorted(requested) or 'nothing'} names no selection. " @@ -1044,7 +1105,9 @@ def unselect( # apart would drop a patch the whole selection never held. matched = self if not stated else self.select(_attrs=stated) return self._select_channels( - channels, complement=True, applies_to=matched._catalog._ordered_ids() + _stated_channels(channels), + complement=True, + applies_to=matched._catalog._ordered_ids(), ) def _index_names(self) -> tuple[set[str], set[str]]: @@ -1086,9 +1149,12 @@ def _channel_query( "channels rather than whole patches." ) raise InvalidSpoolQueryError(msg) - return {name: stated[name] for name in wanted} + # Sorted so a query built from a set does not order itself by hash: + # the masks are combined with AND either way, but which selector a + # message complains about first should not move between runs. + return {name: stated[name] for name in sorted(wanted)} - def _split_inventory_query(self, _attrs, _coords, kwargs, samples): + def _split_inventory_query(self, _attrs, _coords, kwargs, samples, relative=False): """ Split selectors into the ones the index answers and the rest. @@ -1104,15 +1170,22 @@ def _split_inventory_query(self, _attrs, _coords, kwargs, samples): channels = self._channel_query( requested, names, known_attrs, known_coords, _coords, kwargs ) - if channels and samples: + # Neither keyword has anything to mean about a value the fiber + # states: it has no sample numbering of its own -- the channels it + # describes are the patch's -- and no endpoints to be relative to, + # since what it says varies from one acquisition to the next. + # Ignoring either quietly would answer a question nobody asked. + for flag, label in ((samples, "samples"), (relative, "relative")): + if not (channels and flag): + continue msg = ( f"{sorted(channels)} name coordinates the inventory defines " - "along the fiber, which have no sample numbering of their " - "own: the channels they describe are the patch's, and " - "samples=True asks about that axis instead." + f"along the fiber, which {label}=True cannot describe: it " + "asks about the patch's own axis, and these say what is " + "attached to each channel of it." ) raise InvalidSpoolQueryError(msg) - kwargs = {k: v for k, v in kwargs.items() if k not in channels} + kwargs = _without_keys(kwargs, channels) _coords = _without_names(_coords, channels) # A name the index already uses for a coordinate keeps its meaning; # bare names resolve to attrs first, and an inventory must not @@ -1177,12 +1250,7 @@ def _select_channels(self, query: dict, *, complement=False, applies_to=None): name, pieces, reasons = resolve_channel_pieces( self._inventory, contexts, working, query, complement=complement ) - _refuse_rows( - source_rows, - reasons, - "are described by the inventory but cannot have their channels " - "placed along the fiber", - ) + _refuse_rows(source_rows, reasons, _UNPLACEABLE) if name is None: # No row has a fiber to be judged along, so the query matched # nothing: an empty spool, or the whole of it complemented. @@ -1458,21 +1526,25 @@ def split_by( Most often an annotation group. Every kind of group splits: a categorical one by each of its strings, a membership group into the channels it includes and those it does not, and a numeric one - by each distinct measurement. Groups are allowed to overlap, so a - channel may appear in more than one output patch, and a patch - whose channels take several values becomes several patches — this - can greatly expand the spool. + by each distinct measurement. Intervals of one group may overlap, + but a channel still holds only one of its values, so the outputs + of one call divide the fiber rather than share it. A patch whose + channels take several values becomes several patches — this can + greatly expand the spool. Parameters ---------- name The inventory-derived coordinate to split on. include, exclude - Glob patterns matched against each value written as a string, - so `"trench_?a"` and `"hole_*"` both work and `True` is - matched as `"True"`. With `include`, only values matching one - of them are kept; `exclude` drops the values it matches, and - wins where both match. + Glob patterns matched against each value *written as a + string*, which is what lets one spelling cover all three + kinds of group: `"hole_*"` reads a categorical one, `"Tru*"` + a membership one, and `"1.*"` a numeric one. Selecting on the + stamp afterwards compares typed values instead, so the two + are not interchangeable. With `include`, only values matching + one of them are kept; `exclude` drops the values it matches, + and wins where both match. stamp Whether to record the value on each output patch as an attr named after the coordinate, so overlapping siblings stay @@ -1505,17 +1577,23 @@ def split_by( "fiber. Attach one with Spool.attach_inventory." ) raise ParameterError(msg) + # A name the inventory could not contribute has no values to split + # into, so it would quietly give an empty spool. Selection refuses + # a name it does not know, and a misspelling is no more meaningful + # here than it is there. + if name not in set(self._inventory.get_names().coords): + msg = ( + f"{name!r} is not a coordinate the attached inventory defines " + "along the fiber, so there is nothing to split on. " + "Inventory.get_names().coords lists the names it could." + ) + raise InvalidSpoolQueryError(msg) source_rows, working = self._plan_frames() contexts = self._plan_contexts(working) dim, rows, reasons = resolve_split_pieces( self._inventory, contexts, working, name, _glob_filter(include, exclude) ) - _refuse_rows( - source_rows, - reasons, - "are described by the inventory but cannot have their channels " - "placed along the fiber", - ) + _refuse_rows(source_rows, reasons, _UNPLACEABLE) if dim is None: # nothing to split: no row has a fiber to split on return self._restrict_to_rows([]) pieces = [[piece for _, piece in row] for row in rows] @@ -1655,6 +1733,11 @@ def _subdivided(self, sources, rows, pieces, name: str, stamp=None) -> Self: exact — and loading goes through the machinery which already trims a member on extraction. `stamp` names an attr to record on each output, in the order the pieces were given. + + A plan whose pieces do not cover their rows is marked lossy, so + that re-planning the same dimension nests rather than collapsing + onto the sources — which would load back the samples the pieces + left out. """ from dascore.io.index.planned import derived_catalog # noqa: PLC0415 @@ -1672,6 +1755,7 @@ def _subdivided(self, sources, rows, pieces, name: str, stamp=None) -> Self: mode="chunk", origin_path=self.spool_path, stamped=stamped, + lossy=_drops_samples(rows, pieces, name), ) return self._new_from_catalog(catalog) diff --git a/dascore/io/index/planned.py b/dascore/io/index/planned.py index d98dadfe..0d86dd5d 100644 --- a/dascore/io/index/planned.py +++ b/dascore/io/index/planned.py @@ -368,6 +368,7 @@ def __init__( check_behavior: WARN_LEVELS = "warn", origin_path=None, stamped: tuple[str, ...] = (), + lossy: bool = False, ): if "output_id" not in member_rows.columns: msg = "member_rows must carry an output_id column." @@ -385,6 +386,11 @@ def __init__( self.origin_path = origin_path # attrs the outputs state about themselves rather than inherit self.stamped = tuple(stamped) + # Whether the outputs leave samples of their sources out. A lossy + # plan must never be collapsed: its members do not cover their + # sources, so re-planning over them would load back what it + # dropped. See `collapse_working_df`. + self.lossy = bool(lossy) def live_entries(self) -> dict[str, dc.Patch]: """Expose the loader's live registry (for absorption/transfer).""" @@ -463,20 +469,22 @@ def resolve(self, row: Mapping, **trim) -> dc.Patch: if self.mode == "identity": # one untouched member per output; residuals apply at load assert len(members) == 1 - return self._stamp(self._load_member(members.iloc[0].to_dict()), row) - if self.mode == "concat": - patches = [ + patch = self._load_member(members.iloc[0].to_dict()) + elif self.mode == "concat": + loaded = [ self._load_member(kwargs) for kwargs in members.to_dict("records") ] out = concatenate_patches( - patches, check_behavior=self.check_behavior, **{self.dim: None} + loaded, check_behavior=self.check_behavior, **{self.dim: None} ) assert len(out) == 1 - return self._stamp(out[0], row) - joined = members.assign(current_index=output_id) - patches = self._assembler()._patch_from_instruction_df(joined) - assert len(patches) == 1 - return self._stamp(patches[0], row) + patch = out[0] + else: + joined = members.assign(current_index=output_id) + assembled = self._assembler()._patch_from_instruction_df(joined) + assert len(assembled) == 1 + patch = assembled[0] + return self._stamp(patch, row) def _stamp(self, patch: dc.Patch, row: Mapping) -> dc.Patch: """ @@ -525,6 +533,7 @@ def derived_catalog( check_behavior: WARN_LEVELS = "warn", origin_path=None, stamped: tuple[str, ...] = (), + lossy: bool = False, ) -> PatchCatalog: """ Materialize a plan into a fresh in-memory catalog. @@ -594,6 +603,7 @@ def derived_catalog( check_behavior=check_behavior, origin_path=origin_path, stamped=stamped, + lossy=lossy, ) backend = get_backend(":memory:") coord_dims_map = {} if parent is None else parent.backend.coord_dims_map() @@ -627,9 +637,16 @@ def collapse_working_df(catalog: PatchCatalog) -> pd.DataFrame | None: applied to the envelopes. (Planning a different dimension must keep the assembled boundaries, so its caller plans over the output rows instead and never collapses.) + + A *lossy* plan is the exception and never collapses. Collapsing is + sound because the members of a chunk or a subdivision together cover + their sources, so a re-plan which merges them back is entitled to + load a source whole. A plan which drops samples — channel selection + keeping some channels of a patch and not others — breaks exactly + that, and collapsing it would quietly load back what it removed. """ resolver = catalog.resolver - if not isinstance(resolver, PlanResolver): + if not isinstance(resolver, PlanResolver) or resolver.lossy: return None members = resolver.member_rows if catalog.is_view: diff --git a/dascore/proc/coords.py b/dascore/proc/coords.py index e60de253..ff4af15e 100644 --- a/dascore/proc/coords.py +++ b/dascore/proc/coords.py @@ -442,6 +442,17 @@ def coords_from_df( return out +def _check_coord_names(patch: PatchType, kwargs) -> None: + """Refuse a name the patch has no coordinate for, naming what it has.""" + if not (invalid := set(kwargs) - set(patch.coords.coord_map)): + return + valid_list = sorted(patch.coords.coord_map) + msg = ( + f"Coordinate(s) {sorted(invalid)} not found in patch coordinates: {valid_list}" + ) + raise PatchCoordinateError(msg) + + @patch_function(history=None) @compose_docstring(select_params=select_values_description) def select( @@ -543,15 +554,7 @@ def select( 1999 """ - # Check for and raise on invalid kwargs. - if invalid_coords := set(kwargs) - set(patch.coords.coord_map): - invalid_list = sorted(invalid_coords) - valid_list = sorted(patch.coords.coord_map) - msg = ( - f"Coordinate(s) {invalid_list} not found in patch coordinates: {valid_list}" - ) - raise PatchCoordinateError(msg) - + _check_coord_names(patch, kwargs) new_coords, data = patch.coords.select( **kwargs, array=patch.data, @@ -615,28 +618,37 @@ def unselect( - Each named coordinate is complemented on its own. Selecting on two coordinates keeps the samples in both ranges, and everything - outside that is an L rather than a block, which no array can hold — - so `unselect` removes each named range instead, which is the part - of the complement that is expressible. - """ - if invalid_coords := set(kwargs) - set(patch.coords.coord_map): - invalid_list = sorted(invalid_coords) - valid_list = sorted(patch.coords.coord_map) - msg = ( - f"Coordinate(s) {invalid_list} not found in patch coordinates: {valid_list}" - ) - raise PatchCoordinateError(msg) - complements = {} + outside that is a frame around them rather than a block, which no + array can hold — so `unselect` removes each named range instead, + which is the part of the complement that is expressible. Two + coordinates along one dimension therefore both take their range + out of it, leaving what neither removed. + """ + _check_coord_names(patch, kwargs) + keep: dict[str, np.ndarray] = {} for name, value in kwargs.items(): coord = patch.coords.coord_map[name] + dims = patch.coords.dim_map[name] + if len(dims) != 1: + msg = ( + f"Coordinate {name!r} spans {list(dims)}, so removing a range " + "of it does not name samples of one dimension to drop." + ) + raise PatchCoordinateError(msg) # Asking select itself which samples it would keep is what stops # the two from drifting: one selector cannot come to mean # different things in select and its complement. _, indexer = coord.select(value, relative=relative, samples=samples) - kept = np.zeros(len(coord), dtype=bool) - kept[indexer] = True - complements[name] = ~kept - return patch.select(**complements, copy=copy) + selected = np.zeros(len(coord), dtype=bool) + selected[indexer] = True + keep[dims[0]] = keep.get(dims[0], True) & ~selected + # Kept as sample numbers along each dimension rather than as a mask + # per coordinate: coordinates sharing a dimension are applied in + # separate passes, so the second mask would meet an already trimmed + # axis, and a dimension carrying no values of its own takes samples + # where it would refuse an array. + trims = {dim: np.flatnonzero(mask) for dim, mask in keep.items()} + return patch.select(**trims, samples=True, copy=copy) @patch_function(history=None) diff --git a/dascore/proc/inventory.py b/dascore/proc/inventory.py index e61bf3a6..7c26654b 100644 --- a/dascore/proc/inventory.py +++ b/dascore/proc/inventory.py @@ -543,14 +543,34 @@ def _apply_conflicts(patch, new_attrs, conflicts) -> tuple[dict, list]: assert set(_AXIS_COORDS) == set(DISTANCE_MAP_AXES) -def _get_channel_axes(patch, acquisition) -> list[tuple[str, str]]: +def _map_axis_coords(dist_map, names) -> list[tuple[str, str]]: """ - Return every map axis the patch can be read on, with its coordinate. - - The map states its control points in whichever coordinates were - measured, so the patch decides which of them is read. Guessing instead - would be wrong by the channel spacing, silently. + Return each map axis paired with the name it can be read on. + + One name per axis, the first of them present: the map states its + control points in whichever coordinates were measured, so what is + carried decides which is read. Guessing instead would be wrong by the + channel spacing, silently. Stated once because a patch and an index + row must agree about it — selection answering differently than + enrichment would is the failure this whole path exists to avoid. No + two axes share a name in `_AXIS_COORDS`, so a name appears once. """ + out = [] + for axis in dist_map.axes: + for name in _AXIS_COORDS[axis]: + if name in names: + out.append((axis, name)) + break + return out + + +def _readable_on(dist_map) -> list[str]: + """The names a map could be read on, whichever axis answers.""" + return sorted({x for axis in dist_map.axes for x in _AXIS_COORDS[axis]}) + + +def _get_channel_axes(patch, acquisition) -> list[tuple[str, str]]: + """Return every map axis the patch can be read on, with its coordinate.""" dist_map = acquisition.distance_map if dist_map is None: msg = ( @@ -558,22 +578,15 @@ def _get_channel_axes(patch, acquisition) -> list[tuple[str, str]]: "its channels cannot be placed on the optical path." ) raise PatchError(msg) - out = [] - for axis in dist_map.axes: - for name in _AXIS_COORDS[axis]: - if name in patch.coords.coord_map: - out.append((axis, name)) - break - if out: + if out := _map_axis_coords(dist_map, patch.coords.coord_map): return out - wanted = sorted({x for axis in dist_map.axes for x in _AXIS_COORDS[axis]}) msg = ( f"Acquisition {acquisition.code!r} maps {list(dist_map.axes)} onto " - f"path distance, so it needs one of the {wanted} coordinates, and " - f"this patch has {sorted(patch.coords.coord_map)}. An acquisition " - "whose patches carry interrogator meters is calibrated with a " - "distance_map on the instrument_distance axis, one control point " - "being enough to state an origin." + f"path distance, so it needs one of the {_readable_on(dist_map)} " + f"coordinates, and this patch has {sorted(patch.coords.coord_map)}. " + "An acquisition whose patches carry interrogator meters is " + "calibrated with a distance_map on the instrument_distance axis, " + "one control point being enough to state an origin." ) raise PatchError(msg) @@ -795,27 +808,19 @@ def _channel_placement(dims: set[str], acquisition) -> tuple: f"{acquisition.code!r} defines no distance_map, so its channels " "cannot be placed on the optical path" ) - found = {} - for axis in dist_map.axes: - # One coordinate per axis, the same preference `_get_channel_axes` - # applies; two axes landing on one dimension is that dimension - # stating both, which the map's own validator has already agreed. - for name in _AXIS_COORDS[axis]: - if name in dims: - found.setdefault(name, axis) - break + found = _map_axis_coords(dist_map, dims) if not found: - wanted = sorted({x for axis in dist_map.axes for x in _AXIS_COORDS[axis]}) return None, ( f"has dimensions {sorted(dims)}, and {acquisition.code!r} places " - f"channels by one of {wanted}" + f"channels by one of {_readable_on(dist_map)}" ) - if len(found) > 1: + if len({name for _, name in found}) > 1: return None, ( - f"carries {sorted(found)} as separate dimensions, so which of " - f"them {acquisition.code!r} places its channels by is ambiguous" + f"carries {sorted(name for _, name in found)} as separate " + f"dimensions, so which of them {acquisition.code!r} places its " + "channels by is ambiguous" ) - name, axis = next(iter(found.items())) + axis, name = found[0] return name, axis @@ -893,13 +898,15 @@ def resolve_channel_pieces( Returns ------- A `(name, pieces, reasons)` triple. `reasons` holds a refusal per row - and None elsewhere; when any row is refused the other two are None, - since the caller raises rather than selecting. + and None elsewhere; when any row is refused `pieces` is None, since + the caller raises rather than selecting, and `name` is whatever the + rows which placed fine agreed on — possibly None. """ name, placements, reasons = channel_placements(contexts, frame) if any(x is not None for x in reasons) or name is None: - # No dimension to trim along, so nothing was judged and there are - # no pieces to report; the caller reads that off `name`. + # Nothing was judged -- either a row must be refused, or no row + # has a fiber at all -- so there are no pieces to report; the + # caller reads which of the two off `reasons` and `name`. return name, None, reasons out, unusable = [], [] for row in _placed_rows(contexts, placements, frame, name): @@ -989,6 +996,11 @@ def _placed_rows(contexts, placements, frame, name: str): reason = "states no channel spacing to place its channels on" yield PlacedRow(context, low, high, None, None, reason) continue + # Envelopes are value-ordered whatever the coordinate's orientation, + # so the grid is walked by the step's magnitude -- a reverse-sorted + # patch states a negative one, and counting samples with it would + # give none at all. + step = abs(step) key = (id(context), axis, low, high, step) if (placed := cache.get(key)) is None: grid = np.arange(round((high - low) / step) + 1) * step + low @@ -1004,9 +1016,10 @@ def resolve_split_pieces(inventory, contexts, frame, name, keep) -> tuple: Splitting expands the spool into one patch per value a group takes along each row, so a row is answered with a list of `(value, piece)` - pairs rather than pieces alone. A group whose values overlap gives - some channels to more than one output, which is the whole point of a - membership group and why the pieces of a row need not be disjoint. + pairs rather than pieces alone. A value the group states in two + places gives that value two pieces, so a row's pieces are disjoint + but neither contiguous nor in envelope order — they arrive grouped by + value, and the values are sorted rather than the envelopes. Parameters ---------- @@ -1036,12 +1049,11 @@ def resolve_split_pieces(inventory, contexts, frame, name, keep) -> tuple: if row.grid is None: continue path = row.context.optical_path - values = _get_coord_values(inventory, path, name, row.distances) + values, _ = _channel_values(inventory, path, name, row.distances) if values is None: # A group this path defines nowhere puts none of its channels # anywhere, so the row contributes no output at all. continue - values = values.values if isinstance(values, BaseCoord) else values for value in _split_values(values): if not keep(value): continue @@ -1067,24 +1079,48 @@ def _split_values(values) -> list: array = np.asarray(values) if array.dtype == bool: return sorted(set(array.tolist())) - keep = ~np.isnan(array) if np.issubdtype(array.dtype, np.number) else array != "" - return sorted(set(array[keep].tolist())) + return sorted(set(array[~_undefined_mask(array)].tolist())) + + +def _channel_values(inventory, path, name, distances) -> tuple: + """ + Return one name's value per channel, with the units it carries. + + A path is what states anything along the fiber, so an acquisition + without one -- a valid inventory, describing a system which simply + projects nothing -- has no values, exactly as a path which defines + the name nowhere has none. + """ + values = ( + None if path is None else _get_coord_values(inventory, path, name, distances) + ) + if values is None: + return None, None + if not isinstance(values, BaseCoord): + return values, None + # The projection carries the units the inventory documents for the + # field, and dropping them would refuse the unit-bearing selectors + # the index accepts against a stated attr of the same kind. + units = None if values.units is None else {"num": get_quantity_str(values.units)} + return values.values, units def _channel_matches(inventory, path, name, selector, distances) -> np.ndarray: """Return the channels one selector matches, as the index would.""" from dascore.io.index.query import evaluate_attr_predicate # noqa: PLC0415 - values = _get_coord_values(inventory, path, name, distances) + values, units = _channel_values(inventory, path, name, distances) if values is None: # A name this path defines nowhere states nothing about any of its # channels, so it matches none of them -- listing a name is not # promising a value for it. return np.zeros(len(distances), dtype=bool) - values = values.values if isinstance(values, BaseCoord) else values if selector is None: + # `None` is the query spelling of the undefined marker, which is a + # value here rather than the "select everything" a bare None means + # of an attr: a channel the track says nothing about is a channel. return _undefined_mask(values) - return evaluate_attr_predicate(list(values), name, selector) + return evaluate_attr_predicate(list(values), name, selector, units) def _coords_equal(existing, values) -> bool: diff --git a/dascore/utils/chunk_plan.py b/dascore/utils/chunk_plan.py index efa7025f..d23a5712 100644 --- a/dascore/utils/chunk_plan.py +++ b/dascore/utils/chunk_plan.py @@ -1145,9 +1145,10 @@ def build_subdivision_plan(df: pd.DataFrame, pieces, name: str) -> ChunkPlan: pieces One sequence of inclusive `(low, high)` envelopes per row, in the row's own units and already on its sample grid. They must be - disjoint and ascending within a row; `subdivision_pieces` builds - them from cut values, and a mask over the row's samples gives - them directly. + disjoint within a row but need not be in envelope order — + `Spool.split_by` emits them grouped by value. `subdivision_pieces` + builds them from cut values, and a mask over the row's samples + gives them directly. name The dimension being subdivided. """ diff --git a/docs/tutorial/patch.qmd b/docs/tutorial/patch.qmd index dc969689..c9f55ae2 100644 --- a/docs/tutorial/patch.qmd +++ b/docs/tutorial/patch.qmd @@ -462,7 +462,7 @@ assert patch.get_coord("distance").step is not None assert outside.get_coord("distance").step is None ``` -That is why [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses coordinates. A patch can have samples removed from its middle; at spool level the complement of a range would be a hole in every patch rather than a choice between patches. When several coordinates are named, each is complemented on its own — the true complement of a block is an L, which no array can hold. +That is why [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses coordinates. A patch can have samples removed from its middle; at spool level the complement of a range would be a hole in every patch rather than a choice between patches. When several coordinates are named, each is complemented on its own — the true complement of a block is a frame around it, which no array can hold. ## Order Order is similar to [`Patch.select`](`dascore.Patch.select`), but will re-arrange data to the order specified by a value array. This may also cause parts of the patch to be duplicated. diff --git a/tests/test_proc/test_proc_coords.py b/tests/test_proc/test_proc_coords.py index 4f26bccb..1514ed9d 100644 --- a/tests/test_proc/test_proc_coords.py +++ b/tests/test_proc/test_proc_coords.py @@ -704,8 +704,8 @@ def test_each_named_coordinate_is_complemented(self, random_patch): """ Two names remove two ranges rather than the one intersection. - The complement of a block is an L, which no array can hold, so - unselect removes the part of it which is expressible. + The complement of a block is a frame around it, which no array + can hold, so unselect removes the part which is expressible. """ time = random_patch.get_array("time") window = (time[0], time[4]) @@ -720,6 +720,19 @@ def test_unknown_coordinate_raises(self, random_patch): with pytest.raises(PatchCoordinateError, match="not found in patch"): random_patch.unselect(not_a_coord=(1, 2)) + def test_a_multidimensional_coordinate_raises(self, random_patch): + """ + A range of one names no samples of a single dimension to drop. + + Its complement is a shape spanning both, which is the same reason + two coordinates cannot be complemented jointly. + """ + size = random_patch.shape + grid = np.arange(size[0] * size[1]).reshape(size) + patch = random_patch.update_coords(quality=(("distance", "time"), grid)) + with pytest.raises(PatchCoordinateError, match="spans"): + patch.unselect(quality=(0, 10)) + def test_non_dimensional_coordinate(self, random_patch): """A coordinate along a dimension trims that dimension.""" size = random_patch.coord_shapes["distance"][0] diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index ff0a05df..014ca344 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -24,6 +24,7 @@ from dascore.core.inventory import ( Acquisition, CoordinateReferenceSystem, + CouplingCondition, DistanceMap, FiberArray, FiberSegment, @@ -46,6 +47,7 @@ UnresolvedPatchError, ) from dascore.proc.inventory import resolve_row_epochs +from dascore.units import cm, m @pytest.fixture(scope="module") @@ -2586,14 +2588,6 @@ def test_trims_to_the_matching_channels(self, patch, inventory): assert out[0].get_coord("distance").min() == 0 assert out[0].get_coord("distance").max() == 150 - def test_len_is_exact_before_anything_loads(self, patch, inventory): - """Selection is metadata work: the count is final without a read.""" - spool = dc.spool(patch).attach_inventory(inventory) - out = spool.select(coupling="trench") - contents = out.get_contents() - assert len(out) == len(contents) == 1 - assert contents["distance_max"].iloc[0] == 150.0 - def test_the_data_is_the_channels_it_names(self, patch, inventory): """ The rows kept are the rows the coordinate names. @@ -2625,7 +2619,9 @@ def test_a_disjoint_match_keeps_each_piece_whole(self, patch, two_zones): """Each piece holds exactly the rows its own envelope names.""" spool = dc.spool(patch).attach_inventory(two_zones) whole = patch.get_array("distance") - for piece in spool.select(hole="a"): + pieces = list(spool.select(hole="a")) + assert len(pieces) == 2 + for piece in pieces: rows = np.searchsorted(whole, piece.get_array("distance")) assert np.array_equal(piece.data, patch.data[rows]) @@ -2637,7 +2633,10 @@ def test_selection_agrees_with_enrichment(self, patch, two_zones): makes that worth doing rather than a coincidence to maintain. """ spool = dc.spool(patch).attach_inventory(two_zones) - for piece in spool.select(hole="a").enrich(): + pieces = list(spool.select(hole="a").enrich()) + assert len(pieces) == 2 + for piece in pieces: + assert len(piece.get_array("hole")) assert set(piece.get_array("hole")) == {"a"} def test_a_value_nothing_holds_keeps_nothing(self, patch, inventory): @@ -2687,6 +2686,7 @@ def test_composes_with_a_time_split(self, patch, inventory): split = _split_epochs(inventory, middle, second={"name": "moved"}) spool = dc.spool(patch).attach_inventory(split) out = spool.conform_to_inventory().select(coupling="trench") + assert len(out) == 2 times, dists = patch.get_array("time"), patch.get_array("distance") for piece in out: rows = np.searchsorted(dists, piece.get_array("distance")) @@ -2718,6 +2718,10 @@ def test_a_whole_match_changes_nothing(self, patch, inventory): """Keeping every channel needs no plan, so the spool is unchanged.""" spool = dc.spool(patch).attach_inventory(inventory) out = spool.select(geometry="trench") + # The same resolver, not merely an equal spool: building a plan + # whose single piece covered its row would load the same patch + # and pass every other assertion here. + assert out._catalog.resolver is spool._catalog.resolver assert len(out) == 1 assert out[0].shape == patch.shape @@ -2793,11 +2797,18 @@ def test_the_patch_axis_keeps_its_own_meaning(self, patch, inventory): out = spool.select(**form, **kwargs) assert out.get_contents()["distance_max"].tolist() == [100.0] - def test_samples_with_a_channel_name_raises(self, patch, inventory): - """A fiber coordinate has no sample numbering of its own.""" + @pytest.mark.parametrize("flag", ["samples", "relative"]) + def test_axis_keywords_with_a_channel_name_raise(self, patch, inventory, flag): + """ + Neither keyword has anything to say about a fiber coordinate. + + Both describe the patch's own axis, while these say what is + attached to each channel of it. Ignoring one quietly would answer + a different question than the caller asked. + """ spool = dc.spool(patch).attach_inventory(inventory) - with pytest.raises(InvalidSpoolQueryError, match="no sample numbering"): - spool.select(coupling="trench", samples=True) + with pytest.raises(InvalidSpoolQueryError, match=f"{flag}=True cannot"): + spool.select(coupling="trench", **{flag: True}) def test_a_channel_name_through_attrs_raises(self, patch, inventory): """It describes channels, so the attrs namespace is the wrong one.""" @@ -2828,6 +2839,10 @@ def test_complements_the_selection(self, patch, two_zones): spool = dc.spool(patch).attach_inventory(two_zones) kept = _channels(spool.select(hole="a")) dropped = _channels(spool.unselect(hole="a")) + # Stated rather than merely complementary: kept=nothing and + # dropped=everything satisfies a partition too, and is exactly + # what the complement plumbing could produce by mistake. + assert len(kept) == 82 # 10-50 and 200-240 inclusive assert not set(kept) & set(dropped) assert sorted([*kept, *dropped]) == sorted(patch.get_array("distance")) @@ -2884,7 +2899,9 @@ def test_the_pieces_hold_their_own_channels(self, patch, inventory): """Each output holds the rows its value covers, and no others.""" spool = dc.spool(patch).attach_inventory(inventory) whole = patch.get_array("distance") - for piece in spool.split_by("zone"): + pieces = list(spool.split_by("zone")) + assert len(pieces) == 2 + for piece in pieces: rows = np.searchsorted(whole, piece.get_array("distance")) assert np.array_equal(piece.data, patch.data[rows]) @@ -2897,14 +2914,23 @@ def test_stamp_false_leaves_the_attrs_alone(self, patch, inventory): """Which is what a nested split wants of the second one.""" spool = dc.spool(patch).attach_inventory(inventory) out = spool.split_by("zone", stamp=False) + assert len(out) == 2 # the split still happened; only the stamp is off assert "zone" not in out.get_contents().columns assert not any(dict(x.attrs).get("zone") for x in out) - def test_the_stamp_is_selectable(self, patch, inventory): - """A stamped value is an ordinary attr, so it filters like one.""" + def test_the_stamp_shadows_the_inventory_name(self, patch, inventory): + """ + A stamped value is an ordinary attr, and attrs win over the fiber. + + Both paths would keep one patch here, so the patch's shape is + what tells them apart: the attr path keeps the north piece whole, + while resolving `zone` along the fiber again would trim it. + """ spool = dc.spool(patch).attach_inventory(inventory) out = spool.split_by("zone") - assert len(out.select(zone="north")) == 1 + north = out.select(zone="north") + assert len(north) == 1 + assert north[0].shape == out[0].shape def test_a_membership_group_splits_in_two(self, patch, inventory): """Both sides come out: the channels included and those not.""" @@ -2912,18 +2938,25 @@ def test_a_membership_group_splits_in_two(self, patch, inventory): out = spool.split_by("noisy") assert sorted(out.get_contents()["noisy"].tolist()) == [False, False, True] - def test_overlapping_values_share_channels(self, patch, inventory): + def test_one_split_partitions_the_channels(self, patch, inventory): """ - A group may overlap another, so a channel can land in two outputs. + A channel holds one value of a group, so one split cannot share it. - `noisy` runs from 150 to 300 m and the zones meet at 200, so the - two split differently over the same fiber. + Overlapping intervals of a group resolve to a single value per + channel — the projection `Patch.enrich` uses — so the outputs of + one call divide the fiber rather than covering it twice. Two + *different* groups may still cut it differently, which is what + makes a nested split worth doing. """ spool = dc.spool(patch).attach_inventory(inventory) - zones = _channels(spool.split_by("zone")) - noisy = _channels(spool.split_by("noisy")) - assert sorted(zones) == sorted(patch.get_array("distance")) - assert sorted(noisy) == sorted(patch.get_array("distance")) + for group in ("zone", "noisy"): + channels = _channels(spool.split_by(group)) + assert sorted(channels) == sorted(patch.get_array("distance")) + assert len(channels) == len(set(channels.tolist())) + # The two groups do not agree about where the fiber divides. + zoned = {tuple(x.get_array("distance")) for x in spool.split_by("zone")} + noisy = {tuple(x.get_array("distance")) for x in spool.split_by("noisy")} + assert zoned != noisy def test_a_disjoint_value_becomes_several_patches(self, patch, two_zones): """One value covering two stretches keeps them apart.""" @@ -2944,10 +2977,17 @@ def test_exclude_wins_over_include(self, patch, inventory): out = spool.split_by("zone", include=("nor*", "sou*"), exclude="north") assert out.get_contents()["zone"].tolist() == ["south"] - def test_a_group_no_path_defines_yields_nothing(self, patch, inventory): - """There is no value to split on, so there is no output.""" + def test_a_name_the_inventory_lacks_raises(self, patch, inventory): + """ + A misspelling has no values to split into, so it says so. + + Returning an empty spool would be indistinguishable from a group + which happens to cover nothing, and selection refuses a name it + does not know for the same reason. + """ spool = dc.spool(patch).attach_inventory(inventory) - assert len(spool.split_by("not_a_group")) == 0 + with pytest.raises(InvalidSpoolQueryError, match="not a coordinate"): + spool.split_by("not_a_group") def test_needs_an_inventory(self, patch): """The values it expands into are ones an inventory states.""" @@ -3027,7 +3067,11 @@ def test_patches_on_different_channel_dimensions_refuse(self, patch, inventory): def test_none_on_a_membership_group(self, patch, inventory): """A membership group says something about every channel: False.""" spool = dc.spool(patch).attach_inventory(inventory) + # noisy runs from 150 to 300 m, which is patch distance 50 to 200, + # so the channels it says nothing about are those on either side. undefined = spool.select(noisy=None).get_contents() + assert undefined["distance_min"].tolist() == [0.0, 201.0] + assert undefined["distance_max"].tolist() == [49.0, 299.0] assert undefined["distance_min"].tolist() == ( spool.select(noisy=False).get_contents()["distance_min"].tolist() ) @@ -3115,3 +3159,224 @@ def test_splitting_a_lag_time_patch_yields_nothing(self, patch, inventory): lags = patch.get_coord("time").values - patch.get_coord("time").min() spool = dc.spool(patch.update_coords(time=lags)).attach_inventory(inventory) assert len(spool.split_by("zone")) == 0 + + +@pytest.fixture(scope="module") +def uneven_spool(patch, inventory): + """ + Three patches whose fiber answers differ, in one spool. + + The first resolves to a path putting one group in two stretches, the + second to a path putting it in one, and the third to nothing at all. + """ + array = inventory.networks[0].fiber_arrays[0] + acquisition, path = array.acquisitions[0], array.optical_paths[0] + holes = ( + OpticalPathAnnotation( + start_distance=110.0, end_distance=150.0, group="hole", value="a" + ), + OpticalPathAnnotation( + start_distance=300.0, end_distance=340.0, group="hole", value="a" + ), + ) + both = inventory.replace( + array, + array.new( + acquisitions=(acquisition, acquisition.new(location_code="01")), + optical_paths=( + path.new(annotations=(*path.annotations, *holes)), + path.new(location_code="01", annotations=(*path.annotations, holes[0])), + ), + ), + ) + patches = [ + patch.update_attrs(tag="split"), + patch.update_attrs(acquisition_key="DAS.R2D1.01.RAW", tag="solid"), + patch.update_attrs(acquisition_key="DAS.R2D1..NOPE", tag="unknown"), + ] + return dc.spool(patches).attach_inventory(both), patches + + +class TestChannelSelectAlignment: + """ + Rows subdividing by different amounts must not swap answers. + + One patch becoming two while its neighbour becomes one and a third + leaves entirely is where a plan built by position rather than by + identity would quietly hand a patch someone else's envelope. + """ + + def test_each_patch_keeps_its_own_pieces(self, uneven_spool): + """The split patch gets two, the solid one gets one, the third none.""" + spool, _ = uneven_spool + contents = spool.select(hole="a").get_contents() + by_tag = contents.groupby("tag")["distance_min"].apply(sorted) + assert by_tag["split"] == [10.0, 200.0] + assert by_tag["solid"] == [10.0] + assert "unknown" not in by_tag + + def test_each_piece_holds_its_own_patch_data(self, uneven_spool): + """And the array behind each envelope is that patch's, not a neighbour's.""" + spool, patches = uneven_spool + sources = {x.attrs.tag: x for x in patches} + whole = patches[0].get_array("distance") + for piece in spool.select(hole="a"): + source = sources[piece.attrs.tag] + rows = np.searchsorted(whole, piece.get_array("distance")) + assert np.array_equal(piece.data, source.data[rows]) + + def test_the_complement_is_exact_for_every_patch(self, uneven_spool): + """Including the one the inventory says nothing about, kept whole.""" + spool, patches = uneven_spool + whole = sorted(patches[0].get_array("distance")) + kept, dropped = {}, {} + for store, out in ( + (kept, spool.select(hole="a")), + (dropped, spool.unselect(hole="a")), + ): + for piece in out: + store.setdefault(piece.attrs.tag, []).extend( + piece.get_array("distance") + ) + for tag in ("split", "solid"): + assert not set(kept[tag]) & set(dropped[tag]) + assert sorted([*kept[tag], *dropped[tag]]) == whole + assert sorted(dropped["unknown"]) == whole + + def test_splitting_keeps_each_patch_with_its_value(self, uneven_spool): + """split_by builds the same plan, so it aligns the same way.""" + spool, patches = uneven_spool + sources = {x.attrs.tag: x for x in patches} + whole = patches[0].get_array("distance") + out = spool.split_by("hole") + assert out.get_contents()["hole"].tolist() == ["a", "a", "a"] + for piece in out: + source = sources[piece.attrs.tag] + rows = np.searchsorted(whole, piece.get_array("distance")) + assert np.array_equal(piece.data, source.data[rows]) + + +class TestChannelReviewFindings: + """Defects the review pipeline found, each with its own scenario.""" + + def test_a_reverse_sorted_patch_is_still_placed(self, patch, inventory): + """ + A descending coordinate states a negative step, and envelopes + do not: counting samples with the signed one gives none at all, + so every fiber query would silently drop the patch. + """ + reversed_patch = patch.sort_coords("distance", reverse=True) + spool = dc.spool(reversed_patch).attach_inventory(inventory) + assert spool.get_contents()["distance_step"].iloc[0] < 0 + out = spool.select(coupling="trench") + assert len(out) == 1 + assert out.get_contents()["distance_max"].tolist() == [150.0] + assert out[0].shape == (151, patch.shape[1]) + + def test_an_acquisition_without_a_path_states_nothing(self, patch, inventory): + """ + A pathless acquisition is valid, and simply projects nothing. + + The CRS labels stay listed whatever the paths say, so a geometry + axis is still a name a query may use — and must answer with no + channel rather than by dereferencing the absent path. + """ + array = inventory.networks[0].fiber_arrays[0] + pathless = inventory.replace(array, array.new(optical_paths=())) + spool = dc.spool(patch).attach_inventory(pathless) + assert "x" in pathless.get_names().coords + assert len(spool.select(x=(-1e9, 1e9))) == 0 + assert len(spool.split_by("x")) == 0 + + def test_a_unit_bearing_selector_is_converted(self, patch, inventory): + """ + The inventory documents the units of the fields it projects, so + dropping them would refuse selectors the index accepts against a + stated attr of the same kind. + """ + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + deep = inventory.replace( + path, + path.new( + coupling=( + CouplingCondition( + start_distance=100.0, + end_distance=250.0, + coupling_type="trench", + medium="soil", + depth=2.0, + ), + ) + ), + ) + spool = dc.spool(patch).attach_inventory(deep) + bare = spool.select(**{"coupling.depth": (1.0, 3.0)}) + assert bare.get_contents()["distance_max"].tolist() == [150.0] + for equivalent in ((1 * m, 3 * m), (100 * cm, 300 * cm)): + out = spool.select(**{"coupling.depth": equivalent}) + assert out.get_contents()["distance_max"].tolist() == [150.0] + + def test_a_bare_ellipsis_selects_everything(self, patch, inventory): + """ + As it does everywhere else, so attaching an inventory does not + turn a no-op selector into an error. `None` is the exception, + and means the undefined marker rather than "everything". + """ + spool = dc.spool(patch).attach_inventory(inventory) + assert len(spool.select(coupling=...)) == 1 + assert spool.select(coupling=...)[0].shape == patch.shape + + def test_unselect_takes_a_channel_name_through_coords(self, patch, inventory): + """ + Both spellings `select` accepts, since `_coords` is the only way + to name a group which shares an acquisition field's name. + """ + spool = dc.spool(patch).attach_inventory(inventory) + expected = [151.0] + assert ( + spool.unselect(coupling="trench").get_contents()["distance_min"].tolist() + == expected + ) + assert ( + spool.unselect(_coords={"coupling": "trench"}) + .get_contents()["distance_min"] + .tolist() + == expected + ) + assert ( + spool.unselect(_coords="coupling", coupling="trench") + .get_contents()["distance_min"] + .tolist() + == expected + ) + + def test_rechunking_a_selection_keeps_it(self, patch, inventory): + """ + Re-planning the same dimension collapses onto the sources, which + is sound only while the plan's pieces cover them. A selection's + do not, so collapsing would load back the channels it removed — + and the contents would go on describing the ones it kept. + """ + spool = dc.spool(patch).attach_inventory(inventory) + selected = spool.select(coupling="trench") + assert selected[0].shape == (151, patch.shape[1]) + rechunked = selected.chunk(distance=None) + assert rechunked.get_contents()["distance_max"].tolist() == [150.0] + assert rechunked[0].shape == (151, patch.shape[1]) + + def test_only_a_plan_which_drops_samples_is_lossy( + self, patch, inventory, path_epochs + ): + """ + The other side of it: conform's pieces do cover their row, so it + stays collapsible and only a selection does not. + + Checked on the plans themselves rather than by re-chunking each, + because chunking a derived catalog is broken on dev already + (#871) and would fail here for a reason of its own. + """ + conformed = dc.spool(patch).conform_to_inventory(path_epochs) + assert len(conformed) == 2 + assert not conformed._catalog.resolver.lossy + selected = dc.spool(patch).attach_inventory(inventory).select(coupling="trench") + assert selected._catalog.resolver.lossy diff --git a/tests/test_utils/test_chunk.py b/tests/test_utils/test_chunk.py index e5c867ff..2a52d712 100644 --- a/tests/test_utils/test_chunk.py +++ b/tests/test_utils/test_chunk.py @@ -585,11 +585,12 @@ def test_an_on_grid_cut_lands_on_its_own_sample(self, index): """ A cut sitting exactly on a sample opens the piece at that sample. - Dividing the native types keeps this exact for a time - coordinate — datetime64 arithmetic is integer nanoseconds — where - converting each operand to float seconds first rounds three - times and puts the boundary sample one piece too early at these - step and index combinations. + Dividing the native types rounds once — datetime64 arithmetic + is integer nanoseconds, so only the quotient rounds — where + converting each operand to float seconds first rounds three times + and puts the boundary sample one piece too early at these step + and index combinations. The grid comparison settles it either + way, which is what these pin. """ start, step = STARTTIME, np.timedelta64(17060962, "ns") df = _one_row_df(start, step, index + 10) @@ -618,10 +619,10 @@ def test_an_on_grid_float_cut_is_not_pushed_past_its_sample(self): `1.0 + 0.1` is the textbook case: the difference back out is `0.10000000000000009`, so the ratio is a hair above 1 and `ceil` - answers 2 — the sample the cut lands exactly on would open the - piece *before* it. Comparing against the grid is what settles - it, and a distance axis makes this the ordinary case rather than - the exotic one. + answers 2 — the sample the cut lands exactly on would be left in + the piece *before* the cut instead of opening the one after it. + Comparing against the grid is what settles it, and a distance + axis makes this the ordinary case rather than the exotic one. """ start, step = 1.0, 0.1 df = _one_row_df(start, step, 10, name="distance") @@ -693,3 +694,48 @@ def test_one_piece_sequence_per_row(self): df = _one_row_df(0.0, 1.0, 10, name="distance") with pytest.raises(AssertionError): build_subdivision_plan(pd.concat([df, df]), [[]], "distance") + + +class TestSnappedCutExactness: + """The property the two correction loops exist to guarantee.""" + + @pytest.mark.parametrize("seed", [0, 1, 2]) + def test_an_on_grid_cut_always_lands_on_its_own_sample(self, seed): + """ + Over many float grids, a cut on a sample opens the piece there. + + The ratio alone gets this wrong on a few percent of inputs in + either direction, and which few depends on the values, so a + handful of hand-picked cases cannot stand for it. The pieces are + read straight off `subdivision_pieces` rather than a plan so the + grid arithmetic is what is under test. + """ + rng = np.random.default_rng(seed) + for _ in range(300): + start = float(rng.uniform(-1e4, 1e4)) + step = float(rng.uniform(1e-6, 10)) + index = int(rng.integers(1, 200_000)) + cut = start + index * step + df = _one_row_df(start, step, index + 2, name="distance") + pieces = subdivision_pieces(df, [(cut,)], "distance")[0] + assert len(pieces) == 2 + # The second piece opens exactly on the cut, and the first + # ends one step short of it, whichever way the ratio erred. + assert pieces[1][0] == cut + assert pieces[0][1] == cut - step + + def test_a_cut_between_samples_opens_at_the_next(self): + """The other direction, over the same spread of grids.""" + rng = np.random.default_rng(3) + for _ in range(300): + start = float(rng.uniform(-1e4, 1e4)) + step = float(rng.uniform(1e-6, 10)) + index = int(rng.integers(1, 200_000)) + on_grid = start + index * step + cut = np.nextafter(on_grid, np.inf) + df = _one_row_df(start, step, index + 2, name="distance") + pieces = subdivision_pieces(df, [(cut,)], "distance")[0] + # A cut past a sample leaves that sample behind, so the piece + # it opens starts at the next one. + assert pieces[1][0] > on_grid + assert pieces[1][0] >= cut From fa28959e4d9dd8e29c9561d5dec0216022c84808 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 23:44:18 +0200 Subject: [PATCH 7/8] Address the PR bots Codex found one which corrupts a catalog: an annotation group may be named anything the inventory does not reserve, and the stamp is assigned straight onto the plan's outputs, so a group called output_id replaced the column binding each output to the data it came from. split_by now refuses to stamp over the spool's own bookkeeping, and stamp=False still splits on such a group. It also found that a bare ... vetoed samples and relative for the whole query, though it asks nothing of the coordinate it names, and that scalar equality validated a unit conversion without keeping it -- so a geometry axis stored in degrees never matched a selector pint bases in radians, where the range form did. That last one is in the index's own predicate and applies equally to a stated attr; it was reachable before this branch and is much more reachable now. CodeRabbit found three pieces of prose the branch made false: attach_inventory still said fiber coordinates were not selectable, and Patch.unselect described itself as an exact complement without the qualification its own note goes on to make. --- dascore/core/spool.py | 43 +++++++++++++++++--- dascore/io/index/query.py | 15 +++++-- dascore/proc/coords.py | 14 ++++--- docs/tutorial/patch.qmd | 4 +- tests/test_proc/test_proc_inventory.py | 55 +++++++++++++++++++++++++- 5 files changed, 114 insertions(+), 17 deletions(-) diff --git a/dascore/core/spool.py b/dascore/core/spool.py index 17194e2d..d583c326 100644 --- a/dascore/core/spool.py +++ b/dascore/core/spool.py @@ -269,6 +269,31 @@ def _drops_samples(rows: pd.DataFrame, pieces, name: str) -> bool: return False +def _check_stampable(name: str, rows: pd.DataFrame) -> None: + """ + Refuse a stamp which would overwrite the plan's own bookkeeping. + + An annotation group may be named anything the inventory does not + reserve, and the stamp is assigned onto the outputs — so a group + called `output_id` would replace the column binding each output to + its members, and one called `time_min` an envelope. Overwriting a + carried attr is fine and is how re-splitting restamps; these are not + attrs. + """ + envelopes = { + x for x in rows.columns if x.rsplit("_", 1)[-1] in {"min", "max", "step"} + } + if name not in {"output_id", "dims", *envelopes} and not name.startswith("_"): + return + msg = ( + f"{name!r} is how the spool itself describes a patch, so stamping " + "it would overwrite what binds each output to the data it came " + "from. Rename the group, or pass stamp=False to split on it " + "without recording the value." + ) + raise InvalidSpoolQueryError(msg) + + def _stated_channels(channels: dict) -> dict: """ Return the channel selectors which actually select something. @@ -1175,12 +1200,16 @@ def _split_inventory_query(self, _attrs, _coords, kwargs, samples, relative=Fals # describes are the patch's -- and no endpoints to be relative to, # since what it says varies from one acquisition to the next. # Ignoring either quietly would answer a question nobody asked. + # Judged against the selectors which actually select something: a + # bare `...` names a fiber coordinate without asking anything of + # it, so it must not veto a flag the rest of the query needs. + stated_channels = _stated_channels(channels) for flag, label in ((samples, "samples"), (relative, "relative")): - if not (channels and flag): + if not (stated_channels and flag): continue msg = ( - f"{sorted(channels)} name coordinates the inventory defines " - f"along the fiber, which {label}=True cannot describe: it " + f"{sorted(stated_channels)} name coordinates the inventory " + f"defines along the fiber, which {label}=True cannot describe: it " "asks about the patch's own axis, and these say what is " "attached to each channel of it." ) @@ -1389,8 +1418,10 @@ def attach_inventory(self, inventory) -> Self: Attaching promises nothing about the spool matching the inventory; [`conform_to_inventory`](`dascore.core.spool.Spool.conform_to_inventory`) - is what makes it so. Selecting on the coordinates an inventory - defines along the fiber is not implemented yet. + is what makes it so. Once attached, the coordinates the + inventory defines along the fiber become selectable, and + [`split_by`](`dascore.core.spool.Spool.split_by`) can expand the + spool by the values of one. """ if not isinstance(inventory, Inventory): msg = f"attach_inventory needs an Inventory, got {type(inventory)}." @@ -1589,6 +1620,8 @@ def split_by( ) raise InvalidSpoolQueryError(msg) source_rows, working = self._plan_frames() + if stamp: + _check_stampable(name, working) contexts = self._plan_contexts(working) dim, rows, reasons = resolve_split_pieces( self._inventory, contexts, working, name, _glob_filter(include, exclude) diff --git a/dascore/io/index/query.py b/dascore/io/index/query.py index e3b9e404..eb4d936c 100644 --- a/dascore/io/index/query.py +++ b/dascore/io/index/query.py @@ -14,7 +14,7 @@ import operator import re from collections.abc import Sequence -from dataclasses import dataclass, field +from dataclasses import dataclass, field, replace from enum import Enum, auto import numpy as np @@ -435,10 +435,17 @@ def evaluate_attr_predicate(values, name: str, value, units=None) -> np.ndarray: units = units or {} def selector(item): - """Type one selector value the way the index would.""" + """ + Type one selector value the way the index would, in its units. + + The converted value is what the comparison needs, not just the + check that it converts: a stored coordinate in degrees against a + selector pint bases in radians would otherwise match nothing, + where the range form — which keeps its converted bounds — does. + """ out = _coerce_scalar(item, kinds) - _to_target_unit(out, units.get(out.kind), name) - return out + value = _to_target_unit(out, units.get(out.kind), name) + return replace(out, value=value) def each(func) -> np.ndarray: return np.array([x is not None and bool(func(x)) for x in typed], dtype=bool) diff --git a/dascore/proc/coords.py b/dascore/proc/coords.py index ff4af15e..7b97bc0f 100644 --- a/dascore/proc/coords.py +++ b/dascore/proc/coords.py @@ -577,8 +577,9 @@ def unselect( Return the patch outside a selection. The complement of [`Patch.select`](`dascore.Patch.select`): it takes - the same selectors and removes exactly the samples that selection - would have kept. + the same selectors and removes the samples that selection would have + kept. With one coordinate named that is exactly the complement; with + several, each is complemented on its own — see the note below. Parameters ---------- @@ -612,9 +613,12 @@ def unselect( - Removing a range from the middle of a coordinate leaves a hole in it, so the result is no longer evenly sampled and the coordinate becomes a monotonic array. That is exactly what - [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses - coordinates for: at spool level the complement of a range is a hole - in every patch rather than a choice between patches. + [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses the + patches' *own* coordinates for: at spool level the complement of a + range is a hole in every patch rather than a choice between + patches. The coordinates an attached inventory defines along the + fiber it does accept, since removing one of those chooses which + channels a patch holds. - Each named coordinate is complemented on its own. Selecting on two coordinates keeps the samples in both ranges, and everything diff --git a/docs/tutorial/patch.qmd b/docs/tutorial/patch.qmd index c9f55ae2..ccafde17 100644 --- a/docs/tutorial/patch.qmd +++ b/docs/tutorial/patch.qmd @@ -440,7 +440,7 @@ assert len(sub_patch.get_array('distance')) == 4 ## Unselect -[`Patch.unselect`](`dascore.Patch.unselect`) is the complement of `select`: it takes the same selectors and removes exactly the samples `select` would have kept. +[`Patch.unselect`](`dascore.Patch.unselect`) is the complement of `select`: it takes the same selectors and removes the samples `select` would have kept. Naming one coordinate gives exactly the complement; naming several complements each on its own, as below. ```{python} import dascore as dc @@ -462,7 +462,7 @@ assert patch.get_coord("distance").step is not None assert outside.get_coord("distance").step is None ``` -That is why [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses coordinates. A patch can have samples removed from its middle; at spool level the complement of a range would be a hole in every patch rather than a choice between patches. When several coordinates are named, each is complemented on its own — the true complement of a block is a frame around it, which no array can hold. +That is why [`Spool.unselect`](`dascore.core.spool.Spool.unselect`) refuses the patches' own coordinates. A patch can have samples removed from its middle; at spool level the complement of a range would be a hole in every patch rather than a choice between patches. The coordinates an attached DASDAE inventory defines along the fiber are a separate case, and are accepted: removing one of those chooses which channels a patch holds. When several coordinates are named, each is complemented on its own — the true complement of a block is a frame around it, which no array can hold. ## Order Order is similar to [`Patch.select`](`dascore.Patch.select`), but will re-arrange data to the order specified by a value array. This may also cause parts of the patch to be duplicated. diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index 014ca344..e15ae969 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -47,7 +47,7 @@ UnresolvedPatchError, ) from dascore.proc.inventory import resolve_row_epochs -from dascore.units import cm, m +from dascore.units import cm, get_quantity, m @pytest.fixture(scope="module") @@ -3380,3 +3380,56 @@ def test_only_a_plan_which_drops_samples_is_lossy( assert not conformed._catalog.resolver.lossy selected = dc.spool(patch).attach_inventory(inventory).select(coupling="trench") assert selected._catalog.resolver.lossy + + def test_a_stamp_cannot_overwrite_the_plan(self, patch, inventory): + """ + A group may be named anything the inventory does not reserve, + and the stamp is assigned onto the plan's outputs — so a group + called `output_id` would replace what binds each output to the + data it came from, and the catalog would be built from nonsense. + """ + path = inventory.networks[0].fiber_arrays[0].optical_paths[0] + clash = inventory.replace( + path, + path.new( + annotations=( + OpticalPathAnnotation( + start_distance=100.0, + end_distance=200.0, + group="output_id", + value="a", + ), + ) + ), + ) + spool = dc.spool(patch).attach_inventory(clash) + with pytest.raises(InvalidSpoolQueryError, match="overwrite"): + spool.split_by("output_id") + # Splitting on it is still legal; only recording the value is not. + assert len(spool.split_by("output_id", stamp=False)) == 1 + + def test_a_no_op_channel_selector_does_not_veto_a_flag(self, patch, inventory): + """ + `...` names a fiber coordinate without asking anything of it, so + it must not refuse a `samples` the rest of the query needs. + """ + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.select(distance=(0, 10), coupling=..., samples=True) + assert len(out) == 1 + assert out[0].shape == (10, patch.shape[1]) + # A selector which does ask something still refuses it. + with pytest.raises(InvalidSpoolQueryError, match="samples=True cannot"): + spool.select(coupling="trench", samples=True) + + def test_a_scalar_quantity_selector_is_converted(self, patch, inventory): + """ + A geometry axis is stored in the CRS's own units while a quantity + selector is typed in its base unit, so comparing the magnitudes + unconverted matched nothing — where the range form matched. + """ + spool = dc.spool(patch).attach_inventory(inventory) + enriched = patch.enrich(inventory, coords=("y",), attrs=False) + stated = float(enriched.get_array("y")[0]) + degrees = get_quantity("degree") + assert len(spool.select(y=stated)) == 1 + assert len(spool.select(y=stated * degrees)) == 1 From 8ea1f93522de667f8ae613ab3130270b1ba8c036 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Wed, 12 Aug 2026 23:49:48 +0200 Subject: [PATCH 8/8] Pin that a rebuilt grid names the channels the patch's own does CodeRabbit reasoned that since the grid is reconstructed from the index envelope rather than read off the patch, and Patch.select does not snap an interior bound to the nearest sample, a float grid which drifted could trim one channel too many. It does not happen over six pathological grids -- the count and the membership agree with what enrichment projects every time -- but the property is worth holding onto, so it is a test rather than a reply. Compared by count and position, not by value: a trimmed CoordRange regenerates its values from the piece's own start, so a plain select(distance=(a, b)) with no inventory anywhere already differs from the original in the last ulp. --- tests/test_proc/test_proc_inventory.py | 86 +++++++++++++++++++++++++- 1 file changed, 85 insertions(+), 1 deletion(-) diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index e15ae969..fb58b787 100644 --- a/tests/test_proc/test_proc_inventory.py +++ b/tests/test_proc/test_proc_inventory.py @@ -35,7 +35,7 @@ OpticalPath, OpticalPathAnnotation, ) -from dascore.examples import inventory_patch_pair +from dascore.examples import get_example_patch, inventory_patch_pair from dascore.exceptions import ( CoordMergeError, InvalidInventoryError, @@ -3256,6 +3256,57 @@ def test_splitting_keeps_each_patch_with_its_value(self, uneven_spool): assert np.array_equal(piece.data, source.data[rows]) +def _float_grid_pair(distance, span): + """A patch on a given float grid, and an inventory annotating its middle.""" + time = get_example_patch().get_array("time") + data = np.random.default_rng(0).random((len(distance), len(time))) + patch = dc.Patch( + data=data, + coords={"distance": distance, "time": time}, + dims=("distance", "time"), + attrs={"acquisition_key": "DAS.R2D1..RAW", "category": "DAS"}, + ) + acquisition = Acquisition( + code="RAW", + location_code="", + data_type="velocity", + data_category="DAS", + gauge_length=10.0, + sample_rate=1.0 / dc.to_float(patch.get_coord("time").step), + distance_map=DistanceMap( + instrument_distance=(float(distance.min()), float(distance.max())), + distance=(0.0, span), + ), + ) + path = OpticalPath( + name="main", + location_code="", + optical_components=(FiberSegment(name="c", optical_length=span + 10),), + annotations=( + OpticalPathAnnotation( + start_distance=span * 0.23, + end_distance=span * 0.61, + group="zone", + value="mid", + ), + ), + ) + return patch, Inventory( + networks=( + Network( + code="DAS", + fiber_arrays=( + FiberArray( + code="R2D1", + acquisitions=(acquisition,), + optical_paths=(path,), + ), + ), + ), + ) + ).check() + + class TestChannelReviewFindings: """Defects the review pipeline found, each with its own scenario.""" @@ -3433,3 +3484,36 @@ def test_a_scalar_quantity_selector_is_converted(self, patch, inventory): degrees = get_quantity("degree") assert len(spool.select(y=stated)) == 1 assert len(spool.select(y=stated * degrees)) == 1 + + @pytest.mark.parametrize( + ("start", "step", "size"), + [ + (1.0, 0.1, 300), + (0.1, 0.25, 500), + (-1234.5678, 0.037, 400), + (0.0, 1 / 3, 350), + ], + ) + def test_the_index_selects_what_enrichment_projects(self, start, step, size): + """ + A rebuilt grid must name the same channels the patch's own does. + + The grid is reconstructed from the index envelope rather than + read off the patch, and `Patch.select` does not snap an interior + bound to the nearest sample — so a float grid which drifted from + the patch's own values could trim one channel too many. Compared + by count and position rather than by value, since a trimmed + CoordRange regenerates its values from the piece's own start and + so differs in the last ulp for any trim, inventory or not. + """ + distance = start + np.arange(size) * step + span = float(distance.max() - distance.min()) + patch, inventory = _float_grid_pair(distance, span) + spool = dc.spool(patch).attach_inventory(inventory) + selected = spool.select(zone="mid") + got = np.concatenate([x.get_array("distance") for x in selected]) + projected = patch.enrich(inventory, coords=("zone",), attrs=False) + wanted = distance[projected.get_array("zone") == "mid"] + assert len(wanted) # the group really does cover part of this fiber + assert len(got) == len(wanted) + assert np.allclose(np.sort(got), np.sort(wanted), rtol=0, atol=1e-9)