Pick the DP score width from bounds on the scores; check it, retry at 32 bits by glennhickey · Pull Request #14 · ComparativeGenomicsToolkit/abPOA

glennhickey · 2026-10-12T01:07:13Z

With large scores, abPOA can stop with Error in cg_backtrack or, worse, return a wrong MSA without any error. Cactus hit this on a CHM13 HPRC build: chr2's trained scoring model reached abPOA as -O 1290,3870 -E 60,2 (match up to 210, mismatch down to −1020), and cactus_consolidated(chr2) died 9 times out of 9. This PR replaces the 16/32-bit choice with one worked out from bounds on the scores. It checks any alignment whose width it cannot prove safe, and redoes a failed 16-bit alignment at 32 bits. Output is unchanged wherever the old code was right.

The bug:

The 16-bit gate (max_score <= INT16_MAX - min_mis - gap_oe1 - gap_oe2) budgets the range as if the "minus infinity" floor sat at INT16_MIN + max(min_mis, oe1, oe2). Upstream Backtrack errors when trying different scoring matrices yangao07/abPOA#79 (add820b, Oct 2024) raised the floor by 512*max(e1, e2) but did not take that out of the gate's budget.
With a gap extension of 60 the floor is +1824, above real scores. Floor values then win maxes, and the backtrack either cannot retrace the DP (Error in cg_backtrack) or retraces a wrong path.
The 512 was not derived. Floor-derived values fall below the floor by amounts that scale with the mismatch and gap-open scores, not the extension. In banded runs, cells no path reaches can also climb above the floor by the match score per diagonal step. So 512*e wastes much of the 16-bit range and can be too small at 32 bits (AVX-512 with trained scores).
Measured with the chr2 model on 2,327 windows (the yeast, chr10 and short windows below): 36110ac aborted on 1,030. It finished 1,297, and only 437 of those had the correct MSA (the result of the same model with every score halved).
This is the fourth abPOA score-overflow fix cactus has needed (after 2021, Jan 2024 and Oct 2024). The earlier ones each adjusted a constant or patched one case; this one derives the bounds.

What changes (src/abpoa_align_simd.c, Makefile):

Width rule: simd_abpoa_score_range() works out, per width and in 64-bit arithmetic, how far floor-derived values can move from the floor:
- down: mismatches along a diagonal, one gap open and one SIMD_SET_F step;
- up, when the band is narrower than the query: matches along a diagonal of unreachable cells, which stay in the band's first SIMD block.
It puts the floor just high enough that nothing wraps. It uses 16 bits only if every real score stays above the zone that floor-derived values can reach.
- The rule covers linear, affine and convex gaps, global, local and extend modes, any band, and path scores.
- M and X are recomputed from the matrix rather than taken from abpt.
Proven or checked: when every row covers the whole query, the lower bound on real scores is rigorous, and the width is proven safe. With a narrower band the bound is only an estimate, so the alignment is checked.
Checks:
- a failed backtrack sets a flag instead of calling err_fatal;
- the best score must lie in the valid range;
- for a width that is not proven, one pass after the DP (out of line) checks every stored H. It fails on a value in the zone just above the floor, or on one outside the range from which a later step could wrap.
If none of these fires, the result is the exact DP.
Retry: a 16-bit alignment that fails a check is redone at 32 bits. Before the retry, res and the per-node band positions (node_id_to_max_pos_left/right) are restored.
32 bits:
- When the bounds prove it, it runs unchecked. Otherwise it runs with every check on, and abPOA stops with a message naming the scores only if a check fails.
- abPOA refuses up front only when the checks themselves cannot fit in 32 bits: 3X + M + 3OE + (pn/2)*EM + 5PS >= 2^31.
Row 0: each gap piece's -(o + e*i) is computed in 64 bits and clamped at the floor. Before this, a first gap extension larger than the match score (346 against 210, say) could wrap the dearer piece in a 16-bit run counted as proven, and give a wrong MSA with no check firing. At 32 bits the same products could overflow int.
Loop alignment (e2aa0f1, a separate commit that can be dropped): on x86 the DP object is built with -falign-loops=32 if the compiler accepts it. On Skylake-class CPUs (Intel's JCC erratum), the DP's speed otherwise depends on where unrelated edits happen to place its loops. Without it, this code took 15% more cycles than 36110ac on a Skylake i7-6700, although it executes only 2% more instructions.
Other: the public API is unchanged. Error messages give the M and X the rule used, and a failure at a 32-bit width the bounds prove says the scores are in range, since it cannot be an overflow. At debug verbosity a score width: N bits line is printed. A test-only build (-DABPOA_TEST_SCORE_RANGE) reads environment variables that force 32 bits or a 16-bit floor that loses.

Results (MSA and consensus vs 36110ac; cactus's settings: convex gaps, global mode, -b 1000 -f 0.1, progressive; AVX2 build, as cactus builds abPOA):

windows	count	score sets	identical to `36110ac`
yeast chrI windows dumped from cactus	462	D, T1, T2h	1,386 / 1,386
yeast chrIV windows dumped from cactus	1,420	D, T1, T2h	4,260 / 4,260
HPRC haplotypes cut at GRCh38 chr10 repeat and indel sites	285	D, T1, T2h	855 / 855
HPRC-like chr9 windows (up to 10 kb, 464 rows)	15	D, T1, T2h	45 / 45
short windows built to drive scores low, plus 10 reproducers	160	D, T1, T2h	480 / 480

D is cactus's default (-O 400,1200 -E 30,1, match ≤ 100). T1 is a trained model with gap extension 30 (-O 1410,4230 -E 30,1, match ≤ 210). T2h is the chr2 model with every score halved (-O 645,1935 -E 30,1).

The chr2 model itself (-E 60,2) on the 2,327 yeast, chr10 and short windows: this PR gives the T2h result on all 2,327. 36110ac aborted on 1,030 and gave a wrong MSA on 860.
Forced widths: forcing 32 bits, or forcing a 16-bit floor that loses (+1, +2000), gives the same output on the yeast, chr10 and short windows. The losing floors are caught and redone.
Other builds: the default (runtime-dispatch), SSE2 and SSE4.1 builds match on the chr10 and short windows.
Fuzzing:
- About 3,000 cases ran on the first version of this branch: 1,000 random cases (all gap and alignment modes, bands, scores up to about 1e7), 800 with forced floors, 250 scale sweeps and 517 real windows. They found the row-0 wrap (fixed here) and refusals of very large scores (below).
- The final version, re-run on the 1,000 random and 517 real-window cases, never differs from or fails where the first version completed, except where the first version was wrong.
Very large scores (about 1e5 and up): the first version refused 28 fuzz inputs at such scores and stopped on 1 more.
- 18 now give the correct result.
- 5 stop when a check fails. Their real scores do not fit in 32 bits; 36110ac happened to be right on 1 of them.
- 6 are refused up front (gap opens of 0.68e9 to 1.7e9, a mismatch of 0.72e9, or input scores beyond 32 bits); 36110ac happened to be right on 3 of them.
Speed: on a 27-sequence HPRC-like window that runs entirely at 32 bits, old and new ran at the same time, pinned to the two hyperthreads of one core, and user cycles were counted. Measured that way, the alignment takes 0.997 of 36110ac's cycles (cactus's AVX2 build, D and T1, 10 pairs each; old against old gives 1.001). It executes 2% more instructions. The runtime-dispatch build measured 0.998.
More 16-bit work: the new rule sends more alignments to 16 bits than the old gate (yeast windows 96% → 97–98%). Runtime barely moves, because nearly all DP cells are in large windows that run at 32 bits either way.
Memory: peak memory is unchanged.

Tests:

tests/run_all.sh has a new regression block: test_data/seq.fa with the chr2 scores must equal the same run with every score halved, under three band settings. 36110ac aborts on it.
tests/test_score_range.sh is a new file, which run_all.sh picks up with the other test_*.sh scripts; it takes about 2 minutes. It builds abPOA with -DABPOA_TEST_SCORE_RANGE and checks:
- that 16 and 32 bits agree in every gap and alignment mode, banded and not, and with amb-strand, path scores and progressive;
- that floors forced to lose are caught and redone;
- both sides of the 16-bit query-length and graph-depth limits;
- the chr2 model against forced 32 bits and against itself halved;
- the row-0 input, which the first version of this branch got wrong;
- that cactus's default scores × 10000 on a banded 1.2 kb window give the unscaled result;
- that a mismatch of 8e8 is refused.
tests/run_all.sh passes 12/12, in about 2 minutes. The new checks fail on 36110ac.
The SSE2, SSE4.1, AVX2 and dispatch builds compile with no warnings. ASan and UBSan are clean on the retry path.
No AVX-512 hardware was available, so the 32-lane path is covered by the formulas only.

Merging: this targets main, whose tip is 36110ac, the commit cactus pins today. The fork's default branch, master, is the July 2021 v1.2.5 merge and is not what cactus uses. Please merge with a merge commit, not a squash or rebase: cactus will pin this branch's tip, 12f99c4, which must stay reachable.

🤖 Generated with Claude Code

https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

… bits The 16-bit gate budgeted the DP's range as if "minus infinity" sat at INT16_MIN + max(mismatch, oe1, oe2), but add820b (yangao07#79) put inf_min 512*max(e1,e2) higher without taking that out of the budget. With large scores the floor reached real scores: a last-train model scaled for cactus (match <= 210, mismatch >= -1020, -O 1290,3870 -E 60,2) put it at +1824, floor values won maxes, and the backtrack aborted with "Error in cg_backtrack" or, when it did not, could retrace a wrong alignment. On 150 short windows with that model 36110ac aborted on 81 and gave a different MSA from the 32-bit DP on 9. simd_abpoa_score_range() now works out, per width and in 64 bits, how far floor-derived values can move from the floor (down: mismatches along a diagonal, a gap open, one SIMD_SET_F step; up, when the band is narrower than the query: matches along a diagonal of unreachable cells, which stay in the band's first SIMD block), puts the floor just high enough that nothing wraps, and uses 16 bits only if every real score stays above the zone where a floor-derived value could reach it. When every row covers the whole query that is proven (no real H is below -(g(D) + g(qlen)): delete down to the row, insert across); with a narrower band it is an estimate. 32 bits uses bounds that need no band geometry. This holds for every gap and alignment mode and any non-negative scores; max_mat/min_mis are recomputed from the matrix. Checks: the backtrack flags a failure instead of aborting; the best score must be in range; and an alignment whose width is only estimated safe (16 bits with a band narrower than the query, or 32 bits on graphs deeper than about 1.5M rows for a 10 kb query with trained scores) has every stored H checked, three SIMD ops per block, for a value in that zone or above qlen*max_match. A run that passes is the exact DP. If a check fails, a 16-bit alignment is redone in 32 bits after restoring the result and the band bounds (node_id_to_max_pos_left/right, through an undo log) it changed; a 32-bit one aborts with a message naming the scores and sizes. No 64-bit path: for trained scores 32 bits is proven up to graphs whose banded DP would need tens of GB. Public API unchanged. -DABPOA_TEST_SCORE_RANGE (tests only) lets the environment force 32 bits or a losing 16-bit floor. Output is unchanged wherever 36110ac was right: byte-identical MSAs with cactus's default and trained (E1=30) scores and abPOA's defaults, and identical to a forced 32-bit DP for SSE2, SSE4.1, AVX2 and dispatch builds and for 16-bit floors forced high enough to lose (all caught and redone). The gate admits more alignments to 16 bits than before (two gaps at the cheaper piece, not one at e1), all of them proven. tests/run_all.sh: test_data/seq.fa with those scores must equal the same run with every score halved, under three band settings (36110ac aborts). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

tests/test_score_range.sh builds abPOA with -DABPOA_TEST_SCORE_RANGE and checks that 16- and 32-bit DP agree in every gap and alignment mode, banded and not; that a 16-bit floor forced to lose is caught and redone at 32 bits; that both sides of the 16-bit query-length and graph-depth limits give the 32-bit result; and that a trained model with gap extension 60 gives the same result as with every score halved. run_all.sh picks it up with the other test_*.sh scripts. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

The zone scan for widths that are not proven safe ran inside the DP row loop. It never runs for cactus's alignments, but its code alone made GCC compile the hot loop differently: on a window that runs entirely at 32 bits, f46f6c6 took about 3-4% more cycles than 36110ac. The scan now runs once, after the DP, in a separate non-inlined function per width. Every row's band is still in memory for the backtrack, so it reads the same stored H values the in-loop scan saw: row 0 and every row computed, skipping rows outside the subgraph and, in extend mode, rows after a z-drop stop. On pinned, simultaneous old/new pairs the cycle count is back to 36110ac's. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

Outside local mode, row 0 holds a gap along the query on each gap piece, -(o_k + e_k*i). The width choice bounds only the cheaper piece (the real score), so with convex gaps whose first extension exceeds the largest match score, the other piece left the 16-bit range in a run counted as proven. The best-score check happened to catch it and redo the alignment at 32 bits (-M 1 -X 4 -O 4,24 -E 100,1 -b -1 with a 400 bp query), but nothing guaranteed that. At 32 bits the same products could overflow int. Each piece is now computed in 64 bits and clamped at inf_min. Row 0's F values are not used after its H is set (the next row recomputes F and the backtrack stops at row 0), and H, the cheaper piece, is above the floor whenever the run is valid, so output does not change for any input the old code handled correctly. Adding a per-piece term to the width test instead would send these alignments to 32 bits for nothing. tests/test_score_range.sh: that input must give the 32-bit result at 16 bits without a retry. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

The score-range error messages printed abpt->max_mat and min_mis, which the width rule does not trust: a library caller that sets its own matrix may leave them stale. Print the M and X that simd_abpoa_score_range() recomputes from the matrix instead. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

…sing them When even the 32-bit bounds did not fit (qlen*(M+X) near 2^31, or qlen*M itself beyond 2^31), abPOA stopped before aligning with "scores too large". A fuzz campaign found 28 such inputs that 36110ac aligned correctly. Now a 32-bit alignment the bounds do not prove is run with every check on, and abPOA stops only if a check fails: - The floor sits as low as the band's first SIMD block allows (pn diagonal steps for floor-derived values, not qlen+1), leaving real scores the most room; floor-derived values that go further fail the checks rather than go unseen. This also applies to 32-bit runs the old rule checked with the qlen+1 allowance, where real scores could reach that higher zone although the result was right (fuzz case 1129). - The scan also requires every stored H to be at least low, the floor for a full band: from there no later step can wrap below INT_MIN. With that, a wrap no longer needs qlen*M + step <= INT_MAX to be seen, so where qlen*M does not fit, the highest H accepted is INT_MAX - M, from which no step can wrap. - Only scores whose checks cannot fit in 32 bits at all are refused up front: 3X + M + 3OE + pn/2*EM + 5PS >= 2^31, in 64 bits. Inputs the bounds do prove are unchanged: same width, floor and output. Of the 28 fuzz inputs (and 1129), 18 now give the correct output (17 equal 36110ac, confirmed by its forced-32-bit run and, where scaled, by the unscaled scores; on 1559, 36110ac's row 0 overflowed int and the new output equals the unscaled one). Five fail a check, all with real scores beyond 32 bits (local best near 3.5e9 and 1e10, in-band scores to -3.7e9 and -7.5e9, H to 1.6e10); on one of them 36110ac's output was still right. Six are refused: a gap open of 0.68e9 to 1.7e9, a mismatch of 0.72e9, or input scores beyond 32 bits. tests/test_score_range.sh: cactus's default scores times 10000 on a banded 1.2 kb window must give the result of the scores themselves, and a mismatch of 8e8 must be refused. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

A fuzz campaign on real windows found that, before 0f6d155, models whose first gap extension exceeds the match score (346 or 400 against 210, trained-like matrices) gave a different MSA from the 32-bit DP with no check firing: row 0's dearer gap piece wrapped in a 16-bit run counted as proven. Divergent generated sequences reproduce it (16 of 60 inputs with -M 210 -X 1020 -O 871,2479 -E 346,2); the one added to tests/test_score_range.sh must now give the 32-bit result. f46f6c6, 7312421 and 8eeb089 fail it, at both SIMD widths. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

With the score-range changes the DP executes the same instructions in its loops as before (about 2% more in total), yet on this Skylake its cycle count moved with unrelated edits: 36110ac and 8eeb089..8d41b30 were level, while 1a682ec, whose only effect on the loops is to shift them by a few bytes, took 14.5% more cycles in the AVX2 build (the one cactus uses) and 5.5% more in the runtime-dispatch build. The cause is Intel's JCC erratum: a conditional jump of the scalar max_in_row loop now crossed a 32-byte boundary, which keeps the loop out of the decoded-uop cache. Built with -Wa,-mbranches-within-32B-boundaries, all three are level. Aligning loops to 32 bytes fixes where each loop sits relative to those boundaries, so its speed no longer depends on the code before it. On pinned, simultaneous pairs on a pure 32-bit window, the DP is now level with 36110ac: cycles 0.998 of it (AVX2 and dispatch, cactus's default and a trained model). The flag is added to the DP object only, on x86, and only if the compiler takes it (it is probed with -Werror). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

When the bounds prove a 32-bit alignment's scores fit, a failed check cannot be an overflow, yet a failed best-score check (with or without a failed backtrack) still printed "Scores this large cannot be aligned exactly: scale them down", even for small scores (-M 89 -X 2279 -O 7100,2696 -E 110,52 -b 0 -f 0 -S on a 51 bp query, which also kills 36110ac with "Error in cg_backtrack"). Any failure at a proven 32-bit width now says the scores are in range, as a backtrack-only failure already did. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV

glennhickey and others added 9 commits October 11, 2026 13:08

This was referenced Oct 12, 2026

Fix abPOA's cg_backtrack crash with large trained scores; divide trained models by their common factor ComparativeGenomicsToolkit/cactus#2040

Merged

Update abPOA to v1.5.7 with the 16/32-bit score width fix ComparativeGenomicsToolkit/taffy#95

Merged

glennhickey merged commit c7bdc84 into main Oct 12, 2026
8 checks passed

Navigation Menu

Repository navigation

Pick the DP score width from bounds on the scores; check it, retry at 32 bits - #14

Pick the DP score width from bounds on the scores; check it, retry at 32 bits#14
glennhickey merged 9 commits into
mainfrom
robust-score-width

glennhickey commented Oct 12, 2026

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

glennhickey commented Oct 12, 2026

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant