r330: unimap heuristic for mz sampling

CommandLines.cpp

@@ -101,6 +101,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->thread_num = 1;
     asm_opt->k_mer_length = 51;
 	asm_opt->mz_win = 51;
+	asm_opt->sample_dist = 200;
 	asm_opt->bf_shift = 37;
 	asm_opt->max_kmer_cnt = 2000;
 	asm_opt->high_factor = 5.0;

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.13-r329-dirty"
+#define HA_VERSION "0.13-r330-dirty"
 
 #define VERBOSE 0
 
@@ -34,6 +34,7 @@ typedef struct {
     int thread_num;
     int k_mer_length;
 	int mz_win;
+	int sample_dist;
 	int bf_shift;
 	int max_kmer_cnt;
 	double high_factor; // coverage cutoff set to high_factor*hom_cov

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, 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.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;

htab.cpp

@@ -556,7 +556,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, 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, 0, 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));

htab.h

@@ -58,7 +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, 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);
 
 static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function

sketch.cpp

@@ -6,7 +6,7 @@
 #include "htab.h"
 #include "ksort.h"
 
-#define MAX_HIGH_OCC     3   // TODO: don't hard code if we need to tune this parameter
+#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
@@ -37,15 +37,13 @@ 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 max_high_occ)
+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
-	if (max_high_occ > MAX_MAX_HIGH_OCC)
-		max_high_occ = MAX_MAX_HIGH_OCC;
 	for (i = 0; i < n; ++i)
 		if (p->a[i].rid != 0) ++m;
 	if (m == 0) return; // no high-frequency k-mers; do nothing
@@ -55,6 +53,9 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ)
 				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
@@ -66,8 +67,7 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ)
 				}
 				//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;
+					p->a[b[j].pos].rid = 0;
 			}
 			last0 = i;
 		}
@@ -89,7 +89,7 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ)
  * @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, 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;
@@ -202,7 +202,7 @@ 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);
-	select_mz(p, len, MAX_HIGH_OCC);
+	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;
 }