diff --git a/lchain.c b/lchain.c index 1404c0a..9c76517 100644 --- a/lchain.c +++ b/lchain.c @@ -365,109 +365,3 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski } return compact_a(km, n_u, u, n_v, v, a); } - -mm128_t *mg_lchain_dp_rmq(int dist_s, int bw_s, int bw_l, int max_skip, int max_iter, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip, - int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km) -{ - int32_t *f, *t, *v, n_u, n_v, mmax_f = 0; - int32_t dist_l = dist_s > bw_l? dist_s : bw_l; - int64_t *p, i, j, max_ii, st_s = 0, i0, st_l = 0; - uint64_t *u; - lc_elem_t *root = 0; - void *mem_mp = 0; - kmp_rmq_t *mp; - - if (_u) *_u = 0, *n_u_ = 0; - if (n == 0 || a == 0) { - kfree(km, a); - return 0; - } - p = Kmalloc(km, int64_t, n); - f = Kmalloc(km, int32_t, n); - v = Kmalloc(km, int32_t, n); - t = Kcalloc(km, int32_t, n); - mem_mp = km_init2(km, 0x10000); - mp = kmp_init_rmq(mem_mp); - - // fill the score and backtrack arrays - for (i = 0, max_ii = -1; i < n; ++i) { - int64_t max_j = -1, end_j; - int32_t max_f = a[i].y>>32&0xff, n_skip = 0; - lc_elem_t s, *q, lo, hi; - // add in-range anchors - if (i0 < i && a[i0].x != a[i].x) { - int64_t j; - for (j = i0; j < i; ++j) { - q = kmp_alloc_rmq(mp); - q->y = (int32_t)a[j].y, q->i = j, q->pri = -(f[j] + 0.25 * chn_pen_gap * ((int32_t)a[j].x + (int32_t)a[j].y)); - krmq_insert(lc_elem, &root, q, 0); - } - i0 = i; - } - // get rid of active chains out of range - while (st_l < i && (a[i].x>>32 != a[st_l].x>>32 || a[i].x > a[st_l].x + dist_l || krmq_size(head, root) > cap_rmq_size)) { - s.y = (int32_t)a[st_l].y, s.i = st_l; - if ((q = krmq_find(lc_elem, root, &s, 0)) != 0) { - q = krmq_erase(lc_elem, &root, q, 0); - kmp_free_rmq(mp, q); - } - ++st_l; - } - // RMQ - lo.i = INT32_MAX, lo.y = (int32_t)a[i].y - dist_l; - hi.i = 0, hi.y = (int32_t)a[i].y; - if ((q = krmq_rmq(lc_elem, root, &lo, &hi)) != 0) { - int32_t sc, exact, width; - int64_t j = q->i; - assert(q->y >= lo.y && q->y <= hi.y); - sc = f[j] + comput_sc_simple(&a[i], &a[j], chn_pen_gap, chn_pen_skip, &exact, &width); - if (width <= bw_l && sc > max_f) max_f = sc, max_j = j; - } - // determine st_s - while (st_s < i && (a[i].x>>32 != a[st_s].x>>32 || a[i].x > a[st_s].x + dist_s)) ++st_s; - if (i - st_s > max_iter) st_s = i - max_iter; - // core dp loop - for (j = i - 1; j >= st_s; --j) { - int32_t sc; - sc = comput_sc(&a[i], &a[j], dist_s, dist_s, bw_s, chn_pen_gap, chn_pen_skip, 0, 1); - if (sc == INT32_MIN) continue; - sc += f[j]; - if (sc > max_f) { - max_f = sc, max_j = j; - if (n_skip > 0) --n_skip; - } else if (t[j] == (int32_t)i) { - if (++n_skip > max_skip) - break; - } - if (p[j] >= 0) t[p[j]] = i; - } - end_j = j; - if (max_ii < 0 || a[i].x - a[max_ii].x > (int64_t)dist_s) { - int32_t max = INT32_MIN; - max_ii = -1; - for (j = i - 1; j >= st_s; --j) - if (max < f[j]) max = f[j], max_ii = j; - } - if (max_ii >= 0 && max_ii < end_j) { - int32_t tmp; - tmp = comput_sc(&a[i], &a[max_ii], dist_s, dist_s, bw_s, chn_pen_gap, chn_pen_skip, 0, 1); - if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) - max_f = tmp + f[max_ii], max_j = max_ii; - } - f[i] = max_f, p[i] = max_j; - v[i] = max_j >= 0 && v[max_j] > max_f? v[max_j] : max_f; // v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak - if (max_ii < 0 || (a[i].x - a[max_ii].x <= (int64_t)dist_s && f[max_ii] < f[i])) - max_ii = i; - if (mmax_f < max_f) mmax_f = max_f; - } - km_destroy(mem_mp); - - u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, bw_l, &n_u, &n_v); - *n_u_ = n_u, *_u = u; // NB: note that u[] may not be sorted by score here - kfree(km, p); kfree(km, f); kfree(km, t); - if (n_u == 0) { - kfree(km, a); kfree(km, v); - return 0; - } - return compact_a(km, n_u, u, n_v, v, a); -} diff --git a/map.c b/map.c index 73ca7ec..038888f 100644 --- a/map.c +++ b/map.c @@ -287,13 +287,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** for (i = 0, n_a = 0; i < n_regs0; ++i) n_a += (int32_t)u[i]; kfree(b->km, u); radix_sort_128x(a, a + n_a); - if (opt->flag & MM_F_RMQ) { - a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score, - chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km); - } else { - a = mg_lchain_dp_rmq(opt->max_gap, opt->bw, opt->bw_long, opt->max_chain_skip, opt->max_chain_iter, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score, - chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km); - } + a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score, + chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km); } } else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) { // re-chain, mostly for short reads int rechain = 0; diff --git a/minimap.h b/minimap.h index 1ed4923..54f8596 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.27-r1203-dirty" +#define MM_VERSION "2.27-r1205-dirty" #define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit #define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name diff --git a/mmpriv.h b/mmpriv.h index b633789..6de9d0c 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -84,8 +84,6 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km); mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km); -mm128_t *mg_lchain_dp_rmq(int dist_s, int bw_s, int bw_l, int max_skip, int max_iter, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip, - int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km); void mm_mark_alt(const mm_idx_t *mi, int n, mm_reg1_t *r); void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a, int is_qstrand); diff --git a/options.c b/options.c index 1a033e1..cd2990d 100644 --- a/options.c +++ b/options.c @@ -118,6 +118,9 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo) io->flag = 0, io->k = 25, io->w = 51; mo->min_mid_occ = 50, mo->max_mid_occ = 500; mo->rmq_inner_dist = 5000; + mo->occ_dist = 200; + mo->best_n = 100; + mo->chain_gap_scale = 5.0f; } else if (strcmp(preset, "map-iclr-prerender") == 0) { io->flag = 0, io->k = 15; mo->b = 6, mo->transition = 1;