Mercurial > repos > siyuan > prada

comparison pyPRADA_1.2/tools/bwa-0.5.7-mh/bwtsw2_core.c @ 0:acc2ca1a3ba4

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