Merge branch 'master' into hifiasm-v0.14

fix hc_edge
debug hic_hits
2026-09-21 04:08:11 +08:00 · 2021-02-09 18:46:29 -05:00 · 2021-02-09 18:44:06 -05:00 · 2021-02-09 18:40:49 -05:00 · 2021-02-09 13:00:23 -05:00 · 2021-02-02 12:45:50 -05:00
28 changed files with 19529 additions and 5162 deletions
@@ -11,9 +11,12 @@
 #include "htab.h"
 #include "kthread.h"

-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);
+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);
+void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);

 All_reads R_INF;
+Debug_reads R_INF_FLAG;

 void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length)
 {
@@ -296,30 +299,6 @@ void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int
    }
 }

-int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
-{
-    long long overlapLen = xEnd - xBeg + 1;
-
-    if(yEnd - yBeg + 1 == overlapLen)
-    {
-        long long i;
-
-        for (i = 0; i < overlapLen; i++)
-        {
-            if(x[xBeg + i] != y[yBeg + i])
-            {
-                break;
-            }
-        }
-
-        if(i == overlapLen)
-        {
-            return 1;
-        }
-    }
-
-    return 0;
-}

 long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag)
 {
@@ -343,6 +322,64 @@ long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_l



+            ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
+            tmp.rev = overlap_list->list[i].y_pos_strand;
+
+
+            /**********************target***************************/
+            tmp.tn = overlap_list->list[i].y_id;
+            if(tmp.rev == 1)
+            {
+                long long y_readLen = R_INF.read_length[overlap_list->list[i].y_id];
+                tmp.ts = y_readLen - overlap_list->list[i].y_pos_e - 1;
+                tmp.te = y_readLen - overlap_list->list[i].y_pos_s - 1;
+            }
+            else
+            {
+                tmp.ts = overlap_list->list[i].y_pos_s;
+                tmp.te = overlap_list->list[i].y_pos_e;
+            }
+            ///the end pos is open
+            tmp.te++;
+            /**********************target***************************/
+
+            tmp.bl = R_INF.read_length[overlap_list->list[i].y_id];
+            tmp.ml = overlap_list->list[i].strong;
+            tmp.no_l_indel = overlap_list->list[i].without_large_indel;
+
+            tmp.el = overlap_list->list[i].shared_seed;
+
+            add_ma_hit_t_alloc(paf, &tmp);
+        }
+    }
+
+    return available_overlaps;
+}
+
+
+
+long long push_final_overlaps_increment(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag)
+{
+    long long i = 0;
+    long long available_overlaps = paf->length;
+    ma_hit_t tmp;
+    ///clear_ma_hit_t_alloc(paf); // paf has been preallocated, so we don't need preallocation
+    for (i = 0; i < (long long)overlap_list->length; i++)
+    {
+        if (overlap_list->list[i].is_match == flag)
+        {
+            available_overlaps++;
+            /**********************query***************************/
+            //the interval of overlap is half-open [start, end) 
+            tmp.qns = overlap_list->list[i].x_id;
+            tmp.qns = tmp.qns << 32;
+            tmp.qns = tmp.qns | (uint64_t)(overlap_list->list[i].x_pos_s);
+            ///the end pos is open
+            tmp.qe = overlap_list->list[i].x_pos_e + 1;
+            /**********************query***************************/
+
+
+
            ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
            tmp.rev = overlap_list->list[i].y_pos_strand;

@@ -383,6 +420,7 @@ typedef struct {
 	UC_Read self_read, ovlp_read;
 	Candidates_list clist;
 	overlap_region_alloc olist;
+    overlap_region_alloc olist_hp;
 	ha_abuf_t *ab;
 	// error correction related buffers
 	int64_t num_read_base, num_correct_base, num_recorrect_base;
@@ -392,6 +430,10 @@ typedef struct {
 	Correct_dumy correct;
 	haplotype_evdience_alloc hap;
 	Round2_alignment round2;
+    kvec_t_u32_warp b_buf;
+    kvec_t_u64_warp r_buf;
+    kvec_t_u8_warp k_flag;
+    overlap_region tmp_region;
 } ha_ovec_buf_t;

 ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
@@ -403,6 +445,11 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
 	init_UC_Read(&b->ovlp_read);
 	init_Candidates_list(&b->clist);
 	init_overlap_region_alloc(&b->olist);
+    init_overlap_region_alloc(&b->olist_hp);
+    init_fake_cigar(&(b->tmp_region.f_cigar));
+    kv_init(b->b_buf.a);
+    kv_init(b->r_buf.a);
+    kv_init(b->k_flag.a);
 	b->ab = ha_abuf_init();
 	if (!b->is_final) {
 		init_Cigar_record(&b->cigar1);
@@ -421,7 +468,12 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
 	destory_UC_Read(&b->ovlp_read);
 	destory_Candidates_list(&b->clist);
 	destory_overlap_region_alloc(&b->olist);
+    destory_overlap_region_alloc(&b->olist_hp);
 	ha_abuf_destroy(b->ab);
+    destory_fake_cigar(&(b->tmp_region.f_cigar));
+    kv_destroy(b->b_buf.a);
+    kv_destroy(b->r_buf.a);
+    kv_destroy(b->k_flag.a);
 	if (!b->is_final) {
 		destory_Cigar_record(&b->cigar1);
 		destory_Graph(&b->POA_Graph);
@@ -451,6 +503,7 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
 {
 	int64_t i, mem = 0, mem_clist, mem_olist;
 	mem_clist = b->clist.size * sizeof(k_mer_hit) + b->clist.chainDP.size * 7 * 4;
+    
 	mem_olist = b->olist.size * sizeof(overlap_region);
 	for (i = 0; i < (int64_t)b->olist.size; ++i) {
 		const overlap_region *r = &b->olist.list[i];
@@ -458,6 +511,14 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
 		mem_olist += r->f_cigar.size * 8;
 		mem_olist += r->boundary_cigars.size * sizeof(window_list);
 	}
+    mem_olist += b->olist_hp.size * sizeof(overlap_region);
+	for (i = 0; i < (int64_t)b->olist_hp.size; ++i) {
+		const overlap_region *r = &b->olist_hp.list[i];
+		mem_olist += r->w_list_size * sizeof(window_list);
+		mem_olist += r->f_cigar.size * 8;
+		mem_olist += r->boundary_cigars.size * sizeof(window_list);
+	}
+
 	mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist;
 	if (!b->is_final) {
 		mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4;
@@ -475,8 +536,8 @@ static void worker_ovec(void *data, long i, int tid)
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
 	int fully_cov, abnormal;

-	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
-    ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);
+    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);

 	clear_Cigar_record(&b->cigar1);
 	clear_Round2_alignment(&b->round2);
@@ -498,37 +559,133 @@ static void worker_ovec(void *data, long i, int tid)
 	}
 	R_INF.paf[i].is_abnormal = abnormal;

-	if (b->save_ov) {
+    R_INF.trio_flag[i] = AMBIGU;
+    
+    ///need to be fixed in r305
+    // if(ha_idx_hp == NULL)
+    // {
+    //     R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), NULL);
+    // }
+
+    if (R_INF.trio_flag[i] != AMBIGU || b->save_ov) {
 		int is_rev = (asm_opt.number_of_round % 2 == 0);
 		push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, is_rev);
 		push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, is_rev);
 	}
 }

+
 static void worker_ovec_related_reads(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
-	int required_read_name_length = strlen(asm_opt.required_read_name);
-	uint64_t k;
-	if (required_read_name_length == (int)Get_NAME_LENGTH((R_INF),i)
-			&&
-			memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
-	{
-		ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
-        ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);

-		fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
-		recover_UC_Read(&b->self_read, &R_INF, i);
-		fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
+    uint64_t k, queryNameLen;
+    for (k = 0; k < R_INF_FLAG.query_num; k++)
+    {
+        queryNameLen = strlen(R_INF_FLAG.read_name[k]);
+        if (queryNameLen != Get_NAME_LENGTH((R_INF),i)) continue;
+        if (memcmp(R_INF_FLAG.read_name[k], Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
+        {
+            break;
+        }
+    }

-		for (k = 0; k < b->olist.length; k++) {
-			fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
-			recover_UC_Read(&b->self_read, &R_INF, b->olist.list[k].y_id);
-			fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
+    if(k < R_INF_FLAG.query_num)
+    {
+        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]));
+
+        clear_Cigar_record(&b->cigar1);
+        clear_Round2_alignment(&b->round2);
+
+        correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
+                &b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal);
+
+        b->num_read_base += b->self_read.length;
+        b->num_correct_base += b->correct.corrected_base;
+        b->num_recorrect_base += b->round2.dumy.corrected_base;
+
+        push_cigar(R_INF.cigars, i, &b->cigar1);
+        push_cigar(R_INF.second_round_cigar, i, &b->round2.cigar);
+
+        R_INF.paf[i].is_fully_corrected = 0;
+        if (fully_cov) {
+            if (get_cigar_errors(&b->cigar1) == 0 && get_cigar_errors(&b->round2.cigar) == 0)
+                R_INF.paf[i].is_fully_corrected = 1;
+        }
+        R_INF.paf[i].is_abnormal = abnormal;
+
+        
+
+        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, "%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);
+
+        fprintf(R_INF_FLAG.fp, "***************************k-mer counts (%d)***************************\n", (int)(R_INF_FLAG.candidate_count[q_idx].a.n));
+        
+        sort_kvec_t_u64_warp(&(R_INF_FLAG.candidate_count[q_idx]), 0);
+        for (k = 0; k < R_INF_FLAG.candidate_count[q_idx].a.n; k++) 
+        {
+            fprintf(R_INF_FLAG.fp, "[%lu] Count(%u): %lu, filtered: %lu\n", k, 
+            (uint32_t)R_INF_FLAG.candidate_count[q_idx].a.a[k], R_INF_FLAG.candidate_count[q_idx].a.a[k]>>33, 
+            (R_INF_FLAG.candidate_count[q_idx].a.a[k]>>32)&(uint64_t)1);
+        }
+
+        
+        fprintf(R_INF_FLAG.fp, "***************************forward ovlp***************************\n");
+        for (k = 0; k < b->olist.length; k++) 
+        {
+            if(b->olist.list[k].is_match != 1) continue;
+			fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
+            fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
+            b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, 
+            b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
+            Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
 		}
-	}
+
+        fprintf(R_INF_FLAG.fp, "***************************reverse ovlp***************************\n");
+        for (k = 0; k < b->olist.length; k++) 
+        {
+            if(b->olist.list[k].is_match != 2) continue;
+			fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
+            fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
+            b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, 
+            b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
+            Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
+		}
+
+        fprintf(R_INF_FLAG.fp, "***************************unmatched ovlp***************************\n");
+        for (k = 0; k < b->olist.length; k++) 
+        {
+            if(b->olist.list[k].is_match == 1) continue;
+            if(b->olist.list[k].is_match == 2) continue;
+			fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
+            fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
+            b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, 
+            b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
+            Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
+		}
+
+        R_INF.trio_flag[i] = AMBIGU;
+
+        ///need to be fixed in r305
+        // if(ha_idx_hp == NULL)
+        // {
+        //     R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), R_INF_FLAG.fp);
+        // }
+        
+        fprintf(R_INF_FLAG.fp, "R_INF.trio_flag[%ld]: %u\n", i, R_INF.trio_flag[i]);
+    
+
+	    pthread_mutex_unlock(&R_INF_FLAG.OutputMutex);
+    }
 }

+
 static inline long long get_N_occ(char* seq, long long length)
 {
 	long long j, N_occ = 0;
@@ -632,25 +789,76 @@ void Output_corrected_reads()
    fclose(output_file);
 }

+void debug_print_pob_regions()
+{
+    uint64_t i, total = 0;
+    for (i = 0; i < R_INF.total_reads; i++)
+    {
+        if(R_INF.trio_flag[i]!=AMBIGU)
+        {
+            total++;
+            fprintf(stderr, "(%lu) %.*s\n", i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
+        }
+    }
+    fprintf(stderr, "total hp reads: %lu, R_INF.total_reads: %lu\n", total, R_INF.total_reads);
+    exit(1);
+}
+
+void rescue_hp_reads(ha_ovec_buf_t **b)
+{
+    int hom_cov, het_cov;
+    ha_flt_tab_hp = ha_idx_hp = NULL;
+    if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
+       ha_flt_tab_hp = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 1);
+    }
+    ha_idx_hp = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 1, &R_INF, &hom_cov, &het_cov);
+
+
+    if (asm_opt.required_read_name)
+		kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
+	else
+		kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);
+
+
+
+
+    ha_ft_destroy(ha_flt_tab_hp); ha_flt_tab_hp = NULL;
+    ha_pt_destroy(ha_idx_hp); ha_idx_hp = NULL;
+}
+
+
+
+
 void ha_overlap_and_correct(int round)
 {
-	int i, hom_cov, het_cov;
+	int i, hom_cov, het_cov, r_out = 0;
 	ha_ovec_buf_t **b;
 	ha_ecsave_buf_t *e;
+    ha_flt_tab_hp = ha_idx_hp = NULL;

+    if((ha_idx == NULL)&&(asm_opt.flag & HA_F_VERBOSE_GFA)&&(round == asm_opt.number_of_round - 1))
+    {
+        r_out = 1;
+    }
+
+    if(asm_opt.required_read_name) init_Debug_reads(&R_INF_FLAG, asm_opt.required_read_name); // for debugging only
 	// overlap and correct reads
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1));
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov, &het_cov); // build the index
-	if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
+    if(ha_idx) hom_cov = asm_opt.hom_cov;
+	if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
+	///debug_adapter(&asm_opt, &R_INF);
+    if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
 		ha_opt_update_cov(&asm_opt, hom_cov);
 	if (asm_opt.required_read_name)
 		kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
 	else
 		kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);
+
+    if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
 	ha_pt_destroy(ha_idx);
-	ha_idx = 0;
+	ha_idx = NULL;

 	// collect statistics
 	for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -662,8 +870,8 @@ void ha_overlap_and_correct(int round)
 	}
 	free(b);

-	if (asm_opt.required_read_name) exit(0); // for debugging only
-
+	if (asm_opt.required_read_name) destory_Debug_reads(&R_INF_FLAG), exit(0); // for debugging only
+    
 	// save corrected reads to R_INF
 	CALLOC(e, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -679,8 +887,10 @@ void ha_overlap_and_correct(int round)
 		free(e[i].second_round_read);
 	}
 	free(e);
+    ///debug_print_pob_regions();
 }

+
 void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, 
 UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
 {
@@ -762,6 +972,116 @@ UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
    }
 }

+
+int check_chain_indels(Fake_Cigar* chain, long long xBeg, long long xEnd, float indel_rate)
+{
+    uint64_t i = 0;
+    long long indels = 0, xOffset;
+    if(chain->length != 0)
+    {
+        indels += abs(get_fake_gap_shift(chain, 0));
+        xOffset = get_fake_gap_pos(chain, 0);
+        if(indels > (xOffset - xBeg + 1) * indel_rate) return 0;
+
+        for (i = 1; i < chain->length; i++)
+        {
+            indels += abs((get_fake_gap_shift(chain, i) - get_fake_gap_shift(chain, i-1)));
+            xOffset = get_fake_gap_pos(chain, i);
+            if(indels > (xOffset - xBeg + 1) * indel_rate) return 0;
+        }
+    }
+
+    if(indels > (xEnd - xBeg + 1) * indel_rate) return 0;
+    return 1;
+}
+
+void update_overlaps_chain_width(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, 
+UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact, float indel_rate)
+{
+
+    uint64_t inner_j = 0;
+    uint64_t j = 0;
+    long long x_overlapLen, y_overlapLen;
+    while (j < overlap_list->length && inner_j < paf->length)
+    {
+        if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn)
+        {
+            j++;
+        }
+        else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn)
+        {
+            inner_j++;
+        }
+        else
+        {
+            if(check_chain_indels(&(overlap_list->list[j].f_cigar), overlap_list->list[j].x_pos_s, 
+            overlap_list->list[j].x_pos_e, indel_rate) == 1)
+            {
+                if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev)
+                {
+                    x_overlapLen = Get_qe(paf->buffer[inner_j]) - Get_qs(paf->buffer[inner_j]) + 1;
+                    y_overlapLen = Get_te(paf->buffer[inner_j]) - Get_ts(paf->buffer[inner_j]) + 1;
+                    if(x_overlapLen < y_overlapLen) x_overlapLen = y_overlapLen;
+                    x_overlapLen = x_overlapLen * 0.1;
+
+                    // if(
+                    // ((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < x_overlapLen)
+                    // && (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < x_overlapLen))
+                    // ||
+                    // ((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < x_overlapLen)
+                    // && (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < x_overlapLen)))
+                    if(
+                    ((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < (uint64_t)x_overlapLen)
+                    && (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < (uint64_t)x_overlapLen))
+                    || 
+                    ((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < (uint64_t)x_overlapLen)
+                    && (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < (uint64_t)x_overlapLen))
+                    )
+                    {
+                        overlap_list->list[j].is_match = is_match;
+                        overlap_list->list[j].strong = paf->buffer[inner_j].ml;
+                        overlap_list->list[j].without_large_indel = paf->buffer[inner_j].no_l_indel;
+                        if(is_exact == 1)
+                        {
+                            if(overlap_list->list[j].y_pos_strand == 0)
+                            {
+                                recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id);
+                            }
+                            else
+                            {
+                                recover_UC_Read_RC(overlap_read, &R_INF, overlap_list->list[j].y_id);
+                            }
+                            if(if_exact_match(g_read->seq, g_read->length, overlap_read->seq, overlap_read->length, 
+                            overlap_list->list[j].x_pos_s, overlap_list->list[j].x_pos_e, 
+                            overlap_list->list[j].y_pos_s, overlap_list->list[j].y_pos_e))
+                            {
+                                overlap_list->list[j].shared_seed = 1;
+                            }
+                            else
+                            {
+                                overlap_list->list[j].shared_seed = 0;
+                            }
+                        }
+                    }
+                    else
+                    {
+                        overlap_list->list[j].is_match = 3;
+                    }
+                }
+                else
+                {
+                    overlap_list->list[j].is_match = 3;
+                }
+            }
+
+            j++;
+            inner_j++;
+        }
+    }
+}
+
+
+
 void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read, UC_Read* overlap_read)
 {
    uint64_t j;
@@ -769,6 +1089,12 @@ void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read,
    {
        if (overlap_list->list[j].is_match != 1)
        {
+            if((overlap_list->list[j].x_pos_e + 1 - overlap_list->list[j].x_pos_s) != 
+               (overlap_list->list[j].y_pos_e + 1 - overlap_list->list[j].y_pos_s))
+            {
+                continue;
+            }
+
            if(overlap_list->list[j].y_pos_strand == 0)
            {
                recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id);
@@ -819,9 +1145,11 @@ void ha_print_ovlp_stat(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long
 void fill_chain(Fake_Cigar* chain, char* x_string, char* y_string, long long xBeg, long long yBeg,
 long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
 {
-    long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, mapScore, zdroped;
+    /**
+    long long i, xOffset, yOffset, xRegionLen, yRegionLen, maxXpos, maxYpos, mapGlobalScore, mapExtentScore, zdroped;
+    long long xBuoundaryScore, yBuoundaryScore;
    ///float band_rate = 0.08;
-    int endbouns;
+    int endbouns,mode;
    if(chain->length <= 0) return;

    kvec_t(uint8_t) x_num;
@@ -850,7 +1178,9 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
        ///text is x, query is y
        afine_gap_alignment(x_string, x_num.a, xRegionLen, y_string, y_num.a, yRegionLen, 
        c2n, BACKWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
-        /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapScore, &zdroped);
+        BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, 
+        &mapExtentScore, &xBuoundaryScore, &yBuoundaryScore, &zdroped);
+
        // fprintf(stderr, "* xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n",
        // xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped);
    }
@@ -891,7 +1221,8 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
        ///text is x, query is y
        afine_gap_alignment(x_string+xOffset, x_num.a, xRegionLen, y_string+yOffset, y_num.a, yRegionLen, 
        c2n, FORWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
-        /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapScore, &zdroped);
+        BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, 
+        &mapExtentScore, &xBuoundaryScore, &yBuoundaryScore, &zdroped);
        // fprintf(stderr, "# xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n",
        // xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped);
    }
@@ -899,6 +1230,7 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)

    kv_destroy(x_num);
    kv_destroy(y_num);
+    **/
 }
 void Final_phasing(overlap_region_alloc* overlap_list, Cigar_record_alloc* cigarline,
 UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n)
@@ -952,21 +1284,10 @@ UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n)
 static void worker_ov_final(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
-    uint8_t c2n[256]; // this may be moved to ha_ovec_buf_t, but it should be fast to populate anyway
-    memset(c2n, 4, 256);
-    c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1;
-	c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table

 	//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
-	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0);
-    ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 0);
-
-	/**
-	  correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
-	  &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
-	  &current_cigar, &hap, &second_round, 0, 0);
-	  push_final_overlaps(&(R_INF.paf[i]), &overlap_list);
-	 **/
+    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);

 	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);
@@ -982,6 +1303,81 @@ static void worker_ov_final(void *data, long i, int tid)
 	push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
 }

+
+
+
+void reset_final_overlaps(overlap_region_alloc *overlap_list)
+{
+    
+    uint64_t i;
+    for (i = 0; i < overlap_list->length; i++)
+    {
+        if (overlap_list->list[i].is_match == 1 || overlap_list->list[i].is_match == 2)
+        {
+            overlap_list->list[i].x_pos_s = overlap_list->list[i].x_pos_e = (uint32_t)-1;
+            overlap_list->list[i].y_pos_s = overlap_list->list[i].y_pos_e = (uint32_t)-1;
+            overlap_list->list[i].is_match = 0;
+        }
+    }
+
+    ha_sort_list_by_anchor(overlap_list);
+}
+
+void debug_affine_gap_alignment(overlap_region_alloc *overlap_list, UC_Read* g_read, UC_Read* overlap_read)
+{
+    uint64_t i;
+    kvec_t(uint8_t) x_num;
+    kvec_t(uint8_t) y_num;
+    kv_init(x_num);
+    kv_init(y_num);
+    for (i = 0; i < overlap_list->length; i++)
+    {
+        if (overlap_list->list[i].is_match == 1 && overlap_list->list[i].shared_seed == 1)
+        {
+
+            kv_resize(uint8_t, x_num, (uint64_t)(Get_READ_LENGTH(R_INF, overlap_list->list[i].x_id)));
+            kv_resize(uint8_t, y_num, (uint64_t)(Get_READ_LENGTH(R_INF, overlap_list->list[i].y_id)));
+
+            get_affine_gap_score(&(overlap_list->list[i]), g_read, overlap_read, x_num.a, y_num.a,
+            overlap_list->list[i].x_pos_e + 1 - overlap_list->list[i].x_pos_s, 
+            overlap_list->list[i].y_pos_e + 1 - overlap_list->list[i].y_pos_s);
+        }
+    }
+
+    kv_destroy(x_num);
+    kv_destroy(y_num);
+}
+
+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);
+
+    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);
+    ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length);
+
+
+    ///update_overlaps(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1);
+    update_overlaps_chain_width(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1, 0.002);
+    update_overlaps(&b->olist, &(R_INF.reverse_paf[i]), &b->self_read, &b->ovlp_read, 2, 0);
+    ///recover missing exact overlaps
+    update_exact_overlaps(&b->olist, &b->self_read, &b->ovlp_read);
+
+    
+    ///Final_phasing(&overlap_list, &cigarline, &g_read, &overlap_read, c2n);
+    push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1);
+    push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
+
+    ///debug_affine_gap_alignment(&b->olist, &b->self_read, &b->ovlp_read);
+
+    reset_final_overlaps(&b->olist);
+    correct_overlap_high_het(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read);
+    push_final_overlaps_increment(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
+}
+
 void Output_PAF()
 {
    fprintf(stderr, "Writing PAF to disk ...... \n");
@@ -1156,7 +1552,7 @@ void hap_recalculate_peaks(char* output_file_name)
    int hom_cov, het_cov;
    // construct hash table for high occurrence k-mers
    if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
-        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
+        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
        ha_opt_update_cov(&asm_opt, hom_cov);
    }
    free(R_INF.read_length);
@@ -1164,7 +1560,7 @@ void hap_recalculate_peaks(char* output_file_name)

    load_All_reads(&R_INF, gfa_name);

-    ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
+    ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
    asm_opt.hom_cov = hom_cov;
    asm_opt.het_cov = het_cov;
    ha_pt_destroy(ha_idx);
@@ -1181,11 +1577,21 @@ void ha_overlap_final(void)
 {
 	int i, hom_cov, het_cov;
 	ha_ovec_buf_t **b;
+    ha_flt_tab_hp = ha_idx_hp = NULL;
+
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
-		b[i] = ha_ovec_init(1, 1);
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
-	kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
+		b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1);
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
+    if(asm_opt.flag & HA_F_HIGH_HET)
+    {
+        kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads);
+    }
+    else
+    {
+        kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
+    }
+    
 	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
 	for (i = 0; i < asm_opt.thread_num; ++i)
@@ -1195,6 +1601,8 @@ void ha_overlap_final(void)
    asm_opt.het_cov = het_cov;
 }

+
+
 int ha_assemble(void)
 {
 	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
@@ -1213,14 +1621,18 @@ int ha_assemble(void)
        if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
 	}
 	if (!ovlp_loaded) {
+        ha_flt_tab = ha_idx = NULL;
+        if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
+
 		// construct hash table for high occurrence k-mers
-		if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
-			ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
+		if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) 
+        {
+			ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
 			ha_opt_update_cov(&asm_opt, hom_cov);
 		}
 		// error correction
 		assert(asm_opt.number_of_round > 0);
-		for (r = 0; r < asm_opt.number_of_round; ++r) {
+		for (r = ha_idx?asm_opt.number_of_round-1:0; r < asm_opt.number_of_round; ++r) {
 			ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
 			ha_overlap_and_correct(r);
 			fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
@@ -8,6 +8,11 @@
 #define Get_Cigar_Type(RECORD) (RECORD&3)
 #define Get_Cigar_Length(RECORD) (RECORD>>2)

+#define RESEED_DP 4
+#define RESEED_PEAK_RATE 0.15
+#define RESEED_LEN 2000
+#define RESEED_HP_RATE 0.9
+
 int ha_assemble(void);

 #endif
@@ -23,6 +23,13 @@ static ko_longopt_t long_options[] = {
 	{ "ex-iter",       ko_required_argument, 308 },
    { "purge-cov",     ko_required_argument, 309 },
    { "pri-range",     ko_required_argument, 310 },
+    { "high-het",      ko_no_argument, 311 },
+    { "lowQ",          ko_required_argument, 312 },
+	{ "min-hist-cnt",  ko_required_argument, 313 },
+    { "h1",      ko_required_argument, 314 },
+    { "h2",      ko_required_argument, 315 },
+    { "enzyme",      ko_required_argument, 316 },
+    { "b-cov",      ko_required_argument, 317 },
 	{ 0, 0, 0 }
 };

@@ -58,6 +65,11 @@ void Print_H(hifiasm_opt_t* asm_opt)
    fprintf(stderr, "    -x FLOAT    max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
    fprintf(stderr, "    -y FLOAT    min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
    fprintf(stderr, "    -u          disable post join contigs step which may improve N50\n");
+    fprintf(stderr, "    --lowQ      INT\n");
+    fprintf(stderr, "                output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
+    fprintf(stderr, "    --b-cov     INT\n");
+    fprintf(stderr, "                break contigs at breakpoints with coverage drop at <INT-fold coverage [%d]\n", asm_opt->break_cov);
+    
 //	fprintf(stderr, "    --pri-range INT1[,INT2]\n");
 //	fprintf(stderr, "                keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");

@@ -77,6 +89,11 @@ void Print_H(hifiasm_opt_t* asm_opt)
                                     asm_opt->purge_overlap_len);
    fprintf(stderr, "    --purge-cov INT\n");
    fprintf(stderr, "                coverage upper bound of Purge-dups [auto]\n");
+    fprintf(stderr, "    --high-het  enable this mode for high heterozygosity sample [experimental, not stable]\n");
+    
+    fprintf(stderr, "  Hi-C-partition [experimental, not stable]:\n");
+    fprintf(stderr, "    --h1 FILEs   file names of Hi-C R1  [r1_1.fq,r1_2.fq,...]\n");
+    fprintf(stderr, "    --h2 FILEs   file names of Hi-C R2  [r2_1.fq,r2_2.fq,...]\n");

    fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
    fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -91,8 +108,12 @@ void init_opt(hifiasm_opt_t* asm_opt)
    asm_opt->read_file_names = NULL;
    asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);
    asm_opt->required_read_name = NULL;
+    asm_opt->hic_enzymes = NULL;
+    asm_opt->hic_reads[0] = NULL;
+    asm_opt->hic_reads[1] = NULL;
    asm_opt->thread_num = 1;
    asm_opt->k_mer_length = 51;
+    asm_opt->hic_mer_length = 31;
 	asm_opt->mz_win = 51;
 	asm_opt->bf_shift = 37;
 	asm_opt->high_factor = 5.0;
@@ -101,8 +122,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
 	asm_opt->hom_cov = 20;
    asm_opt->het_cov = -1024;
 	asm_opt->max_n_chain = 100;
-    asm_opt->k_mer_min_freq = 3;
-    asm_opt->k_mer_max_freq = 66;
+	asm_opt->min_hist_kmer_cnt = 5;
    asm_opt->load_index_from_disk = 1;
    asm_opt->write_index_to_disk = 1;
    asm_opt->number_of_round = 3;
@@ -123,18 +143,40 @@ void init_opt(hifiasm_opt_t* asm_opt)
    asm_opt->purge_level_primary = 2;
    asm_opt->purge_level_trio = 0;
    asm_opt->purge_simi_rate = 0.75;
+    asm_opt->purge_simi_rate_hic = 0.85;
    asm_opt->purge_overlap_len = 1;
+    asm_opt->purge_overlap_len_hic = 50;
    asm_opt->recover_atg_cov_min = -1024;
    asm_opt->recover_atg_cov_max = INT_MAX;
    asm_opt->hom_global_coverage = -1;
+    asm_opt->bed_inconsist_rate = 70;
+    asm_opt->hic_inconsist_rate = 30;
+    ///asm_opt->bub_mer_length = 3;
+    asm_opt->bub_mer_length = 1000000;
+    asm_opt->break_cov = 0;
+}
+
+void destory_enzyme(enzyme* f)
+{
+    int i;
+    if(f != NULL)
+    {
+        for (i = 0; i < f->n; i++)
+        {
+            free(f->a[i]);
+        }
+        free(f->a);
+        free(f->l);
+        free(f);
+    }
 }

 void destory_opt(hifiasm_opt_t* asm_opt)
 {
-    if(asm_opt->read_file_names != NULL)
-    {
-        free(asm_opt->read_file_names);
-    }
+    if(asm_opt->read_file_names != NULL) free(asm_opt->read_file_names);
+    if(asm_opt->hic_enzymes != NULL) destory_enzyme(asm_opt->hic_enzymes);
+    if(asm_opt->hic_reads[0] != NULL) destory_enzyme(asm_opt->hic_reads[0]);
+    if(asm_opt->hic_reads[1] != NULL) destory_enzyme(asm_opt->hic_reads[1]);
 }

 void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
@@ -174,6 +216,16 @@ static int check_file(char* name, const char* opt)
    return 1;
 }

+static int check_hic_reads(enzyme* f, const char* opt)
+{
+    int i;
+    for (i = 0; i < f->n; i++)
+    {
+        if(check_file(f->a[i], opt) == 0) return 0;
+    }
+    return 1;
+}
+
 int check_option(hifiasm_opt_t* asm_opt)
 {
    if(asm_opt->read_file_names == NULL || asm_opt->num_reads == 0)
@@ -316,11 +368,56 @@ int check_option(hifiasm_opt_t* asm_opt)
        return 0;
    }

+    if(asm_opt->bed_inconsist_rate < 0 || asm_opt->bed_inconsist_rate > 100)
+    {
+        fprintf(stderr, "[ERROR] inconsistency rate should be [0, 100] (--lowQ)\n");
+        return 0;
+    }
+

    if(asm_opt->fn_bin_yak[0] != NULL && check_file(asm_opt->fn_bin_yak[0], "YAK1") == 0) return 0;
    if(asm_opt->fn_bin_yak[1] != NULL && check_file(asm_opt->fn_bin_yak[1], "YAK2") == 0) return 0;
    if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0;
    if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0;
+    if(asm_opt->required_read_name != NULL && check_file(asm_opt->required_read_name, "b") == 0) return 0;
+    
+    if(asm_opt->hic_reads[0] != NULL && check_hic_reads(asm_opt->hic_reads[0], "HIC1") == 0) return 0;
+    if(asm_opt->hic_reads[1] != NULL && check_hic_reads(asm_opt->hic_reads[1], "HIC2") == 0) return 0;
+    if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] == NULL)
+    {
+        fprintf(stderr, "[ERROR] lack r2 of HiC reads (--h2)\n");
+        return 0;
+    }
+    if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[0] == NULL)
+    {
+        fprintf(stderr, "[ERROR] lack r1 of HiC reads (--h1)\n");
+        return 0;
+    }
+
+    if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] != NULL && 
+                            asm_opt->hic_reads[0]->n != asm_opt->hic_reads[1]->n)
+    {
+        fprintf(stderr, "[ERROR] wrong r1 and r2 of HiC reads (--h1 && --h2)\n");
+        return 0;
+    }
+
+    if(asm_opt->hic_enzymes != NULL && asm_opt->hic_enzymes->n == 0)
+    {
+        fprintf(stderr, "[ERROR] wrong HiC enzymes (--enzyme)\n");
+        return 0;
+    }
+
+    if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[0]->n == 0)
+    {
+        fprintf(stderr, "[ERROR] wrong r1 of HiC reads (--h1)\n");
+        return 0;
+    }
+
+    if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[1]->n == 0)
+    {
+        fprintf(stderr, "[ERROR] wrong r2 of HiC reads (--h2)\n");
+        return 0;
+    }

    // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
    // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
@@ -371,6 +468,61 @@ void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt)
    }
 }

+void get_hic_enzymes(char *argv, enzyme** x, int check_name)
+{
+    int i, k, pre_i, len = strlen(argv);
+    (*x) = (enzyme*)calloc(1, sizeof(enzyme));
+    if(len == 0)
+    {
+        (*x)->n = 0; (*x)->l = NULL; (*x)->a = NULL;
+        return;
+    }
+
+
+    (*x)->n = 1;
+    for (i = pre_i = 0; i < len; i++)
+    {
+        if(argv[i] == ',')
+        {
+            (*x)->n++;
+            continue;
+        } 
+        
+        if(check_name)
+        {
+            if(argv[i] != 'A' && argv[i] != 'C' && argv[i] != 'G' && argv[i] != 'T' &&
+                argv[i] != 'a' && argv[i] != 'c' && argv[i] != 'g' && argv[i] != 't' && 
+                argv[i] != 'N' && argv[i] != 'n')
+            {
+                (*x)->n = 0;
+                (*x)->l = NULL;
+                (*x)->a = NULL;
+                return;
+            }
+        }
+        
+    }
+    (*x)->l = (int*)calloc((*x)->n, sizeof(int));
+    (*x)->a = (char**)calloc((*x)->n, sizeof(char*));
+
+    for (i = pre_i = k = 0; i < len; i++)
+    {
+        if(argv[i] == ',')
+        {
+            (*x)->l[k] = i - pre_i;
+            (*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
+            memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
+            (*x)->a[k][(*x)->l[k]] = '\0';
+            pre_i = i + 1;
+            k++;
+        }
+    }
+
+    (*x)->l[k] = i - pre_i;
+    (*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
+    memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
+    (*x)->a[k][(*x)->l[k]] = '\0';
+}

 int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 {
@@ -433,7 +585,14 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
            {
                asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1;
            }
-        } 
+        }
+        else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
+        else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg);
+        else if (c == 313) asm_opt->min_hist_kmer_cnt = atoi(opt.arg);
+        else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
+        else if (c == 315) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[1]), 0);
+        else if (c == 316) get_hic_enzymes(opt.arg, &(asm_opt->hic_enzymes), 1);
+        else if (c == 317) asm_opt->break_cov = atoi(opt.arg);
        else if (c == 'l')
        {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
            asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
@@ -452,15 +611,16 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		}
    }

-
+    
    if (argc == opt.ind)
    {
        Print_H(asm_opt);
        return 0;
    }
-    ///fprintf(stderr, "max_ov_diff_ec: %f, max_ov_diff_final: %f\n", asm_opt->max_ov_diff_ec, asm_opt->max_ov_diff_final);

    get_queries(argc, argv, &opt, asm_opt);

+
+
    return check_option(asm_opt);
 }
@@ -3,7 +3,7 @@

 #include <pthread.h>

-#define HA_VERSION "0.9-r289"
+#define HA_VERSION "0.14-r309"

 #define VERBOSE 0

@@ -17,9 +17,15 @@
 #define HA_F_PURGE_JOIN      0x80
 #define HA_F_BAN_POST_JOIN   0x100
 #define HA_F_BAN_ASSEMBLY    0x200
+#define HA_F_HIGH_HET        0x400

 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap

+typedef struct{
+    int *l, n; 
+    char **a;
+}enzyme;
+
 typedef struct {
 	int flag;
    int num_reads;
@@ -29,9 +35,13 @@ typedef struct {
 	char *fn_bin_yak[2];
 	char *fn_bin_list[2];
 	char *extract_list;
+    enzyme *hic_reads[2];
+    enzyme *hic_enzymes;
 	int extract_iter;
    int thread_num;
    int k_mer_length;
+    int hic_mer_length;
+    int bub_mer_length;
 	int mz_win;
 	int bf_shift;
 	double high_factor; // coverage cutoff set to high_factor*hom_cov
@@ -39,9 +49,9 @@ typedef struct {
 	double max_ov_diff_final;
 	int hom_cov;
    int het_cov;
+    int break_cov;
 	int max_n_chain; // fall-back max number of chains to consider
-    int k_mer_min_freq;
-    int k_mer_max_freq;
+	int min_hist_kmer_cnt;
    int load_index_from_disk;
    int write_index_to_disk;
    int number_of_round;
@@ -58,14 +68,18 @@ typedef struct {
    int purge_level_primary;
    int purge_level_trio;
    int purge_overlap_len;
+    int purge_overlap_len_hic;
    int recover_atg_cov_min;
    int recover_atg_cov_max;
    int hom_global_coverage;
+    int bed_inconsist_rate;
+    int hic_inconsist_rate;

    float max_hang_rate;
    float min_drop_rate;
    float max_drop_rate;
    float purge_simi_rate;
+    float purge_simi_rate_hic;

    long long small_pop_bubble_size;
    long long large_pop_bubble_size;
@@ -91,4 +105,9 @@ static inline int ha_opt_triobin(const hifiasm_opt_t *opt)
 	return ((opt->fn_bin_yak[0] && opt->fn_bin_yak[1]) || (opt->fn_bin_list[0] && opt->fn_bin_list[1]));
 }

+static inline int ha_opt_hic(const hifiasm_opt_t *opt)
+{
+    return ((opt->hic_reads[0] && opt->hic_reads[1]));
+}
+
 #endif
@@ -21,8 +21,8 @@

 ///#define FLAG_THRE 0

-#define MAX(x, y) ((x >= y)?(x):(y))
-#define MIN(x, y) ((x <= y)?(x):(y))
+#define MAX(x, y) (((x) >= (y))?(x):(y))
+#define MIN(x, y) (((x) <= (y))?(x):(y))
 #define DIFF(x, y) ((MAX((x), (y))) - (MIN((x), (y))))
 #define OVERLAP(x_start, x_end, y_start, y_end) (MIN(x_end, y_end) - MAX(x_start, y_start) + 1) 
 ///#define OVERLAP(x_start, x_end, y_start, y_end) MIN(x_end, y_end) - MAX(x_start, y_start) + 1 
@@ -1162,10 +1162,32 @@ void init_Cigar_record_alloc(Cigar_record_alloc* x);
 void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size);
 void destory_Cigar_record_alloc(Cigar_record_alloc* x);

-void afine_gap_alignment(const char *tseq, uint8_t* tnum, const int tl, 
-const char *qseq, uint8_t* qnum, const int ql, const uint8_t *c2n, const int strand,
+void afine_gap_alignment(const char *qseq, uint8_t* qnum, const int ql, 
+const char *tseq, uint8_t* tnum, const int tl, const uint8_t *c2n, const int strand,
 int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop, int end_bonus,
-long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped);
+long long* max_q_pos, long long* max_t_pos, long long* global_score, 
+long long* extention_score, long long* q_boundary_score, long long* q_boundary_t_coordinate,
+long long* t_boundary_score, long long* t_boundary_q_coordinate,
+long long* droped, int mode);
+void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
+long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
+uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen);
+int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, float hp_rate, int rlen, FILE* fp);
+
+inline int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
+{
+    long long overlapLen = xEnd - xBeg + 1;
+
+    if(yEnd - yBeg + 1 == overlapLen)
+    {
+        if(memcmp(x+xBeg, y+yBeg, overlapLen)==0) return 1;
+    }
+
+    return 0;
+}
+
+

 #define FORWARD_KSW 0
 #define BACKWARD_KSW 1
@@ -1174,5 +1196,5 @@ long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped)
 #define GAP_OPEN_KSW 4
 #define GAP_EXT_KSW 2
 #define Z_DROP_KSW 400
-#define BAND_KSW 50
+#define BAND_KSW 500
 #endif
@@ -305,6 +305,9 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
    }
 }

+
+
+
 long long get_chainLen(long long x_beg, long long x_end, long long xLen, 
                       long long y_beg, long long y_end, long long yLen)
 {
@@ -337,6 +340,50 @@ long long get_chainLen(long long x_beg, long long x_end, long long xLen,
    return x_end - x_beg + 1;
 }

+
+void debug_chain_single_site(k_mer_hit* a, long long a_n, Chain_Data* dp, int x_readLen, int y_readLen, int s_index)
+{
+    long long j, current_j = s_index;
+    long long selfLen = 0, indels = 0;
+    long long distance_self_pos, distance_pos, distance_gap;
+
+    j = s_index;
+    while (j >= 0)
+    {
+        current_j = j;
+
+        j = dp->pre[j];
+
+        if(j != -1)
+        {
+            distance_self_pos = a[current_j].self_offset - a[j].self_offset;
+            distance_pos = a[current_j].offset - a[j].offset;
+            distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos;
+
+            indels += distance_gap; 
+            selfLen += distance_self_pos;
+        }
+        fprintf(stderr, "j: %lld, score: %lld, occ: %d, pre_j: %lld\n", 
+                    current_j, (long long)dp->score[current_j], dp->occ[current_j], j);
+    }
+
+    fprintf(stderr, "s_self_offset: %u, s_offset: %u, e_self_offset: %u, e_offset: %u, ovlp length: %lld, x_readLen: %d, y_readLen: %d\n", 
+    a[s_index].self_offset, a[s_index].offset, a[current_j].self_offset, a[current_j].offset,
+    get_chainLen(a[s_index].self_offset, a[current_j].self_offset, x_readLen, 
+                                     a[s_index].offset, a[current_j].offset, y_readLen), x_readLen, y_readLen);
+
+    if(indels != dp->indels[s_index])
+    {
+        fprintf(stderr, "indels: %lld, dp->indels[i]: %ld\n", indels, (long)dp->indels[s_index]);
+    }
+
+    if(selfLen != dp->self_length[s_index])
+    {
+        fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %ld\n", selfLen, (long)dp->self_length[s_index]);
+    }
+    fprintf(stderr,"\n");
+}
+
 int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc, double bw_thres)
 {
 	int32_t i, tot_indel = 0, tot_len = 0;
@@ -348,7 +395,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
 	if (i < n_a) return -1;
 	bw_pen = 1.0 / bw_thres;
 	dp->score[0] = a[0].good? min_sc : min_sc>>1;
-	dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0;
+	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;
 		int32_t dx = (int32_t)a[i].offset - (int32_t)a[i-1].offset;
@@ -368,6 +415,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
 		dp->pre[i] = i - 1;
 		dp->indels[i] = tot_indel;
 		dp->self_length[i] = tot_len;
+        dp->occ[i] = i + 1;
 	}
 	if (i < n_a) return -1;
 	return n_a;
@@ -391,6 +439,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
    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;
@@ -435,7 +484,9 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
            ///min distance
            distance_min = distance_pos < distance_self_pos? distance_pos:distance_self_pos;
            score = distance_min < min_score? distance_min : min_score;
-			if (!a[j].good) score >>= 1;
+            ///need to be fixed in r305
+			///if (!a[j].good) score = (score >> 1) + (score & 1);
+            if (!a[j].good) score >>= 1;

            gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
            ///if the gap rate > 0.06, score will be negative
@@ -444,7 +495,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
            score += dp->score[j];

 			///find a new max score
-			if (score > max_score) {
+			if (score > max_score) {///must use > instead of >=
 				max_score = score;
 				max_j = j;
 				max_indels = total_indels;
@@ -466,10 +517,12 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
        dp->pre[i] = max_j;
        dp->indels[i] = max_indels;
        dp->self_length[i] = max_self_length;
+        dp->occ[i] = 1;
+        if(max_j != -1) dp->occ[i] = dp->occ[max_j] + 1;
    }

    ///debug_chain(a, a_n, dp);
-
+    
 skip_dp:

    max_score = -1;
@@ -562,7 +615,7 @@ skip_dp:
    }
 }

-void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
+void calculate_overlap_region_by_chaining_back(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
                                          uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end)
 {
    overlap_region tmp_region;
@@ -627,6 +680,71 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r
    destory_fake_cigar(&(tmp_region.f_cigar));
 }

+
+void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,
+                                          uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar)
+{
+    long long i = 0;
+    uint64_t current_ID;
+    uint64_t current_stand;
+
+    if (candidates->length == 0)
+    {
+        return;
+    }
+
+    long long sub_region_beg;
+    long long sub_region_end;
+
+    clear_fake_cigar(&((*f_cigar).f_cigar));
+
+    i = 0;
+    while (i < candidates->length)
+    {
+        chain_idx->a.n = 0;
+        current_ID = candidates->list[i].readID;
+        current_stand = candidates->list[i].strand;
+
+        ///reference read
+        (*f_cigar).x_id = readID;
+        (*f_cigar).x_pos_strand = current_stand;
+        ///query read
+        (*f_cigar).y_id = current_ID;
+        ///here the strand of query is always 0
+        (*f_cigar).y_pos_strand = 0;  
+
+        sub_region_beg = i;
+        sub_region_end = i;
+        i++;
+
+        while (i < candidates->length 
+        && 
+        current_ID == candidates->list[i].readID
+        &&
+        current_stand == candidates->list[i].strand)
+        {
+            sub_region_end = i;
+            i++;
+        }
+
+        if ((*f_cigar).x_id == (*f_cigar).y_id)
+        {
+            continue;
+        }
+
+		chain_DP(candidates->list + sub_region_beg,
+				sub_region_end - sub_region_beg + 1, &(candidates->chainDP), f_cigar, band_width_threshold,
+				25, Get_READ_LENGTH((*R_INF), (*f_cigar).x_id), Get_READ_LENGTH((*R_INF), (*f_cigar).y_id));
+
+        ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1)
+        if ((*f_cigar).x_id != (*f_cigar).y_id)
+        {
+            append_inexact_overlap_region_alloc(overlap_list, f_cigar, R_INF, add_beg_end);
+        }
+    }
+}
+
+
 void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
                        int extra_begin, int extra_end, int error_threshold)
 {
@@ -704,6 +822,7 @@ void destory_Chain_Data(Chain_Data* x)
    free(x->pre);
    free(x->indels);
    free(x->self_length);
+    free(x->occ);
 	free(x->tmp);
 }

@@ -716,6 +835,7 @@ void resize_Chain_Data(Chain_Data* x, long long size)
 		REALLOC(x->pre, x->size);
 		REALLOC(x->indels, x->size);
 		REALLOC(x->self_length, x->size);
+        REALLOC(x->occ, x->size);
 		REALLOC(x->tmp, x->size);
 	}
 }
@@ -12,6 +12,8 @@
 #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
 #define THRESHOLD  15
 #define OVERLAP_THRESHOLD_FILTER 0.9
+#define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3
+#define HIGH_HET_ERROR_RATE 0.08
 #define THRESHOLD_MAX_SIZE  31

 #define GROUP_SIZE 4
@@ -117,6 +119,7 @@ typedef struct {
 	int64_t *pre;
 	int32_t *indels;
 	int32_t *self_length;
+    int32_t *occ;
 	int64_t *tmp; // MUST BE 64-bit integer
 	int64_t length;
 	int64_t size;
@@ -142,8 +145,8 @@ int extra_begin, int extra_end, int error_threshold);

 void overlap_region_sort_y_id(overlap_region *a, long long n);

-void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
-uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end);
+void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,
+uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar);

 void init_fake_cigar(Fake_Cigar* x);
 void destory_fake_cigar(Fake_Cigar* x);
@@ -1,10 +1,12 @@
 CXX=		g++
+CC=			gcc
 CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
+CFLAGS=		$(CXXFLAGS)
 CPPFLAGS=
 INCLUDES=
 OBJS=		CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
 			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
-			htab.o hist.o sketch.o anchor.o extract.o sys.o
+			htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm

@@ -13,12 +15,15 @@ ifneq ($(asan),)
 	LIBS+=-fsanitize=address
 endif

-.SUFFIXES:.cpp .o
+.SUFFIXES:.cpp .c .o
 .PHONY:all clean depend

 .cpp.o:
 		$(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@

+.c.o:
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
+
 all:$(EXE)

 $(EXE):$(OBJS) main.o
@@ -38,6 +43,7 @@ Assembly.o: kthread.h
 CommandLines.o: CommandLines.h ketopt.h
 Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
 Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
+Correct.o: ksw2.h
 Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
 Hash_Table.o: CommandLines.h ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h
@@ -65,3 +71,4 @@ main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h
 main.o: Levenshtein_distance.h htab.h
 sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
 sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
+hic.o: hic.h
@@ -21,8 +21,11 @@
 #define PRIMARY_LABLE 0
 #define ALTER_LABLE 1
 #define HAP_LABLE 2
+#define FAKE_LABLE 4
 #define TRIO_THRES 0.9
-#define DOUBLE_CHECK_THRES 0.2
+#define DOUBLE_CHECK_THRES 0.1
+#define FINAL_DOUBLE_CHECK_THRES 0.2
+#define CHIMERIC_TRIM_THRES 4
 // #define PRIMARY_LABLE 1
 // #define ALTER_LABLE 2
 // #define HAP_LABLE 4
@@ -104,6 +107,11 @@ typedef struct {
 	uint8_t no_l_indel;
 } asg_arc_t;

+typedef struct {
+	size_t n, m;
+	asg_arc_t* a;
+} kv_asg_arc_t;
+

 typedef struct {
 	uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
@@ -153,6 +161,7 @@ typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;
 typedef struct {
 	ma_utg_v u;
 	asg_t *g;
+	kvec_t(uint64_t) occ;
 } ma_ug_t;

 typedef struct {
@@ -392,7 +401,8 @@ typedef struct {
 }kvec_asg_arc_t_warp;

 void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend);
-
+int asg_arc_del_multi(asg_t *g);
+int asg_arc_del_asymm(asg_t *g);

 typedef struct {
 	uint32_t q_pos;
@@ -440,11 +450,8 @@ long long max_hang_length, long long clean_round, long long gap_fuzz,
 float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
 long long bubble_dist, int read_graph, int write);

-void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, 
-ma_hit_t_alloc* reverse_sources, int id, char* command);
-
+void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, char* command);
 void collect_abnormal_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum);
-
 void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen);
 void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
 void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 
@@ -453,18 +460,6 @@ uint64_t* source_index, long long listLen);
 #define EvaluateLen(U, id) ((U).a[(id)].start)
 #define IsMerge(U, id) ((U).a[(id)].end)
 #define kv_reuse(v, rn, rm, r) ((v).n = (rn), (v).m = (rm), (v).a = (r))
-#define long_tip(U, id, threshold) ((EvaluateLen((U), (id))>=(threshold))&&(!((U).a[(id)].circ)))
-///there are threee cases: 
-///1. if this untig is too long (>maxShortUntig), it must be not short untig/must be a long untig
-///2. if this untig is long (>minLongUntig && EvaluateLen(ug->u, av[i].v>>1) > (EvaluateLen(ug->u, v>>1)*l_untig_rate)), it might be a long tip
-#define check_long_tip(U, id, minLongUntig, maxShortUntig, ShortUntigRate, mainLen) \
-                        ((!((U).a[(id)].circ)) \
-                        && \
-                        ((EvaluateLen((U), (id)) > (maxShortUntig))\
-                        ||\
-                        ((long_tip((U), (id), (minLongUntig)))\
-                        &&\
-                        (EvaluateLen((U), (id)) > (ShortUntigRate)*(mainLen)))))
 #define Get_vis(visit, v, d) (((visit)[(v)>>1])&(((((v)<<(d))&1)+1)))
 #define Set_vis(visit, v, d) (((visit)[(v)>>1])|=(((((v)<<(d))&1)+1)))

@@ -479,14 +474,10 @@ typedef struct {

 void init_R_to_U(R_to_U* x, uint64_t len);
 void destory_R_to_U(R_to_U* x);
-void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig);
+void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
 void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
 void transfor_R_to_U(R_to_U* x);
 void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
-void clean_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
-long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
-R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, 
-uint32_t miniBiGraph, float chimeric_rate, int is_final_clean);
 int asg_pop_bubble_primary(asg_t *g, int max_dist);
 long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);

@@ -741,13 +732,12 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 	stops_threshold, b_0) == LOOP)
 	{
 		return UNAVAILABLE;
-	} 
+	}
 	if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &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;
@@ -768,7 +758,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)

 	uint32_t max_count = 0, min_count = 0;
 	ma_utg_t *node_min = NULL, *node_max = NULL;
-
 	if(ug != NULL)
 	{
 		/*****************************label all unitigs****************************************/
@@ -779,12 +768,11 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 			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);
+				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****************************************/

-
 		///each unitig
 		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
 		{
@@ -818,7 +806,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 				}
 			}
 		}
-
 		/*****************************label all unitigs****************************************/
 		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
 		{
@@ -831,7 +818,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 			}
 		}
 		/*****************************label all unitigs****************************************/
-		
 	}
 	else
 	{
@@ -839,7 +825,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 		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);
+			set_R_to_U(ruIndex, qn, 1, 1, &(read_sg->seq[qn].c));
 		}
 		/*****************************label all reads****************************************/

@@ -894,8 +880,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 	return NON_PLOID;
 }

-
-
 inline uint32_t check_different_haps_naive(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)
@@ -1045,10 +1029,8 @@ uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp*
 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, long long bubble_dist, uint32_t is_bubble_check, 
 uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
-void deduplicate(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 resolve_tangle);
 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);
+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);
 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 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);
@@ -1059,9 +1041,77 @@ uint32_t is_bubble_check, uint32_t is_primary_check);
 uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
 R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
 uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, 
-uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop);
+uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes);
 int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);

+
+typedef struct{
+    double weight;
+    uint32_t uID:31, del:1;
+    uint64_t dis;
+	uint64_t occ;
+	///uint64_t occ:63, scaff:1;
+    ///uint32_t enzyme;
+} hc_edge;
+
+typedef struct{
+    kvec_t(hc_edge) e;
+    kvec_t(hc_edge) f;//forbiden
+} hc_linkeage;
+
+typedef struct{
+	uint64_t beg, end;
+}bed_interval;
+
+typedef struct{
+	size_t n, m;
+	bed_interval* a;
+}bed_in;
+
+typedef struct{
+    kvec_t(hc_linkeage) a;
+    kvec_t(uint64_t) enzymes;
+	kvec_t(bed_in) bed;
+    uint32_t* u_idx;
+	uint64_t r_num;
+} hc_links;
+
+typedef struct{
+    ///kvec_t(hc_edge) a;
+    size_t n, m; 
+    hc_edge *a;
+}hc_edge_warp;
+
+void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
+void destory_hc_links(hc_links* link);
+void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
+long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
+R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, 
+uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop, 
+float drop_ratio, hc_links* link);
+void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
+ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
+long long bubble_dist, 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,
+kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link);
+void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg);
+ma_ug_t* copy_untig_graph(ma_ug_t *src);
+ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, 
+uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, 
+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, int is_bench);
+asg_t* copy_read_graph(asg_t *src);
+ma_ug_t *ma_ug_gen(asg_t *g);
+void ma_ug_destroy(ma_ug_t *ug);
+
+inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e)
+{
+    if(a_s > b_e || b_s > a_e) return 0;
+    if(i_s) (*i_s) = a_s >= b_s? a_s : b_s; ///MAX(a_s, b_s);
+    if(i_e) (*i_e) = a_e <= b_e? a_e : b_e; ///MIN(a_e, b_e);
+    return 1;
+}
+
 #define JUNK_COV 5
 #define DISCARD_RATE 0.8

@@ -46,6 +46,7 @@ void destory_All_reads(All_reads* r)
 		if (r->read_sperate[i]) free(r->read_sperate[i]);
 		if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer);
 		if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer);
+		///if (r->pb_regions) kv_destroy(r->pb_regions[i].a);
 	}
 	free(r->paf);
 	free(r->reverse_paf);
@@ -55,6 +56,7 @@ void destory_All_reads(All_reads* r)
 	free(r->name_index);
 	free(r->read_length);
 	free(r->trio_flag);
+	///if (r->pb_regions) free(r->pb_regions);
 }

 void write_All_reads(All_reads* r, char* read_file_name)
@@ -192,6 +194,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
 		r->second_round_cigar[i].lost_base_length = r->cigars[i].lost_base_length = 0;
 		r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
 	}
+	///r->pb_regions = NULL;

    free(index_name);    
    fclose(fp);
@@ -263,6 +266,8 @@ void malloc_All_reads(All_reads* r)
 	r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
 	r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
 	r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
+	///r->pb_regions = (kvec_t_u64_warp*)malloc(r->total_reads*sizeof(kvec_t_u64_warp));
+
 	for (i = 0; i < (long long)r->total_reads; i++)
 	{
 		r->second_round_cigar[i].size = r->cigars[i].size = 0;
@@ -274,6 +279,7 @@ void malloc_All_reads(All_reads* r)
 		r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
 		init_ma_hit_t_alloc(&(r->paf[i]));
 		init_ma_hit_t_alloc(&(r->reverse_paf[i]));
+		///kv_init(r->pb_regions[i].a);
 	}

 	r->name = (char*)malloc(sizeof(char)*r->total_name_length);
@@ -332,93 +338,6 @@ void init_UC_Read(UC_Read* r)
 	}
 }

-void recover_UC_Read_sub_region_begin_end(char* r, long long start_pos, long long length, uint8_t strand,
-										  All_reads* R_INF, long long ID, int extra_begin, int extra_end)
-{
-	long long readLen = Get_READ_LENGTH((*R_INF), ID);
-	uint8_t* src = Get_READ((*R_INF), ID);
-
-	long long i;
-	long long copyLen;
-	long long end_pos = start_pos + length - 1;
-
-	if (strand == 0)
-	{
-		i = start_pos;
-		copyLen = 0;
-
-		long long initLen = start_pos % 4;
-
-		if (initLen != 0)
-		{
-			memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
-			copyLen = copyLen + 4 - initLen;
-			i = i + copyLen;
-		}
-		while (copyLen < length)
-		{
-			memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
-			copyLen = copyLen + 4;
-			i = i + 4;
-		}
-
-		if (R_INF->N_site[ID])
-		{
-			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
-			{
-				if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos)
-				{
-					r[R_INF->N_site[ID][i] - start_pos] = 'N';
-				}
-				else if((long long)R_INF->N_site[ID][i] > end_pos)
-				{
-					break;
-				}
-			}
-		}
-	}
-	else
-	{
-		start_pos = readLen - start_pos - 1;
-		end_pos = readLen - end_pos - 1;
-
-		///start_pos > end_pos
-		i = start_pos;
-		copyLen = 0;
-		long long initLen = (start_pos + 1) % 4;
-
-		if (initLen != 0)
-		{
-			memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
-			copyLen = copyLen + initLen;
-			i = i - initLen;
-		}
-
-		while (copyLen < length)
-		{
-			memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
-			copyLen = copyLen + 4;
-			i = i - 4;
-		}
-
-		if (R_INF->N_site[ID])
-		{
-			long long offset = readLen - start_pos - 1;
-
-			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
-			{
-				if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
-				{
-					r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
-				}
-				else if((long long)R_INF->N_site[ID][i] > start_pos)
-				{
-					break;
-				}
-			}
-		}
-	}
-}

 void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
 {
@@ -509,6 +428,104 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
 }


+void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
+{
+	i_r->length = length;i_r->RID = ID;
+	if (i_r->length + 8 > i_r->size)
+	{
+		i_r->size = i_r->length + 4;
+		i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size));
+	}
+	char* r = i_r->seq;
+	long long readLen = Get_READ_LENGTH((*R_INF), ID);
+	uint8_t* src = Get_READ((*R_INF), ID);
+
+	long long i;
+	long long copyLen;
+	long long end_pos = start_pos + length - 1;
+
+	if (strand == 0)
+	{
+		i = start_pos;
+		copyLen = 0;
+
+		long long initLen = start_pos % 4;
+
+		if (initLen != 0)
+		{
+			memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
+			copyLen = copyLen + 4 - initLen;
+			i = i + copyLen;
+		}
+
+		while (copyLen < length)
+		{
+			memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
+			copyLen = copyLen + 4;
+			i = i + 4;
+		}
+
+		if (R_INF->N_site[ID])
+		{
+			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
+			{
+				if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos)
+				{
+					r[R_INF->N_site[ID][i] - start_pos] = 'N';
+				}
+				else if((long long)R_INF->N_site[ID][i] > end_pos)
+				{
+					break;
+				}
+			}
+		}
+	}
+	else
+	{
+		start_pos = readLen - start_pos - 1;
+		end_pos = readLen - end_pos - 1;
+
+		///start_pos > end_pos
+		i = start_pos;
+		copyLen = 0;
+		long long initLen = (start_pos + 1) % 4;
+
+		if (initLen != 0)
+		{
+			memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
+			copyLen = copyLen + initLen;
+			i = i - initLen;
+		}
+
+		while (copyLen < length)
+		{
+			memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
+			copyLen = copyLen + 4;
+			i = i - 4;
+		}
+
+		if (R_INF->N_site[ID])
+		{
+			long long offset = readLen - start_pos - 1;
+
+			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
+			{
+				if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
+				{
+					r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
+				}
+				else if((long long)R_INF->N_site[ID][i] > start_pos)
+				{
+					break;
+				}
+			}
+		}
+	}
+}
+
+
+
+
 void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID)
 {
 	r->length = Get_READ_LENGTH((*R_INF), ID);
@@ -667,3 +684,54 @@ void reverse_complement(char* pattern, uint64_t length)
 		pattern[end] = RC_CHAR(pattern[end]);
 	}
 }
+
+
+void init_Debug_reads(Debug_reads* x, const char* file)
+{
+	int nameLen, i, bufLen = 1000;
+	if((uint64_t)(bufLen) < strlen(file) + 50) bufLen = strlen(file) + 50;
+	char* Name_Buffer = (char*)malloc(sizeof(char)*bufLen);
+	fprintf(stderr, "Queried debugging reads at: %s\n", file);
+
+	x->fp = fopen(file,"r");
+	x->query_num = 0;
+	
+	while(fgets(Name_Buffer, bufLen, x->fp))
+	{
+		x->query_num++;
+	}
+	x->read_name = (char**)malloc(sizeof(char*)*x->query_num);
+	x->candidate_count = (kvec_t_u64_warp*)malloc(sizeof(kvec_t_u64_warp)*x->query_num);
+	fseek(x->fp, 0, SEEK_SET);
+
+	i = 0;
+	while(fgets(Name_Buffer, bufLen, x->fp))
+	{
+		nameLen = strlen(Name_Buffer) - 1;
+		x->read_name[i] = (char*)malloc(sizeof(char)*(nameLen+1));
+		memcpy(x->read_name[i], Name_Buffer, sizeof(char)*nameLen);
+		x->read_name[i][nameLen] = '\0';
+		kv_init(x->candidate_count[i].a);
+		i++;
+	}
+	
+    fclose(x->fp);
+
+	sprintf(Name_Buffer, "%s.debug.stdout", file);
+	x->fp = fopen(Name_Buffer,"w");
+	fprintf(stderr, "Print debugging information to: %s\n", Name_Buffer);
+	free(Name_Buffer);
+}
+
+void destory_Debug_reads(Debug_reads* x)
+{
+	uint64_t i;
+	for (i = 0; i < x->query_num; i++)
+	{
+		free(x->read_name[i]);
+		kv_destroy(x->candidate_count[i].a);
+	}
+	
+	free(x->read_name);
+	fclose(x->fp);
+}
@@ -22,7 +22,8 @@
 ///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID
 #define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)]
 #define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)])
-
+#define CHECK_BY_NAME(R_INF, NAME, ID) (Get_NAME_LENGTH((R_INF),(ID))==strlen((NAME)) && \
+                                        memcmp((NAME), Get_NAME((R_INF), (ID)), Get_NAME_LENGTH((R_INF),(ID))) == 0)

 extern uint8_t seq_nt6_table[256];
 extern char bit_t_seq_table[256][4];
@@ -93,6 +94,7 @@ typedef struct
 	uint32_t new_length;
 } Compressed_Cigar_record;

+
 #define AMBIGU 0
 #define FATHER 1
 #define MOTHER 2
@@ -128,6 +130,8 @@ typedef struct

    ma_hit_t_alloc* paf;
    ma_hit_t_alloc* reverse_paf;
+
+    ///kvec_t_u64_warp* pb_regions;
 } All_reads;

 extern All_reads R_INF;
@@ -140,6 +144,15 @@ typedef struct
 	long long RID;
 } UC_Read;

+typedef struct
+{
+    char** read_name;
+    uint64_t query_num;
+    kvec_t_u64_warp* candidate_count;
+    FILE* fp;
+    pthread_mutex_t OutputMutex;
+} Debug_reads;
+
 void init_All_reads(All_reads* r);
 void malloc_All_reads(All_reads* r);
 void ha_insert_read_len(All_reads *r, int read_len, int name_len);
@@ -154,5 +167,7 @@ void write_All_reads(All_reads* r, char* read_file_name);
 int load_All_reads(All_reads* r, char* read_file_name);
 void destory_All_reads(All_reads* r);
 int destory_read_bin(All_reads* r);
-
+void init_Debug_reads(Debug_reads* x, const char* file);
+void destory_Debug_reads(Debug_reads* x);
+void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID);
 #endif
@@ -7,6 +7,7 @@
 #include "Correct.h"
 #include "kthread.h"
 #include "kdq.h"
+#include "hic.h"

 KDQ_INIT(uint64_t)

@@ -3291,7 +3292,7 @@ int asg_pop_bubble_purge_graph(asg_t *purge_g, int max_dist)
 		for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs
 			if (!av[i].del) ++n_arc;
 		if (n_arc > 1)
-			n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1);
+			n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1, NULL, NULL);
 	}
 	free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
 	if (n_pop) asg_cleanup(purge_g);
@@ -3925,10 +3926,57 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg
    
 }

+
+
+void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, hap_overlaps* t)
+{
+    uint32_t i = 0, k = 0, rId_0, rId_1, pre_0, pre_1, b_0 = t->xUid, b_1 = t->yUid;
+    uint64_t d = RC_2;
+    ma_utg_t* u_b_0 = &(ug->u.a[b_0]);
+    ma_utg_t* u_b_1 = &(ug->u.a[b_1]);
+    if(u_b_0->n == 0) return;
+    if(u_b_1->n == 0) return;
+
+    for (i = t->x_beg_id, pre_0 = (uint32_t)-1; i < t->x_end_id; i++)
+    {
+        rId_0 = u_b_0->a[i]>>33;
+        if(link->u_idx[rId_0] == (uint32_t)-1) continue;
+        if(pre_0 == link->u_idx[rId_0]) continue;
+        pre_0 = link->u_idx[rId_0];
+        
+        for (k = t->y_beg_id, pre_1  = (uint32_t)-1; k < t->y_end_id; k++)
+        {
+            rId_1 = u_b_1->a[k]>>33;
+            if(link->u_idx[rId_1] == (uint32_t)-1) continue;
+            if(pre_1 == link->u_idx[rId_1]) continue;
+            pre_1 = link->u_idx[rId_1];
+            push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d);
+            push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d);
+        }
+    }
+
+}
+
+
+void collect_reverse_unitigs_purge(buf_t* b_0, hc_links* link, ma_ug_t *ug, hap_overlaps_list* all_ovlp)
+{
+    if(b_0->b.n <= 1) return;
+    uint32_t k;
+    int index = 0;
+    for (k = 0; k < b_0->b.n - 1; k++)
+    {
+        index = get_specific_hap_overlap(&(all_ovlp->x[b_0->b.a[k]>>1]), b_0->b.a[k]>>1, b_0->b.a[k+1]>>1);
+        if(index == -1) continue;
+        collect_reverse_unitig_pair(link, ug, &(all_ovlp->x[b_0->b.a[k]>>1].a.a[index]));
+    }
+}
+
+
 void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp,
 R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, 
 uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, 
-kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit)
+kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit,
+hc_links* link)
 {
    uint32_t v, n_vtx = purge_g->n_seq * 2, beg, end;
    long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen;
@@ -3948,6 +3996,8 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg
        {
            continue;
        }
+
+        if(link) collect_reverse_unitigs_purge(&b_0, link, ug, all_ovlp);
        purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut, 
        read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit);
    }
@@ -4168,11 +4218,10 @@ uint32_t minLen, double purge_threshold)
    return 0;
 }

-
 void purge_dups(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, kvec_asg_arc_t_warp* edge, float density, 
 uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, 
-uint32_t just_contain, uint32_t just_coverage)
+uint32_t just_contain, uint32_t just_coverage, hc_links* link)
 {
    asg_t *purge_g = NULL;
    purge_g = asg_init();
@@ -4229,7 +4278,7 @@ uint32_t just_contain, uint32_t just_coverage)
        for (i = 0, offset = 0; i < reads->n; i++)
        {
            rId = reads->a[i]>>33;
-            set_R_to_U(ruIndex, rId, uId, 1);
+            set_R_to_U(ruIndex, rId, uId, 1, &(read_g->seq[rId].c));

            position_index[rId] = offset;
            position_index[rId] = position_index[rId] << 32;
@@ -4301,6 +4350,7 @@ uint32_t just_contain, uint32_t just_coverage)
                purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE;
                purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del = 1;
                all_ovlp.x[uId].a.a[i].status = DELETE;
+                if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
            }

            if(all_ovlp.x[uId].a.a[i].type == XCY)
@@ -4309,6 +4359,7 @@ uint32_t just_contain, uint32_t just_coverage)
                purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE;
                purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del = 1;
                all_ovlp.x[uId].a.a[i].status = DELETE;
+                if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
            }
            ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
        }
@@ -4360,7 +4411,7 @@ uint32_t just_contain, uint32_t just_coverage)

        link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
        &(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex),
-        max_hang, min_ovlp, edge, hap_buf.buf[0].visit);
+        max_hang, min_ovlp, edge, hap_buf.buf[0].visit, link);
    }

    for (v = 0; v < all_ovlp.num; v++)
@@ -15,7 +15,7 @@
 void purge_dups(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, kvec_asg_arc_t_warp* edge, float density, 
 uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, 
-uint32_t just_contain, uint32_t just_coverage);
+uint32_t just_contain, uint32_t just_coverage, hc_links* link);
 void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
 uint32_t is_circle, uint64_t* rLen);
 void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);
@@ -8,7 +8,7 @@ cd hifiasm && make
 # Run on test data (use -f0 for small datasets)
 wget https://github.com/chhylp123/hifiasm/releases/download/v0.7/chr11-2M.fa.gz
 ./hifiasm -o test -t4 -f0 chr11-2M.fa.gz 2> test.log
-awk '/^S/{print ">"$1;print $2}' test.p_ctg.gfa > test.p_ctg.fa  # get primary contigs in FASTA
+awk '/^S/{print ">"$2;print $3}' test.p_ctg.gfa > test.p_ctg.fa  # get primary contigs in FASTA

 # Assemble inbred/homozygous genomes (-l0 disables duplication purging)
 hifiasm -o CHM13.asm -t32 -l0 CHM13-HiFi.fa.gz 2> CHM13.asm.log
@@ -30,6 +30,29 @@ produce primary/alternate assemblies of quality competitive with the best
 assemblers. It also introduces a new graph binning algorithm and achieves
 the best haplotype-resolved assembly given trio data.

+## Why Hifiasm?
+
+* Hifiasm delivers high-quality assemblies. It tends to generate longer contigs
+  and resolve more segmental duplications than other assemblers.
+
+* Given sequence reads from the parents, hifiasm can produce overall the best
+  haplotype-resolved assembly so far. It is the assembler of choice by the
+  [Human Pangenome Project][hpp] for the first batch of samples.
+
+* Hifiasm can purge duplications between haplotigs without relying on
+  third-party tools such as purge\_dups. Hifiasm does not need polishing tools
+  like pilon or racon, either. This simplifies the assembly pipeline and saves
+  running time.
+
+* Hifiasm is fast. It can assemble a human genome in half a day and assemble a
+  ~30Gb redwood genome in three days. No genome is too large for hifiasm.
+
+* Hifiasm is trivial to install and easy to use. It does not required python,
+  R or C++11 compilers and can be compiled into a single executable. The
+  default setting works well with a variety of genomes.
+
+[hpp]: https://humanpangenome.org
+
 ## Usage

 A typical hifiasm command line looks like:
@@ -144,12 +167,12 @@ redwood genome in a few days on a single machine. For trio binning assembly:
 [NA12891-data]: https://www.ebi.ac.uk/ena/data/view/ERR194160
 [NA12892-data]: https://www.ebi.ac.uk/ena/data/view/ERR194161

-Except NA12878, the assemblies above were produced by hifiasm v0.7 and can be
+Except NA12878, the assemblies above were produced by hifiasm v0.12 and can be
 downloaded at
 ```txt
-ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/submission/v0.7/
+ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/submission/hifiasm-0.12/
 ```
-NA12878 was assembled with a more recent version of hifiasm and is available at
+NA12878 was assembled with an older version of hifiasm and is available at
 ```txt
 ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/NA12878-r253/
 ```
@@ -169,4 +192,4 @@ page](https://github.com/chhylp123/hifiasm/issues).

 ## Limitations

-1. Purging haplotig duplications may introduce misassemblies.
+1. Purging haplotig duplications may introduce misassemblies.
@@ -67,7 +67,7 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ..
 {
 	va_list ap;
 	FILE *fp;
-	uint32_t t[3];
+	uint32_t t[3], f_tmp = 0;
 	char magic[4];
 	int i, j, absent, min_cnt = 0, mid_cnt = 0, mode_err = 0;
 	uint64_t mask = (1ULL<<YAK_COUNTER_BITS) - 1, n_ins = 0, n_new = 0;
@@ -92,33 +92,36 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ..
 		fclose(fp);
 		return 0;
 	}
-	fread(t, 4, 3, fp);
+	f_tmp += fread(t, 4, 3, fp);
 	if (t[2] != YAK_COUNTER_BITS) {
 		fprintf(stderr, "ERROR: saved counter bits: %d; compile-time counter bits: %d\n", t[2], YAK_COUNTER_BITS);
 		fclose(fp);
 		return 0;
 	}
-
+	///t[0] = k; t[1] = pre, t[2] = YAK_COUNTER_BITS;
 	ch = ch0 == 0? yak_ch_init(t[0], t[1]) : ch0;
 	assert((int)t[0] == ch->k && (int)t[1] == ch->pre);
 	for (i = 0; i < 1<<ch->pre; ++i) {
 		yak_ht_t *h = ch->h[i].h;
-		fread(t, 4, 2, fp);
+		f_tmp += fread(t, 4, 2, fp);
+		///t[0] = kh_capacity(h), t[1] = kh_size(h);
 		if (ch0 == 0) yak_ht_resize(h, t[0]);
 		for (j = 0; j < (int)t[1]; ++j) {
 			uint64_t key;
-			fread(&key, 8, 1, fp);
+			f_tmp += fread(&key, 8, 1, fp);
 			if (mode == YAK_LOAD_ALL) {
 				++n_ins;
 				yak_ht_put(h, key, &absent);
 				if (absent) ++n_new;
 			} else if (mode == YAK_LOAD_TRIOBIN1 || mode == YAK_LOAD_TRIOBIN2) {
 				int cnt = key & mask, x, shift = mode == YAK_LOAD_TRIOBIN1? 0 : 2;
+				//1. filter singleton k-mer; 2. label non-repeat and repeat
 				if (cnt >= mid_cnt) x = 2<<shift;
 				else if (cnt >= min_cnt) x = 1<<shift;
 				else x = -1;
 				if (x >= 0) {
 					khint_t k;
+					///no need cnt at all
 					key = (key & ~mask) | x;
 					++n_ins;
 					k = yak_ht_put(h, key, &absent);
@@ -28,8 +28,8 @@ typedef struct {
 } seed1_t;

 struct ha_abuf_s {
-	uint64_t n_a, m_a;
-	uint32_t old_mz_m;
+	uint64_t n_a, m_a;///number of anchors (seed positions)
+	uint32_t old_mz_m;///number of seeds
 	ha_mz1_v mz;
 	seed1_t *seed;
 	anchor1_t *a;
@@ -50,17 +50,16 @@ uint64_t ha_abuf_mem(const ha_abuf_t *ab)
 	return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t);
 }

-static int ha_ov_type(const overlap_region *r, uint32_t len)
+int ha_ov_type(const overlap_region *r, uint32_t len)
 {
 	if (r->x_pos_s == 0 && r->x_pos_e == len - 1) return 2; // contained in a longer read
 	else if (r->x_pos_s > 0 && r->x_pos_e < len - 1) return 3; // containing a shorter read
 	else return r->x_pos_s == 0? 0 : 1;
 }

-void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain)
+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)
 {
-	extern void *ha_flt_tab;
-	extern ha_pt_t *ha_idx;
 	uint32_t i, rlen;
 	uint64_t k, l;
 	double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
@@ -74,7 +73,8 @@ 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);
+	ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct);
+	// minimizer of queried read
 	if (ab->mz.m > ab->old_mz_m) {
 		ab->old_mz_m = ab->mz.m;
 		REALLOC(ab->seed, ab->old_mz_m);
@@ -93,6 +93,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	}
 	for (i = 0, k = 0; i < ab->mz.n; ++i) {
 		int j;
+		///z is one of the minimizer
 		ha_mz1_t *z = &ab->mz.a[i];
 		seed1_t *s = &ab->seed[i];
 		for (j = 0; j < s->n; ++j) {
@@ -116,6 +117,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		}
 	}

+
 	// copy over to _cl_
 	if (ab->m_a >= (uint64_t)cl->size) {
 		cl->size = ab->m_a;
@@ -131,7 +133,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	}
 	cl->length = ab->n_a;

-	calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain);
+	calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain, f_cigar);

 	#if 0
 	if (overlap_list->length > 0) {
@@ -172,5 +174,132 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		}
 	}

+	///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+}
+
+
+void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
+{
+	uint64_t j = 0, inner_j = 0;
+	while (j < overlap_list->length && inner_j < paf->length)
+    {
+        if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn)
+        {
+            j++;
+        }
+        else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn)
+        {
+            inner_j++;
+        }
+        else
+        {
+            if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev)
+            {
+				overlap_list->list[j].is_match = 1;
+            }
+            j++;
+            inner_j++;
+        }
+    }
+}
+
+
+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)
+{
+	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);
+
+	if(ha_idx_hp)
+	{
+		uint32_t i, k, y_id, overlapLen, max_i;
+		int shared_seed;
+		overlap_region t;
+		overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
+		ma_hit_sort_tn(paf->buffer, paf->length);
+		ma_hit_sort_tn(rev_paf->buffer, rev_paf->length);
+		lable_matched_ovlp(overlap_list, paf);
+		lable_matched_ovlp(overlap_list, rev_paf);
+
+		for (i = 0, k = 0; i < overlap_list->length; ++i) 
+		{
+			if(overlap_list->list[i].is_match == 1)
+			{
+				if(k != i)
+				{
+					t = overlap_list->list[k];
+					overlap_list->list[k] = overlap_list->list[i];
+					overlap_list->list[i] = t;
+					overlap_list->list[k].is_match = 0;
+				}
+				k++;
+			}
+		}
+		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);	
+		
+		if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
+		{
+			overlap_list->list = (overlap_region*)realloc(overlap_list->list, 
+				sizeof(overlap_region)*(overlap_list->length + overlap_list_hp->length));
+			memset(overlap_list->list + overlap_list->size, 0, sizeof(overlap_region)*
+				(overlap_list->length + overlap_list_hp->length - overlap_list->size));
+			overlap_list->size = overlap_list->length + overlap_list_hp->length;
+		}
+		
+		for (i = 0, k = overlap_list->length; i < overlap_list_hp->length; i++, k++)
+		{
+			t = overlap_list->list[k];
+			overlap_list->list[k] = overlap_list_hp->list[i];
+			overlap_list_hp->list[i] = t;
+		}
+		overlap_list->length = k;
+		
+		overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
+
+		i = k = 0;
+		while (i < overlap_list->length)
+		{
+			y_id = overlap_list->list[i].y_id;
+			shared_seed = overlap_list->list[i].shared_seed;
+			overlapLen = overlap_list->list[i].overlapLen;
+			max_i = i;
+			i++;
+			while (i < overlap_list->length && overlap_list->list[i].y_id == y_id)
+			{
+				if((overlap_list->list[i].shared_seed > shared_seed) || 
+				  ((overlap_list->list[i].shared_seed == shared_seed) && (overlap_list->list[i].overlapLen <= overlapLen)))
+				{
+					y_id = overlap_list->list[i].y_id;
+					shared_seed = overlap_list->list[i].shared_seed;
+					overlapLen = overlap_list->list[i].overlapLen;
+					max_i = i;
+				}
+				i++;
+			}
+
+			if(k != max_i)
+			{
+				t = overlap_list->list[k];
+				overlap_list->list[k] = overlap_list->list[max_i];
+				overlap_list->list[max_i] = t;
+			}
+			k++;
+		}
+
+		overlap_list->length = k;
+	}
+	
 	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
 }
+
+void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
+{
+	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+}
@@ -0,0 +1,72 @@
+#ifndef __HIC__
+#define __HIC__
+#include <stdint.h>
+#include "Overlaps.h"
+
+#define kdq_clear(q) ((q)->count = (q)->front = 0)
+#define kv_malloc(v, s) ((v).n = 0, (v).m = (s), MALLOC((v).a, (s)))
+#define RC_0 0
+#define RC_1 1
+#define RC_2 2
+
+hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir);
+void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d);
+void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link);
+void hic_benchmark(ma_ug_t *ug, asg_t* read_g);
+
+typedef struct {
+    double w;
+    uint32_t id, occ;
+    ///uint32_t *bid, bid_n;
+    ma_utg_t *u;
+    uint64_t l_d, r_d;
+}chain_hic_w_type;
+
+typedef struct {
+    size_t n, m;
+    chain_hic_w_type* a;
+    uint32_t max_bub_id;
+    uint32_t *chain_idx, u_n;
+}chain_hic_warp;
+
+typedef struct {
+    long long g_occ, b_occ;
+    uint64_t id;
+    uint8_t del;
+}chain_w_type;
+
+typedef struct {
+    uint32_t* index, round_id, n_round;
+    ma_ug_t* ug;
+    kvec_t(uint32_t) list;
+    kvec_t(uint32_t) num;
+    kvec_t(uint64_t) pathLen;
+    kvec_t(uint64_t) b_s_idx;
+    uint64_t s_bub, f_bub, b_bub, b_end_bub, tangle_bub, cross_bub, mess_bub;
+    uint32_t check_het;
+    asg_t *b_g;
+    ma_ug_t* b_ug;
+    kvec_t(chain_w_type) chain_weight;
+    chain_hic_warp c_w;
+} bubble_type;
+#define P_het(B) ((B).num.n)
+#define M_het(B) ((B).num.n + 1)
+// #define IF_BUB(ID, B) ((B).index[(ID)] < (B).num.n)
+// #define IF_HET(ID, B) ((B).index[(ID)] == (B).num.n)
+// #define IF_HOM(ID, B) ((B).index[(ID)] > (B).num.n)
+#define IF_BUB(ID, B) ((B).index[(ID)] < (B).f_bub+1)
+#define IF_HET(ID, B) ((B).index[(ID)] == (B).f_bub+1)
+#define IF_HOM(ID, B) ((B).index[(ID)] > (B).f_bub+1)
+#define Get_bub_num(RECORD) ((RECORD).num.n-1)
+void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** a, uint32_t* n, uint64_t* pathBase);
+int load_hc_links(hc_links* link, const char *fn);
+void write_hc_links(hc_links* link, const char *fn);
+void destory_bubbles(bubble_type* bub);
+void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link);
+void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub);
+uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het);
+void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het);
+void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end);
+void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ);
+
+#endif
@@ -111,6 +111,11 @@ assembly.
 .BI -r \ INT
 Rounds of haplotype-aware error corrections [3]. This option affects all outputs of hifiasm.

+.TP
+.BI --min-hist-cnt \ INT
+When analyzing the k-mer spectrum, ignore counts below
+.IR INT .
+
 .SS Assembly options

 .TP
@@ -198,7 +203,14 @@ Min and max coverage cutoff of primary contigs.
 Keep contigs with coverage in this range at p_ctg.gfa. 
 Inferred automatically in default.
 If INT2 is not specified, it is set to infinity. 
-Set -1 to disable 
+Set -1 to disable.
+
+.TP
+.BI --lowQ \ INT
+Output contig regions with >=INT% inconsistency to the bed file 
+with suffix
+.B lowQ.bed
+[70]. Set 0 to disable. 

 .SS Trio-partition options

@@ -264,6 +276,12 @@ Min number of overlapped reads for duplicate haplotigs that should be purged [1]
 Coverage upper bound of Purge-dups, which is inferred automatically in default.
 If the coverage of a contig is higher than this bound, don't apply Purge-dups.

+.TP
+.BI --high-het \ INT
+Enable this mode for high heterozygosity sample, which will increase running time. 
+For ordinary samples, no need to enable this mode [experimental, not stable].
+
+
 .SS Debugging options

 .TP 10
@@ -1,4 +1,5 @@
 #include <stdio.h>
+#include <assert.h>
 #include "htab.h"

 static void ha_hist_line(int c, int x, int exceed, int64_t cnt)
@@ -11,16 +12,19 @@ static void ha_hist_line(int c, int x, int exceed, int64_t cnt)
 	fprintf(stderr, " %lld\n", (long long)cnt);
 }

-int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het)
+int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het)
 {
 	const int hist_max = 100;
 	int i, start, low_i, max_i, max2_i, max3_i;
 	int64_t max, max2, max3, min;

-	// find the low point from the left
+	// determine the start point
+	assert(n_cnt > start_cnt);
 	*peak_het = -1;
 	start = cnt[1] > 0? 1 : 2;
-	low_i = start;
+
+	// find the low point from the left
+	low_i = start > start_cnt? start : start_cnt;
 	for (i = low_i + 1; i < n_cnt; ++i)
 		if (cnt[i] > cnt[i-1]) break;
 	low_i = i - 1;
@@ -35,6 +35,9 @@ const unsigned char seq_nt4_table[256] = { // translate ACGT to 0123

 void *ha_flt_tab;
 ha_pt_t *ha_idx;
+void *ha_flt_tab_hp;
+ha_pt_t *ha_idx_hp;
+void *ha_ct_table;

 /***************************
 * Yak specific parameters *
@@ -46,6 +49,7 @@ typedef struct {
 	int32_t pre;
 	int32_t n_thread;
 	int64_t chunk_size;
+	int adaLen;
 } yak_copt_t;

 void yak_copt_init(yak_copt_t *o)
@@ -71,7 +75,7 @@ typedef struct {
 	int n_shift, n_hashes;
 	uint8_t *b;
 } yak_bf_t;
-
+///in most cases, n_shift = 25, n_hashes = 4
 yak_bf_t *yak_bf_init(int n_shift, int n_hashes)
 {
 	yak_bf_t *b;
@@ -126,10 +130,12 @@ typedef struct {

 typedef struct {
 	int k, pre, n_hash, n_shift;
-	uint64_t tot;
+	uint64_t tot;  ///number of distinct k-mers
 	ha_ct1_t *h;
 } ha_ct_t;

+///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 37
+///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
 static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
 {
 	ha_ct_t *h;
@@ -138,12 +144,15 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
 	CALLOC(h, 1);
 	h->k = k, h->pre = pre;
 	CALLOC(h->h, 1<<h->pre);
+	///i<h->pre = 4096
+	///it seems there is a large hash table h, consisting 4096 small hash tables
 	for (i = 0; i < 1<<h->pre; ++i)
 		h->h[i].h = yak_ct_init();
+	///for 0-th counting, enter here; used for bloom filter
 	if (n_hash > 0 && n_shift > h->pre) {
 		h->n_hash = n_hash, h->n_shift = n_shift;
 		for (i = 0; i < 1<<h->pre; ++i)
-			h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash);
+			h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); ///h->n_shift = 37, h->pre = 12, h->n_hash = 4 
 	}
 	return h;
 }
@@ -173,15 +182,23 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *
 	int j, mask = (1<<h->pre) - 1, n_ins = 0;
 	ha_ct1_t *g;
 	if (n == 0) return 0;
+	///corresponding small hash index
 	g = &h->h[a[0]&mask];
 	for (j = 0; j < n; ++j) {
 		int ins = 1, absent;
+		///x is a 64-bit word, h->pre=12
+		///all elements at a have the same low 12 bits
+		///so low 12 bits are not useful
 		uint64_t x = a[j] >> h->pre;
 		khint_t k;
 		if ((a[j]&mask) != (a[0]&mask)) continue;
 		if (create_new) {
+			///for 0-th counting, g->b = NULL
 			if (g->b)
 				ins = (yak_bf_insert(g->b, x) == h->n_hash);
+			///for 0-th counting, g->b = NULL
+			///x = the high 52 bits of a[j] + low 12 bits 0
+			///the low 12 bits are used for counting
 			if (ins) {
 				k = yak_ct_put(g->h, x << YAK_COUNTER_BITS | (g->b? 1 : 0), &absent);
 				if (absent) ++n_ins;
@@ -219,6 +236,8 @@ static void worker_ct_hist(void *data, long i, int tid) // callback for kt_for()
 			++cnt[kh_key(g, k)&YAK_MAX_COUNT];
 }

+///YAK_N_COUNTS is also 4096
+///used for calculating k-mer histogram
 static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread)
 {
 	hist_aux_t a;
@@ -226,6 +245,7 @@ static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread
 	a.h = h;
 	memset(cnt, 0, YAK_N_COUNTS * sizeof(uint64_t));
 	CALLOC(a.cnt, n_thread);
+	///start 4096 threads
 	kt_for(n_thread, worker_ct_hist, &a, 1<<h->pre);
 	for (i = 0; i < YAK_N_COUNTS; ++i) cnt[i] = 0;
 	for (j = 0; j < n_thread; ++j)
@@ -265,6 +285,7 @@ static void ha_ct_shrink(ha_ct_t *h, int min, int max, int n_thread)
 	int i;
 	shrink_aux_t a;
 	a.h = h, a.min = min, a.max = max;
+	///still start 4096 threads
 	kt_for(n_thread, worker_ct_shrink, &a, 1<<h->pre);
 	for (i = 0, h->tot = 0; i < 1<<h->pre; ++i)
 		h->tot += kh_size(h->h[i].h);
@@ -306,6 +327,7 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for()
 			int absent;
 			khint_t l;
 			l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent);
+			///this should be the start index of kh_key's corresponding pos at ha_idxpos_t* a
 			kh_val(b->h, l) = b->n;
 			b->n += kh_key(g, k) & YAK_MAX_COUNT;
 		}
@@ -412,8 +434,12 @@ typedef struct {
 	ha_mz1_t *b;
 } ch_buf_t;

+///p = 12
 static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
 {
+	///assign k-mer to one of the 4096 bins
+	///using low 12 bits for assigning
+	///so all elements at b have the same low 12 bits
 	int pre = y & ((1<<p) - 1);
 	ch_buf_t *b = &buf[pre];
 	if (b->n == b->m) {
@@ -425,6 +451,7 @@ static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k

 static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
 {
+	///assign minimizer to one of 4096 bins by low 12 bits
 	int pre = y->x & ((1<<p) - 1);
 	ch_buf_t *b = &buf[pre];
 	if (b->n == b->m) {
@@ -434,13 +461,17 @@ static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
 	b->b[b->n++] = *y;
 }

+///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read
 static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf
 {
 	int i, l;
 	uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;
 	for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
 		int c = seq_nt4_table[(uint8_t)seq[i]];
+		///c = 00, 01, 10, 11
 		if (c < 4) { // not an "N" base
+			///x[0] & x[1] are the forward k-mer
+			///x[2] & x[3] are the reverse complementary k-mer
 			x[0] = (x[0] << 1 | (c&1))  & mask;
 			x[1] = (x[1] << 1 | (c>>1)) & mask;
 			x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;
@@ -483,12 +514,13 @@ KSEQ_INIT(gzFile, gzread)
 #define HAF_RS_WRITE_SEQ 0x8
 #define HAF_RS_READ      0x10
 #define HAF_CREATE_NEW   0x20
+#define HAF_SKIP_READ    0x40

 typedef struct { // global data structure for kt_pipeline()
 	const yak_copt_t *opt;
 	const void *flt_tab;
 	int flag, create_new, is_store;
-	uint64_t n_seq;
+	uint64_t n_seq; ///number of total reads
 	kseq_t *ks;
 	UC_Read ucr;
 	ha_ct_t *ct;
@@ -499,7 +531,8 @@ typedef struct { // global data structure for kt_pipeline()

 typedef struct { // data structure for each step in kt_pipeline()
 	pl_data_t *p;
-	uint64_t n_seq0;
+	uint64_t n_seq0;  ///the start index of current buffer block at R_INF
+	///sum_len = total bases, nk = number of k-mers
 	int n_seq, m_seq, sum_len, nk;
 	int *len;
 	char **seq;
@@ -514,15 +547,17 @@ static void worker_for_insert(void *data, long i, int tid) // callback for kt_fo
 	ch_buf_t *b = &s->buf[i];
 	if (s->p->pt)
 		b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b);
-	else
+	else///for 0-th count, go into here
 		b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);
 }

 static void worker_for_mz(void *data, long i, int tid)
 {
 	st_data_t *s = (st_data_t*)data;
+	///get the corresponding minimzer vector of this read
 	ha_mz1_v *b = &s->mz_buf[tid];
 	s->mz_buf[tid].n = 0;
+	///s->p->opt->w = 51, s->p->opt->k
 	ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab);
 	s->mz[i].n = s->mz[i].m = b->n;
 	MALLOC(s->mz[i].a, b->n);
@@ -540,6 +575,11 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		s->n_seq0 = p->n_seq;
 		if (p->rs_in && (p->flag & HAF_RS_READ)) {
 			while (p->n_seq < p->rs_in->total_reads) {
+				if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
+				{
+					++p->n_seq;
+					continue;
+				}
 				int l;
 				recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
 				l = p->ucr.length;
@@ -559,12 +599,15 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 			}
 		} else {
 			while ((ret = kseq_read(p->ks)) >= 0) {
-				int l = p->ks->seq.l;
+				int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);
+				if(l <= 0) continue;
+
 				if (p->n_seq >= 1<<28) {
 					fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
 					exit(1);
 				}
 				if (p->rs_out) {
+					///for 0-th count, just insert read length to R_INF, instead of read
 					if (p->flag & HAF_RS_WRITE_LEN) {
 						assert(p->n_seq == p->rs_out->total_reads);
 						ha_insert_read_len(p->rs_out, l, p->ks->name.l);
@@ -572,23 +615,25 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 						int i, n_N;
 						assert(l == (int)p->rs_out->read_length[p->n_seq]);
 						for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
-							if (seq_nt4_table[(uint8_t)p->ks->seq.s[i]] >= 4)
+							if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)
 								++n_N;
-						ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s, l, &p->rs_out->N_site[p->n_seq], n_N);
+						ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);
 						memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
 					}
 				}
+				///for 0-th count, insert both seq and length to local block
 				if (s->n_seq == s->m_seq) {
 					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
 					REALLOC(s->len, s->m_seq);
 					REALLOC(s->seq, s->m_seq);
 				}
 				MALLOC(s->seq[s->n_seq], l);
-				memcpy(s->seq[s->n_seq], p->ks->seq.s, l);
+				memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);
 				s->len[s->n_seq++] = l;
 				++p->n_seq;
 				s->sum_len += l;
 				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
+				///p->opt->chunk_size is the block max size
 				if (s->sum_len >= p->opt->chunk_size)
 					break;
 			}
@@ -596,18 +641,24 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		if (s->sum_len == 0) free(s);
 		else return s;
 	} else if (step == 1) { // step 2: extract k-mers
+		///s is the block of reads
 		st_data_t *s = (st_data_t*)in;
+		///for 0-th counting, n_pre = 4096
 		int i, n_pre = 1<<p->opt->pre, m;
 		// allocate the k-mer buffer
 		CALLOC(s->buf, n_pre);
 		m = (int)(s->nk * 1.2 / n_pre) + 1;
+		//pre-allocate memory for each of 4096 buffer
 		for (i = 0; i < n_pre; ++i) {
 			s->buf[i].m = m;
+			///for 0-th counting, p->pt = NULL
 			if (p->pt) MALLOC(s->buf[i].b, m);
 			else MALLOC(s->buf[i].a, m);
 		}
 		// fill the buffer
+		///for 0-th counting, p->opt->w == 1
 		if (p->opt->w == 1) { // enumerate all k-mers
+			///scan all reads
 			for (i = 0; i < s->n_seq; ++i) {
 				if (p->opt->is_HPC)
 					count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
@@ -618,18 +669,21 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		} else { // minimizers only
 			uint32_t j;
 			// compute minimizers
+			// s->n_seq is how many reads at this buffer
+			// s->mz && s->mz_buf are lists of minimzer vectors
 			CALLOC(s->mz, s->n_seq);
 			CALLOC(s->mz_buf, 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)
 				free(s->mz_buf[i].a);
 			free(s->mz_buf);
 			// insert minimizers
-			if (p->pt) {
+			if (p->pt) {///insert whole minimizer
 				for (i = 0; i < s->n_seq; ++i)
 					for (j = 0; j < s->mz[i].n; ++j)
 						pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
-			} else {
+			} else {///just insert the hash key of minimizer
 				for (i = 0; i < s->n_seq; ++i)
 					for (j = 0; j < s->mz[i].n; ++j)
 						ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
@@ -640,6 +694,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 			}
 			free(s->mz);
 		}
+		///just clean seq
 		free(s->seq); free(s->len);
 		s->seq = 0, s->len = 0;
 		return s;
@@ -647,7 +702,9 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		st_data_t *s = (st_data_t*)in;
 		int i, n = 1<<p->opt->pre;
 		uint64_t n_ins = 0;
+		///for 0-th counting, p->pt = NULL
 		kt_for(p->opt->n_thread, worker_for_insert, s, n);
+		///n_ins is number of distinct k-mers
 		for (i = 0; i < n; ++i) {
 			n_ins += s->buf[i].n_ins;
 			if (p->pt) free(s->buf[i].b);
@@ -666,8 +723,51 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 	return 0;
 }

+void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs)
+{
+	int ret;
+    uint32_t i, m, pass, unpass;
+    gzFile fp = 0;
+    kseq_t *ks = NULL;
+	UC_Read ucr;
+	init_UC_Read(&ucr);
+
+    for (i = m = pass = unpass = 0; i < (uint32_t)asm_opt->num_reads; ++i)
+    {
+		if ((fp = gzopen(asm_opt->read_file_names[i], "r")) == 0) continue;
+		ks = kseq_init(fp);
+		while ((ret = kseq_read(ks)) >= 0) 
+		{
+			int l = ks->seq.l;
+			if((l - asm_opt->adapterLen*2) <= 0) continue;
+			recover_UC_Read(&ucr, rs, m);
+			fprintf(stderr, "l: %d, ucr.length: %lld, asm_opt->adapterLen: %d\n", 
+			l, ucr.length, asm_opt->adapterLen);
+			if(memcmp(ucr.seq, ks->seq.s+asm_opt->adapterLen, ucr.length) == 0)
+			{
+				pass++;
+			}
+			else
+			{
+				unpass++;
+			}
+			m++;
+		}
+        kseq_destroy(ks);
+        gzclose(fp);
+		ks = NULL;
+		fp = 0;
+    }
+
+	destory_UC_Read(&ucr);
+
+	fprintf(stderr, "[M::%s::# reads: %u, # pass: %u, # unpass: %u\n]", __func__, m, pass, unpass);
+	exit(1);
+}
+
 static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs, int64_t *n_seq)
 {
+	///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW
 	int read_rs = (rs && (flag & HAF_RS_READ));
 	pl_data_t pl;
 	gzFile fp = 0;
@@ -676,23 +776,29 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 	if (read_rs) {
 		pl.rs_in = rs;
 		init_UC_Read(&pl.ucr);
-	} else {
+	} else {///for 0-th counting, go into here
 		if ((fp = gzopen(fn, "r")) == 0) return 0;
 		pl.ks = kseq_init(fp);
 	}
+	///for 0-th counting, read all reads into pl.rs_out
 	if (rs && (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ)))
 		pl.rs_out = rs;
+	///for 0-th counting, flt_tab = NULL
+	///for 1-th counting, flt_tab = NULL
 	pl.flt_tab = flt_tab;
 	pl.opt = opt;
 	pl.flag = flag;
-	if (p0) {
+	if (p0) {///for 1-th counting, p0 = NULL
 		pl.pt = p0, pl.create_new = 0; // never create new elements in a position table
 		assert(p0->k == opt->k && p0->pre == opt->pre);
 	} else if (c0) {
 		pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW);
 		assert(c0->k == opt->k && c0->pre == opt->pre);
-	} else {
+	} else {///for ft-th counting and 1-th counting, go into here
 		pl.create_new = 1; // alware create new elements if the count table is empty
+		///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
+		///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
+		///building a large hash table consisting of 4096 small hash tables
 		pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift);
 	}
 	kt_pipeline(3, worker_count, &pl, 3);
@@ -713,6 +819,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
 	yak_copt_t opt;
 	ha_ct_t *h = 0;
 	assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
+	///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN
 	if (rs) {
 		if (flag & HAF_RS_WRITE_LEN)
 			init_All_reads(rs);
@@ -721,10 +828,16 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
 	}
 	yak_copt_init(&opt);
 	opt.k = asm_opt->k_mer_length;
+	///always 0
 	opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
+	///for ft-counting, shoud be 1
 	opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
+	///for ft-counting, shoud be 37
+	///for ha_pt_gen, shoud be 0
 	opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
 	opt.n_thread = asm_opt->thread_num;
+	opt.adaLen = asm_opt->adapterLen;
+	///asm_opt->num_reads is the number of fastq files
 	for (i = 0; i < asm_opt->num_reads; ++i)
 		h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
 	if (h && opt.bf_shift > 0)
@@ -771,23 +884,63 @@ void ha_ft_destroy(void *h)
 	if (h) yak_ft_destroy((yak_ft_t*)h);
 }

+
+void debug_ct_index(void* q_ct_idx, void* r_ct_idx)
+{
+	ha_ct_t* ct_idx = (ha_ct_t*)q_ct_idx;
+	yak_ct_t *g = NULL;
+	uint64_t i;
+	khint_t k;
+	for (i = 0; (int)i < 1<<ct_idx->pre; i++)
+	{
+		g = ct_idx->h[i].h;
+		for (k = 0; k < kh_end(g); ++k) 
+		{
+			if (kh_exist(g, k)) 
+			{
+				int c = kh_key(g, k) & YAK_MAX_COUNT;
+				uint64_t hash = ((kh_key(g, k) >> ct_idx->pre)<<ct_idx->pre) | i;
+				int q = query_ct_index(r_ct_idx, hash);
+				if(q!=c)
+				{
+					fprintf(stderr, "ERROR:c: %d, q: %d\n", c, q);
+				}
+			}
+		}
+	}
+}
+
 /*************************
 * High-level interfaces *
 *************************/

-void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode)
 {
 	yak_ft_t *flt_tab;
 	int64_t cnt[YAK_N_COUNTS];
-	int peak_hom, peak_het, cutoff;
+	int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
+	if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
 	ha_ct_t *h;
-	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs);
-	ha_ct_hist(h, cnt, asm_opt->thread_num);
-	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
-	if (hom_cov) *hom_cov = peak_hom;
-	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
-	cutoff = (int)(peak_hom * asm_opt->high_factor);
-	if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
+	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs);
+	if((asm_opt->flag & HA_F_VERBOSE_GFA))
+	{
+		write_ct_index((void*)h, asm_opt->output_file_name);
+		// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
+		// debug_ct_index((void*)h, ha_ct_table);
+		// debug_ct_index(ha_ct_table, (void*)h);
+		// ha_ct_destroy((ha_ct_t *)ha_ct_table);
+	}
+	
+	if(!(ex_flag & HAF_SKIP_READ))
+	{
+		ha_ct_hist(h, cnt, asm_opt->thread_num);
+		peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
+		if (hom_cov) *hom_cov = peak_hom;
+		if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+		///in default, asm_opt->high_factor = 5.0
+		cutoff = (int)(peak_hom * asm_opt->high_factor);
+		if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
+	}
 	ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
 	flt_tab = gen_hh(h);
 	ha_ct_destroy(h);
@@ -796,36 +949,42 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
 	return (void*)flt_tab;
 }

-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov)
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov)
 {
 	int64_t cnt[YAK_N_COUNTS], tot_cnt;
 	int peak_hom, peak_het, i, extra_flag1, extra_flag2;
 	ha_ct_t *ct;
 	ha_pt_t *pt;
-	if (read_from_store) {
+	if (read_from_store) {///if reads have already been read
 		extra_flag1 = extra_flag2 = HAF_RS_READ;
-	} else if (rs->total_reads == 0) {
+	} else if (rs->total_reads == 0) {///if reads & length have not been scanned
 		extra_flag1 = HAF_RS_WRITE_LEN;
 		extra_flag2 = HAF_RS_WRITE_SEQ;
-	} else {
+	} else {///if length has been loaded but reads have not
 		extra_flag1 = HAF_RS_WRITE_SEQ;
 		extra_flag2 = HAF_RS_READ;
 	}
+	if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
+
 	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
 			yak_realtime(), yak_cpu_usage(), (long)ct->tot);
 	ha_ct_hist(ct, cnt, asm_opt->thread_num);
 	fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
-	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
+	peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
 	if (hom_cov) *hom_cov = peak_hom;
 	if (het_cov) *het_cov = peak_het;
 	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+	///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time
 	if (flt_tab == 0) {
 		int cutoff = (int)(peak_hom * asm_opt->high_factor);
 		if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
+		if((extra_flag1 & HAF_SKIP_READ) && (extra_flag2 & HAF_SKIP_READ)) cutoff = YAK_MAX_COUNT - 1;
 		ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num);
 		for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i;
 	} else {
+		///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times
+		///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen
 		ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
 		for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;
 	}
@@ -837,3 +996,239 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 			yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
 	return pt;
 }
+
+int query_ct_index(void* ct_idx, uint64_t hash)
+{
+	ha_ct1_t *g = &(((ha_ct_t*)ct_idx)->h[hash & ((1ULL<<((ha_ct_t*)ct_idx)->pre) - 1)]);
+	khint_t k;
+	k = yak_ct_get(g->h, hash);
+	if (k == kh_end(g->h)) return 0;
+	return kh_key(g->h, k)&YAK_MAX_COUNT;
+}
+
+
+int write_ct_index(void *i_ct_idx, char* file_name)
+{
+	char* gfa_name = (char*)malloc(strlen(file_name)+25);
+    sprintf(gfa_name, "%s.ct_flt", file_name);
+    FILE* fp = fopen(gfa_name, "w");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+	ha_ct_t* ct_idx = (ha_ct_t*)i_ct_idx;
+	int i;
+	ha_ct1_t *g;
+	fwrite(&ct_idx->k, sizeof(ct_idx->k), 1, fp);
+	fwrite(&ct_idx->pre, sizeof(ct_idx->pre), 1, fp);
+	fwrite(&ct_idx->n_hash, sizeof(ct_idx->n_hash), 1, fp);
+	fwrite(&ct_idx->n_shift, sizeof(ct_idx->n_shift), 1, fp);
+	fwrite(&ct_idx->tot, sizeof(ct_idx->tot), 1, fp);
+	for (i = 0; i < 1<<ct_idx->pre; i++)
+	{
+		g = &(ct_idx->h[i]);
+		yak_ct_save(g->h, fp);
+	}
+
+
+	fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
+    free(gfa_name);
+	fclose(fp);
+	return 1;
+}
+
+int load_ct_index(void **i_ct_idx, char* file_name)
+{
+	char* gfa_name = (char*)malloc(strlen(file_name)+25);
+    sprintf(gfa_name, "%s.ct_flt", file_name);
+    FILE* fp = fopen(gfa_name, "r");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+	ha_ct_t** ct_idx = (ha_ct_t**)i_ct_idx;
+	double index_time = 0;
+	uint64_t flag = 0;
+	int i;
+	ha_ct_t *h = 0;
+	ha_ct1_t *g;
+	CALLOC(h, 1);
+
+	flag += fread(&h->k, sizeof(h->k), 1, fp);
+	flag += fread(&h->pre, sizeof(h->pre), 1, fp);
+	flag += fread(&h->n_hash, sizeof(h->n_hash), 1, fp);
+	flag += fread(&h->n_shift, sizeof(h->n_shift), 1, fp);
+	flag += fread(&h->tot, sizeof(h->tot), 1, fp);
+	CALLOC(h->h, 1<<h->pre);
+
+
+	index_time = yak_realtime();
+	for (i = 0; i < 1<<h->pre; ++i) 
+	{
+		g = &(h->h[i]);
+		yak_ct_load(&(g->h), fp);
+	}
+
+	(*ct_idx) = h;
+	fprintf(stderr, "[M::%s::%.3f] ==> Loaded count table\n", __func__, yak_realtime() - index_time);
+	fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
+	free(gfa_name);
+	return 1;
+}
+
+int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
+{
+    char* gfa_name = (char*)malloc(strlen(file_name)+25);
+    sprintf(gfa_name, "%s.pt_flt", file_name);
+    FILE* fp = fopen(gfa_name, "w");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+	yak_ft_t *ha_flt_tab = (yak_ft_t*)flt_tab;
+
+	if(ha_flt_tab)
+	{
+		fwrite("f", 1, 1, fp);
+		yak_ft_save(ha_flt_tab, fp);
+	}
+
+
+	if(ha_idx)
+	{
+		int i;
+		ha_pt1_t *g;
+		fwrite("h", 1, 1, fp);
+		fwrite(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
+		fwrite(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
+		fwrite(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
+		fwrite(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
+
+		for (i = 0; i < 1<<ha_idx->pre; ++i) 
+		{
+			g = &(ha_idx->h[i]);
+			yak_pt_save(g->h, fp);
+			fwrite(&g->n, sizeof(g->n), 1, fp);
+			fwrite(g->a, sizeof(ha_idxpos_t), g->n, fp);
+		}
+	}
+
+	fwrite(&opt->number_of_round, sizeof(opt->number_of_round), 1, fp);
+	fwrite(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
+	fwrite(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
+	fwrite(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
+
+
+	write_All_reads(r, gfa_name);
+
+	fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
+    free(gfa_name);
+	fclose(fp);
+	return 1;
+}
+
+int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
+{
+	char* gfa_name = (char*)malloc(strlen(file_name)+25);
+    sprintf(gfa_name, "%s.pt_flt", file_name);
+    FILE* fp = fopen(gfa_name, "r");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+
+	ha_pt_t *ha_idx = NULL;
+	char mode = 0;
+	int f_flag, absent, i;
+	double index_time, index_s_time, pos_time, pos_s_time;
+
+	
+	
+	f_flag += fread(&mode, 1, 1, fp);
+	if(mode == 'f')
+	{
+		index_time = yak_realtime();
+
+		yak_ft_load((yak_ft_t **)r_flt_tab, fp);
+		
+	    f_flag += fread(&mode, 1, 1, fp);
+
+		fprintf(stderr, "[M::%s::%.3f] ==> Loaded flt table\n", __func__, yak_realtime()-index_time);
+	}
+	///insert using multiple threads???
+	if(mode == 'h')
+	{
+		pos_time = index_time = 0;
+
+		CALLOC(ha_idx, 1);
+		ha_pt1_t *g;
+		f_flag += fread(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
+		f_flag += fread(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
+		f_flag += fread(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
+		f_flag += fread(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
+		CALLOC(ha_idx->h, 1<<ha_idx->pre);
+		for (i = 0; i < 1<<ha_idx->pre; ++i) 
+		{
+			index_s_time = yak_realtime();
+
+			g = &(ha_idx->h[i]);
+			yak_pt_load(&(g->h), fp);
+
+			index_time += yak_realtime() - index_s_time;
+			
+			pos_s_time = yak_realtime();
+
+			f_flag += fread(&g->n, sizeof(g->n), 1, fp);
+			MALLOC(g->a, g->n);
+			f_flag += fread(g->a, sizeof(ha_idxpos_t), g->n, fp);
+
+			pos_time += yak_realtime() - pos_s_time;
+		}
+		(*r_ha_idx) = ha_idx;
+
+		fprintf(stderr, "[M::%s::%.3f(index)/%.3f(pos)] ==> Loaded pos table\n", __func__, index_time, pos_time);
+	}
+
+	if(mode != 'h' && mode != 'f')
+	{
+		free(gfa_name);
+		fclose(fp);
+		return 0;
+	}
+
+
+	f_flag += fread(&absent, sizeof(absent), 1, fp);
+	if(absent != opt->number_of_round)
+	{
+		fprintf(stderr, "ERROR: different number of rounds!\n");
+		exit(1);
+	}
+
+	f_flag += fread(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
+	f_flag += fread(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
+	f_flag += fread(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
+
+
+	fclose(fp);
+	
+	if(!load_All_reads(r, gfa_name))
+	{
+		free(gfa_name);
+		return 0;
+	}
+
+
+	memset(r->trio_flag, AMBIGU, r->total_reads*sizeof(uint8_t));
+	r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
+    r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
+	for (i = 0; i < (long long)r->total_reads; i++)
+    {
+        init_ma_hit_t_alloc(&(r->paf[i]));
+        init_ma_hit_t_alloc(&(r->reverse_paf[i]));
+    }
+
+	fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
+
+	free(gfa_name);
+	return 1;
+}
@@ -6,7 +6,9 @@
 #include "CommandLines.h"

 typedef struct {
-	uint64_t x;
+	uint64_t x; ///x is the hash key
+	///rid is the read id, pos is the end pos of this minimizer, rev is the direction
+	///span is the length of this k-mer. For non-HPC k-mer, span may not be equal to k
 	uint64_t rid:28, pos:27, rev:1, span:8;
 } ha_mz1_t;

@@ -25,15 +27,25 @@ typedef struct ha_abuf_s ha_abuf_t;
 extern const unsigned char seq_nt4_table[256];
 extern void *ha_flt_tab;
 extern ha_pt_t *ha_idx;
+extern void *ha_flt_tab_hp;
+extern ha_pt_t *ha_idx_hp;
+extern void *ha_ct_table;

-void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov);
+
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
 int ha_ft_isflt(const void *hh, uint64_t y);
 void ha_ft_destroy(void *h);

-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov);
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov);
 void ha_pt_destroy(ha_pt_t *h);
 const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);

+int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
+int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
+int write_ct_index(void *ct_idx, char* file_name);
+int load_ct_index(void **ct_idx, char* file_name);
+int query_ct_index(void* ct_idx, uint64_t hash);
+
 ha_abuf_t *ha_abuf_init(void);
 void ha_abuf_destroy(ha_abuf_t *ab);
 uint64_t ha_abuf_mem(const ha_abuf_t *ab);
@@ -48,7 +60,9 @@ 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 ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
+void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct);
+int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het);
+void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);

 static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
 {
@@ -76,6 +90,7 @@ static inline uint64_t yak_hash64_64(uint64_t key)

 static inline uint64_t yak_hash_long(uint64_t x[4])
 {
+	///compare forward k-mer and reverse complementary strand
 	int j = x[1] < x[3]? 0 : 1;
 	return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
 }
@@ -133,6 +133,30 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 			h->count = 0; \
 		} \
 	}
+#define __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { \
+		if (!h) return 0; \
+		khint_t n_buckets = (h->keys? 1U<<h->bits : 0U); \
+		fwrite(&n_buckets, sizeof(n_buckets), 1, fp); \
+		fwrite(&h->bits, sizeof(h->bits), 1, fp); \
+		fwrite(&h->count, sizeof(h->count), 1, fp); \
+		fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \
+		return 1; \
+	} \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \
+		(*h) = prefix##_init(); \
+		khint_t n_buckets; \
+		uint64_t flag = 0;\
+		flag += fread(&n_buckets, sizeof(n_buckets), 1, fp); \
+		flag += fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \
+		flag += fread(&(*h)->count, sizeof((*h)->count), 1, fp); \
+		(*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \
+		(*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \
+		flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
+		return 1; \
+	} \

 #define __KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	SCOPE khint_t prefix##_getp(const HType *h, const khkey_t *key) { \
@@ -245,6 +269,7 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 #define KHASHL_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	__KHASHL_TYPE(HType, khkey_t) \
 	__KHASHL_IMPL_BASIC(SCOPE, HType, prefix) \
+	__KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
 	__KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	__KHASHL_IMPL_RESIZE(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	__KHASHL_IMPL_PUT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
@@ -264,6 +289,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_s, HType##_s_bucket_t, prefix##_s_hash, prefix##_s_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_s_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_s_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_s_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_s_load(h, fp); } \
 	SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_s_resize(h, new_n_buckets); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_s_bucket_t t; t.key = key; return prefix##_s_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_s_del(h, k); } \
@@ -276,6 +303,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_m_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_m_load(h, fp); } \
 	SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_m_resize(h, new_n_buckets); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \
@@ -287,6 +316,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_cs, HType##_cs_bucket_t, __kh_cached_hash, prefix##_cs_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_cs_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_cs_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cs_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cs_load(h, fp); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cs_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cs_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cs_del(h, k); } \
 	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cs_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cs_putp(h, &t, absent); }
@@ -297,6 +328,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_cm, HType##_cm_bucket_t, __kh_cached_hash, prefix##_cm_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_cm_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_cm_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cm_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cm_load(h, fp); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cm_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cm_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cm_del(h, k); } \
 	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cm_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cm_putp(h, &t, absent); }
@@ -0,0 +1,177 @@
+#ifndef KSW2_H_
+#define KSW2_H_
+
+#include <stdint.h>
+
+#define KSW_NEG_INF -0x40000000
+
+#define KSW_EZ_SCORE_ONLY  0x01 // don't record alignment path/cigar
+#define KSW_EZ_RIGHT       0x02 // right-align gaps
+#define KSW_EZ_GENERIC_SC  0x04 // without this flag: match/mismatch only; last symbol is a wildcard
+#define KSW_EZ_APPROX_MAX  0x08 // approximate max; this is faster with sse
+#define KSW_EZ_APPROX_DROP 0x10 // approximate Z-drop; faster with sse
+#define KSW_EZ_EXTZ_ONLY   0x40 // only perform extension
+#define KSW_EZ_REV_CIGAR   0x80 // reverse CIGAR in the output
+#define KSW_EZ_SPLICE_FOR  0x100
+#define KSW_EZ_SPLICE_REV  0x200
+#define KSW_EZ_SPLICE_FLANK 0x400
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct {
+	uint32_t max:31, zdropped:1;
+	int max_q, max_t;      // max extension coordinate
+	int mqe, mqe_t;        // max score when reaching the end of query
+	int mte, mte_q;        // max score when reaching the end of target
+	int score;             // max score reaching both ends; may be KSW_NEG_INF
+	int m_cigar, n_cigar;
+	int reach_end;
+	uint32_t *cigar;
+} ksw_extz_t;
+
+/**
+ * NW-like extension
+ *
+ * @param km        memory pool, when used with kalloc
+ * @param qlen      query length
+ * @param query     query sequence with 0 <= query[i] < m
+ * @param tlen      target length
+ * @param target    target sequence with 0 <= target[i] < m
+ * @param m         number of residue types
+ * @param mat       m*m scoring mattrix in one-dimension array
+ * @param gapo      gap open penalty; a gap of length l cost "-(gapo+l*gape)"
+ * @param gape      gap extension penalty
+ * @param w         band width (<0 to disable)
+ * @param zdrop     off-diagonal drop-off to stop extension (positive; <0 to disable)
+ * @param flag      flag (see KSW_EZ_* macros)
+ * @param ez        (out) scores and cigar
+ */
+void ksw_extz(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+			  int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez);
+
+void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+
+void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+			  int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
+
+void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+
+void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
+
+void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez);
+
+/**
+ * Global alignment
+ *
+ * (first 10 parameters identical to ksw_extz_sse())
+ * @param m_cigar   (modified) max CIGAR length; feed 0 if cigar==0
+ * @param n_cigar   (out) number of CIGAR elements
+ * @param cigar     (out) BAM-encoded CIGAR; caller need to deallocate with kfree(km, )
+ *
+ * @return          score of the alignment
+ */
+int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
+int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
+int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
+
+void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat);
+int ksw_ll_i16(void *q, int tlen, const uint8_t *target, int gapo, int gape, int *qe, int *te);
+
+#ifdef __cplusplus
+}
+#endif
+
+/************************************
+ *** Private macros and functions ***
+ ************************************/
+
+#ifdef HAVE_KALLOC
+#include "kalloc.h"
+#else
+#include <stdlib.h>
+#define kmalloc(km, size) malloc((size))
+#define kcalloc(km, count, size) calloc((count), (size))
+#define krealloc(km, ptr, size) realloc((ptr), (size))
+#define kfree(km, ptr) free((ptr))
+#endif
+
+static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uint32_t *cigar, uint32_t op, int len)
+{
+	if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) {
+		if (*n_cigar == *m_cigar) {
+			*m_cigar = *m_cigar? (*m_cigar)<<1 : 4;
+			cigar = (uint32_t*)krealloc(km, cigar, (*m_cigar) << 2);
+		}
+		cigar[(*n_cigar)++] = len<<4 | op;
+	} else cigar[(*n_cigar)-1] += len<<4;
+	return cigar;
+}
+
+// In the backtrack matrix, value p[] has the following structure:
+//   bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F}
+//   bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E})
+//   bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F})
+static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intron_len, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0,
+								 int *m_cigar_, int *n_cigar_, uint32_t **cigar_)
+{ // p[] - lower 3 bits: which type gets the max; bit
+	int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0;
+	uint32_t *cigar = *cigar_, tmp;
+	while (i >= 0 && j >= 0) { // at the beginning of the loop, _state_ tells us which state to check
+		int force_state = -1;
+		if (is_rot) {
+			r = i + j;
+			if (i < off[r]) force_state = 2;
+			if (off_end && i > off_end[r]) force_state = 1;
+			tmp = force_state < 0? p[(size_t)r * n_col + i - off[r]] : 0;
+		} else {
+			if (j < off[i]) force_state = 2;
+			if (off_end && j > off_end[i]) force_state = 1;
+			tmp = force_state < 0? p[(size_t)i * n_col + j - off[i]] : 0;
+		}
+		if (state == 0) state = tmp & 7; // if requesting the H state, find state one maximizes it.
+		else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H
+		if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
+		if (force_state >= 0) state = force_state;
+		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
+		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
+		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
+		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
+	}
+	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
+	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
+	if (!is_rev)
+		for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
+			tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp;
+	*m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar;
+}
+
+static inline void ksw_reset_extz(ksw_extz_t *ez)
+{
+	ez->max_q = ez->max_t = ez->mqe_t = ez->mte_q = -1;
+	ez->max = 0, ez->score = ez->mqe = ez->mte = KSW_NEG_INF;
+	ez->n_cigar = 0, ez->zdropped = 0, ez->reach_end = 0;
+}
+
+static inline int ksw_apply_zdrop(ksw_extz_t *ez, int is_rot, int32_t H, int a, int b, int zdrop, int8_t e)
+{
+	int r, t;
+	if (is_rot) r = a, t = b;
+	else r = a + b, t = a;
+	if (H > (int32_t)ez->max) {
+		ez->max = H, ez->max_t = t, ez->max_q = r - t;
+	} else if (t >= ez->max_t && r - t >= ez->max_q) {
+		int tl = t - ez->max_t, ql = (r - t) - ez->max_q, l;
+		l = tl > ql? tl - ql : ql - tl;
+		if (zdrop >= 0 && ez->max - H > zdrop + l * e) {
+			ez->zdropped = 1;
+			return 1;
+		}
+	}
+	return 0;
+}
+#endif
@@ -0,0 +1,305 @@
+#include <string.h>
+#include <assert.h>
+#include "ksw2.h"
+
+#ifdef __SSE2__
+#include <emmintrin.h>
+
+#ifdef KSW_SSE2_ONLY
+#undef __SSE4_1__
+#endif
+
+#ifdef __SSE4_1__
+#include <smmintrin.h>
+#endif
+
+#ifdef KSW_CPU_DISPATCH
+#ifdef __SSE4_1__
+void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+#else
+void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+#endif
+#else
+void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+#endif // ~KSW_CPU_DISPATCH
+{
+#define __dp_code_block1 \
+	z = _mm_add_epi8(_mm_load_si128(&s[t]), qe2_); \
+	xt1 = _mm_load_si128(&x[t]);                     /* xt1 <- x[r-1][t..t+15] */ \
+	tmp = _mm_srli_si128(xt1, 15);                   /* tmp <- x[r-1][t+15] */ \
+	xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_); /* xt1 <- x[r-1][t-1..t+14] */ \
+	x1_ = tmp; \
+	vt1 = _mm_load_si128(&v[t]);                     /* vt1 <- v[r-1][t..t+15] */ \
+	tmp = _mm_srli_si128(vt1, 15);                   /* tmp <- v[r-1][t+15] */ \
+	vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); /* vt1 <- v[r-1][t-1..t+14] */ \
+	v1_ = tmp; \
+	a = _mm_add_epi8(xt1, vt1);                      /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \
+	ut = _mm_load_si128(&u[t]);                      /* ut <- u[t..t+15] */ \
+	b = _mm_add_epi8(_mm_load_si128(&y[t]), ut);     /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */
+
+#define __dp_code_block2 \
+	z = _mm_max_epu8(z, b);                          /* z = max(z, b); this works because both are non-negative */ \
+	z = _mm_min_epu8(z, max_sc_); \
+	_mm_store_si128(&u[t], _mm_sub_epi8(z, vt1));    /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \
+	_mm_store_si128(&v[t], _mm_sub_epi8(z, ut));     /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \
+	z = _mm_sub_epi8(z, q_); \
+	a = _mm_sub_epi8(a, z); \
+	b = _mm_sub_epi8(b, z);
+
+	int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc;
+	int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
+	int32_t *H = 0, H0 = 0, last_H0_t = 0;
+	uint8_t *qr, *sf, *mem, *mem2 = 0;
+	__m128i q_, qe2_, zero_, flag1_, flag2_, flag8_, flag16_, sc_mch_, sc_mis_, sc_N_, m1_, max_sc_;
+	__m128i *u, *v, *x, *y, *s, *p = 0;
+
+	ksw_reset_extz(ez);
+	if (m <= 0 || qlen <= 0 || tlen <= 0) return;
+
+	zero_   = _mm_set1_epi8(0);
+	q_      = _mm_set1_epi8(q);
+	qe2_    = _mm_set1_epi8((q + e) * 2);
+	flag1_  = _mm_set1_epi8(1);
+	flag2_  = _mm_set1_epi8(2);
+	flag8_  = _mm_set1_epi8(0x08);
+	flag16_ = _mm_set1_epi8(0x10);
+	sc_mch_ = _mm_set1_epi8(mat[0]);
+	sc_mis_ = _mm_set1_epi8(mat[1]);
+	sc_N_   = mat[m*m-1] == 0? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m*m-1]);
+	m1_     = _mm_set1_epi8(m - 1); // wildcard
+	max_sc_ = _mm_set1_epi8(mat[0] + (q + e) * 2);
+
+	if (w < 0) w = tlen > qlen? tlen : qlen;
+	wl = wr = w;
+	tlen_ = (tlen + 15) / 16;
+	n_col_ = qlen < tlen? qlen : tlen;
+	n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 15) / 16 + 1;
+	qlen_ = (qlen + 15) / 16;
+	for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) {
+		max_sc = max_sc > mat[t]? max_sc : mat[t];
+		min_sc = min_sc < mat[t]? min_sc : mat[t];
+	}
+	if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches
+
+	mem = (uint8_t*)kcalloc(km, tlen_ * 6 + qlen_ + 1, 16);
+	u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned
+	v = u + tlen_, x = v + tlen_, y = x + tlen_, s = y + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
+	if (!approx_max) {
+		H = (int32_t*)kmalloc(km, tlen_ * 16 * 4);
+		for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
+	}
+	if (with_cigar) {
+		mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
+		p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
+		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
+		off_end = off + qlen + tlen - 1;
+	}
+
+	for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t];
+	memcpy(sf, target, tlen);
+
+	for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) {
+		int st = 0, en = tlen - 1, st0, en0, st_, en_;
+		int8_t x1, v1;
+		uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v;
+		__m128i x1_, v1_;
+		// find the boundaries
+		if (st < r - qlen + 1) st = r - qlen + 1;
+		if (en > r) en = r;
+		if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil
+		if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor
+		if (st > en) {
+			ez->zdropped = 1;
+			break;
+		}
+		st0 = st, en0 = en;
+		st = st / 16 * 16, en = (en + 16) / 16 * 16 - 1;
+		// set boundary conditions
+		if (st > 0) {
+			if (st - 1 >= last_st && st - 1 <= last_en)
+				x1 = ((uint8_t*)x)[st - 1], v1 = v8[st - 1]; // (r-1,s-1) calculated in the last round
+			else x1 = v1 = 0; // not calculated; set to zeros
+		} else x1 = 0, v1 = r? q : 0;
+		if (en >= r) ((uint8_t*)y)[r] = 0, u8[r] = r? q : 0;
+		// loop fission: set scores first
+		if (!(flag & KSW_EZ_GENERIC_SC)) {
+			for (t = st0; t <= en0; t += 16) {
+				__m128i sq, st, tmp, mask;
+				sq = _mm_loadu_si128((__m128i*)&sf[t]);
+				st = _mm_loadu_si128((__m128i*)&qrr[t]);
+				mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_));
+				tmp = _mm_cmpeq_epi8(sq, st);
+#ifdef __SSE4_1__
+				tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp);
+				tmp = _mm_blendv_epi8(tmp,     sc_N_,   mask);
+#else
+				tmp = _mm_or_si128(_mm_andnot_si128(tmp,  sc_mis_), _mm_and_si128(tmp,  sc_mch_));
+				tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp),     _mm_and_si128(mask, sc_N_));
+#endif
+				_mm_storeu_si128((__m128i*)((uint8_t*)s + t), tmp);
+			}
+		} else {
+			for (t = st0; t <= en0; ++t)
+				((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]];
+		}
+		// core loop
+		x1_ = _mm_cvtsi32_si128(x1);
+		v1_ = _mm_cvtsi32_si128(v1);
+		st_ = st / 16, en_ = en / 16;
+		assert(en_ - st_ + 1 <= n_col_);
+		if (!with_cigar) { // score only
+			for (t = st_; t <= en_; ++t) {
+				__m128i z, a, b, xt1, vt1, ut, tmp;
+				__dp_code_block1;
+#ifdef __SSE4_1__
+				z = _mm_max_epi8(z, a);                          // z = z > a? z : a (signed)
+#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8()
+				z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_));  // z = z > 0? z : 0;
+				z = _mm_max_epu8(z, a);                          // z = max(z, a); this works because both are non-negative
+#endif
+				__dp_code_block2;
+#ifdef __SSE4_1__
+				_mm_store_si128(&x[t], _mm_max_epi8(a, zero_));
+				_mm_store_si128(&y[t], _mm_max_epi8(b, zero_));
+#else
+				tmp = _mm_cmpgt_epi8(a, zero_);
+				_mm_store_si128(&x[t], _mm_and_si128(a, tmp));
+				tmp = _mm_cmpgt_epi8(b, zero_);
+				_mm_store_si128(&y[t], _mm_and_si128(b, tmp));
+#endif
+			}
+		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
+			off[r] = st, off_end[r] = en;
+			for (t = st_; t <= en_; ++t) {
+				__m128i d, z, a, b, xt1, vt1, ut, tmp;
+				__dp_code_block1;
+				d = _mm_and_si128(_mm_cmpgt_epi8(a, z), flag1_); // d = a > z? 1 : 0
+#ifdef __SSE4_1__
+				z = _mm_max_epi8(z, a);                          // z = z > a? z : a (signed)
+				tmp = _mm_cmpgt_epi8(b, z);
+				d = _mm_blendv_epi8(d, flag2_, tmp);             // d = b > z? 2 : d
+#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
+				z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_));  // z = z > 0? z : 0;
+				z = _mm_max_epu8(z, a);                          // z = max(z, a); this works because both are non-negative
+				tmp = _mm_cmpgt_epi8(b, z);
+				d = _mm_or_si128(_mm_andnot_si128(tmp, d), _mm_and_si128(tmp, flag2_)); // d = b > z? 2 : d; emulating blendv
+#endif
+				__dp_code_block2;
+				tmp = _mm_cmpgt_epi8(a, zero_);
+				_mm_store_si128(&x[t], _mm_and_si128(tmp, a));
+				d = _mm_or_si128(d, _mm_and_si128(tmp, flag8_));  // d = a > 0? 0x08 : 0
+				tmp = _mm_cmpgt_epi8(b, zero_);
+				_mm_store_si128(&y[t], _mm_and_si128(tmp, b));
+				d = _mm_or_si128(d, _mm_and_si128(tmp, flag16_)); // d = b > 0? 0x10 : 0
+				_mm_store_si128(&pr[t], d);
+			}
+		} else { // gap right-alignment
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
+			off[r] = st, off_end[r] = en;
+			for (t = st_; t <= en_; ++t) {
+				__m128i d, z, a, b, xt1, vt1, ut, tmp;
+				__dp_code_block1;
+				d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), flag1_); // d = z > a? 0 : 1
+#ifdef __SSE4_1__
+				z = _mm_max_epi8(z, a);                          // z = z > a? z : a (signed)
+				tmp = _mm_cmpgt_epi8(z, b);
+				d = _mm_blendv_epi8(flag2_, d, tmp);             // d = z > b? d : 2
+#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
+				z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_));  // z = z > 0? z : 0;
+				z = _mm_max_epu8(z, a);                          // z = max(z, a); this works because both are non-negative
+				tmp = _mm_cmpgt_epi8(z, b);
+				d = _mm_or_si128(_mm_andnot_si128(tmp, flag2_), _mm_and_si128(tmp, d)); // d = z > b? d : 2; emulating blendv
+#endif
+				__dp_code_block2;
+				tmp = _mm_cmpgt_epi8(zero_, a);
+				_mm_store_si128(&x[t], _mm_andnot_si128(tmp, a));
+				d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag8_));  // d = 0 > a? 0 : 0x08
+				tmp = _mm_cmpgt_epi8(zero_, b);
+				_mm_store_si128(&y[t], _mm_andnot_si128(tmp, b));
+				d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag16_)); // d = 0 > b? 0 : 0x10
+				_mm_store_si128(&pr[t], d);
+			}
+		}
+		if (!approx_max) { // find the exact max with a 32-bit score array
+			int32_t max_H, max_t;
+			// compute H[], max_H and max_t
+			if (r > 0) {
+				int32_t HH[4], tt[4], en1 = st0 + (en0 - st0) / 4 * 4, i;
+				__m128i max_H_, max_t_, qe_;
+				max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] - qe : H[en0] + v8[en0] - qe; // special casing the last element
+				max_t = en0;
+				max_H_ = _mm_set1_epi32(max_H);
+				max_t_ = _mm_set1_epi32(max_t);
+				qe_    = _mm_set1_epi32(q + e);
+				for (t = st0; t < en1; t += 4) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t;
+					__m128i H1, tmp, t_;
+					H1 = _mm_loadu_si128((__m128i*)&H[t]);
+					t_ = _mm_setr_epi32(v8[t], v8[t+1], v8[t+2], v8[t+3]);
+					H1 = _mm_add_epi32(H1, t_);
+					H1 = _mm_sub_epi32(H1, qe_);
+					_mm_storeu_si128((__m128i*)&H[t], H1);
+					t_ = _mm_set1_epi32(t);
+					tmp = _mm_cmpgt_epi32(H1, max_H_);
+#ifdef __SSE4_1__
+					max_H_ = _mm_blendv_epi8(max_H_, H1, tmp);
+					max_t_ = _mm_blendv_epi8(max_t_, t_, tmp);
+#else
+					max_H_ = _mm_or_si128(_mm_and_si128(tmp, H1), _mm_andnot_si128(tmp, max_H_));
+					max_t_ = _mm_or_si128(_mm_and_si128(tmp, t_), _mm_andnot_si128(tmp, max_t_));
+#endif
+				}
+				_mm_storeu_si128((__m128i*)HH, max_H_);
+				_mm_storeu_si128((__m128i*)tt, max_t_);
+				for (i = 0; i < 4; ++i)
+					if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i;
+				for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE
+					H[t] += (int32_t)v8[t] - qe;
+					if (H[t] > max_H)
+						max_H = H[t], max_t = t;
+				}
+			} else H[0] = v8[0] - qe - qe, max_H = H[0], max_t = 0; // special casing r==0
+			// update ez
+			if (en0 == tlen - 1 && H[en0] > ez->mte)
+				ez->mte = H[en0], ez->mte_q = r - en;
+			if (r - st0 == qlen - 1 && H[st0] > ez->mqe)
+				ez->mqe = H[st0], ez->mqe_t = st0;
+			if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e)) break;
+			if (r == qlen + tlen - 2 && en0 == tlen - 1)
+				ez->score = H[tlen - 1];
+		} else { // find approximate max; Z-drop might be inaccurate, too.
+			if (r > 0) {
+				if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) {
+					int32_t d0 = v8[last_H0_t] - qe;
+					int32_t d1 = u8[last_H0_t + 1] - qe;
+					if (d0 > d1) H0 += d0;
+					else H0 += d1, ++last_H0_t;
+				} else if (last_H0_t >= st0 && last_H0_t <= en0) {
+					H0 += v8[last_H0_t] - qe;
+				} else {
+					++last_H0_t, H0 += u8[last_H0_t] - qe;
+				}
+				if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e)) break;
+			} else H0 = v8[0] - qe - qe, last_H0_t = 0;
+			if (r == qlen + tlen - 2 && en0 == tlen - 1)
+				ez->score = H0;
+		}
+		last_st = st, last_en = en;
+		//for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging
+	}
+	kfree(km, mem);
+	if (!approx_max) kfree(km, H);
+	if (with_cigar) { // backtrack
+		int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
+		if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) {
+			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+		} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) {
+			ez->reach_end = 1;
+			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+		} else if (ez->max_t >= 0 && ez->max_q >= 0) {
+			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+		}
+		kfree(km, mem2); kfree(km, off);
+	}
+}
+#endif // __SSE2__
@@ -37,7 +37,11 @@ static inline int tq_shift(tiny_queue_t *q)
 * @param p      minimizers
 */
 void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf)
-{
+{	///in default, w = 51, k = 51, is_hpc = 1
+	/**
+	 uint64_t x;
+	 uint64_t rid:28, pos:27, rev:1, span:8;
+	 **/
 	static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
 	uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
 	int i, j, l, buf_pos, min_pos, kmer_span = 0;
@@ -45,8 +49,10 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 	tiny_queue_t tq;

 	assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
+	///sizeof(ha_mz1_t) = 16
 	memset(buf, 0xff, w * 16);
 	memset(&tq, 0, sizeof(tiny_queue_t));
+	///len/w is the evaluated minimizer numbers
 	kv_resize(ha_mz1_t, *p, p->n + len/w);

 	for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
@@ -65,7 +71,11 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 				tq_push(&tq, skip_len);
 				kmer_span += skip_len;
 				if (tq.count > k) kmer_span -= tq_shift(&tq);
-			} else kmer_span = l + 1 < k? l + 1 : k;
+			} else kmer_span = l + 1 < k? l + 1 : k; 
+			///kmer_span should be used for HPC k-mer
+			///so for non-HPC k-mer, kmer_span should be k in any case?
+			///kmer_span is used to calculate anchor pos on reverse complementary strand
+
 			kmer[0] = (kmer[0] << 1 | (c&1))  & mask;                  // forward k-mer
 			kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
 			kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1))  << shift1; // reverse k-mer
@@ -80,6 +90,12 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 					info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
 			}
 		} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
+
+
+		//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
+		//i is the real base iterator, while l is the HPC base iterator
+		//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
+		//but even if l < k, infor is still stored into buf
 		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
 		if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
 			for (j = buf_pos + 1; j < w; ++j)
@@ -87,15 +103,26 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 			for (j = 0; j < buf_pos; ++j)
 				if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
 		}
+		/**
+		 * There are three cases:
+		 * 1. info.x <= min.x, means info is a new minimizer
+		 * 2. info.x > min.x, info is not a new minimizer
+		 *    (1) buf_pos != min_pos, do nothing
+		 *    (2) buf_pos == min_pos, means current minimizer has moved outside the window
+		 * **/
+		///three cases: 1. 
 		if (info.x <= min.x) { // a new minimum; then write the old min
 			if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
 			min = info, min_pos = buf_pos;
 		} else if (buf_pos == min_pos) { // old min has moved outside the window
 			if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			///buf_pos == min_pos, means current minimizer has moved outside the window
+			///so for now we need to find a new minimizer at the current window (w k-mers)
 			for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
 				if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
 			for (j = 0; j <= buf_pos; ++j)
 				if (min.x >= buf[j].x) min = buf[j], min_pos = j;
+
 			if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
 				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
 					if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
@@ -108,3 +135,127 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 	if (min.x != UINT64_MAX)
 		kv_push(ha_mz1_t, *p, min);
 }
+
+
+
+void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, 
+kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
+{	///in default, w = 51, k = 51, is_hpc = 1
+	/**
+	 uint64_t x;
+	 uint64_t rid:28, pos:27, rev:1, span:8;
+	 **/
+	extern void *ha_ct_table;
+	if(dbg_ct != NULL) dbg_ct->a.n = 0;
+
+	static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
+	uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0}, filtered;
+	int i, j, l, buf_pos, min_pos, kmer_span = 0;
+	ha_mz1_t buf[256], min = dummy;
+	tiny_queue_t tq;
+	if(k_flag != NULL)
+	{
+		kv_resize(uint8_t, k_flag->a, (uint64_t)len);
+		k_flag->a.n = len; 
+		memset(k_flag->a.a, 0, k_flag->a.n);
+	} 
+
+
+	assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
+	///sizeof(ha_mz1_t) = 16
+	memset(buf, 0xff, w * 16);
+	memset(&tq, 0, sizeof(tiny_queue_t));
+	///len/w is the evaluated minimizer numbers
+	kv_resize(ha_mz1_t, *p, p->n + len/w);
+
+	for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
+		int c = seq_nt4_table[(uint8_t)str[i]];
+		ha_mz1_t info = dummy;
+		if (c < 4) { // not an ambiguous base
+			int z;
+			if (is_hpc) {
+				int skip_len = 1;
+				if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
+					for (skip_len = 2; i + skip_len < len; ++skip_len)
+						if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
+							break;
+					i += skip_len - 1; // put $i at the end of the current homopolymer run
+				}
+				tq_push(&tq, skip_len);
+				kmer_span += skip_len;
+				///how many bases that are covered by this HPC k-mer
+				///kmer_span includes at most k HPC elements
+				if (tq.count > k) kmer_span -= tq_shift(&tq);
+			} else kmer_span = l + 1 < k? l + 1 : k; 
+			///kmer_span should be used for HPC k-mer
+			///non-HPC k-mer, kmer_span should be k
+			///kmer_span is used to calculate anchor pos on reverse complementary strand
+
+			if(k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC
+
+			kmer[0] = (kmer[0] << 1 | (c&1))  & mask;                  // forward k-mer
+			kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
+			kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1))  << shift1; // reverse k-mer
+			kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
+			if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
+			z = kmer[1] < kmer[3]? 0 : 1; // strand
+			++l;
+			if (l >= k && kmer_span < 256) {
+				uint64_t y;
+				y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
+
+				filtered = 0;
+				if(hf != 0) filtered = ha_ft_isflt(hf, y);
+				if(dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
+				///if (hf == 0 || ha_ft_isflt(hf, y) == 0)
+				if(filtered == 0)
+					info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
+				if(k_flag != NULL) k_flag->a.a[i]++;
+				if(k_flag != NULL && filtered == 1) k_flag->a.a[i]++;
+			}
+		} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
+
+
+		//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
+		//i is the real base iterator, while l is the HPC base iterator
+		//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
+		//but even if l < k, infor is still stored into buf
+		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
+		if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
+			for (j = buf_pos + 1; j < w; ++j)
+				if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
+			for (j = 0; j < buf_pos; ++j)
+				if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
+		}
+		/**
+		 * There are three cases:
+		 * 1. info.x <= min.x, means info is a new minimizer
+		 * 2. info.x > min.x, info is not a new minimizer
+		 *    (1) buf_pos != min_pos, do nothing
+		 *    (2) buf_pos == min_pos, means current minimizer has moved outside the window
+		 * **/
+		///three cases: 1. 
+		if (info.x <= min.x) { // a new minimum; then write the old min
+			if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			min = info, min_pos = buf_pos;
+		} else if (buf_pos == min_pos) { // old min has moved outside the window
+			if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			///buf_pos == min_pos, means current minimizer has moved outside the window
+			///so for now we need to find a new minimizer at the current window (w k-mers)
+			for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
+				if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
+			for (j = 0; j <= buf_pos; ++j)
+				if (min.x >= buf[j].x) min = buf[j], min_pos = j;
+
+			if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
+				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
+					if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
+				for (j = 0; j <= buf_pos; ++j)
+					if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
+			}
+		}
+		if (++buf_pos == w) buf_pos = 0;
+	}
+	if (min.x != UINT64_MAX)
+		kv_push(ha_mz1_t, *p, min);
+}
Author	SHA1	Message	Date
chhylp123	284cd0784a	Merge branch 'master' into hifiasm-v0.14	2021-02-09 18:46:29 -05:00
chhylp123	ceeb4562af	fix hc_edge	2021-02-09 18:44:06 -05:00
chhylp123	5bb22f701c	debug hic_hits	2021-02-09 18:40:49 -05:00
chhylp123	8467d58ef1	for v0.14	2021-02-09 13:00:23 -05:00
chhylp123	fda13cb0e4	update -z	2021-02-02 12:45:50 -05:00
chhylp123	02b9a3e109	hic multiple rounds	2021-01-31 18:49:08 -05:00
chhylp123	75ce7cd576	for v0.14-hic	2021-01-31 10:34:21 -05:00
chhylp123	a009026e65	fix bug for broken bubbles	2021-01-15 07:33:13 -05:00
chhylp123	88f5b11f04	bubble graph	2021-01-01 04:50:42 -05:00
chhylp123	c545fd84c7	hic bubble	2020-12-29 09:15:44 -05:00
chhylp123	5d506e82c9	phasing	2020-12-28 03:35:17 -05:00
chhylp123	53a655d802	init partition	2020-12-26 19:01:24 -05:00
chhylp123	f3e390eee8	for hic	2020-12-17 02:53:31 -05:00
chhylp123	73b5ef6769	bubble phasing	2020-11-21 20:50:31 -05:00
chhylp123	3a1a3750e9	Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly	2020-11-16 01:49:07 -05:00
chhylp123	d87314a07e	resolve conflict	2020-11-16 01:44:12 -05:00
chhylp123	ab073aa1fc	hic first	2020-11-14 23:35:32 -05:00
Heng Li	eddb8173da	added --min-hist-cnt for #49	2020-10-26 16:28:30 -04:00
chhylp123	35e3423f93	debug counting	2020-10-26 12:32:55 -04:00
chhylp123	0cc02035c5	HiC full index	2020-10-26 11:01:48 -04:00
chhylp123	b2b4624cea	rerelease v0.13	2020-10-22 14:24:59 -04:00
chhylp123	de4163e881	Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly merge README	2020-10-21 22:03:43 -04:00
chhylp123	53adde2b37	fix conflict	2020-10-21 22:01:07 -04:00
chhylp123	91d4915da7	release v0.13	2020-10-21 21:57:43 -04:00
chhylp123	04f94514ee	fix circle	2020-10-19 19:51:17 -04:00
chhylp123	b1f844801a	fix c->l bug	2020-10-19 16:39:53 -04:00
chhylp123	e5d503ea65	update polishing	2020-10-19 13:22:45 -04:00
chhylp123	9d4e3b3283	update polishing	2020-10-18 17:59:31 -04:00
chhylp123	753ea81194	reomve redundancy	2020-10-16 18:08:12 -04:00
chhylp123	89fcb35a64	for debugging	2020-10-11 14:32:33 -04:00
Heng Li	8cb131dca3	added a "Why hifiasm?" section	2020-10-05 01:01:58 -04:00
Heng Li	00128c3076	Fixed the wrong download link	2020-10-04 23:01:16 -04:00
Heng Li	519d5b1eb1	Updated the download link	2020-10-04 23:00:27 -04:00
chhylp123	7a314d1ef3	fix type in README	2020-08-27 23:57:48 -04:00
chhylp123	55f95c0191	fix typo	2020-08-26 00:25:42 -04:00
chhylp123	b629dd8556	freeze-0.12-r304	2020-08-23 23:38:04 -04:00
chhylp123	0a2559e4a0	freeze-r304	2020-08-23 23:36:45 -04:00
chhylp123	724990cbe7	freeze-r304	2020-08-23 23:34:36 -04:00
chhylp123	03cf0d475c	save bfore multiple overlaps	2020-08-23 08:46:33 -04:00
chhylp123	90d3cb69b9	fix version number	2020-08-05 18:04:47 -04:00
chhylp123	927d66262b	rerelase v0.11	2020-08-05 17:50:53 -04:00
chhylp123	e031c70d0f	Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly	2020-08-05 17:46:17 -04:00
chhylp123	a8aa37d201	fix makefile	2020-08-05 17:45:23 -04:00
Heng Li	ebedae2a21	updated Makefile	2020-08-05 17:33:14 -04:00
chhylp123	320dc58127	fix bug at r_utg	2020-08-04 16:43:50 -04:00
chhylp123	ef5ba0cc2c	update help/man	2020-07-31 21:32:02 -04:00
chhylp123	8bdc21799c	update READMEN	2020-07-31 12:42:27 -04:00
chhylp123	7c385e219a	update version number	2020-07-31 12:36:37 -04:00
chhylp123	d92f9519ae	Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly	2020-07-31 12:34:49 -04:00
chhylp123	6ca9d01f06	update version number	2020-07-31 12:28:30 -04:00
chhylp123	23f0a8b020	Merge pull request #23 from chhylp123/hifiasm_high_het Hifiasm high het	2020-07-31 12:24:43 -04:00
chhylp123	9d659fd2e0	update v0.10	2020-07-31 12:23:25 -04:00
chhylp123	ace34ae2b8	for high-het and coverage	2020-07-28 13:07:35 -04:00