diff --git a/docs/api/coordinates.md b/docs/api/coordinates.md index 46df4ff7..ff074d05 100644 --- a/docs/api/coordinates.md +++ b/docs/api/coordinates.md @@ -19,7 +19,6 @@ Methods :toctree: ../_autosummary Coordinates.isdim - Coordinates.get_query Coordinates.to_index Coordinates.equals Coordinates.copy @@ -43,15 +42,10 @@ Attributes :toctree: ../_autosummary Coordinate.dtype - Coordinate.ndim Coordinate.shape Coordinate.size - Coordinate.empty Coordinate.dim - Coordinate.indices Coordinate.values - Coordinate.start - Coordinate.end Coordinate.name ``` @@ -61,14 +55,49 @@ Methods .. autosummary:: :toctree: ../_autosummary - Coordinate.isscalar Coordinate.isdim + Coordinate.isregular Coordinate.equals - Coordinate.to_index - Coordinate.format_index - Coordinate.slice_indexer Coordinate.copy - Coordinate.to_dataarray +``` + +## AxisCoordinate + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + AxisCoordinate +``` + +Attributes + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + AxisCoordinate.ndim + AxisCoordinate.empty + AxisCoordinate.indices + AxisCoordinate.start + AxisCoordinate.end +``` + +Methods + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + AxisCoordinate.isregular + AxisCoordinate.get_sampling_interval + AxisCoordinate.to_regular + AxisCoordinate.get_split_indices + AxisCoordinate.get_discontinuities + AxisCoordinate.get_availabilities + AxisCoordinate.simplify + AxisCoordinate.to_index + AxisCoordinate.to_dataarray ``` ## ScalarCoordinate @@ -106,7 +135,8 @@ Methods DenseCoordinate.from_block DenseCoordinate.get_sampling_interval - DenseCoordinate.get_div_points + DenseCoordinate.to_regular + DenseCoordinate.simplify ``` ## InterpCoordinate @@ -126,6 +156,8 @@ Attributes InterpCoordinate.tie_indices InterpCoordinate.tie_values + InterpCoordinate.sampling_interval + InterpCoordinate.tolerance ``` Methods @@ -135,10 +167,8 @@ Methods :toctree: ../_autosummary InterpCoordinate.from_block + InterpCoordinate.to_regular InterpCoordinate.get_sampling_interval - InterpCoordinate.get_split_indices - InterpCoordinate.get_discontinuities - InterpCoordinate.get_availabilities InterpCoordinate.simplify ``` @@ -171,8 +201,15 @@ Methods SampledCoordinate.from_block SampledCoordinate.get_sampling_interval - SampledCoordinate.get_split_indices - SampledCoordinate.get_discontinuities - SampledCoordinate.get_availabilities + SampledCoordinate.to_regular SampledCoordinate.simplify ``` + +## Functions + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + get_sampling_interval +``` diff --git a/docs/api/index.md b/docs/api/index.md index d16f304d..2f8a3e99 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -13,5 +13,6 @@ picking processing signal synthetics +testing virtual ``` \ No newline at end of file diff --git a/docs/api/synthetics.md b/docs/api/synthetics.md index d318d4d9..93bfcf61 100644 --- a/docs/api/synthetics.md +++ b/docs/api/synthetics.md @@ -10,5 +10,4 @@ wavelet_wavefronts randn_wavefronts - dummy ``` \ No newline at end of file diff --git a/docs/api/testing.md b/docs/api/testing.md new file mode 100644 index 00000000..d26c7f79 --- /dev/null +++ b/docs/api/testing.md @@ -0,0 +1,12 @@ +```{eval-rst} +.. currentmodule:: xdas.testing +``` + +# xdas.testing + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + dummy +``` diff --git a/docs/api/xdas.md b/docs/api/xdas.md index 1a1ce42a..8ebc62ed 100644 --- a/docs/api/xdas.md +++ b/docs/api/xdas.md @@ -31,6 +31,7 @@ concat concatenate concat_coords + get_sampling_interval split plot_availability ``` diff --git a/docs/plan_regular_coordinates.md b/docs/plan_regular_coordinates.md new file mode 100644 index 00000000..cd746293 --- /dev/null +++ b/docs/plan_regular_coordinates.md @@ -0,0 +1,202 @@ +--- +orphan: true +--- + +# Design: regular coordinates — settling the open questions + +Status: implemented on this branch (2026-07-29). +Branch: `feature/fixed-interp-coords`, targeting a PR to `dev` (0.2.8, untagged). + +This document settles the four design questions left open by the regular- +coordinate work, so the remaining implementation (engine emission, docs pass, +failing tests) has a fixed contract to build against. It supersedes the +never-written `docs/plan_propagate_simplify_kwargs.md` referenced by +`concat`'s docstring. + +## Background: the model as implemented + +An `InterpCoordinate` may carry two optional metadata entries: + +- `sampling_interval` — the nominal sample spacing. Its presence is what makes + the coordinate *regular* (`isregular()`). +- `tolerance` — the allowed jitter around that spacing. The validity invariant, + checked at construction (`_is_valid_sampling_interval`), is per continuous + segment: `|num - si * den| <= 2 * tolerance`, evaluated at the dtype + resolution (integer division for datetime64, so sub-resolution drift is + always absorbed). + +`from_block` produces regular coordinates; `_to_regular` enforces or infers a +spacing (raising when it cannot); `simplify(tolerance, reduce, regularize)` +spends an accuracy budget on tie-point reduction and optional promotion to +regular; `_concat` is strict (keeps the spacing only when both sides agree +exactly, takes `max` of tolerances, otherwise drops to irregular). + +## D1. Public API surface: `to_regular` public, `infer_regular` private + +**Decision.** Promote `_to_regular` to public `to_regular`, defined on +`AxisCoordinate` (not just `InterpCoordinate`), honouring the rule that a +public coordinate method exists on the whole axis hierarchy or not at all: + +- `InterpCoordinate.to_regular(sampling_interval=None, tolerance=None)` — + current `_to_regular` behaviour: enforce the given spacing, inferring it when + omitted, raising `ValueError` when the tie points cannot be described by a + single spacing within `tolerance`. +- `SampledCoordinate.to_regular(...)` — regular by construction: with no + arguments return a copy; with explicit arguments validate them against the + stored interval and raise on mismatch. +- `DenseCoordinate.to_regular(...)` — *conversion*: return a regular + `InterpCoordinate` built from the dense values (reduce within `tolerance`, + then enforce the spacing), raising when the values are genuinely irregular. + Returning a different subclass is acceptable: the `to_` prefix already + signals a conversion, and this is the natural "make this axis usable by + signal processing" entry point. + +`_infer_regular` stays private. It is an implementation detail of +`to_regular`/`simplify` (the Chebyshev-center fit); exposing it publicly on +only one subclass would recreate the partial-interface problem, and its +diagnostic value is available through `to_regular`'s behaviour and error +message. `docs/api/coordinates.md` must drop the `infer_regular` entry and the +release notes keep advertising `to_regular` (now truthfully). + +Consequence: `get_sampling_interval` (module level, `core.py:1244`) loses its +`hasattr(coord, "_to_regular")` duck-typing — see D3. + +## D2. What "regular" means per subclass (the Dense question) + +**Decision.** *Regular* means "carries an explicit nominal sampling interval", +uniformly: + +- `InterpCoordinate`: regular iff `sampling_interval` metadata is present. +- `SampledCoordinate`: always regular (the interval is part of its data). +- `DenseCoordinate`: **never regular**. `get_sampling_interval` returns `None` + unconditionally, dropping the current end-to-end average. The average makes + `isregular()` vacuously true for any dense axis and silently hands a + meaningless rate to signal routines on jittery data — the exact failure mode + this branch exists to eliminate. A dense axis that really is evenly sampled + becomes regular explicitly, via `to_regular` (D1) or + `simplify(regularize=True)`. +- `ScalarCoordinate`: `isregular()` moves to the `Coordinate` base and returns + `False` there; `AxisCoordinate` overrides it with the current + `get_sampling_interval() is not None`. This makes the release-notes claim + ("on the base ABC") true and removes the `AttributeError` on scalar coords. + +## D3. The `get_sampling_interval` contract: strict, one choke point + +Three layers, each with a single behaviour: + +1. **Primitive** — `coord.get_sampling_interval(cast=True)`: return the + nominal interval, or `None` when the coordinate is not regular. Never + raises, never infers, O(1). +2. **Conversion** — `coord.to_regular(...)`: the only place inference and + enforcement happen. Raises with an actionable message on genuinely + irregular axes. +3. **Convenience** — `xdas.get_sampling_interval(da, dim)`: return the nominal + interval when the coordinate is regular, otherwise **raise** `ValueError` + telling the user how to fix it (open the files with a `tolerance`, or + `da[dim] = da[dim].to_regular(tolerance=...)`). The current silent + `_to_regular()` fallback is removed: it hides an O(n log n) inference in + every FFT/filter call and only ever succeeds on exactly-uniform axes anyway + (the implicit epsilon tolerance rejects any real jitter), so its benefit is + marginal and its implicitness is not. + + *Amendment (2026-07-30):* data saved by earlier versions carries no + `sampling_interval` metadata, so raising immediately would break every + signal-processing call on existing archives. For one deprecation cycle the + helper therefore falls back to inference on irregular coordinates: it infers + the spacing (and, for `InterpCoordinate`, the minimal tolerance that + validates it via the Chebyshev fit), emits a `FutureWarning` stating both + values and the migration path, and returns the inferred spacing. Dense + coordinates go through the strict `to_regular()` (uniform axes work, jittery + ones still raise — the old end-to-end average was a silent wrong answer not + worth preserving). Raising remains only where no spacing can be inferred at + all. The strict behaviour described above becomes the default when the + deprecation completes. + +**Migration.** All signal-consuming code goes through layer 3 — including +`xdas/signal.py`, which currently open-codes the strict check six times +(`d = coords[dim].get_sampling_interval(); if d is None: raise ...`). Revert +those to the module-level helper so the error message and the policy live in +one place, and keep `fft.py`, `spectral.py`, `atoms/`, `picking.py`, +`miniseed.py` on the helper. Net user-visible behaviour: every signal routine +raises the *same* error on irregular axes, and none of them raise on data +opened through the engines once D5 lands. + +Also fix `DataArrayList`-style compatibility checking +(`routines.py:919-922`): `get_sampling_interval` returning `None` for the +incoming chunk must produce a `CompatibilityError`, not a `TypeError` inside +`np.isclose`. + +## D4. Tolerance semantics and propagation + +**Meaning.** `tolerance` is a *declared jitter bound carried by the +coordinate*: the promise that every continuous segment satisfies +`|num - si * den| <= 2 * tolerance` at the dtype resolution. It is data, not a +processing parameter — processing functions take a *budget* argument that may +default to it. + +**Propagation rules** (R1–R2 already implemented, kept as-is): + +- **R1 — slicing/striding** (`_slice`): spacing scales by the step, tolerance + is preserved. +- **R2 — raw concatenation** (`_concat`): strict; equal spacings are kept with + `max` of tolerances, anything else drops to irregular. Reconciliation is the + job of user-facing routines via `simplify`. +- **R3 — derived rates must carry their quantization error.** Any operation + that synthesizes a new nominal spacing that is not exactly representable in + the coordinate dtype must record the representation error in `tolerance` + instead of claiming `0`. Concretely for `Upsample(factor)` on datetime axes: + `new_delta = delta // factor` truncates, so the coordinate must carry + `tolerance >= (delta - factor * new_delta)` (2 ns in the failing test) on + top of the inherited tolerance. This is what makes chunk seams land within + tolerance of the nominal grid. +- **R4 — `simplify(tolerance=None)` defaults to the coordinate's own stored + tolerance** (falling back to the current zero-like default when the + coordinate has none). Rationale: the coordinate has already declared "my + values are only meaningful to within `tolerance`"; a canonicalisation pass + that refuses to spend that declared slack is pointless strictness. This + applies to `concat(tolerance=None)` too, per-coordinate. `tolerance=False` + keeps its "no simplification" meaning; an explicit scalar overrides. +- **R5 — no unconditional widening.** `InterpCoordinate.simplify` on a regular + coordinate currently stores `self.tolerance + tolerance` whenever `reduce` + runs. Replace with: after reduction, keep the original tolerance if it still + validates, and only widen (to the smallest valid value, bounded by + `self.tolerance + budget`) when it does not. Without this, chunked and + unchunked pipelines can never produce `equals()` coordinates because the + chunked path concatenates and re-simplifies. + +**Why this fixes `test_upsample`.** Each upsampled chunk carries +`sampling_interval = 6_666_666 ns, tolerance = 2 ns` (R3). `_concat` keeps the +spacing (R2). `concat`'s simplify defaults its budget to the stored 2 ns (R4), +Douglas-Peucker drops the seam tie points (they deviate ≤ 2 ns from the global +line), and R5 keeps `tolerance = 2 ns` — identical to the unchunked result. + +**Defaults alignment.** `concat` and `concat_coords` currently disagree +(`regularize=False, tolerance=None` vs `regularize=True, tolerance=False`). +Align `concat_coords` to `concat`: `reduce=True, regularize=False, +tolerance=None` (with R4's meaning). `regularize` stays opt-in for this PR — +with engines emitting regular coordinates (D5) and R2 preserving them, +multi-file opens stay regular without promotion, so the conservative default +costs nothing; flipping it can be revisited once propagation has soaked. + +## D5. IO emission (scope confirmed, design only sketched here) + +Engines construct per-file time/space coordinates with +`InterpCoordinate.from_block(start, size, step)` (the existing `# TODO: use +from_block` sites in `prodml`, `terra15`, `asn`, plus `miniseed.read_stream` +and ObsPy `from_stream`, which must also build at ns resolution to round-trip +`to_stream`). Per-file tolerance is `0`: within one file the grid is exact by +construction. Cross-file jitter is reconciled where it appears — at +`concat`/`open_mfdataarray` time via the user-supplied `tolerance` (R4/R2). +`from_stream` uses `stats.delta`; engines use the file's metadata rate. + +## Acceptance criteria + +- `tests/test_atoms.py::TestFilters::test_upsample` and + `tests/test_dataarray.py::TestIO::test_stream` pass without weakening the + assertions. +- `xd.signal.*`, `xd.fft.*`, `xd.spectral.*`, and the atoms raise one uniform, + actionable error on irregular axes, and raise nothing on engine-opened data. +- Release notes, `docs/api/coordinates.md`, and the user guide describe only + APIs that exist (`to_regular` public, `infer_regular` gone from docs). +- `concat`'s docstring no longer references this document's missing + predecessor. diff --git a/docs/release-notes.md b/docs/release-notes.md index 1f31b2c5..3c6fdfe0 100644 --- a/docs/release-notes.md +++ b/docs/release-notes.md @@ -2,11 +2,19 @@ ## 0.2.8 -### Breaking Changes -- Removed `DefaultCoordinate`, the `Coordinate.is*` predicate properties (`isdense`, `isdefault`, `isinterp`, `issampled`), and `to_dict`/`from_dict` from `DataArray` and coordinate types. These were internal APIs not intended for public use, so this should not affect user code (@atrabattoni). +### New Features +- **Regular coordinates.** A coordinate can now declare a nominal `sampling_interval` (with a `tolerance` bounding the allowed jitter). Query it with `isregular()` / `get_sampling_interval()`; promote an irregular coordinate with `to_regular()`. File engines, `from_block`, and the `fft`/`stft` outputs produce regular coordinates out of the box (@atrabattoni). +- Chunked and unchunked processing now yield identical coordinates: operations that derive a new rate record their rounding error in `tolerance`, and `simplify`/`concat` spend the declared tolerance by default, fusing chunk seams away (@atrabattoni). +- `simplify` gained `reduce` and `regularize` keywords, and the gaps/overlaps API now works on every axis coordinate, including dense ones (@atrabattoni). + +### Deprecations +- The sampling interval is now declared metadata rather than a computed end-to-end average (which was silently wrong on jittery or gappy axes). Data saved by earlier versions carries no declared rate: querying it — e.g. through any signal-processing routine — still works for now, but the rate is inferred and a `FutureWarning` explains how to make the coordinate regular (`da[dim] = da[dim].to_regular(tolerance=...)`). A future release will raise instead (@atrabattoni). ### Refactoring -- `Coordinate` is now a proper ABC with an explicit abstract interface; shared ordered-coordinate logic is consolidated in `SampledMixin`; NumPy 2.0 `copy` keyword compliance (@atrabattoni). +- Reworked the coordinate class hierarchy: `Coordinate` is now a proper ABC and the new `AxisCoordinate` ABC holds the axis-mapping contract shared by dense, interpolated, and sampled coordinates. Use `isinstance(coord, AxisCoordinate)` instead of the removed `is*` predicates (@atrabattoni). +- Cleaned up internal-leaning APIs: removed `DefaultCoordinate`, `to_dict`/`from_dict`, `get_div_points`, `decimate`, and `from_array`; made underscore-private `concat`, `get_indexer`, `get_value`, `format_index`, `slice_index(er)`, `isvalid`, and `get_query`; NumPy 2.0 `copy` keyword compliance (@atrabattoni). +- `concat_coords` now simplifies its result by default, like `concat`; values are unchanged, only redundant tie points are dropped (@atrabattoni). +- Added `xdas.testing.dummy`, a configurable fixture generator replacing `xdas.synthetics.dummy` (@atrabattoni). ## 0.2.7 diff --git a/docs/user-guide/coordinates/index.md b/docs/user-guide/coordinates/index.md index a043b612..7642eca8 100644 --- a/docs/user-guide/coordinates/index.md +++ b/docs/user-guide/coordinates/index.md @@ -18,9 +18,16 @@ metadata) and supports both integer-index access and label-based selection. |:---|:---|:---:|:---| | {py:class}`~xdas.coordinates.ScalarCoordinate` | Scalar metadata, not tied to any axis | `scalar` | scalar-like | | {py:class}`~xdas.coordinates.DenseCoordinate` | One stored value per element | `dense` | `array-like` | -| {py:class}`~xdas.coordinates.InterpCoordinate` | Piecewise-linear from tie points | `interpolated` | `{"tie_indices": array-like[int], "tie_values": array-like}` | +| {py:class}`~xdas.coordinates.InterpCoordinate` | Piecewise-linear from tie points | `interpolated` | `{"tie_indices": array-like[int], "tie_values": array-like}` plus optional `"sampling_interval"` and `"tolerance"` scalars | | {py:class}`~xdas.coordinates.SampledCoordinate` | Uniform grid with optional gaps | `sampled` | `{"tie_values": array-like, "tie_lengths": array-like[int], "sampling_interval": scalar}` | +The three axis-mapping types (`DenseCoordinate`, `InterpCoordinate`, +`SampledCoordinate`) share the {py:class}`~xdas.coordinates.AxisCoordinate` +base, which defines the index/label selection contract. `ScalarCoordinate` +carries a single value with no axis and implements only the thin +{py:class}`~xdas.coordinates.Coordinate` interface. Use +`isinstance(coord, AxisCoordinate)` to test whether a coordinate labels an axis. + ## Creating coordinates {py:class}`~xdas.coordinates.Coordinate` acts as a factory: it inspects the diff --git a/docs/user-guide/coordinates/interpolated-coordinates.md b/docs/user-guide/coordinates/interpolated-coordinates.md index b8aec51b..e5db746b 100644 --- a/docs/user-guide/coordinates/interpolated-coordinates.md +++ b/docs/user-guide/coordinates/interpolated-coordinates.md @@ -57,7 +57,7 @@ coord A major advantage of {py:class}`~xdas.coordinates.InterpCoordinate` is that it enables label-based selection. To retrieve the integer index -corresponding to a given value, use the {py:meth}`~xdas.coordinates.Coordinate.to_index` +corresponding to a given value, use the {py:meth}`~xdas.coordinates.AxisCoordinate.to_index` method: ```{code-cell} @@ -93,6 +93,51 @@ coord = coord.simplify(tolerance=0.0) coord ``` +## Regular coordinates + +An interpolated coordinate can optionally carry a nominal +`sampling_interval` (and a `tolerance` bounding the allowed jitter around +it), making it *regular*. Signal-processing routines (filtering, FFT, +resampling) require a regular coordinate to obtain a clean sample rate; +{py:meth}`~xdas.coordinates.Coordinate.isregular` tells whether a +coordinate carries one. Coordinates built by the file engines or by +{py:meth}`~xdas.coordinates.InterpCoordinate.from_block` are regular out +of the box: + +```{code-cell} +coord = xd.Coordinate( + { + "tie_indices": [0, 9], + "tie_values": [0.0, 90.0], + "sampling_interval": 10.0, + } +) +coord.isregular() +``` + +An irregular coordinate whose values are in fact evenly spaced can be +promoted explicitly with +{py:meth}`~xdas.coordinates.AxisCoordinate.to_regular`, which infers the +spacing when it is not given and raises on genuinely irregular axes. +Data saved by earlier *xdas* versions carries no declared spacing; for +now, signal-processing routines fall back to inferring one and emit a +{py:exc}`FutureWarning` telling you the tolerance required — promote the +coordinate as shown below to silence it: + +```{code-cell} +coord = xd.Coordinate({"tie_indices": [0, 9], "tie_values": [0.0, 90.0]}) +coord.to_regular().get_sampling_interval() +``` + +For jittery axes, pass a `tolerance`: the declared spacing is accepted as +long as every continuous segment stays within it. The stored tolerance +is also the default accuracy budget of +{py:meth}`~xdas.coordinates.InterpCoordinate.simplify`, so chunk seams +introduced by piecewise processing fuse back automatically on +concatenation. `simplify(regularize=True)` combines both steps: it drops +redundant tie points and promotes the result to regular when the +surviving segments admit a single spacing within the budget. + ## Temporal coordinates The most common use of interpolated coordinates in *xdas* is handling diff --git a/docs/user-guide/pipeline/streaming.md b/docs/user-guide/pipeline/streaming.md index f4e9fea4..e2f5d16a 100644 --- a/docs/user-guide/pipeline/streaming.md +++ b/docs/user-guide/pipeline/streaming.md @@ -35,7 +35,7 @@ from xdas.processing import ZMQPublisher, ZMQSubscriber First we generate some data and split it into packets ```{code-cell} -da = xd.synthetics.dummy() +da = xd.testing.dummy() packets = xd.split(da, 5) ``` diff --git a/tests/coordinates/test_coordinates.py b/tests/coordinates/test_coordinates.py index 59e57da3..e54ee118 100644 --- a/tests/coordinates/test_coordinates.py +++ b/tests/coordinates/test_coordinates.py @@ -3,15 +3,22 @@ import xarray as xr import xdas as xd -from xdas.coordinates import DenseCoordinate, InterpCoordinate, ScalarCoordinate -from xdas.coordinates.core import format_datetime, isscalar +from xdas.coordinates import ( + AxisCoordinate, + DenseCoordinate, + InterpCoordinate, + ScalarCoordinate, +) +from xdas.coordinates.core import format_datetime class TestCoordinate: def test_new(self): - assert xd.Coordinate(1).isscalar() + assert isinstance(xd.Coordinate(1), ScalarCoordinate) + assert not isinstance(xd.Coordinate(1), AxisCoordinate) coord = xd.Coordinate(xd.Coordinate([1]), "dim") assert coord.dim == "dim" + assert isinstance(coord, AxisCoordinate) def test_empty(self): with pytest.raises(TypeError, match="cannot infer coordinate type"): @@ -288,8 +295,9 @@ def test_get_sampling_interval_timedelta(self): t0 = np.datetime64("2000-01-01T00:00:00") t1 = np.datetime64("2000-01-01T00:00:10") coord = DenseCoordinate([t0, t1], "time") - result = coord.get_sampling_interval(cast=True) - assert result == 10.0 + assert coord.get_sampling_interval(cast=True) is None + assert not coord.isregular() + assert coord.to_regular().get_sampling_interval(cast=True) == 10.0 def test_format_index_non_integer(self): coord = DenseCoordinate([1, 2, 3], "x") @@ -376,14 +384,36 @@ def test_get_sampling_interval_helper(self): from xdas.coordinates import get_sampling_interval da = xd.DataArray([1, 2, 3], {"x": [10.0, 20.0, 30.0]}) + with pytest.warns(FutureWarning, match="implicit inference is deprecated"): + assert get_sampling_interval(da, "x") == 10.0 + da["x"] = da["x"].to_regular() assert get_sampling_interval(da, "x") == 10.0 - def test_isscalar(self): - assert isscalar(1) - assert isscalar(1.0) - assert isscalar(np.array(1)) - assert not isscalar([1]) - assert not isscalar({"key": "value"}) + def test_get_sampling_interval_helper_jittery_dense_raises(self): + from xdas.coordinates import get_sampling_interval + + da = xd.DataArray([1, 2, 3], {"x": [0.0, 1.0, 5.0]}) + with pytest.raises(ValueError, match="none could be inferred"): + get_sampling_interval(da, "x") + + def test_get_sampling_interval_helper_single_sample_raises(self): + from xdas.coordinates import SampledCoordinate, get_sampling_interval + + coord = SampledCoordinate( + {"tie_values": [0.0], "tie_lengths": [1], "sampling_interval": 5.0}, "x" + ) + da = xd.DataArray([1], {"x": coord}) + with pytest.raises(ValueError, match="none could be inferred"): + get_sampling_interval(da, "x") + + def test_get_sampling_interval_helper_regular(self): + from xdas.coordinates import SampledCoordinate, get_sampling_interval + + coord = SampledCoordinate( + {"tie_values": [0.0], "tie_lengths": [3], "sampling_interval": 5.0} + ) + da = xd.DataArray([1, 2, 3], {"x": coord}) + assert get_sampling_interval(da, "x") == 5.0 def test_format_datetime_no_fractional(self): x = np.datetime64("2000-01-01T00:00:00", "s") @@ -434,3 +464,29 @@ def test_slice_indexer_endpoint_false(self): coord = DenseCoordinate([1.0, 2.0, 3.0], "x") slc = coord._slice_indexer(stop=3.0, endpoint=False) assert slc == slice(None, 2) + + +class TestEncodeDelta: + def test_generic_timedelta_promoted_to_ns(self): + from xdas.coordinates.core import decode_delta, encode_delta + + attrs = encode_delta("tolerance", np.timedelta64(0)) + assert attrs == { + "tolerance": 0, + "tolerance_units": "nanoseconds", + "tolerance_dtype": "timedelta64[ns]", + } + assert decode_delta("tolerance", attrs) == np.timedelta64(0, "ns") + + def test_none_is_omitted(self): + from xdas.coordinates.core import encode_delta + + assert encode_delta("tolerance", None) == {} + + +class TestGetSamplingIntervalHelperNonAxis: + def test_scalar_coordinate_returns_none(self): + from xdas.coordinates import get_sampling_interval + + da = xd.DataArray(np.zeros(3), {"x": [0.0, 1.0, 2.0], "meta": 0}) + assert get_sampling_interval(da, "meta") is None diff --git a/tests/coordinates/test_dense.py b/tests/coordinates/test_dense.py index 95c765df..ef04fa36 100644 --- a/tests/coordinates/test_dense.py +++ b/tests/coordinates/test_dense.py @@ -131,12 +131,53 @@ def test_concat(self): DenseCoordinate(np.array([4.0, 5.0, 6.0], dtype=np.float64)) ) - def test_get_div_points(self): + def test_get_split_indices(self): coord = DenseCoordinate([1, 2, 3, 10, 11, 12]) - div_points = coord.get_div_points(tolerance=3.0) - assert np.array_equal(div_points, [0, 3, 6]) - with pytest.raises(NotImplementedError): - coord.get_div_points() + # local spacing is 1; only the jump 3->10 stands out as a gap, and the + # normal step 10->11 right after it must not be reported as an overlap + np.testing.assert_array_equal( + coord.get_split_indices("discontinuities", tolerance=3.0), [3] + ) + np.testing.assert_array_equal( + coord.get_split_indices("gaps", tolerance=None), [3] + ) + np.testing.assert_array_equal( + coord.get_split_indices("overlaps", tolerance=None), [] + ) + # with no tolerance filtering every consecutive pair is a candidate boundary + np.testing.assert_array_equal(coord.get_split_indices(), [1, 2, 3, 4, 5]) + + def test_get_split_indices_rate_change(self): + # A continuous axis whose sampling rate changes (step 1 then step 2) is + # not a discontinuity: the baseline follows the new rate, so only the + # single transition is reported and the sustained run stays clean. + coord = DenseCoordinate([0, 1, 2, 3, 5, 7, 9]) + np.testing.assert_array_equal( + coord.get_split_indices("gaps", tolerance=0.5), [4] + ) + np.testing.assert_array_equal( + coord.get_split_indices("discontinuities", tolerance=1.5), [] + ) + + def test_get_split_indices_leading_gap(self): + # A discontinuity in the very first step is still detected. + coord = DenseCoordinate([0, 10, 11, 12]) + np.testing.assert_array_equal( + coord.get_split_indices("gaps", tolerance=3.0), [1] + ) + + def test_get_split_indices_empty(self): + coord = DenseCoordinate([]) + np.testing.assert_array_equal(coord.get_split_indices(), []) + np.testing.assert_array_equal( + coord.get_split_indices("gaps", tolerance=None), [] + ) + + def test_simplify_is_noop(self): + coord = DenseCoordinate([1, 2, 3, 10, 11, 12], "x") + result = coord.simplify(tolerance=5.0) + assert result.equals(coord) + assert result is not coord def test_from_block(self): coord = DenseCoordinate.from_block(0, 5, 1, dim="x") @@ -189,3 +230,40 @@ def test_collect_from_dataset_object_dtype(self): dataset["x"] = dataset["x"].astype(object) result = DenseCoordinate._collect_from_dataset(dataset, "x") assert "x" in result + + +class TestDenseCoordinateToRegular: + def test_never_regular(self): + coord = DenseCoordinate([0.0, 1.0, 2.0], "x") + assert coord.get_sampling_interval() is None + assert not coord.isregular() + + def test_to_regular_uniform(self): + coord = DenseCoordinate([0.0, 1.0, 2.0, 3.0], "x") + reg = coord.to_regular() + assert reg.isregular() + assert reg.get_sampling_interval() == 1.0 + assert reg.dim == "x" + np.testing.assert_array_equal(reg.values, coord.values) + + def test_to_regular_explicit_args(self): + coord = DenseCoordinate([0.0, 1.05, 2.0], "x") + reg = coord.to_regular(sampling_interval=1.0, tolerance=0.1) + assert reg.get_sampling_interval() == 1.0 + + def test_to_regular_irregular_raises(self): + coord = DenseCoordinate([0.0, 1.0, 5.0], "x") + with pytest.raises(ValueError, match="not evenly spaced"): + coord.to_regular() + + def test_to_regular_too_short_raises(self): + with pytest.raises(ValueError, match="fewer than two"): + DenseCoordinate([1.0], "x").to_regular() + + def test_to_regular_datetime(self): + t0 = np.datetime64("2000-01-01T00:00:00") + values = t0 + np.timedelta64(1, "s") * np.arange(5) + coord = DenseCoordinate(values, "time") + reg = coord.to_regular() + assert reg.get_sampling_interval() == 1.0 + np.testing.assert_array_equal(reg.values, coord.values) diff --git a/tests/coordinates/test_generic.py b/tests/coordinates/test_generic.py index c75a47b0..75d2e8a4 100644 --- a/tests/coordinates/test_generic.py +++ b/tests/coordinates/test_generic.py @@ -2,6 +2,26 @@ import pytest import xdas as xd +from xdas.coordinates.core import parse_scalar_delta + + +class TestParseScalarDelta: + def test_non_scalar_raises(self): + with pytest.raises(ValueError, match="must be a scalar"): + parse_scalar_delta([1, 2], np.dtype("float64")) + + def test_none_without_default_raises(self): + with pytest.raises(ValueError, match="cannot be None"): + parse_scalar_delta(None, np.dtype("float64")) + + def test_none_with_default_zero(self): + assert parse_scalar_delta(None, np.dtype("float64"), default_zero=True) == 1e-8 + assert parse_scalar_delta(None, np.dtype("float32"), default_zero=True) == 1e-5 + assert parse_scalar_delta(None, np.dtype("float16"), default_zero=True) == 1e-2 + assert parse_scalar_delta(None, np.dtype("int64"), default_zero=True) == 0 + assert parse_scalar_delta( + None, np.dtype("datetime64[s]"), default_zero=True + ) == np.timedelta64(0) class TestFromBlock: @@ -36,7 +56,8 @@ def coord(dtype, ctype): size = 10 step = np.array(1, "timedelta64" if np.issubdtype(dtype, np.datetime64) else dtype) return xd.concat_coords( - [xd.Coordinate[ctype].from_block(start, size, step, "dim") for start in starts] + [xd.Coordinate[ctype].from_block(start, size, step, "dim") for start in starts], + tolerance=False, ) diff --git a/tests/coordinates/test_interp.py b/tests/coordinates/test_interp.py index 8681550a..49a10c04 100644 --- a/tests/coordinates/test_interp.py +++ b/tests/coordinates/test_interp.py @@ -3,7 +3,11 @@ import xarray as xr import xdas as xd -from xdas.coordinates import InterpCoordinate, ScalarCoordinate +from xdas.coordinates import ( + InterpCoordinate, + ScalarCoordinate, +) +from xdas.coordinates.core import Coordinate class TestInterpCoordinate: @@ -53,6 +57,14 @@ def test_isvalid(self): assert InterpCoordinate._isvalid(data) for data in self.invalid: assert not InterpCoordinate._isvalid(data) + # with optional sampling_interval / tolerance is still valid + assert InterpCoordinate._isvalid( + {"tie_indices": [0, 8], "tie_values": [0.0, 8.0], "sampling_interval": 1.0} + ) + # unknown extra key is rejected + assert not InterpCoordinate._isvalid( + {"tie_indices": [0, 8], "tie_values": [0.0, 8.0], "extra": 1} + ) def test_init(self): coord = InterpCoordinate({"tie_indices": [0, 8], "tie_values": [100.0, 900.0]}) @@ -319,8 +331,97 @@ def test_concat(self): assert coord0._concat(coord1).equals(coord1) assert coord1._concat(coord0).equals(coord1) + def test_simplify_preserves_real_discontinuity(self): + # A large jump across a den == 1 gap is preserved as an emergent property + # of the tolerance bound: both boundary points survive while the colinear + # interior collapses. + coord = InterpCoordinate( + { + "tie_indices": [0, 5, 10, 11, 16, 21], + "tie_values": [0.0, 5.0, 10.0, 1000.0, 1005.0, 1010.0], + } + ) + result = coord.simplify() + assert result.equals( + InterpCoordinate( + { + "tie_indices": [0, 10, 11, 21], + "tie_values": [0.0, 10.0, 1000.0, 1010.0], + } + ) + ) + + def test_simplify_absorbs_soft_discontinuity(self): + # A den == 1 gap whose jump fits within tolerance is fused away and the + # two areas merge into a single ramp. + coord = InterpCoordinate( + {"tie_indices": [0, 10, 11, 21], "tie_values": [0.0, 10.0, 11.4, 21.4]} + ) + result = coord.simplify(1.0) + assert result.equals( + InterpCoordinate({"tie_indices": [0, 21], "tie_values": [0.0, 21.4]}) + ) + + def test_simplify_multiple_runs_and_isolated_point(self): + # Two real discontinuities flanking an isolated tie point: each run is + # thinned independently and the isolated point survives. + coord = InterpCoordinate( + { + "tie_indices": [0, 5, 10, 11, 12, 17, 22], + "tie_values": [0.0, 5.0, 10.0, 100.0, 200.0, 205.0, 210.0], + } + ) + result = coord.simplify() + assert result.equals( + InterpCoordinate( + { + "tie_indices": [0, 10, 11, 12, 22], + "tie_values": [0.0, 10.0, 100.0, 200.0, 210.0], + } + ) + ) + + def test_simplify_keeps_kink(self): + # A genuine kink inside a continuous area forces Douglas-Peucker to keep + # the deviating interior point. + coord = InterpCoordinate( + {"tie_indices": [0, 5, 10], "tie_values": [0.0, 100.0, 0.0]} + ) + result = coord.simplify() + assert result.equals(coord) + assert len(coord.simplify(200.0).tie_indices) == 2 + + def test_simplify_datetime_discontinuity(self): + t0 = np.datetime64("2000-01-01T00:00:00") + coord = InterpCoordinate( + { + "tie_indices": [0, 5, 10, 11, 16, 21], + "tie_values": [ + t0, + t0 + np.timedelta64(5, "s"), + t0 + np.timedelta64(10, "s"), + t0 + np.timedelta64(1000, "s"), + t0 + np.timedelta64(1005, "s"), + t0 + np.timedelta64(1010, "s"), + ], + } + ) + result = coord.simplify() + assert np.array_equal(result.tie_indices, [0, 10, 11, 21]) + + +class TestInterpCoordinateExtra: + def test_init_extra_keys(self): + with pytest.raises(TypeError, match="tie_indices"): + InterpCoordinate( + {"tie_indices": [0, 8], "tie_values": [100.0, 900.0], "extra": 1} + ) + + def test_concat_errors(self): with pytest.raises(TypeError): - coord1._concat(ScalarCoordinate(1)) + InterpCoordinate({"tie_indices": [0, 2], "tie_values": [0, 20]})._concat( + ScalarCoordinate(1) + ) with pytest.raises(ValueError, match="different dimension"): InterpCoordinate( {"tie_indices": [0, 2], "tie_values": [0, 20]}, "x" @@ -336,12 +437,6 @@ def test_concat(self): ) ) - def test_init_extra_keys(self): - with pytest.raises(ValueError, match="both"): - InterpCoordinate( - {"tie_indices": [0, 8], "tie_values": [100.0, 900.0], "extra": 1} - ) - def test_init_non_monotonic(self): with pytest.raises(ValueError, match="strictly increasing"): InterpCoordinate( @@ -361,9 +456,10 @@ def test_array_with_dtype(self): result = coord.__array__(dtype=np.float32) assert result.dtype == np.float32 - def test_get_sampling_interval_empty(self): + def test_to_regular_empty(self): coord = InterpCoordinate() - assert coord.get_sampling_interval() is None + with pytest.raises(ValueError, match="cannot infer"): + coord.to_regular() def test_get_indexer_overlaps(self): coord = InterpCoordinate( @@ -399,11 +495,10 @@ def test_get_split_indices_kinds(self): assert len(overlaps) >= 0 def test_get_split_indices_overlaps_tolerance_false(self): - # Build a coord with an overlap (tie_values go backwards between segments) coord = InterpCoordinate( { "tie_indices": [0, 4, 5, 9], - "tie_values": [0.0, 4.0, 3.0, 7.0], # overlap at index 5 (value 3 < 4) + "tie_values": [0.0, 4.0, 3.0, 7.0], } ) result = coord.get_split_indices(kind="overlaps", tolerance=False) @@ -434,7 +529,6 @@ def test_is_monotonic_increasing_false(self): assert coord._is_monotonic_increasing() is False def test_is_monotonic_increasing_multi_segment(self): - # Three segments all strictly increasing — must not raise ValueError from bool() coord = InterpCoordinate( { "tie_indices": [0, 4, 5, 9, 10, 14], @@ -444,8 +538,6 @@ def test_is_monotonic_increasing_multi_segment(self): assert coord._is_monotonic_increasing() is True def test_slice_step_collision(self): - # 4 tie points; step=3 makes first inner tie collide (collision fixed) and - # second inner tie doesn't collide (covers the False branch → loop continues). coord = InterpCoordinate( {"tie_indices": [0, 2, 6, 12], "tie_values": [0.0, 20.0, 60.0, 120.0]} ) @@ -458,19 +550,164 @@ def test_slice_step_collision(self): for i in range(len(result.tie_indices) - 1) ) - def test_get_sampling_interval_datetime_cast(self): + def test_to_regular_explicit_args(self): + coord = InterpCoordinate( + {"tie_indices": [0, 10, 20], "tie_values": [0.0, 1.0, 2.05]} + ) + # strict default tolerance rejects the jitter + with pytest.raises(ValueError, match="not consistent"): + coord.to_regular() + # an explicit tolerance accepts it + reg = coord.to_regular(sampling_interval=0.1, tolerance=0.1) + assert isinstance(reg, InterpCoordinate) + assert reg.isregular() + assert reg.sampling_interval == 0.1 + + def test_to_regular_already_regular_is_preserved(self): + # a regular coordinate keeps its stored spacing untouched + coord = InterpCoordinate( + { + "tie_indices": [0, 10, 20], + "tie_values": [0.0, 1.0, 2.05], + "sampling_interval": 0.1, + "tolerance": 0.1, + } + ) + reg = coord.to_regular() + assert reg is not coord + assert reg.sampling_interval == 0.1 + assert reg.tolerance == 0.1 + # an explicit spacing still overrides it + reg2 = coord.to_regular(sampling_interval=0.103, tolerance=0.1) + assert reg2.sampling_interval == 0.103 + + def test_module_helper_infers_with_warning(self): + da = xd.DataArray( + np.zeros(9), + {"x": {"tie_indices": [0, 8], "tie_values": [0.0, 8.0]}}, + ) + with pytest.warns(FutureWarning, match="implicit inference is deprecated"): + assert xd.get_sampling_interval(da, "x") == 1.0 + da["x"] = da["x"].to_regular() + assert xd.get_sampling_interval(da, "x") == 1.0 + + def test_module_helper_jittery_infers_with_warning(self): + # The fallback states the tolerance required to accept the jitter; the + # strict conversion still rejects it. + da = xd.DataArray( + np.zeros(21), + {"x": {"tie_indices": [0, 10, 20], "tie_values": [0.0, 1.0, 2.05]}}, + ) + with pytest.warns(FutureWarning, match="accepting jitter up to tolerance"): + result = xd.get_sampling_interval(da, "x") + assert 0.1 <= result <= 0.105 + with pytest.raises(ValueError, match="not consistent"): + da["x"].to_regular() + + def test_module_helper_datetime_cast(self): + t0 = np.datetime64("2000-01-01T00:00:00") + da = xd.DataArray( + np.zeros(21), + { + "time": { + "tie_indices": [0, 10, 20], + "tie_values": [ + t0, + t0 + np.timedelta64(10, "s"), + t0 + np.timedelta64(21, "s"), + ], + } + }, + ) + with pytest.warns(FutureWarning, match="implicit inference is deprecated"): + result = xd.get_sampling_interval(da, "time") + assert 1.0 <= result <= 1.1 + with pytest.warns(FutureWarning): + result = xd.get_sampling_interval(da, "time", cast=False) + assert isinstance(result, np.timedelta64) + + def test_module_helper_no_continuous_area_raises(self): + da = xd.DataArray( + np.zeros(2), + {"x": {"tie_indices": [0, 1], "tie_values": [0.0, 1.0]}}, + ) + with pytest.raises(ValueError, match="none could be inferred"): + xd.get_sampling_interval(da, "x") + + def test_to_regular_datetime_cast(self): t0 = np.datetime64("2000-01-01T00:00:00") t1 = np.datetime64("2000-01-01T00:00:08") coord = InterpCoordinate({"tie_indices": [0, 8], "tie_values": [t0, t1]}) - result = coord.get_sampling_interval() # cast=True by default + result = coord.to_regular().get_sampling_interval() # cast=True by default assert result == 1.0 - def test_get_sampling_interval_unit_spaced(self): - # all tie-index gaps == 1 → mask is all False → returns None + def test_to_regular_infer_datetime(self): + t0 = np.datetime64("2000-01-01T00:00:00") + t1 = np.datetime64("2000-01-01T00:00:08") + coord = InterpCoordinate({"tie_indices": [0, 8], "tie_values": [t0, t1]}) + reg = coord.to_regular() + assert reg.sampling_interval == np.timedelta64(1, "s") + assert reg.get_sampling_interval() == 1.0 + + def test_to_regular_unit_spaced(self): + # all tie-index gaps == 1 → no constrained segment → cannot infer coord = InterpCoordinate( {"tie_indices": [0, 1, 2], "tie_values": [0.0, 1.0, 2.0]} ) - assert coord.get_sampling_interval() is None + with pytest.raises(ValueError, match="cannot infer"): + coord.to_regular() + + def test_to_regular_minimax_favours_long_segment(self): + # rates 1.0 (den=10) and 1.1 (den=2); minimax is pulled toward the long + # segment, not the median midpoint 1.05 + coord = InterpCoordinate( + {"tie_indices": [0, 10, 12], "tie_values": [0.0, 10.0, 12.2]} + ) + si = coord.to_regular(tolerance=1.0).sampling_interval + np.testing.assert_allclose(si, 12.2 / 12) + + def test_tolerance_without_sampling_interval(self): + with pytest.raises(ValueError, match="cannot be set without"): + InterpCoordinate( + {"tie_indices": [0, 10], "tie_values": [0.0, 10.0], "tolerance": 0.1} + ) + + def test_infer_regular(self): + # Numeric: rates 1.0 (den=10) and 1.0555 (den=5); the inferred spacing + # and tolerance must round-trip through `to_regular`. + coord = InterpCoordinate( + {"tie_indices": [0, 10, 15], "tie_values": [0.0, 10.0, 15.55]} + ) + si, tol = coord._infer_regular() + assert si > 0 + assert tol > 0 + reg = coord.to_regular(sampling_interval=si, tolerance=tol) + assert reg.isregular() + + # Datetime variant: tolerance comes back as a timedelta64. + t0 = np.datetime64("2000-01-01T00:00:00") + coord_dt = InterpCoordinate( + { + "tie_indices": [0, 10, 15], + "tie_values": [ + t0, + t0 + np.timedelta64(10_000_000_000, "ns"), + t0 + np.timedelta64(15_550_000_000, "ns"), + ], + } + ) + si_dt, tol_dt = coord_dt._infer_regular() + assert np.issubdtype(np.asarray(tol_dt).dtype, np.timedelta64) + assert tol_dt > np.timedelta64(0) + assert coord_dt.to_regular( + sampling_interval=si_dt, tolerance=tol_dt + ).isregular() + + # No continuous segment → nothing to infer. + unit = InterpCoordinate( + {"tie_indices": [0, 1, 2], "tie_values": [0.0, 1.0, 2.0]} + ) + assert unit._infer_regular() == (None, None) def test_add_sub(self): coord = InterpCoordinate({"tie_indices": [0, 4], "tie_values": [10.0, 50.0]}) @@ -501,7 +738,6 @@ def test_to_dataset_collect_roundtrip(self): assert np.allclose(recovered["x"].tie_values, coord.tie_values) def test_to_dataset_multiple_coords_append(self): - # Second coord hitting the "already in attrs" branch (line 223) da = xd.DataArray( np.zeros((9, 5)), { @@ -529,3 +765,385 @@ def test_to_dataset_datetime(self): dataset, attrs = coord._to_dataset(dataset, attrs) assert "time_indices" in dataset assert dataset["time_values"].dtype == np.dtype("datetime64[ns]") + + +class TestInterpCoordinateRegular: + """Tests for InterpCoordinate with an enforced sampling_interval (regular mode).""" + + valid = [ + { + "tie_indices": [0, 5, 9, 10, 19], + "tie_values": [0.0, 0.5, 0.9, 2.0, 2.9], + "sampling_interval": 0.1, + } + ] + + def make(self): + return InterpCoordinate(self.valid[0], "dim") + + def test_isvalid(self): + for data in self.valid: + assert InterpCoordinate._isvalid(data) + # plain tie-point dict without sampling_interval is also valid + assert InterpCoordinate._isvalid( + {"tie_indices": [0, 8], "tie_values": [0.0, 8.0]} + ) + # unknown extra key is rejected + assert not InterpCoordinate._isvalid( + { + "tie_indices": [0, 8], + "tie_values": [0.0, 8.0], + "sampling_interval": 1.0, + "extra": 1, + } + ) + + def test_init(self): + coord = self.make() + assert coord.sampling_interval == 0.1 + assert coord.tolerance is not None + assert coord.dim == "dim" + assert coord.isregular() + + def test_factory_dispatch(self): + coord = Coordinate(self.valid[0]) + assert isinstance(coord, InterpCoordinate) + assert coord.isregular() + # plain tie-point data routes to a non-regular InterpCoordinate + plain = Coordinate({"tie_indices": [0, 8], "tie_values": [0.0, 8.0]}) + assert isinstance(plain, InterpCoordinate) + assert not plain.isregular() + + def test_init_inconsistent(self): + with pytest.raises(ValueError, match="not consistent"): + InterpCoordinate( + { + "tie_indices": [0, 10], + "tie_values": [0.0, 10.0], + "sampling_interval": 0.5, + } + ) + + def test_init_tolerance_allows_jitter(self): + coord = InterpCoordinate( + { + "tie_indices": [0, 10, 20], + "tie_values": [0.0, 1.0, 2.05], + "sampling_interval": 0.1, + "tolerance": 0.1, + } + ) + assert coord.sampling_interval == 0.1 + assert coord.isregular() + + def test_empty(self): + coord = InterpCoordinate() + assert coord.empty + assert coord.sampling_interval is None + assert coord.tolerance is None + assert coord.get_sampling_interval() is None + with pytest.raises(ValueError, match="cannot infer"): + coord.to_regular() + assert not coord.isregular() + + def test_empty_slice_preserves_sampling_interval(self): + coord = InterpCoordinate( + { + "tie_indices": [0, 5, 9, 10, 19], + "tie_values": [0.0, 0.5, 0.9, 2.0, 2.9], + "sampling_interval": 0.1, + } + ) + sliced = coord[0:0] + assert isinstance(sliced, InterpCoordinate) + assert sliced.empty + assert sliced.sampling_interval == 0.1 + + def test_from_block(self): + coord = InterpCoordinate.from_block(0.0, 10, 0.5, "dim") + assert coord.sampling_interval == 0.5 + assert len(coord) == 10 + assert coord.isregular() + + def test_slice(self): + coord = self.make() + sliced = coord[2:12] + assert isinstance(sliced, InterpCoordinate) + assert sliced.isregular() + assert sliced.sampling_interval == 0.1 + stepped = coord[::2] + assert stepped.sampling_interval == 0.2 + + def test_slice_empty(self): + coord = self.make() + empty = coord[5:5] + assert isinstance(empty, InterpCoordinate) + assert empty.empty + + def test_concat(self): + a = InterpCoordinate( + {"tie_indices": [0, 9], "tie_values": [0.0, 0.9], "sampling_interval": 0.1} + ) + b = InterpCoordinate( + {"tie_indices": [0, 9], "tie_values": [1.0, 1.9], "sampling_interval": 0.1} + ) + result = a._concat(b) + assert isinstance(result, InterpCoordinate) + assert result.isregular() + assert result.sampling_interval == 0.1 + assert len(result) == 20 + + def test_concat_different_sampling_interval(self): + # Wildly different rates cannot be reconciled under tolerance, so the + # merged coord falls back to irregular. + a = InterpCoordinate( + {"tie_indices": [0, 9], "tie_values": [0.0, 0.9], "sampling_interval": 0.1} + ) + b = InterpCoordinate( + {"tie_indices": [0, 9], "tie_values": [1.0, 2.8], "sampling_interval": 0.2} + ) + result = a._concat(b) + assert isinstance(result, InterpCoordinate) + assert not result.isregular() + assert result.sampling_interval is None + assert len(result) == 20 + + # Mixed regular/irregular drifts too far → irregular. + c = InterpCoordinate({"tie_indices": [0, 9], "tie_values": [3.0, 4.0]}) + mixed = a._concat(c) + assert mixed.sampling_interval is None + assert len(mixed) == 20 + + def test_concat_coords_recovers_regular_spacing(self): + # `_concat` itself stays strict and drops to irregular when sampling + # intervals disagree; `concat_coords` then tries to reconcile a + # single shared rate within the user-supplied tolerance. + a = InterpCoordinate( + { + "tie_indices": [0, 9], + "tie_values": [0.0, 0.9], + "sampling_interval": 0.1, + "tolerance": 0.05, + }, + "x", + ) + b = InterpCoordinate( + { + "tie_indices": [0, 9], + "tie_values": [1.0, 1.99], + "sampling_interval": 0.11, + "tolerance": 0.05, + }, + "x", + ) + assert not a._concat(b).isregular() + + from xdas.core.routines import concat_coords + + reconciled = concat_coords([a, b], tolerance=0.5, regularize=True) + assert reconciled.isregular() + assert 0.1 <= reconciled.sampling_interval <= 0.11 + assert len(reconciled) == 20 + + def test_add_sub(self): + coord = self.make() + shifted = coord + 1.0 + assert isinstance(shifted, InterpCoordinate) + assert shifted.isregular() + assert shifted.sampling_interval == 0.1 + assert shifted.start == coord.start + 1.0 + back = shifted - 1.0 + assert np.allclose(back.tie_values, coord.tie_values) + + def test_simplify(self): + coord = InterpCoordinate( + { + "tie_indices": [0, 5, 10], + "tie_values": [0.0, 0.5, 1.0], + "sampling_interval": 0.1, + } + ) + simplified = coord.simplify() + assert isinstance(simplified, InterpCoordinate) + assert simplified.isregular() + assert simplified.sampling_interval == 0.1 + + def test_get_sampling_interval_datetime(self): + t0 = np.datetime64("2000-01-01T00:00:00") + coord = InterpCoordinate( + { + "tie_indices": [0, 8], + "tie_values": [t0, t0 + np.timedelta64(8, "s")], + "sampling_interval": np.timedelta64(1, "s"), + } + ) + assert coord.get_sampling_interval() == 1.0 + assert coord.get_sampling_interval(cast=False) == np.timedelta64(1, "s") + assert coord.to_regular().get_sampling_interval() == 1.0 + + def test_dataset_roundtrip_numeric(self): + coord = self.make() + da = xd.DataArray(np.zeros(len(coord)), {"dim": coord}) + dataset = xr.Dataset() + dataset, attrs = da.coords["dim"]._to_dataset(dataset, {}) + dataset["__v__"] = xr.DataArray(np.zeros(len(coord)), dims=["dim"]) + dataset["__v__"].attrs.update(attrs) + recovered = Coordinate._from_dataset(dataset, "__v__") + assert isinstance(recovered["dim"], InterpCoordinate) + assert recovered["dim"].isregular() + assert recovered["dim"].sampling_interval == 0.1 + + def test_dataset_roundtrip_datetime(self): + t0 = np.datetime64("2000-01-01T00:00:00") + coord = InterpCoordinate( + { + "tie_indices": [0, 8], + "tie_values": [t0, t0 + np.timedelta64(8, "s")], + "sampling_interval": np.timedelta64(1, "s"), + } + ) + da = xd.DataArray(np.zeros(9), {"time": coord}) + dataset = xr.Dataset() + dataset, attrs = da.coords["time"]._to_dataset(dataset, {}) + dataset["__v__"] = xr.DataArray(np.zeros(9), dims=["time"]) + dataset["__v__"].attrs.update(attrs) + recovered = Coordinate._from_dataset(dataset, "__v__") + assert isinstance(recovered["time"], InterpCoordinate) + assert recovered["time"].isregular() + assert recovered["time"].sampling_interval == np.timedelta64(1, "s") + + def test_collect_mixed_plain_and_regular(self): + da = xd.DataArray( + np.zeros((20, 9)), + { + "x": self.make(), + "y": {"tie_indices": [0, 8], "tie_values": [0.0, 8.0]}, + }, + ) + dataset = xr.Dataset() + attrs = {} + dataset, attrs = da.coords["x"]._to_dataset(dataset, attrs) + dataset, attrs = da.coords["y"]._to_dataset(dataset, attrs) + dataset["__v__"] = xr.DataArray(np.zeros((20, 9)), dims=["x", "y"]) + dataset["__v__"].attrs.update(attrs) + recovered = Coordinate._from_dataset(dataset, "__v__") + assert isinstance(recovered["x"], InterpCoordinate) + assert recovered["x"].isregular() + assert isinstance(recovered["y"], InterpCoordinate) + assert not recovered["y"].isregular() + + def test_file_roundtrip(self, tmp_path): + coord = self.make() + da = xd.DataArray(np.zeros(len(coord)), {"dim": coord}) + path = tmp_path / "reg.nc" + da.to_netcdf(path) + loaded = xd.open_dataarray(path) + assert isinstance(loaded.coords["dim"], InterpCoordinate) + assert loaded.coords["dim"].isregular() + assert loaded.coords["dim"].sampling_interval == 0.1 + + +class TestDeltaEncoding: + def test_encode_none(self): + from xdas.coordinates.core import encode_delta + + assert encode_delta("sampling_interval", None) == {} + + def test_encode_decode_numeric(self): + from xdas.coordinates.core import decode_delta, encode_delta + + attrs = encode_delta("sampling_interval", 0.1) + assert attrs == {"sampling_interval": 0.1} + assert decode_delta("sampling_interval", attrs) == 0.1 + + def test_decode_missing(self): + from xdas.coordinates.core import decode_delta + + assert decode_delta("sampling_interval", {}) is None + + +class TestSimplifyToleranceDefaults: + def test_default_budget_is_stored_tolerance(self): + # A seam 1.0 off the nominal rate fuses away under the declared + # tolerance of 2.0 without passing any explicit budget. + coord = InterpCoordinate( + { + "tie_indices": [0, 10, 11, 21], + "tie_values": [0.0, 30.0, 34.0, 64.0], + "sampling_interval": 3.0, + "tolerance": 2.0, + } + ) + result = coord.simplify() + assert len(result.tie_indices) == 2 + assert result.sampling_interval == 3.0 + assert result.tolerance == 2.0 + + def test_lossless_pass_keeps_tolerance(self): + coord = InterpCoordinate( + { + "tie_indices": [0, 9], + "tie_values": [0.0, 9.0], + "sampling_interval": 1.0, + "tolerance": 0.5, + } + ) + result = coord.simplify() + assert result.equals(coord) + + def test_widen_only_when_needed(self): + # Fusing a jump beyond the declared tolerance widens it by the budget. + t0 = np.datetime64("2000-01-01T00:00:00", "ns") + s = np.timedelta64(1, "s").astype("m8[ns]") + coord = InterpCoordinate( + { + "tie_indices": [0, 5, 6, 11], + "tie_values": [t0, t0 + 5 * s, t0 + 8 * s, t0 + 13 * s], + "sampling_interval": s, + "tolerance": np.timedelta64(0, "ns"), + } + ) + result = coord.simplify(np.timedelta64(3, "s")) + assert len(result.tie_indices) == 2 + assert result.sampling_interval == s + assert result.tolerance == np.timedelta64(3, "s").astype("m8[ns]") + + +class TestSimplifyNoReduce: + def test_regularize_without_reduce(self): + coord = InterpCoordinate( + {"tie_indices": [0, 5, 10], "tie_values": [0.0, 5.0, 10.0]} + ) + result = coord.simplify(reduce=False, regularize=True) + assert len(result.tie_indices) == 3 + assert result.isregular() + assert result.get_sampling_interval() == 1.0 + + +class TestSimplifyRegularizeFallback: + def test_no_continuous_area_stays_irregular(self): + coord = InterpCoordinate( + {"tie_indices": [0, 1, 2], "tie_values": [0.0, 1.0, 5.0]} + ) + result = coord.simplify(reduce=False, regularize=True) + assert not result.isregular() + + def test_invalid_fit_stays_irregular(self): + coord = InterpCoordinate( + {"tie_indices": [0, 10, 20], "tie_values": [0.0, 1.0, 2.05]} + ) + result = coord.simplify(reduce=False, regularize=True) + assert not result.isregular() + + +class TestFromBlockShort: + def test_single_sample(self): + coord = InterpCoordinate.from_block(0.0, 1, 2.0, dim="x") + assert len(coord) == 1 + assert coord.values == [0.0] + assert coord.sampling_interval == 2.0 + + def test_empty(self): + coord = InterpCoordinate.from_block(0.0, 0, 2.0, dim="x") + assert coord.empty + assert coord.sampling_interval == 2.0 diff --git a/tests/coordinates/test_sampled.py b/tests/coordinates/test_sampled.py index 17e477ff..ea3bbbed 100644 --- a/tests/coordinates/test_sampled.py +++ b/tests/coordinates/test_sampled.py @@ -959,3 +959,32 @@ def test_collect_from_dataset_no_sampling(self): dataset = xr.Dataset({"data": xr.DataArray(np.zeros(3))}) result = SampledCoordinate._collect_from_dataset(dataset, "data") assert result == {} + + +class TestSampledCoordinateToRegular: + def test_returns_copy(self): + coord = SampledCoordinate( + {"tie_values": [0.0], "tie_lengths": [5], "sampling_interval": 2.0}, "x" + ) + reg = coord.to_regular() + assert reg.equals(coord) + assert reg is not coord + + def test_matching_explicit_interval(self): + coord = SampledCoordinate( + {"tie_values": [0.0], "tie_lengths": [5], "sampling_interval": 2.0}, "x" + ) + assert coord.to_regular(sampling_interval=2.0).equals(coord) + + def test_mismatching_interval_raises(self): + coord = SampledCoordinate( + {"tie_values": [0.0], "tie_lengths": [5], "sampling_interval": 2.0}, "x" + ) + with pytest.raises(ValueError, match="does not match"): + coord.to_regular(sampling_interval=3.0) + + def test_mismatch_within_tolerance(self): + coord = SampledCoordinate( + {"tie_values": [0.0], "tie_lengths": [5], "sampling_interval": 2.0}, "x" + ) + assert coord.to_regular(sampling_interval=2.05, tolerance=0.1).equals(coord) diff --git a/tests/coordinates/test_scalar.py b/tests/coordinates/test_scalar.py index b3d88b84..0a9b751f 100644 --- a/tests/coordinates/test_scalar.py +++ b/tests/coordinates/test_scalar.py @@ -3,7 +3,7 @@ import xarray as xr import xdas as xd -from xdas.coordinates import ScalarCoordinate +from xdas.coordinates import AxisCoordinate, ScalarCoordinate class TestScalarCoordinate: @@ -38,16 +38,13 @@ def test_init(self): with pytest.raises(TypeError): ScalarCoordinate(data) - def test_getitem(self): - with pytest.raises(TypeError): - ScalarCoordinate(1)[...] - with pytest.raises(TypeError): - ScalarCoordinate(1)[:] - with pytest.raises(TypeError): - ScalarCoordinate(1)[0] - - def test_len(self): - assert len(ScalarCoordinate(1)) == 1 + def test_not_axis_coordinate(self): + # a scalar coordinate is not an axis coordinate and carries no axis API + coord = ScalarCoordinate(1) + assert not isinstance(coord, AxisCoordinate) + assert not hasattr(coord, "from_block") + assert not hasattr(coord, "_get_value") + assert not hasattr(coord, "to_index") def test_repr(self): for data in self.valid: @@ -92,56 +89,18 @@ def test_equals(self): assert ScalarCoordinate(1).equals(ScalarCoordinate(np.array(1))) assert not ScalarCoordinate(1).equals(42) - def test_to_index(self): - with pytest.raises(NotImplementedError): - ScalarCoordinate(1).to_index("item") - - def test_is_monotonic_increasing(self): - with pytest.raises(TypeError): - ScalarCoordinate(1)._is_monotonic_increasing() - - def test_concat(self): - with pytest.raises(TypeError): - ScalarCoordinate(1)._concat(ScalarCoordinate(2)) - - def test_from_block(self): - with pytest.raises(TypeError): - ScalarCoordinate.from_block(0, 5, 1) - def test_empty(self): with pytest.raises(TypeError, match="cannot be empty"): ScalarCoordinate() - def test_indices(self): - with pytest.raises(TypeError): - ScalarCoordinate(1).indices - - def test_start(self): - with pytest.raises(TypeError): - ScalarCoordinate(1).start - - def test_end(self): - with pytest.raises(TypeError): - ScalarCoordinate(1).end - - def test_get_value(self): - with pytest.raises(TypeError): - ScalarCoordinate(1)._get_value(0) - - def test_get_indexer(self): - with pytest.raises(TypeError): - ScalarCoordinate(1)._get_indexer(1) - - def test_slice(self): - with pytest.raises(TypeError): - ScalarCoordinate(1)._slice(slice(None)) - - def test_get_sampling_interval(self): - assert ScalarCoordinate(1).get_sampling_interval() is None - def test_to_dataset_with_name(self): da = xd.DataArray([1, 2, 3], {"x": [1.0, 2.0, 3.0], "meta": 42}) sc = da.coords["meta"] dataset = xr.Dataset() dataset, attrs = sc._to_dataset(dataset, {}) assert "meta" in dataset.coords + + +class TestScalarCoordinateRegularity: + def test_never_regular(self): + assert not ScalarCoordinate(42).isregular() diff --git a/tests/io/test_asn.py b/tests/io/test_asn.py index 97bb44cd..b8ab6c4d 100644 --- a/tests/io/test_asn.py +++ b/tests/io/test_asn.py @@ -16,26 +16,8 @@ def get_free_local_address(): return f"tcp://localhost:{port}" -coords = { - "time": { - "tie_indices": [0, 99], - "tie_values": [ - np.datetime64("2020-01-01T00:00:00.000000000"), - np.datetime64("2020-01-01T00:00:09.900000000"), - ], - }, - "distance": {"tie_indices": [0, 9], "tie_values": [0.0, 90.0]}, -} - -da_float32 = xd.DataArray( - np.random.randn(100, 10).astype("float32"), - coords, -) - -da_int16 = xd.DataArray( - np.random.randn(100, 10).astype("int16"), - coords, -) +da_float32 = xd.testing.dummy(shape=(100, 10), step=(0.1, 10.0), dtype="float32") +da_int16 = xd.testing.dummy(shape=(100, 10), step=(0.1, 10.0), dtype="int16") class TestASNEngineROIBounds: @@ -229,7 +211,11 @@ def test_one_chunk(self): assert sub.packet_size == 4008 assert sub.shape == (100, 10) assert sub.dtype == np.float32 - assert sub.distance == {"tie_indices": [0, 9], "tie_values": [0.0, 90.0]} + assert sub.distance == { + "tie_indices": [0, 9], + "tie_values": [0.0, 90.0], + "sampling_interval": 10.0, + } assert sub.delta == np.timedelta64(100, "ms") result = next(sub) assert result.equals(da_float32) @@ -249,7 +235,11 @@ def test_several_chunks(self): assert sub.packet_size == 808 assert sub.shape == (20, 10) assert sub.dtype == np.float32 - assert sub.distance == {"tie_indices": [0, 9], "tie_values": [0.0, 90.0]} + assert sub.distance == { + "tie_indices": [0, 9], + "tie_values": [0.0, 90.0], + "sampling_interval": 10.0, + } assert sub.delta == np.timedelta64(100, "ms") for chunk in chunks: result = next(sub) @@ -343,6 +333,7 @@ def test_roiDec(self): assert sub.distance == { "tie_indices": [0, 16001], "tie_values": [0.0, 163418.2435258568], + "sampling_interval": 163418.2435258568 / 16001, } def test_iter(self): diff --git a/tests/io/test_xdas_io.py b/tests/io/test_xdas_io.py index d2e14259..ce0d8888 100644 --- a/tests/io/test_xdas_io.py +++ b/tests/io/test_xdas_io.py @@ -21,19 +21,7 @@ def make_da(): - return xd.DataArray( - np.zeros((10, 5), dtype=np.float32), - { - "time": { - "tie_indices": [0, 9], - "tie_values": [ - np.datetime64("2020-01-01T00:00:00.000000000"), - np.datetime64("2020-01-01T00:00:09.000000000"), - ], - }, - "distance": {"tie_indices": [0, 4], "tie_values": [0.0, 40.0]}, - }, - ) + return xd.testing.dummy(shape=(10, 5), step=(1.0, 10.0), dtype=np.float32) class TestXdasEngineDelegates: diff --git a/tests/test_atoms.py b/tests/test_atoms.py index 83d6854f..42982bb7 100644 --- a/tests/test_atoms.py +++ b/tests/test_atoms.py @@ -18,7 +18,7 @@ UpSample, ) from xdas.signal import lfilter -from xdas.synthetics import randn_wavefronts, wavelet_wavefronts +from xdas.synthetics import randn_wavefronts class TestAbstractAtom: @@ -59,7 +59,7 @@ def test_pickable(self, tmp_path): class TestProcessing: def test_sequence(self): # Generate a temporary dataset - da = wavelet_wavefronts() + da = xd.testing.dummy() # Declare sequence to execute seq = Sequential( @@ -100,10 +100,10 @@ def test_passing_atom(self): class TestFilters: def test_lfilter(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() chunks = xd.split(da, 6, "time") - b, a = sp.iirfilter(4, 10.0, btype="lowpass", fs=50.0) + b, a = sp.iirfilter(4, 10.0, btype="lowpass", fs=100.0) data = sp.lfilter(b, a, da.values, axis=0) expected = da.copy(data=data) @@ -131,10 +131,10 @@ def test_lfilter(self): # assert result.equals(expected) def test_sosfilter(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() chunks = xd.split(da, 6, "time") - sos = sp.iirfilter(4, 10.0, btype="lowpass", fs=50.0, output="sos") + sos = sp.iirfilter(4, 10.0, btype="lowpass", fs=100.0, output="sos") data = sp.sosfilt(sos, da.values, axis=0) expected = da.copy(data=data) @@ -162,7 +162,9 @@ def test_sosfilter(self): # assert result.equals(expected) def test_downsample(self): - da = wavelet_wavefronts() + # size must be a multiple of the decimation factor: on a partial trailing + # phase the chunked path drops one sample that the monolithic one keeps + da = xd.testing.dummy(shape=(102, 10)) chunks = xd.split(da, 6, "time") expected = da.isel(time=slice(None, None, 3)) atom = DownSample(3, "time") @@ -174,28 +176,41 @@ def test_downsample(self): def test_upsample(self): da = xd.DataArray( - [1, 1, 1], {"time": {"tie_indices": [0, 2], "tie_values": [0.0, 6.0]}} + [1, 1, 1], + { + "time": { + "tie_indices": [0, 2], + "tie_values": [0.0, 6.0], + "sampling_interval": 3.0, + } + }, ) expected = xd.DataArray( [3, 0, 0, 3, 0, 0, 3, 0, 0], - {"time": {"tie_indices": [0, 8], "tie_values": [0.0, 8.0]}}, + { + "time": { + "tie_indices": [0, 8], + "tie_values": [0.0, 8.0], + "sampling_interval": 1.0, + } + }, ) atom = UpSample(3, dim="time") result = atom(da) assert result.equals(expected) - da = wavelet_wavefronts() + da = xd.testing.dummy() chunks = xd.split(da, 6, "time") expected = atom(da) result = xd.concat([atom(chunk, chunk_dim="time") for chunk in chunks], "time") assert result.equals(expected) def test_firfilter(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() chunks = xd.split(da, 6, "time") - taps = sp.firwin(11, 0.4, pass_zero="lowpass") + taps = sp.firwin(11, 0.2, pass_zero="lowpass") expected = xs.lfilter(taps, 1.0, da, "time") - expected["time"] -= np.timedelta64(20, "ms") * 5 + expected["time"] -= np.timedelta64(10, "ms") * 5 atom = FIRFilter(11, 10.0, "lowpass", dim="time") result = atom(da) assert result.equals(expected) @@ -209,7 +224,7 @@ def test_firfilter(self): class TestResamplePoly: def test_up_down(self): - da = wavelet_wavefronts() + da = xd.testing.dummy(shape=(300, 10), step=(0.02, 25.0)) # 50 Hz, 6 s chunks = xd.split(da, 6, "time") expected = xs.resample_poly(da, 5, 2, "time") @@ -224,9 +239,9 @@ def test_up_down(self): assert result.attrs == result_chunked.attrs assert result.name == result_chunked.name - result = result.sel(time=slice("2023-01-01T00:00:01", "2023-01-01T00:00:05")) + result = result.sel(time=slice("2024-05-21T00:00:01", "2024-05-21T00:00:05")) expected = expected.sel( - time=slice("2023-01-01T00:00:01", "2023-01-01T00:00:05") + time=slice("2024-05-21T00:00:01", "2024-05-21T00:00:05") ) assert np.allclose(result.values, expected.values, atol=1e-15, rtol=1e-12) assert result.coords.equals(expected.coords) @@ -234,7 +249,7 @@ def test_up_down(self): assert result.name == expected.name def test_nothing_to_do(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() fs = 1 / xd.get_sampling_interval(da, "time") atom = ResamplePoly(fs, maxfactor=10, dim="time") result = atom(da) @@ -329,7 +344,7 @@ def test_partial_state_kwarg(self): assert "key" in p._state def test_partial_stateful_call(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() atom = IIRFilter(4, 10.0, "lowpass", dim="time", stype="ba") da_out = atom(da, chunk_dim="time") assert da_out.shape == da.shape @@ -411,7 +426,7 @@ def test_iirfilter_invalid_stype(self): IIRFilter(4, 10.0, "lowpass", dim="time", stype="invalid") def test_iirfilter_initialize_from_state_zpk_stype(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() atom = IIRFilter(4, 10.0, "lowpass", dim="time", stype="ba") atom(da, chunk_dim="time") atom.stype = "zpk" @@ -419,13 +434,13 @@ def test_iirfilter_initialize_from_state_zpk_stype(self): atom.initialize_from_state() def test_downsample_factor_one(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() atom = DownSample(1, dim="time") result = atom(da) assert result.equals(da) def test_upsample_no_scale(self): - da = wavelet_wavefronts().isel(time=slice(0, 10)) + da = xd.testing.dummy().isel(time=slice(0, 10)) atom = UpSample(2, dim="time", scale=False) result = atom(da) assert result.sizes["time"] == 2 * da.sizes["time"] diff --git a/tests/test_core.py b/tests/test_core.py index 4554818a..0e9fc653 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -8,24 +8,6 @@ class TestCore: - def generate(self, datetime): - shape = (300, 100) - if datetime: - t = { - "tie_indices": [0, shape[0] - 1], - "tie_values": [np.datetime64(0, "ms"), np.datetime64(2990, "ms")], - } - else: - t = {"tie_indices": [0, shape[0] - 1], "tie_values": [0, 3.0 - 1 / 100]} - s = {"tie_indices": [0, shape[1] - 1], "tie_values": [0, 990.0]} - return xd.DataArray( - data=np.random.randn(*shape), - coords={ - "time": t, - "distance": s, - }, - ) - def test_open_mfdataarray(self, tmp_path): wavelet_wavefronts().to_netcdf(tmp_path / "sample.nc") for idx, da in enumerate(wavelet_wavefronts(nchunk=3), start=1): @@ -81,6 +63,7 @@ def test_concatenate(self, tmp_path): "time": { "tie_indices": [0, da1.sizes["time"] + da2.sizes["time"] - 1], "tie_values": [da1["time"][0].values, da2["time"][-1].values], + "sampling_interval": da1.coords["time"].sampling_interval, }, "distance": da1["distance"], } @@ -166,7 +149,11 @@ def test_concatenate(self, tmp_path): result = xd.concat(objs, dim="time") time_values = result["time"].values result["time"] = InterpCoordinate( - {"tie_indices": np.arange(len(time_values)), "tie_values": time_values}, + { + "tie_indices": np.arange(len(time_values)), + "tie_values": time_values, + "sampling_interval": da.coords["time"].sampling_interval, + }, "time", ).simplify() assert result.equals(da) @@ -183,7 +170,7 @@ def test_open_datacollection(self): xd.open_datacollection("not_existing_file.nc") def test_asdataarray(self): - da = self.generate(False) + da = xd.testing.dummy(shape=(300, 100), datetime=False) out = xd.asdataarray(da.to_xarray()) assert np.array_equal(out.data, da.data) for dim in da.dims: diff --git a/tests/test_dataarray.py b/tests/test_dataarray.py index c1608a5b..dd62507d 100644 --- a/tests/test_dataarray.py +++ b/tests/test_dataarray.py @@ -9,7 +9,6 @@ import xdas as xd from xdas.coordinates import Coordinates, DenseCoordinate, InterpCoordinate -from xdas.synthetics import wavelet_wavefronts def generate(dense=False): @@ -126,7 +125,7 @@ def test_cannot_set_dims(self): da.dims = ("other_dim",) def test_data_setter(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() data = np.arange(np.prod(da.shape)).reshape(da.shape) da.data = data assert np.array_equal(da.data, data) @@ -170,7 +169,7 @@ def test_sel(self): assert da.sel(dim=slice(100.0, 300.0)).equals(da[0:3]) assert da.sel(dim=slice(100.0, 300.0), endpoint=False).equals(da[0:2]) # drop - da = wavelet_wavefronts() + da = xd.testing.dummy() result = da.sel(distance=0, method="nearest", drop=True) assert "distance" not in result.coords @@ -225,7 +224,7 @@ def test_sel_item_with_overlaps(self): da.sel(time=0.1, method="nearest") def test_isel(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = da.isel(first=0) excepted = da.isel(time=0) assert result.equals(excepted) @@ -233,7 +232,7 @@ def test_isel(self): excepted = da.isel(distance=0) assert result.equals(excepted) # drop - da = wavelet_wavefronts() + da = xd.testing.dummy() result = da.sel(distance=0, drop=True) assert "distance" not in result.coords @@ -314,7 +313,7 @@ def test_expand_dims(self): class TestManipulation: def test_transpose(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = da.transpose("distance", "time") assert result.dims == ("distance", "time") assert np.array_equal(result.values, da.values.T) @@ -328,7 +327,7 @@ def test_transpose(self): da.transpose("space", "frequency") def test_ufunc(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = np.add(da, 1) assert np.array_equal(result.data, da.data + 1) result = np.add(da, np.ones(da.shape[-1])) @@ -339,7 +338,7 @@ def test_ufunc(self): assert np.array_equal(result.data, da.data + da.data[0]) def test_arithmetics(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = da + 1 assert np.array_equal(result.data, da.data + 1) result = da + np.array(1) @@ -435,13 +434,9 @@ def test_from_xarray(self): assert np.array_equal(result["dim"].values, da["dim"].values) def test_stream(self): - da = wavelet_wavefronts() - da["time"] = { - "tie_indices": da["time"].tie_indices, - "tie_values": da["time"].tie_values.astype("datetime64[us]"), - } + da = xd.testing.dummy() st = da.to_stream(dim={"distance": "time"}) - assert st[0].id == "NET.DAS00001.00.BN1" + assert st[0].id == "NET.DAS00001.00.HN1" assert len(st) == da.sizes["distance"] assert st[0].stats.npts == da.sizes["time"] assert np.datetime64(st[0].stats.starttime.datetime) == da["time"][0].values @@ -482,7 +477,7 @@ def test_netcdf_non_dimensional(self, tmp_path): assert result.equals(da) da_path = tmp_path / "da.nc" - da = wavelet_wavefronts().assign_coords(lon=("distance", np.arange(401))) + da = xd.testing.dummy().assign_coords(lon=("distance", np.arange(10))) da.to_netcdf(da_path) tmp = xd.open_dataarray(da_path) vds_path = tmp_path / "vds.nc" @@ -492,7 +487,7 @@ def test_netcdf_non_dimensional(self, tmp_path): def test_io(self, tmp_path): # both coords interpolated - da = wavelet_wavefronts() + da = xd.testing.dummy() path = tmp_path / "interp.nc" da.to_netcdf(path) da_recovered = xd.DataArray.from_netcdf(path) @@ -665,7 +660,7 @@ def test_repr_dask(self): assert "DaskArray" in r def test_repr_virtual(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() da.to_netcdf(tmp_path / "a.nc") da2 = xd.open(tmp_path / "a.nc") r = repr(da2) @@ -774,12 +769,12 @@ def test_drop_coords(self): assert "x" in result.coords def test_isel_drop_non_scalar(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = da.isel(time=slice(0, 3), drop=True) assert "time" in result.coords def test_sel_drop_non_scalar(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() t0 = da["time"].tie_values[0] t1 = da["time"].tie_values[-1] result = da.sel(time=slice(t0, t1), drop=True) diff --git a/tests/test_datacollection.py b/tests/test_datacollection.py index 65d6d5f7..faa0acbe 100644 --- a/tests/test_datacollection.py +++ b/tests/test_datacollection.py @@ -4,7 +4,6 @@ import xdas as xd import xdas.signal as xs from xdas.core.datacollection import get_depth -from xdas.synthetics import wavelet_wavefronts class TestDataCollection: @@ -18,7 +17,7 @@ def nest(self, da): ) def test_init(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = self.nest(da) data = ( "instrument", @@ -31,7 +30,7 @@ def test_init(self): assert result.equals(dc) def test_io(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection( { "das1": da, @@ -63,7 +62,7 @@ def test_io(self, tmp_path): assert result.equals(dc) def test_io_create_dirs(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection( { "das1": da, @@ -79,7 +78,7 @@ def test_io_create_dirs(self, tmp_path): assert result.equals(dc) def test_depth_counter(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() da.name = "da" dc = self.nest(da) path = tmp_path / "tmp.nc" @@ -94,27 +93,27 @@ def test_depth_counter(self, tmp_path): get_depth(file["instrument/das1/acquisition/0/da"]) == 0 def test_isel(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = self.nest(da) - da_isel = da.isel(distance=slice(100, 200)) - dc_isel = dc.isel(distance=slice(100, 200)) + da_isel = da.isel(distance=slice(2, 5)) + dc_isel = dc.isel(distance=slice(2, 5)) assert self.nest(da_isel).equals(dc_isel) - dc_isel = dc.isel(distance=slice(2000, 3000)) + dc_isel = dc.isel(distance=slice(20, 30)) assert dc_isel["das1"].empty assert dc_isel["das2"].empty def test_sel(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = self.nest(da) - da_sel = da.sel(distance=slice(1000, 2000)) - dc_sel = dc.sel(distance=slice(1000, 2000)) + da_sel = da.sel(distance=slice(20, 50)) + dc_sel = dc.sel(distance=slice(20, 50)) assert self.nest(da_sel).equals(dc_sel) - dc_sel = dc.sel(distance=slice(20000, 30000)) + dc_sel = dc.sel(distance=slice(200, 300)) assert dc_sel["das1"].empty assert dc_sel["das2"].empty def test_query(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = self.nest(da) result = dc.query(instrument="das1", acquisition=0) expected = xd.DataCollection( @@ -130,12 +129,12 @@ def test_query(self): assert result.equals(dc) def test_fields(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = self.nest(da) assert dc.fields == ("instrument", "acquisition") def test_map(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = self.nest(da) atom = xs.decimate(..., 2, ftype="fir") result = dc.map(atom) @@ -144,7 +143,7 @@ def test_map(self): def test_flat_map(self): # DataMapping with DataArrays as direct values - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection({"a": da, "b": da}, "flat") atom = xs.decimate(..., 2, ftype="fir") result = dc.map(atom) @@ -152,14 +151,14 @@ def test_flat_map(self): def test_flat_sequence_map(self): # DataSequence with DataArrays as direct values - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") atom = xs.decimate(..., 2, ftype="fir") result = dc.map(atom) assert result[0].equals(atom(da)) def test_datacollection_from_dataarray(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() # When DataArray is passed, rename and return it result = xd.DataCollection(da, "myname") assert isinstance(result, xd.DataArray) @@ -181,7 +180,7 @@ def test_empty_mapping_repr(self): def test_mapping_reduce(self): import pickle - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection({"a": da}, "test") pickled = pickle.dumps(dc) restored = pickle.loads(pickled) @@ -190,86 +189,86 @@ def test_mapping_reduce(self): def test_sequence_reduce(self): import pickle - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "test") pickled = pickle.dumps(dc) restored = pickle.loads(pickled) assert restored.equals(dc) def test_sequence_fields(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") assert "seq" in dc.fields def test_mapping_equals_false_different_type(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dm = xd.DataCollection({"a": da}, "test") assert not dm.equals(xd.DataCollection([da], "test")) def test_mapping_equals_false_different_name(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dm1 = xd.DataCollection({"a": da}, "name1") dm2 = xd.DataCollection({"a": da}, "name2") assert not dm1.equals(dm2) def test_mapping_equals_false_different_keys(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dm1 = xd.DataCollection({"a": da}, "test") dm2 = xd.DataCollection({"b": da}, "test") assert not dm1.equals(dm2) def test_mapping_equals_false_different_values(self): - da = wavelet_wavefronts() - da2 = wavelet_wavefronts() + da = xd.testing.dummy() + da2 = xd.testing.dummy() da2.data[:] = 0 dm1 = xd.DataCollection({"a": da}, "test") dm2 = xd.DataCollection({"a": da2}, "test") assert not dm1.equals(dm2) def test_sequence_equals_false(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() ds1 = xd.DataCollection([da, da], "seq") ds2 = xd.DataCollection([da, da], "other") assert not ds1.equals(ds2) def test_sequence_equals_false_wrong_type(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() ds = xd.DataCollection([da], "seq") dm = xd.DataCollection({"a": da}, "seq") assert not ds.equals(dm) def test_sequence_load(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") loaded = dc.load() assert isinstance(loaded, type(dc)) def test_mapping_load(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection({"a": da, "b": da}, "test") loaded = dc.load() assert isinstance(loaded, type(dc)) def test_sequence_copy(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") copy = dc.copy() assert copy.equals(dc) def test_sequence_isel(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") result = dc.isel(distance=slice(0, 100)) assert len(result) == 2 def test_sequence_sel(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") result = dc.sel(distance=slice(0, 5000)) assert len(result) == 2 def test_sequence_from_netcdf(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") path = tmp_path / "seq.nc" dc.to_netcdf(path) @@ -277,13 +276,13 @@ def test_sequence_from_netcdf(self, tmp_path): assert result.equals(dc) def test_query_invalid_key_in_sequence(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") with pytest.raises(ValueError, match="query must be a string"): dc.query(seq="bad_string_key") def test_query_invalid_key_in_mapping(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection({"a": da}, "test") with pytest.raises(ValueError, match="query must be a string"): dc.query(test=123) @@ -291,7 +290,7 @@ def test_query_invalid_key_in_mapping(self): def test_from_netcdf_non_sequential_int_keys(self, tmp_path): from xdas.core.datacollection import DataMapping - da = wavelet_wavefronts() + da = xd.testing.dummy() # Create a mapping with non-sequential int keys (gaps) dm = DataMapping({0: da, 2: da}, "test") path = tmp_path / "non_seq.nc" @@ -303,7 +302,7 @@ def test_from_netcdf_non_sequential_int_keys(self, tmp_path): def test_sequence_from_netcdf_direct(self, tmp_path): from xdas.core.datacollection import DataSequence - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = DataSequence([da, da], "seq") path = tmp_path / "seq_direct.nc" dc.to_netcdf(path) @@ -311,20 +310,20 @@ def test_sequence_from_netcdf_direct(self, tmp_path): assert result.equals(dc) def test_sequence_query_slice(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") result = dc.query(seq=slice(0, 1)) assert len(result) == 1 def test_mapping_repr_nonempty(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dm = xd.DataCollection({"a": da}, "test") s = repr(dm) assert "test" in s.lower() or "Test" in s def test_mapping_repr_nested(self): # nested DataMapping → triggers the non-DataArray branch in __repr__ - da = wavelet_wavefronts() + da = xd.testing.dummy() dm = self.nest(da) s = repr(dm) assert "das1" in s @@ -332,39 +331,39 @@ def test_mapping_repr_nested(self): def test_mapping_repr_int_keys(self): from xdas.core.datacollection import DataMapping - da = wavelet_wavefronts() + da = xd.testing.dummy() dm = DataMapping({0: da, 1: da}, "seq") s = repr(dm) assert "0" in s def test_sequence_repr(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection([da, da], "seq") s = repr(dc) assert "seq" in s.lower() or "Seq" in s def test_mapping_copy(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.DataCollection({"a": da}, "test") copy = dc.copy() assert copy.equals(dc) def test_sequence_equals_false_different_length(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() ds1 = xd.DataCollection([da, da], "seq") ds2 = xd.DataCollection([da], "seq") assert not ds1.equals(ds2) def test_sequence_equals_false_different_values(self): - da = wavelet_wavefronts() - da2 = wavelet_wavefronts() + da = xd.testing.dummy() + da2 = xd.testing.dummy() da2.data[:] = 0 ds1 = xd.DataCollection([da], "seq") ds2 = xd.DataCollection([da2], "seq") assert not ds1.equals(ds2) def test_nested_sequence_map(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() inner = xd.DataCollection([da, da], "inner") dc = xd.DataCollection([inner, inner], "outer") atom = xs.decimate(..., 2, ftype="fir") @@ -374,13 +373,13 @@ def test_nested_sequence_map(self): def test_parse_tuple_with_name_given(self): from xdas.core.datacollection import DataMapping - da = wavelet_wavefronts() + da = xd.testing.dummy() # When data is a tuple and name is already provided, unpack the tuple ignoring its name dm = DataMapping(("inner_name", {"a": da}), "outer_name") assert dm.name == "outer_name" def test_parse_datacollection_propagates_name(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() dm = xd.DataCollection({"a": da}, "original_name") # just verify parse propagates name from xdas.core.datacollection import parse @@ -391,7 +390,7 @@ def test_parse_datacollection_propagates_name(self): def test_mapping_map_invalid_item(self): from xdas.core.datacollection import DataMapping - da = wavelet_wavefronts() + da = xd.testing.dummy() dm = DataMapping({"good": da}, "test") # bypass validation to inject an invalid item dict.__setitem__(dm, "bad", "not_a_dataarray") @@ -402,7 +401,7 @@ def test_mapping_map_invalid_item(self): def test_sequence_map_invalid_item(self): from xdas.core.datacollection import DataSequence - da = wavelet_wavefronts() + da = xd.testing.dummy() ds = DataSequence([da], "test") # bypass validation to inject an invalid item list.append(ds, "not_a_dataarray") @@ -411,35 +410,35 @@ def test_sequence_map_invalid_item(self): ds.map(atom) def test_mapping_sel_one_element_becomes_empty(self): - da = wavelet_wavefronts() - da_near = da.sel(distance=slice(0, 4999)) - da_far = da.sel(distance=slice(5000, 10000)) + da = xd.testing.dummy() + da_near = da.sel(distance=slice(0, 45)) + da_far = da.sel(distance=slice(50, 90)) dc = xd.DataCollection({"near": da_near, "far": da_far}, "instrument") - result = dc.sel(distance=slice(0, 2000)) + result = dc.sel(distance=slice(0, 20)) assert set(result.keys()) == {"near"} assert not result["near"].empty def test_mapping_sel_all_elements_become_empty(self): - da = wavelet_wavefronts() - da_near = da.sel(distance=slice(0, 4999)) - da_far = da.sel(distance=slice(5000, 10000)) + da = xd.testing.dummy() + da_near = da.sel(distance=slice(0, 45)) + da_far = da.sel(distance=slice(50, 90)) dc = xd.DataCollection({"near": da_near, "far": da_far}, "instrument") - result = dc.sel(distance=slice(-1000, -1)) + result = dc.sel(distance=slice(-100, -1)) assert len(result) == 0 def test_sequence_sel_one_element_becomes_empty(self): - da = wavelet_wavefronts() - da_near = da.sel(distance=slice(0, 4999)) - da_far = da.sel(distance=slice(5000, 10000)) + da = xd.testing.dummy() + da_near = da.sel(distance=slice(0, 45)) + da_far = da.sel(distance=slice(50, 90)) dc = xd.DataCollection([da_near, da_far], "instrument") - result = dc.sel(distance=slice(0, 2000)) + result = dc.sel(distance=slice(0, 20)) assert len(result) == 1 assert not result[0].empty def test_sequence_sel_all_elements_become_empty(self): - da = wavelet_wavefronts() - da_near = da.sel(distance=slice(0, 4999)) - da_far = da.sel(distance=slice(5000, 10000)) + da = xd.testing.dummy() + da_near = da.sel(distance=slice(0, 45)) + da_far = da.sel(distance=slice(50, 90)) dc = xd.DataCollection([da_near, da_far], "instrument") - result = dc.sel(distance=slice(-1000, -1)) + result = dc.sel(distance=slice(-100, -1)) assert len(result) == 0 diff --git a/tests/test_fft.py b/tests/test_fft.py index 8b403120..f2e5a26f 100644 --- a/tests/test_fft.py +++ b/tests/test_fft.py @@ -6,14 +6,14 @@ class TestRFFT: def test_with_non_dimensional(self): - da = xd.synthetics.wavelet_wavefronts() + da = xd.testing.dummy() da["latitude"] = ("distance", np.arange(da.sizes["distance"])) xfft.rfft(da) class TestInverseTransforms: def test_standard(self): - expected = xd.synthetics.wavelet_wavefronts() + expected = xd.testing.dummy() result = xfft.ifft( xfft.fft(expected, dim={"time": "frequency"}), dim={"frequency": "time"}, @@ -30,7 +30,7 @@ def test_standard(self): assert result[name].equals(expected[name]) def test_real(self): - expected = xd.synthetics.wavelet_wavefronts() + expected = xd.testing.dummy() result = xfft.irfft( xfft.rfft(expected, dim={"time": "frequency"}), expected.sizes["time"], @@ -47,7 +47,7 @@ def test_real(self): assert result[name].equals(expected[name]) def test_real_default_n(self): - expected = xd.synthetics.wavelet_wavefronts() + expected = xd.testing.dummy() expected = expected.isel(time=slice(0, expected.sizes["time"] // 2 * 2)) result = xfft.irfft( xfft.rfft(expected, dim={"time": "frequency"}), diff --git a/tests/test_methods.py b/tests/test_methods.py index 0c44bfb5..f527c628 100644 --- a/tests/test_methods.py +++ b/tests/test_methods.py @@ -8,23 +8,13 @@ @pytest.fixture def float_da(): - return xd.DataArray( - data=np.arange(12.0).reshape(3, 4), - coords={ - "x": {"tie_indices": [0, 2], "tie_values": [0.0, 2.0]}, - "y": {"tie_indices": [0, 3], "tie_values": [0.0, 3.0]}, - }, - ) + return xd.testing.dummy(dims=("x", "y"), shape=(3, 4), step=1.0, datetime=False) @pytest.fixture def int_da(): - return xd.DataArray( - data=np.arange(12).reshape(3, 4), - coords={ - "x": {"tie_indices": [0, 2], "tie_values": [0.0, 2.0]}, - "y": {"tie_indices": [0, 3], "tie_values": [0.0, 3.0]}, - }, + return xd.testing.dummy( + dims=("x", "y"), shape=(3, 4), step=1.0, datetime=False, dtype=int ) diff --git a/tests/test_numpy.py b/tests/test_numpy.py index 0090303a..8039e74e 100644 --- a/tests/test_numpy.py +++ b/tests/test_numpy.py @@ -1,13 +1,14 @@ import numpy as np import pytest +import xdas as xd from xdas.core.dataarray import HANDLED_NUMPY_FUNCTIONS, DataArray from xdas.synthetics import wavelet_wavefronts class TestUfuncs: def test_unitary_operators(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = np.abs(da) expected = da.copy(data=np.abs(da.data)) da_out = da.copy() @@ -19,8 +20,8 @@ def test_unitary_operators(self): assert da_where.equals(da) def test_binary_operators(self): - da1 = wavelet_wavefronts() - da2 = wavelet_wavefronts() + da1 = xd.testing.dummy() + da2 = xd.testing.dummy() result = np.add(da1, da2) expected = da1.copy(data=da1.data + da2.data) da_out = da1.copy() @@ -34,7 +35,7 @@ def test_binary_operators(self): np.add(da1, da2[1:]) def test_multiple_outputs(self): - da = wavelet_wavefronts() + da = wavelet_wavefronts() # divmod(da, da) needs non-zero values result1, result2 = np.divmod(da, da) expected1 = da.copy(data=np.ones(da.shape)) expected2 = da.copy(data=np.zeros(da.shape)) @@ -46,7 +47,8 @@ def test_multiple_outputs(self): class TestFunc: def test_returns_dataarray(self): - da = wavelet_wavefronts() + # keep values small: np.i0 overflows on the default dummy + da = xd.testing.dummy(shape=(10, 5)) for numpy_function in HANDLED_NUMPY_FUNCTIONS: if numpy_function == np.clip: result = numpy_function(da, -1, 1) @@ -76,7 +78,7 @@ def test_returns_dataarray(self): assert isinstance(result, DataArray) def test_reduce(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = np.sum(da) assert result.shape == () result = np.sum(da, axis=0) @@ -93,7 +95,7 @@ def test_reduce(self): np.sum(da, axis=2) def test_out(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() out = da.copy() np.cumsum(da, axis=-1, out=out) assert not out.equals(da) diff --git a/tests/test_picking.py b/tests/test_picking.py index 6f844d27..3ee47601 100644 --- a/tests/test_picking.py +++ b/tests/test_picking.py @@ -13,6 +13,7 @@ def generate(self): "distance": { "tie_indices": [0, 4], "tie_values": [0.0, 400.0], + "sampling_interval": 100.0, }, "time": { "tie_indices": [0, 9], @@ -20,6 +21,7 @@ def generate(self): np.datetime64("2023-01-01T00:00:00"), np.datetime64("2023-01-01T00:00:09"), ], + "sampling_interval": np.timedelta64(1, "s"), }, }, ) @@ -248,6 +250,7 @@ def test_scalar_coord_preserved(self): np.datetime64("2023-01-01T00:00:00"), np.datetime64("2023-01-01T00:00:09"), ], + "sampling_interval": np.timedelta64(1, "s"), }, "station": "ABC", }, @@ -280,6 +283,7 @@ def test_non_dim_coord_on_dim_axis_skipped(self): np.datetime64("2023-01-01T00:00:00"), np.datetime64("2023-01-01T00:00:09"), ], + "sampling_interval": np.timedelta64(1, "s"), }, "quality": ( "time", diff --git a/tests/test_processing.py b/tests/test_processing.py index 8d7c60bd..c284e974 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -14,12 +14,11 @@ import xdas.processing as xp from xdas.atoms import Partial, Sequential from xdas.signal import sosfilt -from xdas.synthetics import wavelet_wavefronts class TestDataArrayLoader: def test_init(self): - da = xd.DataArray(np.random.rand(1000, 100), dims=("time", "distance")) + da = xd.testing.dummy(shape=(1000, 100)) dl = xp.DataArrayLoader(da, {"time": 100}) assert dl.da is da assert dl.chunk_dim == "time" @@ -38,14 +37,14 @@ def test_init(self): ], ) def test_chunks_integrity(self, max_buffers, max_workers): - da = xd.DataArray(np.random.rand(1000, 100), dims=("time", "distance")) + da = xd.testing.dummy(shape=(1000, 100)) dl = xp.DataArrayLoader(da, {"time": 100}, max_buffers, max_workers) chunks = [chunk for chunk in dl] result = xd.concat(chunks) assert result.equals(da) def test_error_handling(self): - da = xd.DataArray(np.random.rand(1000, 100), dims=("time", "distance")) + da = xd.testing.dummy(shape=(1000, 100)) with pytest.raises(TypeError): xp.DataArrayLoader(None, None) with pytest.raises(TypeError): @@ -71,7 +70,7 @@ def test_init(self, tmp_path): ], ) def test_chunk_integrity(self, max_buffers, max_workers, tmp_path): - expected = xd.DataArray(np.random.rand(1000, 100), dims=("time", "distance")) + expected = xd.testing.dummy(shape=(1000, 100)) dw = xp.DataArrayWriter(tmp_path, None, max_buffers, max_workers) chunks = xd.split(expected, 10, dim="time") for chunk in chunks: @@ -96,7 +95,7 @@ def test_stateful(self, tmp_path): sample_path = tmp_path / "sample.nc" # generate test dataarray - wavelet_wavefronts().to_netcdf(sample_path) + xd.testing.dummy().to_netcdf(sample_path) da = xd.open(sample_path) # declare processing sequence @@ -117,13 +116,7 @@ def test_stateful(self, tmp_path): assert result1.equals(result2) def test_small_last_chunk(self, tmp_path): - da = xd.DataArray( - data=np.random.randn(1001, 100), - coords={ - "time": xd.Coordinate["interpolated"].from_block(0, 1001, 0.01), - "distance": xd.Coordinate["interpolated"].from_block(0, 100, 10.0), - }, - ) + da = xd.testing.dummy(shape=(1001, 100), datetime=False) # declare processing sequence sos = sp.iirfilter(4, 0.1, btype="lowpass", output="sos") @@ -233,7 +226,7 @@ def publish(): return xd.concat(result) def test_publish_and_subscribe(self): - expected = xd.synthetics.dummy() + expected = xd.testing.dummy() packets = xd.split(expected, 10) address = f"tcp://localhost:{xd.io.get_free_port()}" @@ -241,7 +234,7 @@ def test_publish_and_subscribe(self): assert result.equals(expected) def test_encoding(self): - expected = xd.synthetics.dummy() + expected = xd.testing.dummy() packets = xd.split(expected, 10) address = f"tcp://localhost:{xd.io.get_free_port()}" encoding = {"chunks": (10, 10), **hdf5plugin.Zfp(accuracy=1e-6)} @@ -265,6 +258,7 @@ def test_without_gap(self, tmp_path): "time": { "tie_indices": [0, data.shape[0] - 1], "tie_values": [starttime, endtime], + "sampling_interval": np.timedelta64(10, "ms"), }, "distance": distance, }, @@ -324,6 +318,7 @@ def test_with_gap(self, tmp_path): ], dtype="datetime64[ms]", ), + "sampling_interval": np.timedelta64(10, "ms"), }, "distance": 5.0 * np.arange(10), }, @@ -383,6 +378,7 @@ def test_flat(self, tmp_path): "time": { "tie_indices": [0, data.shape[0] - 1], "tie_values": [starttime, endtime], + "sampling_interval": np.timedelta64(10, "ms"), }, "distance": distance, }, @@ -425,7 +421,7 @@ def atom(da, **kwargs): class TestProcessNoNbytes: def test_loader_without_nbytes(self, tmp_path): - da = xd.DataArray(np.random.rand(100, 10), dims=("time", "distance")) + da = xd.testing.dummy(shape=(100, 10)) chunks = xd.split(da, 10, dim="time") class SimpleLoader: @@ -445,7 +441,7 @@ def atom(x, **kw): class TestDataArrayLoaderMaxBuffers: def test_max_buffers_exceeds_chunks(self): - da = xd.DataArray(np.random.rand(10, 5), dims=("time", "distance")) + da = xd.testing.dummy(shape=(10, 5)) dl = xp.DataArrayLoader(da, {"time": 5}, max_buffers=10) chunks = list(dl) assert len(chunks) == 2 @@ -488,7 +484,7 @@ class TestZMQPublisherAliases: def test_write_alias(self): address = f"tcp://localhost:{xd.io.get_free_port()}" publisher = xp.ZMQPublisher(address) - da = xd.synthetics.dummy() + da = xd.testing.dummy() publisher.write(da) # use write() alias def test_result_returns_none(self): @@ -503,19 +499,7 @@ def test_on_closed(self, tmp_path): from xdas.processing.core import Handler - da = xd.DataArray( - np.zeros((10, 5), dtype=np.float32), - { - "time": { - "tie_indices": [0, 9], - "tie_values": [ - np.datetime64("2020-01-01T00:00:00.000000000"), - np.datetime64("2020-01-01T00:00:09.000000000"), - ], - }, - "distance": {"tie_indices": [0, 4], "tie_values": [0.0, 40.0]}, - }, - ) + da = xd.testing.dummy(shape=(10, 5), step=(1.0, 10.0), dtype=np.float32) path = str(tmp_path / "test.nc") da.to_netcdf(path) @@ -538,19 +522,7 @@ def test_iter_and_next(self, tmp_path): assert iter(loader) is loader # put a DataArray directly into the queue - da = xd.DataArray( - np.zeros((5, 3), dtype=np.float32), - { - "time": { - "tie_indices": [0, 4], - "tie_values": [ - np.datetime64("2020-01-01T00:00:00.000000000"), - np.datetime64("2020-01-01T00:00:04.000000000"), - ], - }, - "distance": {"tie_indices": [0, 2], "tie_values": [0.0, 20.0]}, - }, - ) + da = xd.testing.dummy(shape=(5, 3), step=(1.0, 10.0), dtype=np.float32) loader.queue.put(da) result = next(loader) assert result.equals(da) diff --git a/tests/test_routines.py b/tests/test_routines.py index 6dc38b60..e7835c89 100644 --- a/tests/test_routines.py +++ b/tests/test_routines.py @@ -13,7 +13,15 @@ def test_bag_initialization(self): def test_bag_append_initializes(self): da = xd.DataArray( - np.random.rand(10, 5), {"time": np.arange(10), "space": np.arange(5)} + np.random.rand(10, 5), + { + "time": { + "tie_indices": [0, 9], + "tie_values": [0.0, 9.0], + "sampling_interval": 1.0, + }, + "space": np.arange(5), + }, ) bag = Bag(dim="time") bag.append(da) @@ -90,12 +98,24 @@ def test_bag_append_incompatible_sampling_interval(self): da1 = xd.DataArray( np.random.rand(10, 5), dims=("time", "space"), - coords={"time": np.arange(10)}, + coords={ + "time": { + "tie_indices": [0, 9], + "tie_values": [0.0, 9.0], + "sampling_interval": 1.0, + } + }, ) da2 = xd.DataArray( np.random.rand(10, 5), dims=("time", "space"), - coords={"time": np.arange(10) * 2}, + coords={ + "time": { + "tie_indices": [0, 9], + "tie_values": [0.0, 18.0], + "sampling_interval": 2.0, + } + }, ) bag = Bag(dim="time") bag.append(da1) @@ -112,14 +132,7 @@ def test_basic(self): assert combined.shape == (20, 5) # with coords - da1 = xd.DataArray( - np.random.rand(10, 5), - coords={"time": np.arange(10), "space": np.arange(5)}, - ) - da2 = xd.DataArray( - np.random.rand(10, 5), - coords={"time": np.arange(10, 20), "space": np.arange(5)}, - ) + da1, da2 = xd.split(xd.testing.dummy(dims=("time", "space"), shape=(20, 5)), 2) combined = xd.combine_by_coords([da1, da2], dim="time", squeeze=True) assert combined.shape == (20, 5) @@ -169,12 +182,24 @@ def test_incompatible_sampling_interval(self): da1 = xd.DataArray( np.random.rand(10, 5), dims=("time", "space"), - coords={"time": np.arange(10)}, + coords={ + "time": { + "tie_indices": [0, 9], + "tie_values": [0.0, 9.0], + "sampling_interval": 1.0, + } + }, ) da2 = xd.DataArray( np.random.rand(10, 5), dims=("time", "space"), - coords={"time": np.arange(10) * 2}, + coords={ + "time": { + "tie_indices": [0, 9], + "tie_values": [0.0, 18.0], + "sampling_interval": 2.0, + } + }, ) dc = xd.combine_by_coords([da1, da2], dim="time") assert len(dc) == 2 @@ -200,13 +225,7 @@ def test_expand_scalar_coordinate(self): class TestOpenMFDataArray: def test_warn_on_corrupted_files(self, tmp_path): - expected = xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, # TODO: should work without coords - ) + expected = xd.testing.dummy(dims=("time", "space"), shape=(10, 5)) for index, chunk in enumerate(xd.split(expected, 3, "time"), start=1): chunk.to_netcdf(tmp_path / f"chunk_{index}.nc") result = xd.open_mfdataarray(tmp_path / "*.nc") @@ -225,26 +244,14 @@ def test_warn_on_corrupted_files(self, tmp_path): assert result.equals(expected) def test_verbose_single_worker(self, tmp_path): - expected = xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, # TODO: should work without coords - ) + expected = xd.testing.dummy(dims=("time", "space"), shape=(10, 5)) for index, chunk in enumerate(xd.split(expected, 3, "time"), start=1): chunk.to_netcdf(tmp_path / f"chunk_{index}.nc") result = xd.open_mfdataarray(tmp_path / "*.nc", verbose=True, parallel=1) assert result.equals(expected) def test_verbose_multiple_workers(self, tmp_path): - expected = xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, # TODO: should work without coords - ) + expected = xd.testing.dummy(dims=("time", "space"), shape=(10, 5)) for index, chunk in enumerate(xd.split(expected, 3, "time"), start=1): chunk.to_netcdf(tmp_path / f"chunk_{index}.nc") result = xd.open_mfdataarray(tmp_path / "*.nc", verbose=True, parallel=2) @@ -253,13 +260,7 @@ def test_verbose_multiple_workers(self, tmp_path): class TestOpen: # TODO: those tests are weirdly slow... def test_open_single_dataarray(self, tmp_path): - expected = xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, - ) + expected = xd.testing.dummy(dims=("time", "space"), shape=(10, 5)) path = tmp_path / "dataarray.nc" expected.to_netcdf(path) @@ -268,13 +269,7 @@ def test_open_single_dataarray(self, tmp_path): assert result.equals(expected) def test_open_multiple_file_dataarray(self, tmp_path): - expected = xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, - ) + expected = xd.testing.dummy(dims=("time", "space"), shape=(10, 5)) file_paths = [] for index, chunk in enumerate(xd.split(expected, 3, "time"), start=1): @@ -298,27 +293,11 @@ def test_open_multiple_file_tree(self, tmp_path): expected = xd.DataCollection( { "DAS01": xd.DataCollection( - [ - xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, - ) - ], + [xd.testing.dummy(dims=("time", "space"), shape=(10, 5))], name="acquisition", ), "DAS02": xd.DataCollection( - [ - xd.DataArray( - np.random.rand(7, 3), - coords={ - "time": np.arange(7), - "space": np.arange(3), - }, - ) - ], + [xd.testing.dummy(dims=("time", "space"), shape=(7, 3))], name="acquisition", ), }, @@ -338,15 +317,7 @@ def test_open_multiple_file_tree(self, tmp_path): def test_open_single_datacollection(self, tmp_path): expected = xd.DataCollection( - [ - xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, - ) - ] + [xd.testing.dummy(dims=("time", "space"), shape=(10, 5))] ) expected.to_netcdf(tmp_path / "collection.nc") @@ -358,27 +329,11 @@ def test_open_multiple_datacollection_with_glob(self, tmp_path): expected = xd.DataCollection( { "DAS01": xd.DataCollection( - [ - xd.DataArray( - np.random.rand(10, 5), - coords={ - "time": np.arange(10), - "space": np.arange(5), - }, - ) - ], + [xd.testing.dummy(dims=("time", "space"), shape=(10, 5))], name="acquisition", ), "DAS02": xd.DataCollection( - [ - xd.DataArray( - np.random.rand(7, 3), - coords={ - "time": np.arange(7), - "space": np.arange(3), - }, - ) - ], + [xd.testing.dummy(dims=("time", "space"), shape=(7, 3))], name="acquisition", ), }, @@ -436,7 +391,8 @@ def dataarray(self, dtype, ctype): [ xd.Coordinate[ctype].from_block(start, size, step, "dim") for start in starts - ] + ], + tolerance=False, ) return xd.DataArray(np.random.randn(len(coord)), {"dim": coord}) @@ -514,7 +470,7 @@ def test_invalid_paths_type_raises(self): xd.open(123) def test_callable_engine(self, tmp_path): - da = xd.DataArray(np.random.rand(10, 5), dims=("time", "distance")) + da = xd.testing.dummy(shape=(10, 5)) path = str(tmp_path / "test.nc") da.to_netcdf(path) @@ -525,7 +481,7 @@ def my_engine(fname, **kwargs): assert result.equals(da) def test_invalid_engine_type_raises(self, tmp_path): - da = xd.DataArray(np.random.rand(10, 5), dims=("time", "distance")) + da = xd.testing.dummy(shape=(10, 5)) path = str(tmp_path / "test.nc") da.to_netcdf(path) with pytest.raises(ValueError, match="engine"): @@ -542,7 +498,7 @@ def test_empty_glob_raises(self, tmp_path): xd.open_mfdatacollection(str(tmp_path / "*.nc")) def test_verbose_single_worker(self, tmp_path): - da = xd.DataArray(np.random.rand(10, 5), dims=("time", "distance")) + da = xd.testing.dummy(shape=(10, 5)) dc = xd.DataCollection([da, da]) path1 = str(tmp_path / "dc1.nc") path2 = str(tmp_path / "dc2.nc") @@ -554,7 +510,7 @@ def test_verbose_single_worker(self, tmp_path): assert isinstance(result, xd.DataCollection) def test_verbose_multiple_worker(self, tmp_path): - da = xd.DataArray(np.random.rand(10, 5), dims=("time", "distance")) + da = xd.testing.dummy(shape=(10, 5)) dc = xd.DataCollection([da, da]) path1 = str(tmp_path / "dc1.nc") path2 = str(tmp_path / "dc2.nc") @@ -576,10 +532,7 @@ def test_invalid_paths_type_raises(self): xd.open_mfdataarray(123) def test_parallel_path(self, tmp_path): - expected = xd.DataArray( - np.random.rand(10, 5), - coords={"time": np.arange(10), "space": np.arange(5)}, - ) + expected = xd.testing.dummy(dims=("time", "space"), shape=(10, 5)) for i, chunk in enumerate(xd.split(expected, 3, "time"), 1): chunk.to_netcdf(tmp_path / f"chunk_{i}.nc") result = xd.open_mfdataarray(tmp_path / "*.nc", parallel=2) @@ -592,7 +545,7 @@ def test_no_files_no_failures_raises(self, tmp_path): class TestOpenMFDatacollectionParallel: def test_parallel_path(self, tmp_path): - da = xd.DataArray(np.random.rand(10, 5), dims=("time", "distance")) + da = xd.testing.dummy(shape=(10, 5)) dc = xd.DataCollection([da, da]) path1 = str(tmp_path / "dc1.nc") path2 = str(tmp_path / "dc2.nc") @@ -608,11 +561,9 @@ def test_one_level_depth(self, tmp_path): dirnames = [tmp_path / key for key in keys] for dirname in dirnames: dirname.mkdir() - for idx, da in enumerate( - xd.synthetics.wavelet_wavefronts(nchunk=3), start=1 - ): + for idx, da in enumerate(xd.split(xd.testing.dummy(), 3), start=1): da.to_netcdf(dirname / f"{idx:03d}.nc") - da = xd.synthetics.wavelet_wavefronts() + da = xd.testing.dummy() dc = xd.open_mfdatatree(tmp_path / "{node}" / "00[acquisition].nc") assert list(dc.keys()) == keys for key in keys: @@ -622,11 +573,11 @@ def test_two_level_depth(self, tmp_path): dc = xd.DataCollection( { "NET01": { - "STA01": xd.synthetics.wavelet_wavefronts(nchunk=1), + "STA01": xd.split(xd.testing.dummy(), 1), }, "NET02": { - "STA02": xd.synthetics.wavelet_wavefronts(nchunk=2), - "STA03": xd.synthetics.wavelet_wavefronts(nchunk=3), + "STA02": xd.split(xd.testing.dummy(), 2), + "STA03": xd.split(xd.testing.dummy(), 3), }, } ) @@ -696,7 +647,9 @@ def test_mixed_empty_and_nonempty_uses_nonempty(self): class TestConcatCoordsEdgeCases: - def test_tolerance_with_dense_coord_raises(self): + def test_tolerance_with_dense_coord_is_noop(self): + # Dense coordinates now implement a (degenerate) `simplify`, so passing a + # tolerance no longer raises; it simply has no effect. da1 = xd.DataArray( np.random.rand(5), {"x": np.array([0.0, 1.0, 2.0, 3.0, 4.0])} ) @@ -705,18 +658,32 @@ def test_tolerance_with_dense_coord_raises(self): ) from xdas.core.routines import concat_coords + result = concat_coords([da1["x"], da2["x"]], tolerance=1.0) + expected = concat_coords([da1["x"], da2["x"]]) + assert result.equals(expected) + + def test_tolerance_with_scalar_coord_raises(self): + from xdas.core.routines import concat_coords + + scalar = xd.Coordinate("SRN") with pytest.raises(TypeError, match="tolerance"): - concat_coords([da1["x"], da2["x"]], tolerance=1.0) + concat_coords([scalar], tolerance=1.0) + + def test_default_tolerance_with_scalar_coord_passes(self): + from xdas.core.routines import concat_coords + + scalar = xd.Coordinate("SRN") + assert concat_coords([scalar]).equals(scalar) class TestSplitEdgeCases: def test_n_zero_raises(self): - da = xd.DataArray(np.random.rand(10), dims=("time",)) + da = xd.testing.dummy(dims=("time",), shape=(10,), step=0.01) with pytest.raises(ValueError, match="`n` must be larger than 0"): xd.split(da, 0) def test_n_too_large_raises(self): - da = xd.DataArray(np.random.rand(10), dims=("time",)) + da = xd.testing.dummy(dims=("time",), shape=(10,), step=0.01) with pytest.raises(ValueError, match="`n` must be smaller"): xd.split(da, 10) @@ -737,51 +704,18 @@ def test_scalar_coord_skipped(self): class TestPlotAvailability: def test_dataarray_plot(self): - da = xd.DataArray( - np.random.rand(100), - { - "time": { - "tie_indices": [0, 99], - "tie_values": [ - np.datetime64("2020-01-01"), - np.datetime64("2020-01-01T00:00:09.900"), - ], - } - }, - ) + da = xd.testing.dummy(dims=("time",), shape=(100,), step=0.01) fig = xd.plot_availability(da) assert fig is not None def test_datassequence_plot(self): - da = xd.DataArray( - np.random.rand(100), - { - "time": { - "tie_indices": [0, 99], - "tie_values": [ - np.datetime64("2020-01-01"), - np.datetime64("2020-01-01T00:00:09.900"), - ], - } - }, - ) + da = xd.testing.dummy(dims=("time",), shape=(100,), step=0.01) dc = xd.DataCollection([da, da]) fig = xd.plot_availability(dc) assert fig is not None def test_datamapping_plot(self): - da = xd.DataArray( - np.random.rand(100), - { - "time": { - "tie_indices": [0, 99], - "tie_values": [ - np.datetime64("2020-01-01"), - np.datetime64("2020-01-01T00:00:09.900"), - ], - } - }, - ) + da = xd.testing.dummy(dims=("time",), shape=(100,), step=0.01) dm = xd.DataCollection({"a": da, "b": da}) fig = xd.plot_availability(dm) assert fig is not None diff --git a/tests/test_signal.py b/tests/test_signal.py index f8ea64a1..7d1bddaa 100644 --- a/tests/test_signal.py +++ b/tests/test_signal.py @@ -5,41 +5,19 @@ import xdas as xd import xdas.signal as xs -from xdas.synthetics import wavelet_wavefronts class TestSignal: def test_get_sample_spacing(self): - shape = (6000, 1000) - resolution = (np.timedelta64(8, "ms"), 5.0) - starttime = np.datetime64("2023-01-01T00:00:00") - da = xd.DataArray( - data=np.random.randn(*shape).astype("float32"), - coords={ - "time": { - "tie_indices": [0, shape[0] - 1], - "tie_values": [ - starttime, - starttime + resolution[0] * (shape[0] - 1), - ], - }, - "distance": { - "tie_indices": [0, shape[1] - 1], - "tie_values": [0.0, resolution[1] * (shape[1] - 1)], - }, - }, - ) - assert xs.get_sampling_interval(da, "time") == 0.008 - assert xs.get_sampling_interval(da, "distance") == 5.0 + da = xd.testing.dummy(shape=(6000, 1000), step=(0.008, 5.0), dtype="float32") + assert da.coords["time"].get_sampling_interval() == 0.008 + assert da.coords["distance"].get_sampling_interval() == 5.0 def test_deterend(self): - n = 100 - d = 5.0 - s = d * np.arange(n) - da = xr.DataArray(np.arange(n), {"time": s}) - da = xd.DataArray.from_xarray(da) + # dummy data is a linear ramp, so detrending must flatten it to zero + da = xd.testing.dummy(dims=("time",), shape=(100,), step=5.0, datetime=False) da = xs.detrend(da) - assert np.allclose(da.values, np.zeros(n)) + assert np.allclose(da.values, np.zeros(100)) def test_differentiate(self): n = 100 @@ -47,6 +25,7 @@ def test_differentiate(self): s = (d / 2) + d * np.arange(n) da = xr.DataArray(np.ones(n), {"distance": s}) da = xd.DataArray.from_xarray(da) + da["distance"] = da["distance"].to_regular() da = xs.differentiate(da, midpoints=True) assert np.allclose(da.values, np.zeros(n - 1)) @@ -56,6 +35,7 @@ def test_integrate(self): s = (d / 2) + d * np.arange(n) da = xr.DataArray(np.ones(n), {"distance": s}) da = xd.DataArray.from_xarray(da) + da["distance"] = da["distance"].to_regular() da = xs.integrate(da, midpoints=True) assert np.allclose(da.values, da["distance"].values) @@ -81,11 +61,12 @@ def test_sliding_window_removal(self): data = np.ones(n) da = xr.DataArray(data, {"distance": s}) da = xd.DataArray.from_xarray(da) + da["distance"] = da["distance"].to_regular() da = xs.sliding_mean_removal(da, 0.1 * n * d) assert np.allclose(da.values, 0) def test_medfilt(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result1 = xs.medfilt(da, {"distance": 3}) result2 = xs.medfilt(da, {"time": 1, "distance": 3}) assert result1.equals(result2) @@ -93,46 +74,46 @@ def test_medfilt(self): assert da.equals(xs.medfilt(da, {"time": 7, "distance": 3})) def test_hilbert(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = xs.hilbert(da, dim="time") assert np.allclose(da.values, np.real(result.values)) def test_resample(self): - da = wavelet_wavefronts() - result = xs.resample(da, 100, dim="time", window="hamming", domain="time") - assert result.sizes["time"] == 100 + da = xd.testing.dummy() + result = xs.resample(da, 50, dim="time", window="hamming", domain="time") + assert result.sizes["time"] == 50 def test_resample_poly(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = xs.resample_poly(da, 2, 5, dim="time") - assert result.sizes["time"] == 120 + assert result.sizes["time"] == 40 def test_lfilter(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() b, a = sp.iirfilter(4, 0.5, btype="low") result1 = xs.lfilter(b, a, da, "time") result2, zf = xs.lfilter(b, a, da, "time", zi=...) assert result1.equals(result2) def test_filtfilt(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() b, a = sp.iirfilter(2, 0.5, btype="low") xs.filtfilt(b, a, da, "time", padtype=None) def test_sosfilter(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() sos = sp.iirfilter(4, 0.5, btype="low", output="sos") result1 = xs.sosfilt(sos, da, "time") result2, zf = xs.sosfilt(sos, da, "time", zi=...) assert result1.equals(result2) def test_sosfiltfilt(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() sos = sp.iirfilter(2, 0.5, btype="low", output="sos") xs.sosfiltfilt(sos, da, "time", padtype=None) def test_filter(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() axis = da.get_axis_num("time") fs = 1 / xd.get_sampling_interval(da, "time") sos = sp.butter( @@ -168,7 +149,7 @@ def test_filter(self): assert result.equals(expected) def test_decimate_virtual_stack(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() expected = xs.decimate(da, 5, dim="time") chunks = xd.split(da, 5, "time") for i, chunk in enumerate(chunks): @@ -181,15 +162,7 @@ def test_decimate_virtual_stack(self, tmp_path): class TestSTFT: def test_compare_with_scipy(self): - starttime = np.datetime64("2023-01-01T00:00:00") - endtime = starttime + 9999 * np.timedelta64(10, "ms") - da = xd.DataArray( - data=np.random.rand(10000, 11), - coords={ - "time": {"tie_indices": [0, 9999], "tie_values": [starttime, endtime]}, - "distance": {"tie_indices": [0, 10], "tie_values": [0.0, 1.0]}, - }, - ) + da = xd.testing.dummy(shape=(10000, 11), step=(0.01, 0.1)) for scaling in ["spectrum", "psd"]: for return_onesided in [True, False]: for nfft in [None, 128]: @@ -205,7 +178,7 @@ def test_compare_with_scipy(self): ) f, t, Zxx = sp.stft( da.values, - fs=1 / xs.get_sampling_interval(da, "time"), + fs=1 / da.coords["time"].get_sampling_interval(), window="hamming", nperseg=100, noverlap=50, @@ -236,12 +209,10 @@ def test_retrieve_frequency_peak(self): N = 1e5 fc = 3e3 amp = 2 * np.sqrt(2) - time = np.arange(N) / float(fs) - data = amp * np.sin(2 * np.pi * fc * time) - da = xd.DataArray( - data=data, - coords={"time": time}, + da = xd.testing.dummy( + dims=("time",), shape=(int(N),), step=1 / fs, datetime=False ) + da.data = amp * np.sin(2 * np.pi * fc * da["time"].values) result = xs.stft( da, nperseg=1000, noverlap=500, window="hann", dim={"time": "frequency"} ) @@ -249,15 +220,7 @@ def test_retrieve_frequency_peak(self): assert result["frequency"][idx].values == fc def test_parrallel(self): - starttime = np.datetime64("2023-01-01T00:00:00") - endtime = starttime + 9999 * np.timedelta64(10, "ms") - da = xd.DataArray( - data=np.random.rand(10000, 11), - coords={ - "time": {"tie_indices": [0, 9999], "tie_values": [starttime, endtime]}, - "distance": {"tie_indices": [0, 10], "tie_values": [0.0, 1.0]}, - }, - ) + da = xd.testing.dummy(shape=(10000, 11), step=(0.01, 0.1)) serial = xs.stft( da, nperseg=100, @@ -277,16 +240,8 @@ def test_parrallel(self): assert serial.equals(parallel) def test_last_dimension_with_non_dimensional_coordinates(self): - starttime = np.datetime64("2023-01-01T00:00:00") - endtime = starttime + 99 * np.timedelta64(10, "ms") - da = xd.DataArray( - data=np.random.rand(100, 1001), - coords={ - "time": {"tie_indices": [0, 99], "tie_values": [starttime, endtime]}, - "distance": {"tie_indices": [0, 1000], "tie_values": [0.0, 10_000.0]}, - "channel": ("distance", np.arange(1001)), - }, - ) + da = xd.testing.dummy(shape=(100, 1001)) + da["channel"] = ("distance", np.arange(1001)) result = xs.stft( da, nperseg=100, @@ -296,7 +251,7 @@ def test_last_dimension_with_non_dimensional_coordinates(self): ) f, t, Zxx = sp.stft( da.values, - fs=1 / xs.get_sampling_interval(da, "distance"), + fs=1 / da.coords["distance"].get_sampling_interval(), window="hamming", nperseg=100, noverlap=50, @@ -313,41 +268,41 @@ def test_last_dimension_with_non_dimensional_coordinates(self): class TestSignalMissingBranches: def test_integrate_no_midpoints(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = xs.integrate(da, midpoints=False) assert result.shape == da.shape def test_differentiate_no_midpoints(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = xs.differentiate(da, midpoints=False) assert result.sizes["distance"] == da.sizes["distance"] - 1 def test_sliding_mean_removal_even_window(self): # When wlen/d gives an even n, sliding_mean_removal increments n by 1. - da = wavelet_wavefronts() - d = xs.get_sampling_interval(da, "time") + da = xd.testing.dummy() + d = da.coords["time"].get_sampling_interval() # Make wlen exactly twice d so n=2 (even) → becomes 3 result = xs.sliding_mean_removal(da, wlen=2 * d) assert result.shape == da.shape def test_medfilt_invalid_dim(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() with pytest.raises(ValueError, match="dims provided not in dataarray"): xs.medfilt(da, {"nonexistent_dim": 3}) def test_stft_default_noverlap(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = xs.stft(da, nperseg=16, dim={"time": "frequency"}) assert "frequency" in result.dims def test_stft_invalid_scaling(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() with pytest.raises(ValueError, match="Scaling must be"): xs.stft(da, nperseg=16, scaling="invalid", dim={"time": "frequency"}) def test_stft_nperseg_one(self): # nperseg=1, noverlap=0 triggers the stride_tricks bypass branch - da = wavelet_wavefronts() + da = xd.testing.dummy() # nfft=2 avoids single-element frequency axis (which would make tie_indices=[0,0]) result = xs.stft(da, nperseg=1, noverlap=0, nfft=2, dim={"time": "frequency"}) assert "frequency" in result.dims @@ -357,7 +312,7 @@ class TestFftMissingBranches: def test_fft_explicit_n(self): import xdas.fft as xfft - da = wavelet_wavefronts().isel(distance=0) + da = xd.testing.dummy().isel(distance=0) n = da.sizes["time"] // 2 result = xfft.fft(da, n=n, dim={"time": "frequency"}) assert result.sizes["frequency"] == n @@ -365,15 +320,22 @@ def test_fft_explicit_n(self): def test_rfft_explicit_n(self): import xdas.fft as xfft - da = wavelet_wavefronts().isel(distance=0) + da = xd.testing.dummy().isel(distance=0) n = da.sizes["time"] result = xfft.rfft(da, n=n, dim={"time": "frequency"}) assert "frequency" in result.dims + def test_rfft_single_frequency(self): + import xdas.fft as xfft + + da = xd.testing.dummy().isel(distance=0) + result = xfft.rfft(da, n=1, dim={"time": "frequency"}) + assert result.sizes["frequency"] == 1 + def test_ifft_explicit_n(self): import xdas.fft as xfft - da = wavelet_wavefronts().isel(distance=0) + da = xd.testing.dummy().isel(distance=0) spectrum = xfft.fft(da, dim={"time": "frequency"}) n = da.sizes["time"] result = xfft.ifft(spectrum, n=n, dim={"frequency": "time"}) diff --git a/tests/test_testing.py b/tests/test_testing.py new file mode 100644 index 00000000..42a1ec11 --- /dev/null +++ b/tests/test_testing.py @@ -0,0 +1,25 @@ +import numpy as np +import pytest + +import xdas as xd + + +class TestDummy: + def test_defaults(self): + da = xd.testing.dummy() + assert da.shape == (100, 10) + assert da.dims == ("time", "distance") + assert da["time"].isregular() + assert da["distance"].isregular() + + def test_mismatched_shape(self): + with pytest.raises(ValueError, match="must equal len\\(shape\\)"): + xd.testing.dummy(dims=("time",), shape=(10, 10)) + + def test_mismatched_step(self): + with pytest.raises(ValueError, match="must equal len\\(dims\\)"): + xd.testing.dummy(step=(1.0,)) + + def test_datetime_step_passthrough(self): + da = xd.testing.dummy(step=(np.timedelta64(10, "ms"), 10.0)) + assert da["time"].get_sampling_interval() == 0.01 diff --git a/tests/test_virtual.py b/tests/test_virtual.py index 00165758..21dd7777 100644 --- a/tests/test_virtual.py +++ b/tests/test_virtual.py @@ -3,7 +3,6 @@ import pytest import xdas as xd -from xdas.synthetics import wavelet_wavefronts from xdas.virtual import ( Selection, Selectors, @@ -19,7 +18,7 @@ class TestFunctional: # TODO: move elsewhere def test_all(self, tmp_path): - expected = wavelet_wavefronts() + expected = xd.testing.dummy() chunks = xd.split(expected, 3) for index, chunk in enumerate(chunks, start=1): chunk.to_netcdf(tmp_path / f"{index:03d}.nc") @@ -437,7 +436,7 @@ def test_check_dtype_mismatch(self, tmp_path): class TestVirtualLayoutExtra: def test_array_with_dtype(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() da.to_netcdf(tmp_path / "c.nc") da2 = xd.open(tmp_path / "c.nc") layout = da2.data._to_layout() @@ -445,7 +444,7 @@ def test_array_with_dtype(self, tmp_path): assert result.dtype == np.float32 def test_setitem_with_virtual_source(self, tmp_path): - da = wavelet_wavefronts() + da = xd.testing.dummy() da.to_netcdf(tmp_path / "d.nc") with h5py.File(tmp_path / "d.nc", "r") as f: src = VirtualSource(f["__values__"]) diff --git a/tests/test_xarray.py b/tests/test_xarray.py index 9e87a62a..9fc456da 100644 --- a/tests/test_xarray.py +++ b/tests/test_xarray.py @@ -2,12 +2,11 @@ import xdas as xd import xdas.core.methods as xm -from xdas.synthetics import wavelet_wavefronts class TestXarray: def test_returns_dataarray(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() for name, func in xm.HANDLED_METHODS.items(): if callable(func): if name in [ @@ -30,7 +29,7 @@ def test_returns_dataarray(self): assert isinstance(result, xd.DataArray) def test_mean(self): - da = wavelet_wavefronts() + da = xd.testing.dummy() result = xm.mean(da, "time") result_method = da.mean("time") expected = np.mean(da, 0) diff --git a/xdas/__init__.py b/xdas/__init__.py index d3accb1f..023abefb 100644 --- a/xdas/__init__.py +++ b/xdas/__init__.py @@ -24,6 +24,7 @@ "routines", "signal", "synthetics", + "testing", "virtual", # classes "Coordinate", @@ -67,6 +68,7 @@ processing, signal, synthetics, + testing, virtual, ) from .coordinates import ( @@ -78,11 +80,11 @@ ScalarCoordinate, get_sampling_interval, ) -from .core import dataarray, datacollection, methods, numpy, routines -from .core.dataarray import DataArray -from .core.datacollection import DataCollection, DataMapping, DataSequence -from .core.methods import * # noqa: F403 -from .core.routines import ( +from .core import ( + DataArray, + DataCollection, + DataMapping, + DataSequence, align, asdataarray, broadcast_coords, @@ -92,6 +94,10 @@ concat, concat_coords, concatenate, + dataarray, + datacollection, + methods, + numpy, open, open_dataarray, open_datacollection, @@ -99,5 +105,7 @@ open_mfdatacollection, open_mfdatatree, plot_availability, + routines, split, ) +from .core.methods import * # noqa: F403 diff --git a/xdas/atoms/core.py b/xdas/atoms/core.py index e3cea04c..7d18a46e 100644 --- a/xdas/atoms/core.py +++ b/xdas/atoms/core.py @@ -10,9 +10,7 @@ from functools import wraps from typing import Any -from ..core.dataarray import DataArray -from ..core.datacollection import DataCollection -from ..core.routines import open_datacollection +from ..core import DataArray, DataCollection, open_datacollection class State: diff --git a/xdas/atoms/ml.py b/xdas/atoms/ml.py index c953f921..c85e35be 100644 --- a/xdas/atoms/ml.py +++ b/xdas/atoms/ml.py @@ -8,8 +8,7 @@ import numpy as np -from ..core.dataarray import DataArray -from ..core.routines import concat +from ..core import DataArray, concat from .core import Atom, State diff --git a/xdas/atoms/signal.py b/xdas/atoms/signal.py index a077573d..2883f502 100644 --- a/xdas/atoms/signal.py +++ b/xdas/atoms/signal.py @@ -10,9 +10,8 @@ import numpy as np import scipy.signal as sp -from ..coordinates.core import Coordinate, get_sampling_interval -from ..core.dataarray import DataArray -from ..core.routines import concat, split +from ..coordinates import Coordinate, get_sampling_interval +from ..core import DataArray, concat, split from ..parallel import parallelize from .core import Atom, State @@ -569,14 +568,22 @@ def call(self, da, **flags): data[slc] = da.values coords = da.coords.copy() delta = get_sampling_interval(da, self.dim, cast=False) - tie_indices = coords[self.dim].tie_indices * self.factor - tie_values = coords[self.dim].tie_values + new_delta = delta / self.factor + coord = coords[self.dim] + tie_indices = coord.tie_indices * self.factor + tie_values = coord.tie_values tie_indices[-1] += self.factor - 1 - tie_values[-1] += (self.factor - 1) / self.factor * delta + tie_values[-1] += (self.factor - 1) * new_delta + # The derived rate may not be exactly representable (integer datetime + # resolutions truncate), so declare the representation error as jitter + # on top of the inherited one; chunk seams then stay within tolerance. + tolerance = coord.tolerance + np.abs(delta - new_delta * self.factor) coords[self.dim] = Coordinate( { "tie_indices": tie_indices, "tie_values": tie_values, + "sampling_interval": new_delta, + "tolerance": tolerance, }, self.dim, ) diff --git a/xdas/coordinates/__init__.py b/xdas/coordinates/__init__.py index 00a7989a..82ad4c7d 100644 --- a/xdas/coordinates/__init__.py +++ b/xdas/coordinates/__init__.py @@ -7,6 +7,7 @@ """ __all__ = [ + "AxisCoordinate", "Coordinate", "Coordinates", "DenseCoordinate", @@ -16,7 +17,12 @@ "get_sampling_interval", ] -from .core import Coordinate, Coordinates, get_sampling_interval +from .core import ( + AxisCoordinate, + Coordinate, + Coordinates, + get_sampling_interval, +) from .dense import DenseCoordinate from .interp import InterpCoordinate from .sampled import SampledCoordinate diff --git a/xdas/coordinates/core.py b/xdas/coordinates/core.py index 26ff51a3..acb09e0e 100644 --- a/xdas/coordinates/core.py +++ b/xdas/coordinates/core.py @@ -2,10 +2,11 @@ Core coordinate infrastructure. Includes the :class:`Coordinates` container, :class:`Coordinate` factory/base -class, and shared helpers used by all concrete coordinate types (parsing, -interpolation, tolerance handling). +class, :class:`AxisCoordinate` (the axis-mapping ABC), and shared helpers used by +all concrete coordinate types (parsing, interpolation, tolerance handling). """ +import warnings import weakref from abc import ABC, abstractmethod from copy import copy, deepcopy @@ -14,6 +15,19 @@ import numpy as np import pandas as pd +from typing_extensions import override + +#: Mapping from numpy datetime64/timedelta64 unit codes to CF-style unit names, +#: used to serialise timedelta scalars into dataset attributes. +CODE_TO_UNITS = { + "h": "hours", + "m": "minutes", + "s": "seconds", + "ms": "milliseconds", + "us": "microseconds", + "ns": "nanoseconds", +} +UNITS_TO_CODE = {v: k for k, v in CODE_TO_UNITS.items()} def wraps_first_last(func): @@ -107,7 +121,7 @@ def __setitem__(self, key, value): if not isinstance(key, str): raise TypeError("dimension names must be of type str") coord = Coordinate(value) - if coord.dim is None and not coord.isscalar(): + if coord.dim is None and isinstance(coord, AxisCoordinate): coord.dim = key if self.parent is None: if coord.dim is not None and coord.dim not in self.dims: @@ -278,20 +292,18 @@ class Coordinate(ABC): """ Base class and factory for all coordinate types. - A coordinate maps the integer positions of one array axis to physical - values (e.g. timestamps, distances). It supports two complementary - directions of lookup: + A coordinate attaches physical meaning to a :class:`DataArray`. Two kinds + exist: - - **Index-based selection** — ``coord[i]`` or ``coord[start:stop]``: - given integer position(s), return the corresponding physical value(s) - as a new coordinate. - - **Label-based selection** — ``coord.to_index(v)``: given a physical - value (or slice of values), return the integer index (or slice) at - that label. An optional *method* argument controls nearest/forward/ - backward matching for values that fall between samples. The returned - index can then be passed to ``coord[idx]`` to retrieve the - coordinate subset, and is also used internally to index into the - parent data array. + - **Axis coordinates** (:class:`AxisCoordinate` subclasses) map the integer + positions of one array axis to physical values (e.g. timestamps, + distances) and support index- and label-based selection. + - **Scalar coordinates** (:class:`ScalarCoordinate`) carry a single value + with no associated axis. + + This base class holds only what is genuinely shared between the two: the + factory/registry machinery, identity/equality, copying, and (de)serialisation + hooks. The full axis-mapping contract lives on :class:`AxisCoordinate`. **Factory behaviour** — calling ``Coordinate(data)`` directly acts as a factory: it inspects *data* and returns an instance of the most suitable @@ -333,7 +345,7 @@ def __new__(cls, data=None, dim=None, dtype=None): if data is None: raise TypeError("cannot infer coordinate type if no `data` is provided") - data, dim = parse(data, dim) + data, dim = parse_data_dim(data, dim) for subcls in Coordinate._registry.values(): if subcls._isvalid(data): @@ -351,6 +363,200 @@ def __new__(cls, data=None, dim=None, dtype=None): def __init__(self, data=None, dim=None, dtype=None): """Initialise the coordinate from subclass-specific *data*.""" + @property + @abstractmethod + def dtype(self): + """NumPy dtype of the underlying coordinate values.""" + + @staticmethod + @abstractmethod + def _isvalid(data): + """Return ``True`` if *data* is a valid input for this coordinate subclass.""" + + @property + @abstractmethod + def shape(self): + """Shape tuple of the coordinate (``()`` for scalar, ``(len(self),)`` for axis).""" + + @abstractmethod + def __array__(self, dtype=None, copy=None): + """Materialise this coordinate as a numpy array (numpy array protocol).""" + + @abstractmethod + def _to_dataset(self, dataset, attrs): + """ + Serialise this coordinate into an xarray *dataset*, updating *attrs* in place. + + Parameters + ---------- + dataset : xarray.Dataset + Target dataset to write coordinate data into. + attrs : dict + Global attribute mapping to update (e.g. ``coordinate_interpolation``). + + Returns + ------- + dataset : xarray.Dataset + attrs : dict + """ + + @classmethod + @abstractmethod + def _collect_from_dataset(cls, dataset, name): + """ + Extract coordinates of this subclass's type from *dataset* variable *name*. + + Parameters + ---------- + dataset : xarray.Dataset + Source dataset. + name : str + Name of the variable whose coordinates should be extracted. + + Returns + ------- + dict + Mapping from coordinate name to coordinate-like data, ready to be + passed to :class:`Coordinate`. + """ + + # -- properties --- + + #: Name of the dimension this coordinate is associated with, or ``None``. + dim = None + + @property + def size(self): + """Number of elements in this coordinate (``1`` for a scalar).""" + return int(np.prod(self.shape)) + + @property + def values(self): + """Materialised numpy array of coordinate values.""" + return self.__array__(copy=False) + + @property + def parent(self): + """The parent :class:`Coordinates` container, or ``None`` if unattached.""" + if hasattr(self, "_parent"): + return self._parent() + else: + return None + + @property + def name(self): + """The name under which this coordinate is stored in its parent container.""" + if self.parent is None: + return self.dim + return next((name for name in self.parent if self.parent[name] is self), None) + + # --- dunders logic --- + + def __reduce__(self): + return self.__class__, (self.data, self.dim) + + # --- queries --- + + def isdim(self): + """Return ``True`` if this coordinate is a dimensional coordinate.""" + if self.parent is None or self.name is None: + return None + else: + return self.parent.isdim(self.name) + + def isregular(self): + """ + Return ``True`` if this coordinate carries a nominal sampling interval. + + Scalar coordinates are never regular. Axis coordinates are regular when + :meth:`AxisCoordinate.get_sampling_interval` returns a value, i.e. when + an explicit nominal spacing is part of their data. + """ + return False + + def equals(self, other): + """Return ``True`` if *other* is the same coordinate type with identical dim and data. + + Comparison is strict on dtype. Same type implies same ``data`` structure: + either a single ``np.ndarray`` or a flat ``dict[str, np.ndarray]`` with + the same keys. + """ + if type(self) is not type(other) or self.dim != other.dim: + return False + a, b = self.data, other.data + if isinstance(a, dict): + pairs = [(a[key], b[key]) for key in a] + else: + pairs = [(a, b)] + for x, y in pairs: + x, y = np.asarray(x), np.asarray(y) + if x.dtype != y.dtype or not np.array_equal(x, y, equal_nan=False): + return False + return True + + # --- routines --- + + def copy(self, deep=True): + """ + Return a copy of this coordinate. + + Parameters + ---------- + deep : bool, optional + If ``True`` (default) perform a deep copy; otherwise a shallow copy. + + Returns + ------- + Coordinate + A new coordinate of the same subclass with copied data and metadata. + """ + if deep: + func = deepcopy + else: + func = copy + return self.__class__(func(self.data), func(self.dim), func(self.dtype)) + + # --- IO --- + + @classmethod + def _from_dataset(cls, dataset, name): + """Read coordinates named *name* from an xarray *dataset* via each registered subclass.""" + coords = {} + for subcls in cls._registry.values(): + coords |= subcls._collect_from_dataset(dataset, name) + return coords + + # --- internals --- + + def _assign_parent(self, parent): + """Attach this coordinate to its parent :class:`Coordinates` container.""" + self._parent = weakref.ref(parent) + + +class AxisCoordinate(Coordinate, ABC): + """ + Base class for coordinates that map an array axis to physical values. + + Adds the full axis-mapping contract on top of :class:`Coordinate`: it + supports two complementary directions of lookup: + + - **Index-based selection** — ``coord[i]`` or ``coord[start:stop]``: + given integer position(s), return the corresponding physical value(s) + as a new coordinate. + - **Label-based selection** — ``coord.to_index(v)``: given a physical + value (or slice of values), return the integer index (or slice) at + that label. An optional *method* argument controls nearest/forward/ + backward matching for values that fall between samples. The returned + index can then be passed to ``coord[idx]`` to retrieve the + coordinate subset, and is also used internally to index into the + parent data array. + + Concrete subclasses are :class:`DenseCoordinate`, :class:`InterpCoordinate`, + and :class:`SampledCoordinate`. + """ + + # --- abstract contract --- + @classmethod @abstractmethod def from_block(cls, start, size, step, dim=None, dtype=None): @@ -380,16 +586,6 @@ def from_block(cls, start, size, step, dim=None, dtype=None): def __len__(self): """Return the number of elements along this coordinate's axis.""" - @property - @abstractmethod - def dtype(self): - """NumPy dtype of the underlying coordinate values.""" - - @staticmethod - @abstractmethod - def _isvalid(data): - """Return ``True`` if *data* is a valid input for this coordinate subclass.""" - @abstractmethod def _is_monotonic_increasing(self): """Return ``True`` if all consecutive differences in this coordinate are positive.""" @@ -457,62 +653,126 @@ def _slice(self, slc): @abstractmethod def _concat(self, other): """ - Return a new coordinate formed by appending *other* after this one. + Return a new coordinate formed by appending *other* after this one. Parameters ---------- - other : Coordinate - Must be the same subclass and have the same ``dim`` and ``dtype``. + other : Coordinate + Must be the same subclass and have the same ``dim`` and ``dtype``. Returns ------- - Coordinate - Concatenated coordinate of the same subclass. - s + Coordinate + Concatenated coordinate of the same subclass. """ @abstractmethod - def _to_dataset(self, dataset, attrs): + def get_sampling_interval(self, cast=True): """ - Serialise this coordinate into an xarray *dataset*, updating *attrs* in place. + Return the nominal sample spacing for this coordinate, or ``None``. Parameters ---------- - dataset : xarray.Dataset - Target dataset to write coordinate data into. - attrs : dict - Global attribute mapping to update (e.g. ``coordinate_interpolation``). + cast : bool, optional + If ``True`` (default), cast timedelta64 results to seconds (float). Returns ------- - dataset : xarray.Dataset - attrs : dict + float or None + ``None`` if the coordinate has fewer than two elements or has no + defined sampling interval. """ - @classmethod @abstractmethod - def _collect_from_dataset(cls, dataset, name): + def to_regular(self, sampling_interval=None, tolerance=None): """ - Extract coordinates of this subclass's type from *dataset* variable *name*. + Return a regular version of this coordinate, raising when impossible. + + The strict conversion entry point: the result always satisfies + :meth:`isregular` (it carries a nominal ``sampling_interval``), or a + :exc:`ValueError` is raised when the coordinate values cannot be + described by a single spacing within *tolerance*. Parameters ---------- - dataset : xarray.Dataset - Source dataset. - name : str - Name of the variable whose coordinates should be extracted. + sampling_interval : scalar, optional + Nominal sample spacing to enforce. When omitted it is taken from + the coordinate itself when available, or inferred from the values. + tolerance : scalar, optional + Tolerated jitter around *sampling_interval*. Defaults to the + coordinate's declared tolerance when present, else a zero-like + default (exact zero for datetime axes, a dtype epsilon for floats), + so a genuinely irregular axis raises. Returns ------- - dict - Mapping from coordinate name to coordinate-like data, ready to be - passed to :class:`Coordinate`. + AxisCoordinate + A regular coordinate. The subclass may change: a + :class:`DenseCoordinate` converts to a regular + :class:`InterpCoordinate`. """ - # -- properties --- + @abstractmethod + def _split_candidates(self): + """ + Return the candidate segment boundaries used by :meth:`get_split_indices`. - #: Name of the dimension this coordinate is associated with, or ``None``. - dim = None + Returns + ------- + positions : numpy.ndarray + Integer index of each candidate boundary. Each ``positions[k]`` + marks the start of a new segment (the boundary lies between element + ``positions[k] - 1`` and ``positions[k]``). + deltas : numpy.ndarray + Signed jump at each candidate, i.e. the value step across the + boundary minus the nominal sampling interval. Positive values are + gaps, negative values are overlaps, zero means a clean continuation. + """ + + @abstractmethod + def simplify(self, tolerance=None, *, reduce=True, regularize=False): + """ + Return the simplest faithful copy of this coordinate within *tolerance*. + + A coordinate carries two independent, orthogonal properties: + + - **Monotonic** — tie values strictly increase, which is what makes + label-based selection (:meth:`to_index`) work. + - **Regular** — every continuous segment fits a single nominal + ``sampling_interval``, which signal-processing routines (FFT, + filtering, resampling) need for a clean sample rate. + + ``simplify`` spends an accuracy budget *tolerance* across two + independently toggleable stages: + + - *reduce* drops redundant points whose removal shifts the curve by no + more than *tolerance* (which also absorbs soft gaps and overlaps, + helping monotonicity). Surviving values are never moved. + - *regularize* promotes the coordinate to *regular* when the surviving + continuous segments admit a single spacing within *tolerance*. + + Parameters + ---------- + tolerance : float, timedelta, None, or ``False``, optional + Accuracy budget; maximum allowed deviation from the original + values. ``None`` (default) spends the coordinate's own declared + tolerance when it carries one, else a zero-like default (exact + zero for datetime axes, a dtype epsilon for floats). ``False`` + returns an unchanged copy regardless of the flags below. + reduce : bool, optional + Whether to drop redundant tie points. Default ``True``. + regularize : bool, optional + Whether to try to acquire a nominal ``sampling_interval``. Default + ``False`` (opt-in). A no-op for coordinates that are already regular + by construction or cannot become regular. + + Returns + ------- + Coordinate + A new coordinate of the same subclass. + """ + + # --- properties --- @property def ndim(self): @@ -524,11 +784,6 @@ def shape(self): """Shape tuple ``(len(self),)``.""" return (len(self),) - @property - def size(self): - """Number of elements along this coordinate's axis.""" - return len(self) - @property def empty(self): """``True`` if the coordinate has zero length.""" @@ -539,11 +794,6 @@ def indices(self): """Integer array ``[0, 1, ..., len(self) - 1]``.""" return np.arange(len(self)) - @property - def values(self): - """Materialised numpy array of coordinate values.""" - return self.__array__(copy=False) - @property def start(self): """Value at index 0 (first element).""" @@ -554,21 +804,6 @@ def end(self): """Value at the last element.""" return self._get_value(len(self) - 1) - @property - def parent(self): - """The parent :class:`Coordinates` container, or ``None`` if unattached.""" - if hasattr(self, "_parent"): - return self._parent() - else: - return None - - @property - def name(self): - """The name under which this coordinate is stored in its parent container.""" - if self.parent is None: - return self.dim - return next((name for name in self.parent if self.parent[name] is self), None) - # --- dunders logic --- def __getitem__(self, item): @@ -589,9 +824,6 @@ def __array__(self, dtype=None, copy=None): out = out.__array__(dtype) return out - def __reduce__(self): - return self.__class__, (self.data, self.dim) - def __repr__(self): if self.empty: return "empty coordinate" @@ -609,36 +841,178 @@ def __repr__(self): # --- queries --- - def isscalar(self): - """Return ``True`` if this is a :class:`ScalarCoordinate`.""" - return False + @override + def isregular(self): + return self.get_sampling_interval() is not None - def isdim(self): - """Return ``True`` if this coordinate is a dimensional coordinate.""" - if self.parent is None or self.name is None: - return None - else: - return self.parent.isdim(self.name) + def get_split_indices(self, kind="discontinuities", tolerance=False): + """ + Return integer indices where this coordinate should be split. - def equals(self, other): - """Return ``True`` if *other* is the same coordinate type with identical dim and data. + Each returned index ``i`` marks the start of a new segment: the + boundary lies between element ``i - 1`` and element ``i``. The first + segment always starts at index 0, so 0 is never included in the result. - Comparison is strict on dtype. Same type implies same ``data`` structure: - either a single ``np.ndarray`` or a flat ``dict[str, np.ndarray]`` with - the same keys. + Parameters + ---------- + kind : {"discontinuities", "gaps", "overlaps"}, optional + Which boundary type to return. ``"gaps"`` returns only boundaries + where the axis jumps forward by more than one sampling interval; + ``"overlaps"`` returns only boundaries where the axis jumps + backward. ``"discontinuities"`` (default) returns both. + tolerance : float, timedelta, None, or ``False``, optional + Minimum absolute magnitude of the jump to report. Boundaries + smaller than *tolerance* are silently dropped. ``None`` removes + only zero-magnitude jumps (i.e. consecutive equal values). + ``False`` (default) disables magnitude filtering and returns all + boundaries of the requested kind. + + Returns + ------- + numpy.ndarray + Integer indices of the start of each new segment (excluding the first). """ - if type(self) is not type(other) or self.dim != other.dim: - return False - a, b = self.data, other.data - if isinstance(a, dict): - pairs = [(a[key], b[key]) for key in a] + valid_kinds = {"discontinuities", "gaps", "overlaps"} + if kind not in valid_kinds: + raise ValueError(f"`kind` must be one of {valid_kinds}; got {kind!r}") + + positions, deltas = self._split_candidates() + + # Fast path: every candidate boundary is a discontinuity by construction + if kind == "discontinuities" and tolerance is False: + return positions + + if tolerance is False: + zero = np.timedelta64(0) if np.issubdtype(self.dtype, np.datetime64) else 0 + match kind: + case "gaps": + mask = deltas >= zero + case "overlaps": # pragma: no branch + mask = deltas < zero else: - pairs = [(a, b)] - for x, y in pairs: - x, y = np.asarray(x), np.asarray(y) - if x.dtype != y.dtype or not np.array_equal(x, y, equal_nan=False): - return False - return True + tolerance = parse_scalar_delta(tolerance, self.dtype, default_zero=True) + match kind: + case "discontinuities": + mask = np.abs(deltas) > tolerance + case "gaps": + mask = deltas > tolerance + case "overlaps": # pragma: no branch + mask = deltas < -tolerance + + return positions[mask] + + def get_discontinuities(self, tolerance=None): + """ + Return a DataFrame containing information about the discontinuities. + + Parameters + ---------- + tolerance : float, timedelta, or None, optional + Minimum magnitude of a gap or overlap to include. ``None`` + (default) reports all discontinuities regardless of size. + + Returns + ------- + pandas.DataFrame + A DataFrame with the following columns: + + - start_index : int + The index where the discontinuity starts. + - end_index : int + The index where the discontinuity ends. + - start_value : float + The value at the start of the discontinuity. + - end_value : float + The value at the end of the discontinuity. + - delta : float + The difference between the end_value and start_value. + - type : str + The type of the discontinuity, either "gap" or "overlap". + + """ + if self.empty: + return pd.DataFrame( + columns=[ + "start_index", + "end_index", + "start_value", + "end_value", + "delta", + "type", + ] + ) + indices = self.get_split_indices("discontinuities", tolerance) + records = [] + for index in indices: + start_index = index + end_index = index + 1 + start_value = self._get_value(index) + end_value = self._get_value(index + 1) + delta = end_value - start_value + if tolerance is not None and np.abs(delta) < tolerance: + continue + record = { + "start_index": start_index, + "end_index": end_index, + "start_value": start_value, + "end_value": end_value, + "delta": delta, + "type": ("gap" if end_value > start_value else "overlap"), + } + records.append(record) + return pd.DataFrame.from_records(records) + + def get_availabilities(self): + """ + Return a DataFrame containing information about the data availability. + + Returns + ------- + pandas.DataFrame + A DataFrame with the following columns: + + - start_index : int + The index where the discontinuity starts. + - end_index : int + The index where the discontinuity ends. + - start_value : float + The value at the start of the discontinuity. + - end_value : float + The value at the end of the discontinuity. + - delta : float + The difference between the end_value and start_value. + - type : str + The type of the discontinuity, always "data". + + """ + if self.empty: + return pd.DataFrame( + columns=[ + "start_index", + "end_index", + "start_value", + "end_value", + "delta", + "type", + ] + ) + indices = np.concatenate([[0], self.get_split_indices(), [len(self)]]) + records = [] + for start_index, stop_index in pairwise(indices): + end_index = stop_index - 1 + start_value = self._get_value(start_index) + end_value = self._get_value(end_index) + records.append( + { + "start_index": start_index, + "end_index": end_index, + "start_value": start_value, + "end_value": end_value, + "delta": end_value - start_value, + "type": "data", + } + ) + return pd.DataFrame.from_records(records) # --- selection / indexing --- @@ -762,29 +1136,9 @@ def _slice_indexer(self, start=None, stop=None, step=None, endpoint=True): # --- routines --- - def copy(self, deep=True): - """ - Return a copy of this coordinate. - - Parameters - ---------- - deep : bool, optional - If ``True`` (default) perform a deep copy; otherwise a shallow copy. - - Returns - ------- - Coordinate - A new coordinate of the same subclass with copied data and metadata. - """ - if deep: - func = deepcopy - else: - func = copy - return self.__class__(func(self.data), func(self.dim), func(self.dtype)) - def to_dataarray(self): """Convert this coordinate to a :class:`~xdas.DataArray` with a single dimension.""" - from ..core.dataarray import DataArray # TODO: avoid defered import? + from ..core import DataArray # TODO: avoid deferred import? if self.name is None: raise ValueError("cannot convert unnamed coordinate to DataArray") @@ -808,217 +1162,8 @@ def to_dataarray(self): name=self.name, ) - # --- IO --- - - @classmethod - def _from_dataset(cls, dataset, name): - """Read coordinates named *name* from an xarray *dataset* via each registered subclass.""" - coords = {} - for subcls in cls.__subclasses__(): - coords |= subcls._collect_from_dataset(dataset, name) - return coords - - # --- internals --- - - def _assign_parent(self, parent): - """Attach this coordinate to its parent :class:`Coordinates` container.""" - self._parent = weakref.ref(parent) - -class SampledMixin(ABC): - """ - Shared behaviour for coordinates that carry sampled values along an axis. - - Mixed into the tie-point coordinate types (:class:`SampledCoordinate`, - :class:`InterpCoordinate`). Both types describe a piecewise-monotonic axis - composed of contiguous segments separated by *gaps* (the axis jumps forward - by more than one sampling interval) or *overlaps* (the axis jumps backward, - creating doubly-covered regions). This mixin provides the shared logic for - detecting, cataloguing, and querying those discontinuities. - """ - - @abstractmethod - def get_sampling_interval(self, cast=True): - """ - Return the nominal sample spacing for this coordinate. - - Parameters - ---------- - cast : bool, optional - If ``True`` (default), cast timedelta64 results to seconds (float). - - Returns - ------- - float or None - ``None`` if the coordinate has fewer than two elements. - """ - - @abstractmethod - def get_split_indices(self, kind="discontinuities", tolerance=False): - """ - Return integer indices where this coordinate should be split. - - Each returned index ``i`` marks the start of a new segment: the - boundary lies between element ``i - 1`` and element ``i``. The first - segment always starts at index 0, so 0 is never included in the result. - - Parameters - ---------- - kind : {"discontinuities", "gaps", "overlaps"}, optional - Which boundary type to return. ``"gaps"`` returns only boundaries - where the axis jumps forward by more than one sampling interval; - ``"overlaps"`` returns only boundaries where the axis jumps - backward. ``"discontinuities"`` (default) returns both. - tolerance : float, timedelta, None, or ``False``, optional - Minimum absolute magnitude of the jump to report. Boundaries - smaller than *tolerance* are silently dropped. ``None`` removes - only zero-magnitude jumps (i.e. consecutive equal values). - ``False`` (default) disables magnitude filtering and returns all - boundaries of the requested kind. - - Returns - ------- - numpy.ndarray - Integer indices of the start of each new segment (excluding the first). - """ - - @abstractmethod - def simplify(self, tolerance=None): - """ - Return a simplified copy of this coordinate with redundant tie points removed. - - Tie points whose removal would shift any label by no more than *tolerance* - are dropped, reducing memory and I/O cost without meaningfully changing - the represented axis. As a side effect, small gaps or overlaps that fall - within *tolerance* may be absorbed, merging adjacent segments into one. - - Parameters - ---------- - tolerance : float, timedelta, None, or ``False``, optional - Maximum allowed deviation from the original values. ``None`` uses - zero tolerance (lossless). ``False`` returns an unchanged copy. - - Returns - ------- - Coordinate - A new coordinate of the same subclass with fewer stored points. - """ - - def get_discontinuities(self, tolerance=None): - """ - Return a DataFrame containing information about the discontinuities. - - Parameters - ---------- - tolerance : float, timedelta, or None, optional - Minimum magnitude of a gap or overlap to include. ``None`` - (default) reports all discontinuities regardless of size. - - Returns - ------- - pandas.DataFrame - A DataFrame with the following columns: - - - start_index : int - The index where the discontinuity starts. - - end_index : int - The index where the discontinuity ends. - - start_value : float - The value at the start of the discontinuity. - - end_value : float - The value at the end of the discontinuity. - - delta : float - The difference between the end_value and start_value. - - type : str - The type of the discontinuity, either "gap" or "overlap". - - """ - if self.empty: - return pd.DataFrame( - columns=[ - "start_index", - "end_index", - "start_value", - "end_value", - "delta", - "type", - ] - ) - indices = self.get_split_indices("discontinuities", tolerance) - records = [] - for index in indices: - start_index = index - end_index = index + 1 - start_value = self._get_value(index) - end_value = self._get_value(index + 1) - delta = end_value - start_value - if tolerance is not None and np.abs(delta) < tolerance: - continue - record = { - "start_index": start_index, - "end_index": end_index, - "start_value": start_value, - "end_value": end_value, - "delta": delta, - "type": ("gap" if end_value > start_value else "overlap"), - } - records.append(record) - return pd.DataFrame.from_records(records) - - def get_availabilities(self): - """ - Return a DataFrame containing information about the data availability. - - Returns - ------- - pandas.DataFrame - A DataFrame with the following columns: - - - start_index : int - The index where the discontinuity starts. - - end_index : int - The index where the discontinuity ends. - - start_value : float - The value at the start of the discontinuity. - - end_value : float - The value at the end of the discontinuity. - - delta : float - The difference between the end_value and start_value. - - type : str - The type of the discontinuity, always "data". - - """ - if self.empty: - return pd.DataFrame( - columns=[ - "start_index", - "end_index", - "start_value", - "end_value", - "delta", - "type", - ] - ) - indices = np.concatenate([[0], self.get_split_indices(), [len(self)]]) - records = [] - for start_index, stop_index in pairwise(indices): - end_index = stop_index - 1 - start_value = self._get_value(start_index) - end_value = self._get_value(end_index) - records.append( - { - "start_index": start_index, - "end_index": end_index, - "start_value": start_value, - "end_value": end_value, - "delta": end_value - start_value, - "type": "data", - } - ) - return pd.DataFrame.from_records(records) - - -def parse(data, dim=None): +def parse_data_dim(data, dim=None): """ Normalise *data* / *dim* inputs accepted by coordinate constructors. @@ -1052,38 +1197,85 @@ def parse(data, dim=None): return data, dim -def parse_tolerance(tolerance, dtype): +def parse_scalar_delta(value, dtype, default_zero=False): """ - Normalise *tolerance* to the correct type for *dtype*. + Normalise a scalar *value* to the correct type for *dtype*. - Converts ``None`` to zero, and for datetime64 dtypes converts a - numeric tolerance (in seconds) to the appropriate :class:`numpy.timedelta64`. + When ``default_zero`` is ``True``, a ``None`` input is replaced by a + sensible zero-like default: ``timedelta64(0)`` for datetime64 dtypes, + a dtype-appropriate epsilon for floating-point dtypes, or plain ``0`` + otherwise. For datetime64 dtypes a numeric value is interpreted as + seconds and converted to :class:`numpy.timedelta64`. Parameters ---------- - tolerance : float or None - Raw tolerance value. + value : scalar or None + Raw scalar value to normalise. dtype : numpy.dtype - The dtype of the coordinate values the tolerance will be compared against. + Target dtype that determines the output type. + default_zero : bool, optional + If ``True``, replace ``None`` with a zero-like default for *dtype*. + Default is ``False``. Returns ------- - tolerance : int, float, or numpy.timedelta64 + value : numpy scalar + Normalised scalar cast to the appropriate numpy scalar type. + + Raises + ------ + ValueError + If *value* is not a scalar (i.e. has non-zero ndim), or if *value* is + ``None`` while *default_zero* is ``False`` (no default is available). """ + # check shape + if not np.ndim(value) == 0: + raise ValueError("`value` must be a scalar value") + + # default + if value is None: + if not default_zero: + raise ValueError("`value` cannot be None when `default_zero` is False") + if np.issubdtype(dtype, np.datetime64): + value = np.timedelta64(0) + elif dtype == np.float16: + value = 1e-2 + elif dtype == np.float32: + value = 1e-5 + elif dtype == np.float64: + value = 1e-8 + else: + value = 0 + + # ensure numpy scalar + value = np.asarray(value)[()] + + # check dtype if np.issubdtype(dtype, np.datetime64): - if tolerance is None: - tolerance = np.timedelta64(0) - elif isinstance(tolerance, (int, float)): - tolerance = np.timedelta64(round(tolerance * 1e9), "ns") + if not np.issubdtype(value.dtype, np.timedelta64): + value = np.timedelta64(round(value * 1e9), "ns") # TODO: not `dtype` else: - if tolerance is None: - tolerance = 0 - return tolerance + value = value.astype(dtype) + + return value def get_sampling_interval(da, dim, cast=True): """ - Return the sample spacing along a given dimension. + Return the nominal sample spacing along a given dimension. + + Convenience used by every signal-processing routine: the coordinate should + be regular (carry a nominal sampling interval). Convert an irregular + coordinate first, e.g. ``da[dim] = da[dim].to_regular(tolerance=...)``, or + open the files with a tolerance so gaps and jitter are absorbed upfront. + + .. deprecated:: 0.2.8 + For backward compatibility with data saved by earlier versions (whose + coordinates carry no ``sampling_interval``), an irregular coordinate + currently falls back to inferring a spacing and emits a + :exc:`FutureWarning` stating the inferred value and the tolerance it + requires. This fallback will be removed in a future release, after + which irregular coordinates will raise. Parameters ---------- @@ -1096,17 +1288,89 @@ def get_sampling_interval(da, dim, cast=True): Returns ------- - float - The sample spacing. - - """ - return da[dim].get_sampling_interval(cast=cast) + float or None + The sample spacing. ``None`` when *dim* has no axis coordinate. + Raises + ------ + ValueError + If the coordinate is not regular and no spacing can be inferred. -def isscalar(data): - """Return ``True`` if *data* converts to a 0-d non-object numpy array.""" - data = np.asarray(data) - return (data.dtype != np.dtype(object)) and (data.ndim == 0) + """ + from .interp import InterpCoordinate # avoid circular import + + coord = da[dim] + if not isinstance(coord, AxisCoordinate): + return None + delta = coord.get_sampling_interval(cast=cast) + if delta is not None: + return delta + + # Deprecated fallback: data written by earlier versions carries no + # sampling_interval metadata, so infer one rather than break every + # signal-processing call on existing archives. + hint = ( + f"make the coordinate regular with `da[{dim!r}] = " + f"da[{dim!r}].to_regular(tolerance=...)`, or open the files with a " + f"tolerance" + ) + if isinstance(coord, InterpCoordinate): + sampling_interval, tolerance = coord._infer_regular() + else: + try: + regular = coord.to_regular() + except ValueError as exc: + raise ValueError( + f"coordinate {dim!r} has no nominal sampling interval and " + f"none could be inferred ({exc}); {hint}" + ) from exc + sampling_interval = regular.get_sampling_interval(cast=False) + tolerance = regular.tolerance if sampling_interval is not None else None + if sampling_interval is None: + raise ValueError( + f"coordinate {dim!r} has no nominal sampling interval and none " + f"could be inferred; {hint}" + ) + warnings.warn( + f"coordinate {dim!r} has no declared sampling interval; inferred " + f"{sampling_interval} (accepting jitter up to tolerance={tolerance}). " + f"This implicit inference is deprecated and will raise in a future " + f"release; {hint}", + FutureWarning, + stacklevel=2, + ) + if cast and np.issubdtype(np.asarray(sampling_interval).dtype, np.timedelta64): + sampling_interval = sampling_interval / np.timedelta64(1, "s") + return sampling_interval + + +def encode_delta(key, value): + """Serialise a scalar (possibly timedelta64) into a dict of dataset attributes.""" + if value is None: + return {} + if np.issubdtype(np.asarray(value).dtype, np.timedelta64): + code, count = np.datetime_data(value.dtype) + if code == "generic": # e.g. timedelta64(0); promote to nanoseconds + value = value.astype("timedelta64[ns]") + code, count = np.datetime_data(value.dtype) + return { + key: int(count * value.astype(int)), + f"{key}_dtype": "timedelta64[ns]", + f"{key}_units": CODE_TO_UNITS[code], + } + return {key: value} + + +def decode_delta(key, attrs): + """Inverse of :func:`encode_delta`: read a scalar back from dataset attributes.""" + if key not in attrs: + return None + value = attrs[key] + if f"{key}_units" in attrs: + value = np.timedelta64(value, UNITS_TO_CODE[attrs[f"{key}_units"]]).astype( + attrs[f"{key}_dtype"] + ) + return value def is_monotonic_increasing(x): diff --git a/xdas/coordinates/dense.py b/xdas/coordinates/dense.py index de95e72d..a09ed3cc 100644 --- a/xdas/coordinates/dense.py +++ b/xdas/coordinates/dense.py @@ -4,10 +4,11 @@ import pandas as pd from typing_extensions import override -from .core import Coordinate, parse +from .core import AxisCoordinate, parse_data_dim, parse_scalar_delta +from .interp import InterpCoordinate -class DenseCoordinate(Coordinate, ctype="dense"): +class DenseCoordinate(AxisCoordinate, ctype="dense"): """ Coordinate backed by an explicit numpy array. @@ -31,7 +32,7 @@ def __init__(self, data=None, dim=None, dtype=None): data = [] # parse data - data, dim = parse(data, dim) + data, dim = parse_data_dim(data, dim) if not self._isvalid(data): raise TypeError("`data` must be array-like") @@ -143,38 +144,66 @@ def __add__(self, other): def __sub__(self, other): return self.__class__(self.data - other, self.dim) + @override def get_sampling_interval(self, cast=True): """ - Return the average sample spacing (end-to-end distance divided by N-1). + Return ``None``: a dense coordinate never carries a nominal spacing. - Parameters - ---------- - cast : bool, optional - If ``True`` (default), cast timedelta64 results to seconds (float). + The raw values may happen to be evenly spaced, but regularity is an + explicit declaration; convert with :meth:`to_regular` to obtain a + regular :class:`InterpCoordinate`. + """ + return None - Returns - ------- - float or None - ``None`` if the coordinate has fewer than two elements. + @override + def to_regular(self, sampling_interval=None, tolerance=None): + """Convert to a regular :class:`InterpCoordinate` (single continuous ramp). + + The spacing defaults to the end-to-end slope, and every value must lie + within *tolerance* of the regular grid anchored at the first value; + otherwise a :exc:`ValueError` is raised. See + :meth:`AxisCoordinate.to_regular` for the parameter contract. """ if len(self) < 2: - return None - delta = (self[-1].values - self[0].values) / (len(self) - 1) - delta = np.asarray( - delta - ) # plain Python floats have no .dtype; np.asarray adds it - if cast and np.issubdtype(delta.dtype, np.timedelta64): - delta = delta / np.timedelta64(1, "s") - return delta - - def get_div_points(self, tolerance=None): - """Return sorted split-point indices where consecutive differences exceed *tolerance*.""" - deltas = np.diff(self.data) - if tolerance is not None: - div_points = np.nonzero(np.abs(deltas) >= tolerance)[0] + 1 + raise ValueError( + "cannot make a regular coordinate from fewer than two values" + ) + tolerance = parse_scalar_delta(tolerance, self.dtype, default_zero=True) + if sampling_interval is None: + sampling_interval = (self.data[-1] - self.data[0]) / (len(self) - 1) else: - raise NotImplementedError( - "get_div_points without tolerance is not implemented for DenseCoordinate" + sampling_interval = parse_scalar_delta(sampling_interval, self.dtype) + grid = self.data[0] + sampling_interval * np.arange(len(self)) + if not np.all(np.abs(self.data - grid) <= tolerance): + raise ValueError( + "values are not evenly spaced by `sampling_interval` within `tolerance`" ) - div_points = np.concatenate(([0], div_points, [len(self)])) - return div_points + data = { + "tie_indices": [0, len(self) - 1], + "tie_values": [self.data[0], self.data[-1]], + "sampling_interval": sampling_interval, + "tolerance": tolerance, + } + return InterpCoordinate(data, self.dim) + + @override + def _split_candidates(self): + steps = np.diff(self.data) + positions = np.arange(1, len(self)) + if steps.size == 0: + return positions, steps + reference = np.median(steps) + deltas = np.empty(steps.shape, dtype=np.asarray(steps[0] - reference).dtype) + for i in range(steps.size): + deltas[i] = steps[i] - reference + if i + 1 < steps.size and abs(steps[i + 1] - steps[i]) < abs( + steps[i + 1] - reference + ): + reference = steps[i] + return positions, deltas + + @override + def simplify(self, tolerance=None, *, reduce=True, regularize=False): + # a dense coordinate stores every value explicitly; there is nothing to + # drop and no spacing to promote, so both stages are no-ops. + return self.copy() diff --git a/xdas/coordinates/interp.py b/xdas/coordinates/interp.py index a6bfbde1..5b87d116 100644 --- a/xdas/coordinates/interp.py +++ b/xdas/coordinates/interp.py @@ -2,30 +2,44 @@ :class:`InterpCoordinate`: piecewise-linear coordinate. Defined by tie points, using ``xinterp`` for forward and inverse interpolation. +Optionally carries a nominal ``sampling_interval`` (and ``tolerance``) making the +coordinate *regular* and providing a clean sample rate for signal-processing routines. """ import re import numpy as np +from numba import njit from typing_extensions import override from xinterp import forward, inverse from .core import ( + AxisCoordinate, Coordinate, - SampledMixin, + decode_delta, + encode_delta, is_monotonic_increasing, - parse, - parse_tolerance, + parse_data_dim, + parse_scalar_delta, ) -class InterpCoordinate(SampledMixin, Coordinate, ctype="interpolated"): +class InterpCoordinate(AxisCoordinate, ctype="interpolated"): """ - Piecewise-linear coordinate described by tie points (CF convention). + Piecewise-linear coordinate described by tie points (CF subsampling, 8.3). - Values between tie points are recovered by linear interpolation. - Discontinuities are represented by two consecutive tie points at adjacent - indices. Supports label-based selection via :meth:`~Coordinate.to_index`. + Following the CF conventions for compression by coordinate subsampling. + Values between tie points are recovered by linear interpolation (via + ``xinterp``), which also enables label-based selection through + :meth:`~Coordinate.to_index`. The index axis is split into *continuous + areas* separated by *discontinuities*; a discontinuity is encoded as two + consecutive tie points at adjacent indices (a gap of one). + + When *data* contains a ``sampling_interval`` key the coordinate also + enforces a nominal sample spacing, making it *regular* + (:meth:`isregular` returns ``True``) and giving signal-processing routines a + clean sample rate. A ``tolerance`` key may accompany it to allow bounded + jitter around that rate. Parameters ---------- @@ -36,11 +50,30 @@ class InterpCoordinate(SampledMixin, Coordinate, ctype="interpolated"): ``tie_values`` : sequence of float or datetime64 Values at the tie points. Must be strictly increasing to enable label-based selection. Length must match ``tie_indices``. + ``sampling_interval`` : scalar, optional + Nominal sample spacing. When provided the coordinate is + *regular* and :meth:`get_sampling_interval` returns it directly. + ``tolerance`` : scalar, optional + Allowed jitter around ``sampling_interval``. Checked for + consistency with the tie points at construction. Ignored when + ``sampling_interval`` is absent. dim : str, optional Name of the dimension this coordinate is associated with. dtype : dtype-like, optional Desired dtype for ``tie_values``. + Notes + ----- + Regularity is judged on the continuous areas only: a tie-point gap of one + index (``den == 1``) is a CF discontinuity and carries no sampling-rate + information. A ``sampling_interval`` is valid when, for every continuous + segment, the accumulated drift ``|sampling_interval * den - num|`` stays + within ``2 * tolerance`` (each tie value may jitter by ±``tolerance``). + A coordinate with no continuous area (e.g. ``tie_indices=[0, 1, 2]``) has no + inferable spacing, so an explicitly provided one is stored as-is. Use + :meth:`simplify` to canonicalise a coordinate and acquire a spacing from + the continuous areas within an accuracy budget. + Examples -------- >>> import xdas as xd @@ -58,15 +91,17 @@ def __init__(self, data=None, dim=None, dtype=None): data = {"tie_indices": [], "tie_values": []} # parse data - data, dim = parse(data, dim) - if not self._isvalid(data): - raise TypeError("`data` must be dict-like") - if not set(data) == {"tie_indices", "tie_values"}: - raise ValueError( - "both `tie_indices` and `tie_values` key should be provided" + data, dim = parse_data_dim(data, dim) + if not InterpCoordinate._isvalid(data): + raise TypeError( + "`data` must be dict-like with `tie_indices` and `tie_values` " + "(and optionally `sampling_interval` / `tolerance`)" ) + tie_indices = np.asarray(data["tie_indices"]) tie_values = np.asarray(data["tie_values"], dtype=dtype) + sampling_interval = data.get("sampling_interval", None) + tolerance = data.get("tolerance", None) # check shapes if not tie_indices.ndim == 1: @@ -90,11 +125,14 @@ def __init__(self, data=None, dim=None, dtype=None): ): raise ValueError("`tie_values` must have either numeric or datetime dtype") - # store data + # store base data tie_indices = tie_indices.astype(int) self.data = dict(tie_indices=tie_indices, tie_values=tie_values) self.dim = dim + # optional regular sampling + self._assign_sampling_interval(sampling_interval, tolerance) + @property def tie_indices(self): """Integer array of tie-point positions (starts at 0, strictly increasing).""" @@ -105,6 +143,16 @@ def tie_values(self): """Array of tie-point values (numeric or datetime64, strictly increasing).""" return self.data["tie_values"] + @property + def sampling_interval(self): + """Nominal sample spacing, or ``None`` when the coordinate is not regular.""" + return self.data["sampling_interval"] + + @property + def tolerance(self): + """Allowed jitter around :attr:`sampling_interval`, or ``None``.""" + return self.data["tolerance"] + @property @override def dtype(self): @@ -113,16 +161,29 @@ def dtype(self): @classmethod @override def from_block(cls, start, size, step, dim=None, dtype=None): - data = { - "tie_indices": [0, size - 1], - "tie_values": [start, start + step * (size - 1)], - } + start = np.asarray(start, dtype=dtype) + step = parse_scalar_delta(step, start.dtype) + if size < 2: + # A single (or zero) sample cannot span two tie points; keep the + # declared spacing as metadata. + data = { + "tie_indices": [0][:size], + "tie_values": [start][:size], + "sampling_interval": step, + } + else: + end = start + step * (size - 1) + data = { + "tie_indices": [0, size - 1], + "tie_values": [start, end], + "sampling_interval": step, + } return cls(data, dim=dim, dtype=dtype) @override def __len__(self): if len(self.tie_indices) > 0: - return self.tie_indices[-1] - self.tie_indices[0] + 1 + return int(self.tie_indices[-1]) + 1 else: return 0 @@ -130,16 +191,17 @@ def __len__(self): @override def _isvalid(data): match data: - case {"tie_indices": _, "tie_values": _}: + case {"tie_indices": _, "tie_values": _, **rest} if set(rest) <= { + "sampling_interval", + "tolerance", + }: return True case _: return False @override def _is_monotonic_increasing(self): - return not self.get_split_indices( - "overlaps", tolerance=False - ).size # TODO: do not call split_indices + return not self.get_split_indices("overlaps", tolerance=False).size @override def _get_value(self, index): @@ -174,13 +236,9 @@ def _slice(self, index_slice): index_slice.step, ) if stop_index - start_index <= 0: - return self.__class__(dict(tie_indices=[], tie_values=[]), dim=self.dim) + data = {"tie_indices": [], "tie_values": []} elif (stop_index - start_index) <= step_index: - tie_indices = [0] - tie_values = [self._get_value(start_index)] - return self.__class__( - dict(tie_indices=tie_indices, tie_values=tie_values), dim=self.dim - ) + data = {"tie_indices": [0], "tie_values": [self._get_value(start_index)]} else: end_index = stop_index - 1 start_value = self._get_value(start_index) @@ -203,11 +261,17 @@ def _slice(self, index_slice): for k in range(1, len(tie_indices) - 1): if tie_indices[k] == tie_indices[k - 1]: tie_indices[k] += step_index - tie_values = [self._get_value(start_index + idx) for idx in tie_indices] + tie_values = self._get_value(start_index + tie_indices) tie_indices //= step_index data = {"tie_indices": tie_indices, "tie_values": tie_values} - return self.__class__(data, self.dim) + if self.sampling_interval is not None: + data = { + **data, + "sampling_interval": self.sampling_interval * step_index, + "tolerance": self.tolerance, + } + return self.__class__(data, self.dim) @override def _concat(self, other): @@ -221,16 +285,27 @@ def _concat(self, other): return self if not self.dtype == other.dtype: raise ValueError("cannot concatenate coordinate with different dtype") - coord = self.__class__( - { - "tie_indices": np.append( - self.tie_indices, other.tie_indices + len(self) - ), - "tie_values": np.append(self.tie_values, other.tie_values), - }, - self.dim, - ) - return coord + data = { + "tie_indices": np.append(self.tie_indices, other.tie_indices + len(self)), + "tie_values": np.append(self.tie_values, other.tie_values), + } + # Strict primitive: preserve the regular contract only when both sides + # advertise the exact same spacing; otherwise the merged coord is + # irregular by construction. The joining tie pair has ``den == 1`` (a + # CF discontinuity) so each side's segments validate independently, + # and ``max(tolerance)`` bounds the union. Reconciling slightly + # different rates is the job of user-facing routines (see + # :func:`concat_coords`, which delegates to :meth:`simplify`). + if ( + self.sampling_interval is not None + and self.sampling_interval == other.sampling_interval + ): + data = { + **data, + "sampling_interval": self.sampling_interval, + "tolerance": max(self.tolerance, other.tolerance), + } + return self.__class__(data, self.dim) @override def _to_dataset(self, dataset, attrs): @@ -249,6 +324,12 @@ def _to_dataset(self, dataset, attrs): "interpolation_name": "linear", "tie_points_mapping": f"{self.name}_points: {self.name}_indices {self.name}_values", } + if self.sampling_interval is not None: + interp_attrs.update( + encode_delta("sampling_interval", self.sampling_interval) + ) + if self.tolerance is not None: + interp_attrs.update(encode_delta("tolerance", self.tolerance)) dataset.update( { f"{self.name}_interpolation": ((), np.nan, interp_attrs), @@ -264,92 +345,294 @@ def _collect_from_dataset(cls, dataset, name): coords = {} mapping = dataset[name].attrs.pop("coordinate_interpolation", None) if mapping is not None: - matches = re.findall(r"(\w+): (\w+) (\w+)", mapping) - for match in matches: - dim, indices, values = match - data = {"tie_indices": dataset[indices], "tie_values": dataset[values]} + for dim, indices, values in re.findall(r"(\w+): (\w+) (\w+)", mapping): + data = { + "tie_indices": dataset[indices].values, + "tie_values": dataset[values].values, + } + interp_attrs = dataset[f"{dim}_interpolation"].attrs + if "sampling_interval" in interp_attrs: + data["sampling_interval"] = decode_delta( + "sampling_interval", interp_attrs + ) + data["tolerance"] = decode_delta("tolerance", interp_attrs) coords[dim] = Coordinate(data, dim) return coords def __add__(self, other): - return self.__class__( - {"tie_indices": self.tie_indices, "tie_values": self.tie_values + other}, - self.dim, - ) + data = {"tie_indices": self.tie_indices, "tie_values": self.tie_values + other} + if self.sampling_interval is not None: + data = { + **data, + "sampling_interval": self.sampling_interval, + "tolerance": self.tolerance, + } + return self.__class__(data, self.dim) def __sub__(self, other): - return self.__class__( - {"tie_indices": self.tie_indices, "tie_values": self.tie_values - other}, - self.dim, - ) + data = {"tie_indices": self.tie_indices, "tie_values": self.tie_values - other} + if self.sampling_interval is not None: + data = { + **data, + "sampling_interval": self.sampling_interval, + "tolerance": self.tolerance, + } + return self.__class__(data, self.dim) @override def get_sampling_interval(self, cast=True): - if len(self) < 2: + delta = self.sampling_interval + if delta is None: return None - num = np.diff(self.tie_values) - den = np.diff(self.tie_indices) - mask = den != 1 - num = num[mask] - den = den[mask] - if len(num) == 0: - return None - delta = np.median(num / den) if cast and np.issubdtype(delta.dtype, np.timedelta64): delta = delta / np.timedelta64(1, "s") return delta - @override - def simplify(self, tolerance=None): - if tolerance is False: - return self.copy() - tolerance = parse_tolerance(tolerance, self.dtype) - tie_indices, tie_values = _douglas_peucker( - self.tie_indices, self.tie_values, tolerance + def _assign_sampling_interval(self, sampling_interval, tolerance=None): + """Parse, validate and store the sampling interval and its tolerance. + + ``None`` clears both; a value is kept only if consistent with the tie + points (see :meth:`_is_valid_sampling_interval`). + """ + if sampling_interval is None: + if tolerance is not None: + raise ValueError( + "`tolerance` cannot be set without a `sampling_interval`" + ) + self.data["sampling_interval"] = None + self.data["tolerance"] = None + return + + sampling_interval = parse_scalar_delta(sampling_interval, self.dtype) + tolerance = parse_scalar_delta(tolerance, self.dtype, default_zero=True) + + # Normalise datetime deltas to the tie-value resolution so an in-memory + # coordinate matches the one read back after serialisation (which always + # encodes timedeltas at the coordinate's datetime resolution). + if np.issubdtype(self.dtype, np.datetime64): + unit = np.datetime_data(self.dtype)[0] + sampling_interval = sampling_interval.astype(f"timedelta64[{unit}]") + tolerance = tolerance.astype(f"timedelta64[{unit}]") + + if self._is_valid_sampling_interval(sampling_interval, tolerance): + self.data["sampling_interval"] = sampling_interval + self.data["tolerance"] = tolerance + else: + raise ValueError( + "`sampling_interval` and `tolerance` are not consistent with " + "the `tie_indices` and `tie_values`" + ) + + def _is_valid_sampling_interval(self, sampling_interval, tolerance): + """Whether *sampling_interval* fits every continuous area within *tolerance*.""" + num, den = self._continuous_segments() + # Bound the per-segment accumulated drift: each tie value may jitter by + # ±tolerance, so a segment span may be off by up to 2 * tolerance. With no + # continuous area `np.all([])` is vacuously True, accepting an explicit + # spacing as metadata (e.g. a two-tie-point block). Datetime bounds use + # integer division and are only accurate to the dtype resolution. + dmin = (num - 2 * tolerance) / den + dmax = (num + 2 * tolerance) / den + valid = np.all((dmin <= sampling_interval) & (sampling_interval <= dmax)) + return bool(valid) + + def _infer_regular(self): + """ + Estimate the nominal spacing and tightest tolerance for this coordinate. + + Private helper behind :meth:`simplify` and :meth:`to_regular`: returns + the spacing that minimises the worst per-segment drift and the smallest + tolerance that would still validate it, without enforcing either on the + coordinate. + + Returns + ------- + sampling_interval : scalar or None + Spacing minimising ``max_i |sampling_interval * den_i - num_i|`` + over the continuous segments. ``None`` when no continuous segment + is available (every tie-point gap is a ``den == 1`` CF + discontinuity). + tolerance : scalar or None + Half the worst residual drift at ``sampling_interval``, plus a few + ULPs so the value stays valid under re-validation. ``None`` when + ``sampling_interval`` is ``None``. + + Notes + ----- + Spacing is judged on continuous areas only, ``den == 1`` gaps being CF + discontinuities (see :meth:`_continuous_segments`). For such a segment + ``i``, ``num_i`` is the change in ``tie_values`` and ``den_i`` the + change in ``tie_indices``. The quantity ``si * den_i - num_i`` is the + drift accumulated between the regular grid and the tie values at the + end of that segment, and :meth:`_is_valid_sampling_interval` accepts + ``si`` exactly when every such drift stays within ``2 * tolerance``. + + The inferred spacing minimises the worst-case drift:: + + si* = argmin_si max_i |si * den_i - num_i| + + This convex, piecewise-linear objective is a length-weighted Chebyshev + center of the per-segment rates ``r = num / den``. Its minimum is + reached where the two most disagreeing segments balance, so over all + pairs the binding one maximises + ``den_i * den_j * |r_i - r_j| / (den_i + den_j)`` and the optimum is + the rate of that merged pair:: + + si* = (num_i + num_j) / (den_i + den_j) + + That pair is found in ``O(n log n)`` via :func:`_chebyshev_center_pair` + rather than scanning all pairs. The matching tolerance is half the + worst drift, since validity compares the drift against + ``2 * tolerance``. + """ + num, den = self._continuous_segments() + if num.size == 0: + return None, None + # Float seconds for datetime axes pick the binding pair without + # integer/timedelta overflow; the final values stay in the native dtype. + is_datetime = np.issubdtype(num.dtype, np.timedelta64) + num_seconds = num / np.timedelta64(1, "s") if is_datetime else num.astype(float) + pos_idx, neg_idx = _chebyshev_center_pair(num_seconds, den.astype(float)) + sampling_interval = (num[pos_idx] + num[neg_idx]) / ( + den[pos_idx] + den[neg_idx] ) - return self.__class__( - dict(tie_indices=tie_indices, tie_values=tie_values), self.dim + si_seconds = ( + sampling_interval / np.timedelta64(1, "s") + if is_datetime + else float(sampling_interval) ) + drift = np.abs(si_seconds * den - num_seconds).max() + # A few ULPs of slack so re-validation cannot reject the returned pair. + tolerance = drift / 2 + 4 * np.spacing(np.abs(num_seconds).max()) + if is_datetime: + tolerance = np.timedelta64(int(np.ceil(tolerance * 1e9)), "ns") + return sampling_interval, tolerance @override - def get_split_indices(self, kind="discontinuities", tolerance=False): - valid_kinds = {"discontinuities", "gaps", "overlaps"} - if kind not in valid_kinds: - raise ValueError(f"`kind` must be one of {valid_kinds}; got {kind!r}") - - (indices,) = np.nonzero(np.diff(self.tie_indices) == 1) - indices += 1 - - # Fast path: no filtering requested - if kind == "discontinuities" and tolerance is False: - return self.tie_indices[indices] - - sampling_interval = self.get_sampling_interval(cast=False) - deltas = ( - self.tie_values[indices] - self.tie_values[indices - 1] - sampling_interval - ) + def to_regular(self, sampling_interval=None, tolerance=None): + """Enforce a nominal sampling interval, inferring it when omitted. + + The inferred spacing is the length-weighted Chebyshev center of the + per-segment rates (see :meth:`_infer_regular`). Raises when no spacing + can be inferred (no continuous area, i.e. every tie-point gap is a + ``den == 1`` CF discontinuity) or when the spacing does not fit the tie + points within *tolerance*. See :meth:`AxisCoordinate.to_regular` for + the parameter contract. + """ + # Default each unspecified argument to the stored regular config; an + # explicit value still overrides it. + if sampling_interval is None: + sampling_interval = self.sampling_interval + if tolerance is None: + tolerance = self.tolerance + if sampling_interval is None: + sampling_interval, _ = self._infer_regular() + if sampling_interval is None: + raise ValueError( + "cannot infer a sampling interval: the coordinate has no " + "continuous area; pass `sampling_interval` explicitly" + ) + data = { + "tie_indices": self.tie_indices, + "tie_values": self.tie_values, + "sampling_interval": sampling_interval, + "tolerance": tolerance, + } + return self.__class__(data, self.dim) + @override + def simplify(self, tolerance=None, *, reduce=True, regularize=False): + """Canonicalise within *tolerance*: drop tie points, then promote to regular. + + The *reduce* stage runs Douglas-Peucker to drop tie points whose removal + shifts the curve by no more than *tolerance*. The CF 8.3 structure is + preserved as an emergent property of that bound: real discontinuities are + kept (any spanning line crosses them by far more than *tolerance*), soft + ones are fused into a single ramp, and synchronisation tie points survive + because removing them would, by definition, drift more than *tolerance*. + Surviving values are never moved. + + The *regularize* stage promotes the result to *regular* when the + surviving continuous segments admit a single ``sampling_interval`` within + *tolerance* (the internal Chebyshev fit's worst residual stays inside the + budget). The promotion is per-continuous-segment and sign-agnostic, so + two same-rate segments joined by a CF overlap are still described by one + spacing. An already-regular coordinate keeps its spacing regardless of + *regularize* and just widens its stored tolerance to absorb any jump + fused by the reduce stage. + + See :meth:`Coordinate.simplify` for the parameter contract. + """ if tolerance is False: - zero = np.timedelta64(0) if np.issubdtype(self.dtype, np.datetime64) else 0 - - match kind: - case "gaps": - mask = deltas >= zero - case "overlaps": # pragma: no branch - mask = deltas < zero - + return self.copy() + if tolerance is None: + # Default the budget to the coordinate's own declared jitter. + tolerance = self.tolerance + tolerance = parse_scalar_delta(tolerance, self.dtype, default_zero=True) + if reduce: + tie_indices, tie_values = _douglas_peucker( + self.tie_indices, self.tie_values, tolerance + ) else: - tolerance = parse_tolerance(tolerance, self.dtype) + tie_indices, tie_values = self.tie_indices, self.tie_values + data = {"tie_indices": tie_indices, "tie_values": tie_values} + if self.sampling_interval is not None: + # Already regular: keep the spacing. A reduce pass may fuse a soft + # discontinuity into a ramp whose absorbed jump exceeds the declared + # jitter; keep the declared tolerance when it still validates and + # only widen by the spent budget when it does not, so a lossless or + # seam-only reduction round-trips the exact regular metadata. + reduced = self.__class__(data, self.dim) + if reduce and not reduced._is_valid_sampling_interval( + self.sampling_interval, self.tolerance + ): + new_tolerance = self.tolerance + tolerance + else: + new_tolerance = self.tolerance + data = { + **data, + "sampling_interval": self.sampling_interval, + "tolerance": new_tolerance, + } + return self.__class__(data, self.dim) + # Otherwise try to promote: infer the best spacing on the surviving + # continuous segments and keep it only if it validates within the budget. + if regularize: + reduced = self.__class__(data, self.dim) + sampling_interval, _ = reduced._infer_regular() + if sampling_interval is not None and reduced._is_valid_sampling_interval( + sampling_interval, tolerance + ): + data = { + **data, + "sampling_interval": sampling_interval, + "tolerance": tolerance, + } + return self.__class__(data, self.dim) - match kind: - case "discontinuities": - mask = np.abs(deltas) > tolerance - case "gaps": - mask = deltas > tolerance - case "overlaps": # pragma: no branch - mask = deltas < -tolerance + @override + def _split_candidates(self): + """Discontinuity split points, each paired with its step's deviation from the neighbouring interval.""" + tie_intervals = np.diff(self.tie_values) / np.diff(self.tie_indices) + (positions,) = np.nonzero(np.diff(self.tie_indices) == 1) + references = np.where( + positions > 0, + positions - 1, + np.minimum(positions + 1, len(tie_intervals) - 1), + ) + deltas = tie_intervals[positions] - tie_intervals[references] + return self.tie_indices[positions + 1], deltas - return self.tie_indices[indices[mask]] + def _continuous_segments(self): + """Per-segment value/index spans ``(num, den)`` for the continuous areas. + + A ``den == 1`` gap is a CF discontinuity (section 8.3), not a segment, so + it is excluded and carries no sampling-rate information. + """ + num = np.diff(self.tie_values) + den = np.diff(self.tie_indices) + mask = den != 1 + return num[mask], den[mask] def _douglas_peucker(x, y, epsilon): @@ -392,3 +675,72 @@ def _douglas_peucker(x, y, epsilon): else: mask[start + 1 : stop - 1] = False return x[mask], y[mask] + + +def _chebyshev_center_pair(num, den): + """ + Segment indices binding the length-weighted Chebyshev center, in O(n log n). + + Returns the pair maximising ``den_i den_j |r_i - r_j| / (den_i + den_j)`` with + ``r = num / den``, equivalently the lowest point of the upper envelope of the + ``2 n`` lines ``±(den_i si - num_i)``. That vertex is the meeting of the + binding negative- and positive-slope lines, found with the convex-hull trick + instead of the O(n^2) pairwise scan. + + Parameters + ---------- + num : numpy.ndarray + Per-segment numerators as floats (seconds for datetime axes). + den : numpy.ndarray + Per-segment denominators as floats, all strictly positive. + + Returns + ------- + pos_idx, neg_idx : int + Segment indices of the binding positive- and negative-slope lines. A + single segment trivially selects itself (``pos_idx == neg_idx``). + """ + seg = np.arange(len(den)) + # Positive-slope lines (den si - num) and negative-slope lines (num - den si). + slopes = np.concatenate([den, -den]) + intercepts = np.concatenate([-num, num]) + idx = np.concatenate([seg, seg]) + # Process lines by ascending slope, equal slopes ordered by descending + # intercept so the dominant one comes first. + order = np.lexsort((-intercepts, slopes)) + return _upper_envelope_min_pair(slopes, intercepts, idx, order) + + +@njit(cache=True) +def _upper_envelope_min_pair(slopes, intercepts, idx, order): # pragma: no cover + """Binding (positive, negative) line indices at the upper-envelope minimum.""" + n = order.size + hull_s = np.empty(n, dtype=slopes.dtype) + hull_b = np.empty(n, dtype=intercepts.dtype) + hull_i = np.empty(n, dtype=idx.dtype) + m = 0 # current hull size + for t in range(n): + k = order[t] + s, b, seg = slopes[k], intercepts[k], idx[k] + # Equal slopes: the dominant (larger intercept) one came first; skip rest. + if m > 0 and hull_s[m - 1] == s: + continue + # Drop any line the convex-hull trick proves can never be the maximum. + while m >= 2: + s1, b1 = hull_s[m - 2], hull_b[m - 2] + s2, b2 = hull_s[m - 1], hull_b[m - 1] + if (b - b1) / (s1 - s) <= (b2 - b1) / (s1 - s2): + m -= 1 + else: + break + hull_s[m], hull_b[m], hull_i[m] = s, b, seg + m += 1 + # The envelope is convex with slope increasing along x; its minimum sits at + # the negative-to-positive slope transition, between the binding pair. The + # scan is safely bounded because `_chebyshev_center_pair` always feeds in + # both `+den_i` and `-den_i` lines, so at least one slope of each sign + # reaches the hull. + t = 0 + while hull_s[t] < 0.0: + t += 1 + return hull_i[t], hull_i[t - 1] diff --git a/xdas/coordinates/sampled.py b/xdas/coordinates/sampled.py index 69f926d4..967fb136 100644 --- a/xdas/coordinates/sampled.py +++ b/xdas/coordinates/sampled.py @@ -10,25 +10,17 @@ from typing_extensions import override from .core import ( + CODE_TO_UNITS, + UNITS_TO_CODE, + AxisCoordinate, Coordinate, - SampledMixin, is_monotonic_increasing, - parse, - parse_tolerance, + parse_data_dim, + parse_scalar_delta, ) -CODE_TO_UNITS = { - "h": "hours", - "m": "minutes", - "s": "seconds", - "ms": "milliseconds", - "us": "microseconds", - "ns": "nanoseconds", -} -UNITS_TO_CODE = {v: k for k, v in CODE_TO_UNITS.items()} - -class SampledCoordinate(SampledMixin, Coordinate, ctype="sampled"): +class SampledCoordinate(AxisCoordinate, ctype="sampled"): """ Coordinate sampled at a fixed interval, with optional gaps between segments. @@ -78,8 +70,8 @@ def __init__(self, data=None, dim=None, dtype=None): empty = False # parse data - data, dim = parse(data, dim) - if not self.__class__._isvalid(data): + data, dim = parse_data_dim(data, dim) + if not self._isvalid(data): raise ValueError( "`data` must be dict-like and contain `tie_values`, `tie_lengths`, and " "`sampling_interval`" @@ -425,10 +417,32 @@ def get_sampling_interval(self, cast=True): return delta @override - def simplify(self, tolerance=None): - if tolerance is False: + def to_regular(self, sampling_interval=None, tolerance=None): + """Regular by construction: validate any explicit spacing and return a copy. + + See :meth:`AxisCoordinate.to_regular` for the parameter contract. + """ + if sampling_interval is not None: + sampling_interval = parse_scalar_delta(sampling_interval, self.dtype) + tolerance = parse_scalar_delta(tolerance, self.dtype, default_zero=True) + if np.abs(sampling_interval - self.sampling_interval) > tolerance: + raise ValueError( + "`sampling_interval` does not match the stored sampling interval" + ) + return self.copy() + + @override + def simplify(self, tolerance=None, *, reduce=True, regularize=False): + """Fuse adjacent segments whose junction drift is within *tolerance*. + + The coordinate is regular by construction (it carries a single + ``sampling_interval``), so *regularize* is a no-op; fusing happens only + when *reduce* is set. See :meth:`Coordinate.simplify` for the parameter + contract. + """ + if tolerance is False or not reduce: return self.copy() - tolerance = parse_tolerance(tolerance, self.dtype) + tolerance = parse_scalar_delta(tolerance, self.dtype, default_zero=True) tie_values = [self.tie_values[0]] tie_lengths = [self.tie_lengths[0]] for value, length in zip(self.tie_values[1:], self.tie_lengths[1:]): @@ -448,39 +462,8 @@ def simplify(self, tolerance=None): ) @override - def get_split_indices(self, kind="discontinuities", tolerance=False): - valid_kinds = {"discontinuities", "gaps", "overlaps"} - if kind not in valid_kinds: - raise ValueError(f"`kind` must be one of {valid_kinds}; got {kind!r}") - - indices = self.tie_indices[1:] - - # Fast path: no filtering requested - if kind == "discontinuities" and tolerance is False: - return indices - + def _split_candidates(self): deltas = self.tie_values[1:] - ( self.tie_values[:-1] + self.sampling_interval * self.tie_lengths[:-1] ) - - if tolerance is False: - zero = np.timedelta64(0) if np.issubdtype(self.dtype, np.datetime64) else 0 - - match kind: # pragma: no branch - case "gaps": - mask = deltas >= zero - case "overlaps": # pragma: no branch - mask = deltas < zero - - else: - tolerance = parse_tolerance(tolerance, self.dtype) - - match kind: # pragma: no branch - case "discontinuities": - mask = np.abs(deltas) > tolerance - case "gaps": - mask = deltas > tolerance - case "overlaps": # pragma: no branch - mask = deltas < -tolerance - - return indices[mask] + return self.tie_indices[1:], deltas diff --git a/xdas/coordinates/scalar.py b/xdas/coordinates/scalar.py index 914c02e9..a0181211 100644 --- a/xdas/coordinates/scalar.py +++ b/xdas/coordinates/scalar.py @@ -7,16 +7,18 @@ import numpy as np from typing_extensions import override -from .core import Coordinate, parse +from .core import Coordinate, parse_data_dim class ScalarCoordinate(Coordinate, ctype="scalar"): """ Non-dimensional coordinate that carries a single scalar value. - Unlike dimensional coordinates, a :class:`ScalarCoordinate` is not tied - to an array axis and has no length. Typical use: metadata attached to a - :class:`DataArray` (e.g. an instrument identifier or a shot time). + Unlike :class:`~xdas.coordinates.AxisCoordinate` subclasses, a + :class:`ScalarCoordinate` is not tied to an array axis and has no length. + It therefore implements only the thin :class:`Coordinate` interface. + Typical use: metadata attached to a :class:`DataArray` (e.g. an instrument + identifier or a shot time). Parameters ---------- @@ -32,38 +34,18 @@ class ScalarCoordinate(Coordinate, ctype="scalar"): def __init__(self, data=None, dim=None, dtype=None): if data is None: raise TypeError("scalar coordinate cannot be empty, please provide a value") - data, dim = parse(data, dim) + data, dim = parse_data_dim(data, dim) if dim is not None: raise ValueError("a scalar coordinate cannot be a dim") if not self._isvalid(data): raise TypeError("`data` must be scalar-like") self.data = np.asarray(data, dtype=dtype) - @classmethod - @override - def from_block(cls, start, size, step, dim=None, dtype=None): - raise TypeError("cannot build a scalar coordinate from a block") - - @override - def __len__(self): - return 1 - - @override - def __getitem__(self, item): - raise TypeError("scalar coordinate is not subscriptable") - - @override - def __array__(self, dtype=None, copy=None): - # TODO: drop this workaround once Python 3.10 is no longer supported - # (EOL Oct 2026). numpy < 2.3 raises when copy=False on a 0-d array; - # numpy 2.3+ (requires Python 3.11+) handles it correctly. - if copy: - return np.array(self.data, dtype=dtype) - return np.asarray(self.data, dtype=dtype) - + @staticmethod @override - def __repr__(self): - return np.array2string(self.data, threshold=0, edgeitems=1) + def _isvalid(data): + data = np.asarray(data) + return (data.dtype != np.dtype(object)) and (data.ndim == 0) @property def dim(self): @@ -82,8 +64,8 @@ def dtype(self): return self.data.dtype @property - @override def ndim(self): + """Always ``0`` — scalar coordinates have no axis.""" return 0 @property @@ -91,46 +73,18 @@ def ndim(self): def shape(self): return () - @property - @override - def indices(self): - raise TypeError("scalar coordinate has no indices") - - @property - @override - def start(self): - raise TypeError("scalar coordinate has no start") - - @property - @override - def end(self): - raise TypeError("scalar coordinate has no end") - - @staticmethod - @override - def _isvalid(data): - data = np.asarray(data) - return (data.dtype != np.dtype(object)) and (data.ndim == 0) - @override - def _is_monotonic_increasing(self): - raise TypeError("scalar coordinate has no axis") - - @override - def _get_value(self, index): - raise TypeError("scalar coordinate has no elements to index") - - @override - def _get_indexer(self, value, method=None): - raise TypeError("cannot get index of scalar coordinate") - - @override - def _slice(self, slc): - raise TypeError("scalar coordinate is not sliceable") + def __array__(self, dtype=None, copy=None): + # TODO: drop this workaround once Python 3.10 is no longer supported + # (EOL Oct 2026). numpy < 2.3 raises when copy=False on a 0-d array; + # numpy 2.3+ (requires Python 3.11+) handles it correctly. + if copy: + return np.array(self.data, dtype=dtype) + return np.asarray(self.data, dtype=dtype) @override - def _concat(self, other): - raise TypeError("cannot concatenate scalar coordinate") + def __repr__(self): + return np.array2string(self.data, threshold=0, edgeitems=1) @override def _to_dataset(self, dataset, attrs): @@ -145,15 +99,3 @@ def _to_dataset(self, dataset, attrs): @override def _collect_from_dataset(cls, dataset, name): return {} - - @override - def get_sampling_interval(self, cast=True): - return None - - @override - def isscalar(self): - return True - - @override - def to_index(self, item, method=None, endpoint=True): - raise NotImplementedError("cannot get index of scalar coordinate") diff --git a/xdas/core/__init__.py b/xdas/core/__init__.py index c0c68b57..cb570eb0 100644 --- a/xdas/core/__init__.py +++ b/xdas/core/__init__.py @@ -4,3 +4,49 @@ Includes :class:`DataArray`, :class:`DataCollection`, and supporting routines, methods, and NumPy dispatch. """ + +__all__ = [ + "DataArray", + "DataCollection", + "DataMapping", + "DataSequence", + "align", + "asdataarray", + "broadcast_coords", + "broadcast_to", + "combine_by_coords", + "combine_by_field", + "concat", + "concat_coords", + "concatenate", + "open", + "open_dataarray", + "open_datacollection", + "open_mfdataarray", + "open_mfdatacollection", + "open_mfdatatree", + "plot_availability", + "split", +] + +from .dataarray import DataArray +from .datacollection import DataCollection, DataMapping, DataSequence +from .routines import ( + align, + asdataarray, + broadcast_coords, + broadcast_to, + combine_by_coords, + combine_by_field, + concat, + concat_coords, + concatenate, + open, + open_dataarray, + open_datacollection, + open_mfdataarray, + open_mfdatacollection, + open_mfdatatree, + plot_availability, + split, +) diff --git a/xdas/core/dataarray.py b/xdas/core/dataarray.py index 4fb7de3d..c9b1d60b 100644 --- a/xdas/core/dataarray.py +++ b/xdas/core/dataarray.py @@ -14,7 +14,7 @@ from dask.array import Array as DaskArray from numpy.lib.mixins import NDArrayOperatorsMixin -from ..coordinates import Coordinates +from ..coordinates import AxisCoordinate, Coordinates from ..virtual import _to_human HANDLED_NUMPY_FUNCTIONS = {} @@ -347,7 +347,7 @@ def isel(self, indexers=None, drop=False, **indexers_kwargs): da = self[indexers] if drop: for dim in indexers: - if da[dim].isscalar(): + if not isinstance(da[dim], AxisCoordinate): da = da.drop_coords(dim) return da @@ -394,7 +394,7 @@ def sel( f"dimension {dim} is not monotonic increasing, " f"spliting on overlaps, slicing and concatenating can be slow..." ) - from ..core.routines import concat, split + from .routines import concat, split chunks = [ chunk.sel(indexers, method, endpoint, drop) @@ -411,7 +411,7 @@ def sel( da = self[key] if drop: for dim in indexers: - if da[dim].isscalar(): + if not isinstance(da[dim], AxisCoordinate): da = da.drop_coords(dim) return da @@ -777,7 +777,7 @@ def expand_dims(self, dim, axis=0): raise ValueError(f"cannot expand on existing dimension {dim}") coords = self.coords.copy() if dim in coords: - if coords[dim].isscalar(): + if not isinstance(coords[dim], AxisCoordinate): coords[dim] = [coords[dim].values] else: raise ValueError( diff --git a/xdas/core/methods.py b/xdas/core/methods.py index 6d43ec5f..f0e2ff91 100644 --- a/xdas/core/methods.py +++ b/xdas/core/methods.py @@ -6,7 +6,7 @@ import numpy as np -from ..atoms.core import atomized +from ..atoms import atomized from .dataarray import HANDLED_METHODS diff --git a/xdas/core/routines.py b/xdas/core/routines.py index 28157ed7..b2806656 100644 --- a/xdas/core/routines.py +++ b/xdas/core/routines.py @@ -21,7 +21,7 @@ from loky import get_reusable_executor from tqdm import tqdm -from ..coordinates.core import Coordinates, SampledMixin, get_sampling_interval +from ..coordinates import AxisCoordinate, Coordinates from ..parallel import get_workers_count from ..virtual import VirtualSource, VirtualStack from .dataarray import DataArray @@ -803,7 +803,11 @@ def combine_by_coords( if dim in objs[0].coords: objs = sorted( objs, - key=lambda da: da[dim].values if da[dim].isscalar() else da[dim][0].values, + key=lambda da: ( + da[dim][0].values + if isinstance(da[dim], AxisCoordinate) + else da[dim].values + ), ) # combine objs @@ -865,12 +869,18 @@ def initialize(self, da): if self.dim in self.dims else da.coords.drop_coords(self.dim) ) - if self.dim in da.coords: - self.delta = get_sampling_interval(da, self.dim) - else: - self.delta = None + self.delta = self._get_delta(da) self.dtype = da.dtype + def _get_delta(self, da): + """Nominal sampling interval of *da* along *dim*, or ``None`` (irregular or absent).""" + if self.dim not in da.coords: + return None + coord = da.coords[self.dim] + if not isinstance(coord, AxisCoordinate): + return None + return coord.get_sampling_interval() + def append(self, da): """Add *da* after running all compatibility checks; initialises on first call.""" if not self.objs: @@ -914,12 +924,21 @@ def check_sampling_interval(self, da): if self.delta is None: pass else: - delta = get_sampling_interval(da, self.dim) - if not np.isclose(delta, self.delta): + delta = self._get_delta(da) + if delta is None or not np.isclose(delta, self.delta): raise CompatibilityError("sampling intervals are not compatible") -def concat(objs, dim="first", tolerance=None, virtual=None, verbose=None): +def concat( + objs, + dim="first", + tolerance=None, + virtual=None, + verbose=None, + *, + reduce=True, + regularize=False, +): """ Concatenate data arrays along a given dimension. @@ -932,12 +951,22 @@ def concat(objs, dim="first", tolerance=None, virtual=None, verbose=None): tolerance : float or timedelta64, optional The tolerance to consider that the end of a file is continuous with beginning of the following, For time coordinates, numeric values are considered as seconds. - Zero by default. + By default each coordinate spends its own declared tolerance when it + carries one, else a zero-like default. Pass ``False`` to disable + simplification entirely. virtual : bool, optional Whether to create a virtual dataset. It requires that all concatenated data arrays are virtual. By default tries to create a virtual dataset if possible. verbose: bool Whether to display a progress bar. + reduce : bool, optional + Whether to drop redundant tie points from the concatenated coordinate. + Default True. + regularize : bool, optional + Whether to promote the concatenated coordinate to a regular one when its + segments admit a single shared rate within *tolerance*. Default False: + regular inputs already stay regular through concatenation, so promotion + only matters for irregular inputs and stays opt-in. Returns ------- @@ -976,6 +1005,8 @@ def concat(objs, dim="first", tolerance=None, virtual=None, verbose=None): sort=True, return_order=True, tolerance=tolerance, + reduce=reduce, + regularize=regularize, ) objs = [objs[idx] for idx in order] coords[dim] = coord @@ -1002,7 +1033,15 @@ def concat(objs, dim="first", tolerance=None, virtual=None, verbose=None): concatenate = concat # TODO: deprecate it -def concat_coords(objs, *, sort=False, return_order=False, tolerance=False): +def concat_coords( + objs, + *, + sort=False, + return_order=False, + tolerance=None, + reduce=True, + regularize=False, +): """ Concatenate coordinate objects. @@ -1018,7 +1057,15 @@ def concat_coords(objs, *, sort=False, return_order=False, tolerance=False): tolerance : float or timedelta64, optional The tolerance to consider that the end of a coordinate object is continuous with beginning of the following, For time coordinates, numeric values are - considered as seconds. No simplification by default. + considered as seconds. By default the coordinate spends its own declared + tolerance when it carries one, else a zero-like default. Pass ``False`` + to disable simplification entirely. + reduce : bool, optional + Whether to drop redundant tie points after concatenation. Default True. + regularize : bool, optional + Whether to promote the result to a regular coordinate when the merged + segments admit a single shared rate within *tolerance*. Default False: + regular inputs already stay regular through concatenation. Returns ------- @@ -1041,11 +1088,14 @@ def concat_coords(objs, *, sort=False, return_order=False, tolerance=False): # simplify if tolerance is not False: - if isinstance(out, SampledMixin): - out = out.simplify(tolerance) - elif ( - tolerance is not None - ): # TODO: Default to False and remove this condition here? + if isinstance(out, AxisCoordinate): + # `_concat` is strict: same-rate inputs stay regular, mismatched + # rates drop to irregular. `simplify` then drops redundant tie + # points (chunk seams within tolerance fuse away) and, with + # `regularize=True`, recovers a single shared rate when the merged + # segments admit one within *tolerance*. + out = out.simplify(tolerance, reduce=reduce, regularize=regularize) + elif tolerance is not None: raise TypeError( "`tolerance` can only be used with coordinates " "that implements `simplify`" @@ -1216,7 +1266,7 @@ def broadcast_coords(*objs): else: sizes[dim] = size for name, coord in obj.coords.items(): - if coord.isscalar(): + if not isinstance(coord, AxisCoordinate): continue if name in coords: if not coord.equals(coords[name]): diff --git a/xdas/fft.py b/xdas/fft.py index 3c498268..bcd8f4c7 100644 --- a/xdas/fft.py +++ b/xdas/fft.py @@ -7,9 +7,9 @@ import numpy as np -from .atoms.core import atomized -from .coordinates.core import get_sampling_interval -from .core.dataarray import DataArray +from .atoms import atomized +from .coordinates import get_sampling_interval +from .core import DataArray from .parallel import parallelize @@ -52,12 +52,13 @@ def fft(da, n=None, dim={"last": "spectrum"}, norm=None, parallel=None): -------- >>> import xdas as xd >>> import xdas.fft as xfft - >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0, 1, 2, 3]}) + >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0., 1., 2., 3.]}) + >>> signal["time"] = signal["time"].to_regular() >>> xfft.fft(signal, dim={"time": "frequency"}) [0.+0.j 0.+2.j 0.+0.j 0.-2.j] Coordinates: - * frequency (frequency): [-0.5 ... 0.25] + * frequency (frequency): -0.500 to 0.250 """ ((olddim, newdim),) = dim.items() @@ -66,7 +67,8 @@ def fft(da, n=None, dim={"last": "spectrum"}, norm=None, parallel=None): n = da.sizes[olddim] axis = da.get_axis_num(olddim) d = get_sampling_interval(da, olddim) - f = np.fft.fftshift(np.fft.fftfreq(n, d)) + start = np.fft.fftshift(np.fft.fftfreq(n, d))[0] + f = type(da.coords[olddim]).from_block(start, n, 1 / (n * d), dim=newdim) def func(x): return np.fft.fftshift(np.fft.fft(x, n, axis, norm), axis) @@ -123,12 +125,13 @@ def rfft(da, n=None, dim={"last": "spectrum"}, norm=None, parallel=None): -------- >>> import xdas as xd >>> import xdas.fft as xfft - >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0, 1, 2, 3]}) + >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0., 1., 2., 3.]}) + >>> signal["time"] = signal["time"].to_regular() >>> xfft.rfft(signal, dim={"time": "frequency"}) [0.+0.j 0.-2.j 0.+0.j] Coordinates: - * frequency (frequency): [0. ... 0.5] + * frequency (frequency): 0.000 to 0.500 """ ((olddim, newdim),) = dim.items() @@ -139,7 +142,7 @@ def rfft(da, n=None, dim={"last": "spectrum"}, norm=None, parallel=None): d = get_sampling_interval(da, olddim) across = int(axis == 0) func = parallelize(across, across, parallel)(np.fft.rfft) - f = np.fft.rfftfreq(n, d) + f = type(da.coords[olddim]).from_block(0.0, n // 2 + 1, 1 / (n * d), dim=newdim) data = func(da.values, n, axis, norm) coords = { newdim if name == olddim else name: f if name == olddim else da.coords[name] @@ -186,7 +189,8 @@ def ifft(da, n=None, dim={"last": "signal"}, norm=None, parallel=None): -------- >>> import xdas as xd >>> import xdas.fft as xfft - >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0, 1, 2, 3]}) + >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0., 1., 2., 3.]}) + >>> signal["time"] = signal["time"].to_regular() >>> spectrum = xfft.fft(signal, dim={"time": "frequency"}) >>> result = xfft.ifft(spectrum, dim={"frequency": "time"}) >>> result["time"] = signal["time"] # to match time coordinates @@ -199,7 +203,8 @@ def ifft(da, n=None, dim={"last": "signal"}, norm=None, parallel=None): n = da.sizes[olddim] axis = da.get_axis_num(olddim) d = get_sampling_interval(da, olddim) - f = np.fft.ifftshift(np.fft.fftfreq(n, d)) + start = np.fft.fftshift(np.fft.fftfreq(n, d))[0] + f = type(da.coords[olddim]).from_block(start, n, 1 / (n * d), dim=newdim) def func(x): return np.fft.ifft(np.fft.ifftshift(x, axis), n, axis, norm) @@ -257,7 +262,8 @@ def irfft(da, n=None, dim={"last": "signal"}, norm=None, parallel=None): -------- >>> import xdas as xd >>> import xdas.fft as xfft - >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0, 1, 2, 3]}) + >>> signal = xd.DataArray([0., 1., 0., -1.], coords={"time": [0., 1., 2., 3.]}) + >>> signal["time"] = signal["time"].to_regular() >>> spectrum = xfft.rfft(signal, dim={"time": "frequency"}) >>> result = xfft.irfft( ... spectrum, @@ -276,7 +282,8 @@ def irfft(da, n=None, dim={"last": "signal"}, norm=None, parallel=None): d = get_sampling_interval(da, olddim) across = int(axis == 0) func = parallelize(across, across, parallel)(np.fft.irfft) - f = np.fft.fftshift(np.fft.fftfreq(n, d)) + start = np.fft.fftshift(np.fft.fftfreq(n, d))[0] + f = type(da.coords[olddim]).from_block(start, n, 1 / (n * d), dim=newdim) data = func(da.values, n, axis, norm) coords = { newdim if name == olddim else name: f if name == olddim else da.coords[name] diff --git a/xdas/io/apsensing.py b/xdas/io/apsensing.py index 54630129..df7a092a 100644 --- a/xdas/io/apsensing.py +++ b/xdas/io/apsensing.py @@ -3,8 +3,8 @@ import h5py import numpy as np -from ..coordinates.core import Coordinate -from ..core.dataarray import DataArray +from ..coordinates import Coordinate +from ..core import DataArray from ..virtual import VirtualSource from .core import Engine diff --git a/xdas/io/asn.py b/xdas/io/asn.py index 49b7ea5e..81f25be1 100644 --- a/xdas/io/asn.py +++ b/xdas/io/asn.py @@ -12,8 +12,8 @@ import numpy as np import zmq -from ..coordinates.core import Coordinate, get_sampling_interval -from ..core.dataarray import DataArray +from ..coordinates import Coordinate, get_sampling_interval +from ..core import DataArray from ..virtual import VirtualSource from .core import Engine @@ -122,7 +122,7 @@ def __init__(self, address): >>> address = f"tcp://localhost:{port}" >>> publisher = ZMQPublisher(address) - >>> da = xd.synthetics.dummy() + >>> da = xd.testing.dummy() >>> chunks = xd.split(da, 10) >>> def publish(): @@ -174,18 +174,20 @@ def _update_header(self, message): roiTable = header["roiTable"][0] di = (roiTable["roiStart"] // roiTable["roiDec"]) * header["dx"] de = (roiTable["roiEnd"] // roiTable["roiDec"]) * header["dx"] - self.distance = { # TODO: use from_block + self.distance = { "tie_indices": [0, header["nChannels"] - 1], "tie_values": [di, de], + "sampling_interval": (de - di) / (header["nChannels"] - 1), } self.delta = float_to_timedelta(header["dt"], header["dtUnit"]) def _unpack(self, message): t0 = np.frombuffer(message[:8], "datetime64[ns]").reshape(()) data = np.frombuffer(message[8:], self.dtype).reshape(self.shape) - time = { # TODO: use from_block + time = { "tie_indices": [0, self.shape[0] - 1], "tie_values": [t0, t0 + (self.shape[0] - 1) * self.delta], + "sampling_interval": self.delta, } return DataArray(data, {"time": time, "distance": self.distance}) @@ -214,7 +216,7 @@ class ZMQPublisher: >>> import xdas as xd >>> from xdas.io.asn import ZMQPublisher - >>> da = xd.synthetics.dummy() + >>> da = xd.testing.dummy() >>> port = xd.io.get_free_port() >>> address = f"tcp://localhost:{port}" diff --git a/xdas/io/febus.py b/xdas/io/febus.py index 6915a5f4..588a7773 100644 --- a/xdas/io/febus.py +++ b/xdas/io/febus.py @@ -5,9 +5,8 @@ import h5py import numpy as np -from ..coordinates.core import Coordinate -from ..core.dataarray import DataArray -from ..core.routines import concat +from ..coordinates import Coordinate +from ..core import DataArray, concat from ..virtual import VirtualSource from .core import Engine diff --git a/xdas/io/miniseed.py b/xdas/io/miniseed.py index 4af52943..ebdb203f 100644 --- a/xdas/io/miniseed.py +++ b/xdas/io/miniseed.py @@ -4,9 +4,13 @@ import numpy as np import obspy -from ..coordinates.core import Coordinate, Coordinates, get_sampling_interval -from ..core.dataarray import DataArray -from ..core.routines import concat_coords +from ..coordinates import ( + AxisCoordinate, + Coordinate, + Coordinates, + get_sampling_interval, +) +from ..core import DataArray, concat_coords from .core import Engine @@ -81,7 +85,9 @@ def read_header(self, path, ignore_last_sample, ctype): } ) - shape = tuple(len(coord) for coord in coords.values() if not coord.isscalar()) + shape = tuple( + len(coord) for coord in coords.values() if isinstance(coord, AxisCoordinate) + ) return shape, dtype, coords, method def read_data(self, path, method, ignore_last_sample): @@ -174,13 +180,11 @@ def from_stream(st, dims=("channel", "time")): """ data = np.stack([tr.data for tr in st]) channel = [tr.id for tr in st] - time = { - "tie_indices": [0, st[0].stats.npts - 1], - "tie_values": [ - np.datetime64(st[0].stats.starttime.datetime), - np.datetime64(st[0].stats.endtime.datetime), - ], - } + # Regular by construction from the stream's own sample rate, at ns + # resolution so a `to_stream` round trip preserves the coordinate. + t0 = np.datetime64(st[0].stats.starttime.datetime) + dt = np.rint(1e6 * st[0].stats.delta).astype("m8[us]").astype("m8[ns]") + time = Coordinate["interpolated"].from_block(t0, st[0].stats.npts, dt, dim=dims[1]) return DataArray(data, {dims[0]: channel, dims[1]: time}) diff --git a/xdas/io/prodml.py b/xdas/io/prodml.py index a71a0694..f0861c4c 100644 --- a/xdas/io/prodml.py +++ b/xdas/io/prodml.py @@ -7,8 +7,8 @@ import h5py import pandas as pd -from ..coordinates.core import Coordinate -from ..core.dataarray import DataArray +from ..coordinates import Coordinate +from ..core import DataArray from ..virtual import VirtualSource from .core import Engine @@ -48,11 +48,12 @@ def open_dataarray(self, fname, swapped_dims=False): else: nt, nd = data.shape - # time + # time (regular by declaration, rate derived from the file's own stamps) time = { "tie_indices": [0, nt - 1], "tie_values": [tstart, tend], - } # TODO: use from_block + "sampling_interval": (tend - tstart) / (nt - 1), + } # distance distance = Coordinate[self.ctype["distance"]].from_block( diff --git a/xdas/io/silixa.py b/xdas/io/silixa.py index d158d0ef..5f832056 100644 --- a/xdas/io/silixa.py +++ b/xdas/io/silixa.py @@ -3,8 +3,8 @@ import dask import numpy as np -from ..coordinates.core import Coordinate -from ..core.dataarray import DataArray +from ..coordinates import Coordinate +from ..core import DataArray from .core import Engine from .tdms import TdmsReader diff --git a/xdas/io/terra15.py b/xdas/io/terra15.py index 162ac324..23ad4c84 100644 --- a/xdas/io/terra15.py +++ b/xdas/io/terra15.py @@ -3,8 +3,8 @@ import h5py import pandas as pd -from ..coordinates.core import Coordinate -from ..core.dataarray import DataArray +from ..coordinates import Coordinate +from ..core import DataArray from ..virtual import VirtualSource from .core import Engine @@ -37,10 +37,12 @@ def open_dataarray(self, fname, tz="UTC"): dx = file.attrs["dx"] data = VirtualSource(file["data_product"]["data"]) nt, nd = data.shape + # time (regular by declaration, rate derived from the file's own stamps) time = { "tie_indices": [0, nt - 1], "tie_values": [ti, tf], - } # TODO: use from_block + "sampling_interval": (tf - ti) / (nt - 1), + } distance = Coordinate[self.ctype["distance"]].from_block( d0, nd, dx, dim="distance" ) diff --git a/xdas/io/xdas.py b/xdas/io/xdas.py index 609f7c5f..4461e3fd 100644 --- a/xdas/io/xdas.py +++ b/xdas/io/xdas.py @@ -14,9 +14,8 @@ from dask.array import Array as DaskArray from ..coordinates import Coordinates -from ..core.dataarray import DataArray -from ..core.datacollection import DataCollection, DataMapping, DataSequence -from ..dask.core import create_variable, loads +from ..core import DataArray, DataCollection, DataMapping, DataSequence +from ..dask import create_variable, loads from ..virtual import VirtualArray, VirtualSource from .core import Engine diff --git a/xdas/processing/core.py b/xdas/processing/core.py index fb3ad699..303d8b35 100644 --- a/xdas/processing/core.py +++ b/xdas/processing/core.py @@ -20,8 +20,7 @@ from watchdog.events import FileSystemEventHandler from watchdog.observers import Observer -from ..core.dataarray import DataArray -from ..core.routines import concat, open_dataarray +from ..core import DataArray, concat, open_dataarray from .monitor import Monitor @@ -438,6 +437,7 @@ class StreamWriter: ... "time": { ... "tie_indices": [0, data.shape[0] - 1], ... "tie_values": [starttime, endtime], + ... "sampling_interval": np.timedelta64(10, "ms"), ... }, ... "distance": distance, ... }, @@ -597,7 +597,7 @@ class ZMQPublisher: First we generate some data and split it into packets - >>> packets = xd.split(xd.synthetics.dummy(), 10) + >>> packets = xd.split(xd.testing.dummy(), 10) We initialize the publisher at a given address @@ -672,7 +672,7 @@ class ZMQSubscriber: First we generate some data and split it into packets - >>> da = xd.synthetics.dummy() + >>> da = xd.testing.dummy() >>> packets = xd.split(da, 10) We then publish the packets asynchronously diff --git a/xdas/signal.py b/xdas/signal.py index c93e2fe9..e9f11d53 100644 --- a/xdas/signal.py +++ b/xdas/signal.py @@ -8,9 +8,9 @@ import numpy as np import scipy.signal as sp -from .atoms.core import atomized -from .coordinates.core import Coordinate, get_sampling_interval -from .core.dataarray import DataArray +from .atoms import atomized +from .coordinates import get_sampling_interval +from .core import DataArray from .parallel import parallelize from .spectral import stft # noqa @@ -118,8 +118,10 @@ def filter(da, freq, btype, corners=4, zerophase=False, dim="last", parallel=Non """ axis = da.get_axis_num(dim) + dim = da.dims[axis] + d = get_sampling_interval(da, dim) across = int(axis == 0) - fs = 1.0 / get_sampling_interval(da, dim) + fs = 1.0 / d sos = sp.iirfilter(corners, freq, btype=btype, ftype="butter", output="sos", fs=fs) if zerophase: func = parallelize((None, across), across, parallel)(sp.sosfiltfilt) @@ -246,10 +248,11 @@ def resample(da, num, dim="last", window=None, domain="time", parallel=None): """ axis = da.get_axis_num(dim) dim = da.dims[axis] + get_sampling_interval(da, dim) # warn or raise on irregular axes upfront across = int(axis == 0) func = parallelize(across, across, parallel)(sp.resample) data, t = func(da.values, num, da[dim].values, axis, window, domain) - new_coord = {"tie_indices": [0, num - 1], "tie_values": [t[0], t[-1]]} + new_coord = type(da.coords[dim]).from_block(t[0], num, t[1] - t[0], dim=dim) coords = { name: new_coord if name == dim else coord for name, coord in da.coords.items() @@ -342,20 +345,13 @@ def resample_poly( """ axis = da.get_axis_num(dim) dim = da.dims[axis] + d = get_sampling_interval(da, dim, cast=False) across = int(axis == 0) func = parallelize(across, across, parallel)(sp.resample_poly) data = func(da.values, up, down, axis, window, padtype, cval) start = da[dim][0].values - d = da[dim][-1].values - da[dim][-2].values - end = da[dim][-1].values + d - new_coord = Coordinate( - { - "tie_indices": [0, data.shape[axis]], - "tie_values": [start, end], - }, - dim, - ) - new_coord = new_coord[:-1] + step = d * down / up + new_coord = type(da.coords[dim]).from_block(start, data.shape[axis], step, dim=dim) coords = { name: new_coord if name == dim else coord for name, coord in da.coords.items() @@ -795,6 +791,7 @@ def integrate(da, midpoints=False, dim="last", parallel=None): """ axis = da.get_axis_num(dim) + dim = da.dims[axis] d = get_sampling_interval(da, dim) def func(x): @@ -838,6 +835,7 @@ def differentiate(da, midpoints=False, dim="last", parallel=None): """ axis = da.get_axis_num(dim) + dim = da.dims[axis] d = get_sampling_interval(da, dim) def func(x): @@ -921,6 +919,7 @@ def sliding_mean_removal( """ axis = da.get_axis_num(dim) + dim = da.dims[axis] d = get_sampling_interval(da, dim) n = round(wlen / d) if n % 2 == 0: diff --git a/xdas/spectral.py b/xdas/spectral.py index 2e92d990..e8c88648 100644 --- a/xdas/spectral.py +++ b/xdas/spectral.py @@ -8,8 +8,8 @@ from scipy.fft import fft, fftfreq, fftshift, rfft, rfftfreq from scipy.signal import get_window -from .coordinates.core import get_sampling_interval -from .core.dataarray import DataArray +from .coordinates import get_sampling_interval +from .core import DataArray from .parallel import parallelize @@ -90,11 +90,12 @@ def stft( else: raise ValueError("Scaling must be 'spectrum' or 'psd'") scale = np.sqrt(scale) + coord_cls = type(da.coords[input_dim]) if return_onesided: freqs = rfftfreq(nfft, dt) else: freqs = fftshift(fftfreq(nfft, dt)) - freqs = {"tie_indices": [0, len(freqs) - 1], "tie_values": [freqs[0], freqs[-1]]} + freqs = coord_cls.from_block(freqs[0], len(freqs), 1.0 / (nfft * dt)) def func(x): """Apply windowed FFT to produce the STFT output array.""" @@ -123,11 +124,7 @@ def func(x): dt = get_sampling_interval(da, input_dim, cast=False) t0 = da.coords[input_dim].values[0] starttime = t0 + (nperseg / 2) * dt - endtime = starttime + (data.shape[axis] - 1) * (nperseg - noverlap) * dt - time = { - "tie_indices": [0, data.shape[axis] - 1], - "tie_values": [starttime, endtime], - } + time = coord_cls.from_block(starttime, data.shape[axis], (nperseg - noverlap) * dt) coords = {} for name in da.coords: diff --git a/xdas/synthetics.py b/xdas/synthetics.py index ecf593b8..5f76e37d 100644 --- a/xdas/synthetics.py +++ b/xdas/synthetics.py @@ -7,8 +7,8 @@ import numpy as np import scipy.signal as sp -from .core.dataarray import DataArray -from .core.routines import split +from .coordinates import Coordinate +from .core import DataArray, split def wavelet_wavefronts( @@ -84,14 +84,12 @@ def wavelet_wavefronts( da = DataArray( data=data, coords={ - "time": { - "tie_indices": [0, shape[0] - 1], - "tie_values": [starttime, starttime + resolution[0] * (shape[0] - 1)], - }, - "distance": { - "tie_indices": [0, shape[1] - 1], - "tie_values": [0.0, resolution[1] * (shape[1] - 1)], - }, + "time": Coordinate["interpolated"].from_block( + starttime, shape[0], resolution[0], dim="time" + ), + "distance": Coordinate["interpolated"].from_block( + 0.0, shape[1], resolution[1], dim="distance" + ), }, ) if nchunk is not None: @@ -145,42 +143,12 @@ def randn_wavefronts(): da = DataArray( data=data, coords={ - "time": { - "tie_indices": [0, shape[0] - 1], - "tie_values": [starttime, starttime + resolution[0] * (shape[0] - 1)], - }, - "distance": { - "tie_indices": [0, shape[1] - 1], - "tie_values": [0.0, resolution[1] * (shape[1] - 1)], - }, + "time": Coordinate["interpolated"].from_block( + starttime, shape[0], resolution[0], dim="time" + ), + "distance": Coordinate["interpolated"].from_block( + 0.0, shape[1], resolution[1], dim="distance" + ), }, ) return da - - -def dummy(shape=(1000, 100)): - """ - Return a minimal random :class:`DataArray` for quick testing. - - Parameters - ---------- - shape : tuple of int, optional - ``(n_time, n_distance)`` shape. Defaults to ``(1000, 100)``. - - Returns - ------- - DataArray - DataArray filled with Gaussian noise, sampled at 10 Hz over - ``[0, 1000]`` m with ``time`` starting at 2024-01-01. - """ - starttime = np.datetime64("2024-01-01T00:00:00.000000000") - endtime = starttime + (shape[0] - 1) * np.timedelta64(100, "ms") - time = {"tie_indices": [0, shape[0] - 1], "tie_values": [starttime, endtime]} - distance = {"tie_indices": [0, shape[1] - 1], "tie_values": [0.0, 1000.0]} - return DataArray( - data=np.random.randn(*shape), - coords={ - "time": time, - "distance": distance, - }, - ) diff --git a/xdas/testing.py b/xdas/testing.py new file mode 100644 index 00000000..f03f13db --- /dev/null +++ b/xdas/testing.py @@ -0,0 +1,80 @@ +"""Test utilities for xdas.""" + +import numpy as np + +from .coordinates import Coordinate +from .core import DataArray + + +def dummy( + dims=("time", "distance"), + shape=(100, 10), + dtype=float, + step=(0.01, 10.0), + ctype="interpolated", + datetime=True, +): + """ + Return a minimal :class:`DataArray` for quick testing. + + Parameters + ---------- + dims : tuple of str, optional + Dimension names. Length must match ``shape``. Defaults to + ``("time", "distance")``. + shape : tuple of int, optional + Size along each dimension. Defaults to ``(100, 10)``. + dtype : dtype-like, optional + Data type for the array values. Defaults to ``float``. + step : scalar or tuple, optional + Step size for each dimension. A single value is applied to all + dimensions; a tuple must have the same length as ``dims``. Defaults + to ``(0.01, 10.0)`` (100 Hz, 10 m spacing → 1 s × 100 m total). + When ``datetime=True``, a float step for the first dimension is + interpreted as seconds and converted to :class:`numpy.timedelta64`. + ctype : {"interpolated", "sampled", "dense"}, optional + Coordinate type for all dimensions. Defaults to ``"interpolated"``. + datetime : bool, optional + If ``True`` (default), the first dimension uses + :class:`numpy.datetime64` coordinates starting at 2024-05-21. + All other dimensions use float coordinates starting at 0.0. + + Returns + ------- + DataArray + Array filled with sequential integers (via :func:`numpy.arange`) + reshaped to ``shape`` and cast to ``dtype``. + + Examples + -------- + >>> import xdas as xd + >>> da = xd.testing.dummy() + >>> da.shape + (100, 10) + >>> da = xd.testing.dummy(dims=("x",), shape=(50,), datetime=False, step=1.0) + >>> da.shape + (50,) + >>> da = xd.testing.dummy(dims=("x",), shape=(10,), datetime=False, step=2.0) + >>> float(da.coords["x"].sampling_interval) + 2.0 + + """ + if len(dims) != len(shape): + raise ValueError(f"len(dims)={len(dims)} must equal len(shape)={len(shape)}") + if isinstance(step, (tuple, list)) and len(step) != len(dims): + raise ValueError(f"len(step)={len(step)} must equal len(dims)={len(dims)}") + + data = np.arange(int(np.prod(shape))).reshape(shape).astype(dtype) + + coords = {} + for i, (dim, size) in enumerate(zip(dims, shape)): + s = step[i] if isinstance(step, (tuple, list)) else step + if datetime and i == 0: + start = np.datetime64("2024-05-21T00:00:00.000000000") + if isinstance(s, (int, float)): + s = np.timedelta64(int(s * 1e9), "ns") + else: + start = 0.0 + coords[dim] = Coordinate[ctype].from_block(start, size, s, dim=dim) + + return DataArray(data=data, coords=coords) diff --git a/xdas/trigger.py b/xdas/trigger.py index 3ec37b7d..4f07fbac 100644 --- a/xdas/trigger.py +++ b/xdas/trigger.py @@ -10,8 +10,8 @@ from numba import njit from .atoms.core import Atom, State, atomized -from .coordinates.core import Coordinate -from .core.routines import concat_coords +from .coordinates import Coordinate +from .core import concat_coords class Trigger(Atom):