Merge branch 'master' into dev-lh3

CommandLines.cpp

@@ -122,7 +122,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->k_mer_length = 51;
     asm_opt->hic_mer_length = 31;
 	asm_opt->mz_win = 51;
+	asm_opt->mz_sample_dist = 500;
 	asm_opt->bf_shift = 37;
+	asm_opt->max_kmer_cnt = 2000;
 	asm_opt->high_factor = 5.0;
 	asm_opt->max_ov_diff_ec = 0.04;
 	asm_opt->max_ov_diff_final = 0.03;

CommandLines.h

@@ -43,7 +43,9 @@ typedef struct {
     int hic_mer_length;
     int bub_mer_length;
 	int mz_win;
+	int mz_sample_dist;
 	int bf_shift;
+	int max_kmer_cnt;
 	double high_factor; // coverage cutoff set to high_factor*hom_cov
 	double max_ov_diff_ec;
 	double max_ov_diff_final;

Overlaps.cpp

@@ -9449,8 +9449,7 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
 
 
 
-void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources, 
-ma_hit_t_alloc* reverse_sources, int id, const char* command)
+void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, const char* command)
 {
     long long i, j, Len;
     uint32_t tn;
@@ -9474,7 +9473,7 @@ ma_hit_t_alloc* reverse_sources, int id, const char* command)
             fprintf(stderr, "\n\n\nafter %s\n", command);
 
             fprintf(stderr, "****************ma_hit_t (%lld)ref_read: %.*s, len: %lu****************\n", 
-            i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), Get_READ_LENGTH(R_INF, i));
+					i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), (unsigned long)Get_READ_LENGTH(R_INF, i));
 
 
             fprintf(stderr, "sources Len: %d, is_fully_corrected: %d\n", 

anchor.cpp

@@ -57,8 +57,8 @@ int ha_ov_type(const overlap_region *r, uint32_t len)
 	else return r->x_pos_s == 0? 0 : 1;
 }
 
-void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag,
-kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
+void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain,
+						   kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
 {
 	uint32_t i, rlen;
 	uint64_t k, l;
@@ -73,7 +73,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
 	rlen = Get_READ_LENGTH(R_INF, rid); // read length
 
 	// get the list of anchors
-	ha_sketch_query(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, 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);
 	// minimizer of queried read
 	if (ab->mz.m > ab->old_mz_m) {
 		ab->old_mz_m = ab->mz.m;
@@ -205,7 +205,8 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
 
 
 void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
-int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
+								 int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar,
+								 kvec_t_u64_warp* dbg_ct)
 {
 	extern void *ha_flt_tab;
 	extern ha_pt_t *ha_idx;
@@ -302,4 +303,4 @@ int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp*
 void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
 {
 	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
-}
+}

htab.cpp

@@ -1,4 +1,3 @@
-#include <stdint.h>
 #include <zlib.h>
 #include <stdio.h>
 #include <string.h>
@@ -191,7 +190,7 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *
 		///so low 12 bits are not useful
 		uint64_t x = a[j] >> h->pre;
 		khint_t k;
-		if ((a[j]&mask) != (a[0]&mask)) continue;
+		assert((a[j]&mask) == (a[0]&mask));
 		if (create_new) {
 			///for 0-th counting, g->b = NULL
 			if (g->b)
@@ -367,9 +366,9 @@ int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a)
 		khint_t k;
 		int n;
 		ha_idxpos_t *p;
-		if ((a[j].x&mask) != (a[0].x&mask)) continue;
+		assert((a[j].x&mask) == (a[0].x&mask));
 		k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS);
-		if (k == kh_end(g->h)) continue;
+		if (k == kh_end(g->h)) continue; // TODO: understand why we sometimes come here
 		n = kh_key(g->h, k) & YAK_MAX_COUNT;
 		assert(n < YAK_MAX_COUNT);
 		p = &g->a[kh_val(g->h, k) + n];
@@ -520,7 +519,7 @@ typedef struct { // global data structure for kt_pipeline()
 	const yak_copt_t *opt;
 	const void *flt_tab;
 	int flag, create_new, is_store;
-	uint64_t n_seq; ///number of total reads
+	uint64_t n_mz, n_seq; ///number of total reads
 	kseq_t *ks;
 	UC_Read ucr;
 	ha_ct_t *ct;
@@ -558,7 +557,7 @@ static void worker_for_mz(void *data, long i, int tid)
 	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);
+	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->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));
@@ -575,8 +574,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		s->n_seq0 = p->n_seq;
 		if (p->rs_in && (p->flag & HAF_RS_READ)) {
 			while (p->n_seq < p->rs_in->total_reads) {
-				if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
-				{
+				if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) {
 					++p->n_seq;
 					continue;
 				}
@@ -689,6 +687,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 						ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
 			}
 			for (i = 0; i < s->n_seq; ++i) {
+				p->n_mz += s->mz[i].n;
 				free(s->mz[i].a);
 				if (!p->is_store) free(s->seq[i]);
 			}
@@ -809,6 +808,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 		gzclose(fp);
 	}
 	*n_seq = pl.n_seq;
+	if (pl.opt->w > 1) fprintf(stderr, "[M::%s] collected %ld minimizers\n", __func__, (long)pl.n_mz);
 	return pl.ct;
 }
 
@@ -849,34 +849,41 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
  * High count filter table *
  ***************************/
 
-KHASHL_SET_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, kh_hash_dummy, kh_eq_generic)
+// Warning: the max count is 32767
+KHASHL_MAP_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, int16_t, kh_hash_dummy, kh_eq_generic)
 
-static yak_ft_t *gen_hh(const ha_ct_t *h)
+static yak_ft_t *gen_hh(const ha_ct_t *h, int max_cnt)
 {
 	int i;
 	yak_ft_t *hh;
+	if (max_cnt > YAK_MAX_COUNT - 1) max_cnt = YAK_MAX_COUNT - 1;
+	if (max_cnt > INT16_MAX - 1) max_cnt = INT16_MAX - 1;
 	hh = yak_ft_init();
 	yak_ft_resize(hh, h->tot * 2);
 	for (i = 0; i < 1<<h->pre; ++i) {
 		yak_ct_t *ht = h->h[i].h;
-		khint_t k;
+		khint_t k, l;
 		for (k = 0; k < kh_end(ht); ++k) {
 			if (kh_exist(ht, k)) {
 				uint64_t y = kh_key(ht, k) >> h->pre << YAK_COUNTER_BITS | i;
 				int absent;
-				yak_ft_put(hh, y, &absent);
+				l = yak_ft_put(hh, y, &absent);
+				if (absent) {
+					int cnt = kh_key(ht, k) & YAK_MAX_COUNT;
+					kh_val(hh, l) = cnt > max_cnt? INT16_MAX : cnt;
+				}
 			}
 		}
 	}
 	return hh;
 }
 
-int ha_ft_isflt(const void *hh, uint64_t y)
+int32_t ha_ft_cnt(const void *hh, uint64_t y)
 {
 	yak_ft_t *h = (yak_ft_t*)hh;
 	khint_t k;
 	k = yak_ft_get(h, y);
-	return k == kh_end(h)? 0 : 1;
+	return k == kh_end(h)? 0 : kh_val(h, k) == INT16_MAX? INT32_MAX : kh_val(h, k);
 }
 
 void ha_ft_destroy(void *h)
@@ -942,7 +949,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
 		if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
 	}
 	ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
-	flt_tab = gen_hh(h);
+	flt_tab = gen_hh(h, 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);
@@ -1006,7 +1013,6 @@ int query_ct_index(void* ct_idx, uint64_t hash)
 	return kh_key(g->h, k)&YAK_MAX_COUNT;
 }
 
-
 int write_ct_index(void *i_ct_idx, char* file_name)
 {
 	char* gfa_name = (char*)malloc(strlen(file_name)+25);
@@ -1139,7 +1145,7 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_op
 
 	ha_pt_t *ha_idx = NULL;
 	char mode = 0;
-	int f_flag, absent, i;
+	int f_flag = 0, absent, i;
 	double index_time, index_s_time, pos_time, pos_s_time;
 
 	

htab.h

@@ -31,9 +31,8 @@ extern void *ha_flt_tab_hp;
 extern ha_pt_t *ha_idx_hp;
 extern void *ha_ct_table;
 
-
 void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
-int ha_ft_isflt(const void *hh, uint64_t y);
+int32_t ha_ft_cnt(const void *hh, uint64_t y);
 void ha_ft_destroy(void *h);
 
 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);
@@ -59,8 +58,7 @@ 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);
-void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, 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);
 int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het);
 void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);
 

ksort.h

@@ -40,7 +40,32 @@ typedef struct {
 
 #define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
 
-#define KSORT_INIT(name, type_t, __sort_lt)								\
+#define KSORT_INIT(name, type_t, __sort_lt) \
+	void ks_heapdown_##name(size_t i, size_t n, type_t l[]) \
+	{ \
+		size_t k = i; \
+		type_t tmp = l[i]; \
+		while ((k = (k << 1) + 1) < n) { \
+			if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k; \
+			if (__sort_lt(l[k], tmp)) break; \
+			l[i] = l[k]; i = k; \
+		} \
+		l[i] = tmp; \
+	} \
+	void ks_heapmake_##name(size_t lsize, type_t l[]) \
+	{ \
+		size_t i; \
+		for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \
+			ks_heapdown_##name(i, lsize, l); \
+	} \
+	void ks_heapsort_##name(size_t lsize, type_t l[]) \
+	{ \
+		size_t i; \
+		for (i = lsize - 1; i > 0; --i) { \
+			type_t tmp; \
+			tmp = *l; *l = l[i]; l[i] = tmp; ks_heapdown_##name(0, i, l); \
+		} \
+	} \
 	static inline void __ks_insertsort_##name(type_t *s, type_t *t)		\
 	{																	\
 		type_t *i, *j, swap_tmp;										\

sketch.cpp

@@ -4,6 +4,10 @@
 #include <string.h>
 #include "kvec.h"
 #include "htab.h"
+#include "ksort.h"
+
+#define MAX_HIGH_OCC     8   // TODO: don't hard code if we need to tune this parameter
+#define MAX_MAX_HIGH_OCC 16
 
 typedef struct { // a simplified version of kdq
 	int front, count;
@@ -25,6 +29,56 @@ static inline int tq_shift(tiny_queue_t *q)
 	return x;
 }
 
+static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b)
+{
+	return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x));
+}
+
+#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)
+{ // 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;
+	ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation
+
+	if (n == 0 || n == 1) return;
+	assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers
+	for (i = 0; i < n; ++i)
+		if (p->a[i].rid != 0) ++m;
+	if (m == 0) return; // no high-frequency k-mers; do nothing
+	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;
+			}
+			last0 = i;
+		}
+	}
+	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];
+	p->n = n;
+}
+
 /**
  * Find symmetric (w,k)-minimizers on a DNA sequence
  *
@@ -36,134 +90,28 @@ static inline int tq_shift(tiny_queue_t *q)
  * @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)
-{	///in default, w = 51, k = 51, is_hpc = 1
-	/**
-	 uint64_t x;
-	 uint64_t rid:28, pos:27, rev:1, span:8;
-	 **/
-	static const ha_mz1_t dummy = { UINT64_MAX, 0, 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));
-	///sizeof(ha_mz1_t) = 16
-	memset(buf, 0xff, w * 16);
-	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;
-				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
-			///so for non-HPC k-mer, kmer_span should be k in any case?
-			///kmer_span is used to calculate anchor pos on reverse complementary strand
-
-			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;
-				y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
-				if (hf == 0 || ha_ft_isflt(hf, y) == 0)
-					info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
-			}
-		} 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 (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
-			for (j = 0; j < buf_pos; ++j)
-				if (min.x == buf[j].x && 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 (min.x >= buf[j].x) 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 (min.x >= buf[j].x) 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 (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
-				for (j = 0; j <= buf_pos; ++j)
-					if (min.x == buf[j].x && 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);
-}
-
-
-
-void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, 
-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)
 {	///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;
-	if(dbg_ct != NULL) dbg_ct->a.n = 0;
-
-	static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
-	uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0}, filtered;
+	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;
-	if(k_flag != NULL)
-	{
+
+	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);
 	} 
 
-
-	assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
-	///sizeof(ha_mz1_t) = 16
-	memset(buf, 0xff, w * 16);
+	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);
@@ -191,7 +139,7 @@ kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
 			///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
+			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;
@@ -202,20 +150,17 @@ kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
 			++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]);
-
-				filtered = 0;
-				if(hf != 0) filtered = ha_ft_isflt(hf, y);
-				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 (hf == 0 || ha_ft_isflt(hf, y) == 0)
-				if(filtered == 0)
-					info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
-				if(k_flag != NULL) k_flag->a.a[i]++;
-				if(k_flag != NULL && filtered == 1) k_flag->a.a[i]++;
+				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 = 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)
@@ -223,9 +168,9 @@ kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
 		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 (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[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 (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
+				if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
 		}
 		/**
 		 * There are three cases:
@@ -243,19 +188,22 @@ kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
 			///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 (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
+				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 (min.x >= buf[j].x) min = buf[j], min_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 (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
+					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 (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[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 (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;
+}