seperate phase dp to ont/hifi

fix bubble issue

Assembly.cpp

@@ -12,6 +12,7 @@
 #include "kthread.h"
 #include "rcut.h"
 #include "kalloc.h"
+#include "ecovlp.h"
 
 void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
 int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp);
@@ -595,9 +596,9 @@ 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;
-    // if(i != 33) return;
+    // if(i != 12578) return;
     // fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
-    // if (memcmp("m64012_190920_173625/88015004/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
+    // if (memcmp("7897e875-76e5-42c8-bc37-94b370c4cc8d", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
     //     fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
     // } else {
     //     return;
@@ -606,6 +607,9 @@ static void worker_ovec(void *data, long i, int tid)
     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, NULL/**&(b->k_flag)**/, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
 
+    // prt_chain(&b->olist);
+    // return;
+
 	clear_Cigar_record(&b->cigar1);
 	clear_Round2_alignment(&b->round2);
 
@@ -649,13 +653,6 @@ static void worker_ovec_cal0(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
 	int fully_cov, abnormal;
-    // if(i != 33) return;
-    // fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
-    // if (memcmp("m64012_190920_173625/88015004/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
-    //     fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
-    // } else {
-    //     return;
-    // }
 
     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, NULL/**&(b->k_flag)**/, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
@@ -664,18 +661,25 @@ static void worker_ovec_cal0(void *data, long i, int tid)
 	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, &b->r_buf, &(b->tmp_region.w_list), 0, 0, &fully_cov, &abnormal);
+			&b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 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;
 
 
-
+	// R_INF.paf[i].is_fully_corrected = 0;
-	R_INF.paf[i].is_fully_corrected = 0;
+	// R_INF.paf[i].is_abnormal = abnormal;
-	R_INF.paf[i].is_abnormal = abnormal;
     R_INF.trio_flag[i] = AMBIGU;
 
+    // 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;
+    // R_INF.trio_flag[i] = AMBIGU;
+    
     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);
 
@@ -879,16 +883,14 @@ static void worker_ec_save(void *data, long i, int tid)
 
 void Output_corrected_reads()
 {
-    long long i;
+    uint64_t i; UC_Read g_read;
-    UC_Read g_read;
     init_UC_Read(&g_read);
     char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35);
     sprintf(gfa_name, "%s.ec.fa", asm_opt.output_file_name);
-    FILE* output_file = fopen(gfa_name, "w");
+    FILE *output_file = fopen(gfa_name, "w");
     free(gfa_name);
 
-    for (i = 0; i < (long long)R_INF.total_reads; i++)
+    for (i = 0; i < R_INF.total_reads; i++) {
-    {
         recover_UC_Read(&g_read, &R_INF, i);
         fwrite(">", 1, 1, output_file);
         fwrite(Get_NAME(R_INF, i), 1, Get_NAME_LENGTH(R_INF, i), output_file);
@@ -900,6 +902,40 @@ void Output_corrected_reads()
     fclose(output_file);
 }
 
+void Output_corrected_fastq()
+{
+    uint64_t i, k;
+    UC_Read g_read; asg8_v dv;
+    init_UC_Read(&g_read); kv_init(dv);
+    char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35);
+    sprintf(gfa_name, "%s.ec.fq", asm_opt.output_file_name);
+    FILE* fp = fopen(gfa_name, "w");
+    free(gfa_name);
+
+    for (i = 0; i < R_INF.tqn; i++) {
+        recover_UC_Read(&g_read, &R_INF, i);
+        fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
+        fprintf(fp, "%.*s\n", (int32_t)g_read.length, g_read.seq);        
+        fprintf(fp, "+\n");
+        retrive_bqual(&dv, NULL, i, -1, -1, 0, sc_bn);
+        for (k = 0; k < dv.n; k++) fprintf(fp, "%c", (char)(sc_tb[dv.a[k]] + 33 - 1));
+        fprintf(fp, "\n");
+    }
+
+    for (; i < R_INF.total_reads; i++) {
+        recover_UC_Read(&g_read, &R_INF, i);
+        fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
+        fprintf(fp, "%.*s\n", (int32_t)g_read.length, g_read.seq);        
+        fprintf(fp, "+\n");
+        // retrive_bqual(&dv, NULL, i, -1, -1, 0, sc_bn);
+        // for (k = 0; k < dv.n; k++) fprintf(fp, "%c", (char)(sc_tb[dv.a[k]] + 33 - 1));
+        for (k = 0; k < (uint64_t)g_read.length; k++) fprintf(fp, "%c", (char)(3 + 33 - 1));
+        fprintf(fp, "\n");
+    }
+    destory_UC_Read(&g_read); kv_destroy(dv);
+    fclose(fp);
+}
+
 void debug_print_pob_regions()
 {
     uint64_t i, total = 0;
@@ -920,7 +956,7 @@ 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_flt_tab_hp = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 1, 0);
     }
     ha_idx_hp = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 1, &R_INF, &hom_cov, &het_cov);
 
@@ -966,6 +1002,81 @@ void prt_dbg_rs(FILE *fp, Debug_reads* x, uint64_t round)
     destory_UC_Read(&g_read);
 }
 
+void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t *tot_e, uint64_t w_tmp)
+{
+	int hom_cov, het_cov, r_out = 0;
+    ha_flt_tab_hp = ha_idx_hp = NULL; (*tot_b) = (*tot_e) = 0;
+
+    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
+    
+    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
+        asm_opt.hom_cov = hom_cov; asm_opt.het_cov = het_cov;
+    }
+	///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);
+    het_cnt = NULL;
+    if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
+
+    if (r_out) {
+        write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
+        if((asm_opt.flag & HA_F_VERBOSE_GFA) && (asm_opt.bin_only == 1)) exit(1);///just for debug
+    }
+    if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0);
+
+    // Output_corrected_fastq();
+
+
+    cal_ec_r(asm_opt.thread_num, round, num_pround, R_INF.total_reads, (round == (asm_opt.number_of_round-1))?1:0, tot_b, tot_e);
+
+    // exit(1);    
+
+    // if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
+    if(des_idx) {
+        ha_pt_destroy(ha_idx); ha_idx = NULL;
+    }
+	
+
+    if(het_cnt) {
+        print_het_cnt_log(het_cnt); free(het_cnt); het_cnt = NULL;
+    }
+
+    // exit(1);
+
+
+    if (asm_opt.required_read_name) prt_dbg_rs(R_INF_FLAG.fp_r0, &R_INF_FLAG, 0); // for debugging only
+    
+	// save corrected reads to R_INF
+	// sl_ec_r(asm_opt.thread_num, R_INF.total_reads);
+
+    if (asm_opt.required_read_name) prt_dbg_rs(R_INF_FLAG.fp_r1, &R_INF_FLAG, 1); // for debugging only
+    if (asm_opt.required_read_name) destory_Debug_reads(&R_INF_FLAG), exit(0); // for debugging only
+    ///debug_print_pob_regions();
+
+    // Output_corrected_reads();
+
+    // exit(1);
+}
+
+
+int ha_ec_dbg(void)
+{
+	int hom_cov, het_cov;
+
+    ha_idx = ha_pt_gen(&asm_opt, 0, 0, 0, &R_INF, &hom_cov, &het_cov); // build the index
+    asm_opt.hom_cov = hom_cov; asm_opt.het_cov = het_cov;
+	ha_opt_update_cov(&asm_opt, hom_cov);
+
+    cal_ec_r_dbg(asm_opt.thread_num, R_INF.total_reads);
+
+    ha_pt_destroy(ha_idx); ha_idx = NULL;
+
+    return 0;
+}
 
 void ha_overlap_and_correct(int round)
 {
@@ -992,11 +1103,14 @@ void ha_overlap_and_correct(int round)
     het_cnt = NULL;
     if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
     // fprintf(stderr, "[M::%s-start]\n", __func__);
+    // double tt0 = yak_realtime_0();
 	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);///debug_for_fix
     // fprintf(stderr, "[M::%s-end]\n", __func__);
+    // fprintf(stderr, "[M::%s::%.3f] ==> chaining\n", __func__, yak_realtime_0()-tt0);
+    // exit(1);
     
     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);
@@ -1038,7 +1152,7 @@ void ha_overlap_and_correct(int round)
     ///debug_print_pob_regions();
 }
 
-void ha_overlap_cal(int round)
+void ha_overlap_cal(int round, int read_from_store)
 {
 	int i, hom_cov, het_cov;
 	ha_ovec_buf_t **b;
@@ -1049,9 +1163,9 @@ void ha_overlap_cal(int round)
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1),0);
     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
+	if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, ((round == 0)&&(read_from_store == 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
+    if (/**round == 0 &&**/ ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
 		ha_opt_update_cov(&asm_opt, hom_cov);
     het_cnt = NULL;
     if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
@@ -1075,6 +1189,8 @@ void ha_overlap_cal(int round)
 		ha_ovec_destroy(b[i]);
 	}
 	free(b);
+    asm_opt.hom_cov = hom_cov;
+    asm_opt.het_cov = het_cov;
 }
 
 
@@ -1787,7 +1903,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, 0);
+        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0, 0);
         ha_opt_update_cov(&asm_opt, hom_cov);
     }
     free(R_INF.read_length);
@@ -1836,6 +1952,30 @@ void ha_overlap_final(void)
     asm_opt.het_cov = het_cov;
 }
 
+void ha_ec_ff(int renew_idx)
+{
+	int hom_cov, het_cov;
+    ha_flt_tab_hp = ha_idx_hp = NULL;
+
+    if(ha_idx && renew_idx) {
+        ha_pt_destroy(ha_idx); ha_idx = NULL;
+    }
+
+    if(!ha_idx) {
+        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;
+    }
+
+    cal_ov_r(asm_opt.thread_num, R_INF.total_reads, renew_idx);
+
+    if(asm_opt.write_pos_idx) {
+        refresh_pt_idx(&ha_flt_tab, &ha_idx, NULL, &asm_opt, asm_opt.output_file_name, 1);
+        // write_pt_index(ha_flt_tab, ha_idx, NULL, &asm_opt, asm_opt.output_file_name);
+    } else {
+        ha_pt_destroy(ha_idx); ha_idx = NULL;
+    }	
+}
+
 static void worker_ov_utg(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
@@ -1896,15 +2036,14 @@ int ha_assemble_ovec(void)
     ha_flt_tab = ha_idx = NULL;
 
     // construct hash table for high occurrence k-mers
-    if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) 
+    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, 0);
-        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);
     ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
-    ha_overlap_cal(r);
+    ha_overlap_cal(r, 0);
     fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
     
     
@@ -1919,13 +2058,107 @@ int ha_assemble_ovec(void)
 	return 0;
 }
 
+int ha_assemble_ovec_cc(void)
+{
+    ha_idx = NULL;
+
+    ha_opt_reset_to_round(&asm_opt, 0); // this update asm_opt.roundID and a few other fields
+    ha_ec_dbg();
+
+    // Output_PAF0(R_INF.paf, "0");
+	destory_All_reads(&R_INF);
+	return 0;
+}
+
 int ha_assemble(void)
 {
     // debug_mc_g_t(MC_NAME);
     // debug_mc_gg_t(MC_NAME, 0, 0);
     // quick_debug_phasing(MC_NAME);
 	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
-	int r, hom_cov = -1, ovlp_loaded = 0;
+	int r, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e;
+	if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
+		ovlp_loaded = 1;
+		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
+		if (asm_opt.extract_list) {
+			ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list);
+			exit(0);
+		}
+		if (asm_opt.flag & HA_F_WRITE_EC) {
+            if(asm_opt.is_sc) Output_corrected_fastq();
+            else Output_corrected_reads();
+        }
+		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
+        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);
+        r = ha_idx?asm_opt.number_of_round-1:0;
+        if((!ha_idx) && (asm_opt.restart)) {
+            for (r = asm_opt.number_of_round - 1; r >= 0; --r) {
+                if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1)) {
+                    load_ct_index(&ha_ct_table, asm_opt.output_file_name); r0 = r;
+                    break;
+                }
+            }
+            if(r < 0) r = 0;
+        }
+
+		// construct hash table for high occurrence k-mers
+		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, 0);
+			ha_opt_update_cov(&asm_opt, hom_cov);
+		}
+		// error correction
+		assert(asm_opt.number_of_round > 0);
+		for (; 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
+            tot_b = tot_e = 0;
+			// ha_overlap_and_correct(r);
+            ha_ec(r, asm_opt.number_of_pround, (r<asm_opt.number_of_round-1)?1:0, &tot_b, &tot_e, ((r > r0) && (asm_opt.restart))?1:0);
+			fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
+					yak_cpu_usage(), yak_peakrss_in_gb(), r + 1);
+            fprintf(stderr, "[M::%s] # bases: %lu; # corrected bases: %lu\n", __func__, tot_b, tot_e);
+			// fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__,
+			// 		asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases);
+			// fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
+		}
+		if (asm_opt.flag & HA_F_WRITE_EC) {
+            if(asm_opt.is_sc) Output_corrected_fastq();
+            else Output_corrected_reads();
+        }
+		// overlap between corrected reads
+		ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round);
+		// ha_overlap_final();
+        ha_ec_ff(1/**0**/);
+        fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb());
+		// fprintf(stderr, "\n[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb());
+		// ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
+		if(!(asm_opt.write_pos_idx)) {
+            ha_ft_destroy(ha_flt_tab); ha_flt_tab = NULL;
+        }
+		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
+		ha_triobin(&asm_opt);
+
+        // exit(1);
+	}
+    if(ovlp_loaded == 2) ovlp_loaded = 0;
+    ha_opt_update_cov_min(&asm_opt, asm_opt.hom_cov, MIN_N_CHAIN);
+
+    build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
+        R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 
+        asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
+	destory_All_reads(&R_INF); if(asm_opt.dbg_bam) destroy_cc_v(&scb);
+	return 0;
+}
+
+
+int ha_assemble_pair(void)
+{
+	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
+	int r = 0, hom_cov = -1, ovlp_loaded = 0; memset((&R_INF), 0, sizeof(R_INF));
+
 	if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
 		ovlp_loaded = 1;
 		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
@@ -1937,37 +2170,44 @@ int ha_assemble(void)
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
         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(!append_All_reads(&R_INF, asm_opt.output_file_name, 0)) {
+            fprintf(stderr, "[M::%s::] Cannot load %s.0\n", __func__, asm_opt.output_file_name); 
+            exit(1);
+        }
+
+        if(!append_All_reads(&R_INF, asm_opt.output_file_name, 1)) {
+            fprintf(stderr, "[M::%s::] Cannot load %s.1\n", __func__, asm_opt.output_file_name); 
+            exit(1);
+        }
+        // Output_corrected_reads(); exit(0);
+
+        ha_flt_tab = ha_idx = NULL; r = asm_opt.number_of_round - 1;
         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 == NULL) 
+        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, 1);
-			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 = 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);
+        ha_overlap_cal(r, 1);
-			fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
-					yak_cpu_usage(), yak_peakrss_in_gb(), r + 1);
-			fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__,
-					asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases);
         fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
-		}
-		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
-		// overlap between corrected reads
-		ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round);
-		ha_overlap_final();
         fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(),
 				yak_cpu_usage(), yak_peakrss_in_gb());
+        
+
+        if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
         ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
         ha_ft_destroy(ha_flt_tab);
         if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
         ha_triobin(&asm_opt);
+
 	}
     if(ovlp_loaded == 2) ovlp_loaded = 0;
     ha_opt_update_cov_min(&asm_opt, asm_opt.hom_cov, MIN_N_CHAIN);

Assembly.h

@@ -44,11 +44,13 @@ typedef struct {
 } ha_ovec_buf_t;
 
 int ha_assemble(void);
+int ha_assemble_pair(void);
 void ug_idx_build(ma_ug_t *ug, int hap_n);
 ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug);
 ha_ovec_buf_t *ha_ovec_buf_init(void *km, int is_final, int save_ov, int is_ug);
 void ha_ovec_destroy(ha_ovec_buf_t *b);
 int64_t ha_ovec_mem(const ha_ovec_buf_t *b, int64_t *mem_a);
 int ha_assemble_ovec(void);
+int ha_ec_dbg(void);
 
 #endif

CommandLines.cpp

@@ -7,6 +7,9 @@
 #include <sys/time.h>
 #include "CommandLines.h"
 #include "ketopt.h"
+#include "kseq.h"
+
+KSEQ_INIT(gzFile, gzread)
 
 #define DEFAULT_OUTPUT "hifiasm.asm"
 
@@ -18,8 +21,8 @@ static ko_longopt_t long_options[] = {
 	{ "write-paf",     ko_no_argument, 302 },
 	{ "write-ec",      ko_no_argument, 303 },
 	{ "skip-triobin",  ko_no_argument, 304 },
-	{ "max-od-ec",     ko_no_argument, 305 },
+	{ "max-od-ec",     ko_required_argument, 305 },
-	{ "max-od-final",  ko_no_argument, 306 },
+	{ "max-od-final",  ko_required_argument, 306 },
 	{ "ex-list",       ko_required_argument, 307 },
 	{ "ex-iter",       ko_required_argument, 308 },
     { "hom-cov",     ko_required_argument, 309 },
@@ -64,6 +67,33 @@ static ko_longopt_t long_options[] = {
     { "ul-cut",     ko_required_argument, 349},
     { "dual-scaf",     ko_no_argument, 350},
     { "scaf-gap",     ko_required_argument, 351},
+    { "sec-in",     ko_required_argument, 352},
+    { "somatic-cov",     ko_required_argument, 353},
+    { "telo-m",     ko_required_argument, 354},
+    { "telo-p",     ko_required_argument, 355},
+    { "telo-d",     ko_required_argument, 356},
+    { "telo-s",     ko_required_argument, 357},
+    { "ctg-n",     ko_required_argument, 358},
+    { "ont",       ko_no_argument, 359},
+    // { "sc-n",       ko_no_argument, 360},
+    { "chem-c",     ko_required_argument, 361},
+    { "chem-f",     ko_required_argument, 362},
+    { "ul-m",     ko_required_argument, 363},
+    { "rl-cut",     ko_required_argument, 364},
+    { "sc-cut",     ko_required_argument, 365},
+    { "hf",     ko_required_argument, 366},
+    { "cb",     ko_required_argument, 367},
+    { "gpath",       ko_no_argument, 368},
+    { "het-cov",     ko_required_argument, 369},
+    { "resume",     ko_no_argument, 370},
+    { "flt-kocc",     ko_required_argument, 371},
+    { "chn-occ",     ko_required_argument, 372},
+    { "dbg-in1",     ko_required_argument, 373},
+    { "dbg-in2",     ko_required_argument, 374},
+    { "ec-only",     ko_no_argument, 375},
+    { "hyb-syn",     ko_required_argument, 376},
+    { "simd-m",     ko_required_argument, 377},
+    { "del-hf",     ko_no_argument, 378},
     // { "path-round",     ko_required_argument, 348},
 	{ 0, 0, 0 }
 };
@@ -84,11 +114,15 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    -t INT       number of threads [%d]\n", asm_opt->thread_num);
     fprintf(stderr, "    -h           show help information\n");
     fprintf(stderr, "    --version    show version number\n");
+    fprintf(stderr, "  Preset options:\n");
+    fprintf(stderr, "    --ont        assemble Oxford Nanopore reads\n");
 	fprintf(stderr, "  Overlap/Error correction:\n");
     fprintf(stderr, "    -k INT       k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
 	fprintf(stderr, "    -w INT       minimizer window size [%d]\n", asm_opt->mz_win);
 	fprintf(stderr, "    -f INT       number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
-	fprintf(stderr, "    -D FLOAT     drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
+	fprintf(stderr, "    -D FLOAT     drop k-mers occurring >FLOAT*coverage times [%.1f]; work with --flt-kocc or -N\n", asm_opt->high_factor);
+    fprintf(stderr, "    --flt-kocc   INT\n");
+    fprintf(stderr, "                 drop k-mers occurring >max(-D*coverage,--flt-kocc) times [%ld]\n", asm_opt->hf_cutoff);    
 	fprintf(stderr, "    -N INT       consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
     fprintf(stderr, "    -r INT       round of correction [%d]\n", asm_opt->number_of_round);
     fprintf(stderr, "    -z INT       length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
@@ -96,6 +130,16 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "                 employ k-mers occurring <INT times to rescue repetitive overlaps [%d]\n", asm_opt->max_kmer_cnt);
     fprintf(stderr, "    --hg-size    INT(k, m or g)\n");
     fprintf(stderr, "                 estimated haploid genome size used for inferring read coverage [auto]\n");
+    fprintf(stderr, "    --resume     resume from the previously incomplete assembly [%ld]\n", asm_opt->restart);
+    fprintf(stderr, "    --het-cov    INT\n");
+    fprintf(stderr, "                 heterozygous read coverage [auto]; used for error correction and assembly; manual value overrides auto\n");
+    fprintf(stderr, "    --hom-cov    INT\n");
+    fprintf(stderr, "                 homozygous read coverage [auto]; used for error correction and assembly; manual value overrides auto\n");
+    fprintf(stderr, "    --chn-occ    INT\n");
+    fprintf(stderr, "                 discard overlaps supported by <INT minimizers [%ld]\n", asm_opt->chn_occ);
+    fprintf(stderr, "    --ec-only    error correction only; disable overlapping and assembly\n");
+    fprintf(stderr, "    --simd-m     use SIMD acceleration when supported: AVX-512 (2), AVX2 (1), or non-SIMD (0)\n");
+    
     fprintf(stderr, "  Assembly:\n");
     fprintf(stderr, "    -a INT       round of assembly cleaning [%d]\n", asm_opt->clean_round);
     fprintf(stderr, "    -m INT       pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
@@ -107,8 +151,6 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    -u           post-join step for contigs which may improve N50; 0 to disable; 1 to enable\n");
     fprintf(stderr, "                 [%u] and [%u] in default for the UL+HiFi assembly and the HiFi assembly, respectively\n",
                                       asm_opt->ul_pst_join, asm_opt->hifi_pst_join);
-    fprintf(stderr, "    --hom-cov    INT\n");
-    fprintf(stderr, "                 homozygous read coverage [auto]\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");
@@ -119,6 +161,10 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "                 break contigs at positions with <=FLOAT*coverage exact overlaps;\n");
     fprintf(stderr, "                 only work with '--b-cov' or '--h-cov'[%.2f]\n", asm_opt->m_rate);
     fprintf(stderr, "    --primary    output a primary assembly and an alternate assembly\n");
+    fprintf(stderr, "    --ctg-n      INT\n");
+    fprintf(stderr, "                 remove tip contigs composed of <=INT reads [%d]\n", asm_opt->max_contig_tip);
+    fprintf(stderr, "    --gpath      output the corresponding path of each contig (p_ctg) within the assembly graph (d_utg.noseq.gfa)\n");
+    
     
 //	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");
@@ -167,6 +213,11 @@ void Print_H(hifiasm_opt_t* asm_opt)
 
     fprintf(stderr, "  Ultra-Long-integration:\n");
     fprintf(stderr, "    --ul FILEs   file names of Ultra-Long reads [r1.fq,r2.fq,...]\n");
+    ///pending for integration
+    /**
+    fprintf(stderr, "    --ul-m       INT\n");
+    fprintf(stderr, "                 hybrid assembly mode. 0: fast and memory efficent; 1: may produce better assembly with ONT R10 [%d]\n", asm_opt->ul_mod);
+    **/
     fprintf(stderr, "    --ul-rate    FLOAT\n");
     fprintf(stderr, "                 error rate of Ultra-Long reads [%.3g]\n", asm_opt->ul_error_rate);
     fprintf(stderr, "    --ul-tip     INT\n");
@@ -186,6 +237,37 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    --scaf-gap   INT\n");
     fprintf(stderr, "                 max gap size for scaffolding [%ld]\n", asm_opt->self_scaf_gap_max);
 
+    fprintf(stderr, "  Telomere-identification:\n");
+    fprintf(stderr, "    --telo-m     STR\n");
+    fprintf(stderr, "                 telomere motif at 5'-end; CCCTAA for human [%s]\n", ((asm_opt->telo_motif)?(asm_opt->telo_motif):("NULL")));///5'-end, check CCCTAA
+    fprintf(stderr, "    --telo-p     INT\n");
+    fprintf(stderr, "                 non-telomeric penalty [%ld]\n", asm_opt->telo_pen);
+    fprintf(stderr, "    --telo-d     INT\n");
+    fprintf(stderr, "                 max drop [%ld]\n", asm_opt->telo_drop);
+    fprintf(stderr, "    --telo-s     INT\n");
+    fprintf(stderr, "                 min score for telomere reads [%ld]\n", asm_opt->telo_mic_sc);
+
+    fprintf(stderr, "  ONT Simplex assembly (beta):\n");
+    fprintf(stderr, "    --ont        assemble ONT Simplex reads in fastq format\n");
+    // fprintf(stderr, "    --sc-n       consider base qual value for assembly\n");
+    fprintf(stderr, "    --chem-c     INT\n");
+    // fprintf(stderr, "                 detect chimeric reads with <=INT other reads support [%lu]\n", asm_opt->chemical_cov);
+    fprintf(stderr, "                 detect chimeric reads with <=INT other reads support [auto]\n");
+    fprintf(stderr, "    --chem-f     INT\n");
+    fprintf(stderr, "                 length of flanking regions for chimeric read detection [%lu]\n", asm_opt->chemical_flank);
+    fprintf(stderr, "    --rl-cut     INT\n");
+    fprintf(stderr, "                 filter out ONT Simplex reads shorter than <INT> for assembly [%ld]\n", asm_opt->rl_cut);
+    fprintf(stderr, "    --sc-cut     INT\n");
+    fprintf(stderr, "                 filter out ONT Simplex reads with a mean base quality score below <INT> [%ld]\n", asm_opt->sc_cut);
+    fprintf(stderr, "    --hf FILEs    HiFi read file(s)\n");
+    fprintf(stderr, "    --hyb-syn    INT\n");
+    fprintf(stderr, "                 hybrid correction mode (requires --hf) [%d]:\n", asm_opt->hyb_syn);
+    fprintf(stderr, "                 1: all-vs-all (ONT<-all, HiFi<-all)\n");
+    fprintf(stderr, "                 2: ONT<-all, HiFi<-HiFi\n");
+    fprintf(stderr, "                 3: mode 2 + ONT/HiFi sync to reduce bias\n");
+
+
+
     fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
     fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n");
 }
@@ -204,6 +286,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->hic_reads[0] = NULL;
     asm_opt->hic_reads[1] = NULL;
     asm_opt->fn_bin_poy = NULL;
+    asm_opt->fn_chr_bin = NULL;
     asm_opt->ar = NULL;
     asm_opt->thread_num = 1;
     asm_opt->k_mer_length = 51;
@@ -220,6 +303,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
 	asm_opt->max_kmer_cnt = 2000;
 	asm_opt->high_factor = 5.0;
 	asm_opt->max_ov_diff_ec = 0.04;
+    asm_opt->max_ov_diff_ec_sec = 0.04;
 	asm_opt->max_ov_diff_final = 0.03;
 	asm_opt->hom_cov = 20;
     asm_opt->het_cov = -1024;
@@ -228,6 +312,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->load_index_from_disk = 1;
     asm_opt->write_index_to_disk = 1;
     asm_opt->number_of_round = 3;
+    asm_opt->number_of_pround = 0/**3**/;
     asm_opt->adapterLen = 0;
     asm_opt->clean_round = 4;
     ///asm_opt->small_pop_bubble_size = 100000;
@@ -242,6 +327,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->min_overlap_coverage = 0;
     asm_opt->max_short_tip = 3;
     asm_opt->max_short_ul_tip = 6;
+    asm_opt->max_contig_tip = 3;
     asm_opt->min_cnt = 2;
     asm_opt->mid_cnt = 5;
     asm_opt->purge_level_primary = 3;
@@ -305,6 +391,59 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->self_scaf_min = 250000;
     asm_opt->self_scaf_reliable_min = 5000000;
     asm_opt->self_scaf_gap_max = 3000000;
+    asm_opt->sec_in = NULL;
+    asm_opt->somatic_cov = -1;
+
+    asm_opt->telo_motif = NULL;
+    asm_opt->telo_pen = 1;
+    asm_opt->telo_drop = 2000;
+    asm_opt->telo_mic_sc = 500;
+
+    asm_opt->is_ont = 0;
+    asm_opt->is_sc = 0;
+    asm_opt->chemical_cov = -1/**1**/;
+    asm_opt->chemical_flank = 256;
+    asm_opt->ul_mod = 0;
+
+    asm_opt->rl_cut = 1000;
+    asm_opt->sc_cut = 10;
+
+    asm_opt->hf = NULL;
+
+    asm_opt->gpath = 0;
+
+    asm_opt->hf_rate = 4;
+    asm_opt->ont_rate = 1;///must be 1 or 0
+    asm_opt->hf_rate_max = 4;
+
+    asm_opt->het_cov_set = -1;
+    asm_opt->restart = 0;
+
+    asm_opt->hf_cutoff = -1;
+
+    asm_opt->write_pos_idx = 0/**1**/;
+
+    asm_opt->hom_cov_0 = -1;
+    asm_opt->het_cov_0 = -1;
+	asm_opt->max_n_chain_0 = -1;
+
+
+    asm_opt->hmo_cov_ss = -1;
+    asm_opt->het_cov_ss = -1;
+    asm_opt->chn_occ = 2;
+
+    asm_opt->dbg_run_1 = NULL;
+    asm_opt->dbg_run_2 = NULL;
+
+    asm_opt->ec_only = 0;
+    asm_opt->hyb_syn = 1;
+    asm_opt->step_rd = -1/**128**/;
+
+    asm_opt->dbg_bam = 0;
+
+    asm_opt->simd_mm = -1;
+
+    asm_opt->del_hf = 0;
 }   
 
 void destory_enzyme(enzyme* f)
@@ -378,14 +517,43 @@ static int check_file(char* name, const char* opt)
 
 static int check_hic_reads(enzyme* f, const char* opt)
 {
-    int i;
+    int32_t i;
-    for (i = 0; i < f->n; i++)
+    for (i = 0; i < f->n; i++) {
-    {
         if(check_file(f->a[i], opt) == 0) return 0;
     }
     return 1;
 }
 
+static int check_fq_files(enzyme* f, const char* opt, int32_t is_fq)
+{
+    int32_t i, ret; gzFile dfp; kseq_t *ks = NULL;
+    for (i = 0; i < f->n; i++) {
+        if(!(f->a[i])) {
+            fprintf(stderr, "[ERROR] input file does not exist (%s)\n", opt);
+            return 0;
+        }
+
+        dfp = gzopen(f->a[i], "r");
+        if (dfp == 0) {
+            fprintf(stderr, "[ERROR] Cannot find the input file: %s (%s)\n", f->a[i], opt);
+            return 0;
+        } else if(is_fq){
+            ks = kseq_init(dfp);
+            while (((ret = kseq_read(ks)) >= 0)) {
+                if((ks->qual.l == 0) || (ks->qual.s == NULL)) {
+                    fprintf(stderr, "[ERROR] %s is in fasta format rather than fastq format (%s)\n", f->a[i], opt);
+                    fprintf(stderr, "[ERROR] set --ul-m 0 for fasta files\n");
+                    return 0;
+                }
+                break;
+            }
+            kseq_destroy(ks); ks = NULL;
+        }
+        gzclose(dfp);
+    }
+    return 1;
+}
+
 int check_option(hifiasm_opt_t* asm_opt)
 {
     if(asm_opt->read_file_names == NULL || asm_opt->num_reads == 0)
@@ -585,7 +753,7 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
-    if(asm_opt->ar != NULL && check_hic_reads(asm_opt->ar, "UL") == 0) return 0;
+    if(asm_opt->ar != NULL && check_fq_files(asm_opt->ar, "--ul", asm_opt->ul_mod) == 0) return 0;
     if(asm_opt->ar != NULL && asm_opt->ar->n == 0)
     {
         fprintf(stderr, "[ERROR] wrong UL reads (--ul)\n");
@@ -634,30 +802,70 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+    if(asm_opt->ul_mod != 0 && asm_opt->ul_mod != 1) {
+        fprintf(stderr, "[ERROR] must be 0 or 1 (--ul-m)\n");
+        return 0;
+    }
+
+    if(asm_opt->telo_motif) {
+        uint64_t k, tlen = strlen((asm_opt->telo_motif)); char c;
+        if(tlen > 32) {
+            fprintf(stderr, "[ERROR] [--telo-m] must be no longer than 32\n");
+            return 0;
+        }
+        for (k = 0; k < tlen; k++) {
+            c = asm_opt->telo_motif[k];
+            if(c != 'A' && c != 'C' && c != 'G' && c != 'T' &&
+                c != 'a' && c != 'c' && c != 'g' && c != 't') {
+                fprintf(stderr, "[ERROR] [--telo-m] must be A/C/G/T\n");
+                return 0;
+            }
+        }
+    }
+
+
+    if((asm_opt->hf) && (!(asm_opt->is_ont))) {
+        fprintf(stderr, "[ERROR] [--hf] must work with [--ont]\n");
+        return 0;
+    }
+
+    if((asm_opt->hyb_syn != 1) && (asm_opt->hyb_syn != 2) && (asm_opt->hyb_syn != 3)) {
+        fprintf(stderr, "[ERROR] [--hyb-syn] must be 1/2/3\n");
+        return 0;
+    }
+
     return 1;
 }
 
 void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt)
 {
-    if(opt->ind == argc)
+    if(opt->ind == argc) {
-    {
         return;
     }
 
     asm_opt->num_reads = argc - opt->ind;
     asm_opt->read_file_names = (char**)malloc(sizeof(char*)*asm_opt->num_reads);
     
-    long long i;
+    long long i; int ret;
-    gzFile dfp;
+    gzFile dfp; kseq_t *ks = NULL;
-    for (i = 0; i < asm_opt->num_reads; i++)
+    for (i = 0; i < asm_opt->num_reads; i++) {
-    {
         asm_opt->read_file_names[i] = argv[i + opt->ind];
         dfp = gzopen(asm_opt->read_file_names[i], "r");
-        if (dfp == 0)
+        if (dfp == 0) {
-        {
             fprintf(stderr, "[ERROR] Cannot find the input read file: %s\n", 
                     asm_opt->read_file_names[i]);
 		    exit(0);
+        } else if(asm_opt->is_sc){
+            ks = kseq_init(dfp);
+            while (((ret = kseq_read(ks)) >= 0)) {
+                if((ks->qual.l == 0) || (ks->qual.s == NULL)) {
+                    fprintf(stderr, "[ERROR] %s is in fasta format rather than fastq format\n", asm_opt->read_file_names[i]);
+                    asm_opt->is_sc = 0;
+                    exit(0);
+                }
+                break;
+            }
+            kseq_destroy(ks); ks = NULL;
         }
         gzclose(dfp);
     }
@@ -794,12 +1002,14 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg);
 		else if (c == 307) asm_opt->extract_list = opt.arg;
 		else if (c == 308) asm_opt->extract_iter = atoi(opt.arg);
-        else if (c == 309)
+        else if (c == 309) {
-        {
+            asm_opt->hom_global_coverage = asm_opt->hmo_cov_ss = atoi(opt.arg);
-            asm_opt->hom_global_coverage = atoi(opt.arg);
             asm_opt->hom_global_coverage_set = 1;
+            if(asm_opt->hmo_cov_ss <= 0) {
+                fprintf(stderr, "[ERROR] homozygous read coverage should be > 0 (--hom-cov)");
+			    return 1;
             }
-        else if (c == 310)
+        } else if (c == 310)
         {
             char* s = NULL;
             asm_opt->recover_atg_cov_min = strtol(opt.arg, &s, 10);
@@ -853,7 +1063,65 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 349) asm_opt->ul_min_base = atol(opt.arg);
         else if (c == 350) asm_opt->self_scaf = 1;
         else if (c == 351) asm_opt->self_scaf_gap_max = atol(opt.arg);
-        else if (c == 'l') {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
+        else if (c == 352) get_hic_enzymes(opt.arg, &(asm_opt->sec_in), 0);
+        else if (c == 353) asm_opt->somatic_cov = atol(opt.arg);
+        else if (c == 354) asm_opt->telo_motif = opt.arg;
+        else if (c == 355) asm_opt->telo_pen = atol(opt.arg);
+        else if (c == 356) asm_opt->telo_drop = atol(opt.arg);
+        else if (c == 357) asm_opt->telo_mic_sc = atol(opt.arg);
+        else if (c == 358) asm_opt->max_contig_tip = atol(opt.arg);
+        else if (c == 359) {
+            asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; asm_opt->is_sc = 1; ///asm_opt->mz_win = 37; asm_opt->k_mer_length = 37;
+        } /**else if (c == 360) {
+            asm_opt->is_sc = 1; 
+        }**/ else if (c == 361) {
+            asm_opt->chemical_cov = atol(opt.arg);
+        } else if (c == 362) {
+            asm_opt->chemical_flank = atol(opt.arg);
+        ///pending for integration
+        /**        
+        } else if (c == 363) {
+            asm_opt->ul_mod = atol(opt.arg); 
+        **/
+        } else if (c == 364) {
+            asm_opt->rl_cut = atol(opt.arg);
+        } else if (c == 365) {
+            asm_opt->sc_cut = atol(opt.arg);
+        } else if (c == 366) {
+            get_hic_enzymes(opt.arg, &(asm_opt->hf), 0);
+        } else if (c == 367) {
+            asm_opt->fn_chr_bin = opt.arg;
+        } else if (c == 368) {
+            asm_opt->gpath = 1;
+        } else if (c == 369) {
+            asm_opt->het_cov_set = asm_opt->het_cov_ss = atoi(opt.arg);
+            if(asm_opt->het_cov_ss <= 0) {
+                fprintf(stderr, "[ERROR] heterozygous read coverage should be > 0 (--het-cov)");
+			    return 1;
+            }
+        } else if (c == 370) {
+            asm_opt->restart = 1;
+        } else if (c == 371) {
+            asm_opt->hf_cutoff = atoi(opt.arg);
+        } else if (c == 372) {
+            asm_opt->chn_occ = atoi(opt.arg);
+            if(asm_opt->chn_occ <= 0) {
+                fprintf(stderr, "[ERROR] chain cutoff should be > 0 (--chn-occ)");
+			    return 1;
+            }
+        } else if (c == 373) {
+            asm_opt->dbg_run_1 = opt.arg;
+        } else if (c == 374) {
+            asm_opt->dbg_run_2 = opt.arg;
+        } else if (c == 375) {
+            asm_opt->ec_only = 1;
+        } else if (c == 376) {
+            asm_opt->hyb_syn = atoi(opt.arg);
+        } else if (c == 377) {
+            asm_opt->simd_mm = atoi(opt.arg);
+        } else if (c == 378) {
+            asm_opt->del_hf = 1;
+        } 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);
         }
         else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg);
@@ -890,6 +1158,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 
     // fprintf(stderr, "[M::%s::] post join::%u\n", __func__, (uint32_t)(!(asm_opt->flag & HA_F_BAN_POST_JOIN)));
     // exit(1);
+    if(!(asm_opt->is_ont)) {
+        asm_opt->rl_cut = -1; asm_opt->sc_cut = 1;
+    }
 
     return check_option(asm_opt);
 }

CommandLines.h

@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.19.7-r598"
+#define HA_VERSION "0.25.1-r910"
 
 #define VERBOSE 0
 
@@ -41,10 +41,13 @@ typedef struct {
 	char *fn_bin_yak[2];
 	char *fn_bin_list[2];
     char *fn_bin_poy;
+    char *fn_chr_bin;
 	char *extract_list;
     enzyme *hic_reads[2];
     enzyme *hic_enzymes;
     enzyme *ar;
+    enzyme *hf;
+    enzyme *sec_in;
 	int extract_iter;
     int thread_num;
     int k_mer_length;
@@ -62,6 +65,7 @@ typedef struct {
 	int max_kmer_cnt;
 	double high_factor; // coverage cutoff set to high_factor*hom_cov
 	double max_ov_diff_ec;
+    double max_ov_diff_ec_sec;
 	double max_ov_diff_final;
 	int hom_cov;
     int het_cov;
@@ -73,6 +77,7 @@ typedef struct {
     int load_index_from_disk;
     int write_index_to_disk;
     int number_of_round;
+    int number_of_pround;
     int adapterLen;
     int clean_round;
     int roundID;
@@ -82,6 +87,7 @@ typedef struct {
     int min_overlap_coverage;
     int max_short_tip;
     int max_short_ul_tip;
+    int max_contig_tip;
     int min_cnt;
     int mid_cnt;
     int purge_level_primary;
@@ -134,6 +140,7 @@ typedef struct {
     int64_t infor_cov, s_hap_cov, trio_cov_het_ovlp;
     double ul_error_rate, ul_error_rate_low, ul_error_rate_hpc;
     int32_t ul_ec_round;
+    int32_t ul_mod;
     uint8_t is_dbg_het_cnt;
     uint8_t is_low_het_ul;
     uint8_t is_base_trans;
@@ -141,6 +148,7 @@ typedef struct {
     uint8_t is_topo_trans;
     uint8_t is_bub_trans;
     uint8_t bin_only;
+    uint8_t ec_only;
     int32_t ul_clean_round;
     int32_t prt_dbg_gfa;
     int32_t integer_correct_round;
@@ -151,6 +159,51 @@ typedef struct {
     uint64_t self_scaf_min;
     uint64_t self_scaf_reliable_min;
     int64_t self_scaf_gap_max;
+    int64_t somatic_cov;
+
+    char *telo_motif;
+    int64_t telo_pen;
+    int64_t telo_drop;
+    int64_t telo_mic_sc;
+    
+    uint64_t is_ont;
+    uint64_t is_sc;
+    int64_t chemical_cov;
+    uint64_t chemical_flank;
+
+    int64_t rl_cut;
+    int64_t sc_cut;
+    uint8_t gpath;
+
+    uint64_t hf_rate;///cannot be larger than 128?
+    uint64_t hf_rate_max;///cannot be larger than 128?
+    uint64_t ont_rate;///cannot be 0, should be 1 in anyway
+
+    int64_t het_cov_set;
+    int64_t restart;
+
+    int64_t hf_cutoff;
+
+    uint8_t write_pos_idx;
+
+    int hom_cov_0;
+    int het_cov_0;
+	int max_n_chain_0; // fall-back max number of chains to consider
+
+    int64_t hmo_cov_ss;
+    int64_t het_cov_ss;
+    int64_t chn_occ;
+
+    char *dbg_run_1, *dbg_run_2;
+
+    int32_t hyb_syn;
+
+    int64_t step_rd;
+
+    uint8_t dbg_bam;
+
+    int8_t simd_mm;
+    int8_t del_hf;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;

Correct.h

@@ -144,7 +144,10 @@ typedef struct
     char misBase;
 }haplotype_evdience;
 
-
+#define hh_tp(z) (((z).type&1))
+#define hh_hp(z) ((((z).type>>1)&1))
+#define hh_bq(z) ((((z).type>>2))&sc_bm)
+#define hh_wq(z) ((((z).type>>(sc_bn+2)))&sc_bm)
 
 typedef struct
 {
@@ -1384,8 +1387,119 @@ typedef struct {
 	int64_t k, q[2], t[2], cq[2], ct[2], ci[2], werr, werr0, cerr;
 	int64_t qoff, f, toff, coff, cur_qoff;
 } rtrace_iter;
+
+typedef struct
+{
+    overlap_region_alloc *ol;
+    Candidates_list *cl;
+    All_reads *rref;
+    UC_Read *qu; 
+    UC_Read *tu;
+    bit_extz_t *exz;
+    overlap_region *aux_o, *rse_o;
+    double e_rate[2];
+    int64_t wl[2];
+    int64_t rid; 
+    int64_t khit;
+    int64_t move_gap;
+    asg16_v *buf; 
+    asg64_v *srt; 
+    asg8_v *hpz;
+    ma_hit_t_alloc *in;
+
+    int8_t chem_drop[2]; 
+    double align_gap_rate[2]; 
+    int64_t align_gap_max[2];
+
+    uint64_t sec_aln_win; 
+    uint64_t sec_aln_cov; 
+    double sec_aln_err_rate; 
+    double sec_aln_max; 
+    asg64_v *kp;
+
+    asg32_v *v32; 
+    asg64_v *bp; 
+    ha_abuf_t *ab; 
+    uint64_t max_n_chain; 
+    uint64_t max_n_chain_f; 
+    uint64_t chain_cutoff; 
+    uint64_t ave_cov_min;
+    uint64_t ocw;
+
+    uint64_t t_cut;
+
+    uint64_t hom_cov_a;
+} gen_hc_aln_t;
 int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *trace, rtrace_t *tc,
 const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o, 
 bit_extz_t *exz, double e_rate, int64_t qs);
 
+uint64_t gen_hc_r_alin_adp_smp_ff_ec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, int64_t rid, asg64_v *sp, uint64_t ocw, uint32_t *ocn, uint32_t *osc, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min);
+uint64_t gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, 
+    asg64_v *sp, uint64_t ocw, uint8_t *hpf, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match, uint8_t is_dedup);
+uint64_t gen_hc_r_alin_adp_mmp_1(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, 
+    asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match);
+uint64_t gen_hc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, 
+    asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match);
+uint64_t gen_gc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, ul_idx_t *uref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, 
+    asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match);
+void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, uint8_t set_match);
+void gen_hc_r_alin_adv_adp_smp_0(gen_hc_aln_t *ez, uint8_t set_match);
+void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match);
+void pp_chn_a(overlap_region *z, Candidates_list *cl, uint8_t is_raw);
+uint64_t gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, uint8_t *hpf);
+void gen_hc_r_alin_flt(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, 
+    asg64_v *kp, asg64_v *sp, uint64_t ocw, uint8_t *hpf, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min);
+void gen_hc_r_alin_adp(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, 
+    asg64_v *kp, asg64_v *sp, uint64_t ocw, uint8_t *hpf, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min);
+void gen_hc_r_alin_adv(gen_hc_aln_t *ez);
+uint64_t gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, uint8_t *hpf);
+void gen_hc_r_alin_nec_adv(gen_hc_aln_t *ez);
+uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
+void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32, 
+    int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate);
+void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te);
+void push_alnw(overlap_region *aux_o, bit_extz_t *exz);
+void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez);
+void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *va, uint64_t scw, uint64_t *tqual, uint64_t *wqual);
+void gen_reseed_re(overlap_region_alloc *ol, Candidates_list *cl, overlap_region *aux_o, overlap_region *rse_o, All_reads *rref, UC_Read* qu, UC_Read *tu, bit_extz_t *exz, kv_ul_ov_t *c_idx, asg64_v *idx, asg64_v *res, int64_t bd, int64_t mzw, int64_t kl, int64_t rid, double err_h, double err_l, asg16_v *b16, uint64_t tqn, uint8_t *hpf);
+inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl, int64_t qs, int64_t qe, int64_t ts, int64_t te)
+{
+    if(qe - qs != te - ts) return 0;
+    if(memcmp(qstr + qs, tstr + ts, qe - qs) == 0) return 1;
+    return 0;
+}
+// void est_rep_err_rate(overlap_region_alloc* ol, asg64_v *ix, kv_ul_ov_t *c_idx, int64_t ql, int64_t wl, uint64_t *ou_a, uint64_t min_dp, int64_t ph_cov, uint8_t flg_ov, double flg_ov_sec_rate, double flg_cov_rate, uint64_t *ave_e, uint64_t *bd_e, uint64_t *tot_cov);
+void est_rep_err_rate(overlap_region_alloc* ol, asg64_v *ix, kv_ul_ov_t *c_idx, int64_t ql, int64_t wl, uint64_t *ou_a, uint64_t min_dp, int64_t ph_cov, uint8_t flg_ov, double flg_ov_sec_rate, double flg_cov_rate, uint64_t *ave_e, uint64_t *bd_e, uint64_t *tot_cov);
+#define ovlp_id(x) ((x).tn)
+#define ovlp_min_wid(x) ((x).ts)
+#define ovlp_max_wid(x) ((x).te)
+#define ovlp_cur_wid(x) ((x).qn)
+#define ovlp_cur_xoff(x) ((x).qs)
+#define ovlp_cur_yoff(x) ((x).ts)
+#define ovlp_cur_ylen(x) ((x).te)
+#define ovlp_cur_coff(x) ((x).qe)
+#define ovlp_bd(x) ((x).sec)
+#define ovlp_um(x) ((x).sec)
+#define ovlp_hf(x) ((x).el)
+
+#define HPC_PL 12
+#define HPC_RR 4
+#define HPC_CC 2
+#define HPC_RR_Q 5
+#define HPC_CC_Q 3
+#define HC_MF_R 0.5
+#define HC_AV_MIN 0.7
+#define GC_MF_N 12
+
+// #define FORCE_CUT 1
+
+// overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, 
+
+// double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max,
+//     uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp
+
+// uint64_t gen_hc_r_alin_ea_hybrid(overlap_region_alloc* ol, uint64_t bi, uint64_t bn, uint64_t tk, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, uint8_t ec_filter, ma_hit_t_alloc *in, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max,
+//     uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp)
+
 #endif

Hash_Table.cpp

@@ -1512,6 +1512,34 @@ inline int32_t comput_sc_ch(const k_mer_hit *ai, const k_mer_hit *aj, double bw_
     return sc;
 }
 
+inline int32_t comput_sc_ch_ec(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol)
+{
+    ///ai is the suffix of aj
+    int32_t dq, dr, dd, dg, q_span, sc; 
+    dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
+    if(dq <= 0) return INT32_MIN;
+    dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
+    if(dr <= 0) return INT32_MIN;
+    dd = dr > dq? dr - dq : dq - dr;//gap
+    if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN;
+    dg = dr < dq? dr : dq;//len
+    q_span = ai->cnt&(0xffu); 
+    sc = q_span < dg? q_span : dg;
+    sc = normal_w(sc, ((int32_t)(ai->cnt>>8)));
+    if (dd || (dg > q_span && dg > 0)) {
+        double lin_pen, a_pen;
+        lin_pen = (chn_pen_gap*(double)dd);
+        a_pen = ((double)(sc))*((((double)dd)/((double)dg))/bw_rate);
+        ///for long gap
+        // if(lin_pen > a_pen) lin_pen = a_pen;
+        if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
+        else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
+        lin_pen += (chn_pen_skip*(double)dg);
+        sc -= (int32_t)lin_pen;
+    }
+    return sc;
+}
+
 inline int32_t comput_sc_ff(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol)
 {
     ///ai is the suffix of aj
@@ -1721,18 +1749,7 @@ uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, o
     return cL;
 }
 
-
+void push_ovlp_chain_qgen(overlap_region* o, uint32_t xid, int64_t xl, int64_t yl, int64_t sc, k_mer_hit *beg, k_mer_hit *end)
-#define kv_pushp_ol(type, v, p) do {									\
-		if ((v).length == (v).size) {										\
-			(v).list = (type*)realloc((v).list, sizeof(type)*((v).size?((v).size<<1):(2)));	\
-            memset((v).list+(v).size, 0, sizeof(overlap_region)*(((v).size?((v).size<<1):2)-(v).size));\
-            (v).size = (v).size?((v).size<<1):(2);							\
-		}															\
-		*(p) = &((v).list[(v).length++]); \
-	} while (0)
-
-void push_ovlp_chain_qgen(overlap_region* o, uint32_t xid, int64_t xl, int64_t yl, int64_t sc, 
-k_mer_hit *beg, k_mer_hit *end)
 {
     int64_t xr, yr; 
     o->x_id = xid; o->y_id = beg->readID;
@@ -1986,6 +2003,284 @@ uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64
     return cL;
 }
 
+void quick_ck_lchain(k_mer_hit* a, int64_t a_n, int64_t xl, int64_t yl, double chn_pen_gap, double chn_pen_skip, double bw_rate, 
+int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int64_t *plus, int64_t *msc, int64_t *msc_i, int64_t *movl, int64_t *si, int64_t *ei)
+{
+    if(a_n <= 0) return;
+    int64_t l, k, is_srt = 1, z; k_mer_hit *ai, *aj;
+    int64_t dq, dr, dd, dg, q_span, sc, csc, ddt; 
+    int64_t plus0, msc0, msc_i0, movl0; double lin_pen, a_pen;
+
+    *plus = 0; *msc = *msc_i = INT32_MIN; *movl = INT32_MAX; *si = 0; *ei = a_n; 
+
+    for (k = 1, l = 0; k <= a_n; k++) {
+        if(k == a_n || a[k].strand != a[l].strand) {
+            t[k-1] = 0; ii[k-1] = 0;
+            // if(a_n && a[0].readID == 3125488) {
+            //     fprintf(stderr, "[M::%s::] ii::[%ld,%ld)(%c), is_srt::%ld, chn_pen_gap::%f, chn_pen_skip::%f, bw_rate::%f\n", __func__, l, k, "+-"[a[l].strand], is_srt, chn_pen_gap, chn_pen_skip, bw_rate);
+            // }
+            if(is_srt) {
+                plus0 = 0; msc0 = msc_i0 = INT32_MIN; movl0 = INT32_MAX; ddt = 0;
+
+
+                p[l] = -1; f[l] = a[l].cnt&(0xffu);
+                if(f[l] >= msc0) {msc0 = f[l]; msc_i0 = l;}///difference
+                if(f[l] < plus0) plus0 = f[l];
+
+
+                for (z = l + 1; z < k; z++) {
+                    ///roughly same to comput_sc_ch(&a[z], &a[z-1])
+                    ai = &a[z]; aj = &a[z-1];
+                    dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
+                    if(dq <= 0) break;
+                    dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
+                    if(dr <= 0) break;
+                    dd = dr > dq? dr - dq : dq - dr;//gap
+                    // if(a_n && a[0].readID == 3125488) {
+                    //     fprintf(stderr, "%ld,", dd);
+                    // }
+                    if((dd > 16) && (dd > cal_bw(&(a[z]), &(a[z-1]), bw_rate, xl, yl))) break;
+                    dg = dr < dq? dr : dq;//len
+                    q_span = ai->cnt&(0xffu); 
+                    sc = q_span < dg? q_span : dg;
+                    sc = normal_w(sc, ((int32_t)(ai->cnt>>8)));
+                    if (dd || (dg > q_span && dg > 0)) {
+                        lin_pen = (chn_pen_gap*(double)dd);
+                        a_pen = ((double)(sc))*((((double)dd)/((double)dg))/bw_rate);
+                        ///for long gap
+                        // if(lin_pen > a_pen) lin_pen = a_pen;
+                        if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
+                        else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
+                        lin_pen += (chn_pen_skip*(double)dg);
+                        sc -= (int32_t)lin_pen;
+                    }
+
+                    sc += f[z-1]; csc = a[z].cnt&(0xffu); if(sc < csc) break;
+                    p[z] = z - 1; f[z] = sc; ddt += dd; 
+
+                    if(f[z] >= msc0) {msc0 = f[z]; msc_i0 = z;}///difference
+                    if(f[z] < plus0) plus0 = f[z];
+                }
+
+                // if(a_n && a[0].readID == 3125488) {
+                //     fprintf(stderr, "\n");
+                //     fprintf(stderr, "[M::%s::] msc0::%ld, msc_i0::%ld, (%c)\n", __func__, msc0, msc_i0, "+-"[a[l].strand]);
+                // }
+                if((z >= k) && (msc_i0 == (k - 1))) {
+                    if((k - l >= 2) && (ddt > 16) && (ddt > cal_bw(&(a[k-1]), &(a[l]), bw_rate, xl, yl))) msc_i0 = INT32_MIN;
+                    if(msc_i0 == (k - 1)) {
+                        if(msc0 >= (*msc)) {
+                            movl0 = get_chainLen(a[msc_i0].self_offset, a[msc_i0].self_offset, xl, a[msc_i0].offset, a[msc_i0].offset, yl);
+                            if(msc0 > (*msc) || movl0 < (*movl)) {
+                                *msc = msc0; *msc_i = msc_i0; *movl = movl0;
+                            }
+                        }
+                        if(plus0 < (*plus)) *plus = plus0;
+                        if((*ei) > k) {
+                            (*si) = k;
+                        } else {
+                            (*ei) = l;
+                        }
+                    }
+                }
+            }
+            l = k; is_srt = 1;
+        } else {
+            if((a[k].self_offset <= a[k-1].self_offset) || (a[k].offset <= a[k-1].offset)) is_srt = 0;
+            t[k-1] = 0; ii[k-1] = 0;
+        }
+    }
+}
+
+
+uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx, 
+              Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter, 
+              int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, 
+              uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be, 
+              int64_t gen_cigar, int64_t mcopy_num, double mcopy_rate, int64_t mcopy_khit_cutoff, 
+              int64_t khit_n)
+{
+    if(a_n <= 0) return 0;
+    int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, max_ii, ovl, movl, plus = 0, min_sc, ch_n, si, ei; 
+    int32_t *f, max, tmp, *ii; int64_t i, k, j, cL = 0; k_mer_hit* a; k_mer_hit* des; k_mer_hit *swap; overlap_region *z;
+    resize_Chain_Data(dp, a_n, NULL); ch_n = 1; // int64_t bw; bw = ((xl < yl)?xl:yl); bw *= bw_rate;
+    t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ;
+    
+    a = cl->list + a_idx; des = cl->list + des_idx;
+    // if(a_n && (a[0].readID == 27105 || a[0].readID == 7603)) {///r833
+    //     fprintf(stderr, "---[M::%s::rid->%u::%c]\ta_n::%ld\n", 
+    //                     __func__, a[0].readID, "+-"[a[0].strand], a_n);
+    // }
+    if(quick_check) {
+        quick_ck_lchain(a, a_n, xl, yl, chn_pen_gap, chn_pen_skip, bw_rate, p, t, f, ii, &plus, &msc, &msc_i, &movl, &si, &ei);
+    } else {
+        msc = msc_i = INT32_MIN; movl = INT32_MAX; plus = 0; si = 0; ei = a_n;
+        memset(t, 0, (a_n*sizeof((*t))));
+    }   
+    // if(a_n && a[0].readID == 4412344) {
+    //     fprintf(stderr, "[M::%s::] si::%ld, ei::%ld, a_n::%ld\n", __func__, si, ei, a_n);
+    // }
+    for (i = st = si, max_ii = -1; i < ei; ++i) {
+        max_f = a[i].cnt&(0xffu); 
+        n_skip = 0; max_j = end_j = -1;
+        if ((i-st) > max_iter) st = i-max_iter;
+        while (a[i].strand != a[st].strand) ++st;
+
+        for (j = i - 1; j >= st; --j) {
+            sc = comput_sc_ch_ec(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
+            if (sc == INT32_MIN) continue;
+            sc += f[j];
+            if (sc > max_f) {
+                max_f = sc, max_j = j;
+                if (n_skip > 0) --n_skip;
+            } else if (t[j] == (int32_t)i) {
+                if (++n_skip > max_skip)
+                    break;
+            }
+            if (p[j] >= 0) t[p[j]] = i;
+        }
+        end_j = j;
+
+        if ((max_ii<0) || (a[i].self_offset>a[max_ii].self_offset+max_dis) || (a[i].strand!=a[max_ii].strand)) {
+            max = INT32_MIN; max_ii = -1;
+            for (j=i-1; (j>=st) && (a[i].self_offset<=max_dis+a[j].self_offset)&&(a[i].strand==a[j].strand); --j) {
+                if (max < f[j]) {
+                    max = f[j], max_ii = j;
+                }
+            }
+        }
+
+        if ((max_ii >= 0) && (max_ii < end_j) && (a[i].strand == a[max_ii].strand)) {///just have a try with a[i]<->a[max_ii]
+            tmp = comput_sc_ch_ec(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
+            if (tmp != INT32_MIN && max_f < tmp + f[max_ii])
+                max_f = tmp + f[max_ii], max_j = max_ii;
+        }
+        f[i] = max_f; p[i] = max_j;
+        if ((max_ii < 0) || ((a[i].self_offset<=max_dis+a[max_ii].self_offset)&&(a[i].strand==a[max_ii].strand)&&(f[max_ii]<f[i]))) {
+            max_ii = i;
+        }
+        if(f[i] >= msc) {
+            ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl);
+            if(f[i] > msc || ovl < movl) {
+                msc = f[i]; msc_i = i; movl = ovl;
+            }
+        }
+        if(f[i] < plus) plus = f[i];
+        ii[i] = 0;///for mcopy, not here
+        // if(a_n && (a[0].readID == 27105 || a[0].readID == 7603)) {///r833
+        //     fprintf(stderr, "i::%ld[M::%s::rid->%u::%c] q::%u, t::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n", 
+        //                     i, __func__, a[i].readID, "+-"[a[i].strand], 
+        //                     a[i].self_offset, a[i].offset, st, max_ii, f[i], p[i], msc_i, msc, movl);
+        // }
+    }
+
+    for (i = msc_i, cL = 0; i >= 0; i = p[i]) { ii[i] = 1; t[cL++] = i;}///label the best chain
+
+    if(mcopy_num > 1) {
+        // if(a[0].readID == 4412344) {
+        //     fprintf(stderr, "[M::%s::] msc::%ld, cL::%ld\n", __func__, msc, cL);
+        // }
+        if(cL >= mcopy_khit_cutoff) {///if there are too few k-mers, disable mcopy
+            msc -= plus; min_sc = msc*mcopy_rate/**0.2**/; ii[msc_i] = 0;
+            for (i = ch_n = 0; i < a_n; ++i) {///make all f[] positive
+                f[i] -= plus; if(i >= ch_n) t[i] = 0;
+                if((!(ii[i])) && (f[i] >= min_sc)) {///!(ii[i]): skip the best chain
+                    t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++;
+                }
+            }
+            // if(a[0].readID == 4412344) {
+            //     fprintf(stderr, "[M::%s::] msc::%ld, min_sc::%ld, cL::%ld, ch_n::%ld, mcopy_num::%ld\n", __func__, msc, min_sc, cL, ch_n, mcopy_num);
+            // }
+            if(ch_n > 1) {
+                int64_t n_v, n_v0, ni, n_u, n_u0 = res->length; 
+                radix_sort_hc64i(t, t + ch_n);
+                for (k = ch_n-1, n_v = n_u = 0; k >= 0 && n_u < mcopy_num; --k) {
+                    n_v0 = n_v;
+                    for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) {
+                        ii[n_v++] = i; t[i] |= 1; i = p[i];
+                    }
+                    if(n_v0 == n_v) continue;
+                    sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i]));
+                    // if(a[0].readID == 4412344) {
+                    //     fprintf(stderr, "+[M::%s::] sc::%ld, n_a::%ld\n", __func__, sc, n_v-n_v0);
+                    // }
+                    if(sc >= min_sc) {
+                        kv_pushp_ol(overlap_region, (*res), &z);
+                        push_ovlp_chain_qgen(z, xid, xl, yl, sc+plus, &(a[ii[n_v-1]]), &(a[ii[n_v0]]));
+                        // if(a[0].readID == 4412344) {
+                        //     fprintf(stderr, "-[M::%s::] sc::%ld, n_a::%ld, q::[%u,%u), t::[%u,%u), %c\n", __func__, sc, n_v-n_v0, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, "+-"[z->y_pos_strand]);
+                        // }
+                        ///mcopy_khit_cutoff <= 1: disable the mcopy_khit_cutoff filtering, for the realignment
+                        // if((mcopy_khit_cutoff <= 1) || ((z->x_pos_e+1-z->x_pos_s) <= (movl<<2))) {
+                        if((!n_u) || (n_v - n_v0 > 1)) {
+                            z->align_length = n_v-n_v0; z->x_id = n_v0;
+                            n_u++;
+                        } else {///non-best is tiny
+                            res->length--; n_v = n_v0;
+                        }
+                    } else {
+                        n_v = n_v0;
+                    }
+                }
+
+                // if(n_u > 1) ks_introsort_or_sss(n_u, res->list + n_u0); 
+                // res->length = n_u0 + filter_non_ovlp_xchains(res->list + n_u0, n_u, &n_v);
+                n_u = res->length;
+                if(n_u > n_u0 + 1) {
+                    kv_resize_cl(k_mer_hit, (*cl), (n_v+cl->length));
+                    a = cl->list + a_idx; des = cl->list + des_idx; swap = cl->list + cl->length; 
+                    for (k = n_u0, i = n_v0 = n_v = 0; k < n_u; k++) {
+                        z = &(res->list[k]);
+                        z->non_homopolymer_errors = des_idx + i;
+                        n_v0 = z->x_id; ni = z->align_length;
+                        for (j = 0; j < ni; j++, i++) {
+                            ///k0 + (ni - j - 1)
+                            swap[i] = a[ii[n_v0 + (ni- j - 1)]]; 
+                            swap[i].readID = k; 
+                        }
+                        z->x_id = xid; 
+                        if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, swap+i-ni, ni);
+                        if(!khit_n) z->align_length = 0; 
+                    }
+                    memcpy(des, swap, i*sizeof((*swap))); //assert(i == ch_n);
+                    
+                    // fprintf(stderr, "[M::%s::msc->%ld] msc_k_hits::%u, cL::%ld, min_sc::%ld, best_sc::%ld, n_u0_sc::%d, mcopy_rate::%f, # chains::%ld\n", 
+                    // __func__, msc, res->list[n_u0].align_length, cL, min_sc, msc+plus, res->list[n_u0].shared_seed,
+                    // mcopy_rate, n_u-n_u0);
+                } else if(n_u == n_u0 + 1) {
+                    z = &(res->list[n_u0]); k = n_u0; i = 0;
+                    z->non_homopolymer_errors = des_idx + i;
+                    n_v0 = z->x_id; ni = z->align_length;
+                    for (j = 0; j < ni; j++, i++) {
+                        ///k0 + (ni - j - 1)
+                        des[i] = a[ii[n_v0 + (ni- j - 1)]]; 
+                        des[i].readID = k; 
+                    }
+                    z->x_id = xid; 
+                    if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des+i-ni, ni);
+                    if(!khit_n) z->align_length = 0; 
+                }
+                return i;
+            } else {
+                msc += plus; i = msc_i; cL = 0; 
+                while (i >= 0) {t[cL++] = i; i = p[i];}
+            }
+        }
+    }
+
+
+
+    ///a[] has been sorted by self_offset
+    // i = msc_i; cL = 0; 
+    // while (i >= 0) {t[cL++] = i; i = p[i];}
+    kv_pushp_ol(overlap_region, (*res), &z);
+    push_ovlp_chain_qgen(z, xid, xl, yl, msc, &(a[t[cL-1]]), &(a[t[0]]));
+    for (i = 0; i < cL; i++) {des[i] = a[t[cL-i-1]]; des[i].readID = res->length-1;}
+    z->non_homopolymer_errors = des_idx;
+    if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des, cL);
+    if(khit_n) z->align_length = cL;
+    return cL;
+}
 
 
 #define rev_khit(an, xl, yl) do {	\
@@ -2308,6 +2603,57 @@ uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp
     return cL;
 }
 
+uint64_t lchain_simple0(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter)
+{
+    if(a_n <= 0) return 0;
+    int64_t *p, *t, max_f, n_skip, st, max_j, sc, msc, msc_i; 
+    int32_t *f; int64_t i, j, cL = 0; 
+    resize_Chain_Data(dp, a_n, NULL);
+    t = dp->tmp; f = dp->score; p = dp->pre; msc = msc_i = -1;
+    
+    memset(t, 0, (a_n*sizeof((*t))));
+    f[0]=a[0].cnt; p[0]=-1; msc = f[0]; msc_i = 0;
+
+    for (i = 1, st = 0; i < a_n; ++i) {
+        max_f = INT32_MIN; n_skip = 0; max_j = -1;
+        if ((i-st) > max_iter) st = i-max_iter;
+        ///[st, i-2]
+        for (j=i-1; j >= st; --j) {
+            if((a[i].self_offset > a[j].self_offset)&&(a[i].offset > a[j].offset)) {
+                sc = f[j]+a[i].cnt;
+                if (sc > max_f) {
+                    max_f = sc, max_j = j;
+                    if (n_skip > 0) --n_skip;
+                } else if (t[j] == (int32_t)i) {
+                    if (++n_skip > max_skip)
+                        break;
+                }
+                if (p[j] >= 0) t[p[j]] = i;
+            }
+        }
+        f[i] = max_f; p[i] = max_j;
+        if(f[i] > msc) {
+            msc = f[i]; msc_i = i;
+        }
+    }
+
+    ///a[] has been sorted by self_offset
+    i = msc_i; 
+    cL = 0; 
+    while (i >= 0) {
+        t[cL++] = i; i = p[i];
+    }
+
+    n_skip = cL>>1;
+    for (i = 0; i < n_skip; i++) {
+        msc_i = t[i]; t[i] = t[cL-i-1]; t[cL-i-1] = msc_i;
+    }
+    if(des) {
+        for (i = 0; i < cL; i++) des[i] = a[t[i]];
+    }
+    return cL;
+}
+
 inline int64_t hit_long_gap(k_mer_hit *a, k_mer_hit *b, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw)
 {
     int64_t dq, dr, dd, dm;

Hash_Table.h

@@ -8,10 +8,17 @@
 
 #define WINDOW 375
 #define WINDOW_BOUNDARY 375
+#define WINDOW_HC 775
+///ONT high error
+// #define WINDOW_OHC 475
+#define WINDOW_OHC 375
+#define WINDOW_HC_FAST 512
 ///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected
 #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
 #define THRESHOLD  15
 #define OVERLAP_THRESHOLD_HIFI_FILTER 0.9
+#define OVERLAP_THRESHOLD_HIFI_FF_FILTER 0.6
+#define OVERLAP_THRESHOLD_HIFI_FF_DE_FILTER 0.5
 #define OVERLAP_THRESHOLD_NOSI_FILTER 0.7
 #define OVERLAP_THRESHOLD_FILTER_HPC 0.75
 #define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3
@@ -229,6 +236,9 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s
                 int64_t left_fix, int64_t right_fix);
 uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, 
                                                         int64_t max_skip, int64_t max_iter);
+uint64_t lchain_simple0(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter);
+void push_ovlp_chain_qgen(overlap_region* o, uint32_t xid, int64_t xl, int64_t yl, int64_t sc, k_mer_hit *beg, k_mer_hit *end);
+
 uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx, 
               Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter, 
               int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, 
@@ -236,4 +246,20 @@ uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64
               int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff, 
               int64_t khit_n);
 
+uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx, 
+              Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter, 
+              int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, 
+              uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be, 
+              int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff, 
+              int64_t khit_n);
+
+#define kv_pushp_ol(type, v, p) do {									\
+		if ((v).length == (v).size) {										\
+			(v).list = (type*)realloc((v).list, sizeof(type)*((v).size?((v).size<<1):(2)));	\
+            memset((v).list+(v).size, 0, sizeof(overlap_region)*(((v).size?((v).size<<1):2)-(v).size));\
+            (v).size = (v).size?((v).size<<1):(2);							\
+		}															\
+		*(p) = &((v).list[(v).length++]); \
+	} while (0)
+
 #endif

Levenshtein_avx2.cpp

@@ -0,0 +1,162 @@
+#include "Levenshtein_distance.h"
+#include <immintrin.h>
+
+
+#define init_simd_ed4(PSA, PNA, THRE, ABS_DIAG, R_ERR, R_PE, SI, TN, CUT, BD, I, MM, PEQ_MM, LZ, IBD) {\
+	(R_ERR)[(SI)] = INT32_MAX; (R_PE)[(SI)] = -1; (IBD)[(SI)] = ((THRE)<<1) - (ABS_DIAG)[(SI)];\
+	if(((PNA)[(SI)] <= (TN) + (CUT)) && ((TN) <= (PNA)[(SI)] + (CUT))) {\
+		(BD) = (((THRE)<<1)+1)-(ABS_DIAG)[(SI)]; (BD) = (((BD)<=(PNA)[(SI)])?(BD):(PNA)[(SI)]); (LZ) |= (((int32_t)1u) << (SI));\
+		for ((I) = 0, (MM) = (((Word)1)<<((ABS_DIAG)[(SI)])); (I) < (BD); (I)++) {\
+			(PEQ_MM)[seq_nt4_table[(uint8_t)(PSA)[(SI)][(I)]]][(SI)] |= (MM); (MM) <<= 1;\
+		}\
+	}\
+}
+
+#define ed_core_64x4(PEQz, VPz, VNz, Xz, D0z, HNz, HPz) {	\
+	/**(X) = (Peq)|(VN);**/\
+	(Xz) = _mm256_or_si256((PEQz), (VNz)); \
+	/**(D0) = (((VP) + ((X)&(VP))) ^ (VP)) | (X);**/\
+	(D0z) = _mm256_or_si256(_mm256_xor_si256(_mm256_add_epi64((VPz), _mm256_and_si256((Xz), (VPz))), (VPz)), (Xz)); \
+	/**(HN) = (VP)&(D0);**/\
+	(HNz) = _mm256_and_si256((VPz), (D0z)); \
+	/**(HP) = (VN) | ~((VP) | (D0));**/\
+	(HPz) = _mm256_or_si256((VNz), _mm256_andnot_si256(_mm256_or_si256((VPz), (D0z)), _mm256_set1_epi64x(-1))); \
+	/**(X) = (D0) >> 1;**/\
+	(Xz) = _mm256_srli_epi64((D0z), 1); \
+	/**(VN) = (X)&(HP);**/\
+	(VNz) = _mm256_and_si256((Xz), (HPz)); \
+	/**(VP) = (HN) | ~((X) | (HP));**/\
+	(VPz) = _mm256_or_si256((HNz), _mm256_andnot_si256(_mm256_or_si256((Xz), (HPz)), _mm256_set1_epi64x(-1))); \
+}
+
+#define ed_core_upx4(PEQz, PSA, PNA, IBD, HT, CC, MMK, SI) {	\
+	if((HT) & (((int32_t)1u) << (SI))) {\
+		(IBD)[(SI)]++;\
+		if((IBD)[(SI)] < (PNA)[(SI)]) {\
+			(CC) = seq_nt4_table[(uint8_t)(PSA)[(SI)][(IBD)[(SI)]]];\
+			if((CC) < 4) (PEQz)[(CC)] = _mm256_or_si256((PEQz)[(CC)], (MMK)[(SI)]);\
+		}\
+	}\
+}
+
+#define ed_tail_upx4(HT, SI, ST, AI, PNA, ABS_DIAG, K, ERR_MM, VP_MM, VN_MM, THRE, R_ERR, R_PE, BD, I) {\
+	if((HT) & (((int32_t)1u) << (SI))) {\
+		(ST)[(SI)] -= (ABS_DIAG)[(SI)]; (AI)[(SI)] += (PNA)[(SI)] + (ABS_DIAG)[(SI)];\
+		for ((K)[(SI)] = 0; (ST)[(SI)] < 0 && (K)[(SI)] < (AI)[(SI)]; (K)[(SI)]++, (ST)[(SI)]++) {\
+			(ERR_MM)[(SI)] += ((VP_MM)[(SI)]&(1ULL)); (VP_MM)[(SI)]>>=1;\
+			(ERR_MM)[(SI)] -= ((VN_MM)[(SI)]&(1ULL)); (VN_MM)[(SI)]>>=1;\
+		}\
+		if (((ERR_MM)[(SI)] <= (THRE)) && ((ERR_MM)[(SI)] <= (R_ERR)[(SI)])) {\
+			(R_ERR)[(SI)] = (ERR_MM)[(SI)]; (R_PE)[(SI)] = (ST)[(SI)];\
+		}\
+		(ST)[(SI)] -= (K)[(SI)]; (BD)++; (I) = (SI);\
+	}\
+}
+
+#define ED_TAIL_LANE(K) do { \
+    if ((ht) & (((int32_t)1u) << (K))) { \
+        st = tn - 1 - abs_diag_a[(K)]; \
+        ai = pna[(K)] - tn + abs_diag_a[(K)]; \
+        for (i = 0, uge = INT64_MAX; st < 0 && i < ai; i++, st++) { \
+            err_mm[(K)] += ((VP_mm[(K)] >> i) & 1ULL); \
+            err_mm[(K)] -= ((VN_mm[(K)] >> i) & 1ULL); \
+        } \
+        if ((err_mm[(K)] <= thre) && (err_mm[(K)] <= r_err[(K)])) { \
+            r_err[(K)] = err_mm[(K)]; \
+            r_pe[(K)] = st; \
+        } \
+        st -= i; \
+        while (i < ai) { \
+            err_mm[(K)] += ((VP_mm[(K)] >> i) & 1ULL); \
+            err_mm[(K)] -= ((VN_mm[(K)] >> i) & 1ULL); \
+            ++i; \
+            if ((err_mm[(K)] <= thre) && (err_mm[(K)] <= r_err[(K)])) { \
+                r_err[(K)] = err_mm[(K)]; \
+                r_pe[(K)] = st + i; \
+            } \
+            if (i == thre) uge = err_mm[(K)]; \
+        } \
+        if ((uge <= thre) && (uge == r_err[(K)])) r_pe[(K)] = st + thre; \
+    } \
+} while (0)
+
+void ed_band_cal_semi_64_w_absent_diag_avx4(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe)
+{
+	// r_err[0] = r_err[1] = r_err[2] = r_err[3] = r_err[4] = r_err[5] = r_err[6] = r_err[7] = thre+1;
+	// r_pe[0] = r_pe[1] = r_pe[2] = r_pe[3] = r_pe[4] = r_pe[5] = r_pe[6] = r_pe[7] = -1;   
+	
+	Word mm, Peq_mm[5][AVX_GS2] = {{0}}, *VN_mm = NULL, *VP_mm = NULL, c = 0; __m256i Peq[5], VP, VN, X, D0, HN, HP, lone, E, C, mmk[AVX_GS2]; 
+	int32_t lz = 0, ht = (((int32_t)1u)<<AVX_GS2)-1; int32_t bd, ibd[AVX_GS2], i, last_high = (thre<<1), tn0 = tn - 1, cut = thre+last_high;
+
+	lone = _mm256_set1_epi64x(1);
+	VP = _mm256_setzero_si256();
+
+	VN_mm = Peq_mm[0];
+	VN_mm[0] = (((Word)1)<<(abs_diag_a[0]))-1; VN_mm[1] = (((Word)1)<<(abs_diag_a[1]))-1; VN_mm[2] = (((Word)1)<<(abs_diag_a[2]))-1; VN_mm[3] = (((Word)1)<<(abs_diag_a[3]))-1; 
+	VN = _mm256_loadu_si256((const __m256i *)VN_mm);
+
+	VN_mm[0] = abs_diag_a[0]; VN_mm[1] = abs_diag_a[1]; VN_mm[2] = abs_diag_a[2]; VN_mm[3] = abs_diag_a[3]; 
+	E = _mm256_loadu_si256((const __m256i *)VN_mm);
+
+	memset(VN_mm, 0, (sizeof((*VN_mm))*AVX_GS2)); VN_mm = NULL;///reset
+
+	init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 0, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 1, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 2, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 3, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+
+	ht &= lz;
+	if(ht == 0) return;
+
+	Peq[0] = _mm256_loadu_si256((const __m256i *)Peq_mm[0]); 
+	Peq[1] = _mm256_loadu_si256((const __m256i *)Peq_mm[1]); 
+	Peq[2] = _mm256_loadu_si256((const __m256i *)Peq_mm[2]); 
+	Peq[3] = _mm256_loadu_si256((const __m256i *)Peq_mm[3]); 
+	Peq[4] = _mm256_setzero_si256();
+	
+	C = _mm256_set1_epi64x(cut + 1);///_mm512_set1_epi64x(cut);
+
+	mm = ((Word)1 << (thre<<1));///for the incoming char/last char**
+	mmk[0] = _mm256_set_epi64x(0, 0, 0, mm);
+	mmk[1] = _mm256_set_epi64x(0, 0, mm, 0);
+	mmk[2] = _mm256_set_epi64x(0, mm, 0, 0);
+	mmk[3] = _mm256_set_epi64x(mm, 0, 0, 0);
+
+	i = 0; 
+
+	while (i < tn0) {
+		ed_core_64x4(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
+		E = _mm256_add_epi64(_mm256_xor_si256(lone, _mm256_and_si256(D0, lone)), E);
+		// ht = _mm512_cmple_epi64_mask(E, C);
+		ht =_mm256_movemask_pd(_mm256_castsi256_pd(_mm256_cmpgt_epi64(C, E)));
+		ht &= lz;
+		if(ht == 0) return;
+
+		Peq[0] = _mm256_srli_epi64(Peq[0], 1);
+		Peq[1] = _mm256_srli_epi64(Peq[1], 1);
+		Peq[2] = _mm256_srli_epi64(Peq[2], 1);
+		Peq[3] = _mm256_srli_epi64(Peq[3], 1);
+		i++; ///c = 4;
+
+		ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 0);
+		ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 1);
+		ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 2);
+		ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 3);
+	}
+
+	ed_core_64x4(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
+	E = _mm256_add_epi64(_mm256_xor_si256(lone, _mm256_and_si256(D0, lone)), E);
+	// ht = _mm512_cmple_epi64_mask(E, C);
+	ht =_mm256_movemask_pd(_mm256_castsi256_pd(_mm256_cmpgt_epi64(C, E)));
+	ht &= lz;
+	if(ht == 0) return;
+
+	int32_t st, ai; int64_t err_mm[AVX_GS2], uge;
+	VN_mm = Peq_mm[0]; VP_mm = Peq_mm[1];
+	_mm256_storeu_si256((__m256i *)VN_mm, VN); _mm256_storeu_si256((__m256i *)VP_mm, VP); _mm256_storeu_si256((__m256i *)err_mm, E);
+
+	ED_TAIL_LANE(0);
+	ED_TAIL_LANE(1);
+	ED_TAIL_LANE(2);
+	ED_TAIL_LANE(3);
+}

Levenshtein_avx512.cpp

@@ -0,0 +1,276 @@
+#include "Levenshtein_distance.h"
+#include <immintrin.h>
+
+#define init_simd_ed(PSA, PNA, THRE, ABS_DIAG, R_ERR, R_PE, SI, TN, CUT, BD, I, MM, PEQ_MM, LZ, IBD) {\
+	(R_ERR)[(SI)] = INT32_MAX; (R_PE)[(SI)] = -1; (IBD)[(SI)] = ((THRE)<<1) - (ABS_DIAG)[(SI)];\
+	if(((PNA)[(SI)] <= (TN) + (CUT)) && ((TN) <= (PNA)[(SI)] + (CUT))) {\
+		(BD) = (((THRE)<<1)+1)-(ABS_DIAG)[(SI)]; (BD) = (((BD)<=(PNA)[(SI)])?(BD):(PNA)[(SI)]); (LZ) |= (((__mmask8)1u) << (SI));\
+		for ((I) = 0, (MM) = (((Word)1)<<((ABS_DIAG)[(SI)])); (I) < (BD); (I)++) {\
+			(PEQ_MM)[seq_nt4_table[(uint8_t)(PSA)[(SI)][(I)]]][(SI)] |= (MM); (MM) <<= 1;\
+		}\
+	}\
+}
+
+#define ed_core_64x8(PEQz, VPz, VNz, Xz, D0z, HNz, HPz) {	\
+	/**(X) = (Peq)|(VN);**/\
+	(Xz) = _mm512_or_si512((PEQz), (VNz));\
+	/**(D0) = (((VP) + ((X)&(VP))) ^ (VP)) | (X);**/\
+	(D0z) = _mm512_or_si512(_mm512_xor_si512(_mm512_add_epi64((VPz), _mm512_and_si512((Xz), (VPz))), (VPz)), (Xz));\
+	/**(HN) = (VP)&(D0);**/\
+	(HNz) = _mm512_and_si512((VPz), (D0z));\
+	/**(HP) = (VN) | ~((VP) | (D0));**/\
+	(HPz) = _mm512_or_si512((VNz), _mm512_andnot_si512(_mm512_or_si512((VPz), (D0z)), _mm512_set1_epi64(-1)));\
+	/**(X) = (D0) >> 1;**/\
+	(Xz)  = _mm512_srli_epi64((D0z), 1);\
+	/**(VN) = (X)&(HP);**/\
+	(VNz) = _mm512_and_si512((Xz), (HPz));\
+	/**(VP) = (HN) | ~((X) | (HP));**/\
+	(VPz) = _mm512_or_si512((HNz), _mm512_andnot_si512(_mm512_or_si512((Xz), (HPz)), _mm512_set1_epi64(-1)));\
+}
+
+#define ed_core_upx8(PEQz, PSA, PNA, IBD, HT, CC, MMK, SI) {	\
+	if((HT) & (((__mmask8)1u) << (SI))) {\
+		(IBD)[(SI)]++;\
+		if((IBD)[(SI)] < (PNA)[(SI)]) {\
+			(CC) = seq_nt4_table[(uint8_t)(PSA)[(SI)][(IBD)[(SI)]]];\
+			if((CC) < 4) (PEQz)[(CC)] = _mm512_or_si512((PEQz)[(CC)], (MMK)[(SI)]);\
+		}\
+	}\
+}
+
+#define ed_tail_upx8(HT, SI, ST, AI, PNA, ABS_DIAG, K, ERR_MM, VP_MM, VN_MM, THRE, R_ERR, R_PE, BD, I) {\
+	if((HT) & (((__mmask8)1u) << (SI))) {\
+		(ST)[(SI)] -= (ABS_DIAG)[(SI)]; (AI)[(SI)] += (PNA)[(SI)] + (ABS_DIAG)[(SI)];\
+		for ((K)[(SI)] = 0; (ST)[(SI)] < 0 && (K)[(SI)] < (AI)[(SI)]; (K)[(SI)]++, (ST)[(SI)]++) {\
+			(ERR_MM)[(SI)] += ((VP_MM)[(SI)]&(1ULL)); (VP_MM)[(SI)]>>=1;\
+			(ERR_MM)[(SI)] -= ((VN_MM)[(SI)]&(1ULL)); (VN_MM)[(SI)]>>=1;\
+		}\
+		if (((ERR_MM)[(SI)] <= (THRE)) && ((ERR_MM)[(SI)] <= (R_ERR)[(SI)])) {\
+			(R_ERR)[(SI)] = (ERR_MM)[(SI)]; (R_PE)[(SI)] = (ST)[(SI)];\
+		}\
+		(ST)[(SI)] -= (K)[(SI)]; (BD)++; (I) = (SI);\
+	}\
+}
+
+#define ed_tail_ck8(MBEST, SI, R_PE, ST, K, THRE, UGE_MM, ERR_MM, AI, HT) {\
+	(K)[(SI)]++;\
+	if((MBEST) & (((__mmask8)1u) << (SI))) {\
+		(R_PE)[(SI)] = (ST)[(SI)] + (K)[(SI)];\
+	}\
+	if((K)[(SI)] >= (AI)[(SI)]) (HT) &= ~(((__mmask8)1u) << (SI));\
+	if((K)[(SI)] == (THRE)) (UGE_MM)[(SI)] = (ERR_MM)[(SI)];\
+}
+
+
+void ed_band_cal_semi_64_w_absent_diag_avx8(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe)
+{
+	// r_err[0] = r_err[1] = r_err[2] = r_err[3] = r_err[4] = r_err[5] = r_err[6] = r_err[7] = thre+1;
+	// r_pe[0] = r_pe[1] = r_pe[2] = r_pe[3] = r_pe[4] = r_pe[5] = r_pe[6] = r_pe[7] = -1;
+	/**   
+	ed_band_cal_semi_64_w_absent_diag_avx4(psa, pna, tstr, tn, thre, abs_diag_a, r_err, r_pe);
+	ed_band_cal_semi_64_w_absent_diag_avx4(psa + 4, pna + 4, tstr, tn, thre, abs_diag_a + 4, r_err + 4, r_pe + 4);
+	return;
+	**/
+
+	
+	Word mm, Peq_mm[5][AVX_GS] = {{0}}, *VN_mm = NULL, *VP_mm = NULL, c = 0; __m512i Peq[5], VP, VN, X, D0, HN, HP, lone, E, C, bestE, bestPE, curPE, cutPE, threPE, ugE, mmk[AVX_GS]; 
+	__mmask8 lz = ((__mmask8)0u), ht = (((__mmask8)1u)<<AVX_GS)-1, mtf, mt, mbest; int32_t bd, ibd[AVX_GS], i, last_high = (thre<<1), tn0 = tn - 1, cut = thre+last_high;
+
+	lone = _mm512_set1_epi64(1);
+	VP = _mm512_setzero_si512();
+
+	VN_mm = Peq_mm[0];
+	VN_mm[0] = (((Word)1)<<(abs_diag_a[0]))-1; VN_mm[1] = (((Word)1)<<(abs_diag_a[1]))-1; VN_mm[2] = (((Word)1)<<(abs_diag_a[2]))-1; VN_mm[3] = (((Word)1)<<(abs_diag_a[3]))-1; 
+	VN_mm[4] = (((Word)1)<<(abs_diag_a[4]))-1; VN_mm[5] = (((Word)1)<<(abs_diag_a[5]))-1; VN_mm[6] = (((Word)1)<<(abs_diag_a[6]))-1; VN_mm[7] = (((Word)1)<<(abs_diag_a[7]))-1; 
+	VN = _mm512_loadu_si512(VN_mm);
+
+	VN_mm[0] = abs_diag_a[0]; VN_mm[1] = abs_diag_a[1]; VN_mm[2] = abs_diag_a[2]; VN_mm[3] = abs_diag_a[3]; 
+	VN_mm[4] = abs_diag_a[4]; VN_mm[5] = abs_diag_a[5]; VN_mm[6] = abs_diag_a[6]; VN_mm[7] = abs_diag_a[7];
+	E = _mm512_loadu_si512(VN_mm);
+
+	memset(VN_mm, 0, (sizeof((*VN_mm))*AVX_GS)); VN_mm = NULL;///reset
+
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 0, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 1, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 2, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 3, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 4, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 5, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 6, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+	init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 7, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
+
+	ht &= lz;
+	if(ht == 0) return;
+
+	Peq[0] = _mm512_loadu_si512(Peq_mm[0]); 
+	Peq[1] = _mm512_loadu_si512(Peq_mm[1]); 
+	Peq[2] = _mm512_loadu_si512(Peq_mm[2]); 
+	Peq[3] = _mm512_loadu_si512(Peq_mm[3]); 
+	Peq[4] = _mm512_setzero_si512();
+
+	C = _mm512_set1_epi64(cut);
+
+	mm = ((Word)1 << (thre<<1));///for the incoming char/last char**
+	mmk[0] = _mm512_mask_set1_epi64(VP, 1, mm);
+	mmk[1] = _mm512_mask_set1_epi64(VP, 2, mm);
+	mmk[2] = _mm512_mask_set1_epi64(VP, 4, mm);
+	mmk[3] = _mm512_mask_set1_epi64(VP, 8, mm);
+	mmk[4] = _mm512_mask_set1_epi64(VP, 16, mm);
+	mmk[5] = _mm512_mask_set1_epi64(VP, 32, mm);
+	mmk[6] = _mm512_mask_set1_epi64(VP, 64, mm);
+	mmk[7] = _mm512_mask_set1_epi64(VP, 128, mm);
+
+	i = 0; 
+
+	while (i < tn0) {
+		ed_core_64x8(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
+		E = _mm512_add_epi64(_mm512_xor_si512(lone, _mm512_and_si512(D0, lone)), E);
+		ht = _mm512_cmple_epi64_mask(E, C);
+		ht &= lz;
+		if(ht == 0) return;
+
+		Peq[0] = _mm512_srli_epi64(Peq[0], 1);
+		Peq[1] = _mm512_srli_epi64(Peq[1], 1);
+		Peq[2] = _mm512_srli_epi64(Peq[2], 1);
+		Peq[3] = _mm512_srli_epi64(Peq[3], 1);
+		i++; ///c = 4;
+
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 0);
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 1);
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 2);
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 3);
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 4);
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 5);
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 6);
+		ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 7);
+	}
+
+	ed_core_64x8(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
+	E = _mm512_add_epi64(_mm512_xor_si512(lone, _mm512_and_si512(D0, lone)), E);
+	ht = _mm512_cmple_epi64_mask(E, C);
+	ht &= lz;
+	if(ht == 0) return;
+
+
+	// site = tn - 1 - abs_diag;/**up bound**/
+	// ai = pn - tn + abs_diag; /**in most cases, ai = (thre<<1)**/
+	int32_t st[AVX_GS] = {tn-1, tn-1, tn-1, tn-1, tn-1, tn-1, tn-1, tn-1}; 
+	int32_t ai[AVX_GS] = {-tn, -tn, -tn, -tn, -tn, -tn, -tn, -tn};
+	int32_t k[AVX_GS] = {0}; i = -1; bd = 0;
+	int64_t err_mm[AVX_GS], uge_mm[AVX_GS] = {INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX}; 
+	VN_mm = Peq_mm[0]; VP_mm = Peq_mm[1];
+	_mm512_storeu_si512(VN_mm, VN); _mm512_storeu_si512(VP_mm, VP); _mm512_storeu_si512(err_mm, E);
+
+	ed_tail_upx8(ht, 0, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+	ed_tail_upx8(ht, 1, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+	ed_tail_upx8(ht, 2, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+	ed_tail_upx8(ht, 3, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+	ed_tail_upx8(ht, 4, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+	ed_tail_upx8(ht, 5, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+	ed_tail_upx8(ht, 6, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+	ed_tail_upx8(ht, 7, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
+
+	if(bd <= 0) return;
+
+	if(bd > 1) {
+		VN = _mm512_loadu_si512(VN_mm); VP = _mm512_loadu_si512(VP_mm); E = _mm512_loadu_si512(err_mm); i = 0;
+
+		bestE = _mm512_loadu_si512(r_err); ///threE = _mm512_set1_epi64(thre);
+		if(k[0] >= ai[0]) ht &= ((__mmask8)(255-1));
+		if(k[1] >= ai[1]) ht &= ((__mmask8)(255-2));
+		if(k[2] >= ai[2]) ht &= ((__mmask8)(255-4));
+		if(k[3] >= ai[3]) ht &= ((__mmask8)(255-8));
+		if(k[4] >= ai[4]) ht &= ((__mmask8)(255-16));
+		if(k[5] >= ai[5]) ht &= ((__mmask8)(255-32));
+		if(k[6] >= ai[6]) ht &= ((__mmask8)(255-64));
+		if(k[7] >= ai[7]) ht &= ((__mmask8)(255-128));
+
+		err_mm[0] = r_pe[0]; err_mm[1] = r_pe[1]; err_mm[2] = r_pe[2]; err_mm[3] = r_pe[3]; 
+		err_mm[4] = r_pe[4]; err_mm[5] = r_pe[5]; err_mm[6] = r_pe[6]; err_mm[7] = r_pe[7]; 
+		bestPE = _mm512_loadu_si512(err_mm);
+		err_mm[0] = st[0] + k[0]; err_mm[1] = st[1] + k[1]; err_mm[2] = st[2] + k[2]; err_mm[3] = st[3] + k[3];
+		err_mm[4] = st[4] + k[4]; err_mm[5] = st[5] + k[5]; err_mm[6] = st[6] + k[6]; err_mm[7] = st[7] + k[7];
+		curPE = _mm512_loadu_si512(err_mm);	
+		err_mm[0] = st[0] + thre; err_mm[1] = st[1] + thre; err_mm[2] = st[2] + thre; err_mm[3] = st[3] + thre;
+		err_mm[4] = st[4] + thre; err_mm[5] = st[5] + thre; err_mm[6] = st[6] + thre; err_mm[7] = st[7] + thre;
+		threPE = _mm512_loadu_si512(err_mm);
+		err_mm[0] = st[0] + ai[0]; err_mm[1] = st[1] + ai[1]; err_mm[2] = st[2] + ai[2]; err_mm[3] = st[3] + ai[3];
+		err_mm[4] = st[4] + ai[4]; err_mm[5] = st[5] + ai[5]; err_mm[6] = st[6] + ai[6]; err_mm[7] = st[7] + ai[7];
+		cutPE = _mm512_loadu_si512(err_mm);
+
+		ugE = _mm512_loadu_si512(uge_mm);
+
+		// mtf = _mm512_cmpge_epi64_mask(curPE, threPE) | ((__mmask8)(~ht));
+		mtf = _mm512_cmpge_epi64_mask(curPE, threPE);
+
+		while ((ht != 0) && ((mtf|((__mmask8)(~ht))) != (__mmask8)255)) {
+			E = _mm512_add_epi64(E, _mm512_and_si512(VP, lone)); VP = _mm512_srli_epi64(VP, 1);
+			E = _mm512_sub_epi64(E, _mm512_and_si512(VN, lone)); VN = _mm512_srli_epi64(VN, 1);
+			// i++;
+			
+			curPE = _mm512_add_epi64(curPE, lone);
+			ht &= _mm512_cmple_epi64_mask(curPE, cutPE);
+			if (ht == 0) break;
+
+			mbest = _mm512_cmple_epi64_mask(E, bestE) & ht;
+			
+			bestE = _mm512_mask_mov_epi64(bestE, mbest, E);
+			bestPE = _mm512_mask_mov_epi64(bestPE, mbest, curPE);
+
+			mt = _mm512_cmpeq_epi64_mask(curPE, threPE);
+			ugE = _mm512_mask_mov_epi64(ugE, mt&ht, E);
+
+			mtf |= mt;
+	
+			// if(mbest && i < thre) _mm512_storeu_si512(err_mm, E);
+
+			// ed_tail_ck8(mbest, 0, r_pe, st, k, thre, uge_mm, err_mm, ai, ht);
+			// ed_tail_ck8(mbest, 1, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); 
+			// ed_tail_ck8(mbest, 2, r_pe, st, k, thre, uge_mm, err_mm, ai, ht);
+			// ed_tail_ck8(mbest, 3, r_pe, st, k, thre, uge_mm, err_mm, ai, ht);
+			// ed_tail_ck8(mbest, 4, r_pe, st, k, thre, uge_mm, err_mm, ai, ht);
+			// ed_tail_ck8(mbest, 5, r_pe, st, k, thre, uge_mm, err_mm, ai, ht);
+			// ed_tail_ck8(mbest, 6, r_pe, st, k, thre, uge_mm, err_mm, ai, ht);
+			// ed_tail_ck8(mbest, 7, r_pe, st, k, thre, uge_mm, err_mm, ai, ht);
+		}
+		
+
+		while (ht != 0) {
+			E = _mm512_add_epi64(E, _mm512_and_si512(VP, lone)); VP = _mm512_srli_epi64(VP, 1);
+			E = _mm512_sub_epi64(E, _mm512_and_si512(VN, lone)); VN = _mm512_srli_epi64(VN, 1);
+			// i++;
+			
+			curPE = _mm512_add_epi64(curPE, lone);
+			ht &= _mm512_cmple_epi64_mask(curPE, cutPE);
+			if (ht == 0) break;
+
+			mbest = _mm512_cmple_epi64_mask(E, bestE) & ht;
+			
+			bestE = _mm512_mask_mov_epi64(bestE, mbest, E);
+			bestPE = _mm512_mask_mov_epi64(bestPE, mbest, curPE);
+		}
+
+		cutPE = _mm512_set1_epi64(thre);
+		ht = _mm512_cmpgt_epi64_mask(bestE, cutPE);
+		bestE = _mm512_mask_set1_epi64(bestE, ht, INT32_MAX);
+		bestPE = _mm512_mask_set1_epi64(bestPE, ht, -1);
+
+		ht = _mm512_cmple_epi64_mask(ugE, cutPE) & _mm512_cmpeq_epi64_mask(ugE, bestE);
+		bestPE = _mm512_mask_mov_epi64(bestPE, ht, threPE);
+
+		_mm512_storeu_si512(r_err, bestE);
+		_mm512_storeu_si512(r_pe, bestPE);
+		_mm512_storeu_si512(uge_mm, ugE);
+	} else {///bd == 1
+		while (k[i] < ai[i]) {
+			err_mm[i] += (VP_mm[i]&(1ULL)); VP_mm[i]>>=1; 
+			err_mm[i] -= (VN_mm[i]&(1ULL)); VN_mm[i]>>=1;
+			++k[i];
+			if ((err_mm[i] <= thre) && (err_mm[i] <= r_err[i])) {
+				r_err[i] = err_mm[i]; r_pe[i] = st[i] + k[i];
+			}
+			if(k[i] == thre) uge_mm[i] = err_mm[i];
+		}
+		if((uge_mm[i] <= thre) && (uge_mm[i] == r_err[i])) r_pe[i] = st[i] + thre;
+	}
+}

Levenshtein_distance.h

@@ -12,6 +12,9 @@
 #include <stdio.h>
 #include "kvec.h"
 
+#define AVX_GS 8
+#define AVX_GS2 4
+
 extern const unsigned char seq_nt4_table[256];
 typedef uint64_t Word;
 typedef uint32_t Word_32;
@@ -548,6 +551,201 @@ inline int32_t pop_trace_back(asg16_v *res, int32_t i, uint16_t *c, uint32_t *le
 	return i;
 }
 
+///compact functions
+#define pop_trac_bpc(in, rc, rb, rl) do {	\
+		(rc) = ((in)>>14);\
+		if((rc) == 1 || (rc) == 2) {(rb) = (((in)>>12)&3); (rl) = ((in)&(0xfff));}\
+		else {(rl) = ((in)&(0x3fff));}\
+	} while (0)	
+
+inline void push_trace_bp(asg16_v *res, uint16_t c, uint16_t b, uint32_t len, uint32_t is_append)
+{
+    uint16_t p, c0, b0, len0, mm; 
+	if((is_append) && (res->n)) {
+		b0 = b;
+		pop_trac_bpc(res->a[res->n-1], c0, b0, len0);
+		if((c == c0) && (b == b0)) {
+			res->n--; len += len0;
+		}
+	}
+	
+	mm = (0x3fff); c0 = c; c <<= 14;
+	if(c0 == 1 || c0 == 2) {
+		mm = (0xfff); c += ((b&3) << 12);
+	}
+
+	
+	while (len >= mm) {
+		p = (c + mm); kv_push(uint16_t, *res, p); len -= mm;
+	}
+	// fprintf(stderr, "[M::%s] c::%u, len::%u\n", __func__, c, len);
+	if(len) {
+		p = (c + len); kv_push(uint16_t, *res, p);
+	}
+}
+
+inline uint32_t pop_trace_bp(asg16_v *res, uint32_t i, uint16_t *c, uint16_t *b, uint32_t *len)
+{
+	(*c) = (res->a[i]>>14); 
+	if((*c) == 1 || (*c) == 2) {
+		(*b) = ((res->a[i]>>12)&3); 
+		(*len) = (res->a[i]&(0xfff));
+	} else {
+		(*b) = (uint16_t)-1; 
+		(*len) = (res->a[i]&(0x3fff));
+	}
+
+	uint32_t sl; uint16_t sb;
+	for (i++; (i < res->n) && ((*c) == (res->a[i]>>14)); i++) {
+		if((*c) == 1 || (*c) == 2) {
+			sb = ((res->a[i]>>12)&3); sl = (res->a[i]&(0xfff));
+		} else {
+			sb = (uint16_t)-1; sl = (res->a[i]&(0x3fff));
+		}
+		if((*b) != sb) break;
+		(*len) += sl;
+	}
+	return i;
+}
+
+inline int64_t pop_trace_bp_rev(asg16_v *res, int64_t i, uint16_t *c, uint16_t *b, uint32_t *len)
+{
+	(*c) = (res->a[i]>>14); 
+	if((*c) == 1 || (*c) == 2) {
+		(*b) = ((res->a[i]>>12)&3); 
+		(*len) = (res->a[i]&(0xfff));
+	} else {
+		(*b) = (uint16_t)-1; 
+		(*len) = (res->a[i]&(0x3fff));
+	}
+
+	uint32_t sl; uint16_t sb;
+	for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) {
+		if((*c) == 1 || (*c) == 2) {
+			sb = ((res->a[i]>>12)&3); sl = (res->a[i]&(0xfff));
+		} else {
+			sb = (uint16_t)-1; sl = (res->a[i]&(0x3fff));
+		}
+		if((*b) != sb) break;
+		(*len) += sl;
+	}
+	return i;
+}
+
+///full functions
+#define pop_trac_bpc_f(in, rc, rbq, rbt, rl) do {	\
+		(rc) = ((in)>>14);\
+		if((rc) == 2 || (rc) == 3) {(rbt) = (((in)>>12)&3); (rl) = ((in)&(0xfff));}\
+		else if((rc) == 1) {(rbt) = (((in)>>12)&3); (rbq) = (((in)>>10)&3); (rl) = ((in)&(0x3ff));}\
+		else {(rl) = ((in)&(0x3fff));}\
+	} while (0)	
+
+inline void push_trace_bp_f(asg16_v *res, uint16_t c, uint16_t bq, uint16_t bt, uint32_t len, uint32_t is_append)
+{
+    uint16_t p, c0 = c, bq0, bt0, len0, mm; 
+	if(c == 3) {
+		bt = bq; bq = (uint16_t)-1;
+	}
+	if((is_append) && (res->n)) {
+		bq0 = bq; bt0 = bt;
+		pop_trac_bpc_f(res->a[res->n-1], c0, bq0, bt0, len0);
+		if((c == c0) && (bq == bq0) && (bt == bt0)) {
+			res->n--; len += len0;
+		}
+	}
+
+
+	c0 = c; c <<= 14;
+	if(c0 == 2 || c0 == 3) {
+		mm = (0xfff); c += ((bt&3) << 12);
+	} else if(c0 == 1) {
+		mm = (0x3ff); c += ((bt&3) << 12); c += ((bq&3) << 10);
+	} else {
+		mm = (0x3fff);
+	}
+
+	
+	while (len >= mm) {
+		p = (c + mm); kv_push(uint16_t, *res, p); len -= mm;
+	}
+	// fprintf(stderr, "[M::%s] c::%u, len::%u\n", __func__, c, len);
+	if(len) {
+		p = (c + len); kv_push(uint16_t, *res, p);
+	}
+}
+
+inline uint32_t pop_trace_bp_f(asg16_v *res, uint32_t i, uint16_t *c, uint16_t *bq, uint16_t *bt, uint32_t *len)
+{
+	(*c) = (res->a[i]>>14); (*bq) = (*bt) = (uint16_t)-1;
+	if((*c) == 2 || (*c) == 3) {
+		(*bt) = ((res->a[i]>>12)&3); 
+		(*len) = (res->a[i]&(0xfff));
+	} else if((*c) == 1) {
+		(*bt) = ((res->a[i]>>12)&3);
+		(*bq) = ((res->a[i]>>10)&3);
+		(*len) = (res->a[i]&(0x3ff));
+	} else {
+		(*len) = (res->a[i]&(0x3fff));
+	}
+
+	uint32_t sl; uint16_t sbq, sbt;
+	for (i++; (i < res->n) && ((*c) == (res->a[i]>>14)); i++) {
+		sbq = sbt = (uint16_t)-1;
+		if((*c) == 2 || (*c) == 3) {
+			sbt = ((res->a[i]>>12)&3); 
+			sl = (res->a[i]&(0xfff));
+		} else if((*c) == 1) {
+			sbt = ((res->a[i]>>12)&3);
+			sbq = ((res->a[i]>>10)&3);
+			sl = (res->a[i]&(0x3ff));
+		} else {
+			sl = (res->a[i]&(0x3fff));
+		}
+		if((*bq) != sbq || (*bt) != sbt) break;
+		(*len) += sl;
+	}
+	if((*c) == 3) {
+		(*bq) = (*bt); (*bt) = (uint16_t)-1;
+	}
+	return i;
+}
+
+inline int64_t pop_trace_bp_rev_f(asg16_v *res, int64_t i, uint16_t *c, uint16_t *bq, uint16_t *bt, uint32_t *len)
+{
+	(*c) = (res->a[i]>>14); (*bq) = (*bt) = (uint16_t)-1;
+	if((*c) == 2 || (*c) == 3) {
+		(*bt) = ((res->a[i]>>12)&3); 
+		(*len) = (res->a[i]&(0xfff));
+	} else if((*c) == 1) {
+		(*bt) = ((res->a[i]>>12)&3);
+		(*bq) = ((res->a[i]>>10)&3);
+		(*len) = (res->a[i]&(0x3ff));
+	} else {
+		(*len) = (res->a[i]&(0x3fff));
+	}
+
+	uint32_t sl; uint16_t sbq, sbt;
+	for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) {
+		sbq = sbt = (uint16_t)-1;
+		if((*c) == 2 || (*c) == 3) {
+			sbt = ((res->a[i]>>12)&3); 
+			sl = (res->a[i]&(0xfff));
+		} else if((*c) == 1) {
+			sbt = ((res->a[i]>>12)&3);
+			sbq = ((res->a[i]>>10)&3);
+			sl = (res->a[i]&(0x3ff));
+		} else {
+			sl = (res->a[i]&(0x3fff));
+		}
+		if((*bq) != sbq || (*bt) != sbt) break;
+		(*len) += sl;
+	}
+	if((*c) == 3) {
+		(*bq) = (*bt); (*bt) = (uint16_t)-1;
+	}
+	return i;
+}
+
 ///511 -> 16 64-bits
 // #define MAX_E 511
 // #define MAX_L 2500
@@ -601,7 +799,7 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
 		if(c == 0) {
 			for (k=0;(k<cl)&&(pstr[pi]==tstr[ti]);k++,pi++,ti++);
 			if(k!=cl) {
-				fprintf(stderr, "ERROR-d-0\n"); 
+				fprintf(stderr, "ERROR-d-0, pi::%d, ti::%d, ci::%u\n", pi, ti, ci); 
 				return 0;
 			}
 		} else {
@@ -609,7 +807,7 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
 			if(c == 1) {
 				for (k=0;(k<cl)&&(pstr[pi]!=tstr[ti]);k++,pi++,ti++);
 				if(k!=cl) {
-					fprintf(stderr, "ERROR-d-1\n"); 
+					fprintf(stderr, "ERROR-d-1, pi::%d, ti::%d, ci::%u\n", pi, ti, ci);
 					return 0;
 				}
 			} else if(c == 2) {///more p
@@ -619,10 +817,20 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
 			}
 		}
 	}
+	
 	if(err != ez->err) {
 		fprintf(stderr, "ERROR-err\n"); 
 		return 0;
 	}
+	if(pi != ez->pe + 1) {
+		fprintf(stderr, "ERROR-pi\n"); 
+		return 0;
+	}
+	if(ti != ez->te + 1) {
+		fprintf(stderr, "ERROR-ti\n"); 
+		return 0;
+	}
+	
 	return 1;
 }
 
@@ -3529,6 +3737,10 @@ inline void ed_band_cal_extension_64_1_w_trace(char *pstr, int32_t pn, char *tst
     return;
 }
 
+void ed_band_cal_semi_64_w_absent_diag_avx4(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe);
+
+void ed_band_cal_semi_64_w_absent_diag_avx8(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe);
+
 inline void ed_band_cal_semi_64_w_absent_diag(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t abs_diag, bit_extz_t *ez)
 {
 	init_base_ed(*ez, thre, pn, tn); ez->ps = ez->pe = -1; ez->ts = 0; ez->te = tn-1;

Makefile

@@ -5,8 +5,8 @@ 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 \
+			POA.o Correct.o Levenshtein_distance.o Levenshtein_avx2.o Levenshtein_avx512.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
-			htab.o hist.o sketch.o anchor.o extract.o sys.o hic.o rcut.o horder.o \
+			htab.o hist.o sketch.o anchor.o extract.o sys.o hic.o rcut.o horder.o ecovlp.o\
 			tovlp.o inter.o kalloc.o gfa_ut.o gchain_map.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
@@ -19,13 +19,22 @@ endif
 .SUFFIXES:.cpp .c .o
 .PHONY:all clean depend
 
+all:$(EXE)
+
 .cpp.o:
 		$(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
 
 .c.o:
 		$(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
 
-all:$(EXE)
+# compiled only with AVX2
+Levenshtein_avx2.o: Levenshtein_avx2.cpp Levenshtein_distance.h
+		$(CXX) -c $(CXXFLAGS) -mavx2 $(CPPFLAGS) $(INCLUDES) $< -o $@
+
+# compiled only with AVX512
+Levenshtein_avx512.o: Levenshtein_avx512.cpp Levenshtein_distance.h
+		$(CXX) -c $(CXXFLAGS) -mavx512f $(CPPFLAGS) $(INCLUDES) $< -o $@
+
 
 $(EXE):$(OBJS) main.o
 		$(CXX) $(CXXFLAGS) $^ -o $@ $(LIBS)
@@ -40,7 +49,7 @@ depend:
 
 Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h
 Assembly.o: Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h
-Assembly.o: kthread.h
+Assembly.o: kthread.h ecovlp.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
@@ -58,6 +67,7 @@ Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 Purge_Dups.o: ksort.h Purge_Dups.h kvec.h kdq.h Overlaps.h Hash_Table.h
 Purge_Dups.o: htab.h Process_Read.h CommandLines.h Correct.h
 Purge_Dups.o: Levenshtein_distance.h POA.h kthread.h
+ecovlp.o: Hash_Table.h Process_Read.h Overlaps.h kthread.h
 Trio.o: khashl.h kthread.h kseq.h Process_Read.h Overlaps.h kvec.h kdq.h
 Trio.o: CommandLines.h htab.h
 anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h

Overlaps.h

@@ -35,6 +35,7 @@
 // #define ALTER_LABLE 2
 // #define HAP_LABLE 4
 #define ug_ext_len 75000
+#define UL_COV_THRES 2
 
 #define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
 #define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
@@ -85,6 +86,12 @@ typedef struct {
 	uint32_t n;
 } idx_emask_t;
 
+typedef struct {
+    uint64_t n, mask;
+    uint8_t *hh;
+    uint64_t tlen, tm;
+} telo_end_t;
+
 typedef struct {
 	///off: start idx in mg128_t * a[];
 	///cnt: how many eles in this chain
@@ -133,8 +140,8 @@ void add_ma_hit_t_alloc(ma_hit_t_alloc* x, ma_hit_t* element);
 void ma_hit_sort_tn(ma_hit_t *a, long long n);
 void ma_hit_sort_qns(ma_hit_t *a, long long n);
 
-int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, 
+int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, char* output_file_name);
-char* output_file_name);
+
 
 typedef struct {
 	uint32_t s:31, del:1, e;
@@ -256,6 +263,7 @@ typedef struct { size_t n, m; uint64_t *a; } asg64_v;
 typedef struct { size_t n, m; uint32_t *a; } asg32_v;
 typedef struct { size_t n, m; ma_utg_t *a;} ma_utg_v;
 typedef struct { asg64_v idx; kv_ul_ov_t srt;} mask_ul_ov_t;
+typedef struct { uint64_t n, m; char *a; } asgchr_v;
 
 typedef struct {
 	ma_utg_v u;
@@ -625,7 +633,7 @@ static inline int count_out_without_del(const asg_t *g, uint32_t v)
 
 void build_string_graph_without_clean(
 int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
-long long n_read, uint64_t* readLen, long long mini_overlap_length, 
+uint64_t n_read, uint64_t* readLen, long long mini_overlap_length, 
 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);
@@ -920,7 +928,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int m
 void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, 
 float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, 
-bub_label_t* b_mask_t, long long no_trio_recover);
+bub_label_t* b_mask_t, long long no_trio_recover, uint8_t *cmk);
 
 typedef struct{
     double weight;
@@ -1115,6 +1123,7 @@ typedef struct{
 	int64_t min_dp;
     bub_label_t* b_mask_t;
 	uint64_t* readLen;
+	telo_end_t *te;
 }ug_opt_t;
 
 typedef struct{
@@ -1130,7 +1139,6 @@ typedef struct{
 	int64_t mini_ovlp;
 }ul_renew_t;
 
-
 void adjust_utg_by_trio(ma_ug_t **ug, asg_t* read_g, uint8_t flag, float drop_rate,
 ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
@@ -1235,5 +1243,12 @@ uint64_t infer_mmhap_copy(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *
 uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_rate, uint64_t bd, ma_ug_t *ug);
 void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src, ug_rid_cov_t *in);
 void gen_ug_rid_cov_t_by_ovlp(kv_u_trans_t *ta, ug_rid_cov_t *cc);
+void rescue_chimeric_reads_aggressive(ma_ug_t *i_ug, asg_t *rg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t chainLenThres, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, 
+kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t, uint8_t *cmk);
+
+#define UC_Read_resize(v, s) do {\
+		if ((v).size<(s)) {REALLOC((v).seq,(s));(v).size=(s);}\
+	} while (0)
 
 #endif

Process_Read.cpp

@@ -46,13 +46,20 @@ void init_All_reads(All_reads* r)
 void destory_All_reads(All_reads* r)
 {
 	uint64_t i = 0;
-	for (i = 0; i < r->total_reads; i++) {
+	for (i = 0; i < r->tqn; i++) {
 		if (r->N_site[i]) free(r->N_site[i]);
 		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);
+		if(r->rsc && r->rsc[i]) free(r->rsc[i]);
 	}
+	for (; i < r->total_reads; i++) {
+		if (r->N_site[i]) free(r->N_site[i]);
+		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);
+	}
+
 	free(r->paf);
 	free(r->reverse_paf);
 	free(r->N_site);
@@ -61,6 +68,7 @@ void destory_All_reads(All_reads* r)
 	free(r->name_index);
 	free(r->read_length);
 	free(r->trio_flag);
+	free(r->rsc);
 	///if (r->pb_regions) free(r->pb_regions);
 }
 
@@ -108,6 +116,15 @@ void write_All_reads(All_reads* r, char* read_file_name)
 	fwrite(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
 	fwrite(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
 
+	uint64_t mm = 2;///1;
+	if(asm_opt.is_sc) {
+		fwrite(&mm, sizeof(mm), 1, fp);
+		fwrite(&(r->tqn), sizeof(r->tqn), 1, fp);
+		for (i = 0; i < r->tqn; i++) {
+			fwrite(r->rsc[i], sizeof(uint8_t), ((r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0)), fp);
+		}
+	}
+
     free(index_name);    
 	fflush(fp);
     fclose(fp);
@@ -123,6 +140,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
 		free(index_name);
         return 0;
     }
+	// fprintf(stderr, "[M::%s]\tindex_name::%s\n", __func__, index_name);
 	int local_adapterLen;
 	int f_flag;
     f_flag = fread(&local_adapterLen, sizeof(local_adapterLen), 1, fp);
@@ -201,6 +219,27 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	}
 	///r->pb_regions = NULL;
 
+	uint64_t mm = 0;
+	if (!feof(fp)) {
+		if((fread(&mm, sizeof(mm), 1, fp)) && (mm == 1 || mm == 2)) {
+			if(mm == 1) {
+				mm = r->total_reads;
+			} else {
+				assert(mm == 2);
+				fread(&mm, sizeof(mm), 1, fp);
+			}
+			r->tqn = mm;
+
+			MALLOC(r->rsc, r->tqn);
+			for (i = 0; i < r->tqn; i++) {
+				MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
+				f_flag += fread(r->rsc[i], sizeof(uint8_t), (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0), fp);
+			}
+		}
+	}
+
+
+
     free(index_name);    
     fclose(fp);
     fprintf(stderr, "Reads has been loaded.\n");
@@ -209,6 +248,206 @@ int load_All_reads(All_reads* r, char* read_file_name)
 }
 
 
+void write_cc_v(cc_v* r, char* read_file_name)
+{
+    fprintf(stderr, "Writing raw reads to disk... \n");
+    char* index_name = (char*)malloc(strlen(read_file_name)+32);
+    sprintf(index_name, "%s.bin", read_file_name);
+    FILE* fp = fopen(index_name, "w"); free(index_name); 
+
+	///typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v;
+	asg16_v *z;
+	uint64_t k, rn = r->n; uint32_t zn; fwrite(&rn, sizeof(rn), 1, fp);
+	for (k = 0; k < r->n; k++) {
+		z = &(r->a[k]);
+		zn = z->n;
+		fwrite(&zn, sizeof(zn), 1, fp);
+		fwrite(z->a, sizeof((*(z->a))), zn, fp);
+	}
+	   
+	fflush(fp); fclose(fp);
+    fprintf(stderr, "Raw reads has been written.\n");
+}
+
+uint8_t load_cc_v(cc_v* r, char* read_file_name)
+{
+    fprintf(stderr, "Loading raw reads... \n");
+    char* index_name = (char*)malloc(strlen(read_file_name)+32);
+    sprintf(index_name, "%s.bin", read_file_name);
+    FILE* fp = fopen(index_name, "r"); free(index_name); 
+	if (!fp) {
+		fprintf(stderr, "No raw read bin.\n");
+        return 0;
+    }
+
+	///typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v;
+	asg16_v *z; int f_flag = 0;
+	uint64_t k, rn; uint32_t zn; f_flag += fread(&rn, sizeof(rn), 1, fp);
+	r->n = r->m = rn; MALLOC(r->a, r->n);
+	for (k = 0; k < r->n; k++) {
+		z = &(r->a[k]);
+		f_flag += fread(&zn, sizeof(zn), 1, fp);
+		z->n = z->m = zn; MALLOC(z->a, z->n);
+		f_flag += fread(z->a, sizeof((*(z->a))), zn, fp);
+	}
+	   
+	fflush(fp); fclose(fp);
+    fprintf(stderr, "Raw reads has been loaded.\n");
+	return 1;
+}
+
+void read_ma(ma_hit_t* x, FILE* fp)
+{
+    int f_flag;
+    f_flag = fread(&(x->qns), sizeof(x->qns), 1, fp);
+    f_flag += fread(&(x->qe), sizeof(x->qe), 1, fp);
+    f_flag += fread(&(x->tn), sizeof(x->tn), 1, fp);
+    f_flag += fread(&(x->ts), sizeof(x->ts), 1, fp);
+    f_flag += fread(&(x->te), sizeof(x->te), 1, fp);
+    f_flag += fread(&(x->el), sizeof(x->el), 1, fp);
+    f_flag += fread(&(x->no_l_indel), sizeof(x->no_l_indel), 1, fp);
+
+    uint32_t t;
+    f_flag += fread(&(t), sizeof(t), 1, fp);
+    x->ml = t;
+
+    f_flag += fread(&(t), sizeof(t), 1, fp);
+    x->rev = t;
+    
+    f_flag += fread(&(t), sizeof(t), 1, fp);
+    x->bl = t;
+
+    f_flag += fread(&(t), sizeof(t), 1, fp);
+    x->del = t;
+}
+
+
+int append_All_reads(All_reads* r, char *idx, uint32_t id)
+{
+	char *gfa_name = NULL; MALLOC(gfa_name, (strlen(idx)+100));
+	FILE *fp = NULL, *fpo = NULL;
+
+	sprintf(gfa_name, "%s.ec%u.bin", idx, id);
+    fp = fopen(gfa_name, "r");
+	if (!fp) {free(gfa_name); return 0;}
+
+	sprintf(gfa_name, "%s.ovlp%u.source.bin", idx, id);
+	fpo = fopen(gfa_name, "r");
+	if (!fpo) {free(gfa_name); fclose(fp); return 0;}
+	free(gfa_name);
+
+	int local_adapterLen, f_flag;
+    f_flag = fread(&local_adapterLen, sizeof(local_adapterLen), 1, fp);
+    if(local_adapterLen != asm_opt.adapterLen) {
+        fprintf(stderr, "[M::%s] the adapterLen of index is: %d, but the adapterLen set by user is: %d\n", 
+		__func__, local_adapterLen, asm_opt.adapterLen);
+		exit(1);
+    }
+	uint64_t index_size0, name_index_size0, total_reads0, total_reads_bases0, total_name_length0;
+	index_size0 = r->index_size; 
+	total_reads0 = r->total_reads; 
+	total_reads_bases0 = r->total_reads_bases;
+	total_name_length0 = r->total_name_length;
+	name_index_size0 = ((total_reads0)?(total_reads0+1):(0));///r->name_index_size;
+
+    f_flag += fread(&r->index_size, sizeof(r->index_size), 1, fp);
+	f_flag += fread(&r->name_index_size, sizeof(r->name_index_size), 1, fp);
+	f_flag += fread(&r->total_reads, sizeof(r->total_reads), 1, fp);
+	f_flag += fread(&r->total_reads_bases, sizeof(r->total_reads_bases), 1, fp);
+	f_flag += fread(&r->total_name_length, sizeof(r->total_name_length), 1, fp);
+
+	r->index_size += index_size0;
+	r->name_index_size += name_index_size0; if(name_index_size0) r->name_index_size--;
+	r->total_reads += total_reads0;
+	r->total_reads_bases += total_reads_bases0;
+	r->total_name_length += total_name_length0;
+
+
+	uint64_t i = 0, zero = 0, k;
+	REALLOC(r->N_site, r->total_reads);
+
+	for (i = total_reads0; i < r->total_reads; i++) {
+		f_flag += fread(&zero, sizeof(zero), 1, fp);
+
+		if (zero) {
+			r->N_site[i] = (uint64_t*)malloc(sizeof(uint64_t)*(zero + 1));
+			r->N_site[i][0] = zero;
+			if (r->N_site[i][0]) {
+				f_flag += fread(r->N_site[i]+1, sizeof(r->N_site[i][0]), r->N_site[i][0], fp);
+			}
+		} else {
+			r->N_site[i] = NULL;
+		}
+	}
+
+	REALLOC(r->read_length, r->total_reads);
+	f_flag += fread(r->read_length + total_reads0, sizeof(uint64_t), r->total_reads-total_reads0, fp);
+
+	REALLOC(r->read_size, r->total_reads);
+	memcpy (r->read_size + total_reads0, r->read_length + total_reads0, sizeof(uint64_t)*(r->total_reads-total_reads0));
+
+	REALLOC(r->read_sperate, r->total_reads);
+	for (i = total_reads0; i < r->total_reads; i++) {
+		r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
+		f_flag += fread(r->read_sperate[i], sizeof(uint8_t), r->read_length[i]/4+1, fp);
+	}
+
+
+	REALLOC(r->name, r->total_name_length);
+	f_flag += fread(r->name + total_name_length0, sizeof(char), r->total_name_length - total_name_length0, fp);
+
+	REALLOC(r->name_index, r->name_index_size); uint64_t sft = 0;
+	if(name_index_size0) sft = r->name_index[--name_index_size0];
+	// fprintf(stderr, "name_index_size0::%lu, total_reads0::%lu, sft::%lu\n", name_index_size0, total_reads0, sft);
+	f_flag += fread(r->name_index + name_index_size0, sizeof(uint64_t), r->name_index_size-name_index_size0, fp);
+	if(sft) {
+		for (i = name_index_size0; i < r->name_index_size; i++) r->name_index[i] += sft;
+	}
+
+	/****************************may have bugs********************************/
+	REALLOC(r->trio_flag, r->total_reads);
+	f_flag += fread(r->trio_flag + total_reads0, sizeof(uint8_t), r->total_reads - total_reads0, fp);
+
+	int hom_cov0 = asm_opt.hom_cov, het_cov0 = asm_opt.het_cov;
+	f_flag += fread(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp); asm_opt.hom_cov += hom_cov0;
+    f_flag += fread(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp); asm_opt.het_cov += het_cov0;
+	/****************************may have bugs********************************/
+
+	REALLOC(r->cigars, r->total_reads);
+	REALLOC(r->second_round_cigar, r->total_reads);
+	for (i = total_reads0; i < r->total_reads; i++) {
+		r->second_round_cigar[i].size = r->cigars[i].size = 0;
+		r->second_round_cigar[i].length = r->cigars[i].length = 0;
+		r->second_round_cigar[i].record = r->cigars[i].record = NULL;
+
+		r->second_round_cigar[i].lost_base_size = r->cigars[i].lost_base_size = 0;
+		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;
+
+	REALLOC(r->paf, r->total_reads);
+	memset(r->paf+total_reads0, 0, sizeof((*(r->paf)))*(r->total_reads - total_reads0));
+	long long n_read; f_flag = fread(&n_read, sizeof(n_read), 1, fpo); ma_hit_t t;
+	for (i = total_reads0; i < r->total_reads; i++) {
+		f_flag += fread(&(r->paf[i].is_fully_corrected), sizeof(r->paf[i].is_fully_corrected), 1, fpo);
+        f_flag += fread(&(r->paf[i].is_abnormal), sizeof(r->paf[i].is_abnormal), 1, fpo);
+        f_flag += fread(&(r->paf[i].length), sizeof(r->paf[i].length), 1, fpo);
+
+        if(r->paf[i].length == 0) continue;
+        for (k = 0; k < r->paf[i].length; k++) read_ma(&t, fpo);
+		r->paf[i].length = 0;
+	}
+
+	REALLOC(r->reverse_paf, r->total_reads); 
+	memset(r->reverse_paf+total_reads0, 0, sizeof((*(r->reverse_paf)))*(r->total_reads - total_reads0));
+
+    fclose(fp); fclose(fpo);
+    fprintf(stderr, "Reads has been loaded.\n");
+	return 1;
+}
+
+
 int destory_read_bin(All_reads* r)
 {
 
@@ -255,16 +494,24 @@ void ha_insert_read_len(All_reads *r, int read_len, int name_len)
 	r->name_index[r->total_reads] = r->name_index[r->total_reads - 1] + name_len;
 }
 
-void malloc_All_reads(All_reads* r)
+void malloc_All_reads(All_reads *r)
 {
 	r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
 	memcpy(r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
 
 	r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads);
-	long long i = 0;
+	if(asm_opt.is_sc && r->tqn) MALLOC(r->rsc, r->tqn);
-	for (i = 0; i < (long long)r->total_reads; i++)
+	assert(r->tqn <= r->total_reads);
-	{
+
-		r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
+	uint64_t i = 0;
+	if(r->rsc) {
+		for (i = 0; i < r->tqn; i++) {
+			MALLOC(r->read_sperate[i], (r->read_length[i]/4+1));
+			MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
+		}
+	}
+	for (; i < r->total_reads; i++) {
+		MALLOC(r->read_sperate[i], (r->read_length[i]/4+1));
 	}
 
 	r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
@@ -273,8 +520,7 @@ void malloc_All_reads(All_reads* r)
 	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++)
+	for (i = 0; i < r->total_reads; i++) {
-	{
 		r->second_round_cigar[i].size = r->cigars[i].size = 0;
 		r->second_round_cigar[i].length = r->cigars[i].length = 0;
 		r->second_round_cigar[i].record = r->cigars[i].record = NULL;
@@ -672,6 +918,130 @@ void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_sit
 	}
 }
 
+void convert_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitu, uint64_t rev, uint64_t sc_off)
+{
+	uint64_t i = 0; uint8_t c = 0, sc;
+	// fprintf(stderr, "\n[M::%s]\n", __func__);
+	for (i = 0; i < src_l; i++) {
+		for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+		if(c >= bitu) c = bitu - 1;
+		dest[(rev?(src_l-i-1):(i))] = c;
+		// fprintf(stderr, "%u->%u\n", sc, c);
+	}
+}
+
+void ha_compress_qual_bit(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn)
+{
+
+	uint64_t i = 0, k, bit_r = 8/bitn, dest_i = 0;
+	uint8_t tmp = 0, c = 0;
+
+	for (i = 0; i + bit_r <= src_l;) {
+		for (k = tmp = 0; k < bit_r; k++) {
+			c = ((uint8_t)src[i]);
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		dest[dest_i++] = tmp;
+	}
+
+	if(i < src_l) {
+		for (k = tmp = 0; i < src_l; k++) {
+			c = ((uint8_t)src[i]);
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		 
+		dest[dest_i++] = (tmp<<(8-(bitn*k)));
+	}
+}
+
+void ha_compress_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn, uint64_t sc_off)
+{
+
+	uint64_t i = 0, k, bit_r = 8/bitn, dest_i = 0, bitu = (1<<bitn);
+	uint8_t tmp = 0, c = 0, sc;
+
+	for (i = 0; i + bit_r <= src_l;) {
+		for (k = tmp = 0; k < bit_r; k++) {
+			for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+			if(c >= bitu) c = bitu - 1;
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		dest[dest_i++] = tmp;
+	}
+
+	if(i < src_l) {
+		for (k = tmp = 0; i < src_l; k++) {
+			for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+			if(c >= bitu) c = bitu - 1;
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		 
+		dest[dest_i++] = (tmp<<(8-(bitn*k)));
+	}
+}
+
+///[s, e)
+int64_t retrive_bqual(asg8_v *dv, uint8_t *ds, uint64_t id, int64_t s, int64_t e, uint8_t rev, int64_t bitn)
+{
+	int64_t rl = Get_READ_LENGTH(R_INF, id), l; 
+	if(s < 0) {s = 0;} if(e < 0) {e = rl;}
+	if(s >= e || e > rl) return -1;
+
+	uint8_t *da = NULL, *src = Get_QUAL(R_INF, id), mm = (((uint8_t)1)<<bitn)-1, mlf = 8 - bitn, mrf;
+	int64_t bitr = 8/bitn, dk, sk, swk;
+	l = e - s;
+	if(dv) {
+		kv_resize(uint8_t, *dv, ((uint64_t)l)); da = dv->a; dv->n = l;
+	} else {
+		da = ds;
+	}
+
+	if(!rev) {
+		dk = 0; sk = s;
+
+		mrf = ((s%bitr)*bitn);
+		// if(s == 21519 && e == 22332) {
+		// 	fprintf(stderr, "+[M::%s] id::%lu, in::[%ld, %ld), rev::%u, bitn::%ld, bitr::%ld, mrf::%u\n", __func__, id, s, e, rev, bitn, bitr, mrf);
+		// }
+		if(mrf) {
+			for (swk = sk/bitr; mrf < 8 && sk < e; mrf += bitn, sk++) da[dk++] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		
+		for (swk = sk/bitr; (sk + bitr) <= e; sk += bitr, swk++) {
+			for (mrf = 0; mrf < 8; mrf += bitn) da[dk++] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+
+		if(sk < e) {
+			for (mrf = 0; sk < e; mrf += bitn, sk++) da[dk++] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		// if(dk != l) {
+		// 	fprintf(stderr, "+[M::%s] id::%lu, in::[%ld, %ld), rev::%u, bitn::%ld, bitr::%ld\n", __func__, id, s, e, rev, bitn, bitr);
+		// }
+		assert(dk == l);
+	} else {
+		sk = s; s = e; e = sk;
+		s = rl - s; e = rl - e;
+		dk = l; sk = s; 
+
+		mrf = ((s%bitr)*bitn);
+		if(mrf) {
+			for (swk = sk/bitr; mrf < 8 && sk < e; mrf += bitn, sk++) da[--dk] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		
+		for (swk = sk/bitr; (sk + bitr) <= e; sk += bitr, swk++) {
+			for (mrf = 0; mrf < 8; mrf += bitn) da[--dk] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+
+		if(sk < e) {
+			for (mrf = 0; sk < e; mrf += bitn, sk++) da[--dk] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		
+		assert(dk == 0);
+	}
+
+	return l;
+}
+
 void reverse_complement(char* pattern, uint64_t length)
 {
 	uint64_t i = 0;
@@ -694,6 +1064,24 @@ void reverse_complement(char* pattern, uint64_t length)
 	}
 }
 
+void print_fastq(FILE *fp, char *id, char *bs, char *qual, uint64_t bitu, uint64_t sc_off)
+{
+	uint64_t i = 0, ql = strlen(qual); uint8_t c = 0, sc;
+	
+	if(fp) fprintf(fp, "@%s\n%s\n+\n", id, bs);
+	else fprintf(stdout, "@%s\n%s\n+\n", id, bs);
+
+	for (i = 0; i < ql; i++) {
+		for (c = 0, sc = ((uint8_t)qual[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+		if(c >= bitu) c = bitu - 1;
+		if(fp) fprintf(fp, "%u", c);
+		else fprintf(stdout, "%u", c);
+	}
+
+	if(fp) fprintf(fp, "\n");
+	else fprintf(stdout, "\n");
+}
+
 
 void init_Debug_reads(Debug_reads* x, const char* file)
 {
@@ -1179,7 +1567,7 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
 		p->rlen = str_l; 
 		// fprintf(stderr, "str_l->%ld, str->%u\n", str_l, str?1:0);
 		
-		if(o == NULL || on == 0) on = 0; en = 0;
+		if(o == NULL || on == 0) {on = 0;} en = 0;
 		for (i = on-1, st = et = str_l; i >= 0; i--) {
 			z = &(o[i]);
 			if(z->el) {
@@ -1455,6 +1843,62 @@ void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_
 	}
 }
 
+void retrieve_u_seq_fast(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l, void *km)
+{
+	if(u->m == 0 || u->n == 0) return;
+	if(l < 0) l = u->len;
+	char *r = NULL, *a = NULL;
+	int64_t e = s + l, ssp, sep, rs, re, des_i;
+	uint64_t k, rId, ori, r_l;
+	if(i_r) {
+        i_r->length = l; i_r->RID = 0;
+        if(i_r->length > i_r->size) {
+            i_r->size = i_r->length;
+			if(!km) REALLOC(i_r->seq, i_r->size);
+			else KREALLOC(km, i_r->seq, i_r->size);
+            // i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size));
+        }
+        r = i_r->seq;
+    }
+    if(i_s) r = i_s;
+
+	if(u->s) {
+		memcpy(r, u->s + s, l * sizeof((*(u->s))));
+	} else {
+		if(strand == 1) {
+			sep = u->len - s; 
+			ssp = u->len - e;
+			s = ssp; e = sep; 
+		}
+		for (k = l = des_i = 0; k < u->n; k++) {
+			rId = u->a[k]>>33;
+			ori = u->a[k]>>32&1;
+			r_l = (uint32_t)u->a[k];
+			if(r_l == 0) continue;
+			ssp = l; sep = l + r_l;
+			l += r_l;
+			if(sep <= s) continue;
+			if(ssp >= e) break;
+			rs = MAX(ssp, s); re = MIN(sep, e); 
+			a = r + des_i; des_i += re - rs;
+			recover_UC_Read_sub_region(a, rs-ssp, re-rs, ori, &R_INF, rId);
+		}
+	}
+
+	if(strand == 1) {
+		char t;
+		re = (e - s);
+		l = re>>1; 
+		for (k = 0; k < (uint64_t)l; k++) {
+			des_i = re - k - 1;
+			t = r[des_i];
+			r[des_i] = RC_CHAR(r[k]);
+			r[k] = RC_CHAR(t);
+		}
+		if(re&1) r[l] = RC_CHAR(r[l]);
+	}
+}
+
 uint32_t retrieve_u_cov(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t pos, uint8_t dir, int64_t *pi)
 {
 	uint64_t *a = ul->cc->interval.a + ul->cc->idx[id], cc = 0, ff = 0;

Process_Read.h

@@ -8,6 +8,7 @@
 #include <zlib.h>
 #include "Overlaps.h"
 #include "CommandLines.h"
+#include "Levenshtein_distance.h"
 ///#include "Hash_Table.h"
 
 #define READ_INIT_NUMBER 1000
@@ -22,6 +23,7 @@
 #define Get_NAME_LENGTH(R_INF, ID) ((R_INF).name_index[(ID)+1] - (R_INF).name_index[(ID)])
 ///#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_QUAL(R_INF, ID) (R_INF).rsc[(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)
@@ -38,6 +40,8 @@ extern char rc_Table[6];
 
 void init_aux_table();
 
+typedef struct { size_t n, m; uint8_t *a; } asg8_v;
+
 typedef struct
 {
     uint64_t x_id;
@@ -107,6 +111,8 @@ typedef struct
 #define CHAIN_MATCH 1
 #define CHAIN_UNMATCH 0.334
 
+#define NEC 1
+
 typedef struct
 {
 	uint64_t** N_site;
@@ -117,6 +123,7 @@ typedef struct
 	uint64_t* read_length;
 	uint64_t* read_size;
 	uint8_t* trio_flag;
+    uint8_t** rsc;
 
 	///seq start pos in uint8_t* read
 	///do not need it
@@ -127,8 +134,10 @@ typedef struct
 	uint64_t* name_index;
 	uint64_t name_index_size;
 	uint64_t total_reads;
+    uint64_t tqn;
 	uint64_t total_reads_bases;
 	uint64_t total_name_length;
+    uint64_t tr[2];
 
 	Compressed_Cigar_record* cigars; 
 	Compressed_Cigar_record* second_round_cigar;
@@ -136,11 +145,16 @@ typedef struct
     ma_hit_t_alloc* paf;
     ma_hit_t_alloc* reverse_paf;
 
+    uint8_t is_syn;
+
     ///kvec_t_u64_warp* pb_regions;
 } All_reads;
 
 extern All_reads R_INF;
 
+typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v;
+extern cc_v scb;
+
 typedef struct
 {
 	char* seq;
@@ -221,6 +235,7 @@ 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);
 void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ);
+void ha_compress_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn, uint64_t sc_off);
 void init_UC_Read(UC_Read* r);
 void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID);
 void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID);
@@ -229,6 +244,9 @@ void destory_UC_Read(UC_Read* r);
 void reverse_complement(char* pattern, uint64_t length);
 void write_All_reads(All_reads* r, char* read_file_name);
 int load_All_reads(All_reads* r, char* read_file_name);
+uint8_t load_cc_v(cc_v* r, char* read_file_name);
+void write_cc_v(cc_v* r, char* read_file_name);
+int append_All_reads(All_reads* r, char *idx, uint32_t id);
 void destory_All_reads(All_reads* r);
 int destory_read_bin(All_reads* r);
 void init_Debug_reads(Debug_reads* x, const char* file);
@@ -249,5 +267,20 @@ void write_compress_base_disk(FILE *fp, uint64_t ul_rid, char *str, uint32_t len
 int64_t load_compress_base_disk(FILE *fp, uint64_t *ul_rid, char *dest, uint32_t *len, ul_vec_t *buf);
 scaf_res_t *init_scaf_res_t(uint32_t n);
 void destroy_scaf_res_t(scaf_res_t *p);
+void read_ma(ma_hit_t* x, FILE* fp);
+void retrieve_u_seq_fast(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l, void *km);
+
+const uint64_t sc_tb[8] = { 
+	10, 20, 30, 40, 50, 60, 70, 80
+};
+
+#define sc_bn 2
+#define sc_bm ((((uint64_t)1)<<sc_bn)-1)
+#define sc_wn 5
+
+void convert_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitu, uint64_t rev, uint64_t sc_off);
+int64_t retrive_bqual(asg8_v *dv, uint8_t *ds, uint64_t id, int64_t s, int64_t e, uint8_t rev, int64_t bitn);
+void print_fastq(FILE *fp, char *id, char *bs, char *qual, uint64_t bitu, uint64_t sc_off);
+void ha_compress_qual_bit(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn);
 
 #endif

README.md

@@ -15,6 +15,9 @@ hifiasm -o CHM13.asm -t32 -l0 CHM13-HiFi.fa.gz 2> CHM13.asm.log
 # Assemble heterozygous genomes with built-in duplication purging
 hifiasm -o HG002.asm -t32 HG002-file1.fq.gz HG002-file2.fq.gz
 
+# Assemble genomes with ONT R10 reads rather than PacBio HiFi reads using the latest release of hifiasm (>0.21.0-r686)
+hifiasm -o HG002.asm --ont -t32 HG002-ont.fq.gz
+
 # Hi-C phasing with paired-end short reads in two FASTQ files
 hifiasm -o HG002.asm --h1 read1.fq.gz --h2 read2.fq.gz HG002-HiFi.fq.gz
 
@@ -23,8 +26,18 @@ yak count -b37 -t16 -o pat.yak <(cat pat_1.fq.gz pat_2.fq.gz) <(cat pat_1.fq.gz
 yak count -b37 -t16 -o mat.yak <(cat mat_1.fq.gz mat_2.fq.gz) <(cat mat_1.fq.gz mat_2.fq.gz)
 hifiasm -o HG002.asm -t32 -1 pat.yak -2 mat.yak HG002-HiFi.fa.gz
 
-# Single-sample telomere-to-telomere assembly with HiFi, ultralong and Hi-C reads
+# Improve contiguity for diploid genome assembly by self-scaffolding (`--dual-scaf`)
+hifiasm -o HG002.asm --dual-scaf --h1 read1.fq.gz --h2 read2.fq.gz HG002-HiFi.fq.gz
+
+# Preserve more telomeres for human genomes (`--telo-m CCCTAA`)
+hifiasm -o HG002.asm --telo-m CCCTAA --h1 read1.fq.gz --h2 read2.fq.gz HG002-HiFi.fq.gz
+
+# Hybrid assembly with HiFi, ultralong and Hi-C reads
 hifiasm -o HG002.asm --h1 read1.fq.gz --h2 read2.fq.gz --ul ul.fq.gz HG002-HiFi.fq.gz
+
+# Single-sample telomere-to-telomere assembly for diploid human genomes
+hifiasm -o HG002.asm --dual-scaf --telo-m CCCTAA --h1 read1.fq.gz --h2 read2.fq.gz --ul ul.fq.gz HG002-HiFi.fq.gz
+
 ```
 See [tutorial][tutorial] for more details. 
 
@@ -35,6 +48,7 @@ See [tutorial][tutorial] for more details.
 - [Why Hifiasm?](#why)
 - [Usage](#use)
   - [Assembling HiFi reads without additional data types](#hifionly)
+  - [Assembling ONT reads](#ontonly)
   - [Hi-C integration](#hic)
   - [Trio binning](#trio)
   - [Ultra-long ONT integration](#ul)
@@ -106,6 +120,16 @@ bloom filter which takes 16GB memory at the beginning. For genomes much larger
 than human, applying `-f38` or even `-f39` is preferred to save memory on k-mer
 counting.
 
+### <a name="ontonly"></a>Assembling ONT reads
+
+Since version 0.21.0 (r686), hifiasm can support ONT assembly using ONT simplex R10 reads. 
+To enable this feature, add the `--ont` option as shown below:
+```sh
+hifiasm -t64 --ont -o ONT.asm ONT.read.fastq.gz
+```
+Please note that this module requires input reads in FASTQ format.
+
+
 ### <a name="hic"></a>Hi-C integration
 
 Hifiasm can generate a pair of haplotype-resolved assemblies with paired-end
@@ -154,11 +178,29 @@ For the single-sample telomere-to-telomere assembly with Hi-C reads:
 ```sh
 hifiasm -o NA12878.asm -t32 --ul ul.fq.gz --h1 read1.fq.gz --h2 read2.fq.gz HiFi-reads.fq.gz
 ```
-For the trio-binning telomere-to-telomere assembly；
+For the trio-binning telomere-to-telomere assembly:
 ```sh
 hifiasm -o NA12878.asm -t32 --ul ul.fq.gz -1 pat.yak -2 mat.yak HiFi-reads.fq.gz
 ```
 
+### <a name="ul"></a>Self-scaffolding
+
+For diploid haplotype-resolved genome assembly, hifiasm can further enhance assembly contiguity 
+by introducing scaffolding. It leverages the assemblies of the two haplotypes to scaffold each other. 
+Specifically, if there is a gap within the haplotype 1 assembly, hifiasm will use the corresponding 
+homologous region in haplotype 2 to scaffold haplotype 1. Below is an example using the `--dual-scaf` option.
+```sh
+hifiasm -o NA12878.asm -t32 --dual-scaf HiFi-reads.fq.gz
+```
+
+### <a name="ul"></a>Preserve more telomeres for T2T assemblies
+
+Hifiasm can preserve more telomeres by specifying the telomere motif using the `--telo-m` option. 
+Below is an example applied to human genome assembly.
+```sh
+hifiasm -o NA12878.asm -t32 --telo-m CCCTAA HiFi-reads.fq.gz
+```
+
 ### <a name="output"></a>Output files
 
 Hifiasm generates different types of assemblies based on the input data. 
@@ -247,3 +289,8 @@ If you use hifiasm in your work, please cite:
 > Haplotype-resolved assembly of diploid genomes without parental data. 
 > *Nature Biotechnology*, **40**:1332–1335.
 > https://doi.org/10.1038/s41587-022-01261-x
+
+> Cheng, H., Asri, M., Lucas, J., Koren, S., Li, H. (2024)
+> Scalable telomere-to-telomere assembly for diploid and polyploid genomes with double graph. 
+> *Nat Methods*, **21**:967-970.
+> https://doi.org/10.1038/s41592-024-02269-8

Trio.cpp

@@ -388,6 +388,31 @@ inline void phrase_hstatus(char *s, char **rname, uint32_t *hid)
 	*hid = atoi(id);
 }
 
+inline void phrase_hchar(char *s, char **rname, uint32_t *hid)
+{
+	*rname = NULL; *hid = (uint32_t)-1;
+	uint32_t tot = 0, k, l, z, sl = strlen(s); 
+	for (k = 1, l = 0; k <= sl; k++) {
+		if((k == sl) || (s[k] == '\t') || (s[k] == ' ')) {
+			if(k > l) {
+				s[k] = 0; 
+				if(tot == 0) {
+					*rname = s + l;
+				} else if(tot == 1) {
+					for (z = l; (z < k) && (s[z] >= '0') && (s[z] <= '9'); ++z);
+					if(z < k) {
+						fprintf(stderr, "ERROR: wrong hap id\n");
+						return;
+					}
+					*hid = atoi(s + l);
+				}
+				tot++;
+			}
+			l = k + 1;
+		}
+	}
+}
+
 uint32_t *ha_polybin_list(const hifiasm_opt_t *opt)
 {
 	int64_t i;
@@ -447,6 +472,74 @@ uint32_t *ha_polybin_list(const hifiasm_opt_t *opt)
 	return ss;
 }
 
+uint32_t *ha_charbin_list(const hifiasm_opt_t *opt, uint8_t **idx, uint32_t *idx_n)
+{
+	int64_t i;
+	khint_t k;
+	cstr_ht_t *h; (*idx) = NULL; *idx_n = 0;
+	assert(R_INF.total_reads < (uint32_t)-1);
+	h = cstr_ht_init();
+	for (i = 0; i < (int64_t)R_INF.total_reads; ++i) {
+		int absent;
+		char *str = (char*)calloc(Get_NAME_LENGTH(R_INF, i) + 1, 1);
+		strncpy(str, Get_NAME(R_INF, i), Get_NAME_LENGTH(R_INF, i));
+		k = cstr_ht_put(h, str, &absent);
+		if (absent) kh_val(h, k) = i;
+	}
+	fprintf(stderr, "[M::%s::%.3f*%.2f] created the hash table for read names\n", __func__, yak_realtime(), yak_cpu_usage());
+
+	gzFile fp;
+	kstream_t *ks;
+	kstring_t str = {0,0,0};
+	char *rname = NULL;
+	uint32_t hid, mhid = 0, *ss = NULL;
+	int dret;
+	int64_t n_tot = 0, n_bin = 0;
+	fp = gzopen(opt->fn_chr_bin, "r");
+	if (fp == 0) {
+		fprintf(stderr, "ERROR: failed to open file '%s'\n", opt->fn_chr_bin);
+		for (k = 0; k < kh_end(h); ++k)
+			if (kh_exist(h, k))
+				free((char*)kh_key(h, k));
+		cstr_ht_destroy(h);
+		return NULL;
+	}
+	MALLOC(ss, R_INF.total_reads); memset(ss, -1, sizeof((*ss))*R_INF.total_reads);
+	ks = ks_init(fp);
+	while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
+		khint_t k; ++n_tot;
+		phrase_hchar(str.s, &rname, &hid);
+		if((!(*rname)) || hid == (uint32_t)-1) {
+			fprintf(stderr, "ERROR: wrong hap id\n");
+			continue;
+		}
+		k = cstr_ht_get(h, rname);
+		if (k != kh_end(h)) {
+			ss[kh_val(h, k)] = hid;
+			if(hid > mhid) mhid = hid;
+			++n_bin;
+			// fprintf(stderr, "%s\t%u\trid::%ld\n", rname, hid, kh_val(h, k));
+		}
+	}
+	free(str.s);
+	ks_destroy(ks);
+	gzclose(fp);
+	
+	for (k = 0; k < kh_end(h); ++k)
+		if (kh_exist(h, k))
+			free((char*)kh_key(h, k));
+	cstr_ht_destroy(h);
+
+	mhid++; CALLOC((*idx), mhid); 
+	for (i = 0; i < (int64_t)R_INF.total_reads; ++i) {
+		if(ss[i] >= mhid) continue;
+		(*idx)[ss[i]] = 1;
+	}
+	*idx_n = mhid;
+	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads with external lists\n", __func__, yak_realtime(), yak_cpu_usage());
+	return ss;
+}
+
 void ha_triobin(const hifiasm_opt_t *opt)
 {
 	memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads * sizeof(uint8_t));

ecovlp.h

@@ -0,0 +1,24 @@
+#ifndef __ECOVLP_PARSER__
+#define __ECOVLP_PARSER__
+
+#define __STDC_LIMIT_MACROS
+#include <stdint.h>
+#include "Hash_Table.h"
+#include "Process_Read.h"
+#include "kdq.h"
+
+
+void prt_chain(overlap_region_alloc *o);
+void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, uint64_t is_sv, uint64_t *tot_b, uint64_t *tot_e);
+void sl_ec_r(uint64_t n_thre, uint64_t n_a);
+void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx);
+void handle_chemical_r(uint64_t n_thre, uint64_t n_a);
+void handle_chemical_arc(uint64_t n_thre, uint64_t n_a);
+uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a);
+void cal_ec_r_dbg(uint64_t n_thre, uint64_t n_a);
+void write_ec_reads(const char *suffix_ou, cc_v *cvt, uint8_t is_rev);
+void destroy_cc_v(cc_v *z);
+void gen_gfa_bam(ma_ug_t *ug, uint64_t n_a);
+void clean_arc_rf(uint64_t n_thre, uint64_t n_a);
+
+#endif

gfa_ut.h

@@ -16,10 +16,10 @@ typedef struct {
 } sset_aux;
 
 void ul_clean_gfa(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, 
-double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file);
+double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, uint8_t *cmk, char *o_file);
-uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_ou, R_to_U *ru);
+uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_ou, R_to_U *ru, telo_end_t *te);
-void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t normal_len, uint32_t pop_chimer, asg64_v *dbg);
+void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t normal_len, uint32_t pop_chimer, asg64_v *dbg, telo_end_t *te);
-void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t ou_thres);
+void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t ou_thres, telo_end_t *te);
 void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio, uint32_t min_diff, float ou_rat/**, asg64_v *dbg**/);
 void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, 
 uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len);
@@ -33,7 +33,7 @@ void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t m
 void normalize_gou(asg_t *g);
 void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);
 asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len);
-void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover);
+void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover, uint8_t *cmk);
 // void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail);
 bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het, uint8_t avoid_het);
 void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt);

hic.cpp

@@ -6569,7 +6569,7 @@ min_cut_t* m, hc_links* link, G_partition* x)
     if(res->full_bub == 0)
     {
         res->a.n = 0;
-        uint32_t v, u = 0, uv, k_n, pre_n = x->n;
+        uint32_t v, u = 0, uv = UINT32_MAX, k_n, pre_n = x->n;
         hc_linkeage* t = NULL;
         x->n--;
         for (i = 0; i < n; i++)

horder.cpp

@@ -32,7 +32,7 @@ KRADIX_SORT_INIT(osg, osg_arc_t, osg_arc_key, member_size(osg_arc_t, u))
 #define BREAK_CUTOFF 0.1
 #define BREAK_BOUNDARY 0.015
 void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
-int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub);
+int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub, uint32_t max_ext);
 
 typedef struct {
 	uint64_t ruid;
@@ -714,7 +714,7 @@ void update_u_hits(kvec_pe_hit *u_hits, kvec_pe_hit *r_hits, ma_ug_t* ug, asg_t*
                 kv_pushp(hit_aux_t, x, &p);
                 p->ruid = u->a[i]>>32; 
                 p->ruid <<= 32;
-                p->ruid |= v;
+                p->ruid |= v;///rid|rev|uid
                 p->off = offset;
                 offset += (uint32_t)u->a[i];
             }
@@ -725,8 +725,8 @@ void update_u_hits(kvec_pe_hit *u_hits, kvec_pe_hit *r_hits, ma_ug_t* ug, asg_t*
         }
     }
 
-    radix_sort_hit_aux_ruid(x.a, x.a + x.n);
+    radix_sort_hit_aux_ruid(x.a, x.a + x.n);///sort by (rid|rev|uid)
-    x.idx.n = x.idx.m = (x.n?(x.a[x.n-1].ruid>>33)+1:0);///how many unitigs
+    x.idx.n = x.idx.m = (x.n?(x.a[x.n-1].ruid>>33)+1:0);///how many reads?
     CALLOC(x.idx.a, x.idx.n);
     for (k = 1, l = 0; k <= x.n; ++k) 
     {   
@@ -2185,7 +2185,7 @@ h_covs *res, h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique
         span_e = MIN(span_e, len-1) + 1;
         //if(span_e - span_s <= ulen*BREAK_CUTOFF)//need it or not?
         {
-            if(span_s >= sPos && span_e <= ePos)
+            if(span_s >= sPos && span_e <= ePos)///test the density of this local region; bug -> should use any HiC pairs that are overlapped with [sPos, sPoe), instead of fully covered by [sPos, sPoe)
             {
                 occ++;
                 if(unique_only && hit->a.a[i].id == 0) continue;
@@ -2590,18 +2590,18 @@ void update_h_w(h_w_t *e, dens_idx_t *idx, double *max_div)
             ii = 0;
             for (pi = l; pi < k; pi++)
             {
-                pos = e->a[pi].d>>32;///
+                pos = e->a[pi].d>>32;///loc of 3'-end
                 while (ii < idn)
                 {
                     if(((uint32_t)id[ii]) == pos)
                     {
-                        if(ori)
+                        if(ori)///the most left one
                         {
                             break;
                         }
                         else
                         {
-                            while (ii < idn && (((uint32_t)id[ii]) == pos))
+                            while (ii < idn && (((uint32_t)id[ii]) == pos))///the most right one
                             {
                                 ii++;
                             }
@@ -2612,7 +2612,7 @@ void update_h_w(h_w_t *e, dens_idx_t *idx, double *max_div)
                     ii++;
                 }
                 if(ii >= idn) fprintf(stderr, "ERROR-1\n");
-                e->a[pi].w += (ori? idn-ii: ii+1);
+                e->a[pi].w += (ori? idn-ii: ii+1);///the smaller the better
                 if(max_div) (*max_div) = MAX((*max_div), e->a[pi].w);
             }
             l = k;
@@ -2785,7 +2785,7 @@ void update_scg(horder_t *h, trans_col_t *t_idx)
         if(!hits->a.a[i].id) continue;//hom hits
         suid = get_hit_suid(*hits, i);
         euid = get_hit_euid(*hits, i);
-        if(suid == euid) continue;
+        if(suid == euid) continue;//same unitig
         slen = ug->u.a[suid].len;
         elen = ug->u.a[euid].len;
 
@@ -2987,7 +2987,7 @@ void get_backbone_layout(horder_t *h, sc_lay_t *sl, osg_t *lg, uint8_t *vis)
     {
         ///I guess this should be (!!(asg_arc_n(lg, k<<1)))^(!!(asg_arc_n(lg, (k<<1)+1)))?
         ///no, since asg_arc_n is at most 1
-        if((asg_arc_n(lg, k<<1))^(asg_arc_n(lg, (k<<1)+1)))
+        if((asg_arc_n(lg, k<<1))^(asg_arc_n(lg, (k<<1)+1)))///end scaffolding
         {
             v = (asg_arc_n(lg, k<<1)?(k<<1):((k<<1)+1));
             if(vis[k<<1] || vis[(k<<1)+1]) continue;
@@ -3512,7 +3512,7 @@ void generate_scaffold(ma_utg_t *su, lay_t *ly, ma_ug_t *pug, asg_t *rg)
         }
         if(i < ly->n - 2) kv_push(uint64_t, *su, (uint64_t)-1);
     }
-    if(ly->n != 2) is_circle = 0;
+    if(ly->n != 2) is_circle = 0;///ly-> == 2: single circle
 
     for (i = 0, totalLen = 0; i < su->n-1; i++)
     {
@@ -3936,7 +3936,7 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt,
     // output_hic_rtg(i_ug, h->r_g, opt, asm_opt.output_file_name);
 
     // reduce_hamming_error(h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz);
-    reduce_hamming_error_adv(NULL, h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz, opt->ruIndex, NULL);
+    reduce_hamming_error_adv(NULL, h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz, opt->ruIndex, NULL, (asm_opt.max_short_tip*2));
     /**
     scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, FATHER);
     scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, MOTHER);

htab.cpp

@@ -385,7 +385,7 @@ ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread, int is_l)
 	ha_ct_destroy_bf(ct);
 	CALLOC(pt, 1);
 	pt->k = ct->k, pt->pre = ct->pre, pt->tot = ct->tot;
-	CALLOC(pt->h, 1<<pt->pre);
+	CALLOC(pt->h, (((uint64_t)1)<<pt->pre));
 	for (i = 0; i < 1<<pt->pre; ++i) {
 		pt->h[i].h = yak_pt_init();
 		yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h));
@@ -422,7 +422,7 @@ ha_pt_t *ha_pt_gen_count(ha_ct_t *ct, int n_thread)
 	ha_ct_destroy_bf(ct);
 	CALLOC(pt, 1);
 	pt->k = ct->k, pt->pre = ct->pre, pt->tot = ct->tot;
-	CALLOC(pt->h, 1<<pt->pre);
+	CALLOC(pt->h, (((uint64_t)1)<<pt->pre));
 	for (i = 0; i < 1<<pt->pre; ++i) {
 		pt->h[i].h = yak_pt_init();
 		yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h));
@@ -546,6 +546,21 @@ const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash)
 	return kh_key(g->h, k) & YAK_MAX_COUNT;
 }
 
+inline uint64_t flt_quals(char *sc_a, uint64_t sc_l, uint64_t sc_off, int64_t sc_cut)
+{
+	int64_t sc_min = sc_l * sc_cut, sc_tot; uint64_t k;
+	for (k = sc_tot = 0; (k < sc_l) && (sc_tot < sc_min); k++) {
+		sc_tot += (((uint8_t)sc_a[k]) - sc_off);
+	}
+
+	// if(sc_tot < sc_min) {
+	// 	fprintf(stderr, "[M::%s] sc_tot::%ld, sc_min::%ld, sc_l::%lu\n", __func__, sc_tot, sc_min, sc_l);
+	// }
+
+	if(sc_tot < sc_min) return 0;
+	return 1;
+}
+
 /**********************************
  * Buffer for counting all k-mers *
  **********************************/
@@ -563,7 +578,7 @@ KSEQ_INIT(gzFile, gzread)
 typedef struct { // global data structure for kt_pipeline()
 	const yak_copt_t *opt;
 	const void *flt_tab;
-	int flag, create_new, is_store, uq;
+	int flag, create_new, is_store, uq, ifq;
 	uint64_t n_mz, n_seq; ///number of total reads
 	kseq_t *ks;
 	UC_Read ucr;
@@ -692,6 +707,7 @@ static inline void sf##_pt_insert_buf(sf##_ch_buf_t *buf, int p, const HType *y)
 static void *sf##_worker_count(void *data, int step, void *in) /** callback for kt_pipeline()**/\
 {\
 	pl_data_t *p = (pl_data_t*)data;\
+	/**uint8_t src_a[1000000], des_a[1000000];**/\
 	if (step == 0) { /** step 1: read a block of sequences**/\
 		int ret;\
 		sf##_st_data_t *s;\
@@ -744,7 +760,8 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
 		} else {\
 			while ((ret = kseq_read(p->ks)) >= 0) {\
 				int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);\
-				if(l <= 0) continue;\
+				if((l <= 0) || (l < asm_opt.rl_cut)) continue;\
+				if((p->ifq) && (asm_opt.sc_cut > 0) && (!flt_quals(p->ks->qual.s+p->opt->adaLen, l, 33, asm_opt.sc_cut))) continue;\
 				if (p->n_seq >= 1<<28) {\
 					fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);\
 					exit(1);\
@@ -762,6 +779,16 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
 								++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);\
+						if(p->ifq) {\
+							ha_compress_qual(Get_QUAL(*p->rs_out, p->n_seq), p->ks->qual.s+p->opt->adaLen, l, sc_bn, 33);\
+							/**print_fastq(NULL, p->ks->name.s, p->ks->seq.s, p->ks->qual.s, (1<<sc_bn), 33);**/\
+							/**if(l <= 1000000) {\
+								convert_qual(src_a, p->ks->qual.s+p->opt->adaLen, l, (1<<sc_bn), 0, 33);\
+								retrive_bqual(NULL, des_a, p->n_seq, -1, -1, 0, sc_bn);\
+								if(memcmp(src_a, des_a, l)!=0) fprintf(stderr, "ERROR: incorrect qual values\n");\
+								else fprintf(stderr, "Correct: correct qual values\n");\
+							}**/\
+						}\
 					}\
 				}\
 				if (s->n_seq == s->m_seq) {\
@@ -809,7 +836,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
 			uint32_t j;\
 			/**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), CALLOC(s->mt, p->opt->n_thread);\
+			CALLOC(s->mz, s->n_seq); CALLOC(s->mz_buf, p->opt->n_thread); CALLOC(s->mt, p->opt->n_thread);\
 			/**calculate minimzers for each read, each read corresponds to one thread**/\
 			kt_for(p->opt->n_thread, sf##_worker_for_mz, s, s->n_seq);\
 			for (i = 0; i < p->opt->n_thread; ++i) free(s->mt[i].a), free(s->mz_buf[i].a);\
@@ -900,7 +927,7 @@ void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs)
 	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, ma_utg_v *us, int64_t *n_seq)
+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, ma_utg_v *us, int64_t *n_seq, uint8_t ifq)
 {
 	///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW
 	int read_rs = (rs && (flag & HAF_RS_READ));
@@ -908,7 +935,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 	pl_data_t pl;
 	gzFile fp = 0;
 	memset(&pl, 0, sizeof(pl_data_t));
-	pl.n_seq = *n_seq; 
+	pl.n_seq = *n_seq; pl.ifq = ifq;
 	if(ug_rs) {
 		pl.us_in = us;
 	} else if (read_rs) {
@@ -980,11 +1007,38 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, h
 	opt.adaLen = (keep_adapter? asm_o->adapterLen : 0);
 	opt.min_rcnt = (low_freq?*low_freq:-1);
 	opt.uq = (unique_only?1:0);
-	///asm_opt->num_reads is the number of fastq files
+	
-	for (i = n_bs = 0; i < (us?1:asm_o->num_reads); ++i){
+	n_bs = 0;
-		h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq);
+	///pending for integration
+	/**
+	if(rs && asm_o->ar && asm_o->ul_mod) {
+		for (i = 0; i < asm_o->ar->n; ++i){
+			h = yak_count(&opt, asm_o->ar->a[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq);
 			if(h) n_bs += h->bs;
 		}
+	}**/
+	///asm_opt->num_reads is the number of fastq files
+	for (i = 0; i < (us?1:asm_o->num_reads); ++i){
+		h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq, asm_opt.is_sc);
+		if(h) n_bs += h->bs;
+	}
+
+	if((rs) && (flag & HAF_RS_WRITE_LEN) && (asm_opt.is_sc)) {
+		rs->tqn = rs->total_reads; rs->tr[0] = rs->total_reads_bases;
+	}
+
+	if((asm_o->hf) && (!us)) {
+		for (i = 0; i < asm_o->hf->n; ++i) {
+			h = yak_count(&opt, asm_o->hf->a[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq, 0);
+			if(h) n_bs += h->bs;
+		}
+	}
+	if((rs) && (flag & HAF_RS_WRITE_LEN) && (asm_opt.is_sc)) {
+		rs->tr[1] = rs->total_reads_bases - rs->tr[0];
+	}
+
+	// fprintf(stderr, "[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]);
+
 	if(h) h->bs = n_bs;	
 	if (h && opt.bf_shift > 0)
 		ha_ct_destroy_bf(h);
@@ -1066,13 +1120,26 @@ void debug_ct_index(void* q_ct_idx, void* r_ct_idx)
 /*************************
  * High-level interfaces *
  *************************/
-void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff)
+
+ int64_t ha_ct_ug_cutoff(ha_ct_t *h, int64_t num_thre, double cut_rate)
 {
-	yak_ft_t *flt_tab;
+	int64_t cnt[YAK_N_COUNTS], k, tot_n = 0, tot_cutn = 0;
-    ha_ct_t *h;
+	ha_ct_hist(h, cnt, num_thre);
+	for (k = tot_n = 0; k < YAK_N_COUNTS; k++) tot_n += cnt[k]; 
+	tot_cutn = tot_n - (tot_n*cut_rate);
+
+	for (k = tot_n = 0; k < YAK_N_COUNTS && tot_n < tot_cutn; k++) tot_n += cnt[k]; 
+	return k;
+}
+
+void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff, int max_cutoff)
+{
+	yak_ft_t *flt_tab; ha_ct_t *h;
 	///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer
 	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ, !(asm_opt->flag&HA_F_NO_HPC), k, w, NULL, NULL, NULL, us, 0, NULL, 0);
 
+	if(cutoff < 0) cutoff = ha_ct_ug_cutoff(h, asm_opt->thread_num, 0.0002);
+	if(cutoff > max_cutoff) cutoff = max_cutoff;
 	// cutoff = (int)(asm_opt->hom_cov * asm_opt->high_factor);
     if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
 	// fprintf(stderr, "[M::%s::] cutoff->%d\n\n", __func__, cutoff);
@@ -1084,27 +1151,35 @@ void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int
 	return (void*)flt_tab;
 }
 
-void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq)
+void *ha_ft_ug_gen(hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq)
 {
 	yak_ft_t *flt_tab; 
-	ha_ct_t *h;
+	ha_ct_t *h; ///int32_t b0; 
+	// if(min_freq < 0 || max_freq < 0) {
+	// 	b0 = asm_opt->bf_shift; asm_opt->bf_shift = 0;
+	// }
 	///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer
 	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL, 0);
+	// if(min_freq < 0 || max_freq < 0) {
+	// 	asm_opt->bf_shift = b0;
+		
+	// }
+
 	ha_ct_shrink(h, min_freq, max_freq>YAK_MAX_COUNT-1?YAK_MAX_COUNT-1:max_freq, asm_opt->thread_num);
 	flt_tab = gen_hh(h, YAK_MAX_COUNT);
 	ha_ct_destroy(h);
 	return (void*)flt_tab;
 }
 
-void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode)
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode, int read_from_store)
 {
 	yak_ft_t *flt_tab;
 	int64_t cnt[YAK_N_COUNTS];
 	int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
 	if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
 	ha_ct_t *h; 
-	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL, 0);
+	h = ha_count(asm_opt, HAF_COUNT_ALL|ex_flag|((read_from_store)?(HAF_RS_READ):(HAF_RS_WRITE_LEN)), !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL, 0);
-	if((asm_opt->flag & HA_F_VERBOSE_GFA))
+	if((asm_opt->flag & HA_F_VERBOSE_GFA) || (asm_opt->restart))
 	{
 		write_ct_index((void*)h, asm_opt->output_file_name);
 		// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
@@ -1117,10 +1192,20 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
 	{
 		ha_ct_hist(h, cnt, asm_opt->thread_num);
 		peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, asm_opt->hg_size>0?(h->bs/asm_opt->hg_size):(-1), cnt, &peak_het);
+		///r850
+		fprintf(stderr, "[M::%s::auto] inferred peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+		if((asm_opt->het_cov_ss > 0) || (asm_opt->hmo_cov_ss > 0)) {
+			fprintf(stderr, "[M::%s::user] override requested peak_hom: %ld; peak_het: %ld\n", __func__, asm_opt->hmo_cov_ss, asm_opt->het_cov_ss);
+			if(asm_opt->hmo_cov_ss > 0) peak_hom = asm_opt->hmo_cov_ss;
+			if(asm_opt->het_cov_ss > 0) peak_het = asm_opt->het_cov_ss;
+		}
 		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);
+		if (peak_hom > 0) fprintf(stderr, "[M::%s::final] using peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+
 		///in default, asm_opt->high_factor = 5.0
+		///r833
 		cutoff = (int)(peak_hom * asm_opt->high_factor);
+		if(cutoff < asm_opt->hf_cutoff) cutoff = asm_opt->hf_cutoff;
 		if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
 	}
 	ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
@@ -1215,9 +1300,16 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	ha_ct_hist(ct, cnt, asm_opt->thread_num);
 	fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
 	peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, asm_opt->hg_size>0?(ct->bs/asm_opt->hg_size):(-1), cnt, &peak_het);
+	///r850
+	fprintf(stderr, "[M::%s::auto] inferred peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+	if((asm_opt->het_cov_ss > 0) || (asm_opt->hmo_cov_ss > 0)) {
+		fprintf(stderr, "[M::%s::user] override requested peak_hom: %ld; peak_het: %ld\n", __func__, asm_opt->hmo_cov_ss, asm_opt->het_cov_ss);
+		if(asm_opt->hmo_cov_ss > 0) peak_hom = asm_opt->hmo_cov_ss;
+		if(asm_opt->het_cov_ss > 0) peak_het = asm_opt->het_cov_ss;
+	}
 	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);
+	if (peak_hom > 0) fprintf(stderr, "[M::%s::final] using 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);
@@ -1329,8 +1421,9 @@ int load_ct_index(void **i_ct_idx, char* file_name)
 
 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);
+    char* gfa_name = (char*)malloc(strlen(file_name)+64);
-    sprintf(gfa_name, "%s.pt_flt", file_name);
+    if(r) sprintf(gfa_name, "%s.pt_flt", file_name);
+	else sprintf(gfa_name, "%s.pt_flt.bin", file_name);
     FILE* fp = fopen(gfa_name, "w");
 	if (!fp) {
 		free(gfa_name);
@@ -1369,19 +1462,36 @@ int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t*
 	fwrite(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
 	fwrite(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
 
-
+	if(r) {
 		write_All_reads(r, gfa_name);
 
+		sprintf(gfa_name, "%s.pt_flt.paf.bin", file_name); 
+		fclose(fp); fp = fopen(gfa_name, "w"); uint64_t k;
+		if (!fp) {
+			free(gfa_name);
+			return 0;
+		}
+		fwrite(&(r->total_reads), sizeof(r->total_reads), 1, fp);
+		for (k = 0; k < r->total_reads; k++) {
+			fwrite(&(r->paf[k].is_fully_corrected), sizeof(r->paf[k].is_fully_corrected), 1, fp);
+			fwrite(&(r->paf[k].is_abnormal), sizeof(r->paf[k].is_abnormal), 1, fp);
+			fwrite(&(r->paf[k].length), sizeof(r->paf[k].length), 1, fp);
+			fwrite(r->paf[k].buffer, sizeof((*(r->paf[k].buffer))), r->paf[k].length, fp);
+		}
+	}
+
 	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)
+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);
+	char* gfa_name = (char*)malloc(strlen(file_name)+64);
-    sprintf(gfa_name, "%s.pt_flt", file_name);
+	if(r) sprintf(gfa_name, "%s.pt_flt", file_name);
+	else sprintf(gfa_name, "%s.pt_flt.bin", file_name);
+	// fprintf(stderr, "[M::%s]\tgfa_name::%s\n", __func__, gfa_name);
     FILE* fp = fopen(gfa_name, "r");
 	if (!fp) {
 		free(gfa_name);
@@ -1460,25 +1570,107 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_op
 	f_flag += fread(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
 
 
+	// fclose(fp);
+	if(r) {
+		if(!load_All_reads(r, gfa_name)) {
+			free(gfa_name);
+			return 0;
+		}
+
+		memset(r->trio_flag, AMBIGU, r->total_reads*sizeof(uint8_t));
+
+		sprintf(gfa_name, "%s.pt_flt.paf.bin", file_name);
+		fclose(fp); fp = fopen(gfa_name, "r"); uint64_t k;
+		if (!fp) {
+			free(gfa_name);
+			return 0;
+		}
+		f_flag += fread(&(r->total_reads), sizeof(r->total_reads), 1, fp);
+		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 (k = 0; k < r->total_reads; k++) {
+			// init_ma_hit_t_alloc(&(r->paf[k]));
+			init_ma_hit_t_alloc(&(r->reverse_paf[k]));
+
+			f_flag += fread(&(r->paf[k].is_fully_corrected), sizeof(r->paf[k].is_fully_corrected), 1, fp);
+			f_flag += fread(&(r->paf[k].is_abnormal), sizeof(r->paf[k].is_abnormal), 1, fp);
+			f_flag += fread(&(r->paf[k].length), sizeof(r->paf[k].length), 1, fp);
+			r->paf[k].size = r->paf[k].length;
+
+			r->paf[k].buffer = NULL;
+			if(r->paf[k].length == 0) continue;
+			
+			r->paf[k].buffer = (ma_hit_t*)malloc(sizeof(ma_hit_t)*r->paf[k].length);
+			fread(r->paf[k].buffer, sizeof((*(r->paf[k].buffer))), r->paf[k].length, fp);
+		}
+	}
 	fclose(fp);
 
-	if(!load_All_reads(r, gfa_name))
+	fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
-	{
+
+	free(gfa_name);
+	return 1;
+}
+
+void refresh_pt_idx(void **flt_tab, ha_pt_t **ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint8_t is_w)
+{
+	char* gfa_name = (char*)malloc(strlen(file_name)+64);
+    sprintf(gfa_name, "%s.ad", file_name);
+
+	if(is_w) {
+		write_pt_index(*flt_tab, *ha_idx, NULL, opt, gfa_name);
+	} else {
+		load_pt_index(flt_tab, ha_idx, NULL, opt, gfa_name);
+	}
+
+	free(gfa_name);
+}
+
+
+uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load)
+{
+	char* gfa_name = (char*)malloc(strlen(file_name)+64);
+	FILE *fp = NULL; int f_flag = 0; uint64_t rr0 = (uint64_t)-1, tot_rr0 = (uint64_t)-1;
+
+
+	if(is_load) {
+		sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr);
+		fp = fopen(gfa_name, "r");
+		if (!fp) {
+			free(gfa_name);
+			return 0;
+		}
+		f_flag += fread(&rr0, sizeof(rr0), 1, fp);
+		f_flag += fread(&tot_rr0, sizeof(tot_rr0), 1, fp);
+		fclose(fp);
+		if(rr0 != rr || tot_rr0 != tot_rr) {
 			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__);
+		sprintf(gfa_name, "%s.r%lu", file_name, rr);
+		if(!load_pt_index(r_flt_tab, r_ha_idx, r, opt, gfa_name)) {
+			free(gfa_name); 
+			return 0;
+		}
+	} else {
+		sprintf(gfa_name, "%s.r%lu", file_name, rr);
+		write_pt_index(*r_flt_tab, *r_ha_idx, r, opt, gfa_name);
+
+
+
+		sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr);
+		fp = fopen(gfa_name, "w");
+		if (!fp) {
+			free(gfa_name);
+			return 0;
+		}
+		fwrite(&rr, sizeof(rr), 1, fp);
+		fwrite(&tot_rr, sizeof(tot_rr), 1, fp);
+		fclose(fp);
+	}
 
 	free(gfa_name); 
 	return 1;

htab.h

@@ -72,9 +72,9 @@ extern void *ha_flt_tab_hp;
 extern ha_pt_t *ha_idx_hp;
 extern void *ha_ct_table;
 
-void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff);
+void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff, int max_cutoff);
-void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq);
+void *ha_ft_ug_gen(hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq);
-void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode, int read_from_store);
 int32_t ha_ft_cnt(const void *hh, uint64_t y);
 void ha_ft_destroy(void *h);
 
@@ -88,11 +88,13 @@ const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash);
 
 int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
 int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
+void refresh_pt_idx(void **flt_tab, ha_pt_t **ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint8_t is_w);
 int uidx_write(void *flt_tab, ha_pt_t *ha_idx, char* file_name, ma_ug_t *ug);
 int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name, ma_ug_t *ug);
 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);
+uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load);
 
 ha_abuf_t *ha_abuf_init_buf(void *km);
 ha_abufl_t *ha_abufl_init_buf(void *km);
@@ -108,6 +110,7 @@ uint64_t ha_abufl_mem(const ha_abufl_t *ab);
 
 double yak_cputime(void);
 void yak_reset_realtime(void);
+double yak_realtime_0(void);
 double yak_realtime(void);
 long yak_peakrss(void);
 double yak_peakrss_in_gb(void);
@@ -116,6 +119,7 @@ double yak_cpu_usage(void);
 void ha_triobin(const hifiasm_opt_t *opt);
 uint32_t test_yak_binning(char* fn, char *cmd);
 uint32_t *ha_polybin_list(const hifiasm_opt_t *opt);
+uint32_t *ha_charbin_list(const hifiasm_opt_t *opt, uint8_t **idx, uint32_t *idx_n);
 
 void mz1_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km);
 void mz2_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mzl_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km);
@@ -124,6 +128,7 @@ int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_
 void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt);
 void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);
 
+
 inline int mz_low_b(int peak_hom, int peak_het)
 {
 	int low_freq = 2;
@@ -164,7 +169,7 @@ static inline uint64_t yak_hash_long(uint64_t x[4])
 	return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
 }
 
-#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr))))
+#define CALLOC(ptr, len) ((ptr) = ((((len)*sizeof(*(ptr))) <= 9223372036854775807)?((__typeof__(ptr))calloc((len), sizeof(*(ptr)))):(NULL)))
 #define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
 #define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
 #define MEMCPY(dest, src, len) (memcpy((dest), (src), (len) * sizeof(*(src))))

inter.cpp

@@ -412,6 +412,7 @@ typedef struct { // global data structure for kt_pipeline()
 	kv_u_trans_t *ta;
 	uint64_t mm;
 	uint64_t soff;
+	uint64_t is_exact;
 } ctdat_t;
 
 
@@ -513,7 +514,7 @@ double diff_ec_ul, double diff_ec_ul_low, double diff_ec_ul_hpc, int ec_ul_round
 
 void uidx_l_build(ma_ug_t *ug, mg_idxopt_t *opt, int cutoff)
 {
-    ha_flt_tab = ha_ft_ul_gen(&asm_opt, &(ug->u), opt->k, opt->w, cutoff);
+    ha_flt_tab = ha_ft_ul_gen(&asm_opt, &(ug->u), opt->k, opt->w, cutoff, -1);
     ha_idx = ha_pt_ul_gen(&asm_opt, ha_flt_tab, &(ug->u), opt->k, opt->w, cutoff);	
 	fprintf(stderr, "[M::%s] Index has been built.\n", __func__);
 }
@@ -1380,8 +1381,8 @@ const int64_t qlen, const int64_t rlen, int32_t *r_qs, int32_t *r_qe, int32_t *r
         re = rlen - 1; qe += rtail; 
     }
 
-	if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe + 1;
+	if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe + 1;}
-	if(r_rs) (*r_rs) = rs; if(r_re) (*r_re) = re + 1;
+	if(r_rs) {(*r_rs) = rs;} if(r_re) {(*r_re) = re + 1;}
     if(rev) {
         if(r_rs) (*r_rs) = rlen - re - 1; 
 		if(r_re) (*r_re) = rlen - rs;
@@ -1412,8 +1413,8 @@ int32_t *r_qs, int32_t *r_qe, int32_t *r_rs, int32_t *r_re)
         re = rlen - 1; qe += rtail; 
     }
 
-	if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe + 1;
+	if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe + 1;}
-	if(r_rs) (*r_rs) = rs; if(r_re) (*r_re) = re + 1;
+	if(r_rs) {(*r_rs) = rs;} if(r_re) {(*r_re) = re + 1;}
     if(ri->v&1) {
         if(r_rs) (*r_rs) = rlen - re - 1; 
 		if(r_re) (*r_re) = rlen - rs;
@@ -4499,7 +4500,7 @@ int64_t g_adjacent_dis(const asg_t *g, uint32_t v, uint32_t w)
 
 void get_r_offset(ma_ug_t *ug, mg_lchain_t *x, int64_t *rs, int64_t *re, int64_t *qs, int64_t *qe)
 {
-	if(qs) *qs = x->qs; if(qe) *qe = x->qe;
+	if(qs) {*qs = x->qs;} if(qe) {*qe = x->qe;}
 	if(!(x->score&1)) {
 		if(rs) *rs = x->rs + x->off;
 		if(re) *re = x->re + x->off;
@@ -4511,7 +4512,7 @@ void get_r_offset(ma_ug_t *ug, mg_lchain_t *x, int64_t *rs, int64_t *re, int64_t
 
 void get_u_offset(ma_ug_t *ug, mg_lchain_t *x, int64_t *rs, int64_t *re, int64_t *qs, int64_t *qe)
 {
-	if(qs) *qs = x->qs; if(qe) *qe = x->qe;
+	if(qs) {*qs = x->qs;} if(qe) {*qe = x->qe;}
 	if(!(x->v&1)) {
 		if(rs) *rs = x->rs + x->off;
 		if(re) *re = x->re + x->off;
@@ -4966,6 +4967,10 @@ void fill_unaligned_alignments(ma_ug_t *ug, mg_lchain_t *a, int64_t a_n, int64_t
 	}
 }
 
+void sort_uc_block_qe(uc_block_t* a, uint64_t a_n) {
+	radix_sort_uc_block_t_qe_srt(a, a + a_n); 
+}
+
 void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc, int64_t ulid)
 {
 	int64_t k, ucn = uc->n, a_n, m, l, lk; ma_ug_t *ug = uref->ug; mg_lchain_t *ix, *a; uc_block_t *z;
@@ -6725,13 +6730,13 @@ int64_t comput_sc_partial_cigar(int64_t sc, int64_t ol, double err_sc_r, overlap
 	// pk = (*wi); pe = (*werr);
 	if(k < 0) {
 		k = 0; err = 0;
-		if(wi) (*wi) = k; if(werr) (*werr) = err;
+		if(wi) {(*wi) = k;} if(werr) {(*werr) = err;}
 		// assert(e <= z->w_list.a[0].x_start);
 	} else {
 		if(z->w_list.a[k].y_end != -1) {
 			err -= z->w_list.a[k].error;
 		}
-		if(wi) (*wi) = k; if(werr) (*werr) = err;
+		if(wi) {(*wi) = k;} if(werr) {(*werr) = err;}
 	
 		// if(!(err >= 0 && k >= 0 && k < wn && z->w_list.a[k].x_start < e && z->w_list.a[k].x_end + 1 >= e)){
 		// 	fprintf(stderr, "[M::%s] ol::%ld, e::%ld, z::[%u, %u], k::%ld, wn::%ld, w::[%d, %d], err::%ld\n", __func__, 
@@ -8812,7 +8817,7 @@ double e_rate, kv_rtrace_t *trace, uint32_t rechain_w, ul_ov_t *res, ul_ov_t *rr
 	for (i = k = 0; i < nv; i++) {
 		if(res[i].qn == (uint32_t)-1) continue;
 		res[k] = res[i]; 
-		if(res[k].sec&mm) res[k].sec -= mm; k++;
+		if(res[k].sec&mm) {res[k].sec -= mm;} k++;
 		// res[k] = res[i]; 
 		// id[k] = res[k].sec; 
 		// if(id[k]&mm) id[k] -= mm; 
@@ -10094,7 +10099,7 @@ uint32_t test_het_aln(ma_ug_t *ug, uint64_t rid, u_trans_t *a, uint64_t a_n, st_
 
 uint32_t is_mmhom_node(uint64_t *ca, ma_utg_t *u, asg_t *sg, uint64_t cov_bd, double cut_rate)
 {
-	if(cut_rate < 0) cut_rate = 0; if(cut_rate > 1.0) cut_rate = 1.0;
+	if(cut_rate < 0) {cut_rate = 0;} if(cut_rate > 1.0) {cut_rate = 1.0;}
 	uint64_t k, a, na, a_cut = u->n*cut_rate, na_cut = u->n*(1.0-cut_rate);
 	for (k = a = na = 0; k < u->n; k++) {
 		if(ca[k] > (cov_bd*((uint64_t)sg->seq[u->a[k]>>33].len))) {
@@ -13632,8 +13637,8 @@ void gen_end_coord(ul_ov_t *z, int64_t qlen, int64_t tlen, int64_t *r_qs, int64_
         te = tlen; qe += ttail; 
     }
 
-	if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe;
+	if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe;}
-	if(r_ts) (*r_ts) = ts; if(r_te) (*r_te) = te;
+	if(r_ts) {(*r_ts) = ts;} if(r_te) {(*r_te) = te;}
     if(z->rev) {
         if(r_ts) (*r_ts) = tlen - te; 
 		if(r_te) (*r_te) = tlen - ts;
@@ -13817,8 +13822,8 @@ void extend_end_coord(mg_lchain_t *li, ul_ov_t *ui, const int64_t qlen, const in
         re = rlen; qe += rtail; 
     }
 
-	if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe;
+	if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe;}
-	if(r_rs) (*r_rs) = rs; if(r_re) (*r_re) = re;
+	if(r_rs) {(*r_rs) = rs;} if(r_re) {(*r_re) = re;}
     if(rev) {
         if(r_rs) (*r_rs) = rlen - re; 
 		if(r_re) (*r_re) = rlen - rs;
@@ -14322,7 +14327,7 @@ uint64_t detect_mul_way, float len_dif)
 int64_t hc_chain_backtrack(int64_t n, const int64_t *f, const uint64_t *p, uint64_t *srt, uint64_t *u, uint64_t *v, 
 int64_t *n_u_, int64_t *n_v_)
 {
-	if(n_u_) *n_u_ = 0; if(n_v_) *n_v_ = 0;
+	if(n_u_) {*n_u_ = 0;} if(n_v_) {*n_v_ = 0;}
 	int64_t i, k, n_v, n_srt, n_v0, n_u, sc;
 	if (n == 0) return 0;
 	// v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak
@@ -14346,7 +14351,7 @@ int64_t *n_u_, int64_t *n_v_)
 		u[n_u++] = (((uint64_t)sc)<<32) | ((uint64_t)(n_v-n_v0));
 	}
 
-	if(n_u_) *n_u_ = n_u; if(n_v_) *n_v_ = n_v;
+	if(n_u_) {*n_u_ = n_u;} if(n_v_) {*n_v_ = n_v;}
 	return n_u;
 }
 
@@ -15527,7 +15532,7 @@ void update_rovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
 		for (i = sidx+1; i < eidx; i++) {
 			get_r_offset(ug, &(a[i]), &rs, &re, NULL, NULL);
 			a[i].qe = cal_qext_coor(left_r[1], re, left_q[1], left_q[1] + re - left_r[1], re);
-			if(a[i].qe < 0) a[i].qe = 0; if(a[i].qe > qlen) a[i].qe = qlen;
+			if(a[i].qe < 0) {a[i].qe = 0;} if(a[i].qe > qlen) {a[i].qe = qlen;}
 			a[i].qs = cal_qext_coor(left_r[0], (a[i].qe<=left_q[1])?re:left_r[1], 
 											left_q[0], (a[i].qe<=left_q[1])?a[i].qe:left_q[1], rs);
 			assert(a[i].qs >= 0 && a[i].qs <= qlen);
@@ -15548,7 +15553,7 @@ void update_rovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
 			a[i].qe = cal_qext_coor(right_r[0], right_r[1], right_q[0], right_q[1], re);
 			assert(a[i].qe >= 0 && a[i].qe <= qlen);
 			a[i].qs = cal_qext_coor(rs, right_r[0], right_q[0]-(right_r[0]-rs), right_q[0], rs);
-			if(a[i].qs < 0) a[i].qs = 0; if(a[i].qs > qlen) a[i].qs = qlen;
+			if(a[i].qs < 0) {a[i].qs = 0;} if(a[i].qs > qlen) {a[i].qs = qlen;}
 			if(a[i].qs > a[i].qe) {
 				tt = a[i].qs; a[i].qs = a[i].qe; a[i].qe = tt;
 			}
@@ -16340,8 +16345,8 @@ void ctg_rg2ug_gen(ma_utg_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint6
 				for (i = l; i < k; i++) {///could be merged
 					m = get_unique_rctg_aln(uref, r_cl->a[i]>>32, i, ls, r_cl->len, &kp, NULL, ((uint32_t)-1)); 
 					assert(m == 1); assert(kp.tn == lp.tn); assert(kp.qn == lp.qn + i - l); assert(kp.sec == lp.sec + i - l);
-					if(kp.qs < lp.qs) lp.qs = kp.qs; if(kp.qe > lp.qe) lp.qe = kp.qe;
+					if(kp.qs < lp.qs) {lp.qs = kp.qs;} if(kp.qe > lp.qe) {lp.qe = kp.qe;}
-					if(kp.ts < lp.ts) lp.ts = kp.ts; if(kp.te > lp.te) lp.te = kp.te;
+					if(kp.ts < lp.ts) {lp.ts = kp.ts;} if(kp.te > lp.te) {lp.te = kp.te;}
 					ls += (uint32_t)(r_cl->a[i]); 
 
 					// if(id == 123) {
@@ -16802,7 +16807,7 @@ uint64_t gen_trans_ovlp_scaf(ma_ug_t *gfa, u_trans_t *map, uc_block_t *qin, uc_b
 	return 1;
 }
 
-void extract_trans_scaf_res_t(uc_block_t *qin, const ul_idx_t *uref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint32_t avid, kv_ul_ov_t *res)
+void extract_trans_scaf_res_t(uc_block_t *qin, const ul_idx_t *uref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint32_t avid, uint32_t is_exact, kv_ul_ov_t *res)
 {
 	utg_rid_dt *a; uint64_t i, a_n; uc_block_t *rin; 
 	u_trans_t *u; uint64_t k, u_n;
@@ -16812,6 +16817,7 @@ void extract_trans_scaf_res_t(uc_block_t *qin, const ul_idx_t *uref, scaf_res_t
 	u = u_trans_a(*ta, qin->hid); 
 	u_n = u_trans_n(*ta, qin->hid);
 	for (k = 0; k < u_n; k++) {
+		if(is_exact && u[k].f != RC_0 && u[k].f != RC_1) continue;
 		// fprintf(stderr, "+[M::%s] utg%.6u%c->utg%.6u%c, q::[%u, %u), %c, t::[%u, %u)\n", __func__, qin->hid+1, "lc"[gfa->u.a[qin->hid].circ], u[k].tn+1, "lc"[gfa->u.a[u[k].tn].circ], u[k].qs, u[k].qe, "+-"[u[k].rev], u[k].ts, u[k].te);
 		a = get_r_ug_region(uref->r_ug, &a_n, u[k].tn);
     	for (i = 0; i < a_n; i++) {
@@ -16827,7 +16833,7 @@ void extract_trans_scaf_res_t(uc_block_t *qin, const ul_idx_t *uref, scaf_res_t
 	}
 }
 
-void cl_trans_gen(uint32_t qid, uint32_t qlen, ul_vec_t *qstr, kv_ul_ov_t *res, const ul_idx_t *uref, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, bubble_type *bub, kv_u_trans_t *ta, uint32_t is_self)
+void cl_trans_gen(uint32_t qid, uint32_t qlen, ul_vec_t *qstr, kv_ul_ov_t *res, const ul_idx_t *uref, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, bubble_type *bub, kv_u_trans_t *ta, uint32_t is_self, uint32_t is_exact)
 {
     uint64_t k, l, m = 0, mi, z, zn, a_n, bid, nw; uc_block_t *a; uint32_t avid = ((is_self)?(qid):((uint32_t)-1));
     res->n = 0;
@@ -16882,7 +16888,7 @@ void cl_trans_gen(uint32_t qid, uint32_t qlen, ul_vec_t *qstr, kv_ul_ov_t *res,
 			if(extract_scaf_res_t(&(qstr->bb.a[z]), uref, ref_sc, avid, res)) continue;
 
 			///trans match
-			extract_trans_scaf_res_t(&(qstr->bb.a[z]), uref, ref_sc, gfa, ta, avid, res);
+			extract_trans_scaf_res_t(&(qstr->bb.a[z]), uref, ref_sc, gfa, ta, avid, is_exact, res);
 		}
 	}
 
@@ -17159,7 +17165,7 @@ double diff_ec_ul, double ovlp_max, double unmatch_max, int64_t qlen, int64_t tl
     return 1;
 }
 
-void ctg_trans_gp_chain(uint32_t qid, kv_ul_ov_t *res, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, int64_t qlen, Chain_Data* dp, ma_ug_t *ref, uint32_t is_self, asg64_v *b)
+void ctg_trans_gp_chain(uint32_t qid, kv_ul_ov_t *res, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, int64_t qlen, Chain_Data* dp, ma_ug_t *ref, uint32_t is_self, asg64_v *b, double unmatch_max)
 {
 	uint64_t k, l, m = 0, rn = res->n; ul_ov_t rr; uint32_t avid = ((is_self)?(qid):((uint32_t)-1));
 	for (k = 0; k < rn; k++) {
@@ -17170,7 +17176,7 @@ void ctg_trans_gp_chain(uint32_t qid, kv_ul_ov_t *res, const ul_idx_t *uref, con
 	for (k = 1, l = 0, b->n = 0; k <= rn; k++) {
 		if((k == rn) || ((res->a[l].tn>>1) != (res->a[k].tn>>1))) {
 			// fprintf(stderr, "+[M::%s]\tutg%.6ul\n", __func__, (res->a[l].tn>>1) + 1);
-			if(((res->a[l].tn>>1) != avid) && (linear_ctg_trans_chain_dp(qid, res->a+l, k-l, uref, uopt, bw, diff_ec_ul, /**0.333333**/0.4, /**0.333333**/0.666666, qlen, ref->u.a[(res->a[l].tn>>1)].len, UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, &rr, b))) {
+			if(((res->a[l].tn>>1) != avid) && (linear_ctg_trans_chain_dp(qid, res->a+l, k-l, uref, uopt, bw, diff_ec_ul, /**0.333333**/0.4, unmatch_max, qlen, ref->u.a[(res->a[l].tn>>1)].len, UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, &rr, b))) {
 				res->a[m++] = rr; 
 				// fprintf(stderr, "-[M::%s]\tutg%.6ul\n", __func__, rr.tn + 1);
 			}
@@ -17345,18 +17351,18 @@ void push_ctg_trans_res(uint32_t id, kv_ul_ov_t *in, kv_ul_ov_t *ou)
 }
 
 uint32_t direct_ctg_trans_chain(mg_tbuf_t *b, uint32_t id, glchain_t *ll, gdpchain_t *gdp, st_mt_t *sps, haplotype_evdience_alloc *hap, const ul_idx_t *uref, const ug_opt_t *uopt,
-int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, bubble_type *bub, kv_u_trans_t *ta, const asg_t *rg, uint64_t soff)
+int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, bubble_type *bub, kv_u_trans_t *ta, const asg_t *rg, uint64_t soff, uint64_t is_exact)
 {
 	// res->bb.n = 0;
 	// if(ulid != 86660) return 0;
 	kv_ul_ov_t *idx = &(ll->lo), *init = &(ll->tk); ///int64_t max_idx;
 	asg64_v b0, b1; uint32_t is_self = (qry?0:1); idx->n = 0;
 	
-	cl_trans_gen(id, qry?qry->u.a[id].len:ref->u.a[id].len, qry_sc?&(qry_sc->a[id]):&(ref_sc->a[id]), idx, uref, ref, ref_sc, gfa, bub, ta, is_self);
+	cl_trans_gen(id, qry?qry->u.a[id].len:ref->u.a[id].len, qry_sc?&(qry_sc->a[id]):&(ref_sc->a[id]), idx, uref, ref, ref_sc, gfa, bub, ta, is_self, is_exact);
 	if(idx->n == 0) return 0;
 
 	copy_asg_arr(b0, (*sps)); 
-	ctg_trans_gp_chain(id, idx, uref, uopt, bw, diff_ec_ul, qry?qry->u.a[id].len:ref->u.a[id].len, dp, ref, is_self, &b0);
+	ctg_trans_gp_chain(id, idx, uref, uopt, bw, diff_ec_ul, qry?qry->u.a[id].len:ref->u.a[id].len, dp, ref, is_self, &b0, ((is_exact)?(0.5):(0.666666)));
 	copy_asg_arr((*sps), b0); 
 	if(idx->n == 0) return 0;
 
@@ -17471,7 +17477,7 @@ static void worker_for_ctg_trans_alignment(void *data, long i, int tid)
 
 	s->hab[tid]->num_read_base++; 
 	s->hab[tid]->num_correct_base += direct_ctg_trans_chain(s->buf[tid], i, &(s->ll[tid]), &(s->gdp[tid]), &(s->sps[tid]), &(s->hab[tid]->hap), s->uu, s->uopt, G_CHAIN_BW, s->opt->diff_ec_ul, UG_SKIP, i, &(s->hab[tid]->clist.chainDP), 
-	c->qry, c->qry_sc, c->ref, c->ref_sc, c->gfa, c->bub, c->ta, s->rg, c->soff);
+	c->qry, c->qry_sc, c->ref, c->ref_sc, c->gfa, c->bub, c->ta, s->rg, c->soff, c->is_exact);
 }
 
 void rm_dup_aln(u_trans_t *u, uint64_t u_n, asg64_v *b, double dup_cut)
@@ -18178,7 +18184,9 @@ const ug_opt_t *uopt, const asg_t *rg, R_to_U *ri, int64_t bw, double diff_ec_ul
 {
 	uint64_t dn = dump->n, i, mqs, mqe; int64_t iqs, iqe, its, ite; idx->n = 0;
 	collect_nc_ovlps(rch, rch_i, idx, uopt, rg, ri, ulid, &mqs, &mqe);
-	assert(idx->n > 1); assert((mqs != ((uint64_t)-1)) && (mqe != ((uint64_t)-1)) && (mqe > mqs));
+	// assert(idx->n > 1); assert((mqs != ((uint64_t)-1)) && (mqe != ((uint64_t)-1)) && (mqe > mqs));
+	if((idx->n <= 1) || (mqs == ((uint64_t)-1)) || (mqe == ((uint64_t)-1)) || (mqe <= mqs)) return;
+	
 	gen_linear_rchains(idx, dump, rg, uopt, bw, diff_ec_ul, rch->rlen, dp, ulid);
 	assert(idx->n);
 	idx->n = gen_cns_chain_linear(idx->a, idx->n, /**dump->a+dn,**/ rg, ri, rch->rlen, bw, diff_ec_ul, dp, UG_SKIP_N, UG_ITER_N, UG_DIS_N, mqs, mqe, ulid);
@@ -18582,7 +18590,7 @@ void collect_pp_ovlps(ul_ov_t *a, uint64_t a_n, kv_ul_ov_t *res, const ug_opt_t
         if(((*mqe) == ((uint64_t)-1)) || ((*mqe) < a[i].qe)) (*mqe) = a[i].qe;
 		
 	}
-	if(st) (*mqs) = st->qs; if(et) (*mqe) = et->qe;
+	if(st) {(*mqs) = st->qs;} if(et) {(*mqe) = et->qe;}
 }
 
 uint64_t gen_cns_chain_linear_hard(ul_ov_t *a, int64_t a_n, const asg_t *rg, int64_t qlen, int64_t bw, double diff_thre, Chain_Data* dp, 
@@ -18972,7 +18980,7 @@ void work_ctg_path_trans(uldat_t *sl, asg_t *sg, ma_ug_t *qry, scaf_res_t *qry_s
 	fprintf(stderr, "[M::%s::] # try:%d, # done:%d\n", __func__, s.sum_len, s.n); 
 }
 
-void work_ctg_path_trans_self(uldat_t *sl, asg_t *sg, ma_ug_t *db, scaf_res_t *db_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint64_t soff, bubble_type *bu, kv_u_trans_t *res)
+void work_ctg_path_trans_self(uldat_t *sl, asg_t *sg, ma_ug_t *db, scaf_res_t *db_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint64_t soff, uint64_t is_exact, bubble_type *bu, kv_u_trans_t *res)
 {
 	utepdat_t s; uint64_t i; memset(&s, 0, sizeof(s)); 
 	s.id = 0; s.opt = sl->opt; s.ug = sl->ug; s.uopt = sl->uopt; s.rg = sl->rg; s.uu = sl->uu; 
@@ -18983,7 +18991,7 @@ void work_ctg_path_trans_self(uldat_t *sl, asg_t *sg, ma_ug_t *db, scaf_res_t *d
 		s.hab[i] = ha_ovec_init(0, 0, 1); s.buf[i] = mg_tbuf_init();
 	}
 
-	ctdat_t c; memset(&c, 0, sizeof(c)); c.soff = soff; c.bub = bu;
+	ctdat_t c; memset(&c, 0, sizeof(c)); c.soff = soff; c.bub = bu; c.is_exact = is_exact;
 	c.s = &(s); c.qry = NULL; c.qry_sc = NULL; c.ref = db; c.ref_sc = db_sc; c.gfa = gfa; c.ta = ta;
 
 	// detect_outlier_len("+++work_ul_gchains");
@@ -19002,7 +19010,6 @@ void work_ctg_path_trans_self(uldat_t *sl, asg_t *sg, ma_ug_t *db, scaf_res_t *d
 	fprintf(stderr, "[M::%s::] # try:%d, # done:%d\n", __func__, s.sum_len, s.n); 
 }
 
-
 uint64_t work_ul_gchains_consensus(uldat_t *sl)
 {
 	utepdat_t s; uint64_t i; memset(&s, 0, sizeof(s));
@@ -21761,7 +21768,7 @@ void asg_arc_push_contain_trans(asg_t *g, ma_hit_t_alloc* ov, int64_t min_ovlp,
 		nv = asg_arc_n(g, v); av = asg_arc_a(g, v); avi = g->idx[v]>>32;
 		while(nv) {
 			for (i = nc = cc = 0; i < nv; ++i) {
-				if(av[i].del) continue; w = av[i].v;
+				if(av[i].del) {continue;} w = av[i].v;
 				if(!(is_contain_r((*ri), (w>>1)))) continue;///new arcs must be bridged by contained reads
 				assert(!(g->seq[w>>1].del));
 				nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); awi = g->idx[w]>>32;
@@ -21931,17 +21938,18 @@ void gen_contig_trans(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *qry, scaf_res_t
 	destroy_ul_idx_t(uu);
 }
 
-void gen_contig_self(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *db, scaf_res_t *db_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint64_t soff, bubble_type *bu, kv_u_trans_t *res)
+void gen_contig_self(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *db, scaf_res_t *db_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint64_t soff, bubble_type *bu, kv_u_trans_t *res, uint32_t is_exact)
 {
 	mg_idxopt_t opt; uldat_t sl; int32_t cutoff; 
 	init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
 	cutoff = asm_opt.max_n_chain;
-	init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round);
+	init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, ((is_exact)?(0.333333):(asm_opt.ul_error_rate)), asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round);
 	ul_idx_t *uu = gen_ul_idx_t_sc(db, sg, db_sc, gfa->u.n);
 	init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, uu); sl.rg = sg; sl.ug = db;
 
-	work_ctg_path_trans_self(&sl, sg, db, db_sc, gfa, ta, soff, bu, res); uu->ug = NULL;
+	work_ctg_path_trans_self(&sl, sg, db, db_sc, gfa, ta, soff, is_exact, bu, res); 
-	destroy_ul_idx_t(uu);
+
+	uu->ug = NULL; destroy_ul_idx_t(uu);
 }
 
 void order_contig_trans(kv_u_trans_t *in)

inter.h

@@ -129,7 +129,8 @@ void dedup_contain_g(const ug_opt_t *uopt, asg_t *sg);
 void trans_base_mmhap_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res);
 scaf_res_t *gen_contig_path(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *ctg, ma_ug_t *ref);
 void gen_contig_trans(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint32_t qoff, uint32_t toff, bubble_type *bu, kv_u_trans_t *res);
-void gen_contig_self(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *db, scaf_res_t *db_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint64_t soff, bubble_type *bu, kv_u_trans_t *res);
+void gen_contig_self(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *db, scaf_res_t *db_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint64_t soff, bubble_type *bu, kv_u_trans_t *res, uint32_t is_exact);
 void order_contig_trans(kv_u_trans_t *in);
+void sort_uc_block_qe(uc_block_t* a, uint64_t a_n);
 
 #endif

main.cpp

@@ -8,6 +8,7 @@
 
 int main(int argc, char *argv[])
 {
+	// setvbuf(stderr, NULL, _IONBF, 0);
 	int i, ret;
 	yak_reset_realtime();
     init_opt(&asm_opt);
@@ -61,14 +62,41 @@ int main(int argc, char *argv[])
 	// ed_band_cal_global((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3),
 	// ed_band_cal_global_128bit((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3));
 	// exit(1);
-	if(!(asm_opt.dbg_ovec_cal)) ret = ha_assemble();
+	
-	else ret = ha_assemble_ovec();
+	int8_t simd_auto = 0;
+	#if defined(__x86_64__) || defined(__i386__)
+		__builtin_cpu_init();
+		if (__builtin_cpu_supports("avx512f")) simd_auto = 2;
+		else if (__builtin_cpu_supports("avx2")) simd_auto = 1;
+		else simd_auto = 0;
+	#endif
+
+	if (simd_auto) {
+		fprintf(stderr, "[M::%s::auto] detected CPU support for %s\n", __func__, (simd_auto == 2) ? "AVX-512" : "AVX2");
+	} else {
+		fprintf(stderr, "[M::%s::auto] no supported SIMD extension detected; falling back to non-SIMD\n", __func__);
+	}
+	
+	if (asm_opt.simd_mm == 0 || asm_opt.simd_mm == 1 || asm_opt.simd_mm == 2) {
+		fprintf(stderr, "[M::%s::user] user requested %s mode\n", __func__, (asm_opt.simd_mm == 2) ? "AVX-512" : ((asm_opt.simd_mm == 1) ? "AVX2" : "non-SIMD"));
+		if(asm_opt.simd_mm > simd_auto) asm_opt.simd_mm = simd_auto;
+	} else {
+		asm_opt.simd_mm = simd_auto; 
+		if(asm_opt.simd_mm >= 1) asm_opt.simd_mm = 1; ///use avx2 rather than avx512, looks like avx512 still has issues right now
+	}
+	
+	fprintf(stderr, "[M::%s::final] using %s mode\n", __func__, (asm_opt.simd_mm == 2) ? "AVX-512" : ((asm_opt.simd_mm == 1) ? "AVX2" : "non-SIMD"));
+
+
+	if(asm_opt.sec_in) ret = ha_assemble_pair();
+	else if(asm_opt.dbg_ovec_cal) ret = ha_ec_dbg();
+	else ret = ha_assemble();
 	
     destory_opt(&asm_opt);
 	fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION);
 	fprintf(stderr, "[M::%s] CMD:", __func__);
 	for (i = 0; i < argc; ++i)
 		fprintf(stderr, " %s", argv[i]);
-	fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss_in_gb());
+	fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB; SIMD: %u\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss_in_gb(), asm_opt.simd_mm);
     return ret;
 }

rcut.cpp

@@ -806,8 +806,19 @@ uint32_t mc_edges_symm(mc_match_t *ma)
 
 static void normalize_mb_edge(mb_edge_t *a, mb_edge_t *b)
 {
-	if(a->w >= b->w)
+	uint8_t f = 0;
-	{
+	if(a->w[0] > b->w[0]) {
+		f = 1;
+	} else if(a->w[0] == b->w[0] && a->w[1] > b->w[1]) {
+		f = 1;
+	} else if(a->w[0] == b->w[0] && a->w[1] == b->w[1] && a->w[2] > b->w[2]) {
+		f = 1;
+	} else if(a->w[0] == b->w[0] && a->w[1] == b->w[1] && a->w[2] == b->w[2] && a->w[3] > b->w[3]) {
+		f = 1;
+	}
+
+	// if(a->w >= b->w) {
+	if(f) {
 		b->x = (uint32_t)a->x;
 		b->x <<= 32;
 		b->x |= (a->x>>32);
@@ -815,9 +826,7 @@ static void normalize_mb_edge(mb_edge_t *a, mb_edge_t *b)
 		b->w[3] = a->w[3];
 		b->w[1] = a->w[2];
 		b->w[2] = a->w[1];
-	}
+	} else {
-	else
-	{
 		a->x = (uint32_t)b->x;
 		a->x <<= 32;
 		a->x |= (b->x>>32);
@@ -3346,7 +3355,7 @@ mc_clus_t *init_mc_clus_t(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub, uin
 	mc_clus_t *p; CALLOC(p, 1); 
 	p->bub = bub; p->opt = opt; p->mg = mg; p->n = bub->ug->g->n_seq;
 	CALLOC(p->lock, p->n);
-	if(n_thread > 64) n_thread = 64; p->n_thread = n_thread;
+	if(n_thread > 64) {n_thread = 64;} p->n_thread = n_thread;
 	CALLOC(p->aux, p->n_thread);
 	uint32_t k, ss = (p->n>>3)+(!!(p->n&7));
 	for (k = 0; k < p->n_thread; k++) {

sys.cpp

@@ -31,6 +31,12 @@ double yak_realtime(void)
 	return yak_realtime_core() - yak_realtime0;
 }
 
+double yak_realtime_0(void)
+{
+	return yak_realtime_core();
+}
+
+
 long yak_peakrss(void)
 {
 	struct rusage r;