fixed a couple of leaks; buggy atm

bwamem.c

@@ -10,7 +10,7 @@
 void mem_fill_scmat(int a, int b, int8_t mat[25])
 {
 	int i, j, k;
-	for (i = k = 0; i < 5; ++i) {
+	for (i = k = 0; i < 4; ++i) {
 		for (j = 0; j < 4; ++j)
 			mat[k++] = i == j? a : -b;
 		mat[k++] = 0; // ambiguous base
@@ -233,11 +233,28 @@ mem_aln_t mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac,
 	s = &c->seeds[0];
 	if (s->qbeg + s->len != l_query) { // right extension of the first seed
 		int qle, tle, qe, re;
+		int16_t *qw = 0;
 		qe = s->qbeg + s->len; re = s->rbeg + s->len - rmax[0];
 		for (j = 0; j < l_query - qe; ++j) putchar("ACGTN"[(int)query[j+qe]]); putchar('\n');
 		for (j = 0; j < rmax[1] - rmax[0] - re; ++j) putchar("ACGTN"[(int)rseq[j+re]]); putchar('\n');
-		a.score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, opt->w, a.score, 0, &qle, &tle);
+		if (c->n > 1) { // generate $qw
+			int l = rmax[1] - (s->rbeg + s->len);
+			assert(l >= 0 && l < 1000);
+			qw = malloc(l * 2);
+			for (i = 0; i < l; ++i) qw[i] = -1; // no constraint by default
+			for (i = 1; i < c->n; ++i) {
+				const mem_seed_t *t = &c->seeds[i];
+				for (j = 0; j < t->len; ++j) {
+					int x = t->rbeg + j - (s->rbeg + s->len), y = t->qbeg + j - (s->qbeg + s->len);
+					assert(x < l);
+					if (qw[x] == -1) qw[x] = x > y? x - y : y - x;
+					else if (qw[x] >= 0) qw[x] = -2; // in a seed overlap, do not set any constraint
+				}
+			}
+		}
+		a.score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, opt->w, a.score, qw, &qle, &tle);
 		a.qe = qe + qle; a.re = rmax[0] + re + tle;
+		free(qw);
 	} else a.qe = l_query, a.re = s->rbeg + s->len;
 
 	a.is_all = 1;

fastmap.c

@@ -17,7 +17,7 @@ int main_mem(int argc, char *argv[])
 	bwt_t *bwt;
 	bntseq_t *bns;
 	int i, j, c;
-	gzFile *fp;
+	gzFile fp;
 	kseq_t *seq;
 	uint8_t *pac = 0;
 
@@ -66,9 +66,15 @@ int main_mem(int argc, char *argv[])
 			putchar('\n');
 		}
 		puts("//");
+		for (i = 0; i < chain.n; ++i) free(chain.chains[i].seeds);
+		free(chain.chains);
 	}
 
-	free(opt);
+	free(pac); free(opt);
+	bns_destroy(bns);
+	bwt_destroy(bwt);
+	kseq_destroy(seq);
+	gzclose(fp);
 	return 0;
 }
 

ksw.c

@@ -315,7 +315,7 @@ typedef struct {
 	int32_t h, e;
 } eh_t;
 
-int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int *qw, int *_qle, int *_tle)
+int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int16_t *qw, int *_qle, int *_tle)
 {
 	eh_t *eh;
 	int8_t *qp;
@@ -350,7 +350,7 @@ int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target,
 		h1 = h0 - (gapo + gape * (i + 1));
 		if (h1 < 0) h1 = 0;
 		// apply the band and the constraint (if provided)
-		t = (qw && qw[i] < w)? qw[i] : w; // this is the band width at $i
+		t = (qw && qw[i] >= 0 && qw[i] < w)? qw[i] : w; // this is the band width at $i
 		if (beg < i - t) beg = i - t;
 		if (end > i + t + 1) end = i + t + 1;
 		if (end > qlen) end = qlen;

ksw.h

@@ -49,7 +49,7 @@ extern "C" {
 	/** Unified interface for ksw_sse2_8() and ksw_sse2_16() */
 	int ksw_sse2(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a);
 
-	int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int *qw, int *_qle, int *_tle);
+	int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int16_t *qw, int *_qle, int *_tle);
 
 #ifdef __cplusplus
 }