ML feature substrate: bounded, reproducible per-locus feature stream (rosalind features)

README.md

@@ -121,6 +121,23 @@ Drop the `rosalind-budget` Action into any pipeline to make a declared memory bu
 
 It runs `plan` (predicts the peak), then `variants --index --enforce` (honors the budget), and uploads the BLAKE3 receipt as a build artifact. This is the one thing a `--max-mem` flag on another caller can't give you: a portable, declarative, **verifiable** memory budget that fails a stranger's build loudly — the contract, enforced where your pipeline already lives. (Available once a release is published; see Quickstart.)
 
+## A reproducible feature substrate for ML
+
+The same bounded streaming engine that calls variants can emit **per-locus features** instead — one tabular row per callable position, ready for a model:
+
+```sh
+rosalind features --index ref.idx --alignments sorted.bam -o features.tsv
+# columns: contig, pos, ref, depth, raw_depth, A/C/G/T counts,
+#          per-allele fwd/rev strand counts, mean base-qual, mean mapq
+```
+
+```python
+import pandas as pd
+df = pd.read_csv("features.tsv", sep="\t")   # one line; ready for sklearn/PyTorch/JAX
+```
+
+Two properties no other pileup gives you together: it is **bounded** (the whole-genome table streams to disk; peak memory tracks coverage, not genome size — a 1 Mbp toy genome's ~983k-row table is produced in ~6 MiB), and it is **byte-identical run-to-run**, with a BLAKE3 receipt over the output. That means **bit-reproducible training inputs**: hash your feature file, and you can prove this quarter's model saw exactly the same data as last quarter's. `features` honors the same `plan`/`--enforce`/`verify` memory contract as `variants`. *(TSV today; an Arrow/Parquet egress and a zero-copy `pyarrow` Python binding are on the roadmap.)*
+
 ## Roadmap
 
 The core primitive is a streaming, CIGAR-aware pileup column stream; variant calling and custom plugins consume it. Performance work deliberately *follows* the unique capability — the target user needs "it fits and is predictable" before "it's fastest."

docs/superpowers/specs/2026-06-02-ml-feature-substrate-design.md

@@ -0,0 +1,110 @@
+# ML Feature Substrate — Egress Core (design)
+
+**Status:** DESIGN SPEC — 2026-06-02. **Branch:** `rosalind/ml-feature-substrate` (off `main`
+`9506d45`). From the reflection audit (angle 3): the kernel already produces a bounded, deterministic
+per-locus feature stream (`PileupColumn`); it is gated from ML builders only by the lack of a feature
+**egress** (the `on_row` sink emits germline calls, not features) and a dead-end Python stub.
+
+## 1. Goal
+
+Expose the kernel as a **bounded, deterministic, byte-identical per-locus feature stream** with a
+verifiable hash receipt — landing the novel claim no other tool makes: **byte-identical features →
+bit-reproducible training inputs.** Reuse the shipped, accuracy-validated, memory-contracted call path.
+
+## 2. Scope (this increment)
+
+The **durable, fully-testable core**: a `rosalind features` CLI + a `call::features` egress module
+streaming a per-locus **TSV** under the same memory contract (`plan`/`--enforce`/`verify`), with a
+BLAKE3 receipt. **TSV-first** (no new deps, universally readable by pandas/polars/R, trivially
+byte-deterministic) — which fully delivers the reproducibility claim; the receipt hashes the feature
+file. **Deferred to a follow-up:** Arrow/Parquet egress (efficiency/zero-copy), the real `pyarrow`
+Python boundary (replacing the stub), and a reference model demo.
+
+## 3. What exists (reuse)
+
+- `PileupColumn { locus, ref_base, raw_depth, obs }` + `depth()`, `allele_counts() -> [u32;4]`
+  (`[A,C,G,T]`), `strand_counts() -> [[u32;2];4]` (`[allele][0=fwd,1=rev]`); `Obs { allele, base_qual,
+  mapq, reverse }`. The engine yields one column per **callable** position (obs non-empty), in
+  deterministic canonical-obs order (the keystone fix).
+- `call_germline_whole_genome` (`src/call/whole_genome.rs`) — the bounded per-contig driver
+  (`PerContig` partition + per-contig `decode_window` + `PileupEngine`), returning
+  `(WorkingSet, SkipCounts)`. Make `PerContig` `pub(crate)` and reuse it.
+- The memory contract: `estimate_variants_working_set` / `--enforce` / `RunManifest` / `peak_rss_bytes`.
+
+## 4. Deliverables
+
+### 4a. `src/call/features.rs` — the egress module
+
+- **`FeatureRow`** (one per callable locus), TSV columns (tab-separated, `\n`-terminated):
+  `contig  pos  ref  depth  raw_depth  a  c  g  t  a_fwd  a_rev  c_fwd  c_rev  g_fwd  g_rev  t_fwd
+  t_rev  mean_bq  mean_mapq`
+  where `pos` is **1-based** (matches VCF POS = `locus.pos.0 + 1`), `contig` is the contig **name**,
+  `ref` is the ASCII ref base, `depth` = callable obs, `raw_depth` = covering reads, the 4 `a/c/g/t`
+  are `allele_counts`, the 8 strand columns are `strand_counts`, and `mean_bq`/`mean_mapq` are the
+  integer means over `obs` formatted `{:.2}` (byte-stable). Header line written once:
+  `#contig\tpos\tref\t…`.
+- **`feature_row_fields(col: &PileupColumn, contig_name: &str) -> impl Iterator<…>`** (or a
+  `write_feature_row<W: Write>(w, contig_name, col)` + `write_feature_header<W: Write>(w)`), mirroring
+  `src/io/vcf.rs`'s streaming writer shape. Means: `sum(base_qual as u64)/depth` etc. → `{:.2}`.
+- **`stream_features_region<S: ReadSource>(source, reference, contig, region, pileup_params, on_row:
+  &mut dyn FnMut(&PileupColumn) -> Result<(), CoreError>) -> Result<(WorkingSet, SkipCounts),
+  CoreError>`** — drive `PileupEngine`, call `on_row` per emitted column, track max working set, read
+  `skip_counts()` after. (No buffer; bounded.)
+- **`stream_features_whole_genome<S>(source, ref_view, contigs, pileup_params, on_row) ->
+  Result<(WorkingSet, SkipCounts), CoreError>`** — per-contig loop reusing `PerContig` + per-contig
+  `decode_window`, summing `SkipCounts` (mirrors `call_germline_whole_genome`). The sink gets
+  `(&PileupColumn, contig_name)` — pass the contig name through (or the sink resolves it).
+- **Unit tests:** `feature_row_from_a_known_column` (exact field values for a hand-built column);
+  `feature_stream_is_bounded_by_coverage_not_read_count` (working set flat in read count, like the
+  caller); `feature_rows_match_per_contig`.
+
+### 4b. `rosalind features` CLI — `src/main.rs`
+
+`Commands::Features { index, alignments, mapq_threshold, max_depth (1000), max_read_len (250),
+memory_budget_mb, enforce, output, manifest }` → `run_features`, structurally mirroring
+`run_variants_index`:
+- `StreamingBamSource` + `stream_features_whole_genome`, streaming each row to a `BufWriter` (TSV
+  header once, then one row per column) — no genome-wide buffer.
+- The **same** `--enforce` gate (exit 3 refuse / exit 4 breach; the working-set model is identical —
+  it is the same pileup engine), the same realized-peak receipt, the same `over_max_depth`/skip
+  surfacing. Receipt `subcommand = "features"`, plus a `feature_rows` param (count emitted).
+- `plan --index` already predicts this peak (same engine) — no change needed.
+
+### 4c. Tests + docs
+
+- **Integration test** (`tests/features.rs`): run `rosalind features` on the `build_sorted_bam_fixture`
+  (or a local fixture); assert the header + ≥1 row with the expected columns; run **twice** and assert
+  the two `features.tsv` are **byte-identical** (the reproducibility claim); assert the receipt records
+  `feature_rows` and that `verify` passes.
+- **README/CONTRACT** short section: "a bounded, deterministic, byte-identical per-locus feature stream
+  — bit-reproducible ML training inputs with a verifiable receipt," with the one-liner
+  `rosalind features --index ref.idx --alignments sorted.bam -o features.tsv` and a note that pandas
+  reads it in one line. Mark Arrow/pyarrow/model as the roadmap follow-up.
+
+## 5. Cross-cutting
+
+- **Determinism is the headline.** `features.tsv` must be byte-identical run-to-run (canonical obs
+  order + integer-sum means formatted `{:.2}` + deterministic column stream). The byte-identical
+  integration test is the proof.
+- **Bounded memory.** Rows stream to disk; one contig's reference resident; active set depth-capped —
+  the same guarantee as the caller, inherited by construction.
+- Per-item commit; 0 warnings (debug+release); fmt clean; full `cargo test` green.
+
+## 6. Out of scope (follow-up increment)
+
+- **Arrow/Parquet** egress (zero-copy, columnar efficiency) behind a feature flag.
+- **`pyarrow` Python boundary** replacing the `python_bindings` stub (`rosalind.features(index, bam,
+  region, budget_mb)` yielding RecordBatches) — needs maturin/pyO3 build verification.
+- **Reference model** demo (a notebook/script training a tiny SNV classifier on the stream + showing
+  the receipts match → bit-reproducible eval).
+
+## 7. Self-review
+
+- **Coverage:** egress module (4a), CLI (4b), tests + docs (4c). ✓
+- **Type consistency:** the feature sink is `&mut dyn FnMut(&PileupColumn) -> Result<(), CoreError>`;
+  `stream_features_*` return `(WorkingSet, SkipCounts)` like the germline drivers; `PerContig` becomes
+  `pub(crate)`. ✓
+- **No placeholders:** the TSV schema, means, and 1-based pos are concrete. ✓
+- **Reproducibility:** integer means `{:.2}` + canonical obs order + byte-identical test → the claim is
+  proven, not asserted. ✓
+- **Scope:** TSV core only; Arrow/pyarrow/model explicitly deferred. ✓

src/call/features.rs

@@ -0,0 +1,231 @@
+//! Per-locus FEATURE egress: the bounded, deterministic `PileupColumn` stream
+//! exposed as a tabular ML feature source. Same streaming engine, same memory
+//! contract as germline calling — but every callable locus is emitted as a
+//! feature row (not just variant sites). The output is byte-identical run-to-run
+//! (canonical obs order + integer-mean formatting), so a BLAKE3 receipt over it
+//! gives **bit-reproducible training inputs**.
+
+use std::io::{self, Write};
+use std::ops::Range;
+use std::sync::Arc;
+
+use crate::call::whole_genome::PerContig;
+use crate::core::{AlignedRead, ContigSet, CoreError, WorkingSet};
+use crate::genomics::ReferenceView;
+use crate::pileup::{PileupColumn, PileupEngine, PileupParams, ReadSource, SkipCounts};
+
+/// The tab-separated header for the per-locus feature table (written once).
+pub const FEATURE_HEADER: &str = "#contig\tpos\tref\tdepth\traw_depth\ta\tc\tg\tt\t\
+a_fwd\ta_rev\tc_fwd\tc_rev\tg_fwd\tg_rev\tt_fwd\tt_rev\tmean_bq\tmean_mapq";
+
+/// Write the feature header line.
+pub fn write_feature_header<W: Write>(out: &mut W) -> io::Result<()> {
+    writeln!(out, "{FEATURE_HEADER}")
+}
+
+/// Write one feature row for a callable pileup column. `pos` is 1-based (VCF POS).
+/// Means are integer sums over observations divided by callable depth, formatted
+/// to 2 decimals so the row is byte-stable across runs.
+pub fn write_feature_row<W: Write>(
+    out: &mut W,
+    contig_name: &str,
+    col: &PileupColumn,
+) -> io::Result<()> {
+    let depth = col.depth();
+    let ac = col.allele_counts();
+    let sc = col.strand_counts();
+    let (sum_bq, sum_mapq) = col.obs.iter().fold((0u64, 0u64), |(b, m), o| {
+        (b + o.base_qual as u64, m + o.mapq as u64)
+    });
+    let (mean_bq, mean_mapq) = if depth > 0 {
+        (sum_bq as f64 / depth as f64, sum_mapq as f64 / depth as f64)
+    } else {
+        (0.0, 0.0)
+    };
+    writeln!(
+        out,
+        "{contig}\t{pos}\t{refb}\t{depth}\t{raw}\t\
+{a}\t{c}\t{g}\t{t}\t\
+{afwd}\t{arev}\t{cfwd}\t{crev}\t{gfwd}\t{grev}\t{tfwd}\t{trev}\t\
+{mean_bq:.2}\t{mean_mapq:.2}",
+        contig = contig_name,
+        pos = col.locus.pos.0 as u64 + 1,
+        refb = col.ref_base as char,
+        depth = depth,
+        raw = col.raw_depth,
+        a = ac[0],
+        c = ac[1],
+        g = ac[2],
+        t = ac[3],
+        afwd = sc[0][0],
+        arev = sc[0][1],
+        cfwd = sc[1][0],
+        crev = sc[1][1],
+        gfwd = sc[2][0],
+        grev = sc[2][1],
+        tfwd = sc[3][0],
+        trev = sc[3][1],
+    )
+}
+
+/// Stream feature columns over `region` of `contig` to a sink, returning the max
+/// pileup working set + skip counts. The sink receives each callable column in
+/// ascending position order; no genome-wide buffer accumulates.
+pub fn stream_features_region<S: ReadSource>(
+    source: S,
+    reference: Arc<[u8]>,
+    contig: u32,
+    region: Range<u32>,
+    pileup_params: PileupParams,
+    on_row: &mut dyn FnMut(&PileupColumn) -> Result<(), CoreError>,
+) -> Result<(WorkingSet, SkipCounts), CoreError> {
+    let mut engine = PileupEngine::new(source, reference, contig, region, pileup_params);
+    let mut max_ws = WorkingSet { bytes: 0 };
+    while let Some(column) = engine.next() {
+        let column = column?;
+        let ws = engine.current_working_set();
+        if ws.bytes > max_ws.bytes {
+            max_ws = ws;
+        }
+        on_row(&column)?;
+    }
+    Ok((max_ws, engine.skip_counts()))
+}
+
+/// Stream feature columns across every contig (in id order), reading each contig's
+/// reference from `ref_view`. The sink receives `(column, contig_name)`. Peak
+/// memory ≈ the largest contig's reference + the pileup working set, independent
+/// of input size. `source` MUST be `(contig, pos)`-ordered (a `StreamingBamSource`
+/// guards this). Mirrors `call_germline_whole_genome`.
+pub fn stream_features_whole_genome<S: ReadSource>(
+    mut source: S,
+    ref_view: &ReferenceView,
+    contigs: &ContigSet,
+    pileup_params: PileupParams,
+    on_row: &mut dyn FnMut(&PileupColumn, &str) -> Result<(), CoreError>,
+) -> Result<(WorkingSet, SkipCounts), CoreError> {
+    let mut max_ws = WorkingSet { bytes: 0 };
+    let mut skips = SkipCounts::default();
+    let mut peeked: Option<AlignedRead> = None;
+
+    for c in contigs.iter() {
+        let start = c.global_offset as usize;
+        let end = c.global_offset as usize + c.length as usize;
+        let mut decoded = Vec::new();
+        ref_view.decode_window(start, end, &mut decoded);
+        let reference: Arc<[u8]> = Arc::from(decoded);
+
+        let per = PerContig {
+            source: &mut source,
+            contig: c.id,
+            peeked: &mut peeked,
+        };
+        let region: Range<u32> = 0..c.length;
+        let (ws, contig_skips) = stream_features_region(
+            per,
+            reference,
+            c.id,
+            region,
+            pileup_params.clone(),
+            &mut |col| on_row(col, &c.name),
+        )?;
+        if ws.bytes > max_ws.bytes {
+            max_ws = ws;
+        }
+        skips.accumulate(&contig_skips);
+    }
+
+    Ok((max_ws, skips))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::core::{CigarOp, CigarOpKind, Position, SamFlags};
+    use crate::pileup::SliceSource;
+
+    fn mread(pos: u32, seq: &[u8], reverse: bool) -> AlignedRead {
+        let flags = if reverse {
+            SamFlags(SamFlags::REVERSE)
+        } else {
+            SamFlags::default()
+        };
+        AlignedRead {
+            contig: 0,
+            pos: Position(pos),
+            mapq: 60,
+            flags,
+            cigar: vec![CigarOp::new(CigarOpKind::Match, seq.len() as u32)],
+            seq: Arc::from(seq.to_vec().into_boxed_slice()),
+            qual: Arc::from(vec![30u8; seq.len()].into_boxed_slice()),
+        }
+    }
+
+    fn rows_for(reads: Vec<AlignedRead>, reference: &[u8]) -> String {
+        let mut out: Vec<u8> = Vec::new();
+        write_feature_header(&mut out).unwrap();
+        stream_features_region(
+            SliceSource::new(reads),
+            Arc::from(reference.to_vec().into_boxed_slice()),
+            0,
+            0..reference.len() as u32,
+            PileupParams::default(),
+            &mut |col| write_feature_row(&mut out, "chr1", col).map_err(CoreError::from),
+        )
+        .unwrap();
+        String::from_utf8(out).unwrap()
+    }
+
+    #[test]
+    fn feature_row_has_the_expected_fields() {
+        // Reference AAAA; at pos 0, two reads: A (ref, fwd) and C (alt, rev).
+        let reference = b"AAAA";
+        let reads = vec![mread(0, b"AAAA", false), mread(0, b"CAAA", true)];
+        let text = rows_for(reads, reference);
+        let line0 = text.lines().nth(1).expect("a data row"); // after header
+                                                              // contig pos ref depth raw a c g t a_fwd a_rev c_fwd c_rev g_fwd g_rev t_fwd t_rev mean_bq mean_mapq
+                                                              // pos 1 (1-based), ref A, depth 2, raw 2, a=1 c=1, a_fwd=1 c_rev=1, bq30 mapq60.
+        let f: Vec<&str> = line0.split('\t').collect();
+        assert_eq!(f[0], "chr1");
+        assert_eq!(f[1], "1"); // 1-based
+        assert_eq!(f[2], "A");
+        assert_eq!(f[3], "2"); // depth
+        assert_eq!(f[5], "1"); // a count
+        assert_eq!(f[6], "1"); // c count
+        assert_eq!(f[9], "1"); // a_fwd
+        assert_eq!(f[12], "1"); // c_rev
+        assert_eq!(f[17], "30.00"); // mean_bq
+        assert_eq!(f[18], "60.00"); // mean_mapq
+        assert_eq!(f.len(), 19);
+    }
+
+    #[test]
+    fn feature_stream_is_deterministic_and_bounded() {
+        use crate::core::MemoryBudget;
+        // 50k single-base reads tiled at depth ~1: identical output twice, tiny WS.
+        let reference = vec![b'A'; 50_000];
+        let reads: Vec<AlignedRead> = (0..50_000u32).map(|p| mread(p, b"C", false)).collect();
+        let a = rows_for(reads.clone(), &reference);
+        let b = rows_for(reads, &reference);
+        assert_eq!(a, b, "feature output must be byte-identical run-to-run");
+
+        let reference2 = vec![b'A'; 50_000];
+        let reads2: Vec<AlignedRead> = (0..50_000u32).map(|p| mread(p, b"C", false)).collect();
+        let mut max_ws = 0u64;
+        let (ws, _skips) = stream_features_region(
+            SliceSource::new(reads2),
+            Arc::from(reference2.into_boxed_slice()),
+            0,
+            0..50_000,
+            PileupParams::default(),
+            &mut |_col| Ok(()),
+        )
+        .unwrap();
+        max_ws = max_ws.max(ws.bytes);
+        // Bounded by coverage, not the 50k read count (reference ~50k + a tiny active set).
+        assert!(
+            ws.fits(MemoryBudget::from_mb(1)),
+            "feature working set {max_ws} should fit 1 MiB"
+        );
+    }
+}

src/call/mod.rs

@@ -2,13 +2,17 @@
 //! abstention-aware variant calls. Built on `crate::core` + `crate::pileup`
 //! only; no VCF writing or CLI wiring (those are later phases).
 
+pub mod features;
 pub mod germline;
 pub mod pipeline;
 pub mod plan;
 pub mod somatic;
 pub mod types;
 pub mod whole_genome;
 
+pub use features::{
+    stream_features_region, stream_features_whole_genome, write_feature_header, write_feature_row,
+};
 pub use germline::call_germline;
 pub use pipeline::{
     call_germline_region, call_germline_region_streaming, call_germline_region_tracked,