forward-backward search seems working

Makefile

@@ -34,6 +34,7 @@ bwt1away.o:bwt.h bwtaln.h
 bwt2fmv.o:bwt.h
 bntseq.o:bntseq.h
 bwtgap.o:bwtgap.h bwtaln.h bwt.h
+fastmap:bwt.h
 
 bwtsw2_core.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h
 bwtsw2_aux.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h

bwt.c

@@ -249,13 +249,13 @@ void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_b
 	int i;
 	bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl);
 	for (i = 0; i != 4; ++i) {
-		ok[i].x[!is_back] = bwt->L2[i] + tk[i] + 1;
-		ok[i].x[2] = (tl[i] -= tk[i]);
+		ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i];
+		ok[i].x[2] = tl[i] - tk[i];
 	}
-	ok[3].x[is_back] = ik->x[is_back];
-	ok[2].x[is_back] = ok[3].x[is_back] + tl[3];
-	ok[1].x[is_back] = ok[2].x[is_back] + tl[2];
-	ok[0].x[is_back] = ok[1].x[is_back] + tl[1];
+	ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary);
+	ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
+	ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
+	ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];
 }
 
 static void bwt_reverse_intvs(bwtintv_v *p)
@@ -298,7 +298,6 @@ int bwt_smem1(const bwt_t *bwt, int len, const uint8_t *q, int x, bwtintv_v *mem
 	ret = curr->a[0].info; // this will be the returned value
 	swap = curr; curr = prev; prev = swap;
 
-	if (x == 40) printf("[%lld,%lld,%lld]\n", prev->a[0].x[0], prev->a[0].x[1], prev->a[0].x[2]);
 	for (i = x - 1; i >= -1; --i) { // backward search for MEMs
 		if (q[i] > 3) break;
 		c = i < 0? 0 : q[i];

bwt.h

@@ -81,7 +81,7 @@ typedef struct { size_t n, m; bwtintv_t *a; } bwtintv_v;
 	 (bwt)->L2[bwt_B0(bwt, k)] + bwt_occ(bwt, k, bwt_B0(bwt, k))		\
 	 : (bwt)->L2[bwt_B0(bwt, (k)-1)] + bwt_occ(bwt, k, bwt_B0(bwt, (k)-1)))
 
-#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)], (ik).x[2] = (bwt)->L2[(int)(c)+1] - (bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3-(c)], (ik).info = 0)
+#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)]+1, (ik).x[2] = (bwt)->L2[(int)(c)+1]-(bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3-(c)]+1, (ik).info = 0)
 
 #ifdef __cplusplus
 extern "C" {

fastmap.c

@@ -40,7 +40,9 @@ int main_fastmap(int argc, char *argv[])
 		for (i = 0; i < seq->seq.l; ++i)
 			seq->seq.s[i] = nst_nt4_table[(int)seq->seq.s[i]];
 		{
-			int beg = 98;
+			int beg = 97;
+			uint8_t str[3];
+			bwtint_t b, e;
 			bwtintv_t ik, ok[4];
 			bwt_set_intv(bwt, seq->seq.s[seq->seq.l - 1], ik);
 			for (i = seq->seq.l - 2; i >= beg; --i) {
@@ -50,6 +52,7 @@ int main_fastmap(int argc, char *argv[])
 				if (ik.x[2] == 0) break;
 			}
 			printf("[%lld,%lld,%lld] @ %d\n", ik.x[0], ik.x[1], ik.x[2], i+1);
+			//str[0] = '\1'; str[1] = '\3'; bwt_match_exact(bwt, 2, str, &b, &e); printf("%lld,%lld,%lld\n", b, e, e-b+1);
 			printf("======================== %lld, [%lld,%lld,%lld,%lld]\n", bwt->primary, bwt->L2[1], bwt->L2[2]-bwt->L2[1], bwt->L2[3]-bwt->L2[2], bwt->L2[4]-bwt->L2[3]);
 			bwt_set_intv(bwt, seq->seq.s[beg], ik);
 			for (i = beg + 1; i < seq->seq.l; ++i) {
@@ -59,6 +62,7 @@ int main_fastmap(int argc, char *argv[])
 				if (ik.x[2] == 0) break;
 			}
 			printf("[%lld,%lld,%lld] @ %d\n", ik.x[0], ik.x[1], ik.x[2], i-1);
+			//str[0] = '\0'; str[1] = '\2'; bwt_match_exact(bwt, 2, str, &b, &e); printf("%lld,%lld,%lld\n", b, e, e-b+1);
 		}
 		bwt_smem(bwt, seq->seq.l, (uint8_t*)seq->seq.s, &mem, tvec);
 		printf(">%s\t%ld\n", seq->name.s, mem.n);