fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-15 13:07:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

backup; DON'T USE!!!

Browse Source
This commit is contained in:
Heng Li
2024-03-20 10:21:10 -04:00
parent c881b14ba5
commit 27b254ed6f
3 changed files with 115 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

106
lchain.c
View File

@@ -365,3 +365,109 @@ 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);
}

9
map.c
View File

@@ -287,8 +287,13 @@ 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);
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);
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);
}
}
} 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;

4
mmpriv.h
View File

@@ -80,12 +80,12 @@ int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, ui
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, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a, int is_qstrand);
mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float gap_scale,
int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
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);