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 diff --git a/src/abpoa_align_simd.c b/src/abpoa_align_simd.c index 7ee8917..1b96ad7 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);*/ \ } @@ -598,11 +604,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 +624,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 +645,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 { \ @@ -649,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) { \ @@ -660,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); \ } \ } \ } @@ -677,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]; \ } \ } \ } @@ -697,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]); \ } \ } \ @@ -1140,6 +1158,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,6 +1168,65 @@ } \ } +/* 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; \ + } \ +} + +#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, 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], 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, 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)); 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) 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 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->top || best_score <= r->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, \ @@ -1178,7 +1258,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); \ } @@ -1211,7 +1292,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); \ } @@ -1244,10 +1326,83 @@ 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 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. + */ +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); + } + // 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.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 + 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; +} + // 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 +1454,49 @@ 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(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, (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: + r = bits == 16 ? &r16 : &r32; + if (!r->fits) + 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) { + 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 +1521,28 @@ 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 && 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]. " + "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_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; + 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 diff --git a/tests/test_score_range.sh b/tests/test_score_range.sh new file mode 100755 index 0000000..4f6233a --- /dev/null +++ b/tests/test_score_range.sh @@ -0,0 +1,170 @@ +#!/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 +# 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 + +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 + +# 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 +# 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 +# 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 +fi +echo "Score range checks passed" +exit 0