diff --git a/.gitea/workflows/ci.yml b/.gitea/workflows/ci.yml index cfd9624..06b42e3 100644 --- a/.gitea/workflows/ci.yml +++ b/.gitea/workflows/ci.yml @@ -38,7 +38,9 @@ jobs: # (clawhdf5-tools) interop tests compare against. apt-get install -y --no-install-recommends python3 python3-venv cmake hdf5-tools python3 -m venv /opt/interop - /opt/interop/bin/pip install --no-cache-dir h5py numpy netCDF4 xarray hdf5plugin + # maturin + pytest: ci-test.sh builds the Python package + # (crates/clawhdf5-py) and runs its tests against h5py. + /opt/interop/bin/pip install --no-cache-dir h5py numpy netCDF4 xarray hdf5plugin maturin pytest echo "/opt/interop/bin" >> "$GITHUB_PATH" - name: Show interop library versions # h5dump's version too: the h5rs dump test requires its exact output diff --git a/.gitignore b/.gitignore index 029ea7a..cf90c2c 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,5 @@ benchmarks/longmemeval/*.json # Local model weights (MiniLM etc.) — large, not committed weights/ .venv +__pycache__/ +.pytest_cache/ diff --git a/CHANGELOG.md b/CHANGELOG.md index ce4b7f2..29804fc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -230,6 +230,119 @@ are not byte-identical to earlier versions). Regression test: `crates/clawhdf5/tests/writer_groups_interop.rs`. +### Python bindings (2026-09-26) +- **Panic: selections of v4 implicit-index chunked datasets** (pre-existing, + facade `Dataset::read_selection`, Rust callers too). A hyperslab whose + bounding box covered more than half of a chunked dataset with the implicit + index (`libver='latest'`, early allocation, no filters) panicked with + "index out of bounds" in `generate_implicit_chunks`: the fallback in + `data_read::read_raw_data_selection` passed the layout's chunk dimensions, + element-size dimension included, and then decoded the whole dataset + anyway. That arm now decodes and extracts directly, for every chunk index. + `crates/clawhdf5/tests/v4_chunk_index_selection.rs` reads small and large + hyperslabs of all five v4 indexes (single chunk, implicit, fixed array, + extensible array, B-tree v2) and compares them with h5py; it panicked + before. The Python bindings made this easy to reach (`ds[0:3]` on + libhdf5's `h5fc_ext*.h5` test files). +- **`pip install` / `maturin develop` now gives `import clawhdf5`.** The + distribution in `crates/clawhdf5-py/pyproject.toml` was still called + `rustyhdf5` while the extension module was `clawhdf5`, and the package's + tests imported `rustyhdf5`, so they failed at collection. Distribution, + module and tests now all say `clawhdf5`, and the module has + `__version__`. +- **h5py-style reads that read the selection, not the dataset.** `ds[...]` + used to read the whole dataset and slice it in numpy, and knew six + dtypes. Now + integers (negative from the end), slices with positive steps, `...`, one + increasing list of integers per key and compound field names map onto the + facade's hyperslab selection (a list is read one group of neighbouring + chunks at a time and picked from in memory), with h5py's results (numpy scalar for an all-integer + key, 0-d array for `scalar[...]`) and h5py's errors for everything else + (negative steps, `None`, boolean masks, out-of-range indices). + `Dataset.dtype` is the numpy dtype h5py reports, for every integer and + IEEE float width (incl. `float16`) in either byte order, `bool`, enums + (base integer with `metadata['enum']`), complex (`r`/`i` compounds), + fixed strings (`S`), variable-length strings (`object` of `bytes`, as + h5py), variable-length sequences (`object` of arrays), opaque (`V`), + HDF5 array types and compounds (numpy structured, offsets and padding + kept, nested). The bytes the library returns become the numpy array's + buffer without a copy. Types the mapping cannot describe exactly + (references, bitfields, time, non-IEEE floats, integers with padding + bits, variable-length members inside compounds) raise `TypeError` rather + than return guessed data. Attributes come back as h5py returns them + (numpy scalars and arrays with the stored dtype, `str` for + variable-length strings, `numpy.bytes_` for fixed ones — **a change**: + string attributes written by this package are fixed-length and used to + come back as `str` — and `clawhdf5.Empty` for a null dataspace, which + datasets return too). `Group`/`File` gain `get`, `values`, `items`, + iteration, `len`, `name`, absolute and relative paths (`g['/a/b']`, + `g['c/d']`, `f['/']`); `Dataset` gains `ndim`, `size`, `maxshape`, + `name`, `len()` and `numpy.asarray(ds)`. File access and decoding run + with the GIL released, so Python threads read in parallel. + `crates/clawhdf5-py/tests/test_read_vs_h5py.py` compares every read with + h5py 3.16 (HDF5 2.0) on a file h5py writes. One difference is h5py's: + it returns variable-length sequences of big-endian floats unswapped; this + package returns the stored values. +- **A panic in the library is an ordinary Python exception.** PyO3 turns a + Rust panic into `PanicException`, a `BaseException` that `except + Exception` does not catch. Every call from the bindings into the library + is now guarded and a panic becomes `clawhdf5.InternalError` (a + `RuntimeError`) naming the object; with the implicit-index panic above + restored, `ds[0:30]` raises it. +- **Wrong data: uninitialised padding in compound results of index lists.** + `ds[[0, 3, 6]]` joined one read per run with `np.concatenate`, which + copies structured dtypes field by field into an `np.empty` result, so the + padding bytes held whatever was in memory (pointers were seen) and leaked + through `tobytes()`, hashes and write-backs. The runs' bytes are now joined + in Rust, whole elements at a time, so the result carries the bytes read + from the file (h5py's, zero for files it wrote) and stays zero-copy. + The h5py comparisons now also compare every byte of structured values + (`test_compound_padding_bytes_match_h5py` and `assert_same`). +- **Index lists no longer decode the same chunks once per run.** A list + index was one uncached hyperslab read per run of consecutive indices, so + on a chunked, compressed dataset every run decoded its chunk again: + `d[list(range(0, 200000, 40))]` over 20 gzip chunks took 8 s (h5py: + 0.014 s). The list is now read in groups — for a chunked dataset a group + ends only where a whole chunk holds no selected index, so each chunk is + decoded once; otherwise at a gap of more than 64 KiB — and the selected + rows are picked from each group in Rust. The same read now takes 3.8 ms + (h5py 4.1 ms; release build on tank, best of 5). + `test_a_long_index_list_decodes_each_chunk_once` compares 1-D, 2-D and + contiguous cases with h5py under a 2 s bound (5.8 s before, debug build). +- **Groups and datasets remember where they are.** Every `ds[...]`, and + every `g[k]`, resolved its path from the root again (two or three times + per open), and in a large group each resolution scans the group's links, + so visiting a group was quadratic: 4000 scalar datasets in one group took + 39 s (`libver='earliest'`) and 131 s (`'latest'`) to list, read and + re-read in `test_big_groups_are_not_quadratic`; now 0.3 s each (debug + build). A `Dataset` keeps its object's address, and a `Group` (and the + file's root) its address and, once listed, its link table. New facade + API: `File::dataset_at(address)` opens a dataset without resolving a + path. libhdf5's `h5stat_newgrat.h5` (35001 members in the root): listing + takes 0.03 s and 2000 opens 1 ms (h5py: 0.022 s). +- **`ds[np.array(1)]` is an integer index**, as in h5py; a 0-d integer + array went down the index-list path and raised a confusing `TypeError`. + The h5py comparison keys now include 0-d arrays on every axis. +- **Tests that would notice a held GIL, and our extra errors.** + `test_reads_release_the_gil` times a Python thread spinning while another + reads: with the read made to hold the GIL it stalls for the whole read + (0.062 s of a 0.064 s read) and the test fails; released, its longest + stall is about 3 ms. (The existing threads test only checked values.) + `test_errors_match_h5py` now also requires that every key h5py reads + reads here too, with the same result, and covers more keys (0-d arrays, + repeated and empty lists, `()`, `...`). +- **Docs say when a selection reads more than itself.** The README and + the package README said `ds[...]` reads only the selected elements, + without condition. The library decodes the whole dataset when the + selection's bounding box covers more than half of it, and for compact, + virtual, unwritten and non-default-fill chunked datasets; the READMEs, + the facade's `read_selection` docs and `docs/known-issues.md` now say so. +- **CI builds and tests the Python package.** It was excluded from CI. + `scripts/ci-test.sh` now lints `clawhdf5-py`, builds the wheel with + maturin, unpacks it under `target/` and runs the pytest suite; skipped + without maturin/pytest in `$CLAWHDF5_PYTHON`, a failure then under + `CLAWHDF5_REQUIRE_INTEROP=1`. The CI interop venv installs both. + ### Plugin filters (2026-09-26) - **LZF, bitshuffle, bzip2 and Blosc read and write, in pure Rust.** Files written by h5py with `compression="lzf"`, or with hdf5plugin's diff --git a/README.md b/README.md index 20294ec..4ae4d1c 100644 --- a/README.md +++ b/README.md @@ -73,8 +73,9 @@ breaking change, are in [CHANGELOG.md](CHANGELOG.md). - Default fusion weights are now the measured 0.4 / 0.6 (v2.5.0). Re-ranking had been discarding the retrieval score, costing the Markdown backend 40.6pp of Hit@1; fixed in v2.6.0. -- Selection reads decode only the chunks they touch (a 64×64 window: 105 ms to - 0.39 ms), and full reads are 1.2–1.9× faster (v2.5.0). +- Selection reads whose bounding box covers at most half the dataset decode + only the chunks they touch (a 64×64 window: 105 ms to 0.39 ms), and full + reads are 1.2–1.9× faster (v2.5.0). **Memory** - A loaded store holds ~30% less (embeddings stored once, v2.6.0), and the @@ -432,6 +433,54 @@ b.write("groups.h5")?; A group holds at most 65 535 links; more is an error, as is a link over 65 515 bytes (a very long soft-link target) in a group of more than 8 links. +### Python + +`crates/clawhdf5-py` is a Python package (PyO3 + numpy) that reads HDF5 with +an h5py-shaped API and no libhdf5. It is not on PyPI; build it with +[maturin](https://www.maturin.rs) into a virtualenv: + +```bash +python -m venv .venv && . .venv/bin/activate +pip install maturin numpy +maturin develop --release -m crates/clawhdf5-py/Cargo.toml +python -c "import clawhdf5; print(clawhdf5.__version__)" +``` + +```python +import numpy as np +import clawhdf5 + +with clawhdf5.File("data.h5", "r") as f: + print(list(f.keys())) # sorted member names, like h5py + ds = f["group/temperatures"] # relative or absolute ("/group/...") paths + print(ds.shape, ds.dtype) # dtype is the numpy dtype h5py reports + block = ds[100:200, ::4] # a small selection reads only its chunks + row = ds[-1] # integers drop the axis + picked = ds[[1, 5, 9], :] # one increasing index list per key + units = ds.attrs["units"] # attributes come back as h5py returns them + everything = np.asarray(ds) + + records = f["table"] # compound -> numpy structured array + ids = records["id"] # one field +``` + +Reads cover integers and IEEE floats of every width in either byte order, +`bool`, enums, complex, fixed and variable-length strings, variable-length +sequences, opaque, HDF5 array types and compounds; other types (references, +bitfields, ...) raise `TypeError` instead of returning guessed data. Keys +follow h5py (negative steps, `None` and boolean masks are refused). The +read itself runs with the GIL released, so Python threads read in parallel. +A selection whose bounding box covers at most half the dataset decodes only +the chunks (or contiguous rows) that box overlaps; a larger one — including +a strided slice across the whole dataset — decodes the whole dataset, as +do datasets that are compact, virtual, unwritten, or chunked with a +non-default fill value (`docs/known-issues.md`). An index list is read one +group of neighbouring chunks at a time. +Writing (`File(path, "w")`, `create_dataset`, `create_group`, `attrs[...] =`) +covers `float64`, `float32`, `int64`, `int32` and `uint8` arrays. The tests +in `crates/clawhdf5-py/tests` compare every read with h5py; run them with +`pip install pytest h5py && pytest crates/clawhdf5-py/tests`. + ### Agent Memory ```rust diff --git a/crates/clawhdf5-format/src/chunked_read.rs b/crates/clawhdf5-format/src/chunked_read.rs index ffc1379..70176e2 100644 --- a/crates/clawhdf5-format/src/chunked_read.rs +++ b/crates/clawhdf5-format/src/chunked_read.rs @@ -513,6 +513,10 @@ fn collect_chunk_info_inner( /// /// Chunks are stored contiguously starting at `base_address`. No stored index; /// addresses are computed from the chunk position. +/// +/// `chunk_dimensions` are the spatial chunk dimensions, one per entry of +/// `dataset_dims` — not the layout message's list, which carries the element +/// size as an extra last dimension. pub fn generate_implicit_chunks( base_address: u64, dataset_dims: &[u64], diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 633ecd7..8aec190 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -288,9 +288,11 @@ pub fn read_raw_data_indexed( /// Read raw bytes for only the selected elements of a dataset. /// -/// For chunked layouts, only chunks that intersect the selection are read -/// and decompressed. For compact/contiguous layouts, the full data is read -/// and then the selection is extracted. +/// When the selection's bounding box covers at most half the dataset, only +/// that box is materialised — the overlapping rows of a contiguous dataset, +/// the overlapping chunks of a chunked one, whatever its chunk index (see +/// [`crate::partial_read`]). Otherwise, and for compact and virtual +/// layouts, the whole dataset is decoded and the selection extracted. #[allow(clippy::too_many_arguments)] pub fn read_raw_data_selection( file_data: &[u8], @@ -358,85 +360,17 @@ pub fn read_raw_data_selection( } DataLayout::Chunked { chunk_dimensions, - btree_address, version, - chunk_index_type, .. } => { + // `partial_read` declined (a bounding box covering most of the + // dataset, or a selection it doesn't box), so decode every chunk + // and pick the selection out, whatever the chunk index. This arm + // used to enumerate the chunks first — passing the layout's + // chunk dimensions, element-size dimension included, to the + // implicit-index generator, which then indexed past the rank and + // panicked — only to decode the full dataset anyway. crate::chunked_read::chunk_geometry(chunk_dimensions, *version, dataspace, elem_size)?; - // For chunked data, only read chunks that intersect the selection - let chunk_dims: Vec = chunk_dimensions.iter().map(|&d| d as u64).collect(); - let rank = dims.len(); - - // Collect chunk info from B-tree - let chunks = if *version == 4 { - match chunk_index_type { - Some(2) => { - // Implicit index - crate::chunked_read::generate_implicit_chunks( - btree_address.unwrap_or(0), - dims, - chunk_dimensions, - elem_size as u32, - ) - } - _ => { - if let Some(_addr) = btree_address { - // Use extensible array or fixed array - // Fall back to full read for complex v4 index types - let full_data = read_raw_data_full( - file_data, - layout, - dataspace, - datatype, - pipeline, - offset_size, - length_size, - )?; - return extract_selection_from_buffer( - &full_data, dims, elem_size, selection, - ); - } else { - return Ok(Vec::new()); - } - } - } - } else { - // v3: B-tree v1 - if let Some(addr) = btree_address { - crate::chunked_read::collect_chunk_info_checked( - file_data, - *addr, - chunk_dimensions, - offset_size, - length_size, - )? - } else { - return Ok(Vec::new()); - } - }; - - // Filter chunks to only those that intersect the selection - let intersecting: Vec<_> = chunks - .iter() - .filter(|ci| { - let offsets: Vec = ci.offsets.iter().take(rank).copied().collect(); - selection.intersects_chunk(&offsets, &chunk_dims[..rank]) - }) - .collect(); - - if intersecting.is_empty() { - return Ok(Vec::new()); - } - - // Decompress only the intersecting chunks - let _chunk_total_bytes: usize = - chunk_dims.iter().map(|&d| d as usize).product::() * elem_size; - let _element_size_u32 = elem_size as u32; - - // First, assemble only the intersecting chunks into a partial buffer, - // then extract the selection. For simplicity, we assemble into a full - // dataset buffer and extract (same as contiguous path). let full_data = read_raw_data_full( file_data, layout, diff --git a/crates/clawhdf5-py/Cargo.toml b/crates/clawhdf5-py/Cargo.toml index 0d54a4f..3513a62 100644 --- a/crates/clawhdf5-py/Cargo.toml +++ b/crates/clawhdf5-py/Cargo.toml @@ -3,7 +3,7 @@ name = "clawhdf5-py" version = "2.7.0" edition = "2024" rust-version.workspace = true -description = "Python bindings for rustyhdf5 — a pure-Rust HDF5 library" +description = "Python bindings for clawhdf5 — a pure-Rust HDF5 library" license = "MIT" repository = "https://git.redclaw.dev/quantumclaw/clawhdf5" readme = "README.md" diff --git a/crates/clawhdf5-py/README.md b/crates/clawhdf5-py/README.md index 6554f32..95cae12 100644 --- a/crates/clawhdf5-py/README.md +++ b/crates/clawhdf5-py/README.md @@ -3,23 +3,81 @@ [![crates.io](https://img.shields.io/crates/v/clawhdf5-py.svg)](https://crates.io/crates/clawhdf5-py) [![docs.rs](https://docs.rs/clawhdf5-py/badge.svg)](https://docs.rs/clawhdf5-py) -Python bindings for clawhdf5 — a pure-Rust HDF5 library. +Python bindings for clawhdf5 — a pure-Rust HDF5 library. The package is +`clawhdf5` (`import clawhdf5`); it needs numpy and no libhdf5. -## Features +## Install -- h5py-compatible API (`File`, `Group`, `Dataset`) -- NumPy array integration -- Read and write HDF5 files from Python with no C dependencies +Not on PyPI yet. Build it into a virtualenv with [maturin](https://www.maturin.rs): -## Usage +```bash +pip install maturin numpy +cd crates/clawhdf5-py +maturin develop --release +python -c "import clawhdf5; print(clawhdf5.__version__)" +``` + +## Reading + +The read API follows h5py: ```python +import numpy as np import clawhdf5 -with clawhdf5.File('data.h5', 'r') as f: - data = f['/dataset'][:] +with clawhdf5.File("data.h5", "r") as f: + f.keys(), f["group"].items(), "group/data" in f + ds = f["group/data"] # or f["/group/data"], f["group"]["data"] + ds.shape, ds.dtype, ds.attrs["units"] + ds[10:20, ::2] # a small selection reads only its chunks + ds[-1], ds[..., 0], ds[[1, 4, 7]] + np.asarray(ds) + f["table"]["id"] # a compound field ``` +- `Dataset.dtype` is the numpy dtype h5py reports: integers and IEEE floats + of every width in either byte order, `bool`, enums (with + `dtype.metadata['enum']`), complex, `S` fixed strings, `object` for + variable-length strings (`bytes` values) and sequences (array values), + `V` opaque, array types, and compounds as structured dtypes. + Other types raise `TypeError`. +- Keys are h5py's: integers, slices with a positive step, `...`, one + increasing list of integers, compound field names. Each maps onto a + hyperslab selection. `None`, negative steps and boolean masks are refused + with h5py's errors. +- What is read from the file: a selection whose bounding box covers at + most half the dataset decodes only the chunks (or contiguous rows) the box + overlaps. The library decodes the whole dataset for a larger box + (including a strided slice such as `ds[::100]` across a chunked dataset), + and for compact, virtual and unwritten datasets and chunked ones with a + non-default fill value. An index list is read one group of neighbouring + chunks at a time (a new group only past a chunk with no selected index), + so each chunk is decoded once. `ds[()]`, `ds[...]` and `np.asarray(ds)` + use the file's chunk cache; other selections do not. +- The bytes the library reads become the numpy array's buffer without a + copy, and the read runs with the GIL released, so threads read in + parallel. A bug in the library (a Rust panic) raises + `clawhdf5.InternalError`, a `RuntimeError`. +- Attributes return what h5py returns; `clawhdf5.Empty` stands for a null + dataspace (h5py's `Empty`). + +## Writing + +`clawhdf5.File(path, "w")` with `create_dataset(name, data=array, +chunks=..., compression="gzip")`, `create_group` and `attrs[...] = ...` +writes `float64`, `float32`, `int64`, `int32` and `uint8` arrays; the file is +written on `close()`. + +## Tests + +```bash +pip install pytest h5py +pytest crates/clawhdf5-py/tests +``` + +`tests/test_read_vs_h5py.py` compares every read with h5py on a file h5py +writes. `scripts/ci-test.sh` builds the wheel and runs these in CI. + ## License MIT diff --git a/crates/clawhdf5-py/pyproject.toml b/crates/clawhdf5-py/pyproject.toml index 87a6f6e..6788adc 100644 --- a/crates/clawhdf5-py/pyproject.toml +++ b/crates/clawhdf5-py/pyproject.toml @@ -3,12 +3,15 @@ requires = ["maturin>=1.0,<2.0"] build-backend = "maturin" [project] -name = "rustyhdf5" +name = "clawhdf5" version = "2.7.0" -description = "Python bindings for rustyhdf5 — a pure-Rust HDF5 library" +description = "Python bindings for clawhdf5 — a pure-Rust HDF5 library" requires-python = ">=3.8" license = { text = "MIT" } dependencies = ["numpy"] [tool.maturin] features = ["extension-module"] +# The extension module is `clawhdf5` (the cdylib's [lib] name): the +# distribution, the import name and the #[pymodule] all agree. +module-name = "clawhdf5" diff --git a/crates/clawhdf5-py/src/attrs.rs b/crates/clawhdf5-py/src/attrs.rs index a0d270e..fd3f0b3 100644 --- a/crates/clawhdf5-py/src/attrs.rs +++ b/crates/clawhdf5-py/src/attrs.rs @@ -1,24 +1,31 @@ //! PyAttrs — dict-like access to HDF5 attributes. -use std::collections::HashMap; use std::sync::{Arc, Mutex}; +use clawhdf5_format::attribute::AttributeMessage; +use pyo3::exceptions::{PyKeyError, PyTypeError, PyValueError}; use pyo3::prelude::*; -use pyo3::types::PyList; +use pyo3::types::{PyList, PyTuple}; -use crate::{OwnedAttrValue, attr_value_to_py, py_to_attr_value}; +use crate::convert::{Converter, Elements, resolve_vl}; +use crate::{OwnedAttrValue, PyEmpty, attr_value_to_py, node, py_to_attr_value}; /// Backing storage for attributes. enum AttrsInner { - /// Read-only attributes from an existing HDF5 object. - Read(HashMap), + /// Attributes of an object in a file opened for reading, sorted by name. + Read { + file: Arc, + attrs: Vec, + }, /// Writable attribute list shared with a parent (PyFile or PyGroup). Write(Arc>>), } /// Dict-like access to HDF5 attributes. /// -/// In read mode, provides immutable access to attribute key/value pairs. +/// In read mode, values are what h5py returns: numpy scalars for scalar +/// attributes, numpy arrays otherwise, `str` for variable-length strings, +/// `numpy.bytes_` for fixed-length ones, and `Empty` for a null dataspace. /// In write mode, attributes set here are accumulated and written when /// the parent file is closed. #[pyclass(name = "Attrs")] @@ -27,11 +34,13 @@ pub struct PyAttrs { } impl PyAttrs { - /// Create a read-only attrs from an existing attribute map. - pub(crate) fn from_read(map: HashMap) -> Self { - Self { - inner: AttrsInner::Read(map), - } + /// The attributes of the object at `addr` (whose path is `path`) in a + /// file opened for reading. + pub(crate) fn read(file: Arc, addr: u64, path: &str) -> PyResult { + let attrs = node::attributes(&file, addr, path)?; + Ok(Self { + inner: AttrsInner::Read { file, attrs }, + }) } /// Create a writable attrs that shares storage with a parent object. @@ -46,11 +55,11 @@ impl PyAttrs { impl PyAttrs { fn __getitem__(&self, py: Python<'_>, key: &str) -> PyResult> { match &self.inner { - AttrsInner::Read(map) => match map.get(key) { - Some(val) => Ok(attr_value_to_py(py, val)), - None => Err(PyErr::new::( - key.to_string(), - )), + AttrsInner::Read { file, attrs } => match attrs.iter().find(|a| a.name == key) { + Some(attr) => Ok(attr_to_py(py, file, attr)?.unbind()), + None => Err(PyKeyError::new_err(format!( + "Can't open attribute (can't locate attribute: '{key}')" + ))), }, AttrsInner::Write(store) => { let guard = store.lock().unwrap(); @@ -60,16 +69,14 @@ impl PyAttrs { return Ok(attr_value_to_py(py, &attr_val)); } } - Err(PyErr::new::( - key.to_string(), - )) + Err(PyKeyError::new_err(key.to_string())) } } } fn __setitem__(&self, key: &str, value: &Bound<'_, PyAny>) -> PyResult<()> { match &self.inner { - AttrsInner::Read(_) => Err(PyErr::new::( + AttrsInner::Read { .. } => Err(PyErr::new::( "cannot set attributes on a read-only file", )), AttrsInner::Write(store) => { @@ -88,14 +95,14 @@ impl PyAttrs { fn __len__(&self) -> usize { match &self.inner { - AttrsInner::Read(map) => map.len(), + AttrsInner::Read { attrs, .. } => attrs.len(), AttrsInner::Write(store) => store.lock().unwrap().len(), } } fn __contains__(&self, key: &str) -> bool { match &self.inner { - AttrsInner::Read(map) => map.contains_key(key), + AttrsInner::Read { attrs, .. } => attrs.iter().any(|a| a.name == key), AttrsInner::Write(store) => store.lock().unwrap().iter().any(|(k, _)| k == key), } } @@ -111,10 +118,20 @@ impl PyAttrs { format!("") } + /// The value of `key`, or `default` if there is no such attribute. + #[pyo3(signature = (key, default=None))] + fn get(&self, py: Python<'_>, key: &str, default: Option>) -> PyResult> { + if self.__contains__(key) { + self.__getitem__(py, key) + } else { + Ok(default.unwrap_or_else(|| py.None())) + } + } + /// Return attribute names as a list. fn keys(&self, py: Python<'_>) -> PyResult> { let names: Vec = match &self.inner { - AttrsInner::Read(map) => map.keys().cloned().collect(), + AttrsInner::Read { attrs, .. } => attrs.iter().map(|a| a.name.clone()).collect(), AttrsInner::Write(store) => store .lock() .unwrap() @@ -129,7 +146,10 @@ impl PyAttrs { /// Return attribute values as a list. fn values(&self, py: Python<'_>) -> PyResult> { let vals: Vec> = match &self.inner { - AttrsInner::Read(map) => map.values().map(|v| attr_value_to_py(py, v)).collect(), + AttrsInner::Read { file, attrs } => attrs + .iter() + .map(|a| attr_to_py(py, file, a).map(Bound::unbind)) + .collect::>()?, AttrsInner::Write(store) => store .lock() .unwrap() @@ -147,10 +167,10 @@ impl PyAttrs { /// Return attribute (key, value) pairs as a list of tuples. fn items(&self, py: Python<'_>) -> PyResult> { let pairs: Vec<(String, Py)> = match &self.inner { - AttrsInner::Read(map) => map + AttrsInner::Read { file, attrs } => attrs .iter() - .map(|(k, v)| (k.clone(), attr_value_to_py(py, v))) - .collect(), + .map(|a| Ok((a.name.clone(), attr_to_py(py, file, a)?.unbind()))) + .collect::>()?, AttrsInner::Write(store) => store .lock() .unwrap() @@ -166,28 +186,69 @@ impl PyAttrs { } } +/// An attribute's value as h5py returns it. +fn attr_to_py<'py>( + py: Python<'py>, + file: &clawhdf5_rs::File, + attr: &AttributeMessage, +) -> PyResult> { + crate::no_panic(|| { + let sb = file.superblock(); + let conv = Converter::new(py, &attr.datatype, sb.offset_size) + .map_err(|e| prefix_err(py, &attr.name, e))?; + if node::is_null(&attr.dataspace) { + return Ok(PyEmpty::new(conv.dtype).into_pyobject(py)?.into_any()); + } + let shape: Vec = attr + .dataspace + .dimensions + .iter() + .map(|&d| d as usize) + .collect(); + let n: usize = shape.iter().product(); + let data = if conv.is_vl() { + let want = n * conv.elem_size; + if attr.raw_data.len() < want { + return Err(PyValueError::new_err(format!( + "attribute {}: {} bytes of variable-length references, expected {want}", + attr.name, + attr.raw_data.len(), + ))); + } + let raw = &attr.raw_data[..want]; + let file_data = file.as_bytes(); + let (osz, lsz, unit) = (sb.offset_size, sb.length_size, conv.vl_unit); + Elements::Vl( + py.detach(|| resolve_vl(file_data, raw, n, osz, lsz, unit)) + .map_err(|e| PyValueError::new_err(format!("attribute {}: {e}", attr.name)))?, + ) + } else { + Elements::Bytes(attr.raw_data.clone()) + }; + let arr = conv + .to_array(py, data, &shape, true) + .map_err(|e| prefix_err(py, &attr.name, e))?; + if shape.is_empty() { + // A scalar dataspace: h5py returns the element itself. + return arr.get_item(PyTuple::empty(py)); + } + Ok(arr) + }) +} + +fn prefix_err(py: Python<'_>, name: &str, e: PyErr) -> PyErr { + let msg = format!("attribute {name}: {}", e.value(py)); + if e.is_instance_of::(py) { + PyTypeError::new_err(msg) + } else { + PyValueError::new_err(msg) + } +} + #[cfg(test)] mod tests { use super::*; - #[test] - fn read_attrs_len() { - let mut map = HashMap::new(); - map.insert("a".into(), clawhdf5_rs::AttrValue::I64(1)); - map.insert("b".into(), clawhdf5_rs::AttrValue::F64(2.0)); - let attrs = PyAttrs::from_read(map); - assert_eq!(attrs.__len__(), 2); - } - - #[test] - fn read_attrs_contains() { - let mut map = HashMap::new(); - map.insert("x".into(), clawhdf5_rs::AttrValue::String("hello".into())); - let attrs = PyAttrs::from_read(map); - assert!(attrs.__contains__("x")); - assert!(!attrs.__contains__("y")); - } - #[test] fn write_attrs_len() { let store = Arc::new(Mutex::new(Vec::new())); diff --git a/crates/clawhdf5-py/src/convert.rs b/crates/clawhdf5-py/src/convert.rs new file mode 100644 index 0000000..4ca7546 --- /dev/null +++ b/crates/clawhdf5-py/src/convert.rs @@ -0,0 +1,600 @@ +//! HDF5 datatypes as numpy dtypes, and element bytes as numpy arrays. +//! +//! The dtype a file's datatype maps to is the one h5py reports for it +//! (byte order kept, compound offsets and padding kept, `r`/`i` compounds as +//! complex, the `FALSE`/`TRUE` enum as `bool`, fixed strings as `S`, +//! variable-length data as `object`). For every fixed-size type that dtype +//! describes the file's element bytes exactly, so the bytes the library +//! returns become the array's buffer as they are: the `Vec` is handed to +//! numpy without a copy and viewed as the dtype. +//! +//! Anything this mapping cannot describe exactly — non-IEEE floats, integers +//! with padding bits, VAX byte order, references, bitfields, time, and +//! variable-length members inside compounds or arrays — is a `TypeError`, +//! never a best-effort guess. + +use std::collections::HashMap; + +use clawhdf5_format::datatype::{CharacterSet, Datatype, DatatypeByteOrder}; +use clawhdf5_format::global_heap::GlobalHeapCollection; +use numpy::PyArray1; +use pyo3::exceptions::{PyTypeError, PyValueError}; +use pyo3::prelude::*; +use pyo3::types::{PyBytes, PyDict, PyList, PyString, PyTuple}; + +/// How the elements of a datatype become Python values. +#[derive(Clone, Debug, PartialEq)] +pub(crate) enum Layout { + /// Fixed-size elements numpy reads as they are. + Fixed, + /// A top-level HDF5 array type: elements are viewed as the base dtype and + /// the array's dimensions are appended to the shape (as h5py does). + Subarray(Vec), + /// Variable-length string: a global heap reference per element. + VlString { utf8: bool }, + /// Variable-length sequence of a fixed-size base type. + VlSequence, +} + +/// Everything needed to turn a dataset's or attribute's bytes into numpy. +pub(crate) struct Converter { + /// The dtype reported to Python (`Dataset.dtype`). + pub dtype: Py, + /// The dtype the element bytes are viewed as: `dtype` itself, the base + /// of a subarray, or the base of a variable-length sequence. + pub view: Py, + pub layout: Layout, + /// Bytes per element in the raw buffer the library returns. + pub elem_size: usize, + /// For variable-length data, bytes per unit of an element's stored + /// length: 1 for strings, the base type's size for sequences. + pub vl_unit: usize, +} + +fn unsupported(what: impl std::fmt::Display) -> PyErr { + PyTypeError::new_err(format!( + "clawhdf5 cannot read this datatype into numpy: {what}" + )) +} + +fn byte_order_char(order: &DatatypeByteOrder, size: u32) -> PyResult<&'static str> { + if size == 1 { + return Ok("|"); + } + match order { + DatatypeByteOrder::LittleEndian => Ok("<"), + DatatypeByteOrder::BigEndian => Ok(">"), + DatatypeByteOrder::Vax => Err(unsupported("VAX byte order")), + } +} + +/// The numpy format string of an integer type, if it is a plain one. +fn int_format(dt: &Datatype) -> PyResult { + match dt { + Datatype::FixedPoint { + size, + byte_order, + signed, + bit_offset, + bit_precision, + } => { + if !matches!(size, 1 | 2 | 4 | 8) { + return Err(unsupported(format!("{size}-byte integer"))); + } + if *bit_offset != 0 || u32::from(*bit_precision) != size * 8 { + return Err(unsupported(format!( + "integer with {bit_precision} significant bits at offset {bit_offset} in {size} bytes" + ))); + } + let kind = if *signed { 'i' } else { 'u' }; + Ok(format!( + "{}{kind}{size}", + byte_order_char(byte_order, *size)? + )) + } + other => Err(unsupported(format!("{other:?} is not an integer"))), + } +} + +/// The numpy format string of an IEEE 754 binary16/32/64 type. +fn float_format(dt: &Datatype) -> PyResult { + let Datatype::FloatingPoint { + size, + byte_order, + bit_offset, + bit_precision, + exponent_location, + exponent_size, + mantissa_location, + mantissa_size, + exponent_bias, + } = dt + else { + return Err(unsupported(format!("{dt:?} is not a float"))); + }; + // (exponent location, exponent size, mantissa size, bias) of IEEE 754. + let ieee = match size { + 2 => (10, 5, 10, 15), + 4 => (23, 8, 23, 127), + 8 => (52, 11, 52, 1023), + _ => return Err(unsupported(format!("{size}-byte float"))), + }; + let layout = ( + *exponent_location, + *exponent_size, + *mantissa_size, + *exponent_bias, + ); + if *bit_offset != 0 + || u32::from(*bit_precision) != size * 8 + || *mantissa_location != 0 + || layout != ieee + { + return Err(unsupported(format!( + "non-IEEE {size}-byte float (exponent {exponent_size} bits at {exponent_location}, \ + mantissa {mantissa_size} bits at {mantissa_location}, bias {exponent_bias})" + ))); + } + Ok(format!("{}f{size}", byte_order_char(byte_order, *size)?)) +} + +/// `r`/`i` compounds of two identical IEEE floats are complex numbers in h5py. +fn complex_format( + size: u32, + members: &[clawhdf5_format::datatype::CompoundMember], +) -> Option { + let [re, im] = members else { return None }; + if re.name != "r" || im.name != "i" || re.datatype != im.datatype { + return None; + } + let Datatype::FloatingPoint { + size: fsize, + byte_order, + .. + } = &re.datatype + else { + return None; + }; + if !matches!(fsize, 4 | 8) + || re.byte_offset != 0 + || im.byte_offset != u64::from(*fsize) + || size != 2 * fsize + { + return None; + } + float_format(&re.datatype).ok()?; + let order = byte_order_char(byte_order, *fsize).ok()?; + Some(format!("{order}c{}", 2 * fsize)) +} + +/// The members of an enum as `{name: value}`. +fn enum_members<'py>( + py: Python<'py>, + base: &Datatype, + members: &[clawhdf5_format::datatype::EnumMember], +) -> PyResult> { + let signed = matches!(base, Datatype::FixedPoint { signed: true, .. }); + let dict = PyDict::new(py); + for m in members { + let value: Py = if signed { + let v = clawhdf5_format::data_read::read_as_i64(&m.value, base) + .map_err(|e| PyValueError::new_err(format!("enum member {}: {e}", m.name)))?; + let v = *v.first().ok_or_else(|| { + PyValueError::new_err(format!("enum member {} has no value", m.name)) + })?; + v.into_pyobject(py)?.into_any().unbind() + } else { + let v = clawhdf5_format::data_read::read_as_u64(&m.value, base) + .map_err(|e| PyValueError::new_err(format!("enum member {}: {e}", m.name)))?; + let v = *v.first().ok_or_else(|| { + PyValueError::new_err(format!("enum member {} has no value", m.name)) + })?; + v.into_pyobject(py)?.into_any().unbind() + }; + dict.set_item(&m.name, value)?; + } + Ok(dict) +} + +/// Whether an enum is h5py's boolean: a one-byte integer with exactly the +/// members `FALSE` = 0 and `TRUE` = 1. +fn is_h5py_bool(base: &Datatype, members: &[clawhdf5_format::datatype::EnumMember]) -> bool { + if base.type_size() != 1 || members.len() != 2 { + return false; + } + let value = |name: &str| { + members + .iter() + .find(|m| m.name == name) + .and_then(|m| m.value.first().copied()) + }; + value("FALSE") == Some(0) && value("TRUE") == Some(1) +} + +fn np_dtype<'py>(py: Python<'py>, spec: impl IntoPyObject<'py>) -> PyResult> { + py.import("numpy")?.getattr("dtype")?.call1((spec,)) +} + +fn np_dtype_with_metadata<'py>( + py: Python<'py>, + spec: impl IntoPyObject<'py>, + metadata: Bound<'py, PyDict>, +) -> PyResult> { + let kwargs = PyDict::new(py); + kwargs.set_item("metadata", metadata)?; + py.import("numpy")? + .getattr("dtype")? + .call((spec,), Some(&kwargs)) +} + +/// The numpy dtype of a fixed-size datatype, whose element bytes numpy can +/// read as they are. +pub(crate) fn fixed_dtype<'py>(py: Python<'py>, dt: &Datatype) -> PyResult> { + match dt { + Datatype::FixedPoint { .. } => np_dtype(py, int_format(dt)?), + Datatype::FloatingPoint { .. } => np_dtype(py, float_format(dt)?), + Datatype::String { size, charset, .. } => { + if *size == 0 { + return Err(unsupported("zero-length fixed string")); + } + let meta = PyDict::new(py); + let enc = match charset { + CharacterSet::Ascii => "ascii", + CharacterSet::Utf8 => "utf-8", + }; + meta.set_item("h5py_encoding", enc)?; + np_dtype_with_metadata(py, format!("S{size}"), meta) + } + Datatype::Opaque { size, .. } => { + if *size == 0 { + return Err(unsupported("zero-length opaque type")); + } + np_dtype(py, format!("V{size}")) + } + Datatype::Enumeration { + base_type, members, .. + } => { + let base = int_format(base_type)?; + if is_h5py_bool(base_type, members) { + return np_dtype(py, "?"); + } + let meta = PyDict::new(py); + meta.set_item("enum", enum_members(py, base_type, members)?)?; + np_dtype_with_metadata(py, base, meta) + } + Datatype::Compound { size, members } => { + if let Some(c) = complex_format(*size, members) { + return np_dtype(py, c); + } + let names = PyList::empty(py); + let formats = PyList::empty(py); + let offsets = PyList::empty(py); + for m in members { + let end = m.byte_offset.checked_add(u64::from(m.datatype.type_size())); + if end.is_none_or(|end| end > u64::from(*size)) { + return Err(PyValueError::new_err(format!( + "compound member {} lies outside the {size}-byte compound", + m.name + ))); + } + names.append(&m.name)?; + formats.append(fixed_dtype(py, &m.datatype).map_err(|e| { + unsupported(format!("compound member {}: {}", m.name, e.value(py))) + })?)?; + offsets.append(m.byte_offset)?; + } + let spec = PyDict::new(py); + spec.set_item("names", names)?; + spec.set_item("formats", formats)?; + spec.set_item("offsets", offsets)?; + spec.set_item("itemsize", size)?; + np_dtype(py, spec) + } + Datatype::Array { + base_type, + dimensions, + } => { + let base = fixed_dtype(py, base_type)?; + let dims = PyTuple::new(py, dimensions)?; + np_dtype(py, (base, dims)) + } + Datatype::VariableLength { is_string, .. } => Err(unsupported(if *is_string { + "variable-length string inside a compound or array type" + } else { + "variable-length sequence inside a compound or array type" + })), + Datatype::Reference { .. } => Err(unsupported("object/region references")), + Datatype::BitField { .. } => Err(unsupported("bitfield")), + Datatype::Time { .. } => Err(unsupported("time")), + } +} + +impl Converter { + /// The converter for a dataset's or attribute's datatype. + pub(crate) fn new(py: Python<'_>, dt: &Datatype, offset_size: u8) -> PyResult { + match dt { + Datatype::VariableLength { + is_string: true, + charset, + .. + } => { + let utf8 = matches!(charset, Some(CharacterSet::Utf8)); + let meta = PyDict::new(py); + if utf8 { + meta.set_item("vlen", py.get_type::())?; + } else { + meta.set_item("vlen", py.get_type::())?; + } + let dtype = np_dtype_with_metadata(py, "O", meta)?; + Ok(Self { + view: dtype.clone().unbind(), + dtype: dtype.unbind(), + layout: Layout::VlString { utf8 }, + elem_size: vl_ref_size(offset_size)?, + vl_unit: 1, + }) + } + Datatype::VariableLength { + is_string: false, + base_type, + .. + } => { + let base = fixed_dtype(py, base_type)?; + let meta = PyDict::new(py); + meta.set_item("vlen", &base)?; + let dtype = np_dtype_with_metadata(py, "O", meta)?; + Ok(Self { + dtype: dtype.unbind(), + view: base.unbind(), + layout: Layout::VlSequence, + elem_size: vl_ref_size(offset_size)?, + vl_unit: base_type.type_size() as usize, + }) + } + Datatype::Array { + base_type, + dimensions, + } => { + let dtype = fixed_dtype(py, dt)?; + let base = fixed_dtype(py, base_type)?; + Ok(Self { + dtype: dtype.unbind(), + view: base.unbind(), + layout: Layout::Subarray(dimensions.iter().map(|&d| d as usize).collect()), + elem_size: dt.type_size() as usize, + vl_unit: 0, + }) + } + _ => { + let dtype = fixed_dtype(py, dt)?; + Ok(Self { + view: dtype.clone().unbind(), + dtype: dtype.unbind(), + layout: Layout::Fixed, + elem_size: dt.type_size() as usize, + vl_unit: 0, + }) + } + } + } + + pub(crate) fn is_vl(&self) -> bool { + matches!(self.layout, Layout::VlString { .. } | Layout::VlSequence) + } + + /// An empty array of `shape` (some dimension is zero). + pub(crate) fn empty<'py>( + &self, + py: Python<'py>, + shape: &[usize], + ) -> PyResult> { + let np = py.import("numpy")?; + match &self.layout { + Layout::Subarray(dims) => { + let mut full = shape.to_vec(); + full.extend_from_slice(dims); + np.call_method1("empty", (PyTuple::new(py, full)?, self.view.bind(py))) + } + _ => np.call_method1("empty", (PyTuple::new(py, shape)?, self.dtype.bind(py))), + } + } + + /// Turn decoded element data into a numpy array of `shape`. + /// + /// `str_values` decodes variable-length strings to `str` (what h5py + /// does for attributes) instead of `bytes` (what it does for datasets). + pub(crate) fn to_array<'py>( + &self, + py: Python<'py>, + data: Elements, + shape: &[usize], + str_values: bool, + ) -> PyResult> { + let n: usize = shape.iter().product(); + match (data, &self.layout) { + (Elements::Bytes(bytes), Layout::Fixed) => { + bytes_as_array(py, bytes, self.view.bind(py), shape) + } + (Elements::Bytes(bytes), Layout::Subarray(dims)) => { + let mut full = shape.to_vec(); + full.extend_from_slice(dims); + bytes_as_array(py, bytes, self.view.bind(py), &full) + } + (Elements::Vl(items), Layout::VlString { .. }) => { + check_count(items.len(), n)?; + let mut objs: Vec> = Vec::with_capacity(items.len()); + for item in items { + let obj = if str_values { + PyBytes::new(py, &item) + .call_method1("decode", ("utf-8", "surrogateescape"))? + .unbind() + } else { + PyBytes::new(py, &item).into_any().unbind() + }; + objs.push(obj); + } + object_array(py, objs, shape) + } + (Elements::Vl(items), Layout::VlSequence) => { + check_count(items.len(), n)?; + let base = self.view.bind(py); + let itemsize: usize = base.getattr("itemsize")?.extract()?; + let mut objs: Vec> = Vec::with_capacity(items.len()); + for item in items { + if item.len() % itemsize != 0 { + return Err(PyValueError::new_err(format!( + "variable-length element of {} bytes is not a whole number of {itemsize}-byte values", + item.len() + ))); + } + let len = item.len() / itemsize; + objs.push(bytes_as_array(py, item, base, &[len])?.unbind()); + } + object_array(py, objs, shape) + } + _ => Err(PyValueError::new_err( + "internal error: element data does not match the datatype", + )), + } + } +} + +/// Element data as read, before it becomes numpy. +pub(crate) enum Elements { + /// The elements' bytes, back to back. + Bytes(Vec), + /// Each variable-length element's bytes, resolved from the global heap. + Vl(Vec>), +} + +fn vl_ref_size(offset_size: u8) -> PyResult { + // The library sizes a variable-length element as 16 bytes (a length, an + // 8-byte heap address and an index) whatever the file's offset size. + // Refuse the other sizes rather than read misaligned references. + if offset_size != 8 { + return Err(unsupported(format!( + "variable-length data in a file with {offset_size}-byte offsets" + ))); + } + Ok(4 + usize::from(offset_size) + 4) +} + +fn check_count(got: usize, want: usize) -> PyResult<()> { + if got != want { + return Err(PyValueError::new_err(format!( + "read {got} elements, expected {want}" + ))); + } + Ok(()) +} + +/// A numpy array over `bytes` without copying them: the `Vec` becomes the +/// array's buffer and is viewed as `dtype` with `shape`. +pub(crate) fn bytes_as_array<'py>( + py: Python<'py>, + bytes: Vec, + dtype: &Bound<'py, PyAny>, + shape: &[usize], +) -> PyResult> { + let itemsize: usize = dtype.getattr("itemsize")?.extract()?; + let n: usize = shape.iter().product(); + if n.checked_mul(itemsize) != Some(bytes.len()) { + return Err(PyValueError::new_err(format!( + "read {} bytes, expected {n} elements of {itemsize} bytes", + bytes.len() + ))); + } + let shape = PyTuple::new(py, shape)?; + if n == 0 { + return py.import("numpy")?.call_method1("empty", (shape, dtype)); + } + let raw = PyArray1::from_vec(py, bytes); + let arr = raw + .call_method1("view", (dtype,))? + .call_method1("reshape", (shape,))?; + // A `Vec` carries no alignment promise. numpy copes with unaligned + // arrays, but slowly and not in every routine, so hand out an aligned + // copy in the (allocator-dependent, rare) case the buffer is not. + if !arr + .getattr("flags")? + .getattr("aligned")? + .extract::()? + { + return arr.call_method0("copy"); + } + Ok(arr) +} + +fn object_array<'py>( + py: Python<'py>, + objs: Vec>, + shape: &[usize], +) -> PyResult> { + let arr = PyArray1::from_vec(py, objs); + arr.call_method1("reshape", (PyTuple::new(py, shape)?,)) +} + +/// Resolve variable-length elements (global heap references in `raw`) to +/// their bytes: each element's stored length times `unit` (1 for strings, +/// the base type's size for sequences). Pure Rust, so it runs without the +/// GIL. +pub(crate) fn resolve_vl( + file_data: &[u8], + raw: &[u8], + count: usize, + offset_size: u8, + length_size: u8, + unit: usize, +) -> Result>, String> { + let refs = clawhdf5_format::vl_data::parse_vl_references(raw, count as u64, offset_size) + .map_err(|e| e.to_string())?; + let undefined = match offset_size { + 2 => 0xFFFF, + 4 => 0xFFFF_FFFF, + _ => u64::MAX, + }; + let mut collections: HashMap = HashMap::new(); + let mut out = Vec::with_capacity(refs.len()); + for vl in &refs { + if vl.collection_address == 0 || vl.collection_address == undefined { + if vl.length != 0 { + return Err(format!( + "variable-length element of length {} has no heap address", + vl.length + )); + } + out.push(Vec::new()); + continue; + } + let coll = match collections.entry(vl.collection_address) { + std::collections::hash_map::Entry::Occupied(e) => e.into_mut(), + std::collections::hash_map::Entry::Vacant(e) => { + let addr = usize::try_from(vl.collection_address) + .map_err(|_| "global heap address out of range".to_string())?; + e.insert( + GlobalHeapCollection::parse(file_data, addr, length_size) + .map_err(|e| e.to_string())?, + ) + } + }; + let index = u16::try_from(vl.object_index) + .map_err(|_| format!("global heap object index {} out of range", vl.object_index))?; + let obj = coll.get_object(index).ok_or_else(|| { + format!( + "global heap object {index} not found in the collection at {}", + vl.collection_address + ) + })?; + let need = (vl.length as usize) + .checked_mul(unit) + .ok_or("variable-length element too long")?; + if need > obj.data.len() { + return Err(format!( + "variable-length element of {need} bytes in a {}-byte heap object", + obj.data.len() + )); + } + out.push(obj.data[..need].to_vec()); + } + Ok(out) +} diff --git a/crates/clawhdf5-py/src/dataset.rs b/crates/clawhdf5-py/src/dataset.rs index e1b7a38..bcfd852 100644 --- a/crates/clawhdf5-py/src/dataset.rs +++ b/crates/clawhdf5-py/src/dataset.rs @@ -1,241 +1,373 @@ -//! PyDataset — read access to HDF5 datasets with numpy integration. +//! PyDataset — h5py-style read access to HDF5 datasets. +//! +//! `ds[key]` parses the key into hyperslab selections (see `select`) and +//! reads them through the facade's `read_selection`, which decodes only the +//! chunks a small selection touches (see its docs for when it decodes the +//! whole dataset instead); the +//! bytes it returns become the numpy array's buffer without a copy (see +//! `convert`). All file access and decoding runs with the GIL released, so +//! Python threads reading the same or different datasets run in parallel. use std::sync::Arc; -use numpy::PyArrayDyn; -use numpy::ndarray::{ArrayD, IxDyn}; +use clawhdf5_format::datatype::Datatype; +use clawhdf5_format::object_header::ObjectHeader; +use pyo3::exceptions::{PyTypeError, PyValueError}; use pyo3::prelude::*; -use pyo3::types::PyList; - -use clawhdf5_rs::DType; +use pyo3::types::{PyList, PyTuple}; use crate::attrs::PyAttrs; -use crate::to_py_err; +use crate::convert::{Converter, Elements, resolve_vl}; +use crate::select::{self, Plan}; +use crate::{PyEmpty, node, to_py_err}; -/// A handle to an HDF5 dataset (read mode). +/// A dataset in a file opened for reading. /// -/// Supports numpy-style indexing via `__getitem__`: /// ```python -/// ds = f['dataset_name'] -/// data = ds[:] # read all data as numpy array -/// shape = ds.shape -/// dtype = ds.dtype +/// ds = f['group/dataset'] +/// ds.shape, ds.dtype, ds.attrs['units'] +/// block = ds[10:20, ::2] # a small selection reads only its chunks /// ``` #[pyclass(name = "Dataset")] pub struct PyDataset { file: Arc, path: String, - cached_shape: Vec, - cached_dtype: DType, + /// Where the dataset's object header is: reads open it from here rather + /// than resolve `path` again. + addr: u64, + /// `None` for a dataset with a null dataspace (h5py's `Empty`). + shape: Option>, + /// The chunk shape, for a chunked dataset. + chunks: Option>, + datatype: Datatype, + /// Why the datatype cannot be read into numpy, if it cannot. + conv: Result, } impl PyDataset { - pub fn new(file: Arc, path: String) -> PyResult { - let ds = file.dataset(&path).map_err(to_py_err)?; - let cached_shape = ds.shape().map_err(to_py_err)?; - let cached_dtype = ds.dtype().map_err(to_py_err)?; - Ok(Self { - file, - path, - cached_shape, - cached_dtype, + pub(crate) fn open( + py: Python<'_>, + file: Arc, + path: String, + addr: u64, + hdr: &ObjectHeader, + ) -> PyResult { + crate::no_panic(|| { + let null = node::is_null(&node::dataspace(&file, hdr)?); + let (shape, datatype) = { + let ds = file.dataset_at(addr).map_err(to_py_err)?; + let shape = if null { + None + } else { + Some(ds.shape().map_err(to_py_err)?) + }; + (shape, ds.raw_datatype().map_err(to_py_err)?) + }; + let conv = Converter::new(py, &datatype, file.superblock().offset_size) + .map_err(|e| e.value(py).to_string()); + let chunks = shape + .as_ref() + .and_then(|s| node::chunk_shape(&file, hdr, s.len())); + Ok(Self { + file, + path, + addr, + shape, + chunks, + datatype, + conv, + }) }) } -} -/// Map a `DType` to a numpy dtype string. -fn dtype_to_numpy_str(dt: &DType) -> &'static str { - match dt { - DType::F64 => "float64", - DType::F32 => "float32", - DType::I64 => "int64", - DType::I32 => "int32", - DType::I16 => "int16", - DType::I8 => "int8", - DType::U64 => "uint64", - DType::U32 => "uint32", - DType::U16 => "uint16", - DType::U8 => "uint8", - DType::String | DType::VariableLengthString => "object", - _ => "object", - } -} - -#[pymethods] -impl PyDataset { - /// The shape of the dataset as a tuple. - #[getter] - fn shape(&self, py: Python<'_>) -> PyResult> { - let tuple = pyo3::types::PyTuple::new(py, self.cached_shape.iter().map(|&d| d as usize))?; - Ok(tuple.into_any().unbind()) + fn converter(&self) -> PyResult<&Converter> { + self.conv + .as_ref() + .map_err(|msg| PyTypeError::new_err(format!("{}: {msg}", node::name(&self.path)))) } - /// The numpy dtype string of the dataset. - #[getter] - fn dtype(&self) -> &'static str { - dtype_to_numpy_str(&self.cached_dtype) - } + /// Read the selection described by `plan` into a numpy array. + fn read_plan<'py>(&self, py: Python<'py>, plan: &Plan) -> PyResult> { + let conv = self.converter()?; + let dims = self.shape.as_deref().unwrap_or(&[]); + let out_shape = plan.out_shape(); - /// Attribute access (read-only). - #[getter] - fn attrs(&self) -> PyResult { - let ds = self.file.dataset(&self.path).map_err(to_py_err)?; - let map = ds.attrs().map_err(to_py_err)?; - Ok(PyAttrs::from_read(map)) - } - - /// Read data via indexing. Supports `ds[:]`, `ds[0]`, `ds[0:5]`, etc. - /// - /// The full dataset is always read from the underlying file; the index - /// is then applied on the resulting numpy array. - fn __getitem__<'py>(&self, py: Python<'py>, key: &Bound<'py, PyAny>) -> PyResult> { - let arr = self.read_as_numpy(py)?; - let indexed = arr.get_item(key)?; - Ok(indexed.unbind()) - } - - fn __repr__(&self) -> String { - format!( - "", - self.path, - self.cached_shape, - dtype_to_numpy_str(&self.cached_dtype), - ) - } - - fn __len__(&self) -> usize { - self.cached_shape.first().copied().unwrap_or(0) as usize - } -} - -impl PyDataset { - /// Read the full dataset and return it as a numpy array (or list for strings). - /// - /// For numeric types, the Rust I/O (file reading + decompression) is - /// performed inside `py.detach()` so that the GIL is released - /// during the potentially expensive operation. The numpy array - /// construction still happens with the GIL held. - fn read_as_numpy<'py>(&self, py: Python<'py>) -> PyResult> { - let file = &self.file; - let path = &self.path; - let shape: Vec = self.cached_shape.iter().map(|&d| d as usize).collect(); - - match &self.cached_dtype { - DType::F64 => { - let data = py - .detach(|| file.dataset(path).and_then(|ds| ds.read_f64())) - .map_err(to_py_err)?; - let nd = ArrayD::from_shape_vec(IxDyn(&shape), data) - .map_err(|e| PyErr::new::(e.to_string()))?; - let arr = PyArrayDyn::from_owned_array(py, nd); - Ok(arr.into_any()) - } - DType::F32 => { - let data = py - .detach(|| file.dataset(path).and_then(|ds| ds.read_f32())) - .map_err(to_py_err)?; - let nd = ArrayD::from_shape_vec(IxDyn(&shape), data) - .map_err(|e| PyErr::new::(e.to_string()))?; - let arr = PyArrayDyn::from_owned_array(py, nd); - Ok(arr.into_any()) - } - DType::I32 => { - let data = py - .detach(|| file.dataset(path).and_then(|ds| ds.read_i32())) - .map_err(to_py_err)?; - let nd = ArrayD::from_shape_vec(IxDyn(&shape), data) - .map_err(|e| PyErr::new::(e.to_string()))?; - let arr = PyArrayDyn::from_owned_array(py, nd); - Ok(arr.into_any()) - } - DType::I64 => { - let data = py - .detach(|| file.dataset(path).and_then(|ds| ds.read_i64())) - .map_err(to_py_err)?; - let nd = ArrayD::from_shape_vec(IxDyn(&shape), data) - .map_err(|e| PyErr::new::(e.to_string()))?; - let arr = PyArrayDyn::from_owned_array(py, nd); - Ok(arr.into_any()) - } - DType::U8 => { - // Try zero-copy first (contiguous layout), fall back to - // read_u64 + cast for chunked/compact datasets. - let data: Vec = py - .detach(|| { - let ds = file.dataset(path)?; - match ds.read_u8_zerocopy() { - Ok(slice) => Ok(slice.to_vec()), - Err(_) => { - let raw = ds.read_u64()?; - Ok(raw.iter().map(|&v| v as u8).collect()) - } + let arr = if plan.is_empty() { + conv.empty(py, &out_shape)? + } else { + let (vl, elem_size, unit) = (conv.is_vl(), conv.elem_size, conv.vl_unit); + let list_axis = plan.list_axis(); + let chunk_len = match (&self.chunks, list_axis) { + (Some(c), Some(a)) => c.get(a).copied(), + _ => None, + }; + let (reads, list_axis) = plan.reads(dims, chunk_len, elem_size); + let read_shape = plan.read_shape(); + let file = &*self.file; + let addr = self.addr; + // Everything below touches only Rust data: release the GIL. + let read = || -> Result { + let ds = file.dataset_at(addr)?; + let mut blocks = Vec::with_capacity(reads.len()); + for read in reads { + let raw = ds.read_selection(&read.sel)?; + let mut shape = read.shape; + let want = shape.iter().product::() * elem_size; + if raw.len() != want { + return Err(ReadError::Other(format!( + "read {} bytes, expected {want}", + raw.len() + ))); + } + let raw = match (&read.pick, list_axis) { + (Some(pick), Some(axis)) => { + let kept = select::gather_along(&raw, &shape, axis, pick, elem_size); + shape[axis] = pick.len(); + kept } - }) - .map_err(to_py_err)?; - let nd = ArrayD::from_shape_vec(IxDyn(&shape), data) - .map_err(|e| PyErr::new::(e.to_string()))?; - let arr = PyArrayDyn::from_owned_array(py, nd); - Ok(arr.into_any()) + _ => raw, + }; + blocks.push((raw, shape)); + } + // Several blocks only for a list index: join their bytes + // (every byte of every element, padding included) along + // that axis before anything becomes numpy. + let raw = match (blocks.len(), list_axis) { + (1, _) => blocks.pop().expect("one block").0, + (_, Some(axis)) => select::join_along(&blocks, axis, elem_size), + _ => { + return Err(ReadError::Other( + "several reads without an index list".into(), + )); + } + }; + if !vl { + return Ok(Elements::Bytes(raw)); + } + let sb = file.superblock(); + let n = read_shape.iter().product(); + resolve_vl( + file.as_bytes(), + &raw, + n, + sb.offset_size, + sb.length_size, + unit, + ) + .map(Elements::Vl) + .map_err(ReadError::Other) + }; + let data = py + .detach(|| { + std::panic::catch_unwind(std::panic::AssertUnwindSafe(read)) + .unwrap_or_else(|p| Err(ReadError::Panic(crate::panic_text(&*p)))) + }) + .map_err(|e| e.into_py(&self.path))?; + let joined = conv.to_array(py, data, &read_shape, false)?; + // Drop the axes indexed by an integer (length 1 in the blocks). + let mut shape = out_shape.clone(); + if let crate::convert::Layout::Subarray(sub) = &conv.layout { + shape.extend_from_slice(sub); } - DType::U64 => { - let data = py - .detach(|| file.dataset(path).and_then(|ds| ds.read_u64())) - .map_err(to_py_err)?; - let nd = ArrayD::from_shape_vec(IxDyn(&shape), data) - .map_err(|e| PyErr::new::(e.to_string()))?; - let arr = PyArrayDyn::from_owned_array(py, nd); - Ok(arr.into_any()) - } - DType::String | DType::VariableLengthString => { - // String reads need the GIL for PyList construction, but we - // release it during the Rust I/O portion. - let data = py - .detach(|| file.dataset(path).and_then(|ds| ds.read_string())) - .map_err(to_py_err)?; - let list = PyList::new(py, &data)?; - Ok(list.into_any()) - } - other => Err(PyErr::new::(format!( - "unsupported dataset dtype for reading: {other}" - ))), + joined.call_method1("reshape", (PyTuple::new(py, shape)?,))? + }; + + let arr = select_fields(py, arr, &plan.fields)?; + if plan.scalar { + return arr.get_item(PyTuple::empty(py)); + } + Ok(arr) + } +} + +/// An error from the read closure, turned into a Python error with the GIL. +enum ReadError { + Lib(clawhdf5_rs::Error), + Other(String), + Panic(String), +} + +impl From for ReadError { + fn from(e: clawhdf5_rs::Error) -> Self { + ReadError::Lib(e) + } +} + +impl ReadError { + fn into_py(self, path: &str) -> PyErr { + match self { + ReadError::Lib(e) => to_py_err(e), + ReadError::Other(msg) => PyValueError::new_err(format!("{}: {msg}", node::name(path))), + ReadError::Panic(msg) => crate::InternalError::new_err(format!( + "{}: clawhdf5 internal error (please report it): {msg}", + node::name(path) + )), } } } -#[cfg(test)] -mod tests { - use super::*; +/// Keep only the named compound fields, as h5py's `ds['x']` / `ds['x', 'y']`. +fn select_fields<'py>( + py: Python<'py>, + arr: Bound<'py, PyAny>, + fields: &[String], +) -> PyResult> { + if fields.is_empty() { + return Ok(arr); + } + let names = arr.getattr("dtype")?.getattr("names")?; + if names.is_none() { + return Err(PyValueError::new_err( + "Field names only allowed for compound types", + )); + } + let names: Vec = names.extract()?; + for f in fields { + if !names.contains(f) { + return Err(PyValueError::new_err(format!( + "Field {f} does not appear in this type." + ))); + } + } + let np = py.import("numpy")?; + if let [one] = fields { + return np.call_method1("ascontiguousarray", (arr.get_item(one)?,)); + } + let picked = arr.get_item(PyList::new(py, fields)?)?; + py.import("numpy.lib.recfunctions")? + .call_method1("repack_fields", (picked,)) +} - #[test] - fn dtype_mapping() { - assert_eq!(dtype_to_numpy_str(&DType::F64), "float64"); - assert_eq!(dtype_to_numpy_str(&DType::F32), "float32"); - assert_eq!(dtype_to_numpy_str(&DType::I32), "int32"); - assert_eq!(dtype_to_numpy_str(&DType::I64), "int64"); - assert_eq!(dtype_to_numpy_str(&DType::U8), "uint8"); - assert_eq!(dtype_to_numpy_str(&DType::String), "object"); +#[pymethods] +impl PyDataset { + /// The shape of the dataset (`None` for an empty/null dataspace). + #[getter] + fn shape<'py>(&self, py: Python<'py>) -> PyResult> { + match &self.shape { + Some(s) => Ok(PyTuple::new(py, s)?.into_any()), + None => Ok(py.None().into_bound(py)), + } } - #[test] - fn dataset_from_file() { - let mut b = clawhdf5_rs::FileBuilder::new(); - b.create_dataset("vals").with_f64_data(&[1.0, 2.0, 3.0]); - let bytes = b.finish().unwrap(); - let file = Arc::new(clawhdf5_rs::File::from_bytes(bytes).unwrap()); - let ds = PyDataset::new(file, "vals".into()).unwrap(); - assert_eq!(ds.cached_shape, vec![3]); - assert_eq!(ds.cached_dtype, DType::F64); + /// The maximum shape (`None` per unlimited dimension), like h5py. + #[getter] + fn maxshape<'py>(&self, py: Python<'py>) -> PyResult> { + crate::no_panic(|| { + let Some(shape) = &self.shape else { + return Ok(py.None().into_bound(py)); + }; + let max = self + .file + .dataset_at(self.addr) + .and_then(|ds| ds.max_dimensions()) + .map_err(to_py_err)? + .unwrap_or_else(|| shape.clone()); + let items: Vec> = max + .into_iter() + .map(|d| (d != u64::MAX).then_some(d)) + .collect(); + Ok(PyTuple::new(py, items)?.into_any()) + }) } - #[test] - fn dataset_len() { - let mut b = clawhdf5_rs::FileBuilder::new(); - b.create_dataset("data") - .with_i32_data(&[10, 20, 30, 40]) - .with_shape(&[2, 2]); - let bytes = b.finish().unwrap(); - let file = Arc::new(clawhdf5_rs::File::from_bytes(bytes).unwrap()); - let ds = PyDataset::new(file, "data".into()).unwrap(); - assert_eq!(ds.__len__(), 2); + /// The dataset's numpy dtype, as h5py reports it. + #[getter] + fn dtype<'py>(&self, py: Python<'py>) -> PyResult> { + Ok(self.converter()?.dtype.bind(py).clone()) + } + + #[getter] + fn ndim(&self) -> usize { + self.shape.as_ref().map_or(0, Vec::len) + } + + /// Number of elements (`None` for an empty/null dataspace, as h5py). + #[getter] + fn size(&self) -> Option { + self.shape.as_ref().map(|s| s.iter().product()) + } + + /// The dataset's full name, e.g. `/group/data`. + #[getter] + fn name(&self) -> String { + node::name(&self.path) + } + + /// The dataset's attributes (read-only, dict-like). + #[getter] + fn attrs(&self) -> PyResult { + PyAttrs::read(Arc::clone(&self.file), self.addr, &self.path) + } + + /// Read with h5py indexing: integers, slices with positive steps, + /// `...`, one increasing list of integers, and compound field names. + /// A selection whose bounding box covers at most half the dataset reads + /// only the chunks (or contiguous rows) it overlaps. + fn __getitem__<'py>( + &self, + py: Python<'py>, + key: &Bound<'py, PyAny>, + ) -> PyResult> { + let Some(dims) = &self.shape else { + let is_empty_tuple = key.cast::().is_ok_and(|t| t.is_empty()); + let is_ellipsis = key.is_instance_of::(); + if is_empty_tuple || is_ellipsis { + let empty = PyEmpty::new(self.converter()?.dtype.clone_ref(py)); + return Ok(empty.into_pyobject(py)?.into_any()); + } + return Err(PyValueError::new_err("Empty datasets cannot be sliced")); + }; + let plan = select::parse(key, dims)?; + self.read_plan(py, &plan) + } + + /// `numpy.asarray(ds)` reads the whole dataset. + #[pyo3(signature = (dtype=None, copy=None))] + fn __array__<'py>( + &self, + py: Python<'py>, + dtype: Option<&Bound<'py, PyAny>>, + copy: Option, + ) -> PyResult> { + let _ = copy; // every read is a fresh array + let Some(dims) = &self.shape else { + return Err(PyValueError::new_err("an empty dataset has no array value")); + }; + let ellipsis = pyo3::types::PyEllipsis::get(py).to_owned().into_any(); + let plan = select::parse(&ellipsis, dims)?; + let arr = self.read_plan(py, &plan)?; + match dtype { + Some(dt) => arr.call_method1("astype", (dt,)), + None => Ok(arr), + } + } + + fn __len__(&self) -> PyResult { + match self.shape.as_deref() { + Some([first, ..]) => Ok(*first as usize), + _ => Err(PyTypeError::new_err( + "Attempt to take len() of scalar dataset", + )), + } + } + + fn __repr__(&self, py: Python<'_>) -> String { + let dtype = match &self.conv { + Ok(c) => c + .dtype + .bind(py) + .str() + .map(|s| s.to_string()) + .unwrap_or_default(), + Err(_) => format!("{:?}", self.datatype), + }; + let shape = match &self.shape { + Some(s) => format!("{s:?}"), + None => "None".to_string(), + }; + format!( + "", + node::name(&self.path) + ) } } diff --git a/crates/clawhdf5-py/src/file.rs b/crates/clawhdf5-py/src/file.rs index 23e5f95..5f7f51b 100644 --- a/crates/clawhdf5-py/src/file.rs +++ b/crates/clawhdf5-py/src/file.rs @@ -4,10 +4,10 @@ use std::path::PathBuf; use std::sync::{Arc, Mutex}; use pyo3::prelude::*; +use pyo3::types::PyList; use crate::attrs::PyAttrs; -use crate::dataset::PyDataset; -use crate::group::{PyGroup, WriteGroupState, finalize_write_group}; +use crate::group::{PyGroup, ReadGroup, WriteGroupState, finalize_write_group}; use crate::{DatasetSpec, OwnedAttrValue, apply_dataset_spec, extract_numpy_data, to_py_err}; /// Internal state for write mode. @@ -35,10 +35,12 @@ struct WriteState { #[pyclass(name = "File")] pub struct PyFile { inner: Option, + filename: String, } enum FileInner { - Read(Arc), + /// The root group; it holds the file. + Read(ReadGroup), Write(WriteState), } @@ -51,15 +53,20 @@ impl PyFile { /// mode: 'r' for read (default), 'w' for write #[new] #[pyo3(signature = (path, mode="r"))] - fn new(path: &str, mode: &str) -> PyResult { + fn new(py: Python<'_>, path: &str, mode: &str) -> PyResult { + let filename = path.to_string(); match mode { "r" => { - let file = clawhdf5_rs::File::open(path).map_err(to_py_err)?; + let file = py.detach(|| { + crate::no_panic(|| clawhdf5_rs::File::open(path).map_err(to_py_err)) + })?; Ok(Self { - inner: Some(FileInner::Read(Arc::new(file))), + inner: Some(FileInner::Read(root_group(Arc::new(file)))), + filename, }) } "w" => Ok(Self { + filename, inner: Some(FileInner::Write(WriteState { path: PathBuf::from(path), root_datasets: Vec::new(), @@ -101,44 +108,51 @@ impl PyFile { Ok(false) // don't suppress exceptions } - /// Get a child object (dataset or group) by path. + /// Get a child object (dataset or group) by path; `f['/']` is the root. fn __getitem__(&self, py: Python<'_>, key: &str) -> PyResult> { - let file = self.read_file()?; - // Try dataset first - match file.dataset(key) { - Ok(_) => { - let ds = PyDataset::new(Arc::clone(file), key.to_string())?; - Ok(ds.into_pyobject(py)?.into_any().unbind()) - } - Err(clawhdf5_rs::Error::NotADataset(_)) => { - let grp = PyGroup::from_read(Arc::clone(file), key.to_string()); - Ok(grp.into_pyobject(py)?.into_any().unbind()) - } - Err(_) => { - // Could be a group (no DataLayout message, no error) - match file.group(key) { - Ok(_) => { - let grp = PyGroup::from_read(Arc::clone(file), key.to_string()); - Ok(grp.into_pyobject(py)?.into_any().unbind()) - } - Err(e) => Err(PyErr::new::(format!( - "{key}: {e}" - ))), - } - } - } + self.read_file()?.get_item(py, key) + } + + /// `f.get(key, default=None)`. + #[pyo3(signature = (key, default=None))] + fn get(&self, py: Python<'_>, key: &str, default: Option>) -> PyResult> { + self.read_file()?.get(py, key, default) } /// List the names of all children in the root group. fn keys(&self, py: Python<'_>) -> PyResult> { - let file = self.read_file()?; - let root = file.root(); - let mut names = root.datasets().map_err(to_py_err)?; - let groups = root.groups().map_err(to_py_err)?; - names.extend(groups); - names.sort(); - let list = pyo3::types::PyList::new(py, &names)?; - Ok(list.into_any().unbind()) + let names = self.read_file()?.member_names()?; + Ok(PyList::new(py, names)?.into_any().unbind()) + } + + fn values(&self, py: Python<'_>) -> PyResult> { + let vals = self.read_file()?.values(py)?; + Ok(PyList::new(py, vals)?.into_any().unbind()) + } + + fn items(&self, py: Python<'_>) -> PyResult> { + let items = self.read_file()?.items(py)?; + Ok(PyList::new(py, items)?.into_any().unbind()) + } + + fn __iter__(&self, py: Python<'_>) -> PyResult> { + self.keys(py)?.call_method0(py, "__iter__") + } + + fn __len__(&self) -> PyResult { + Ok(self.read_file()?.member_names()?.len()) + } + + /// The root group's name, `/`. + #[getter] + fn name(&self) -> &'static str { + "/" + } + + /// The path the file was opened with. + #[getter] + fn filename(&self) -> &str { + &self.filename } /// Create a dataset in the root group (write mode only). @@ -192,10 +206,7 @@ impl PyFile { #[getter] fn attrs(&self) -> PyResult { match self.inner.as_ref() { - Some(FileInner::Read(file)) => { - let map = file.root().attrs().map_err(to_py_err)?; - Ok(PyAttrs::from_read(map)) - } + Some(FileInner::Read(root)) => root.attrs(), Some(FileInner::Write(state)) => Ok(PyAttrs::from_write(Arc::clone(&state.root_attrs))), None => Err(PyErr::new::( "file is closed", @@ -205,8 +216,8 @@ impl PyFile { fn __repr__(&self) -> String { match &self.inner { - Some(FileInner::Read(f)) => { - format!("", f.as_bytes().len()) + Some(FileInner::Read(root)) => { + format!("", root.file.as_bytes().len()) } Some(FileInner::Write(s)) => { format!("", s.path.display()) @@ -216,13 +227,13 @@ impl PyFile { } fn __contains__(&self, key: &str) -> PyResult { - let file = self.read_file()?; - Ok(file.dataset(key).is_ok() || file.group(key).is_ok()) + Ok(self.read_file()?.contains(key)) } } impl PyFile { - fn read_file(&self) -> PyResult<&Arc> { + /// The root group of a file opened for reading. + fn read_file(&self) -> PyResult<&ReadGroup> { match &self.inner { Some(FileInner::Read(f)) => Ok(f), Some(FileInner::Write(_)) => Err(PyErr::new::( @@ -260,31 +271,38 @@ fn parse_compression( } } +fn root_group(file: Arc) -> ReadGroup { + let root = file.superblock().root_group_address; + ReadGroup::new(file, String::new(), root) +} + /// Build and write the HDF5 file from accumulated write state. fn finalize_write(state: WriteState) -> PyResult<()> { - let mut builder = clawhdf5_rs::FileBuilder::new(); + crate::no_panic(|| { + let mut builder = clawhdf5_rs::FileBuilder::new(); - // Root attributes - let root_attrs = state.root_attrs.lock().unwrap_or_else(|e| e.into_inner()); - for (name, val) in root_attrs.iter() { - builder.set_attr(name, val.clone().into()); - } - drop(root_attrs); + // Root attributes + let root_attrs = state.root_attrs.lock().unwrap_or_else(|e| e.into_inner()); + for (name, val) in root_attrs.iter() { + builder.set_attr(name, val.clone().into()); + } + drop(root_attrs); - // Root datasets - for spec in &state.root_datasets { - let db = builder.create_dataset(&spec.name); - apply_dataset_spec(db, spec); - } + // Root datasets + for spec in &state.root_datasets { + let db = builder.create_dataset(&spec.name); + apply_dataset_spec(db, spec); + } - // Groups - for group_arc in &state.groups { - let guard = group_arc.lock().unwrap(); - finalize_write_group(&mut builder, &guard); - } + // Groups + for group_arc in &state.groups { + let guard = group_arc.lock().unwrap(); + finalize_write_group(&mut builder, &guard); + } - builder.write(&state.path).map_err(to_py_err)?; - Ok(()) + builder.write(&state.path).map_err(to_py_err)?; + Ok(()) + }) } #[cfg(test)] diff --git a/crates/clawhdf5-py/src/group.rs b/crates/clawhdf5-py/src/group.rs index 9f995f8..7585e57 100644 --- a/crates/clawhdf5-py/src/group.rs +++ b/crates/clawhdf5-py/src/group.rs @@ -1,13 +1,14 @@ //! PyGroup — navigable HDF5 group with read and write support. -use std::sync::{Arc, Mutex}; +use std::collections::HashMap; +use std::sync::{Arc, Mutex, OnceLock}; +use pyo3::exceptions::{PyIOError, PyKeyError, PyValueError}; use pyo3::prelude::*; use pyo3::types::PyList; use crate::attrs::PyAttrs; -use crate::dataset::PyDataset; -use crate::{DatasetSpec, OwnedAttrValue, apply_dataset_spec, extract_numpy_data, to_py_err}; +use crate::{DatasetSpec, OwnedAttrValue, apply_dataset_spec, extract_numpy_data, node}; /// Shared state for a group being written. pub(crate) struct WriteGroupState { @@ -18,32 +19,24 @@ pub(crate) struct WriteGroupState { /// An HDF5 group. /// -/// In read mode, provides `__getitem__` navigation and child listing. -/// In write mode, supports `create_dataset` and `create_group` and -/// attribute setting. -/// -/// ```python -/// grp = f['group_name'] -/// grp.keys() -/// ds = grp['dataset'] -/// ``` +/// In read mode it behaves like an h5py group: `grp['name']`, +/// `grp['sub/path']` and `grp['/absolute/path']`, `keys()`, `values()`, +/// `items()`, iteration, `len()`, `in`, `get()`, `name` and `attrs`. +/// In write mode, supports `create_dataset` and attribute setting. #[pyclass(name = "Group")] pub struct PyGroup { inner: GroupInner, } enum GroupInner { - Read { - file: Arc, - path: String, - }, + Read(ReadGroup), Write(Arc>), } impl PyGroup { - pub(crate) fn from_read(file: Arc, path: String) -> Self { + pub(crate) fn from_read(file: Arc, path: String, addr: u64) -> Self { Self { - inner: GroupInner::Read { file, path }, + inner: GroupInner::Read(ReadGroup::new(file, path, addr)), } } @@ -52,63 +45,172 @@ impl PyGroup { inner: GroupInner::Write(state), } } + + fn read_group(&self, what: &str) -> PyResult<&ReadGroup> { + match &self.inner { + GroupInner::Read(g) => Ok(g), + GroupInner::Write(_) => Err(PyIOError::new_err(format!( + "cannot {what} a group opened for writing" + ))), + } + } +} + +/// A group in a file opened for reading (a file is its root group, as in +/// h5py). It keeps its own address and, once listed, its links, so looking +/// up a child neither resolves the path from the root nor scans the group's +/// links again: visiting every member of a large group is linear, not +/// quadratic. +pub(crate) struct ReadGroup { + pub file: Arc, + pub path: String, + pub addr: u64, + /// Link name -> object address (soft links resolved), filled on first use. + links: OnceLock>, + /// Names of the datasets and subgroups, sorted (h5py's order). + members: OnceLock>, +} + +impl ReadGroup { + pub(crate) fn new(file: Arc, path: String, addr: u64) -> Self { + Self { + file, + path, + addr, + links: OnceLock::new(), + members: OnceLock::new(), + } + } + + fn links(&self) -> PyResult<&HashMap> { + if let Some(links) = self.links.get() { + return Ok(links); + } + let entries = crate::no_panic(|| { + clawhdf5_format::group_v2::resolve_group_children( + self.file.as_bytes(), + self.file.superblock(), + self.addr, + ) + .map_err(|e| PyValueError::new_err(format!("{}: {e}", node::name(&self.path)))) + })?; + let map = entries + .into_iter() + .map(|e| (e.name, e.object_header_address)) + .collect(); + Ok(self.links.get_or_init(|| map)) + } + + /// The path and address of `key` (a name, a relative or an absolute path). + fn locate(&self, key: &str) -> PyResult<(String, u64)> { + let path = node::join(&self.path, key); + let rel = if self.path.is_empty() { + Some(path.as_str()) + } else if path == self.path { + Some("") + } else { + path.strip_prefix(self.path.as_str()) + .and_then(|r| r.strip_prefix('/')) + }; + let addr = match rel { + // A direct child: the link table, when it has the name. + Some(name) if !name.is_empty() && !name.contains('/') => { + match self.links()?.get(name) { + Some(&a) => a, + None => node::resolve_from(&self.file, self.addr, name, &path)?, + } + } + Some(rel) => node::resolve_from(&self.file, self.addr, rel, &path)?, + None => node::address(&self.file, &path)?, + }; + Ok((path, addr)) + } + + /// `group[key]`. + pub(crate) fn get_item(&self, py: Python<'_>, key: &str) -> PyResult> { + let (path, addr) = self.locate(key)?; + node::open(py, &self.file, path, addr) + } + + /// `group.get(key, default)`. + pub(crate) fn get( + &self, + py: Python<'_>, + key: &str, + default: Option>, + ) -> PyResult> { + match self.get_item(py, key) { + Err(e) if e.is_instance_of::(py) => { + Ok(default.unwrap_or_else(|| py.None())) + } + other => other, + } + } + + /// Names of the group's datasets and subgroups, sorted (h5py's order). + pub(crate) fn member_names(&self) -> PyResult<&[String]> { + if let Some(m) = self.members.get() { + return Ok(m); + } + let mut names = Vec::new(); + for (name, &addr) in self.links()? { + let hdr = node::header_at(&self.file, addr, &node::join(&self.path, name))?; + if matches!( + node::kind(&hdr), + Some(node::Kind::Dataset | node::Kind::Group) + ) { + names.push(name.clone()); + } + } + names.sort_by(|a, b| a.as_bytes().cmp(b.as_bytes())); + Ok(self.members.get_or_init(|| names)) + } + + pub(crate) fn contains(&self, key: &str) -> bool { + self.locate(key) + .and_then(|(path, addr)| node::header_at(&self.file, addr, &path)) + .ok() + .and_then(|h| node::kind(&h)) + .is_some_and(|k| k != node::Kind::Datatype) + } + + pub(crate) fn values(&self, py: Python<'_>) -> PyResult>> { + self.member_names()? + .iter() + .map(|n| self.get_item(py, n)) + .collect() + } + + pub(crate) fn items(&self, py: Python<'_>) -> PyResult)>> { + self.member_names()? + .iter() + .map(|n| Ok((n.clone(), self.get_item(py, n)?))) + .collect() + } + + pub(crate) fn attrs(&self) -> PyResult { + PyAttrs::read(Arc::clone(&self.file), self.addr, &self.path) + } } #[pymethods] impl PyGroup { /// Get a child object (dataset or subgroup) by name or path. fn __getitem__(&self, py: Python<'_>, key: &str) -> PyResult> { - match &self.inner { - GroupInner::Read { file, path } => { - let full_path = if path.is_empty() { - key.to_string() - } else { - format!("{path}/{key}") - }; - // Try dataset first - match file.dataset(&full_path) { - Ok(_) => { - let ds = PyDataset::new(Arc::clone(file), full_path)?; - Ok(ds.into_pyobject(py)?.into_any().unbind()) - } - Err(clawhdf5_rs::Error::NotADataset(_)) => { - let grp = PyGroup::from_read(Arc::clone(file), full_path); - Ok(grp.into_pyobject(py)?.into_any().unbind()) - } - Err(e) => { - // Could be a group without a DataLayout message - match file.group(&full_path) { - Ok(_) => { - let grp = PyGroup::from_read(Arc::clone(file), full_path); - Ok(grp.into_pyobject(py)?.into_any().unbind()) - } - Err(_) => Err(PyErr::new::(format!( - "{key}: {e}" - ))), - } - } - } - } - GroupInner::Write(_) => Err(PyErr::new::( - "cannot read children from a group opened for writing", - )), - } + self.read_group("read children from")?.get_item(py, key) + } + + /// `group.get(key, default=None)`. + #[pyo3(signature = (key, default=None))] + fn get(&self, py: Python<'_>, key: &str, default: Option>) -> PyResult> { + self.read_group("read children from")?.get(py, key, default) } /// List the names of all children (datasets and subgroups). fn keys(&self, py: Python<'_>) -> PyResult> { match &self.inner { - GroupInner::Read { file, path } => { - let group = if path.is_empty() { - file.root() - } else { - file.group(path).map_err(to_py_err)? - }; - let mut names = group.datasets().map_err(to_py_err)?; - let groups = group.groups().map_err(to_py_err)?; - names.extend(groups); - names.sort(); - let list = PyList::new(py, &names)?; + GroupInner::Read(g) => { + let list = PyList::new(py, g.member_names()?)?; Ok(list.into_any().unbind()) } GroupInner::Write(state) => { @@ -120,6 +222,36 @@ impl PyGroup { } } + fn values(&self, py: Python<'_>) -> PyResult> { + let g = self.read_group("read children from")?; + Ok(PyList::new(py, g.values(py)?)?.into_any().unbind()) + } + + fn items(&self, py: Python<'_>) -> PyResult> { + let g = self.read_group("read children from")?; + Ok(PyList::new(py, g.items(py)?)?.into_any().unbind()) + } + + fn __iter__(&self, py: Python<'_>) -> PyResult> { + self.keys(py)?.call_method0(py, "__iter__") + } + + fn __len__(&self) -> PyResult { + match &self.inner { + GroupInner::Read(g) => Ok(g.member_names()?.len()), + GroupInner::Write(state) => Ok(state.lock().unwrap().datasets.len()), + } + } + + /// The group's full name, e.g. `/sensors`. + #[getter] + fn name(&self) -> String { + match &self.inner { + GroupInner::Read(g) => node::name(&g.path), + GroupInner::Write(state) => node::name(&state.lock().unwrap().name), + } + } + /// Create a dataset inside this group (write mode only). /// /// Parameters: @@ -161,7 +293,7 @@ impl PyGroup { state.lock().unwrap().datasets.push(spec); Ok(()) } - GroupInner::Read { .. } => Err(PyErr::new::( + GroupInner::Read { .. } => Err(PyIOError::new_err( "cannot create datasets on a read-only group", )), } @@ -171,15 +303,7 @@ impl PyGroup { #[getter] fn attrs(&self) -> PyResult { match &self.inner { - GroupInner::Read { file, path } => { - let group = if path.is_empty() { - file.root() - } else { - file.group(path).map_err(to_py_err)? - }; - let map = group.attrs().map_err(to_py_err)?; - Ok(PyAttrs::from_read(map)) - } + GroupInner::Read(g) => g.attrs(), GroupInner::Write(state) => { let store = Arc::clone(&state.lock().unwrap().attrs); Ok(PyAttrs::from_write(store)) @@ -189,12 +313,9 @@ impl PyGroup { fn __repr__(&self) -> String { match &self.inner { - GroupInner::Read { path, .. } => { - if path.is_empty() { - "".to_string() - } else { - format!("") - } + GroupInner::Read(g) => { + let n = g.member_names().map_or(0, |m| m.len()); + format!("", node::name(&g.path)) } GroupInner::Write(state) => { let name = &state.lock().unwrap().name; @@ -205,14 +326,7 @@ impl PyGroup { fn __contains__(&self, key: &str) -> PyResult { match &self.inner { - GroupInner::Read { file, path } => { - let full_path = if path.is_empty() { - key.to_string() - } else { - format!("{path}/{key}") - }; - Ok(file.dataset(&full_path).is_ok() || file.group(&full_path).is_ok()) - } + GroupInner::Read(g) => Ok(g.contains(key)), GroupInner::Write(state) => { let guard = state.lock().unwrap(); Ok(guard.datasets.iter().any(|d| d.name == key)) @@ -244,26 +358,28 @@ mod tests { use super::*; #[test] - fn read_group_construction() { + fn member_names_are_sorted() { let mut b = clawhdf5_rs::FileBuilder::new(); - let mut g = b.create_group("grp"); + b.create_dataset("zeta").with_f64_data(&[1.0]); + b.create_dataset("alpha").with_f64_data(&[1.0]); + let mut g = b.create_group("mid"); g.create_dataset("x").with_f64_data(&[1.0]); let finished = g.finish(); b.add_group(finished); let bytes = b.finish().unwrap(); let file = Arc::new(clawhdf5_rs::File::from_bytes(bytes).unwrap()); - let _grp = PyGroup::from_read(file, "grp".into()); - } - - #[test] - fn write_group_state() { - let state = WriteGroupState { - name: "test".into(), - datasets: vec![], - attrs: Arc::new(Mutex::new(vec![])), - }; - let arc = Arc::new(Mutex::new(state)); - let _grp = PyGroup::from_write(arc); + let root = file.superblock().root_group_address; + let top = ReadGroup::new(Arc::clone(&file), String::new(), root); + assert_eq!(top.member_names().unwrap(), ["alpha", "mid", "zeta"]); + let (path, addr) = top.locate("mid").unwrap(); + assert_eq!(path, "mid"); + let mid = ReadGroup::new(Arc::clone(&file), path, addr); + assert_eq!(mid.member_names().unwrap(), ["x"]); + assert!(top.contains("mid/x")); + assert!(mid.contains("/alpha")); + assert!(mid.contains("x") && mid.contains("./x")); + assert!(!top.contains("nope")); + assert!(!mid.contains("alpha")); } #[test] diff --git a/crates/clawhdf5-py/src/lib.rs b/crates/clawhdf5-py/src/lib.rs index 5721964..e3b5619 100644 --- a/crates/clawhdf5-py/src/lib.rs +++ b/crates/clawhdf5-py/src/lib.rs @@ -10,9 +10,12 @@ //! ``` mod attrs; +mod convert; mod dataset; mod file; mod group; +mod node; +mod select; use pyo3::prelude::*; @@ -21,6 +24,42 @@ pub(crate) use dataset::PyDataset; pub(crate) use file::PyFile; pub(crate) use group::PyGroup; +pyo3::create_exception!( + clawhdf5, + InternalError, + pyo3::exceptions::PyRuntimeError, + "A bug in clawhdf5 met while reading or writing a file (a Rust panic, \ + caught). Derived from RuntimeError, so `except Exception` handles it." +); + +/// The text of a caught panic. +pub(crate) fn panic_text(payload: &(dyn std::any::Any + Send)) -> String { + payload + .downcast_ref::<&str>() + .map(|s| (*s).to_string()) + .or_else(|| payload.downcast_ref::().cloned()) + .unwrap_or_else(|| "unknown panic".to_string()) +} + +/// Run `f`, turning a panic in the library into [`InternalError`] instead of +/// PyO3's `PanicException` (a `BaseException`, which `except Exception` +/// does not catch). Wraps every call into the library. +pub(crate) fn no_panic(f: impl FnOnce() -> PyResult) -> PyResult { + std::panic::catch_unwind(std::panic::AssertUnwindSafe(f)).unwrap_or_else(|p| { + Err(InternalError::new_err(format!( + "clawhdf5 internal error (please report it): {}", + panic_text(&*p) + ))) + }) +} + +/// A test hook: panics inside [`no_panic`], so the tests can check that a +/// library panic reaches Python as an ordinary exception. +#[pyfunction] +fn _panic_for_test() -> PyResult<()> { + no_panic(|| panic!("deliberate panic for the test suite")) +} + /// Convert a `clawhdf5_rs::Error` into a `PyErr`. /// /// Maps different error variants to more specific Python exception types: @@ -46,6 +85,54 @@ pub(crate) fn to_py_err(e: clawhdf5_rs::Error) -> PyErr { } } +/// The value of a dataset or attribute with a null dataspace: a type but no +/// data. Mirrors `h5py.Empty`. +#[pyclass(name = "Empty", frozen)] +pub struct PyEmpty { + dtype: Py, +} + +impl PyEmpty { + pub(crate) fn new(dtype: Py) -> Self { + Self { dtype } + } +} + +#[pymethods] +impl PyEmpty { + #[new] + fn py_new(py: Python<'_>, dtype: &Bound<'_, PyAny>) -> PyResult { + let dtype = py.import("numpy")?.getattr("dtype")?.call1((dtype,))?; + Ok(Self::new(dtype.unbind())) + } + + #[getter] + fn dtype(&self, py: Python<'_>) -> Py { + self.dtype.clone_ref(py) + } + + #[getter] + fn shape(&self, py: Python<'_>) -> Py { + py.None() + } + + #[getter] + fn size(&self, py: Python<'_>) -> Py { + py.None() + } + + fn __eq__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult { + match other.cast::() { + Ok(o) => self.dtype.bind(py).eq(o.get().dtype.bind(py)), + Err(_) => Ok(false), + } + } + + fn __repr__(&self, py: Python<'_>) -> PyResult { + Ok(format!("Empty(dtype={})", self.dtype.bind(py).repr()?)) + } +} + /// The data payload for a dataset being written. #[derive(Clone)] pub(crate) enum DatasetData { @@ -219,10 +306,14 @@ pub(crate) fn extract_numpy_data( /// The clawhdf5 Python module. #[pymodule] fn clawhdf5(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add("__version__", env!("CARGO_PKG_VERSION"))?; m.add_class::()?; m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; + m.add("InternalError", m.py().get_type::())?; + m.add_function(wrap_pyfunction!(_panic_for_test, m)?)?; Ok(()) } @@ -232,9 +323,9 @@ mod tests { #[test] fn owned_attr_value_roundtrip() { - let val = OwnedAttrValue::F64(3.14); + let val = OwnedAttrValue::F64(2.5); let attr: clawhdf5_rs::AttrValue = val.into(); - assert!(matches!(attr, clawhdf5_rs::AttrValue::F64(v) if (v - 3.14).abs() < 1e-10)); + assert!(matches!(attr, clawhdf5_rs::AttrValue::F64(v) if (v - 2.5).abs() < 1e-10)); } #[test] diff --git a/crates/clawhdf5-py/src/node.rs b/crates/clawhdf5-py/src/node.rs new file mode 100644 index 0000000..fa9ede7 --- /dev/null +++ b/crates/clawhdf5-py/src/node.rs @@ -0,0 +1,215 @@ +//! Resolving paths to objects in a file opened for reading. + +use std::sync::Arc; + +use clawhdf5_format::attribute::AttributeMessage; +use clawhdf5_format::dataspace::{Dataspace, DataspaceType}; +use clawhdf5_format::message_type::MessageType; +use clawhdf5_format::object_header::ObjectHeader; +use pyo3::exceptions::{PyKeyError, PyTypeError, PyValueError}; +use pyo3::prelude::*; + +use crate::dataset::PyDataset; +use crate::group::PyGroup; + +/// Join `key` onto the group path `base` the way h5py does: an absolute key +/// starts from the root, a relative one from `base`. Paths are kept without +/// a leading `/`; the root is `""`. +pub(crate) fn join(base: &str, key: &str) -> String { + let parts = if key.starts_with('/') { + key.split('/').collect::>() + } else { + base.split('/').chain(key.split('/')).collect() + }; + parts + .into_iter() + .filter(|p| !p.is_empty() && *p != ".") + .collect::>() + .join("/") +} + +/// The HDF5 name (`/a/b`) of a path. +pub(crate) fn name(path: &str) -> String { + format!("/{path}") +} + +/// The address of the object at `path`, resolved from the root group. +pub(crate) fn address(file: &clawhdf5_rs::File, path: &str) -> PyResult { + resolve_from(file, file.superblock().root_group_address, path, path) +} + +/// The address of `rel` resolved from the group at `group` (`full` is the +/// resulting path, for the error message). +pub(crate) fn resolve_from( + file: &clawhdf5_rs::File, + group: u64, + rel: &str, + full: &str, +) -> PyResult { + if rel.is_empty() { + return Ok(group); + } + crate::no_panic(|| { + clawhdf5_format::group_v2::resolve_path_from(file.as_bytes(), file.superblock(), group, rel) + .map_err(|e| { + PyKeyError::new_err(format!( + "Unable to open object (object '{}' doesn't exist): {e}", + name(full) + )) + }) + }) +} + +/// The object header at `addr` (the object at `path`). +pub(crate) fn header_at(file: &clawhdf5_rs::File, addr: u64, path: &str) -> PyResult { + crate::no_panic(|| { + let sb = file.superblock(); + let at = usize::try_from(addr) + .map_err(|_| PyValueError::new_err(format!("{}: address out of range", name(path))))?; + ObjectHeader::parse(file.as_bytes(), at, sb.offset_size, sb.length_size) + .map_err(|e| PyValueError::new_err(format!("{}: {e}", name(path)))) + }) +} + +/// What an object header describes. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Kind { + Dataset, + Group, + Datatype, +} + +pub(crate) fn kind(hdr: &ObjectHeader) -> Option { + let has = |t: MessageType| hdr.messages.iter().any(|m| m.msg_type == t); + if has(MessageType::DataLayout) { + Some(Kind::Dataset) + } else if has(MessageType::LinkInfo) + || has(MessageType::Link) + || has(MessageType::SymbolTable) + || has(MessageType::GroupInfo) + { + Some(Kind::Group) + } else if has(MessageType::Datatype) { + Some(Kind::Datatype) + } else { + None + } +} + +/// Open the object at `addr` (whose path is `path`) as a `Dataset` or +/// `Group`. Both keep the address, so later reads resolve nothing. +pub(crate) fn open( + py: Python<'_>, + file: &Arc, + path: String, + addr: u64, +) -> PyResult> { + let hdr = header_at(file, addr, &path)?; + match kind(&hdr) { + Some(Kind::Dataset) => Ok(PyDataset::open(py, Arc::clone(file), path, addr, &hdr)? + .into_pyobject(py)? + .into_any() + .unbind()), + Some(Kind::Group) => Ok(PyGroup::from_read(Arc::clone(file), path, addr) + .into_pyobject(py)? + .into_any() + .unbind()), + Some(Kind::Datatype) => Err(PyTypeError::new_err(format!( + "{}: committed (named) datatypes are not supported by clawhdf5", + name(&path) + ))), + None => Err(PyValueError::new_err(format!( + "{}: not a dataset, group or datatype", + name(&path) + ))), + } +} + +/// The dataspace message of an object header. +pub(crate) fn dataspace(file: &clawhdf5_rs::File, hdr: &ObjectHeader) -> PyResult { + crate::no_panic(|| { + let sb = file.superblock(); + let msg = hdr + .messages + .iter() + .find(|m| m.msg_type == MessageType::Dataspace) + .ok_or_else(|| PyValueError::new_err("object has no dataspace message"))?; + let data = clawhdf5_format::shared_message::message_data( + file.as_bytes(), + msg, + sb.offset_size, + sb.length_size, + ) + .map_err(|e| PyValueError::new_err(e.to_string()))?; + Dataspace::parse(&data, sb.length_size).map_err(|e| PyValueError::new_err(e.to_string())) + }) +} + +/// The chunk shape of a chunked dataset (one entry per dataset dimension), +/// or `None` for other layouts or a layout message that does not parse. +pub(crate) fn chunk_shape( + file: &clawhdf5_rs::File, + hdr: &ObjectHeader, + rank: usize, +) -> Option> { + let sb = file.superblock(); + let msg = hdr + .messages + .iter() + .find(|m| m.msg_type == MessageType::DataLayout)?; + match clawhdf5_format::data_layout::DataLayout::parse(&msg.data, sb.offset_size, sb.length_size) + .ok()? + { + clawhdf5_format::data_layout::DataLayout::Chunked { + chunk_dimensions, .. + } if chunk_dimensions.len() >= rank => Some( + chunk_dimensions[..rank] + .iter() + .map(|&d| u64::from(d)) + .collect(), + ), + _ => None, + } +} + +pub(crate) fn is_null(space: &Dataspace) -> bool { + space.space_type == DataspaceType::Null +} + +/// The attributes of the object at `addr` (whose path is `path`), sorted by +/// name (h5py's order). Attributes whose messages cannot be parsed are left +/// out, as the facade's `attrs()` does. +pub(crate) fn attributes( + file: &clawhdf5_rs::File, + addr: u64, + path: &str, +) -> PyResult> { + let hdr = header_at(file, addr, path)?; + crate::no_panic(|| { + let sb = file.superblock(); + let (mut attrs, _errors) = clawhdf5_format::attribute::extract_attributes_tolerant( + file.as_bytes(), + &hdr, + sb.offset_size, + sb.length_size, + ) + .map_err(|e| PyValueError::new_err(format!("{}: {e}", name(path))))?; + attrs.sort_by(|a, b| a.name.as_bytes().cmp(b.name.as_bytes())); + Ok(attrs) + }) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn join_paths() { + assert_eq!(join("", "a"), "a"); + assert_eq!(join("a", "b/c"), "a/b/c"); + assert_eq!(join("a/b", "/x"), "x"); + assert_eq!(join("a", "/"), ""); + assert_eq!(join("", "/a//b/"), "a/b"); + assert_eq!(join("a", "./b"), "a/b"); + } +} diff --git a/crates/clawhdf5-py/src/select.rs b/crates/clawhdf5-py/src/select.rs new file mode 100644 index 0000000..7d11f3c --- /dev/null +++ b/crates/clawhdf5-py/src/select.rs @@ -0,0 +1,554 @@ +//! h5py-style indexing (`ds[1, 2:10:3, ...]`) mapped onto hyperslab +//! selections, so the library reads the selection rather than the whole +//! dataset (it still decodes everything for large selections; see the +//! facade's `Dataset::read_selection`). +//! +//! The rules and error messages follow h5py's `selections.py`: integers +//! (negative from the end) drop their axis, slices must have a positive +//! step, one `Ellipsis` fills the unmentioned axes, a single increasing list +//! of integers may index one axis, and strings name compound fields. +//! Everything else (`None`/`np.newaxis`, boolean masks, several index lists) +//! is refused with the error h5py gives. + +use clawhdf5_format::selection::Selection; +use pyo3::exceptions::{PyIndexError, PyTypeError, PyValueError}; +use pyo3::prelude::*; +use pyo3::types::{PyEllipsis, PySlice, PyString, PyTuple}; + +/// The selection along one axis. +#[derive(Clone, Debug, PartialEq)] +pub(crate) enum Axis { + /// A single index: the axis is dropped from the result. + Index(u64), + /// `start, start + step, ...`, `count` of them. + Slice { start: u64, step: u64, count: u64 }, + /// Increasing, distinct indices. + List(Vec), +} + +impl Axis { + fn len(&self) -> u64 { + match self { + Axis::Index(_) => 1, + Axis::Slice { count, .. } => *count, + Axis::List(v) => v.len() as u64, + } + } +} + +/// A parsed index expression. +#[derive(Clone, Debug, PartialEq)] +pub(crate) struct Plan { + /// One entry per dataset axis. + pub axes: Vec, + /// Compound field names to keep (empty: all). + pub fields: Vec, + /// For a scalar dataset: `ds[()]` gives a scalar, `ds[...]` a 0-d array. + /// For other datasets: every axis was an integer, so h5py gives a scalar. + pub scalar: bool, +} + +impl Plan { + /// The shape of the result. + pub fn out_shape(&self) -> Vec { + self.axes + .iter() + .filter(|a| !matches!(a, Axis::Index(_))) + .map(|a| a.len() as usize) + .collect() + } + + /// The shape of the result before the integer-indexed axes are dropped + /// (they have length 1 here): the shape of the joined reads. + pub fn read_shape(&self) -> Vec { + self.axes.iter().map(|a| a.len() as usize).collect() + } + + /// The axis indexed by a list, if any. + pub fn list_axis(&self) -> Option { + self.axes.iter().position(|a| matches!(a, Axis::List(_))) + } + + /// Whether the selection is empty. + pub fn is_empty(&self) -> bool { + self.axes.iter().any(|a| a.len() == 0) + } + + /// The hyperslab reads that make up this selection, each with the shape + /// of its block (index axes kept at length 1). More than one only when an + /// axis is indexed by a list; those are joined along `list_axis` + /// afterwards (`join_along`), after keeping each read's `pick` rows. + /// + /// A list is read in groups, each one hyperslab over a stretch of the + /// axis, not once per index: every read decodes the chunks it touches + /// (and lists the dataset's chunks), so a read per run of indices decoded + /// the same chunk again and again. For a chunked dataset (`chunk_len` is + /// the chunk's length along the list axis) a group ends only where a + /// whole chunk holds no selected index, so no chunk is decoded twice or + /// without need. Otherwise a group ends at a gap of more than + /// [`MAX_GAP_BYTES`] of unselected data. + pub fn reads( + &self, + dims: &[u64], + chunk_len: Option, + elem_size: usize, + ) -> (Vec, Option) { + let list_axis = self.list_axis(); + let groups: Vec<&[u64]> = match list_axis.map(|i| &self.axes[i]) { + Some(Axis::List(idx)) => { + let row_bytes = self.row_bytes(elem_size); + group_indices(idx, |last, next| match chunk_len { + Some(c) if c > 0 => next / c <= last / c + 1, + _ => (next - last - 1).saturating_mul(row_bytes) <= MAX_GAP_BYTES, + }) + } + _ => vec![&[]], + }; + let mut out = Vec::with_capacity(groups.len()); + for group in groups { + let (first, span) = match (group.first(), group.last()) { + (Some(&f), Some(&l)) => (f, l - f + 1), + _ => (0, 0), + }; + let pick = (span != group.len() as u64) + .then(|| group.iter().map(|&i| (i - first) as usize).collect()); + let mut start = Vec::with_capacity(dims.len()); + let mut stride = Vec::with_capacity(dims.len()); + let mut count = Vec::with_capacity(dims.len()); + for axis in &self.axes { + let (s, st, c) = match axis { + Axis::Index(i) => (*i, 1, 1), + Axis::Slice { start, step, count } => (*start, *step, *count), + Axis::List(_) => (first, 1, span), + }; + start.push(s); + // A stride only matters between blocks; keep it >= 1. + stride.push(if c <= 1 { 1 } else { st }); + count.push(c); + } + let shape: Vec = count.iter().map(|&c| c as usize).collect(); + let whole = start.iter().all(|&s| s == 0) + && stride.iter().all(|&s| s == 1) + && count.as_slice() == dims; + let sel = if whole { + Selection::All + } else { + let block = vec![1; dims.len()]; + Selection::Hyperslab { + start, + stride, + count, + block, + } + }; + out.push(Read { sel, shape, pick }); + } + (out, list_axis) + } + + /// Bytes of one step along the list axis within a read's bounding box. + fn row_bytes(&self, elem_size: usize) -> u64 { + self.axes + .iter() + .map(|a| match a { + Axis::Slice { step, count, .. } if *count > 0 => (count - 1) * step + 1, + _ => 1, + }) + .fold(elem_size as u64, u64::saturating_mul) + } +} + +/// Unselected data a read of a non-chunked dataset copies through rather +/// than start another read. +pub(crate) const MAX_GAP_BYTES: u64 = 64 * 1024; + +/// One hyperslab read of a selection. +#[derive(Clone, Debug, PartialEq)] +pub(crate) struct Read { + pub sel: Selection, + /// The block's shape (index axes at length 1). + pub shape: Vec, + /// For a list: the positions along the list axis, within the block, to + /// keep (`None`: all of them). + pub pick: Option>, +} + +/// Split increasing indices into groups; `joins(last, next)` says whether +/// `next` extends the group whose last index is `last`. +fn group_indices(idx: &[u64], joins: impl Fn(u64, u64) -> bool) -> Vec<&[u64]> { + let mut groups = Vec::new(); + let mut from = 0; + for k in 1..idx.len() { + if !joins(idx[k - 1], idx[k]) { + groups.push(&idx[from..k]); + from = k; + } + } + if from < idx.len() { + groups.push(&idx[from..]); + } + groups +} + +/// Keep the elements at positions `pick` along `axis` of a row-major block. +pub(crate) fn gather_along( + bytes: &[u8], + shape: &[usize], + axis: usize, + pick: &[usize], + elem_size: usize, +) -> Vec { + let outer: usize = shape[..axis].iter().product(); + let inner: usize = shape[axis + 1..].iter().product::() * elem_size; + let len = shape[axis]; + let mut out = Vec::with_capacity(outer * pick.len() * inner); + for o in 0..outer { + for &p in pick { + let at = (o * len + p) * inner; + out.extend_from_slice(&bytes[at..at + inner]); + } + } + out +} + +/// Join row-major blocks of `elem_size`-byte elements whose shapes differ +/// only along `axis` into one buffer, in order along that axis. Whole +/// elements are copied, so compound padding keeps the bytes that were read. +pub(crate) fn join_along( + blocks: &[(Vec, Vec)], + axis: usize, + elem_size: usize, +) -> Vec { + let Some((_, first)) = blocks.first() else { + return Vec::new(); + }; + let outer: usize = first[..axis].iter().product(); + let inner: usize = first[axis + 1..].iter().product::() * elem_size; + let total: usize = blocks.iter().map(|(_, s)| s[axis]).sum(); + let mut out = vec![0u8; outer * total * inner]; + let mut at = 0; + for (bytes, shape) in blocks { + let len = shape[axis] * inner; + for o in 0..outer { + let dst = (o * total) * inner + at; + out[dst..dst + len].copy_from_slice(&bytes[o * len..(o + 1) * len]); + } + at += len; + } + out +} + +/// Parse `key` for a dataset of shape `dims`. +pub(crate) fn parse(key: &Bound<'_, PyAny>, dims: &[u64]) -> PyResult { + let items: Vec> = match key.cast::() { + Ok(t) => t.iter().collect(), + Err(_) => vec![key.clone()], + }; + let mut fields = Vec::new(); + let mut args = Vec::new(); + for item in items { + if let Ok(s) = item.cast::() { + fields.push(s.to_str()?.to_owned()); + } else { + args.push(item); + } + } + + if args.iter().any(|a| a.is_none()) { + return Err(PyTypeError::new_err( + "Indexing with None (or np.newaxis) is not supported", + )); + } + + let rank = dims.len(); + if rank == 0 { + return match args.as_slice() { + [] => Ok(Plan { + axes: vec![], + fields, + scalar: true, + }), + [a] if a.is_instance_of::() => Ok(Plan { + axes: vec![], + fields, + scalar: false, + }), + _ => Err(PyValueError::new_err( + "Illegal slicing argument for scalar dataspace", + )), + }; + } + + // Expand the ellipsis (at most one) to full slices. + let n_ellipsis = args + .iter() + .filter(|a| a.is_instance_of::()) + .count(); + if n_ellipsis > 1 { + return Err(PyValueError::new_err("Only one ellipsis may be used.")); + } + let explicit = args.len() - n_ellipsis; + if explicit > rank { + return Err(PyValueError::new_err(format!( + "{explicit} indexing arguments for {rank} dimensions" + ))); + } + let py = key.py(); + let mut expanded: Vec>> = Vec::with_capacity(rank); + for a in args { + if a.is_instance_of::() { + for _ in 0..(rank - explicit) { + expanded.push(None); + } + } else { + expanded.push(Some(a)); + } + } + while expanded.len() < rank { + expanded.push(None); + } + + let mut axes = Vec::with_capacity(rank); + for (arg, &n) in expanded.iter().zip(dims) { + axes.push(match arg { + None => Axis::Slice { + start: 0, + step: 1, + count: n, + }, + Some(a) => parse_axis(py, a, n)?, + }); + } + if axes.iter().filter(|a| matches!(a, Axis::List(_))).count() > 1 { + return Err(PyTypeError::new_err( + "Only one indexing vector or array is currently allowed for fancy indexing", + )); + } + let scalar = axes.iter().all(|a| matches!(a, Axis::Index(_))); + Ok(Plan { + axes, + fields, + scalar, + }) +} + +fn parse_axis(py: Python<'_>, a: &Bound<'_, PyAny>, n: u64) -> PyResult { + if a.is_none() { + return Err(PyTypeError::new_err( + "Indexing with None (or np.newaxis) is not supported", + )); + } + if let Ok(s) = a.cast::() { + let n_isize = isize::try_from(n) + .map_err(|_| PyValueError::new_err("dimension too large to slice"))?; + let ind = s.indices(n_isize)?; + if ind.step < 1 { + return Err(PyValueError::new_err(format!( + "Step must be >= 1 (got {})", + ind.step + ))); + } + // `slicelength` is the number of elements selected, >= 0. + let count = ind.slicelength as u64; + let start = if count == 0 { 0 } else { ind.start as u64 }; + return Ok(Axis::Slice { + start, + step: ind.step as u64, + count, + }); + } + let np = py.import("numpy")?; + let is_bool = + a.is_instance_of::() || a.is_instance(&np.getattr("bool_")?)?; + let is_array_like = a.is_instance(&np.getattr("ndarray")?)? + || a.is_instance_of::() + || a.is_instance_of::(); + // A 0-d integer array (`ds[np.array(1)]`) is an integer index, as in h5py. + if a.is_instance(&np.getattr("ndarray")?)? && a.getattr("ndim")?.extract::()? == 0 { + let kind: String = a.getattr("dtype")?.getattr("kind")?.extract()?; + if kind == "i" || kind == "u" { + let i: i128 = a.call_method0("item")?.extract()?; + return Ok(Axis::Index(normalize(i, n)?)); + } + } + if !is_bool && !is_array_like && a.hasattr("__index__")? { + let i: i128 = a.call_method0("__index__")?.extract()?; + return Ok(Axis::Index(normalize(i, n)?)); + } + if is_array_like { + let arr = np.call_method1("asarray", (a,))?; + let kind: String = arr.getattr("dtype")?.getattr("kind")?.extract()?; + if kind == "b" { + return Err(PyTypeError::new_err( + "Boolean mask indexing is not supported by clawhdf5", + )); + } + let ndim: usize = arr.getattr("ndim")?.extract()?; + let size: usize = arr.getattr("size")?.extract()?; + if size > 0 && kind != "i" && kind != "u" { + return Err(PyTypeError::new_err( + "Indexing arrays must have integer dtypes", + )); + } + if ndim > 1 { + return Err(PyTypeError::new_err( + "Only 1-D integer lists or arrays can be used for fancy indexing", + )); + } + let vals: Vec = arr.call_method0("tolist")?.extract()?; + let mut idx = Vec::with_capacity(vals.len()); + for v in vals { + idx.push(normalize(v, n)?); + } + if idx.windows(2).any(|w| w[0] >= w[1]) { + return Err(PyTypeError::new_err( + "Indexing elements must be in increasing order", + )); + } + return Ok(Axis::List(idx)); + } + Err(PyTypeError::new_err(format!( + "Illegal index type for clawhdf5 datasets: {}", + a.get_type().name()? + ))) +} + +fn normalize(i: i128, n: u64) -> PyResult { + let n_i = i128::from(n); + let j = if i < 0 { i + n_i } else { i }; + if j < 0 || j >= n_i { + let hi = n_i - 1; + return Err(PyIndexError::new_err(format!( + "Index ({i}) out of range for (0-{hi})" + ))); + } + Ok(j as u64) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn indices_group_by_chunk() { + let chunked = |c: u64| move |last: u64, next: u64| next / c <= last / c + 1; + // Chunks of 10: 3, 5 and 15 are in neighbouring chunks; 42 skips two. + let idx = [3, 5, 15, 42, 43, 99]; + assert_eq!( + group_indices(&idx, chunked(10)), + vec![&[3, 5, 15][..], &[42, 43], &[99]] + ); + assert_eq!(group_indices(&[], chunked(10)), Vec::<&[u64]>::new()); + } + + #[test] + fn a_list_reads_once_per_group() { + let plan = Plan { + axes: vec![ + Axis::List(vec![0, 2, 3, 40]), + Axis::Slice { + start: 0, + step: 1, + count: 5, + }, + ], + fields: vec![], + scalar: false, + }; + // Chunks of 8 rows: rows 0-3 are one read, row 40 another. + let (reads, axis) = plan.reads(&[50, 5], Some(8), 4); + assert_eq!(axis, Some(0)); + assert_eq!(reads.len(), 2); + assert_eq!(reads[0].shape, vec![4, 5]); + assert_eq!(reads[0].pick, Some(vec![0, 2, 3])); + assert_eq!(reads[1].shape, vec![1, 5]); + assert_eq!(reads[1].pick, None); + // Not chunked: a gap under MAX_GAP_BYTES is read through. + let (reads, _) = plan.reads(&[50, 5], None, 4); + assert_eq!(reads.len(), 1); + assert_eq!(reads[0].pick, Some(vec![0, 2, 3, 40])); + } + + #[test] + fn gather_keeps_picked_rows() { + // A 2x3 block of 1-byte elements; keep columns 0 and 2. + let block = [1, 2, 3, 4, 5, 6]; + assert_eq!( + gather_along(&block, &[2, 3], 1, &[0, 2], 1), + vec![1, 3, 4, 6] + ); + } + + #[test] + fn blocks_join_along_the_list_axis() { + // Two 2x1 and 2x2 blocks of 1-byte elements, joined along axis 1. + let a = (vec![1, 2], vec![2, 1]); + let b = (vec![3, 4, 5, 6], vec![2, 2]); + assert_eq!(join_along(&[a, b], 1, 1), vec![1, 3, 4, 2, 5, 6]); + // Along axis 0 it is concatenation; 2-byte elements stay whole. + let a = (vec![1, 2, 3, 4], vec![1, 2]); + let b = (vec![5, 6, 7, 8], vec![1, 2]); + assert_eq!(join_along(&[a, b], 0, 2), vec![1, 2, 3, 4, 5, 6, 7, 8]); + } + + #[test] + fn full_selection_reads_everything() { + let plan = Plan { + axes: vec![ + Axis::Slice { + start: 0, + step: 1, + count: 4, + }, + Axis::Slice { + start: 0, + step: 1, + count: 3, + }, + ], + fields: vec![], + scalar: false, + }; + let (reads, list) = plan.reads(&[4, 3], None, 8); + assert_eq!(list, None); + assert_eq!( + reads, + vec![Read { + sel: Selection::All, + shape: vec![4, 3], + pick: None + }] + ); + } + + #[test] + fn index_and_step_map_to_a_hyperslab() { + let plan = Plan { + axes: vec![ + Axis::Index(2), + Axis::Slice { + start: 1, + step: 3, + count: 2, + }, + ], + fields: vec![], + scalar: false, + }; + let (reads, _) = plan.reads(&[4, 8], None, 8); + assert_eq!( + reads, + vec![Read { + sel: Selection::Hyperslab { + start: vec![2, 1], + stride: vec![1, 3], + count: vec![1, 2], + block: vec![1, 1], + }, + shape: vec![1, 2], + pick: None + }] + ); + assert_eq!(plan.out_shape(), vec![2]); + } +} diff --git a/crates/clawhdf5-py/tests/conftest.py b/crates/clawhdf5-py/tests/conftest.py new file mode 100644 index 0000000..967f1e0 --- /dev/null +++ b/crates/clawhdf5-py/tests/conftest.py @@ -0,0 +1,18 @@ +"""Shared fixtures for the clawhdf5 Python binding tests.""" + +import os + +import pytest + + +@pytest.fixture(scope="session") +def h5py(): + """h5py, or a skip — unless CLAWHDF5_REQUIRE_INTEROP=1, which makes a + missing h5py a failure (as for the Rust interop suites).""" + try: + import h5py as mod + except ImportError: + if os.environ.get("CLAWHDF5_REQUIRE_INTEROP") == "1": + pytest.fail("h5py is required (CLAWHDF5_REQUIRE_INTEROP=1) but not importable") + pytest.skip("h5py not installed") + return mod diff --git a/crates/clawhdf5-py/tests/test_read_vs_h5py.py b/crates/clawhdf5-py/tests/test_read_vs_h5py.py new file mode 100644 index 0000000..6a3b6cb --- /dev/null +++ b/crates/clawhdf5-py/tests/test_read_vs_h5py.py @@ -0,0 +1,653 @@ +"""Every read through clawhdf5 compared against h5py (libhdf5) on files h5py +writes: dtypes, shapes, values and the type of what comes back (array, +numpy scalar, bytes, str, Empty), for every datatype the bindings map and a +spread of index expressions; plus the errors h5py gives for the same keys.""" + +import threading +from concurrent.futures import ThreadPoolExecutor + +import numpy as np +import pytest + +import clawhdf5 + +# --------------------------------------------------------------------------- +# The generated file +# --------------------------------------------------------------------------- + +NUMERIC = [ + "i2", ">i4", ">i8", ">u2", ">u4", ">u8", + "f2", ">f4", ">f8", +] + + +def _values(dtype, shape, seed): + rng = np.random.default_rng(seed) + dt = np.dtype(dtype) + n = int(np.prod(shape)) + if dt.kind == "f": + return rng.standard_normal(n).astype(dt).reshape(shape) + info = np.iinfo(dt) + return rng.integers(info.min, info.max, size=n, dtype=dt.newbyteorder("=")).astype(dt).reshape(shape) + + +def _compound_dtype(): + return np.dtype( + { + "names": ["id", "pos", "label", "flag", "vec"], + "formats": ["f8", "S6", "u1", ("", "be_") + f.create_dataset(f"num/{name}_1d", data=_values(dt, (37,), i)) + f.create_dataset( + f"num/{name}_2d_gzip", + data=_values(dt, (13, 11), 100 + i), + chunks=(4, 5), + compression="gzip", + ) + f.create_dataset("num/f8_3d", data=_values("f8")) + fseqs = np.empty(3, dtype=object) + fseqs[:] = [np.linspace(0, 1, 5).astype(">f8"), np.array([2.5], dtype=">f8"), np.array([], dtype=">f8")] + f.create_dataset("vlen/f8_be", data=fseqs, dtype=vl_f) + + # Compounds. + cdt = _compound_dtype() + rec = np.zeros(10, dtype=cdt) + rec["id"] = np.arange(10) * 3 + rec["pos"] = np.linspace(-1, 1, 10) + rec["label"] = [f"n{i}".encode() for i in range(10)] + rec["flag"] = np.arange(10) % 2 + rec["vec"] = np.arange(30, dtype=" 1 else [0], (Ellipsis,)] + if n0 > 5: + keys += [[1, 2, 3, 5], [0, 4, 5]] + if len(shape) >= 2: + n1 = shape[1] + keys += [ + (0, 0), + (-1, -1), + (slice(1, 3), slice(None, None, 2)), + (Ellipsis, 1), + (1, Ellipsis), + (slice(None), [0, n1 - 1] if n1 > 1 else [0]), + (slice(None, None, 2), 1), + (slice(0, 2), slice(3, 1)), + (np.array(1) if n0 > 1 else np.array(0), slice(None)), + (slice(None), np.array(n1 - 1, dtype="u2")), + ] + if len(shape) >= 3: + keys += [(0, slice(None), -1), (slice(1, None, 2), 2, slice(None, None, 3)), (Ellipsis, 0, 0), (0, Ellipsis, 1)] + return keys + + +# h5py 3.16 (HDF5 2.0) returns the elements of a variable-length sequence of +# big-endian floats unswapped (0.25 comes back as 2.6e-319), so these are +# checked against the values written instead (test_vlen_big_endian). +H5PY_MISREADS = {"vlen/f8_be"} + +ERROR_KEYS_1D = [ + slice(None, None, -1), 10**6, -(10**6), None, (0, 0, 0, 0, 0), [3, 1], (Ellipsis, Ellipsis), 1.5, "nope", + np.array(1.0), np.array(True), np.array(10**6), [0, 0], [], (), Ellipsis, (0,), [-1], np.array([1, 2]), +] + + +# --------------------------------------------------------------------------- +# Tests +# --------------------------------------------------------------------------- + + +def test_every_dataset_matches_h5py(h5py, pair): + ours, theirs, _ = pair + checked = 0 + for name in _all_datasets(h5py, theirs): + if name.startswith("unsupported/") or name in H5PY_MISREADS or name == "cmp/with_vlen": + continue + t = theirs[name] + o = ours[name] + assert o.shape == t.shape, name + assert o.dtype == t.dtype, f"{name}: {o.dtype} vs {t.dtype}" + assert dict(o.dtype.metadata or {}) == dict(t.dtype.metadata or {}), name + assert o.ndim == t.ndim and o.size == t.size, name + assert o.maxshape == t.maxshape, name + assert o.name == t.name, name + if t.shape is None: + assert_same(o[()], t[()], name) + continue + for key in keys_for(t.shape): + what = f"{name}[{key!r}]" + try: + expected = t[key] + except Exception as e: # noqa: BLE001 - h5py refuses: so must we + with pytest.raises(type(e)): + o[key] + continue + assert_same(o[key], expected, what) + checked += 1 + assert checked > 500 + + +def test_vlen_big_endian(pair): + ours, _, _ = pair + ds = ours["vlen/f8_be"] + assert ds.dtype.metadata["vlen"] == np.dtype(">f8") + got = ds[()] + expected = [np.linspace(0, 1, 5), np.array([2.5]), np.array([])] + assert got.shape == (3,) + for g, e in zip(got, expected): + assert g.dtype == np.dtype(">f8") + np.testing.assert_array_equal(g, e) + np.testing.assert_array_equal(ds[1], [2.5]) + + +def test_errors_match_h5py(h5py, pair): + ours, theirs, _ = pair + for name in ["num/le_i4_1d", "num/be_f8_2d_gzip", "str/vlen", "num/scalar_f8"]: + for key in ERROR_KEYS_1D: + try: + expected = theirs[name][key] + except Exception as e: # noqa: BLE001 + with pytest.raises(type(e)): + ours[name][key] + else: + # h5py reads it, so must we (and the same values). + assert_same(ours[name][key], expected, f"{name}[{key!r}]") + + +def test_compound_fields_match_h5py(pair): + ours, theirs, _ = pair + for name in ["cmp/padded", "cmp/padded_chunked", "cmp/nested_2d"]: + t, o = theirs[name], ours[name] + for field in t.dtype.names: + assert_same(o[field], t[field], f"{name}[{field}]") + assert_same(o[field, 1:3], t[field, 1:3], f"{name}[{field}, 1:3]") + two = list(t.dtype.names[:2]) + expected = t[tuple(two)] + got = o[tuple(two)] + assert got.dtype.names == expected.dtype.names + for field in two: + np.testing.assert_array_equal(got[field], expected[field]) + with pytest.raises(ValueError): + ours["cmp/padded"]["no_such_field"] + with pytest.raises(ValueError): + ours["num/le_i4_1d"]["id"] + + +def test_numpy_asarray_and_len(pair): + ours, theirs, _ = pair + assert_same(np.asarray(ours["num/f8_3d"]), np.asarray(theirs["num/f8_3d"])) + assert len(ours["num/f8_3d"]) == len(theirs["num/f8_3d"]) + with pytest.raises(TypeError): + len(ours["num/scalar_f8"]) + + +def test_attributes_match_h5py(h5py, pair): + ours, theirs, _ = pair + for path in ["/", "num/le_f8_1d", "deep/er/est"]: + t = theirs[path].attrs + o = ours[path].attrs + assert list(o.keys()) == sorted(t.keys()), path + assert len(o) == len(t) + for k in t.keys(): + assert k in o + assert_same(o[k], t[k], f"{path}.attrs[{k}]") + assert o.get("missing", 5) == 5 + with pytest.raises(KeyError): + o["missing"] + assert [k for k, _ in o.items()] == list(o.keys()) + + +def test_groups_match_h5py(h5py, pair): + ours, theirs, _ = pair + assert list(ours.keys()) == list(theirs.keys()) + assert len(ours) == len(theirs) + for path in ["num", "deep", "deep/er", "deep/er/est", "empty_group"]: + o, t = ours[path], theirs[path] + assert isinstance(o, clawhdf5.Group) + assert list(o.keys()) == list(t.keys()), path + assert list(o) == list(t), path + assert len(o) == len(t), path + assert o.name == t.name + g = ours["deep/er"] + assert isinstance(g["est"], clawhdf5.Group) + assert isinstance(g["est/leaf"], clawhdf5.Dataset) + assert g["/deep/er/est/leaf"].name == "/deep/er/est/leaf" + assert "est/leaf" in g and "/num" in g and "nope" not in g + assert ours["/"].name == "/" + assert "deep/er/est/leaf" in ours + assert ours.get("nope") is None + with pytest.raises(KeyError): + ours["deep/nope"] + names = [k for k, v in ours["deep/er/est"].items()] + assert names == ["leaf"] + assert isinstance(ours["deep/er/est"].values()[0], clawhdf5.Dataset) + + +def test_unsupported_types_are_errors_not_data(pair): + ours, _, _ = pair + ds = ours["unsupported/refs"] # opening works + with pytest.raises(TypeError): + ds.dtype + with pytest.raises(TypeError): + ds[()] + with pytest.raises(TypeError): + ours["cmp/with_vlen"][()] + + +def test_boolean_masks_are_refused(pair): + ours, _, _ = pair + with pytest.raises(TypeError): + ours["num/le_i4_1d"][np.ones(37, dtype=bool)] + + +def test_reads_hand_numpy_the_rust_buffer(pair): + """A fixed-size read is a view over the buffer the library filled, not a + copy of it.""" + ours, _, _ = pair + arr = ours["num/le_f8_2d_gzip"][2:9, 1:4] + assert not arr.flags.owndata + assert arr.base is not None + assert arr.flags.aligned and arr.flags.c_contiguous + + +def test_only_the_selected_chunks_are_read(h5py, tmp_path): + """Damage one chunk: a selection that avoids it still reads, one that + touches it fails. Reading everything and slicing afterwards (what the + bindings did before) failed both. (The library reads the whole dataset + anyway when a selection's bounding box covers more than half of it, so + the selections here stay below that.)""" + path = str(tmp_path / "damaged.h5") + data = np.arange(1000, dtype=" 0.03, f"a read took only {one_read:.3f} s; too short to measure" + + old = sys.getswitchinterval() + sys.setswitchinterval(0.001) + stop = threading.Event() + gaps = [] + + def spin(): + last = time.perf_counter() + worst = 0.0 + while not stop.is_set(): + now = time.perf_counter() + worst = max(worst, now - last) + last = now + gaps.append(worst) + + try: + t = threading.Thread(target=spin) + t.start() + time.sleep(0.01) + t0 = time.perf_counter() + for _ in range(2): + ds[key] + reading = time.perf_counter() - t0 + stop.set() + t.join() + finally: + sys.setswitchinterval(old) + np.testing.assert_array_equal(ds[key], data[:1000]) + # Held, the spinner would stall for about one read. + assert gaps[0] < one_read / 3, f"spinner stalled {gaps[0]:.3f} s during reads of {one_read:.3f} s ({reading:.3f} s)" + + +def _v4_index_fixture(h5py, path): + """One 2-D dataset per v4 chunk index (HDF5 1.10+ layout, libver='latest').""" + data = (np.arange(37 * 23, dtype="= 2 def test_attrs_contains(sample_read_file): - with rustyhdf5.File(sample_read_file, "r") as f: + with clawhdf5.File(sample_read_file, "r") as f: assert "version" in f.attrs assert "nonexistent" not in f.attrs def test_attrs_keys(sample_read_file): - with rustyhdf5.File(sample_read_file, "r") as f: + with clawhdf5.File(sample_read_file, "r") as f: keys = f.attrs.keys() assert "version" in keys assert "description" in keys @@ -114,7 +114,7 @@ def test_attrs_keys(sample_read_file): def test_read_group_keys(grouped_read_file): - with rustyhdf5.File(grouped_read_file, "r") as f: + with clawhdf5.File(grouped_read_file, "r") as f: keys = f.keys() assert "sensors" in keys assert "metadata" in keys @@ -122,7 +122,7 @@ def test_read_group_keys(grouped_read_file): def test_read_group_dataset(grouped_read_file): - with rustyhdf5.File(grouped_read_file, "r") as f: + with clawhdf5.File(grouped_read_file, "r") as f: grp = f["sensors"] ds = grp["temperature"] data = ds[:] @@ -130,14 +130,14 @@ def test_read_group_dataset(grouped_read_file): def test_read_group_attrs(grouped_read_file): - with rustyhdf5.File(grouped_read_file, "r") as f: + with clawhdf5.File(grouped_read_file, "r") as f: grp = f["sensors"] - assert grp.attrs["location"] == "lab" + assert grp.attrs["location"] == b"lab" def test_nested_path_access(grouped_read_file): """Test f['group/dataset'] path navigation.""" - with rustyhdf5.File(grouped_read_file, "r") as f: + with clawhdf5.File(grouped_read_file, "r") as f: ds = f["sensors/temperature"] data = ds[:] np.testing.assert_array_almost_equal(data, [22.5, 23.1, 21.8]) @@ -149,7 +149,7 @@ def test_nested_path_access(grouped_read_file): def test_context_manager(sample_read_file): - with rustyhdf5.File(sample_read_file, "r") as f: + with clawhdf5.File(sample_read_file, "r") as f: data = f["temperatures"][:] np.testing.assert_array_almost_equal(data, [22.5, 23.1, 21.8]) # File should be closed after with block @@ -162,30 +162,30 @@ def test_context_manager(sample_read_file): def test_write_simple(tmp_h5): - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("data", data=np.array([1.0, 2.0, 3.0])) # Verify by reading back - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: data = f["data"][:] np.testing.assert_array_almost_equal(data, [1.0, 2.0, 3.0]) def test_write_with_attrs(tmp_h5): - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("values", data=np.array([10, 20], dtype=np.int32)) f.attrs["author"] = "test" f.attrs["count"] = 42 - with rustyhdf5.File(tmp_h5, "r") as f: - assert f.attrs["author"] == "test" + with clawhdf5.File(tmp_h5, "r") as f: + assert f.attrs["author"] == b"test" assert f.attrs["count"] == 42 def test_write_with_group(tmp_h5): - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: grp = f.create_group("experiment") grp.create_dataset("results", data=np.array([3.14, 2.72])) grp.attrs["version"] = 1 - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: ds = f["experiment/results"] np.testing.assert_array_almost_equal(ds[:], [3.14, 2.72]) grp = f["experiment"] @@ -199,9 +199,9 @@ def test_write_with_group(tmp_h5): def test_roundtrip_float64(tmp_h5): original = np.array([1.1, 2.2, 3.3], dtype=np.float64) - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("data", data=original) - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: result = f["data"][:] np.testing.assert_array_almost_equal(result, original) assert result.dtype == np.float64 @@ -209,9 +209,9 @@ def test_roundtrip_float64(tmp_h5): def test_roundtrip_float32(tmp_h5): original = np.array([1.5, 2.5, 3.5], dtype=np.float32) - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("data", data=original) - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: result = f["data"][:] np.testing.assert_array_almost_equal(result, original) assert result.dtype == np.float32 @@ -219,9 +219,9 @@ def test_roundtrip_float32(tmp_h5): def test_roundtrip_int32(tmp_h5): original = np.array([-10, 0, 10, 100], dtype=np.int32) - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("data", data=original) - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: result = f["data"][:] np.testing.assert_array_equal(result, original) assert result.dtype == np.int32 @@ -229,9 +229,9 @@ def test_roundtrip_int32(tmp_h5): def test_roundtrip_int64(tmp_h5): original = np.array([-1, 0, 1, 2**40], dtype=np.int64) - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("data", data=original) - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: result = f["data"][:] np.testing.assert_array_equal(result, original) assert result.dtype == np.int64 @@ -239,9 +239,9 @@ def test_roundtrip_int64(tmp_h5): def test_roundtrip_uint8(tmp_h5): original = np.array([0, 127, 255], dtype=np.uint8) - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("data", data=original) - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: result = f["data"][:] np.testing.assert_array_equal(result, original) assert result.dtype == np.uint8 @@ -254,7 +254,7 @@ def test_roundtrip_uint8(tmp_h5): def test_chunked_gzip(tmp_h5): original = np.arange(100, dtype=np.float64) - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset( "compressed", data=original, @@ -262,7 +262,7 @@ def test_chunked_gzip(tmp_h5): compression="gzip", compression_opts=6, ) - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: result = f["compressed"][:] np.testing.assert_array_equal(result, original) @@ -276,7 +276,7 @@ def test_h5py_can_read_our_file(tmp_h5): """Verify that h5py can read files we create.""" import h5py - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("values", data=np.array([1.0, 2.0, 3.0])) f.attrs["meta"] = "hello" with h5py.File(tmp_h5, "r") as f: @@ -292,7 +292,7 @@ def test_we_can_read_h5py_file(tmp_h5): with h5py.File(tmp_h5, "w") as f: f.create_dataset("data", data=np.array([10.0, 20.0, 30.0])) f.attrs["version"] = 2 - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: data = f["data"][:] np.testing.assert_array_equal(data, [10.0, 20.0, 30.0]) assert f.attrs["version"] == 2 @@ -305,9 +305,9 @@ def test_we_can_read_h5py_file(tmp_h5): def test_2d_array_roundtrip(tmp_h5): original = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], dtype=np.float64) - with rustyhdf5.File(tmp_h5, "w") as f: + with clawhdf5.File(tmp_h5, "w") as f: f.create_dataset("matrix", data=original) - with rustyhdf5.File(tmp_h5, "r") as f: + with clawhdf5.File(tmp_h5, "r") as f: ds = f["matrix"] assert ds.shape == (2, 3) result = ds[:] @@ -321,15 +321,15 @@ def test_2d_array_roundtrip(tmp_h5): def test_open_nonexistent_file(): with pytest.raises(OSError): - rustyhdf5.File("/nonexistent/path.h5", "r") + clawhdf5.File("/nonexistent/path.h5", "r") def test_invalid_mode(tmp_h5): with pytest.raises(ValueError): - rustyhdf5.File(tmp_h5, "x") + clawhdf5.File(tmp_h5, "x") def test_key_error_on_missing_dataset(sample_read_file): - with rustyhdf5.File(sample_read_file, "r") as f: + with clawhdf5.File(sample_read_file, "r") as f: with pytest.raises(KeyError): f["nonexistent"] diff --git a/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index 635735d..90faf11 100644 --- a/crates/clawhdf5/src/reader.rs +++ b/crates/clawhdf5/src/reader.rs @@ -176,6 +176,22 @@ impl File { }) } + /// A `Dataset` handle for the object header at `address` (an address + /// from a group listing, or one kept from an earlier lookup), without + /// resolving a path. Resolving a path walks every group on it, which in + /// a large group costs a scan of its links; keep the address instead to + /// open the same dataset repeatedly. + pub fn dataset_at(&self, address: u64) -> Result, Error> { + let hdr = self.parse_header(address)?; + if !has_message(&hdr, MessageType::DataLayout) { + return Err(Error::NotADataset(format!("object at address {address}"))); + } + Ok(Dataset { + file: self, + header: hdr, + }) + } + /// Resolve a path and return a `Group` handle. /// /// The path uses `/` separators (e.g., `"sensors"`). @@ -606,8 +622,16 @@ impl<'f> Dataset<'f> { /// Read selected elements as raw bytes. /// - /// Only the elements matching the [`clawhdf5_format::selection::Selection`] are returned. For chunked - /// datasets, only intersecting chunks are decompressed. + /// Only the elements matching the [`clawhdf5_format::selection::Selection`] are returned. + /// + /// What is read to get them: when the selection's bounding box covers at + /// most half the dataset, only that box — the chunks overlapping it, or + /// the rows of a contiguous dataset. The whole dataset is decoded instead + /// when the box covers more than half (a strided selection spanning the + /// dataset does), for compact and virtual layouts, for a dataset with no + /// storage, and for a chunked dataset with a non-default fill value. + /// [`Selection::All`](clawhdf5_format::selection::Selection::All) goes + /// through the file's chunk cache; other selections do not. pub fn read_selection( &self, selection: &clawhdf5_format::selection::Selection, diff --git a/crates/clawhdf5/tests/integration_tests.rs b/crates/clawhdf5/tests/integration_tests.rs index e8ddab9..47798e8 100644 --- a/crates/clawhdf5/tests/integration_tests.rs +++ b/crates/clawhdf5/tests/integration_tests.rs @@ -986,3 +986,35 @@ fn u64_data_roundtrip() { values ); } + +// --------------------------------------------------------------------------- +// Opening a dataset by address +// --------------------------------------------------------------------------- + +#[test] +fn dataset_at_opens_the_same_dataset_as_its_path() { + let mut b = FileBuilder::new(); + let mut g = b.create_group("grp"); + g.create_dataset("vals").with_f64_data(&[1.0, 2.5, -3.0]); + b.add_group(g.finish()); + let file = File::from_bytes(b.finish().unwrap()).unwrap(); + + let addr = + clawhdf5_format::group_v2::resolve_path_any(file.as_bytes(), file.superblock(), "grp/vals") + .unwrap(); + let by_addr = file.dataset_at(addr).unwrap(); + assert_eq!(by_addr.read_f64().unwrap(), vec![1.0, 2.5, -3.0]); + assert_eq!( + by_addr.shape().unwrap(), + file.dataset("grp/vals").unwrap().shape().unwrap() + ); + + // The group's own header is not a dataset. + let group_addr = + clawhdf5_format::group_v2::resolve_path_any(file.as_bytes(), file.superblock(), "grp") + .unwrap(); + assert!(matches!( + file.dataset_at(group_addr), + Err(clawhdf5::Error::NotADataset(_)) + )); +} diff --git a/crates/clawhdf5/tests/v4_chunk_index_selection.rs b/crates/clawhdf5/tests/v4_chunk_index_selection.rs new file mode 100644 index 0000000..23786bd --- /dev/null +++ b/crates/clawhdf5/tests/v4_chunk_index_selection.rs @@ -0,0 +1,241 @@ +//! Partial hyperslab reads of every v4 (`libver='latest'`) chunk index type, +//! compared element for element with h5py. +//! +//! `Dataset::read_selection` takes two routes: `partial_read` materialises +//! the selection's bounding box when it covers at most half the dataset, and +//! `data_read::read_raw_data_selection` handles the rest. The second route +//! once passed the layout's full chunk dimensions (which carry the element +//! size as an extra, last dimension) to the implicit-index chunk generator, +//! which then indexed past the dataset's rank and panicked. So every case +//! below reads both a small window and one covering most of the dataset. +//! +//! Each case asserts which chunk index the file actually uses (parsed from +//! the layout message), so a change in how h5py lays the file out can't turn +//! this into a test of the wrong index. +//! +//! Skipped when python3 with h5py is unavailable, unless +//! `CLAWHDF5_REQUIRE_INTEROP=1`. + +use std::process::Command; + +use clawhdf5::File; +use clawhdf5_format::data_layout::DataLayout; +use clawhdf5_format::message_type::MessageType; +use clawhdf5_format::object_header::ObjectHeader; +use clawhdf5_format::selection::Selection; +use clawhdf5_format::superblock::Superblock; + +fn python() -> String { + std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string()) +} + +fn interop_required() -> bool { + std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1") +} + +fn python_available() -> bool { + Command::new(python()) + .args(["-c", "import h5py"]) + .output() + .map(|o| o.status.success()) + .unwrap_or(false) +} + +fn run_python(script: &str) -> String { + let output = Command::new(python()) + .args(["-c", script]) + .output() + .expect("failed to run python"); + assert!( + output.status.success(), + "Python script failed:\nSTDOUT: {}\nSTDERR: {}", + String::from_utf8_lossy(&output.stdout), + String::from_utf8_lossy(&output.stderr) + ); + String::from_utf8_lossy(&output.stdout).trim().to_string() +} + +/// The v4 chunk index type recorded in `name`'s layout message +/// (1 single chunk, 2 implicit, 3 fixed array, 4 extensible array, 5 B-tree v2). +fn chunk_index_type(path: &std::path::Path, name: &str) -> u8 { + let data = std::fs::read(path).unwrap(); + let sb = Superblock::parse(&data, 0).unwrap(); + let addr = clawhdf5_format::group_v2::resolve_path_any(&data, &sb, name).unwrap(); + let hdr = ObjectHeader::parse(&data, addr as usize, sb.offset_size, sb.length_size).unwrap(); + let msg = hdr + .messages + .iter() + .find(|m| m.msg_type == MessageType::DataLayout) + .expect("layout message"); + match DataLayout::parse(&msg.data, sb.offset_size, sb.length_size).unwrap() { + DataLayout::Chunked { + version: 4, + chunk_index_type: Some(t), + .. + } => t, + other => panic!("{name}: expected a v4 chunked layout, got {other:?}"), + } +} + +struct Case { + name: &'static str, + /// Python keyword arguments to `create_dataset` besides `data`. + kwargs: &'static str, + index_type: u8, +} + +const SHAPE: [u64; 2] = [37, 23]; + +const CASES: &[Case] = &[ + Case { + name: "implicit", + // Early allocation, no filters, fixed maximum: the implicit index. + kwargs: "chunks=(5, 4), dcpl=early()", + index_type: 2, + }, + Case { + name: "fixed_array", + kwargs: "chunks=(5, 4), compression='gzip'", + index_type: 3, + }, + Case { + name: "extensible_array", + kwargs: "chunks=(5, 4), maxshape=(None, 23), compression='gzip'", + index_type: 4, + }, + Case { + name: "btree2", + kwargs: "chunks=(5, 4), maxshape=(None, None), compression='gzip'", + index_type: 5, + }, + Case { + name: "single_chunk", + kwargs: "chunks=(37, 23), compression='gzip'", + index_type: 1, + }, + Case { + name: "single_chunk_unfiltered", + kwargs: "chunks=(37, 23)", + index_type: 1, + }, +]; + +/// `(start, stride, count, block)` per dimension; the first few stay below +/// half the dataset (bounding-box path), the rest exceed it (full path). +/// `(start, stride, count, block)` of a 2-D hyperslab. +type Hyperslab2 = ([u64; 2], [u64; 2], [u64; 2], [u64; 2]); + +fn selections() -> Vec { + vec![ + ([0, 0], [1, 1], [3, 23], [1, 1]), // ds[0:3] + ([7, 3], [1, 1], [9, 6], [1, 1]), // interior window across chunks + ([36, 22], [1, 1], [1, 1], [1, 1]), // last element (edge chunk) + ([2, 1], [3, 4], [4, 3], [1, 1]), // strided, small + ([0, 0], [1, 1], [30, 23], [1, 1]), // most rows + ([1, 0], [2, 1], [18, 23], [1, 1]), // every other row, spanning all + ([0, 2], [1, 1], [37, 20], [1, 1]), // columns 2..22 of every row + ([3, 1], [5, 3], [7, 7], [2, 2]), // strided blocks over everything + ] +} + +fn py_slice(start: u64, stride: u64, count: u64, block: u64) -> String { + // Each case is expressible as a numpy index when block == 1; with a block + // the selected indices are listed explicitly. + let idx: Vec = (0..count) + .flat_map(|c| (0..block).map(move |b| start + c * stride + b)) + .map(|i| i.to_string()) + .collect(); + format!("[{}]", idx.join(",")) +} + +#[test] +fn partial_hyperslabs_of_every_v4_chunk_index_match_h5py() { + if !python_available() { + assert!( + !interop_required(), + "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available" + ); + eprintln!("SKIP: python3 with h5py not available"); + return; + } + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("v4_index_selection.h5"); + let path_str = path.display().to_string(); + + // h5py writes the file, then reads every selection back and prints the + // values, one line per (case, selection). + let mut script = format!( + "import h5py, numpy as np\n\ + def early():\n\ + \x20 p = h5py.h5p.create(h5py.h5p.DATASET_CREATE)\n\ + \x20 p.set_alloc_time(h5py.h5d.ALLOC_TIME_EARLY)\n\ + \x20 return p\n\ + data = (np.arange({n}, dtype=' = expected + .next() + .expect("h5py printed too few lines") + .split_whitespace() + .map(|v| v.parse().unwrap()) + .collect(); + let raw = ds + .read_selection(&sel) + .unwrap_or_else(|e| panic!("{}: read_selection {sel:?} failed: {e}", case.name)); + let got: Vec = raw + .as_chunks::<4>() + .0 + .iter() + .map(|b| i32::from_le_bytes(*b)) + .collect(); + assert_eq!(got, want, "{}: selection {sel:?}", case.name); + assert_eq!( + ds.read_i32_selection(&sel).unwrap(), + want, + "{}: read_i32_selection {sel:?}", + case.name + ); + } + } + assert!(expected.next().is_none(), "h5py printed extra lines"); +} diff --git a/docs/QUICKSTART.md b/docs/QUICKSTART.md index 2ede177..34b5325 100644 --- a/docs/QUICKSTART.md +++ b/docs/QUICKSTART.md @@ -395,19 +395,22 @@ clawhdf5 --path agent.h5 snapshot backup_2026-03-19.h5 Read HDF5 files from Python without libhdf5: ```bash -pip install clawhdf5 # coming soon — build from source for now -cd crates/clawhdf5-py && maturin develop +# Not on PyPI yet: build from source into a virtualenv +pip install maturin numpy +cd crates/clawhdf5-py && maturin develop --release ``` ```python import clawhdf5 -# Read -f = clawhdf5.open("data.h5") -temps = f.read_f64("temperatures") -print(temps) # [22.5, 23.1, 21.8] +# Read (h5py-style) +with clawhdf5.File("data.h5", "r") as f: + temps = f["temperatures"][:] + print(temps) # [22.5 23.1 21.8] ``` +See `crates/clawhdf5-py/README.md` for the supported types and indexing. + --- ## Common Patterns diff --git a/docs/known-issues.md b/docs/known-issues.md index 3c12ff9..63a5fa2 100644 --- a/docs/known-issues.md +++ b/docs/known-issues.md @@ -7,6 +7,23 @@ deleting it. --- +## Selection reads that decode more than the selection + +**Status:** open (documented 2026-09-26). `Dataset::read_selection` (and so +the Python `ds[...]`) materialises only the selection's bounding box when +that box covers at most half the dataset (`partial_read`). It decodes the +whole dataset and extracts the selection instead when: +- the bounding box covers more than half the dataset — which a strided + selection across a chunked dataset (`ds[::100]`) always does, although + it may touch few chunks; +- the dataset is compact or virtual, or has no storage; +- it is chunked with a non-default fill value (the box path does not fill + unallocated chunks, so the fill-aware full read is used). +Values are correct in every case; this is cost only. Selections other than +`Selection::All` also bypass the file's chunk cache. The bounding-box +heuristic's other cost is measured under "Concurrent and contiguous read +performance" below. + ## Concurrent and contiguous read performance (measured 2026-09-26) **Status:** open for chunked full reads (one cause fixed 2026-09-26); the diff --git a/scripts/ci-test.sh b/scripts/ci-test.sh index e99cbf1..be25825 100755 --- a/scripts/ci-test.sh +++ b/scripts/ci-test.sh @@ -7,7 +7,9 @@ # # Environment: # CLAWHDF5_REQUIRE_INTEROP=1 Fail (instead of skip) when python3 with -# h5py/netCDF4/xarray is missing. CI sets this. +# h5py/netCDF4/xarray is missing, or without +# maturin/pytest for the Python package step. +# CI sets this. # Unset locally, the interop steps are skipped # if python3+h5py is not importable. # CLAWHDF5_FUZZ_SECONDS=N Run each cargo-fuzz target for N seconds @@ -55,7 +57,8 @@ run_step "cargo fmt --check" cargo fmt --check # 2. Clippy over every target (lib, bins, tests, benches, examples). Without # --all-targets, test and bench code is never linted. clawhdf5-py is -# excluded because it needs PyO3/Python headers. +# excluded here: PyO3 needs a Python interpreter to build, so it is +# linted in the Python package step (5b) instead. run_step "cargo clippy --all-targets" cargo clippy \ --workspace \ --exclude clawhdf5-py \ @@ -216,6 +219,36 @@ else STEPS+=("SKIP: h5py interop (format, ignored tests)") fi +# 5b. The Python package (crates/clawhdf5-py): lint it, build the wheel with +# maturin and run its pytest suite, which compares every read with h5py. +# The wheel is unpacked under target/ and put on PYTHONPATH, so the +# interpreter's environment is left as it was. Needs maturin and pytest +# in $PYTHON (CI installs both); skipped without them, and a failure +# instead when CLAWHDF5_REQUIRE_INTEROP=1. +python_package() { + local root="$SCRIPT_DIR/.." out + out="${CARGO_TARGET_DIR:-$root/target}/py-package" + rm -rf "$out" && mkdir -p "$out/wheel" "$out/site" || return 1 + cargo clippy -p clawhdf5-py --all-targets -- -D warnings || return 1 + "$PYTHON" -m maturin build \ + -m "$root/crates/clawhdf5-py/Cargo.toml" \ + -i "$PYTHON" \ + --out "$out/wheel" || return 1 + "$PYTHON" -m pip install --quiet --no-deps --target "$out/site" "$out"/wheel/*.whl || return 1 + PYTHONPATH="$out/site" "$PYTHON" -m pytest -q -p no:cacheprovider \ + "$root/crates/clawhdf5-py/tests" +} +if "$PYTHON" -m maturin --version >/dev/null 2>&1 && "$PYTHON" -c "import pytest" >/dev/null 2>&1; then + run_step "Python package (maturin build + pytest vs h5py)" python_package +elif [ "${CLAWHDF5_REQUIRE_INTEROP:-0}" = "1" ]; then + run_step "Python package (maturin build + pytest vs h5py)" \ + bash -c "echo \"maturin and pytest are required in $PYTHON (CLAWHDF5_REQUIRE_INTEROP=1)\"; exit 1" +else + echo "" + echo "==> [Python package] SKIPPED: needs maturin and pytest in $PYTHON" + STEPS+=("SKIP: Python package (maturin build + pytest)") +fi + # 6. Benches must keep compiling (they are not run). run_step "cargo bench --no-run" cargo bench \ --workspace \