Merge branch 'fix/p0-reader-numeric' into fix/phase0-correctness

crates/clawhdf5-format/src/data_read.rs

@@ -773,14 +773,7 @@ pub fn read_as_f64_zerocopy<'a>(raw: &'a [u8], datatype: &Datatype) -> Option<&'
     // Only native LE f64 is eligible
     #[cfg(target_endian = "little")]
     {
-        if !matches!(
+        if !is_native_le_float(datatype, FloatFormat::Double) {
-            datatype,
-            Datatype::FloatingPoint {
-                size: 8,
-                byte_order: DatatypeByteOrder::LittleEndian,
-                ..
-            }
-        ) {
             return None;
         }
         if !raw.len().is_multiple_of(8) {
@@ -809,14 +802,7 @@ pub fn read_as_f64_zerocopy<'a>(raw: &'a [u8], datatype: &Datatype) -> Option<&'
 pub fn read_as_f32_zerocopy<'a>(raw: &'a [u8], datatype: &Datatype) -> Option<&'a [f32]> {
     #[cfg(target_endian = "little")]
     {
-        if !matches!(
+        if !is_native_le_float(datatype, FloatFormat::Single) {
-            datatype,
-            Datatype::FloatingPoint {
-                size: 4,
-                byte_order: DatatypeByteOrder::LittleEndian,
-                ..
-            }
-        ) {
             return None;
         }
         if !raw.len().is_multiple_of(4) {
@@ -902,9 +888,9 @@ fn native_le_to_vec<T: Copy>(raw: &[u8], count: usize) -> Vec<T> {
 
 /// Convert raw bytes to `f64` values.
 pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result<Vec<f64>, FormatError> {
-    // Array datatypes (e.g. an array-typed compound member) are read as a flat
+    // Array datatypes read as a flat sequence of their base elements, and
-    // sequence of their base elements.
+    // enumerations (h5py's bool among them) as their integer values.
-    if let Datatype::Array { base_type, .. } = datatype {
+    if let Datatype::Array { base_type, .. } | Datatype::Enumeration { base_type, .. } = datatype {
         return read_as_f64(raw, base_type);
     }
     ensure_numeric(datatype, "FloatingPoint or FixedPoint")?;
@@ -919,20 +905,19 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result<Vec<f64>, FormatEr
 
     // Fast path: native-endian f64 — single bulk memcpy
     #[cfg(target_endian = "little")]
-    if matches!(
+    if is_native_le_float(datatype, FloatFormat::Double) {
-        datatype,
-        Datatype::FloatingPoint {
-            size: 8,
-            byte_order: DatatypeByteOrder::LittleEndian,
-            ..
-        }
-    ) {
         return Ok(native_le_to_vec::<f64>(raw, count));
     }
 
     let order = get_byte_order(datatype);
     let mut result = Vec::with_capacity(count);
-
+    if let Datatype::FloatingPoint { .. } = datatype {
+        let format = FloatFormat::of(datatype)?;
+        for chunk in raw.chunks_exact(elem_size) {
+            result.push(format.decode(chunk, &order));
+        }
+        return Ok(result);
+    }
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
         let val = convert_to_f64(chunk, datatype, &order)?;
@@ -947,18 +932,7 @@ fn convert_to_f64(
     order: &DatatypeByteOrder,
 ) -> Result<f64, FormatError> {
     match dt {
-        Datatype::FloatingPoint { size, .. } => match size {
+        Datatype::FloatingPoint { .. } => Ok(FloatFormat::of(dt)?.decode(bytes, order)),
-            4 => {
-                let v = read_f32_bytes(bytes, order);
-                Ok(v as f64)
-            }
-            8 => Ok(read_f64_bytes(bytes, order)),
-            2 => Ok(read_f16_bytes(bytes, order) as f64),
-            _ => Err(FormatError::DataSizeMismatch {
-                expected: 8,
-                actual: *size as usize,
-            }),
-        },
         Datatype::FixedPoint {
             size,
             signed,
@@ -982,9 +956,83 @@ fn convert_to_f64(
     }
 }
 
+/// One numeric element as stored, before conversion to the caller's type.
+#[derive(Debug, Clone, Copy, PartialEq)]
+enum Scalar {
+    Signed(i64),
+    Unsigned(u64),
+    Float(f64),
+}
+
+impl Scalar {
+    // Every conversion follows libhdf5's default (hard) conversions: a value
+    // outside the target type's range saturates to its minimum or maximum —
+    // including a negative value read as unsigned, which reads as 0 — rather
+    // than being truncated to its low bits. Floats truncate toward zero; NaN
+    // converts to 0 (libhdf5 leaves that case to the C cast, whose result is
+    // platform-dependent).
+
+    fn to_i64(self) -> i64 {
+        match self {
+            Scalar::Signed(v) => v,
+            Scalar::Unsigned(v) => i64::try_from(v).unwrap_or(i64::MAX),
+            Scalar::Float(v) => v as i64,
+        }
+    }
+
+    fn to_u64(self) -> u64 {
+        match self {
+            Scalar::Signed(v) => u64::try_from(v).unwrap_or(0),
+            Scalar::Unsigned(v) => v,
+            Scalar::Float(v) => v as u64,
+        }
+    }
+
+    fn to_i32(self) -> i32 {
+        match self {
+            Scalar::Signed(v) => v.clamp(i32::MIN.into(), i32::MAX.into()) as i32,
+            Scalar::Unsigned(v) => i32::try_from(v).unwrap_or(i32::MAX),
+            Scalar::Float(v) => v as i32,
+        }
+    }
+}
+
+/// Decode one element of a numeric datatype.
+fn decode_scalar(
+    bytes: &[u8],
+    dt: &Datatype,
+    order: &DatatypeByteOrder,
+) -> Result<Scalar, FormatError> {
+    match dt {
+        Datatype::FixedPoint {
+            size,
+            signed,
+            bit_offset,
+            bit_precision,
+            ..
+        } => {
+            let full = read_unsigned_int(bytes, *size as usize, order);
+            let (off, prec) = effective_bits(*size as usize, *bit_offset, *bit_precision);
+            Ok(if *signed {
+                Scalar::Signed(extract_signed(full, off, prec))
+            } else {
+                Scalar::Unsigned(extract_unsigned(full, off, prec))
+            })
+        }
+        _ => convert_to_f64(bytes, dt, order).map(Scalar::Float),
+    }
+}
+
 /// Convert raw bytes to `i64` values.
+///
+/// Values are converted the way libhdf5 converts them: integers outside the
+/// target range saturate at its minimum or maximum (a negative value read as
+/// unsigned is 0), and floating-point data is truncated toward zero and
+/// saturated, with NaN read as 0.
 pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result<Vec<i64>, FormatError> {
-    if let Datatype::Array { base_type, .. } = datatype {
+    // Array datatypes read as a flat sequence of their base elements, and
+    // enumerations (h5py's bool among them) as their integer values.
+    if let Datatype::Array { base_type, .. } | Datatype::Enumeration { base_type, .. } = datatype {
         return read_as_i64(raw, base_type);
     }
     ensure_numeric(datatype, "FixedPoint (signed)")?;
@@ -1014,19 +1062,24 @@ pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result<Vec<i64>, FormatEr
     }
 
     let order = get_byte_order(datatype);
-    let (off, prec) = fixed_bits(datatype);
     let mut result = Vec::with_capacity(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
-        let full = read_unsigned_int(chunk, elem_size, &order);
+        result.push(decode_scalar(chunk, datatype, &order)?.to_i64());
-        result.push(extract_signed(full, off, prec));
     }
     Ok(result)
 }
 
 /// Convert raw bytes to `u64` values.
+///
+/// Values are converted the way libhdf5 converts them: integers outside the
+/// target range saturate at its minimum or maximum (a negative value read as
+/// unsigned is 0), and floating-point data is truncated toward zero and
+/// saturated, with NaN read as 0.
 pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result<Vec<u64>, FormatError> {
-    if let Datatype::Array { base_type, .. } = datatype {
+    // Array datatypes read as a flat sequence of their base elements, and
+    // enumerations (h5py's bool among them) as their integer values.
+    if let Datatype::Array { base_type, .. } | Datatype::Enumeration { base_type, .. } = datatype {
         return read_as_u64(raw, base_type);
     }
     ensure_numeric(datatype, "FixedPoint (unsigned)")?;
@@ -1039,19 +1092,19 @@ pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result<Vec<u64>, FormatEr
     }
     let count = raw.len() / elem_size;
     let order = get_byte_order(datatype);
-    let (off, prec) = fixed_bits(datatype);
     let mut result = Vec::with_capacity(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
-        let full = read_unsigned_int(chunk, elem_size, &order);
+        result.push(decode_scalar(chunk, datatype, &order)?.to_u64());
-        result.push(extract_unsigned(full, off, prec));
     }
     Ok(result)
 }
 
 /// Convert raw bytes to `f32` values.
 pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatError> {
-    if let Datatype::Array { base_type, .. } = datatype {
+    // Array datatypes read as a flat sequence of their base elements, and
+    // enumerations (h5py's bool among them) as their integer values.
+    if let Datatype::Array { base_type, .. } | Datatype::Enumeration { base_type, .. } = datatype {
         return read_as_f32(raw, base_type);
     }
     ensure_numeric(datatype, "FloatingPoint")?;
@@ -1066,25 +1119,11 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatEr
 
     // Fast path: native-endian f32 — single bulk memcpy
     #[cfg(target_endian = "little")]
-    if matches!(
+    if is_native_le_float(datatype, FloatFormat::Single) {
-        datatype,
-        Datatype::FloatingPoint {
-            size: 4,
-            byte_order: DatatypeByteOrder::LittleEndian,
-            ..
-        }
-    ) {
         return Ok(native_le_to_vec::<f32>(raw, count));
     }
-    // Little-endian half precision (numpy float16): widen directly.
+    // Little-endian IEEE half precision (numpy float16): widen directly.
-    if matches!(
+    if is_native_le_float(datatype, FloatFormat::Half) {
-        datatype,
-        Datatype::FloatingPoint {
-            size: 2,
-            byte_order: DatatypeByteOrder::LittleEndian,
-            ..
-        }
-    ) {
         let (halves, _) = raw[..count * 2].as_chunks::<2>();
         return Ok(halves
             .iter()
@@ -1094,18 +1133,22 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatEr
 
     let order = get_byte_order(datatype);
     let mut result = Vec::with_capacity(count);
+    if let Datatype::FloatingPoint { .. } = datatype {
+        let format = FloatFormat::of(datatype)?;
+        for chunk in raw.chunks_exact(elem_size) {
+            result.push(match format {
+                FloatFormat::Single => read_f32_bytes(chunk, &order),
+                FloatFormat::Half => read_f16_bytes(chunk, &order),
+                // Double rounds; every other supported layout (bfloat16, FP8)
+                // is exact in f32.
+                _ => format.decode(chunk, &order) as f32,
+            });
+        }
+        return Ok(result);
+    }
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
         match datatype {
-            Datatype::FloatingPoint { size: 4, .. } => {
-                result.push(read_f32_bytes(chunk, &order));
-            }
-            Datatype::FloatingPoint { size: 8, .. } => {
-                result.push(read_f64_bytes(chunk, &order) as f32);
-            }
-            Datatype::FloatingPoint { size: 2, .. } => {
-                result.push(read_f16_bytes(chunk, &order));
-            }
             Datatype::FixedPoint {
                 signed: true,
                 size,
@@ -1140,8 +1183,15 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatEr
 }
 
 /// Convert raw bytes to `i32` values.
+///
+/// Values are converted the way libhdf5 converts them: integers outside the
+/// target range saturate at its minimum or maximum (a negative value read as
+/// unsigned is 0), and floating-point data is truncated toward zero and
+/// saturated, with NaN read as 0.
 pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result<Vec<i32>, FormatError> {
-    if let Datatype::Array { base_type, .. } = datatype {
+    // Array datatypes read as a flat sequence of their base elements, and
+    // enumerations (h5py's bool among them) as their integer values.
+    if let Datatype::Array { base_type, .. } | Datatype::Enumeration { base_type, .. } = datatype {
         return read_as_i32(raw, base_type);
     }
     ensure_numeric(datatype, "FixedPoint")?;
@@ -1162,6 +1212,7 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result<Vec<i32>, FormatEr
             datatype,
             Datatype::FixedPoint {
                 byte_order: DatatypeByteOrder::LittleEndian,
+                signed: true,
                 ..
             }
         )
@@ -1170,12 +1221,10 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result<Vec<i32>, FormatEr
     }
 
     let order = get_byte_order(datatype);
-    let (off, prec) = fixed_bits(datatype);
     let mut result = Vec::with_capacity(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
-        let full = read_unsigned_int(chunk, elem_size, &order);
+        result.push(decode_scalar(chunk, datatype, &order)?.to_i32());
-        result.push(extract_signed(full, off, prec) as i32);
     }
     Ok(result)
 }
@@ -1616,6 +1665,174 @@ fn reorder_bytes(bytes: &[u8], order: &DatatypeByteOrder) -> [u8; 8] {
     buf
 }
 
+/// How the bits of a floating-point datatype are laid out, read from the
+/// datatype message's fields rather than assumed from its size (a 2-byte
+/// float may be IEEE half or bfloat16).
+#[derive(Debug, Clone, Copy, PartialEq)]
+enum FloatFormat {
+    /// IEEE-754 binary16.
+    Half,
+    /// IEEE-754 binary32.
+    Single,
+    /// IEEE-754 binary64.
+    Double,
+    /// Any other IEEE-style layout (implied leading mantissa bit, all-ones
+    /// exponent for infinity/NaN) whose values are all exact in `f64`:
+    /// bfloat16, the FP8 formats, and similar.
+    Other(FloatLayout),
+}
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+struct FloatLayout {
+    exponent_location: u32,
+    exponent_size: u32,
+    mantissa_location: u32,
+    mantissa_size: u32,
+    exponent_bias: u32,
+}
+
+impl FloatFormat {
+    fn of(dt: &Datatype) -> Result<FloatFormat, FormatError> {
+        let Datatype::FloatingPoint {
+            size,
+            exponent_location,
+            exponent_size,
+            mantissa_location,
+            mantissa_size,
+            exponent_bias,
+            ..
+        } = dt
+        else {
+            return Err(FormatError::TypeMismatch {
+                expected: "FloatingPoint",
+                actual: datatype_name(dt),
+            });
+        };
+        let layout = FloatLayout {
+            exponent_location: u32::from(*exponent_location),
+            exponent_size: u32::from(*exponent_size),
+            mantissa_location: u32::from(*mantissa_location),
+            mantissa_size: u32::from(*mantissa_size),
+            exponent_bias: *exponent_bias,
+        };
+        let fields = (
+            layout.exponent_location,
+            layout.exponent_size,
+            layout.mantissa_location,
+            layout.mantissa_size,
+            layout.exponent_bias,
+        );
+        let bits = size.saturating_mul(8);
+        // The sign bit is not kept in `Datatype`; every standard layout has it
+        // directly above the exponent, with the mantissa below.
+        let well_formed = layout.exponent_size > 0
+            && layout.mantissa_size > 0
+            && layout.mantissa_location + layout.mantissa_size <= layout.exponent_location
+            && layout.exponent_location + layout.exponent_size < bits;
+        match (size, fields) {
+            (2, (10, 5, 0, 10, 15)) => Ok(FloatFormat::Half),
+            (4, (23, 8, 0, 23, 127)) => Ok(FloatFormat::Single),
+            (8, (52, 11, 0, 52, 1023)) => Ok(FloatFormat::Double),
+            _ if well_formed
+                && *size <= 8
+                && layout.exponent_size <= 11
+                && layout.mantissa_size <= 52 =>
+            {
+                Ok(FloatFormat::Other(layout))
+            }
+            // Fields that cannot describe any float (e.g. left zeroed by a
+            // hand-built datatype): fall back to the IEEE type of that size.
+            (2, _) if !well_formed => Ok(FloatFormat::Half),
+            (4, _) if !well_formed => Ok(FloatFormat::Single),
+            (8, _) if !well_formed => Ok(FloatFormat::Double),
+            // x87 80-bit extended, binary128, ...: not representable in f64.
+            _ => Err(FormatError::TypeMismatch {
+                expected: "floating point of at most 64 bits (IEEE-style layout)",
+                actual: "FloatingPoint",
+            }),
+        }
+    }
+
+    fn decode(self, bytes: &[u8], order: &DatatypeByteOrder) -> f64 {
+        match self {
+            FloatFormat::Half => f64::from(read_f16_bytes(bytes, order)),
+            FloatFormat::Single => f64::from(read_f32_bytes(bytes, order)),
+            FloatFormat::Double => read_f64_bytes(bytes, order),
+            FloatFormat::Other(layout) => {
+                layout.decode(read_unsigned_int(bytes, bytes.len(), order))
+            }
+        }
+    }
+}
+
+impl FloatLayout {
+    /// Decode the value held in the low `size * 8` bits of `bits`.
+    fn decode(self, bits: u64) -> f64 {
+        let field = |location: u32, size: u32| (bits >> location) & ((1u64 << size) - 1);
+        let exponent = field(self.exponent_location, self.exponent_size);
+        let mantissa = field(self.mantissa_location, self.mantissa_size);
+        let negative = field(self.exponent_location + self.exponent_size, 1) == 1;
+        let max_exponent = (1u64 << self.exponent_size) - 1;
+        let magnitude = if exponent == max_exponent {
+            if mantissa == 0 {
+                f64::INFINITY
+            } else {
+                f64::NAN
+            }
+        } else {
+            let bias = i64::from(self.exponent_bias);
+            let msize = i64::from(self.mantissa_size);
+            // value = significand * 2^power, with an implied leading 1 unless
+            // the number is subnormal (exponent field 0).
+            let (significand, power) = if exponent == 0 {
+                (mantissa, 1 - bias - msize)
+            } else {
+                (
+                    mantissa | (1u64 << self.mantissa_size),
+                    exponent as i64 - bias - msize,
+                )
+            };
+            scale_by_pow2(significand as f64, power)
+        };
+        if negative { -magnitude } else { magnitude }
+    }
+}
+
+/// `x * 2^power` without `std` (no `powi`/`libm`). `x` is a non-negative
+/// integer below 2^53, so it is exact.
+fn scale_by_pow2(x: f64, power: i64) -> f64 {
+    if x == 0.0 || power < -1200 {
+        return 0.0;
+    }
+    if power > 1100 {
+        return f64::INFINITY;
+    }
+    let pow2 = |p: i64| f64::from_bits(((p + 1023) as u64) << 52);
+    let mut x = x;
+    let mut power = power;
+    while power > 1023 {
+        x *= pow2(1023);
+        power -= 1023;
+    }
+    while power < -1022 {
+        x *= pow2(-1022);
+        power += 1022;
+    }
+    x * pow2(power)
+}
+
+/// Whether `datatype` is the little-endian IEEE float `format`, whose bytes
+/// can be copied straight into native values on a little-endian target.
+fn is_native_le_float(datatype: &Datatype, format: FloatFormat) -> bool {
+    matches!(
+        datatype,
+        Datatype::FloatingPoint {
+            byte_order: DatatypeByteOrder::LittleEndian,
+            ..
+        }
+    ) && FloatFormat::of(datatype).is_ok_and(|f| f == format)
+}
+
 fn read_f64_bytes(bytes: &[u8], order: &DatatypeByteOrder) -> f64 {
     let buf = reorder_bytes(bytes, order);
     f64::from_le_bytes(buf)
@@ -1666,20 +1883,6 @@ fn effective_bits(size: usize, bit_offset: u16, bit_precision: u16) -> (u32, u32
     (bit_offset as u32, prec)
 }
 
-/// `(bit_offset, bit_precision)` for a fixed-point datatype, full width for
-/// other types.
-fn fixed_bits(datatype: &Datatype) -> (u32, u32) {
-    match datatype {
-        Datatype::FixedPoint {
-            size,
-            bit_offset,
-            bit_precision,
-            ..
-        } => effective_bits(*size as usize, *bit_offset, *bit_precision),
-        _ => (0, 0),
-    }
-}
-
 /// Whether a datatype occupies its full storage width (bit offset 0, precision
 /// == size·8), in which case the bulk-copy fast read paths apply. Non
 /// fixed-point types are treated as full width.
@@ -1892,6 +2095,67 @@ mod tests {
         assert_eq!(read_as_u64(&raw, &dt).unwrap(), vec![4095, 1, 2048]);
     }
 
+    #[test]
+    fn float_to_int_truncates_and_saturates() {
+        // Values libhdf5 hands to an undefined C cast: NaN reads as 0 and
+        // exactly 2^63 saturates instead of wrapping to i64::MIN.
+        let dt = make_f64_le_type();
+        let vals = [f64::NAN, 2f64.powi(63), -2.5, 2.0f64.powi(64)];
+        let raw: Vec<u8> = vals.iter().flat_map(|v| v.to_le_bytes()).collect();
+        assert_eq!(
+            read_as_i64(&raw, &dt).unwrap(),
+            vec![0, i64::MAX, -2, i64::MAX]
+        );
+        assert_eq!(
+            read_as_u64(&raw, &dt).unwrap(),
+            vec![0, 1 << 63, 0, u64::MAX]
+        );
+        assert_eq!(
+            read_as_i32(&raw, &dt).unwrap(),
+            vec![0, i32::MAX, -2, i32::MAX]
+        );
+    }
+
+    #[test]
+    fn bfloat16_and_fp8_decode_by_fields() {
+        // bfloat16 is a 2-byte float that is not IEEE half.
+        let bf16 = Datatype::FloatingPoint {
+            size: 2,
+            byte_order: DatatypeByteOrder::LittleEndian,
+            bit_offset: 0,
+            bit_precision: 16,
+            exponent_location: 7,
+            exponent_size: 8,
+            mantissa_location: 0,
+            mantissa_size: 7,
+            exponent_bias: 127,
+        };
+        let raw: Vec<u8> = [0x3FC0u16, 0xC010, 0x7F80, 0x0001]
+            .iter()
+            .flat_map(|v| v.to_le_bytes())
+            .collect();
+        let got = read_as_f64(&raw, &bf16).unwrap();
+        assert_eq!(&got[..3], &[1.5, -2.25, f64::INFINITY]);
+        assert_eq!(got[3], 2f64.powi(-133)); // smallest subnormal
+        assert_eq!(read_as_f32(&raw, &bf16).unwrap()[..2], [1.5, -2.25]);
+
+        // FP8 E4M3: 1, -1, 2, 0, NaN (IEEE-style, as libhdf5 treats it).
+        let e4m3 = Datatype::FloatingPoint {
+            size: 1,
+            byte_order: DatatypeByteOrder::LittleEndian,
+            bit_offset: 0,
+            bit_precision: 8,
+            exponent_location: 3,
+            exponent_size: 4,
+            mantissa_location: 0,
+            mantissa_size: 3,
+            exponent_bias: 7,
+        };
+        let got = read_as_f64(&[0x38, 0xB8, 0x40, 0x00, 0x7E], &e4m3).unwrap();
+        assert_eq!(&got[..4], &[1.0, -1.0, 2.0, 0.0]);
+        assert!(got[4].is_nan());
+    }
+
     #[test]
     fn full_width_signed_unchanged() {
         // Regression: full-width 32-bit signed must be unaffected.

crates/clawhdf5-format/src/vl_data.rs

@@ -148,7 +148,12 @@ pub fn read_vl_strings(
     Ok(result)
 }
 
-/// Resolve VL byte sequences from raw data.
+/// Resolve VL sequences from raw data, returning each element's bytes.
+///
+/// Each element is the sequence's full encoding — element count × base type
+/// size bytes, in the base type's byte order — so a sequence of `i32` yields
+/// four bytes per value. Decode it with the base type (e.g.
+/// [`crate::data_read::read_as_i64`]).
 pub fn read_vl_bytes(
     file_data: &[u8],
     raw_data: &[u8],
@@ -177,8 +182,10 @@ pub fn read_vl_bytes(
             },
         )?;
 
-        let len = (vl.length as usize).min(obj.data.len());
+        // The heap object holds the whole sequence. `vl.length` counts
-        result.push(obj.data[..len].to_vec());
+        // elements, not bytes, so it is only the byte length when the base
+        // type is one byte wide.
+        result.push(obj.data.clone());
     }
 
     Ok(result)

crates/clawhdf5/tests/numeric_conversion_interop.rs

@@ -0,0 +1,374 @@
+//! Numeric conversions on read, checked against h5py/libhdf5.
+//!
+//! h5py writes each file and prints what libhdf5 converts the data to
+//! (`Dataset.astype`); the typed readers must return the same values.
+//! Skipped when python3 with h5py is unavailable, unless
+//! `CLAWHDF5_REQUIRE_INTEROP=1`.
+
+use std::collections::HashMap;
+use std::process::Command;
+
+use clawhdf5::File;
+
+fn python() -> String {
+    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
+}
+
+fn interop_required() -> bool {
+    std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1")
+}
+
+fn python_available() -> bool {
+    Command::new(python())
+        .args(["-c", "import h5py, numpy"])
+        .output()
+        .map(|o| o.status.success())
+        .unwrap_or(false)
+}
+
+macro_rules! skip_if_no_python {
+    () => {
+        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;
+        }
+    };
+}
+
+/// Run `script` (which writes the file at `path`) and return its stdout as
+/// `key -> values`, one `key v1 v2 ...` line per key.
+fn run_python(script: &str) -> HashMap<String, Vec<String>> {
+    let output = Command::new(python())
+        .args(["-c", script])
+        .output()
+        .expect("failed to run python");
+    assert!(
+        output.status.success(),
+        "python failed:\n{}",
+        String::from_utf8_lossy(&output.stderr)
+    );
+    String::from_utf8_lossy(&output.stdout)
+        .lines()
+        .filter_map(|line| {
+            let mut words = line.split_whitespace().map(str::to_string);
+            Some((words.next()?, words.collect()))
+        })
+        .collect()
+}
+
+fn parse<T: std::str::FromStr>(values: &[String]) -> Vec<T>
+where
+    T::Err: std::fmt::Debug,
+{
+    values.iter().map(|v| v.parse().unwrap()).collect()
+}
+
+/// Python prelude: `emit(key, array)` prints one line of integers.
+const PRELUDE: &str = r#"
+import h5py, numpy as np
+def emit(key, arr):
+    print(key, *[int(v) for v in np.asarray(arr).ravel()])
+"#;
+
+#[test]
+fn float_dataset_read_as_integers_converts_like_libhdf5() {
+    // read_i32/read_i64/read_u64 on a float dataset used to return the raw
+    // IEEE bit patterns (1.5 read as i64 was 4609434218613702656).
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("float_to_int.h5");
+    let script = format!(
+        r#"{PRELUDE}
+vals = [1.5, -2.75, 3e9, 1e300, -1e300, -0.5, 0.0, 7.99, np.inf, -np.inf, 1e19, -1e19]
+with h5py.File("{path}", "w") as f:
+    for name, dt in (("f8", "<f8"), ("f8be", ">f8"), ("f4", "<f4")):
+        f.create_dataset(name, data=np.array(vals).astype(dt))
+    # libhdf5's half conversions are not saturating (an infinite half becomes
+    # INT_MIN whatever its sign, a negative one wraps as u64), so the half
+    # case stays finite and its u64 read is checked separately below.
+    f.create_dataset("f2", data=np.array([1.5, -2.75, -0.5, 0.0, 7.99, 65504, -65504], "<f2"))
+with h5py.File("{path}", "r") as f:
+    for name in ("f8", "f8be", "f4", "f2"):
+        d = f[name]
+        emit(name + ":i32", d.astype("<i4")[()])
+        emit(name + ":i64", d.astype("<i8")[()])
+        emit(name + ":u64", d.astype("<u8")[()])
+"#,
+        path = path.display()
+    );
+    let expected = run_python(&script);
+
+    let file = File::open(&path).unwrap();
+    for name in ["f8", "f8be", "f4", "f2"] {
+        let ds = file.dataset(name).unwrap();
+        assert_eq!(
+            ds.read_i32().unwrap(),
+            parse::<i32>(&expected[&format!("{name}:i32")]),
+            "{name} as i32"
+        );
+        assert_eq!(
+            ds.read_i64().unwrap(),
+            parse::<i64>(&expected[&format!("{name}:i64")]),
+            "{name} as i64"
+        );
+        if name != "f2" {
+            assert_eq!(
+                ds.read_u64().unwrap(),
+                parse::<u64>(&expected[&format!("{name}:u64")]),
+                "{name} as u64"
+            );
+        }
+    }
+    // Negative values saturate at 0 rather than wrapping.
+    let f2 = file.dataset("f2").unwrap();
+    assert_eq!(f2.read_u64().unwrap(), vec![1, 0, 0, 0, 7, 65504, 0]);
+}
+
+#[test]
+fn integer_reads_saturate_out_of_range_values_like_libhdf5() {
+    // Narrowing reads used to keep the low bits (i64 2^40+5 read as i32 was
+    // 5, u64::MAX read as i64 was -1) and signed-to-unsigned reads wrapped
+    // (-1 read as u64 was 4294967295).
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("int_narrowing.h5");
+    let script = format!(
+        r#"{PRELUDE}
+data = {{
+    "i8": np.array([2**40 + 5, -(2**35), 7, -1, 2**63 - 1, -(2**63)], "<i8"),
+    "i8be": np.array([2**40 + 5, -(2**35), 7, -1], ">i8"),
+    "u8": np.array([2**64 - 1, 2**63, 5, 0], "<u8"),
+    "u4": np.array([2**32 - 1, 2**31, 2**31 - 1, 3], "<u4"),
+    "i4": np.array([-1, 5, -(2**31), 2**31 - 1], "<i4"),
+    "i2be": np.array([-300, 300, -1], ">i2"),
+    "u1": np.array([255, 0, 128], "u1"),
+}}
+with h5py.File("{path}", "w") as f:
+    for name, arr in data.items():
+        f.create_dataset(name, data=arr)
+with h5py.File("{path}", "r") as f:
+    for name in data:
+        d = f[name]
+        emit(name + ":i32", d.astype("<i4")[()])
+        emit(name + ":i64", d.astype("<i8")[()])
+        emit(name + ":u64", d.astype("<u8")[()])
+"#,
+        path = path.display()
+    );
+    let expected = run_python(&script);
+
+    let file = File::open(&path).unwrap();
+    for name in ["i8", "i8be", "u8", "u4", "i4", "i2be", "u1"] {
+        let ds = file.dataset(name).unwrap();
+        assert_eq!(
+            ds.read_i32().unwrap(),
+            parse::<i32>(&expected[&format!("{name}:i32")]),
+            "{name} as i32"
+        );
+        assert_eq!(
+            ds.read_i64().unwrap(),
+            parse::<i64>(&expected[&format!("{name}:i64")]),
+            "{name} as i64"
+        );
+        // libhdf5 wraps a negative big-endian i64 read as little-endian u64
+        // (it only byte-swaps when the sizes match and the order differs);
+        // every other signed-to-unsigned read saturates at 0, so do that.
+        if name != "i8be" {
+            assert_eq!(
+                ds.read_u64().unwrap(),
+                parse::<u64>(&expected[&format!("{name}:u64")]),
+                "{name} as u64"
+            );
+        }
+    }
+    let i8be = file.dataset("i8be").unwrap();
+    assert_eq!(i8be.read_u64().unwrap(), vec![(1 << 40) + 5, 0, 7, 0]);
+}
+
+#[test]
+fn floats_decode_by_their_datatype_fields() {
+    // Every 2-byte float used to decode as IEEE half, so bfloat16 1.5 read as
+    // 1.9375 and +inf as NaN; 1-byte FP8 floats were refused.
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("float_layouts.h5");
+    let script = format!(
+        r#"{PRELUDE}
+def custom(base, fields, bias, size):
+    # fields: (sign pos, exponent pos, exponent size, mantissa pos, mantissa size)
+    t = base.copy()
+    t.set_fields(*fields)
+    t.set_ebias(bias)
+    t.set_precision(size * 8)
+    t.set_size(size)
+    return t
+
+def write(f, name, ftype, raw):
+    raw = np.ascontiguousarray(raw)
+    space = h5py.h5s.create_simple(raw.shape)
+    ds = h5py.h5d.create(f.id, name.encode(), ftype, space)
+    ds.write(h5py.h5s.ALL, h5py.h5s.ALL, raw, mtype=ftype)
+
+# bfloat16: 1.5, -2.25, +inf, 0, 3.140625, 1, -0, smallest subnormal,
+# largest finite, NaN
+bf16 = np.array([0x3FC0, 0xC010, 0x7F80, 0x0000, 0x4049, 0x3F80, 0x8000, 0x0001,
+                 0x7F7F, 0x7FC1], "<u2")
+# 8-bit patterns, all 256 of them
+fp8 = np.arange(256, dtype="u1")
+types = {{
+    "bf16_le": (custom(h5py.h5t.IEEE_F32LE, (15, 7, 8, 0, 7), 127, 2), bf16),
+    "bf16_be": (custom(h5py.h5t.IEEE_F32BE, (15, 7, 8, 0, 7), 127, 2), bf16.byteswap()),
+    "e4m3": (custom(h5py.h5t.IEEE_F32LE, (7, 3, 4, 0, 3), 7, 1), fp8),
+    "e5m2": (custom(h5py.h5t.IEEE_F32LE, (7, 2, 5, 0, 2), 15, 1), fp8),
+}}
+with h5py.File("{path}", "w") as f:
+    for name, (t, raw) in types.items():
+        write(f, name, t, raw)
+    vals = np.array([1.5, -2.25, np.inf, -np.inf, 0.0, -0.0, 6e-8, 65504, 1e-40, 1e300])
+    f.create_dataset("f2_be", data=vals.astype(">f2"))
+    f.create_dataset("f4_be", data=vals.astype(">f4"))
+    f.create_dataset("f8_le", data=vals.astype("<f8"))
+with h5py.File("{path}", "r") as f:
+    for name in list(types) + ["f2_be", "f4_be", "f8_le"]:
+        v = f[name].astype("<f8")[()]
+        # NaN payloads differ between converters; compare NaN as one value.
+        v[np.isnan(v)] = np.nan
+        emit(name, v.view("<u8"))
+"#,
+        path = path.display()
+    );
+    let expected = run_python(&script);
+
+    let canonical = |v: f64| if v.is_nan() { f64::NAN } else { v };
+    let file = File::open(&path).unwrap();
+    for name in [
+        "bf16_le", "bf16_be", "e4m3", "e5m2", "f2_be", "f4_be", "f8_le",
+    ] {
+        let ds = file.dataset(name).unwrap();
+        let want: Vec<u64> = parse(&expected[name]);
+        let got: Vec<u64> = ds
+            .read_f64()
+            .unwrap()
+            .into_iter()
+            .map(|v| canonical(v).to_bits())
+            .collect();
+        assert_eq!(got, want, "{name} as f64");
+        // f32 reads agree too (every value here is exact in f32 except the
+        // f64 dataset, which rounds like `as f32`).
+        let got32: Vec<u32> = ds
+            .read_f32()
+            .unwrap()
+            .into_iter()
+            .map(|v| if v.is_nan() { f32::NAN } else { v }.to_bits())
+            .collect();
+        let want32: Vec<u32> = want
+            .iter()
+            .map(|&b| canonical(f64::from_bits(b)) as f32)
+            .map(|v| if v.is_nan() { f32::NAN } else { v }.to_bits())
+            .collect();
+        assert_eq!(got32, want32, "{name} as f32");
+    }
+}
+
+#[test]
+fn enum_and_bool_datasets_read_as_their_integer_values() {
+    // Enumerations (h5py stores bool as an enum of int8) were refused by the
+    // numeric readers with a type mismatch.
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("enums.h5");
+    let script = format!(
+        r#"{PRELUDE}
+with h5py.File("{path}", "w") as f:
+    f.create_dataset("bool", data=np.array([True, False, True]))
+    e = h5py.enum_dtype({{"RED": 0, "GREEN": 7, "BLUE": -3}}, basetype=">i2")
+    f.create_dataset("enum_i2be", data=np.array([0, 7, -3, 7], ">i2"), dtype=e)
+    e = h5py.enum_dtype({{"LOW": 0, "HIGH": 200}}, basetype="u1")
+    f.create_dataset("enum_u1", data=np.array([200, 0, 200], "u1"), dtype=e)
+    e = h5py.enum_dtype({{"A": -(2**40), "B": 2**40}}, basetype="<i8")
+    f.create_dataset("enum_i8", data=np.array([2**40, -(2**40)], "<i8"), dtype=e)
+with h5py.File("{path}", "r") as f:
+    for name in ("bool", "enum_i2be", "enum_u1", "enum_i8"):
+        emit(name, np.asarray(f[name][()]).astype(np.int64))
+"#,
+        path = path.display()
+    );
+    let expected = run_python(&script);
+
+    let file = File::open(&path).unwrap();
+    for name in ["bool", "enum_i2be", "enum_u1", "enum_i8"] {
+        let ds = file.dataset(name).unwrap();
+        let want: Vec<i64> = parse(&expected[name]);
+        assert_eq!(ds.read_i64().unwrap(), want, "{name} as i64");
+        let want_f64: Vec<f64> = want.iter().map(|&v| v as f64).collect();
+        assert_eq!(ds.read_f64().unwrap(), want_f64, "{name} as f64");
+    }
+    let bools = file.dataset("bool").unwrap();
+    assert_eq!(bools.read_u64().unwrap(), vec![1, 0, 1]);
+    assert_eq!(bools.read_i32().unwrap(), vec![1, 0, 1]);
+}
+
+#[test]
+fn vl_sequences_of_wide_base_types_read_whole() {
+    // read_vl_bytes took the sequence's element count as its byte length, so
+    // [1, 2, 3] as VL int32 came back as 3 bytes instead of 12.
+    use clawhdf5::Selection;
+    use clawhdf5_format::datatype::{Datatype, DatatypeByteOrder};
+    use clawhdf5_format::vl_data::read_vl_bytes;
+
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("vlen.h5");
+    let script = format!(
+        r#"{PRELUDE}
+data = {{
+    "i4": (h5py.vlen_dtype("<i4"), [[1, 2, 3], [], [-5], list(range(40))]),
+    "f8": (h5py.vlen_dtype("<f8"), [[1.5, -2.0], [0.25]]),
+    "u1": (h5py.vlen_dtype("u1"), [[1, 2, 255], []]),
+}}
+with h5py.File("{path}", "w") as f:
+    for name, (dt, seqs) in data.items():
+        ds = f.create_dataset(name, (len(seqs),), dtype=dt)
+        for i, s in enumerate(seqs):
+            ds[i] = s
+with h5py.File("{path}", "r") as f:
+    for name in data:
+        for i, s in enumerate(f[name][()]):
+            emit(f"{{name}}:{{i}}", np.frombuffer(np.asarray(s).tobytes(), "u1"))
+"#,
+        path = path.display()
+    );
+    let expected = run_python(&script);
+
+    let file = File::open(&path).unwrap();
+    let sb = file.superblock();
+    for (name, count) in [("i4", 4), ("f8", 2), ("u1", 2)] {
+        let raw = file
+            .dataset(name)
+            .unwrap()
+            .read_selection(&Selection::All)
+            .unwrap();
+        let got =
+            read_vl_bytes(file.as_bytes(), &raw, count, sb.offset_size, sb.length_size).unwrap();
+        let want: Vec<Vec<u8>> = (0..count)
+            .map(|i| parse(expected.get(&format!("{name}:{i}")).unwrap()))
+            .collect();
+        assert_eq!(got, want, "{name}");
+        if name == "i4" {
+            let i32_le = Datatype::FixedPoint {
+                size: 4,
+                byte_order: DatatypeByteOrder::LittleEndian,
+                signed: true,
+                bit_offset: 0,
+                bit_precision: 32,
+            };
+            let values = clawhdf5_format::data_read::read_as_i64(&got[0], &i32_le).unwrap();
+            assert_eq!(values, vec![1, 2, 3]);
+            assert_eq!(got[3].len(), 40 * 4);
+        }
+    }
+}