dp opt

CommandLines.cpp

@@ -40,6 +40,9 @@ static ko_longopt_t long_options[] = {
     { "n-weight",      ko_required_argument, 326 },
     { "l-msjoin",      ko_required_argument, 327 },
     { "purge-max",     ko_required_argument, 328 },
+    { "fast",     ko_no_argument, 329 },
+    { "dp-er",     ko_required_argument, 330},
+    { "max-kocc",     ko_required_argument, 331},
 	{ 0, 0, 0 }
 };
 
@@ -207,6 +210,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->is_alt = 0;
     asm_opt->misjoin_len = 500000;
     asm_opt->scffold = 0;
+    asm_opt->dp_min_len = 2000;
+    asm_opt->dp_e = 0.0025;
 }
 
 void destory_enzyme(enzyme* f)
@@ -677,6 +682,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 326) asm_opt->n_weight = atoi(opt.arg);
         else if (c == 327) asm_opt->misjoin_len = atol(opt.arg);
         else if (c == 328) asm_opt->pur_global_coverage = atoi(opt.arg);
+        else if (c == 329) asm_opt->flag |= HA_F_FAST;
+        else if (c == 330) asm_opt->dp_e = atof(opt.arg);
+        else if (c == 331) asm_opt->max_kmer_cnt = atol(opt.arg);
         else if (c == 'l')
         {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
             asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);

CommandLines.h

@@ -4,7 +4,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.15.5-r352"
+#define HA_VERSION "0.15.5-r358"
 
 #define VERBOSE 0
 
@@ -20,6 +20,7 @@
 #define HA_F_BAN_ASSEMBLY    0x200
 #define HA_F_HIGH_HET        0x400
 #define HA_F_PARTITION       0x800
+#define HA_F_FAST            0x1000
 
 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 
@@ -109,6 +110,8 @@ typedef struct {
     uint32_t is_alt;
     uint64_t misjoin_len;
     uint64_t scffold;
+    int32_t dp_min_len;
+    float dp_e;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;

Hash_Table.cpp

@@ -778,70 +778,6 @@ skip_dp:
     return chainLen;
 }
 
-void calculate_overlap_region_by_chaining_back(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
-                                          uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end)
-{
-    overlap_region tmp_region;
-    long long i = 0;
-    uint64_t current_ID;
-    uint64_t current_stand;
-
-    if (candidates->length == 0)
-    {
-        return;
-    }
-
-    long long sub_region_beg;
-    long long sub_region_end;
-
-    init_fake_cigar(&(tmp_region.f_cigar));
-
-    i = 0;
-    while (i < candidates->length)
-    {
-        current_ID = candidates->list[i].readID;
-        current_stand = candidates->list[i].strand;
-
-        ///reference read
-        tmp_region.x_id = readID;
-        tmp_region.x_pos_strand = current_stand;
-        ///query read
-        tmp_region.y_id = current_ID;
-        ///here the strand of query is always 0
-        tmp_region.y_pos_strand = 0;  
-
-        sub_region_beg = i;
-        sub_region_end = i;
-        i++;
-
-        while (i < candidates->length 
-        && 
-        current_ID == candidates->list[i].readID
-        &&
-        current_stand == candidates->list[i].strand)
-        {
-            sub_region_end = i;
-            i++;
-        }
-
-        if (tmp_region.x_id == tmp_region.y_id)
-        {
-            continue;
-        }
-
-		chain_DP(candidates->list + sub_region_beg,
-				sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold,
-				25, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id));
-
-        ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1)
-        if (tmp_region.x_id != tmp_region.y_id)
-        {
-            append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF, add_beg_end);
-        }
-    }
-
-    destory_fake_cigar(&(tmp_region.f_cigar));
-}
 
 
 void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,

Overlaps.cpp

@@ -16473,7 +16473,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
             ///need this line
             if(nsg->seq[uId].c == ALTER_LABLE) continue;
             pre = uId; 
-            pre = pre | (0x100000000);
+            pre = pre | (uint64_t)(0x100000000);
             kv_push(uint64_t, u_vecs->a, pre);
         }
     }
@@ -16483,7 +16483,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
     {
         if(u_vecs->a.a[k] == u_vecs->a.a[k-1])
         {
-            u_vecs->a.a[m] += (0x100000000);
+            u_vecs->a.a[m] += (uint64_t)(0x100000000);
         }
         else
         {

anchor.cpp

@@ -73,7 +73,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	rlen = Get_READ_LENGTH(R_INF, rid); // read length
 
 	// get the list of anchors
-	ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct);
+	ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, NULL);
 	// minimizer of queried read
 	if (ab->mz.m > ab->old_mz_m) {
 		ab->old_mz_m = ab->mz.m;
@@ -145,7 +145,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		}
 	}
 	#endif
-
+	
 	if ((int)overlap_list->length > max_n_chain) {
 		int32_t w, n[4], s[4];
 		n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
@@ -157,9 +157,12 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 			if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
 		}
 		if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
+			// n[0] = n[1] = n[2] = n[3] = 0;
 			for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
 				overlap_region *r = &overlap_list->list[i];
 				w = ha_ov_type(r, rlen);
+				// ++n[w];
+				// if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) {
 				if (r->shared_seed >= s[w]) {
 					if ((uint32_t)k != i) {
 						overlap_region t;

htab.cpp

@@ -48,7 +48,7 @@ typedef struct {
 	int32_t pre;
 	int32_t n_thread;
 	int64_t chunk_size;
-	int adaLen;
+	int adaLen, min_rcnt;
 } yak_copt_t;
 
 void yak_copt_init(yak_copt_t *o)
@@ -286,6 +286,8 @@ static void ha_ct_shrink(ha_ct_t *h, int min, int max, int n_thread)
 	kt_for(n_thread, worker_ct_shrink, &a, 1<<h->pre);
 	for (i = 0, h->tot = 0; i < 1<<h->pre; ++i)
 		h->tot += kh_size(h->h[i].h);
+	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
+			yak_realtime(), yak_cpu_usage(), (long)h->tot);
 }
 
 /***********************
@@ -313,6 +315,42 @@ typedef struct {
 	ha_pt_t *pt;
 } pt_gen_aux_t;
 
+
+static void worker_pt_shrink(void *data, long i, int tid) // callback for kt_for()
+{
+    ha_pt_t *h = (ha_pt_t*)data;
+	ha_pt1_t *b = &h->h[i];
+    yak_pt_t *f = NULL;
+    khint_t k;
+    f = yak_pt_init();
+    for (k = 0, b->n = 0; k < kh_end(b->h); ++k) {
+        if (kh_exist(b->h, k)) {
+			if(kh_val(b->h, k) <= 0) continue;
+			int absent; khint_t l;
+			l = yak_pt_put(f, (kh_key(b->h, k) >> h->pre) << YAK_COUNTER_BITS, &absent);
+			kh_val(f, l) = b->n; 
+			b->n += kh_key(b->h, k) & YAK_MAX_COUNT;
+        }
+    }
+    yak_pt_destroy(b->h);
+    h->h[i].h = f;
+	CALLOC(b->a, b->n);///need fix
+}
+
+static uint64_t ha_pt_shrink(ha_pt_t *h, int n_thread)
+{
+    int i;
+	uint64_t occ;
+    ///still start 4096 threads
+    kt_for(n_thread, worker_pt_shrink, h, 1<<h->pre);
+    for (i = 0, occ = 0, h->tot = 0; i < 1<<h->pre; ++i)
+	{
+		h->tot += kh_size(h->h[i].h);
+		occ += h->h[i].n;
+	}
+	return occ;
+}
+
 static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for()
 {
 	pt_gen_aux_t *a = (pt_gen_aux_t*)data;
@@ -353,6 +391,44 @@ ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread)
 	return pt;
 }
 
+static void worker_pt_gen_count(void *data, long i, int tid) // callback for kt_for()
+{
+	pt_gen_aux_t *a = (pt_gen_aux_t*)data;
+	ha_pt1_t *b = &a->pt->h[i];
+	yak_ct_t *g = a->ct->h[i].h;
+	khint_t k;
+	for (k = 0, b->n = 0; k != kh_end(g); ++k) {
+		if (kh_exist(g, k)) {
+			int absent;
+			khint_t l;
+			l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent);
+			kh_val(b->h, l) = 0; kh_key(b->h, l) |= kh_key(g, k) & YAK_MAX_COUNT;
+		}
+	}
+	yak_ct_destroy(g);
+	a->ct->h[i].h = 0;
+}
+
+ha_pt_t *ha_pt_gen_count(ha_ct_t *ct, int n_thread)
+{
+	pt_gen_aux_t a;
+	int i;
+	ha_pt_t *pt;
+	ha_ct_destroy_bf(ct);
+	CALLOC(pt, 1);
+	pt->k = ct->k, pt->pre = ct->pre, pt->tot = ct->tot;
+	CALLOC(pt->h, 1<<pt->pre);
+	for (i = 0; i < 1<<pt->pre; ++i) {
+		pt->h[i].h = yak_pt_init();
+		yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h));
+	}
+	a.ct = ct, a.pt = pt;
+	kt_for(n_thread, worker_pt_gen_count, &a, 1<<pt->pre);
+	free(ct->h); free(ct);
+	return pt;
+}
+
+
 int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a)
 {
 	int j, mask = (1<<h->pre) - 1, n_ins = 0;
@@ -377,6 +453,26 @@ int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a)
 	}
 	return n_ins;
 }
+
+
+int ha_pt_cnt_insert_list(ha_pt_t *h, int n, const uint64_t *a)
+{
+	int j, mask = (1<<h->pre) - 1, n_ins = 0;
+	ha_pt1_t *g;
+	if (n == 0) return 0;
+	g = &h->h[a[0]&mask];
+	for (j = 0; j < n; ++j) {
+		uint64_t x = a[j] >> h->pre;
+		khint_t k;
+		assert((a[j]&mask) == (a[0]&mask));
+		k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS);
+		if (k == kh_end(g->h)) continue; // TODO: understand why we sometimes come here
+		++kh_val(g->h, k);
+		++n_ins;
+	}
+	// fprintf(stderr, "n: %d, n_ins: %d\n", n, n_ins);
+	return n_ins;
+}
 /*
 static void worker_pt_sort(void *data, long i, int tid)
 {
@@ -420,6 +516,15 @@ const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n)
 	return &g->a[kh_val(g->h, k)];
 }
 
+const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash)
+{
+	khint_t k;
+	const ha_pt1_t *g = &h->h[hash & ((1ULL<<h->pre) - 1)];
+	k = yak_pt_get(g->h, hash >> h->pre << YAK_COUNTER_BITS);
+	if (k == kh_end(g->h)) return 0;
+	return kh_key(g->h, k) & YAK_MAX_COUNT;
+}
+
 /**********************************
  * Buffer for counting all k-mers *
  **********************************/
@@ -513,6 +618,7 @@ KSEQ_INIT(gzFile, gzread)
 #define HAF_CREATE_NEW   0x20
 #define HAF_SKIP_READ    0x40
 #define HAF_UG_READ      0x80
+#define HAF_COUNT_REFINE 0x100
 
 typedef struct { // global data structure for kt_pipeline()
 	const yak_copt_t *opt;
@@ -538,6 +644,7 @@ typedef struct { // data structure for each step in kt_pipeline()
 	ha_mz1_v *mz_buf;
 	ha_mz1_v *mz;
 	ch_buf_t *buf;
+	st_mt_t *mt;
 } st_data_t;
 
 static void worker_for_insert(void *data, long i, int tid) // callback for kt_for()
@@ -545,9 +652,20 @@ static void worker_for_insert(void *data, long i, int tid) // callback for kt_fo
 	st_data_t *s = (st_data_t*)data;
 	ch_buf_t *b = &s->buf[i];
 	if (s->p->pt)
-		b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b);
+	{
+		if(s->p->flag&HAF_COUNT_REFINE)
+		{
+			b->n_ins += ha_pt_cnt_insert_list(s->p->pt, b->n, b->a);
+		}
+		else
+		{
+			b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b);
+		}
+	}
 	else///for 0-th count, go into here
+	{
 		b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);
+	}		
 }
 
 static void worker_for_mz(void *data, long i, int tid)
@@ -556,8 +674,8 @@ static void worker_for_mz(void *data, long i, int tid)
 	///get the corresponding minimzer vector of this read
 	ha_mz1_v *b = &s->mz_buf[tid];
 	s->mz_buf[tid].n = 0;
-	///s->p->opt->w = 51, s->p->opt->k
-	ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0);
+	ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0, 
+	(s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, s->mt?&(s->mt[tid]):NULL);
 	s->mz[i].n = s->mz[i].m = b->n;
 	MALLOC(s->mz[i].a, b->n);
 	memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t));
@@ -668,7 +786,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		for (i = 0; i < n_pre; ++i) {
 			s->buf[i].m = m;
 			///for 0-th counting, p->pt = NULL
-			if (p->pt) MALLOC(s->buf[i].b, m);
+			if (p->pt && !(p->flag&HAF_COUNT_REFINE)) MALLOC(s->buf[i].b, m);
 			else MALLOC(s->buf[i].a, m);
 		}
 		// fill the buffer
@@ -689,13 +807,20 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 			// s->mz && s->mz_buf are lists of minimzer vectors
 			CALLOC(s->mz, s->n_seq);
 			CALLOC(s->mz_buf, p->opt->n_thread);
+			s->mt = NULL;
+			if(s->p->pt&&(s->p->flag&HAF_COUNT_REFINE)) CALLOC(s->mt, p->opt->n_thread);						
 			///calculate minimzers for each read, each read corresponds to one thread
 			kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
+
 			for (i = 0; i < p->opt->n_thread; ++i)
+			{
+				if(s->mt) free(s->mt[i].a);
 				free(s->mz_buf[i].a);
+			}
+			if(s->mt) free(s->mt);
 			free(s->mz_buf);
 			// insert minimizers
-			if (p->pt) {///insert whole minimizer
+			if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {///insert whole minimizer
 				for (i = 0; i < s->n_seq; ++i)
 					for (j = 0; j < s->mz[i].n; ++j)
 						pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
@@ -724,7 +849,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		///n_ins is number of distinct k-mers
 		for (i = 0; i < n; ++i) {
 			n_ins += s->buf[i].n_ins;
-			if (p->pt) free(s->buf[i].b);
+			if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b);
 			else free(s->buf[i].a);
 		}
 		if (p->ct) p->ct->tot += n_ins;
@@ -833,7 +958,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 	return pl.ct;
 }
 
-ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter)
+ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq)
 {
 	int i;
 	int64_t n_seq = 0;
@@ -858,6 +983,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
 	opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
 	opt.n_thread = asm_opt->thread_num;
 	opt.adaLen = (keep_adapter? asm_opt->adapterLen : 0);
+	opt.min_rcnt = (low_freq?*low_freq:-1);
 	///asm_opt->num_reads is the number of fastq files
 	for (i = 0; i < asm_opt->num_reads; ++i)
 		h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq);
@@ -948,7 +1074,7 @@ void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int hap_n)
 	int64_t cnt[YAK_N_COUNTS];
 	int cutoff = hap_n + 1;
 	ha_ct_t *h;
-	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, NULL, NULL, NULL, us, 0);
+	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, NULL, NULL, NULL, us, 0, NULL);
 
 	ha_ct_hist(h, cnt, asm_opt->thread_num);
 	print_hist_lines(YAK_N_COUNTS, 1, cnt);
@@ -969,7 +1095,7 @@ ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_
 	ha_ct_t *ct;
 	ha_pt_t *pt;
 	///HAF_COUNT_EXACT: no bf
-	ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, NULL, flt_tab, NULL, us, 0);
+	ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, NULL, flt_tab, NULL, us, 0, NULL);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
 			yak_realtime(), yak_cpu_usage(), (long)ct->tot);
 
@@ -987,7 +1113,7 @@ ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_
 		for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;
 	}
 	pt = ha_pt_gen(ct, asm_opt->thread_num);
-	ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, pt, flt_tab, NULL, us, 0);
+	ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, pt, flt_tab, NULL, us, 0, NULL);
 	assert((uint64_t)tot_cnt == pt->tot_pos);
 	//ha_pt_sort(pt, asm_opt->thread_num);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
@@ -1002,7 +1128,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
 	int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
 	if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
 	ha_ct_t *h;
-	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs, NULL, 1);
+	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs, NULL, 1, NULL);
 	if((asm_opt->flag & HA_F_VERBOSE_GFA))
 	{
 		write_ct_index((void*)h, asm_opt->output_file_name);
@@ -1030,6 +1156,90 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
 	return (void*)flt_tab;
 }
 
+void *ha_ft_gen_worse(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode)
+{
+	yak_ft_t *flt_tab;
+	int64_t cnt[YAK_N_COUNTS];
+	int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
+	if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
+	ha_ct_t *h;
+	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs, NULL, 1, NULL);
+	if((asm_opt->flag & HA_F_VERBOSE_GFA))
+	{
+		write_ct_index((void*)h, asm_opt->output_file_name);
+		// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
+		// debug_ct_index((void*)h, ha_ct_table);
+		// debug_ct_index(ha_ct_table, (void*)h);
+		// ha_ct_destroy((ha_ct_t *)ha_ct_table);
+	}
+	
+	if(!(ex_flag & HAF_SKIP_READ))
+	{
+		ha_ct_hist(h, cnt, asm_opt->thread_num);
+		peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
+		if (hom_cov) *hom_cov = peak_hom;
+		if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+	}
+	ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
+	flt_tab = gen_hh(h, 1/**asm_opt->max_kmer_cnt**/);
+	ha_ct_destroy(h);
+	fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__,
+			yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff);
+	return (void*)flt_tab;
+}
+
+ha_pt_t *ha_pt_gen_dp(const hifiasm_opt_t *asm_opt, ha_ct_t *ct, int flag, int n_thread, const void *flt_tab, All_reads *rs, int peak_hom, int peak_het)
+{
+	int low_freq = mz_low_b(peak_hom, peak_het);
+	ha_pt_t *pt = ha_pt_gen_count(ct, n_thread); ///key = cnt, val = 0
+	ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, pt, flt_tab, rs, NULL, 1, &low_freq);
+	uint64_t occ = ha_pt_shrink(pt, n_thread);
+	if(flag&HAF_RS_WRITE_LEN) flag -= HAF_RS_WRITE_LEN;
+	if(flag&HAF_RS_WRITE_SEQ) flag -= HAF_RS_WRITE_SEQ;
+	flag |= HAF_RS_READ; pt->tot_pos = 0;
+	ha_count(asm_opt, HAF_COUNT_EXACT|flag, pt, flt_tab, rs, NULL, 1, NULL);
+	// fprintf(stderr, "[M::%s::] counted %lu distinct minimizer k-mers\n", __func__, pt->tot);
+	// fprintf(stderr, "[M::%s::] collected %lu minimizers\n\n\n", __func__, pt->tot_pos);
+	assert(occ == pt->tot_pos);
+	return pt;
+}
+
+ha_pt_t *ha_pt_gen_worst(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov)
+{
+	int64_t cnt[YAK_N_COUNTS];
+	int peak_hom, peak_het, extra_flag1, extra_flag2;
+	ha_ct_t *ct;
+	ha_pt_t *pt;
+	if (read_from_store) {///if reads have already been read
+		extra_flag1 = extra_flag2 = HAF_RS_READ;
+	} else if (rs->total_reads == 0) {///if reads & length have not been scanned
+		extra_flag1 = HAF_RS_WRITE_LEN;
+		extra_flag2 = HAF_RS_WRITE_SEQ;
+	} else {///if length has been loaded but reads have not
+		extra_flag1 = HAF_RS_WRITE_SEQ;
+		extra_flag2 = HAF_RS_READ;
+	}
+	if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
+
+	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs, NULL, 1, NULL);
+	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
+			yak_realtime(), yak_cpu_usage(), (long)ct->tot);
+	ha_ct_hist(ct, cnt, asm_opt->thread_num);
+	fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
+	peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
+	if (hom_cov) *hom_cov = peak_hom;
+	if (het_cov) *het_cov = peak_het;
+	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+
+	///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time
+	ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
+	pt = ha_pt_gen_dp(asm_opt, ct, HAF_COUNT_EXACT|extra_flag2, asm_opt->thread_num, flt_tab, rs, peak_hom, peak_het);
+	//ha_pt_sort(pt, asm_opt->thread_num);
+	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
+			yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
+	return pt;
+}
+
 ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov)
 {
 	int64_t cnt[YAK_N_COUNTS], tot_cnt;
@@ -1047,7 +1257,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	}
 	if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
 
-	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs, NULL, 1);
+	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs, NULL, 1, NULL);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
 			yak_realtime(), yak_cpu_usage(), (long)ct->tot);
 	ha_ct_hist(ct, cnt, asm_opt->thread_num);
@@ -1069,12 +1279,21 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 		ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
 		for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;
 	}
-	pt = ha_pt_gen(ct, asm_opt->thread_num);
-	ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs, NULL, 1);
-	assert((uint64_t)tot_cnt == pt->tot_pos);
+	if(!(asm_opt->flag & HA_F_FAST))
+	{
+		fprintf(stderr, "[M::%s::] counting in normal mode\n", __func__);
+		pt = ha_pt_gen(ct, asm_opt->thread_num);
+		ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs, NULL, 1, NULL);
+		assert((uint64_t)tot_cnt == pt->tot_pos);
+	}
+	else
+	{
+		fprintf(stderr, "[M::%s::] counting in fast mode\n", __func__);
+		pt = ha_pt_gen_dp(asm_opt, ct, HAF_COUNT_EXACT|extra_flag2, asm_opt->thread_num, flt_tab, rs, peak_hom, peak_het);
+	}
 	//ha_pt_sort(pt, asm_opt->thread_num);
-	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
-			yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
+	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions, counted %ld distinct minimizer k-mers\n", __func__,
+			yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos, (long)pt->tot);
 	return pt;
 }
 
@@ -1084,7 +1303,7 @@ int query_ct_index(void* ct_idx, uint64_t hash)
 	khint_t k;
 	k = yak_ct_get(g->h, hash);
 	if (k == kh_end(g->h)) return 0;
-	return kh_key(g->h, k)&YAK_MAX_COUNT;
+	return ((kh_key(g->h, k)&YAK_MAX_COUNT)==YAK_MAX_COUNT)?-1:(kh_key(g->h, k)&YAK_MAX_COUNT);
 }
 
 int write_ct_index(void *i_ct_idx, char* file_name)

htab.h

@@ -5,6 +5,11 @@
 #include "Process_Read.h"
 #include "CommandLines.h"
 
+typedef struct {
+	int n, m;
+	uint64_t *a;
+} st_mt_t;
+
 typedef struct {
 	uint64_t x; ///x is the hash key
 	///rid is the read id, pos is the end pos of this minimizer, rev is the direction
@@ -75,6 +80,7 @@ ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_
 ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov);
 void ha_pt_destroy(ha_pt_t *h);
 const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
+const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash);
 
 int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
 int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
@@ -95,11 +101,20 @@ double yak_cpu_usage(void);
 
 void ha_triobin(const hifiasm_opt_t *opt);
 
-void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct);
+void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt);
 int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het);
 void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt);
 void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);
 
+inline int mz_low_b(int peak_hom, int peak_het)
+{
+	int low_freq = 2;
+	if(peak_het > 0) low_freq = peak_het/2;
+	else if(peak_hom > 0) low_freq = peak_hom/4;
+	if(low_freq < 2) low_freq = 2;
+	return low_freq;
+}
+
 static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
 {
 	key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1;

sketch.cpp

@@ -5,10 +5,11 @@
 #include "kvec.h"
 #include "htab.h"
 #include "ksort.h"
+#include "Correct.h"
 
 #define MAX_HIGH_OCC     8   // TODO: don't hard code if we need to tune this parameter
 #define MAX_MAX_HIGH_OCC 16
-
+#define GMC(a, x,y,xn) ((a)[(x)*(xn)+(y)])
 
 static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b)
 {
@@ -18,9 +19,219 @@ static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b)
 #define mz_lt(a, b) (mzcmp(&(a), &(b)) < 0)
 KSORT_INIT(mz, ha_mz1_t, mz_lt)
 
-static void select_mz(ha_mz1_v *p, int len, int sample_dist)
+void debug_refine(ha_mz1_t *ma, uint64_t *mmt, int32_t sn, int32_t n, int32_t m, int32_t end)
+{
+	uint64_t ks = end;
+	int64_t t = 0, i, k, sp = -1, ep = -1, ovlp, tot = mmt[end]&0xffffffff, nt = 0;;
+	while (ks != 0xffffffff)
+	{
+		i = ks/m; k = ks%m;
+		ks = mmt[ks]>>32;
+		if(ks == 0xffffffff || (int32_t)(ks/m) == (i-1))
+		{
+			t++;
+			ovlp = ((MIN(ep, (int64_t)ma[k].pos) >= MAX(sp, (int64_t)(ma[k].pos+1-ma[k].span)))? 
+						MIN(ep, (int64_t)ma[k].pos) - MAX(sp, (int64_t)(ma[k].pos+1-ma[k].span)) + 1:0);
+			if(ovlp != 0) fprintf(stderr, "ERROR-OVLP\n");
+			if(sp == -1 || sp > (ma[k].pos+1-ma[k].span)) sp = ma[k].pos+1-ma[k].span;
+			if(ep == -1 || ep < ma[k].pos) ep = ma[k].pos;
+			nt += (ma[k].rid);
+		}
+	}
+	if(t != sn) fprintf(stderr, "ERROR-TN, t: %ld, sn: %d\n", t, sn);
+	if(nt != tot) fprintf(stderr, "ERROR-TOT, nt: %ld, tot: %ld\n", nt, tot);
+}
+
+void refine_select(ha_mz1_v *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, 
+int32_t *rsi, int32_t *rei)
+{
+	int32_t n = sn, m = eidx + 1 - sidx, i, k, t, mk=-1;
+	uint64_t ix, kx, ks;
+	kv_resize(uint64_t, *mm, mm->n+n*m);
+	ha_mz1_t *ma = mz->a + sidx;
+	uint64_t *mmt = mm->a + mm->n;
+	// fprintf(stderr, "[M::%s::] ==> +n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);
+
+	for (i = 0; i < n; i++) ///how many selected minimizers 
+	{
+		for (k = 0, mk = -1; k < m; k++) ///how many minimizers in total
+		{
+			if((int32_t)(ma[k].rid)<min_freq) continue;
+			
+			ks = ma[k].pos + 1 - ma[k].span; t = -1;
+			if(i > 0)
+			{
+				for (t = k-1; t >= 0 && (ma[t].pos >= ks||(int32_t)(ma[t].rid)<min_freq); t--);
+			} 
+			
+			ix = (i <= 0?0:(t<0?0xffffffff:(GMC(mmt, i-1,t,m)&0xffffffff)));
+			if(ix < 0xffffffff) ix += (ma[k].rid);
+			kx = (mk < 0?0xffffffff:(GMC(mmt, i, mk,m)&0xffffffff));
+
+			
+
+			ks = MIN(ix, kx);
+			// fprintf(stderr, "ks: %lu, i: %d (n-%d), k: %d (m-%d), ix: %lu, kx: %lu, t: %d, mk: %d\n", 
+			// ks, i, n, k, m, ix, kx, t, mk);
+			if((ks&0xffffffff) == 0xffffffff) ks |= ((uint64_t)0xffffffff)<<32;
+			else if(ks == ix) ks |= (uint64_t)(i>0?(i-1)*m+t:0xffffffff)<<32;
+			else if(ks == kx) ks |= (uint64_t)(mk>=0?i*m+mk:0xffffffff)<<32;
+			
+			GMC(mmt, i,k,m) = ks;	
+			mk = k;
+		}
+	}
+	// fprintf(stderr, "[M::%s::] ==> ++n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);
+
+	ks = (n-1)*m + mk; ix = (uint64_t)-1; kx = 0;
+	while (ks != 0xffffffff)
+	{
+		i = ks/m; k = ks%m;
+		ks = mmt[ks]>>32;
+		// fprintf(stderr, "i: %d, k: %d, ks: %lu\n", i, k, ks);
+		if(ks == 0xffffffff || (int32_t)(ks/m) == (i-1))
+		{
+			mm->a[sidx+k] = 1;
+			ix = MIN((uint64_t)k, ix); kx = MAX((uint64_t)k, kx);
+		}
+	}
+	///debug
+	// debug_refine(ma, mmt, sn, n, m, (n-1)*m + mk);
+
+	if(rsi) (*rsi) = ix + sidx;
+	if(rei) (*rei) = kx + sidx;	
+}
+
+void refine_sketch(ha_mz1_v *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, float er, int32_t min_freq, st_mt_t *mt)
+{
+	// fprintf(stderr, "[M::%s::] ==> #########10#########, rlen: %d\n", __func__, rlen);
+
+	int32_t i, n = p->n, bd, len = MIN(rlen, dp_min_len), sublen, cnt, ei, li, ri;
+	int32_t sn =  len*er + 1;
+	kv_resize(uint64_t, *mt, (int64_t)p->n); 
+	mt->n = p->n; memset(mt->a, 0, sizeof(uint64_t)*p->n);
+	for (i = 0; i < n; i++) p->a[i].rid = ha_pt_cnt(pt, p->a[i].x);
+
+	for (i = cnt = 0, bd = -1, ei = -1; i < n; i++)
+	{
+		if((int32_t)(p->a[i].rid)<min_freq) continue;
+		sublen = p->a[i].pos + 1;
+		if(sublen > len) break;
+		else ei = i;
+		
+		if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd)
+		{
+			bd = p->a[i].pos;
+			cnt++;
+		}
+	}
+
+	// fprintf(stderr, "[M::%s::] ==> +cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);
+
+	
+	if(cnt >= sn) refine_select(p, 0, ei, sn, min_freq, mt, NULL, &li);
+	else
+	{
+		li = i-1; 
+		for (i = 0; i <= li; i++) mt->a[i] = 1;
+	} 
+	
+	
+	if(len < rlen)
+	{
+		for (i = n-1, cnt = 0, bd = rlen+1, ei = -1; i >= 0; i--)
+		{
+			if((int32_t)(p->a[i].rid)<min_freq) continue;
+			sublen = rlen - (p->a[i].pos + 1 - p->a[i].span);
+			if(sublen > len) break;
+			else ei = i;
+
+			if((int32_t)(p->a[i].pos) < bd)
+			{
+				bd = p->a[i].pos + 1 - p->a[i].span;
+				cnt++;
+			}
+		}
+
+		// fprintf(stderr, "[M::%s::] ==> -cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);
+
+		
+		if(cnt >= sn) refine_select(p, ei, n-1, sn, min_freq, mt, &ri, NULL);
+		else 
+		{
+			ri = i+1;
+			for (i = ri; i <= n-1; i++) mt->a[i] = 1;
+		} 
+		
+		// fprintf(stderr, "[M::%s::] ==> --cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);
+
+		if(ri - li >= 2)
+		{
+			li++; ri--;
+			sn =  (p->a[ri].pos - p->a[li].pos + p->a[li].span)*er + 1;
+			for (i = li, cnt = 0, bd = -1; i <= ri; i++)
+			{
+				if((int32_t)(p->a[i].rid)<min_freq) continue;
+				if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd)
+				{
+					bd = p->a[i].pos;
+					cnt++;
+					if(cnt >= sn) break;
+				}
+			}
+
+			if(cnt >= sn) refine_select(p, li, ri, sn, min_freq, mt, NULL, NULL);
+			else for (i = li; i <= ri; i++) mt->a[i] = 1;
+		}
+	}
+
+	// fprintf(stderr, "[M::%s::] ==> #########20#########, p->n: %u, n: %d\n", __func__, p->n, n);
+	for (i = sn = 0; i < n; i++)
+	{
+		if(mt->a[i])
+		{
+			p->a[sn] = p->a[i];
+			sn++;
+		}
+	}
+	// if(p->n != sn) fprintf(stderr, "[M::%s::] ==> #########21#########, p->n: %u, sn: %d\n", __func__, p->n, sn);
+	p->n = sn;
+	
+}
+
+inline int hf_dp(ha_mz1_v *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, 
+int32_t *rsi, int32_t *rei)
+{
+	return 0;
+}
+
+inline void hf_select(ha_mz1_v *p, int32_t si, int32_t ei, int32_t n, int32_t len, int32_t sample_dist, ha_mz1_t *b, int32_t force)
+{
+	if(ei - si <= 1) return;
+	int32_t ps = si < 0? 0 : p->a[si].pos;
+	int32_t pe = ei == n? len : p->a[ei].pos;
+	int32_t j, k, st = si + 1, en = ei;
+	int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499);
+	if (max_high_occ > MAX_MAX_HIGH_OCC)
+		max_high_occ = MAX_MAX_HIGH_OCC;
+	for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
+		b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[]
+	ks_heapmake_mz(k, b); // initialize the binomial heap
+	for (; j < en; ++j) { // if there are more, choose top max_high_occ
+		if (mz_lt(p->a[j], b[0])) { // then update the heap
+			b[0] = p->a[j], b[0].pos = j;
+			ks_heapdown_mz(0, k, b);
+		}
+	}
+	//ks_heapsort_mz(k, b); // sorting is not needed for now
+	for (j = 0; j < k; ++j)
+		if (b[j].rid < pe - ps || force)
+			p->a[b[j].pos].rid = 0;
+}
+
+static void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
 { // for high-occ minimizers, choose up to max_high_occ in each high-occ streak
-	int32_t i, last0 = -1, n = (int32_t)p->n, m = 0;
+	int32_t i, last0 = -1, n = (int32_t)p->n, m = 0, nw[2], min_len;
 	ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation
 
 	if (n == 0 || n == 1) return;
@@ -31,29 +242,47 @@ static void select_mz(ha_mz1_v *p, int len, int sample_dist)
 	for (i = 0; i <= n; ++i) {
 		if (i == n || p->a[i].rid == 0) {
 			if (i - last0 > 1) {
-				int32_t ps = last0 < 0? 0 : p->a[last0].pos;
-				int32_t pe = i == n? len : p->a[i].pos;
-				int32_t j, k, st = last0 + 1, en = i;
-				int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499);
-				if (max_high_occ > MAX_MAX_HIGH_OCC)
-					max_high_occ = MAX_MAX_HIGH_OCC;
-				for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
-					b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[]
-				ks_heapmake_mz(k, b); // initialize the binomial heap
-				for (; j < en; ++j) { // if there are more, choose top max_high_occ
-					if (mz_lt(p->a[j], b[0])) { // then update the heap
-						b[0] = p->a[j], b[0].pos = j;
-						ks_heapdown_mz(0, k, b);
-					}
-				}
-				//ks_heapsort_mz(k, b); // sorting is not needed for now
-				for (j = 0; j < k; ++j)
-					if (b[j].rid < pe - ps)
-						p->a[b[j].pos].rid = 0;
+				hf_select(p, last0, i, n, len, sample_dist, b, 0);
+				// int32_t ps = last0 < 0? 0 : p->a[last0].pos;
+				// int32_t pe = i == n? len : p->a[i].pos;
+				// int32_t j, k, st = last0 + 1, en = i;
+				// int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499);
+				// if (max_high_occ > MAX_MAX_HIGH_OCC)
+				// 	max_high_occ = MAX_MAX_HIGH_OCC;
+				// for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
+				// 	b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[]
+				// ks_heapmake_mz(k, b); // initialize the binomial heap
+				// for (; j < en; ++j) { // if there are more, choose top max_high_occ
+				// 	if (mz_lt(p->a[j], b[0])) { // then update the heap
+				// 		b[0] = p->a[j], b[0].pos = j;
+				// 		ks_heapdown_mz(0, k, b);
+				// 	}
+				// }
+				// //ks_heapsort_mz(k, b); // sorting is not needed for now
+				// for (j = 0; j < k; ++j)
+				// 	if (b[j].rid < pe - ps)
+				// 		p->a[b[j].pos].rid = 0;
 			}
 			last0 = i;
 		}
 	}
+
+	min_len = MAX(dp_min_len, (p->a[0].pos+1)+sample_dist);
+	for (i = 0, nw[0] = nw[1] = 0; i < n; i++)
+    {
+		nw[(p->a[i].rid!=0)]++;
+        if((p->a[i].pos + 1) > min_len) break;
+    }
+	if(nw[0]==0 && nw[1]>0) hf_select(p, -1, i, n, len, sample_dist, b, 1);
+
+	min_len = MAX(dp_min_len, (len - (p->a[n-1].pos + 1 - p->a[n-1].span))+sample_dist);
+	for (i = n-1, nw[0] = nw[1] = 0; i >= 0; i--)
+	{
+		nw[(p->a[i].rid!=0)]++;
+		if((len - (p->a[i].pos + 1 - p->a[i].span)) > min_len) break;
+	}
+	if(nw[0]==0 && nw[1]>0) hf_select(p, i, n, n, len, sample_dist, b, 1);
+
 	for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers
 		if (p->a[i].rid == 0)
 			p->a[n++] = p->a[i];
@@ -71,7 +300,7 @@ static void select_mz(ha_mz1_v *p, int len, int sample_dist)
  * @param is_hpc homopolymer-compressed or not
  * @param p      minimizers
  */
-void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
+void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt)
 {	///in default, w = 51, k = 51, is_hpc = 1
 	/**
 	 uint64_t x;
@@ -184,7 +413,127 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 	}
 	if (min.x != UINT64_MAX)
 		kv_push(ha_mz1_t, *p, min);
-	if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC);
+	if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC, dp_min_len);
+	if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);
+	for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts
+		p->a[i].rid = rid;
+}
+
+void ha_sketch_worse(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt)
+{	///in default, w = 51, k = 51, is_hpc = 1
+	/**
+	 uint64_t x;
+	 uint64_t rid:28, pos:27, rev:1, span:8;
+	 **/
+	extern void *ha_ct_table;
+	static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0 };
+	uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
+	int i, j, l, buf_pos, min_pos, kmer_span = 0;
+	ha_mz1_t buf[256], min = dummy;
+	tiny_queue_t tq;
+
+	assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
+	if (dbg_ct != NULL) dbg_ct->a.n = 0;
+	if (k_flag != NULL) {
+		kv_resize(uint8_t, k_flag->a, (uint64_t)len);
+		k_flag->a.n = len; 
+		memset(k_flag->a.a, 0, k_flag->a.n);
+	} 
+
+	memset(buf, 0xff, w * sizeof(ha_mz1_t));
+	memset(&tq, 0, sizeof(tiny_queue_t));
+	///len/w is the evaluated minimizer numbers
+	kv_resize(ha_mz1_t, *p, p->n + len/w);
+
+	for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
+		int c = seq_nt4_table[(uint8_t)str[i]];
+		ha_mz1_t info = dummy;
+		if (c < 4) { // not an ambiguous base
+			int z;
+			if (is_hpc) {
+				int skip_len = 1;
+				if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
+					for (skip_len = 2; i + skip_len < len; ++skip_len)
+						if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
+							break;
+					i += skip_len - 1; // put $i at the end of the current homopolymer run
+				}
+				tq_push(&tq, skip_len);
+				kmer_span += skip_len;
+				///how many bases that are covered by this HPC k-mer
+				///kmer_span includes at most k HPC elements
+				if (tq.count > k) kmer_span -= tq_shift(&tq);
+			} else kmer_span = l + 1 < k? l + 1 : k; 
+			///kmer_span should be used for HPC k-mer
+			///non-HPC k-mer, kmer_span should be k
+			///kmer_span is used to calculate anchor pos on reverse complementary strand
+
+			if (k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC
+
+			kmer[0] = (kmer[0] << 1 | (c&1))  & mask;                  // forward k-mer
+			kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
+			kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1))  << shift1; // reverse k-mer
+			kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
+			if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
+			z = kmer[1] < kmer[3]? 0 : 1; // strand
+			++l;
+			if (l >= k && kmer_span < 256) {
+				uint64_t y;
+				int32_t cnt, filtered;
+				y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
+				cnt = hf? ha_ft_cnt(hf, y) : 0;
+				filtered = (cnt >= 1<<28);
+				if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
+				if (!filtered) info.x = y, info.rid = pt? ha_pt_cnt(pt, y):cnt, info.pos = i, info.rev = z, info.span = kmer_span; // initially ha_mz1_t::rid keeps the k-mer count
+				if (k_flag != NULL) k_flag->a.a[i]++;
+				if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;
+			}
+		} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
+
+		//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
+		//i is the real base iterator, while l is the HPC base iterator
+		//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
+		//but even if l < k, infor is still stored into buf
+		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
+		if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
+			for (j = buf_pos + 1; j < w; ++j)
+				if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
+			for (j = 0; j < buf_pos; ++j)
+				if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
+		}
+		/**
+		 * There are three cases:
+		 * 1. info.x <= min.x, means info is a new minimizer
+		 * 2. info.x > min.x, info is not a new minimizer
+		 *    (1) buf_pos != min_pos, do nothing
+		 *    (2) buf_pos == min_pos, means current minimizer has moved outside the window
+		 * **/
+		///three cases: 1. 
+		if (info.x <= min.x) { // a new minimum; then write the old min
+			if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			min = info, min_pos = buf_pos;
+		} else if (buf_pos == min_pos) { // old min has moved outside the window
+			if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			///buf_pos == min_pos, means current minimizer has moved outside the window
+			///so for now we need to find a new minimizer at the current window (w k-mers)
+			for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
+				if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
+			for (j = 0; j <= buf_pos; ++j)
+				if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j;
+
+			if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
+				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
+					if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
+				for (j = 0; j <= buf_pos; ++j)
+					if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
+			}
+		}
+		if (++buf_pos == w) buf_pos = 0;
+	}
+	if (min.x != UINT64_MAX)
+		kv_push(ha_mz1_t, *p, min);
+	if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);
+	// if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC);
 	for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts
 		p->a[i].rid = rid;
 }