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/core/spool.py b/dascore/core/spool.py index 827e692d..d583c326 100644 --- a/dascore/core/spool.py +++ b/dascore/core/spool.py @@ -6,7 +6,9 @@ import inspect import warnings 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 @@ -49,12 +51,14 @@ 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 from dascore.utils.misc import ( _spool_map, deep_equality_check, + iterate, ) from dascore.utils.namespace import NamespaceOwner from dascore.utils.patch import ( @@ -95,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' " @@ -154,56 +165,55 @@ 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: + """ + 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: """ - Refuse the patches whose acquisition changes partway through. + 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: @@ -232,6 +242,125 @@ 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 _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. + + 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. + + 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_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. + + 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 _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] + return kept or None + + def _namespace_names(spec: namespace_select_type) -> set[str]: """ Return the names an `_attrs`/`_coords` argument designates. @@ -607,12 +736,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, @@ -626,7 +760,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. @@ -919,8 +1055,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, relative) ) catalog = self._catalog.select( _attrs=_attrs, @@ -930,8 +1066,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 := _stated_channels(channel_query): + out = out._select_channels(channel_query) return out @compose_docstring(doc=get_docstring(BaseSpool.unselect)) @@ -945,36 +1083,35 @@ 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( known_attrs | selectable, known_coords, _attrs=_attrs, - _coords=_coords, - kwargs=kwargs, + _coords=_without_names(_coords, channels), + kwargs=_without_keys(kwargs, 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 _stated_channels(channels): msg = ( "unselect needs something to remove; " f"{sorted(requested) or 'nothing'} names no selection. " @@ -984,33 +1121,65 @@ 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( + _stated_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. + 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. """ - candidates = (requested - known) | (wanted_coords & coords) - if channel_level := sorted(candidates & coords): + 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) + # 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. @@ -1019,16 +1188,34 @@ 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 ) + # 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. + # 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 (stated_channels and flag): + continue + msg = ( + 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." + ) + raise InvalidSpoolQueryError(msg) + 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 # quietly move one out of the namespace it has always been in. @@ -1036,7 +1223,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, @@ -1052,7 +1239,64 @@ 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 resolve_channel_pieces # noqa: PLC0415 + + source_rows, working = self._plan_frames() + if not len(working): + return self + 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[~judged] = None + name, pieces, reasons = resolve_channel_pieces( + self._inventory, contexts, working, query, complement=complement + ) + _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. + 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: """ @@ -1174,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)}." @@ -1297,6 +1543,107 @@ 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. 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*, 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 + 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) + # 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() + 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) + ) + _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] + 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, @@ -1364,7 +1711,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 +1733,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 +1747,50 @@ 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, 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. `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 + + 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(kept, cuts, "time"), - parent=new._catalog, + plan=plan, + parent=self._catalog, merge_kwargs={}, mode="chunk", - origin_path=new.spool_path, + origin_path=self.spool_path, + stamped=stamped, + lossy=_drops_samples(rows, pieces, name), ) - 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/io/index/planned.py b/dascore/io/index/planned.py index 52ebfb74..0d86dd5d 100644 --- a/dascore/io/index/planned.py +++ b/dascore/io/index/planned.py @@ -367,6 +367,8 @@ def __init__( mode: str = "chunk", 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." @@ -382,6 +384,13 @@ 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) + # 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).""" @@ -460,20 +469,37 @@ 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()) - 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 out[0] - joined = members.assign(current_index=output_id) - patches = self._assembler()._patch_from_instruction_df(joined) - assert len(patches) == 1 - return patches[0] + 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: + """ + 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 +532,8 @@ def derived_catalog( mode: str = "chunk", check_behavior: WARN_LEVELS = "warn", origin_path=None, + stamped: tuple[str, ...] = (), + lossy: bool = False, ) -> PatchCatalog: """ Materialize a plan into a fresh in-memory catalog. @@ -574,6 +602,8 @@ def derived_catalog( mode=mode, 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() @@ -607,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/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 4cbb2573..7b97bc0f 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, @@ -566,6 +569,92 @@ 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 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 + ---------- + 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 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 + 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) + 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) def order( patch: PatchType, *, copy=False, relative=False, samples=False, **kwargs diff --git a/dascore/proc/inventory.py b/dascore/proc/inventory.py index 2cd8ccae..7c26654b 100644 --- a/dascore/proc/inventory.py +++ b/dascore/proc/inventory.py @@ -30,6 +30,7 @@ ) from dascore.exceptions import ( InvalidInventoryError, + InvalidSpoolQueryError, ParameterError, PatchError, UnresolvedPatchError, @@ -542,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 = ( @@ -557,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) @@ -770,6 +784,345 @@ 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 = _map_axis_coords(dist_map, dims) + if not found: + return None, ( + f"has dimensions {sorted(dims)}, and {acquisition.code!r} places " + f"channels by one of {_readable_on(dist_map)}" + ) + if len({name for _, name in found}) > 1: + return None, ( + 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" + ) + axis, name = found[0] + 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 `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: + # 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): + 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) + 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. + 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)) + msg = ( + "The patches of this spool place their channels along different " + f"dimensions ({joined}), so an operation along the fiber has no " + "one dimension to work on. Select them apart first." + ) + raise InvalidSpoolQueryError(msg) + return next(iter(named), None), placements, reasons + + +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 _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, tuple] = {} + for context, (dim, axis), (low, high), step in zip( + contexts, placements, bounds, steps, strict=True + ): + if context is None or dim is None: + yield PlacedRow(context, low, high, None, None, None) + continue + if pd.isnull(step) or not step: + # 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 + # 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 + placed = (grid, context.acquisition.channel_to_distance(grid, axis=axis)) + cache[key] = placed + grid, distances = placed + yield PlacedRow(context, low, high, grid, distances, 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 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 + ---------- + 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, _ = _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 + 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. + + 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 array.dtype == bool: + return sorted(set(array.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, 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) + 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, units) + + 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/dascore/utils/chunk_plan.py b/dascore/utils/chunk_plan.py index 92936a74..d23a5712 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,68 @@ 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 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. + """ + 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/docs/tutorial/patch.qmd b/docs/tutorial/patch.qmd index cf9eea52..ccafde17 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 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 + +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 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_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_coords.py b/tests/test_proc/test_proc_coords.py index deff814e..1514ed9d 100644 --- a/tests/test_proc/test_proc_coords.py +++ b/tests/test_proc/test_proc_coords.py @@ -644,6 +644,103 @@ 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 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]) + 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_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] + 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.""" diff --git a/tests/test_proc/test_proc_inventory.py b/tests/test_proc/test_proc_inventory.py index 3562ab21..fb58b787 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, @@ -34,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, @@ -46,6 +47,7 @@ UnresolvedPatchError, ) from dascore.proc.inventory import resolve_row_epochs +from dascore.units import cm, get_quantity, m @pytest.fixture(scope="module") @@ -722,11 +724,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 +1458,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 +1564,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 +1873,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 +1895,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): """ @@ -2358,7 +2375,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) @@ -2527,3 +2544,976 @@ 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_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") + 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]) + + 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) + 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): + """ + 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") + 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")) + 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") + # 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 + + 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] + + @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=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.""" + 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")) + # 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")) + + 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") + 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]) + + 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 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_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") + 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.""" + spool = dc.spool(patch).attach_inventory(inventory) + out = spool.split_by("noisy") + assert sorted(out.get_contents()["noisy"].tolist()) == [False, False, True] + + def test_one_split_partitions_the_channels(self, patch, inventory): + """ + A channel holds one value of a group, so one split cannot share it. + + 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) + 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.""" + 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_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) + 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.""" + 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) + # 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() + ) + + 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 + + +@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]) + + +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.""" + + 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 + + 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 + + @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) diff --git a/tests/test_utils/test_chunk.py b/tests/test_utils/test_chunk.py index 32b384ec..2a52d712 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,17 @@ 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 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) 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 +610,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 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") + 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 +659,83 @@ 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): + """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) + 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") + + +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