feat(concat): add gvl.concat for merging datasets along regions or samples - #346
Merged
Merged
Conversation
Design for gvl.concat over the (region, sample) grid, plus the finding that mcvickerlab#334's first ask needs no feature: gvl.write already supports `samples=` against the parent PGEN, so pre-splitting with plink2 is what forces the per-cohort index rebuild. Measured on NFSv3 rather than assumed: memmap reads are fastest (316 MB/s) but memmap writes are ~21x slower than buffered writes (15.2 vs 316 MB/s), so the merge streams memmap-read -> buffered-write in 16 MiB chunks, destination-ordered, single-threaded. Keeps the variants.arrow hardlink and adds a fingerprint instead of copying: gvl never mmaps that file (memory_map=False) and genoray rewrites its indexes via atomic_write_path, so the SIGBUS concern does not apply. Mirrors the existing svar_link/svar2_link precedent. Relates to mcvickerlab#334. Files mcvickerlab#337 (fingerprint guard for write()'s hardlink) and mcvickerlab#338 (buffered writes for bulk output on NFS). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…low on NFS too The first benchmark did not evict the page cache and its sources had just been written, so its "memmap read -> buffered write = 316 MB/s" row was measuring RAM, not NFS. Re-measured with posix_fadvise(DONTNEED) on both NFSv3 and node-local XFS, 3 reps each. Corrected finding: on NFSv3 *any* memmap in the path costs ~5-6x (buffered/buffered 93 MB/s vs ~14-17 MB/s for every combination involving a memmap), because faults go out as 4 KiB RPCs instead of using the mount's 1 MiB rsize/wsize. On local XFS all four patterns are within noise of ~100 MB/s, so memmap is not harmful there and buffered IO is a safe unconditional choice. The design accordingly uses buffered seek()+read() for sources as well as buffered write() for destinations — no memmap in the bulk data path. Runs are already coalesced into large contiguous ranges, so this loses nothing relative to a memmap gather. Cost model updated from ~1h to ~3h per merged TB. Adds an explicit scope limit: these are sequential-streaming numbers and say nothing about the random fancy-indexed access in gvl's read path. Relates to mcvickerlab#334. Updates mcvickerlab#338. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Eight tasks over four new modules (_concat_plan, _concat_io, _concat_validate, _concat). Tasks 2-4 are independent and dispatch in parallel; 5-8 are sequential. Carries the spec's measured constraints into every task: buffered IO on both sides in 16 MiB chunks with no np.memmap in the bulk path, single-threaded, destination-ordered, and neither pass materializing a full array. Test oracle is asymmetric by design — byte-identity on the region axis, read-equality only on the sample axis, since extend_to_length sizes each region's window to the cohort present at write time. Relates to mcvickerlab#334. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Pure planning layer for gvl.concat. Both axes reduce to one provenance map over the (R, S[, P]) C-order grid; coalesce() compresses it into maximal contiguous runs so the IO layer moves large byte ranges. Relates to mcvickerlab#334.
Region- and sample-sharded datasets plus a single-shot whole dataset, which later tasks compare against. Relates to mcvickerlab#334.
All checks run before any bytes move. Variant-source identity reuses the bounded blake2b fingerprint idiom from _fasta_cache rather than hashing a multi-GB Arrow file. Relates to mcvickerlab#334.
Buffered read/write on both sides in 16 MiB chunks, no np.memmap in the bulk path: on NFSv3 any memmap in the path costs ~5-6x, while on local XFS the patterns are equivalent. Relates to mcvickerlab#334.
test_link_or_copy_produces_identical_bytes now asserts the destination is a genuine hardlink (same st_dev/st_ino), not just byte-identical, so a regression that silently fell through to shutil.copyfile would be caught. Adds coverage for the EXDEV fallback path (proves a real copy landed, distinct inode) and for non-EXDEV OSErrors propagating instead of being swallowed into a copy. Relates to mcvickerlab#334.
…orce one variant source has_dosages checked <dataset>/genotypes/dosages.npy, a path gvl.write never writes -- dosages only ever live inside the externally-linked svar/svar2 store. Resolve it the same way _haps.py does (via the svar/svar2 link), falling back to False for pgen/vcf and tracks-only datasets. validate_concat also only compared the coarse backend string, so two unrelated pgen_vcf datasets with different variant tables would pass. Add a ConcatInput.fingerprint (defaulting to None), populated from the linked store's existing fingerprint for svar/svar2 or via variants_fingerprint for pgen_vcf, and raise on a mismatch when both sides have one -- enforcing the "one variant source" global constraint. Relates to mcvickerlab#334.
Merges genotypes, regions, and metadata along either axis. variants.arrow is hardlinked from the first input and guarded by a recorded fingerprint. Relates to mcvickerlab#334.
…nsistency Two integration-coverage gaps from review: nothing byte-compared merged regions.npy for axis="regions" (the Correction-3 gather path), and the sample-axis genotype merge was only unit-tested for the S_d failure mode. Adds a regions.npy byte-identity check against a single-shot write for axis="regions", and a self-consistency check (merged variant_idxs slice must equal the contributing shard's own slice, per (region, sample, ploid)) for axis="samples" -- self-consistency rather than byte-identity against `whole`, since extend_to_length legitimately makes a shard's stored variant set differ from the full write's. Also adds interleaved_sample_shards (shard0=[s0,s2], shard1=[s1]) so the sample-axis self-consistency check runs against a merge that genuinely interleaves datasets after sorting, not one that reduces to block concatenation -- exercising _sample_order/_merged_bed on the case old block-concatenation code would get wrong. It passes. Relates to mcvickerlab#334.
Completes the per-store merge table: per-sample tracks gather over (R, S), annot tracks are sample-independent, and the svar/svar2 range arrays copy verbatim as fixed-stride gathers since they hold absolute ranges into an external store. Relates to mcvickerlab#334.
…axis Finding 1: axis="samples" linked annot_intervals from input #0 without the fingerprint-compare the design spec requires. Since each input's annot data is computed against its own regions.npy, and extend_to_length can diverge chromEnd per input (the same case Correction 3's elementwise-max exists for), copying input #0 verbatim could silently under-cover the merged dataset's annotation reads. Added _assert_annot_track_matches, a bounded fingerprint-compare (reusing the existing _fasta_cache.fingerprint idiom) that raises naming the offending input and track instead of copying divergent data. Finding 2: adds test coverage for annot tracks (previously untested on both axes) and region-axis per-sample tracks (previously untested; only the sample axis had coverage). Relates to mcvickerlab#334.
Converts an out-of-band variant-index rewrite from silent wrong data into a loud error. Absent fingerprints skip the check so pre-existing datasets still open. Relates to mcvickerlab#334, mcvickerlab#337.
Adds the api.md autodoc entry required to keep docs in sync with __all__, the concat workflow and cost model, and the parent.pgen + samples= guidance that answers the first half of mcvickerlab#334. Closes mcvickerlab#334.
Final whole-branch review fixes for gvl.concat, applied in one pass: - Reject within-input duplicate (chrom, chromStart, chromEnd) rows in validate_concat: the merged union sort is unstable and doesn't include strand, so a within-input tie can silently attribute the wrong stored region to a merged r_idx_map slot. Cross-input coordinate collisions remain legal. - Replace every `f.write(<array>.tobytes())` with `<array>.tofile(f)` in _concat.py, and write _gather_svar_offsets's two offset planes sequentially instead of np.stack-ing them, removing several full-size copies of exactly the arrays the design spec flags as the OOM risk at biobank scale. - Mark format.md's changelog row for the still-unreleased variants_fingerprint format bump as "(unreleased)" instead of claiming 0.41.0 as already shipped. - Raise a clear ValueError naming the offending input when inputs' bed schemas disagree (e.g. a missing `name` column), instead of surfacing polars' raw SchemaError from pl.concat. - Raise in Dataset.open when metadata records a variants_fingerprint but the hardlinked genotypes/variants.arrow is missing, instead of silently skipping verification. - Reuse input #0's already-stored r_idx_map on axis="samples" instead of recomputing it via a second unstable sort pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The Fix-1 duplicate-(chrom,chromStart,chromEnd) guard in validate_concat was applied to both axes, but axis="samples" has no tie-order exposure: _merged_bed reuses input 0's stored r_idx_map verbatim and regions.npy is never re-sorted on that axis. Applying the guard there blocked a realistic input shape (same window on both strands) that gvl.write accepts and that is the standard per-shard cohort-build-by-sample workflow. Narrow the guard to axis="regions", where it remains load-bearing against the unstable sp.bed.sort tie-order seam.
…at-spec # Conflicts: # docs/source/write.md # skills/genvarloader/SKILL.md
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #334.
Adds
gvl.concat(path, datasets, axis, *, overwrite=False, max_mem="4g"), merging on-disk GVL datasets along the region or sample axis without re-extracting genotypes.Implements
docs/superpowers/specs/2026-07-31-dataset-concat-design.md.Approach
A GVL dataset is a set of parallel ragged stores over an
(R, S[, P])C-order grid. The merge builds a provenance map (merged flat slot → source dataset + source flat slot), coalesces it into maximal contiguous runs, and streams each run with buffered IO. Both axes use the same gather;axisonly selects the provenance map._concat_plan.py_concat_io.py_concat_validate.py_concat.pyconcat()entry point, orchestration._open.pyvariants.arrowfingerprint at open.Backends: PGEN/VCF,
.svar,.svar2. Per-sample tracks and annotation tracks are merged too.Design constraints held throughout: buffered IO on both sides in 16 MiB chunks with no
np.memmapin the bulk path (a memmap there costs ~5-6x on NFSv3), single-threaded, destination-ordered iteration, and no full-array materialization — offsets alone reach 16-32 GB at biobank scale.Notable corrections found during implementation
Several defects in the original plan were caught by review and fixed. Recording them since they are the substance of the work:
_write.py:267sorts samples unconditionally, and all stores are written in sorted region order. These coincide only when shards don't interleave, which is true of the obvious fixtures and false in general (e.g. sample shards[s0,s2]+[s1]).provenance()gained anorder=parameter; all six call sites pass it.genotypes/dosages.npy, a pathgvl.writenever creates. Real dosages live in the externally linked.svarstore and are preserved by the link — no merging needed. Block removed.regions.npyis an elementwise max, soextend_to_lengthdivergence across shards silently under-covered the merged windows. The spec's required fingerprint compare was missing; added.sp.bed.sortis an unstable polars sort that excludesstrand, so an input containing duplicate(chrom, chromStart, chromEnd)rows could produce anr_idx_mapdisagreeing with the payload — a silent wrong read. Rejected in validation onaxis="regions"; deliberately not applied toaxis="samples", which has no such exposure and where the same window on+/-strands is a legitimate, common input.f.write(arr.tobytes())allocated a third full-size copy of exactly the array the spec names as the OOM risk; nowarr.tofile(f).Testing
Full tree green: 1178 passed, 56 skipped, 4 xfailed.
Tests are built to be discriminating rather than merely present — the region and sample fixtures interleave across datasets, so they genuinely fail under the pre-correction ordering. Sample-axis genotypes are checked for self-consistency against the contributing shard (not byte-identity against a single-shot write, which
extend_to_lengthmakes an invalid oracle); region-axis output is asserted byte-identical.Follow-ups filed
.svar2concat has no test coverage (verified consistent with the on-disk format by inspection, but never executed; needs ansvar2_pathfixture). Carries three deferred minors.ulimit -nwhen merging many hundreds of shards. Deferred rather than reworking the tested IO hot path late.Docs
api.md(sync gate passes),write.md(merge workflow + the parent-PGENsamples=guidance answering #334's first question),faq.md,format.md, andskills/genvarloader/SKILL.md.🤖 Generated with Claude Code