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) <[email protected]>
This commit is contained in:
osobh
2026-09-28 21:25:03 -05:00
co-authored by Claude Opus 5.5
parent bce07e9cb9
commit d0d5347cd9
21 changed files with 790 additions and 24 deletions
+4
View File
@@ -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<Vec<CompoundFieldData>, 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 {
+176 -17
View File
@@ -140,6 +140,19 @@ pub enum Datatype {
base_type: Box<Datatype>,
dimensions: Vec<u32>,
},
/// 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<Datatype> },
}
/// 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);
@@ -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<T: Copy, const N: usize>(
&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,
@@ -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<u8> = [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() {