diff --git a/Assembly.cpp b/Assembly.cpp
index 358bbe2..ba5aa97 100644
--- a/Assembly.cpp
+++ b/Assembly.cpp
@@ -13,7 +13,7 @@
 #include "rcut.h"
 
 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, st_mt_t *sp);
 void ha_get_ug_candidates(ha_abuf_t *ab, int64_t rid, ma_utg_t *u, ma_utg_v *ua, 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, double chain_match_rate);
 void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
@@ -438,6 +438,7 @@ typedef struct {
     kvec_t_u8_warp k_flag;
     overlap_region tmp_region;
     ma_utg_v *ua;
+    st_mt_t sp;
 } ha_ovec_buf_t;
 
 ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
@@ -454,6 +455,7 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
     kv_init(b->b_buf.a);
     kv_init(b->r_buf.a);
     kv_init(b->k_flag.a);
+    kv_init(b->sp);
 	b->ab = ha_abuf_init();
 	if (!b->is_final) {
 		init_Cigar_record(&b->cigar1);
@@ -478,6 +480,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
     kv_destroy(b->b_buf.a);
     kv_destroy(b->r_buf.a);
     kv_destroy(b->k_flag.a);
+    kv_destroy(b->sp);
 	if (!b->is_final) {
 		destory_Cigar_record(&b->cigar1);
 		destory_Graph(&b->POA_Graph);
@@ -541,7 +544,7 @@ static void worker_ovec(void *data, long i, int tid)
 	int fully_cov, abnormal;
 
     ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 
-    0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
+    0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
 
 	clear_Cigar_record(&b->cigar1);
 	clear_Round2_alignment(&b->round2);
@@ -599,7 +602,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
         int fully_cov, abnormal, q_idx = k;
 
         ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 
-        0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), &(R_INF_FLAG.candidate_count[q_idx]));
+        0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), &(R_INF_FLAG.candidate_count[q_idx]), &(b->sp));
 
         clear_Cigar_record(&b->cigar1);
         clear_Round2_alignment(&b->round2);
@@ -625,7 +628,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
 
         pthread_mutex_lock(&R_INF_FLAG.OutputMutex);
 
-        fprintf(R_INF_FLAG.fp, "\n>%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
+        fprintf(R_INF_FLAG.fp, "\n>%.*s, rid: %ld\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i), i);
         fprintf(R_INF_FLAG.fp, "%d-th round, len: %lu, hom_cov: %d, max_n_chain: %d\n", 
             asm_opt.number_of_round, Get_READ_LENGTH(R_INF, i), asm_opt.hom_cov, asm_opt.max_n_chain);
 
@@ -1291,7 +1294,7 @@ static void worker_ov_final(void *data, long i, int tid)
 
 	//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
     ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 0.001, 
-    asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
+    asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
 
 	overlap_region_sort_y_id(b->olist.list, b->olist.length);
 	ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
@@ -1357,7 +1360,7 @@ static void worker_ov_final_high_het(void *data, long i, int tid)
     ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
 
     ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, HIGH_HET_ERROR_RATE, 
-    asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
+    asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
 
     overlap_region_sort_y_id(b->olist.list, b->olist.length);
     ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
diff --git a/CommandLines.cpp b/CommandLines.cpp
index 9ddae12..7f92410 100644
--- a/CommandLines.cpp
+++ b/CommandLines.cpp
@@ -154,6 +154,7 @@ 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_rewin = 1000;
 	asm_opt->mz_sample_dist = 500;
 	asm_opt->bf_shift = 37;
 	asm_opt->max_kmer_cnt = 2000;
diff --git a/CommandLines.h b/CommandLines.h
index f682237..db780f7 100644
--- a/CommandLines.h
+++ b/CommandLines.h
@@ -4,7 +4,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.15.5-r358"
+#define HA_VERSION "0.15.5-r366"
 
 #define VERBOSE 0
 
@@ -46,6 +46,7 @@ typedef struct {
     int hic_mer_length;
     int bub_mer_length;
 	int mz_win;
+    int mz_rewin;
 	int mz_sample_dist;
 	int bf_shift;
 	int max_kmer_cnt;
diff --git a/Hash_Table.cpp b/Hash_Table.cpp
index 421a596..c1b7dbd 100644
--- a/Hash_Table.cpp
+++ b/Hash_Table.cpp
@@ -9,6 +9,8 @@ pthread_mutex_t output_mutex;
 #define overlap_region_key(a) ((a).y_id)
 KRADIX_SORT_INIT(overlap_region_sort, overlap_region, overlap_region_key, member_size(overlap_region, y_id))
 
+#define normal_w(x, y) ((x)>=(y)?(x)/(y):1)
+
 void overlap_region_sort_y_id(overlap_region *a, long long n)
 {
 	radix_sort_overlap_region_sort(a, a + n);
@@ -467,7 +469,46 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
     }
 }
 
+void print_chain(k_mer_hit* a, long long a_n, Chain_Data* dp, long long topN)
+{
+    fprintf(stderr, "topN: %lld\n", topN);
+    long long max_score = -1, max_i = -1, max_n = 0;;
+    long long ss, i, j, current_j;
+    kvec_t(long long) si; kv_init(si);
+    for (ss = 0; ss < topN && ss < a_n; ss++){
+        for (i = 0, max_i = -1, max_score = -1; i < a_n; ++i) {
+            for (j = 0; j < (long long)si.n; j++){
+                if(i == si.a[j]) break;
+            }
+            if(j < (long long)si.n) continue;
+            if(dp->score[i] > max_score) max_score = dp->score[i], max_i = i;
+        }
+        if(max_i < 0) continue;
+        j = max_i; max_n = 0;
+        while (j >= 0)
+        {
+            current_j = j;
+            if(current_j == -1) continue;
+            j = dp->pre[j];
+            max_n++;
+        }
 
+        fprintf(stderr, "\nmax_i: %lld, max_score: %lld, max_n: %lld\n", max_i, max_score, max_n);
+        
+        j = max_i; 
+        while (j >= 0)
+        {
+            current_j = j;
+            if(current_j == -1) continue;
+
+            kv_push(long long, si, current_j);
+            j = dp->pre[j];
+            fprintf(stderr, "self_offset: %u, offset: %u, cnt: %u, score: %d\n", 
+            a[current_j].self_offset, a[current_j].offset, a[current_j].cnt, dp->score[current_j]);
+        }
+    }
+    kv_destroy(si);
+}
 
 
 long long get_chainLen(long long x_beg, long long x_end, long long xLen, 
@@ -556,7 +597,8 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
 			break;
 	if (i < n_a) return -1;
 	bw_pen = 1.0 / bw_thres;
-	dp->score[0] = a[0].good? min_sc : min_sc>>1;
+	// dp->score[0] = a[0].good? min_sc : min_sc>>1;
+    dp->score[0] = normal_w(min_sc, (int64_t)a[0].cnt);
 	dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0, dp->occ[0] = 1;
 	for (i = 1; i < n_a; ++i) {
 		int32_t score, dg;
@@ -570,7 +612,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
 		dg = dx < dy? dx : dy;
 		if (dd > THRESHOLD_MAX_SIZE && dd > dg * bw_thres) break;
 		score = dg < min_sc? dg : min_sc;
-		if (!a[i].good) score >>= 1;
+        score = normal_w(score, (int64_t)a[i].cnt);
 		gap_rate = (double)tot_indel / tot_len;
 		score -= (int)(gap_rate * score * bw_pen);
 		dp->score[i] = dp->score[i-1] + score;
@@ -601,7 +643,6 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
     resize_Chain_Data(dp, a_n);
 
 	ret = ha_chain_check(a, a_n, dp, min_score, band_width_threshold);
-
 	if (ret > 0) {
 		a_n = ret;
 		goto skip_dp;
@@ -617,7 +658,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
         pos = a[i].offset;
         self_pos = a[i].self_offset;
         max_j = -1;
-        max_score = a[i].good? min_score : min_score>>1;
+        max_score = normal_w(min_score, (int64_t)a[i].cnt);
         max_indels = 0;
         max_self_length = 0;
 
@@ -648,7 +689,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
             score = distance_min < min_score? distance_min : min_score;
             ///need to be fixed in r305
 			///if (!a[j].good) score = (score >> 1) + (score & 1);
-            if (!a[j].good) score >>= 1;
+            score = normal_w(score, (int64_t)a[j].cnt);
 
             gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
             ///if the gap rate > 0.06, score will be negative
@@ -684,6 +725,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
     }
 
     ///debug_chain(a, a_n, dp);
+    // if((*result).x_id == 2162668 && (*result).y_id == 182804) print_chain(a, a_n, dp, 10);
     
 skip_dp:
 
diff --git a/Hash_Table.h b/Hash_Table.h
index 39f6d6e..2948309 100644
--- a/Hash_Table.h
+++ b/Hash_Table.h
@@ -110,8 +110,8 @@ typedef struct
 
 typedef struct
 {
-	uint32_t readID:30, strand:1, good:1;
-	uint32_t offset, self_offset;
+	uint32_t readID:31, strand:1;
+	uint32_t offset, self_offset, cnt;
 } k_mer_hit;
 
 typedef struct {
diff --git a/Overlaps.cpp b/Overlaps.cpp
index 82a05f9..a0dbade 100644
--- a/Overlaps.cpp
+++ b/Overlaps.cpp
@@ -153,7 +153,7 @@ void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend)
 ///if ug == NULL, nsg should be equal to read_sg
 inline uint32_t check_different_haps(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
 uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
-R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
+R_to_U* ruIndex, uint8_t* is_r_het, uint32_t min_edge_length, uint32_t stops_threshold)
 {
 	uint32_t vEnd, qn, tn, j, is_Unitig, uId;
 	long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;
@@ -187,7 +187,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
         b_max.b_0 = b_0;
     }
 
-	uint32_t max_count = 0, min_count = 0;
+	uint32_t max_count = 0, min_count = 0, n_het = 0, n_hom = 0;
 	ma_utg_t *node_min = NULL, *node_max = NULL;
 	if(ug != NULL)
 	{
@@ -217,7 +217,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 
 				/************************BUG: don't forget****************************/
 				if(reverse_sources[qn].length > 0) min_count++;
-				///if(reverse_sources[qn].length >= 0) min_count++;
+                if((is_r_het[qn] & C_HET) || (is_r_het[qn] & P_HET)) n_het++;
+                n_hom++;
 				/************************BUG: don't forget****************************/
 				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
 				{
@@ -268,7 +269,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 
 			/************************BUG: don't forget****************************/
 			if(reverse_sources[qn].length > 0) min_count++;
-			///if(reverse_sources[qn].length >= 0) min_count++;
+			if((is_r_het[qn] & C_HET) || (is_r_het[qn] & P_HET)) n_het++;
+            n_hom++;
 			/************************BUG: don't forget****************************/
 
 			for (j = 0; j < (long long)reverse_sources[qn].length; j++)
@@ -300,7 +302,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 	}
 
 	if(min_count == 0) return UNAVAILABLE;
-	if(max_count > min_count*asm_opt.purge_simi_thres/**DIFF_HAP_RATE**/) return PLOID;
+	if(max_count > min_count*asm_opt.purge_simi_thres && n_het >= n_hom*HET_HOM_RATE) return PLOID;
 	return NON_PLOID;
 }
 
@@ -476,107 +478,6 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t *min_count, uint32_t
     (*max_count) = inp_match;
     (*min_count) = hap_match;
 }
-inline uint32_t check_different_haps_base(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
-uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
-R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
-{
-	uint32_t vEnd, qn;
-	long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;
-
-	b_0->b.n = b_1->b.n = 0;
-	if(get_unitig(nsg, ug, v_0, &vEnd, &tmp, &ELen_0, &max_stop_nodeLen, &max_stop_baseLen, 
-	stops_threshold, b_0) == LOOP)
-	{
-		return UNAVAILABLE;
-	}
-	if(get_unitig(nsg, ug, v_1, &vEnd, &tmp, &ELen_1, &max_stop_nodeLen, &max_stop_baseLen, 
-	stops_threshold, b_1) == LOOP)
-	{
-		return UNAVAILABLE;
-	}
-	if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;
-
-	rIdContig b_max, b_min;
-	b_max.b_0 = b_min.b_0 = NULL;
-	b_max.offset = b_max.readI = b_max.untigI = 0;
-	b_min.offset = b_min.readI = b_min.untigI = 0;
-
-	if(ELen_0<=ELen_1)
-	{
-        b_min.b_0 = b_0;
-        b_max.b_0 = b_1;
-    }
-    else
-    {
-        b_min.b_0 = b_1;
-        b_max.b_0 = b_0;
-    }
-
-	uint32_t max_count = 0, min_count = 0;
-	ma_utg_t *node_max = NULL;
-	if(ug != NULL)
-	{
-		/*****************************label all unitigs****************************************/
-		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
-		{
-			node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
-			///each read
-			for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
-			{
-				qn = (node_max->a[b_max.readI]>>33);
-				set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1, &(read_sg->seq[qn].c));
-			}
-		}
-		/*****************************label all unitigs****************************************/
-
-        calculate_match_cover(b_min.b_0->b.a, b_min.b_0->b.n, nsg, ug, read_sg,
-        reverse_sources, ruIndex, &min_count, &max_count);
-
-		/*****************************label all unitigs****************************************/
-		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
-		{
-			node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
-			///each read
-			for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
-			{
-				qn = (node_max->a[b_max.readI]>>33);
-				ruIndex->index[qn] = (uint32_t)-1;
-			}
-		}
-		/*****************************label all unitigs****************************************/
-	}
-	else
-	{
-		/*****************************label all reads****************************************/
-		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
-		{
-			qn = (b_max.b_0->b.a[b_max.untigI]>>1);
-			set_R_to_U(ruIndex, qn, 1, 1, &(read_sg->seq[qn].c));
-		}
-		/*****************************label all reads****************************************/
-
-		calculate_match_cover(b_min.b_0->b.a, b_min.b_0->b.n, nsg, NULL, read_sg,
-        reverse_sources, ruIndex, &min_count, &max_count);
-
-		/*****************************label all reads****************************************/
-		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
-		{
-			qn = (b_max.b_0->b.a[b_max.untigI]>>1);
-			ruIndex->index[qn] = (uint32_t)-1;
-		}
-		/*****************************label all reads****************************************/
-	}
-
-    // if(v_0 == 67 && v_1 == 510)
-    // {
-    //     fprintf(stderr, "v_0-%u, v_1-%u, min_count-%u, max_count-%u\n", v_0, v_1, min_count, max_count);
-    // }
-
-	if(min_count == 0) return UNAVAILABLE;
-	if(max_count > min_count*asm_opt.purge_simi_thres/**DIFF_HAP_RATE**/) return PLOID;
-	return NON_PLOID;
-}
-
 
 asg_t *asg_init(void)
 {
@@ -12504,7 +12405,7 @@ ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov)
                 rid = (ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33));
                 uCov += cov->cov[rid];
 
-                if(t_ch) t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
+                if(t_ch) t_ch->ir_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
             }
         }
 
@@ -12520,7 +12421,7 @@ ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov)
                 rid = (ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33));
                 uLen += read_sg->seq[rid].len;
 
-                if(t_ch) t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
+                if(t_ch) t_ch->ir_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
             }
             if(occ >= thre_pri) break;
         }
@@ -12693,7 +12594,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp)
 
 void set_ug_coverage_aggressive(ma_ug_t *ug, uint32_t uID, asg_t* read_g, 
 const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag,
-trans_chain* t_ch, long long het_cov_thres)
+uint8_t* is_r_het, long long het_cov_thres)
 {
     ma_utg_t *u = &(ug->u.a[uID]);
     uint32_t k, j, rId, tn, is_Unitig;
@@ -12791,7 +12692,7 @@ trans_chain* t_ch, long long het_cov_thres)
 
         if((R_bases <= 0) || ((C_bases/R_bases) <= het_cov_thres))
         {
-            t_ch->is_r_het[rId] |= C_HET;
+            is_r_het[rId] |= C_HET;
         }
     }
 
@@ -12820,7 +12721,7 @@ trans_chain* t_ch, long long het_cov_thres)
 }
 
 
-void set_r_het_flag(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, trans_chain* t_ch)
+void set_r_het_flag(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* is_r_het)
 {
     uint64_t m, dip_thre_max, dip_thres;
     uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t));
@@ -12842,7 +12743,7 @@ void set_r_het_flag(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_all
     {
         dip_thres = dip_thre_max;
         ///if(ug->u.a[m].n <= dip_thre_max) dip_thres = dip_thre_max * 1.1;
-        set_ug_coverage_aggressive(ug, m, sg, coverage_cut, sources, ruIndex, primary_flag, t_ch, dip_thres);
+        set_ug_coverage_aggressive(ug, m, sg, coverage_cut, sources, ruIndex, primary_flag, is_r_het, dip_thres);
     }
     free(primary_flag); 
 }
@@ -12859,7 +12760,6 @@ trans_chain* init_trans_chain(ma_ug_t *ug, uint64_t r_num)
     memset(x->rUidx, -1, x->r_num*sizeof(uint32_t));
     MALLOC(x->rUpos, r_num);
     memset(x->rUpos, -1, x->r_num*sizeof(uint64_t));
-    CALLOC(x->is_r_het, x->r_num);
     memset(&(x->b_buf_0), 0, sizeof(buf_t));
     memset(&(x->b_buf_1), 0, sizeof(buf_t));
     kv_init(x->topo_buf);
@@ -12939,7 +12839,6 @@ void destory_trans_chain(trans_chain **x)
         kv_destroy((*x)->k_t_b);
         free((*x)->rUidx);
         free((*x)->rUpos);
-        free((*x)->is_r_het);
         uint32_t k;
         for (k = 0; k < (*x)->bed.n; k++) kv_destroy((*x)->bed.a[k]);
         kv_destroy((*x)->bed);
@@ -13046,7 +12945,7 @@ void write_trans_chain(trans_chain* t_ch, const char *fn)
     FILE* fp = fopen(buf, "w");
 
     fwrite(&t_ch->r_num, sizeof(t_ch->r_num), 1, fp);
-    fwrite(t_ch->is_r_het, sizeof(uint8_t), t_ch->r_num, fp);
+    fwrite(t_ch->ir_het, sizeof(uint8_t), t_ch->r_num, fp);
     
     uint32_t i;
     fwrite(&t_ch->bed.n, sizeof(t_ch->bed.n), 1, fp);
@@ -13082,8 +12981,8 @@ trans_chain* load_hc_trans(const char *fn)
     CALLOC(t_ch, 1);
 
     flag += fread(&t_ch->r_num, sizeof(t_ch->r_num), 1, fp);
-    MALLOC(t_ch->is_r_het, t_ch->r_num);
-    flag += fread(t_ch->is_r_het, sizeof(uint8_t), t_ch->r_num, fp);
+    MALLOC(t_ch->ir_het, t_ch->r_num);
+    flag += fread(t_ch->ir_het, sizeof(uint8_t), t_ch->r_num, fp);
 
     uint32_t i;
     flag += fread(&t_ch->bed.n, sizeof(t_ch->bed.n), 1, fp);
@@ -13351,7 +13250,12 @@ int max_hang, int min_ovlp, R_to_U* ruIndex, bub_label_t* b_mask_t)
     kv_init(new_rtg_edges.a); kv_init(d_edges.a);
 
     ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
-    adjust_utg_advance(sg, ug, reverse_sources, ruIndex, b_mask_t);
+
+    uint8_t* is_r_het = NULL;
+    CALLOC(is_r_het, sg->n_seq);
+    set_r_het_flag(ug, sg, coverage_cut, sources, ruIndex, is_r_het);
+
+    adjust_utg_advance(sg, ug, reverse_sources, ruIndex, b_mask_t, is_r_het);
     asg_t* nsg = (*ug).g;
     uint32_t v, n_vtx = nsg->n_seq;
     for (v = 0; v < n_vtx; ++v)
@@ -13366,6 +13270,7 @@ int max_hang, int min_ovlp, R_to_U* ruIndex, bub_label_t* b_mask_t)
 
     kv_destroy(new_rtg_edges.a);  kv_destroy(d_edges.a);
     horder_clean_sg_by_utg(sg, ug);
+    free(is_r_het);
     return ug;
 }
 
@@ -13377,8 +13282,8 @@ int max_hang, int min_ovlp, R_to_U* ruIndex, bub_label_t* b_mask_t)
     p->r_num = sg->n_seq; p->u_num = ug->u.n;
     kv_malloc(p->bed, p->u_num); p->bed.n = p->u_num;
     for (k = 0; k < p->bed.n; k++) kv_init(p->bed.a[k]);
-    CALLOC(p->is_r_het, p->r_num);
-    kv_u_trans_t *ta = get_utg_ovlp(&ug, sg, sources, reverse_sources, coverage_cut, ruIndex, max_hang, min_ovlp, NULL, b_mask_t, p->is_r_het);
+    CALLOC(p->ir_het, p->r_num);
+    kv_u_trans_t *ta = get_utg_ovlp(&ug, sg, sources, reverse_sources, coverage_cut, ruIndex, max_hang, min_ovlp, NULL, b_mask_t, p->ir_het);
     p->k_trans = *ta; free(ta);
     return p;
 }
@@ -13539,7 +13444,7 @@ void set_trio_flag_by_cov(ma_ug_t *ug, asg_t *read_g, hap_cov_t *cov)
         for (k = 0; k < u->n; k++)
         {
             if(R_INF.trio_flag[u->a[k]>>33]&SET_TRIO) continue;
-            if(cov->t_ch->is_r_het[u->a[k]>>33] == N_HET) continue;
+            if(cov->t_ch->ir_het[u->a[k]>>33] == N_HET) continue;
             R_INF.trio_flag[u->a[k]>>33] |= flag;
         }
     }
@@ -13639,23 +13544,6 @@ const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t*
 	return R_bases == 0? 0 : C_bases/R_bases;
 }
 
-
-void print_r_het(hap_cov_t *cov, uint8_t* trio_flag, const char* cmd)
-{
-    if(cov && cov->t_ch)
-    {
-        fprintf(stderr, "\n+%s-is_r_het[1369536]=%u\n", cmd, cov->t_ch->is_r_het[1369536]);
-        fprintf(stderr, "+%s-is_r_het[5097804]=%u\n", cmd, cov->t_ch->is_r_het[5097804]);
-        fprintf(stderr, "+%s-is_r_het[603738]=%u\n", cmd, cov->t_ch->is_r_het[603738]);
-    }   
-    if(trio_flag)
-    {
-        fprintf(stderr, "+%s-trio_flag[1369536]=%u\n", cmd, trio_flag[1369536]);
-        fprintf(stderr, "+%s-trio_flag[5097804]=%u\n", cmd, trio_flag[5097804]);
-        fprintf(stderr, "+%s-trio_flag[603738]=%u\n", cmd, trio_flag[603738]);
-    }
-}
-
 void kt_u_trans_t_idx(kv_u_trans_t *ta, uint32_t n)
 {
     radix_sort_u_trans(ta->a, ta->a + ta->n);
@@ -14774,7 +14662,7 @@ uint32_t cal_trio_vec(buf_t* b, ma_ug_t *ug, float thres)
 }
 
 int cut_trio_tip_primary(asg_t *g, ma_ug_t *ug, uint32_t max_ext, uint32_t trio_flag, uint32_t keep_out_node,
-asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t min_edge_length)
+asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het, uint32_t min_edge_length)
 {
     double startTime = Get_T();
     uint32_t n_vtx = g->n_seq * 2, v, w, i, cnt = 0, tipEvaluateLen, flag, inner_flag, operation, tip_trio_flag;
@@ -14856,8 +14744,8 @@ asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t min_e
                 if(operation == CUT) break;
                 if(aw[i].del) continue;
                 if(aw[i].v == (b.b.a[b.b.n-1]^1)) continue;
-                inner_flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, b.b.a[b.b.n-1]^1, aw[i].v,
-                reverse_sources, &b_0, &b_1, ruIndex, min_edge_length, 1);
+                inner_flag = check_different_haps(g, ug, read_sg, b.b.a[b.b.n-1]^1, aw[i].v,
+                reverse_sources, &b_0, &b_1, ruIndex, is_r_het, min_edge_length, 1);
                 if(inner_flag == NON_PLOID) operation = CUT;
             }
         }
@@ -14989,8 +14877,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg
                 n_reduced++;
 
                 operation = TRIM;
-                flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, av[v_maxLen_i].v, av[i].v,
-                reverse_sources, &b_0, &b_1, ruIndex, min_edge_length, 1);
+                flag = check_different_haps(g, ug, read_sg, av[v_maxLen_i].v, av[i].v, reverse_sources, 
+                &b_0, &b_1, ruIndex, cov->is_r_het, min_edge_length, 1);
                 // #define UNAVAILABLE (uint32_t)-1
                 // #define PLOID 0
                 // #define NON_PLOID 1
@@ -15096,8 +14984,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre
             {
                 n_reduced++;
                 operation = TRIM;
-                flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
-                reverse_sources, &b_0, &b_1, ruIndex, min_edge_length, stops_threshold);
+                flag = check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
+                reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, min_edge_length, stops_threshold);
                 // #define UNAVAILABLE (uint32_t)-1
                 // #define PLOID 0
                 // #define NON_PLOID 1
@@ -15290,8 +15178,8 @@ hap_cov_t *cov, utg_trans_t *o)
 
             if(return_flag != END_TIPS) continue;
 
-            flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, av[base_maxLen_i].v, av[i].v,
-            reverse_sources, &b_0, &b_1, ruIndex, miniedgeLen, 1);
+            flag = check_different_haps(g, ug, read_sg, av[base_maxLen_i].v, av[i].v,
+            reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, miniedgeLen, 1);
 
             // if((av[i].v>>1) == 255 && (av[base_maxLen_i].v>>1) == 33)
             // if((av[i].v>>1) == 1852 && (av[base_maxLen_i].v>>1) == 2441)
@@ -15616,10 +15504,8 @@ hap_cov_t *cov, utg_trans_t *o)
 
             if(ll>convexLen  && max_stop_baseLen>=ll*MAX_STOP_RATE)
             {
-                // flag = check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
-                // reverse_sources, &b_0, &b_1, ruIndex, miniedgeLen, stops_threshold);
-                flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
-                reverse_sources, &b_0, &b_1, ruIndex, miniedgeLen, stops_threshold);
+                flag = check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
+                reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, miniedgeLen, stops_threshold);
                 // #define UNAVAILABLE (uint32_t)-1
                 // #define PLOID 0
                 // #define NON_PLOID 1
@@ -15682,7 +15568,7 @@ hap_cov_t *cov, utg_trans_t *o)
 
 int detect_chimeric_by_topo(asg_t *g, ma_ug_t *ug, asg_t *read_sg, 
 ma_hit_t_alloc* reverse_sources, long long miniedgeLen, uint32_t stops_threshold, float drop_rate,
-R_to_U* ruIndex, utg_trans_t *o)
+R_to_U* ruIndex, utg_trans_t *o, uint8_t* is_r_het)
 {
     double startTime = Get_T();
     uint32_t i, k, v_i, v_beg, v_end, selfLen, w1, w2, wv, nw, n_vtx = g->n_seq * 2, n_reduced = 0, convex, convex_T, read_num;
@@ -15762,8 +15648,8 @@ R_to_U* ruIndex, utg_trans_t *o)
             }
             if(k != b_0.b.n) break;
             
-            if(/**check_different_haps_base**/check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1, 
-            ruIndex, miniedgeLen, stops_threshold)==PLOID)
+            if(check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1, 
+            ruIndex, is_r_het, miniedgeLen, stops_threshold)==PLOID)
             {
                 break;
             }
@@ -15808,8 +15694,8 @@ R_to_U* ruIndex, utg_trans_t *o)
             }
             if(k != b_0.b.n) break;
 
-            if(/**check_different_haps_base**/check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1, 
-            ruIndex, miniedgeLen, stops_threshold)==PLOID)
+            if(check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1, 
+            ruIndex, is_r_het, miniedgeLen, stops_threshold)==PLOID)
             {
                 break;
             }
@@ -16055,7 +15941,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
     /**********debug**********/
     if(just_bubble_pop == 0)
     {
-        cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, 2);
+        cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
     }
     /**********debug**********/
     long long pre_cons = get_graph_statistic(g);
@@ -16073,7 +15959,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
             asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL);            
             asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL);
             asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL);
-            detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL);
+            detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het);
             ///need consider tangles
             ///note we need both the read graph and the untig graph
         }
@@ -16082,14 +15968,14 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
     }
     if(just_bubble_pop == 0)
     {   
-        cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, 2);
+        cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
     }
 
     ///print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
 
-    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio);    
-    drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
-    all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres);
+    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, trio_flag, drop_ratio);    
+    drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het);
+    all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, cov->is_r_het, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres);
     if(is_first)
     {
         is_first = 0;
@@ -16133,24 +16019,23 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
     asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, NULL, 1);
     if(just_bubble_pop == 0)
     {
-        cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2);
+        cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
     }
     // print_debug_gfa(read_g, ug, coverage_cut, "debug_init", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
     long long pre_cons = get_graph_statistic(g);
     long long cur_cons = 0;
     while(pre_cons != cur_cons)
-    {    
+    {   
         pre_cons = get_graph_statistic(g);
         asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, NULL, 1);
         if(just_bubble_pop == 0)
         {
             ///need consider tangles  
-            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL);              
+            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL); 
             asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, NULL);   
             asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL);   
             asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL);
-            detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL);    
-
+            detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het);    
             if(round != T_ROUND)
             {
                 unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, 
@@ -16158,14 +16043,13 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
             }
         }
         cur_cons = get_graph_statistic(g);
-    }    
+    }
     if(just_bubble_pop == 0)
     {   
-        cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 
-        2);
+        cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
     }
-    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio); 
-    drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); 
+    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, (uint32_t)-1, drop_ratio); 
+    drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het); 
     // print_debug_gfa(read_g, ug, coverage_cut, "debug_clean_end", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
     unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1);    
     ///print_graph_statistic(g, "end");
@@ -16196,7 +16080,7 @@ int min_ovlp, hap_cov_t *cov)
 
     redo:
     asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, o, 1);
-    cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2);
+    cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
 
     long long pre_cons = get_graph_statistic(g);
     long long cur_cons = 0;
@@ -16214,7 +16098,7 @@ int min_ovlp, hap_cov_t *cov)
                 asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, o);   
                 asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, o);   
                 asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, o);
-                detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, o);    
+                detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, o, cov->is_r_het);    
                 
                 cur_cons = get_graph_statistic(g);
             }    
@@ -16236,11 +16120,10 @@ int min_ovlp, hap_cov_t *cov)
     }
 
     
-    cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 
-    2);
+    cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
     
-    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio); 
-    drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); 
+    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, (uint32_t)-1, drop_ratio); 
+    drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het); 
     print_debug_gfa(read_g, ug, coverage_cut, "debug_clean_end", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
     unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1);    
     
@@ -16278,7 +16161,7 @@ void set_drop_trio_flag(ma_ug_t *ug)
 }
 
 void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
-ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, 
+ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het,
 uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
 {
     asg_t* nsg = ug->g;
@@ -16286,7 +16169,7 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
     ma_utg_t *u;
     uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t));
 
-    drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex);
+    drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex, is_r_het);
 
     while (n_reduce)
     {
@@ -16365,7 +16248,7 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
         }
     }
 
-    drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex);
+    drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex, is_r_het);
     asg_cleanup(nsg);
     free(primary_flag);
 }
@@ -16496,9 +16379,9 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
 }
 
 uint32_t unitig_simi(uint32_t x, uint32_t y, ma_ug_t* ug, ma_hit_t_alloc* reverse_sources, 
-R_to_U* ruIndex)
+R_to_U* ruIndex, uint8_t* is_r_het)
 {
-    uint32_t k, j, uId, tn, is_Unitig, rId, ref_unitig, min_count, max_count;
+    uint32_t k, j, uId, tn, is_Unitig, rId, ref_unitig, min_count, max_count, n_het, n_hom;
     ma_utg_t *nsu_x = NULL, *nsu_y = NULL, *nsu_query = NULL;
     nsu_x = &(ug->u.a[x]);
     nsu_y = &(ug->u.a[y]);
@@ -16516,11 +16399,13 @@ R_to_U* ruIndex)
         ref_unitig = x;
     }
      
-    min_count = max_count = 0;
+    min_count = max_count = n_het = n_hom = 0;
     for (k = 0; k < nsu_query->n; k++)
     {
         rId = nsu_query->a[k]>>33;
         if(reverse_sources[rId].length >= 0) min_count++;
+        if((is_r_het[rId] & C_HET) || (is_r_het[rId] & P_HET)) n_het++;
+        n_hom++;
 
         for (j = 0; j < reverse_sources[rId].length; j++)
         {
@@ -16539,7 +16424,7 @@ R_to_U* ruIndex)
     }
     
     if(min_count == 0) return UNAVAILABLE;
-	if(max_count > min_count*asm_opt.purge_simi_thres/**DIFF_HAP_RATE**/) return PLOID;
+    if(max_count > min_count*asm_opt.purge_simi_thres && n_het >= n_hom*HET_HOM_RATE) return PLOID;
 	return NON_PLOID;
 }
 
@@ -16573,7 +16458,7 @@ uint32_t* non_require, uint32_t* ambigious)
 ///note: to use this function, don't renew unitig graph!!!!!!!!!
 void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
 ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, 
-ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ)
+ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het, float double_check_rate, int non_tig_occ)
 {
     
 
@@ -16679,7 +16564,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int n
                     if(flag_occ <= ((non_flag_occ+flag_occ)*drop_rate)) continue;
                     if((flag_occ+non_flag_occ) == 0) continue;
                 }**/
-                if(unitig_simi(uId, (uint32_t)(u_vecs.a.a[k]), ug, reverse_sources, ruIndex)==PLOID)
+                if(unitig_simi(uId, (uint32_t)(u_vecs.a.a[k]), ug, reverse_sources, ruIndex, is_r_het)==PLOID)
                 {
                     break;
                 }
@@ -16810,7 +16695,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex)
 
 
 
-void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
+void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het)
 {
     uint32_t v, n_vtx = nsg->n_seq*2, convex_f, convex_b, i, nv;
     long long ll, tmp, max_stop_nodeLen, max_stop_baseLen;
@@ -16845,8 +16730,8 @@ void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* re
             }
             get_real_length(nsg, convex_f, &convex_f);
             if(convex_f != convex_b) continue;
-            if(/**check_different_haps_base**/check_different_haps(nsg, ug, read_g, v^1, av[i].v, 
-            reverse_sources, &b_0, &b_1, ruIndex, 2, 1) == PLOID)
+            if(check_different_haps(nsg, ug, read_g, v^1, av[i].v, 
+            reverse_sources, &b_0, &b_1, ruIndex, is_r_het, 2, 1) == PLOID)
             {
                 av[i].del = 1;
 		        asg_arc_del(nsg, av[i].v^1, v^1, 1);
@@ -16987,10 +16872,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
     uint32_t v, n_vtx = nsg->n_seq;    
     hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, 
                                                     coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_trio>0?1:0);
-    if(cov->t_ch)
-    {
-        set_r_het_flag(*ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
-    }
+    if(cov->t_ch) cov->t_ch->ir_het = cov->is_r_het;
 
     if(asm_opt.recover_atg_cov_min == -1024)
     {
@@ -17011,9 +16893,9 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
     }
     
     
-    adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t);
+    adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t, cov->is_r_het);
     ///primary_flag = get_utg_attributes(*ug, read_g, coverage_cut, sources, ruIndex);
-    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, 
+    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0, 
     DOUBLE_CHECK_THRES, flag, drop_rate);
 
     nsg = (*ug)->g;
@@ -17033,7 +16915,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
 
     update_hap_label(*ug, read_g);
 
-    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, 
+    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0, 
     DOUBLE_CHECK_THRES, flag, drop_rate);
 
     force_trio_clean((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, flag, 0.55, 0.01, 5);
@@ -17059,7 +16941,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
     ///if(flag == MOTHER) print_untig_by_read(*ug, "m64043_200627_000137/124716590/ccs", 2789716, NULL, NULL, "beg");
     
     
-    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 1, 
+    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 1, 
     FINAL_DOUBLE_CHECK_THRES, flag, drop_rate);
 
     update_hap_label(NULL, read_g);
@@ -17730,7 +17612,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
                             ori = t_ch->b_buf_0.b.a[i]&1;
                             for (k = 0; k < p->n; k++)
                             {
-                                t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
+                                t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
                             }
                         }
                         
@@ -17744,7 +17626,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
                             ori = t_ch->b_buf_1.b.a[i]&1;
                             for (k = 0; k < p->n; k++)
                             {
-                                t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
+                                t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
                             }
                         }
 
@@ -17776,7 +17658,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
             ori = cov->t_ch->topo_res.a[i]&1;
             for (k = 0; k < p->n; k++)
             {
-                t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
+                t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
             }
             /***********************x***********************/
 
@@ -17789,7 +17671,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
                 ori = t_ch->b_buf_0.b.a[k_i]&1;
                 for (k = 0; k < p->n; k++)
                 {
-                    t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
+                    t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
                 }
             }
             /***********************y***********************/
@@ -20382,7 +20264,7 @@ uint32_t type)
 inline uint32_t walk_through(asg_t *read_g, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
 long long maxShortUntig, float l_untig_rate, float max_node_threshold, buf_t* b_0, buf_t* b_1,
 kvec_t_u32_warp* u_vecs, uint8_t* visit, uint32_t v, uint32_t* r_beg, uint32_t* r_end, 
-uint32_t* r_next_uID, R_to_U* ruIndex)
+uint32_t* r_next_uID, R_to_U* ruIndex, uint8_t* is_r_het)
 {
     (*r_beg) = (*r_end) = (uint32_t)-1;
     asg_t* nsg = ug->g;
@@ -20478,8 +20360,8 @@ uint32_t* r_next_uID, R_to_U* ruIndex)
             // #define UNAVAILABLE (uint32_t)-1
             // #define PLOID 0
             // #define NON_PLOID 1
-            if(returnFlag == 1 && /**check_different_haps_base**/check_different_haps(nsg, ug, read_g, beg, next_uID, 
-                reverse_sources, b_0, b_1, ruIndex, minLongUntig-1, 1) == PLOID)
+            if(returnFlag == 1 && check_different_haps(nsg, ug, read_g, beg, next_uID, 
+                reverse_sources, b_0, b_1, ruIndex, is_r_het, minLongUntig-1, 1) == PLOID)
             {
                 ///output_tangles(beg, next_uID, u_vecs->a.a, u_vecs->a.n, (char*)("???"));
                 returnFlag = 0;
@@ -22858,8 +22740,8 @@ float drop_ratio)
 }
 
 void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
-long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t trio_flag, 
-float drop_ratio)
+long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint8_t* is_r_het, 
+uint32_t trio_flag, float drop_ratio)
 {
     buf_t b_0, b_1;
     memset(&b_0, 0, sizeof(buf_t));
@@ -22910,7 +22792,7 @@ float drop_ratio)
             {
                 flag = walk_through(read_g, ug, reverse_sources, minLongUntig, 
                     maxShortUntig, l_untig_rate, max_node_threshold, &b_0, &b_1, 
-                    &u_vecs, visit, sv, &beg, &end, &next_uID, ruIndex);
+                    &u_vecs, visit, sv, &beg, &end, &next_uID, ruIndex, is_r_het);
                 n_reduce += flag;
                 if(flag != UNROLL_M)
                 {
@@ -23604,7 +23486,7 @@ int get_arc_t(Edge_iter* x, asg_arc_t* get)
 }
 
 
-void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t, double dupLenThres)
+void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t, uint8_t* is_r_het, double dupLenThres)
 {
     asg_t* nsg = ug->g;
     uint32_t v, n_vtx = nsg->n_seq * 2, rnw, nw, beg, end, i;
@@ -23752,8 +23634,8 @@ void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reve
                     continue;
                 }
 
-                if(/**check_different_haps_base**/check_different_haps(nsg, ug, read_g, beg, end, reverse_sources, &b_0, &b_1,
-                    ruIndex, 2, 1) == PLOID)
+                if(check_different_haps(nsg, ug, read_g, beg, end, reverse_sources, &b_0, &b_1,
+                    ruIndex, is_r_het, 2, 1) == PLOID)
                 {
                     continue;
                 }
@@ -23806,12 +23688,12 @@ void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reve
 
 
 
-void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t)
+void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t, uint8_t* is_r_het)
 {
     double startTime = Get_T();
     asg_t* nsg = ug->g;
-    unroll_simple_case_advance(ug, sg, reverse_sources, ruIndex, b_mask_t, 2.5);
-    drop_semi_circle(ug, ug->g, sg, reverse_sources, ruIndex);
+    unroll_simple_case_advance(ug, sg, reverse_sources, ruIndex, b_mask_t, is_r_het, 2.5);
+    drop_semi_circle(ug, ug->g, sg, reverse_sources, ruIndex, is_r_het);
     asg_cleanup(nsg);
     asg_symm(nsg);
     ///debug_utg_graph(ug, sg, 0, 0);
@@ -24455,6 +24337,8 @@ uint32_t is_collect_trans)
     }
 
     if(set) free(set);
+    CALLOC(x->is_r_het, read_g->n_seq);
+    set_r_het_flag(ug, read_g, coverage_cut, sources, ruIndex, x->is_r_het);
 
     x->t_ch = NULL;
     if(is_collect_trans) x->t_ch = init_trans_chain(ug, read_g->n_seq);
@@ -24468,6 +24352,7 @@ void destory_hap_cov_t(hap_cov_t **x)
     {
         free((*x)->cov);
         free((*x)->pos_idx);
+        free((*x)->is_r_het);
         kv_destroy((*x)->u_buffer.a);
         kv_destroy((*x)->tailIndex.a);
         kv_destroy((*x)->prevIndex.a);
@@ -24514,7 +24399,7 @@ void reset_trans_chain(trans_chain* t_ch, ma_utg_t *u)
 {
     uint32_t k = 0;
     if(u->n == 0 || u->m == 0) return;
-    for (k = 0; k < u->n; k++) t_ch->is_r_het[u->a[k]>>33] = N_HET;
+    for (k = 0; k < u->n; k++) t_ch->ir_het[u->a[k]>>33] = N_HET;
 }
 
 void append_utg(ma_ug_t* ptg, ma_ug_t* atg, trans_chain* t_ch)
@@ -24733,8 +24618,8 @@ uint32_t collect_p_trans, uint32_t collect_p_trans_f)
     ma_utg_t* u = NULL;
     hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, 
                             coverage_cut, max_hang, min_ovlp, (asm_opt.purge_level_primary>0||i_cov)?1:0);
-    if(cov->t_ch) set_r_het_flag(*ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
-    adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t);
+    if(cov->t_ch) cov->t_ch->ir_het = cov->is_r_het;
+    adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t, cov->is_r_het);
 
     nsg = (*ug)->g;
     n_vtx = nsg->n_seq;
@@ -28170,7 +28055,7 @@ void reset_bub(bubble_type* bub, ma_ug_t *ug, trans_chain* back_ug_chain, kvec_a
 
     new_rtg_edges->a.n = 0;
     ///classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp);
-    identify_bubbles(ug, bub, back_ug_chain->is_r_het, NULL);
+    identify_bubbles(ug, bub, back_ug_chain->ir_het, NULL);
     // update_bubble_chain(ug, bub, 0, 1);
     // resolve_bubble_chain_tangle(ug, bub);
     // fprintf(stderr, "bub.f_bub: %lu, bub.b_bub: %lu, bub.b_end_bub: %lu, bub.tangle_bub: %lu, bub.cross_bub: %lu\n", 
@@ -28461,7 +28346,7 @@ int max_hang, int min_ovlp, bubble_type* bub, long long gap_fuzz)
 }
 
 
-uint8_t *rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
+void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, 
 float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, 
 bub_label_t* b_mask_t)
@@ -28505,14 +28390,11 @@ bub_label_t* b_mask_t)
         rescue_missing_hap_ovlp(ug, sg, sources, coverage_cut, max_hang, min_ovlp, &bub, gap_fuzz);
     }
 
-    uint8_t *het_flag = cov->t_ch->is_r_het;
-    cov->t_ch->is_r_het = NULL;
     destory_bubbles(&bub);
     destory_hap_cov_t(&cov);
     ma_ug_destroy(ug);
     kv_destroy(new_rtg_edges.a);
     ma_ug_destroy(copy_ug); copy_ug = NULL;
-    return het_flag;
 }
 
 void update_unitig(long long step, long long init, ma_utg_t* nsu, asg_t *r_g, 
@@ -30978,12 +30860,12 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, 
     max_hang_length, mini_overlap_length);
 
-    ruIndex->is_het = rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, 
+    rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, 
     ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz, &b_mask_t);
 
     output_unitig_graph(sg, coverage_cut, o_file, sources, ruIndex, max_hang_length, mini_overlap_length);
     // flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
-    flat_soma_v(sg, sources, ruIndex); free(ruIndex->is_het); ruIndex->is_het = NULL;
+    flat_soma_v(sg, sources, ruIndex);
 
     output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, 
         asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length);
diff --git a/Overlaps.h b/Overlaps.h
index dc690ad..b4eec32 100644
--- a/Overlaps.h
+++ b/Overlaps.h
@@ -715,6 +715,7 @@ uint32_t stops_threshold, buf_t* b)
 #define TRIO_DROP_LENGTH_THRES 0.8
 #define MAX_STOP_RATE 0.6
 #define TANGLE_MISSED_THRES 0.6
+#define HET_HOM_RATE 0.7
 
 typedef struct {
 	uint32_t father_occ;
@@ -747,17 +748,17 @@ typedef struct{
 } bub_label_t;
 
 void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
-long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t trio_flag, 
-float drop_ratio);
-void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t);
+long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint8_t* is_r_het, 
+uint32_t trio_flag, float drop_ratio);
+void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t, uint8_t* is_r_het);
 void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
 R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t chainLenThres, uint32_t is_bubble_check, 
 uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t);
 void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
 R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
 void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
-ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ);
-void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het, float double_check_rate, int non_tig_occ);
+void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het);
 void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges, bub_label_t* b_mask_t);
 void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
@@ -822,7 +823,7 @@ typedef struct {
 typedef struct{
 	uint32_t* rUidx;
 	uint64_t* rUpos;
-	uint8_t* is_r_het;
+	uint8_t* ir_het;
 	uint32_t r_num, u_num;
 	kvec_t(bed_in) bed;
 	kvec_t(uint32_t) topo_buf;
@@ -848,7 +849,7 @@ typedef struct {
     kvec_asg_arc_t_offset u_buffer; 
     kvec_t_i32_warp tailIndex;
     kvec_t_i32_warp prevIndex;
-	///hc_links* link;
+	uint8_t* is_r_het;
 	trans_chain* t_ch;
 }hap_cov_t;
 
diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp
index 00c1887..e84bf4e 100644
--- a/Purge_Dups.cpp
+++ b/Purge_Dups.cpp
@@ -176,6 +176,7 @@ long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* cove
 		///if the second peak is not significant
 		if(max2 < max * 0.05 || min > max2 * 0.95) max2 = max2_i = -1;
 	}
+    if(max2 < max*0.0075) max2 = max2_i = -1;
 		
 
     // look for smaller peak on the high end
@@ -205,7 +206,7 @@ long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* cove
 			
 		if (max3 < max * 0.05 || min > max3 * 0.95 || max3_i > max_i * 3) max3 = max3_i = -1;
 	}
-
+    if (max3 < max*0.0075) max3 = max3_i = -1;
 
     
 
@@ -216,6 +217,13 @@ long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* cove
         topo_peak = cov_buf[topo_peak_i];
         if (topo_peak <= max * 0.05) topo_peak_i = topo_peak = -1;
     }
+
+    if(asm_opt.purge_level_primary == 0) 
+    {
+        (*hom_peak) = max_i;
+        return;
+    }
+    
     
     long long k_mer_het, k_mer_hom, coverage_het, coverage_hom, alter_peak;
     k_mer_het = k_mer_hom = coverage_het = coverage_hom = alter_peak = -1;
@@ -345,12 +353,6 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
             }
         }
 
-        // if(qn == 1893151 || qn == 1929038)
-        // {
-        //     fprintf(stderr, "qn: %lu, C_bases_primary: %lld, C_bases_alter: %lld, C_bases: %lld\n",
-        //     qn, C_bases_primary, C_bases_alter, C_bases);
-        // }
-
         C_bases = C_bases_primary + C_bases_alter;
         if(C_bases_alter < C_bases * ALTER_COV_THRES) continue;
 
@@ -369,8 +371,6 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
         }
     }
 
-    ///fprintf(stderr, "alter max_i: %lld, max: %lld\n", max_i, max);
-    ///if(max_i < 5) max_i = max = -1;
     if(alter_bases < primary_bases * REAL_ALTER_THRES) max_i = max = -1;
 
     free(cov_buf);
@@ -4981,9 +4981,9 @@ void debug_p_g_t(p_g_t* pg, hap_cov_t *cov, asg_t *read_g)
             offset += (uint32_t)u->a[i];
             if(i >= sid && i <= eid)
             {
-                if(cov->t_ch->is_r_het[u->a[i]>>33] != t->h_status)
+                if(cov->t_ch->ir_het[u->a[i]>>33] != t->h_status)
                 {
-                    fprintf(stderr, "ERROR-(-3): is_r_het: %u, h_status: %u\n", cov->t_ch->is_r_het[u->a[i]>>33], t->h_status);
+                    fprintf(stderr, "ERROR-(-3): is_r_het: %u, h_status: %u\n", cov->t_ch->ir_het[u->a[i]>>33], t->h_status);
                 }
             }
         }
@@ -5065,11 +5065,11 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
         for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) 
         {   
             ///if (k == u->n || (!!cov->t_ch->is_r_het[u->a[k]>>33]) != (!!cov->t_ch->is_r_het[u->a[l]>>33]))
-            if (k == u->n || cov->t_ch->is_r_het[u->a[k]>>33] != cov->t_ch->is_r_het[u->a[l]>>33])
+            if (k == u->n || cov->t_ch->ir_het[u->a[k]>>33] != cov->t_ch->ir_het[u->a[l]>>33])
             {
                 kv_pushp(p_node_t, pg->pg_het_node, &t);
                 t->c_ug_id = uId;
-                t->h_status = cov->t_ch->is_r_het[u->a[l]>>33];
+                t->h_status = cov->t_ch->ir_het[u->a[l]>>33];
                 t->baseBeg = l_pos;
                 t->baseEnd = offset + read_g->seq[u->a[k-1]>>33].len - 1;
                 t->nodeBeg = l;
@@ -5261,7 +5261,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
     init_hap_overlaps_list(&all_ovlp, nsg->n_seq);
     hap_overlaps_list back_all_ovlp;
     init_hap_overlaps_list(&back_all_ovlp, nsg->n_seq);
-    ///uint32_t junk_cov, hap_cov, dip_cov, junk_occ, repeat_occ, single_cov;
     asg_arc_t t, *p = NULL;
     int r;
     hap_alignment_struct_pip hap_buf;
@@ -5382,12 +5381,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
                 r = get_hap_arch(&(all_ovlp.x[uId].a.a[i]), ug->u.a[all_ovlp.x[uId].a.a[i].xUid].len, 
                 ug->u.a[all_ovlp.x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
                 
-                // if(all_ovlp.x[uId].a.a[i].xUid == 118 && all_ovlp.x[uId].a.a[i].yUid == 82)
-                // {
-                //     fprintf(stderr, "r: %d\n", r);
-                //     print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
-                // }
-                
                 if(r < 0) continue;
                 p = asg_arc_pushp(pg->pg_h_lev);
                 *p = t;
diff --git a/Purge_Dups.h b/Purge_Dups.h
index 55107c4..7b16bd2 100644
--- a/Purge_Dups.h
+++ b/Purge_Dups.h
@@ -10,7 +10,7 @@
 #define HOM_PEAK_RATE 1.25
 #define HET_PEAK_RATE (HOM_PEAK_RATE*2)
 #define ALTER_COV_THRES 0.9
-#define REAL_ALTER_THRES 0.1
+#define REAL_ALTER_THRES 0.25
 #define CHAIN_FILTER_RATE 0.7
 
 #define SELF_EXIST 0
diff --git a/anchor.cpp b/anchor.cpp
index 072abad..08b9c31 100644
--- a/anchor.cpp
+++ b/anchor.cpp
@@ -1,4 +1,5 @@
 #include <stdio.h>
+#include <math.h>
 #include "htab.h"
 #include "ksort.h"
 #include "Hash_Table.h"
@@ -7,8 +8,9 @@
 
 typedef struct { // this struct is not strictly necessary; we can use k_mer_pos instead, with modifications
 	uint64_t srt;
-	uint32_t self_off:31, good:1;
+	uint32_t self_off;
 	uint32_t other_off;
+	uint32_t cnt;
 } anchor1_t;
 
 #define an_key1(a) ((a).srt)
@@ -23,7 +25,7 @@ KSORT_INIT(or_xs, overlap_region, oreg_xs_lt)
 KSORT_INIT(or_ss, overlap_region, oreg_ss_lt)
 
 typedef struct {
-	int n, good;
+	int n;
 	const ha_idxpos_t *a;
 } seed1_t;
 
@@ -58,12 +60,13 @@ int ha_ov_type(const overlap_region *r, uint32_t len)
 }
 
 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)
+						   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, st_mt_t *sp)
 {
 	uint32_t i, rlen;
 	uint64_t k, l;
-	double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
-	double high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
+	uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
+	uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
+	if(low_occ < 2) low_occ = 2;
 
 	// prepare
     clear_Candidates_list(cl);
@@ -73,7 +76,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, NULL, -1, asm_opt.dp_min_len, -1, NULL);
+	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, sp, asm_opt.mz_rewin);
 	// minimizer of queried read
 	if (ab->mz.m > ab->old_mz_m) {
 		ab->old_mz_m = ab->mz.m;
@@ -83,7 +86,6 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		int n;
 		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
 		ab->seed[i].n = n;
-		ab->seed[i].good = (n > low_occ && n < high_occ);
 		ab->n_a += n;
 	}
 	if (ab->n_a > ab->m_a) {
@@ -102,7 +104,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 			uint8_t rev = z->rev == y->rev? 0 : 1;
 			an->other_off = y->pos;
 			an->self_off = rev? ucr->length - 1 - (z->pos + 1 - z->span) : z->pos;
-			an->good = s->good;
+			an->cnt = s->n;
 			an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
 		}
 	}
@@ -129,7 +131,16 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		p->strand = ab->a[k].srt >> 32 & 1;
 		p->offset = ab->a[k].other_off;
 		p->self_offset = ab->a[k].self_off;
-		p->good = ab->a[k].good;
+		if(ab->a[k].cnt > low_occ && ab->a[k].cnt < high_occ){
+			p->cnt = 1;
+		}
+		else if(ab->a[k].cnt <= low_occ){
+			p->cnt = 2;
+		}
+		else{
+			p->cnt = 1 + ((ab->a[k].cnt + (high_occ<<1) - 1)/(high_occ<<1));
+			p->cnt = pow(p->cnt, 1.1);
+		}
 	}
 	cl->length = ab->n_a;
 
@@ -269,7 +280,6 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
 		int n;
 		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
 		ab->seed[i].n = n;
-		ab->seed[i].good = 0;
 		ab->n_a += n;
 	}
 	if (ab->n_a > ab->m_a) {
@@ -288,7 +298,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
 			uint8_t rev = z->rev == y->rev? 0 : 1;
 			an->other_off = y->pos;
 			an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos;
-			an->good = s->good;
+			an->cnt = 1;
 			an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
 		}
 	}
@@ -315,7 +325,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
 		p->strand = ab->a[k].srt >> 32 & 1;
 		p->offset = ab->a[k].other_off;
 		p->self_offset = ab->a[k].self_off;
-		p->good = ab->a[k].good;
+		p->cnt = 1;
 	}
 	cl->length = ab->n_a;
 
@@ -391,14 +401,14 @@ 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)
+								 kvec_t_u64_warp* dbg_ct, st_mt_t *sp)
 {
 	extern void *ha_flt_tab;
 	extern ha_pt_t *ha_idx;
 	extern void *ha_flt_tab_hp;
 	extern ha_pt_t *ha_idx_hp;
 
-	ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct);
+	ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct, sp);
 
 	if(ha_idx_hp)
 	{
@@ -428,7 +438,7 @@ void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overl
 		overlap_list->length = k;
 
 
-		ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar, dbg_ct);	
+		ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar, dbg_ct, sp);	
 		
 		if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
 		{
diff --git a/hic.cpp b/hic.cpp
index ab44a93..309a741 100644
--- a/hic.cpp
+++ b/hic.cpp
@@ -14786,7 +14786,7 @@ void debug_gfa_space(ha_ug_index* idx, ma_ug_t* ug, trans_chain* t_ch, kv_u_tran
     memset(&bub, 0, sizeof(bubble_type));
     bub.round_id = 0; bub.n_round = 2;
 
-    identify_bubbles(ug, &bub, t_ch->is_r_het, ref);
+    identify_bubbles(ug, &bub, t_ch->ir_het, ref);
 
     hc_links link;
     init_hc_links(&link, ug->g->n_seq, t_ch);
@@ -15687,7 +15687,7 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits,
     pe_hit *h_a = NULL;
     u_trans_t *p = NULL;
 
-    identify_bubbles(idx->ug, bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans));
+    identify_bubbles(idx->ug, bub, idx->t_ch->ir_het, &(idx->t_ch->k_trans));
 
     ta->idx.n = ta->n = 0;
     if(hits->idx.n == 0) idx_hc_links(hits, idx, NULL);
diff --git a/htab.cpp b/htab.cpp
index 025bf24..04c02e4 100644
--- a/htab.cpp
+++ b/htab.cpp
@@ -675,7 +675,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;
 	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->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[tid]), asm_opt.mz_rewin);
 	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));
@@ -807,17 +807,16 @@ 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);						
+			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->mt[i].a);
 				free(s->mz_buf[i].a);
 			}
-			if(s->mt) free(s->mt);
+			free(s->mt);
 			free(s->mz_buf);
 			// insert minimizers
 			if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {///insert whole minimizer
diff --git a/htab.h b/htab.h
index 0e4f9c3..9e2caa0 100644
--- a/htab.h
+++ b/htab.h
@@ -101,7 +101,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, 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);
+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, int32_t ws);
 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);
diff --git a/rcut.cpp b/rcut.cpp
index dbe80a5..22932e5 100644
--- a/rcut.cpp
+++ b/rcut.cpp
@@ -889,9 +889,9 @@ trans_chain* t_ch)
             offset += (uint32_t)u->a[i];
             if(i >= sid && i <= eid)
             {
-                if(t_ch->is_r_het[u->a[i]>>33] != t->hs)
+                if(t_ch->ir_het[u->a[i]>>33] != t->hs)
                 {
-                    fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", t_ch->is_r_het[u->a[i]>>33], t->hs);
+                    fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", t_ch->ir_het[u->a[i]>>33], t->hs);
                 }
             }
         }
@@ -939,11 +939,11 @@ void update_mc_edges(mc_g_t *mg, hap_overlaps_list* ha, kv_u_trans_t *ta, trans_
 			u = &(mg->ug->u.a[v]);
 			for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) 
 			{   
-				if (k == u->n || t_ch->is_r_het[u->a[k]>>33] != t_ch->is_r_het[u->a[l]>>33])
+				if (k == u->n || t_ch->ir_het[u->a[k]>>33] != t_ch->ir_het[u->a[l]>>33])
 				{
 					kv_pushp(mc_interval_t, p, &t);
 					t->uID = v;
-					t->hs = t_ch->is_r_het[u->a[l]>>33];
+					t->hs = t_ch->ir_het[u->a[l]>>33];
 
 					t->bS = l_pos;
 					t->bE = offset + mg->rg->seq[u->a[k-1]>>33].len - 1;
@@ -2680,10 +2680,10 @@ void set_p_flag(mc_g_t *mg, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch,
 	for (i = 0; i < u->n; i++)
 	{
 		trio_flag[u->a[i]>>33] |= SET_TRIO;
-		if(t_ch->is_r_het[u->a[i]>>33] == N_HET) continue;
+		if(t_ch->ir_het[u->a[i]>>33] == N_HET) continue;
 		if(s == 0)
 		{
-			if(t_ch->is_r_het[u->a[i]>>33]&P_HET)//special case
+			if(t_ch->ir_het[u->a[i]>>33]&P_HET)//special case
 			{
 				trio_flag[u->a[i]>>33] |= FATHER;
 			}
diff --git a/sketch.cpp b/sketch.cpp
index 7bca350..c9a92ab 100644
--- a/sketch.cpp
+++ b/sketch.cpp
@@ -229,7 +229,7 @@ inline void hf_select(ha_mz1_v *p, int32_t si, int32_t ei, int32_t n, int32_t le
 			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)
+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, nw[2], min_len;
 	ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation
@@ -289,6 +289,288 @@ static void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
 	p->n = n;
 }
 
+static inline int mzcmp_l(const ha_mz1_v *p, int32_t ai, int32_t bi)
+{
+	if(ai >= 0 && bi >= 0){
+		ha_mz1_t *a = &(p->a[ai]), *b = &(p->a[bi]);
+		if(a->rid > 0 && b->rid > 0) return mzcmp(a, b);
+		return (a->rid == 0) - (b->rid == 0);
+	}
+	return (ai < 0) - (bi < 0);
+}
+
+#define GL(x, i) ((int64_t)((uint32_t)((x).a[(i)])))
+#define A_M(p, i) ((i) >= 0 && (p).a[(i)].rid > 0)
+int32_t qfw(ha_mz1_v *p, st_mt_t *mt, int32_t n, int32_t tot_l, int32_t ws, int32_t i, int32_t *mi)
+{
+	int32_t m, si;
+	for (si = i, (*mi) = -1; i < n; i++){
+		if(GL(*mt, i) >= ws || (i+1 < n && GL(*mt, i) < ws && GL(*mt, i+1) > ws) || 
+				(i+1 == n && tot_l >= ws && GL(*mt, i) < ws)){
+			for (m = si; m <= i; m++){
+				if(!A_M(*p, m)) continue;
+				if(mzcmp_l(p, *mi, m) >= 0) (*mi) = m;
+			}
+			if((*mi) >= 0 && A_M(*p, *mi)){
+				for (m = si; m <= i; m++){
+					if(!A_M(*p, m)) continue;
+					if(mzcmp_l(p, *mi, m) == 0) mt->a[m] |= 0x100000000;
+				}
+			}
+			break;
+		}
+	}
+	return i;
+}
+
+void dbg_boundary(ha_mz1_v *p, st_mt_t *mt, int32_t w, int32_t k, int32_t tot_l)
+{
+	if(tot_l < w + k -1) return;
+	int32_t i, m, n = p->n, s, a;
+	for (i = 0; i < n; i++){
+		if(GL(*mt, i) >= w+k-1){
+			for (m = s = a = 0; m <= i; m++){
+				if(!A_M(*p, m)) continue;
+				if(GL(*mt, m) <= w+k-1){
+					a++;
+					if(mt->a[m]&0x100000000) s++;
+				}
+			}
+			if(a > 0 && s == 0){
+				fprintf(stderr, "\nERROR1, s: %d, n: %d, tot_l: %d, end_l: %ld\n", s, n, tot_l, GL(*mt, i));
+				for (m = s = a = 0; m <= i; m++){
+					if(!A_M(*p, m)) continue;
+					if(GL(*mt, m) <= w+k-1){
+						fprintf(stderr, "lp: %ld\n", GL(*mt, m));
+						a++;
+						if(mt->a[m]&0x100000000) s++;
+					}
+				}
+				
+			} 
+			
+			break;
+		}
+	}
+	if(i == n){
+		for (m = s = a = 0; m < n; m++){
+			if(!A_M(*p, m)) continue;
+			if(GL(*mt, m) <= w+k-1){
+				a++;
+				if(mt->a[m]&0x100000000) s++;
+			}
+		}
+		if(a > 0 && s == 0) fprintf(stderr, "ERROR2\n");
+	}
+
+	for (i = n-1; i >= 0; i--)
+	{
+		if (GL(*mt, i) + w <= tot_l + 1) {
+			for (m = i, s = a = 0; m < n; m++){
+				if(!A_M(*p, m)) continue;
+				if(GL(*mt, m) + w >= tot_l + 1){
+					a++;
+					if(mt->a[m]&0x100000000) s++;
+				} 				
+			}
+			if(a > 0 && s == 0) {
+				fprintf(stderr, "\nERROR3, s: %d, n: %d, tot_l: %d, end_l: %ld\n", s, n, tot_l, GL(*mt, i));
+				for (m = i, s = a = 0; m < n; m++){
+					if(!A_M(*p, m)) continue;
+					if(GL(*mt, m) + w >= tot_l + 1){
+						fprintf(stderr, "lp: %ld\n", GL(*mt, m));
+						a++;
+						if(mt->a[m]&0x100000000) s++;
+					} 				
+				}
+			}
+			
+			break;
+		}
+	}
+	
+	if(i < 0){
+		for (m = s = a = 0; m < n; m++){
+			if(!A_M(*p, m)) continue;
+			if(GL(*mt, m) + w >= tot_l + 1){
+				a++;
+				if(mt->a[m]&0x100000000) s++;
+			} 				
+		}
+		if(a > 0 && s == 0) fprintf(stderr, "ERROR4\n");
+	}
+}
+
+static void select_mz_h(ha_mz1_v *p, st_mt_t *mt, int len, int sample_dist, int32_t w, int32_t k, int32_t tot_l)
+{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak
+	int32_t i, mi = -1, si, last0 = -1, n = (int32_t)p->n, m = 0, ws = w + k - 1;
+
+	if (n == 0) 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 = m = 0, last0 = -1; 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;
+				if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){
+					last0 = -2;
+					m++;
+					break;
+				}
+			}
+			last0 = i;
+		}
+	}
+	if (m == 0) return; // no high-frequency k-mers; do nothing
+	if(last0 >= -1) goto ff;
+	i = 0;
+	i = qfw(p, mt, n, tot_l, ws, i, &mi);
+	if(i == n) goto ff;
+	
+	for (si = 0, i++; i < n; i++){
+		for (; si < i; si++){
+			if(GL(*mt, si) + w > GL(*mt, i)) break;
+		}
+
+		// a new minimum; then write the old min
+		if(mzcmp_l(p, i, mi) <= 0) {
+			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
+			mi = i;
+		}// old min has moved outside the window
+		else if(si > mi){
+			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
+			for (m = si, mi = -1; m <= i; m++){
+				if(mzcmp_l(p, mi, m) >= 0) mi = m;
+            }
+			if(A_M(*p, mi)){
+				for (m = si; m <= i; m++){
+					if(!A_M(*p, m)) continue;
+					if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;
+				}
+			}
+		}
+	}	
+
+	if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
+
+	for (i = n - 1; si < n && GL(*mt, si) + w <= tot_l + 1; si++){
+		if(si > mi){
+			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
+			for (m = si, mi = -1; m <= i; m++){
+				if(mzcmp_l(p, mi, m) >= 0) mi = m;
+            }
+			if(A_M(*p, mi)){
+				for (m = si; m <= i; m++){
+					if(!A_M(*p, m)) continue;
+					if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;
+				}
+			}
+		}
+	}
+	
+	/**
+	dbg_boundary(p, mt, w, k, tot_l);
+	fprintf(stderr, "\n");
+	for (i = 0; i < (int32_t)p->n; ++i){
+		if(p->a[i].rid == 0) continue;
+		fprintf(stderr, "%cl: %u, pos: %lu, cnt: %lu, key: %lu, i: %d\n", "+-"[!!(mt->a[i]&0x100000000)], 
+													(uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid, p->a[i].x, i);
+		// if (mt->a[i]&0x100000000){
+		// 	fprintf(stderr, "+l: %u, pos: %lu, cnt: %lu\n", (uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid);
+		// }
+	}
+	**/
+	ha_mz1_t b[MAX_MAX_HIGH_OCC];
+	for (i = 0, last0 = -1; 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;
+				if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){
+					for (m = last0 + 1, mi = 0; m < i; ++m){
+						if(mt->a[m]&0x100000000) p->a[m].rid = 0, mi++;
+					}
+					if(mi == 0) hf_select(p, last0, i, n, len, sample_dist, b, 0);
+				} 
+			}
+			last0 = i;
+		}
+	}
+
+	ff:
+	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;
+}
+
+void debug_pl(const char *str, int len, int w, int k, int is_hpc, ha_mz1_v *p, const void *hf, st_mt_t *mt)
+{
+	int i, l, dbi, dbcnt = 0, kmer_span = 0;
+	tiny_queue_t tq;
+	memset(&tq, 0, sizeof(tiny_queue_t));
+	uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
+	
+    for (i = l = dbi = 0; i < len; ++i) {
+        int c = seq_nt4_table[(uint8_t)str[i]];
+        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
+
+            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;
+                y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
+                cnt = hf? ha_ft_cnt(hf, y) : 0;
+
+				for (dbi = 0; dbi < mt->n; dbi++)
+				{
+					if(p->a[dbi].x == y && p->a[dbi].rid == cnt && p->a[dbi].pos == i && p->a[dbi].rev == z && p->a[dbi].span == kmer_span)
+					{
+						if(l != (int)mt->a[dbi]) fprintf(stderr, "ERROR\n");
+						dbcnt++;
+					}
+				}
+            }
+        } else l = 0, tq.count = tq.front = 0, kmer_span = 0;
+	}
+
+	if(dbcnt != mt->n) fprintf(stderr, "ERROR\n");
+	if(mt->n != (int)p->n) fprintf(stderr, "ERROR\n");
+	for (dbi = 1; dbi < mt->n; dbi++)
+	{
+		if(p->a[dbi].pos <= p->a[dbi-1].pos || (int)mt->a[dbi] <= (int)mt->a[dbi-1])
+		{
+			fprintf(stderr, "ERROR\n");
+		}
+	}
+}
+
 /**
  * Find symmetric (w,k)-minimizers on a DNA sequence
  *
@@ -300,123 +582,140 @@ static void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
  * @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, 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;
+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, int32_t ws)
+{   ///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, 0};
+    uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
+    int i, j, l, tl = 0, buf_pos, min_pos, kmer_span = 0;
+    ha_mz1_t buf[256], min = dummy;
+    uint32_t buf_p[256], min_s = (uint32_t)-1;
+    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);
-	} 
+    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);
+    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);
+    kv_resize(uint64_t, *mt, (int64_t)p->m); mt->n = p->n;
 
-	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
+    for (i = l = tl = 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
+            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 = 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;
+            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; tl++;
+            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 = 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;
+        buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
+        buf_p[buf_pos] = l;
+        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]); kv_push(uint64_t, *mt, buf_p[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]); kv_push(uint64_t, *mt, buf_p[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 (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;
+        /**
+         * 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 (mzcmp(&min, &info) >= 0) { // a new minimum; then write the old min
+            if (l >= w + k && min.x != UINT64_MAX){
+                kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s);
+            } 
+            min = info, min_pos = buf_pos, min_s = buf_p[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); kv_push(uint64_t, *mt, min_s);
+            } 
+            ///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 = dummy; 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, min_s = buf_p[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, min_s = buf_p[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]); kv_push(uint64_t, *mt, buf_p[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]); kv_push(uint64_t, *mt, buf_p[j]);
+                    } 
+            }
+        }
+        if (++buf_pos == w) buf_pos = 0;
+    }
+    if (min.x != UINT64_MAX){
+        kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s); 
+    }
+	// debug_pl(str, len, w, k, is_hpc, p, hf, mt);
+    // if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC, dp_min_len);
+    select_mz_h(p, mt, len, sample_dist, ws, k, tl);
+    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)