1 #include <stdlib.h>
2 #include <string.h>
3 #include <stdio.h>
4 #include <sys/resource.h>
5 #include <assert.h>
6 #include "bwt_lite.h"
7 #include "bwtsw2.h"
8 #include "bwt.h"
9 #include "kvec.h"
10 
11 #ifdef USE_MALLOC_WRAPPERS
12 #  include "malloc_wrap.h"
13 #endif
14 
15 typedef struct {
16 	bwtint_t k, l;
17 } qintv_t;
18 
19 #define qintv_eq(a, b) ((a).k == (b).k && (a).l == (b).l)
20 #define qintv_hash(a) ((a).k>>7^(a).l<<17)
21 
22 #include "khash.h"
23 KHASH_INIT(qintv, qintv_t, uint64_t, 1, qintv_hash, qintv_eq)
24 KHASH_MAP_INIT_INT64(64, uint64_t)
25 
26 #define MINUS_INF -0x3fffffff
27 #define MASK_LEVEL 0.90f
28 
29 struct __mempool_t;
30 static void mp_destroy(struct __mempool_t*);
31 typedef struct {
32 	bwtint_t qk, ql;
33 	int I, D, G;
34 	uint32_t pj:2, qlen:30;
35 	int tlen;
36 	int ppos, upos;
37 	int cpos[4];
38 } bsw2cell_t;
39 
40 #include "ksort.h"
41 KSORT_INIT_GENERIC(int)
42 #define __hitG_lt(a, b) (((a).G + ((int)(a).n_seeds<<2)) > (b).G + ((int)(b).n_seeds<<2))
43 KSORT_INIT(hitG, bsw2hit_t, __hitG_lt)
44 
45 static const bsw2cell_t g_default_cell = { 0, 0, MINUS_INF, MINUS_INF, MINUS_INF, 0, 0, 0, -1, -1, {-1, -1, -1, -1} };
46 
47 typedef struct {
48 	int n, max;
49 	uint32_t tk, tl; // this is fine
50 	bsw2cell_t *array;
51 } bsw2entry_t, *bsw2entry_p;
52 
53 /* --- BEGIN: Stack operations --- */
54 typedef struct {
55 	int n_pending;
56 	kvec_t(bsw2entry_p) stack0, pending;
57 	struct __mempool_t *pool;
58 } bsw2stack_t;
59 
60 #define stack_isempty(s) (kv_size(s->stack0) == 0 && s->n_pending == 0)
stack_destroy(bsw2stack_t * s)61 static void stack_destroy(bsw2stack_t *s) { mp_destroy(s->pool); kv_destroy(s->stack0); kv_destroy(s->pending); free(s); }
stack_push0(bsw2stack_t * s,bsw2entry_p e)62 inline static void stack_push0(bsw2stack_t *s, bsw2entry_p e) { kv_push(bsw2entry_p, s->stack0, e); }
stack_pop(bsw2stack_t * s)63 inline static bsw2entry_p stack_pop(bsw2stack_t *s)
64 {
65 	assert(!(kv_size(s->stack0) == 0 && s->n_pending != 0));
66 	return kv_pop(s->stack0);
67 }
68 /* --- END: Stack operations --- */
69 
70 /* --- BEGIN: memory pool --- */
71 typedef struct __mempool_t {
72 	int cnt; // if cnt!=0, then there must be memory leak
73 	kvec_t(bsw2entry_p) pool;
74 } mempool_t;
mp_alloc(mempool_t * mp)75 inline static bsw2entry_p mp_alloc(mempool_t *mp)
76 {
77 	++mp->cnt;
78 	if (kv_size(mp->pool) == 0) return (bsw2entry_t*)calloc(1, sizeof(bsw2entry_t));
79 	else return kv_pop(mp->pool);
80 }
mp_free(mempool_t * mp,bsw2entry_p e)81 inline static void mp_free(mempool_t *mp, bsw2entry_p e)
82 {
83 	--mp->cnt; e->n = 0;
84 	kv_push(bsw2entry_p, mp->pool, e);
85 }
mp_destroy(struct __mempool_t * mp)86 static void mp_destroy(struct __mempool_t *mp)
87 {
88 	int i;
89 	for (i = 0; i != kv_size(mp->pool); ++i) {
90 		free(kv_A(mp->pool, i)->array);
91 		free(kv_A(mp->pool, i));
92 	}
93 	kv_destroy(mp->pool);
94 	free(mp);
95 }
96 /* --- END: memory pool --- */
97 
98 /* --- BEGIN: utilities --- */
bsw2_connectivity(const bwtl_t * b)99 static khash_t(64) *bsw2_connectivity(const bwtl_t *b)
100 {
101 	khash_t(64) *h;
102 	uint32_t k, l, cntk[4], cntl[4]; // this is fine
103 	uint64_t x;
104 	khiter_t iter;
105 	int j, ret;
106 	kvec_t(uint64_t) stack;
107 
108 	kv_init(stack);
109 	h = kh_init(64);
110 	kh_resize(64, h, b->seq_len * 4);
111 	x = b->seq_len;
112 	kv_push(uint64_t, stack, x);
113 	while (kv_size(stack)) {
114 		x = kv_pop(stack);
115 		k = x>>32; l = (uint32_t)x;
116 		bwtl_2occ4(b, k-1, l, cntk, cntl);
117 		for (j = 0; j != 4; ++j) {
118 			k = b->L2[j] + cntk[j] + 1;
119 			l = b->L2[j] + cntl[j];
120 			if (k > l) continue;
121 			x = (uint64_t)k << 32 | l;
122 			iter = kh_put(64, h, x, &ret);
123 			if (ret) { // if not present
124 				kh_value(h, iter) = 1;
125 				kv_push(uint64_t, stack, x);
126 			} else ++kh_value(h, iter);
127 		}
128 	}
129 	kv_destroy(stack);
130 	//fprintf(stderr, "[bsw2_connectivity] %u nodes in the DAG\n", kh_size(h));
131 	return h;
132 }
133 // pick up top T matches at a node
cut_tail(bsw2entry_t * u,int T,bsw2entry_t * aux)134 static void cut_tail(bsw2entry_t *u, int T, bsw2entry_t *aux)
135 {
136 	int i, *a, n, x;
137 	if (u->n <= T) return;
138 	if (aux->max < u->n) {
139 		aux->max = u->n;
140 		aux->array = (bsw2cell_t*)realloc(aux->array, aux->max * sizeof(bsw2cell_t));
141 	}
142 	a = (int*)aux->array;
143 	for (i = n = 0; i != u->n; ++i)
144 		if (u->array[i].ql && u->array[i].G > 0)
145 			a[n++] = -u->array[i].G;
146 	if (n <= T) return;
147 	x = -ks_ksmall(int, n, a, T);
148 	n = 0;
149 	for (i = 0; i < u->n; ++i) {
150 		bsw2cell_t *p = u->array + i;
151 		if (p->G == x) ++n;
152 		if (p->G < x || (p->G == x && n >= T)) {
153 			p->qk = p->ql = 0; p->G = 0;
154 			if (p->ppos >= 0) u->array[p->ppos].cpos[p->pj] = -1;
155 		}
156 	}
157 }
158 // remove duplicated cells
remove_duplicate(bsw2entry_t * u,khash_t (qintv)* hash)159 static inline void remove_duplicate(bsw2entry_t *u, khash_t(qintv) *hash)
160 {
161 	int i, ret, j;
162 	khiter_t k;
163 	qintv_t key;
164 	kh_clear(qintv, hash);
165 	for (i = 0; i != u->n; ++i) {
166 		bsw2cell_t *p = u->array + i;
167 		if (p->ql == 0) continue;
168 		key.k = p->qk; key.l = p->ql;
169 		k = kh_put(qintv, hash, key, &ret);
170 		j = -1;
171 		if (ret == 0) {
172 			if ((uint32_t)kh_value(hash, k) >= p->G) j = i;
173 			else {
174 				j = kh_value(hash, k)>>32;
175 				kh_value(hash, k) = (uint64_t)i<<32 | p->G;
176 			}
177 		} else kh_value(hash, k) = (uint64_t)i<<32 | p->G;
178 		if (j >= 0) {
179 			p = u->array + j;
180 			p->qk = p->ql = 0; p->G = 0;
181 			if (p->ppos >= 0) u->array[p->ppos].cpos[p->pj] = -3;
182 		}
183 	}
184 }
185 // merge two entries
merge_entry(const bsw2opt_t * __restrict opt,bsw2entry_t * u,bsw2entry_t * v,bwtsw2_t * b)186 static void merge_entry(const bsw2opt_t * __restrict opt, bsw2entry_t *u, bsw2entry_t *v, bwtsw2_t *b)
187 {
188 	int i;
189 	if (u->n + v->n >= u->max) {
190 		u->max = u->n + v->n;
191 		u->array = (bsw2cell_t*)realloc(u->array, u->max * sizeof(bsw2cell_t));
192 	}
193 	for (i = 0; i != v->n; ++i) {
194 		bsw2cell_t *p = v->array + i;
195 		if (p->ppos >= 0) p->ppos += u->n;
196 		if (p->cpos[0] >= 0) p->cpos[0] += u->n;
197 		if (p->cpos[1] >= 0) p->cpos[1] += u->n;
198 		if (p->cpos[2] >= 0) p->cpos[2] += u->n;
199 		if (p->cpos[3] >= 0) p->cpos[3] += u->n;
200 	}
201 	memcpy(u->array + u->n, v->array, v->n * sizeof(bsw2cell_t));
202 	u->n += v->n;
203 }
204 
push_array_p(bsw2entry_t * e)205 static inline bsw2cell_t *push_array_p(bsw2entry_t *e)
206 {
207 	if (e->n == e->max) {
208 		e->max = e->max? e->max<<1 : 256;
209 		e->array = (bsw2cell_t*)realloc(e->array, sizeof(bsw2cell_t) * e->max);
210 	}
211 	return e->array + e->n;
212 }
213 
time_elapse(const struct rusage * curr,const struct rusage * last)214 static inline double time_elapse(const struct rusage *curr, const struct rusage *last)
215 {
216 	long t1 = (curr->ru_utime.tv_sec - last->ru_utime.tv_sec) + (curr->ru_stime.tv_sec - last->ru_stime.tv_sec);
217 	long t2 = (curr->ru_utime.tv_usec - last->ru_utime.tv_usec) + (curr->ru_stime.tv_usec - last->ru_stime.tv_usec);
218 	return (double)t1 + t2 * 1e-6;
219 }
220 /* --- END: utilities --- */
221 
222 /* --- BEGIN: processing partial hits --- */
save_hits(const bwtl_t * bwt,int thres,bsw2hit_t * hits,bsw2entry_t * u)223 static void save_hits(const bwtl_t *bwt, int thres, bsw2hit_t *hits, bsw2entry_t *u)
224 {
225 	int i;
226 	uint32_t k; // this is fine
227 	for (i = 0; i < u->n; ++i) {
228 		bsw2cell_t *p = u->array + i;
229 		if (p->G < thres) continue;
230 		for (k = u->tk; k <= u->tl; ++k) {
231 			int beg, end;
232 			bsw2hit_t *q = 0;
233 			beg = bwt->sa[k]; end = beg + p->tlen;
234 			if (p->G > hits[beg*2].G) {
235 				hits[beg*2+1] = hits[beg*2];
236 				q = hits + beg * 2;
237 			} else if (p->G > hits[beg*2+1].G) q = hits + beg * 2 + 1;
238 			if (q) {
239 				q->k = p->qk; q->l = p->ql; q->len = p->qlen; q->G = p->G;
240 				q->beg = beg; q->end = end; q->G2 = q->k == q->l? 0 : q->G;
241 				q->flag = q->n_seeds = 0;
242 			}
243 		}
244 	}
245 }
246 /* "narrow hits" are node-to-node hits that have a high score and
247  * are not so repetitive (|SA interval|<=IS). */
save_narrow_hits(const bwtl_t * bwtl,bsw2entry_t * u,bwtsw2_t * b1,int t,int IS)248 static void save_narrow_hits(const bwtl_t *bwtl, bsw2entry_t *u, bwtsw2_t *b1, int t, int IS)
249 {
250 	int i;
251 	for (i = 0; i < u->n; ++i) {
252 		bsw2hit_t *q;
253 		bsw2cell_t *p = u->array + i;
254 		if (p->G >= t && p->ql - p->qk + 1 <= IS) { // good narrow hit
255 			if (b1->max == b1->n) {
256 				b1->max = b1->max? b1->max<<1 : 4;
257 				b1->hits = realloc(b1->hits, b1->max * sizeof(bsw2hit_t));
258 			}
259 			q = &b1->hits[b1->n++];
260 			q->k = p->qk; q->l = p->ql;
261 			q->len = p->qlen;
262 			q->G = p->G; q->G2 = 0;
263 			q->beg = bwtl->sa[u->tk]; q->end = q->beg + p->tlen;
264 			q->flag = 0;
265 			// delete p
266 			p->qk = p->ql = 0; p->G = 0;
267 			if (p->ppos >= 0) u->array[p->ppos].cpos[p->pj] = -3;
268 		}
269 	}
270 }
271 /* after this, "narrow SA hits" will be expanded and the coordinates
272  * will be obtained and stored in b->hits[*].k. */
bsw2_resolve_duphits(const bntseq_t * bns,const bwt_t * bwt,bwtsw2_t * b,int IS)273 int bsw2_resolve_duphits(const bntseq_t *bns, const bwt_t *bwt, bwtsw2_t *b, int IS)
274 {
275 	int i, j, n, is_rev;
276 	if (b->n == 0) return 0;
277 	if (bwt && bns) { // convert to chromosomal coordinates if requested
278 		int old_n = b->n;
279 		bsw2hit_t *old_hits = b->hits;
280 		for (i = n = 0; i < b->n; ++i) { // compute the memory to allocated
281 			bsw2hit_t *p = old_hits + i;
282 			if (p->l - p->k + 1 <= IS) n += p->l - p->k + 1;
283 			else if (p->G > 0) ++n;
284 		}
285 		b->n = b->max = n;
286 		b->hits = calloc(b->max, sizeof(bsw2hit_t));
287 		for (i = j = 0; i < old_n; ++i) {
288 			bsw2hit_t *p = old_hits + i;
289 			if (p->l - p->k + 1 <= IS) { // the hit is no so repetitive
290 				bwtint_t k;
291 				if (p->G == 0 && p->k == 0 && p->l == 0 && p->len == 0) continue;
292 				for (k = p->k; k <= p->l; ++k) {
293 					b->hits[j] = *p;
294 					b->hits[j].k = bns_depos(bns, bwt_sa(bwt, k), &is_rev);
295 					b->hits[j].l = 0;
296 					b->hits[j].is_rev = is_rev;
297 					if (is_rev) b->hits[j].k -= p->len - 1;
298 					++j;
299 				}
300 			} else if (p->G > 0) {
301 				b->hits[j] = *p;
302 				b->hits[j].k = bns_depos(bns, bwt_sa(bwt, p->k), &is_rev);
303 				b->hits[j].l = 0;
304 				b->hits[j].flag |= 1;
305 				b->hits[j].is_rev = is_rev;
306 				if (is_rev) b->hits[j].k -= p->len - 1;
307 				++j;
308 			}
309 		}
310 		free(old_hits);
311 	}
312 	for (i = j = 0; i < b->n; ++i) // squeeze out empty elements
313 		if (b->hits[i].G) b->hits[j++] = b->hits[i];
314 	b->n = j;
315 	ks_introsort(hitG, b->n, b->hits);
316 	for (i = 1; i < b->n; ++i) {
317 		bsw2hit_t *p = b->hits + i;
318 		for (j = 0; j < i; ++j) {
319 			bsw2hit_t *q = b->hits + j;
320 			int compatible = 1;
321 			if (p->is_rev != q->is_rev) continue; // hits from opposite strands are not duplicates
322 			if (p->l == 0 && q->l == 0) {
323 				int qol = (p->end < q->end? p->end : q->end) - (p->beg > q->beg? p->beg : q->beg); // length of query overlap
324 				if (qol < 0) qol = 0;
325 				if ((float)qol / (p->end - p->beg) > MASK_LEVEL || (float)qol / (q->end - q->beg) > MASK_LEVEL) {
326 					int64_t tol = (int64_t)(p->k + p->len < q->k + q->len? p->k + p->len : q->k + q->len)
327 						- (int64_t)(p->k > q->k? p->k : q->k); // length of target overlap
328 					if ((double)tol / p->len > MASK_LEVEL || (double)tol / q->len > MASK_LEVEL)
329 						compatible = 0;
330 				}
331 			}
332 			if (!compatible) {
333 				p->G = 0;
334 				if (q->G2 < p->G2) q->G2 = p->G2;
335 				break;
336 			}
337 		}
338 	}
339 	n = i;
340 	for (i = j = 0; i < n; ++i) {
341 		if (b->hits[i].G == 0) continue;
342 		if (i != j) b->hits[j++] = b->hits[i];
343 		else ++j;
344 	}
345 	b->n = j;
346 	return b->n;
347 }
348 
bsw2_resolve_query_overlaps(bwtsw2_t * b,float mask_level)349 int bsw2_resolve_query_overlaps(bwtsw2_t *b, float mask_level)
350 {
351 	int i, j, n;
352 	if (b->n == 0) return 0;
353 	ks_introsort(hitG, b->n, b->hits);
354 	{ // choose a random one
355 		int G0 = b->hits[0].G;
356 		for (i = 1; i < b->n; ++i)
357 			if (b->hits[i].G != G0) break;
358 		j = (int)(i * drand48());
359 		if (j) {
360 			bsw2hit_t tmp;
361 			tmp = b->hits[0]; b->hits[0] = b->hits[j]; b->hits[j] = tmp;
362 		}
363 	}
364 	for (i = 1; i < b->n; ++i) {
365 		bsw2hit_t *p = b->hits + i;
366 		int all_compatible = 1;
367 		if (p->G == 0) break;
368 		for (j = 0; j < i; ++j) {
369 			bsw2hit_t *q = b->hits + j;
370 			int64_t tol = 0;
371 			int qol, compatible = 0;
372 			float fol;
373 			if (q->G == 0) continue;
374 			qol = (p->end < q->end? p->end : q->end) - (p->beg > q->beg? p->beg : q->beg);
375 			if (qol < 0) qol = 0;
376 			if (p->l == 0 && q->l == 0) {
377 				tol = (int64_t)(p->k + p->len < q->k + q->len? p->k + p->len : q->k + q->len)
378 					- (p->k > q->k? p->k : q->k);
379 				if (tol < 0) tol = 0;
380 			}
381 			fol = (float)qol / (p->end - p->beg < q->end - q->beg? p->end - p->beg : q->end - q->beg);
382 			if (fol < mask_level || (tol > 0 && qol < p->end - p->beg && qol < q->end - q->beg)) compatible = 1;
383 			if (!compatible) {
384 				if (q->G2 < p->G) q->G2 = p->G;
385 				all_compatible = 0;
386 			}
387 		}
388 		if (!all_compatible) p->G = 0;
389 	}
390 	n = i;
391 	for (i = j = 0; i < n; ++i) {
392 		if (b->hits[i].G == 0) continue;
393 		if (i != j) b->hits[j++] = b->hits[i];
394 		else ++j;
395 	}
396 	b->n = j;
397 	return j;
398 }
399 /* --- END: processing partial hits --- */
400 
401 /* --- BEGIN: global mem pool --- */
bsw2_global_init()402 bsw2global_t *bsw2_global_init()
403 {
404 	bsw2global_t *pool;
405 	bsw2stack_t *stack;
406 	pool = calloc(1, sizeof(bsw2global_t));
407 	stack = calloc(1, sizeof(bsw2stack_t));
408 	stack->pool = (mempool_t*)calloc(1, sizeof(mempool_t));
409 	pool->stack = (void*)stack;
410 	return pool;
411 }
412 
bsw2_global_destroy(bsw2global_t * pool)413 void bsw2_global_destroy(bsw2global_t *pool)
414 {
415 	stack_destroy((bsw2stack_t*)pool->stack);
416 	free(pool->aln_mem);
417 	free(pool);
418 }
419 /* --- END: global mem pool --- */
420 
fill_cell(const bsw2opt_t * o,int match_score,bsw2cell_t * c[4])421 static inline int fill_cell(const bsw2opt_t *o, int match_score, bsw2cell_t *c[4])
422 {
423 	int G = c[3]? c[3]->G + match_score : MINUS_INF;
424 	if (c[1]) {
425 		c[0]->I = c[1]->I > c[1]->G - o->q? c[1]->I - o->r : c[1]->G - o->qr;
426 		if (c[0]->I > G) G = c[0]->I;
427 	} else c[0]->I = MINUS_INF;
428 	if (c[2]) {
429 		c[0]->D = c[2]->D > c[2]->G - o->q? c[2]->D - o->r : c[2]->G - o->qr;
430 		if (c[0]->D > G) G = c[0]->D;
431 	} else c[0]->D = MINUS_INF;
432 	return(c[0]->G = G);
433 }
434 
init_bwtsw2(const bwtl_t * target,const bwt_t * query,bsw2stack_t * s)435 static void init_bwtsw2(const bwtl_t *target, const bwt_t *query, bsw2stack_t *s)
436 {
437 	bsw2entry_t *u;
438 	bsw2cell_t *x;
439 
440 	u = mp_alloc(s->pool);
441 	u->tk = 0; u->tl = target->seq_len;
442 	x = push_array_p(u);
443 	*x = g_default_cell;
444 	x->G = 0; x->qk = 0; x->ql = query->seq_len;
445 	u->n++;
446 	stack_push0(s, u);
447 }
448 /* On return, ret[1] keeps not-so-repetitive hits (narrow SA hits); ret[0] keeps all hits (right?) */
bsw2_core(const bntseq_t * bns,const bsw2opt_t * opt,const bwtl_t * target,const bwt_t * query,bsw2global_t * pool)449 bwtsw2_t **bsw2_core(const bntseq_t *bns, const bsw2opt_t *opt, const bwtl_t *target, const bwt_t *query, bsw2global_t *pool)
450 {
451 	bsw2stack_t *stack = (bsw2stack_t*)pool->stack;
452 	bwtsw2_t *b, *b1, **b_ret;
453 	int i, j, score_mat[16], *heap, heap_size, n_tot = 0;
454 	struct rusage curr, last;
455 	khash_t(qintv) *rhash;
456 	khash_t(64) *chash;
457 
458 	// initialize connectivity hash (chash)
459 	chash = bsw2_connectivity(target);
460 	// calculate score matrix
461 	for (i = 0; i != 4; ++i)
462 		for (j = 0; j != 4; ++j)
463 			score_mat[i<<2|j] = (i == j)? opt->a : -opt->b;
464 	// initialize other variables
465 	rhash = kh_init(qintv);
466 	init_bwtsw2(target, query, stack);
467 	heap_size = opt->z;
468 	heap = calloc(heap_size, sizeof(int));
469 	// initialize the return struct
470 	b = (bwtsw2_t*)calloc(1, sizeof(bwtsw2_t));
471 	b->n = b->max = target->seq_len * 2;
472 	b->hits = calloc(b->max, sizeof(bsw2hit_t));
473 	b1 = (bwtsw2_t*)calloc(1, sizeof(bwtsw2_t));
474 	b_ret = calloc(2, sizeof(void*));
475 	b_ret[0] = b; b_ret[1] = b1;
476 	// initialize timer
477 	getrusage(0, &last);
478 	// the main loop: traversal of the DAG
479 	while (!stack_isempty(stack)) {
480 		int old_n, tj;
481 		bsw2entry_t *v;
482 		uint32_t tcntk[4], tcntl[4];
483 		bwtint_t k, l;
484 
485 		v = stack_pop(stack); old_n = v->n;
486 		n_tot += v->n;
487 
488 		for (i = 0; i < v->n; ++i) { // test max depth and band width
489 			bsw2cell_t *p = v->array + i;
490 			if (p->ql == 0) continue;
491 			if (p->tlen - (int)p->qlen > opt->bw || (int)p->qlen - p->tlen > opt->bw) {
492 				p->qk = p->ql = 0;
493 				if (p->ppos >= 0) v->array[p->ppos].cpos[p->pj] = -5;
494 			}
495 		}
496 
497 		// get Occ for the DAG
498 		bwtl_2occ4(target, v->tk - 1, v->tl, tcntk, tcntl);
499 		for (tj = 0; tj != 4; ++tj) { // descend to the children
500 			bwtint_t qcntk[4], qcntl[4];
501 			int qj, *curr_score_mat = score_mat + tj * 4;
502 			khiter_t iter;
503 			bsw2entry_t *u;
504 
505 			k = target->L2[tj] + tcntk[tj] + 1;
506 			l = target->L2[tj] + tcntl[tj];
507 			if (k > l) continue;
508 			// update counter
509 			iter = kh_get(64, chash, (uint64_t)k<<32 | l);
510 			--kh_value(chash, iter);
511 			// initialization
512 			u = mp_alloc(stack->pool);
513 			u->tk = k; u->tl = l;
514 			memset(heap, 0, sizeof(int) * opt->z);
515 			// loop through all the nodes in v
516 		    for (i = 0; i < v->n; ++i) {
517 				bsw2cell_t *p = v->array + i, *x, *c[4]; // c[0]=>current, c[1]=>I, c[2]=>D, c[3]=>G
518 				int is_added = 0;
519 				if (p->ql == 0) continue; // deleted node
520 				c[0] = x = push_array_p(u);
521 				x->G = MINUS_INF;
522 				p->upos = x->upos = -1;
523 				if (p->ppos >= 0) { // parent has been visited
524 					c[1] = (v->array[p->ppos].upos >= 0)? u->array + v->array[p->ppos].upos : 0;
525 					c[3] = v->array + p->ppos; c[2] = p;
526 					if (fill_cell(opt, curr_score_mat[p->pj], c) > 0) { // then update topology at p and x
527 						x->ppos = v->array[p->ppos].upos; // the parent pos in u
528 						p->upos = u->n++; // the current pos in u
529 						if (x->ppos >= 0) u->array[x->ppos].cpos[p->pj] = p->upos; // the child pos of its parent in u
530 						is_added = 1;
531 					}
532 				} else {
533 					x->D = p->D > p->G - opt->q? p->D - opt->r : p->G - opt->qr;
534 					if (x->D > 0) {
535 						x->G = x->D;
536 						x->I = MINUS_INF; x->ppos = -1;
537 						p->upos = u->n++;
538 						is_added = 1;
539 					}
540 				}
541 				if (is_added) { // x has been added to u->array. fill the remaining variables
542 					x->cpos[0] = x->cpos[1] = x->cpos[2] = x->cpos[3] = -1;
543 					x->pj = p->pj; x->qk = p->qk; x->ql = p->ql; x->qlen = p->qlen; x->tlen = p->tlen + 1;
544 					if (x->G > -heap[0]) {
545 						heap[0] = -x->G;
546 						ks_heapadjust(int, 0, heap_size, heap);
547 					}
548 				}
549 				if ((x->G > opt->qr && x->G >= -heap[0]) || i < old_n) { // good node in u, or in v
550 					if (p->cpos[0] == -1 || p->cpos[1] == -1 || p->cpos[2] == -1 || p->cpos[3] == -1) {
551 						bwt_2occ4(query, p->qk - 1, p->ql, qcntk, qcntl);
552 						for (qj = 0; qj != 4; ++qj) { // descend to the prefix trie
553 							if (p->cpos[qj] != -1) continue; // this node will be visited later
554 							k = query->L2[qj] + qcntk[qj] + 1;
555 							l = query->L2[qj] + qcntl[qj];
556 							if (k > l) { p->cpos[qj] = -2; continue; }
557 							x = push_array_p(v);
558 							p = v->array + i; // p may not point to the correct position after realloc
559 							x->G = x->I = x->D = MINUS_INF;
560 							x->qk = k; x->ql = l; x->pj = qj; x->qlen = p->qlen + 1; x->ppos = i; x->tlen = p->tlen;
561 							x->cpos[0] = x->cpos[1] = x->cpos[2] = x->cpos[3] = -1;
562 							p->cpos[qj] = v->n++;
563 						} // ~for(qj)
564 					} // ~if(p->cpos[])
565 				} // ~if
566 			} // ~for(i)
567 			if (u->n) save_hits(target, opt->t, b->hits, u);
568 			{ // push u to the stack (or to the pending array)
569 				uint32_t cnt, pos;
570 				cnt = (uint32_t)kh_value(chash, iter);
571 				pos = kh_value(chash, iter)>>32;
572 				if (pos) { // something in the pending array, then merge
573 					bsw2entry_t *w = kv_A(stack->pending, pos-1);
574 					if (u->n) {
575 						if (w->n < u->n) { // swap
576 							w = u; u = kv_A(stack->pending, pos-1); kv_A(stack->pending, pos-1) = w;
577 						}
578 						merge_entry(opt, w, u, b);
579 					}
580 					if (cnt == 0) { // move from pending to stack0
581 						remove_duplicate(w, rhash);
582 						save_narrow_hits(target, w, b1, opt->t, opt->is);
583 						cut_tail(w, opt->z, u);
584 						stack_push0(stack, w);
585 						kv_A(stack->pending, pos-1) = 0;
586 						--stack->n_pending;
587 					}
588 					mp_free(stack->pool, u);
589 				} else if (cnt) { // the first time
590 					if (u->n) { // push to the pending queue
591 						++stack->n_pending;
592 						kv_push(bsw2entry_p, stack->pending, u);
593 						kh_value(chash, iter) = (uint64_t)kv_size(stack->pending)<<32 | cnt;
594 					} else mp_free(stack->pool, u);
595 				} else { // cnt == 0, then push to the stack
596 					bsw2entry_t *w = mp_alloc(stack->pool);
597 					save_narrow_hits(target, u, b1, opt->t, opt->is);
598 					cut_tail(u, opt->z, w);
599 					mp_free(stack->pool, w);
600 					stack_push0(stack, u);
601 				}
602 			}
603 		} // ~for(tj)
604 		mp_free(stack->pool, v);
605 	} // while(top)
606 	getrusage(0, &curr);
607 	for (i = 0; i < 2; ++i)
608 		for (j = 0; j < b_ret[i]->n; ++j)
609 			b_ret[i]->hits[j].n_seeds = 0;
610 	bsw2_resolve_duphits(bns, query, b, opt->is);
611 	bsw2_resolve_duphits(bns, query, b1, opt->is);
612 	//fprintf(stderr, "stats: %.3lf sec; %d elems\n", time_elapse(&curr, &last), n_tot);
613 	// free
614 	free(heap);
615 	kh_destroy(qintv, rhash);
616 	kh_destroy(64, chash);
617 	stack->pending.n = stack->stack0.n = 0;
618 	return b_ret;
619 }
620