From f46f6c6f3139f5a6019424845f6abf421dc8f742 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 12:52:12 -0400 Subject: [PATCH 1/9] Pick the score width from bounds on the scores; check it, retry in 32 bits The 16-bit gate budgeted the DP's range as if "minus infinity" sat at INT16_MIN + max(mismatch, oe1, oe2), but add820b (#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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- src/abpoa_align_simd.c | 209 +++++++++++++++++++++++++++++++++++++---- tests/run_all.sh | 15 +++ 2 files changed, 205 insertions(+), 19 deletions(-) diff --git a/src/abpoa_align_simd.c b/src/abpoa_align_simd.c index 7ee8917..e8bcfce 100644 --- a/src/abpoa_align_simd.c +++ b/src/abpoa_align_simd.c @@ -181,15 +181,17 @@ } \ } \ } \ - if (hit == 0) err_fatal_simple("Error in lg_backtrack."); \ + if (hit == 0) { aln_fail |= ABPOA_FAIL_BT; break; } /* see simd_abpoa_check_scores */ \ simd_output_pre_nodes(pre_index[i], pre_n[i], i, j, 0, abpt->verbose); \ } \ + if (aln_fail) { free(cigar); } else { \ if (j > 0) cigar = abpoa_push_cigar(&n_c, &m_c, cigar, ABPOA_CINS, j, -1, j-1); \ /* reverse cigar */ \ res->graph_cigar = abpt->rev_cigar ? cigar : abpoa_reverse_cigar(n_c, cigar); \ res->n_cigar = n_c; res->m_cigar = m_c; \ res->node_e = abpoa_graph_index_to_node_id(graph, best_i+beg_index), res->query_e=best_j-1; /*0-based*/ \ res->node_s = abpoa_graph_index_to_node_id(graph, _start_i+beg_index), res->query_s=_start_j-1; \ + } \ /*abpoa_print_cigar(n_c, *graph_cigar, graph);*/ \ } @@ -294,15 +296,17 @@ } \ } \ } \ - if (hit == 0) err_fatal_simple("Error in ag_backtrack."); \ + if (hit == 0) { aln_fail |= ABPOA_FAIL_BT; break; } /* see simd_abpoa_check_scores */ \ simd_output_pre_nodes(pre_index[i], pre_n[i], i, j, cur_op, abpt->verbose); \ } \ + if (aln_fail) { free(cigar); } else { \ if (j > 0) cigar = abpoa_push_cigar(&n_c, &m_c, cigar, ABPOA_CINS, j, -1, j-1); \ /* reverse cigar */ \ res->graph_cigar = abpt->rev_cigar ? cigar : abpoa_reverse_cigar(n_c, cigar); \ res->n_cigar = n_c; res->m_cigar = m_c; \ res->node_e = abpoa_graph_index_to_node_id(graph, best_i+beg_index), res->query_e=best_j-1; /*0-based*/ \ res->node_s = abpoa_graph_index_to_node_id(graph, _start_i+beg_index), res->query_s=_start_j-1; \ + } \ /*abpoa_print_cigar(n_c, *graph_cigar, graph);*/ \ } @@ -445,15 +449,17 @@ } \ } \ } \ - if (hit == 0) err_fatal_simple("Error in cg_backtrack."); \ + if (hit == 0) { aln_fail |= ABPOA_FAIL_BT; break; } /* see simd_abpoa_check_scores */ \ simd_output_pre_nodes(pre_index[i], pre_n[i], i, j, cur_op, abpt->verbose); \ } \ + if (aln_fail) { free(cigar); } else { \ if (j > 0) cigar = abpoa_push_cigar(&n_c, &m_c, cigar, ABPOA_CINS, j, -1, j-1); \ /* reverse cigar */ \ res->graph_cigar = abpt->rev_cigar ? cigar : abpoa_reverse_cigar(n_c, cigar); \ res->n_cigar = n_c; res->m_cigar = m_c; \ res->node_e = abpoa_graph_index_to_node_id(graph, best_i+beg_index), res->query_e=best_j-1; /*0-based*/ \ res->node_s = abpoa_graph_index_to_node_id(graph, _start_i+beg_index), res->query_s=_start_j-1; \ + } \ /*abpoa_print_cigar(n_c, *graph_cigar, graph);*/ \ } @@ -475,6 +481,8 @@ int w = abpt->wb < 0 ? qlen : abpt->wb+(int)(abpt->wf*qlen); /* when w < 0, do whole global */ \ int best_i = 0, best_j = 0, best_id = 0, max, left_max_i=-1, right_max_i=-1; \ SIMDi zero = SIMDSetZeroi(), SIMD_INF_MIN = SIMDSetOne(inf_min); \ + SIMDi ZONE_FMAX = SIMDSetOne((score_t)chk_fmax), ZONE_MIN = SIMDSetOne((score_t)chk_hi); \ + SIMDi H_MAX = SIMD_INF_MIN; /* simd_abpoa_zone_scan */ \ pn = sp.num_of_value; qp_sn = dp_sn = (matrix_col_n + pn - 1) / pn; \ log_n = sp.log_num, size = sp.size; qp = abm->s_mem; \ int set_num; SIMDi *PRE_MASK, *SUF_MIN, *PRE_MIN; \ @@ -598,11 +606,14 @@ #define simd_abpoa_lg_first_row { \ /* fill the first row */ \ if (abpt->wb >= 0) { \ + simd_abpoa_save_max_pos(beg_node_id); \ graph->node_id_to_max_pos_left[beg_node_id] = graph->node_id_to_max_pos_right[beg_node_id] = 0; \ for (i = 0; i < graph->node[beg_node_id].out_edge_n; ++i) { /* set max pos for out_id */ \ int out_id = graph->node[beg_node_id].out_id[i]; \ - if (index_map[abpoa_graph_node_id_to_index(graph, out_id)]) \ + if (index_map[abpoa_graph_node_id_to_index(graph, out_id)]) { \ + simd_abpoa_save_max_pos(out_id); \ graph->node_id_to_max_pos_left[out_id] = graph->node_id_to_max_pos_right[out_id] = 1; \ + } \ } \ dp_beg[0] = 0, dp_end[0] = GET_AD_DP_END(graph, w, beg_node_id, end_node_id, qlen); \ } else { \ @@ -615,11 +626,14 @@ #define simd_abpoa_ag_first_row { \ /* fill the first row */ \ if (abpt->wb >= 0) { \ + simd_abpoa_save_max_pos(beg_node_id); \ graph->node_id_to_max_pos_left[beg_node_id] = graph->node_id_to_max_pos_right[beg_node_id] = 0; \ for (i = 0; i < graph->node[beg_node_id].out_edge_n; ++i) { /* set max pos for out_id */ \ int out_id = graph->node[beg_node_id].out_id[i]; \ - if (index_map[abpoa_graph_node_id_to_index(graph, out_id)]) \ + if (index_map[abpoa_graph_node_id_to_index(graph, out_id)]) { \ + simd_abpoa_save_max_pos(out_id); \ graph->node_id_to_max_pos_left[out_id] = graph->node_id_to_max_pos_right[out_id] = 1; \ + } \ } \ dp_beg[0] = 0, dp_end[0] = GET_AD_DP_END(graph, w, beg_node_id, end_node_id, qlen); \ } else { \ @@ -633,11 +647,14 @@ #define simd_abpoa_cg_first_row { \ /* fill the first row */ \ if (abpt->wb >= 0) { \ + simd_abpoa_save_max_pos(beg_node_id); \ graph->node_id_to_max_pos_left[beg_node_id] = graph->node_id_to_max_pos_right[beg_node_id] = 0; \ for (i = 0; i < graph->node[beg_node_id].out_edge_n; ++i) { /* set max pos for out_id */ \ int out_id = graph->node[beg_node_id].out_id[i]; \ - if (index_map[abpoa_graph_node_id_to_index(graph, out_id)]) \ + if (index_map[abpoa_graph_node_id_to_index(graph, out_id)]) { \ + simd_abpoa_save_max_pos(out_id); \ graph->node_id_to_max_pos_left[out_id] = graph->node_id_to_max_pos_right[out_id] = 1; \ + } \ } \ dp_beg[0] = 0, dp_end[0] = GET_AD_DP_END(graph, w, beg_node_id, end_node_id, qlen); \ } else { \ @@ -1140,6 +1157,9 @@ /* set max_pos_left/right for next nodes */ \ for (i = 0; i < graph->node[node_id].out_edge_n; ++i) { \ int out_node_id = graph->node[node_id].out_id[i]; \ + if (right_max_i+1 > graph->node_id_to_max_pos_right[out_node_id] || \ + left_max_i+1 < graph->node_id_to_max_pos_left[out_node_id]) \ + simd_abpoa_save_max_pos(out_node_id); \ if (right_max_i+1 > graph->node_id_to_max_pos_right[out_node_id]) \ graph->node_id_to_max_pos_right[out_node_id] = right_max_i+1; \ if (left_max_i+1 < graph->node_id_to_max_pos_left[out_node_id]) \ @@ -1147,15 +1167,53 @@ } \ } +/* band bounds a 16-bit attempt changes, logged so that a retry at 32 bits starts from the same state */ +#define simd_abpoa_save_max_pos(id) { \ + if (mpos_log_on) { \ + if (mpos_n == mpos_m) { \ + mpos_m = mpos_m ? mpos_m << 1 : 256; \ + mpos_log = (int*)_err_realloc(mpos_log, mpos_m * 3 * sizeof(int)); \ + } \ + mpos_log[mpos_n*3] = (id); mpos_log[mpos_n*3+1] = graph->node_id_to_max_pos_left[id]; \ + mpos_log[mpos_n*3+2] = graph->node_id_to_max_pos_right[id]; ++mpos_n; \ + } \ +} + +/* only when the width is not proven safe (see simd_abpoa_score_range): every stored H of the row, + * min over those above fmax and max over all, tested once per alignment by simd_abpoa_check_scores */ +#define simd_abpoa_zone_scan(first_sn, last_sn, SIMDSetIfGreater, SIMDMax) { \ + SIMDi _t; \ + for (sn_i = first_sn; sn_i <= last_sn; ++sn_i) { \ + SIMDSetIfGreater(_t, dp_h[sn_i], ZONE_FMAX, dp_h[sn_i], ZONE_MIN); \ + SIMDSetIfGreater(ZONE_MIN, ZONE_MIN, _t, _t, ZONE_MIN); \ + H_MAX = SIMDMax(H_MAX, dp_h[sn_i]); \ + } \ +} + +/* a real score in (fmax, ztop] can lose to, or be confused with, a floor-derived one, and a score above + * u can only be a wrap: either way the alignment is redone at 32 bits (simd_abpoa_align_sequence_to_subgraph) */ +#define simd_abpoa_check_scores(score_t) { \ + if (chk_zone) { \ + score_t *_zmin = (score_t*)&ZONE_MIN, *_hmax = (score_t*)&H_MAX; \ + for (i = 0; i < pn; ++i) { \ + if (_zmin[i] <= chk_ztop || _hmax[i] > chk_u) aln_fail |= ABPOA_FAIL_ZONE; \ + } \ + } \ + if ((best_i > 0 || abpt->align_mode == ABPOA_GLOBAL_MODE) && (best_score > chk_u || best_score <= chk_ztop)) \ + aln_fail |= ABPOA_FAIL_BEST; \ +} + // TODO end_bonus for extension // linear gap penalty: gap_open1 == 0 #define simd_abpoa_lg_align_sequence_to_graph_core(score_t, sp, SIMDSetOne, SIMDMax, SIMDAdd, \ SIMDSub, SIMDShiftOneN, SIMDSetIfGreater, SIMDGetIfGreater) { \ simd_abpoa_lg_var(score_t, sp, SIMDSetOne, SIMDShiftOneN, SIMDAdd); \ simd_abpoa_lg_first_dp(score_t); \ + if (chk_zone) simd_abpoa_zone_scan(0, _end_sn, SIMDSetIfGreater, SIMDMax); \ for (index_i=beg_index+1, dp_i=1; index_ialign_mode == ABPOA_LOCAL_MODE) { \ simd_abpoa_max_in_row(score_t, SIMDSetIfGreater, SIMDGetIfGreater); \ set_global_max_score(max, dp_i, left_max_i); \ @@ -1178,7 +1236,8 @@ simd_abpoa_print_lg_matrix(score_t, beg_index, end_index); \ fprintf(stderr, "best_score: (%d, %d) -> %d\n", best_i, best_j, best_score); \ } \ - if (abpt->ret_cigar) simd_abpoa_lg_backtrack(score_t); \ + simd_abpoa_check_scores(score_t); \ + if (abpt->ret_cigar && aln_fail == 0) simd_abpoa_lg_backtrack(score_t); \ simd_abpoa_free_var; SIMDFree(GAP_E1S); \ } @@ -1187,9 +1246,11 @@ SIMDSub, SIMDShiftOneN, SIMDSetIfGreater, SIMDGetIfGreater, SIMDSetIfEqual) { \ simd_abpoa_ag_var(score_t, sp, SIMDSetOne, SIMDShiftOneN, SIMDAdd); \ simd_abpoa_ag_first_dp(score_t); \ + if (chk_zone) simd_abpoa_zone_scan(0, _end_sn, SIMDSetIfGreater, SIMDMax); \ for (index_i=beg_index+1, dp_i=1; index_ialign_mode == ABPOA_LOCAL_MODE) { \ simd_abpoa_max_in_row(score_t, SIMDSetIfGreater, SIMDGetIfGreater); \ set_global_max_score(max, dp_i, left_max_i); \ @@ -1211,7 +1272,8 @@ simd_abpoa_print_ag_matrix(score_t, beg_index, end_index); \ fprintf(stderr, "best_score: (%d, %d) -> %d\n", best_i, best_j, best_score); \ } \ - if (abpt->ret_cigar) simd_abpoa_ag_backtrack(score_t); \ + simd_abpoa_check_scores(score_t); \ + if (abpt->ret_cigar && aln_fail == 0) simd_abpoa_ag_backtrack(score_t); \ simd_abpoa_free_var; SIMDFree(GAP_E1S); \ } @@ -1220,9 +1282,11 @@ SIMDSub, SIMDShiftOneN, SIMDSetIfGreater, SIMDGetIfGreater, SIMDSetIfEqual) { \ simd_abpoa_cg_var(score_t, sp, SIMDSetOne, SIMDShiftOneN, SIMDAdd); \ simd_abpoa_cg_first_dp(score_t); \ + if (chk_zone) simd_abpoa_zone_scan(0, _end_sn, SIMDSetIfGreater, SIMDMax); \ for (index_i=beg_index+1, dp_i=1; index_ialign_mode == ABPOA_LOCAL_MODE) { \ simd_abpoa_max_in_row(score_t, SIMDSetIfGreater, SIMDGetIfGreater); \ set_global_max_score(max, dp_i, left_max_i); \ @@ -1244,10 +1308,72 @@ simd_abpoa_print_cg_matrix(score_t, beg_index, end_index); \ fprintf(stderr,"best_score: (%d, %d) -> %d\n",best_i, best_j, best_score); \ } \ - if (abpt->ret_cigar) simd_abpoa_cg_backtrack(score_t); \ + simd_abpoa_check_scores(score_t); \ + if (abpt->ret_cigar && aln_fail == 0) simd_abpoa_cg_backtrack(score_t); \ simd_abpoa_free_var; SIMDFree(GAP_E1S); SIMDFree(GAP_E2S); \ } +/* + * Score width. Cells no path reaches start at a "minus infinity" floor, and a floor-derived value moves away + * from it: down by at most `decay' (mismatches along a diagonal, a gap open, one log-step of SIMD_SET_F) and, + * when the band is narrower than the query, up by at most `growth' (matches along a diagonal of unreachable + * cells, which only exist in the band's first SIMD block). With + * floor = MIN + decay, fmax = floor + growth, ztop = fmax + M + 2X + 2OE + 3PS, u = qlen*M, + * nothing wraps as long as u plus the largest single subtraction fits under MAX, and no floor-derived value + * reaches a real one as long as every real H stays above ztop: the result is then the exact DP, identical at + * either width. Real H are at most u and, when every row covers the whole query, at least + * -(g(D) + g(qlen) + PS*D) (delete down to the row, insert across, each gap on its cheaper piece): this + * "proves" a width. For narrower bands that value is only an estimate (a rigorous bound is far lower), so a + * width the estimate allows is "feasible" and is used with every stored H checked: a real H in (fmax, ztop] + * is the first sign of the floor winning, and a value above u is a wrap, so a run that passes is exact and + * one that fails is redone at 32 bits. 32 bits takes bounds that need no band geometry. + * M/X: largest match score/mismatch penalty, OE/EM: largest gap open+extension/extension, PS: largest path + * score penalty (20, or 0 without inc_path_score), D: rows on a path (gn bounds it), g(n): cheapest gap of n. + */ +#define ABPOA_FAIL_BT 1 /* the backtrack could not retrace the DP */ +#define ABPOA_FAIL_ZONE 2 /* a stored H in (fmax, ztop] or above u */ +#define ABPOA_FAIL_BEST 4 /* the best score in (-inf, ztop] or above u */ + +typedef struct { int64_t floor, fmax, ztop, u; int fits, feasible, proven; } abpoa_score_range_t; + +static int64_t simd_abpoa_gap_cost(abpoa_para_t *abpt, int64_t n) { + if (n <= 0) return 0; + int64_t c = (abpt->gap_mode == ABPOA_LINEAR_GAP ? 0 : abpt->gap_open1) + (int64_t)abpt->gap_ext1 * n; + if (abpt->gap_mode == ABPOA_CONVEX_GAP) c = MIN_OF_TWO(c, abpt->gap_open2 + (int64_t)abpt->gap_ext2 * n); + return c; +} + +// geometry: 0: every row covers the whole query, 1: band of w >= 1, 2: any band +static abpoa_score_range_t simd_abpoa_score_range(abpoa_para_t *abpt, int bits, int pn, int qlen, int gn, int geometry) { + abpoa_score_range_t r; int i; + int64_t lo = bits == 16 ? INT16_MIN : INT32_MIN, hi = bits == 16 ? INT16_MAX : INT32_MAX; + int64_t M = 0, X = 0, Q = qlen, D = gn, PS = abpt->inc_path_score ? 20 : 0; + for (i = 0; i < abpt->m * abpt->m; ++i) { + M = MAX_OF_TWO(M, abpt->mat[i]); X = MAX_OF_TWO(X, -(int64_t)abpt->mat[i]); + } + int64_t OE = (abpt->gap_mode == ABPOA_LINEAR_GAP ? 0 : abpt->gap_open1) + (int64_t)abpt->gap_ext1; + int64_t OEmin = OE, EM = abpt->gap_ext1; + if (abpt->gap_mode == ABPOA_CONVEX_GAP) { + int64_t oe2 = (int64_t)abpt->gap_open2 + abpt->gap_ext2; + OE = MAX_OF_TWO(OE, oe2); OEmin = MIN_OF_TWO(OEmin, oe2); EM = MAX_OF_TWO(EM, abpt->gap_ext2); + } + int64_t chain = geometry == 0 ? X + PS : geometry == 1 ? (pn + 1) * (X + PS) : (Q + 2) * (X + PS); + int64_t growth = geometry == 0 ? 0 : geometry == 1 ? pn * M : (Q + 1) * M; + int64_t step = MAX_OF_THREE(X + PS, OE + PS, pn / 2 * EM); // largest single subtraction + r.u = Q * M; + r.floor = lo + chain + OE + pn / 2 * EM + PS; + r.fmax = r.floor + growth; + r.ztop = r.fmax + M + 2 * X + 2 * OE + 3 * PS; + int64_t low_est = -(simd_abpoa_gap_cost(abpt, D) + simd_abpoa_gap_cost(abpt, Q) + PS * D); + int64_t low = geometry == 0 ? low_est : -(MAX_OF_TWO(X, OEmin) * (Q + D) + PS * D); + if (abpt->align_mode == ABPOA_LOCAL_MODE) low_est = low = 0; + // ztop < 0: the checks can only miss the floor winning if the source's score of 0 is in or under the zone + r.fits = pn / 2 * EM <= hi && OE <= hi && M <= hi && X <= hi && Q < hi && r.ztop < 0 && r.u + step <= hi; + r.feasible = r.fits && low_est > r.ztop; + r.proven = r.feasible && low > r.ztop; + return r; +} + // align query to subgraph between beg_node_id and end_node_id (both are excluded) // generally: beg/end are the SRC/SINK_node #ifdef ABPOA_SIMD_DISPATCH @@ -1299,26 +1425,51 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b // fprintf(stderr, "reg_n: %d\n", simd_p16.reg_n); - int32_t max_score, bits, mem_ret=0, gap_ext1 = abpt->gap_ext1, gap_ext2 = abpt->gap_ext2; - int32_t gap_oe1 = abpt->gap_open1+gap_ext1, gap_oe2 = abpt->gap_open2+gap_ext2; - int len = qlen > gn ? qlen : gn; + int32_t bits, mem_ret=0, gap_ext1 = abpt->gap_ext1, gap_ext2 = abpt->gap_ext2; #ifdef __SIMD_DEBUG__ + int32_t gap_oe1 = abpt->gap_open1+gap_ext1, gap_oe2 = abpt->gap_open2+gap_ext2; simd_p32.inf_min = MAX_OF_THREE(INT32_MIN + abpt->min_mis, INT32_MIN + gap_oe1, INT32_MIN + gap_oe2) + 512 * MAX_OF_TWO(gap_ext1, gap_ext2); if (simd_abpoa_realloc(ab, gn, qlen, abpt, simd_p32)) return 0; if (abpt->gap_mode == ABPOA_CONVEX_GAP) abpoa_cg_global_align_sequence_to_graph_core(ab, beg_node_id, beg_index, end_node_id, end_index, index_map, qlen, query, abpt, _simd_p32, res); #else - max_score = MAX_OF_TWO(qlen * abpt->max_mat, len * abpt->gap_ext1 + abpt->gap_open1); - if (max_score <= INT16_MAX - abpt->min_mis - gap_oe1 - gap_oe2) { - simd_p16.inf_min = MAX_OF_THREE(INT16_MIN + abpt->min_mis, INT16_MIN + gap_oe1, INT16_MIN + gap_oe2) + 512 * MAX_OF_TWO(gap_ext1, gap_ext2); + // 16 bits where the scores provably (or, for narrow bands, probably) fit, else 32; see simd_abpoa_score_range + int bw = abpt->wb < 0 ? qlen : abpt->wb + (int)(abpt->wf * qlen); + int geometry = (abpt->wb < 0 || bw > qlen) ? 0 : (bw >= 1 ? 1 : 2); + abpoa_score_range_t r16 = simd_abpoa_score_range(abpt, 16, simd_p16.num_of_value, qlen, gn, geometry); + abpoa_score_range_t r32 = simd_abpoa_score_range(abpt, 32, simd_p32.num_of_value, qlen, gn, geometry == 0 ? 0 : 2), *r; +#ifdef ABPOA_TEST_SCORE_RANGE + { // tests only (tests/test_score_range.sh): force 32 bits, or 16 bits at a floor that can lose + char *e = getenv("ABPOA_TEST_BITS"); + if (e && atoi(e) == 32) r16.feasible = 0; + if ((e = getenv("ABPOA_TEST_FLOOR16")) && r16.fits) { + int64_t d = MIN_OF_TWO(atoi(e) - r16.floor, -1 - r16.ztop); // the checks need ztop < 0 + r16.floor += d; r16.fmax += d; r16.ztop += d; + r16.feasible = 1; r16.proven = 0; + } + } +#endif + bits = r16.feasible ? 16 : 32; + if (abpt->verbose >= ABPOA_DEBUG_VERBOSE) fprintf(stderr, "score width: %d bits%s\n", bits, r16.feasible && !r16.proven ? " (checked)" : ""); + // what a failed 16-bit attempt changes, restored before it is redone at 32 bits + abpoa_res_t res0 = *res; int *mpos_log = NULL, mpos_n = 0, mpos_m = 0, mpos_log_on, aln_fail, chk_zone; + int64_t chk_fmax, chk_ztop, chk_u, chk_hi; +retry: + r = bits == 16 ? &r16 : &r32; + if (!r->fits) + err_fatal(__func__, "scores too large for 32-bit alignment (qlen=%d, graph rows=%d, match<=%d, mismatch>=%d, " + "gap open=%d,%d, gap extension=%d,%d): scale them down.", qlen, gn, abpt->max_mat, -abpt->min_mis, + abpt->gap_open1, abpt->gap_open2, gap_ext1, gap_ext2); + aln_fail = 0; chk_zone = !r->proven; mpos_log_on = bits == 16; mpos_n = 0; + chk_fmax = r->fmax, chk_ztop = r->ztop, chk_u = r->u, chk_hi = bits == 16 ? INT16_MAX : INT32_MAX; + if (bits == 16) { + simd_p16.inf_min = (int)r->floor; mem_ret = simd_abpoa_realloc(ab, gn, qlen, abpt, simd_p16); - bits = 16; } else { - simd_p32.inf_min = MAX_OF_THREE(INT32_MIN + abpt->min_mis, INT32_MIN + gap_oe1, INT32_MIN + gap_oe2) + 512 * MAX_OF_TWO(gap_ext1, gap_ext2); + simd_p32.inf_min = (int)r->floor; mem_ret = simd_abpoa_realloc(ab, gn, qlen, abpt, simd_p32); - bits = 32; } - if (mem_ret) return 0; + if (mem_ret) { free(mpos_log); return 0; } if (bits == 16) { if (abpt->gap_mode == ABPOA_LINEAR_GAP) { @@ -1343,6 +1494,26 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b SIMDAddi32, SIMDSubi32, SIMDShiftOneNi32, SIMDSetIfGreateri32, SIMDGetIfGreateri32, SIMDSetIfEquali32); } } + if (aln_fail) { + if (bits == 32 && aln_fail == ABPOA_FAIL_BT && r->proven) // the scores provably fit: not an overflow + err_fatal(__func__, "Error in backtrack (32-bit DP, scores in range): qlen=%d, graph rows=%d.", qlen, gn); + if (bits == 32) + err_fatal(__func__, "alignment failed its score checks at 32 bits (%s%s%s): qlen=%d, graph rows=%d, match<=%d, " + "mismatch>=%d, gap open=%d,%d, gap extension=%d,%d, floor=%lld, zone=(%lld,%lld], max=%lld. " + "Scores this large cannot be aligned exactly: scale them down.", + aln_fail & ABPOA_FAIL_BT ? "backtrack " : "", aln_fail & ABPOA_FAIL_ZONE ? "zone " : "", + aln_fail & ABPOA_FAIL_BEST ? "best" : "", qlen, gn, abpt->max_mat, -abpt->min_mis, abpt->gap_open1, + abpt->gap_open2, gap_ext1, gap_ext2, (long long)r->floor, (long long)r->fmax, (long long)r->ztop, (long long)r->u); + if (abpt->verbose >= ABPOA_DEBUG_VERBOSE) + fprintf(stderr, "16-bit alignment failed its score checks (%d), redone at 32 bits: qlen=%d, graph rows=%d\n", aln_fail, qlen, gn); + *res = res0; + for (i = mpos_n - 1; i >= 0; --i) { + ab->abg->node_id_to_max_pos_left[mpos_log[i*3]] = mpos_log[i*3+1]; + ab->abg->node_id_to_max_pos_right[mpos_log[i*3]] = mpos_log[i*3+2]; + } + bits = 32; goto retry; + } + free(mpos_log); #endif free(index_map); abpt->verbose = old_verbose; diff --git a/tests/run_all.sh b/tests/run_all.sh index 087f1a6..c2e9178 100755 --- a/tests/run_all.sh +++ b/tests/run_all.sh @@ -49,6 +49,21 @@ if [ -f bin/abpoa ]; then else echo "FAIL"; FAIL=$((FAIL + 1)); ERRORS="$ERRORS\n - heter regression" fi + + # a trained model with gap extension 60, which the 16-bit DP could not hold ("Error in + # cg_backtrack"): every score halved is the same model, so the output must not change + printf " %-45s " "large scores match the same scores halved" + ok=1 + for opts in "-r1" "-r1 -b -1" "-r1 -b 1000 -f 0.1 -p"; do + a=$(./bin/abpoa test_data/seq.fa $opts -M 210 -X 1020 -O 1290,3870 -E 60,2 2>/dev/null) || ok=0 + b=$(./bin/abpoa test_data/seq.fa $opts -M 105 -X 510 -O 645,1935 -E 30,1 2>/dev/null) || ok=0 + [ -n "$a" ] && [ "$a" = "$b" ] || ok=0 + done + if [ $ok = 1 ]; then + echo "PASS"; PASS=$((PASS + 1)) + else + echo "FAIL"; FAIL=$((FAIL + 1)); ERRORS="$ERRORS\n - large score regression" + fi echo "" fi From 731242172a6ca9067d301b9bb4fbcd53707dde94 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 14:00:28 -0400 Subject: [PATCH 2/9] Test the score width choice, its checks and the 32-bit retry 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- tests/test_score_range.sh | 128 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100755 tests/test_score_range.sh diff --git a/tests/test_score_range.sh b/tests/test_score_range.sh new file mode 100755 index 0000000..2676104 --- /dev/null +++ b/tests/test_score_range.sh @@ -0,0 +1,128 @@ +#!/bin/bash +# Score width (simd_abpoa_score_range in src/abpoa_align_simd.c): the choice between 16- and 32-bit +# DP, the checks on every 16-bit alignment that is not proven safe, and the retry at 32 bits. +# +# abPOA is built once more with -DABPOA_TEST_SCORE_RANGE, under which ABPOA_TEST_BITS=32 forces +# 32 bits and ABPOA_TEST_FLOOR16= forces 16 bits wherever the scores fit, with "minus infinity" +# at v, high enough to lose to real scores (capped at the highest floor the checks can cover). +# Every output must equal the plain build's: +# 1. 16 and 32 bits agree, in every gap and alignment mode, banded and not +# 2. a floor that loses is caught, and the alignment redone at 32 bits gives the same result +# 3. on both sides of the 16-bit limits (query length, graph depth) the result is the 32-bit one +# 4. scores that the previous 16-bit choice could not hold (a trained model with gap extension +# 60, which aborted with "Error in cg_backtrack") give the 32-bit result, and the same as with +# every score halved +set -e + +SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)" +PROJECT_DIR="$(cd "$SCRIPT_DIR/.." && pwd)" +cd "$PROJECT_DIR" + +BUILD_DIR="${TMPDIR:-/tmp}/abpoa_score_range_$$" +trap "rm -rf $BUILD_DIR" EXIT +ERRORS=0 +fail() { echo "FAIL: $*"; ERRORS=$((ERRORS + 1)); } + +# test builds: the default one, and on x86 also SSE2, whose vectors hold half as many scores +TEST_FLAGS="-Wall -Wno-unused-function -Wno-misleading-indentation -DUSE_SIMDE -DSIMDE_ENABLE_NATIVE_ALIASES -DABPOA_TEST_SCORE_RANGE" +BINS="" +build() { # name, make args... + local d="$BUILD_DIR/$1"; shift + mkdir -p "$d" && cp -r src Makefile "$d/" && ln -s "$PROJECT_DIR/include" "$d/include" + make -C "$d" clean >/dev/null 2>&1 + if make -C "$d" -j2 EXTRA_FLAGS="$TEST_FLAGS" "$@" >"$d/make.log" 2>&1; then + BINS="$BINS $(ls "$d"/bin/abpoa*)" + else + tail -20 "$d/make.log"; fail "test build ($*) failed" + fi +} +build default +case "$(uname -m)" in x86_64|amd64) build sse2 sse2=1 ;; esac +[ -n "$BINS" ] || { echo "$ERRORS score range errors"; exit 1; } + +# deterministic random sequences: randseq [, applied to the seed's own sequence] +randseq() { + awk -v seed="$1" -v n="$2" -v mu="${3:-0}" 'BEGIN { srand(seed); s = ""; + for (i = 0; i < n; ++i) s = s substr("ACGT", int(rand() * 4) + 1, 1); + if (mu > 0) { srand(seed + 7919); t = ""; + for (i = 1; i <= n; ++i) { r = rand(); c = substr(s, i, 1); + if (r < mu / 3) continue; else if (r < 2 * mu / 3) c = substr("ACGT", int(rand() * 4) + 1, 1); + else if (r < mu) c = c substr("ACGT", int(rand() * 4) + 1, 1); t = t c } + s = t } + print s }' +} + +# every way of running bin against in.fa with opts must give the plain run's output +same_as_plain() { # bin, in.fa, opts... + local bin="$1" fa="$2"; shift 2 + local ref out f + if ! ref=$("$bin" "$@" "$fa" 2>/dev/null); then fail "$(basename "$bin") $* $fa exits non-zero"; return; fi + out=$(ABPOA_TEST_BITS=32 "$bin" "$@" "$fa" 2>/dev/null) || true + [ "$out" = "$ref" ] || fail "$(basename "$bin") $* $fa: 32 bits differ from the plain run" + for f in -20000 -6000 -2000 -500 1 2000; do + out=$(ABPOA_TEST_FLOOR16=$f "$bin" "$@" "$fa" 2>/dev/null) || true + [ "$out" = "$ref" ] || fail "$(basename "$bin") $* $fa: 16 bits at floor $f differ from the plain run" + done +} + +# 1, 2: gap modes x alignment modes x banding, on the shipped test data +for bin in $BINS; do + for fa in test_data/seq.fa test_data/heter.fa; do + for gap in "" "-O 4,0 -E 2,0" "-O 0,0 -E 2,0"; do + for mode in 0 1 2; do + for band in "" "-b -1"; do + same_as_plain "$bin" "$fa" -r1 -m $mode $gap $band + done + done + done + for extra in "-s -r2" "-G -r1" "-p -r1"; do + same_as_plain "$bin" "$fa" $extra + done + done + # the losing floors must actually have been caught and retried + n=$(ABPOA_TEST_FLOOR16=2000 "$bin" -V 2 -r1 test_data/heter.fa 2>&1 >/dev/null | grep -c "redone at 32 bits" || true) + [ "$n" -gt 0 ] || fail "$(basename "$bin"): a 16-bit floor of 2000 was never caught" +done + +# 3: across the 16-bit limits, with cactus-like trained scores and band (every row covers the query) +T1="-M 210 -X 1080 -O 1410,4230 -E 30,1 -b 1000 -f 0.1 -r1" +widths() { "$@" -V 2 2>&1 >/dev/null | grep -o "score width: [0-9]* bits" | sort -u | tr '\n' ' '; } +for bin in $BINS; do + seen="" + for len in $(awk 'BEGIN { for (i = 125; i <= 145; ++i) print i }'); do # query length: 135 is the last 16-bit one + fa="$BUILD_DIR/q$len.fa" + printf ">a\n%s\n>b\n%s\n>c\n%s\n" "$(randseq $len $len)" "$(randseq $len $len 0.1)" "$(randseq $((len + 1)) $len 0.05)" > "$fa" + same_as_plain "$bin" "$fa" $T1 + seen="$seen$(widths "$bin" $T1 "$fa")" + done + for w in 16 32; do echo "$seen" | grep -q "$w bits" || fail "$(basename "$bin"): query length sweep never ran at $w bits"; done + seen="" + for len in $(awk 'BEGIN { for (i = 7300; i <= 8100; i += 50) print i }'); do # graph depth for a 130 bp query + fa="$BUILD_DIR/d$len.fa" + printf ">a\n%s\n>b\n%s\n" "$(randseq $len $len)" "$(randseq 1 130)" > "$fa" + same_as_plain "$bin" "$fa" $T1 + seen="$seen$(widths "$bin" $T1 "$fa")" + done + for w in 16 32; do echo "$seen" | grep -q "$w bits" || fail "$(basename "$bin"): graph depth sweep never ran at $w bits"; done +done + +# 4: the HPRC oct9 chr2 model, scaled 30x, and the same model halved +T2="-M 210 -X 1020 -O 1290,3870 -E 60,2" +T2H="-M 105 -X 510 -O 645,1935 -E 30,1" +for bin in $BINS; do + for fa in test_data/seq.fa $BUILD_DIR/q125.fa $BUILD_DIR/q140.fa; do + for opts in "-r1" "-r1 -b -1" "-r1 -b 1000 -f 0.1 -p" "-r1 -m 1" "-r1 -m 2 -b -1"; do + same_as_plain "$bin" "$fa" $T2 $opts + a=$("$bin" $T2 $opts "$fa" 2>/dev/null) || true + b=$("$bin" $T2H $opts "$fa" 2>/dev/null) || true + [ -n "$a" ] && [ "$a" = "$b" ] || fail "$(basename "$bin") $opts $fa: chr2 model differs from it halved" + done + done +done + +if [ $ERRORS -gt 0 ]; then + echo "$ERRORS score range errors" + exit 1 +fi +echo "Score range checks passed" +exit 0 From 8eeb08929254b9c43e164df51be9078e2f774719 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 15:57:11 -0400 Subject: [PATCH 3/9] Scan for floor-zone scores after the DP, out of line 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- src/abpoa_align_simd.c | 67 ++++++++++++++++++++++-------------------- 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/src/abpoa_align_simd.c b/src/abpoa_align_simd.c index e8bcfce..c0a6e1c 100644 --- a/src/abpoa_align_simd.c +++ b/src/abpoa_align_simd.c @@ -481,8 +481,6 @@ int w = abpt->wb < 0 ? qlen : abpt->wb+(int)(abpt->wf*qlen); /* when w < 0, do whole global */ \ int best_i = 0, best_j = 0, best_id = 0, max, left_max_i=-1, right_max_i=-1; \ SIMDi zero = SIMDSetZeroi(), SIMD_INF_MIN = SIMDSetOne(inf_min); \ - SIMDi ZONE_FMAX = SIMDSetOne((score_t)chk_fmax), ZONE_MIN = SIMDSetOne((score_t)chk_hi); \ - SIMDi H_MAX = SIMD_INF_MIN; /* simd_abpoa_zone_scan */ \ pn = sp.num_of_value; qp_sn = dp_sn = (matrix_col_n + pn - 1) / pn; \ log_n = sp.log_num, size = sp.size; qp = abm->s_mem; \ int set_num; SIMDi *PRE_MASK, *SUF_MIN, *PRE_MIN; \ @@ -1179,27 +1177,46 @@ } \ } -/* only when the width is not proven safe (see simd_abpoa_score_range): every stored H of the row, - * min over those above fmax and max over all, tested once per alignment by simd_abpoa_check_scores */ -#define simd_abpoa_zone_scan(first_sn, last_sn, SIMDSetIfGreater, SIMDMax) { \ - SIMDi _t; \ - for (sn_i = first_sn; sn_i <= last_sn; ++sn_i) { \ - SIMDSetIfGreater(_t, dp_h[sn_i], ZONE_FMAX, dp_h[sn_i], ZONE_MIN); \ - SIMDSetIfGreater(ZONE_MIN, ZONE_MIN, _t, _t, ZONE_MIN); \ - H_MAX = SIMDMax(H_MAX, dp_h[sn_i]); \ - } \ +#define ABPOA_FAIL_BT 1 /* the backtrack could not retrace the DP */ +#define ABPOA_FAIL_ZONE 2 /* a stored H in (fmax, ztop] or above u */ +#define ABPOA_FAIL_BEST 4 /* the best score in (-inf, ztop] or above u */ + +typedef struct { int64_t floor, fmax, ztop, u; int fits, feasible, proven; } abpoa_score_range_t; /* simd_abpoa_score_range */ + +/* only when the width is not proven safe (see simd_abpoa_score_range), once the DP is done: every stored H of + * every row computed (row 0 included, none after an extend-mode z-drop), for a value in (fmax, ztop] or above u. + * Out of line: inside the DP function, the scan's code alone made the compiler build a slower DP loop */ +#define simd_abpoa_zone_scan_def(name, score_t, SIMDSetOne, SIMDSetIfGreater, SIMDMax) \ +static ATTRIBUTE((noinline)) int name(SIMDi **dp_row, int *dp_beg_sn, int *dp_end_sn, uint8_t *row_map, \ + int last_row, abpoa_score_range_t *r) { \ + union { SIMDi v; score_t s[sizeof(SIMDi) / sizeof(score_t)]; } zmin, hmax; /* lanes read at the end */ \ + int row, sn_i, i, fail = 0; \ + SIMDi FMAX = SIMDSetOne((score_t)r->fmax), t, *dp_h; \ + zmin.v = SIMDSetOne((score_t)(r->ztop + 1)); hmax.v = FMAX; \ + for (row = 0; row <= last_row; ++row) { \ + if (row_map[row] == 0) continue; /* not part of the subgraph: never computed */ \ + dp_h = dp_row[row]; \ + for (sn_i = dp_beg_sn[row]; sn_i <= dp_end_sn[row]; ++sn_i) { \ + SIMDSetIfGreater(t, dp_h[sn_i], FMAX, dp_h[sn_i], zmin.v); /* t = H > fmax ? H : ztop+1 */ \ + SIMDSetIfGreater(zmin.v, zmin.v, t, t, zmin.v); /* min(zmin, t) */ \ + hmax.v = SIMDMax(hmax.v, dp_h[sn_i]); \ + } \ + } \ + for (i = 0; i < (int)(sizeof(SIMDi) / sizeof(score_t)); ++i) { \ + if (zmin.s[i] <= r->ztop || hmax.s[i] > r->u) fail = ABPOA_FAIL_ZONE; \ + } \ + return fail; \ } +simd_abpoa_zone_scan_def(simd_abpoa_zone_scan16, int16_t, SIMDSetOnei16, SIMDSetIfGreateri16, SIMDMaxi16) +simd_abpoa_zone_scan_def(simd_abpoa_zone_scan32, int32_t, SIMDSetOnei32, SIMDSetIfGreateri32, SIMDMaxi32) /* a real score in (fmax, ztop] can lose to, or be confused with, a floor-derived one, and a score above * u can only be a wrap: either way the alignment is redone at 32 bits (simd_abpoa_align_sequence_to_subgraph) */ #define simd_abpoa_check_scores(score_t) { \ - if (chk_zone) { \ - score_t *_zmin = (score_t*)&ZONE_MIN, *_hmax = (score_t*)&H_MAX; \ - for (i = 0; i < pn; ++i) { \ - if (_zmin[i] <= chk_ztop || _hmax[i] > chk_u) aln_fail |= ABPOA_FAIL_ZONE; \ - } \ - } \ - if ((best_i > 0 || abpt->align_mode == ABPOA_GLOBAL_MODE) && (best_score > chk_u || best_score <= chk_ztop)) \ + if (chk_zone) /* rows 0 to the last computed: the loop ran to its end, or stopped at a z-drop */ \ + aln_fail |= (sizeof(score_t) == 2 ? simd_abpoa_zone_scan16 : simd_abpoa_zone_scan32)(dp_row, dp_beg_sn, \ + dp_end_sn, index_map + beg_index, index_i < end_index ? dp_i : dp_i - 1, r); \ + if ((best_i > 0 || abpt->align_mode == ABPOA_GLOBAL_MODE) && (best_score > r->u || best_score <= r->ztop)) \ aln_fail |= ABPOA_FAIL_BEST; \ } @@ -1209,11 +1226,9 @@ SIMDSub, SIMDShiftOneN, SIMDSetIfGreater, SIMDGetIfGreater) { \ simd_abpoa_lg_var(score_t, sp, SIMDSetOne, SIMDShiftOneN, SIMDAdd); \ simd_abpoa_lg_first_dp(score_t); \ - if (chk_zone) simd_abpoa_zone_scan(0, _end_sn, SIMDSetIfGreater, SIMDMax); \ for (index_i=beg_index+1, dp_i=1; index_ialign_mode == ABPOA_LOCAL_MODE) { \ simd_abpoa_max_in_row(score_t, SIMDSetIfGreater, SIMDGetIfGreater); \ set_global_max_score(max, dp_i, left_max_i); \ @@ -1246,11 +1261,9 @@ SIMDSub, SIMDShiftOneN, SIMDSetIfGreater, SIMDGetIfGreater, SIMDSetIfEqual) { \ simd_abpoa_ag_var(score_t, sp, SIMDSetOne, SIMDShiftOneN, SIMDAdd); \ simd_abpoa_ag_first_dp(score_t); \ - if (chk_zone) simd_abpoa_zone_scan(0, _end_sn, SIMDSetIfGreater, SIMDMax); \ for (index_i=beg_index+1, dp_i=1; index_ialign_mode == ABPOA_LOCAL_MODE) { \ simd_abpoa_max_in_row(score_t, SIMDSetIfGreater, SIMDGetIfGreater); \ set_global_max_score(max, dp_i, left_max_i); \ @@ -1282,11 +1295,9 @@ SIMDSub, SIMDShiftOneN, SIMDSetIfGreater, SIMDGetIfGreater, SIMDSetIfEqual) { \ simd_abpoa_cg_var(score_t, sp, SIMDSetOne, SIMDShiftOneN, SIMDAdd); \ simd_abpoa_cg_first_dp(score_t); \ - if (chk_zone) simd_abpoa_zone_scan(0, _end_sn, SIMDSetIfGreater, SIMDMax); \ for (index_i=beg_index+1, dp_i=1; index_ialign_mode == ABPOA_LOCAL_MODE) { \ simd_abpoa_max_in_row(score_t, SIMDSetIfGreater, SIMDGetIfGreater); \ set_global_max_score(max, dp_i, left_max_i); \ @@ -1330,12 +1341,6 @@ * M/X: largest match score/mismatch penalty, OE/EM: largest gap open+extension/extension, PS: largest path * score penalty (20, or 0 without inc_path_score), D: rows on a path (gn bounds it), g(n): cheapest gap of n. */ -#define ABPOA_FAIL_BT 1 /* the backtrack could not retrace the DP */ -#define ABPOA_FAIL_ZONE 2 /* a stored H in (fmax, ztop] or above u */ -#define ABPOA_FAIL_BEST 4 /* the best score in (-inf, ztop] or above u */ - -typedef struct { int64_t floor, fmax, ztop, u; int fits, feasible, proven; } abpoa_score_range_t; - static int64_t simd_abpoa_gap_cost(abpoa_para_t *abpt, int64_t n) { if (n <= 0) return 0; int64_t c = (abpt->gap_mode == ABPOA_LINEAR_GAP ? 0 : abpt->gap_open1) + (int64_t)abpt->gap_ext1 * n; @@ -1453,7 +1458,6 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b if (abpt->verbose >= ABPOA_DEBUG_VERBOSE) fprintf(stderr, "score width: %d bits%s\n", bits, r16.feasible && !r16.proven ? " (checked)" : ""); // what a failed 16-bit attempt changes, restored before it is redone at 32 bits abpoa_res_t res0 = *res; int *mpos_log = NULL, mpos_n = 0, mpos_m = 0, mpos_log_on, aln_fail, chk_zone; - int64_t chk_fmax, chk_ztop, chk_u, chk_hi; retry: r = bits == 16 ? &r16 : &r32; if (!r->fits) @@ -1461,7 +1465,6 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b "gap open=%d,%d, gap extension=%d,%d): scale them down.", qlen, gn, abpt->max_mat, -abpt->min_mis, abpt->gap_open1, abpt->gap_open2, gap_ext1, gap_ext2); aln_fail = 0; chk_zone = !r->proven; mpos_log_on = bits == 16; mpos_n = 0; - chk_fmax = r->fmax, chk_ztop = r->ztop, chk_u = r->u, chk_hi = bits == 16 ? INT16_MAX : INT32_MAX; if (bits == 16) { simd_p16.inf_min = (int)r->floor; mem_ret = simd_abpoa_realloc(ab, gn, qlen, abpt, simd_p16); From 0f6d1559a9f5e563f321a321edab52678b43361d Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 16:06:13 -0400 Subject: [PATCH 4/9] Clamp row 0's gap scores at the floor, in 64 bits 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- src/abpoa_align_simd.c | 13 ++++++++----- tests/test_score_range.sh | 15 +++++++++++++++ 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/src/abpoa_align_simd.c b/src/abpoa_align_simd.c index c0a6e1c..6ca3bfc 100644 --- a/src/abpoa_align_simd.c +++ b/src/abpoa_align_simd.c @@ -664,6 +664,9 @@ _end_sn = MIN_OF_TWO(dp_end_sn[0]+1, dp_sn-1); \ } +/* row 0, outside local mode: a gap along the query, each gap piece's score taken in 64 bits and clamped at + * inf_min. The width choice bounds only the cheaper piece, and the other one's -(o + e*i) would wrap; clamped, + * it still loses to the cheaper one, and a real score that low fails the checks (simd_abpoa_check_scores) */ #define simd_abpoa_lg_first_dp(score_t) { \ simd_abpoa_lg_first_row; \ if (abpt->align_mode == ABPOA_LOCAL_MODE) { \ @@ -675,7 +678,7 @@ } \ _dp_h = (score_t*)dp_h; \ for (i = 0; i <= dp_end[0]; ++i) { /* no SIMD parallelization */ \ - _dp_h[i] = -gap_ext1 * i; \ + _dp_h[i] = MAX_OF_TWO(-(int64_t)gap_ext1 * i, inf_min); \ } \ } \ } @@ -692,8 +695,8 @@ _dp_h=(score_t*)dp_h,_dp_e1=(score_t*)dp_e1,_dp_f1=(score_t*)dp_f1; \ _dp_h[0] = 0; _dp_e1[0] = -(gap_oe1), _dp_f1[0] = inf_min; \ for (i = 1; i <= dp_end[0]; ++i) { /* no SIMD parallelization */ \ - _dp_f1[i] = -gap_open1 - gap_ext1 * i; \ - _dp_h[i] = -gap_open1 - gap_ext1 * i; \ + _dp_f1[i] = MAX_OF_TWO(-gap_open1 - (int64_t)gap_ext1 * i, inf_min); \ + _dp_h[i] = _dp_f1[i]; \ } \ } \ } @@ -712,8 +715,8 @@ _dp_h[0] = 0; _dp_e1[0] = -(gap_oe1); _dp_e2[0] = -(gap_oe2); \ _dp_f1[0] = _dp_f2[0] = inf_min; \ for (i = 1; i <= dp_end[0]; ++i) { /* no SIMD parallelization */ \ - _dp_f1[i] = -gap_open1 - gap_ext1 * i; \ - _dp_f2[i] = -gap_open2 - gap_ext2 * i; \ + _dp_f1[i] = MAX_OF_TWO(-gap_open1 - (int64_t)gap_ext1 * i, inf_min); \ + _dp_f2[i] = MAX_OF_TWO(-gap_open2 - (int64_t)gap_ext2 * i, inf_min); \ _dp_h[i] = MAX_OF_TWO(_dp_f1[i], _dp_f2[i]); \ } \ } \ diff --git a/tests/test_score_range.sh b/tests/test_score_range.sh index 2676104..15316a8 100755 --- a/tests/test_score_range.sh +++ b/tests/test_score_range.sh @@ -12,6 +12,8 @@ # 4. scores that the previous 16-bit choice could not hold (a trained model with gap extension # 60, which aborted with "Error in cg_backtrack") give the 32-bit result, and the same as with # every score halved +# 5. a 16-bit alignment proven safe needs no retry, even where row 0's gap on the dearer piece +# leaves the 16-bit range set -e SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)" @@ -120,6 +122,19 @@ for bin in $BINS; do done done +# 5: convex gaps with a gap extension above the match score: the width choice bounds the cheaper +# piece's gap along row 0, -(24 + i), while the other piece's -(4 + 100*i) is out of 16-bit range +# past column 327; that must neither change the result nor make a proven alignment need its retry +fa="$BUILD_DIR/row0.fa" +printf ">a\n%s\n>b\n%s\n" "$(randseq 5 500)" "$(randseq 6 400)" > "$fa" +R0="-M 1 -X 4 -O 4,24 -E 100,1 -b -1 -r1" +for bin in $BINS; do + same_as_plain "$bin" "$fa" $R0 + [ "$(widths "$bin" $R0 "$fa")" = "score width: 16 bits " ] || fail "$(basename "$bin") $R0: not a 16-bit alignment" + n=$("$bin" $R0 -V 2 "$fa" 2>&1 >/dev/null | grep -c "redone at 32 bits" || true) + [ "$n" = 0 ] || fail "$(basename "$bin") $R0: a proven 16-bit alignment was redone at 32 bits" +done + if [ $ERRORS -gt 0 ]; then echo "$ERRORS score range errors" exit 1 From 8d41b3073f995d5dbeca12f2be26e3f790e55ffd Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 16:13:05 -0400 Subject: [PATCH 5/9] Print the match and mismatch scores the width rule used 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- src/abpoa_align_simd.c | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/abpoa_align_simd.c b/src/abpoa_align_simd.c index 6ca3bfc..cf1f9b4 100644 --- a/src/abpoa_align_simd.c +++ b/src/abpoa_align_simd.c @@ -1184,7 +1184,7 @@ #define ABPOA_FAIL_ZONE 2 /* a stored H in (fmax, ztop] or above u */ #define ABPOA_FAIL_BEST 4 /* the best score in (-inf, ztop] or above u */ -typedef struct { int64_t floor, fmax, ztop, u; int fits, feasible, proven; } abpoa_score_range_t; /* simd_abpoa_score_range */ +typedef struct { int64_t M, X, floor, fmax, ztop, u; int fits, feasible, proven; } abpoa_score_range_t; /* simd_abpoa_score_range */ /* only when the width is not proven safe (see simd_abpoa_score_range), once the DP is done: every stored H of * every row computed (row 0 included, none after an extend-mode z-drop), for a value in (fmax, ztop] or above u. @@ -1368,7 +1368,7 @@ static abpoa_score_range_t simd_abpoa_score_range(abpoa_para_t *abpt, int bits, int64_t chain = geometry == 0 ? X + PS : geometry == 1 ? (pn + 1) * (X + PS) : (Q + 2) * (X + PS); int64_t growth = geometry == 0 ? 0 : geometry == 1 ? pn * M : (Q + 1) * M; int64_t step = MAX_OF_THREE(X + PS, OE + PS, pn / 2 * EM); // largest single subtraction - r.u = Q * M; + r.M = M; r.X = X; r.u = Q * M; r.floor = lo + chain + OE + pn / 2 * EM + PS; r.fmax = r.floor + growth; r.ztop = r.fmax + M + 2 * X + 2 * OE + 3 * PS; @@ -1464,8 +1464,8 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b retry: r = bits == 16 ? &r16 : &r32; if (!r->fits) - err_fatal(__func__, "scores too large for 32-bit alignment (qlen=%d, graph rows=%d, match<=%d, mismatch>=%d, " - "gap open=%d,%d, gap extension=%d,%d): scale them down.", qlen, gn, abpt->max_mat, -abpt->min_mis, + err_fatal(__func__, "scores too large for 32-bit alignment (qlen=%d, graph rows=%d, match<=%lld, mismatch>=%lld, " + "gap open=%d,%d, gap extension=%d,%d): scale them down.", qlen, gn, (long long)r->M, -(long long)r->X, abpt->gap_open1, abpt->gap_open2, gap_ext1, gap_ext2); aln_fail = 0; chk_zone = !r->proven; mpos_log_on = bits == 16; mpos_n = 0; if (bits == 16) { @@ -1504,11 +1504,11 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b if (bits == 32 && aln_fail == ABPOA_FAIL_BT && r->proven) // the scores provably fit: not an overflow err_fatal(__func__, "Error in backtrack (32-bit DP, scores in range): qlen=%d, graph rows=%d.", qlen, gn); if (bits == 32) - err_fatal(__func__, "alignment failed its score checks at 32 bits (%s%s%s): qlen=%d, graph rows=%d, match<=%d, " - "mismatch>=%d, gap open=%d,%d, gap extension=%d,%d, floor=%lld, zone=(%lld,%lld], max=%lld. " + err_fatal(__func__, "alignment failed its score checks at 32 bits (%s%s%s): qlen=%d, graph rows=%d, match<=%lld, " + "mismatch>=%lld, gap open=%d,%d, gap extension=%d,%d, floor=%lld, zone=(%lld,%lld], max=%lld. " "Scores this large cannot be aligned exactly: scale them down.", aln_fail & ABPOA_FAIL_BT ? "backtrack " : "", aln_fail & ABPOA_FAIL_ZONE ? "zone " : "", - aln_fail & ABPOA_FAIL_BEST ? "best" : "", qlen, gn, abpt->max_mat, -abpt->min_mis, abpt->gap_open1, + aln_fail & ABPOA_FAIL_BEST ? "best" : "", qlen, gn, (long long)r->M, -(long long)r->X, abpt->gap_open1, abpt->gap_open2, gap_ext1, gap_ext2, (long long)r->floor, (long long)r->fmax, (long long)r->ztop, (long long)r->u); if (abpt->verbose >= ABPOA_DEBUG_VERBOSE) fprintf(stderr, "16-bit alignment failed its score checks (%d), redone at 32 bits: qlen=%d, graph rows=%d\n", aln_fail, qlen, gn); From 1a682ec73c7c99a322394654b31c01a8823d8ce0 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 16:41:02 -0400 Subject: [PATCH 6/9] Align scores the 32-bit bounds cannot prove, checked, instead of refusing 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- src/abpoa_align_simd.c | 76 +++++++++++++++++++++++++-------------- tests/test_score_range.sh | 19 ++++++++++ 2 files changed, 68 insertions(+), 27 deletions(-) diff --git a/src/abpoa_align_simd.c b/src/abpoa_align_simd.c index cf1f9b4..2de328f 100644 --- a/src/abpoa_align_simd.c +++ b/src/abpoa_align_simd.c @@ -1180,46 +1180,50 @@ } \ } -#define ABPOA_FAIL_BT 1 /* the backtrack could not retrace the DP */ -#define ABPOA_FAIL_ZONE 2 /* a stored H in (fmax, ztop] or above u */ -#define ABPOA_FAIL_BEST 4 /* the best score in (-inf, ztop] or above u */ +#define ABPOA_FAIL_BT 1 /* the backtrack could not retrace the DP */ +#define ABPOA_FAIL_ZONE 2 /* a stored H in (fmax, ztop] */ +#define ABPOA_FAIL_BEST 4 /* the best score in (-inf, ztop] or above top */ +#define ABPOA_FAIL_RANGE 8 /* a stored H below low or above top: a later step could wrap */ -typedef struct { int64_t M, X, floor, fmax, ztop, u; int fits, feasible, proven; } abpoa_score_range_t; /* simd_abpoa_score_range */ +typedef struct { int64_t M, X, low, floor, fmax, ztop, top; int fits, feasible, proven; } abpoa_score_range_t; /* simd_abpoa_score_range */ /* only when the width is not proven safe (see simd_abpoa_score_range), once the DP is done: every stored H of - * every row computed (row 0 included, none after an extend-mode z-drop), for a value in (fmax, ztop] or above u. - * Out of line: inside the DP function, the scan's code alone made the compiler build a slower DP loop */ + * every row computed (row 0 included, none after an extend-mode z-drop), for a value in (fmax, ztop], below low + * or above top. Out of line: inside the DP function, the scan's code alone made the compiler build a slower DP loop */ #define simd_abpoa_zone_scan_def(name, score_t, SIMDSetOne, SIMDSetIfGreater, SIMDMax) \ static ATTRIBUTE((noinline)) int name(SIMDi **dp_row, int *dp_beg_sn, int *dp_end_sn, uint8_t *row_map, \ int last_row, abpoa_score_range_t *r) { \ - union { SIMDi v; score_t s[sizeof(SIMDi) / sizeof(score_t)]; } zmin, hmax; /* lanes read at the end */ \ + union { SIMDi v; score_t s[sizeof(SIMDi) / sizeof(score_t)]; } zmin, hmin, hmax; /* lanes read at the end */ \ int row, sn_i, i, fail = 0; \ SIMDi FMAX = SIMDSetOne((score_t)r->fmax), t, *dp_h; \ - zmin.v = SIMDSetOne((score_t)(r->ztop + 1)); hmax.v = FMAX; \ + zmin.v = SIMDSetOne((score_t)(r->ztop + 1)); hmin.v = hmax.v = FMAX; \ for (row = 0; row <= last_row; ++row) { \ if (row_map[row] == 0) continue; /* not part of the subgraph: never computed */ \ dp_h = dp_row[row]; \ for (sn_i = dp_beg_sn[row]; sn_i <= dp_end_sn[row]; ++sn_i) { \ SIMDSetIfGreater(t, dp_h[sn_i], FMAX, dp_h[sn_i], zmin.v); /* t = H > fmax ? H : ztop+1 */ \ SIMDSetIfGreater(zmin.v, zmin.v, t, t, zmin.v); /* min(zmin, t) */ \ + SIMDSetIfGreater(hmin.v, hmin.v, dp_h[sn_i], dp_h[sn_i], hmin.v); \ hmax.v = SIMDMax(hmax.v, dp_h[sn_i]); \ } \ } \ for (i = 0; i < (int)(sizeof(SIMDi) / sizeof(score_t)); ++i) { \ - if (zmin.s[i] <= r->ztop || hmax.s[i] > r->u) fail = ABPOA_FAIL_ZONE; \ + if (zmin.s[i] <= r->ztop) fail |= ABPOA_FAIL_ZONE; \ + if (hmin.s[i] < r->low || hmax.s[i] > r->top) fail |= ABPOA_FAIL_RANGE; \ } \ return fail; \ } simd_abpoa_zone_scan_def(simd_abpoa_zone_scan16, int16_t, SIMDSetOnei16, SIMDSetIfGreateri16, SIMDMaxi16) simd_abpoa_zone_scan_def(simd_abpoa_zone_scan32, int32_t, SIMDSetOnei32, SIMDSetIfGreateri32, SIMDMaxi32) -/* a real score in (fmax, ztop] can lose to, or be confused with, a floor-derived one, and a score above - * u can only be a wrap: either way the alignment is redone at 32 bits (simd_abpoa_align_sequence_to_subgraph) */ +/* a real score in (fmax, ztop] can lose to, or be confused with, a floor-derived one, and from a score below low + * or above top a later step could wrap: either way a 16-bit alignment is redone at 32 bits, and a 32-bit one + * stops abPOA (simd_abpoa_align_sequence_to_subgraph) */ #define simd_abpoa_check_scores(score_t) { \ if (chk_zone) /* rows 0 to the last computed: the loop ran to its end, or stopped at a z-drop */ \ aln_fail |= (sizeof(score_t) == 2 ? simd_abpoa_zone_scan16 : simd_abpoa_zone_scan32)(dp_row, dp_beg_sn, \ dp_end_sn, index_map + beg_index, index_i < end_index ? dp_i : dp_i - 1, r); \ - if ((best_i > 0 || abpt->align_mode == ABPOA_GLOBAL_MODE) && (best_score > r->u || best_score <= r->ztop)) \ + if ((best_i > 0 || abpt->align_mode == ABPOA_GLOBAL_MODE) && (best_score > r->top || best_score <= r->ztop))\ aln_fail |= ABPOA_FAIL_BEST; \ } @@ -1339,8 +1343,14 @@ simd_abpoa_zone_scan_def(simd_abpoa_zone_scan32, int32_t, SIMDSetOnei32, SIMDSet * -(g(D) + g(qlen) + PS*D) (delete down to the row, insert across, each gap on its cheaper piece): this * "proves" a width. For narrower bands that value is only an estimate (a rigorous bound is far lower), so a * width the estimate allows is "feasible" and is used with every stored H checked: a real H in (fmax, ztop] - * is the first sign of the floor winning, and a value above u is a wrap, so a run that passes is exact and - * one that fails is redone at 32 bits. 32 bits takes bounds that need no band geometry. + * is the first sign of the floor winning, and an H below low (the floor for a full band) or above top (u) is + * where a later step could wrap, so a run that passes is exact and one that fails is redone at 32 bits. + * 32 bits takes bounds that need no band geometry. A 32-bit run they do not prove is checked instead, with the + * floor as low as the band's first SIMD block allows (floor-derived values outside it fail the checks, they do not + * go unseen), which leaves real scores the most room; where u does not fit either, top = min(u, MAX - M), from + * which no step can wrap. Such a run stops abPOA only if a check fails. Only scores whose zone does not fit + * below 0 at all, 3X + M + 3OE + pn/2*EM + 5PS >= 2^31 in 64 bits (or one of them beyond 32 bits), are refused + * before aligning: 32 bits cannot hold the range the checks need, so no result could be verified. * M/X: largest match score/mismatch penalty, OE/EM: largest gap open+extension/extension, PS: largest path * score penalty (20, or 0 without inc_path_score), D: rows on a path (gn bounds it), g(n): cheapest gap of n. */ @@ -1365,20 +1375,31 @@ static abpoa_score_range_t simd_abpoa_score_range(abpoa_para_t *abpt, int bits, int64_t oe2 = (int64_t)abpt->gap_open2 + abpt->gap_ext2; OE = MAX_OF_TWO(OE, oe2); OEmin = MIN_OF_TWO(OEmin, oe2); EM = MAX_OF_TWO(EM, abpt->gap_ext2); } - int64_t chain = geometry == 0 ? X + PS : geometry == 1 ? (pn + 1) * (X + PS) : (Q + 2) * (X + PS); - int64_t growth = geometry == 0 ? 0 : geometry == 1 ? pn * M : (Q + 1) * M; + // n: cells on a diagonal of floor-derived values: none, the band's first SIMD block, or any column + int64_t n = geometry == 0 ? 0 : geometry == 1 ? pn : Q + 1, zone = M + 2 * X + 2 * OE + 3 * PS; int64_t step = MAX_OF_THREE(X + PS, OE + PS, pn / 2 * EM); // largest single subtraction - r.M = M; r.X = X; r.u = Q * M; - r.floor = lo + chain + OE + pn / 2 * EM + PS; - r.fmax = r.floor + growth; - r.ztop = r.fmax + M + 2 * X + 2 * OE + 3 * PS; + r.M = M; r.X = X; r.top = Q * M; + r.low = lo + X + 2 * PS + OE + pn / 2 * EM; + r.floor = r.low + n * (X + PS); + r.fmax = r.floor + n * M; + r.ztop = r.fmax + zone; int64_t low_est = -(simd_abpoa_gap_cost(abpt, D) + simd_abpoa_gap_cost(abpt, Q) + PS * D); int64_t low = geometry == 0 ? low_est : -(MAX_OF_TWO(X, OEmin) * (Q + D) + PS * D); if (abpt->align_mode == ABPOA_LOCAL_MODE) low_est = low = 0; // ztop < 0: the checks can only miss the floor winning if the source's score of 0 is in or under the zone - r.fits = pn / 2 * EM <= hi && OE <= hi && M <= hi && X <= hi && Q < hi && r.ztop < 0 && r.u + step <= hi; + int consts = pn / 2 * EM <= hi && OE <= hi && M <= hi && X <= hi && Q < hi; + r.fits = consts && r.ztop < 0 && r.top + step <= hi; r.feasible = r.fits && low_est > r.ztop; r.proven = r.feasible && low > r.ztop; + if (bits == 32 && !r.proven && consts && r.low + zone < 0) { // checked: the floor as low as the band allows + n = MIN_OF_TWO(n, pn); + if (r.low + n * (X + PS + M) + zone >= 0) n = (-1 - r.low - zone) / (X + PS + M); + r.floor = r.low + n * (X + PS); + r.fmax = r.floor + n * M; + r.ztop = r.fmax + zone; + if (r.top + step > hi) r.top = MIN_OF_TWO(r.top, hi - M); + r.fits = 1; r.feasible = low_est > r.ztop; + } return r; } @@ -1451,14 +1472,14 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b char *e = getenv("ABPOA_TEST_BITS"); if (e && atoi(e) == 32) r16.feasible = 0; if ((e = getenv("ABPOA_TEST_FLOOR16")) && r16.fits) { - int64_t d = MIN_OF_TWO(atoi(e) - r16.floor, -1 - r16.ztop); // the checks need ztop < 0 + int64_t d = MIN_OF_TWO(MAX_OF_TWO(atoi(e), r16.low) - r16.floor, -1 - r16.ztop); // checks need ztop < 0, floor >= low r16.floor += d; r16.fmax += d; r16.ztop += d; r16.feasible = 1; r16.proven = 0; } } #endif bits = r16.feasible ? 16 : 32; - if (abpt->verbose >= ABPOA_DEBUG_VERBOSE) fprintf(stderr, "score width: %d bits%s\n", bits, r16.feasible && !r16.proven ? " (checked)" : ""); + if (abpt->verbose >= ABPOA_DEBUG_VERBOSE) fprintf(stderr, "score width: %d bits%s\n", bits, (bits == 16 ? r16 : r32).proven ? "" : " (checked)"); // what a failed 16-bit attempt changes, restored before it is redone at 32 bits abpoa_res_t res0 = *res; int *mpos_log = NULL, mpos_n = 0, mpos_m = 0, mpos_log_on, aln_fail, chk_zone; retry: @@ -1504,12 +1525,13 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b if (bits == 32 && aln_fail == ABPOA_FAIL_BT && r->proven) // the scores provably fit: not an overflow err_fatal(__func__, "Error in backtrack (32-bit DP, scores in range): qlen=%d, graph rows=%d.", qlen, gn); if (bits == 32) - err_fatal(__func__, "alignment failed its score checks at 32 bits (%s%s%s): qlen=%d, graph rows=%d, match<=%lld, " - "mismatch>=%lld, gap open=%d,%d, gap extension=%d,%d, floor=%lld, zone=(%lld,%lld], max=%lld. " + err_fatal(__func__, "alignment failed its score checks at 32 bits (%s%s%s%s): qlen=%d, graph rows=%d, match<=%lld, " + "mismatch>=%lld, gap open=%d,%d, gap extension=%d,%d; floor=%lld, zone=(%lld,%lld], range=[%lld,%lld]. " "Scores this large cannot be aligned exactly: scale them down.", aln_fail & ABPOA_FAIL_BT ? "backtrack " : "", aln_fail & ABPOA_FAIL_ZONE ? "zone " : "", - aln_fail & ABPOA_FAIL_BEST ? "best" : "", qlen, gn, (long long)r->M, -(long long)r->X, abpt->gap_open1, - abpt->gap_open2, gap_ext1, gap_ext2, (long long)r->floor, (long long)r->fmax, (long long)r->ztop, (long long)r->u); + aln_fail & ABPOA_FAIL_RANGE ? "range " : "", aln_fail & ABPOA_FAIL_BEST ? "best" : "", qlen, gn, + (long long)r->M, -(long long)r->X, abpt->gap_open1, abpt->gap_open2, gap_ext1, gap_ext2, (long long)r->floor, + (long long)r->fmax, (long long)r->ztop, (long long)r->low, (long long)r->top); if (abpt->verbose >= ABPOA_DEBUG_VERBOSE) fprintf(stderr, "16-bit alignment failed its score checks (%d), redone at 32 bits: qlen=%d, graph rows=%d\n", aln_fail, qlen, gn); *res = res0; diff --git a/tests/test_score_range.sh b/tests/test_score_range.sh index 15316a8..074e6eb 100755 --- a/tests/test_score_range.sh +++ b/tests/test_score_range.sh @@ -14,6 +14,8 @@ # every score halved # 5. a 16-bit alignment proven safe needs no retry, even where row 0's gap on the dearer piece # leaves the 16-bit range +# 6. scores too large for the 32-bit bounds to prove an alignment are aligned anyway, checked, and +# give the same result as the scores scaled down; only scores 32 bits cannot check are refused set -e SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)" @@ -135,6 +137,23 @@ for bin in $BINS; do [ "$n" = 0 ] || fail "$(basename "$bin") $R0: a proven 16-bit alignment was redone at 32 bits" done +# 6: cactus's default scores times 10000 on a banded 1.2 kb window: beyond what the 32-bit bounds can +# prove (qlen*(M+X) is above 2^31), so every alignment is checked, and the result must be that of the +# scores themselves. A mismatch of 8e8 leaves no room for the checks in 32 bits: refused up front. +fa="$BUILD_DIR/big.fa" +printf ">a\n%s\n>b\n%s\n>c\n%s\n" "$(randseq 11 1200)" "$(randseq 11 1200 0.1)" "$(randseq 11 1200 0.2)" > "$fa" +for bin in $BINS; do + for opts in "-r1" "-r1 -m 2"; do + a=$("$bin" $opts -M 100 -X 125 -O 400,1200 -E 30,1 "$fa" 2>/dev/null) || true + b=$("$bin" $opts -M 1000000 -X 1250000 -O 4000000,12000000 -E 300000,10000 "$fa" 2>/dev/null) || true + [ -n "$a" ] && [ "$a" = "$b" ] || fail "$(basename "$bin") $opts: cactus's scores times 10000 differ from the scores" + done + if "$bin" -M 1 -X 800000000 -O 4,24 -E 2,1 -r1 test_data/seq.fa >/dev/null 2>"$BUILD_DIR/large.err" || + ! grep -q "scores too large" "$BUILD_DIR/large.err"; then + fail "$(basename "$bin"): a mismatch of 8e8 was not refused" + fi +done + if [ $ERRORS -gt 0 ]; then echo "$ERRORS score range errors" exit 1 From a41d204492ffacc631ed035eb6f7374f2fe7979f Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 17:10:15 -0400 Subject: [PATCH 7/9] Test row 0 on an input where its wrap gave a wrong result 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- tests/test_score_range.sh | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/tests/test_score_range.sh b/tests/test_score_range.sh index 074e6eb..4f6233a 100755 --- a/tests/test_score_range.sh +++ b/tests/test_score_range.sh @@ -12,8 +12,8 @@ # 4. scores that the previous 16-bit choice could not hold (a trained model with gap extension # 60, which aborted with "Error in cg_backtrack") give the 32-bit result, and the same as with # every score halved -# 5. a 16-bit alignment proven safe needs no retry, even where row 0's gap on the dearer piece -# leaves the 16-bit range +# 5. a 16-bit alignment proven safe needs no retry, and gives the 32-bit result, even where row 0's +# gap on the dearer piece leaves the 16-bit range # 6. scores too large for the 32-bit bounds to prove an alignment are aligned anyway, checked, and # give the same result as the scores scaled down; only scores 32 bits cannot check are refused set -e @@ -136,6 +136,14 @@ for bin in $BINS; do n=$("$bin" $R0 -V 2 "$fa" 2>&1 >/dev/null | grep -c "redone at 32 bits" || true) [ "$n" = 0 ] || fail "$(basename "$bin") $R0: a proven 16-bit alignment was redone at 32 bits" done +# the same with a trained-like model whose first gap extension (346) is above the match score (210), +# on divergent sequences: before row 0 was clamped, no check fired here and the result was wrong +fa="$BUILD_DIR/row0b.fa"; : > "$fa" +for j in $(seq 1 12); do printf ">m%d\n%s\n" $j "$(randseq 10 150 "$(awk -v j=$j 'BEGIN { print 0.05 * j }')")" >> "$fa"; done +for j in 1 2 3; do printf ">u%d\n%s\n" $j "$(randseq $((1000 + j)) $((20 * j)))" >> "$fa"; done +for bin in $BINS; do + same_as_plain "$bin" "$fa" -M 210 -X 1020 -O 871,2479 -E 346,2 -b 1000 -f 0.1 -p -r1 +done # 6: cactus's default scores times 10000 on a banded 1.2 kb window: beyond what the 32-bit bounds can # prove (qlen*(M+X) is above 2^31), so every alignment is checked, and the result must be that of the From e2aa0f16579dee0795db0eccdf0af9562ed37a94 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 18:00:33 -0400 Subject: [PATCH 8/9] Align the DP's loops to 32 bytes on x86 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- Makefile | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/Makefile b/Makefile index f5d7aa1..0132b9e 100644 --- a/Makefile +++ b/Makefile @@ -93,6 +93,13 @@ endif endif endif +# Loops of the DP aligned to 32 bytes on x86: on CPUs with Intel's JCC erratum (Skylake and its derivatives) the +# speed of the DP otherwise depends on where unrelated code changes happen to put its loops, by up to 15%. +# Used only if the compiler takes the flag. +ifeq ($(ARCH), $(filter $(ARCH), x86_64 amd64)) + LOOP_ALIGN := $(shell $(CC) -falign-loops=32 -Werror -x c -c /dev/null -o /dev/null 2>/dev/null && echo -falign-loops=32) +endif + # some more possible overrides FLAG_SSE2 = -msse2 FLAG_SSE41 = -msse4.1 @@ -182,16 +189,16 @@ $(SRC_DIR)/simd_check.o:$(SRC_DIR)/simd_check.c $(SRC_DIR)/simd_instruction.h $(CC) -c $(CFLAGS) $(SIMD_FLAG) -I$(INC_DIR) $< -o $@ $(SRC_DIR)/abpoa_align_simd.o:$(SRC_DIR)/abpoa_align_simd.c $(SRC_DIR)/abpoa_graph.h $(SRC_DIR)/abpoa_align.h $(SRC_DIR)/simd_instruction.h $(SRC_DIR)/utils.h - $(CC) -c $(CFLAGS) $(SIMD_FLAG) -I$(INC_DIR) $< -o $@ + $(CC) -c $(CFLAGS) $(LOOP_ALIGN) $(SIMD_FLAG) -I$(INC_DIR) $< -o $@ $(SRC_DIR)/abpoa_align_simd_sse2.o:$(SRC_DIR)/abpoa_align_simd.c $(SRC_DIR)/abpoa_graph.h $(SRC_DIR)/abpoa_align.h $(SRC_DIR)/simd_instruction.h $(SRC_DIR)/utils.h - $(CC) -c $(CFLAGS) -DABPOA_SIMD_DISPATCH -msse2 -U__SSE4_1__ -I$(INC_DIR) $< -o $@ + $(CC) -c $(CFLAGS) $(LOOP_ALIGN) -DABPOA_SIMD_DISPATCH -msse2 -U__SSE4_1__ -I$(INC_DIR) $< -o $@ $(SRC_DIR)/abpoa_align_simd_sse41.o:$(SRC_DIR)/abpoa_align_simd.c $(SRC_DIR)/abpoa_graph.h $(SRC_DIR)/abpoa_align.h $(SRC_DIR)/simd_instruction.h $(SRC_DIR)/utils.h - $(CC) -c $(CFLAGS) -DABPOA_SIMD_DISPATCH -msse4.1 -I$(INC_DIR) $< -o $@ + $(CC) -c $(CFLAGS) $(LOOP_ALIGN) -DABPOA_SIMD_DISPATCH -msse4.1 -I$(INC_DIR) $< -o $@ $(SRC_DIR)/abpoa_align_simd_avx2.o:$(SRC_DIR)/abpoa_align_simd.c $(SRC_DIR)/abpoa_graph.h $(SRC_DIR)/abpoa_align.h $(SRC_DIR)/simd_instruction.h $(SRC_DIR)/utils.h - $(CC) -c $(CFLAGS) -DABPOA_SIMD_DISPATCH -mavx2 -I$(INC_DIR) $< -o $@ + $(CC) -c $(CFLAGS) $(LOOP_ALIGN) -DABPOA_SIMD_DISPATCH -mavx2 -I$(INC_DIR) $< -o $@ $(SRC_DIR)/abpoa_align_simd_avx512bw.o:$(SRC_DIR)/abpoa_align_simd.c $(SRC_DIR)/abpoa_graph.h $(SRC_DIR)/abpoa_align.h $(SRC_DIR)/simd_instruction.h $(SRC_DIR)/utils.h - $(CC) -c $(CFLAGS) -DABPOA_SIMD_DISPATCH -mavx512bw -I$(INC_DIR) $< -o $@ + $(CC) -c $(CFLAGS) $(LOOP_ALIGN) -DABPOA_SIMD_DISPATCH -mavx512bw -I$(INC_DIR) $< -o $@ $(SRC_DIR)/abpoa_dispatch_simd_sse2.o:$(SRC_DIR)/abpoa_dispatch_simd.c $(SRC_DIR)/abpoa.h $(CC) -c $(CFLAGS) -DABPOA_SIMD_DISPATCH -msse2 -U__SSE4_1__ -I$(INC_DIR) $< -o $@ $(SRC_DIR)/abpoa_dispatch_simd_sse41.o:$(SRC_DIR)/abpoa_dispatch_simd.c $(SRC_DIR)/abpoa.h From 12f99c4bc2a9f92546d583e111bd15245a21da23 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 11 Oct 2026 21:00:34 -0400 Subject: [PATCH 9/9] Say the scores are in range when a proven 32-bit alignment fails 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 Claude-Session: https://claude.ai/code/session_013YWGkqE3ND862JF8QRZERV --- src/abpoa_align_simd.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/abpoa_align_simd.c b/src/abpoa_align_simd.c index 2de328f..1b96ad7 100644 --- a/src/abpoa_align_simd.c +++ b/src/abpoa_align_simd.c @@ -1522,8 +1522,9 @@ int simd_abpoa_align_sequence_to_subgraph(abpoa_t *ab, abpoa_para_t *abpt, int b } } if (aln_fail) { - if (bits == 32 && aln_fail == ABPOA_FAIL_BT && r->proven) // the scores provably fit: not an overflow - err_fatal(__func__, "Error in backtrack (32-bit DP, scores in range): qlen=%d, graph rows=%d.", qlen, gn); + if (bits == 32 && r->proven) // the scores provably fit: not an overflow + err_fatal(__func__, "Error in %s (32-bit DP, scores in range): qlen=%d, graph rows=%d.", + aln_fail == ABPOA_FAIL_BT ? "backtrack" : "alignment, its score checks failed", qlen, gn); if (bits == 32) err_fatal(__func__, "alignment failed its score checks at 32 bits (%s%s%s%s): qlen=%d, graph rows=%d, match<=%lld, " "mismatch>=%lld, gap open=%d,%d, gap extension=%d,%d; floor=%lld, zone=(%lld,%lld], range=[%lld,%lld]. "