diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml index b67a3372..2170a002 100644 --- a/.github/workflows/tests.yaml +++ b/.github/workflows/tests.yaml @@ -8,7 +8,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] + python-version: ["3.11", "3.12", "3.13", "3.14"] steps: - uses: actions/checkout@v4 diff --git a/docs/api/index.md b/docs/api/index.md index 2f8a3e99..e2179dc9 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -14,5 +14,6 @@ processing signal synthetics testing +tiles virtual ``` \ No newline at end of file diff --git a/docs/api/tiles.md b/docs/api/tiles.md new file mode 100644 index 00000000..dec6c782 --- /dev/null +++ b/docs/api/tiles.md @@ -0,0 +1,65 @@ +```{eval-rst} +.. currentmodule:: xdas.virtual.tiles +``` + +# xdas.virtual.tiles + +Lazy tile-backed virtual arrays: the only backend of the formats that +HDF5 virtual datasets cannot serve (Silixa TDMS, MiniSEED), the default +one for Febus, and available on request from every other engine +(`vtype="tiles"`). + +## TileArray + +A dense rectilinear grid of file-backed tiles exposed as one lazy +numpy-like array. + +Attributes + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + TileArray.shape + TileArray.dtype + TileArray.ndim + TileArray.size + TileArray.chunks + TileArray.ntiles + TileArray.engine + TileArray.attrs +``` + +Methods + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + TileArray.from_tiles + TileArray.from_variable + TileArray.to_dataset + TileArray.create_variable + TileArray.finalize_save + TileArray.concat + TileArray.expand_dims + TileArray.squeeze + TileArray.transpose + TileArray.astype + TileArray.equals +``` + +## Engine lookup + +Tiles are decoded by the ``load_tile`` half of the +{class}`xdas.io.Engine` format plugins; the engine names stored in tile +manifests resolve on that registry (``Engine[name]``). + +```{eval-rst} +.. currentmodule:: xdas.io + +.. autosummary:: + :toctree: ../_autosummary + + Engine.load_tile +``` diff --git a/docs/api/virtual.md b/docs/api/virtual.md index 3c1e40a1..93104ba2 100644 --- a/docs/api/virtual.md +++ b/docs/api/virtual.md @@ -4,9 +4,20 @@ # xdas.virtual +## VirtualBackend + +Marker base and `vtype` registry for the virtual array backends. + +```{eval-rst} +.. autosummary:: + :toctree: ../_autosummary + + VirtualBackend +``` + ## VirtualArray -Base class for all virtual array types. +Base class of the HDF5 virtual dataset backend (`vtype="hdf5"`). Attributes @@ -28,6 +39,7 @@ Methods .. autosummary:: :toctree: ../_autosummary + VirtualArray.from_variable VirtualArray.to_dataset ``` diff --git a/docs/api/xdas.md b/docs/api/xdas.md index 8ebc62ed..8012fb94 100644 --- a/docs/api/xdas.md +++ b/docs/api/xdas.md @@ -32,6 +32,7 @@ concatenate concat_coords get_sampling_interval + sortby split plot_availability ``` diff --git a/docs/release-notes.md b/docs/release-notes.md index 6285f66b..b5cfc0e4 100644 --- a/docs/release-notes.md +++ b/docs/release-notes.md @@ -1,5 +1,28 @@ # Release notes +## 0.2.9 (unreleased) + +### New Features +- **Tile-backed virtual arrays.** The new `xdas.virtual.tiles` module exposes file archives as one lazy `TileArray`. Slicing (any step, including negative), integer indexing, `np.newaxis`, concatenation, and the numpy manipulation routines (the `transpose`, `flip`, `split`, `stack` and `atleast` families, `expand_dims`, `squeeze`, `roll`, `tile`, `delete`, `append`/`insert`) all stay lazy; whole-array reductions (`sum`, `mean`, `min`, `max`, …) stream one tile row at a time; reads touch only the tiles the selection overlaps (@atrabattoni). +- **`vtype="tiles"` on every HDF5 engine.** The open functions with `vtype="tiles"` return tile-backed arrays for the asn, febus, terra15, apsensing, prodml and native xdas engines. Silixa and MiniSEED always emit them now (replacing the serialized-dask-graph fallback, with time-axis push-down for Silixa), and Febus defaults to them — one tile per file, where the HDF5 backing needed one virtual mapping per data block. Custom engines add support by implementing `Engine.load_tile(path, selection, **params)` (@atrabattoni). +- Tile-backed arrays round-trip through the native xdas netCDF format: the manifest is stored as a compact `__tiles__` sibling group, relocatable by editing its single root path and directly readable by the 0.3 line (@atrabattoni). +- **Explicit engine configuration.** The open functions declare `engine`, `vtype` and `ctype`, and `engine` also accepts a configured `xdas.io.Engine` instance. Format-specific parameters (`overlaps`/`offset` for febus, `ignore_last_sample` for miniseed, `swapped_dims` for prodml, `tz` for terra15, `group` for the native format) are engine constructor parameters, validated up front (@atrabattoni). +- `open_mfdataarray` no longer caps the number of files when the resolved vtype consolidates its scan results, which `tiles` does; the 100 000 ceiling remains for `hdf5`, which builds one virtual mapping per file and so cannot be fused into anything smaller (@atrabattoni). +- **Streamed multi-file combining.** `open_mfdataarray` now fuses scan results every 100 000 files instead of holding one data array per file until the end, so memory no longer grows with the archive: results are accumulated without coordinate simplification (lossless in any arrival order) and sorted once at the end, giving the same result as before whatever the file naming. Acquisitions interleaved in time now group by compatibility, one array per acquisition, instead of splitting at each alternation. Since the batch size is also the ceiling, anything that opened in one call before still takes the single-batch path unchanged (@atrabattoni). +- **`xdas.sortby`.** Sort a tile- or stack-backed data array along a dimension by coordinate value, lazily: the blocks are permuted through the manifest without reading any data. This is how the streamed combine orders shuffled archives, exposed for standalone use (@atrabattoni). +- `simplify` runs in linear time whatever the number of gaps: the reduce stage is now a one-pass sleeve instead of Douglas-Peucker, which degenerated quadratically on gap-rich coordinates (a 100 000-file gappy archive simplified in minutes; now milliseconds). The deviation guarantee is unchanged — dropped tie points stay within `tolerance` of the curve, surviving values never move — though the surviving tie-point selection may differ slightly on jittery axes (@atrabattoni). + +### Breaking Changes +- Python 3.10 support is dropped and the numpy requirement is raised to 2.3: the tile manifests use `np.strings` routines introduced in numpy 2.3, which itself requires Python 3.11+. Python 3.10 reaches end of life in October 2026 (@atrabattoni). +- Passing a bare read function as `engine` now raises a `TypeError`: subclass `xdas.io.Engine` instead (see the data-formats documentation) (@atrabattoni). +- Misspelled or unsupported keyword arguments passed next to an engine name now raise a `TypeError` instead of being silently ignored, and combining `vtype`, `ctype` or engine keywords with an already configured engine instance raises a `ValueError` (@atrabattoni). + +### Deprecations +- Writing dask-backed virtual arrays is deprecated and emits a `FutureWarning`; existing files still open, but no engine emits them any more (@atrabattoni). + +### Bug Fixes +- The miniseed `ctype` argument is now honored: it previously routed to an unused attribute and the reader always built interpolated time coordinates. The default is unchanged (@atrabattoni). + ## 0.2.8 ### New Features diff --git a/docs/user-guide/coordinates/interpolated-coordinates.md b/docs/user-guide/coordinates/interpolated-coordinates.md index e5db746b..3dddbed6 100644 --- a/docs/user-guide/coordinates/interpolated-coordinates.md +++ b/docs/user-guide/coordinates/interpolated-coordinates.md @@ -84,8 +84,10 @@ Gaps represent missing data and are generally not problematic; overlaps usually arise from labelling errors and should be resolved. Using the {py:meth}`~xdas.coordinates.InterpCoordinate.simplify` method, -the coordinate can be compressed with controlled accuracy using the -[Ramer–Douglas–Peucker algorithm][RDP]. In the example below, the +the coordinate can be compressed with controlled accuracy using a +one-pass [sleeve algorithm][SDT] (tie points are dropped as long as the +segment joining the surviving neighbours passes within `tolerance` of +them, in a single left-to-right walk). In the example below, the second tie point carries no additional information and is safely discarded: ```{code-cell} @@ -162,4 +164,4 @@ coord.to_index(slice("2023-01-01T00:10:00", "2023-01-01T00:20:00")) ``` [CF]: -[RDP]: +[SDT]: diff --git a/docs/user-guide/io/data-formats.md b/docs/user-guide/io/data-formats.md index e32b327c..29f01c8f 100644 --- a/docs/user-guide/io/data-formats.md +++ b/docs/user-guide/io/data-formats.md @@ -20,75 +20,87 @@ os.chdir("../../_data") ## Implemented file formats -Here below the list of formats that are currently implemented. All HDF5 based formats support native virtualization. Other formats support Dask virtualization. Please refer to the [](virtual-datasets) section. Xdas should automatically detect the correct file format. You can still specify which one you want in the `engine` argument in {py:func}`xdas.open`. +Here below the list of formats that are currently implemented. Every format supports tile virtualization; HDF5 based formats also support native HDF5 virtualization, and the tables give the backing each engine uses when you do not ask for one. Pass `vtype` to {py:func}`xdas.open` to choose the other, and see [](virtual-datasets) for how to pick. Xdas should automatically detect the correct file format. You can still specify which one you want in the `engine` argument. Xdas support the following DAS formats: -| Constructor | Instrument | `engine` argument | Virtualization | -|:-----------------:|:-----------------:|:-----------------:|:-----------------:| -| AP Sensing | DAS N5* | `"apsensing"` | HDF5 | -| ASN | OptoDAS | `"asn"` | HDF5 | -| FEBUS | A1 | `"febus"` | HDF5 | -| OptaSense | OLA, ODH*, ... | `"optasense"` | HDF5 | -| Silixa | iDAS | `"silixa"` | Dask | -| SINTELA | ONYX | `"sintela"` | HDF5 | -| Terra15 | Treble | `"terra15"` | HDF5 | +| Constructor | Instrument | `engine` argument | Virtualization | Default | +|:-----------------:|:-----------------:|:-----------------:|:-----------------:|:---------:| +| AP Sensing | DAS N5* | `"apsensing"` | HDF5, tiles | `hdf5` | +| ASN | OptoDAS | `"asn"` | HDF5, tiles | `hdf5` | +| FEBUS | A1 | `"febus"` | HDF5, tiles | `tiles` | +| OptaSense | OLA, ODH*, ... | `"optasense"` | HDF5, tiles | `hdf5` | +| Silixa | iDAS | `"silixa"` | tiles | `tiles` | +| SINTELA | ONYX | `"sintela"` | HDF5, tiles | `hdf5` | +| Terra15 | Treble | `"terra15"` | HDF5, tiles | `hdf5` | It also implements its own format and support ProdML and miniSEED: -| Format | `engine` argument | Virtualization | -|:-----------------:|:-----------------:|:-----------------:| -| Xdas | `None` | HDF5 | -| ProdML | `"prodml"` | HDF5 | -| miniSEED | `"miniseed"` | Dask | +| Format | `engine` argument | Virtualization | Default | +|:-----------------:|:-----------------:|:-----------------:|:---------:| +| Xdas | `None` | HDF5, tiles | `hdf5` | +| ProdML | `"prodml"` | HDF5, tiles | `hdf5` | +| miniSEED | `"miniseed"` | tiles | `tiles` | + +```{note} +A Febus file stores a stack of overlapping blocks rather than one contiguous array. The +HDF5 backing needs a separate mapping for every block, so the manifest of a Febus dataset +grows with the number of blocks as well as the number of files; a tile array needs one +tile per file whatever the block count. That is why `febus` defaults to `tiles`. +``` ```{warning} Due to poor documentation of the various version of the Febus format, it is recommended to manually provide the required trimming and the position of the timestamps within each block. For example to trim 100 samples on both side of each block and to set the timestamp location at the center of the block for a block of 2000 samples: `xdas.open("path.h5", engine="febus", overlaps=(100, 100), offset=1000)` ``` -## Extending *xdas* with your file format - -*xdas* insists on its extensibility, the power is in the hands of the users. Extending *xdas* usually consists of writing few-line-of-code-long functions. The process consists in dealing with the two main aspects of a {py:class}`xarray.DataArray`: unpacking the data and coordinates objects, eventually processing them and packing them back into a Database object. +### Engine parameters -### Function-based solution +Every open function ({py:func}`xdas.open`, {py:func}`xdas.open_dataarray`, +{py:func}`xdas.open_mfdataarray`, {py:func}`xdas.open_mfdatatree`) takes the same +engine-related arguments. `engine` selects the file format; `vtype` and `ctype` +select the virtualization backend and the coordinate types, and exist for every +engine. Some formats take additional parameters (the trimming of Febus blocks +shown above, the timezone of Terra15 timestamps, ...). Those are engine +constructor parameters: when `engine` is given by name, any extra keyword +argument is forwarded to the engine constructor; alternatively you can configure +an engine instance yourself and pass it as `engine`. The three calls below are +equivalent: -To add a new file format the user can specify a function that read one file and outputs a {py:class}`xarray.DataArray`. This function can then be passed as an engine keyword argument to the {py:func}`xdas.open` function. The reading function must fetch and parse the data and coordinates information. - -Adding the support for a new file format generally consists in providing the path to the data array and parsing the start time and spatial and temporal spacing as in the example below. - -```{code-cell} -import h5py -import numpy as np +```python import xdas as xd -from xdas import DataArray -from xdas.virtual import VirtualSource +from xdas.io.febus import FebusEngine + +da = xd.open("path.h5", engine="febus", vtype="tiles", overlaps=(100, 100), offset=1000) +da = xd.open( + "path.h5", engine=FebusEngine(vtype="tiles", overlaps=(100, 100), offset=1000) +) +engine = FebusEngine(vtype="tiles", overlaps=(100, 100), offset=1000) # reusable +da = xd.open("path.h5", engine=engine) +``` -def open_dataarray(fname): - with h5py.File(fname, "r") as file: - t0 = np.datetime64(file["dataset"].attrs["t0"]).astype("datetime64[ms]") - dt = np.timedelta64(int(file["dataset"].attrs["dt"]*1e3), "ms") - dx = file["dataset"].attrs["dx"][()] - data = VirtualSource(file["dataset"]) - nt, nx = data.shape - t = {"tie_indices": [0, nt - 1], "tie_values": [t0, t0 + (nt - 1) * dt]} - x = {"tie_indices": [0, nx - 1], "tie_values": [0.0, (nx - 1) * dx]} - return DataArray(data, {"time": t, "distance": x}) +Misspelled or unsupported parameters raise a `TypeError` from the engine +constructor. A configured instance is a complete specification: combining it +with `vtype`, `ctype` or extra keyword arguments raises an error. Format +auto-detection (`engine=None`) accepts `vtype` and `ctype` but no +format-specific parameters, since those require knowing the format. -# Replace "other_format.hdf5" by the path of your file -da = xd.open("other_format.hdf5", engine=open_dataarray) -da -``` +## Extending *xdas* with your file format -This example is for one file. For multi-file datasets please indicate the path of your files with a '*' before the file format if all your files are in the same folder or pass a list of paths. +*xdas* insists on its extensibility, the power is in the hands of the users. Extending *xdas* usually consists of writing a few-line-of-code-long engine class. The process consists in dealing with the two main aspects of a {py:class}`xarray.DataArray`: unpacking the data and coordinates objects, eventually processing them and packing them back into a Database object. -### Class-based solution +### Writing an engine -To add support in a more complete way, you can also create your own engine by inheriting from the `xdas.io.Engine` abstract class. Note that when the class is defined, the `name` keyword argument allows to register the new engine along with the `aliases` one that is useful when several instruments share the same data format. This allows to add your engine to the `Engine._registry` and to retrieve it by doing `Engine[name]`. The `_supported_vtypes` and `_supported_ctypes` class attributes allow to determine which kind of virtualization backend and type of coordinates can be used with this file format. When you open any file, you can additionally provide the `vtype` and `ctype` keyword arguments to specify which backends to use. The `Engine` class defines the `__init__` method that checks those passed kwargs and stores in `self.vtype` and `self.ctype` the chosen backends. +To add a new file format, create your own engine by inheriting from the `xdas.io.Engine` abstract class. Note that when the class is defined, the `name` keyword argument allows to register the new engine along with the `aliases` one that is useful when several instruments share the same data format. This allows to add your engine to the `Engine._registry` and to retrieve it by doing `Engine[name]`. The `_supported_vtypes` and `_supported_ctypes` class attributes allow to determine which kind of virtualization backend and type of coordinates can be used with this file format. When you open any file, you can additionally provide the `vtype` and `ctype` keyword arguments to specify which backends to use. The `Engine` class defines the `__init__` method that checks those passed kwargs and stores in `self.vtype` and `self.ctype` the chosen backends. If your format needs parameters of its own, define an `__init__` taking them after `vtype` and `ctype` and calling `super().__init__(vtype, ctype)`: they then become available next to the engine name in the open functions, like the built-in ones described above. ```{code-cell} -from xdas.io import Engine +import h5py +import numpy as np +import xdas as xd +from xdas import DataArray from xdas.coordinates import Coordinate +from xdas.io import Engine +from xdas.virtual import VirtualSource class MyEngine(Engine, name="my_engine", aliases=["other_engine"]): _supported_vtypes = ["hdf5"] @@ -117,3 +129,59 @@ Once the class is created and instanciated you can then use it : da = xd.open("other_format.hdf5", engine="my_engine", ctype="sampled") da ``` + +### Tile-backed engines + +Beside the `hdf5` vtype shown above (an HDF5 virtual source), an engine can +offer the `tiles` vtype: `open_dataarray` then backs the data array with a lazy +{py:class}`xdas.virtual.TileArray` describing the file, and the engine implements +the decoding half as a `load_tile` static method — called once per tile +touched, with exactly one source-local slice per source axis, in source order +(whatever transposes or inserted axes the tile array presents), and the +manifest's engine specification as keyword arguments, returning exactly the +selected sub-box: + +```{code-cell} +from xdas.virtual import TileArray + +class MyTileEngine(Engine, name="my_tile_engine"): + _supported_vtypes = ["hdf5", "tiles"] + _supported_ctypes = { + "distance": ["interpolated", "sampled", "dense"], + "time": ["interpolated", "sampled", "dense"], + } + + def open_dataarray(self, fname): + with h5py.File(fname, "r") as file: + t0 = np.datetime64(file["dataset"].attrs["t0"]).astype("datetime64[ms]") + dt = np.timedelta64(int(file["dataset"].attrs["dt"]*1e3), "ms") + x0 = file["dataset"].attrs["x0"][()] + dx = file["dataset"].attrs["dx"][()] + if self.vtype == "tiles": + data = TileArray.from_tiles( + str(fname), + file["dataset"].shape, + file["dataset"].dtype, + "my_tile_engine", + ) + else: + data = VirtualSource(file["dataset"]) + nt, nx = data.shape + t = Coordinate[self.ctype["time"]].from_block(t0, nt, dt, dim="time") + x = Coordinate[self.ctype["distance"]].from_block(x0, nx, dx, dim="distance") + return DataArray(data, {"time": t, "distance": x}) + + @staticmethod + def load_tile(path, selection): + with h5py.File(path, "r") as file: + return file["dataset"][selection] +``` + +`load_tile` must depend only on its arguments — never on engine instance +state — so that saved tile views decode identically everywhere. + +This is the backing used by default for the formats HDF5 virtual datasets cannot +serve (Silixa TDMS, MiniSEED) and for Febus, whose files hold many blocks each; it +is available on request from every other built-in engine. The order of +`_supported_vtypes` decides the default, so listing `"tiles"` first is all it +takes for a new engine to prefer it. diff --git a/docs/user-guide/io/virtual-datasets.md b/docs/user-guide/io/virtual-datasets.md index 6830a05c..2a205161 100644 --- a/docs/user-guide/io/virtual-datasets.md +++ b/docs/user-guide/io/virtual-datasets.md @@ -16,10 +16,14 @@ os.chdir("../../_data") To deal with large multi-file dataset, *Xdas* uses the concept of virtual datasets. A virtual dataset is a file that contains pointers towards an arbitrary number of files that can then be accessed seamlessly as a single, contiguous dataset. -*Xdas* uses two types of virtualization: +*Xdas* uses several types of virtualization, selected with the `vtype` argument: -- For HDF5 based format, it leverages the performance offered by the [virtual datasets](https://docs.h5py.org/en/stable/vds.html) native capabilities of netCDF4/HDF5 which comes with almost no overhead (C compiled). -- For other type of files, it leverage the flexibility of [Dask arrays](https://docs.Dask.org/en/stable/array.html). +- `hdf5`: for HDF5 based formats, it leverages the [virtual datasets](https://docs.h5py.org/en/stable/vds.html) native capabilities of netCDF4/HDF5, which resolve the mapping in compiled C with no per-file Python overhead. The cost of that mapping does however grow with the number of linked files. +- `tiles`: a manifest of file-backed tiles stored as a plain array, decoded by the engine itself. It works with any format and keeps the file mapping inspectable. + +Which types an engine offers is declared by its `_supported_vtypes` attribute; the first one listed is the default. See [](#choosing-a-virtualization-backend) for how to pick, and [](data-formats.md) for what each engine supports and defaults to. + +A third backing, [Dask arrays](https://docs.Dask.org/en/stable/array.html), is deprecated and no longer used by any engine — see [](#dask-virtualization). ## HDF5 Virtualization @@ -83,10 +87,135 @@ A virtual dataset can point to another virtual dataset. This can be beneficial f When loading large part of a virtual dataset, you might end up with nan values. This normally happens when linked files are missing. But due to a [known limitation](https://forum.hdfgroup.org/t/virtual-datasets-and-open-file-limit/6757) of the HDF5 C library it can be due to the opening of too many files. Try increasing the number of possible file to open with the `ulimit` command. Or load smaller chunk of data. ``` -## Dask Virtualization +(tile-virtualization)= +## Tile Virtualization + +With the `tiles` vtype, the mapping is not delegated to HDF5. *Xdas* stores it as a +{py:class}`xdas.virtual.TileArray`: a plain array manifest that records, for each tile, which +file it comes from and which part of that file it contributes. Reading a region resolves +which tiles it touches and asks the engine to decode each of them through its `load_tile` +method. The manifest is ordinary data, so it can be inspected, sliced and concatenated +like any other array, and it is stored as such inside the *Xdas* netCDF format. + +(choosing-a-virtualization-backend)= +## Choosing a virtualization backend + +For formats that HDF5 virtual datasets cannot serve, the choice is made for you. When an +engine supports both, the trade-off is essentially *who resolves the mapping*: the HDF5 C +library, or *Xdas* itself. + +That choice decides how each cost scales with the size of the archive. Writing and +reopening an HDF5 virtual dataset both cost one operation per linked file, so both grow +with the file count; a tile manifest is written and read back as an array, so neither +does. Reading inverts the expectation one might have: resolving a region inside the C +library involves no Python at all, but its cost grows with how many mappings the dataset +*contains*, while a tile manifest is searched, so its cost grows only with how many tiles +the read *touches*. Modest file counts therefore favour HDF5, and the advantage moves to +tiles as the archive grows. + +### HDF5 virtualization + +**Advantages** + +- Resolution happens inside the HDF5 C library, so reading involves no per-file Python + call. On modest file counts this makes it the faster of the two to read. +- Any HDF5-aware tool can read the result, not only *Xdas*. +- Virtual datasets can point at other virtual datasets, so a growing archive can be + linked in batches without relinking everything. +- A subset saved from a virtual dataset is very compact, because it refers to the dataset + it was cut from rather than restating the underlying file list. + +**Limitations** + +- Building the mapping costs one HDF5 call per source file, so both the time and the + memory needed to write a manifest grow in proportion to the number of files. Beyond + some point, writing a single flat manifest stops being practical. +- Reopening a virtual dataset reads its whole mapping table, so opening cost also grows + with the number of linked files. Deep archives therefore tend to require a pyramid of + virtual datasets, which shifts that cost to read time and multiplies the number of + manifest files to keep track of. The top of such a pyramid opens quickly precisely + because it defers the work: the first read of a region then has to open the level below + it, a toll that a short-lived process pays on every run. +- Read latency grows with the number of mappings the dataset holds, not only with the + amount of data asked for, so the same request gets slower as the archive it lives in + gets bigger. Past a large enough file count this outweighs the advantage of resolving + in C, and reads become slower than the tile equivalent. +- Once written, the mapping is opaque: HDF5 presents a virtual dataset as a regular + dataset, so the list of linked files can no longer be inspected or edited. +- Because a saved subset refers to its parent, extracts are not self-contained. Moving or + deleting the parent breaks them, and each extract adds one more level of indirection. +- Strided (decimating) selection along the concatenation axis is not supported. +- Missing files are read as NaN rather than raising, and exceeding the C library's + open-file limit produces the same symptom, which makes such problems easy to miss. + +### Tile virtualization + +**Advantages** + +- Writing a manifest is an array write rather than a per-file operation, so it stays fast + and light as the number of files grows. A single flat manifest remains workable at + scales where an HDF5 one does not. +- Opening loads only the tile geometry, so open time stays low even for very large + manifests, and no cost is deferred to the first read. +- Read latency depends on how much of the manifest a request touches, not on how large + the manifest is, so reads do not get slower as the archive grows. +- The file list is data: it can be inspected, modified and saved again. +- Concatenation fuses manifests without reading any values. +- A saved subset names the data files it needs directly, so extracts are self-contained + and never gain an extra level of indirection. +- Decimating along the stacking axis is supported: the stride is folded into the tile + geometry and stays lazy. +- The per-file footprint of the manifest is smaller. + +**Limitations** + +- Decoding goes through Python, one call per tile touched. A request spread over a great + many tiles therefore carries a per-tile overhead that the C library avoids, which is + what makes HDF5 the quicker reader while file counts stay modest. +- For very small manifests the result can be larger than the HDF5 equivalent, since paths + and geometry are written out explicitly instead of referring to a parent dataset. +- The manifest is only meaningful to *Xdas*. +- The engine must implement `load_tile`. + +### Considerations that apply to both + +- Building either manifest starts with reading the metadata of every file. That scan is + dominated by disk access and is usually the bulk of the total build time, so it is not a + criterion for choosing between the two. +- A scan holds at most 100 000 results at once. Tiles combines them every 100 000 files + into a manifest whose per-file cost is negligible, so the batch is freed and a single + call can scan any number of files. Fusing changes nothing for `hdf5`, which keeps one + HDF5 virtual mapping per file whatever it does, so 100 000 is a ceiling for it instead + of a batch size. +- The combined result does not depend on the scan order: files are sorted by their + coordinate values, not by their names. +- Neither backend helps when a coordinate is not monotonic — for instance when files + overlap in time. Label-based selection then falls back to a slow path in both cases, and + is better addressed in the data itself. + +```{hint} +As a rule of thumb, prefer `hdf5` when the dataset is modest in file count or when other +HDF5-based tooling has to read it, and `tiles` when the file count is large, when the +mapping needs to remain inspectable, when extracts must stand on their own, or when +decimated reads matter. The larger the archive, the stronger the case for `tiles`: it is +the only one of the two whose write, open and read costs do not all grow with the number +of files. +``` -Other type of formats will be loaded as Dask arrays. Those latter are a N-dimensional stack of chunks. At each chunk is associated a task to complete to get the values of that chunk. It results in a computation graph that Xdas is capable to serialize and store within its native NetCDF format. To be able to serialize the graph, it must only contain xdas or Dask functions. +(dask-virtualization)= +## Dask Virtualization (deprecated) -From an user point of view the use of this type of virtualization is very similar to HDF5 one. +```{deprecated} 0.2.9 +Dask virtualization is no longer used by any engine and writing it is deprecated. The +formats that once relied on it — those HDF5 virtual datasets cannot serve — now use +[tile virtualization](#tile-virtualization) instead. Existing files that store a Dask +graph can still be read, so nothing on disk is lost, but new datasets should not be +written this way. +``` -The main difference is that when opening a dataset with Dask virtualization, the entire graph of pointers to the files is loaded, can be modified and saved again. In the HDF5 case, opening a virtual dataset is handled the same way as if it is a regular file meaning that the underlying mapping of pointers is hidden and cannot be modified. Dask graph can be slow when they start to become very big (more than one million tasks). +Formats that HDF5 could not virtualize used to be loaded as Dask arrays: an +N-dimensional stack of chunks, each with a task attached that produces its values, +serialized into the native *Xdas* netCDF format as a computation graph. Tiles replace it +with a manifest that describes the same mapping as plain array data, which is both more +compact and far quicker to build, and which does not go sluggish once the graph reaches +millions of tasks. diff --git a/pyproject.toml b/pyproject.toml index d908894d..321a3ec6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "setuptools.build_meta" [project] name = "xdas" dynamic = ["version"] -requires-python = ">= 3.10" +requires-python = ">= 3.11" authors = [ { name = "Alister Trabattoni", email = "alister.trabattoni@gmail.com" }, ] @@ -17,7 +17,7 @@ dependencies = [ "loky", "msgpack", "numba", - "numpy", + "numpy>=2.3", "obspy", "pandas", "plotly", diff --git a/tests/coordinates/test_interp.py b/tests/coordinates/test_interp.py index 562badc4..610925e9 100644 --- a/tests/coordinates/test_interp.py +++ b/tests/coordinates/test_interp.py @@ -392,7 +392,7 @@ def test_simplify_multiple_runs_and_isolated_point(self): ) def test_simplify_keeps_kink(self): - # A genuine kink inside a continuous area forces Douglas-Peucker to keep + # A genuine kink inside a continuous area forces the reduction to keep # the deviating interior point. coord = InterpCoordinate( {"tie_indices": [0, 5, 10], "tie_values": [0.0, 100.0, 0.0]} @@ -1122,7 +1122,7 @@ def test_widen_only_when_needed(self): assert result._is_valid_sampling_interval(s, result.tolerance) def test_widening_beyond_the_budget_never_raises(self): - # Douglas-Peucker bounds how far values move, not how much drift fusing + # The reduction bounds how far values move, not how much drift fusing # a discontinuity exposes, so the required tolerance can exceed the # budget. Real OptoDAS seams: 2 ms late every 10 s at 125 Hz. t0 = np.datetime64("2021-10-27T15:44:10.721999872", "ns") @@ -1215,3 +1215,74 @@ def test_empty(self): coord = InterpCoordinate.from_block(0.0, 0, 2.0, dim="x") assert coord.empty assert coord.sampling_interval == 2.0 + + +class TestSleeve: + """The one-pass reduce stage: same deviation guarantee as the former + Douglas-Peucker, O(n) whatever survives.""" + + def test_deviation_bound_holds_on_jitter(self): + rng = np.random.default_rng(0) + n = 200 + starts = np.arange(n, dtype="int64") * 10 + jitter = rng.integers(-1_500_000, 1_500_000, n) + t0 = np.datetime64("2026-01-01", "ns").astype("i8") + seg = t0 + np.arange(n) * 100_000_000 + jitter + tie_indices = np.empty(2 * n, dtype="int64") + tie_indices[0::2] = starts + tie_indices[1::2] = starts + 9 + tie_values = np.empty(2 * n, dtype="i8") + tie_values[0::2] = seg + tie_values[1::2] = seg + 90_000_000 + coord = InterpCoordinate( + {"tie_indices": tie_indices, "tie_values": tie_values.astype("M8[ns]")}, + "time", + ) + tolerance = np.timedelta64(1_000_000, "ns") + result = coord.simplify(tolerance) + # dropped points stay within tolerance of the simplified curve, and + # surviving values are never moved + deviation = np.abs(result._get_value(coord.tie_indices) - coord.tie_values) + assert deviation.max() <= tolerance + kept = np.isin(result.tie_indices, coord.tie_indices) + assert kept.all() + + def test_every_gap_survives(self): + n = 50 + tie_indices = np.arange(2 * n, dtype="int64") + tie_indices[1::2] = tie_indices[0::2] + 1 + tie_indices = np.cumsum(np.where(np.arange(2 * n) % 2, 1, 9)) + tie_indices -= tie_indices[0] + values = np.arange(2 * n) * 1_000_000_000 + coord = InterpCoordinate( + {"tie_indices": tie_indices, "tie_values": values.astype("M8[ns]")}, + "time", + ) + result = coord.simplify(np.timedelta64(1, "ms")) + assert len(result.tie_indices) == len(coord.tie_indices) + + def test_subunit_tolerance_is_not_truncated(self): + # microsecond values with a nanosecond tolerance: a 1 us seam jitter + # needs a sub-microsecond budget to survive (400 ns keeps it) and a + # 1100 ns one to fuse — both unrepresentable in truncated us + values = np.array([0, 999, 1001, 2000], dtype="M8[us]") + coord = InterpCoordinate( + {"tie_indices": [0, 999, 1000, 1999], "tie_values": values}, "time" + ) + kept = coord.simplify(np.timedelta64(400, "ns")) + fused = coord.simplify(np.timedelta64(1100, "ns")) + assert len(kept.tie_indices) == 4 + assert len(fused.tie_indices) == 2 + + def test_float_values(self): + coord = InterpCoordinate( + {"tie_indices": [0, 9, 10, 19], "tie_values": [0.0, 9.0, 10.5, 19.5]}, + "x", + ) + assert len(coord.simplify(1.0).tie_indices) == 2 + assert len(coord.simplify(0.1).tie_indices) == 4 + + def test_two_ties_pass_through(self): + coord = InterpCoordinate({"tie_indices": [0, 9], "tie_values": [0.0, 9.0]}) + result = coord.simplify(0.0) + assert len(result.tie_indices) == 2 diff --git a/tests/io/test_febus.py b/tests/io/test_febus.py index 5defdd01..85502e68 100644 --- a/tests/io/test_febus.py +++ b/tests/io/test_febus.py @@ -45,17 +45,13 @@ def test_open_with_freqres_attr(self, tmp_path): da = xd.open(str(path), engine="febus", overlaps=(1, 1), offset=0) assert isinstance(da, xd.DataArray) - def test_invalid_overlaps_raises(self, tmp_path): - path = tmp_path / "febus.h5" - make_febus_file(path) + def test_invalid_overlaps_raises(self): with pytest.raises(ValueError, match="overlaps must be"): - FebusEngine().open_dataarray(str(path), overlaps="bad") + FebusEngine(overlaps="bad") - def test_invalid_offset_raises(self, tmp_path): - path = tmp_path / "febus.h5" - make_febus_file(path) + def test_invalid_offset_raises(self): with pytest.raises(ValueError, match="offset must be an integer"): - FebusEngine().open_dataarray(str(path), overlaps=(1, 1), offset="bad") + FebusEngine(overlaps=(1, 1), offset="bad") def test_missing_block_rate_raises(self, tmp_path): path = tmp_path / "febus_no_blockrate.h5" @@ -70,4 +66,4 @@ def test_missing_block_rate_raises(self, tmp_path): zone.attrs["Origin"] = np.array([0.0, 0.0]) zone.create_dataset("Data", data=np.zeros((nchunks, nt, nx))) with pytest.raises(KeyError, match="Could not find the block size"): - FebusEngine().open_dataarray(str(path), overlaps=(0, 0), offset=0) + FebusEngine(overlaps=(0, 0), offset=0).open_dataarray(str(path)) diff --git a/tests/io/test_generic.py b/tests/io/test_generic.py index 0861e3d7..28f44a95 100644 --- a/tests/io/test_generic.py +++ b/tests/io/test_generic.py @@ -9,13 +9,25 @@ class TestEngineRegistry: def test_unknown_engine_raises_key_error(self): - with pytest.raises(KeyError, match="not found"): + with pytest.raises(KeyError, match="no engine registered"): Engine["nonexistent_engine_xyz"] def test_invalid_vtype_raises_value_error(self): with pytest.raises(ValueError, match="vtype must be None or a string"): Engine["asn"](vtype=42) + def test_unregistered_vtype_raises_key_error(self): + with pytest.raises(KeyError, match="no virtual backend registered"): + Engine["asn"](vtype="netcdf") + + def test_auto_engine_validates_vtype_upfront(self): + with pytest.raises(KeyError, match="no virtual backend registered"): + AutoEngine(vtype="netcdf") + + def test_registered_but_unsupported_vtype_still_engine_checked(self): + with pytest.raises(NotImplementedError, match="not supported by"): + Engine["silixa"](vtype="hdf5") + def test_dict_ctype_fills_missing_keys(self): engine = Engine["asn"](ctype={"time": "interpolated"}) assert engine.ctype["time"] == "interpolated" diff --git a/tests/io/test_miniseed.py b/tests/io/test_miniseed.py index 34321d98..71c10f7a 100644 --- a/tests/io/test_miniseed.py +++ b/tests/io/test_miniseed.py @@ -1,9 +1,12 @@ import numpy as np +import numpy.testing as npt import obspy import pytest import xdas as xd +from xdas.coordinates import Coordinate from xdas.io.miniseed import MiniSEEDEngine, get_band_code, to_stream +from xdas.virtual import TileArray def make_network(dirpath, gap=False, samples=100): @@ -79,6 +82,10 @@ def test_miniseed(tmp_path): assert da.coords["location"].values == "00" assert da.coords["channel"].values.tolist() == ["HHZ", "HHN", "HHE"] + # the ctype engine parameter drives the time coordinate flavor + da = xd.open(paths[0], engine="miniseed", ctype="dense") + assert isinstance(da.coords["time"], Coordinate["dense"]) + # read one file with gaps make_network(tmp_path, gap=True, samples=100) paths = sorted(tmp_path.glob("*_gap.mseed")) @@ -176,6 +183,42 @@ def test_miniseed(tmp_path): assert values_gap_trimmed.shape == (3, 89) +def test_miniseed_tile_backend(tmp_path): + make_network(tmp_path, samples=100) + + # single files open lazily, tile-backed + paths = sorted(tmp_path.glob("*00.mseed")) + da = xd.open(paths[0], engine="miniseed") + assert isinstance(da.data, TileArray) + assert da.data.engine == { + "name": "miniseed", + "method": "synchronized", + "ignore_last_sample": False, + } + + # concatenation along a new dimension stays lazy and reads correctly + objs = [xd.open(path, engine="miniseed") for path in paths] + stacked = xd.concat(objs, "station") + assert isinstance(stacked.data, TileArray) + npt.assert_array_equal(stacked.values, np.stack([obj.values for obj in objs])) + + # single-trace files fold the channel axis to a scalar coordinate + single = tmp_path / "single.mseed" + st = obspy.Stream([obspy.Trace(np.random.rand(50), make_header(1, "Z", 0))]) + st.write(str(single)) + da = xd.open(single, engine="miniseed") + assert da.dims == ("time",) + assert da.shape == (50,) + npt.assert_allclose(da.values, st[0].data, rtol=1e-6) + + # round trip through the native format + da = xd.open(paths[0], engine="miniseed") + da.to_netcdf(tmp_path / "view.nc") + reopened = xd.open_dataarray(tmp_path / "view.nc") + assert isinstance(reopened.data, TileArray) + npt.assert_array_equal(reopened.values, da.values) + + def test_miniseed_helpers(tmp_path): # get_band_code with out-of-range sampling rate assert get_band_code(0.0) == "X" @@ -204,4 +247,4 @@ def test_miniseed_unsynchronized_traces(tmp_path): ) st.write(str(path), format="MSEED") with pytest.raises(ValueError, match="synchronized"): - MiniSEEDEngine().read_header(str(path), False, "interpolated") + MiniSEEDEngine().read_header(str(path)) diff --git a/tests/io/test_silixa.py b/tests/io/test_silixa.py new file mode 100644 index 00000000..68b55eaf --- /dev/null +++ b/tests/io/test_silixa.py @@ -0,0 +1,65 @@ +import numpy as np +import numpy.testing as npt + +from xdas.io import silixa +from xdas.virtual import TileArray + + +class FakeTdms: + """In-memory stand-in for :class:`~xdas.io.tdms.TdmsReader`.""" + + data = np.arange(20.0 * 4).reshape(20, 4) + + channel_length = 20 + fileinfo = {"n_channels": 4} + _data_type = np.dtype("float64") + + def __init__(self, path): + self.path = path + + def __enter__(self): + return self + + def __exit__(self, *args): + return False + + def get_properties(self): + return { + "GPSTimeStamp": "2020-01-01T00:00:00", + "SamplingFrequency[Hz]": 1000.0, + "Start Distance (m)": 0.0, + "Fibre Length Multiplier": 1.0, + "SpatialResolution[m]": 4.0, + } + + def get_data(self, first_s=None, last_s=None): + first_s = 0 if first_s is None else first_s + last_s = len(self.data) - 1 if last_s is None else last_s + return self.data[first_s : last_s + 1] + + +def test_tile_load(monkeypatch): + monkeypatch.setattr(silixa, "TdmsReader", FakeTdms) + expected = FakeTdms.data + manifest = TileArray.from_tiles("fake.tdms", (20, 4), "float64", {"name": "silixa"}) + npt.assert_array_equal(np.asarray(manifest), expected) + npt.assert_array_equal(np.asarray(manifest[3:15:2, 1:3]), expected[3:15:2, 1:3]) + expanded = np.expand_dims(manifest, 0) + assert isinstance(expanded, TileArray) + npt.assert_array_equal(np.asarray(expanded), expected[None]) + + +def test_read_data(monkeypatch): + monkeypatch.setattr(silixa, "TdmsReader", FakeTdms) + npt.assert_array_equal(silixa.SilixaEngine().read_data("fake.tdms"), FakeTdms.data) + + +def test_open_dataarray(monkeypatch): + monkeypatch.setattr(silixa, "TdmsReader", FakeTdms) + da = silixa.SilixaEngine().open_dataarray("fake.tdms") + assert isinstance(da.data, TileArray) + assert da.dims == ("time", "distance") + assert da.shape == (20, 4) + assert da.coords["time"][0].values == np.datetime64("2020-01-01T00:00:00") + npt.assert_allclose(da.coords["distance"].values, 4.0 * np.arange(4)) + npt.assert_array_equal(da.values, FakeTdms.data) diff --git a/tests/io/test_tiles_vtype.py b/tests/io/test_tiles_vtype.py new file mode 100644 index 00000000..f6e9c3ff --- /dev/null +++ b/tests/io/test_tiles_vtype.py @@ -0,0 +1,294 @@ +"""Cross-format tests of the optional ``tiles`` vtype of the io engines.""" + +import h5py +import numpy as np +import numpy.testing as npt +import pytest + +import xdas as xd +from xdas.virtual import TileArray + + +def ramp(shape, dtype="float32"): + """Distinct values everywhere, so misplaced reads cannot cancel out.""" + return np.arange(np.prod(shape), dtype=dtype).reshape(shape) + + +def make_asn_file(path, nt=20, nd=4, t0=0.0): + data = ramp((nt, nd)) + with h5py.File(path, "w") as file: + header = file.create_group("header") + header["time"] = t0 + header["dt"] = 0.1 + header["dx"] = 10.0 + file.create_dataset("data", data=data) + cable_spec = file.create_group("cableSpec") + cable_spec["sensorDistances"] = 10.0 * np.arange(nd) + demod_spec = file.create_group("demodSpec") + demod_spec["roiStart"] = np.array([0]) + demod_spec["roiEnd"] = np.array([nd]) + return data + + +def make_febus_file(path, nchunks=3, nt=12, nx=5): + data = ramp((nchunks, nt, nx)) + times = np.arange(nchunks, dtype=np.float64) * 0.01 + with h5py.File(path, "w") as file: + source = file.create_group("DeviceName").create_group("Source1") + source.create_dataset("time", data=times) + zone = source.create_group("Zone1") + zone.attrs["BlockRate"] = np.array([100.0]) + zone.attrs["Spacing"] = np.array([5.0, 1.0]) + zone.attrs["Extent"] = np.array([0.0, (nx - 1) * 5.0]) + zone.attrs["Origin"] = np.array([0.0, 0.0]) + zone.create_dataset("StrainRate", data=data) + return data + + +def make_terra15_file(path, nt=15, nd=6): + data = ramp((nt, nd)) + with h5py.File(path, "w") as file: + product = file.create_group("data_product") + # small epoch offsets stay exact in float64 down to the nanosecond + product.create_dataset("gps_time", data=0.001 * np.arange(nt)) + product.create_dataset("data", data=data) + file.attrs["sensing_range_start"] = 12.0 + file.attrs["dx"] = 2.0 + return data + + +def make_apsensing_file(path, nt=10, nd=8): + data = ramp((nt, nd)) + with h5py.File(path, "w") as file: + file.create_dataset("DAS", data=data) + meta = file.create_group("Metadata") + meta.create_dataset("Timestamp", data=np.bytes_(b"2020-01-01T00:00:00.000Z")) + proc = file.create_group("ProcessingServer") + proc["DataRate"] = 1000.0 + proc["SpatialSampling"] = 2.0 + file.create_group("DAQ")["PositionStart"] = 0.0 + return data + + +def make_prodml_file(path, nt=10, nd=5, swapped=False): + data = ramp((nd, nt) if swapped else (nt, nd)) + with h5py.File(path, "w") as file: + acquisition = file.create_group("Acquisition") + acquisition.attrs["SpatialSamplingInterval"] = 2.0 + acquisition.attrs["StartLocusIndex"] = 0 + rawdata = acquisition.create_group("Raw[0]").create_dataset( + "RawData", data=data + ) + rawdata.attrs["PartStartTime"] = np.bytes_(b"2020-01-01T00:00:00.000+00:00") + rawdata.attrs["PartEndTime"] = np.bytes_(b"2020-01-01T00:00:00.900+00:00") + return data + + +MAKERS = { + "asn": (make_asn_file, {}), + "febus": (make_febus_file, {"overlaps": (1, 1), "offset": 0}), + "terra15": (make_terra15_file, {}), + "apsensing": (make_apsensing_file, {}), + "prodml": (make_prodml_file, {}), +} + + +@pytest.mark.parametrize("fmt", sorted(MAKERS)) +def test_tiles_vtype_matches_hdf5(tmp_path, fmt): + """The opt-in tiles backing yields the exact same array as the VDS one.""" + maker, kwargs = MAKERS[fmt] + path = str(tmp_path / f"{fmt}.h5") + maker(path) + expected = xd.open_dataarray(path, engine=fmt, **kwargs) + result = xd.open_dataarray(path, engine=fmt, vtype="tiles", **kwargs) + assert isinstance(result.data, TileArray) + assert result.data.engine["name"] == fmt + assert result.equals(expected) + sliced = result[3:9:2, 1:3] + assert isinstance(sliced.data, TileArray) + npt.assert_array_equal(sliced.values, expected.values[3:9:2, 1:3]) + + +def test_febus_block_crossing_reads(tmp_path): + """Row ranges spanning block boundaries fuse the right trimmed parts.""" + path = str(tmp_path / "febus.h5") + make_febus_file(path) + kwargs = {"overlaps": (1, 1), "offset": 0} + # the hdf5 vtype is the independent reference (tiles is the febus default) + expected = xd.open_dataarray(path, engine="febus", vtype="hdf5", **kwargs).values + result = xd.open_dataarray(path, engine="febus", vtype="tiles", **kwargs) + assert result.data.engine["block_size"] == 12 + assert result.data.engine["overlaps"] == [1, 1] + npt.assert_array_equal(result[8:12].values, expected[8:12]) + npt.assert_array_equal(result[::3].values, expected[::3]) + npt.assert_array_equal(result[25:].values, expected[25:]) + + +def test_prodml_transpose_param(tmp_path): + """The shared contract reads distance-major files time-major on request. + + The 0.2 engine never writes ``transpose`` (its manifests keep the + on-disk layout), but manifests written by the 0.3 line use it. + """ + path = str(tmp_path / "prodml_swapped.h5") + data = make_prodml_file(path, swapped=True) + manifest = TileArray.from_tiles( + path, data.T.shape, data.dtype, {"name": "prodml", "transpose": True} + ) + npt.assert_array_equal(np.asarray(manifest), data.T) + npt.assert_array_equal(np.asarray(manifest[2:7:2, 1:4]), data.T[2:7:2, 1:4]) + + +def test_tiles_view_roundtrip(tmp_path): + """A tile-backed view persists (spec params included) and reopens lazily.""" + path = str(tmp_path / "febus.h5") + make_febus_file(path) + da = xd.open_dataarray( + path, engine="febus", vtype="tiles", overlaps=(1, 1), offset=0 + ) + out = str(tmp_path / "view.nc") + da.to_netcdf(out) + reopened = xd.open_dataarray(out) + assert isinstance(reopened.data, TileArray) + assert reopened.data.engine == da.data.engine + assert reopened.equals(da) + + +def test_open_mfdataarray_fuses_tiles(tmp_path): + """Multi-file opening fuses at the manifest level and stays lazy.""" + paths, parts = [], [] + for k in range(2): + path = str(tmp_path / f"asn{k}.h5") + parts.append(make_asn_file(path, t0=k * 2.0)) + paths.append(path) + da = xd.open_mfdataarray(paths, engine="asn", vtype="tiles", parallel=False) + assert isinstance(da, xd.DataArray) + assert isinstance(da.data, TileArray) + assert da.data.ntiles == 2 + npt.assert_array_equal(da.values, np.concatenate(parts)) + + +def test_xdas_engine_tiles_vtype(tmp_path): + """Materialized native files reopen lazily as tile arrays on request.""" + da = xd.testing.dummy(shape=(10, 5), step=(1.0, 10.0), dtype=np.float32) + path = str(tmp_path / "da.nc") + da.to_netcdf(path) + result = xd.open_dataarray(path, engine="xdas", vtype="tiles") + assert isinstance(result.data, TileArray) + assert result.data.engine == {"name": "xdas", "dataset": "/__values__"} + assert result.equals(da) + + +def test_manifest_strings_stored_as_char_arrays(tmp_path): + """Manifest strings land on disk as fixed-width char arrays, not vlen.""" + da = xd.testing.dummy(shape=(10, 5), step=(1.0, 10.0), dtype=np.float32) + path = str(tmp_path / "da.nc") + da.to_netcdf(path) + tiled = xd.open_dataarray(path, engine="xdas", vtype="tiles") + out = str(tmp_path / "view.nc") + tiled.to_netcdf(out) + with h5py.File(out, "r") as file: + assert file["__tiles__/paths"].dtype == np.dtype("S1") + reopened = xd.open_dataarray(out) + assert isinstance(reopened.data, TileArray) + # and reopen as fixed-width bytes: no per-string decode, no heap + assert reopened.data.dataset["paths"].dtype.kind == "S" + npt.assert_array_equal(reopened.values, da.values) + + +def test_reopened_tile_file_stays_writable(tmp_path): + """Opening loads the manifest and closes the file: it accepts appends.""" + import h5netcdf + + da = xd.testing.dummy(shape=(10, 5), step=(1.0, 10.0), dtype=np.float32) + path = str(tmp_path / "da.nc") + da.to_netcdf(path) + tiled = xd.open_dataarray(path, engine="xdas", vtype="tiles") + out = str(tmp_path / "view.nc") + tiled.to_netcdf(out) + reopened = xd.open_dataarray(out) + with h5netcdf.File(out, "a") as file: + file.attrs["appended"] = 1 + npt.assert_array_equal(reopened.values, da.values) + + +def test_legacy_vlen_manifest_reopens(tmp_path): + """Manifests stored with variable-length strings (pre char-array) reopen.""" + da = xd.testing.dummy(shape=(10, 5), step=(1.0, 10.0), dtype=np.float32) + path = str(tmp_path / "da.nc") + da.to_netcdf(path) + tiled = xd.open_dataarray(path, engine="xdas", vtype="tiles") + out = str(tmp_path / "legacy.nc") + tiled.to_netcdf(out) + # rewrite the manifest group the way the old writer did: vlen strings + manifest = tiled.data.to_dataset() + for name in list(manifest.variables): + manifest[name].encoding.clear() + if manifest[name].dtype.kind == "S": + manifest[name] = manifest[name].astype(str) + with h5py.File(out, "a") as file: + del file["__tiles__"] + manifest.to_netcdf(out, mode="a", group="__tiles__", engine="h5netcdf") + with h5py.File(out, "r") as file: + assert h5py.check_string_dtype(file["__tiles__/paths"].dtype) is not None + reopened = xd.open_dataarray(out) + assert isinstance(reopened.data, TileArray) + assert reopened.data.dataset["paths"].dtype.kind == "S" + npt.assert_array_equal(reopened.values, da.values) + + +def test_tiles_datacollection_roundtrip(tmp_path): + """Collections of tile-backed arrays reopen as collections, not as errors. + + The manifest lives in a sibling group of the data array's variables, which + used to make the array look like a nested collection to the reader. + """ + das = [ + xd.testing.dummy(shape=(10, 5), step=(1.0, 10.0), dtype=np.float32) + for _ in range(2) + ] + paths = [] + for k, da in enumerate(das): + path = str(tmp_path / f"da{k}.nc") + da.to_netcdf(path) + paths.append(path) + tiled = [xd.open_dataarray(path, engine="xdas", vtype="tiles") for path in paths] + + sequence = xd.DataCollection(tiled, name="acquisition") + fname = str(tmp_path / "sequence.nc") + sequence.to_netcdf(fname, virtual=True) + result = xd.open_datacollection(fname) + assert len(result) == 2 + for expected, actual in zip(das, result): + assert isinstance(actual.data, TileArray) + npt.assert_array_equal(actual.values, expected.values) + + mapping = xd.DataCollection({"a": sequence, "b": sequence}, name="node") + fname = str(tmp_path / "mapping.nc") + mapping.to_netcdf(fname, virtual=True) + result = xd.open_datacollection(fname) + assert sorted(result) == ["a", "b"] + assert isinstance(result["a"][1].data, TileArray) + npt.assert_array_equal(result["b"][0].values, das[0].values) + + +def test_default_vtypes(): + """The backing each engine picks when the caller does not say. + + Pinned because it is what the I/O guide documents: formats that store + several blocks per file, or that HDF5 virtual datasets cannot serve at + all, default to tiles. + """ + from xdas.io import Engine + + expected = { + "apsensing": "hdf5", + "asn": "hdf5", + "febus": "tiles", + "miniseed": "tiles", + "prodml": "hdf5", + "silixa": "tiles", + "terra15": "hdf5", + "xdas": "hdf5", + } + assert {name: Engine[name]().vtype for name in expected} == expected diff --git a/tests/test_core.py b/tests/test_core.py index b2ea9ba3..83bc8864 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -25,6 +25,35 @@ def test_open_mfdataarray(self, tmp_path): with pytest.raises(FileNotFoundError): xd.open_mfdataarray(["not_existing_file.nc"]) + def test_open_mfdataarray_file_limit(self, tmp_path, monkeypatch): + from xdas.core import routines + + for idx, da in enumerate(wavelet_wavefronts(nchunk=3), start=1): + da.to_netcdf(tmp_path / f"{idx:03}.nc") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + with pytest.raises(NotImplementedError, match="the limit is 2"): + xd.open_mfdataarray(tmp_path / "00*.nc") + + def test_open_mfdataarray_no_file_limit_for_tiles(self, tmp_path, monkeypatch): + from xdas.core import routines + + for idx, da in enumerate(wavelet_wavefronts(nchunk=3), start=1): + da.to_netcdf(tmp_path / f"{idx:03}.nc") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + da = xd.open_mfdataarray(tmp_path / "00*.nc", engine="xdas", vtype="tiles") + assert da.shape == wavelet_wavefronts().shape + + def test_open_mfdataarray_file_limit_engine_instance(self, tmp_path, monkeypatch): + from xdas.core import routines + from xdas.io.xdas import XdasEngine + + for idx, da in enumerate(wavelet_wavefronts(nchunk=3), start=1): + da.to_netcdf(tmp_path / f"{idx:03}.nc") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + # the configured instance carries the vtype the limit is keyed on + with pytest.raises(NotImplementedError, match="the limit is 2"): + xd.open_mfdataarray(tmp_path / "00*.nc", engine=XdasEngine(vtype="hdf5")) + def test_open_mfdataarray_grouping(self, tmp_path): acqs = [ { diff --git a/tests/test_routines.py b/tests/test_routines.py index 8fafc7b8..8d46f041 100644 --- a/tests/test_routines.py +++ b/tests/test_routines.py @@ -468,22 +468,42 @@ def test_invalid_paths_type_raises(self): with pytest.raises(Exception, match="paths"): xd.open(123) - def test_callable_engine(self, tmp_path): + def test_engine_instance(self, tmp_path): + from xdas.io import Engine + da = xd.testing.dummy(shape=(10, 5)) path = str(tmp_path / "test.nc") da.to_netcdf(path) + result = xd.open_dataarray(path, engine=Engine["xdas"]()) + assert result.equals(da) - def my_engine(fname, **kwargs): - return xd.open_dataarray(fname) + def test_engine_instance_rejects_extra_config(self, tmp_path): + from xdas.io import Engine - result = xd.open_dataarray(path, engine=my_engine) - assert result.equals(da) + da = xd.testing.dummy(shape=(10, 5)) + path = str(tmp_path / "test.nc") + da.to_netcdf(path) + engine = Engine["xdas"]() + with pytest.raises(ValueError, match="configured engine instance"): + xd.open_dataarray(path, engine=engine, vtype="tiles") + with pytest.raises(ValueError, match="configured engine instance"): + xd.open_dataarray(path, engine=engine, group="somegroup") + + def test_unknown_engine_kwarg_raises(self, tmp_path): + da = xd.testing.dummy(shape=(10, 5)) + path = str(tmp_path / "test.nc") + da.to_netcdf(path) + with pytest.raises(TypeError, match="overlpas"): + xd.open_dataarray(path, engine="febus", overlpas=(1, 1)) + # auto-detection accepts no format-specific parameters + with pytest.raises(TypeError, match="overlaps"): + xd.open_dataarray(path, overlaps=(1, 1)) def test_invalid_engine_type_raises(self, tmp_path): da = xd.testing.dummy(shape=(10, 5)) path = str(tmp_path / "test.nc") da.to_netcdf(path) - with pytest.raises(ValueError, match="engine"): + with pytest.raises(TypeError, match="engine must be"): xd.open_dataarray(path, engine=42) @@ -724,3 +744,246 @@ def test_invalid_type_raises(self): with pytest.raises(TypeError, match="DataCollection"): _get_timeline_dataframe("not_valid") + + +class TestSortby: + def make_archive(self, tmp_path, vtype, pairs=((0, 2), (1, 3))): + """Save 4 time chunks and fuse them losslessly as two runs. + + `concat` sorts whatever it is given, so tile-level disorder is + built the way streamed scans produce it: runs that are internally + ordered but interleave each other. The default pairing yields the + tile order 0, 2, 1, 3. + """ + expected = xd.testing.dummy(dims=("time", "space"), shape=(20, 5)) + chunks = xd.split(expected, 4, "time") + parts = [] + for index, chunk in enumerate(chunks): + path = tmp_path / f"chunk_{index}.nc" + chunk.to_netcdf(path) + parts.append(xd.open_dataarray(path, engine="xdas", vtype=vtype)) + runs = [ + xd.concat([parts[i] for i in pair], "time", tolerance=False) + for pair in pairs + ] + return expected, xd.concat(runs, "time", tolerance=False) + + def test_sorts_tiles_lazily(self, tmp_path): + from xdas.virtual import TileArray + + expected, shuffled = self.make_archive(tmp_path, "tiles") + result = xd.sortby(shuffled, "time") + assert isinstance(result.data, TileArray) + assert result.equals(expected) + assert result["time"].equals(expected["time"]) + + def test_sorts_virtual_stack(self, tmp_path): + from xdas.virtual import VirtualStack + + expected, shuffled = self.make_archive(tmp_path, "hdf5") + result = xd.sortby(shuffled, "time") + assert isinstance(result.data, VirtualStack) + assert result.equals(expected) + + def test_already_sorted_fast_path(self, tmp_path): + expected, arranged = self.make_archive( + tmp_path, "tiles", pairs=((0, 1), (2, 3)) + ) + result = xd.sortby(arranged, "time") + assert result.equals(expected) + # a second sort is a no-op even though the coordinate is simplified + assert xd.sortby(result, "time").equals(expected) + + def test_tolerance_false_skips_simplification(self, tmp_path): + _, shuffled = self.make_archive(tmp_path, "tiles") + result = xd.sortby(shuffled, "time", tolerance=False) + # sorted but not simplified: one tie pair per chunk remains + assert len(result["time"].tie_indices) == 8 + assert bool(np.all(np.diff(result["time"].tie_values.astype("i8")) > 0)) + + def test_eager_data_raises(self): + da = xd.testing.dummy(dims=("time", "space"), shape=(10, 5)) + with pytest.raises(NotImplementedError, match="TileArray or a VirtualStack"): + xd.sortby(da, "time") + + def test_dense_coordinate_raises(self, tmp_path): + _, shuffled = self.make_archive(tmp_path, "tiles") + shuffled["time"] = shuffled["time"].values + with pytest.raises(NotImplementedError, match="interpolated"): + xd.sortby(shuffled, "time") + + def test_simplified_unsorted_raises(self, tmp_path): + # a coordinate whose ties span the first two tiles as one segment + # (the state a prior simplification leaves) while the tile order + # still needs fixing: the exact blockwise gather is impossible + from xdas.coordinates import InterpCoordinate + + _, misaligned = self.make_archive(tmp_path, "tiles") + coord = misaligned["time"] + misaligned["time"] = InterpCoordinate( + { + "tie_indices": np.array([0, 9, 10, 14, 15, 19]), + "tie_values": coord.tie_values[[0, 3, 4, 5, 6, 7]], + }, + "time", + ) + with pytest.raises(NotImplementedError, match="align"): + xd.sortby(misaligned, "time") + + +class TestStreamingCombine: + def save_shuffled(self, tmp_path, nchunk=6): + """Save chunks under names whose lexicographic order shuffles time.""" + expected = xd.testing.dummy(dims=("time", "space"), shape=(30, 5)) + names = ["e", "b", "f", "a", "d", "c"][:nchunk] + for chunk, name in zip(xd.split(expected, nchunk, "time"), names): + chunk.to_netcdf(tmp_path / f"{name}.nc") + return expected + + def test_matches_monolithic(self, tmp_path, monkeypatch): + from xdas.core import routines + + expected = self.save_shuffled(tmp_path) + mono = xd.open_mfdataarray( + tmp_path / "*.nc", engine="xdas", vtype="tiles", parallel=False + ) + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + streamed = xd.open_mfdataarray( + tmp_path / "*.nc", engine="xdas", vtype="tiles", parallel=False + ) + assert streamed.equals(expected) + assert streamed["time"].equals(mono["time"]) + np.testing.assert_array_equal(np.asarray(streamed.data), np.asarray(mono.data)) + + def test_non_consolidating_vtype_raises_instead_of_streaming( + self, tmp_path, monkeypatch + ): + from xdas.core import routines + + # the batch size is the ceiling, so a vtype that cannot consolidate + # never reaches the streaming path: it raises at the first batch + self.save_shuffled(tmp_path) + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + with pytest.raises(NotImplementedError, match="cannot be consolidated"): + xd.open_mfdataarray( + tmp_path / "*.nc", engine="xdas", vtype="hdf5", parallel=False + ) + + def test_warns_and_recovers_on_corrupted_file(self, tmp_path, monkeypatch): + from xdas.core import routines + + expected = self.save_shuffled(tmp_path) + with (tmp_path / "ba.nc").open("wb") as file: + file.write(b"corrupted") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + with pytest.warns(RuntimeWarning): + streamed = xd.open_mfdataarray( + tmp_path / "*.nc", engine="xdas", vtype="tiles", parallel=False + ) + assert streamed.equals(expected) + + def test_groups_interleaved_acquisitions_by_signature(self, tmp_path, monkeypatch): + from xdas.core import routines + + # acquisition A (5 channels) at t0 and t2, B (3 channels) at t1: + # signature grouping fuses A whole where the monolithic time-ordered + # walk would split it around B + wide = xd.testing.dummy(dims=("time", "space"), shape=(20, 5)) + chunks = xd.split(wide, 2, "time") + narrow = xd.testing.dummy(dims=("time", "space"), shape=(10, 3)) + narrow["time"] = narrow["time"] + ( + chunks[1]["time"][0].values - narrow["time"][0].values + ) + chunks[0].to_netcdf(tmp_path / "a.nc") + narrow.to_netcdf(tmp_path / "b.nc") + chunks[1].to_netcdf(tmp_path / "c.nc") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + streamed = xd.open_mfdataarray( + tmp_path / "*.nc", engine="xdas", vtype="tiles", parallel=False + ) + assert isinstance(streamed, xd.DataCollection) + assert len(streamed) == 2 + + def test_single_run_squeezes(self, tmp_path, monkeypatch): + from xdas.core import routines + + expected = self.save_shuffled(tmp_path) + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + collection = xd.open_mfdataarray( + tmp_path / "*.nc", + engine="xdas", + vtype="tiles", + parallel=False, + squeeze=False, + ) + assert isinstance(collection, xd.DataCollection) + assert len(collection) == 1 + assert collection[0].equals(expected) + + +class TestStreamingCombineFallbacks: + def test_dim_last_and_plain_name(self, tmp_path, monkeypatch): + from xdas.core import routines + + expected = xd.testing.dummy(dims=("time",), shape=(30,), step=0.01) + names = ["c", "a", "b"] + for chunk, name in zip(xd.split(expected, 3, "time"), names): + chunk.to_netcdf(tmp_path / f"{name}.nc") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + for dim in ("last", "time"): + result = xd.open_mfdataarray( + tmp_path / "*.nc", + dim=dim, + engine="xdas", + vtype="tiles", + parallel=False, + ) + assert result.equals(expected) + + def test_unsortable_group_falls_back_to_plain_concat(self, tmp_path, monkeypatch): + from xdas.core import routines + + # dense time coordinates: sortby cannot permute them, the group is + # concatenated with the tolerance directly (runs sorted by start) + expected = xd.testing.dummy( + dims=("time", "space"), shape=(30, 5), ctype="dense" + ) + for index, chunk in enumerate(xd.split(expected, 3, "time")): + chunk.to_netcdf(tmp_path / f"chunk_{index}.nc") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + result = xd.open_mfdataarray( + tmp_path / "*.nc", engine="xdas", vtype="tiles", parallel=False + ) + assert result.equals(expected) + + def test_no_dim_coordinate(self, tmp_path, monkeypatch): + from xdas.core import routines + + da = xd.DataArray( + np.arange(30.0 * 5).reshape(30, 5), + coords={"space": {"tie_indices": [0, 4], "tie_values": [0.0, 40.0]}}, + dims=("time", "space"), + ) + for index in range(3): + da[10 * index : 10 * (index + 1)].to_netcdf(tmp_path / f"chunk_{index}.nc") + monkeypatch.setattr(routines, "MAX_OPEN_FILES", 2) + result = xd.open_mfdataarray( + tmp_path / "*.nc", engine="xdas", vtype="tiles", parallel=False + ) + assert result.shape == (30, 5) + + +class TestSortbyMetadataFree: + def test_permutes_without_declared_sampling_interval(self, tmp_path): + from xdas.coordinates import InterpCoordinate + + helper = TestSortby() + expected, shuffled = helper.make_archive(tmp_path, "tiles") + coord = shuffled["time"] + shuffled["time"] = InterpCoordinate( + {"tie_indices": coord.tie_indices, "tie_values": coord.tie_values}, + "time", + ) + result = xd.sortby(shuffled, "time") + assert np.array_equal(result["time"].values, expected["time"].values) + np.testing.assert_array_equal(np.asarray(result.data), expected.values) diff --git a/tests/test_xdas.py b/tests/test_xdas.py index d3fcf89d..775815f8 100644 --- a/tests/test_xdas.py +++ b/tests/test_xdas.py @@ -2,8 +2,8 @@ import xdas as xd -# Release segment, optionally followed by a PEP 440 pre-release marker (e.g. 0.2.8rc0). -VERSION_PATTERN = re.compile(r"^\d+(\.\d+)*((a|b|rc)\d+)?$") +# Release segment, plus the optional PEP 440 pre/post/dev markers (e.g. 0.2.9.dev0). +VERSION_PATTERN = re.compile(r"^\d+(\.\d+)*((a|b|rc)\d+)?(\.post\d+)?(\.dev\d+)?$") def test_version(): diff --git a/tests/virtual/test_core.py b/tests/virtual/test_core.py new file mode 100644 index 00000000..a63b071a --- /dev/null +++ b/tests/virtual/test_core.py @@ -0,0 +1,72 @@ +import h5py +import numpy as np +import pytest + +from xdas.virtual import ( + TileArray, + VirtualArray, + VirtualBackend, + VirtualSource, + VirtualStack, +) + + +class TestVirtualBackend: + def test_lookup(self): + assert VirtualBackend["hdf5"] is VirtualArray + assert VirtualBackend["tiles"] is TileArray + + def test_unknown_vtype_raises_key_error(self): + with pytest.raises(KeyError, match="no virtual backend registered"): + VirtualBackend["netcdf"] + + def test_registry_holds_only_named_backends(self): + assert set(VirtualBackend._registry) == {"hdf5", "tiles"} + + def test_subclasses_inherit_vtype_without_reregistering(self): + assert VirtualSource.vtype == "hdf5" + assert VirtualStack.vtype == "hdf5" + assert VirtualBackend["hdf5"] is VirtualArray + + def test_consolidates(self): + assert TileArray.consolidates + assert not VirtualArray.consolidates + assert not VirtualBackend.consolidates + + def test_isinstance_covers_both_backends(self): + source = VirtualSource("path.h5", "data", (2, 3), "f8") + assert isinstance(source, VirtualBackend) + assert isinstance(VirtualStack([source]), VirtualBackend) + + def test_base_is_abstract(self): + with pytest.raises(TypeError, match="abstract"): + VirtualBackend() + + def test_finalize_save_defaults_to_nothing(self): + source = VirtualSource("path.h5", "data", (2, 3), np.dtype("f8")) + assert source.finalize_save("path.nc") is None + + def test_derived_properties_shared_by_both_backends(self): + source = VirtualSource("path.h5", "data", (2, 3), np.dtype("f8")) + assert source.ndim == 2 + assert source.size == 6 + assert source.nbytes == 48 + assert not source.empty + tiled = TileArray.from_tiles("path.h5", (2, 3), "f8", "xdas") + assert tiled.ndim == 2 + assert tiled.size == 6 + assert tiled.nbytes == 48 + assert not tiled.empty + + def test_from_variable_dispatches_to_each_backend(self, tmp_path): + data = np.arange(6.0).reshape(2, 3) + path = tmp_path / "source.h5" + with h5py.File(path, "w") as file: + file.create_dataset("data", data=data) + with h5py.File(path) as file: + source = VirtualBackend["hdf5"].from_variable(file["data"]) + tiled = VirtualBackend["tiles"].from_variable(file["data"]) + assert isinstance(source, VirtualSource) + assert isinstance(tiled, TileArray) + np.testing.assert_array_equal(np.asarray(source), data) + np.testing.assert_array_equal(np.asarray(tiled), data) diff --git a/tests/test_virtual.py b/tests/virtual/test_hdf5.py similarity index 97% rename from tests/test_virtual.py rename to tests/virtual/test_hdf5.py index 503e4c8d..2d9596a0 100644 --- a/tests/test_virtual.py +++ b/tests/virtual/test_hdf5.py @@ -3,7 +3,7 @@ import pytest import xdas as xd -from xdas.virtual import ( +from xdas.virtual.hdf5 import ( Selection, Selectors, SingleSelector, @@ -389,18 +389,13 @@ def test_get_indexer(self): class TestVirtualArrayAbstract: - def test_abstract_stubs(self): - va = VirtualArray() - with pytest.raises(NotImplementedError): - _ = va[0] - with pytest.raises(NotImplementedError): - va.__array__() - with pytest.raises(NotImplementedError): - _ = va.shape - with pytest.raises(NotImplementedError): - _ = va.dtype + def test_cannot_instantiate(self): + with pytest.raises(TypeError, match="abstract"): + VirtualArray() + + def test_to_dataset_stub(self): with pytest.raises(NotImplementedError): - va.to_dataset(None, None) + VirtualArray.to_dataset(None, None, None) class TestVirtualStackExtra: diff --git a/tests/virtual/test_tiles.py b/tests/virtual/test_tiles.py new file mode 100644 index 00000000..12de5034 --- /dev/null +++ b/tests/virtual/test_tiles.py @@ -0,0 +1,2171 @@ +"""The tile-backed virtual array and its integration in the DataArray and native format.""" + +import math +import os + +import dask.array as da_ +import h5py +import numpy as np +import numpy.testing as npt +import pytest +import xarray as xr + +import xdas as xd +from xdas.io import Engine +from xdas.virtual.tiles import TileArray + +NX = 5 + +DIMS = ("time", "distance") + +ENGINE = {"name": "h5py", "dataset": "data"} + + +class H5pyEngine(Engine, name="h5py"): + """Read any HDF5 dataset — the engine of the synthetic test files. + + The format engines each read their own layout; test files belong to + no format, so they are described by this generic load-only engine + (its opening half stays abstract). The selection always has one + slice per source axis, in source order, whatever the virtual + arrangement. + """ + + @staticmethod + def load_tile(path, selection, *, dataset): + with h5py.File(path, "r") as file: + return file[dataset][selection] + + +@pytest.fixture +def stack(tmp_path): + """Three gzip-compressed HDF5 files with junk edge rows to trim. + + Emulates overlap trimming: each file carries one junk row at its start + and end that the tile's start row (plus the row ``size``) cuts out. + Returns the manifest and the expected stacked values. + """ + paths = [] + sizes = [] + parts = [] + row = 0 + for k, raw_nt in enumerate([12, 9, 14]): + path = str(tmp_path / f"src{k}.h5") + useful = raw_nt - 2 + data = np.full((raw_nt, NX), -999.0) + data[1:-1] = (row + np.arange(useful))[:, None] + np.arange(NX) / 10 + with h5py.File(path, "w") as file: + file.create_dataset("data", data=data, chunks=(4, NX), compression="gzip") + paths.append(path) + sizes.append(useful) + parts.append(data[1:-1]) + row += useful + manifest = TileArray.from_tiles( + paths, (sizes, NX), "float64", ENGINE, attrs={"units": "strain"} + ) + # per-tile source origins are view state: assigned through the manifest + manifest = TileArray( + manifest.dataset.assign(starts_0=("tile_0", np.array([1, 1, 1]))), + manifest.dtype, + manifest.engine, + ) + return manifest, np.concatenate(parts) + + +@pytest.fixture +def windowed(tmp_path): + """Three files whose rows contribute blob-local windows via ``starts_0``. + + Each file holds junk rows around the useful window; the manifest + exposes blob rows ``[start, start + size)``. The middle file has a + zero start (window at the top of the blob). + """ + paths, sizes, starts, parts = [], [], [], [] + row = 0 + for k, raw_nt in enumerate([12, 9, 14]): + path = str(tmp_path / f"win{k}.h5") + useful = raw_nt - 4 + first = 0 if k == 1 else 2 + data = np.full((raw_nt, NX), -999.0) + good = (row + np.arange(useful))[:, None] + np.arange(NX) / 10 + data[first : first + useful] = good + with h5py.File(path, "w") as file: + file.create_dataset("data", data=data) + paths.append(path) + sizes.append(useful) + starts.append(first) + parts.append(good) + row += useful + manifest = TileArray.from_tiles( + paths, (sizes, NX), "float64", {"name": "h5py", "dataset": "data"} + ) + manifest = TileArray( + manifest.dataset.assign(starts_0=("tile_0", np.array(starts))), + manifest.dtype, + manifest.engine, + ) + return manifest, np.concatenate(parts) + + +@pytest.fixture +def engine_calls(monkeypatch): + """Record the path of every h5py engine read, delegating to the real one.""" + calls = [] + original = Engine["h5py"].load_tile + + def counting(path, selection, **params): + calls.append(path) + return original(path, selection, **params) + + monkeypatch.setattr(Engine["h5py"], "load_tile", counting) + return calls + + +def _tile_file(path, data, **kwargs): + """Write *data* to an HDF5 file at *path*.""" + with h5py.File(path, "w") as file: + file.create_dataset("data", data=data, **kwargs) + + +def _with_starts(manifest, *starts): + """Rebuild *manifest* with per-axis tile origins inside their sources. + + ``starts_k`` is view state the :meth:`TileArray.from_tiles` encoder + does not take: windowed manifests assign it through the dataset and + the canonical constructor. + """ + assign = { + f"starts_{k}": (f"tile_{k}", np.asarray(entry, dtype=np.int64)) + for k, entry in enumerate(starts) + if entry is not None + } + return TileArray(manifest.dataset.assign(assign), manifest.dtype, manifest.engine) + + +def _random_key(rng, shape, max_step=1): + """A random non-empty positive-step slice per axis.""" + key = [] + for extent in shape: + a = int(rng.integers(0, extent)) + b = int(rng.integers(a + 1, extent + 1)) + s = int(rng.integers(1, max_step + 1)) + key.append(slice(a, b, s)) + return tuple(key) + + +def _random_grid(tmp_path, rng, ndim): + """A random rectilinear grid of junk-padded files, one per tile. + + Per-axis margins model source cropping (axis-separable, as the grid + requires): every tile at index ``i`` along axis ``k`` starts at + ``margins[k][i]`` inside its own file. + """ + counts = tuple(int(rng.integers(1, 4)) for _ in range(ndim)) + sizes = [rng.integers(2, 6, count).astype(np.int64) for count in counts] + margins = [rng.integers(0, 3, count).astype(np.int64) for count in counts] + shape = tuple(int(entry.sum()) for entry in sizes) + edges = [np.concatenate(([0], np.cumsum(entry))) for entry in sizes] + reference = np.empty(shape) + paths = np.empty(counts, dtype=object) + for number, index in enumerate(np.ndindex(counts)): + extents = tuple(int(sizes[k][i]) for k, i in enumerate(index)) + raw = tuple( + int(margins[k][i]) + extent + int(rng.integers(0, 2)) + for (k, i), extent in zip(enumerate(index), extents) + ) + data = np.full(raw, -1.0) + block = 1000.0 * number + np.arange(math.prod(extents)).reshape(extents) + inner = tuple( + slice(int(margins[k][i]), int(margins[k][i]) + extent) + for (k, i), extent in zip(enumerate(index), extents) + ) + data[inner] = block + placed = tuple( + slice(int(edges[k][i]), int(edges[k][i + 1])) for k, i in enumerate(index) + ) + reference[placed] = block + path = str(tmp_path / f"grid{number}.h5") + _tile_file(path, data) + paths[index] = path + manifest = TileArray.from_tiles( + paths, sizes, "float64", {"name": "h5py", "dataset": "data"} + ) + return _with_starts(manifest, *margins), reference + + +class TestManifest: + def test_shape_and_geometry(self, stack): + manifest, reference = stack + assert manifest.shape == reference.shape + assert manifest.ntiles == 3 + npt.assert_array_equal(manifest._edges[0], [0, 10, 17, 29]) + + def test_reads_across_sources(self, stack): + manifest, reference = stack + npt.assert_array_equal(np.asarray(manifest), reference) + npt.assert_array_equal(np.asarray(manifest[9:13]), reference[9:13]) + npt.assert_array_equal(np.asarray(manifest[3:5]), reference[3:5]) + + def test_dataset_model(self, stack, tmp_path): + manifest, _ = stack + dataset = manifest.dataset + assert tuple(dataset["sizes_0"].dims) == ("tile_0",) + assert tuple(dataset["sizes_1"].dims) == ("tile_1",) + # per-file paths vary along tile_0 only: the trailing axis folds + assert tuple(dataset["paths"].dims) == ("tile_0",) + # the common directory splits off: 0-d root, root-relative paths + assert dataset["root"].ndim == 0 + assert os.fsdecode(dataset["root"].values[()]) == str(tmp_path) + # strings are held as fixed-width bytes, not str objects + assert dataset["paths"].dtype.kind == "S" + assert dataset["paths"].values.tolist() == [b"src0.h5", b"src1.h5", b"src2.h5"] + npt.assert_array_equal(dataset["starts_0"].values, [1, 1, 1]) + # all-default geometry columns are not stored + assert "starts_1" not in dataset and "steps_0" not in dataset + + def test_param_folding(self, stack): + manifest, _ = stack + path = manifest._full_paths().item(0) + uniform = TileArray.from_tiles( + path, ([10, 10, 10], NX), "float64", ENGINE, record=0, nbytes=80 + ) + # one path everywhere: 0-d; uniform per-tile params: 0-d + assert uniform.dataset["paths"].ndim == 0 + assert uniform.dataset["record"].ndim == 0 + assert uniform.shape == (30, NX) + varying = TileArray.from_tiles( + path, ([10, 10], NX), "float64", ENGINE, record=[[0], [80]] + ) + assert tuple(varying.dataset["record"].dims) == ("tile_0",) + + def test_validation(self, stack): + manifest, _ = stack + with pytest.raises(ValueError, match="at least one axis"): + TileArray.from_tiles("a", (), "f8", ENGINE) + with pytest.raises(ValueError, match="little-endian"): + TileArray.from_tiles("a", (5, NX), ">f8", ENGINE) + with pytest.raises(ValueError, match="strictly positive"): + TileArray.from_tiles("a", (0, NX), "f8", ENGINE) + with pytest.raises(ValueError, match="does not match the grid"): + TileArray.from_tiles( + np.array(["a", "b"], dtype=object), ([1, 2, 3], NX), "f8", ENGINE + ) + with pytest.raises(ValueError, match="reserved"): + TileArray.from_tiles("a", (5, NX), "f8", ENGINE, sizes_0=[5]) + with pytest.raises(ValueError, match="reserved"): + TileArray.from_tiles("a", (5, NX), "f8", ENGINE, starts_0=[0]) + with pytest.raises(ValueError, match="reserved"): + TileArray.from_tiles("a", (5, NX), "f8", ENGINE, root=["r"]) + bad_root = manifest.dataset.copy() + bad_root["root"] = (("tile_0",), np.array(["a", "b", "c"], dtype=object)) + with pytest.raises(ValueError, match="0-d"): + TileArray(bad_root, manifest.dtype, manifest.engine) + dataset = manifest.dataset.copy() + with pytest.raises(ValueError, match="`sizes_0`"): + TileArray( + dataset.drop_vars(["sizes_0", "sizes_1"]), + manifest.dtype, + manifest.engine, + ) + with pytest.raises(ValueError, match="`paths`"): + TileArray(dataset.drop_vars("paths"), manifest.dtype, manifest.engine) + bad_starts = dataset.assign(starts_0=("tile_0", np.array([-1, 0, 0]))) + with pytest.raises(ValueError, match="non-negative"): + TileArray(bad_starts, manifest.dtype, manifest.engine) + + def test_extra_variables_are_params(self, stack): + """Any non-geometry manifest variable is a per-tile engine parameter.""" + manifest, _ = stack + arr = TileArray( + manifest.dataset.assign(record=(("tile_0",), np.arange(3))), + manifest.dtype, + manifest.engine, + ) + assert arr._params == ("record",) + + def test_string_params_decode_to_str(self, tmp_path): + """Per-tile string parameters store as bytes but reach the engine as str.""" + paths, parts = [], [] + for k in range(2): + path = str(tmp_path / f"named{k}.h5") + data = 100.0 * k + np.arange(3.0 * NX).reshape(3, NX) + with h5py.File(path, "w") as file: + file.create_dataset(f"data{k}", data=data) + paths.append(path) + parts.append(data) + manifest = TileArray.from_tiles( + paths, ([3, 3], NX), "float64", "h5py", dataset=["data0", "data1"] + ) + assert manifest.dataset["dataset"].dtype.kind == "S" + npt.assert_array_equal(np.asarray(manifest), np.concatenate(parts)) + + def test_engine_validation(self): + with pytest.raises(KeyError, match="no engine registered"): + TileArray.from_tiles("a", (5, NX), "f8", {"name": "bogus"}) + with pytest.raises(ValueError, match="`name` key"): + TileArray.from_tiles("a", (5, NX), "f8", {"dataset": "data"}) + with pytest.raises(ValueError, match="`name` key"): + TileArray.from_tiles("a", (5, NX), "f8", None) + + def test_engine_string_shorthand(self): + arr = TileArray.from_tiles("a", (5, NX), "f8", "h5py") + assert arr.engine == {"name": "h5py"} + + def test_engine_registration(self): + class DummyEngine(Engine, name="dummy"): + @staticmethod + def load_tile(path, selection): + return np.zeros((1, 1)) + + try: + assert Engine["dummy"] is DummyEngine + finally: + del Engine._registry["dummy"] + + def test_engine_without_tile_loader(self): + # a registered engine that predates the tiles machinery resolves + # but fails loudly when a manifest asks it to decode + class NoTilesEngine(Engine, name="notiles"): + pass + + try: + arr = TileArray.from_tiles("a", (5, NX), "f8", {"name": "notiles"}) + with pytest.raises(NotImplementedError): + np.asarray(arr) + finally: + del Engine._registry["notiles"] + + def test_repr(self, stack): + manifest, _ = stack + assert repr(manifest) == "TileArray[h5py] 1kB (float64) 3 tiles" + assert manifest._repr_inline_(40) == "TileArray[h5py] (3 tiles)" + assert manifest._repr_inline_(10) == "TileArray" + + def test_repr_of_a_single_tile(self, tmp_path): + """One tile reads as one tile, and the volume scales with the array.""" + path = str(tmp_path / "one.h5") + _tile_file(path, np.zeros((250, NX))) + arr = TileArray.from_tiles(path, (250, NX), "float64", ENGINE) + assert repr(arr) == "TileArray[h5py] 10kB (float64) 1 tile" + assert arr._repr_inline_(40) == "TileArray[h5py] (1 tile)" + + def test_relative_paths_are_anchored(self, tmp_path, monkeypatch): + """Relative paths absolutize at construction and survive a chdir.""" + data = np.arange(4.0 * NX).reshape(4, NX) + _tile_file(tmp_path / "rel.h5", data) + monkeypatch.chdir(tmp_path) + manifest = TileArray.from_tiles("rel.h5", (4, NX), "f8", ENGINE) + assert manifest.root == str(tmp_path) + assert os.path.isabs(manifest._full_paths().item(0)) + monkeypatch.chdir(tmp_path.parent) + npt.assert_array_equal(np.asarray(manifest), data) + + def test_attrs(self, stack): + manifest, _ = stack + assert manifest.attrs == {"units": "strain"} + + +class TestSourcePaths: + """Paths are stored split: a common 0-d root and root-relative values.""" + + def make(self, path): + # the file's first row is skipped: sliced away, as views are made + return TileArray.from_tiles([str(path)], ([5], NX), "