diff --git a/BENCHMARKS.md b/BENCHMARKS.md index 969d1ed..a672bbe 100644 --- a/BENCHMARKS.md +++ b/BENCHMARKS.md @@ -149,6 +149,40 @@ vectors, and recall is measured against brute-force ground truth rather than against the f32 index, whose own approximation errors a re-scored search is entitled to get right. +### float16 embedding storage (`MemoryConfig::float16`) + +Measured 2026-09-23 on tank (AMD Ryzen 7 7800X3D). The same clustered +384-dim data in an `f32` store and a `float16` store, both with the default +int8 index and Hebbian boosting off (so every query sees the same store). +Recall is vector-only `hybrid_search` against an exact scan of the original +`f32` vectors, 200 queries. Six runs, three with each store going first; +medians. Nothing depended on the order. + +```bash +cargo run --release -p clawhdf5-bench --bin search_harness -- --float16-study --full +cargo run --release -p clawhdf5-bench --bin search_harness -- --float16-study --full --f16-first +``` + +| N | embeddings | file MiB | checkpoint ms | open ms | recall@10 | top-10 overlap | hybrid p50 ms | +|---:|---|---:|---:|---:|---:|---:|---:| +| 1 000 | f32 | 1.6 | 8 | 1.6 | 1.0000 | | 0.072 | +| 1 000 | float16 | 0.8 | 6 | 1.9 | 0.9980 | 0.9980 | 0.072 | +| 10 000 | f32 | 15.4 | 66 | 15.2 | 1.0000 | | 0.495 | +| 10 000 | float16 | 8.2 | 45 | 18.1 | 1.0000 | 1.0000 | 0.494 | +| 100 000 | f32 | 154.0 | 752 | 299.8 | 0.9940 | | 4.676 | +| 100 000 | float16 | **80.8** | **512** | **252.1** | 0.9990 | 0.9940 | 4.654 | + +The file is 48% smaller, checkpoints write less and open reads less. At +small N opening is slightly slower (widening halves costs more than the I/O it +saves: +3 ms at 10K). Recall does not move: half precision keeps about three +significant digits, far finer than the gaps between neighbours on unit-length +embeddings. Recall was identical in every run; the 0.999 against 0.994 at +100K is two slightly different HNSW graphs, not an improvement to claim. + +The in-memory cache holds the half-rounded values, so the store searches the +same before and after a reopen; RAM use is unchanged (the cache is still +`f32`). What `float16` saves is disk, and the I/O that goes with it. + ### Opening a store (`read_from_disk`) `HDF5Memory::open` memory-mapped the file, copied the whole mapping into a diff --git a/CHANGELOG.md b/CHANGELOG.md index 5f8dc5b..7d2a4f1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,23 @@ ## Unreleased ### Upgrade Notes +- **Files written by clawhdf5 now open in h5py and libhdf5.** Every `f32` + dataset we wrote — including every agent store's embeddings — was refused + with "sign bit position out of bounds", and every empty dataset with + "invalid dataset size". Both were write-side bugs present in every release; + clawhdf5's own reader was unaffected. An agent store is rewritten in full at + each checkpoint, so it becomes readable at its next checkpoint on this + version; other files with `f32` or empty datasets need rewriting. Details in + `docs/known-issues.md`. +- **`MemoryConfig::float16` now does what it says.** It was persisted and + otherwise ignored; embeddings were always stored as `f32`. A store created + with it on now writes half-precision embeddings (48% smaller files) and + rounds embeddings to half precision as they are saved. A store that already + had `float16 = true` rounds its embeddings when next opened and writes them + as `float16` at its next checkpoint. Off by default. +- **Breaking:** `MemoryError` gained `InvalidEntry`, returned when a + `float16` store is given an embedding value beyond ±65504. Exhaustive + matches need the new arm. - **The default build no longer compiles any C.** Deflate now defaults to the pure-Rust zlib-rs instead of zlib-ng, so building the core crates needs neither cmake nor a C compiler. Speed on HDF5 reads and writes is within 6% @@ -23,6 +40,40 @@ `quantized_index = false`, or pass `create --f32-index` to the CLI, to opt out. The CLI's `--quantized-index` is still accepted but is now a no-op. +### Interop +- `clawhdf5-format`: **every `f32` dataset was unreadable by h5py and + libhdf5.** The float datatype encoder hard-coded the sign bit's position to + 63, correct only for `f64`; libhdf5 validates it and refused the dataset. It + is now derived from the type (15 / 31 / 63). Our reader ignores the field, + and the interop suites only wrote `f64`, which is how it went unnoticed. +- `clawhdf5-format`: **every empty dataset was unreadable by h5py and + libhdf5.** It was written with a real address and zero bytes, which trips + libhdf5's `addr + size <= addr` overflow check. An empty contiguous dataset + now gets the undefined address, as libhdf5 writes it. This affected every + agent store without sessions or a knowledge graph. +- New interop tests: `f32` and `float16` datasets in both directions (our + `float16` rounding matches numpy's bit for bit on 4 020 probe values, + including ties, subnormals and the overflow boundary), and an agent store — + `f32` and `float16` — opened by h5py with every dataset decoded. + +### Storage +- `clawhdf5-format`: **half-precision datasets.** + `DatasetBuilder::with_f16_data` writes IEEE binary16 (numpy `float16`), + rounding to nearest-even; `make_f16_type`, and `clawhdf5_format::float16` + with the conversions, which are checked against the `half` crate on 16.7M + values and round-trip all 65 536 half values. Reading `float16` as `f32` + gained a little-endian fast path. +- `clawhdf5-agent`: **`MemoryConfig::float16` stores embeddings as half + precision.** At 100K x 384 the file goes from 154.0 to 80.8 MiB (−48%), a + checkpoint from 752 to 512 ms and open from 300 to 252 ms, with the same + vector recall@10 against an exact scan (0.999 vs 0.994) and the same + `hybrid_search` latency; at 10K open is 3 ms slower. The cache rounds each + embedding as it is saved, so memory and file agree bit for bit and a store + returns the same results before and after a reopen (tested). Out-of-range + values are refused with `MemoryError::InvalidEntry` rather than stored as + infinity; batches are all or nothing. CLI: `create --float16`. See + `BENCHMARKS.md`, "float16 embedding storage". + ### Build - **Pure-Rust default.** `clawhdf5-format`, `clawhdf5-filters` and the `clawhdf5` facade default to the `zlib-rs` deflate backend; `fast-deflate` diff --git a/CLAUDE.md b/CLAUDE.md index cb10c9e..5423c38 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -87,6 +87,16 @@ Cargo workspace with 16 crates under `crates/` (plus `libaec-sys`, an internal F `export` do). An unreadable WAL (torn header, bad magic) is quarantined to `.h5.wal.corrupt-` rather than blocking `open()`; a WAL with an unknown *newer* version still fails and is left untouched. +- `MemoryConfig::float16` (off by default, persisted; CLI `create --float16`) + writes `/memory/embeddings` as IEEE half precision (48% smaller file at + 100K, same recall). `MemoryCache::half_precision` rounds each embedding as + it enters the cache (push, update, WAL replay, and on load of a store still + `f32` on disk), so memory and file agree bit for bit; the conversions live + in `clawhdf5_format::float16` and must stay the single implementation. + Values beyond ±65504 are `MemoryError::InvalidEntry`. Interop: every file + must open in h5py — `f32` datasets and empty datasets did not until + 2026-09-23 (see `docs/known-issues.md`); the agent's `h5py_interop` test + guards a whole store. - `MemoryConfig::compression` is off by default; when on, embeddings are deflate-compressed, or Zstd with the agent's `zstd` feature (links libzstd). - `Dataset::verify_provenance()` (clawhdf5 facade, `provenance` feature, on by diff --git a/README.md b/README.md index 58799af..f661c4b 100644 --- a/README.md +++ b/README.md @@ -82,6 +82,15 @@ breaking change, are in [CHANGELOG.md](CHANGELOG.md). 100K × 384 store to 1.74× the raw vectors. At equal recall it is also faster than `f32`: 1.63× QPS on AVX2, 1.18× on a Raspberry Pi 5 (NEON `SDOT`). +**Interop (unreleased)** +- **Files we write now open in h5py and libhdf5.** Every `f32` dataset — + including every agent store's embeddings — and every empty dataset was + refused by libhdf5. Both were write-side bugs in every release; agent stores + fix themselves at their next checkpoint. See + [docs/known-issues.md](docs/known-issues.md). +- `MemoryConfig::float16` now stores half-precision embeddings (it was + ignored): 48% smaller files at the same recall. + **Tooling** - CI now runs the h5py/netCDF4 interop suites for real (they had been skipping silently) and runs an aarch64 job for the NEON kernels. @@ -244,6 +253,9 @@ retrieval recall reported as QA accuracy typically overstates by 20–30 points. | 10K | 17.0 MB | 1.7 KB | 2.7 MB (6.2x) | | 100K | 169.8 MB | 1.7 KB | 26.9 MB (6.2x) | +With `MemoryConfig::float16` the embeddings take half the space: an agent +store of 100K × 384 records is 80.8 MiB instead of 154.0. + **In memory** — a store reopened from disk, 384-dim `f32`, measured with a counting allocator ([BENCHMARKS.md § Memory footprint](BENCHMARKS.md#memory-footprint)): @@ -538,7 +550,7 @@ ClawhDF5's agent memory design draws from 15+ recent papers: | Flag | Default | Description | |------|---------|-------------| -| `float16` | **yes** | Half-precision cosine kernel (`cosine_similarity_f16`). The store itself always writes `f32` embeddings; `MemoryConfig::float16` is recorded in `/meta` but not yet applied | +| `float16` | **yes** | Half-precision cosine kernel (`cosine_similarity_f16`). Half-precision *storage* is the `MemoryConfig::float16` setting below, and needs no feature | | `hnsw` | **yes** | HNSW approximate vector index for `hybrid_search` (via `clawhdf5-ann`); disable for an exact linear scan | | `parallel` | **yes** | Parallel HNSW bulk build (same graph, ~3× faster on 16 cores) and Rayon brute-force search strategies | | `zstd` | no | Compress embeddings with Zstd instead of deflate when `MemoryConfig::compression` is on (links libzstd) | @@ -568,6 +580,14 @@ setting existed keep their `f32` index; opt out for new stores with `quantized_index = false` or `clawhdf5-cli create --f32-index`. See [BENCHMARKS.md § Quantising the index copy](BENCHMARKS.md#quantising-the-index-copy-quantized_index). +`MemoryConfig::float16` (off by default; CLI `create --float16`) stores the +embeddings on disk as IEEE half precision (numpy `float16`): at 100K × 384 the +file drops from 154 to 81 MiB, checkpoints and opens get faster, and vector +recall and search latency do not change. Embeddings are rounded as they are +saved, so the store searches the same before and after a reopen; values must +lie within ±65504. See +[BENCHMARKS.md § float16 embedding storage](BENCHMARKS.md#float16-embedding-storage-memoryconfigfloat16). + ### `clawhdf5-format` | Flag | Default | Description | @@ -656,8 +676,9 @@ agent_memory.h5 │ └── ann_generation (ties the .ann sidecar to this checkpoint) ├── /memory │ ├── chunks: string[N] -│ ├── embeddings: f32[N × D] (chunked; deflate, or Zstd with the -│ │ `zstd` feature, when compression is on) +│ ├── embeddings: f32[N × D], or f16 for a `float16` store +│ │ (chunked; deflate, or Zstd with the `zstd` +│ │ feature, when compression is on) │ ├── source_channel: string[N] │ ├── timestamps: f64[N] │ ├── session_ids: string[N] diff --git a/crates/clawhdf5-agent/src/cache.rs b/crates/clawhdf5-agent/src/cache.rs index e69da1b..37cd001 100644 --- a/crates/clawhdf5-agent/src/cache.rs +++ b/crates/clawhdf5-agent/src/cache.rs @@ -1,6 +1,7 @@ //! In-memory cache for memory entries, sessions, and knowledge graph. use crate::vector_search; +use clawhdf5_format::float16::round_to_f16; /// Every entry's embedding, in one contiguous `[N x dim]` buffer. /// @@ -149,6 +150,11 @@ pub struct MemoryCache { pub norms: Vec, /// Hebbian activation weights (default 1.0 per entry). pub activation_weights: Vec, + /// Round every embedding to IEEE half precision as it enters the cache, + /// so the cache holds exactly what a `float16` store writes to disk. Set + /// it with [`MemoryCache::set_half_precision`], which also rounds the + /// rows already held. + pub half_precision: bool, } impl MemoryCache { @@ -164,9 +170,44 @@ impl MemoryCache { embedding_dim, norms: Vec::new(), activation_weights: Vec::new(), + half_precision: false, } } + /// Switch half-precision rounding on or off. Turning it on rounds every + /// embedding already held (and recomputes norms where one changed) — + /// e.g. a `float16` store whose last checkpoint predates half-precision + /// storage and so is still `f32` on disk. + pub fn set_half_precision(&mut self, on: bool) { + self.half_precision = on; + if !on { + return; + } + for i in 0..self.embeddings.len() { + let row = &self.embeddings[i]; + if row + .iter() + .all(|&v| round_to_f16(v).to_bits() == v.to_bits()) + { + continue; + } + let rounded: Vec = row.iter().map(|&v| round_to_f16(v)).collect(); + self.norms[i] = vector_search::compute_norm(&rounded); + self.embeddings.set(i, &rounded); + } + } + + /// The embedding as the cache will hold it: rounded to half precision + /// when [`Self::half_precision`] is on, otherwise unchanged. + fn stored_form(&self, mut embedding: Vec) -> Vec { + if self.half_precision { + for v in &mut embedding { + *v = round_to_f16(*v); + } + } + embedding + } + /// Kept for callers that used to have to re-flatten after a bulk load. /// The buffer is always flat now, so there is nothing to rebuild. #[deprecated(note = "embeddings are stored flat; this is a no-op")] @@ -202,6 +243,7 @@ impl MemoryCache { tags: String, ) -> usize { let idx = self.chunks.len(); + let embedding = self.stored_form(embedding); let norm = vector_search::compute_norm(&embedding); self.chunks.push(chunk); self.embeddings.push(&embedding); @@ -240,6 +282,7 @@ impl MemoryCache { session_id: String, ) { if idx < self.chunks.len() { + let embedding = self.stored_form(embedding); let norm = vector_search::compute_norm(&embedding); self.chunks[idx] = chunk; self.embeddings.set(idx, &embedding); @@ -439,4 +482,60 @@ mod tests { .reset_from(2, vec![vec![1.0, 2.0], vec![3.0, 4.0]]); assert_eq!(cache.embeddings.as_flat(), vec![1.0, 2.0, 3.0, 4.0]); } + + #[test] + fn set_half_precision_rounds_existing_rows_and_their_norms() { + // A store with float16 set whose checkpoint is still f32 on disk + // loads full-precision rows; switching rounding on must bring them to + // exactly what the next checkpoint will write. + let mut cache = MemoryCache::new(3); + cache.push( + "a".into(), + vec![0.1, 0.2, 0.3], + "c".into(), + 0.0, + "s".into(), + "".into(), + ); + cache.push( + "b".into(), + vec![0.5, 0.25, 1.0], + "c".into(), + 0.0, + "s".into(), + "".into(), + ); + let exact_norm = cache.norms[0]; + + cache.set_half_precision(true); + let row0: Vec = [0.1f32, 0.2, 0.3] + .iter() + .map(|&v| round_to_f16(v)) + .collect(); + assert_eq!(&cache.embeddings[0], row0.as_slice()); + assert_eq!(cache.norms[0], vector_search::compute_norm(&row0)); + assert_ne!(cache.norms[0], exact_norm); + // Already representable: untouched. + assert_eq!(&cache.embeddings[1], &[0.5, 0.25, 1.0]); + + // New rows are rounded as they arrive, and updates too. + cache.push( + "c".into(), + vec![0.1, 0.0, 0.0], + "c".into(), + 0.0, + "s".into(), + "".into(), + ); + assert_eq!(cache.embeddings[2][0], round_to_f16(0.1)); + cache.update( + 2, + "c".into(), + vec![0.3, 0.0, 0.0], + "c".into(), + 0.0, + "s".into(), + ); + assert_eq!(cache.embeddings[2][0], round_to_f16(0.3)); + } } diff --git a/crates/clawhdf5-agent/src/lib.rs b/crates/clawhdf5-agent/src/lib.rs index 5d660ea..4607449 100644 --- a/crates/clawhdf5-agent/src/lib.rs +++ b/crates/clawhdf5-agent/src/lib.rs @@ -63,6 +63,7 @@ use std::path::{Path, PathBuf}; use cache::MemoryCache; #[cfg(feature = "hnsw")] use clawhdf5_ann::{DistanceMetric, HnswIndex, Storage}; +use clawhdf5_format::float16::round_to_f16; use ephemeral::{EphemeralConfig, EphemeralStore}; // EphemeralEntry and EphemeralStats are part of the crate public API via @@ -83,6 +84,9 @@ pub enum MemoryError { NotFound(String), /// Another `HDF5Memory` (in this or another process) has the store open. Locked(String), + /// A record the store cannot hold as given, e.g. an embedding value + /// outside the half-precision range of a `float16` store. + InvalidEntry(String), } impl std::fmt::Display for MemoryError { @@ -93,6 +97,7 @@ impl std::fmt::Display for MemoryError { MemoryError::Schema(e) => write!(f, "schema error: {e}"), MemoryError::NotFound(e) => write!(f, "not found: {e}"), MemoryError::Locked(e) => write!(f, "store is locked: {e}"), + MemoryError::InvalidEntry(e) => write!(f, "invalid entry: {e}"), } } } @@ -124,6 +129,12 @@ pub struct MemoryConfig { pub embedding_dim: usize, pub chunk_size: usize, pub overlap: usize, + /// Store embeddings as IEEE half precision (numpy `float16`): half the + /// bytes of the embeddings dataset on disk. Every embedding is rounded to + /// the nearest half as it enters the store, in memory as well as on disk, + /// so search results are the same before and after a reopen. Values must + /// lie within ±65504; a save outside that is `MemoryError::InvalidEntry`. + /// Fixed when the store is created (persisted in `/meta`). pub float16: bool, pub compression: bool, pub compression_level: u32, @@ -317,7 +328,8 @@ impl HDF5Memory { /// Create a new HDF5 memory file with the given configuration. pub fn create(config: MemoryConfig) -> Result { let lock = store_lock::StoreLock::acquire(&config.path)?; - let cache = MemoryCache::new(config.embedding_dim); + let mut cache = MemoryCache::new(config.embedding_dim); + cache.set_half_precision(config.float16); let sessions = SessionCache::new(); let knowledge = KnowledgeCache::new(); @@ -1070,7 +1082,29 @@ impl HDF5Memory { /// Upsert: if an active entry with the same tags (key) exists, update it in-place. /// Otherwise append a new entry. Use this for key-based memory stores where /// the same key should not create duplicates. + /// A `float16` store holds embeddings as IEEE half precision, which has no + /// finite value beyond ±65504. Refuse such an embedding rather than + /// silently store infinity. (Values that are already infinite or NaN are + /// stored as they are, as in an `f32` store.) + fn check_embedding(&self, embedding: &[f32]) -> Result<()> { + if !self.config.float16 { + return Ok(()); + } + let overflow = embedding + .iter() + .enumerate() + .find(|&(_, &v)| v.is_finite() && round_to_f16(v).is_infinite()); + match overflow { + None => Ok(()), + Some((i, v)) => Err(MemoryError::InvalidEntry(format!( + "embedding[{i}] = {v} is outside the half-precision range (±65504) \ + of this float16 store" + ))), + } + } + pub fn save_or_update(&mut self, entry: MemoryEntry) -> Result { + self.check_embedding(&entry.embedding)?; if let Some(existing_idx) = self.cache.find_by_tags(&entry.tags) { if let Some(ref mut w) = self.wal { let wal_entry = wal::WalEntry { @@ -1126,6 +1160,7 @@ impl HDF5Memory { impl AgentMemory for HDF5Memory { fn save(&mut self, entry: MemoryEntry) -> Result { + self.check_embedding(&entry.embedding)?; if let Some(ref mut w) = self.wal { let wal_entry = wal::WalEntry { entry_type: wal::WalEntryType::Save, @@ -1167,6 +1202,10 @@ impl AgentMemory for HDF5Memory { } fn save_batch(&mut self, entries: Vec) -> Result> { + // All or nothing: check every entry before storing any. + for entry in &entries { + self.check_embedding(&entry.embedding)?; + } let mut indices = Vec::with_capacity(entries.len()); for entry in entries { let idx = self.cache.push( @@ -1327,6 +1366,9 @@ impl HDF5Memory { })?; let view = memory_strategy::CacheStoreView::new(&self.cache, &self.knowledge); let output = strat.evaluate(&exchange, &view); + for e in &output.entries { + self.check_embedding(&e.embedding)?; + } for e in &output.entries { self.cache.push( e.chunk.clone(), @@ -1411,6 +1453,16 @@ impl HDF5Memory { let mut promoted = 0; for key in candidates { + // Check before taking, so a rejected entry stays in the ephemeral + // tier rather than being lost. + if let Some(emb) = self + .ephemeral + .as_ref() + .and_then(|s| s.get_entry(&key)) + .and_then(|e| e.embedding.as_deref()) + { + self.check_embedding(emb)?; + } let entry = match self .ephemeral .as_mut() diff --git a/crates/clawhdf5-agent/src/schema.rs b/crates/clawhdf5-agent/src/schema.rs index a2fd15a..3555986 100644 --- a/crates/clawhdf5-agent/src/schema.rs +++ b/crates/clawhdf5-agent/src/schema.rs @@ -159,20 +159,27 @@ fn build_memory_group( // chunks: fixed-length string array write_string_dataset(&mut group, "chunks", &cache.chunks); - // embeddings: f32 [N x D] + // embeddings: [N x D], f32 — or IEEE half precision for a `float16` + // store. The cache already holds half-rounded values then, so this + // conversion is exact and a reopened store sees the same numbers. let n = cache.embeddings.len() as u64; let d = cache.embedding_dim as u64; let flat = cache.flat_embeddings(); { - let ds = group - .create_dataset("embeddings") - .with_f32_data(flat) - .with_shape(&[n, d]); + let ds = group.create_dataset("embeddings"); + let elem_bytes: u64 = if config.float16 { + ds.with_f16_data(flat); + 2 + } else { + ds.with_f32_data(flat); + 4 + }; + ds.with_shape(&[n, d]); // Chunk size tuning: target ~256KB per chunk for optimal I/O if n > 0 && d > 0 { let target_chunk_bytes: u64 = 256 * 1024; - let rows_per_chunk = (target_chunk_bytes / (d * 4)).max(1).min(n); + let rows_per_chunk = (target_chunk_bytes / (d * elem_bytes)).max(1).min(n); ds.with_chunks(&[rows_per_chunk, d]); // Compression. Shuffle is applied automatically (auto-shuffle @@ -513,7 +520,16 @@ pub fn validate_and_load( }; // Load /memory group - let memory_cache = load_memory_group(file, embedding_dim)?; + let mut memory_cache = load_memory_group(file, embedding_dim)?; + // A float16 store's cache holds half-rounded embeddings. Embeddings read + // from an f16 dataset already are; a float16 store whose last checkpoint + // predates half-precision storage is still f32 on disk and is rounded + // here. + if config.float16 && embeddings_are_f16(file) { + memory_cache.half_precision = true; + } else { + memory_cache.set_half_precision(config.float16); + } // Load /sessions group let session_cache = load_sessions_group(file)?; @@ -756,6 +772,13 @@ fn read_string_dataset_from_group( .map_err(|e| MemoryError::Hdf5(format!("cannot read strings from {name}: {e}"))) } +/// Whether `/memory/embeddings` is stored as IEEE half precision. +fn embeddings_are_f16(file: &clawhdf5::File) -> bool { + file.dataset("memory/embeddings") + .and_then(|ds| ds.dtype()) + .is_ok_and(|dt| matches!(dt, clawhdf5::DType::Other(ref s) if s == "float16")) +} + fn read_f32_dataset(group: &clawhdf5::Group<'_>, name: &str) -> Result, MemoryError> { let ds = group .dataset(name) diff --git a/crates/clawhdf5-agent/tests/float16_store.rs b/crates/clawhdf5-agent/tests/float16_store.rs new file mode 100644 index 0000000..5c83107 --- /dev/null +++ b/crates/clawhdf5-agent/tests/float16_store.rs @@ -0,0 +1,208 @@ +//! `MemoryConfig::float16`: embeddings stored as IEEE half precision. +//! +//! The setting used to be recorded in `/meta` and otherwise ignored — the +//! embeddings dataset was always `f32`. These tests pin what it now does: the +//! dataset is `float16`, the in-memory cache holds exactly the values the file +//! holds (so search results survive a reopen bit for bit), and a value half +//! precision cannot represent is refused rather than stored as infinity. + +use std::path::{Path, PathBuf}; + +use clawhdf5_agent::{AgentMemory, HDF5Memory, MemoryConfig, MemoryEntry, MemoryError}; +use clawhdf5_format::float16::round_to_f16; +use tempfile::TempDir; + +const DIM: usize = 64; + +/// Deterministic, embedding-like unit vectors. +fn embedding(seed: u64) -> Vec { + let mut x = seed.wrapping_mul(0x9E37_79B9_7F4A_7C15) | 1; + let v: Vec = (0..DIM) + .map(|_| { + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + (x >> 40) as f32 / (1u64 << 24) as f32 - 0.5 + }) + .collect(); + let norm = v.iter().map(|a| a * a).sum::().sqrt(); + v.iter().map(|a| a / norm).collect() +} + +fn entry(i: u64) -> MemoryEntry { + MemoryEntry { + chunk: format!("memory number {i} about topic {}", i % 7), + embedding: embedding(i), + source_channel: "test".into(), + timestamp: i as f64, + session_id: "s".into(), + tags: format!("t{i}"), + } +} + +fn config(dir: &TempDir, name: &str, float16: bool) -> MemoryConfig { + let mut c = MemoryConfig::new(dir.path().join(name), "agent", DIM); + c.float16 = float16; + c +} + +fn embeddings_dtype_and_values(path: &Path) -> (String, Vec) { + let file = clawhdf5::File::open(path).unwrap(); + let ds = file.dataset("memory/embeddings").unwrap(); + (format!("{:?}", ds.dtype().unwrap()), ds.read_f32().unwrap()) +} + +fn search_bits(m: &mut HDF5Memory, q: u64) -> Vec<(usize, u32)> { + m.hybrid_search(&embedding(q), "memory topic 3", 0.4, 0.6, 10) + .iter() + .map(|r| (r.index, r.score.to_bits())) + .collect() +} + +#[test] +fn float16_store_writes_half_precision_and_reopens_identically() { + let dir = TempDir::new().unwrap(); + // Two identical stores. Search is not read-only (it boosts the Hebbian + // activation of what it returns, and checkpoints persist that), so each + // is queried exactly once: one live, one after a checkpoint and reopen. + let live_cfg = config(&dir, "live.h5", true); + let cfg = config(&dir, "f16.h5", true); + let path: PathBuf = cfg.path.clone(); + + let mut live = HDF5Memory::create(live_cfg).unwrap(); + live.save_batch((0..200).map(entry).collect()).unwrap(); + let mut m = HDF5Memory::create(cfg).unwrap(); + m.save_batch((0..200).map(entry).collect()).unwrap(); + drop(m); + + // On disk: a genuine float16 dataset holding the rounded inputs. + let (dtype, values) = embeddings_dtype_and_values(&path); + assert_eq!(dtype, "Other(\"float16\")"); + let expected: Vec = (0..200) + .flat_map(|i| embedding(i).into_iter().map(|v| round_to_f16(v).to_bits())) + .collect(); + let got: Vec = values.iter().map(|v| v.to_bits()).collect(); + assert_eq!(got, expected); + + // Reopened, the store answers exactly as the live one does: the cache + // held the half-rounded values before the checkpoint. + let mut reopened = HDF5Memory::open(&path).unwrap(); + for q in 0..5 { + assert_eq!( + search_bits(&mut live, 1000 + q), + search_bits(&mut reopened, 1000 + q), + "query {q}" + ); + } +} + +#[test] +fn float16_halves_the_embeddings_on_disk() { + let dir = TempDir::new().unwrap(); + let mut sizes = Vec::new(); + for float16 in [false, true] { + let cfg = config(&dir, &format!("s{float16}.h5"), float16); + let path = cfg.path.clone(); + let mut m = HDF5Memory::create(cfg).unwrap(); + m.save_batch((0..2000).map(entry).collect()).unwrap(); + drop(m); + sizes.push(std::fs::metadata(&path).unwrap().len()); + } + let embedding_bytes_f32 = (2000 * DIM * 4) as u64; + let saved = sizes[0] - sizes[1]; + // Half of the f32 embeddings, give or take metadata and alignment. + assert!( + saved.abs_diff(embedding_bytes_f32 / 2) < 16 * 1024, + "f32 {} B, f16 {} B, saved {saved} B, expected ~{} B", + sizes[0], + sizes[1], + embedding_bytes_f32 / 2 + ); +} + +#[test] +fn f32_store_is_unchanged() { + let dir = TempDir::new().unwrap(); + let cfg = config(&dir, "f32.h5", false); + let path = cfg.path.clone(); + let mut m = HDF5Memory::create(cfg).unwrap(); + m.save_batch((0..50).map(entry).collect()).unwrap(); + drop(m); + let (dtype, values) = embeddings_dtype_and_values(&path); + assert_eq!(dtype, "F32"); + let expected: Vec = (0..50).flat_map(embedding).collect(); + assert_eq!(values, expected); +} + +#[test] +fn out_of_range_values_are_refused_not_stored_as_infinity() { + let dir = TempDir::new().unwrap(); + let mut cfg = config(&dir, "range.h5", true); + cfg.wal_enabled = true; + let path = cfg.path.clone(); + let mut m = HDF5Memory::create(cfg).unwrap(); + m.save(entry(1)).unwrap(); + + let mut bad = entry(2); + bad.embedding[5] = 70_000.0; + match m.save(bad.clone()) { + Err(MemoryError::InvalidEntry(msg)) => assert!(msg.contains("embedding[5]"), "{msg}"), + other => panic!("expected InvalidEntry, got {other:?}"), + } + assert!(matches!( + m.save_or_update(bad.clone()), + Err(MemoryError::InvalidEntry(_)) + )); + // A batch is all or nothing. + assert!(matches!( + m.save_batch(vec![entry(3), bad.clone(), entry(4)]), + Err(MemoryError::InvalidEntry(_)) + )); + assert_eq!(m.count(), 1); + + // The largest finite half, and values that round down to it, are fine. + let mut edge = entry(5); + edge.embedding[0] = 65504.0; + edge.embedding[1] = -65519.0; + m.save(edge).unwrap(); + assert_eq!(m.count(), 2); + drop(m); + + // Nothing rejected reached the WAL or the file. + let m = HDF5Memory::open(&path).unwrap(); + assert_eq!(m.count(), 2); + + // An f32 store takes the same value as it always did. + let mut m32 = HDF5Memory::create(config(&dir, "range32.h5", false)).unwrap(); + m32.save(bad).unwrap(); +} + +#[test] +fn wal_replay_rounds_like_a_live_save() { + let dir = TempDir::new().unwrap(); + let mut cfg = config(&dir, "wal.h5", true); + cfg.wal_enabled = true; + cfg.wal_max_entries = 10_000; // keep everything in the WAL + let path = cfg.path.clone(); + let mut m = HDF5Memory::create(cfg).unwrap(); + for i in 0..30 { + m.save(entry(i)).unwrap(); + } + let live = search_bits(&mut m, 77); + + // Crash image: the .h5 is still the empty checkpoint; everything is in + // the WAL, which holds the caller's f32 values. + let crash = TempDir::new().unwrap(); + let image = crash.path().join("image.h5"); + std::fs::copy(&path, &image).unwrap(); + std::fs::copy( + path.with_extension("h5.wal"), + image.with_extension("h5.wal"), + ) + .unwrap(); + drop(m); + + let mut recovered = HDF5Memory::open(&image).unwrap(); + assert_eq!(recovered.count(), 30); + assert_eq!(search_bits(&mut recovered, 77), live); +} diff --git a/crates/clawhdf5-agent/tests/h5py_interop.rs b/crates/clawhdf5-agent/tests/h5py_interop.rs new file mode 100644 index 0000000..d82ed46 --- /dev/null +++ b/crates/clawhdf5-agent/tests/h5py_interop.rs @@ -0,0 +1,94 @@ +//! An agent store is a standard HDF5 file: h5py can open it and read every +//! dataset. +//! +//! It could not: the float datatype's sign-bit position was hard-coded for +//! f64, so every f32 dataset (embeddings, norms, activation weights) made +//! libhdf5 refuse the file with "sign bit position out of bounds". + +use std::process::Command; + +use clawhdf5_agent::{AgentMemory, HDF5Memory, MemoryConfig, MemoryEntry}; + +fn python() -> String { + std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string()) +} + +fn h5py_available() -> bool { + Command::new(python()) + .args(["-c", "import h5py"]) + .output() + .map(|o| o.status.success()) + .unwrap_or(false) +} + +#[test] +fn h5py_reads_every_dataset_of_an_agent_store() { + if !h5py_available() { + assert!( + std::env::var("CLAWHDF5_REQUIRE_INTEROP").as_deref() != Ok("1"), + "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available" + ); + eprintln!("SKIP: python3 with h5py not available"); + return; + } + let dir = tempfile::tempdir().unwrap(); + for float16 in [false, true] { + let path = dir.path().join(format!("store_{float16}.h5")); + let mut cfg = MemoryConfig::new(path.clone(), "agent", 8); + cfg.float16 = float16; + let mut m = HDF5Memory::create(cfg).unwrap(); + // save_batch checkpoints, so the records are in the .h5, not the WAL. + m.save_batch( + (0..20) + .map(|i| MemoryEntry { + chunk: format!("memory {i}"), + embedding: (0..8).map(|j| ((i * 8 + j) as f32).sin()).collect(), + source_channel: "test".into(), + timestamp: i as f64, + session_id: "s".into(), + tags: String::new(), + }) + .collect(), + ) + .unwrap(); + drop(m); + + // Exact expected values, as bits: numpy's sin need not match Rust's + // to the last place. + let bits = (0..160) + .map(|k| (k as f32).sin().to_bits().to_string()) + .collect::>() + .join(","); + let script = format!( + r#" +import h5py, numpy as np +want = np.float16 if {py_bool} else np.float32 +with h5py.File("{path}", "r") as f: + names = [] + f.visititems(lambda n, o: names.append(n) if isinstance(o, h5py.Dataset) else None) + for n in names: + f[n][()] # every dataset must decode + e = f["memory/embeddings"] + assert e.dtype == want, e.dtype + assert e.shape == (20, 8), e.shape + ref = np.array([{bits}], dtype=np.uint32).view(np.float32).astype(want).reshape(20, 8) + assert (e[()] == ref).all() + assert f["memory/norms"].dtype == np.float32 +print(len(names)) +"#, + py_bool = if float16 { "True" } else { "False" }, + path = path.display() + ); + let out = Command::new(python()) + .args(["-c", &script]) + .output() + .unwrap(); + assert!( + out.status.success(), + "float16={float16}: {}", + String::from_utf8_lossy(&out.stderr) + ); + let n: usize = String::from_utf8_lossy(&out.stdout).trim().parse().unwrap(); + assert!(n >= 10, "only {n} datasets"); + } +} diff --git a/crates/clawhdf5-bench/src/bin/search_harness.rs b/crates/clawhdf5-bench/src/bin/search_harness.rs index 951f0de..1df37a9 100644 --- a/crates/clawhdf5-bench/src/bin/search_harness.rs +++ b/crates/clawhdf5-bench/src/bin/search_harness.rs @@ -19,6 +19,7 @@ //! cargo run --release -p clawhdf5-bench --bin search_harness -- --full # + 100K //! cargo run --release -p clawhdf5-bench --bin search_harness -- --json out.json //! cargo run --release -p clawhdf5-bench --bin search_harness -- --ann-only --uniform +//! cargo run --release -p clawhdf5-bench --bin search_harness -- --float16-study --full //! ``` use std::time::{Duration, Instant}; @@ -88,6 +89,9 @@ static UNIFORM: std::sync::atomic::AtomicBool = std::sync::atomic::AtomicBool::n /// the memory) instead of f32, to price the recall it costs. static INT8: std::sync::atomic::AtomicBool = std::sync::atomic::AtomicBool::new(false); +/// `--f16-first`: in `--float16-study`, run the float16 store first. +static F16_FIRST: std::sync::atomic::AtomicBool = std::sync::atomic::AtomicBool::new(false); + /// `--rerank`: re-score the candidate pool against the exact vectors before /// taking the top K. static RERANK: std::sync::atomic::AtomicBool = std::sync::atomic::AtomicBool::new(false); @@ -483,6 +487,148 @@ fn bench_end_to_end(n: usize, json: &mut Vec) { })); } +// --------------------------------------------------------------------------- +// float16 study: what does half-precision embedding storage cost? +// --------------------------------------------------------------------------- + +/// `--float16-study`: the same data in an `f32` store and a `float16` store. +/// Reports file size, checkpoint and open time, vector-search recall@10 +/// against an exact scan of the *original* f32 vectors, how often the two +/// stores return the same top 10, and `hybrid_search` latency. Hebbian +/// boosting is off, so every query sees the same store. +fn float16_study(n: usize) { + let data = make_dataset(n, 0xF16 ^ n as u64); + let mut rng = Rng(11); + let query_texts: Vec = data + .query_cluster + .iter() + .enumerate() + .map(|(i, c)| text_for(*c, i, &mut rng)) + .collect(); + + // Exact top K by cosine (the vectors are unit length) on the f32 inputs. + let exact: Vec> = data + .queries + .iter() + .map(|q| { + let mut scored: Vec<(usize, f32)> = data + .vectors + .iter() + .enumerate() + .map(|(i, v)| (i, v.iter().zip(q).map(|(a, b)| a * b).sum())) + .collect(); + scored.sort_by(|a, b| b.1.total_cmp(&a.1).then(a.0.cmp(&b.0))); + scored.into_iter().take(K).map(|(i, _)| i).collect() + }) + .collect(); + + let dir = tempfile::tempdir().unwrap(); + let mut per_variant: Vec<(bool, Vec>)> = Vec::new(); + // `--f16-first` swaps the order, to check the numbers do not depend on + // which store runs first (page cache, allocator, CPU frequency). + let order = if F16_FIRST.load(std::sync::atomic::Ordering::Relaxed) { + [true, false] + } else { + [false, true] + }; + for float16 in order { + let path = dir.path().join(format!("f16study_{float16}.h5")); + let mut rng = Rng(3); + let entries: Vec = data + .vectors + .iter() + .enumerate() + .map(|(i, v)| MemoryEntry { + chunk: text_for(data.cluster_of[i], i, &mut rng), + embedding: v.clone(), + source_channel: "bench".into(), + timestamp: i as f64, + session_id: format!("s{}", i % 50), + tags: format!("t{i}"), + }) + .collect(); + let mut config = MemoryConfig::new(path.clone(), "bench", DIM); + config.float16 = float16; + config.hebbian_boost = 0.0; + let mut mem = HDF5Memory::create(config).unwrap(); + mem.save_batch(entries).unwrap(); + // Build the indexes, then time a checkpoint that writes everything. + std::hint::black_box(mem.hybrid_search(&data.queries[0], "", 1.0, 0.0, K)); + let t = Instant::now(); + mem.flush_wal().unwrap(); + let checkpoint = t.elapsed(); + drop(mem); + let file_bytes = std::fs::metadata(&path).unwrap().len(); + + // Median of three opens. + let mut opens: Vec = (0..3) + .map(|_| { + let t = Instant::now(); + let m = HDF5Memory::open(&path).unwrap(); + let d = t.elapsed(); + drop(m); + d + }) + .collect(); + opens.sort(); + let mut mem = HDF5Memory::open(&path).unwrap(); + + // Vector-only search: empty text, all weight on the vector stage. + let results: Vec> = data + .queries + .iter() + .map(|q| { + mem.hybrid_search(q, "", 1.0, 0.0, K) + .iter() + .map(|r| r.index) + .collect() + }) + .collect(); + let hits: usize = results + .iter() + .zip(&exact) + .map(|(got, want)| got.iter().filter(|i| want.contains(i)).count()) + .sum(); + let recall = hits as f64 / (K * data.queries.len()) as f64; + + let latency = summarize( + (0..N_QUERIES) + .map(|i| { + let t = Instant::now(); + std::hint::black_box(mem.hybrid_search( + &data.queries[i], + &query_texts[i], + 0.4, + 0.6, + K, + )); + t.elapsed() + }) + .collect(), + ); + let overlap = match per_variant.first() { + Some((_, other)) => { + let same: usize = results + .iter() + .zip(other) + .map(|(a, b)| a.iter().filter(|i| b.contains(i)).count()) + .sum(); + format!("{:.4}", same as f64 / (K * data.queries.len()) as f64) + } + None => "—".into(), + }; + println!( + "| {n} | {} | {:.1} | {:.0} | {:.1} | {recall:.4} | {overlap} | {:.3} |", + if float16 { "float16" } else { "f32" }, + mib(file_bytes), + millis(checkpoint), + millis(opens[1]), + millis(latency.p50), + ); + per_variant.push((float16, results)); + } +} + // --------------------------------------------------------------------------- // Fusion study: does capping the keyword candidate pool change the ranking? // --------------------------------------------------------------------------- @@ -642,6 +788,24 @@ fn main() { } return; } + if args.iter().any(|a| a == "--f16-first") { + F16_FIRST.store(true, std::sync::atomic::Ordering::Relaxed); + } + if args.iter().any(|a| a == "--float16-study") { + println!("## float16 embedding storage ({DIM}-dim, int8 index, Hebbian boost off)\n"); + println!( + "| N | embeddings | file MiB | checkpoint ms | open ms | recall@10 | top-10 overlap with the other | hybrid p50 ms |" + ); + println!("|---:|---|---:|---:|---:|---:|---:|---:|"); + for &n in if full { + &[1_000, 10_000, 100_000][..] + } else { + &[1_000, 10_000][..] + } { + float16_study(n); + } + return; + } if args.iter().any(|a| a == "--int8") { INT8.store(true, std::sync::atomic::Ordering::Relaxed); println!("(int8-quantised index vectors)"); diff --git a/crates/clawhdf5-cli/src/main.rs b/crates/clawhdf5-cli/src/main.rs index ed02190..c3f15e3 100644 --- a/crates/clawhdf5-cli/src/main.rs +++ b/crates/clawhdf5-cli/src/main.rs @@ -36,6 +36,11 @@ enum Commands { /// Accepted for compatibility; int8 is now the default #[arg(long, hide = true, conflicts_with = "f32_index")] quantized_index: bool, + /// Store embeddings on disk as IEEE half precision (float16): half + /// the bytes, about three significant digits; values must lie within + /// ±65504 + #[arg(long)] + float16: bool, }, /// Save a memory entry (reads JSON from stdin or --json) Save { @@ -102,9 +107,11 @@ fn run(cli: Cli) -> Result<(), Box> { wal, f32_index, quantized_index: _, + float16, } => { let mut config = MemoryConfig::new(cli.path.clone(), &agent_id, dim); config.wal_enabled = wal; + config.float16 = float16; // Only ever switch *off* the library default: assigning the flag // outright would force every CLI-created store back to f32 unless // the caller knew to ask for int8. @@ -120,6 +127,7 @@ fn run(cli: Cli) -> Result<(), Box> { "embedding_dim": dim, "wal_enabled": wal, "quantized_index": config_quantized, + "float16": float16, "count": mem.count(), }); println!("{}", serde_json::to_string_pretty(&j)?); diff --git a/crates/clawhdf5-format/Cargo.toml b/crates/clawhdf5-format/Cargo.toml index 8567a7f..a88b734 100644 --- a/crates/clawhdf5-format/Cargo.toml +++ b/crates/clawhdf5-format/Cargo.toml @@ -24,6 +24,7 @@ libaec-sys = { path = "../libaec-sys", version = "0.1", optional = true } pco = { version = "1.0", optional = true } [dev-dependencies] +half = { workspace = true } serde_json = "1" criterion = { workspace = true } clawhdf5-derive = { path = "../clawhdf5-derive", version = "2.7.0" } diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 0e6773f..e3c1bbd 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -1076,6 +1076,21 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr ) { return Ok(native_le_to_vec::(raw, count)); } + // Little-endian half precision (numpy float16): widen directly. + if matches!( + datatype, + Datatype::FloatingPoint { + size: 2, + byte_order: DatatypeByteOrder::LittleEndian, + .. + } + ) { + let (halves, _) = raw[..count * 2].as_chunks::<2>(); + return Ok(halves + .iter() + .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b))) + .collect()); + } let order = get_byte_order(datatype); let mut result = Vec::with_capacity(count); @@ -1622,36 +1637,7 @@ fn read_f16_bytes(bytes: &[u8], order: &DatatypeByteOrder) -> f32 { f16_bits_to_f32(u16::from_le_bytes(buf)) } -/// Convert the bit pattern of an IEEE-754 half (binary16) to an `f32`. -fn f16_bits_to_f32(h: u16) -> f32 { - let h = h as u32; - let sign = (h & 0x8000) << 16; - let exp = (h >> 10) & 0x1f; - let mant = h & 0x3ff; - let bits = if exp == 0 { - if mant == 0 { - sign // signed zero - } else { - // Subnormal: normalize into an f32 normal. - let mut e: i32 = -1; - let mut m = mant; - loop { - e += 1; - m <<= 1; - if m & 0x400 != 0 { - break; - } - } - let m = m & 0x3ff; - sign | (((127 - 15 - e) as u32) << 23) | (m << 13) - } - } else if exp == 0x1f { - sign | 0x7f80_0000 | (mant << 13) // inf / NaN - } else { - sign | ((exp + (127 - 15)) << 23) | (mant << 13) - }; - f32::from_bits(bits) -} +use crate::float16::f16_bits_to_f32; fn read_f32_bytes(bytes: &[u8], order: &DatatypeByteOrder) -> f32 { let mut buf = [0u8; 4]; diff --git a/crates/clawhdf5-format/src/datatype.rs b/crates/clawhdf5-format/src/datatype.rs index b59652f..4a1ba26 100644 --- a/crates/clawhdf5-format/src/datatype.rs +++ b/crates/clawhdf5-format/src/datatype.rs @@ -829,8 +829,12 @@ mod tests { // The HDF5 library rejects a float whose sign position is not inside // its precision; this was hard-coded to 63, so every f32 we wrote was // unreadable by h5py. Byte 2 of the message is the sign position. - use crate::type_builders::{make_f32_type, make_f64_type}; - for (dt, sign) in [(make_f32_type(), 31), (make_f64_type(), 63)] { + use crate::type_builders::{make_f16_type, make_f32_type, make_f64_type}; + for (dt, sign) in [ + (make_f16_type(), 15), + (make_f32_type(), 31), + (make_f64_type(), 63), + ] { let bytes = dt.serialize(); assert_eq!(bytes[2], sign, "{dt:?}"); let (parsed, _) = Datatype::parse(&bytes).unwrap(); diff --git a/crates/clawhdf5-format/src/float16.rs b/crates/clawhdf5-format/src/float16.rs new file mode 100644 index 0000000..1039035 --- /dev/null +++ b/crates/clawhdf5-format/src/float16.rs @@ -0,0 +1,155 @@ +//! IEEE-754 half precision (binary16) conversions. +//! +//! Pure integer bit manipulation, so it works under `no_std` and needs no +//! `libm`. The writer ([`crate::type_builders::DatasetBuilder::with_f16_data`]), +//! the reader and `clawhdf5-agent`'s half-precision embedding store all use +//! these two functions, so a value rounded in memory is bit-for-bit the value +//! that reads back from the file. + +/// Largest finite half-precision value. Anything larger in magnitude rounds +/// to infinity. +pub const F16_MAX: f32 = 65504.0; + +/// Convert an `f32` to the bit pattern of the nearest half-precision value, +/// rounding ties to even (the IEEE default, and what numpy and the `half` +/// crate do). +/// +/// Values beyond ±[`F16_MAX`] become ±infinity, values too small for a +/// subnormal become signed zero, and NaN stays NaN (quiet, payload +/// truncated). +pub fn f32_to_f16_bits(value: f32) -> u16 { + let x = value.to_bits(); + let sign = (x >> 16) & 0x8000; + let exp = x & 0x7F80_0000; + let man = x & 0x007F_FFFF; + + // Infinity and NaN. + if exp == 0x7F80_0000 { + let quiet_nan = if man == 0 { 0 } else { 0x0200 }; + return (sign | 0x7C00 | quiet_nan | (man >> 13)) as u16; + } + + let half_exp = ((exp >> 23) as i32) - 127 + 15; + + // Too large: infinity. + if half_exp >= 0x1F { + return (sign | 0x7C00) as u16; + } + + // Subnormal half, or zero. + if half_exp <= 0 { + if 14 - half_exp > 24 { + return sign as u16; + } + let man = man | 0x0080_0000; // implicit leading bit + let shift = (14 - half_exp) as u32; + let mut half_man = man >> shift; + let round_bit = 1u32 << (shift - 1); + // Round half to even: up if above half, or exactly half and odd. + if (man & round_bit) != 0 && (man & (3 * round_bit - 1)) != 0 { + half_man += 1; + } + return (sign | half_man) as u16; + } + + // Normal half. A mantissa carry correctly rolls into the exponent (and + // from the largest finite value into infinity). + let half = sign | ((half_exp as u32) << 10) | (man >> 13); + let round_bit = 0x0000_1000; + if (man & round_bit) != 0 && (man & (3 * round_bit - 1)) != 0 { + (half + 1) as u16 + } else { + half as u16 + } +} + +/// Convert the bit pattern of a half-precision value to `f32` (exact: every +/// half value is representable as an `f32`). +pub fn f16_bits_to_f32(h: u16) -> f32 { + let h = h as u32; + let sign = (h & 0x8000) << 16; + let exp = (h >> 10) & 0x1f; + let mant = h & 0x3ff; + let bits = if exp == 0 { + if mant == 0 { + sign // signed zero + } else { + // Subnormal: normalize into an f32 normal. + let mut e: i32 = -1; + let mut m = mant; + loop { + e += 1; + m <<= 1; + if m & 0x400 != 0 { + break; + } + } + let m = m & 0x3ff; + sign | (((127 - 15 - e) as u32) << 23) | (m << 13) + } + } else if exp == 0x1f { + sign | 0x7f80_0000 | (mant << 13) // inf / NaN + } else { + sign | ((exp + 127 - 15) << 23) | (mant << 13) + }; + f32::from_bits(bits) +} + +/// Round an `f32` to the nearest half-precision value, returned as `f32`. +pub fn round_to_f16(value: f32) -> f32 { + f16_bits_to_f32(f32_to_f16_bits(value)) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn every_half_value_round_trips() { + for bits in 0..=u16::MAX { + let v = f16_bits_to_f32(bits); + if v.is_nan() { + assert!(f16_bits_to_f32(f32_to_f16_bits(v)).is_nan(), "{bits:#06x}"); + } else { + assert_eq!(f32_to_f16_bits(v), bits, "{bits:#06x} -> {v}"); + } + } + } + + #[test] + fn matches_the_half_crate() { + // Every 257th f32 bit pattern (~16.7M values) covers every exponent, + // the subnormal range, both signs, ties and the overflow boundary. + let mut bits: u32 = 0; + loop { + let v = f32::from_bits(bits); + let ours = f32_to_f16_bits(v); + let theirs = half::f16::from_f32(v); + if v.is_nan() { + assert!(theirs.is_nan() && f16_bits_to_f32(ours).is_nan()); + } else { + assert_eq!(ours, theirs.to_bits(), "{bits:#010x} ({v:e})"); + assert_eq!(f16_bits_to_f32(ours).to_bits(), theirs.to_f32().to_bits()); + } + match bits.checked_add(257) { + Some(b) => bits = b, + None => break, + } + } + } + + #[test] + fn rounds_ties_to_even_and_saturates_to_infinity() { + // 1 + 2^-11 is exactly halfway between 1.0 and the next half (1 + 2^-10). + assert_eq!(round_to_f16(1.0 + 2f32.powi(-11)), 1.0); + assert_eq!( + round_to_f16(1.0 + 3.0 * 2f32.powi(-11)), + 1.0 + 2.0 * 2f32.powi(-10) + ); + assert_eq!(round_to_f16(F16_MAX), F16_MAX); + assert_eq!(round_to_f16(65520.0), f32::INFINITY); // halfway to 2^16 rounds up + assert_eq!(round_to_f16(-1e9), f32::NEG_INFINITY); + assert_eq!(round_to_f16(1e-9).to_bits(), 0); + assert_eq!(round_to_f16(-1e-9).to_bits(), (-0.0f32).to_bits()); + } +} diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index 83e9740..54bd326 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -72,6 +72,7 @@ pub mod filter_pipeline; pub mod filters; mod filters_szip; pub mod fixed_array; +pub mod float16; pub mod fractal_heap; pub mod global_heap; pub mod group_info; diff --git a/crates/clawhdf5-format/src/type_builders.rs b/crates/clawhdf5-format/src/type_builders.rs index 270cd48..e698899 100644 --- a/crates/clawhdf5-format/src/type_builders.rs +++ b/crates/clawhdf5-format/src/type_builders.rs @@ -56,6 +56,21 @@ pub fn make_f64_type() -> Datatype { } } +/// IEEE-754 half precision (binary16), little-endian — numpy's `float16`. +pub fn make_f16_type() -> Datatype { + Datatype::FloatingPoint { + size: 2, + byte_order: DatatypeByteOrder::LittleEndian, + bit_offset: 0, + bit_precision: 16, + exponent_location: 10, + exponent_size: 5, + mantissa_location: 0, + mantissa_size: 10, + exponent_bias: 15, + } +} + pub fn make_f32_type() -> Datatype { Datatype::FloatingPoint { size: 4, @@ -478,6 +493,24 @@ impl DatasetBuilder { self } + /// Store `data` as IEEE half precision (numpy `float16`), rounding each + /// value to the nearest half ([`crate::float16::f32_to_f16_bits`]). + /// Half the bytes of [`Self::with_f32_data`], at about three significant + /// decimal digits; values beyond ±65504 become ±infinity. Reading it back + /// with `read_f32` yields the rounded values exactly. + pub fn with_f16_data(&mut self, data: &[f32]) -> &mut Self { + self.datatype = Some(make_f16_type()); + let mut b = Vec::with_capacity(data.len() * 2); + for &v in data { + b.extend_from_slice(&crate::float16::f32_to_f16_bits(v).to_le_bytes()); + } + self.data = Some(b); + if self.shape.is_none() { + self.shape = Some(vec![data.len() as u64]); + } + self + } + pub fn with_i32_data(&mut self, data: &[i32]) -> &mut Self { self.datatype = Some(make_i32_type()); let mut b = Vec::with_capacity(data.len() * 4); diff --git a/crates/clawhdf5/tests/h5py_interop_tests.rs b/crates/clawhdf5/tests/h5py_interop_tests.rs index acac2dc..042940f 100644 --- a/crates/clawhdf5/tests/h5py_interop_tests.rs +++ b/crates/clawhdf5/tests/h5py_interop_tests.rs @@ -1388,6 +1388,130 @@ with h5py.File("{path_str}", "w", libver="latest") as f: } } +// --------------------------------------------------------------------------- +// Half precision (float16) in both directions +// --------------------------------------------------------------------------- + +/// Values that exercise rounding: ties, subnormals, the overflow boundary and +/// ordinary embedding-sized components. +fn f16_probe_values() -> Vec { + let mut v = vec![ + 0.0, + -0.0, + 1.0, + -1.0, + 0.5, + 1.0 + 2f32.powi(-11), + 1.0 + 3.0 * 2f32.powi(-11), + 65504.0, + 65519.0, + 65520.0, + -70000.0, + 6.0e-8, + 3.0e-8, + 1.0e-9, + 1.0e-5, + 0.1, + 0.333_333, + 1234.567, + f32::INFINITY, + f32::NEG_INFINITY, + ]; + // A deterministic spread of embedding-like values. + let mut x = 0x2545_F491u32; + for _ in 0..4000 { + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + v.push((x as f32 / u32::MAX as f32 - 0.5) * 0.4); + } + v +} + +#[test] +fn clawhdf5_writes_f16_h5py_reads() { + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("ours_f16.h5"); + let path_str = path.display().to_string(); + let values = f16_probe_values(); + + let mut fb = FileBuilder::new(); + fb.create_dataset("plain").with_f16_data(&values); + fb.create_dataset("chunked") + .with_f16_data(&values) + .with_shape(&[values.len() as u64]) + .with_chunks(&[512]) + .with_deflate(6); + fb.write(&path).unwrap(); + + // h5py must see a genuine float16 dataset, and our rounding must agree + // with numpy's own float32 -> float16 conversion bit for bit. + let input = values + .iter() + .map(|v| format!("{:?}", v.to_bits())) + .collect::>() + .join(","); + let script = format!( + r#" +import h5py, numpy as np +src = np.array([{input}], dtype=np.uint32).view(np.float32) +expected = src.astype(np.float16).view(np.uint16) +with h5py.File("{path_str}", "r") as f: + for name in ("plain", "chunked"): + d = f[name] + assert d.dtype == np.float16, (name, d.dtype) + got = d[:].view(np.uint16) + bad = np.nonzero(got != expected)[0] + assert bad.size == 0, (name, bad[:5], got[bad[:5]], expected[bad[:5]]) +print("ok") +"# + ); + assert_eq!(run_python_output(&script), "ok"); +} + +#[test] +fn h5py_writes_f16_clawhdf5_reads() { + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("h5py_f16.h5"); + let path_str = path.display().to_string(); + let values = f16_probe_values(); + let input = values + .iter() + .map(|v| format!("{:?}", v.to_bits())) + .collect::>() + .join(","); + + let script = format!( + r#" +import h5py, numpy as np +src = np.array([{input}], dtype=np.uint32).view(np.float32).astype(np.float16) +with h5py.File("{path_str}", "w") as f: + f.create_dataset("plain", data=src) + f.create_dataset("chunked", data=src, chunks=(512,), compression="gzip", shuffle=True) + f.create_dataset("big_endian", data=src.astype(">f2")) +"# + ); + run_python(&script); + + let expected: Vec = values + .iter() + .map(|&v| clawhdf5_format::float16::round_to_f16(v).to_bits()) + .collect(); + let file = File::open(&path).unwrap(); + for name in ["plain", "chunked", "big_endian"] { + let ds = file.dataset(name).unwrap(); + assert_eq!( + ds.dtype().unwrap(), + DType::Other("float16".into()), + "{name}" + ); + let got: Vec = ds.read_f32().unwrap().iter().map(|v| v.to_bits()).collect(); + assert_eq!(got, expected, "{name}"); + } +} + #[test] fn clawhdf5_writes_f32_h5py_reads() { // Every f32 dataset used to be unreadable by h5py ("sign bit position out