diff --git a/Assembly.cpp b/Assembly.cpp
index f7f65ab..26f3c01 100644
--- a/Assembly.cpp
+++ b/Assembly.cpp
@@ -10,6 +10,7 @@
 #include "Correct.h"
 #include "htab.h"
 #include "kthread.h"
+#include "rcut.h"
 
 void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
 int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct);
@@ -1605,6 +1606,7 @@ void ha_overlap_final(void)
 
 int ha_assemble(void)
 {
+    // debug_mc_g_t(MC_NAME);
 	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
 	int r, hom_cov = -1, ovlp_loaded = 0;
 	if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
diff --git a/Overlaps.cpp b/Overlaps.cpp
index 86220fb..aa5f028 100644
--- a/Overlaps.cpp
+++ b/Overlaps.cpp
@@ -12736,7 +12736,7 @@ bub_label_t* b_mask_t)
         ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic;
         adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, 
         tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-        max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1);
+        max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1/**0**/);
 
         print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, 
         min_ovlp, &new_rtg_edges);
@@ -13466,14 +13466,12 @@ bub_label_t* b_mask_t)
 
     hap_cov_t *cov = NULL;
     asg_t *copy_sg = copy_read_graph(sg);
-    ma_ug_t *copy_ug = copy_untig_graph(ug);    
+    ma_ug_t *copy_ug = copy_untig_graph(ug);
     adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, 
     tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-    max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1);
-
+    max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1);
     print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, 
     min_ovlp, &new_rtg_edges);
-
     ma_ug_destroy(copy_ug);
     asg_destroy(copy_sg);
 
@@ -15318,7 +15316,6 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
     {    
         pre_cons = get_graph_statistic(g);
         asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, 1);
-
         if(just_bubble_pop == 0)
         {
             ///need consider tangles  
@@ -16088,7 +16085,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
     {
 		purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, 
         asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 1, 0, 
-        cov, 0);
+        cov, 0, 0);
 		///delete_useless_nodes(ug);
         delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
 	}
@@ -16388,6 +16385,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
     {
         b->b.n--;
         cur = b->b.a[b->b.n];
+        if(flag[cur>>1] == 0 && (cur>>1) != (v>>1)) continue;
         flag[cur>>1] = 0;
 
         ncur = asg_arc_n(g, cur);
@@ -16396,6 +16394,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
         {
             if(acur[i].del) continue;
             if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue;
+            if(flag[acur[i].v>>1] == 0) continue;
             kv_push(uint32_t, b->b, acur[i].v);
         }
     }
@@ -16406,6 +16405,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
     {
         b->b.n--;
         cur = b->b.a[b->b.n];
+        if(flag[cur>>1] == 0 && (cur>>1) != (v>>1)) continue;
         flag[cur>>1] = 0;
 
         ncur = asg_arc_n(g, cur);
@@ -16414,6 +16414,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
         {
             if(acur[i].del) continue;
             if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue;
+            if(flag[acur[i].v>>1] == 0) continue;
             kv_push(uint32_t, b->b, acur[i].v);
         }
     }
@@ -16641,14 +16642,14 @@ void chain_origin_trans_uid_c_bubble(uint32_t query, buf_t *target, buf_t *idx,
 // in a resolved bubble, mark unused vertices and arcs as "reduced"
 static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, hap_cov_t *cov, uint32_t is_update_chain)
 {
-	uint32_t i, k, k_i, v, u, uLen = 0, uCov = 0, uId, rId, ori;
+    uint32_t i, k, k_i, v, u, uLen = 0, uCov = 0, uId, rId, ori;
     ma_utg_t* p = NULL;
     trans_chain* t_ch = (is_update_chain?cov->t_ch:NULL);
     ///b->S.a[0] is the sink of this bubble
 
-	///assert(b->S.n == 1);
-	///first remove all nodes in this bubble
-	for (i = 0; i < b->b.n; ++i)
+    ///assert(b->S.n == 1);
+    ///first remove all nodes in this bubble
+    for (i = 0; i < b->b.n; ++i)
     {
         uId = b->b.a[i]>>1;
         if(uId == (b->S.a[0]>>1)) continue;
@@ -16664,10 +16665,10 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
 
     
     ///v is the sink of this bubble
-	v = b->S.a[0];
-	///recover node
-	do {
-		u = b->a[v].p; // u->v
+    v = b->S.a[0];
+    ///recover node
+    do {
+        u = b->a[v].p; // u->v
         if(v != b->S.a[0])
         {
             uId = v>>1;
@@ -16680,16 +16681,16 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
                 uLen += cov->read_g->seq[rId].len;
             }
         }
-		v = u;
-	} while (v != v0);
+        v = u;
+    } while (v != v0);
 
     uCov = (uLen == 0? 0 : uCov / uLen);
 
     ///v is the sink of this bubble
-	v = b->S.a[0];
-	///recover node
-	do {
-		u = b->a[v].p; // u->v
+    v = b->S.a[0];
+    ///recover node
+    do {
+        u = b->a[v].p; // u->v
         if(v != b->S.a[0])
         {
             uId = v>>1;
@@ -16701,8 +16702,8 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
                 cov->cov[rId] += (uCov * cov->read_g->seq[rId].len);
             }
         }
-		v = u;
-	} while (v != v0);
+        v = u;
+    } while (v != v0);
 
 
     if(t_ch)
@@ -16763,7 +16764,6 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
             }
         }
 
-
         topologicalSortUtil(ug->g, cov, v0, b->S.a[0]);
         ///if(cov->t_ch->topo_res.n != b->b.n - 1) fprintf(stderr, "ERROR-4\n");
         if(cov->t_ch->topo_res.n == 0) return;
@@ -16773,8 +16773,10 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
             if(uId == (b->S.a[0]>>1)) continue;
 
             dfs_trans_chain_bub(ug->g, cov, uId, v0>>1, b->S.a[0]>>1);
+
             if(cov->t_ch->b_buf_0.b.n == 0) continue;
             chain_origin_trans_uid_c_bubble(cov->t_ch->topo_res.a[i], &(t_ch->b_buf_0), b, ug, cov);
+
             /***********************x***********************/
             uId = cov->t_ch->topo_res.a[i]>>1;
             p = &(ug->u.a[uId]);
@@ -17528,7 +17530,6 @@ hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d)
     /****************************may have bugs********************************/
     ///if(keep_d != 0) debug_asg_bub_pop1_primary_trio(g, utg, v0, max_dist, b, positive_flag, negative_flag, 1);
     /****************************may have bugs********************************/
-
     if(cov && utg) asg_bub_backtrack_primary_cov(utg, v0, b, cov, is_update_chain);
 	if(is_pop) asg_bub_backtrack_primary(g, v0, b);
     if(path_base_len || path_nodes) asg_bub_backtrack_primary_length(g, utg, v0, b, path_base_len, path_nodes);
@@ -23678,7 +23679,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
 R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
 kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, 
-uint32_t collect_p_trans)
+uint32_t collect_p_trans, uint32_t collect_p_trans_f)
 {
     asg_t* nsg = (*ug)->g;
     uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain;
@@ -23686,9 +23687,8 @@ uint32_t collect_p_trans)
     hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, 
                             coverage_cut, max_hang, min_ovlp, (asm_opt.purge_level_primary>0||i_cov)?1:0);
     if(cov->t_ch) set_r_het_flag(*ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
-    
     adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t);
-    
+
     nsg = (*ug)->g;
     n_vtx = nsg->n_seq;
     for (v = 0; v < n_vtx; ++v) 
@@ -23698,14 +23698,11 @@ uint32_t collect_p_trans)
         EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n;
     }
 
-
     clean_primary_untig_graph(*ug, read_g, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, 
     stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, cov);
     delete_useless_nodes(ug);
     renew_utg(ug, read_g, new_rtg_edges);
-    
     if(i_cov && collect_p_trans == 0) goto skip_purge;
-
     if(asm_opt.purge_level_primary > 0)
     {
         ///print_debug_gfa(read_g, *ug, coverage_cut, "debug_purge", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
@@ -23713,7 +23710,7 @@ uint32_t collect_p_trans)
         if(asm_opt.purge_level_primary == 1) just_contain = 1;
         purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, 
         asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 
-        just_contain, 0, cov, !!(cov->t_ch&&collect_p_trans));
+        just_contain, 0, cov, !!(cov->t_ch&&collect_p_trans), collect_p_trans_f);
         delete_useless_nodes(ug);
         renew_utg(ug, read_g, new_rtg_edges);
     }
@@ -23855,26 +23852,21 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov
     kvec_asg_arc_t_warp new_rtg_edges;
     kv_init(new_rtg_edges.a);
 
-
     adjust_utg_by_primary(&ug, sg, TRIO_THRES, sources, reverse_sources, coverage_cut, 
     tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-    max_hang, min_ovlp, &new_rtg_edges, NULL, b_mask_t, 0);
-
+    max_hang, min_ovlp, &new_rtg_edges, NULL, b_mask_t, 0, 0);
 
     if(asm_opt.b_low_cov > 0)
     {
         break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp, 
         &asm_opt.b_low_cov, NULL, asm_opt.m_rate);
     }
-
     if(asm_opt.b_high_cov > 0)
     {
         break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp, 
         NULL, &asm_opt.b_high_cov, asm_opt.m_rate);
     }
-
     ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0);
-
     
     
     fprintf(stderr, "Writing primary contig GFA to disk... \n");
@@ -27096,8 +27088,8 @@ void reset_bub(bubble_type* bub, ma_ug_t *ug, trans_chain* back_ug_chain, kvec_a
     new_rtg_edges->a.n = 0;
     ///classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp);
     identify_bubbles(ug, bub, back_ug_chain->is_r_het, NULL);
-    update_bubble_chain(ug, bub, 0, 1);
-    resolve_bubble_chain_tangle(ug, bub);
+    // update_bubble_chain(ug, bub, 0, 1);
+    // resolve_bubble_chain_tangle(ug, bub);
     // fprintf(stderr, "bub.f_bub: %lu, bub.b_bub: %lu, bub.b_end_bub: %lu, bub.tangle_bub: %lu, bub.cross_bub: %lu\n", 
     // bub->f_bub, bub->b_bub, bub->b_end_bub, bub->tangle_bub, bub->cross_bub);
 }
@@ -27396,16 +27388,20 @@ bub_label_t* b_mask_t)
     ma_ug_t *ug = NULL;
     ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
 
+    // FILE* output_file = fopen("straw-debug.noseq.gfa", "w");
+    // ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file);
+    // fclose(output_file);
+
+
     hap_cov_t *cov = NULL;
     asg_t *copy_sg = copy_read_graph(sg);
     ma_ug_t *copy_ug = copy_untig_graph(ug);    
     adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, 
     tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-    max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 0);
+    max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 0, 0);
     ma_ug_destroy(copy_ug); copy_ug = NULL;
     asg_destroy(copy_sg); copy_sg = NULL;
 
-
     uint32_t beg_idx, occ;
     bubble_type bub; 
     memset(&bub, 0, sizeof(bubble_type));
@@ -29276,7 +29272,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int m
 
         hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0);
         purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, NULL, 
-        opt->purge_simi_thres, opt->purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0);
+        opt->purge_simi_thres, opt->purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0, 0);
         destory_hap_cov_t(&cov);
 
         ma_ug_destroy(ug);
@@ -29562,7 +29558,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
         init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq);
         goto debug_gfa;
     }
-    
 
     ///just for debug
     renew_graph_init(sources, reverse_sources, sg, coverage_cut, ruIndex, n_read);
@@ -29571,7 +29566,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     ///normalize_ma_hit_t_single_side(sources, n_read);
     normalize_ma_hit_t_single_side_advance(sources, n_read);
     normalize_ma_hit_t_single_side_advance(reverse_sources, n_read);
-
     if (ha_opt_triobin(&asm_opt))
     {
         drop_edges_by_trio(sources, n_read);
@@ -29580,7 +29574,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     {
         memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
     }
-    
     ///print_binned_reads(sources, n_read, coverage_cut);
 
     clean_weak_ma_hit_t(sources, reverse_sources, n_read);
@@ -29592,7 +29585,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut);
     ///print_binned_reads(sources, n_read, coverage_cut);
     ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
-
     ///fix_binned_reads(sources, n_read, coverage_cut);
     ///just need to deal with trio here
     ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
@@ -29601,7 +29593,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     ///debug_info_of_specfic_node((char*)"m64043_200504_050026/93784180/ccs", sg, ruIndex, (char*)"sbsbsb");
     init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq);
 
-
     asg_arc_del_trans(sg, gap_fuzz);
 
     asm_opt.coverage = get_coverage(sources, coverage_cut, n_read);
@@ -29617,7 +29608,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     asg_cut_tip(sg, asm_opt.max_short_tip);
     ///debug_info_of_specfic_node("m64043_200505_112554/8849050/ccs", sg, "inner_1");
     ///drop_inexact_edegs_at_bubbles(sg, bubble_dist);
-
     
     if(clean_round > 0)
     {
@@ -29696,7 +29686,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
             asg_cut_tip(sg, asm_opt.max_short_tip);
         }
     }
-
     if(VERBOSE >= 1)
     {
         fprintf(stderr, "\n\n**********final clean**********\n");
@@ -29715,7 +29704,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     asg_arc_del_orthology_multiple_way(sg, reverse_sources, 0.4, asm_opt.max_short_tip, ruIndex);
     asg_cut_tip(sg, asm_opt.max_short_tip);
 
-    
 
 
 
@@ -29725,7 +29713,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     asg_cut_tip(sg, asm_opt.max_short_tip);
 
     asg_arc_del_simple_circle_untig(sources, coverage_cut, sg, 100, 0);
-
     ///note: don't apply asg_arc_del_too_short_overlaps() after this function!!!!
     rescue_contained_reads_aggressive(NULL, sg, sources, coverage_cut, ruIndex, max_hang_length, 
     mini_overlap_length, 10, 1, 0, NULL, NULL, &b_mask_t);
@@ -29738,13 +29725,13 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     // rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length, 
     // mini_overlap_length, bubble_dist, 10);
 
+
     set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, 
     max_hang_length, mini_overlap_length);
 
     ruIndex->is_het = rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, 
     ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz, &b_mask_t);
 
-    
     output_unitig_graph(sg, coverage_cut, o_file, sources, ruIndex, max_hang_length, mini_overlap_length);
     flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
     
diff --git a/Overlaps.h b/Overlaps.h
index 33d20c9..761b26e 100644
--- a/Overlaps.h
+++ b/Overlaps.h
@@ -1160,7 +1160,7 @@ void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
 ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
 R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
-kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans);
+kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans, uint32_t collect_p_trans_f);
 ma_ug_t* copy_untig_graph(ma_ug_t *src);
 ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, 
 uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp
index 8e0beef..9068cdd 100644
--- a/Purge_Dups.cpp
+++ b/Purge_Dups.cpp
@@ -5156,7 +5156,7 @@ void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, hap_cov_t *cov,
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
 uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, 
-uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
+uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t collect_p_trans_f)
 {
     p_g_t *pg = NULL;
     asg_t* nsg = ug->g;
@@ -5235,16 +5235,20 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
 
     normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
 
+    if(collect_p_trans && collect_p_trans_f == 0)
+    {
+        collect_purge_trans_cov(ug, &all_ovlp, cov, position_index);
+    } 
+    
     if(asm_opt.polyploidy <= 2)
     {
         mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL);
         ///pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag);
     } 
     
-    if(collect_p_trans)
+    if(collect_p_trans && collect_p_trans_f == 1)
     {
         collect_purge_trans_cov(ug, &all_ovlp, cov, position_index);
-        ///goto end_coverage;
     } 
     
     pg = init_p_g_t(ug, cov, read_g);
diff --git a/Purge_Dups.h b/Purge_Dups.h
index 6d8d100..337f590 100644
--- a/Purge_Dups.h
+++ b/Purge_Dups.h
@@ -66,7 +66,7 @@ typedef struct {
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
 uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, 
-uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans);
+uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t collect_p_trans_f);
 void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
 uint32_t is_circle, uint64_t* rLen);
 void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);
diff --git a/hic.cpp b/hic.cpp
index bad0f86..f674fbd 100644
--- a/hic.cpp
+++ b/hic.cpp
@@ -2267,7 +2267,7 @@ void dfs_bubble(asg_t *g, kvec_t_u32_warp* stack, kvec_t_u32_warp* result, uint3
 uint32_t get_unitig_het_arb(ma_ug_t* ug, uint32_t uid, uint8_t *r_het_flag, kv_u_trans_t *ref, uint32_t m_het_occ,
 uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label)
 {
-    if(u_trans_n(*ref, uid) > 0) return m_het_label;
+    if(ref && u_trans_n(*ref, uid) > 0) return m_het_label;
     ma_utg_t *u = &(ug->u.a[uid]);
     uint32_t k, rId;
     uint32_t het_occ, hom_occ;
diff --git a/rcut.cpp b/rcut.cpp
index 66384ec..66657a0 100644
--- a/rcut.cpp
+++ b/rcut.cpp
@@ -28,6 +28,8 @@ uint8_t bit_filed[8] = {1, 2, 4, 8, 16, 32, 64, 128};
 typedef struct {
 	int32_t max_iter;
     int32_t n_perturb;
+	int32_t n_b_perturb;
+	int32_t n_s_perturb;
     double f_perturb;
     uint64_t seed;
 } mc_opt_t;
@@ -111,13 +113,12 @@ typedef struct {
 void mc_opt_init(mc_opt_t *opt, int32_t n_perturb, double f_perturb, uint64_t seed)
 {
 	memset(opt, 0, sizeof(mc_opt_t));
-	// opt->n_perturb = 50000;
 	opt->n_perturb = n_perturb;
-	// opt->f_perturb = 0.1;
 	opt->f_perturb = f_perturb;
 	opt->max_iter = 1000;
-	// opt->seed = 11;
 	opt->seed = seed;
+	opt->n_s_perturb = n_perturb;
+	opt->n_b_perturb = n_perturb*0.5;
 }
 
 void mc_merge_dup(mc_g_t *mg) // MUST BE sorted
@@ -1774,6 +1775,7 @@ static void mc_perturb_node(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_
 {
 	uint32_t i, k, n_bfs = 0;
 	k = (uint32_t)(kr_drand_r(&b->x) * b->cc_size + .499);
+	if(k >= b->cc_size) k = b->cc_size - 1;
 	k = (uint32_t)ma->cc[b->cc_off + k];///node id
 	n_bfs = mc_bfs(ma, b, k, bfs_round, (int32_t)(b->cc_size * opt->f_perturb));
 	for (i = 0; i < n_bfs; ++i)
@@ -1808,6 +1810,7 @@ static void mb_perturb_node(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b, int
 {
 	uint32_t i, k, n_bfs = 0;
 	k = (uint32_t)(kr_drand_r(&b->x) * b->cc_size + .499);
+	if(k >= b->cc_size) k = b->cc_size - 1;
 	k = (uint32_t)mbg->e->cc[b->cc_off + k];///node id
 	n_bfs = mb_bfs(mbg->e, b, k, bfs_round, (int32_t)(b->cc_size * opt->f_perturb));
 	for (i = 0; i < n_bfs; ++i)
@@ -2052,11 +2055,51 @@ void mc_solve_bp(mc_bp_t *bp)
 	fprintf(stderr, "[M::%s::%.3f] ==> round %u\n", __func__, yak_realtime()-index_time, r);
 }
 
-void print_sc(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b, t_w_t sc_opt, uint32_t n_iter)
+t_w_t mc_score_all_advance(const mc_match_t *ma, int8_t *s)
 {
-	t_w_t w = mc_score(ma, b);
-	if(w != sc_opt) fprintf(stderr, "ERROR\n");
-	fprintf(stderr, "# iter: %u, sc_opt: %f, sc-local: %f, sc-global: %f\n", n_iter, sc_opt, w, mc_score_all(ma, b));
+	uint32_t k;
+	t_w_t z[2], zt = 0;
+	for (k = 0; k < ma->n_seq; ++k) 
+	{
+		uint32_t o = ma->idx.a[k] >> 32;
+        uint32_t j, n = (uint32_t)ma->idx.a[k];
+		z[0] = z[1] = 0;
+		for (j = 0; j < n; ++j) {
+            const mc_edge_t *e = &ma->ma.a[o + j];
+            uint32_t t = ma_y(*e);
+            if (s[t] > 0) z[0] += e->w;
+            else if (s[t] < 0) z[1] += e->w;
+        }
+		zt += -((t_w_t)(s[k])) * (z[0] - z[1]);
+	}
+	return zt;
+}
+
+t_w_t mb_score_all_advance(const mc_match_t *ma, mb_g_t *mbg)
+{
+	uint32_t k;
+	t_w_t z[2], zt = 0;
+	for (k = 0; k < ma->n_seq; ++k) 
+	{
+		uint32_t o = ma->idx.a[k] >> 32;
+        uint32_t j, n = (uint32_t)ma->idx.a[k];
+		z[0] = z[1] = 0;
+		for (j = 0; j < n; ++j) {
+            const mc_edge_t *e = &ma->ma.a[o + j];
+            uint32_t t = ma_y(*e);
+            if (mbg->u->u.a[mbg->u->idx.a[t]>>1].s[mbg->u->idx.a[t]&1] > 0) z[0] += e->w;
+            else if (mbg->u->u.a[mbg->u->idx.a[t]>>1].s[mbg->u->idx.a[t]&1] < 0) z[1] += e->w;
+        }
+		zt += -((t_w_t)(mbg->u->u.a[mbg->u->idx.a[k]>>1].s[mbg->u->idx.a[k]&1])) * (z[0] - z[1]);
+	}
+	return zt;
+}
+
+void print_sc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, t_w_t sc_opt, uint32_t n_iter)
+{
+	t_w_t w = mc_score(mg->e, b);
+	// if(w != sc_opt) fprintf(stderr, "ERROR\n");
+	fprintf(stderr, "# iter: %u, sc_opt: %f, sc-local: %f, sc-global: %f\n", n_iter, sc_opt, w, mc_score_all_advance(mg->e, mg->s.a));
 }
 
 uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint32_t cc_off, uint32_t cc_size)
@@ -2065,8 +2108,8 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
 	t_w_t sc_opt = -(1<<30), sc;///problem-w
 	b->cc_off = cc_off, b->cc_size = cc_size;
 	if (b->cc_size < 2) return 0;
-	// print_sc(opt, mg->e, b, sc_opt, (uint32_t)-1);
 	sc_opt = mc_init_spin(mg->e, b);
+	// print_sc(opt, mg, b, sc_opt, n_iter);
 	if (b->cc_size == 2) return 0;
 	for (j = 0; j < b->cc_size; ++j) {///backup s and z in s_opt and z_opt
 		b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; ///hap status of each unitig
@@ -2086,6 +2129,7 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
 			b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
 		}
 	}
+	// print_sc(opt, mg, b, sc_opt, n_iter);
 	// mc_reset_z_debug(mg->e, b);
 	// print_sc(opt, mg->e, b, sc_opt, n_iter);
 	// fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off);
@@ -2100,6 +2144,7 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
 				b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
 			}
 			sc_opt = sc;
+			// print_sc(opt, mg, b, sc_opt, n_iter);
 		} else {
 			for (j = 0; j < b->cc_size; ++j) {
 				b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
@@ -2117,11 +2162,14 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
 			}
 			sc_opt = sc;
 		}
-		
-		// print_sc(opt, mg->e, b, sc_opt, n_iter);
 	}
+
 	for (j = 0; j < b->cc_size; ++j)
+	{
 		b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
+		b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
+	}
+		
 	return n_iter;
 }
 
@@ -2291,6 +2339,7 @@ void mc_init_spin_all(const mc_opt_t *opt, mc_g_t *mg, mb_g_t *mbg, mc_svaux_t *
 	}
 
 	if(!mbg) return;
+	// memcpy(b->s_opt, b->s, sizeof(int8_t)*mg->e->n_seq);
 
 	///adjust by block
 	kvec_t(uint32_t) s; kv_init(s);
@@ -2385,6 +2434,7 @@ void mc_init_spin_all(const mc_opt_t *opt, mc_g_t *mg, mb_g_t *mbg, mc_svaux_t *
 	/*******************************for debug************************************/
 	// debug_mbg(mbg, mg, b);
 	/*******************************for debug************************************/
+	// memcpy(b->s, b->s_opt, sizeof(int8_t)*mg->e->n_seq);
 }
 
 
@@ -2448,6 +2498,9 @@ void mb_g_cc(mb_g_t *mbg)
 
 void mc_set_by_mbg(mc_g_t *mg, mb_g_t *mbg)
 {
+	// t_w_t z0 = mc_score_all_advance(mg->e, mg->s.a);
+	// t_w_t z1 = mb_score_all_advance(mg->e, mbg);
+	// if(z0 >= z1) return;
 	uint32_t i, k, qn, *a[2], a_n[2];
 	int8_t s[2];
 	for (i = 0; i < mbg->u->u.n; i++)
@@ -2502,27 +2555,9 @@ void print_mb_g_blcok(mb_g_t *mbg)
 	}
 }
 
-t_w_t mc_score_all_advance(const mc_match_t *ma, int8_t *s)
-{
-	uint32_t k;
-	t_w_t z[2], zt = 0;
-	for (k = 0; k < ma->n_seq; ++k) 
-	{
-		uint32_t o = ma->idx.a[k] >> 32;
-        uint32_t j, n = (uint32_t)ma->idx.a[k];
-		z[0] = z[1] = 0;
-		for (j = 0; j < n; ++j) {
-            const mc_edge_t *e = &ma->ma.a[o + j];
-            uint32_t t = ma_y(*e);
-            if (s[t] > 0) z[0] += e->w;
-            else if (s[t] < 0) z[1] += e->w;
-        }
-		zt += -((t_w_t)(s[k])) * (z[0] - z[1]);
-	}
-	return zt;
-}
 
-void mb_solve_core(const mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys)
+
+void mb_solve_core(mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys)
 {
 	if(!ref) return;
 	double index_time = yak_realtime();
@@ -2530,11 +2565,11 @@ void mb_solve_core(const mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t
 	mb_g_t *mbg = init_mb_g_t(mg, ref, is_sys);
 	mb_svaux_t *bb;
 	/**************************init**************************/
+	fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
 	mc_svaux_t *b;
 	mc_g_cc(mg->e);
 	b = mc_svaux_init(mg, opt->seed);
 	mc_init_spin_all(opt, mg, mbg, b);
-	
 	mc_svaux_destroy(b);
 	free(mg->e->cc); 
 	mg->e->cc = NULL;
@@ -2543,18 +2578,21 @@ void mb_solve_core(const mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t
 	bb = mb_svaux_init(mbg, opt->seed);
 	/*******************************for debug************************************/
 	// print_mb_g_blcok(mbg);
+	fprintf(stderr, "*********before-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
+	fprintf(stderr, "*********before-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg));
 	/*******************************for debug************************************/
-
-	fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
+	opt->n_perturb = opt->n_b_perturb;
 	for (st = 0, i = 1; i <= mbg->e->n_seq; ++i) {
 		if (i == mbg->e->n_seq || mbg->e->cc[st]>>32 != mbg->e->cc[i]>>32) {
 			mb_solve_cc(opt, mbg, bb, st, i - st);
 			st = i;
 		}
 	}
-	
+	opt->n_perturb = opt->n_s_perturb - opt->n_b_perturb;
 	/*******************************for debug************************************/
 	// debug_mb_solve_core(mbg);
+	fprintf(stderr, "*********after-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
+	fprintf(stderr, "*********after-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg));
 	/*******************************for debug************************************/
 	mc_set_by_mbg(mg, mbg);
 	fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
@@ -2583,7 +2621,7 @@ void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
 			st = i;
 		}
 	}
-	fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
+	fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b));
 	if(bp) mc_solve_bp(bp);	
 	///mc_write_info(g, b);
 	mc_svaux_destroy(b);
@@ -2694,6 +2732,60 @@ void p_nodes(mc_g_t *mg, trans_chain* t_ch, uint8_t* trio_flag)
 	}
 }
 
+void write_mc_g_t(mc_opt_t *opt, mc_g_t *mg, const char *name)
+{
+	FILE* fp = fopen(name, "w");
+
+    fwrite(opt, sizeof(mc_opt_t), 1, fp);
+	fwrite(&(mg->s.n), sizeof(mg->s.n), 1, fp);
+	fwrite(mg->s.a, sizeof(mc_node_t), mg->s.n, fp);
+	fwrite(&(mg->e->n_seq), sizeof(mg->e->n_seq), 1, fp);
+	fwrite(&(mg->e->ma.n), sizeof(mg->e->ma.n), 1, fp);
+	fwrite(mg->e->ma.a, sizeof(mc_edge_t), mg->e->ma.n, fp);
+	fwrite(&(mg->e->idx.n), sizeof(mg->e->idx.n), 1, fp);
+	fwrite(mg->e->idx.a, sizeof(uint64_t), mg->e->idx.n, fp);
+
+    fclose(fp);
+}
+
+mc_g_t* load_mc_g_t(mc_opt_t *opt, const char *name)
+{
+	FILE* fp = NULL; 
+    fp = fopen(name, "r"); 
+    if(!fp) return NULL;
+
+	uint64_t flag = 0;
+	mc_g_t *mg = NULL; CALLOC(mg, 1);
+	kv_init(mg->s); CALLOC(mg->e, 1);
+
+    flag += fread(opt, sizeof(mc_opt_t), 1, fp);
+
+	flag += fread(&(mg->s.n), sizeof(mg->s.n), 1, fp);
+	mg->s.m = mg->s.n; MALLOC(mg->s.a, mg->s.n);
+	flag += fread(mg->s.a, sizeof(mc_node_t), mg->s.n, fp);
+
+	flag += fread(&(mg->e->n_seq), sizeof(mg->e->n_seq), 1, fp);
+
+	flag += fread(&(mg->e->ma.n), sizeof(mg->e->ma.n), 1, fp);
+	mg->e->ma.m = mg->e->ma.n; MALLOC(mg->e->ma.a, mg->e->ma.n);
+	flag += fread(mg->e->ma.a, sizeof(mc_edge_t), mg->e->ma.n, fp);
+
+	flag += fread(&(mg->e->idx.n), sizeof(mg->e->idx.n), 1, fp);
+	mg->e->idx.m = mg->e->idx.n; MALLOC(mg->e->idx.a, mg->e->idx.n);
+	flag += fread(mg->e->idx.a, sizeof(uint64_t), mg->e->idx.n, fp);
+
+    fclose(fp);
+	return mg;
+}
+
+void debug_mc_g_t(const char* name)
+{
+	mc_opt_t opt;
+	mc_g_t *mg = load_mc_g_t(&opt, name);
+	mc_solve_core(&opt, mg, NULL);
+	destory_mc_g_t(&mg);
+	exit(1);
+}
 
 void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref)
 {
@@ -2704,6 +2796,7 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u
 	
 	mb_solve_core(&opt, mg, ref, is_sys);
 	///debug_mc_g_t(mg);
+	if(renew_s == 0) write_mc_g_t(&opt, mg, MC_NAME);
 	mc_solve_core(&opt, mg, bub);
 
 	if((asm_opt.flag & HA_F_PARTITION) && t_ch)
@@ -2713,6 +2806,5 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u
 
 	if(ovlp) clean_ovlp_by_mc(mg, ovlp);
 
-	destory_mc_g_t(&mg);
-	
+	destory_mc_g_t(&mg);	
 }
\ No newline at end of file
diff --git a/rcut.h b/rcut.h
index 42a2bc0..64ee4e6 100644
--- a/rcut.h
+++ b/rcut.h
@@ -22,7 +22,7 @@ typedef struct {
 #define w_t double
 #define t_w_t double
 #define w_cast(x) ((t_w_t)((x)))
-
+#define MC_NAME "debug_mc.bin"
 
 typedef struct {
 	uint64_t x; ///(uint64_t)nid1 << 32 | nid2;
@@ -89,4 +89,5 @@ static inline double kr_drand_r(uint64_t *x)
 }
 
 void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref);
+void debug_mc_g_t(const char* name);
 #endif
\ No newline at end of file