From 8053bfc4e299f02f9b6e6681f870dfa98cfd46ba Mon Sep 17 00:00:00 2001 From: unavailable-2374 Date: Wed, 28 Jan 2026 18:23:51 +0000 Subject: [PATCH 1/3] Fix wave algorithm boundary conditions to prevent invalid alignments Root cause analysis: - Wave algorithm uses diagonal coordinate system where k = x - y - For valid coordinates: x = (anti+k)/2 >= 0 and y = (anti-k)/2 >= 0 - This requires k to be in range [-anti, anti] - Original code only checked hgh <= anti (y >= 0), not low >= -anti (x >= 0) Fixes applied: 1. Entry point validation (Local_Alignment, Wrap_Around_Alignment, Find_Extension): - Add check: while (((anti+low)>>1) < 0) low += 1 - Constrain minp to -anti (prevents wave expansion into x < 0 region) - Constrain maxp to anti (prevents wave expansion into y < 0 region) - Return error if no valid diagonals remain (low > hgh) 2. Internal bounds checking (all 6 wave functions): - Add check: if (x < 0 || x < k) mark diagonal as dead - Applied to both 0-wave initialization and main wave loops 3. Safety checks (forward_wave, reverse_wave): - Add check: if (avail == 0) return error - Catches case where all diagonals fail bounds checking 4. Early termination (all wave expansion loops): - Add check: if (hgh < low) break - Prevents infinite loops when all diagonals are trimmed This multi-layer fix prevents invalid coordinates from being computed, rather than just detecting and skipping them after the fact. Co-Authored-By: Claude Opus 4.5 --- align.c | 169 +++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 161 insertions(+), 8 deletions(-) diff --git a/align.c b/align.c index 7f6f16e..db4d279 100644 --- a/align.c +++ b/align.c @@ -432,6 +432,16 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { V[k] = -2; // Mark dead + T[k] = PATH_INT; + M[k] = PATH_LEN; + HA[k] = -1; + NA[k] = 0; + bs += 1; + continue; + } if (avail >= cmax-1) { cmax = ((int) (avail*1.2)) + 10000; @@ -536,6 +546,11 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, bclip = -INT32_MAX; } + // Safety check: if no valid diagonals were processed in 0-wave, the alignment is invalid + // This can happen when all diagonals fail bounds checking (x < 0 || x < k) + if (avail == 0) + return (2); // Return error code 2 = no valid alignment found + #ifdef DEBUG_WAVE printf("\nFORWARD WAVE:\n"); print_wave(V,M,low,hgh,besta); @@ -679,6 +694,15 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { t = T[k]; + n = M[k]; + ua = HA[k]; + V[k] = -2; // Mark dead, will be trimmed + bs += 1; + continue; + } while (1) { c = bs[x]; if (c == 4) @@ -789,6 +813,10 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, break; } + // If all diagonals were trimmed (hgh < low), terminate the wave + if (hgh < low) + break; + #ifdef WAVE_STATS k = (hgh-low)+1; if (k > MAX) @@ -956,6 +984,16 @@ static int reverse_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { V[k] = -2; // Mark dead + T[k] = PATH_INT; + M[k] = PATH_LEN; + HA[k] = -1; + NA[k] = 0; + bs -= 1; + continue; + } if (avail >= cmax-1) { cmax = ((int) (avail*1.2)) + 10000; @@ -1199,6 +1237,15 @@ static int reverse_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { t = T[k]; + n = M[k]; + ua = HA[k]; + V[k] = -2; // Mark dead, will be trimmed + bs += 1; + continue; + } while (1) { c = bs[x]; if (c == 4) @@ -1460,14 +1507,27 @@ int Local_Alignment(Alignment *align, Work_Data *ework, Align_Spec *espec, selfie = (align->aseq == align->bseq); + // Ensure valid diagonal range: k must be in [-anti, anti] + // For k outside this range, x = (anti+k)/2 or y = (anti-k)/2 becomes negative + + // Ensure hgh <= anti (guarantees y >= 0 at k = hgh) while (((anti-hgh)>>1) < 0) hgh -= 1; - + + // Ensure low >= -anti (guarantees x >= 0 at k = low) + while (((anti+low)>>1) < 0) + low += 1; + + // If no valid diagonals remain, return error + if (low > hgh) + return (1); + if (lbord < 0) { if (selfie && low >= 0) minp = 1; else - minp = -INT32_MAX; + // Constrain minp to -anti to prevent wave expansion into invalid x < 0 region + minp = -anti; } else minp = low-lbord; @@ -1475,7 +1535,8 @@ int Local_Alignment(Alignment *align, Work_Data *ework, Align_Spec *espec, { if (selfie && hgh <= 0) maxp = -1; else - maxp = INT32_MAX; + // Constrain maxp to anti to prevent wave expansion into invalid y < 0 region + maxp = anti; } else maxp = hgh+hbord; @@ -1662,6 +1723,16 @@ static int forward_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { V[k] = -2; // Mark dead + T[k] = PATH_INT; + M[k] = PATH_LEN; + HA[k] = -1; + NA[k] = 0; + bs += 1; + continue; + } if (avail >= cmax-1) { cmax = ((int) (avail*1.2)) + 10000; @@ -1992,6 +2063,10 @@ static int forward_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, break; } + // If all diagonals were trimmed (hgh < low), terminate the wave + if (hgh < low) + break; + #ifdef WAVE_STATS k = (hgh-low)+1; if (k > MAX) @@ -2154,6 +2229,16 @@ static int reverse_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { V[k] = -2; // Mark dead + T[k] = PATH_INT; + M[k] = PATH_LEN; + HA[k] = -1; + NA[k] = 0; + bs -= 1; + continue; + } if (avail >= cmax-1) { cmax = ((int) (avail*1.2)) + 10000; @@ -2627,15 +2712,29 @@ int Wrap_Around_Alignment(Alignment *align, Work_Data *ework, Align_Spec *espec, printf("\n"); #endif + // Ensure valid diagonal range: k must be in [-anti, anti] + // For k outside this range, x = (anti+k)/2 or y = (anti-k)/2 becomes negative + + // Ensure hgh <= anti (guarantees y >= 0 at k = hgh) while (((anti-hgh)>>1) < 0) hgh -= 1; - + + // Ensure low >= -anti (guarantees x >= 0 at k = low) + while (((anti+low)>>1) < 0) + low += 1; + + // If no valid diagonals remain, return error + if (low > hgh) + return (1); + if (lbord < 0) - minp = -INT32_MAX; + // Constrain minp to -anti to prevent wave expansion into invalid x < 0 region + minp = -anti; else minp = low-lbord; if (hbord < 0) - maxp = INT32_MAX; + // Constrain maxp to anti to prevent wave expansion into invalid y < 0 region + maxp = anti; else maxp = hgh+hbord; @@ -2793,6 +2892,16 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { V[k] = -2; // Mark dead + T[k] = PATH_INT; + M[k] = PATH_LEN; + HA[k] = -1; + NA[k] = 0; + bs += 1; + continue; + } if (avail >= cmax-1) { cmax = ((int) (avail*1.2)) + 10000; @@ -2897,6 +3006,11 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, bclip = -INT32_MAX; } + // Safety check: if no valid diagonals were processed in 0-wave, the alignment is invalid + // This can happen when all diagonals fail bounds checking (x < 0 || x < k) + if (avail == 0) + return (2); // Return error code 2 = no valid alignment found + #ifdef DEBUG_WAVE printf("\nFORWARD WAVE:\n"); print_wave(V,M,low,hgh,besta); @@ -3040,6 +3154,15 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { t = T[k]; + n = M[k]; + ua = HA[k]; + V[k] = -2; // Mark dead, will be trimmed + bs += 1; + continue; + } while (1) { c = bs[x]; if (c == 4) @@ -3150,6 +3273,10 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, break; } + // If all diagonals were trimmed (hgh < low), terminate the wave + if (hgh < low) + break; + #ifdef WAVE_STATS k = (hgh-low)+1; if (k > MAX) @@ -3311,6 +3438,16 @@ static int reverse_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { V[k] = -2; // Mark dead + T[k] = PATH_INT; + M[k] = PATH_LEN; + HA[k] = -1; + NA[k] = 0; + bs -= 1; + continue; + } if (avail >= cmax-1) { cmax = ((int) (avail*1.2)) + 10000; @@ -3554,6 +3691,15 @@ static int reverse_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + if (x < 0 || x < k) + { t = T[k]; + n = M[k]; + ua = HA[k]; + V[k] = -2; // Mark dead, will be trimmed + bs += 1; + continue; + } while (1) { c = bs[x]; if (c == 4) @@ -3807,12 +3953,19 @@ int Find_Extension(Alignment *align, Work_Data *ework, Align_Spec *espec, printf("\n"); #endif + // Validate diagonal is within valid range [-anti, anti] + // For diag outside this range, x = (anti+diag)/2 or y = (anti-diag)/2 becomes negative + if (((anti+diag)>>1) < 0 || ((anti-diag)>>1) < 0) + return (1); // Invalid diagonal for given anti-diagonal + if (lbord < 0) - minp = -INT32_MAX; + // Constrain minp to -anti to prevent wave expansion into invalid x < 0 region + minp = -anti; else minp = diag-lbord; if (hbord < 0) - maxp = INT32_MAX; + // Constrain maxp to anti to prevent wave expansion into invalid y < 0 region + maxp = anti; else maxp = diag+hbord; From 0001aaf6397128be7f032b4fbc1e2b49e7b3b211 Mon Sep 17 00:00:00 2001 From: unavailable-2374 Date: Mon, 2 Feb 2026 17:06:30 +0000 Subject: [PATCH 2/3] Add coordinate validation in ALNtoPAF to skip invalid alignments Skip alignments where coordinates exceed sequence bounds: - path->abpos < 0 or path->aepos > aln->alen - path->bbpos < 0 or path->bepos > aln->blen - start >= end for either sequence This handles edge cases where wave algorithm produces out-of-bounds coordinates, preventing "Subrange out of bounds" errors in Get_Contig_Piece. Co-Authored-By: Claude Opus 4.5 --- ALNtoPAF.c | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/ALNtoPAF.c b/ALNtoPAF.c index c37cb56..819e9ca 100644 --- a/ALNtoPAF.c +++ b/ALNtoPAF.c @@ -181,6 +181,14 @@ void *gen_paf(void *args) ascaff = contigs1[acontig].scaf; bscaff = contigs2[bcontig].scaf; + // Validate alignment coordinates are within bounds + // Skip alignments with invalid coordinates (can occur from edge cases in wave algorithm) + if (path->abpos < 0 || path->aepos > aln->alen || path->abpos >= path->aepos || + path->bbpos < 0 || path->bepos > aln->blen || path->bbpos >= path->bepos) + { alast = acontig; + continue; // Skip this invalid alignment + } + aoff = contigs1[acontig].sbeg; if (SWAP_G) From 7d8e0e882aa6b903e1098491a442a11393e59143 Mon Sep 17 00:00:00 2001 From: unavailable-2374 Date: Wed, 4 Feb 2026 15:59:10 +0000 Subject: [PATCH 3/3] Fix wavefront boundary with dynamic sequence length constraints During each wavefront expansion step, calculate the effective diagonal range based on the current anti-diagonal position: Forward wave: curr_anti = mida + 2*(dif+1) Reverse wave: curr_anti = mida - 2*(dif+1) Sequence length constraints: - seq_minp = curr_anti - 2*blen (ensures y <= blen) - seq_maxp = 2*alen - curr_anti (ensures x <= alen) Take the stricter of user constraints and sequence constraints: - eff_minp = max(minp, seq_minp) - eff_maxp = min(maxp, seq_maxp) Applied to all 6 wavefront functions: - forward_wave, reverse_wave - forward_wrap, reverse_wrap - forward_extend, reverse_extend Also restored simple boundary check: if (x < 0 || x < k) This dynamic approach is more precise than static constraints because the valid diagonal range changes as the wavefront expands. Co-Authored-By: Claude Opus 4.5 --- align.c | 230 +++++++++++++++++++++++++++++++++++--------------------- 1 file changed, 146 insertions(+), 84 deletions(-) diff --git a/align.c b/align.c index db4d279..16b1281 100644 --- a/align.c +++ b/align.c @@ -353,6 +353,8 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, int *mind, int maxd, int mida, int minp, int maxp, int aoff) { char *aseq = align->aseq; char *bseq = align->bseq; + int alen = align->alen; + int blen = align->blen; Path *apath = align->path; int hgh, low, dif; @@ -432,7 +434,7 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { V[k] = -2; // Mark dead T[k] = PATH_INT; @@ -630,19 +632,27 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, T = _T-vmin; } - if (low >= minp) - { NA[low] = NA[low+1]; - V[low] = -1; - } - else - low += 1; + // Dynamic sequence length constraint for wave expansion + { int curr_anti = mida + 2*(dif+1); + int seq_minp = curr_anti - 2*blen; + int seq_maxp = 2*alen - curr_anti; + int eff_minp = (minp > seq_minp) ? minp : seq_minp; + int eff_maxp = (maxp < seq_maxp) ? maxp : seq_maxp; - if (hgh <= maxp) - { NA[hgh] = NA[hgh-1]; - V[hgh] = am = -1; - } - else - am = V[--hgh]; + if (low >= eff_minp) + { NA[low] = NA[low+1]; + V[low] = -1; + } + else + low += 1; + + if (hgh <= eff_maxp) + { NA[hgh] = NA[hgh-1]; + V[hgh] = am = -1; + } + else + am = V[--hgh]; + } dif += 1; @@ -694,7 +704,7 @@ static int forward_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { t = T[k]; n = M[k]; @@ -907,6 +917,8 @@ static int reverse_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, int mind, int maxd, int mida, int minp, int maxp, int aoff) { char *aseq = align->aseq - 1; char *bseq = align->bseq - 1; + int alen = align->alen; + int blen = align->blen; Path *apath = align->path; int hgh, low, dif; @@ -984,7 +996,7 @@ static int reverse_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { V[k] = -2; // Mark dead T[k] = PATH_INT; @@ -1173,19 +1185,27 @@ static int reverse_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, T = _T-vmin; } - if (low >= minp) - { NA[low] = NA[low+1]; - V[low] = ap = INT32_MAX; - } - else - ap = V[++low]; + // Dynamic sequence length constraint for wave expansion + { int curr_anti = mida - 2*(dif+1); + int seq_minp = curr_anti - 2*blen; + int seq_maxp = 2*alen - curr_anti; + int eff_minp = (minp > seq_minp) ? minp : seq_minp; + int eff_maxp = (maxp < seq_maxp) ? maxp : seq_maxp; - if (hgh <= maxp) - { NA[hgh] = NA[hgh-1]; - V[hgh] = INT32_MAX; - } - else - hgh -= 1; + if (low >= eff_minp) + { NA[low] = NA[low+1]; + V[low] = ap = INT32_MAX; + } + else + ap = V[++low]; + + if (hgh <= eff_maxp) + { NA[hgh] = NA[hgh-1]; + V[hgh] = INT32_MAX; + } + else + hgh -= 1; + } dif += 1; @@ -1237,7 +1257,7 @@ static int reverse_wave(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { t = T[k]; n = M[k]; @@ -1478,10 +1498,12 @@ int Local_Alignment(Alignment *align, Work_Data *ework, Align_Spec *espec, int selfie; int fshort, rshort; - { int alen; + { int alen, blen; int maxtp, wsize; alen = align->alen; + blen = align->blen; + (void) blen; // Used by wave functions via align->blen if (hgh-low >= 7500) wsize = VectorEl*(hgh-low+1); @@ -1647,6 +1669,8 @@ static int forward_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, int *mind, int maxd, int mida, int minp, int maxp, int tspace) { char *aseq = align->aseq; char *bseq = align->bseq; + int alen = align->alen; + int blen = align->blen; Path *apath = align->path; int hgh, low, dif; @@ -1723,7 +1747,7 @@ static int forward_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { V[k] = -2; // Mark dead T[k] = PATH_INT; @@ -1903,19 +1927,27 @@ static int forward_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, T = _T-vmin; } - if (low >= minp) - { NA[low] = NA[low+1]; - V[low] = -1; - } - else - low += 1; + // Dynamic sequence length constraint for wave expansion + { int curr_anti = mida + 2*(dif+1); + int seq_minp = curr_anti - 2*blen; + int seq_maxp = 2*alen - curr_anti; + int eff_minp = (minp > seq_minp) ? minp : seq_minp; + int eff_maxp = (maxp < seq_maxp) ? maxp : seq_maxp; - if (hgh <= maxp) - { NA[hgh] = NA[hgh-1]; - V[hgh] = am = -1; - } - else - am = V[--hgh]; + if (low >= eff_minp) + { NA[low] = NA[low+1]; + V[low] = -1; + } + else + low += 1; + + if (hgh <= eff_maxp) + { NA[hgh] = NA[hgh-1]; + V[hgh] = am = -1; + } + else + am = V[--hgh]; + } dif += 1; @@ -2151,10 +2183,12 @@ static int forward_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, return (0); } -static int reverse_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, +static int reverse_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, int mind, int maxd, int mida, int minp, int maxp, int tspace) { char *aseq = align->aseq - 1; char *bseq = align->bseq - 1; + int alen = align->alen; + int blen = align->blen; Path *apath = align->path; int hgh, low, dif; @@ -2229,7 +2263,7 @@ static int reverse_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { V[k] = -2; // Mark dead T[k] = PATH_INT; @@ -2405,19 +2439,27 @@ static int reverse_wrap(_Work_Data *work, _Align_Spec *spec, Alignment *align, T = _T-vmin; } - if (low >= minp) - { NA[low] = NA[low+1]; - V[low] = ap = INT32_MAX; - } - else - ap = V[++low]; + // Dynamic sequence length constraint for wave expansion + { int curr_anti = mida - 2*(dif+1); + int seq_minp = curr_anti - 2*blen; + int seq_maxp = 2*alen - curr_anti; + int eff_minp = (minp > seq_minp) ? minp : seq_minp; + int eff_maxp = (maxp < seq_maxp) ? maxp : seq_maxp; - if (hgh <= maxp) - { NA[hgh] = NA[hgh-1]; - V[hgh] = INT32_MAX; - } - else - hgh -= 1; + if (low >= eff_minp) + { NA[low] = NA[low+1]; + V[low] = ap = INT32_MAX; + } + else + ap = V[++low]; + + if (hgh <= eff_maxp) + { NA[hgh] = NA[hgh-1]; + V[hgh] = INT32_MAX; + } + else + hgh -= 1; + } dif += 1; @@ -2814,6 +2856,8 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, int midd, int mida, int minp, int maxp) { char *aseq = align->aseq; char *bseq = align->bseq; + int alen = align->alen; + int blen = align->blen; Path *apath = align->path; int hgh, low, dif; @@ -2892,7 +2936,7 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { V[k] = -2; // Mark dead T[k] = PATH_INT; @@ -3090,19 +3134,27 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, T = _T-vmin; } - if (low >= minp) - { NA[low] = NA[low+1]; - V[low] = -1; - } - else - low += 1; + // Dynamic sequence length constraint for wave expansion + { int curr_anti = mida + 2*(dif+1); + int seq_minp = curr_anti - 2*blen; + int seq_maxp = 2*alen - curr_anti; + int eff_minp = (minp > seq_minp) ? minp : seq_minp; + int eff_maxp = (maxp < seq_maxp) ? maxp : seq_maxp; - if (hgh <= maxp) - { NA[hgh] = NA[hgh-1]; - V[hgh] = am = -1; - } - else - am = V[--hgh]; + if (low >= eff_minp) + { NA[low] = NA[low+1]; + V[low] = -1; + } + else + low += 1; + + if (hgh <= eff_maxp) + { NA[hgh] = NA[hgh-1]; + V[hgh] = am = -1; + } + else + am = V[--hgh]; + } dif += 1; @@ -3154,7 +3206,7 @@ static int forward_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { t = T[k]; n = M[k]; @@ -3362,6 +3414,8 @@ static int reverse_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, int midd, int mida, int minp, int maxp) { char *aseq = align->aseq - 1; char *bseq = align->bseq - 1; + int alen = align->alen; + int blen = align->blen; Path *apath = align->path; int hgh, low, dif; @@ -3438,7 +3492,7 @@ static int reverse_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, Pebble *pb; x = (mida+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { V[k] = -2; // Mark dead T[k] = PATH_INT; @@ -3627,19 +3681,27 @@ static int reverse_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, T = _T-vmin; } - if (low >= minp) - { NA[low] = NA[low+1]; - V[low] = ap = INT32_MAX; - } - else - ap = V[++low]; + // Dynamic sequence length constraint for wave expansion + { int curr_anti = mida - 2*(dif+1); + int seq_minp = curr_anti - 2*blen; + int seq_maxp = 2*alen - curr_anti; + int eff_minp = (minp > seq_minp) ? minp : seq_minp; + int eff_maxp = (maxp < seq_maxp) ? maxp : seq_maxp; - if (hgh <= maxp) - { NA[hgh] = NA[hgh-1]; - V[hgh] = INT32_MAX; - } - else - hgh -= 1; + if (low >= eff_minp) + { NA[low] = NA[low+1]; + V[low] = ap = INT32_MAX; + } + else + ap = V[++low]; + + if (hgh <= eff_maxp) + { NA[hgh] = NA[hgh-1]; + V[hgh] = INT32_MAX; + } + else + hgh -= 1; + } dif += 1; @@ -3691,7 +3753,7 @@ static int reverse_extend(_Work_Data *work, _Align_Spec *spec, Alignment *align, b <<= 1; x = (c+k)>>1; - // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]) + // Bounds check: x must be >= 0 (for aseq[x]) and >= k (for bs[x]=bseq[y]) if (x < 0 || x < k) { t = T[k]; n = M[k];