diff --git a/chain.c b/chain.c index 1ad83e3..53e24b4 100644 --- a/chain.c +++ b/chain.c @@ -19,7 +19,7 @@ static inline int ilog2_32(uint32_t v) return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v]; } -int mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, uint64_t **_u, void *km) +mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km) { // TODO: make sure this works when n has more than 32 bits int32_t st = 0, k, *f, *p, *t, *v, n_u, n_v; int64_t i, j; @@ -148,6 +148,10 @@ int mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cn k += n; } memcpy(u, u2, n_u * 8); - kfree(km, b); kfree(km, w); kfree(km, u2); - return n_u; + + // write _a_ to _b_ and deallocate _a_ because _a_ is oversized, sometimes a lot + memcpy(b, a, n_v * sizeof(mm128_t)); + *n_u_ = n_u; + kfree(km, a); kfree(km, w); kfree(km, u2); + return b; } diff --git a/map.c b/map.c index 4778b07..262d1a2 100644 --- a/map.c +++ b/map.c @@ -252,7 +252,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k void mm_map_multi(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname) { - int i, j, n_u, max_gap_ref, rep_len, qlen_sum, n_regs0; + int i, j, max_gap_ref, rep_len, qlen_sum, n_regs0; int64_t n_a; uint64_t *u; mm128_t *a; @@ -275,23 +275,11 @@ void mm_map_multi(const mm_idx_t *mi, int n_segs, const int *qlens, const char * } max_gap_ref = opt->max_gap_ref >= 0? opt->max_gap_ref : opt->max_gap; - n_u = mm_chain_dp(max_gap_ref, opt->max_gap, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, !!(opt->flag&MM_F_SPLICE), n_segs, n_a, a, &u, b->km); - if (n_u > 0) { // shrink _a_ because after chaining, the size of _a_ may be much reduced - int64_t n_a0 = n_a; - for (j = 0, n_a = 0; j < n_u; ++j) - n_a += (int32_t)u[j]; - if (n_a < n_a0>>1) { - mm128_t *a0 = a; - a = (mm128_t*)kmalloc(b->km, n_a * sizeof(mm128_t)); - memcpy(a, a0, n_a * sizeof(mm128_t)); // anything beyond n_a is not used - kfree(b->km, a0); - } - } - regs0 = mm_gen_regs(b->km, qlen_sum, n_u, u, a); - n_regs0 = n_u; + a = mm_chain_dp(max_gap_ref, opt->max_gap, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, !!(opt->flag&MM_F_SPLICE), n_segs, n_a, a, &n_regs0, &u, b->km); + regs0 = mm_gen_regs(b->km, qlen_sum, n_regs0, u, a); if (mm_dbg_flag & MM_DBG_PRINT_SEED) - for (j = 0; j < n_u; ++j) + for (j = 0; j < n_regs0; ++j) for (i = regs0[j].as; i < regs0[j].as + regs0[j].cnt; ++i) fprintf(stderr, "CN\t%d\t%s\t%d\t%c\t%d\t%d\t%d\n", j, mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff), i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x)); diff --git a/mmpriv.h b/mmpriv.h index 0afb432..922066d 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -63,7 +63,7 @@ void mm_idxopt_init(mm_idxopt_t *opt); const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n); int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq); int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f); -int mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, uint64_t **_u, void *km); +mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km); mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a); mm_reg1_t *mm_gen_regs(void *km, int qlen, int n_u, uint64_t *u, mm128_t *a);