diff --git a/Assembly.cpp b/Assembly.cpp
index f24d9c4..7798f24 100644
--- a/Assembly.cpp
+++ b/Assembly.cpp
@@ -23,6 +23,7 @@ All_reads R_INF;
 Debug_reads R_INF_FLAG;
 all_ul_t UL_INF, ULG_INF;
 uint32_t *het_cnt = NULL;
+// uint32_t debug_out = 0;
 
 void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length)
 {
diff --git a/CommandLines.cpp b/CommandLines.cpp
index f002fe7..761a18b 100644
--- a/CommandLines.cpp
+++ b/CommandLines.cpp
@@ -56,6 +56,7 @@ static ko_longopt_t long_options[] = {
     { "bin-only",     ko_no_argument, 341},
     { "ul-round",     ko_required_argument, 342},
     { "prt-raw",     ko_no_argument, 343},
+    { "integer-correct",     ko_required_argument, 344},
 	{ 0, 0, 0 }
 };
 
@@ -267,6 +268,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->bin_only = 0;
     asm_opt->ul_clean_round = 1;
     asm_opt->prt_dbg_gfa = 0;
+    asm_opt->integer_correct_round = 0;
 }
 
 void destory_enzyme(enzyme* f)
@@ -801,6 +803,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 341) asm_opt->bin_only = 1;
         else if (c == 342) asm_opt->ul_clean_round = atol(opt.arg);
         else if (c == 343) asm_opt->prt_dbg_gfa = 1;
+        else if (c == 344) asm_opt->integer_correct_round = atol(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 afedfde..e737492 100644
--- a/CommandLines.h
+++ b/CommandLines.h
@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.19.0-r534"
+#define HA_VERSION "0.19.0-r554"
 
 #define VERBOSE 0
 
@@ -140,6 +140,7 @@ typedef struct {
     uint8_t bin_only;
     int32_t ul_clean_round;
     int32_t prt_dbg_gfa;
+    int32_t integer_correct_round;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;
diff --git a/Overlaps.cpp b/Overlaps.cpp
index ea320dc..92cf3bf 100644
--- a/Overlaps.cpp
+++ b/Overlaps.cpp
@@ -2770,6 +2770,7 @@ R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ)
     }
 
     asg_cleanup(g);
+    asg_symm(g);
     g->r_seq = g->n_seq;
 	return g;
 }
@@ -10308,16 +10309,14 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
     long long i, j, index;
     uint32_t qn, tn, ou;
 
-    for (i = 0; i < num_sources; i++)
-    {
-        for (j = 0; j < sources[i].length; j++)
-        {
+    for (i = 0; i < num_sources; i++) {
+        for (j = 0; j < sources[i].length; j++) {
             qn = Get_qn(sources[i].buffer[j]);
             tn = Get_tn(sources[i].buffer[j]);
 
             if(sources[i].buffer[j].del) continue;
             ou = (sources[i].buffer[j].bl&((uint32_t)0x3fffffff));
-            //if this is a weak overlap
+            //if this is a weak overlap; ml == 0 -> weak overlap
             if((sources[i].buffer[j].ml == 0) && ((ou_thres==((uint32_t)-1)) || (ou < ou_thres)))
             {   
                 if(
@@ -33680,11 +33679,19 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t,
         if(asm_opt.prt_dbg_gfa) prt_dbg_gfa(sg, "raw", *cov, src, ruIndex, max_hang_length, mini_overlap_length);        
         asg_arc_del_trans(sg, gap_fuzz);
     } else {
+        ug_opt_t uopt; 
         sg = ma_sg_gen_ul(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length, UL_COV_THRES);
         if(asm_opt.prt_dbg_gfa) prt_dbg_gfa(sg, "raw", *cov, src, ruIndex, max_hang_length, mini_overlap_length);
+        gen_ug_opt_t(&uopt, src, NULL, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, *cov, ruIndex, -1, -1, -1, -1, -1, b_mask_t);
+        
+        ///debug
+        // asg_symm(sg);
+        // dedup_contain_g(&uopt, sg);
+        
+        if(clean_contain_g(&uopt, sg, 0)) update_sg_uo(sg, src);
         // prt_specfic_sge(sg, 10498, 10505, "--*--");
         asg_arc_del_trans_ul(sg, gap_fuzz);
-        // prt_specfic_sge(sg, 10498, 10505, "--#--");
+        // prt_specfic_sge(sg, 10498, 10505, "--#--");        
     }
     
     init_bub_label_t(b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq);
@@ -33727,7 +33734,7 @@ bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, con
     int32_t k, strl = strlen(bin_file)+1, kt, cl, sl; char *id = NULL;
     renew_g(src, rev_src, n_read, readLen, cov, ruIndex, sg, mini_overlap_length, max_hang_length,
         uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, b_mask_t, 
-        is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, ((is_trio)?(0):(1))/**1**/);
+        is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, /**((is_trio)?(0):(1))**/1);
     gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen, 
             *cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t);
     ug_ext_gfa(uopt, *sg, ug_ext_len);
diff --git a/Overlaps.h b/Overlaps.h
index f147f3d..5d0107f 100644
--- a/Overlaps.h
+++ b/Overlaps.h
@@ -1206,6 +1206,7 @@ void ma_hit_contained_advance(ma_hit_t_alloc* sources, long long n_read, ma_sub_
 R_to_U* ruIndex, int max_hang, int min_ovlp);
 void hic_clean_adv(asg_t *sg, ug_opt_t *uopt);
 void update_ug_ou(ma_ug_t *ug, asg_t *sg);
+int asg_arc_del_trans_ul(asg_t *g, int fuzz);
 
 #define JUNK_COV 5
 #define DISCARD_RATE 0.8
diff --git a/Process_Read.h b/Process_Read.h
index 88362af..3d0fa09 100644
--- a/Process_Read.h
+++ b/Process_Read.h
@@ -206,6 +206,7 @@ typedef struct
 extern all_ul_t UL_INF;
 extern all_ul_t ULG_INF;
 // extern uint32_t *het_cnt;
+// extern uint32_t debug_out;
 
 void init_All_reads(All_reads* r);
 void malloc_All_reads(All_reads* r);
diff --git a/gfa_ut.cpp b/gfa_ut.cpp
index e109fc3..f20c097 100644
--- a/gfa_ut.cpp
+++ b/gfa_ut.cpp
@@ -22,7 +22,6 @@ KRADIX_SORT_INIT(srt64, uint64_t, generic_key, 8)
 #define ASG_ET_MULTI_NEI 3
 #define UL_TRAV_HERATE 0.2
 #define UL_TRAV_FT_RATE 0.8
-#define is_contain_r(ri, z) (((z)<(ri).len)&&((ri).index[(z)]!=(uint32_t)(-1))&&(!((ri).index[(z)]>>31)))
 
 KDQ_INIT(uint64_t)
 
@@ -98,6 +97,8 @@ typedef struct{
     double min_ovlp_drop_ratio;
     double max_ovlp_drop_ratio;
     double hom_check_drop_rate;
+    double min_path_drop_ratio;
+    double max_path_drop_ratio;
     int64_t max_tip, max_tip_hifi;
     uint32_t is_trio;
 }ulg_opt_t;
@@ -114,6 +115,7 @@ typedef struct {
     size_t n, m;
     uint32_t id:31, is_del:1;
     uint32_t cn;
+    // uint8_t is_consist;
 } ul2ul_item_t;
 
 typedef struct {
@@ -1921,6 +1923,8 @@ void flex_asg_t_cleanup(flex_asg_t *fg)
         fg->g->is_srt = 0;
     }
     asg_cleanup(fg->g);
+    // asg_symm(fg->g);
+    // asg_arc_del_trans_ul(fg->g, fg->gap_fuzz);
 }
 
 void asg_arc_cut_contain(flex_asg_t *fg, asg64_v *in, asg64_v *in0, R_to_U* rI, float ou_rat, uint32_t only_trans_nn)
@@ -2747,7 +2751,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
     double step = 
         (clean_round==1?max_ovlp_drop_ratio:((max_ovlp_drop_ratio-min_ovlp_drop_ratio)/(clean_round-1)));
     double drop = min_ovlp_drop_ratio;
-    int64_t i; asg64_v bu = {0,0,0}, ba = {0,0,0}; uint32_t l_drop = 2000; flex_asg_t *fg = NULL; uint32_t min_diff = (is_ou?2000:0);
+    int64_t i; asg64_v bu = {0,0,0}, ba = {0,0,0}; uint32_t l_drop = 2000; flex_asg_t *fg = NULL; uint32_t min_diff = 0, step_diff = 2000;
     if(is_ou) fg = init_flex_asg_t(sg, uopt->sources, uopt->min_ovlp, uopt->max_hang, asm_opt.max_hang_rate, gap_fuzz);
     // if(is_ou) update_sg_uo(sg, src);///do not do it here
     // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
@@ -2769,9 +2773,13 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
 
 	asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
     // fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847));
-    
+    // if(is_ou) dedup_contain_g(uopt, sg);
     for (i = 0; i < clean_round; i++, drop += step) {
         if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio;
+        if(is_ou) {
+            if(drop <= 0.500001) min_diff = step_diff>>1;
+            else min_diff = step_diff;
+        }
         // fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop);
         // prt_specfic_sge(sg, 10531, 10519, "--0--");
 
@@ -2795,7 +2803,10 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
 
         // prt_specfic_sge(sg, 10531, 10519, "--3--");
         // if(is_ou) asg_arc_cut_contain(fg, &bu, &ba, rI, ((i+1)<clean_round)?ou_drop_rate:-1);
-        if(is_ou) asg_arc_cut_contain(fg, &bu, &ba, rI, ou_drop_rate, 0);
+        if(is_ou) {
+            asg_arc_cut_contain(fg, &bu, &ba, rI, ou_drop_rate, 0);
+            // dedup_contain_g(uopt, sg);
+        }
 
         asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1);
         asg_arc_cut_bub_links(sg, &bu, HARD_OL_DROP, HARD_OL_SEC_DROP, HARD_OU_DROP, is_ou, asm_opt.large_pop_bubble_size, rev, rI, max_tip);
@@ -2813,14 +2824,19 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
 
         if(is_ou) {
             if(ul_refine_alignment(uopt, sg)) update_sg_uo(sg, src);
+            if(clean_contain_g(uopt, sg, 1)) update_sg_uo(sg, src);
         }
     }
+
+    if(is_ou) min_diff = step_diff;
     // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty4.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
     // debug_info_of_specfic_node("m64012_190921_234837/111673711/ccs", sg, rI, "end");
     // debug_info_of_specfic_node("m64011_190830_220126/95028102/ccs", sg, rI, "end");
     if(is_ou) {
         asg_arc_cut_contain(fg, &bu, &ba, rI, ou_drop_rate, 0);
         asg_arc_cut_contain(fg, &bu, &ba, rI, -1, 1);
+        // dedup_contain_g(uopt, sg);
+        if(clean_contain_g(uopt, sg, 1)) update_sg_uo(sg, src);
     }
     if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1);
 
@@ -2863,7 +2879,11 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
     
     ug_ext_gfa(uopt, sg, ug_ext_len);
 
+    // if(is_ou) dedup_contain_g(uopt, sg);
+    // exit(1)
+
     output_unitig_graph(sg, uopt->coverage_cut, o_file, src, rI, uopt->max_hang, uopt->min_ovlp);
+    // exit(1);
     // flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
     flat_soma_v(sg, src, rI);
 
@@ -6262,8 +6282,8 @@ void integer_gen_ovlp(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, ul2ul_it
 {
     ul_str_idx_t *str_idx = &(uidx->pstr); ma_ug_t *ug = uidx->l1_ug; ul2ul_t res;
     ul_str_t *str = &(str_idx->str.a[qid]); integer_aln_t *p; ul_chain_t sc;
-    uint64_t k, z, *hid_a, hid_n, sck, b_n, m; uint32_t vk, vz; uc_block_t *xi;
-    o->n = o->cn = 0; o->id = qid; o->is_del = 0;
+    uint64_t k, z, *hid_a, hid_n, sck, b_n, m/**, ovn = 0**/; uint32_t vk, vz; uc_block_t *xi;
+    o->n = o->cn = 0; o->id = qid; o->is_del = 0; ///o->is_consist = 0;
     // if(qid == 95 || qid == 36) print_integer_seq(ug, str_idx->str.a, qid, 1);
     if(str->cn < 2) return;
     for (k = 0, buf->b.n = 0; k < str->cn; k++) {
@@ -6297,6 +6317,14 @@ void integer_gen_ovlp(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, ul2ul_it
     }
 
     radix_sort_integer_aln_t_srt(buf->b.a, buf->b.a + buf->b.n); 
+
+    // b_n = buf->b.n;
+    // for (k = 1, z = ovn = 0; k <= b_n; k++) {
+    //     if(k == b_n || (buf->b.a[z].tn_rev_qk>>33) != (buf->b.a[k].tn_rev_qk>>33)) {
+    //         ovn++; z = k;
+    //     }
+    // }
+
     b_n = buf->b.n; buf->sc.n = 0;
     for (k = 1, z = 0; k <= b_n; k++) {
         if(k == b_n || (buf->b.a[z].tn_rev_qk>>32) != (buf->b.a[k].tn_rev_qk>>32)) {
@@ -6341,7 +6369,8 @@ void integer_gen_ovlp(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, ul2ul_it
     //     // v = ((uint32_t)str->a[str->cn-1]);
     //     // nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
     // }
-    o->cn = o->n;
+    o->cn = o->n; 
+    // if(ovn == o->cn) o->is_consist = 1;
     if(o->n <= 0) return;
 }
 
@@ -6641,6 +6670,32 @@ void gen_integer_normalize(ul_resolve_t *uidx)
     }
 }
 
+uint64_t dd_path_connect(asg_t *g, ul_str_t *str)
+{
+    if(str->cn < 2) return 1;///actually should return 1, doesn't matter
+    uint64_t i, v, w, nv, k; asg_arc_t *av;
+    v = ((uint32_t)str->a[0]);
+    for (i = 1; i < str->cn; i++) {
+        w = ((uint32_t)str->a[i]);
+
+        nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+        for (k = 0; k < nv; k++) {
+            if(av[k].del) continue;
+            if(av[k].v == w) break;
+        }
+        if(k >= nv) return 0;
+
+        nv = asg_arc_n(g, w^1); av = asg_arc_a(g, w^1);
+        for (k = 0; k < nv; k++) {
+            if(av[k].del) continue;
+            if(av[k].v == (v^1)) break;
+        }
+        if(k >= nv) return 0;
+
+        v = w;
+    }
+    return 1;
+}
 
 void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul2ul_idx_t *ul2, integer_t *buf, int64_t min_dp)
 {
@@ -6713,6 +6768,14 @@ void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul
         if(is_left == 0 && is_right == 0) is_del = 1;
         if(is_left && is_right && is_middle) is_del = 1;
     }
+
+    // if(is_del && str->cn > 1 && o->is_consist && dd_path_connect(uidx->l1_ug->g, str)) {
+    //     is_del = 0;
+    // }
+
+    // fprintf(stderr, "[M::%s] qid::%u, o->is_consist::%u, str->cn::%u, is_connect::%lu\n", 
+    //     __func__, qid, o->is_consist, str->cn, dd_path_connect(uidx->l1_ug->g, str));
+
     /**
     if(!is_del) {
         for (k = 0; k < str->cn; k++) {
@@ -7720,6 +7783,7 @@ void ma_integer_ug_print0(const ma_ug_t *ug, ul_resolve_t *uidx, int print_seq,
     for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA
         p = &ug->u.a[i];
         if(p->m == 0) continue;
+        if(ug->g && ug->g->seq[i].del) continue;
         sprintf(name, "%s%.6d%c", prefix, i + 1, "lc"[p->circ]);
         if(is_seq) {
             gen_ul_seq(uidx, i, &t);
@@ -9369,6 +9433,9 @@ uint32_t ulg_arc_cut_supports(ul_resolve_t *uidx, ma_ug_t *ug, int32_t max_ext,
 float len_rat, uint32_t is_trio, uint32_t topo_level, uint32_t skip_hom, uint32_t *max_drop_len, uint32_t collapse_check,
 asg64_v *in, asg64_v *ib)
 {
+    // fprintf(stderr, "\n[M::%s::] max_ext::%d, max_ext_hifi::%d, len_rat::%f, is_trio::%u, topo_level::%u, skip_hom::%u, collapse_check::%u\n", 
+    //         __func__, max_ext, max_ext_hifi, len_rat, is_trio, topo_level, skip_hom, collapse_check);
+
     asg64_v tx = {0,0,0}, tb = {0,0,0}, *b = NULL, *ub = NULL; asg_t *g = ug->g;
     uint32_t v, w, i, k, kv, nv, kw, nw, cnt = 0, n_vtx = g->n_seq<<1, to_del, collapse;
     asg_arc_t *av, *aw, *ve, *we; uint64_t w_q, w_t, pb;
@@ -9398,7 +9465,8 @@ asg64_v *in, asg64_v *ib)
     kt_for(uidx->str_b.n_thread, worker_update_ul_arc_supports, uidx, b->n);///all ul + ug
     uidx->uovl.iug_tra = NULL;
     // fprintf(stderr, "#[M::%s::] Done\n", __func__);
-
+    // fprintf(stderr, "#[M::%s::] collapse_check::%u, len_rat::%f\n", __func__, collapse_check, len_rat);
+    
     radix_sort_srt64(b->a, b->a + b->n);
     for (k = 0; k < b->n; k++) {
         if(g->arc[(uint32_t)b->a[k]].del) continue;
@@ -9434,7 +9502,7 @@ asg64_v *in, asg64_v *ib)
             }
         }
 
-        // if((v>>1) == 409 && (w>>1) == 407) {
+        // if(((v>>1) == 6788 && (w>>1) == 17213) || ((w>>1) == 6788 && (v>>1) == 17213)) {
         //     fprintf(stderr, "#[M::%s::] v>>1::%u, v&1::%u, kv::%u, w>>1::%u, w&1::%u, kw::%u, collapse::%u\n", 
         //     __func__, v>>1, v&1, kv, w>>1, w&1, kw, collapse);
         // }
@@ -9448,6 +9516,10 @@ asg64_v *in, asg64_v *ib)
                 //     fprintf(stderr, "+[M::%s::] v>>1::%u, v&1::%u, kv::%u, w>>1::%u, w&1::%u, kw::%u, w_q::%lu, w_t::%lu\n", 
                 //     __func__, v>>1, v&1, kv, w>>1, w&1, kw, w_q, w_t);
                 // }
+                // if(((v>>1) == 6788 && (w>>1) == 17213) || ((w>>1) == 6788 && (v>>1) == 17213)) {
+                //     fprintf(stderr, "+[M::%s::] v>>1::%u, v&1::%u, w_q::%lu, w>>1::%u, w&1::%u, w_t::%lu\n", 
+                //     __func__, v>>1, v&1, w_q, w>>1, w&1, w_t);
+                // }
                 if(w_q == (uint64_t)-1) continue;
                 if(w_q > w_t*len_rat) continue;
             }
@@ -9460,13 +9532,21 @@ asg64_v *in, asg64_v *ib)
                 //     fprintf(stderr, "-[M::%s::] v>>1::%u, v&1::%u, kv::%u, w>>1::%u, w&1::%u, kw::%u, w_q::%lu, w_t::%lu\n", 
                 //     __func__, v>>1, v&1, kv, w>>1, w&1, kw, w_q, w_t);
                 // }
+                // if(((v>>1) == 6788 && (w>>1) == 17213) || ((w>>1) == 6788 && (v>>1) == 17213)) {
+                //     fprintf(stderr, "-[M::%s::] v>>1::%u, v&1::%u, w_q::%lu, w>>1::%u, w&1::%u, w_t::%lu\n", 
+                //     __func__, v>>1, v&1, w_q, w>>1, w&1, w_t);
+                // }
                 if(w_q == (uint64_t)-1) continue;
                 if(w_q > w_t*len_rat) continue;
             }
         }
 
         to_del = check_ulg_to_del(uidx, ug, v, w, kv, kw, max_ext, max_ext_hifi, topo_level, collapse, b, ub);
-
+        
+        // if(((v>>1) == 6788 && (w>>1) == 17213) || ((w>>1) == 6788 && (v>>1) == 17213)) {
+        //     fprintf(stderr, "#[M::%s::] v>>1::%u, v&1::%u, kv::%u, w>>1::%u, w&1::%u, kw::%u, to_del::%u\n", 
+        //     __func__, v>>1, v&1, kv, w>>1, w&1, kw, to_del);
+        // }
         if (to_del) {
             ve->del = we->del = 1, ++cnt;
         }
@@ -12623,6 +12703,7 @@ typedef struct {
     uint64_t int_an;
     uint64_t *ridx_a; 
     uint64_t *ridx;
+    // uint8_t is_double_check;
 } unique_bridge_check_t;
 
 uint64_t get_arc_support(ul_resolve_t *uidx, uint64_t v, uint64_t w)
@@ -12731,6 +12812,53 @@ uint64_t get_arc_support_chain(ul_resolve_t *uidx, uint64_t *a, uint64_t a_n, us
     return occ;
 }
 
+// uint64_t is_consist_ul(ul_resolve_t *uidx, uint64_t *a, uint64_t a_n, usg_t *ng)
+// {
+//     if(a_n <= 0) return 0;
+//     uint64_t k, v, w, z, *hid_a, hid_n, ulid, i; ul2ul_item_t *it;;
+//     ul_str_idx_t *str_idx = &(uidx->pstr);
+//     for (k = 0; k < a_n; k++) {
+//         v = (ng->a[a[k]>>1].mm<<1)|(a[k]&1);
+//         hid_a = str_idx->occ.a + str_idx->idx.a[v>>1];
+//         hid_n = str_idx->idx.a[(v>>1)+1] - str_idx->idx.a[v>>1];
+//         for (z = 0; z < hid_n; z++) {
+//             ulid = hid_a[z]>>32;
+//             if(ulid < uidx->uovl.uln) {
+//                 it = get_ul_ovlp(&(uidx->uovl), ulid, 1);
+//                 if(!it) continue;
+//                 assert(uidx->pstr.str.a[ulid].cn > 1);
+//                 if(it->is_consist == 0) return 0;
+//             }
+//         }
+//     }
+
+//     uint64_t nv; asg_arc_t *av;
+//     v = (ng->a[a[0]>>1].mm<<1)|(a[0]&1);
+//     for (i = 1; i < a_n; i++) {
+//         w = (ng->a[a[i]>>1].mm<<1)|(a[i]&1);
+
+//         nv = asg_arc_n(uidx->l1_ug->g, v); 
+//         av = asg_arc_a(uidx->l1_ug->g, v);
+//         for (k = 0; k < nv; k++) {
+//             if(av[k].del) continue;
+//             if(av[k].v == w) break;
+//         }
+//         if(k >= nv) return 0;
+
+//         nv = asg_arc_n(uidx->l1_ug->g, w^1); 
+//         av = asg_arc_a(uidx->l1_ug->g, w^1);
+//         for (k = 0; k < nv; k++) {
+//             if(av[k].del) continue;
+//             if(av[k].v == (v^1)) break;
+//         }
+//         if(k >= nv) return 0;
+
+//         v = w;
+//     }
+
+//     return 1;
+// }
+
 static void worker_unique_bridge_check_s(void *data, long i, int tid) // callback for kt_for()
 {
     unique_bridge_check_t *uaux = (unique_bridge_check_t *)data;
@@ -12747,7 +12875,9 @@ static void worker_unique_bridge_check_s(void *data, long i, int tid) // callbac
             // (ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1));
             nse = get_arc_support_chain(uaux->uidx, integer_seq+k-1, 2, ng);
             if(nse < unique_bridge_occ) {
-                gidx[i] |= ((uint64_t)0x8000000000000000); return;
+                // if((!(uaux->is_double_check)) || (!is_consist_ul(uaux->uidx, integer_seq+k-1, 2, ng))) {
+                    gidx[i] |= ((uint64_t)0x8000000000000000); return;
+                // }
             }
         }
         if(IF_HOM((integer_seq[k]>>1), *bub)) continue;
@@ -12757,7 +12887,9 @@ static void worker_unique_bridge_check_s(void *data, long i, int tid) // callbac
             // nse = get_arc_support(uaux->uidx, v, w);
             nse = get_arc_support_chain(uaux->uidx, integer_seq+vk, wk+1-vk, ng);
             if(nse < unique_bridge_occ) {
-                gidx[i] |= ((uint64_t)0x8000000000000000); return;
+                // if((!(uaux->is_double_check)) || (!is_consist_ul(uaux->uidx, integer_seq+vk, wk+1-vk, ng))) {
+                    gidx[i] |= ((uint64_t)0x8000000000000000); return;
+                // }
             }
         }
         vk = wk;
@@ -12772,6 +12904,10 @@ uint32_t ava_pass_unique_bridge_cov(ul_resolve_t *uidx, usg_t *g, asg64_v *b64,
     uaux.uidx = uidx; uaux.integer_seq = integer_seq; 
     uaux.gidx = b64->a + g_s; uaux.gidx_n = g_e - g_s;
     uaux.interval_idx = b64->a; uaux.ng = g;
+    ///if tip_l == 0, it is more likely to be right, so give more chance by double checking
+    // if(!tip_l) uaux.is_double_check = 1;
+    // else uaux.is_double_check = 0;
+
     kt_for(uidx->str_b.n_thread, worker_unique_bridge_check_s, (&uaux), g_e-g_s);///seq->n > 1
 
     for (i = g_s; i < g_e; i++) {///available intervals within the same cluster
@@ -12813,8 +12949,7 @@ uint8_t *ff, uint32_t *ng_occ, uint32_t max_ext)
         if(k == b64->n || (b64->a[k]>>32) != (b64->a[i]>>32)) {
             tip_l = ava_pass_unique_bridge_tips(ng, b64, i, k, int_a, ff, max_ext);
             if(tip_l < max_ext) {///no long tip
-                if(ava_pass_unique_bridge_cov(uidx, ng, b64, i, k, int_a, tip_l)) 
-                {
+                if(/**(!tip_l) || (**/ava_pass_unique_bridge_cov(uidx, ng, b64, i, k, int_a, tip_l)) {
                     for (z = i; z < k; z++) {
                         b64->a[mm++] = b64->a[z];
                     }
@@ -15694,11 +15829,11 @@ void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, ui
         for (z = 0; z < ub->n; z++) {///all unreliable arcs
             i = ub->a[z]; tm = 1; p = get_usg_arc(ng, seq->a[i].v, seq->a[i+1].v);
             if(p) {
-                if(p->ou < (uint64_t)tm) p->ou = tm;
+                if(p->ou < (uint64_t)tm) p->ou = tm;///since the initial ou is 0, p->ou = 1
                 pushp_usg_arc_mm(ng, seq->a[i].v, seq->a[i+1].v, k, i); 
 
                 p = get_usg_arc(ng, seq->a[i+1].v^1, seq->a[i].v^1);
-                if(p->ou < (uint64_t)tm) p->ou = tm;
+                if(p->ou < (uint64_t)tm) p->ou = tm;///since the initial ou is 0, p->ou = 1
                 pushp_usg_arc_mm(ng, seq->a[i+1].v^1, seq->a[i].v^1, k, i); 
             } 
             ///give up unreliable arcs if they are not adjacent
@@ -15779,14 +15914,18 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg, uint
 {
     ul2ul_idx_t *idx = &(uidx->uovl); asg64_v bu = {0,0,0}, uu = {0,0,0}; int64_t max_tip_hifi0 = ulopt->max_tip_hifi;
     ma_ug_t *iug = idx->i_ug; int64_t i, mm_tip = ulopt->max_tip; uint64_t cnt = 1, topo_level, ss = 0;
-    double step = (ulopt->clean_round==1?ulopt->max_ovlp_drop_ratio:
-                        ((ulopt->max_ovlp_drop_ratio-ulopt->min_ovlp_drop_ratio)/(ulopt->clean_round-1)));
-    double drop = ulopt->min_ovlp_drop_ratio; CALLOC(iug->g->seq_vis, iug->g->n_seq*2);
+    double step = (ulopt->clean_round==1?ulopt->max_path_drop_ratio:
+                        ((ulopt->max_path_drop_ratio-ulopt->min_path_drop_ratio)/(ulopt->clean_round-1)));
+    double drop = ulopt->min_path_drop_ratio; CALLOC(iug->g->seq_vis, iug->g->n_seq*2);
+
+
+    // char sb[1000]; 
+    // output_integer_graph(uidx, iug, "ig_h0", 0);
 
 
     for (ss = 0; ss < 2; ss++) {
-        for (i = 0, drop = ulopt->min_ovlp_drop_ratio; i < ulopt->clean_round; i++, drop += step) {
-            if(drop > ulopt->max_ovlp_drop_ratio) drop = ulopt->max_ovlp_drop_ratio;
+        for (i = 0, drop = ulopt->min_path_drop_ratio; i < ulopt->clean_round; i++, drop += step) {
+            if(drop > ulopt->max_path_drop_ratio) drop = ulopt->max_path_drop_ratio;
             // fprintf(stderr, "\n[M::%s::] Starting round-%ld, drop::%f\n", __func__, i, drop);
             cnt = 1; topo_level = 2; mm_tip = ulopt->max_tip;
             while (cnt) {
@@ -15812,16 +15951,23 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg, uint
                 cnt += ulg_arc_cut_tips(uidx, iug, mm_tip, ulopt->max_tip_hifi, 0, &bu, &uu);
             }
             // fprintf(stderr, "[M::%s::] Done round-%ld, drop::%f\n", __func__, i, drop);
+            // sprintf(sb, "ig_ss_%lu_i_%ld", ss, i);
+            // output_integer_graph(uidx, iug, sb, 0);
         }
 
         ulg_pop_bubble(uidx, iug, NULL, ((int64_t)0x7fffffff), ulopt->max_tip_hifi, 1, &bu, &uu);
 
         while(ulg_arc_cut_z(uidx, iug, ((int64_t)0x7fffffff), ulopt->max_tip_hifi, 1.5, 0.15, 100, 0.8, ulopt->is_trio, 1, NULL, &bu, &uu));
 
+        // sprintf(sb, "ig_ss_%lu_i_%ld", ss, i);
+        // output_integer_graph(uidx, iug, sb, 0);
+
         ulopt->max_tip_hifi <<= 1;
     }
     ulopt->max_tip_hifi = max_tip_hifi0;
 
+    // output_integer_graph(uidx, iug, "ig_h1", 0);
+
     u2g_threading(uidx, ulopt, 3, is_bridg, &bu, &uu);
 
 
@@ -16132,7 +16278,12 @@ void renew_u2g_bg(ul_resolve_t *uidx)
     bg->bg = asg_init();
     bg->bg->n_seq = 0; bg->bg->m_seq = raw->g->n_seq; MALLOC(bg->bg->seq, bg->bg->m_seq);
     for (k = 0; k < raw->g->n_seq; k++) {
-        raw_id = k; if(IF_HOM(raw_id, *bub)) continue;
+        raw_id = k; 
+        if(IF_HOM(raw_id, *bub)) {///updated-line
+            // asg_seq_set(bg->bg, k, raw->g->seq[k].len, 1);
+            // bg->bg->seq[k].c = 0;
+            continue;
+        }
         raw_a = idx->cc.iug_b + idx->cc.iug_idx[raw_id];
         raw_n = idx->cc.iug_idx[raw_id+1] - idx->cc.iug_idx[raw_id];
         for (z = 0, buf.n = 0; z < raw_n; z++) {
@@ -16322,7 +16473,7 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin
 
     kt_for(uidx->str_b.n_thread, worker_integert_debug_sym, uidx, z->tot);///all ul + ug
     print_integert_ovlp_stat(z);
-    print_uls_seq(uidx, asm_opt.output_file_name);
+    // print_uls_seq(uidx, asm_opt.output_file_name);
     // print_uls_ovs(uidx, asm_opt.output_file_name);
     
     z->i_g = integer_sg_gen(uidx, uopt->min_ovlp);
@@ -16369,8 +16520,10 @@ void gen_cul_g_t(ul_resolve_t *uidx)
     }
 }
 
-void init_ulg_opt_t(ulg_opt_t *z, ug_opt_t *uopt, int64_t clean_round, double min_ovlp_drop_ratio, 
-double max_ovlp_drop_ratio, double hom_check_drop_rate, int64_t max_tip, int64_t max_tip_hifi, bub_label_t *b_mask_t, uint32_t is_trio)
+void init_ulg_opt_t(ulg_opt_t *z, ug_opt_t *uopt, int64_t clean_round, 
+double min_path_drop_ratio, double max_path_drop_ratio,
+double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double hom_check_drop_rate, 
+int64_t max_tip, int64_t max_tip_hifi, bub_label_t *b_mask_t, uint32_t is_trio)
 {
     z->tipsLen = uopt->tipsLen;
     z->tip_drop_ratio = uopt->tip_drop_ratio;
@@ -16381,6 +16534,8 @@ double max_ovlp_drop_ratio, double hom_check_drop_rate, int64_t max_tip, int64_t
 
     z->b_mask_t = b_mask_t;
     z->clean_round = clean_round;
+    z->min_path_drop_ratio = min_path_drop_ratio;
+    z->max_path_drop_ratio = max_path_drop_ratio;
     z->min_ovlp_drop_ratio = min_ovlp_drop_ratio;
     z->max_ovlp_drop_ratio = max_ovlp_drop_ratio;
     z->hom_check_drop_rate = hom_check_drop_rate;
@@ -17070,7 +17225,7 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, ui
     if(deep_clean) {
         // print_debug_gfa(sg, NULL, uopt->coverage_cut, "bclean", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
         deep_graph_clean(uopt, sg, 1, is_trio, asm_opt.max_short_tip, asm_opt.min_drop_rate, 
-        MIN(asm_opt.max_drop_rate, 0.7), 0.75, 1, asm_opt.clean_round, 20);
+        MIN(asm_opt.max_drop_rate, 0.7), 0.75, 1, asm_opt.clean_round, 10/**20**/);
         // print_debug_gfa(sg, NULL, uopt->coverage_cut, "aclean", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
     }
     
@@ -17079,11 +17234,11 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, ui
     // fprintf(stderr, "2[M::%s]\n", __func__);
     // exit(1);
 
-    char* gfa_name = NULL; MALLOC(gfa_name, strlen(o_file)+strlen(bin_file)+50);
-    sprintf(gfa_name, "%s.%s", o_file, bin_file);
-    print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
-    print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
-    free(gfa_name);
+    // char* gfa_name = NULL; MALLOC(gfa_name, strlen(o_file)+strlen(bin_file)+50);
+    // sprintf(gfa_name, "%s.%s", o_file, bin_file);
+    // print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+    // print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
+    // free(gfa_name);
     
     
     filter_sg_by_ug(sg, init_ug, uopt);
@@ -17102,12 +17257,12 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, ui
     // print_ul_alignment(init_ug, &UL_INF, 47072, "after-2");
     // exit(1);
     // if(free_uld) {
-        print_raw_uls_seq(uidx, asm_opt.output_file_name);
-    //     print_raw_uls_aln(uidx, asm_opt.output_file_name);
+        // print_raw_uls_seq(uidx, asm_opt.output_file_name);
+        // print_raw_uls_aln(uidx, asm_opt.output_file_name);
     // }
     
-    ul_re_correct(uidx, 3); 
-    init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_ul_tip, b_mask_t, is_trio);
+    ul_re_correct(uidx, asm_opt.integer_correct_round/**3**/); 
+    init_ulg_opt_t(&uu, uopt, clean_round, 0.2, 0.6, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_ul_tip, b_mask_t, is_trio);
     // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
     /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0, is_bridg);
     // print_ul_alignment(init_ug, &UL_INF, 47072, "after-3");
@@ -17304,12 +17459,9 @@ uint32_t usg_topocut_aux_unambi1(asg_t *g, uint32_t v)
 int usg_topocut_aux_del(ma_ug_t *ug, uint32_t v, int max_ext, asg64_v *b)
 {
 	int32_t n_ext; 
-	for (n_ext = 1; n_ext < max_ext && v != (uint32_t)-1; ++n_ext) {
-		if (usg_topocut_aux_unambi1(ug->g, v^1) == (uint32_t)-1) {
-			--n_ext;
-			break;
-		}
-        kv_push(uint64_t, *b, v);
+	for (n_ext = 0; n_ext < max_ext && v != (uint32_t)-1; ) {
+		if (usg_topocut_aux_unambi1(ug->g, v^1) == (uint32_t)-1) break;
+        if(b) kv_push(uint64_t, *b, v); n_ext += ug->u.a[v>>1].n;
 		v = usg_topocut_aux_unambi1(ug->g, v);
 	}
 	return n_ext;
@@ -17335,10 +17487,28 @@ inline void asg_arc_del_by_ug(asg_t *sg, ma_ug_t *ug, uint32_t uv, uint32_t uw,
     av[i].del = del;
 }
 
-void cal_bub_best_by_len(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t is_trio, uint32_t is_ou, double len_rat, double ou_rat, uint32_t min_ou)
+uint32_t cal_utg_occ(ma_ug_t *ug, uint32_t begNode)
+{
+    uint32_t v = begNode, w, kv, occ = 0; ma_utg_t *u;
+
+    while (1) {
+        kv = get_real_length(ug->g, v, NULL);
+        u = &(ug->u.a[v>>1]); occ += u->n; 
+        if(kv!=1) break;
+        ///kv must be 1 here
+        kv = get_real_length(ug->g, v, &w);
+        if(get_real_length(ug->g, w^1, NULL)!=1) break;
+        v = w;
+        if(v == begNode) break;
+    }
+    return occ;
+}
+
+void cal_bub_best_by_len(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t is_trio, uint32_t is_ou, double len_rat, double ou_rat, uint32_t min_ou, 
+uint32_t min_node)
 {
     // fprintf(stderr, "[M::%s]\tStart\n", __func__);
-    ma_ug_t *ug = sl->ug; asg_t *g = sl->ug->g; uint32_t ol_max, ou_max;
+    ma_ug_t *ug = sl->ug; asg_t *g = sl->ug->g; uint32_t ol_max, ou_max, lnid;
     uint32_t v, w, n_vtx = (g->n_seq<<1), nv, nw, i, k, z, kv, kw, bb, to_del, bn, tip;
     asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av, *aw, *ve, *we; ma_utg_t *u;
     uint32_t trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, mm_ol, mm_ou, cnt = 0, del_v, del_w;
@@ -17391,7 +17561,7 @@ void cal_bub_best_by_len(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t
             }
         }
         ///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted
-        mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou);
+        mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou); lnid = 0;
 
         for (i = kv = ol_max = ou_max = 0; i < nv; ++i) {
             if(av[i].del) continue;
@@ -17399,6 +17569,9 @@ void cal_bub_best_by_len(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t
             if(is_trio && get_ug_tip_trio_infor(ug, av[i].v) == ntrioF) continue;
             if(ol_max < av[i].ol) ol_max = av[i].ol; 
             if(ou_max < av[i].ou) ou_max = av[i].ou;
+            if((!lnid) && ((is_best_arc((*sl), v, av[i].v)))) {
+                if((cal_utg_occ(ug, v^1) >= min_node) || (cal_utg_occ(ug, av[i].v) >= min_node)) lnid = 1;
+            }
         }
         if (kv < 1) continue;
         if (kv >= 2) {
@@ -17413,6 +17586,9 @@ void cal_bub_best_by_len(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t
             if(is_trio && get_ug_tip_trio_infor(ug, aw[i].v) == ntrioF) continue;
             if(ol_max < aw[i].ol) ol_max = aw[i].ol; 
             if(ou_max < aw[i].ou) ou_max = aw[i].ou;
+            if((!lnid) && ((is_best_arc((*sl), w, aw[i].v)))) {
+                if((cal_utg_occ(ug, w^1) >= min_node) || (cal_utg_occ(ug, aw[i].v) >= min_node)) lnid = 1;
+            }
         }
         if (kw < 1) continue;
         if (kw >= 2) {
@@ -17421,22 +17597,26 @@ void cal_bub_best_by_len(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t
         }
 
         if (kv <= 1 && kw <= 1) continue;
+        if(!lnid) continue;
 
         to_del = 0; del_v = del_w = (uint32_t)-1; tip = 0;
         if (kv > 1 && kw > 1) {
             to_del = 1;
         } else if (kw == 1) {
-            tip = asg_topocut_aux(g, w^1, max_ext);
+            // tip = asg_topocut_aux(g, w^1, max_ext);
+            tip = usg_topocut_aux_del(ug, w^1, max_ext, NULL);
             if (tip < max_ext) {
                 to_del = 1; del_w = w^1;
             }
         } else if (kv == 1) {
-            tip = asg_topocut_aux(g, v^1, max_ext);
+            // tip = asg_topocut_aux(g, v^1, max_ext);
+            tip = usg_topocut_aux_del(ug, v^1, max_ext, NULL);
             if (tip < max_ext) {
                 to_del = 1; del_v = v^1;
             }
         }
 
+        
         if (to_del) {
             bn = b->n;
             if(del_v != ((uint32_t)-1)) usg_topocut_aux_del(ug, del_v, max_ext, b);
@@ -17630,7 +17810,7 @@ double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double ou_rat, int64_t m
         if(sl->len_rat > max_ovlp_drop_ratio) sl->len_rat = max_ovlp_drop_ratio;
         update_ug_clean_t(sl);
         kt_for(asm_opt.thread_num, cal_bub_best, sl, sl->ug->g->n_seq);
-        cal_bub_best_by_len(sl, &b, max_ext, is_trio, sl->is_ou, sl->len_rat, sl->ou_rat, min_ou);
+        cal_bub_best_by_len(sl, &b, max_ext, is_trio, sl->is_ou, sl->len_rat, sl->ou_rat, min_ou, long_tip);
     }
 
     // if(is_ou) {
diff --git a/gfa_ut.h b/gfa_ut.h
index 4640e33..6b47493 100644
--- a/gfa_ut.h
+++ b/gfa_ut.h
@@ -3,6 +3,7 @@
 #include "Overlaps.h"
 #include "hic.h"
 
+#define is_contain_r(ri, z) (((z)<(ri).len)&&((ri).index[(z)]!=(uint32_t)(-1))&&(!((ri).index[(z)]>>31)))
 
 typedef struct {
 	asg_t *g; 
@@ -37,5 +38,7 @@ void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc
 bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het);
 void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt);
 void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len);
+void update_sg_uo(asg_t *g, ma_hit_t_alloc *src);
+uint32_t get_arcs(asg_t *g, uint32_t v, uint32_t* idx, uint32_t idx_n);
 
 #endif
diff --git a/inter.cpp b/inter.cpp
index ff689e0..4b54085 100644
--- a/inter.cpp
+++ b/inter.cpp
@@ -15097,6 +15097,7 @@ void dd_ul_vec_t(const ul_idx_t *uref, mg_lchain_t *a, int64_t a_n, ul_vec_t *rc
 					z->qs = a[k].qs; z->qe = a[k].qe;
 					z->te = a[k].re; z->ts = a[k].rs;
 					z->pidx = k;
+					z->aidx = z->pdis = (uint32_t)-1;
 				} else {
 					z->qs = a[k].qs; z->qe = a[k].qe;
 					z->te = a[k].re; z->ts = a[k].rs;
@@ -15110,6 +15111,7 @@ void dd_ul_vec_t(const ul_idx_t *uref, mg_lchain_t *a, int64_t a_n, ul_vec_t *rc
 			z->qs = a[k].qs; z->qe = a[k].qe;
 			z->te = a[k].re; z->ts = a[k].rs;
 			z->pidx = k;
+			z->aidx = z->pdis = (uint32_t)-1;
 		}
 	}
 
@@ -15191,6 +15193,42 @@ void dd_ul_vec_t(const ul_idx_t *uref, mg_lchain_t *a, int64_t a_n, ul_vec_t *rc
 	if(sp != (uint32_t)-1) l += ep - sp;
 	if(l == (int64_t)rch->rlen) rch->dd = 1;
 
+
+	///debug_ssb
+	// for (k = 0; (uint32_t)k < rch->bb.n; k++) {
+	// 	if(rch->bb.a[k].pidx != (uint32_t)-1) {
+	// 		if((rch->bb.a[k].pidx >= 0) && (rch->bb.a[k].pidx < rch->bb.n) && 
+	// 					(rch->bb.a[rch->bb.a[k].pidx].aidx == k)) {
+	// 			;
+	// 		} else {
+	// 			fprintf(stderr, "[M::%s::k->%ld] +rch->bb.n::%u\n", __func__, k, (uint32_t)rch->bb.n);
+	// 			for (l = 0; (uint32_t)l < rch->bb.n; l++) {
+	// 				fprintf(stderr, "[M::%.*s::k->%lu] q::[%u, %u), t::[%u, %u), pidx::%u, aidx::%u\n", 
+	// 				(int)Get_NAME_LENGTH(R_INF, rch->bb.a[l].hid), Get_NAME(R_INF, rch->bb.a[l].hid), l, 
+	// 				rch->bb.a[l].qs, rch->bb.a[l].qe, rch->bb.a[l].ts, rch->bb.a[l].te, 
+	// 				rch->bb.a[l].pidx, rch->bb.a[l].aidx);
+	// 			}
+	// 			exit(1);
+	// 		}
+	// 	}
+
+	// 	if(rch->bb.a[k].aidx != (uint32_t)-1) {
+	// 		if((rch->bb.a[k].aidx >= 0) && (rch->bb.a[k].aidx < rch->bb.n) && 
+	// 					(rch->bb.a[rch->bb.a[k].aidx].pidx == k)) {
+	// 			;
+	// 		} else {
+	// 			fprintf(stderr, "[M::%s::k->%ld] -rch->bb.n::%u\n", __func__, k, (uint32_t)rch->bb.n);
+	// 			for (l = 0; (uint32_t)l < rch->bb.n; l++) {
+	// 				fprintf(stderr, "[M::%.*s::k->%lu] q::[%u, %u), t::[%u, %u), pidx::%u, aidx::%u\n", 
+	// 				(int)Get_NAME_LENGTH(R_INF, rch->bb.a[l].hid), Get_NAME(R_INF, rch->bb.a[l].hid), l, 
+	// 				rch->bb.a[l].qs, rch->bb.a[l].qe, rch->bb.a[l].ts, rch->bb.a[l].te, 
+	// 				rch->bb.a[l].pidx, rch->bb.a[l].aidx);
+	// 			}
+	// 			exit(1);
+	// 		}
+	// 	}
+	// }
+
 	// if(ulid == 292) {
 	// 	for (i = 0; i < rch->bb.n; i++) {
 	// 		fprintf(stderr, "(%lu) qs:%u, qe:%u, ts:%u, te:%u, pidx:%u\n", i, 
@@ -15428,6 +15466,1215 @@ static void worker_for_ul_gchains_alignment(void *data, long i, int tid)
 	// gl_chain_refine_advance(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, km);
 }
 
+uint32_t inline is_rg_connect(const asg_t *rg, uint32_t v, uint32_t w)
+{
+	uint32_t nv = asg_arc_n(rg, v), i; 
+	asg_arc_t *av = asg_arc_a(rg, v);
+    for (i = 0; i < nv; i++) {
+		if((!(av[i].del)) && (av[i].v == w)) return 1;
+    }
+	return 0;
+}
+
+uint32_t refine_contain_consensus_chain(const asg_t *rg, ul_vec_t *rch, R_to_U *ri, asg64_v *bu, uint64_t ulid)
+{
+	bu->n = 0; 
+	if(rch->bb.n == 1 && rch->bb.a[0].base) return bu->n;///no alignment
+	if(rch->bb.n == 0) return bu->n;///no alignment
+	uint64_t i, m, nc, occ, is_end; int64_t k, kn;
+	kv_resize(uint64_t, *bu, rch->bb.n); 
+	memset(bu->a, -1, (sizeof((*(bu->a)))*rch->bb.n));
+	for (k = rch->bb.n-1; k >= 0; k--) {
+		// if(ulid == 2 || ulid == 46 || ulid == 475 || ulid == 3274 || ulid == 3360) {
+		// 	fprintf(stderr, "[M::%s::] ulid::%lu, k::%ld, q::[%u, %u), t::[%u, %u), del::%u, pidx::%u, pconn::%u, aidx::%u, aconn::%u\n", 
+		// 		__func__, ulid, k, 
+		// 		rch->bb.a[k].qs, rch->bb.a[k].qe, rch->bb.a[k].ts, rch->bb.a[k].te, rg->seq[rch->bb.a[k].hid].del,
+		// 		rch->bb.a[k].pidx, 
+		// 		((rch->bb.a[k].pidx!=((uint32_t)-1))&&is_rg_connect(rg, (uint64_t)((rch->bb.a[rch->bb.a[k].pidx].hid<<1)|(rch->bb.a[rch->bb.a[k].pidx].rev)), (uint64_t)((rch->bb.a[k].hid<<1)|(rch->bb.a[k].rev))))?1:0,
+		// 		rch->bb.a[k].aidx,
+		// 		((rch->bb.a[k].aidx!=((uint32_t)-1))&&is_rg_connect(rg, (uint64_t)((rch->bb.a[k].hid<<1)|(rch->bb.a[k].rev)),(uint64_t)((rch->bb.a[rch->bb.a[k].aidx].hid<<1)|(rch->bb.a[rch->bb.a[k].aidx].rev))))?1:0);
+		// }
+
+		
+
+		if(bu->a[k] != ((uint64_t)-1)) continue;
+		is_end = 0;
+		if((rch->bb.a[k].aidx == ((uint32_t)-1)) && (rch->bb.a[k].pidx != ((uint32_t)-1))) {///the end of a chain
+			is_end = 1;
+		} else if((rch->bb.a[k].aidx != ((uint32_t)-1)) && (rch->bb.a[k].pidx != ((uint32_t)-1))) {
+			if(!is_rg_connect(rg, (uint64_t)((rch->bb.a[k].hid<<1)|(rch->bb.a[k].rev)), 
+						(uint64_t)((rch->bb.a[rch->bb.a[k].aidx].hid<<1)|(rch->bb.a[rch->bb.a[k].aidx].rev)))) {
+				is_end = 1;
+			}
+		}
+		if(!is_end) continue;
+		i = m = k; occ = 0;
+		while (i != ((uint32_t)-1)) {
+			assert(bu->a[i] == ((uint64_t)-1));
+			if(rg->seq[rch->bb.a[i].hid].del) {
+				if(occ > 1) {
+					bu->a[m] = occ;
+					// if(debug_out) {
+					// 	fprintf(stderr, "+[M::%s::] ulid::%lu, k::%ld, m::%lu, occ::%lu\n", 
+					// 	__func__, ulid, k, m, occ);
+					// }
+				}
+				// assert(i!=m);
+				bu->a[i] = 0; m = (uint32_t)-1; occ = 0;
+			} else {
+				if(m == (uint32_t)-1) m = i;
+				bu->a[i] = 0; occ++;
+			}
+			i = rch->bb.a[i].pidx;
+			if(i != ((uint32_t)-1)) {
+				if(!is_rg_connect(rg, (uint64_t)((rch->bb.a[i].hid<<1)|(rch->bb.a[i].rev)), 
+						(uint64_t)((rch->bb.a[rch->bb.a[i].aidx].hid<<1)|(rch->bb.a[rch->bb.a[i].aidx].rev)))) {
+					i = ((uint32_t)-1);
+				}
+			}
+		}
+
+		if(occ > 1) {
+			bu->a[m] = occ;
+			// if(debug_out) {
+			// 	fprintf(stderr, "+[M::%s::] ulid::%lu, k::%ld, m::%lu, occ::%lu\n", 
+			// 	__func__, ulid, k, m, occ);
+			// }
+		}
+	}
+
+	kn = rch->bb.n;
+	for (k = bu->n = 0; k < kn; k++) {
+		if((bu->a[k] == ((uint64_t)-1)) || (bu->a[k] == 0)) continue;
+		i = k; occ = nc = 0;
+		while (i != (uint32_t)-1) {
+			if(rg->seq[rch->bb.a[i].hid].del) break;
+
+			occ++; if(is_contain_r((*ri), rch->bb.a[i].hid)) nc++; 
+
+			i = rch->bb.a[i].pidx;
+			if(i != ((uint32_t)-1)) {
+				if(!is_rg_connect(rg, (uint64_t)((rch->bb.a[i].hid<<1)|(rch->bb.a[i].rev)), 
+						(uint64_t)((rch->bb.a[rch->bb.a[i].aidx].hid<<1)|(rch->bb.a[rch->bb.a[i].aidx].rev)))) {
+					i = ((uint32_t)-1);
+				}
+			}
+		}
+		// for (i = k, occ = nc = 0; (i != (uint32_t)-1) && (!rg->seq[rch->bb.a[i].hid].del); i = rch->bb.a[i].pidx) {
+		// 	occ++; if(is_contain_r((*ri), rch->bb.a[i].hid)) nc++; 
+		// }
+		assert(occ == bu->a[k]);
+		if(nc) bu->a[bu->n++] = k;
+		// if(debug_out) {
+			// if(occ != rch->bb.n) {
+			// 	fprintf(stderr, "-[M::%s::] ulid::%lu, nc::%lu, occ::%lu, k::%ld, rch->bb.n::%u\n", 
+			// 	__func__, ulid, nc, occ, k, (uint32_t)rch->bb.n);
+			// }
+		// }
+	}
+	return bu->n;///# chains have contained reads
+}
+
+uint32_t extract_ccov(ma_hit_t *in, uc_block_t *rovlp, const asg_t *rg, ul_ov_t *res, uint32_t adjust_rev, int64_t min_ovlp, int64_t max_hang)
+{
+	uint64_t qn, tn; int32_t r = 1; asg_arc_t e; 
+	qn = Get_qn((*in)); tn = Get_tn((*in));
+	if(rg->seq[tn].del) return 0;	
+	
+	r = ma_hit2arc(in, Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), max_hang, asm_opt.max_hang_rate, min_ovlp, &e);
+	if(r != MA_HT_QCONT) return 0; ///qn is contained in tn
+
+	uint64_t ori = in->rev, ts = rovlp->ts, te = rovlp->te, tl;
+	if(ori) {
+		ts = rg->seq[qn].len - rovlp->te;
+        te = rg->seq[qn].len - rovlp->ts;
+	}
+
+	ts += in->ts; te += in->ts; 
+	if(te <= ts) return 0;
+	tl = te - ts; tl = tl*0.01; if(tl > 8) tl = 8;
+	if((ts >= (rg->seq[tn].len+tl)) || (te >= (rg->seq[tn].len+tl))) return 0;
+	if(ts > rg->seq[tn].len) ts = rg->seq[tn].len;
+	if(te > rg->seq[tn].len) te = rg->seq[tn].len;
+	if(te <= ts) return 0;
+
+	memset(res, 0, sizeof(*res));
+	res->qn = 0; res->qs = rovlp->qs; res->qe = rovlp->qe; 
+	res->tn = tn; res->ts = ts; res->te = te; 
+	res->el = rovlp->el; res->rev = (rovlp->rev == ori?0:1);
+	if(adjust_rev && res->rev) {///for linear chaining
+		res->ts = rg->seq[tn].len - te;
+		res->te = rg->seq[tn].len - ts;
+	}
+	return 1;
+}
+
+uint32_t extract_nccov(ma_hit_t *in, uc_block_t *rovlp, const asg_t *rg, ul_ov_t *res, uint32_t adjust_rev, int64_t min_ovlp, int64_t max_hang)
+{
+	uint64_t tn; int64_t os, oe, s_shift, e_shift, tt, qs, qe, ts, te; 
+	tn = Get_tn((*in)); if(rg->seq[tn].del) return 0;	
+	os = MAX(rovlp->ts, Get_qs((*in)));
+	oe = MIN(rovlp->te, Get_qe((*in)));
+	if(oe <= os) return 0;
+
+	///[os, oe) -> rovlp->t*
+	s_shift = get_offset_adjust(os-rovlp->ts, rovlp->te-rovlp->ts, rovlp->qe-rovlp->qs);
+    e_shift = get_offset_adjust(rovlp->te-oe, rovlp->te-rovlp->ts, rovlp->qe-rovlp->qs);
+	if(rovlp->rev) {
+		tt = s_shift; s_shift = e_shift; e_shift = tt;
+	}
+	qs = rovlp->qs + s_shift; qe = ((int64_t)rovlp->qe)-e_shift; 
+	if(qs >= qe) return 0;
+
+	///[os, oe) -> in->q*
+	s_shift = get_offset_adjust(os-Get_qs((*in)), Get_qe((*in))-Get_qs((*in)), Get_te((*in))-Get_ts((*in)));
+	e_shift = get_offset_adjust(Get_qe((*in))-oe, Get_qe((*in))-Get_qs((*in)), Get_te((*in))-Get_ts((*in)));
+	if(in->rev) {
+		tt = s_shift; s_shift = e_shift; e_shift = tt;
+	}
+	ts = Get_ts((*in)) + s_shift; te = ((int64_t)Get_te((*in)))-e_shift;
+	if(ts >= te) return 0;
+
+	memset(res, 0, sizeof(*res));
+	res->qn = 0; res->qs = qs; res->qe = qe; 
+	res->tn = tn; res->ts = ts; res->te = te; 
+	res->el = rovlp->el; res->rev = ((rovlp->rev == in->rev)?0:1);
+	if(adjust_rev && res->rev) {///for linear chaining
+		res->ts = rg->seq[tn].len - te;
+		res->te = rg->seq[tn].len - ts;
+	}
+
+	// fprintf(stderr, "\n[M::%s::id->%u::%c] q::[%u, %u), t::[%u, %u)\n", 
+	// 			__func__, rovlp->hid, "+-"[rovlp->rev], rovlp->qs, rovlp->qe, rovlp->ts, rovlp->te);
+	// fprintf(stderr, "+[M::%s::] qn::%u, q::[%u, %u), %c, tn::%u, t::[%u, %u)\n", 
+	// 			__func__, Get_qn((*in)), Get_qs((*in)), in->qe, "+-"[in->rev], in->tn, in->ts, in->te);
+	// fprintf(stderr, "-[M::%s::] qn::%u, q::[%u, %u), %c, tn::%u, t::[%u, %u)\n", 
+	// 			__func__, res->qn, res->qs, res->qe, "+-"[res->rev], res->tn, res->ts, res->te);
+	return 1;
+}
+
+
+void collect_nc_ovlps(ul_vec_t *rch, uint32_t rch_i, kv_ul_ov_t *res, const ug_opt_t *uopt, const asg_t *rg, R_to_U *ri, int64_t ulid, uint64_t *mqs, uint64_t *mqe)
+{
+	uint64_t i, occ, nc, k; ma_hit_t_alloc* src; 
+	int64_t min_ovlp = uopt->min_ovlp, max_hang = uopt->max_hang; ul_ov_t p;
+	res->n = occ = nc = 0; (*mqs) = (*mqe) = (uint64_t)-1;
+	i = rch_i;
+	// for (i = rch_i; (i != (uint32_t)-1) && (!rg->seq[rch->bb.a[i].hid].del); i = rch->bb.a[i].pidx) {
+	while(i != (uint32_t)-1) {
+		if(rg->seq[rch->bb.a[i].hid].del) break;
+		// if(ulid == 442462) {		
+		// 	fprintf(stderr, "\n>>>[M::%s::i->%lu::id->%u::%.*s]\n", __func__, i, rch->bb.a[i].hid, 
+		// 		(int)Get_NAME_LENGTH(R_INF, rch->bb.a[i].hid), Get_NAME(R_INF, rch->bb.a[i].hid));
+		// 	fprintf(stderr, "*[M::%s::id->%u::%.*s]\tq::[%u, %u),\tt::[%u, %u),\tis_rev::%u,\tis_cr::%u\n", 
+		// 			__func__, rch->bb.a[i].hid, 
+		// 			(int)Get_NAME_LENGTH(R_INF, rch->bb.a[i].hid), Get_NAME(R_INF, rch->bb.a[i].hid), 
+		// 			rch->bb.a[i].qs, rch->bb.a[i].qe, rch->bb.a[i].ts, rch->bb.a[i].te, rch->bb.a[i].rev,
+		// 			!!(is_contain_r((*(ri)), rch->bb.a[i].hid)));
+		// }
+		
+		if(is_contain_r((*ri), rch->bb.a[i].hid)) {
+			src = &(uopt->sources[rch->bb.a[i].hid]);
+			for (k = 0; k < src->length; k++) {
+				if(rg->seq[Get_tn(src->buffer[k])].del) continue;///tn must exist
+				// if(!extract_ccov(&(src->buffer[k]), &(rch->bb.a[i]), rg, &p, 1, min_ovlp, max_hang)) continue;
+				// if(ulid == 442462) {	
+				// 	fprintf(stderr, "+[M::%s::id->%u::%.*s]\tq::[%u, %u),\tql::%u,\tt::[%u, %u),\ttl::%u,\tis_rev::%u,\tis_cr::%u\n", 
+				// 	__func__, src->buffer[k].tn, 
+				// 	(int)Get_NAME_LENGTH(R_INF, src->buffer[k].tn), Get_NAME(R_INF, src->buffer[k].tn), 
+				// 	Get_qs(src->buffer[k]), Get_qe(src->buffer[k]), rg->seq[Get_qn(src->buffer[k])].len,
+				// 	Get_ts(src->buffer[k]), Get_te(src->buffer[k]), rg->seq[Get_tn(src->buffer[k])].len,
+				// 	src->buffer[k].rev, !!(is_contain_r((*(ri)), src->buffer[k].tn)));
+				// }
+				if(!extract_nccov(&(src->buffer[k]), &(rch->bb.a[i]), rg, &p, 1, min_ovlp, max_hang)) continue;
+				// if(ulid == 442462) {	
+				// 	fprintf(stderr, "-[M::%s::id->%u::%.*s]\tq::[%u, %u),\tt::[%u, %u),\tis_rev::%u,\tis_cr::%u\n", 
+				// 	__func__, p.tn, (int)Get_NAME_LENGTH(R_INF, p.tn), Get_NAME(R_INF, p.tn), 
+				// 	p.qs, p.qe, p.ts, p.te, p.rev, !!(is_contain_r((*(ri)), p.tn)));
+				// 	// fprintf(stderr, "cc[M::%s::id->%u::%.*s] q::[%u, %u), t::[%u, %u), is_cr::%u, del::%u\n", 
+				// 	// __func__, p.tn, (int)Get_NAME_LENGTH(R_INF, p.tn), Get_NAME(R_INF, p.tn), 
+				// 	// p.qs, p.qe, p.ts, p.te, !!(is_contain_r((*(ri)), p.tn)), rg->seq[p.tn].del);
+				// }
+				p.el = 0; p.tn <<= 1; p.tn |= p.rev; p.qn = i;//for linear chain
+            	kv_push(ul_ov_t, *res, p);
+			}
+			nc++;
+		}
+		///push itself into the chain
+		memset(&p, 0, sizeof(p));
+		p.qn = 0; p.qs = rch->bb.a[i].qs; p.qe = rch->bb.a[i].qe; 
+    	p.tn = rch->bb.a[i].hid; p.ts = rch->bb.a[i].ts; p.te = rch->bb.a[i].te; 
+    	p.el = 1; p.rev = rch->bb.a[i].rev;
+		if(p.rev) {///for linear chaining
+			p.ts = rg->seq[p.tn].len - rch->bb.a[i].te;
+			p.te = rg->seq[p.tn].len - rch->bb.a[i].ts;
+		}
+		p.el = 1; p.tn <<= 1; p.tn |= p.rev; p.qn = i;//for linear chain
+		kv_push(ul_ov_t, *res, p);
+		if(((*mqs) == ((uint64_t)-1)) || ((*mqs) > rch->bb.a[i].qs)) (*mqs) = rch->bb.a[i].qs;
+		if(((*mqe) == ((uint64_t)-1)) || ((*mqe) < rch->bb.a[i].qe)) (*mqe) = rch->bb.a[i].qe;
+		occ++; 
+
+		i = rch->bb.a[i].pidx;
+		if(i != ((uint32_t)-1)) {
+			if(!is_rg_connect(rg, (uint64_t)((rch->bb.a[i].hid<<1)|(rch->bb.a[i].rev)), 
+					(uint64_t)((rch->bb.a[rch->bb.a[i].aidx].hid<<1)|(rch->bb.a[rch->bb.a[i].aidx].rev)))) {
+				i = ((uint32_t)-1);
+			}
+		}
+	}
+	assert(occ > 1 && nc > 0);
+}
+
+
+inline int64_t comput_rlinear_sc(ul_ov_t *li, ul_ov_t *lj, int64_t jidx, int32_t *bq, int32_t *bt, double diff_ec_ul, int64_t bw)
+{	///li is the suffix of lj; sorted by qe, so li->qe >= lj->qe, so li->te >= lj->te
+	if(li->te < lj->te) return INT32_MIN;
+	int64_t dq, dt, dd, mm, os, oe;
+	oe = MIN((li->qe), (lj->qe)); os = MAX((li->qs), (lj->qs));
+	if(oe <= os) return INT32_MIN;
+	oe = MIN((li->te), (lj->te)); os = MAX((li->ts), (lj->ts));
+	if(oe <= os) return INT32_MIN;
+
+	dq = li->qe - lj->qs; dt = li->te - lj->ts;
+	dd = (dq>dt? dq-dt:dt-dq);
+	mm = MAX(dq, dt); mm *= diff_ec_ul; if(mm < bw) mm = bw;
+	if(dd > mm) return INT32_MIN;
+	dq = (int64_t)li->qe - bq[jidx];
+	dt = (int64_t)li->te - bt[jidx];
+	if(dq < ((int64_t)(li->qe - li->qs))) dq = li->qe - li->qs;
+	if(dt < ((int64_t)(li->te - li->ts))) dt = li->te - li->ts;
+	return MIN(dq, dt);
+}
+
+uint64_t linear_rchain_dp_adv(ul_ov_t *ch, int64_t ch_n, ul_ov_t *sv, const ug_opt_t *uopt, int64_t bw, 
+double diff_ec_ul, int64_t qlen, int64_t max_skip, int64_t max_iter, int64_t max_dis, Chain_Data* dp, const asg_t *rg, int64_t ulid)
+{	///all in[].el must be 1
+	if(ch_n == 0) return 0;
+    int64_t i, j, k, sc, csc, mm_sc, mm_idx, its, ite, max; int32_t *f, *bq, *bt;
+    ul_ov_t *li = NULL, *lj = NULL; int64_t *p, *t, st, plus, max_ii, n_skip, end_j; 
+	resize_Chain_Data(dp, ch_n, NULL);
+	t = dp->tmp; f = dp->score; p = dp->pre; bq = dp->indels; bt = dp->self_length;
+
+    radix_sort_ul_ov_srt_qe(ch, ch + ch_n);
+	for (i = 1, j = 0; i <= ch_n; i++) {
+        if (i == ch_n || ch[i].qe != ch[j].qe) {
+            if(i - j > 1) radix_sort_ul_ov_srt_qs(ch+j, ch+i);
+            j = i;
+        }
+    }
+
+	memset(t, 0, (ch_n*sizeof((*t))));
+	for (i = st = plus = 0, max_ii = -1; i < ch_n; ++i) {
+		li = &(ch[i]); csc = MIN((li->qe-li->qs), (li->te-li->ts));
+		mm_sc = /**csc**/-1; mm_idx = -1; n_skip = 0; end_j = -1; 
+		st = (i<max_iter)?(0):(i-max_iter);
+		for (j = i-1; j >= st; --j) {
+			lj = &(ch[j]); 
+			if(lj->qe <= li->qs) break;
+			sc = comput_rlinear_sc(li, lj, j, bq, bt, diff_ec_ul, bw); ///should allow contain
+			if(sc == INT32_MIN) continue;
+			if(sc > mm_sc) {
+				mm_sc = sc, mm_idx = j;
+				if (n_skip > 0) --n_skip;
+			} else if (t[j] == i) {
+				if (++n_skip > max_skip)
+					break;
+			}
+			if (p[j] >= 0) t[p[j]] = i;
+		}
+		end_j = j;
+		if (max_ii < 0 || (ch[i].qe>(ch[max_ii].qe+max_dis))) {//too long
+			max = INT32_MIN; max_ii = -1;
+            for (j = i - 1; (j >= st) && (ch[i].qe<=(max_dis+ch[j].qe)); --j) {
+                if (max < f[j]) {
+                    max = f[j], max_ii = j;
+                }
+            }
+		}
+		
+        if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii]
+			lj = &(ch[max_ii]); 
+			if(lj->qe > li->qs) {
+				sc = comput_rlinear_sc(li, lj, max_ii, bq, bt, diff_ec_ul, bw); ///should allow contain
+				if(sc != INT32_MIN) {
+					if(sc > mm_sc) {
+						mm_sc = sc; mm_idx = max_ii;
+					}
+				}
+			}
+		}
+
+		if(mm_idx < 0) mm_sc = csc;
+		f[i] = mm_sc; p[i] = mm_idx; 
+		bq[i] = li->qs; bt[i] = li->ts;
+		if(mm_idx >= 0) {
+			if(bq[i] > bq[mm_idx]) bq[i] = bq[mm_idx];
+			if(bt[i] > bt[mm_idx]) bt[i] = bt[mm_idx];
+		} 
+
+        if ((max_ii < 0) || ((ch[i].qe<=max_dis+ch[max_ii].qe) && (f[max_ii]<f[i]))) {
+            max_ii = i;
+        }
+		li->sec = (mm_idx<0?0x3FFFFFFF:i-mm_idx); sv[i] = *li;
+		// if(ulid == 442462) {    
+		// 	fprintf(stderr, "[M::%s::id->%u::%.*s]\tq::[%u, %u),\tt::[%u, %u),\tis_rev::%u\n", 
+		// 	__func__, li->tn, (int)Get_NAME_LENGTH(R_INF, li->tn), Get_NAME(R_INF, li->tn), 
+		// 	li->qs, li->qe, li->ts, li->te, li->rev);
+		// }
+
+	}
+	// radix_sort_gfa64i
+	for (i = 0; i < ch_n; ++i) {
+		t[i] = ((uint64_t)f[i])<<32; t[i] += i; bt[i] = 0;
+	}
+	radix_sort_gfa64i(t, t+ch_n);
+	int64_t n_u, z;
+	for (z = ch_n-1, n_u = 0; z >= 0; --z) {
+		k = (uint32_t)t[z];
+		if(bt[k]) continue;
+		i = k; ch[n_u]=sv[i]; sc = f[i];
+		for (;i>=0;) {
+			if(sv[i].qs < ch[n_u].qs) ch[n_u].qs = sv[i].qs;
+			if(sv[i].ts < ch[n_u].ts) ch[n_u].ts = sv[i].ts;
+			if(sv[i].qe > ch[n_u].qe) ch[n_u].qe = sv[i].qe;
+			if(sv[i].te > ch[n_u].te) ch[n_u].te = sv[i].te;
+			// ch[n_u].qn = i;//start idx of read alignment in chain
+			bt[i] = 1; i = p[i];
+		}
+		adjust_rev_tse(&(ch[n_u]), rg->seq[ch[n_u].tn].len, &its, &ite);
+		ch[n_u].ts = its; ch[n_u].te = ite; ch[n_u].sec = (sc>0x3FFFFFFF?0x3FFFFFFF:sc);
+		ch[n_u].qn = 0; 
+		n_u++;
+	}
+	return n_u;
+}
+
+
+void gen_linear_rchains(kv_ul_ov_t *res, kv_ul_ov_t *buf, const asg_t *rg, const ug_opt_t *uopt, int64_t bw, 
+double diff_ec_ul, int64_t qlen, Chain_Data* dp, int64_t ulid)
+{
+	uint64_t k, l, z, an, m, bn = buf->n; 
+    radix_sort_ul_ov_srt_tn(res->a, res->a + res->n);
+    ///after this function, res keeps unitig alignment, while buf keeps read alignments
+    for (k = 1, l = m = 0; k <= res->n; k++) {
+        if(k == res->n || res->a[k].tn != res->a[l].tn) {///qn <- (tn|rev)
+			for (z = l; z < k; z++) res->a[z].tn>>=1;
+			kv_resize(ul_ov_t, *buf, bn+k-l); 
+			an = l + linear_rchain_dp_adv(res->a+l, k-l, buf->a+bn, uopt, bw, diff_ec_ul, qlen, 
+			UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, rg, ulid);
+			for (z = l; z < an; z++) res->a[m++] = res->a[z];
+			l = k;
+		}
+    }
+    res->n = m; 
+	// kv_resize(ul_ov_t, *buf, bn+res->n);
+	// int64_t iqs, iqe, its, ite;
+	// for (k = 0, z = bn; k < res->n; k++) {
+	// 	buf->a[z] = res->a[k]; res->a[k].qn = z; 
+	// 	extend_end_coord(NULL, &(res->a[k]), qlen, rg->seq[res->a[k].tn].len, &iqs, &iqe, &its, &ite);
+	// 	res->a[k].qs = iqs; res->a[k].qe = iqe; res->a[k].ts = its; res->a[k].te = ite;
+	// 	z++;
+	// }
+}
+
+int64_t gconnect_test(const asg_t *g, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq)
+{
+	int64_t dt = -1, dif, mm; 
+	uint32_t nv, i; asg_arc_t *av = NULL; ///ma_hit_t *x = NULL;
+	nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+	for (i = 0; i < nv; i++) {
+		if(av[i].del || av[i].v != w) continue;
+		dt = av[i].ol; 
+		break;
+	}
+
+	if(dt < 0) return 0;
+	dif = (dq>dt? dq-dt:dt-dq);
+	mm = MAX(dq, dt); mm *= diff_ec_ul; if(mm < bw) mm = bw;
+	if(dif <= mm) return 1;
+	return 0;
+}
+
+uint64_t gen_cns_chain_linear(ul_ov_t *a, int64_t a_n, /**ul_ov_t *ab,**/ const asg_t *rg, R_to_U *ri, int64_t qlen, int64_t bw, double diff_thre, Chain_Data* dp, 
+int64_t max_skip, int64_t max_iter, int64_t max_dis, uint64_t mqs, uint64_t mqe, uint64_t ulid)
+{
+	if(a_n == 0) return 0;
+	uint32_t li_v, lj_v; int32_t *f, *c_n, *len; int64_t *p, *t, st, max_ii, max, qo, cL, sn, ln, csn, mm_sn, cln, mm_ln; 
+	int64_t mm_ovlp, x, i, j, sc, csc, mm_sc, mm_idx, n_skip, end_j, ch_sc, cn_sn, ch_ln, ch_i; ul_ov_t *li, *lj; 
+	resize_Chain_Data(dp, a_n, NULL); 
+	t = dp->tmp; f = dp->score; p = dp->pre; c_n = dp->occ; len = dp->indels;
+        
+	radix_sort_ul_ov_srt_qe(a, a + a_n);
+    for (i = 1, j = 0; i <= a_n; i++) {
+        if (i == a_n || a[i].qe != a[j].qe) {
+            if(i - j > 1) radix_sort_ul_ov_srt_qs(a+j, a+i);
+            j = i;
+        }
+    }
+
+	memset(t, 0, (a_n*sizeof((*t)))); 
+	ch_sc = ch_i = cn_sn = INT32_MIN; ch_ln = INT32_MAX;
+	for (i = st = 0, max_ii = -1; i < a_n; ++i) {
+        li = &(a[i]); 
+        mm_ovlp = max_ovlp(rg, ((li->tn<<1)|li->rev)^1);        
+        x = (li->qs + mm_ovlp)*diff_thre; 
+        if(x < bw) x = bw; 
+        x += li->qs + mm_ovlp;
+        if (x > qlen+1) x = qlen+1;
+        x = find_ul_ov_max(i, a, x+G_CHAIN_INDEL); 
+
+		csc = 0; csn = 1; cln = 0; 
+		if(is_contain_r((*ri), li->tn)) {csc = -1; csn = 0; cln = rg->seq[li->tn].len;}
+        mm_sc = INT32_MIN+1; mm_sn = INT32_MIN+1; mm_ln = INT32_MAX; 
+		mm_idx = -1; n_skip = 0; end_j = -1; 
+
+		li_v = (li->tn<<1)|li->rev; li_v ^= 1;
+        if ((x-st) > max_iter) st = x-max_iter;
+		for (j = x; j >= st; --j) { // collect potential destination vertices
+			lj = &(a[j]); 
+			if(lj->qe+G_CHAIN_INDEL <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore
+			if(f[j] == INT32_MIN) continue;///could not reach the left end
+			sc = sn = ln = INT32_MIN;
+			if(lj->qs <= li->qs && lj->qe <= li->qe) {///not contained
+				lj_v = (lj->tn<<1)|lj->rev; lj_v ^= 1;
+				qo = infer_rovlp(li, lj, NULL, NULL, NULL, NULL);
+				if(gconnect_test(rg, li_v, lj_v, bw, diff_thre, qo)) {
+					sc = f[j] + csc; sn = c_n[j] + csn; ln = len[j] + cln;
+				}
+			}
+			if(sc == INT32_MIN) continue;
+			if((sc > mm_sc) || ((sc == mm_sc) && (sn > mm_sn)) || 
+				((sc == mm_sc) && (sn == mm_sn) && (ln < mm_ln))) {
+				mm_sc = sc, mm_idx = j, mm_sn = sn, mm_ln = ln;
+				if (n_skip > 0) --n_skip;
+			} else if (t[j] == i) {
+				if (++n_skip > max_skip) break;
+			}
+			if (p[j] >= 0) t[p[j]] = i;
+		}
+        
+        end_j = j;
+        if (max_ii < 0 || (a[i].qe > (a[max_ii].qe+max_dis))) {//too long
+            max = INT32_MIN; max_ii = -1;
+            for (j = i - 1; (j >= st) && (a[i].qe<=(max_dis+a[j].qe)); --j) {
+                if (max < f[j]) {
+                    max = f[j], max_ii = j;
+                }
+            }
+        }
+                
+        if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii]
+			lj = &(a[max_ii]); 
+			if(((lj->qe+G_CHAIN_INDEL)>li->qs) && (lj->qs<=li->qs) && (lj->qe<=li->qe) && (f[max_ii]!=INT32_MIN)) {
+				lj_v = (lj->tn<<1)|lj->rev; lj_v ^= 1; sc = sn = ln = INT32_MIN;
+				qo = infer_rovlp(li, lj, NULL, NULL, NULL, NULL);
+				if(gconnect_test(rg, li_v, lj_v, bw, diff_thre, qo)) {
+					sc = f[max_ii] + csc; sn = c_n[max_ii] + csn; ln = len[max_ii] + cln;
+				}
+				if(sc != INT32_MIN) {
+					if((sc > mm_sc) || ((sc == mm_sc) && (sn > mm_sn)) || 
+						((sc == mm_sc) && (sn == mm_sn) && (ln < mm_ln))) {
+						mm_sc = sc, mm_idx = max_ii, mm_sn = sn, mm_ln = ln;
+					}
+				}
+			}
+		}
+
+		if((mm_idx == -1) && (((int64_t)li->qs) > ((int64_t)(mqs+bw)))) {
+			mm_sc = mm_idx = mm_sn = INT32_MIN; mm_ln = INT32_MAX;
+		}
+		if((mm_sc==(INT32_MIN+1))) {
+			mm_sc = csc; mm_sn = csn; mm_ln = cln;
+		}
+
+        
+        f[i] = mm_sc; p[i] = mm_idx; c_n[i] = mm_sn; len[i] = mm_ln;
+        if ((max_ii < 0) || ((a[i].qe<=max_dis+a[max_ii].qe) && (f[max_ii]<f[i]))) {
+            max_ii = i;
+        }
+
+		if((mm_sc!=INT32_MIN) && (((int64_t)(a[i].qe+bw))>=((int64_t)(mqe)))) {
+			if((mm_sc > ch_sc) || ((mm_sc == ch_sc) && (mm_sn > cn_sn)) || 
+				((mm_sc == ch_sc) && (mm_sn == cn_sn) && (mm_ln < ch_ln))) {
+				ch_sc = mm_sc; ch_i = i; cn_sn = mm_sn; ch_ln = mm_ln;
+			}
+		}
+		// if(ulid == 442462) {
+		// 	fprintf(stderr, "[M::%s::id->%u::%.*s] i::%ld, q::[%u, %u), t::[%u, %u), is_cr::%u, mm_idx::%ld, end::%u\n", 
+		// 	__func__, li->tn, (int)Get_NAME_LENGTH(R_INF, li->tn), Get_NAME(R_INF, li->tn), i, 
+		// 	li->qs, li->qe, li->ts, li->te, !!(is_contain_r((*ri), li->tn)), mm_idx, !!(((int64_t)(a[i].qe+bw))>=((int64_t)(mqe))));
+		// }
+    }
+	// if(ulid == 442462) {
+	// 	fprintf(stderr, "[M::%s::] ch_i::%ld, mqs::%lu, mqe::%lu\n", 
+	// 	__func__, ch_i, mqs, mqe);
+	// }
+	if(ch_i == INT32_MIN) return 0;
+
+	i = ch_i; cL = 0;
+    while (i >= 0) {t[cL++] = i; i = p[i];}
+	for (i = 0; i < cL; i++) a[i] = a[t[cL-i-1]];
+	return cL;
+}
+
+void gen_contain_consensus_chain(ul_vec_t *rch, uint32_t rch_i, kv_ul_ov_t *idx, kv_ul_ov_t *dump,
+const ug_opt_t *uopt, const asg_t *rg, R_to_U *ri, int64_t bw, double diff_ec_ul, int64_t ulid, Chain_Data* dp)
+{
+	uint64_t dn = dump->n, i, mqs, mqe; int64_t iqs, iqe, its, ite; idx->n = 0;
+	collect_nc_ovlps(rch, rch_i, idx, uopt, rg, ri, ulid, &mqs, &mqe);
+	assert(idx->n > 1); assert((mqs != ((uint64_t)-1)) && (mqe != ((uint64_t)-1)) && (mqe > mqs));
+	gen_linear_rchains(idx, dump, rg, uopt, bw, diff_ec_ul, rch->rlen, dp, ulid);
+	assert(idx->n);
+	idx->n = gen_cns_chain_linear(idx->a, idx->n, /**dump->a+dn,**/ rg, ri, rch->rlen, bw, diff_ec_ul, dp, UG_SKIP_N, UG_ITER_N, UG_DIS_N, mqs, mqe, ulid);
+	// fprintf(stderr, "[M::%s::] idx->n::%u, mqs::%lu, mqe::%lu\n", __func__, (uint32_t)idx->n, mqs, mqe);
+	if(idx->n) {
+		kv_resize(ul_ov_t, *dump, dn + idx->n);
+		///mask existing overlaps
+		for (i = rch_i; (i != (uint32_t)-1) && (!rg->seq[rch->bb.a[i].hid].del); i = rch->bb.a[i].pidx) {
+			dump->a[i].tn = dump->a[i].qn = (uint32_t)-1;
+		}
+		for (i = 0; i < idx->n; i++) {
+			dump->a[dn] = idx->a[i];
+			extend_end_coord(NULL, &(dump->a[dn]), rch->rlen, rg->seq[dump->a[dn].tn].len, &iqs, &iqe, &its, &ite);
+			dump->a[dn].qs = iqs; dump->a[dn].qe = iqe; dump->a[dn].ts = its; dump->a[dn].te = ite;
+			dump->a[dn].qn = ((i>0)?(dn-1):((uint32_t)-1)); dn++;
+		}
+		dump->n = dn;
+	}
+}
+
+uint32_t update_consensus_chain(const ug_opt_t *uopt, const asg_t *rg, kv_ul_ov_t *idx, kv_ul_ov_t *dump, ul_vec_t *rch, asg64_v *b, R_to_U *ri)
+{
+	uint64_t *idm, k, l, kn, dn, cc = 0;
+	assert(dump->n >= rch->bb.n);
+	kv_resize(uint64_t, *b, dump->n); idm = b->a;
+
+	for (k = kn = 0; k < rch->bb.n; k++) {
+		idm[k] = (uint64_t)-1;
+		if((dump->a[k].qn == ((uint32_t)-1)) && (dump->a[k].tn == ((uint32_t)-1))) continue;
+		dump->a[k].qn = rch->bb.a[dump->a[k].qn].pidx; idm[k] = kn; kn++;
+	}
+	for (; k < dump->n; k++) {
+		idm[k] = kn; kn++;
+	}
+	// fprintf(stderr, "[M::%s::kn->%lu] dump->n::%u\n", __func__, kn, (uint32_t)dump->n);
+	
+	dn = dump->n; dump->n = 0; 
+	kv_resize(ul_ov_t, *idx, kn); idx->n = 0;
+	for (k = 0; k < dn; k++) {
+		if(idm[k] == ((uint64_t)-1)) continue;
+		dump->a[idm[k]] = dump->a[k]; kn--; dump->n++;
+		if(dump->a[idm[k]].qn != ((uint32_t)-1)) {
+			dump->a[idm[k]].qn = (uint32_t)idm[dump->a[idm[k]].qn];
+		}
+		kv_push(ul_ov_t, *idx, dump->a[idm[k]]); idx->a[idx->n-1].qn = idm[k];
+	}
+	assert(kn == 0); 
+
+	radix_sort_ul_ov_srt_qe(idx->a, idx->a + idx->n);
+	for (k = 1, l = 0; k <= idx->n; k++) {
+        if (k == idx->n || idx->a[l].qe != idx->a[k].qe) {
+            if(k - l > 1) radix_sort_ul_ov_srt_qs(idx->a+l, idx->a+k);
+            l = k;
+        }
+    }
+	for (k = 0; k < idx->n; k++) dump->a[idx->a[k].qn].tn = k;
+
+	for (k = 0; k < idx->n; k++) {
+		idx->a[k].qn = ((dump->a[idx->a[k].qn].qn!=((uint32_t)-1))?
+									(dump->a[dump->a[idx->a[k].qn].qn].tn):((uint32_t)-1));
+		idx->a[k].el = 1;
+	}
+
+	uc_block_t *z, *p; int64_t tt;
+	kv_resize(uc_block_t, rch->bb, idx->n); rch->bb.n = 0;
+	for (k = 0; k < idx->n; k++) {
+		// fprintf(stderr, "+k::%lu[M::%s::id->%u] q::[%u, %u), t::[%u, %u), is_cr::%u\n", 
+		// k, __func__, idx->a[k].tn, idx->a[k].qs, idx->a[k].qe, idx->a[k].ts, idx->a[k].te, 
+		// !!(is_contain_r((*ri), idx->a[k].tn)));
+		kv_pushp(uc_block_t, rch->bb, &z); memset(z, 0, sizeof((*z)));
+		z->hid = idx->a[k].tn; z->rev = idx->a[k].rev; 
+    	z->pchain = 1; z->base = 0; z->el = 1;
+		z->qs = idx->a[k].qs; z->qe = idx->a[k].qe;
+		z->te = idx->a[k].te; z->ts = idx->a[k].ts;
+        z->pidx = idx->a[k].qn; z->pdis = z->aidx = (uint32_t)-1;
+		
+	}
+	for (k = cc = 0; k < rch->bb.n; k++) {
+		if(is_contain_r((*ri), rch->bb.a[k].hid)) cc++;
+		// fprintf(stderr, "-k::%lu[M::%s::id->%u] q::[%u, %u), t::[%u, %u), is_cr::%u\n", 
+		// k, __func__, rch->bb.a[k].hid, rch->bb.a[k].qs, rch->bb.a[k].qe, rch->bb.a[k].ts, rch->bb.a[k].te, 
+		// !!(is_contain_r((*ri), rch->bb.a[k].hid)));
+		if(rch->bb.a[k].pidx == (uint32_t)-1) continue;
+		z = &(rch->bb.a[k]); p = &(rch->bb.a[z->pidx]);
+		assert(p->aidx == (uint32_t)-1); p->aidx = k;
+		tt = g_adjacent_dis_mul(NULL, uopt->sources, uopt->max_hang, uopt->min_ovlp, ((z->hid<<1)|((uint32_t)z->rev))^1, ((p->hid<<1)|((uint32_t)p->rev))^1);
+		if(tt >= 0) z->pdis = tt; ///assert(tt >= 0);
+	}
+
+	///debug_ssb
+	// for (k = 0; k < rch->bb.n; k++) {
+	// 	if(rch->bb.a[k].pidx != (uint32_t)-1) {
+	// 		if((rch->bb.a[k].pidx >= 0) && (rch->bb.a[k].pidx < rch->bb.n) && 
+	// 					(rch->bb.a[rch->bb.a[k].pidx].aidx == k)) {
+	// 			;
+	// 		} else {
+	// 			fprintf(stderr, "[M::%s::k->%lu] +rch->bb.n::%u\n", __func__, k, (uint32_t)rch->bb.n);
+	// 			for (l = 0; l < rch->bb.n; l++) {
+	// 				fprintf(stderr, "[M::%.*s::k->%lu] q::[%u, %u), t::[%u, %u), is_cr::%u, pidx::%u, aidx::%u\n", 
+	// 				(int)Get_NAME_LENGTH(R_INF, rch->bb.a[l].hid), Get_NAME(R_INF, rch->bb.a[l].hid), l, 
+	// 				rch->bb.a[l].qs, rch->bb.a[l].qe, rch->bb.a[l].ts, rch->bb.a[l].te, 
+	// 				!!(is_contain_r((*ri), rch->bb.a[l].hid)), 
+	// 				rch->bb.a[l].pidx, rch->bb.a[l].aidx);
+	// 			}
+	// 			exit(1);
+	// 		}
+	// 	}
+
+	// 	if(rch->bb.a[k].aidx != (uint32_t)-1) {
+	// 		if((rch->bb.a[k].aidx >= 0) && (rch->bb.a[k].aidx < rch->bb.n) && 
+	// 					(rch->bb.a[rch->bb.a[k].aidx].pidx == k)) {
+	// 			;
+	// 		} else {
+	// 			fprintf(stderr, "[M::%s::k->%lu] -rch->bb.n::%u\n", __func__, k, (uint32_t)rch->bb.n);
+	// 			for (l = 0; l < rch->bb.n; l++) {
+	// 				fprintf(stderr, "[M::%.*s::k->%lu] q::[%u, %u), t::[%u, %u), is_cr::%u, pidx::%u, aidx::%u\n", 
+	// 				(int)Get_NAME_LENGTH(R_INF, rch->bb.a[l].hid), Get_NAME(R_INF, rch->bb.a[l].hid), l, 
+	// 				rch->bb.a[l].qs, rch->bb.a[l].qe, rch->bb.a[l].ts, rch->bb.a[l].te, 
+	// 				!!(is_contain_r((*ri), rch->bb.a[l].hid)), 
+	// 				rch->bb.a[l].pidx, rch->bb.a[l].aidx);
+	// 			}
+	// 			exit(1);
+	// 		}
+	// 	}
+	// }
+
+	return ((cc==0)?1:0);
+}
+
+static void worker_for_contain_consensus(void *data, long i, int tid)
+{
+    ul_vec_t *p = &(UL_INF.a[i]); asg64_v b0; 
+    utepdat_t *s = (utepdat_t*)data; uint32_t ff, k; 
+
+	// if(i != 304) return;
+    copy_asg_arr(b0, s->ll[tid].srt.a); 
+    ff = refine_contain_consensus_chain(s->rg, p, s->uopt->ruIndex, &b0, i);
+    copy_asg_arr(s->ll[tid].srt.a, b0);
+    // fprintf(stderr, "[M::%s::%.*s(id:%ld), len:%u] ff:%u\n", __func__, 
+    // UL_INF.nid.a[i].n, UL_INF.nid.a[i].a, i, p->rlen, ff);
+	// if(debug_out) {
+	// 	for (k = 0; k < p->bb.n; k++) {
+	// 		fprintf(stderr, "[M::%.*s] q::[%u, %u), t::[%u, %u), is_cr::%u, pidx::%u, aidx::%u\n", 
+	// 		(int)Get_NAME_LENGTH(R_INF, p->bb.a[k].hid), Get_NAME(R_INF, p->bb.a[k].hid), 
+	// 		p->bb.a[k].qs, p->bb.a[k].qe, p->bb.a[k].ts, p->bb.a[k].te, 
+	// 		!!(is_contain_r((*(s->uopt->ruIndex)), p->bb.a[k].hid)), p->bb.a[k].pidx, p->bb.a[k].aidx);
+    // 	}
+	// 	fprintf(stderr, "***[M::%s::%.*s(id:%ld), len:%u] ff:%u, sp_chn::%u, p->bb.n::%u\n\n", __func__, 
+    // 	UL_INF.nid.a[i].n, UL_INF.nid.a[i].a, i, p->rlen, ff, (uint32_t)s->ll[tid].srt.a.n, (uint32_t)p->bb.n);
+	// }
+    if(!ff) return;
+	
+	// char *as = NULL; 
+	// asprintf(&as, "\n[M::%s]\trid::%ld\tlen::%lu\tname::%.*s\tb0.n::%u\n", 
+	// __func__, s->id+i, s->len[i], (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a, (uint32_t)b0.n);
+    // push_vlog(&(overall_zdbg->a[s->id+i]), as); free(as); as = NULL;
+
+
+
+    kv_ul_ov_t *idx = &(s->ll[tid].lo), *dump = &(s->ll[tid].tk); ul_ov_t *z;
+    // if(p->dd == 1) return; //fully aligned
+    // if(p->bb.n == 1 && p->bb.a[0].base) return;///no alignment
+    // if(p->bb.n == 0) return;///no alignment
+    s->hab[tid]->num_read_base++; 
+
+	// if(i == 442462) {
+	// 	fprintf(stderr, "\n[M::%s::%.*s(id:%ld), len:%u] ff:%u, sp_chn::%u, p->bb.n::%u\n", __func__, 
+    // 	UL_INF.nid.a[i].n, UL_INF.nid.a[i].a, i, p->rlen, ff, (uint32_t)s->ll[tid].srt.a.n, (uint32_t)p->bb.n);
+	// }
+
+	// fprintf(stderr, "[M::%s]\tp->bb.n::%u\n", __func__, (uint32_t)p->bb.n);
+    kv_resize(ul_ov_t, *dump, p->bb.n);
+    for (k = dump->n = 0; k < p->bb.n; k++) {
+        z = &(dump->a[dump->n++]); ///memset(z, 0, sizeof((*z)));
+        z->qn = k; z->qs = p->bb.a[k].qs; z->qe = p->bb.a[k].qe; 
+        z->tn = p->bb.a[k].hid; z->ts = p->bb.a[k].ts; z->te = p->bb.a[k].te; 
+        z->el = 1; z->rev = p->bb.a[k].rev; z->sec = 0;
+		// if(i == 442462) {
+		// 	fprintf(stderr, "[M::%s::id->%u::%.*s] q::[%u, %u), t::[%u, %u), is_rev::%u, is_cr::%u\n", 
+		// 	__func__, z->tn, (int)Get_NAME_LENGTH(R_INF, z->tn), Get_NAME(R_INF, z->tn), 
+		// 	z->qs, z->qe, z->ts, z->te, z->rev, !!(is_contain_r((*(s->uopt->ruIndex)), z->tn)));
+		// }
+    }
+    
+    for (k = 0; k < s->ll[tid].srt.a.n; k++) {
+        gen_contain_consensus_chain(p, s->ll[tid].srt.a.a[k], idx, dump, s->uopt, s->rg, s->uopt->ruIndex, G_CHAIN_BW, s->opt->diff_ec_ul, i, &(s->hab[tid]->clist.chainDP));
+    }
+	
+    copy_asg_arr(b0, s->ll[tid].srt.a);
+	ff = update_consensus_chain(s->uopt, s->rg, idx, dump, p, &b0, s->uopt->ruIndex);
+    copy_asg_arr(s->ll[tid].srt.a, b0);
+	s->hab[tid]->num_correct_base += ff;
+	// fprintf(stderr, "[M::%s::%.*s(id:%ld), len:%u] ffa->%u, sp_chn::%u\n", __func__, 
+    // UL_INF.nid.a[i].n, UL_INF.nid.a[i].a, i, p->rlen, ff, (uint32_t)s->ll[tid].srt.a.n);
+
+    // s->hab[tid]->num_correct_base += direct_gchain(s->buf[tid], p, &(s->ll[tid]), &(s->gdp[tid]), &(s->sps[tid]), &(s->hab[tid]->hap), s->uu, s->uopt, G_CHAIN_BW, s->opt->diff_ec_ul, UG_SKIP, i, &(s->hab[tid]->clist.chainDP), s->rg, ((ff==2)?1:0));
+    // gl_chain_refine_advance(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, km);
+}
+
+uint32_t extract_contain_tig(asg_t *g, uint64_t v0, R_to_U *ri, uint64_t offset, uint64_t must_tip, kv_ul_ov_t *res)
+{
+	uint32_t v = v0, w, kv, kw, is_tip = 0, l = offset;
+	ul_ov_t p; memset(&p, 0, sizeof(p));
+
+	while (1) {
+		if(!is_contain_r((*ri), (v>>1))) return 0;
+		kv = get_arcs(g, v, &w, 1);        
+		if(kv > 1) return 0;
+
+		p.qn = 0; p.qs = l; p.qe = l + g->seq[v>>1].len; 
+		p.tn = v>>1; p.ts = 0; p.te = g->seq[v>>1].len;
+		p.el = 0; p.rev = v&1; kv_push(ul_ov_t, *res, p);
+
+		if(kv == 0) {is_tip = 1; break;}
+
+		l += asg_arc_len(g->arc[w]);
+		w = g->arc[w].v;
+        ///kv must be 1 here
+        kw = get_arcs(g, w^1, NULL, 0);
+		assert(kw >= 1);
+		if(kw > 1) {
+			p.qn = 0; p.qs = l; p.qe = l + g->seq[w>>1].len; 
+			p.tn = w>>1; p.ts = 0; p.te = g->seq[w>>1].len;
+			p.el = 1; p.rev = w&1; kv_push(ul_ov_t, *res, p);
+			break;
+		}
+        v = w;
+        if(v == v0) return 0;
+    }
+	if((must_tip) && (!is_tip)) return 0;
+	return 1;
+}
+
+uint32_t gen_dup_path(const asg_t *g, R_to_U *ri, kv_ul_ov_t *res, uint64_t v)
+{
+    asg_t *rg = (asg_t *)g; res->n = 0;
+	if(rg->seq[v>>1].del) return res->n;
+	uint64_t nv = asg_arc_n(rg, v); asg_arc_t *av = asg_arc_a(rg, v);
+	uint64_t k, kv, w = (uint64_t)-1, kw, l, rn, z; ul_ov_t p; memset(&p, 0, sizeof(p));
+
+	// for (k = 0; k < nv; k++) {
+	// 	if(av[k].del) continue;
+	// 	asg_arc_t *aw = asg_arc_a(rg, av[k].v);
+	// 	uint64_t nw = asg_arc_n(rg, av[k].v);
+	// 	for (l = 0; l < nw; l++) {
+	// 		if(aw[l].del) continue;
+	// 		for (rn = 0; rn < nv; rn++) {
+	// 			if(av[rn].del) continue;
+	// 			if(av[rn].v == aw[l].v) break;
+	// 		}
+	// 		if(rn < nv) {
+	// 			fprintf(stderr, "[M::%s] v0::%lu, w::%u, v1::%u\n", __func__, v, av[k].v, aw[l].v);
+	// 		}
+	// 	}
+	// }
+	
+
+
+
+	for (k = kv = 0; k < nv && kv < 2; k++) {
+		if(!(av[k].del)) kv++; 
+	}
+	// if(v == 7723) {
+	// 	fprintf(stderr, "[M::%s] v0::%lu, kv::%lu\n", __func__, v, kv);
+	// }
+	if(kv > 1) {
+		for (k = 0; k < nv; k++) {
+			if(av[k].del) continue;
+			w = av[k].v; kw = get_arcs(rg, w^1, NULL, 0);
+			// if(v == 7723) {
+			// 	fprintf(stderr, "*[M::%s] v0::%lu, av[%lu].v::%u, kw::%lu,\n", __func__, v, k, av[k].v, kw);
+			// }
+			if(kw > 1) {///kv > 1 && kw > 1
+				if((v>>1) != (w>>1)) {///cannot handle if the prefix and suffix nodes are the same
+					l = 0;
+					p.qn = 0; p.qs = l; p.qe = l + rg->seq[v>>1].len; 
+					p.tn = v>>1; p.ts = 0; p.te = rg->seq[v>>1].len;
+					p.el = 1; p.rev = v&1; kv_push(ul_ov_t, *res, p);
+					// if(v == 7723) {
+					// 	fprintf(stderr, "*[M::%s] v0::%lu, p.tn::%u, p.rev::%u\n", __func__, v, p.tn, p.rev);
+					// }
+
+					l += asg_arc_len(av[k]);
+					p.qn = 0; p.qs = l; p.qe = l + rg->seq[w>>1].len; 
+					p.tn = w>>1; p.ts = 0; p.te = rg->seq[w>>1].len;
+					p.el = 1; p.rev = w&1; kv_push(ul_ov_t, *res, p);
+
+					// if(v == 7723) {
+					// 	fprintf(stderr, "*[M::%s] v0::%lu, p.tn::%u, p.rev::%u\n", __func__, v, p.tn, p.rev);
+					// }
+
+					p.tn = (uint32_t)-1; p.qn = ((uint64_t)(av-g->arc+k));
+					kv_push(ul_ov_t, *res, p);
+				}
+			} else if((kw == 1) && (is_contain_r((*ri), (w>>1)))) {///kv > 1 && kw == 1 && w is contained
+				rn = res->n; l = 0;
+				p.qn = 0; p.qs = l; p.qe = l + rg->seq[v>>1].len; 
+				p.tn = v>>1; p.ts = 0; p.te = rg->seq[v>>1].len;
+				p.el = 1; p.rev = v&1; kv_push(ul_ov_t, *res, p);
+				l += asg_arc_len(av[k]);
+				if(!extract_contain_tig(rg, w, ri, l, 0, res)) {
+					res->n = rn;
+				} else {
+					for (z = rn+1; z < res->n; z++) {
+						if(res->a[z].tn == (v>>1)) break;
+					}
+					if(z >= res->n) {///cannot handle circles
+						p.tn = (uint32_t)-1; p.qn = ((uint64_t)(av-g->arc+k));
+						kv_push(ul_ov_t, *res, p);
+					} else {
+						res->n = rn;
+					}
+				}
+			}
+		}
+	} else if(((kv == 1) || (kv == 0)) && (!get_arcs(rg, v^1, NULL, 0)) && (is_contain_r((*ri), (v>>1)))) {
+		///1. kv == 1 && kv^ == 0 && v is contained
+		///2. kv == 0 && kv^ == 0 && v is contained; an isolated single contained read
+		rn = res->n; l = 0;
+		if(!extract_contain_tig(rg, v, ri, l, 1, res)) {
+			res->n = rn;
+		} else {
+			for (z = rn+1; z < res->n; z++) {
+				if(res->a[z].tn == (v>>1)) break;
+			}
+			if(z >= res->n) {///cannot handle circles
+				p.tn = (uint32_t)-1; p.qn = (uint32_t)-1;
+				if(kv > 0) {
+					for (k = 0; k < nv; k++) {
+						if(!(av[k].del)) break; 
+					}
+					assert(k < nv);
+					p.qn = ((uint64_t)(av-g->arc+k));
+				}
+
+				kv_push(ul_ov_t, *res, p);
+			} else {
+				res->n = rn;
+			}
+		}
+	} 
+	// if(v == 7723) {
+	// 	fprintf(stderr, "[M::%s] v0::%lu, res->n::%lu,\n", __func__, v, (uint64_t)res->n);
+	// }
+    return res->n;///# chains have contained reads
+}
+
+void collect_pp_ovlps(ul_ov_t *a, uint64_t a_n, kv_ul_ov_t *res, const ug_opt_t *uopt, const asg_t *rg, uint64_t rlen, int64_t ulid, uint64_t *mqs, uint64_t *mqe)
+{
+	uint64_t i, k; ma_hit_t_alloc* src; uc_block_t rch; ul_ov_t p, *st, *et;
+    int64_t min_ovlp = uopt->min_ovlp, max_hang = uopt->max_hang, iqs, iqe; st = et = NULL;
+    res->n = 0; (*mqs) = (*mqe) = (uint64_t)-1; memset(&rch, 0, sizeof(rch));
+	if(a_n && a[0].el) st = &(a[0]);
+	if(a_n > 1 && a[a_n-1].el) et = &(a[a_n-1]);
+	// if(ulid == 7723 && a_n == 2 && a[0].tn == (7723>>1) && a[a_n-1].tn == (7829>>1)) {
+	// 	fprintf(stderr, "[M::%s] ulid::%ld, st::%ld, st_q::[%ld, %ld), et::%ld, et_q::[%ld, %ld)\n", 
+	// 	__func__, ulid, st?(int64_t)st->tn:-1, st?(int64_t)st->qs:-1, st?(int64_t)st->qe:-1, 
+	// 	et?(int64_t)et->tn:-1, et?(int64_t)et->qs:-1, et?(int64_t)et->qe:-1);
+	// }
+	for (i = 0; i < a_n; i++) {
+		src = &(uopt->sources[a[i].tn]);
+		rch.hid = a[i].tn; rch.rev = a[i].rev;
+		rch.qs = a[i].qs; rch.qe = a[i].qe;
+		rch.ts = a[i].ts; rch.te = a[i].te;
+		for (k = 0; k < src->length; k++) {
+            if(rg->seq[Get_tn(src->buffer[k])].del) continue;///tn must exist
+			if(st && st->tn == Get_tn(src->buffer[k])) continue;//not useful
+			if(et && et->tn == Get_tn(src->buffer[k])) continue;//not useful
+			if(!extract_nccov(&(src->buffer[k]), &rch, rg, &p, 0, min_ovlp, max_hang)) continue;
+			
+			extend_end_coord(NULL, &p, rlen, rg->seq[p.tn].len, &iqs, &iqe, NULL, NULL);
+			if(st && iqs <= (int64_t)st->qe && iqe <= (int64_t)st->qe) continue;//not useful
+			if(et && iqs >= (int64_t)et->qs && iqe >= (int64_t)et->qs) continue;//not useful
+			
+			if(p.rev) {
+				iqs = rg->seq[p.tn].len - p.te;
+				iqe = rg->seq[p.tn].len - p.ts;
+				p.ts = iqs; p.te = iqe;
+			}
+			// if(ulid == 7723 && a_n == 2 && a[0].tn == (7723>>1) && a[a_n-1].tn == (7829>>1)) {
+			// 	fprintf(stderr, "[M::%s] ulid::%ld, p.tn::%u, pq::[%u, %u), p.el::%u\n", 
+			// 	__func__, ulid, p.tn, p.qs, p.qe, p.el);
+			// }
+			p.el = 0; p.tn <<= 1; p.tn |= p.rev; p.qn = i;//for linear chain
+            kv_push(ul_ov_t, *res, p);
+		}
+
+		///it is unnecessary to push any non-end reads; 
+		///since the read graph will remove any old read; each read/node is unique within the read graph
+		if(a[i].el) {
+			///push itself into the chain
+			p = a[i]; 
+			// if(ulid == 7723 && a_n == 2 && a[0].tn == (7723>>1) && a[a_n-1].tn == (7829>>1)) {
+			// 	fprintf(stderr, "[M::%s] ulid::%ld, p.tn::%u, pq::[%u, %u), p.el::%u\n", 
+			// 	__func__, ulid, p.tn, p.qs, p.qe, p.el);
+			// }
+			p.tn <<= 1; p.tn |= p.rev; p.qn = i;//for linear chain
+			kv_push(ul_ov_t, *res, p);
+		}
+        if(((*mqs) == ((uint64_t)-1)) || ((*mqs) > a[i].qs)) (*mqs) = a[i].qs;
+        if(((*mqe) == ((uint64_t)-1)) || ((*mqe) < a[i].qe)) (*mqe) = a[i].qe;
+		
+	}
+	if(st) (*mqs) = st->qs; if(et) (*mqe) = et->qe;
+}
+
+uint64_t gen_cns_chain_linear_hard(ul_ov_t *a, int64_t a_n, const asg_t *rg, int64_t qlen, int64_t bw, double diff_thre, Chain_Data* dp, 
+int64_t max_skip, int64_t max_iter, int64_t max_dis, uint64_t mqs, uint64_t mqe, uint64_t ulid, uint64_t ltn, uint64_t rtn)
+{
+    if(a_n == 0) return 0;
+    uint32_t li_v, lj_v; int32_t *f; int64_t *p, *t, st, max_ii, max, qo, cL; ul_ov_t *li, *lj; 
+    int64_t mm_ovlp, x, i, j, sc, csc, mm_sc, mm_idx, n_skip, end_j, ch_sc, ch_i; 
+    resize_Chain_Data(dp, a_n, NULL); t = dp->tmp; f = dp->score; p = dp->pre;
+        
+    radix_sort_ul_ov_srt_qe(a, a + a_n);
+    for (i = 1, j = 0; i <= a_n; i++) {
+        if (i == a_n || a[i].qe != a[j].qe) {
+            if(i - j > 1) radix_sort_ul_ov_srt_qs(a+j, a+i);
+            j = i;
+        }
+    }
+
+    memset(t, 0, (a_n*sizeof((*t)))); 
+    ch_sc = ch_i = INT32_MIN; 
+    for (i = st = 0, max_ii = -1; i < a_n; ++i) {
+        li = &(a[i]); 
+        mm_ovlp = max_ovlp(rg, ((li->tn<<1)|li->rev)^1);        
+        x = (li->qs + mm_ovlp)*diff_thre; 
+        if(x < bw) x = bw; 
+        x += li->qs + mm_ovlp;
+        if (x > qlen+1) x = qlen+1;
+        x = find_ul_ov_max(i, a, x+G_CHAIN_INDEL); 
+
+        csc = 1; if(li->el) csc = 0;
+        mm_sc = INT32_MIN+1; mm_idx = -1; n_skip = 0; end_j = -1; 
+
+        li_v = (li->tn<<1)|li->rev; li_v ^= 1;
+        if ((x-st) > max_iter) st = x-max_iter;
+        for (j = x; j >= st; --j) { // collect potential destination vertices
+            lj = &(a[j]); 
+            if(lj->qe+G_CHAIN_INDEL <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore
+            if(f[j] == INT32_MIN) continue;///could not reach the left end
+            sc = INT32_MIN;
+            if(lj->qs <= li->qs && lj->qe <= li->qe) {///not contained
+                lj_v = (lj->tn<<1)|lj->rev; lj_v ^= 1;
+                qo = infer_rovlp(li, lj, NULL, NULL, NULL, NULL);
+                if(gconnect_test(rg, li_v, lj_v, bw, diff_thre, qo)) {
+                    sc = f[j] + csc; 
+                }
+            }
+            if(sc == INT32_MIN) continue;
+            if(sc > mm_sc) {
+                mm_sc = sc, mm_idx = j;
+                if (n_skip > 0) --n_skip;
+            } else if (t[j] == i) {
+                if (++n_skip > max_skip) break;
+            }
+            if (p[j] >= 0) t[p[j]] = i;
+        }
+        
+        end_j = j;
+        if (max_ii < 0 || (a[i].qe > (a[max_ii].qe+max_dis))) {//too long
+            max = INT32_MIN; max_ii = -1;
+            for (j = i - 1; (j >= st) && (a[i].qe<=(max_dis+a[j].qe)); --j) {
+                if (max < f[j]) {
+                    max = f[j], max_ii = j;
+                }
+            }
+        }
+                
+        if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii]
+            lj = &(a[max_ii]); 
+            if(((lj->qe+G_CHAIN_INDEL)>li->qs) && (lj->qs<=li->qs) && (lj->qe<=li->qe) && (f[max_ii]!=INT32_MIN)) {
+                lj_v = (lj->tn<<1)|lj->rev; lj_v ^= 1; sc = INT32_MIN;
+                qo = infer_rovlp(li, lj, NULL, NULL, NULL, NULL);
+                if(gconnect_test(rg, li_v, lj_v, bw, diff_thre, qo)) {
+                    sc = f[max_ii] + csc; 
+                }
+                if(sc != INT32_MIN) {
+                    if(sc > mm_sc) {
+                        mm_sc = sc, mm_idx = max_ii;
+                    }
+                }
+            }
+        }
+
+		if(mm_idx == -1) {
+			if((ltn) != ((uint64_t)-1)) {
+				if(ltn != li->tn) mm_sc = mm_idx = INT32_MIN;
+			} else {
+				if(((int64_t)li->qs) > ((int64_t)(mqs+bw))) mm_sc = mm_idx = INT32_MIN;
+			}
+		}
+        
+        if((mm_sc==(INT32_MIN+1))) {
+            mm_sc = csc; 
+        }
+
+        
+        f[i] = mm_sc; p[i] = mm_idx; 
+
+		// if(ulid == 7723 && ltn == (7723>>1) && rtn == (7829>>1)) {
+		// 	fprintf(stderr, "[M::%s::id->%u::%.*s] i::%ld, q::[%u, %u), t::[%u, %u), mm_idx::%ld, mm_sc::%ld\n", 
+		// 	__func__, li->tn, (int)Get_NAME_LENGTH(R_INF, li->tn), Get_NAME(R_INF, li->tn), i, 
+		// 	li->qs, li->qe, li->ts, li->te, mm_idx, mm_sc);
+        // }
+
+		if(mm_sc == INT32_MIN) continue;
+		if((rtn != ((uint64_t)-1)) && (rtn != li->tn)) continue;
+		if((rtn == ((uint64_t)-1)) && (((int64_t)(a[i].qe+bw))<((int64_t)(mqe)))) continue;
+        if(mm_sc >= ch_sc) {
+			ch_sc = mm_sc; ch_i = i;
+		}
+		if ((max_ii < 0) || (f[max_ii]<f[i])) {
+            max_ii = i;
+        }
+    }
+    // if(ulid == 442462) {
+    //  fprintf(stderr, "[M::%s::] ch_i::%ld, mqs::%lu, mqe::%lu\n", 
+    //  __func__, ch_i, mqs, mqe);
+    // }
+    if(ch_i == INT32_MIN) return 0;
+
+    i = ch_i; cL = 0;
+    while (i >= 0) {t[cL++] = i; i = p[i];}
+    for (i = 0; i < cL; i++) a[i] = a[t[cL-i-1]];
+    return cL;
+}
+
+void gen_linear_rchains_dedup(kv_ul_ov_t *res, kv_ul_ov_t *buf, const asg_t *rg, const ug_opt_t *uopt, int64_t bw, 
+double diff_ec_ul, int64_t qlen, Chain_Data* dp, int64_t ulid, uint64_t *flt, uint64_t flt_n)
+{
+	uint64_t k, l, z, an, m, bn = buf->n, fi, tn; 
+    radix_sort_ul_ov_srt_tn(res->a, res->a + res->n);
+    ///after this function, res keeps unitig alignment, while buf keeps read alignments
+    for (k = 1, l = m = fi = 0; k <= res->n; k++) {
+        if(k == res->n || res->a[k].tn != res->a[l].tn) {///qn <- (tn|rev)
+			tn = res->a[l].tn>>1;
+			for (z = l; z < k; z++) res->a[z].tn>>=1;
+			for (; (fi<flt_n) && (flt[fi]<tn); fi++);
+			if((fi>=flt_n) || (flt[fi]!=tn)) {
+				kv_resize(ul_ov_t, *buf, bn+k-l); 
+				an = l + linear_rchain_dp_adv(res->a+l, k-l, buf->a+bn, uopt, bw, diff_ec_ul, qlen, 
+				UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, rg, ulid);
+				for (z = l; z < an; z++) res->a[m++] = res->a[z];
+			}
+			l = k;
+		}
+    }
+    res->n = m; 
+}
+
+void renew_consensus_chain(kv_ul_ov_t *in, uint32_t in_s, uint32_t in_e, kv_ul_ov_t *idx, asg64_v *ou, R_to_U *ri,
+const ug_opt_t *uopt, const asg_t *rg, int64_t bw, double diff_ec_ul, uint32_t arc_id, Chain_Data* dp, uint64_t rlen, uint32_t ulid)
+{	
+	if(in_e <= in_s) return;
+	uint64_t mqs, mqe, k, l, ou_n = ou->n, stn, etn;
+	collect_pp_ovlps(in->a + in_s, in_e - in_s, idx, uopt, rg, rlen, ulid, &mqs, &mqe);
+	// if(ulid == 7723 && in_e-in_s == 2 && in->a[in_s].tn == (7723>>1) && in->a[in_e-1].tn == (7829>>1)) {
+	// 	fprintf(stderr, "-0-[M::%s] ulid::%u, idx->n::%u\n", __func__, ulid, (uint32_t)idx->n);
+	// }
+	if(!idx->n) return;///it is possible
+	kv_resize(uint64_t, *ou, (ou_n+(in_e-in_s)));
+	for (k = in_s; k < in_e; k++) {
+		if(in->a[k].el) continue;
+		ou->a[ou_n++] = in->a[k].tn;
+	}
+	assert(ou_n <= ou->m);
+	radix_sort_gfa64(ou->a+ou->n, ou->a+ou_n);
+	gen_linear_rchains_dedup(idx, in, rg, uopt, bw, diff_ec_ul, rlen, dp, ulid, ou->a+ou->n, ou_n-ou->n);
+	// if(ulid == 7723 && in_e-in_s == 2 && in->a[in_s].tn == (7723>>1) && in->a[in_e-1].tn == (7829>>1)) {
+	// 	fprintf(stderr, "-1-[M::%s] ulid::%u, idx->n::%u\n", __func__, ulid, (uint32_t)idx->n);
+	// }
+	if(!idx->n) return;///it is possible
+	stn = etn = (uint64_t)-1;
+	if(in_e > in_s && in->a[in_s].el) stn = in->a[in_s].tn;
+	if(in_e-in_s>1 && in->a[in_e-1].el) etn = in->a[in_e-1].tn;
+	idx->n = gen_cns_chain_linear_hard(idx->a, idx->n, rg, rlen, bw, diff_ec_ul, dp, UG_SKIP_N, UG_ITER_N, UG_DIS_N, mqs, mqe, ulid, stn, etn);
+	// if(ulid == 7723 && in_e-in_s == 2 && in->a[in_s].tn == (7723>>1) && in->a[in_e-1].tn == (7829>>1)) {
+	// 	fprintf(stderr, "-2-[M::%s] ulid::%u, idx->n::%u\n", __func__, ulid, (uint32_t)idx->n);
+	// }
+	if(!idx->n) return;
+	assert((stn==((uint64_t)-1))||((idx->a[0].tn == stn) && (idx->a[0].el)));
+	assert((etn==((uint64_t)-1))||((idx->a[idx->n-1].tn == etn) && (idx->a[idx->n-1].el)));
+	if((stn != ((uint64_t)-1)) && (etn != ((uint64_t)-1)) && idx->n < 2) return;///could happen for circle
+	if(((stn != ((uint64_t)-1)) || (etn != ((uint64_t)-1))) && idx->n < 1) return;///not sure if it will happen
+	ou_n = ou->n;
+	for (k = in_s; k < in_e; k++) {
+		if(in->a[k].el) continue;
+		kv_push(uint64_t, *ou, in->a[k].tn);
+	}
+	for (k = 0; k < idx->n; k++) {
+		if(idx->a[k].el) continue;
+		kv_push(uint64_t, *ou, idx->a[k].tn);
+	}
+	if(ou->n == ou_n) return;
+	radix_sort_gfa64(ou->a+ou_n, ou->a+ou->n);
+	for (k = ou_n+1, l = ou_n; k <= ou->n; k++) {
+		if(k == ou->n || ou->a[l] != ou->a[k]) {
+			if(k-l > 1) break;
+			l = k;
+		}
+	}
+	assert(k > ou->n);
+	ou->n = ou_n;
+	l = arc_id; 
+	if(arc_id == (uint32_t)-1) {///this is an isloated node
+		l = idx->a[0].tn; l |= ((uint64_t)0x80000000);
+	}
+	l <<= 32; l |= ((uint64_t)((uint32_t)-1));
+	kv_push(uint64_t, *ou, l);
+	for (k = 0; k < idx->n; k++) {
+		l = ((uint64_t)(idx->a[k].tn<<1))|((uint64_t)idx->a[k].rev);
+		if(idx->a[k].el) l += ((uint64_t)0x100000000);
+		kv_push(uint64_t, *ou, l);
+	}
+	kv_push(uint64_t, *ou, ((uint64_t)-1));
+	for (k = in_s; k < in_e; k++) {
+		l = ((uint64_t)(in->a[k].tn<<1))|((uint64_t)in->a[k].rev);
+		if(in->a[k].el) l += ((uint64_t)0x100000000);
+		kv_push(uint64_t, *ou, l);
+	}
+}
+
+static void worker_for_contain_dedup(void *data, long i, int tid)
+{
+    utepdat_t *s = (utepdat_t*)data; asg64_v ou; 
+	kv_ul_ov_t *idx = &(s->ll[tid].lo), *dump = &(s->ll[tid].tk); 
+	uint32_t k, v = i, l, dump_n; 
+
+	if(!gen_dup_path(s->rg, s->uopt->ruIndex, dump, v)) return;
+	// if(v == 7723) {
+	// 	fprintf(stderr, "[M::%s] v::%u, dump->n::%u\n", __func__, v, (uint32_t)dump->n);
+	// }
+
+
+	copy_asg_arr(ou, s->ll[tid].srt.a); 
+	for (k = 1, l = 0, dump_n = dump->n; k <= dump_n; k++) {
+		if(k == dump_n || (dump->a[k].tn == ((uint32_t)-1))) {
+			assert(k > l || k == dump_n); 
+			if(k > l) {
+				assert(dump->a[k].tn == ((uint32_t)-1));
+				renew_consensus_chain(dump, l, k, idx, &ou, s->uopt->ruIndex, s->uopt, s->rg, G_CHAIN_BW, 
+							0.02/**s->opt->diff_ec_ul**/, dump->a[k].qn, &(s->hab[tid]->clist.chainDP), dump->a[k-1].qe, v);
+			}
+			l = k + 1;
+		}
+	}
+	copy_asg_arr(s->ll[tid].srt.a, ou);
+}
+
 void detect_outlier_len(const char* cmd)
 {
 	uint64_t k;
@@ -15467,6 +16714,460 @@ uint64_t work_ul_gchains(uldat_t *sl)
 	return s.n;
 }
 
+
+uint64_t work_ul_gchains_consensus(uldat_t *sl)
+{
+	utepdat_t s; uint64_t i; memset(&s, 0, sizeof(s));
+	s.id = 0; s.uopt = sl->uopt; s.rg = sl->rg; s.opt = sl->opt; ///s.ug = sl->ug; s.uu = sl->uu; 
+	//  CALLOC(s.buf, sl->n_thread); CALLOC(s.gdp, sl->n_thread); CALLOC(s.mzs, sl->n_thread); 
+	CALLOC(s.hab, sl->n_thread); CALLOC(s.ll, sl->n_thread); ///CALLOC(s.sps, sl->n_thread); 
+
+	for (i = 0; i < sl->n_thread; ++i) {
+		s.hab[i] = ha_ovec_init(0, 0, 1); ///s.buf[i] = mg_tbuf_init();
+	}
+
+	// detect_outlier_len("+++work_ul_gchains");
+
+	//debug
+    // overall_zdbg = init_mul_debug_prt_t(UL_INF.n);
+
+	kt_for(sl->n_thread, worker_for_contain_consensus, &s, UL_INF.n);
+
+	// detect_outlier_len("---work_ul_gchains");
+
+	for (i = 0; i < sl->n_thread; ++i) {
+		s.sum_len += s.hab[i]->num_read_base; s.n += s.hab[i]->num_correct_base;
+		// hc_gdpchain_destroy(&(s.gdp[i])); kv_destroy(s.mzs[i]); mg_tbuf_destroy(s.buf[i]); 
+		ha_ovec_destroy(s.hab[i]); hc_glchain_destroy(&(s.ll[i])); ///kv_destroy(s.sps[i]);
+	}
+
+	// free(s.buf); free(s.gdp); free(s.mzs); 
+	free(s.hab); free(s.ll); ///free(s.sps);
+	fprintf(stderr, "[M::%s::] # try:%d, # done:%d\n", __func__, s.sum_len, s.n); 
+	return s.sum_len;
+}
+
+
+void update_contain_dedup_arr(asg64_v *in, uint64_t prefix, asg64_v *res)
+{
+	if(in->n == 0) return;
+	uint64_t i, m;
+	for (i = 0; i < in->n; i++) {
+		if(in->a[i] == ((uint64_t)-1)) continue;
+		if(((uint32_t)in->a[i]) == ((uint32_t)-1)) {
+			m = prefix; m <<= 32; m += i;
+			kv_push(uint64_t, *res, m);
+		}
+	}
+}
+
+uint64_t path_del_topo(asg_t *g, uint64_t v)
+{
+	uint64_t nv, nw, kw, k, z, w; asg_arc_t *av, *aw;
+	nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+	for (k = 0; k < nv; k++){
+		if(av[k].del) continue;
+		w = av[k].v^1;
+		nw = asg_arc_n(g, w); aw = asg_arc_a(g, w);
+		for (z = kw = 0; z < nw; z++) {
+			if(aw[z].del) continue;
+			if((aw[z].v>>1) != (v>>1)) continue;
+			kw++;
+		}
+		if(kw <= 0) return 0; 
+	}
+	return 1;
+}
+
+uint64_t tes_path_connec(asg_t *g, uint64_t sn, uint64_t en, uint64_t *pa, uint64_t pn, uint64_t is_single_path, uint64_t *min_ou)
+{
+	// fprintf(stderr, "-0-[M::%s::]\n", __func__); 
+	if(min_ou) (*min_ou) = 0;
+	if((sn != (uint64_t)-1) && ((g->seq[sn>>1].del)||(get_arcs(g, sn, NULL, 0)<2))) return 0;
+	if((en != (uint64_t)-1) && ((g->seq[en>>1].del)||((get_arcs(g, en^1, NULL, 0)<2)))) return 0;
+	uint64_t nv, k, v, w, i, ou = (uint64_t)-1; asg_arc_t *av;
+	if((sn != (uint64_t)-1) && (en != (uint64_t)-1) && (!pn)) {
+		v = sn; w = en;
+		nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+		for (k = 0; k < nv; k++) {
+			if(av[k].del) continue;
+			if(av[k].v == w) {
+				if(((uint64_t)av[k].ou) < ou) ou = av[k].ou;
+				break;
+			}
+		}
+		if(k >= nv) return 0;
+		
+		v = en^1; w = sn^1;
+		nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+		for (k = 0; k < nv; k++) {
+			if(av[k].del) continue;
+			if(av[k].v == w) {
+				if(((uint64_t)av[k].ou) < ou) ou = av[k].ou;
+				break;
+			}
+		}
+		if(k >= nv) return 0;
+
+		// if(get_arcs(g, sn, NULL, 0) < 2) return 0;
+		// if(get_arcs(g, en^1, NULL, 0) < 2) return 0;
+		if(min_ou) (*min_ou) = ou;
+		return 1;
+	}
+
+	// fprintf(stderr, "-1-[M::%s::]\n", __func__); 
+	if(sn != (uint64_t)-1) {
+		v = sn; i = 0;
+	} else {
+		v = pa[0]; i = 1; assert(!(pa[0]>>32));
+	}
+	// fprintf(stderr, "-2-[M::%s::]\n", __func__); 
+	for (; i < pn; i++) {
+		w = pa[i]; assert(!(pa[i]>>32));
+		if(g->seq[v>>1].del || g->seq[w>>1].del) return 0;
+
+		nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+		for (k = 0; k < nv; k++) {
+			if(av[k].del) continue;
+			if(av[k].v == w) {
+				if(((uint64_t)av[k].ou) < ou) ou = av[k].ou;
+				break;
+			}
+		}
+		if(k >= nv) return 0;
+
+		nv = asg_arc_n(g, w^1); av = asg_arc_a(g, w^1);
+		for (k = 0; k < nv; k++) {
+			if(av[k].del) continue;
+			if(av[k].v == (v^1)) {
+				if(((uint64_t)av[k].ou) < ou) ou = av[k].ou;
+				break;
+			}
+		}
+		if(k >= nv) return 0;
+
+		if(is_single_path) {
+			if(v != sn && v != en && get_arcs(g, v, NULL, 0) != 1) return 0;
+			if(w != sn && w != en && get_arcs(g, w^1, NULL, 0) != 1) return 0;
+		}
+
+		v = w;
+	}
+	// fprintf(stderr, "-3-[M::%s::]\n", __func__); 
+	if(en != (uint64_t)-1) {
+		w = en;
+		if(g->seq[v>>1].del || g->seq[w>>1].del) return 0;
+		// fprintf(stderr, "-3a-[M::%s::]\n", __func__); 
+		nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+		for (k = 0; k < nv; k++) {
+			if(av[k].del) continue;
+			if(av[k].v == w) {
+				if(((uint64_t)av[k].ou) < ou) ou = av[k].ou;
+				break;
+			}
+		}
+		if(k >= nv) return 0;
+		// fprintf(stderr, "-3b-[M::%s::]\n", __func__); 
+		nv = asg_arc_n(g, w^1); av = asg_arc_a(g, w^1);
+		for (k = 0; k < nv; k++) {
+			if(av[k].del) continue;
+			if(av[k].v == (v^1)) {
+				if(((uint64_t)av[k].ou) < ou) ou = av[k].ou;
+				break;
+			}
+		}
+		if(k >= nv) return 0;
+
+		if(is_single_path) {
+			if(v != sn && v != en && get_arcs(g, v, NULL, 0) != 1) return 0;
+			if(w != sn && w != en && get_arcs(g, w^1, NULL, 0) != 1) return 0;
+		}
+	}
+
+
+
+	// fprintf(stderr, "-4-[M::%s::]\n", __func__); 
+	if(!path_del_topo(g, pa[0]^1)) return 0;
+	// if(sn != ((uint64_t)-1)) {
+	// 	v = sn; w = pa[0];
+	// 	nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+	// 	for (k = 0; k < nv; k++) {
+	// 		if(av[k].del) continue;
+	// 		if(av[k].v == w) break;
+	// 	}
+	// 	if(k >= nv) return 0;
+
+	// 	v = pa[0]^1; w = sn^1;
+	// 	nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+	// 	for (k = 0; k < nv; k++) {
+	// 		if(av[k].del) continue;
+	// 		if(av[k].v == w) break;
+	// 	}
+	// 	if(k >= nv) return 0;
+	// }
+
+
+
+	if(!path_del_topo(g, pa[pn-1])) return 0;
+	// if(en != ((uint64_t)-1)) {
+	// 	v = pa[pn-1]; w = en;
+	// 	nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+	// 	for (k = 0; k < nv; k++) {
+	// 		if(av[k].del) continue;
+	// 		if(av[k].v == w) break;
+	// 	}
+	// 	if(k >= nv) return 0;
+
+	// 	v = en^1; w = pa[pn-1]^1;
+	// 	nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+	// 	for (k = 0; k < nv; k++) {
+	// 		if(av[k].del) continue;
+	// 		if(av[k].v == w) break;
+	// 	}
+	// 	if(k >= nv) return 0;
+	// }
+
+	if(min_ou) (*min_ou) = ou;
+	return 1;
+}
+
+void append_path_connec(asg_t *g, uint64_t sn, uint64_t en, uint64_t *pa, uint64_t pn, uint64_t min_ou)
+{
+	uint64_t nv, k, v, w, i, ou; asg_arc_t *av;
+    if(sn != (uint64_t)-1) {
+        v = sn; i = 0;
+    } else {
+        v = pa[0]; i = 1; 
+    }
+
+    for (; i < pn; i++) {
+        w = pa[i]; 
+
+        nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+        for (k = 0; k < nv; k++) {
+            if(av[k].del) continue;
+            if(av[k].v == w) {
+				ou = av[k].ou; ou += min_ou; 
+				if(ou > OU_MASK) ou = OU_MASK;
+				av[k].ou = ou;
+                break;
+            }
+        }
+		assert(k < nv);
+
+        nv = asg_arc_n(g, w^1); av = asg_arc_a(g, w^1);
+        for (k = 0; k < nv; k++) {
+            if(av[k].del) continue;
+            if(av[k].v == (v^1)) {
+                ou = av[k].ou; ou += min_ou; 
+				if(ou > OU_MASK) ou = OU_MASK;
+				av[k].ou = ou;
+                break;
+            }
+        }
+		assert(k < nv);
+
+        v = w;
+    }
+    // fprintf(stderr, "-3-[M::%s::]\n", __func__); 
+    if(en != (uint64_t)-1) {
+        w = en;
+        // fprintf(stderr, "-3a-[M::%s::]\n", __func__); 
+        nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+        for (k = 0; k < nv; k++) {
+            if(av[k].del) continue;
+            if(av[k].v == w) {
+                ou = av[k].ou; ou += min_ou; 
+				if(ou > OU_MASK) ou = OU_MASK;
+				av[k].ou = ou;
+                break;
+            }
+        }
+		assert(k < nv);
+
+        // fprintf(stderr, "-3b-[M::%s::]\n", __func__); 
+        nv = asg_arc_n(g, w^1); av = asg_arc_a(g, w^1);
+        for (k = 0; k < nv; k++) {
+            if(av[k].del) continue;
+            if(av[k].v == (v^1)) {
+                ou = av[k].ou; ou += min_ou; 
+				if(ou > OU_MASK) ou = OU_MASK;
+				av[k].ou = ou;
+                break;
+            }
+        }
+		assert(k < nv);
+    }
+}
+
+
+uint64_t asg_clean_idx_tig(asg_t *g, uint64_t *idx, uint64_t idx_n, uint64_t *srt, uint64_t srt_n, utepdat_t *s)
+{
+	uint64_t i, k, *a, an, sp, *del_a, del_n, *ref_a, ref_n, *ka, kk, cnt = 0, min_ou;
+	uint64_t ref_s, ref_e, del_s, del_e;
+	for (i = 0; i < srt_n; i++) {
+		a = s->ll[idx[(uint32_t)srt[i]]>>32].srt.a.a; 
+		an = s->ll[idx[(uint32_t)srt[i]]>>32].srt.a.n;
+		sp = (uint32_t)idx[(uint32_t)srt[i]];
+		assert((((uint32_t)a[sp])==((uint32_t)-1)) && ((a[sp]>>32)!=((uint32_t)-1)));
+
+		// fprintf(stderr, "\n+[M::%s::] ((uint32_t)(a[sp]>>32))::%u, ((uint32_t)a[sp])::%u, sp::%lu, an::%lu\n", __func__, 
+		// ((uint32_t)(a[sp]>>32)), ((uint32_t)a[sp]), sp, an); 
+
+		k = sp + 1; ref_a = a + k;
+		for (sp = k; k < an; k++) {
+			if(a[k] == (uint64_t)-1) break;
+		}
+		assert(k < an);
+		ref_n = k - sp;
+		// fprintf(stderr, "+[M::%s::] k::%lu\n", __func__, k); 
+		// fprintf(stderr, "+[M::%s::] k::%lu, an::%lu\n", __func__, k, an); 
+
+		k++; del_a = a + k;
+		for (sp = k; k < an; k++) {
+			if((a[k] != (uint64_t)-1) && (((uint32_t)a[k]) == ((uint32_t)-1))) break;
+		}
+		
+		// assert(k < an);//for the last batch of reads, there is no end marker, so k == an
+		del_n = k - sp;
+
+		ref_s = ref_e = del_s = del_e = (uint64_t)-1;
+		if(ref_n > 0 && (ref_a[0]>>32)) ref_s = (uint32_t)ref_a[0];
+		if(ref_n > 1 && (ref_a[ref_n-1]>>32)) ref_e = (uint32_t)ref_a[ref_n-1];
+
+		if(del_n > 0 && (del_a[0]>>32)) del_s = (uint32_t)del_a[0];
+		if(del_n > 1 && (del_a[del_n-1]>>32)) del_e = (uint32_t)del_a[del_n-1];
+
+		// if((ref_e != del_e) || (ref_s != del_s)) {
+		// 	fprintf(stderr, "-[M::%s::] idx_n::%lu, idx_n::%lu, ref_n::%lu, del_n::%lu, ref_s::%lu, ref_e::%lu, del_s::%lu, del_e::%lu\n", 
+		// 	__func__, idx_n, idx_n, ref_n, del_n, ref_s, ref_e, del_s, del_e); 
+		// 	for (k = 0; k < ref_n; k++) {
+		// 		fprintf(stderr, "-[M::ref::] ref_a[k]>>1::%lu\n", ref_a[k]>>1); 
+		// 	}
+		// 	for (k = 0; k < del_n; k++) {
+		// 		fprintf(stderr, "-[M::del::] del_a[k]>>1::%lu\n", del_a[k]>>1); 
+		// 	}
+		// }
+		assert(ref_s == del_s); assert(ref_e == del_e); 
+
+		if(ref_s != (uint64_t)-1) {
+			ref_a = ref_a + 1; ref_n--;
+		}
+		if(ref_e != (uint64_t)-1) ref_n--;
+
+		if(del_s != (uint64_t)-1) {
+			del_a = del_a + 1; del_n--;
+		}
+		if(del_e != (uint64_t)-1) del_n--;
+
+		if(ref_n == 0 && del_n > 0) {
+			ka = ref_a; ref_a = del_a; del_a = ka;
+			kk = ref_n; ref_n = del_n; del_n = kk;
+			kk = ref_s; ref_s = del_s; del_s = kk;
+			kk = ref_e; ref_e = del_e; del_e = kk;
+		}
+
+		// fprintf(stderr, "-[M::%s::] idx_n::%lu, idx_n::%lu, ref_n::%lu, del_n::%lu, ref_s::%lu, ref_e::%lu\n", 
+		// __func__, idx_n, idx_n, ref_n, del_n, ref_s, ref_e); 
+		// for (k = 0; k < ref_n; k++) {
+		// 	fprintf(stderr, "-[M::ref::] ref_a[k]>>1::%lu\n", ref_a[k]>>1); 
+		// }
+		// for (k = 0; k < del_n; k++) {
+		// 	fprintf(stderr, "-[M::del::] del_a[k]>>1::%lu\n", del_a[k]>>1); 
+		// }
+		
+
+
+		///just for debug
+		// if(ref_s != (uint64_t)-1) {
+		// 	assert(get_arcs(g, ref_s, NULL, 0) > 1);
+		// }
+		// if(ref_e != (uint64_t)-1) {
+		// 	assert(get_arcs(g, ref_e^1, NULL, 0) > 1);
+		// }
+		// assert(tes_path_connec(g, ref_s, ref_e, ref_a, ref_n, 0, NULL));
+		// assert(tes_path_connec(g, del_s, del_e, del_a, del_n, 1, NULL));
+		
+		min_ou = 0;
+		if((ref_s!=(uint64_t)-1) && (get_arcs(g, ref_s, NULL, 0)<=1)) continue;
+		if((ref_e!=(uint64_t)-1) && (get_arcs(g, ref_e^1, NULL, 0)<=1)) continue;
+		if(!tes_path_connec(g, ref_s, ref_e, ref_a, ref_n, 0, NULL)) continue;
+		if(!tes_path_connec(g, del_s, del_e, del_a, del_n, 1, &min_ou)) continue;
+
+		if(del_n > 0) {
+			for (k = 0; k < del_n; k++) asg_seq_del(g, del_a[k]>>1);
+			cnt += del_n;
+		} else {
+			asg_arc_del(g, del_s, del_e, 1);
+			asg_arc_del(g, del_e^1, del_s^1, 1);
+			cnt++;
+		}
+		if(min_ou > 0) append_path_connec(g, ref_s, ref_e, ref_a, ref_n, min_ou);
+	}
+	if(cnt) {
+		asg_cleanup(g);
+	}
+	return cnt;
+}
+
+uint64_t work_rg_contain_dedup(uldat_t *sl)
+{
+	utepdat_t s; uint64_t i, occ, ou_n, m, *srt, srt_n, *idx, idx_n, tot; 
+	asg64_v ou, bu; kv_init(ou); memset(&s, 0, sizeof(s));
+	s.id = 0; s.uopt = sl->uopt; s.rg = sl->rg; s.opt = sl->opt; ///s.ug = sl->ug; s.uu = sl->uu; 
+	//  CALLOC(s.buf, sl->n_thread); CALLOC(s.gdp, sl->n_thread); CALLOC(s.mzs, sl->n_thread); 
+	CALLOC(s.hab, sl->n_thread); CALLOC(s.ll, sl->n_thread); ///CALLOC(s.sps, sl->n_thread); 
+
+	for (i = 0; i < sl->n_thread; ++i) {
+		s.hab[i] = ha_ovec_init(0, 0, 1); ///s.buf[i] = mg_tbuf_init();
+	}
+
+	// detect_outlier_len("+++work_ul_gchains");
+
+	//debug
+    // overall_zdbg = init_mul_debug_prt_t(UL_INF.n);
+	occ = 1; tot = 0;
+	while (occ) {
+		kt_for(sl->n_thread, worker_for_contain_dedup, &s, (s.rg->n_seq<<1));
+		for (i = ou.n = 0; i < sl->n_thread; ++i) {
+			copy_asg_arr(bu, s.ll[i].srt.a);
+			update_contain_dedup_arr(&bu, i, &ou);
+			copy_asg_arr(s.ll[i].srt.a, bu);
+		}
+		ou_n = ou.n;
+		for (i = 0; i < ou_n; i++) {
+			assert(((uint32_t)s.ll[ou.a[i]>>32].srt.a.a[(uint32_t)ou.a[i]]) == ((uint32_t)-1));
+			m = s.ll[ou.a[i]>>32].srt.a.a[(uint32_t)ou.a[i]]>>32; m <<= 32; m += i;
+			kv_push(uint64_t, ou, m);
+		}
+		// fprintf(stderr, "[M::%s::] ou_n::%lu, ou.n::%u\n", __func__, ou_n, (uint32_t)ou.n); 
+		idx = ou.a; idx_n = ou_n; srt = ou.a + ou_n; srt_n = ou.n - ou_n;
+		radix_sort_gfa64(srt, srt + srt_n);
+		// print_debug_gfa((asg_t *)s.rg, NULL, s.uopt->coverage_cut, "UL.dirty.debug", s.uopt->sources, s.uopt->ruIndex, s.uopt->max_hang, s.uopt->min_ovlp, 0, 0, 0);
+		occ = asg_clean_idx_tig((asg_t *)s.rg, idx, idx_n, srt, srt_n, &s);
+		for (i = 0; i < sl->n_thread; ++i) s.ll[i].srt.a.n = 0;
+		tot += occ;
+		// fprintf(stderr, "[M::%s::] idx_n::%lu, srt_n::%lu, occ::%lu\n", __func__, idx_n, srt_n, occ); 
+	}
+	
+
+	
+
+	// detect_outlier_len("---work_ul_gchains");
+
+	for (i = 0; i < sl->n_thread; ++i) {
+		s.sum_len += s.hab[i]->num_read_base; s.n += s.hab[i]->num_correct_base;
+		// hc_gdpchain_destroy(&(s.gdp[i])); kv_destroy(s.mzs[i]); mg_tbuf_destroy(s.buf[i]); 
+		ha_ovec_destroy(s.hab[i]); hc_glchain_destroy(&(s.ll[i])); ///kv_destroy(s.sps[i]);
+	}
+
+	// free(s.buf); free(s.gdp); free(s.mzs); 
+	free(s.hab); free(s.ll); kv_destroy(ou); ///free(s.sps);
+	fprintf(stderr, "[M::%s::] # duplicated reads::%lu\n", __func__, tot); 
+	// exit(1);
+	return tot;
+}
+
 void print_ul_ovlps(all_ul_t *x, int32_t prt_ovlp)
 {
 	uint64_t k, i, ucov_occ = 0, cov_occ = 0, ucov_len = 0, cov_len = 0, unaligned_len = 0, unaligned_occ = 0, aligned_occ = 0;
@@ -15847,7 +17548,7 @@ void determine_connective_adv(all_ul_t *m, const ug_opt_t *uopt, int64_t bw, dou
 	}
 }
 
-void determine_connective_backtrack(all_ul_t *m, const ug_opt_t *uopt, ul_vec_t *p, uint32_t ii, uint64_t rid)
+void determine_connective_backtrack(all_ul_t *m, const ug_opt_t *uopt, ul_vec_t *p, uint32_t ii, uint64_t rid, uint64_t ulid)
 {
 	assert((!p->bb.a[ii].base)&&(p->bb.a[ii].hid == rid)&&(p->bb.a[ii].el));
 	if(ii <= 0) return;
@@ -15876,7 +17577,13 @@ void determine_connective_backtrack(all_ul_t *m, const ug_opt_t *uopt, ul_vec_t
 					break;
 				}
 			}
-			if(z < x->length) x->buffer[z].bl++;
+			if(z < x->length) {
+				x->buffer[z].bl++;
+				// if(((li_v>>1) == 1656459) || ((li_v>>1) == 4179430)|| 
+				// 	((lk_v>>1) == 1656459) || ((lk_v>>1) == 4179430)) {
+				// 	fprintf(stderr, "[M::%s::]\tulid::%lu\n", __func__, ulid);
+				// }
+			}
 		}
 
 		lk = ((lk->pidx==(uint32_t)-1)?NULL:&(p->bb.a[lk->pidx]));
@@ -15895,7 +17602,7 @@ static void update_ovlp_src(void *data, long i, int tid) // callback for kt_for(
 	for (k = 0; k < a_n; k++) {
 		///note: we only label reliable chains
 		// if(i == 4217) fprintf(stderr, "\nul_id: %lu\n", a[k]>>32);
-		determine_connective_backtrack(&UL_INF, sl->uopt, &(UL_INF.a[a[k]>>32]), (uint32_t)(a[k]), i);
+		determine_connective_backtrack(&UL_INF, sl->uopt, &(UL_INF.a[a[k]>>32]), (uint32_t)(a[k]), i, a[k]>>32);
 		// determine_connective_adv(&UL_INF, sl->uopt, G_CHAIN_BW, sl->opt->diff_ec_ul, &(UL_INF.a[a[k]>>32]), (uint32_t)(a[k]), i);
 		// determine_connective(&UL_INF, sl->uopt, G_CHAIN_BW, sl->opt->diff_ec_ul, 
 		// 										&(UL_INF.a[a[k]>>32]), (uint32_t)(a[k]), i);
@@ -15904,6 +17611,38 @@ static void update_ovlp_src(void *data, long i, int tid) // callback for kt_for(
 }
 
 
+static void clean_contained_chg(void *data, long i, int tid) // callback for kt_for()
+{
+	uldat_t *sl = (uldat_t *)data;
+	sl->rg->seq_vis[i] = 0;
+	if(sl->rg->seq[i].del) return;
+	if(!(is_contain_r((*(sl->uopt->ruIndex)), ((uint64_t)i)))) return;
+	ma_hit_t_alloc* src = sl->uopt->sources;
+	uint64_t z;
+	for (z = 0; z < src[i].length; z++) {
+		if(src[i].buffer[z].bl) break;
+	}
+	if(z >= src[i].length) sl->rg->seq_vis[i] = 1;
+}
+
+static void label_contained_chg(void *data, long i, int tid) // callback for kt_for()
+{
+	uldat_t *sl = (uldat_t *)data; uint32_t v = i;
+	sl->rg->seq_vis[v] = 0;
+	if(sl->rg->seq[v>>1].del) return;
+	if(is_contain_r((*(sl->uopt->ruIndex)), (v>>1))) {
+		sl->rg->seq_vis[v] = 1; return;
+	}
+	asg_arc_t *av = asg_arc_a(sl->rg, v);
+	uint32_t nv = asg_arc_n(sl->rg, v), k;
+	for (k = 0; k < nv; ++k) {
+		if(av[k].del) continue;
+		if(is_contain_r((*(sl->uopt->ruIndex)), (av[k].v>>1))) break;
+	}
+	if(k < nv) sl->rg->seq_vis[v] = 1;
+}
+
+
 uint64_t* get_hifi2ul_list(all_ul_t *x, uint64_t hid, uint64_t* a_n)
 {
 	(*a_n) = x->ridx.idx.a[hid+1] - x->ridx.idx.a[hid];
@@ -17599,6 +19338,173 @@ void ul_load(const ug_opt_t *uopt)
 	// destory_all_ul_t(&UL_INF); 
 }
 
+void clean_contain_g0(uldat_t *sl)
+{
+	kt_for(sl->n_thread, clean_contained_chg, sl, R_INF.total_reads);
+	asg_t *sg = (asg_t *)sl->rg; uint32_t k, cnt;
+	for (k = cnt = 0; k < sg->n_seq; k++) {
+		// if(is_contain_r((*(sl->uopt->ruIndex)), k)) {
+		// 	fprintf(stderr, "[M::%s::]\t%.*s\n", __func__, (int)Get_NAME_LENGTH(R_INF, k), Get_NAME(R_INF, k));
+		// }
+		if(sg->seq_vis[k]) {
+			asg_seq_del(sg, k); cnt++;
+		}
+		sg->seq_vis[k] = 0; 
+	}
+	if(cnt) asg_cleanup(sg);
+	fprintf(stderr, "[M::%s::] # discard cread::%u\n", __func__, cnt); 
+	// exit(1);
+}
+
+void asg_arc_push_contain_trans(asg_t *g, ma_hit_t_alloc* ov, int64_t min_ovlp, int64_t max_hang, double max_hang_rate, int64_t fuzz, R_to_U *ri, uint8_t *mark)
+{
+	uint32_t n_vtx = g->n_seq<<1, cnt = 0, nc, cc, qn, tn, avi, awi;
+	uint32_t v, w, i, k, nv, nw; asg_arc_t *av, *aw; int32_t r = 1;
+	ma_hit_t_alloc *z; asg_arc_t p, *t;
+	uint32_t *dis, *idx; CALLOC(dis, n_vtx); MALLOC(idx, n_vtx);
+	for (v = 0; v < n_vtx; ++v) {
+		if(g->seq[v>>1].del || (!mark[v])) continue;
+        nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
+        if (nv == 0) continue; // no hits    
+		for (i = 0; i < nv; ++i) {
+			if((!(av[i].del)) && (is_contain_r((*ri), (av[i].v>>1)))) break;
+		}
+		if(i >= nv) continue; 
+		z = &(ov[v>>1]);
+		for (i = 0; i < z->length; i++) {
+        	qn = Get_qn(z->buffer[i]); tn = Get_tn(z->buffer[i]);
+			if(z->buffer[i].del) continue;
+			r = ma_hit2arc(&(z->buffer[i]), g->seq[qn].len, g->seq[tn].len, max_hang, max_hang_rate, min_ovlp, &p);
+			if((r < 0) || ((p.ul>>32) != v)) continue;
+			if(g->seq[p.v>>1].del) continue;
+			if(mark[p.v^1]) {
+				dis[p.v] = asg_arc_len(p) + fuzz; idx[p.v] = i;
+			}
+		}
+
+		nv = asg_arc_n(g, v); av = asg_arc_a(g, v); avi = g->idx[v]>>32;
+		while(nv) {
+			for (i = nc = cc = 0; i < nv; ++i) {
+				if(av[i].del) continue; w = av[i].v;
+				if(!(is_contain_r((*ri), (w>>1)))) continue;///new arcs must be bridged by contained reads
+				assert(!(g->seq[w>>1].del));
+				nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); awi = g->idx[w]>>32;
+				for (k = 0; k < nw; k++) {
+					if(aw[k].del) continue;
+					if((dis[aw[k].v] == ((uint32_t)-1)) || (dis[aw[k].v] == 0)) continue;
+					assert(!(g->seq[aw[k].v>>1].del));
+					if((asg_arc_len(av[i]) + asg_arc_len(aw[k])) <= dis[aw[k].v]) {
+						dis[aw[k].v] = ((uint32_t)-1); assert(!(z->buffer[idx[aw[k].v]].del));
+						qn = Get_qn(z->buffer[idx[aw[k].v]]); tn = Get_tn(z->buffer[idx[aw[k].v]]);
+						r = ma_hit2arc(&(z->buffer[idx[aw[k].v]]), g->seq[qn].len, g->seq[tn].len, max_hang, max_hang_rate, min_ovlp, &p);
+						assert(r>=0); assert((p.ul>>32)==v); assert(p.v==aw[k].v);
+						cnt++; p.ou = 0; t = asg_arc_pushp(g); *t = p; 
+
+						// fprintf(stderr, "[M::%s::]\t%.*s(id::%lu::%c)(is_c::%u)\t%.*s(id::%u::%c)(is_c::%u)\n", __func__, 
+						// (int)Get_NAME_LENGTH(R_INF, (g->arc[g->n_arc-1].ul>>33)), 
+						// Get_NAME(R_INF, (g->arc[g->n_arc-1].ul>>33)), 
+						// g->arc[g->n_arc-1].ul>>33, "+-"[(g->arc[g->n_arc-1].ul>>32)&1], 
+						// (is_contain_r((*ri), (g->arc[g->n_arc-1].ul>>33))),
+						// (int)Get_NAME_LENGTH(R_INF, (g->arc[g->n_arc-1].v>>1)),
+						// Get_NAME(R_INF, (g->arc[g->n_arc-1].v>>1)), 
+						// g->arc[g->n_arc-1].v>>1, "+-"[(g->arc[g->n_arc-1].v)&1],
+						// (is_contain_r((*ri), (g->arc[g->n_arc-1].v>>1)))
+						// );
+						// fprintf(stderr, "[M::%s::]\tmiddle::%.*s(id::%u::%c)(is_c::%u)\n", __func__, 
+						// (int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), 
+						// w>>1, "+-"[w&1], (is_contain_r((*ri), (w>>1))));
+
+
+
+						av = g->arc + avi; aw = g->arc + awi;///renew av and aw since asg_arc_pushp
+						if((is_contain_r((*ri), (g->arc[g->n_arc-1].v>>1)))) {
+							cc++;
+						} else {
+							if(g->n_arc!=(nc+1)) {		
+								p = g->arc[g->n_arc-1]; 						
+								g->arc[g->n_arc-1] = g->arc[nc];
+								g->arc[nc] = p;
+							}
+							nc++;
+						}
+					}
+				}
+			}
+
+			if(cc) {///new arcs must be bridged by contained reads
+				nv = cc; av = g->arc + g->n_arc - cc; avi = g->n_arc - cc;
+			} else {
+				nv = 0; av = NULL; avi = ((uint32_t)-1);
+			}
+		}
+
+		z = &(ov[v>>1]);
+		for (i = 0; i < z->length; i++) {
+        	qn = Get_qn(z->buffer[i]); tn = Get_tn(z->buffer[i]);
+			if(z->buffer[i].del) continue;
+			r = ma_hit2arc(&(z->buffer[i]), g->seq[qn].len, g->seq[tn].len, max_hang, max_hang_rate, min_ovlp, &p);
+			if((r < 0) || ((p.ul>>32) != v)) continue;
+			dis[p.v] = 0;
+		}
+	}
+
+	free(dis); free(idx);
+	if(cnt) {
+		free(g->idx);
+        g->idx = 0;
+        g->is_srt = 0;
+        asg_cleanup(g);
+		asg_symm(g);
+	}
+}
+
+void repush_contain_trans_archs(uldat_t *sl)
+{
+	asg_t *sg = (asg_t *)sl->rg;
+	kt_for(sl->n_thread, label_contained_chg, sl, (sg->n_seq<<1));
+
+	// print_debug_gfa(sg, NULL, sl->uopt->coverage_cut, "UL.dirty.debug0", sl->uopt->sources, sl->uopt->ruIndex, sl->uopt->max_hang, sl->uopt->min_ovlp, 0, 0, 0);
+	
+	asg_arc_push_contain_trans(sg, sl->uopt->sources, sl->uopt->min_ovlp, sl->uopt->max_hang, asm_opt.max_hang_rate, sl->uopt->gap_fuzz, sl->uopt->ruIndex, sg->seq_vis);
+
+	// print_debug_gfa(sg, NULL, sl->uopt->coverage_cut, "UL.dirty.debug1", sl->uopt->sources, sl->uopt->ruIndex, sl->uopt->max_hang, sl->uopt->min_ovlp, 0, 0, 0);
+	// exit(1);
+}
+
+uint32_t clean_contain_g(const ug_opt_t *uopt, asg_t *sg, uint32_t push_trans)
+{
+	mg_idxopt_t opt; uldat_t sl; int32_t cutoff, f = 0;
+	init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
+	cutoff = asm_opt.max_n_chain;
+	init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round);
+	init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, NULL); sl.rg = sg;
+	
+	if(push_trans) repush_contain_trans_archs(&sl);
+	if(work_ul_gchains_consensus(&sl)) {
+		free(UL_INF.ridx.idx.a); free(UL_INF.ridx.occ.a); memset(&(UL_INF.ridx), 0, sizeof(UL_INF.ridx));
+		gen_ul_vec_rid_t(&UL_INF, &R_INF, NULL);
+		// fprintf(stderr, "+[M::%s::]\n", __func__);
+		kt_for(sl.n_thread, update_ovlp_src, &sl, R_INF.total_reads);
+		// fprintf(stderr, "-[M::%s::]\n", __func__);
+		kt_for(sl.n_thread, update_ovlp_src_bl, &sl, R_INF.total_reads);
+		clean_contain_g0(&sl);
+		f = 1;
+	}
+	if(push_trans) asg_arc_del_trans_ul(sg, sl.uopt->gap_fuzz);
+	return f;
+}
+
+
+void dedup_contain_g(const ug_opt_t *uopt, asg_t *sg)
+{
+	mg_idxopt_t opt; uldat_t sl; int32_t cutoff;
+	init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
+	cutoff = asm_opt.max_n_chain;
+	init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round);
+	init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, NULL); sl.rg = sg;
+	work_rg_contain_dedup(&sl);
+}
+
 
 uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg)
 {
diff --git a/inter.h b/inter.h
index 1732391..e5b74b5 100644
--- a/inter.h
+++ b/inter.h
@@ -124,5 +124,7 @@ uint32_t gen_src_shared_interval_simple(uint32_t src, ma_ug_t *ug, uint64_t *flt
 uint64_t check_ul_ov_t_consist(ul_ov_t *x, ul_ov_t *y, int64_t ql, int64_t tl, double diff);
 uint32_t infer_se(uint32_t qs, uint32_t qe, uint32_t ts, uint32_t te, uint32_t rev, 
 uint32_t rqs, uint32_t rqe, uint32_t *rts, uint32_t *rte);
+uint32_t clean_contain_g(const ug_opt_t *uopt, asg_t *sg, uint32_t push_trans);
+void dedup_contain_g(const ug_opt_t *uopt, asg_t *sg);
 
 #endif