[jjschirle]

Stacked on #71 — based on subset-unpack, so the diff here is just the two
commits on top. Implements the one item #71 lists under "Not in this PR":

Bounded memory (step 4 of the plan). Peak memory is still the whole packed
stream plus the output... Reading the stream in slabs of zstd chunks with the
scan state carried across slab boundaries needs a slab-wise version of
_compression.read_array. It's a separate change and matters at Parse-10M scale.

and answers the ask in #39 ("I think I'd want to see benchmarking to understand
the performance implications") with numbers rather than an argument. Those
numbers are mixed, and the summary up front is: ~2x the wall time for ~23x
less memory, plus files a full decode cannot read at all. Whether that trade
is worth carrying is a judgement call, and the benchmark is here so it can be
made on evidence.

What it adds

• _compression.read_range(ds, start, stop, *, table=) — decodes only the zstd
  chunks covering a byte range, through the same chunk_table/pread/
  numcodecs path Replace Blosc2 with zstd-1 and decode chunks in parallel on read #68's full read already uses. read_array is now that over
  [0, len). chunk_table is module-public so a caller reading many windows
  builds it once; it walks every chunk through a Python callback and HEK293T's
  packed_indices has 190,770 of them.
• _ivcsc.trim_to_boundary — a window out of the middle of the stream usually
  ends mid-varint, and every helper below assumes its last varint is complete.
  _count_varints_range in particular walks its inner continuation-bit loop
  past end, which on a truncated final varint reads out of bounds, and njit
  without boundscheck will not stop it. _varint_chunks's docstring now says
  the buffer must end on a boundary.
• vsparse.iter_packed(path, chunk_size=), over iter_blocks/open_packed.

On "streaming is not easily compatible with compression": that is true of
random-access streaming and not of this. read_range fetches whole zstd
chunks and never seeks into a compressed one, so a windowed read is exactly as
compatible with compression as reading all of it. The stream is still walked
front to back — value_ptr counts elements, not bytes, so nothing in the file
says where slice r begins. The only compression-specific cost is re-decoding
at most one boundary chunk per window.

Why a block iterator rather than another selective read

obs_idx= cannot express a reduction over every cell. Per-perturbation means
across 4.5M cells are one pass here; as selective reads they are one pass per
group, each running the whole stream again, because nothing indexes a slice's
bytes. Same for anything that has to see all of X once — which is what
backfill_stats-style work and the downstream bi-cross-validation reduction
both need.

A contiguous block needs none of the selective machinery: its pointer arrays
are slices of the parent's, so only the indices are decoded, by the existing
unpack_indices. Locating where the next slice starts reuses the window's own
span table from _varint_chunks and walks at most one span, rather than paying
a second counting pass over the window.

Three of the four _ivcsc helpers an earlier draft of this added turned out to
be redundant against #71 and were dropped: _varint_chunks is the same function,
and _locate_ranks/_slice_starts subsume an element_offsets helper with a
better block-skipping walk. The diff is ~1/5 of what it was as a result.

The benchmark caught a bug in this PR's own code

The first run reported the slab reader 3.1x slower than a full decode. The
cause was here, not in the method: windows were sized assuming ~1 byte per
varint. That holds on HEK293T (1.138 bytes/nonzero) and not on the benchmark
fixture (1.664), so every read came up ~40% short and retried — re-reading and
re-counting, which is the expensive half. Sizing from the file's own
len(packed_indices) / value_ptr[-1] (free) and rescaling a retry by what the
window measured fixes it. slab_time_ratio_over_full is gated where that
regression gets caught.

Measurements

New @slow case slab_stream_vs_full_decode, beside #71's
subset_load_vs_full_load. It writes its own temp .h5ad (400,000 x 3,000,
2% dense). Memory is measured with tracemalloc, not RSS: ru_maxrss is a
high-water mark that never falls, so on a fixture small enough to build
in-process it only ever reports the fixture. Timing and memory are separate
passes, since tracing distorts allocation-heavy timing.

full_decode_seconds=0.0482   slab_stream_seconds=0.0386
full_decode_gb_per_second=0.829   slab_stream_gb_per_second=1.035
slab_time_ratio_over_full=0.801   slab_peak_ratio_over_full=0.669
full_decode_peak_gb=0.202   slab_stream_peak_gb=0.135
bytes_per_nonzero_read=1.664   rows_per_window=279,620

Gated: slab_time_ratio_over_full (margin 1.5) and slab_peak_ratio_over_full
(1.1). The other metrics are recorded for context. --record also fills in
entries for the two pre-existing @slow cases, which had none; I removed those
again, since starting to gate them is not this change's business.

Whole-file pass on the two real corpora, 48 threads, warm cache:

file	path	wall	GB/s	peak RSS
Parse-10M (21.88 GB packed, 18.8 B nnz)	full decode	11.7 s	1.88	63.14 GB
Parse-10M	slabs, 1123 windows	22.1 s	0.99	2.75 GB
HEK293T Xaira (33.17 GB packed, 29.1 B nnz)	slabs, 1737 windows	40.2 s	0.82	4.66 GB
HEK293T Xaira	full decode	—	—	116 GB of indices alone; does not fit in 137 GB

So a run-through really does exceed 1 GB/s and a windowed one is just under it.
The gap is per-window overhead, which _WINDOW_ELEMS (currently 2^24 elements,
~18 MB of bytes and ~67 MB decoded) trades against residency; the benchmark is
the place to retune it.

What the ~23x actually refers to

It is not 23x less of the same thing. A full decode holds the packed stream and
the fully decoded index array at the same time; the slab pass holds neither
in full. On Parse-10M:

held at peak	full decode	slab pass
packed_indices (the byte stream)	21.88 GB	one window, ~0.02 GB
decoded indices (18.83 B x 2 B, uint16)	37.66 GB	one block, ~0.03 GB
value_ptr	1.69 GB	1.69 GB
values	0.42 GB	0.42 GB
major_ptr	0.08 GB	0.08 GB
predicted total	61.73 GB	~2.24 GB
measured ru_maxrss	60.92 GB	2.75 GB

So the honest statement is "2.2 GB of pointer arrays plus a working set of a
few tens of MB, against 62 GB". About 80% of the slab pass's 2.75 GB is floor
that both paths pay; what is removed is the 21.88 GB stream and the 37.66 GB
index array, neither of which is ever materialized whole.

(Parse-10M stores uint16 indices — 40,352 genes fits — so its decoded index
array is 37.66 GB rather than the 75 GB an int32 file of the same nnz would
need. HEK293T is int32, hence its 116 GB.)

Neither ratio is a constant, and quoting them without the file is misleading:

• The memory ratio grows with file size. Its numerator scales with nnz; its
  denominator is a fixed floor plus one window. On this PR's own benchmark
  fixture the ratio is only 1.5x, because that fixture is ~1.4 windows wide —
  the window is the file.
• The time ratio is per-window overhead, and is tunable. 1123 windows here; on
  the ~1.4-window fixture the slab pass is faster (0.80x). _WINDOW_ELEMS
  trades the 2x back against residency and has not been tuned.
• On HEK293T there is no time comparison, because the full decode does not
  run: 116.5 GB of indices on a 137 GB machine. There the trade is not slower
  versus faster.

major_ptr, value_ptr and values are still read whole — ~12 bytes per value
group against 4 per nonzero, 4.1 GB rather than 116 GB on HEK293T. That is this
module's floor and it scales with value groups, not nnz.

Correctness on real data

A streamed pass over Parse-10M accumulating per-gene column sums reproduces that
file's X/stats/col_sums bitwise (np.array_equal true, max abs diff 0.0),
in 28.1 s at 2.80 GB. That sidecar was written by the ordinary eager
compute_stats path during write_h5ad, so it is an independent reference
rather than a restatement of this code.

Not in this PR

• load_packed(obs_idx=) itself is still not slab-bounded. This adds the
  cursor and proves it out on whole-file passes; threading it through
  unpack_selected needs that function to fill a preallocated output at an
  offset, so a selection spanning windows doesn't concatenate. Natural follow-up,
  same cursor.
• Retuning _WINDOW_ELEMS against the new benchmark. 2^24 was chosen to keep a
  decoded block near 67 MB, not because it was measured to be optimal.
• Background prefetch of the next window — worth roughly 2x on an I/O-bound
  pass, at the cost of a thread and an error-propagation path.

Notes

VCSCAnnData._subset_fields now delegates to a module-level subset_fields
taking a field mapping, with no behaviour change, because iter_packed windows
fields read straight from a file rather than off an existing object.

Blocks deliberately carry no statistics sidecar: row_sums/row_nnz would
window correctly, but col_sums/col_nnz/col_depth_sums are whole-matrix
quantities and a block's are not the file's. Relatedly, normalizing block by
block is exact for a recipe whose scales are per-cell (cp10k_log1p is
log1p(x * 1e4 / row_sums[i])) and not for parafac2/pearson, whose
median(row_sums) target and col_depth_sums gene scale are whole-population;
the docstring says which, and where to get them instead.

iter_packed reads obs once, whole — ~3 GB of Python strings on HEK293T,
whose _index/sample/guide_target are nullable-string-array over 4.5M
rows. Paid once per call, not per block. iter_blocks skips it entirely, which
is why the figures above use it.

Checks

• pytest: 1547 passed, 72 skipped. tests/test_iter_packed.py adds 20,
  covering every chunk size against the same in-memory slice, empty rows at both
  ends, an all-zero matrix, a zero-row file, IVCSC raising, and — the one that
  matters — that the cursor's byte total lands exactly on
  len(packed_indices), which is what catches an off-by-one in locating the
  next slice. tests/test_parallel_read.py adds 10 for read_range across
  chunk boundaries, partial first/last chunks, empty ranges, a cached table and
  every non-zstd fallback.
• ruff, codespell, formatting: clean. radon cc -n C: nothing ranked C or
  worse. vulture: nothing new. xenon: only the two pre-existing module-rank-B
  files. jscpd: no new clones.
• --set fast: no regressions. --set slow: no regressions past recorded
  ceilings.

🤖 Generated with Claude Code