Fix the background null model: the sampler was discarding the tail it measures

[lucapinello]

Nine commits. The headline: the null's ceiling was wrong by construction, not merely noisy — and fixing it needed no GPU.

The defect

percentile = min(rank/denominator, 1.0) clamps once an effect passes the largest sampled background value. That had been patched three times (re-anchoring, union-at-2N, a read-side exceedance ratio). None was the cause.

merge_effect_shards.py called ReservoirSampler.from_flat_samples(*parts) with no capacity, silently inheriting DEFAULT_CAPACITY = 50_000. Every AlphaGenome gene_expression track offers 148,367 effect values, so the union was subsampled 2.97×. A uniform m-of-N subsample keeps the population maximum with probability exactly m/N: 50,000/148,367 = 0.3370, and 33.9% of the 667 RNA rows were measured to have kept theirs.

AlphaGenome effect null	tracks	max ratio (true/shipped)	p99 ratio
gene_expression	667	median 1.332, p90 3.18, worst 8.34	1.006
every other layer	4,501	1.0000	1.0000

The tail was wrong by up to 8.3× while p99 was right to 0.02%. That asymmetry is why this survived every calibration gate for months: reservoir sampling is unbiased for the body, and every gate measures the body.

Neither existing guard could have caught it:

• No shard exceeded its capacity (~18.5k offered against 20,000), so the builders' capacity was never the issue. The merge is a second, independent site.
• cdf_grid_violations is fed the offered count while checking geometry set by the retained count, then does if n >= n_points: continue. Offered is always ≥ 10,000, so it skipped every thinned row by construction.

Fixed by exact retention. No GPU — the 8 shards hold every offered value, so the union recovers all 148,367 in 23 seconds. capacity is now keyword-only with no default, making this specific defect a TypeError.

Verified against the raw samples, not against itself

• all 4,501 never-thinned rows come back bit-identical
• 0 of 5,168 tracks thinned; shipped maxima now equal the true population maxima
• release gates hold: CAGE 0.6250, RNA 0.7770, 0% saturated
• Borzoi/Enformer verified bit-exact and deliberately not rebuilt (their unions are under capacity, so the merge never takes the rng.choice path) — pinned as a test so the omission is a recorded negative result rather than an oversight

One gate I had specified wrongly: "p99 within ±10% for every track". 3 of 667 rows exceed it, 2 above / 1 below — unbiased noise in the old 50,000-sample estimate, not movement in the new one, which is now exact. Body invariance is verified at p50/p90, where 0.0% of rows fall outside 10%. At p99.9 the ratio shifts systematically up (median 1.022): that is the recovered tail mass, i.e. the point.

Two more defects found on the way

Sei was entirely dark. sei.py set assay_type = info.name for its sequence classes; classify_track_layer dispatches on the literal "sequence-class", so it returned "other", whose LAYER_CONFIGS entry is None, and every one of Sei's 40 tracks scored raw_score=None. It had a built, verified, zero-degeneracy null that no query could reach — and its absence from every committed example read as "we didn't include Sei" rather than "Sei can't be scored". One line. All 40 now score.

Anchored positions were being clamped onto the contig margin. usable selects whole chromosomes long enough for the margin, but the tss/junction/gene-body populations were filtered only by chrom in usable, never by the margin interval. 2,515 of 20,083 protein-coding TSS (12.5%) sit within 5 Mb of a contig end and were _clamped onto the boundary coordinate — up to 5 Mb from the TSS they were labelled as being within 1 kb of.

	before	after
positions exactly on a boundary	12.1–14.6% per stratum	0.0–0.1%
distinct positions (of 6,000)	5,265	6,000
worst duplicate coordinate	64× (chr16:5,000,000)	1×

The mislabelling matters; the duplication matters more. Identical positions give identical effect values, so they pad the sample count without adding information and manufacture tied CDF runs — the degeneracy _rank_with_tie_breaking exists to compensate for, injected by the sampler rather than present in the biology. This affected the shipped backgrounds.

Also in here

• DHS stratum, additive at scaled N (12,000 → 18,000), every pre-existing stratum keeping its exact absolute count. Justified by tail width and max(union) = max(max_a, max_b), not by calibration — and it does not fix motif-creation saturation: the DHS-anchored null ChromBPNet has always used still pins on rs12740374 CEBPA at 1.11×.
• A landmine disarmed. The stratum dispatch ended in a bare else drawing a uniformly random position, which doubled as the random handler, the empty-pool fallback and the catch-all for unrecognised names. Adding "dhs" without a branch would have shipped 6,000 uniformly random positions tagged, tallied and stamped as DHS, with entirely plausible-looking numbers. Both samplers now raise; the decisive guard is an annotation round-trip (≥99% of dhs within 150 bp of a Meuleman summit, ccre inside a cCRE, tss_near within 1 kb) — which is what found the clamping bug.
• effect_exceedance: above the ceiling the percentile is clamped and says nothing, so report the ratio to it. 1.011× / 1.109× / 2.973× / 10.109× now separate four effects that all read 1.000000. Deliberately not an extrapolated percentile — measured, a GPD overshoots the far tail by 3.8× and an exponential undershoots by 0.27×, while the plain empirical max is within 13%.
• HF_HOME no longer relocated — HF_HUB_CACHE is. HF_HOME is the parent of both the blob store and the huggingface-cli login token, so setting it orphaned the credential and every gated model failed telling users to run a command they had already run. Self-inflicted in the data-dir change; found by a measurement that needed AlphaGenome, not by a test.
• data dir defaults to the install tree rather than $HOME; LegNet and ChromBPNet list_cell_types fixes (each had hidden cell types that shipped background rows for).

Tests

1,181 fast tests pass. Every new guard was checked to fail on the pre-fix code.

Two tests corrected their own wrong premises while being written: unsigned CDF rows start at the smallest sampled magnitude (8.4e-6), not a structural zero; and the two support ends are asymmetric because searchsorted(side="right") counts values ≤ x, so exactly the maximum pins while exactly the minimum does not.

One was rewritten rather than retuned: a TSS-proximity test asserted an aggregate median, which additive union does not preserve — the TSS counts were untouched but the mixture median rose because DHS is legitimately distal. It now asserts per stratum (every tss_near position within 1 kb: measured 100%, median 413 bp) plus aggregate enrichment (13.6% vs 1.4% uniform).

Not in here

The unified GPU rebuild (effect + summary + perbin for all 8 oracles under one build_id) is planned but not started; it is gated on a sampler rewrite so thinning cannot recur.

Still open, and the same defect in other layers: perbin is thinned 16–43× on five oracles, and summary up to 5.2× on three.

🤖 Generated with Claude Code