From d0d5347cd9478935af7f26cea7dd6117f4d8c733 Mon Sep 17 00:00:00 2001 From: osobh Date: Mon, 28 Sep 2026 21:16:17 -0500 Subject: [PATCH] feat: write complex numbers, incl. HDF5 2.0 native complex (class 11) - Datatype::Complex serializes class 11 version 5 byte-identically to libhdf5 2.2.0; containers holding it are written as version 5. - DatasetBuilder::with_complex_f32/f64_data (h5py's {r, i} compound, default) and with_native_complex_f32/f64_data (class 11, opt-in); make_(native_)complex_f32/f64_type for attributes. - Dataset::read_complex_f64/f32 read either form. - Python create_dataset accepts complex64/complex128 (compound form). - Parsing unchanged: class 11 still surfaces as {r, i}. - Tests vs h5py 3.16 / libhdf5 2.0.0 and h5dump 2.2.0; docs. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 39 ++++ README.md | 4 +- crates/clawhdf5-format/src/data_read.rs | 4 + crates/clawhdf5-format/src/datatype.rs | 193 ++++++++++++++++-- crates/clawhdf5-format/src/type_builders.rs | 95 +++++++++ .../tests/writer_h5py_tests.rs | 163 +++++++++++++++ crates/clawhdf5-py/README.md | 7 +- crates/clawhdf5-py/src/convert.rs | 3 + crates/clawhdf5-py/src/edit.rs | 5 +- crates/clawhdf5-py/src/lib.rs | 24 ++- crates/clawhdf5-py/tests/test_write_read.py | 45 ++++ crates/clawhdf5-tools/src/dtype.rs | 10 + crates/clawhdf5-tools/src/value.rs | 5 + crates/clawhdf5-wasm/src/core.rs | 3 + crates/clawhdf5/src/lib.rs | 5 +- crates/clawhdf5/src/reader.rs | 54 +++++ crates/clawhdf5/src/vlen.rs | 1 + crates/clawhdf5/tests/h5py_interop_tests.rs | 48 +++++ crates/clawhdf5/tests/integration_tests.rs | 94 +++++++++ docs/QUICKSTART.md | 3 +- docs/known-issues.md | 9 + 21 files changed, 790 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b2a578c..6590714 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -59,6 +59,45 @@ `'r+'` raises `NotImplementedError` before anything is written; inside a compound or array type they still raise `TypeError`. +### Writing complex numbers, including HDF5 2.0's native complex type (2026-09-28) +- **h5py's form (default):** `DatasetBuilder::with_complex_f32_data` / + `with_complex_f64_data` take `[re, im]` pairs and write the compound + `{r, i}` h5py writes for numpy `complex64`/`complex128` (h5py 3.16 still + writes this by default); `make_complex_f32_type`/`make_complex_f64_type` + give the datatype for `AttrValue::Raw` attributes. +- **Native complex (opt-in):** `with_native_complex_f32_data` / + `with_native_complex_f64_data` and `make_native_complex_f32_type` / + `make_native_complex_f64_type` write `H5T_COMPLEX_IEEE_F32LE`/`F64LE` + (datatype class 11, version 5), through the new `Datatype::Complex` + variant. The encoding is byte-identical to libhdf5 2.2.0's + (`H5Odtype.c`: homogeneous, rectangular, base type follows), and a + compound, array or variable-length type holding one is written as + version 5, as libhdf5 raises it. No file-level version bound is needed: + the superblock and object headers we write already open in libhdf5 2.0. + Only libhdf5 2.0+ reads class 11 (Debian's h5dump 1.14 fails on the + object), hence opt-in. Checked on tank: h5py 3.16 (libhdf5 2.0.0) reads + datasets (contiguous and chunked+deflate), attributes, a compound member + and an array of them as numpy `complex64`/`complex128` with class 11; + h5dump 2.2.0 prints them as `H5T_COMPLEX_IEEE_F*LE`; `h5rs check --data` + finds no problems. +- **Reading:** `Dataset::read_complex_f64`/`read_complex_f32` return + `[re, im]` pairs from either form (and from h5py's files, native or + not). Parsing is unchanged: class 11 still surfaces as the `{r, i}` + compound (`Datatype::complex_as_compound`), so `Datatype::parse` never + returns `Datatype::Complex` (see `docs/known-issues.md`). +- **Breaking for exhaustive matches:** `Datatype` gained the `Complex` + variant; downstream `match`es over `Datatype` without a wildcard need an + arm (`Datatype::complex_as_compound(size, base)` gives the compound view). +- **Python:** `create_dataset` accepts `complex64` and `complex128` arrays, + written as h5py's compound (native class 11 is Rust-only). +- Tests: `datatype.rs` (byte equality with libhdf5's encoding), + `integration_tests::complex_datasets_and_attributes_round_trip`, + `writer_h5py_tests::h5py_reads_our_complex_datasets_and_attributes` + (runs h5dump 2.x when `CLAWHDF5_H5DUMP2` names one), + `h5py_interop_tests::h5py_complex_datasets_read_as_complex`, and + `test_write_read.py::test_roundtrip_complex` / + `test_read_native_complex_from_h5py`. + ### `ObjectHeader::parse` back at its pre-M2/M3 speed (2026-09-27) - Parsing a version-1 object header was 4% slower than before range-read M2/M3 (`docs/known-issues.md`). The cause was the call to the per-chunk diff --git a/README.md b/README.md index 7c34ff5..bbef7ac 100644 --- a/README.md +++ b/README.md @@ -115,11 +115,13 @@ Limits and open issues, with dates, are in | **File format** | Superblock v0–v3, user blocks, v1/v2 object headers | Metadata cache images | Writing files HDF5 1.8 can read | | **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 (HDF5 2.0 class 11) | 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 | + +| **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 | Writing variable-length data; 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) | | **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) | | **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 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) | +| **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) | Plugin filters other than LZF are cargo features (`bitshuffle`, `bzip2`, `blosc`, `blosc2`, `zfp`, or `plugin-filters` for all of them), all pure diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 588d568..47e0c58 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -815,6 +815,7 @@ fn datatype_name(dt: &Datatype) -> &'static str { Datatype::Enumeration { .. } => "Enumeration", Datatype::VariableLength { .. } => "VariableLength", Datatype::Array { .. } => "Array", + Datatype::Complex { .. } => "Complex", } } @@ -1566,6 +1567,9 @@ pub fn read_compound_fields( datatype: &Datatype, ) -> Result, FormatError> { match datatype { + Datatype::Complex { size, base_type } => { + read_compound_fields(raw, &Datatype::complex_as_compound(*size, base_type)) + } Datatype::Compound { size, members } => { let elem_size = *size as usize; if elem_size == 0 { diff --git a/crates/clawhdf5-format/src/datatype.rs b/crates/clawhdf5-format/src/datatype.rs index 974c249..8e2120d 100644 --- a/crates/clawhdf5-format/src/datatype.rs +++ b/crates/clawhdf5-format/src/datatype.rs @@ -140,6 +140,19 @@ pub enum Datatype { base_type: Box, dimensions: Vec, }, + /// Class 11: HDF5 2.0 native complex number (`H5T_COMPLEX`, datatype + /// message version 5): two consecutive `base_type` values, real then + /// imaginary, in rectangular form. `size` is twice the base size and the + /// base is an IEEE float. + /// + /// This variant exists for **writing** (see + /// `type_builders::make_native_complex_f64_type`): only libhdf5 2.0 and + /// newer can read class 11, so it is opt-in and h5py's compound `{r, i}` + /// stays the default complex encoding. [`Datatype::parse`] still + /// surfaces a class-11 message as that equivalent `{r, i}` compound, so + /// every compound reader handles both encodings; parsing what this + /// variant serializes therefore yields a `Compound`, not a `Complex`. + Complex { size: u32, base_type: Box }, } /// Longest opaque tag that can be stored: its NUL-padded length must fit @@ -862,19 +875,7 @@ impl Datatype { actual: size as usize, }); } - let members = vec![ - CompoundMember { - name: String::from("r"), - byte_offset: 0, - datatype: base_type.clone(), - }, - CompoundMember { - name: String::from("i"), - byte_offset: base_size as u64, - datatype: base_type, - }, - ]; - Ok((Datatype::Compound { size, members }, pos)) + Ok((Self::complex_as_compound(size, &base_type), pos)) } _ => Err(FormatError::InvalidDatatypeClass(class_id)), } @@ -939,7 +940,8 @@ impl Datatype { .try_for_each(|m| m.datatype.check_unused_bits()), Datatype::Enumeration { base_type, .. } | Datatype::VariableLength { base_type, .. } - | Datatype::Array { base_type, .. } => base_type.check_unused_bits(), + | Datatype::Array { base_type, .. } + | Datatype::Complex { base_type, .. } => base_type.check_unused_bits(), _ => Ok(()), } } @@ -1046,7 +1048,8 @@ impl Datatype { } else { 0 }; - let mut buf = Self::build_header(9, 1, [bf0, bf1, 0], *size); + let version = base_type.min_parent_version().max(1); + let mut buf = Self::build_header(9, version, [bf0, bf1, 0], *size); buf.extend_from_slice(&base_type.serialize()); buf } @@ -1054,7 +1057,11 @@ impl Datatype { let num = members.len() as u16; let bf0 = (num & 0xFF) as u8; let bf1 = ((num >> 8) & 0xFF) as u8; - let mut buf = Self::build_header(6, 3, [bf0, bf1, 0], *size); + let version = members + .iter() + .map(|m| m.datatype.min_parent_version()) + .fold(3, u8::max); + let mut buf = Self::build_header(6, version, [bf0, bf1, 0], *size); let ob = offset_bytes_for_size(*size); for m in members { // Null-terminated name @@ -1097,7 +1104,8 @@ impl Datatype { base_type, dimensions, } => { - let mut buf = Self::build_header(10, 3, [0, 0, 0], self.type_size()); + let version = base_type.min_parent_version().max(3); + let mut buf = Self::build_header(10, version, [0, 0, 0], self.type_size()); buf.push(dimensions.len() as u8); for &d in dimensions { buf.extend_from_slice(&d.to_le_bytes()); @@ -1152,6 +1160,59 @@ impl Datatype { }; Self::build_header(7, version, [bf0, 0, 0], *size) } + Datatype::Complex { size, base_type } => { + // Version 5 (HDF5 2.0), as libhdf5's `H5O__dtype_encode_helper` + // writes it: bit 0 = homogeneous (the only kind libhdf5 + // supports), bits 1-2 = form (0, rectangular); the base + // datatype message follows. + let mut buf = Self::build_header(11, 5, [0x01, 0, 0], *size); + buf.extend_from_slice(&base_type.serialize()); + buf + } + } + } + + /// The `{r, i}` compound equivalent to a native complex type of `size` + /// bytes over `base_type`: `r` at offset 0, `i` right after it — the + /// shape h5py writes for numpy complex dtypes, and what [`Self::parse`] + /// returns for a class-11 message. Readers that meet a + /// [`Datatype::Complex`] handle it through this view. + pub fn complex_as_compound(size: u32, base_type: &Datatype) -> Datatype { + let base_size = base_type.type_size(); + Datatype::Compound { + size, + members: vec![ + CompoundMember { + name: String::from("r"), + byte_offset: 0, + datatype: base_type.clone(), + }, + CompoundMember { + name: String::from("i"), + byte_offset: u64::from(base_size), + datatype: base_type.clone(), + }, + ], + } + } + + /// The lowest datatype message version a type that contains this one + /// may be encoded with. libhdf5 raises a compound, array, variable-length + /// or enum type to the version of its members (`H5O_DTYPE_CHECK_VERSION` + /// in `H5Odtype.c`), so a type holding a native complex (version 5) is + /// itself written as version 5; everything else we write keeps the + /// container's own version. + fn min_parent_version(&self) -> u8 { + match self { + Datatype::Complex { .. } => 5, + Datatype::Compound { members, .. } => members + .iter() + .map(|m| m.datatype.min_parent_version()) + .fold(0, u8::max), + Datatype::Enumeration { base_type, .. } + | Datatype::VariableLength { base_type, .. } + | Datatype::Array { base_type, .. } => base_type.min_parent_version(), + _ => 0, } } @@ -1189,6 +1250,20 @@ impl Datatype { Datatype::Enumeration { base_type, .. } | Datatype::VariableLength { base_type, .. } | Datatype::Array { base_type, .. } => base_type.check_encodable_parts(), + // libhdf5 only builds complex types over IEEE floats + // (`H5Tcomplex_create`), always twice the base size. + Datatype::Complex { size, base_type } => match base_type.as_ref() { + Datatype::FloatingPoint { size: b, .. } if b.checked_mul(2) == Some(*size) => { + Ok(()) + } + Datatype::FloatingPoint { .. } => Err(FormatError::SerializationError(format!( + "complex datatype of size {size} is not twice its base size {}", + base_type.type_size() + ))), + _ => Err(FormatError::SerializationError( + "complex datatype base must be a floating-point type".into(), + )), + }, _ => Ok(()), } } @@ -1216,6 +1291,7 @@ impl Datatype { Datatype::Reference { size, .. } => *size, Datatype::Enumeration { size, .. } => *size, Datatype::VariableLength { size, .. } => *size, + Datatype::Complex { size, .. } => *size, Datatype::Array { base_type, dimensions, @@ -1815,6 +1891,89 @@ mod tests { )); } + #[test] + fn native_complex_serializes_as_libhdf5_2_0_does() { + use crate::type_builders::{make_native_complex_f32_type, make_native_complex_f64_type}; + let dt = make_native_complex_f64_type(); + assert_eq!(dt.serialize(), COMPLEX_F64_HDF5_2_0); + assert_eq!(dt.type_size(), 16); + dt.check_encodable().unwrap(); + // Parsing surfaces class 11 as the equivalent `{r, i}` compound. + let (parsed, _) = Datatype::parse(&dt.serialize()).unwrap(); + assert_eq!( + parsed, + Datatype::complex_as_compound(16, &crate::type_builders::make_f64_type()) + ); + + let f32c = make_native_complex_f32_type().serialize(); + assert_eq!( + &f32c[..8], + &[0x5b, 0x01, 0x00, 0x00, 0x08, 0x00, 0x00, 0x00] + ); + assert_eq!( + &f32c[8..], + &crate::type_builders::make_f32_type().serialize()[..] + ); + } + + #[test] + fn compound_holding_native_complex_serializes_as_libhdf5_2_0_does() { + // The same type as `test_compound_with_complex_member_from_hdf5_2_0`: + // libhdf5 raises the compound to version 5 for its complex member. + let dt = Datatype::Compound { + size: 24, + members: vec![ + CompoundMember { + name: "z".into(), + byte_offset: 0, + datatype: crate::type_builders::make_native_complex_f64_type(), + }, + CompoundMember { + name: "k".into(), + byte_offset: 16, + datatype: crate::type_builders::make_i64_type(), + }, + ], + }; + let mut want = vec![ + 0x56, 0x02, 0x00, 0x00, 0x18, 0x00, 0x00, 0x00, b'z', 0x00, 0x00, + ]; + want.extend_from_slice(&COMPLEX_F64_HDF5_2_0); + want.extend_from_slice(&[b'k', 0x00, 0x10]); + want.extend_from_slice(&[ + 0x10, 0x08, 0x00, 0x00, 0x08, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40, 0x00, + ]); + assert_eq!(dt.serialize(), want); + + // An array of complex is raised to version 5 as well; one without + // stays at version 3. + let arr = Datatype::Array { + base_type: Box::new(crate::type_builders::make_native_complex_f32_type()), + dimensions: vec![2], + }; + assert_eq!(arr.serialize()[0], 0x5a); + arr.check_encodable().unwrap(); + let plain = Datatype::Array { + base_type: Box::new(crate::type_builders::make_f32_type()), + dimensions: vec![2], + }; + assert_eq!(plain.serialize()[0], 0x3a); + } + + #[test] + fn native_complex_must_be_twice_an_ieee_float() { + let bad_size = Datatype::Complex { + size: 12, + base_type: Box::new(crate::type_builders::make_f64_type()), + }; + assert!(bad_size.check_encodable().is_err()); + let int_base = Datatype::Complex { + size: 8, + base_type: Box::new(crate::type_builders::make_i32_type()), + }; + assert!(int_base.check_encodable().is_err()); + } + #[test] fn test_reference_object() { let buf = build_dt_header(7, 1, [0, 0, 0], 8); diff --git a/crates/clawhdf5-format/src/type_builders.rs b/crates/clawhdf5-format/src/type_builders.rs index 9d32e5c..ba8c292 100644 --- a/crates/clawhdf5-format/src/type_builders.rs +++ b/crates/clawhdf5-format/src/type_builders.rs @@ -137,6 +137,37 @@ pub fn make_f32_type() -> Datatype { } } +/// numpy `complex64` the way h5py stores it: a compound `{r: f32, i: f32}`. +/// Every HDF5 reader opens it; h5py reads it back as `complex64`. +pub fn make_complex_f32_type() -> Datatype { + Datatype::complex_as_compound(8, &make_f32_type()) +} + +/// numpy `complex128` the way h5py stores it: a compound `{r: f64, i: f64}`. +pub fn make_complex_f64_type() -> Datatype { + Datatype::complex_as_compound(16, &make_f64_type()) +} + +/// HDF5 2.0's native complex type `H5T_COMPLEX_IEEE_F32LE` (datatype class +/// 11). Only libhdf5 2.0 and newer (h5py built on it) can read a file that +/// uses it; older libhdf5, including h5dump 1.14, refuses the object. +/// Prefer [`make_complex_f32_type`] unless the consumer wants class 11. +pub fn make_native_complex_f32_type() -> Datatype { + Datatype::Complex { + size: 8, + base_type: Box::new(make_f32_type()), + } +} + +/// HDF5 2.0's native complex type `H5T_COMPLEX_IEEE_F64LE` (datatype class +/// 11); see [`make_native_complex_f32_type`] for who can read it. +pub fn make_native_complex_f64_type() -> Datatype { + Datatype::Complex { + size: 16, + base_type: Box::new(make_f64_type()), + } +} + pub fn make_i32_type() -> Datatype { Datatype::FixedPoint { size: 4, @@ -640,6 +671,70 @@ impl DatasetBuilder { self } + /// Store complex numbers, each `[re, im]`, as h5py does for numpy + /// `complex64`: a compound `{r, i}` of `f32` ([`make_complex_f32_type`]), + /// readable by every HDF5 library. For HDF5 2.0's native complex type use + /// [`Self::with_native_complex_f32_data`]. + pub fn with_complex_f32_data(&mut self, data: &[[f32; 2]]) -> &mut Self { + self.set_complex( + make_complex_f32_type(), + data.as_flattened(), + f32::to_le_bytes, + ) + } + + /// Store complex numbers, each `[re, im]`, as h5py does for numpy + /// `complex128`: a compound `{r, i}` of `f64` ([`make_complex_f64_type`]). + pub fn with_complex_f64_data(&mut self, data: &[[f64; 2]]) -> &mut Self { + self.set_complex( + make_complex_f64_type(), + data.as_flattened(), + f64::to_le_bytes, + ) + } + + /// Store complex numbers, each `[re, im]`, as HDF5 2.0's native complex + /// type `H5T_COMPLEX_IEEE_F32LE` (datatype class 11). The bytes are the + /// same as [`Self::with_complex_f32_data`]; only the datatype differs. + /// h5py on libhdf5 2.0+ reads it as `complex64`; libhdf5 1.x cannot open + /// the dataset at all, so this is opt-in. + pub fn with_native_complex_f32_data(&mut self, data: &[[f32; 2]]) -> &mut Self { + self.set_complex( + make_native_complex_f32_type(), + data.as_flattened(), + f32::to_le_bytes, + ) + } + + /// Store complex numbers, each `[re, im]`, as HDF5 2.0's native complex + /// type `H5T_COMPLEX_IEEE_F64LE` (class 11); see + /// [`Self::with_native_complex_f32_data`]. + pub fn with_native_complex_f64_data(&mut self, data: &[[f64; 2]]) -> &mut Self { + self.set_complex( + make_native_complex_f64_type(), + data.as_flattened(), + f64::to_le_bytes, + ) + } + + fn set_complex( + &mut self, + datatype: Datatype, + parts: &[T], + le: fn(T) -> [u8; N], + ) -> &mut Self { + self.datatype = Some(datatype); + let mut b = Vec::with_capacity(parts.len() * N); + for &v in parts { + b.extend_from_slice(&le(v)); + } + self.data = Some(b); + if self.shape.is_none() { + self.shape = Some(vec![(parts.len() / 2) as u64]); + } + self + } + /// Write a compound (struct) dataset. pub fn with_compound_data( &mut self, diff --git a/crates/clawhdf5-format/tests/writer_h5py_tests.rs b/crates/clawhdf5-format/tests/writer_h5py_tests.rs index cfbaaa8..a4a0d33 100644 --- a/crates/clawhdf5-format/tests/writer_h5py_tests.rs +++ b/crates/clawhdf5-format/tests/writer_h5py_tests.rs @@ -383,6 +383,169 @@ else: assert_eq!((fields[1].name.as_str(), im), ("i", vec![2.0, 4.0])); } +/// Complex data both ways we write it: h5py's compound `{r, i}` (the +/// default, readable everywhere) and HDF5 2.0's native complex type (class +/// 11, opt-in). h5py on libhdf5 2.0+ must read the native datasets and +/// attributes as numpy `complex64`/`complex128`, and a compound holding a +/// native complex member (written as datatype version 5, as libhdf5 does). +/// Skips the native checks when h5py's libhdf5 predates 2.0; also runs +/// h5dump 2.x when `CLAWHDF5_H5DUMP2` names one. +#[test] +#[ignore = "requires Python h5py module"] +fn h5py_reads_our_complex_datasets_and_attributes() { + use clawhdf5_format::type_builders::{ + make_complex_f64_type, make_native_complex_f32_type, make_native_complex_f64_type, + }; + let path = std::env::temp_dir().join("clawhdf5_test_complex.h5"); + let native_ok = h5py_read( + &path, + "import h5py; print(int(getattr(h5py.get_config(), 'has_native_complex', False)))", + ) == "1"; + + let z64 = [[1.5f32, -2.0], [0.0, 3.25], [-7.0, 1.0e-3]]; + let z128 = [[1.0f64, 2.0], [-3.5, 4.0e300], [0.25, -0.0]]; + let big: Vec<[f64; 2]> = (0..600).map(|k| [k as f64, -(k as f64) / 4.0]).collect(); + let c128 = |re: f64, im: f64| [re.to_le_bytes(), im.to_le_bytes()].concat(); + + let mut fw = FileWriter::new(); + fw.create_dataset("compound128") + .with_complex_f64_data(&z128) + .set_attr( + "c", + AttrValue::Raw { + datatype: make_complex_f64_type(), + shape: vec![], + data: c128(1.0, -1.0), + }, + ); + if native_ok { + fw.create_dataset("native64") + .with_native_complex_f32_data(&z64) + .set_attr( + "c", + AttrValue::Raw { + datatype: make_native_complex_f64_type(), + shape: vec![2], + data: [c128(0.5, -1.5), c128(2.0, 3.0)].concat(), + }, + ); + fw.create_dataset("native128") + .with_native_complex_f64_data(&z128); + fw.create_dataset("native_chunked") + .with_native_complex_f64_data(&big) + .with_shape(&[20, 30]) + .with_chunks(&[7, 16]) + .with_deflate(4); + // A compound with a native complex member, and an array of them. + let rec = CompoundTypeBuilder::new() + .field("z", make_native_complex_f64_type()) + .i64_field("k") + .build(); + let raw = [ + c128(1.0, 2.0), + 7i64.to_le_bytes().to_vec(), + c128(-1.0, 0.5), + (-8i64).to_le_bytes().to_vec(), + ] + .concat(); + fw.create_dataset("records").with_compound_data(rec, raw, 2); + let arr: Vec = [1.0f32, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0] + .iter() + .flat_map(|v| v.to_le_bytes()) + .collect(); + fw.create_dataset("arrays") + .with_array_data(make_native_complex_f32_type(), &[2], arr, 2); + fw.set_root_attr( + "zroot", + AttrValue::Raw { + datatype: make_native_complex_f32_type(), + shape: vec![], + data: [4.0f32.to_le_bytes(), (-4.0f32).to_le_bytes()].concat(), + }, + ); + } + std::fs::write(&path, fw.finish().unwrap()).unwrap(); + + let script = format!( + r#" +import h5py, json, numpy as np +f = h5py.File('{}', 'r') +def z(a): return [[float(np.real(v)), float(np.imag(v))] for v in np.asarray(a).ravel()] +out = {{}} +d = f['compound128'] +out['compound128'] = [str(d.dtype), z(d[()]), str(d.attrs['c'].dtype), z(d.attrs['c'])] +if {native}: + for n in ['native64', 'native128', 'native_chunked']: + d = f[n] + out[n] = [str(d.dtype), list(d.shape), z(d[()]), d.id.get_type().get_class()] + a = f['native64'].attrs['c'] + out['attr'] = [str(a.dtype), z(a)] + r = f['records'][()] + out['records'] = [str(r.dtype['z']), z(r['z']), r['k'].tolist()] + a = f['arrays'][()] + out['arrays'] = [str(a.dtype), list(a.shape), z(a)] + a = f.attrs['zroot'] + out['zroot'] = [str(a.dtype), z(a)] + out['CLASS'] = h5py.h5t.COMPLEX +print(json.dumps(out)) +"#, + path.display(), + native = if native_ok { "True" } else { "False" }, + ); + let v: serde_json::Value = serde_json::from_str(&h5py_read(&path, &script)).unwrap(); + let pairs = |p: &[[f64; 2]]| serde_json::json!(p); + assert_eq!(v["compound128"][0], "complex128"); + assert_eq!(v["compound128"][1], pairs(&z128)); + assert_eq!(v["compound128"][2], "complex128"); + assert_eq!(v["compound128"][3], pairs(&[[1.0, -1.0]])); + if !native_ok { + eprintln!("HDF5 < 2.0: native complex not checked"); + return; + } + let z64_wide: Vec<[f64; 2]> = z64.iter().map(|p| [p[0].into(), p[1].into()]).collect(); + let class_complex = v["CLASS"].clone(); + for (name, dtype, shape, values) in [ + ("native64", "complex64", vec![3], z64_wide.clone()), + ("native128", "complex128", vec![3], z128.to_vec()), + ("native_chunked", "complex128", vec![20, 30], big.clone()), + ] { + assert_eq!(v[name][0], dtype, "{name}"); + assert_eq!(v[name][1], serde_json::json!(shape), "{name}"); + assert_eq!(v[name][2], pairs(&values), "{name}"); + // Stored as class 11, not converted from a compound. + assert_eq!(v[name][3], class_complex, "{name}"); + } + assert_eq!(v["attr"][0], "complex128"); + assert_eq!(v["attr"][1], pairs(&[[0.5, -1.5], [2.0, 3.0]])); + assert_eq!(v["records"][0], "complex128"); + assert_eq!(v["records"][1], pairs(&[[1.0, 2.0], [-1.0, 0.5]])); + assert_eq!(v["records"][2], serde_json::json!([7, -8])); + assert_eq!(v["arrays"][0], "complex64"); + assert_eq!(v["arrays"][1], serde_json::json!([2, 2])); + assert_eq!( + v["arrays"][2], + pairs(&[[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]]) + ); + assert_eq!(v["zroot"][0], "complex64"); + assert_eq!(v["zroot"][1], pairs(&[[4.0, -4.0]])); + + // h5dump from libhdf5 2.x, when available (Debian's 1.14 cannot read + // class 11 at all). + if let Ok(h5dump) = std::env::var("CLAWHDF5_H5DUMP2") { + let o = std::process::Command::new(&h5dump) + .arg(&path) + .output() + .expect("run h5dump"); + let text = String::from_utf8_lossy(&o.stdout); + assert!( + o.status.success(), + "h5dump: {}", + String::from_utf8_lossy(&o.stderr) + ); + assert!(text.contains("H5T_COMPLEX"), "{text}"); + } +} + #[test] #[ignore = "requires Python h5py module"] fn read_h5py_generated_enum() { diff --git a/crates/clawhdf5-py/README.md b/crates/clawhdf5-py/README.md index 3a62261..c05f97e 100644 --- a/crates/clawhdf5-py/README.md +++ b/crates/clawhdf5-py/README.md @@ -100,8 +100,11 @@ f.remote_stats # {'requests': ..., 'bytes_fetched': ..., 'hits': ..., ...} `clawhdf5.File(path, "w")` with `create_dataset(name, data=array, chunks=..., compression="gzip")`, `create_group` and `attrs[...] = ...` -writes `float64`, `float32`, `int64`, `int32` and `uint8` arrays; the file is -written on `close()`. +writes `float64`, `float32`, `int64`, `int32`, `uint8`, `complex64` and +`complex128` arrays; the file is written on `close()`. Complex arrays are +stored as h5py stores them, a compound `{r, i}` that every libhdf5 reads +(not HDF5 2.0's native complex type, which only libhdf5 2.0+ reads; the +Rust API writes that on request). ## Editing a file in place diff --git a/crates/clawhdf5-py/src/convert.rs b/crates/clawhdf5-py/src/convert.rs index b93675b..8e14e6d 100644 --- a/crates/clawhdf5-py/src/convert.rs +++ b/crates/clawhdf5-py/src/convert.rs @@ -313,6 +313,9 @@ fn np_dtype_with_metadata<'py>( /// read as they are. pub(crate) fn fixed_dtype<'py>(py: Python<'py>, dt: &Datatype) -> PyResult> { match dt { + Datatype::Complex { size, base_type } => { + fixed_dtype(py, &Datatype::complex_as_compound(*size, base_type)) + } Datatype::FixedPoint { .. } => np_dtype(py, int_format(dt)?), Datatype::FloatingPoint { .. } => np_dtype(py, float_format(dt)?), Datatype::String { size, charset, .. } => { diff --git a/crates/clawhdf5-py/src/edit.rs b/crates/clawhdf5-py/src/edit.rs index 6646a4d..8c77d4b 100644 --- a/crates/clawhdf5-py/src/edit.rs +++ b/crates/clawhdf5-py/src/edit.rs @@ -47,6 +47,7 @@ fn not_implemented(what: impl std::fmt::Display) -> PyErr { /// `edit_helpers._convert_array`), or why they cannot be written. pub(crate) fn category(dt: &Datatype) -> PyResult<&'static str> { match dt { + Datatype::Complex { .. } => Ok("complex"), Datatype::FixedPoint { .. } => Ok("int"), Datatype::FloatingPoint { .. } if crate::convert::is_ieee_float(dt) => Ok("float"), // Read as a wider IEEE float; writing would need the reverse @@ -107,7 +108,9 @@ fn check_exact(dt: &Datatype) -> PyResult<()> { Datatype::Compound { members, .. } => { members.iter().try_for_each(|m| check_exact(&m.datatype)) } - Datatype::Array { base_type, .. } => check_exact(base_type), + Datatype::Array { base_type, .. } | Datatype::Complex { base_type, .. } => { + check_exact(base_type) + } Datatype::String { padding: StringPadding::NullPad, .. diff --git a/crates/clawhdf5-py/src/lib.rs b/crates/clawhdf5-py/src/lib.rs index 4ae81f7..3ffb73b 100644 --- a/crates/clawhdf5-py/src/lib.rs +++ b/crates/clawhdf5-py/src/lib.rs @@ -161,6 +161,10 @@ pub(crate) enum DatasetData { I64(Vec), I32(Vec), U8(Vec), + /// numpy `complex64`, `[re, im]` pairs, written as h5py does. + C64(Vec<[f32; 2]>), + /// numpy `complex128`, `[re, im]` pairs, written as h5py does. + C128(Vec<[f64; 2]>), } /// Specification for a dataset to be written. @@ -281,6 +285,14 @@ pub(crate) fn apply_dataset_spec( DatasetData::U8(v) => { db.with_u8_data(v); } + // h5py's compound `{r, i}`, not HDF5 2.0's native complex type: it + // is what h5py writes (3.16 included) and every libhdf5 can read it. + DatasetData::C64(v) => { + db.with_complex_f32_data(v); + } + DatasetData::C128(v) => { + db.with_complex_f64_data(v); + } } if !spec.shape.is_empty() { db.with_shape(&spec.shape); @@ -314,9 +326,19 @@ pub(crate) fn extract_numpy_data( "int64" => DatasetData::I64(flat.extract::>()?), "int32" => DatasetData::I32(flat.extract::>()?), "uint8" => DatasetData::U8(flat.extract::>()?), + // A 1-D complex array viewed as floats is its (re, im) parts in order. + "complex64" => { + let parts: Vec = flat.call_method1("view", ("().0.to_vec()) + } + "complex128" => { + let parts: Vec = flat.call_method1("view", ("().0.to_vec()) + } _ => { return Err(PyErr::new::(format!( - "unsupported numpy dtype: {dtype_str}; expected float64, float32, int64, int32, or uint8" + "unsupported numpy dtype: {dtype_str}; expected float64, float32, int64, int32, \ + uint8, complex64 or complex128" ))); } }; diff --git a/crates/clawhdf5-py/tests/test_write_read.py b/crates/clawhdf5-py/tests/test_write_read.py index d5b8672..280496a 100644 --- a/crates/clawhdf5-py/tests/test_write_read.py +++ b/crates/clawhdf5-py/tests/test_write_read.py @@ -247,6 +247,51 @@ def test_roundtrip_uint8(tmp_h5): assert result.dtype == np.uint8 +@pytest.mark.parametrize("dtype", [np.complex64, np.complex128]) +def test_roundtrip_complex(tmp_h5, dtype): + """Complex arrays are written as h5py writes them (a compound {r, i}); + h5py and clawhdf5 both read them back as the same numpy complex dtype.""" + import h5py + + original = (np.arange(12).reshape(3, 4) * (1.5 - 0.25j)).astype(dtype) + with clawhdf5.File(tmp_h5, "w") as f: + f.create_dataset("z", data=original) + f.create_dataset("zc", data=original, chunks=(2, 2), compression="gzip") + with clawhdf5.File(tmp_h5, "r") as f: + for name in ["z", "zc"]: + result = f[name][:] + assert result.dtype == dtype + np.testing.assert_array_equal(result, original) + with h5py.File(tmp_h5, "r") as f: + for name in ["z", "zc"]: + assert f[name].dtype == dtype + assert f[name].id.get_type().get_class() == h5py.h5t.COMPOUND + np.testing.assert_array_equal(f[name][:], original) + + +def test_read_native_complex_from_h5py(tmp_h5): + """HDF5 2.0's native complex type (class 11), written through h5py's + low-level API, reads as numpy complex.""" + import h5py + from h5py import h5s, h5t + + if not getattr(h5py.get_config(), "has_native_complex", False): + pytest.skip("h5py's libhdf5 predates 2.0") + original = np.array([1 + 2j, -3.5 + 0j, 0 - 1e-3j]) + with h5py.File(tmp_h5, "w") as f: + for name, t, dt in [ + (b"n64", h5t.COMPLEX_IEEE_F32LE, np.complex64), + (b"n128", h5t.COMPLEX_IEEE_F64LE, np.complex128), + ]: + d = h5py.h5d.create(f.id, name, t, h5s.create_simple((3,))) + d.write(h5s.ALL, h5s.ALL, original.astype(dt), mtype=t) + with clawhdf5.File(tmp_h5, "r") as f: + for name, dt in [("n64", np.complex64), ("n128", np.complex128)]: + result = f[name][:] + assert result.dtype == dt + np.testing.assert_array_equal(result, original.astype(dt)) + + # --------------------------------------------------------------------------- # Test: chunked + compressed datasets # --------------------------------------------------------------------------- diff --git a/crates/clawhdf5-tools/src/dtype.rs b/crates/clawhdf5-tools/src/dtype.rs index 2c1a1fd..f7cea2e 100644 --- a/crates/clawhdf5-tools/src/dtype.rs +++ b/crates/clawhdf5-tools/src/dtype.rs @@ -135,6 +135,9 @@ pub fn short(dt: &Datatype) -> String { return f.short.into(); } match dt { + Datatype::Complex { size, base_type } => { + short(&Datatype::complex_as_compound(*size, base_type)) + } Datatype::FixedPoint { size, signed, @@ -213,6 +216,9 @@ pub fn long(dt: &Datatype) -> String { return f.long.into(); } match dt { + Datatype::Complex { size, base_type } => { + long(&Datatype::complex_as_compound(*size, base_type)) + } Datatype::FixedPoint { size, signed, @@ -492,6 +498,9 @@ fn string_ddl( /// hdf5-json type object. pub fn json(dt: &Datatype) -> J { match dt { + Datatype::Complex { size, base_type } => { + json(&Datatype::complex_as_compound(*size, base_type)) + } Datatype::FixedPoint { .. } => { json!({"class": "H5T_INTEGER", "base": atomic_ddl(dt)}) } @@ -589,6 +598,7 @@ fn pad_json(p: &StringPadding) -> &'static str { /// Class name used to decide whether two datatypes can be compared. pub fn class(dt: &Datatype) -> &'static str { match dt { + Datatype::Complex { .. } => "compound", Datatype::FixedPoint { .. } => "integer", Datatype::FloatingPoint { .. } => "float", Datatype::Time { .. } => "time", diff --git a/crates/clawhdf5-tools/src/value.rs b/crates/clawhdf5-tools/src/value.rs index c4f4af6..8412e8a 100644 --- a/crates/clawhdf5-tools/src/value.rs +++ b/crates/clawhdf5-tools/src/value.rs @@ -183,6 +183,11 @@ impl<'a> Decoder<'a> { return Value::Error("short element".into()); }; match dt { + Datatype::Complex { size, base_type } => self.decode( + &Datatype::complex_as_compound(*size, base_type), + b, + depth + 1, + ), Datatype::FixedPoint { .. } => match decode_int(dt, b) { Some(v) => Value::Int(v), None => Value::Bytes(b.to_vec()), diff --git a/crates/clawhdf5-wasm/src/core.rs b/crates/clawhdf5-wasm/src/core.rs index 0b5f78e..ea92fb5 100644 --- a/crates/clawhdf5-wasm/src/core.rs +++ b/crates/clawhdf5-wasm/src/core.rs @@ -487,6 +487,9 @@ pub fn describe(dt: &Datatype) -> String { } } match dt { + Datatype::Complex { size, base_type } => { + describe(&Datatype::complex_as_compound(*size, base_type)) + } Datatype::FixedPoint { size, signed, diff --git a/crates/clawhdf5/src/lib.rs b/crates/clawhdf5/src/lib.rs index 1d763d9..24bd435 100644 --- a/crates/clawhdf5/src/lib.rs +++ b/crates/clawhdf5/src/lib.rs @@ -76,7 +76,10 @@ pub use clawhdf5_format::provenance; pub use clawhdf5_format::selection::Selection; pub use clawhdf5_format::storage::Storage; pub use clawhdf5_format::superblock::swmr_flags; -pub use clawhdf5_format::type_builders::{CompoundTypeBuilder, EnumTypeBuilder, FillTime}; +pub use clawhdf5_format::type_builders::{ + CompoundTypeBuilder, EnumTypeBuilder, FillTime, make_complex_f32_type, make_complex_f64_type, + make_native_complex_f32_type, make_native_complex_f64_type, +}; #[cfg(test)] mod tests { diff --git a/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index 7e2ceb5..4456d66 100644 --- a/crates/clawhdf5/src/reader.rs +++ b/crates/clawhdf5/src/reader.rs @@ -1317,6 +1317,60 @@ impl<'f> Dataset<'f> { }) } + /// Read a complex dataset as `[re, im]` pairs of `f64`. Accepts both + /// encodings: h5py's compound `{r, i}` (`with_complex_f64_data`) and + /// HDF5 2.0's native complex type, class 11 + /// (`with_native_complex_f64_data`), over any float base (`f32` and + /// `f16` parts are widened exactly). Any other datatype is a + /// `TypeMismatch` error. + pub fn read_complex_f64(&self) -> Result, Error> { + let (re, im) = self.read_complex_parts(data_read::read_as_f64)?; + Ok(re.into_iter().zip(im).map(|(r, i)| [r, i]).collect()) + } + + /// Read a complex dataset as `[re, im]` pairs of `f32`; see + /// [`read_complex_f64`](Self::read_complex_f64). `f64` parts are + /// rounded to `f32`. + pub fn read_complex_f32(&self) -> Result, Error> { + let (re, im) = self.read_complex_parts(data_read::read_as_f32)?; + Ok(re.into_iter().zip(im).map(|(r, i)| [r, i]).collect()) + } + + fn read_complex_parts( + &self, + convert: fn(&[u8], &Datatype) -> Result, FormatError>, + ) -> Result<(Vec, Vec), Error> { + let dt = self.datatype()?; + let is_complex = match &dt { + // Class 11 parses to this same `{r, i}` compound. + Datatype::Compound { size, members } => { + matches!(members.as_slice(), [r, i] + if r.name == "r" && i.name == "i" + && r.datatype == i.datatype + && matches!(r.datatype, Datatype::FloatingPoint { .. }) + && r.byte_offset == 0 + && i.byte_offset == u64::from(r.datatype.type_size()) + && *size == 2 * r.datatype.type_size()) + } + _ => false, + }; + if !is_complex { + return Err(Error::Format(FormatError::TypeMismatch { + expected: "complex ({r, i} compound or HDF5 2.0 native complex)", + actual: "another datatype", + })); + } + let raw = self.read_raw()?; + let mut fields = data_read::read_compound_fields(&raw, &dt)?.into_iter(); + let (Some(r), Some(i)) = (fields.next(), fields.next()) else { + return Ok((Vec::new(), Vec::new())); + }; + Ok(( + convert(&r.raw_data, &r.datatype)?, + convert(&i.raw_data, &i.datatype)?, + )) + } + /// Read all data as `i32` values. pub fn read_i32(&self) -> Result, Error> { self.file.retry(|| { diff --git a/crates/clawhdf5/src/vlen.rs b/crates/clawhdf5/src/vlen.rs index f471b24..f92f21e 100644 --- a/crates/clawhdf5/src/vlen.rs +++ b/crates/clawhdf5/src/vlen.rs @@ -64,6 +64,7 @@ fn class_name(dt: &Datatype) -> &'static str { } => "variable-length string", Datatype::VariableLength { .. } => "variable-length sequence", Datatype::Array { .. } => "array", + Datatype::Complex { .. } => "complex", } } diff --git a/crates/clawhdf5/tests/h5py_interop_tests.rs b/crates/clawhdf5/tests/h5py_interop_tests.rs index 042940f..c0e5b1b 100644 --- a/crates/clawhdf5/tests/h5py_interop_tests.rs +++ b/crates/clawhdf5/tests/h5py_interop_tests.rs @@ -1550,3 +1550,51 @@ print("ok") ); assert_eq!(run_python_output(&script), "ok"); } + +/// h5py writes complex numbers two ways: its default compound `{r, i}` (the +/// only form its high-level API writes, h5py 3.16 included) and, through the +/// low-level API on libhdf5 2.0+, the native complex type (class 11). Both +/// read back through `read_complex_f32`/`read_complex_f64`. +#[test] +fn h5py_complex_datasets_read_as_complex() { + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("complex.h5"); + let path_str = path.display().to_string(); + let native = run_python_output(&format!( + r#" +import h5py, numpy as np +from h5py import h5t, h5s, h5d, h5p +z = np.arange(24).reshape(4, 6) * (1 - 0.5j) +with h5py.File("{path_str}", "w") as f: + f["c64"] = z.astype(np.complex64) + f.create_dataset("c128", data=z, chunks=(2, 4), compression="gzip") + native = getattr(h5py.get_config(), "has_native_complex", False) + if native: + for name, t, dt in [(b"n64", h5t.COMPLEX_IEEE_F32LE, np.complex64), + (b"n128", h5t.COMPLEX_IEEE_F64LE, np.complex128)]: + dcpl = h5p.create(h5p.DATASET_CREATE) + dcpl.set_chunk((2, 4)) + dcpl.set_deflate(4) + d = h5d.create(f.id, name, t, h5s.create_simple((4, 6)), dcpl=dcpl) + d.write(h5s.ALL, h5s.ALL, np.ascontiguousarray(z.astype(dt)), mtype=t) +print(int(native)) +"# + )); + let want: Vec<[f64; 2]> = (0..24).map(|k| [k as f64, -(k as f64) / 2.0]).collect(); + let want32: Vec<[f32; 2]> = want.iter().map(|p| [p[0] as f32, p[1] as f32]).collect(); + + let file = File::open(&path).unwrap(); + let mut names = vec!["c64", "c128"]; + if native == "1" { + names.extend(["n64", "n128"]); + } else { + eprintln!("HDF5 < 2.0: native complex not checked"); + } + for name in names { + let ds = file.dataset(name).unwrap(); + assert_eq!(ds.shape().unwrap(), vec![4, 6], "{name}"); + assert_eq!(ds.read_complex_f64().unwrap(), want, "{name}"); + assert_eq!(ds.read_complex_f32().unwrap(), want32, "{name}"); + } +} diff --git a/crates/clawhdf5/tests/integration_tests.rs b/crates/clawhdf5/tests/integration_tests.rs index 47798e8..f978a3a 100644 --- a/crates/clawhdf5/tests/integration_tests.rs +++ b/crates/clawhdf5/tests/integration_tests.rs @@ -1018,3 +1018,97 @@ fn dataset_at_opens_the_same_dataset_as_its_path() { Err(clawhdf5::Error::NotADataset(_)) )); } + +// --------------------------------------------------------------------------- +// Complex numbers: h5py's compound {r, i} and HDF5 2.0 native (class 11) +// --------------------------------------------------------------------------- + +#[test] +fn complex_datasets_and_attributes_round_trip() { + let z64 = [[1.5f32, -2.0], [0.0, 3.25], [f32::MAX, f32::MIN_POSITIVE]]; + let z128 = [[1.0f64, 2.0], [-3.5, 4.0e300], [f64::NAN, -0.0]]; + let big: Vec<[f64; 2]> = (0..1000).map(|k| [k as f64, -(k as f64) / 3.0]).collect(); + + let mut b = FileBuilder::new(); + b.create_dataset("compound64").with_complex_f32_data(&z64); + b.create_dataset("compound128").with_complex_f64_data(&z128); + b.create_dataset("native64") + .with_native_complex_f32_data(&z64) + .set_attr( + "c", + AttrValue::Raw { + datatype: clawhdf5::make_native_complex_f64_type(), + shape: vec![], + data: [0.5f64.to_le_bytes(), (-1.5f64).to_le_bytes()].concat(), + }, + ); + b.create_dataset("native128") + .with_native_complex_f64_data(&z128); + b.create_dataset("native_chunked") + .with_native_complex_f64_data(&big) + .with_shape(&[10, 100]) + .with_chunks(&[5, 30]) + .with_deflate(4); + b.create_dataset("empty").with_native_complex_f64_data(&[]); + let file = File::from_bytes(b.finish().unwrap()).unwrap(); + + let same = |got: Vec<[f64; 2]>, want: &[[f64; 2]]| { + assert_eq!(got.len(), want.len()); + for (g, w) in got.iter().zip(want) { + for k in 0..2 { + assert!(g[k].to_bits() == w[k].to_bits(), "{g:?} != {w:?}"); + } + } + }; + for name in ["compound64", "native64"] { + let ds = file.dataset(name).unwrap(); + assert_eq!(ds.read_complex_f32().unwrap(), z64, "{name}"); + assert_eq!(ds.shape().unwrap(), vec![3]); + // Both encodings read as the same {r, i} compound. + assert_eq!( + ds.dtype().unwrap(), + DType::Compound(vec![("r".into(), DType::F32), ("i".into(), DType::F32)]) + ); + } + for name in ["compound128", "native128"] { + same( + file.dataset(name).unwrap().read_complex_f64().unwrap(), + &z128, + ); + } + let chunked = file.dataset("native_chunked").unwrap(); + assert_eq!(chunked.shape().unwrap(), vec![10, 100]); + same(chunked.read_complex_f64().unwrap(), &big); + assert!( + file.dataset("empty") + .unwrap() + .read_complex_f64() + .unwrap() + .is_empty() + ); + + match file.dataset("native64").unwrap().attr("c").unwrap() { + Some(AttrValue::Raw { + datatype, + shape, + data, + }) => { + assert!(shape.is_empty()); + assert_eq!(datatype.type_size(), 16); + let fields = + clawhdf5_format::data_read::read_compound_fields(&data, &datatype).unwrap(); + let part = |k: usize| { + clawhdf5_format::data_read::read_as_f64(&fields[k].raw_data, &fields[k].datatype) + .unwrap() + }; + assert_eq!((part(0), part(1)), (vec![0.5], vec![-1.5])); + } + other => panic!("expected a raw complex attribute, got {other:?}"), + } + + // Not complex: a clear error, not garbage. + let mut b = FileBuilder::new(); + b.create_dataset("x").with_f64_data(&[1.0, 2.0]); + let file = File::from_bytes(b.finish().unwrap()).unwrap(); + assert!(file.dataset("x").unwrap().read_complex_f64().is_err()); +} diff --git a/docs/QUICKSTART.md b/docs/QUICKSTART.md index f634a85..a288613 100644 --- a/docs/QUICKSTART.md +++ b/docs/QUICKSTART.md @@ -218,7 +218,8 @@ with clawhdf5.File("data.h5", "r+") as f: `'r+'` cannot create or delete datasets and groups, or delete attributes (`NotImplementedError`, nothing written). New files (`'w'`) take numeric -arrays (`float64`, `float32`, `int64`, `int32`, `uint8`): +arrays (`float64`, `float32`, `int64`, `int32`, `uint8`, and `complex64`/ +`complex128`, stored as h5py's compound `{r, i}`): ```python with clawhdf5.File("new.h5", "w") as f: diff --git a/docs/known-issues.md b/docs/known-issues.md index 947c356..46446d0 100644 --- a/docs/known-issues.md +++ b/docs/known-issues.md @@ -174,6 +174,15 @@ wrong data. decoded, and external references are an error (object references decode). A multi-dimensional numeric attribute is returned as a flat array (its shape is not reported; `AttrValue::Raw` carries the shape). + HDF5 2.0's native complex type (class 11) is read as the equivalent + `{r, i}` compound: values are right (`Dataset::read_complex_f64`, and + numpy complex in Python), but `raw_datatype()`, `dtype()`, `h5rs + ls`/`dump` and the browser reader show a compound where h5dump 2.x + prints `H5T_COMPLEX_IEEE_F64LE`, so `h5rs dump` of such a file does not + match h5dump 2.x (h5dump 1.14 cannot read it at all). Writing class 11 + (added 2026-09-28) is opt-in (`with_native_complex_f64_data`, + `make_native_complex_f64_type`); the Python bindings write complex + arrays as h5py's compound only. - **Metadata cache images** (read since 2026-09-26) differ from libhdf5 in that: libhdf5 fails only the first metadata read of an image it cannot load and then reads the file's own (possibly stale) metadata, where we