diff --git a/.gitea/workflows/ci.yml b/.gitea/workflows/ci.yml index 06b42e3..fbca967 100644 --- a/.gitea/workflows/ci.yml +++ b/.gitea/workflows/ci.yml @@ -40,7 +40,7 @@ jobs: python3 -m venv /opt/interop # 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 + /opt/interop/bin/pip install --no-cache-dir h5py numpy netCDF4 xarray h5netcdf 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 cf90c2c..d0a22c0 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,6 @@ weights/ .venv __pycache__/ .pytest_cache/ + +# Scratch files the heavy tests generate (huge_chunks_interop) +crates/*/tests/scratch/ diff --git a/BENCHMARKS.md b/BENCHMARKS.md index 75024aa..c81101a 100644 --- a/BENCHMARKS.md +++ b/BENCHMARKS.md @@ -536,6 +536,116 @@ The rows and columns of the uncompressed layouts are within 20% (chunked column 0.45 -> 0.49 ms, contiguous column 2.55 -> 2.61 ms). This run does not explain the slower windows. +### HDF5 1.8 format: version-1 B-tree chunk indexes (2026-09-28, tank, loaded) + +**Superseded** by the idle re-run below; kept as the record the default +was first decided on. + +Measured 2026-09-28 on tank (AMD Ryzen 7 7800X3D), branch `feat/libver-v18` +at `3c61635`, to decide whether the writer's default should become the +HDF5 1.8 format (`libver_bounds(LibVer::V18, LibVer::V18)`: version-1 B-tree +chunk indexes) instead of the 1.10 format (Fixed Array indexes for these +datasets). **Not an idle machine:** two other agents were building; the +1-minute load average was 7.0 to 7.6 throughout (the rule is below 2), so +treat differences under about 20% as noise. Default and `--v18` runs +alternated, three of each per chunk size; medians of the three. + +> **Run:** `cargo run --release -p clawhdf5-bench --bin read_harness -- --chunk N [--v18]` +> with N = 256 (128 chunks per dataset) and N = 32 (8192 chunks per dataset) + +Write (the whole 3-dataset file) and file size: + +| chunks | 1.10 write ms | 1.8 write ms | 1.10 file bytes | 1.8 file bytes | size | +|---|---:|---:|---:|---:|---:| +| 256 x 256 | 467-500 | 467-483 | 134 916 080 | 134 933 848 | +0.013% | +| 32 x 32 | 490-495 | 489-506 | 138 879 336 | 139 465 112 | +0.42% | + +Reads (chunked datasets; the contiguous one does not change), ms. The +window labels are the harness's, which count 256 x 256 chunks: with 32 x 32 +chunks the 64 x 64 window covers 9 chunks and the 512 x 512 one 289. + +| chunks | layout | read | 1.10 | 1.8 | 1.8 / 1.10 | +|---|---|---|---:|---:|---:| +| 256 | chunked + deflate | full (first) | 6.70 | 7.30 | 1.09x | +| 256 | chunked + deflate | full (repeat) | 4.80 | 3.80 | 0.79x | +| 256 | chunked + deflate | 64 x 64 window (1 chunk) | 0.15 | 0.16 | 1.07x | +| 256 | chunked + deflate | 512 x 512 window (4-9 chunks) | 2.40 | 1.25 | 0.52x | +| 256 | chunked + deflate | one row | 0.95 | 0.95 | 1.00x | +| 256 | chunked + deflate | one column | 1.93 | 1.94 | 1.01x | +| 256 | chunked | full (first) | 11.30 | 11.80 | 1.04x | +| 256 | chunked | full (repeat) | 6.20 | 5.60 | 0.90x | +| 256 | chunked | 64 x 64 window (1 chunk) | 0.04 | 0.04 | 1.00x | +| 256 | chunked | 512 x 512 window (4-9 chunks) | 1.50 | 0.37 | 0.25x | +| 256 | chunked | one row | 0.05 | 0.05 | 1.00x | +| 256 | chunked | one column | 0.44 | 0.45 | 1.02x | +| 32 | chunked + deflate | full (first) | 11.50 | 12.20 | 1.06x | +| 32 | chunked + deflate | full (repeat) | 9.30 | 9.30 | 1.00x | +| 32 | chunked + deflate | 64 x 64 window (1 chunk) | 0.45 | 0.89 | 1.98x | +| 32 | chunked + deflate | 512 x 512 window (4-9 chunks) | 1.97 | 2.41 | 1.22x | +| 32 | chunked + deflate | one row | 0.69 | 1.13 | 1.64x | +| 32 | chunked + deflate | one column | 1.26 | 1.70 | 1.35x | +| 32 | chunked | full (first) | 11.80 | 13.00 | 1.10x | +| 32 | chunked | full (repeat) | 6.30 | 6.20 | 0.98x | +| 32 | chunked | 64 x 64 window (1 chunk) | 0.36 | 0.83 | 2.31x | +| 32 | chunked | 512 x 512 window (4-9 chunks) | 0.70 | 1.16 | 1.66x | +| 32 | chunked | one row | 0.38 | 0.84 | 2.21x | +| 32 | chunked | one column | 1.05 | 1.21 | 1.15x | + +With 128 chunks per dataset the two formats read and write alike (the 512 x +512 windows' 0.25x and 0.52x are not explained by the index and are likely +the load). With 8192 chunks, full reads stay within 10%, but a selection on +a freshly opened file costs about 0.4 to 0.5 ms more through the version-1 +B-tree (1.2x to 2.3x). Each timed selection opens the file anew, so the +likely cause (not profiled) is walking the B-tree's nodes (2.6 KB each, +about 150 per dataset here) against a Fixed Array's few blocks. Writing costs the +same; files grow by about 36 bytes per chunk. The default therefore stays +the 1.10 format; the 1.8 format is opt-in. + +### HDF5 1.8 format: version-1 B-tree chunk indexes, idle re-run (2026-09-28, tank) + +Measured 2026-09-28 on tank (AMD Ryzen 7 7800X3D), idle: 1-minute load +average 1.66 to 1.84 at the start of each run. Stacked branch +`feat/huge-chunks` at `b55768f` (contains `feat/libver-v18`). Default and +`--v18` runs alternated, five of each per chunk size; median (min-max), ms. + +> **Run:** `cargo run --release -p clawhdf5-bench --bin read_harness -- --chunk N [--v18]` +> with N = 256 (128 chunks per dataset) and N = 32 (8192 chunks per dataset) + +| chunks | layout | read | 1.10 | 1.8 | 1.8 / 1.10 | +|---|---|---|---:|---:|---:| +| 256 | chunked + deflate | full (first) | 5.9 (5.5-6.2) | 6.2 (6.0-6.7) | 1.05x | +| 256 | chunked + deflate | full (repeat) | 4.4 (4.1-4.9) | 4.3 (4.2-4.4) | 0.98x | +| 256 | chunked + deflate | 64 x 64 window (1 chunk) | 0.15 (0.15-0.17) | 0.16 (0.15-0.16) | 1.07x | +| 256 | chunked + deflate | 512 x 512 window (4-9 chunks) | 1.23 (1.23-1.36) | 1.22 (1.22-1.25) | 0.99x | +| 256 | chunked + deflate | one row | 0.95 (0.94-1.04) | 0.94 (0.93-0.98) | 0.99x | +| 256 | chunked + deflate | one column | 1.92 (1.91-2.12) | 1.90 (1.90-1.91) | 0.99x | +| 256 | chunked | full (first) | 11.1 (10.7-11.6) | 10.9 (10.6-11.9) | 0.98x | +| 256 | chunked | full (repeat) | 5.1 (4.9-5.2) | 5.2 (4.7-5.6) | 1.02x | +| 256 | chunked | 64 x 64 window (1 chunk) | 0.03 (0.03-0.03) | 0.04 (0.04-0.04) | 1.33x | +| 256 | chunked | 512 x 512 window (4-9 chunks) | 0.36 (0.36-0.38) | 0.37 (0.36-0.38) | 1.03x | +| 256 | chunked | one row | 0.05 (0.04-0.05) | 0.05 (0.05-0.05) | 1.00x | +| 256 | chunked | one column | 0.45 (0.43-0.45) | 0.44 (0.44-0.45) | 0.98x | +| 32 | chunked + deflate | full (first) | 11.7 (11.2-12.4) | 12.5 (11.8-13.2) | 1.07x | +| 32 | chunked + deflate | full (repeat) | 9.5 (9.3-9.8) | 9.4 (9.3-11.2) | 0.99x | +| 32 | chunked + deflate | 64 x 64 window (1 chunk) | 0.46 (0.45-0.49) | 0.89 (0.89-0.90) | 1.93x | +| 32 | chunked + deflate | 512 x 512 window (4-9 chunks) | 1.95 (1.94-2.04) | 2.40 (2.39-2.51) | 1.23x | +| 32 | chunked + deflate | one row | 0.69 (0.67-0.72) | 1.15 (1.13-1.18) | 1.67x | +| 32 | chunked + deflate | one column | 1.26 (1.24-1.28) | 1.72 (1.69-1.73) | 1.37x | +| 32 | chunked | full (first) | 12.5 (11.9-12.7) | 12.6 (12.5-13.4) | 1.01x | +| 32 | chunked | full (repeat) | 5.9 (5.8-6.0) | 6.2 (6.0-6.4) | 1.05x | +| 32 | chunked | 64 x 64 window (1 chunk) | 0.36 (0.35-0.36) | 0.84 (0.83-0.86) | 2.33x | +| 32 | chunked | 512 x 512 window (4-9 chunks) | 0.70 (0.69-0.74) | 1.17 (1.16-1.17) | 1.67x | +| 32 | chunked | one row | 0.38 (0.38-0.41) | 0.85 (0.84-0.86) | 2.24x | +| 32 | chunked | one column | 1.03 (1.02-1.21) | 1.22 (1.20-1.25) | 1.18x | + +Writes: 377 vs 374 ms (128 chunks), 399 vs 409 ms (8192 chunks); file +bytes as in the loaded run. The loaded run's 0.25x and 0.52x windows at 128 +chunks were the load: idle, every 128-chunk read is within 7% except the +0.03 ms single-chunk window (one timer tick). The 8192-chunk result stands: +a selection on a freshly opened file costs 0.2 to 0.5 ms more through the +version-1 B-tree (1.2x to 2.3x), full reads are within 7%. The default stays +the 1.10 format. + ## Local file speed after range reads ### `ObjectHeader::parse` back at 8f59b2e's speed (2026-09-27, tank) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6590714..88e255b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,189 @@ ## Unreleased +### NetCDF-4: variables' dimensions come from the file (2026-09-28) +- `clawhdf5-netcdf4` gave each variable the first unused dimension of + equal size (else an anonymous `dim_`), so a variable on an unlimited + dimension it had written fewer records of got `dim_`, and dimensions + of one size could be swapped. It now resolves them as netCDF-C does + (`libhdf5/hdf5open.c`): the ids in the variable's `_Netcdf4Coordinates` + (each scale's `_Netcdf4Dimid`), else the scales its `DIMENSION_LIST` + references (object references, also the revised `H5T_STD_REF`; of + several scales attached to one axis, the last, as netCDF-C's + `dimscale_visitor` keeps), looked up in its group and then each parent + group; a coordinate variable is on + its own scale. Only an axis with neither (a file not written by a netCDF + library) is still matched by size. A variable may use one dimension + twice (`(p, p)`). +- `variables()` and `variable_names()` leave out the dimension scales that + are only dimensions (`NAME` "This is a netCDF dimension but not a netCDF + variable."), as netCDF-C does, and `variable()` refuses them + (`VariableNotFound`). A variable stored as `_nc4_non_coord_` + (netCDF-C's name for a variable sharing a dimension's name without being + its coordinate variable) is listed and found as ``. + `NetCDF4File::variable_names` is new; `NetCDF4Group::variable_names` used + to list every dataset. +- A variable along an unlimited dimension has the dimension's length, as in + netCDF: `Variable::shape` is that length and the reads return that many + values, the records it has not written as the fill value (`_FillValue`, + else netCDF's `NC_FILL_*` for the type; `""` for strings; NaN from + `read_f64`). It was the HDF5 extent. `Variable::stored_shape` (new) is the + HDF5 extent. Where the unlimited dimension is not a variable's first, + unwritten values are placed per row, as netCDF-C's element and row reads + return them; a whole-variable read through netCDF-C 4.9.3 + (netCDF4-python 1.7.4) instead returns the written values first, then the + fill. +- A variable's attributes are read when it is opened (they were read on + first use). +- Tests, compared with netCDF4-python 1.7.4 (netCDF-C 4.9.3) variable by + variable (dimensions, shape, every value): the known-issues reproducer; + two dimensions of one size in either order, one dimension used twice, a + scalar, a non-coordinate variable named like a dimension, a subgroup and + a sub-subgroup on their ancestors' dimensions; unwritten records of + `i1`/`i4`/`u8`/`f4`/`f8`/string variables, with and without + `_FillValue`; a file of h5py dimension scales (no netCDF attributes); + h5netcdf 1.8.1 and xarray 2026.7.0 files (engines netcdf4 and h5netcdf). + The h5netcdf cases skip when h5netcdf is not installed; CI now installs + it. A one-off comparison over the 78 conformance-corpus files + netCDF4-python opens (tank, 2026-09-28; a throwaway test, not committed) + found 52 files with the same variables, dimensions and shapes, and the + same values in every numeric variable of up to 5000 elements (enum + variables were not compared). The other 26 have no dimension scales + (netCDF-C names their axes `phony_dim_`; this crate still matches by + size or names them `dim_`) or hold datasets of types netCDF-C + skips (opaque, references). Affected v2.1.0 to v2.7.0. + `docs/known-issues.md`. + +### Writing files HDF5 1.8 can read (2026-09-28) +- New `LibVer` (`V18`, `V110`, `V112`, `V114`, `V200`, `Latest`; re-exported + as `clawhdf5::LibVer`) and `FileBuilder::libver_bounds(low, high)` + (`FileWriter::libver_bounds` in the format crate), as libhdf5's + `H5Pset_libver_bounds` and h5py's `libver=(low, high)`. The default, + `(V110, Latest)`, writes exactly what was written before (the HDF5 1.10 + format, which HDF5 1.8.23 refuses to open). +- A low bound of `V18` writes what libhdf5 2.x writes for h5py's + `libver=('v108', 'latest')`: a version-2 superblock (version 3 only for a + paged file), version-3 layout messages for contiguous, compact and chunked + datasets, and a version-1 B-tree chunk index for every chunked dataset, + resizable and multi-unlimited ones included, instead of the single chunk, + Fixed Array, Extensible Array and version-2 B-tree indexes. The B-tree + writer (`btree_v1_write`) replays `H5B_insert` with the chunk callbacks of + `H5Dbtree.c` for chunks arriving in row-major order: libhdf5's split + ratios, right keys moved exactly when `H5D__btree_cmp3` moves them, the + root kept at its address. Its trees are libhdf5's node for node (levels, + child counts, keys) for 1-D, 2-D and 3-D datasets with two- and + three-level trees, deflated or not, against libhdf5 2.0 writing without a + chunk cache (`chunk_btrees_match_libhdf5`). An empty chunked dataset has + no tree (undefined address), as in libhdf5. +- A high bound refuses, with the new `FormatError::LibverBound` and before + anything is written, what needs a newer format: virtual datasets and the + paged file-space strategy (1.10), the 1.12 reference types (datatype + version 4), native complex numbers (datatype version 5, HDF5 2.0), and a + low bound above the high one. `(V18, V18)` therefore writes a file HDF5 + 1.8 reads, or fails; `(V18, Latest)` writes such objects in their newer + format, as libhdf5 does. `Datatype::max_encoded_version` reports the + version a type needs. +- Checked against a real HDF5 1.8: `scripts/build-hdf5-1.8.sh` builds + 1.8.23 (the last 1.8 release) with its tools. `clawhdf5-tools`' + `tests/libver_v18.rs` writes every writer feature under `(V18, V18)` + (contiguous, empty, scalar, compact, f16/f32/f64/i64, chunked with + deflate + shuffle + Fletcher-32 and edge chunks, resizable with one and + two unlimited dimensions, finite maxshape, 100 000 one-element chunks for + a three-level tree, a 100 x 100 grid of chunks, fill values, fixed-length + strings, compound, enum, array, compact/dense/creation-ordered groups, + soft, hard and external links, dense attributes); HDF5 1.8.23's h5dump + dumps the whole file exactly as h5dump 1.14 does and returns our bytes + for every numeric dataset (`-b LE`), and h5py, clawhdf5 and + `h5rs check --data` read it. Then `FileEditor` appends to the resizable + datasets 300 times (splitting B-tree nodes), grows the three-level one, + rewrites a deflated row, grows a 2-D dataset and sets compact and dense + attributes; h5py appends too; and every reader checks again. The 1.8 + checks are skipped where no HDF5 1.8 is found (`CLAWHDF5_H5DUMP18` or + `~/.cache/hdf5-1.8.23`), as in CI. +- The default stays the 1.10 format: with 8192 chunks per dataset, small + selections on a freshly opened file read 1.2x to 2.3x slower through a + version-1 B-tree (read harness, tank, 2026-09-28, under load; full reads + and writes within 10%, files about 36 bytes per chunk larger). + `read_harness` gains `--v18` and `--chunk N`. `BENCHMARKS.md`, "HDF5 1.8 + format". +- Not yet: the Python bindings' `'w'` mode has no `libver` argument, and + the pre-1.8 format (version-0 superblock, symbol-table groups) cannot be + written. `docs/known-issues.md`. + +### Chunks of 4 GiB or more (HDF5 2.0's layout message version 5) (2026-09-28) +- **What libhdf5 does** (read in libhdf5 2.2.0's `H5Dchunk.c`, `H5Dfarray.c`, + `H5Dearray.c`, `H5Dbtree2.c`, `H5Dbtree.c`, `H5Olayout.c`): a chunk of + more than 0xFFFFFFFF bytes requires layout message version 5, so a high + bound of `H5F_LIBVER_V200` ("chunk size > 4GB requires H5F_LIBVER_V200"), + whatever the low bound, and version 5 always uses the newer indexes. + Under version 5 a filtered Fixed Array, Extensible Array or v2 B-tree + element stores the chunk's size in "size of lengths" bytes (8); under + version 4 in one byte more than the chunk needs. A filtered Single Chunk + stores it in "size of lengths" bytes under both. Unfiltered elements + store no size, and the Implicit index none at all. A v1 B-tree key keeps + a 32-bit size: libhdf5 never writes a larger chunk there and refuses to + open one ("chunk size must be < 4GB with v1 b-tree index"). +- **Reading** such chunks works in every index: `ChunkInfo::chunk_size` and + `ChunkMapping::file_size` are now `u64` (**public API change**). The + Single Chunk, Implicit, Fixed and Extensible Array readers truncated an + unfiltered chunk's size (the read failed), and the v2 B-tree reader + refused any chunk over 4 GiB. Checked against a fixture libhdf5 2.0.0 + wrote (`tests/fixtures/huge_chunks_filtered.h5`, 188 KB: a 4 GiB + 8 + byte chunk in each filtered index, deflated twice) and a 44 GiB sparse + file h5py writes at test time (every unfiltered index). +- **Selections read what they touch:** a selection of a chunked dataset with + a non-default fill value is read over a box of fill values from the + chunks it touches (new `partial_read::read_selection_filled_in`); it + used to decode the whole dataset. An unfiltered chunk of a file that is + not in memory (`File::open_storage`) is read row by row instead of whole. + A filtered chunk is still decoded whole. +- **Writing:** a chunk of more than `u32::MAX` bytes gets layout message + version 5 and libhdf5's index element widths; the Fixed and Extensible + Array structures match libhdf5 2.0.0's byte for byte. New public helpers + `chunked_write::{layout_version_for, chunk_size_len, MAX_V4_CHUNK_BYTES}`. + Refused before anything is written: chunk dimensions of 2^32 or more + (they were cut to 32 bits) and, for such chunks, LZF, bitshuffle, bzip2, + Blosc and pcodec. Chunks are extracted row by row and one at a time + (compressed in parallel only up to 64 MiB), and a filter pipeline no + longer copies its input first. +- **Deflate** compressed only the first 4 GiB - 1 bytes of a larger chunk + (zlib takes at most that much per call, and `Finish` ended the stream + there; since the one-pass deflate of 2026-09-23, in `clawhdf5-format` and + `clawhdf5-filters`), and it reserved the compressor's worst-case bound, + which flate2's Rust backends zero on every call: 4 GiB kept per 4 GiB + chunk. Inputs over 64 MiB now start at 1/16 of the bound and grow, and a + result more than 64 MiB too large is shrunk. On the read side, an + intermediate deflate stage (one followed by a filter other than shuffle + or Fletcher-32) no longer reserves the chunk's whole bound: reading the + double-deflated fixture peaked at 8.5 GiB, now 4.0 GiB. +- **LZ4:** a chunk decoding to more than 256 MiB was refused ("lz4: + declared size exceeds limit") though the chunk size bounded it; the + ceiling now applies only to a decode of unknown size. A chunk of 4 GiB + or more is read as the registered HDF5 framing, whose 64-bit size no + longer starts with four zero bytes. +- **`FileEditor`** refuses, before writing anything, to write values into + or to prune, fill or allocate chunks of 4 GiB or more; growing the extent + (late allocation) and setting attributes work. It sizes a version-5 + filtered index element from the file's size of lengths (it assumed 8). +- **32-bit targets:** reading or writing such a chunk is + `FormatError::Overflow`; `scripts/check-32bit-casts.sh` passes, and the + wasm package test reads the fixture and gets that error in every index. +- Tests: `crates/clawhdf5/tests/huge_chunks_interop.rs` (index listing, + byte-for-byte array indexes and the editor always; with + `CLAWHDF5_HUGE_CHUNKS=1`, reads of every index filtered and unfiltered, + and a write of every index read back by clawhdf5, h5py 3.16 and — layout + and storage size — h5dump 2.2.0 via `CLAWHDF5_H5DUMP2`), + `partial_read_equivalence::fill_value_selections_match_full_reads`, unit + tests of the encodings (`chunked_write::tests`) and of LZ4. The opt-in + run (`cargo test --release -p clawhdf5 --test huge_chunks_interop -- + --test-threads=1`) peaked at 8.06 GiB resident, h5py included (tank, + 2026-09-28, commit `9143737`). libhdf5 2.2.0's h5dump cannot print these + chunks' values (its deflate filter fails on any chunk over 4 GiB, + libhdf5's own included); h5py 3.16 reads them. `docs/known-issues.md`: three fixed + entries and the open "Chunks of 4 GiB or more: limits" (chunk dimensions + of 2^32 or more, the refused filters, memory, the untested unfiltered + write). + ### A dropped `FileEditor` releases its lock at once (2026-09-28) - `FileEditor`'s `flock` could outlive the editor for a moment when another thread forked to spawn a process: the child shared the locked diff --git a/README.md b/README.md index 87758ad..d7f61da 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,8 @@ It reads superblocks v0–3, every group and chunk-index structure libhdf5 writes, the standard filters and the common plugin filters, variable-length data and virtual datasets, and follows files a SWMR writer is appending to. It reads files from libhdf5, h5py and netCDF-4 and writes -files they read. The same library opens files over HTTP and in object +files they read, in the HDF5 1.10 format or, on request, in one HDF5 1.8 +reads. The same library opens files over HTTP and in object stores by range requests, runs in the browser as WebAssembly, and has Python bindings with an h5py-shaped API. @@ -112,12 +113,12 @@ Limits and open issues, with dates, are in | Area | Supported | Read only | Not supported | |---|---|---|---| -| **File format** | Superblock v0–v3, user blocks, v1/v2 object headers | Metadata cache images | Writing files HDF5 1.8 can read | +| **File format** | Superblock v0–v3, user blocks, v1/v2 object headers; writing the HDF5 1.10 format (default) or, with `libver_bounds(LibVer::V18, LibVer::V18)`, files HDF5 1.8 reads (checked with HDF5 1.8.23) | Metadata cache images | Writing the pre-1.8 format (version-0 superblock, symbol-table groups) | | **Groups and links** | Symbol-table, compact and dense groups (tested to 100 000 links), creation order, soft and hard links; writing external links | | Following external links (explicit error); user-defined links are skipped | | **Datatypes** | Integers and IEEE floats of every width and byte order (incl. `f16`), enums, compounds (every version, incl. HDF5 2.0's v5), arrays, fixed-length strings, opaque, complex: h5py's `{r, i}` compound (`with_complex_f64_data`) and HDF5 2.0's native class 11 (`with_native_complex_f64_data`, opt-in: only libhdf5 2.0+ reads it; reads surface it as `{r, i}`) | Variable-length strings and sequences, object references; HDF5 2.x's small floats (bfloat16, FP8 E4M3/E5M2, FP6 E2M3/E3M2, FP4 E2M1: every bit pattern decoded as libhdf5 2.2.0 decodes it) and other non-IEEE floats up to 64 bits | Writing variable-length data; writing non-IEEE floats; decoding region and attribute references; x87 long double and binary128 | -| **Layouts and chunk indexes** | Compact, contiguous and chunked; chunk indexes single chunk, Fixed Array, Extensible Array and v2 B-tree (the writer picks one as libhdf5 does); fill values; resizable datasets; virtual datasets (read limits in known-issues) | Chunk indexes v1 B-tree and implicit (the editor also changes them) | External raw data files (explicit error) | +| **Layouts and chunk indexes** | Compact, contiguous and chunked; chunk indexes single chunk, Fixed Array, Extensible Array and v2 B-tree (the writer picks one as libhdf5 does), and v1 B-tree (every chunked dataset under the 1.8 bound, split as libhdf5 splits it); chunks of 4 GiB or more (HDF5 2.0's layout message version 5); fill values; resizable datasets; virtual datasets (read limits in known-issues) | The implicit chunk index (the editor also changes it) | External raw data files (explicit error); chunk dimensions of 2^32 or more | | **Filters** | deflate (pure-Rust zlib-rs), shuffle, Fletcher-32, LZ4 (opt-in), Zstd (C, opt-in); plugins LZF, bitshuffle, bzip2, Blosc 1 | N-Bit, scale-offset, SZIP (C, opt-in); plugins Blosc2 and ZFP | Other filter IDs, unless you register a codec (`filter_registry::register_filter`) | -| **Editing in place** | `FileEditor`: overwrite values, grow and shrink chunked datasets (every index), set attributes (compact and dense), in files from h5py or clawhdf5 | | Creating or deleting objects in an existing file; deleting attributes; new chunks in implicit indexes; VL data; filters this build cannot encode (refused before any write) | +| **Editing in place** | `FileEditor`: overwrite values, grow and shrink chunked datasets (every index), set attributes (compact and dense), in files from h5py or clawhdf5 | | Creating or deleting objects in an existing file; deleting attributes; new chunks in implicit indexes; VL data; rewriting chunks of 4 GiB or more; filters this build cannot encode (refused before any write) | | **Access** | Local files (mmap or buffered), bytes in memory, any `Storage` backend, HTTP(S) and S3/GCS/Azure via `clawhdf5-remote`, SWMR reading (`File::open_swmr`, `Dataset::refresh`) | Remote files and the browser are read-only | SWMR writing; remote SWMR; MPI collective I/O (`clawhdf5-io`'s `mpi-io` reads on one rank and broadcasts) | | **Bindings** | Python (read, `'w'` for numeric and complex arrays, `'r+'` editing, URLs), NetCDF-4 (CF scale/offset/fill) | WebAssembly (`open(bytes)`, `openUrl`); no Zstd/SZIP/pcodec, no compound, reference, opaque, bitfield, time or VL-sequence datasets | Node.js (the package does not work; see known-issues) | @@ -411,7 +412,8 @@ scripts/ci-test.sh # what CI runs: fmt, clippy matrix, tests, conformance/run.sh # the conformance report (needs h5py, hdf5plugin, h5dump) ``` -The interop suites need a Python with h5py (and netCDF4, xarray); on a +The interop suites need a Python with h5py (and netCDF4, xarray; the +NetCDF-4 tests of h5netcdf-written files skip without h5netcdf); on a PEP 668 system that has to be a virtualenv, which `ci-test.sh` finds as `.venv` or through `CLAWHDF5_PYTHON`. Without one they skip; set `CLAWHDF5_REQUIRE_INTEROP=1` to make that a failure, as CI does: diff --git a/crates/clawhdf5-bench/src/bin/read_harness.rs b/crates/clawhdf5-bench/src/bin/read_harness.rs index b09879a..0a517d6 100644 --- a/crates/clawhdf5-bench/src/bin/read_harness.rs +++ b/crates/clawhdf5-bench/src/bin/read_harness.rs @@ -7,15 +7,19 @@ //! ```text //! cargo run --release -p clawhdf5-bench --bin read_harness //! cargo run --release -p clawhdf5-bench --bin read_harness -- --large # 512 MB +//! cargo run --release -p clawhdf5-bench --bin read_harness -- --v18 # HDF5 1.8 format +//! cargo run --release -p clawhdf5-bench --bin read_harness -- --chunk 32 # 32 x 32 chunks //! ``` +//! +//! `--v18` writes the file with `libver_bounds(V18, V18)` (version-1 B-tree +//! chunk indexes) instead of the default 1.10 format (Fixed Array indexes +//! here), to compare the two. use std::time::{Duration, Instant}; -use clawhdf5::{File, FileBuilder}; +use clawhdf5::{File, FileBuilder, LibVer}; use clawhdf5_format::selection::Selection; -const CHUNK: u64 = 256; - struct Layout { name: &'static str, chunked: bool, @@ -46,16 +50,19 @@ fn value(row: u64, col: u64) -> f64 { (row * 100_003 + col) as f64 * 0.5 } -fn write_file(path: &std::path::Path, rows: u64, cols: u64) { +fn write_file(path: &std::path::Path, rows: u64, cols: u64, chunk: u64, v18: bool) { let data: Vec = (0..rows) .flat_map(|r| (0..cols).map(move |c| value(r, c))) .collect(); let mut builder = FileBuilder::new(); + if v18 { + builder.libver_bounds(LibVer::V18, LibVer::V18); + } for (i, layout) in LAYOUTS.iter().enumerate() { let ds = builder.create_dataset(&format!("d{i}")); ds.with_f64_data(&data).with_shape(&[rows, cols]); if layout.chunked { - ds.with_chunks(&[CHUNK, CHUNK]); + ds.with_chunks(&[chunk, chunk]); } if layout.deflate { ds.with_deflate(4); @@ -91,7 +98,14 @@ fn slab(start: [u64; 2], count: [u64; 2]) -> Selection { } fn main() { - let large = std::env::args().any(|a| a == "--large"); + let args: Vec = std::env::args().collect(); + let large = args.iter().any(|a| a == "--large"); + let v18 = args.iter().any(|a| a == "--v18"); + let chunk: u64 = args + .iter() + .position(|a| a == "--chunk") + .and_then(|i| args.get(i + 1)) + .map_or(256, |c| c.parse().expect("--chunk N")); let (rows, cols) = if large { (8192, 8192) } else { (4096, 2048) }; let total_mb = (rows * cols * 8) as f64 / (1 << 20) as f64; if cfg!(debug_assertions) { @@ -100,12 +114,21 @@ fn main() { let dir = tempfile::TempDir::new().unwrap(); let path = dir.path().join("read_harness.h5"); - write_file(&path, rows, cols); - let file_mb = std::fs::metadata(&path).unwrap().len() as f64 / (1 << 20) as f64; + let t = Instant::now(); + write_file(&path, rows, cols, chunk, v18); + let write_ms = t.elapsed().as_secs_f64() * 1e3; + let file_bytes = std::fs::metadata(&path).unwrap().len(); + let file_mb = file_bytes as f64 / (1 << 20) as f64; println!("## Read harness"); println!( - "\n{rows} x {cols} f64 ({total_mb:.0} MB per dataset), chunks {CHUNK} x {CHUNK}, file {file_mb:.0} MB\n" + "\n{rows} x {cols} f64 ({total_mb:.0} MB per dataset), chunks {chunk} x {chunk}, \ + format {}, file {file_mb:.0} MB ({file_bytes} bytes), written in {write_ms:.0} ms\n", + if v18 { + "1.8 (v1 B-tree)" + } else { + "1.10 (default)" + } ); // (label, selection, elements selected) diff --git a/crates/clawhdf5-filters/src/fast_deflate.rs b/crates/clawhdf5-filters/src/fast_deflate.rs index 6ce268f..2b4e5f0 100644 --- a/crates/clawhdf5-filters/src/fast_deflate.rs +++ b/crates/clawhdf5-filters/src/fast_deflate.rs @@ -327,24 +327,53 @@ fn inflate_bounded(data: &[u8], size_hint: usize, limit: usize) -> Result Result, String> { use flate2::{Compress, Compression, FlushCompress, Status}; // zlib's compressBound, plus the zlib header and trailer. let bound = data.len() + (data.len() >> 12) + (data.len() >> 14) + (data.len() >> 25) + 13 + 6; + // flate2's Rust backends (zlib-rs, miniz_oxide) zero the whole spare + // capacity on each call, so a large input's worst-case bound would be + // memory held for nothing (4 GiB for a 4 GiB chunk that deflates to a + // few MiB): past 64 MiB the output starts at 1/16 of the bound and + // doubles as needed. + let first = if bound <= DEFLATE_EXACT_BOUND { + bound + } else { + bound / 16 + }; let mut out = Vec::new(); - out.try_reserve_exact(bound) + out.try_reserve_exact(first) .map_err(|e| format!("deflate: cannot allocate output: {e}"))?; let mut deflater = Compress::new(Compression::new(level), true); loop { let (in_before, out_before) = (deflater.total_in(), deflater.total_out()); + let rest = &data[in_before as usize..]; + // zlib takes at most u32::MAX input bytes per call, and `Finish` + // ends the stream after the bytes it took: input of 4 GiB or more + // was cut at 4 GiB - 1. Finish only once the rest fits one call. + let flush = if rest.len() > u32::MAX as usize { + FlushCompress::None + } else { + FlushCompress::Finish + }; let status = deflater - .compress_vec(&data[in_before as usize..], &mut out, FlushCompress::Finish) + .compress_vec(rest, &mut out, flush) .map_err(|e| format!("deflate: {e}"))?; match status { - Status::StreamEnd => return Ok(out), + Status::StreamEnd => { + if out.capacity() - out.len() > DEFLATE_EXACT_BOUND { + out.shrink_to_fit(); + } + return Ok(out); + } + // Out of room (the bound makes it unreachable below + // `DEFLATE_EXACT_BOUND`): grow rather than fail. Status::Ok | Status::BufError if out.len() == out.capacity() => out .try_reserve(out.capacity().max(4096)) .map_err(|e| format!("deflate: cannot allocate output: {e}"))?, diff --git a/crates/clawhdf5-format/README.md b/crates/clawhdf5-format/README.md index 2223b55..9b3923a 100644 --- a/crates/clawhdf5-format/README.md +++ b/crates/clawhdf5-format/README.md @@ -36,7 +36,9 @@ clawhdf5-format = { git = "https://git.redclaw.dev/quantumclaw/clawhdf5" } `type_builders` (datasets, groups, attributes, compound and enum types, links, virtual datasets, creation-order tracking); chunk indexes and dense-storage B-trees of any size (`chunked_write`, `btree_v2_write`, - `ea_writer`). Output is read by h5py and h5dump. + `ea_writer`, and the version-1 chunk B-tree of `btree_v1_write`). Output + is read by h5py and h5dump; `FileWriter::libver_bounds` (`libver`) picks + the format: HDF5 1.10 by default, or one HDF5 1.8 reads. - **Filters:** `filter_pipeline` and `filter_registry` (look up by ID; other IDs can be registered at run time with `register_filter`). Built in: deflate, shuffle, Fletcher-32, N-Bit, scale-offset; behind features LZ4, diff --git a/crates/clawhdf5-format/src/btree_v1_write.rs b/crates/clawhdf5-format/src/btree_v1_write.rs new file mode 100644 index 0000000..c0524d8 --- /dev/null +++ b/crates/clawhdf5-format/src/btree_v1_write.rs @@ -0,0 +1,470 @@ +//! Writing a version-1 B-tree chunk index (node type 1): the chunk index of +//! layout message versions 1-3, and the only one HDF5 1.8 reads. +//! +//! The tree is built the way libhdf5 builds it when the chunks reach it one +//! after another in row-major order (a whole-dataset `H5Dwrite` of a 1-D +//! dataset, or of any dataset without a chunk cache; with one, libhdf5 +//! inserts the small chunks of a multi-dimensional dataset in the order its +//! cache evicts them, which fills the nodes differently): each +//! chunk goes through the same steps as `H5B_insert` (`H5B.c`) with the +//! chunk callbacks of `H5Dbtree.c`, so nodes split where libhdf5's split, +//! with its default split ratios (a full right-most node keeps 90% of its +//! children, a left-most one 10%, any other half), and keys hold what +//! libhdf5's hold: +//! +//! - a chunk's key is its size in the file, its filter mask and its offsets +//! (the element-size coordinate 0); +//! - a node's final key is the zero-size key one chunk past the chunk that +//! last moved it (every scaled coordinate plus one, `H5D__btree_new_node`), +//! which libhdf5 moves only when a new chunk is not below it +//! (`H5D__btree_cmp3`) — so after an even number of appends in one +//! dimension it lies on the last chunk itself; +//! - a full root is copied to a new node and becomes the parent of the copy +//! and its new sibling, so the root's address (the layout message's) never +//! changes. +//! +//! Nodes are laid out in the order libhdf5 allocates them (the root first, +//! then each new node as a split creates it), all of the full node size, the +//! unused slots zero. + +#[cfg(not(feature = "std"))] +use alloc::{format, vec, vec::Vec}; + +use core::cmp::Ordering; + +use crate::error::FormatError; + +/// libhdf5's default chunk B-tree K (`HDF5_BTREE_CHUNK_IK_DEF`): nodes hold +/// up to 2K = 64 children. Superblocks of version 2 cannot record another +/// value without a superblock extension, which this writer does not emit. +pub(crate) const CHUNK_BTREE_K: u16 = 32; + +/// libhdf5's default split ratios (`H5D_XFER_BTREE_SPLIT_RATIO_DEF`) for a +/// left-most, middle and right-most node. +const SPLIT_RATIOS: [f64; 3] = [0.1, 0.5, 0.9]; + +/// A chunk to index: scaled coordinates (offset / chunk dimension) in each +/// dataset dimension, stored size, filter mask and address. +pub(crate) struct ChunkEntry { + pub(crate) scaled: Vec, + pub(crate) nbytes: u64, + pub(crate) filter_mask: u32, + pub(crate) address: u64, +} + +#[derive(Debug, Clone, PartialEq, Eq)] +struct Key { + nbytes: u32, + mask: u32, + /// Scaled coordinates, the element-size one (0 or 1) last. + scaled: Vec, +} + +impl Key { + /// `H5D__btree_new_node`'s right key: one chunk past `self` in every + /// dimension, with no storage. + fn right_of(&self) -> Key { + Key { + nbytes: 0, + mask: 0, + scaled: self.scaled.iter().map(|s| s + 1).collect(), + } + } + + fn cmp_scaled(&self, other: &Key) -> Ordering { + self.scaled.cmp(&other.scaled) + } +} + +#[derive(Debug, Clone)] +struct Node { + level: u8, + left: Option, + right: Option, + /// `children.len() + 1` keys once the node holds a child. + keys: Vec, + /// Chunk addresses in a leaf, node indexes above. + children: Vec, +} + +/// What an insertion below a node did (`H5B__insert_helper`'s outputs). +#[derive(Default)] +struct Ret { + /// The node's new left key (`lt_key_changed`). + lt: Option, + /// The node's new right key (`rt_key_changed`). + rt: Option, + /// The node split: the key shared by the halves and the new right node. + split: Option<(Key, usize)>, +} + +struct Tree { + nodes: Vec, + two_k: usize, +} + +fn bad(why: &str) -> FormatError { + FormatError::SerializationError(format!("version-1 B-tree chunk index: {why}")) +} + +impl Tree { + fn new(k: u16) -> Self { + Self { + nodes: vec![Node { + level: 0, + left: None, + right: None, + keys: Vec::new(), + children: Vec::new(), + }], + two_k: 2 * usize::from(k), + } + } + + /// `H5B_insert` of `key` (a chunk after every chunk already inserted). + fn insert(&mut self, key: &Key, addr: u64) -> Result<(), FormatError> { + let r = self.insert_helper(0, key, addr, 64)?; + let Some((md, split)) = r.split else { + return Ok(()); + }; + // The root split: copy it to a new node and make the root the + // parent of the copy and its new right sibling. + let lt = r.lt.unwrap_or_else(|| self.nodes[0].keys[0].clone()); + let rt = match r.rt { + Some(rt) => rt, + None => self.nodes[split] + .keys + .last() + .cloned() + .ok_or_else(|| bad("empty node"))?, + }; + let moved = self.nodes[0].clone(); + let level = moved.level; + let moved_id = self.nodes.len(); + self.nodes.push(moved); + self.nodes[split].left = Some(moved_id); + self.nodes[0] = Node { + level: level + 1, + left: None, + right: None, + keys: vec![lt, md, rt], + children: vec![moved_id as u64, split as u64], + }; + Ok(()) + } + + fn insert_helper( + &mut self, + id: usize, + key: &Key, + addr: u64, + depth: u8, + ) -> Result { + if depth == 0 { + return Err(bad("tree too deep")); + } + let n = self.nodes[id].children.len(); + let level = self.nodes[id].level; + let mut ret = Ret::default(); + if n == 0 { + // The first chunk (H5B_INS_FIRST): its key and the right key. + let node = &mut self.nodes[id]; + node.keys = vec![key.clone(), key.right_of()]; + node.children = vec![addr]; + return Ok(ret); + } + // Binary search with H5D__btree_cmp3: 1 when the chunk is not below + // the right key, -1 when below the left key, else 0. + let (mut lo, mut hi, mut idx) = (0usize, n, 0usize); + let mut cmp = Ordering::Less; + while lo < hi && cmp != Ordering::Equal { + idx = (lo + hi) / 2; + let node = &self.nodes[id]; + cmp = if key.cmp_scaled(&node.keys[idx + 1]) != Ordering::Less { + Ordering::Greater + } else if key.cmp_scaled(&node.keys[idx]) == Ordering::Less { + Ordering::Less + } else { + Ordering::Equal + }; + if cmp == Ordering::Less { + hi = idx; + } else { + lo = idx + 1; + } + } + let (mut lt_changed, mut rt_changed) = (false, false); + // The child to add after child `idx`, with its left key. + let mut new_child: Option<(Key, u64)> = None; + match cmp { + Ordering::Less => return Err(bad("chunks out of order")), + Ordering::Greater if idx + 1 < n => { + return Err(bad("cannot place chunk")); + } + Ordering::Greater if level == 0 => { + // Past every chunk of the right-most leaf: a new maximum + // (H5B_INS_RIGHT through `new_node`), which moves the right + // key one chunk past it. + idx = n - 1; + self.nodes[id].keys[idx + 1] = key.right_of(); + rt_changed = true; + new_child = Some((key.clone(), addr)); + } + Ordering::Equal if level == 0 => { + // Inside the last chunk's range: H5D__btree_insert adds it + // to the right of that chunk; the right key stays. + if key.scaled == self.nodes[id].keys[idx].scaled { + return Err(bad("duplicate chunk")); + } + new_child = Some((key.clone(), addr)); + } + _ => { + if cmp == Ordering::Greater { + idx = n - 1; + } + let child = usize::try_from(self.nodes[id].children[idx]) + .map_err(|_| bad("bad node index"))?; + let r = self.insert_helper(child, key, addr, depth - 1)?; + if let Some(lt) = r.lt { + self.nodes[id].keys[idx] = lt; + lt_changed = true; + } + if let Some(rt) = r.rt { + self.nodes[id].keys[idx + 1] = rt; + rt_changed = true; + } + if let Some((md, split)) = r.split { + new_child = Some((md, split as u64)); + } + } + } + // Pass the node's changed end keys up, as H5B__insert_helper does. + if lt_changed && idx == 0 { + ret.lt = Some(self.nodes[id].keys[0].clone()); + } + if rt_changed && idx + 1 >= n { + ret.rt = Some(self.nodes[id].keys[idx + 1].clone()); + } + if let Some((md, child)) = new_child { + // A full node splits first; the child goes to the half that + // holds child `idx`. + let (mut target, mut split) = (id, None); + if n == self.two_k { + let s = self.split(id, idx); + let nleft = self.nodes[id].children.len(); + if idx >= nleft { + idx -= nleft; + target = s; + } + split = Some(s); + } + // H5B__insert_child (H5B_INS_RIGHT): the new child after child + // `idx`, its left key after that child's. + let node = &mut self.nodes[target]; + node.keys.insert(idx + 1, md); + node.children.insert(idx + 1, child); + ret.split = split.map(|s| (self.nodes[s].keys[0].clone(), s)); + } + Ok(ret) + } + + /// `H5B__split` of the full node `id`, the insertion going after child + /// `idx`; returns the new right node. + fn split(&mut self, id: usize, idx: usize) -> usize { + let node = &self.nodes[id]; + let ratio = if node.right.is_none() { + SPLIT_RATIOS[2] + } else if node.left.is_none() { + SPLIT_RATIOS[0] + } else { + SPLIT_RATIOS[1] + }; + let mut nleft = (self.two_k as f64 * ratio) as usize; + if idx < nleft && nleft == self.two_k { + nleft -= 1; + } else if idx >= nleft && nleft == 0 { + nleft += 1; + } + let new_id = self.nodes.len(); + let right = Node { + level: node.level, + left: Some(id), + right: node.right, + keys: node.keys[nleft..].to_vec(), + children: node.children[nleft..].to_vec(), + }; + let old_right = node.right; + self.nodes.push(right); + if let Some(r) = old_right { + self.nodes[r].left = Some(new_id); + } + let node = &mut self.nodes[id]; + node.keys.truncate(nleft + 1); + node.children.truncate(nleft); + node.right = Some(new_id); + new_id + } +} + +/// Bytes of one node of a chunk B-tree with `ndims` key dimensions (the +/// dataset's rank plus the element-size one). +fn node_size(two_k: usize, ndims: usize, offset_size: usize) -> usize { + let key = 8 + 8 * ndims; + 8 + 2 * offset_size + (two_k + 1) * key + two_k * offset_size +} + +/// Build the chunk B-tree for `chunks`, given in row-major order of their +/// scaled coordinates, with nodes laid out from `base_address`. `chunk_dims` +/// are the chunk's dimensions (the dataset's rank of them) and `elem_size` +/// the element size, the key's last dimension. Returns the nodes' bytes; the +/// root is at `base_address`. `chunks` must not be empty: an index without +/// chunks has no tree (its address is undefined). +pub(crate) fn build_chunk_btree_v1_at( + chunks: &[ChunkEntry], + chunk_dims: &[u64], + elem_size: u32, + base_address: u64, + offset_size: u8, +) -> Result, FormatError> { + if chunks.is_empty() { + return Err(bad("no chunks")); + } + let rank = chunk_dims.len(); + let mut tree = Tree::new(CHUNK_BTREE_K); + for c in chunks { + if c.scaled.len() != rank { + return Err(bad("chunk rank differs from the dataset's")); + } + let nbytes = u32::try_from(c.nbytes).map_err(|_| { + FormatError::SerializationError(format!( + "a chunk of {} bytes cannot be indexed by a version-1 B-tree \ + (HDF5 1.8 chunks are under 4 GiB)", + c.nbytes + )) + })?; + let mut scaled = c.scaled.clone(); + scaled.push(0); + let key = Key { + nbytes, + mask: c.filter_mask, + scaled, + }; + tree.insert(&key, c.address)?; + } + + let os = usize::from(offset_size); + let ndims = rank + 1; + let nsize = node_size(tree.two_k, ndims, os); + let addr_of = |id: usize| base_address + (id * nsize) as u64; + let mut dims: Vec = chunk_dims.to_vec(); + dims.push(u64::from(elem_size)); + let mut out = vec![0u8; tree.nodes.len() * nsize]; + for (i, node) in tree.nodes.iter().enumerate() { + let d = &mut out[i * nsize..(i + 1) * nsize]; + d[0..4].copy_from_slice(b"TREE"); + d[4] = 1; // node type: raw data chunks + d[5] = node.level; + let n = u16::try_from(node.children.len()).map_err(|_| bad("node too large"))?; + d[6..8].copy_from_slice(&n.to_le_bytes()); + let undef = u64::MAX; + put_addr(&mut d[8..], node.left.map_or(undef, addr_of), os); + put_addr(&mut d[8 + os..], node.right.map_or(undef, addr_of), os); + let mut p = 8 + 2 * os; + for (k, key) in node.keys.iter().enumerate() { + d[p..p + 4].copy_from_slice(&key.nbytes.to_le_bytes()); + d[p + 4..p + 8].copy_from_slice(&key.mask.to_le_bytes()); + for (j, (&s, &dim)) in key.scaled.iter().zip(&dims).enumerate() { + let off = s + .checked_mul(dim) + .ok_or_else(|| FormatError::Overflow("chunk key offset".into()))?; + d[p + 8 + 8 * j..p + 16 + 8 * j].copy_from_slice(&off.to_le_bytes()); + } + p += 8 + 8 * ndims; + if let Some(&child) = node.children.get(k) { + let a = if node.level == 0 { + child + } else { + addr_of(usize::try_from(child).map_err(|_| bad("bad node index"))?) + }; + put_addr(&mut d[p..], a, os); + p += os; + } + } + } + Ok(out) +} + +fn put_addr(d: &mut [u8], v: u64, os: usize) { + d[..os].copy_from_slice(&v.to_le_bytes()[..os]); +} + +#[cfg(test)] +mod tests { + use super::*; + + fn build(n: u64) -> Tree { + let mut t = Tree::new(CHUNK_BTREE_K); + for i in 0..n { + let key = Key { + nbytes: 80, + mask: 0, + scaled: vec![i, 0], + }; + t.insert(&key, 1000 + i).unwrap(); + } + t + } + + /// Leaves in order from the root, with their child counts. + fn leaves(t: &Tree, id: usize, out: &mut Vec) { + let n = &t.nodes[id]; + if n.level == 0 { + out.push(n.children.len()); + } else { + for &c in &n.children { + leaves(t, c as usize, out); + } + } + } + + #[test] + fn sequential_appends_split_as_libhdf5_does() { + // libhdf5 2.0 (h5py, libver=('v108', 'latest')) writes 1000 chunks + // as a root over 17 leaves of 57 chunks and one of 31, with the + // root's right key on the last chunk (9990, 8 for 10-element f8 + // chunks). + let t = build(1000); + assert_eq!(t.nodes[0].level, 1); + let mut l = Vec::new(); + leaves(&t, 0, &mut l); + let mut want = vec![57; 17]; + want.push(31); + assert_eq!(l, want); + assert_eq!(t.nodes[0].keys.last().unwrap().scaled, vec![999, 1]); + // 100 000 chunks: three levels, a root of 31 children. + let t = build(100_000); + assert_eq!(t.nodes[0].level, 2); + assert_eq!(t.nodes[0].children.len(), 31); + } + + #[test] + fn right_key_moves_every_other_append() { + let t = build(5); + assert_eq!(t.nodes[0].keys.last().unwrap().scaled, vec![5, 1]); + let t = build(6); + assert_eq!(t.nodes[0].keys.last().unwrap().scaled, vec![5, 1]); + } + + #[test] + fn keys_and_siblings_are_consistent() { + let t = build(5000); + for (i, n) in t.nodes.iter().enumerate() { + assert!(n.children.len() <= t.two_k); + assert_eq!(n.keys.len(), n.children.len() + 1); + if let Some(r) = n.right { + assert_eq!(t.nodes[r].left, Some(i)); + assert_eq!(n.keys.last(), t.nodes[r].keys.first()); + } + } + } +} diff --git a/crates/clawhdf5-format/src/chunk_cache.rs b/crates/clawhdf5-format/src/chunk_cache.rs index 6e3fee8..cca133d 100644 --- a/crates/clawhdf5-format/src/chunk_cache.rs +++ b/crates/clawhdf5-format/src/chunk_cache.rs @@ -899,7 +899,7 @@ mod tests { fn make_chunk(offsets: Vec, address: u64, size: u32) -> ChunkInfo { ChunkInfo { - chunk_size: size, + chunk_size: u64::from(size), filter_mask: 0, offsets, address, diff --git a/crates/clawhdf5-format/src/chunk_index.rs b/crates/clawhdf5-format/src/chunk_index.rs index 706bad0..1e083a0 100644 --- a/crates/clawhdf5-format/src/chunk_index.rs +++ b/crates/clawhdf5-format/src/chunk_index.rs @@ -109,7 +109,7 @@ pub struct ChunkMapping { /// File byte address of the compressed chunk. pub file_offset: u64, /// Size of the compressed chunk in the file. - pub file_size: u32, + pub file_size: u64, /// Filter mask (0 = all filters applied). pub filter_mask: u32, /// Pre-computed row-copy operations for assembling this chunk into output. @@ -338,7 +338,7 @@ mod tests { fn make_chunk(offsets: Vec, address: u64, size: u32) -> ChunkInfo { ChunkInfo { - chunk_size: size, + chunk_size: u64::from(size), filter_mask: 0, offsets, address, diff --git a/crates/clawhdf5-format/src/chunked_read.rs b/crates/clawhdf5-format/src/chunked_read.rs index 4d97adb..5f0528a 100644 --- a/crates/clawhdf5-format/src/chunked_read.rs +++ b/crates/clawhdf5-format/src/chunked_read.rs @@ -314,7 +314,7 @@ pub(crate) fn chunk_req( chunk_bytes: usize, wanted: bool, ) -> ExtentReq { - let len = c.chunk_size as usize; + let len = crate::addr::saturating_usize(c.chunk_size); ExtentReq { addr: c.address, len, @@ -521,7 +521,7 @@ pub fn decompress_all_chunks_with_stats_in( #[derive(Debug, Clone)] pub struct ChunkInfo { /// Size of chunk data in the file (after compression). - pub chunk_size: u32, + pub chunk_size: u64, /// Bitmask of filters that were NOT applied (0 = all applied). pub filter_mask: u32, /// N-dimensional offset of this chunk in dataset space. @@ -610,7 +610,7 @@ fn stored_element_size(dt: &Datatype, offset_size: u8) -> u64 { /// (`H5D__chunk_set_sizes`: "stored datatype size in chunk layout does not /// match datatype description"). Reading it anyway laid the chunks out with /// the wrong element size. -pub(crate) fn check_chunk_element_size( +pub fn check_chunk_element_size( layout: &DataLayout, datatype: &Datatype, offset_size: u8, @@ -1028,7 +1028,7 @@ fn parse_chunk_node( if node_level == 0 { chunks.push(stored.len()); stored.push(ChunkInfo { - chunk_size, + chunk_size: u64::from(chunk_size), filter_mask, offsets: keys[k..].to_vec(), address, @@ -1131,7 +1131,7 @@ pub fn generate_implicit_chunks_in_grid( } chunks.push(ChunkInfo { - chunk_size: chunk_byte_size as u32, + chunk_size: chunk_byte_size, filter_mask: 0, offsets, address: base_address.saturating_add(grid_idx.saturating_mul(chunk_byte_size)), @@ -1190,9 +1190,15 @@ fn read_btree_v2_chunks( } _ => return Err(bad("tree is not a chunk index")), }; - let unfiltered_bytes = checked_chunk_byte_len(chunk_dims, elem_size)?; - let unfiltered_bytes = - u32::try_from(unfiltered_bytes).map_err(|_| bad("chunk larger than 4 GiB"))?; + // A u64 whatever the platform: an unfiltered chunk's size is only + // recorded here; a chunk this platform cannot address fails when read. + let unfiltered_bytes = u64::try_from( + chunk_dims + .iter() + .try_fold(elem_size as u128, |acc, &c| acc.checked_mul(c as u128)) + .ok_or_else(|| bad("chunk size overflows"))?, + ) + .map_err(|_| bad("chunk size overflows"))?; let records = collect_btree_v2_records_in(file_data, &header, offset_size, length_size)?; let mut chunks = Vec::with_capacity(records.len()); @@ -1213,10 +1219,7 @@ fn read_btree_v2_chunks( pos += size_len; let mask = u32::from_le_bytes([data[pos], data[pos + 1], data[pos + 2], data[pos + 3]]); pos += 4; - ( - u32::try_from(size).map_err(|_| bad("stored chunk larger than 4 GiB"))?, - mask, - ) + (size, mask) }; let mut offsets = Vec::with_capacity(rank); for &dim in chunk_dims { @@ -1334,9 +1337,9 @@ pub fn list_chunks_in( // Single chunk — one chunk covering the entire dataset let chunk_byte_size = checked_chunk_byte_len(&chunk_dims, elem_size)?; let (csize, fmask) = if let Some(fs) = single_filtered_size { - (fs as u32, single_filter_mask.unwrap_or(0)) + (fs, single_filter_mask.unwrap_or(0)) } else { - (chunk_byte_size as u32, 0) + (chunk_byte_size as u64, 0) }; vec![ChunkInfo { chunk_size: csize, @@ -1488,7 +1491,7 @@ pub fn list_chunks_for_read_in( let chunk_bytes = checked_chunk_byte_len(&chunk_dims, elem_size)?; if let Some(c) = chunks .iter() - .find(|c| c.address != u64::MAX && c.chunk_size as usize != chunk_bytes) + .find(|c| c.address != u64::MAX && c.chunk_size != chunk_bytes as u64) { return Err(FormatError::ChunkedReadError(format!( "incorrect chunk size returned from index for unfiltered chunk at {:?}: \ @@ -2124,7 +2127,7 @@ pub fn read_chunked_data_indexed_in( .iter() .zip(&hits) .map(|(m, hit)| { - let len = m.file_size as usize; + let len = crate::addr::saturating_usize(m.file_size); ExtentReq { addr: m.file_offset, len, @@ -2467,7 +2470,7 @@ mod tests { // Entries: key[i], child[i] pairs, then final key for chunk in chunks { // Key: chunk_size(4) + filter_mask(4) + ndims offsets - buf.extend_from_slice(&chunk.chunk_size.to_le_bytes()); + buf.extend_from_slice(&(chunk.chunk_size as u32).to_le_bytes()); buf.extend_from_slice(&chunk.filter_mask.to_le_bytes()); for d in 0..ndims { let off = if d < chunk.offsets.len() { @@ -2757,7 +2760,7 @@ mod tests { } chunk_infos.push(ChunkInfo { - chunk_size: chunk_bytes as u32, + chunk_size: chunk_bytes as u64, filter_mask: 0, offsets: vec![start as u64, 0], address: data_offset as u64, @@ -2869,7 +2872,7 @@ mod tests { .collect(); let stored = crate::filters::compress_chunk(&chunk, &pipeline, 4).unwrap(); chunks.push(ChunkInfo { - chunk_size: stored.len() as u32, + chunk_size: stored.len() as u64, filter_mask: 0, offsets: vec![r0 as u64, c0 as u64, 0], address: file.len() as u64, @@ -2943,7 +2946,7 @@ mod tests { let short = crate::filters::compress_chunk(&[1u8; 64], &pipeline, 4).unwrap(); for bad in [5usize, 11, 40] { chunks[bad].address = file.len() as u64; - chunks[bad].chunk_size = short.len() as u32; + chunks[bad].chunk_size = short.len() as u64; file.extend_from_slice(&short); } for _ in 0..20 { @@ -3133,7 +3136,7 @@ mod tests { file_data[data_offset..data_offset + compressed.len()].copy_from_slice(&compressed); chunk_infos.push(ChunkInfo { - chunk_size: compressed.len() as u32, + chunk_size: compressed.len() as u64, filter_mask: 0, offsets: vec![start as u64, 0], address: data_offset as u64, @@ -3215,7 +3218,7 @@ mod tests { file_data[data_offset..data_offset + chunk_size].copy_from_slice(&chunk_bytes); chunk_infos.push(ChunkInfo { - chunk_size: chunk_size as u32, + chunk_size: chunk_size as u64, filter_mask: 0, offsets: vec![row_start as u64, col_start as u64, 0], address: data_offset as u64, @@ -3340,7 +3343,7 @@ mod tests { assert_eq!(c.address, 0x1000 + i as u64 * chunk_byte_size as u64); assert_eq!(c.offsets, vec![i as u64 * 20]); assert_eq!(c.filter_mask, 0); - assert_eq!(c.chunk_size, chunk_byte_size as u32); + assert_eq!(c.chunk_size, chunk_byte_size as u64); } } diff --git a/crates/clawhdf5-format/src/chunked_write.rs b/crates/clawhdf5-format/src/chunked_write.rs index 78c4161..d750fd8 100644 --- a/crates/clawhdf5-format/src/chunked_write.rs +++ b/crates/clawhdf5-format/src/chunked_write.rs @@ -7,6 +7,7 @@ use crate::addr::saturating_usize; #[cfg(not(feature = "std"))] use alloc::{format, vec, vec::Vec}; +use crate::btree_v1_write; use crate::btree_v2_write::{BTreeV2Params, build_btree_v2}; use crate::checksum::jenkins_lookup3; use crate::chunk_cache::{CACHE_LINE_SIZE, align_to_cache_line}; @@ -19,6 +20,7 @@ use crate::filter_pipeline::{ FilterPipeline, }; use crate::filters::compress_chunk_masked; +use crate::libver::LibVer; /// Round a file offset up to the next cache-line boundary. /// /// This ensures chunk data starts at an address that is a multiple of the @@ -400,86 +402,83 @@ pub fn split_into_chunks( chunk_dims: &[u64], element_size: usize, ) -> Vec<(Vec, Vec)> { - let rank = shape.len(); - if rank == 0 { + if shape.is_empty() { return vec![(vec![], raw_data.to_vec())]; } + (0..chunk_count(shape, chunk_dims)) + .map(|i| extract_chunk(raw_data, shape, chunk_dims, element_size, i)) + .collect() +} - // Compute number of chunks per dimension - let mut num_chunks_per_dim = Vec::with_capacity(rank); - for d in 0..rank { - num_chunks_per_dim.push(shape[d].div_ceil(chunk_dims[d])); +/// Number of chunks of the current extent `shape`. +fn chunk_count(shape: &[u64], chunk_dims: &[u64]) -> u64 { + shape + .iter() + .zip(chunk_dims) + .map(|(&s, &c)| s.div_ceil(c)) + .product() +} + +/// The `linear_idx`-th chunk (row-major over the chunks of the current +/// extent) of the row-major dataset `raw_data`: its offset in dataset space +/// and its bytes, a whole chunk with the part past the dataset's edge zero. +/// Copied one row (a run along the last dimension) at a time. +fn extract_chunk( + raw_data: &[u8], + shape: &[u64], + chunk_dims: &[u64], + element_size: usize, + linear_idx: u64, +) -> (Vec, Vec) { + let rank = shape.len(); + let mut offsets = vec![0u64; rank]; + let mut remaining = linear_idx; + for d in (0..rank).rev() { + let n = shape[d].div_ceil(chunk_dims[d]); + offsets[d] = (remaining % n) * chunk_dims[d]; + remaining /= n; } - let total_chunks: u64 = num_chunks_per_dim.iter().product(); + let chunk_elements: usize = chunk_dims.iter().map(|&d| saturating_usize(d)).product(); + let mut chunk = vec![0u8; chunk_elements.saturating_mul(element_size)]; - // Dataset strides (row-major) + // Elements of the chunk inside the dataset, per dimension. + let valid: Vec = (0..rank) + .map(|d| saturating_usize(shape[d].saturating_sub(offsets[d]).min(chunk_dims[d]))) + .collect(); + if valid.contains(&0) { + return (offsets, chunk); + } let mut ds_strides = vec![1usize; rank]; - for i in (0..rank.saturating_sub(1)).rev() { - ds_strides[i] = ds_strides[i + 1] * saturating_usize(shape[i + 1]); - } - - // Chunk strides let mut chunk_strides = vec![1usize; rank]; - for i in (0..rank.saturating_sub(1)).rev() { - chunk_strides[i] = chunk_strides[i + 1] * saturating_usize(chunk_dims[i + 1]); + for d in (0..rank - 1).rev() { + ds_strides[d] = ds_strides[d + 1] * saturating_usize(shape[d + 1]); + chunk_strides[d] = chunk_strides[d + 1] * saturating_usize(chunk_dims[d + 1]); } - - let chunk_total_elements: usize = chunk_dims.iter().map(|&d| saturating_usize(d)).product(); - - let mut result = Vec::with_capacity(saturating_usize(total_chunks)); - - for linear_idx in 0..total_chunks { - // Convert linear index to chunk grid coordinates - let mut chunk_grid_coords = vec![0u64; rank]; - let mut remaining = linear_idx; - for d in (0..rank).rev() { - chunk_grid_coords[d] = remaining % num_chunks_per_dim[d]; - remaining /= num_chunks_per_dim[d]; + let row = valid[rank - 1] * element_size; + let mut idx = vec![0usize; rank]; + loop { + let src: usize = (0..rank) + .map(|d| (saturating_usize(offsets[d]) + idx[d]) * ds_strides[d]) + .sum::() + * element_size; + let dst: usize = (0..rank).map(|d| idx[d] * chunk_strides[d]).sum::() * element_size; + // Whole elements only, as far as `raw_data` reaches. + let n = row.min(raw_data.len().saturating_sub(src)) / element_size * element_size; + chunk[dst..dst + n].copy_from_slice(&raw_data[src..src + n]); + // Next row: advance every dimension but the last. + let mut d = rank - 1; + loop { + if d == 0 { + return (offsets, chunk); + } + d -= 1; + idx[d] += 1; + if idx[d] < valid[d] { + break; + } + idx[d] = 0; } - - // Chunk offset in dataset space - let offsets: Vec = (0..rank) - .map(|d| chunk_grid_coords[d] * chunk_dims[d]) - .collect(); - - // Extract chunk data - let mut chunk_bytes = vec![0u8; chunk_total_elements * element_size]; - - for flat_idx in 0..chunk_total_elements { - let mut remaining_idx = flat_idx; - let mut ds_flat = 0usize; - let mut out_of_bounds = false; - - for d in 0..rank { - let coord_in_chunk = remaining_idx / chunk_strides[d]; - remaining_idx %= chunk_strides[d]; - - let global_coord = saturating_usize(offsets[d]) + coord_in_chunk; - if global_coord >= saturating_usize(shape[d]) { - out_of_bounds = true; - break; - } - ds_flat += global_coord * ds_strides[d]; - } - - if out_of_bounds { - // Zero-filled (already initialized) - continue; - } - - let src_start = ds_flat * element_size; - let dst_start = flat_idx * element_size; - - if src_start + element_size <= raw_data.len() { - chunk_bytes[dst_start..dst_start + element_size] - .copy_from_slice(&raw_data[src_start..src_start + element_size]); - } - } - - result.push((offsets, chunk_bytes)); } - - result } /// Parallel compression threshold: use rayon when chunk count exceeds this. @@ -490,46 +489,68 @@ pub fn split_into_chunks( #[cfg(feature = "parallel")] const PARALLEL_COMPRESS_THRESHOLD: usize = 2; -/// Compress all chunks, using parallel compression when beneficial, and -/// return each chunk's stored bytes with its filter mask. +/// Largest chunk compressed in parallel: every thread holds a chunk and its +/// compressed copy at once, so chunks larger than this (up to 4 GiB and +/// more) are compressed one after another. +#[cfg(feature = "parallel")] +const PARALLEL_COMPRESS_MAX_CHUNK_BYTES: u64 = 64 << 20; + +/// Extract and compress every chunk of the dataset, returning each chunk's +/// raw size, stored bytes and filter mask, in chunk order. /// /// Chunks run through the pipeline as libhdf5 runs them /// ([`compress_chunk_masked`]): an optional filter that fails — LZF or Blosc /// output no smaller than its input — is skipped and its mask bit set. /// -/// With the `parallel` feature and more than [`PARALLEL_COMPRESS_THRESHOLD`] -/// filtered chunks, compression runs across rayon threads; otherwise it is -/// sequential. Output order matches input order, so per-chunk bytes are -/// identical to the sequential path. +/// Each chunk is extracted just before it is compressed and dropped after, +/// so at most one raw chunk per thread is held. With the `parallel` feature, +/// more than [`PARALLEL_COMPRESS_THRESHOLD`] filtered chunks, and chunks of +/// at most [`PARALLEL_COMPRESS_MAX_CHUNK_BYTES`], compression runs across +/// rayon threads; otherwise it is sequential. Output order matches chunk +/// order, so per-chunk bytes are identical to the sequential path. fn compress_all_chunks( - chunks: &[(Vec, Vec)], + raw_data: &[u8], + shape: &[u64], + chunk_dims: &[u64], + element_size: usize, + chunk_bytes: u64, pipeline: &Option, - element_size: u32, -) -> Result, u32)>, FormatError> { +) -> Result, u32)>, FormatError> { + let one = |i: u64| -> Result<(u64, Vec, u32), FormatError> { + let (_, raw) = if shape.is_empty() { + (Vec::new(), raw_data.to_vec()) + } else { + extract_chunk(raw_data, shape, chunk_dims, element_size, i) + }; + let raw_size = raw.len() as u64; + match pipeline { + Some(pl) => { + let (stored, mask) = compress_chunk_masked(&raw, pl, element_size as u32)?; + Ok((raw_size, stored, mask)) + } + None => Ok((raw_size, raw, 0)), + } + }; + let n = if shape.is_empty() { + 1 + } else { + chunk_count(shape, chunk_dims) + }; #[cfg(feature = "parallel")] { - if let Some(pl) = pipeline - && chunks.len() > PARALLEL_COMPRESS_THRESHOLD + if pipeline.is_some() + && n > PARALLEL_COMPRESS_THRESHOLD as u64 + && chunk_bytes <= PARALLEL_COMPRESS_MAX_CHUNK_BYTES { use rayon::prelude::*; - return chunks - .par_iter() - .map(|(_offsets, chunk_bytes)| compress_chunk_masked(chunk_bytes, pl, element_size)) - .collect(); + return (0..n).into_par_iter().map(one).collect(); } } + #[cfg(not(feature = "parallel"))] + let _ = chunk_bytes; // Sequential fallback - chunks - .iter() - .map(|(_offsets, chunk_bytes)| { - if let Some(pl) = pipeline { - compress_chunk_masked(chunk_bytes, pl, element_size) - } else { - Ok((chunk_bytes.clone(), 0)) - } - }) - .collect() + (0..n).map(one).collect() } /// Build the complete chunked dataset blob (chunk data + index) and return @@ -550,10 +571,12 @@ pub fn serialize_v4_single_chunk_pub( filter_mask, offset_size, element_size, + 4, ) } -/// Serialize a v4 single chunk layout message. +/// Serialize a v4 (or, for a chunk of 4 GiB or more, v5) single chunk +/// layout message. fn serialize_v4_single_chunk( chunk_dims: &[u32], chunk_address: u64, @@ -561,9 +584,10 @@ fn serialize_v4_single_chunk( filter_mask: Option, offset_size: u8, element_size: u32, + version: u8, ) -> Vec { let mut buf = Vec::new(); - buf.push(4); // version + buf.push(version); buf.push(2); // class = chunked // flags: bit 0 = unknown meaning in some files, bit 1 = filters for single chunk @@ -603,8 +627,9 @@ fn serialize_v4_fixed_array( offset_size: u8, element_size: u32, max_bits: u8, + version: u8, ) -> Vec { - let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size); + let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size, version); // chunk index type = 3 (Fixed Array) buf.push(3); @@ -643,9 +668,9 @@ pub(crate) fn push_v4_chunk_dims(buf: &mut Vec, chunk_dims: &[u32], element_ } } -fn layout_v4_chunked_prefix(chunk_dims: &[u32], element_size: u32) -> Vec { +fn layout_v4_chunked_prefix(chunk_dims: &[u32], element_size: u32, version: u8) -> Vec { let mut buf = Vec::new(); - buf.push(4); // version + buf.push(version); buf.push(2); // class = chunked let flags: u8 = 0x00; @@ -671,21 +696,53 @@ pub(crate) fn push_addr(buf: &mut Vec, addr: u64, offset_size: u8) { /// Width of the chunk-size field of a filtered chunk index element. Must /// match the library's `H5D_FARRAY_FILT_COMPUTE_CHUNK_SIZE_LEN` (the EA and -/// B-tree v2 indexes use the same formula): -/// `1 + ((log2(unfiltered chunk bytes) + 8) / 8)`, capped at 8. -pub(crate) fn filtered_chunk_size_len(slots: &[Option]) -> usize { +/// B-tree v2 indexes use the same formula): see [`chunk_size_len`]. Chunks +/// of more than `u32::MAX` bytes are written with layout version 5 +/// ([`layout_version_for`]). +pub(crate) fn filtered_chunk_size_len(slots: &[Option], length_size: u8) -> usize { let max_raw = slots .iter() .flatten() .map(|c| c.raw_size) .max() .unwrap_or(1); - let log2_val = if max_raw <= 1 { + chunk_size_len(max_raw, layout_version_for(max_raw), length_size) +} + +/// Largest chunk, in bytes, a layout message of version 4 or lower may +/// describe: libhdf5 writes a larger one with version 5 +/// (`H5D__chunk_construct`: "chunk size > 4GB requires H5F_LIBVER_V200"), +/// which libhdf5 before 2.0 cannot read. +pub const MAX_V4_CHUNK_BYTES: u64 = u32::MAX as u64; + +/// The layout message version clawhdf5 writes for chunks of `chunk_bytes` +/// bytes: 4, or 5 for a chunk larger than [`MAX_V4_CHUNK_BYTES`] (what +/// libhdf5 2.x writes for it; the chunk index is chosen as for version 4). +pub fn layout_version_for(chunk_bytes: u64) -> u8 { + if chunk_bytes > MAX_V4_CHUNK_BYTES { + 5 + } else { + 4 + } +} + +/// Width libhdf5 gives the stored-size field of a filtered chunk index +/// element (Fixed Array, Extensible Array, v2 B-tree) for chunks of +/// `chunk_bytes` bytes under layout message `layout_version` +/// (`H5D_FARRAY_FILT_COMPUTE_CHUNK_SIZE_LEN` and its EA and B-tree twins): +/// up to version 4, one byte more than the chunk's size needs, +/// `1 + ((log2(chunk_bytes) + 8) / 8)` capped at 8; from version 5 (HDF5 +/// 2.0), the file's size of lengths (`length_size`), whatever the chunk. +pub fn chunk_size_len(chunk_bytes: u64, layout_version: u8, length_size: u8) -> usize { + if layout_version >= 5 { + return usize::from(length_size); + } + let log2 = if chunk_bytes <= 1 { 0 } else { - 63 - max_raw.leading_zeros() + 63 - chunk_bytes.leading_zeros() }; - (1 + ((log2_val + 8) / 8) as usize).min(8) + (1 + ((log2 + 8) / 8) as usize).min(8) } /// Append one chunk index element: the chunk's address, plus its stored size @@ -732,7 +789,7 @@ pub fn build_fixed_array_at( let os = offset_size as usize; let num_elements = slots.len(); - let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots)); + let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots, length_size)); let elem_size = os + chunk_size_bytes.map_or(0, |n| n + 4); let client_id: u8 = if has_filters { 1 } else { 0 }; @@ -827,23 +884,23 @@ pub fn precompress_chunks( element_size: usize, options: &ChunkOptions, ) -> Result { - let chunk_bytes = chunk_dims - .iter() - .try_fold(element_size as u64, |acc, &d| acc.checked_mul(d)) - .and_then(|b| u32::try_from(b).ok()) - .unwrap_or(0); - let pipeline = options.build_pipeline_for_chunk(element_size as u32, chunk_bytes); + let (_, chunk_bytes) = checked_chunk_dims(chunk_dims, element_size)?; + if chunk_bytes > MAX_V4_CHUNK_BYTES { + check_huge_chunk_filters(options, chunk_bytes)?; + } + let pipeline = options + .build_pipeline_for_chunk(element_size as u32, u32::try_from(chunk_bytes).unwrap_or(0)); let has_filters = pipeline.is_some(); let pipeline_message = pipeline.as_ref().map(|pl| pl.serialize()); - let raw_chunks = split_into_chunks(raw_data, shape, chunk_dims, element_size); - let compressed = compress_all_chunks(&raw_chunks, &pipeline, element_size as u32)?; - - let chunks = raw_chunks - .into_iter() - .zip(compressed) - .map(|((_offsets, raw_bytes), (c, mask))| (raw_bytes.len() as u64, c, mask)) - .collect(); + let chunks = compress_all_chunks( + raw_data, + shape, + chunk_dims, + element_size, + chunk_bytes, + &pipeline, + )?; Ok(PrecompressedChunks { chunks, @@ -855,6 +912,67 @@ pub fn precompress_chunks( }) } +/// The chunk dimensions as the layout message stores them (each below +/// 2^32), and one chunk's size in bytes. A chunk dimension of 2^32 or more +/// (which HDF5 2.0 can store, in wider fields) is refused, as it is when +/// read: clawhdf5 holds chunk dimensions as `u32`. So is a chunk whose size +/// overflows 64 bits, or that this platform cannot hold in memory (a chunk +/// of 4 GiB or more on a 32-bit target). +fn checked_chunk_dims( + chunk_dims: &[u64], + element_size: usize, +) -> Result<(Vec, u64), FormatError> { + let dims = chunk_dims + .iter() + .map(|&d| { + u32::try_from(d).map_err(|_| { + FormatError::InvalidChunkDimensions(format!( + "chunk dimension {d} is 2^32 or more, which clawhdf5 does not support" + )) + }) + }) + .collect::, _>>()?; + let bytes = chunk_dims + .iter() + .try_fold(element_size as u64, |acc, &d| acc.checked_mul(d)) + .filter(|&b| usize::try_from(b).is_ok()) + .ok_or_else(|| { + FormatError::Overflow(format!( + "a chunk of {chunk_dims:?} x {element_size} bytes exceeds this platform's address space" + )) + })?; + Ok((dims, bytes)) +} + +/// The filters clawhdf5 can apply to a chunk of more than `u32::MAX` bytes: +/// shuffle, deflate, Zstandard, LZ4 (whose HDF5 framing records the size in +/// 64 bits and splits the chunk into blocks) and Fletcher32. The others +/// record the chunk size or their block lengths in 32 bits, or cannot take +/// a buffer that large (h5py's LZF, bitshuffle, bzip2, Blosc), and pcodec is +/// clawhdf5's own; they are refused rather than written into a chunk +/// libhdf5 could not decode. +fn check_huge_chunk_filters(options: &ChunkOptions, chunk_bytes: u64) -> Result<(), FormatError> { + let refused = if let Some(plugin) = &options.plugin { + Some(match plugin { + PluginFilter::Lzf => "LZF", + PluginFilter::Bitshuffle { .. } => "bitshuffle", + PluginFilter::Bzip2 { .. } => "bzip2", + PluginFilter::Blosc { .. } => "Blosc", + }) + } else if options.pcodec { + Some("pcodec") + } else { + None + }; + match refused { + Some(name) => Err(FormatError::FilterError(format!( + "{name} cannot compress a chunk of {chunk_bytes} bytes (4 GiB or more); \ + use smaller chunks, or deflate, Zstandard or LZ4" + ))), + None => Ok(()), + } +} + /// Lay out precompressed chunks at `base_address` and build index structures. /// /// This is the address-dependent half of chunk writing. Call it in Pass 1 @@ -866,6 +984,44 @@ pub fn build_chunked_data_from_precompressed( base_address: u64, maxshape: Option<&[u64]>, ) -> Result { + build_chunked_data_from_precompressed_libver( + pre, + base_address, + maxshape, + LibVer::Latest, + LibVer::Latest, + ) +} + +/// [`build_chunked_data_from_precompressed`] for a file whose low library +/// version bound is `low`: below [`LibVer::V110`] (that is, for HDF5 1.8) +/// every chunked dataset gets a version-3 layout message and a version-1 +/// B-tree chunk index, whatever its maximum shape, as libhdf5 writes it; +/// otherwise the version-4 layout and the index libhdf5 picks for it. +/// +/// A chunk of 4 GiB or more (over [`MAX_V4_CHUNK_BYTES`]) takes layout +/// message version 5 whatever `low` is, as in libhdf5 (a version-1 B-tree +/// key holds a 32-bit size), and needs a `high` bound of at least +/// [`LibVer::V200`] ([`FormatError::LibverBound`] otherwise). +pub fn build_chunked_data_from_precompressed_libver( + pre: &PrecompressedChunks, + base_address: u64, + maxshape: Option<&[u64]>, + low: LibVer, + high: LibVer, +) -> Result { + let (_, chunk_bytes) = checked_chunk_dims(&pre.chunk_dims, pre.element_size)?; + if chunk_bytes > MAX_V4_CHUNK_BYTES { + if high < LibVer::V200 { + return Err(FormatError::LibverBound { + what: format!("a chunk of {chunk_bytes} bytes (4 GiB or more)"), + needs: LibVer::V200, + high, + }); + } + } else if low < LibVer::V110 { + return build_btree_v1_chunked_data(pre, base_address, maxshape); + } let index = ChunkIndexPlan::new(&pre.shape, maxshape, &pre.chunk_dims)?; let offset_size: u8 = 8; let length_size: u8 = 8; @@ -891,7 +1047,8 @@ pub fn build_chunked_data_from_precompressed( }); } - let chunk_dims_u32: Vec = pre.chunk_dims.iter().map(|&d| d as u32).collect(); + let (chunk_dims_u32, chunk_bytes) = checked_chunk_dims(&pre.chunk_dims, element_size)?; + let version = layout_version_for(chunk_bytes); let aligned_idx = align_to_cache_line(data_buf.len()); if aligned_idx > data_buf.len() { @@ -915,6 +1072,7 @@ pub fn build_chunked_data_from_precompressed( ea_address, offset_size, element_size as u32, + version, ) } ChunkIndexPlan::SingleChunk => { @@ -932,6 +1090,7 @@ pub fn build_chunked_data_from_precompressed( filter_mask, offset_size, element_size as u32, + version, ) } ChunkIndexPlan::FixedArray(grid, nslots) => { @@ -957,6 +1116,7 @@ pub fn build_chunked_data_from_precompressed( offset_size, element_size as u32, FA_PAGE_BITS, + version, ) } ChunkIndexPlan::BTreeV2 => { @@ -981,6 +1141,7 @@ pub fn build_chunked_data_from_precompressed( offset_size, element_size as u32, node_size, + version, ) } }; @@ -992,6 +1153,92 @@ pub fn build_chunked_data_from_precompressed( }) } +/// Lay out precompressed chunks at `base_address` followed by a version-1 +/// B-tree chunk index, with a version-3 layout message: what libhdf5 writes +/// for a chunked dataset under a low bound of 1.8. +fn build_btree_v1_chunked_data( + pre: &PrecompressedChunks, + base_address: u64, + maxshape: Option<&[u64]>, +) -> Result { + if let Some(ms) = maxshape { + let bad = |what: &str| FormatError::ChunkedReadError(format!("maxshape: {what}")); + if ms.len() != pre.shape.len() { + return Err(bad("rank differs from the shape")); + } + if ms.iter().zip(&pre.shape).any(|(&m, &s)| m < s) { + return Err(bad("smaller than the shape")); + } + } + let offset_size: u8 = 8; + let mut data_buf = Vec::new(); + let mut entries = Vec::with_capacity(pre.chunks.len()); + for (i, (_raw_size, stored, filter_mask)) in pre.chunks.iter().enumerate() { + let aligned_offset = align_to_cache_line(data_buf.len()); + if aligned_offset > data_buf.len() { + data_buf.resize(aligned_offset, 0u8); + } + entries.push(btree_v1_write::ChunkEntry { + scaled: scaled_coords(&pre.shape, &pre.chunk_dims, i), + nbytes: stored.len() as u64, + filter_mask: *filter_mask, + address: base_address + data_buf.len() as u64, + }); + data_buf.extend_from_slice(stored); + } + let element_size = u32::try_from(pre.element_size) + .map_err(|_| FormatError::Overflow("element size".into()))?; + // A dataset with no chunks has no tree: its address is undefined, as + // libhdf5 leaves it until the first chunk is written. + let btree_address = if entries.is_empty() { + u64::MAX + } else { + let aligned_idx = align_to_cache_line(data_buf.len()); + if aligned_idx > data_buf.len() { + data_buf.resize(aligned_idx, 0u8); + } + let addr = base_address + data_buf.len() as u64; + let tree = btree_v1_write::build_chunk_btree_v1_at( + &entries, + &pre.chunk_dims, + element_size, + addr, + offset_size, + )?; + data_buf.extend_from_slice(&tree); + addr + }; + let layout_message = + serialize_v3_chunked(&pre.chunk_dims, btree_address, offset_size, element_size)?; + Ok(ChunkedDataResult { + data_bytes: data_buf, + layout_message, + pipeline_message: pre.pipeline_message.clone(), + }) +} + +/// A version-3 layout message for a chunked dataset: dimensionality (the +/// rank plus one), the B-tree's address, then each chunk dimension and the +/// element size, four bytes each. +fn serialize_v3_chunked( + chunk_dims: &[u64], + btree_address: u64, + offset_size: u8, + element_size: u32, +) -> Result, FormatError> { + let ndims = u8::try_from(chunk_dims.len() + 1) + .map_err(|_| FormatError::Overflow("chunked layout rank".into()))?; + let mut buf = vec![3u8, 2, ndims]; + push_addr(&mut buf, btree_address, offset_size); + for &d in chunk_dims { + let d = + u32::try_from(d).map_err(|_| FormatError::Overflow(format!("chunk dimension {d}")))?; + buf.extend_from_slice(&d.to_le_bytes()); + } + buf.extend_from_slice(&element_size.to_le_bytes()); + Ok(buf) +} + /// Most slots a Fixed Array index may have before we refuse to build it: its /// data block holds one element per chunk of the *maximum* extent, so a huge /// finite maxshape with small chunks would otherwise exhaust memory. @@ -1132,7 +1379,7 @@ fn build_btree_v2_chunk_index_at( let chunk_size_bytes = has_filters.then(|| { let slots: Vec> = records.iter().map(|(_, c)| Some((*c).clone())).collect(); - filtered_chunk_size_len(&slots) + filtered_chunk_size_len(&slots, length_size) }); let record_size = os + chunk_size_bytes.map_or(0, |n| n + 4) + 8 * rank; let record_size_u16 = u16::try_from(record_size) @@ -1182,8 +1429,9 @@ fn serialize_v4_btree_v2( offset_size: u8, element_size: u32, node_size: u32, + version: u8, ) -> Vec { - let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size); + let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size, version); buf.push(5); // chunk index type = 5 (version-2 B-tree) buf.extend_from_slice(&node_size.to_le_bytes()); buf.push(BT2_SPLIT_PERCENT); @@ -1763,7 +2011,7 @@ mod tests { #[test] fn serialize_v4_single_chunk_no_filters_roundtrip() { - let msg = serialize_v4_single_chunk(&[20], 0x1000, None, None, 8, 8); + let msg = serialize_v4_single_chunk(&[20], 0x1000, None, None, 8, 8, 4); let layout = DataLayout::parse(&msg, 8, 8).unwrap(); match layout { DataLayout::Chunked { @@ -1788,7 +2036,7 @@ mod tests { #[test] fn serialize_v4_single_chunk_with_filters_roundtrip() { - let msg = serialize_v4_single_chunk(&[100], 0x2000, Some(500), Some(0), 8, 8); + let msg = serialize_v4_single_chunk(&[100], 0x2000, Some(500), Some(0), 8, 8, 4); let layout = DataLayout::parse(&msg, 8, 8).unwrap(); match layout { DataLayout::Chunked { @@ -1807,7 +2055,7 @@ mod tests { #[test] fn serialize_v4_fixed_array_roundtrip() { - let msg = serialize_v4_fixed_array(&[20], 0x3000, 8, 8, 4); + let msg = serialize_v4_fixed_array(&[20], 0x3000, 8, 8, 4, 4); let layout = DataLayout::parse(&msg, 8, 8).unwrap(); match layout { DataLayout::Chunked { @@ -1851,11 +2099,245 @@ mod tests { assert_eq!(&fa[28..32], b"FADB"); } + // ---- Chunks of 4 GiB or more (layout message version 5) ---- + + /// 2^29 + 1 `f64`: 4 GiB + 8 bytes, the smallest `f64` chunk past + /// `u32::MAX`. + const HUGE_DIM: u64 = (1 << 29) + 1; + const HUGE_BYTES: u64 = HUGE_DIM * 8; + + fn huge_chunk(address: u64, compressed_size: u64) -> WrittenChunk { + WrittenChunk { + address, + compressed_size, + raw_size: HUGE_BYTES, + filter_mask: 0, + } + } + + #[test] + fn chunks_past_u32_max_take_layout_version_5() { + assert_eq!(layout_version_for(0), 4); + assert_eq!(layout_version_for(u64::from(u32::MAX)), 4); + assert_eq!(layout_version_for(u64::from(u32::MAX) + 1), 5); + assert_eq!(layout_version_for(HUGE_BYTES), 5); + } + + /// A chunk of 4 GiB or more never goes into a version-1 B-tree (its key + /// holds a 32-bit size): under a 1.8 low bound it still takes layout + /// version 5, and a high bound below 2.0 refuses it, as in libhdf5. + #[test] + fn huge_chunks_ignore_the_v18_low_bound_and_need_v200() { + // One filtered chunk; the stored bytes stand in for its compression. + let pre = PrecompressedChunks { + chunks: vec![(HUGE_BYTES, vec![0u8; 16], 0)], + has_filters: true, + element_size: 8, + shape: vec![HUGE_DIM], + chunk_dims: vec![HUGE_DIM], + pipeline_message: None, + }; + let r = build_chunked_data_from_precompressed_libver( + &pre, + 4096, + None, + LibVer::V18, + LibVer::Latest, + ) + .unwrap(); + assert_eq!(r.layout_message[0], 5, "layout message version"); + for high in [LibVer::V18, LibVer::V114] { + match build_chunked_data_from_precompressed_libver(&pre, 4096, None, LibVer::V18, high) + { + Err(FormatError::LibverBound { needs, .. }) => assert_eq!(needs, LibVer::V200), + other => panic!("high bound {high}: {:?}", other.map(|r| r.layout_message)), + } + } + } + + /// `H5D_FARRAY_FILT_COMPUTE_CHUNK_SIZE_LEN` (and its EA and v2 B-tree + /// twins) in libhdf5 2.2.0: one byte more than the chunk size needs up to + /// layout version 4, the size of lengths from version 5. + #[test] + fn chunk_size_len_follows_libhdf5() { + assert_eq!(chunk_size_len(1, 4, 8), 2); + assert_eq!(chunk_size_len(160, 4, 8), 2); + assert_eq!(chunk_size_len(255, 4, 8), 2); + assert_eq!(chunk_size_len(256, 4, 8), 3); + assert_eq!(chunk_size_len(u64::from(u32::MAX), 4, 8), 5); + assert_eq!(chunk_size_len(1 << 32, 4, 8), 6); + assert_eq!(chunk_size_len(u64::MAX, 4, 8), 8); + assert_eq!(chunk_size_len(HUGE_BYTES, 5, 8), 8); + assert_eq!(chunk_size_len(160, 5, 8), 8); + assert_eq!(chunk_size_len(HUGE_BYTES, 5, 4), 4); + let slots = [Some(huge_chunk(0x1000, 20_000)), None]; + assert_eq!(filtered_chunk_size_len(&slots, 8), 8); + let small = [Some(WrittenChunk { + raw_size: 160, + ..huge_chunk(0x1000, 100) + })]; + assert_eq!(filtered_chunk_size_len(&small, 8), 2); + } + + /// A 4 GiB + 8 byte chunk: every layout message is version 5 with the + /// dimensions libhdf5 writes (4 bytes each: 0x20000001 and 8), and a + /// filtered index element stores the chunk's size in 8 bytes. + #[test] + fn huge_chunk_layout_messages_and_index_elements() { + let dims = [HUGE_DIM as u32]; + let parsed = |msg: &[u8]| { + assert_eq!(msg[0], 5, "layout message version"); + // Class chunked, then (after the flags) 2 dimensions of 4 bytes. + assert_eq!(&msg[1..2], &[2]); + assert_eq!(&msg[3..5], &[2, 4]); + assert_eq!(&msg[5..13], &[1, 0, 0, 0x20, 8, 0, 0, 0]); + match DataLayout::parse(msg, 8, 8).unwrap() { + DataLayout::Chunked { + chunk_dimensions, + chunk_index_type, + .. + } => { + assert_eq!(chunk_dimensions, vec![HUGE_DIM as u32, 8]); + chunk_index_type.unwrap() + } + other => panic!("{other:?}"), + } + }; + let single = serialize_v4_single_chunk(&dims, 0x800, Some(20_000), Some(0), 8, 8, 5); + assert_eq!(parsed(&single), 1); + // Filtered size (8 bytes), filter mask, address. + assert_eq!(single.len(), 13 + 1 + 8 + 4 + 8); + assert_eq!( + parsed(&serialize_v4_fixed_array(&dims, 0x800, 8, 8, 10, 5)), + 3 + ); + let ea = ea_writer::serialize_v4_extensible_array(&dims, 0x800, 8, 8, 5); + assert_eq!(parsed(&ea), 4); + assert_eq!( + parsed(&serialize_v4_btree_v2(&dims, 0x800, 8, 8, 2048, 5)), + 5 + ); + + let slots = [Some(huge_chunk(0x1000, 20_000)), None]; + // FAHD: element size (address 8 + size 8 + mask 4) at byte 6. + let fa = build_fixed_array_at(&slots, 8, 8, true, 0x2000); + assert_eq!(fa[6], 20); + // FADB element 0: address, then the stored size in 8 bytes. + let fadb = 28; + let prefix = 4 + 1 + 1 + 8; + assert_eq!( + &fa[fadb + prefix..fadb + prefix + 8], + &0x1000u64.to_le_bytes() + ); + assert_eq!( + &fa[fadb + prefix + 8..fadb + prefix + 16], + &20_000u64.to_le_bytes() + ); + // AEHD: element size at byte 6 as well. + let ea = ea_writer::build_extensible_array_at(&slots, 8, 8, true, 0x2000); + assert_eq!(&ea[..4], b"EAHD"); + assert_eq!(ea[6], 20); + // BTHD: record size (element + 8-byte scaled offset) at bytes 10-11. + let chunk = huge_chunk(0x1000, 20_000); + let (bt, _) = + build_btree_v2_chunk_index_at(1, &[(vec![0], &chunk)], 8, 8, true, 0x2000).unwrap(); + assert_eq!(&bt[..4], b"BTHD"); + assert_eq!(u16::from_le_bytes([bt[10], bt[11]]), 28); + } + + #[test] + fn huge_chunks_refused_where_unsupported() { + // A chunk dimension of 2^32 or more. + assert!(matches!( + checked_chunk_dims(&[1 << 32], 1), + Err(FormatError::InvalidChunkDimensions(m)) if m.contains("2^32") + )); + // A chunk size that overflows u64. + assert!(matches!( + checked_chunk_dims(&[u32::MAX.into(), u32::MAX.into(), 2], 8), + Err(FormatError::Overflow(_)) + )); + assert_eq!( + checked_chunk_dims(&[HUGE_DIM], 8).unwrap(), + (vec![HUGE_DIM as u32], HUGE_BYTES) + ); + // Filters that cannot take a chunk that large. + for plugin in [ + PluginFilter::Lzf, + PluginFilter::Bzip2 { level: 9 }, + PluginFilter::Blosc { + codec: BloscCodec::Lz4, + level: 5, + shuffle: BloscShuffle::Byte, + }, + PluginFilter::Bitshuffle { + block_size: 0, + compression: BitshuffleCompression::Lz4, + }, + ] { + let options = ChunkOptions { + plugin: Some(plugin), + ..Default::default() + }; + assert!(matches!( + check_huge_chunk_filters(&options, HUGE_BYTES), + Err(FormatError::FilterError(m)) if m.contains("4 GiB") + )); + } + for options in [ + ChunkOptions { + deflate_level: Some(6), + shuffle: true, + fletcher32: true, + ..Default::default() + }, + ChunkOptions { + zstd_level: Some(3), + ..Default::default() + }, + ChunkOptions { + lz4: true, + ..Default::default() + }, + ] { + check_huge_chunk_filters(&options, HUGE_BYTES).unwrap(); + } + } + + /// Chunks are extracted row by row: the result is the element-by-element + /// split, edge padding included. + #[test] + fn extract_chunk_matches_elementwise_split() { + let shape = [5u64, 7, 3]; + let chunks = [2u64, 3, 2]; + let data: Vec = (0..5 * 7 * 3 * 2).map(|i| i as u8).collect(); + let n = chunk_count(&shape, &chunks); + assert_eq!(n, 3 * 3 * 2); + for i in 0..n { + let (offsets, chunk) = extract_chunk(&data, &shape, &chunks, 2, i); + assert_eq!(chunk.len(), 2 * 3 * 2 * 2); + for (e, pair) in chunk.as_chunks::<2>().0.iter().enumerate() { + let c = [e / 6, (e / 2) % 3, e % 2]; + let g: Vec = (0..3).map(|d| offsets[d] + c[d] as u64).collect(); + let expect = if (0..3).all(|d| g[d] < shape[d]) { + let flat = ((g[0] * 7 + g[1]) * 3 + g[2]) as usize; + [data[2 * flat], data[2 * flat + 1]] + } else { + [0, 0] + }; + assert_eq!(*pair, expect, "chunk {i} element {e}"); + } + } + // Data shorter than the shape: the missing elements stay zero. + let (_, chunk) = extract_chunk(&data[..5], &[4], &[4], 2, 0); + assert_eq!(chunk, [0, 1, 2, 3, 0, 0, 0, 0]); + } + // ---- Extensible Array tests ---- #[test] fn serialize_v4_extensible_array_roundtrip() { - let msg = ea_writer::serialize_v4_extensible_array(&[10], 0x4000, 8, 8); + let msg = ea_writer::serialize_v4_extensible_array(&[10], 0x4000, 8, 8, 4); let layout = DataLayout::parse(&msg, 8, 8).unwrap(); match layout { DataLayout::Chunked { @@ -2020,7 +2502,7 @@ mod tests { for info in &infos { // Skipped chunks are stored at the chunk's size (shuffled). assert_eq!( - info.chunk_size == (c * 8) as u32, + info.chunk_size == (c * 8) as u64, info.filter_mask != 0, "{info:?}" ); diff --git a/crates/clawhdf5-format/src/datatype.rs b/crates/clawhdf5-format/src/datatype.rs index 8e2120d..5e9a678 100644 --- a/crates/clawhdf5-format/src/datatype.rs +++ b/crates/clawhdf5-format/src/datatype.rs @@ -1216,6 +1216,26 @@ impl Datatype { } } + /// The highest datatype message version in this type's encoding, its + /// members' and base types' included (the version decides which HDF5 + /// releases can read it: 1-3 HDF5 1.8, 4 HDF5 1.12, 5 HDF5 2.0). + pub fn max_encoded_version(&self) -> u8 { + let own = self.serialize().first().map_or(0, |b| b >> 4); + let inner = match self { + Datatype::Compound { members, .. } => members + .iter() + .map(|m| m.datatype.max_encoded_version()) + .max() + .unwrap_or(0), + Datatype::Enumeration { base_type, .. } + | Datatype::VariableLength { base_type, .. } + | Datatype::Array { base_type, .. } + | Datatype::Complex { base_type, .. } => base_type.max_encoded_version(), + _ => 0, + }; + own.max(inner) + } + /// Check that this datatype can be written: every part of it has an /// on-disk encoding, and the encoding is one the reader (and libhdf5) /// accepts. [`Self::serialize`] cannot report errors, so the writer calls diff --git a/crates/clawhdf5-format/src/ea_writer.rs b/crates/clawhdf5-format/src/ea_writer.rs index de7859c..5951c50 100644 --- a/crates/clawhdf5-format/src/ea_writer.rs +++ b/crates/clawhdf5-format/src/ea_writer.rs @@ -18,9 +18,10 @@ pub(crate) fn serialize_v4_extensible_array( ea_address: u64, offset_size: u8, element_size: u32, + version: u8, ) -> Vec { let mut buf = Vec::new(); - buf.push(4); // version + buf.push(version); buf.push(2); // class = chunked buf.push(0x00); // flags @@ -87,7 +88,7 @@ pub fn build_extensible_array_at( ea_base_address: u64, ) -> Vec { let os = offset_size as usize; - let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots)); + let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots, length_size)); let elem_size = os + chunk_size_bytes.map_or(0, |n| n + 4); let client_id: u8 = if has_filters { 1 } else { 0 }; let arr_off_size = (MAX_NELMTS_BITS as usize).div_ceil(8); diff --git a/crates/clawhdf5-format/src/error.rs b/crates/clawhdf5-format/src/error.rs index 449ff13..9dd22b4 100644 --- a/crates/clawhdf5-format/src/error.rs +++ b/crates/clawhdf5-format/src/error.rs @@ -167,6 +167,19 @@ pub enum FormatError { VlDataError(String), /// Serialization error. SerializationError(String), + /// The file's library version bounds + /// ([`FileWriter::libver_bounds`](crate::file_writer::FileWriter::libver_bounds)) + /// do not allow what was asked for: `what` needs the format of HDF5 + /// `needs` or later, and the high bound is `high` (or the low bound is + /// above the high one, with `needs` the low bound). + LibverBound { + /// What cannot be written. + what: String, + /// The oldest release whose format holds it. + needs: crate::libver::LibVer, + /// The file's high bound. + high: crate::libver::LibVer, + }, /// Dataset is missing data. DatasetMissingData, /// Dataset is missing shape. @@ -450,6 +463,13 @@ impl fmt::Display for FormatError { FormatError::SerializationError(msg) => { write!(f, "serialization error: {msg}") } + FormatError::LibverBound { what, needs, high } => { + write!( + f, + "{what} needs the HDF5 {needs} file format, above the high \ + library version bound ({high})" + ) + } FormatError::DatasetMissingData => { write!(f, "dataset is missing data") } diff --git a/crates/clawhdf5-format/src/extensible_array.rs b/crates/clawhdf5-format/src/extensible_array.rs index 42cf595..889f2cf 100644 --- a/crates/clawhdf5-format/src/extensible_array.rs +++ b/crates/clawhdf5-format/src/extensible_array.rs @@ -224,7 +224,7 @@ fn read_element( }; Ok(( Some(ChunkInfo { - chunk_size: chunk_byte_size as u32, + chunk_size: chunk_byte_size, filter_mask: 0, offsets, address, @@ -259,7 +259,7 @@ fn read_element( }; Ok(( Some(ChunkInfo { - chunk_size: chunk_size as u32, + chunk_size, filter_mask, offsets, address, @@ -946,7 +946,7 @@ mod tests { assert_eq!(chunks.len(), 2); assert_eq!(chunks[0].address, base_addr); assert_eq!(chunks[0].offsets, vec![0]); - assert_eq!(chunks[0].chunk_size, chunk_byte_size as u32); + assert_eq!(chunks[0].chunk_size, chunk_byte_size); assert_eq!(chunks[1].address, base_addr + chunk_byte_size); assert_eq!(chunks[1].offsets, vec![20]); } diff --git a/crates/clawhdf5-format/src/file_writer.rs b/crates/clawhdf5-format/src/file_writer.rs index c64e1cf..4eb59ee 100644 --- a/crates/clawhdf5-format/src/file_writer.rs +++ b/crates/clawhdf5-format/src/file_writer.rs @@ -5,12 +5,13 @@ use crate::addr::saturating_usize; #[cfg(not(feature = "std"))] -use alloc::{format, vec, vec::Vec}; +use alloc::{format, string::String, vec, vec::Vec}; use crate::attribute::AttributeMessage; use crate::btree_v2_write::{BTreeV2Params, build_btree_v2}; use crate::chunked_write::{ - ChunkOptions, PrecompressedChunks, build_chunked_data_from_precompressed, precompress_chunks, + ChunkOptions, PrecompressedChunks, build_chunked_data_from_precompressed_libver, + precompress_chunks, }; use crate::data_layout::VdsMapping; use crate::dataspace::{Dataspace, DataspaceType}; @@ -31,6 +32,7 @@ pub use crate::type_builders::ProvenanceConfig; pub use crate::type_builders::{AttrValue, CompoundTypeBuilder, EnumTypeBuilder}; use crate::datatype::{CharacterSet, Datatype}; +use crate::libver::LibVer; pub(crate) const OFFSET_SIZE: u8 = 8; pub(crate) const LENGTH_SIZE: u8 = 8; @@ -168,13 +170,15 @@ pub(crate) fn build_dataset_oh( attrs: AttrStorage<'_>, fill_message: &[u8], refcount: u32, + layout_version: u8, ) -> Result, FormatError> { let mut w = ObjectHeaderWriter::new(); w.add_message_with_flags(MessageType::Datatype, dt.serialize(), 0x01); w.add_message(MessageType::Dataspace, ds.serialize(LENGTH_SIZE)); w.add_message_with_flags(MessageType::FillValue, fill_message.to_vec(), 0x01); + // Versions 3 and 4 encode a contiguous layout the same way. let mut dl = Vec::new(); - dl.push(4); // version + dl.push(layout_version); dl.push(1); // class = contiguous // An empty dataset has no storage: its address must be the undefined // address, as libhdf5 writes it. A real address with size 0 trips @@ -198,14 +202,16 @@ pub(crate) fn build_compact_dataset_oh( attrs: AttrStorage<'_>, fill_message: &[u8], refcount: u32, + layout_version: u8, ) -> Result, FormatError> { let mut w = ObjectHeaderWriter::new(); w.add_message_with_flags(MessageType::Datatype, dt.serialize(), 0x01); w.add_message(MessageType::Dataspace, ds.serialize(LENGTH_SIZE)); w.add_message_with_flags(MessageType::FillValue, fill_message.to_vec(), 0x01); - // Compact layout message: version=4, class=0, u16 size, inline data + // Compact layout message: version (3 and 4 are the same here), class=0, + // u16 size, inline data let mut dl = Vec::new(); - dl.push(4); // version + dl.push(layout_version); dl.push(0); // class = compact dl.extend_from_slice(&(data.len() as u16).to_le_bytes()); dl.extend_from_slice(data); @@ -1371,6 +1377,10 @@ pub struct FileWriter { /// file-space strategy (a File Space Info message in the superblock /// extension). page_size: Option, + /// Library version bounds: the low bound picks the format versions + /// written, the high bound limits the features allowed. + low: LibVer, + high: LibVer, } impl Default for FileWriter { @@ -1485,9 +1495,39 @@ impl FileWriter { alignment_threshold: 0, alignment_bytes: 0, page_size: None, + low: LibVer::V110, + high: LibVer::Latest, } } + /// Set the library version bounds, as libhdf5's `H5Pset_libver_bounds` + /// (h5py's `libver=(low, high)`): the oldest HDF5 release whose format + /// the file uses (`low`), and the newest whose features it may use + /// (`high`). See [`crate::libver`] for what each bound changes. + /// + /// The default, `(LibVer::V110, LibVer::Latest)`, is what clawhdf5 has + /// always written: the HDF5 1.10 format (version-3 superblock, version-4 + /// layouts with the 1.10 chunk indexes), readable by HDF5 1.10 and later. + /// + /// `(LibVer::V18, LibVer::V18)` writes a file HDF5 1.8 can read — the + /// low bound libhdf5 2.0 uses by default: a version-2 superblock, + /// version-3 layouts, and a version-1 B-tree for every chunked dataset, + /// resizable ones included; [`Self::finish`] then fails with + /// [`FormatError::LibverBound`] for anything HDF5 1.8 cannot read + /// (virtual datasets, a paged file, the 1.12 reference types, native + /// complex numbers). With a low bound of 1.8 and a later high bound + /// such objects are written in the newer format, as libhdf5 writes them; + /// the rest of the file stays readable by 1.8. A chunk of 4 GiB or more + /// always takes HDF5 2.0's version-5 layout (never a version-1 B-tree) + /// and so needs a high bound of at least [`LibVer::V200`]. + /// + /// A low bound above the high bound makes [`Self::finish`] fail. + pub fn libver_bounds(&mut self, low: LibVer, high: LibVer) -> &mut Self { + self.low = low; + self.high = high; + self + } + /// Set global file alignment: datasets with raw data >= `threshold` bytes /// will have their data aligned to `bytes` boundary. /// @@ -1583,6 +1623,33 @@ impl FileWriter { ))); } + let (low, high) = (self.low, self.high); + let within_bounds = |what: &dyn Fn() -> String, needs: LibVer| { + if needs > high { + Err(FormatError::LibverBound { + what: what(), + needs, + high, + }) + } else { + Ok(()) + } + }; + within_bounds(&|| format!("a low library version bound of {low}"), low)?; + if page_size.is_some() { + within_bounds(&|| "the paged file-space strategy".into(), LibVer::V110)?; + } + // Versions 3 of the layout message and 2 of the superblock are what + // HDF5 1.8 reads; 1.10 added version 4 (with its chunk indexes) and + // version 3. A paged file needs the version-3 superblock whatever + // the low bound (libhdf5 raises it as far as the high bound allows). + let layout_version: u8 = if low < LibVer::V110 { 3 } else { 4 }; + let superblock_version: u8 = if low < LibVer::V110 && page_size.is_none() { + 2 + } else { + 3 + }; + // The group tree, in layout order: groups depth-first from the root, // then every group's datasets in the same order. let tree = writer_tree::build(self.root, self.track_order)?; @@ -1621,9 +1688,20 @@ impl FileWriter { let ds_attrs = all_ds.iter().flat_map(|d| &d.attrs); for a in group_attrs.chain(ds_attrs) { a.datatype.check_encodable()?; + within_bounds( + &|| format!("the datatype of attribute {:?}", a.name), + LibVer::for_datatype_version(a.datatype.max_encoded_version()), + )?; } for d in &all_ds { d.dt.check_encodable()?; + within_bounds( + &|| "a dataset's datatype".into(), + LibVer::for_datatype_version(d.dt.max_encoded_version()), + )?; + if d.virtual_sources.is_some() { + within_bounds(&|| "a virtual dataset".into(), LibVer::V110)?; + } } let is_vds: Vec = all_ds.iter().map(|d| d.virtual_sources.is_some()).collect(); @@ -1749,10 +1827,12 @@ impl FileWriter { elem_size, &d.chunk_options, )?; - let result = build_chunked_data_from_precompressed( + let result = build_chunked_data_from_precompressed_libver( &pre, dummy_cursor, d.maxshape.as_deref(), + low, + high, )?; dummy_cursor += result.data_bytes.len() as u64; let oh = build_chunked_dataset_oh( @@ -1785,6 +1865,7 @@ impl FileWriter { }, &d.fill_message, d.refcount, + layout_version, )?; dummy_blobs.push(DataBlob { data: vec![], @@ -1804,6 +1885,7 @@ impl FileWriter { }, &d.fill_message, d.refcount, + layout_version, )?; dummy_blobs.push(DataBlob { data: vec![], @@ -1904,13 +1986,15 @@ impl FileWriter { let base_address = cursor2 as u64; // Reuse precompressed chunks from Pass 1 — avoids re-compressing // the same data a second time. - let result = build_chunked_data_from_precompressed( + let result = build_chunked_data_from_precompressed_libver( dummy_blobs[i] .precompressed .as_ref() .expect("chunked dataset missing precompressed cache"), base_address, d.maxshape.as_deref(), + low, + high, )?; cursor2 += result.data_bytes.len(); let oh = build_chunked_dataset_oh( @@ -1944,6 +2028,7 @@ impl FileWriter { }, &d.fill_message, d.refcount, + layout_version, )?; ds_blobs2.push(DataBlob { data: vec![], @@ -1973,6 +2058,7 @@ impl FileWriter { }, &d.fill_message, d.refcount, + layout_version, )?; let mut data = vec![0u8; padding]; data.extend_from_slice(&d.raw); @@ -1997,7 +2083,7 @@ impl FileWriter { let mut buf = Vec::with_capacity(cursor2); let sb = Superblock { - version: 3, + version: superblock_version, offset_size: OFFSET_SIZE, length_size: LENGTH_SIZE, base_address: 0, @@ -2783,4 +2869,150 @@ mod tests { assert_eq!(sb.version, 3); assert_eq!(sb.page_size, None); } + + fn layout_of(bytes: &[u8], name: &str) -> Vec { + let sb = Superblock::parse(bytes, 0).unwrap(); + let addr = resolve_path_any(bytes, &sb, name).unwrap(); + let hdr = ObjectHeader::parse(bytes, addr as usize, 8, 8).unwrap(); + hdr.messages + .iter() + .find(|m| m.msg_type == MessageType::DataLayout) + .unwrap() + .data + .clone() + } + + #[test] + fn libver_v18_writes_the_1_8_format() { + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::V18); + fw.create_dataset("contig").with_f64_data(&[1.0, 2.0]); + fw.create_dataset("compact").with_f64_data(&[3.0]).compact(); + fw.create_dataset("grow") + .with_f64_data(&[1.0, 2.0, 3.0]) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[2]); + fw.create_dataset("none") + .with_f64_data(&[]) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[2]); + let bytes = fw.finish().unwrap(); + assert_eq!(Superblock::parse(&bytes, 0).unwrap().version, 2); + assert_eq!(layout_of(&bytes, "contig")[..2], [3, 1]); + assert_eq!(layout_of(&bytes, "compact")[..2], [3, 0]); + let grow = layout_of(&bytes, "grow"); + // Version 3, chunked, 2 dimensions (the element size is the last), + // B-tree address, chunk dims 2 and 8. + assert_eq!(grow[..3], [3, 2, 2]); + assert_eq!(grow[11..], [2, 0, 0, 0, 8, 0, 0, 0]); + let root = u64::from_le_bytes(grow[3..11].try_into().unwrap()) as usize; + assert_eq!(&bytes[root..root + 5], b"TREE\x01"); + // No chunks, no tree. + assert_eq!(layout_of(&bytes, "none")[3..11], [0xff; 8]); + assert_eq!(read_dataset_f64(&bytes, "grow"), vec![1.0, 2.0, 3.0]); + assert_eq!(read_dataset_f64(&bytes, "contig"), vec![1.0, 2.0]); + assert_eq!(read_dataset_f64(&bytes, "compact"), vec![3.0]); + } + + #[test] + fn default_libver_bounds_keep_the_1_10_format() { + let mut fw = FileWriter::new(); + fw.create_dataset("contig").with_f64_data(&[1.0, 2.0]); + fw.create_dataset("grow") + .with_f64_data(&[1.0, 2.0, 3.0]) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[2]); + let default = fw.finish().unwrap(); + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V110, LibVer::Latest); + fw.create_dataset("contig").with_f64_data(&[1.0, 2.0]); + fw.create_dataset("grow") + .with_f64_data(&[1.0, 2.0, 3.0]) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[2]); + assert_eq!(fw.finish().unwrap(), default); + assert_eq!(layout_of(&default, "contig")[0], 4); + assert_eq!(layout_of(&default, "grow")[..2], [4, 2]); + } + + #[test] + fn libver_high_bound_refuses_newer_features() { + let bound = |r: Result, FormatError>, needs: LibVer| match r { + Err(FormatError::LibverBound { needs: n, high, .. }) => { + assert_eq!((n, high), (needs, LibVer::V18)); + } + other => panic!("expected a bound error, got {other:?}"), + }; + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::V18); + fw.create_dataset("z") + .with_native_complex_f64_data(&[[1.0, 2.0]]); + bound(fw.finish(), LibVer::V200); + + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::V18); + fw.create_dataset("x").with_f64_data(&[1.0]).set_attr( + "z", + AttrValue::Raw { + datatype: crate::type_builders::make_native_complex_f64_type(), + shape: vec![], + data: vec![0; 16], + }, + ); + bound(fw.finish(), LibVer::V200); + + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::V18); + fw.create_dataset("r").with_compound_data( + Datatype::Reference { + size: 16, + ref_type: crate::datatype::ReferenceType::Object2, + }, + vec![0; 16], + 1, + ); + bound(fw.finish(), LibVer::V112); + + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::V18); + fw.create_dataset("src").with_f64_data(&[1.0, 2.0]); + fw.create_dataset("vds") + .with_shape(&[2]) + .with_f64_data(&[]) + .with_virtual_sources(vec![VdsMapping { + source_file: ".".into(), + source_dataset: "src".into(), + source_selection: sel_all(), + virtual_selection: sel_hyper_1d(0, 2), + }]); + bound(fw.finish(), LibVer::V110); + + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::V18) + .with_page_size(4096); + bound(fw.finish(), LibVer::V110); + + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V110, LibVer::V18); + bound(fw.finish(), LibVer::V110); + } + + #[test] + fn libver_low_v18_high_latest_allows_newer_objects() { + // As libhdf5 does: the object that needs a newer format gets it, + // the rest of the file keeps the 1.8 format. + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::Latest); + fw.create_dataset("z") + .with_native_complex_f64_data(&[[1.0, 2.0]]); + let bytes = fw.finish().unwrap(); + assert_eq!(Superblock::parse(&bytes, 0).unwrap().version, 2); + let mut fw = FileWriter::new(); + fw.libver_bounds(LibVer::V18, LibVer::Latest) + .with_page_size(4096); + fw.create_dataset("x").with_f64_data(&[1.0]); + let bytes = fw.finish().unwrap(); + assert_eq!(Superblock::parse(&bytes, 0).unwrap().version, 3); + assert_eq!(layout_of(&bytes, "x")[0], 3); + } } diff --git a/crates/clawhdf5-format/src/filters.rs b/crates/clawhdf5-format/src/filters.rs index 28070d4..d902041 100644 --- a/crates/clawhdf5-format/src/filters.rs +++ b/crates/clawhdf5-format/src/filters.rs @@ -340,10 +340,23 @@ pub fn decompress_chunk_exact_with<'s>( } else { MAX_DECOMPRESS_SIZE }; - let size_hint = if ctx.max_output != 0 { + // The bound is the output's size when only size-preserving + // filters (shuffle, Fletcher32) remain to be undone; before + // any other filter (a second deflate, N-Bit, ...) it is only + // a ceiling, and the output starts smaller and grows. + // Reserving the bound of a 4 GiB chunk for a stage a few MiB + // long doubled the peak memory of reading it (8.5 GiB, now + // 4.0, for `huge_chunks_filtered.h5`'s double deflate). + let exact = pipeline.filters[..i].iter().enumerate().all(|(j, f)| { + filter_skipped(filter_mask, j) + || matches!(f.filter_id, FILTER_SHUFFLE | FILTER_FLETCHER32) + }); + let size_hint = if ctx.max_output == 0 { + input.len().saturating_mul(4).min(1 << 20) + } else if exact { ctx.max_output } else { - input.len().saturating_mul(4).min(1 << 20) + input.len().saturating_mul(4).min(ctx.max_output) }; let inflater = scratch .inflater @@ -429,7 +442,10 @@ pub fn compress_chunk_masked( "more than 32 filters in a pipeline".into(), )); } - let mut result = data.to_vec(); + // The input is not copied: the first filter reads it where it is (a + // chunk may be 4 GiB or more), and each filter's output replaces the + // previous one. + let mut owned: Option> = None; let mut mask = 0u32; for (i, filter) in pipeline.filters.iter().enumerate() { let ctx = FilterContext { @@ -437,7 +453,8 @@ pub fn compress_chunk_masked( element_size: element_size as usize, max_output: 0, }; - let out = match filter_registry::encode(&result, &ctx) { + let result: &[u8] = owned.as_deref().unwrap_or(data); + let out = match filter_registry::encode(result, &ctx) { Ok(out) if FAIL_UNLESS_SMALLER.contains(&filter.filter_id) && out.len() >= result.len() => { @@ -449,13 +466,13 @@ pub fn compress_chunk_masked( r => r, }; match out { - Ok(out) => result = out, + Ok(out) => owned = Some(out), Err(e @ FormatError::UnsupportedFilter(_)) => return Err(e), Err(_) if filter.flags & FILTER_FLAG_OPTIONAL != 0 => mask |= 1 << i, Err(e) => return Err(e), } } - Ok((result, mask)) + Ok((owned.unwrap_or_else(|| data.to_vec()), mask)) } /// The filters compiled into this build, sorted by ID (see @@ -1339,6 +1356,9 @@ fn deflate_compress(data: &[u8], level: u32) -> Result, FormatError> { deflate_bounded(data, level).map_err(FormatError::CompressionError) } +/// Largest compression output reserved at its worst-case size up front. +const DEFLATE_EXACT_BOUND: usize = 64 << 20; + /// Deflate `data` into a zlib stream in one pass, into a buffer sized for the /// worst case up front (the same reasoning as [`inflate_bounded`]). #[cfg(feature = "deflate")] @@ -1347,24 +1367,44 @@ pub(crate) fn deflate_bounded(data: &[u8], level: u32) -> Result, String // zlib's compressBound, plus the zlib header and trailer. let bound = data.len() + (data.len() >> 12) + (data.len() >> 14) + (data.len() >> 25) + 13 + 6; + // flate2's Rust backends (zlib-rs, miniz_oxide) zero the whole spare + // capacity on each call, so a large input's worst-case bound would be + // memory held for nothing (4 GiB for a 4 GiB chunk that deflates to a + // few MiB): past 64 MiB the output starts at 1/16 of the bound and + // doubles as needed. + let first = if bound <= DEFLATE_EXACT_BOUND { + bound + } else { + bound / 16 + }; let mut out = Vec::new(); - out.try_reserve_exact(bound) + out.try_reserve_exact(first) .map_err(|e| format!("deflate: cannot allocate output: {e}"))?; let mut deflater = Compress::new(Compression::new(level), true); loop { let (in_before, out_before) = (deflater.total_in(), deflater.total_out()); + let rest = &data[saturating_usize(in_before)..]; + // zlib takes at most u32::MAX input bytes per call, and `Finish` + // ends the stream after the bytes it took: a chunk of 4 GiB or more + // was cut at 4 GiB - 1. Finish only once the rest fits one call. + let flush = if rest.len() > u32::MAX as usize { + FlushCompress::None + } else { + FlushCompress::Finish + }; let status = deflater - .compress_vec( - &data[saturating_usize(in_before)..], - &mut out, - FlushCompress::Finish, - ) + .compress_vec(rest, &mut out, flush) .map_err(|e| format!("deflate: {e}"))?; match status { - Status::StreamEnd => return Ok(out), - // The bound should make running out of room unreachable; grow - // rather than fail if it happens. + Status::StreamEnd => { + if out.capacity() - out.len() > DEFLATE_EXACT_BOUND { + out.shrink_to_fit(); + } + return Ok(out); + } + // Out of room (the bound makes it unreachable below + // `DEFLATE_EXACT_BOUND`): grow rather than fail. Status::Ok | Status::BufError if out.len() == out.capacity() => out .try_reserve(out.capacity().max(4096)) .map_err(|e| format!("deflate: cannot allocate output: {e}"))?, @@ -1395,10 +1435,12 @@ const LZ4_DEFAULT_BLOCK_SIZE: usize = 1 << 30; /// * The legacy clawhdf5 framing (up to 2.7.0): a 4-byte little-endian size /// followed by one raw LZ4 block. libhdf5 cannot read it. /// -/// They are told apart unambiguously: an HDF5 chunk is smaller than 4 GiB, so -/// the registered format's big-endian `u64` size always starts with four zero -/// bytes and the whole chunk is at least 12 bytes; a legacy chunk starts with -/// four zero bytes only when it is empty, and is then 5 bytes long. +/// They are told apart unambiguously: for a chunk under 4 GiB the registered +/// format's big-endian `u64` size starts with four zero bytes and the whole +/// chunk is at least 12 bytes; a legacy chunk starts with four zero bytes +/// only when it is empty, and is then 5 bytes long. A chunk of 4 GiB or more +/// (HDF5 2.0) is always the registered format: clawhdf5 never wrote legacy +/// chunks that large. /// /// Every size read from the payload is bounded against `expected_bytes` (the /// pipeline's declared chunk size) before it sizes an allocation, so a crafted @@ -1416,14 +1458,18 @@ fn lz4_decompress(data: &[u8], expected_bytes: usize) -> Result, FormatE "lz4: declared size exceeds chunk size".into(), )); } - if size > MAX_DECOMPRESS_SIZE { + // Without a chunk size, a ceiling; with one, the chunk size is the + // bound (a chunk may be 4 GiB or more, like a deflated one). + if expected_bytes == 0 && size > MAX_DECOMPRESS_SIZE { return Err(FormatError::DecompressionError( "lz4: declared size exceeds limit".into(), )); } Ok(()) }; - if data.len() >= 12 && data[..4] == [0, 0, 0, 0] { + // A chunk of 4 GiB or more cannot be legacy (its size field is 32 + // bits), and its registered-format size does not start with zeros. + if data.len() >= 12 && (data[..4] == [0, 0, 0, 0] || expected_bytes > u32::MAX as usize) { return lz4_decompress_hdf5(data, check_size); } // Legacy clawhdf5 framing: 4-byte LE size + one LZ4 block. @@ -1441,9 +1487,10 @@ fn lz4_decompress_hdf5( ) -> Result, FormatError> { let err = |m: &str| FormatError::DecompressionError(format!("lz4: {m}")); let be32 = |b: &[u8]| u32::from_be_bytes([b[0], b[1], b[2], b[3]]) as usize; - // The first four bytes are zero (checked by the caller), so the size is - // the low 32 bits of the big-endian u64. - let orig_size = be32(&data[4..8]); + let orig_size = u64::from_be_bytes([ + data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], + ]); + let orig_size = usize::try_from(orig_size).map_err(|_| err("chunk too large"))?; check_size(orig_size)?; let block_size = be32(&data[8..12]).min(orig_size); if block_size == 0 && orig_size != 0 { @@ -2717,6 +2764,39 @@ mod tests { assert_eq!(decompressed, data); } + /// A chunk's size bounds an LZ4 chunk, not the 256 MiB ceiling for an + /// unknown size: a 300 MiB chunk was refused ("declared size exceeds + /// limit"). + #[test] + #[cfg(feature = "lz4")] + fn lz4_chunks_over_256_mib_decode() { + let mut data = vec![0u8; 300 << 20]; + data[12345] = 7; + let compressed = lz4_compress(&data, &[]).unwrap(); + let decompressed = lz4_decompress(&compressed, data.len()).unwrap(); + assert!(decompressed == data); + // Without a chunk size the ceiling still applies. + assert!(lz4_decompress(&compressed, 0).is_err()); + } + + /// A chunk of 4 GiB or more is always in the registered framing, whose + /// big-endian size then does not start with four zero bytes: its size + /// is read whole (here larger than the chunk, so refused before any + /// allocation), not taken for a legacy 4-byte size. + #[test] + #[cfg(all(feature = "lz4", target_pointer_width = "64"))] + fn lz4_chunks_of_4_gib_use_the_registered_framing() { + let chunk = (1usize << 32) + 8; + let mut data = ((chunk + 8) as u64).to_be_bytes().to_vec(); + data.extend_from_slice(&(1u32 << 30).to_be_bytes()); + data.extend_from_slice(&[0; 8]); + let err = lz4_decompress(&data, chunk).unwrap_err(); + assert!( + matches!(&err, FormatError::DecompressionError(m) if m.contains("exceeds chunk size")), + "{err:?}" + ); + } + #[test] #[cfg(feature = "lz4")] fn pipeline_lz4_only() { diff --git a/crates/clawhdf5-format/src/fixed_array.rs b/crates/clawhdf5-format/src/fixed_array.rs index de80b7a..dbe00ae 100644 --- a/crates/clawhdf5-format/src/fixed_array.rs +++ b/crates/clawhdf5-format/src/fixed_array.rs @@ -383,7 +383,7 @@ fn parse_fa_element( offset_size: u8, element_size: u8, chunk_byte_size: u64, -) -> Result, FormatError> { +) -> Result, FormatError> { let os = offset_size as usize; if client_id == 0 { // Non-filtered: element is just the chunk address. @@ -393,7 +393,7 @@ fn parse_fa_element( return Ok(None); } let address = read_offset(file_data, abs, offset_size)?; - Ok(Some((address, chunk_byte_size as u32, 0))) + Ok(Some((address, chunk_byte_size, 0))) } else { // Filtered: address(offset_size) + chunk_size(variable) + filter_mask(4) let es = element_size as usize; @@ -418,7 +418,7 @@ fn parse_fa_element( file_data[fm_off + 2], file_data[fm_off + 3], ]); - Ok(Some((address, chunk_size as u32, filter_mask))) + Ok(Some((address, chunk_size, filter_mask))) } } @@ -688,7 +688,7 @@ mod tests { assert_eq!(c.address, base_addr + i as u64 * chunk_byte_size as u64); assert_eq!(c.offsets, vec![i as u64 * 20]); assert_eq!(c.filter_mask, 0); - assert_eq!(c.chunk_size, chunk_byte_size as u32); + assert_eq!(c.chunk_size, chunk_byte_size as u64); } } diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index 9ce3ff9..32c3834 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -61,6 +61,7 @@ pub mod addr; pub mod attribute; pub mod attribute_info; pub mod btree_v1; +mod btree_v1_write; pub mod btree_v2; mod btree_v2_write; mod bulk_alloc; @@ -107,6 +108,7 @@ pub mod group_v1; pub mod group_v2; #[cfg(feature = "parallel")] pub mod lane_partition; +pub mod libver; pub mod link_info; pub mod link_message; pub mod local_heap; diff --git a/crates/clawhdf5-format/src/libver.rs b/crates/clawhdf5-format/src/libver.rs new file mode 100644 index 0000000..42fe4fa --- /dev/null +++ b/crates/clawhdf5-format/src/libver.rs @@ -0,0 +1,93 @@ +//! Library version bounds for writing: which HDF5 releases can read a file. +//! +//! libhdf5 picks the version of every object it writes from the file's +//! *low* bound (`H5Pset_libver_bounds`; h5py's `libver=`): the oldest +//! format version that holds the object, but never older than the one the +//! low bound names. The *high* bound caps it: a feature that needs a newer +//! format than the high bound is an error. [`LibVer`] names the same +//! releases, and [`crate::file_writer::FileWriter::libver_bounds`] sets them. +//! +//! What the low bound changes in what clawhdf5 writes: +//! +//! | | low [`LibVer::V18`] | low [`LibVer::V110`] or later (the default) | +//! |---|---|---| +//! | superblock | version 2 | version 3 | +//! | data layout message | version 3 | version 4 | +//! | chunk index | version-1 B-tree (every chunked dataset) | single chunk, Fixed Array, Extensible Array or version-2 B-tree, as libhdf5 picks | +//! +//! Everything else (version-2 object headers, link and group-info messages, +//! dense storage in fractal heaps with version-2 B-trees, filter pipeline +//! version 2, fill value version 3, datatype versions up to 3) is the same +//! and already readable by HDF5 1.8. +//! +//! What the high bound refuses: anything that needs 1.10 (virtual datasets, +//! the paged file-space strategy) above [`LibVer::V18`], the 1.12 reference +//! types (datatype version 4) above [`LibVer::V110`], and HDF5 2.0's native +//! complex numbers (datatype version 5) above [`LibVer::V114`]. +//! `libver_bounds(LibVer::V18, LibVer::V18)` therefore writes a file HDF5 +//! 1.8 can read, or fails. + +use core::fmt; + +/// An HDF5 library release, as a bound on the file format versions a writer +/// may use (libhdf5's `H5F_libver_t`). Ordered oldest first. +/// +/// There is no `Earliest`: clawhdf5 cannot write the pre-1.8 format +/// (symbol-table groups, version-1 object headers). +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +#[non_exhaustive] +pub enum LibVer { + /// HDF5 1.8 (`H5F_LIBVER_V18`, h5py `'v108'`). + V18, + /// HDF5 1.10 (`H5F_LIBVER_V110`, h5py `'v110'`). + V110, + /// HDF5 1.12 (`H5F_LIBVER_V112`, h5py `'v112'`). + V112, + /// HDF5 1.14 (`H5F_LIBVER_V114`, h5py `'v114'`). + V114, + /// HDF5 2.0 (`H5F_LIBVER_V200`). + V200, + /// The newest format this build of clawhdf5 writes + /// (`H5F_LIBVER_LATEST`, h5py `'latest'`). + Latest, +} + +impl LibVer { + /// The release a datatype message of this version first appeared in: + /// versions 1-3 are readable by HDF5 1.8, 4 needs 1.12, 5 needs 2.0. + pub(crate) fn for_datatype_version(version: u8) -> Self { + match version { + 0..=3 => LibVer::V18, + 4 => LibVer::V112, + _ => LibVer::V200, + } + } +} + +impl fmt::Display for LibVer { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + f.write_str(match self { + LibVer::V18 => "1.8", + LibVer::V110 => "1.10", + LibVer::V112 => "1.12", + LibVer::V114 => "1.14", + LibVer::V200 => "2.0", + LibVer::Latest => "latest", + }) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn ordered_oldest_first() { + assert!(LibVer::V18 < LibVer::V110); + assert!(LibVer::V114 < LibVer::V200); + assert!(LibVer::V200 < LibVer::Latest); + assert_eq!(LibVer::for_datatype_version(3), LibVer::V18); + assert_eq!(LibVer::for_datatype_version(4), LibVer::V112); + assert_eq!(LibVer::for_datatype_version(5), LibVer::V200); + } +} diff --git a/crates/clawhdf5-format/src/parallel_read.rs b/crates/clawhdf5-format/src/parallel_read.rs index 606173c..38288b4 100644 --- a/crates/clawhdf5-format/src/parallel_read.rs +++ b/crates/clawhdf5-format/src/parallel_read.rs @@ -255,7 +255,7 @@ pub fn decompress_chunks_lane_partitioned_in( for &local in &indices { let index = batch.start + local; let chunk_info = &chunks[index]; - let size = chunk_info.chunk_size as usize; + let size = crate::addr::saturating_usize(chunk_info.chunk_size); let raw_chunk = raw_bytes.get(index, &reqs[index])?; let decompressed = decompress_chunk_exact( @@ -426,7 +426,7 @@ mod tests { for i in 0..8u64 { let len = if short && i == 5 { 16 } else { 32 }; infos.push(ChunkInfo { - chunk_size: len as u32, + chunk_size: len as u64, filter_mask: 0, offsets: vec![i * 8], address: file.len() as u64, diff --git a/crates/clawhdf5-format/src/partial_read.rs b/crates/clawhdf5-format/src/partial_read.rs index 3361485..f9896de 100644 --- a/crates/clawhdf5-format/src/partial_read.rs +++ b/crates/clawhdf5-format/src/partial_read.rs @@ -243,6 +243,58 @@ fn copy_overlap( } } +/// Copy the part of unfiltered chunk `chunk` (`chunk_bytes` long, shape +/// `chunk_shape`) that overlaps the box into `out`, reading only the runs of +/// the overlap from the file. The whole chunk must still lie inside the +/// file, as it must when it is fetched whole. +#[allow(clippy::too_many_arguments)] +fn read_unfiltered_overlap( + file_data: &S, + chunk: &crate::chunked_read::ChunkInfo, + chunk_bytes: usize, + chunk_shape: &[u64], + elem_size: usize, + out: &mut [u8], + box_start: &[u64], + box_extent: &[u64], +) -> Result<(), FormatError> { + let rank = chunk_shape.len(); + let origin = &chunk.offsets[..rank]; + let (lo, extent): (Vec, Vec) = (0..rank) + .map(|d| { + let lo = origin[d].max(box_start[d]); + let hi = origin[d] + .saturating_add(chunk_shape[d]) + .min(box_start[d] + box_extent[d]); + (lo, hi.saturating_sub(lo)) + }) + .unzip(); + let file_len = crate::storage::len_usize(file_data); + let base = crate::addr::to_usize(chunk.address)?; + if base > file_len || chunk_bytes > file_len - base { + return Err(FormatError::UnexpectedEof { + expected: base.saturating_add(chunk_bytes), + available: file_len, + }); + } + let overlap = Selection::Hyperslab { + start: lo.iter().zip(origin).map(|(l, o)| l - o).collect(), + stride: vec![1; rank], + count: extent.clone(), + block: vec![1; rank], + }; + let rows = crate::gather::gather_storage( + file_data, + chunk.address, + chunk_bytes, + chunk_shape, + elem_size, + &overlap, + )?; + copy_overlap(&rows, &lo, &extent, out, box_start, box_extent, elem_size); + Ok(()) +} + /// Read `selection` without materialising the whole dataset, when that is /// possible and worthwhile. `Ok(None)` means "use the full-read path": an /// `All`/`None`/invalid selection, a layout this doesn't handle (compact, @@ -281,9 +333,45 @@ pub fn read_selection_in( offset_size: u8, length_size: u8, selection: &Selection, +) -> Result>, FormatError> { + read_selection_filled_in( + file_data, + layout, + dataspace, + elem_size, + pipeline, + offset_size, + length_size, + selection, + None, + ) +} + +/// [`read_selection_in`] for a dataset whose fill value is `fill` (one +/// element's bytes; `None` or all zeros is the default fill): the elements +/// of a chunked dataset's selection that lie in chunks never written read as +/// `fill`, as they do in a full read. Only the chunks the selection's +/// bounding box overlaps are read, so a selection of a few elements of a +/// dataset whose chunks are 4 GiB or more costs one decoded chunk (a +/// filtered chunk has to be decoded whole) or, unfiltered, only the bytes +/// it selects. +#[allow(clippy::too_many_arguments)] +pub fn read_selection_filled_in( + file_data: &S, + layout: &DataLayout, + dataspace: &Dataspace, + elem_size: usize, + pipeline: Option<&FilterPipeline>, + offset_size: u8, + length_size: u8, + selection: &Selection, + fill: Option<&[u8]>, ) -> Result>, FormatError> { let dims = &dataspace.dimensions; - if dims.is_empty() || elem_size == 0 { + // A fill value that is not one element's bytes is the full path's to + // interpret. + let odd_fill = fill.is_some_and(|f| !f.is_empty() && f.len() != elem_size); + if dims.is_empty() || elem_size == 0 || odd_fill { return Ok(None); } let total = dataspace.checked_num_elements()?; @@ -341,8 +429,18 @@ pub fn read_selection_in( return Ok(None); } let mut boxed = alloc_output(checked_byte_len(box_elements, elem_size)?)?; + if let Some(fill) = fill.filter(|f| <[u8]>::len(f) == elem_size && f.iter().any(|&b| b != 0)) { + for element in boxed.chunks_exact_mut(elem_size) { + element.copy_from_slice(fill); + } + } match layout { + // No chunk was ever written: every element is the fill value. + DataLayout::Chunked { + btree_address: None, + .. + } if fill.is_some() => {} DataLayout::Chunked { btree_address: Some(_), .. @@ -373,6 +471,25 @@ pub fn read_selection_in( }) }) .collect(); + // An unfiltered chunk of a file that is not in memory: fetch + // only the rows the box needs, not the whole chunk (which may be + // 4 GiB or more). + let (direct, wanted): (Vec<_>, Vec<_>) = wanted.into_iter().partition(|c| { + file_data.as_contiguous().is_none() + && pipeline.is_none_or(|pl| all_filters_skipped(pl, c.filter_mask)) + }); + for chunk in direct { + read_unfiltered_overlap( + file_data, + chunk, + chunk_bytes, + &chunk_shape, + elem_size, + &mut boxed, + &box_start, + &box_extent, + )?; + } // Their stored bytes, batch by batch when the file is not in // memory; each batch's chunks are decoded into this thread's // reusable buffers before the next batch is fetched. diff --git a/crates/clawhdf5-format/tests/raw_fetch_bounds.rs b/crates/clawhdf5-format/tests/raw_fetch_bounds.rs index 8217323..c8b2471 100644 --- a/crates/clawhdf5-format/tests/raw_fetch_bounds.rs +++ b/crates/clawhdf5-format/tests/raw_fetch_bounds.rs @@ -151,7 +151,7 @@ fn crafted() -> (Vec, Chunked, Vec) { // v1 B-tree key (size, filter mask, offsets + 0) then the child // address. let mut pat = Vec::new(); - pat.extend_from_slice(&c.chunk_size.to_le_bytes()); + pat.extend_from_slice(&(c.chunk_size as u32).to_le_bytes()); pat.extend_from_slice(&c.filter_mask.to_le_bytes()); // The key holds one offset per dimension plus the element offset // (0); `offsets` may or may not list that last one. @@ -173,7 +173,7 @@ fn crafted() -> (Vec, Chunked, Vec) { assert!( chunks .iter() - .all(|c| c.chunk_size == HUGE && c.address == blob) + .all(|c| c.chunk_size == u64::from(HUGE) && c.address == blob) ); (bytes, ds, chunks) } @@ -313,7 +313,7 @@ fn large_reads_are_fetched_in_batches() { let data: Vec = (0..2 * CHUNK).map(|i| (i % 251) as u8).collect(); let chunks: Vec = (0..40u64) .map(|i| ChunkInfo { - chunk_size: CHUNK as u32, + chunk_size: CHUNK as u64, filter_mask: 0, offsets: vec![i * CHUNK as u64], address: (i % 2) * CHUNK as u64, diff --git a/crates/clawhdf5-netcdf4/README.md b/crates/clawhdf5-netcdf4/README.md index a5f63f0..426e823 100644 --- a/crates/clawhdf5-netcdf4/README.md +++ b/crates/clawhdf5-netcdf4/README.md @@ -36,17 +36,35 @@ let values: Vec = temp.read_f64()?; | Item | What | |---|---| -| `NetCDF4File` | `open`, `from_bytes`, `dimensions`, `variables`, `variable`, `global_attrs`, `group`, `group_names`, `nc_properties`, and `hdf5_file` for the underlying `clawhdf5::File` | -| `NetCDF4Group` | the same for a sub-group (`dimensions`, `variables`, `attrs`, nested `group`) | -| `Variable` | `name`, `shape`, `dimensions`, `nc_type`, `is_coordinate`, `attrs`, `cf_attributes`; `read_f64` (CF scale/offset and fill applied), `read_raw_f32`/`_f64`/`_i32`/`_i64`/`_u64`, `read_string`, `read_raw` | +| `NetCDF4File` | `open`, `from_bytes`, `dimensions`, `variables`, `variable_names`, `variable`, `global_attrs`, `group`, `group_names`, `nc_properties`, and `hdf5_file` for the underlying `clawhdf5::File` | +| `NetCDF4Group` | the same for a sub-group (`dimensions`, `variables`, `variable_names`, `attrs`, nested `group`) | +| `Variable` | `name`, `shape`, `stored_shape`, `dimensions`, `nc_type`, `is_coordinate`, `attrs`, `cf_attributes`; `read_f64` (CF scale/offset and fill applied), `read_raw_f32`/`_f64`/`_i32`/`_i64`/`_u64`, `read_string`, `read_raw` | | `Dimension` | `name`, `size`, `is_unlimited` (an unlimited dimension's `size` is its current length as netCDF-C reports it: the largest extent of the variables using it) | | `CfAttributes` | CF convention attributes: `units`, `long_name`, `standard_name`, `fill_value` (`_FillValue`), `missing_value`, `scale_factor`, `add_offset`, `valid_range`, `calendar`, `axis` | | `NcType` | the NetCDF type of a variable | -No cargo features. Tests compare against files written by netCDF4-python -(`tests/interop_tests.rs`; the CI job requires them with -`CLAWHDF5_REQUIRE_INTEROP=1`). What the HDF5 reader underneath cannot -read is listed in [`docs/known-issues.md`](../../docs/known-issues.md). +Variables and dimensions follow netCDF-C: + +- A variable's dimensions are the ones the file names: the ids in its + `_Netcdf4Coordinates` attribute, else the dimension scales its + `DIMENSION_LIST` references, found in its group or a parent group. Only + an axis the file names no dimension for (an HDF5 file not written by a + netCDF library) gets the first dimension of the group of the same size, + else an anonymous `dim_`. +- Dimension scales that are only dimensions are not variables; a dataset + `_nc4_non_coord_` is the variable ``. +- A variable along an unlimited dimension has the dimension's length: + `shape` is that length, and the reads return that many values, the + records the variable has not written as its fill value (`_FillValue`, + else netCDF's default for the type; NaN from `read_f64`). + `stored_shape` is the HDF5 dataset's extent. + +No cargo features. Tests compare against files written by netCDF4-python, +h5py dimension scales, h5netcdf and xarray, variable by variable with what +netCDF4-python reads (`tests/interop_tests.rs`; the CI job requires them +with `CLAWHDF5_REQUIRE_INTEROP=1`; the h5netcdf cases skip when h5netcdf is +not installed). What the HDF5 reader underneath cannot read is listed in +[`docs/known-issues.md`](../../docs/known-issues.md). ## License diff --git a/crates/clawhdf5-netcdf4/src/dimension.rs b/crates/clawhdf5-netcdf4/src/dimension.rs index 3e20163..4f88819 100644 --- a/crates/clawhdf5-netcdf4/src/dimension.rs +++ b/crates/clawhdf5-netcdf4/src/dimension.rs @@ -22,73 +22,123 @@ pub struct Dimension { pub is_unlimited: bool, } -/// Extract dimensions from an HDF5 group (root or subgroup). +/// A dimension scale of one group: the dataset that defines a dimension. +#[derive(Debug, Clone)] +pub(crate) struct Scale { + /// Object header address of the scale's dataset (what a variable's + /// `DIMENSION_LIST` references). + pub address: u64, + /// Its `_Netcdf4Dimid` (what a variable's `_Netcdf4Coordinates` lists). + pub dimid: Option, + /// Index of its dimension in [`GroupDims::dims`]. + pub dim: usize, +} + +/// The dimensions a group defines, with the scales that define them. +#[derive(Debug, Clone, Default)] +pub(crate) struct GroupDims { + /// The group's dimensions, in `_Netcdf4Dimid` order (then discovery order). + pub dims: Vec, + /// The dimension scales behind `dims`; empty when the group has no + /// dimension scale and `dims` were inferred from 1-D datasets. + pub scales: Vec, +} + +impl GroupDims { + /// The dimension defined by the scale at `address`. + pub fn by_address(&self, address: u64) -> Option<&Dimension> { + self.scales + .iter() + .find(|s| s.address == address) + .map(|s| &self.dims[s.dim]) + } + + /// The dimension whose scale has `_Netcdf4Dimid` `id`. + pub fn by_dimid(&self, id: i64) -> Option<&Dimension> { + self.scales + .iter() + .find(|s| s.dimid == Some(id)) + .map(|s| &self.dims[s.dim]) + } +} + +/// The dimensions of an HDF5 group (root or subgroup). /// /// NetCDF-4 stores dimensions as datasets with `CLASS=DIMENSION_SCALE`. A fixed /// dimension's size is the dataset's first (and typically only) shape extent. /// Unlimited dimensions have `max_dimensions[0] == u64::MAX` in the HDF5 dataspace; -/// their size is computed by `unlimited_len`. -pub(crate) fn extract_dimensions( +/// their size is computed by `unlimited_len`. A group with no dimension +/// scale at all (not written by a netCDF library) gets one dimension per +/// 1-D dataset instead. +pub(crate) fn group_dims( file: &clawhdf5::File, group: &clawhdf5::Group<'_>, -) -> Result, Error> { +) -> Result { + let addresses: HashMap = group.entries()?.into_iter().collect(); let dataset_names = group.datasets()?; - let mut dims = Vec::new(); - let mut seen_dimids: HashMap = HashMap::new(); + // (dimid, dimension, scale address), in discovery order. + let mut found: Vec<(Option, Dimension, u64)> = Vec::new(); for ds_name in &dataset_names { let ds = group.dataset(ds_name)?; let attrs = ds.attrs()?; - - // Check if this is a dimension scale if !is_dimension_scale(&attrs) { continue; } - + let Some(&address) = addresses.get(ds_name) else { + continue; + }; let shape = ds.shape()?; - let is_unlimited = check_unlimited(file, group, ds_name); + let is_unlimited = is_unlimited(&ds); let size = if is_unlimited { unlimited_len(file, &attrs, &shape) } else { shape.first().copied().unwrap_or(0) }; - - let dimid = get_dimid(&attrs); - let dim = Dimension { name: ds_name.clone(), size, is_unlimited, }; - - if let Some(id) = dimid { - seen_dimids.insert(id, dims.len()); - } - dims.push(dim); + found.push((get_dimid(&attrs), dim, address)); } - // Sort by dimid if available, otherwise keep discovery order - if !seen_dimids.is_empty() { - let mut pairs: Vec<(i64, Dimension)> = Vec::new(); - let mut unordered = Vec::new(); - - for (i, dim) in dims.into_iter().enumerate() { - let id = seen_dimids - .iter() - .find(|(_, idx)| **idx == i) - .map(|(k, _)| *k); - if let Some(id) = id { - pairs.push((id, dim)); - } else { - unordered.push(dim); + if found.is_empty() { + // Fallback: infer dimensions from dataset shapes and names. + // In NetCDF-4, coordinate variables are datasets whose name matches + // a dimension name. If there are no explicit DIMENSION_SCALE attributes, + // we look for 1-D datasets that might be coordinate variables. + let mut dims = Vec::new(); + for ds_name in &dataset_names { + let ds = group.dataset(ds_name)?; + let shape = ds.shape()?; + if shape.len() == 1 { + dims.push(Dimension { + name: ds_name.clone(), + size: shape[0], + is_unlimited: is_unlimited(&ds), + }); } } - pairs.sort_by_key(|(id, _)| *id); - dims = pairs.into_iter().map(|(_, d)| d).collect(); - dims.extend(unordered); + return Ok(GroupDims { + dims, + scales: Vec::new(), + }); } - Ok(dims) + // By dimid; scales without one keep their discovery order after those + // with one (the sort is stable). + found.sort_by_key(|(id, ..)| (id.is_none(), id.unwrap_or(0))); + let mut out = GroupDims::default(); + for (i, (dimid, dim, address)) in found.into_iter().enumerate() { + out.dims.push(dim); + out.scales.push(Scale { + address, + dimid, + dim: i, + }); + } + Ok(out) } /// The start of the `NAME` attribute netCDF-C gives a dimension scale that @@ -107,10 +157,7 @@ const PURE_DIMENSION_NAME: &str = "This is a netCDF dimension but not a netCDF v /// scale's own extent, as before. fn unlimited_len(file: &clawhdf5::File, attrs: &HashMap, shape: &[u64]) -> u64 { let own = shape.first().copied().unwrap_or(0); - let is_variable = !matches!( - attrs.get("NAME"), - Some(AttrValue::String(n)) if n.starts_with(PURE_DIMENSION_NAME) - ); + let is_variable = !is_pure_dimension(attrs); let Some(refs) = reference_list(file, attrs) else { return own; }; @@ -171,8 +218,59 @@ fn reference_list( Some(addresses.into_iter().map(|r| r.address).zip(axes).collect()) } +/// The dimension scale attached to each axis of a variable, from its +/// `DIMENSION_LIST` attribute (HDF5 dimension scales: one variable-length +/// sequence of object references per axis) — the address of the scale +/// netCDF-C takes for the axis, or `None` for an axis with none. netCDF-C's +/// `dimscale_visitor` lets `H5DSiterate_scales` visit every scale attached +/// to the axis and keeps the last, so with several (h5py's `attach_scale` +/// twice) the last one is the axis's dimension. `None` overall when the +/// attribute is missing or not in that form. +pub(crate) fn dimension_list( + file: &clawhdf5::File, + attrs: &HashMap, +) -> Option>> { + use clawhdf5_format::data_read::read_object_references; + use clawhdf5_format::datatype::Datatype; + use clawhdf5_format::vl_data::VlResolver; + let Some(AttrValue::Raw { datatype, data, .. }) = attrs.get("DIMENSION_LIST") else { + return None; + }; + let Datatype::VariableLength { + is_string: false, + base_type, + .. + } = datatype + else { + return None; + }; + let sb = file.superblock(); + let base_size = usize::try_from(base_type.type_size()).ok()?; + let sequences = VlResolver::new_in(file.storage(), sb.offset_size, sb.length_size) + .sequences(data, base_size) + .ok()?; + sequences + .iter() + .map(|refs| { + let refs = read_object_references(refs, base_type, sb.offset_size).ok()?; + Some(refs.iter().rev().find(|r| !r.is_null()).map(|r| r.address)) + }) + .collect() +} + +/// Whether a dataset is a dimension scale that is only a dimension, not a +/// netCDF variable: netCDF-C and h5netcdf give it this `NAME`, and netCDF-C +/// does not list it among the variables. +pub(crate) fn is_pure_dimension(attrs: &HashMap) -> bool { + is_dimension_scale(attrs) + && matches!( + attrs.get("NAME"), + Some(AttrValue::String(n)) if n.starts_with(PURE_DIMENSION_NAME) + ) +} + /// Check if a dataset's attributes mark it as a dimension scale. -fn is_dimension_scale(attrs: &HashMap) -> bool { +pub(crate) fn is_dimension_scale(attrs: &HashMap) -> bool { if let Some(AttrValue::String(class)) = attrs.get("CLASS") { return class == "DIMENSION_SCALE"; } @@ -180,7 +278,7 @@ fn is_dimension_scale(attrs: &HashMap) -> bool { } /// Get the _Netcdf4Dimid attribute value if present. -fn get_dimid(attrs: &HashMap) -> Option { +pub(crate) fn get_dimid(attrs: &HashMap) -> Option { match attrs.get("_Netcdf4Dimid") { Some(AttrValue::I64(id)) => Some(*id), Some(AttrValue::U64(id)) => Some(*id as i64), @@ -188,56 +286,8 @@ fn get_dimid(attrs: &HashMap) -> Option { } } -/// Check if a dimension is unlimited by inspecting the HDF5 dataspace max_dimensions. -/// -/// A dimension is unlimited when `max_dimensions[0] == u64::MAX` in the HDF5 dataspace. -fn check_unlimited(_file: &clawhdf5::File, group: &clawhdf5::Group<'_>, ds_name: &str) -> bool { - let ds = match group.dataset(ds_name) { - Ok(ds) => ds, - Err(_) => return false, - }; - - match ds.max_dimensions() { - Ok(Some(max_dims)) => max_dims.first().copied() == Some(u64::MAX), - _ => false, - } -} - -/// Extract dimensions from an HDF5 group using both dimension scale attributes -/// and variable DIMENSION_LIST references. -/// -/// This is a more robust approach that also discovers dimensions from variables -/// that reference them, even when dimension scales aren't explicitly set. -pub(crate) fn extract_dimensions_from_datasets( - group: &clawhdf5::Group<'_>, - file: &clawhdf5::File, -) -> Result, Error> { - // First try the standard approach with DIMENSION_SCALE - let mut dims = extract_dimensions(file, group)?; - - // If we found dimensions, return them - if !dims.is_empty() { - return Ok(dims); - } - - // Fallback: infer dimensions from dataset shapes and names. - // In NetCDF-4, coordinate variables are datasets whose name matches - // a dimension name. If there are no explicit DIMENSION_SCALE attributes, - // we look for 1-D datasets that might be coordinate variables. - let dataset_names = group.datasets()?; - for ds_name in &dataset_names { - let ds = group.dataset(ds_name)?; - let shape = ds.shape()?; - if shape.len() == 1 { - // This 1-D dataset could be a coordinate variable / dimension - let is_unlimited = check_unlimited(file, group, ds_name); - dims.push(Dimension { - name: ds_name.clone(), - size: shape[0], - is_unlimited, - }); - } - } - - Ok(dims) +/// Whether a dataset's first axis is unlimited (`max_dimensions[0] == +/// u64::MAX` in its dataspace). +fn is_unlimited(ds: &clawhdf5::Dataset<'_>) -> bool { + matches!(ds.max_dimensions(), Ok(Some(max_dims)) if max_dims.first() == Some(&u64::MAX)) } diff --git a/crates/clawhdf5-netcdf4/src/group.rs b/crates/clawhdf5-netcdf4/src/group.rs index a280435..4178ee5 100644 --- a/crates/clawhdf5-netcdf4/src/group.rs +++ b/crates/clawhdf5-netcdf4/src/group.rs @@ -9,12 +9,15 @@ use clawhdf5::AttrValue; use crate::dimension::{self, Dimension}; use crate::error::Error; -use crate::variable::{self, Variable}; +use crate::scope::{self, Scope}; +use crate::variable::Variable; /// A NetCDF-4 group corresponding to an HDF5 group. pub struct NetCDF4Group<'f> { /// Group name. name: String, + /// Path of the group from the root (`/`-separated). + path: String, /// Underlying HDF5 file. file: &'f clawhdf5::File, /// Underlying HDF5 group. @@ -25,11 +28,13 @@ impl<'f> NetCDF4Group<'f> { /// Create a new NetCDF4Group from an HDF5 group. pub(crate) fn new( name: String, + path: String, file: &'f clawhdf5::File, hdf5_group: clawhdf5::Group<'f>, ) -> Self { Self { name, + path, file, hdf5_group, } @@ -40,27 +45,22 @@ impl<'f> NetCDF4Group<'f> { &self.name } - /// List dimensions defined in this group. + /// List dimensions defined in this group (not those of its parent + /// groups, which its variables can also use). pub fn dimensions(&self) -> Result, Error> { - dimension::extract_dimensions_from_datasets(&self.hdf5_group, self.file) + Ok(dimension::group_dims(self.file, &self.hdf5_group)?.dims) } - /// List variables in this group. + /// List variables in this group: its datasets, except the dimension + /// scales that are only dimensions. Their dimensions can be defined in + /// this group or a parent group. pub fn variables(&self) -> Result>, Error> { - let dims = self.dimensions()?; - variable::build_variables(&self.hdf5_group, &dims) + Scope::new(self.file, &self.path)?.variables() } /// Get a specific variable by name. pub fn variable(&self, name: &str) -> Result, Error> { - let dims = self.dimensions()?; - let ds = self - .hdf5_group - .dataset(name) - .map_err(|_| Error::VariableNotFound(name.to_string()))?; - let shape = ds.shape()?; - let var_dims = crate::variable::match_dimensions_to_variable(&shape, &dims); - Ok(Variable::new(name.to_string(), ds, var_dims)) + scope::variable_at(self.file, &self.path, name) } /// Read all attributes of this group. @@ -79,12 +79,18 @@ impl<'f> NetCDF4Group<'f> { .hdf5_group .group(name) .map_err(|_| Error::GroupNotFound(name.to_string()))?; - Ok(NetCDF4Group::new(name.to_string(), self.file, hdf5_group)) + Ok(NetCDF4Group::new( + name.to_string(), + format!("{}/{name}", self.path), + self.file, + hdf5_group, + )) } - /// List dataset (variable) names in this group. + /// The names of this group's variables (see + /// [`variables`](Self::variables)). pub fn variable_names(&self) -> Result, Error> { - Ok(self.hdf5_group.datasets()?) + Scope::new(self.file, &self.path)?.variable_names() } } diff --git a/crates/clawhdf5-netcdf4/src/lib.rs b/crates/clawhdf5-netcdf4/src/lib.rs index 90ed40f..0e22773 100644 --- a/crates/clawhdf5-netcdf4/src/lib.rs +++ b/crates/clawhdf5-netcdf4/src/lib.rs @@ -27,6 +27,7 @@ pub mod cf; pub mod dimension; pub mod error; pub mod group; +mod scope; pub mod types; pub mod variable; @@ -75,25 +76,25 @@ impl NetCDF4File { /// List dimensions defined in the root group. pub fn dimensions(&self) -> Result, Error> { - dimension::extract_dimensions_from_datasets(&self.hdf5.root(), &self.hdf5) + Ok(dimension::group_dims(&self.hdf5, &self.hdf5.root())?.dims) } - /// List all variables in the root group. + /// List all variables in the root group: its datasets, except the + /// dimension scales that are only dimensions (netCDF-C does not list + /// them either). pub fn variables(&self) -> Result>, Error> { - let dims = self.dimensions()?; - variable::build_variables(&self.hdf5.root(), &dims) + scope::Scope::new(&self.hdf5, "/")?.variables() + } + + /// The names of the root group's variables (see + /// [`variables`](Self::variables)). + pub fn variable_names(&self) -> Result, Error> { + scope::Scope::new(&self.hdf5, "/")?.variable_names() } /// Get a specific variable by name from the root group. pub fn variable(&self, name: &str) -> Result, Error> { - let dims = self.dimensions()?; - let ds = self - .hdf5 - .dataset(name) - .map_err(|_| Error::VariableNotFound(name.to_string()))?; - let shape = ds.shape()?; - let var_dims = variable::match_dimensions_to_variable(&shape, &dims); - Ok(Variable::new(name.to_string(), ds, var_dims)) + scope::variable_at(&self.hdf5, "", name) } /// Read all global (root group) attributes. @@ -112,7 +113,12 @@ impl NetCDF4File { .hdf5 .group(name) .map_err(|_| Error::GroupNotFound(name.to_string()))?; - Ok(NetCDF4Group::new(name.to_string(), &self.hdf5, hdf5_group)) + Ok(NetCDF4Group::new( + name.to_string(), + name.to_string(), + &self.hdf5, + hdf5_group, + )) } /// Access the underlying HDF5 file for advanced operations. diff --git a/crates/clawhdf5-netcdf4/src/scope.rs b/crates/clawhdf5-netcdf4/src/scope.rs new file mode 100644 index 0000000..d74213b --- /dev/null +++ b/crates/clawhdf5-netcdf4/src/scope.rs @@ -0,0 +1,231 @@ +//! A group's variables and the dimensions they are defined on. +//! +//! netCDF-C (`libhdf5/hdf5open.c`) gives a variable its dimensions from the +//! file, never by size: the dimension ids in its `_Netcdf4Coordinates` +//! attribute (each dimension scale's `_Netcdf4Dimid`), else the dimension +//! scales its `DIMENSION_LIST` attribute references, looked up in the +//! variable's group and then each parent group up to the root. Only an axis +//! with neither (a file not written by a netCDF library) gets a dimension +//! by size. Dimension scales that are only dimensions are not variables, and +//! a variable stored as `_nc4_non_coord_` (a variable sharing a +//! dimension's name without being its coordinate variable) is ``. + +use std::collections::{HashMap, HashSet}; + +use clawhdf5::AttrValue; + +use crate::dimension::{self, Dimension, GroupDims}; +use crate::error::Error; +use crate::variable::Variable; + +/// The prefix netCDF-C gives the dataset of a variable that has a +/// dimension's name but is not that dimension's coordinate variable (the +/// dimension's scale holds the name). +const NON_COORD_PREFIX: &str = "_nc4_non_coord_"; + +/// A group, with the dimensions visible from it. +pub(crate) struct Scope<'f> { + file: &'f clawhdf5::File, + group: clawhdf5::Group<'f>, + /// This group's dimensions, then its parent's, and so on to the root's. + levels: Vec, +} + +impl<'f> Scope<'f> { + /// The group at `path` (`/`-separated from the root; `""` or `"/"` is + /// the root). + pub fn new(file: &'f clawhdf5::File, path: &str) -> Result { + let parts: Vec<&str> = path.split('/').filter(|p| !p.is_empty()).collect(); + let mut levels = Vec::with_capacity(parts.len() + 1); + for n in (0..=parts.len()).rev() { + let group = file.group(&parts[..n].join("/"))?; + levels.push(dimension::group_dims(file, &group)?); + } + let group = file.group(&parts.join("/"))?; + Ok(Self { + file, + group, + levels, + }) + } + + /// The group's datasets, as `(dataset name, object header address)` in + /// listing order. + fn datasets(&self) -> Result, Error> { + let datasets: HashSet = self.group.datasets()?.into_iter().collect(); + Ok(self + .group + .entries()? + .into_iter() + .filter(|(name, _)| datasets.contains(name)) + .collect()) + } + + /// The group's variables: every dataset but the dimension scales that + /// are only dimensions. + pub fn variables(&self) -> Result>, Error> { + let mut variables = Vec::new(); + for (ds_name, address) in self.datasets()? { + let ds = self.file.dataset_at(address)?; + let attrs = ds.attrs()?; + if dimension::is_pure_dimension(&attrs) { + continue; + } + variables.push(self.variable_from(nc_name(&ds_name), address, ds, attrs)?); + } + Ok(variables) + } + + /// The names of the group's variables. + pub fn variable_names(&self) -> Result, Error> { + let mut names = Vec::new(); + for (ds_name, address) in self.datasets()? { + let attrs = self.file.dataset_at(address)?.attrs()?; + if !dimension::is_pure_dimension(&attrs) { + names.push(nc_name(&ds_name)); + } + } + Ok(names) + } + + /// The variable called `name`: the dataset `_nc4_non_coord_` if + /// there is one, else the dataset `` unless it is only a + /// dimension. + pub fn variable(&self, name: &str) -> Result, Error> { + let not_found = || Error::VariableNotFound(name.to_string()); + let datasets = self.datasets()?; + let prefixed = format!("{NON_COORD_PREFIX}{name}"); + let address = datasets + .iter() + .find(|(n, _)| *n == prefixed) + .or_else(|| datasets.iter().find(|(n, _)| n == name)) + .map(|&(_, address)| address) + .ok_or_else(not_found)?; + let ds = self.file.dataset_at(address)?; + let attrs = ds.attrs()?; + if dimension::is_pure_dimension(&attrs) { + return Err(not_found()); + } + self.variable_from(nc_name(name), address, ds, attrs) + } + + fn variable_from( + &self, + name: String, + address: u64, + ds: clawhdf5::Dataset<'f>, + attrs: HashMap, + ) -> Result, Error> { + let shape = ds.shape()?; + let dims = self.variable_dims(address, &attrs, &shape); + Ok(Variable::new(name, ds, dims, attrs)) + } + + /// The first dimension, searching this group and then its ancestors, + /// that `find` picks. + fn find<'a>( + &'a self, + find: impl Fn(&'a GroupDims) -> Option<&'a Dimension>, + ) -> Option { + self.levels.iter().find_map(find).cloned() + } + + /// The dimensions of the dataset at `address`, one per axis of `shape`, + /// as netCDF-C resolves them (see the module docs). + fn variable_dims( + &self, + address: u64, + attrs: &HashMap, + shape: &[u64], + ) -> Vec { + let rank = shape.len(); + let mut dims: Vec> = vec![None; rank]; + if rank == 0 { + return Vec::new(); + } + // A coordinate variable is the scale of its (first) dimension. + dims[0] = self.levels[0].by_address(address).cloned(); + + if let Some(ids) = coordinates(attrs).filter(|ids| ids.len() == rank) { + for (slot, id) in dims.iter_mut().zip(ids) { + if slot.is_none() { + *slot = self.find(|level| level.by_dimid(id)); + } + } + } + if dims.iter().any(Option::is_none) + && let Some(scales) = + dimension::dimension_list(self.file, attrs).filter(|s| s.len() == rank) + { + for (slot, scale) in dims.iter_mut().zip(scales) { + if slot.is_none() + && let Some(scale) = scale + { + *slot = self.find(|level| level.by_address(scale)); + } + } + } + + // Neither: the first dimension of this group of the same size not + // already taken by another such axis, else an anonymous one. + let own = &self.levels[0].dims; + let mut used = vec![false; own.len()]; + dims.into_iter() + .zip(shape) + .map(|(dim, &size)| { + dim.unwrap_or_else(|| { + match own + .iter() + .enumerate() + .find(|&(i, d)| !used[i] && d.size == size) + { + Some((i, d)) => { + used[i] = true; + d.clone() + } + None => Dimension { + name: format!("dim_{size}"), + size, + is_unlimited: false, + }, + } + }) + }) + .collect() + } +} + +/// The variable `name` of the group at `group_path`; `name` may itself be +/// a path (`"sub/var"`), relative to that group. +pub(crate) fn variable_at<'f>( + file: &'f clawhdf5::File, + group_path: &str, + name: &str, +) -> Result, Error> { + match name.trim_start_matches('/').rsplit_once('/') { + Some((dir, leaf)) => Scope::new(file, &format!("{group_path}/{dir}")) + .map_err(|_| Error::VariableNotFound(name.to_string()))? + .variable(leaf), + None => Scope::new(file, group_path)?.variable(name.trim_start_matches('/')), + } +} + +/// The netCDF name of the dataset `ds_name`. +fn nc_name(ds_name: &str) -> String { + ds_name + .strip_prefix(NON_COORD_PREFIX) + .unwrap_or(ds_name) + .to_string() +} + +/// A variable's `_Netcdf4Coordinates`: the `_Netcdf4Dimid` of the dimension +/// of each axis. +fn coordinates(attrs: &HashMap) -> Option> { + match attrs.get("_Netcdf4Coordinates")? { + AttrValue::I64Array(ids) => Some(ids.clone()), + AttrValue::I64(id) => Some(vec![*id]), + AttrValue::U64Array(ids) => ids.iter().map(|&id| i64::try_from(id).ok()).collect(), + AttrValue::U64(id) => Some(vec![i64::try_from(*id).ok()?]), + _ => None, + } +} diff --git a/crates/clawhdf5-netcdf4/src/variable.rs b/crates/clawhdf5-netcdf4/src/variable.rs index ddb0b68..51d529e 100644 --- a/crates/clawhdf5-netcdf4/src/variable.rs +++ b/crates/clawhdf5-netcdf4/src/variable.rs @@ -2,12 +2,18 @@ //! //! Variables in NetCDF-4 are HDF5 datasets. This module wraps them with //! dimension associations and CF attribute support. +//! +//! A variable along an unlimited dimension has that dimension's length in +//! netCDF, even when fewer records of it have been written (its HDF5 dataset +//! is shorter): [`Variable::shape`] is the netCDF shape and the reads return +//! that many values, the unwritten ones as the fill value, as netCDF-C does. +//! [`Variable::stored_shape`] is the dataset's extent. use std::collections::HashMap; use clawhdf5::AttrValue; -use crate::cf::{self, CfAttributes}; +use crate::cf::{self, CfAttributes, FillValue}; use crate::dimension::Dimension; use crate::error::Error; use crate::types::{NcType, dtype_to_nctype}; @@ -20,18 +26,23 @@ pub struct Variable<'f> { dataset: clawhdf5::Dataset<'f>, /// Dimensions associated with this variable. dims: Vec, - /// Cached attributes. - attrs_cache: Option>, + /// The dataset's attributes. + attrs: HashMap, } impl<'f> Variable<'f> { /// Create a new Variable wrapping an HDF5 dataset. - pub(crate) fn new(name: String, dataset: clawhdf5::Dataset<'f>, dims: Vec) -> Self { + pub(crate) fn new( + name: String, + dataset: clawhdf5::Dataset<'f>, + dims: Vec, + attrs: HashMap, + ) -> Self { Self { name, dataset, dims, - attrs_cache: None, + attrs, } } @@ -40,13 +51,27 @@ impl<'f> Variable<'f> { &self.name } - /// The dimensions of this variable. + /// The dimensions of this variable, one per axis: the ones the file + /// gives it (`_Netcdf4Coordinates`, else `DIMENSION_LIST`), found in its + /// group or a parent group. An axis the file gives no dimension (a file + /// not written by a netCDF library) gets the first dimension of the + /// variable's group of the same size, else an anonymous `dim_`. pub fn dimensions(&self) -> &[Dimension] { &self.dims } - /// The shape of this variable (dimension sizes). + /// The shape of this variable as netCDF reports it: along an unlimited + /// dimension, the dimension's current length (the longest variable on + /// it), even if fewer records of this variable have been written; + /// otherwise the dataset's extent. The reads return this many values. pub fn shape(&self) -> Result, Error> { + Ok(nc_shape(&self.dataset.shape()?, &self.dims)) + } + + /// The extent of the HDF5 dataset: what has been written. It differs + /// from [`shape`](Self::shape) only along an unlimited dimension that + /// another variable has more records of. + pub fn stored_shape(&self) -> Result, Error> { Ok(self.dataset.shape()?) } @@ -58,59 +83,88 @@ impl<'f> Variable<'f> { /// Read all attributes as a HashMap. pub fn attrs(&mut self) -> Result<&HashMap, Error> { - if self.attrs_cache.is_none() { - self.attrs_cache = Some(self.dataset.attrs()?); - } - Ok(self - .attrs_cache - .as_ref() - .expect("invariant: attrs_cache is Some after initialization")) + Ok(&self.attrs) } /// Extract CF convention attributes. pub fn cf_attributes(&mut self) -> Result { - let attrs = self.attrs()?; - Ok(cf::extract_cf_attributes(attrs)) + Ok(cf::extract_cf_attributes(&self.attrs)) } /// Read data as f64 with scale_factor/add_offset applied. /// /// Missing values (matching `_FillValue` or `missing_value`) become NaN. /// If no scale_factor or add_offset attributes exist, returns the raw f64 data. + /// Records along an unlimited dimension that this variable has not + /// written (see [`shape`](Self::shape)) are NaN. pub fn read_f64(&mut self) -> Result, Error> { let raw = self.dataset.read_f64()?; - let cf = self.cf_attributes()?; - Ok(cf::apply_scale_offset(&raw, &cf)) + let cf = cf::extract_cf_attributes(&self.attrs); + self.padded(cf::apply_scale_offset(&raw, &cf), || Ok(f64::NAN)) } /// Read raw data as f64 without any scale/offset transformation. + /// + /// Unwritten records along an unlimited dimension read as the fill + /// value (`_FillValue`, else netCDF's default for the type), as in the + /// other `read_raw_*` methods and [`read_string`](Self::read_string). pub fn read_raw_f64(&self) -> Result, Error> { - Ok(self.dataset.read_f64()?) + self.padded_read(self.dataset.read_f64()?) } /// Read raw data as f32 without any scale/offset transformation. pub fn read_raw_f32(&self) -> Result, Error> { - Ok(self.dataset.read_f32()?) + self.padded_read(self.dataset.read_f32()?) } /// Read raw data as i32 without any scale/offset transformation. pub fn read_raw_i32(&self) -> Result, Error> { - Ok(self.dataset.read_i32()?) + self.padded_read(self.dataset.read_i32()?) } /// Read raw data as i64 without any scale/offset transformation. pub fn read_raw_i64(&self) -> Result, Error> { - Ok(self.dataset.read_i64()?) + self.padded_read(self.dataset.read_i64()?) } /// Read raw data as u64 without any scale/offset transformation. pub fn read_raw_u64(&self) -> Result, Error> { - Ok(self.dataset.read_u64()?) + self.padded_read(self.dataset.read_u64()?) } /// Read raw data as strings. pub fn read_string(&self) -> Result, Error> { - Ok(self.dataset.read_string()?) + self.padded_read(self.dataset.read_string()?) + } + + /// The fill value netCDF-C gives the variable's unwritten values: its + /// `_FillValue`, else the default fill value of its type (`NC_FILL_*`). + fn fill_value(&self) -> Result { + if let Some(fill) = cf::extract_cf_attributes(&self.attrs).fill_value { + return Ok(fill); + } + Ok(default_fill(self.nc_type()?)) + } + + /// `data`, read in the dataset's extent, laid out in the variable's + /// netCDF shape with the fill value in the positions not written. + fn padded_read(&self, data: Vec) -> Result, Error> { + self.padded(data, || Ok(T::from_fill(&self.fill_value()?))) + } + + /// Like [`padded_read`](Self::padded_read), padding with what `fill` + /// returns (called only when there is something to pad). + fn padded( + &self, + data: Vec, + fill: impl FnOnce() -> Result, + ) -> Result, Error> { + let extent = self.dataset.shape()?; + let shape = nc_shape(&extent, &self.dims); + if shape == extent { + return Ok(data); + } + pad(data, &extent, &shape, fill()?) } /// Read raw bytes without any type conversion. @@ -122,23 +176,23 @@ impl<'f> Variable<'f> { let dtype = self.dataset.dtype()?; match dtype { clawhdf5::DType::F64 => { - let vals = self.dataset.read_f64()?; + let vals = self.read_raw_f64()?; Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect()) } clawhdf5::DType::F32 => { - let vals = self.dataset.read_f32()?; + let vals = self.read_raw_f32()?; Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect()) } clawhdf5::DType::I32 => { - let vals = self.dataset.read_i32()?; + let vals = self.read_raw_i32()?; Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect()) } clawhdf5::DType::I64 => { - let vals = self.dataset.read_i64()?; + let vals = self.read_raw_i64()?; Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect()) } clawhdf5::DType::U64 => { - let vals = self.dataset.read_u64()?; + let vals = self.read_raw_u64()?; Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect()) } other => { @@ -169,64 +223,137 @@ impl std::fmt::Debug for Variable<'_> { } } -/// Build variables from a group's datasets and associated dimensions. -pub(crate) fn build_variables<'f>( - group: &clawhdf5::Group<'f>, - available_dims: &[Dimension], -) -> Result>, Error> { - let dataset_names = group.datasets()?; - let mut variables = Vec::new(); - - for ds_name in &dataset_names { - let ds = group.dataset(ds_name)?; - let shape = ds.shape()?; - - // Associate dimensions with this variable. - // First try DIMENSION_LIST attribute, then fall back to shape matching. - let var_dims = match_dimensions_to_variable(&shape, available_dims); - - variables.push(Variable::new(ds_name.clone(), ds, var_dims)); +/// The netCDF shape of a variable whose dataset has `extent`: along an +/// unlimited dimension the dimension's length, which is at least the extent. +fn nc_shape(extent: &[u64], dims: &[Dimension]) -> Vec { + if dims.len() != extent.len() { + return extent.to_vec(); } - - Ok(variables) + extent + .iter() + .zip(dims) + .map(|(&e, d)| if d.is_unlimited { e.max(d.size) } else { e }) + .collect() } -/// Match dimensions to a variable based on shape. -/// -/// For each axis of the variable, find a dimension with matching size. -/// If multiple dimensions have the same size, prefer exact name matching -/// from the convention order. -pub(crate) fn match_dimensions_to_variable( - shape: &[u64], - available_dims: &[Dimension], -) -> Vec { - let mut result = Vec::with_capacity(shape.len()); - - // Track which dimensions have been used to avoid duplicates - let mut used = vec![false; available_dims.len()]; - - for &dim_size in shape { - let mut matched = false; - - // Find a dimension with matching size that hasn't been used yet - for (i, dim) in available_dims.iter().enumerate() { - if !used[i] && dim.size == dim_size { - result.push(dim.clone()); - used[i] = true; - matched = true; +/// `data`, row-major in `extent`, placed in a row-major array of `shape` +/// (as many axes, each at least as long) filled with `fill`. +fn pad(data: Vec, extent: &[u64], shape: &[u64], fill: T) -> Result, Error> { + let too_big = || Error::TypeError(format!("variable of shape {shape:?} is too large")); + let to_usize = |dims: &[u64]| -> Result, Error> { + dims.iter() + .map(|&d| usize::try_from(d).map_err(|_| too_big())) + .collect() + }; + let (extent, shape) = (to_usize(extent)?, to_usize(shape)?); + let total = shape + .iter() + .try_fold(1usize, |n, &d| n.checked_mul(d)) + .ok_or_else(too_big)?; + if extent.len() != shape.len() + || extent.iter().zip(&shape).any(|(e, s)| e > s) + || extent.iter().product::() != data.len() + { + return Err(Error::TypeError(format!( + "{} values of extent {extent:?} do not fit shape {shape:?}", + data.len() + ))); + } + let (Some((&row, outer)), Some(&row_stride)) = (extent.split_last(), shape.last()) else { + return Ok(data); + }; + let mut out = vec![fill; total]; + if row == 0 { + return Ok(out); + } + // The position of the current row along each outer axis. + let mut index = vec![0usize; outer.len()]; + for chunk in data.chunks_exact(row) { + let offset = index.iter().zip(&shape).fold(0, |o, (&i, &n)| o * n + i); + out[offset * row_stride..][..row].clone_from_slice(chunk); + for (i, &n) in index.iter_mut().zip(outer).rev() { + *i += 1; + if *i < n { break; } - } - - if !matched { - // Create an anonymous dimension for unmatched sizes - result.push(Dimension { - name: format!("dim_{dim_size}"), - size: dim_size, - is_unlimited: false, - }); + *i = 0; } } + Ok(out) +} - result +/// netCDF's default fill value for a type (`NC_FILL_*` in `netcdf.h`). +fn default_fill(nc_type: NcType) -> FillValue { + match nc_type { + NcType::Byte => FillValue::Int(-127), + NcType::UByte => FillValue::UInt(255), + NcType::Short => FillValue::Int(-32767), + NcType::UShort => FillValue::UInt(65535), + NcType::Int => FillValue::Int(-2_147_483_647), + NcType::UInt => FillValue::UInt(4_294_967_295), + NcType::Int64 => FillValue::Int(-9_223_372_036_854_775_806), + NcType::UInt64 => FillValue::UInt(18_446_744_073_709_551_614), + NcType::Float => FillValue::Float(f64::from(9.969_21e36_f32)), + NcType::Double => FillValue::Float(9.969_209_968_386_869e36), + NcType::String => FillValue::String(String::new()), + NcType::Char => FillValue::Int(0), + } +} + +/// A fill value converted to the element type of a read, as the read +/// converts the stored values. +trait FromFill { + fn from_fill(fill: &FillValue) -> Self; +} + +macro_rules! numeric_from_fill { + ($($t:ty),*) => {$( + impl FromFill for $t { + fn from_fill(fill: &FillValue) -> Self { + match fill { + FillValue::Float(v) => *v as $t, + FillValue::Int(v) => *v as $t, + FillValue::UInt(v) => *v as $t, + FillValue::String(_) => <$t>::default(), + } + } + } + )*}; +} +numeric_from_fill!(f64, f32, i32, i64, u64); + +impl FromFill for String { + fn from_fill(fill: &FillValue) -> Self { + match fill { + FillValue::String(s) => s.clone(), + _ => String::new(), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn pad_places_rows() { + // (2, 1) written of (2, 4): each row padded, not the tail. + let out = pad(vec![1, 2], &[2, 1], &[2, 4], 0).unwrap(); + assert_eq!(out, vec![1, 0, 0, 0, 2, 0, 0, 0]); + // Leading axis short. + let out = pad(vec![1, 2, 3, 4], &[2, 2], &[3, 2], -1).unwrap(); + assert_eq!(out, vec![1, 2, 3, 4, -1, -1]); + // Nothing written. + let out = pad(Vec::::new(), &[0, 3], &[2, 3], 7).unwrap(); + assert_eq!(out, vec![7; 6]); + // 3-D, middle axis short. + let out = pad(vec![1, 2, 3, 4], &[2, 1, 2], &[2, 2, 2], 0).unwrap(); + assert_eq!(out, vec![1, 2, 0, 0, 3, 4, 0, 0]); + } + + #[test] + fn pad_rejects_wrong_length() { + assert!(pad(vec![1, 2, 3], &[2, 2], &[3, 2], 0).is_err()); + assert!(pad(vec![1, 2, 3, 4], &[2, 2], &[1, 4], 0).is_err()); + } } diff --git a/crates/clawhdf5-netcdf4/tests/interop_tests.rs b/crates/clawhdf5-netcdf4/tests/interop_tests.rs index ce3d5d0..2beffca 100644 --- a/crates/clawhdf5-netcdf4/tests/interop_tests.rs +++ b/crates/clawhdf5-netcdf4/tests/interop_tests.rs @@ -4,7 +4,7 @@ use std::process::Command; -use clawhdf5_netcdf4::{AttrValue, NetCDF4File}; +use clawhdf5_netcdf4::{AttrValue, NcType, NetCDF4File}; // --------------------------------------------------------------------------- // Helpers @@ -461,3 +461,384 @@ with nc.Dataset({path:?}) as f: } assert_eq!(got, expected); } + +// =========================================================================== +// Variables' dimensions, shapes and values as netCDF4-python reports them +// =========================================================================== + +/// Whether python can import `module`. +fn python_has(module: &str) -> bool { + Command::new(python()) + .args(["-c", &format!("import {module}")]) + .output() + .map(|o| o.status.success()) + .unwrap_or(false) +} + +/// h5netcdf is not in every interop environment (CI installs it; a local +/// `.venv` may not have it), so its tests skip without it even under +/// `CLAWHDF5_REQUIRE_INTEROP=1`. +macro_rules! skip_if_no_h5netcdf { + () => { + if !python_has("h5netcdf") { + eprintln!("SKIP: python3 with h5netcdf not available"); + return; + } + }; +} + +/// Every variable of the file at `path`, in every group, as netCDF4-python +/// reports it: `" () ()"` and its values +/// (numeric variables; element by element with masking off, so unwritten +/// records are the fill value), sorted by the description. +/// +/// Values are read one element at a time because netCDF-C 4.9.3 lays out a +/// whole-variable read of a variable shorter than an unlimited dimension +/// that is not its first wrongly (the written values first, then the fill); +/// element reads, and reads of one index of the leading axis, are right. +fn netcdf4_view(path: &std::path::Path) -> Vec<(String, Vec)> { + let script = r#" +import sys +import numpy as np +import netCDF4 as nc +def walk(g): + for name, v in g.variables.items(): + v.set_auto_mask(False) + head = "%s %s (%s) (%s)" % (g.path, name, ",".join(v.dimensions), ",".join(map(str, v.shape))) + vals = [] + if v.dtype != str and v.dtype.kind in "iuf": + vals = [repr(float(v[i])) for i in np.ndindex(v.shape)] + print(head + "|" + " ".join(vals)) + for sub in g.groups.values(): + walk(sub) +with nc.Dataset(sys.argv[1]) as f: + walk(f) +"#; + let out = Command::new(python()) + .args(["-c", script, &path.display().to_string()]) + .output() + .expect("failed to run python3"); + assert!( + out.status.success(), + "{}", + String::from_utf8_lossy(&out.stderr) + ); + let mut view: Vec<(String, Vec)> = String::from_utf8(out.stdout) + .unwrap() + .lines() + .map(|line| { + let (head, vals) = line.split_once('|').unwrap(); + let vals = vals + .split_whitespace() + .map(|v| v.parse().unwrap()) + .collect(); + (head.to_string(), vals) + }) + .collect(); + view.sort_by(|a, b| a.0.cmp(&b.0)); + view +} + +/// The same view of the file through clawhdf5-netcdf4. +fn clawhdf5_view(path: &std::path::Path) -> Vec<(String, Vec)> { + fn describe( + group_path: &str, + vars: Vec>, + ) -> Vec<(String, Vec)> { + vars.into_iter() + .map(|v| { + let dims: Vec<&str> = v.dimensions().iter().map(|d| d.name.as_str()).collect(); + let shape: Vec = v.shape().unwrap().iter().map(u64::to_string).collect(); + let head = format!( + "{group_path} {} ({}) ({})", + v.name(), + dims.join(","), + shape.join(",") + ); + let vals = match v.nc_type().unwrap() { + NcType::String | NcType::Char => Vec::new(), + _ => v.read_raw_f64().unwrap(), + }; + (head, vals) + }) + .collect() + } + fn walk( + group_path: &str, + group: &clawhdf5_netcdf4::NetCDF4Group<'_>, + out: &mut Vec<(String, Vec)>, + ) { + out.extend(describe(group_path, group.variables().unwrap())); + for name in group.group_names().unwrap() { + walk( + &format!("{group_path}/{name}"), + &group.group(&name).unwrap(), + out, + ); + } + } + let file = NetCDF4File::open(path).unwrap(); + let mut view = describe("/", file.variables().unwrap()); + for name in file.group_names().unwrap() { + walk(&format!("/{name}"), &file.group(&name).unwrap(), &mut view); + } + view.sort_by(|a, b| a.0.cmp(&b.0)); + view +} + +/// clawhdf5-netcdf4 reports the same variables, dimensions, shapes and +/// values (bit for bit, NaN equal to NaN) as netCDF4-python. +fn assert_same_view(path: &std::path::Path) { + let want = netcdf4_view(path); + let got = clawhdf5_view(path); + let heads = |v: &[(String, Vec)]| v.iter().map(|(h, _)| h.clone()).collect::>(); + assert_eq!(heads(&got), heads(&want), "variables differ from netCDF4's"); + for ((head, got), (_, want)) in got.iter().zip(&want) { + let same = got.len() == want.len() + && got + .iter() + .zip(want) + .all(|(a, b)| a.to_bits() == b.to_bits() || (a.is_nan() && b.is_nan())); + assert!(same, "{head}: got {got:?}, netCDF4 reads {want:?}"); + } +} + +/// The reproducer of the known-issues entry: `a` is on the unlimited `time` +/// (5 long through `b`) with 2 records, not on an anonymous `dim_2`; the +/// pure dimension scales `time` and `empty` are not variables; `a` has +/// shape (5,) and reads its 3 unwritten records as the fill value. +#[test] +fn variable_dimensions_come_from_the_file() { + skip_if_no_netcdf4!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("repro.nc"); + run_python(&format!( + r#" +import netCDF4 as nc +import numpy as np +with nc.Dataset({path:?}, "w") as f: + f.createDimension("time", None) + f.createDimension("empty", None) + f.createDimension("x", 3) + f.createVariable("a", "i4", ("time",))[0:2] = [1, 2] + f.createVariable("b", "f4", ("time", "x"))[0:5, :] = np.arange(15).reshape(5, 3) + f.createVariable("e", "i4", ("empty",)) + f.createVariable("c", "i4", ("x",))[:] = [7, 8, 9] +"#, + path = path.display().to_string() + )); + assert_same_view(&path); + + let file = NetCDF4File::open(&path).unwrap(); + let mut names = file.variable_names().unwrap(); + names.sort(); + assert_eq!(names, ["a", "b", "c", "e"]); + assert!(matches!( + file.variable("time"), + Err(clawhdf5_netcdf4::Error::VariableNotFound(_)) + )); + let a = file.variable("a").unwrap(); + assert_eq!(a.dimensions()[0].name, "time"); + assert_eq!(a.shape().unwrap(), [5]); + assert_eq!(a.stored_shape().unwrap(), [2]); + assert_eq!( + a.read_raw_i32().unwrap(), + [1, 2, -2_147_483_647, -2_147_483_647, -2_147_483_647] + ); +} + +/// Dimensions of one size are told apart by the file, not by order: `p` +/// and `q` are both 2 long, and `v(q, p)`, `same(p, p)` (one dimension +/// twice), a scalar, `q`'s coordinate variable, a variable called `p` that +/// is not `p`'s coordinate variable (stored as `_nc4_non_coord_p`), and +/// variables in a subgroup and a sub-subgroup on dimensions of their +/// ancestors. +#[test] +fn equal_size_and_inherited_dimensions_match_netcdf4_python() { + skip_if_no_netcdf4!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("dims.nc"); + run_python(&format!( + r#" +import netCDF4 as nc +import numpy as np +with nc.Dataset({path:?}, "w") as f: + f.createDimension("p", 2) + f.createDimension("q", 2) + f.createVariable("v", "i4", ("q", "p"))[:] = np.array([[1, 2], [3, 4]]) + f.createVariable("same", "i4", ("p", "p"))[:] = np.array([[5, 6], [7, 8]]) + f.createVariable("s", "f8", ())[...] = 3.5 + f.createVariable("q", "f4", ("q",))[:] = [0, 1] + f.createVariable("p", "f4", ("q", "p"))[:] = np.array([[0, 1], [2, 3]]) + g = f.createGroup("g") + g.createDimension("r", 2) + g.createVariable("w", "i4", ("r", "q", "p"))[:] = np.arange(8).reshape(2, 2, 2) + h = g.createGroup("h") + h.createVariable("z", "i4", ("p", "r"))[:] = np.array([[1, 2], [3, 4]]) +"#, + path = path.display().to_string() + )); + assert_same_view(&path); + + let file = NetCDF4File::open(&path).unwrap(); + let v = file.variable("v").unwrap(); + let dims: Vec<&str> = v.dimensions().iter().map(|d| d.name.as_str()).collect(); + assert_eq!(dims, ["q", "p"]); + let p = file.variable("p").unwrap(); + assert!(!p.is_coordinate()); + assert!(file.variable("q").unwrap().is_coordinate()); + let s = file.variable("s").unwrap(); + assert!(s.dimensions().is_empty()); + assert_eq!(s.shape().unwrap(), Vec::::new()); + let z = file + .group("g") + .unwrap() + .group("h") + .unwrap() + .variable("z") + .unwrap(); + let dims: Vec<&str> = z.dimensions().iter().map(|d| d.name.as_str()).collect(); + assert_eq!(dims, ["p", "r"]); +} + +/// Variables shorter than their unlimited dimension have its length and +/// read the fill value (`_FillValue`, else netCDF's default for the type) +/// where nothing was written — also when the unlimited dimension is not +/// the first; `read_f64` gives NaN there. +#[test] +fn unwritten_records_read_as_fill_like_netcdf4_python() { + skip_if_no_netcdf4!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("pad.nc"); + run_python(&format!( + r#" +import netCDF4 as nc +import numpy as np +with nc.Dataset({path:?}, "w") as f: + f.createDimension("t", None) + f.createDimension("x", 2) + f.createVariable("a", "i4", ("t",))[0:2] = [1, 2] + f.createVariable("f", "f4", ("x", "t"), fill_value=-5.0)[:, 0:1] = np.array([[1], [2]]) + f.createVariable("d", "f8", ("t",))[0:4] = [1, 2, 3, 4] + f.createVariable("u", "u8", ("t",))[0:1] = [1] + f.createVariable("b", "i1", ("t", "x"))[0:3, :] = np.ones((3, 2)) + f.createVariable("st", str, ("t",))[0] = "hi" + g = f.createGroup("g") + g.createVariable("k", "f4", ("t",))[0:1] = [9] +"#, + path = path.display().to_string() + )); + assert_same_view(&path); + + let file = NetCDF4File::open(&path).unwrap(); + let mut f = file.variable("f").unwrap(); + assert_eq!(f.shape().unwrap(), [2, 4]); + assert_eq!(f.stored_shape().unwrap(), [2, 1]); + assert_eq!( + f.read_raw_f32().unwrap(), + [1.0, -5.0, -5.0, -5.0, 2.0, -5.0, -5.0, -5.0] + ); + let read = f.read_f64().unwrap(); + assert_eq!(read[0], 1.0); + assert!(read[1].is_nan() && read[7].is_nan()); + let st = file.variable("st").unwrap(); + assert_eq!(st.read_string().unwrap(), ["hi", "", "", ""]); + assert_eq!(st.shape().unwrap(), [4]); +} + +/// A file with HDF5 dimension scales but none of netCDF's own attributes +/// (h5py's `dims` API): the dimensions come from `DIMENSION_LIST`, so +/// `v(q, p)` is not `v(p, q)` although both are 2 long; with two scales +/// attached to one axis (`w`), netCDF-C takes the last. +#[test] +fn h5py_dimension_scales_match_netcdf4_python() { + skip_if_no_netcdf4!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("scales.h5"); + run_python(&format!( + r#" +import h5py +import numpy as np +with h5py.File({path:?}, "w") as f: + f["p"] = np.arange(2.0) + f["q"] = np.arange(2.0) + 10 + f["p"].make_scale("p") + f["q"].make_scale("q") + f["v"] = np.arange(4).reshape(2, 2) + f["v"].dims[0].attach_scale(f["q"]) + f["v"].dims[1].attach_scale(f["p"]) + f["w"] = np.arange(2) + f["w"].dims[0].attach_scale(f["p"]) + f["w"].dims[0].attach_scale(f["q"]) +"#, + path = path.display().to_string() + )); + assert_same_view(&path); +} + +/// Files h5netcdf writes (its own implementation of the netCDF-4 +/// conventions over h5py): an unlimited dimension, equal sizes, a subgroup +/// on inherited dimensions, a scalar. +#[test] +fn h5netcdf_file_matches_netcdf4_python() { + skip_if_no_netcdf4!(); + skip_if_no_h5netcdf!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("h5netcdf.nc"); + run_python(&format!( + r#" +import h5netcdf +import numpy as np +with h5netcdf.File({path:?}, "w") as f: + f.dimensions = {{"p": 2, "q": 2, "t": None}} + f.create_variable("v", ("q", "p"), "i4")[...] = np.array([[1, 2], [3, 4]]) + f.create_variable("q", ("q",), "f4")[...] = [0, 1] + f.create_variable("same", ("p", "p"), "i4")[...] = np.array([[5, 6], [7, 8]]) + a = f.create_variable("a", ("t", "p"), "f8") + f.resize_dimension("t", 3) + a[...] = np.ones((3, 2)) + f.create_variable("short", ("t",), "i4") + g = f.create_group("g") + g.dimensions = {{"r": 2}} + g.create_variable("w", ("r", "q", "p"), "i4")[...] = np.arange(8).reshape(2, 2, 2) + g.create_variable("s", (), "f8")[...] = 2.5 +"#, + path = path.display().to_string() + )); + assert_same_view(&path); +} + +/// Files xarray writes, through netCDF4 and (when installed) h5netcdf: +/// coordinates, two dimensions of one size, an unlimited dimension. +#[test] +fn xarray_files_match_netcdf4_python() { + skip_if_no_netcdf4!(); + skip_if_no_xarray!(); + let dir = tempfile::tempdir().unwrap(); + let mut engines = vec!["netcdf4"]; + if python_has("h5netcdf") { + engines.push("h5netcdf"); + } else { + eprintln!("SKIP: xarray with engine h5netcdf (h5netcdf not available)"); + } + for engine in engines { + let path = dir.path().join(format!("xarray_{engine}.nc")); + run_python(&format!( + r#" +import numpy as np +import xarray as xr +ds = xr.Dataset( + {{ + "temp": (("time", "lat", "lon"), np.arange(12.0).reshape(3, 2, 2)), + "grid": (("lon", "lat"), np.array([[1, 2], [3, 4]], dtype="i4")), + "scalar": ((), 1.5), + }}, + coords={{"time": [0.0, 6.0, 12.0], "lat": [10.0, 20.0], "lon": [5.0, 6.0]}}, +) +ds.to_netcdf({path:?}, engine={engine:?}, unlimited_dims=["time"]) +"#, + path = path.display().to_string() + )); + assert_same_view(&path); + } +} diff --git a/crates/clawhdf5-netcdf4/tests/netcdf4_tests.rs b/crates/clawhdf5-netcdf4/tests/netcdf4_tests.rs index 877214e..a093759 100644 --- a/crates/clawhdf5-netcdf4/tests/netcdf4_tests.rs +++ b/crates/clawhdf5-netcdf4/tests/netcdf4_tests.rs @@ -754,3 +754,53 @@ fn test_dimension_struct_equality() { }; assert_ne!(d1, d3); } + +/// A dimension scale that is only a dimension (netCDF-C's `NAME`) is not a +/// variable, and `_nc4_non_coord_` is the variable ``, found in +/// place of the scale of the same name. +#[test] +fn test_pure_dimensions_hidden_and_non_coord_names() { + let pure = "This is a netCDF dimension but not a netCDF variable. 2"; + let mut b = FileBuilder::new(); + b.create_dataset("x") + .with_f32_data(&[0.0, 0.0]) + .with_shape(&[2]) + .set_attr("CLASS", AttrValue::String("DIMENSION_SCALE".into())) + .set_attr("NAME", AttrValue::String(pure.into())) + .set_attr("_Netcdf4Dimid", AttrValue::I64(0)); + b.create_dataset("_nc4_non_coord_x") + .with_f64_data(&[1.0, 2.0, 3.0]) + .with_shape(&[3]); + b.create_dataset("v") + .with_f64_data(&[5.0, 6.0]) + .with_shape(&[2]); + let file = NetCDF4File::from_bytes(b.finish().unwrap()).unwrap(); + + let dims = file.dimensions().unwrap(); + assert_eq!(dims.len(), 1); + assert_eq!(dims[0].name, "x"); + let mut names = file.variable_names().unwrap(); + names.sort(); + assert_eq!(names, ["v", "x"]); + let x = file.variable("x").unwrap(); + assert_eq!(x.name(), "x"); + assert_eq!(x.read_raw_f64().unwrap(), [1.0, 2.0, 3.0]); + assert!(!x.is_coordinate()); + // No DIMENSION_LIST: `v` gets `x` by size, as before. + assert_eq!(file.variable("v").unwrap().dimensions()[0].name, "x"); +} + +/// `variable` still takes a path relative to the group, as it did when it +/// opened the dataset by path. +#[test] +fn test_variable_by_path() { + let file = NetCDF4File::from_bytes(make_grouped_netcdf4()).unwrap(); + let pressure = file.variable("surface/pressure").unwrap(); + assert_eq!(pressure.name(), "pressure"); + assert_eq!(pressure.read_raw_f64().unwrap(), [1013.25, 1012.0, 1011.5]); + assert!(file.variable("/time").is_ok()); + assert!(matches!( + file.variable("nowhere/pressure"), + Err(clawhdf5_netcdf4::Error::VariableNotFound(_)) + )); +} diff --git a/crates/clawhdf5-tools/src/check.rs b/crates/clawhdf5-tools/src/check.rs index 4a36554..1ba8239 100644 --- a/crates/clawhdf5-tools/src/check.rs +++ b/crates/clawhdf5-tools/src/check.rs @@ -916,7 +916,7 @@ impl Checker<'_> { bad.push("has size 0".into()); } else if !filtered && let Some(cb) = chunk_bytes - && u64::from(c.chunk_size) != cb + && c.chunk_size != cb { bad.push(format!( "is {} bytes; an unfiltered chunk is {cb}", @@ -933,7 +933,7 @@ impl Checker<'_> { ); } } - self.extent(c.address, u64::from(c.chunk_size), path); + self.extent(c.address, c.chunk_size, path); } if reported > 50 { self.problem( diff --git a/crates/clawhdf5-tools/src/info.rs b/crates/clawhdf5-tools/src/info.rs index d169f73..b0239f1 100644 --- a/crates/clawhdf5-tools/src/info.rs +++ b/crates/clawhdf5-tools/src/info.rs @@ -224,7 +224,7 @@ pub fn allocated_bytes(h5: &H5, info: &DsInfo) -> Result { let ds = info.ds.as_ref().map_err(Clone::clone)?; chunks(h5, layout, ds, dt)? .iter() - .map(|c| u64::from(c.chunk_size)) + .map(|c| c.chunk_size) .sum() } DataLayout::Virtual { .. } => 0, diff --git a/crates/clawhdf5-tools/tests/libver_v18.rs b/crates/clawhdf5-tools/tests/libver_v18.rs new file mode 100644 index 0000000..0d27a39 --- /dev/null +++ b/crates/clawhdf5-tools/tests/libver_v18.rs @@ -0,0 +1,736 @@ +//! Files written with `FileBuilder::libver_bounds(LibVer::V18, LibVer::V18)` +//! must be readable by HDF5 1.8. Every writer feature is written under that +//! bound and read back by HDF5 1.8.23's h5dump (values dumped in binary and +//! compared), by h5py (libhdf5 2.x), by h5dump 1.14, by clawhdf5 and by +//! `h5rs check --data`; then `FileEditor` grows, appends to and annotates +//! the file (splitting version-1 B-tree nodes) and every reader checks it +//! again. The version-1 B-trees we write are compared node by node with +//! the ones libhdf5 writes for the same data under h5py's +//! `libver=('v108', 'latest')`. +//! +//! HDF5 1.8 is found through `CLAWHDF5_H5DUMP18` (the path to its h5dump) +//! or at `~/.cache/hdf5-1.8.23/bin/h5dump`, where +//! `scripts/build-hdf5-1.8.sh` builds it; without it the 1.8 checks are +//! skipped (CI has no HDF5 1.8), even with `CLAWHDF5_REQUIRE_INTEROP=1`. +//! h5py/numpy (`CLAWHDF5_PYTHON`) and h5dump are needed otherwise; they +//! skip when missing unless `CLAWHDF5_REQUIRE_INTEROP=1`. + +use std::path::{Path, PathBuf}; +use std::process::{Command, Output}; + +use clawhdf5::{AttrValue, File, FileBuilder, FileEditor, LibVer, Selection}; +use clawhdf5_format::datatype::{CharacterSet, Datatype, StringPadding}; +use clawhdf5_format::file_writer::{CompoundTypeBuilder, EnumTypeBuilder}; +use clawhdf5_format::type_builders::make_i32_type; + +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 available(cmd: &str, args: &[&str]) -> bool { + Command::new(cmd) + .args(args) + .output() + .map(|o| o.status.success()) + .unwrap_or(false) +} + +fn tools_ok() -> bool { + let ok = + available(&python(), &["-c", "import h5py, numpy"]) && available("h5dump", &["--version"]); + if !ok { + assert!( + !interop_required(), + "CLAWHDF5_REQUIRE_INTEROP=1 but h5py/numpy or h5dump is not available" + ); + eprintln!("SKIP: h5py/numpy or h5dump not available"); + } + ok +} + +/// HDF5 1.8's h5dump, when there is one. +fn h5dump18() -> Option { + let p = match std::env::var_os("CLAWHDF5_H5DUMP18") { + Some(p) => PathBuf::from(p), + None => PathBuf::from(std::env::var_os("HOME")?).join(".cache/hdf5-1.8.23/bin/h5dump"), + }; + let o = Command::new(&p).arg("--version").output().ok()?; + let v = String::from_utf8_lossy(&o.stdout).to_string(); + if !v.contains("1.8.") { + eprintln!("SKIP 1.8 checks: {} is not HDF5 1.8 ({v})", p.display()); + return None; + } + Some(p) +} + +fn py(script: &str) -> String { + let o = Command::new(python()) + .args(["-c", script]) + .output() + .expect("run python"); + assert!( + o.status.success(), + "python failed:\n{script}\nSTDOUT: {}\nSTDERR: {}", + String::from_utf8_lossy(&o.stdout), + String::from_utf8_lossy(&o.stderr) + ); + String::from_utf8_lossy(&o.stdout).trim().to_string() +} + +fn text(o: &Output) -> String { + format!( + "{}{}", + String::from_utf8_lossy(&o.stdout), + String::from_utf8_lossy(&o.stderr) + ) +} + +fn tmpdir() -> tempfile::TempDir { + tempfile::TempDir::new_in(env!("CARGO_TARGET_TMPDIR")).unwrap() +} + +/// The hyperslab of `count` elements from `start`. +fn block(start: &[u64], count: &[u64]) -> Selection { + Selection::Hyperslab { + start: start.to_vec(), + stride: vec![1; start.len()], + count: count.to_vec(), + block: vec![1; start.len()], + } +} + +fn le(v: &[T], f: impl Fn(T) -> [u8; N]) -> Vec { + v.iter().flat_map(|&x| f(x)).collect() +} + +/// The datasets of the test file and the bytes each holds (little-endian, +/// row-major, as h5py's `tobytes()` and h5dump's `-b LE` give them). +struct Expect { + datasets: Vec<(String, Vec)>, +} + +impl Expect { + fn set(&mut self, name: &str, bytes: Vec) { + match self.datasets.iter_mut().find(|(n, _)| n == name) { + Some(e) => e.1 = bytes, + None => self.datasets.push((name.to_string(), bytes)), + } + } +} + +const MANY: usize = 100_000; + +/// Write every feature under the 1.8 bound. +fn write_file(path: &Path) -> Expect { + let mut e = Expect { + datasets: Vec::new(), + }; + let mut b = FileBuilder::new(); + b.libver_bounds(LibVer::V18, LibVer::V18); + + // Contiguous, with dense attributes (more than 8). + let v: Vec = (0..1000).map(|i| i as f64 * 0.5).collect(); + let d = b.create_dataset("contig").with_f64_data(&v); + for i in 0..12 { + d.set_attr(&format!("a{i:02}"), AttrValue::I64(i)); + } + e.set("/contig", le(&v, f64::to_le_bytes)); + b.create_dataset("empty").with_f64_data(&[]); + e.set("/empty", vec![]); + b.create_dataset("scalar") + .with_f64_data(&[2.5]) + .with_shape(&[]); + e.set("/scalar", 2.5f64.to_le_bytes().to_vec()); + let v: Vec = (0..16).map(|i| i * 3 - 7).collect(); + b.create_dataset("compact").with_i32_data(&v).compact(); + e.set("/compact", le(&v, i32::to_le_bytes)); + let v: Vec = (0..20).map(|i| i as f32 / 3.0).collect(); + b.create_dataset("f32").with_f32_data(&v); + e.set("/f32", le(&v, f32::to_le_bytes)); + let v: Vec = vec![0.5, -2.0, 1024.0, 0.0]; + b.create_dataset("f16").with_f16_data(&v); + e.set( + "/f16", + le(&v, |x| { + clawhdf5_format::float16::f32_to_f16_bits(x).to_le_bytes() + }), + ); + + // Chunked with every built-in filter HDF5 1.8 has, and edge chunks. + let v: Vec = (0..60 * 70).map(|i| (i % 97) as f32 * 1.25).collect(); + b.create_dataset("chunked") + .with_f32_data(&v) + .with_shape(&[60, 70]) + .with_chunks(&[16, 16]) + .with_deflate(6) + .with_shuffle() + .with_fletcher32(); + e.set("/chunked", le(&v, f32::to_le_bytes)); + // What the 1.10 indexes would be: Extensible Array, version-2 B-tree, + // Fixed Array, single chunk. All become version-1 B-trees. + let v: Vec = (0..25).collect(); + b.create_dataset("resizable") + .with_i32_data(&v) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[4]); + e.set("/resizable", le(&v, i32::to_le_bytes)); + let v: Vec = (0..30).map(|i| i * 1_000_000_007).collect(); + b.create_dataset("resizable2") + .with_i64_data(&v) + .with_shape(&[5, 6]) + .with_maxshape(&[u64::MAX, u64::MAX]) + .with_chunks(&[2, 4]); + e.set("/resizable2", le(&v, i64::to_le_bytes)); + let v: Vec = (0..10).map(|i| -i).collect(); + b.create_dataset("fixedmax") + .with_i32_data(&v) + .with_maxshape(&[100]) + .with_chunks(&[3]) + .with_deflate(1); + e.set("/fixedmax", le(&v, i32::to_le_bytes)); + let v: Vec = (0..8).map(|i| i * i).collect(); + b.create_dataset("single") + .with_i32_data(&v) + .with_chunks(&[8]); + e.set("/single", le(&v, i32::to_le_bytes)); + // Enough chunks for a three-level tree, and a 2-D two-level one. + let v: Vec = (0..MANY as i32).map(|i| i ^ 0x5a5a).collect(); + b.create_dataset("many") + .with_i32_data(&v) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[1]); + e.set("/many", le(&v, i32::to_le_bytes)); + let v: Vec = (0..100 * 100).map(|i| i as f64).collect(); + b.create_dataset("grid") + .with_f64_data(&v) + .with_shape(&[100, 100]) + .with_chunks(&[1, 1]); + e.set("/grid", le(&v, f64::to_le_bytes)); + b.create_dataset("empty_chunked") + .with_f64_data(&[]) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[10]); + e.set("/empty_chunked", vec![]); + let v: Vec = (0..12).collect(); + b.create_dataset("filled") + .with_i32_data(&v) + .with_maxshape(&[u64::MAX]) + .with_chunks(&[5]) + .with_fill_value(&(-9i32).to_le_bytes()); + e.set("/filled", le(&v, i32::to_le_bytes)); + + // Datatypes. + let raw = b"abcdehello\0\0\0\0\0".to_vec(); + b.create_dataset("strings").with_compound_data( + Datatype::String { + size: 5, + padding: StringPadding::NullPad, + charset: CharacterSet::Ascii, + }, + raw.clone(), + 3, + ); + e.set("/strings", raw); + let ct = CompoundTypeBuilder::new() + .i32_field("a") + .f64_field("b") + .build(); + let mut raw = Vec::new(); + for i in 0..5i32 { + raw.extend_from_slice(&i.to_le_bytes()); + raw.extend_from_slice(&(f64::from(i) * 1.5).to_le_bytes()); + } + b.create_dataset("compound") + .with_compound_data(ct, raw.clone(), 5); + e.set("/compound", raw); + let et = EnumTypeBuilder::i32_based() + .value("RED", 0) + .value("GREEN", 1) + .value("BLUE", 7) + .build(); + let v = [0, 7, 1, 1, 0]; + b.create_dataset("enum").with_enum_i32_data(et, &v); + e.set("/enum", le(&v, i32::to_le_bytes)); + let v: Vec = (0..24).collect(); + b.create_dataset("array").with_array_data( + make_i32_type(), + &[2, 3], + le(&v, i32::to_le_bytes), + 4, + ); + e.set("/array", le(&v, i32::to_le_bytes)); + let v: Vec = vec![-1, 0, i64::MAX]; + b.create_dataset("i64").with_i64_data(&v); + e.set("/i64", le(&v, i64::to_le_bytes)); + + // Groups: compact with links of every kind, dense (links and + // attributes), and tracking creation order. + let mut g = b.create_group("g_compact"); + g.create_dataset("x").with_i32_data(&[1, 2, 3]); + e.set("/g_compact/x", le(&[1i32, 2, 3], i32::to_le_bytes)); + g.add_soft_link("soft", "/contig"); + g.add_external_link("ext", "other.h5", "/data"); + g.set_attr("title", AttrValue::String("compact group".into())); + b.add_group(g.finish()); + b.add_hard_link("hard", "/g_compact/x"); + let mut g = b.create_group("g_dense"); + for i in 0..20 { + let v = [i, i + 1]; + g.create_dataset(&format!("d{i:02}")).with_i32_data(&v); + e.set(&format!("/g_dense/d{i:02}"), le(&v, i32::to_le_bytes)); + } + for i in 0..12 { + g.set_attr(&format!("attr{i:02}"), AttrValue::F64(f64::from(i) / 4.0)); + } + b.add_group(g.finish()); + let mut g = b.create_group("g_order"); + g.track_order(true); + for i in 0..10 { + let n = format!("z{}", 9 - i); + g.create_dataset(&n).with_i32_data(&[i]); + e.set(&format!("/g_order/{n}"), i.to_le_bytes().to_vec()); + } + for i in 0..10 { + g.set_attr(&format!("b{}", 9 - i), AttrValue::I64(i)); + } + b.add_group(g.finish()); + b.create_dataset("deep/er/path").with_i32_data(&[42]); + e.set("/deep/er/path", 42i32.to_le_bytes().to_vec()); + + b.set_attr("version", AttrValue::F64(1.8)); + b.set_attr("ints", AttrValue::I64Array(vec![1, -2, 3])); + b.set_attr("text", AttrValue::String("readable by 1.8".into())); + b.set_attr( + "texts", + AttrValue::StringArray(vec!["a".into(), "bb".into(), "ccc".into()]), + ); + b.set_attr("big", AttrValue::U64(u64::MAX)); + b.write(path).unwrap(); + e +} + +/// Structural facts of a file: superblock and layout message versions. +fn check_versions(path: &Path) { + use clawhdf5_format::message_type::MessageType; + use clawhdf5_format::object_header::ObjectHeader; + let bytes = std::fs::read(path).unwrap(); + assert_eq!(bytes[8], 2, "superblock version"); + let f = File::open(path).unwrap(); + let sb = f.superblock().clone(); + for name in [ + "contig", + "compact", + "chunked", + "resizable", + "resizable2", + "many", + ] { + let addr = clawhdf5_format::group_v2::resolve_path_any(&bytes, &sb, name).unwrap(); + let oh = ObjectHeader::parse(&bytes, addr as usize, 8, 8).unwrap(); + let layout = oh + .messages + .iter() + .find(|m| m.msg_type == MessageType::DataLayout) + .unwrap(); + assert_eq!(layout.data[0], 3, "layout version of {name}"); + } +} + +/// clawhdf5 reads every dataset's bytes back. +fn check_ours(path: &Path, e: &Expect) { + let f = File::open(path).unwrap(); + for (name, want) in &e.datasets { + let got = f + .dataset(name) + .unwrap() + .read_selection(&Selection::All) + .unwrap(); + assert!(got == *want, "our read of {name} differs"); + } + let attrs = f.dataset("contig").unwrap().attrs().unwrap(); + assert!((12..=13).contains(&attrs.len()), "{attrs:?}"); +} + +/// h5py (libhdf5 2.x) reads every dataset's bytes back. +fn check_h5py(path: &Path, e: &Expect, dir: &Path) { + let mut script = format!( + "import h5py, numpy as np\nf = h5py.File({p:?}, 'r')\n", + p = path.to_str().unwrap() + ); + for (i, (name, want)) in e.datasets.iter().enumerate() { + let exp = dir.join(format!("expect{i}.bin")); + std::fs::write(&exp, want).unwrap(); + script.push_str(&format!( + "a = np.ascontiguousarray(f[{name:?}][()]).tobytes()\n\ + assert a == open({x:?}, 'rb').read(), {name:?}\n", + x = exp.to_str().unwrap() + )); + } + script.push_str( + "assert f.attrs['text'] in (b'readable by 1.8', 'readable by 1.8')\n\ + assert list(f.attrs['ints']) == [1, -2, 3]\n\ + assert len(f['g_dense'].attrs) == 12 and len(f['g_dense']) == 20\n\ + assert list(f['g_order']) == ['z9', 'z8', 'z7', 'z6', 'z5', 'z4', 'z3', 'z2', 'z1', 'z0']\n\ + assert f['g_compact/soft'].shape == (1000,)\n\ + assert f['resizable'].maxshape == (None,)\n\ + assert f['filled'].fillvalue == -9\n\ + print('ok')\n", + ); + assert_eq!(py(&script), "ok"); +} + +/// h5dump (1.14) and `h5rs check --data` accept the file. +fn check_tools(path: &Path) { + let p = path.to_str().unwrap(); + let o = Command::new(env!("CARGO_BIN_EXE_h5rs")) + .args(["check", "--data", "-q", p]) + .output() + .unwrap(); + assert!(o.status.success(), "h5rs check --data {p}:\n{}", text(&o)); + let o = Command::new("h5dump") + .args(["-o", "/dev/null", p]) + .output() + .unwrap(); + assert!( + o.status.success() && o.stderr.is_empty(), + "h5dump {p}:\n{}", + text(&o) + ); +} + +/// HDF5 1.8's h5dump reads the whole file exactly as h5dump 1.14 does, and +/// dumps each numeric dataset's values as the bytes we wrote. +fn check_18(h5dump: &Path, path: &Path, e: &Expect, dir: &Path) { + let p = path.to_str().unwrap(); + let o = Command::new(h5dump).arg(p).output().unwrap(); + assert!( + o.status.success() && o.stderr.is_empty(), + "h5dump 1.8 {p}:\n{}", + text(&o) + ); + let out18 = String::from_utf8_lossy(&o.stdout).to_string(); + assert!(out18.contains("EXTERNAL_LINK \"ext\""), "{out18}"); + let o = Command::new("h5dump").arg(p).output().unwrap(); + assert!(o.status.success(), "h5dump {p}:\n{}", text(&o)); + // HDF5 1.8 has no name for IEEE half floats. + let out = String::from_utf8_lossy(&o.stdout).replace( + "H5T_IEEE_F16LE", + "16-bit little-endian floating-point 16-bit precision", + ); + if let Some((n, (a, b))) = out18 + .lines() + .zip(out.lines()) + .enumerate() + .find(|(_, (a, b))| a != b) + { + panic!("h5dump 1.8 and h5dump differ at line {}:\n{a}\n{b}", n + 1); + } + assert_eq!(out18.lines().count(), out.lines().count()); + // `-b` writes nothing for compounds, enums, arrays, strings and half + // floats (HDF5 1.8 + // has no native half float), for h5py's files too: those are covered + // by the text above. + let skip = ["/f16", "/strings", "/compound", "/enum", "/array"]; + for (i, (name, want)) in e.datasets.iter().enumerate() { + if want.is_empty() || skip.contains(&name.as_str()) { + continue; + } + let bin = dir.join(format!("dump18_{i}.bin")); + let o = Command::new(h5dump) + .args(["-d", name, "-b", "LE", "-o"]) + .arg(&bin) + .arg(p) + .output() + .unwrap(); + assert!(o.status.success(), "h5dump 1.8 -d {name}:\n{}", text(&o)); + let got = std::fs::read(&bin).unwrap(); + assert!(got == *want, "HDF5 1.8 read of {name} differs"); + } +} + +fn check_all(path: &Path, e: &Expect, dir: &Path, h5dump: Option<&Path>) { + check_ours(path, e); + check_h5py(path, e, dir); + check_tools(path); + if let Some(h) = h5dump { + check_18(h, path, e, dir); + } +} + +#[test] +fn every_feature_reads_in_hdf5_1_8() { + if !tools_ok() { + return; + } + let h18 = h5dump18(); + let dir = tmpdir(); + let path = dir.path().join("v18.h5"); + let mut e = write_file(&path); + check_versions(&path); + check_all(&path, &e, dir.path(), h18.as_deref()); + + // FileEditor: grow the unlimited datasets (splitting B-tree nodes), + // overwrite values in filtered and unfiltered chunks, set attributes + // (compact and dense). + let mut ed = FileEditor::open(&path).unwrap(); + let mut res: Vec = (0..25).collect(); + for round in 0..300usize { + let n = res.len() as u64; + let add = 1 + (round % 9) as u64; + ed.resize("resizable", &[n + add]).unwrap(); + let vals: Vec = (0..add).map(|k| (n + k) as i32 * 7 - 3).collect(); + ed.write_values("resizable", &block(&[n], &[add]), &vals) + .unwrap(); + res.extend(&vals); + } + e.set("/resizable", le(&res, i32::to_le_bytes)); + let mut many: Vec = (0..MANY as i32).map(|i| i ^ 0x5a5a).collect(); + ed.resize("many", &[MANY as u64 + 500]).unwrap(); + let vals: Vec = (0..500).collect(); + ed.write_values("many", &block(&[MANY as u64], &[500]), &vals) + .unwrap(); + many.extend(&vals); + many[12_345] = -1; + ed.write_values("many", &block(&[12_345], &[1]), &[-1i32]) + .unwrap(); + e.set("/many", le(&many, i32::to_le_bytes)); + let mut chunked: Vec = (0..60 * 70).map(|i| (i % 97) as f32 * 1.25).collect(); + let row: Vec = (0..70).map(|i| -(i as f32)).collect(); + ed.write_values("chunked", &block(&[31, 0], &[1, 70]), &row) + .unwrap(); + chunked[31 * 70..32 * 70].copy_from_slice(&row); + e.set("/chunked", le(&chunked, f32::to_le_bytes)); + ed.resize("resizable2", &[9, 6]).unwrap(); + let mut r2: Vec = (0..30).map(|i| i * 1_000_000_007).collect(); + r2.extend(std::iter::repeat_n(0, 24)); + e.set("/resizable2", le(&r2, i64::to_le_bytes)); + ed.set_attr("contig", "added", &AttrValue::F64(3.25)) + .unwrap(); + ed.set_attr("g_dense", "attr03", &AttrValue::F64(-1.0)) + .unwrap(); + ed.set_attr("/", "text", &AttrValue::String("edited".into())) + .unwrap(); + drop(ed); + + check_ours(&path, &e); + let f = File::open(&path).unwrap(); + assert!(matches!( + f.dataset("contig").unwrap().attr("added").unwrap(), + Some(AttrValue::F64(v)) if v == 3.25 + )); + drop(f); + // Everything but the root's "text" attribute is as check_h5py expects. + py(&format!( + "import h5py\nf = h5py.File({p:?}, 'r')\n\ + assert f.attrs['text'] in (b'edited', 'edited')\n\ + assert f['contig'].attrs['added'] == 3.25\n\ + assert f['g_dense'].attrs['attr03'] == -1.0\n", + p = path.to_str().unwrap() + )); + let o = Command::new(env!("CARGO_BIN_EXE_h5rs")) + .args(["check", "--data", "-q", path.to_str().unwrap()]) + .output() + .unwrap(); + assert!(o.status.success(), "h5rs check after edits:\n{}", text(&o)); + if let Some(h) = h18.as_deref() { + check_18(h, &path, &e, dir.path()); + } + + // h5py (libhdf5 2.x) goes on appending to the 1.8-format file, and 1.8 + // still reads it. + py(&format!( + "import h5py, numpy as np\n\ + with h5py.File({p:?}, 'r+') as f:\n\ + \x20 d = f['resizable']\n\ + \x20 n = d.shape[0]\n\ + \x20 d.resize((n + 100,))\n\ + \x20 d[n:] = np.arange(100, dtype=')>, + children: Vec, +} + +fn read_tree(bytes: &[u8], addr: u64, ndims: usize, sizes: bool) -> TreeNode { + let a = addr as usize; + assert_eq!(&bytes[a..a + 4], b"TREE"); + let level = bytes[a + 5]; + let n = u16::from_le_bytes([bytes[a + 6], bytes[a + 7]]) as usize; + let mut p = a + 24; + let mut keys = Vec::new(); + let mut kids = Vec::new(); + for i in 0..=n { + let size = u32::from_le_bytes(bytes[p..p + 4].try_into().unwrap()); + let offs = (0..ndims) + .map(|d| u64::from_le_bytes(bytes[p + 8 + 8 * d..p + 16 + 8 * d].try_into().unwrap())) + .collect(); + keys.push((if sizes { size } else { 0 }, offs)); + p += 8 + 8 * ndims; + if i < n { + kids.push(u64::from_le_bytes(bytes[p..p + 8].try_into().unwrap())); + p += 8; + } + } + let children = if level > 0 { + kids.iter() + .map(|&c| read_tree(bytes, c, ndims, sizes)) + .collect() + } else { + Vec::new() + }; + TreeNode { + level, + keys, + children, + } +} + +/// Where two trees first differ (libhdf5's first). +fn first_difference(a: &TreeNode, b: &TreeNode, at: &str) -> Option { + if a.level != b.level || a.keys.len() != b.keys.len() { + return Some(format!( + "{at}: level {} with {} keys vs level {} with {} keys", + a.level, + a.keys.len(), + b.level, + b.keys.len() + )); + } + if let Some(i) = (0..a.keys.len()).find(|&i| a.keys[i] != b.keys[i]) { + return Some(format!("{at} key {i}: {:?} vs {:?}", a.keys[i], b.keys[i])); + } + a.children + .iter() + .zip(&b.children) + .enumerate() + .find_map(|(i, (x, y))| first_difference(x, y, &format!("{at}/{i}"))) +} + +/// The chunk B-tree of dataset `name`: (tree, ndims) from its layout. +fn tree_of(path: &Path, name: &str, sizes: bool) -> TreeNode { + use clawhdf5_format::message_type::MessageType; + use clawhdf5_format::object_header::ObjectHeader; + let bytes = std::fs::read(path).unwrap(); + let f = File::open(path).unwrap(); + let sb = f.superblock().clone(); + let addr = clawhdf5_format::group_v2::resolve_path_any(&bytes, &sb, name).unwrap(); + let oh = ObjectHeader::parse(&bytes, addr as usize, 8, 8).unwrap(); + let l = &oh + .messages + .iter() + .find(|m| m.msg_type == MessageType::DataLayout) + .unwrap() + .data; + assert_eq!((l[0], l[1]), (3, 2), "{name}: layout v3, chunked"); + let ndims = l[2] as usize; + let root = u64::from_le_bytes(l[3..11].try_into().unwrap()); + read_tree(&bytes, root, ndims, sizes) +} + +/// Our version-1 chunk B-trees are libhdf5's, node for node (levels, child +/// counts, every key's offsets, and its chunk size where the chunks are +/// the same bytes), for 1-D, 2-D and 3-D datasets with two- and three-level +/// trees, filtered or not. +/// +/// libhdf5 inserts each chunk into the tree when it leaves its chunk cache. +/// A whole-dataset write with no cache (`rdcc_nbytes=0`, or chunks larger +/// than the cache) inserts them in row-major order, as we build the tree; +/// with the default cache small chunks of a 1-D dataset still arrive in +/// order, but those of a multi-dimensional one arrive in the order the +/// cache's hash evicts them, which gives the same keys in differently +/// filled nodes. Both trees index the same chunks; we do not model the +/// cache. +#[test] +fn chunk_btrees_match_libhdf5() { + if !tools_ok() { + return; + } + let dir = tmpdir(); + let theirs = dir.path().join("libhdf5.h5"); + py(&format!( + "import h5py, numpy as np\n\ + with h5py.File({p:?}, 'w', libver=('v108', 'latest'), rdcc_nbytes=0) as f:\n\ + \x20 f.create_dataset('d1000', data=np.arange(10000.0), chunks=(10,))\n\ + \x20 f.create_dataset('d999', data=np.arange(9990.0), chunks=(10,))\n\ + \x20 f.create_dataset('big', data=np.arange(100000, dtype=' = (0..10000).map(f64::from).collect(); + b.create_dataset("d1000") + .with_f64_data(&v) + .with_chunks(&[10]); + b.create_dataset("d999") + .with_f64_data(&v[..9990]) + .with_chunks(&[10]); + let v: Vec = (0..100_000).collect(); + b.create_dataset("big") + .with_i32_data(&v) + .with_chunks(&[1]) + .with_maxshape(&[u64::MAX]); + let v: Vec = (0..10000).map(f64::from).collect(); + b.create_dataset("grid") + .with_f64_data(&v) + .with_shape(&[100, 100]) + .with_chunks(&[1, 1]); + let raw: Vec = (0..27000i16).flat_map(|x| x.to_le_bytes()).collect(); + b.create_dataset("cube") + .with_compound_data( + clawhdf5_format::datatype::Datatype::FixedPoint { + size: 2, + byte_order: clawhdf5_format::datatype::DatatypeByteOrder::LittleEndian, + signed: true, + bit_offset: 0, + bit_precision: 16, + }, + raw, + 27000, + ) + .with_shape(&[30, 30, 30]) + .with_chunks(&[2, 3, 5]); + let v: Vec = (0..20000).map(|i| i % 13).collect(); + b.create_dataset("gz") + .with_i64_data(&v) + .with_chunks(&[7]) + .with_deflate(4); + b.write(&ours).unwrap(); + for (name, sizes) in [ + ("d1000", true), + ("d999", true), + ("big", true), + ("grid", true), + ("cube", true), + // Compressed sizes differ between zlib-rs and zlib. + ("gz", false), + ] { + let a = tree_of(&theirs, name, sizes); + let b = tree_of(&ours, name, sizes); + if let Some(d) = first_difference(&a, &b, "root") { + panic!("{name}: our chunk B-tree differs from libhdf5's: {d}"); + } + } +} diff --git a/crates/clawhdf5/src/edit/earray.rs b/crates/clawhdf5/src/edit/earray.rs index bae196f..6c777dd 100644 --- a/crates/clawhdf5/src/edit/earray.rs +++ b/crates/clawhdf5/src/edit/earray.rs @@ -49,17 +49,9 @@ pub(crate) fn encode_elem( /// element for chunks of `chunk_bytes` bytes (`H5D__earray_idx_create`, /// `H5D__farray_idx_create`): one byte more than the nominal size needs — /// except under layout message version 5 (HDF5 2.0's own format), which -/// always uses 8 bytes. -pub(crate) fn chunk_size_len(chunk_bytes: u64, layout_version: u8) -> usize { - if layout_version >= 5 { - return 8; - } - let log2 = if chunk_bytes <= 1 { - 0 - } else { - 63 - chunk_bytes.leading_zeros() - }; - (1 + ((log2 + 8) / 8) as usize).min(8) +/// uses the file's size of lengths (`length_size`). +pub(crate) fn chunk_size_len(chunk_bytes: u64, layout_version: u8, length_size: u8) -> usize { + clawhdf5_format::chunked_write::chunk_size_len(chunk_bytes, layout_version, length_size) } /// Creation parameters, in the layout message's order. @@ -226,7 +218,7 @@ impl Ea { ) -> Result { let os = img.os; let elem_size = if filtered { - os as usize + chunk_size_len(chunk_bytes, layout_version) + 4 + os as usize + chunk_size_len(chunk_bytes, layout_version, img.ls) + 4 } else { os as usize }; diff --git a/crates/clawhdf5/src/edit/farray.rs b/crates/clawhdf5/src/edit/farray.rs index e7a1696..8e737fc 100644 --- a/crates/clawhdf5/src/edit/farray.rs +++ b/crates/clawhdf5/src/edit/farray.rs @@ -83,7 +83,7 @@ impl Fa { let os = img.os; let osz = os as usize; let elem_size = if filtered { - osz + chunk_size_len(chunk_bytes, layout_version) + 4 + osz + chunk_size_len(chunk_bytes, layout_version, img.ls) + 4 } else { osz }; diff --git a/crates/clawhdf5/src/edit/mod.rs b/crates/clawhdf5/src/edit/mod.rs index 28e1fa6..169c6d3 100644 --- a/crates/clawhdf5/src/edit/mod.rs +++ b/crates/clawhdf5/src/edit/mod.rs @@ -412,12 +412,24 @@ impl<'t> ChunkedEdit<'t> { .iter() .map(|&c| u64::from(c)) .collect(); + // Chunks of 4 GiB or more (HDF5 2.0 writes them with layout message + // version 5) are read, but not rewritten: the editor holds a chunk + // it rewrites in memory, compressed and not. Writing values, and a + // resize that prunes or allocates chunks, are refused here, before + // anything is written; growing the extent (late allocation) and + // setting attributes do not come here. let chunk_bytes = cd .iter() .try_fold(t.es as u64, |a, &c| a.checked_mul(c)) + .filter(|&b| b <= u64::from(u32::MAX)) .and_then(|b| usize::try_from(b).ok()) - .filter(|&b| b <= u32::MAX as usize) - .ok_or_else(|| Error::Unsupported("chunk larger than 4 GiB".into()))?; + .ok_or_else(|| { + Error::Unsupported( + "rewriting chunks of 4 GiB or more (writing values, or a resize that \ + prunes or allocates chunks)" + .into(), + ) + })?; let hdr = Header::load(img, t.addr)?; let layout_msg = hdr.find(MSG_LAYOUT).ok_or(Error::MissingMessage( clawhdf5_format::message_type::MessageType::DataLayout, @@ -600,7 +612,8 @@ impl<'t> ChunkedEdit<'t> { let d = self.hdr.data(img, self.layout_msg)?; let node_size = u32::from_le_bytes([d[at], d[at + 1], d[at + 2], d[at + 3]]); let size_len = if filtered { - earray::chunk_size_len(self.chunk_bytes as u64, self.lpos.version) + 4 + earray::chunk_size_len(self.chunk_bytes as u64, self.lpos.version, img.ls) + + 4 } else { 0 }; @@ -671,7 +684,7 @@ impl<'t> ChunkedEdit<'t> { } _ => return Err(Error::Unsupported("chunk index missing".into())), } - img.free(info.address, u64::from(info.chunk_size)); + img.free(info.address, info.chunk_size); Ok(()) } @@ -1625,7 +1638,7 @@ fn store_chunk( let len = bytes.len() as u64; let placed = match existing { Some(info) if t.pipeline.is_none() => { - if u64::from(info.chunk_size) != len { + if info.chunk_size != len { return Err(Error::Unsupported( "unfiltered chunk stored at an unexpected size".into(), )); @@ -1637,11 +1650,10 @@ fn store_chunk( // thing in the file (the chunk an append keeps rewriting usually is) // and can grow there. Some(info) - if len <= u64::from(info.chunk_size) - || img.grow_tail(info.address, u64::from(info.chunk_size), len)? => + if len <= info.chunk_size || img.grow_tail(info.address, info.chunk_size, len)? => { img.write(info.address, &bytes)?; - (len != u64::from(info.chunk_size) || info.filter_mask != mask).then_some(Elem { + (len != info.chunk_size || info.filter_mask != mask).then_some(Elem { addr: info.address, size: len, mask, @@ -1651,7 +1663,7 @@ fn store_chunk( let a = img.alloc(len)?; img.write(a, &bytes)?; if let Some(info) = existing { - img.free(info.address, u64::from(info.chunk_size)); + img.free(info.address, info.chunk_size); } Some(Elem { addr: a, @@ -1669,14 +1681,16 @@ fn store_chunk( fn img_read<'a>(f: &'a File, info: &ChunkInfo) -> Result<&'a [u8], Error> { let start = usize::try_from(info.address) .map_err(|_| Error::Unsupported("chunk address out of range".into()))?; - f.as_bytes() - .get(start..start + info.chunk_size as usize) - .ok_or_else(|| { - Error::Format(clawhdf5_format::error::FormatError::UnexpectedEof { - expected: start + info.chunk_size as usize, - available: f.as_bytes().len(), - }) + let end = usize::try_from(info.chunk_size) + .ok() + .and_then(|n| start.checked_add(n)) + .unwrap_or(usize::MAX); + f.as_bytes().get(start..end).ok_or_else(|| { + Error::Format(clawhdf5_format::error::FormatError::UnexpectedEof { + expected: end, + available: f.as_bytes().len(), }) + }) } fn decode_chunk( diff --git a/crates/clawhdf5/src/lib.rs b/crates/clawhdf5/src/lib.rs index 24bd435..9678c29 100644 --- a/crates/clawhdf5/src/lib.rs +++ b/crates/clawhdf5/src/lib.rs @@ -68,6 +68,7 @@ pub use writer::{DatasetSpec, create_datasets_parallel}; // Re-export useful types from clawhdf5-format for advanced users pub use clawhdf5_format::data_layout::VdsMapping; pub use clawhdf5_format::dict_encoding::{DictEncoded, DictionaryEncoder}; +pub use clawhdf5_format::libver::LibVer; pub use clawhdf5_format::property_list::{ DatasetCreateProps, FileAccessProps, FileCreateProps, lib_version, }; diff --git a/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index 4456d66..4426d39 100644 --- a/crates/clawhdf5/src/reader.rs +++ b/crates/clawhdf5/src/reader.rs @@ -1528,11 +1528,11 @@ impl<'f> Dataset<'f> { let ds = self.dataspace()?; let dl = self.data_layout()?; let pipeline = self.filter_pipeline()?; - // The selection reader knows nothing about fill values. When they - // matter — no storage at all, or a non-zero fill on a chunked (possibly - // sparse) dataset — select from a fill-aware full read instead. (The - // selection reader currently decodes the full dataset too, so this - // costs nothing extra.) + // Fill values matter when there is no storage at all, or a non-zero + // fill on a chunked (possibly sparse) dataset: a chunked dataset's + // selection is then read over a box of fill values; anything else + // (and a selection whose box is most of the dataset) is selected + // from a fill-aware full read. let fill = clawhdf5_format::fill_value::dataset_fill_value_from_storage( &self.file.data, &self.header.messages, @@ -1545,6 +1545,29 @@ impl<'f> Dataset<'f> { && !clawhdf5_format::fill_value::is_default(fill.as_deref())); if fill_matters { clawhdf5_format::partial_read::validate(selection, &ds.dimensions)?; + // A chunked dataset: read the chunks the selection touches over + // a box of fill values, not the whole dataset (which may be many + // GiB when its chunks are large). + if matches!(dl, DataLayout::Chunked { .. }) { + clawhdf5_format::chunked_read::check_chunk_element_size( + &dl, + &dt, + self.file.offset_size(), + )?; + if let Some(selected) = clawhdf5_format::partial_read::read_selection_filled_in( + &self.file.data, + &dl, + &ds, + dt.type_size() as usize, + pipeline.as_ref(), + self.file.offset_size(), + self.file.length_size(), + selection, + Some(fill.as_deref().unwrap_or(&[])), + )? { + return Ok(selected); + } + } let full = self.read_raw()?; return Ok(data_read::extract_selection_from_buffer( &full, diff --git a/crates/clawhdf5/src/writer.rs b/crates/clawhdf5/src/writer.rs index ea48bc9..9fcf6fc 100644 --- a/crates/clawhdf5/src/writer.rs +++ b/crates/clawhdf5/src/writer.rs @@ -1,6 +1,7 @@ //! Writing API: FileBuilder and GroupBuilder for creating HDF5 files. use clawhdf5_format::file_writer::FileWriter as FormatWriter; +use clawhdf5_format::libver::LibVer; use clawhdf5_format::type_builders::{ AttrValue, DatasetBuilder as FormatDatasetBuilder, FinishedGroup, GroupBuilder as FormatGroupBuilder, @@ -98,6 +99,34 @@ impl FileBuilder { self } + /// Set the library version bounds, as h5py's `libver=(low, high)`: the + /// oldest HDF5 release whose file format is used, and the newest whose + /// features are allowed. The default is `(LibVer::V110, + /// LibVer::Latest)`, the HDF5 1.10 format clawhdf5 has always written. + /// + /// `libver_bounds(LibVer::V18, LibVer::V18)` writes a file HDF5 1.8 can + /// read (version-1 B-tree chunk indexes, version-2 superblock, the low + /// bound libhdf5 2.0 uses by default), or fails with + /// [`Error::Format`] (`FormatError::LibverBound`) for anything 1.8 + /// cannot read: virtual datasets, the 1.12 reference types, native + /// complex numbers. See `clawhdf5_format::libver` for the details. + /// + /// ``` + /// use clawhdf5::{FileBuilder, LibVer}; + /// + /// let mut b = FileBuilder::new(); + /// b.libver_bounds(LibVer::V18, LibVer::V18); + /// b.create_dataset("x") + /// .with_f64_data(&[1.0, 2.0, 3.0]) + /// .with_maxshape(&[u64::MAX]); + /// let bytes = b.finish().unwrap(); + /// assert_eq!(bytes[8], 2); // superblock version 2 + /// ``` + pub fn libver_bounds(&mut self, low: LibVer, high: LibVer) -> &mut Self { + self.writer.libver_bounds(low, high); + self + } + /// Set an attribute on the root group. pub fn set_attr(&mut self, name: &str, value: AttrValue) { self.writer.set_root_attr(name, value); diff --git a/crates/clawhdf5/tests/fixtures/gen_huge_chunks.py b/crates/clawhdf5/tests/fixtures/gen_huge_chunks.py new file mode 100644 index 0000000..f901fd5 --- /dev/null +++ b/crates/clawhdf5/tests/fixtures/gen_huge_chunks.py @@ -0,0 +1,110 @@ +"""Write HDF5 files whose chunks are 4 GiB or more, one dataset per chunk index. + + python gen_huge_chunks.py filtered OUT.h5 # huge_chunks_filtered.h5 + python gen_huge_chunks.py unfiltered OUT.h5 # generated at test time + +libhdf5 2.x writes a chunk of more than 0xFFFFFFFF bytes with layout message +version 5 (`H5D__chunk_construct`: "chunk size > 4GB requires +H5F_LIBVER_V200"), so never with a version-1 B-tree; a filtered chunk index +element of a version-5 layout stores the chunk's size in "size of lengths" +bytes (8). Every dataset is ` N: + ds[N : N + 10] = np.arange(100.0, 110.0) + if index == "implicit": + ds[2 * N - 10 : 2 * N] = np.arange(500.0, 510.0) + dsid.close() + + +def main(): + mode, out = sys.argv[1], sys.argv[2] + filtered = {"filtered": True, "unfiltered": False}[mode] + fapl = h5py.h5p.create(h5py.h5p.FILE_ACCESS) + fapl.set_libver_bounds(h5py.h5f.LIBVER_EARLIEST, h5py.h5f.LIBVER_V200) + fid = h5py.h5f.create(out.encode(), h5py.h5f.ACC_TRUNC, fapl=fapl) + indexes = ["single", "farray", "earray", "btree2"] + if not filtered: + indexes.insert(1, "implicit") + for index in indexes: + make(fid, index, index, filtered) + fid.close() + + +if __name__ == "__main__": + main() diff --git a/crates/clawhdf5/tests/fixtures/huge_chunks_filtered.h5 b/crates/clawhdf5/tests/fixtures/huge_chunks_filtered.h5 new file mode 100644 index 0000000..60a8931 Binary files /dev/null and b/crates/clawhdf5/tests/fixtures/huge_chunks_filtered.h5 differ diff --git a/crates/clawhdf5/tests/huge_chunks_interop.rs b/crates/clawhdf5/tests/huge_chunks_interop.rs new file mode 100644 index 0000000..52885e7 --- /dev/null +++ b/crates/clawhdf5/tests/huge_chunks_interop.rs @@ -0,0 +1,518 @@ +//! Chunks of 4 GiB or more, which HDF5 2.0 writes (layout message version 5, +//! `H5F_LIBVER_V200`), in every chunk index libhdf5 uses for them. +//! +//! `fixtures/huge_chunks_filtered.h5` (written by libhdf5 2.0.0 through h5py +//! 3.16, `fixtures/gen_huge_chunks.py filtered`) holds one dataset per +//! filtered index — Single Chunk, Fixed Array, Extensible Array, v2 B-tree — +//! whose chunks are 2^29 + 1 `f64` (4 GiB + 8 bytes), stored deflated twice +//! so each takes about 20 KiB. Listing their chunks is cheap and always runs; +//! decoding one inflates 4 GiB, so those tests are opt-in: +//! `CLAWHDF5_HUGE_CHUNKS=1` (run them one at a time, `--test-threads=1`: +//! each needs about 4.5 GiB of memory). The same variable enables the +//! unfiltered tests, which have h5py write a sparse file (about 44 GiB long, +//! a few blocks on disk) under `tests/scratch/`, and the writer tests. +//! +//! libhdf5 never writes such a chunk with a version-1 B-tree (a chunk of more +//! than 0xFFFFFFFF bytes forces layout version 5 and with it the newer +//! indexes) and refuses to open one; `chunked_read`'s unit tests cover that. + +use std::path::{Path, PathBuf}; +use std::process::Command; + +use clawhdf5::File; +use clawhdf5_format::chunked_read::list_chunks; +use clawhdf5_format::data_layout::DataLayout; +use clawhdf5_format::dataspace::Dataspace; +use clawhdf5_format::message_type::MessageType; +use clawhdf5_format::object_header::ObjectHeader; +use clawhdf5_format::selection::Selection; +use clawhdf5_format::superblock::Superblock; + +/// Elements per chunk along the chunked axis: 4 GiB + 8 bytes of `f64`. +const N: u64 = (1 << 29) + 1; + +fn fixture() -> PathBuf { + Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures/huge_chunks_filtered.h5") +} + +fn heavy() -> bool { + if std::env::var("CLAWHDF5_HUGE_CHUNKS").is_ok_and(|v| v == "1") { + return true; + } + eprintln!("SKIP: set CLAWHDF5_HUGE_CHUNKS=1 to decode 4 GiB chunks"); + false +} + +fn python() -> String { + std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string()) +} + +fn python_available() -> bool { + Command::new(python()) + .args(["-c", "import h5py"]) + .output() + .is_ok_and(|o| o.status.success()) +} + +fn interop_required() -> bool { + std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1") +} + +/// The raw layout message version, the parsed layout, the dataspace and the +/// chunks of `name` in the in-memory file `data`. +fn layout_of( + data: &[u8], + name: &str, +) -> ( + u8, + DataLayout, + Dataspace, + Vec, +) { + 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 = |t| { + hdr.messages + .iter() + .find(|m| m.msg_type == t) + .unwrap_or_else(|| panic!("{name}: no {t:?} message")) + }; + let lm = msg(MessageType::DataLayout); + let layout = DataLayout::parse(&lm.data, sb.offset_size, sb.length_size).unwrap(); + let space = Dataspace::parse(&msg(MessageType::Dataspace).data, sb.length_size).unwrap(); + let (chunks, _) = + list_chunks(data, &layout, &space, 8, sb.offset_size, sb.length_size).unwrap(); + (lm.data[0], layout, space, chunks) +} + +/// The fixture's datasets: name, chunk index type, shape, and the scaled +/// origins of the chunks libhdf5 wrote. +type Case = (&'static str, u8, &'static [u64], &'static [&'static [u64]]); + +const FILTERED: &[Case] = &[ + ("single", 1, &[N], &[&[0]]), + ("farray", 3, &[N + 10], &[&[0], &[N]]), + ("earray", 4, &[N + 10], &[&[0], &[N]]), + ( + "btree2", + 5, + &[2, N + 10], + &[&[0, 0], &[0, N], &[1, 0], &[1, N]], + ), +]; + +/// Every index of the fixture: layout version 5, its chunks listed at the +/// right origins with their stored (deflated) sizes. The index elements +/// store a size in 8 bytes (libhdf5's `H5F_SIZEOF_SIZE` under layout +/// version 5), where version 4 would use 6 for a chunk this size. +#[test] +fn filtered_huge_chunk_indexes_list() { + let data = std::fs::read(fixture()).unwrap(); + for &(name, index, shape, origins) in FILTERED { + let (version, layout, space, mut chunks) = layout_of(&data, name); + assert_eq!(version, 5, "{name}: layout message version"); + let DataLayout::Chunked { + chunk_dimensions, + chunk_index_type, + .. + } = &layout + else { + panic!("{name}: not chunked: {layout:?}"); + }; + assert_eq!(*chunk_index_type, Some(index), "{name}"); + assert_eq!(chunk_dimensions.last(), Some(&8), "{name}"); + assert_eq!(space.dimensions, shape, "{name}"); + chunks.sort_by(|a, b| a.offsets.cmp(&b.offsets)); + let got: Vec<&[u64]> = chunks.iter().map(|c| &c.offsets[..shape.len()]).collect(); + assert_eq!(got, origins, "{name}"); + for c in &chunks { + assert!( + (10_000..40_000).contains(&c.chunk_size), + "{name}: stored size {} of chunk {:?}", + c.chunk_size, + c.offsets + ); + assert_eq!(c.filter_mask, 0); + } + } +} + +/// `f64` values of `sel` in dataset `name` of `file`. +fn sel(file: &File, name: &str, start: &[u64], count: &[u64]) -> Vec { + let ds = file.dataset(name).unwrap(); + let s = Selection::Hyperslab { + start: start.to_vec(), + stride: vec![1; start.len()], + count: count.to_vec(), + block: vec![1; start.len()], + }; + ds.read_f64_selection(&s).unwrap() +} + +fn f(range: std::ops::Range) -> Vec { + range.map(f64::from).collect() +} + +/// Reads that touch a few elements of each 4 GiB chunk: written values, the +/// fill value (-1) next to them, and the edge of the dataset. +#[test] +fn filtered_huge_chunks_read() { + if !heavy() { + return; + } + let file = File::open(fixture()).unwrap(); + let mut first = f(0..10); + first.extend([-1.0; 2]); + for name in ["single", "farray", "earray"] { + assert_eq!(sel(&file, name, &[0], &[12]), first, "{name}"); + } + assert_eq!(sel(&file, "single", &[N - 2], &[2]), [-1.0; 2]); + let mut edge = vec![-1.0; 2]; + edge.extend(f(100..110)); + for name in ["farray", "earray"] { + assert_eq!(sel(&file, name, &[N - 2], &[12]), edge, "{name}"); + } + assert_eq!(sel(&file, "btree2", &[0, 0], &[1, 10]), f(0..10)); + assert_eq!(sel(&file, "btree2", &[1, N], &[1, 10]), f(300..310)); + assert_eq!( + sel(&file, "btree2", &[0, N - 1], &[2, 2]), + [-1.0, 100.0, -1.0, 300.0] + ); +} + +/// A directory under `tests/scratch/` (on disk: the sparse files must not +/// land on a tmpfs `/tmp`), removed when dropped. +fn scratch() -> tempfile::TempDir { + let root = Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/scratch"); + std::fs::create_dir_all(&root).unwrap(); + tempfile::tempdir_in(root).unwrap() +} + +/// Have h5py (libhdf5 2.x) write `fixtures/gen_huge_chunks.py`'s file for +/// `mode` into `dir`; `None` when there is no h5py (a failure under +/// `CLAWHDF5_REQUIRE_INTEROP=1`). +fn generate(dir: &Path, mode: &str) -> Option { + 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 None; + } + let script = Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures/gen_huge_chunks.py"); + let path = dir.join(format!("{mode}.h5")); + let out = Command::new(python()) + .arg(&script) + .arg(mode) + .arg(&path) + .output() + .unwrap(); + assert!( + out.status.success(), + "gen_huge_chunks.py {mode} failed:\n{}", + String::from_utf8_lossy(&out.stderr) + ); + Some(path) +} + +/// Positioned reads of a file (no mmap), counting the bytes read. +struct Counting { + file: std::fs::File, + len: u64, + read: std::sync::atomic::AtomicU64, +} + +impl clawhdf5_format::storage::Storage for Counting { + fn read_at( + &self, + offset: u64, + len: usize, + ) -> Result, clawhdf5_format::error::FormatError> { + use std::os::unix::fs::FileExt; + let len = len.min(usize::try_from(self.len.saturating_sub(offset)).unwrap_or(usize::MAX)); + let mut buf = vec![0; len]; + self.file.read_exact_at(&mut buf, offset).unwrap(); + self.read + .fetch_add(len as u64, std::sync::atomic::Ordering::Relaxed); + Ok(buf.into()) + } + + fn len(&self) -> u64 { + self.len + } +} + +fn check_unfiltered(file: &File) { + let mut first = f(0..10); + first.extend([0.0; 2]); + for name in ["single", "implicit", "farray", "earray"] { + assert_eq!(sel(file, name, &[0], &[12]), first, "{name}"); + } + assert_eq!(sel(file, "single", &[N - 2], &[2]), [0.0; 2]); + let mut edge = vec![0.0; 2]; + edge.extend(f(100..110)); + for name in ["implicit", "farray", "earray"] { + assert_eq!(sel(file, name, &[N - 2], &[12]), edge, "{name}"); + } + let mut last = vec![0.0; 2]; + last.extend(f(500..510)); + assert_eq!(sel(file, "implicit", &[2 * N - 12], &[12]), last); + assert_eq!(sel(file, "btree2", &[0, 0], &[1, 10]), f(0..10)); + assert_eq!(sel(file, "btree2", &[1, N], &[1, 10]), f(300..310)); + assert_eq!( + sel(file, "btree2", &[0, N - 1], &[2, 2]), + [0.0, 100.0, 0.0, 300.0] + ); +} + +/// Unfiltered chunks of 4 GiB + 8 bytes in every index libhdf5 gives them +/// (Single Chunk, Implicit, Fixed Array, Extensible Array, v2 B-tree), in a +/// sparse file h5py writes: read through a memory map, and through +/// positioned reads, where a selection reads only the rows it needs. +#[test] +fn unfiltered_huge_chunks_read() { + if !heavy() { + return; + } + let dir = scratch(); + let Some(path) = generate(dir.path(), "unfiltered") else { + return; + }; + check_unfiltered(&File::open(&path).unwrap()); + + let file = std::fs::File::open(&path).unwrap(); + let len = file.metadata().unwrap().len(); + assert!(len > 40 << 30, "{len}"); + let storage = std::sync::Arc::new(Counting { + file, + len, + read: 0.into(), + }); + let positioned = File::open_storage(storage.clone()).unwrap(); + check_unfiltered(&positioned); + // Every read above together: metadata and a few rows, not 4 GiB chunks. + let read = storage.read.load(std::sync::atomic::Ordering::Relaxed); + assert!(read < 1 << 20, "{read} bytes read"); +} + +/// `FileEditor` does not rewrite chunks of 4 GiB or more: writing values, +/// or a resize that prunes or allocates chunks, is refused before anything +/// is written. Growing the extent and setting attributes still work. +#[test] +fn editor_refuses_rewriting_huge_chunks() { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("huge.h5"); + std::fs::copy(fixture(), &path).unwrap(); + let before = std::fs::read(&path).unwrap(); + let mut ed = clawhdf5::FileEditor::open(&path).unwrap(); + let one = Selection::Hyperslab { + start: vec![0], + stride: vec![1], + count: vec![1], + block: vec![1], + }; + for name in ["single", "farray", "earray"] { + let err = ed.write_values(name, &one, &[5.0f64]).unwrap_err(); + assert!( + matches!(&err, clawhdf5::Error::Unsupported(m) if m.contains("4 GiB")), + "{name}: {err:?}" + ); + } + // Shrinking prunes and fills chunks: refused. + for (name, shape) in [("earray", vec![10]), ("btree2", vec![1, N + 10])] { + let err = ed.resize(name, &shape).unwrap_err(); + assert!( + matches!(&err, clawhdf5::Error::Unsupported(m) if m.contains("4 GiB")), + "{name}: {err:?}" + ); + } + assert_eq!( + std::fs::read(&path).unwrap(), + before, + "a refused edit wrote" + ); + + // Growing without early allocation touches no chunk: the dataspace + // changes, as libhdf5's `H5Dset_extent` changes it. + ed.resize("earray", &[N + 20]).unwrap(); + ed.resize("btree2", &[3, N + 10]).unwrap(); + ed.set_attr("earray", "note", &clawhdf5::AttrValue::F64(1.5)) + .unwrap(); + let file = ed.reader().unwrap(); + assert_eq!(file.dataset("earray").unwrap().shape().unwrap(), [N + 20]); + assert_eq!( + file.dataset("btree2").unwrap().shape().unwrap(), + [3, N + 10] + ); + assert!(matches!( + file.dataset("earray").unwrap().attr("note").unwrap(), + Some(clawhdf5::AttrValue::F64(v)) if v == 1.5 + )); +} + +/// The Fixed and Extensible Array structures clawhdf5 builds for 4 GiB +/// chunks are libhdf5's byte for byte: built at the fixture's addresses +/// from the fixture's chunks, they match what libhdf5 2.0.0 wrote (the +/// header's own address fields aside, which point where each writer put +/// the next block). +#[test] +fn huge_chunk_array_indexes_match_libhdf5() { + use clawhdf5_format::chunked_write::WrittenChunk; + let data = std::fs::read(fixture()).unwrap(); + let u64_at = |at: usize| u64::from_le_bytes(data[at..at + 8].try_into().unwrap()); + for name in ["farray", "earray"] { + let (_, layout, _, mut chunks) = layout_of(&data, name); + let DataLayout::Chunked { + btree_address: Some(hdr), + .. + } = layout + else { + panic!("{name}: {layout:?}"); + }; + let hdr = hdr as usize; + chunks.sort_by_key(|c| c.offsets[0]); + let slots: Vec> = chunks + .iter() + .map(|c| { + Some(WrittenChunk { + address: c.address, + compressed_size: c.chunk_size, + raw_size: N * 8, + filter_mask: c.filter_mask, + }) + }) + .collect(); + if name == "farray" { + let ours = clawhdf5_format::chunked_write::build_fixed_array_at( + &slots, 8, 8, true, hdr as u64, + ); + // FAHD up to (not including) the data block address. + assert_eq!(&ours[..16], &data[hdr..hdr + 16], "FAHD"); + let dblk = u64_at(hdr + 16) as usize; + let fadb = &ours[28..]; + assert_eq!(fadb, &data[dblk..dblk + fadb.len()], "FADB"); + } else { + let ours = clawhdf5_format::ea_writer::build_extensible_array_at( + &slots, 8, 8, true, hdr as u64, + ); + // EAHD up to (not including) the index block address. + assert_eq!(&ours[..60], &data[hdr..hdr + 60], "EAHD"); + let iblk = u64_at(hdr + 60) as usize; + let eaib = &ours[72..]; + assert_eq!(eaib, &data[iblk..iblk + eaib.len()], "EAIB"); + } + } +} + +/// h5dump from libhdf5 2.x, when `CLAWHDF5_H5DUMP2` names one (Debian's +/// h5dump is 1.14, which cannot read layout message version 5). +fn h5dump2() -> Option { + std::env::var_os("CLAWHDF5_H5DUMP2").map(PathBuf::from) +} + +/// clawhdf5 writes chunks of 4 GiB + 8 bytes (deflated) in every index it +/// uses — Single Chunk, Fixed Array, Extensible Array, v2 B-tree — with +/// layout message version 5 and 8-byte stored sizes, as libhdf5 2.x does; +/// h5py (libhdf5 2.x) and h5dump 2.x read them. Each dataset holds ten +/// values in one 4 GiB chunk, so writing it holds one 4 GiB chunk. +#[test] +fn writer_huge_chunks_round_trip() { + if !heavy() { + return; + } + const UNLIM: u64 = u64::MAX; + let dir = scratch(); + let path = dir.path().join("ours.h5"); + let values: Vec = (0..10).map(f64::from).collect(); + let mut b = clawhdf5::FileBuilder::new(); + for (name, shape, max, chunk) in [ + ("single", vec![10], vec![N], vec![N]), + ("farray", vec![10], vec![N + 10], vec![N]), + ("earray", vec![10], vec![UNLIM], vec![N]), + ("btree2", vec![1, 10], vec![UNLIM, UNLIM], vec![1, N]), + ] { + b.create_dataset(name) + .with_f64_data(&values) + .with_shape(&shape) + .with_maxshape(&max) + .with_chunks(&chunk) + .with_deflate(6); + } + b.write(&path).unwrap(); + + let data = std::fs::read(&path).unwrap(); + assert!(data.len() < 64 << 20, "{} bytes", data.len()); + let mut stored = Vec::new(); + for (name, index) in [("single", 1), ("farray", 3), ("earray", 4), ("btree2", 5)] { + let (version, layout, _, chunks) = layout_of(&data, name); + assert_eq!(version, 5, "{name}"); + assert!( + matches!(layout, DataLayout::Chunked { chunk_index_type: Some(t), .. } if t == index), + "{name}: {layout:?}" + ); + assert_eq!(chunks.len(), 1, "{name}"); + stored.push((name, chunks[0].chunk_size)); + } + + let file = File::open(&path).unwrap(); + let mut expect = values.clone(); + expect.extend([0.0; 2]); + for name in ["single", "farray", "earray"] { + let got = file.dataset(name).unwrap().read_f64().unwrap(); + assert_eq!(got, values, "{name}"); + assert_eq!(sel(&file, name, &[0], &[10]), values, "{name}"); + } + assert_eq!(sel(&file, "btree2", &[0, 3], &[1, 7]), f(3..10)); + + if python_available() { + let out = Command::new(python()) + .arg("-c") + .arg( + r#" +import sys, h5py, numpy as np +with h5py.File(sys.argv[1], "r") as f: + for name in ("single", "farray", "earray", "btree2"): + d = f[name] + assert d.chunks[-1] == 2**29 + 1, (name, d.chunks) + got = d[0] if d.ndim == 2 else d[:] + assert list(got) == list(np.arange(10.0)), (name, got) +print("ok") +"#, + ) + .arg(&path) + .output() + .unwrap(); + assert!( + out.status.success(), + "h5py:\n{}", + String::from_utf8_lossy(&out.stderr) + ); + } else { + assert!(!interop_required(), "h5py not available"); + } + + // h5dump 2.2.0 reads the layouts and walks every index: the storage + // size it reports is the chunk's. (It cannot print the values: its + // deflate filter fails on any chunk over 4 GiB, libhdf5's own included, + // with "memory allocation failed for deflate uncompression".) + if let Some(h5dump) = h5dump2() { + for (name, size) in stored { + let out = Command::new(&h5dump) + .args(["-H", "-p", "-d", name]) + .arg(&path) + .output() + .unwrap(); + let text = format!( + "{}{}", + String::from_utf8_lossy(&out.stdout), + String::from_utf8_lossy(&out.stderr) + ); + assert!(out.status.success(), "h5dump {name}:\n{text}"); + assert!(text.contains("536870913 )"), "{name}:\n{text}"); + assert!(text.contains(&format!("SIZE {size} ")), "{name}:\n{text}"); + assert!(!text.contains("rror"), "{name}:\n{text}"); + } + } +} diff --git a/crates/clawhdf5/tests/parallel_integration.rs b/crates/clawhdf5/tests/parallel_integration.rs index bdcf4f0..964d3c8 100644 --- a/crates/clawhdf5/tests/parallel_integration.rs +++ b/crates/clawhdf5/tests/parallel_integration.rs @@ -283,7 +283,7 @@ mod parallel_tests { for (i, chunk) in compressed_chunks.iter().enumerate() { file_data[offset..offset + chunk.len()].copy_from_slice(chunk); chunk_infos.push(ChunkInfo { - chunk_size: chunk.len() as u32, + chunk_size: chunk.len() as u64, filter_mask: 0, offsets: vec![(i * chunk_elems) as u64, 0], address: offset as u64, diff --git a/crates/clawhdf5/tests/partial_read_equivalence.rs b/crates/clawhdf5/tests/partial_read_equivalence.rs index a84b220..86fa568 100644 --- a/crates/clawhdf5/tests/partial_read_equivalence.rs +++ b/crates/clawhdf5/tests/partial_read_equivalence.rs @@ -195,3 +195,90 @@ fn out_of_bounds_selections_are_errors() { [99] ); } + +/// A chunked dataset with a fill value and chunks never written: a +/// selection is read over a box of fill values from the chunks it touches +/// (it used to be picked out of a full read), and through positioned reads +/// an unfiltered chunk is read row by row. Both equal the full read, in +/// every chunk index h5py writes. +#[test] +fn fill_value_selections_match_full_reads() { + let python = std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".into()); + let has_h5py = std::process::Command::new(&python) + .args(["-c", "import h5py"]) + .output() + .is_ok_and(|o| o.status.success()); + if !has_h5py { + assert!( + std::env::var("CLAWHDF5_REQUIRE_INTEROP").as_deref() != Ok("1"), + "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("fill.h5"); + let script = r#" +import sys, h5py, numpy as np +with h5py.File(sys.argv[1], "w", libver="latest") as f: + for name, maxshape, gz in [ + ("farray", None, False), ("farray_gz", None, True), + ("earray", (None, 23), False), ("earray_gz", (None, 23), True), + ("btree2", (None, None), False), ("btree2_gz", (None, None), True), + ]: + d = f.create_dataset(name, shape=(37, 23), maxshape=maxshape, chunks=(5, 4), + dtype=" 4GB requires H5F_LIBVER_V200"), which +also makes a filtered chunk index element store the chunk's size in "size +of lengths" bytes (8) instead of one byte more than the chunk needs. +clawhdf5 reads and writes such chunks in every index libhdf5 gives them +(Single Chunk, Implicit, Fixed Array, Extensible Array, v2 B-tree; libhdf5 +never puts one in a v1 B-tree and refuses to open one, as clawhdf5 does). +The [fix](#chunks-of-4-gib-or-more-could-not-be-read) is tested by +`crates/clawhdf5/tests/huge_chunks_interop.rs`; its decoding tests are +opt-in (`CLAWHDF5_HUGE_CHUNKS=1`, one at a time: `--test-threads=1`). +What remains: +- **Chunk dimensions of 2^32 or more** (which libhdf5 2.x also allows under + `H5F_LIBVER_V200`) are refused, when read (`InvalidChunkDimensions`) + and when written: `DataLayout::Chunked::chunk_dimensions` is `Vec`. + A 4 GiB chunk of 1-byte elements needs one. +- **Filters the writer refuses for such a chunk** (`FilterError`, before + anything is written): LZF, bitshuffle, bzip2 and Blosc (their HDF5 + filters record sizes or block lengths in 32 bits, or cannot take such a + buffer) and pcodec. Deflate, shuffle, Fletcher-32, LZ4 and Zstd are + written; only shuffle + deflate is tested end to end (LZ4's framing and + the 256 MiB limit it had are unit-tested; Zstd is not tested at this + size). +- **Unfiltered chunks this size are written but not tested end to end:** + `FileBuilder` assembles the whole file in memory, several copies of the + chunk (well over the 12 GiB the tests may use). Their layout messages + are unit-tested (`chunked_write::tests::huge_chunk_layout_messages_and_index_elements`). +- **`FileEditor`** refuses to rewrite such chunks (see + [its limits](#in-place-modification-fileeditor-limits)). +- **Memory:** decoding a filtered chunk holds all of it (4 GiB and more), + twice when it is shuffled (the inflated and the unshuffled copy, as + libhdf5's filters do too); writing one holds the chunk and its shuffled + copy. A selection decodes only the chunks it touches; one of an + unfiltered chunk reads only the rows it selects, through a memory map or + positioned reads (`File::open_storage`): every selection of + `unfiltered_huge_chunks_read` together read under 1 MiB. Peak resident + memory of `cargo test --release -p clawhdf5 --test huge_chunks_interop + -- --test-threads=1` with `CLAWHDF5_HUGE_CHUNKS=1` (all six tests, h5py + included) was 8.06 GiB (tank, 2026-09-28, commit `9143737`); selection + reads of the double-deflated fixture alone, 4.0 GiB. +- **32-bit targets (wasm32):** such a chunk cannot be held in memory: + reading or writing one is `FormatError::Overflow` ("exceeds the + addressable size" / "address space"); the file still opens and lists + (`examples/wasm-viewer/test/test.mjs`, "huge chunk"). +- **libhdf5 2.2.0's h5dump cannot print these values:** built here from tag + `2.2.0`, it fails to inflate any deflated chunk over 4 GiB, those + libhdf5 2.0.0 writes included ("memory allocation failed for deflate + uncompression", `H5Zdeflate.c`), while h5py 3.16 (libhdf5 2.0.0) reads + them. The tests check h5dump 2.2.0 on layouts and storage sizes only + (`CLAWHDF5_H5DUMP2`). +- libhdf5 2.0.0 itself (h5py 3.16) drops a write to an unallocated, + unfiltered Single Chunk this large when the fill time is "never" (the + dataset stays unallocated); `fixtures/gen_huge_chunks.py` allocates it + early instead. Not a clawhdf5 issue; recorded because the generator + depends on it. + ## HDF5 features still unsupported **Status:** open. What remains of the gaps the @@ -200,8 +266,16 @@ wrong data. `filter_registry::register_filter`. - **Writer:** in dense storage (more than 8 attributes on an object, or more than 8 links in a group) one attribute or link message over 65 515 - bytes is an error (no huge fractal-heap objects). The writer does not - produce output that HDF5 1.8 can read. + bytes is an error (no huge fractal-heap objects). Files HDF5 1.8 reads + are opt-in (`libver_bounds(LibVer::V18, LibVer::V18)`, since 2026-09-28, + [history](#hdf5-18-could-not-read-the-files-we-wrote)): by default the + writer uses the HDF5 1.10 format. The pre-1.8 format (version-0 + superblock, symbol-table groups; h5py's `libver='earliest'`) cannot be + written, and the Python bindings' `'w'` mode has no `libver` argument. + Under the 1.8 bound, chunk B-trees equal libhdf5's node for node when + libhdf5 inserts the chunks in row-major order; libhdf5 with its chunk + cache inserts the small chunks of a multi-dimensional dataset in eviction + order, which fills the nodes differently (same chunks, same keys). - **Checks we deliberately do not make:** - a float sign bit position outside the type, and a size-0 string type: clawhdf5 up to v2.7.0 wrote them; @@ -413,25 +487,6 @@ always expose it. A fix belongs in `conformance/ref_bugs.py` (more or more varied reads for this object) or in documenting the file as a known refusal; neither is done. -## NetCDF-4: variables' dimensions are guessed from sizes - -**Status:** open (found 2026-09-28 while fixing unlimited dimension sizes). -`clawhdf5-netcdf4` gives each variable the dimensions it finds by size -(`match_dimensions_to_variable`: the first unused dimension of equal size, -else an anonymous `dim_`), not the ones its `DIMENSION_LIST` names, and -`variables()` also lists the dimension scales that are only dimensions -(netCDF-C hides them). With netCDF4-python: unlimited `time` and `empty`, -`x` (3), `a(time)` with 2 records, `b(time, x)` with 5, `e(empty)` — netCDF4 -reports `time` = 5, `a` on `time`, and variables `a`, `b`, `e` and a `c(x)`; -clawhdf5-netcdf4 reports `time` = 5 (right) but `a` on `dim_2`, and also -variables `time` and `empty` (the scales, put on `empty`). Two dimensions of -one size can be swapped the same way. Related: netCDF4 gives `a` the shape -(5,) (a variable along an unlimited dimension has the dimension's length, -unwritten records read as fill); `Variable::shape` is the HDF5 extent, -`[2]`, and reads return those 2 values. `dimensions()` is right. -Workaround: read the variable's `_Netcdf4Coordinates` attribute (dimension -ids, matching each scale's `_Netcdf4Dimid`). - ## Small floats decode as libhdf5 does, not as the OCP MX specification **Status:** open, deliberate (documented 2026-09-28). HDF5 2.x predefines @@ -484,6 +539,164 @@ Newest first. "Before any release" means no tagged release (v2.7.0 and earlier) contains the bug. Full detail is in `CHANGELOG.md` under the date given. +## NetCDF-4: variables' dimensions are guessed from sizes + + +**Status:** fixed 2026-09-28 (#27). Affected +every release (v2.1.0 to v2.7.0: size matching dates from the crate's +first version). Wrong metadata only: stored values were always read right. +Users who read `_Netcdf4Coordinates` themselves can use +`Variable::dimensions` again; code that relied on `Variable::shape` being +the HDF5 extent, or on the reads returning only the written records, should +use `Variable::stored_shape` (new) — `shape` and the reads now follow +netCDF (below). `variables()` no longer lists pure dimension scales. + +Found 2026-09-28 while fixing unlimited dimension sizes. +`clawhdf5-netcdf4` gave each variable the dimensions it found by size +(`match_dimensions_to_variable`: the first unused dimension of equal size, +else an anonymous `dim_`), not the ones its `DIMENSION_LIST` names, and +`variables()` also listed the dimension scales that are only dimensions +(netCDF-C hides them). With netCDF4-python: unlimited `time` and `empty`, +`x` (3), `a(time)` with 2 records, `b(time, x)` with 5, `e(empty)` — netCDF4 +reports `time` = 5, `a` on `time`, and variables `a`, `b`, `e` and a `c(x)`; +clawhdf5-netcdf4 reported `time` = 5 (right) but `a` on `dim_2`, and also +variables `time` and `empty` (the scales, put on `empty`). Two dimensions of +one size could be swapped the same way. Related: netCDF4 gives `a` the shape +(5,) (a variable along an unlimited dimension has the dimension's length, +unwritten records read as fill); `Variable::shape` was the HDF5 extent, +`[2]`, and reads returned those 2 values. `dimensions()` was right. The +workaround was to read the variable's `_Netcdf4Coordinates` attribute +(dimension ids, matching each scale's `_Netcdf4Dimid`). + +Variables now get their dimensions as netCDF-C resolves them +(`libhdf5/hdf5open.c`): `_Netcdf4Coordinates` ids, else the scales +`DIMENSION_LIST` references, looked up in the variable's group and its +parents; size matching remains only for axes with neither (files not +written by a netCDF library). Pure dimension scales are not variables, and +`_nc4_non_coord_` datasets are the variables ``. A variable +along an unlimited dimension has the dimension's length, and its unwritten +records read as the fill value. One difference from netCDF-C 4.9.3 is +deliberate: when the unlimited dimension is not a variable's first +(`f(x, t)` with 1 of 4 records), a whole-variable read through netCDF-C +returns the written values first and then the fill +(`[[1, 2, fill, fill], [fill, ...]]` for rows `[1]` and `[2]`), while its +element and row reads — and clawhdf5-netcdf4 — place each row's values in +their row (`[[1, fill, fill, fill], [2, fill, fill, fill]]`). +`interop_tests` compares variables, dimensions, shapes and every value +with netCDF4-python 1.7.4 for the reproducer, dimensions of equal size, +one dimension used twice, a scalar, subgroups on their parents' +dimensions, unwritten records with and without `_FillValue`, h5py +dimension scales, and h5netcdf 1.8.1 and xarray files (tank, 2026-09-28, +`CLAWHDF5_PYTHON= cargo test -p clawhdf5-netcdf4`). +Files without dimension scales still get dimensions by size (netCDF-C +gives them `phony_dim_`), as before; see `CHANGELOG.md` for a +comparison over the conformance corpus's netCDF-readable files. + +## HDF5 1.8 could not read the files we wrote + +**Status:** fixed 2026-09-28 (#28), as an opt-in. +Affected every release (v2.1.0 to v2.7.0): the writer only ever wrote the +HDF5 1.10 format. Users who need HDF5 1.8 to read their files call +`FileBuilder::libver_bounds(LibVer::V18, LibVer::V18)` (format crate: +`FileWriter::libver_bounds`) and write them again; the default output is +unchanged. Listed until now under +[HDF5 features still unsupported](#hdf5-features-still-unsupported). + +HDF5 1.8.23's h5dump refused a default clawhdf5 file outright ("unable to +open file": its superblock is version 3; tank, 2026-09-28). The files +also use version-4 layout messages and the 1.10 chunk indexes (single +chunk, Fixed Array, Extensible Array, version-2 B-tree). libhdf5 2.0 made +the 1.8 format its default low bound (`H5F_LIBVER_V18`). + +With a low bound of 1.8 the writer now writes what libhdf5 2.x writes for +h5py's `libver=('v108', 'latest')`: a version-2 superblock, version-3 +layout messages, and a version-1 B-tree for every chunked dataset, +resizable ones included, built the way `H5B_insert` builds it (checked node +for node against libhdf5's trees). Everything else the writer emits was +already 1.8's (version-2 object headers, link messages, dense storage in +fractal heaps and version-2 B-trees, filter pipeline version 2, fill value +version 3, datatypes up to version 3). With a high bound of 1.8, what 1.8 +cannot read is `FormatError::LibverBound` before anything is written: +virtual datasets, the paged file-space strategy, the 1.12 reference types, +native complex numbers. + +`crates/clawhdf5-tools/tests/libver_v18.rs` writes every writer feature +under the 1.8 bound; HDF5 1.8.23's h5dump (built by +`scripts/build-hdf5-1.8.sh`; skipped where it is missing, as in CI) dumps +the whole file exactly as h5dump 1.14 does and returns our bytes for each +numeric dataset, and h5py, clawhdf5 and `h5rs check --data` agree; again +after `FileEditor` grows and appends to it (splitting B-tree nodes) and +sets attributes, and after h5py appends. + +The default stays the 1.10 format: on the read harness (tank, 2026-09-28, +under load from other builds) a freshly opened file with 8192 chunks +per dataset reads small selections 1.2x to 2.3x slower through a version-1 +B-tree (see `BENCHMARKS.md`, "HDF5 1.8 format"). + +## Chunks of 4 GiB or more could not be read + +**Status:** fixed 2026-09-28 (#29); affected every +release (v2.1.0 to v2.7.0). Errors, and cost; no wrong values were +returned. Nothing for users to do but upgrade. + +HDF5 2.0 writes chunks of more than 4 GiB - 1 bytes (layout message +version 5). `ChunkInfo::chunk_size` was a `u32`, so the Single Chunk, +Implicit, Fixed Array and Extensible Array readers truncated an +unfiltered chunk's size and the read failed ("incorrect chunk size +returned from index for unfiltered chunk"); a v2 B-tree index refused any +such chunk ("chunk larger than 4 GiB"), filtered or not. Filtered chunks +in the other indexes read, because their stored sizes were small. A +selection of a chunked dataset with a non-default fill value decoded the +whole dataset (8 GiB of output for the fixture's 2-D dataset), and a +selection of an unfiltered chunk fetched the whole chunk from a file that +is not in memory. `ChunkInfo::chunk_size` and `ChunkMapping::file_size` are now +`u64`; the selection path fills its box with the fill value and reads an +unfiltered chunk's rows only. Tests: `huge_chunks_interop` +(`filtered_huge_chunk_indexes_list` always; with `CLAWHDF5_HUGE_CHUNKS=1` +`filtered_huge_chunks_read`, over a fixture libhdf5 2.0.0 wrote, and +`unfiltered_huge_chunks_read`, over a 44 GiB sparse file h5py writes at +test time). Limits that remain: +[Chunks of 4 GiB or more](#chunks-of-4-gib-or-more-limits). + +## Chunks of 4 GiB or more were written unreadable + +**Status:** fixed 2026-09-28 (#29). The deflate +truncation was before any release (the one-pass deflate dates from +2026-09-23); the rest affected every release (v2.1.0 to v2.7.0). Files +clawhdf5 wrote with a chunk of 4 GiB or more should be written again. + +The writer gave such a chunk layout message version 4 (which libhdf5 +before 2.0 cannot read, and which libhdf5 2.x never writes for it; whether +2.x reads what clawhdf5 wrote was not checked), cut a chunk dimension of +2^32 or more to 32 bits, and, since 2026-09-23, deflated only the first +4 GiB - 1 bytes of the chunk (zlib takes at most that much per call and +`Finish` ended the stream there): the chunk failed to decode ("decoded to +4294967295 bytes, expected 4294967304"). It now writes layout version 5 +with libhdf5's index element widths (the Fixed and Extensible Array +structures match libhdf5 2.0.0's byte for byte: +`huge_chunk_array_indexes_match_libhdf5`), refuses chunk dimensions of +2^32 or more and the filters that cannot take such a chunk, deflates the +whole chunk, and no longer keeps the compressor's worst-case bound (4 GiB +of zeroed memory for a chunk that deflates to 4 MiB): writing the four +datasets of `writer_huge_chunks_round_trip` went past the tests' 12 GiB +cap, and compressing one shuffled 4 GiB chunk now peaks at 4.3 GiB +(tank, 2026-09-28). h5py 3.16 reads what it writes. + +## LZ4 chunks larger than 256 MiB were refused + +**Status:** fixed 2026-09-28 (#29); affected every +release (v2.1.0 to v2.7.0). An error, never wrong data. Nothing for users +to do but upgrade. + +The LZ4 decoder refused a chunk that decodes to more than 256 MiB +("lz4: declared size exceeds limit") even when the dataset's chunk size +bounded it; that ceiling is meant for a decode whose size is unknown, and +now applies only then (deflate never had it). A chunk of 4 GiB or more is +always read as the registered HDF5 framing, whose 64-bit size then no +longer starts with four zero bytes. Tests: +`filters::tests::lz4_chunks_over_256_mib_decode`, +`lz4_chunks_of_4_gib_use_the_registered_framing`. + ## A dropped `FileEditor` could keep its file locked for a moment **Status:** fixed 2026-09-28 (#23), before any release diff --git a/examples/wasm-viewer/test/test.mjs b/examples/wasm-viewer/test/test.mjs index e76ecdc..ee2adee 100644 --- a/examples/wasm-viewer/test/test.mjs +++ b/examples/wasm-viewer/test/test.mjs @@ -500,6 +500,20 @@ async function limitTests() { await fails(() => pkg.openUrl("http://huge.invalid/x.h5", { fetch: mockFetch(farBytes, { total }) }), total > 2 ** 53 ? /2\^53 - 1/ : /4 GiB/, `length ${total}`); } + + // Nor can it hold a chunk of 4 GiB or more (HDF5 2.0, layout message + // version 5): the file opens and lists, and reading such a chunk is an + // error naming the size, in every chunk index. + const hugeChunks = join(import.meta.dirname, "..", "..", "..", + "crates/clawhdf5/tests/fixtures/huge_chunks_filtered.h5"); + const hc = pkg.open(new Uint8Array(readFileSync(hugeChunks))); + eq(hc.list("/").map((e) => e.name), ["btree2", "earray", "farray", "single"], "huge chunks: list"); + for (const name of ["single", "farray", "earray", "btree2"]) { + const [start, count] = name === "btree2" ? [[0, 0], [1, 4]] : [[0], [4]]; + await fails(() => hc.readHyperslab(`/${name}`, start, count), /exceeds the addressable size/, + `huge chunk: ${name}`); + } + hc.free(); } // A body of `total` bytes in 64 KiB pieces, made as they are read; `pulled()` diff --git a/scripts/build-hdf5-1.8.sh b/scripts/build-hdf5-1.8.sh new file mode 100644 index 0000000..cc75fc0 --- /dev/null +++ b/scripts/build-hdf5-1.8.sh @@ -0,0 +1,43 @@ +#!/usr/bin/env bash +# Build libhdf5 1.8.23 (the last 1.8 release) with its command-line tools, as +# the oracle for files written with `LibVer::V18` (the `libver_v18` test in +# clawhdf5-tools finds h5dump through CLAWHDF5_H5DUMP18 or this default prefix). +# +# bash scripts/build-hdf5-1.8.sh [PREFIX] +# +# PREFIX defaults to ~/.cache/hdf5-1.8.23; sources go to PREFIX-src and the +# build tree to PREFIX-build. Reuses an existing install. Needs git, cmake, a C +# compiler and zlib headers. +set -euo pipefail +PREFIX="${1:-$HOME/.cache/hdf5-1.8.23}" +SRC="$PREFIX-src" +BUILD="$PREFIX-build" +if [ -x "$PREFIX/bin/h5dump" ]; then + echo "reusing $PREFIX/bin/h5dump" + "$PREFIX/bin/h5dump" --version + exit 0 +fi +if [ ! -d "$SRC" ]; then + git clone --depth 1 --branch hdf5-1_8_23 https://github.com/HDFGroup/hdf5.git "$SRC" +fi +# 1.8 predates current compilers (GCC 14 turns its pointer-type mismatches in +# the tools into errors): pin gnu99 and demote those errors to warnings. +CFLAGS18="-w -std=gnu99 -Wno-error=incompatible-pointer-types" +CFLAGS18="$CFLAGS18 -Wno-error=implicit-function-declaration -Wno-error=int-conversion" +cmake -S "$SRC" -B "$BUILD" \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_INSTALL_PREFIX="$PREFIX" \ + -DCMAKE_C_FLAGS="$CFLAGS18" \ + -DBUILD_SHARED_LIBS=ON \ + -DBUILD_TESTING=OFF \ + -DHDF5_BUILD_TOOLS=ON \ + -DHDF5_BUILD_EXAMPLES=OFF \ + -DHDF5_BUILD_CPP_LIB=OFF \ + -DHDF5_BUILD_FORTRAN=OFF \ + -DHDF5_BUILD_HL_LIB=OFF \ + -DHDF5_BUILD_JAVA=OFF \ + -DHDF5_ENABLE_Z_LIB_SUPPORT=ON \ + -DHDF5_ENABLE_SZIP_SUPPORT=OFF +cmake --build "$BUILD" -j "${JOBS:-6}" +cmake --install "$BUILD" +"$PREFIX/bin/h5dump" --version