Repository navigation
Integrate opt-in M4 position-profile exploration - #48
Merged
Merged
Conversation
A derived matrix is fitted as costs, and the Gram transform from similarities is lossy, so round-tripping one through a VDJAM-style similarity table would not reproduce what was fitted. _substitution_costs now takes a cost-table path directly, and reference-rule-scores threads --matrix through to the kernel and the V-loop tables. Built-in names are unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A 190-entry cost table carries no information about B, Z, X or *, and those rows are not inert: the V-loop templates are padded with characters outside the standard twenty, so charging padding the maximum cost silently rewrote the V-loop distances feeding the kernel. Measured effect: the same BLOSUM62 scored 0.1635790 through the built-in name and 0.1617818 through a cost-table path, a 0.0018 macro AP gap on human TRB that has nothing to do with the matrix under study. The padding rows are now taken from BLOSUM62's own Gram penalties, rescaled to the table's range, so every candidate handles padding identically and a table written from BLOSUM62 reads back as BLOSUM62 across all 24 symbols. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…profile
`vdjmatch --matrix` consumes a position-independent cost table: `match/geometry.py`
builds one 24x24 matrix and the native engine applies it at every column. A positional
hypothesis — the 21 October structural contact mask, which scores +0.0008926 macro AP
over a flat matrix on a fixed candidate tensor, 2.7x the shipped geometry's +0.0003280
and 13x BLOSUM62's +0.0000672 — therefore could not be scored on the real one-vs-all
retrieval metric. This closes that gap.
ROUTE: the native path carries it; no rescoring layer, no seqtree change.
`seqtree.gapblock` already takes `position_weights_by_length`, a mapping from full
sequence length to one nonnegative integer weight per position, indexed in the longer
sequence's frame and applied in `src/gapblock.cpp:56,65` as `pen[a][b] * weights[j]`.
It is plumbed through `grouped_count_batch` and `grouped_paired_sum_count_batch`
(`_alpha`/`_beta` suffixes) in released seqtree 1.3.1, the pinned floor, so no C++
rebuild is involved. vdjmatch already drives that knob: `_total_distance_options`
builds the uniform table `[0]*3 + [100]*(L-3-ctrim) + [0]*ctrim`, and
`grouped_total_distance_count_batch` overrides it for `significance-trimmed`. Weight
100 is unmodified, so a multiplier `m` is `round(100*m)` and an all-ones profile is
bit-for-bit the shipped geometry. A rescoring layer over `--kernel-only` output was
the fallback and is not needed.
WHAT THE MATRIX MAXIMUM MEANS ONCE COSTS VARY BY POSITION
`distance_units(costs)` returns `(max(costs), gcd(costs))` and the search radius
`22.5*scale`, the `ladder_cutoffs(scale)` bin edges and the CDR3 gap cost `600*scale`
all derive from it. A profile does not change `max(costs)`; it changes the ceiling at
one *column*, to `max(costs)*max(multiplier)`. Handled explicitly rather than silently
compensated:
- gaps stay unweighted at `600*scale`, so multipliers above 1 make gaps relatively
cheaper than the substitutions they compete with;
- radius and ladder stay matrix-derived, so a profile with mean multiplier above 1
retrieves strictly fewer neighbours at the same radius. The contact mask (x64 above
propensity 0.203) has mean multiplier ~13.6 and must be normalised to mean 1 over
scored positions in the file to stay radius-comparable with the flat arm. The code
applies the file as written and never normalises for you;
- `score_quantum` is deliberately not set, unlike `significance-trimmed` which floors
to 600: flooring would break the all-ones identity;
- the manifest records the profile hash and its per-locus multiplier min/max/mean, so
an arm's effective scale is recoverable from its own output.
Native rejects `longest * max(pen) * max(weight) > INT32_MAX`; headroom on BLOSUM62 at
junction length 25 is ~4.7e3 on the multiplier, so x64 is far inside it.
PROFILE FILE
TSV, header `locus length position multiplier`. `length` is the full junction length
(anchors included, the frame the search keys use) and `position` is 0-based in it.
Every position of every `(locus, length)` present in the batch must appear or the run
fails naming the gap — a partially covered profile would otherwise silently mix two
geometries. The trimmed flanks (first 3, last `ctrim`) carry weight 0 whatever the
profile says there, exactly as the uniform table does.
SCOPE
`--position-profile` is accepted only with `--kernel-only` and defaults to None; no
shipped geometry, default or calibration changes. Scoring one arm with a profile:
python -m vdjmatch.cli reference-rule-scores \
--vdjdb reference.zip --species human --locus TRB \
--kernel-only --matrix blosum62 \
--position-profile contact-mask.tsv \
--reference-population curated-benchmark --min-records 30 \
--threads 8 --output-prefix arm-contact
COST OF ONE ARM
Recorded full human TRB kernel-only arms on this cohort (108,390 queries, 252 eligible
peptide targets, `runs/vdjam-20261021-v1/cli/*.manifest.json` in the analysis repo):
619.665 s at 8 threads (BLOSUM62), and at 5 threads 716.156 s (VDJam-B), 716.850 s
(BLOSUM62), 786.872 s (clipped), 791.824 s (Miyazawa-Jernigan), 1059.973 s (PAM250).
The flag adds no traversal and no second pass — one integer row per junction length,
built once before the native call — so a profiled arm costs what its matrix costs,
less whatever the smaller ball saves when the mean multiplier exceeds 1.
Tests: `tests/unit/test_position_profile.py` (10) covers the all-ones bitwise identity
on single and paired axes at 1 and 4 threads, the uniform-table equality of the weight
builder, an exact directional check — multipliers 0, 2 and 4 on the one mismatched
column of a known pair give exactly 0x, 2x and 4x its 1800-unit native distance — and
the coverage, header, duplicate and sign rejections. The CLI case in
`test_reference_rule_scores_cli.py` runs the real command three ways and asserts the
all-ones arm reproduces the unprofiled scores while a heavy profile does not.
1547 non-slow unit tests pass in 92.48 s; strict Sphinx builds with no warnings.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…n lengths The per-(length, position) profile format was not estimable at the tails and was incomplete by construction. Valid human TRB carries 105,677 source records over 29 distinct junction lengths; lengths 13-16 hold 81,944 of them (77.5%), while 12 lengths hold fewer than 30 records each (50 records, 0.047%) and six -- 27, 31, 33, 34, 38 and 126 -- hold exactly one. A dense table asks for 126 multipliers for a junction seen once, and `profile_weights` hard-failed on any length the file did not enumerate, so those singletons would have had to be filled in or the run would die. The file is now end-anchored, like the `position_significance.tsv` the repository already ships for the same reason: rows are `locus side offset multiplier`, with `side` either `cys` (offsets from the conserved Cys104) or `fw` (offsets counted back from Phe/Trp118), each side contiguous from offset 0 and both required for a declared locus. The last declared offset is the pooled core and every deeper position reuses it. Each junction position is owned by its NEARER anchor, so the two counts meet in the middle and never overlap; once `length-1` exceeds the two declared depths, every interior position carries its nearer side's pooled core, which is the sense in which the two counts converge. The exact centre of an odd-length junction is equidistant and averages the two sides. Every length is then covered without enumeration, including the ones no cohort can estimate. If the two sides should instead combine by `min`, matching the `1 - max(p_same_V, p_same_J)` convention of `significance_weights`, that is the one branch in `profile_multipliers`. Default unchanged: `--position-profile` stays opt-in and off, and no shipped geometry, default or calibration moves. 1548 non-slow unit tests pass in 92.77 s; strict Sphinx builds with no warnings. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`significance_weights` already resolves a doubly-claimed column with `omega = 1 - max(p_same_V, p_same_J)`: a position is uninformative if it sits near *either* germline anchor, and neither anchor needs the other's agreement to say so. Stated about costs rather than proximities that is `min`, so `profile_multipliers` now takes the minimum of the two sides instead of handing the column to the nearer anchor. This is not a corner case. With the nine offsets per side the shipped `position_significance.tsv` declares, both anchors reach every junction of length <= 17, which is 101,093 of the 105,677 valid human TRB source records (95.7%), so the rule governs almost all of the cohort. The two counts still converge: past the two declared depths every interior column carries the minimum of the two pooled cores, so no length is special and none has to be enumerated. A contact mask is the shape this rule suits -- low at both anchors, high core on both sides -- because `min` then keeps the flanks cheap and lifts only the columns both anchors call deep. Default unchanged; `--position-profile` stays opt-in and off. 1548 non-slow unit tests pass in 90.17 s; strict Sphinx builds with no warnings. 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
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.
Add opt-in M4 exploration through
reference-rule-scores --kernel-only --position-profile PROFILE.tsv, using end-anchored multipliers in the existing native batch search. Include its explicit cost-table prerequisites from the matrix branch. No fitted profiles/matrices, scoring defaults, calibration changes or version bump are introduced.The integration preserves the released 0.6.2 alias fix. It corrects stale profile-format help, freezes and hash-binds explicit matrix files alongside profiles, rejects matrix overrides outside kernel-only mode, and preserves built-in matrix precedence over same-named local files. Documentation and the packaged skill describe the experimental scope. Mean-one normalization preserves nominal scaling; only an all-ones profile guarantees exact flat-geometry equivalence.
Validation on the merged tree: 1,552 non-slow tests passed in 75.63 s under the 120-second limit; 75 focused geometry/profile/cost-table/CLI tests passed; strict Sphinx passed. Coverage includes all-ones equivalence, paired/threaded native counts, matrix round trips, file-mutation rejection and unchanged built-in defaults. The final documentation/manifest wording also passed 31 focused tests in 2.81 s and strict Sphinx. No scientific analyses were rerun.