Repository navigation
v4.4.0 — the per-item-loop audit
Every module in src/mir read for one pattern: Python iterating over data where a batched call, a
polars expression or a table gather does the same work. Two sites had it. Everything else
iterates over structure — loci, samples, donors, prototypes, PCs, channels, columns, epochs,
minibatches — which is a different thing, and each step is a real numpy or torch call.
Requires vdjtools >= 4.6.0 as a hard floor: strip_allele_values does not exist before it.
shm_penalty_batch said "vectorised" and was a per-row loop
It called the scalar shm_penalty once per clonotype, which re-parsed the mutation string and made
one seqtree penalty call per substitution — so a BCR repertoire of 10^5 clonotypes carrying
~15 mutations each made ~1.5 million calls to answer 676 distinct questions. It sits on
TCREmp.embed, so it ran for every frame carrying v_mutations or v_identity.
| cohort | was | now | gain |
|---|---|---|---|
| 10,000 x 15 mutations | 32.7 ms | 10.0 ms | 3.3x |
| 100,000 x 15 mutations | 324.3 ms | 83.2 ms | 3.9x |
| 200,000 x 18 mutations | 748.2 ms | 210.0 ms | 3.6x |
| 200,000, identity only | 87.2 ms | 2.0 ms | 43.4x |
shm_penalty stays the scalar reference and the new tests/test_shm.py asserts the two agree with
==, not a tolerance: every BLOSUM62 penalty is a small integer, so a float64 sum is exact however
it is ordered, and a tolerance would hide a mis-indexed table — the one bug this rewrite could
plausibly introduce. This module had no tests at all, which is the other half of why its
docstring and its body disagreed for as long as they did.
Three behaviours the rewrite had to keep, now pinned: J, O and U are letters BLOSUM62 does not
define and must raise in both paths (a table with a silent default would return a finite,
plausible penalty instead); a non-letter token is dropped in silence (23V, A23, AB); and
truthiness is read off the raw value, so an int 0 takes the identity path while 3 does not —
casting to string first makes "0" truthy and scores that row from the wrong branch.
The isotype channels spent their time stripping allele suffixes, not matching bands
isotype_masks and isotype_shares each materialised the allele-stripped c_call column into a
Python list and ran one in test per clonotype per band. The membership tests were not the
cost. Profiled on 200,000 rows: strip_allele 100.7 ms, the three is_in tests 3.7 ms, the
to_list() 5.0 ms, the Python in tests 25.9 ms. Replacing only the comprehensions moved
132.2 ms to 134.8 — nothing, which is how the real cause surfaced.
| clonotypes | isotype_masks |
isotype_shares |
|---|---|---|
| 10,000 | 8.63 -> 1.20 ms (7.2x) | 7.66 -> 1.32 ms (5.8x) |
| 50,000 | 30.85 -> 2.29 ms (13.5x) | 32.80 -> 2.88 ms (11.4x) |
| 200,000 | 132.23 -> 6.10 ms (21.7x) | 124.89 -> 8.39 ms (14.9x) |
Output is identical on plain calls, nulls, comma ambiguity, surrounding whitespace, a single
occupied band, a band under the floor, and a frame with no c_call. isotype_shares takes its
masks from isotype_masks now rather than from a second copy of the same comprehension — two
copies is how the two came to be able to disagree about a null.
On reaching for C++
mirpy is pure Python by design (hatchling, no extension), so "make it native" means reaching into
seqtree's or vdjtools' C++. Both fixes went the other way and the numbers are the reason: a C++
mutation-spec parser would answer 676 distinct questions, so the (26, 26) lookup table is the
native code; and a native strip_allele over every row would still be a pass over every row — the
win there is doing the work 261 times instead of a million.
No API change. 380 tests pass; ruff and sphinx -W clean. The seven sites measured and
deliberately left alone are recorded in the repo's working notes, with their numbers, so nobody
re-derives them. Full detail in
CHANGELOG.md.