From 5e202afb5f7f181c3f6ae6978b278661247fcec8 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@me.com>
Date: Wed, 28 Jun 2017 18:17:57 -0400
Subject: [PATCH] fixed invalid backtracking, a temporary solution

---
 ksw2.h           |  1 -
 ksw2_extz2_sse.c | 17 ++++++++++++++++-
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/ksw2.h b/ksw2.h
index 14a8c53..035b5f6 100644
--- a/ksw2.h
+++ b/ksw2.h
@@ -99,7 +99,6 @@ static inline void ksw_backtrack(void *km, int is_rot, int is_rev, const uint8_t
 		if (is_rot) r = i + j, tmp = p[r * n_col + i - off[r]];
 		else tmp = p[i * n_col + j - off[i]];
 		which = tmp >> (which << 1) & 3;
-		if (which == 0 && tmp>>6) break;
 		if (which == 0) which = tmp & 3;
 		if (which == 0)      cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
 		else if (which == 1) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i;      // deletion
diff --git a/ksw2_extz2_sse.c b/ksw2_extz2_sse.c
index 1fedb40..af5d7a3 100644
--- a/ksw2_extz2_sse.c
+++ b/ksw2_extz2_sse.c
@@ -95,7 +95,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) {
 		int st = 0, en = tlen - 1, st0, en0, st_, en_;
 		int8_t x1, v1;
-		uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v;
+		uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v, *x8 = (uint8_t*)x, p_en0 = 0;
 		__m128i x1_, v1_;
 		// find the boundaries
 		if (st < r - qlen + 1) st = r - qlen + 1;
@@ -163,6 +163,13 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
 			__m128i *pr = p + r * n_col_ - st_;
 			off[r] = st;
+			if (en0 < r) { // to avoid backtracking out of the band; this assumes a fixed band
+				int8_t a, z = ((uint8_t*)s)[en0] + 2 * qe;
+				a = x8[en0-1] + v8[en0-1];
+				p_en0 = a > z? 1 : 0;
+				z = a > z? a : z;
+				p_en0 |= a - (z - q) > 0? 1<<2 : 0;
+			}
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, xt1, vt1, ut, tmp;
 				__dp_code_block1;
@@ -189,6 +196,13 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		} else { // gap right-alignment
 			__m128i *pr = p + r * n_col_ - st_;
 			off[r] = st;
+			if (en0 < r) {
+				int8_t a, z = ((uint8_t*)s)[en0] + 2 * qe;
+				a = x8[en0-1] + v8[en0-1];
+				p_en0 = a >= z? 1 : 0;
+				z = a >= z? a : z;
+				p_en0 |= a - (z - q) >= 0? 1<<2 : 0;
+			}
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, xt1, vt1, ut, tmp;
 				__dp_code_block1;
@@ -213,6 +227,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				_mm_store_si128(&pr[t], d);
 			}
 		}
+		if (with_cigar && en0 < r) ((uint8_t*)(p + r * n_col_))[en0 - st] = p_en0;
 		if (!approx_max) { // find the exact max with a 32-bit score array
 			int32_t max_H, max_t;
 			// compute H[], max_H and max_t