diff --git a/CommandLines.cpp b/CommandLines.cpp
index 9c6833e..06610b1 100644
--- a/CommandLines.cpp
+++ b/CommandLines.cpp
@@ -90,6 +90,8 @@ static ko_longopt_t long_options[] = {
     { "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},
     // { "path-round",     ko_required_argument, 348},
 	{ 0, 0, 0 }
 };
@@ -133,6 +135,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
     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, "  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);
@@ -251,7 +254,13 @@ void Print_H(hifiasm_opt_t* asm_opt)
     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   file names of HiFi reads\n");
+    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");
@@ -420,6 +429,10 @@ void init_opt(hifiasm_opt_t* asm_opt)
 
     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 = 128;
 }   
 
 void destory_enzyme(enzyme* f)
@@ -805,6 +818,11 @@ int check_option(hifiasm_opt_t* asm_opt)
         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;
 }
 
@@ -1084,6 +1102,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
             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 == 'l') {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
             asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
         }
diff --git a/CommandLines.h b/CommandLines.h
index a3a9893..cd8f027 100644
--- a/CommandLines.h
+++ b/CommandLines.h
@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.25.0-r861"
+#define HA_VERSION "0.25.0-r866"
 
 #define VERBOSE 0
 
@@ -148,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;
@@ -195,6 +196,9 @@ typedef struct {
 
     char *dbg_run_1, *dbg_run_2;
 
+    int32_t hyb_syn;
+
+    int64_t step_rd;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;
diff --git a/Correct.cpp b/Correct.cpp
index 0ce20c7..0f0ff6e 100644
--- a/Correct.cpp
+++ b/Correct.cpp
@@ -8845,7 +8845,8 @@ void prt_sub_cigar(overlap_region* z, uint64_t str_l, uint64_t site, uint64_t wi
 
 #define is_st_bs(s, rr, mm) (((mm) != ((uint64_t)-1)) && (((s).overlap_num + mm) >= ((s).occ_0)) && ((((s).occ_0*(rr) + (s).overlap_num)) >= ((s).occ_0)))
 
-void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, uint64_t multi_check, double st_rate, uint64_t st_max)
+void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, uint64_t multi_check, double st_rate, uint64_t st_max, 
+    int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     // fprintf(stderr, "[M::%s::] Done\n", __func__);
     if(hap->length == 0) return;
@@ -8877,18 +8878,6 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
     hap->snp_stat.n = m_snp_stat; hap->length = m_list;
     if(hap->snp_stat.n == 0 || hap->length == 0) return;
 
-    // for (k = 0; k < hap->snp_stat.n; k++) {
-    //     s = &(hap->snp_stat.a[k]); if(s->site != 1502) continue;
-    //     fprintf(stderr, "[M::%s] site::%u, occ_0::%u, occ_1::%u, occ_2::%u, is_homopolymer::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_2, s->is_homopolymer);
-    //     prt_sub_read(g_read->seq, g_read->length, s->site, 25);
-    //     for (i = 0; i < overlap_list->length; i++) {
-    //         if(/**overlap_list->list[i].is_match == 1 &&**/ overlap_list->list[i].x_pos_s <= s->site && s->site <= overlap_list->list[i].x_pos_e) {
-    //             fprintf(stderr, "%.*s\tis_match::%u\tid::%u\n", (int)Get_NAME_LENGTH(R_INF, overlap_list->list[i].y_id), Get_NAME(R_INF, overlap_list->list[i].y_id), overlap_list->list[i].is_match, overlap_list->list[i].y_id);
-    //         }   
-    //     }
-    // }
-    
-
     hap->snp_srt.n = 0;
     radix_sort_haplotype_evdience_id_srt(hap->list, hap->list + hap->length);
     for (k = 1, l = 0; k <= hap->length; ++k) {  
@@ -8899,53 +8888,11 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
                 assert(s->site == hap->list[i].site);
                 if(s->occ_0 < 2 || s->occ_1 < 2) continue;
                 if(is_st_bs((*s), st_rate, st_max)) continue;
-                if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
-                    // fprintf(stderr, "[M::%s] site::%u, occ_0::%u, occ_1::%u, occ_2::%u, is_homopolymer::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_1, s->is_homopolymer);
+                if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) {
                     o++;///allels must be real
                 }
             }
 
-            // if(overlap_list->list[hap->list[l].overlapID].y_id == 1740) {
-            //     fprintf(stderr, "str1074[M::%s-id::%u] o->%lu(%c), l::%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o, "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], l);
-            //     for (i = l; i < k; i++) {
-            //         if(hh_tp(hap->list[i])!=1) continue;///mismatch
-            //         s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
-            //         assert(s->site == hap->list[i].site);
-            //         if(s->occ_0 < 2 || s->occ_1 < 2) continue;
-            //         if(is_st_bs((*s), st_rate, st_max)) continue;
-            //         if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
-            //             fprintf(stderr, "[M::%s] site::%u, occ_0::%u, occ_1::%u, occ_2::%u, is_homopolymer::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_1, s->is_homopolymer);
-            //             // prt_sub_read(g_read->seq, g_read->length, s->site, 50);
-            //         }
-            //     }
-            // } 
-            // else {
-            //     fprintf(stderr, "[M::%s-id::%u] o->%lu(%c), l::%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o, "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], l);
-            // } 
-            // else {
-            //     for (i = l, o = 0; i < k; i++) {
-            //         if(hh_tp(hap->list[i])!=0) continue;///mismatch
-            //         s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
-            //         assert(s->site == hap->list[i].site);
-            //         if(s->occ_0 < 2 || s->occ_1 < 2) continue;
-            //         if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
-            //             if(s->site == 7878 || s->site == 9682) {
-                        
-            //                 fprintf(stderr, "***[M::%s] site::%u, rid::%u\t%.*s\n", __func__, s->site, overlap_list->list[hap->list[l].overlapID].y_id, (int)Get_NAME_LENGTH(R_INF, overlap_list->list[hap->list[l].overlapID].y_id), 
-            //                     Get_NAME(R_INF, overlap_list->list[hap->list[l].overlapID].y_id));
-            //                 prt_sub_cigar(&(overlap_list->list[hap->list[l].overlapID]), g_read->length, s->site, 50);
-            //             }
-            //             // o++;///allels must be real
-            //             // fprintf(stderr, "[M::%s] site::%u, occ_0::%u, occ_1::%u, occ_2::%u, is_homopolymer::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_1, s->is_homopolymer);
-            //             // prt_sub_read(g_read->seq, g_read->length, s->site, 50);
-            //         }
-            //     }
-            // }
-
-            // if(overlap_list->list[hap->list[l].overlapID].y_id == 3626) {
-            //     fprintf(stderr, "***0***[M::%s-id::%u] o->%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o);
-            // }
-
             if(o > 0) {
                 o = ((uint32_t)-1) - o;
                 o <<= 32; o += l; 
@@ -8954,8 +8901,6 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
             l = k;
         }
     }
-
-    // fprintf(stderr, "[M::%s] snp_srt.n->%lu\n", __func__, ((uint64_t)hap->snp_srt.n));
     
     if (hap->snp_srt.n > 0) {
         radix_sort_bc64(hap->snp_srt.a, hap->snp_srt.a + hap->snp_srt.n);///sort by how many snps in one overlap
@@ -8966,16 +8911,10 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
                 s = &(hap->snp_stat.a[hap->list[i].overlapSite]);
                 if(s->occ_0 < 2 || s->occ_1 < 2) continue;
                 if(is_st_bs((*s), st_rate, st_max)) continue;
-                if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
+                if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) {
                     o++;
-                    // if(overlap_list->list[hap->list[l].overlapID].y_id == 20835) {
-                    //     fprintf(stderr, "[M::%s-id::%u] occ_0->%u, occ_1->%u, occ_2->%u, site->%u\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, s->occ_0, s->occ_1, s->occ_2, s->site);
-                    // }
                 }
             }
-            // if(overlap_list->list[hap->list[l].overlapID].y_id == 1740) {
-            //  fprintf(stderr, "srt[M::%s-id::%u] o->%lu, l->%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o, l);
-            // }
             if(o == 0) continue;
 
             ii = hap->list[l].overlapID;
@@ -8991,9 +8930,6 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
                     for (z = hap->list[i].overlapSite; z >= 0; z--) {
                         t = &(hap->snp_stat.a[z]);
                         if(s->site!=t->site) break;
-                        // if(t->site == 14217) {
-                        //     fprintf(stderr, "[M::%s]\tsite::%u\tn0::%u\tn1::%u\to::%lu\t%.*s\n", __func__, t->site, t->occ_0, t->occ_1, o, (int)Get_NAME_LENGTH(R_INF, overlap_list->list[ii].y_id), Get_NAME(R_INF, overlap_list->list[ii].y_id));
-                        // }
                         t->occ_0 -= hap->list[i].cov;
                         assert(t->occ_0 >= 1);
                         if((st_max != ((uint64_t)-1)) && (overlap_list->list[ii].y_pos_strand == 0)) {
@@ -9015,10 +8951,6 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
                 if(s->score == 1) o++;
             }
             ii = hap->list[l].overlapID;
-            ///for HiFi, do not flip trans to cis
-            // if(overlap_list->list[ii].is_match == 2 && o == 0) {
-            //     overlap_list->list[ii].is_match = 1;
-            // }
             if(overlap_list->list[ii].is_match == 1 && o > 0) {
                 overlap_list->list[ii].is_match = 2;
             }
@@ -9054,16 +8986,12 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
                     s = &(hap->snp_stat.a[hap->list[i].overlapSite]);
                     if(s->occ_0 < 2 || s->occ_1 < 2) continue;
                     if(is_st_bs((*s), st_rate, st_max)) continue;
-                    if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) continue;
+                    if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) continue;
                     if(s->score == 1) continue;
                     o++; 
                     kv_push(uint64_t, hap->snp_srt, hap->list[i].overlapSite);
                 }
-                // if(overlap_list->list[hap->list[l].overlapID].y_id == 317 || overlap_list->list[hap->list[l].overlapID].y_id == 287) {
-                //     fprintf(stderr, "***2***[M::%s-id::%u] o->%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o);
-                // }
                 hap->snp_srt.n -= o;
-                ///there are at least two variants at one read
                 if(o>=(overlap_list->list[hap->list[l].overlapID].align_length*up)) {
                     radix_sort_bc64(hap->snp_srt.a + hap->snp_srt.n, hap->snp_srt.a + hap->snp_srt.n + o);
                     a = hap->snp_srt.a + hap->snp_srt.n;
@@ -9102,9 +9030,6 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
                 for (i = l; i < k; i++) {
                     if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) {
                         s = &(hap->snp_stat.a[hap->list[i].overlapSite]);
-                        if(s->site == 267) {
-                            fprintf(stderr, "[M::%s::]\tsite::%u\tocc0::%u\tocc1::%u\tocc2::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_2);
-                        }
                         if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) {
                             overlap_list->list[ii].strong = 1;
                             if(hh_tp(hap->list[i])==1) {
@@ -9123,6 +9048,7 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio
     }    
 }
 
+
 typedef struct {
     uint64_t *ack;
     uint64_t ack_w64n;///how many words(64) required
@@ -9499,24 +9425,24 @@ void adjust_srt_on_t_obs(srt_on_t *z, uint64_t oid)
 #define g_sn_az(za, zi) (((za).b32->a[(zi)+1]-(za).b32->a[(zi)]))
 #define g_sa_az(za, zi) (((za).b32->a+(za).sn_tot)+(za).b32->a[(zi)])
 
-inline int8_t is_plus_sc_obs(SnpStats *s, double st_rate, uint64_t st_max)
+inline int8_t is_plus_sc_obs(SnpStats *s, double st_rate, uint64_t st_max, uint64_t hap_cov_match, uint64_t hap_cov_unmatch)
 {
     ///if(s->score < 0) return 0;
-    if(s->occ_0 < 2 || s->occ_1 < 2 || s->occ_0 < asm_opt.s_hap_cov || s->occ_1 < asm_opt.infor_cov) return -1;
+    if(s->occ_0 < 2 || s->occ_1 < 2 || s->occ_0 < hap_cov_match || s->occ_1 < hap_cov_unmatch) return -1;
     if(is_st_bs((*s), st_rate, st_max)) return 0;
     return 1;
 }
 
-inline uint8_t is_plus_sc(SnpStats *s, double st_rate, uint64_t st_max)
+inline uint8_t is_plus_sc(SnpStats *s, double st_rate, uint64_t st_max, uint64_t hap_cov_match, uint64_t hap_cov_unmatch)
 {
     ///if(s->score < 0) return 0;
-    if(s->occ_0 < 2 || s->occ_1 < 2 || s->occ_0 < asm_opt.s_hap_cov || s->occ_1 < asm_opt.infor_cov) return 0;
+    if(s->occ_0 < 2 || s->occ_1 < 2 || s->occ_0 < hap_cov_match || s->occ_1 < hap_cov_unmatch) return 0;
     if(is_st_bs((*s), st_rate, st_max)) return 0;
     return 1;
 }
 
 ///sn_tot: how many snps in total
-inline void insert_srt_on_t_i32(srt_on_t *z, haplotype_evdience *a, uint64_t an, asg32_v *bu, SnpStats *sn_a, uint64_t sn_tot, double st_rate, uint64_t st_max)
+inline void insert_srt_on_t_i32(srt_on_t *z, haplotype_evdience *a, uint64_t an, asg32_v *bu, SnpStats *sn_a, uint64_t sn_tot, double st_rate, uint64_t st_max, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     uint64_t k, oid, sid, sa_tot0, ss; 
     uint32_t *snp_a, snp_n, snp_k;
@@ -9534,7 +9460,7 @@ inline void insert_srt_on_t_i32(srt_on_t *z, haplotype_evdience *a, uint64_t an,
         if(g_w_az((*z), oid) == 0) continue;
 
         sid = a[k].overlapSite;
-        if(!is_plus_sc(&sn_a[sid], st_rate, st_max)) continue;
+        if(!is_plus_sc(&sn_a[sid], st_rate, st_max, hap_cov_match, hap_cov_unmatch)) continue;
 
         sn[sid]++; z->sa_tot++;
     }
@@ -9559,7 +9485,7 @@ inline void insert_srt_on_t_i32(srt_on_t *z, haplotype_evdience *a, uint64_t an,
         oid = a[k].overlapID;
         if(g_w_az((*z), oid) == 0) continue;
         sid = a[k].overlapSite; 
-        if(!is_plus_sc(&sn_a[sid], st_rate, st_max)) continue;
+        if(!is_plus_sc(&sn_a[sid], st_rate, st_max, hap_cov_match, hap_cov_unmatch)) continue;
         
         sa_tot0++;
         snp_a = g_sa_az(*z, sid);
@@ -9574,7 +9500,7 @@ inline void insert_srt_on_t_i32(srt_on_t *z, haplotype_evdience *a, uint64_t an,
 }
 
 ///sn_tot: how many snps in total
-inline void insert_srt_on_t_i32_obs(srt_on_t *z, haplotype_evdience *a, uint64_t an, asg32_v *bu, SnpStats *sn_a, uint64_t sn_tot, double st_rate, uint64_t st_max)
+inline void insert_srt_on_t_i32_obs(srt_on_t *z, haplotype_evdience *a, uint64_t an, asg32_v *bu, SnpStats *sn_a, uint64_t sn_tot, double st_rate, uint64_t st_max, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     uint64_t k, oid, sid, sa_tot0, ss; 
     uint32_t *snp_a, snp_n, snp_k;
@@ -9591,7 +9517,7 @@ inline void insert_srt_on_t_i32_obs(srt_on_t *z, haplotype_evdience *a, uint64_t
         oid = a[k].overlapID;
         // if(g_w_az((*z), oid) == 0) continue;
         sid = a[k].overlapSite;
-        zf = is_plus_sc_obs(&sn_a[sid], st_rate, st_max);
+        zf = is_plus_sc_obs(&sn_a[sid], st_rate, st_max, hap_cov_match, hap_cov_unmatch);
         if(zf < 0) continue;
 
         sn[sid]++; z->sa_tot++;
@@ -9617,7 +9543,7 @@ inline void insert_srt_on_t_i32_obs(srt_on_t *z, haplotype_evdience *a, uint64_t
         oid = a[k].overlapID;
         // if(g_w_az((*z), oid) == 0) continue;
         sid = a[k].overlapSite; 
-        zf = is_plus_sc_obs(&sn_a[sid], st_rate, st_max);
+        zf = is_plus_sc_obs(&sn_a[sid], st_rate, st_max, hap_cov_match, hap_cov_unmatch);
         if(zf < 0) continue;
         
         sa_tot0++;
@@ -9763,7 +9689,7 @@ void update_srt_on_t_oid_obs(srt_on_t *z)
 
 
 
-void dbg_srt_ont(srt_on_t *z, haplotype_evdience_alloc* hap, double st_rate, uint64_t st_max, uint64_t oid0, uint64_t w0, uint64_t rid)
+void dbg_srt_ont(srt_on_t *z, haplotype_evdience_alloc* hap, double st_rate, uint64_t st_max, uint64_t oid0, uint64_t w0, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     uint64_t k, l, i, oi, oim = (uint64_t)-1, om = 0, o; SnpStats *s = NULL;
     for (k = 1, l = 0; k <= hap->length; ++k) {  
@@ -9778,7 +9704,7 @@ void dbg_srt_ont(srt_on_t *z, haplotype_evdience_alloc* hap, double st_rate, uin
                 if(hh_tp(hap->list[i])!=1) continue;///mismatch
                 s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
                 assert(s->site == hap->list[i].site);
-                if(!is_plus_sc(s, st_rate, st_max)) continue;
+                if(!is_plus_sc(s, st_rate, st_max, hap_cov_match, hap_cov_unmatch)) continue;
                 o++;///allels must be real
                 // fprintf(stderr, "[M::%s] s->site::%u, s->occ_0::%u, s->occ_1::%u, s->overlap_num::%u, oi::%lu\n", __func__, s->site, s->occ_0, s->occ_1, s->overlap_num, oi);
             }
@@ -9827,7 +9753,7 @@ void dbg_srt_ont_maxk(srt_on_t *z, const char *cmd)
     }
 }
 
-void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid)
+void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid, uint64_t hap_cov_match, uint64_t hap_cov_unmatch)
 {
     if(hap->length == 0) return;
     uint64_t k, l, i, o, ii, m_snp_stat, m_list, m_off;
@@ -9873,7 +9799,7 @@ void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_r
                 if(hh_tp(hap->list[i])!=1) continue;///mismatch
                 s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
                 assert(s->site == hap->list[i].site);
-                if(!is_plus_sc(s, st_rate, st_max)) continue;
+                if(!is_plus_sc(s, st_rate, st_max, hap_cov_match, hap_cov_unmatch)) continue;
                 // if(hap->list[l].overlapID == 35) {
                 //     fprintf(stderr, "[M::%s] site::%u, occ_0::%u, occ_1::%u, occ_2::%u, is_homopolymer::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_1, s->is_homopolymer);
                 // }
@@ -9897,7 +9823,7 @@ void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_r
     if (bs.kan) {
         // fprintf(stderr, "\n-a-[M::%s]\n\n", __func__);
         ///insert SNP -> overlap id
-        insert_srt_on_t_i32(&bs, hap->list, hap->length, b32, hap->snp_stat.a, hap->snp_stat.n, st_rate, st_max);
+        insert_srt_on_t_i32(&bs, hap->list, hap->length, b32, hap->snp_stat.a, hap->snp_stat.n, st_rate, st_max, hap_cov_match, hap_cov_unmatch);
         // fprintf(stderr, "\n-b-[M::%s]\n\n", __func__);
 
         while (pop_srt_k(&bs, &ii, &l, &o)) {///pop one read/overlap
@@ -9938,7 +9864,7 @@ void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_r
                             assert(t->overlap_num >= 1);
                         }
 
-                        if((t->score > 0) && (!is_plus_sc(t, st_rate, st_max))) {
+                        if((t->score > 0) && (!is_plus_sc(t, st_rate, st_max, hap_cov_match, hap_cov_unmatch))) {
                             t->score = -1;
                             ///mark read that need to be updated
                             // fprintf(stderr, "-1-[M::%s]\tg_w_az(bs, 35)::%u\tt->site::%u\n",  __func__, g_w_az(bs, 35), t->site);
@@ -9999,7 +9925,7 @@ void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_r
 
 
 
-void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid)
+void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     if(hap->length == 0) return;
     uint64_t k, l, i, o, obs, ii, m_snp_stat, m_list, m_off;
@@ -10027,7 +9953,7 @@ void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overla
         
             for (; l < k; l++) {
                 s = &(hap->snp_stat.a[l]); assert(s->score > 0);
-                zf = is_plus_sc_obs(s, st_rate, st_max);
+                zf = is_plus_sc_obs(s, st_rate, st_max, hap_cov_match, hap_cov_unmatch);
                 if(zf <= 0) s->score = -s->score;
 
                 hap->snp_stat.a[m_snp_stat++] = hap->snp_stat.a[l];
@@ -10056,7 +9982,7 @@ void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overla
                 if(hh_tp(hap->list[i])!=1) continue;///mismatch
                 s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
                 assert(s->site == hap->list[i].site);
-                zf = is_plus_sc_obs(s, st_rate, st_max);
+                zf = is_plus_sc_obs(s, st_rate, st_max, hap_cov_match, hap_cov_unmatch);
                 // if(s->site == 11114) {
                 //     fprintf(stderr, "[M::%s] s->site::%u, s->occ_0::%u, s->occ_1::%u, s->overlap_num::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->overlap_num);
                 // }
@@ -10088,7 +10014,7 @@ void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overla
     if (bs.kan) {
         // fprintf(stderr, "\n-a-[M::%s]\n\n", __func__);
         ///insert SNP -> overlap id
-        insert_srt_on_t_i32_obs(&bs, hap->list, hap->length, b32, hap->snp_stat.a, hap->snp_stat.n, st_rate, st_max);
+        insert_srt_on_t_i32_obs(&bs, hap->list, hap->length, b32, hap->snp_stat.a, hap->snp_stat.n, st_rate, st_max, hap_cov_match, hap_cov_unmatch);
         // fprintf(stderr, "\n-b-[M::%s]\n\n", __func__);
 
         // fprintf(stderr, "+4+[M::%s] rid::%lu\n", __func__, rid);
@@ -10133,7 +10059,7 @@ void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overla
                             assert(t->overlap_num >= 1);
                         }
 
-                        zf = is_plus_sc(t, st_rate, st_max);
+                        zf = is_plus_sc(t, st_rate, st_max, hap_cov_match, hap_cov_unmatch);
                         // if(s->site == 11114) {
                         //     fprintf(stderr, "[M::%s] s->site::%u, s->occ_0::%u, s->occ_1::%u, s->overlap_num::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->overlap_num);
                         //     fprintf(stderr, "[M::%s]\tzf::%d\tt->score::%d\n",  __func__, zf, t->score);
@@ -10392,7 +10318,10 @@ void generate_haplotypes_naive_HiFi_adv_back(haplotype_evdience_alloc* hap, over
 }
 
 
-void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, uint64_t multi_check, double st_rate, uint64_t st_max, uint64_t snp_dis, int64_t snp_cut)
+
+
+void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, uint64_t multi_check, double st_rate, uint64_t st_max, uint64_t snp_dis, int64_t snp_cut, 
+    int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     // fprintf(stderr, "-0-[M::%s::] Done\n", __func__);
     if(hap->length == 0) return;
@@ -10474,7 +10403,7 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al
                 assert(s->site == hap->list[i].site);
                 if(s->occ_0 < 2 || s->occ_1 < 2) continue;
                 if(is_st_bs((*s), st_rate, st_max)) continue;
-                if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
+                if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) {
                     assert(s->score > 0);
                     o += s->score;// o++;///allels must be real
                     // if(overlap_list->list[hap->list[l].overlapID].y_id == 1740) {
@@ -10484,46 +10413,6 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al
             }
             if(o > omax) o = omax;
 
-            // if(overlap_list->list[hap->list[l].overlapID].y_id == 1740) {
-            //     fprintf(stderr, "str1074[M::%s-id::%u] o->%lu(%c), l::%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o, "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], l);
-            //     for (i = l; i < k; i++) {
-            //         if(hh_tp(hap->list[i])!=1) continue;///mismatch
-            //         s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
-            //         assert(s->site == hap->list[i].site);
-            //         if(s->occ_0 < 2 || s->occ_1 < 2) continue;
-            //         if(is_st_bs((*s), st_rate, st_max)) continue;
-            //         if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
-            //             fprintf(stderr, "[M::%s] site::%u, occ_0::%u, occ_1::%u, occ_2::%u, is_homopolymer::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_1, s->is_homopolymer);
-            //             // prt_sub_read(g_read->seq, g_read->length, s->site, 50);
-            //         }
-            //     }
-            // } 
-            // else {
-                // fprintf(stderr, "[M::%s-id::%u] o->%lu(%c), l::%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o, "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], l);
-            // } 
-            // else {
-            //     for (i = l, o = 0; i < k; i++) {
-            //         if(hh_tp(hap->list[i])!=0) continue;///mismatch
-            //         s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
-            //         assert(s->site == hap->list[i].site);
-            //         if(s->occ_0 < 2 || s->occ_1 < 2) continue;
-            //         if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
-            //             if(s->site == 7878 || s->site == 9682) {
-                        
-            //                 fprintf(stderr, "***[M::%s] site::%u, rid::%u\t%.*s\n", __func__, s->site, overlap_list->list[hap->list[l].overlapID].y_id, (int)Get_NAME_LENGTH(R_INF, overlap_list->list[hap->list[l].overlapID].y_id), 
-            //                     Get_NAME(R_INF, overlap_list->list[hap->list[l].overlapID].y_id));
-            //                 prt_sub_cigar(&(overlap_list->list[hap->list[l].overlapID]), g_read->length, s->site, 50);
-            //             }
-            //             // o++;///allels must be real
-            //             // fprintf(stderr, "[M::%s] site::%u, occ_0::%u, occ_1::%u, occ_2::%u, is_homopolymer::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_1, s->is_homopolymer);
-            //             // prt_sub_read(g_read->seq, g_read->length, s->site, 50);
-            //         }
-            //     }
-            // }
-
-            // if(overlap_list->list[hap->list[l].overlapID].y_id == 3626) {
-            //     fprintf(stderr, "***0***[M::%s-id::%u] o->%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o);
-            // }
 
             if(o > 0) {
                 o = omax - o;
@@ -10548,7 +10437,7 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al
                 s = &(hap->snp_stat.a[hap->list[i].overlapSite]);
                 if(s->occ_0 < 2 || s->occ_1 < 2) continue;
                 if(is_st_bs((*s), st_rate, st_max)) continue;
-                if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
+                if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) {
                     o++;
                     // if(overlap_list->list[hap->list[l].overlapID].y_id == 1740) {
                         // fprintf(stderr, "srt[M::%s] site->%u, occ_0::%u, occ_1::%u\n", __func__, s->site, s->occ_0, s->occ_1);
@@ -10567,7 +10456,7 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al
                     s = &(hap->snp_stat.a[hap->list[i].overlapSite]);
                     if(s->occ_0 < 2 || s->occ_1 < 2) continue;
                     if(is_st_bs((*s), st_rate, st_max)) continue;
-                    if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
+                    if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) {
                         s->id = UINT32_MAX; ///s->score = 1; 
                         // fprintf(stderr, "set[M::%s-site::%u] occ_0::%u, occ_1::%u\n", __func__, s->site, s->occ_0, s->occ_1);
                     }
@@ -10636,7 +10525,7 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al
                     s = &(hap->snp_stat.a[hap->list[i].overlapSite]);
                     if(s->occ_0 < 2 || s->occ_1 < 2) continue;
                     if(is_st_bs((*s), st_rate, st_max)) continue;
-                    if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) continue;
+                    if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) continue;
                     if(s->id == UINT32_MAX) continue; //if(s->score == 1) continue;
                     o++; 
                     kv_push(uint64_t, hap->snp_srt, hap->list[i].overlapSite);
@@ -10713,9 +10602,7 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al
 }
 
 
-
-
-void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, uint64_t rid)
+void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     uint64_t k, l, i, o, ii;
     int64_t z; SnpStats *s = NULL;
@@ -10731,7 +10618,7 @@ void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc*
                 s = &(hap->snp_stat.a[hap->list[i].overlapSite]); 
                 assert(s->site == hap->list[i].site);
                 if(s->occ_0 < 2 || s->occ_1 < 2) continue;
-                if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++;///allels must be real
+                if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) o++;///allels must be real
             }
 
             if(o > 0) {
@@ -10752,7 +10639,7 @@ void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc*
                 if(hh_tp(hap->list[i])!=1) continue;
                 s = &(hap->snp_stat.a[hap->list[i].overlapSite]);
                 if(s->occ_0 < 2 || s->occ_1 < 2) continue;
-                if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) {
+                if(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch) {
                     o++;
                 }
             }
@@ -10844,6 +10731,7 @@ void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc*
 }
 
 
+
 inline int64_t comput_sc_rphase(SnpStats *ai, uint64_t id, SnpStats *aj, uint64_t jd, haplotype_evdience *za, uint64_t occ0_cut)
 {
     if(ai->site == aj->site) return INT64_MIN;
@@ -11032,7 +10920,7 @@ void gen_rphase_dp0_multiple_path(SnpStats *a, int64_t an, haplotype_evdience *z
     
 }
 
-int64_t is_hpc_vec(SnpStats *ai, uint64_t id, haplotype_evdience *za)
+int64_t is_hpc_vec(SnpStats *ai, uint64_t id, haplotype_evdience *za, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     haplotype_evdience *iz = NULL; int64_t in, ik, n0 = ai->occ_0, n1 = ai->occ_1, f = 0; 
     iz = za + ai->non_homopolymer_num; in = ai->homopolymer_num - ai->non_homopolymer_num;
@@ -11045,7 +10933,7 @@ int64_t is_hpc_vec(SnpStats *ai, uint64_t id, haplotype_evdience *za)
             ai->occ_1 -= iz[ik].cov;
         }
     }
-    if((ai->occ_0 < 2 || ai->occ_1 < 2) || (!(ai->occ_0 >= asm_opt.s_hap_cov && ai->occ_1 >= asm_opt.infor_cov))) f = 1;
+    if((ai->occ_0 < 2 || ai->occ_1 < 2) || (!(ai->occ_0 >= hap_cov_match && ai->occ_1 >= hap_cov_unmatch))) f = 1;
     ai->occ_0 = n0; ai->occ_1 = n1;
     // if((sec_check) && (((n0)<=((n0+n1)*0.333333)) || ((n1)<=((n0+n1)*0.333333)))) f = 1;
     return f;
@@ -11095,7 +10983,7 @@ int32_t get_hq_value(SnpStats *ai, uint64_t id, haplotype_evdience *za, int64_t
     }
 }
 
-void gen_rphase_dp0_single_path(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, uint8_t site_sc)
+void gen_rphase_dp0_single_path(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, uint8_t site_sc, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     if(an <= 0) return;
     int64_t *p, i, k, j, st, max_f, max_j, sc, ch_n, plus = 0, rn, rn0; int32_t *f, *ii; uint64_t m, cc0, cc = 0, cci, cc_min; int64_t b0l, b0h, b1l, b1h, krn;
@@ -11158,7 +11046,7 @@ void gen_rphase_dp0_single_path(SnpStats *a, int64_t an, haplotype_evdience *za,
             if(rn > 1) {
                 plus = 1;
             } else {
-                if((!is_hpc_vec(&(a[res->a[rn0]]), res->a[rn0], za)) && (a[res->a[rn0]].occ_0 >= cc)) plus = 1;
+                if((!is_hpc_vec(&(a[res->a[rn0]]), res->a[rn0], za, hap_cov_match, hap_cov_unmatch)) && (a[res->a[rn0]].occ_0 >= cc)) plus = 1;
             }
 
             for (i = 0; i < rn; i++) {
@@ -11271,7 +11159,7 @@ inline uint64_t get_occ0(SnpStats *ai, uint64_t id, haplotype_evdience *za, over
 }
 
 void gen_rphase_dp0_single_path_hybrid_0(SnpStats *a, int64_t an, haplotype_evdience *za, int32_t *f, int64_t *p, int32_t *ii, uint8_t *qual_a, uint64_t *ra, asg64_v *idx, int64_t plus, uint64_t cc, uint64_t cci, uint64_t cc_min,
-    double h_rate, overlap_region_alloc *ol, uint64_t rid, uint64_t tcut, uint64_t hf_only, int64_t n_hap, double cut_rate, uint64_t cut_bd, uint64_t site_sc)
+    double h_rate, overlap_region_alloc *ol, uint64_t rid, uint64_t tcut, uint64_t hf_only, int64_t n_hap, double cut_rate, uint64_t cut_bd, uint64_t site_sc, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     int64_t i, j, k, rn, rn0, ch_n; uint64_t m, cc0 = cc; idx->n = 0; int64_t b0l, b0h, b1l, b1h, krn;
     for (i = 0; i < an; i++) {
@@ -11310,7 +11198,7 @@ void gen_rphase_dp0_single_path_hybrid_0(SnpStats *a, int64_t an, haplotype_evdi
                 plus = 1;
             } else {
                 // if((!is_hpc_vec(&(a[ra[rn0]]), ra[rn0], za)) && (a[ra[rn0]].occ_0 >= cc)) plus = 1;
-                if((!is_hpc_vec(&(a[ra[rn0]]), ra[rn0], za)) && (get_occ0((&a[ra[rn0]]), ra[rn0], za, ol, rid, tcut) >= cc)) plus = 1;
+                if((!is_hpc_vec(&(a[ra[rn0]]), ra[rn0], za, hap_cov_match, hap_cov_unmatch)) && (get_occ0((&a[ra[rn0]]), ra[rn0], za, ol, rid, tcut) >= cc)) plus = 1;
             }
 
             if(hf_only && plus == -1) continue;
@@ -11389,8 +11277,9 @@ void gen_rphase_dp0_single_path_hybrid_0(SnpStats *a, int64_t an, haplotype_evdi
 }
 
 
+
 void gen_rphase_dp0_single_path_hybrid_0_multi(SnpStats *a, int64_t an, haplotype_evdience *za, int32_t *f, int64_t *p, int32_t *ii, uint8_t *qual_a, uint64_t *rz, asg64_v *idx, int64_t plus, uint64_t cc, uint64_t cci, uint64_t cc_min,
-    overlap_region_alloc *ol, uint64_t rid, uint64_t tcut, uint64_t hf_only, int64_t n_hap, double cut_rate, uint64_t cut_bd, uint64_t site_sc)
+    overlap_region_alloc *ol, uint64_t rid, uint64_t tcut, uint64_t hf_only, int64_t n_hap, double cut_rate, uint64_t cut_bd, uint64_t site_sc, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     int64_t ri, rj, rk, rn; uint64_t m, cc0 = cc, occ0; idx->n = 0; int64_t b0l, b0h, b1l, b1h, krn;
     for (ri = 0; ri < an; ri++) {
@@ -11430,14 +11319,14 @@ void gen_rphase_dp0_single_path_hybrid_0_multi(SnpStats *a, int64_t an, haplotyp
                 if(cc > cc0) cc = cc0;
 
                 if(krn > 1) {
-                    if(b0h > b0l && b0h >= asm_opt.s_hap_cov && b1h > b1l && b1h >= asm_opt.infor_cov) {
+                    if(b0h > b0l && b0h >= hap_cov_match && b1h > b1l && b1h >= hap_cov_unmatch) {
                         a[rz[rk]].score = plus;
                     } else {
                         a[rz[rk]].score = -1;
                     }
                 } else {
-                    if(((b0h > b0l) && (b0h > ((b0h + b0l)*0.7)) && (b0h >= asm_opt.s_hap_cov) && (b0h >= (int64_t)cc)) && 
-                            ((b1h > b1l) && (b1h > ((b1h + b1l)*0.7)) && (b1h >= asm_opt.infor_cov) && (b1h >= (int64_t)cc))) {
+                    if(((b0h > b0l) && (b0h > ((b0h + b0l)*0.7)) && (b0h >= hap_cov_match) && (b0h >= (int64_t)cc)) && 
+                            ((b1h > b1l) && (b1h > ((b1h + b1l)*0.7)) && (b1h >= hap_cov_unmatch) && (b1h >= (int64_t)cc))) {
                         a[rz[rk]].score = plus;
                     } else {
                         a[rz[rk]].score = -1;
@@ -11484,7 +11373,7 @@ void gen_rphase_dp0_single_path_hybrid_0_multi(SnpStats *a, int64_t an, haplotyp
                 if(rn > 1) {
                     plus = krn/**1**/;
                 } else {
-                    if(!is_hpc_vec(&(a[rz[rk]]), rz[rk], za)) {
+                    if(!is_hpc_vec(&(a[rz[rk]]), rz[rk], za, hap_cov_match, hap_cov_unmatch)) {
                         occ0 = get_occ0((&a[rz[rk]]), rz[rk], za, ol, rid, tcut);
                         if(occ0 >= cc) plus = krn/**1**/;
                     }
@@ -11503,7 +11392,8 @@ void gen_rphase_dp0_single_path_hybrid_0_multi(SnpStats *a, int64_t an, haplotyp
     }
 }
 
-void gen_rphase_dp0_single_path_multi(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, uint8_t site_sc)
+
+void gen_rphase_dp0_single_path_multi(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, uint8_t site_sc, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     if(an <= 0) return;
     int64_t *p, ri, rj, st, max_f, max_j, sc, plus = 0; int32_t *f, *ii; uint64_t cc = 0, cci, cc_min;
@@ -11533,10 +11423,11 @@ void gen_rphase_dp0_single_path_multi(SnpStats *a, int64_t an, haplotype_evdienc
     }
 
     kv_resize(uint64_t, *res, ((uint64_t)an));
-    gen_rphase_dp0_single_path_hybrid_0_multi(a, an, za, f, p, ii, qual_a, res->a, idx, plus, cc, cci, cc_min, NULL, ((uint64_t)-1), ((uint64_t)-1), 0, n_hap, cut_rate, cut_bd, site_sc);
+    gen_rphase_dp0_single_path_hybrid_0_multi(a, an, za, f, p, ii, qual_a, res->a, idx, plus, cc, cci, cc_min, NULL, ((uint64_t)-1), ((uint64_t)-1), 0, n_hap, cut_rate, cut_bd, site_sc, hap_cov_match, hap_cov_unmatch);
 }
 
-void gen_rphase_dp0_single_path_hybrid(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, double cut_rate_hf, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, overlap_region_alloc *ol, uint64_t rid, uint64_t tcut, uint64_t site_sc)
+void gen_rphase_dp0_single_path_hybrid(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, double cut_rate_hf, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, overlap_region_alloc *ol, uint64_t rid, uint64_t tcut, uint64_t site_sc, 
+    int64_t hap_cov_match, int64_t hap_cov_unmatch, double hf_rate)
 {
     if(an <= 0) return;
     int64_t *p_h, *p_a, i, j, st, max_f_h, max_f_a, max_j_h, max_j_a, sc_h, sc_a, plus_h = 0, plus_a = 0; int32_t *f_h, *f_a, *ii_h, *ii_a; uint64_t cc = 0, cc0 = 0, cci, cc_min;
@@ -11583,13 +11474,13 @@ void gen_rphase_dp0_single_path_hybrid(SnpStats *a, int64_t an, haplotype_evdien
    
 
     kv_resize(uint64_t, *res, ((uint64_t)an));
-    gen_rphase_dp0_single_path_hybrid_0_multi(a, an, za, f_a, p_a, ii_a, qual_a, res->a, idx, plus_a, cc, cci, cc_min, ol, rid, (uint64_t)-1, 0, n_hap, cut_rate, cut_bd, site_sc);
+    gen_rphase_dp0_single_path_hybrid_0_multi(a, an, za, f_a, p_a, ii_a, qual_a, res->a, idx, plus_a, cc, cci, cc_min, ol, rid, (uint64_t)-1, 0, n_hap, cut_rate, cut_bd, site_sc, hap_cov_match, hap_cov_unmatch);
 
-    cc0 = cc; cc = ((het_cov > 0)?(het_cov):(hom_cov/n_hap)); cc *= (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1]))); cc *= cut_rate_hf; 
-    if(cc < ((uint64_t)(MAX(asm_opt.s_hap_cov, asm_opt.infor_cov) + 1))) cc = MAX(asm_opt.s_hap_cov, asm_opt.infor_cov) + 1; 
+    cc0 = cc; cc = ((het_cov > 0)?(het_cov):(hom_cov/n_hap)); cc *= hf_rate; cc *= cut_rate_hf; 
+    if(cc < ((uint64_t)(MAX(hap_cov_match, hap_cov_unmatch) + 1))) cc = MAX(hap_cov_match, hap_cov_unmatch) + 1; 
     if(cc <= 0) cc = 1;
     if(cc < cc0) cc = cc0;
-    /**if(rid >= tcut)**/ gen_rphase_dp0_single_path_hybrid_0_multi(a, an, za, f_h, p_h, ii_h, qual_a, res->a, idx, plus_h, cc, cci, cc_min, ol, rid, tcut, 1, n_hap, cut_rate, cut_bd, site_sc);///additional SNPs
+    /**if(rid >= tcut)**/ gen_rphase_dp0_single_path_hybrid_0_multi(a, an, za, f_h, p_h, ii_h, qual_a, res->a, idx, plus_h, cc, cci, cc_min, ol, rid, tcut, 1, n_hap, cut_rate, cut_bd, site_sc, hap_cov_match, hap_cov_unmatch);///additional SNPs
 }
 
 inline void fill_incom(asg64_v *om, uint64_t oid, uint64_t pe, uint64_t *idx_a, int64_t idx_n, uint64_t qid)
@@ -11643,7 +11534,8 @@ void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *v
     // fprintf(stderr, "[M::%s] s::%lu, e::%lu, wk::%lu, scw::%lu, msc/wk::%lu, wqual::%lu\n", __func__, s, e, wk, scw, msc/wk, (*wqual));
 }
 
-void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, double st_rate, uint64_t st_max, asg64_v *idx, asg64_v *res, uint64_t rid, uint8_t *qa, uint64_t tcut)
+
+void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, double st_rate, uint64_t st_max, asg64_v *idx, asg64_v *res, uint64_t rid, uint8_t *qa, uint64_t tcut, int64_t hap_cov_match, int64_t hap_cov_unmatch)
 {
     if(hl->length <= 0) return;
 
@@ -11662,7 +11554,7 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub
             for (o = l, m1 = 0; o < k; o++) {
                 s = &(hl->snp_stat.a[o]);
                 // fprintf(stderr, "+[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1);
-                if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov))) {
+                if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch))) {
                     continue;
                 }
                 // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1);
@@ -11725,7 +11617,7 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub
                 if((hq[4] >= hq_cut) && (hq[0] >= hq_cut || hq[1] >= hq_cut || hq[2] >= hq_cut || hq[3] >= hq_cut)) {
                     for (o = l, m_snp_stat0 = m_snp_stat; o < k; o++) {
                         s = &(hl->snp_stat.a[o]);
-                        if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov))) {
+                        if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch))) {
                             continue;
                         }
                         
@@ -11845,6 +11737,7 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub
     }
 }
 
+
 /**
  * //r829
 void recal_rphase0(overlap_region *z, SnpStats *sa, int64_t sn, haplotype_evdience_alloc *hl, int64_t hn, uint64_t *idx, int64_t idx_n)
@@ -11916,26 +11809,23 @@ void recal_rphase(All_reads *rref, haplotype_evdience_alloc *hl, overlap_region_
 }
 **/
 
-void gen_rphase_dp_adv(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Read* g_read, double st_rate, uint64_t st_max, Chain_Data *dp, asg64_v *idx, asg64_v *res, uint64_t rid, uint8_t *qa, uint64_t tcut, uint64_t site_sc)
+
+void gen_rphase_dp_adv(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Read* g_read, double st_rate, uint64_t st_max, Chain_Data *dp, asg64_v *idx, asg64_v *res, uint64_t rid, uint8_t *qa, uint64_t tcut, uint64_t site_sc,
+    int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t n_hap, int64_t het_c, int64_t hom_c, double hf_rate)
 {
     if(hl->length <= 0) return;
     uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list /**m_off**/, m1, c0, c1, rev_n; SnpStats *s; 
-    int64_t het_a, hom_a;
 
-    call_rphase_sc(hl, ol, st_rate, st_max, idx, res, rid, qa, tcut);
+    call_rphase_sc(hl, ol, st_rate, st_max, idx, res, rid, qa, tcut, hap_cov_match, hap_cov_unmatch);
 
-    het_a = asm_opt.het_cov; hom_a = asm_opt.hom_cov;
-    if(asm_opt.het_cov_set >= 0) het_a = asm_opt.het_cov_set;
-    if(het_a < 0) {
-        het_a = hom_a/asm_opt.polyploidy;
-    } else if(het_a > (hom_a/asm_opt.polyploidy)) {
-        het_a = hom_a/asm_opt.polyploidy;
-    }
+    // fprintf(stderr, "+[M::%s]\tyid::%u\thap_cov_match::%ld\thap_cov_unmatch::%ld\tn_hap::%ld\thet_c::%ld\thom_c::%ld\thf_rate::%f\n",
+    //          __func__, ol->list[0].y_id, hap_cov_match, hap_cov_unmatch, n_hap, het_c, hom_c, hf_rate);
+    
     if(tcut == ((uint64_t)-1)) {
-        // gen_rphase_dp0_single_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, asm_opt.polyploidy, 0.7, 6, res, qv->a, site_sc);
-        gen_rphase_dp0_single_path_multi(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, asm_opt.polyploidy, 0.6/**0.7**/, 6, res, qa, site_sc);
+        // gen_rphase_dp0_single_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, n_hap, 0.7, 6, res, qv->a, site_sc);
+        gen_rphase_dp0_single_path_multi(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.6/**0.7**/, 6, res, qa, site_sc, hap_cov_match, hap_cov_unmatch);
     } else {
-        gen_rphase_dp0_single_path_hybrid(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, asm_opt.polyploidy, 0.6/**0.7**/, 0.4, 6, res, qa, ol, rid, tcut, site_sc);
+        gen_rphase_dp0_single_path_hybrid(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.6/**0.7**/, 0.4, 6, res, qa, ol, rid, tcut, site_sc, hap_cov_match, hap_cov_unmatch, hf_rate);
     }
 
     for (k = 1, l = 0, i = m_snp_stat = m_list = 0; k <= hl->snp_stat.n; ++k) {///filter snps
@@ -12003,12 +11893,13 @@ void gen_rphase_dp_adv(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, U
     // fprintf(stderr, "[M::%s] m_snp_stat->%lu\n", __func__, m_snp_stat);
 }
 
-void gen_rphase_dp(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Read* g_read, double st_rate, uint64_t st_max, Chain_Data *dp, asg64_v *idx, asg64_v *res, uint64_t rid, asg8_v *qv, uint64_t tcut, uint64_t site_sc)
+
+void gen_rphase_dp(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Read* g_read, double st_rate, uint64_t st_max, Chain_Data *dp, asg64_v *idx, asg64_v *res, uint64_t rid, asg8_v *qv, uint64_t tcut, uint64_t site_sc,
+    int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t n_hap, int64_t het_c, int64_t hom_c, double hf_rate)
 {
     if(hl->length <= 0) return;
 
     uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list /**m_off**/, m1, c0, c1, rev_n, tqual, wqual, hq_cut = 2; uint32_t hq[5], hp[4], is_st; SnpStats *s; haplotype_evdience ev; char mc;
-    int64_t het_a, hom_a;
     if(rid < tcut) {
         retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn);
     } else {
@@ -12026,7 +11917,7 @@ void gen_rphase_dp(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Re
             for (o = l, m1 = 0; o < k; o++) {
                 s = &(hl->snp_stat.a[o]);
                 // fprintf(stderr, "+[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1);
-                if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov))) {
+                if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch))) {
                     continue;
                 }
                 // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1);
@@ -12089,7 +11980,7 @@ void gen_rphase_dp(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Re
                 if((hq[4] >= hq_cut) && (hq[0] >= hq_cut || hq[1] >= hq_cut || hq[2] >= hq_cut || hq[3] >= hq_cut)) {
                     for (o = l, m_snp_stat0 = m_snp_stat; o < k; o++) {
                         s = &(hl->snp_stat.a[o]);
-                        if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov))) {
+                        if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch))) {
                             continue;
                         }
                         
@@ -12209,18 +12100,12 @@ void gen_rphase_dp(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Re
     }
     // fprintf(stderr, "-[M::%s]\tres->n::%lu\tol->length::%lu\tsnp_stat.n::%lu\n", __func__, (uint64_t)res->n, ol->length, (uint64_t)hl->snp_stat.n);
     // gen_rphase_dp0_multiple_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, res);
-    het_a = asm_opt.het_cov; hom_a = asm_opt.hom_cov;
-    if(asm_opt.het_cov_set >= 0) het_a = asm_opt.het_cov_set;
-    if(het_a < 0) {
-        het_a = hom_a/asm_opt.polyploidy;
-    } else if(het_a > (hom_a/asm_opt.polyploidy)) {
-        het_a = hom_a/asm_opt.polyploidy;
-    }
+   
     if(tcut == ((uint64_t)-1)) {
-        // gen_rphase_dp0_single_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, asm_opt.polyploidy, 0.7, 6, res, qv->a, site_sc);
-        gen_rphase_dp0_single_path_multi(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, asm_opt.polyploidy, 0.7, 6, res, qv->a, site_sc);
+        // gen_rphase_dp0_single_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, n_hap, 0.7, 6, res, qv->a, site_sc);
+        gen_rphase_dp0_single_path_multi(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.7, 6, res, qv->a, site_sc, hap_cov_match, hap_cov_unmatch);
     } else {
-        gen_rphase_dp0_single_path_hybrid(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, asm_opt.polyploidy, 0.7, 0.4, 6, res, qv->a, ol, rid, tcut, site_sc);
+        gen_rphase_dp0_single_path_hybrid(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.7, 0.4, 6, res, qv->a, ol, rid, tcut, site_sc, hap_cov_match, hap_cov_unmatch, hf_rate);
     }
 
 
@@ -13071,7 +12956,7 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R
     hap->length = m;
 
     // generate_haplotypes_naive_advance(hap, overlap_list, NULL);
-    generate_haplotypes_naive_HiFi(hap, overlap_list, 0.04, g_read, 1, 0, ((uint64_t)-1));
+    generate_haplotypes_naive_HiFi(hap, overlap_list, 0.04, g_read, 1, 0, ((uint64_t)-1), asm_opt.s_hap_cov, asm_opt.infor_cov);
     // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat);
     // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat);
 
@@ -21620,6 +21505,62 @@ int64_t extract_sub_cigar_hc(overlap_region *z, All_reads *rref, haplotype_evdie
     return 1;
 }
 
+
+int64_t extract_sub_err_hc(overlap_region *z, uint64_t ql, int64_t s, int64_t e, ul_ov_t *p, int64_t *rerr)
+{
+    int64_t wk = ovlp_cur_wid(*p), xk = ovlp_cur_xoff(*p), yk = ovlp_cur_yoff(*p), ck = ovlp_cur_coff(*p), os, oe;
+    bit_extz_t ez; int64_t s0, e0; (*rerr) = INT64_MAX;
+    s0 = ((int64_t)(z->w_list.a[wk].x_start)); 
+    e0 = ((int64_t)(z->w_list.a[wk].x_end)) + 1;
+    if(s < s0) {s = s0;} if(e > e0) {e = e0;}///exclude boundary
+    if(s >= e) return -1;
+    os = MAX(s, s0); oe = MIN(e, e0);
+    if(oe <= os) return -1;
+
+    set_bit_extz_t(ez, (*z), wk);
+    if(!ez.cigar.n) return -1;
+    int64_t cn = ez.cigar.n, op; int64_t ws, we, ovlp; 
+    if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed
+        ck = 0; xk = ez.ts; yk = ez.ps;
+    }
+    
+    while (ck > 0 && xk > s) {///x -> t; y -> p
+        --ck;
+        op = ez.cigar.a[ck]>>14;
+        if(op!=2) xk -= (ez.cigar.a[ck]&(0x3fff));
+        if(op!=3) yk -= (ez.cigar.a[ck]&(0x3fff)); 
+    }
+    
+    //some cigar will span s or e
+    (*rerr) = 0;
+    while (ck < cn && xk < e) {//[s, e)
+        ws = xk; 
+        op = ez.cigar.a[ck]>>14;
+        if(op!=2) xk += (ez.cigar.a[ck]&(0x3fff)); ///op == 3
+        if(op!=3) yk += (ez.cigar.a[ck]&(0x3fff)); ///op == 2
+        ck++; we = xk;
+
+        os = MAX(s, ws); oe = MIN(e, we);
+        ovlp = ((oe>os)? (oe-os):0);
+
+        if(op != 2) {
+            if(!ovlp) continue;
+        } else {///ws == we
+            if(ws < s || ws >= e) continue;
+        }
+
+        if(op == 1 || op == 3) {
+            (*rerr) += ovlp;
+        } else if(op == 2) {
+            (*rerr) += (ez.cigar.a[ck-1]&(0x3fff));
+        }
+    }
+
+    ovlp_cur_xoff(*p) = xk; ovlp_cur_yoff(*p) = yk; ovlp_cur_coff(*p) = ck; ovlp_cur_ylen(*p) = 0;
+    return 1;
+}
+
+
 inline char* cal_hpc_len(All_reads *rref, uint8_t rev, uint32_t id, int64_t s0, int64_t e0, int64_t cs, int64_t ce, int64_t *rs, int64_t *re, UC_Read *z)
 {
     int64_t os, oe, ol, cl; (*rs) = (*re) = -1;
@@ -22241,6 +22182,176 @@ uint64_t hc_phase_robust_rr(overlap_region* ol, All_reads *rref, haplotype_evdie
     return rr;
 }
 
+inline uint64_t median_inplace(uint64_t *a, uint64_t m, uint8_t is_srt) 
+{
+    if(is_srt) radix_sort_bc64(a, a + m);
+    if (m == 0) return 0;
+    if (m&1) return a[m>>1];
+    return (a[(m>>1)-1]+a[m>>1])>>1;
+}
+
+inline uint64_t mad_from_sorted(uint64_t *a, uint64_t m, uint64_t med) {
+    // reuse array: convert to abs deviations, median them
+    uint64_t k;
+    for (k = 0; k < m; ++k) {
+        a[k] = ((a[k]>=med)?(a[k]-med):(med-a[k])); 
+    }
+    return median_inplace(a, m, 1);
+}
+
+uint64_t hc_est_robust_rr_diff_gap(overlap_region* ol, uint64_t ql, uint64_t *id_a, uint64_t id_n, uint64_t *srt_a, uint64_t *srt_b, uint64_t s, uint64_t e, uint64_t win_g, double k_mad, 
+    ul_ov_t *c_idx, int64_t *est_err, uint8_t is_dbg)
+{
+    uint64_t k, q[2], rr = 0, os, oe, srt_n = 0; int64_t rerr, es[2], eg; ul_ov_t *p; overlap_region *z;
+    *est_err = INT64_MAX;
+    if(is_dbg) fprintf(stderr, "\n[M::%s]\tpos::[%lu,\t%lu)\n", __func__, s, e);
+    for (k = 0; k < id_n; k++) {
+        p = &(c_idx[id_a[k]]); z = &(ol[ovlp_id(*p)]); 
+        q[0] = z->w_list.a[ovlp_cur_wid(*p)].x_start;
+        q[1] = z->w_list.a[ovlp_cur_wid(*p)].x_end+1;
+        if(q[1] <= e) rr = 1;
+        os = MAX(q[0], s); oe = MIN(q[1], e);
+        if((oe > os) && (oe - os == e - s)) {
+            extract_sub_err_hc(z, ql, os, oe, p, &rerr);
+            // if(is_dbg) fprintf(stderr, "[M::%s]\ttn::%u\t%c\tq::[%lu,\t%lu)\to::[%lu,\t%lu)\trerr::%ld\n", 
+            //     __func__, z->y_id, "+-"[z->y_pos_strand], q[0], q[1], os, oe, rerr);
+            srt_a[srt_n] = rerr; srt_a[srt_n] <<= 32; srt_a[srt_n] |= k; srt_n++;
+            // extract_sub_cigar_hc_hpc(z, rref, hp, qstr, ql, tu, os, oe, p, set_f, hp->flag + os - s, occ_thres, hpc_len, h0_w);
+        }
+    }
+
+    if(srt_n > 0) {
+        double eps = 1.0, sc, med, mad, sigma, thr; uint64_t bk = srt_n, wgn, med_f, mad_f;
+        radix_sort_bc64(srt_a, srt_a + srt_n);
+        if(is_dbg) {
+            for (k = 0; k < srt_n; k++) {
+                p = &(c_idx[id_a[(uint32_t)srt_a[k]]]); z = &(ol[ovlp_id(*p)]); 
+                fprintf(stderr, "[M::%s]\ttn::%u\t%c\tq::[%d,\t%d)\tlocal_err::%lu\tglobal_err::%u\n", 
+                    __func__, z->y_id, "+-"[z->y_pos_strand], z->w_list.a[ovlp_cur_wid(*p)].x_start, z->w_list.a[ovlp_cur_wid(*p)].x_end+1, srt_a[k]>>32, z->non_homopolymer_errors);
+            }
+        }
+        
+        if((srt_n > 1) && ((srt_a[0]>>32) != (srt_a[srt_n-1]>>32))) {
+            wgn = MIN((srt_n-1), win_g);
+
+            for (k = 1; k < srt_n; k++) {
+                es[0] = srt_a[k-1]>>32;
+                es[1] = srt_a[k]>>32;
+                eg = es[1] - es[0];
+
+                sc = ((double)eg)/(((double)es[0]) + eps);
+                rr = sc * 1000; if(eg > 0 && rr == 0) rr = 1;
+                srt_b[k-1] = rr; 
+
+                if(is_dbg) {
+                    fprintf(stderr, "[M::%s]\tes[0]::%ld,\tes[1]::%ld,\teg::%ld,\tsc::%f\trr::%lu\n", 
+                        __func__, es[0], es[1], eg, sc, rr);
+                }
+            }
+            med_f = median_inplace(srt_b, wgn, 1);
+            mad_f = mad_from_sorted(srt_b, wgn, med_f);
+            med = ((double)med_f)/1000.0;
+            mad = ((double)mad_f)/1000.0;
+            sigma = mad * 1.4826;
+            thr = (sigma > 1e-12) ? (med + k_mad * sigma) : (med + 0.5);
+            if(is_dbg) {
+                fprintf(stderr, "[M::%s]\tmed_f::%lu,\tmed::%f,\tmad_f::%lu,\tmad::%f,\tsigma::%f,\tthr::%f\n", 
+                    __func__, med_f, med, mad_f, mad, sigma, thr);
+            }
+
+            for (k = 1; k < srt_n; k++) {
+                es[0] = srt_a[k-1]>>32;
+                es[1] = srt_a[k]>>32;
+                eg = es[1] - es[0];
+                if(eg <= 0) continue;
+                sc = ((double)eg)/(((double)es[0]) + eps);
+                if(sc >= thr) {
+                    bk = k; break;
+                }
+            }                
+        } 
+        *est_err = (srt_a[bk-1]>>32); ///e_cut_n = bk;
+    }
+    if(is_dbg) fprintf(stderr, "[M::%s]\t*est_err::%ld\n", __func__, *est_err);
+
+    return rr;
+}
+
+
+uint64_t hc_est_robust_rr(overlap_region* ol, uint64_t ql, uint64_t *id_a, uint64_t id_n, uint64_t *srt_a, uint64_t *srt_b, uint64_t s, uint64_t e, double k_mad, ul_ov_t *c_idx, 
+    int64_t *max_est_err, int64_t *ave_est_err, uint8_t is_dbg)
+{
+    uint64_t k, q[2], rr = 0, os, oe, srt_n = 0, cut; int64_t rerr, mm_cut = 0; ul_ov_t *p; overlap_region *z;
+    *max_est_err = *ave_est_err = INT64_MAX;
+    // if(is_dbg) fprintf(stderr, "\n[M::%s]\tpos::[%lu,\t%lu)\n", __func__, s, e);
+    for (k = 0; k < id_n; k++) {
+        // fprintf(stderr, "[M::%s]\tid_a[%lu]::%lu\n", __func__, k, id_a[k]);
+        p = &(c_idx[id_a[k]]); z = &(ol[ovlp_id(*p)]); 
+        q[0] = z->w_list.a[ovlp_cur_wid(*p)].x_start;
+        q[1] = z->w_list.a[ovlp_cur_wid(*p)].x_end+1;
+        if(q[1] <= e) rr = 1;
+        os = MAX(q[0], s); oe = MIN(q[1], e);
+        if((oe > os) && (oe - os == e - s)) {
+            extract_sub_err_hc(z, ql, os, oe, p, &rerr);
+            // if(is_dbg) fprintf(stderr, "[M::%s]\ttn::%u\t%c\tq::[%lu,\t%lu)\to::[%lu,\t%lu)\trerr::%ld\n", 
+            //     __func__, z->y_id, "+-"[z->y_pos_strand], q[0], q[1], os, oe, rerr);
+            srt_a[srt_n] = srt_b[srt_n] = rerr; srt_b[srt_n] <<= 32; srt_b[srt_n] |= k;
+            srt_n++;
+            if(rerr > mm_cut) mm_cut = rerr;
+            // if(is_dbg) fprintf(stderr, "[M::%s]\ttn::%u\t%c\tq::[%d,\t%d)\tlocal_err::%ld\tglobal_err::%u\n", 
+            //     __func__, z->y_id, "+-"[z->y_pos_strand], z->w_list.a[ovlp_cur_wid(*p)].x_start, z->w_list.a[ovlp_cur_wid(*p)].x_end+1, rerr, z->non_homopolymer_errors);
+        }
+    }
+
+    if(srt_n > 0) {
+        double mad; uint64_t med_f, mad_f, thr;
+        radix_sort_bc64(srt_a, srt_a + srt_n);
+        // if(is_dbg) {
+        //     for (k = 0; k < srt_n; k++) {
+        //         fprintf(stderr, "[M::%s]\tsorted_local_err::%lu\n", 
+        //             __func__, srt_a[k]);
+        //     }
+        // }
+        
+        if((srt_n > 1) && (srt_a[0] != srt_a[srt_n-1])) {
+            med_f = median_inplace(srt_a, srt_n, 0);
+            mad_f = mad_from_sorted(srt_a, srt_n, med_f);
+            mad = ((double)mad_f) * 1.4826;
+            if (mad < 1e-12) mad = 1e-12;
+            thr = med_f + k_mad * mad;
+            *max_est_err = thr;
+            *ave_est_err = med_f + mad_f;
+            // for (k = 0; k < srt_n; ++k) {
+            //     if((srt_a[k]>>32) <= thr) {
+
+            //     }
+            // }
+            // if(is_dbg) {
+            //     fprintf(stderr, "[M::%s]\tmed_f::%lu,\tmad_f::%lu,\tmad::%f,\tthr::%lu\n", 
+            //         __func__, med_f, mad_f, mad, thr);
+            // }
+                         
+        } else {
+            *max_est_err = *ave_est_err = (srt_a[srt_n-1]); ///e_cut_n = bk;
+        } 
+
+        if((*max_est_err) < mm_cut) {
+            cut = (*max_est_err);
+            for (k = 0; k < srt_n; k++) {
+                if((srt_b[k]>>32) > cut) {
+                    p = &(c_idx[id_a[((uint32_t)srt_b[k])]]);
+                    ovlp_um(*p) += (oe - os);
+                }
+            }
+        }
+    }
+    // if(is_dbg) fprintf(stderr, "[M::%s]\t*max_est_err::%ld\t*ave_est_err::%ld\n", __func__, *max_est_err, *ave_est_err);
+
+    return rr;
+}
+
+
+
 
 int64_t infer_rovlp(ul_ov_t *li, ul_ov_t *lj, uc_block_t *bi, uc_block_t *bj, All_reads *ridx, ma_ug_t *ug)
 {
@@ -23023,7 +23134,7 @@ void gen_ov_lidel_variant(kv_ul_ov_t *cz, overlap_region_alloc *ol, asg64_v *idx
 }
 
 
-uint64_t rcall_lidel_variant(kv_ul_ov_t *cz, overlap_region_alloc *ol, haplotype_evdience_alloc *hp, asg64_v *idx, asg64_v *idz, uint64_t rid)
+uint64_t rcall_lidel_variant(kv_ul_ov_t *cz, overlap_region_alloc *ol, haplotype_evdience_alloc *hp, asg64_v *idx, asg64_v *idz, uint64_t rid, int64_t hom_cov_a, int64_t het_cov_a, int64_t n_hap)
 {
     // fprintf(stderr, "-0-[M::%s] %.*s\tis_match::%u\n", __func__, (int)Get_NAME_LENGTH(R_INF, ol->list[48].y_id), Get_NAME(R_INF, ol->list[48].y_id), ol->list[48].is_match); 
 
@@ -23135,7 +23246,7 @@ uint64_t rcall_lidel_variant(kv_ul_ov_t *cz, overlap_region_alloc *ol, haplotype
             }
         }
 
-        push_idel_info(hp, s, e, svi[k].qn, ov, ovn, sv + svi[k].ts, svi[k].te - svi[k].ts, ol->list, idx->a + ovn, idx->n - ovn, rid, asm_opt.het_cov, asm_opt.hom_cov, asm_opt.polyploidy, 0.333333, 5);
+        push_idel_info(hp, s, e, svi[k].qn, ov, ovn, sv + svi[k].ts, svi[k].te - svi[k].ts, ol->list, idx->a + ovn, idx->n - ovn, rid, het_cov_a, hom_cov_a, n_hap, 0.333333, 5);
 
         ///relabel matched overlaps
         for (p = svi[k].ts; p < svi[k].te; p++) {
@@ -23148,7 +23259,7 @@ uint64_t rcall_lidel_variant(kv_ul_ov_t *cz, overlap_region_alloc *ol, haplotype
     return hp->length;
 }
 
-uint64_t rphase_lidel(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, uint64_t rid, uint64_t hpc_len, uint64_t std_bs)
+uint64_t rphase_lidel(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, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, int64_t hom_cov_a, int64_t het_cov_a, int64_t n_hap)
 {
     hp->length = hp->snp_stat.n = 0; 
     int64_t on = ol->length, k, i, zwn, q[2], t[2]; 
@@ -23180,7 +23291,7 @@ uint64_t rphase_lidel(overlap_region_alloc* ol, All_reads *rref, haplotype_evdie
     on = rphase_lidel_cc(ol, c_idx->a, on, 0.500001, 3, 0.25, 3, rid, buf, idx);
     c_idx->n = on; 
     // fprintf(stderr, "-1-[M::%s]\n", __func__); 
-    return rcall_lidel_variant(c_idx, ol, hp, idx, buf, rid);
+    return rcall_lidel_variant(c_idx, ol, hp, idx, buf, rid, hom_cov_a, het_cov_a, n_hap);
 }
 
 void dbg_cp_cov(overlap_region_alloc* ol, kv_ul_ov_t *c_idx, int64_t ps, int64_t ic)
@@ -23222,8 +23333,207 @@ void gen_qvec_hvec(All_reads *rref, asg8_v *t, uint8_t **hf, uint8_t **qual, uin
     }
 }
 
+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)
+{
+    overlap_region *z; uint64_t m, *wea = NULL, iin = 0, ewn = (ql/wl) + (((ql%wl)>0)?1:0), ewk; ul_ov_t *cp; uint8_t fm;
+    if(ou_a) {
+        wea = ou_a;
+    } else {
+        iin = ewn;
+        ix->n += iin; kv_resize(uint64_t, *ix, ix->n);
+    }
+    int64_t k, i, ixn0 = ix->n, on = ol->length, zwn, q[2], est_e, est_bd; uint64_t tot_e_av = 0, tot_e_bd = 0, tot_cov_l = 0;
+    for (k = 0; k < on; k++) {
+        z = &(ol->list[k]); zwn = z->w_list.n; 
+        if((!zwn) || ((z->is_match != 1) && (z->is_match != 2))) continue;
 
-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)
+        for (i = 0; i < zwn; i++) {
+            if(is_ualn_win(z->w_list.a[i])) continue;
+            q[0] = z->w_list.a[i].x_start; q[1] = z->w_list.a[i].x_end;
+            if(q[1] >= q[0]) {
+                m = ((uint64_t)q[0]); m <<= 32; 
+                m += c_idx->n; kv_push(uint64_t, *ix, m);
+
+                kv_pushp(ul_ov_t, *c_idx, &cp);
+                ovlp_id(*cp) = k; ///ovlp id
+                ovlp_cur_wid(*cp) = i; ///cur id of windows
+                ovlp_cur_xoff(*cp) = z->w_list.a[i].x_start; ///cur xpos
+                ovlp_cur_yoff(*cp) = z->w_list.a[i].y_start; ///cur xpos
+                ovlp_cur_ylen(*cp) = 0;
+                ovlp_cur_coff(*cp) = 0; ///cur cigar off in cur window
+                ovlp_um(*cp) = 0;
+            }
+        }
+    }
+
+    int64_t srt_n = ix->n, s, e, t, os, oe, rm_n, rr; i = 0;
+    radix_sort_bc64(ix->a + ixn0, ix->a + ix->n);
+    for (k = ixn0 + 1, i = ixn0; k <= srt_n; k++) {
+        if (k == srt_n || (ix->a[k]>>32) != (ix->a[i]>>32)) {
+            if(k - i > 1) {
+                for (t = i; t < k; t++) {
+                    cp = &(c_idx->a[(uint32_t)ix->a[t]]);
+                    // s = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start+ovlp_bd(*cp);
+                    // assert(s == (int64_t)(idx->a[i]>>32));
+                    m = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end + 1;
+                    m <<= 32; m += ((uint32_t)ix->a[t]); ix->a[t] = m;
+                    // fprintf(stderr, "[M::%s] s::%ld\tsi::%lu\n", __func__, s, (idx->a[i]>>32));
+                }
+                radix_sort_bc64(ix->a + i, ix->a + k);
+            }
+            i = k;
+        }
+    }
+
+
+    i = ixn0; s = 0; e = wl; e = ((e<=ql)?e:ql); rr = 0; ewk = 0;
+    for (; s < ql; ) {///[s, e)
+        if(rr) {
+            // rr = 0;
+            for (m = rm_n = srt_n; m < ix->n; m++) {
+                cp = &(c_idx->a[ix->a[m]]);
+                q[0] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start;
+                q[1] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1;
+                ///[s, e) && [q[0], q[1])
+                os = MAX(q[0], s); oe = MIN(q[1], e);
+                if(oe > os) {
+                    ix->a[rm_n++] = ix->a[m];
+                    // if(q[1] <= e) rr = 1;
+                }
+            }
+            ix->n = rm_n;
+        }
+
+
+        for (; i < srt_n; ++i) {
+            cp = &(c_idx->a[(uint32_t)ix->a[i]]);
+            q[0] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start;
+            q[1] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1;
+            if(q[0] >= e) break;
+            os = MAX(q[0], s); oe = MIN(q[1], e);
+            if(oe > os) {
+                kv_push(uint64_t, *ix, ((uint32_t)ix->a[i]));
+                // if(q[1] <= e) rr = 1;
+            }
+        }
+
+        m = ix->n + ix->n - srt_n + ix->n - srt_n; kv_resize(uint64_t, *ix, m);
+        // rr = hc_est_robust_rr_diff_gap(ol->list, ql, ix->a + srt_n, ix->n - srt_n, ix->a + ix->n, ix->a + ix->n + ix->n - srt_n, s, e, 64, 5.0, c_idx->a, &est_e, 1);
+        rr = hc_est_robust_rr(ol->list, ql, ix->a + srt_n, ix->n - srt_n, ix->a + ix->n, ix->a + ix->n + ix->n - srt_n, s, e, 2.0, c_idx->a, &est_bd, &est_e, 0);
+        // fprintf(stderr, "[M::%s-0-]\tq::[%ld,%ld)\test_bd::%ld\test_e::%ld\n", __func__, s, e, est_bd, est_e);
+
+        if(!ou_a) wea = ix->a + ixn0 - iin; 
+        if(est_bd == INT64_MAX || est_e == INT64_MAX) {
+            wea[ewk] = (uint64_t)-1;
+        } else {
+            wea[ewk] = (uint32_t)est_bd;  wea[ewk] <<= 32; wea[ewk] |= (uint32_t)est_e; 
+        }
+        ewk++;
+
+        s += wl; e += wl; e = ((e<=ql)?e:ql);
+    }
+    assert(ewn == ewk);
+
+    if(!ou_a) wea = ix->a + ixn0 - iin; 
+    for (ewk = tot_e_av = tot_e_bd = tot_cov_l = 0; ewk < ewn; ewk++) {
+        if(wea[ewk] == ((uint64_t)-1)) continue;
+        s = ewk*wl; e = s + wl;
+        if((int64_t)s >= ql) s = ql; 
+        if((int64_t)e >= ql) e = ql;
+        tot_cov_l += e - s;
+        tot_e_av += (uint32_t)wea[ewk];
+        tot_e_bd += wea[ewk]>>32;
+    }
+    fprintf(stderr, "-a-[M::%s]\test_err::%lu(%f),\tmax_err::%lu(%f),\tql::%ld,\ttot_cov_l::%lu\n", 
+        __func__, tot_e_av, tot_cov_l?((double)tot_e_av)/((double)tot_cov_l):0, tot_e_bd, tot_cov_l?((double)tot_e_bd)/((double)tot_cov_l):0, ql, tot_cov_l);
+
+    // fprintf(stderr, "[M::%s]\tn0::%ld\tc_idx->n::%u\n", __func__, srt_n - ixn0, (uint32_t)c_idx->n);
+
+    for (i = k = ixn0, fm = 0; i < srt_n; i++) {
+        cp = &(c_idx->a[(uint32_t)ix->a[i]]);
+        q[0] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start;
+        q[1] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1;
+        m = q[1] - q[0];
+        if(ovlp_um(*cp) > 0 && ovlp_um(*cp) > (m*0.75)) {
+            fm = 1; continue;
+        }
+        ovlp_um(*cp) = 0;
+        ix->a[k++] = ix->a[i];
+    }
+    srt_n = ix->n = k;
+    // fprintf(stderr, "[M::%s]\tn1::%ld\tc_idx->n::%u\n", __func__, srt_n - ixn0, (uint32_t)c_idx->n);
+
+    if(srt_n > ixn0 && fm) {
+        i = ixn0; s = 0; e = wl; e = ((e<=ql)?e:ql); rr = 0; ewk = 0;
+        for (; s < ql; ) {///[s, e)
+            // fprintf(stderr, "[M::%s-0-]\tq::[%ld,%ld)\tix->n::%u\tsrt_n::%ld\n", __func__, s, e, (uint32_t)ix->n, srt_n);
+            if(rr) {
+                // rr = 0;
+                for (m = rm_n = srt_n; m < ix->n; m++) {
+                    cp = &(c_idx->a[ix->a[m]]);
+                    q[0] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start;
+                    q[1] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1;
+                    ///[s, e) && [q[0], q[1])
+                    os = MAX(q[0], s); oe = MIN(q[1], e);
+                    if(oe > os) {
+                        ix->a[rm_n++] = ix->a[m];
+                        // if(q[1] <= e) rr = 1;
+                    }
+                }
+                ix->n = rm_n;
+            }
+            // fprintf(stderr, "[M::%s-1-]\tq::[%ld,%ld)\tix->n::%u\tsrt_n::%ld\n", __func__, s, e, (uint32_t)ix->n, srt_n);
+
+            for (; i < srt_n; ++i) {
+                cp = &(c_idx->a[(uint32_t)ix->a[i]]);
+                q[0] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start;
+                q[1] = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1;
+                if(q[0] >= e) break;
+                os = MAX(q[0], s); oe = MIN(q[1], e);
+                if(oe > os) {
+                    kv_push(uint64_t, *ix, ((uint32_t)ix->a[i]));
+                    // if(q[1] <= e) rr = 1;
+                }
+            }
+            // fprintf(stderr, "[M::%s-2-]\tq::[%ld,%ld)\tix->n::%u\tsrt_n::%ld\n", __func__, s, e, (uint32_t)ix->n, srt_n);
+
+            m = ix->n + ix->n - srt_n + ix->n - srt_n; kv_resize(uint64_t, *ix, m);
+            // rr = hc_est_robust_rr_diff_gap(ol->list, ql, ix->a + srt_n, ix->n - srt_n, ix->a + ix->n, ix->a + ix->n + ix->n - srt_n, s, e, 64, 5.0, c_idx->a, &est_e, 1);
+            rr = hc_est_robust_rr(ol->list, ql, ix->a + srt_n, ix->n - srt_n, ix->a + ix->n, ix->a + ix->n + ix->n - srt_n, s, e, 2.0, c_idx->a, &est_bd, &est_e, 0);
+
+            if(!ou_a) wea = ix->a + ixn0 - iin; 
+            if(est_bd == INT64_MAX || est_e == INT64_MAX) {
+                wea[ewk] = (uint64_t)-1;
+            } else {
+                wea[ewk] = (uint32_t)est_bd;  wea[ewk] <<= 32; wea[ewk] |= (uint32_t)est_e; 
+            }
+            ewk++;
+
+            s += wl; e += wl; e = ((e<=ql)?e:ql);
+        }
+        assert(ewn == ewk);
+        if(!ou_a) wea = ix->a + ixn0 - iin; 
+        for (ewk = tot_e_av = tot_e_bd = tot_cov_l = 0; ewk < ewn; ewk++) {
+            if(wea[ewk] == ((uint64_t)-1)) continue;
+            s = ewk*wl; e = s + wl;
+            if((int64_t)s >= ql) s = ql; 
+            if((int64_t)e >= ql) e = ql;
+            tot_cov_l += e - s;
+            tot_e_av += (uint32_t)wea[ewk];
+            tot_e_bd += wea[ewk]>>32;
+        }
+        fprintf(stderr, "-b-[M::%s]\test_err::%lu(%f),\tmax_err::%lu(%f),\tql::%ld,\ttot_cov_l::%lu\n", 
+            __func__, tot_e_av, tot_cov_l?((double)tot_e_av)/((double)tot_cov_l):0, tot_e_bd, tot_cov_l?((double)tot_e_bd)/((double)tot_cov_l):0, ql, tot_cov_l);
+    }
+
+    
+
+    ix->n = ixn0 - iin; 
+
+}
+
+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)
 {   
     int64_t on = ol->length, k, i, zwn, q[2]/**, ndp, odp, ms, me, hfs, hfe, hfi**/; 
     uint64_t m, l0, wi, wl0, si, ei, fi; overlap_region *z; ul_ov_t *cp; uint8_t *qhf = NULL, *qual = NULL;
@@ -23397,20 +23707,20 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all
 
     if(!dp) {
         // generate_haplotypes_naive_advance(hap, overlap_list, NULL);
-        generate_haplotypes_naive_HiFi(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)));
+        generate_haplotypes_naive_HiFi(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), hap_cov_match, hap_cov_unmatch);
         // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat);
         // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat);
     } else {
         // gen_rphase_dp(hp, ol, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), dp, idx, buf, rid, q8, tcut, site_sc);
-        gen_rphase_dp_adv(hp, ol, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), dp, idx, buf, rid, qual, tcut, site_sc);
+        gen_rphase_dp_adv(hp, ol, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), dp, idx, buf, rid, qual, tcut, site_sc, hap_cov_match, hap_cov_unmatch, n_hap, het_cov_a, hom_cov_a, hf_rate);
         // if(!site_sc) generate_haplotypes_naive_HiFi_adv(hp, ol, 0.04, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid);//r835
-        if(!site_sc) generate_haplotypes_naive_HiFi_adv_hc(hp, ol, 0.04, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid);///r835
-        else generate_haplotypes_weight(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), 32, 2);
+        if(!site_sc) generate_haplotypes_naive_HiFi_adv_hc(hp, ol, 0.04, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid, hap_cov_match, hap_cov_unmatch);///r835
+        else generate_haplotypes_weight(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), 32, 2, hap_cov_match, hap_cov_unmatch);
     }
 
     if(lindel) {
-        if(rphase_lidel(ol, rref, hp, qu, tu, c_idx, idx, buf, bd, wl, ql, occ_thres, rid, hpc_len, std_bs)) {
-            generate_haplotypes_sv(hp, ol, rid);
+        if(rphase_lidel(ol, rref, hp, qu, tu, c_idx, idx, buf, bd, wl, ql, occ_thres, rid, hpc_len, std_bs, hom_cov_a, het_cov_a, n_hap)) {
+            generate_haplotypes_sv(hp, ol, rid, hap_cov_match, hap_cov_unmatch);
         }
     }
 }
@@ -33423,6 +33733,10 @@ void rescue_cu_aln_adv(gen_hc_aln_t *ez, int64_t ql, uint64_t *wcut, uint64_t wc
     }
 }
 
+// void est_hc_r_alin_err(overlap_region_alloc* ol, uint32_t *werr, uint64_t werr_n, uint64_t wl)
+// {
+//     memset(werr, -1, sizeof((*werr))*werr_n);
+// }
 
 
 ///need to consider coverage, this information is missing right now (currently only use numbers)
diff --git a/Correct.h b/Correct.h
index 83268bf..1e36554 100644
--- a/Correct.h
+++ b/Correct.h
@@ -1427,6 +1427,8 @@ typedef struct
     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, 
@@ -1446,12 +1448,14 @@ 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);
+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);
+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);
 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;
@@ -1468,6 +1472,7 @@ inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl,
 #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
diff --git a/Process_Read.h b/Process_Read.h
index c728f71..2205242 100644
--- a/Process_Read.h
+++ b/Process_Read.h
@@ -144,6 +144,8 @@ 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;
 
diff --git a/anchor.cpp b/anchor.cpp
index 6a69b87..73784de 100644
--- a/anchor.cpp
+++ b/anchor.cpp
@@ -985,7 +985,7 @@ uint32_t *low_occ)
 
 
 uint64_t minimizers_qgen0(ha_abuf_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag, 
-void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint64_t ti_cut)
+void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint64_t ti_cut, uint8_t flt_chh)
 {
 	// fprintf(stderr, "+[M::%s]\n", __func__);
 	uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0; int n, j; ha_mz1_t *z; seed1_t *s; 
@@ -1009,7 +1009,6 @@ void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_m
 	}
 
 	for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
-
 		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
 		ab->seed[i].n = n;
 		ab->n_a += n;
@@ -1025,6 +1024,7 @@ void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_m
 		z = &ab->mz.a[i]; s = &ab->seed[i];
 		for (j = 0; j < s->n; ++j) {
 			const ha_idxpos_t *y = &s->a[j];
+			if((flt_chh) && (y->rid < ti_cut)) continue;///ONT
 			anchor1_t *an = &ab->a[k++];
 			uint8_t rev = z->rev == y->rev? 0 : 1;
 			an->other_off = rev?((uint32_t)-1)-1-(y->pos+1-y->span):y->pos;
@@ -1035,6 +1035,7 @@ void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_m
 			an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off;
 		}
 	}
+	ab->n_a = k;
 
 	// copy over to _cl_
 	if (ab->m_a >= (uint64_t)cl->size) {
@@ -2682,7 +2683,7 @@ void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz
 	int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip;
 	set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
 	// minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ);
-	minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ, ((uint64_t)-1));
+	minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ, ((uint64_t)-1), 0);
 	// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
 	// lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
 	///no need to sort here, overlap_list has been sorted at lchain_gen
@@ -2690,14 +2691,14 @@ void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz
 }
 
 void h_ec_lchain_hybrid(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres_h, double bw_thres_l,  
-	int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint64_t ti_cut, uint8_t is_raw_chain)
+	int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint64_t ti_cut, uint8_t flt_chh, uint8_t is_raw_chain)
 {
 	extern void *ha_flt_tab;
 	extern ha_pt_t *ha_idx;
 	int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip; uint64_t tcut_n = 0;
 	set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
 	// minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ);
-	tcut_n = minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ, ti_cut);
+	tcut_n = minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ, ti_cut, flt_chh);
 	// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
 	// lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
 	///no need to sort here, overlap_list has been sorted at lchain_gen
diff --git a/ecovlp.cpp b/ecovlp.cpp
index c402145..41be62e 100644
--- a/ecovlp.cpp
+++ b/ecovlp.cpp
@@ -125,13 +125,18 @@ cc_v scc = {0, 0, NULL, NULL};
 cc_v scb = {0, 0, NULL, NULL};
 cc_v sca = {0, 0, NULL, NULL};
 
+typedef struct {uint64_t p, pn, pm, tov, tov_size, tqn; asg64_v *idx; ma_hit_t_alloc *pf;} tsrt_v_buf;
+typedef struct {uint64_t p, pn, pm, rid, tot, chunk_size, tqn, n_thr; uint64_t n_ov, n_bl; ma_hit_t_alloc *pf;} tsrt_v_m;
+
+
+
 typedef struct {size_t n, m; char *a; UC_Read z; asg8_v q;} sl_v;
 
 
 void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, 
 								 int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint8_t is_raw_chain);
 void h_ec_lchain_hybrid(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres_h, double bw_thres_l,  
-                                int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint64_t ti_cut, uint8_t is_raw_chain);
+                                int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint64_t ti_cut, uint8_t flt_chh, uint8_t is_raw_chain);
 
 void h_ec_lchain_amz(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, 
 								 int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w);
@@ -2407,14 +2412,15 @@ void print_debug_ovlp_cigar(overlap_region_alloc* ol, asg64_v* idx, kv_ul_ov_t *
     }    
 }
 
-uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, uint64_t wl, int64_t ql, uint64_t occ_tot, double occ_exact, overlap_region *aux_o, asg32_v* b32, cns_gfa *cns, uint64_t cns_g_wl, uint32_t rid, uint64_t tcut, asg64_v *hf_idx)
+uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, uint64_t wl, int64_t ql, uint64_t occ_tot, double occ_exact, overlap_region *aux_o, asg32_v* b32, cns_gfa *cns, uint64_t cns_g_wl, uint32_t rid, uint64_t tcut, 
+    uint64_t tot_ont_b, uint64_t tot_hf_b, uint64_t ont_rate_w, uint64_t hf_rate_w, uint64_t hf_rate_w_max, asg64_v *hf_idx)
 {
     int64_t on = ol->length, k, i, zwn, q[2]; cns->cns_g_wl = cns_g_wl;
     uint64_t m, *ra, rn, nec = 0, n_id, l_nid, p[2], li; uint64_t o_rate = ((uint64_t)-1), h_rate = ((uint64_t)-1); overlap_region *z; ul_ov_t *cp; 
     bit_extz_t ez; uint64_t ci; uint32_t cl; uint16_t c, hf;
     if(hf_idx) {
-        o_rate = asm_opt.ont_rate;
-        h_rate = ceil(((double)rref->tr[0])/((double)rref->tr[1]))*asm_opt.hf_rate; h_rate = MAX(h_rate, asm_opt.hf_rate_max);
+        o_rate = ont_rate_w;
+        h_rate = ceil(((double)tot_ont_b)/((double)tot_hf_b))*hf_rate_w; h_rate = MAX(h_rate, hf_rate_w_max);
         if(o_rate == 0) h_rate = 1;
     }
     
@@ -2710,6 +2716,132 @@ uint64_t extract_max_exact(overlap_region *z, asg16_v *ec, /**UC_Read *qu, UC_Re
     return 0;
 }
 
+int64_t extract_max_exact_smp(asg16_v *in, int64_t xs0, int64_t xe0, int64_t ys0, int64_t ye0, uint32_t *rxs, uint32_t *rxe, uint32_t *rys, uint32_t *rye)
+{
+    *rxs = *rxe = *rys = *rye = 0; 
+    if((xe0 <= xs0) || (ye0 <= ys0)) return 0;
+    int64_t ok = 0, nk = 0, ck = 0, cn = in->n, wo[2], wn[2], os, oe; uint16_t op, bq, bt; uint32_t cl; uint64_t ovlp;
+    while (ck < cn && ok < xe0) {
+        wo[0] = ok; wn[0] = nk; 
+        ck = pop_trace_bp_f(in, ck, &op, &bq, &bt, &cl);
+        if(op != 2) ok += cl;
+        if(op != 3) nk += cl;
+        wo[1] = ok; wn[1] = nk;  
+
+        if(op == 0) {
+            os = MAX(xs0, wo[0]); oe = MIN(xe0, wo[1]);
+            ovlp = ((oe>os)? (oe-os):0);
+            if((ovlp > 0) && (ovlp > (*rxe) - (*rxs))) {
+                // fprintf(stderr, "[M::%s]\to::[%ld,%ld)\n", __func__, os, oe);
+                (*rxs) = wn[0] + os - wo[0]; 
+                (*rxe) = wn[0] + oe - wo[0]; 
+
+                (*rys) = ys0 + os - xs0; 
+                (*rye) = ys0 + oe - xs0; 
+            }
+        }
+    }
+
+    return (*rxe) - (*rxs);
+}
+
+void push_ne_ovlp_flt(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, All_reads* R_INF, int64_t ql, asg16_v *ec, uint64_t tqn, uint16_t hf_only)
+{
+    uint64_t k, n; uint8_t el; ma_hit_t *z; uint32_t rxs, rxe, rys, rye; int64_t tl, qs, qe, ts, te, qr, tr;
+    if(hf_only) {
+        for (k = n = 0; k < paf->length; k++) {
+            if(paf->buffer[k].tn >= tqn) continue;///no existing HiFi-to-HiFi, get them from ov
+            z = &(paf->buffer[k]); el = z->el; z->el = 0;
+            if(el) {
+                if((ec) && (extract_max_exact_smp(ec, (uint32_t)(z->qns), z->qe, z->ts, z->te, &rxs, &rxe, &rys, &rye) > 0)) {
+                    z->qns = ov->list[k].x_id; z->qns = z->qns << 32; z->qns = z->qns | (uint64_t)(rxs); z->qe = rxe;
+                    z->ts = rys; z->te = rye;
+                    z->el = 1;
+                }
+
+                if(z->el == 0) {///extend to normal
+                    tl = Get_READ_LENGTH((*R_INF), Get_tn(*z));
+                    qs = ((uint32_t)(z->qns)); qe = z->qe; ts = z->ts; te = z->te;
+                    if(qs >= ql) {qs = ql;} if(qe > ql) {qe = ql;} 
+                    if(ts >= tl) {ts = tl;} if(te > tl) {te = tl;} 
+
+                    if(qs <= ts) {
+                        ts -= qs; qs = 0;
+                    } else {
+                        qs -= ts; ts = 0;
+                    }
+                
+                    qr = ql - qe; tr = tl - te;
+                    if(qr <= tr) {
+                        qe = ql; te += qr;        
+                    } else {
+                        te = tl; qe += tr; 
+                    }
+                
+                    z->qns >>= 32; z->qns <<= 32; z->qns |= ((uint64_t)(qs)); z->qe = qe;
+                    z->ts = ts; z->te = te;
+                }
+            }
+            if((z->qe <= ((uint32_t)(z->qns))) || (z->te <= z->ts)) continue;
+            paf->buffer[n++] = paf->buffer[k];
+        }
+        paf->length = n;
+    } else {
+        paf->length = 0;
+    }
+
+    for (k = 0, n = paf->length; k < ov->length; k++) {
+        if(ov->list[k].is_match == flag) n++;
+    }
+
+    if(n > paf->size) {
+        paf->size = n;
+        REALLOC(paf->buffer, paf->size);
+    }
+
+
+    for (k = 0; k < ov->length; k++) {
+        if(ov->list[k].is_match == flag) {
+            // fprintf(stderr, "@%s\tSN:%.*s(id::%u)\terr::%u\n", flag==1?"SQ":"RQ", (int32_t)Get_NAME_LENGTH((*R_INF), ov->list[k].y_id), Get_NAME((*R_INF), ov->list[k].y_id), ov->list[k].y_id, ov->list[k].non_homopolymer_errors);
+
+            z = &(paf->buffer[paf->length++]);
+
+            z->qns = ov->list[k].x_id;
+            z->qns = z->qns << 32;
+            z->tn = ov->list[k].y_id;
+
+            z->qns = z->qns | (uint64_t)(ov->list[k].x_pos_s);
+            z->qe = ov->list[k].x_pos_e + 1;
+            z->ts = ov->list[k].y_pos_s;
+            z->te = ov->list[k].y_pos_e + 1;
+
+            ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
+            z->rev = ov->list[k].y_pos_strand;
+
+            z->bl = Get_READ_LENGTH((*R_INF), ov->list[k].y_id);
+            z->ml = ov->list[k].strong;
+            z->no_l_indel = ov->list[k].without_large_indel;
+
+            z->el = 0;
+            if(ec) {
+                extract_max_exact(&ov->list[k], ec, /**qu, tu,**/ &rxs, &rxe, &rys, &rye);
+                // z->el = 0;
+                // fprintf(stderr, "[M::%s]\tq::[%u,\t%u)\tt::[%u,\t%u)\teq::[%u,\t%u)\tet::[%u,\t%u)\n", __func__, ov->list[k].x_pos_s, ov->list[k].x_pos_e + 1, ov->list[k].y_pos_s, ov->list[k].y_pos_e + 1, rxs, rxe, rys, rye);
+                if(rxe > rxs) {
+                    z->qns = ov->list[k].x_id;
+                    z->qns = z->qns << 32;
+                    z->qns = z->qns | (uint64_t)(rxs);
+                    z->qe = rxe;
+                    z->ts = rys;
+                    z->te = rye;
+
+                    z->el = 1;
+                }
+            }
+        }
+    }
+}
+
 void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, All_reads* R_INF, asg16_v *ec/**, uint64_t qid, UC_Read *qu, UC_Read *tu**/)
 {
     // if(qu && tu) {
@@ -2725,6 +2857,114 @@ void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag,
         REALLOC(paf->buffer, paf->size);
     }
 
+    for (k = paf->length = 0; k < ov->length; k++) {
+        if(ov->list[k].is_match == flag) {
+            // fprintf(stderr, "@%s\tSN:%.*s(id::%u)\terr::%u\n", flag==1?"SQ":"RQ", (int32_t)Get_NAME_LENGTH((*R_INF), ov->list[k].y_id), Get_NAME((*R_INF), ov->list[k].y_id), ov->list[k].y_id, ov->list[k].non_homopolymer_errors);
+
+            z = &(paf->buffer[paf->length++]);
+
+            z->qns = ov->list[k].x_id;
+            z->qns = z->qns << 32;
+            z->tn = ov->list[k].y_id;
+
+            z->qns = z->qns | (uint64_t)(ov->list[k].x_pos_s);
+            z->qe = ov->list[k].x_pos_e + 1;
+            z->ts = ov->list[k].y_pos_s;
+            z->te = ov->list[k].y_pos_e + 1;
+
+            ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
+            z->rev = ov->list[k].y_pos_strand;
+
+            z->bl = Get_READ_LENGTH((*R_INF), ov->list[k].y_id);
+            z->ml = ov->list[k].strong;
+            z->no_l_indel = ov->list[k].without_large_indel;
+
+            z->el = 0;
+            if(ec) {
+                extract_max_exact(&ov->list[k], ec, /**qu, tu,**/ &rxs, &rxe, &rys, &rye);
+                // z->el = 0;
+                // fprintf(stderr, "[M::%s]\tq::[%u,\t%u)\tt::[%u,\t%u)\teq::[%u,\t%u)\tet::[%u,\t%u)\n", __func__, ov->list[k].x_pos_s, ov->list[k].x_pos_e + 1, ov->list[k].y_pos_s, ov->list[k].y_pos_e + 1, rxs, rxe, rys, rye);
+                if(rxe > rxs) {
+                    z->qns = ov->list[k].x_id;
+                    z->qns = z->qns << 32;
+                    z->qns = z->qns | (uint64_t)(rxs);
+                    z->qe = rxe;
+                    z->ts = rys;
+                    z->te = rye;
+
+                    z->el = 1;
+                }
+            }
+        }
+    }
+}
+
+
+void pull_ovlp_syn(ma_hit_t_alloc* paf, uint64_t tqn)
+{
+    uint64_t k, m;
+    for (k = m = 0; k < paf->length; k++) {
+        if(paf->buffer[k].tn < tqn) continue; ///ONT
+        paf->buffer[m++] = paf->buffer[k];
+    }
+    paf->length = m;
+}
+
+uint64_t gen_ne_ovlp_hf_region(asg64_v *idx, int64_t het_cov, double hf_rate)
+{
+    if(idx->n == 0) return 0;///no HiFi reads
+
+    radix_sort_ec64(idx->a, idx->a + idx->n); 
+
+    int64_t min_dp = ceil(((double)het_cov)*hf_rate*0.6); if(min_dp < 6) min_dp = 6;
+    int64_t k, dp, old_dp, st, ed, s0, e0, n = idx->n; uint64_t m;
+    for (k = dp = old_dp = st = ed = 0; k < n; k++) {
+        old_dp = dp;
+        if (idx->a[k]&1) {--dp;}///if a[j] is qe
+        else {++dp;}
+
+        ed = idx->a[k]>>1;
+        if((ed > st) && (old_dp >= min_dp)) {
+            s0 = e0 = -1;
+            if(idx->n > (uint64_t)n) {
+                s0 = idx->a[idx->n-1]>>32;
+                e0 = (uint32_t)idx->a[idx->n-1];
+            }
+            if(e0 == st) {
+                m = s0; m <<= 32; m |= (uint64_t)ed;
+                idx->a[idx->n-1] = m;
+            } else {
+                assert(st > e0);
+                m = st; m <<= 32; m |= (uint64_t)ed;
+                kv_push(uint64_t, *idx, m);
+            }
+        }
+        st = ed;
+    }
+
+    return idx->n - n;
+}
+
+void push_ne_ovlp_syn(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, All_reads* R_INF, asg16_v *ec, asg64_v *idx, uint64_t tqn, int64_t het_cov, double hf_rate) 
+{
+    uint64_t k, s, e, *sa = NULL, sn, sk; ma_hit_t *z; uint32_t rxs, rxe, rys, rye;
+    for (k = idx->n = sn = 0; k < ov->length; k++) {
+        if(ov->list[k].is_match == flag) {
+            if(ov->list[k].y_id >= tqn) {///HiFi
+                s = ov->list[k].x_pos_s; s = s<<1; kv_push(uint64_t, *idx, s);
+                e = ov->list[k].x_pos_e + 1; e = (e<<1)|1; kv_push(uint64_t, *idx, e);                
+                // idx->n += 2;
+            }
+            sn++;
+        }
+    }
+
+    if(sn > paf->size) {
+        paf->size = sn;
+        REALLOC(paf->buffer, paf->size);
+    }
+    sn = gen_ne_ovlp_hf_region(idx, het_cov, hf_rate); sa = idx->a + idx->n - sn;
+
     for (k = paf->length = 0; k < ov->length; k++) {
         if(ov->list[k].is_match == flag) {
             // fprintf(stderr, "@%s\tSN:%.*s(id::%u)\terr::%u\n", flag==1?"SQ":"RQ", (int32_t)Get_NAME_LENGTH((*R_INF), ov->list[k].y_id), Get_NAME((*R_INF), ov->list[k].y_id), ov->list[k].y_id, ov->list[k].non_homopolymer_errors);
@@ -2764,8 +3004,25 @@ void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag,
             }
         }
     }
+    if(sn <= 0) return;///no HiFi
+    
+    uint64_t s0, e0, st, ed;
+    for (k = 0; k < paf->length; k++) {
+        z = &(paf->buffer[k]); 
+        if(z->tn >= tqn) {///HiFi
+            s0 = (uint32_t)z->qns; e0 = z->qe;
+            for (sk = 0; sk < sn; sk++) {
+                st = sa[sk]>>32; ed = (uint32_t)sa[sk];
+                if(st <= s0 && ed >= e0) {
+                    z->bl = 0x7FFFFFFF; break;
+                }
+                if(st >= s0) break;
+            }
+        }
+    }
 }
 
+
 void push_ff_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, All_reads* R_INF, uint64_t *cnt)
 {
     // if(qu && tu) {
@@ -3763,7 +4020,7 @@ void init_gen_hc_aln_t(gen_hc_aln_t *ez, overlap_region_alloc *ol, Candidates_li
     asg8_v *hpz, double e_rate_l, double e_rate_h, int64_t wl_l, int64_t wl_h, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in,
     int8_t chem_drop_l, int8_t chem_drop_h, double align_gap_rate_l, double align_gap_rate_h, int64_t align_gap_max_l, int64_t align_gap_max_h,
     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 ave_cov_min, uint64_t ocw, uint64_t t_cut, uint64_t hom_cov_a)
 {
     ez->ol = ol;
     ez->cl = cl;
@@ -3808,6 +4065,8 @@ void init_gen_hc_aln_t(gen_hc_aln_t *ez, overlap_region_alloc *ol, Candidates_li
 
     ez->t_cut = t_cut;
     ez->hpz = hpz;
+
+    ez->hom_cov_a = hom_cov_a;
 }
 
 uint64_t cal_aln_bs(overlap_region_alloc *ol)
@@ -3819,13 +4078,22 @@ uint64_t cal_aln_bs(overlap_region_alloc *ol)
     return tot;
 }
 
+
+#define set_ec_cov(het_i, hom_i, het_s, n_hap, het_r, hom_r) do {\
+    (het_r) = (het_i); (hom_r) = (hom_i);\
+    if((het_s) >= 0) (het_r) = (het_s);\
+    if(((het_r) < 0) || ((het_r) > ((hom_r)/(n_hap)))) {(het_r) = (hom_r)/(n_hap);}\
+} while (0)
+
+
+
 static void worker_hap_ec(void *data, long i, int tid)
 {
 	ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
     uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
-    uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; ///gen_hc_aln_t ez;
+    uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; int64_t het_a, hom_a; ///gen_hc_aln_t ez; 
     overlap_region *aux_o = NULL/**, *rse_o = NULL**/; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t tot_b = 0; double tt0 = 0, tt1 = 0;
-    b->v8q.n = b->v8t.n = 0;
+    b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a);
 
     // if((i != 733166) && (i != 858708) && (i != 858732) && (i != 859819) && (i != 859899) && (i != 863486) && (i != 872165) && (i != 899887) && (i != 902298) && 
     //    (i != 906808) && (i != 946173) && (i != 952685) && (i != 983977) && (i != 1000227) && (i != 1011228) && (i != 1042858) && (i != 1045860) && (i != 1118558) && 
@@ -3874,14 +4142,11 @@ static void worker_hap_ec(void *data, long i, int tid)
 
     // if(i != 339646) return;
     // if(i!=854835) return;
-    
+    // if(i != 533) return;
     
     
     // if(i % 100000 == 0) fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i);
-    // if (memcmp("c42804f3-0e13-43a0-8a71-b91b40accf9a", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
-    // if (memcmp("b2e68ecf-381a-439c-b676-c1e6831d6acf", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
-    // if (memcmp("e3f3f43a-e200-4cac-8acd-3f85428f3811", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
-    // if (memcmp("4e144e93-4653-4ebf-8920-7943e378cf9a", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
+    // if (memcmp("6c55c5f1-e86d-4065-bbf7-68b24a995bee", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
     //     fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i);
     // } else {
     //     return;
@@ -3983,11 +4248,14 @@ static void worker_hap_ec(void *data, long i, int tid)
 
     copy_asg_arr(buf0, b->sp); 
     //site_sc: r765 -> r766: 1 -> 0
-    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32);
+    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
+    asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0);
     copy_asg_arr(b->sp, buf0);
     ///for debug indel
     // stderr_phase_ovlp(&b->olist);
 
+    // est_rep_err_rate(&b->olist, &b->v64, &b->pidx, qlen, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), NULL);
+
     if(DBG_TIME && dbg_a) {
         tt1 = yak_realtime_0();
         dbg_a[i].phs_tm = tt1 - tt0;
@@ -4000,7 +4268,8 @@ static void worker_hap_ec(void *data, long i, int tid)
 
 
     copy_asg_arr(buf0, b->sp);
-    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), NULL);
+    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), 
+    R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL);
     copy_asg_arr(b->sp, buf0);
 
     if(DBG_TIME && dbg_a) {
@@ -4104,9 +4373,9 @@ static void worker_hap_ec_ss(void *data, long i, int tid)
 {
 	ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
     uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
-    uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; ///gen_hc_aln_t ez;
+    uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; int64_t het_a, hom_a; ///gen_hc_aln_t ez;
     overlap_region *aux_o = NULL/**, *rse_o = NULL**/; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t /**tot_b = 0,**/ i0 = i, prt_n0; 
-    b->v8q.n = b->v8t.n = 0;
+    b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a);
 
     // if(((uint64_t)i) < dbgss->fn) return;
     i = (uint32_t)dbgss->fa[i0]; prt_n0 = dbgss->spt_mul[tid].n;
@@ -4193,7 +4462,8 @@ static void worker_hap_ec_ss(void *data, long i, int tid)
 
     copy_asg_arr(buf0, b->sp); 
     //site_sc: r765 -> r766: 1 -> 0
-    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32);
+    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
+    asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0);
     copy_asg_arr(b->sp, buf0);
     ///for debug indel
     // stderr_phase_ovlp(&b->olist);
@@ -4204,7 +4474,8 @@ static void worker_hap_ec_ss(void *data, long i, int tid)
 
 
     copy_asg_arr(buf0, b->sp);
-    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), NULL);
+    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), 
+    R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL);
     copy_asg_arr(b->sp, buf0);
 
 
@@ -4227,15 +4498,19 @@ static void worker_hap_ec_ss(void *data, long i, int tid)
 static void worker_hap_ec_hybrid(void *data, long i, int tid)
 {
 	ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
-    uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
+    uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); int64_t het_a, hom_a;
     uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; double bw_h, bw_l, e_h, e_l;
     gen_hc_aln_t ez; overlap_region *aux_o = NULL, *rse_o = NULL; asg64_v buf0, buf1; uint64_t qlen = 0, qw = 0, qid = i; //uint64_t sk[2], ek[2], fn, qid = i, nec;
     if(qid < R_INF.tqn) {///ont
         bw_h = 0.05; bw_l = 0.035; e_h = asm_opt.max_ov_diff_ec; e_l = (asm_opt.max_ov_diff_ec + asm_opt.max_ov_diff_ec_sec)/2;
     } else { ///HiFi
         bw_h = 0.035; bw_l = 0.02; e_h = (asm_opt.max_ov_diff_ec + asm_opt.max_ov_diff_ec_sec)/2; e_l = asm_opt.max_ov_diff_ec_sec;
+        if(R_INF.is_syn) {///remove all HiFi-to-ONT overlaps as we will get them later
+            pull_ovlp_syn(&(R_INF.paf[i]), R_INF.tqn); pull_ovlp_syn(&(R_INF.reverse_paf[i]), R_INF.tqn);
+            return;
+        }
     }
-    b->v8q.n = b->v8t.n = 0;
+    b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a);
 
     // if(i != 5966) return;
 
@@ -4260,7 +4535,7 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
     // if(qlen <= 0) return;
 
     h_ec_lchain_hybrid(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, bw_h, bw_l,
-        /**((asm_opt.is_ont)?(0.05):(0.02)),**/ asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, R_INF.tqn, 1);///ONT high error
+        /**((asm_opt.is_ont)?(0.05):(0.02)),**/ asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, R_INF.tqn, 0, 1);///ONT high error
 
     // fprintf(stderr, "-b-[M::%s] rid::%ld\n", __func__, i);
 
@@ -4275,7 +4550,7 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
         ((qid < R_INF.tqn)?(1):(0)), 1, (qid < R_INF.tqn)?(0.006):(-1), 0.006, (qid < R_INF.tqn)?(64):(-1), 64, 
         (qid < R_INF.tqn)?(512):(0), (qid < R_INF.tqn)?(6):(0), (qid < R_INF.tqn)?(1.5):(-1), (qid < R_INF.tqn)?(0.1):(-1), NULL, 
         &(b->v32), &buf0, b->ab, (asm_opt.max_n_chain>0)?(asm_opt.max_n_chain):(1), ((asm_opt.max_n_chain*HC_MF_R)>0)?(asm_opt.max_n_chain*HC_MF_R):1, 
-        asm_opt.chn_occ, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), qw, R_INF.tqn); 
+        asm_opt.chn_occ, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), qw, R_INF.tqn, asm_opt.hom_cov); 
         // gen_hc_r_alin_ea_adv(&ez);
     gen_hc_r_alin_ea_adv_flt(&ez);
     copy_asg_arr(b->sp, buf0);
@@ -4317,15 +4592,17 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
     // prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length);
 
     copy_asg_arr(buf0, b->sp); 
-    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32);
+    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32,
+    asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1]))));
     copy_asg_arr(b->sp, buf0);
     ///for debug indel
-    // stderr_phase_ovlp(&b->olist);
+    // if(i == 23863) stderr_phase_ovlp(&b->olist);
 
     dedup_chains(&b->olist);
 
     copy_asg_arr(buf0, b->sp); copy_asg_arr(buf1, b->hap.snp_srt);
-    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, R_INF.tqn, &buf1);
+    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, R_INF.tqn, 
+        R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, &buf1);
     copy_asg_arr(b->sp, buf0); copy_asg_arr(b->hap.snp_srt, buf1);
 
     push_nec_re(aux_o, &(scc.a[i]));
@@ -4335,8 +4612,8 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
     //     // b->olist.length = 0;
     //     fprintf(stderr, "[M::%s] rid::%ld\t%.*s\n\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
     // }
-
-    push_ne_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, &(scc.a[i])/**, i, &b->self_read, &b->ovlp_read**/);
+    if((qid < R_INF.tqn) && (R_INF.is_syn == 1)) push_ne_ovlp_syn(&(R_INF.paf[i]), &b->olist, 1, &R_INF, &(scc.a[i]), &b->v64, R_INF.tqn, het_a, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1])))/**, i, &b->self_read, &b->ovlp_read**/);
+    else push_ne_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, &(scc.a[i])/**, i, &b->self_read, &b->ovlp_read**/);
     push_ne_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, NULL/**, i, NULL, NULL**/);
 
 
@@ -4352,6 +4629,172 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
     // exit(0);
 }
 
+uint32_t is_HiFi_only_gen(uint64_t rid, ma_hit_t_alloc *paf, asg64_v *b, uint64_t tqn, int64_t ql, double cut)
+{
+    uint64_t k, s, e; ma_hit_t *z; int64_t tl, qs, qe, ts, te, qr, tr;
+    for (k = b->n = 0; k < paf->length; k++) {
+        z = &(paf->buffer[k]); 
+        if(z->tn >= tqn) continue; ///no HiFi-to-HiFi
+        s = (uint32_t)z->qns; e = z->qe;
+        if(z->el) {
+            tl = Get_READ_LENGTH(R_INF, Get_tn(*z));
+            qs = ((uint32_t)(z->qns)); qe = z->qe; 
+            ts = z->ts; te = z->te;
+            if(qs >= ql) {qs = ql;} if(qe > ql) {qe = ql;} 
+            if(ts >= tl) {ts = tl;} if(te > tl) {te = tl;} 
+
+            if(qs <= ts) {
+                ts -= qs; qs = 0;
+            } else {
+                qs -= ts; ts = 0;
+            }
+
+            qr = ql - qe; tr = tl - te;
+            if(qr <= tr) {
+                qe = ql; te += qr;        
+            } else {
+                te = tl; qe += tr; 
+            }
+            s = qs; e = qe;
+            if(e <= s) continue;
+        }
+        // if(rid == 23863) {
+        //     fprintf(stderr, "[M::%s]\ts::%lu(%u)\te::%lu(%u)\tbl::%u\tz->el::%u\n", __func__, s, (uint32_t)z->qns, e, z->qe, z->bl, z->el);
+        // }
+        s <<= 2; if(z->bl == 0x7FFFFFFF) s |= 2;
+        e <<= 2; if(z->bl == 0x7FFFFFFF) e |= 2; e |= 1;
+        kv_push(uint64_t, *b, s); kv_push(uint64_t, *b, e);
+    }
+    radix_sort_ec64(b->a, b->a + b->n);
+
+    int64_t dp, old_dp, dp_h, old_dp_h, st, ed;
+    for (k = dp = old_dp = dp_h = old_dp_h = st = ed = 0; k < b->n; k++) {
+        old_dp = dp; old_dp_h = dp_h;
+        if (b->a[k]&1) {
+            --dp;///if a[j] is qe
+            if(b->a[k]&2) --dp_h;
+        } else {
+            ++dp;
+            if(b->a[k]&2) ++dp_h;
+        }
+
+        ed = b->a[k]>>2;
+        if(ed > st) {
+            if((old_dp <= 0) || (old_dp_h <= 0)) return 0;
+            if(old_dp_h < (old_dp*cut)) return 0;
+
+        }
+        st = ed;
+    }
+
+
+    ed = ql; old_dp = dp; old_dp_h = dp_h;
+    if(ed > st) {
+        if((old_dp <= 0) || (old_dp_h <= 0)) return 0;
+        if(old_dp_h < (old_dp*cut)) return 0;
+    }
+    
+    return 1;
+}
+
+static void worker_hap_ec_hybrid_sync(void *data, long i, int tid)
+{
+	ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); i += R_INF.tqn; 
+    uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); int64_t max_n_chain_a = asm_opt.max_n_chain;
+    uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; double bw_h, bw_l, e_h, e_l, hf_rate; uint8_t hf_only = 0; int64_t het_a, hom_a;
+    gen_hc_aln_t ez; overlap_region *aux_o = NULL, *rse_o = NULL; asg64_v buf0, buf1; uint64_t qlen = 0, qw = 0/**, qid = i**/; //uint64_t sk[2], ek[2], fn, qid = i, nec;
+
+    hf_only = is_HiFi_only_gen(i, &(R_INF.paf[i]), &b->v64, R_INF.tqn, Get_READ_LENGTH(R_INF, i), 0.666666);
+    // if(i == 23863) {
+    //     fprintf(stderr, "[M::%s]\tqid::%ld\thf_only::%u\trl::%lu\n", __func__, i, hf_only, Get_READ_LENGTH(R_INF, i));
+    //     // fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%lu\tq::[%u,%u)\t%c\t%.*s(tid::%u)\t\ttl::%lu\tt::[%u,%u)\tel::%u\n", __func__,
+    //     //     (int32_t)Get_NAME_LENGTH(R_INF, (z->buffer[k].qns>>32)), Get_NAME(R_INF, (z->buffer[k].qns>>32)), (uint32_t)(z->buffer[k].qns>>32), Get_READ_LENGTH(R_INF, (z->buffer[k].qns>>32)), 
+    //     //     (uint32_t)z->buffer[k].qns, z->buffer[k].qe, "+-"[z->buffer[k].rev], 
+    //     //     (int32_t)Get_NAME_LENGTH(R_INF, z->buffer[k].tn), Get_NAME(R_INF, z->buffer[k].tn), z->buffer[k].tn, Get_READ_LENGTH(R_INF, (z->buffer[k].tn)),  
+    //     //     z->buffer[k].ts, z->buffer[k].te, z->buffer[k].el?1:0);
+    // }
+    bw_h = 0.035; bw_l = 0.02; e_h = (asm_opt.max_ov_diff_ec + asm_opt.max_ov_diff_ec_sec)/2; e_l = asm_opt.max_ov_diff_ec_sec;
+    b->v8q.n = b->v8t.n = 0; hf_rate = ((double)R_INF.tr[1])/(((double)R_INF.tr[1]) + ((double)R_INF.tr[0]));
+    set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a);
+    if(hf_only) {
+        max_n_chain_a = ceil(hf_rate*max_n_chain_a); 
+        if(max_n_chain_a < 6) max_n_chain_a = 6; 
+        if(max_n_chain_a > asm_opt.max_n_chain) max_n_chain_a = asm_opt.max_n_chain;
+
+        het_a = ceil(hf_rate*het_a); 
+        if(het_a < 5) het_a = 5;
+        if(het_a > asm_opt.hom_cov) het_a = asm_opt.hom_cov;
+
+        hom_a = ceil(hf_rate*hom_a); 
+        if(hom_a < 5) hom_a = 5;
+        if(hom_a > asm_opt.het_cov) hom_a = asm_opt.het_cov;
+    }
+    if(max_n_chain_a <= 0) max_n_chain_a = 1;
+
+    // if(i != 5966) return;    
+
+    // debug_retrive_bqual(D, &b->v8t, i, 256); return;
+
+    recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; 
+    qw = ((qlen < (COV_W_AC<<1))?(qlen>>1):(COV_W_AC)); if(!qw) qw = 1;
+
+    h_ec_lchain_hybrid(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, bw_h, bw_l,
+        max_n_chain_a, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, R_INF.tqn, hf_only, 1);///ONT high error
+
+    // fprintf(stderr, "-b-[M::%s] rid::%ld\n", __func__, i);
+
+    // b->num_read_base += b->olist.length;
+    b->cnt[0] += b->self_read.length;
+
+    aux_o = fetch_aux_ovlp(&b->olist, NULL/**&rse_o**/);///must be here
+
+    copy_asg_arr(buf0, b->sp); 
+    init_gen_hc_aln_t(&ez, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, rse_o, &b->v8q,
+        e_l, e_h, WINDOW_HC, WINDOW_OHC, i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), 0, 1, -1, 0.006, -1, 64, 
+        0, 0, -1, -1, NULL, &(b->v32), &buf0, b->ab, max_n_chain_a, ((max_n_chain_a*HC_MF_R)>0)?(max_n_chain_a*HC_MF_R):1, 
+        asm_opt.chn_occ, ((hom_a*HC_AV_MIN)>0)?(hom_a*HC_AV_MIN):(1), qw, R_INF.tqn, hom_a); 
+        // gen_hc_r_alin_ea_adv(&ez);
+    gen_hc_r_alin_ea_adv_flt(&ez);
+    copy_asg_arr(b->sp, buf0);
+
+
+    copy_asg_arr(buf0, b->sp); 
+    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32, 
+    asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, hf_rate);
+    copy_asg_arr(b->sp, buf0);
+    ///for debug indel
+    // if(i == 23863) stderr_phase_ovlp(&b->olist);
+
+    dedup_chains(&b->olist);
+
+    copy_asg_arr(buf0, b->sp); copy_asg_arr(buf1, b->hap.snp_srt);
+    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, R_INF.tqn, 
+        R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, hf_only?(NULL):(&buf1));
+    copy_asg_arr(b->sp, buf0); copy_asg_arr(b->hap.snp_srt, buf1);
+
+    push_nec_re(aux_o, &(scc.a[i]));
+    push_nec_re(aux_o, &(scb.a[i]));
+
+    // if((asm_opt.is_ont) && is_chemical_r_qual(&b->olist, &b->v64, qlen, 1, 16, &(b->v8q), i)/**(is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))**/) {
+    //     // b->olist.length = 0;
+    //     fprintf(stderr, "[M::%s] rid::%ld\t%.*s\n\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
+    // }
+    push_ne_ovlp_flt(&(R_INF.paf[i]), &b->olist, 1, &R_INF, b->self_read.length, &(scc.a[i]), R_INF.tqn, hf_only);
+    push_ne_ovlp_flt(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, b->self_read.length, NULL, R_INF.tqn, hf_only);
+
+
+    check_well_cal(&(scc.a[i]), &b->v64, &(R_INF.paf[i].is_fully_corrected), &(R_INF.paf[i].is_abnormal), qlen, (MIN_COVERAGE_THRESHOLD*2), &(R_INF.paf[i]));
+    R_INF.trio_flag[i] = AMBIGU;
+
+    // prt_chain(&b->olist);
+
+    // ul_map_lchain(b->abl, (uint32_t)-1, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->clist, s->opt->bw_thres, 
+    //         s->opt->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2/**0.75**/, 2, 3);
+    // exit(1);
+    refresh_ec_ovec_buf_t0(b, REFRESH_N);
+    // exit(0);
+}
+
 
 
 static void worker_hap_ec_dbg_paf(void *data, long i, int tid)
@@ -4480,8 +4923,9 @@ uint32_t quick_exact_match(ma_hit_t *z, All_reads *rref, UC_Read* qu, UC_Read* t
 {
     uint64_t rts, rte, rqs, rqe, f = 0; int64_t ql, tl, qr, tr, qs, qe, ts, te;
 
-    // fprintf(stderr, "-0-[M::%s]\tf::%lu\n", __func__, f);
-    if(adjust_exact_match(&(sc->a[z->tn]), z->ts, z->te, ((uint32_t)(z->qns)), z->qe, &rts, &rte, &rqs, &rqe, z->rev)) {
+    if((R_INF.is_syn == 1) && ((z->qns>>32) >= R_INF.tqn) && (z->tn < R_INF.tqn)) {///HiFi-to-ONT
+        f = 1;
+    } else if(adjust_exact_match(&(sc->a[z->tn]), z->ts, z->te, ((uint32_t)(z->qns)), z->qe, &rts, &rte, &rqs, &rqe, z->rev)) {
         z->ts = rts; z->te = rte; f = 1;
         z->qns >>= 32; z->qns <<= 32; z->qns |= ((uint64_t)(rqs)); z->qe = rqe;
 
@@ -4789,6 +5233,13 @@ static void worker_update_dc_ec(void *data, long i, int tid)
     recover_UC_Read(&b->self_read, &R_INF, i);
     for (k = 0; k < R_INF.paf[i].length; k++) {
         z = &(R_INF.paf[i].buffer[k]);
+        // if(((z->qns>>32) == 23863 && z->tn == 4) || ((z->qns>>32) == 4 && z->tn == 23863)) {
+        //     fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%lu\tq::[%u,%u)\t%c\t%.*s(tid::%u)\t\ttl::%lu\tt::[%u,%u)\tel::%u\n", __func__, 
+        //         (int32_t)Get_NAME_LENGTH(R_INF, (z->qns>>32)), Get_NAME(R_INF, (z->qns>>32)), (uint32_t)(z->qns>>32), Get_READ_LENGTH(R_INF, (z->qns>>32)), 
+        //         (uint32_t)z->qns, z->qe, "+-"[z->rev], 
+        //         (int32_t)Get_NAME_LENGTH(R_INF, z->tn), Get_NAME(R_INF, z->tn), z->tn, Get_READ_LENGTH(R_INF, (z->tn)),  
+        //         z->ts, z->te, z->el?1:0);
+        // }
         if((z->el) && (quick_exact_match(z, &R_INF, &b->self_read, &b->ovlp_read, &scc))) {
             z->el = 1; b->cnt[0]++;
         } else {
@@ -7170,7 +7621,8 @@ static void worker_hap_dc_ec0(void *data, long i, int tid)
     ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
     uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
     uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
-    asg64_v buf0; overlap_region *aux_o = NULL; uint32_t qlen = 0;
+    asg64_v buf0; overlap_region *aux_o = NULL; uint32_t qlen = 0; int64_t het_a, hom_a;
+    set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a);
     // overlap_region *aux_o = NULL; asg64_v buf0;
 
     // gen_ovlst_paf(&(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->v64));
@@ -7203,11 +7655,13 @@ static void worker_hap_dc_ec0(void *data, long i, int tid)
     b->cnt[0] += b->self_read.length;
 
     copy_asg_arr(buf0, b->sp); 
-    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 1**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32);
+    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 1**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, 
+    asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1);
     copy_asg_arr(b->sp, buf0); 
 
     copy_asg_arr(buf0, b->sp);
-    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), NULL);
+    b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), 
+    R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL);
     copy_asg_arr(b->sp, buf0);
 
     push_nec_re(aux_o, &(scc.a[i]));
@@ -7603,6 +8057,215 @@ void cal_ec_multiple_stat_cmp(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, c
     exit(1);
 }
 
+static void ff_ihyb_syn_worker_insert(void *data, long i, int tid) /** callback for kt_for()**/
+{
+    tsrt_v_buf *s = ((tsrt_v_buf*)data);
+    asg64_v *za = &(s->idx[i]); int64_t n0;
+    uint64_t k, qn, tn, ok, ss; ma_hit_t *fv = NULL, *rv = NULL; ma_hit_t_alloc *rva = NULL; 
+    for (k = 0; k < za->n; k++) {
+        qn = za->a[k]>>32; ok = ((uint32_t)za->a[k])>>1; ss = za->a[k]&1;
+        fv = &(s->pf[qn].buffer[ok]); tn = fv->tn;
+        assert((!ss) || (fv->bl == 0x7FFFFFFF));
+        // fv->bl = Get_READ_LENGTH(R_INF, tn);
+        rva = &(s->pf[tn]);
+        if(rva->length >= rva->size) {
+            rva->length++;
+        } else {
+            rv = &(rva->buffer[rva->length++]);
+            rv->qns = (uint64_t)-1;
+        }
+    }	
+
+
+    for (k = 0; k < za->n; k++) {
+        qn = za->a[k]>>32; ok = ((uint32_t)za->a[k])>>1; ss = za->a[k]&1;
+        fv = &(s->pf[qn].buffer[ok]); tn = fv->tn;
+        assert((!ss) || (fv->bl == 0x7FFFFFFF));
+        fv->bl = Get_READ_LENGTH(R_INF, tn);
+        rva = &(s->pf[tn]);
+        n0 = MIN(rva->length, rva->size);
+        if(rva->length > rva->size) {
+            rva->size = rva->length;
+            REALLOC(rva->buffer, rva->size);
+        }
+
+        for (n0--; (n0 >= 0) && (rva->buffer[n0].qns == ((uint64_t)-1)); n0--);
+        rva->length = n0 + 1;
+        rv = &(rva->buffer[rva->length++]);
+
+        rv->qns = Get_tn(*fv); rv->qns = rv->qns << 32; rv->qns = rv->qns | Get_ts(*fv);
+        rv->qe = Get_te(*fv);
+
+        rv->tn = Get_qn(*fv);
+        rv->ts = Get_qs(*fv);
+        rv->te = Get_qe(*fv);
+
+        rv->rev = fv->rev;
+        rv->el = fv->el;
+
+        rv->ml = fv->ml;
+        rv->no_l_indel = fv->no_l_indel;
+        rv->bl = ((ss)?(0x7FFFFFFF):(Get_READ_LENGTH(R_INF, rv->tn)));
+        /**
+        if(rv->el) {///HiFi->ONT: this should extend
+            ql = Get_READ_LENGTH(R_INF, Get_qn(*rv));
+            tl = Get_READ_LENGTH(R_INF, Get_tn(*rv));
+            qs = ((uint32_t)(rv->qns)); qe = rv->qe; ts = rv->ts; te = rv->te;
+            if(qs >= ql) {qs = ql;} if(qe > ql) {qe = ql;} 
+            if(ts >= tl) {ts = tl;} if(te > tl) {te = tl;} 
+
+            if(qs <= ts) {
+                ts -= qs; qs = 0;
+            } else {
+                qs -= ts; ts = 0;
+            }
+
+            qr = ql - qe; tr = tl - te;
+            if(qr <= tr) {
+                qe = ql; te += qr;        
+            } else {
+                te = tl; qe += tr; 
+            }
+
+            rv->qns >>= 32; rv->qns <<= 32; rv->qns |= ((uint64_t)(qs)); rv->qe = qe;
+            rv->ts = ts; rv->te = te;
+            if((qe <= qs) || (te <= ts)) rva->length--;
+        }
+        **/
+        // if(((rv->qns>>32) == 23863 && rv->tn == 4) || ((rv->qns>>32) == 4 && rv->tn == 23863)) {
+        //     fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%lu\tq::[%u,%u)\t%c\t%.*s(tid::%u)\t\ttl::%lu\tt::[%u,%u)\tel::%u\n", __func__, 
+        //     (int32_t)Get_NAME_LENGTH(R_INF, (rv->qns>>32)), Get_NAME(R_INF, (rv->qns>>32)), (uint32_t)(rv->qns>>32), Get_READ_LENGTH(R_INF, (rv->qns>>32)), 
+        //     (uint32_t)rv->qns, rv->qe, "+-"[rv->rev], 
+        //     (int32_t)Get_NAME_LENGTH(R_INF, rv->tn), Get_NAME(R_INF, rv->tn), rv->tn, Get_READ_LENGTH(R_INF, (rv->tn)),  
+        //     rv->ts, rv->te, rv->el?1:0);
+        // }
+    }
+}
+
+static void *ff_ihyb_syn_worker_count(void *data, int step, void *in) 
+{
+    tsrt_v_m *p = (tsrt_v_m*)data;
+    if (step == 0) {
+        ma_hit_t_alloc *z = NULL; uint64_t k, qn, tn, m; asg64_v *zi = NULL;
+        tsrt_v_buf *s = NULL; CALLOC(s, 1);
+        s->p = p->p; s->pm = p->pm; s->pn = p->pn; s->tqn = p->tqn; s->pf = p->pf; CALLOC(s->idx, s->pn);
+
+        while (p->rid < p->tqn) {
+            // fprintf(stderr, "\n[M::%s] p->rid::%lu, p->tot::%lu, R_INF.tqn::%lu\n", __func__, p->rid, p->tot, R_INF.tqn);
+            z = &(p->pf[p->rid++]);
+            // fprintf(stderr, "[M::%s] R_INF.paf[23856].length::%u\n", __func__, R_INF.paf[23856].length);
+            // fprintf(stderr, "[M::%s] s->pf[23856].lengthlength::%u\n", __func__, s->pf[23856].length);
+            // fprintf(stderr, "[M::%s] z->length::%u\n", __func__, z->length);
+            for (k = 0; k < z->length; k++) {
+                tn = z->buffer[k].tn;
+                if(tn < p->tqn) continue;///no ont-2-ont
+                qn = z->buffer[k].qns>>32;
+                m = qn<<=32; m |= (k<<1);
+                if(z->buffer[k].bl == 0x7FFFFFFF) {
+                    m |= 1; p->n_bl++;
+                    // z->buffer[k].bl = Get_READ_LENGTH(R_INF, tn);
+                }
+                zi = &(s->idx[(tn-p->tqn)&s->pm]);
+                s->tov_size -= zi->m;
+                kv_push(uint64_t, *zi, m); s->tov++;
+                s->tov_size += zi->m;
+                p->n_ov++;
+                // if(((z->buffer[k].qns>>32) == 23863 && z->buffer[k].tn == 4) || ((z->buffer[k].qns>>32) == 4 && z->buffer[k].tn == 23863)) {
+                //     fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%lu\tq::[%u,%u)\t%c\t%.*s(tid::%u)\t\ttl::%lu\tt::[%u,%u)\tel::%u\n", __func__,
+                //         (int32_t)Get_NAME_LENGTH(R_INF, (z->buffer[k].qns>>32)), Get_NAME(R_INF, (z->buffer[k].qns>>32)), (uint32_t)(z->buffer[k].qns>>32), Get_READ_LENGTH(R_INF, (z->buffer[k].qns>>32)), 
+                //         (uint32_t)z->buffer[k].qns, z->buffer[k].qe, "+-"[z->buffer[k].rev], 
+                //         (int32_t)Get_NAME_LENGTH(R_INF, z->buffer[k].tn), Get_NAME(R_INF, z->buffer[k].tn), z->buffer[k].tn, Get_READ_LENGTH(R_INF, (z->buffer[k].tn)),  
+                //         z->buffer[k].ts, z->buffer[k].te, z->buffer[k].el?1:0);
+                // }
+            }
+            if(s->tov_size >= p->chunk_size) break;
+        }
+        if (s->tov == 0) {
+            for (k = 0; k < s->pn; k++) {
+                free(s->idx[k].a);
+            }
+            free(s->idx); free(s);
+        } else {
+            return s;
+        }
+    } else if (step == 1) {
+        tsrt_v_buf *s  = (tsrt_v_buf* )in; uint64_t k;
+        kt_for(p->n_thr, ff_ihyb_syn_worker_insert, s, s->pn);
+        for (k = 0; k < s->pn; k++) {
+            free(s->idx[k].a);
+        }
+        free(s->idx); free(s);
+    }
+    
+    return 0;
+}
+
+void ff_ihyb_syn_tid(ec_ovec_buf_t *b, uint64_t pre, uint64_t n_a, uint64_t n_thre)
+{
+    tsrt_v_m sp = {0, 0, 0, 0, 0, 0};
+    sp.p = pre; sp.pn = ((uint64_t)1) << pre; sp.pm = (((uint64_t)1) << pre) - 1; sp.n_ov = sp.n_bl = 0;
+    sp.rid = 0; sp.tot = n_a; sp.chunk_size = 10000000; sp.tqn = R_INF.tqn; sp.n_thr = n_thre;
+
+    sp.pf = R_INF.paf; sp.n_ov = sp.n_bl = 0; sp.rid = 0; 
+    kt_pipeline(n_thre, ff_ihyb_syn_worker_count, &sp, 2);
+    fprintf(stderr, "[M::%s::cis-paf] # syn overlaps::%lu, # syn informative overlaps::%lu\n", __func__, sp.n_ov, sp.n_bl);
+
+
+    sp.pf = R_INF.reverse_paf; sp.n_ov = sp.n_bl = 0; sp.rid = 0; 
+    kt_pipeline(n_thre, ff_ihyb_syn_worker_count, &sp, 2);
+    fprintf(stderr, "[M::%s::trans-paf] # syn overlaps::%lu, # syn informative overlaps::%lu\n", __func__, sp.n_ov, sp.n_bl);
+
+    kt_for(n_thre, worker_hap_ec_hybrid_sync, b, n_a - R_INF.tqn);///HiFi-only
+
+    
+
+
+
+    /**
+    tsrt_v_buf bsrt_v = {0, 0, 0, NULL};
+    bsrt_v.p = pre; bsrt_v.pn = ((uint64_t)1) << pre; bsrt_v.pm = (((uint64_t)1) << pre) - 1;
+    CALLOC(bsrt_v.idx, bsrt_v.pn);
+
+    kt_pipeline(3, ff_ihyb_syn_worker_count, &bsrt_v, 2);
+
+    for (k = 0; k < R_INF.tqn; k++) {
+        z = &(R_INF.paf[k]);
+        for (i = 0; i < z->length; i++) {
+            tn = z->buffer[i].tn;
+            qn = z->buffer[i].qns>>32;
+            m = qn<<=32; m |= (i<<1);
+            if(z->buffer[i].bl == 0x7FFFFFFF) {
+                m |= 1; 
+                z->buffer[i].bl = Get_READ_LENGTH(R_INF, tn);
+            }
+            zi = &(bsrt_v.idx[(tn-R_INF.tqn)&bsrt_v.pm]);
+            kv_push(uint64_t, *zi, m);
+        }
+    }
+    
+    for (k = 0; k < bsrt_v.pn; k++) {
+        zi = &(bsrt_v.idx[k]);
+        radix_sort_ec64(zi->a, zi->a + zi->n);
+    }
+    
+
+
+
+    for (k = 0; k < bsrt_v.pn; k++) free(bsrt_v.idx[k].a);
+    free(bsrt_v.idx); bsrt_v = {0, 0, 0, NULL};
+    **/
+}
+
+void gen_ihyb_syn(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a)
+{
+    // R_INF.is_syn = 1;
+    kt_for(n_thre, worker_hap_ec_hybrid, b, n_a);
+    // R_INF.is_syn = 1;
+
+    ff_ihyb_syn_tid(b, 10, n_a, n_thre);
+    
+}
+
 uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64_t *r_base)
 {
     double tt0 = yak_realtime_0();
@@ -7621,8 +8284,16 @@ uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64
     // fprintf(stderr, "[M::%s] n_thre->%lu\n", __func__, n_thre);
     if(asm_opt.dbg_run_1 && asm_opt.dbg_run_2) cal_ec_multiple_stat_cmp(b, n_thre, n_a, asm_opt.dbg_run_1, asm_opt.dbg_run_2, 0.1, 0.1);
 
-    if(!(asm_opt.hf)) kt_for(n_thre, worker_hap_ec, b, n_a);///debug_for_fix
-    else kt_for(n_thre, worker_hap_ec_hybrid, b, n_a);///debug_for_fix
+    R_INF.is_syn = 0;
+    if(!(asm_opt.hf)) {
+        kt_for(n_thre, worker_hap_ec, b, n_a);///debug_for_fix
+        // exit(1);
+    } else if(asm_opt.hyb_syn == 1) {///all-to-all 
+        kt_for(n_thre, worker_hap_ec_hybrid, b, n_a);///debug_for_fix
+    } else {
+        R_INF.is_syn = 1;
+        gen_ihyb_syn(b, n_thre, n_a);
+    }
 
     for (k = 0; k < n_thre; ++k) {
         num_base += b->a[k].cnt[0];
@@ -7817,6 +8488,20 @@ dbg_cnt_ss* gen_dbg_cnt_ss(uint64_t n_a)
     return p;
 }
 
+void prt_nel_ovlp(ma_hit_t_alloc *pa, uint64_t p_n)
+{
+    uint64_t k, t; ma_hit_t *z;
+    for (k = 0; k < p_n; k++) {
+        for (t = 0; t < pa[k].length; t++) {
+            z = &(pa[k].buffer[t]);
+            fprintf(stderr, "%.*s(qid::%u)\t%c\t%.*s(tid::%u)\tel::%u\n", 
+                (int32_t)Get_NAME_LENGTH(R_INF, (z->qns>>32)), Get_NAME(R_INF, (z->qns>>32)), (uint32_t)(z->qns>>32), "+-"[z->rev], (int32_t)Get_NAME_LENGTH(R_INF, z->tn), Get_NAME(R_INF, z->tn), z->tn,
+            z->el?1:0);
+        }
+    }
+    exit(1);
+}
+
 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)
 {
     // write_ec_reads("ec0.fa");
@@ -7850,6 +8535,9 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u
     // if(is_sv) kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps
     fprintf(stderr, "-2-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]);
 
+    prt_nel_ovlp(R_INF.paf, n_a);
+    exit(1);
+
     if((!is_sv) || (is_sv && is_cr)) {
         kt_for(n_thre, worker_hap_post_rev, b, n_a);
     }