diff --git a/chain.c b/chain.c index 45801b7..7f3ce28 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 match_len, int max_dist, int bw, int max_skip, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km) +int mm_chain_dp(int max_dist, int bw, int max_skip, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km) { int32_t st = 0, i, j, k, *p, *f, *t, *v, n_u, n_v; uint64_t *u; @@ -34,7 +34,7 @@ int mm_chain_dp(int match_len, int max_dist, int bw, int max_skip, int min_sc, i // fill the score and backtrack arrays for (i = 0; i < n; ++i) { uint64_t ri = a[i].x; - int32_t qi = (int32_t)a[i].y; + int32_t qi = (int32_t)a[i].y, q_span = a[i].y>>32; int32_t max_f = -INT32_MAX, max_j = -1, n_skip = 0; while (st < i && ri - a[st].x > max_dist) ++st; for (j = i - 1; j >= st; --j) { @@ -48,12 +48,12 @@ int mm_chain_dp(int match_len, int max_dist, int bw, int max_skip, int min_sc, i } dd = dr > dq? dr - dq : dq - dr; if (dd > bw) continue; - sc = dq > match_len && dr > match_len? match_len : dq < dr? dq : dr; - sc = f[j] + sc - ilog2_32(dd); + sc = dq > q_span && dr > q_span? q_span : dq < dr? dq : dr; + sc = f[j] + sc - (dd? ilog2_32(dd) : 0); if (sc > max_f) max_f = sc, max_j = j; } if (max_j >= 0) f[i] = max_f, p[i] = max_j; - else f[i] = match_len, p[i] = -1; + else f[i] = q_span, p[i] = -1; } // find the ending positions of chains diff --git a/map.c b/map.c index 8f826ac..ea89c24 100644 --- a/map.c +++ b/map.c @@ -227,7 +227,7 @@ void mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b, uint3 kfree(b->km_fixed, m); //for (i = 0; i < n_a; ++i) printf("%c\t%s\t%d\t%d\n", "+-"[a[i].x>>63], mi->seq[a[i].x<<1>>33].name, (uint32_t)a[i].x>>1, (uint32_t)a[i].y); - n_u = mm_chain_dp(mi->k, opt->max_gap, opt->bw, opt->max_skip, opt->min_score, n_a, a, &u, b->km_fixed); + n_u = mm_chain_dp(opt->max_gap, opt->bw, opt->max_skip, opt->min_score, n_a, a, &u, b->km_fixed); printf("%s\t%d", qname, n_u); for (i = j = 0; i < n_u; ++i) { int n = (uint32_t)u[i]; diff --git a/mmpriv.h b/mmpriv.h index 13d64aa..82684b3 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -14,7 +14,7 @@ void radix_sort_128x(mm128_t *beg, mm128_t *end); void radix_sort_64(uint64_t *beg, uint64_t *end); uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk); -int mm_chain_dp(int match_len, int max_dist, int bw, int max_skip, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km); +int mm_chain_dp(int max_dist, int bw, int max_skip, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km); #ifdef __cplusplus }