diff --git a/lchain.c b/lchain.c index 7df5cab..0bb1e73 100644 --- a/lchain.c +++ b/lchain.c @@ -149,7 +149,7 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int int is_cdna, int n_seg, 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 *f, *t, *v, n_u, n_v, mmax_f = 0, max_drop = bw; - int64_t *p, i, j, max_ii, st = 0, n_iter = 0; + int64_t *p, i, j, max_ii, st = 0; uint64_t *u; if (_u) *_u = 0, *n_u_ = 0; @@ -174,7 +174,6 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int for (j = i - 1; j >= st; --j) { int32_t sc; sc = comput_sc(&a[i], &a[j], max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg); - ++n_iter; if (sc == INT32_MIN) continue; sc += f[j]; if (sc > max_f) { @@ -204,6 +203,7 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int if (max_ii < 0 || (a[i].x - a[max_ii].x <= (int64_t)max_dist_x && f[max_ii] < f[i])) max_ii = i; if (mmax_f < max_f) mmax_f = max_f; + //fprintf(stderr, "X1\t%ld\t%ld:%d\t%ld\t%ld:%d\t%ld\t%ld\t%ld\n", (long)i, (long)(a[i].x>>32), (int32_t)a[i].x, (long)max_j, max_j<0?-1L:(long)(a[max_j].x>>32), max_j<0?-1:(int32_t)a[max_j].x, (long)max_f, (long)v[i], (long)mmax_f); } u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, max_drop, &n_u, &n_v); @@ -325,12 +325,11 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski krmq_interval(lc_elem, root_inner, &s, &lo, &hi); if (lo) { const lc_elem_t *q; - int32_t width, n_rmq_iter = 0; + int32_t width; krmq_itr_t(lc_elem) itr; krmq_itr_find(lc_elem, root_inner, lo, &itr); while ((q = krmq_at(&itr)) != 0) { if (q->y < (int32_t)a[i].y - max_dist_inner) break; - ++n_rmq_iter; j = q->i; sc = f[j] + comput_sc_simple(&a[i], &a[j], chn_pen_gap, chn_pen_skip, 0, &width); if (width <= bw) { diff --git a/main.c b/main.c index 4da93bc..5b6aacc 100644 --- a/main.c +++ b/main.c @@ -78,6 +78,8 @@ static ko_longopt_t long_options[] = { { "no-hash-name", ko_no_argument, 353 }, { "secondary-seq", ko_no_argument, 354 }, { "ds", ko_no_argument, 355 }, + { "rmq-inner", ko_required_argument, 356 }, + { "dbg-seed-freq", ko_no_argument, 501 }, { "help", ko_no_argument, 'h' }, { "max-intron-len", ko_required_argument, 'G' }, { "version", ko_no_argument, 'V' }, @@ -245,6 +247,8 @@ int main(int argc, char *argv[]) else if (c == 353) opt.flag |= MM_F_NO_HASH_NAME; // --no-hash-name else if (c == 354) opt.flag |= MM_F_SECONDARY_SEQ; // --secondary-seq else if (c == 355) opt.flag |= MM_F_OUT_DS; // --ds + else if (c == 356) opt.rmq_inner_dist = mm_parse_num(o.arg); // --rmq-inner + else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-freq else if (c == 330) { fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n"); } else if (c == 314) { // --frag diff --git a/minimap.h b/minimap.h index b944199..884d42d 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.27-r1193" +#define MM_VERSION "2.27-r1200-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 2f5034b..fd82021 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -14,6 +14,7 @@ #define MM_DBG_PRINT_SEED 0x4 #define MM_DBG_PRINT_ALN_SEQ 0x8 #define MM_DBG_PRINT_CHAIN 0x10 +#define MM_DBG_SEED_FREQ 0x20 #define MM_SEED_LONG_JOIN (1ULL<<40) #define MM_SEED_IGNORE (1ULL<<41) diff --git a/seed.c b/seed.c index 76a67ae..fc5c14e 100644 --- a/seed.c +++ b/seed.c @@ -112,7 +112,8 @@ mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int ma } for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < n_m0; ++i) { mm_seed_t *q = &m[i]; - //fprintf(stderr, "X\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt); + if (mm_dbg_flag & MM_DBG_SEED_FREQ) + fprintf(stderr, "SF\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt); if (q->flt) { int en = (q->q_pos >> 1) + 1, st = en - q->q_span; if (st > rep_en) {