fix(chrombpnet): invert log-counts with expm1, not exp - #113
Merged
Conversation
ChromBPNet's count head is trained against log(1 + count) — upstream batchgen_generator.py feeds np.log(1+batch_cts.sum(-1, keepdims=True)) as the count target — but chorus inverted it with np.exp, so every recovered count was high by exactly +1. That is ~0.1% at a 1,000-count peak and up to 100% at a low-activity site, which is precisely the regime the activity CDFs are built from. The bug was self-consistent: the oracle and the CDF builder made the same error, so ChromBPNet's percentiles stayed internally valid, which is why it went unnoticed. Raw counts and cross-oracle comparability were not valid. Cherimoya already does this correctly and is unaffected (cherimoya_source/scoring.py). Three count-inversion sites, all changed: - oracles/chrombpnet.py:579 _transform_predictions_to_tracks - oracles/chrombpnet.py:802 predict_sliding - scripts/build_backgrounds_chrombpnet.py:348 predict_profiles_batch The profile-softmax np.exp calls at :577, :801 and :347 are a different transform (softmax over profile logits, not a log1p inverse) and are deliberately untouched. tests/test_chrombpnet_counts.py covers all three: behavioural tests on _transform_predictions_to_tracks across the dynamic range including a 1-count site, plus a builder/oracle consistency check so the CDF builder and query path cannot drift apart again. Verified to fail 9/9 against the pre-fix code and pass 9/9 after. chrombpnet_pertrack.npz must be rebuilt against this transform; the CDFs shift and are being regenerated separately. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The KLF1 enhancer-profile notebook does not go through ChromBPNetOracle._transform_predictions_to_tracks — it calls model.predict_on_batch directly and folds the count head in itself, in two cells (s4-cbp, 5f10e19e). Both used np.exp(count_logits), so they carried the same +1-count inflation the oracle had. Source only. The committed cell OUTPUTS were computed with the old inversion and are still stale; they need GPU re-execution along with the rest of the ChromBPNet-affected examples. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
lucapinello
added a commit
that referenced
this pull request
Jul 31, 2026
…defects
Re-checked my own findings against README.md and examples/notebooks/README.md
after Luca asked whether the docs were actually wrong or whether I just
hadn't followed them. Four claims were unfair to the docs and are
demoted:
- 'README assumes mamba exists' is FALSE: README.md:19 lists Miniforge
as prerequisite #1 in a 'Before you start' block.
- The chorus setup HF-token gate is documented BY DESIGN at
README.md:312, with the rationale ('so the other oracles aren't built
for nothing'), and README.md:40 front-loads the token as required.
Not a P1. (Nit that stands: :312 says '5 oracles', now 7.)
- The MCP e2e failure was a documented trap: README.md:1004-1011
describes the two-mamba-root problem AND gives the fix
(MAMBA_ROOT_PREFIX). With it set, the test passes in 287 s. Not
following documented troubleshooting was my error. The only real nit
is in the test's hardcoded ~/.local/share/mamba fallback.
- The notebook kernelspec P1 drops to P2: the ipykernel registration
step IS documented at examples/notebooks/README.md:45, so this is
AUDIT_CHECKLIST §6 omitting a prerequisite, not a product defect.
- Concurrent per-oracle setup colliding on GPU 0 is demoted too: no doc
suggests parallelising setup; that pattern was invented here.
Also updates the report to the finished state: 9/9 oracles installed and
Healthy, fast suite 475 passed / 4 skipped / 0 errors, AlphaGenome
backend-equivalence and MCP e2e both passing, and the reason the NPZ
upload is deliberately held until #113 merges.
The genuine findings are unchanged: the three P0s (Enformer track-0
substitution, Borzoi prediction_interval, Borzoi direct mode), the
_get_denominator percentile inflation, merge_shards discarding a rebuild,
the TF-env-direct CPU-only gap, device='cuda:N' escaping the CVD mask,
MCP drift, and the 26 stale AUDIT_CHECKLIST gates.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
lucapinello
added a commit
that referenced
this pull request
Jul 31, 2026
…114) * audit(2026-07-31): 8xH100 fresh full install + AUDIT_CHECKLIST pass Fresh install from nothing on Linux/8xH100 (driver 595.71.05, CUDA 13.2): bootstrapped Miniforge (the box had no conda at all), built the base env plus every per-oracle env that needs no credential, downloaded weights and hg38, ran the fast and integration suites, and worked all 18 AUDIT_CHECKLIST sections. Every failure and every claimed stale gate was independently re-run before being recorded. Headline results: - 7/9 oracles installed and Healthy; the AlphaGenome pair is blocked on an HF token, and `chorus setup` hard-gates ALL nine on it. - Cherimoya verified on the real Triton/CUDA path (HAS_TRITON=True, torch 2.13.0+cu130, triton 3.7.1, sm_90) — not the macOS CPU fallback. Integration 9 passed / 1 skipped. - Fast suite 474 passed / 4 skipped / 1 error (the error is the AlphaGenome smoke test, i.e. the same token gap). - Determinism green: all 7 installed oracles bitwise identical. Three P0s found that are unrelated to any current work and survive green CI because the direct (use_environment=False) predict path is untested: Enformer silently returns track index 0 for any non-ENCFF track id (all 638 FANTOM CAGE tracks); Borzoi variant scores silently return None from a 2.67x-too-wide prediction_interval; Borzoi is entirely broken in direct mode. Also P1: effect percentiles inflated ~4.1% by _get_denominator, and Enformer/ChromBPNet run silently CPU-only whenever their env is used directly — the exact invocation CLAUDE.md prescribes. 26 AUDIT_CHECKLIST gates are themselves stale (6-oracle/22-tool counts, ChromBPNet "24 CDFs", ">=334 pass / <=1 skip", selenium method, and more) and manufacture failures if run verbatim. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * audit(2026-07-31): concrete Phase 6 inventory, split by what blocks it 33 files under examples/ reference ChromBPNet. Separates the artefacts that can be regenerated as soon as the rebuilt NPZ lands from the multioracle ones that cannot, because regenerate_multioracle.py iterates ('chrombpnet','legnet','alphagenome') and no *_variant_report.pkl is committed. Also records that there are four regeneration entry points, not two, and that the SORT1_chrombpnet output is already stale from the 2026-06-17 windowing fix so that drift is not misattributed to expm1. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * audit(2026-07-31): rebuild complete — 786 tracks, all gates pass, per-family impact Both CHIP shards finished 372/372; merge-shards reported '42 existing + 744 new', so the append-and-dedup trap did not fire. Rebuilt NPZ is sha256 cf30e17f..., 83.0 MB, with track_id set AND order identical to shipped, counts uniformly 18672/34004/1088128, every CDF row monotone, and zero all-zero rows. Records the per-family impact, which the ATAC/DNASE half alone hid: the 744 CHIP/BPNet rows carried a ~30% median count inflation (median true count 3.3) versus ~1.2% for ATAC/DNASE (median 82.2). Effect-CDF rankings are preserved (row correlation >= 0.973) while raw counts are corrected. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * audit(2026-07-31): correct findings that were my process errors, not defects Re-checked my own findings against README.md and examples/notebooks/README.md after Luca asked whether the docs were actually wrong or whether I just hadn't followed them. Four claims were unfair to the docs and are demoted: - 'README assumes mamba exists' is FALSE: README.md:19 lists Miniforge as prerequisite #1 in a 'Before you start' block. - The chorus setup HF-token gate is documented BY DESIGN at README.md:312, with the rationale ('so the other oracles aren't built for nothing'), and README.md:40 front-loads the token as required. Not a P1. (Nit that stands: :312 says '5 oracles', now 7.) - The MCP e2e failure was a documented trap: README.md:1004-1011 describes the two-mamba-root problem AND gives the fix (MAMBA_ROOT_PREFIX). With it set, the test passes in 287 s. Not following documented troubleshooting was my error. The only real nit is in the test's hardcoded ~/.local/share/mamba fallback. - The notebook kernelspec P1 drops to P2: the ipykernel registration step IS documented at examples/notebooks/README.md:45, so this is AUDIT_CHECKLIST §6 omitting a prerequisite, not a product defect. - Concurrent per-oracle setup colliding on GPU 0 is demoted too: no doc suggests parallelising setup; that pattern was invented here. Also updates the report to the finished state: 9/9 oracles installed and Healthy, fast suite 475 passed / 4 skipped / 0 errors, AlphaGenome backend-equivalence and MCP e2e both passing, and the reason the NPZ upload is deliberately held until #113 merges. The genuine findings are unchanged: the three P0s (Enformer track-0 substitution, Borzoi prediction_interval, Borzoi direct mode), the _get_denominator percentile inflation, merge_shards discarding a rebuild, the TF-env-direct CPU-only gap, device='cuda:N' escaping the CVD mask, MCP drift, and the 26 stale AUDIT_CHECKLIST gates. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
This was referenced Jul 31, 2026
lucapinello
added a commit
that referenced
this pull request
Jul 31, 2026
…bias shapes (#120) Three separate defects in the BPNet/CHIP path, all of which had to be settled together because each changes what the 744 CHIP CDF rows must be built from. The convention question is answered from the vendored TRAINING code, not by analogy: chorus/oracles/chrombpnet_source/templates/BPNet/ is upstream kundajelab/bpnet-refactor (losses.py byte-identical; arch.py and custommodel.py differ only by import shims). 1. JOINT strand softmax, not per-strand. custommodel.py:93-104 (the orig_multi_loss=False branch) reshapes the (B,1000,2) profile to (B,2000) and applies ONE multinomial_nll, and losses.py normalizes over the last axis. So the trained profile parameterization IS a single softmax spanning both strands. The oracle softmaxed each strand separately and scaled EACH by the full count total, making the two emitted tracks sum to 2.00x what the model predicted (measured exactly 2.00x). The builder did a third thing again - summed the strand LOGITS before one softmax - whose 501bp window sum drifted 0.84-1.30x versus per-strand across five loci, i.e. sequence-dependently, so the two could never be reconciled by rescaling. Which branch trained these weights is provable, not assumed: orig_multi_loss=False gives a (2,1) logcounts kernel, True gives (3,2) and REFUSES to load the real h5. All 744 cached BPNet models have (2,1). 2. The count inverse is exp(C) - n_tracks, not expm1(C). bpnet-refactor's generator builds the count target PER TRACK as log1p - np.log(np.sum(profile_predictions, axis=1) + 1) over an array shaped (batch, positions, tracks), so the sum is over positions - and custommodel.py:57 pools a task's tracks with reduce_logsumexp. The trained target is therefore log(n_tracks + total), so the inverse is exp(C) - n_tracks. For the 42 single-track ATAC/DNASE models n_tracks == 1 and this IS expm1(C): verified BIT-IDENTICAL (max diff exactly 0.000e+00) against the pre-change formula, so PR #113 stands and those 42 CDF rows stay valid. For the 744 two-track CHIP models expm1 left exactly ONE read of inflation - measured 1.12-1.95x on real loci, worse at quiet sites, median 1.78x on background 501bp window sums. Kept expm1 for n_tracks == 1 since it is more accurate for small C. (Upstream's own bpnet/cli/predict.py uses a bare exp(C), i.e. it is off by +n_tracks in the other direction. The target construction is the authority.) 3. The count-bias input shape was hardcoded (N, 1) for a (None, 2) input. BPNet logsumexp-reduces that bias before the final Dense (main_logsumexp_counts_bias_0), so the wrong width - which Keras silently broadcasts rather than rejecting - made the term log(1)=0 instead of log(2)=0.6931 and shifted every predicted log-count down by a constant w*log(2) = 0.588506 with w = 0.849035 (the Dense weight; matches the measured shift to all printed digits). Counts came out 1.80x too low at a peak and up to 3.04x at a quiet site. Profile logits were unaffected. The background builder derived the shape correctly and was NOT affected, so the oracle silently disagreed with its own CDFs. Fixed at all four call sites via _zero_bias_inputs, which derives the shapes from model.inputs: _predict_direct, predict_sliding, the env-path predict_template.py, and the public one-shot helper in oracles/bpnet.py (whose docstring also prescribed the superseded per-strand exp(counts) recipe). Also fixed: predict_sliding used SUM-LOGITS for CHIP while _predict used per-strand, and reported the result under a ':+' track_id. It now takes the joint softmax and emits the PLUS strand, matching the id it reports and the pooled per-strand CDF row it is looked up against. Emitting a both-strand total there would put 1.53-2.93x the _predict value behind the same id. Verified end-to-end on real weights: builder and oracle now agree to 2.0e-07 on every strand/locus combination, the two strands conserve exp(C)-2 exactly, and DNASE:K562 is bit-identical to before. CONSEQUENCE: the 744 CHIP CDF rows must be rebuilt against this transform; the 42 ATAC/DNASE rows must NOT (they are unchanged, and merge_shards appends and de-dups by track_id, so the base must contain only those 42 at merge time). Not yet rebuilt - see the PR. tests/test_chrombpnet_counts.py grew to 23 tests, including a helper that constructs the count target the way training does (per-track log1p pooled with logsumexp) so the expectations cannot drift from the loss, an explicit bit-identity test for the single-track inverse, and a check that expm1 overstates a two-track model by exactly one read. Fast suite: 491 passed, 4 skipped, 0 errors. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
lucapinello
added a commit
that referenced
this pull request
Jul 31, 2026
…#117) The ChromBPNet per-track CDFs changed (#113 + the HF background rebuild), so every committed ChromBPNet artefact drifted. Regenerating them also surfaced three artefacts stale against EARLIER merged changes that were never re-executed. All four effects are separated so nothing is misattributed to expm1. 1. expm1 (the intended change). ChromBPNet raw counts come down by exactly the window softmax mass w; percentile RANKINGS unchanged: SORT1 DNASE:HepG2 ref 287.857 -> 287.237 (-0.620) alt 747.653 -> 747.014 (-0.639) log2FC +1.37394 -> +1.37581 effect %ile 0.9995, activity %ile 0.9059 (unchanged) klf1 locus profile max 127.333 -> 127.307, mean 0.462 -> 0.461 2. SORT1_chrombpnet was stale since 2026-05-09, predating the 2026-06-17 windowing fix (#94). That is the large jump (ref 48.53 -> 287.24, log2FC +0.318 -> +1.376), NOT expm1 — the post-#94 multioracle artefact already carried +1.37394, which the regenerated value matches to within the expm1 residual. 3. klf1's EPInformer-seq cells were stale against the 2026-06-04 "widewin-roadmap retrain (sparse H3K27ac)" weights swap on lucapinello/chorus-epinformerseq-v2 (HF commit 2654c05b7c) — committed alongside it but never re-executed. Hence K562 H3K27ac mean 1.771 -> 0.132. Verified NOT device-dependent: cpu vs cuda agree to 0.006 on a 3.78 range (0.16%). epinformerseq.py:416 calls hf_hub_download with no `revision`, which is what let this drift silently. 4. The AlphaGenome report was stale against the CAGE minus-strand dedup (#99/#101); it still carried both strands. all_scores 61 -> 32 is entirely CAGE:HepG2 58 -> 29 — the duplicate #101 exists to remove. No track lost: same 4 unique descriptions before and after. Notebooks re-executed on H100, zero error outputs: klf1 14/14 cells, advanced_multi_oracle 57/57, comprehensive_showcase 38/38. Registering the documented `chorus` kernelspec (examples/notebooks/README.md:45) was required first. Re-execution also dropped a foreign-home path leak (/Users/jieconglin) from klf1's outputs. DELIBERATELY NOT REGENERATED — blocked, needs a decision: rs12740374_SORT1_legnet_report.html and rs12740374_SORT1_multioracle_report.html are left at their origin/main bytes. #99 (2026-06-18) added "tile LegNet across the locus for the IGV browser", but the legnet artefact was last committed 2026-06-17 (#95, single-window), so the tiling has never been applied. Applying it embeds a locus-wide 1-bp array and the report becomes 137 MB / the consolidated one 145 MB — both above GitHub's 100 MiB file limit, i.e. unpushable. Dropping legnet's .pkl instead yields a shippable 9.55 MB but silently loses legnet's 3 IGV tracks (21 -> 18). The committed 9.93 MB / 21-track artefact came from a pre-#99 single-window .pkl and is not reproducible with current code. So #99 has been unshippable since it merged. Options: bin/downsample the embedded IGV arrays, revert the tiling, or move these HTMLs to Git LFS. Consequence: the consolidated example_output.{md,json} here carry the corrected ChromBPNet numbers while the consolidated HTML does not, until that is resolved. Fast suite: 477 passed, 4 skipped, 0 errors. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
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.
What
ChromBPNet's count head is trained against
log(1 + count)— upstreamchrombpnet/training/data_generators/batchgen_generator.pyfeedsnp.log(1+batch_cts.sum(-1, keepdims=True))as the count target — but chorus inverted it withnp.exp. Every recovered count was therefore high by exactly+1: negligible at a peak (~0.1 % at 1,000 counts), up to 100 % at a low-activity site, which is precisely the regime the activity CDFs are built from.The bug was self-consistent — the oracle and the CDF builder made the same error — so ChromBPNet's percentiles stayed internally valid, which is why it went unnoticed. Raw counts and cross-oracle comparability were not.
Cherimoya already did this correctly (
cherimoya_source/scoring.py) and is untouched; its docstring is what flagged this.Three count-inversion sites
chorus/oracles/chrombpnet.py_transform_predictions_to_trackschorus/oracles/chrombpnet.pypredict_slidingscripts/build_backgrounds_chrombpnet.pypredict_profiles_batchThe profile-softmax
np.expcalls at:577,:801and:347are a different transform (softmax over profile logits, not alog1pinverse) and are deliberately unchanged. The diff is exactly three functional lines.Tests
New
tests/test_chrombpnet_counts.py(9 tests), mirroringtest_cherimoya.py::test_expected_counts_uses_expm1_not_exp:_transform_predictions_to_trackswith a stubself, so no TF/weights/env needed — alog(101)count head must yield 100, not 101expm1(log1p(c)) == cround-trips across0 … 10,000predict_sliding, which needs a live TF model to reach behaviourallycherimoya_source/scoring.pywarns aboutVerified the suite fails 9/9 against the pre-fix code and passes 9/9 after (bug reintroduced by hand, then reverted).
Fast suite on 8×H100 / Linux: 474 passed, 4 skipped, 1 error. The single error is
TestSmokeAlphagenome::test_predict, which needs thechorus-alphagenomeenv — absent on this box becausechorus setuphard-gates on an HF token that is not yet available here. Unrelated to this change.Docs corrected in the same pass
These asserted the buggy behaviour and become false with this commit:
chorus/oracles/cherimoya_source/scoring.py— said chorus "recovers its counts withnp.expin three places" (and cited:800; the count line is:802)docs/NORMALIZATION_GUIDE.md— same claimscripts/build_backgrounds_chrombpnet.py:327—predict_profiles_batchdocstring still saidsoftmax × exp(counts)tests/test_cherimoya.py:226— docstring, now past tenseREADME.md:1228— drive-by factual fix: the ChromBPNet CDF sample sizes were listed as 10,000/31,500; the shipped NPZ carries 18,672/34,004 (confirmed by reading it, and byNORMALIZATION_GUIDE.md:212-218)CHANGELOG.md—[Unreleased] / FixedentryFollow-up (not in this PR)
chrombpnet_pertrack.npzmust be rebuilt against the corrected transform; all three CDFs shift. In progress separately — do not merge-and-forget, or the shipped CDFs andpredict()will disagree.CHIP:{cell}:{TF}:+/:-tracks that each carry the full count mass. So 744/786 CDF rows measure a different quantity than any trackpredict()returns — a separate correctness issue from this one, worth its own decision.🤖 Generated with Claude Code