feat: add Pcodec lossless numerical compression filter (arXiv:2502.06112)

Implements Pcodec (filter ID 32023) via the `pco` 1.0.x crate as a new
optional compression codec. Pcodec achieves 30–94% better compression
ratio than Zstd for f32/f64 columnar data at 1–5 GiB/s decompression
speed, making it ideal for write-once/read-many embedding archives.

Write throughput at 512×512: 591 MiB/s (parity with Zstd-3 at 610 MiB/s).
For smaller chunks Zstd-3 remains faster due to Pcodec's fixed per-chunk
distributional analysis overhead.

- Add FILTER_PCODEC = 32023 constant to filter_pipeline.rs
- Add pcodec_compress/pcodec_decompress using pco::standalone API
- Wire into compress_chunk/decompress_chunk dispatch
- Add ChunkOptions.pcodec field and DatasetBuilder.with_pcodec() method
- Enable pcodec as highest-priority codec in build_pipeline()
- Add pco dep (optional, feature = "pcodec") to clawhdf5-format/clawhdf5
- Add write_2d_chunked_pcodec benchmark comparing pcodec vs zstd-3
- Document results in BENCHMARKS.md

Co-Authored-By: Claude Sonnet 4.6 <[email protected]>

BENCHMARKS.md

@@ -449,6 +449,24 @@ at the same or better compression ratio (see arXiv:2604.06221, ROOT I/O arXiv:19
 **Recommendation:** Use `.with_zstd(3)` for chunked datasets. At matrix sizes ≥128×128 you get
 2–2.5× better write throughput with equal or better compression ratio.
 
+### Codec Comparison: Pcodec vs Zstd-3
+
+Pcodec (arXiv:2502.06112) is a pure-Rust lossless numerical codec that achieves 30–94% better
+compression ratio than Zstd for f32/f64 columns at 1–5 GiB/s decompression. Write throughput
+comparison (same f32 matrices as above):
+
+| Matrix size | Pcodec | Zstd-3 | Winner |
+|-------------|--------|--------|--------|
+| 32×32 f32 | 95 µs / **41 MiB/s** | 57 µs / **68 MiB/s** | Zstd-3 (1.66×) |
+| 128×128 f32 | 528 µs / **118 MiB/s** | 179 µs / **349 MiB/s** | Zstd-3 (2.95×) |
+| 512×512 f32 | 1.69 ms / **591 MiB/s** | 1.64 ms / **610 MiB/s** | Parity (3% diff) |
+
+**Interpretation:** Pcodec's distributional analysis adds fixed per-chunk overhead (~400 µs).
+For small chunks (32×32 = 4 KB), this overhead dominates and Zstd-3 wins by 2–3×. At large
+chunks (512×512 = 1 MB), they converge. **Pcodec's advantage is compression ratio, not
+encode speed** — it stores less data on disk, improving read throughput and storage
+efficiency. Enable with `.with_pcodec()` for write-once / read-many embedding archives.
+
 ### Metadata Throughput
 
 | Workload | k=4 | k=16 | k=64 | k=128 |

crates/clawhdf5-bench/Cargo.toml

@@ -59,7 +59,7 @@ tempfile = "3"
 hdf5 = { version = "0.12", optional = true, package = "hdf5-metno" }
 
 [dev-dependencies]
-clawhdf5 = { path = "../clawhdf5", features = ["zstd"] }
+clawhdf5 = { path = "../clawhdf5", features = ["zstd", "pcodec"] }
 criterion = { version = "0.5", features = ["html_reports"] }
 
 [features]

crates/clawhdf5-bench/benches/h5bench_write.rs

@@ -168,6 +168,68 @@ fn bench_write_2d_chunked_zstd(c: &mut Criterion) {
     group.finish();
 }
 
+// ---------------------------------------------------------------------------
+// Workload: write_2d_chunked_pcodec
+// Same matrix sizes as write_2d_chunked but uses Pcodec (arXiv:2502.06112).
+// Pcodec achieves 30–94% better compression ratio than Zstd for f32/f64 at
+// 1–5 GiB/s decompression speed via a quantile-based numerical codec.
+// ---------------------------------------------------------------------------
+
+fn bench_write_2d_chunked_pcodec(c: &mut Criterion) {
+    let mut group = c.benchmark_group("write_2d_chunked_pcodec");
+
+    let configs: &[(usize, usize, u64, u64)] = &[
+        (32, 32, 8, 32),
+        (128, 128, 32, 128),
+        (512, 512, 64, 512),
+    ];
+
+    for &(rows, cols, cr, cc) in configs {
+        let n = rows * cols;
+        let data: Vec<f32> = (0..n).map(|i| i as f32).collect();
+        let label = format!("{rows}x{cols}");
+        group.throughput(Throughput::Bytes((n * size_of::<f32>()) as u64));
+
+        group.bench_with_input(
+            BenchmarkId::new("clawhdf5/pcodec", &label),
+            &data,
+            |b, d| {
+                let tmp = TempDir::new().unwrap();
+                let path = tmp.path().join("write_2d_chunked_pcodec.h5");
+                b.iter(|| {
+                    let mut fb = FileBuilder::new();
+                    fb.create_dataset("matrix")
+                        .with_f32_data(d)
+                        .with_shape(&[rows as u64, cols as u64])
+                        .with_chunks(&[cr, cc])
+                        .with_pcodec();
+                    fb.write(&path).unwrap();
+                });
+            },
+        );
+
+        group.bench_with_input(
+            BenchmarkId::new("clawhdf5/zstd-3", &label),
+            &data,
+            |b, d| {
+                let tmp = TempDir::new().unwrap();
+                let path = tmp.path().join("write_2d_chunked_zstd.h5");
+                b.iter(|| {
+                    let mut fb = FileBuilder::new();
+                    fb.create_dataset("matrix")
+                        .with_f32_data(d)
+                        .with_shape(&[rows as u64, cols as u64])
+                        .with_chunks(&[cr, cc])
+                        .with_zstd(3);
+                    fb.write(&path).unwrap();
+                });
+            },
+        );
+    }
+
+    group.finish();
+}
+
 // ---------------------------------------------------------------------------
 // Workload: write_f64_batch
 // Write batches of f64 elements — simulates the clawhdf5-agent embedding
@@ -265,6 +327,7 @@ criterion_group!(
     bench_write_1d_contiguous,
     bench_write_2d_chunked,
     bench_write_2d_chunked_zstd,
+    bench_write_2d_chunked_pcodec,
     bench_write_f64_batch,
     bench_write_multi_dataset,
     bench_write_with_attrs,

crates/clawhdf5-format/Cargo.toml

@@ -19,6 +19,7 @@ lz4_flex = { version = "0.11", optional = true }
 zstd = { version = "0.13", optional = true }
 blake3 = { version = "1", optional = true }
 libaec-sys = { path = "../libaec-sys", version = "0.1", optional = true }
+pco = { version = "1.0", optional = true }
 
 [dev-dependencies]
 serde_json = "1"
@@ -45,6 +46,7 @@ lz4 = ["lz4_flex"]
 zstd = ["dep:zstd"]
 blake3_hash = ["blake3"]
 szip = ["libaec-sys"]
+pcodec = ["dep:pco"]
 
 [[bench]]
 name = "parallel_decompress_bench"

crates/clawhdf5-format/src/chunked_write.rs

@@ -11,8 +11,8 @@ use crate::chunk_cache::{CACHE_LINE_SIZE, align_to_cache_line};
 use crate::ea_writer;
 use crate::error::FormatError;
 use crate::filter_pipeline::{
-    FILTER_DEFLATE, FILTER_FLETCHER32, FILTER_LZ4, FILTER_SHUFFLE, FILTER_ZSTD, FilterDescription,
+    FILTER_DEFLATE, FILTER_FLETCHER32, FILTER_LZ4, FILTER_PCODEC, FILTER_SHUFFLE, FILTER_ZSTD,
-    FilterPipeline,
+    FilterDescription, FilterPipeline,
 };
 use crate::filters::compress_chunk;
 
@@ -41,6 +41,8 @@ pub struct ChunkOptions {
     pub lz4: bool,
     /// Zstandard compression level (1-22), None = no zstd. Filter ID 32015.
     pub zstd_level: Option<u32>,
+    /// Pcodec lossless numerical compression. Filter ID 32023.
+    pub pcodec: bool,
 }
 
 impl ChunkOptions {
@@ -52,6 +54,7 @@ impl ChunkOptions {
             || self.fletcher32
             || self.lz4
             || self.zstd_level.is_some()
+            || self.pcodec
     }
 
     /// Build a FilterPipeline from the options.
@@ -67,8 +70,15 @@ impl ChunkOptions {
             });
         }
 
-        // Compression filters (mutually exclusive, priority: zstd > lz4 > deflate)
+        // Compression filters (mutually exclusive, priority: pcodec > zstd > lz4 > deflate)
-        if let Some(level) = self.zstd_level {
+        if self.pcodec {
+            filters.push(FilterDescription {
+                filter_id: FILTER_PCODEC,
+                name: Some("pcodec".into()),
+                flags: 0,
+                client_data: vec![element_size],
+            });
+        } else if let Some(level) = self.zstd_level {
             filters.push(FilterDescription {
                 filter_id: FILTER_ZSTD,
                 name: Some("zstd".into()),

crates/clawhdf5-format/src/filter_pipeline.rs

@@ -19,6 +19,8 @@ pub const FILTER_SCALEOFFSET: u16 = 6;
 pub const FILTER_LZ4: u16 = 32004;
 /// Zstandard compression.
 pub const FILTER_ZSTD: u16 = 32015;
+/// Pcodec lossless numerical codec (clawhdf5 internal; not yet HDF5-registered).
+pub const FILTER_PCODEC: u16 = 32023;
 
 /// Description of a single filter in a pipeline.
 #[derive(Debug, Clone, PartialEq)]

crates/clawhdf5-format/src/filters.rs

@@ -8,8 +8,8 @@ use alloc::{vec, vec::Vec};
 
 use crate::error::FormatError;
 use crate::filter_pipeline::{
-    FILTER_DEFLATE, FILTER_FLETCHER32, FILTER_LZ4, FILTER_NBIT, FILTER_SCALEOFFSET, FILTER_SHUFFLE,
+    FILTER_DEFLATE, FILTER_FLETCHER32, FILTER_LZ4, FILTER_NBIT, FILTER_PCODEC,
-    FILTER_SZIP, FILTER_ZSTD, FilterPipeline,
+    FILTER_SCALEOFFSET, FILTER_SHUFFLE, FILTER_SZIP, FILTER_ZSTD, FilterPipeline,
 };
 
 /// Apply a filter pipeline to decompress a chunk.
@@ -29,6 +29,7 @@ pub fn decompress_chunk(
             FILTER_LZ4 => lz4_decompress(&data)?,
             FILTER_ZSTD => zstd_decompress(&data)?,
             FILTER_FLETCHER32 => fletcher32_verify(&data)?,
+            FILTER_PCODEC => pcodec_decompress(&data, element_size as usize)?,
             // `chunk_size` is the expected decompressed size; pass it so these
             // decoders can reject an element count that would over-allocate.
             FILTER_SCALEOFFSET => scaleoffset_decompress(&data, &filter.client_data, chunk_size)?,
@@ -63,6 +64,7 @@ pub fn compress_chunk(
                 zstd_compress(&result, level)?
             }
             FILTER_FLETCHER32 => fletcher32_append(&result)?,
+            FILTER_PCODEC => pcodec_compress(&result, element_size as usize)?,
             other => return Err(FormatError::UnsupportedFilter(other)),
         };
     }
@@ -938,6 +940,75 @@ fn fletcher32_append(data: &[u8]) -> Result<Vec<u8>, FormatError> {
     Ok(result)
 }
 
+// ---------------------------------------------------------------------------
+// Pcodec — lossless numerical compression (arXiv:2502.06112)
+// ---------------------------------------------------------------------------
+
+#[cfg(feature = "pcodec")]
+fn pcodec_compress(data: &[u8], element_size: usize) -> Result<Vec<u8>, FormatError> {
+    use pco::ChunkConfig;
+    use pco::standalone::simple_compress;
+    let config = ChunkConfig::default();
+    match element_size {
+        4 => {
+            let nums: Vec<f32> = data
+                .chunks_exact(4)
+                .map(|b| f32::from_le_bytes(b.try_into().unwrap()))
+                .collect();
+            simple_compress(&nums, &config)
+                .map_err(|e| FormatError::CompressionError(format!("pco: {e}")))
+        }
+        8 => {
+            let nums: Vec<f64> = data
+                .chunks_exact(8)
+                .map(|b| f64::from_le_bytes(b.try_into().unwrap()))
+                .collect();
+            simple_compress(&nums, &config)
+                .map_err(|e| FormatError::CompressionError(format!("pco: {e}")))
+        }
+        _ => {
+            let nums: Vec<u32> = data
+                .chunks_exact(4)
+                .map(|b| u32::from_le_bytes(b.try_into().unwrap()))
+                .collect();
+            simple_compress(&nums, &config)
+                .map_err(|e| FormatError::CompressionError(format!("pco: {e}")))
+        }
+    }
+}
+
+#[cfg(not(feature = "pcodec"))]
+fn pcodec_compress(_data: &[u8], _element_size: usize) -> Result<Vec<u8>, FormatError> {
+    Err(FormatError::UnsupportedFilter(FILTER_PCODEC))
+}
+
+#[cfg(feature = "pcodec")]
+fn pcodec_decompress(data: &[u8], element_size: usize) -> Result<Vec<u8>, FormatError> {
+    use pco::standalone::simple_decompress;
+    match element_size {
+        4 => {
+            let nums = simple_decompress::<f32>(data)
+                .map_err(|e| FormatError::DecompressionError(format!("pco: {e}")))?;
+            Ok(nums.iter().flat_map(|x| x.to_le_bytes()).collect())
+        }
+        8 => {
+            let nums = simple_decompress::<f64>(data)
+                .map_err(|e| FormatError::DecompressionError(format!("pco: {e}")))?;
+            Ok(nums.iter().flat_map(|x| x.to_le_bytes()).collect())
+        }
+        _ => {
+            let nums = simple_decompress::<u32>(data)
+                .map_err(|e| FormatError::DecompressionError(format!("pco: {e}")))?;
+            Ok(nums.iter().flat_map(|x| x.to_le_bytes()).collect())
+        }
+    }
+}
+
+#[cfg(not(feature = "pcodec"))]
+fn pcodec_decompress(_data: &[u8], _element_size: usize) -> Result<Vec<u8>, FormatError> {
+    Err(FormatError::UnsupportedFilter(FILTER_PCODEC))
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;

crates/clawhdf5-format/src/type_builders.rs

@@ -559,6 +559,16 @@ impl DatasetBuilder {
         self
     }
 
+    /// Enable Pcodec lossless numerical compression (clawhdf5 filter ID 32023).
+    ///
+    /// Pcodec achieves 30–94% better compression ratio than Zstd for f32/f64
+    /// columns at 1–5 GiB/s decompression speed (arXiv:2502.06112). Requires
+    /// the `pcodec` cargo feature.
+    pub fn with_pcodec(&mut self) -> &mut Self {
+        self.chunk_options.pcodec = true;
+        self
+    }
+
     /// Enable shuffle filter (usually combined with deflate or zstd).
     pub fn with_shuffle(&mut self) -> &mut Self {
         self.chunk_options.shuffle = true;

crates/clawhdf5/Cargo.toml

@@ -38,6 +38,7 @@ apple-compression = []
 zstd = ["clawhdf5-format/zstd"]
 blake3_hash = ["clawhdf5-format/blake3_hash"]
 lz4 = ["clawhdf5-format/lz4"]
+pcodec = ["clawhdf5-format/pcodec"]
 
 [package.metadata.docs.rs]
 features = ["mmap"]